# A Numerical Study of the Forces on Two Tandem Cylinders Exerted by Internal Solitary Waves

Authors: Yin Wang; Lingling Wang; Hai Zhu; Hongwu Tang; Gang Wei
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A three-dimensional numerical wave flume is employed to investigate the forces exerted by internal solitary waves (ISWs) on a pair of circular cylinders in tandem arrangement, using large-eddy simulation (LES) model. The effect of the centre-to-centre distance (L) ranging from 1.5 to 5 diameters (D) is studied for various ISWs amplitudes (η
a) in the two-layer fluid system. Vertical-averaged vorticity distribution and vertical-averaged pressure gradient distribution in each layer are presented to investigate the different hydrodynamic interference between cylinders and the ISWs forces on each cylinder at variousL. Furthermore, the force behaviors of the two cylinders are also compared with that of an isolated cylinder in the same environment. The interaction between the two piles occurs in both layers, and it is found that, for 1.5
≤
L
/
D
<
3.5, strong mutual interference appears between two cylinders; for 3.5
≤
L
/
D
<
5, the two cylinders continue to influence each other in a weak-interference state; for L
/
D
≥
5, the interaction gradually decreases into a noninteracting state. This paper tries to provide some references to structural arrangement of double-cylinder structure and grouped-cylinder structure in stratified flow environment.


---

## Body

## 1. Introduction

In oceans, estuaries, and lakes, the stable stratification of density happens while the fluid density changes along with the depth due to the variations of temperature, salinity, and other environmental factors. Internal solitary waves (ISWs) with different amplitudes may be produced by a tiny or weak perturbation in such stable stratified environment [1]. Based on the monitoring and measured data collected from the South China Sea, strong underwater currents caused by internal waves could be a serious threat to underwater structures, such as oil drilling platforms or supporting cylinders [2].Thus more and more researchers got involved in the studies of IWs action on such structures. Laboratory experiment is an important approach to investigate the ISWs loading on cylinders. Ermanyuk and Gavrilov [3] experimentally studied the hydrodynamic loads exerted by ISW on a submerged circular cylinder in a two-layer system and discovered the locations of the maximum and minimum horizontal loads on the cylinder. Wei et al. [4] manufactured a new wave-maker to excite the ISWs and developed a measurement technique of hydrodynamic load to determine the interaction characteristics between IWs and an isolated cylindrical body in the tank.Along with rapid development of computer and CFD technology, numerical methods have gotten more and more extensive application to investigate the loading induced by ISWs. The Morison empirical method [5], modal separation [6], and regression analyses [7] were applied to estimate the forces exerted by internal soliton on cylindrical piles. After that, a simplified estimation method of the force exerted by ISWs involving only the first mode internal wave was proposed [8]. The method can be used to estimate the force without observational current data and simplify the calculation procedure. For better solving the issues related to depression ISWs, an extended KdV model (EKdV) was employed to calculate the forces on a pile exerted by ISWs propagation from deep sea to shallow water [9]. Afterwards, Si et al. [10] used a generalized KdV model (GKdV) to obtain the vertical distribution of horizontal velocity of large-amplitude ISWs and discovered that the shear force was the largest at the turning point of the horizontal velocity. Based on water-wave theory, Linton and McIver [11] and Cadby and Linton [12] built a two-dimensional and three-dimensional numerical model, respectively, to study the interaction of waves with structures in two-layer fluids, and multipole expansions were used to solve the problems of wave radiation and scattering by a submerged structure in either the upper or lower layer. Then Sturova [13] conducted a study of radiation loads on interface piercing cylinder in a two-layer fluid of finite depth by a coupled element technique. By comparing with the ISWs force and the surface wave force, Du et al. [14] and Song et al. [15] found that the order magnitude of total force exerted by internal soliton was the same as that exerted by a surface wave, and the maximum total horizontal force caused by an internal wave amounted to 37.7% of that exerted by a surface wave [16]. Considering the above-mentioned, the ISWs force on underwater structures cannot be neglected.Nevertheless, most researches before just focused on the cases of ISWs forces on an isolated cylindrical pile. In practical engineering, pile-supported structures are generally arranged in tandem or side by side with various centre-to-centre distance (L) to form the so-called multiple slender structures, such as the bundle of risers linking the seabed to the offshore platforms and the piles being used for supporting bridges [17]. In general, two cylinders of the structures submerged in water behave in a similar manner to a single cylinder when the two cylinders are sufficiently apart. In some situations the cylinders need to be placed at a close proximity [18], and the interference between the two bodies significantly changes the flow around them [19]. Different flow patterns can be characterized by the behavior of the wake region [20], and unexpected flow structures and forces can be generated as the spacing between two circular cylinders changes [21].Up to date, very few numerical investigations about the internal solitary waves (ISWs) action on a pair of circular cylinders in tandem arrangements can be found in the literatures. Thus, this paper aims to study numerically ISWs force behaviors of two tandem cylinders (placed parallel to wave direction). To begin with, large-eddy simulation model (LES) is employed to simulate the generation and propagation of depression ISWs in a three-dimensional numerical wave flume. With respect to two-phase stratified flow system, the horizontal current induced by ISWs reverses the flow direction between upper and lower layers. Therefore, the hydrodynamic interference occurring in each layer at diverseL ranging from 1.5D to 5D will be investigated, respectively, from the perspective of vorticity distribution and pressure gradient distribution. The changes in vorticity fields and pressure gradient fields induced by different pile-to-pile interactions will be employed to study and explain the ISWs force behaviors of the two tandem cylinders. Finally, the results are also compared with that of an isolated cylinder in the same environment.


## 2. Theoretical Foundation

### 2.1. Korteweg-de Vries (KdV) Theory

Many nonlinear equations, represented by KdV equation, are widely applied to describe the propagation of the internal solitary waves in horizontal direction. KdV equation can be written as(1)
∂
η
∂
t
+
c
0
∂
η
∂
x
+
3
2
h
1
-
h
2
h
1
h
2
c
0
η
∂
η
∂
x
+
1
6
h
1
h
2
c
0
∂
3
η
∂
x
3
=
0
,where t is the time; h
1 and h
2 represent the thicknesses of the upper layer and the lower layer; the linear velocity of linear IWs c
o and the interfacial vertical displacement η can be, respectively, expressed as (2)
c
o
=
g
Δ
ρ
h
1
h
2
ρ
2
h
1
+
h
2
,
η
=
-
η
′
sech
2
⁡
x
-
c
p
t
l
,where η
′ is the amplitude of incident wave; c
p is the phase speed; l is the characteristic wavelengths; Δ
ρ is the density difference between two layers; g is the gravity acceleration. c
p and l read (3)
c
p
=
c
o
1
-
η
′
h
1
-
h
2
2
h
1
h
2
,
l
=
2
h
1
h
2
3
η
′
h
2
-
h
1
.


## 2.1. Korteweg-de Vries (KdV) Theory

Many nonlinear equations, represented by KdV equation, are widely applied to describe the propagation of the internal solitary waves in horizontal direction. KdV equation can be written as(1)
∂
η
∂
t
+
c
0
∂
η
∂
x
+
3
2
h
1
-
h
2
h
1
h
2
c
0
η
∂
η
∂
x
+
1
6
h
1
h
2
c
0
∂
3
η
∂
x
3
=
0
,where t is the time; h
1 and h
2 represent the thicknesses of the upper layer and the lower layer; the linear velocity of linear IWs c
o and the interfacial vertical displacement η can be, respectively, expressed as (2)
c
o
=
g
Δ
ρ
h
1
h
2
ρ
2
h
1
+
h
2
,
η
=
-
η
′
sech
2
⁡
x
-
c
p
t
l
,where η
′ is the amplitude of incident wave; c
p is the phase speed; l is the characteristic wavelengths; Δ
ρ is the density difference between two layers; g is the gravity acceleration. c
p and l read (3)
c
p
=
c
o
1
-
η
′
h
1
-
h
2
2
h
1
h
2
,
l
=
2
h
1
h
2
3
η
′
h
2
-
h
1
.


## 3. Numerical Model

### 3.1. Navier-Stokes Equations

The process of three-dimensional unsteady incompressible fluid motion governed by Navier-Stokes (N-S) equations can be described by(4)
∂
u
j
∂
x
j
=
0
,
∂
u
i
∂
t
+
u
j
∂
u
i
∂
x
j
=
-
1
ρ
∂
p
∂
x
i
+
v
∂
2
u
i
∂
x
j
∂
x
j
+
f
i
,where t is the time; u
j is the velocity component; x
j is the Cartesian coordinate; p is the pressure; v is the kinematic viscosity; and f
i is the body force that equals the gravity acceleration g in the vertical direction.


### 3.2. Scalar Transport Equation

The mass transfer between the two-layer water system has been taken into account. The scalar transport equation that governs the advection-diffusion effect is as follows:(5)
∂
C
∂
t
+
u
j
∂
C
∂
x
j
=
k
∂
2
C
∂
x
j
∂
x
j
,where C is the volume concentration of brine water in the lower layer of the fluid system;k is the molecular diffusivity coefficient.


### 3.3. LES Governing Equations

