Generation of Water Waves Using Momentum Source Wave-Maker Applied to a RANS Solver

Nowadays, as the development of Computational Fluid Dynamics (CFD) and the numerical wave tank (NWT) has advanced, numerical analysis has become increasingly useful and powerful for the ship designing and ship hydrodynamics. In this study, a momentum source wave-maker and an analytical relaxation wave absorber were embedded into 2D RANS equation model with RSM turbulence closure scheme to establish the NWT for ship designing and hydrodynamics. The VOF (volume-of-fluid) method was applied to accurately capture the water free surface.The body force-weighted scheme is chosen for pressure interpolation and the second order upwind scheme for discretization of themomentumequation. In order to calculate convection and diffusion fluxes through the control volume faces, PISO algorithm is adopted for pressure-velocity coupling. The momentum source function for wave generation and the analytical relaxation function for wave absorption were deduced for constructing the NWT (numerical wave tank). The proposed NWT was then validated by the laboratory measurements of Umeyama and the analytical solution, indicating that the constructed NWT is effective and accurate.


Introduction
Recently, the computer technology and numerical methods used in marine science and technology have evolved a lot including these aspects: ocean wave power generation, marine resources exploitation, and marine transportation.All these are related to the numerical simulation of ocean waves.The ocean waves can be categorized into linear waves, high order stokes waves, solitary waves, irregular waves, and so on.The physical methods are very limited due to the aspects of high cost (for establishing the experimental wave tank and for flow measuring), inconvenience for maintenance and management, and hardness for monitoring of some specific flow data leading to the ignorance of some important details [1][2][3][4].As the development of computer technology and numerical methods based on Computational Fluid Dynamics (CFD) has advanced, numerical simulation has become increasingly useful and powerful for showing the complicated flow fields encountered in developing the marine science and technology under water wave and current conditions.First of all, it is essential to establish a numerical wave tank accurately.Numerical wave tanks are computer codes whose main goal is to reproduce physical wave basins as closely as possible, which have solved many complex water wave problems as an effective tool.As an important influencing aspect of the marine science and technology, the ocean waves have been analyzed a lot recently based on the numerical wave tanks (NWT).At present, there are many models used to establish the numerical wave tanks, and the Navier-Stokes or Reynolds Averaged Navier-Stokes (RANS) equation models have been widely used.The Navier-Stokes equations can accurately describe water wave problems, but the direct numerical simulation (DNS) of the Navier-Stokes equations which solves the entire spectrum of motions (ranging from large eddy motions to the smallest turbulence scale motions) as well as nonlinear terms, viscous and turbulent stresses, has low computational efficiency and is always restricted to the hardware limitations.Reynolds-Averaged Navier-Stokes (RANS) equations [5,6], which were firstly applied in marine technology about 20 years ago, are alternatives that hugely reduce computational load compared to that of DNS.The studies on establishing the numerical wave tank using the RANS equations have been increasing, among which three types of wave-makers have been developed: static boundary wave-maker, moving boundary wavemaker, and internal wave-maker [7].
With the application of static boundary wave-maker, waves are generated at the inlet boundary given specified velocities and the water surface elevation derived from any wave theory formulation [8].As pointed out by Pengzhi Lin [9], the static boundary wave-maker will produce obvious numerical dissipation due to the strong wave reflection.With the application of moving boundary wave-maker, the motion of a realistic wave-maker (flap-type wave-maker or pistontype wave-maker) is numerically simulated by a moving boundary.For example, Wood et al. [10] analyzed the runup of steep nonbreaking waves using the piston-type wavemaker based on the FLUENT software, and Finnegan and Goggins [11] simulated the linear water waves and wavestructure interactions using the flap-type wave-maker based on the Ansys CFX commercial software.Since the computational domain changes as the solid body moves toward or away from the fluid, a remeshing is inevitable in each time step or after a large distortion of the generated grid, which should be avoided as suggested by most researchers.The moving boundary wave-maker has also been applied with the meshless method (smoothed particle hydrodynamics) for wave generation.For example, Jia-wen Sun [12] simulated the wave impact on a horizontal deck by establishing a numerical wave tank with the moving boundary wave-maker method based on an improved incompressible smoothed particle hydrodynamics.The application of mesh-less method in establishing the numerical wave tanks is still in an early stage, and there are many technical obstacles, such as the lack of particles nearby the wave-maker, that need to be conquered.
With the application of internal wave-maker, waves are generated by applying a source function or a source line within a designated region inside the computational zone not interacting with reflected waves and being able to be combined with various wave absorber methods.For the source line method, waves are generated at a single point in the wave propagation direction, by adding, at each time step, an incremental water surface elevation computed by the resolution of the model equations.Larsen and Dancy [13] were the first to use the method based on the Boussinesq equations and suggested that the phase velocity is related to the incremental water surface elevation.For the application of the source line method based on the mild-slope equations model, Lee and Suh [14] successfully generated directional monochromatic and random waves.This method was developed and applied to the study of wave energy converters by Beels et al. [15].The source function method may be subdivided into analytical relaxation method, mass source function method, and momentum source function method.Mayer et al. [16] proposed the analytical relaxation method for generating and absorbing waves simultaneously in the relaxation domain.The analytical relaxation method is usually linked to the numerical sponge layer method, and thus the relaxation zone should be enlarged.The mass source function method and the momentum source function method were first proposed by Ge Wei [17] based on Boussinesq-type equations and were developed based on Navier-Stokes equations or RANS equations.The mass source function method based on the Navier-Stokes equations model was first proposed by Lin and Liu by using the finite difference method (FDM) on a staggered grid system.The VOF scheme was adopted to track the free surface.The proposed numerical wave tank was applied to generate various types of wave, and the rules of thumb for designing the source region have become a general principle.J.-S.Zhang et al. [18] applied the mass source function method and the rules of thumb for designing the source region in analyzing the wave-current interactions by using the finite volume method (FVM).Nevertheless, a lot of numerical tests are needed for the mass source function designing for generating a particular wave.The momentum source function method proposed firstly by Junwoo Choi [19] is more convenient than the mass source function method because the width of the momentum source domain is the only parameter to be specified, whereas the mass source function method has three parameters (i.e., height, width, and location of the mass source domain), but the application based on Navier-Stokes equations or RANS equations for simulating the regular and irregular wave propagation over constant water depth has been reported in few studies.
In this study, we deduced the momentum source functions of the RANS equations for the internal wave-maker and the wave absorber based on analytical relaxation method mentioned above.FLUENT software was chosen as the base solver due to its popularity and applicability to water wave problems.All the numerical simulations were run in parallel using series of Intel XEON E5-2640V3 processors (2.60 GHz).To evaluate the application of this method to generate a target wave train in a vertical 2-dimensional channel of constant depth, the numerical results were compared to the laboratory measurements by Umeyama [20] and the analytical solutions.

