Fluid Dynamics & Materials Processing

DOI: 10.32604/fdmp.2021.014652

ARTICLE

A Numerical Study on the Mechanisms Producing Forces on Cylinders Interacting with Stratified Shear Environments

Yin Wang1,*, Lingling Wang2, Yong Ji1, Zhicheng Xi1 and Wenwen Zhang1

1College of Water Conservancy and Ecological Engineering, Nanchang Institute of Technology, Nanchang, 330099, China
2College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing, 210098, China
*Corresponding Author: Yin Wang. Email: 2008984059@nit.edu.cn
Received: 16 October 2020; Accepted: 05 February 2021

1  Introduction

Internal waves (IWs) usually occur in fluids with stable density stratification. Stable stratified flow exists in oceans, estuaries, and lakes, while fluid density varies along the vertical direction induced by temperature and salinity differences and other environmental conditions [1]. In such an environment, IWs carrying enormous energy can be produced by a small perturbation at the density interface (pycnocline) [2]. Internal solitary waves (ISWs) are a particular type of nonlinear IW with a large amplitude, short period, and considerable energy [3]. Large amplitude ISWs significantly impact the safety and stability of underwater structures in coasts and estuaries with complex hydrodynamic environments [4].

During the propagation of ISWs, fluid in the upper and lower layers, bounded by the pycnocline, flows in opposite directions. The enormous energy carried by the strong shear flow will lead to an unusually high velocity. The maximum wave-induced force on a cylinder exerted by an ISW with a velocity of 2.1 m/s is equivalent to the force exerted by a surface wave with a wavelength of 300 m and a wave height of 18 m [5].This strong force could be a severe threat to underwater structures, such as oil drilling platforms or supporting cylinders [6]. It is of great significance to carry out in-depth research on the mechanism of ISW dynamics and ISW-induced forces on underwater structures.

The horizontal force on cylinders can be divided into differential pressure resistance and frictional resistance(viscous force) [7]. The differential pressure resistance is the horizontal force generated by the pressure difference between the windward side (upstream surface) and leeside (downstream surface) of the cylinder, while the viscous force is the horizontal force generated by fluid viscosity. Most previous studies have focused on the overall force on a cylinder exerted by ISWs and very few studies have precisely obtained the force distribution on different parts of the cylinder in the upper layer or in the lower layer [8–10]. The detailed mechanical characteristics of the cylinders are still not fully studied.

In the current research, a three-dimensional (3D) numerical wave flume is used to investigate the ISW forces acting on cylinders meticulously through a large-eddy Simulation (LES) approach. By extracting the pressure distribution and tangential velocity of the cylinder based on the numerical simulation results, the infinitesimal method is adopted to obtain the differential pressure resistance and the viscous force on the cylinder in each layer by the method of integration. Finally, the vertical distribution of the IWS forces along depth can be obtained to reveal the force mechanism of the cylinder in the stratified flow environment.

2  Mathematical Model

2.1 Governing Equations

Based on Navier-Stokes (N-S) equations, the moving course of 3D unsteady incompressible fluid can be expressed by:

∂ui∂xi=0 (1)

∂ui∂t+∂(uiuj)∂xj=−1ρ∂p∂xi+ν∂2ui∂xj∂xj+fi (2)

where ρ is the density; t is the time; ui is the velocity component; p is the pressure; ν is the kinetic viscosity and fi is the body force.

2.2 Scalar Transport Equation

In current investigation, the internal solitary wave is generated by the salinity variation in a two-phase flow system. The convective-diffusive effect can be described by the following equation:

∂C∂t+∂(ujC)∂xj=∂∂xj(k∂C∂xj)+S (3)

where C represents a scalar volume concentration ranging from 0 to 1, and the fluid density difference between the two layers can be dominated by the equation: ρ = Cρlow + (1 – C)ρup, (ρlow and ρup are the designed fluid density in the two-layer system, respectively); k is the molecular diffusivity coefficient; and S is the source or sink term.

2.3 Turbulence Model

The turbulence is simulated by Large-eddy simulation (LES) in this paper by using a spatial filter to the two governing equations above [11]. Therefore, the space filtered Eqs. (1)–(3) can be expressed (denoted by an overbar) as:

∂uj¯∂xj=0 (4)