Applying a spatial filter to Navier-Stokes equations, the filtered momentum and mass equations can be written as(6)
∂
u
j
¯
∂
x
j
=
0
,
∂
u
i
¯
∂
t
+
u
j
¯
∂
u
i
¯
∂
x
j
=
-
1
ρ
∂
p
¯
∂
x
i
+
v
∂
2
u
i
¯
∂
x
j
∂
x
j
+
∂
τ
i
j
∂
x
i
+
f
i
¯
,
∂
C
¯
∂
t
+
u
j
¯
∂
C
¯
∂
x
j
=
k
∂
2
C
¯
∂
x
j
∂
x
j
+
∂
χ
j
∂
x
j
,where the overbar notation denotes the application of top-hat filter; the subgrid scale (SGS) stress tensor term τ
i
j
=
u
i
¯
u
j
¯
-
u
i
u
j
¯ in (9) is responsible for the momentum exchange between the subgrid scale and the resolved scale; χ
j
=
C
¯
u
j
¯
-
C
u
j
¯ is the subgrid scalar flux responsible for the scalar flux exchange between the subgrid scale and the resolved scale. Hence, subsequent modeling is required to determine the turbulent viscosity ν
t. As a function of the filter size and strain rate tensor, ν
t can be defined as (7)
ν
t
=
C
s
Δ
2
2
S
i
j
¯
S
i
j
¯
1
/
2
,where S
i
j
¯
=
1
/
2
∂
2
u
i
¯
/
∂
x
j
+
∂
2
u
j
¯
/
∂
x
j is the strain rate tensor; C
s is the Smagorinsky constant; Δ is the filtered width. In general, the model coefficient C
s varies in both time and space due to different flow conditions and even can be negative for a long time. Thus a dynamic procedure developed by Germano et al. [22] is employed to determine C
s.


### 3.4. Numerical Method and Boundary Condition

Large-eddy simulation model (LES) is employed to simulate the generation and propagation of ISWs of depression type. Velocity-pressure term is solved by SIMPLE algorithm to enforce mass conservation and to obtain the pressure field, while second-order centred differencing scheme is adopted for the spatial discretization, and the time step is discretized by second-order implicit scheme.The left boundary, belonging to the wave generation area where gravity collapse happens, and the two sidewalls and bottom of the wave tank are specified as a rigid wall with no-slip condition. Sommerfeld radiation type to avoid wave reflection is adopted to specify the right boundary, and the “rigid lid” approximation is used here to filter the free surface mode to ignore the influence of the surface wave [23].


### 3.5. Building of Numerical Wave Flume

A three-dimensional numerical wave flume established in current study is illustrated in Figure1. The numerical wave tank in this paper has a dimension of 12 m × 0.5 m × 0.4 m in the stream-wise (X), span-wise (Z), and vertical (Y) direction, which represents the length, width, and height, respectively. This flume is divided into two parts, for one is the wave generation area which locates L
o = 0.3 m from the left boundary in the X direction and the remaining is the wave propagation area. The two piles are placed in tandem arrangement. The bottom centre of the upstream cylindrical pile (Pile1) locates at (
x
,
y
,
z
)
=
(
5,0
,
0.25
) from the coordinate origin, whileL is the centre-to-centre distance between Pile1 and the downstream cylindrical pile (Pile2).D is the diameter of two piles. The density is ρ, and the thickness of each layer ish, of which the subscripts 1 and 2 represent the upper layer and the lower layer. The upper-layer fluid density ρ
1 is set to be 1000 kg/m3, the lower-layer fluid density ρ
2 is set to be 1030 kg/m3, and the volumetric concentrationC of the brine water in the lower layer is around 3%. The step height Δh we called here is the height difference of the pycnocline. The whole two-layer system keeps quiescent at the initial time.Figure 1
Schematic diagram showing the numerical model layout and the gravity collapse mechanism for ISW generation.Gravity collapse in a two-layer stratified fluid system is adopted to excite the ISWs [24]. Compared with the method of oscillating boundary or the movable body, the ISWs generation method of gravity collapse is simpler and widely employed in the laboratory experiments and numerical simulations [25]. An ISW of depression type forms here by making the height of the pycnocline in the wave propagation area higher than the middle depth, as shown in Figure 1.


### 3.6. Verification of the Numerical Model

The KdV equation is employed to validate the waveform firstly, then followed by a grid-size convergence analysis. The numerical simulation results of ISWs forces on an isolated pile are verified finally by comparing with the experiments in laboratory and the Morison equation.


#### 3.6.1. Verification of Wave Profiles

Two cases of waveform, with different wave amplitudesη
o, simulated by the numerical model have been compared with the KdV solutions. The corresponding parameters are shown in Table 1, where Δh is the step height, η
o is the amplitude of an ISW, and H
=
0.4 m is the total water depth in the tank.Table 1
Corresponding parameters of waveform verification cases.
Simulation number
Wave type
Δ
h
h
1 (m)
h
2 (m)
η
o/H
A1
Depression
0.05
0.1
0.3
0.062
A2
Depression
0.1
0.1
0.3
0.0985The total simulation time is 50 s for both verification cases and wave profiles are extracted att
=
33 s. Verification results show that the numerical simulations are in good agreement with the KdV equation, and waveform of smaller amplitude fits better with the solutions (see Figure 2).Figure 2
Comparison of numerical results with KdV solutions for wave profiles: (a) A1; (b) A2.
(a)
(b)


#### 3.6.2. Grid-Size Convergence Analysis

To assess the reliability of the proposed results, a grid independence study carried out on three different grids for an isolated circular cylinder is performed below. The difference in the grids is principally based on the number of nodes on the cylinder circumference and the total nodes in theY direction. Details of the grid independence test are given in Table 2, where η
a/H is the nondimension amplitude, Δt is the time step, and C
F
x
m
a
x represents the dimensionless global horizontal force amplitude exerted by internal solitary waves. The time step Δt can be chosen by the following expression [26]: C
r
=
c
p
Δ
t
/
Δ
x
<
1, Δ
t
<
Δ
x
/
c
p
≈
0.02, where Cr is the Courant number, c
p is the phase speed of internal solitary wave, and Δx is the minimum grid length in the X direction. Finally, we selected the time step Δ
t
=
0.01 s in this paper. Close-up views of the triangular unstructured mesh around an isolated cylinder of the three meshes are illustrated in Figure 3.Table 2
Details of the grid independence test carried out atη
a/H = 0.215.
Test number
Number of nodes on the cylinder circumference
Number of nodes in the vertical direction
Δ
t (s)
C
F
x
m
ax
T1 (low)
30
33
0.02
1.899
T2 (moderate)
60
65
0.01
2.258
T3 (high)
90
98
0.006
2.298Figure 3
Close-up view of the triangular unstructured mesh around a single circular cylinder: (a) T1: low density; (b) T2: moderate density; (c) T3: high density.
(a)
(b)
(c)Figure4 shows the computational results of three different mesh densities. It is observed that an increase in resolution from case T1 (low) to case T2 (moderate) gave an obvious difference in the magnitude of C
F
x; however, by comparing case T2 with case T3 (high), the curves of the nondimensional ISWs forces show that the difference between the two computations is very small. Hence, the computations are considered grid independent.Figure 4
Computational results of three different mesh densities, carried out atη
a
/
H = 0.21.


#### 3.6.3. Verification of Forces Acting on an Isolated Circular Pile

