Wake Control of Flow Past Twin Cylinders via Small Cylinders

Abstract

The drag and lift force of a twin-cylinder structure are often greater than those of a single cylinder, causing serious structural safety problems. However, there are few studies on the passive control of twin cylinders. The study aimed to investigate the performance of passive drag reduction measures using small cylinders on twin cylinders at a Reynolds number of 100. The effects of small cylinder height (H_D/D = 0~1.0, D is the side length of the twin cylinder) and cross-sectional shape on fluid force and flow structures were studied by direct numerical simulations. The control mechanism was analyzed using high-order dynamic mode decomposition (HODMD). The results showed that significant drag reduction occurred in the co-shedding state, particularly when the gap length of the twin cylinders L/D = 6.0. The small control cylinders with H_D = 0.6, by contrast, showed the best performance in reducing the mean drag and fluctuating lift of the twin cylinders. It reduced the mean drag of the upstream cylinder (UC) by 2.58% and the downstream cylinder (DC) by more than 62.97%. The fluctuating lift coefficient for UC (DC) was also decreased by more than 70.41% (59.74%). The flow structures showed that when the flow hit UC under the action of small control cylinders, a virtual missile-like aerodynamic shape was formed at the leading edge of UC. In this way, the gap vortex consisted of two asymmetric steady vortices and the vortex length significantly increased. This was also confirmed by HODMD. The coherence modes in the gap were suppressed and thus the interaction between gap flow and wake flow was mitigated, which resulted in the fluid force reduction.

1. Introduction

In recent years, there has been rapid development in bridge construction, particularly long-span bridges, many of which are built in typhoon-prone areas. Consequently, the wind resistance of bridges is of paramount importance. Similarly, square cylindrical structures are commonly used in various practical engineering projects, such as high-rise buildings, marine engineering, water conservancy, and hydropower engineering. When airflow moves around a bluff body such as a square cylinder, vortex shedding occurs, leading to periodic fluctuations [1,2,3]. Consequently, the cylinder itself produces fluctuating lift and drag forces, causing periodic vibrations known as flow-induced vibration (FIV). These vibrations [4] encompass a range of complex flow phenomena, such as vortex-induced vibration (VIV) and galloping [5,6]. The phenomenon of fluid force on a cylindrical structure can severely damage the structure and impact structural safety. Previous research has extensively explored methods to reduce fluid force on a single square cylinder. However, in practical engineering, multiple cylindrical structures are commonly used instead of isolated square cylinders. Unfortunately, little research has been conducted on drag reduction for this type of structure. Cylindrical structures experience larger fluid forces than single cylinders, leading to the production of shear layers, wakes, and shed vortices [6]. These phenomena interact with one another, resulting in significantly more complex flow behavior than are observed for single cylinders [7,8]. It is imperative to conduct further research on reducing the fluid force on two or more cylindrical structures, including analyzing the mechanism responsible for large oscillation amplitudes in FIV [9,10]. Doing so will enable us to find effective methods to improve the wind resistance of similar structures, such as double-box girder bridges.

Scholars have conducted extensive research on the hydrodynamic performance of two tandem square cylinders, but the two tandem square cylinder model is more frequently used in practical engineering. For this structure, the gap spacing between cylinders is a critical factor that affects the flow characteristics. Tatsutani et al. [11] proposed that there are three flow regimes for two tandem square cylinders: (i) the single slender-body regime (L/D < 1.5), (2) the reattachment regime (1.5 < L/D < 4), and (3) the co-shedding regime (L/D > 4). This conclusion has been confirmed by scholars in their studies [12]. Sohankar et al. [13] conducted two-dimensional and three-dimensional unsteady simulations to study changes in physical quantities such as lift and drag under different vortex patterns and continuous changes in Reynolds number (Re) and spacing. This study explained the hysteresis phenomenon of physical quantities and pointed out that the hysteresis limit depends on Re. The vortex shedding of the upstream cylinder is obtained by experiment, where there are two modes in which the shedding is suppressed, and the upper branch is related to co-shedding or the binary vortex regime. The results of the upstream cylinder in the tandem arrangement for different Re (L/D = 4) are similar to those of a single cylinder, while those of the downstream cylinder differ.

Flow-induced vibration (FIV) and vortex-induced vibration (VIV) are crucial research subjects for the flow characteristics of two tandem square cylinders. Qiu et al. [14] used numerical simulations to study the influence of variables on the VIV of two tandem square cylinders under the condition of Re = 150 and mass ratios (m*) of 3, 10, and 20, respectively. They reported that the mass ratio has an important influence on the vortex-induced vibration response of the double cylinders, and the mass ratio has a great influence on the distance between the two cylinders. The distance in the synchronous region decreased sharply at m* = 3, but it remained almost unchanged at m* = 10 and 20, indicating that the mass ratio may change the flow pattern. Nepali et al. [15] investigated the impact of Reynolds number (Re) and reduced velocity (Ur) on the vortex-induced vibration (VIV) characteristics of two tandem square cylinders with low mass ratio (m* = 2) and free oscillation. The study was conducted using the characteristics-based-split (CBS) finite element method under the conditions of L/D = 5, Re = 40–200, and Ur = 3–13. The research revealed that the locked region covered a wide range of Ur at lower Re, and a typical 8-shaped trajectory could be observed when Re was low. However, under the influence of wake interaction, the trajectory became irregular as Re increased. Han and his colleagues [16] simulated two tandem square cylinders with Re = 40–200, L/D = 5, and a degree of freedom of 2, observing a stable mode at smaller Re (Re = 40 and 80). With the increase in Re, different states of vortices gradually appeared, and when Re increased to 160 or 200, it was observed that the vortex mode in the locked region changed from 2S to P + S and 2T. The study highlighted the need to accelerate research on drag reduction methods for two tandem square cylinders in order to ensure the stability of this structure.

