Development and implementation of a Direct Surface Description method for free surface flows in OpenFOAM

The solution procedure for two phases in OpenFOAM suffers from unphysical velocity oscillations at the free surface between the two phases. It is likely that this problem also exists in other two-phase Computational Fluid Dynamics (CFD) codes. We aim to solve this by imposing boundary conditions directly on the free surface. We have taken the first step towards a new two-phase solution method by first addressing the water phase alone. It is a free surface modelling method based on merging concepts from two existing methods: (1) A single-phase free surface method and (2) the solution method used in OpenFOAM. The underlying motivation is to enable more accurate estimation of wave induced load distributions from wave crest impacts on offshore structures. This first and foremost requires an accurate prediction of the kinematics near the free surface. We present a solution method with boundary conditions directly on the free surface, thereby the name: Direct Surface Description (DSD). Additionally it is the first time that the isoAdvector algorithm is combined with a single phase free surface method. The implementation is made in OpenFOAM, but may also relevant to other codes as well. First a still water level simulation is presented to illustrate the unphysical behaviour of the existing solvers and validate the behaviour of the DSD method. The second test case is a moderately steep stream function wave in intermediate water depth. The DSD method is validated and compared to the existing solution methods of OpenFOAM: interFoam and interIsoFoam . We present a detailed comparison of surface elevations and velocity profiles. This is followed by a convergence study including wave height, velocity and phase shift. Additionally the influence of the Courant–Friedrichs–Lewy (CFL) number is studied. The stream function wave case demonstrates that the DSD method accurately predicts the free surface elevation and velocity fields without free surface undulations or oscillatory velocity fields. The convergence study underlines an increased accuracy of the DSD method. Finally, a 2D and a 3D showcase with breaking waves are presented to show that the DSD method is capable of simulating more complex and realistic cases.


Introduction
The observation of extreme breaking waves from offshore structures at intermediate water depth in the Danish North sea has made the industry question the safety level and failure probability of the existing structures.Safe operations are important within the offshore industry, which deals with decommissioning and maintenance of oil and gas platforms as well as construction and maintenance of offshore wind turbines.Detailed analyses of wave loads from both breaking and nonbreaking waves are key to reduce the uncertainty of the wave load on individual beams in offshore structures.
Much research have investigated wave loads on monopiles, oriented both vertically as well as horizontally and substantial knowledge is available for the load on fully submerged cylinders.For semisubmerged structures like monopiles or jacket structures it is often  25 years (OpenFOAM-v1912, 2019.We chose to implement our work in OpenFOAM, because (1) it is popular in both academia and industry, (2) it has an active community and (3) the code is open-source and maintained with official releases every six months.The official release of OpenFOAM-v1912 provides two solvers called interFoam and interIsoFoam which are the common choices for simulation of water waves.
The numerical method in interFoam and interIsoFoam solves the governing equations with smooth density differences near the free surface.This is a well accepted method for smoothly varying densities, such as stratified ocean flows.This method was successfully used to study mixing of saline layers due to bridge piers by Jensen et al. (2018).However, the density jump at the interface between water and air is very large and the approach of solving the flow directly with density differences might not be adequate, as the density gradient is unresolved.The white cells are air cells.Wroniszewski et al. (2014) compared the performance of several codes including interFoam in a test case with a progressive solitary wave and highlighted that interFoam over-predicts the velocity at the wave crest relative to the analytical solution.The same problem was also highlighted by Roenby et al. (2017), Amini Afshar (2010), Tomaselli (2016) and Larsen et al. (2019).Amini Afshar (2010) and Larsen et al. (2019) noted a second problem in interFoam, which is generation of wiggles in the interface between air and water.A third problem is the development of spurious velocities in the air region above the free surface, which may adversely affect the velocities in the water near the crest.The spurious air velocity problem has received most attention as indicated by Larsen et al. (2019), who mentioned nine articles published in the time span from 1999 to 2008.All these studies attribute the development of spurious velocities to the surface tension model, however more recently it was shown that spurious air velocities can be generated without surface tension modelling (Wemmenhove et al., 2015;Vukčević, 2016;Vukčević et al., 2017).Larsen et al. (2019) followed a stream function wave for 100 wave periods in a cyclic domain.The study provided a detailed description of the OpenFOAM setup and aimed to establish a best practice for the interFoam solver.They showed that interFoam can propagate waves quite accurately over long distances using specific settings for the numerical schemes, which resulted in a numerical diffusive balance.However, even in the optimised simulation oscillations in the velocity field close the free surface are present.The study also highlighted that interFoam is very sensitive to the settings of the numerical schemes and the applied CFL number, which they recommend to be less than 0.05 to obtain reasonable predictions of the wave kinematics.
The solver interIsoFoam (Roenby et al., 2018) was based on the interFoam solver, where the algebraic advection algorithm Multidimensional Universal Limited Explicit Solver (MULES) (Ubbink, 1997;Rusche, 2002) was substituted with the geometrical advection algorithm isoAdvector developed by Roenby et al. (2016).isoAdvector has been further improved by Scheufler and Roenby (2018), who improved the prediction of the orientation of the iso surfaces on both structured and unstructured grids.The method was introduced in OpenFOAM-v2006 and we have implemented our work in OpenFOAM-v1912.Therefore, we have not been able to study the impact of the new interface reconstruction method in the present work.The improved interface reconstruction is however expected to improve the free surface advection and thereby also improve the overall performance of the DSD solvers.interIsoFoam has been used to propagate a stream function wave by Roenby et al. (2017Roenby et al. ( , 2018) ) and Larsen et al. (2019).The studies showed that isoAdvector improved the prediction of the surface elevation and eliminated the wiggly interface observed with interFoam.However, the velocity prediction by interIsoFoam has severe over-and underprediction, which is much larger than with interFoam.Roenby et al. (2017) used a fine resolution of 20 cells per wave height to get acceptable results, where significant errors could still be observed in both the surface elevation and velocity profile.Vukčević et al. (2017) showed that the continuous density field creates an imbalance in the momentum equation, which causes spurious velocities in absence of surface tension.The discontinuity of the density field was captured with a newly developed Ghost Fluid Method (GFM).Originally the GFM method was developed by Fedkiw et al. (1999), Fedkiw and Liu (2001) and Fedkiw (2001).The GFM method in combination with a novel level set method was shown to accurately propagate waves over long distances by Vukčević (2016), however no clear validation, e.g.velocity profiles, was given for the velocity field.
More recently, GFM was combined with isoAdvector by Vukčević et al. (2018) where they studied a progressive stream function wave train with intermediate water depth and moderate steepness.The study only presented results from a mesh with 26 cells per wave height.This is a relatively fine resolution compared to the literature referenced earlier, where resolutions from 2 to 30 cells per wave height have been seen.They compared different components of the stream function solution, and reported a first-order surface elevation error of 3.81% after eight wave periods and a first-order horizontal velocity error of 4.2% after seven wave periods.The GFM method have also been used for prediction of wave loads by Jasak et al. (2015), Gatin et al. (2018), Liu et al. (2020) and Gatin et al. (2020).More recently the GFM method have also been implemented by Peltonen et al. (2020), who suggests an alternative formulation for reconstructing the free surface position from the volume fraction field of water.Other related work on the GFM method includes Meyer et al. (2016), Wemmenhove et al. (2015), Queutey and Visonneau (2007), Huang et al. (2007), Desjardins et al. (2008) and Egan and Gibou (2020).
The original Marker and Cell method (MAC) and later Surface Markers methods (Harlow and Welch, 1965;Nichols and Hirt, 1975;Chen et al., 1991Chen et al., , 1995;;Christensen and Deigaard, 2001;Raad and Bidoae, 2005) only resolved the fluid domain, which was also the case for the original Volume of Fluid Method (VOF) by Hirt and Nichols (1981).Many of these methods were developed for structured grids.Nielsen (2003) and Christensen (2006) developed and used the CFD code NS3, which only solves the governing equations in the water region and models the behaviour at the interface.The interest in a single fluid free surface method for OpenFOAM was already mentioned by Jacobsen et al. (2012), because the air flow induced by the wave motion is irrelevant in many engineering applications.Jacobsen et al. (2012) observed that the CFL criterion is sometimes dominated by velocities in the air, which restricts the time step and thereby increases the simulation time.Furthermore, it is also mentioned that only solving the governing equations in the water region would reduce the cell count significantly and thereby substantially reduce the simulation time.
The Direct Surface Description (DSD) method has been inspired by Nielsen (2003) and Christensen (2006), however the methods differ in: (1) the discretisation of the Laplacian operator in the pressure equation, (2) how the continuity equation and momentum equation is coupled and (3) how the fields are extrapolated near the free surface.
The DSD method has been implemented in the general polyhedral FVM framework of OpenFOAM, whereas NS3 is limited structured hexahedral meshes.The general polyhedral discretisation in OpenFOAM limits the implicit discretisation to the two immediate neighbours of each cell face, which gives a compact stencil for the Laplacian operator.NS3 uses a larger, higher-order stencil for the Laplacian operator, which is valid on structured hexahedral meshes.Another difference between NS3 and the DSD method is the advection scheme, where NS3 applies the algebraic Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) scheme (Ubbink and Issa, 1999) and the DSD method uses isoAdvector.
The DSD method eliminates the air phase, which is treated like a void phase with constant pressure.In the present work simple cases with a well defined free surface and theoretical solution are considered for validation and comparison to the existing solvers.The ultimate goal is to use the DSD method to predict wave loads on structures from breaking wave impacts.During such simulations pockets of air may be occluded by the water phase, leading to a compression of the air phase and impulsive loading.It should be emphasised that this cannot be captured by the method in its present form, because the air phase is treated as a void phase with constant pressure.In order to capture these effects the method must be extended to simulate the flow in both the water and air phase.Furthermore, the assumption of incompressibility of the air phase needs to be relaxed meaning that the air phase should be modelled as an compressible phase.This is beyond the scope of the present paper, where the main focus is to solve the issue with the spurious velocities near the free surface.However, it is certainly an aspect which needs to be studied in the future.
The objective of the present work is to develop an efficient single phase method for free surface flows in OpenFOAM based on Christensen (2006) and Vukčević et al. (2017), which can accurately predict the free surface elevation, pressure fields and the wave kinematics without resolving the boundary layer near the free surface.Furthermore, we aim to study: (1) the influence of different temporal integration methods (2nd order Adams-Moulton, 2nd order explicit Runge-Kutta and 4th order explicit Runge-Kutta), (2) the convergence, (3) the influence of the CFL number and (4) the computational speed.

