# Generating Water Waves on Steady Currents Using Mass Source Wave Maker Coupled with Analytical Relaxation Wave Absorber

Authors: Xing Feng; Jia Liu; Ruina Ma
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In order to numerically simulate the wave-current interaction problems frequently encountered by aquaculture structures, a two-dimensional numerical wave-current tank model was established here based on a mass source wave maker coupled with an analytical relaxation wave absorber. The wave-maker model and the wave-absorber model were embedded into a two-dimensional RANS solver, which was closed with RSM turbulence scheme. The volume of fluid (VOF) method was adapted to accurately capture the free surface between water and air. To generate a steady uniform current flow, the uniform current flow velocity was calculated at the left-hand-side (LHS) and right-had-side (RHS) outflow boundaries, respectively. Once the steady uniform current flow was generated over the whole computational domain, the target water wave was marked within a specified region by embedding the mass source function based on wave theory into the mass conservation equation and then propagated on the generated uniform current flow. To verify the accuracy of the numerical wave-current tank established here, some of the obtained numerical results were then compared with the experimental results and the analytical solutions, and they agreed well with each other, indicating that the model developed here has great ability in simulating water waves on uniform currents over constant water depth. The established numerical wave-current tank was then used to study the optimal layout of the mass source region and the effects of water current velocity on water surface wave parameters during regular wave coupling with uniform water currents. Meanwhile, the established model was extended to generate steep wave and apply in deep water conditions. Finally, the proposed methods were applied to investigate the wave-current-structure interaction problems and the propagation of solitary waves traveling with coplanar/counter currents. Model-data comparisons show that the developed model here is potentially useful and efficient for investigating the inevitable wave-current-structure interaction problems in aquaculture technologies.


---

## Body

## 1. Introduction

Recently, marine aquaculture, also known as aquafarming or fish farming, research and technologies are developing rapidly. In order to ensure the safety of such activities and the corresponding structures, it is important to accurately calculate the ambient flow field when making designs, and for the marine environments the ambient flow field is all the more important as maritime activities take place where the water waves and ocean currents always coexist. The interaction of waves and currents, which can significantly alter the wave parameters and induce wave breaking and scattering, plays an important role in response to floating structures and other aspects related to the development of marine science and technology. The interaction between wave and current, and its effect on marine structures has attracted the attention of ocean and coastal engineers and researchers, and the mechanisms have been extensively investigated in the past. Among these, many physical laboratory experiments have been carried out to examine the wave-current interaction problem and its effect on marine structures. Umeyama carried out a series of physical experiments on the wave-current interactions in a laboratory flume with a constant water depth. The turbulence properties of pure-current conditions, pure-wave conditions, and waves traveling with following/opposing currents conditions were measured, analyzed, and compared (by varying the wave height and/or wave period in different uniform current flow fields) [1]. Toffoli and Waseda present two independent sets of physical experiments performed in the experimental wave tank of Plymouth University and the Ocean Engineering Tank of the University of Tokyo, which measure and analyze the flow field properties of rogue waves in an opposite current in detail. These two independent sets of laboratory experiments confirmed a recent conjecture using a current-modified cubic nonlinear Schrodinger (NLS) equation, which establishes that a stable wave traveling with a counter current may trigger rogue waves [2]. Kristiansen and Faltinsen experimentally measured and analyzed the mooring loads on an aquaculture net cage considering the wave-current interactions [3]. Tambroni and Figueiredo da Silva performed a lot of experimental investigations in the Total Environment Simulator (TES) recirculating flume belonging to the University of Hull (UK) for determining the impacts of macroalgal mats of Ulva intestinalis on flow resistance and bed stability considering wave-current interactions when waves and currents coexisted [4]. de Jesus Henriques et al. experimentally measured and analyzed the resultant effects of the wave-current interactions on the power and thrust performance of a model-scale horizontal marine current turbine by employing the high-speed recirculating water channel at the university of Liverpool [5].Although experimental investigation can be straightforward and accurate, physical experiments are usually utilized for calibration cases because of the high cost of establishing the physical wave tank and flow field measuring, inconvenience for maintenance and management, and difficulties in monitoring some specific flow data (may lead to the ignorance of some important details). The numerical investigations based on the advanced computer technologies, efficacious numerical methods, and high-level Computational Fluid Dynamics theory have emerged and sprung forth. Numerical methods that can faithfully reproduce the experimental process have become increasingly useful and powerful for measuring and analyzing the intricate flow fields, especially when water wave and current coexists. In recent years, the numerical application of Navier-Stokes or RANS equations model for analyzing wave-current interaction problems have been explored to a great extent. Compared with other models, the Navier-Stokes or RANS equations model can directly solve the coupled fields of pressure and velocity, which can show the wave pattern and the current state simultaneously. Sung-Yong Kim carried out numerical calculation of the hydrodynamic load of fixed offshore substructures in waves and uniform currents based on the numerical wave tank technique, which was established by using the commercial CFD software-FLUENT [6]. For generating the waves on the uniform currents, the static boundary wave-maker method was used. As pointed out by many researchers, this kind of method will unavoidably bring evident numerical dissipation due to the simultaneous specification of wave and current at boundaries. Besides this kind of wave generation methods, there are moving boundary wave-maker and internal wave-maker methods. For the former method, the moving boundary coupled with dynamic mesh technology is employed to numerically simulate the motion of a realistic physical flap-type wave-maker or piston-type wave-maker. For the posteriori method, water wave generation, including wave absorption, is achieved by designating a source function deduced based on the target wave theory within the specified region inside the computational domain.There are plenty of applications of both the methods. Based on the moving boundary wave-maker method, Wood et al. numerically analyzed the run-up of steep nonbreaking waves by utilizing the commercial CFD software-FLUENT [7], Finnegan and Goggins successfully simulated the linear water waves in both deep water and water with finite depth and developed the model to investigate wave-structure interactions by the application of the commercial CFD software-Ansys CFX [8], Anbarsooz et al. numerically investigated the intricate propagation physics of the fully nonlinear viscous waves by utilizing both the piston-type and the flap-type moving boundary wave-makers [9], Telesetal et al. developed a second-order, piston-type wave-maker to simulate wave-current interactions at a local scale by employing a RANS-equations-based CFD solver [10], Markus et al. established a numerical wave channel by adopting a CFD solver for simulating wave-current interactions focusing on evaluating the flow field of a nonlinear wave combining with a nonuniform current [11]. In all these models with the utilization of a moving boundary wave-maker, the computational domain varies as the moving boundary (solid body) steps toward or away from the fluid domain, making a re-meshing technology inevitable. During the numerical calculation, the computational grids should be updated in each calculation time step usually leading to a large grid distortion, which will make simulation unstable, inaccurate, and inefficacious. For the internal wave-maker method, there are mainly three subbranches including analytical relaxation method, mass source function method, and momentum source function method. Xing and Wu adopted the analytical relaxation method proposed by Madsen et al. to generate regular waves on uniform currents over constant finite water depth based on the RANS equations model [12, 13]. As pointed out by many researchers, the analytical relaxation method has low computational efficiency due to the maximum time step limitation. The mass source function method and the momentum source function method were first adopted and embedded in Boussinesq-type equations by Wei et al. and were developed by others to be applied based on Navier-Stokes equations or RANS equations [14, 15]. As firstly reported by Choi and Yoon [16], the momentum source function method based on the RANS equation model was adopted to generate various water waves such as linear waves, Stokes waves, solitary waves, and random waves in finite water depth. As reported by Chen and Hsiao, it is difficult to generate stable water wave profile with deep water depth by employing the momentum source function as it is derived from linear Boussinesq equations under the shallow water assumption [17]. Meanwhile, the momentum source function method achieved by adding momentum source items in the momentum equations makes it also challenging to generate waves on currents. The application of mass source function method based on the Navier-Stokes or RANS equations model was firstly reported by Lin and Liu with the utilization of the finite difference method (FDM) and VOF algorithm on a staggered computational grid system [18]. Recently, this method has been applied to a great extent, as shown by Vasarmidis et al. [19]. This method is advantageous compared to the static boundary, especially when waves are reflected back to the numerical boundary (wave-structure interaction problems), and the rules of thumb presented by Wu and Wu for specifying the source region have become a general and useful principle [20, 21]. Zhang et al. adopted this method based on a RANS solver and the rules of thumb for specifying the source region in investigating the wave-current interactions with the employment of the VOF algorithm and the finite volume method (FVM) [22]. Nevertheless, it is essential to carry out a lot of numerical tests with the mass source function method for designing the source region for generating a particularly designated wave. Although the mass source function method does not have an inherent limitation on the relative water depth and can be applied to the wave-current interaction simulation, the simulation of the regular and/or irregular wave propagation on uniform current over constant water depth from shallow water conditions to deep water conditions based on Navier-Stokes equations or RANS equations is still scarce. Therefore, an approach for designing the mass source region for generating waves on uniform current over constant water depth from shallow water conditions to deep water conditions should be developed and applied in solving the wave-current-structure interaction problems.In this study, the mass source function method and the analytical relaxation wave absorption method were embedded into a two-dimensional RANS equations model for generating regular waves on uniform currents over constant water depth from shallow water conditions to deep water conditions. The licensed FLUENT software, which was owned by the research center of Marine Safety and Pollution Prevention (MSPP) at Dalian Maritime University, was chosen as the base solver due to its popularity and applicability to numerical calculation of water wave problems. All the cases presented in the following of this context were run in parallel on series of Intel XEON E5-2640V3 processors (2.60 GHz). To verify the accuracy of the model proposed and developed here, the laboratory measurements by Umeyama and pertinent analytical solutions were employed as comparisons [23]. An optimal layout for designing the mass source region was brought out. Finally, the mass source function method and its designing approach proposed here are applied to generate steep wave and investigate the influence of the current velocity on the wave parameters of the regular and solitary waves from shallow water to deep water conditions and simulate the wave-current-structure interactions.


## 2. Mathematical Models

### 2.1. Governing Equations

In order to numerically simulate the water wave propagation process categorized as incompressible fluid motion, the following mass conservation equation, momentum conservation equation, and RANS equation are employed.(1)∂ui∂xi=Sm,(2)∂ρui∂t+∂ρuiuj∂xj=∂∂xjμ∂ui∂xj+∂uj∂xi−∂p∂xi+∂∂xj−ρui′uj′¯+ρgi+Si,where t represents time and ρ represents the fluid density, xirepresents the abscissa and ordinate, the subscripts i and j represent the indices used in Einstein’s sum convention, ui denotes the ensemble mean velocity component, p represents the fluid pressure, μ represents the dynamic viscosity, and g denotes the gravitational acceleration. Smrepresents the additional mass source term specified for target wave generation, Si represents the additional momentum source term designated for wave absorption. −ρui′uj′¯ denotes the Reynolds stress term. The solution about the Reynolds stress term can be found in the Manual Guide of FLUENT in detail [24].The VOF algorithm is adopted to track the free surface between the two fluids, namely, water and air. In VOF algorithm, which employs a volume function to represent the fractional volume of each fluid in the computational domain, the main principle can be expressed as follows:(3)∂α∂t+∂uiα∂xi=0.α=1 indicates that the computational cell is fully occupied by water phase, while α=0 corresponds to the computational cell full of air phase. Computational cells with volume function value of 0<α<1 represent cells containing both phases and a free surface.


### 2.2. Wave Generation by Employing the Mass Source Function Method

The mass source function method based on the RANS equations model developed by Lin and Liu was adopted in the present study with the utilization of FVM discretization scheme, which protrudes by avoiding generation evident numerical dissipation (particularly commonly detected with the simultaneous specification of wave and current at inflow boundary) and disturbing the flow field [25, 26]. Within the specified source function region, the mass source function deduced based on the corresponding wave theory was added into the mass conservation equation (1):(4)Sm=2CAη.in whichC is the water wave phase velocity of the target wave. A is the area of the source region. η is the water wave surface elevation; for linear monochromatic wave, it can be represented as η=asinωt, and for solitary wave, it can be represented as the following formula:(5)ηx,t=Hsech23H4h3x−Ct.


### 2.3. Wave Absorption by Employing the Momentum Source Function Method

