A Three-Dimensional Numerical Model with an L-Type Wave-Maker System for Water Wave Simulations by the Moving Boundary Method

: A three-dimensional numerical wave tank was developed based on Reynolds averaged Navier–Stokes equations and the volume of fluid method. The moving boundary method is adopted in this model to generate water waves. Piston-type wave-makers are mimicked for the total replication of the physical wave tank conditions. Two-dimensional regular and irregular waves are simulated, with the capability to trigger the active wave absorption algorithm. The two-sided wave-maker system with L-type arrangement is adopted in this model to expand the effective wave areas for three-dimensional waves. Oblique regular waves and multidirectional random waves are simulated, yielding a good agreement with theoretical solutions. The results indicate that this numerical model is an effective tool to provide finer details or complement data unavailable due to the physical setting of a tank experiment.


Introduction
Within the field of coastal engineering, physical modeling, and numerical simulation are the two main experimental approaches. Physical modeling can exactly reproduce the wave propagation and deformation process in the laboratory, but this method is costly and limited by the setting of the experimental facilities in which they are conducted. With the rapid development of computer technology, numerical modeling has become an effective tool to study water waves and their action on marine structures. Numerical simulation is less expensive and enables finer details difficult to measure via physical experiments-such as three-dimensional vortex in turbulent flow-to be shown.
Over the past few decades, various numerical models have been developed to simulate wave propagation, deformation and the water-structure interaction. The Boussinesq equation models and shallow water equation models are widely used to simulate the wave propagation process in coastal zones. These equations can be seen as simplified Navier-Stokes equations, averaged in the vertical direction. However, the simplified equations are subject to many limitations, among which the impossibility of simulating the wave breaking process without an additional dissipation scheme. The Navier-Stokes equation models present an advantage in simulation accuracy, especially in the case of the direct numerical simulation method, but their huge computational cost remains a big challenge. For this reason, Reynolds-averaged Navier-Stokes (RANS) equation models have been widely used to simulate water waves. Lin and Liu [1] developed COBRAS model, which solves the RANS equations and k-ε equations to simulate the breaking waves in the surf zone. Other RANS equation models have also been developed for water wave simulations, including the works of Troch and Rouck [2], Choi et al. [3], Cheng et al. [4], Del Jesus et al. [5], Jacobsen et al. [6], Higuera et al. [7] and Zhang et al. [8].
The treatment of the wave generation and active absorption boundaries for a RANS numerical model is a critical issue. One commonly used approach to generate waves is the velocity-inlet method, in which the theoretical values of surface elevations and velocities are assigned to the boundary cells. However, this method requires additional corrections to be taken in order to maintain the conservation of total water mass. This was the approach adopted by Del Jesus et al. [5] to generate waves in the IH-3VOF code. Another widely used approach is the internal wave-maker method, which generates waves via a mass source or momentum source. Lin and Liu [9] proposed several mass source functions, which are added in the continuity equation to generate different kinds of waves. Similar functions have been adopted in other numerical models to simulate two-dimensional waves. The internal wave-maker method has also been expanded to three dimensions by Ha et al. [10]. Another method, the relaxation zone method was adopted by Jacobsen et al. [6] to generate and absorb waves in the waves2Foam code. Furthermore, as various types of waves can be generated by piston-type or flap-type wave-makers in physical wave tanks, a straightforward numerical wave generation method consists in replicating the movement of physical wave-maker paddles by the moving boundary method.
The advantage of the moving boundary method is that it enables the thorough replication of a physical wave tank, including the influence of the wave-maker paddle movement. The combination of Navier-Stokes equations and moving boundary method makes the numerical simulation closer to the actual experimental environment. With an effective numerical model, the simulation results can be used to fill the data-void zones of the physical experiment results, especially at the points inaccessible for physical measurements. However, the moving boundary method has not been widely used until now because of the requirement of dynamic meshes. Orszaghova et al. [11] adopted this method in a potential flow model. Higuera et al. [12] applied the dynamic mesh in the OpenFOAM (Open Source Field Operation and Manipulation) code as a wave generation boundary, in which the movement signals were obtained from the physical wave-makers to generate two-dimensional and three-dimensional focused waves. The features of the waves simulated by the Navier-Stokes equations with the moving boundary method, especially for irregular waves and three-dimensional waves or multidirectional waves, have seldom been investigated in previous studies. In addition, in three-dimensional numerical wave tanks, preventing the influence of the gaps between adjacent paddles is an unsolved problem. In this paper, a three-dimensional wave tank is established based on the Navier-Stokes equations and the moving boundary method. The two-sided wave-maker system (L-type multidirectional wave-maker system) is replicated by smooth boundaries to prevent deviations due to the gaps between adjacent paddles. The features of two-dimensional regular waves and irregular waves generated by the moving boundary method in the RANS equation model are discussed, and the method is expanded to three dimensions to simulate oblique waves and multidirectional waves.

