Highfidelity simulation of regular waves based on multimoment finite volume formulation and THINC method
Abstract
The performance of interFoam (a widely used solver within OpenFOAM package) in simulating the propagation of water waves has been reported to be sensitive to the temporal and spatial resolution. To facilitate more accurate simulations, a numerical wave tank is built based on a highorder accurate NavierStokes model, which employs the VPM (volumeaverage/pointvalue multimoment) scheme as the fluid solver and the THINC/QQ method (THINC method with quadratic surface representation and Gaussian quadrature) for the freesurface capturing. Simulations of regular waves in an intermediate water depth are conducted and the results are assessed via comparing with the analytical solutions. The performance of the present model and interFoam solver in simulating the wave propagation is systematically compared in this work. The results clearly demonstrate that compared with interFoam solver, the present model significantly improves the dissipation properties of the propagating wave, where the waveforms as well as the velocity distribution can be substantially maintained while the waves propagating over long distances even with large time steps and coarse grids. It is also shown that the present model requires much less computation time to reach a given error level in comparison with interFoam solver.
keywords:
regular waves, OpenFOAM, highorder accurate model, finite volume method, numerical dissipation1 Introduction
The finite volume method (FVM) has gained great popularity in numerical simulations of water waves for its better conservativeness and flexibility to adapt to both structured and unstructured grids. Conventional FVM solvers typically employ secondorder accurate spatial discretizations, among which the industrystandard solver STARCCM+ and the opensource code interFoam (or other solvers build on interFoam, e.g. waves2Foam jacobsen2012wave (), IHFOAM higuera2013realistic () and interFlow roenby2017new ()) are becoming increasingly popular and have been reported to deal successfully with oceanographic and coastal engineering applications jacobsen2014formation (); jacobsen2014formation02 (); brown2016evaluation (); higuera2013simulating (). However, wave tanks based on secondorder FVM are expected to suffer from numerical dissipation, especially along the direction of wave propagation. Simulations of regular waves using STARCCM+ software peric2015generation () and interFoam solver cha2011numerical (); larsen2018performance () demonstrate that the simulated surface elevations are sensitive to the temporal and spatial resolution, decreasing wave heights and phase shifts are observed in those cases with large time steps and coarse grids. The velocity profiles along the vertical crosssection are also examined by some authors larsen2018performance (); wroniszewski2014benchmarking (), severely overestimated velocities are found near the free surface not only in interFoam simulations but also in the solutions with other NavierStokes solvers, including Gerris and Thétis.
Refining the grids and time steps is a common strategy to minimize those undesirable effects, commonly hundreds of cells are used within a wavelength jacobsen2012wave (); larsen2018performance (); jin2014numerical (); hu2016numerical (); windt2017assessment (), whereas refinement will indisputably lead to an increase in computational cost. The cost could be extremely expensive for some cases, for instance, those involving shortperiod waves, since the relatively higher frequencies cause greater numerical dissipation. In view of this, it is natural to pursue highorder accurate models for the simulations of water waves. Conventional FVM requires wide cell stencil to generate highorder reconstructions, the choice of the stencil, however, is not straightforward on unstructured grids. The numerical schemes with a compact stencil, on the other hand, add degrees of freedoms (DOFs) locally on each cell for highorder reconstructions, and are more flexible and easier to implement on unstructured grids. The representative methods of this sort are the discontinuous Galerkin (DG) method cockburn1989tvb (), the spectral finite volume (SV) method wang2002spectral () and the CIP (constrained interpolation profile) type multimoment finite volume methods (CIP/MM FVM) ii20054th (); xiao2006unified (); chen2008shallow (). The VPM scheme xie2014multi () is an extension of the CIP/MM FVM to incompressible NavierStokes equations on unstructured grids with triangular and tetrahedral elements, it is reported to be a superior tradeoff regarding the numerical accuracy and computational cost. In this work, the VPM scheme is utilized to solve the fluid dynamics and the THINC/QQ scheme xie2017toward () is used for the freesurface capturing, where the latter is claimed to have comparable solution quality to other existing VOF methods with PLIC geometric interface reconstruction. A numerical framework for incompressible interfacial multiphase flows xie2017unstructured () is employed to combine the VPM and THINC/QQ schemes. For wave generation, a mass source function lin1999internal () is applied in this work. By positioning the wavemaker inside the solution domain, wave damping wei1995time () can be straightfowardly applied to all domain boundaries, wave reflections from the tank walls (or structures of interest in future studies) can thus be eliminated in simulations. Comparisons will be conducted between the present model and interFoam solver in simulating the propagation of regular waves within a twodimensional (2D) domain.
The rest of this paper is organized as follows. After the description of governing equations and solution procedures for both the present model and interFoam solver in section 2, the computation domain along with the wave generation and absorption methods is demonstrated in section 3. Section 4 presents the simulation results, followed by a discussion on the numerical schemes in section 5. This paper is ended with conclusion remarks in section 6.
2 Governing equations and solution methods
In this work, the incompressible twofluid flows with moving interface are solved based on the onefluid model prosperetti2009computational (). The NavierStokes equations containing the effects of surface tension and gravity are used for both fluids in the same form,
(1) 
(2) 
(3) 
where is the velocity vector with components and in and directions respectively, is the pressure in excess of the hydrostatic part, is the gravity acceleration, the density and the dynamic viscosity coefficient. is the surface tension force formulated by with being the surface tension coefficient and the interface curvature. A mass source term and a momentum source term have been introduced into the equations for wave generation and absorption respectively, the specific forms of which will be given in section 3. The introduction of the mass source term will be embodied in the computation of Poisson equation and the divergence term in the right hand side of equation 3. A volumeoffluid (VOF) function is used with an indicator function distinguishing two kinds of fluids,
(4) 
In the onefluid model, it is assumed that the intrinsic fluid properties, such as density and viscosity, are updated based on the VOF function as,
(5) 
where and are the densities, and are the dynamic viscosity coefficients of water and air respectively.
2.1 The present model
The present model combines two newly developed numerical schemes, namely, the VPM scheme and the THINC/QQ scheme. The fluid dynamic equations are discretized by the multimoment finite volume method, i.e. VPM scheme, while the VOF transport equation is solved by the THINC/QQ scheme. A numerical framework for interfacial multiphase flows on unstructured grids xie2017unstructured () is employed to combine these two schemes. In this framework, the fractionstep method chorin1968numerical () is used to update the numerical solution from time level () to () and the solution procedure is summarized as follows.

