Research Article | Open Access

Ying Gou, Xin-jia Chen, Teng Bin, "A Time-Domain Boundary Element Method for Wave Diffraction in a Two-Layer Fluid", *Journal of Applied Mathematics*, vol. 2012, Article ID 686824, 15 pages, 2012. https://doi.org/10.1155/2012/686824

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

**Academic Editor:**Ioannis K. Chatjigeorgiou

#### 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 radiationdiffraction 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.

##### 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: where denote the upper and lower fluid domains. The linearized boundary conditions for free surface and interface are described as: 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: Inserting (2.4) into boundary conditions of (2.2)~(2.3): where . The solutions of (2.5), (2.6), and (2.8) can be written as the following forms: 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 The wave elevations and have the relations as following: Assuming that , the coefficients can be obtained based on (2.12): we can then obtain the incident wave potentials and wave elevations as follows: 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 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: 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.

##### 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:

For the diffraction potential, it must satisfy the following governing equations and the boundary conditions: For simplicity, we define a function on , 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 , and denotes the damping coefficient is given by 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: 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: where Greenโs function 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 where 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: 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: 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: 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 , 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 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 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 , . The densities of the fluids in the upper and the lower layers are kg/m^{3} and kg/m^{3}. That means . The cylinder has a radius of , and a draft of . 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 .

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 ( 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 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 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.

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 and . The dimensionless factor for wave forces are , and for moment is *.* The results are shown in Figures 9, 10, and 11.

Three types of density ratio , 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).

#### References

- Y. L. Cheng, J. C. Li, and Y. F. Liu, โThe induced flow field by internal solitary wave and it action on cylindrical piles in the stratified ocean,โ in
*Proceedings of the 4th International Conference on Fluid Mechanics*, pp. 296โ299, 2004. View at: Google Scholar - Z. J. Song, B. Teng, Y. Gou et al., โComparisons of internal solitary wave and surface wave actions on marine structures and their responses,โ
*Applied Ocean Research*, vol. 33, no. 2, pp. 120โ129, 2011. View at: Google Scholar - H. Lamb,
*Hydrodynamics*, Cambridge University Press, Cambridge, UK, 1993. View at: Zentralblatt MATH - L. D. Landau and F. M. Lifshitz,
*Fluid Mechanics*, Pergamon Press, Oxford, UK, 1959, (translation into Japanese by H. Takeuchi, published in 1970 by Tokyo-Tosho). - R. W. Yeung and T. Nguyen, โRadiation and diffraction of waves in a two-layer fluid,โ in
*Proceedings of the 22nd Symposium on Naval Hydrodynamics*, pp. 875โ891, Washington, DC, USA, 1999. View at: Google Scholar - I. Ten and M. Kashiwagi, โHydrodynamics of a body floating in a two-layer fluid of finite depth. Part 1. Radiation problem,โ
*Journal of Marine Science and Technology*, vol. 9, pp. 127โ141, 2004. View at: Google Scholar - M. Kashiwagi, I. Ten, and M. Yasunaga, โHydrodynamics of a body floating in a two-layer fluid of finite depth. Part 2. Diffraction problem,โ
*Journal of Marine Science and Technology*, vol. 11, pp. 150โ164, 2006. View at: Google Scholar - C. M. Linton and M. McIver, โThe interaction of waves with horizontal cylinders in two-layer fluids,โ
*Journal of Fluid Mechanics*, vol. 304, pp. 213โ229, 1995. View at: Publisher Site | Google Scholar | Zentralblatt MATH - I. V. Sturova, โProblems of radiation and diffraction for a circular cylinder in a stratified fluid,โ
*Fluid Dynamics*, vol. 34, pp. 512โ533, 1999. View at: Google Scholar - D. Das and B. N. Mandal, โWater wave radiation by a sphere submerged in water with an ice-cover,โ
*Archive of Applied Mechanics*, vol. 78, no. 8, pp. 649โ661, 2008. View at: Publisher Site | Google Scholar - Y. X. You, Q. Shi, G. Wei, and G. P. Miao, โThe interaction of water waves with a vertically floating cylinder in a two-layer fluid,โ
*Journal of Shanghai Jiao Tong University*, vol. 41, no. 9, pp. 1460โ1464, 2007. View at: Google Scholar - B. Teng, G. Ying, N. Dezhi et al., โA higher-order BEM for wave-current action on structures direct computation of free-term coefficient and CPV integrals,โ
*China Ocean Engineering*, vol. 20, no. 3, pp. 395โ410, 2006. View at: Google Scholar - B. Teng and R. Eatock Taylor, โNew higher-order boundary element methods for wave diffraction/radiation,โ
*Applied Ocean Research*, vol. 17, pp. 71โ77, 1995. View at: Google Scholar - T. C. Nguyen and R. W. Yeung, โUnsteady three-dimensional sources for a two-layer fluid of finite depth and their applications,โ
*Journal of Engineering Mathematics*, vol. 70, no. 1–3, pp. 67โ91, 2011. View at: Publisher Site | Google Scholar

#### Copyright

Copyright © 2012 Ying Gou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.