Two-Phase Flow Model Description
OpenFOAM is an open source code package for computational fluid dynamics and the release version OpenFOAM-2.3.0 was adopted as our code development platform in this work. The most significant advantage of OpenFOAM is the freedom for users to extend the existing basic solvers by new codes, and the compatibility with some pre-processing and post-processing tools.
The basic solver used in this work is interDyMFoam, which solves the three-dimensional RANS equations for two phases of incompressible, immiscible fluids. The governing equations are continuity and mass conservation equations, written as: where u is the velocity vector, ρ is the density of the fluid, μ is the dynamic viscosity, pd is the pseudo-dynamic pressure, g is the gravitational acceleration vector, x is the position vector and fc is an additional source term for wave damping, whose detailed expression is given in Section 2.4. The volume of fluid (VOF) method is adopted to deal with the interface between two fluids. The equation for the volume fraction function is: in which the third term is an artificial compression term set to reduce the smearing of the interface. In Equation (3), α is the volume fraction of water. cα is a coefficient to control the strength of compression, taking cα = 1 in this work. It means that the value of α is 1 for water and 0 for air and in the case where the value of α is between 0 and 1, the cell contains the interface. It is convenient to calculate the physical properties of the fluid at the interface by volume fraction, as follows: (1 ) water air

Solving Procedure
An algorithm PIMPLE, the combination of the SIMPLE [13] and PISO [14] algorithms, is adopted in the OpenFOAM code. The detailed algorithms and explanations for PIMPLE can be found in [15]. In the simulation cases of this study, the number of outer circulation were set to one, i.e., the algorithm can be seen as the PISO algorithm. The maximum Courant number was set as 0.25 by the adjustable time step method to guarantee the accuracy. The volume fraction α was calculated via the multidimensional universal limiter with explicit solution method [16]. This method ensures that the volume fraction α is a value bounded by a limiter, similarly with the flux corrected transport technique [17].
In short, the simulation procedures during a time step are as follows: (1) Calculation of the time step based on Courant number; (2) Update of the dynamic meshes; (3) Update of the fluid properties (i.e., density ρ and dynamic viscosity μ); (4) Solving of Equation (3) by the MULES method to update the volume fraction function α; (5) Calculation of the velocity and pressure fields by the PISO algorithm.

