Numerical Simulation of the Interaction between Solitary Waves and Underwater Barriers Using a VPM–THINC/QQ-Coupled Model

Abstract

The interaction between solitary waves and underwater barriers is investigated using our in-house code, entitled VPM (volume-average/point-value multi-moment)–THINC/QQ (THINC method with quadratic surface representation and Gaussian quadrature)-coupled model. The stability and accuracy of the proposed model are validated by comparing the numerical results with those of the well established two-phase flow solver interFoam. All the results indicate that the presented coupled model has the advantage of high fidelity in simulating solitary wave propagation. Subsequently, solitary waves passing over a single underwater barrier are simulated by the present model. Numerical results are compared with experimental results in terms of the free surface elevation, velocity profile, vorticity field, and wave forces. Great agreements are obtained. In the end, the interactions between solitary waves and double underwater barriers are investigated numerically. The results reveal that the reflection coefficient increases first, and then decreases, with the increasing space between the two barriers. For cases with different wave heights, the transmission coefficient decreases monotonically, and the dissipation coefficient is opposed to the transmission coefficient.

1. Introduction

Due to their hydrodynamic similarities, solitary waves are often considered when studying the behavior characteristics of tsunamis on coastal effects [1,2,3]. Many scholars have studied the run-up height of solitary waves. Maiti et al. [4] indicated that the run-up height of solitary waves is closely related to the wave’s steepness and the slope of the plane. Wu et al. [5] found that the run-up height of breaking solitary waves increases with increases in the slope of the beach, but under the conditions of non-breaking waves, this relationship is the opposite. Li and Raichlen [6] obtained a nonlinear analytical solution to the classical shallow water equation using a hodograph transformation to describe the characteristics of solitary waves on a uniform plane beach. Boussinesq and shallow-water equation-based models are compared regarding their ability to predict solitary wave run-up height on plane beaches in [7]. Studying the behavior characteristics of tsunamis on coastal effects, Hsiao et al. [8] and Chang et al. [9] discussed the run-up process of solitary waves in a 300 m long flume. Moreover, Smith [10] used the boundary integral technique [11], combined with a boundary layer model, to compute the run-up and fluid flow for non-breaking waves and the velocity profiles at the lower region of the beach, which were well predicted.

In addition to the run-up height of solitary waves, the interaction between non-beach-type structures and solitary waves has also attracted a great deal of attention. The evolution of a solitary wave passing over a shelf is experimentally discussed by Seabra et al. [12]. The classic experiment of Seabra et al. [12] is a numerical reconstruction by different numerical models [13,14]. According to their results, it is noted that the surface elevation measured at the monitoring points can be well simulated, despite a certain degree of delay or advance in time. The deformation of solitary waves has been observed in detail in their studies under different physical conditions, including wave height and step size. Similar studies on different step sizes were carried out by Lin et al. [15]. Javier et al. [16] also explored the effect of a porous step on solitary wave breaking. Experiments and numerical simulations of solitary waves passing over a submerged permeable breakwater were implemented by Wu et al. [17]. Chang et al. [18] described vortex generation and evolution in water waves propagating over a submerged rectangular obstacle. Furthermore, Zhou et al. [19] used two submerged obstacles to explore the hydraulic characteristics of solitary waves. Wang et al. [20] showed the three-dimensional scattering of solitary waves passing over a vertical cylinder. A circular cylinder array was subsequently used to study solitary wave scattering in detail by [21]. Francesco et al. [22] analysed solitary wave-induced forces on horizontal circular cylinders by laboratory experiments and SPH simulations.

Since the development of numerical simulations, many numerical methods or models have been developed to predict the characteristics of solitary waves. These methods and models include a semi-implicit finite-volume shallow-water Boussinesq model, the SL-VOF method, the SPH method, and a Lagrangian finite-element Boussinesq wave model, among others [23,24,25,26,27,28,29,30,31,32]. Significant progress has been made in these areas. However, accurately simulating complex nonlinear phenomena, such as free surface deformation and wave breaking, is still a challenge for some numerical models.

