Fully nonlinear analysis of second-order gap resonance between two floating barges

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Fully nonlinear analysis of second-order gap resonance between two floating barges Yajie Li Research Center of Fluid Machinery Engineering and Technology, Jiangsu University, Zhenjiang 212013, China

a r t i c l e

i n f o

Keywords: Second-order gap resonance Boundary element method Fully nonlinear numerical wave tank

a b s t r a c t A two-dimensional boundary element method (BEM) based on the fully nonlinear potential flow theory is adopted to study second-order wave resonance in a gap formed by twin barges. Second-order resonance is related to second-order effects of wave elevations and hydrodynamic forces, which have components of the sum or difference frequency. Firstly, forced body motions (heave, sway, and roll mode) excited second-order resonance is studied with regard to different gap width and resonant mode. Then incident wave induced second-order resonance in the gap is simulated. It is found that second-order resonance is more easily provoked and significant when the gap width over draught ratio is large.

1. Introduction and background Parallel placed multiple floating bodies are common configurations in naval architecture and ocean engineering. There is water confined between the floating bodies. The water motions in the gaps resemble liquid sloshing in a tank. Thus, the terms used in tank sloshing problems like natural frequencies, natural modes and resonance are adopted to describe the liquid motions in the gap. The resonance of liquid motion in the gap is commonly briefed as gap resonance. The natural frequency 𝜔n in a gap is defined for a fixed configuration, i.e. gap width and depth. The natural frequency of the lowest resonance mode is called piston mode, denoted as 𝜔0 . The natural frequency of other sloshing mode is denoted as 𝜔i , i = 1, 2,…. The order i is determined by the corresponding mode shape in the gap. The mode shapes are characterised by the presence of standing waves. First-order resonance is the classic resonance, which can occur at a natural frequency 𝜔n in the gap. Second-order resonance depends on second-order effects, which have components of the sum or difference frequencies. The liquid motion in the gap at resonance will also result in a large wave elevation and high wave load on the surrounding structures. Considerable studies have been carried out focusing on the resonance phenomenon via different methods within potential flow theory. For researches in the frequency domain, Molin [1] studied the resonant behaviour of piston and sloshing modes in moonpools for 2D and 3D cases. The results were given for both the natural frequencies and the associated mode shapes for a range of moonpool draught to width ratios. Newman [2] analysed the low-frequency piston mode for moonpools with slowly varying cross-sections. Three different configurations were examined, namely a simple torus with one moonpool, a structure

with two separated moonpools, and a structure with two concentric annular moonpools. Yeung and Seah [3] considered the piston and other higher-order symmetric modes in a gap between two heaving rectangular cylinders. The natural frequencies and free surface modal shapes in the gap were presented together with the hydrodynamic coefficients around resonance. They pointed out that higher-order mode resonance occurred at fairly regular intervals of square of frequency. Faltinsen et al. [4] studied similar problem as in Yeung and Seah [3] by linear potential flow theory and model tests. They focused on the piston mode and found that the wave amplitudes inside the moonpool evaluated through linear theory were in general larger than that of model tests, especially for the experimental series with larger heaving amplitude. Sun et al. [5] and [6] investigated the first-order and second-order resonance behaviour between twin adjacent barges and a FLNG and LNG carrier, respectively. Zhang and Bandyk [7] studied the two-dimensional moonpool resonances for interface and surface-piercing twin rectangular bodies in a two-layer fluid. For studies in the time domain, Wang and Wu [8] gave detailed analysis of first- and second-order resonance of the liquid trapped in a gap formed by twin rectangular cylinders. They considered the two cylinders in forced vertical motion and horizontal motion in both opposite and the same directions. When first- and second-order resonance occurred, the wave runups in the gap grew quickly with time. The gap free surface elevations could be amplified up to 15 times of the heaving amplitude. Second-order resonance was clearly observed at frequencies equalling half of the natural frequencies. Wang et al. [9] further conducted firstand second-order resonance analysis of twin wedged cylinders and semielliptic cylinders in forced motions. They pointed out that when the body had flare or curvature at the water line, the resonant effect became less

E-mail address: [email protected] https://doi.org/10.1016/j.enganabound.2019.05.002 Received 4 December 2018; Received in revised form 26 March 2019; Accepted 2 May 2019 Available online 10 May 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

significant. The resonant frequency was affected by the body shape at the water line as well. Wang et al. [10] investigated the first- and secondorder gap resonance between an array of equally spaced wedge-shaped cylinders and rectangular cylinders in vertical motions. Feng and Bai [11] considered the wave resonances in the gap between side-by-side barges in beam seas. The incident waves were prescribed as the fifthorder Stokes waves with a range of wave steepness to examine the nonlinear effects on gap resonances. The term ‘response amplitude operator (RAO)’ is employed to express the amplitude ratio of the response to excitation. The gap free-surface RAOs and the resonant mode shapes were presented. They found that the resonance frequency was slightly shifted to higher values as the incident wave steepness increases. Lu et al. [12] studied the fluid resonance in narrow gaps of three identical barges fixed in incident waves. The numerical wave tank with viscous fluid flow theory was employed. Ning et al. [13] investigated the effect of the number of barges on the resonant frequency in narrow gaps. The fluid flow problem was solved based on the 2D fully-nonlinear potential flow theory. Ning et al. [14] did experimental and numerical study on wave response at the gap between two barges of different draughts. It is found that most of the researches concentrate on first-order gap resonance, while the study on second-order gap resonance is limited (Sun et al. [5] and [6]; Wang and Wu [8]; Wang et al. [9]; Wang et al. [10]). The several existing studies only confine to a fixed gap configuration. Thus, this paper aims to explore the effect of gap width on second-order gap resonance, which can be excited by forced body motions (radiation problems) and by incident waves (diffraction problems), respectively. When the radiation problems are considered, the different motion mode is studied separately. The content of this paper is organised as follows. The mathematical model based on the fully nonlinear potential flow theory and numerical procedures based on a 2D boundary element method are provided in Section 2. Section 3 gives the convergence study and validation of the present methodology. In Section 4, secondorder resonance in the gap liquid motion is studied first for heave, sway and roll mode, respectively. Then, incident wave induced second-order gap resonance is investigated. Conclusions are given in Section 5, which highlights the main contributions of the present study.

breadth of both bodies is denoted as B. The water depth is large compared to both the body dimension and wavelength and thus its effect on this problem can be neglected. A Cartesian coordinate system fixed in space is chosen such that the x-axis is along the undisturbed free surface and z-axis is directed vertically upwards. The distance between the inner sides of the two bodies is L and its middle point is taken as the origin. The initial draught of both bodies is set as D. Since the fluid motion is irrotational, a velocity potential 𝜙(x, z, t) exists which satisfies Laplace’s equation in the fluid domain: ∇2 𝜙 = 0.

(1)