The analytical relaxation method proposed by Madsen et al. has been widely adopted by many numerical hydrodynamics researchers for establishing NWT (numerical wave tank) in pure-wave conditions. Here, the analytical relaxation method was applied for wave absorption in wave and current conditions and the elimination of wave reflection without varying the water depth and the current velocity at outflow boundary. The basic mechanism of the application of analytical relaxation method for wave absorption can be presented as following: Within the specified relaxation domain (usually with the length of one to two wavelengths), the water particle velocity and pressure will be updated at each calculational time step by the added momentum source function in the momentum equations to approach the target current flow velocity and the corresponding pressure. The relevant relaxation algorithm can be expressed as:(6)uM=CuJ+1−CU,vM=CvJ,pM=CpJ,where the subscript M represents the renewed physical quantity in the designated domain—the wave absorption domain, the subscript J represents the computed physical quantity at the previous calculational time step, U denotes the target current velocity, C=Cx represents the relaxation function, Cxmin=1, Cxmax=0, and more details about the analytical relaxation method and the relaxation function can be found in Madsen et al. [13].Then, the momentum source function adopted here for the analytical relaxation wave absorption can be reversely deduced from the Euler equations ignoring the water viscosity by substituting the aforementioned relaxation algorithm. The different discretization forms of the corresponding momentum equations with and without the additive momentum source can be expressed as the following:(7)ρuMn+1−uJnΔt+ρuJ∂uJ∂x+ρvJ∂uJ∂z=−∂pJ∂x+Sx,(8)ρvMn+1−vJnΔt+ρuJ∂vJ∂x+ρvJ∂vJ∂z=−∂pJ∂z+Sz−ρg,(9)ρuMn+1−uMnΔt+ρuM∂uM∂x+ρvM∂uM∂z=−∂pM∂x,(10)ρvMn+1−vMnΔt+ρuM∂vM∂x+ρvM∂vM∂z=−∂pM∂z−ρg.By subtracting equation (7) from (9) and equation (8) from (10), the corresponding momentum source function in abscissa and ordinate can be deduced as(11)Sx=1−C∂pJ∂x−ρuJ−UΔt+ρ1−C2uJ∂uJ∂x+vJ∂uJ∂z−ρ∗C∗1−C∗U∗uJ,(12)Sz=1−C∂pJ∂z−ρvJΔt+ρ1−C2uJ∂vJ∂x+vJ∂vJ∂z−ρ∗C∗1−C∗U∗vJ.


## 2.1. Governing Equations

In order to numerically simulate the water wave propagation process categorized as incompressible fluid motion, the following mass conservation equation, momentum conservation equation, and RANS equation are employed.(1)∂ui∂xi=Sm,(2)∂ρui∂t+∂ρuiuj∂xj=∂∂xjμ∂ui∂xj+∂uj∂xi−∂p∂xi+∂∂xj−ρui′uj′¯+ρgi+Si,where t represents time and ρ represents the fluid density, xirepresents the abscissa and ordinate, the subscripts i and j represent the indices used in Einstein’s sum convention, ui denotes the ensemble mean velocity component, p represents the fluid pressure, μ represents the dynamic viscosity, and g denotes the gravitational acceleration. Smrepresents the additional mass source term specified for target wave generation, Si represents the additional momentum source term designated for wave absorption. −ρui′uj′¯ denotes the Reynolds stress term. The solution about the Reynolds stress term can be found in the Manual Guide of FLUENT in detail [24].The VOF algorithm is adopted to track the free surface between the two fluids, namely, water and air. In VOF algorithm, which employs a volume function to represent the fractional volume of each fluid in the computational domain, the main principle can be expressed as follows:(3)∂α∂t+∂uiα∂xi=0.α=1 indicates that the computational cell is fully occupied by water phase, while α=0 corresponds to the computational cell full of air phase. Computational cells with volume function value of 0<α<1 represent cells containing both phases and a free surface.


## 2.2. Wave Generation by Employing the Mass Source Function Method

The mass source function method based on the RANS equations model developed by Lin and Liu was adopted in the present study with the utilization of FVM discretization scheme, which protrudes by avoiding generation evident numerical dissipation (particularly commonly detected with the simultaneous specification of wave and current at inflow boundary) and disturbing the flow field [25, 26]. Within the specified source function region, the mass source function deduced based on the corresponding wave theory was added into the mass conservation equation (1):(4)Sm=2CAη.in whichC is the water wave phase velocity of the target wave. A is the area of the source region. η is the water wave surface elevation; for linear monochromatic wave, it can be represented as η=asinωt, and for solitary wave, it can be represented as the following formula:(5)ηx,t=Hsech23H4h3x−Ct.


## 2.3. Wave Absorption by Employing the Momentum Source Function Method

The analytical relaxation method proposed by Madsen et al. has been widely adopted by many numerical hydrodynamics researchers for establishing NWT (numerical wave tank) in pure-wave conditions. Here, the analytical relaxation method was applied for wave absorption in wave and current conditions and the elimination of wave reflection without varying the water depth and the current velocity at outflow boundary. The basic mechanism of the application of analytical relaxation method for wave absorption can be presented as following: Within the specified relaxation domain (usually with the length of one to two wavelengths), the water particle velocity and pressure will be updated at each calculational time step by the added momentum source function in the momentum equations to approach the target current flow velocity and the corresponding pressure. The relevant relaxation algorithm can be expressed as:(6)uM=CuJ+1−CU,vM=CvJ,pM=CpJ,where the subscript M represents the renewed physical quantity in the designated domain—the wave absorption domain, the subscript J represents the computed physical quantity at the previous calculational time step, U denotes the target current velocity, C=Cx represents the relaxation function, Cxmin=1, Cxmax=0, and more details about the analytical relaxation method and the relaxation function can be found in Madsen et al. [13].Then, the momentum source function adopted here for the analytical relaxation wave absorption can be reversely deduced from the Euler equations ignoring the water viscosity by substituting the aforementioned relaxation algorithm. The different discretization forms of the corresponding momentum equations with and without the additive momentum source can be expressed as the following:(7)ρuMn+1−uJnΔt+ρuJ∂uJ∂x+ρvJ∂uJ∂z=−∂pJ∂x+Sx,(8)ρvMn+1−vJnΔt+ρuJ∂vJ∂x+ρvJ∂vJ∂z=−∂pJ∂z+Sz−ρg,(9)ρuMn+1−uMnΔt+ρuM∂uM∂x+ρvM∂uM∂z=−∂pM∂x,(10)ρvMn+1−vMnΔt+ρuM∂vM∂x+ρvM∂vM∂z=−∂pM∂z−ρg.By subtracting equation (7) from (9) and equation (8) from (10), the corresponding momentum source function in abscissa and ordinate can be deduced as(11)Sx=1−C∂pJ∂x−ρuJ−UΔt+ρ1−C2uJ∂uJ∂x+vJ∂uJ∂z−ρ∗C∗1−C∗U∗uJ,(12)Sz=1−C∂pJ∂z−ρvJΔt+ρ1−C2uJ∂vJ∂x+vJ∂vJ∂z−ρ∗C∗1−C∗U∗vJ.


## 3. Numerical Methods and Validation

### 3.1. Numerical Methods

In order to implement the aforementioned numerical algorithms, the commercial Computational Fluid Dynamics (CFD) software-FLUENT 13.0, which has been corroborated to be viable for solving two-dimensional or three-dimensional wave hydrodynamics problems based on the N-S equations or RANS equations, is chosen as the base solver. For clarity and conciseness, only the principal numerical methods employed here are briefly described. The FLUENT solver employing the RANS equations as the governing equation is based on the FVM discretization scheme on the staggered grid system, which has been reported in research papers. For modelling turbulent flows, there are various turbulence models such as the one equation model, DNS model, two equations model, and the LES model available in the software. When deciding the most suitable turbulence model for implementing the aforementioned algorithms, the Reynolds Stress Model (RSM) presented better convergence property and higher accuracy than the other turbulence models. Therefore, the RSM utilizing a linear pressure-strain model is adopted in conjunction with the standard wall function method for handling near wall turbulence flow for all the following simulations. The volume-of-fluid method (VOF) algorithm coupled with a geometric reconstruction scheme is applied for free surface movement tracking. Fluid pressure and volume fraction function are arranged at the cell center of the computational grid system, and fluid velocity components are arranged at the face center of the corresponding cells. The body force-weighted scheme and the second-order upwind scheme are adopted for pressure interpolation and the discretization of the momentum equations, respectively, and the PISO algorithm is employed for pressure-velocity coupling when calculating convection and diffusion fluxes through the control-volume faces. The convergence criterion of all the variables is set to be 0.00001 for all cases.Figure1 shows the illustrative sketch of the computational domain representing the numerical wave and current tank (NWCT) established by utilizing the aforementioned internal mass source function wave-maker and the analytical relaxation wave absorber. The no-slip boundary in conjunction with the smooth wall function algorithm is set for bottom boundary, and the top boundary of the NWCT is set to be the pressure-outlet boundary. To generate a steady uniform current flow, the uniform current flow velocity is provided at the left-hand-side (LHS) and the right-hand-side (RHS) outflow boundary. After a stable uniform current flow is achieved in the whole computational domain, the target wave is generated within a designated region located in the middle of the abscissa of the NWCT by embedding the above-deduced mass source function Sm into the mass conservation equation (1) with the utilization of the internal macro of the UDF (User Defined Function)-DEFINE_SOURCE (mom_source, cell, thread, ds, eqn) defined in the software FLUENT. The generated wave in the specified region propagating to the left-hand-side (LHS) are opposing the pre-generated steady uniform current flow, while propagating to the right-hand-side (RHS) are following the pre-generated steady uniform current flows. At both ends of the abscissa with a length of approximately three wavelengths, the analytical relaxation wave absorber was employed to deplete the boundary reflection and actively assimilate the wave energy by embedding the above-deduced momentum source function presented by equations (11) and (12) into the momentum equations. When the wave and the current travel into the wave absorption regions, the water particle velocities are gradually and artificially updated by the analytical relaxation function to approach the current flow velocity, and the wave energy is absorbed. The calculational grid of the computational domain was generated with the utilization of the software GAMBIT by using the structured grids. The calculational grid step size and the temporal step size will affect the computational results to some extent, which will be comprehensively analyzed later.Figure 1
An illustrative sketch of the computational domain representing the NWCT.


### 3.2. Model Verification