Update the velocity field from step () to by solving the VOF function 3 and the advection term 6 simultaneously in the time integration framework of the thirdorder TVD RungeKutta gottlieb1998total (),
(6) 
Update the velocity field from to by solving the diffusion terms,
(7) 
Update the velocity field from to by adding the effects of surface tension, gravity force and damping force,
(8) 
To make the intermediate velocity field satisfy the mass eqaution 1, it must be corrected by the following projection step. First, the pressure field at step is obtained by solving Poisson equation,
(9) then the velocity field is corrected by projecting the pressure field,
(10)
The numerical schemes will be explained with a 2D triangular element for brevity, more details of the numerical schemes and the 3D formulations can be found in xie2014multi (); xie2017toward (). Shown in Fig. 1(left), the vertices are denoted by located at , the boundary segments are denoted by where subscript denotes the index of the th surface .
VPM scheme for incompressible NavierStokes equations
Different from conventional FVM where only the volume integrated average (VIA) variable is memorized and updated at each time step, the VPM scheme treats both VIA and point value (PV) as the computational variables, such that it realizes a highorder reconstruction with a compact stencil. As shown in Fig.1(right), in the multimoment finite volume method, two kinds of moments of the velocity field are used, namely VIA, and PVs, at cell vertices, which are defined by
(11) 
where denotes the volume of cell . For other physical fields, i.e. pressure field and volume fraction , only VIA is used as the dicrete moment for simplicity. Given both VIA and PV, a quadratic polynomial for velocity field on cell can be constructed in the local cordinate system as,
(12) 
where the unknown coefficients are determined by the corresponding constraint conditions in terms of the multiple moments and are identical for all mesh cells. Both the PVs at the vertices, , and the VIA, , are updated at every time step through the projection procedure shown above. The derivatives at the cell center, and , can be computed from the PVs and VIAs of the neighboring cells. By constructing highorder polynomial, the accuracy is significantly improved in the computation of the fluid dynamics.
THINC/QQ scheme for interface capturing
The discrete form of volumeoffluid transport equation can be obtained from the finite volume formulation 3 over control volume as,
(13) 
where the volumeoffluid of each discrete grid element is defined by
(14) 
and denotes the normal velocity on surface , stands for the reconstruction function in the upwinding cell of surface . To compute the surface flux in equation 13, a hyperbolic tangent function in the local cordinate is used to approximate the indicator function for the target cell at each time step,
(15) 
where parameter determines the steepness of the jump in the interpolation function, which has a constant value of 1.5 in this work. represents the interface surface in the standard element of local element, which is approximated as a curved surface by using fully quadratic polynomial that includes the geometric information,
(16) 
where coefficients are evaluated from the interface normal vector and curvature tensor with leastsquare method.
The only unknown in equation 15 that indicates the location of the interface is determined from VOF values using constrained condition 14,
(17) 
as no general analytical expression is available for integration of multidimensional hyperbolic tangent function 15, a fully multidimensional Gaussian quadrature is applied to approximate the integration,
(18) 
where and denote the cordinates and weights of Gaussian points in element . To balance the quadrature accuracy and computational cost, the number of quadrature points in each dimension is specified as 3 in this work. Once the indicator function is determined, the numerical fluxes on each surface can be calculated by Gaussian quadrature and then used to update the VOF value by equation 13.
2.2 interFoam solver
The solution of the momentum equation in interFoam is performed by constructing a predicted velocity field and then correcting it using the ”Pressure Implicit with Splitting of Operators”(PISO) procedures to time advance the pressure and velocity fields. In particular, to evaluate the fluxes on cell boundaries in the computation of advection, a central differencing flux and a nonoscillatory upwind flux are adopted and a TVD conforming flux limiter of is applied to switch between these two flux schemes when using a socalled limitedLinearV scheme. The variable is an indicator function of consecutive gradients and is specified as 1 in this work. As central differencing schemes are secondorder accurate and upwind schemes are firstorder accurate, this solver has its upper limit in solving the advection term as secondorder accurate. For interface capturing, the ”Multidimensional universal limiter with explicit solution” (MULES) limiter is applied to limit the phase fluxes in solving VOF transport equation 3. And an artificial compression term is added to obtain numerical interface compression, such that it attains the form
(19) 
The rest of the numerical settings will be those found in the popular dambreak tutorial and interested readers are referred to jasak1996error () deshpande2012evaluating () for more details.
3 Numerical wave tank
For this work, a numerical wave tank that has a constant water depth is built with a rectangular source region located at m and an elevation of about from the still water level (SWL), a sponge layer that has a length of is set at each end of the tank. The computational domain along with its coordinate system is presented in Fig. 2. The dimensions of the domain are given in Table 1. It is noted that the present wave tank has a relatively large span in horizontal direction (50m), which is designed to facilitate simulations of longdistance wave propagation.
50.0m  0.5m  0.35m  0.12m  0.0667m  0.04m  3.0m 
As mentioned in section 2, two source functions are embedded in the continuity equation 1 and the momentum equation 2 for wave generation and absorption respectively. The mass source function will be determined by the surface displacement of the target wave lin1999internal (), for a linear monochromatic wave with a wave height and a wave frequency , , it will be given as,
(19) 
where stands for the phase velocity and is the area of the source region.
To prevent undesirable wave reflections from the domain boundaries, absorbing regions are employed by adding damping force into the momentum equation, which has a form as
(20) 
where is the absorbing coefficient, a commonly used form wei1995time () is applied in this work,
(21) 
where and are the starting position and length of the absorbing regions, respectively. and are the empirical damping coefficients to be determined via the numerical tests, which in this work are specified as and . It is noted that, the steepness parameter from the THINC function is reduced to 0.5 within the absorbing regions to increase the viscosity and consequently to enhance the absorption efficiency.
4 Numerical results
Linear monochromatic waves in a constant water depth are simulated with both the present model and interFoam solver. The quality of simulated waves will be assessed in terms of surface elevations and velocity distribution. The generated waves have a moderate nonlinearity due to the intermediate water depth condition (), in this paper, results derived from the secondorder Stokes wave theory are referred to as the reference solutions.
4.1 The base case
First, a base case is presented, in what follows unless stated the numerical settings will remain in accord with this base case. The parameters of the target wave used to define the mass source function are collected in Table 2 and a set of grid system discretizes the entire wave tank with roughly 43 cells per wavelength and 10 cells per wave height () in the vicinity of the SWL. The grid remains uniform in the horizontal direction but nonuniform in the vertical direction, which is coarsened gradually from the SWL to the tank bottom. The time step for the base case is set as , which corresponds to 250 time steps per wave period.
1.0s  0.03m  1.425m  1.425m/s  6.2832rad/s  4.4092/m  0.021 
Simulations are performed for 60 ( is the nominal wave period). Figs.4 and 4 show the time series of normalized surface elevations at different positions, i.e. , where refers to the distance between the center of the source region and the wave gauge, is the nominal wavelength. Profiles recorded at (Figs.4a and 4a) fit well the reference solution, indicating that both wave makers can produce target waves constantly without significant reflections from the boundaries. However, results from the distant gauges show that the two solvers have utterly different performance in propagating the wave trains. The results with the present model compare well with the reference solutions (Figs.4bf), the wave shapes are well maintained and no significant phase errors can be observed even after traveling a long distance of . While with interFoam solver, generated waves have a noticeable amplitude decay after propagating a distance of , in the meantime the simulated profiles start to lead the reference solutions (Figs.4bf). Furthermore, a snapshot of the wave tank at time in Fig.5 intuitively shows that with interFoam solver, the wave decay aggravates as the wave trains propagate along the tank and there is also a change in the wavelength, resulting in phase lags over location. The wave trains simulated with the present model, on the other hand, are well established inside the wave tank, the wave height and wavelength are well maintained and agree well with the reference solutions throughout the working area. Further validation of the wave generation and wave absorption will not be implemented in this work, main attention will be paid to the numerical dissipation in wave propagation.
4.2 Time step study
After a preliminary validation of the wave tank, a time step study and a grid study are carried out with the aim of investigating the sensitivity of the solvers to the temporal and spatial resolution. It is noteworthy that in this section, both the time and spatial series of the surface elevations will be plotted to demonstrate the phase consistency with regard to temporal and spatial distribution. Despite the working area of the tank has a length over 20, the time series recorded at will be assessed in the view of engineering applications.
For the time step study, simulations are performed with different time steps (0.008s, 0.004s and 0.002s) and with otherwise the identical setup to the base case. An enlarged view of the surface elevations recorded at is presented in Fig.7 and a snapshot of the wave tank at is shown in Fig.7. It is demonstrated that the simulations with the present model are basically insensitive to the temporal resolution, the results are substantially consistent with the reference solutions even using a large time step of 0.008s. The interFoam simulations, on the other hand, show significant wave decay and phase shifts in the 0.008s case and the wave profiles approach the reference solutions as the time step reduces, but noticeable deviation exits even with an extremely small time step of 0.0002s that roughly corresponds to a maximum Courant number of 0.01.
The phase errors in the propagation of regular waves simulated with interFoam solver have also been reported in larsen2018performance () and paulsen2014forcing (), they are supposed to arise from the changes in wave periods and wavelengths, so we collect the normalized wave period and the normalized wavelength in Table 3. The simulated wave period is computed by averaging that of the last 30 waves recorded at , while the predicted is calculated from the surface elevations along the wave tank at . It can be seen from the table that the present model yields a constant wave period and wavelength that are highly consistent with the nominal values. While reduced wave periods and expanded wavelengths are found in the interFoam simulations, the disparities increase with the time step and have the maximum values in the 0.008s case, including a decrease of 0.4 in wave periods and a wavelength increase of 2.5. These changes in wave periods and wavelengths account for the phase shifts found in the surface elevations, though the changes are not significant, noticeable phase shifts are produced due to the cumulative effects. The normalized wave height is collected in the table as well, where the strategy for evaluating is utilized again to compute the simulated wave height . It can be seen that the wave decay rates reach the maximum values when a large time step of 0.008s is applied, which are 2.6 and 26.8 for the present model and interFoam solver, respectively.
Present model  interFoam solver  