On the free surface Sf , there exist both kinematic and dynamic conditions. Physically, the kinematic condition means that a fluid particle on the free surface will always stay on the free surface. Mathematically, it can be expressed as 𝜕𝜂 𝜕𝜙 𝜕𝜂 𝜕𝜙 + − =0 on 𝑆𝑓 , (2) 𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑧 where 𝜂(x, t) is the free surface elevation. The dynamic condition on the free surface requires that the water pressure on Sf is equal to the constant atmospheric pressure. Thus from Bernoulli’s equation, it becomes 𝜕𝜙 1 + |∇𝜙|2 + 𝑔𝜂 = 0 on 𝑆𝑓 . (3) 𝜕𝑡 2 On the wetted body-surfaces S0 and S1 , the impermeability condition should be satisfied because no fluid particle can penetrate or leave the rigid body’s surface. Mathematically, it yields ) ( ) ( ) 𝜕𝜙 (⃖⃖⃖⃖⃗ ⃖⃖⃗ ⃖⃖⃖⃖⃗0 ⋅ 𝑛⃗ = 𝑈0 − Ω0 𝑍0 𝑛𝑥 + 𝑉0 + Ω0 𝑋0 𝑛𝑧 = 𝑈 0 + Ω0 × 𝑋 𝜕𝑛

on 𝑆0 , (4)

( ) ( ) 𝜕𝜙 (⃖⃖⃖⃖⃗ ⃖⃖⃖⃖⃗ ⃖⃖⃖⃖⃗) = 𝑈1 + Ω1 × 𝑋1 ⋅ 𝑛⃗ = 𝑈1 − Ω1 𝑍1 𝑛𝑥 + 𝑉1 + Ω1 𝑋1 𝑛𝑧 𝜕𝑛

on 𝑆1 , (5)

where 𝜕 /𝜕 n is the partial derivative along the normal direction to the body-surface, pointing out of the fluid domain; 𝑛⃗ = (𝑛𝑥 , 𝑛𝑧 ) is the unit ⃖⃖⃖⃖⃗0 = (𝑈0 , 𝑉0 ) and 𝑈 ⃖⃖⃖⃖⃗1 = (𝑈1 , 𝑉1 ) represent the translanormal vector; 𝑈 tional velocities of each body in the x and z direction, respectively. And Ω0 , Ω1 is the rotational velocity about the rotational centre of each ⃖⃖⃖⃖⃗0 = (𝑋0 , 𝑍0 ) and 𝑋 ⃖⃖⃖⃖⃗1 = (𝑋1 , 𝑍1 ) are the position vectors in body. Also 𝑋 the body-fixed system relative to each rotational centre. The impermeability condition is applied also on the sea bottom Sb . Since the sea bottom is stationary, the condition has the form

2. Mathematical model and numerical procedures This paper mainly studies the second-order wave resonance in the gap between two floating barges. The water oscillations that take place in the gap are either excited by forced body motions or induced by incident waves. They are both covered in this study. The gap resonance problem is studied through the fully nonlinear velocity potential theory. This implies that the fluid itself is inviscid, incompressible and the fluid velocity field is irrotational.

𝜕𝜙 = 0 on 𝑆𝑏 . (6) 𝜕𝑛 A truncation boundary Sc is placed at a distance away from each body to make the fluid domain finite. A damping layer is added in front of each truncated boundary to absorb the outgoing waves generated by the motion or presence of the body. Ideally, the disturbed outgoing waves can be damped out completely by the artificial damping sponge. Thus, ‘non-reflecting boundary conditions’ can be adopted on Sc [15]. So we have 𝜕𝜙 = 0, on 𝑆𝑐 (7) 𝜕𝑛 when no incoming waves are present, or

2.1. Governing equation and boundary conditions We consider the wave motion of the liquid confined between twin floating barges, as shown in Fig. 1. For convenience, we call the body on the left hand side ‘Body-0’, and the one on the right ‘Body-1’. The

𝜕 𝜙𝐼 𝜕𝜙 = , on 𝑆𝑐 𝜕𝑛 𝜕𝑛 when there exist incoming waves, denoted by 𝜙I .

(8)

2.2. BVP involving incident waves For the diffraction problems, there are some different features. Foremost, the total velocity potential 𝜙 and free surface elevation 𝜂 can both be decomposed into two parts [16]: 𝜙 = 𝜙𝐼 + 𝜙𝐷 ,

Fig. 1. Sketch of the problem of wave interactions with two floating bodies. 2

𝜂 = 𝜂𝐼 + 𝜂𝐷 ,

(9)

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

where 𝜙I , 𝜂 I and 𝜙D , 𝜂 D denote the velocity potential and free surface elevation due to incident waves and disturbed waves, respectively. Physically, 𝜙D and 𝜂 D include the nonlinear diffraction effects of body presence on the incident waves. This split of the velocity potential and free surface elevation results in a new scheme of solving the BVP, which is to solve the disturbed velocity potential 𝜙D instead of the total velocity potential 𝜙. Firstly, the disturbed velocity potential, 𝜙D , satisfies Laplace’s equation in the fluid domain ∇2 𝜙𝐷 = 0.

potential should also be about 𝜙D . From Eq. (15) and (9), we have [19] 𝐷𝜙𝐷 𝜕 𝜙𝐷 𝜕𝜙 1 1 2 2 = + (∇𝜙 ⋅ ∇)𝜙𝐷 = − 𝐼 + ||∇𝜙𝐷 || − ||∇𝜙𝐼 || − 𝑔𝜂 on 𝑆𝑓 . 𝐷𝑡 𝜕𝑡 𝜕𝑡 2 2 (16) The updating of the free surface position is the same as in Eq. (14). It has been mentioned above that to satisfy the radiation condition, an artificial damping zone is added to minimize the reflection of the truncated boundary Sc . This is achieved by adding a damping term in Eqs. (14) and (15) artificially. So they become

(10)

𝜕𝜙 𝐷𝑥 = , 𝐷𝑡 𝜕𝑥

On the free surface Sf , the dynamic boundary condition for 𝜙D can be derived from Eqs. (3) and (9): 𝜕 𝜙𝐷 𝜕 𝜙𝐼 𝜕𝜙 𝜕 𝜙𝐼 1 2 = − = − ||∇𝜙𝐼 + ∇𝜙𝐷 || − 𝑔 𝜂 − 𝜕𝑡 𝜕𝑡 𝜕𝑡 2 𝜕𝑡

𝐷𝜙 1 = |∇𝜙|2 − 𝑔𝜂 − 𝜈(𝑥)𝜙 𝐷𝑡 2

(11)

𝜕𝜙 𝐷𝑧 = 𝐷𝑡 𝜕𝑧

on

𝑆𝑓 ,

𝐷𝜙𝐷 𝐷𝑡

(12)

𝐷𝑧 𝐷𝑡

(18)

(19)

=−

=

𝜕𝜙 𝜕𝑧

𝜕 𝜙𝐼 𝜕𝑡

2 2 + 12 ||∇𝜙𝐷 || − 12 ||∇𝜙𝐼 || − 𝑔𝜂 − 𝜈(𝑥)𝜙𝐷 ,

− 𝜈(𝑥)𝜂𝐷 ,

on 𝑆𝑓

(20)

The performance of the damping zone depends on the combination of 𝛼 and 𝛽. Quantitative study of the wave absorption efficiency is needed for each case study. For the radiation problems, the criterion for appropriate 𝛼 and 𝛽 is to make sure that the radiated outgoing waves are damped out gradually inside the damping zone and completely vanish at the truncated boundary. For the diffraction problems, the criterion is to enable the disturbed waves to be absorbed gradually with only incident wave left at the truncated boundary. A detailed example of the selection of 𝛼 and 𝛽 for a single floating body in incident waves is given by Tanizawa’s [23] study.

(13)

2.3. Hydrodynamic forces and moments calculation The hydrodynamic forces on the barges can be obtained by direct integration of the pressure over the instantaneous wetted body-surfaces S0 and S1 . The pressure in the fluid can be determined through Bernoulli’s equation ) ( 𝜕𝜙 1 𝑃 = −𝜌 + |∇𝜙|2 + 𝑔𝑧 . (21) 𝜕𝑡 2

(14)

⃖⃖⃗ and moments 𝑀 ⃖⃖⃖⃗ acting on the bodies The hydrodynamic forces 𝐹 are then obtained through

and 𝐷𝜙 1 = |∇𝜙|2 − 𝑔𝜂 𝐷𝑡 2

on 𝑆𝑓 ,

𝜔 and 𝜆 in v(x) are angular wave frequency and linear wavelength, respectively. The two nondimensional parameter 𝛼 and 𝛽 control the strength and the length of the damping zone, respectively. When there are incoming waves, the damping term should be reflected on the updating of the disturbed velocity potential 𝜙D and disturbed free surface elevation 𝜂 D as

Eqs. (10)–(13) completely define the disturbed flow involving incident waves. When dealing with this kind of wave interaction problems, there are generally two options. One can choose either to solve the BVP for the total velocity potential 𝜙 directly based on Eqs. (1)–(7) or to solve the BVP for the disturbed wave 𝜙D defined by Eqs. (10)–(13) since the incident wave is given. Both schemes are widely used in the previous researches (Sun et al. [17]; Sun et al. [18]; Zhou and Wu [19]; Zhou et al. [20]). In the present study, the former scheme is used for the radiation problems and the latter scheme is adopted for the diffraction problems. For simulations in the time domain, Lagrangian specification should be employed to track the flow quantities from time to time [21]. The kinematic and dynamic boundary conditions on the free surface can be expressed in their Lagrangian form as: 𝜕𝜙 𝐷𝑥 = , 𝐷𝑡 𝜕𝑥

(17)

in which v(x) is the damping coefficient [22] { ( ) |𝑥|−𝑥0 2 𝛼𝜔 when 𝑥0 ≤ |𝑥| ≤ 𝑥1 = 𝑥0 + 𝛽𝜆, 𝜆 𝑣(𝑥) = 0 when |𝑥| < 𝑥0 ;

Similarly on the truncation boundaries and sea bottom, the impermeability condition leads to 𝜕 𝜙𝐷 𝜕𝜙 𝜕 𝜙𝐼 = − = −∇𝜙𝐼 ⋅ 𝑛⃗ on 𝑆𝑏 , 𝑆𝑐 . 𝜕𝑛 𝜕𝑛 𝜕𝑛

on 𝑆𝑓 ,

and on 𝑆𝑓 .

It should be mentioned that the incident wave potential 𝜙I is defined in the region below the pure incident wave elevation and above the bed, which is z ≤ 𝜂 I . The total free surface elevation 𝜂 is different from 𝜂 I as a result of the disturbance of the body. Thus in Eq. (11), it may be difficult to understand the physical meaning of 𝜕 𝜙I /𝜕 t and |∇𝜙I | at positions beyond the region z ≤ 𝜂 I . But mathematically, 𝜕 𝜙I /𝜕 t and |∇𝜙I | can be readily evaluated at any position based on the expressions of Stokes wave velocity potential and free surface elevation. From Eqs. (4), (5) and (9), the boundary conditions of 𝜙D on the wetted body-surfaces S0 , S1 must satisfy ( ) 𝜕 𝜙𝐷 𝜕 𝜙𝐼 𝜕𝜙 𝜕 𝜙𝐼 ⃖⃖⃖⃗𝑖 × 𝑋 ⃖⃖⃖⃗𝑖 ⋅ 𝑛⃗ − ⃖⃖⃖⃗𝑖 + Ω = − = 𝑈 𝜕𝑛 𝜕𝑛 𝜕𝑛 𝜕𝑛 ( ) ⃖⃖⃖⃗𝑖 × 𝑋 ⃖⃖⃖⃗𝑖 − ∇𝜙𝐼 ⋅ 𝑛⃗ on 𝑆𝑖 , 𝑖 = 0, 1. ⃖⃖⃖⃗𝑖 + Ω = 𝑈

𝜕𝜙 𝐷𝑧 = − 𝜈(𝑥)𝑧 𝐷𝑡 𝜕𝑧

on 𝑆𝑓 ,

(15)

⃖⃖⃖⃗𝑖 = 𝐹

where D/Dt = 𝜕 /𝜕 t + ∇𝜙 · ∇ is the material derivative following a given fluid particle. As pointed out by Longuet-Higgins and Cokelet [21], we are able to follow the free surface elevation and velocity potential from one position on the free surface to the succeeding position using Eqs. (14) and (15). Numerically, the updating of the free surface and velocity potential can be realised by integration of the total differentials with small time steps. Since we solve the disturbed wave potential 𝜙D instead of total velocity potential 𝜙 for the diffraction problems, the updating of the velocity

∫𝑆𝑖

𝑃 𝑛⃗𝑑𝑙,

⃖⃖⃖⃖⃗𝑖 = 𝑀

∫𝑆𝑖

( ) ⃖⃖⃖⃗𝑖 × 𝑛⃗ 𝑑𝑙. 𝑖 = 0, 1 𝑃 𝑋

(22)

In Eq. (21), the time derivative of velocity potential, 𝜕 𝜙/𝜕 t, is not explicitly given even when 𝜙 has been found. If we look into the 𝜙t terms in Eq. (22), we can see that the products of 𝜙t and normal vector components are always put together as the integrands. An auxiliary function approach is proposed by Wu and Eatock Taylor [24] to circumvent the need of calculating 𝜙t so that the hydrodynamic forces can be obtained directly. One of the drawbacks of this approach is that the hydrodynamic pressure cannot be calculated because 𝜙t itself is still not known. 3

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

2.4. Numerical procedures

are unknowns and to be determined through solving Eq. (26). Applying the boundary conditions to Eq. (26) and moving all the known terms to the right hand side and the unknown terms to the left, Eq. (26) can be transformed explicitly into the matrix form as [17] [ ]{ } [ ]{( ) } 𝜙𝑛 𝑆 (𝜙)𝑆𝑁 𝐻𝑆𝑁 ,𝑆𝑁 −𝐺𝑆𝑁 ,𝑆𝑓 𝐺𝑆𝑁 ,𝑆𝑁 −𝐻𝑆𝑁 ,𝑆𝑓 ( ) 𝑁 = , 𝜙𝑛 𝑆 𝐻𝑆𝑓 ,𝑆𝑁 −𝐺𝑆𝑓 ,𝑆𝑓 𝐺 − 𝐻 𝜙) ( 𝑆 , 𝑆 𝑆 , 𝑆 𝑆𝑓 𝑓 𝑁 𝑓 𝑓 𝑓

The above boundary value problems of 𝜙 (or 𝜙D when incident waves are present) and 𝜒 are solved by the boundary element method. Take the solution of 𝜙 for example to show the numerical implementation of the BEM. Based on Green’s second identity, Laplace’s equation can be reformed as [ ] 𝜕 𝜙(𝑞 ) 𝜕 ln 𝑟𝑝𝑞 𝐴(𝑝)𝜙(𝑝) = ln 𝑟𝑝𝑞 − 𝜙(𝑞) 𝑑𝑠, (23) ∫𝜕Ω 𝜕 𝑛𝑞 𝜕 𝑛𝑞

(28) where (𝜙)𝑆𝑓 , (𝜙)𝑆𝑁 , (𝜙𝑛 )𝑆𝑓 and (𝜙𝑛 )𝑆𝑁 denote the column vectors of discrete values of 𝜙 and 𝜙n on the corresponding boundaries, respectively. In Eq. (28), 𝐻𝑆𝑁 ,𝑆𝑁 is a partition matrix of the coefficient matrix, whose elements are composed of H(m, i) and nodes m, i both belong to SN . The other seven partition matrices of the coefficient matrices in Eq. (28) are constructed similarly. The accuracy of the solution of the above matrix equation determines the accuracy of the numerical results for a certain mesh. For time domain simulations, the fluid velocity at the free surface and wetted body nodes should be determined in order to update the information from one time to next. Since 𝜙 and 𝜙n are known at these nodes, it is straightforward to compute the normal and tangential velocity at the mid-point of each element through linear interpolation. Let v𝜏 and vn denote tangential and normal components of the fluid velocity at mid-points, we have

where rpq is the distance between points p and q, A(p) is the solid angle at point p. And 𝜕Ω is the boundary of the fluid domain described anticlockwise. This integral representation suggests that the unknown function 𝜙 at any point can be calculated, if all the values of 𝜙 and 𝜕 𝜙/𝜕 n on the boundary 𝜕Ω are known. Eq. (23) is usually referred to as the boundary integral equation when p has been moved to the boundary. The solid angle A(p) in Eq. (23) depends on the relative position of point p and the fluid boundary. We first discretise the whole fluid boundary 𝜕Ω into a large number N of straight-line elements. Along each element, a linear interpolation function is adopted to approximate the velocity potential 𝜙 and its normal derivative 𝜙n between two nodes. Introducing a local coordinate 𝜉 along element i, denoted as 𝜕Ωi , i = 1, 2,…, N, the velocity potential 𝜙 and its normal derivative 𝜙n on the whole discretised boundary can be written in the form 𝜙=

𝑁+1 ∑ 𝑖=1

𝑓𝑖 (𝜉, 𝑥⃗)𝜙𝑖 ,

𝜙𝑛 =

𝑁+1 ∑ 𝑖=1

𝑓𝑖 (𝜉, 𝑥⃗)𝜙𝑛,𝑖 ,

𝑣𝜏,𝑖 = (𝜙𝑖+1 − 𝜙𝑖 )∕𝐿𝑖 ,

⎧(1 − 𝜉)∕2, 𝑥⃗ ∈ 𝜕 Ω 𝑖 ⎪ 𝑓𝑖 (𝜉, 𝑥⃗) = ⎨(1 + 𝜉)∕2, 𝑥⃗ ∈ 𝜕 Ω𝑖−1 ⎪0.𝑥⃗ ∉ 𝜕 Ω𝑖 , 𝜕 Ω𝑖−1 ⎩

𝑢mid,𝑖 = 𝑣𝑛,𝑖 𝑛𝑥 − 𝑣𝜏,𝑖 𝑛𝑧 ,

𝐴(𝑗 )𝜙𝑗 =

𝑖=1

∫𝜕 Ω𝑖

] 𝜕 ln 𝑟 𝜙𝑛 ln 𝑟 − 𝜙 𝑑 𝑠. 𝜕𝑛

Finally, substituting Eq. (24) into the discretised integral Eq. (25) and rearranging the expression according to velocity potential and its normal derivative, we obtain (26)

𝑖=1

in which 𝐺(𝑗, 𝑖) =

∫𝜕Ω

𝑓𝑖 (𝑥⃗, 𝜉) ln 𝑟𝑑𝑠(𝑥⃗) = 𝐺2 (𝑗, 𝑖 − 1) + 𝐺1 (𝑗, 𝑖);

⎧ 𝑖+1 = 𝑌 𝑖 + 𝐾2 Δ𝑡, ⎪𝑌 ⎪𝐾 = 𝑓 (𝑥𝑖 , 𝑧𝑖 , 𝜙𝑖 ), ⎨ 1 ( ⎪ Δ𝑡 𝑖 𝑖 ⎪𝐾 2 = 𝑓 𝑥 + 2 𝐾 1 , 𝑧 + ⎩

𝜕 ln 𝑟 | 𝐻(𝑗, 𝑖)|𝑗≠𝑖 = 𝑓 (𝑥⃗, 𝜉) 𝑑𝑠(𝑥⃗) = 𝐻2 (𝑗, 𝑖 − 1) + 𝐻1 (𝑗, 𝑖); | ∫𝜕Ω 𝑖 𝜕𝑛 𝐺1 (𝑗, 𝑖) = 𝐻1 (𝑗, 𝑖) =

∫𝜕 Ω𝑖

(1 − 𝜉) (1 + 𝜉) ln 𝑟𝑑𝑠; 𝐺2 (𝑗, 𝑖 − 1) = ln 𝑟𝑑𝑠; ∫𝜕 Ω𝑖−1 2 2

∫𝜕 Ω𝑖

(1 − 𝜉) 2

𝜕 ln 𝑟 𝑑𝑠; 𝐻2 (𝑗, 𝑖 𝜕𝑛

− 1) =

∫𝜕 Ω𝑖−1

(30)

Special treatment is needed for the intersection points, which belong to both free surface and wetted body-surface. The velocities at these points are essential because they decide the accuracy of the wave runup on the body and affect the evaluation of the hydrodynamic loads on structures. Thus, a four-point Lagrangian interpolation, which gives higher accuracy, is used to calculate the velocities at intersection points instead of linear interpolation. Having obtained the fluid velocity at free surface nodes, a secondorder Runge-Kutta time integration method is adopted to update the free-surface position and the corresponding velocity potential on it. Let Y = (x, z, 𝜙) be expressed by a function DY/Dt = f(x, z, 𝜙), then second order Runge-Kutta scheme gives

(25)

𝑁 ∑ ( ) 𝐺(𝑗 , 𝑖)𝜙𝑛,𝑖 − 𝐻(𝑗 , 𝑖)𝜙𝑖 = 0, 𝑗 = 1, 2, … , 𝑁

𝑣mid,𝑖 = 𝑣𝑛,𝑖 𝑛𝑧 + 𝑣𝜏,𝑖 𝑛𝑥 .

Finally, the Cartesian components (u, v) of the fluid velocity at node i are calculated through the interpolation ( ) ( ) 𝑢𝑖 = 𝑢mid,𝑖 𝐿𝑖−1 + 𝑢𝑚𝑖𝑑,𝑖−1 𝐿𝑖 ∕ 𝐿𝑖 + 𝐿𝑖−1 , ( ) ( ) (31) 𝑣𝑖 = 𝑣mid,𝑖 𝐿𝑖−1 + 𝑣𝑚𝑖𝑑,𝑖−1 𝐿𝑖 ∕ 𝐿𝑖 + 𝐿𝑖−1 .

and 𝑥⃗ = (𝑥, 𝑧) is a vector in the global Cartesian coordinate system, representing the position of point 𝜉. Then, the boundary integral Eq. (23) should be discretised according to the newly discretised boundary, which has the form [

(29)

Then the horizontal and vertical components of velocity at the midpoint of element i can be calculated by vector decomposition of the normal and tangential velocity v𝜏,i and vn,i . Hence the velocity components umid,i and vmid,i can be expressed as

(24)

where

𝑁 ∑

𝑣𝑛,𝑖 = (𝜙𝑛,𝑖+1 + 𝜙𝑛,𝑖 )∕2.

Δ𝑡 𝐾 , 𝜙𝑖 2 1

+

Δ𝑡 𝐾 2 1

)

(32) ,

where Yi = (xi ,zi ,𝜙i ) are known vectors corresponding to time t at the ith step and Yi + 1 are unknown terms at time t + Δt with step i + 1. For long time simulations, some over-distortion of elements may occur due to accumulated numerical error. Thus, special treatments, such as remeshing, smoothing, jet cutting and so on, should be applied to make the simulation run successfully (refer to [25]) for more details about how these treatments are performed).

(1 + 𝜉) 𝜕 ln 𝑟 𝑑𝑠. (27) 𝜕𝑛 2

The integrations over the element 𝜕Ωi for G(j, i) and H(j, i) can be performed analytically. There are usually three types of boundary conditions: Dirichlet condition, Neumann condition and mixed conditions. The present BVP for 𝜙 has mixed boundary conditions because on Sf the velocity potential itself is known and on SN its normal derivative is known. Here, SN stands for the whole boundary, excluding the free surface for convenience. The subscript N of SN means that Neumann condition is applied on the corresponding boundaries. The normal derivative of 𝜙 on Sf and 𝜙 itself on SN

3. Convergence study and comparison In this section, we verify the convergence and accuracy of the numerical method. The twin barges undergoing forced heaving oscillation 4

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 2. Mesh convergence study. (a) Wave runup on the left side of Body-1; (b) vertical hydrodynamic forces on Body-1.

Fig. 3. Time convergence study. (a) Wave runup on the left side of Body-1; (b) vertical hydrodynamic forces on Body-1.

are simulated. The barges are subjected to the following motion 𝑠𝑧 (𝑡) = 𝐴 sin 𝜔𝑧 𝑡

the results related to Body-1 are compared. Note that Fz excludes the contribution from the initial buoyancy forces. It is observed that both the wave runup and vertical force results from the latter two meshes agree well with each other. The mesh with the element size Δd = 𝜆/114 can ensure the convergence of the result, which will be adopted in the following simulations. The convergence of the time step is also studied. Three different time steps Δt = T/60, T/120 and T/240 are applied. Fig. 3 gives the comparison of time histories of wave runup on the left side of Body-1 and vertical hydrodynamic force on Body-1. Results from Δt = T/120 and T/240 have shown an excellent agreement. The comparison suggests that the time step Δt ≤ T/120 can guarantee the time convergence of results. The time step Δt = T/120 is chosen for future simulations. As the accuracy study, the added-mass in heave and heave damping coefficients obtained from the present numerical method are compared with the experimental and theoretical results in Lee et al. [26]. The heave added-mass coefficient 𝜇̄ and heave damping coefficient 𝜆 are determined from the force history F(t) through

(33)

where sz , A and 𝜔z are the displacement, amplitude and angular frequency of the prescribed vertical motion, respectively. In order to compare with the results in Lee et al. [26], the fluid density 𝜌, gravitational acceleration g and half breadth B/2 are chosen to perform the nondimensionalisation. Then, √ √ the time and frequency can be scaled as 𝑡 → 𝑡 𝑔 ∕(𝐵 ∕2) and 𝜔 → 𝜔𝑧 (𝐵∕2)∕𝑔 . The case of B = 2, D = 1, L = 6 and A = 1/24, 𝜔2 = 1.1 is chosen to do the convergence study. We could know the heave period T ≈ 6 and the wavelength 𝜆 ≈ 5.7 from the dispersion relation 𝜔2 = 2𝜋∕𝜆.

(34)

The truncated boundaries are placed at 5𝜆 away from the body on each side. Damping parameters are chosen as 𝛼 = 0.5 and 𝛽 = 2. Initially, all elements on the near region of the free surface and wetted bodysurface are distributed equally with size Δd. On the truncated boundaries and sea bottom, large elements are used because they are far away from the bodies. And they are not included in the mesh convergence study. During the simulation, the wetted body-surface elements always have the same size Δd at every time step. Three meshes with the element size Δd = 𝜆/57, 𝜆/114 and 𝜆/228 are considered. The time step is chosen temporarily as Δt = T/120. The comparison of the corresponding results is shown in Fig. 2. Since the problem is symmetric about the z-axis, only

𝜇̄ =

1 4𝐴𝜔2 ∫𝑡

𝜆=

1 4𝐴𝜔2 ∫𝑡

𝑡+𝑇

𝑡+𝑇

𝐹 (𝑡) sin 𝜔𝑡𝑑𝑡

(35)

𝐹 (𝑡) cos 𝜔𝑡𝑑𝑡

(36)

Here, 4 in the denominator is the non-dimensional mass of the two barges. The comparison of heave added-mass and damping coefficients 5

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 4. Comparison of present results with previous results. (a) Heave added-mass coefficients; (b) heave damping coefficients.

Fig. 5. Wave elevation history in the gap along Body-1 for L = 1, 𝜔z = 0.5𝜔0 , A = 0.05 and its corresponding FFT analysis.

Fig. 6. Wave elevation history in the gap along Body-1 for L = 1, 𝜔z = 0.5𝜔0 , A = 0.2 and its corresponding FFT analysis.

4. Numerical results and discussions

can be found in Fig. 4(a) and (b), respectively. It can be seen that, for both cases, the present results agree well with the theoretical ones. Minor discrepancies between the numerical and experimental results are also observed, which might due to physical factors such as the fluid viscosity. In general, the present numerical method could provide reliable results.

Previous studies have shown that a classic resonance or first-order resonance can occur when the external excitation frequency is at one of the natural frequencies of the liquid confined in the gap. One would speculate that even if the external excitation frequency is not at a 6

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 7. FFT analysis of wave elevations along the left side of Body-1 for L = 2. (a) 𝜔z = 0.5𝜔0 ; (b) 𝜔z = 0.5𝜔2 .

Fig. 8. Wave elevation history in the gap along Body-1 for L = 4, 𝜔z = 0.5𝜔2 , A = 0.2 and its corresponding FFT analysis.

Fig. 9. Wave elevation history in the gap along Body-1 for L = 8 , 𝜔z = 0.5𝜔2 , A = 0.2 and its corresponding FFT analysis.

natural frequency in the gap, second-order resonance may occur if its double frequency equals one of them. This speculation has been proved mathematically by Wu [27] for the sloshing in a 2D rectangular tank. Second-order resonance is related to second-order effects, in which the amplitudes of the second-order terms (instead of the first-order terms) grow with time. Wang and Wu [8] numerically confirmed that this dis-

tinct type of resonance exists in the liquid confined in the gap. The physical parameters are nondimensionalised using the density of the fluid 𝜌, the body width B and the acceleration due to gravity g. The draught of the bodies is set at D = 1 for all the following cases. The forces presented below are normalised by mA, where the non-dimensional mass m of each body equals the initial mass displacement. 7

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 10. Wave elevation results for L = 4, 𝜔x = 𝜔1 . (a) Wave profiles in the gap; (b) wave elevation history in the gap along Body-1 and its corresponding power spectrum analysis.

Fig. 11. Wave elevation results for L = 4, 𝜔x = 𝜔2 . (a) Wave profiles in the gap; (b) wave elevation history in the gap along Body-1 and its corresponding power spectrum analysis.

Fig. 12. Wave elevation results for L = 4, 𝜔x = 𝜔3 . (a) Wave profiles in the gap; (b) wave elevation history in the gap along Body-1 and its corresponding power spectrum analysis.

8

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 13. Wave elevation history in the gap along Body-1 for L = 2, 𝜔x = 0.5𝜔1 and its corresponding FFT analysis.

Fig. 14. Wave elevation history at the left side of Body-1 for L = 2, 𝜔x = 0.5𝜔2 and its corresponding FFT analysis.

Fig. 15. Wave elevation histories at the left side of Body-1 for L = 4, 𝜔x = 0.5𝜔2 and its corresponding FFT analysis.

configurations can be calculated through the disturbance approach [28].

4.1. Second-order gap resonance excited by forced body motions Second-order resonance behaviour in the gap liquid motion is studied, which is excited by heave, sway and roll mode of forced body motions, respectively. The forces on the surrounding bodies are calculated. The lowest several natural frequencies of different gap

4.1.1. Heave motions The twin bodies are excited harmonically as sz = Asin (𝜔z t). As suggested by Wang and Wu [8], the typical condition of the second-order 9

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 16. The hydrodynamic forces on Body-1 and their corresponding FFT analysis when L = 4, 𝜔x = 0.5𝜔2 , A = 0.25.

Fig. 17. Wave elevation histories at the left side of Body-1 for L = 8, 𝜔x = 0.5𝜔2 and its corresponding FFT analysis.

resonance is 𝜔z = 0.5𝜔2i , i = 0, 1, 2,…. The situation of 𝜔z = 0.5𝜔0 for L = 1 is considered first. The wave elevation history along the left side of Body-1 for A = 0.05 is shown in Fig. 5(a). The FFT analysis is performed for two time intervals of the elevation history, as presented in Fig. 5(b). The amplitude at the frequency component 𝜔0 in this case reflects the level of second-order effect. The increase of its amplitude with time will suggest the occurrence of second-order resonance, according to its definition. The wave elevation history shows the second-order effect is minor. Meanwhile, the increase of the amplitude at 𝜔0 is

hardly observable. We now amplify the amplitude A to 0.2 to see the second-order effect more clearly, because the second-order terms are proportional to A2 . The wave elevation for A = 0.2 is provided in Fig. 6(a), which shows the second-order effect is evident. However, the FFT analysis performed for three time intervals still demonstrates very little increase of it, see Fig. 6(b). This may suggest the second-order resonance will not be triggered easily in the case of L = 1. The cases of L = 2 are considered next. The FFT analysis of wave elevations along the left side of Body-1 for both 𝜔z = 0.5𝜔0 and 𝜔z = 0.5𝜔2 10

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 18. The hydrodynamic forces on Body-1 and their corresponding FFT analysis when L = 8, 𝜔x = 0.5𝜔2 , A = 0.25.

Fig. 19. Wave elevation histories at the left side of Body-1 for L = 2, 𝜔r = 𝜔2 and its corresponding FFT analysis.

is presented in Fig. 7(a) and (b), respectively. The component of 𝜔0 is small and its amplitude is almost unchanged with time. The increase of the amplitude at 𝜔2 is subtle, although the second-order effect is pronounced at 𝜔z = 0.5𝜔2 . Now we consider the case of larger gap width L = 4. Its wave elevation history along the left side of Body-1 is given in Fig. 8(a). The nonlinear effects are obvious with complex components, as detailed in Fig. 8(b). The amplitude at the double frequency component is larger than that at the excitation frequency, and it clearly increases with time. Second-order resonance is evident in this case.

The above analysis suggests that the second-order resonance phenomenon is largely dependent on the gap width over depth ratio. The gap depth is defined as the initial body draught, regardless of the wave elevations in the gap during oscillations. Specifically, the bigger the value of the ratio is, the more evident the second-order resonance becomes. To test if this speculation is correct, an even wider gap, L = 8 is further simulated. The wave elevation along the left side of Body1 is presented in Fig. 9(a). Its shape reveals stronger second-order effect. The corresponding FFT analysis performed on three time intervals, which is given in Fig. 9(b), shows clearly the occurrence of second-order 11

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 20. Wave elevation histories at the left side of Body-1 for L = 2, 𝜔r = 0.5𝜔2 and its corresponding FFT analysis.

Fig. 21. Wave elevation histories at the left side of Body-1 for L = 4, 𝜔r = 0.5𝜔2 and its corresponding FFT analysis.

Fig. 22. Wave elevation histories at the left side of Body-1 for L = 4, 𝜔r = 0.5𝜔3 and its corresponding FFT analysis.

resonance. The amplitude at 𝜔2 is even several times larger than that at 0.5𝜔2 . Therefore, it is correct to say that the second-order resonance is more pronounced when the gap is wider.

for L = 4 have been calculated through the disturbance approach, namely 𝜔1 = 0.9188, 𝜔2 = 1.2531, 𝜔3 = 1.5038. Before investigating the second-order resonance, the first-order resonance is studied first for comparison. The simulations are run at each of the three natural frequencies. The wave elevation results are provided in Figs. 10–12, respectively. For the case of 𝜔x = 𝜔1 , the wave profiles in the gap show that it is indeed the first sloshing-type mode. Standing waves are formed with a node at the middle of the gap and two anti-nodes at

4.1.2. Sway motions The twin bodies are undergoing the same horizontal motion, which we prescribe as sx = Asin (𝜔x t). The first three natural frequencies 12

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 23. Wave elevation histories at the left side of Body-1 for L = 8, 𝜔r = 0.5𝜔2 and its corresponding FFT analysis.

Fig. 24. Wave elevation along the left side of Body-1 and its corresponding FFT analysis for L = 1, 𝜔I = 0.5𝜔0 .

the gap walls. The free surface elevation along the left side of Body-1 demonstrates a clearly increasing trend for the first 50 periods. It then stops growing due to radiation damping in the gap. The corresponding power spectrum analysis on the wave elevation also shows that there is only one frequency dominating the liquid motion in the gap, which coincides with the excitation frequency. This confirms that first-order resonance has taken place when 𝜔x = 𝜔1 . For the case of 𝜔x = 𝜔2 , the free surface profiles in the gap demonstrate that the anti-nodes move 𝜆/4 towards the centre of the gap, which suggests no resonance has occurred. Additionally, the first-order resonance cannot be observed from the wave elevation history in Fig. 11(b). This particular irregular shape of the wave elevation implies multiple dominant frequencies. The corresponding power spectrum analysis verifies this. The first natural frequency 𝜔1 is triggered and the power spectral density of it overtakes that of the excitation frequency 𝜔2 . As a matter of fact, 𝜔2 is associated with a symmetric mode, while the sway motion of the two bodies in the same direction results in an anti-symmetric motion of the liquid. The first-order resonance has also occurred in the case of 𝜔x = 𝜔3 . The wave elevation history and power spectral distribution are similar to that of 𝜔x = 𝜔1 . The wave profiles show it is the third sloshing-type mode based on the observation that there are standing waves formed in the gap with three nodes and two anti-nodes at the gap walls. It is interesting to notice that first-order resonance does not occur at 𝜔x = 𝜔2 . This is because the flow induced by the same horizon-

tal motion of twin bodies is anti-symmetric. Thus, the first-order resonance is expected at an odd mode, or rather 𝜔x = 𝜔2i − 1 , i = 1, 2,…. Furthermore, Wu [27] has proved this mathematically based on liquid sloshing in a tank. In terms of the second-order resonance, he has pointed out that it will not be excited when the double frequency is equal to one of the natural frequencies of an odd mode. When the double frequency equals one of the even modes, such resonance does occur. Fig. 13 gives the wave elevation history in the gap along Body-1 and its FFT analysis for three time intervals for L = 2, 𝜔x = 0.5𝜔1 . The amplitude at its double frequency 𝜔1 is large. However, it decreases with time, which means no second-order resonance in the gap liquid motion. The FFT analysis of wave elevation history along the same point for L = 2, 𝜔x = 0.5𝜔2 , which is given in Fig. 14(b), has shown a different pattern. It is found that the amplitude at its double frequency 𝜔2 is in general small but does not drop with time. It is hard to tell if the second-order resonance occurs. Considering that the second-order resonance associated with heave motion is more easily triggered in wide gaps, we consequently consider the case of L = 4, 𝜔x = 0.5𝜔2 . The wave elevation along the left side of body-1 is given in Fig. 15(a). Its FFT analysis is provided in Fig. 15(b). Three dominant frequencies are observed in the gap motion, namely the excitation frequency 0.5𝜔2 , first natural frequency 𝜔1 and double frequency 𝜔2 . The amplitude at 𝜔1 clearly decreases with time, while the amplitude at 𝜔2 slightly increases with time. It suggests the occurrence of second-order resonance. The 13

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 25. Horizontal forces on both bodies and their corresponding FFT analysis for L = 1, 𝜔I = 0.5𝜔0 . (a) and (b) Upwave body Body-0; (c) and (d) leeside body Body-1.

Fig. 26. Wave elevation on the right side of the gap and corresponding FFT analysis for L = 4, 𝜔I = 0.5𝜔1 .

hydrodynamic forces on Body-1 and their corresponding FFT analysis are presented in Fig. 16. The second-order resonance has been reflected in the forces. In order to observe the second-order resonance more clearly, we consider the case of even larger gap width L = 8. The FFT analysis of wave elevation along the left side of Body-1 is provided in Fig. 17(b). Evidently, the amplitude at the double frequency grows with time in this case. It exceeds the amplitude at the excitation frequency after a certain

period. The hydrodynamic forces on Body-1 are given in Fig. 18. The same trend as wave elevation is shown. In a word, the second-order resonance is more easily provoked in a wider gap in sway motion mode. This is consistent with the observations in horizontal liquid sloshing in a tank [29,30]. Regarding twin barges undergoing horizontal motion in opposite direction prescribed as s0x = Asin (𝜔x t), s1x = −Asin (𝜔x t), the flow will be symmetric as in the vertical motion. As pointed out by Wang and Wu 14

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 27. Horizontal forces on both bodies and their corresponding FFT analysis for L = 4, 𝜔I = 0.5𝜔1 . (a) and (b) Upwave body Body-0; (c) and (d) leeside body Body-1.

Fig. 28. Wave elevation on the right side of the gap and corresponding FFT analysis for L = 4, 𝜔I = 0.5𝜔2 .

[8], the resonant behaviour is expected to be similar to that in the vertical motion. Specifically, the first-order resonance will take place at natural frequency of an even mode and the second-order resonance at half of the natural frequency of an even mode.

rotating about its own mass centre, which is located at half of the draught below the still water. The case of L = 2 is simulated. The wave runup on the left side of body-1 for 𝜔r = 𝜔2 is plotted in Fig. 19(a). Its irregular form implies multiple frequencies existing in the system, and the general decreasing trend indicates that no first-order resonance takes place. The corresponding FFT analysis of wave elevation confirms there are three dominant frequencies. The amplitude at the excitation frequency 𝜔2 does not increase with time. The case of 𝜔r = 0.5𝜔2 is then tested to see if there exists second-order resonance. It is observed in

4.1.3. Roll motions The angular displacements of the twin barges are both expressed as 𝜃 = 𝜃 0 sin (𝜔r t). The positive direction is anticlockwise. Each body is 15

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 29. Horizontal forces on both bodies and their corresponding FFT analysis for L = 4, 𝜔I = 0.5𝜔2 . (a) and (b) Upwave body Body-0; (c) and (d) leeside body Body-1.

Fig. 20(a) that the wave elevation is nearly linear. Its FFT analysis shows directly that the second-order effect is too small in this case. The cases of L = 4 are simulated next. For 𝜔r = 0.5𝜔2 , the wave elevation in the gap along Body-1 is presented in Fig. 21. It can be seen clearly that the second-order effect is strong. Moreover, the amplitude at its double frequency grows obviously with time, which proves the occurrence of second-order resonance. For 𝜔r = 0.5𝜔3 , the amplitude of wave elevation at double frequency can be neglected as shown in Fig. 22(b). The dominant frequencies in the system are excitation frequency 0.5𝜔3 and the first natural frequency 𝜔1 . No second-order resonance can be observed in this case. The case of larger gap width L = 8 is simulated. Clearer second-order resonance can be observed in Fig. 23, which gives the wave elevation results for 𝜔r = 0.5𝜔2 . The wave elevation history in the gap demonstrates an increasing trend because of the resonance. The amplitude at its double frequency is several times larger than that at the excitation frequency. If we compare the above cases of roll motion carefully, the following conclusions can be drawn. Firstly, the first-order resonance cannot occur at the natural frequency of an even mode. Secondly, the second-order resonance can take place at half of the natural frequency of an even mode. Thirdly, for the same excitation amplitude, second-order resonance is more pronounced in a wider gap. Therefore, the flow produced by roll motion resembles that of sway motion in terms of the feature of first- and second-order resonance.

4.2. Second-order gap resonance induced by incident waves In terms of the resonance effects on the surrounding bodies, there should be some differences between radiation problems and diffraction problems. In the previous radiation problems, the flow is either symmetric or anti-symmetric, while in the diffraction problems, the flow is asymmetric. Moreover, the two bodies are distinguished by upwave body and leeside body depending on the propagating direction of the incident wave. This subsection will study the incident wave induced second-order resonance in the gap formed by two fixed bodies. Second-order resonance condition is the same as in the radiation problems, which gives that double frequency of the incident wave 𝜔I is equal to one of the natural frequencies. We take the fifth-order Stokes wave as the incoming wave, the corresponding incident velocity potential and the wave elevation are provided as [31]: √ 𝜙𝐼 =

(( ) 1 1 37 5 𝑘𝑧 1 𝜀 − 𝜀3 − 𝜀 𝑒 sin 𝜑 + 𝜀4 𝑒2𝑘𝑧 sin 2𝜑 2 24 2 𝑘3 ) 1 5 3𝑘𝑧 + 𝜀 𝑒 sin 3𝜑 , 12

) ( ) 3 2 211 4 1 1 𝜀 − 𝜀 cos 𝜑 + 𝜀 + 𝜀3 cos 2𝜑 8 192 2 3 ( ) ) 3 2 99 4 1 3 125 4 + 𝜀 + 𝜀 cos 3𝜑 + 𝜀 cos 4𝜑 + 𝜀 cos 5𝜑 , 8 128 3 384