In recent years, the open source software OpenFORM has been given widespread attention due to its convenient operation. InterFoam, which is a widely used solver within the OpenFOAM package, is based on the FVM (finite volume method)–VOF (volume of fluid) method. Deshpande et al. [33] evaluated the performance of the two-phase flow solver interFoam in OpenFORM and found that numerical dissipation is a major problem faced by interForm. For most existing FVM solvers, the difficulties in building higher-order stencils on unstructured grids limit the second-order spatial discretization, which leads to numerical decay and phase shifts in wave propagation [34] and over-dissipation in wave–structure interactions. To solve the above problem, a variant of muli-moment FVM, the VPM method, was proposed and developed in the OpenFOAM package by [35,36]. Unlike conventional high-order schemes that require a wide stencil, the VPM method provides quadratic reconstructions to the velocity field within a compact stencil by introducing the point values and partial derivatives of every control volume as additional constraints. Hence, it achieves a theoretically third-order accuracy on unstructured tetrahedral (or triangular) grids. Numerical decay was greatly improved with this method. The differences between the interFoam and VPM models are also introduced in [37].

In this work, the VPM scheme and THINC/QQ [38] method are employed to solve the flow field and interface capture, where the latter is claimed to have comparable solution quality to other existing VOF methods with PLIC geometric interface reconstruction. The numerical framework is named the VPM–THINC/QQ model. The classic experimental results of [39] are used, including the free surface, velocity profile, vorticity fields, and wave forces, to verify the reliability of the present model. Meanwhile, the performance of the present model is assessed by comparing the numerical results with the interFoam solver’s simulations. The generation of solitary waves and the influence of underwater structures on RTD coefficients of solitary waves are investigated in this study.

The main purpose of this paper is to develop a high-precision numerical model to achieve the high-fidelity propagation of solitary waves. The effects of double underwater barriers on the energy dissipation of solitary waves are also investigated. The rest of this article is arranged as follows. Governing equations, the flow solver of the present model, and a relaxation wavemaker method are described in Section 2. In Section 3, solitary waves are simulated by the present model and compared with the interFoam solver in the numerical wave tank without structure. Meanwhile, the numerical results of the present model are compared with the experimental and numerical results of [39,40] in terms of the free surface, velocity profile, vorticity field, and wave forces. The wave characteristics, including wave reflection and transmission and dissipation coefficients, are discussed in Section 4. Finally, the paper contains some conclusions in Section 5.

2. Numerical Model

2.1. Governing Equations

In this study, the incompressible two-fluid flows with moving interface are solved based on the one-fluid model [41]. The Navier-Stokes equations, containing the effects of surface tension and gravity, are used for both fluids in the same form. The large eddy simulation (LES) is adopted to deal with the turbulent flow in the vicinity of the submerged barriers. The governing equation are shown as:

$$\nabla ·\mathbf{u}=0,$$

(1)

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla ·(\rho \mathbf{u}\otimes \mathbf{u})=-\nabla p+(\nabla ·({\mu }_{\mathrm{eff}}\nabla \mathbf{u})+\nabla \mathbf{u}·\nabla {\mu }_{\mathrm{eff}})+{\mathbf{F}}_s-\mathbf{g}·\mathbf{x}\nabla \rho ,$$

(2)

$$\frac{\partial \varphi }{\partial t}+\nabla ·\left(\mathbf{u}\varphi \right)=\varphi \nabla ·\mathbf{u},$$

(3)

where $\mathbf{u}=(u,v)$ is the velocity vector with components u and v in x and y directions, respectively, p is the pressure in excess of the hydrostatic part, $\mathbf{g}$ is the gravity acceleration, $\rho$ the density, ${\mathbf{F}}_s$ is the surface tension force formulated by ${\mathbf{F}}_s=\sigma \kappa \nabla \varphi$ with $\sigma$ being the surface tension coefficient and $\kappa$ the interface curvature. ${\mu }_{\mathrm{eff}}$ is the effective viscosity coefficient, calculated as ${\mu }_{\mathrm{eff}}$ = $\mu$ + ${\mu }_t$. In this work, $\mu$ is the dynamic viscosity of the fluids, and ${\mu }_t$ is calculated by the standard Smagorinsky model [42] as:

$${\mu }_t=\rho {\left(C_s\mathsf{\Delta }\right)}^2\left|\mathbf{S}\right|,$$

(4)

in which $\mathsf{\Delta }$ is the dimension of grid, $C_s$ is the Smagorinsky constant (evaluated as 0.20), and S is the tensor of fluid strain rate.

A volume-of-fluid (VOF) function is used with an indicator function $\varphi (\mathbf{x}, t)$, distinguishing two kinds of fluids, shown as:

$$\varphi (\mathbf{x}, t)=\{\begin{array}{cc}1\hfill & \mathrm{water}\hfill \\ 0\hfill & \mathrm{air}\hfill \\ 0<\varphi <1\hfill & \mathrm{interface}.\hfill \end{array}$$

(5)

In the one-fluid model, it is assumed that the intrinsic fluid properties, such as density and viscosity, are updated based on the VOF function as,

$$\rho ={\rho }_1\varphi +{\rho }_2(1-\varphi ) \mathrm{and} \mu ={\mu }_1\varphi +{\mu }_2(1-\varphi ),$$

(6)

where ${\rho }_1$ and ${\rho }_2$ are the densities, and ${\mu }_1$ and ${\mu }_2$ are the dynamic viscosity coefficients of water and air, respectively.

2.2. VPM–THINC/QQ Model

The VPM–THINC/QQ model combines two newly developed numerical schemes, namely, the VPM scheme [35,36] and the THINC/QQ scheme [38]. The fluid dynamic equations are discretized by the multi-moment finite volume method, i.e., VPM scheme, and the VOF transport equation is solved by the THINC/QQ scheme. A numerical framework for interfacial multiphase flows on unstructured grids is employed to combine the two schemes. In this framework, the fraction-step method is used to update the numerical solution from time level n($t=t^n$) to $n+1$($t=t^n+\mathsf{\Delta }t$), and the solution procedure is summarized as follows.

• Update the velocity field from step n(${\mathbf{u}}^n$) to ${\mathbf{u}}^{*}$ by solving the VOF function (Equation (3)) and the advection term (Equation (7)) simultaneously in the time integration framework of the third-order TVD Runge-Kutta [43], shown as:

  $$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}=-\nabla ·(\rho {\overline{\mathbf{u}}}^n\otimes {\overline{\mathbf{u}}}^n).$$

  (7)
• Update the velocity field from ${\mathbf{u}}^{*}$ to ${\mathbf{u}}^{**}$ by solving the diffusion terms,

  $$\frac{{\mathbf{u}}^{**}-{\mathbf{u}}^{*}}{\mathsf{\Delta }t}=(\nabla ·(\mu \nabla {\mathbf{u}}^{*})+\nabla {\mathbf{u}}^{*}·\nabla \mu ).$$

  (8)
• Update the velocity field from ${\mathbf{u}}^{**}$ to ${\mathbf{u}}^{***}$ by adding the effects of surface tension and gravity force,

  $$\frac{{\mathbf{u}}^{***}-{\mathbf{u}}^{**}}{\mathsf{\Delta }t}=\frac{1}{\rho }(\sigma \kappa \nabla \varphi -\mathbf{g}·\mathbf{x}\nabla \rho ).$$

  (9)
• To make the intermediate velocity field ${\mathbf{u}}^{***}$ satisfy the mass conservation in Equation (1), it must be corrected by the following projection step. First, the pressure field at step $n+1$ is obtained by solving the Poisson equation,

  $$\nabla ·(\frac{1}{\rho }\nabla p^{n+1})=\frac{\nabla ·{\mathbf{u}}^{***}}{\mathsf{\Delta }t},$$

  (10)

  then, the velocity field is corrected by projecting the pressure field,

  $$\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^{***}}{\mathsf{\Delta }t}=-\frac{1}{\rho }(\nabla p^{n+1}).$$

  (11)

The numerical framework has already been built, and the high-fidelity propagation of regular waves has been realized. More detailed information about the numerical framework can be found in [37]. Next, the generation theory of solitary waves is introduced.

2.3. Wave-Maker with a Relaxation Region

For viscous numerical wave tanks, two kinds of wave makers have been widely used, including the piston wave maker [44,45,46] and the mass source wave maker [47,48]. In order to solve the secondary reflected waves generated during the interaction between solitary waves and structures, a relaxation wavemaker [49] is employed in this study. In the numerical wave tank, the free surface elevation and the velocity components of solitary wave are applied to the entrance boundary to generate a solitary wave. The expressions shown in the research of [14,50] are as follows.

$$\eta =H{\mathrm{sech}}^2\left(\sqrt{\frac{3H}{4h^3}}X\right),$$

(12)

$$u=\frac{\eta }{h}\sqrt{gh}\left[1-\frac{1}{4}\frac{\eta }{h}+\frac{h}{3}\left(\frac{h}{\eta }\right)\left(1-\frac{3}{2}\frac{y^2}{h^2}\right)\frac{d^2\eta }{d{X}^2}\right],$$

(13)