To validate the aforementioned developed model, the laboratory measurements of Umeyama [1] are chosen for the comparison of the numerical results. Umeyama investigated the wave-current interactions by performing a series of laboratory experiments in a physical wave channel with a length of 25m, width of 0.7m, and maximum depth of 1.0m. Waves were generated with the utilization of a piston-type wave-maker placed at one end of the wave channel and were absorbed by a wave absorber installed at the other end. A pipe under the wave channel was utilized to recirculate the water flow pumped by a centrifugal pump to achieve the generation of a steady coplanar current with a depth-averaged velocity of U=0.08m/s. For all the investigations performed in the channel, the water depth (h) was kept at 0.3m and the target wave period was chosen as 1.0s; the experimental parameters for all the investigations are listed in Table 1. The free water surface elevation was measured with the utilization of a resistance-type probe, located 14.0m from the piston-type wave-maker, and the horizontal velocities distribution of water particles were measured by two nonintrusive measuring technique-PIV (particle image velocimetry) and PTV (particle tracking velocimetry), respectively. More details about the laboratory introduction can be found in Umeyama [1]. Considering the generated wave traveling with coplanar and counter current at the same time in the NWCT, a computational domain with a range of −20m to 20m in abscissa and0m to 0.5min ordinate is set in all the numerical simulations, and the designing of the source region differs from case to case, which will be further discussed in detail in Section 4.1. The calculational grid step size and the temporal grid size are determined by employing the sensitivity analysis. A grid refinement including spatial and temporal grid step study is carried out to ensure the accuracy of the NWCT established here. A linear wave with wavelength λ=1.37m, wave amplitude a=0.515cm , and wave period T=1sis generated in the NWCT at a constant water depth of h=0.3m. All numerical simulations are performed for a duration of 25 wave periods for ensuring the calculational stability. The temporal grid sensitivity analysis is performed first by utilizing the medium special grid system (containing 1000000 quadrilateral cells withΔx=0.01m, Δy=0.002m) and three different temporal grid sizes (0.01 s (T/100), 0.005 s (T/200), and 0.002 s (T/500)). The numerical results are presented in Figure 2(a). It is observed that the calculated time series of surface elevation at location x=10mfor all three temporal grid size cases conform well with each other revealing that the numerical model achieves convergence at the temporal grid. Consequently, for computational efficiency, the temporal grid step size is fixed at T/200 for all the following cases. Subsequently, the spatial grid sensitivity analysis is carried out based on three different special grid systems including the coarse (containing 250000 quadrilateral cells with Δx=0.02m, Δy=0.004m), medium (aforementioned grid system), and fine grid systems (containing 1440000 quadrilateral cells withΔx=0.005m, nonuniform in ordinate with Δy=0.005∼0.001m). The numerical results are shown in Figure 2(b). It is observed that the numerically calculated surface elevation converges better as the spatial grid is refined, but the improvement in the numerical results based on spatial grid refinement from the medium grid system to the fine grid system is unnoticeable. Considering the compromising of the accuracy of capturing the fluid interface and the computational efficiency, the medium grid system is utilized for all of the following simulations in this study.Table 1
Experimental cases.
Experimental parametersW1W2W3WC1WC2WC3Current velocity (cm/s)8.08.08.0Wave height (cm)1.032.343.610.912.023.09Wave period (s)1.01.01.01.01.01.0Wave steepness (ak)0.02360.05360.08270.02080.04630.0708Figure 2
Computational grid and time step convergence by comparing the time series of surface elevation atx = 10 m with three different calculational time steps and three different grid systems.
(a)(b)Figure3 shows the comparison of numerically calculated, physically measured, and analytically deduced phase-averaged surface displacement within one wave period in both pure-wave condition and wave-current interacting condition (following the current case). As expected, the simulated results agree well with the physical experimental results and the theoretical solution. It can be seen that the physical experimental results in case-WC1 depart from the theoretical solution and the numerical results (mainly generated during the physical laboratory testing and data processing procedures). By comparing the right figure (numerical results in wave-current interacting condition) with the left figure (numerical results in pure-wave condition) in the same row, it can be observed that the original wave pattern is varied by the current due to the wave-current interaction and as the wave is propagating in the uniform coplanar current, the wave height is decreased. The numerically calculated discrepancies of wave amplitude between W1-WC1, W2-WC2, and W3-WC3 are 0.06cm, 0.13cm and 0.21cm, respectively. For the waves with the same period propagating in the same following current, a larger wave height will result in a larger reduction of wave amplitude.Figure 3
Comparison of simulated and measured, analytical relative free surface elevation within one wave period under both pure-wave and wave-current interacting conditions. (a) W1. (b) WC1. (c) W2. (d) WC2. (e) W3. (f) WC3.
(a)(b)(c)(d)(e)(f)Figure4 displays the comparison of numerically and theoretically calculated free surface elevation distribution along the wave tank for wave-alone cases. Firstly, it can be observed that the wave elevation profiles in the working domain (three to five wavelengths away from the wave generation region) for all simulations fit well with the desired wave elevation profiles. In the two ends of the NWCT, namely, the wave absorption domains, the wave elevation approaches zero indicating that the analytical relaxation wave absorber can absorb the wave energy effectively to deplete the reflected waves generated from the left and right wall boundaries which will affect the wave profile in the working domain.Figure 4
Comparison between numerically simulated and theoretically calculated free surface elevation for wave-alone cases. (a) W1. (b) W2. (c) W3.
(a)(b)(c)Figure5 displays the comparison of simulated and theoretical free surface elevation distribution along the wave tank for wave-current cases (from x=−20 to x=0 for opposing current domain and from x=0 to x=20 for coplanar current domain). The theoretical free surface elevation is estimated by η=aω−k′U/ωcosk′x−ωt, in which k′can be solved by ω−k′U2=gk′tanhk′h. As can be found in Figure 5, the simulated results fit well with the theoretical results. Then, the effects of current on the water surface profile are investigated by choosing cases W3 & WC3 as examples. The averaged free surface elevations in the region from x=4L to 9L (for following current condition and pure-wave condition), x=−4L∼−9L (for opposing current condition) at t=20∼25T were used to compute the averaged wavelength and wave amplitude. By analyzing the effects of currents on the water surface profile as shown in Figure 6, it can be found that the wavelength is increased in the following current conditions and decreased in the opposing current conditions. This is due to the Doppler shift proposed by Wolf and Prandle [27]. Meanwhile, the wave amplitude gets smaller when waves travel with a following current and gets bigger with an opposing current. As a result, the wave steepness becomes larger due to shorter wavelength and increased wave amplitude in an opposing current and becomes smaller in a following current. The quantitative analysis about the influence of current flow velocity on the wave parameters will be discussed in Section 4.2.Figure 5
Comparison between numerically simulated and theoretically calculated free surface elevation for wave-current interacting cases. (a) WC1. (b) WC2. (c) WC3.
(a)(b)(c)Figure 6
Effects of currents on the water surface profile.Figures7–9 depict the comparison of numerically and theoretically calculated horizontal- and vertical- velocity profiles under pure-wave conditions (cases W1, W2, and W3), showing that they compare well with each other indicating the accuracy of the established numerical wave tank in calculating the flow velocity field during the generated wave propagating. Meanwhile, it can be inferred by comparing the flow velocity field under different wave conditions that the water wave particle velocity profiles are mainly affected by the surface wave motion under pure-wave conditions. As shown in the upper part of Figures 7–9, the absolute value of the vertical velocity of the water wave particles above the bottom boundary layer increases as moving toward the free surface (as depicted in the velocity profile at time t=0.25s and/or t=0.75s) when the zero wave elevation arrives. On the contrary, the absolute value of horizontal velocity of the water wave particles above the bottom boundary layer varies slightly as moving toward the free surface (as depicted in the velocity profile at time t=0.0s and/or t=0.5s) when the wave trough and/or wave crest arrive, as shown in the lower part of Figures 7–9. As a compensation step due to the lack of justifiable theoretical models for obtaining the velocity profiles under wave-current interacting conditions, the numerically calculated horizontal velocity profiles under waves interacting with coplanar current conditions was compared with the available experimental data measured with the utilization of PIV technology by professor Umeyama [1]. Figure 10 displays the comparison of numerically and theoretically calculated horizontal-velocity profiles under wave-current interacting conditions (cases WC1, WC2 and WC3 considering coplanar current conditions only). The numerically calculated horizontal-velocity profiles fit well with the relevant measured results and the relatively larger discrepancy between simulation and measurement is mainly caused by the scattered data of measurement (mainly generated during the physical laboratory testing and data processing procedures). By comparing the horizontal velocity profiles under pure-wave conditions depicted in Figures 7–9 with the corresponding results under wave-current interacting conditions displayed in Figure 10, it can be inferred that the wave-current interaction has obvious influence on water wave particle velocity profiles including the current flow and the surface wave motion. For the wave with small wave amplitude (WC1), the horizontal velocities varied slightly with the surface wave motion (the dominant influencing parameter is the current velocity). When the wave amplitude increased (WC3), the horizontal velocities in wave-current cases varied a lot with the surface wave motion.Figure 7
Comparison of numerically and theoretically calculated velocity profiles for case-W1. (a) Vertical velocity profiles. (b) Horizontal velocity profiles.
(a)(b)Figure 8
Comparison of numerically and theoretically calculated velocity profiles for case-W2. (a) Vertical velocity profiles. (b) Horizontal velocity profiles.
(a)(b)Figure 9
Comparison of numerically and theoretically calculated velocity profiles for case-W3. (a) Vertical velocity profiles. (b) Horizontal velocity profiles.
(a)(b)Figure 10
Comparison of numerically and theoretically calculated velocity profiles for wave-current interacting cases. (a) WC1. (b) WC2. (c) WC3.
(a)(b)(c)


## 3.1. Numerical Methods

In order to implement the aforementioned numerical algorithms, the commercial Computational Fluid Dynamics (CFD) software-FLUENT 13.0, which has been corroborated to be viable for solving two-dimensional or three-dimensional wave hydrodynamics problems based on the N-S equations or RANS equations, is chosen as the base solver. For clarity and conciseness, only the principal numerical methods employed here are briefly described. The FLUENT solver employing the RANS equations as the governing equation is based on the FVM discretization scheme on the staggered grid system, which has been reported in research papers. For modelling turbulent flows, there are various turbulence models such as the one equation model, DNS model, two equations model, and the LES model available in the software. When deciding the most suitable turbulence model for implementing the aforementioned algorithms, the Reynolds Stress Model (RSM) presented better convergence property and higher accuracy than the other turbulence models. Therefore, the RSM utilizing a linear pressure-strain model is adopted in conjunction with the standard wall function method for handling near wall turbulence flow for all the following simulations. The volume-of-fluid method (VOF) algorithm coupled with a geometric reconstruction scheme is applied for free surface movement tracking. Fluid pressure and volume fraction function are arranged at the cell center of the computational grid system, and fluid velocity components are arranged at the face center of the corresponding cells. The body force-weighted scheme and the second-order upwind scheme are adopted for pressure interpolation and the discretization of the momentum equations, respectively, and the PISO algorithm is employed for pressure-velocity coupling when calculating convection and diffusion fluxes through the control-volume faces. The convergence criterion of all the variables is set to be 0.00001 for all cases.Figure1 shows the illustrative sketch of the computational domain representing the numerical wave and current tank (NWCT) established by utilizing the aforementioned internal mass source function wave-maker and the analytical relaxation wave absorber. The no-slip boundary in conjunction with the smooth wall function algorithm is set for bottom boundary, and the top boundary of the NWCT is set to be the pressure-outlet boundary. To generate a steady uniform current flow, the uniform current flow velocity is provided at the left-hand-side (LHS) and the right-hand-side (RHS) outflow boundary. After a stable uniform current flow is achieved in the whole computational domain, the target wave is generated within a designated region located in the middle of the abscissa of the NWCT by embedding the above-deduced mass source function Sm into the mass conservation equation (1) with the utilization of the internal macro of the UDF (User Defined Function)-DEFINE_SOURCE (mom_source, cell, thread, ds, eqn) defined in the software FLUENT. The generated wave in the specified region propagating to the left-hand-side (LHS) are opposing the pre-generated steady uniform current flow, while propagating to the right-hand-side (RHS) are following the pre-generated steady uniform current flows. At both ends of the abscissa with a length of approximately three wavelengths, the analytical relaxation wave absorber was employed to deplete the boundary reflection and actively assimilate the wave energy by embedding the above-deduced momentum source function presented by equations (11) and (12) into the momentum equations. When the wave and the current travel into the wave absorption regions, the water particle velocities are gradually and artificially updated by the analytical relaxation function to approach the current flow velocity, and the wave energy is absorbed. The calculational grid of the computational domain was generated with the utilization of the software GAMBIT by using the structured grids. The calculational grid step size and the temporal step size will affect the computational results to some extent, which will be comprehensively analyzed later.Figure 1
An illustrative sketch of the computational domain representing the NWCT.


## 3.2. Model Verification