Mathematical model
This section introduces the mathematical background for the numerical methods presented.A nomenclature is provided in Appendix A.

Governing equations
The present study analyses a flow of two incompressible Newtonian fluids in a gravitational field, where the fluids are separated by a sharp interface.The density, , is eliminated from the general form of the continuity equation, when the fluids are incompressible and  is piece-wise constant.The continuity equation for an incompressible fluid is where  is the velocity vector.In our application for ocean waves the heavy phase is water with a density of   = 1000 kg∕m 3 and the light phase is air with a density of   = 1 kg∕m 3 .Conservation of momentum for each fluid is described by the incompressible Navier-Stokes equations: where  is time,  is the kinematic viscosity,  is the gravitational acceleration vector, and  is the total pressure.The gravity term  from Eq. ( 2) is reformulated to 1  ∇ ( ⋅ ) assuming a constant density field.The momentum equation is rewritten to The density is isolated in a single term accounting for both pressure and gravity.The expression inside the gradient on the right hand side of Eq. ( 3) is denoted  ℎ =  −  ⋅  and called the relative pressure.It is the total pressure minus the hydrostatic pressure.Sometimes this expression is also called the dynamic pressure, however this is only valid if  is determined with respect to the still water level.The final form of Eq. (2) becomes The equation for  ℎ is presented later in Section 3.1.4after the discretisation has been introduced.

Interface capturing
The interface between the two phases is captured with the VOF method originally introduced by Hirt and Nichols (1981).The volume fraction of the water phase is computed by Hence,  = 1 indicates that a cell is filled with water and  = 0 that a cell is filled with air.Cells with 0 <  < 1 represent the interface.The VOF advection equation for the water phase is derived from the continuity equation Eq. (1) and reads According to Hirt and Nichols (1981), a VOF method needs to provide: 1.A numerical description of the location and shape of the free surface 2.An algorithm for the time evolution of the free surface 3. A scheme for imposing the desired free surface boundary conditions in the discretisation of the governing equations interFoam is based on a diffusive interface description, which is not in agreement with the original VOF method.This is handled with a ''work-around'' solution introducing the compressibility algorithm MULES.The interIsoFoam method uses an advection algorithm which ensures that the volume fraction field remains sharp and well defined.However, neither interFoam nor interIsoFoam provide a scheme for imposing free surface boundary conditions, the third condition above.
The interIsoFoam solver is missing a scheme for imposing boundary conditions at the free surface.To implement boundary conditions in the discretisation of the governing equations the location of the free surface is needed.Cells with  > 0.5 are denoted wet and cells with  < 0.5 are denoted dry.The free surface location    along the vector   spanning from the face owner cell  to the face neighbour cell  is defined as where  is the normalised location of the free surface between  and , see defined  as whereas Christensen (2006) defined  as a conditional expression where   is the face normal vector and   is the face tangential vector.The gradient of the volume fraction normal and tangential to the cell face is defined as Fig. 1(a) illustrates the projections performed in Eq. (10).