$$v=-\frac{y}{h}\sqrt{gh}\left[\left(1-\frac{1}{2}\frac{\eta }{h}\right)\frac{d\eta }{dX}+\frac{1}{3}{h}^2\left(1-\frac{1}{2}\frac{y^2}{h^2}\right)\frac{d^3\eta }{d{X}^3}\right],$$

(14)

where

$$X=-c\left(t,-,\frac{\tau }{2}\right),c=\sqrt{g(h+H)},\tau =\frac{2}{c}\sqrt{\frac{4h^3}{3H}}\left({\tanh }^{-1}\left(0.999\right)+\frac{H}{h}\right),$$

(15)

in which $\eta$ is the free surface elevation, u denotes the horizontal velocity, and v is the vertical velocity. h is the still-water depth, and H is wave height. c is wave celerity, and $99.9\%$ of the solitary wave can be generated during the duration $\tau$.

Here, a relaxation zone is implemented at the entrance to correct the solitary waves generated by the wave-maker boundary. By this way, the reflected wave in the computational domain to interfere with the wave-maker boundary could be avoided internally [49]. Another relaxation zone is also set up at the wave outlet to avoid wave reflection from the outlet boundary. The diagram of the wave flume is shown in Figure 1. A relaxation function

$$\begin{array}{cc}{\alpha }_R\left({\chi }_R\right)=1-\frac{\exp \left({\chi }_R^{3.5}\right)-1}{\exp \left(1\right)-1}& {\chi }_R\in \left[0,1\right],\end{array}$$

(16)

is applied inside the relaxation zone in the following way,

$$\psi ={\alpha }_R{\psi }_{computed}+\left(1-{\alpha }_R\right){\psi }_{target},$$

(17)

where $\psi$ is either $\mathbf{u}$ or $\varphi$. ${\psi }_{computed}$ represents the solution of the flow field, and ${\psi }_{target}$ denotes the value we expect. More details about the wave making method could be found in [49].

3. Assessment of the VPM–THINC/QQ

3.1. Comparison between the VPM–THINC/QQ Model and interFoam Solver

In this part, the performance of the VPM–THINC/QQ model in simulating the generation and propagation of solitary waves is compared with the interFoam solver. The setup of the computational domain is shown in Figure 2. The length of the empty wave tank from $x_0$ to $x_1$ is 18.0 m, and the height is 0.3 m. The still water depth is h = 0.14 m. A solitary wave with wave height of H = 0.07 m is generated. The wave-maker region with 1.5 m wide is at the leftmost part of the wave tank, and the wave damping region with 2.0 m wide is at the rightmost part of the wave tank. The same grid (grid1) is used in the VPM–THINC/QQ model and interFoam solver. The longitudinal grid is uniform along the wave tank, with a size $dx$ of 0.03 m. The minimum vertical grid $d{y}_{min}$ is 0.003 m. Besides, a finer mesh (grid2) is also used in interFoam solver. The mesh is divided into 1800 grids, evenly distributed in the longitudinal direction, with the smallest vertical grid size $d{y}_{min}$ of 0.002 m. The initial time step is set as $\mathsf{\Delta }t$ = 0.001 s, which would be automatically updated according to the maximum Courant-Friendrichs-Lewy (CFL) number, 0.5. The monitoring points a, c, and e are placed at x = 6, 8, and 10 m, respectively, to measure the variation of water surface. The total simulation time $t_{total}$ is 12 s.

The snapshots of solitary wave propagation at different times are shown in Figure 3. The simulated results of the VPM–THINC/QQ model, based on grid 1, are compared with those of the InterFoam solver simulations with grid 1 and grid 2. As shown in Figure 3, it can be seen that, under the same grid size, the wave height of solitary waves using the VPM–THINC/QQ model remains constant with increasing time at the same grid size, while the wave height of solitary waves simulated by interFoam solver decrease along the water tank. The wave height of solitary waves at t = 11 s is 0.721 times of the incident wave height. It indicates that the VPM–THINC/QQ model has the advantage to maintain more waveform than the interFoam solver in the solitary wave propagation process. The time series of solitary waves at three monitoring points are shown in Figure 4. The wave heights at x = 6, 8, and 10 m are 0.833, 0.794, and 0.767 times of the original wave height in interFoam solver, respectively. For the VPM–THINC/QQ model, the computed wave heights agree excellently with the original wave heights, which shows that the VPM–THINC/QQ model can effectively eliminate the numerical consumption and achieve high-fidelity propagation of solitary waves. Grid 2 is also used to simulate the solitary wave in empty wave tank by interFoam solver. From Figure 3 and Figure 4, the solitary wave heights are attenuated to about 0.915 times of the original wave heights at t = 11 s. The attenuation of interFoam solver has been significantly improved due to mesh refinement. Moreover, the wave heights at x = 6, 8, and 10 m become 0.962, 0.957, and 0.954 times the original wave height, respectively. However, there is still numerical dissipation compared with the results obtained by VPM–THINC/QQ model, which even uses coarse mesh (grid1). Obviously, the VPM–THINC/QQ model has high fidelity in simulating solitary waves.