Governing Equations.
The incompressible fluid motion due to the water wave propagation can be described by the mass conservation equation, momentum conservation equations, and RANS equations: in which  is time,  is the density,   denotes the Cartesian coordinates,   is an ensemble mean velocity component,  is the fluid pressure,  is the dynamic viscosity, and  is the gravitational acceleration.  is the additional mass source term, which is equal to 0 here for the momentum wave-maker method.  is the additional momentum source term which will be deduced in Sections 2.2 and 2.3 for wave generation and absorption, respectively.−      is the Reynolds stress term.
The VOF method is used to track the free water surface.The main idea of VOF method is to define a function  to represent the fractional volume of water fluid.
= 1 indicates that the cell is full of water, while  = 0 corresponds to a cell fully occupied by air.Cells with value of 0 <  < 1 contain a water free surface.

Wave Generation.
As shown above, the water free surface displacement is not a variable of the RANS equations; it is impossible to get the exact relation expression between the water free surface displacement and its momentum source function of the RANS equations for generating the target water waves.However, the water free surface displacement is an unknown variable of the depth-integrated equations and Wei et al. (1999) derived the mass source function based on the Boussinesq approximation for generating water waves.Their idea can be embedded into the RANS equation, and the momentum source function can be calculated via the appropriate numerical technique.In other words, the mass source function of the continuity equation, in which the water free surface is a key variable, will be transformed into the momentum source function of the RANS equations for generating the target waves in this study.
Wei et al. used the following linear Boussinesq equations to derive the mass source function to generate target water waves over constant water depth: in which ℎ is the water depth,  the gravitational acceleration,  the surface elevation, and  ⇀  the velocity vector.The parameters  1 and  used her are  1 =  + 1/3 and The velocity vector can be expressed as  ⇀  = ∇, in which  is the velocity potential.The model equations (( 4) and ( 5)) become as follows.
For internal wave generation the mass source function (, , ) in (4) has been proposed as , where angular frequency , wave number , and wave height  0 are the wave parameters used to obtain a target wave train.And  1 = √/ exp(− 2 /4), where  = 80/() 2 ,  is the wavelength, and  is a parameter for the width of internal generation region.
In order to transform the mass source function derived above into a momentum source function for internal wave generation, the linear Boussinesq equations ( 4) and ( 5) will be transformed into the following equations by deleting mass source function in the continuity equation and the momentum source function will be added.
By demanding (11) on the  partial derivative, we get in which   can be expressed as follows according to (10).
By substituting the   expression (13) into (12), we get the following.
By integrating the above equation, the momentum source function can be expressed as follows.
Using the relation above and the source mass function derived in Wei et al., the x-directional component of the momentum source function can be derived as follows.
Similarly, we can get the y-directional component of the momentum source function.
According to the Fourier transform, the momentum source function for a linear monochromatic wave can be simplified as follows.
For 2D simulation of the propagation of the linear waves over constant water depth, the wave direction angle  = 0 and   = 0.