Free surface boundary conditions
At the interface between two phases the kinematic and dynamic free surface boundary conditions must be considered.The kinematic boundary condition (Batchelor, 1967) states that the velocity field is continuous across the free surface: Here the light phase is air and the heavy phase is water.The boundary layer around the free surface due to the kinematic condition is difficult to resolve in numerical simulations of ocean waves due to the needed spatial resolution.The kinematic boundary condition also makes the continuity equation valid for two incompressible phases.The dynamic boundary condition follows from momentum conservation, where the forces acting on the fluid at the interface are in equilibrium.Following Tuković and Jasak (2012) the tangential stress balance is where    is the unit normal vector on the interface pointing from heavy towards light fluid, ∇  = ∇ −       ⋅ ∇ is the surface gradient operator,  ,  = ( −       ) ⋅    is the tangential velocity component,  is the identity tensor and  is the surface tension coefficient.The normal stress balance yields a relative pressure jump across the free surface that reads where  = −∇  ⋅    is twice the mean curvature of the interface.
The first term models a pressure jump due to surface tension.The second term −2(  −   )∇  ⋅    models a pressure jump due to the normal viscous force at the interface expressed in terms of the surface divergence of the interface velocity according to Chen et al. (2000).
The last term (  −   ) ⋅    models a pressure jump caused by the density jump at the free surface.The free surface boundary conditions are simplified through assumptions (Vukčević et al., 2017).The tangential stress balance in Eq. ( 12) is neglected, as it is of minor importance for free surface flows at high Reynolds numbers (Huang et al., 2007).The normal stress balance in Eq. ( 13) simplifies to assuming a continuous viscosity field and neglecting the effect of surface tension, which is valid for high Weber number flows.The simplified free surface conditions yields an additional condition, which follows from inspection of the terms in the momentum equation in  Eq. ( 4) (Vukčević, 2016).Since the temporal, convective and diffusive terms are continuous, the pressure gradient term on the right hand side must also be continuous: As the DSD method only models the flow in the water, we need to reformulate the boundary conditions at the free surface.The dynamic boundary condition from Eq. ( 14) is reformulated and total pressure in the air is approximated with a constant total pressure:  = 0.The flow is driven by pressure gradients, and assuming  = 0 in the air  The total pressure at the free surface on the air side is set to     = 0.The final relative pressure boundary condition is then which defines an inhomogeneous Dirichlet condition that must be enforced at the location of the free surface.It is not possible to derive a free surface boundary condition for the velocity at the free surface when only the water phase is considered.The movement of the free surface should correspond to the velocity at the free surface, but that does not provide a boundary condition.Instead we use linear spatial extrapolation from internal cell centres to compute velocities at face and cell centres outside the water domain needed in the discretisation.Fig. 2 illustrates the definition of a cell centre and a face centre and the extrapolation procedure is described in Section 3.2.4.The air faces of the surface cells receives an extrapolated velocity, which is used to evaluate the divergence term on right hand side of the pressure equation, which is presented later in Section 3.1.4.The air cell faces are those cell faces that have an air cell as its face neighbour.The predicted pressure is thereby directly affected by the extrapolated velocities through the continuity constraint.

Numerical method
During our review of the literature around OpenFOAM, we have noticed that the discretisation of the governing equations is often presented in a very compact form.For the sake of completeness, our description includes a broader introduction to the discretisation in OpenFOAM followed with specific details on our development.The first part describes the fundamental interpolation and discretisation of the momentum equations, the continuity equation and the derivation of the pressure equation.The numerical discretisation is presented with the implicit time integration used in OpenFOAM.Especially the details of the derivation of the pressure equation is often hidden in the compact notation.The knowledge about these details is necessary to understand the underlying numerical method and the developed new method.The second part is dedicated to aspects of the numerical method and discretisation of the free surface boundary specific to the DSD method.The first section describes the discretisation of the pressure gradient at the free surface.The second section describes the movement of the free surface.The third section defines the extrapolation practice used at the free surface boundary in the DSD method.The fourth section describes how the least squares gradient calculation has been modified in OpenFOAM for cells at the free surface.Lastly, the fifth section describes the explicit Runge-Kutta time integration of the momentum and pressure equations.Flowcharts for the solvers dsdFoam, dsdrkFoam and interIsoFoam can be found in Appendix B.

Interpolation
In OpenFOAM linear interpolation determines face values from cell values by where  is a general dependent variable and the overbar indicates linear interpolation.assumed to intersect the face centre.It is possible to apply a skewness correction, 1 when ⃖⃖⃖⃖⃖⃖ ⃗   does not intersect the face centre, but it is not used in this work.

Momentum equation
The time integration in OpenFOAM is performed according to Section 6.2.4 (Ferziger and Perić, 2002).The momentum equation in Eq. ( 4) is integrated over a time step  centred around : The integral is approximated with the mid-point rule, which yields The equation is divided by  and integrated over a control volume  : The mesh is assumed to be static in the present derivation.
Temporal term.The implicit time derivative is discretised by secondorder Adams-Moulton scheme. 2 OpenFOAM uses a finite difference formulation of the temporal discretisation, which is described in Section 13.4.2(Moukalled et al., 2016).The velocity at − and −− 0 is denoted  0 and  00 .The scheme is derived by expressing  0 and  00 with a Taylor series expansion around .The discretised form of the acceleration term for non-uniform time steps reads where the algebraic coefficients are defined as The first time step in a simulation is discretised with the first-order forward Euler scheme.
1 Specified by the keyword skewCorrected in OpenFOAM.
2 Known by the keyword backward in OpenFOAM.
Convection term.The volume integral of the convection term is transformed to a surface integral with Gauss' theorem.The surface integral is approximated with a summation over the cell faces, where the face integral is approximated with the mid-point rule: The summation over cell faces is split in two sums for the owned and neighbouring cell faces to account for the face normal vector always pointing from owner cell P towards neighbour cell N. Hence the above cell face summation assumes the following split, which is valid for all face summations.
The volumetric flux defined at the cell faces is defined as  *  =   ⋅  *  .The surface area normal vector is given by   =     , where   is the face area.When the convection term is treated implicitly, the term is linearised by computing  *  explicitly from the previous time step or iteration.This is indicated with superscript * .The nonlinearity is resolved by the outer iteration loop in the PIMPLE algorithm in OpenFOAM.The PIMPLE algorithm consists of an outer and inner loop.The outer loop is a series of Picard iterations and the inner loop is the Pressure-Implicit with Splitting of Operators (PISO) algorithm (Issa, 1986;Tuković and Jasak, 2012).
Diffusion term.The discretisation of the diffusion term in Eq. ( 26) follows Jasak (1996) and Muzaferija (1994). is the surface of the control volume.
The surface normal gradient   ⋅ (∇)  is approximated by where the orthogonal face area vector is computed as . This is called the over-relaxed approach (Jasak, 1996).The magnitude of the surface area vector   is denoted by two vertical lines as |  |.The operation is defined as Pressure gradient term.The volume integral of the pressure gradient term is transformed to a surface integral with Gauss' theorem.The surface integral is approximated by a sum over the cell faces, where the integral of each face is given by the mid-point approximation: When the pressure gradient term is left undiscretised in the momentum equations the volume integral yields

Continuity equation
The continuity equation from Eq. ( 1) is integrated in time and the integral is approximated using the mid-point rule: Eq. ( 30) is integrated over  via Gauss' theorem as a sum over the cell faces, and the integral over each cell face is approximated using the mid-point rule:

Pressure equation
The PISO algorithm implemented in OpenFOAM is described by Tuković and Jasak (2012), Tuković et al. (2018) and Uroić (2019).The following derivation is presented for static meshes.The pressure equation is derived from the discretised continuity equation Eq. ( 31) and a partially discretised momentum equation, where the pressure gradient is not discretised.This form of the momentum equation is commonly known as the semi-discretised form in the literature.The semidiscretised momentum equation presented in Eq. ( 32) is derived from Eq. ( 21) after substituting the discretisation presented in Eqs. ( 18), ( 22), ( 24), ( 27) and ( 29).Finally the semi-discretised momentum equation is divided by  to normalise the equation.This leads to where ∑  is the off-diagonal contributions from the neighbouring cells.The diagonal matrix coefficient for cell  (    ) corresponds to  in Eq. ( 23) and    is The off-diagonal matrix coefficient from each cell neighbour  is given by and the source term is given by The cell centre velocity,   , is expressed from the momentum equation in Eq. ( 32), which leads to The essence of the momentum interpolation method is to express the momentum equation at the cell faces to mimic a staggered grid.This is obtained by interpolating the collocated grid coefficients in Eq. ( 36) to the cell faces.The resulting expression for the cell face velocity    using the momentum interpolation method is The above expression is needed in Eq. ( 31), the discretised continuity equation.The terms ( )  are found from linear interpolation of neighbouring cells to cell faces indicated by the overbar.The pressure gradient term The face velocity of the previous time steps  0  and  00  are described by Tuković and Jasak (2012), which explains why small time steps may lead to pressure oscillations and how the solution of the pressure equation is made independent of the time step size using a method developed by Yu et al. (2002).The cell face flux   satisfies the discretised continuity equation, whereas the cell centre velocity interpolated to the face centres (   )  does not satisfy the discretised continuity equation.The idea is to correct by the difference between the face normal where   is the face normal vector and   is the face area.The pressure equation is obtained by substituting the face velocity   from Eq. ( 37) in the discretised continuity equation Eq. ( 31): After the pressure equation is solved, the volumetric face flux is computed by and the cell centre velocity is computed by the expression in Eq. ( 36).
To the authors impression, understanding the connection between the theoretical background and the source code in relation to the pressure-velocity coupling is a challenge for a majority of the Open-FOAM users.We hope that linking the presented theory to the variable names used in OpenFOAM will be useful especially for those who are new to or not familiar with the OpenFOAM source code.To highlight the connection between the presented theory and the implementation in OpenFOAM-v1912,  0  and  00  from Eq. ( 38) are inserted in Eq. ( 40), which yields The former equation is rearranged to obtain: in the third line is a correction to the flux from the discretisation of acceleration in the second line, which ensures that the solution is independent of time step size and eliminates pressure oscillations in the solution for very small time steps (Tuković and Jasak, 2012).
The fifth line describes the flux from the pressure gradient.Since   ⋅   =   the third and fourth line can be further simplified and the expression is rearranged to obtain as given in Box I.The OpenFOAM code terminology is indicated with the braces, where phiHbyA is the uncorrected volumetric flux computed by linear interpolation of HbyA in Eq. ( 36) and rAUf* fvc::ddtCorr(U, phi) is the correction from the source terms of the discretised acceleration term.

Discretisation of cell centred pressure gradient
The cell-centred gradient in cell  from Eq. ( 28) forms a sum over values at the cell faces, which is obtained by linear interpolation: As previously mentioned the sum needs to be separated in two summations, due to the sign convention of the surface area normal vector   .The sum ∑  includes all the cell faces where  is the owner of cell face  .Likewise, ∑ ℎ includes all the cell faces where  is the neighbour of cell face  .The difference in sign comes from the positive sign convention of the surface area normal vector   from owner cell  towards neighbour cell .However, the Gauss discretisation of the cell gradient assumes that the cell face normal vectors point away from the discretised cell .Hence, the contribution from face  to the discretised gradient in  is positive, whereas the contribution from face  to the discretised gradient in  is negative.
The discretisation includes the pressure boundary condition at the free surface using extrapolated pressures in the neighbouring air cells.When the owner cell  is water the extrapolated value   ℎ, at neighbour cell  is given by linear extrapolation: The contribution from face  to the computed gradient in cell  can be written as follows after inserting  ℎ, =   ℎ, in Eq. ( 44): When the neighbour cell  is water the extrapolated value   ℎ, at owner cell  is given by: The contribution from face  to the computed gradient in cell  can be written as follows after inserting  ℎ, =   ℎ, in Eq. ( 44): The negative sign comes from the aforementioned sign convention of the surface area normal vector   from  towards .

Discretisation of face normal pressure gradient
The pressure gradient term from Eq. ( 39) at a fully submerged face  is discretised with the same scheme used for the diffusion term in Eq. ( 27).The non-orthogonal correction term is omitted in the following derivation, since the only change is the application of the pressure gradient scheme from the previous section.The orthogonal part of the face normal pressure gradient is: At free surface faces where cell  is wet and cell  is dry, the Dirichlet condition for the relative pressure from Eq. ( 17) gives the following expression for the pressure gradient term: Similarly when cell  is dry and cell  is wet, the Dirichlet condition  ℎ,  is implemented by the expression: The density at the cell face is that of the wet phase, i.e   =   .

Free surface advection
The governing advection equation presented in Eq. ( 6) is solved using the isoAdvector algorithm (Roenby et al., 2016).It geometrically captures the change of the VOF field  over a time step  based on a fixed velocity translation  of an iso-surface reconstructed from the volume fraction field .In brief, the isoAdvector algorithm updates the volume fraction field  by Eq. ( 52).
is a list of the cell faces for a single cell.  is the total volume of water flowing from the current cell during one time step  through face  into the neighbouring cell.  is approximated from ,  and  in Eq. ( 53), where   is the volumetric cell face flux through face ,   is the surface area vector at face  and   () is the cell face area of face  that is submerged in water as function of time.
Furthermore, it is assumed that the iso-surface, reconstructed from , moves with a constant normal velocity that is representative for the movement of the iso-surface during the time step from   to  +1 .
When isoAdvector is used in combination with an implicit time integration the solver consists of two loops: An outer and an inner loop.The outer loop iterates over the movement of the free surface and the non-linearity in the convection term in the momentum equations.The inner loop is an implementation of the PISO algorithm for a given number of correction steps.The implicit solver employs an explicit Euler extrapolation in time for the velocity field used to advect the free surface, when the solver passes the first outer corrector.If only a single outer corrector is used the advection of the free surface is expected to be first-order accurate.If the solver employs more than a single outer corrector the velocity field advecting the free surface is found by a central difference approximation between the previous velocity field and the newest estimate of the velocity field to the new time step.This should lead to second-order accuracy.For the implicit method DSD solver and the explicit DSD solvers ,  and  used to advect the free surface are estimated from the previous time instance in the first outer iteration: In all the subsequent outer iterations in the implicit DSD solver  and  are estimated with a central difference approximation, and  is equal to the newest prediction from the previous outer iteration: After the prediction of the new  field in each outer iteration, the velocity and volumetric flux fields are reset to the value of the previous outer iteration: The DSD Runge-Kutta solvers use an explicit first-order Adams-Bashforth prediction of the input parameters  and  to isoAdvector.

Extrapolation
Whenever values are required in face centres or cell centres at the other side of the free surface in the air region, a linear extrapolation based on known cell centre values and gradients is performed.The extrapolation is given in Eq. ( 57).Here subscript  indicates the target coordinate point of the extrapolation.The subscript  indicates the points, that the extrapolation is based on.It is assumed that the gradient is known at the points , where the extrapolation is performed from.
The expression is a weighted linear extrapolation based on the field value and field gradient in the selected extrapolation cells .The linear extrapolation performed from each extrapolation point is given by the expression . The weight 1∕|  −   | increases the importance of the extrapolation points close to the target point and decreases the importance of the extrapolation point far away from the target point.The extrapolation is combined with different conditions for the selection of the points included in  to yield different extrapolation practices.These conditions are constructed to ensure that extrapolation is only performed from cells in the water region.The selection of extrapolation points  can be divided in four cases: 1. Extrapolation to new water cell centres, which have entered the water region during the advection step, where the cell shares at least one cell face with the water region at the former time instance.2. Extrapolation to new water cell centres, which have entered the water region during the advection step, where the cell only shares an edge line or a vertex point with the water region at the former time instance.3. Extrapolation to new water cell faces, which have entered the water region during the advection step.4. Extrapolation of fields at previous time instances used in the discretisation of ∕.
Fig. 3 illustrates the first three extrapolation cases (1-3) from the list above.The first case represents the case where an air cell has entered the fluid region during the advection of the free surface.The target  1 is the cell centre of the new fluid cell and the extrapolated value is based on the field value and the field in gradient the cells marked with 1.These cells shares a cell face with the new water cell  1.
The second case also represents the case where an air cell has entered the fluid region during the advection of the free surface.The target  2 is the cell centre of the new fluid cell and the extrapolated value is based on the field value and the field gradient in the cell marked with 2.This cell only shares a vertex with the new water cell  2.
The third case is an example of the extrapolation procedure, when new face values are needed in new water cells.In this case cells marked by 3 are used to estimate the value at the new cell face marked by  3. The fourth extrapolation case is used in the discretisation of the acceleration term in the momentum equation, when cells with no previous time history are encountered.The velocity field at the previous time steps are extrapolated according to the new location of the free surface.This means that extrapolated values will be assigned to cells, which have changed from air to water during the time interval from the previous time instance to the new time instance.In relation to this it is assumed that cells, which have entered the water region, are limited to a region of a single layer of cells from the old time to the new time.The single layer of cells is illustrated in Fig. 3, where the green cells mark the single layer next to the fluid region marked with blue cells.This restricts the CFL number to   < 0.5, however this is not considered a problem, since decreasing the CFL number is a more effective way to get more accurate results compared to increasing the mesh resolution (Larsen et al., 2019).

