Figures 15–20 show the mean and fluctuating aerodynamic force coefficients of the twin-box girder under the aerodynamic interference by leading circular cylinders at different spacing ratios. Figure 15 shows that the time-averaged drag force coefficient of the twin-box girder with cylinder interference decreased significantly compared with the original force coefficient without interference, because the shedding vortices in the wake of the circular cylinder destroyed the original laminar boundary layers; hence, the laminar boundary layers transition to turbulence at the windward corner, and the turbulent boundary layers can be maintained throughout the girder section. Thus, the boundary layer is composed of many small vortices, which effectively reduce the contact area between the fluid boundary layer and the twin-box girder, leading to a decrease in the friction between the turbulent boundary layer and the model. Meanwhile, when the vortices contact the girder, the velocity direction is generally opposite to the incoming flow. In this way, the vortices generate a force that is opposite to the streamwise direction, and the force generated by the interaction between the vortices and the girder reduces the total drag force of the bridge deck. Meanwhile, when the diameter of the cylinder is 3 mm, it clearly shows that the time-averaged drag force coefficient exhibits a periodic trend at different locations from the model, as shown in Figure 15, which indicates that the vortex scale is closer to the diameter. In addition, the time-averaged drag force with cylinder interference at low Re achieves a similar characteristic of undisturbed drag force at high Re.

**Figure 15.** (**a**–**d**) The relationship between the time-averaged drag force coefficient and the Reynolds number at different spacing ratios and cylinder diameters.

**Figure 16.** (**a**–**d**) The relationship between the fluctuating drag force coefficient and the Reynolds number at different spacing ratios and cylinder diameters.

**Figure 17.** (**a**–**d**) The relationship between the time-averaged lift force coefficient and the Reynolds number at different spacing ratios and cylinder diameters.

**Figure 18.** (**a**–**d**) The relationship between the fluctuating lift force coefficient and the Reynolds number at different spacing ratios and cylinder diameters.

**Figure 19.** (**a**–**d**) The relationship between the time-averaged moment force coefficient and the Reynolds number at different spacing ratios and cylinder diameters.

**Figure 20.** (**a**–**d**) The relationship between the fluctuating moment force coefficient and the Reynolds number at different spacing ratios and cylinder diameters.

Figure 16 shows that the fluctuating drag force coefficient of the twin-box girder with cylinder interference is much smaller than that of the twin-box girder without interference at Re ≤ 9.0 × 103. Moreover, the fluctuating drag force is insensitive to the Reynolds number, and the disturbed fluctuating drag force at low Re is closer to that of the undisturbed fluctuating force at high Re. Furthermore, the fluctuating drag force is independent of the spacing ratio *ε*.

Figure 17 shows that the time-averaged lift force coefficient of the girder with cylinder interference also decreased sharply compared with that of the bare twin-box girder. It should be noted that with the increase in distance between the cylinder and the twin-box girder, the *CL mean* shows less sensitivity to the Reynolds number, because the turbulence in the wake of cylinder has been fully developed, which causes the boundary layer of the body surface to transition to turbulence. When the cylinder is close to the twin-box girder, the wake of the cylinder cannot fully develop into turbulence to change the boundary layer, which explains why the Reynolds number effect still exists. In conclusion, when the spacing ratio *ε* > 2, the turbulence generated by the wake of the cylinder effectively eliminates the Reynolds number effects on the time-averaged lift force.

Figure 18 shows that the fluctuating lift force coefficient of the twin-box girder with cylinder interference is independent of the Reynolds number and spacing ratio *ε*. Moreover, the fluctuating lift force is much smaller than that of the undisturbed force at Re ≤ 9.0 × 103. This may be because the vertical pulsation component in the cylindrical wake is very small and the turbulence generated by the cylinder suppresses the vortex shedding of the twinbox girder to reduce the fluctuation of the lift force.

Compared with the undisturbed time-averaged moment force, Figure 19 shows that the mean moment force coefficient with cylinder interference is large, and is dependent on spacing ratio and diameter. With the decrease in diameter, the time-averaged moment force coefficient *Cm mean* is slightly enhanced. Meanwhile, with the increase in the spacing ratio, the absolute value of *Cm mean* decreases gradually. When the spacing ratio *ε =* 5, the *Cm mean* is smaller than the undisturbed mean moment force. This indicates that the fully developed turbulence wake generated by a large spacing ratio can effectively suppress the time-averaged moment force.

Figure 20 shows the fluctuating moment force coefficient *Cm rms* of the twin-box girder with cylinder interference. At Re ≤ 9.0 × 103, the interfered *Cm rms* is greater than the undisturbed *Cm rms*. When the spacing ratio *ε* ≤ 3, the *Cm rms* is insensitive to the diameter, and the *Cm rms* decreases slightly with the increase in the Reynolds number.

#### **4. Conclusions**

In the present work, the effects of two cutting-edge aerodynamic interference measures on the pressure distribution and aerodynamic force of a twin-box girder were investigated. We used the leading body-height grid and leading circular cylinder to increase the turbulence intensity of the incoming flow. The flows gained more energy before reaching the separation point, and the entrainment rate of the flow was enhanced, which not only affects the separation bubbles and reattachment points, but also changes the characteristics of the boundary layer. The conclusions are summarized as follows:

**Author Contributions:** Conceptualization, B.W. and S.L.; methodology, B.W. and S.L.; experiments, B.W., J.F. and S.L.; software, B.W.; validation, B.W.; formal analysis, B.W. and S.L.; investigation, B.W., G.X. and S.L.; resources, S.L.; data curation, B.W. and S.L.; writing—original draft preparation, B.W.; writing—review and editing, B.W. and S.L.; supervision, S.L.; project administration, S.L.; funding acquisition, S.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the NSFC under Grant No. 51878230,the Natural Science Foundation of Heilongjiang Province under Grant No. YQ2021E033, and the Postdoctoral Scientific Research Development Fund of Heilongjiang Province under Grant No. LBH-Q20021, along with support by the Heilongjiang Touyan Team and Fundamental Research Funds for the Central Universities.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **An Experimental Investigation of Passive Jet Control Method on Bridge Tower Wake**

**Yewei Huang 1,2 and Wenli Chen 1,2,\***

**Abstract:** In this study, we employed a four-hole cobra probe to measure the wake characteristics of a rec-tangular bridge tower model in a wind tunnel. The scale of the model was 1:30, and the Reynolds number varied from 1.38 × 105 to 2.27 × 105 by changing the yaw angle. A measurement plane with 9 × 19 measurement points was horizontally set at the middle height behind the model. The wake characteristics of the test model without control, i.e., the baseline case, was first tested in the yaw angle range from 0◦ to 90◦; then, four kinds of passive jet control cases were tested to study their control effects on the bridge tower wake. To evaluate the wake characteristics, three main aspects, i.e., mean velocity, turbulence intensity, and velocity frequency, were investigated. The meas-urement results indicate that the passive jet control method can achieve an effect in suppressing the turbulence of the wake but can slightly modify the mean velocity distribution. The dominant frequency distribution region was eliminated when the yaw angle was small but slightly expanded at a large angle. The differences between cases show a trend that the larger the suction coefficient is, the better the control effects are.

**Keywords:** passive jet control; tower wake characteristics; cobra probe

#### **1. Introduction**

The cable system bridge is the most widely used large-span bridge system. With the increasing requirement of the span, a growing number of extra-large cable system bridges have come into service, and the wind resistance performance is the primary limitation to the span of the bridge. As the bridge span grows, the wind effect appears more complex, and unprecedented wind-induced phenomena could occur. The wake flow characteristics of bridge towers is an interesting research branch that could lead to the oscillation of the cable behind the tower. For cable-stayed bridges, some parts are under the influence of the tower wake, particularly the long stay cable. In suspender bridges, the suspenders near the tower are wholly in the influence region of the tower wake. The oscillation motion of these suspenders is much more complex due to the unstable wake flow.

The characteristics of cylinder wake have been sufficiently investigated. Knisely [1] reviewed and further researched the Strouhal numbers of the wake of rectangular cylinders; Norberg [2] also investigated the wake frequencies by using a hot wire. Shimada and Ishihara [3] used a two-layer k-ε model to numerically study the aerodynamic characteristics of an infinite-length rectangular cylinder and gained a great agreement in the distribution of mean pressure, but the fluctuating pressure distribution was underestimated compared with the real model. Assi et al. [4] investigated the wake-induced vibration (WIV) of tandem circular cylinders and considered the vibration as a consequence of the interaction between the wake and the structure, which differs from vortex-induced vibration (VIV). The experiment confirmed the wake displacement theory proposed by Zdravkovich [5]. Carmo et al. [6] performed a numerical study on the phenomenon of WIV to a cylinder at a low Reynolds

**Citation:** Huang, Y.; Chen, W. An Experimental Investigation of Passive Jet Control Method on Bridge Tower Wake. *Appl. Sci.* **2022**, *12*, 4691. https://doi.org/10.3390/app12094691

Academic Editor: Giovanni Bernardini

Received: 2 March 2022 Accepted: 1 May 2022 Published: 6 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

number (*Re*). Bearman [7] reviewed the research of circular cylinder wakes, regarding wakeinduced vibration as a form of VIV. The downstream cylinder is under the collective effect of vortex shedding from both upstream structures and itself. Wang et al. [8] performed a two-dimensional numerical simulation of two tandem cylinders at low *Re* flow and found that the WIV phenomenon comes to aggravation when the size of the cylinder upstream is larger than that downstream. Li et al. [9] conducted a wind tunnel test and a computational fluid dynamics (CFD) simulation to investigate the WIV of the suspenders in wake of a bridge tower. The experiment found that the suspender in the wake of the tower would carry out an elliptical orbit motion at a specific wind speed when the frequency of the wake fits the natural frequency of the suspenders. In another part, a tower with four corners cut decreases the RMS aerodynamic coefficient and facilitates the suppression of the WIV of suspenders.

On the other hand, many studies have been conducted to investigate the suppression of VIV of cylinders. A technology based on suction and jet flow was confirmed to be effective. Amitay et al. [10,11] modified the aerodynamic characteristics of cylinders by using fluidic actuators. The closed recirculating flow regimes generated by jet flow could be equivalent to a virtual surface which can modify the characteristics of aerodynamics. Crook et al. [12] used synthetic jets to delay the separation point on cylinders, with the jet placed near the separation line. Many relevant studies were performed in the 21st century [13–15] that also indicated that the turbulence of wake also can be decreased, and even the vortex street could be eliminated. Hao Feng et al. [16,17] investigated the mode of wake by proper orthogonal decomposition, as the synthetic jet is applied to control the vortex shedding. The particle image velocimetry (PIV) results show that a synthetic jet could generate a vortex pair further downstream, and the proper orthogonal decomposition (POD) analysis indicates that this vortex pair, together with the vorticity shear layers near the cylinder, forms the vortex structures further downstream. Chen et al. [18] performed an experimental investigation on the suction-based flow control method, investigated the pressure distribution on the surface of a circular cylinder with suction slits, and found that the steady suction method exhibits excellent control effectiveness and can distinctly suppress the VIV, the amplitudes of cylinder vibrations, fluctuating pressure coefficients and lift coefficients of the circular cylinder model. They also applied PIV for further study [19], and the results show that the separation points on the test model in the PIV measurement plane are pushed further downstream, and the wake closure length also varies significantly along the span-wise direction of the test model for the cases with the suction flow control. Chen et al. [20,21] first proposed a passive jet apparatus, performed an experimental and numerical study to investigate a passive control method that has passive windward suction and a leeward jet over a circular cylinder, and found that the mean drag and the dynamic wind loads acting on the cylinder are suppressed. This apparatus could manipulate the shedding process of the wake vortices, and the jets could modify the wake stability by displacing the region of absolute instability further downstream. This method was also investigated by Zhang et al. [22] to apply to a bridge girder to control the flow separation. Chen et al. [23] further developed the passive jet method on a rectangular cylinder, taking a tower column as a prototype. The passive jet also has adequate control effects on the aerodynamic forces acting on the column.

The WIV of suspenders was found to be caused by the vortices from the bridge tower shedding into the tower wake [9,24]. Therefore, in this paper, we aim to further apply this passive jet apparatus to the bridge tower to modify the characteristics of the wake. We conducted an experiment to measure the time histories of wind speed at each measurement point. The experiment setup is introduced in Section 2, the results of tower wake characteristics are shown in Section 3, and finally, a discussion and conclusions are presented in Section 4.

#### **2. Experimental Setup**

The experiment was conducted in a closed-circuit boundary layer wind tunnel affiliated with the Joint Laboratory of Wind Tunnel and Wave Flume (Joint Lab WTWF), Harbin Institute of Technology. The wind tunnel contains two test sections. The size (width × height × length) of the larger one is 6.0 m × 3.6 m × 50 m, and the smaller one is 4.0 m × 3.0 m × 25 m. Our experiment was performed in the small section, where the wind speed can be adjusted continuously from 3 m/s to 30 m/s, the turbulence intensity of incoming wind speed is no more than 0.46%. The small section equips two turntable systems and a three-dimensional automated traverse measurement system.

#### *2.1. Bridge Tower Model*