Wave Generation Based on the Moving Boundary Method
In three-dimensional physical model tests, oblique waves and multi-directional waves are usually generated by a segmented wave-maker system. In this paper, to expand the effective area of numerically simulated multidirectional waves, the L-type wave-maker system was considered. The arrangement of a typical L-type segmented wave-maker system in a physical wave basin is shown in Figure 1a. The fixed wall, as generally taken in real laboratory wave basin, was set at the corner of the two rows of wave-makers to prevent the collision between them. Theoretically, ordinary segmented wave-makers on one side (for example just on the left side) are capable of generating oblique waves with large incident angles, but in fact, the wave generation capability is related to the target wavelengths and the width of paddles. Besides, the bigger the generated wave angle (θ in Figure 1) is, the smaller the effective test area is. Therefore in the laboratory, the two-sided segmented wave-maker system (L-type arrangement) is a feasible method to expand the effective test area. Since the wave-maker paddle can be replicated by the moving boundary, a straightforward way to simulate the three-dimensional waves is to replicate the two-sided wave-maker system shown in Figure 1a.
It is noteworthy that in the numerical model, the gaps between the paddle with fixed wall and adjacent paddles should be solid boundaries to avoid water leakage. In three-dimensional cases, the grid cells at these gaps will be stretched into parallelograms due to the phase difference between the adjacent paddles. Besides, the grid cells near the moving boundary are also stretched or compressed during the boundary movement. Therefore, in the numerical model, the mesh deformation should be controlled by artificial restraint to avoid the grid destruction (topological deformations) of the grid cells. In this work, the meshes are formed by the straight lines connecting adjacent grid points. Therefore, the mesh deformation can be restrained by controlling the grid point movements. The grid points include boundary mesh points and the internal mesh points, which can be controlled by different equations to generate target wave and avoid the grid destruction. A sketch map of the moving boundaries when simulating oblique waves by this model is shown in Figure 1b. The movements of the boundary mesh points should mimic the movements of a physical segmented multidirectional wave-maker system. The displacement or velocity of wave-maker paddles is required to control the boundary mesh point movements at each time step. The signals to control the boundary movement can be obtained from the displacement series of the physical wave-maker during the experiment. However, as stated above, the time steps in this model will be changed according to the calculated Courant number, so only the approximate displacement can be obtained by interpolation. In this paper, the displacements of the boundary mesh points are calculated at every step by theoretical formulas, which is more effective and accurate. The specific formulas for boundary mesh point movements for different kinds of waves will be discussed in Sections 3-5.
Further, a key issue for the boundary mesh movement is how to deal with the gaps and fixed wall in the physical wave basin. A straightforward approach is to thoroughly replicate the three-dimensional wave basin shown in Figure 1a. The transition regions for boundary mesh movement should be added at the gaps to avoid grid destruction. So, the grid point distribution on the boundary from point a to point b shown in Figure 1b will be treated as displayed in Figure 1c. The displacement of straight line b2-b3 and b4-b are identical to the wave-maker paddle movements, which can be calculated by wave generation theory. The movements of the grid points in the transition region b1-b2 and b3-b4 can be calculated by the cosine interpolation method [12]. However, the simulation error increases with the length of the transition region. Since the length of the transition region proposed by Higuera et al. [12] is proportional to the displacement difference of the adjacent paddles, this method will cause some deviations for three-dimensional wave simulation cases with large paddle displacement differences, for example, multi-directional irregular wave simulations. Further, the fixed wall at the corner also influences the simulation results. For instance, the oblique wave in the shadow area of Figure 1a will be deviated from the target wave.
In this work, the two-sided wave-makers are simulated as continuously deforming boundary to prevent the deviations due to the fixed wall and gaps, as shown in Figure 1d. The accurate displacements of the grid points on the moving boundary can be calculated via the wave-maker theory, so the grid lines on the wave generation boundary are formed by the straight lines joining two adjacent points. To generate oblique regular wave, the moving boundary on the x-axis will be a cosine curve of wavelength L' = L/cos θ, in which L is the wavelength of the target wave. The mesh convergence discussed in Section 3.2 shows that the grid number per wavelength is no less than 80, so the grid point number within the length L' is more than 80/cos θ, which ensures that the moving boundary is an almost completely smooth cosine curve. Under such conditions, every point on the smooth moving boundary can be seen as an infinitesimal wave-maker paddle. The deviations due to the fixed wall and the gaps can be avoided by this method.
The internal mesh points near the boundary meshes also moves during the wave generation to counteract the cell deformation due to the boundary movements. The governing equations of the internal mesh points should work as artificial restraint to avoid mesh destruction. A straightforward approach for artificial restraint is to treat every single grid cell as a solid element. If the deformation causes no destruction of the solid elements, then the meshes are still in good condition. The stress-strain equations can be used as the governing equations to control the deformation, however this will introduce new variables and increase simulation time. Some simplified approaches such as the Arbitrary Lagrangian-Eulerian method [18], the spring analogy method [19] or Laplace equations [20] can be used. In this work, the Laplace equation is adopted to control the mesh movements: where γ is the diffusion coefficient, up is the motion velocity of grid point, x is the position of mesh points and the superscripts n and o represent new and old time step, respectively. The diffusion coefficient γ can be defined as a field, whose values are the reciprocals of the distances with the nearest wave generation boundary. With this approach, the movements of the mesh points near the dynamic boundary will be larger than the internal mesh points. Take cell b-f-g-h in Figure 1d as an example, the cell moves to b'-f'-g'-h' at the next time step, as shown in Figure 2a. The movement of the internal mesh point is smaller than the boundary mesh point movement (∆ζg < ∆ζh, ∆ζf < ∆ζb). The longer the distance with the moving boundary, the smaller the displacement is. So, the movements of the boundary mesh points pass to the internal mesh points decreasingly, and the deformation of boundary grid cell can be balanced by the internal grid cells. Further, point b and point h replicate the movements of wave-maker, i.e., move perpendicular to the x-axis. Since the time step is very small under the control of Courant number, the displacement differences (∆ζ) between two time steps are small and the cell has no topological deformation. In two-dimensional wave simulations, the grid destruction can be avoided by the above methods. However, in three-dimensional cases, the only exception is the cell at the corner (cell a-c-d-e in Figure 1d), which movement sketch is displayed in Figure 2b. Since point a is the joint point of the two moving boundaries, its displacement is a vector ζa = (ζax, ζay). If ∆ζax ≥ ∆xac or ∆ζay ≥ ∆yae, the grid cell will be destroyed. In the numerical model, the displacement of a grid point is defined as the displacement difference with the initial position when t = 0 s. So, the above inequalities can be written as ζax ≥ in which θmax and θmin are the maximum and minimum interior angles of the quadrangle a-c-d-e in degrees. The aspect ratio is the ratio of the longest and shortest side lengths of the quadrangle. If the skewness is smaller than 0.8 and the aspect ratio is smaller than 40 [21], the grid cell is in good condition. Otherwise, the boundary mesh points should be redistributed by solving Laplace equations. For the boundary a-c, the grid points can be redistributed as: where upx n is the component of up n on x-axis, x is the coordinate of mesh points, the superscripts n and o represent new and old time step, respectively. The diffusion coefficient γ can be calculated as the reciprocals of the coordinate differences with point a on x-axis. It is noteworthy that the coordinates used here are the coordinates at the end of last time step. For instance, the diffusion coefficient of point c is calculated as: For the grid points on boundary a-b, the redistribution algorithm only changes the x-coordinate, while the displacements at new positions can be calculated by wave generation theory. The grid points on boundary a-e are also redistributed by applying this method on y-axis. After the boundary mesh point redistribution, the positions of internal mesh points can be updated by Equations (6) and (7). These procedures can be repeated several times until the grid cell a-c-d-e fulfills the above conditions. Generally, one cycle is enough to redistribute the grid points well, so the increased simulation time by this procedure is acceptable. The flow chart of the dynamic mesh algorithm is shown in Figure 3. In general, the movements of the boundary mesh points are calculated by theoretical formulas to mimic the physical wave-maker, while the movements of the internal mesh points are controlled by Equations (6) and (7) (6) and (7) Three-dimensional wave simulation? simulations, the grid cell at the joint point of the two moving boundaries should be checked at every time step to decide whether the redistribution procedures are needed. Since OpenFOAM is written by C++, a class-based object-oriented language, the moving boundary in this work is created as a derived class from a base class named 'fixedValuePointPatchField' in OpenFOAM code. The movements of the boundary mesh points can be calculated and assigned to the grid points by the new boundary. Even though the new boundary is capable of replicating the physical wave-makers, the computational cost also increases due to the calculation of dynamic meshes. The simulation time for two-dimensional wave cases increased about 12%-16% compared to waves2Foam code [6].