Numerous researchers have explored methods for controlling the flow field and improving the aerodynamic/hydrodynamic performance of single control cylinders. Such methods can be categorized as active or passive, depending on whether they involve an external power supply. Passive control techniques do not require power, but may necessitate the use of control objects or surface modification. For instance, Gupta et al. [17] adopted passive control techniques and investigated three Re values (75, 100, and 125) by placing a small control cylinder near the wake of single control cylinders of different sizes. The researchers found that this approach could entirely suppress symmetric and asymmetric vortex shedding in the “effective zone.” The effective zone expanded with increasing cylinder size and decreased with increasing Re. Even when Re was large, the drag reduction effect of the cylinder was significant within the complete suppression zone, with the coefficient of drag (C_D) reduced by approximately 10–15% compared to the baseline scenario. Chauhan et al. [18] examined the flow control of a control cylinder tail with a small control rod and discovered that at Re = 485, the highest drag coefficient reduction (22.11%) was achieved when the control rod was located at x/D = 0.5, y/D = ±0.6, while secondary vortices formed behind the control rod caused the drag coefficient to increase when it was positioned at x/D = 1, y/D = ±1. Alonzo-Garcia et al. [19] varied the angle range of the small control rod and concluded that the presence of the rod would result in drag reduction and a decrease in the pressure deficiency behind the cylinder. When the rod was placed on the upper surface (θ = 45–135°), the flow separation did not change due to the gap flow’s impact. However, when the rod’s angle was θ = 135–180°, the reduction in the average recirculation zone increased the pressure loss.

Islam and Manzoor et al. [20] conducted a study on the unsteady laminar flow of a control cylinder using a two-dimensional numerical simulation and passive control of multiple small control cylinders with different gaps and angles around the main cylinder. The study revealed that vortex shedding disappeared completely when Re = 100 and (θ, g) = (450, 0.5), (450, 0.75), (450, 1), and (450, 1.5). Similarly, vortex shedding also disappeared at Re = 160 when (θ, g) = (300, 0.5). It is worth noting that compared to Re = 160, Re = 100 had a more significant inhibitory effect on vortex shedding. The authors [21,22] also utilized a two-dimensional multi-relaxation-time Lattice Boltzmann Method (LBM) with Re = 80–200, cylinder spacing g = 1–7, control plate width w = 0.1–1 d, and under the action of multiple control cylinders, the gap ratio was 0.5–8 and Re was 160, to study the flow of a cylinder with a control plate upstream. The experimental data revealed that the average drag coefficients of the downstream, upstream, and dual configurations were reduced by 8.3%, 51%, and 50.8%, respectively, while the CLrms (root-mean-square value of lift coefficient) values of all three configurations were reduced by 84.4%, 58.2%, and 86.4%, respectively. Using the same approach, Islam et al. [23] studied the effects of l/w and g/D on force reduction, flow separation, and vortex shedding frequency under the conditions of 0.1 ≤ l/w ≤ 1, 1 ≤ g/D ≤ 10, and Re = 150. The study found that when g/D ≥ 6 and l/w ≥ 0.6, the drag reduction effect of the upstream control panel was not evident. However, as the length of the control plate l/w increased, the upstream control plate played a role in slowly suppressing vortex shedding at the downstream position of the cylinder.

Similarly, the research of Zhou and Cheng et al. [24] is akin to that of Islam et al. [21,23]. Their study focused on a Re of 250 and the height of the control panel varied from 10% to 100% of the square side width. The study concluded that the total drag coefficients of the cylinder and the control plate generally decreased with the increase in the height of the control plate. If the height was high enough, the cylinder experienced negative drag. If the height was appropriate, the fluctuating lift could be entirely suppressed.

Coating the outer surface of the main cylinder is also a common passive control method. Hasegawa and Sakaue [25] applied a microfiber coating to the flow-separated region of the cylinder, which was 8% of the cylinder’s length. They studied the influence of the coating on the wake resistance and eddy current under the subcritical condition of Re of 6.1 × 10⁴, taking the different angle positions of the coating at the leeward side of the cylinder as variables. They found that the coating angular position affected the vortex formation length and resistance and the main wake frequency corresponding to the vortex formation. The experimental data indicated that passive drag reduction measures could reduce resistance by 16%. Moreover, the formation length of the vortex was extended and the main wake frequency increased. Additionally, Du and Zhang et al. [26] conducted wind tunnel experiments at subcritical Re and applied porous materials to the trailing edge of a smooth cylinder for force measurement based on aerodynamic balance. Through PIV flow field visualization and transient vorticity field DMD analysis, they demonstrated that this passive method could effectively suppress wake vortex strength and vortex shedding, thereby reducing wake drag and enhancing vortex shedding stability. The drag reduction effect increased with increases in Re and the permeability coefficient, and DMD analysis showed that the porous material coating could effectively weaken the energy of different modes of vortices and suppress the generation of large-scale vortices.

The flow characteristics of a cylinder with slits are also an area of focus in electronic engineering technology. The shear stress plays a leading role in fluctuating lift while the pressure on the outer surface dominates drag force. Hus and Chen [27] classified the cylinder into two flow patterns based on the inclination angle: injection and suction. The tangent slit had a higher vortex-shedding frequency than the oblique slit, with the frequency increasing as the distance between the vortex and the cylinder decreased. In addition, in flow states with high Re, the drag increased with the increase in slit width, whereas the drag value was small in the flow of an inclined slit cylinder with an inclined angle less than 45° and low Re.