Gradient calculation at the free surface
The least squares gradient calculation in OpenFOAM is programmed to include all cell neighbours in the calculation (Jasak and Weller, 2000).However, the DSD method does not compute any values in the air, hence the dry neighbouring cells next to the free surface should not contribute to the gradient calculation.Therefore, we have modified the least squares gradient calculation in order to only account for the neighbouring cells categorised as water.
The normal least squares gradient is determined by minimising the error of a plane fitted to a cloud of points, which include cell  and its immediate neighbours .The velocity gradient in cell  is computed by where  is a symmetric 3 × 3 matrix for the three dimensional case.
The  × 3 matrices   and   , the  × 3 distance matrix  and the  ×  weighting matrix  are defined as where  is the number of neighbour cells, ⃖ ⃗   is the 1 × 3 velocity vector in cell centre , and ⃖⃖⃖⃖⃖⃖ ⃗    is the 1 × 3 vector spanning from cell centre  to cell neighbour centre   .
The idea behind this scheme is to evaluate the gradient with the standard least squares method in all cells and then correct the estimated gradients in the free surface cells.The free surface cells are the cells containing one or more free surface cell faces.This is illustrated in Fig. 1(b), where the cell faces along the red line are free surface faces.
For each free surface cell in the water the immediate neighbours inside the water region are identified.A distance weighted average of the neighbouring gradients is computed and assigned to the free surface cell using the extrapolation procedure presented in Section 3.2.4.If cell  has one or more fully submerged cells, i.e. cells with no free surface faces, then these cells will be used to extrapolate from.This is extrapolation case 1 in Fig. 3 where the cells marked with 1 are used to estimate value in  1.For some surface cells there will be no immediate neighbour cells fully submerged in water.This is extrapolation case 2 in Fig. 3, where the cell marked with 2 is used to estimate value in  2. The extrapolation of the gradient fields corresponds to a weighted zero gradient condition on the velocity gradient field.In other words, the velocity field is expected to continue to change in the same way in cell P as it does in the immediate neighbouring cells in water.The scheme algorithm can be summarised in the following steps: 1. Compute standard least squares gradient in all cells.

Table 1
Butcher tableau for n-stages, RK2 and RK4. n-stages 2. Extrapolate ∇ according to extrapolation procedure from Section 3.2.4,where  is simply replaced by ∇. 3. Set gradient to zero in all air cells.This scheme has proven to be very robust when applied to complex free surface flow simulations.The scheme extracts as much information as possible from the nearest cells, when limiting the search radius to one layer of cells from the target cells.

Explicit Runge-Kutta time integration
The explicit Runge-Kutta scheme can be defined by a Butcher tableau, which enables us to write an algorithm that can handle any explicit Runge-Kutta time integration.The Butcher tableau for an explicit Runge-Kutta scheme with  stages, the RK2 scheme and the RK4 scheme are given in Table 1.
Each row in the table represents a single stage in the Runge-Kutta method.After evaluating all stages the final prediction is made from a linear combination of the  stages.For explicit Runge-Kutta methods the diagonal elements are zero, which means that the first stage simply evaluates the function  (  +  1 ,  1 ) =  (  ,   ) from the solution at the last time step.The coefficients  and  are defined according to the chosen Runge-Kutta scheme, while the coefficients  are computed by where  is the row index and  is the column index in the Butcher tableau.
The implementation of the explicit Runge-Kutta time integration together with an equation for the pressure follows (Sanderse and Koren, 2012).We solve an equation of the form The acceleration term is discretised with the forward Euler scheme, and the pressure gradient term is left undiscretised to form an equation for the pressure solved at each Runge-Kutta stage.After the pressure has been found its contribution is added to the intermediate velocity and volumetric flux at the end of each Runge-Kutta stage.The function  (, ) is evaluated from the solution at the previous Runge-Kutta stages or previous time step.The general form of the equation at stage  is The unknown velocity is isolated, leading to The intermediate velocity, ũ , is corrected by the pressure gradient after the pressure equation has been solved.The continuity equation must be satisfied at each stage, which yields the pressure equation: The constant   works like a Lagrangian multiplier, which can be neglected without affecting the result, as long as it is also removed from the pressure gradient updating the velocity (Sanderse and Koren, 2012).
After the Runge-Kutta stages have been completed, the final prediction of the velocity to the new time step is found by The pressure equation derived from Eq. ( 65) is After the pressure has been found, the velocity and flux are corrected to give the final prediction to the new time  +1 .The discretisation of the pressure equation follows the same discretisation procedures presented in Sections 3.1.4and 3.2.2.

Results and discussion
The DSD method has been implemented in two solvers.The first solver is an implicit solver dsdFoam using the second-order Adams-Moulton scheme to discretise the acceleration.The second solver is an explicit solver dsdrkFoam with a Butcher tableau based explicit Runge-Kutta time integration of the pressure-velocity coupling and a first-order Adams-Bashforth prediction of the input variables to the free surface advection algorithm isoAdvector.The solver dsdrkFoam is tested with a second and fourth order Runge-Kutta scheme.The second order scheme is Heun's method and the short name used throughout the paper is DSD-RK2.The fourth order scheme is the low cost RK4 scheme and will be known by the name DSD-RK4.The short name for the implicit solver   is DSD.This leads to the three different numerical methods which are compared to interFoam and interIsoFoam.
The solvers are used to simulate two 2D test cases: (1) A still water level and (2) A progressive wave train described by stream function theory.

Still water level test
It may seem trivial to simulate a still water level, but the case is in fact quite challenging since it exposes any imbalances in the numerical model.The test case is a still water level in an inclined square of 1 × 1 m with 50 × 50 cells.We simulate 1 s using a constant time step of  = 0.001 s.The case is simulated with interFoam (IF), interIsoFoam (IIF), DSD, DSD-RK2 and DSD-RK4.
Before we developed the DSD method the starting point for our research was the GFM method (Vukčević et al., 2017).The solver GFM is our implementation of the GFM method in OpenFOAM combined with isoAdvector (Vukčević et al., 2018).The only modification is the formulation of the inverse distance from Eq. ( 9).The solver GFM-RK4 is our implementation of the GFM method combined with an explicit fourth-order Runge-Kutta time integration for the pressure-velocity coupling.
Fig. 4 presents the simulated velocity magnitude fields at  = 1 s for the seven solution algorithms.The free surface contour is given by a red line.The DSD and GFM methods provide very low velocities caused by the discrepancies between the true position of the free surface and the reconstructed position from the volume fraction field.When the inverse distance is specified based on the known location of the free surface the resulting velocities are of the order 10 −11 m∕s for the DSD solvers and 10 −12 m∕s for the GFM solvers.The maximum velocity in the air is 4.32 m∕s for interIsoFoam and 1.81 m∕s for interFoam.
The air region from the interFoam and interIsoFoam simulations have been clipped out to highlight the effect in the water cells.The imbalance is caused by the linear interpolation of the density field in interIsoFoam and interFoam (Vukčević et al., 2017).However, after the air has been accelerated the linear interpolation of the velocity field accelerates the momentum transfer to the water column through the convection term, as the boundary layer is not resolved.This test shows a general problem with the interFoam and interIsoFoam methodologies.

