Analysis of the Force Characteristics of Two Tandem Cylinders by Internal Waves over Slope Topography

Abstract

Large-amplitude internal waves (IWs) with strong lateral shear forces can cause destructive effects on marine engineering structures. In this study, a large eddy simulation (LES) method was employed to simulate the generation and propagation of IWs in a three-dimensional numerical wave tank, and the pressure distribution, flow field characteristics, and force behaviors of two tandem cylinders under the coupling effect of the IWs and slope terrain were studied. The influence mechanism of the normalized value of center-to-center spacing (L) and the diameter of the cylinder (D), i.e., (L/D), on the strength of the vortex disturbance between cylinders was studied by comparing the simulation results of two tandem cylinders with those of a single cylinder (SC) to further explore the mechanical response characteristics of the upstream cylinder (P₁) and downstream cylinder (P₂). The simulation results showed that the vortex interaction between cylinders is the critical factor that affects the forces acting on the cylinders. The strength of the vortex disturbance could be distinguished by the dimensionless critical center-to-center spacing between cylinders (L_c/D = 3.0). When L/D ≤ L_c/D, the vortex disturbance was severe, causing P₁ and P₂ to experience significant horizontal positive forces and negative forces, respectively. In the case of L/D > L_c/D, the forces acting on both cylinders gradually returned to those on a single cylinder.

1. Introduction

Internal waves (IWs) are the density stratification phenomena caused by vertical density differences in fluid salinity or temperature in various waters, which then cause fluctuations in the density interface (pycnocline) [1]. IWs have the characteristics of large amplitude, long period, and tremendous energy [2]. The IWs generated in stratified strong shear environments pose a serious threat to the supporting structures of risers, drilling platforms, or cross-sea bridges in the complex hydrodynamic environments of the near shore or estuaries [3]. Therefore, IWs are an important consideration in the design of underwater supporting structures.

Many scholars have conducted research on the characteristics of IW forces acting on underwater structures. Cai et al. [4] adopted Morison’s empirical formula to study the characteristics of the forces and moments exerted by IWs on a cylinder in a shear flow environment. Yang Fan et al. [5] established a mathematical model simulating the interaction between IWs and underwater vehicles in a two-layer fluid system, and analyzed the changes in horizontal forces, vertical forces, and moment of the submarine under the action of IWs. Wang Lingling et al. [6,7] employed physical experiments and numerical simulations to compare the force characteristics of cylinders with different shapes in IW environments. Helfrich [8] conducted a series of physical investigations to explore the propagation characteristics of IWs over slope terrains, including the phenomena of climbing and shallow fragmentation. The experimental results indicated that the IWs generated significant vertical mixing, and approximately 10–20% of the wave energy was consumed during vertical mixing. During the propagation of IWs over the coastal slope, continuous changes in the waveform had a continuous impact on the hydrodynamic characteristics, leading to shear instability, turbulent mixing [9], and exacerbating the shear effect of the density thermocline [10]. When the slope terrain causes the thickness of the upper and lower layers to be equal, the polarity of the waveform may even flip [11], transitioning from a depression wave to an elevation wave [12], thereby causing flow and vortex direction reversal. Accordingly, the topography has a significant impact on the propagation and evolution of IWs, greatly altering the hydrodynamic characteristics.

Moreover, Zdravkovich [13] noted that when underwater cylinders are arranged at a certain center-to-center distance, the resultant forces on multiple cylinders, as well as the surrounding flow field, are significantly different from those of a single cylinder (SC) at the same Reynolds number. Gopalan [14] numerically studied the flow field characteristics around two tandem cylinders and found that the center-to-center spacing between the cylinders has a significant influence on the flow field. The forces on the cylinders undergo sudden changes at a critical spacing. When the spacing (L) is less than the critical spacing (L_c) [15], the vortex disturbance between the cylinders will greatly affect the flow field evolution, leading to unpredictable changes in the flow field, pressure distribution, and force behaviors [16]. Wang [17] numerically studied the influence of the spacing between two tandem cylinders on the intensity of the vortex disturbance, exploring the mechanical properties of P₁ and P₂.