Chauhan et al. [28] employed an attached splitter plate with a length of cylinder width of L/D = 0–6 to regulate the flow of a control cylinder at a fixed Re = 485. It was observed that the time-averaged resistance coefficient continued to decrease when L/D = 2–4, increased slightly when L/D = 5, but decreased again when L/D = 6. Furthermore, the formation time of the vortex increased as the length of the plate increased. When L/D = 2, the maximum reduction of the Strouhal number was observed to reach 33.56%. However, when L/D = 3, a strong secondary vortex was generated, and the Strouhal number increased suddenly until vortex shedding began to be suppressed when L/D > 6. This conclusion was further confirmed in subsequent experiments [29]. Sharma et al. [30] found that as the flexibility of the flexible foil attached to the aluminum square cylinder increased under the same length, the liquid–solid coupling became complicated, and it was proven that the swing amplitude of the foil had no obvious relationship with the drag coefficient, while the wake frequency and the locking of the inherent mode would affect it.

Active control has a promising future as it can adjust the flow characteristics based on the environment. Saha et al. [31] conducted numerical simulations on a control cylinder with Re = 100 and pointed out that the control performance decreased with the increase in blowing in the form of a jet. Huang et al. [32] employed rear jet injection technology to actively control the flow and aerodynamic characteristics of the cross-flow of a control cylinder. The article mentioned that the drag coefficient decreased with the increase in the injection ratio, and four characteristic flow modes, namely wake-dominated, transitional, critical, and jet-dominated, were produced at low injection ratio (IR), medium-low IR, moderate-high IR, and large IR, respectively.

Qu et al. [33] employed a slot synthetic jet to investigate its control effect on the front surface of a square cylinder. The results indicated that the high excitation frequency of the near-wake flow (f_e/f₀ = 3.4) enabled the synthetic jet deflection to effectively suppress airflow separation, which was not possible with low excitation frequency. Additionally, the authors proposed that the synthetic jet placed on the front surface could more effectively suppress the wake shedding of the square cylinder compared to the continuous jet on the front surface, based on experimental comparisons. Sohankar et al. [34] also found that at an incident angle α = 45°, Re = 150, and |Γ| = 0.4, the drag coefficient was reduced by approximately 39% compared to |Γ| = 0 and the vortex shedding phenomenon was suppressed, achieving the best experimental outcome.

Previous studies mainly focused on the drag reduction of a single cylinder, while very few studies concentrated on twin cylinders [35]. Moreover, the flow mechanism of a single cylinder is very different from that of twin cylinders, and the control methods for a single cylinder may not be suitable for twin cylinders. As far as we know, only suction control has been applied to perform drag reduction for twin cylinders at present. This method is an active control, which requires the input of external energy, and it requires the help of blow inhalation equipment or an air pump. This study explored the lift and drag of twin cylinders in series under the condition of co-stripping. The suction flow was set on both sides near the upper cylinder (UC). The results showed that the drag reduction of the lower cylinder was greater than that of the upper cylinder, and the lift values of the twin cylinders were similar [36]. So far, no publication has applied the passive method to two-cylinder drag reduction; therefore, this paper placed a small control cylinder above each twin cylinder to explore its drag reduction effect. The advantage of this passive control is that there is no need for external energy input, just putting the small control cylinder at a certain height can reduce the drag and fluctuating lift of the twin cylinders. This paper aimed to study the details of this passive control method. It is worth noting that the results are beneficial for the control of fluid-induced vibration (FIV), as this method can reduce vortex shedding, which is the main reason for introducing FIV [16,37].

This paper aimed to explore the drag reduction effect of passive control on two tandem square cylinders. The study utilized a small control cylinder (d = 0.1) placed perpendicular to the longitudinal direction of the main cylinders at a certain height above them. The height variation range was locked at h/D = [0–1.0]. When Re = 100, the gap between the two main cylinders was varied as L/D = 2.0, 3.0, 4.0, 5.0, and 6.0 for numerical simulations in order to investigate changes in lift, drag, and pressure cloud diagrams under different flow conditions. A set of working conditions with the smallest fluctuating lift value and average drag value was then selected. Then, the small cylinders with a circular section shape were modified into square, octagonal, and decagonal shapes, respectively, and numerical simulations were carried out to determine the influence of the control column with different sections on the convective dynamic and fluid structure at different heights, as well as reveal certain laws.

It is worth noting that at present, no publication has reported passive drag reduction experiments applying small control cylinders above the twin cylinders, and this study also added the variables of small control cylinder height and cross-sectional shape on this basis. More importantly, this paper has high practical engineering application value and can provide a theoretical basis for mitigation of wind-induced effects. For example, because the main cylinders resemble a twin-box girder bridge and the cross-section of the small control cylinder resembles the auxiliary facility railing on the bridge, this research explores the influence of railing height and cross-sectional shape on the wind resistance performance of the bridge. These findings are highly relevant to practical engineering applications, such as developing strategies to mitigate the adverse effects of wind load on bridge structures by altering the flow field, improving the selection of railing cross-sectional shapes for future wind-resistant bridge design, and enhancing measures to control wind loads [38,39].

2. Materials and Methods

2.1. Control Equation

Assuming the system is two-dimensional, unsteady, incompressible, laminar, and Newtonian flow, the control of such fluid flow is represented by the two-dimensional Navier–Stokes equation. The governing equation for flow past two tandem cylinders can be expressed by:

$$\nabla \cdot \mathit{u}=0$$

(1)

$$\frac{\partial u}{\partial t}=-(u\cdot \nabla )u-\nabla p+\frac{1}{Re}{\nabla }^{2}u$$

(2)