∂ui¯∂t+uj¯∂ui¯∂xj=−1ρ∂p¯∂xi+v∂2ui¯∂xj∂xj+∂τij∂xi+fi¯ (5)

∂C¯∂t+uj¯∂C¯∂xj=k∂2C¯∂xj∂xj+∂χj∂xj (6)

where τij is the sub-grid stress tensor and χj is the Sub-Grid Scale (SGS) flux. The former is used for the momentum exchange between the SGS and the resolved scale (RS), and the latter is applied for the scalar transport between the SGS and the RS. The two parameters can be defined by the formulas below:

τij=ui¯uj¯−uiuj¯ (7)

χj=C¯uj¯−Cuj¯ (8)

Actually, SGS tensor τij can be decomposed into an off-diagonal tensor τtij and a diagonal tensor τtkk yielding these formulas:

τijt=τij−δij3τkkt=−2νtSij¯ (9)

χj=−kt∂C¯∂xj (10)

where νt is the turbulent eddy viscosity; Sij¯ is the strain rate tensor; kt is the eddy diffusion coefficient. νt and Sij¯ can be solved by:

Sij¯=12(∂2ui¯∂xj+∂2uj¯∂xj) (11)

νt=(CsΔ)2(2Sij¯Sij¯)1/2 (12)

where Cs is the Smagorinsky constant and varing in time and space because of the stratified strong shear environment, which can be determined by a dynamic procedure developed by Germano et al. [12]; ∆ is the filtered scale.

2.4 Differential Pressure Resistance on Cylinder

2.4.1 Circular Cylinder

As shown in Fig. 1, P is the pressure on a certain point of the circular cylinder circumference at a depth y, and the horizontal differential pressure resistance on the circular cylinder fc can be defined as:

fc=−∫02πP⋅cos⁡θdl (13)

where θ is the angle between the point pressure P and the X axis, (°), counterclockwise is positive; l is the arc length, cm.

Based on the infinitesimal method, Eq. (13) can be discretized as:

fc=−Δli∑i=1nPi⋅cos⁡θi (14)

where Pi is the pressure acting on the i-th infinitesimal, Pa; θi is the angle between the extension cord of the normal line of the i-th infinitesimal and the X axis, (°), counterclockwise is positive; n is the number of division on the circumference, in this paper, n is taken as 360; △li=2πr /n = 0.087 cm is the arc length of the i-th infinitesimal (r is the cylindrical radius, r is taken as 5 cm in this paper).

Figure 1: Calculation diagram of differential pressure resistance on circular cylinder

2.4.2 Square Cylinder

As shown in Fig. 2, P is the pressure on a certain point of the square cylinder edges at a depth y, and the horizontal differential pressure resistance on the square cylinder fs can be defined as:

fs=fs3−fs1=∑i=190Pi⋅Δsi−∑i=190Pi⋅Δsi (15)

where fsj is the pressure acting on the j-th side of the square cylinder, N , j = 1, 2, 3, 4, as shown in Fig. 2; Pi is the pressure of the i-th infinitesimal on each side of the square cylinder, N; △si=s/4⋅n=0.11 cm is the side length of the i-th infinitesimal (s is the perimeter of the square cylinder , s is taken as 10 cm in this paper; m is the number of division on the edges, m is taken as 360).

Figure 2: Calculation diagram of differential pressure resistance on square cylinder

2.5 Viscous Force on Cylinder

2.5.1 Circular Cylinder

As shown in Fig. 3, the tangential velocity uτ and radial velocity uγ of a certain infinitesimal on the grid nearest to the circumference can be defined as:

uτ=uzcos⁡θ−uxsin⁡θ (16)

uγ=uzsin⁡θ+uxcos⁡θ (17)

where ux and uz are the x-direction and z-direction velocity of the infinitesimal, cm/s; θ is the angle between the infinitesimal a and the X axis, (°), counterclockwise is positive. Then the viscous force Δfcτ along the horizontal direction on a certain infinitesimal i can be expressed as:

Δfcτxi=−∂uτi∂γ⋅μ⋅sin⁡θi⋅Δli (18)

where uτi is the tangential velocity of the i-th infinitesimal, ∂uτi∂γ is the normal gradient of the tangential velocity, cm/s; μ is the dynamic viscosity, Pa⋅s , μ=ν⋅ρ , which ν is the kinematic viscosity.

Figure 3: Calculation diagram of viscous force on circular cylinder