Wave Damping
To avoid wave reflections, a damping zone is set up near the outlets of the wave tank. In this study, the linear form damping coefficient is adopted based on the comparison results of Han [22]. The additional source term for wave damping in Equation (2) is denoted as: where β is the empirical coefficient and takes 3.0 in the current study, xs is the starting coordinate of the damping zone and xe is the end coordinate. Equation (13) is used for two-dimensional wave flumes, while for three-dimensional wave tanks, the expression changes to Equation (14):

Boundary Conditions
The basic numerical model being established, the looping algorithm can start running when the initial fields and boundary conditions are given. In the two-phase flow model, the appropriate boundary conditions should be assigned for the velocity field, the pressure field and the volume fraction field. The top boundary of the numerical wave tank is simulated as an open boundary. All the other boundaries can be seen as moving walls and static walls. OpenFOAM provides some common boundary conditions for users. In the present model, the boundaries of the volume fraction field are set as zero gradient boundaries. The pressure field uses the boundary called fixedFluxPressure, which guarantees that the flux at the boundary is consistent with the velocity boundary conditions to ensure mass conservation. To make sure the moving boundary works as a physical wave-maker, a boundary called movingWallVelocity is applied for the velocity boundary. The flux at this boundary is zero, so that the boundary movement can drive the fluid before it moving precisely.

Wave Generation
Piston-type and flap-type wave-makers are widely used in physical wave tanks to generate waves, and the present work focuses on the implementation of the piston-type wave-maker in a numerical model. The flap-type wave-maker can then also be simulated via a similar methodology. The relationship between the stroke and wave parameters can be obtained via the potential flow theory [23], written as: in which H is the wave height of the target regular wave in m, S is the stroke of the wave-maker paddles in m, k0 is wave number, d is water depth and T0 is the transfer function of progressive wave.
In accordance with physical wave-makers, the moving boundary can generate regular waves by sinusoidal motion: where ζ is the displacement of the moving boundary in m. During the first two wave periods of simulation, the displacement function should be multiplied by a damping coefficient δ to increase the stability of wave generation: in which t is the simulation time, t = 0 at the beginning of simulation and T is the wave period.