Note that this equation serves as a crucial tool for investigating the control of fluid flow in a two-dimensional, unsteady, incompressible, laminar, and Newtonian setting, allowing for the analysis of various quantities such as drag and lift coefficients, pressure coefficients, and Strouhal number. Please note that the following variables are non-dimensional. Where u = (u_x,u_y) is the two-dimensional flow velocity vector, and p is the fluid pressure. U represents the free flow velocity, and D represents the width of the cylinder, which is used to represent all the velocity components and geometric length scaling. Additionally, the static pressure and time increase coefficients are adjusted to ρU² (ρ represents fluid density) and D/U, respectively. The overall dimensionless quantities in this formula include the drag coefficient (C_D), lift coefficient (C_L), pressure coefficient (Cp), and Strouhal number.

$$\sigma (p,u)=-pI+\frac{1}{Re}\left(\nabla u+\nabla u^{\mathrm{T}}\right)$$

(3)

$$F_D=2{\displaystyle {\oint }_{}\{\sigma (}p,u)\cdot n\}\cdot e_{x}dl$$

(4)

$$F_L=2{\displaystyle {\oint }_{}\{\sigma (}p,u)\cdot n\}\cdot e_{y}dl$$

(5)

$$C_L=\frac{F_L}{0.5\rho U^{2}A}$$

(6)

$$C_D=\frac{F_D}{0.5\rho U^{2}A}$$

(7)

$$St=\frac{fD}{U}$$

(8)

$$C\mathrm{p}=\frac{P-P_{\infty }}{0.5\rho U^2}$$

(9)

The above formulas represent the total drag and lift in the x and y directions as F_d and F_L, respectively. P_∞ represents the free pressure, σ is the stress tensor, and n is the outward-facing unit normal at the cylinder surface. The frequency of associated vortex shedding and Strouhal number can be determined by the fast Fourier transform of fluctuating lift (f, St). The formulas for calculating the average drag coefficient ($\overline{C_D}$), average lift coefficient ($\overline{C_L}$), root-mean-square (RMS) lift coefficient (C_L’), average pressure coefficient ($\overline{C\mathrm{p}}$), and root-mean-square pressure coefficient Cp’ are given below:

$$\overline{C_D}=\frac{1}{N}{\displaystyle \sum _1^N{C}_D}$$

(10)

$$\overline{C_L}=\frac{1}{N}{\displaystyle \sum _1^N{C}_L}$$

(11)

$$C_L{}^{\prime }=\sqrt{\frac{1}{N}{\displaystyle \sum _1^N{(C_L-\overline{C_L})}^2}}$$

(12)

$$\overline{C\mathrm{p}}=\frac{1}{N}{\displaystyle \sum _1^{N}C\mathrm{p}}$$

(13)

$$C{\mathrm{p}}^{\prime }=\sqrt{\frac{1}{N}{\displaystyle \sum _1^N{(C\mathrm{p}-\overline{C\mathrm{p}})}^2}}$$

(14)

Note: N is the number of vortex shedding cycles, greater than or equal to 40.

2.2. Computational Domain and Boundary Conditions

The computational domain in this study consisted of twin cylinders in series, with Re = UD/ν, where D is the cylinder width, ν is the viscosity, and Re is set to 100. To vary the Reynolds number, the viscosity was changed while keeping U and D constant. Five double cylinders with different dimensionless gap ratios (L/D = 2, 3, 4, 5, and 6, where L is the distance between the center points of the twin cylinders) were studied under uniform flow (U). In addition, two small control cylinders were applied at a certain height above the two cylinders, with the dimensionless small circular cylinder height (H_D) defined as H_D = h/L, (H_D = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0), where h is the height from the center of the small control cylinder section to the upper edge of the cylinder. The study investigated the influence of changing the cross-sectional shape of the small circular cylinders on drag reduction under the optimal drag reduction conditions of the small circular cylinders. The variable of cross-sectional side length is defined as Edge, with the circular cross-section approximated by Edge = 32, as shown in Figure 1b. ANSYS-Fluent 17.0 was used to solve the simulated N–S equation based on the finite volume method, with the third-order MUSCL scheme used to discretize the convective component in order to ensure second-order time accuracy.

The computational domain for the small circular cylinders is depicted in Figure 1a. The domain is defined in the Cartesian coordinate system, with the origin located at the center of the twin cylinders. The distance from the inlet to the origin was 23.5D, while the distance from the outlet to the origin was 76.5D. The upper and lower boundaries were symmetrically placed at 23.5D from the centerline. Figure 1c provides a two-dimensional view of the non-structured grid distribution around the cylinder throughout the computational domain. The low exponential rate stretching grid was employed in the far field to optimize and improve the computational efficiency, and the grid density around the cylinder was increased to enhance the calculation accuracy. For the isogonal skew, according to the histogram of the isogonal skew value automatically generated by ICEM, we found that 72.227% of the mesh mass was 0.90–1.0, the lower the mesh quality, the smaller the proportion, and all of the mesh mass was above 0.45. The laminar model was used for the overall viscosity. At the inlet, Dirichlet boundary conditions (u_x = U, u_y = 0, ∂p/∂x = 0) were applied, while zero gradient velocity and the pressure outlet (∂u/∂x = 0, p = 0) were set at the outlet of the computational domain. The upper and lower sides were the same as the inlet. No-slip conditions (u = 0, ▽p • n = 0, n is the normal vector of the wall) were applied on the wall of the cylinder.

2.3. Higher-Order Dynamic Mode Decomposition

The Koopman operator is an infinite-dimensional linear operator used to represent the time evolution of nonlinear dynamic systems [40,41]. By decomposing complex systems into series modes, it can be used to represent the periodic components of complex nonlinear systems [41]. Dynamic mode decomposition (DMD) is a technique that utilizes the Koopman theory to linearly approximate high-dimensional dynamic systems. There have been many improved versions of DMD, with HODMD being one such improved version that has achieved satisfactory results in handling complex flow data and noise [41].

Under the Koopman assumption, the matrix can be mapped to another matrix by introducing a linear operator a, which can be expressed as: Q_N−1, Q_N.