Based on the infinitesimal method, the horizontal viscous force on the circular cylinder fcτ can be expressed as:

fcτ=∑i=1nΔfcτi (19)

2.5.2 Square Cylinder

Following the calculation method of the viscous force on circular cylinder in Section 2.5.1, the viscous force fsτ along the horizontal direction on square cylinder at a depth y can be obtained:

fsτ=fsτ2+fsτ4=∑i=190Δfsτi+∑i=190Δfsτi (20)

where Δfsτi=∂uxi∂γ⋅μ⋅Δsi , fsτj is the viscous force acting on the j-th edge of the square cylinder, N; as shown in Fig. 4, uxi is the velocity along the X direction of the i-th infinitesimal, cm/s; ∂uxi∂γ is the normal gradient of the velocity along the X direction, cm/(s·cm).

Figure 4: Calculation diagram of viscous force on square cylinder

3  Numerical Models

A 3D numerical wave flume is developed in this study, as illustrated in Fig. 5, with a size of 12.0 m in length (X), 0.8 m in height (Y) and 0.6 m in width (Z). Schematic diagrams of the unstructured triangular mesh around the cylinder and the structured mesh of the tank in the vertical direction are shown in Figs. 6 and 7. The numbers of circumference nodes on the cylinder and vertical nodes along the tank are 60 and 85, respectively.

Figure 5: Schematic diagram of the 3D numerical flume, units: meters

Figure 6: Close-up view of the triangular unstructured mesh around a cylinder

Figure 7: Close-up view of the structured mesh of the tank in vertical direction

In the current investigation, the dimensionless horizontal total force F′ on the cylinder can be expressed as [13]:

F′=FρgAH (21)

where F is the total horizontal force on the circular cylinder and square cylinder calculated in the model, N; A is the windward area of the two types of cylinders, cm2; H is the total depth of the numerical flume, cm; and g is the gravitational acceleration, m/s2.

3.1 Cases and Parameters Setting

In current investigations, ISWs of depression type can be produced by a method of gravity collapse in a two-phase flow system [14]. The fluid densities of the upper layer and lower layers are ρ1 = 0.998 g/cm3 and ρ2 = 1.017 g/cm3, and the depths of the upper layer and lower layer are h1 = 0.2 m and h2 = 0.6 m, respectively. Consequently, H is equal to 0.8 m in this study. ηo/H represents the dimensionless wave amplitude when the ISW form is stable after propagating 4 m.

As illustrated in Fig. 8, the whole flume can be divided into a wave generation area and a wave propagation area. The former starts from the left boundary and is 0.7 m in length, and the latter is 11.3 m in length. ho, the height difference of the pycnocline between the two areas, and is called the step height, and the step length Lo, as called in this study is the generation area length.

Figure 8: Schematic diagram of the gravity collapse to produce the ISWs in a two-phase flow system

3.2 Numerical Methods and Boundary Conditions

The generation and propagation of depression-type ISWs are simulated by a large-eddy simulation (LES) model. The governing equations are discretized by the finite-volume method [15]. The SIMPLE algorithm is applied to couple the velocity-pressure term and enforce mass conservation, and then the pressure field can be obtained [16,17]. A bounded central differencing scheme is adopted for solving momentum equation. The diffusion term is solved by the second-order central difference scheme, while the advection term is solved by the second order upwind scheme. The second-order implicit scheme is adopted to discretize the time step.

The left boundary, the two sidewalls, and the bottom of the wave flume are set to be rigid walls with no-slip conditions as illustrated in Fig. 9. To avoid wave reflection, the right boundary is specified by a “Sommerfeld radiation” type [18]. The top boundary is specified by the “rigid lid” approximation [19], because the surface waves are much smaller than the internal solitary waves [20,21]. Therefore, the “rigid lid” approximation is reasonable to use in the current research [22].

Figure 9: Instantaneous three-dimensional vorticity contours around cylinders, when F′ = F′max, ηo/H = 0.121 (a) Circular cylinder (b) Square cylinder

3.3 Numerical Model Verification

The validation of the established numerical model is executed by consulting the same procedure performed in reference [23]. The size of the tank, the meshing method and the boundary conditions used previously in reference [23] are chosen for the present numerical simulations.

4  Simulation Results and Analysis

