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

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.

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 system 𝑂𝑥𝑦𝑧 is adopted with the origin 𝑂 located at the undisturbed free surface, where the 𝑧-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 𝜌₁, ℎ₁ and 𝜌₂, ℎ₂, respectively. The ratio between the upper and lower layer density can be defined as 𝛾=𝜌₁/𝜌₂ and the total water depth isℎ=ℎ₁+ℎ₂.

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 Φ⁽¹⁾(𝑥,𝑦,𝑧,𝑡) and Φ⁽²⁾(𝑥,𝑦,𝑧,𝑡) in the fluid domain Ω₁ and Ω₂ satisfy Laplace equation: ∇²Φ^(𝑚)=0,(2.1) where 𝑚=1,2 denote the upper and lower fluid domains. The linearized boundary conditions for free surface and interface are described as: ΦΦ^(1)𝑡𝑡+𝑔Φ^𝑧(1)=0,𝑧=0,_𝑧(1)𝛾Φ_(1)𝑡𝑡=Φ^𝑧(2),𝑧=−ℎ₁,+𝑔Φ^𝑧(1)=Φ^(2)𝑡𝑡+𝑔Φ^𝑧(2),𝑧=−ℎ₁.(2.2) A rigid sea bottom satisfies the no-penetration condition in the lower domain is given as: Φ^𝑧(2)=0,𝑧=−ℎ₁−ℎ₂.(2.3) The velocity potential is defined as the following equation: Φ^(𝑚)(𝑥,𝑦,𝑧,𝑡)=Re𝜙^(𝑚)(𝑥,𝑦,𝑧)𝑒^{−𝑖𝜔𝑡}.(2.4) Inserting (2.4) into boundary conditions of (2.2)~(2.3): 𝜕𝜙⁽¹⁾(𝑥,𝑦,𝑧)𝜕𝑧=𝜔𝑔²𝜙⁽¹⁾𝜕𝜙₍₁₎(𝑥,𝑧)=𝑘₀𝜙⁽¹⁾(𝑥,𝑦,𝑧),𝑧=0,(2.5)(𝑥,𝑧)𝜕𝑧=𝜕𝜙⁽²⁾(𝑥,𝑧)𝜕𝑧,𝑧=−ℎ¹𝛾𝜕𝜙₍₁₎,(2.6)(𝑥,𝑦,𝑧)𝜕𝑧−𝑘⁰𝜙⁽¹⁾(𝑥,𝑦,𝑧)=𝜕𝜙⁽²⁾(𝑥,𝑦,𝑧)𝜕𝑧−𝑘⁰𝜕𝜙_𝑧(2)𝜙⁽²⁾(𝑥,𝑦,𝑧),𝑧=−ℎ₁,(2.7)(𝑥,𝑧)𝜕𝑧=0,𝑧=−ℎ¹−ℎ₂,(2.8) where 𝑘₀=𝜔²/𝑔. The solutions of (2.5), (2.6), and (2.8) can be written as the following forms: 𝜙⁽¹⁾(𝑥,𝑦,𝑧)=𝑍₁(𝑧)𝑒^{𝑖𝑘(𝑥cos𝛽+𝑦sin𝛽)},𝑍₁(𝑧)=𝑎₁cosh𝑘(𝑧+ℎ)cosh𝑘ℎ+𝑏¹sinh𝑘𝑧sinh𝑘ℎ₁𝜙⁽²⁾(𝑥,𝑦,𝑧)=𝑍₂(𝑧)𝑒^{𝑖𝑘(𝑥cos𝛽+𝑦sin𝛽)},𝑍₂(𝑧)=𝑎₂cosh𝑘(𝑧+ℎ)cosh𝑘ℎ,(2.9) where 𝑘 denotes the wave number, 𝑎₁, 𝑏₁, and 𝑎₂ are the unknown coefficients. Inserting the expression of 𝜙^(𝑚)(𝑥,𝑦,𝑧) into the free surface boundary conditions (2.5), (2.6), and (2.8), we can get 𝑍₁(𝑧)=𝑎₁𝑍coscho𝑘sh(𝑧𝑘+ℎℎ)+𝑘𝑘⁰−tanh𝑘ℎsinh𝑘𝑧₂(𝑧)=𝑎₁1+𝑘𝑘⁰−tanh𝑘ℎcosh𝑘ℎ¹cosh𝑘ℎsinh𝑘ℎ₂coscho𝑘sh(𝑧𝑘+ℎℎ).(2.10) The wave elevations 𝜂⁽¹⁾ and 𝜂⁽²⁾ have the relations as following: ΦΦ^𝑡(1)=−𝑔𝜂⁽¹⁾,𝑧=0,(2.11)_𝑧(2)=𝜂^𝑡(2),𝑧=−ℎ₁.(2.12) Assuming that 𝜂⁽²⁾=𝐴𝑒^{𝑖(𝑘𝑥cos𝛽+𝑘𝑦sin𝛽−𝜔𝑡)}, the coefficients 𝑎₁ can be obtained based on (2.12): 𝑎₁=−𝑖𝑤𝐴cosh𝑘ℎ𝑘sinh𝑘ℎ₂1+𝑘₀/𝑘−tanh𝑘ℎcosh𝑘ℎ₁cosh𝑘ℎ/sinh𝑘ℎ₂,(2.13) we can then obtain the incident wave potentials and wave elevations as follows: Φ^𝑖(1)=𝜔𝐴𝑘⁰/𝑘sinh𝑘𝑧+cosh𝑘𝑧𝑘₀cosh𝑘ℎ₁−𝑘sinh𝑘ℎ₁Φsin(𝑘𝑥cos𝛽+𝑦sin𝛽−𝜔𝑡)intheupperdomain_𝑖(2)=𝜔𝐴cosh𝑘(𝑧+ℎ)𝑘sinh𝑘ℎ₂𝜂sin(𝑘𝑥cos𝛽+𝑦sin𝛽−𝜔𝑡)inthelowerdomain_𝑖(1)=𝐴𝜂cosh𝑘ℎ₁1−𝑘/𝑘₀tanh𝑘ℎ₁cos(𝑘(𝑥cos𝛽+𝑦sin𝛽)−𝜔𝑡)._𝑖(2)=𝐴cos(𝑘(𝑥cos𝛽+𝑦sin𝛽)−𝜔𝑡),(2.14) where 𝐴 denotes the incident wave amplitude on the interface.