The numerical simulation results of ISWs forces on an isolated pile are verified by Morison equation and the experiments in laboratory.(1
) Verification by Morison Equation.  The Morison equation [27] can be adopted if the diameter of a pile,D, is supposed to be less than the wave length,L
o,D/L
o ≤ 0.15. In this investigation, the diameter of pile, D
=
0.05 m, is much smaller than the leading soliton length,L
o ≈ 2 m, makingD/L
o ≤ 0.15 become satisfied.The global forceF
m, made up of a drag component F
D and an inertial component F
I, on a cylindrical pile can be expressed as(8)
F
m
=
F
D
+
F
I
.F
D and F
I are determined by the following equations, respectively [14]:(9)
F
D
=
1
2
∫
0
H
ρ
g
C
D
D
u
u
d
z
,
(10)
F
I
=
∫
0
H
ρ
g
C
M
π
D
2
4
∂
u
∂
t
d
z
,where ρ is water density, g is gravity acceleration,D is the diameter of pile,u is the horizontal component of the current velocity,t is the time,z is the depth, and H is the total water depth in the tank. C
D is the drag-force coefficient, and C
M is the inertia-force coefficient. The values of C
D and C
M over a wave cycle vary with the current intensity and the size of the cylinders. The two coefficients, correlated with the Keulegan-Carpenter number (K
C) and the Reynolds number (Re), can be determined experimentally [28]. In this paper, K
C
=
u
m
a
x
T
/
D
=
26.5, R
e
=
u
m
a
x
D
/
ν
=
12000, where u
m
a
x is the maximum internal wave-induced velocity,T is the period of the internal wave,D is the diameter of the pile, and ν is the kinematic viscosity of fluid. According to the experimental results of Sarpkaya and Michael [29], C
D = 0.6, C
M = 1.8 are chosen in the calculation.(2
) Verification by Experiments.  Beyond the analysis performed by Morison equation, further validation of the numerical method is added by comparison with laboratory experiments.The experiments of ISWs force measurement are conducted in a long stratified fluid tank at the PLA University of Science and Technology, which is made of steel frame and glass materials with the dimension of 12 m length, 0.5 m width, and 0.5 m depth, as shown in Figure5. The experiment is implemented with initial conditions as close to those in the numerical setup as possible: the total water depth H in the tank is 0.4 m, while h
1
=
0.1 m, h
2
=
0.3 m, ρ
1
=
1000 kg/m3, ρ
2
=
1030 kg/m3, and the density difference Δ
ρ between the two layers is set to be 30 kg/m3. Referring to the wave-generating method used in numerical simulation, gravity collapse in a two-layer stratified fluid system is adopted to excite the ISWs in the experiments. The initial time of gravity collapse in the tank is shown in Figure 6.Figure 5
Schematic diagram of the stratified fluid test tank.Figure 6
Snapshot of the gravity collapse mechanism for ISW generation at the initial time.A slender cylinder is of 0.4 m length and 0.05 m diameter, the size of which is the same as the one used in the numerical simulation. The cylinder is placed at the centre of the tank in the span-wise direction and a distance of 5 meters from the right side. A force-measuring sensor with a small range and high precision and sensitivity is used to measure the force exerted by internal waves on the model. The sensor is connected to the test cylinder through an aluminum slender rod functioning as force transferring portion of the system; see Figure7. Then records of ISWs forces F
e versus time of different step heights and ISWs amplitudes can be obtained by the Spider8 processor made in Germany and the special software, which was of good reliability, high precision, and rapid speed. The principle of generating internal solitons and measuring forces in a stratified fluid tank is shown in Figure 8.Figure 7
Photo of the force measurement setup and ISW acting on the cylinder.Figure 8
Schematic diagram of generating internal solitons and measuring forces in a stratified fluid tank.Three nondimensional formulas of force are defined as follows:(11)
C
F
x
=
2
F
x
ρ
u
m
a
x
2
D
H
,
C
F
m
=
2
F
m
ρ
u
m
a
x
2
D
H
,
C
F
e
=
2
F
e
ρ
u
m
a
x
2
D
H
,where F
x denotes the global horizontal force calculated by the numerical model, F
m denotes the global horizontal force calculated by the Morison equation, F
e denotes the global horizontal force obtained by experiments, and u
m
a
x represents the maximum internal wave-induced velocity (positive in the direction of wave propagation from left to right when facing the wave flume). The computing time for each verification case is 50 s, and the corresponding parameters of two force verification cases are listed in Table 3, where η
a is the amplitude when the leading ISW just reaches the pile.Table 3
Corresponding parameters of force verification cases.
Simulation number
Wave type
Δ
h (m)
h
1 (m)
h
2 (m)
η
a/H
B1
Depression
0.15
0.1
0.3
0.183
B2
Depression
0.2
0.1
0.3
0.215The comparisons of numerical simulations with experiments and Morison equation for ISWs forcesC
F on the isolated pile versus time are given in Figure 9. A similar overall trend among the three approaches can be observed in the curves, and the experimental results C
F
e agree better with the numerical results C
F
x than that of Morison equation C
F
m. It should be pointed out that the positive direction of C
F defined in this paper is along the direction of ISWs propagation from left to right.Figure 9
Comparison of numerical simulations with experiments and Morison equation for ISWs forces on an isolated pile: (a) B1:η
a/H = 0.1825; (b) B2: η
a/H = 0.215.
(a)
(b)From the above three aspects (waveform, grid-size convergence, and wave force) of verification, it can be concluded that the numerical model is reliable to simulate the generation and propagation of ISWs of depression type and capable of calculating the force on a cylinder exerted by ISWs.


## 3.1. Navier-Stokes Equations

The process of three-dimensional unsteady incompressible fluid motion governed by Navier-Stokes (N-S) equations can be described by(4)
∂
u
j
∂
x
j
=
0
,
∂
u
i
∂
t
+
u
j
∂
u
i
∂
x
j
=
-
1
ρ
∂
p
∂
x
i
+
v
∂
2
u
i
∂
x
j
∂
x
j
+
f
i
,where t is the time; u
j is the velocity component; x
j is the Cartesian coordinate; p is the pressure; v is the kinematic viscosity; and f
i is the body force that equals the gravity acceleration g in the vertical direction.


## 3.2. Scalar Transport Equation

The mass transfer between the two-layer water system has been taken into account. The scalar transport equation that governs the advection-diffusion effect is as follows:(5)
∂
C
∂
t
+
u
j
∂
C
∂
x
j
=
k
∂
2
C
∂
x
j
∂
x
j
,where C is the volume concentration of brine water in the lower layer of the fluid system;k is the molecular diffusivity coefficient.


## 3.3. LES Governing Equations

Applying a spatial filter to Navier-Stokes equations, the filtered momentum and mass equations can be written as(6)
∂
u
j
¯
∂
x
j
=
0
,
∂
u
i
¯
∂
t
+
u
j
¯
∂
u
i
¯
∂
x
j
=
-
1
ρ
∂
p
¯
∂
x
i
+
v
∂
2
u
i
¯
∂
x
j
∂
x
j
+
∂
τ
i
j
∂
x
i
+
f
i
¯
,
∂
C
¯
∂
t
+
u
j
¯
∂
C
¯
∂
x
j
=
k
∂
2
C
¯
∂
x
j
∂
x
j
+
∂
χ
j
∂
x
j
,where the overbar notation denotes the application of top-hat filter; the subgrid scale (SGS) stress tensor term τ
i
j
=
u
i
¯
u
j
¯
-
u
i
u
j
¯ in (9) is responsible for the momentum exchange between the subgrid scale and the resolved scale; χ
j
=
C
¯
u
j
¯
-
C
u
j
¯ is the subgrid scalar flux responsible for the scalar flux exchange between the subgrid scale and the resolved scale. Hence, subsequent modeling is required to determine the turbulent viscosity ν
t. As a function of the filter size and strain rate tensor, ν
t can be defined as (7)
ν
t
=
C
s
Δ
2
2
S
i
j
¯
S
i
j
¯
1
/
2
,where S
i
j
¯
=
1
/
2
∂
2
u
i
¯
/
∂
x
j
+
∂
2
u
j
¯
/
∂
x
j is the strain rate tensor; C
s is the Smagorinsky constant; Δ is the filtered width. In general, the model coefficient C
s varies in both time and space due to different flow conditions and even can be negative for a long time. Thus a dynamic procedure developed by Germano et al. [22] is employed to determine C
s.


## 3.4. Numerical Method and Boundary Condition

Large-eddy simulation model (LES) is employed to simulate the generation and propagation of ISWs of depression type. Velocity-pressure term is solved by SIMPLE algorithm to enforce mass conservation and to obtain the pressure field, while second-order centred differencing scheme is adopted for the spatial discretization, and the time step is discretized by second-order implicit scheme.The left boundary, belonging to the wave generation area where gravity collapse happens, and the two sidewalls and bottom of the wave tank are specified as a rigid wall with no-slip condition. Sommerfeld radiation type to avoid wave reflection is adopted to specify the right boundary, and the “rigid lid” approximation is used here to filter the free surface mode to ignore the influence of the surface wave [23].


## 3.5. Building of Numerical Wave Flume

A three-dimensional numerical wave flume established in current study is illustrated in Figure1. The numerical wave tank in this paper has a dimension of 12 m × 0.5 m × 0.4 m in the stream-wise (X), span-wise (Z), and vertical (Y) direction, which represents the length, width, and height, respectively. This flume is divided into two parts, for one is the wave generation area which locates L
o = 0.3 m from the left boundary in the X direction and the remaining is the wave propagation area. The two piles are placed in tandem arrangement. The bottom centre of the upstream cylindrical pile (Pile1) locates at (
x
,
y
,
z
)
=
(
5,0
,
0.25
) from the coordinate origin, whileL is the centre-to-centre distance between Pile1 and the downstream cylindrical pile (Pile2).D is the diameter of two piles. The density is ρ, and the thickness of each layer ish, of which the subscripts 1 and 2 represent the upper layer and the lower layer. The upper-layer fluid density ρ
1 is set to be 1000 kg/m3, the lower-layer fluid density ρ
2 is set to be 1030 kg/m3, and the volumetric concentrationC of the brine water in the lower layer is around 3%. The step height Δh we called here is the height difference of the pycnocline. The whole two-layer system keeps quiescent at the initial time.Figure 1
Schematic diagram showing the numerical model layout and the gravity collapse mechanism for ISW generation.Gravity collapse in a two-layer stratified fluid system is adopted to excite the ISWs [24]. Compared with the method of oscillating boundary or the movable body, the ISWs generation method of gravity collapse is simpler and widely employed in the laboratory experiments and numerical simulations [25]. An ISW of depression type forms here by making the height of the pycnocline in the wave propagation area higher than the middle depth, as shown in Figure 1.


## 3.6. Verification of the Numerical Model

The KdV equation is employed to validate the waveform firstly, then followed by a grid-size convergence analysis. The numerical simulation results of ISWs forces on an isolated pile are verified finally by comparing with the experiments in laboratory and the Morison equation.