The bridge tower model was designed using the Xihoumen Bridge as prototype, which is a steel box girder suspension bridge with a main span of 1650 m, located in Zhejiang Province, China. The prototype bridge tower is an overhead gantry type, which can be simplified as two separate independent rectangle cross-section bridge towers when only the wake flow of tower is investigated. We investigate only one tower model, without considering the interaction of the wake of two towers. The cross-section of the model is a 283 mm (*B*) × 217 mm (*D*) rectangle, with a side length of 23 mm square corner cut at each corner, and the scale ratio is set to 1:30, as shown in Figure 1. The data presented in this paper correspond to dimensionless distances of *B* = 283 mm and *D* = 217 mm in cross-wind and along-wind directions. The height of the column model is 2.9 m, and thus, the remaining 0.1 m to the ceiling can spare space for the clamp device. The model is assembled with thick planks, which provide enough stiffness to limit deformation. The surface of the model is flat and glossy by polishing. A jack system, mounted between the tower model and ceiling, firmly immobilizes the model. The tower model always maintains rigidity and keeps static during the experiment. At the bottom of the tower model, a 20 mm thick wooden turn plate is arranged to adjust the yaw angle to the incoming flow when the jack system is inoperative.

**Figure 1.** Schematic diagram of test layout.

The yaw angle of model *α* is set from 0 degrees to 90 degrees, with an interval of 10◦. The blocking ratio of the model varies by the yaw angle, reaching a minimum of 5.4% when the yaw angle is 0◦ and a maximum of 8.9% when the yaw angle is 50◦.

#### *2.2. Passive Jet Rings*

The present study adopts four kinds of passive jet rings to take control of the tower wake flow, which refers to the design in Chen et al. [23]. The ring is an annular flow channel with slits made by a 1 mm thick plastic plate, as Figure 2 shows. The top view of the ring is a rounded rectangular box, of which the outer size is 323 mm × 257 mm, with a 43 mm radius rounded corner in each corner, and the inner size is 283 mm × 217 mm; thus, it can closely fit the tower model. The intermediate part of 18 mm thickness, considering the thickness of the plate, constitutes the hollow flow channel. The height of the passive jet ring is 100 mm. Several rectangular slits are distributed around the side walls of the ring. The major differences between these four categories of jet ring are the size and the distribution of these slits. The size of slits in the first and third categories is 80 mm × 30 mm, and that in the second and fourth categories is 60 mm × 30 mm. For the first category, there are nine and seven slits on each long and short side, respectively, and for the second category, there are seven and five slits. For the third and fourth categories, the slits of the first and second categories on the long side are sealed, seven and five slits are distributed on each short side, respectively, and no slits are arranged on the long side. The clear spacing between two slits on each side is 30 mm. All the slits adopt axisymmetric distribution on both two principal axes.

In this study, we calculate the dimensionless coefficient of suction (Chen et al., 2015):

$$\mathcal{C}\_{\text{suc}} = \frac{\mathcal{S}\_{\text{suc}}}{H l} = \frac{h d (n\_{\text{a}} \sin \alpha + n\_{\text{b}} \cos \alpha)}{H (a \sin \alpha + b \cos \alpha)}$$

where *S*suc is the summary area of the suction slits on one passive jet ring projection to the plane perpendicular to the direction of flow; *l=a*sin *α* + *b*cos *α* is the total projection length of the ring perpendicular to the incoming flow; *H* is the height of the passive jet ring; *h* is the length of the suction slit; *d* is the width of the suction slit; *a* and *b* are the length of the long and short side of the ring, respectively; *n*a and *n*b are the number of the slits distributed on each long and short side, respectively, and *α* is the yaw angle of the ring.

#### *2.3. Experiment Details and Measurement Point Arrangement*

In this experiment, a four-hole pressure probe (Series 100 Cobra Probe, Turbulent Flow Instrumentation Pty Ltd., Victoria, Australia) was used to measure the characteristics of the bridge tower wake flow. There are 172 measurement points in the zone of tower wake, distributed on a 2.4 m × 1.8 m horizontal plane 1.5 m above the ground (half of the wind tunnel section height), as shown in Figure 1. In this plane, nine rows perpendicular to the incoming flow distribute at the interval of 300 mm, and nineteen columns parallel to the incoming flow distribute at the interval of 100 mm. The middle measurement point in the first row is 200 mm away from the back of the tower model. Taking the same scale ratio as the tower model, the realistic locations of four suspenders behind the bridge tower correspond to the middle points in rows 3, 5, 7, and 9. The measurement points distribution

can be seen in Figure 3. At each measurement point, the sampling frequency is 625 Hz, and the sampling time is 20 s.

**Figure 3.** Measurement points arrangement sketch map.

The cobra probe is fixed on the bracket on the three-dimensional traverse measurement system, with its probe head straight towards the incoming flow. The traverse measurement system has a built-in coordinate system, by which we can precisely change the position of the probe. It also has a jack system to firmly fix the whole traverse measurement system when data are collected. The position of the probe is changed to the measurement points mentioned above one by one during the experiment. The incoming wind speed *U*0 is set to 10 m/s, corresponding to the range of Reynolds numbers (*Re* = *ρU*0*L*/*μ*, where *ρ* is the air density, *μ* is the dynamic viscosity coefficient of air, *L* is the projection length of the tower cross-section) from 1.38 × 105 to 2.27 × 105, which varies by the projection length (*L=B*sin *α* + *D*cos *α*) when the yaw angle *α* changes. As the yaw angle changes, the corresponding positions of the suspenders vary, revolving the same degrees around the bridge tower. The locations of the first four imaginary suspenders behind the tower are labeled A1 to A4. When the yaw angle continues to increase, the imaginary suspenders behind the other side tower come into the wake zone of the tower model, the locations of which are labeled B1 to B4. There are five test conditions in this study, recorded as Case 0–Case 4. Case 0 is the baseline condition, in which only a tower model is tested, without passive jet control. Case 1–Case 4, corresponding to the category of the passive jet ring from 1 to 4, respectively, are controlled conditions, in which fourteen passive jet rings are installed on the tower model, distributed at the interval of 100 mm. The ring in the middle is installed 1.5 m from the ground, the same height as the probe. The yaw angle changes from 0 to 90 in Case 0 to Case 2, and from 0 to 30 in Case 3 to Case 4, because few suction slits work at a large angle in Case 3 and 4. The test conditions are shown in Table 1, and the value of suction coefficients at different yaw angles in each case are listed in Table 2.

**Table 2.** Suction coefficients for each case in the experiment.

#### **3. Results and Discussions**

*3.1. Mean Velocity Distribution*

The cobra probe has a natural property called the acceptance ratio. The wind speed data are collected only when the local wind speed is in the range from 2 to 50 m/s; moreover, the direction of the speed vector is contained by the front cone with a 45-degree half-angle [25]. Data history would be placed as zero when the acceptance condition is not satisfied. Zero data history is rejected when the mean velocity distribution is calculated. Only u and v velocity components are taken into consideration due to the two-dimensional model, and the velocity is expressed in non-dimensional form as the local wind speed divides the incoming wind speed. The fourth-order spline interpolation is adopted for each measurement point.

Figure 4 shows the mean dimensionless speed (*U*/*U*0) distribution in each test condition, where the black arrows refer to the local wind velocity vector. The general distribution trend of wind velocity is almost the same, and it can be divided into three regions. A low-value zone is located near the position (0, 0), which can be regarded as the low-speed region. This region is created due to the bridge tower shielding effects and the recirculation zone. Due to the limitation of the cobra probe, the mean velocity is always downstream, which is not tally with the fact; thus, the low-speed region is defined as the location where the data acceptance ratio (the proportion of non-zero data) of the probe is less than 50%. A slow diffusion region is located downstream of the low-speed region, where the speed is approximately half of the incoming wind speed. The high-speed region is distributed on two sides of the Y-axis, where the mean speed is almost not affected by the bridge tower model.

**Figure 4.** *Cont.*

**Figure 4.** Dimensionless mean speed (*U*/*U*0) distribution of tower wake.

The non-dimensional area (*S*1/*BD*) of the low-speed region is shown in Figure 5, where *S*1 is the area where the data acceptance ratio is less than 50% in the tower wake. In Case 0, there is a notable minimum area at 10 degrees equal to 0.386, and the second minimum area is equal to 0.503 at 80 degrees. This indicates that the low-speed area of the rectangular column with corner cuts minimizes not at the position perpendicular to the incoming flow, but at the position of little degree rotation. In Case 1 and Case 2, the low-speed region extends downstream compared with Case 0, every 0–90 degrees; in Case 3 and Case 4, an obvious expanding low-speed region can be observed at 10 and 20 degrees. These phenomena can be explicated by the previous study (Chen et al., 2015), where the PIV results show that the wake vortex is pushed by jet flow and stretched, which would induce the expansion of the low-speed region.

**Figure 5.** Dimensionless area of low-speed region.

To show the distribution of the velocity on the *Y*-axis, Figure 6 adopts the velocity of measurement points in the final row (*X*/*B* = 8.48) in different cases. The trend of velocity is consistent in all cases. The distribution on the *Y*-axis shows the shape of a saddle, and the velocity increases first and then decreases as the yaw angle increases. It can be indicated that the passive jet ring significantly modifies the velocity in the near wake but slightly changes the velocity distribution further downstream.

**Figure 6.** Dimensionless speed distribution at *X*/*B* = 8.48.

#### *3.2. Turbulence Intensity Distribution*

The turbulence intensity distribution in the wake region is achieved. To the measurement points with zero data, the zeros are rejected in the calculation, the same as velocity. The processing method may lead to some deviation from the true value when the acceptance ratio is low.

The fourth-order spline interpolation is also used to gain the final distribution. Figure 7 shows the distribution of the turbulence intensity distribution. All pictures reveal a universal distribution mode in the tower wake. A core area of the turbulence is located at the back of the tower model, which possesses a high turbulence intensity attaining about 60%. Then, a diffusion zone of the turbulence appears downstream, where the turbulence generally decreases, and the width of the zone increases. Two main differences between uncontrolled and controlled cases must be noticed. First, the core zone of turbulence is enlarged in some controlled cases. The explanation for this phenomenon is that the low-speed region is enlarged in controlled cases as shown in the velocity part, which directly leads to this consequence. Second, there is a low turbulence region inside the core zone in some controlled cases, which may be due to the jet flow from the passive jet ring, as the jet flow would suppress the shedding vortex at the back of the tower model.

**Figure 7.** Turbulence intensity distribution of tower wake.

To quantitatively compare the turbulence intensity, Figure 8 shows the turbulence intensity distribution in the row at X/B = 8.48, where the turbulence is totally diffused, and the acceptance ratio is 100%; thus, the deviation can be averted. At 0 degrees, Case 3 shows the best control effect, the peak value reduces by 19.3%, and the width of the wake decreases evidently; other cases show more or fewer control effects. At 10 and 20 degrees, all cases control the turbulence effectively. Case 3 and Case 4 attain the best consequence, which reduces 32.5% and 30.2% at 10 degrees and 17.9% and 20.7% at 20 degrees, respectively. When the yaw angle comes to 30 degrees, Case 3 has no control effect but amplifies the peak value for 7.4%, while Case 1 gains the best effect for a 32.5% reduction. The suction coefficient and flow transmission loss can elucidate this variation trend rationally. When the yaw angle is small (0–20 degrees), the differences in the suction coefficients between Cases 1 and 2 and Cases 3 and 4 are small, but the flow would leak from the side slits in Cases 1 and 2, which induces the abating of the jet flow. Hence, Cases 3 and 4 perform better in small yaw angles. As the yaw angle continuously increases, Cases 1 and 2 obtain larger suction coefficients than Cases 3 and 4; on the other hand, because each side functions as a suction or jet wall, no conception of leakage would exist, so Cases 1 and 2 can perform better at a large yaw angle. This is the reason that Cases 3 and 4 are only conducted at 0–30 degrees. A similar phenomenon also appears in the study of the application to a circular cylinder (Chen et al., 2015), where five pairs of slits perform better in wind fluctuating load suppression than 13 pairs of slits. At 40–90 degrees, Cases 1 and 2 gain a control effect in most degrees, and generally, Case 1 performs better. Case 1 attains the control reduction of [34.1, 26.8, 27.7, 30.2, 0.8, 27.0] percent at 40–90 degrees and Case 2 attains the control reduction of [6.6, 15.9, 10.9, 19.7, −18.5, 18.0] percent. Notably, at 80 degrees, Cases 1 and 2 have no control effects—Case 2 even has a negative reduction. The turbulence of Case 0 has an obvious decrease at 80 degrees, which needs to be further explained. The general tendency indicates that Case 1 would have considerable control effects in different yaw angles because the suction coefficient of Case 1 is the largest. The larger the suction coefficient, the better control effects.

**Figure 8.** *Cont.*

**Figure 8.** Turbulence intensity distribution at *X*/*B* = 8.48.

#### *3.3. Dominant Frequency Distribution*

Here, we processed the wind speed time history data to gain the power spectrum density. In this part, if there is zero data history, this is also taken into consideration, which means the total length of time history data is the same. We gathered the dominant frequencies of the spectrum in all measurement points, then transformed them into dimensionless frequencies (Strouhal number, *St* = *f* s*L*/*U*0, where *f* s is the shedding frequency, i.e., the dominant frequency in the wake). The distribution of dimensionless frequencies is obtained, as shown in Figure 9. To eliminate interference from dominant frequencies with small peak values, and for the comparison between cases, a threshold value of power spectral density (PSD) was set. Only the frequencies whose peak value was over the threshold value were counted, and the remaining zones were filled with zeros. If no PSD threshold was set, there is always a dominant frequency anywhere in the wake field, no matter how weak the energy level corresponding to that frequency. The threshold value of PSD in Figure 9 is set as 0.5. In this figure, most St numbers of wake in each case are uniform, but some

divergences are mostly located in the middle of wake, especially in Case 0. Comparing the cases, we can find that divergences decline in Cases 1 and 2; especially in Case 1, the St numbers are almost uniform, and the distribution appears clean.

**Figure 9.** Dimensionless dominant frequency (*fD*/*U*0) distribution of tower wake.