Mesh Convergence Study
Considering the numerical dissipation effect, it is important to choose the appropriate mesh size to simulate propagating waves. A two-dimensional numerical wave flume was established to perform the mesh convergence study. The computation domain is divided into hexahedral meshes, as shown in Figure 4. The wave flume is 40 m long, for a water depth d = 0.6 m. The mesh size along the x-axis is uniform, while that along the z-axis is refined near the water surface. A 9 m long (about 3 times the wavelength) damping zone is set at the end side to absorb the incoming waves.  Figure 5 shows the numerical results of the water surface profiles simulated by different grid dimensions. The target regular wave height is H = 0.04 m, and wave period is T = 1.5 s. The wave energy was absorbed thoroughly in the damping zone. The label X40Z10, as an example, represents the number of mesh per wavelength (n.p.w.l. for short), 40 in this case; and the number of mesh per wave height (n.p.w.h. for short) is 10 here. The spatial resolution has an obviously influence on the stability of the generated wave profiles, as can be seen from Figure 5. With a fixed n.p.w.h. of 10, obvious wave height attenuation can be observed when the n.p.w.l. is 40. However, when the n.p.w.l. increases to 60 and 80, almost no attenuation remains. Similar results to X80Z10 are obtained if the n.p.w.l. increases to 100 or 120 (not shown). In terms of vertical resolution, at a fixed n.p.w.l. of 80, steady results are obtained for an n.p.w.h. greater than 6. Obviously, the n.p.w.l. affects the stability of wave height while n.p.w.h. affects both the wave height and phase. In general, the grid dimensions necessary to simulate regular waves should meet the following requirements: the n.p.w.l. should be no less than 80, and the n.p.w.h. no less than 6. Figure 5 also shows that residual attenuation still exists even with sufficient mesh numbers. This attenuation is caused by numerical dissipation and is inevitable even with high order numerical schemes. Further, some non-linear features appeared in the generated wave profiles. As shown in Figure  5, the wave crests and wave troughs were asymmetric, similarly with high order Stokes waves. This characteristic was consistent with the waves generated in the laboratory, which makes the moving boundary method a better choice to replicate the physical wave flume in a numerical model. In Figure 6, the numerical results for four different waves were compared with the analytical solutions at the positions 10 m, 10 m, 10 m and 20 m away from the wave generation boundary, respectively. The wave parameters and mesh sizes are listed in Table 1 for the four cases. The wave simulated in case 1 was the deep water wave and case 2, 3 and 4 were the waves in intermediate depths. The water steepness (H/L) of case 3 was 0.08, which aimed to test the ability to generate and absorb the wave with large wave steepness. The analytical solutions shown in Figure 6 were calculated by fifth-order Stokes wave theory [24] since the wave parameters were in the range of Stokes wave [25] and wavelengths were shorter than ten times the water depth [24]. The root-mean-square error [26] between the numerical results and analytical solutions during the last fifteen wave periods of the four cases are listed in Table 2. In general, the numerical results are shown to be in good agreement with the theoretical solution. The biggest root-mean-square error and relatively significant attenuation occurred in the large wave steepness case (case 3) because of the inevitable numerical dissipation. Besides, the deviations between the numerical results and theoretical solutions were also due to the influences of spurious free-waves with first order wave generation theory. Since this work aimed to replicate the ordinary physical wave-maker, the second-order wave-maker theory will be included in the future work.     Table 3 shows the reflection coefficients of right boundary for the four cases, which were calculated by the method proposed by Goda and Suzuki [27]. The lengths of the damping zone in these cases were three times the corresponding wavelength. The reflection coefficients were under or about 10%, which indicates that the wave damping zone was effective, even for short wave and steep wave. The biggest value of reflection coefficient for case 3 took place due to its biggest wave steepness and nonlinearity. Another feasible method was to absorb the waves by a moving boundary at the end of the wave flume. The movement of the absorption boundary can be controlled by the equation proposed by Schaffer and Klopman [28]. Even though this method can reduce the computational cost (no need to extend the computational domain as damping zone), the absorption efficiency decreases obviously for short waves, as tested by Higuera et al. [12]. Besides, for three-dimensional cases, this method is impracticable due to the wave velocity component parallel to the moving absorption boundary. Since the three-dimensional waves are also included in this work, the damping zone method was adopted. Except for the short wave cases, the reflection coefficients listed in Table 3 were also slightly lower than the results by the moving absorption boundary method [12].

Active Wave Absorption
A secondary reflected wave occurs when an incident wave is reflected by the structure and then reflected again by the wave-maker paddles. For a long term physical model test, for example, irregular wave interaction with structures, the multiple reflections between structure and wave-maker paddles will make experimental results unreasonable. An effective mitigation method is the adaptive wave absorption, i.e., the correction of the paddle positions at every time step to avoid the occurrence of secondary reflected waves. The position correction can be calculated by different feedback signals, including the surface elevation at or before the paddles and the wave force on the paddles. This method is listed among the different approaches for wave absorption delineated by Schaffer and Klopman [28]. Even though these methods were originally developed for physical model testing, the absorbing wave-makers can be mimicked by the moving boundary method in our numerical model.
In addition, the influence of secondary reflected waves is ignored at times in physical experiments, especially for the structures with high transmission coefficients. Therefore, the active wave absorption function in the wave generation boundary is controlled by a switch. The correction algorithm for paddle positions will only be activated when the switch is turned on.
In the present work, the capability of active absorption by the moving boundary method is discussed. The piston-type absorbing wave-maker with a front-mounted wave gauge was proposed by Hirakuchi et al. [29] and validated in the laboratory by Liu et al. [30]. The velocity of the absorbing wave-maker can be expressed as: in which ζ is the displacement of the moving boundary, ω is the angular frequency of the target regular wave, ηp and η are the target and measured surface elevations at the wave paddle, respectively. The parameter D is calculated by: 1 n n D T ∞ = =  (19) in which Tn is the transfer function of the evanescent mode: where kn is the wave number of the evanescent mode, also called an imaginary wave number. It can be calculated by the dispersion relation for the evanescent mode: The value of kn is in the interval from 0.5π (2n + 1)/d to 0.5π (2n + 3)/d and can be calculated in the program via the binary search method. Zhang et al. [31] gives the relationship between the sums of the evanescent-mode transfer functions Tn with different n. The results show that the curves of Tn match well with each other when n = 20 and n = 200, especially under the condition that Therefore, n is set to 20 in the present work.
Equation (18) listed above is a differential equation for the unknown variable ζ. Considering the very short time step in the numerical model, the term ζ(t) on the right side of Equation (18) can be replaced by the displacement in the last time step. In addition, the differential can be simplified as difference, so the final expressions become: n o ζ ζ ζ = +Δ (23) in which ∆t is the time step in s.
To validate the performance of the active wave absorption algorithm, a 9.01 m-long wave flume is established without damping zone. The incident wave height is H = 0.04 m, the wave period is T = 1.5 s, and the water depth in the flume is d = 0.6 m. Two cases are simulated with a total reflection boundary at the end of the flume. First, using the normal regular wave generation method (no adaptive absorption) and second, using the adaptive wave absorption method. Figure 7 shows the time histories of surface elevations at x = 9.0 m for the two cases. With the normal wave generation method, the water surface is unstable due to the superposition of incident waves, reflected waves and secondary reflected waves. When the active absorption is activated, the surface fluctuation becomes stable, but the mean water level gradually declines. This is because the wave flume is stretched due to the slow drift movement of the moving boundary. The slow drift comes from the zero frequency components in the displacement signals [31], and high-pass filtering can be used to circumvent this issue. In this work, the moving boundary displacement was corrected by an additional slow uniform motion, that is, to push the boundary to the initial position slowly. Figure 8 shows the displacements of moving boundaries before and after the low frequency correction, while the comparison between the simulated results after drift correction and the analytical solution are shown in Figure 9. The analytical solution is the superposition of the incident wave and reflected wave, under the assumption that the secondary reflection is totally absorbed. This result demonstrates that a stable standing wave is formed in the flume and the active absorption method is effective in the numerical model.