### 3.6.1. Verification of Wave Profiles

Two cases of waveform, with different wave amplitudesη
o, simulated by the numerical model have been compared with the KdV solutions. The corresponding parameters are shown in Table 1, where Δh is the step height, η
o is the amplitude of an ISW, and H
=
0.4 m is the total water depth in the tank.Table 1
Corresponding parameters of waveform verification cases.
Simulation number
Wave type
Δ
h
h
1 (m)
h
2 (m)
η
o/H
A1
Depression
0.05
0.1
0.3
0.062
A2
Depression
0.1
0.1
0.3
0.0985The total simulation time is 50 s for both verification cases and wave profiles are extracted att
=
33 s. Verification results show that the numerical simulations are in good agreement with the KdV equation, and waveform of smaller amplitude fits better with the solutions (see Figure 2).Figure 2
Comparison of numerical results with KdV solutions for wave profiles: (a) A1; (b) A2.
(a)
(b)


### 3.6.2. Grid-Size Convergence Analysis

To assess the reliability of the proposed results, a grid independence study carried out on three different grids for an isolated circular cylinder is performed below. The difference in the grids is principally based on the number of nodes on the cylinder circumference and the total nodes in theY direction. Details of the grid independence test are given in Table 2, where η
a/H is the nondimension amplitude, Δt is the time step, and C
F
x
m
a
x represents the dimensionless global horizontal force amplitude exerted by internal solitary waves. The time step Δt can be chosen by the following expression [26]: C
r
=
c
p
Δ
t
/
Δ
x
<
1, Δ
t
<
Δ
x
/
c
p
≈
0.02, where Cr is the Courant number, c
p is the phase speed of internal solitary wave, and Δx is the minimum grid length in the X direction. Finally, we selected the time step Δ
t
=
0.01 s in this paper. Close-up views of the triangular unstructured mesh around an isolated cylinder of the three meshes are illustrated in Figure 3.Table 2
Details of the grid independence test carried out atη
a/H = 0.215.
Test number
Number of nodes on the cylinder circumference
Number of nodes in the vertical direction
Δ
t (s)
C
F
x
m
ax
T1 (low)
30
33
0.02
1.899
T2 (moderate)
60
65
0.01
2.258
T3 (high)
90
98
0.006
2.298Figure 3
Close-up view of the triangular unstructured mesh around a single circular cylinder: (a) T1: low density; (b) T2: moderate density; (c) T3: high density.
(a)
(b)
(c)Figure4 shows the computational results of three different mesh densities. It is observed that an increase in resolution from case T1 (low) to case T2 (moderate) gave an obvious difference in the magnitude of C
F
x; however, by comparing case T2 with case T3 (high), the curves of the nondimensional ISWs forces show that the difference between the two computations is very small. Hence, the computations are considered grid independent.Figure 4
Computational results of three different mesh densities, carried out atη
a
/
H = 0.21.


### 3.6.3. Verification of Forces Acting on an Isolated Circular Pile

The numerical simulation results of ISWs forces on an isolated pile are verified by Morison equation and the experiments in laboratory.(1
) Verification by Morison Equation.  The Morison equation [27] can be adopted if the diameter of a pile,D, is supposed to be less than the wave length,L
o,D/L
o ≤ 0.15. In this investigation, the diameter of pile, D
=
0.05 m, is much smaller than the leading soliton length,L
o ≈ 2 m, makingD/L
o ≤ 0.15 become satisfied.The global forceF
m, made up of a drag component F
D and an inertial component F
I, on a cylindrical pile can be expressed as(8)
F
m
=
F
D
+
F
I
.F
D and F
I are determined by the following equations, respectively [14]:(9)
F
D
=
1
2
∫
0
H
ρ
g
C
D
D
u
u
d
z
,
(10)
F
I
=
∫
0
H
ρ
g
C
M
π
D
2
4
∂
u
∂
t
d
z
,where ρ is water density, g is gravity acceleration,D is the diameter of pile,u is the horizontal component of the current velocity,t is the time,z is the depth, and H is the total water depth in the tank. C
D is the drag-force coefficient, and C
M is the inertia-force coefficient. The values of C
D and C
M over a wave cycle vary with the current intensity and the size of the cylinders. The two coefficients, correlated with the Keulegan-Carpenter number (K
C) and the Reynolds number (Re), can be determined experimentally [28]. In this paper, K
C
=
u
m
a
x
T
/
D
=
26.5, R
e
=
u
m
a
x
D
/
ν
=
12000, where u
m
a
x is the maximum internal wave-induced velocity,T is the period of the internal wave,D is the diameter of the pile, and ν is the kinematic viscosity of fluid. According to the experimental results of Sarpkaya and Michael [29], C
D = 0.6, C
M = 1.8 are chosen in the calculation.(2
) Verification by Experiments.  Beyond the analysis performed by Morison equation, further validation of the numerical method is added by comparison with laboratory experiments.The experiments of ISWs force measurement are conducted in a long stratified fluid tank at the PLA University of Science and Technology, which is made of steel frame and glass materials with the dimension of 12 m length, 0.5 m width, and 0.5 m depth, as shown in Figure5. The experiment is implemented with initial conditions as close to those in the numerical setup as possible: the total water depth H in the tank is 0.4 m, while h
1
=
0.1 m, h
2
=
0.3 m, ρ
1
=
1000 kg/m3, ρ
2
=
1030 kg/m3, and the density difference Δ
ρ between the two layers is set to be 30 kg/m3. Referring to the wave-generating method used in numerical simulation, gravity collapse in a two-layer stratified fluid system is adopted to excite the ISWs in the experiments. The initial time of gravity collapse in the tank is shown in Figure 6.Figure 5
Schematic diagram of the stratified fluid test tank.Figure 6
Snapshot of the gravity collapse mechanism for ISW generation at the initial time.A slender cylinder is of 0.4 m length and 0.05 m diameter, the size of which is the same as the one used in the numerical simulation. The cylinder is placed at the centre of the tank in the span-wise direction and a distance of 5 meters from the right side. A force-measuring sensor with a small range and high precision and sensitivity is used to measure the force exerted by internal waves on the model. The sensor is connected to the test cylinder through an aluminum slender rod functioning as force transferring portion of the system; see Figure7. Then records of ISWs forces F
e versus time of different step heights and ISWs amplitudes can be obtained by the Spider8 processor made in Germany and the special software, which was of good reliability, high precision, and rapid speed. The principle of generating internal solitons and measuring forces in a stratified fluid tank is shown in Figure 8.Figure 7
Photo of the force measurement setup and ISW acting on the cylinder.Figure 8
Schematic diagram of generating internal solitons and measuring forces in a stratified fluid tank.Three nondimensional formulas of force are defined as follows:(11)
C
F
x
=
2
F
x
ρ
u
m
a
x
2
D
H
,
C
F
m
=
2
F
m
ρ
u
m
a
x
2
D
H
,
C
F
e
=
2
F
e
ρ
u
m
a
x
2
D
H
,where F
x denotes the global horizontal force calculated by the numerical model, F
m denotes the global horizontal force calculated by the Morison equation, F
e denotes the global horizontal force obtained by experiments, and u
m
a
x represents the maximum internal wave-induced velocity (positive in the direction of wave propagation from left to right when facing the wave flume). The computing time for each verification case is 50 s, and the corresponding parameters of two force verification cases are listed in Table 3, where η
a is the amplitude when the leading ISW just reaches the pile.Table 3
Corresponding parameters of force verification cases.
Simulation number
Wave type
Δ
h (m)
h
1 (m)
h
2 (m)
η
a/H
B1
Depression
0.15
0.1
0.3
0.183
B2
Depression
0.2
0.1
0.3
0.215The comparisons of numerical simulations with experiments and Morison equation for ISWs forcesC
F on the isolated pile versus time are given in Figure 9. A similar overall trend among the three approaches can be observed in the curves, and the experimental results C
F
e agree better with the numerical results C
F
x than that of Morison equation C
F
m. It should be pointed out that the positive direction of C
F defined in this paper is along the direction of ISWs propagation from left to right.Figure 9
Comparison of numerical simulations with experiments and Morison equation for ISWs forces on an isolated pile: (a) B1:η
a/H = 0.1825; (b) B2: η
a/H = 0.215.
(a)
(b)From the above three aspects (waveform, grid-size convergence, and wave force) of verification, it can be concluded that the numerical model is reliable to simulate the generation and propagation of ISWs of depression type and capable of calculating the force on a cylinder exerted by ISWs.


## 3.6.1. Verification of Wave Profiles

Two cases of waveform, with different wave amplitudesη
o, simulated by the numerical model have been compared with the KdV solutions. The corresponding parameters are shown in Table 1, where Δh is the step height, η
o is the amplitude of an ISW, and H
=
0.4 m is the total water depth in the tank.Table 1
Corresponding parameters of waveform verification cases.
Simulation number
Wave type
Δ
h
h
1 (m)
h
2 (m)
η
o/H
A1
Depression
0.05
0.1
0.3
0.062
A2
Depression
0.1
0.1
0.3
0.0985The total simulation time is 50 s for both verification cases and wave profiles are extracted att
=
33 s. Verification results show that the numerical simulations are in good agreement with the KdV equation, and waveform of smaller amplitude fits better with the solutions (see Figure 2).Figure 2
Comparison of numerical results with KdV solutions for wave profiles: (a) A1; (b) A2.
(a)
(b)