To validate the aforementioned developed model, the laboratory measurements of Umeyama [1] are chosen for the comparison of the numerical results. Umeyama investigated the wave-current interactions by performing a series of laboratory experiments in a physical wave channel with a length of 25m, width of 0.7m, and maximum depth of 1.0m. Waves were generated with the utilization of a piston-type wave-maker placed at one end of the wave channel and were absorbed by a wave absorber installed at the other end. A pipe under the wave channel was utilized to recirculate the water flow pumped by a centrifugal pump to achieve the generation of a steady coplanar current with a depth-averaged velocity of U=0.08m/s. For all the investigations performed in the channel, the water depth (h) was kept at 0.3m and the target wave period was chosen as 1.0s; the experimental parameters for all the investigations are listed in Table 1. The free water surface elevation was measured with the utilization of a resistance-type probe, located 14.0m from the piston-type wave-maker, and the horizontal velocities distribution of water particles were measured by two nonintrusive measuring technique-PIV (particle image velocimetry) and PTV (particle tracking velocimetry), respectively. More details about the laboratory introduction can be found in Umeyama [1]. Considering the generated wave traveling with coplanar and counter current at the same time in the NWCT, a computational domain with a range of −20m to 20m in abscissa and0m to 0.5min ordinate is set in all the numerical simulations, and the designing of the source region differs from case to case, which will be further discussed in detail in Section 4.1. The calculational grid step size and the temporal grid size are determined by employing the sensitivity analysis. A grid refinement including spatial and temporal grid step study is carried out to ensure the accuracy of the NWCT established here. A linear wave with wavelength λ=1.37m, wave amplitude a=0.515cm , and wave period T=1sis generated in the NWCT at a constant water depth of h=0.3m. All numerical simulations are performed for a duration of 25 wave periods for ensuring the calculational stability. The temporal grid sensitivity analysis is performed first by utilizing the medium special grid system (containing 1000000 quadrilateral cells withΔx=0.01m, Δy=0.002m) and three different temporal grid sizes (0.01 s (T/100), 0.005 s (T/200), and 0.002 s (T/500)). The numerical results are presented in Figure 2(a). It is observed that the calculated time series of surface elevation at location x=10mfor all three temporal grid size cases conform well with each other revealing that the numerical model achieves convergence at the temporal grid. Consequently, for computational efficiency, the temporal grid step size is fixed at T/200 for all the following cases. Subsequently, the spatial grid sensitivity analysis is carried out based on three different special grid systems including the coarse (containing 250000 quadrilateral cells with Δx=0.02m, Δy=0.004m), medium (aforementioned grid system), and fine grid systems (containing 1440000 quadrilateral cells withΔx=0.005m, nonuniform in ordinate with Δy=0.005∼0.001m). The numerical results are shown in Figure 2(b). It is observed that the numerically calculated surface elevation converges better as the spatial grid is refined, but the improvement in the numerical results based on spatial grid refinement from the medium grid system to the fine grid system is unnoticeable. Considering the compromising of the accuracy of capturing the fluid interface and the computational efficiency, the medium grid system is utilized for all of the following simulations in this study.Table 1
Experimental cases.
Experimental parametersW1W2W3WC1WC2WC3Current velocity (cm/s)8.08.08.0Wave height (cm)1.032.343.610.912.023.09Wave period (s)1.01.01.01.01.01.0Wave steepness (ak)0.02360.05360.08270.02080.04630.0708Figure 2
Computational grid and time step convergence by comparing the time series of surface elevation atx = 10 m with three different calculational time steps and three different grid systems.
(a)(b)Figure3 shows the comparison of numerically calculated, physically measured, and analytically deduced phase-averaged surface displacement within one wave period in both pure-wave condition and wave-current interacting condition (following the current case). As expected, the simulated results agree well with the physical experimental results and the theoretical solution. It can be seen that the physical experimental results in case-WC1 depart from the theoretical solution and the numerical results (mainly generated during the physical laboratory testing and data processing procedures). By comparing the right figure (numerical results in wave-current interacting condition) with the left figure (numerical results in pure-wave condition) in the same row, it can be observed that the original wave pattern is varied by the current due to the wave-current interaction and as the wave is propagating in the uniform coplanar current, the wave height is decreased. The numerically calculated discrepancies of wave amplitude between W1-WC1, W2-WC2, and W3-WC3 are 0.06cm, 0.13cm and 0.21cm, respectively. For the waves with the same period propagating in the same following current, a larger wave height will result in a larger reduction of wave amplitude.Figure 3
Comparison of simulated and measured, analytical relative free surface elevation within one wave period under both pure-wave and wave-current interacting conditions. (a) W1. (b) WC1. (c) W2. (d) WC2. (e) W3. (f) WC3.
(a)(b)(c)(d)(e)(f)Figure4 displays the comparison of numerically and theoretically calculated free surface elevation distribution along the wave tank for wave-alone cases. Firstly, it can be observed that the wave elevation profiles in the working domain (three to five wavelengths away from the wave generation region) for all simulations fit well with the desired wave elevation profiles. In the two ends of the NWCT, namely, the wave absorption domains, the wave elevation approaches zero indicating that the analytical relaxation wave absorber can absorb the wave energy effectively to deplete the reflected waves generated from the left and right wall boundaries which will affect the wave profile in the working domain.Figure 4
Comparison between numerically simulated and theoretically calculated free surface elevation for wave-alone cases. (a) W1. (b) W2. (c) W3.
(a)(b)(c)Figure5 displays the comparison of simulated and theoretical free surface elevation distribution along the wave tank for wave-current cases (from x=−20 to x=0 for opposing current domain and from x=0 to x=20 for coplanar current domain). The theoretical free surface elevation is estimated by η=aω−k′U/ωcosk′x−ωt, in which k′can be solved by ω−k′U2=gk′tanhk′h. As can be found in Figure 5, the simulated results fit well with the theoretical results. Then, the effects of current on the water surface profile are investigated by choosing cases W3 & WC3 as examples. The averaged free surface elevations in the region from x=4L to 9L (for following current condition and pure-wave condition), x=−4L∼−9L (for opposing current condition) at t=20∼25T were used to compute the averaged wavelength and wave amplitude. By analyzing the effects of currents on the water surface profile as shown in Figure 6, it can be found that the wavelength is increased in the following current conditions and decreased in the opposing current conditions. This is due to the Doppler shift proposed by Wolf and Prandle [27]. Meanwhile, the wave amplitude gets smaller when waves travel with a following current and gets bigger with an opposing current. As a result, the wave steepness becomes larger due to shorter wavelength and increased wave amplitude in an opposing current and becomes smaller in a following current. The quantitative analysis about the influence of current flow velocity on the wave parameters will be discussed in Section 4.2.Figure 5
Comparison between numerically simulated and theoretically calculated free surface elevation for wave-current interacting cases. (a) WC1. (b) WC2. (c) WC3.
(a)(b)(c)Figure 6
Effects of currents on the water surface profile.Figures7–9 depict the comparison of numerically and theoretically calculated horizontal- and vertical- velocity profiles under pure-wave conditions (cases W1, W2, and W3), showing that they compare well with each other indicating the accuracy of the established numerical wave tank in calculating the flow velocity field during the generated wave propagating. Meanwhile, it can be inferred by comparing the flow velocity field under different wave conditions that the water wave particle velocity profiles are mainly affected by the surface wave motion under pure-wave conditions. As shown in the upper part of Figures 7–9, the absolute value of the vertical velocity of the water wave particles above the bottom boundary layer increases as moving toward the free surface (as depicted in the velocity profile at time t=0.25s and/or t=0.75s) when the zero wave elevation arrives. On the contrary, the absolute value of horizontal velocity of the water wave particles above the bottom boundary layer varies slightly as moving toward the free surface (as depicted in the velocity profile at time t=0.0s and/or t=0.5s) when the wave trough and/or wave crest arrive, as shown in the lower part of Figures 7–9. As a compensation step due to the lack of justifiable theoretical models for obtaining the velocity profiles under wave-current interacting conditions, the numerically calculated horizontal velocity profiles under waves interacting with coplanar current conditions was compared with the available experimental data measured with the utilization of PIV technology by professor Umeyama [1]. Figure 10 displays the comparison of numerically and theoretically calculated horizontal-velocity profiles under wave-current interacting conditions (cases WC1, WC2 and WC3 considering coplanar current conditions only). The numerically calculated horizontal-velocity profiles fit well with the relevant measured results and the relatively larger discrepancy between simulation and measurement is mainly caused by the scattered data of measurement (mainly generated during the physical laboratory testing and data processing procedures). By comparing the horizontal velocity profiles under pure-wave conditions depicted in Figures 7–9 with the corresponding results under wave-current interacting conditions displayed in Figure 10, it can be inferred that the wave-current interaction has obvious influence on water wave particle velocity profiles including the current flow and the surface wave motion. For the wave with small wave amplitude (WC1), the horizontal velocities varied slightly with the surface wave motion (the dominant influencing parameter is the current velocity). When the wave amplitude increased (WC3), the horizontal velocities in wave-current cases varied a lot with the surface wave motion.Figure 7
Comparison of numerically and theoretically calculated velocity profiles for case-W1. (a) Vertical velocity profiles. (b) Horizontal velocity profiles.
(a)(b)Figure 8
Comparison of numerically and theoretically calculated velocity profiles for case-W2. (a) Vertical velocity profiles. (b) Horizontal velocity profiles.
(a)(b)Figure 9
Comparison of numerically and theoretically calculated velocity profiles for case-W3. (a) Vertical velocity profiles. (b) Horizontal velocity profiles.
(a)(b)Figure 10
Comparison of numerically and theoretically calculated velocity profiles for wave-current interacting cases. (a) WC1. (b) WC2. (c) WC3.
(a)(b)(c)


## 4. Model Discussion and Application

### 4.1. Region Design for the Mass Source Function Method

As mentioned by Chen et al., there are basically three principal parameters that will affect the accuracy of the wave generation with the utilization of the mass source function wave-maker, i.e., the width, height, and location of the mass source region, and should be specified by judicious study. Firstly, the influence of the location of the mass source region on wave generation was investigated by conducting a series of numerical simulations with the width of the source regionws fixed at ws=Δx (Δx-the spatial grid step in the abscissa), the height hs of the source region fixed aths=0.1h, and with the center of the region varied. The distance d between the center of the mass source region and the bottom boundary was employed to demonstrate the position of the mass source. The normalized wave amplitude (a′/a, where a′ is the numerically calculated wave amplitude and a is the desired wave amplitude) was used to represent the accuracy of the wave generation. As shown in Figure 11, the accuracy for all cases is acceptable. The closer the region moves to the bottom a′<a or to the free surface a′<a, the higher is the inaccuracy. For d/h<0.4, the inaccuracy of the wave-current cases is bigger than the wave-only cases. For d/h>0.65, the inaccuracy of the wave-current cases is less than the wave-only cases. In general, the center of the mass source region should be located at d/h∈0.3,0.65, i.e., the vertical distance d between the center of the mass source wave-maker and the bottom boundary of the NWT should be in the range of 0.3∼0.65 of the steady water depth, and the larger the wave amplitude the smaller the d/h.Figure 11
Simulated/measured amplitude relative ratio for different locations of the mass source.
(a)(b)(c)Then, the influence of the widthws of the source wave-maker on the wave generation will be discussed in detail. The W1-case was chosen as the example, and six different source widths (ws=Δx, ws=2Δx, ws=4Δx, ws=8Δx, ws=12Δx and ws=16Δx) were considered. The time series of surface elevation at x = 10 m was used to judge the quality of the generated waves. As shown in Figure 12, the width of the source region ws has little influence on wave generation. During the simulation process, it is noted that a smaller width of the source region ws will impose a larger mass source input at each time step and result in a smaller time step, and the smaller width is preferred for wave-current conditions. As shown in equation (4), the mass source input at each time step is related to the product of the width ws and the height hs of the mass source region, i.e., the area of the mass source region. Additional numerical tests show that that the generated wave is insensitive to the height hs of the source region as long as it is larger than the wave height, and the optimal height hs of the source region is about 10% of the water depth.Figure 12
Comparison of the generated waves with the utilization of different source widths.Overall, the designing of the mass source region will affect the accuracy and the convergence of the simulation and vary with the wave parameters (wave nonlinearity), current velocity, and water depth (relative water depth). For the present model, the quantitative analysis of the designing of the mass source region cannot be carried out because the opposing and coplanar current cases are simulated at the same time, and it will be accomplished by further work.


### 4.2. Influence of Current Flow Velocity on the Wave Parameters