In Cases 3 and 4, the distribution region shrinks to 0 and 10 degrees. For a better comparison of the area change in the distribution region, Figure 10 shows the dimensionless area (*S*2/*BD*) of the region in each case. The threshold value of PSD is also set as 0.5. Two trends in different sections of the yaw angle appear. When the yaw angle is small (0–10 degrees), Cases 1, 3, and 4 all reveal a great control effect. Cases 1–4 attain [36.0, −47.2, 59.1, 50.6] and [63.3, 26.9, 78.0, 61.9] percent reduction of area in 0 and 10 degrees, respectively. However, when the yaw angle increases, the area expands, where Case 1 shows the least expansion, except at 80 degrees, with no more than 4% of the area expanding. As in the turbulence intensity distribution, a drop of

the area appears at 80 degrees in Case 0, which causes the amplification of the negative control effects of all cases at 80 degrees. Case 2 has the same trend of the area variation by degrees as Case 1, but the effect is worse—29% amplification of area is the maximum at 30 degrees, without considering 80 degrees. Cases 1 and 3 have larger suction coefficients, so these cases perform better than others. If only a small yaw angle is taken into consideration, Case 3 would be recommended, because the position of the suspenders would come out of the wake region at large yaw angles. For engineering applications, Case 1 would be better for various yaw angles.

**Figure 10.** Area of PSD > 0.5 dominant frequency distribution region.

#### **4. Conclusions**

In the present study, the characteristics of the wake of a bridge tower were tested with a four-hole cobra probe using the traverse system moving to each measurement point one by one. The cross-section of the tower model is rectangular with corner cuts, and the yaw angle was thus set from 0 to 90 degrees. Overall, the wake region is a diffusion zone, the width increases downstream, and a low-speed region is located at the near back of the tower. As the yaw angle increases, the area of the low-speed region first increases and then decreases, and the maximum appears at 50 degrees. The performance of the cobra probe is evidently restrained in the low-speed region because the acceptance ratio would be extremely low, and deviation is not avoidable. However, the ratio would rapidly reach 100% out of that region. Four controlled cases were set for comparison of the effects of modification in the wake of the model. Two of them were tested at 0–90 degrees, with suction and jet slits distributed on four sides, and the rest were designed for small yaw angles, with slits only distributed on the windward and leeward sides. The suction coefficients of Cases 1 and 2 vary in a small region, from [0.65, 0.35] to [0.67, 0.39], respectively, when the yaw angle changes. The suction coefficients of Cases 3 and 4 sharply decrease when the degree increases. The experimental results indicate that passive jet flow would enlarge the low-speed region; however, velocity distribution further downstream has no significant change. For the turbulence of the wake, the passive jet control method shows the control effects in several cases, where the peak value is suppressed and the width of the high turbulence region is also limited. An imperative trend should be noted: the higher the suction coefficient is, the better the control effects are. From the frequency results, we can conclude that the passive jet method, in some cases, can conspicuously shrink the area of the dominant frequency distribution region with a small yaw angle. However, this conclusion cannot fit the situation with a large yaw angle. The influence of the flow transmission loss is also reflected in the difference between Cases 1 and 3 or Cases 2 and 4. The jet slits distributed at the profile would cause the leak of flow, thus weakening the control effects, but it can contribute to keeping the suction coefficient at a steady level when

the yaw angle grows, and thus, a control effect on the turbulence can be maintained. Whether the side jet slits are adopted must take actual demand into consideration.

Similar passive jet rings have been proven to have a good control effect on the aerodynamic force of the bridge tower itself and can effectively suppress the vortex-induced vibration of the bridge tower [23]. In this study, we mainly focused on the effects of the passive jet method on wake characteristics. The results show that the passive jet can change the velocity, turbulence, and main frequency characteristics in the wake region of the bridge tower. It has been proven that the passive jet rings for the bridge tower will effectively control both the around flow field and the wake flow region; however, the control capability on the WIV of the suspenders remains to be further investigated for practical engineering.

**Author Contributions:** Conceptualization, Y.H. and W.C.; methodology, Y.H. and W.C.; software, Y.H.; validation, W.C.; formal analysis, Y.H.; investigation, Y.H.; resources, W.C.; data curation, Y.H.; writing—original draft preparation, Y.H.; writing—review and editing, W.C.; visualization, Y.H.; supervision, W.C.; project administration, W.C.; funding acquisition, W.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China grant number U2106222, 52008140, and 51978222.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work was funded by the National Natural Science Foundation of China through Grants U2106222, 52008140, and 51978222.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Aerodynamic Characteristics of a Square Cylinder with Vertical-Axis Wind Turbines at Corners**

**Zhuoran Wang 1, Gang Hu 1,2,3,\*, Dongqin Zhang 1,\*, Bubryur Kim 4, Feng Xu 1 and Yiqing Xiao 1,2,3**

**Abstract:** A preliminary study is carried out to investigate the aerodynamic characteristics of a square cylinder with Savonius wind turbines and to explain the reason why this kind of structure can suppress wind-induced vibrations. A series of computational fluid dynamics simulations are performed for the square cylinders with stationary and rotating wind turbines at the cylinder corners. The turbine orientation and the turbine rotation speed are two key factors that affect aerodynamic characteristics of the cylinder for the stationary and rotating turbine cases, respectively. The numerical simulation results show that the presence of either the stationary or rotating wind turbines has a significant effect on wind forces acting on the square cylinder. For the stationary wind turbine cases, the mean drag and fluctuating lift coefficients decrease by 37.7% and 90.7%, respectively, when the turbine orientation angle is 45◦. For the rotating wind turbine cases, the mean drag and fluctuating lift coefficients decrease by 34.2% and 86.0%, respectively, when the rotation speed is 0.2 times of vortex shedding frequency. Wind turbines installed at the corners of the square cylinder not only enhance structural safety but also exploit wind energy simultaneously.

**Keywords:** square cylinder; wind turbines; aerodynamic characteristics; vortex shedding

#### **1. Introduction**

Techniques to reduce wind-induced vibration of high-rise buildings have practical significance for ensuring structure safety and occupant comfort. In general, there are two traditional strategies, mechanical strategies and aerodynamic modifications, to reduce wind-induced vibrations. The tuned mass damper (TMD) installed inside tall buildings is one of the most effective mechanical strategies. In order to achieve an ideal vibration control effect, the mass of the damper normally needs to reach about 0.5–3% of the building mass. As a result, the damper not only occupies a large amount of building space, but also requires a stronger structure system to support such a heavy device. Furthermore, aerodynamic modification strategies, including rounding corners [1–3] and recessing corners [4–6], have also been studied. Wind tunnel model testing results demonstrate that slotted corners and chamfered corners can cause significant reductions in both along-wind and cross-wind responses [7,8]. Drag forces on a cylinder can also be reduced by adding a splitter plate behind it [9]. Although these strategies are very effective, they sacrifice the valuable space of tall buildings for the sole purpose of reduction in wind loading and wind-induced vibrations.

To further achieve the purpose of reducing drag, many researchers use the moving surface boundary layer control method (MSBC) to inject momentum into the flow field to prevent the formation of boundary layers. In [10], experimental research was conducted

**Citation:** Wang, Z.; Hu, G.; Zhang, D.; Kim, B.; Xu, F.; Xiao, Y. Aerodynamic Characteristics of a Square Cylinder with Vertical-Axis Wind Turbines at Corners. *Appl. Sci.* **2022**, *12*, 3515. https://doi.org/10.3390/ app12073515

Academic Editor: Mohsen Soltani

Received: 21 February 2022 Accepted: 28 March 2022 Published: 30 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

on the aerodynamic characteristic of a square cylinder under the momentum injection method. The results show that in the presence of MSBC, the drag of the square cylinder is considerably reduced at all angles of attack. In research reported in [11], a rectangular prism with two rotating cylinders was studied in a wind tunnel, and the effects of different rotation speed and modes on the aerodynamic characteristic were analyzed. Previous scholars mainly relied on the rotation of the small cylinder to inject momentum into the flow field. We consider replacing the small cylinder with wind turbines to increase its possibility and value in engineering applications.

Compared to the sole purpose of aerodynamic modification strategies to suppress wind-induced vibrations of tall buildings described in previous studies, it is desirable to develop effective aerodynamic devices for not only suppressing the wind-induced vibrations but also harvesting wind energy. Extensive studies aimed at making wind energy more efficient have also been carried out. The authors of [12] added two small diameter cylindrical rods on both sides of the cylinder to significantly improve the efficiency of a vortex-induced vibration energy collection system. Wind turbines have been installed in various positions of tall buildings, such as rooftop and through opening and channel between two buildings. The authors of [13] used the scale model to evaluate the wind speed amplifications in the tunnels for wind-power generation through the installation of wind turbines. To obtain optimal environmental benefits, the optimizing design of wind turbines installed on the tall buildings is also discussed [14,15]. Because previous studies have only focused on the energy-saving benefits of installing wind turbines on structures, we consider that installation of the wind turbine at the corners of the square cylinder will not only enhance structural safety but also exploit wind energy simultaneously.

This is a preliminary study to investigate the aerodynamic characteristics of a square cylinder with Savonius wind turbines and attempts to explain the reason why this kind of structure can suppress wind-induced vibrations. Numerical models, including numerical method, computational configuration and mesh arrangement, grid size and time step configuration, and simulation cases are illustrated in Section 2. Mean pressure coefficients, force coefficient and vortex shedding, and flow pattern around the cylinders for both stationary and rotating wind turbines are shown in Section 3. Finally, conclusions are summarized in Section 4.

#### **2. Numerical Simulation Configurations**

#### *2.1. Numerical Method*

RANS [16–18] treats the vortices of different scales equally in the flow field and smoothens the details of the temporal and spatial variations in pulsation motion through average operation. The turbulence model we used is the SST *k*–*ω*.

The continuity equation and momentum equation for impressible flow are expressed as:

$$\nabla \cdot \mathbf{U} = 0 \tag{1}$$

$$\frac{\partial \mathcal{U}}{\partial t} + \nabla \cdot (\mathcal{U}\mathcal{U}) = -\frac{1}{\rho} \nabla P + \nabla \cdot (\nu \nabla \mathcal{U}) - \nabla \cdot \tau \tag{2}$$

where *U* and *P* are the averaged velocity and pressure, respectively; *ρ* shows air density and *ν* refers to eddy viscosity; *τ* is the Reynolds stress which can be expressed as *τij* = *νt*(*∂ui*/*∂xj* + *∂uj*/*∂xi*) − 2/3*ρkδij* and *δij* is the Kronecker delta.

According to [18,19], the control equations of SST *k*–*ω* are expressed as follows:

$$\frac{\partial \rho k}{\partial t} + \mu\_i \frac{\partial \rho k}{\partial x\_i} = \frac{\partial}{\partial x\_j} \left[ \Gamma\_k \frac{\partial k}{\partial x\_j} \right] + G\_k - \Upsilon\_k \tag{3}$$

$$\frac{\partial \rho \omega}{\partial t} + u\_l \frac{\partial \rho k \omega}{\partial x\_i} = \frac{\partial}{\partial x\_j} \left[ \Gamma\_\omega \frac{\partial \omega}{\partial x\_j} \right] + G\_\omega - Y\_\omega + D\_\omega \tag{4}$$

where *ρ* = 1.225 kg/m3; *k* is turbulent kinetic energy; *ω* is the dissipation rate of turbulent kinetic energy; *ui* is the velocity component and subscripts *i*, *j* = *x*, *y*, *z*; *Gw* is the term of turbulent kinetic energy generation caused by the velocity gradient; Γ*k* and Γ*ω* are the convection terms for *k* and *ω*, respectively; *Yk* and *Yw* are the effective diffusion terms of *k* and *w* caused by turbulence, respectively; and *Dω* is the cross-convection term 1 *∂k ∂ω*

$$(D\_{\omega} = 2(1 - F1)\rho \sigma\_{\omega 2} \frac{1}{\omega} \frac{\cdots}{\partial \mathbf{x}\_j} \frac{\cdots}{\partial \mathbf{x}\_j}).$$

The coefficient of eddy viscosity *νt* can be obtained by:

$$\nu\_t = a\_1 \frac{k}{\max(a\_1 \omega, \Omega F\_2)}\tag{5}$$

$$F\_2 = \tanh(\arg\_2^2), \arg\_2 = \max(2\frac{\sqrt{k}}{0.09\omega y'}, \frac{500\nu}{y^2\omega})\tag{6}$$

where Ω is the absolute value of vorticity, *F*2 is a function that is one for boundary-layer flows and zero for free shear layers, and *a*1 = 0.31 is the empirical coefficient.

The transition from the *k*–*ω* model near the wall to the *k*–*ε* model far from the wall is controlled by mixed function *F*1:

$$F\_1 = \tanh(\arg\_1^4) \tag{7}$$

$$\text{arg}\_{1} = \min(\max(2\frac{\sqrt{k}}{0.09\omega y'}, \frac{500\nu}{y^2\omega}); \frac{4\rho\sigma\_{\omega 2}k}{\text{CD}\_{k\omega}y^2})\tag{8}$$

$$\mathbb{C}D\_{k\omega} = \max(2\rho\sigma\_{\omega 2}\frac{1}{\omega}\frac{\partial k}{\partial x\_{\rangle}}\frac{\partial \omega}{\partial x\_{\rangle}}; 10^{-20})\tag{9}$$

where *y* is the distance to the next surface, *νt* is the eddy viscosity, and *CDkw* is the positive portion of the cross-diffusion term. Model coefficient *σω*2 is 0.856.

The open-source computational fluid dynamics package OpenFOAM is used for numerical simulations in this study. The finite volume method is adopted to solve the unsteady incompressible Navier–Stokes equation. The PIMPLE scheme is used to solve the pressure–velocity coupling. The second order implicit scheme is selected for time discretization. Gauss linear upwind scheme is used for spatial discretization. The time step is set as 1.0 × 10−6 s. In terms of time step selection, to assess how the size of the time steps may affect the results, three time steps of 2 × 10−5, 1.0 × 10−6, and 1.0 × 10−7 s are tested. At last, 1.0 × 10−6 s is selected to balance the simulation time and the result accuracy. For all simulations, more than 1000 nondimensional time steps, corresponding to approximately 50 vortex-shedding cycles, which is much larger than the number of vortex-shedding cycles adopted in the previous studies [20,21], are taken for assuring reliable results.

#### *2.2. Computational Configuration and Mesh Arrangement*

The 2D computational domain, 92 *W* long and 25 *W* wide, is illustrated in Figure 1a in which the width (*W*) of the cross section is 0.1 m and the cylinder centroid is defined as the origin of the coordinate. The lateral boundary is 12.5 *W* away from the cylinder centroid in terms of the previous studies [22]. In order to allow flow redevelopment behind the wake region, the distances between the cylinder centroid and the inlet and the outlet boundary are 30.5 *W* and 61.5 *W,* respectively.

**Figure 1.** Views of grids for a square cylinder with wind turbines at two leading corners: (**a**) Grid distribution in x–y plane; (**b**) Grid surrounding square cylinder; (**c**) Grid surrounding upper wind turbine; (**d**) Grid surrounding lower wind turbine.

As for the boundary conditions, a smooth flow condition with a wind speed of 9 m/s is used at the inlet (*u* = 9*m*/*s*, *∂p*/*∂n* = 0) and the Reynolds number is 6.0 × 104, which is consistent with that in [22]. At the outlet boundary, a zero static pressure is adopted (*p* = *p*0 = 1.013 × 105, *∂u*/*∂n* = 0). Two lateral boundaries are defined as a symmetry boundary condition (*v* = 0, *∂u*/*∂n* = 0), which is used to model zero-shear slip walls in viscous flows. An empty condition is enforced on the top and bottom surfaces of the computational domain, which is an ordinary method to solve the 2D problem in OpenFOAM. Two vertical-axis wind turbines are installed at two leading corners of the square cylinder and the nonslip wall boundary condition (*ui* = 0, *∂p*/*∂n* = 0) is applied to surfaces of the cylinder and the wind turbines, where *n* refers to the normal direction, *u* and *v* represent the velocity components in x and y directions, respectively.