3.2. Verification of the VPM–THINC/QQ Model

In this part, the performance of the present model in simulating solitary waves, interacting with the underwater barrier, is verified according to the experiment of [39]. The available experimental data and the other numerical results in [39,40] are employed to compare with the results simulated by the present model. The numerical computation domain is shown in Figure 5, which sets the same conditions as Figure 2. The difference is that an underwater barrier is placed in the numerical wave tank. Its central position is at x = 8.0 m, width is W = 0.02 m, and height is d = 0.10 m, respectively. At the same time, five wave monitoring points are placed in the wave flume, and the position of the five wave gauges are at $x_a$ = 6.0 m, $x_b$ = 7.333 m, $x_c$ = 8.0 m, $x_d$ = 8.347 m, and $x_e$ = 10 m, respectively.

In this case, the solitary wave with wave height of H = 0.07 m is simulated. The free surface variation, the distribution of velocity, and the vortex field are discussed. Following this, the wave forces on the structure subjected to solitary waves with different wave heights are analyzed. In order to expediently compare with the results of [39,40], the time t = 0 s is defined when the crest of the incident wave reaches the position of wave gauge $x_b$.

3.2.1. Free Surface

The time series of free surface variation during solitary wave passing the underwater barrier are shown in Figure 6. The free surface variation in x = 7.333 m, 8.0 m, and 8.347 m are depicted as Figure 6a–c, respectively. The simulated results are compared with the experimental data [39], and the numerical results are simulated by the Reynolds-Averaged Navier–Stokes (RANS) based model [39], as well as the constrained interpolation profile (CIP) based model [40]. From Figure 6, it can be seen that the free surface variation simulated by VPM–THINC/QQ model shows a satisfactory agreement with the experimental measurements and other numerical results on the whole. It should be noted that data observed by particle image velocimetry (PIV) equipment is only plotted in Figure 6b because the view field of PIV system is focused near the underwater barrier. In Figure 6a, the second soliton is clearly observed due to wave reflection. Besides, wave fluctuations caused by wave breaking are found in Figure 6b,c. A similar phenomenon is also observed by [12] when solitary waves pass over the underwater step.

The free surface elevations during solitary wave passing the underwater barrier at t = 0.46, 0.60, 0.74, 0.88, and 1.02 s are shown in Figure 7. The numerical results of the present model (solid line) are compared with the experimental data measured by PIV equipment in [39], which are shown in circle. In Figure 7, it can be seen that the numerical results of the present model keep consistent with PIV’s observations. In their entirety, there are some stages during solitary wave propagating over the underwater barrier, including waves’ “crest-crest exchange”, becoming steep, rolled up, and breaking, respectively. These stages have also been explicitly described in [39,40]. In Figure 7d, some air is entrapped into water, so some water and air mix together (Figure 7e). It is worth mentioning that only one bubble is observed in the present model and the RANS model, which is shown in Figure 7 of [39], while multiple bubbles are observed in Figure 3e of [40]. The explanation for this phenomenon may be that the bubbles burst easily because the surface tension is not taken into account in the CIP-based model. Overall, the present model could simulate the free surface during solitary wave passing the underwater barrier primely in both space and time.

3.2.2. Velocity Distribution and Vorticity Field

In Figure 8, the velocity distributions at t = 0.46, 0.60, 0.74, 0.88, and 1.02 s are shown as Figure 8a–e, respectively. The velocity profiles from x = 8.06 m to x = 8.20 m, with an identical interval of 0.02m, are depicted. The simulated results are compared with PIV’s measured data, as well as numerical results [39]. From Figure 8, it can be seen that there are velocities above the water surface only found by the present model. The present model is a multiphase flow model in which air is included in fluid calculation. While the RANS model in [39] solves single-phase flow problem and PIV system obtains velocity distributions by monitoring the tracer particles in water, so the air velocity is ignored in both methods. In Figure 8e, the horizontal velocity profile of the present model at x = 8.18 m shows a satisfactory agreement with experimental data, while there are some deviations for the RANS model in [39], expecially at x = 8.14 m and x = 8.18 m. In general, the VPM–THINC/QQ model has an advantage of high accuracy in terms of the velocity distributions.

