# A Time-Domain Boundary Element Method for Wave Diffraction in a Two-Layer Fluid

Authors: Ying Gou; Xin-jia Chen; Teng Bin
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

A time-domain numerical model is established based on the higher-order boundary element method (HOBEM) to simulate wave diffraction problem in a two-layer fluid of finite depth. There are two possible incident wave modes (surface-wave mode and internal-wave mode) exist in the incident wave for a prescribed frequency in a two-layer fluid. For surface-wave mode, the hydrodynamic characters of fluid particles are similar to single-layer fluid. For the internal-wave mode, through the definition of a new function respected to velocity potentials of upper and lower fluid on the interface by using matching condition, a single set of linear equations is set up to compute the time histories of wave forces and wave profiles by using a fourth-order Runge-Kutta method. An artificial damping layer is adopted both on the free surface and interface to avoid the wave reflection. Examinations of the accuracy of this time-domain algorithm are carried out for a truncated cylinder and a rectangular barge, and the results demonstrate the effectiveness of this method.


---

## Body

## 1.  Introduction

In the classical view of ocean hydrodynamics, the fluid under consideration is usually assumed to be of constant density. In the real ocean, however, density of sea water changes due to the variations in temperature and salinity in the water depth direction. In the ocean with density stratification, we usually simplify this complex case as a two-layer fluid model. In this model there is a thin layer called the pycnocline, the density of fluid above and below this pycnocline is approximately constant. In this paper we also focus on this two-layer fluid.The internal waves will be generated on the interface between the two fluid layers with different types such as incident solitary wave and periodic wave. For internal solitary wave, the wave length is much longer than the characteristic length of structure. So the internal solitary wave forces acting on structures can be simulated only by the Morison formula [1, 2]. For the other typical harmonic internal waves, the wave length is over a wide range. The diffraction/radiation theory should be adopted when the characteristic length of structure is relative large.Lamb [3] and Landau and Lifshitz [4] discussed some of the types of linear wave motion in a two-layer fluid. It was shown that incident waves in a two-layer fluid can propagate with two different wave numbers for a given frequency, corresponding to the surface-wave mode and internal-wave mode, respectively. Yeung and Nguyen [5] derived the Green functions in a two-layer fluid of finite depth for 3D problems. Ten and Kashiwagi [6] and Kashiwagi et al. [7] studied a two-dimensional radiation∖diffraction problem for a general body floating in a two-layer fluid of finite depth by means of a boundary integral-equation method. Linton and McIver [8] studied the horizontal cylinders interacted with waves in two-layer fluids by multipole expansion method. By using the same method, Sturova [9] solved the problem of wave motions of a heavy fluid driven by the oscillations of a cylinder (radiation problem) and the scattering of an internal wave by a fixed cylinder (diffraction problem). Das and Mandal [10] also used this method to study the problems of radiation of water waves by a submerged sphere in deep water as well as in finite depth water with an ice-cover. You et al. [11] calculated the radiation and diffraction of water waves by a floating circular cylinder in a two-layer fluid of finite depth by using analytical method, the wave exciting forces for a floating circular cylinder due to incident waves of both surface- and internal-wave modes were presented.To authors knowledge, most of the studies of the interaction between internal waves and structures are based on analytic method or frequency-domain approach based on the use of wave Green function. Because of the limitation of analytic method, it couldnot be used in practical engineering. For the wave Green function approach whether in frequency-domain or in time-domain approach, the complex wave Green function in two-layer fluid should be calculated. Therefore, developing a time-domain approach easy to implement to study the wave interactions between internal waves and structures has important significance.In this present work, a time-domain higher-order boundary element method (THOBEM) was developed for internal wave diffraction from a 3D body located in the upper layer fluid. Integral equations in the upper and lower fluid domains are derived by applying the Green’s second identity to simple Green function and velocity potential in each layer, respectively. By the construction of a function on the interface, a single set of linear equations is set up to compute the time histories of wave forces and wave profiles by using a fourth-order Runge-Kutta method. An artificial damping layer is adopted to dissipate the scattering waves on both free surface and interface. The internal wave force and moment on a floating body are calculated and compared with published frequency-domain solution and analytic method, and a relatively good agreement was found.


## 2. Numerical Simulations

A fixed Cartesian coordinate systemOxyz is adopted with the origin O located at the undisturbed free surface, where the z-axis is positive upwards and the body is located in the upper layer fluid, as shown in Figure 1. Denote the densities and depths of the upper and lower fluid layers by ρ1, h1 and ρ2, h2, respectively. The ratio between the upper and lower layer density can be defined as γ=ρ1/ρ2 and the total water depth ish=h1+h2.Figure 1
The coordinate system.