Progressive surface gravity wave based on stream function wave theory
The stream function wave test case from Vukčević et al. ( 2018) is adopted to enable a qualitative comparison to the performance of their model, which combines the GFM method with isoAdvector.
The domain is 13 corresponding to 70.32 m and the height is 2 m.The water depth is ℎ = 1 m and the still water level is located at  = 0 m.The regular wave train is initialised with the stream function wave theory (Fenton, 1988), generated by the waves2Foam toolbox (Jacobsen et al., 2012).The first and last 1.5 are relaxation zones, where the theoretical solution is gradually blended with the numerical solution at the inlet and vice versa at the outlet.
The wave height is  = 0.3 m, the wave period  = 2 s and the initial phase shift is  =  rad.The stream function wave is based on a zero Eulerian time mean current.The corresponding OpenFOAM settings are specifyEuler true; and eulerVelocity 0;.The gravitational constant is  = 9.81 m∕s 2 .The stream function theory is based on 32 coefficients and the solution is found with 10 iterations.
The derived wave properties are the wave length  = 5.409 m, the wave number  = 1.162 m −1 , the angular frequency  = 3.142 s −1 and the wave celerity  = 2.704 m∕s.The steepness of the wave is moderate with a value of ∕ = 0.05546.The relative depth of ℎ = 1.162 tells that the wave is travelling on an intermediate water depth between the shallow and deep water limits: 0.25 < ℎ < 3.14.
The case is simulated on five different meshes with a constant time step that decreases when the mesh is refined to keep the initial CFL number constant at approximately CFL = 0.09.The solvers are capable of handling variable time steps.The time step is kept constant to enable direct comparison between simulations from different solvers.The simulated duration of time is 40 s corresponding to 20 wave periods.The properties of each mesh case is given in Table 2, which also gives the resolution of wave with the number of cells per wave height ∕ and the number of cells per wave length ∕.The cell aspect ratio is  = 1.An overview of the applied boundary conditions is provided in Table 3.
The acceleration term ∕ can be discretised with different discretisation schemes.For DSD and interIsoFoam, it is discretised with the second-order Adams-Moulton scheme.The second-order Adams-Moulton scheme cannot be selected in OpenFOAM when using the The boundary condition on the top patch for velocity and pressure are not used in the solvers that only consider the water region, instead we apply the free surface boundary conditions at location of the free surface.b There is no boundary condition on the fluid volume fraction field alpha.water,which is defined in the entire domain.
interFoam.As an alternative a 50∕50% blend of the implicit Euler and Crank-Nicolson schemes is employed.The DSD-RK2 and DSD-RK4 employs a different time integration method using an explicit Runge-Kutta scheme, which is described in Section 3.2.6.The OpenFOAM case settings for interFoam and interIsoFoam can be found at the git repository: https://gitlab.gbar.dtu.dk/jrkp/Article1.git.

Surface elevation
Fig. 5 presents the surface elevation centred around the last wave crest identified within  = 10−11 at ∕ = 20.The focus is the shape of the wave, whereas the wave phase error is examined in Section 4.2.4.The wave is propagating from left to right.The black line is the numerical surface elevation and the red line is the theoretical surface elevation.Selected subplots in the figure contain a 300% zoom on the wave crest in a small box.The DSD solvers are better at maintaining the wave height and form compared to interFoam and interIsoFoam.
Local undulations on the free surface are observed in the solution by interFoam.The undulations and the wave length of the undulations decrease in magnitude when the mesh is refined.The cause for the undulations seems to be the MULES advection algorithm, since the undulations cannot be observed in the surface elevation predicted by interIsoFoam.The only difference between the two solvers is the advection algorithm.
The surface elevation predicted by the DSD solvers is smooth without oscillations.It is observed that the explicit DSD solvers DSD-RK2 and DSD-RK4 are slightly more accurate than the implicit solver DSD.The two explicit Runge-Kutta solvers show almost identical results.However, later in Section 4.2.4 a convergence study of the errors is presented, which reveals that DSD-RK4 is the most accurate solver with respect to the surface elevation.
Fig. 6 presents the surface elevation in the interval 9 − 11 at ∕ = 20 for all cases and solvers.The last cycle of identified crests and troughs have been plotted with dots in colours matching the legend.The theoretical crest and trough envelope are shown as black lines and the theoretical crest and trough positions are indicated with red and blue circles.This qualitative comparison clearly illustrates the deficiency of interFoam.It consistently produces a wiggly surface elevation in the cases ∕ = [3, 6, 12, 15], and halves the wave height for ∕ = 1.5.The coarse cases illustrate that the DSD solvers are better at maintaining the wave height also on coarse meshes.In terms of phase shift interFoam and interIsoFoam are also outperformed by the solvers using the DSD method.
In case ∕ = 6 for the interIsoFoam solver the wave crest envelope is observed to oscillate.We have studied this further and found that all cases and solvers have some degree of wave reflection with a period of  ∕2.Fig. 7 shows the envelope of the surface elevation with zoomed plots of the crest region at the top and the trough region at the bottom.The oscillations at the crest and trough are in phase with each other and they have a wave length of ∕2.This shows that the phenomenon is a first harmonic reflection.The magnitude of the oscillations in the wave height is around 0.33% of the specified wave height.This is well below 1%, which is the reported amount of reflection generated by waves2Foam for a relaxation zone of 1.5 (Jacobsen et al., 2012).

Velocity fields
The velocity field given by the stream function wave theory is shown in Fig. 8. Figs.9-13 presents the 2D velocity field from case ∕ = 6 at  = [10.5− 11.5] and ∕ = 19.5 for the five solution methods.
The velocity field predicted by interFoam appears acceptable in the water.Still, the velocity field in the air shows that interFoam has some fundamental instability issues when simulating surface gravity waves.As illustrated with the still water level test, the air velocities can significantly effect the fluid velocities near the free surface during a simulation over a longer period.The velocity field in the water predicted by interIsoFoam is clearly faster in the upper region of the wave compared to the theoretical solution.Furthermore, comparing the velocity vectors to interFoam they are pointing more forwards at the front of the wave.Very large velocities can be observed in the air showing that accurate advection of the free surface alone cannot eliminate the problems in the air.On the contrary, the problem seems to have worsened.A possible explanation for the intensified velocities is as follows: The geometrical advection algorithm does not smear the volume fraction field as the MULES algorithm.The combination of a sharp advection algorithm with a continuous pressure and velocity field over the free surface seems to make the numerical solution more unstable and incorrect.The reason may be that the two modelling philosophies are not consistent, one part of the method in interIsoFoam keeps the surface sharp, whereas the other tends to smear the solution.
The velocities in the air is obviously not a problem for the DSD solvers, where the air phase is eliminated.For the relatively coarse mesh with 6 cells per wave height small differences between the DSD solvers are observed.