$$A{Q}_{N-1}=Q_N+R$$

(15)

$${\mathrm{Q}}_{\mathrm{N}-1}=\left[{\mathrm{q}}_1,{\mathrm{q}}_2,{\mathrm{q}}_3,\cdots ,{\mathrm{q}}_{\mathrm{N}-1}\right]$$

(16)

$${\mathrm{Q}}_{\mathrm{N}}=\left[{\mathrm{q}}_2,{\mathrm{q}}_3,{\mathrm{q}}_4,\cdots ,{\mathrm{q}}_{\mathrm{N}}\right]$$

(17)

Linear operators are commonly used in mathematics to describe the relationship between two vector spaces. When working with linear operators, it can be useful to approximate their behavior in order to gain insights into their properties. One way to do this is to use the singular value decomposition method (SVD). The SVD method allows us to obtain an approximation of linear operator A by decomposing it into three matrices: U, Σ, and V*. Here, R denotes the residual matrix, which represents the difference between the original operator A and its approximation. Using SVD, we can express the eigenvalues of A in logarithmic form, which can provide useful information about the behavior of the operator. By analyzing the eigenvalues, we can determine important properties such as the rank and null space of the operator [42].

$${\lambda }_{\mathrm{i}}=\frac{\log {\sigma }_{\mathrm{i}}}{\mathsf{\Delta }\mathrm{t}}$$

(18)

where ${\sigma }_i$ represents the eigenvalues of the matrix M. For the developed flow state, ${\lambda }_i$ represents an imaginary part, which denotes the frequency of the mode. The corresponding DMD mode ${\psi }_i$ is defined as:

$${\psi }_{\mathrm{i}}={\mathrm{Ly}}_{\mathrm{i}}$$

(19)

where ${\mathrm{y}}_{\mathrm{i}}$ is the eigenvector of the matrix M. $a_i$ is the mode amplitude, which is defined as:

$${\mathrm{a}}_{\mathrm{i}}={\left({\psi }^{-1}\right)}_{\mathrm{ij}}{\left({\mathrm{q}}_1\right)}_{\mathrm{j}}$$

(20)

The mode energy $E_i$ can be expressed as:

$$E_i=\left|{\mathrm{a}}_i{\psi }_{\mathrm{i}}\right|$$

(21)

The mth flow field can be reproduced by using:

$${\mathrm{q}}_m={\displaystyle \sum _{i=1}^n{\mathrm{a}}_i}{\psi }_i{\lambda }^{m-1}$$

(22)

The standard DMD may not be suitable for highly nonlinear systems, particularly when there is a mismatch between spectral and spatial complexity. To address this limitation, HODMD was developed as a supplement for highly nonlinearized systems [41]. The underlying theory of HODMD is Taken’s delay embedding theorem, which augments the dimensionality of observed variables. With the assistance of HODMD, the observation original matrix ${\mathrm{V}}_1^{\mathrm{K}}$ is extended to matrix ${\mathrm{V}}_1^{\mathrm{K}-\mathrm{d}+1}$ with a desirable dimension by adopting d index-lagged snapshots. Note that ${ \mathrm{V}}_1^{\mathrm{K}-\mathrm{d}+1}$ consists of sub-matrices formed with the first and last K-d columns of the modified snapshot matrix. Then, ${\mathrm{V}}_1^{\mathrm{K}-\mathrm{d}+1}$ is used as input data, while the remaining steps are consistent with DMD, which is referred to as HODMD in this study.

$${\mathrm{V}}_1^{\mathrm{K}}=\left[{\mathrm{v}}_1,{\mathrm{v}}_2,\cdots ,{\mathrm{v}}_{\mathrm{k}},{\mathrm{v}}_{\mathrm{k}+1},\cdots ,{\mathrm{v}}_{\mathrm{k}-1},{\mathrm{v}}_{\mathrm{k}}\right]$$

(23)

$${\overline{\mathrm{V}}}_1^{\mathrm{K}-\mathrm{d}+1}\equiv \left[\begin{array}{c}{\mathrm{V}}_1^{\mathrm{K}-\mathrm{d}+1}\\ {\mathrm{V}}_2^{\mathrm{K}-\mathrm{d}+2}\\ \cdots \\ {\mathrm{V}}_{\mathrm{d}-1}^{\mathrm{K}-1}\\ {\mathrm{V}}_{\mathrm{d}}^{\mathrm{K}}\end{array}\right]$$

(24)

It is worth noting that while orthogonal decomposition (POD) is a conventional and widely used technique for extracting coherent modes, it has some known limitations. Specifically, POD is a second-order statistical mode decomposition method, and the main direction in a dataset may not necessarily correspond to the direction of dynamic importance, leading to potentially inaccurate coherent mode extraction [42]. In contrast, DMD and its variant HODMD offer richer physical insight and can more accurately capture coherent modes. Thus, for the analysis of coherent modes under complex flow structures, this study utilized HODMD.

2.4. Grid Independence and Verification

After investigating grid independence and calculation time steps, an appropriate number of grids and calculation times were chosen. The following study focused on calculating and comparing the average drag coefficient (C_D), fluctuating lift coefficient (RMS) (C_L’), and Strouhal number (St) of twin cylinders with previously published literature. The study was conducted using approximately 50,000 control grids, with the assumption that Re = 100 and Δt = 0.005. Additionally, a control group consisting of 80,000 grids was created with Re = 100 and Δt = 0.002.

Compared to the literature, the fine grid ensured that the influence of the number of grids on the calculation results was minimal, which improved the efficiency and accuracy of the calculation, as shown in

Table 1. Test of the independence of the mesh and time step for twin cylinders (Re = 100, L/D = 2).