Wave Generation
The generation of irregular waves can be achieved by the superposition of a series of Airy waves with different wave periods and random initial phases: where ai, ki, ωi and εi are the amplitude, wave number, angular frequency and initial phase of the ith wave component, respectively. εi is a random number uniformly distributed in the 0-2π range. This method, also called random phase spectrum method, can generate irregular waves based on the target spectrum. In the following calculation, the JONSWAP spectrum suggested by Goda [32] is adopted: in which f is the frequency, H1/3 and T1/3 are the significant wave height and significant wave period, Tp and fp are the spectral peak period and spectral peak frequency, and γ = 3.3 is the peak enhancement factor.
Since most of the spectral energy is concentrated in a small frequency interval, the cut-off frequencies in this work are chosen as: The frequency range is divided equally into N frequency intervals: Note that the time history of surface elevation will repeat with the period 2π/∆ω. Therefore, the frequency number N is set as 500 to avoid repetitions in the normal simulation time-scale. The bandwidth of the ith frequency interval can also be calculated by: In each frequency interval, the midpoint is chosen as the representative frequency, i.e., The corresponding wave amplitude is: The displacement of the moving boundary can then be expressed as: Since the surface of irregular waves is much more complex than that of regular waves, the grid used to simulate irregular waves should be denser due to the high frequency components of the irregular waves. The mesh convergence study shows that the number of meshes per wavelength (Lp) should be no less than 150 and the number of meshes per significant wave height (HS) no less than 10. The wavelength used here is the one corresponding to the spectral peak period Tp.
To validate the irregular wave generation performance by the moving boundary method, the target irregular wave is simulated in a 30 m-long wave flume similar to the flume shown in Figure 4. The water depth was 0.55 m, the significant wave height was Hs = 0.04 m and the spectral peak period was Tp = 1.5 s. A 9 m-long (about three times that of Lp) damping zone was set at the end of the flume for wave absorption. The grid sizes used in this case were Lp/∆x = 200, and Hs/∆y = 15. The irregular wave was simulated from 0 to 190 s, and 8192 data points were collected at a sampling rate of 50 Hz (from 26.18 to 190 s). The wave spectra at the positions 5 m, 10 m, 15 m and 20 m from the moving boundary were analyzed and the comparison between the analyzed spectra and the targets is shown in Figure 10. The figure shows that the simulated wave spectra were in good agreement with the target spectrum. Further, the wave heights along the wave flume were analyzed and non-dimensionalized with the inputted significant wave height. The comparison of the simulated non-dimensionalized significant wave heights and the theoretically calculated results is given in Figure 11. The theoretical solutions used here were calculated by the same superposition method described in Equation (24), and with the same parameters as for the calculation of wave generation boundary. Although the numerically calculated wave height tended to decrease and some differences existed between the numerical results and theoretical solutions due to numerical dissipation, the errors between the numerical results and the target wave height were less than 5% and deemed acceptable. It should be mentioned that the numerically calculated wave height increased again several times along the flume. This might be caused by the non-homogeneity in space existing in the superposition method because of inexact randomly uniformity of the initial phases in the range of [0, 2π].
To further analyze the sources of error, the time histories of the simulated surface elevations at the positions 1 m, 10 m and 20 m from the wave generating boundary were compared with the theoretical solutions in Figure 12. The corresponding root-mean-square error between numerical results and analytical solutions during 26.18-190 s at these positions are listed in Table 4. The simulated result at the position 1 m away from the moving boundary was almost exactly the same as the theoretical result, which indicates that the moving boundary worked well. This result also indicates that the damping zone method was effective for the irregular wave absorption with a long simulation time. A small difference between the numerical and the theoretical result appeared at x = 10 m, and the differences became more obvious at farther locations. These differences might be attributed to the numerical dissipation and the grid size limitation in the simulations. The results provided in Section 3.2 indicate that the selection of grid size should be based on the wave parameters. Obvious wave height attenuation occurred in the cases with less than 60 cells per wavelength. However, only one grid size could be applied for the irregular wave simulation although many different wave components with different wavelengths exist. We could infer that the high frequency wave components will decay along the wave flume, which then results in the increased differences related to the distance from the moving boundary. As the theoretical wave surface was also calculated by linear superposition, the nonlinear wave generated during the simulated wave propagation can be another reason for the errors. However, the errors of the statistical wave parameters were within an acceptable limit even at positions relatively far away from the moving boundary, and they could meet the requirement for the numerical simulation.    Position (m) 1 10 20 Root-mean-square error 0.00092 0.00371 0.00491

