CFD Simulations of the Propagation of Free-Surface Waves Past Two Side-By-Side Fixed Squares with a Narrow Gap

Abstract

Two-dimensional computational fluid dynamics (CFD) simulations are carried out to investigate the gap resonance phenomenon that occurs when free-surface waves travel past twin squares in tandem. The volume of fluid method is used to capture the free surface. Validation studies of the present numerical model are conducted for different incident wave frequencies. The numerical results agree well with the published experimental data in terms of the free surface elevations in the gap. The hydrodynamic characteristics of the water column in the gap are investigated at different incident wave frequencies and gap widths. It is found that the free surface elevation in the gap increases and then decreases with the increasing incident wave frequency. The horizontal force on the weather side square structure (the structure in front of the gap) reaches the peak value at a larger frequency than the gap resonance frequency, whereas the variation of the horizontal force on the lee side structure (the structure behind the gap) is in-phase with the free surface elevation in the gap. Moreover, the added mass of the water column in the gap increases with the increasing gap width, which results in the decrease of resonance amplitude and frequency in the gap. However, this does not necessarily reduce the peak value of the horizontal forces on the structures.

1. Introduction

In a multiple floating body system, such as the ship-by-ship offloading system and the assembly of very large floating structures (VLFS), the free surface elevations in the narrow gap between structures can be much higher than the incident wave heights. This phenomenon is called gap resonance, and it is caused by the proximity of incident wave frequencies to the natural frequency of the fluid vibrations in the gap. The gap resonance may cause higher wave forces, lead to violent body motions, and influence the stability of the system.

The sketch of the gap resonance motion is shown in Figure 1. The water column in the gap can have more than five times higher free surface elevations as compared to the incident waves, according to the experiments of Saitoh et al. [1]. The water column in the gap can be simplified as a rigid body (Molin, [2]). The motion of this rigid body then can be solved by considering the gravity force G, buoyancy force $F_b$, wave excitation force $F_w$, and radiation force $F_r$. Therefore, the oscillation of the water column in the gap can be influenced by these forces, which are relevant to the dimension of the gap and the incident wave conditions.

Extensive theoretical analyses, numerical simulations, and experiments have been carried out to investigate the generation mechanism and hydrodynamic characters of this phenomenon. The potential flow theory is an efficient method to investigate gap resonance. Theoretical explanations have been made to figure out the generation mechanism of the narrow gap, and simplified solutions have been proposed to estimate the natural frequency of the gap when the gap width is relatively small (Saitoh et al., [1]; Molin, [2]; Liu and Li, [3]; Tan et al., [4]). The estimated two-dimensional (2-D) gap resonance frequencies agree well with experiments (Saitoh et al., [1]; Iwata et al., [5]; Cong et al., [6]). However, since the linear potential flow theory over-predicts experimentally determined resonant responses in the gap under the excitation of sinusoidal waves (Faltinsen et al., [7]; Kristiansen et al., [8]), most of the potential flow models are semi-analytical or combined with empirical terms (Huijsmans et al., [9]; Newman et al., [10]; Chen, [11]).

As one of the sources of the damping terms in the potential flow model, the viscous effect is inevitable while solving the hydrodynamic forces and the resonance amplitudes at the gap between structures. Faltinsen and Timokha [7] and Kristiansen and Faltinsen [12] systematically investigated the damping in the gap both by experimental measurement and the potential flow method. The damping effects of the free surface nonlinearity, the flow separation, and the boundary layer on the hull were evaluated. It was concluded that the flow separation from the entrance of the gap is the main source of the discrepancy between the potential solver and their experimental results. Furthermore, Kristiansen and Faltinsen [13] investigated the coupled ship motion and the piston-mode flow in a gap. The flow separation at the ship bilge was modeled by an inviscid vortex tracking method. In this manner, the potential flow results were effectively suppressed. They further put forward a combined potential and viscous solver to deal with the flow separation from the barge bilge (Kristiansen and Faltinsen, [13]). Feng and Bai [14] investigated the gap resonance under nonlinear waves by a mixed Eulerian–Lagrangian scheme by the high-order boundary element method (HOBEM). They demonstrated that the main source of the overestimated response is viscosity instead of the free surface nonlinearity. The effect of free surface nonlinearity is negligible except for slightly increasing the resonance frequency. Another truth that proves the importance of viscosity is the influence of the hull bilge shape on the resonance amplitude. Moradi et al. [15] investigated a 2-D gap resonance through numerical simulations including viscosity. They found that the resonance amplitude became smaller with the increasing bilge radius. Tan et al. [16] carried out large quantities of experiments to investigate the effects of the bilge shape. Their results suggested that the energy dissipations of sharp bilge cases are larger than round bilge cases. This was likely due to the different flow separation behavior at the bilge.

Several computational fluid dynamics (CFD) simulations have been carried out to evaluate the viscous dissipation of narrow gap resonance problem. Lu et al. [17,18] investigated the 2-D gap resonance both by the modified potential model [11] and viscous numerical methods. By optimizing the artificial damping of the potential solver, the numerical results based on the potential theory agreed well with the turbulent viscous solver in terms of the resonance amplitude and the hydrodynamic forces. Tan et al. [16] figured out the relationship between the damping coefficients of a theoretical dynamic model ($\epsilon$) and that used in the modified potential model [11,17] (${\mu }_p$): ${\mu }_p=3\pi \epsilon {\omega }_n/8$, where ${\omega }_n$ is the natural frequency, and the damping coefficient $\epsilon$ could be calculated according to the surface elevation results of the numerical viscous solver. Zhao et al. [19] experimentally investigated the 3-D gap resonance under new wave-type transient waves. They found that each gap resonance mode had a characteristic damping which was somewhat larger than the damping calculated using the linear potential flow theory alone. Later, Wang et al. [20] employed OpenFOAM-based Navier–Stokes (N–S) equations to reproduce the experiment results of Zhao et al. [19]. The grid resolution, mesh topology, domain size, and boundary conditions were systematically optimized. The transient wave group was considered to be a better choice for investigating the 3-D gap resonance phenomenon as compared to the regular incident waves. Based on a CFD solver, Chua et al. [21,22] developed a framework to evaluate the damping coefficients of a 3-D gap resonance problem. The energy dissipations regarding wave scattering, frictional force, flow separation, appendages, and hull motions were investigated thoroughly. It could be concluded that the modified potential model considering the damping effects is able to predict well regarding the gap resonance problem. The value of damping, however, still relies on CFD simulations, which are less expensive compared to experiments.