Wave Absorption.
The analytical relaxation method proposed by Mayer and Madson was used here for wave absorption, and the mechanism of the analytical relaxation method can be described as follows: Within the relaxation domain, the velocity and the pressure will be updated at every time step by the added source.As application to the Navier-Stokes equations, the relaxation algorithm for velocity and pressure can be renewed as (20) in which the subscript  stands for the renewed value of the physical quantity in the specified zone, the subscript  stands for the computed value at previous time step, and  = () is the relaxation function,  min = 1,  max = 0.
Then, the source function for the analytical relaxation wave absorbing method can be deduced from the Euler equations by ignoring the water viscosity.The difference forms of momentum equations with and without the additive source can be written as follows.By subtracting ( 21) with ( 23) and ( 22) with ( 24), the source function can be expressed as follows.

𝜌 𝑢
More details about the source function can be found in [16].

Setup of the Numerical Model. The Computational Fluid
Dynamics code, FLUENT 13.0 [21], which proved to be viable for solving compressible and incompressible flows based on the two-dimensional or three-dimensional N-S equations or RANS equations, was chosen as the base solver.In the following, only the numerical methods used here are briefly described.For turbulent flows, the FLUENT solver supports various turbulence models.When deciding the exact turbulence model for wave generation using the momentum source function method, we found that the Reynolds Stress Model (RSM) converged faster and gives better results than the other turbulence models.Therefore, the RSM using a linear pressure-strain model was used in conjunction with the standard wall function method for a smooth wall for all the present simulations.With the above models, the FLUENT solver, employing the finite volume method (FVM) for the discretion of the governing equations on the basis of staggered grids, generally solved the RANS and employed the volume-of-fluid method (VOF) employing a geometric reconstruction scheme to track free surface movement.Pressure, turbulent kinetic energy, turbulent dissipation rate, and water volume fraction function are arranged at the central point of grid; fluid velocity components are arranged at the central point of corresponding grid boundary.The body force-weighted scheme was chosen for pressure interpolation and the second order upwind scheme for discretization of the momentum equation.In order to calculate convection and diffusion fluxes through the control volume faces, PISO algorithm is adopted for pressure-velocity coupling.The convergence criterion is set to 0.001 for all simulations.For unidirectional wave generation in two-dimensional numerical wave tank, only the momentum source function was embedded into the x-momentum equations without any z-directional variation by using the macro of the UDF-DEFINE SOURCE (mom source, cell, thread, ds, eqn) in software FLUENT.At both ends in the x-direction with thickness of approximately three wavelengths, the analytical relaxation wave absorber, by using the source function expressed by ( 25) and ( 26) embedding into the momentum equations, was used for absorbing the wave energy also achieved by the UDF-DEFINE SOURCE (mom source, cell, thread, ds, eqn).
Figure 1 shows the sketch of the numerical wave tank established by using the momentum source function wave generation method and the analytical relaxation wave absorption method.The no-slip boundary with the smooth wall function was set for left, right, and bottom boundary conditions, and the pressure-outlet boundary was set for the top boundary.The computational zones were discreted by the structured grids using GAMBIT.The grid size and the temporal grid size have huge influence on the computational results, which will be discussed later in detail.