#### *2.3. Grid Size and Time Step Configuration*

The hybrid grid scheme is used for discretization, as shown in Figure 1. Because of the geometric irregularity of the wind turbines, an unstructured grid is adopted around the turbines in the vicinity of the cylinder. Structural meshes are used for the whole computational domain except the region containing wind turbines. As shown in Figure 2, the diameter of wind turbines is 5/100 *W*, *d* is 7/300 *W*, and *e* is 1/240 *W*, which leads to a ratio (*e/d*) of 1/6. The height of the minimum grid (δ) is 0.004 mm, the maximum *y*+ values for all cases is less than 1. The rotation of the wind turbine is implemented via the dynamic grid model, and a pair of coupling interfaces is set between the rotating region and the stationary region.

**Figure 2.** Sketch of Savonius wind turbines.

In order to consider the influence of grid size on numerical simulations, five grids (i.e., G1, G2, G3, G4, and G5) of the plain square cylinder with corner turbines are tested. The drag, lift coefficients, and Strouhal number of the square cylinder obtained from the five grid schemes are compared in Table 1. The mean drag, lift coefficients, and the Strouhal number are defined as:

$$\mathcal{C}\_d = \frac{F\_d}{1/2\rho l I^2 W} \tag{10}$$

$$C\_l = \frac{F\_l}{1/2\rho Ul^2 W} \tag{11}$$

$$S\_l = \frac{f\_s W}{U} \tag{12}$$

where *Fd* and *Fl* are the drag and lift forces of the cylinder, respectively; *fs* is the vortex shedding frequency; *W* is the width of the cross section; and *U* is the oncoming flow velocity set at the inlet boundary condition. The specifications of the grids used are given in Table 1 and wind turbine rotation speed is 11.88 rev/s. The Strouhal number and *RMS* lift coefficients *CRMS l* of G4 and G5 are identical. Considering accuracy and calculation speed, G4 is selected as the final grid scheme.

**Table 1.** Grid schemes for grid independence tests.

In addition, the simulation result of a plain cylinder without corner wind turbines is a baseline case to compare to the cases with corner wind turbines. As shown in Figure 3, the mean pressure coefficients on the four faces have favorable agreement with the results in the previous studies [22,23]. Table 2 gives the comparisons of force coefficients between previous studies and present numerical simulations. The definition of δcd is the relative error of *Cd*. As shown in Table 2, the drag coefficient *Cd* for the plain cylinder is 1.99, which agrees well with that of 2.05 for a square cylinder reported in [24]. The relative error δcd is nearly 5% when compared with all the experimental results. The *RMS* lift coefficient for the plain cylinder is 1.29, which matches well with the *RMS* lift coefficients of 1.37 for a square cylinder reported in [25].

**Figure 3.** Comparisons of mean pressure distribution between previous studies and present numerical simulations.

**Table 2.** Comparisons of force coefficients between previous studies and present numerical simulations.

#### *2.4. Simulation Cases*

As shown in Figure 4, two types of simulation cases, stationary and rotating wind turbines, are considered. For the stationary wind turbines, their orientation related to the wind flow direction may have a significant effect on aerodynamic characteristics of the cylinder. Similarly, for the rotating wind turbines, it is anticipated that the rotation speed of the turbines is a key factor which affects the aerodynamic characteristics of the cylinder. In order to study the effects of these two factors, 4 different orientations of the stationary wind turbines and 9 different rotation speeds are studied, as shown in Tables 3 and 4. The orientation of the stationary turbine is denoted by *θ,* which represents the angle of rotation of the wind turbines. Two wind turbines rotate inward relative to the square cylinder. *Vt* is the rotation speed of the wind turbines and *fs* is the vortex shedding frequency of the square cylinder calculated according to Strouhal number of 0.132 for the square cylinder provided by [22]. The rotation speed ranges from 0 to 4 times that of *fs*, and a rotation speed of 0 refers to a stationary case. In this study, the turbine rotates independently from the wind speed. According to [30], the rotation speed of the Savonius wind turbine is about 19.5 (r/s), which is closest to the rotation speed of R6. Thus, R6 is considered as an actual case. In order to conduct a comparative study on the influence of different rotational speeds on the aerodynamic characteristics of the square cylinder, we created 8 other working conditions.

**Figure 4.** Simulation cases: (**a**) Stationary cases; (**b**) Grid surrounding square cylinder.

**Table 3.** Rotating test cases.

**Table 4.** Stationary test cases.

#### **3. Results and Discussions**

*3.1. Mean Pressure Coefficients*

3.1.1. Stationary Wind Turbines at Cylinder Corners

Comparisons of mean pressure coefficients (*Cp*) around the cylinder surface among the four cases with stationary wind turbines at the cylinder corners and the plain cylinder are provided in Figure 5.

On the windward face, *Cp* at the stagnation point is approximately 1.0 for all cases. *Cp* gradually decreases from the stagnation point to two corners and *Cp* on two side faces are symmetrically distributed as expected. Because of the flow separation, the value of *Cp* on two side faces is negative. Compared with the cases with turbines at corners, *Cp* on the side faces of the plain cylinder has the lowest value. It means that the existence of stationary wind turbines in the corners reduces the suction on the side faces. Furthermore, *Cp* on the side faces has the largest value when *θ* is ±45◦ (*S2*), which means that the suction on the side faces of *S2* is lower than the other cases.

**Figure 5.** Mean pressure coefficients around the cylinder for the four stationary cases and plain cylinder.

On the leeward face, the suctions for cases with stationary wind turbines are significantly lower than that of the plain cylinder. Similar to the observations on the side face, the stationary wind turbines weaken the suction on the leeward face and it becomes more obvious when *θ* is ±45◦, i.e., *S2.*

#### 3.1.2. Rotating Wind Turbines at Cylinder Corners

Mean pressure coefficients (*Cp*) around the cylinder surfaces of nine rotating cases and the plain cylinder are compared in Figure 6. Figure 6a presents the cases when *Vt*/*fs* is less than 1 while Figure 6b shows the cases when *Vt/fv* is equal to or larger than 1.

**Figure 6.** Mean pressure coefficients around the cylinder for the nine rotating cases and the plain cylinder: (**a**) *Vt*/*fs* < 1; (**b**) *Vt*/*fs* ≥ 1.

On the two side faces, *Cp* of the plain cylinder has the lowest value among the nine cases. Therefore, the existence of rotating wind turbines in the corner reduces the suction on the two side faces. As shown in Figure 6a, the suction on the two sides increases as the rotation speed increases when the rotation speed of wind turbines is less than *fs*. As shown in Figure 6b, the suction on the two sides has an obvious reduction when the wind turbine rotation speed increases from *f* s (*R5*) to 2*f* s (*R6*). Lastly, the suction will not be lower than *R6* as the rotation speed further increases.

On the leeward face, the suction of cases with the rotating wind turbines is lower than the plain cylinder. Similar to the observations on the side faces, it is believed that the rotating wind turbine is able to weaken the suction on the leeward face. It is also noted that the suction on the leeward face has the lowest value, as shown in Figure 6b, when the rotating speed is 2fs (*R6*).

*3.2. Force Coefficients and Vortex Shedding Characteristics*

3.2.1. Stationary Wind Turbines at Cylinder Corners

*Cd*, *CRMS l* , and *St* of four stationary cases and the plain cylinder are listed in Table 5. The percentage of the values relative to the plain cylinder is shown in Figure 7.

**Table 5.** Force coefficient for four stationary cases and plain cylinder.

**Figure 7.** Comparisons of *RMS* lift coefficients, mean drag coefficients, and Strouhal number for the four stationary cases and plain cylinder.

As shown in Table 5 and Figure 7, four cases with stationary wind turbines have lower *Cd* than the plain cylinder since the suction on the two sides and leeward face is lower than the plain cylinder, as shown in Figure 5.

As shown in Figure 7, *Cd* is reduced 37.7% compared to the plain cylinder when *θ* is 45◦ (*S2*), which is the most significant for all cases. The lowest suction on the leeward face shown in Figure 5 contributes to the lowest *Cd* of *S2*. The variation in *CRMS l* is similar to that of *Cd* and the *RMS* lift coefficients has a 90.7% reduction when *θ* is 45◦ (*S2*). That is to say, the fluctuation of *Cl* was suppressed effectively. There is still a drop of 66.7% of *CRMS l* compared with the plain cylinder when *θ* is 135◦ (*S4*). Figure 8a,b presents the comparison of time history of *Cl* for *S2* and *S4*. The red curve represents the time history of *Cl* of the plain cylinder and the black curves represent the time history of the typical cases. It can be found that the fluctuation of *Cl* for *S2* has a greater inhibition than *S4.*

**Figure 8.** Comparison of time history of *Cl*: (**a**) *S2*; (**b**) *S4*; (**c**) *R1*; (**d**) *R6*; (**e**) *R8*.

Figure 9 illustrates the power spectrum of the lift force for the four stationary cases and plain cylinder; the vortex shedding energy is represented by the peak value of the spectrum. The frequency corresponding to the maximum value of the power spectrum is the vortex shedding frequency, and then *St* can be calculated by *fsW*/*U* (*W* is the width of the cylinder and *U* is the wind speed). As shown in Figure 9, the spectrum of the lift coefficient of the plain cylinder fits well with that in [31]. Four stationary cases have a much larger Strouhal number, approximately 20% larger than that of the plain cylinder. The vortex shedding energy presented by the peak value of the spectrum of four stationary cases is evidently smaller than that of the plain cylinder. Even though the stationary wind turbines installed at two leading corners increase the frequency of vortex shedding of the square cylinder, they can effectively reduce the energy of vortex shedding. Furthermore, the vortex shedding energy of the cylinder with corner wind turbines exhibits the lowest value when *θ* is 45◦ (*S2*).

**Figure 9.** Power spectra of lift forces for the four stationary cases and plain cylinder.

#### 3.2.2. Rotating Wind Turbines at Cylinder Corners