The vortex generation and evolution during the strong interaction between solitary waves and the underwater barrier are shown in Figure 9. The experimental data of [39] (placed at left) and the VPM–THINC/QQ model results (placed at right) at t = 0.46, 0.60, 0.74, 0.88, and 1.02 s are shown as Figure 9a–e, respectively. The vorticity fields are computed by the operation of ($\partial v/\partial x-\partial u/\partial y$). The negative and positive values denote the clockwise and anticlockwise vortices, respectively. From Figure 9a–e, it is obvious that the vorticity fields above the free surface are observed in the present multiphase flow model, but not in experiment [39]. The vorticity fields in air at t = 0.46, 0.60 s are anticlockwise near water surface, but, at t = 0.74, 0.88, and 1.02 s, they appear partly clockwise. The reason for this phenomenon is that water surface begins to move downstream, and the part of the water surface breaks up gradually due to the underwater barrier’s effect. The vortices in water of the present model are basically consistent with the measurements of PIV equipment [39], and the vorticity field of bubble in water is anticlockwise, which also agrees well with the experiment.

3.2.3. Wave forces

The horizontal force acting on the left and right sides of the underwater barrier by solitary waves is shown in Figure 10a. The simulated results are compared with other numerical results computed by the RANS model and the CIP-based model. As shown in Figure 10a, the maximum force on the left side of the barrier is earlier than that on the right side. Meanwhile, the force on the left side is greater than that on the right side. The total horizontal force on the underwater barrier is calculated by subtracting the right force from the left force. The total horizontal force is shown in Figure 10b. In Figure 10, a discrepancy is observed in the tailing stage of the time series due to the complex flow field after wave breaking. In general, the time series of wave forces fit well with other published numerical results.

In order to obtain the universality of the present model, the cases of $H/h$ = 0.1, 0.15, 0.2, 0.35, and 0.5 are simulated by the present model, where H is the solitary wave height, and h is the still water depth. The time series of the total horizontal forces acting on the underwater barrier are plotted in Figure 11a. It can be seen that, for the cases of small wave height, the maximum force appears later because of the low phase velocity. The maximum horizontal forces under different wave heights are compared with that of [39,40], which are shown in Figure 11b. In Figure 11b, the crest of total horizontal force on the barrier increases linearly with the increase in wave height. The phenomenon is also referenced by [51].

In this part, the interaction between solitary waves and the underwater barrier is simulated by the VPM–THINC/QQ model. In terms of the free surface variation, the distribution of velocity, vortex and wave forces, the simulated results of the present model are in good agreement with experimental and other numerical results. On the whole, the ability of the VPM–THINC/QQ model to simulate the interaction between solitary waves and the underwater barrier has been fully verified and assessed.

4. Application of the VPM–THINC/QQ Model

How to effectively reduce wave energy damage is an important issue in coastal protection. From the above study, it can be found that an underwater barrier can consume a portion of the energy of solitary waves to reduce their impact on the coast. In this section, a newly added barrier (see Figure 12) is placed in the downstream of the original barrier in order to achieve more wave energy dissipation. Note that the newly added barrier is the same as the original barrier. The width of the two barrieres are W = 0.02 m, and height are d = 0.10 m, respectively. In Figure 12, $\delta$ is the space between the two underwater barriers, h is the still water depth. At the same time, two wave monitoring points, a and e, are placed in the wave flume.

In order to quantitatively study the influence of two underwater barriers on wave energy dissipation, the reflection coefficient, transmission coefficient, and dissipation coefficient of solitary waves are analyzed. Shao [52] defines the coefficient of transmission as $C_t=H_t/H_i$, reflection as $C_r=H_r/H_i$, and dissipation as $C_d=H_d/H_i$, where $H_t$, $H_r$, $H_d$, and $H_i$ are the transmitted, reflected, dissipated, and incident wave heights, respectively. As described in [40], the incident wave height $H_i$ and the reflected wave height $H_r$ correspond to the maximum value and the second largest value of the wave surface at x = 6.0 m, respectively. In addition, the largest value of the wave surface at x = 10.0 m denotes the transmitted wave height $H_t$. In this study, $C_r$, $C_t$, and $C_d$ are totally named as the RTD coefficients, and the relationship between them can be established according to energy conservation law in Equation (18). It should be noted that $C_r$ and $C_t$ are computed by $H_i$, $H_r$, and $H_t$, but $C_d$ is obtained according to the energy conservation law.