More quantitative analysis was carried out to investigate the impact of current velocity on wave parameters. The input uniform current flow velocity was chosen as0m/s, ±0.08m/s, ±0.12m/s, ±0.16m/s, ±0.2m/s, and ±0.24m/s respectively, and the input wave parameters were chosen as: wave height H=0.0234m, wave period T1=1s,T2=1.5s,T3=2s, and the average water depth was fixed at h=0.3m. The averaged wave profiles generated in the working domain, x=2L to 5L (wave propagating with coplanar current) and x=−2L to −5L (wave propagating with opposing current) at t=20T to 25T were used to compute the averaged wavelength L′ and wave amplitude A′. The statistical wavelength L′, wave amplitude A′, and the corresponding change rate (ΔL/L and ΔA/A) are presented in Table 2, in which ΔL=L′−L, ΔA=A′−A (L and A represent the simulated wavelength and wave amplitude for pure-wave conditions). Figure 13 shows the rate of change of wave amplitude and wavelength against current velocity for fixed waves, and the horizontal axis is nondimensionalized by dividing the wave phase velocity and the vertical axis is nondimensionalized by dividing the wave parameters of the pure-wave conditions. It can be observed that the wavelength of the wave profile will be enlarged and the wave amplitude will be reduced when the generated wave is propagating with coplanar current and vice versa for opposing current condition. As the velocity of the current flow increases, the wave parameters of the wave profile vary more obviously, and the calculated rate of change is compared with the theoretical results (represented by the solid lines in Figure 13) deduced based on the potential flow theory proposed by Zou [28]. As shown in Figure 13, the numerically calculated change rate of the wave parameters of the wave profile in wave-current interacting conditions fit well with the relevant theoretical results, especially for comparison of the wavelength. The simulated change rate of the wave amplitude is smaller than the theoretical results for all cases mainly due to the potential flow theory ignoring the water viscosity, and the larger the current velocity the bigger the discrepancy.Table 2
Numerically calculated wave parameters under various wave-current interacting conditions.
U (m/s)Wave-T1Wave-T2Wave-T3L′ (m)A′ (m)L′ (m)A′ (m)L′ (m)A′ (m)ΔL/L (%)ΔA/A (%)ΔL/L (%)ΔA/A (%)ΔL/L (%)ΔA/A (%)−0.240.8190.01971.540.01772.180.0170−40.268.0−34.451.0−33.045.0−0.20.9260.01731.680.01612.380.0157−32.448.0−28.038.0−27.034.0−0.161.020.01581.830.01502.580.0149−25.235.0−22.028.0−21.027.0−0.121.170.01441.970.01382.770.0136−18.523.0−16.018.0−15.016.0−0.081.200.01332.080.01312.930.0129−12.114.0−11.012.0−10.010.001.370.01172.340.1173.260.1170.000.000.000.000.000.000.081.520.01042.570.01053.590.010511.0−11.110.0−10.010.0−10.00.121.600.009952.690.01023.720.010116.8−15.015.0−13.014.0−14.00.161.670.009362.810.009713.880.0097122.0−20.020.0−17.019.0−17.00.21.750.008892.900.009104.010.0093627.5−24.024.0−22.023.0−20.00.241.820.008543.000.008784.170.0086632.7−27.028.0−25.027.8−26.0Figure 13
Change of wave parameters against current speedL′A′ and LA represent wavelength (amplitude) with and without current effect, respectively; U represents the current velocity, C0 represents wave phase velocity for pure-wave condition, and the solid line ( ) represents the estimated results by the formula proposed by Zou.


### 4.3. Model Application

#### 4.3.1. Wave and Current Propagating over Submerged Breakwater

Water waves propagating over a submerged trapezoid breakwater are very useful and popular for showing the applicability of an established NWT in solving the wave-structure interaction problems [29]. Here, the famous experiment of Beji and Battjes is numerically simulated [30]. Firstly, a regular wave with wave height H=0.02m and wave period T=2.02s is generated at a constant still water depth h=0.4m by utilizing the model established above. Figure 14 depicts the illustrative sketch of the breakwater in the numerical wave tank, which is 49 m long (x=−21m to x=28m) and up to 0.5 m deep. The internal mass source wave-maker is located in the region with abscissa ranging from x=−6.6m to x=−6m and ordinate ranging fromy=0.1m to y=0.14m. The analytical wave absorbers with a length of about two wavelengths are placed at the two ends of the numerical wave tank for absorbing the wave energy and depleting the possible reflected waves caused by the wall boundary. The wave surface elevation around the breakwater was measured with the utilization of six wave gauges located at x=10.5m, x=12.5m, x=13.5m, x=14.5m, x=15.7m, and x=17.3m, respectively. The entire domain is discretized using nonuniform structured grid system, with Δx in the range of 0.002∼0.01m and Δy in the range of 0.0003∼0.002m (137 MB, containing 1470000 quadrilateral cells).Figure 14
Sketch of the established numerical wave tank for reproducing Beji’s experiment.Prior to simulating the wave-structure interaction, it is necessary to validate the applicability of the established NWT in simulating the wave generation and propagation over a flat bottom by removing the structure-breakwater. As shown in Figure15, the computed water surface elevation at the first wave gauge was observed to fit well with the desired water surface elevation, which confirms the accuracy of the NWT established hereunder specified wave conditions.Figure 15
Comparison of the numerically and theoretically calculated water surface elevation without submerged breakwater atx = 10.5 m.The formal simulations about wave propagating over the submerged trapezoid breakwater were then carried out, and the total computation time was 1 day 7 hours 20 minutes for a simulation time of35s. The desired wave was generated in the specified mass source region for propagating toward the breakwater. Figure 15 shows the comparison of the water surface elevation between the experimental data and the simulated results. By analyzing the wave surface elevation at different wave gauges, it can be found that the wave amplitude is larger than the target wave amplitude before reaching the top of the breakwater due to the effect of wave shoaling as shown in Figures 16(a) and 16(b). On the top water area of the submerged breakwater, a high-order harmonic wave starts to form and evolve. Meanwhile, a secondary wave appears as depicted in Figures 16(c) and 16(d). After propagating over the breakwater, the secondary wave mode obtains more energy from the main wave mode and the harmonic effect gets stronger as displayed in Figures 16(e) and 16(f). The prediction of wave transformation ranging in the region represented by wave gauges (e) and (f) is the toughest task due to the intricate flow separations and nonlinear wave energy transfers. It can be found that the numerically calculated wave surface elevations at six gauges with the utilization of the model developed here fit well with the experimental results, indicating that the model developed here has considerable applicability in studying the wave-structure interaction problems in the coastal/ocean engineering fields.Figure 16
Water surface elevation for regular wave propagation over submerged breakwater obtained from experimental data and simulated results. (a)x = 10.5 m. (b) x = 12.5 m. (c) x = 13.5 m. (d) x = 14.5 m. (e) x = 15.7 m. (f) x = 17.3 m.
(a)(b)(c)(d)(e)(f)In order to further present the transformation due to shoaling and breaking, the turbulent intensity distributions before and after wave breaking at four typical temporal grid nodes starting fromt0=32sare analyzed as depicted in Figure 17. Obviously, the turbulent intensity increased slightly as the nonlinearity of the generated waves became stronger and rapidly as wave breaking occurred. It can also be found that the turbulent intensity of the wave crest increased as the wave propagated toward the breakwater due to the increase of nonlinearity caused by wave shoaling (without wave breaking). However, the turbulent intensity numerically calculated in this paper with the utilization of the RSM turbulence model might be overestimated, especially at the crest of highly nonlinear waves [18].Figure 17
Contour map of numerically computed turbulent intensity at timet=t0, t=t0+T/4, t=t0+T/2, and t=t0+3T/4 (from top to bottom).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Then, the water wave interacting with coplanar current propagating over the breakwater was investigated. The current velocity was chosen asU=0.08m/s, Figure 18 depicts the comparison of the water surface elevation under wave-only condition and wave interacting with coplanar current condition. It is noticeable that the water surface elevation in these two conditions differs with each other obviously before reaching the top of the breakwater x=10.5m due to the wave scattering caused by the coexisting current, especially when on the top of the breakwater (x=12.5m, x=13.5m, x=14.5m). After propagating over the breakwater, the properties of the water surface elevation show similar characteristics. Overall, this numerical method based on the model developed here will play an important role in revealing the mechanism of the wave-current-breakwater interaction, which has not been explained so far.Figure 18
Water surface elevation for regular wave propagation with coplanar uniform current flow over a submerged breakwater. (a)x = 10.5 m. (b) x = 12.5 m. (c) x = 13.5 m. (d) x = 14.5 m. (e) x = 15.7 m. (f) x = 17.3 m.
(a)(b)(c)(d)(e)(f)


#### 4.3.2. Steep Water Wave Generation

The aforementioned analysis concerns only the tiny regular waves; whether the established NWCT is applicable or not in steeper wave conditions is still in doubt. To this end, we considered the simulation of steeper sea state waves in large scale. As the wave got steeper, the linear wave theory was not suitable anymore. In order to apply the presented model to generate steeper waves, the second-order Stokes wave theoretical solution was chosen as the progressive wave. Accordingly, the mass source function within the specified wave generation zone was added into the mass conservation equation (1), depicted by equation (13) as the following:(13)Sm=HCA∗−ak2sh2kh+coskx−wt+ak4chkh2ch2kh+1sh3khcos2kx−wt.A validation computation of a Stokes-wave train was carried out based on the following chosen parameters: the incoming wave height being set to beH=1m, the corresponding wave period being T=5.8s, the length of the numerical tank ranging from −40m to 200m, the total depth being 8m, still water depth chosen as h=6.2m, and the length of wave absorbing area equalling to 40 m at both ends. Figure 19 depicts the numerically computed time series of water surface elevation of a wave gauge located at x=50m and the instantaneous phase volume fraction distribution at time t=48s, respectively. From the left part, it can be observed that the numerically calculated time series of water wave surface elevation compared well with the analytical results. From the right part, the generated wave propagated well with almost nil wave attenuation and the momentum-source wave absorption model can absorb the wave energy effectively and deplete the possible reflected wave generated by the right wall boundary influencing the wave profile in the working zone.Figure 19
Computed time series of surface elevation atx=50m and the instantaneous phase volume fraction distribution at time t=48s.
(a)(b)Meanwhile, the pressure field of the generated wave was analyzed. By referring to the second-order Stokes wave theory, the pressure distribution in the wave field can be calculated as the following:(14)P=ρgηchkz+hchkh−ρgkA22sh2khsh2kz+h+ρg3kA221sh2khch2kz+hsh2kh−13cos2kx−ωt.Figure20 shows the comparison of numerically and theoretically calculated pressure distribution within one typical wavelength. It is obvious that the numerically calculated results show fairly good agreement with the theoretical results obtained from equation (14), indicating that the applicability of the established model for generating steep wave is extended.Figure 20
Comparisons of pressure contour between numerical results and theoretical results.


#### 4.3.3. Wave-Current Interactions under Deep Water Conditions

As stated in the first part, wave propagation on currents in deep water has been found in few studies. Here, the numerical wave and current tank (NWCT) was established according to the main scale of the multipurpose deep water tank of Dalian Maritime University. The abscissa of the present NWT with a total depth of11m ranges from −100m to 100m, and the still water depth is chosen and kept at h=10m. The target wavelength and the wave height are equal to 8mand 0.2m, respectively, and the current velocity is chosen to be U=0.1m/s. The internal mass source wave-maker is located in the region with the abscissa ranging from tox=−0.1mto x=0.1m , and the ordinate ranging from y−8.85m to y=9.35m The analytical wave absorbers with a length of 20 m are provided at both the ends of the tank for depleting the wave energy of the transmitted and reflected waves. The wave surface profiles were measured by four wave gauges in the experiment located at x=±32mx/L=±4, x=±48mx/L=±6, and the comparisons of monitored and analytical time series of water surface profile are displayed in Figure 21 (pure-wave case) and Figure 22 (wave-current case).Figure 21
Comparison of numerically and theoretically calculated time series of water wave surface profile in wave-alone case. (a)x/L=4, (b) x/L=−4, (c) x/L=6, and (d) x/L=−6.
(a)(b)(c)(d)Figure 22
Comparison of numerically and theoretically calculated water wave surface profiles in wave-current interacting cases. (a)x/L=4, (b) x/L=−4, (c) x/L=6, and (d) x/L=−6.
(a)(b)(c)(d)As shown in Figures21 and 22, the numerically and theoretically calculated time series of water wave surface profile in both pure-wave and wave-current interacting conditions compare well with each other, which proves that the numerical wave and current tank model established here can be developed and applied in solving wave-current interaction problems in deep water conditions. By comparing the calculated water surface profile in (a)–(d) in Figure 21, no wave attenuation was found among these four wave gauges located at x=±32mx/L=±4, x=±48mx/L=±6. Meanwhile, the wave propagates faster in coplanar-current direction x≥0 than in counter-current x≤0 direction as comparing (a) with (b) or (c) with (d) in Figure 22.Then, the variation of the instantaneous water wave surface elevation ranging in the whole computational domain was analyzed. As shown in Figure23, the numerically calculated water wave surface elevations under both wave-alone and wave-current interacting conditions fit well with the analytical water wave surface elevations. By analyzing the water wave free surface elevation distribution in the entire computational domain shown in the upper part of Figure 23 (wave-alone condition), the wave amplitudes generated are basically the same with very small attenuation, indicating that the model established here can generate accurate and stable repeated target waves. By analyzing the lower part of Figure 23 (wave-current condition), the water wave free surface elevation distribution in the entire computational domain can be divided into two parts: left part (counter-current direction) and right part (coplanar-current condition). In both counter-current direction and coplanar-current condition, the generated waveform repeated well and no attenuation was found, which means that the model established here can be employed to numerically solve the wave-current interaction problems in deep water conditions and that the wave-current interactions will not distort the waveform. However, by comparing the left part with the right part, as the wave propagated against the current (left part), the wave amplitude was increased and the wavelength was reduced. On the contrary, the wave amplitude was decreased and the wavelength was increased as the wave propagated with the current (right part).Figure 23
Comparison between numerically and theoretically calculated water wave free surface elevation for pure wave (a) and wave-current interacting cases (b).
(a)(b)