## 3.6.2. Grid-Size Convergence Analysis

To assess the reliability of the proposed results, a grid independence study carried out on three different grids for an isolated circular cylinder is performed below. The difference in the grids is principally based on the number of nodes on the cylinder circumference and the total nodes in theY direction. Details of the grid independence test are given in Table 2, where η
a/H is the nondimension amplitude, Δt is the time step, and C
F
x
m
a
x represents the dimensionless global horizontal force amplitude exerted by internal solitary waves. The time step Δt can be chosen by the following expression [26]: C
r
=
c
p
Δ
t
/
Δ
x
<
1, Δ
t
<
Δ
x
/
c
p
≈
0.02, where Cr is the Courant number, c
p is the phase speed of internal solitary wave, and Δx is the minimum grid length in the X direction. Finally, we selected the time step Δ
t
=
0.01 s in this paper. Close-up views of the triangular unstructured mesh around an isolated cylinder of the three meshes are illustrated in Figure 3.Table 2
Details of the grid independence test carried out atη
a/H = 0.215.
Test number
Number of nodes on the cylinder circumference
Number of nodes in the vertical direction
Δ
t (s)
C
F
x
m
ax
T1 (low)
30
33
0.02
1.899
T2 (moderate)
60
65
0.01
2.258
T3 (high)
90
98
0.006
2.298Figure 3
Close-up view of the triangular unstructured mesh around a single circular cylinder: (a) T1: low density; (b) T2: moderate density; (c) T3: high density.
(a)
(b)
(c)Figure4 shows the computational results of three different mesh densities. It is observed that an increase in resolution from case T1 (low) to case T2 (moderate) gave an obvious difference in the magnitude of C
F
x; however, by comparing case T2 with case T3 (high), the curves of the nondimensional ISWs forces show that the difference between the two computations is very small. Hence, the computations are considered grid independent.Figure 4
Computational results of three different mesh densities, carried out atη
a
/
H = 0.21.


## 3.6.3. Verification of Forces Acting on an Isolated Circular Pile

The numerical simulation results of ISWs forces on an isolated pile are verified by Morison equation and the experiments in laboratory.(1
) Verification by Morison Equation.  The Morison equation [27] can be adopted if the diameter of a pile,D, is supposed to be less than the wave length,L
o,D/L
o ≤ 0.15. In this investigation, the diameter of pile, D
=
0.05 m, is much smaller than the leading soliton length,L
o ≈ 2 m, makingD/L
o ≤ 0.15 become satisfied.The global forceF
m, made up of a drag component F
D and an inertial component F
I, on a cylindrical pile can be expressed as(8)
F
m
=
F
D
+
F
I
.F
D and F
I are determined by the following equations, respectively [14]:(9)
F
D
=
1
2
∫
0
H
ρ
g
C
D
D
u
u
d
z
,
(10)
F
I
=
∫
0
H
ρ
g
C
M
π
D
2
4
∂
u
∂
t
d
z
,where ρ is water density, g is gravity acceleration,D is the diameter of pile,u is the horizontal component of the current velocity,t is the time,z is the depth, and H is the total water depth in the tank. C
D is the drag-force coefficient, and C
M is the inertia-force coefficient. The values of C
D and C
M over a wave cycle vary with the current intensity and the size of the cylinders. The two coefficients, correlated with the Keulegan-Carpenter number (K
C) and the Reynolds number (Re), can be determined experimentally [28]. In this paper, K
C
=
u
m
a
x
T
/
D
=
26.5, R
e
=
u
m
a
x
D
/
ν
=
12000, where u
m
a
x is the maximum internal wave-induced velocity,T is the period of the internal wave,D is the diameter of the pile, and ν is the kinematic viscosity of fluid. According to the experimental results of Sarpkaya and Michael [29], C
D = 0.6, C
M = 1.8 are chosen in the calculation.(2
) Verification by Experiments.  Beyond the analysis performed by Morison equation, further validation of the numerical method is added by comparison with laboratory experiments.The experiments of ISWs force measurement are conducted in a long stratified fluid tank at the PLA University of Science and Technology, which is made of steel frame and glass materials with the dimension of 12 m length, 0.5 m width, and 0.5 m depth, as shown in Figure5. The experiment is implemented with initial conditions as close to those in the numerical setup as possible: the total water depth H in the tank is 0.4 m, while h
1
=
0.1 m, h
2
=
0.3 m, ρ
1
=
1000 kg/m3, ρ
2
=
1030 kg/m3, and the density difference Δ
ρ between the two layers is set to be 30 kg/m3. Referring to the wave-generating method used in numerical simulation, gravity collapse in a two-layer stratified fluid system is adopted to excite the ISWs in the experiments. The initial time of gravity collapse in the tank is shown in Figure 6.Figure 5
Schematic diagram of the stratified fluid test tank.Figure 6
Snapshot of the gravity collapse mechanism for ISW generation at the initial time.A slender cylinder is of 0.4 m length and 0.05 m diameter, the size of which is the same as the one used in the numerical simulation. The cylinder is placed at the centre of the tank in the span-wise direction and a distance of 5 meters from the right side. A force-measuring sensor with a small range and high precision and sensitivity is used to measure the force exerted by internal waves on the model. The sensor is connected to the test cylinder through an aluminum slender rod functioning as force transferring portion of the system; see Figure7. Then records of ISWs forces F
e versus time of different step heights and ISWs amplitudes can be obtained by the Spider8 processor made in Germany and the special software, which was of good reliability, high precision, and rapid speed. The principle of generating internal solitons and measuring forces in a stratified fluid tank is shown in Figure 8.Figure 7
Photo of the force measurement setup and ISW acting on the cylinder.Figure 8
Schematic diagram of generating internal solitons and measuring forces in a stratified fluid tank.Three nondimensional formulas of force are defined as follows:(11)
C
F
x
=
2
F
x
ρ
u
m
a
x
2
D
H
,
C
F
m
=
2
F
m
ρ
u
m
a
x
2
D
H
,
C
F
e
=
2
F
e
ρ
u
m
a
x
2
D
H
,where F
x denotes the global horizontal force calculated by the numerical model, F
m denotes the global horizontal force calculated by the Morison equation, F
e denotes the global horizontal force obtained by experiments, and u
m
a
x represents the maximum internal wave-induced velocity (positive in the direction of wave propagation from left to right when facing the wave flume). The computing time for each verification case is 50 s, and the corresponding parameters of two force verification cases are listed in Table 3, where η
a is the amplitude when the leading ISW just reaches the pile.Table 3
Corresponding parameters of force verification cases.
Simulation number
Wave type
Δ
h (m)
h
1 (m)
h
2 (m)
η
a/H
B1
Depression
0.15
0.1
0.3
0.183
B2
Depression
0.2
0.1
0.3
0.215The comparisons of numerical simulations with experiments and Morison equation for ISWs forcesC
F on the isolated pile versus time are given in Figure 9. A similar overall trend among the three approaches can be observed in the curves, and the experimental results C
F
e agree better with the numerical results C
F
x than that of Morison equation C
F
m. It should be pointed out that the positive direction of C
F defined in this paper is along the direction of ISWs propagation from left to right.Figure 9
Comparison of numerical simulations with experiments and Morison equation for ISWs forces on an isolated pile: (a) B1:η
a/H = 0.1825; (b) B2: η
a/H = 0.215.
(a)
(b)From the above three aspects (waveform, grid-size convergence, and wave force) of verification, it can be concluded that the numerical model is reliable to simulate the generation and propagation of ISWs of depression type and capable of calculating the force on a cylinder exerted by ISWs.


## 4. Results and Discussion

The flow field largely depends on Reynolds number, range of which covers7800
≤
R
e
=
u
m
a
x
D
/
v
≤
12000 for different dimensionless ISWs amplitudes η
a/H ranging from 0.12 to 0.215, when it comes to flow around a single circular cylinder [30], while as to the case of two tandem cylinders, the dimensionless centre-to-centre distance L
/
D also plays a decisive role [31, 32]. In order to investigate the forces on a multiple slender structure and the hydrodynamic interference between the cylinders, simulations have been carried out for two cylinders in tandem arrangements.To begin with, note that in single-layer flow environment, flows around two circular cylinders in tandem arrangements can be grouped into two main categories: with or without mutual interference [22]. Furthermore, different interference regions in the “mutual interference” category were identified by the classification of Igarashi [33]. Then the pressure and forces on two tandem cylinders in these regions can be analyzed by employing the vortex method [34].With respect to two-phase stratified flow system, as seen in Figure10, while the direction of the ISWs propagation is from left to right, the horizontal current induced by internal waves reverses the flow direction between upper and lower layers. The interaction between Pile1 and Pile2 occurs in both layers. Therefore, the physical mechanisms of the interference between the two cylinders should be studied separately in each layer, which acts as an important method contributing to investigating the direction and the effect of ISWs loading on a multiple slender structure.Figure 10
Instantaneous streamlines profiles showing the wave direction of two layers.Totally 45 cases are simulated to investigate the ISWs forces on two tandem piles for the range of Δh from 0.1 m to 0.2 m and L
/
D from 1.5 to 5. The typical simulation conditions of Δh = 0.2 are listed in Table 4, where η
a is the wave amplitude when the leading ISW just reaches Pile1. The hydrodynamic interference induced by ISWs occurring in each layer at diverseL will be, respectively, discussed in detail below from the perspective of vorticity distribution and pressure gradient distribution [21]. The analysis of flow field characteristic and the forces of Pile1 will be discussed in Sections 4.1 and 4.2 and then followed by those of Pile2 in Sections 4.3 and 4.4.Table 4
Typical simulation conditions.
Simulation number
Wave type
Δ
h (m)
η
a
/
H
L
/
D
C1
Depression
0.1
0.12
1.5
C2
Depression
0.2
0.215
1.5
C3
Depression
0.2
0.215
2
C4
Depression
0.2
0.215
2.5
C5
Depression
0.2
0.215
3.5
C6
Depression
0.2
0.215
5