The velocity field and pressure field in the ISW environment are different from those in the pure current environment, because the three-dimensional characteristics of the flow field induced by ISWs around a cylinder are more obvious as illustrated in Fig. 9. Instantaneous three-dimensional vorticity contours around cylinders clearly show that the vortex sheds in opposite directions in a stratified flow environment, and the vortex shedding from the square cylinder is more intense and erratic. Also, the flow structure and turbulence effect around the square cylinder has a stronger strength three-dimensional character than that around the circular cylinder.

The comparison of the instantaneous streamlines around the two types of cylinders at depths y = 0.3 m(Section P0.3) and y = 0.7 m (Section P0.7) is presented in Figs. 10a and 10b. The vortex around the square cylinder is more asymmetric and is larger in the recirculation zone, which could induce the more uneven pressure distribution on the cylinder.

Figure 10: Instantaneous streamlines around cylinders on section P0.3 and P0.7, when F′ = F′max, ηo/H = 0.121 (a) Circular cylinder (b) Square cylinder

As a result, the ISW-induced force on the cylinder differs significantly from the force induced in a one-phase flow environment. Hence, it is of great significance to carry out research on the hydrodynamic characteristics around cylinders and the pressure distributions on them, which promotes the force behaviors for cylinders to be obtained in strong shear stratified environment.

4.1 Total Dimensionless Horizontal Force on a Cylinder

Based on the numerical simulation results, the trend curves of F′ with time t under the cases of ηo/H of 0.034, 0.075, 0.103 and 0.121 are illustrated in Fig. 11. F′ of each case increases gradually with increasing spacing ηo/H. The ISW forces on the square cylinder are larger than those on the circular cylinder because of the more uneven pressure distribution andr the higher differential pressure between the front and rear of the square cylinder.

Figure 11: The trend curves of F′ with time t for different wave amplitudes

Meanwhile, it can be found that the propagation of ISWs is a transient motion. When passing through the cylinder, the ISW force on the cylinder is continuously changing. Therefore, the most unfavorable moment can be selected as the calculation feature time when F′ reaches Fmax′ . Fmax′ is the ISW-induced dimensionless maximum horizontal force as shown in Fig. 11.

4.2 Differential Pressure Resistance and Viscous Force Distribution on the Cylinder along the Depth

The horizontal force on the cylinder can be divided into differential pressure resistance and viscous force.The basic ideas for solving these two kinds of forces are as follows: ① intercept a cross section at regular intervals on the cylinder along the vertical (Y) direction (9 sections in total), selecting 360 points along the circumference of each cross section. By extracting the pressure and velocity values at these points, the pressure and velocity distribution on circumference can be obtained; ② the infinitesimal method is adopted. By integrating the pressure and velocity distribution on circumferences at various depths, the differential pressure resistance f (Eqs. (14)–(15)) and viscous force fτ (Eqs. (19)–(20)) acting on the circular cylinder and square cylinder can be obtained respectively.

4.2.1 Differential Pressure Resistance

The dimensionless horizontal differential pressure resistance f′ on the cylinder at depth y is defined as:

f′=fρgAH (22)

f′ on circular cylinder fc′ and square cylinder fs′ at various depths y are shown in Tab. 1 , taking the case of ηo/H = 0.098 as an example, and their vertical distribution can be subsequently obtained. With an ISW of the same amplitude, bounded by the pycnocline, the differential pressure resistance fs′ on the square cylinder in both the upper and lower layers are greater than fc′ on the circular cylinder: in the upper layer fs′ is approximately 1.5 times that of fc′ , while in the lower layer fs′ is approximately 3.5 times that of fc′ .

Table 1: the dimensionless horizontal forces on circular and square cylinders, ηo/H = 0.121

4.2.2 Viscous Force

The dimensionless horizontal viscous force fτ′ on the cylinder at depth y is defined as:

fτ′=fτρgAH (23)

fτ′ on circular cylinder fcτ′ and square cylinder fsτ′ at various depths y are shown in Tab. 1, taking the case of ηo/H = 0.098 as an example, and their vertical distribution can be obtained. Comparing f′ with fτ′ , it can be found that the viscous forces on the two types of cylinders are 1–2 orders of magnitude smaller than the differential pressure resistance.

