In recent years, with the fast development of ocean engineering, the hydrodynamic research about the narrow gap between two floating structures attracts more and more attentions because of the possible violent resonance phenomenon of the fluid within the narrow gap. In reality, the fluid resonance phenomenon may exist in the narrow gap between two offshore structures subject to wave loading, such as the ship-by-ship operations, a LNG carrier alongside a terminal and the moonpools of the FPSO (Floating production storage and offloading facilities).

Many existing studies are focused on two-dimensional configurations by assuming a infinite long gap. Newman and Sclavounos (1988) modeled wave interaction with two close rectangular barges by a panel method and reported unusually high wave elevations within the narrow gap between the barges as well as large hydrodynamic forces. Miao et al. (2000, 2001) studied the influence of the gap between two stationary rectangular caissons on wave forces using the boundary element method. It was found that the resonant wave forces on each caisson could be ten times of the forces on an isolated caisson. Miao et al.(2000, 2001) also conducted theoretical analyses to demonstrate that the fluid resonance takes place when kL=nπ (n=1, 2, 3, …, ∞), where k is the incident wave number and L is the width of caisson. Iwata et al. (2007) investigated the hydrodynamic resonance of three identical rectangular caissons with two narrow gaps by laboratory tests. It was found that the fluid resonance characteristics of three identical bodies with two narrow gaps were different from those observed in the case of twin bodies with a single narrow gap. The fluid resonance in narrow gaps may take place at two distinct frequencies for three bodies with two gaps, while the resonance was observed to occur only at one frequency for two structures with one gap. Kristiansen and Faltinsen (2010) investigated the gap resonance between a ship and a vertical wall in a two-dimensional setting. The ship was allowed to oscillate in three degrees of freedom, i.e., Sway, heave and roll. The investigation was carried out by means of model tests as well as by time-domain simulations based on a Boundary Element Method and the mixed Eulerian-Lagrangians cheme. Flow separation was visualized in the simulation by means of an inviscid vortex tracking method. The numerical results and model test results agree well with each other. Lu et al. (2010) studied the fluid resonance of the three identical bodies with two narrow gaps by a numerical method based on the Navier-Stokes equations to account for the viscous dissipation of the energy. In the numerical model, the free surface was captured using the CLEAR-Volume of Fluid (CLEAR-VOF) method and the results were in good agreement with the experiental data in Iwata et al. (2007). To improve the simulation accuracy (mainly the response amplitude in the gap) and maintain the efficiency of the potential flow model, Lu et al. (2011) introduced an artificial damping term into the potential flow model to account for the viscous dissipation effect. In the numerical model by Kristiansen and Faltinsen (2012), the computation domain of Navier-Stokes solver was fully submerged in the fluid so that it could capture viscous effects and the free surface was simulated using the potential flowtheory. Ning et al. (2016) numerically investigated double solitary waves interaction with a two-boxes system with a narrow gap and found the spoon-shaped variation of the second wave run-up with the time interval between two solitary waves.

Two-dimensional models, however, are limited in capturing real three-dimensional characteristics of the fluid dynamics in the gap. Some recent experimental and numerical simulations have been focused on three-dimensional structures. Molin et al. (2009) carried out laboratory experiments about narrow gap resonance of two fixed rectangular barges inwaves. Wave elevations along the centerline between the two barges were recorded, and the experimental data have been widely used to validate numerical models. Sun et al. (2010) simulated wave interaction with two barges with a narrow gap using the second-order diffraction theory based on the quadratic boundary element method, and examined the behavior of the wave motion in the gap and the influence of spacing between the two barges on the gap resonance. Hong et al. (2013) investigated the hydrodynamic forces and wave elevations between two floating side-by-side bodies using two methods, i.e., a nine-node discontinuous higher order boundary element method based on the conventional two-body formulation and a constant boundary element method based on the boundary matching formulation. The results from the two methods were found to agree with each other. Yan et al. (2009) used the quasi arbitrary Lagrangian-Eulerian finite element method to investigate the relationship of resonance frequency and the gap width for nonlinear hydrodynamic interaction between two floating structures. Feng and Bai (2015) applied fully nonlinear potential flow model to investigate the wave resonance in the gap between two side-by-side barges underbeam sea conditions. They focused their study on the nonlinear behaviors and the stiff/soft spring behaviors of the gap resonances. Watai et al. (2015) carried out experiments and numerical simulations using potential flow theory of two different bodies of canonical geometries (a barge and a geosim) with two different side-by-side configurations. The comparison between measurements and numerical results demonstrated that an external damping factor is required to improve the accuracy of the potential flow model.