The previous research concerns on the gap resonance problem were usually the viscous effects and the free surface elevations in the gap. There is still lack of research on the characteristics of hydrodynamic loads on structures with a narrow gap. In the present study, the hydrodynamic loads, viscous effects, and free surface elevations related to the narrow gap resonance were investigated thoroughly. The forces on the floating structures are more important as compared to the surface elevations in the gap, while the horizontal wave loads are more sensitive to the gap resonance, as compared to the vertical wave loads [18]. The present study is organized as follows. Firstly, a mesh and time-step refinement study is conducted for the case with the smallest gap width, which is considered to be the most severe case. Secondly, the free surface elevations in the gap, the viscous dissipation, and the horizontal wave loads on the structures are studied at different incident wave frequencies and different gap widths between structures. Meanwhile, viscous dissipation and horizontal wave loads versus the gap width are also discussed. The present research can provide a reasonable reference to the engineering design process in term of nap gap damping.

2. Numerical Modeling and Setup

2.1. Governing Equations

The open source CFD software OpenFOAM [23] was used to solve the N–S equations numerically. The toolbox waves2Foam (Jacobsen et al., [23]) was applied for the numerical wave tank. The governing equations for the incompressible viscous flow are given as follows [24]:

$$\{\begin{array}{c}\nabla ·\mathit{v}=0\\ \frac{\partial \left(\rho \mathit{v}\right)}{\partial t}+\nabla ·\left(\rho \mathit{v}\mathit{v}\right)=\nabla ·\mu \nabla \mathit{v}-\nabla p_{rgh}-\mathit{g}·\left(\mathit{x}-{\mathit{x}}_r\right)\nabla \rho \end{array}$$

(1)

where $\nabla$ denotes the Hamiltonian operator, $\mathit{v}$ is the velocity field, $\rho$ is the fluid density, $\mu$ is the dynamic viscosity coefficient, $t$ is the time, $p_{rgh}$ is the pressure in excess of the hydrostatic pressure $p_{rgh}=\rho -\rho \mathit{g}·\left(\mathit{x}-{\mathit{x}}_r\right)$, $\mathit{g}$ is the gravity vector, $\mathit{x}$ is the Cartesian coordinate vector, and ${\mathit{x}}_r$ is the reference coordinate vector defined at sea level. The contribution of surface tension effect is less than 1% of the inertial force when T > 0.35 s, H > 0.02 m [24], where T is the deep-water wave period and H is the deep-water wave height. Therefore, the surface tension effect is considered negligible in the present study (see the numerical set-up in Section 3).

2.2. Free Surface Capturing

The volume of fluid (VOF) method was applied to capture the free surface. The water volume fraction is defined as follows:

$$\mathsf{\alpha }\left(\mathit{x},t\right)=\{\begin{array}{cc}0& air\\ 0<{\alpha }_s<1& freesurface\\ 1& water\end{array}$$

(2)

where α is the volume fraction of the water phase. The density $\rho$ and the dynamic viscosity $\mu$ are calculated as follows:

$$\{\begin{array}{c}\rho =\alpha {\rho }_{water}+\left(1-\alpha \right){\rho }_{air}\\ \mu =\alpha {\mu }_{water}+\left(1-\alpha \right){\mu }_{air}\end{array}$$

(3)

The transport equation of the water volume fraction is:

$$\frac{\partial \alpha }{\partial t}+\nabla ·\left(\alpha \mathit{v}\right)+\nabla ·\left[\alpha \left(1-\alpha \right){\mathit{v}}_r\right]=0$$

(4)

where $t$ is the time; $\alpha =\alpha \left(\mathit{x},t\right)$ denotes the cell-based water volume fraction at coordinate $\mathit{x}$ in time $t$; $\mathit{v}$ is the local fluid velocity; and ${\mathit{v}}_r$ is the compress velocity at the interface [25]:

$${\mathit{v}}_r={\mathit{v}}_{water}-{\mathit{v}}_{air}$$

(5)

2.3. The Numerical Wave Flume

The relaxation zone technique [24] was applied for making and absorbing waves. In each time-step, this technique corrects fluid fields in relaxation zones by Equation (6). As a consequence, fluid fields in these zones change gradually from the computed fields (that are obtained by theories) to the target fields (according to the selected wave theory):

$$\Phi =\left(1-{\mathsf{\omega }}_R\right){\Phi }_{target}+{\mathsf{\omega }}_R{\Phi }_{computed}$$

(6)

where $\Phi$ denotes the corrected fields; ${\Phi }_{target}$ and ${\Phi }_{computed}$ are the target and the computed fields, respectively; and ${\mathsf{\omega }}_R={\mathsf{\omega }}_R\left(\sigma \right)\in \left[0, 1\right]$ denotes the weighting function which varies from 1 to 0 as the relaxation zone local coordinate $\sigma$ increases from 0 to 1. The default weighting functions used in waves2Foam is exponential:

$${\mathsf{\omega }}_R=1-\frac{\exp \left({\sigma }^{3.5}\right)-1}{\exp \left(1\right)-1}$$

(7)

2.4. The Post-Processing

The nondimensionalized surface elevations inside the gap and the forces on the floating structures are defined in Equations (8)–(13), which will be used in the post-processing. The free surface elevations at the wave gauges in the wave flume are normalized by the incident wave height $H_i$:

$${\eta }^{*}=\frac{\eta \left(t\right)}{H_i}$$

(8)

where $\eta \left(t\right)$ is the instantaneous free surface elevation relative to the still water level and $H_i$ is the incident wave height. Correspondingly, the normalized surface elevation amplitude is:

$${\eta }_a^{*}=\frac{{\eta }_a}{A_i}$$

(9)

where ${\eta }_a$ is the average value of the elevation amplitudes of several steady-state periods and $A_i$ is the incident wave amplitude.

The wave force $\mathit{f}$ is calculated by integrating over the surface $\mathsf{\Omega }$ of the squares:

$$\mathit{f}={{\displaystyle \int }}_{\mathsf{\Omega }}^{}\left(-\mathbf{n}{p}_{rgh}\left(\mathbf{r},t\right)+\mathbf{n}·\mathit{\tau }\left(\mathbf{r},t\right)\right)d\mathsf{\Omega }$$

(10)

where $\mathbf{n}$ is the unit normal vector pointing into the fluid, $\mathbf{r}$ is the coordinate vector, $p_{rgh}\left(\mathbf{r},t\right)$ is the in excess of the static pressure (see Equation (1)), and $\mathit{\tau }\left(\mathbf{r},t\right)=\mu \nabla \mathit{v}$ is the shear stress tensor.

As a consequence, the instantaneous magnitude $f$ of the wave force $\mathit{f}$ contains both the hydrodynamic components and the hydrostatic components:

$$f=f_{inertial+\mathit{v}iscous}+f_{hydrostatic}$$

(11)

where $f_{hydrostatic}$ is the hydrodynamic component—the hydrostatic component $f_{hydrostatic}$ is the variation of the transient hydrostatic force relative to the initial hydrostatic force in still water.

The normalized parameter $f^{*}$ is introduced to normalize the hydrodynamic force on the squares:

$$f^{*}=\frac{f}{\rho gB{A}_i}$$

(12)

where $B$ is the breadth of the square and $d$ is the draft of the square;

Furthermore, the subscript $x$ was introduced for horizontal force components, e.g., $f_x^{*}$ is the horizontal normalized force and $f_x$ is the horizontal hydrodynamic force.

The normalized horizontal force amplitude is defined as:

$$F_x^{*}=\frac{F_x}{\rho gB{A}_i}$$

(13)

where $F_x$ is the average value of the force amplitude of steady-state periods.

3. Simulation Cases

Twin square bodies are fixed side-by-side in regular incident waves—see SQ1,2 in Figure 2. $B$ is the width of the square, $h$ is the water depth. Several normalized parameters are introduced here: $d^{*}=\frac{d}{B}$; $h^{*}=\frac{h}{d}$; ${\omega }^{*}= \omega /\sqrt{\frac{g}{B}}$, where $\omega$ is the incident wave frequency; $B_g^{*}=\frac{B_g}{B}$, where $B_g$ is the gap width; the Keulegan–Carpenter number (KC number) $\mathrm{KC}=\frac{2\pi A_i}{B_g}$. From hereon, the parameter ${\omega }^{*}$ indicates the nondimensional incident wave frequency of the cases, while the parameter $B_g^{*}$ indicates the nondimensional gap width of the cases.

Two variables were considered in the present study, i.e., the incident wave frequency and the gap width between two squares, which are relevant to the resonance of the water column in the gap. The variations of the gap width can result in different natural frequencies of the water column in the gap. The incident wave frequency can influence the behavior of the gap resonance. Both the characteristics of the hydrodynamic forces and the gap resonance were investigated in the present study. Here $B$, $d^{*},$ and $h^{*}$ are constants which are same as the existing experiments (Saitoh et al., [1]; Tan et al., [16]): $B=0.5$ m; $d^{*}=0.5$; and $h^{*}=2$. The KC numbers of the present cases were small, i.e., less than 3 for all the simulation cases. Therefore, the flows were considered to be inertial dominated [26].

3.1. Numerical Wave Flumes

The reference length for the wave flume design is the resonance wave length ${\lambda }_n$, which is estimated by the first order wave theory:

$$\{\begin{array}{c}{\lambda }_n= 2\pi /k\\ {\omega }_n^2=gk tanh\left(kh\right)\end{array}$$

(14)

where $g$ is the gravitational acceleration, $k$ is the wave number, $h$ is the water depth (see Figure 2); and the resonance frequency ${\omega }_n$ is estimated according to [1]:

$${\omega }_n= \sqrt{\frac{g}{B_{g}B/\left(h-d\right)+d}}$$

(15)

Jacobsen et al. [24] recommended that the length of the relaxation zones should be larger than the incident wave length. Therefore, $L_R$ (Figure 2) was set to $2{\lambda }_n$, $L_W$ was $5{\lambda }_n,$ and the total length $L$ of the numerical wave tank was $9{\lambda }_n$. The $B_g^{*}$ was fixed as 0.1 to investigate the influence of incident wave frequencies, which is the same as the experimental set-up of Saitoh et al. [1] and Tan et al. [4].

3.2. Simulation Cases

(i) The Influences of the Incident Wave Frequencies

In the present simulation cases, different incident wave frequencies were set at the input boundary, which covers the resonance frequency region. However, waves with higher steepness result in larger viscous dissipation in the process of wave propagation [27,28], and the grid resolution needs to vary correspondingly to obtain sufficient numerical accuracy. Therefore, the largest wave frequency was $1.1{\omega }_n$, i.e., the shortest wave length was no smaller than 80% of ${\lambda }_n$. On the other hand, the waves with the smallest frequency pertain to the second-order Stokes wave theory. The incident wave frequencies are listed in

Table 1. Summary of the present numerical simulation cases.

Index	${\mathit{B}}_{\mathit{g}}^{*}$	${\mathit{\omega }}_{\mathbf{n}}^{*}$	KC				${\mathit{\omega }}^{*}$
C1~C9	0.10	1.20	1.51	C1	C2	C3	C4	C5	C6	C7	C8	C9
				0.99	1.07	1.11	1.16	1.17	1.20	1.21	1.24	1.28
C10~C17	0.17	1.09	0.89	C10	C11	C12	C13	C14	C15	C16	C17
				0.89	1.03	1.08	1.09	1.10	1.13	1.14	1.20
C18~C27	0.25	1.00	0.60	C18	C19	C20	C21	C22	C23	C24	C25	C26	C27
				0.77	0.92	0.98	1.00	1.01	1.04	1.06	1.07	1.09	1.11

.

The incident wave amplitude $A_i$ was fixed to $0.012\mathrm{m}$ in accordance with existing experiments (Saitoh et al., [1]; Tan et al., [4]). Besides, four wave gauges (G1–G4) were applied for the free surface elevations and the wave energy dissipation ratio, $R_d=1-R_r{}^2-R_t{}^2$, where $R_r=\frac{A_r}{A_i}$ is the reflection ratio, $A_r$ is the reflected wave amplitude obtained by the two point method [29], $R_t=\frac{A_t}{A_i}$ is the transmission ratio, and $A_t$ is the transmitted wave amplitude measured at G4.

(ii) The Influences of the Gap Width

The $B_g^{*}$ changed from $0.1$ to $0.25$ to include the general range of gap width between two squares in published studies. It should be noted that $B_g^{*}$ is usually small (e.g., $B_g^{*}\approx 0.1$ [1,4,5]). For each $B_g^{*}$, ${\omega }_n$ was calculated by Equation (15) and normalized to ${\omega }_n^{*}$. The selections of incident wave frequencies are shown in Table 1 in detail.

3.3. Boundary Conditions

As shown in Figure 2, the 2-D rectangular computational domain was built around the square structures with inlet, outlet, top, and bottom boundaries. The boundary conditions used in the present study are summarized as follows:

(i) At the inlet boundary, the velocity for the water was given according to the linear wave theory while the air velocity is zero. The normal zero-gradient condition was applied for the pressure.

(ii) At the outlet boundary, the pressure was specified as normal zero-gradient, and the velocities for both water and air were set to zero.

(iii) The top boundary of the computational domain was set as an atmospheric condition, which allowed the air to flow in and out of the domain. The pressure at the top boundary was calculated by $p_T=p_0-0.5u^2$. The velocity u was obtained from the flux at the patch for the inflow and a normal zero-gradient condition for the outflow.

(iv) The no-slip wall boundary condition was used at the bottom and the structures’ surface, where the velocity was zero. The pressure was set as the normal zero-gradient condition.

(v) The “empty” boundary condition was employed at the front and back boundaries due to the present 2-D simulations in OpenFOAM.

4. Mesh and Time-Step Refinement Studies

Mesh and time-step refinement studies were carried out for case C6, which was the case with the smallest gap. The numerical results with different meshes and time-steps were compared to each other in terms of the normalized free surface elevations at G2 (in front of the squares), G3 (in the gap), G4 (behind the squares; see Figure 2), and the hydrodynamic forces on SQ1. As $B_g^{*}$ becomes larger, the resonance amplitude in the gap decreases (Moradi et al., [15]). Therefore, the verified grid settings could be utilized for the other cases with larger $B_g^{*}$. Figure 3 shows the zoom-in view of the mesh around the twin squares of case C6. The vertical grid size over the free surface region is noted as $\Delta {\mathrm{y}}_f$, the horizontal grid size in the free surface region is $\Delta {\mathrm{x}}_f$, and the grid size in the gap and in the vicinity of the square structure boundaries is $\Delta {\mathrm{y}}_{bn}$. The grid size changes gradually at different locations, e.g., between the free surface and the bottom region the grids distribute in the form of hyperbolic cosine function vertically according to the water wave velocity. Details of the grid resolutions are shown in

Table 2. The three different grid resolutions used in the present study.

Mesh Index	$\mathbf{\Delta }{\mathbf{x}}_{\mathit{f}}$	$\mathbf{\Delta }{\mathbf{y}}_{\mathit{f}}$	$\mathbf{\Delta }{\mathbf{y}}_{\mathit{bn}}$	Cell Number
A	${\lambda }_n/80$	$H_i/8$	$B_g/40$	139,270
B	${\lambda }_n/200$	$H_i/20$	$B_g/100$	723,362
C	${\lambda }_n/250$	$H_i/25$	$B_g/125$	1,100,640

. An adaptive time-step scheme was used for simulations with a maximum Courant number of 0.25.

Figure 4 shows the results of the different grid resolutions over five wave periods after the results repeat their cycles. Figure 4a is the normalized surface elevation ${\eta }^{*}$ at G2, G3, and G4 versus the normalized time $t^{*}=\frac{t-t_0}{T_n}$, where $t$ is the simulation time and $t_0$ is the start time of the captured steady-state periods. The relative error of ${\eta }^{*}$ between different grid resolutions is noted as $\mathsf{\epsilon }$, and the numerical model is considered converged when $\mathsf{\epsilon }<5\%$. In Figure 4a, the largest $\mathsf{\epsilon }$ (the $\mathsf{\epsilon }$ at G4) between the numerical results of mesh A and mesh B is 19.79%, and the largest $\mathsf{\epsilon }$ between the numerical results of mesh B and mesh C is 4.95%. Therefore, mesh B is considered to give the sufficient accuracy in terms of the prediction of ${\eta }^{*}$. Figure 4b illustrates $F_{x1}^{*}$ versus $t^{*}$ over the same duration as in Figure 4a, where the largest $\mathsf{\epsilon }$ between the numerical results of mesh A and mesh B is 2.55%, and the largest $\mathsf{\epsilon }$ between the numerical results of mesh B and mesh C is 0.40%. The results show that mesh B is sufficiently fine for predicting both the surface elevations and the horizontal wave forces. Therefore, mesh B was employed for the numerical simulations in the present study.

A time-step refinement study was carried out based on mesh B. A simulation with a maximum Courant number of 0.15 was performed to verify the convergence of the time-step settings. Figure 5 shows the results of C6 using different time-steps over five steady-state periods. The maximum $\mathsf{\epsilon }$ between ${\eta }^{*}$ with two time-steps was 1.91%; the $\mathsf{\epsilon }$ between $F_{x1}^{*}$ with two different time-steps was 2.83%. Therefore, a maximum Courant number of 0.25 was employed for the present numerical simulations.

The verified grid resolution and time-step settings were further applied to the cases with $B_g^{*}=0.17$ and the cases with $B_g^{*}=0.25$. The mesh for cases with $B_g^{*}=0.17$ (cases C10~C17) contained 730,410 cells, while the mesh for the cases with $B_g^{*}=0.25$ (cases C18~C27) contained 899,190 cells; see Table 2 for more information.

5. Results and Discussions

5.1. The Influences of Incident Wave Frequencies

Simulations were performed for cases C1~C9 ($B_g^{*}=0.1$) to investigate the influence of incident wave frequencies. Discussions were carried out in terms of the surface elevation amplitudes in the gap, the viscous dissipation in the gap, the phase-frequency characters of the wave field, and the horizontal wave forces on squares. The subscripts 1 and 2 were introduced to distinguish the forces on the SQ1 and the SQ2, e.g., $f_{x2}$ denotes the horizontal force on SQ2. The present numerical results were compared with the existing experimental and numerical results.

Figure 6 is the comparison of elevation amplitudes inside the gap (${\eta }_a^{*}$) versus the incident wave frequency ${\omega }^{*}$ between the present study and the experiment by Saitoh [1]. The relative difference $\mathsf{\epsilon }$ here was the relative error between the interpolation points of the simulations and the experiment data points. For case C1~C7, $\mathsf{\epsilon }$ was less than 5%, and $\mathsf{\epsilon }$ was slightly larger, i.e., less than 15% for cases C7~C9. The good agreement shows that the present numerical model is able to predict the elevation amplitudes in the gap with reasonable accuracy. As the incident wave frequency ${\omega }^{*}$ became larger, the elevation amplitude inside the gap ${\eta }_a^{*}$ increased rapidly and reached its peak at ${\omega }^{*}=1.195$, before decreasing. This tendency is in accordance with the resonance phenomenon. The mechanism of the resonance can be explained by considering the motion of fluid in the gap as rigid body motion. The motion equation of the fluid vibration in the gap can be expressed as:

$$\left(m+m_a\right)\ddot{\eta }+c\dot{\eta }+k_s\eta =f_{excitation}$$

(16)

where $\eta$ is the fluid displacement or the surface elevation inside the gap, $f_{excitation}$ is the excitation force, $m$ is the mass of fluid inside the gap, $m_a$ is the added mass, $c$ is the damping coefficient, $k_s$ is the stiffness. Based on Equation (16), the natural frequency of the fluid vibration in the gap is:

$${\omega }_n=\sqrt{\frac{k_s}{m+m_a}}$$

(17)

Resonance occurs when the incident wave frequency is equal or nearly equal to ${\omega }_n$. To be specific, the stiffness $k_s$ is the coefficient of $y$ in the buoyancy variation term as the fluid moves, $k_s=\rho g{B}_g$. The first order mass of the fluid inside the gap is the fluid mass inside the gap of still-water state, $m=\rho B_{g}d$. Equation (17) can further be simplified according to these formulas:

$${\omega }_n=\sqrt{\frac{g}{\frac{m_a}{\rho B_g}+d}}$$

(18)

Equations (16)–(18) are the basis of follow-up discussions in the present study. It should be noted that these equations are not limited to the rigid body assumption. It is a general motion equation of the vertical flow motion in the gap. The coefficients $m_a$ and $c$ change significantly between different theoretical models. For example, Equation (15) (derived by Saitoh, [1] using the energy method) is equivalent to take $m_a$ as $\frac{\rho B{B}_g^2}{h-d}$. In Tan et al. [4], $m_a$ was $\frac{\rho B^{*}{B}_g^2}{4\left(h-d\right)}$, where $B^{*}$ is the artificial coefficient.

The viscous dissipation influences the flow field around the squares pronouncedly. In analytical and potential methods, damping is usually obtained from CFD simulations or experimental data. It is generally agreed that the vorticity near the entrance of the gap is the main source of the viscous dissipation. Figure 7 shows the variation of the vorticity contour in the water phase for case C6 (${\omega }^{*}$ equals to ${\omega }_n^{*}$) at $t^{*}=0,0.17, 0.33,0.50,0.67,1$. When $t^{*}=0$, the water level in the gap was at the lowest position, and the average flow velocity in the gap was nearly zero. As the wave gradually travels past the square bodies at $t^{*}=0~0.25$, part of the waves reflected back and overlapped with the incident waves, part of the wave energy was absorbed by the movement of the elevation in the gap, and the rest was the transmitted waves. Due to the large velocity gradient at the corners of the SQ1 and SQ2, two vortices began to grow symmetrically in the gap to dissipate the energy. From $t^{*}=0.25~0.5$, the decreasing flux in the gap resulted in the hydrodynamic pressure gradient, which was at the same direction as the flow direction in the gap, and pushed the vortex back to the gap entrance. From $t^{*}=0.5~1$, the vortices gradually spread to the bottom region below the square structures as the water level in the gap returned back to the lowest position.

Figure 8 and Figure 9 present the variations of the vorticity contour of case C2 (${\omega }^{*}$ smaller than ${\omega }_n^{*}$) and C9 (${\omega }^{*}$ larger than ${\omega }_n^{*}$). In these two cases, the flow in the gap produced smaller vortices as compared to case C6. Furthermore, Figure 10 shows the square of transmission ratio $R_t{}^2$, the square of the reflection ratio $R_r{}^2$ and the dissipation ratio $R_d$ versus the incident wave frequency ${\omega }^{*}$. The present simulation results were compared with the experimental data from Tan et al. [16]. The linear interpolation method was applied for the numerical results according to the sample points from the experimental data. The root mean square errors (RMSE) of the interpolated numerical results relative to the experimental data were calculated. For the $m$th wave frequency, $R_{1m}$ and $R_{2m}$ represent the interpolated numerical results and the experiment data, and the RMSE was calculated by: $\epsilon =\sqrt{{{\displaystyle \sum }}_{m=1}^M[{\left(R_{1m}-R_{2m}\right)}^2/M]}$, where $M$ is the total number of frequencies of the experiments. In Figure 10, the RMSE of $R_r^2$, $R_t^2,$ and $R_d$ are 0.061, 0.04, and 0.053, respectively, which are relatively small values. There was a slight difference of peak frequency in the present simulation and experiments—i.e., $R_t^2$ was approximately 1.17 in the present numerical results, and in the experiments, the peak frequency for $R_t^2$ was 1.15. Though the difference was only around 1.7%, it was able to result in a significant difference when the resonance interval was narrow and $R_t^2$ was nearly zero out of the resonance interval. Another reason for the discrepancies was the uncertainties of the experiments. For example, when the incident wave frequency ${\omega }^{*}$ was approximately 1.02, the $R_d$ in the experiment was negative, Tan et al. (2016) believed that this unphysical phenomenon was caused by the measurement uncertainties. In general, the present numerical results are in good agreement with the experiments.

When gap resonance happens, the wave reflection ratio decreases a lot, which is mainly caused by the increased dissipation ratio. Assuming that the energy loss in the gap is equal to the energy loss of the whole wave field in each period, the value of damping coefficient $c$ (Equation (16)) can be estimated as:

$$\frac{\rho g\omega }{4k}\left[1+\frac{2kh}{\sinh \left(2kh\right)}\right]A_i^2{R}_d=\frac{1}{T}\underset{0}{\overset{T}{{\displaystyle \int }}}c{\dot{\eta }}^{2}dt \stackrel{}{\Rightarrow }\mathrm{c}=\frac{\rho g}{2k\omega }\left[1+\frac{2kh}{\sin \mathrm{h}\left(2kh\right)}\right]\frac{R_d}{{\left({\eta }_a^{*}\right)}^2}$$

(19)

The damping coefficient $\mathrm{c}$ versus ${\omega }^{*}$ is illustrated in Figure 11. As ${\omega }^{*}$ approaches ${\omega }_n^{*}$, the damping coefficient $c$ gradually converges to a constant value. Tan et al. [16] suggested that the linearized damping coefficient is proportional to the excitation frequency $\mathsf{\omega }$: $c=\epsilon \omega M$, where $\epsilon$ is a non-dimensional damping coefficient obtained from the CFD results of the resonance case. In order to make comparisons, the damping coefficient $c$ in present simulations was reformulated as:

$$c=c_{\omega }\omega$$

(20)

where $c_{\omega }$ is the damping coefficient $c$ of the resonance case divided by the resonance frequency ${\omega }_n$. In cases with the gap width $B_g^{*}=0.1$ (cases C1~C9) the value of $c_{\omega }=\frac{5.858}{{\omega }_n}=1.107$, which agrees well with that of the cases with the same gap width $B_g^{*}$ of Tan et al. [16] ($c_{\omega }=1.010$). It should be noted that their results were obtained from CFD calculations including turbulence modeling. This demonstrates that the present laminar model has enough accuracy for predicting the viscous dissipation in the gap, as mentioned in Section 3.

To simplify the problem, in later discussions, it was assumed that the only influence of viscosity is to determine the surface elevation inside the gap, while the flow outside the gap is irrotational. Firstly, as shown in Figure 7, Figure 8 and Figure 9, the vorticity mostly gathered in the bottom area of the gap and had little influence on the flow tendency of the entire flow field. Secondly, several studies [9,10,11,17,18] indicated that the linear potential flow theory based results are reasonable after taking into account the damping inside the gap.

From Figure 7, Figure 8 and Figure 9, the phase retardation of the elevation inside the gap relative to the wave profile in front of SQ1 varied visibly. This interesting character was further investigated. The phase of the elevation inside the gap (${\varphi }_1$) was evaluated from the elevation time history inside the gap ($\eta \left(t\right)$):

$${\varphi }_1=\omega t_{crest}$$

(21)

where $t_{crest}$ is the time corresponding to the peak of $\eta \left(t\right)$ and $\omega$ is the incident wave frequency.

In the present simulations the incident wave phase in front of SQ1 (${\varphi }_2$) was:

$${\varphi }_2=k{x}_{front}+2\pi n_p$$

(22)

where $k$ is the wave number, $x_{front}$ is the x coordinate in front of the SQ1, and $n_p$ is the number of the wave period.

The wave phase in front of the SQ1 may deviate slightly from the result of Equation (22) as a consequence of the superposition between the incident wave and the reflected wave. However, the effect is negligible. It was assumed that the incident wave completely reflected back at the left side of SQ1. Based on the wave superposition theory, this generated a standing wave with the following profile function:

$$\begin{array}{cc}\hfill {\eta }_{sw}& =A_i\cos \left(\omega t-kx\right)+A_i\cos \left(\omega t+kx-2x_{front}\right)\hfill \\ & =2{\mathrm{A}}_{\mathrm{i}}\cos \left(\omega t-k{x}_{front}\right)\cos \left[k\left(x-x_{front}\right)\right]\hfill \end{array}$$

(23)

Figure 12 is the comparison between the surface elevations in front of SQ1 in the present simulations and of the presumed standing waves in front of SQ1. The time history curves of the surface elevation in front of SQ1 were almost coincident with the elevation time histories of the presumed standing wave in front of SQ1, which were calculated based on Equation (22). This indicates that most parts of the incident wave reflected back because of the strong shielding effect of the twin square structures. By comparing Equations (22) and (23), it could readily be drawn that the wave profile in front of SQ1 was almost in-phase with the incident wave profile at the left side surface of SQ1.

Therefore, the phase retardation of the elevation inside the gap relative to the wave phase in front of SQ1 (Hereinafter referred to as the phase leg $\varphi$) was:

$$\varphi ={\varphi }_1-{\varphi }_2$$

(24)

In order to clarify the character of the phase leg $\varphi$, a further investigation on the excitation force ($f_{excitation}$) of the liquid inside the gap was made. As a matter of investigating the excitation force, the liquid inside the gap was considered stationary. The twin squares and the liquid inside the gap could be regarded as an entire rectangular obstacle with a width of $B_t=2B+B_g$. In the control volume below the obstacle (see region CV in Figure 13), the velocity potential $\psi$ followed the impermeable boundary condition at the boundaries $ad$ and $bc$; it could therefore be written as (Cong et al., [6]):

$$\psi =A_{0}x e^{-j\omega t}+{\displaystyle \sum }_{n=1}^{\infty }\left[A_n{e}^{{\mu }_{n}x}+B_n{e}^{{\mu }_n\left(B_t-x\right)}\right]\cos {\mu }_n\left(y+h\right)e^{-j\omega t}$$

(25)

where ${\mu }_n=\frac{n\pi }{h-d}$ is the eigenvalue ($n=1,2,\dots$); $A_0$, $A_n$, and $B_n$ are complex constants; and $j$ is the imaginary unit.

Equation (25) indicates that the flow in the region CV was composed of two parts, i.e., the first part is the horizontal uniform oscillation with the velocity potential of $A_{0}x e^{-j\omega t}$, and the second part is the disturbance of velocity with the velocity potential of ${{\displaystyle \sum }}_{n=1}^{\infty }\left[A_n{e}^{{\mu }_{n}x}+B_n{e}^{{\mu }_n\left(B_t-x\right)}\right]\cos {\mu }_n\left(y+h\right)e^{-j\omega t}$. In fact, the free surface wave velocity on the left side of the boundary $ab$ decays with the water depth. When the relative draft $d/\lambda$ is sufficiently large, the disturbance of velocity at the boundary $ab$ becomes small quantity. The fluid flow in the control volume is mainly the horizontal uniform oscillation. The pressure distribution in the uniform flow field is the same everywhere. This indicates that the pressure inside the whole control volume CV was in-phase with the pressure at the boundary $ab$, which is the pressure of the wave field on the left side of boundary $ab$. When the wave profile in front of SQ1 rose up, the pressure inside the control volume CV increased and vice versa. Therefore, the excitation force ($f_{excitation}$) on the liquid inside the gap synchronized with the wave elevation in front of SQ1. The phase leg $\varphi$ was the same as the phase retardation of the response relative to the excitation of the vibration system.

Liu [30] carried out a series of experiments to investigate the phase leg $\varphi$ and concluded that the change of the phase leg $\varphi$ versus the incident wave frequency satisfied the phase-frequency character of the vibration system under harmonic excitation. For verification of the present simulation results, the phase leg $\varphi$ was theoretically estimated based on Equations (16)–(18) and the structural dynamics theory (relevant theories could be found in textbooks like Humar, [31]):

$$\varphi =\arctan \frac{2\zeta \gamma }{1-{\gamma }^2}$$

(26a)

$$\zeta =\frac{c}{2\left(\rho B_{g}d+m_a\right){\omega }_n}$$

(26b)

$$\gamma =\frac{\omega }{{\omega }_n}$$

(26c)

At present, the parameters $\left(\rho B_{g}d+m_a\right)$, ${\omega }_n$ in Equation (26b,c) are known constants. The damping coefficient $c$ is proportional to the incident wave frequency $\omega$ (Equation (20)). Therefore, the phase leg $\varphi$ versus the incident wave frequency $\omega$ could be plotted immediately based on Equation (26a–c). On the other hand, the phase leg $\varphi$ in the present simulations was evaluated based on Equation (24). Figure 14 shows the comparison of the phase leg $\varphi$ of the theoretical estimation (based on Equation (26a–c) and the simulation results (based on Equation (24)). The theoretical estimation curve and the simulation results are in good agreement with each other. A brief summary could be made from the above results. Firstly, in the present simulations, the excitation force ($f_{excitation}$) on the liquid inside the gap synchronized with the wave elevation in front of SQ1, which was almost in-phase with the incident wave profile in front of SQ1. Secondly, the phase leg $\varphi$ could be estimated by referring to the phase-character of the vibration system, i.e., Equation (26a–c).

Figure 15 shows the horizontal force amplitudes on the twin squares. The results of the present simulations were compared with the results of Lu et al. [18] using a viscous flow solver based on the three-step Taylor–Galerkin finite element method. The normalized horizontal force amplitude $F_x^{*}$ in resonance interval of Lu et al. [18] was smaller than that of the present simulation, whereas the variation tendencies of force amplitudes were similar in both the present solver and Lu et al. [18]’s solvers. It was found that the horizontal force amplitude on SQ2 ($F_{x2}^{*}$) has the same phase with the elevation amplitude in the gap (${\eta }_a^{*}$). However, the peak frequency of the horizontal force amplitude on SQ1 ($F_{x1}^{*}$) was considerably larger than that of SQ2 ($F_{x2}^{*}$).

To explain the generation mechanics of the horizontal forces, first order theoretical estimations of the horizontal force amplitudes were conducted. The first order free surface elevation in the gap was harmonic and can be expressed as:

$$\eta ={\eta }_a\cos \left(\omega t-\varphi \right)$$

(27)

Assuming that the flow in the gap is uniform, the horizontal force produced by the vibration in the gap can be estimated using the Lagrange integral from the free surface to the target location:

$$f_{xg}^{\left(1\right)}={{\displaystyle \int }}_{-d}^0(-\rho \frac{\partial \psi }{\partial t}) dy+{{\displaystyle \int }}_{-d}^0\rho g\eta dy=(\rho gd{\eta }_a-\frac{{\omega }^2{d}^2}{2}{\eta }_a)\cos \left(\omega t-\varphi \right)$$

(28)

where $\psi =\dot{\eta }y$ is the velocity potential of the uniform flow inside the gap.

Using the standing wave theory, the first order wave force on the left side of SQ1 could be estimated by:

$$f_{xw}^{\left(1\right)}={{\displaystyle \int }}_{-\mathrm{d}}^0\rho g{\xi }_a\frac{\cos \mathrm{h}k\left(y+h\right)}{\cos \mathrm{h}\left(kh\right)}\cos \left(\omega t\right) dy=\rho g{\xi }_a\frac{\sinh kh-\sinh k\left(h-d\right)}{k\cos \mathrm{h}\left(kh\right)}\cos \left(\omega t\right)$$

(29)

where ${\xi }_a$ is the equivalent standing wave amplitude in front of SQ1. According to the energy conservation law, ${\xi }_a$ is limited between two times the reflection wave amplitude and is two times the incident wave amplitude. Therefore, it is estimated as $\frac{2A_r+2A_i}{2}$.

The horizontal force on SQ1 includes the effects of both $f_{xg}^{\left(1\right)}$ and $f_{xw}^{\left(1\right)}$:

$$f_{x1}^{\left(1\right)}=f_{xw}^{\left(1\right)}-f_{xg}^{\left(1\right)}=\sqrt{{\left(F_{xw}-F_{xg}\cos \varphi \right)}^2+{\left(F_{xg}\sin \varphi \right)}^2}\cos \left(\omega t-{\varphi }_1\right)$$

(30)

where ${\varphi }_1$ is the phase difference between $f_{x1}^{\left(1\right)}$ and $f_{xw}^{\left(1\right)}$.

The horizontal force on SQ2 is mainly induced by the vibration in the gap:

$$f_{x2}^{\left(1\right)}=f_{xg}^{\left(1\right)}$$

(31)

Figure 16 shows the comparison between the horizontal force amplitudes that were estimated by Equations (30) and (31) and in the present simulation results. $F_{x1}^{*}$ and $F_{x2}^{*}$ ($F^{*}$ and $F$ respectively denote the normalized amplitude and amplitude of corresponded force term (see Section 2.4) estimated by Equations (30) and (31) agree well with the results from the present CFD simulations.

As a summary, firstly, the phase difference between $f_{xw}$ and $f_{xg}$ influenced $F_{x1}$ significantly. To be specific, $f_{xw}$ and $f_{xg}$ counteracted each other when ${\omega }^{*}$ was smaller than ${\omega }_n^{*}$, resulting in a smaller $F_{x1}$ than $F_{xg}$; $f_{xw}$ and $f_{xg}$ mutually promoted when ${\omega }^{*}$ was larger than ${\omega }_n^{*}$; the peak frequency of $F_{x1}$ therefore fell behind the resonance frequency. Secondly, due to the shielding effect of the up-stream structures, the wave forces from the transmitted waves were negligible compared with the hydrodynamic forces in the gap $f_{xg}$. Therefore, the force amplitude on SQ2 $F_{x2}$ was almost equal to the hydrodynamic forces in the gap $F_{xg}$. Aforementioned discussions reveal the mechanism of the horizontal forces on the square structures.

5.2. The Influences of Gap Width

Figure 17 shows the elevation amplitude inside the gap ${\eta }_a^{*}$ versus the incident wave frequency ${\omega }^{*}$ for the cases with different gap width $B_g^{*}$. The resonance frequency and the resonance amplitude became smaller as the gap width increased. A similar relationship has also been reported in the experiments conducted by Saitoh et al. [1] and the numerical simulations of Lu et al. [17] and Moradi et al. [15]. As the gap width $B_g^{*}$ increased, the resonance frequency ${\omega }_n^{*}$ in the present simulations gradually became smaller than the theoretical estimations based on Equation (17) (Saitoh et al., [1]), which is due to the stronger transverse flow in the gap. However, Equations (18)–(20) could still be applied to describe the mean vertical flow in the gap, as described in Section 5.1. For the cases with larger gap width $B_g^{*}$, the value of the resonance frequency ${\omega }_n^{*}$ could be estimated based on the characteristics of the phase leg $\varphi$ of the elevation inside the gap relative to the wave phase in front of SQ1. According to Equation (26a), the phase leg $\varphi$ approached $\mathsf{\pi }/2$ when the incident wave frequency approached to the gap resonance frequency ${\omega }_n$. For example, the phase leg $\varphi$ of cases C23 and C24 were most close to $\mathsf{\pi }/2$. among the cases with the gap width $B_g^{*}=0.25$. Note the value of the phase leg $\varphi$ of case C23 as ${\varphi }_{\mathsf{\pi }/2-}$, which was slightly smaller than $\mathsf{\pi }/2$; the phase leg $\varphi$ of case C24 as ${\varphi }_{\mathsf{\pi }/2+}$, which was slightly larger than $\mathsf{\pi }/2$. The resonance frequency ${\omega }_n^{*}$ of the cases with the gap width $B_g^{*}=0.25$ could be estimated by interpolation:

$${\omega }_n^{*}={\omega }_{C23}^{*}+\frac{\mathsf{\pi }/2-{\varphi }_{\mathsf{\pi }/2-}}{{\varphi }_{\mathsf{\pi }/2+}-{\varphi }_{\mathsf{\pi }/2-}}\left({\omega }_{C24}^{*}-{\omega }_{C23}^{*}\right)$$

(32)

where ${\omega }_{C23}^{*}$, ${\omega }_{C24}^{*}$ are, respectively, the incident wave frequencies of case C23 and case C24.

The resonance frequency ${\omega }_n^{*}$ of the cases with different gap width $B_g^{*}$ and the dimensionless added mass $m_a/\left[\rho B{B}_g^2/\left(h-d\right)\right]$ that were deduced from the resonance frequency ${\omega }_n^{*}$ based on Equation (18) are illustrated in Figure 18. The dissipation ratio $R_d$ and the damping coefficient $\mathrm{c}$ of the cases with different $B_g^{*}$ were calculated and shown in Figure 19, and the phase leg $\varphi$ is illustrated in Figure 20. The resonance frequency ${\omega }_n^{*}$ and the corresponding dimensionless added mass $m_a/\left[\rho B{B}_g^2/\left(h-d\right)\right]$ decreased with the increase of $B_g^{*}$. Meanwhile, the damping coefficient $\mathrm{c}$ increased with the increase of the gap width $B_g^{*}$, which indicates that the viscous dissipation in the gap became stronger when the gap width $B_g^{*}$ became larger. However, in Figure 20, the phase leg $\varphi$ estimated based on Equation (26a–c) shows a big difference from the present simulation results when $B_g^{*}$ was large. To be specific, as the incident wave frequency increasing, the phase leg $\varphi$ estimated based on Equation (26a–c) increased faster than the data points of the present simulations, which indicates that the damping ratio $\zeta$ and the damping coefficient $c$ that are used in Equation (26a–c) were smaller than the corresponding values that the simulations indicate. The reasonable explanation is that when the gap width $B_g^{*}$ was relatively large, the reduction of the piston flow in the gap was not only caused by the viscous dissipation—it was also caused by the sloshing mode flow in the gap. For example, Figure 21 shows the comparison of the flow structures in the gap between cases C6 ($B_g^{*}=0.1$) and C21 ($B_g^{*}=0.25$). For case C6, the flow in the gap was mostly in the vertical direction. However, for case C21, the sloshing mode of the flow became much stronger and extracted a large part of the mechanical energy from the vertical flow in the gap. The damping coefficient $\mathrm{c}$ in Equations (16) and (26) could not simply be estimated according to the dissipation ratio $R_d$ (Equation (19)). The coupling between the piston and the sloshing types of flow motion in the gap was very complex, which could be captured well using the present numerical method.

The horizontal force amplitude versus the incident wave frequency of the cases with $B_g^{*}=0.17$ and $B_g^{*}=0.25$ are presented in Figure 22 and Figure 23, respectively. The variation tendencies of the horizontal force amplitudes were well predicted by the theoretical estimations based on Equations (30) and (31). As the gap width $B_g^{*}$ increased, the decreasing resonance amplitude in the gap did not necessarily reduce the peak value of the horizontal force amplitude $F_x^{*}$. For example, the maximum $F_{x1}^{*}$ increased from 1.85 to 1.92 as $B_g^{*}$ increased from 0.17 to 0.25. The maximum $F_{x2}^{*}$ increased from 1.63 to 1.68 as $B_g^{*}$ increased from 0.1 to 0.17. The above results are reasonable. In Equation (28), the horizontal force amplitude induced by the fluid motion inside the gap ($F_{xg}$) includes two parts, namely the hydrostatic part $\rho gd{\eta }_a$ and the hydrodynamic part $-\frac{{\omega }^2{d}^2}{2}{\eta }_a$. As the resonance amplitude decreases, the hydrostatic part decreases proportionally. However, the hydrodynamic part may increase due to the rapid decrease of the resonance frequency ${\omega }_n$. As a consequence, the horizontal force amplitude on the SQ2 ($F_{x2}$) was likely to increase. On the other hand, due to the effect of the wave field at the left side of SQ1, the probability for a horizontal force amplitude on the SQ1 $F_{x1}$ increased even greater. The aforementioned discussions indicate that the gap width influences the hydrodynamic forces from various aspects. Enlarging gap width can reduce resonance amplitude; however, it does not directly result in the decrease of the maximum horizontal force amplitude on the floating structures.

6. Conclusions

Two-dimensional computational fluid dynamics (CFD) simulations were carried out to investigate the gap resonance phenomenon that occurs when free-surface waves propagate past twin alongside placed squares. The laminar model was used to solve the governing equations, and the volume of fluid method was used for free surface capturing. Both mesh and time-step convergence studies were carried out to ensure sufficient numerical accuracy. The influences of the incident wave frequency, as well as the gap width, were investigated. Detailed discussions were carried out based on the motion equation of the fluid inside the gap in terms of the damping coefficient, the phase retardation of the elevation inside the gap relative to the wave phase in front of the weather side structure, and the horizontal force amplitudes on the square structures. The following conclusions can be drawn from the discussions:

• For the cases of free surface waves past through twin alongside placed 2-D identical rectangular squares, the laminar model with the volume of fluid method was able to predict the gap resonance amplitudes, the viscous dissipation effects, and wave forces on the structures accurately when the KC number was small (e.g., $KC=\frac{2\pi A_i}{B_g}<3$, where $A_i$ is the incident wave amplitude and $B_g$ is the gap width).
• As the incident wave frequency increased, the surface elevation amplitude inside the gap first increased and then decreased, a tendency which is in accordance with resonance phenomena. The horizontal force amplitude on the lee side square structure changed in-phase with the elevation amplitude inside the gap, while the horizontal force amplitude on the weather side structure reached the peak value at a larger frequency than the gap resonance frequency.
• In present simulations, the phase retardation of the elevation inside the gap relative to the wave phase in front of the weather side structure satisfied the phase-frequency characteristics of the vibration system and determined the change tendency of the horizontal force amplitude on the weather side structure. To be specific, the horizontal wave force on the left side of the weather side structure and the horizontal force on the right side of the structure (induced by the fluid inside the gap) counteracted each other when the wave frequency was smaller than the gap resonance frequency. This resulted in a smaller horizontal force amplitude on the structure than the horizontal fluid force inside the gap. The forces on each side of the structure mutually promoted when the incident wave frequency was larger than gap resonance frequency, making the peak frequency of the horizontal force on the structure fall behind the gap resonance frequency.
• The increase of the gap width resulted in the increase of the added mass and reduced the resonance frequency. Moreover, as the gap width increased, the decrease of resonance amplitude in the gap did not necessarily reduce the peak value of the horizontal forces on the squares.
• As the gap width increased, the viscous dissipation and the sloshing mode flow inside the gap both became stronger.