rth Order	Δt	CL’(UC)	CL’(DC)	St
5	0.005	0.319	1.249	0.130
8	0.002	0.316	1.265	0.125
Zhang et al. (2022) [43]	0.00211	0.290	1.230	0.133
Bao et al. (2012) [44]	0.005	0.285	1.196	0.129

. The results demonstrated that the simulation outcomes with Δt = 0.005 and 0.002 showed good consistency under the two different grids compared to the results of Zhang [43] and Bao et al. [44].

In Figure 2, the standard deviations of the lift coefficients for UC and DC were compared with the results from Zhang et al.’s paper [41,44] in 2022. The study found that in a narrow gap (L/D = 2–4), the fluctuating lift coefficient only slightly changed. However, as the flow transitioned from the reattachment state (L/D = 2–4) to the co-shedding state (L/D = 5–6), the fluctuating lift coefficient sharply increased between L/D = 4 and 5. In the co-shedding state, the fluctuating lift coefficient for UC significantly decreased while the lift coefficient for DC significantly increased. These simulation results were in good agreement with Zhang et al.’s findings, which confirmed the accuracy of the numerical simulation in this paper.

To further verify the accuracy of the mesh independence, we conducted tests using both coarse and fine mesh, and the results are listed in Table 1. For a detailed description of the verification and validation of the cylinder’s forced vibration, we recommend referring to Zhang et al.’s previous research.

Table 1 presents the results of the mesh independence test using both coarse and fine grids with Re = 100 and L/D = 2. We compared the Strouhal number (St = fD/U), lift coefficient (C_L), and standard deviation of the lift coefficient (C_L’) with those of previous studies while varying the order r. The force coefficient was reduced to 1/2ρU²D, where ρ represents the air density and U represents the free air velocity.

As demonstrated in Table 1, the results from the two grids showed slight differences, but were still consistent with the results of Zhang et al. [43] and Bao et al. [44].

3. Results and Discussion

3.1. Mean Drag and Fluctuating Lift Forces When H_D Is a Variable

Figure 3 depicts the mean drag and fluctuating lift coefficients of the two cylinders as a function of the height of the small circular cylinders (H_D) with a circular section under the conditions of Re = 100 and L/D = 2 and 6. A shown in Figure 3a, the mean drag coefficient for DC was lower than that for UC due to the shielding effect. Within the range of 0.2 ≤ H_D ≤ 0.6, the drag coefficient for UC decreased as H_D increased, while the drag coefficient for DC increased as H_D increased. Outside this range, the drag coefficient for UC increased as H_D increased, whereas the drag coefficient for DC decreased as H_D increased. The drag of UC decreased by over 2.84%, while the drag of DC decreased by 35.20%. A shown in Figure 3b, the fluctuating lift coefficient for DC was not significantly different from that for UC for the baseline case. Except for 0.4 ≤ H_D ≤ 0.6, the fluctuating lift coefficients for UC and DC increased as H_D increased. The fluctuating lift of UC decreased by over 77.92% and that of DC decreased by over 77.60%. The minimum lift of both cylinders was at H_D = 0.6. In summary, H_D had an insignificant effect on drag but significantly reduced the fluctuating lift.

As shown in panel (c) of Figure 3, under the condition of L/D = 6, the mean drag coefficient for DC was found to be smaller than that for UC due to the shielding effect. In the range of 0.4 ≤ H_D ≤ 1.0, the drag coefficients for both UC and DC followed a similar trend, showing a decrease followed by an increase. When 0.6 ≤ H_D ≤ 1.0, they both stabilized at approximately 1.4 and 0.5, respectively. In the range of 0.0 ≤ H_D ≤ 0.4, the mean drag coefficient for UC first increased and then decreased, while that for DC showed the opposite trend. The drag coefficient for UC decreased by over 2.58%, while that for DC decreased by 62.97%. Panel (d) of Figure 1 shows that, compared to the single cylinder case, the gap flow between the two cylinders promoted severe wake vortex shedding, resulting in a significantly higher fluctuating lift coefficient for DC. In the range of 0.2 ≤ H_D ≤ 1.0, the fluctuating lift coefficients for both UC and DC followed the same trend, with the smallest value occurring at H_D = 0.6. In the range of 0.0 ≤ H_D ≤ 0.2, the trend was the opposite. Compared to the baseline case, the fluctuating lift coefficient for UC was reduced by more than 59.74%, while that for DC was reduced by over 70.41%. These results demonstrated that H_D not only reduced drag but also significantly reduced pulsating lift.

3.2. Edge as a Variable

Figure 4 shows the mean drag and fluctuating lift coefficients generated by twin cylinders under the action of small control cylinders with different sides (the inscribed circle diameter is 1). As shown in (a), the drag coefficient for UC was consistently much larger than that for DC, regardless of the number of cross-sectional edges. Notably, for UC, the drag coefficient remained relatively stable at around 1.4, whereas for DC, it significantly increased from Edge = 4 to 8, and then it stabilized at about 0.4 for larger values of Edge. Additionally, when Edge = 4 in (b), the fluctuating lifts of UC and DC were the smallest, with only a small difference between them. As Edge increased, the fluctuating lift coefficient for UC slightly increased while that for DC significantly increased. Ultimately, they stabilized at around 0.1 and 0.45, respectively.

In summary, when Edge = 4, the mean drag and fluctuating lift coefficients of the two cylinders were the smallest, leading to the best drag reduction effect. When Edge = 8, 16, and 32, the drag reduction effect was not significantly different.

3.3. Re as a Variable