Resolution  0.008s  0.004s  0.002s  0.008s  0.004s  0.002s  0.0002s 
1.000  1.000  1.000  0.996  0.996  0.997  1.000  
0.999  0.999  0.999  1.025  1.013  1.005  1.000  
0.974  0.997  0.997  0.732  0.801  0.849  0.938 
Also of great interest is the velocity distribution beneath the free surface, as the water particle kinematics provide the basis for force calculations on coastal and offshore structures, while also relating to e.g. shear stresses and hence sediment transport predictions. In Fig.8, the velocity profiles along the vertical crosssection are plotted for three characteristic positions, namely, the trough, the node and the crest of the progressive waves. Reference solutions derived from the potential function for a secondorder Stokes wave are shown together for comparison. The velocity profiles for the present model using a time step of 0.004s compare well with the reference solution, only slightly underestimated velocities can be found beneath the wave crest in Fig.8c. While the interFoam simulations display different features. In the 0.004s case, the vertical velocity beneath the wave node is underestimated (Fig.8b), meanwhile the horizontal velocity beneath the wave trough and the wave crest is underestimated in the bulk of the water, especially in the vicinity of the free surface, is severely overestimated (Figs.8a and c). By reducing the time step to 0.0002s, interFoam solver could have comparative solutions with the present model. The presence of spurious currents near the free surface has been reported by other authors in the cases of both surface tension driven flows and gravity driven flows roenby2017new (); larsen2018performance (); francois2006balanced (); wemmenhove2015numerical (). Especially, in a solitary wave simulation by Wroniszewski et al. wroniszewski2014benchmarking (), it is shown that the spurious currents around the free surface were developed not only with OpenFOAM code but also in the cases of Gerris and Thétis. The reason for such behavior is believed to arise from an imbalance in the discretized momentum equation near the free surface and may explain the decreasing wave periods and increasing wavelengths in the interFoam simulation. The exact mechanism, however, is not yet clear, and will be pursued in further study.
4.3 Grid study
For the grid study, simulations are performed on two other sets of grid, namely, a coarser grid with 21 cells per wavelength and 5 cells per wave height () and a refined grid with 86 cells per wavelength and 20 cells per wave height (). The time steps are modified at the same time in order to maintain a constant Courant number. As shown in Figs.10 and 10, with the present model, surface elevations are substantially consistent with the reference solutions when relatively fine grids are applied, while a moderate loss in wave height, as well as phase shifts, can be observed with the coarsest grid. The interFoam simulations evidently show that using higher spatial resolution improves the results, however, significant deviations from the reference solutions still can be found even with the finest grid. It is remarkable that the profile phases show opposite features with two numerical solvers. Compared with the reference solutions, the profiles simulated with the present model lag behind in the time series and have a lead over location, while the interFoam simulations are on the contrary. Following the method utilized in the time step study, we collect the wave parameters in Table 4. It is demonstrated that the case of the present model yields a correct wave period but a shorter wavelength compared with the reference solution, where the latter is obviously the cause of the phase errors. Different to the present model, interFoam solver produces shorter wave periods and longer wavelengths, resulting in the phase lag in the time series and phase lead over location. It is noted that in the case of the present model, the normalized wave height is a bit over 1.000, the reason for this may include two aspects, the high spatial resolution brings scarcely any numerical dissipation, in the meantime the residual reflection may cause an increase in wave height. Corresponding velocity profiles are shown in Fig.11, similar findings with those in time step study can be noticed from the comparison and will not be elaborated in this section.
Present model  interFoam solver  