By inserting the expressions of 𝜙^(𝑚)(𝑥,𝑦,𝑧) into (2.7), the dispersion equation in two-layer fluid is obtained as (𝛾+1)tanh𝑘ℎ₁tanh𝑘ℎ₂𝑘²+𝑘₀tanh𝑘ℎ₁tanh𝑘ℎ₂𝑘₀−𝑘²⁰(𝛾+1)tanh𝑘ℎ₁tanh𝑘ℎ₂=0.(2.15) 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 𝑘₁ and 𝑘₂. In order to find the physical significance of 𝑘₁ and 𝑘₂ more explicitly, on the assumption that the depth of each layer fluid goes to infinite, (2.15) can be simplified to the following version: 𝑘𝑘₀²(1−𝛾)−2𝑘𝑘₀+(𝛾+1)=0.(2.16) The two roots of (2.16) are 𝑘𝑘₀=1,𝑘𝑘₀=(1+𝛾)(1−𝛾).(2.17) 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: Φ𝜂⁽¹⁾₍₁₎=Φ^𝑖(1)+Φ^𝑑(1),Φ⁽²⁾=Φ^𝑖(2)+Φ^𝑑(2),=𝜂^𝑖(1)+𝜂^𝑑(1),𝜂⁽²⁾=𝜂^𝑖(2)+𝜂^𝑑(2).(2.18)