$${C_t}^2+{C_r}^2+{C_d}^2=1.$$

(18)

The space between two underwater barriers $\delta$ and the wave height of solitary waves H are selected as important parameters for the RTD coefficients to explore the interaction between solitary waves and dual underwater barriers. The space $\delta$ = 0.5–6.0 d, with the equal interval of 0.5 d between the two underwater barriers and solitary waves with $H/h$ = 0.1, 0.15, 0.2, 0.35, and 0.5, are chosen for numerical simulation. The total of 60 cases is shown in

Table 1. The characteristics of cases.

Cases	$\mathit{H}/\mathit{h}$ = 0.10	$\mathit{H}/\mathit{h}$ = 0.15	$\mathit{H}/\mathit{h}$ = 0.20	$\mathit{H}/\mathit{h}$ = 0.35	$\mathit{H}/\mathit{h}$ = 0.50
$\delta /d$ = 0.5	case1a	case2a	case3a	case4a	case5a
$\delta /d$ = 1.0	case1b	case2b	case3b	case4b	case5b
$\delta /d$ = 1.5	case1c	case2c	case3c	case4c	case5c
$\delta /d$ = 2.0	case1d	case2d	case3d	case4d	case5d
$\delta /d$ = 2.5	case1e	case2e	case3e	case4e	case5e
$\delta /d$ = 3.0	case1f	case2f	case3f	case4f	case5f
$\delta /d$ = 3.5	case1g	case2g	case3g	case4g	case5g
$\delta /d$ = 4.0	case1h	case2h	case3h	case4h	case5h
$\delta /d$ = 4.5	case1i	case2i	case3i	case4i	case5i
$\delta /d$ = 5.0	case1j	case2j	case3j	case4j	case5j
$\delta /d$ = 5.5	case1k	case2k	case3k	case4k	case5k
$\delta /d$ = 6.0	case1l	case2l	case3l	case4l	case5l

, where h is the still water depth, and d is the height of underwater barrier.

4.1. The RTD Coefficients with Single Underwater Barrier

In Section 3.2, the interaction of solitary waves with a single underwater barrier are simulated by using the VPM–THINC/QQ model. The RTD coefficients during the process of solitary waves passing over a single underwater barrier are calculated and plotted in Figure 13. As shown in Figure 13, the coefficient of wave reflection $C_r$ increases with the increase in wave height for values, ranging from 0.10 to 0.20, and then remains basically unchanged. The wave transmission coefficient $C_t$ decreases as the wave height increases for values, ranging from 0.10 to 0.20, and then slowly increases for non-linearity values of waves, ranging from 0.20 to 0.50. The trend of wave dissipation coefficient $C_d$ is opposite to that of wave transmission coefficient $C_t$. It is worth pointing out that, when the wave height of solitary wave is 0.028 m, the maximum energy dissipation occurs, and the dissipation coefficient is 0.504.

4.2. The RTD Coefficients with Double Underwater Barriers

Figure 14 shows the variation of RTD coefficients during the process of solitary waves interacting with the double underwater barriers. As shown in Figure 14, the changing trends of the RTD coefficients for solitary waves with different wave heights are the same as the cases with different spaces between double underwater barriers in general. The wave reflection coefficient $C_r$ increases rapidly and then decreases slowly as the space between two barriers increases.