Previous numerical studies of gap resonance of two three-dimensional barges were limited to the validation of the numerical models and were mainly focused on beam sea conditions. The influence of the wave direction on the wave elevation in the gap has not been fully understood. The aim of this paper is to study the influence of incident wave direction on the wave resonance in a gap between two fixed barges using a second-order time-domain potential flow model. In the numerical model, a circular computational domain with the structures being placed at its center is employed. A damping zone is placed on the free surface at the external boundary of the domain to absorb the outgoing scattering wave. First, the model is validated by comparing the results with experimental data and other numerical results. Then, the influence of the incident wave direction on the wave resonance is investigated over a range of incident wave frequencies.

A circular fluid domain is used to study the wave-structure interaction problem as shown in Fig. 1. In the numerical simulation, an annulus damping zone is used to disspate the scattered wave energy, in order to satisfy the radiation boundary condition. The origin of a coordinate system Oxyz is placed at the center of the gap on the still water surface, with z-axis pointing upwards and x-direction parallel to the gap. The direction angle of the incident wave (denoted by β) is measured from the positive x-direction in counter clock direction and the damping zone has a width of about one wavelength.

In the second-order time domain model, the fluid is assumed to be incompressible, inviscid and flow is irrotational within the fluid domain, so the flow velocity potential satisfies the Laplace equation and the boundary conditions. The velocity potential is decomposed into the incident potential ϕ_i and the scattered potential ϕ_s. Thus, the scattered potential satisfies the Laplace equation in the domain (Bai and Teng, 2013):

and is subjected to the body boundary condition:

where superscript k (= 1 & 2) denotes the order of approximations in satisfying the nonlinear boundary conditions. First and second-order kinematic and dynamic free surface conditions are given as the following:

where ${\,f' _k}$ and ${\,f'' _k}$ are the forcing terms on the still free surface, and are defined as follows:

We apply a Rankine source and its image about the seabed as Green’s function:

where $\text{x}$₀=(x₀, y₀, z₀) and $\text{x}$=(x, y, z) are the source point and the field point, respectively, and

By applying the second theorem of Green to the scattered potential and Green function, the above boundary value problem can be converted to the following boundary integral equation:

where α is the solid angle coefficient, which depends on the surface shape of the structure. Because both the Rankine source and its image about the seabed are used as the Green function and the scattered waves is absorbed by the damping layer, the integration on the seabed and the outer boundary are not needed. The integration is only needed on the boundary S containing the wetted body S_b and the static water surface S_f. A higher-order boundary element method is used to solve the boundary integration equation and 8-node elements are used to discretize the boundaries.

In the calculation, the normal derivative of the first-order scattered potential on the body surface and the first-order scattered potential on the free surface are known, so the first-order scattered potential on the body surface and the normal derivative of the first-order scattered potential on the free surface can be computed according to the integral Eq. (11). After obtaining the first-order potential and the first-order wave surface elevation, we will construct the second-order terms about the free surface and apply the integral equation to compute the second-order scattered potential of the body surface and the normal derivative of second-order scattered potential of the free surface.

By shifting the terms with unknown to the left-hand side of Eq. (11), we have

for the source point on the body surface, and

for the source point on the free surface.

In the above two equations, the left-hand sides of the integral equations contain the unknowns that are the velocity potential on the body surface and the normal derivative of the velocity potential on the free surface.

On the free surface, the fourth-order Adams-Bashforth predict-correct method is used in time integration. The detailed analysis process is described below.

Step 1: At each time step, the time is t. the wave elevation ${\eta _{\rm s}}(t)$ and velocity potential ${\phi _{\rm s}}(t)$ are known, and the kinetic and dynamic free surface boundary conditions can be written in a general form of function as follows:

Step 2: Δt is used to stand for the time step. The fourth-order Adams-Bashforth method is used to predict the wave elevation and velocity potential on the free surface at the next time step as follows:

Step 3: The predicted free surface condition at t+Δt is substituted into the integral Eqs (12) and (13), we can obtain the normal derivative of the potential of the free surface. Then the fourth-order Adams-Moulton method is used to correct the wave elevation and velocity potential on the free surface as follows:

When $\Delta \phi _{\rm s}^{(k)}$ and $\Delta \eta _{\rm s}^{(k)}$ are small enough, we go back to Step 2 for a new time step. In order to avoid the effect of intitial conditions, a ramp function has also been introduced in the body boundary conditions, which is given as follows:

where T is the wave period.

The above boundary value problem for solving the scattered wave is equivalent to that for the total potential flow if the incident potential satisfies the Laplace equation and all the boundary conditions. For the single directional wave problem, the first and second order incident velocity potentials and wave elevation can be expressed as

where η is the incident wave elevation, k is the wave number, ω is the angular frequency and A is the wave amplitude.

Two identical stationary rectangular barges with square bilges in the experiments by Molin et al. (2009) are simulated. The configuration of the side-by-side barges at model scale is same as that used by Molin et al. (2009). The parameters are: barge length 2.47 m, width 0.6 m, draft 0.18 m and gap width 0.12 m. The water depth is set as 3 m (same as those in the tests). The incident wave angle considered here is 90°, i.e., the beam sea.

In order to demonstrate the mesh convergence of the simulations, three meshes with different densities are used at the two angular frequencies of 5.75 and 6.80 rad/s, respectively. The boundary element sizes on the free surface in the gap and near the barges are much smaller than those far away from the barges, and the boundary element sizes in the gap or near the barges in case of “mesh A”, “mesh B” and “mesh C” are approximately 1/16, 1/20, 1/24 of the shorter incident wavelength, respectively. Utilysing symmetry of the configuration, only a quadrant of the mesh grid is used. Figure 2 shows “Mesh B” with 1 016 points and 350 elements on the barges, 2 512 points and 785 elements on the free water surface per quadrant, respectively.

Figure 3 shows the non-dimensional time history of the wave elevation at the gap centre calculated from the three meshes. Fully developed results are achieved after t/T=30 for all the meshes and the results from different meshes are very close to each other. To quantify the mesh dependency study, Table 1 shows the nondimensional wave crest η_crest/A for ω=5.75 and 6.80 rad/s, where η_crest is the maximum value of the wave elevation. The differences between Meshes B and C for ω=5.75 and 6.80 rad/s are 0.18% and 0.16%, respectively. Therefore Mesh B is chosen for the subsequent simulations as a compromise between computational accuracy and efficiency.

To investigate the resonance frequencies of the wave surface elevation, we simulate the barge system under the beam sea conditions. As stated in Section 2, the second-order Stokes waves are specified as the incoming waves. We used the wave steepness kA=0.003 4 (same as the one in Feng and Bai (2015)) and frequencies in range between 5.0 and 9.0 rad/s. The comparisons of simulated wave crest elevation at two positions of x=0.0 m and 0.3 m with the solutions from the linear theory, fully nonlinear theory and experiment data are shown in Figs 4 and 5, respectively. The overall agreement is good and the wave crest elevation from each method is characterized by three peaks within the simulated range of wave frequencies. The present results are almost identical to the results from the linear wave theory and fully nonlinear numerical model over the range of wave frequency simulated. Among the three numerical results, the second-order results are in between the linear and fully nonlinear results. The fully nonlinear results are closest to the experimental data. That is to say, for the strongly nonlinear wave conditions, the fully nonlinear model is more accurate but time-consuming. On the other hand, the linear model is easily excuted and saves time but it is not applicable for the nonlinear cases. The second-order model is a compromise between linear and fully-nonlinear models. However, discrepancies are observed between the numerical results and experimental data. It is seen that all the potential flow models over-predict the wave elevations near the resonant frequencies. The reason for this is that the energy dissipation due to the viscosity of the flow are not considered by the potential flow models (Lu et al., 2010).

The resonance in the wave motion can be clearly identified by the peaks in the surface elevation. In this case, resonance occurs near three wave frequencies: 5.75 rad/s, 6.80 rad/s and 8.0 rad/s. The corresponding resonance modes are defined as Mode 1, Mode 3 and Mode 5 (Feng and Bai, 2015), respectively, in the case of beam sea. In Modes 1, 3 and 5, there exist 0.5, 1.5 and 2.5 wavelengths along the gap, respectively, as shown in Fig. 6. There exists only one wave crest and none wave trough in the gap in Fig. 6a, but one wave crest and two wave troughs are shown in Fig. 6b, and three wave crests and two wave troughs appear in Fig. 6c. It is also observed that the maximum surface elevation decreases with increasing mode number. The first resonance and third modes, with relatively large wave surface elevations, are the focus in this paper.

The influence of incident wave angle on wave elevations in the gap is investigated in this section.The configuration considered is the same as that used in Section 3, except that the incident wave direction is varied. Figure 7 shows the variations of the dimensionless wave crest with wave angular wave frequency in the gap at x=0.0 m under three wave directions. It can be seen that the largest wave crests in the first, third and fifth resonant modes occur at β being 60°, 30° and 0°, respectively. At the same time, when the wave direction is oblique, the maximum wave elevation at wave frequency instead of resonant frequencies does not necessarily occur at the gap centre. To quantify the location of the maximum wave elevation, the distribution of the dimensionless wave crest along the gap at frequencies close to the resonance frequencies of the first and third resonant modes are presented in Figs 8 and 9 for different incident wave directions. It was found that the wave direction does not affect the resonant frequencies. The resonant frequencies of the first and third modes always occur at 5.75 rad/s and 6.80 rad/s, respectively. The position of the maximum surface elevation along the gap always occurs at the gap center at resonance. From Fig. 8d, it can be seen that the distribution of the wave crests along the gap is symmetric for different wave frequencies as the incident wave angle is equal to 90°. Due to the action of the oblique wave, the location of the maximum wave elevation shifts upstream along the gap for the incident wave frequency that is smaller than the resonant frequency, but shifts downstream for the wave frequency that is larger than the resonant frequency as shown in Figs 8a-c. Such phenomenon is more apparent for smaller incident wave angle. The similar phenomena in the resonant Mode 3 can also be oberseved in Fig. 9.

To have a better understanding of the effects of ω and β on the wave-structure interaction, the contours of the maximum wave surface elevation are plotted on the ω-β plane. It appears that the incident wave angle β does not have effects on the resonance wave frequency. However, it affects the maximum wave surface elevation at resonance. The largest surface elevation appears when the wave direction is 60° for the first resonance mode. The reason may be due to that when the wave direction is oblique, part of the wave enters the gap directly. The oblique wave can be decomposed into the beam wave and head wave, which both can induce the water resonance in the gap as shown in Figs 8 and 9. The interaction of these two waves with different propagation modes may cause a greater surface elevation, which is a quite complicated process and deserved to be further investigated. From Fig. 10, it can also be observed that when the incident wave angles equal to 30° and 45°, the maximum surface elevations reach the largest at the third and fifth resonant modes, respectively. Some of them are even greater than that at the first resonant mode.

The effect of the wave direction on the water resonance in the gap between two side-by-side barges is investigated based on the second-order time-domain potential flow model. Two identical, three-dimensional rectugular barges subjected to different incident wave conditions are simulated. It is found that the distribution of the wave crest along the gap is symmetric about the gap center for all different incident wave directions as the incident wave is perpendicular to the box, i.e., the maximum wave elevations all occur at the gap center. The location of the maximum wave elevation along the gap shifts upstream for the lower freuqnecy waves, but downstream for the higher frequency waves by comparison with the resonant frequency wave. It is also observed that the incident wave angle does not influence the resonant wave frequency, but the corresponding peak surface elevation at resonance varies with the incident wave angle. The maximum wave crest in different resonant mode occurs at different incident wave angle due to the complicated interaction of beam wave and head wave.