#### 4.3.4. Solitary Wave-Current Interaction

In some special marine environments, especially in the South China Sea, solitary waves occur frequently which may propagate toward coastal and near-shore regions and cause huge disasters. By interacting with riding tide current flows (coplanar currents), the wave speed of the propagating solitary wave may be vastly accelerated and the disaster responding time leaving for human being will be tremendously shorten although correspondingly the wave height may be decreased to a certain extent. On the other hand, the ebb tide current flow will enhance the influence of the solitary wave although with a lesser traveling speed. Here, the aforementioned numerical model proposed in this study is applied to investigate how the solitary wave interacts with the rising/ebb tide current flow. As an example, the water depth is chosen ash=8m, the wave height is chosen as H=1m, and the current velocity is chosen as ±1m/s. The numerical flume established here has an entire length of 800m and a total height of 10m. As for the pre-handling, the uniform mesh grid technology Δx=0.2m,Δy=0.004m is applied (containing one million quadrilateral cells). Figure 24 shows the effects of the current flow on the water surface profile at time, t=27s, by utilizing the comparisons of numerically calculated free surface under pure-wave and wave-current interacting (with coplanar and opposing currents) conditions. The numerically calculated wave heights, with the corresponding wave crest occurred at location x=−147.2m, x=−162m and162m for the case with an counter current U=−1m/s, no-current U=0m/s and following current U=1m/s are 1.25 m, 1m and 0.86 m respectively. It can be found that the change of the wave height (0.14 m, −14% of the height of the no-current solitary wave) induced by the following current is smaller than that (0.25 m, 25% of the height of the no-current solitary wave) induced by the opposing current. Meanwhile, the propagating velocity of the solitary wave is also altered by the ambient current flow. As for the opposing current case, the wave crest occurred x=−147.2m later than the current-free solitary wave case x=−162m. As for the following current case, the wave crest occurred x=178.6m earlier than the current-free solitary wave case x=162m.Figure 24
Effect of current velocity on the water surface profile of a solitary wave.


## 4.1. Region Design for the Mass Source Function Method

As mentioned by Chen et al., there are basically three principal parameters that will affect the accuracy of the wave generation with the utilization of the mass source function wave-maker, i.e., the width, height, and location of the mass source region, and should be specified by judicious study. Firstly, the influence of the location of the mass source region on wave generation was investigated by conducting a series of numerical simulations with the width of the source regionws fixed at ws=Δx (Δx-the spatial grid step in the abscissa), the height hs of the source region fixed aths=0.1h, and with the center of the region varied. The distance d between the center of the mass source region and the bottom boundary was employed to demonstrate the position of the mass source. The normalized wave amplitude (a′/a, where a′ is the numerically calculated wave amplitude and a is the desired wave amplitude) was used to represent the accuracy of the wave generation. As shown in Figure 11, the accuracy for all cases is acceptable. The closer the region moves to the bottom a′<a or to the free surface a′<a, the higher is the inaccuracy. For d/h<0.4, the inaccuracy of the wave-current cases is bigger than the wave-only cases. For d/h>0.65, the inaccuracy of the wave-current cases is less than the wave-only cases. In general, the center of the mass source region should be located at d/h∈0.3,0.65, i.e., the vertical distance d between the center of the mass source wave-maker and the bottom boundary of the NWT should be in the range of 0.3∼0.65 of the steady water depth, and the larger the wave amplitude the smaller the d/h.Figure 11
Simulated/measured amplitude relative ratio for different locations of the mass source.
(a)(b)(c)Then, the influence of the widthws of the source wave-maker on the wave generation will be discussed in detail. The W1-case was chosen as the example, and six different source widths (ws=Δx, ws=2Δx, ws=4Δx, ws=8Δx, ws=12Δx and ws=16Δx) were considered. The time series of surface elevation at x = 10 m was used to judge the quality of the generated waves. As shown in Figure 12, the width of the source region ws has little influence on wave generation. During the simulation process, it is noted that a smaller width of the source region ws will impose a larger mass source input at each time step and result in a smaller time step, and the smaller width is preferred for wave-current conditions. As shown in equation (4), the mass source input at each time step is related to the product of the width ws and the height hs of the mass source region, i.e., the area of the mass source region. Additional numerical tests show that that the generated wave is insensitive to the height hs of the source region as long as it is larger than the wave height, and the optimal height hs of the source region is about 10% of the water depth.Figure 12
Comparison of the generated waves with the utilization of different source widths.Overall, the designing of the mass source region will affect the accuracy and the convergence of the simulation and vary with the wave parameters (wave nonlinearity), current velocity, and water depth (relative water depth). For the present model, the quantitative analysis of the designing of the mass source region cannot be carried out because the opposing and coplanar current cases are simulated at the same time, and it will be accomplished by further work.


## 4.2. Influence of Current Flow Velocity on the Wave Parameters

More quantitative analysis was carried out to investigate the impact of current velocity on wave parameters. The input uniform current flow velocity was chosen as0m/s, ±0.08m/s, ±0.12m/s, ±0.16m/s, ±0.2m/s, and ±0.24m/s respectively, and the input wave parameters were chosen as: wave height H=0.0234m, wave period T1=1s,T2=1.5s,T3=2s, and the average water depth was fixed at h=0.3m. The averaged wave profiles generated in the working domain, x=2L to 5L (wave propagating with coplanar current) and x=−2L to −5L (wave propagating with opposing current) at t=20T to 25T were used to compute the averaged wavelength L′ and wave amplitude A′. The statistical wavelength L′, wave amplitude A′, and the corresponding change rate (ΔL/L and ΔA/A) are presented in Table 2, in which ΔL=L′−L, ΔA=A′−A (L and A represent the simulated wavelength and wave amplitude for pure-wave conditions). Figure 13 shows the rate of change of wave amplitude and wavelength against current velocity for fixed waves, and the horizontal axis is nondimensionalized by dividing the wave phase velocity and the vertical axis is nondimensionalized by dividing the wave parameters of the pure-wave conditions. It can be observed that the wavelength of the wave profile will be enlarged and the wave amplitude will be reduced when the generated wave is propagating with coplanar current and vice versa for opposing current condition. As the velocity of the current flow increases, the wave parameters of the wave profile vary more obviously, and the calculated rate of change is compared with the theoretical results (represented by the solid lines in Figure 13) deduced based on the potential flow theory proposed by Zou [28]. As shown in Figure 13, the numerically calculated change rate of the wave parameters of the wave profile in wave-current interacting conditions fit well with the relevant theoretical results, especially for comparison of the wavelength. The simulated change rate of the wave amplitude is smaller than the theoretical results for all cases mainly due to the potential flow theory ignoring the water viscosity, and the larger the current velocity the bigger the discrepancy.Table 2
Numerically calculated wave parameters under various wave-current interacting conditions.
U (m/s)Wave-T1Wave-T2Wave-T3L′ (m)A′ (m)L′ (m)A′ (m)L′ (m)A′ (m)ΔL/L (%)ΔA/A (%)ΔL/L (%)ΔA/A (%)ΔL/L (%)ΔA/A (%)−0.240.8190.01971.540.01772.180.0170−40.268.0−34.451.0−33.045.0−0.20.9260.01731.680.01612.380.0157−32.448.0−28.038.0−27.034.0−0.161.020.01581.830.01502.580.0149−25.235.0−22.028.0−21.027.0−0.121.170.01441.970.01382.770.0136−18.523.0−16.018.0−15.016.0−0.081.200.01332.080.01312.930.0129−12.114.0−11.012.0−10.010.001.370.01172.340.1173.260.1170.000.000.000.000.000.000.081.520.01042.570.01053.590.010511.0−11.110.0−10.010.0−10.00.121.600.009952.690.01023.720.010116.8−15.015.0−13.014.0−14.00.161.670.009362.810.009713.880.0097122.0−20.020.0−17.019.0−17.00.21.750.008892.900.009104.010.0093627.5−24.024.0−22.023.0−20.00.241.820.008543.000.008784.170.0086632.7−27.028.0−25.027.8−26.0Figure 13
Change of wave parameters against current speedL′A′ and LA represent wavelength (amplitude) with and without current effect, respectively; U represents the current velocity, C0 represents wave phase velocity for pure-wave condition, and the solid line ( ) represents the estimated results by the formula proposed by Zou.


## 4.3. Model Application

### 4.3.1. Wave and Current Propagating over Submerged Breakwater

Water waves propagating over a submerged trapezoid breakwater are very useful and popular for showing the applicability of an established NWT in solving the wave-structure interaction problems [29]. Here, the famous experiment of Beji and Battjes is numerically simulated [30]. Firstly, a regular wave with wave height H=0.02m and wave period T=2.02s is generated at a constant still water depth h=0.4m by utilizing the model established above. Figure 14 depicts the illustrative sketch of the breakwater in the numerical wave tank, which is 49 m long (x=−21m to x=28m) and up to 0.5 m deep. The internal mass source wave-maker is located in the region with abscissa ranging from x=−6.6m to x=−6m and ordinate ranging fromy=0.1m to y=0.14m. The analytical wave absorbers with a length of about two wavelengths are placed at the two ends of the numerical wave tank for absorbing the wave energy and depleting the possible reflected waves caused by the wall boundary. The wave surface elevation around the breakwater was measured with the utilization of six wave gauges located at x=10.5m, x=12.5m, x=13.5m, x=14.5m, x=15.7m, and x=17.3m, respectively. The entire domain is discretized using nonuniform structured grid system, with Δx in the range of 0.002∼0.01m and Δy in the range of 0.0003∼0.002m (137 MB, containing 1470000 quadrilateral cells).Figure 14
Sketch of the established numerical wave tank for reproducing Beji’s experiment.Prior to simulating the wave-structure interaction, it is necessary to validate the applicability of the established NWT in simulating the wave generation and propagation over a flat bottom by removing the structure-breakwater. As shown in Figure15, the computed water surface elevation at the first wave gauge was observed to fit well with the desired water surface elevation, which confirms the accuracy of the NWT established hereunder specified wave conditions.Figure 15
Comparison of the numerically and theoretically calculated water surface elevation without submerged breakwater atx = 10.5 m.The formal simulations about wave propagating over the submerged trapezoid breakwater were then carried out, and the total computation time was 1 day 7 hours 20 minutes for a simulation time of35s. The desired wave was generated in the specified mass source region for propagating toward the breakwater. Figure 15 shows the comparison of the water surface elevation between the experimental data and the simulated results. By analyzing the wave surface elevation at different wave gauges, it can be found that the wave amplitude is larger than the target wave amplitude before reaching the top of the breakwater due to the effect of wave shoaling as shown in Figures 16(a) and 16(b). On the top water area of the submerged breakwater, a high-order harmonic wave starts to form and evolve. Meanwhile, a secondary wave appears as depicted in Figures 16(c) and 16(d). After propagating over the breakwater, the secondary wave mode obtains more energy from the main wave mode and the harmonic effect gets stronger as displayed in Figures 16(e) and 16(f). The prediction of wave transformation ranging in the region represented by wave gauges (e) and (f) is the toughest task due to the intricate flow separations and nonlinear wave energy transfers. It can be found that the numerically calculated wave surface elevations at six gauges with the utilization of the model developed here fit well with the experimental results, indicating that the model developed here has considerable applicability in studying the wave-structure interaction problems in the coastal/ocean engineering fields.Figure 16
Water surface elevation for regular wave propagation over submerged breakwater obtained from experimental data and simulated results. (a)x = 10.5 m. (b) x = 12.5 m. (c) x = 13.5 m. (d) x = 14.5 m. (e) x = 15.7 m. (f) x = 17.3 m.
(a)(b)(c)(d)(e)(f)In order to further present the transformation due to shoaling and breaking, the turbulent intensity distributions before and after wave breaking at four typical temporal grid nodes starting fromt0=32sare analyzed as depicted in Figure 17. Obviously, the turbulent intensity increased slightly as the nonlinearity of the generated waves became stronger and rapidly as wave breaking occurred. It can also be found that the turbulent intensity of the wave crest increased as the wave propagated toward the breakwater due to the increase of nonlinearity caused by wave shoaling (without wave breaking). However, the turbulent intensity numerically calculated in this paper with the utilization of the RSM turbulence model might be overestimated, especially at the crest of highly nonlinear waves [18].Figure 17
Contour map of numerically computed turbulent intensity at timet=t0, t=t0+T/4, t=t0+T/2, and t=t0+3T/4 (from top to bottom).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Then, the water wave interacting with coplanar current propagating over the breakwater was investigated. The current velocity was chosen asU=0.08m/s, Figure 18 depicts the comparison of the water surface elevation under wave-only condition and wave interacting with coplanar current condition. It is noticeable that the water surface elevation in these two conditions differs with each other obviously before reaching the top of the breakwater x=10.5m due to the wave scattering caused by the coexisting current, especially when on the top of the breakwater (x=12.5m, x=13.5m, x=14.5m). After propagating over the breakwater, the properties of the water surface elevation show similar characteristics. Overall, this numerical method based on the model developed here will play an important role in revealing the mechanism of the wave-current-breakwater interaction, which has not been explained so far.Figure 18
Water surface elevation for regular wave propagation with coplanar uniform current flow over a submerged breakwater. (a)x = 10.5 m. (b) x = 12.5 m. (c) x = 13.5 m. (d) x = 14.5 m. (e) x = 15.7 m. (f) x = 17.3 m.
(a)(b)(c)(d)(e)(f)