As can be seen from Figure 4, when the cross-section of the small control cylinders was square, the drag reduction effect was the most ideal. Therefore, this part of the study explored the effect of the Reynolds number on the height of the small control cylinders and its effect on the drag reduction of two cylinders. The mean drag and fluctuating lift coefficients for UC and DC at different Re are presented in Figure 5. In the baseline case (a), except for the range of 75 ≤ Re ≤ 125, the drag coefficients for both UC and DC exhibited the same trend, with C_D increasing as Re increased. Under the influence of small control cylinders, the C_D for UC and DC increased in the range of 100 ≤ Re ≤ 175, with the C_D for DC showing a greater change. At Re = 200, C_D significantly decreased. As for fluctuating lift in (b), except for the decrease in the lift coefficient for DC at Re = 125, the C_L’ for UC and DC increased with an increase in Re. Under the influence of small control cylinders, C_L’ significantly decreased, with the C_L’ for DC showing the greatest decrease. The C_L’ for both UC and DC showed an increase with an increase in Re. The change under the influence of small control cylinders was similar to the baseline case, resulting in little change in the drag reduction effect in this Reynolds number range.

3.4. Time History of Fluid Force Coefficients

Figure 6 depicts the time history of the fluid force coefficients for UC and DC under the influence of small circular cylinders at varying heights. It was observed that the primary drag of UC did not show a significant change with the increase in the height of the control cylinder, whereas the drag of DC decreased to varying degrees under the influence of control cylinders of two different heights. Furthermore, the fluctuating lift of both UC and DC decreased under the control effect, with the drag reduction effect of DC being particularly noteworthy. Interestingly, the smallest lift coefficients were observed for both UC and DC at H_D = 0.6.

Figure 7 displays the time history of the fluid force coefficients for UC and DC under different heights of small control cylinders, while Figure 8 shows the power spectral density (PSD) corresponding to the lift coefficient. It was evident that H_D = 0.6 was the optimal condition for drag reduction. With an increase in the height of the small control cylinders, the drag on UC and DC initially decreased and then increased, as compared to the baseline conditions. From Figure 8, it can be observed that the primary peak of PSD in other working conditions decreased while the super harmonic peak remained constant or even increased, as compared to the basic working conditions. Due to the nonlinear nature of the dynamics, the vortex shedding was suppressed, but the vortex of the small-scale structure was enhanced. When H_D = 0.6, the vortex shedding frequency significantly decreased and weak vortex shedding increased. As a result, the lift coefficients of the twin cylinders considerably decreased, with the lift coefficient for UC showing a larger decrease.

3.5. Lift Coefficients and Drag Coefficient Phase Diagram

Figure 9 displays the phase diagrams of C_D and C_L between UC and DC. Under the basic condition, it was observed that the C_L for UC and DC gradually increased with time, while the C_D diagram showed evident limit cycle oscillation, representing a negative correlation state after the oscillation was stable. When H_D = 0.6, the C_D image spiraled upward with time and the spiral amplitude gradually increased. The C_L diagram for UC and DC had two running trajectories, and after adding the railing, the positive correlation trajectory of UC and DC was altered, indicating some instability. For the H_D = 1.0 condition, the image exhibited some periodicity and regularity. Once the C_L image was stable, it displayed an elliptical shape with the left upper right lower diagonal as the long axis, signifying a negative correlation state.

3.6. Instantaneous Flow Structure

Figure 10 displays the instantaneous flow structures for the baseline case and H_D = 0.6 and 1.0. In the absence of small control cylinders, Karman vortex shedding occurred in the gap and wake, forming a 2S mode. Notably, the gap vortex generated by UC interacted strongly with DC, leading to a strong fluctuating lift and reinforcing the wake vortex behind it. This was consistent with the findings of Bhatt R et al. [45]. A small eddy current gradually separated below DC. Under the control of H_D = 0.6 and 1.0, a virtual aerodynamic shape resembling a missile head formed at the flow separation of the UC leading edge, and the flow above UC was smoothed by the action of the small control cylinders. Compared to the baseline case, the vortex at H_D = 0.6 increased but did not fall off over time, forming a steady vortex that remained stable. Additionally, when the flow hit the small control cylinders above DC, the flow trajectory significantly changed, and the flow passed under DC. As a result, the wake became longer and narrower, and the drag was significantly reduced.

In order to provide further insight into the drag reduction mechanism of the small control cylinders, Figure 11 presents a detailed flow field analysis around the cylinder. For the baseline case, the flow was observed to separate at UC, leading to the formation of two small vortices below UC and DC, respectively, which eventually separated over time. By employing H_D = 0.6 small control cylinders, a virtual missile-like aerodynamic shape was formed at the leading edge of UC, which guided the flow at the separation region. The flow width decreased significantly when passing through UC, and the small control cylinders compressed the fluid’s flow, thereby inhibiting its upward trend and resulting in smoother flow. Additionally, the small control cylinders induced a steady vortex to form at the gap, which was significantly increased in comparison to the shedding vortex observed in the baseline case. Furthermore, a stable vortex was formed at the trailing edge of DC. Consequently, the waking edge was longer, the vortex became weaker, the wake became narrower, and the drag was significantly reduced.

3.7. Average Streamwise Velocity Contour Line

Figure 12 presents contour maps of the mean streamwise velocity at different small control cylinder heights at baseline and H_D = 0.6. The vortex formation length of a bluff body strongly depends on the fluid force, as pointed out by Griffin [46]. Specifically, the time-averaged vortex formation length Lf is defined as the distance between the center point of the leeward side and the point with the average velocity u = 0, which is the saddle point of the wake. The presence of the small control cylinders significantly influenced the flow. In the baseline case, the upper and lower vortices near the trailing edge of UC and DC were symmetrically distributed. However, as the height of the small control cylinders increased from H_D = 0.0 to 0.6, the fluid passing through the upper edge of UC and DC became smoother. This led to a significant increase in the length of the vortex at the trailing edge of UC, but the two vortices were not symmetrically distributed. Similarly, the wake after DC became significantly longer and the trailing edge did not show two symmetrical vortices, only forming a certain constant vortex near the trailing edge. As a result, the drag of UC and DC was significantly reduced, as shown in Figure 3.

