Higher-order gap resonance and heave response of two side-by-side barges under Stokes and cnoidal waves

Coupled piston-mode fluid response and the heave motion of two identical barges in side-by-side configuration is studied under finite-depth and shallow-water waves using a two-dimensional fully nonlinear numerical wave tank. To understand possible critical responses of the gap flow and the floating barge, regular-wave conditions which are able to excite up to 5th-order nonlinear gap resonance and also the resonant heave motion of the barge are considered. In shallow-water waves, high-frequency oscillations, featured by secondary peaks in the time histories, are observed for both wave elevation in the gap and the heave motion of the barge. The shallow-water wave-induced 4th- or 5th-order gap resonance can be equally crucial as the 1st- and 2nd-order resonances due to finite-depth waves. At higher-order gap resonance, the higher-harmonic heave motion of the barge is negligibly small, in contrast to the gap-flow response. Compared with fixed barges, the free-heave motion of an upstream barge tends to increase the wave elevation in the gap in most of the resonant conditions, except at 1st-order gap resonance where the gap response is greatly reduced. When the resonant heave motion of a floating barge, either located upstream or downstream, is excited, significant barge motion is observed. However, the relative motion between the gap flow and the floating barge is seen to be very small, ascribed by small phase difference between the two. The present study suggests that the effects of heave motion and water-depth should be carefully considered in the design of side-by-side marine operations, and hiding the small bunkering ships behind the large receiving ships is regarded as a preferred arrangement during the bunkering operations in offshore and coastal environments.