Model Verification.
To validate the proposed model, the numerical results are compared with laboratory measurements of Umeyama.Umeyama performed a series of experiments in a physical wave channel with length of 25m, width of 0.7m, and maximum depth of 1.0m.Waves were generated by placing a piston-type wave-maker at one end of the wave channel and were absorbed by installing a wave absorber at the other end.During all the tests in the channel, the water depth was kept at 0.3m and the wave period was chosen as 1.0s; the parameters for the tests are listed in Table 1.The displacement of the free water surface was measured using a resistance-type probe, located 14.0m from the wave generator, while the horizontal velocities of waterparticles were determined by PIV device and PTV device based on a single-exposure image.For more details about the laboratory measurement, the reader can refer to Umeyama.A computational domain with a range of -20 m to 20 m and 0m to 0.5 m is used in the numerical simulation, and the grid size and the temporal grid size are decided by the following sensitivity analysis.
To check the grid and temporal grid sensitivity of the solutions, the temporal grid sensitivity analysis is carried  out first.The medium grid (93.1MB, containing 1000000 quadrilateral cells) is used and three different time steps are considered, which are 0.01 (T/100), 0.005 (T/200), and 0.002 s (T/500).The wave amplitude is set to 0.00515m, the period is 1 s, and the simulation is carried out for a duration of 25 periods.The momentum source region is located at x=0m to x=0.5m.The momentum source distribution for all the numerical tests is shown in Figure 2. As shown in Figure 2, the amplitude of the source increases with the wave height.
Figure 3(a) shows the time series of surface elevation at x=10m for different time steps, revealing that the results are nearly the same.Consequently, for efficiency, the time step of 0.01s is selected for all cases.Then, three different grid systems of coarse (22.9 MB, containing 250000 quadrilateral cells, Δ = 0.02, Δ = 0.004), medium (93.1MB, containing 1000000 quadrilateral cells, Δ = 0.01, Δ = 0.002), and fine grids (137MB, containing 1440000 quadrilateral cells, Δ = 0.005, nonuniform in  direction with Δ = 0.005 ∼ 0.001) are used to test the grid sensitivity of the solutions.The results are shown in Figure 3(b), where it can be seen that the results are more coincident for the medium grids and fine grids.Considering a further study to capture the fluid interface more exactly and the time consumed, the fine grid system is used for the present simulation.Additional numerical tests show that the numerical results are insensitive to the temporal grid size as long as the temporal grid size is within /100.The grid size Δ in the x-direction was within the range of /50 ∼ /20, the grid size Δ in the z-direction nearby the free surface was within the range of /5 ∼ /10, and the ratio of x-direction grid and z-direction grid was no more than five.
Figure 4 shows the comparison of simulated and measured, analytical water surface profile within one wave period.As expected, there is a good agreement among the simulated results, the physical experimental results, and the theoretical solution.For the numerical results, the error near the wave crest is a little bit larger than the experimental results.As a whole, the numerical results are more close to the theoretical results than to the experimental results.
To demonstrate the effectiveness of the internal wavemaker and the wave absorber, the simulated average wave elevation profiles at time / = 0 in the last five wave periods are compared to the theoretical results as shown in Figure 5.It can be found that the wave elevation profiles in working zone areas for all cases fit well with the theoretical wave elevation profiles, and the analytical relaxation wave absorption method can absorb the wave energy effectively for avoiding the reflected waves generated by the left and right wall boundaries influencing the wave profiles in the working zones.By investigating Figure 4 carefully, a tiny phase discrepancy is observed since the target waves are generated from a narrow region instead of a vertical line.A larger discrepancy is observed in case-W3 (wave steepness  = 0.0827) due to the limitations of momentum source internal wave-maker in generating deep water waves since the momentum source function is derived from linear Boussinesq equations based on the shallow water wave assumption.In general, the momentum source internal wave-maker can generate well-targeted waves by being applied to a RANS solver based on the finite volume method (FVM) with fine spatial and temporal grid.-8 (  =  − 0.3, distance from the free surface) display the comparison of simulated horizontal-and verticalvelocity profiles for the cases-W1, W2, and W3.It can be found that the simulated horizontal-and vertical-velocity profiles agree well with the theoretical results for all cases with tiny discrepancy in cases-W1 and W2 and a larger discrepancy in case-W3 due to the larger wave steepness.The results indicate that the vertical-and horizontal-velocity profiles are mainly affected by the surface wave motion.As shown in the lower parts of Figures 6-8, the absolute value of the vertical velocity above the bottom boundary layer increases toward the water surface (see the velocity profile at time t=0.25s(t=0.75s) of Figures 6-8) when the zero wave elevation arrives.Similarly, the absolute value of horizontal velocity above the bottom boundary layer increases toward the water surface (see the velocity profile at time t=0.0s(t=0.5s) of Figures 6-8) when the wave trough and wave crest arrive, as shown in the upper parts of the Figures 6-8.

