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

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.


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 2 Journal of Applied Mathematics 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 integralequation 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 surfaceand 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 fourthorder 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.

Numerical Simulations
A fixed Cartesian coordinate system Oxyz 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 , h 1 and ρ 2 , h 2 , respectively. The ratio between the upper and lower layer density can be defined as γ ρ 1 /ρ 2 and the total water depth is h h 1 h 2 .

Incident Wave Potential and Dispersion Relation
The expressions of incident wave velocity potential of upper and lower domain in a twolayer 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: 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
A rigid sea bottom satisfies the no-penetration condition in the lower domain is given as: The velocity potential is defined as the following equation: where k 0 ω 2 /g. The solutions of 2.5 , 2.6 , and 2.8 can be written as the following forms: where k denotes the wave number, a 1 , b 1 , and a 2 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
The wave elevations η 1 and η 2 have the relations as following: Assuming that η 2 Ae i kx cos β ky sin β−ωt , the coefficients a 1 can be obtained based on 2.12 : Journal of Applied Mathematics 5 we can then obtain the incident wave potentials and wave elevations as follows: 2.14 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 twolayer fluid is obtained as 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 k 1 and k 2 . In order to find the physical significance of k 1 and k 2 more explicitly, on the assumption that the depth of each layer fluid goes to infinite, 2.15 can be simplified to the following version: The two roots of 2.16 are 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.

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 6 Journal of Applied Mathematics into a known incident wave elevation and an unknown diffractive wave elevation as the following:

2.18
For the diffraction potential, it must satisfy the following governing equations and the boundary conditions:

2.21
For simplicity, we define a function ϕ γΦ , then the boundary condition e of 2.21 can be rewritten as: 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 where m 1, 2, and ν r denotes the damping coefficient is given by where C is the damping coefficient, β is the beach breadth coefficient, and λ is the characteristic wave length. r 0 and r 0 βλ 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: where T is the wave period.

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: where G 2 −1/4πr − 1/4πr 2 includes its images about the seabed.ris the distance between the field and the source point, and r 2 is the distance between the field point and the image of

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: a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33

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: Journal of Applied Mathematics 9 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: All the coefficient matrices a , s , t , and b are constant at all time steps due to the timeindependent 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
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
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.

A Truncated Cylinder in a Two-Layer Fluid of Finite Depth
The water depths of the upper and the lower layers are h 1 /h 0. calculations, the surge and heave exciting forces have been nondimensionalized by ρ 1 gahA, while the dimensionless factor for pitch exciting moment about y-axis is ρ 1 gah 2 A. 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 number kh 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.
To results are shown in Figures 6, 7, and 8. The dimensionless wave frequency varied from 0.1 to 0.4. 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 y-axis wave moment is the combination of two components: due to horizontal wave force and vertical force, so the trend is different.

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    h 2 16 m. The dimensionless factor for wave forces are ρ 1 gLD, and for moment is ρ 1 gL 2 D. The results are shown in Figures 9, 10, and 11.
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.

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 timedomain 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 y-axis wave moment is the combination of two components: horizontal wave force and vertical force, so the trend is different.