Resolution  
1.000  1.000  1.000  0.994  0.996  0.998  
0.992  0.999  0.999  1.029  1.013  1.009  
0.882  0.997  1.001  0.638  0.801  0.903 
4.4 Propagation of waves with different steepness and periods
After the time step study and grid study, the numerical dissipation in simulating regular waves with various steepness and periods is investigated in this section. The parameters used to define the mass source for different cases are shown in Tables 5 and 6, the grid system and the time step of the base case are employed in this section as well, where an identical wave period is applied in the steepness study and vice versa. The maximum wave steepness applied in this section is 0.028, which is much smaller than the breaking point of gravity waves, 0.142, this guarantees that no physical dissipation will be introduced into the simulations. Numerical tests are performed with both the present model and interFoam solver. Numerical dissipation is assessed by looking at the wave heights along the wave tank. The same method mentioned in the time step study is utilized to measure the wave heights, in particular, the last 20 periods instead of 30 periods of the time series will be used to compute the value at as it has a relatively late start.
Wave heights normalized by at various distances to the source region are plotted in Figs.13 and 13. From Fig.13, it can be seen that the two numerical solvers have totally different features in propagating waves with various steepness. With the present model, the wave decay remains at a low level, a moderate wave decay is observed in the case, the decay rate decreases as the steepness increases and practically no dissipation can be found in the cases. This phenomenon may seem surprising, as one would rather expect an expanding wave decay with increasing wave steepness, which is exactly the feature shown in the interFoam simulations. The wave decay in small steepness cases is believed to lie in the relatively low spatial resolution in direction since smaller wave heights are adopted (roughly 3 cells per wave height for the case), while for large steepness, the wave decay is mainly attributed to the nonlinearity effects and the large velocity gradients in the vicinity of the free surface. Apparently the present model is efficient in dealing with the large steepness cases but has a limitation in terms of spatial accuracy. Nevertheless, the present model still has a better performance than interFoam solver with regard to small steepness as it only has a maximum wave decay of in the case comparing to that of with interFoam solver. Fig.13 shows the simulations with various wave periods, the most notable feature is that when simulating a short wave of s, both solvers experience a dramatic wave decay, which is for the present model and for interFoam solver at a distance of to the source region. The low spatial resolution in both and directions is supposed to be the primary cause of the massive wave decay, indicating that simulations with small wave periods need to be performed with care even using highorder accurate models, grid refinement is necessary for this sort of simulations.
(s)  (m)  (m)  