Conclusion
In this study, a momentum source internal wave-maker and an analytical wave absorber were embedded into 2D RANS equation model with RSM turbulence closure scheme to simulate the water wave propagation over constant water depth.By comparing the simulated results with the experimental data and the theoretical results, it can be found that the momentum source internal wave-maker can generate welltargeted waves and the analytical relaxation wave absorption method can absorb the wave energy effectively for avoiding  the reflected waves generated by the left and right wall boundaries influencing the wave profiles in the working zones.The established NWT can be effectively used in the ship designing and ship hydrodynamic analysis.

Figure 1 :
Figure 1: An illustrative sketch of the computational domain.

Figure 2 :Figure 3 :
Figure 2: Momentum source distribution map for all cases.

Figure 4 :
Figure 4: Comparison of simulated and measured, analytical water surface profile within one wave period.

Figure 5 :
Figure 5: Comparison between simulated and theoretical free surface elevation.

Figure 6 :
Figure 6: Comparison of simulated and theoretical horizontal-and vertical-velocity profiles for case-W1.

Figure 7 :
Figure 7: Comparison of simulated and theoretical horizontal-and vertical-velocity profiles for case-W2.

Figures 6
Figures6-8(  =  − 0.3, distance from the free surface) display the comparison of simulated horizontal-and verticalvelocity profiles for the cases-W1, W2, and W3.It can be found that the simulated horizontal-and vertical-velocity profiles agree well with the theoretical results for all cases with tiny discrepancy in cases-W1 and W2 and a larger discrepancy in case-W3 due to the larger wave steepness.The results indicate that the vertical-and horizontal-velocity profiles are mainly affected by the surface wave motion.As shown in the lower parts of Figures6-8, the absolute value of the vertical velocity above the bottom boundary layer increases toward the water surface (see the velocity profile at time t=0.25s(t=0.75s) of Figures6-8) when the zero wave elevation arrives.Similarly, the absolute value of horizontal velocity above the bottom boundary layer increases toward

Figure 8 :
Figure 8: Comparison of simulated and theoretical horizontal-and vertical-velocity profiles for case-W3.