*Cd*, *CRMS l* , and *St* of the nine rotating cases and the plain cylinder are listed in Table 6. Ratios of these values to that for the plain cylinder are shown in Figure 10. All cases of the square cylinder with rotating corner wind turbines have lower mean drag coefficients (*Cd*) than the plain cylinder. Compared with the plain cylinder, when the rotation speed is 0.2*f* s (*R1*), *Cd* of the cylinder with wind turbines dramatically decreases by about 34.2%, and then the drop rate gradually decreases as the rotating speed increases before reaching *f* s. When the rotation speed increases to 2*f* s (*R6*), the drop rate of *Cd* again increases to 34.2%. Similarly, all cases of the square cylinder with rotating corner wind turbines have lower *CRMS l* than the plain cylinder. As rotation speed increases, the variation in the drop rate of *CRMS l* is similar to that of *Cd*, but the reduction in *CRMS l* is more significant than *Cd*. Compared with the plain cylinder, *CRMS l* of the cylinder with corner wind turbines dramatically decreases by about 86.0% as the wind turbine rotation speed is 0.2*f* s(*R1*) and the drop rate becomes 71.3% when the rotating speed increases to 2*f* s (*R6*). When the rotation speed is 4*f* s(*R8*), both *Cd* and *CRMS l* have the least reduction. It is worth noting that even though *Cd* and *CRMS l* have a significant reduction when the wind turbine has a rotating speed of 0.2*fs* and 2*fs*, it is still smaller than the reduction in the case with the stationary wind turbines (*R0*). Figure 8c–e present the comparison of time history of *Cl* of *R1*, *R6*, and *R8*. It can be found that the fluctuation of *Cl* of *R1* is smallest, which corresponds to the largest drop rate of *CRMS l* of *R1.* Figure 8d,e shows that the fluctuation of *Cl* of *R6* is also effectively suppressed, but the effect of fluctuation inhibition of *R8* is not obvious.

**Figure 10.** Comparisons of *RMS* lift coefficients, mean drag coefficients, and Strouhal number for the nine rotating cases and the plain cylinder.