1.0  1.425  0.01  0.007 
1.0  1.425  0.02  0.014 
1.0  1.425  0.03  0.021 
1.0  1.425  0.04  0.028 
(s)  (m)  (m)  

0.7  0.760  0.011  0.014 
1.0  1.425  0.020  0.014 
1.5  2.488  0.035  0.014 
4.5 Computational cost
The present model has shown particular promise in accuracy benefits, many would be concerned with the computational efficiency since highorder schemes are usually more complex and computationally expensive compared to the conventional FVM. So in this section, we will quantify the computational cost of both solvers in terms of computation time on equivalent hardware, as well as quantitative accuracy measurements based on suitable error metrics. All the cases from the grid study in section 4.3 are investigated, which run on a processor of Intel Core i74790 CPU (3.60GHz) and in each case, only a single core is utilized. To evaluate the accuracy, we consider the wave height error instead of the widely used error, where the latter is sensitive to the phase error of the surface elevations that is less concerned in engineering applications. The wave height error is defined as:
(21) 
where the simulated wave height can be found in Table 4. A plot of the wave height error against computation time is shown in Fig.14, where the computation time refers to the execution time spent by the processor and is measured in seconds. It can be seen that with both solvers, the computation time increases while the wave height error is reduced as the resolution increases. The present model outperforms interFoam in terms of accuracy on all three sets of grid, while using the present model is more costly than using interFoam for a specific grid system. For the case, using the present model reduces the error to about compared with interFoam solver, meanwhile it costs more computation time, such a performanceprice ratio is remarkable in view of the fact that interFoam solver takes 96 times longer to achieve a comparative accuracy (see the case of interFoam). In other words, to reach a given error level, for instance, , the present model merely requires a small percentage of computation time that is used by interFoam. Further, as finer grids are applied, the present model demonstrates a faster convergence rate as well as a slower growth in computation time, resulting in an even higher performanceprice ratio, with the case 2 orders of magnitude more accurate than that of interFoam solver at a cost of 37 more computation time.
5 Discussion
Although the exact mechanisms of the waveform changes and the velocity oscillations are not clear at present, some numerical tests are carried out to help better understand the issues behind the behavior and narrow down the dominant factors in the solution procedure. Two interface capturing schemes (MULES and THINC/QQ) and two finite volume methods (the secondorder accurate FVM applied in interFoam and VPM) compose 4 different solvers and all of them are employed to simulate the base case, corresponding surface elevations and velocity profiles are given in Figs.15 and 16.
Both the surface elevations and the velocity profiles demonstrate that simulation results with the same momentum solver are close to each other and the VPMbased solvers produce more accurate results, indicating that the impact of the interface capturing scheme is comparatively limited in this sort of applications. In addition, though not shown herein for brevity, we note that the diffusion terms and the surface tension term make little contribution to the behavior as scarcely any difference occurs in surface elevations and velocity fields when they are turned off. Then it is evident that the advection term is the dominant force in this application. Further, we replace the thirdorder RungeKutta scheme with the firstorder forward Euler scheme to time advance the advection term, the results vary slightly, suggesting that it is the spatial discretization of the advection term that plays a crucial role in the simulations.
6 Conclusion remarks
A 2D wave tank is built in this work, where the VPM scheme and THINC/QQ method are employed to solve the fluid dynamics and capture the free surface, respectively. A mass source function for wave generation and sponge layers for wave absorption are embedded in the present highorder accurate model as well as the opensource code interFoam. The performance of both solvers in simulating the propagation of regular waves is carefully compared in this work.
First, a regular wave in an intermediate depth is simulated as the base case, the surface elevations are compared with the reference solutions given by the secondorder Stokes wave theory. It is demonstrated that both wave makers can produce wave trains constantly without significant reflections from the boundary walls. Wave trains are well reproduced inside the tank with the present model and compare well with the reference solutions even after propagating a distance of 20 wavelengths. With interFoam solver, however, simulated waves show a nonnegligible wave decay and phase shifts after traveling 6 wavelengths, which is believed to be attributed to the numerical dissipation. Then a time step study and a grid study are conducted, special attention has been paid to the numerical dissipation along the propagation direction of the wave trains. By examining the surface elevations, numerical dissipation is evaluated in terms of wave heights and phase consistency. The present model shows a great advantage in simulation fidelity as the simulated results are highly consistent with the reference solutions and shown to be basically insensitive to temporal and spatial resolution, only a moderate wave decay and a slight wavelength decrease are found when a relatively coarse grid is applied. The interFoam simulations, on the other hand, show significant wave decay and phase shifts, where the latter are proved to arise from the reduced wave periods and the expanded wavelengths. The results with interFoam solver can be improved by increasing the temporal and spatial resolution, but deviations from the reference solution still exit even using an extremely small time step of 0.0002s that corresponds to a Courant number of 0.01.
The velocity profiles along the vertical crosssection have also been examined. The present model produces a satisfactory velocity field even with a relatively large time step and a coarse grid. However, the velocities in interFoam simulations are found to be underestimated in the bulk of the water and severely overestimated near the free surface. By reducing the time step to 0.0002s, interFoam solver can have comparative velocity distribution with the present model using a time step of 0.004s. Propagation of regular waves with various steepness and periods is investigated as well. The results indicate that the present model has a much better performance than interFoam solver. While in the short wave case, i.e. s, both solvers experience a dramatic wave decay, which is attributed to the relatively low spatial resolution. This is a reminder that simulations involving shortperiod waves need to be conducted with care, grid refinement is necessary even applying a highorder accurate model. A comparison of the computational cost shows that using the present model causes increase in simulation time for a given grid system, whereas the present model outperforms interFoam solver in terms of accuracy by up to 2 orders of magnitude. Consequently, the present model requires much less computational cost to reach a given error level in comparison with interFoam solver.
Further model analysis suggests that the spatial discretization of the advection term, instead of the interface capturing scheme as many would expect, plays a crucial role in this application. Although the interface capturing scheme THINC/QQ has not demonstrated its advantages over MULES scheme in simulating the propagation of regular waves, it is believed to be promising for applications involving violentlychanging interfaces, for instance, wave breaking and wavestructure interactions.
Acknowledgements
The early stages of this study were conducted during the first author’s stay at Tokyo Institute of Technology. The authors are grateful to the helpful discussions with some members of Xiao Lab. Special acknowledgements should be given to Dr. Feng Xiao. This study was partially supported by the National Natural Science Foundation of China (Grant Nos. 51479175, 51679212), Zhejiang Provincial Natural Science Foundation of China (Grant No. R16E090002).
References
 N. G. Jacobsen, D. R. Fuhrman, and J. Fredsøe, “A wave generation toolbox for the opensource cfd library: Openfoam®,” International Journal for Numerical Methods in Fluids, vol. 70, no. 9, pp. 1073–1088, 2012.
 P. Higuera, J. L. Lara, and I. J. Losada, “Realistic wave generation and active wave absorption for navier–stokes models: Application to openfoam®,” Coastal Engineering, vol. 71, pp. 102–118, 2013.
 J. Roenby, B. E. Larsen, H. Bredmose, and H. Jasak, “A new volumeoffluid method in openfoam,” in VII International Conference on Computational Methods in Marine Engineering, MARINE, vol. 2017, 2017, pp. 1–12.
 N. G. Jacobsen, J. Fredsoe, and J. H. Jensen, “Formation and development of a breaker bar under regular waves. part 1: Model description and hydrodynamics,” Coastal Engineering, vol. 88, pp. 182–193, 2014.
 N. G. Jacobsen and J. Fredsoe, “Formation and development of a breaker bar under regular waves. part 2: Sediment transport and morphology,” Coastal Engineering, vol. 88, pp. 55–68, 2014.
 S. Brown, D. Greaves, V. Magar, and D. Conley, “Evaluation of turbulence closure models under spilling and plunging breakers in the surf zone,” Coastal Engineering, vol. 114, pp. 177–193, 2016.
 P. Higuera, J. L. Lara, and I. J. Losada, “Simulating coastal engineering processes with openfoam®,” Coastal Engineering, vol. 71, pp. 119–134, 2013.
 R. Perić and M. AbdelMaksoud, “Generation of freesurface waves by localized source terms in the continuity equation,” Ocean Engineering, vol. 109, pp. 567–579, 2015.
 J.J. Cha and D.C. Wan, “Numerical wave generation and absorption based on openfoam,” Ocean Engineering(Haiyang Gongcheng), vol. 29, no. 3, pp. 1–12, 2011.
 B. E. Larsen, D. R. Fuhrman, and J. Roenby, “Performance of interfoam on the simulation of progressive waves,” arXiv preprint arXiv:1804.01158, 2018.
 P. A. Wroniszewski, J. C. Verschaeve, and G. K. Pedersen, “Benchmarking of navier–stokes codes for free surface simulations by means of a solitary wave,” Coastal Engineering, vol. 91, pp. 1–17, 2014.
 H. Jin, Y. Liu, S.y. He, H.j. Li et al., “Numerical study on the wave dissipating performance of a submerged horizontal plate breakwater using openfoam,” in The Eleventh ISOPE Pacific/Asia Offshore Mechanics Symposium. International Society of Offshore and Polar Engineers, 2014.
 Z. Z. Hu, D. Greaves, and A. Raby, “Numerical wave tank study of extreme waves and wavestructure interaction using openfoam®,” Ocean Engineering, vol. 126, pp. 329–342, 2016.
 C. Windt, J. Davidson, P. Schmitt, and J. V. Ringwood, “Assessment of numerical wave makers,” in Proceedings of the 12th European Wave and Tidal Energy Conference, 2017.
 B. Cockburn and C.W. Shu, “Tvb rungekutta local projection discontinuous galerkin finite element method for conservation laws. ii. general framework,” Mathematics of computation, vol. 52, no. 186, pp. 411–435, 1989.
 Z. J. Wang, “Spectral (finite) volume method for conservation laws on unstructured grids. basic formulation: Basic formulation,” Journal of computational physics, vol. 178, no. 1, pp. 210–251, 2002.
 S. Ii, M. Shimuta, and F. Xiao, “A 4thorder and singlecellbased advection scheme on unstructured grids using multimoments,” Computer physics communications, vol. 173, no. 12, pp. 17–33, 2005.
 F. Xiao, R. Akoh, and S. Ii, “Unified formulation for compressible and incompressible flows by using multiintegrated moments ii: Multidimensional version for compressible and incompressible flows,” Journal of Computational Physics, vol. 213, no. 1, pp. 31–56, 2006.
 C. Chen and F. Xiao, “Shallow water model on cubedsphere by multimoment finite volume method,” Journal of Computational Physics, vol. 227, no. 10, pp. 5019–5044, 2008.
 B. Xie, S. Ii, A. Ikebata, and F. Xiao, “A multimoment finite volume method for incompressible navier–stokes equations on unstructured grids: volumeaverage/pointvalue formulation,” Journal of Computational Physics, vol. 277, pp. 138–162, 2014.
 B. Xie and F. Xiao, “Toward efficient and accurate interface capturing on arbitrary hybrid unstructured grids: The thinc method with quadratic surface representation and gaussian quadrature,” Journal of Computational Physics, vol. 349, pp. 415–440, 2017.
 B. Xie, P. Jin, and F. Xiao, “An unstructuredgrid numerical model for interfacial multiphase fluids based on multimoment finite volume formulation and thinc method,” International Journal of Multiphase Flow, vol. 89, pp. 375–398, 2017.
 P. Lin and P. L.F. Liu, “Internal wavemaker for navierstokes equations models,” Journal of waterway, port, coastal, and ocean engineering, vol. 125, no. 4, pp. 207–215, 1999.
 G. Wei and J. T. Kirby, “Timedependent numerical code for extended boussinesq equations,” Journal of Waterway, Port, Coastal, and Ocean Engineering, vol. 121, no. 5, pp. 251–261, 1995.
 A. Prosperetti and G. Tryggvason, Computational methods for multiphase flow. Cambridge university press, 2009.
 A. J. Chorin, “Numerical solution of the navierstokes equations,” Mathematics of computation, vol. 22, no. 104, pp. 745–762, 1968.
 S. Gottlieb and C.W. Shu, “Total variation diminishing rungekutta schemes,” Mathematics of computation of the American Mathematical Society, vol. 67, no. 221, pp. 73–85, 1998.
 H. Jasak, “Error analysis and estimation for the finite volume method with applications to fluid flows.” 1996.
 S. S. Deshpande, L. Anumolu, and M. F. Trujillo, “Evaluating the performance of the twophase flow solver interfoam,” Computational science & discovery, vol. 5, no. 1, p. 014016, 2012.
 B. T. Paulsen, H. Bredmose, H. B. Bingham, and N. G. Jacobsen, “Forcing of a bottommounted circular cylinder by steep regular water waves at finite depth,” Journal of fluid mechanics, vol. 755, pp. 1–34, 2014.
 M. M. Francois, S. J. Cummins, E. D. Dendy, D. B. Kothe, J. M. Sicilian, and M. W. Williams, “A balancedforce algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework,” Journal of Computational Physics, vol. 213, no. 1, pp. 141–173, 2006.
 R. Wemmenhove, R. Luppes, A. E. Veldman, and T. Bunnik, “Numerical simulation of hydrodynamic wave loading by a compressible twophase flow method,” Computers & Fluids, vol. 114, pp. 218–231, 2015.