1 VOLUME-OF-FLUID MODEL COMFLOW SIMULATIONS OF WAVE IMPACTS ON A DIKE

ComFLOW is a 3D Volume-of-Fluid (VOF) model to solve the incompressible Navier-Stokes equations including free surface, or to solve the Navier-Stokes equations for two-phase flow problems (two-phase flow: both an incompressible viscous fluid (e.g., water) and a compressible viscous fluid (e.g., air) are present). The problem statement of the present study reads: ‘Is ComFLOW capable of accurate prediction of wave impacts on (impermeable) coastal structures such as dikes? And, if so, what are the preferred model settings and associated computing times?’. In this paper, ComFLOW is validated for this purpose by comparison against pressure data as measured in the Delta flume by pressure sensors at dikes. We have selected three different experiments, with typical dike geometries (slope 1:3.5, with and without berm) at which more than 20 pressure sensors were installed. The results can be summarized as follows. The pressure measurements are reproduced well in the simulations. A grid with about 170 grid cells per wave length in the horizontal, and between 4 and 6 grid cells per wave height in the vertical, proves to be sufficiently fine. At such a grid resolution and with about 450 by 35 grid cells in the computational domain, a typical CPU time is 35 minutes for simulations with a model time of 10 wave periods. For the present application, it is preferable to use the one-phase flow model rather than the two-phase flow model, since the former gives better results in the lower located pressure sensors and consumes less CPU time.


INTRODUCTION
Wave fields in the close vicinity (a few wave lengths) of a coastal structure can be very complex.This is especially the case during extreme storm conditions and for complex geometries.Wave breaking, the generation of super-and sub-harmonics, wave set-down and set-up, shoaling and refraction all significantly transform the wave field nearby the structure.Obviously, these wave transformations affect the subsequent wave impacts on the structure.These impacts may cause damage to the structure or may even lead to complete failure.Also processes related to wave impacts, like wave run-up on and wave overtopping, are strongly influenced by these wave field transformations.Since these transformations are highly nonlinear, analytical models are in most cases inadequate to provide accurate information.Therefore, the majority of studies in this field rely on laboratory experiments in flumes or basins, or on design rules that have resulted from experiments executed at earlier stages.However, laboratory experiments suffer from several deficiencies, such as scale effects, limited reproducibility, and high financial costs.As a result, there is a large interest in a complementary approach: computational models capable of accurate and efficient simulation of the complex wave field near, and the wave impacts on coastal structures.See, to mention just one example, the work on Cobras in Losada et al. (2008) and references therein.In the present paper, Deltares presents results obtained with the Volume-of-Fluid (VOF) simulation model ComFLOW.
The problem statement of the present study is formulated as follows: 'Is ComFLOW capable of accurate prediction of wave impacts on (impermeable) coastal structures such as dikes?And, if so, what are the preferred model settings and associated computing times?'.An answer to the problem statement is obtained by comparing numerical results against laboratory measurement data of pressure sensors located at a dike under severe wave attack.

DESCRIPTION OF COMFLOW
ComFLOW can be run as a one-phase or a two-phase flow model.The one-phase flow model, see Kleefsman et al. (2005), is a Volume-of-Fluid model that solves the incompressible Navier-Stokes equations including free surface motion.The two-phase flow model, see Wemmenhove (2008), is a model to solve the Navier-Stokes equations for problems in which both an incompressible viscous fluid (e.g., water) and a compressible viscous fluid (e.g., air) are present.The two-phase model option is particularly suited to account for the cushioning effect of air entrapment during wave impact.The detailed description of the free surface in the one-phase model, or the water-air interface in the twophase model, allows for accurate simulation of breaking waves.ComFLOW can be applied in 2DV ('flume') and in 3D ('basin') mode.Some essential information on ComFLOW is given in the present paper.More details can be found in the given references.