For the diffraction potential, it must satisfy the following governing equations and the boundary conditions: 𝜕∇Φ²Φ_{𝑑(𝑚)(1)𝑑𝑧}=0,𝑚=1,2𝜕𝑛=−𝜕Φ𝜕𝑛^(1)𝑖𝑧on𝑆^𝐵Φ𝜂^(2)𝑑𝑧=0on𝑧=−ℎ₁−ℎ₂,(2.19)_(1)𝑑𝑡Φ=Φ^(1)𝑑𝑧(a)on𝑧=0,_(1)𝑑𝑡𝜂=−𝑔𝜂^𝑑(1)(b)(2.20)_(2)𝑑𝑡Φ=Φ^(2)𝑑𝑧(c)_(1)𝑑𝑧=Φ^(2)𝑑𝑧(d)on𝑧=−ℎ₁𝛾Φ^(1)𝑑𝑡+𝑔𝜂^𝑑(2)=Φ^(2)𝑑𝑡+𝑔𝜂^𝑑(2)(e).(2.21) For simplicity, we define a function 𝜑=𝛾Φ^𝑑(1)−Φ^𝑑(2) on 𝑧=−ℎ₁, then the boundary condition (e) of (2.21) can be rewritten as: 𝜑_𝑡=(1−𝛾)𝑔𝜂^𝑑(2)on𝑧=−ℎ₁.(2.22) 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 𝜕𝜂^𝑑(𝑚)𝜕𝑡=𝜕𝜙𝜕𝑧^𝑑(𝑚)−𝜈(𝑟)𝜂^𝑑(𝑚)𝜕Φ_𝑑(𝑚),𝑚=1,2,𝜕𝑡=−𝑔𝜂^𝑑(𝑚)𝜕𝜑−𝜈(𝑟)Φ^𝑑(𝑚),𝑚=1,𝜕𝑡=(1−𝛾)𝑔𝜂^𝑑(𝑚)−𝜈(𝑟)𝜑,𝑚=2,(2.23) where 𝑚=1,2, and 𝜈(𝑟) denotes the damping coefficient is given by 𝑣(𝑟)=⎧⎪⎨⎪⎩𝐶𝜔𝑟𝛽−𝜆𝑟⁰²𝑟₀0𝑟<≤𝑟𝑟₀≤𝑟₀+𝛽𝜆,(2.24) where 𝐶 is the damping coefficient, 𝛽 is the beach breadth coefficient, and 𝜆 is the characteristic wave length. 𝑟₀ and 𝑟₀+𝛽𝜆 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: 𝑅_𝑚(𝑡)=1211−cos2𝜋𝑇𝑡𝑡𝑡≤≥22𝑇𝑇,(2.25) where 𝑇 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:𝛼𝜙^1𝑑−_𝑆𝐵𝜙^1𝑑𝜕𝐺₁𝜕𝑛𝑑𝑠+_𝑆𝐹𝐺₁𝜕𝜙^1𝑑𝜕𝑛𝑑𝑠+_𝑆𝐼𝐺₁𝜕𝜙^1𝑑=−_𝑆𝐵𝜕𝑛𝑑𝑠𝐺₁𝜕𝜙^1𝑑𝜕𝑛𝑑𝑠+_𝑆𝐹𝜙^1𝑑𝜕𝐺₁𝜕𝑛𝑑𝑠+_𝑆𝐼𝜙^1𝑑𝜕𝐺₁𝜕𝑛𝑑𝑠,(2.26) where Green’s function 𝐺₁=−1/4𝜋𝑟 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𝑑𝜕𝑛𝑑𝑠=_𝑆𝐼𝜙^2𝑑𝜕𝐺₂𝜕𝑛𝑑𝑠−𝛼𝜙^2𝑑,(2.27) where 𝐺₂=−1/4𝜋𝑟−1/4𝜋𝑟₂ includes its images about the seabed.𝑟is the distance between the field and the source point, and 𝑟₂ 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: 𝛼=⎧⎪⎨⎪⎩011−solidangle((ssoouurrcceeppooiinnttwouittshiidneththeeddoommaainin))4𝜋(sourcepointontheboundary).(2.28) 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: ⎡⎢⎢⎢⎢⎣𝑎₁₁𝑎₁₂𝑎₁₃𝑎₂₁𝑎₂₂𝑎₂₃𝑎₃₁𝑎₃₂𝑎₃₃⎤⎥⎥⎥⎥⎦⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩𝜕Φ^𝑑(1)𝜕𝑛_𝑆𝐹Φ^𝑑(1)_𝑆𝐵𝜕Φ𝜕𝑛^𝑑(1)_𝑆𝐼⎫⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎭=⎡⎢⎢⎢⎢⎣𝑠₁₁𝑠₁₂𝑠₁₃𝑠₂₁𝑠₂₂𝑠₂₃𝑠₃₁𝑠₃₂𝑠₃₃⎤⎥⎥⎥⎥⎦⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩Φ^𝑑(1)_𝑆𝐹𝜕Φ𝜕𝑛^𝑑(1)_𝑆𝐵Φ^𝑑(1)_𝑆𝐼[𝑏]𝜕Φ_𝑑(2)⎫⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎭(2.29)𝜕𝑛_𝑆𝐼=[𝑡]Φ^𝑑(2)_𝑆𝐼.(2.30) 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: 𝜕Φ^𝑑(1)𝜕𝑛=−𝜕Φ𝜕𝑛^𝑑(2).(2.31) 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: ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑎₁₁𝑎₁₂𝑎₁₃+1𝛾𝑠¹³𝑎₂₁[𝑡]⁻¹[𝑏]𝑎₂₂𝑎₂₃+1𝛾𝑠²³𝑎₃₁[𝑡]⁻¹[𝑏]𝑎₃₂𝑎₃₃+1𝛾𝑠³³[𝑡]⁻¹[𝑏]⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩𝜕Φ^𝑑(1)𝜕𝑛_𝑆𝐹Φ^𝑑(1)_𝑆𝐵𝜕Φ𝜕𝑛^𝑑(1)_𝑆𝐼⎫⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎭=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑠₁₁𝑠₁₂1𝛾𝑠¹³𝑠₂₁𝑠₂₂1𝛾𝑠²³𝑠₃₁𝑠₃₂1𝛾𝑠³³⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩Φ^𝑑(1)_𝑆𝐹𝜕Φ𝜕𝑛^𝑑(1)_𝑆𝐵{𝜑}_𝑆𝐼⎫⎪⎪⎪⎪⎬⎪⎪⎪⎪⎭.(2.32) All the coefficient matrices [𝑎],[𝑠],[𝑡], and [𝑏] 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 ΦΦ^𝑑(1)||_𝑆𝐹=0,𝜂^𝑑(1)=0,_𝑑(2)||_𝑆𝐼𝜑=||0,𝜂^𝑑(2)=0,_𝑆𝐼=0.(2.33) 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 𝛿𝜂^𝑚𝑑𝛿𝑡=𝑓Φ^𝑚𝑑,𝜂^𝑚𝑑,𝑡𝑚=1,2𝛿Φ^1𝑑𝛿𝑡=𝑔Φ^1𝑑,𝜂^1𝑑,𝑡𝛿𝜑^2𝑑𝛿𝑡=ℎ𝜑^2𝑑,𝜂^2𝑑,𝑡.(2.34) 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 are ℎ₁/ℎ=0.7, ℎ₂/ℎ=0.3. The densities of the fluids in the upper and the lower layers are 𝜌₁=998.2kg/m³ and 𝜌₂=1027.2kg/m³. That means 𝛾=0.97. The cylinder has a radius of 𝑎/ℎ=0.5, and a draft of 𝑇/ℎ=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 𝜌₁𝑔𝑎ℎ𝐴, while the dimensionless factor for pitch exciting moment about 𝑦-axis is 𝜌₁𝑔𝑎ℎ²𝐴.

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 Δ𝑡=𝑇/100 (𝑇 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 number 𝑘ℎ=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 and 𝑦-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. The 𝑦-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 𝐿, width 𝐵, and draft 𝐷 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 ℎ₁=48m and ℎ₂=16m. The dimensionless factor for wave forces are 𝜌₁𝑔𝐿𝐷, and for moment is 𝜌₁𝑔𝐿²𝐷. 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. The 𝑦-axis wave moment is the combination of two components: horizontal wave force and vertical force, so the trend is different.

Acknowledgments

This work is financially supported by the National Research Program of China (973-Program, Grant no. 2011CB013703). This work is also partially supported by Projects of Natural Science Foundation of China (Grant no. 51009023, 11072052, 50921001, 51109022).