In conclusion, the results of this study demonstrated that the height of the small control cylinders had a significant impact on the vortex formation length and the drag of the small control cylinders. The findings provide useful insight for the design and optimization of bluff body configurations.

3.8. Coherence Modes

We utilized the HODMD technique to examine the flow structures in two scenarios, namely baseline and H_D = 0.6. Our dataset comprised 40 complete cycles of vortex shedding, with each cycle consisting of 25 snapshots. It is worth noting that all of the snapshot data were captured post completion of the flow. Consequently, the primary Ritz values were distributed uniformly on a unit circle. In this paper, we delved into the dominant coherent mode and global mode energy, and we analyzed their relationship with the fluid forces.

Figure 13 depicts the correlation between energy and frequency. This paper computed the vorticity and streamline modes based on the velocity mode. It was evident from Figure 13 that the highest modal energy in both cases was always observed when fD/U = 0, corresponding to the time-averaged mode, where M0 denotes this mode. The first two dynamic modes are determined based on their non-dimensional frequency fD/U (M1: fD/U = St; M2: fD/U = 2St). For the baseline scenario (H_D = 0.0), Figure 13a shows that the wake comprises a dominant mode (M1) and a super harmonic mode (M2). The energy of the super harmonic mode was significantly lower than that of the primary mode, consistent with the PSD results illustrated in Figure 8. As the height of the small control cylinders increased to H_D = 0.6, the average modal energy remained unchanged, while the modal energy of M1 and M2 significantly decreased in comparison to the baseline case.

Concerning the baseline scenario in Figure 14(a1), the primary mode M1 (dynamic mode), resembling Kármán vortices, emerged in the gap and wake and was symmetrical around the cylinder’s centerline. This mode was closely associated with lift. The coherence modes in the gap were suppressed under the controlled case, and thus the interaction between gap flow and wake flow were mitigated. Concerning the super harmonic mode M2 depicted in Figure 14(b2), it appeared to be half the size of the primary mode M1. The coherence mode in the gap was also suppressed. In summary, as the small control cylinders’ height increased, the global energy of the average mode remained unchanged while the local average mode’s pattern changed. Conversely, the dynamic modal energy significantly decreased, whereas the modal shape remained unaltered, which resulted in the fluid force reduction.

4. Conclusions

The drag and lift force of a two-cylinder structure is often greater than that of a single cylinder, causing serious structural safety problems. However, there are few studies on the passive control of twin cylinders. This paper presents a study on the drag reduction performance of twin cylinders via small cylinders in a laminar, two-dimensional, and unsteady flow with Re = 100. The high-order spectral element method was used for the simulations, while high-order dynamic mode decomposition was employed to analyze the dynamic coherent mode and the underlying control mechanism. The study aimed to investigate the impact of small control cylinders with varying heights on the fluid force and flow structure.

Regarding the gap distance, the results showed that the mean drag force of UC was significantly reduced under the influence of the small control cylinders, except for L/D = 2.0. When the flow was in the reattachment state (L/D ≤ 4.0), the drag and fluctuating lift of both UC and DC were negligible, and the control effect was limited compared to other clearance conditions. In the co-shedding state (L/D = 5.0 and 6.0), the drag reduction effect was more significant, and the fluctuating lift of both cylinders was reduced to varying degrees. The optimal drag reduction effect was achieved when L/D = 6.0.

When the small control cylinders were activated, the suction airflow hit the upstream cylinder (UC) and formed a virtual missile-like aerodynamic shape at its leading edge, which guided the flow and reduced flow separation. This resulted in a smoother flow above UC and an increase in the width of flow separation. The gap vortex, which is typically generated by the vorticity from the separated shear layer, was weakened and shrank due to the formation of the aerodynamic shape at the leading edge of UC. Consequently, two asymmetric steady vortices were formed at the gap between the two cylinders, resulting in a significant increase in vortices. The smaller gap vortex only merged with the vortex generated by the downstream cylinder (DC) to a limited extent, reducing the interaction between the gap flow and DC. This caused the wake of UC and DC to become longer, narrower, and weaker, resulting in reduced drag and lift.

The study showed that the fluid force applied to UC and DC was significantly affected by the four small control cylinders’ different heights and cross-sectional shapes. The use of small control cylinders with a height of H_D = 0.6 resulted in a mean drag reduction of more than 2.58% for UC and 62.97% for DC. Furthermore, the fluctuating lift coefficients for UC and DC decreased by more than 70.41% and 59.74%, respectively. It is noteworthy that the control effect remained relatively stable within a wide range of Re (Re = 75–200), indicating that the influence of the Reynolds number on drag reduction was limited.

The results from higher-order dynamic mode decomposition (HODMD) revealed that a main mode (M1) and a super harmonic mode (M2) dominated the wake. The influence of the small control cylinders on the local time-averaged mode became more pronounced compared to its effect on the global time-averaged mode energy. In contrast, the dynamic modal energy was significantly reduced, while the modal shape remained unchanged except for a phase shift. In addition, for the controlled case, the coherence modes in the gap were suppressed and thus the interaction between the gap flow and wake flow were mitigated, which was responsible for the fluid force reduction.

The focus of this study was to reduce the lift and drag of twin cylinders with small control cylinders of variable heights under laminar flow and low Reynolds number conditions. This research result can be applied to the drag reduction of structures such as dual-box girders, twin towers, and neighbor buildings in practical engineering applications. We will conduct future studies on the drag reduction of these engineered structures using the present method. In addition, the drag reduction under high Reynolds number turbulence is also of great interest in our future studies.