### 4.3.2. Steep Water Wave Generation

The aforementioned analysis concerns only the tiny regular waves; whether the established NWCT is applicable or not in steeper wave conditions is still in doubt. To this end, we considered the simulation of steeper sea state waves in large scale. As the wave got steeper, the linear wave theory was not suitable anymore. In order to apply the presented model to generate steeper waves, the second-order Stokes wave theoretical solution was chosen as the progressive wave. Accordingly, the mass source function within the specified wave generation zone was added into the mass conservation equation (1), depicted by equation (13) as the following:(13)Sm=HCA∗−ak2sh2kh+coskx−wt+ak4chkh2ch2kh+1sh3khcos2kx−wt.A validation computation of a Stokes-wave train was carried out based on the following chosen parameters: the incoming wave height being set to beH=1m, the corresponding wave period being T=5.8s, the length of the numerical tank ranging from −40m to 200m, the total depth being 8m, still water depth chosen as h=6.2m, and the length of wave absorbing area equalling to 40 m at both ends. Figure 19 depicts the numerically computed time series of water surface elevation of a wave gauge located at x=50m and the instantaneous phase volume fraction distribution at time t=48s, respectively. From the left part, it can be observed that the numerically calculated time series of water wave surface elevation compared well with the analytical results. From the right part, the generated wave propagated well with almost nil wave attenuation and the momentum-source wave absorption model can absorb the wave energy effectively and deplete the possible reflected wave generated by the right wall boundary influencing the wave profile in the working zone.Figure 19
Computed time series of surface elevation atx=50m and the instantaneous phase volume fraction distribution at time t=48s.
(a)(b)Meanwhile, the pressure field of the generated wave was analyzed. By referring to the second-order Stokes wave theory, the pressure distribution in the wave field can be calculated as the following:(14)P=ρgηchkz+hchkh−ρgkA22sh2khsh2kz+h+ρg3kA221sh2khch2kz+hsh2kh−13cos2kx−ωt.Figure20 shows the comparison of numerically and theoretically calculated pressure distribution within one typical wavelength. It is obvious that the numerically calculated results show fairly good agreement with the theoretical results obtained from equation (14), indicating that the applicability of the established model for generating steep wave is extended.Figure 20
Comparisons of pressure contour between numerical results and theoretical results.


### 4.3.3. Wave-Current Interactions under Deep Water Conditions

As stated in the first part, wave propagation on currents in deep water has been found in few studies. Here, the numerical wave and current tank (NWCT) was established according to the main scale of the multipurpose deep water tank of Dalian Maritime University. The abscissa of the present NWT with a total depth of11m ranges from −100m to 100m, and the still water depth is chosen and kept at h=10m. The target wavelength and the wave height are equal to 8mand 0.2m, respectively, and the current velocity is chosen to be U=0.1m/s. The internal mass source wave-maker is located in the region with the abscissa ranging from tox=−0.1mto x=0.1m , and the ordinate ranging from y−8.85m to y=9.35m The analytical wave absorbers with a length of 20 m are provided at both the ends of the tank for depleting the wave energy of the transmitted and reflected waves. The wave surface profiles were measured by four wave gauges in the experiment located at x=±32mx/L=±4, x=±48mx/L=±6, and the comparisons of monitored and analytical time series of water surface profile are displayed in Figure 21 (pure-wave case) and Figure 22 (wave-current case).Figure 21
Comparison of numerically and theoretically calculated time series of water wave surface profile in wave-alone case. (a)x/L=4, (b) x/L=−4, (c) x/L=6, and (d) x/L=−6.
(a)(b)(c)(d)Figure 22
Comparison of numerically and theoretically calculated water wave surface profiles in wave-current interacting cases. (a)x/L=4, (b) x/L=−4, (c) x/L=6, and (d) x/L=−6.
(a)(b)(c)(d)As shown in Figures21 and 22, the numerically and theoretically calculated time series of water wave surface profile in both pure-wave and wave-current interacting conditions compare well with each other, which proves that the numerical wave and current tank model established here can be developed and applied in solving wave-current interaction problems in deep water conditions. By comparing the calculated water surface profile in (a)–(d) in Figure 21, no wave attenuation was found among these four wave gauges located at x=±32mx/L=±4, x=±48mx/L=±6. Meanwhile, the wave propagates faster in coplanar-current direction x≥0 than in counter-current x≤0 direction as comparing (a) with (b) or (c) with (d) in Figure 22.Then, the variation of the instantaneous water wave surface elevation ranging in the whole computational domain was analyzed. As shown in Figure23, the numerically calculated water wave surface elevations under both wave-alone and wave-current interacting conditions fit well with the analytical water wave surface elevations. By analyzing the water wave free surface elevation distribution in the entire computational domain shown in the upper part of Figure 23 (wave-alone condition), the wave amplitudes generated are basically the same with very small attenuation, indicating that the model established here can generate accurate and stable repeated target waves. By analyzing the lower part of Figure 23 (wave-current condition), the water wave free surface elevation distribution in the entire computational domain can be divided into two parts: left part (counter-current direction) and right part (coplanar-current condition). In both counter-current direction and coplanar-current condition, the generated waveform repeated well and no attenuation was found, which means that the model established here can be employed to numerically solve the wave-current interaction problems in deep water conditions and that the wave-current interactions will not distort the waveform. However, by comparing the left part with the right part, as the wave propagated against the current (left part), the wave amplitude was increased and the wavelength was reduced. On the contrary, the wave amplitude was decreased and the wavelength was increased as the wave propagated with the current (right part).Figure 23
Comparison between numerically and theoretically calculated water wave free surface elevation for pure wave (a) and wave-current interacting cases (b).
(a)(b)


### 4.3.4. Solitary Wave-Current Interaction

In some special marine environments, especially in the South China Sea, solitary waves occur frequently which may propagate toward coastal and near-shore regions and cause huge disasters. By interacting with riding tide current flows (coplanar currents), the wave speed of the propagating solitary wave may be vastly accelerated and the disaster responding time leaving for human being will be tremendously shorten although correspondingly the wave height may be decreased to a certain extent. On the other hand, the ebb tide current flow will enhance the influence of the solitary wave although with a lesser traveling speed. Here, the aforementioned numerical model proposed in this study is applied to investigate how the solitary wave interacts with the rising/ebb tide current flow. As an example, the water depth is chosen ash=8m, the wave height is chosen as H=1m, and the current velocity is chosen as ±1m/s. The numerical flume established here has an entire length of 800m and a total height of 10m. As for the pre-handling, the uniform mesh grid technology Δx=0.2m,Δy=0.004m is applied (containing one million quadrilateral cells). Figure 24 shows the effects of the current flow on the water surface profile at time, t=27s, by utilizing the comparisons of numerically calculated free surface under pure-wave and wave-current interacting (with coplanar and opposing currents) conditions. The numerically calculated wave heights, with the corresponding wave crest occurred at location x=−147.2m, x=−162m and162m for the case with an counter current U=−1m/s, no-current U=0m/s and following current U=1m/s are 1.25 m, 1m and 0.86 m respectively. It can be found that the change of the wave height (0.14 m, −14% of the height of the no-current solitary wave) induced by the following current is smaller than that (0.25 m, 25% of the height of the no-current solitary wave) induced by the opposing current. Meanwhile, the propagating velocity of the solitary wave is also altered by the ambient current flow. As for the opposing current case, the wave crest occurred x=−147.2m later than the current-free solitary wave case x=−162m. As for the following current case, the wave crest occurred x=178.6m earlier than the current-free solitary wave case x=162m.Figure 24
Effect of current velocity on the water surface profile of a solitary wave.


## 4.3.1. Wave and Current Propagating over Submerged Breakwater

Water waves propagating over a submerged trapezoid breakwater are very useful and popular for showing the applicability of an established NWT in solving the wave-structure interaction problems [29]. Here, the famous experiment of Beji and Battjes is numerically simulated [30]. Firstly, a regular wave with wave height H=0.02m and wave period T=2.02s is generated at a constant still water depth h=0.4m by utilizing the model established above. Figure 14 depicts the illustrative sketch of the breakwater in the numerical wave tank, which is 49 m long (x=−21m to x=28m) and up to 0.5 m deep. The internal mass source wave-maker is located in the region with abscissa ranging from x=−6.6m to x=−6m and ordinate ranging fromy=0.1m to y=0.14m. The analytical wave absorbers with a length of about two wavelengths are placed at the two ends of the numerical wave tank for absorbing the wave energy and depleting the possible reflected waves caused by the wall boundary. The wave surface elevation around the breakwater was measured with the utilization of six wave gauges located at x=10.5m, x=12.5m, x=13.5m, x=14.5m, x=15.7m, and x=17.3m, respectively. The entire domain is discretized using nonuniform structured grid system, with Δx in the range of 0.002∼0.01m and Δy in the range of 0.0003∼0.002m (137 MB, containing 1470000 quadrilateral cells).Figure 14
Sketch of the established numerical wave tank for reproducing Beji’s experiment.Prior to simulating the wave-structure interaction, it is necessary to validate the applicability of the established NWT in simulating the wave generation and propagation over a flat bottom by removing the structure-breakwater. As shown in Figure15, the computed water surface elevation at the first wave gauge was observed to fit well with the desired water surface elevation, which confirms the accuracy of the NWT established hereunder specified wave conditions.Figure 15
Comparison of the numerically and theoretically calculated water surface elevation without submerged breakwater atx = 10.5 m.The formal simulations about wave propagating over the submerged trapezoid breakwater were then carried out, and the total computation time was 1 day 7 hours 20 minutes for a simulation time of35s. The desired wave was generated in the specified mass source region for propagating toward the breakwater. Figure 15 shows the comparison of the water surface elevation between the experimental data and the simulated results. By analyzing the wave surface elevation at different wave gauges, it can be found that the wave amplitude is larger than the target wave amplitude before reaching the top of the breakwater due to the effect of wave shoaling as shown in Figures 16(a) and 16(b). On the top water area of the submerged breakwater, a high-order harmonic wave starts to form and evolve. Meanwhile, a secondary wave appears as depicted in Figures 16(c) and 16(d). After propagating over the breakwater, the secondary wave mode obtains more energy from the main wave mode and the harmonic effect gets stronger as displayed in Figures 16(e) and 16(f). The prediction of wave transformation ranging in the region represented by wave gauges (e) and (f) is the toughest task due to the intricate flow separations and nonlinear wave energy transfers. It can be found that the numerically calculated wave surface elevations at six gauges with the utilization of the model developed here fit well with the experimental results, indicating that the model developed here has considerable applicability in studying the wave-structure interaction problems in the coastal/ocean engineering fields.Figure 16
Water surface elevation for regular wave propagation over submerged breakwater obtained from experimental data and simulated results. (a)x = 10.5 m. (b) x = 12.5 m. (c) x = 13.5 m. (d) x = 14.5 m. (e) x = 15.7 m. (f) x = 17.3 m.
(a)(b)(c)(d)(e)(f)In order to further present the transformation due to shoaling and breaking, the turbulent intensity distributions before and after wave breaking at four typical temporal grid nodes starting fromt0=32sare analyzed as depicted in Figure 17. Obviously, the turbulent intensity increased slightly as the nonlinearity of the generated waves became stronger and rapidly as wave breaking occurred. It can also be found that the turbulent intensity of the wave crest increased as the wave propagated toward the breakwater due to the increase of nonlinearity caused by wave shoaling (without wave breaking). However, the turbulent intensity numerically calculated in this paper with the utilization of the RSM turbulence model might be overestimated, especially at the crest of highly nonlinear waves [18].Figure 17
Contour map of numerically computed turbulent intensity at timet=t0, t=t0+T/4, t=t0+T/2, and t=t0+3T/4 (from top to bottom).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Then, the water wave interacting with coplanar current propagating over the breakwater was investigated. The current velocity was chosen asU=0.08m/s, Figure 18 depicts the comparison of the water surface elevation under wave-only condition and wave interacting with coplanar current condition. It is noticeable that the water surface elevation in these two conditions differs with each other obviously before reaching the top of the breakwater x=10.5m due to the wave scattering caused by the coexisting current, especially when on the top of the breakwater (x=12.5m, x=13.5m, x=14.5m). After propagating over the breakwater, the properties of the water surface elevation show similar characteristics. Overall, this numerical method based on the model developed here will play an important role in revealing the mechanism of the wave-current-breakwater interaction, which has not been explained so far.Figure 18
Water surface elevation for regular wave propagation with coplanar uniform current flow over a submerged breakwater. (a)x = 10.5 m. (b) x = 12.5 m. (c) x = 13.5 m. (d) x = 14.5 m. (e) x = 15.7 m. (f) x = 17.3 m.
(a)(b)(c)(d)(e)(f)