### 4.1. Analysis of the Flow Field Characteristic around Pile1 for Various Tandem Arrangements

(a) At the gap of 1.5D: in the upper layer, the plots in Figure 11 show the contours of instantaneous vertical-averaged pressure gradient distribution of the two tandem piles and the isolated pile in the upper layer, where Zone A and Zone B denote the favorable pressure gradient area and the adverse pressure gradient area (all the instantaneous contours shown in this paper are adopted when C
F
x reaches its maximum for each simulation case [4]). By comparing the contours of L
/
D
=
1.5 in Figure 11(a) with that of an isolated pile in Figure 11(b), one can clearly observe how the adverse pressure gradient (Zone B) behind Pile1 is influenced and suppressed by the presence of the downstream body at the gap of 1.5D; in the lower layer, Figure 12(a) shows the instantaneous vertical-averaged vorticity contours with streamlines. Contours are extracted when C
F
x is the maximum. It is important to note that Pile1 is inside the near wake and immersed in a low-pressure region formed by the separated shear layers emanating from Pile2.Figure 11
Instantaneous vertical-averaged pressure gradient contours whenC
F
x is the maximum, η
a
/
H
=
0.215: (a) L
/
D
=
1.5, (b) isolated pile, and (c) L
/
D
=
2, in the upper layer.
(a)
(b)
(c)Figure 12
Instantaneous vertical-averaged vorticity contours and streamlines for diverseL
/
D when C
F
x is the maximum, η
a/H = 0.215: (a) L
/
D
=
1.5, (b) L
/
D
=
2, (c) L
/
D
=
3.5, and (d) L
/
D
=
5, in the lower layer.
(a)
(b)
(c)
(d)(b) At the gap of 2D: in the upper layer, Zone A and Zone B of L
/
D
=
2 shown in Figure 11(c) are very similar to those in Figure 11(b). It indicates that the effect of pressure gradient on Pile1 will not work when the gap is greater than 2D. Thus only the vortex effect in the lower layer is needed to be given consideration for gaps ≥2D for Pile1; in the lower layer, as illustrated in Figure 12(b), the rear of Pile1 begins to be hit by the vortices emanated from Pile2 instead of being immersed in the low-pressure region formed in front of Pile2 [35].(c) At the gap of 3.5D: in the lower layer, vortices generated in front of Pile2 continue to influence the horizontal forces on Pile1 after the strong mutual interference region (spacing ≤ 2D); see Figure 12(c). The impact force exerted by vortices decreases with the increase of the gap.(d) At the gap of 5D: in the lower layer, the interference between the two piles gradually decays into a noninteracting state presenting between the two piles in Figure 12(d).Therefore, on the basis of the above analysis, the interference regions between two cylinders for Pile1 can be preliminarily classified into three interference regions: 1.5
≤
L
/
D
<
2, 2
≤
L
/
D
<
3.5, and 3.5
≤
L
/
D
≤
5. And the changes of vorticity fields and pressure gradient fields induced by different hydrodynamic interactions in these regions will be employed to study and explain the ISWs force behaviors of Pile1 in the following.


### 4.2. ISWs Forces on Pile1

Curves ofC
F
x
m
a
x changing with L
/
D for different η
a/H of Pile1 are shown in Figure 13, where C
F
x
m
a
x represents the dimensionless global horizontal force amplitude for each simulation case, and η
a/H denotes the dimensionless wave amplitude. First of all, it should be noted that the positive direction of all forces defined in this study coincides with the direction of ISWs propagation from left to right, and the negative direction is from right to left. Besides, two physical quantities f
1 and f
2 are defined here to simplify the descriptions of the interaction force between two piles, as illustrated in Figure 14, where f
1 denotes the horizontal force acting on the rear of Pile1 in the negative direction; f
2 denotes the horizontal force acting on the front of Pile2 in the positive direction.Figure 13
Variations ofC
F
x
m
a
x with L
/
D for different η
a/H of Pile1.Figure 14
Schematic diagram showing the interaction forces between two piles.As it can be seen in Figure13, the overall trends of all the five curves of C
F
x
m
a
x on Pile1 are the same in respect to the variations of L
/
D. Below the spacing of 2D, it can be found that C
F
x
m
a
x acting on Pile1 reaches its minimum at the gap of 2D after experiencing the maximum at the gap of 1.5D.At the gap of 1.5D in the upper layer, the drop of the adverse pressure gradient (by comparing Figure 11(a) with Figure 11(b)) reduces f
1, enlarging C
F
x as well as its amplitude C
F
x
m
a
x; meanwhile, in the lower layer, the low-pressure region formed in front of Pile2 also reduces f
1 [20]. So C
F
x
m
a
x increases in both layers and meets its maximum at 1.5D.At the gap of 2D, the influence of the pressure gradient on C
F
x
m
a
x (mentioned in Section 4.1(b)) being ignored, only the vortex effect in the lower layer should be taken into account. As illustrated in Figure 12(b), Pile1 is no longer influenced by the low-pressure region and turns to be hit by the impact force exerted by vortices (in the negative direction), making f
1 larger. So C
F
x
m
a
x will be decreased by being counteracted by f
1 and thus reach its valley value at this gap.In the next stage when the intervals are in the range of2
<
L
/
D
<
3.5, the curves of all the five configurations present a positive slope; see in Figure 13. A continued growth of C
F
x
m
a
x can be found in this region. It is because the impact force of vortices on the rear of Pile1 keeps on reducing as the distance increases further.If the gap is gradually increased to 5D, further increase of L brings little change of  C
F
x
m
a
x above the spacing of 3.5D for all configurations. For a better understanding, C
F
x
m
a
x changing with Δh/H of the two tandem piles is compared with that of an isolated cylinder at this gap, as shown in Figure 15(c), where the forces on the two tandem cylinders are in a similar way to that on the single cylinder. It is indicated that the two-pile interaction will gradually disappear when L
/
D
>
5.Figure 15
Dependence onη
a/H of C
F
x
m
a
x on the tandem piles compared with isolated one: (a) L
/
D
=
1.5, (b) L
/
D
=
3.5, and (c) L
/
D
=
5.
(a)
(b)
(c)A point worth noting is that a shape of “kink” can be found below the gap of 2.5D, as illustrated in Figure 13, and the “kink” becomes more pronounced as η
a/H gets larger. See Line 1 in Figure 15(a); the trend graph of C
F
x
m
a
x with η
a/H shows that the increase of ISWs amplitude leads to a rapid rise in horizontal force on upstream pile at the gap of 1.5D. It can be explained by the fact that the vortex intensity between two piles for the configuration of η
a/H = 0.215 is significantly stronger than that of η
a/H = 0.12, evidenced by comparing the vorticity contours in Figure 16(a) with that in Figure 16(b). The comparison of the two figures also reveals that the case of η
a/H = 0.215 can provide clearer contours diagrams for visualizing when compared with the case of smaller wave amplitude. Since the forces trend curves of Pile1 show a similar behavior, the configuration of Δh = 0.2 m (η
a/H = 0.215) can be taken as a typical condition to elaborate and analyze the force behaviors of piles here. The contours of other η
a/H configurations will not be given in this paper.Figure 16
Instantaneous vertical-averaged vorticity contours with streamlines for diverseη
a
/
H when C
F
x is the maximum, L
/
D
=
1.5: (a) η
a
/
H
=
0.12, (b) η
a
/
H
=
0.215, in the lower layer.
(a)
(b)


### 4.3. Analysis of the Flow Field Characteristic around Pile2 for Various Tandem Arrangements

