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.
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.where represents time and represents the fluid density, represents the abscissa and ordinate, the subscripts i and j represent the indices used in Einstein’s sum convention, denotes the ensemble mean velocity component, represents the fluid pressure, represents the dynamic viscosity, and denotes the gravitational acceleration. represents the additional mass source term specified for target wave generation, represents the additional momentum source term designated for wave absorption. 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: indicates that the computational cell is fully occupied by water phase, while corresponds to the computational cell full of air phase. Computational cells with volume function value of 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):
in which is the water wave phase velocity of the target wave. is the area of the source region. is the water wave surface elevation; for linear monochromatic wave, it can be represented as , and for solitary wave, it can be represented as the following formula:
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:where the subscript represents the renewed physical quantity in the designated domain—the wave absorption domain, the subscript represents the computed physical quantity at the previous calculational time step, denotes the target current velocity, represents the relaxation function, , , 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:
By subtracting equation (7) from (9) and equation (8) from (10), the corresponding momentum source function in abscissa and ordinate can be deduced as
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.
Figure 1 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 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.
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 , width of , and maximum depth of . 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 . For all the investigations performed in the channel, the water depth (h) was kept at and the target wave period was chosen as ; 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 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 − to in abscissa and to in 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 , wave amplitude , and wave period is generated in the NWCT at a constant water depth of . 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, ) and three different temporal grid sizes (0.01 (T/100), 0.005 (T/200), and 0.002 (T/500)). The numerical results are presented in Figure 2(a). It is observed that the calculated time series of surface elevation at location for 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 , ), medium (aforementioned grid system), and fine grid systems (containing 1440000 quadrilateral cells with, nonuniform in ordinate with ). 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.
(a)
(b)
Figure 3 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 , and , 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.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 4 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.
(a)
(b)
(c)
Figure 5 displays the comparison of simulated and theoretical free surface elevation distribution along the wave tank for wave-current cases (from to for opposing current domain and from to for coplanar current domain). The theoretical free surface elevation is estimated by , in which can be solved by . 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 to (for following current condition and pure-wave condition), (for opposing current condition) at 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.
(a)
(b)
(c)
Figures 7–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 and/or ) 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 and/or ) 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.
(a)
(b)
(a)
(b)
(a)
(b)
(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 region fixed at (-the spatial grid step in the abscissa), the height of the source region fixed at, and with the center of the region varied. The distance 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 (, where is the numerically calculated wave amplitude and 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 or to the free surface , the higher is the inaccuracy. For , the inaccuracy of the wave-current cases is bigger than the wave-only cases. For , 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 , i.e., the vertical distance 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 .
(a)
(b)
(c)
Then, the influence of the width 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 (, , , , and ) 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 has little influence on wave generation. During the simulation process, it is noted that a smaller width of the source region 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 and the height 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 of the source region as long as it is larger than the wave height, and the optimal height of the source region is about 10% of the water depth.
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 as , , , , , and respectively, and the input wave parameters were chosen as: wave height , wave period , and the average water depth was fixed at . The averaged wave profiles generated in the working domain, to (wave propagating with coplanar current) and to (wave propagating with opposing current) at to were used to compute the averaged wavelength and wave amplitude . The statistical wavelength , wave amplitude , and the corresponding change rate ( and ) are presented in Table 2, in which , ( and 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.
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 and wave period is generated at a constant still water depth 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 ( to ) and up to 0.5 m deep. The internal mass source wave-maker is located in the region with abscissa ranging from to and ordinate ranging from to . 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 , , , , , and , respectively. The entire domain is discretized using nonuniform structured grid system, with in the range of and in the range of (137 MB, containing 1470000 quadrilateral cells).
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 Figure 15, 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.
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 of . 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.
(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 from are 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].
(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 as , 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 due to the wave scattering caused by the coexisting current, especially when on the top of the breakwater (, , ). 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.
(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:
A validation computation of a Stokes-wave train was carried out based on the following chosen parameters: the incoming wave height being set to be , the corresponding wave period being , the length of the numerical tank ranging from to , the total depth being , still water depth chosen as , 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 and the instantaneous phase volume fraction distribution at time , 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.
(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:
Figure 20 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.
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 of ranges from to , and the still water depth is chosen and kept at . The target wavelength and the wave height are equal to and , respectively, and the current velocity is chosen to be . The internal mass source wave-maker is located in the region with the abscissa ranging from toto , and the ordinate ranging from to 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 , , 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).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
As shown in Figures 21 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 , . Meanwhile, the wave propagates faster in coplanar-current direction than in counter-current 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 Figure 23, 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).
(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 as , the wave height is chosen as , and the current velocity is chosen as . The numerical flume established here has an entire length of and a total height of . As for the pre-handling, the uniform mesh grid technology is applied (containing one million quadrilateral cells). Figure 24 shows the effects of the current flow on the water surface profile at time, , 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 , for the case with an counter current , no-current and following current 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 later than the current-free solitary wave case . As for the following current case, the wave crest occurred earlier than the current-free solitary wave case .
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.
Data Availability
The physical experimental results data supporting this article are from previously reported studies and datasets, which have been cited. The processed data are available from the corresponding author of the cited article [1] upon request. The numerical cases data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors thank Prof. Umeyama Motohiko for sharing the experimental data in Figures 3 and 10. The authors also thank Prof. Yen-Lung Chen for sharing the experimental data in Figure 16. The authors are grateful for the grant of CSC-China Scholarship Council, the Fundamental Research Funds for the Central Universities (China) (Nos. 3132020116 and 31320210123), and NSFC (Project no. 51409032) (National Natural Science Foundation of China).