Velocity profiles under the wave crest
From the references in the introduction we have observed that OpenFOAM solvers are often validated against the surface elevation and less or even no effort is spent on validating velocity fields.The wave kinematics are in many cases very important to derived analyses, e.g.wave forces and wave run-up on structures.
Fig. 14 presents the horizontal velocity profile under the last wave crest identified in  = 10−11 at ∕ = 20.The velocity data are taken from the column of cells closest to the wave crest.The wave crest has been estimated by a parabolic interpolation of the surface elevation.The offset between the location of wave crest and the position of the column of cells with velocity data is accounted for in the theoretical velocity profile in order to match the wave phase of the theoretical profile to the numerical profile.The small boxes show a 250% zoom of the velocity profile just below the free surface.The interFoam solver underpredicts the velocity over the entire water column for the coarse cases, which is also expected since the wave height is significantly reduced.Near the wave crest interFoam has a tendency of overshooting at the crest and underestimating in the region just below the crest.Furthermore, the velocity profile contains small oscillations, which may be caused by the wiggly surface elevation.interIsoFoam overestimates the crest velocity by 20 − 30%, which is maintained with a finer resolution.Just below the overestimation the velocity is underestimated.
The DSD solvers have a tendency to slightly overpredict the velocity near the wave crest.The velocity profile is also underpredicted over the entire depth for coarsest resolution, which is linked to the wave height decrease.The DSD solvers have a tendency to overshoot the velocity at the top and undershoot the velocity at the bottom.The deviations quickly reduces as the mesh is refined and there are almost no visible deviations once the mesh resolution reaches ∕ = 12.Later in Section 4.2.4 the velocity errors are studied in more detail.The small boxes show a 400% zoom of the velocity profile just below the free surface.interFoam and interIsoFoam underpredict the vertical velocity and the velocity profile suffers from oscillations, which attenuates as the mesh is refined.It can be noted from the zoomed plots of the velocity profile near the free surface, that interFoam and interIsoFoam also produce kinks in the profile for the fine mesh resolutions.The zoom on the DSD solvers emphasises how well the solution agree with the theory.To conclude the DSD solvers provide more accurate velocity profiles without unphysical oscillations compared to interFoam and interIsoFoam.Note that a resolution of 15 cells per wave height is in fact only a medium fine resolution with respect to the range used in the literature.

Simulation time and convergence
This section presents a study of the solver performance with respect to accuracy and simulation time.The mean and maximum error of the velocity magnitude profile (((||)), ((||))), the wave height error (((||))) and the wave phase error (((| ℎ  ∕|))) are included in the convergence study.
Data is extracted below the wave crest in the interval  = 10 − 11 during the last 5 s of each simulation.Fig. 16a presents the convergence of the maximum error of the velocity magnitude profile.Fig. 16b presents the convergence of the mean error of the velocity magnitude profile.
Fig. 17a presents the convergence of the mean error of the wave height.Fig. 17b presents the convergence of the mean error of the phase shift.The convergence rate of the interFoam solver is unsteady and in a single case the maximum error of the velocity magnitude profile increases from the second last to the last case.The expected behaviour is a constantly decreasing error.It also seems that there is a tendency for the convergence rate of interFoam to flatten out as the mesh resolution is increased, however a larger study on several different test cases is required before a conclusion can be made.
For interIsoFoam the convergence rate of the maximum velocity magnitude profile error attenuates and reaches an almost constant error, which is undesirable.However, the mean error of the velocity and wave height decreases with a fairly constant rate between −1 and −2.Since the Finite Volume method is based on linear interpolation, it is expected that the error should decrease with a slope of −2 in a double logarithmic plot.
The DSD solvers show a convergence rate closer to −1 than −2 for the maximum velocity error.However, for the mean velocity, wave height and phase shift the convergence rate is −2 and even closing in on −3 for the wave height error with DSD and DSD-RK2.This is better than what could have been expected for the explicit DSD solver, since it uses a first order extrapolation in time for the input variables to isoAdvector.The explicit and implicit DSD solvers are almost equally accurate across the presented mesh resolutions, and they are significantly more accurate than interFoam and interIsoFoam.
For the finest mesh case the wave height error is 18 times lower with DSD-RK2 compared to interFoam for the same simulation.In addition the DSD-RK2 solver quarters the CPU time for the current case in comparison to interFoam as seen from Fig. 18.The phase shift error in Fig. 17b is expressed as a fraction of the wave length.The interFoam solver gives the lowest phase error at the second coarsest mesh, but then the error increases again with further mesh refinement.
The interIsoFoam solver got the largest error among all the solvers and the error converges with a slope around −1, which indicates first order accuracy.The DSD solvers make an accurate prediction for the coarsest case.This is interpreted as a coincidence.The error increases in the second coarsest mesh and then starts to decrease with constant rate with a slope around −2.This indicates that the DSD solvers are second order accurate.Fig. 18 presents the computational time normalised by the computational times from DSD-RK2.In terms of computational   speed interFoam and interIsoFoam are slightly faster on the coarse meshes.However, as the mesh is refined the simulation time relative to DSD-RK2 increases and for the finest mesh interFoam and interIsoFoam are around four times slower than DSD-RK2.The DSD solver gradually performs better as the mesh resolution increases.The DSD-RK4 is roughly 1.5 times slower than DSD-RK2 across all cases.

Effect of CFL number varying the time step
This section presents the effect of the CFL number on the accuracy.Running simulations with higher CFL numbers through a larger time step leads to lower computational cost.However, as the CFL number increases, the errors of the temporal discretisation may increase and change convergence or the numerical method may become unstable.The CFL number is varied by changing the time step size keeping everything else constant.Fig. 19 presents the variation of the maximum velocity error, mean velocity error, wave height error and phase shift error as function of the CFL number.The interIsoFoam solver crashed for the case with an initial   = 0.18, because the velocities had increased enough to make the solution unstable.The interFoam is more stable.It was able run the simulations up to   = 0.38.It should be stressed that the simulations only crashed because a fixed time step was used.It is also observed that the error levels only gradually increases, which indicates that the cause of the crash can be found in the air phase region of the simulation.The implicit DSD solver was able to run up to   = 0.3, where the error suddenly increases very rapidly.The explicit DSD solvers are capable of producing results for all the   numbers up to   = 0.45, however the error increases significantly after   = 0.23.The accuracy of both the implicit and explicit DSD solvers are generally observed to be several times better than interFoam and interIsoFoam up to   = 0.23.
The only exception being the implicit DSD solver at   = 0.23 for the maximum velocity error in Fig. 19a and the phase shift error in Fig. 19d.This suggests that running simulations with a CFL number up to 0.2 is an acceptable choice for the DSD solvers.
Decreasing the CFL number generally reduces the errors for the DSD solvers.The only exception being DSD-RK4 and DSD-RK2 for the maximum velocity error and DSD-RK4 for the phase shift error.In these cases the error is almost constant, however with a small increase as the CFL number decreases.The origin of this behaviour has not been further investigated, because the increase is very small.For interFoam the maximum velocity error decreases when lowering the CFL number as also indicated by Larsen et al. (2019), where   = 0.05 is recommended for simulations performed with interFoam.It is interesting to note that our convergence study shows that the wave height error and the mean velocity error for interFoam do not decrease as the CFL number decreases.The reason for the constant or increasing mean velocity error and wave height error could be that the wave travels a shorter distance of 10 compared to 100 in Larsen et al. (2019).In their work the wave height is observed to be fairly constant in the first 10.