(a) At the gap of 1.5D: Figure 17(a) illustrates the instantaneous vertical-averaged vorticity contours with streamlines in the upper layer. It can be noticed that Pile2 is inside the wake formed behind Pile1, making the downstream body be immersed in a region of low pressure; in the lower layer, by comparing the contours of L
/
D
=
1.5 with a single pile as illustrated in Figures 18(a) and 18(b), we find that the instantaneous vertical-averaged adverse pressure gradient (Zone B) in front of Pile2 is obviously weakened by the existence of Pile1.Figure 17
Instantaneous vertical-averaged vorticity contours and streamlines for diverseL
/
D when C
F
x is the maximum, η
a/H = 0.5: (a) L
/
D
=
1.5, (b) L
/
D
=
2.5, (c) L
/
D
=
3.5, and (d) L
/
D
=
5, in the upper layer.
(a)
(b)
(c)
(d)Figure 18
Instantaneous vertical-averaged pressure gradient contours whenC
F
x is the maximum, η
a
=
0.215: (a) L
/
D
=
1.5, (b) isolated pile, and (c) L
/
D
=
2.5, in the lower layer.
(a)
(b)
(c)(b) At the gap of 2.5D: in the upper layer, it can be seen clearly in Figure 17(b) that Pile2 moves out of the lower-pressure and recirculation area, continuing to be hit by the vortices generated behind Pile1; in the lower layer, Zone A and Zone B of L
/
D
=
2.5 shown in Figure 18(c) are in a similar situation to that of the single pile in Figure 18(b), which reveals that the effect of pressure gradient on Pile2 can be ignored at this gap. Therefore, only the vortex effect in the upper layer should be taken into account for gaps ≥ 2.5.(c) At the gap of 3.5D: in the upper layer, the hydrodynamic interaction still exists between two piles (see Figure 17(c)). But the hitting force exerted by vortices decreases continuously with the spacing increasing accordingly.(d) At the gap of 5D: plot of instantaneous vertical-averaged vorticity contours in Figure 17(d) shows that the interaction between two piles tends to vanish when the spacing is large enough.Similarly, the spacing between two cylinders for Pile2 can also be classified into three interference regions distinguished by different hydrodynamic interactions: 1.5
≤
L
/
D
<
2.5, 2.5
≤
L
/
D
<
3.5, and 3.5
≤
L
/
D
≤
5. The force behaviors of Pile2 in these regions will also be further discussed below.


### 4.4. ISWs Forces on Pile2

The three-dimensional simulation results ofC
F
x
m
a
x on Pile2 are shown in Figure 19. It can be seen that all the five configurations exhibit a similar behavior over L
/
D, except that the systematic differences in the values of C
F
x
m
a
x are observed along the entire L
/
D range investigated. By comparing the following five curves of the measured C
F
x
m
a
x of Pile2, the systematic difference between two adjacent curves diminishes gradually with the increase of η
a/H. It indicates that the influence by wave amplitude on horizontal forces tends to depression. Similar phenomenon can also be found in Figure 13 of  C
F
x
m
a
x on Pile1.Figure 19
C
F
x
m
a
x as a function of L
/
D for different η
a/H of Pile2.For gaps in the range of1.5
≤
L
/
D
≤
2.5, as illustrated in Figure 19, C
F
x
m
a
x of each case increases rapidly with the increasing spacing. First, C
F
x
m
a
x meets its minimum at the gap of 1.5D. This fact can be understood by the weakening of the adverse pressure gradient and the influence of lower-pressure area (studied in Section 4.3(a)), which makes f
2 (defined in Figure 14) on Pile2 decrease in both layers at this gap. When the gap reaches 2.5D (see in Figure 17(b)) Pile2’s moving out of the lower-pressure area and continuing to be hit by vortices induce rapid increase in f
2 and cause C
F
x
m
a
x to increase sharply to its peak. Afterwards, the impact force by vortices continually decreases as the spacing is being enlarged to 3.5D. So the curves for all the configurations present a descending trend between the gap of 2.5D and 3.5D.If the gap is further enlarged to 5D, the variation of C
F
x
m
a
x with cylinder centre-to-centre spacing is very small as shown in Figure 19. Plot of instantaneous vertical-averaged vorticity contours (see in Figure 14(d)) shows that interaction between two cylinders at the gap of 5D gradually decreases into a noninteracting state. Therefore, it is clear that, at this gap, the horizontal force on the multiple slender structure is altered in a similar way to that observed in the isolated cylinder case.


## 4.1. Analysis of the Flow Field Characteristic around Pile1 for Various Tandem Arrangements

(a) At the gap of 1.5D: in the upper layer, the plots in Figure 11 show the contours of instantaneous vertical-averaged pressure gradient distribution of the two tandem piles and the isolated pile in the upper layer, where Zone A and Zone B denote the favorable pressure gradient area and the adverse pressure gradient area (all the instantaneous contours shown in this paper are adopted when C
F
x reaches its maximum for each simulation case [4]). By comparing the contours of L
/
D
=
1.5 in Figure 11(a) with that of an isolated pile in Figure 11(b), one can clearly observe how the adverse pressure gradient (Zone B) behind Pile1 is influenced and suppressed by the presence of the downstream body at the gap of 1.5D; in the lower layer, Figure 12(a) shows the instantaneous vertical-averaged vorticity contours with streamlines. Contours are extracted when C
F
x is the maximum. It is important to note that Pile1 is inside the near wake and immersed in a low-pressure region formed by the separated shear layers emanating from Pile2.Figure 11
Instantaneous vertical-averaged pressure gradient contours whenC
F
x is the maximum, η
a
/
H
=
0.215: (a) L
/
D
=
1.5, (b) isolated pile, and (c) L
/
D
=
2, in the upper layer.
(a)
(b)
(c)Figure 12
Instantaneous vertical-averaged vorticity contours and streamlines for diverseL
/
D when C
F
x is the maximum, η
a/H = 0.215: (a) L
/
D
=
1.5, (b) L
/
D
=
2, (c) L
/
D
=
3.5, and (d) L
/
D
=
5, in the lower layer.
(a)
(b)
(c)
(d)(b) At the gap of 2D: in the upper layer, Zone A and Zone B of L
/
D
=
2 shown in Figure 11(c) are very similar to those in Figure 11(b). It indicates that the effect of pressure gradient on Pile1 will not work when the gap is greater than 2D. Thus only the vortex effect in the lower layer is needed to be given consideration for gaps ≥2D for Pile1; in the lower layer, as illustrated in Figure 12(b), the rear of Pile1 begins to be hit by the vortices emanated from Pile2 instead of being immersed in the low-pressure region formed in front of Pile2 [35].(c) At the gap of 3.5D: in the lower layer, vortices generated in front of Pile2 continue to influence the horizontal forces on Pile1 after the strong mutual interference region (spacing ≤ 2D); see Figure 12(c). The impact force exerted by vortices decreases with the increase of the gap.(d) At the gap of 5D: in the lower layer, the interference between the two piles gradually decays into a noninteracting state presenting between the two piles in Figure 12(d).Therefore, on the basis of the above analysis, the interference regions between two cylinders for Pile1 can be preliminarily classified into three interference regions: 1.5
≤
L
/
D
<
2, 2
≤
L
/
D
<
3.5, and 3.5
≤
L
/
D
≤
5. And the changes of vorticity fields and pressure gradient fields induced by different hydrodynamic interactions in these regions will be employed to study and explain the ISWs force behaviors of Pile1 in the following.


## 4.2. ISWs Forces on Pile1