Introduction
Recently, ship-to-ship bunkering has owned its popularity in fueling the LNG-fueled vessels outside the port. During the bunkering operations from the bunkering vessel to the receiving vessel, a key hydrodynamic issue is related to the occurrence of fluid resonance inside the narrow gap between the two vessels in close proximity under the action of waves. The fluid resonance elevates the water surface within the gap to a large extent and influences the dynamic responses of the vessels due to hydrodynamic interaction. It may also lead to wave impacts on the soft hoses which are typically used to transfer liquid fuels between the two vessels. Two types of resonances may matter in the gap, namely the sloshing-type resonance and the 'piston-mode gap resonance', while the present study focuses only on the latter.
Early studies towards fluid resonance in a confined region with the presence of a free surface started with the theoretical investigation based on linear potential-flow theory (Miao et al., 2000;Iwata et al., 2007;Zhu et al., 2008). Molin (2001) considered rectangular moonpools of large horizontal dimensions and determined the natural modes * Corresponding author. for the oscillation of the inner free-surface flows in both two and three dimensions. Subsequently, physical experiments were conducted to better understand the mechanism behind gap resonance and to provide validation data for numerical analysis Faltinsen, 2008, 2010;Faltinsen et al., 2011;Tan et al., 2014). Saitoh et al. (2006) conducted laboratory experiments to understand the characteristics of fluid resonance in the gap between two fixed barges exposed to waves and concluded that the resonant wave amplitude in the narrow gap is dependent on the body draft and gap width. Perić and Swan (2015) carried out model tests to characterize both the amplitude and nature of the excitation within the gap. Experiments were also carried out by Ning et al. (2018) to investigate the wave response in the gap between two barges of different drafts in incident waves. Some physical tests have also been carried out to study three-dimensional (3D) gap resonance problems (Fournier et al., 2006;Molin et al., 2009;Zhao et al., 2017;Chua et al., 2019).
It has been shown that linear potential-flow theory can predict the resonant frequencies with high accuracy and the resonant modes can https://doi.org/10. 1016/j.oceaneng.2022.112835 Received 19 July 2022; Received in revised form 27 September 2022; Accepted 5 October 2022 also be well captured (Molin, 2001;Faltinsen et al., 2007). However, the wave height in the gap at the resonance is largely overestimated when compared with the experiments, as linear potential-flow theory ignores the energy dissipation caused by the fluid viscosity, flow separation and turbulence. To overcome this, improved potential-flow methods with viscous corrections were applied to avoid non-physically large wave elevation within the gap at the resonance (Tan et al., 2014;Ning et al., 2015a,b;Zhao et al., 2018;Tan et al., 2019). Lu et al. (2011) presented two-dimensional (2D) numerical results of the wave loads on multiple floating bodies in close proximity and also their dependence on the wave frequency, gap width, body draft, body breadth, and the number of bodies based on both a viscous-flow model and a potentialflow model with the inclusion of artificial damping. However, they also pointed out that the artificial damping has to be calibrated by model tests or dedicated Computational Fluid Dynamics (CFD) simulations in order to achieve better prediction, which limits the application of potential-flow models to a certain extent.
With the advances of computation capacity and numerical algorithms, CFD methods based on Navier-Stokes equations have become a popular alternative to deal with the gap resonance problem (Lu et al., 2010;Chua et al., 2018;Jiang and Bai, 2020). Moradi et al. (2016) concluded that potential-flow models not only over-predict the resonant wave heights but also fail to predict the variation trends of the resonant wave height as function of water depth. Piston-mode wave resonance between a ship section and a bottom-mounted terminal was studied by Jiang et al. (2019), where it was also confirmed that the effect of fluid viscosity and the vorticity increases significantly with the incident wave amplitude. Some numerical investigations were also carried out by Gao et al. (2019aGao et al. ( , 2020a to simulate resonant fluid motions inside a narrow gap under various configurations and wave conditions. Besides the viscous dissipation, free-surface nonlinearity is also a key factor that can excite large resonant wave amplitudes in the narrow gap (Kristiansen and Faltinsen, 2008;Liu et al., 2016;Bonfiglio and Brizzolara, 2018;Gao et al., 2020c;Song et al., 2020). In addition to 1st-order resonance, 2nd-order resonances with significant secondary wave elevations and wave loads were also observed when the incident wave frequency ( 0 ) equals to half the natural frequency of the gap flow ( ) (Sun et al., 2010;Gao et al., 2019b;Li, 2019). The existence of 3rdorder gap resonance has been also identified in the literatures (Zhao et al., 2017;Chen et al., 2018;Gao et al., 2021;He et al., 2021;Jiang et al., 2021a;Ding et al., 2022a).
To date, few studies have considered the effect of body motion in the analysis of gap response between two floating structures in close proximity (Li and Teng, 2015;Li and Zhang, 2016;Milne and Graham, 2018;Li, 2019;Lu et al., 2020). However, large relative vertical motions between two vessels may occur during the bunkering operations at sea. Gao et al. (2021) studied the influence of the freeheave motion of an upstream box on the gap resonance under wave actions while the downstream box remains fixed. However, bunkering operations with the bunkering vessel hiding behind the receiving vessel might be a preferable configuration to avoid violent relative motion between the two vessels. As will also be shown later, the coupling between a heaving barge and gap flow can significantly change the flow physics and the responses of the coupled dynamic system.
Meanwhile, a large number of marine operations are carried out in the coastal areas in which appreciable shallow-water-depth effects on the nonlinear wave-structure interactions are expected. The fluid resonances within a narrow gap between a box and a vertical wall were investigated by Jiang et al. (2021b). The incident wave was modeled by cnoidal wave theory to consistently account for shallow water effects and up to 3rd-order gap resonances were identified for finite water depths. Under shallow-water wave conditions, Ding et al. (2022a) has demonstrated the occurrence of even higher-(e.g. 4th-and 5th-) order gap resonances in between two fixed identical barges, where significant wave crest occurred due to the nonlinearities of the incident wave and the nonlinear wave-structure interactions. Large wave crest may lead to wave impact on the soft hoses through which the lowtemperature and high-pressure liquid fuel is transferred, threatening the safety of the operations. Moreover, it has also been observed in Ding et al. (2022a) that the wave elevation in the gap showed non-negligible high-frequency oscillation due to significant higher-order harmonic response.
To the best of the authors' knowledge, there has been no studies in the literature on the ship-to-ship interaction which simultaneously take into account the effects of ship motion, shallow-water depth, freesurface nonlinearity and the viscous flow separation around the barges. With this background information, the present study will extend and apply the 2D fully nonlinear numerical wave tank (NWT) established in Ding et al. (2022a), aiming to answer, among others, the following research questions: (1) how will the free-heave motion of the barges influence the higher-order gap resonance under both finite-depth and shallow-water waves? (2) how much can one gain by deploying the bunkering vessel on the leeside of a receiving vessel? (3) is the violent motion of the gap flow able to excite high-frequency motion of the vessel, which may increase the fatigue damage of the transferring hoses?
In this study, the numerical setup with one fixed square-shape barge and the other identical free-heaving barge is designed to investigate the effects of heave motion and water depth on the gap resonances and wave-structure interactions. The coupled natural frequencies of the gap flow and the heave motion are obtained through free decay tests, and six critical regular-wave conditions are designed to excite the gap resonances and the resonant heave motion of the barge, covering both finite-depth and shallow-water waves. Higher-order gap resonance with consistent consideration of water-depth effect is investigated by analyzing the gap resonance formed by these two identical square barges, as well as the response of the heaving barge, under the actions of Stokes and cnoidal waves.
The remainder of this paper is organized as follows: in Section 2, the numerical model applied in this work is briefly introduced. Then the computational setup of the two-dimensional (2D) numerical wave tank (NWT) is described in Section 3, followed by the validation studies in Section 4 against existing experimental data and numerical solutions. Section 5 shows the determination of the coupled natural frequencies of the gap flow and the heave motion through free decay tests, as well as the design of incident wave conditions. Numerical results and discussions regarding the effects of water depth and free-heave motion on the gap resonance are presented in Section 6. Finally, conclusions are drawn in Section 7.

Numerical model
The Unsteady Reynolds Averaged Navier-Stokes (URANS) turbulence model is applied to consider the viscous dissipation and the Volume of Fluid (VOF) method is used to predict the distribution and the movement of the interface between the air and water phases. The overset meshing technique is employed to track the motion of the heaving barge. Interested readers are referred to the book (Ferziger and Perić, 2002) or the STAR-CCM+ user guide (Siemens, 2020) for more details on the schemes and solvers regarding the numerical model.

Governing equations
The URANS turbulence model in STAR-CCM+ solves the Navier-Stokes equations to describe the motion of the fluid continuum. The conservation of mass and momentum can be expressed as: where is the control volume bounded by the closed surface ; and are the fluid and grid velocities, respectively; is the unit vector normal to and pointing positively outwards; represents time; is the pressure, and is the fluid density; and are the unit vectors in the th and th directions, respectively; is the gravity vector; is the viscous stress tensor. denotes the source term for momentum and is non-zero where wave forcing is applied.
To determine and update the position of the free surface, the VOF method is employed, where the volume of fraction in grid cells is defined as: Then the spatial variation of any fluid property (e.g., the density , and the dynamic viscosity ) can be represented as a weight average of : in which, the subscripts ''water'' and ''air'' denote the property of the corresponding phase, respectively. In STAR-CCM+, the transport equation of the volume fraction is only solved for the first phase, defined as water here, when only two phases are considered, which is represented as: here, = 0.5 is adopted in the present study to capture the free surface. The volume fraction of air is thus adjusted based on that the sum of the volume fractions of both phases is 1. Wave forcing is applied at both ends of the numerical wave flume along the wave propagation to absorb the reflected and transmitted waves, where the solution of the momentum equations (Eq. (2)) is enforced to the solution of the Navier-Stokes equations towards the prescribed nonlinear wave solution by modulating in the source term (Kim et al., 2012): here denotes the forcing coefficient and its optimum value to achieve best wave absorption can be determined by numerical simulations (Perić and Abdel-Maksoud, 2016). represents the current numerical solution of the transport equations and * is the prescribed incidentwave solution. The forcing coefficient varies smoothly within the forcing zone as follows: in which is the maximum value of the forcing coefficient. * is the normalized distance from the cell considered to the boundary in the direction of wave propagation with respect to the length scale of the forcing zone.

Body and mesh motions
The Dynamic Fluid Body Interaction (DFBI) module in STAR-CCM+ is employed to simulate the motion response to forces and moments exerted by the fluid. In this work, the floating barge is confined solely to heave motion and the translation of center of mass is formulated based on Newton's second law.
The mesh motion of the computational domain is calculated by the overset meshing technique, which couples with the DFBI module to update the mesh position around the heaving barge. Meshes in the overset region move rigidly with the barge motion, where the solution is coupled with the flow domain by linearly interpolating donor-cell values to provide data for acceptor cells (Siemens, 2020).

Numerical implementations
The finite volume method is applied to solve the governing equations (Eqs. (1), (2) and (5)) for each cell. The angle factor in the High Resolution Interface Capturing scheme (HRIC) in the VOF multiphase model is set up as 0.15 to reduce the numerical diffusion on the free surface (Muzaferija, 1998). The standard low-Re − turbulence model with all-+ wall treatment is used to prevent excessive growth of turbulent viscosity near the free-surface as the initial values of the turbulent kinetic energy = 10 −5 J/kg and the turbulent dissipation rate = 10 −4 m 2 ∕s 3 are prescribed (Perić, 2018). The applied spatial and temporal integrations are both of second-order accuracy. The implicit unsteady segregated solver is employed to linearize and solve the resulting coupled system of equations by a multigrid method. Ten outer iterations are performed for each time step, which consists of solving the governing equations for the velocity components, the pressurecorrection equation, and the transport equations for the volume fraction of water, as well as and . In addition, the time step is adaptive to satisfy the Courant-Friedrichs-Lewy (CFL) condition, ≤ ∕ , where and are the mesh size and the relative velocity of the fluid to the mesh, respectively. To obtain accurate and stable results, the largest convective Courant number near the free surface is set to around 0.1 in all simulations.

Computational setup
A two-dimensional (2D) fully nonlinear numerical wave tank (NWT) is established to study the effects of heave motion and wave depth on the gap resonances and wave-structure interactions, which is shown in Fig. 1. The height of the wave flume is kept three times of the constant water depth ℎ and the total length of the wave flume adapts to the incident wavelength . A Cartesian coordinate system is applied with the origin set on the still-wave line (SWL) and at the middle of the wave flume. The incident waves propagate along the x-axis and elevate along the z-axis. Two identical square barges are symmetrically located at the center of the flume to the origin with a gap width = 0.1 ℎ. The height and breadth of the barges are both equal in value to the water depth ℎ with an initial draft = 0.5 ℎ. For brevity, the heaving barge is hereafter defined as 'Barge A' and the fixed barge is defined as 'Barge B'. To facilitate the comparison between different setups, the following three systems are defined: • Upstream heaving (UH) system: the upstream barge is free to heave while the downstream barge is fixed; • Downstream heaving (DH) system: the downstream barge is freely heaving with the upstream barge fixed; • Fixed system: both barges are fixed.
It is noted that only the DH system is displayed in Fig. 1 and the model setup of the UH system is almost the same except for a reverse wave propagation (along the x-axis).
The Stokes waves are generated by prescribing the volume fractions and the velocity according to the 5th-order Stokes theory (Fenton, 1988) while a 5th-order periodic solution of Korteweg-de Vries equation (Korteweg and De Vries, 1895) is used to generate the cnoidal waves within the wave forcing zones at both ends of the numerical wave tank along the wave propagation direction, with the lengths of the forcing zones to be twice the incident wavelength . The no-slip boundary conditions are imposed both on the barges and the bottom of the flume. The top boundary of the domain is defined as a pressure outlet. Only one layer of meshes is employed along the transverse direction of the 2D NWT with the symmetry boundary conditions applied on the side-walls of the flume. The overset mesh is set up on the boundaries of the overset region that moves simultaneously with the heaving Barge A. Numerical simulations start from the still-water state with predefined initial conditions of the velocities and the hydrostatic pressures of the incident 5th-order Stokes or cnoidal waves.  and are the breadth and the draft of the barges, respectively. ℎ is the water depth and is the gap width. 'SWL' represents the still-wave line.
Trimmed mesh is employed in the entire domain and a typical mesh configuration is displayed in Fig. 2, with finer resolution at necessary areas. Near the free surface, 20 cells are distributed along the vertical z-direction per wave height and the mesh size along the horizontal x-direction is determined to keep a constant aspect ratio ∕ = 8, which was tested based on careful verification to produce high-quality incident waves. Details for the convergence studies on the mesh and time step can be seen in Ding et al. (2022a). In this study, more than 1,000 time steps are adopted per incident wave period to satisfy the temporal convergence for all the numerical simulations. The meshes are further refined in the narrow gap and close to the barges with adequate prism layers to better resolve the flow field due to flow separation on both barge surfaces in the gap and also at the sharp corners of the barges facing the gap. In addition, meshes are refined around the barges to record the position of the heaving barge, while similar refinement has also been applied in the overset region to achieve accurate data interpolation between the cells in the overset domain and those in the background domain.

Validation of the numerical model
To the best of the authors' knowledge, no experimental results are available in the literature which has similar setup as the present model, thus some widely-used experimental results were utilized to validate the numerical model in two steps, consisting: (I) a fixed system with two identical barges under wave actions; (II) a single free-heaving barge in waves at different wavenumbers.
In the first case, the numerical results obtained by simulations using the same computational setup described in the previous Section 3 were compared with the existing experimental data (Saitoh et al., 2006) and other numerical results (Lu et al., 2008;Feng et al., 2017) regarding the free-surface elevation in the narrow gap. The validation results have been presented in Ding et al. (2022a,b).
For the second case, heave motion amplitude of a floating barge was compared with the experimental measurements by Rodríguez and Spinneken (2016) and the numerical results by Gao et al. (2021). Physical experiments were carried out by Rodríguez and Spinneken (2016) in a 2D wave flume with a length of 63 m and a width of 2.79 m to investigate the free heave motion of a single rectangular barge under regular waves. The water depth was 1.25 m. The barge was located almost in the middle of the wave flume with a breadth of = 0.5 m and a draft of half breadth. Stokes waves with the steepness 0 = 0.10 were considered at various wavenumbers . The numerical wave tank adopted in this case (not shown here) was similar to that in

Free decay tests
According to Faltinsen and Timokha (2015) and Fredriksen et al. (2015), when the fluid bulk motion in the narrow gap is coupled with the floater motion, the natural frequency of the gap flow for the coupled system is different from the uncoupled natural frequency of the gap flow that can be estimated by considering a fixed system (e.g. Molin (2001), Faltinsen et al. (2007)). Therefore, this study will start with free-decay tests in the NWT to obtain the coupled natural frequencies of the water bulk in the narrow gap and the heave motion. Two different free-decay tests will be performed: • Released barge: Barge A is released from an initial vertical position of 0 above its equilibrium position and heaves freely, while Barge B is fixed aside. The fluid motion in the gap is excited by the heaving barge. • Released gap fluid: the water surface in the gap is released from an initial vertical position of 0 above the still-water level. The heave motion of Barge A is excited by the gap flow with the fixed Barge B aside.
For the case of released barge, Figs. 4(a) and 4(b) present the nondimensional time histories of the wave elevation in the gap and the free-decay heave motion of Barge A, respectively. It is observed from Fig. 4(b) that the initial position of the center of gravity of Barge A is not strictly 0 due to a slightly smaller draft in the numerical model compared with the originally specified draft of = 0.5ℎ. This small numerical error is due to the imbalance between the specified mass of the barge and the estimated buoyancy in STAR-CCM+ at this draft, but is found to be negligibly small as 0.001 . Therefore, a downshift of the free-decay result is carried out to get rid of the non-zero mean heave displacement. By using a Fast Fourier transform (FFT), the amplitude components of both time responses are obtained in frequency domain.
From Fig. 5, it can be observed that the coupled natural frequencies of the gap flow ( ) and the heave motion ( ) are roughly 0.9 Hz (corresponding to non-dimensional wavenumber ( ℎ) ≈ 1.732) and 0.6 Hz (corresponding to non-dimensional wavenumber ( ℎ) ≈ 0.967), respectively, which are listed in Table 1. If a fixed system is considered, the corresponding uncoupled natural frequency of the gap flow ( ), as given in Ding et al. (2022a), is 0.84 Hz. Thus, the coupling with the barge motion has increased the natural frequency of the gap flow.
For the case of released gap fluid, the non-dimensional time histories of the free-decay wave elevation in the gap and the heave displacement of Barge A are shown in Fig. 6(a) and Fig. 6(b), respectively. The corresponding normalized amplitude components are obtained in Fig. 7(a) and Fig. 7(b), respectively. It is shown that the natural frequency of the gap flow is also 0.9 Hz, which cross-verifies the prediction of the natural frequency of the gap flow in the free-decay test with released barge. However, the heave motion is hardly excited due to the relative small volume of water bulk in the narrow gap, in which case trivial amplitude component can be observed in Fig. 7(b) at natural frequency of the heaving barge.

Selection of critical wave conditions
Only a few representative critical waves which excite resonant gap responses or resonant heave motion of the barge will be considered is the breadth of the barge. represents the gravitational acceleration and is 9.81 m∕s 2 in the present study. denotes the wave elevation in the gap. and 0 denote the position of the center of gravity and the initial released position of Barge A.

Table 1
Natural frequencies of the coupled heaving system and the uncoupled fixed system. and are the natural frequencies of the gap flow for the heaving systems and the heave motion of Barge A, respectively.
is the natural frequency of the gap flow for the fixed system.

Natural frequency
Value (Hz) 0.9 0.84 0.6 in later analysis. Six incident waves (waves A to E, and M), marked by the stars in Fig. 8, are selected according to the coupled natural frequencies of the gap flow and the heaving barge. In Fig. 8, the thick solid curve ∕ℎ denotes the breaking wave limits and the thin solid curve representing the Ursell number = 0 2 ∕ℎ 3 is a suggested threshold between Stokes and cnoidal theories' areas of application, which is adopted as ≈ 40 according to Fenton (1990).
Among the six wave conditions, four Stokes waves (waves A, B, C and M in Fig. 8) with equal non-dimensional wave height 0 = 0.048 ℎ but different normalized wavelengths ∕ℎ, are considered as incident waves to represent the scenarios at finite water depths. Two cnoidal waves (waves D and E in Fig. 8) with a larger wave height 0 = 0.2 ℎ and longer wavelengths are considered to represent the shallow-water depth conditions. The non-dimensional wave height 0 ∕ℎ, wavelength ∕ℎ, * = ∕ 0 as the ratio between the natural frequency of the gap flow and the primary incident wave frequency 0 , together with the normalized heave natural frequency * = ∕ 0 of Barge A are summarized in Table 2 for the six wave conditions. It has to be noted that waves A, B, C, D and E are designed so that their primary frequencies are equal to , ∕2, ∕3, ∕4 and ∕5, representing 1st-, 2nd-, 3rd-, 4th-and 5th-order resonances of the gap flow in the coupled systems (UH system or DH system), respectively. The wave M is designed to excite the 1st-order resonant heave motion with the primary frequency 0 = = 2∕3 . In addition, since the natural frequency of the gap flow for the fixed system is different from that for the heaving systems, the designed wave parameters also deviate from those listed in Table 2. Table 3 summarizes the wave parameters of the five wave conditions (wave A, B, C, F and G) for the fixed system, which correspond to the incident waves that are able to excite 1st-, 2nd-, 3rd-, 4th-and 5th-order gap resonances for the fixed system, respectively. Further details regarding the wave conditions and the subsequent discussion for the fixed system can be found in Ding et al. (2022a).
The same numerical model as that described in Section 2 is further adopted to investigate the wave elevation within the gap and the heave motion of Barge A under the action of Stokes and cnoidal waves described in Fig. 8 and Table 2. Geometries of the barges and the gap are also the same as those of the validation cases in Ding et al. (2022a,b). The computational setup introduced in Section 3 is applied. The quality of the incident waves has been verified by running another simulation without the barges and comparing the obtained wave time-series and profiles with the corresponding nonlinear wave theories.

Results and discussions
To study the coupled heave motion and gap response excited by waves at different water depths, the free surface elevation in the narrow  represents the upper-limit wave height and represents the Ursell number. and 0 denote the wavelength and the wave height, respectively. ℎ is the water depth.

Table 2
Wave parameters of the six wave conditions for the heaving systems (UH system and DH system). The values of 0 ∕ℎ, ∕ℎ, * = ∕ 0 and * = ∕ 0 are summarized. and are the natural frequencies of the gap flow and the heave motion of Barge A, respectively. 0 is the primary frequency of the incident wave. 0 is the crest-to-trough wave height and is the incident wavelength. ℎ is the water depth.  Table 3 Wave parameters of the five wave conditions for the fixed system. The values of 0 ∕ℎ, ∕ℎ and * = ∕ 0 are summarized. is the natural frequency of the gap flow for the fixed system and 0 is the primary frequency of the incident wave. 0 is the crest-to-trough wave height and is the incident wavelength. ℎ is the water depth. gap is monitored for the UH system and the DH system under the action of incident waves defined in Fig. 8 and Table 2. To obtain the harmonics of the piston-mode response in the gap, Fast Fourier transformation (FFT) is applied to the periodic part of the time-domain results, which will be presented and discussed in the following Section 6.1. Next, the time histories of the heave displacement of Barge A are obtained together with the corresponding harmonic amplitudes in Section 6.2. In addition, analyses are made in Section 6.3 to understand the effect of resonant heave motion for the considered heaving systems. Finally, direct comparisons are made to study the effects of water depth and heave motion on the gap resonances and wave-structure interactions.

Finite-depth results (Stokes waves A, B and C)
The time histories of normalized wave elevatioñ= 2 ∕ 0 in the narrow gap under Stokes waves (A, B and C in Fig. 8) at finite-water depths for the UH system and the DH system are shown in Fig. 9, and the corresponding normalized harmonic amplitudes̃= 2 ∕ 0 with respect to the non-dimensional frequency ∕ 0 are presented in Fig. 10. The non-dimensional natural frequency of the gap flow * = ∕ 0 is denoted by the blue vertical dash lines. The results for the fixed system are also included in Figs. 9 and 10 for comparison purpose.
Under wave condition A (1st-order gap resonance), the largest wave amplitude in the gap is observed in the fixed system, followed by those in the UH system and DH system, respectively. See the time series of wave elevation in Fig. 9(a) for comparison. It is also shown that the wave elevations in the gap for the DH system and fixed system are rather in phase with each other while the UH system has an almost reversed phase to the other two systems. The occurrence of 1st-order gap resonances for all the three systems is confirmed by the corresponding harmonic amplitudes depicted in Fig. 10(a). The 1st harmonics of gap response in UH system and DH system have similar magnitudes, but are much smaller than that in the fixed system.
However, a different observation has been made under wave B (2ndorder gap resonance), where the UH system shows much larger gap response than the other two systems. See the comparison of waveelevation time series in Fig. 9(b). The corresponding harmonic amplitudes of gap response are also shown in Fig. 10(b) for the three systems. Both the 1st and 2nd harmonic amplitudes for the UH system are the largest among the three systems.
Similarly, the 3rd-order resonances are excited by wave C with 8 Y. Ding et al.  is the frequency and 0 is the primary incident-wave frequency.
is the amplitude of the th harmonic of wave elevation in the gap.   Fig. 11 presents the time series of normalized wave elevatioñ= 2 ∕ 0 in the gap under cnoidal waves (D and E in Fig. 8) for the three systems and the corresponding normalized harmonics̃= 2 ∕ 0 are displayed in Fig. 12. Strong higher-order gap resonances, e.g. 4thorder for wave D and 5th-order for wave E, are observed at the troughs of the wave elevation in the gap in Fig. 11 and also at the corresponding harmonic amplitudes in Fig. 12. It is seen from Fig. 12 that the harmonic amplitudes at the resonances are comparable to those of the primary and secondary components for the UH system and the fixed system, while the harmonics at the resonances for the DH system are close to the 2nd-harmonic amplitudes but much smaller than the 1st-harmonic results. Besides, the UH system exhibits the largest wave elevations, followed in consequence by the fixed system and the DH system.

Shallow-water results (cnoidal waves D and E)
In order to understand the relationship between the harmonic components in Fig. 12 and their contributions to the total time-series gap responses shown in Fig. 11, the time-domain results are decomposed as the linear superposition of each harmonic component ( ) = ∑ ( ) = ∑ cos (2 + ). Here, is the time-series gap response of each harmonic component and is the amplitude of the gap response, while and represent the frequency and the phase angle of each harmonic component, respectively.
The phase angles of the harmonic components of wave elevation in the gap under wave conditions D and E for the three systems are compared in Fig. 13. The phase angles of the harmonic responses are almost the same in the two heaving systems when the corresponding frequencies are smaller than the natural frequency of the gap flow (denoted by * in the figure), while they start to deviate from each other when approaching the natural frequency. At the natural frequency, a sudden change of sign of the phase angles are seen for both the UH system and DH system. However, the same behavior has not been observed for the fixed system, where the harmonic components are tuned to be in phase with the incident wave (understood by the 0 • phase angle at * according to the above definition of the phase angle). In fact, for the UH system and DH system, sign change in the phase angles of the harmonic components at the natural frequency of the gap flow have also occurred for Stokes waves A, B and C, which is not shown here due to space limit.
Even though the nonlinear dynamic and coupled system (between barge and gap flow) we are investigating is rather complex, we can make analogy to typical mechanical systems that we are familiar with to understand some of its behavior. For instance, the response of massspring-damper system is in phase with excitation force at low excitation frequency, 180 deg out of phase with excitation at higher frequency, and 90 deg out of phase at resonant frequency (appears as sign change if the phase angle is defined between −180 deg and 180 deg). On the other hand, it is known from a linear theory that the wave excitation force is in phase with the wave elevation at the location of the gap. Typically, the sign change of the phase angles at certain frequency indicates the onset of gap resonance.
The time series of each harmonic component 2 ∕ 0 is shown in Fig. 14(a) for the UH system under wave D (4th-order gap resonance). The maximum crest occurs at where all harmonic components have positive contribution to the total wave crest. However, at the same time instant, as depicted in Fig. 14(b), the 4th-harmonic component has a cancellation with other harmonics in the DH system, and thus the maximum crest occurs at a smaller time instant. Similar behavior has also been observed in the case of 5th-order gap resonance (wave E), and the results will not be shown here due to limited space.

Finite-depth results (Stokes waves A, B and C)
The time series of normalized heave motioñ= 2 ∕ 0 of the floating Barge A excited by incident Stokes waves (A, B and C in Fig. 8) are shown in Fig. 15 for the two heaving systems. The normalized harmonic amplitudes of the heave motioñ= 2 ∕ 0 are presented in Fig. 16 as function of non-dimensional frequency ∕ 0 . The red vertical dash lines in Fig. 16, corresponding to * = ∕ 0 , are used to locate the coupled heave natural frequency.
It is observed from Fig. 15 that the heave motion of Barge A for the UH system is always larger than that for the DH system under the three waves. Furthermore, the heave motion is much smaller under wave A (1st-order gap resonance) than those in waves B (2nd-order gap resonance) and C (3rd-order gap resonance), which have much longer wavelengths. For example, for the UH system, the heave motion under wave B is about 4 times of that under wave A. Based on the parameter given in Table 2, the non-dimensional initial drafts of the barge under waves A, B and C are ≈ 0.85, 0.36 and 0.20, respectively. Here = 2 ∕ is the primary wavenumber. Due to the exponential decay of the incident wave kinematics, the barge will experience smaller excitation force in shorter waves. However, based on our estimation of the wave excitation forces, this can only explain a small part of the differences of the heave motions under different wave conditions. Another explanation could be the motion coupling between barge and gap flow at resonant conditions. Since wave A corresponds to the 1storder resonance of the gap flow, the energy from the incident wave is largely 'trapped' by the gap-flow motion while the barge motion gains relative less energy. The observed large-amplitude gap response at the 1st-order gap resonance illustrated by the results for wave A in Fig. 10(a) seems to support this argument.
As one can observe by comparing Figs. 10 and 16, the corresponding higher-harmonic heave motion of the barge is negligibly small when higher-order gap resonance occurs (marked by the blue vertical lines in the figures), contrasting to the gap-flow response.
As another note for Fig. 15, the heave motion does not show clear secondary peak in any of the three wave conditions as it does in the corresponding gap-flow response depicted in Figs. 9(b) and 9(c). However, it is expected that secondary peaks will still appear if waves with higher nonlinearities are considered.
Different from wave conditions A, B and C, high-frequency oscillations are clearly observed at the troughs of the heave-motion time series in Fig. 17. Those oscillations can potentially increase the fatigue damage of the bunkering system, for instance, at the connection points between the transferring hoses and the vessels. As shown in Fig. 18, since the energy of the incident cnoidal waves is more dispersed into higher-harmonic components, significant higher-order heave responses have been excited in the two shallow-water waves, featured by the secondary peaks in the time series of heave displacements in Fig. 17.
When higher-order gap resonance occurs, despite of significant contribution of higher-harmonic gap responses to the total gap response, the corresponding higher-harmonic heave motion is almost zero. Similar observation has been made for Stokes waves A, B and C corresponding to 1st-, 2nd-and 3rd-order gap resonances, as has been discussed in Section 6.2.1.
The phase angles of the harmonic components of the heave motion of Barge A under cnoidal waves D (4th-order gap resonance) and Y. Ding et al. E (5th-order gap resonance), are displayed in Fig. 19. Here, similar harmonic-decomposition is used and the total heave responses are regarded as the linear superposition of each harmonic component of the heave response cos (2 + ). Here, is the amplitude of the heave displacement of Barge A, and represent the frequency and the phase angle of the harmonic component, respectively. It is shown that the phase angles follow a similar trend as those of wave elevation in the gap for both heaving systems, and show a sign change at the natural frequency of the gap flow indicating the onset of resonant motion of the barge.

Gap and heave responses at natural heave frequency
Wave M, which excites 1st-order heave resonance, has primary frequency equal to the natural heave frequency of the floating barge. Fig. 20 shows the time history of normalized wave elevatioñ= 2 ∕ 0 under wave M for the three systems, with the corresponding harmonic amplitudes presented in Fig. 21. The time history and harmonic amplitudes of the heave response are presented in Fig. 22 and Fig. 23, respectively.
The frequency of Wave M is in between waves A and B, and does not coincide with the natural frequency of the gap flow or an integer multiplier of . Compared with waves A and B whose results have been shown in Figs. 9(a) and 9(b), the overall gap response under wave M is smaller, because no gap resonance is evoked in this case. Comparing the results in Figs. 21 and 10(a), both the 1st-and 2ndharmonic gap responses under wave M are smaller than those under wave A. Comparing Figs. 21 and 10(b), one can also see that even though the 1st-harmonic amplitude under wave M is larger than that in wave B, the corresponding 2nd-harmonic response is negligibly smaller, leading to a lower overall gap response under wave M.
In contrast to gap-flow response, the overall barge motion under wave M is larger than those under waves A and B (results shown in Figs. 15(a) and 15(b)), featured by a higher 1st-harmonic amplitude in Fig. 23 than those under wave conditions A and B (corresponding results in Figs. 16(a) and 16(b)). However, the higher-harmonic heave responses are also found to be negligibly small under wave M.

Summary of the water-depth effects
Direct comparison of wave-elevation response inside the gap and heave response of Barge A for two heaving systems will be made between different wave conditions with finite and shallow-water depths. To shed some light on the engineering application with ship-to-ship operations under offshore and coastal environments, the maximum wave crest in the gap, the maximum barge motion and the maximum relative motion between the gap flow and the floating barge are summarized in this section to assess the effect of water depth. Fig. 24 demonstrates the normalized wave crest̃= 2 ∕ 0 in the gap under the Stokes and cnoidal waves considered in the present study (waves A to E and M) for the two heaving systems. The results for the fixed system are also included for comparison. Table 4 summarizes and harmonic amplitudes of wave elevatioñ= 2 ∕ 0 ( = 1, … , 8) in the gap for the UH system under the six waves considered. Here represents the average value of the maximum wave crest. It is observed from Fig. 24 and Table 4 that, for the UH system, the wave crests in the gap under cnoidal waves D and E which represent 4th-and 5th-order gap resonances respectively, are comparable to those under Stokes waves A (1st-order gap resonance) and B (2nd-order gap   resonance). This is due to the appreciable energy dispersion in higher harmonics for the shallow-water waves. In most cases and except for the 1st-order gap resonance, the largest wave crests in the gap are observed in the UH system, which is a more dangerous situation during the bunkering operations concerning the gap resonance. At the 1st-order gap resonance, the fixed system generates the largest wave elevation in the gap. The gap response in the DH system is the smallest under all wave conditions, indicating that it is a better choice to make use of the sheltering effect by hiding the small bunkering ship (e.g. a carrier ship) behind the large receiving ship (e.g. a container ship) to avoid violent free-surface elevation in the narrow gap at the offshore or coastal bunkering area.

Relative motion between gap flow and the heaving barge
The normalized heave displacements of Barge A under all six wave conditions considered for two heaving systems are demonstrated in Fig. 25 and the corresponding harmonic amplitudes̃= 2 ∕ 0 ( = 1, … , 8) for the UH system are listed in Table 5. represents the time-averaged maximum heave amplitude. As seen from Fig. 25 and Table 5, the largest heave motion of the floating barge is observed under wave M (1st-order heave resonance) for the UH system as a result of the 1st-order heave resonance, which is mainly contributed by the 1st-harmonic response 1 . The maximum heave motion under cnoidal waves D and E are close to those under Stokes waves M (1st-order heave resonance) and B (2nd-order gap resonance), making it necessary to consistently consider the water-depth effect in such operations under coastal environments. The maximum heave motions in waves D and E have been contributed by not only 1st harmonic but also the higherharmonic responses, as one can tell from the considerablẽfor > 1. Nevertheless, since the heave motion is not fully excited by the incident wave A representing the 1st-order gap resonance, much smaller heave crests can be seen for both heaving systems than those by the other wave conditions. In addition, the maximum barge motion for the DH system is always smaller than those for UH system, which roughly equals to 1.0 due to less incident wave transmitted through the upstream fixed barge to the downstream. The non-dimensional heave responses̃are expected to be greater than 1.0 under incident cnoidal waves D and E as a result of the crest-trough asymmetry of the wave elevation.
When wave slamming on transferring hoses is concerned, it is more relevant to study the relative motion ( ) = ( ) − ( ) between the gap flow and the structure. When ( ) is larger than a threshold, e.g. initial         represents the average crest of time-series relative motion between gap and heave responses. distance between the soft hose and the water surface, severe wave slamming may occur on the soft hoses between two vessels if the relative velocity between the gap flow and barge motion is also large. Therefore, the normalized relative motion between gap flow and the heaving Barge A, defined bỹ= 2 ∕ 0 , is presented in Fig. 26 for the two heaving systems under the Stokes and cnoidal waves. Here, denotes the average crest value of ( ). It is interesting to see from Fig. 26 that are very small for both UH system and DH system under wave M, where violent resonant barge motion is excited. It indicates that the gap flow is more in phase with the heave motion of the barge. In term of relative motion, the UH system still represents a more critical configuration than the DH system, and the largest occurs under wave A where 1st-order gap resonance occurs.

Conclusions
A two-dimensional (2D) fully nonlinear numerical wave tank (NWT) is established to numerically investigate the piston-mode gap resonance in the narrow gap between two identical barges under Stokes and cnoidal waves based on the CFD package STAR-CCM+ . Six regular waves, including four Stokes waves and two cnoidal waves, are considered to study the heave motion and the water-depth effects on the fluid response in the gap and also the heave response of the moving barge. Three Stokes waves are designed to excite the 1st-, 2nd-and 3rd-order gap resonances. The fourth Stokes wave is able to excite the 1st-order resonant heave motion of the floating barge. In addition, two cnoidal waves to excite 4th-and 5th-order gap resonances are also considered. Two heaving systems are modeled, namely the upstream heaving (UH) system with an upstream freely-heaving barge and the downstream heaving (DH) system with a downstream freely-heaving barge. The results for the fixed system are also included for comparison. Through this study, some conclusions can be drawn: • The present fully nonlinear numerical wave tank (NWT) is able to capture both the free-surface nonlinearity and viscous dissipation in the gap, which is essential in accurate prediction of resonant gap response. The heave motion of the moving barge can also be well resolved by the present numerical model. • Due to more dispersed energy into higher harmonics, it is easier for the shallow-water wave to trigger important higher-order (e.g. 4th-or 5th-order) gap resonance that might be equally crucial as the 1st-and 2nd-order resonances. • At higher-order gap resonance, the higher-harmonic heave motion of the barge is negligibly small, in contrast to the gap-flow response. • Compared with the fixed system and the DH system, the UH system leads to the largest gap response under almost all the considered resonant conditions, except for the case at 1st-order gap resonance where a significantly larger gap response occurs in the fixed system, 'trapping' the energy of the incident waves. • In general, the heave motion and gap response are significantly smaller in a DH system due to the sheltering effect. One exception is at the 1st-order gap resonance, where the gap response and the barge motion of the DH system are only slightly smaller than those in the UH system. • In the considered shallow-water waves, high-frequency oscillations, featured by secondary peaks in the time histories, are observed not only for the wave elevation in the gap but also the heave motion of the barge in both UH system and DH system, which may potentially increase the fatigue damage of the soft hoses hanging between the two barges. • When the heave resonance of the barge is excited, the heave motion is seen to be the largest among all the six wave conditions considered in this study. However, this heave-resonant condition leads to very little relative motion due to the fact that the gap flow and barge motion tend to be in phase with each other.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.