Active Wave Absorption
As stated above, an irregular wave is generated from the superposition of regular waves. The active wave absorption algorithm for irregular wave can also be achieved via a similar method. In Equation (22), the parameter D should be determined first. However, as it is related to the wave frequency, the representative frequency based weighted average method proposed by Liu et al. [30] was adopted in the present model. The movement of the moving boundary is directly given as: (2 ) (36) in which N is the number of the representative frequencies, ωi is the representative frequency of the ith component, T0i and Di are the corresponding parameters for ωi and λi is the weight coefficient calculated as: in which Si is the value of spectral density at frequency ωi.
Similarly to Section 3.3, a 20.1 m-long wave flume without the damping zone was set up to generate irregular waves by both the normal method (without active absorption) and the adaptive absorption method. The water depth was 0.55 m. The spectral peak period was Tp = 1.5 s and the significant wave height was Hs = 0.04 m. The comparison of the surface elevations at the position x = 1.0 m is shown in Figure 13, while Figure 14 shows the displacement of moving boundary for both cases. Obviously, the surface elevations for the absorption case were much smaller than those of the normal case ( Figure 13). The slow drift of the moving boundary ( Figure 14) only occurred at the initial phase of the correction by the adaptive absorption, from about 20 s (when the first reflected wave reached the moving boundary) to 40 s. This is because when the active absorption switches on, the difference calculation for the reviewed position of the wave generator causes the slow drift due to the reflected wave front. However the following reflected waves become stable, thusly the drift keeps almost unchanged. The total drift was less than 0.02 m, which had no significant influence on the simulation. Figure 15 shows the wave spectra of the separated incident wave obtained by the separation method proposed by Goda and Suzuki [27] according to the measured waves for the two cases. The separated incident wave spectrum in the active absorption case was in good agreement with the target spectrum, which means that the active absorption system was effective in absorbing the reflective wave from the flume.

Oblique Regular Wave Generation
As discussed in Section 2.3, the two-sided wave-maker system (L-type wave-maker) was replicated in the numerical model to expand the effective wave area. The two-dimensional model was extended to mimic a three-dimensional physical wave tank with wave-makers on two sides of the tank and damping zones on the other two sides. The layout of the numerical wave tank is shown in Figure 16.
Theoretically, the oblique wave surface can be written as: in which H is the target wave height, ω is the angular frequency, and k is the wave number. θ is the angle between the incident wave direction and the x-axis.
For the two-sided wave-maker system (L-type wave-maker system), the signal of the moving boundary along the x-axis is: The signal of the moving boundary along the y-axis is: In the above equations, T0 is the 2D transfer function given by Equation (15).  Figure 17 presents the water surface profiles at t = 30 s and 35 s. The simulated wave surface was distributed unevenly in space. This spatial non-uniformity is a key feature of oblique regular waves generated by segmented wave-makers, as reported by Takayama's [33] theoretical analysis. It is mainly attributed to the finite length of the totally wave-maker system and it will cause a fluctuation along the wave maker due to the two ends of the wave maker system. However, the continuous approach versus the segmented approach [12] takes the advantage that it can avoid the effect of the width of the segmented wave board as pointed out by Takayama [33]. Reasonably, the waves generated by L-type wave-maker distributed more evenly than that generated by one-side wave-maker, as discussed by Zhou [34] with the same method as Takayama [33]. The ratios of the calculated wave heights and the target wave height along the lines y = 3 m, y = 6 m and y = 9 m are provided in Figure 18 and show some along-section fluctuations. The trends of the wave height variations were in accordance with the theoretical results of Takayama [33] and Zhou [34].   It should be pointed out, however, that outside the damping zone, the effective simulation area could be expanded to the whole wave tank for oblique wave simulations due to the L-type arrangement of the wave-maker system. For example, the midpoint of the computational domain (6 m, 6 m) will be out of the effective simulation area if the oblique wave is generated only by the moving boundary on the x-axis. Figure 19 shows the comparisons between the time series of surface elevations at the position (6 m, 6 m) and the linear theoretical solution. Good agreement can be observed from the figure.