The coefficient of wave transmission $C_t$ decreases monotonously with increasing space between two barriers. The trend of variation in wave dissipation coefficient $C_d$ is opposite to that of wave transmission coefficient $C_t$. In Figure 14a, the transmission coefficients of double underwater barriers $C_t$ are smaller than that of the single underwater barrier from $\delta /d$ = 1.0 to 6.0 for $H/h$ = 0.1. While in Figure 14b–e, the transmission coefficients of double underwater barriers $C_t$ are larger than that of single underwater barrier between $\delta /d$ = 1.0 and $\delta /d$ = 1.5 for $H/h$ = 0.15, 0.2, 0.35, and 0.5. The reason for this phenomenon is that, as the wave height increases, the transmission coefficient of the double underwater barrier requires more space to reach that of a single underwater barrier. As the wave height increases, there needs to be more space to enhance the dissipation capacity of the double underwater barrier compared to a single barrier. When the distance between the double barriers is small, the wave reflection increases. The reflection coefficient occurs at maximum at $\delta /d$ = 1.0 for $H/h$ = 0.1 in Figure 14a, but for $H/h$ = 0.5 at $\delta /d$ = 2.0 in Figure 14e. This indicates that the higher the wave height, the greater the space between the dual underwater barriers required to achieve maximum wave reflection. As the space between the double submerged barriers increases, the impact of the second barrier on wave reflection is relatively small. For the same wave height, the reflection coefficient gradually becomes the same as that of the single barrier. In Figure 14d, the wave transmission coefficients for $H/h$ = 0.35 is reduced from 0.856 to 0.737 by placing another barrier downstream of the barrier at x = 8.5 m. In Figure 14e, the wave transmission coefficients for $H/h$ = 0.5 is reduced from 0.859 to 0.760 by placing another barrier downstream of the barrier at x = 8.45 m, which could provide a reference for coastal protection.

The vorticity distributions of the maximum reflection coefficient for each wave height are shown in Figure 15. Figure 15a–e depict the $H/h$ = 0.10 and $\delta /d$ = 1.0 (case1b), $H/h$ = 0.15 and $\delta /d$ = 1.5 (case2c), $H/h$ = 0.20 and $\delta /d$ = 1.5 (case3c), $H/h$ = 0.35 and $\delta /d$ = 1.5 (case4c), and $H/h$ = 0.5 and $\delta /d$ = 2.0 (case5d), respectively. As shown in Figure 15, the vortex shedding from the first barrier is blocked by the second barrier, which causes more reflection of the wave surface. The resulting vortices move between barriers and dissipate energy. As a whole, the higher the wave height, the higher the reverse motion of the wave surface. In Figure 15e, the range of vorticity is large because the maximum wave height has a large amount of energy. The increasing space between double underwater barriers can not effectively form wave surface reflection, so there is less energy reflected back. However, it provides more space for vortices in the water to move and dissipate more energy. As the distance between the two barriers increases, the energy of the wave reflection from the second barrier partially dissipates before reaching the first barrier, and it ultimately completely dissipates. At the same time, the wave dissipation between the barriers reaches saturation and does not generate more dissipation. The transmitted wave energy no longer decreases and tends to stabilize.

5. Conclusions

In this study, an in-house numerical framework (the VPM–THINC/QQ model) of the VPM scheme for solving the flow field and the THINC/QQ method for interface capture is established. The propagation of solitary waves and their interactions with underwater barriers are investigated by the present model. The solitary wave is generated by using a relaxed wave-maker method. The well established interFoam solver is employed to comparatively verify the accuracy and applicability of the present model. With the same coarse grid setting (grid1), the wave height at t = 11 s has become 0.721 times of the original wave height in interFoam solver, while remaining almost unchanged in the present model. Compared to the attenuation phenomenon of the simulation results of the interFoam solver, which cannot be eliminated by encrypting the mesh, the present model is able to simulate the generation and propagation of solitary waves stably and accurately.

The solitary wave passing over a single underwater barrier is simulated and studied based on the present model. Good agreements are obtained by comparing with experimental data and other numerical simulation results in the aspects of the free surface, velocity distribution, vorticity field, and structural forces. All comparisons demonstrate the excellent accuracy of the coupled VPM–THINC/QQ model in predicting the interaction characteristics of solitary waves and underwater barriers.

On the basis of the adequately validated and verified VPM–THINC/QQ coupling numerical model, this paper simulates and studies the effects of dual underwater barriers on the reflection, transmission, and dissipation coefficients of solitary waves. The results show that, as the space between two underwater barriers increases, the reflection coefficient increases first and then decreases, the transmission coefficient decreases monotonously, and the dissipation coefficient is opposite to the transmission coefficient for different wave heights. Moreover, it is found that the wave reflection would be enhanced as the number of underwater barriers increases. Besides, the proper space between the double underwater barriers would play an important role in the generation of vortexes in different directions, which is conducive to the dissipation of wave energy. All the above results indicate that the present model can serve as a powerful tool to simulate the generation and propagation of solitary waves and their interaction with ocean structures. In the future, a large number of numerical cases will be carried out by the present model to simulate the real events in coastal engineering.