Experiences from our implementation of the Ghost Fluid Method
A study of the GFM method combined with isoAdvector according to Vukčević et al. (2018), that we implemented in OpenFOAM-v1912 was the initial starting point for the development of the DSD method.This included a study of the stream function wave propagation case.Fig. 20 presents the velocity field from a simulation performed with our implementation of the GFM method combined with isoAdvector, and Fig. 21 presents the velocity field simulated using our implementation of the GFM method solver combined with a fourth-order explicit Runge-Kutta time integration of the pressurevelocity coupling.The velocity field near the free surface cells at the crest is smeared and the velocity is overpredicted under the wave.This observation was confirmed by Vuko Vukčević (personal communication, October 20, 2021).It was suggested that not taking into account the density jump in the viscous term could result in this behaviour.However, the shown simulations do not include diffusion in the numerical implementation, and the velocity field is still smeared.Several other aspects could also cause the problem: (1) the temporal discretisation, (2) the convective discretisation and (3) the assumption of a continuous velocity field at the free surface.
The temporal term, (1), approximates the acceleration using two previous time steps.When the cells near the free surface change from water to air cells or vice versa, the changes in the velocity field can be very sudden.The convective term, (2), is discretised using the linear upwind scheme, which combines the implicit upwind scheme with an explicit correction based on the cell centre velocity gradients.The velocity gradients will for some cells be affected by sudden changes in velocity field direction over the free surface.Furthermore, the convective term uses the volumetric flux, which is computed by linear interpolation of the velocity field.This could also contribute to the error.The assumption of a continuous velocity field, (3), allows a single continuity equation for the entire domain.This means that one should resolve the boundary layer close to the free surface in order to get a correct solution.That would require very fine mesh and it is often not possible in practical engineering applications.
Our results with the GFM method points in the same direction as the results presented by Vukčević et al. (2018), where the first order wave crest decayed with 3.81% after 8 wave lengths for the finest case with 28 cells per wave height and 140 cells per wave length.The first order horizontal velocity component had an error of 4.2% after 7 wave lengths.The horizontal velocity field presented in Fig.
29 (Vukčević et al., 2018) shows that there is a smooth transition of the horizontal velocity at the wave crest, which is also the case for our implementation.
No direct comparisons of velocity profiles were given by Vukčević et al. (2018).Even though the two implementations, (Vukčević et al., 2018) and ours, give similar results we cannot rule out the possibility of differences in the implementation.

Complex free surface flows
The main focus of the present paper has been to validate the velocity fields for simple cases.However, the DSD solvers are indeed capable of simulating complex free surface flows in 2D as well as 3D.The DSD solvers have been used to simulate various 2D cases with breaking waves (spilling and plunging), a dam break and a drop impact in a liquid pool.In 3D the DSD method has been used to simulate regular waves passing cylindrical and square vertical piles and a breaking phase focused wave impacting a vertical cylinder.Fig. 22 shows a snapshot from a simulation of a breaking wave generated with a single piston stroke.The simulated case is based on the case from Wei et al. (2018).The simulation is performed with DSD-RK2.The presented snapshot also shows that the numerical method can handle occluded air pockets without any problems.A movie of the simulation is available from the supplementary material Video S1.The movie shows that DSD-RK2 is capable of simulating a plunging breaking wave event including the post breaking region with a chaotic and violently changing free surface shape.DSD-RK2 was also used to simulate a 3D case of a phase focused breaking wave impact on a vertical cylinder.Fig. 23 presents a single time frame from the simulation, where the focused wave have just impacted the vertical cylinder.A movie of the simulation is available from Video S2 in the supplementary material.Hence, the DSD solvers are fully capable of handling violent and chaotic flows in 2D as well 3D.These simulations are solely intended to show that the DSD method can handle complex free surface simulations.Future work will focus on validating the accuracy of the DSD solvers on cases involving complex free surface flows.

Conclusion
The present study has developed and implemented a method for free surface flows in the OpenFOAM CFD library.The focus was on improving the description of the wave kinematics in the upper part of the water column for instance for analysing wave interaction with local parts of a structure.The problems using an approach with a varying density combined with a VOF-like propagation of the free surface were illustrated with an inclined square box with calm water.The two methods using a varying density, interFoam and interIsoFoam, clearly generated spurious velocities at the free surface in an else calm water.This phenomenon might be the underlying problem in getting reliable wave kinematics in the upper part of the water column.The present study has demonstrated that the new DSD method predicts a smooth surface elevation without any wiggles as seen from interFoam results.Furthermore, the velocity field is smooth over the entire depth.The DSD solvers are more accurate than the existing solvers interFoam and interIsoFoam from OpenFOAM-v1912, with a lower computational cost.For the finest mesh tested DSD-RK2 is 4.5 times faster than interFoam and several times more accurate.The DSD-RK4 solver is 1.5 times slower than DSD-RK2.It is as accurate as DSD-RK2 in some cases and less accurate in other cases.DSD and DSD-RK2 are very equal to each other in accuracy with minor variations.
We have presented an extensive study of the accuracy of the numerical methods for both constant and varying CFL numbers, which underlines the increased accuracy for the DSD solvers.The study focused on a single test case to enable a more extensive study, which can be widened to include other wave conditions with different steepness and relative depth.
From our own implementation of the GFM method, we have found that accurately predicting the pressure is not sufficient to obtain accurate prediction of the velocity kinematics near the wave crest.In the work by Vukčević et al. (2018) the presented results using 26 cells per wave height gave errors of 3.81% for the surface elevation.This is significantly higher than using our DSD-RK2 solver with 15 cells

Fig. 1 .
Fig. 1.Illustration of projection of a scalar gradient on a cell face and discretisation of the free surface.

Fig. 3 .
Fig. 3. Illustration of extrapolation case 1-3.The blue cells indicate water cells before the advection step giving the free position at the new time instance.The green cells illustrate a single layer of air cells next to the fluid region.After the advection step the new cells entering the fluid region is assumed to be within the green layer of cells.The white cells are air cells.

Fig. 4 .
Fig. 4. Velocity magnitude after 1000 time steps.Case: Inclined box with still water level.Mesh: 50 × 50,  = 0.001 s.The air region of the interFoam and interIsoFoam simulation have been clipped out to highlight the effect in the cells filled by water.

Fig. 5 .
Fig. 5. Surface elevation centred around the last wave crest identified in interval  = 10 − 11 at ∕ = 20.  is the position of the theoretical wave crest.

Fig. 7 .
Fig. 7. Envelope of surface elevation from the stream function wave simulated with DSD-RK4 using the finest resolution of ∕ = 15.

Fig. 8 .Fig. 9 .
Fig. 8. Velocity magnitude vector plot around wave crest given by the stream function wave theory.

Fig. 16 .
Fig. 16.Average error under the last wave crest within  = 10 − 11 over the last 5 s of each simulation.(a) Maximum error of the velocity magnitude profile (b) Mean error of the velocity magnitude profile.

Fig. 17 .
Fig. 17.Average error under the last wave crest within  = 10 − 11 over the last 5 s of each simulation.(a) Error of the wave height (b) Phase shift error  ℎ  ∕.

Fig. 19 .
Fig. 19.Error under the last wave crest within  = 10 − 11 averaged over the last 5 s of the simulation.(a)Maximum error of the velocity magnitude profile (b) Mean error of the velocity magnitude profile (c) Error of the wave height (d) Wave crest phase shift error  ℎ  ∕.

Fig. 15
Fig.15presents the vertical velocity profile below the zero down crossing of the last wave identified in  = 10 − 11 at ∕ = 20.The small boxes show a 400% zoom of the velocity profile just below

Table 2
Case properties.

Table 3
Boundary conditions.