Multidirectional Random Wave Generation
The waves observed at sea are multidirectional random waves. The directional spectrum of the multidirectional waves can be written as the product of the frequency spectrum and the directional spreading function, that is: in which the frequency spectrum S(f) is calculated by Equation (25). The directional spreading function G(f, θ) should satisfy the following condition: in which [θmin, θmax] is the wave direction distribution range. The Mitsuyasu-type spreading function was adopted in this paper: where s is the directional spreading parameter, and θ0 is the major wave direction. The function G0(s) can therefore be solved as follows: max min 1 2 0 In the current work, multidirectional random waves are generated by the single direction per frequency method proposed by Yu et al. [35]. The frequency range and the direction range can be divided into I and J parts, respectively. The multidirectional random wave can be calculated by the superposition of I × J wave components: Similarly with the method explained in Section 4.1, the frequency range is divided into I parts equally by Equations (29)- (33). The difference here lies in the representative frequency, which should be calculated as: in which RANij are the random numbers distributed evenly along the interval [0, 1]. By this method, the representative frequencies in a frequency interval are made to be different, i.e., there is a one-to-one correspondence between the representative frequency and the direction. The artificial phase locking phenomenon caused by the double summation method can therefore be prevented. In the numerical model, the two-sided wave-makers are replicated to expand the effective simulation area. The signal of the moving boundary on the x-axis is: and the signal of the moving boundary on the y-axis is: A case of multidirectional irregular wave was simulated in a three-dimensional numerical wave tank, which was 20.0 m long, 20.0 m wide and 0.55 m deep as shown in Figure 17, while the effective simulation area was reduced to 6.0 m by 6.0 m, comprehensively considering the simulated long time duration for multidirectional waves as follows and the endured run time for the case. The significant wave height of the target wave was Hs = 0.04 m and the significant wave period was Ts = 1.5 s. The major wave direction was θ0 = 45°, the wave direction distribution range was [0°, 90°] and the directional spreading parameter was s = 25. In the effective simulation area, the grid sizes were ∆x = ∆y = 0.021 m and in the area near water surface, ∆z = 0.004 m. The total grid number was 7.72 million and this case (from 0 to 190 s) took 78 days to run by nine parallel processes on a workstation powered by an Intel Xeon E5-2690v2 CPU. Figure 20 displays the water surface profiles at t = 90 s and 180 s. Typical short-crested wave surface could be observed. With the same rationale as for the two-sided wave-makers, the effective simulation area was expanded to the whole 6 m by 6 m area. Further, the directional spectrum can be analyzed by the Bayesian approach [36] with a wave gauge array, which consists of five wave gauges as shown in Figure 17. The coordinates of these gauges are (3. Figure 21 shows the time series of surface elevation measured by gauge (3.253 m, 3.244 m). The comparison between the simulated wave directional spectrum and the target directional spectrum is shown in Figure 22. It can be seen that they were basically identical to each other. Nonlinear wave components could also be observed from the analyzed results. Furthermore, the comparison between the simulated and the target frequency spectra and directional functions (corresponding to the peak frequency fp) are shown in Figure 23a,b, respectively. In general, the simulation results were in good agreement with the target values.

Conclusions and Remarks
A three-dimensional numerical wave tank with an L-type wave-maker system was developed based on RANS equations. The piston-type wave-makers were mimicked by the moving boundary method to generate water waves in this model. The advantage of this numerical model is the complete replication of the setting of a laboratory physical wave tank.
The displacement of the moving boundaries was calculated by theoretical method to support the adjustable time step. The two-dimensional regular waves generated in this model were in good agreement with theoretical solutions. As for the irregular wave simulation, the error increased with the propagation distance, although the accuracy still met the requirements of numerical simulation. Furthermore, the active wave absorption algorithm could be activated to prevent the influence of secondary reflection for the two dimensional regular and irregular waves.
The two-sided wave-maker system (L-type arrangement) was adopted in this model to expand the effective wave areas for three-dimensional wave simulations. The wave-makers were simulated as smooth boundaries to prevent the deviations due to the gaps between adjacent paddles. Both the oblique regular wave and multidirectional random wave were generated with designated parameters. The features of the oblique regular waves and multidirectional waves simulated by this model are consistent with their respective theoretical solution, which further illustrates the "replication" efficiency of laboratory wave basin conditions. In particular, the L-type wave-maker system can make the simulated oblique and multidirectional waves effective in the whole basin. This model could thus be applied to reveal or supply finer details or fill up the missing data, which were unavailable due to limitations of the physical instrument performance and setting in a laboratory basin. The disadvantage of this method is the increased computational cost due to the dynamic mesh calculations compared to the relaxation zone method by Jacobsen et al. [6].
While this paper focused on the wave generation by piston-type wave-maker, flap-type and other types of wave-makers could be implemented with similar methods in a future work. The fully three-dimensional active absorption for physical wave-makers is still an open field. Furthermore, as shown in Section 5, the simulation of three-dimensional cases is still computationally expensive, which explains why three-dimensional active absorption simulations were not tested in this work. Attempts at improving the simulation performance will be provided in subsequent works.