𝜂𝐼 = 𝐴

16

(37)

((

1−

(38)

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

Fig. 30. Wave elevation and horizontal forces on both bodies and their corresponding FFT analysis for L = 8, 𝜔I = 0.5𝜔2 . (a) and (b) Wave elevation on the right side of the gap; (c) and (d) upwave body Body-0; (e) and (f) leeside body Body-1.

For numerical time domain simulations of prescribed incident waves, it is possible to allow a gradual development of the disturbed potential. This is achieved by multiplying the incident wave potential 𝜙I and free surface elevation 𝜂 I by a modulation function Mt . The expression of Mt could be written as [32] ] { [ 0.5 1 − cos (𝜋𝑡∕2𝑇 ) , 𝑡 < 2𝑇 𝑀𝑡 = (41) 1, 𝑡 ≥ 2𝑇

where 𝜑 = 𝑘𝑥 − 𝜔𝐼 𝑡, 𝜔𝐼 =

) √ ( 1 1 𝑘 1 + 𝜀2 + 𝜀4 2 8

(39)

(40)

where k = 2𝜋/𝜆 is its corresponding wave number; 𝜔I is the incident wave frequency; A = H/2, where H is the peak to trough incident wave height and 𝜀 = kH/2 = kA. The parameter 𝜀 is known as the wave steepness. The incident wave speed cI = 𝜔I /k increases as 𝜀 increases because of the nonlinear wave relation.

where T = 2𝜋/𝜔I is the period of the incident wave. For narrower gap L = 1, the case of incident wave with 𝜔I = 0.5𝜔0 is considered. The wave steepness is set as 𝜀 = 0.0566 to reveal the 17

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

second-order effect clearer. Fig. 24 gives the wave elevation results. The elevation becomes periodic after about 20 periods. The FFT analysis on the elevation history over three time intervals shows that the amplitude at its double frequency increases slightly with time before the periodic state is reached. It ought to be pointed out that the first 10 periods are excluded from the FFT analysis due to the adoption of modulation function. This rule is applied for all the following FFT analysis. The wave elevation in the gap is much smaller than that of the first-order resonance when 𝜔I = 𝜔0 , where the gap elevation RAO is nearly 6 [28]. In terms of the hydrodynamic forces acting on both bodies when attacked by incident waves, horizontal forces are mainly due to the pressure difference of the two side walls of each body and are therefore considerably affected by the liquid motion in the gap. Vertical forces, however, arise from the pressure distribution on the bottom of each body and may be dominated by the inertia terms. Thus, we focus more on the horizontal forces on the bodies at such resonance. They are given in Fig. 25. The force histories exhibit a strong second-order effect, especially on the leeside body. The amplitude at its double frequency is in fact larger than that at the excitation frequency. The horizontal force on the upwave body is larger than that on the leeside body. Second-order resonance reveals itself more clearly on the horizontal forces. For a wide gap L = 4, the cases of 𝜔I = 0.5𝜔1 and 𝜔I = 0.5𝜔2 are both simulated. The incident wave steepness is 𝜀 = 0.0283. The wave runup on Body-1 in the gap and the corresponding FFT analysis for three different time intervals are presented in Fig. 26 for 𝜔I = 0.5𝜔1 . The increasing trend with time is obvious due to the rapid growth of the amplitude at the double frequency. The amplitude at the excitation frequency 𝜔I remains unchanged. It is expected that the second-order effect of wave elevation in the gap would eventually exceed the first-order effect. Fig. 27 gives the horizontal forces on both bodies and their corresponding FFT analysis at resonance for 𝜔I = 0.5𝜔1 . The horizontal forces become larger and larger as time grows due to the increasing wave elevations in the gap. The wave elevation and horizontal force results of 𝜔I = 0.5𝜔2 are presented in Figs. 28 and 29, respectively. The wave elevation history is complicated because the first three natural modes have all been excited in the gap motion. As time goes on, the amplitude of the second mode increases, while the amplitudes of the first and third modes decrease. The horizontal force on the upwave body is in general larger than that on the leeside body. The second-order effects on the horizontal forces are not that pronounced even though the second-order resonance takes place in the gap motion. For an even wider gap L = 8, the case of 𝜔I = 0.5𝜔2 is tested with the incident wave steepness 𝜀 = 0.0283. The results are presented in Fig. 30. Comparing the FFT analysis of wave elevations on the right side of the gap for L = 4 and L = 8 with 𝜔I = 0.5𝜔2 , see Figs. 28(b) and 30(b), it is clear that second-order resonance is more pronounced for L = 8. This statement holds water for horizontal forces on both barges.

gap and periodic states can be reached. Second-order resonance is more pronounced and easily provoked in a larger gap. It suggests that for sideby-side vessels of small draughts second-order resonance should also be considered during operation. For wave resonance in the gap induced by an incident wave, there are some features different from the above radiation problems. Firstly, the total flow field is now always asymmetric. First-order resonance can be triggered when the incoming wave frequency 𝜔I is at any of the natural frequencies, and second-order resonance when 𝜔I is at half of any natural frequencies. Second-order resonance in the gap can manifest themselves in the wave loads on both bodies. It is found that second-order resonance is more easily provoked and especially significant when the gap width over draught ratio is large. However, it ought to be pointed out that the present work is based on the potential flow theory without considering the fluid viscosity and energy dissipations. The predictions of amplitudes of free surface elevations and hydrodynamic forces at resonant frequencies are not reliable. The viscous flow model will be adopted in future study for accurately predicting the free surfaces elevations and hydrodynamic forces. Acknowledgements I am truly grateful to the Erasmus Mundus Europe Asia (EMEA) scholarship programme, which supported me financially during the early period of this study. I also wish to acknowledge the financial support by Senior Talent Foundation of Jiangsu University (Grant No. 18JDG036). Special thanks go to Prof. Guo Xiong Wu from University College London for his valuable discussions. References [1] Molin B. On the piston and sloshing modes in moonpools. J Fluid Mech 2001;430:27– 50. doi:10.1017/S0022112000002871. [2] Newman JN. Low-frequency resonance of moonpools. In: Proceedings of the eigteenth international workshop on water waves and floating bodies; 2003. p. 9–12. [3] Yeung RW, Seah RKM. On Helmholtz and higher-order resonance of twin floating bodies. J Eng Math 2007;58:251–65. doi:10.1007/s10665-006-9109-3. [4] Faltinsen OM, Rognebakke OF, Timokha AN. Two-dimensional resonant piston-like sloshing in a moonpool. J Fluid Mech 2007;575:359–97. doi:10.1017/S002211200600440X. [5] Sun L, Eatock Taylor R, Taylor PH. First- and second-order analysis of resonant waves between adjacent barges. J Fluids Struct 2010;26:954–78. doi:10.1016/j.jfluidstructs.2010.06.001. [6] Sun L, Eatock Taylor R, Taylor PH. Wave driven free surface motion in the gap between a tanker and an FLNG barge. Appl Ocean Res 2015;51:331–49. doi:10.1016/j.apor.2015.01.011. [7] Zhang X, Bandyk P. Two-dimensional moonpool resonances for interface and surface-piercing twin bodies in a two-layer fluid. Appl Ocean Res 2014;47:204–18. doi:10.1016/j.apor.2014.05.005. [8] Wang CZ, Wu GX. Analysis of second-order resonance in wave interactions with floating bodies through a finite-element method. Ocean Eng 2008;35:717–26. doi:10.1016/j.oceaneng.2008.02.004. [9] Wang CZ, Wu GX, Khoo BC. Fully nonlinear simulation of resonant motion of liquid confined between floating structures. Comput Fluids 2011;44:89–101. doi:10.1016/j.compfluid.2010.12.014. [10] Wang CZ, Meng QC, Huang HC, Khoo BC. Finite element analysis of nonlinear wave resonance by multiple cylinders in vertical motions. Comput Fluids 2013;88:557–68. doi:10.1016/j.compfluid.2013.10.012. [11] Feng X, Bai W. Wave resonances in a narrow gap between two barges using fully nonlinear numerical simulation. Appl Ocean Res 2015;50:119–29. doi:10.1016/j.apor.2015.01.003. [12] Lu L, Cheng L, Teng B, Zhao M. Numerical investigation of fluid resonance in two narrow gaps of three identical rectangular structures. Appl Ocean Res 2010;32:177– 90. doi:10.1016/j.apor.2009.10.003. [13] Ning D, Su X, Zhao M, Teng B. Numerical study of resonance induced by wave action on multiple rectangular boxes with narrow gaps. Acta Oceanol Sin 2015;34:92–102. doi:10.1007/s13131-015-0672-1. [14] Ning D, Zhu Y, Zhang C, Zhao M. Experimental and numerical study on wave response at the gap between two barges of different draughts. Appl Ocean Res 2018;77:14–25. doi:10.1016/j.apor.2018.05.010. [15] Givoli D. Non-reflecting Boundary Conditions. J Comput Phys 1991;94:1–29. doi:10.1016/0021-9991(91)90135-8. [16] Ferrant P. Simulation of strongly nonlinear wave generation and wave-body interactions using a 3D MEL model. In: Proceedings of the twenty-first symposium on naval hydrodynamics; 1996. p. 93–108. [17] Sun SY, Sun SL, Wu GX. Oblique water entry of a wedge into waves with gravity effect. J Fluids Struct 2015;52:49–64. doi:10.1016/j.jfluidstructs.2014.09.011.

5. Conclusions This paper investigates the second-order gap resonance excited by forced body motions and by incident waves through fully nonlinear numerical simulations in the time domain. The associated initial boundary value problem for the velocity potential is solved using boundary element method together with a time stepping scheme. In order to keep the simulations running for a sufficiently long time, some special numerical treatments are also applied to the free surface, such as remeshing, smoothing, jet and thin spray cutting. The methodology is validated through comparing with previous experimental and analytical results. For second-order wave resonance in the gap excited by forced harmonic motions of the surrounding bodies, it depends on the modes of body motions, namely heave, sway and roll. The twin bodies are forced to oscillate with the same parameters, i.e. amplitude, frequency and direction. At such resonance, the wave elevation in the gap cannot grow to infinity due to radiation damping and energy dissipation through the 18

Y. Li

Engineering Analysis with Boundary Elements 106 (2019) 1–19

[18] Sun SY, Sun SL, Ren HL, Wu GX. Splash jet and slamming generated by a rotating flap. Phys Fluids 2015;27:092107. doi:10.1063/1.4931965. [19] Zhou BZ, Wu GX. Resonance of a tension leg platform exited by third-harmonic force in nonlinear regular waves. Philos Trans R Soc A Math Phys Eng Sci 2015;373:20140105. doi:10.1098/rsta.2014.0105. [20] Zhou BZ, Wu GX, Teng B. Fully nonlinear wave interaction with freely floating non-wall-sided structures. Eng Anal Bound Elem 2015;50:117–32. doi:10.1016/j.enganabound.2014.08.003. [21] Longuet-Higgins MS, Cokelet ED. The deformation of steep surface waves on water I. A numerical method of computation. Proc R Soc A Math Phys Eng Sci 1976;350:1– 26. doi:10.1098/rspa.1976.0092. [22] Cointe R, Geyer P, King B, Molin B. Nonlinear and linear motions of a rectangular barge in a perfect fluid. In: Proceedings of the eighteenth symposium on naval hydrodynamics; 1990. [23] Tanizawa K. Long time fully nonlinear simulation of floating body motions with artificial damping zone. J Soc Nav Archit Jpn 1996;180:311–19. doi:10.2534/jjasnaoe1968.1996.180_311. [24] Wu GX, Eatock Taylor R. The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies. Ocean Eng 2003;30:387–400. doi:10.1016/S0029-8018(02)00037-9.

[25] Li Y. Fully nonlinear numerical simulations of wave interactions with multiple structures at resonance. University College London; 2017. [26] Lee CM, Jones H, Bedel JW. Added mass and damping coefficients of heaving twin cylinders in a free surface. NSRDC; 1971. Technical Report No. 3695. [27] Wu GX. Second-order resonance of sloshing in a tank. Ocean Eng 2007;34:2345–9. doi:10.1016/j.oceaneng.2007.05.004. [28] Li Y, Zhang C. Analysis of wave resonance in gap between two heaving barges. Ocean Eng 2016;117:210–20. doi:10.1016/j.oceaneng.2016.03.042. [29] Ning DZ, Song WH, Liu YL, Teng B. A boundary element investigation of liquid sloshing in coupled horizontal and vertical excitation. J Appl Math 2012;2012:340640. doi:10.1155/2012/340640. [30] Zhang C, Li Y, Meng Q. Fully nonlinear analysis of second-order sloshing resonance in a three-dimensional tank. Comput Fluids 2015;116:88–104. doi:10.1016/j.compfluid.2015.04.016. [31] Fenton JD. A fifth-order Stokes theory for steady waves. J Waterw Port Coast Ocean Eng 1985;111:216–34. doi:10.1061/(ASCE)0733-950X(1987)113:4(437). [32] Isaacson M, Cheung K-F. Second order wave diffraction around twodimensional bodies by time-domain method. Appl Ocean Res 1991;13:175–86. doi:10.1016/S0141-1187(05)80073-2.

19