### 2.1. Incident Wave Potential and Dispersion Relation

The expressions of incident wave velocity potential of upper and lower domain in a two-layer fluid had been listed by some researchers [5, 6], and the detail of the present derivative process is listed below.The fluid in each layer is assumed to be inviscid and incompressible, and the flow irrotational. The velocity potentialsΦ(1)(x,y,z,t) and Φ(2)(x,y,z,t) in the fluid domain Ω1 and Ω2 satisfy Laplace equation:
(2.1)∇2Φ(m)=0,
where m=1,2 denote the upper and lower fluid domains. The linearized boundary conditions for free surface and interface are described as:
(2.2)Φ(1)tt+gΦ(1)z=0,z=0,Φ(1)z=Φ(2)z,z=-h1,γ(Φ(1)tt+gΦ(1)z)=Φ(2)tt+gΦ(2)z,z=-h1.
A rigid sea bottom satisfies the no-penetration condition in the lower domain is given as:
(2.3)Φ(2)z=0,z=-h1-h2.
The velocity potential is defined as the following equation:
(2.4)Φ(m)(x,y,z,t)=Re[ϕ(m)(x,y,z)e-iωt].
Inserting (2.4) into boundary conditions of (2.2)~(2.3):
(2.5)∂ϕ(1)(x,y,z)∂z=ω2gϕ(1)(x,z)=k0ϕ(1)(x,y,z),z=0,(2.6)∂ϕ(1)(x,z)∂z=∂ϕ(2)(x,z)∂z,z=-h1,(2.7)γ(∂ϕ(1)(x,y,z)∂z-k0ϕ(1)(x,y,z))=(∂ϕ(2)(x,y,z)∂z-k0ϕ(2)(x,y,z)),z=-h1,(2.8)∂ϕz(2)(x,z)∂z=0,z=-h1-h2,
where k0=ω2/g. The solutions of (2.5), (2.6), and (2.8) can be written as the following forms:
(2.9)ϕ(1)(x,y,z)=Z1(z)eik(xcos⁡β+ysin⁡β),Z1(z)=a1cosh⁡k(z+h)cosh⁡kh+b1sinh⁡kzsinh⁡kh1ϕ(2)(x,y,z)=Z2(z)eik(xcos⁡β+ysin⁡β),Z2(z)=a2cosh⁡k(z+h)cosh⁡kh,
where k denotes the wave number, a1, b1, and a2 are the unknown coefficients. Inserting the expression of ϕ(m)(x,y,z) into the free surface boundary conditions (2.5), (2.6), and (2.8), we can get
(2.10)Z1(z)=a1[cosh⁡k(z+h)cosh⁡kh+(k0k-tanhkh)sinh⁡kz]Z2(z)=a1[1+(k0k-tanhkh)cosh⁡kh1cosh⁡khsinh⁡kh2]cosh⁡k(z+h)cosh⁡kh.
The wave elevations η(1) and η(2) have the relations as following:
(2.11)Φt(1)=-gη(1),z=0,(2.12)Φz(2)=ηt(2),z=-h1.
Assuming that η(2)=Aei(kxcosβ+kysinβ-ωt), the coefficients a1 can be obtained based on (2.12):
(2.13)a1=-iwAcosh⁡khksinh⁡kh2[1+(k0/k-tanhkh)(cosh⁡kh1cosh⁡kh/sinh⁡kh2)],
we can then obtain the incident wave potentials and wave elevations as follows:
(2.14)Φi(1)=ωA((k0/k)sinh⁡kz+cosh⁡kz)k0cosh⁡kh1-ksinh⁡kh1sin⁡(kxcos⁡β+ysin⁡β-ωt)intheupperdomainΦi(2)=ωAcosh⁡k(z+h)ksinh⁡kh2sin⁡(kxcos⁡β+ysin⁡β-ωt)inthelowerdomainηi(1)=Acosh⁡kh1[1-(k/k0)tanhkh1]cos⁡(k(xcos⁡β+ysin⁡β)-ωt).ηi(2)=Acos⁡(k(xcos⁡β+ysin⁡β)-ωt),
where A denotes the incident wave amplitude on the interface.By inserting the expressions ofϕ(m)(x,y,z) into (2.7), the dispersion equation in two-layer fluid is obtained as
(2.15)[(γ+1)tanhkh1tanhkh2]k2+[k0(tanhkh1tanhkh2)]k0-k02(γ+1)tanhkh1tanhkh2=0.
It can be seen that three are two possible roots of (2.15), that is, traveling waves with a prescribed frequency ω in a two-layer fluid can propagate with two possible different wavenumbers k1 and k2. In order to find the physical significance of k1 and k2 more explicitly, on the assumption that the depth of each layer fluid goes to infinite, (2.15) can be simplified to the following version:
(2.16)(kk0)2(1-γ)-2kk0+(γ+1)=0.
The two roots of (2.16) are
(2.17)kk0=1,kk0=(1+γ)(1-γ).
It can be seen obviously that the former root is corresponding to the surface wave mode, and the other is corresponding to the internal wave mode. In this paper we only concern the internal wave mode.