Curves ofC
F
x
m
a
x changing with L
/
D for different η
a/H of Pile1 are shown in Figure 13, where C
F
x
m
a
x represents the dimensionless global horizontal force amplitude for each simulation case, and η
a/H denotes the dimensionless wave amplitude. First of all, it should be noted that the positive direction of all forces defined in this study coincides with the direction of ISWs propagation from left to right, and the negative direction is from right to left. Besides, two physical quantities f
1 and f
2 are defined here to simplify the descriptions of the interaction force between two piles, as illustrated in Figure 14, where f
1 denotes the horizontal force acting on the rear of Pile1 in the negative direction; f
2 denotes the horizontal force acting on the front of Pile2 in the positive direction.Figure 13
Variations ofC
F
x
m
a
x with L
/
D for different η
a/H of Pile1.Figure 14
Schematic diagram showing the interaction forces between two piles.As it can be seen in Figure13, the overall trends of all the five curves of C
F
x
m
a
x on Pile1 are the same in respect to the variations of L
/
D. Below the spacing of 2D, it can be found that C
F
x
m
a
x acting on Pile1 reaches its minimum at the gap of 2D after experiencing the maximum at the gap of 1.5D.At the gap of 1.5D in the upper layer, the drop of the adverse pressure gradient (by comparing Figure 11(a) with Figure 11(b)) reduces f
1, enlarging C
F
x as well as its amplitude C
F
x
m
a
x; meanwhile, in the lower layer, the low-pressure region formed in front of Pile2 also reduces f
1 [20]. So C
F
x
m
a
x increases in both layers and meets its maximum at 1.5D.At the gap of 2D, the influence of the pressure gradient on C
F
x
m
a
x (mentioned in Section 4.1(b)) being ignored, only the vortex effect in the lower layer should be taken into account. As illustrated in Figure 12(b), Pile1 is no longer influenced by the low-pressure region and turns to be hit by the impact force exerted by vortices (in the negative direction), making f
1 larger. So C
F
x
m
a
x will be decreased by being counteracted by f
1 and thus reach its valley value at this gap.In the next stage when the intervals are in the range of2
<
L
/
D
<
3.5, the curves of all the five configurations present a positive slope; see in Figure 13. A continued growth of C
F
x
m
a
x can be found in this region. It is because the impact force of vortices on the rear of Pile1 keeps on reducing as the distance increases further.If the gap is gradually increased to 5D, further increase of L brings little change of  C
F
x
m
a
x above the spacing of 3.5D for all configurations. For a better understanding, C
F
x
m
a
x changing with Δh/H of the two tandem piles is compared with that of an isolated cylinder at this gap, as shown in Figure 15(c), where the forces on the two tandem cylinders are in a similar way to that on the single cylinder. It is indicated that the two-pile interaction will gradually disappear when L
/
D
>
5.Figure 15
Dependence onη
a/H of C
F
x
m
a
x on the tandem piles compared with isolated one: (a) L
/
D
=
1.5, (b) L
/
D
=
3.5, and (c) L
/
D
=
5.
(a)
(b)
(c)A point worth noting is that a shape of “kink” can be found below the gap of 2.5D, as illustrated in Figure 13, and the “kink” becomes more pronounced as η
a/H gets larger. See Line 1 in Figure 15(a); the trend graph of C
F
x
m
a
x with η
a/H shows that the increase of ISWs amplitude leads to a rapid rise in horizontal force on upstream pile at the gap of 1.5D. It can be explained by the fact that the vortex intensity between two piles for the configuration of η
a/H = 0.215 is significantly stronger than that of η
a/H = 0.12, evidenced by comparing the vorticity contours in Figure 16(a) with that in Figure 16(b). The comparison of the two figures also reveals that the case of η
a/H = 0.215 can provide clearer contours diagrams for visualizing when compared with the case of smaller wave amplitude. Since the forces trend curves of Pile1 show a similar behavior, the configuration of Δh = 0.2 m (η
a/H = 0.215) can be taken as a typical condition to elaborate and analyze the force behaviors of piles here. The contours of other η
a/H configurations will not be given in this paper.Figure 16
Instantaneous vertical-averaged vorticity contours with streamlines for diverseη
a
/
H when C
F
x is the maximum, L
/
D
=
1.5: (a) η
a
/
H
=
0.12, (b) η
a
/
H
=
0.215, in the lower layer.
(a)
(b)


## 4.3. Analysis of the Flow Field Characteristic around Pile2 for Various Tandem Arrangements

(a) At the gap of 1.5D: Figure 17(a) illustrates the instantaneous vertical-averaged vorticity contours with streamlines in the upper layer. It can be noticed that Pile2 is inside the wake formed behind Pile1, making the downstream body be immersed in a region of low pressure; in the lower layer, by comparing the contours of L
/
D
=
1.5 with a single pile as illustrated in Figures 18(a) and 18(b), we find that the instantaneous vertical-averaged adverse pressure gradient (Zone B) in front of Pile2 is obviously weakened by the existence of Pile1.Figure 17
Instantaneous vertical-averaged vorticity contours and streamlines for diverseL
/
D when C
F
x is the maximum, η
a/H = 0.5: (a) L
/
D
=
1.5, (b) L
/
D
=
2.5, (c) L
/
D
=
3.5, and (d) L
/
D
=
5, in the upper layer.
(a)
(b)
(c)
(d)Figure 18
Instantaneous vertical-averaged pressure gradient contours whenC
F
x is the maximum, η
a
=
0.215: (a) L
/
D
=
1.5, (b) isolated pile, and (c) L
/
D
=
2.5, in the lower layer.
(a)
(b)
(c)(b) At the gap of 2.5D: in the upper layer, it can be seen clearly in Figure 17(b) that Pile2 moves out of the lower-pressure and recirculation area, continuing to be hit by the vortices generated behind Pile1; in the lower layer, Zone A and Zone B of L
/
D
=
2.5 shown in Figure 18(c) are in a similar situation to that of the single pile in Figure 18(b), which reveals that the effect of pressure gradient on Pile2 can be ignored at this gap. Therefore, only the vortex effect in the upper layer should be taken into account for gaps ≥ 2.5.(c) At the gap of 3.5D: in the upper layer, the hydrodynamic interaction still exists between two piles (see Figure 17(c)). But the hitting force exerted by vortices decreases continuously with the spacing increasing accordingly.(d) At the gap of 5D: plot of instantaneous vertical-averaged vorticity contours in Figure 17(d) shows that the interaction between two piles tends to vanish when the spacing is large enough.Similarly, the spacing between two cylinders for Pile2 can also be classified into three interference regions distinguished by different hydrodynamic interactions: 1.5
≤
L
/
D
<
2.5, 2.5
≤
L
/
D
<
3.5, and 3.5
≤
L
/
D
≤
5. The force behaviors of Pile2 in these regions will also be further discussed below.


## 4.4. ISWs Forces on Pile2

The three-dimensional simulation results ofC
F
x
m
a
x on Pile2 are shown in Figure 19. It can be seen that all the five configurations exhibit a similar behavior over L
/
D, except that the systematic differences in the values of C
F
x
m
a
x are observed along the entire L
/
D range investigated. By comparing the following five curves of the measured C
F
x
m
a
x of Pile2, the systematic difference between two adjacent curves diminishes gradually with the increase of η
a/H. It indicates that the influence by wave amplitude on horizontal forces tends to depression. Similar phenomenon can also be found in Figure 13 of  C
F
x
m
a
x on Pile1.Figure 19
C
F
x
m
a
x as a function of L
/
D for different η
a/H of Pile2.For gaps in the range of1.5
≤
L
/
D
≤
2.5, as illustrated in Figure 19, C
F
x
m
a
x of each case increases rapidly with the increasing spacing. First, C
F
x
m
a
x meets its minimum at the gap of 1.5D. This fact can be understood by the weakening of the adverse pressure gradient and the influence of lower-pressure area (studied in Section 4.3(a)), which makes f
2 (defined in Figure 14) on Pile2 decrease in both layers at this gap. When the gap reaches 2.5D (see in Figure 17(b)) Pile2’s moving out of the lower-pressure area and continuing to be hit by vortices induce rapid increase in f
2 and cause C
F
x
m
a
x to increase sharply to its peak. Afterwards, the impact force by vortices continually decreases as the spacing is being enlarged to 3.5D. So the curves for all the configurations present a descending trend between the gap of 2.5D and 3.5D.If the gap is further enlarged to 5D, the variation of C
F
x
m
a
x with cylinder centre-to-centre spacing is very small as shown in Figure 19. Plot of instantaneous vertical-averaged vorticity contours (see in Figure 14(d)) shows that interaction between two cylinders at the gap of 5D gradually decreases into a noninteracting state. Therefore, it is clear that, at this gap, the horizontal force on the multiple slender structure is altered in a similar way to that observed in the isolated cylinder case.


## 5. Conclusion

A three-dimensional numerical wave flume is adopted to investigate the ISWs force exerted upon two tandem cylinders. Different centre-to-centre spacings (L) between two piles along with various ISWs amplitudes (η
a) are considered in current study, which can provide a better understanding of the hydrodynamic interference and the ISWs force behaviors of a multiple slender structure. The results for two piles are also compared with that of an isolated cylinder present in the same environment. From study conducted in this paper, the conclusions can be made as follows:(1)
Being different from single-layer flow environment, the horizontal currents induced by internal waves are in the opposite direction between upper and lower layers. So the contributing factors (the effects of vortices and pressure gradient) related to the loading on either pile should be studied separately in each layer at different spacings.
(2)
Hydrodynamic interactions regions are sensitive toL
/
D, and the numerical results obtained from these regions are summarized below:
(a)
For1.5
≤
L
/
D
≤
3.5, strong mutual interference appears between two cylinders. Within this interval, the changes in vorticity fields and pressure gradient fields are considered as the key roles in influencing the force behaviors of two piles. However, the effect of pressure gradient can be ignored when L
/
D ≥ 2 for the upstream pile (Pile1) and L
/
D
≥
2.5 for the downstream pile (Pile2). The dimensionless global horizontal force amplitude (C
F
x
m
a
x) on Pile1 reaches its maximum at 1.5D and then touches minimum at 2D, while C
F
x
m
a
x on Pile2 peaks at 2.5D after experiencing its minimum at 1.5D.
(b)
At the range of3.5
<
L
/
D
<
5, the two cylinders keep on influencing each other in a weak-interference state. C
F
x
m
a
x on both piles illustrates slightly declining trends, because the impact force of vortices exerted upon piles decreases gradually as the gap is further enlarged to 5D.
(c)
Beyond the gap of 5D, the interaction between two piles progressively decreases into a noninteracting state, and C
F
x
m
a
x on the multiple slender structure is presented in a similar way to that on the isolated cylinder.
(3)
An important thing to note is that the inflection points inC
F
x
m
a
x trend curves vary for different piles: L
/
D
=
2 for Pile1 and L
/
D
=
2.5 for Pile2. This is because the near wake between two piles in the upper layer covers a wider range than that in the lower layer, evidenced by comparing Figure 12(b) with Figure 14(b). Pile2 moves out of the low-pressure region at around the gap of 2.5D, while for Pile1 at around 2D.



---

*Source: 9086246-2016-06-27.xml*