In summary, IWs can easily cause instability and damage to underwater engineering supporting structures, even leading to loss of control events. IWs have become a catastrophic environmental factor that must be considered in the design process of engineering supporting structures in complex stratified water environments. However, there have been few investigations on the force characteristics of cylinders in the coupled environment of slope terrain and IWs, and it is rare to find related studies on the forces acting on multicylinder structures in such a complex environment. In this study, a series of numerical simulations were carried out to solve the above engineering problems by establishing a three-dimensional numerical wave tank and adopting LES technology to examine the influence mechanism of IWs on the mechanical response mechanism of two tandem cylinders over the slope topography. The research findings can provide a reference for safety and stability research on cylindrical supporting structures in stratified strong shear slope environments. In Section 1, the research background and significance are presented. Section 2 provides the governing equations used in the numerical simulation. In Section 3, the numerical model is established and validated. In Section 4, the mechanical laws and flow field characteristics of P₁ and P₂ are discussed, and the conclusions are given in Section 5.

2. Governing Equation

(1) Momentum equations and continuity equations

In the description of the incompressible viscous fluid three-dimensional transient motion process, the Navier–Stokes equations (NS) are adopted:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0 \left(i=1,2,3\right)$$

(1)

$$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_i}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial u_i}{\partial x_j}\right)+f_i \left(i,j=1,2,3\right)$$

(2)

where ρ is the density of the fluid; $x_i$ (i = 1, 2, 3) are the three direction coordinates of the Cartesian coordinate system; and t is the time. Correspondingly, p is the pressure; u is the velocity term; f is the body force; and µ is the kinematic viscosity.

(2) Scalar transport equation

Scalar transport equations are used to describe the transport processes of water density of different layers in this study:

$$\frac{\partial C}{\partial t}+\frac{\partial \left(u_{i}C\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(k\frac{\partial C}{\partial x_j}\right)+S$$

(3)

$$\rho =C+{\rho }_0$$

(4)

where C is the scalar volume concentration; k is the molecular diffusion coefficient; S is the source term or the sink term; and ${\rho }_0$ is the density of clean water.

(3) Turbulence model

The large eddy simulation (LES) is a spatial low-pass filtering technique employed to capture turbulent fluctuations. It separates large-scale vortexes from small-scale vortexes by a filter function. The large-scale vortexes are directly simulated, and relatively the small-scale vortexes are enclosed by models.

$$\frac{\partial \overline{u_j}}{\partial x_{ij}}=0$$

(5)

$$\frac{\partial \overline{u_i}}{\partial t}+\overline{u_j}\frac{\partial \left(\overline{u_i}\right)}{\partial x_{ij}}=-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(v\frac{\partial \overline{u_i}}{\partial x_j}\right)+\frac{\partial {\tau }_{ij}}{\partial x_i}+\overline{f_i}$$

(6)

$$\frac{\partial \overline{C}}{\partial t}+\overline{u_{ji}}\frac{\partial \overline{C}}{\partial x_j}=k\frac{{\partial }^2\overline{C}}{\partial x_j\partial x_j}+\frac{\partial {\chi }_j}{\partial x_j}$$

(7)

where the overbar represents the filtering function; the subgrid stress tensor ${\tau }_{ij}=\overline{u_i} \overline{u_j}-\overline{u_i{u}_j}$; the subgrid scale flux ${\chi }_j=\overline{C}\overline{u_j}-\overline{C{u}_j}$. As a function that includes the filtering scale and deformation rate tensor, the expression of turbulent eddy viscosity $v_t$ at time t can be described as:

$${\nu }_t={\left(C_s\Delta \right)}^2{\left(2\overline{S_{ij}}\overline{S_{ij}}\right)}^{\frac{1}{2}}$$

(8)

where C_s is the Smagorinsky constant and Δ is the filtered scale. In the IW environment, C_s changes with time and space. In this paper, Germano’s dynamic procedure method is adopted to determine $C_s$ [18]. The deformation rate tensor $\overline{S_{ij}}$ can be expressed as:

$$\overline{S_{ij}}=\frac{1}{2}\left(\frac{{\partial }^2\overline{u_i}}{\partial x_j}+\frac{{\partial }^2\overline{u_j}}{\partial x_j}\right)$$

(9)

3. Numerical Models

In this study, a three-dimensional numerical tank with slope terrain and two tandem cylinders is used to simulate the generation and evolution of depression IWs and the force mechanism of the cylinders (upstream cylinder P₁ and downstream cylinder P₂) over the slope using LES technology. The three-dimensional dimensions of the numerical model are length (X) × width (Z) × height (Y) = 4 m × 0.3 m × 0.3 m, respectively. The dimensions of the slope terrain are length (X) × width (Z) × height (Y) = 2.15 m × 0.3 m × 0.1 m, respectively. P₁ is located at the width center of the slope platform. The coordinates of the bottom center point are located at: (x, y, z) = (2.0, 0.15, 0.15), and the diameter of the cylinder D = 5 cm. The dimensionless center-to-center spacing between P₁ and P₂ is defined as L/D. Figure 1 illustrates the sketch map of the three-dimensional numerical tank.

The formula for calculating the dimensionless horizontal resultant force acting on a cylinder in an IW environment can be defined as:

$$C_{Fn}=\frac{F_n}{\rho gAH}$$

(10)

where $F_n$ is the horizontal resultant force on the cylinder in this numerical simulation; A is cylinder windward (frontal) side; g is the gravitational acceleration; and H is the total depth in the water tank. The stratified water parameters in the numerical tank are set as follows: the total water depth of the tank is H = 30 cm; the upper layer water density is 0.998 g/cm³, and the layer thickness $h_1=7.5 \mathrm{c}\mathrm{m}$; the lower layer water density is 1.020 g/cm³, and the layer thickness $h_2=22.5 \mathrm{c}\mathrm{m}$. The depression IWs are produced by the gravity collapse method [19] in this simulation. The numerical tank is divided into two parts along the X direction: the wave generation area and the propagation area. The height difference $∆h$ between the generation area and the propagation area is the step height, and the length of the generation area $L_0$ is the step length. The slope angle in the model is $\alpha$ = 45°. The specific parameters of the numerical tank are shown in Figure 2.

3.1. Numerical Methods and Boundary Conditions

The large eddy simulation (LES) technique is used to simulate the generation and evolution of the IWs of a depression type in a three-dimensional numerical tank. The LES model of Smagorinsky type is a versatile turbulence model that can be applied to various flow studies, including wall-bounded turbulence as well as free shear flows like jets [20,21] and stratified flows of internal waves [22,23]. The finite-volume method is employed to discretize the control equations. The pressure implicit with splitting of operators (PISO) algorithm is employed to couple velocity and pressure, ensuring mass conservation and thus obtaining the pressure field. The diffusion and convection terms are discretized using the second-order central difference scheme and the second-order upwind scheme, respectively, and the time step is discretized by the second-order implicit scheme.

The left boundary of the wave generation region, the bottom, the tank sidewalls, and the surface of the cylinders are all set as nonslip solid-rigid wall boundaries. To avoid IWs from reflecting, the right end of the tank is defined as a Sommerfeld radiation boundary [24]. The tank top adopts the “rigid-lid” approximation to ignore the influence of surface waves [25].

3.2. Grid Independent Analysis

The corresponding Reynolds number of the internal waves environment ${Re}_{\mathrm{w}}=\frac{{\alpha }_w{c}_w}{\nu } \approx 4500$ [23] is calculated in this paper, where ${\alpha }_w$ is the wave amplitude, $c_w$ is the wave phase speed, and $\nu$ is the kinematic viscosity. $c_w$ can be expressed as:

$$c_{\mathrm{w}}=c_{\mathrm{o}}\left[1-\frac{{\alpha }_w\left(h_1-h_2\right)}{2h_1{h}_2} \right]$$

(11)

$$c_{\mathrm{o}}=\sqrt{\left(g∆\rho h_1{h}_2\right)/{\rho }_2\left(h_1{+h}_2\right)}$$

(12)

where c_o is the linear velocity of IWs.

Grid independence analysis is conducted based on an SC in the environment without terrain, when η_o/H = 0.0575, as shown in Figure 3. Detail views of the grid around the cylinder with different grid densities are illustrated in Figure 4a–c. Figure 5 shows the calculation results of the IW horizontal resultant forces $C_{Fn}$ corresponding to the three grid density cases. When the grid density increases from a low density (L) to a moderate density (M), $C_{Fn}$ on the cylinder clearly changes, with a difference of 7.2% in $C_{Fn-max}$ between the two cases. The difference in $C_{Fn}$ between a moderate density (M) and a high density (H) is very slight, with only a difference of 0.5%. Considering that computing resources can be saved and the IWs’ propagation process and the force behaviors of the cylinder can be accurately simulated simultaneously, it is feasible to choose a moderate density grid to discretize the computational domain. The grid independence analysis details are shown in

Table 1. Details of the grid independence analysis.

Case	∆t (s)	CFn−max	Elements Number
L (low density)	0.02	0.0906	525,454
M (moderate density)	0.01	0.0976	2,384,640
H (high density)	0.006	0.0981	3,318,278

.

3.3. Numerical Simulation Results Verified by the Physical Experimental Results

The physical model experiment conducted by Wang Fei [26] is used to validate the numerical model in this study. The dimensions of the tank are length (X) × width (Z) × height (Y) = 16.0 m × 0.35 m × 0.57 m, respectively. The upper layer of clear water has a thickness of h₁ = 0.95 m and a density of ${\rho }_1$ = 1000 kg/m³, while the lower layer has a thickness of h₂ = 0.475 m and a density of ${\rho }_2$ = 1024 kg/m³. By opening the gate, IWs are generated and they propagate towards the right end under the influence of gravity collapse [19,27], as shown in Figure 6. A vertical cylinder is placed at a distance of 10 m from the left end of the tank. The cylindrical body is divided into six segments, each with a length of 8 cm, and numbered from top to bottom as cylinders 1–6.

The dimensions and the water stratification parameters of the numerical tank are set according to the physical model experiment conducted by Wang Fei [26]. Following the modeling and meshing methods described earlier, a three-dimensional numerical model was established. The model verification results are shown in Figure 7, where the numerical calculations are in good agreement with the experimental measurements, except for the missing experimental data for case cylinder 3. This indicates that the established numerical tank in this study is reliable and capable of accurately simulating the propagation of IWs and the forces on the cylinders.

4. Results and Analysis

4.1. Effects of Slope Terrain on the Forces Acting on SC and Two Tandem Cylinders

Figure 8a–e presents the duration curves of the dimensionless horizontal resultant forces $C_{Fn}$ acting on the single cylinder, upstream cylinder P₁, and downstream cylinder P₂ for different cases, including the no-slope SC model (Case N₁), the slope terrain SC model (Case N₂), the no-slope double cylinder model (Cases N₃ and N₅), and the slope terrain double cylinder model (Cases N₄ and N₆), with a dimensionless wave amplitude of η_o/H = 0.057. The case details are provided in

Table 2. Detailed parameters in cases.

Case	h1/h2	η0/H	L/D	CFn−max (SC)	CFn−max (P1)	CFn−max (P2)
N1	0.33	0.057	/	0.0856	/	/
N2	0.33	0.057	/	−0.0725	/	/
N3	0.33	0.057	1.5	/	0.1084	−0.0575
N4	0.33	0.057	1.5	/	0.1575	−0.1217
N5	0.33	0.057	6.0	/	0.0814	0.0816
N6	0.33	0.057	6.0	/	−0.0738	0.0448

.

From Figure 8a,c,e, when the spacing between the two cylinders is large, the maximum values of the horizontal resultant forces $C_{Fn-max}$ on P₁ and P₂ in the slope terrain cases are negative and slightly smaller in magnitude than those in the no-slope case. When the spacing between the cylinders is small, as shown in Figure 8b,d, there is a strong interaction between P₁ and P₂, resulting in more complex forces on the cylinders. A detailed analysis will be provided in the following sections. To further analyze the force characteristics of the cylinders under different cases, the worst condition of the forces on the cylinders (${\mathrm{when} C}_{Fn}=C_{Fn-max}$) in each case is selected as the research object.

When depression-type IWs propagate, the flow in the upper and lower layers is in opposite directions. In the upper layer, the direction of the flow is the same as the IW propagation direction. In the lower layer, however, the flow direction is opposite to the IW propagation direction. Each cylinder body is divided into 10 segments along the vertical direction, with each segment having a length of 0.02 m, as shown in Figure 9. The portion of the cylinder located in the upper water layer is called the “upper parts”, while the portion in the lower water layer is referred to as the “lower parts”. In the following analysis, the vertical distribution of the dimensionless horizontal force coefficient ($C_f$) along the cylinder sections is calculated. The definition of $C_f$ is based on Equation (10), which is similar to that of $C_{Fn}$.

When the spacing between the cylinders is large, as observed from Figure 10a,c,e, the slope terrain leads to a decrease in the positive horizontal force on the cylinder (in the same direction to that during IW propagation) in the upper layer and an increase in the negative horizontal force on the cylinder (in the opposite direction to that during IW propagation) in the lower layer. This results in a negative $C_{Fn-max}$ on the cylinders in the presence of slope terrain. When the spacing between the cylinders is small, see Figure 10b,d, the vertical distribution of $C_f$ for the cylinders can be either positive or negative. However, the presence of slope terrain increases the magnitudes of $C_{Fn-max}$ for both P₁ and P₂ (as shown in Figure 8b,d).

Further analysis of the differences in forces under different cases is conducted by comparing the pressure distribution along the circumference of the cylinder at two specific sections: Y = 0.16 m (the mid-depth of the lower layer) and Y = 0.26 m (the mid-depth of the upper layer). The schematic diagram of the cylinder circumferential angles is shown in Figure 11.

Figure 12a–e presents the circumferential pressure distribution for different cases at Y = 0.16 m, where C_p is the pressure on the periphery of the cylinder cross section. In the lower layer, the slope terrain has a significant impact on the pressure on the lee side of P₁ and P₂, while the pressure differences on the frontal side are minimal. In the upper layer, as shown in Figure 13, the presence of slope terrain leads to significant changes in the pressure distribution on the frontal side of P₁ and P₂, with minimal differences in the pressure distribution on the lee side. These results indicate that the slope terrain has a substantial influence on the pressure and forces experienced by single and two tandem cylinders. In the following sections, the mechanisms of the combined effects of slope topography and IWs on the flow field, pressure distribution, and forces acting on P₁ and P₂ will be specifically analyzed, considering different L/D cases.

4.2. Effects of L/D on the Force Behaviors and Flow Field Characteristics for P₁

Figure 14 presents the duration curves of $C_{Fn}$ acting on P₁ for eight different L/D cases (1.5 ≤ L/D ≤ 6, cases C₂–C₉) compared with an SC case (case C₁), when η_o/H = 0.057. As shown in the figure, $C_{Fn}$ acting on the cylinder changes with time, and the overall trend of the force curves remains consistent across the nine cases. When L/D < 3.0, the maximum $C_{Fn-max}$ value for P₁ is observed at L/D = 1.5. For 3.0 ≤ L/D ≤ 6, the force duration curves of P₁ exhibit a high degree of similarity, with less variation in the $C_{Fn-max}$ value as L/D changes. These results indicate that L/D = 3.0 can be defined as the dimensionless critical spacing L_c/D. When L/D < L_c/D, there is strong mutual interference between the cylinders. Alternatively, when L/D ≥ L_c/D, the mutual interference gradually diminishes. Moreover, when the spacing between the cylinders is sufficiently large, the force acting on P₁ approaches that on the SC.

Next, we will compare the pressure distribution and flow field of P₁ and SC for different L/D ratios to reveal the mechanical mechanism of P₁ in an IW environment. The specific cases are detailed in

Table 3. Simulation cases for studying the forces acting on P₁.

Case	h1/h2	η0/H	L/D	+CFn–max (P1)
C1 (SC)	0.33	0.057	/	0.0512
C2 (L/D = 1.5)	0.33	0.057	1.5	0.1575
C3 (L/D = 2.0)	0.33	0.057	2.0	0.1556
C4 (L/D = 2.5)	0.33	0.057	2.5	0.0929
C5 (L/D =3.0)	0.33	0.057	3.0	0.0511
C6 (L/D = 3.5)	0.33	0.057	3.5	0.0482
C7 (L/D = 4.0)	0.33	0.057	4.0	0.0470
C8 (L/D = 5.0)	0.33	0.057	5.0	0.0440
C9 (L/D = 6.0)	0.33	0.057	6.0	0.0494

.

4.2.1. Distinction of the Flow Field and Pressure Distribution between P₁ and SC (L/D = 1.5)

Figure 15 shows the vertical distribution of $C_f$ for P₁ for the cases of SC and L/D = 1.5. P₁ and the SC exhibit almost identical $C_f$ distributions in the upper layer. However, in the lower layer, P₁ experiences significantly higher $C_f$ values compared to the SC, resulting in a larger postive $C_{Fn}$ on P₁. Therefore, when L/D = 1.5, the focus is primarily on examining the force characteristics of the lower parts of the cylinder. The following analysis will further investigate the force mechanisms acting on P₁ by comparing the flow field evolution and circumferential pressure distribution when Y = 0.16 m.

Based on the vorticity contour maps shown in Figure 16a,b, it can be observed that when L/D = 1.5, there is an intense disturbance between the two cylinders, and the development of vortices behind P₁ is suppressed. The presence of vortices caused by P₂ clearly contributes to altering the flow field distribution around P₁ in the downstream vortex region.

From the pressure distribution diagram shown in Figure 17, the presence of P₂ causes the lee side of P₁ to be submerged in an area of low pressure. The vortices between the two cylinders are suppressed, significantly reducing the pressure on the lee side of P₁. On the other hand, the pressure on the frontal side of P₁ changes slightly. As a result, $C_f$ acting on P₁ in the positive direction increases. This explains why $C_f$ acting on P₁ is larger than that on the SC (as shown in Figure 15) in the lower layer.

4.2.2. Distinction of the Flow Field and Pressure Distribution between P₁ and SC (2.0 ≤ L/D ≤ 3.0)

As shown in Figure 14, the maximum $C_{Fn-max}$ on P₁ occurs at L/D = 1.5, and it gradually decreases with increasing L/D. When L/D ≥ 3.0, the forces on P₁ approach those on the SC. Therefore, a comparison between P₁ at L/D = 2, L/D = 2.5, L/D = 3.0, and SC is conducted to explore the corresponding mechanical characteristics and flow field.

The comparison of the vertical distribution of $C_f$ on P₁ in the SC, L/D = 2, L/D = 2.5, and L/D = 3.0 is illustrated in Figure 18. For L/D = 2.0, L/D = 2.5, and L/D = 3.0, the upper parts of the cylinder exhibit a slightly smaller $C_f$ than the SC. However, $C_f$ on the lower parts of the cylinder are significantly greater than those on the SC. As a result, the horizontal resultant force $C_{Fn}$ acting on P₁ is greater than that on the SC.

Similarly, the flow field and pressure distribution at Y = 0.16 m are illustrated in Figure 19 and Figure 20, respectively. Compared to the SC, with increasing L/D from 2.0 to 3.0, the suppression of vortices between the two cylinders gradually diminishes (see Figure 19), and the pressure distribution around the circumference of P₁ gradually returns to a state similar to that of the SC (see Figure 20). Therefore, when 2.0 ≤ L/D ≤ 3.0, the differences in forces on the cylinder are also evident in the lower layer, mainly concentrated in the 0°–30° and 330°–360° sections of the cylinder circumference.

4.2.3. Distinction of the Flow Field and Pressure Distribution between P₁ and SC (3.0 < L/D ≤ 6.0)

As shown in Figure 21, when L/D > 3.0, the $C_f$ values on the upper and lower parts of P₁ are both similar to those on the SC. Therefore, if L/D is larger than L_c/D, then the mutual interference between the two cylinders gradually transitions to an undisturbed state, and the flow field and circumferential pressure distribution around the two cylinders gradually return to a state similar to that of SC. Consequently, further discussion on the flow field and pressure distribution characteristics is not necessary.

In conclusion, the size and location of the vortex region are determined by L/D, evidently affecting the forces acting on P₁. Figure 15, Figure 18 and Figure 21 show that the presence of P₂ has a significant impact on the forces acting on P₁ when L/D < L_c/D, with the force differences mainly manifesting in the lower layer.

4.3. Effects of L/D on the Force Behaviors and Flow Field Characteristics for P₂

Similarly, when η_o/H = 0.057, the duration curves of $C_{Fn}$ for P₂ under different L/D cases (1.5 ≤ L/D ≤ 6) are shown in Figure 22. $C_{Fn}$ on P₂ varies significantly over time. When L/D ≥ 3, the overall trend of the $C_{Fn}$ curves for P₂ is similar to that of the SC, although the $C_{Fn-max}$ on the SC is slightly larger than that of L/D ≥ 3. However, when 1.5 ≤ L/D ≤ 2, P₂ experiences a larger negative $C_{Fn}$ (opposite to the wave propagation direction), and the $C_{Fn-max}$ on P₂ is negative. By employing the same analysis method as that for P₁, the comparison of the flow field and pressure distribution between the SC and P₂ at different L/D cases is shown to reveal the force characteristics of P₂ in an IW environment. The specific cases are detailed in

Table 4. Simulation cases for studying the forces acting on P₂.

Case	h1/h2	η0/H	L/D	−CFn−max (P2)
C1 (SC)	0.33	0.057	/	−0.0725
C2 (L/D = 1.5)	0.33	0.057	1.5	−0.1217
C3 (L/D = 2.0)	0.33	0.057	2.0	−0.1124
C4 (L/D = 2.5)	0.33	0.057	2.5	−0.0433
C5 (L/D = 3.0)	0.33	0.057	3.0	−0.0459
C6 (L/D = 3.5)	0.33	0.057	3.5	−0.0458
C7 (L/D = 4.0)	0.33	0.057	4.0	−0.0455
C8 (L/D = 5.0)	0.33	0.057	5.0	−0.0449
C9 (L/D = 6.0)	0.33	0.057	6.0	−0.0445

.

Figure 23 illustrates the vertical distribution of $C_f$ for P₂ at different L/D cases. When 1.5 ≤ L/D < 3.0, the upper parts of P₂ experience significantly smaller forces compared to the SC. However, in the lower layer, the situation is reversed, and the lower parts of P₂ suffer significantly higher forces than those of the SC. When 3 ≤ L/D ≤ 6.0, the forces on the upper parts of P₂ are similar to those on the SC, but the forces on the lower parts are significantly higher. Therefore, the investigation of the forces on P₂ primarily focuses on the flow field and pressure distribution characteristics in the upper layer for 1.5 ≤ L/D < 3.0 and in the lower layer for 1.5 ≤ L/D ≤ 6.

4.3.1. Distinction of the Flow Field and Pressure Distribution in the Upper Layer between P₂ and SC (1.5 ≤ L/D < 3.0)

Figure 24a–c presents the flow field contour maps. It can be observed that when 1.5 ≤ L/D < 3.0, the frontal side of P₂ is immersed in a low-pressure vortex area, significantly reducing the pressure on the frontal side of P₂ (see Figure 25). Clearly, the pressure on the frontal side of P₂ is significantly lower than that on the SC. Consequently, the pressure on the frontal side of P₂ is lower than that on the lee side, compelling the forces on P₂ to be negative.

4.3.2. Distinction of the Flow Field and Pressure Distribution in the Lower Layer between P₂ and SC (1.5 ≤ L/D ≤ 6)

The comparison of the flow field in Figure 26a,b shows that when L/D = 1.5, the frontal side of P₂ is still immersed in vortices, resulting in lower pressure on the frontal side than on the SC (as shown in Figure 27). When L/D = 6.0, the development of vortices in front of P₂ still differs significantly from the SC (as observed by comparing Figure 26a,c). As a result, the pressure on the frontal side of P₂ remains lower than that on the SC (as shown in Figure 27). This explains why, when L/D ≥ 3, the overall force trend of P₂ is similar to that of the SC, but $C_{Fn-max}$ on P₂ is slightly smaller than that on the SC (as shown in Figure 22).

To summarize, the size and position of the vortex region are significantly determined by L/D, directly affecting the forces acting on P₂. Combining the observations from Figure 24, Figure 25, Figure 26 and Figure 27, it can be concluded that for 1.5 ≤ L/D < 3.0, the differences in forces between P₂ and the SC are primarily determined by the pressure differences on the frontal side in the upper layer and the lee side in the lower layer. For 3 ≤ L/D ≤ 6, the differences in forces between P₂ and the SC are primarily determined by the pressure differences on the lee side in the lower layer.

5. Discussion

(1) It should be noted that the LES (large eddy simulation) method used in this study has some level of uncertainties in terms of accuracy compared to the direct numerical simulation (DNS) method, particularly in the context of low to moderate Reynolds numbers. In future research, the DNS method should be employed specifically for cases involving these Reynolds numbers.

(2) We carried out the grid independent analysis, with grids having (low density) 525,454, (moderate density) 2,384,640, and (high density) 3,318,278 grid points. In future research, more thorough grid refinement studies will be conducted, such as adjusting the grid number with medium density to eight times that of low density, and further increasing the grid number difference between high grid density and medium grid density.

(3) The corresponding Reynolds number of the IWs environment ${Re}_{\mathrm{w}}\approx 4500$ is calculated in this paper. The Reynolds number is typically quite small, usually on the order of 10³, which is consistent with the Reynolds numbers found in the relevant studies by Terletska [28], Harnanan [29], Hsieh [30], Arthur [31], and Sutherland [32], all of which are at laboratory scales. Due to limitations in computational resources, a grid quantity of 3,318,278 (chosen for the high-density case of the grid independent analysis) is essentially the computational limit for our current study, requiring approximately 130 h of computation on a quad-core CPU.

In our future research, we will also make every effort to secure funding for the purchase of more powerful computers to further refine our model and grid.

6. Conclusions

This study investigates the force characteristics of two tandem cylinders in the presence of an IW environment over slope terrain using a three-dimensional numerical wave tank. By analyzing the force duration curves, vorticity contour maps, flow field distributions, and pressure distribution maps, the study provides a detailed understanding of the force responses of P₁ and P₂ for different L/D cases. The conclusions are described as follows:

(1) When IWs interact with the slope terrain, the flow field around the cylinders becomes more complex and dynamic, leading to a significant alteration of the force characteristics. The horizontal force acting on the cylinder over the slope is opposite to that action on the cylinder without the slope, with the interaction between IWs and topography playing a crucial role. Bounded by the pycnocline, the impact of the topography on the forces is mainly observed on the lee side of the upper parts and the frontal side of the lower parts of the cylinders.

(2) The dimensionless spacing between the cylinders (L/D) has a significant influence on the forces on P₁ and P₂. A critical spacing, defined as L_c/D = 3.0, is identified to differentiate between strong and weak interference regimes. In the strong interference regime (L/D ≤ L_c/D), P₁ and P₂ experience larger forces in the downstream and upstream directions, respectively. In the weak interference regime (L/D > L_c/D), the force responses of P₁ and P₂ gradually close to those of an SC, but with slightly lower maximum forces. The effects of spacing are predominantly observed in the lower parts of the cylinders.

(3) The pressure distribution and flow field characteristics around the cylinders reveal the mutual disturbance mechanism between P₁ and P₂ under the coupling effect of IWs and slope topography. L/D determines the size and location of the vortex region, directly influencing the force characteristics of the cylinders. The vortex disturbance in the lower layer is more pronounced. When L/D ≤ L_c/D, the influence of the vortexes in the upper layer on the force on the cylinder can be ignored. The presence of P₂ significantly affects the forces on P₁, with a difference from the forces on the SC primarily observed in the flow field and pressure distribution in the lower layer. Similarly, the presence of P₁ also significantly influences the forces on P₂, with differences from the forces on the SC manifested in both the upper and lower layers in terms of the flow field and pressure distribution.