Governing equations
The motion of a homogeneous, incompressible, viscous fluid (e.g.water) is described by the continuity equation (conservation of mass) and the Navier-Stokes equations (conservation of momentum): Here, u is the velocity vector, the fluid density, p the pressure, the dynamic viscosity coefficient, and F the external force vector (e.g. for gravity one has , , 0, 0, At all boundaries in contact with fluid, boundary conditions are required.At solid domain boundaries, the no-slip condition is applied.At the free surface, kinematic and dynamic boundary conditions are imposed.The issue of wave boundary conditions is dealt with further on.
In the two-phase flow model, to close the system of equations the adiabatic equation between pressure and air is used: Extension of ComFLOW with a model to account for wave interaction with permeable structures is discussed in Wellens et al. (2010).

Numerical model
ComFLOW uses simple rectangular (Cartesian) grids, which has several advantages compared to other grid types (e.g.unstructured grids): a more accurate and sharp treatment of the free surface; easier grid generation; simpler and more efficient data structures.The grid may be uniform (all grid cells have the same dimension) or stretched (in one or more dimensions the grid cell dimensions vary, though the cells continue to be rectangular-shaped).The employed cut-cell technique, in which volume and edge apertures indicate which part of the cell volume respectively cell edge are open for fluid, allows for the incorporation of arbitrarily shaped bodies.Cell labeling is introduced to distinguish between cells of different character: B(oundary) cells are completely filled by an object (hence, the volume aperture is zero); E(mpty) cells are empty, but have the possibility of letting fluid in; S(urface) cells are adjacent to E(mpty) cells and contain fluid; the remaining cells are labeled as F(luid) cells.This is illustrated in Fig. 1.A staggered positioning of the primary variables (pressure and velocity) on the grid is used, which avoids the need to apply artificial remedies against unwanted pressure-velocity oscillations (so-called 2 x waves).The continuity and Navier-Stokes equations are discretized using the finite volume method.This ensures conservation of mass and momentum in F-cells; in S-cells, this is not necessarily the case, see Kleefsman et al. (2005).The discretization of the convection term is skew symmetric.This is a physical property ensuring that the discrete convection term does not generate nor destroy energy, but merely redistributes it.Besides the physical kinematic viscosity / , an amount of numerical viscosity equal to that of a first-order upwind scheme is added to the diffusive term.
At some numerical (i.e., open) domain boundaries, incident waves should be able to enter the domain.Simultaneously, outgoing waves should be able to leave the domain through all numerical boundaries without spurious reflections.The Generating and Absorbing Boundary Condition (Wellens et al. (2009)  The equations of motion are discretized in time using the first-order accurate forward Euler method.Convection and diffusion are treated explicit in time.For the pressure, which is treated implicit in time, a Poisson equation needs to be solved.The time step t is automatically adapted during the simulation: it is halved when the CFL-number exceeds a certain threshold cflmax (must be specified by the user; should be <1), or doubled when the CFL-number is smaller than a user-specified value cflmin for ten successive time steps.

DESCRIPTION OF LABORATORY EXPERIMENTS
For validation of ComFLOW for wave impacts on dikes, we selected the experimental data from three laboratory experiments performed between 1997 and 2003 in the Delta flume of Deltares.The selection of the experiments is such that (i) a rather large range of wave steepnesses is considered, and (ii) a large number of pressure sensors was installed at the structure's slope.In other words, the data set, summarized in Table 1, is extensive.
1.17  From each experiment, we have selected four time windows, indicated by t1 to t4.This means that ComFLOW is validated for in total 12 different time windows (4 time windows in 3 experiments).The selection of the time windows was made such that the largest wave crest in the incoming signal (time window t1), the largest pressure head at the structure (time window t2), and the largest block motion (as deduced from another numerical model, ZSTEEN) (time window t3) are included.The duration of these time windows lies between 50 s and 90 s.Time window t4 has a long duration (500 s), and is included to detect whether unwanted numerical effects (e.g.spurious leakage or inflow, which changes the mean water level) accumulate slowly.
The incoming wave signal for the ComFLOW simulations is obtained from measured wave data at three wave gauges located at some distance from the toe of the structure: 62 m for #21o12 and #23o06, and 81 m for #P022.
Experiments #21o12 and #23o06 both have, as shown in Fig. 2, the same dike geometry, with slope 1:3.5 and a berm present.The water level is different for these two experiments.The geometry as used in experiment #P022 also has a slope of 1:3.5, but in this experiment the berm is absent.During the considered experiments, more than 20 pressure sensors are installed roughly between 0.3 m and 1.5 m below the still water level.

COMFLOW MODEL SET-UP AND SIMULATION SERIES
All simulations are done in 2DV, which means a horizontal ( x ) direction and a vertical ( z ) direction are included.The other horizontal ( y ) direction does not play a role, and will be discarded from now on.The horizontal domain length, with 0 x at the wave gauges' location and x positive towards the dike geometry, is 95 m for experiments #21o012 and #23o06, and 110 m for experiment #P022.The vertical domain size, with z pointing upwards from the bottom at zh (hence, 0 z corresponds to the still water level), is 8 m for all experiments.Insertion of the dike geometry in ComFLOW input files turns out to be trivial.
Various computational grids are used, labeled as 'coarse' (short for coarse grid), 'medium' (short for medium grid, i.e. with a grid resolution between that of a coarse and fine grid), 'fine' (short for fine grid) and 'stretched' (short for a stretched grid), see Table 2.The meaning of the symbols is as follows: x and z stand for the grid cell sizes in the horizontal and vertical direction;  The grids 'coarse', 'medium' and 'fine' are uniform, which means that all grid cells have the same size.On uniform grids, the cell aspect ratio is taken equal to one.This means that the horizontal and vertical grid cell sizes are taken equal: xz .The motivation for this is that, at least near the zone where wave impacts occur, the horizontal and vertical length scales of the hydrodynamics have about the same magnitude.To get an idea of the employed grid resolutions as compared to the dimensions of the relevant physical parameters, consider the following.The number of horizontal grid cells per wavelength p L (see Table 1) is roughly 85 for the 'coarse' grids, 170 for the 'medium' grids, and 340 for the 'fine' grids.The number of vertical grid cells per wave height s H (see Table 1) ranges roughly between 2 and 3 for 'coarse' grids', between 4 and 6 for 'medium' grids, and between 7 and 12 for 'fine' grids.
A stretching in both the horizontal and vertical direction is applied for the grids labeled 'stretched'.The reason for studying stretched grids is that they may yield a desired level of accuracy using less computational resources, by choosing the grid resolution based on (expected) complexity of the wave field.For the present application, the wave field is expected to be most complex near the wave impact location at the dike (here we put the 'stretched' grid resolution equal to that of the 'fine' grid), while at the incoming wave boundary and near the bottom the wave field (where we put the 'stretched' grid resolution equal to that of the 'medium' grid) is expected to be less complex.
At An overview of the performed simulations is given in Table 3.Each series contains simulations of the three experiments and for the mentioned time windows.For example, series B contains simulations for the following experiments and time windows: #21o12-t1, #21o12-t2, #21o12-t3, #23o06-t1, #23o06-t2, #23o06-t3, #P022-t1, #P022-t2, #P022-t3.It can be deduced from Table 3 that in total 63 simulations (9 in each of the series B, C, G, H, J and K, and 3 in each of the series L, M and N) are performed.All series but one are one-phase flow simulations; series M contains two-phase flow simulations.The water level is at rest at the beginning of the simulation.As an immediate consequence, no wave impact at the structure occurs during the first part of the simulation.The reason is that the wave field needs time to propagate from the boundary to the structure.
The water density is taken equal to = 1000 kg/m 3 and the dynamic viscosity of water is taken equal to = 1.0 x 10 -3 kg / (m s).Kinematic surface tension is neglected, and the gravity constant is taken equal to g = 9.81 m/s 2 .In the two-phase flow simulations, the air density ref is taken equal to 1.0 kg/m 3 , the dynamic viscosity of air is taken equal to 1.7 x 10 -5 kg / (m s), the atmospheric pressure 0 p is taken equal to 1.0 x 10 5 Pa, and the adiabatic coefficient is put equal to 1.4.Values for cflmin and cflmax for the automatic selection of the time step are taken equal to 0.4 and 0.9 respectively.All simulations are performed on the same PC.The operating system on this PC is Windows XP, version 2002, with Service Pack2.The PC has 2.00 GB of RAM.It has an AMD Athlon™ 64 processor, 4000+, and runs at a clock frequency of 2.41 GHz.

DISCUSSION OF THE RESULTS
Validation is done by comparing numerical results against measured data in the pressure sensors, both qualitatively (by visual inspection of figures) and quantitatively (by studying the root-mean square error of the pressure).Because of space limitations, only a small part of the results is included in the present paper.The discussion of the results is distributed over several topics.

Some overall remarks
Of the 63 simulations, 58 finished successfully, while 5 simulations did not finish well.Of these unsuccessful ones, 3 simulations crashed that were executed with the GABC accurate celerity settings, and 2 simulations ended because of a reduction of the time step to unacceptably small values.The reason for these 5 crashes is yet unknown.Of the 48 simulations done with the GABC with constant celerity settings, 46 finished successfully while 2 simulations crashed.This is considered as an acceptable score for a model still in development.
The overall trends in the pressure measurements of all experiments and in all time windows are reproduced in the simulations.This holds for all grids, even the coarse ones.No slowly growing unwanted numerical effects were detected during the simulations, including those of time window t4.The wave impacts, characterized by a sudden and large increase in pressure, are in most cases accurately captured.This is a rewarding conclusion!For illustration purposes, a snapshot of the computed spatial velocity and pressure distribution is given in Fig. 3, together with the measured and computed pressure head.The pressure head is defined as / p gz .
ComFLOW deals accurately with the process of flooding-and-drying.In some instances, the computed pressure wiggles rapidly between positive and negative values (see Fig. 4, in panel #23o06-t1 dro 31; dro31 means pressure sensor 31).This effect is due to the presence of small droplets, and becomes less with grid refinement.

Grid refinement and grid stretching
Grid refinement improves the quality of the numerical results.The quality obtained at the stretched grid lies, as expected, between that of the medium grid and the fine grid.Some results are shown in Fig. 4, where we omitted the first part of the simulation in which the wave field is not fully developed yet.As mentioned above, the overall trends in the pressure measurements of all experiments and in all time windows are reproduced in the simulations, even on the coarse grids.
The accuracy of the coarse grid solution is for some simulations, and in particular in the higher located pressure sensors, perhaps not good enough; see for example #23o06-t1 dro 31 and #P022-t1 dro 22.The results on the medium grid seems to be good enough.
A quantitative measure to compare (per simulation) the computed and measured pressures is introduced: the root-mean square error (RMSE).The RMSE of the pressure is computed as follows:   It appears, see Table 4, that the reduction in error is less than the reduction in grid size (i.e., no factor 2 in error reduction for a factor 2 in cell size reduction).In other words, in the present study the quality improvement is not very dramatic with grid refinement.It is possible that part of this (relative) insensitivity to grid refinement is due to systematic errors, e.g. in the measurements or in the model.Also for this reason, we think the medium grid solution can be considered as good enough, since further grid refinement leads to only relatively small improvements.

One-phase and two-phase flow simulations
The one-phase and two-phase flow model results show the same trends.In the lower located pressure sensors, the two-phase flow model is less accurate than the one-phase flow model in reproducing the quick water retreat (downwash) before impact, see also Fig. 5, left panel.The CPU time of the two-phase flow model is larger than that of the one-phase flow model, see below.

Boundary conditions
The best boundary condition in terms of robustness and accuracy is formed by the GABC with constant celerity settings.These settings are such that no spatial derivatives in the vertical need to be taken.The GABC with accurate celerity settings turn out not to be robust enough (3 crashed out of 9 simulations).The inflow type boundary condition turns out, as expected, to be less accurate from the moment that wave reflections from the structure reach the incident wave boundary.For reasons of limited space, no figures are included to confirm this claim.

CPU time
The CPU time for a limited set of series is shown in Table 5.Some metadata on the series already given in Table 3 is included again in Table 5 as a service to the reader.The employed time step is an important factor contributing to the CPU time.As discussed earlier, the time step is automatically adapted when the actual CFL number reaches certain thresholds.For all simulations, it turns out that the time steps hardly varies during the simulation (with the crashed simulations being notable exceptions).This is an indication of robustness of the simulations: significant variations in the time step would be an indication for areas with large, possibly spurious, local velocities.The values of the time step as occurring during the major part of the simulations are given in Table 6.As expected, the time step decreases at finer grids.Apparently, the computed velocities are such that the time step is the same for the medium and the coarse grid for #21o12 and #23o06.The time steps at the stretched and fine grid are identical.This is as expected, since the smallest grid cells determine, through the CFL condition, the time step.
Studying Table 5 and Table 6 lead to the conclusion that, for the simulations with GABC and the studied grids, the CPU time scales linearly with the number of grid cells and with the total number of time steps (i.e., the inverse of the time step).This conclusion may seem too good to be true; it is wellknown that the number of operations in a linear solver to solve for the pressure increases more rapidly that the number of unknowns (which is proportional to the number of grid cells).Investigations (not detailed further in this paper) indicate that with a relatively small number of grid cells (say, < 100.000) less than half of the CPU time is spent in the linear solver; the larger part is spent in other activities, such as computation of the matrix elements.Since the number of operations required for these activities scale linearly with the number of grid cells, the (nonlinear) increase in linear solver computational effort appears to be invisible for the grids studied.
Application of the two-phase flow model leads to an increase in CPU time of about a factor 1.5 compared to the one-phase flow model (compare series M and N).This factor is, perhaps not surprisingly, close to the ratio of active grid cells for the two-phase flow model (all grid cells containing water and air) and the one-phase flow model (all grid cells containing water).
Application of the GABC costs about a factor 1.4 more CPU time than application of the inflow type boundary condition (compare series C and H with series B, and compare series G with series N).This is attributed to the different pressure matrix structure required for the GABC, which necessitates the use of another pressure solver.The difference in CPU time between GABC constant celerity settings and GABC accurate celerity settings is negligible (compare series C and H).
Summarizing, sufficiently accurate 2DV simulations on a medium grid resolution and a domain equal to about 450 by 35 grid cells (two wave lengths long) for a duration of about 10 wave periods model time take about 35 minutes on an ordinary PC.

CONCLUSIONS
ComFLOW is a 3D Volume-of-Fluid (VOF) model to solve the incompressible Navier-Stokes equations including free surface, or to solve the Navier-Stokes equations for two-phase flow problems (both an incompressible viscous fluid (e.g., water) and a compressible viscous fluid (e.g., air) are present).The problem statement of the present study reads: 'Is ComFLOW capable of accurate prediction of wave impacts on (impermeable) coastal structures such as dikes?And, if so, what are the preferred model settings and associated computing times?' ComFLOW has been validated for this purpose by comparison against pressure data as measured in the Delta flume by pressure sensors at dikes.We have selected three different experiments, with typical dike geometries (slope 1:3.5, with and without berm) at which more than 20 pressure sensors were installed.From each experiment, we have selected four time windows.Various series, with different model settings (grid resolution, types of boundary conditions, one-phase and two-phase model), have been simulated.Validation is done by comparing numerical results against measured data in the pressure sensors, both qualitatively (by visual inspection of figures) and quantitatively.The results can be summarized as follows.ComFLOW is robust.The overall trends of the pressure measurements of all experiments and in all time windows are reproduced in the simulations.A grid with about 170 grid cells per wave length in the horizontal, and between 4 and 6 grid cells per wave height in the vertical, proves to be sufficiently fine.At such a grid resolution, with about 450 by 35 grid cells in the domain, a typical CPU time is 35 minutes for simulations taking 10 wave periods.The Generating and Absorbing Boundary Conditions with constant celerity settings proves to be sufficiently accurate.For the present application, it is preferable to use the one-phase flow model rather than the two-phase flow model, since the former gives better results in the lower located pressure sensors and consumes less CPU time.

Figure 1 .
Figure 1.Cell labelling: dark grey denotes solid body; light grey is liquid.

FF.
if the open part of the cell is completely filled with fluid, and 01 s F if the open part of the cell is partly filled with fluid.A piecewise constant reconstruction of the free surface is used, where the free surface is displaced by changing the s F value in a cell using the calculated fluxes through cell faces.In ComFLOW, an adapted version of the VOF method is used.The adaptation has been introduced to avoid two main drawbacks of the original VOF method: the occurrence of 'flotsam and jetsam' (loose droplets that separate from the liquid, being a numerical artefact of the liquid displacement), and the loss and gain of water due to rounding off the VOF-By combining the VOF method with a local height function, strict mass conservation is ensured, and almost no flotsam and jetsam occurs.These two properties are a prerequisite for accurate simulation of breaking waves.

N
stand for the number of grid cells in the horizontal and vertical direction;

Figure 3 .
Figure 3. Snapshot of the computed pressure distribution (upper panel; in kN/m 2 ), velocity distribution (middle panel; in m/s) and pressure head (lower panel) at a fine grid (series H) and the same time level, for experiment #21o12-t1.
pressure sensors, and the second summation runs over all time samples in the considered time interval.The considered time interval is chosen such that the first part of the simulation, in which the wave field is not fully developed yet, is omitted.The computed and measured pressures are indicated by comp p at an equidistant time step equal to 0.005 s (200 Hz), which is the sampling time (inverse of the sampling frequency) of the measurements.Interpolation in time of the computed pressures gives the values of comp p at the measurement time steps.

Figure 4 .
Figure 4. Results of the grid refinement study showing the pressure (in kN/m 2 ) as function of time for some experiments (#21o12-t1, #23o06-t1 and #P022-t1) and some pressure sensors (the higher the number after dro, the higher its location on the slope).All simulations are done with GABC constant celerity settings and the onephase flow model.Measurements (black continuous lines) vs computed results: coarse grid (black dashed lines -series J), medium grid (green lines -series G), stretched grid (blue lines -series K) and fine grid (red lines -series H).The horizontal dashed line in black indicates the hydrostatic pressure in the pressure sensor in still water.

Figure 5 .
Figure 5. Results of the one-phase flow and two-phase flow model pressure (in kN/m 2 ) as function of time for experiments #21o12-t1 in some pressure sensors (the higher the number after dro, the higher its location on the slope).All simulations are done with inflow type boundary condition and on a medium grid.Measurements (black continuous lines) vs computed results: one-phase (red lines -series N) and two-phase (blue linesseries M).The horizontal dashed line in black indicates the hydrostatic pressure in the pressure sensor in still water.

Table 1 . Characteristics of selected experiments Test name
The meaning of the symbols in the table is as follows: s H represents the significant wave height; represents the peak period; the wave steepness is defined as / p T

Table 2 . Information with respect to the grids
the incoming wave boundary, three different types of boundary conditions are studied in the present paper: Inflow type boundary condition.As mentioned above, this boundary condition generates incident waves correctly, but does not absorb outgoing waves.GABC with accurate celerity settings.By putting 0 0 c a gh .