### 2.2. Diffraction Problem

Under the assumptions of linear water wave theory, the total velocity potential in upper and lower fluid domain can be divided into a known incident potential and an unknown diffractive potential, while the wave elevation on the free surface and interface can be divided into a known incident wave elevation and an unknown diffractive wave elevation as the following:(2.18)Φ(1)=Φi(1)+Φd(1),Φ(2)=Φi(2)+Φd(2),η(1)=ηi(1)+ηd(1),η(2)=ηi(2)+ηd(2).For the diffraction potential, it must satisfy the following governing equations and the boundary conditions:(2.19)∇2Φd(m)=0,m=1,2∂Φdz(1)∂n=-∂Φiz(1)∂non⁡SBΦdz(2)=0on⁡z=-h1-h2,(2.20)ηdt(1)=Φdz(1)(a)on⁡z=0,Φdt(1)=-gηd(1)(b)(2.21)ηdt(2)=Φdz(2)(c)Φdz(1)=Φdz(2)(d)on⁡z=-h1γ(Φdt(1)+gηd(2))=Φdt(2)+gηd(2)(e).
For simplicity, we define a function φ=γΦd(1)-Φd(2) on z=-h1, then the boundary condition (e) of (2.21) can be rewritten as:
(2.22)φt=(1-γ)gηd(2)on⁡z=-h1.
The velocity potential on interface is discontinuous, but the construction of this function φ makes it easier to solve this boundary value problem. The details are introduced in the next subsection.In order to avoid the reflection of scatter waves on the outer computational boundary, we introduce an artificial damping layer to absorb the scattered wave energy. On the outer part of the free surface and interface, a damping term is added to the free surface boundary conditions(2.23)∂ηd(m)∂t=∂ϕd(m)∂z-ν(r)ηd(m),m=1,2,∂Φd(m)∂t=-gηd(m)-ν(r)Φd(m),m=1,∂φ∂t=(1-γ)gηd(m)-ν(r)φ,m=2,
where m=1,2, and ν(r) denotes the damping coefficient is given by
(2.24)v(r)={Cω(r-r0βλ)2(r0≤r≤r0+βλ)0(r<r0),
where C is the damping coefficient, β is the beach breadth coefficient, and λ is the characteristic wave length. r0 and r0+βλ represent the inner and outer radius of the damping layer, respectively. In this paper values of α and β are equal to 1.0.A ramping function is utilized when the simulation started. Its aim is to make the numerical simulation more stable and reach steady state faster. In this present simulation, a ramp function is applied to the first two waves by the following forms:(2.25)Rm(t)={12(1-cos⁡(πt2T))t≤2T1t≥2T,
where T is the wave period.


### 2.3. Integral Equation and Its Solution by Time-Domain Method

Based on Green’s theorem, the boundary integral equations in the upper and lower layer fluid can be obtained to solve the above boundary value problems. The integral equation in the upper fluid includes the integration over the body surface, the free surface, and the interface:(2.26)αϕd1-∬SBϕd1∂G1∂nds+∬SFG1∂ϕd1∂nds+∬SIG1∂ϕd1∂nds=-∬SBG1∂ϕd1∂nds+∬SFϕd1∂G1∂nds+∬SIϕd1∂G1∂nds,
where Green’s function G1=-1/4πr is the Rankine source. In our case the body located in the upper layer is considered, thus the integral equation in the lower fluid only includes the integration over the interface. The integral equation in the lower layer domain has the forms as
(2.27)∬SIG2∂ϕd2∂nds=∬SIϕd2∂G2∂nds-αϕd2,
where G2=-1/4πr-1/4πr2 includes its images about the seabed.ris the distance between the field and the source point, and r2 is the distance between the field point and the image of the source point about the sea bed. α is the so-called solid angle coefficient and satisfies the following forms:
(2.28)α={1(sourcepointwithinthedomain)0(sourcepointoutsidethedomain)1-solidangle4π(sourcepointon⁡theboundary).
In this paper the direct method is used to calculate the solid angle [12].Compared with constant panel method the higher-order boundary element method [13] is more accurate and efficient. In this numerical method, quadratic isoparametric elements are used to discretize the surface over which the integral is performed. Curved quadrilateral and triangular elements are placed by eight and six nodes, respectively.Through the discretization of (2.26) and (2.27), two sets of linear equations can be obtained as:
(2.29)[a11a12a13a21a22a23a31a32a33]{{∂Φd(1)∂n}SF{Φd(1)}SB{∂Φd(1)∂n}SI}=[s11s12s13s21s22s23s31s32s33]{{Φd(1)}SF{∂Φd(1)∂n}SB{Φd(1)}SI}(2.30)[b]{∂Φd(2)∂n}SI=[t]{Φd(2)}SI.
It is worth noting that the normal velocities of the upper and lower domains on the interface are in the opposite direction. It means that the normal derivatives of the potentials at the interface have the following relation:
(2.31)∂Φd(1)∂n=-∂Φd(2)∂n.
Applying the above relation to (2.30) and taking advantage of the function φ constructed in (2.22), we can combine (2.29) and (2.30) to get a single set of linear equations as follows:
(2.32)[a11a12a13+1γs13[t]-1[b]a21a22a23+1γs23[t]-1[b]a31a32a33+1γs33[t]-1[b]]{{∂Φd(1)∂n}SF{Φd(1)}SB{∂Φd(1)∂n}SI}=[s11s121γs13s21s221γs23s31s321γs33]{{Φd(1)}SF{∂Φd(1)∂n}SB{φ}SI}.
All the coefficient matrices [a],[s],[t], and [b] are constant at all time steps due to the time-independent integral boundary.A time-stepping integration method is employed to obtain the values of wave elevations, the velocity potentials over the body surface, and the derivative of interface. The initial conditions as(2.33)Φd(1)|SF=0,ηd(1)=0,Φd(2)|SI=0,ηd(2)=0,φ|SI=0.
After solving the boundary value problem and obtaining the fluid velocity on the free surface and interface at each time step, the general forms of the dynamic and kinematic boundary conditions (2.23) are used as ordinary differential equations to be marched in time and can be rewritten as
(2.34)δηdmδt=f(Φdm,ηdm,t)m=1,2δΦd1δt=g(Φd1,ηd1,t)δφd2δt=h(φd2,ηd2,t).
A 4th-order Runge-Kutta method is adopted for integrating (2.34) with time.When the unknown velocity potentials over body surface are obtained, the pressure can be derived from the Bernoulli equation, the forces and moments acting on the body can be calculated by integrating the pressure over the mean body surface.


## 2.1. Incident Wave Potential and Dispersion Relation

The expressions of incident wave velocity potential of upper and lower domain in a two-layer fluid had been listed by some researchers [5, 6], and the detail of the present derivative process is listed below.The fluid in each layer is assumed to be inviscid and incompressible, and the flow irrotational. The velocity potentialsΦ(1)(x,y,z,t) and Φ(2)(x,y,z,t) in the fluid domain Ω1 and Ω2 satisfy Laplace equation:
(2.1)∇2Φ(m)=0,
where m=1,2 denote the upper and lower fluid domains. The linearized boundary conditions for free surface and interface are described as:
(2.2)Φ(1)tt+gΦ(1)z=0,z=0,Φ(1)z=Φ(2)z,z=-h1,γ(Φ(1)tt+gΦ(1)z)=Φ(2)tt+gΦ(2)z,z=-h1.
A rigid sea bottom satisfies the no-penetration condition in the lower domain is given as:
(2.3)Φ(2)z=0,z=-h1-h2.
The velocity potential is defined as the following equation:
(2.4)Φ(m)(x,y,z,t)=Re[ϕ(m)(x,y,z)e-iωt].
Inserting (2.4) into boundary conditions of (2.2)~(2.3):
(2.5)∂ϕ(1)(x,y,z)∂z=ω2gϕ(1)(x,z)=k0ϕ(1)(x,y,z),z=0,(2.6)∂ϕ(1)(x,z)∂z=∂ϕ(2)(x,z)∂z,z=-h1,(2.7)γ(∂ϕ(1)(x,y,z)∂z-k0ϕ(1)(x,y,z))=(∂ϕ(2)(x,y,z)∂z-k0ϕ(2)(x,y,z)),z=-h1,(2.8)∂ϕz(2)(x,z)∂z=0,z=-h1-h2,
where k0=ω2/g. The solutions of (2.5), (2.6), and (2.8) can be written as the following forms:
(2.9)ϕ(1)(x,y,z)=Z1(z)eik(xcos⁡β+ysin⁡β),Z1(z)=a1cosh⁡k(z+h)cosh⁡kh+b1sinh⁡kzsinh⁡kh1ϕ(2)(x,y,z)=Z2(z)eik(xcos⁡β+ysin⁡β),Z2(z)=a2cosh⁡k(z+h)cosh⁡kh,
where k denotes the wave number, a1, b1, and a2 are the unknown coefficients. Inserting the expression of ϕ(m)(x,y,z) into the free surface boundary conditions (2.5), (2.6), and (2.8), we can get
(2.10)Z1(z)=a1[cosh⁡k(z+h)cosh⁡kh+(k0k-tanhkh)sinh⁡kz]Z2(z)=a1[1+(k0k-tanhkh)cosh⁡kh1cosh⁡khsinh⁡kh2]cosh⁡k(z+h)cosh⁡kh.
The wave elevations η(1) and η(2) have the relations as following:
(2.11)Φt(1)=-gη(1),z=0,(2.12)Φz(2)=ηt(2),z=-h1.
Assuming that η(2)=Aei(kxcosβ+kysinβ-ωt), the coefficients a1 can be obtained based on (2.12):
(2.13)a1=-iwAcosh⁡khksinh⁡kh2[1+(k0/k-tanhkh)(cosh⁡kh1cosh⁡kh/sinh⁡kh2)],
we can then obtain the incident wave potentials and wave elevations as follows:
(2.14)Φi(1)=ωA((k0/k)sinh⁡kz+cosh⁡kz)k0cosh⁡kh1-ksinh⁡kh1sin⁡(kxcos⁡β+ysin⁡β-ωt)intheupperdomainΦi(2)=ωAcosh⁡k(z+h)ksinh⁡kh2sin⁡(kxcos⁡β+ysin⁡β-ωt)inthelowerdomainηi(1)=Acosh⁡kh1[1-(k/k0)tanhkh1]cos⁡(k(xcos⁡β+ysin⁡β)-ωt).ηi(2)=Acos⁡(k(xcos⁡β+ysin⁡β)-ωt),
where A denotes the incident wave amplitude on the interface.By inserting the expressions ofϕ(m)(x,y,z) into (2.7), the dispersion equation in two-layer fluid is obtained as
(2.15)[(γ+1)tanhkh1tanhkh2]k2+[k0(tanhkh1tanhkh2)]k0-k02(γ+1)tanhkh1tanhkh2=0.
It can be seen that three are two possible roots of (2.15), that is, traveling waves with a prescribed frequency ω in a two-layer fluid can propagate with two possible different wavenumbers k1 and k2. In order to find the physical significance of k1 and k2 more explicitly, on the assumption that the depth of each layer fluid goes to infinite, (2.15) can be simplified to the following version:
(2.16)(kk0)2(1-γ)-2kk0+(γ+1)=0.
The two roots of (2.16) are
(2.17)kk0=1,kk0=(1+γ)(1-γ).
It can be seen obviously that the former root is corresponding to the surface wave mode, and the other is corresponding to the internal wave mode. In this paper we only concern the internal wave mode.


## 2.2. Diffraction Problem

Under the assumptions of linear water wave theory, the total velocity potential in upper and lower fluid domain can be divided into a known incident potential and an unknown diffractive potential, while the wave elevation on the free surface and interface can be divided into a known incident wave elevation and an unknown diffractive wave elevation as the following:(2.18)Φ(1)=Φi(1)+Φd(1),Φ(2)=Φi(2)+Φd(2),η(1)=ηi(1)+ηd(1),η(2)=ηi(2)+ηd(2).For the diffraction potential, it must satisfy the following governing equations and the boundary conditions:(2.19)∇2Φd(m)=0,m=1,2∂Φdz(1)∂n=-∂Φiz(1)∂non⁡SBΦdz(2)=0on⁡z=-h1-h2,(2.20)ηdt(1)=Φdz(1)(a)on⁡z=0,Φdt(1)=-gηd(1)(b)(2.21)ηdt(2)=Φdz(2)(c)Φdz(1)=Φdz(2)(d)on⁡z=-h1γ(Φdt(1)+gηd(2))=Φdt(2)+gηd(2)(e).
For simplicity, we define a function φ=γΦd(1)-Φd(2) on z=-h1, then the boundary condition (e) of (2.21) can be rewritten as:
(2.22)φt=(1-γ)gηd(2)on⁡z=-h1.
The velocity potential on interface is discontinuous, but the construction of this function φ makes it easier to solve this boundary value problem. The details are introduced in the next subsection.In order to avoid the reflection of scatter waves on the outer computational boundary, we introduce an artificial damping layer to absorb the scattered wave energy. On the outer part of the free surface and interface, a damping term is added to the free surface boundary conditions(2.23)∂ηd(m)∂t=∂ϕd(m)∂z-ν(r)ηd(m),m=1,2,∂Φd(m)∂t=-gηd(m)-ν(r)Φd(m),m=1,∂φ∂t=(1-γ)gηd(m)-ν(r)φ,m=2,
where m=1,2, and ν(r) denotes the damping coefficient is given by
(2.24)v(r)={Cω(r-r0βλ)2(r0≤r≤r0+βλ)0(r<r0),
where C is the damping coefficient, β is the beach breadth coefficient, and λ is the characteristic wave length. r0 and r0+βλ represent the inner and outer radius of the damping layer, respectively. In this paper values of α and β are equal to 1.0.A ramping function is utilized when the simulation started. Its aim is to make the numerical simulation more stable and reach steady state faster. In this present simulation, a ramp function is applied to the first two waves by the following forms:(2.25)Rm(t)={12(1-cos⁡(πt2T))t≤2T1t≥2T,
where T is the wave period.


## 2.3. Integral Equation and Its Solution by Time-Domain Method

Based on Green’s theorem, the boundary integral equations in the upper and lower layer fluid can be obtained to solve the above boundary value problems. The integral equation in the upper fluid includes the integration over the body surface, the free surface, and the interface:(2.26)αϕd1-∬SBϕd1∂G1∂nds+∬SFG1∂ϕd1∂nds+∬SIG1∂ϕd1∂nds=-∬SBG1∂ϕd1∂nds+∬SFϕd1∂G1∂nds+∬SIϕd1∂G1∂nds,
where Green’s function G1=-1/4πr is the Rankine source. In our case the body located in the upper layer is considered, thus the integral equation in the lower fluid only includes the integration over the interface. The integral equation in the lower layer domain has the forms as
(2.27)∬SIG2∂ϕd2∂nds=∬SIϕd2∂G2∂nds-αϕd2,
where G2=-1/4πr-1/4πr2 includes its images about the seabed.ris the distance between the field and the source point, and r2 is the distance between the field point and the image of the source point about the sea bed. α is the so-called solid angle coefficient and satisfies the following forms:
(2.28)α={1(sourcepointwithinthedomain)0(sourcepointoutsidethedomain)1-solidangle4π(sourcepointon⁡theboundary).
In this paper the direct method is used to calculate the solid angle [12].Compared with constant panel method the higher-order boundary element method [13] is more accurate and efficient. In this numerical method, quadratic isoparametric elements are used to discretize the surface over which the integral is performed. Curved quadrilateral and triangular elements are placed by eight and six nodes, respectively.Through the discretization of (2.26) and (2.27), two sets of linear equations can be obtained as:
(2.29)[a11a12a13a21a22a23a31a32a33]{{∂Φd(1)∂n}SF{Φd(1)}SB{∂Φd(1)∂n}SI}=[s11s12s13s21s22s23s31s32s33]{{Φd(1)}SF{∂Φd(1)∂n}SB{Φd(1)}SI}(2.30)[b]{∂Φd(2)∂n}SI=[t]{Φd(2)}SI.
It is worth noting that the normal velocities of the upper and lower domains on the interface are in the opposite direction. It means that the normal derivatives of the potentials at the interface have the following relation:
(2.31)∂Φd(1)∂n=-∂Φd(2)∂n.
Applying the above relation to (2.30) and taking advantage of the function φ constructed in (2.22), we can combine (2.29) and (2.30) to get a single set of linear equations as follows:
(2.32)[a11a12a13+1γs13[t]-1[b]a21a22a23+1γs23[t]-1[b]a31a32a33+1γs33[t]-1[b]]{{∂Φd(1)∂n}SF{Φd(1)}SB{∂Φd(1)∂n}SI}=[s11s121γs13s21s221γs23s31s321γs33]{{Φd(1)}SF{∂Φd(1)∂n}SB{φ}SI}.
All the coefficient matrices [a],[s],[t], and [b] are constant at all time steps due to the time-independent integral boundary.A time-stepping integration method is employed to obtain the values of wave elevations, the velocity potentials over the body surface, and the derivative of interface. The initial conditions as(2.33)Φd(1)|SF=0,ηd(1)=0,Φd(2)|SI=0,ηd(2)=0,φ|SI=0.
After solving the boundary value problem and obtaining the fluid velocity on the free surface and interface at each time step, the general forms of the dynamic and kinematic boundary conditions (2.23) are used as ordinary differential equations to be marched in time and can be rewritten as
(2.34)δηdmδt=f(Φdm,ηdm,t)m=1,2δΦd1δt=g(Φd1,ηd1,t)δφd2δt=h(φd2,ηd2,t).
A 4th-order Runge-Kutta method is adopted for integrating (2.34) with time.When the unknown velocity potentials over body surface are obtained, the pressure can be derived from the Bernoulli equation, the forces and moments acting on the body can be calculated by integrating the pressure over the mean body surface.


## 3. Numerical Validation of the Time-Domain Method

### 3.1. A Truncated Cylinder in a Two-Layer Fluid of Finite Depth

The water depths of the upper and the lower layers areh1/h=0.7, h2/h=0.3. The densities of the fluids in the upper and the lower layers are ρ1=998.2kg/m3 and ρ2=1027.2kg/m3. That means γ=0.97. The cylinder has a radius of a/h=0.5, and a draft of T/h=0.5. The truncated cylinder is located in the upper fluid, as shown in Figure 2. In the following calculations, the surge and heave exciting forces have been nondimensionalized by ρ1gahA, while the dimensionless factor for pitch exciting moment about y-axis is ρ1gah2A.Figure 2
Definition of sketch for a truncated cylinder.The structural meshes are used to discretize the body surface, the free surface, and the interface. If the mesh size or the time step is too large, the convergence will be difficult to achieve. Conversely, it will take too much computation time if a large number of meshes are used. The numerical convergence tests are desired to ensure the mesh size, and the time step is suitable. Based on the numerical convergence tests, the radius of the free surface should be two to three times of the wave length. The mesh size on the free surface and interface should be about ten percent of the wave length, and the sizes of elements should be uniformly distributed. The time step for time marching is taken asΔt=T/100 (T is the wave period).The time histories of the wave forces and the wave moments on the cylinder due to internal wave mode at the dimensionless wave numberkh=4.0 are plotted in Figures 3, 4, and 5. It can be seen that very steady time-history results can be obtained by this time-domain method.Figure 3
Time histories of horizontal wave force.Figure 4
Time histories of vertical wave force.Figure 5
Time histories of pitch moment.To examine the accuracy of the present model, the calculation of horizontal wave force, vertical wave force andy-axis wave moment is carried out by the present time-domain method and the analytic solutions of results of  You et al. [11] are also shown to compare. The results are shown in Figures 6, 7, and 8. The dimensionless wave frequency varied from 0.1 to 0.4.Figure 6
Horizontal wave force on cylinder.Figure 7
Vertical wave force on cylinder.Figure 8
Pitch moment on cylinder.From these figures, it is can be seen that the present results are in close agreement with the results obtained by the analytical solutions. The magnitudes of the horizontal and vertical wave forces decrease as dimensionless wave frequency decreases within the calculated range. They-axis wave moment is the combination of two components: due to horizontal wave force and vertical force, so the trend is different.


### 3.2. A Rectangular Box in a Two-Layer Fluid of Finite Depth

To further validate the present method, the wave diffraction of a rectangular box is considered. In order to make a comparison with the numerical results of Nguyen and Yeung [14] by the frequency-domain method, the same box is studied which the length L, width B, and draft D are 90 m, 90 m, and 40 m, respectively. The box is also located in the upper fluid with no forward speed. The water depths of the upper and the lower layers are h1=48m and h2=16m. The dimensionless factor for wave forces are ρ1gLD, and for moment is ρ1gL2D. The results are shown in Figures 9, 10, and 11.Figure 9
Horizontal wave force on box.Figure 10
Vertical wave force on box.Figure 11
Pitch moment on box.Three types of density ratioγ=0.9, 0.7, and 0.1 are adopted. From the comparisons, we can find that the present results by time-domain method also have a good agreement with the frequency-domain method. From these figures, it is clear to see that the density ratio has obvious effects on the hydrodynamic forces of internal-wave mode. As the density ratio decrease, for a given frequencies, the magnitudes of the horizontal and vertical forces increase. As we have mentioned before, the trends of the magnitude of the wave moment is different.


## 3.1. A Truncated Cylinder in a Two-Layer Fluid of Finite Depth

The water depths of the upper and the lower layers areh1/h=0.7, h2/h=0.3. The densities of the fluids in the upper and the lower layers are ρ1=998.2kg/m3 and ρ2=1027.2kg/m3. That means γ=0.97. The cylinder has a radius of a/h=0.5, and a draft of T/h=0.5. The truncated cylinder is located in the upper fluid, as shown in Figure 2. In the following calculations, the surge and heave exciting forces have been nondimensionalized by ρ1gahA, while the dimensionless factor for pitch exciting moment about y-axis is ρ1gah2A.Figure 2
Definition of sketch for a truncated cylinder.The structural meshes are used to discretize the body surface, the free surface, and the interface. If the mesh size or the time step is too large, the convergence will be difficult to achieve. Conversely, it will take too much computation time if a large number of meshes are used. The numerical convergence tests are desired to ensure the mesh size, and the time step is suitable. Based on the numerical convergence tests, the radius of the free surface should be two to three times of the wave length. The mesh size on the free surface and interface should be about ten percent of the wave length, and the sizes of elements should be uniformly distributed. The time step for time marching is taken asΔt=T/100 (T is the wave period).The time histories of the wave forces and the wave moments on the cylinder due to internal wave mode at the dimensionless wave numberkh=4.0 are plotted in Figures 3, 4, and 5. It can be seen that very steady time-history results can be obtained by this time-domain method.Figure 3
Time histories of horizontal wave force.Figure 4
Time histories of vertical wave force.Figure 5
Time histories of pitch moment.To examine the accuracy of the present model, the calculation of horizontal wave force, vertical wave force andy-axis wave moment is carried out by the present time-domain method and the analytic solutions of results of  You et al. [11] are also shown to compare. The results are shown in Figures 6, 7, and 8. The dimensionless wave frequency varied from 0.1 to 0.4.Figure 6
Horizontal wave force on cylinder.Figure 7
Vertical wave force on cylinder.Figure 8
Pitch moment on cylinder.From these figures, it is can be seen that the present results are in close agreement with the results obtained by the analytical solutions. The magnitudes of the horizontal and vertical wave forces decrease as dimensionless wave frequency decreases within the calculated range. They-axis wave moment is the combination of two components: due to horizontal wave force and vertical force, so the trend is different.


## 3.2. A Rectangular Box in a Two-Layer Fluid of Finite Depth

To further validate the present method, the wave diffraction of a rectangular box is considered. In order to make a comparison with the numerical results of Nguyen and Yeung [14] by the frequency-domain method, the same box is studied which the length L, width B, and draft D are 90 m, 90 m, and 40 m, respectively. The box is also located in the upper fluid with no forward speed. The water depths of the upper and the lower layers are h1=48m and h2=16m. The dimensionless factor for wave forces are ρ1gLD, and for moment is ρ1gL2D. The results are shown in Figures 9, 10, and 11.Figure 9
Horizontal wave force on box.Figure 10
Vertical wave force on box.Figure 11
Pitch moment on box.Three types of density ratioγ=0.9, 0.7, and 0.1 are adopted. From the comparisons, we can find that the present results by time-domain method also have a good agreement with the frequency-domain method. From these figures, it is clear to see that the density ratio has obvious effects on the hydrodynamic forces of internal-wave mode. As the density ratio decrease, for a given frequencies, the magnitudes of the horizontal and vertical forces increase. As we have mentioned before, the trends of the magnitude of the wave moment is different.


## 4. Conclusions

A time-domain higher-order boundary element method is developed and it can be used to simulate the wave diffraction problem in a two-layer fluid. Through the construction of a function on the interface, the two integral equations in the upper and lower layers are combined, and a single set of linear equations are set up. The wave diffraction problem are solved with a fourth-order Runge-Kutta method for time marching.In order to verify the validity of the present approach, a truncated cylinder and a rectangular box are taken as examples. Compared with the analytic solution and the numerical results by frequency-domain approach, the computation in the present time-domain approach agrees well with them. So the present time-domain approach can be used to study the wave diffraction problem in two-layer fluid and it can be extend to simulate the interaction between internal wave and floating structure.In the present study, we find that the density ratio has obvious effects on the hydrodynamic forces of internal-wave mode. As the density ratio decrease, for a given frequencies, the magnitudes of the horizontal and vertical forces increase. They-axis wave moment is the combination of two components: horizontal wave force and vertical force, so the trend is different.



---

*Source: 686824-2012-07-01.xml*