The results show that the viscous force has little contribution to the total horizontal force on cylinders; therefore, the fluid viscosity effect can be ignored and only the effect of differential pressure resistance should be considered in the current study. This is similar to the result obtained under the condition of a pure flow environment: differential pressure resistance accounted for more than 95% of the total horizontal force [7], indicating that the viscosity effect is weak in both fluid environments.

4.3 Quantification Analysis of ISW Force on the Cylinder

Numerous studies have shown that forces on cylinders increase with increasing ISW amplitude [24]. In addition to the amplitude factor, the thicknesses of the upper layer h1 and lower layer h2 also determine the propagation characteristics of ISWs [25]. It can be inferred that the water depth ratio h1/h2 can also be a sensitivity parameter of the forces on cylinders. In this study, the influence rule of the relative amplitude ηo/H combined with h1/h2 on ISW forces is comprehensively considered, and the correlation between the maximum dimensionless horizontal force Fmax′ and ηo/H and h1/h2 in the stratified flow environment is quantitatively analyzed. The calculated Fmax′ values on the circular cylinder under all cases are shown in Tab. 2, keeping the density difference Δρ = 0.019 g/cm3 constant, h1/h2 ranging from 0.231 to 0.778, and ηo/H ranging from 0.034 to 0.158.

Table 2: Cases and parameters setting of ISW force quantification analysis

By analyzing the data in Tab. 2, Fmax′ has a strong correlation with ηo/H and h1/h2. Multiple linear regression is used to perform correlation analysis with h1/h2and ηo/H as independent variables and Fmax′ as the dependent variable: Fmax′ shows a positive correlation with ηo/H and a negative correlation with h1/h2 The optimal linear fitting relation among them can be expressed as:

Fmax′=−0.88h1/h2+3.72η0/H+0.34 (24)

The goodness of fit R2 is only 0.3. Taking the calculated Fmax′ as the x-axis and the true value Fmax′ as the y-axis, an exponential behavior is exhibited. Therefore a comprehensive influence coefficient m is defined as in current study:

m=−0.88h1/h2+3.72η0/H+0.34 (25)

By further logarithmically fitting the true value Fmax′ with m, the following equation can be obtained:

Fmax′=0.05e4.18m (26)

R2 is fitted to be 0.957 with a high fitting degree at this time (as shown in Fig. 12). The comprehensive influence coefficient m has its own physical significance: The larger ηo/H and Δh = h1-h2 are, the larger the value of m and the forces on the cylinder are. m provides a quantitative index for the influence of ηo/H and h1/h2 on forces acting on the cylinder. Moreover, m is positively correlated with ηo/H and negatively correlated with h1/h2 (Eq. (25)). In summary, the forces on the cylinder increase with ηo/H under various stratification conditions; with the ISW of the same amplitude, the greater h1/h2 is, the smaller forces on the cylinder are.

Figure 12: The correlation between F′max and m

5  Conclusion

A 3D numerical wave flume is used to meticulously study the ISW forces acting on cylinders through a large-eddy Simulation (LES) approach. By analyzing the hydrodynamic characteristics around cylinders and the pressure distributions on them, the detailed behaviors of differential pressure resistance and viscous force on cylinders can be obtained. The main conclusions from this study are as follows:

a) with the ISW of the same amplitude, the pressure difference resistance on the square cylinder in both the upper and lower layers is greater than that on the circular cylinder: the force on the square cylinder in the upper layer is approximately 1.5 times that on the circular cylinder, while the force on the square cylinder in the lower approximately is about 3.5 times that on the circular cylinder;

b) the viscous force acting on cylinders is 1-2 orders of magnitude smaller than the differential pressure resistance; therefore, the effect of fluid viscosity can be ignored and only the effect of differential pressure resistance should be considered;

c) the maximum dimensionless ISW horizontal force Fmax′ has a strong correlation with the relative amplitude ηo/H and water depth ratio h1/h2. Meanwhile, Fmax′ has a notably exponential relation with the comprehensive influence coefficient m, and the goodness of fit R2 of the two is 0.957. Fmax′ shows a positive correlation with ηo/H and a negative correlation with h1/h2.

Acknowledgement: We gratefully acknowledge the valuable cooperation of Prof. Lingling Wang and the members of her laboratory in preparing this Application note .

Funding Statement: This work was supported by the Jiangxi Natural Science Foundation of China [Grant No. 20202BABL214050]; and the Sci-tech Program by Jiangxi Provincial Education Department of China [Grant No. GJJ190942].

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.