## 4.3.2. Steep Water Wave Generation

The aforementioned analysis concerns only the tiny regular waves; whether the established NWCT is applicable or not in steeper wave conditions is still in doubt. To this end, we considered the simulation of steeper sea state waves in large scale. As the wave got steeper, the linear wave theory was not suitable anymore. In order to apply the presented model to generate steeper waves, the second-order Stokes wave theoretical solution was chosen as the progressive wave. Accordingly, the mass source function within the specified wave generation zone was added into the mass conservation equation (1), depicted by equation (13) as the following:(13)Sm=HCA∗−ak2sh2kh+coskx−wt+ak4chkh2ch2kh+1sh3khcos2kx−wt.A validation computation of a Stokes-wave train was carried out based on the following chosen parameters: the incoming wave height being set to beH=1m, the corresponding wave period being T=5.8s, the length of the numerical tank ranging from −40m to 200m, the total depth being 8m, still water depth chosen as h=6.2m, and the length of wave absorbing area equalling to 40 m at both ends. Figure 19 depicts the numerically computed time series of water surface elevation of a wave gauge located at x=50m and the instantaneous phase volume fraction distribution at time t=48s, respectively. From the left part, it can be observed that the numerically calculated time series of water wave surface elevation compared well with the analytical results. From the right part, the generated wave propagated well with almost nil wave attenuation and the momentum-source wave absorption model can absorb the wave energy effectively and deplete the possible reflected wave generated by the right wall boundary influencing the wave profile in the working zone.Figure 19
Computed time series of surface elevation atx=50m and the instantaneous phase volume fraction distribution at time t=48s.
(a)(b)Meanwhile, the pressure field of the generated wave was analyzed. By referring to the second-order Stokes wave theory, the pressure distribution in the wave field can be calculated as the following:(14)P=ρgηchkz+hchkh−ρgkA22sh2khsh2kz+h+ρg3kA221sh2khch2kz+hsh2kh−13cos2kx−ωt.Figure20 shows the comparison of numerically and theoretically calculated pressure distribution within one typical wavelength. It is obvious that the numerically calculated results show fairly good agreement with the theoretical results obtained from equation (14), indicating that the applicability of the established model for generating steep wave is extended.Figure 20
Comparisons of pressure contour between numerical results and theoretical results.


## 4.3.3. Wave-Current Interactions under Deep Water Conditions

As stated in the first part, wave propagation on currents in deep water has been found in few studies. Here, the numerical wave and current tank (NWCT) was established according to the main scale of the multipurpose deep water tank of Dalian Maritime University. The abscissa of the present NWT with a total depth of11m ranges from −100m to 100m, and the still water depth is chosen and kept at h=10m. The target wavelength and the wave height are equal to 8mand 0.2m, respectively, and the current velocity is chosen to be U=0.1m/s. The internal mass source wave-maker is located in the region with the abscissa ranging from tox=−0.1mto x=0.1m , and the ordinate ranging from y−8.85m to y=9.35m The analytical wave absorbers with a length of 20 m are provided at both the ends of the tank for depleting the wave energy of the transmitted and reflected waves. The wave surface profiles were measured by four wave gauges in the experiment located at x=±32mx/L=±4, x=±48mx/L=±6, and the comparisons of monitored and analytical time series of water surface profile are displayed in Figure 21 (pure-wave case) and Figure 22 (wave-current case).Figure 21
Comparison of numerically and theoretically calculated time series of water wave surface profile in wave-alone case. (a)x/L=4, (b) x/L=−4, (c) x/L=6, and (d) x/L=−6.
(a)(b)(c)(d)Figure 22
Comparison of numerically and theoretically calculated water wave surface profiles in wave-current interacting cases. (a)x/L=4, (b) x/L=−4, (c) x/L=6, and (d) x/L=−6.
(a)(b)(c)(d)As shown in Figures21 and 22, the numerically and theoretically calculated time series of water wave surface profile in both pure-wave and wave-current interacting conditions compare well with each other, which proves that the numerical wave and current tank model established here can be developed and applied in solving wave-current interaction problems in deep water conditions. By comparing the calculated water surface profile in (a)–(d) in Figure 21, no wave attenuation was found among these four wave gauges located at x=±32mx/L=±4, x=±48mx/L=±6. Meanwhile, the wave propagates faster in coplanar-current direction x≥0 than in counter-current x≤0 direction as comparing (a) with (b) or (c) with (d) in Figure 22.Then, the variation of the instantaneous water wave surface elevation ranging in the whole computational domain was analyzed. As shown in Figure23, the numerically calculated water wave surface elevations under both wave-alone and wave-current interacting conditions fit well with the analytical water wave surface elevations. By analyzing the water wave free surface elevation distribution in the entire computational domain shown in the upper part of Figure 23 (wave-alone condition), the wave amplitudes generated are basically the same with very small attenuation, indicating that the model established here can generate accurate and stable repeated target waves. By analyzing the lower part of Figure 23 (wave-current condition), the water wave free surface elevation distribution in the entire computational domain can be divided into two parts: left part (counter-current direction) and right part (coplanar-current condition). In both counter-current direction and coplanar-current condition, the generated waveform repeated well and no attenuation was found, which means that the model established here can be employed to numerically solve the wave-current interaction problems in deep water conditions and that the wave-current interactions will not distort the waveform. However, by comparing the left part with the right part, as the wave propagated against the current (left part), the wave amplitude was increased and the wavelength was reduced. On the contrary, the wave amplitude was decreased and the wavelength was increased as the wave propagated with the current (right part).Figure 23
Comparison between numerically and theoretically calculated water wave free surface elevation for pure wave (a) and wave-current interacting cases (b).
(a)(b)


## 4.3.4. Solitary Wave-Current Interaction

In some special marine environments, especially in the South China Sea, solitary waves occur frequently which may propagate toward coastal and near-shore regions and cause huge disasters. By interacting with riding tide current flows (coplanar currents), the wave speed of the propagating solitary wave may be vastly accelerated and the disaster responding time leaving for human being will be tremendously shorten although correspondingly the wave height may be decreased to a certain extent. On the other hand, the ebb tide current flow will enhance the influence of the solitary wave although with a lesser traveling speed. Here, the aforementioned numerical model proposed in this study is applied to investigate how the solitary wave interacts with the rising/ebb tide current flow. As an example, the water depth is chosen ash=8m, the wave height is chosen as H=1m, and the current velocity is chosen as ±1m/s. The numerical flume established here has an entire length of 800m and a total height of 10m. As for the pre-handling, the uniform mesh grid technology Δx=0.2m,Δy=0.004m is applied (containing one million quadrilateral cells). Figure 24 shows the effects of the current flow on the water surface profile at time, t=27s, by utilizing the comparisons of numerically calculated free surface under pure-wave and wave-current interacting (with coplanar and opposing currents) conditions. The numerically calculated wave heights, with the corresponding wave crest occurred at location x=−147.2m, x=−162m and162m for the case with an counter current U=−1m/s, no-current U=0m/s and following current U=1m/s are 1.25 m, 1m and 0.86 m respectively. It can be found that the change of the wave height (0.14 m, −14% of the height of the no-current solitary wave) induced by the following current is smaller than that (0.25 m, 25% of the height of the no-current solitary wave) induced by the opposing current. Meanwhile, the propagating velocity of the solitary wave is also altered by the ambient current flow. As for the opposing current case, the wave crest occurred x=−147.2m later than the current-free solitary wave case x=−162m. As for the following current case, the wave crest occurred x=178.6m earlier than the current-free solitary wave case x=162m.Figure 24
Effect of current velocity on the water surface profile of a solitary wave.


## 5. Conclusions

A CFD-based two-dimensional numerical wave and current tank (NWCT) with the utilization of an internal mass source wave-maker coupled with the analytical wave absorbing method, which does not interact with waves reflected from inside the computational domain, was used to numerically investigate the wave-current and wave-structure interaction problems. The deduced mass source function for wave generation and the corresponding momentum source function for wave absorption were embedded into the governing equations with RSM turbulence closure scheme and VOF algorithm to achieve wave generation and propagation in uniform current flow over constant water depth. The comparisons between numerically calculated results and the experimental measurements indicate that the approach proposed here is potentially suitable for solving the wave-current-structure interaction problems (except the large deformation problems), and the viability for a lengthy numerical simulation of large domain is achieved. Based on the aforementioned investigation, the following can be drawn:(1)
When wave is traveling with a following (opposing) current, the wave amplitude will decrease (increase) and the wavelength will increase (decrease). The velocity profiles will vary largely with the current flow and the surface wave motion. An increase in the ambient current flow velocity leads to a larger reduction in wave amplitude and a larger increase in wavelength. For the waves with the same period propagating in the same following current, a larger wave height will result in a larger reduction of wave height.(2)
The source region has an obvious influence on the accuracy of the wave generation and should be properly designed for solving the regular wave-current interaction problems’ meanwhile the source region is also affected by the desired wave parameters, ambient current flow velocity, and water depth. The designing method of the source region proposed here is shown to be capable of extending the ability of the wave-maker to simulate the water wave propagation in uniform current over constant water depth from shallow water conditions to deep water conditions.(3)
The model developed here is then applied to solving the problems of wave-current-structure interactions and generating steep waves. The numerically calculated results present fairly good agreement with the corresponding experimental data, indicating that the applicability and viability of the Navier-Stokes or RANS equations model for numerically handling coastal/ocean engineering problems is extended and the developed model here can be a useful and powerful tool for the study of wave-current-structure problems.(4)
The present work indicates that the model developed here with the utilization of coupled mass source wave-maker and analytical relaxation wave absorber can be employed under various wave conditions and that the propagation of a solitary wave on uniform currents over constant water depth can be simulated. As shown in the numerical results, the propagation of the solitary wave can be largely affected by ambient currents. The solitary wave height when propagating with opposing ambient current flow will be enlarged and decreased by the coplanar ambient current flow, and the propagating velocity will be reduced by opposing ambient current flow and enhanced by coplanar ambient current flow.



---

*Source: 8837476-2021-02-02.xml*