Figure 11 presents the power spectrum of the lift force for the nine rotating cases and plain cylinder. Figure 11a refers to the cases with *Vt*/*fs* less than 1 and Figure 11b shows the cases with *Vt/fs* equal to or larger than 1. For Strouhal number (*St*) shown in Figure 11, all of the rotating cases have larger Strouhal numbers than the plain cylinder, except for *R6*. The vortex shedding energy represented by the peak value of the spectrum of the nine rotating cases is evidently smaller than that of the plain cylinder. Moreover, the vortex shedding energy of the cylinder with corner wind turbines has the lowest value when the wind turbine rotation speed is 2*f* s (*R6*).

**Figure 11.** Power spectra of lift forces for the nine rotating cases and plain cylinder: (**a**) *Vt*/*fs* < 1; (**b**) *Vt*/*fs* ≥ 1.

#### *3.3. Flow Pattern around the Cylinders*

The considerable effects of stationary and rotating wind turbines on the aerodynamic characteristics of a square cylinder are explained above. It can be found that there is a significant difference between aerodynamic characteristics of a plain square cylinder and a square cylinder with stationary or rotating wind turbines at the corners. To clarify the related mechanism, the flow patterns are discussed in this section.

#### 3.3.1. Stationary Wind Turbines at Cylinder Corners

Time−averaged pressure coefficient distribution around the plain cylinder and the cylinder of the two stationary cases (*S2* and *S4*) are shown in Figure 12a–c. As for the stationary cases, *S2* has the best effect of drag reduction while *S4* has the worst result, thus the flow pattern of these two cases is discussed herein. As for *S2*, the pressure distribution on the windward face is consistent with the plain cylinder, but the negative pressures coefficient of 0.68 on the leeward face is much smaller than 1.4 of the plain cylinder, which leads to a significant decrease in *Cd*, as presented in Figure 7. The time−averaged pressure coefficient on the leeward face of *S4* is about −0.84, leading to a 25.6% decrease in *Cd*.

**Figure 12.** Time-averaged pressure coefficient distribution: (**a**) plain cylinder; (**b**) *S2*; (**c**) *S4*; (**d**) *R1*; (**e**) *R6*; (**f**) *R8*.

Time-averaged streamlines around the plain cylinder and the cylinder of the two stationary cases (*S2* and *S4*) are shown in Figure 13a–c. For the plain cylinder, two recirculation zones beside the side faces and a pair of recirculation zones which are symmetrically distributed in the near wake can be found. The color of the streamline represents the x component of the velocity (*U*). It can be found that the wind flow has an increasing speed effect at the corner. The length of the recirculation zones is defined as the length between the plain cylinder center and the saddle point of the time-averaged streamtraces along the centerline of the test models [15,32,33]. It has been reported that the drag coefficient is inversely proportional to the length of the recirculation zones in studies on drag reduction

for square cylinders with cut-corners at the front edges [4] and the bluff body fitted with a splitter plate [34–37]. Less fluctuating lift coefficient can also be attributed to the elongated wake recirculation zone [29,38,39]. Apparently, the wake recirculation zone is dramatically elongated when *θ* is 45◦ (*S2*) compared with the plain cylinder and *S4*. Therefore, lower *Cd* and *CRMS l* for *S2* shown in Table 5 can be inferred. By contrast, the elongation of the length of the recirculation zone in *S4* is shorter than *S2,* leading to only a 25.6% decrease in *Cd* compared with the plain cylinder, as presented.

**Figure 13.** Time-averaged streamlines around the cylinder: (**a**) plain cylinder; (**b**) *S2*; (**c**) *S4*; (**d**) *R1*; (**e**) *R6*; (**f**) *R8*.

The suction on the side faces and leeward face depends on the curvature of the separated shear layer, the suction increases as the curvature of the separated shear increases [29]. Figure 14a–c shows the normalized time-averaged vorticity distributions of the plain cylinder and *S2* and *S4*. When *θ* is 45◦ (*S2*), the contour of normalized time-averaged vorticity distributions becomes thinner than that of the plain cylinder, which means the separated shear layer has smaller curvature, resulting in the lower suction that is shown in Figure 5. However, the contour of the normalized time-averaged vorticity distributions of the plain cylinder is stubbier, which means the curvature of the separated layer of the plain cylinder is larger. This may cause the increase in suction on the side faces and leeward faces of the plain cylinder, which may further lead to the increase of *Cd* in the structure.

**Figure 14.** Normalized time-averaged vorticity distributions: (**a**) plain cylinder; (**b**) *S2*; (**c**) *S4*; (**d**) *R1*; (**e**) *R6*; (**f**) *R8*.

Compared with the instantaneous spanwise vorticity distributions of the plain cylinder as shown in Figure 15a, the normalized instantaneous spanwise vorticity distributions in one period of *S2* in Figure 15b indicate that the stationary wind turbines push the separated shear layer away and the vortices form farther downstream compared with the plain cylinder, which weakens its interaction with the side faces [38,40–43]. The weaker interaction causes a lower spectrum peak in the lift force spectrum in Figure 8 and lower *CRMS l* , as shown in Figure 7. Overall, the stationary wind turbines installed in the corner can reduce the mean drag force and the fluctuation of lift remarkably.

**Figure 15.** *Cont*.

**Figure 15.** Normalized instantaneous vorticity distributions during one period: (**a**) plain cylinder; (**b**) *S2*; (**c**) *R1*; (**d**) *R6*.

#### 3.3.2. Rotating Wind Turbines at Cylinder Corner

Time-averaged pressure coefficient distributions of the three rotating cases are shown in Figure 12d–f. (*R1* and *R6* have excellent effects of drag reduction while *R8* has the worst result). The pressure distributions on the windward face of *R1* and *R6* are consistent with the plain cylinder, but *R1* and *R6* have lower negative pressure coefficients (0.7 and 0.65) on the leeward face, leading to about a 35% decrease in drag coefficient, as presented in Table 6.

Time-averaged streamlines around the cylinder of the three rotating cases and the plain cylinder are shown in Figure 13d–f. Apparently, compared with the plain cylinder, the wake recirculation zones are dramatically elongated when the rotation speeds are 0.2*fs*(*R1*) and 2*fs*(*R6*). Thus, 86.0% and 71.3% decrease in CRMS l of *R1* and *R6* can be observed. On the contrary, the wake recirculation zone of *R8* is shorter than that of *R6*, which results in a worse effect on drag reduction.

Figure 14d–f shows the normalized time-averaged vorticity distributions of *R1, R6,* and *R8*. When the rotating speed is 2*fs*(*R6*), the contour of normalized time-averaged vorticity distributions becomes more elongated than the plain cylinder, and the separated shear layer has smaller curvature due to the lower suction, as shown in Figure 6. As for *R1,* the situation is similar. However, the curvature of the separated shear layer of *R8* is distinctly larger, so the suction of *R8* is greater than *R1* and *R6.*

In Figure 15d, the symmetrical vortex shedding is observed in the normalized instantaneous vorticity distributions during one period of *R6*. The typical Karman-type vortex shedding is replaced by the symmetric shedding modes, which can lead to the large suppression of lift fluctuations shown in Figure 9. That is why the *RMS* lift coefficient (*CRMS l* ) of *R6* is significantly lower than that of the plain cylinder. As for *R1,* when the rotation speed is 0.2*f* s, the normalized instantaneous vorticity distribution during one period, shown in Figure 15c, is similar to *S2*, which means that *R1* has a similar mechanism of drag reduction as *S2.*

#### **4. Concluding Remarks**

This is a preliminary study to investigate the aerodynamic characteristics of a square cylinder with Savonius wind turbines and to explain the reason why this kind of structure can suppress wind-induced vibrations. The results indicate that the turbine orientation and the turbine rotation speed are the key factors that affect aerodynamic characteristics of the square cylinder. The numerical simulation results show that presence of either the stationary or the rotating wind turbines has a significant effect on wind forces acting on the square cylinder. The stationary wind turbines and rotating wind turbines installed at the leading corners can reduce both *CRMS l* and *Cd* of the plain cylinder. For the stationary wind turbine cases, when the turbine orientation angle is 45◦, the stationary wind turbines push the separated shear layer away and the vortices form farther away from the square cylinder compared with the plain cylinder, which weakens its interaction with the side faces, thus causing 37.7% and 90.7% reduction in the mean drag and fluctuating lift coefficients, respectively. For the rotating wind turbine cases, the mean drag and fluctuating lift coefficients are reduced by 34.2% and 86.0%, respectively, when the rotation speed is 0.2 times of vortex shedding frequency. When the rotation speed is twice the vortex shedding frequency, the typical Karman-type vortex shedding is replaced by the symmetric shedding modes, which leads to a large suppression of lift fluctuations (71.3%) and mean drag (34.2%).

This study investigated the effect of stationary or rotating Savonius wind turbines with two blades on the aerodynamic characteristics of a square cylinder. For the next stage, other Savonius wind turbines with more blades will be considered to investigate the effects of blades on the aerodynamic performance of a square cylinder.

**Author Contributions:** Conceptualization, G.H. and Y.X.; methodology, G.H., F.X. and D.Z.; software, Z.W.; validation, Z.W., G.H. and D.Z.; formal analysis, Z.W. and D.Z.; investigation, Z.W. and F.X.; resources, G.H. and Y.X.; data curation, Z.W. and B.K.; writing—original draft preparation, Z.W.; writing—review and editing, G.H., D.Z., B.K., F.X. and Y.X.; visualization, Z.W.; supervision, G.H. and Y.X.; project administration, G.H.; funding acquisition, G.H. and Y.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by National Key R&D Program of China (2021YFC3100702), National Natural Science Foundation of China (52108451, 52078175), Shenzhen Science and Technology Innovation Commission (GXWD20201230155427003-20200823230021001, GXWD20201230155427003-<PHONE_NUMBER>4428001), Shenzhen Key Laboratory Launching Project (ZDSYS20200810113601005), and Guangdong–Hong Kong–Macao Joint Laboratory for Data-Driven Fluid Mechanics and Engineering Applications (2020B1212030001).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Effect of the Extended Rigid Flapping Trailing Edge Fringe on an S833 Airfoil**

**Hongtao Yu and Zifeng Yang \***

Department of Mechanical and Materials Engineering, Wright State University, Dayton, OH 45435, USA<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** A 2D numerical simulation was conducted to investigate the effect of an extended rigid trailing edge fringe with a flapping motion on the S833 airfoil and its wake flow field, as an analogy of an owl's wing. This study aims to characterize the influence of the extended flapping fringe on the aerodynamic performance and the wake flow characteristics downstream of the airfoil. The length (*Le*) and flapping frequencies (*fe*) of the fringe are the key parameters that dominate the impact on the airfoil and the flow field, given that the oscillation angular amplitude is fixed at 5◦. The simulation results demonstrated that the airfoil with an extended fringe of 10% of the chord at a flapping frequency of *fe* = 110 Hz showed a substantial effect on the pressure distribution on the airfoil and the flow characteristics downstream of the airfoil. An irregular vortex street was predicted downstream, thus causing attenuations of the vorticities, and shorter streamwise gaps between each pair of vortices. The extended flapping fringe at a lower frequency than the natural shedding vortex frequency can effectively break the large vortex structure up into smaller scales, thus leading to an accelerated attenuation of vorticities in the wake.

**Keywords:** flapping fringe; CFD simulation; vortex attenuation; aerodynamics enhancement

#### **1. Introduction**

For decades, the quiet flight of owls has drawn interest from researchers in the aerodynamic and aeroacoustics fields. Owl-wing-inspired modifications on either leading-edge or trailing-edge designs have attracted plenty of research interest aiming to improve the aerodynamic performance and/or, in the meanwhile, suppress the aerodynamic noise from airfoils or turbine/fan blades. Multiple bio-inspired airfoil designs have been evaluated, including the trailing-edge or leading-edge serrations, trailing-edge flaplets, flapping wings, the compliant surface, and the flexible trailing-edge fringe, etc. [1–9]. It has been found that the soft trailing edge fringe can attenuate the primary noise source generated from the trailing edge through the modification of the coherent scattering mechanism. Consequently, aeroacoustic noise above 2 kHz can be effectively attenuated by the compliant surface covered by the down feathers through a bypass dissipation mechanism [10,11]. In the last decade, the use of an owl-wing-inspired trailing edge design as a passive control approach for an airfoil model to minimize drag and improve lift coefficients has gained attention [12–16].

The brush-like trailing edge has been utilized as the trailing edge extension to reduce both the narrowband bluntness noise and broadband turbulent boundary-layer trailingedge noise, as it contains the characteristics of both porosity and flexibility as the natural fringe [17–21]. A remarkable noise reduction was observed for an airfoil with a flexible trailing edge equipped to the suction side of the airfoil experimentally [21]. The attachment of flaplets on the trailing edge of an airfoil has clearly shown reductions in tonal noise, primarily resulting from the generation of small-scale vortices by the flaplets [22].

Recently, Murayama et al. (2021) carried out an experiment to study the effects of flexible flaps with various gaps on the aerodynamic characteristics of a fixed-wing [23].

**Citation:** Yu, H.; Yang, Z. Effect of the Extended Rigid Flapping Trailing Edge Fringe on an S833 Airfoil. *Appl. Sci.* **2022**, *12*, 444. https://doi.org/ 10.3390/app12010444

Academic Editor: Giangiacomo Minak

Received: 13 October 2021 Accepted: 29 December 2021 Published: 3 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The results indicated that the adjustable flaps enabled the suppression of the large-scale vortex shedding, and thus, reduced the fluctuation of aerodynamic forces. This finding is consistent with our experimental findings of an airfoil with soft feathers installed on the trailing edge [24]. Our previous experiment demonstrated that the real bird feather enabled the large-scale vortex to break up into small-scale turbulence. It is illustrated that the vibrating porous feather fringe produces a significant vertical velocity in the flow field around the trailing edge. As a result, it suppressed the leading-edge vortices approaching the flexible tail, and accelerated the shedding vortices' dissipation. Additionally, Yu et al. (2020) performed computational fluid dynamics (CFD) simulations on an airfoil with an extended trailing-edge fringe [25]. The effects of the fringe lengths and flapping frequencies on shedding vortices have been investigated. The flapping behavior with a certain length and frequency can attenuate the strength of the vortices, and improve the aerodynamics in terms of reducing the drag coefficients, and increasing the lift coefficients. The abovementioned studies aimed at improving the aerodynamics of an airfoil, and reducing the noise. The flapping trailing edge fringe can be applied to wind turbine blades and airplanes to improve noise reduction from the aerodynamic perspective.

In this study, to leverage the influence of the extended flapping trailing edge fringe on the aerodynamic performance of the airfoil, CFD simulations are performed for the S833 airfoil with an extended flapping trailing edge fringe. This numerical simulation model is different from our previous experimental study. First, the experimental model was equipped with real feathers at the trailing edge of the airfoil as a passive flow control without any controlled flapping motion. Due to the complexity of the real feather structure, the porosity, and the biological material properties, such numerical studies haven't been conducted yet. In this work, a simplified bare fringe without any porosity in a 2D fashion, together with specific flapping motions working as an active control, is adopted in the numerical modeling. In the simulation, three fringe lengths (*Le*) at four different flapping frequencies (*fe*) are evaluated thoroughly. The analysis focuses on the flappinginduced changes on the Q-criterion and velocity distributions within the flow field, and the aerodynamic pressure distributions over the airfoil.

#### **2. Materials and Methods**

#### *2.1. Computational Model*

The NREL's S833 airfoil has been extensively used to design wind turbine blades [26]. However, with the rapid development of wind farms, noise emission has become a significant issue. The ultimate goal of this study is to explore a novel airfoil design to reduce the noise generation associated with the shedding vortices from wind turbine blades. Thus, the S833 airfoil is selected to investigate the aerodynamic interactions between the vortex shedding in the wake and the flapping fringe installed at the trailing edge. Two-dimensional CFD simulations are performed for the airfoil model. The dimensions and boundary conditions applied in simulations are shown schematically in Figure 1.

The rectangular area of the outside domain is 20*C* by 10*C* (*C*: the chord length of the airfoil), with 5*C* from the inlet to the leading-edge of the airfoil. The flow boundary conditions applied in the simulation are adopted from the parameters in our previous experimental study [24]. The upstream velocity is set at 6 m/s, and the density and viscosity of the air are assumed as 1.23 kg/m3 and 1.84 × 10−5 Pa s, respectively, to match the air property at room temperature. This boundary condition corresponds to a Reynolds number (Re) of 40,000, which is the same as the Re number evaluated in the experiment. This Re number falls within the range of the Re number of owls' flight. The outlet boundary condition is defined as zero static pressure, and surfaces on the airfoil, as well as the top and bottom walls, are defined as non-slip wall conditions. The airfoil's angle of attack is set at 9 degrees, which is also consistent with the experimental study [24]. The fringe length (*Le*) ranges from 8% to 12% of the chord, with a step increment of 2%. The flapping frequency (*fe*) ranges from 80 Hz to 170 Hz, with a step increment of 30 Hz. The initial *fe* is determined as the shedding vortex frequency (*f*) of the baseline bare airfoil model, which is 140 Hz, obtained via the power spectral density analysis of the velocity in the downstream wake [25]. The dynamic motion of the trailing edge fringe can be expressed as *θ* = *A*sin(2*πft*), where *A* is the amplitude flapping angle (*A* = 5◦). The connection point of the flapping motion is fixed at the tip of the trailing edge.

**Figure 1.** The 2-dimensional computational domain and boundary conditions for the CFD simulation. *C*: chord length of the airfoil; *u*: velocity in the *x*-direction; *v*: velocity in the *y*-direction; *ρ*: density of the fluid; *μ*: viscosity of the fluid; *P*: pressure; *θ*: dynamic motion of the trailing edge fringe; *A*: amplitude of the flapping angle and set at 5 degrees; *f*: frequency of the flapping motion.

#### *2.2. CFD Simulations*

In this work, three fringe lengths at four different flapping frequencies of the extended fringe were evaluated. These variables resulted in 12 simulation cases, and a bare airfoil model was also simulated as the baseline case. Simulations were performed using the CFD software Cradle SC/Tetra 2021 (Dayton, OH, USA). The computational software solves the N–S equations by using the finite volume method with hybrid computational grids, wherein hexahedral grids were generated as prism layers around the airfoil to capture the boundary layer behavior, and tetrahedral grids were generated in the rest of the region. A local mesh refinement technique (grids become coarser along the downstream direction) was adopted to resolve the flow around the flapping fringes installed at the trailing edge to save computational effort.

The turbulence model LKE K-KL-ω was adopted to evaluate the turbulence region over the airfoil. This model is one of the Reynolds-Averaged Navier–Stokes (RANS)-based approaches to predict the transitional flow [25,27,28]. In the LKE model, the energy of the disturbances in the pre-transitional region of a boundary layer is expressed as "Laminar Kinetic Energy", whereas the turbulence energy is defined as k, and the transport equation of KL is solved with two equations of the fully turbulent model. The flow over the airfoil is governed by the Navier–Stokes equation for two-dimensional, viscous, incompressible flow. The continuity equation and the Reynolds-averaged N–S equations can be expressed as:

$$\frac{\partial(\rho \overline{w\_i})}{\partial(x\_i)} = 0 \tag{1}$$

$$\rho \frac{D(\rho \overline{u\_i})}{Dt} = \frac{\partial(\overline{P})}{\partial(\mathbf{x\_i})} + \frac{\partial}{\partial(\mathbf{x\_i})} (\mu \left(\frac{\partial(\overline{u\_i})}{\partial(\mathbf{x\_j})} + \frac{\partial(\overline{u\_i})}{\partial(\mathbf{x\_i})}\right) - \rho \overline{u\_i' u\_j'} \tag{2}$$

where *i* and *j* are coordinate indices; *u* is the velocity; *x* is the flow direction; *ρ* is the density of the fluid; *μ* is the dynamic viscosity of the fluid.

The dynamic motion of the fringe was fulfilled using the stretching mesh technique, as it enables the expansion and shrinking of the mesh elements to accommodate the movement of the solid boundary, which can avoid the pressure or velocity discontinuity/fluctuations at interfaces [25]. The CFD simulation, in terms of the turbulence model, grid sensitivity, and the dynamic moving mesh technique, has been validated statically by comparing

static lift and drag coefficients with experiments, and dynamically, by comparing pressure contours for an oscillatory plunging motion of an airfoil with the literature reported in our previous study [25].

The numerical simulations were performed on a computational cluster containing 3.4 GHz AMD Phenom II X4 965 Quad-cores and 8 GB RAMs. Each case was run with 32 processors. In the current study, a total physical time period of 1.5 s of the flow was simulated with a time step of 5 × 10−5 s. The mean computational time for each simulation case was about 23 h.

#### **3. Results and Discussion**

The impact of the extended flapping fringe on the airfoil's aerodynamics were studied primarily through the alterations in the distribution of the Q-criterion, velocities in the flow field, and pressure coefficients over the airfoil model. The variable Q-criterion quantifying the swirling vortex strength is defined as

$$Q = \frac{1}{2} \left[ \left( \frac{\partial u}{\partial x} \right)^2 + \left( \frac{\partial v}{\partial y} \right)^2 \right] - \frac{\partial u}{\partial y} \frac{\partial v}{\partial x} \tag{3}$$

where *u* is the instantaneous velocity in the *x*-direction, and *v* is the instantaneous velocity in the *y*-direction.

#### *3.1. Effect of the Flapping Motion on the Q-Criterion and Velocity Distributions*

The instantaneous (at t = 1 s) and time-averaged Q-criterion and velocity distributions with streamlines of the bare airfoil model (baseline) are presented in Figure 2. The quasisteady-state of the simulation was reached by checking the periodicity of the flow field. The vortices were shed evenly in the wake, with an identical magnitude of the Q-criterion for each pair, indicating typical Karman vortex street characteristics. The time-averaged Q-criterion shows very weak vortices downstream in the wake. It implies that the rotational direction of two vortices for each pair was opposite, and thus, after averaging, the opposite vorticities counteracts the magnitude of the Q-criterion. A separation bubble was predicted around the middle section on the upper surface of the airfoil, as shown in the instantaneous velocity distribution panel in Figure 2. In the time-averaged flow field, the detachment and reattachment point for the separation bubble can be clearly identified by the streamlines surrounding the circulation region. The separation bubble has a length of about 25% of the chord. In fact, the shedding vortices downstream of the airfoil were initiated by the separation bubble.

The Q-criterion (around the airfoil and in the wake) and velocity distributions for the airfoil model with the fringe length of 8% of the chord at four different *fe* are presented in Figures 3–5. The large-scale vortices initiated over the upper surface of the airfoil were broken up into smaller scale vortices compared to the shed vortices of the bare airfoil model. The flapping motion at the flapping frequency of 80 Hz creates an irregular pattern of vortices compared to the cases with higher frequencies (110–170 Hz), as shown in Figure 3. The time-averaged Q-criterion distributions demonstrate an accelerated decay of vortices by the extended fringe with a low flapping frequency (80 Hz), resulting from the unstructured small-scale vortices being shed irregularly. This observation might be induced by the velocity "offset effect" between induced vortices and the shed vortices from upstream. With the increased flapping frequency, the vortices decayed slowly. This might be attributed to the fact that the contribution from the flapping fringe to the overall swirling strength becomes more dominant with the increased flapping frequency.

**Figure 2.** The simulated instantaneous (time = 1 s) and time-averaged Q-criterion and velocity distributions with streamline of the bare airfoil model.

In the downstream wake, the airfoil model with *fe* = 110, 140, or 170 Hz shows evenly shed vortices, as shown in Figure 4. The increased flapping frequency seems to shorten the gap between each pair of vortices, which is directly correlated to the flapping frequency of the fringe. From the comparison of the four flapping frequencies, it is suggested that when the flapping frequency *fe* is 60 Hz (40%) lower than *f* (natural vortex shedding frequency of the bare airfoil model, 140 Hz), the regular vortex shedding is interfered with, and results in an irregular vortex distribution with non-uniform scales and shorter gaps between each pair. As can be observed, the flapping trailing edge fringe with *fe* = 80 and 100 Hz alters the coherent structure of the large-scale vortices in the wake, and at the same time, reduces the vorticity in the flow field (Figure 4).

However, it is difficult to differentiate the influence of the fringe tip position (fringe tip at the top or bottom positions) on the distribution of the Q-criterion. Thus, the instantaneous and time-averaged flow fields around the airfoil were plotted to illustrate the effects of the flapping fringe on the flow field, as shown in Figure 5. The flapping motion induces movement of the separation bubble towards the leading edge of the airfoil for all cases (Figure 5). The existence of the separation bubble caused fluctuated pressure forces on the upper surface, reflecting the unstable pressure distribution over the airfoil model, as shown in Figure 6. Higher flapping frequencies tend to increase the velocity over the top surface of the airfoil, and, as a consequence, it results in a decrease in the aerodynamic pressure. Subsequently, the pressure drop could lead to an increase in the lift coefficient, and a decrease in drag coefficient. It is indicated that the pressure acting on the top surface tends to be small when the flapping fringe tip is at the bottom position, in contrast to that when the fringe tip is at the top position.

**Figure 3.** Q-criterion distributions around the airfoil model with the fringe length of 8% chord at different flapping frequencies. Left column: instantaneous Q-criterion distributions at the fringe tip moved to the top; Middle column: instantaneous Q-criterion distributions at the fringe tip moved to the bottom; Right column: time-averaged Q-criterion distribution.

**Figure 4.** Instantaneous Q-criterion distributions in the wake of the airfoil model with the fringe length of 8% chord at different flapping frequencies at t = 1.5 s.

**Figure 5.** The velocity distribution around the airfoil model with the fringe length of 8% chord at different flapping frequencies. **Left column**: instantaneous velocity distributions when the fringe tip moved to the top; **Middle column**: instantaneous Q-criterion distributions when the fringe tip moved to the bottom; **Right column**: time-averaged Q-criterion distribution.

**Figure 6.** Pressure coefficient distribution over the airfoil model with the fringe length of 8% chord at four flapping frequencies.

The pressure coefficient distributions over the airfoil model with the fringe length of 8%C at four flapping frequencies are compared, as shown in Figure 6. The instantaneous pressure coefficient distribution shows a significant fluctuation corresponding to the shedding vortex over the upper surface as compared to the averaged pressure distribution over the baseline airfoil model. An obvious phase difference of the pressure oscillations can be observed between the case with the fringe tip at the bottom position and that with the fringe tip at the top position. The phase difference is more complicated in the region of the separation bubble compared to the region downstream of the separation bubble. By comparing the time-averaged pressure distributions, one can observe that higher flapping frequencies tend to move the separation bubble forward (towards the leading ledge). The distance from the leading edge to the take-off point of the separation bubble decreases from 33%C to 24%C as the flapping frequency increases from 80 Hz to 170 Hz. The separation bubble length is about 22%C, 25%C, 20%C, and 20%C for the case of *fe* = 80 Hz, *fe* = 110 Hz, *fe* = 140 Hz, and *fe* = 170 Hz, respectively. Overall, higher flapping frequencies tend to generate a lower time-averaged pressure coefficient over the upper surface, which would result in an increase in the lift coefficient.

Figures 7–10 present the instantaneous Q-criterion distribution over the airfoil, the Qcriterion distribution in the wake, the velocity, and the pressure coefficient distributions over the airfoil model, respectively, with the fringe length of 10% at different flapping frequencies. Figures 11–14 present the instantaneous Q-criterion distribution over the airfoil, the Qcriterion in the wake, the velocity, and the pressure coefficient distributions over the airfoil model, respectively, with the fringe length of 12% chord at various flapping frequencies. Similarly, the vortex shedding characteristics are similar to the airfoil model with the fringe length of 8% chord. The flapping motion with a relatively lower *fe* (<140 Hz) can break the large-scale vortex into small-scale vortices at the trailing edge, and, subsequently, alter the coherent structure of the vortex shedding. The irregularly shed vortices and the increased gap between each pair of vortices for cases with a lower *fe* accelerated the vortex decay compared to those with a higher *fe*. As the flapping frequency goes up to 170 Hz, the shedding frequency of the vortex is dominated by the flapping motion, which results in a regulated pattern of the shedding vortices. The time-averaged Q-criterion distributions show that the resultant time-averaged vortex strength has been elevated by the flapping fringe. It is indicated that the extended flapping fringe did not always facilitate the reduction of vorticity; conversely, it could promote the vortex domination in the wake, its interaction with the blunt trailing edge, and, possibly, the consequent noise emission from the trailing edge of the airfoil.

With the increased fringe length, the scale of the vortex became smaller, as shown in Figures 7 and 11. Correspondingly, the overall Q-criterion distribution over the top surface of the airfoil is decreased, as presented in Figures 8 and 12. These alterations in the flow field around the airfoil led to a decrease in the pressure coefficient over the upper surface of the airfoil model, as shown in Figures 10 and 14. Compared with the flapping frequency, the effect of the fringe length is less significant in terms of slight changes in both the vortex strength and the pressure distribution.

#### *3.2. Effect of the Extended Flapping Fringe on the Lift and Drag Coefficients*

The extended flapping trailing edge fringe not only reduced the swirling strength in the wake of the airfoil, but also decreased the drag coefficient, and increased the lift coefficient, surprisingly, as shown in Figure 15. It should be noted that in the integration calculation of the lift and drag, only the pressure distribution over the airfoil surface, excluding the trailing edge fringe, was considered. The drag coefficient essentially decreased with the increase in fringe length, and always decreased with the increase in flapping frequency, whereas the lift coefficient was essentially proportional to these two factors. The airfoil model with *Le* = 12%C and *fe* = 170 Hz generated the minimum drag coefficient, i.e., 22.5% lower than the baseline, and the maximum lift coefficient, i.e., 63.2% higher than the baseline, as can be observed in Figure 15. However, this combination would not alleviate

the swirling strength of vortices shedding in the wake. For this case, the flapping fringe serves more like a propeller tail, such as the tail of a fish, and thus, reduces the drag, and improves the lift. The drag coefficients of the models with *fe* = 80 Hz are slightly lower than that of the baseline model, whereas the lift coefficients are elevated dramatically. An *fe* = 140 Hz, which is the natural shedding frequency for the bare airfoil, did not generate drastically different results. The difference in the lift and drag generation due to the fringe length is the minimum for this case.

**Figure 7.** Q-criterion distributions around the airfoil model with the fringe length of 10% chord at different flapping frequencies.

**Figure 8.** Instantaneous Q-criterion distributions in the wake of the airfoil model with the fringe length of 10% chord at different flapping frequencies at t = 1.5 s.

**Figure 9.** The velocity distribution around the airfoil model with the fringe length of 10% chord at different flapping frequencies.

**Figure 10.** Pressure coefficient distribution over the airfoil model with the fringe length of 10% chord at different flapping frequencies.

**Figure 11.** Q-criterion distributions around the airfoil model with the fringe length of 12% chord at different flapping frequencies.

**Figure 12.** Instantaneous Q-criterion distributions in the wake of the airfoil model with the fringe length of 12% chord at different flapping frequencies at t = 1.5 s.

**Figure 13.** The velocity distribution around the airfoil model with the fringe length of 12% chord at different flapping frequencies.

**Figure 14.** Pressure coefficient distribution over the airfoil model with the fringe length of 12% chord at different flapping frequencies.

**Figure 15.** Time-averaged force coefficients for airfoil models with different fringe lengths at different flapping frequencies: (**a**) Drag coefficient; (**b**) Lift Coefficient.

In summary, the effects of the fringe length and flapping frequency on the vortex wake and aerodynamic performance are investigated numerically. It is found that a favorable *fe* should be determined according to the nature shedding frequency, *f*, of the bare airfoil model. The structure of the shedding vortices or their strength would be either maintained or elevated by the vibration of the fringe at a higher *fe* (>*f*). A higher flapping frequency (*fe*) tends to serve as a booster, rather than a disturber, to increase the local velocity at the trailing edge, resulting in the elevation of the swirling strength, and a shrinking of the dimension of vortices. However, most importantly, a lower *fe* was able to destruct the shedding vortices' coherent structure due to the disruption of the fringe on the original shedding pattern from the bare airfoil. Consequently, the swirling strength of vortices became weak, and the irregular distribution made it easy to dissipate.

#### *3.3. Limitation of the Current Study*

There are several limitations of this study. First, a RANS turbulence model was used to conduct the simulation. This would cause less accurate simulation results on the structure of shedding vortices. Due to the lack of comparable experimental data, direct validation of the simulation results could not be made, but, qualitatively, the predicted vortex reduction is similar to what we observed in our previous experiment [24]. In addition, to simplify the computational effort, a 2D model geometry was used, and thus, the span-wise porosity effect of feather fringes, as used in the experiment, could not be studied. To improve the understanding, 3D simulations on porous fringes attached at the trailing edge will be investigated in future work. To better mimic the feather fringe of an owl's wing, a soft material with a biological property should be adopted as well. The flapping motion will then be determined by the fluid-structure interaction as a passive flow control.

#### **4. Conclusions**

The airfoil S833 equipped with an extended flapping fringe is used to investigate the effects of two factors, i.e., the trailing fringe length and the flapping frequency, on the vortex shedding characteristics and the aerodynamic performance of the airfoil. The extended fringe with a flapping frequency significantly below the natural shedding frequency (~140 Hz), such as 80 Hz, can alter the coherence structure. Likewise, an irregular vortex structure in the wake can be generated, and thus, cause decrements in the swirling strength, and faster decay of the vortices. The flapping motion can decrease the pressure coefficient over the upper surface of the airfoil, leading to increases in the lift coefficient about 40% relative to the bare airfoil model. Out of the cases in the present study, the model with *Le* = 0.01 m (10% of the airfoil chord) at a flapping frequency of *fe* = 110 Hz outperforms

other cases, and shows an overall substantially better aerodynamic performance of the airfoil, as well as more favorable vortex characteristics downstream of the airfoil. In spite of the limitations of the study, we can observe the aerodynamic benefits of using an extended flapping fringe at the trailing edge. The pressure alteration around the airfoil would potentially reduce the noise generation for wind turbines especially, which will be studied in our future work as well.

**Author Contributions:** Conceptualization, H.Y. and Z.Y.; methodology, H.Y. and Z.Y.; software, H.Y. and Z.Y.; validation, H.Y. and Z.Y.; formal analysis, H.Y. and Z.Y.; investigation, H.Y. and Z.Y.; resources, H.Y. and Z.Y.; data curation, H.Y. and Z.Y.; writing—original draft preparation, H.Y. and Z.Y.; writing—review and editing, H.Y. and Z.Y.; visualization, H.Y. and Z.Y.; supervision, Z.Y.; project administration, Z.Y.; funding acquisition, Z.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The analysis data are available from the corresponding author on reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Study on Traveling Wave Wall Control Method for Suppressing Wake of Flow around a Circular Cylinder at Moderate Reynolds Number**

**Xin Liu 1, Weifeng Bai 2 and Feng Xu 1,\***

**Abstract:** In the present paper, the computational fluid dynamics (CFD) numerical simulation was utilized to investigate the effectiveness of the transverse traveling wave wall (TWW) method with the expectation of inhibiting the vortex shedding from a fixed circular cylinder. We mainly focused on the variations of four kinds of wave propagation directions, five different maximum wave amplitudes and ten different wave velocities for suppressing vortices shedding and aerodynamic forces. The aerodynamic coefficients and vortex structures under different propagation directions, wave amplitudes, wave numbers and wave velocities were investigated in detail. The results demonstrate that the alternate wake behind the cylinder can be effectively eliminated resorting to the "Downstream" propagating TWW. The mean drag coefficient is positively associated with wave velocity. Drag and lift coefficients remain relatively stable at different wave amplitudes. When the velocity ratio (wave velocity divided by incoming velocity) is 1.5, the lift coefficient fluctuation decreases to the minimum. In contrast, the optimal combination of control parameters under the present Reynolds number is concluded with "Downstream" propagating direction, maximum wave amplitude ratio of 0.02, and velocity ratio of 1.5.

**Keywords:** flow control; traveling wave wall; circular cylinder; numerical simulation; CFD

#### **1. Introduction**

Fluid flow around a cylinder is a universal phenomenon in engineering practice, i.e., heat exchanger tubes, marine cables, high-rises, and civil engineering structures. Vortex shedding behind bluff cylinders is the cause of vortex-induced vibration (VIV) which can result in structural damage under certain unfavorable conditions. In addition, there have been a variety of investigations about vortex-induced vibration from different perspectives in recent years [1,2]. Eliminating the alternating shedding vortices in the cylinder wake and suppressing the cylinder vibration are of important significance.

Fluid has a huge influence on bluff body structures, and there are many relative research [3–5]. Studies on flow control of bluff body have specified a new orientation for the flow control problems, and they have gained more and more concerns in recent years. Controlling methods can be mainly divided into two types: passive control and active control. Passive flow control could change the flow conditions to achieve the goal of flow control without consuming any external energy [6–10]. Another flow control method is active flow control, which requires external energy infused into the flow field. Moreover, active flow control could introduce the proper perturbation to change the inner flow mode. Several active flow control methods have been investigated to suppress vibration as well as improve aerodynamic performance, i.e., drag reduction, using wall

**Citation:** Liu, X.; Bai, W.; Xu, F. Study on Traveling Wave Wall Control Method for Suppressing Wake of Flow around a Circular Cylinder at Moderate Reynolds Number. *Appl. Sci.* **2022**, *12*, 3433. https://doi.org/10.3390/

Academic Editor: Artur Tyliszczak

app12073433

Received: 14 February 2022 Accepted: 23 March 2022 Published: 28 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

vibration [11,12], suction and injection [13–16], momentum injection [17–20], bionic control method of traveling wave wall, etc.

The TWW is one of the bionic based flow control methods and it has been developing rapidly in recent years. The boundary layer near the rigid wall comes into existence when the viscous fluid flows pass the wall. However, making the cylinder surface to be flexible or movable may help to weaken the boundary layer.

The "fluid roller bearing" (FRB) effects of the axisymmetric TWW was used by Yang and Wu [21] to sharp into a series of vortex rings. The main flow was separated from the near-wall flow accompanied by friction drag and pressure drag significantly reduced. Wu et al. [22] carried out a numerical simulation to find out the proper parameter (wave amplitude, wavelength and ratio of wave velocity to incoming velocity) for the FRB effects on infinite and finite two-dimensional TWW. Large flow separation was eliminated after applying TWW on airfoils, and the vortex shedding of streamline airfoils with a large wind attack angle was restrained to provide more lift. Wu et al. [23] proposed a transverse traveling wave (TTW) control method to manipulate the unsteady flow around a circular cylinder at Re = 500. A fluid FRB was introduced in this method, and the vortex shedding was eliminated because the global flow remained attached to the surface. Xu et al. [24] simulated the rear section of a cylinder with TWW, and the cylinder was elastically mounted with two degrees of freedom. The simulation took fixed wave amplitude, wave number and ratio of wave velocity to incoming velocity. In addition, the whole process was completely simulated starting with flow around the fixed cylinder, to the oscillating cylinder, and finally to the oscillating cylinder with TWW.

In the present paper, TWW was placed on the rear section of a fixed cylinder, and the CFD numerical simulation was adapted to investigate the influencing factors for suppressing the cylinder oscillating wake. The TWW was activated when the alternating shedding vortex behind the fixed circular could be stably observed in the present numerical simulation. There are several steps we need to comply with as per control variate technique. Firstly, fixed velocity ratio and maximum wave amplitude were used to measure the control effect of four kinds of TWW propagation direction for suppressing the oscillating wake. Then, the most effective propagation direction was selected, based on which we forward studying the controlling effectiveness of 5 different wave amplitudes, 4 different wave numbers and 10 different wave velocities. The control effect under various influencing factors was elaborated by comparing the characteristic values, aerodynamic force time histories and the wake flow patterns.

#### **2. Numerical Model and Validation**

#### *2.1. Governing Equations and TWW*

The governing equations of two-dimensional incompressible flow in a Cartesian coordinate system can be written as follows.

$$\frac{\partial \mu\_i}{\partial \mathbf{x}\_i} = 0 \tag{1}$$

$$\frac{\partial u\_i}{\partial t} + u\_j \frac{\partial u\_i}{\partial x\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial x\_i} + \frac{\mu}{\rho} \frac{\partial^2 u\_i}{\partial x\_j^2} \tag{2}$$

where *ui* is the velocity components, i.e., *u*1 is the velocity component at the inline flow direction and *u*2 is the velocity component at the cross-flow direction, *μ* denotes the kinematic viscosity coefficient, *ρ* is the fluid density, and *p* denotes the pressure in flow field.

As shown in Figure 1, The transverse traveling wave at the rear section of cylinder was established, i.e., the transverse traveling wave can propagate from both point B and C to point A at the same time.

**Figure 1.** Schematic of fixed circular cylinder with TWW.

The origin *O* is positioned at the center of the plane polar coordinate system (*r*, *θ*), and the wave equations of the TWW can be written as,

$$\begin{cases} \quad \text{x = \, \, r\cos\theta \\ \quad y = r\sin\theta \\ \quad r = \, r\_0 + A(l)\cos[k(l-ct)] \end{cases}, -\pi/2 \le \theta \le \pi/2 \tag{3}$$

$$A(l) = \begin{cases} \begin{array}{c} \hat{A}\frac{l}{\lambda}, 0 \le l \le \lambda\\ \hat{A}, \ \lambda \le l \le (N-1)\lambda \end{array} \\\ \hat{A}\frac{(N\lambda - 1)}{\lambda}, (N-1)\lambda \le l \le N\lambda \end{array} \tag{4}$$

where *ro* = *D*/2 is the radius of the standard circular cylinder without any traveling wave, *D* is the diameter of the cylinder, *A*(*l*) is the vibration amplitude of each point on the TWW, *k* = 2*π*/*λ* is the wave number, *l* is the arc length from point A on the cylinder rear edge to any point P on the TWW, *c* is the wave velocity, *t* is the time, *λ* is the wave length, *A*ˆ is the maximum wave amplitude, *N* is the wave number in 1/4 circle. In the first wave, the amplitude increases linearly from 0; In the end wave, the amplitude decreases linearly to 0, while other intermediate waves have the same amplitude as the complete waveform, as shown in Equation (4). The TWW connect smoothly with the non-moving cylinder surface in this configuration.

#### *2.2. Computational Domain and Boundary Conditions*

Figure 2 shows the grid distribution and the computational domain. The computational domain is a rectangle area with a length of 60*D* and a width of 40*D*, and the center of the cylinder is at the coordinate origin. The cylinder diameter is 0.12 m. The upstream inlet is 20*D* ahead of the coordinate origin, and the downstream outlet is 40*D* after the coordinate origin. The distances between either of the upside and downside and the coordinate origin are 20*D*. Unstructured grid is adopted to discretize the computational domain into four layers. The grids in the regions with large gradients of flow parameters are locally refined, i.e., the region near cylinder surface and wake region. The geometric model and meshing are performed by pre-processing software ICEM.

The boundary conditions are set up as follows, the inlet is set as the velocity-inlet with a uniform velocity *U*∞, the outlet is set as the pressure-outlet with the relative pressure of 0, the upper and lower sides are set as the symmetry boundary, and the cylinder surface is set as the no-slip wall.

The numerical calculations are carried out with a Reynold number of 4.1 × 104, and the turbulence flow is calculated by SST *k*-*ω* model. SIMPLE algorithm is adopted to calculate the coupling between the pressure and velocity fields. The pressure interpolation format is set as "Standard". The second-order upwind scheme is used for the momentum discretizing because of its accuracy and stability. During the solution process, the convergence residual standard of continuity equation is 1.0 × 10−6, and the convergence residual standard of momentum equation and turbulence parameters is 3.0 × 10−7.

**Figure 2.** Grid distribution and computational domain near the cylinder surface.

The CFD general software Fluent is used to calculate the flow around a fixed cylinder. Firstly, the governing Equations (1) and (2) are solved to obtain the flow characteristics around the cylinder, and then the aerodynamic coefficient of the cylinder is obtained. When the amplitude of the aerodynamic coefficient is stable (*t* = 1.8 s), TWW is activated. The wave Equation (3) of traveling wave is solved by User Defined Function (UDF) in Fluent, and the radial velocity at any point P in Figure 1 can be obtained according to Equation (3), and then the radial velocity is decomposed into the x direction velocity and the y direction velocity. The DEFINE\_GRID\_MOTION macro of UDF is used to assign the velocity for all grid points on the rear semi-cylindrical surface, and the motion of cylinder surface grids are realized by the dynamic mesh technique. When the grids are updated to the specified position, the flow field is calculated to converge with this shape as the boundary, and then the calculation of the next time step begins. This loop calculation ends until the TWW achieves a significant wake control effect.

#### *2.3. Validity Investigation*

The grid and time step independence verification were carried out because the mesh size and the time step determine the accuracy of CFD computational results. The comparisons of the global parameters for the grid and time step independences studies are shown in Tables 1 and 2. In the tables, *Nc* is the number of nodes for 1/4 circle, *Nmesh* is the total number of the mesh, and Δ*t* is the time step. The global parameters include the fluctuating lift coefficient *C*- *l* , the mean drag coefficient *Cd*, the fluctuating drag coefficient *C*- *d*, the Strouhal number *St*, and the minimum, mean and maximum values *y*+ *min*, *y*+ and *y*+ *max*, respectively. In order to guarantee the correctness of the turbulence simulation results and sufficient amount of mesh in the viscous sublayer, the *y*+ *max* should be less than 11.63 which is the demarcation point between the log-law region and the viscous sublayer recommended by Versteeg and Malalasekera [25].

**Table 1.** Comparisons of the global parameters for the grid independence study.

**Table 2.** Comparisons of the global parameters for the time step independence study.

Through comparison and analyzation of the global parameter results, *Nc* = 250 and Δ*t* = 2.5 × 10−4 s were adopted in the final results discussion. Current results accompanied with some previous simulation and experimental results [26–36] for the flow around a single cylinder are shown in Figure 3. The mean drag coefficient *Cd* versus Reynolds is shown in Figure 3a, and it can be observed that *Cd* is very close to previous results at similar Re range. The *St* in the previous 3D flow experiment is less than that in this paper at similar Re range. In the present 2D numerical simulation, the computational domain height is extraordinarily small, so the 3D flow in spanwise direction cannot be fully developed. Thus, the present *St* results are closer to the 2D flow results [28,34,35], as shown in Figure 3b. From the previous comparison, it can be concluded that the present numerical model and grids are reliable and have good reference value for further research of TWW flow control.

**Figure 3.** Comparisons of aerodynamic coefficient and the Strouhal number with previous results. (**a**) *Cd* versus Re. (**b**) *St* versus Re.

#### **3. Results and Discussion**

#### *3.1. Influence of Different Propagation Directions*

The control effectiveness of 4 types of TWW propagation should be considered first. The four propagation directions include the "Downstream", the "Upstream", the "Corotating-Clockwise" and the "Corotating-Counterclockwise", as shown in Figure 4. During calculation process, the effects of the TWW propagation directions are considered under fixed velocity ratio *c*/*U*∞ = 2.0, wave number *N* = 4, and maximum wave amplitude *A*ˆ/*D* = 0.02. This numerical simulation begins with the flow around a fixed cylinder (FR). The TWW is activated to control the flow field after the flow field is brought into stability with alternating shedding vortices. The aerodynamic coefficients and vortex shedding patterns before and after activating the TWW could serve as a good measurement for the control effect of TWW control method.

**Figure 4.** Four types of propagation direction of TWW. (**a**) Downstream. (**b**) Upstream. (**c**) Corotating-Clockwise. (**d**) Corotating-Counterclockwise.

The time histories of lift and drag coefficient for the TWW cylinder with the four different propagation directions are shown in Figure 5. When *t* < 1.8 s, TWW was not activated, the result of FR was obtained. When *t* ≥ 1.8 s, TWW was activated, the suppression effect of TWW on the vortex shedding and wake under the four different propagation directions can be obtained.

Figure 5 shows that the lift coefficient fluctuation reduces dramatically only when traveling wave propagates in the downstream direction. It indicates that the TWW propagating downstream has the effect of eliminating Vortex Street. In addition, the mean drag coefficient decreases at this propagation direction, while the fluctuating drag coefficient increases dramatically, as shown in Figure 5a. Figure 5b shows that the fluctuating lift coefficient and mean drag coefficient increase significantly when traveling wave propagates in the upstream direction. It illustrates that the traveling wave propagating upstream enhances the wake behind the cylinder. The drag coefficient fluctuation increases, and the mean lift coefficient varies significantly when traveling wave co-propagates along the cylinder rear surface. When traveling wave Co-rotates clockwise, the mean lift coefficient increases significantly to a positive value, as shown in Figure 5c. When traveling wave Co-rotates counterclockwise, the mean lift coefficient decreases significantly to a negative value, as shown in Figure 5d. It indicates that there is a negative pressure zone on the upper or lower surface with these two propagation directions, and the oscillating wake of traveling wave cylinder cannot be eliminated.

**Figure 5.** Lift and drag coefficient time history curves under various propagation directions. (**a**) Downstream. (**b**) Upstream. (**c**) Corotating-Clockwise. (**d**) Corotating-Counterclockwise.

The mean value reflects the variant characteristics within a period, and the fluctuation value (root-mean-square, RMS) reflects the discrete degree of the variant. Before activating the TWW, the mean lift coefficient *Cl* approximates 0 and *C*- *l* is 0.921, the mean drag coefficient *Cd* is 1.267 and *C*- *d* is 0.081. The time history curves of lift and drag coefficients for uncontrolled cylinder are taken for reference in the following parts to measure the influence of different control parameters on the control effect. Figure 6 shows the lift and drag coefficients characteristics with different traveling wave propagation directions, while the dashed curves inside the same figure demonstrate the results of the standard circular cylinder. As shown in Figure 6a, *Cl* approaches to 0 for "Downstream" and "Upstream" propagation direction, while the *Cl* increases dramatically to 3.917 and −4.057 for co-propagations. For the "Downstream" propagation, the *C*- *l* approaches to 0, which indicates that the cylinder oscillating wake is eliminated by TWW. For the "Upstream" propagation, *C*- *l* increases dramatically to 1.825. In addition, *C*- *l* for the co-propagations is close to that of the standard circular cylinder. These illustrate that the other three propagations are not able to suppress the cylinder wake vortex except "Downstream". Figure 6b shows that the fluctuating drag coefficients *C*- *d* vary between 0.172 and 0.763 in various propagation directions, and all of them are larger than that of the standard circular cylinder, which means that TWW leads to a larger amplitude fluctuating drag. For "Downstream" propagation, *Cd* decreases dramatically to −0.152, meaning the cylinder is subjected to reverse thrust due to TWW. For "Upstream" propagation, *Cd* increases dramatically to 2.809, i.e., the inline flow drag increases dramatically. The *Cd* of other two propagation directions equal to 0.905 and 0.902, respectively, which are slightly lower than the standard circular cylinder. Compared with the standard circular cylinder results, *C*- *l* and

*Cd* of "Downstream" propagation decreases by 99.48% and 112%, respectively. It illustrates the cylinder wake oscillation in FR stage can be eliminated by the TWW propagating downstream.

**Figure 6.** Characteristic values of lift and drag coefficient for different propagation directions. (**a**) *C*- *l* & *Cl*. (**b**) *C*- *d* & *Cd*.

The vorticity contours surrounding TWW cylinder with "Downstream" propagation are shown in Figure 7. In the interval *t* < 1.8 s, the traveling wave was not activated, and the vortex sheds from both side of the cylinder alternately, as shown in Figure 7a. Figure 7b shows that the TWW activated shortly has little effects, and the separated free shear layer still rolls into vortices. Form Figure 7c–g, as time advanced, the control effect of TWW became more obvious, and the vortex street appeared in uncontrolled cylinder was eliminated. The shear layer detaches from the cylinder surface, and then enters the wake region without rolling into vortices due to the "fluid roller bearing" (FRB) effect of TWW, and the vorticity is distributed steadily and symmetrically in the cylinder wake region, which means that the TWW cut off the energy source for generating vortex. When *t* = 2.97 s, the TWW eliminates oscillating wakes of the cylinder successfully.