The flapping vortex dynamics of two flexible plates submerged sidebyside in the wake of a square cylinder are investigated through a twoway fluid–structure interaction (FSI) simulation. The gap between the two plates can stabilize wakes, lengthen vortex formation, elongate vortices, suppress vortex shedding, and decrease hydrodynamic forces. The numerical results indicate that the two flexible plates can exhibit four distinct modes of coupled motion: outofphase flapping, inphase flapping, transition flapping, and decoupled flapping, depending on the gap spacing. Additionally, it is discovered that each of the four coupling modes has a unique pattern of vortex development. The findings of this study should proved valuable in the design of FSIbased piezoelectric energy harvesters utilizing cylinder–plate systems.
I. INTRODUCTION
Fluid–structure interaction (FSI) is a complex phenomenon that occurs when a flexible structure is submerged in a fluid flow and experiences deformation and vibration due to fluid forces acting on its surface (Khayyer , 2021). These deflections alter the flow field and result in intricate coupling dynamics that have attracted increased attention over the past few decades. This is due to the widespread use of FSI in various ocean engineering applications, including bionics (Luo , 2020), tidal turbines (Lothode , 2022), schooling fish (Devilliers , 2016), paper processing (Chen , 2019), and piezoelectric energy harvesting (Li , 2018). Despite its prevalence, the dynamics of FSI remain a subject of interest and some confusion.
FSI of flexible structures poses a complex and multidisciplinary challenge, involving the interplay between the elastic structure, fluid flow, and material properties (Bukac, 2016). To gain insight into the dynamics of this interaction, researchers often use flexible filaments (Zhang , 2020) or plates (Zhang , 2022) as simple yet effective models. Over the past few decades, extensive research has been conducted in this field, revealing that flexible plates may suddenly become unstable and exhibit increased flapping motion when fluid flow exceeds a critical speed (Brousseau , 2021; Jeanmonod and Olivier, 2017; Lee and Lee, 2012; 2013; Lee , 2012; Lim and Xiao, 2019; Shin , 2007; and Zhu , 2017). This instability is believed to result from an unstable equilibrium between fluctuating fluid pressure forces and the plate’s flexible nature. Additionally, the interconnected behavior of two cantilevered flexible plates that are aligned with each other has been analyzed. The flapping threshold depends on the distance between the two plates, and the system of two plates can oscillate in both inphase and outofphase modes.
Recently, the demand for power generation for lowpower electronic devices in various ocean engineering applications has increased, leading to a growing use of piezoelectric energy harvesters that leverage the FSI of flexible structures in crossflow. There have been a number of investigations of the flowinduced vibration (FIV) of a cantilevered flexible plate that is equipped with piezoelectric materials with the aim of designing a variety of FSIbased electronic devices (Eugeni , 2020; Goushcha , 2015). Inspired by the energy harvesting mechanism of piezoelectric materials through FIV, it is crucial to gain a comprehensive understanding of the crossflow that occurs around a coupled bluff body, such as a cylinder–plate system. In this system, both square and circular cylinders, which have wellknown vortex shedding mechanisms, have been widely utilized in investigations.
Teksin and Yayla (2016) conducted an experimental investigation of the flow structure in the wake region of a circular cylinder with an attached flexible plate. They discovered that lengthening the splitter plate resulted in a decrease in the peak value of any structure until a certain critical length ratio was reached. Beyond this critical length ratio, no significant variations could be observed. Binyet (2017; 2018; and 2020) conducted a series of studies on a cantilevered flexible plate placed in the wake of a square cylinder. They presented a deflection mechanism that relies on the unequal pressure distribution caused by traveling vortices on either side. In addition, they demonstrated that the presence of the plate has a significant impact on the wake region behind the cylinder. Furthermore, they found that the plate width is a critical factor that greatly affects the interaction between the plate and wake. Sharma and Dutta (2019; 2020) studied the flow over a square cylinder with an attached cambered flexible wake splitter. Their findings indicate that the cambered flexible wake splitter stimulates asymmetric oscillations due to the presence of vortex shedding over the square cylinder. Additionally, they investigated the flow characteristics in the intermediate flow region. Their study demonstrates that the wake frequency and mean drag coefficient exhibit nonmonotonic variations in response to changes in splitter plate length. This length is a critical parameter that can significantly affect wake transition. Dash (2020) explored the potential for drag reduction and wake regime control of a square cylinder through the implementation of an attached dual splitter plate. Their findings indicate that this particular configuration effectively mitigates the occurrence of Kármán vortex shedding and lift force fluctuation, leading to a notable increase in drag reduction. Wang (2022) conducted research on the aerodynamic performance and nearwake flow structure of a square cylinder with a flexible plate attached to its free end. Their findings revealed that the flapping plate tends to transform the alternating spanwise vortex shedding from both sides of the cylinder into symmetrical shedding, particularly in the vicinity of the free end. Zhang (2022) conducted a study of separated flows around a hybrid rigid–flexible thin plate. Their findings demonstrate that the presence of a thin plate results in a significant reduction of drag for the blunt body. Moreover, the vortices shed from the flow exhibit considerable irregularity and instability when subjected to typical flow velocities and blunt body shapes. To successfully apply piezoelectric materials through FIV, it is crucial to have a deep understanding of the flapping vortex dynamics of coupled multiple sidebyside flexible plates behind multiple bluff bodies. Unfortunately, there has been limited research conducted on this topic, specifically the behavior of multiple flexible plates submerged in the wake of square cylinders. As a result, further investigation is necessary to fully comprehend this important area of study.
In the present study, the flow past two flexible plates that are positioned side by side and submerged in the wake of two square cylinders is examined numerically. A comprehensive series of numerical simulations are conducted using a coupled computational fluid dynamics (CFD)–computational solid dynamics (CSD) method. The objective of the study is to examine the flapping vortex dynamics under varying gap spacings between the two flexible plates. Emphasis is put on analyzing the vibrational amplitude, the frequency characteristics, and the dynamic processes underlying the various flapping modes.
II. GOVERNING EQUATIONS
A. Governing equations for CFD and turbulence model
1. Governing equations for CFD
2. Detached eddy simulation (DES)
B. Governing equations for CSD
C. Governing equations for interface
D. Governing equations for new omega vortex identification method
III. NUMERICAL VERIFICATION
To verify the current FSI numerical scheme, the hydrodynamic behavior of a single flexible plate submerged in the wake of a square cylinder is simulated and verified. The computational configuration is depicted schematically in Fig. 1. At the inlet, a velocity profile of uniform magnitude is established, whereas at the outlet, a Neumann boundary condition is imposed to govern the flow behavior. The lateral boundaries of the flow domain are subject to a noslip boundary condition, ensuring that the fluid adheres to the boundaries and experiences zero velocity. The Reynolds number, which quantifies the ratio of inertial forces to viscous forces in a fluid flow, is calculated using the cylinder length D and the flow velocity, resulting in a value of 332.6. The physical parameters for the fluid and the plate are listed in Table I.
.  .  .  .  Young’s .  . 

.  ρ_{f} .  μ_{f} .  ρ_{s} .  modulus E .  Poisson’s . 
Parameter .  (g/cm^{3}) .  [g/(cm s)] .  (g/cm^{3}) .  [g/(cm s^{2})] .  ratio λ . 
Value  1.18 × 10^{−3}  1.82 × 10^{−4}  0.1  2.5 × 10^{6}  0.35 
.  .  .  .  Young’s .  . 

.  ρ_{f} .  μ_{f} .  ρ_{s} .  modulus E .  Poisson’s . 
Parameter .  (g/cm^{3}) .  [g/(cm s)] .  (g/cm^{3}) .  [g/(cm s^{2})] .  ratio λ . 
Value  1.18 × 10^{−3}  1.82 × 10^{−4}  0.1  2.5 × 10^{6}  0.35 
The time step is set to ∆t = 1 ms. The dimensionless parameters used for analysis are the dimensionless displacement y* = y/D, dimensionless vibrational frequency f_{n} = fD/v_{∞}, and dimensionless time T = tv_{∞}/D. The dimensionless lift coefficient is defined as $ C L = F y / 1 2 \rho f v \u221e 2 A $, where F_{y} is the force acting on the plate in the y direction and A is the area.
Both the fluid and solid regions are divided into smaller components using hexahedral elements. For the CSD simulation, a total of 48 416 elements and 219 012 nodes are utilized to accurately represent the solid domains. For the CFD simulation, a meshindependence study is performed on four different fluid mesh systems, designated as G_{1} (coarse), G_{2} (medium), G_{3} (fine), and G_{4} (very fine). The tested mesh parameters and the computed values are listed in Table II. It is determined that the differences in both the tip displacement and frequency between mesh systems G_{3} and G_{4} are insignificant. Consequently, the G_{3} fluid mesh system is selected for the verification test. The meshes for the fluid and solid regions are shown in Fig. 2.
Grid .  Elements .  $ y B * $ .  f_{n} . 

G_{1} (coarse)  71 640  1.05  0.0577 
G_{2} (medium)  118 720  1.06  0.0583 
G_{3} (fine)  142 470  1.06  0.0585 
G_{4} (very fine)  207 407  1.06  0.0585 
Grid .  Elements .  $ y B * $ .  f_{n} . 

G_{1} (coarse)  71 640  1.05  0.0577 
G_{2} (medium)  118 720  1.06  0.0583 
G_{3} (fine)  142 470  1.06  0.0585 
G_{4} (very fine)  207 407  1.06  0.0585 
The illustration in Fig. 3 showcases the evolution of vortex structures and plate deformation over time. On observing the flow around the square cylinder, it becomes evident that the shedding of vortices downstream leads to excitation and vibration of the plate, a phenomenon that has been previously documented by Matthies and Steindorf (2003) and Wood (2010).
The representation in Fig. 4 displays the development of the vertical displacements of the plate tip based on the G_{3} mesh. The results for the vertical displacement and frequency obtained in this study align closely with those in the literature, with a maximum difference of 6% as indicated in Table III. This difference is considered acceptable, indicating that the current numerical method can effectively solve FSI problems in cylinder–plate systems.
Source .  $ y B * $ .  Deviation (%) .  f_{n} .  Deviation (%) . 

Matthies and Steindorf (2003)  1.10  3.7  0.0610  4.3 
Wood (2010)  1.12  5.7  0.0573  2.0 
Present results  1.06  0.0585 
Source .  $ y B * $ .  Deviation (%) .  f_{n} .  Deviation (%) . 

Matthies and Steindorf (2003)  1.10  3.7  0.0610  4.3 
Wood (2010)  1.12  5.7  0.0573  2.0 
Present results  1.06  0.0585 
IV. PROBLEM DESCRIPTION
Figure 5 illustrates the configuration of two sidebyside flexible plates submerged in the wake of a square cylinder. This configuration consists of two fixed, rigid square cylinders separated by a distance S, while the flat plate is designed to be flexible. Both the fluid medium and the flexible plate share the same properties as those previously demonstrated in the verification example. The aim of this study is to investigate the vibrational behavior of a flexible plate located behind two square cylinders, with varying distances between them $ 0.5 \u2264 S / D \u2264 2.5 $. The fluid calculation domain is divided into a mesh, as depicted in Fig. 6, with the minimum mesh size being consistent with the G_{3} mesh used in the verification example. The investigation focuses on exploring the relationship between the plate’s vibration and the varying distance between the two square cylinders.
V. RESULTS AND DISCUSSION
This study explores the interrelated movement of two sidebyside flexible plates located behind square cylinders, along with the associated flapping vortex dynamics, for varying distances between the cylinders. The numerical simulation results reveal that the two flexible plates can exhibit four distinct modes of coupled motion: outofphase flapping, inphase flapping, transition flapping, and decoupled flapping. The following subsections provide a comprehensive analysis and discussion of the four coupled motion modes identified.
A. Outofphase flapping mode
Figure 7(a) displays the dimensionless vertical displacement over time for two monitoring points located at the ends of the flexible plates, under condition S/D = 0.5. It is evident from the figure that when S/D = 0.5, the flow state behind the square cylinders becomes complex owing to the close proximity of the cylinders. This results in a strong interaction between the flexible plates, making it challenging to establish a stable motion pattern. Over time, the vibrational amplitudes of plates 1 and 2 gradually converge to a constant value, with opposite directions of vibration, resulting in an outofphase flapping mode. In this mode, the vibrational amplitude of the plates behind the two square cylinders is significantly smaller than that of a plate behind the single square cylinder. This difference is attributed to the close proximity of the two flat plates, leading to a stronger interaction between them, which restricts the flapping ranges of both plates. In Fig. 7(b), the time histories of the lift coefficients for both plates are depicted. It can be observed that there are variations in the lift coefficient changes between the two plates, but the overall trends are similar. Plate 1 exhibits frequent fluctuations in lift coefficient values ranging from −2 to 4, while plate 2 shows fluctuations between −5 and 2. This discrepancy is a consequence of the close spacing between the plates and the pronounced interaction between them.
A spectral analysis of the dimensionless vertical displacements and lift coefficients of plates 1 and 2 was performed, and the results are depicted in Fig. 8, from which it can be deduced that the dimensionless principal frequency of the vertical displacement for plates 1 and 2 is 0.1267, which diverges from the dimensionless principal frequency of vertical displacement for a single plate (0.0585). This difference highlights that the interaction between the two plates will have an impact on the vibrational behavior of the plates. Over time, the vertical displacements of both plates undergo a transition to a quasiperiodic pattern. This indicates that the outofphase flapping mode gradually transforms into a firstorder bending vibrational mode.
Figure 9 displays the phase trajectories of the two flexible plates. It can be seen that these trajectories are diagonal for both plates, indicating that over time, the vibrations of the two plates approach stability and exhibit periodic coupled motion with a consistent pattern of vibration.
Figure 10 illustrates the distribution of vortex structures in the x–y plane at eight distinct instants (S/D = 0.5). The figure highlights that a fully formed vortex structure develops on either side of each square cylinder, accompanied by periodic shedding. This, in turn, triggers vibration of the flexible plate. As time passes, the vibrating plate affects the progression and shedding of the vortex, leading to a cyclical vortexinduced vibrational process. During the time period from T = 256.5 to T = 264.195, the interaction between the two flexible plates is intense, and a clear state of reverse flapping is observed. During this process, vortices are shed alternately on the upper side of plate 1 and the lower side of plate 2. As these vortices detach, the value of Ω steadily increases, making the vortices more pronounced. Additionally, a vortex with a relatively small Ω value emerges between the two plates, but gradually dissipates over time. This dissipation is attributed to the strong interaction between the closely spaced square cylinders, which hinders the stability of the vortex between the plates. During the time period from T = 266.76 to T = 274.455, both plates 1 and 2 exhibit an outofflapping mode, similar to the previous stage. During this stage, noticeable vortex shedding occurs on the upper side of plate 1 and the lower side of plate 2. Simultaneously, weak vortices are observed between the two plates owing to their interaction, eventually dissipating.
B. Inphase flapping mode
As the distance between the two square cylinders increases to 1D, the two flexible plates become noticeably synchronized in their flapping motion. Figure 11(a) depicts the dimensionless vertical displacements over time of two monitored points located at the tail ends of the flexible plates, under condition S/D = 1. As the distance between the plates gradually increases, the interaction between them begins to weaken, causing the motion of the two plates to tend toward a periodic state. Figure 11(b) displays the time histories of the lift coefficients corresponding to each plate. The data shown in the figure indicate that the lift coefficient of plate 1 fluctuates between 0 and 6, while that of plate 2 fluctuates between −6 and 0. The loads on plates 1 and 2 gradually exhibit symmetry.
To further analyze the relationship between fluid load and plate motion, the dimensionless vertical displacements and lift coefficients of plates 1 and 2 are analyzed using frequency spectra. In Fig. 12(a), the dimensionless main frequency of the vertical vibrational displacement and lift coefficient for plate 1 is 0.115. Similarly, in Fig. 12(b), the vertical vibrational displacement and lift coefficient for plate 2 also share the same dimensionless main frequency of 0.115. This indicates a symmetrical vibrational pattern emerging between the two plates under this specific operating condition.
Figure 13 presents the phase trajectories of the vibration of the two flexible plates for S/D = 1. It can be seen that the two plates swing back and forth under the influence of the fluid load, with the trajectories exhibiting a clear annular shape, indicating that the coupling action between fluid and structure is gradually becoming stable.
Figure 14 illustrates the distribution of vortex structures in the x–y plane at eight distinct instants (S/D = 1.0). As can be seen, at this stage, plates 1 and 2 swing up and down simultaneously, exhibiting an inphase flapping mode. During the time period from T = 256.5 to T = 264.195, the shedding vortex I on the lower side of plate 1 gradually ascends owing to the inertia of the plate’s upward vibration. During this ascent, vortex I gradually collapses and splits into two smaller vortices, significantly affecting the fluid load on plate 1 and causing a variation in the lift coefficient. On collapse of vortex I, a new vortex II forms between the two plates and gradually detaches, flowing downstream. During the time period from T = 266.76 to T = 269.325, vortex II becomes elongated owing to the downward swinging of plate 1. During this process, the value of Ω decreases steadily, resulting in a continuous decrease in vortex intensity. During the time period from T = 269.325 to T = 274.455, a new vortex III appears on the lower side of plate 1. Additionally, under the influence of plate 2’s swinging, vortex III begins to move upward. An examination of the flow field between the two plates during the entire swinging process reveals that the lower side of plate 1 has a greater tendency to form a stable shedding vortex, while the upper side of plate 2 is more susceptible to the creation of a weak attached vortex. This complexity in the vortex formation on plate 1 leads to stronger variations in fluid load and a more chaotic lift coefficient. On the other hand, the flow field of plate 2 more closely resembles the alternating vortex shedding state of a single plate, leading to more consistent fluid load changes and a lift coefficient that tends toward periodic fluctuations.
C. Transition flapping mode
As the distance between the two parallel square cylinders increases gradually to 1.5D, the flexible plates between them start to display a periodic swinging pattern. This marks a transition from aperiodic to periodic motion. Figure 15(a) illustrates the temporal evolution of the dimensionless vertical displacement of two monitoring points at the ends of the flexible plates, under condition S/D = 1.5. Over time, both plates converge toward stable periodic simple harmonic vibration, although there still exist differences compared with the vibration of a single plate. Figure 15(b) depicts the time histories of the lift coefficients for both plates. The lift coefficient of plate 1 oscillates between −2 and 4, and that of plate 2 between −4 and 2. Throughout the process, the lift coefficients of both plates ultimately tend toward a more stable periodic variation, indicating that the fluid loads on the plates become increasingly stable. This phenomenon can be attributed to the increasing separation between the square cylinders, leading to reduced interaction between the two plates and thus a more consistent vibrational pattern.
A spectral analysis of the dimensionless vertical displacement and lift coefficient of plates 1 and 2 was performed, and the results are depicted in Fig. 16. It can be seen that during the transition state, both the vibrational patterns of the plates and the fluid load undergo significant changes, with many highamplitude dimensionless vibrational frequencies being present in the dimensionless vertical displacement and lift coefficient. The spectrogram in Fig. 16(a) reveals dual main frequency characteristics for the dimensionless vibrational vertical displacement and lift coefficient. These include the same dimensionless main frequency as for S/D = 1.0 (0.115) and a secondary main frequency of 0.131. This indicates a shift in the swinging process between the two plates as the distance between the square cylinders increases.
Figure 17 depicts the phase trajectories of the vibration of the two flexible plates for S/D = 1.5. The figure reveals the reciprocal motion of the two plates under the influence of the fluid load, with the formation of a distinct ring shape. In comparison with the previous stage, the phase trajectories exhibit noticeable symmetry, indicating an escalating interaction between the plates as the distance between the square cylinders increases, which then gradually diminishes.
Figure 18 illustrates the distribution of vortex structures in the x–y plane at eight distinct instants (S/D = 1.5). It can be concluded from the figure that as the distance between the square cylinders increases gradually to 1.5D, the upper and lower sides of plates 1 and 2 will create an evident vortex structure during the swinging process. This vortex structure will gradually dissipate downstream, resulting in a flow state similar to that of a Kármán vortex street. In this process, the vortex structure that falls off from both sides of each plate is comparable to the vortex structure generated by a single plate swinging in the wake behind a square cylinder. However, before it detaches from the tail end of the plate, the vortex structure between the inner sides of the two plates is smaller in scale and has a lower Ω value than the vortex outside the plates. This is due to the interaction between the two swinging plates, which significantly restricts the growth of the vortex between their inner sides. Although both sides of the two flexible plates clearly exhibit alternating vortex shedding, the flow field characteristics of the two plates remain asymmetrical in shape owing to the differences in vibrational patterns between the plates.
D. Decoupled flapping mode
As the distance between the two parallel square cylinders increases to 2D, the interaction between the two flexible plates decreases and the vibration of each plate approaches a periodic simple harmonic motion. This resembles the vortexinduced vibration of a single plate and exhibits a state of decoupling. In other words, the motions of the two flexible plates are no longer linked, and they oscillate at their own independent amplitudes and frequencies. Figure 19 displays the time histories of the dimensionless vertical displacements and lift coefficients at two monitoring points located at the tail ends of the flexible plates under condition S/D = 2. According to the overall trend of vibration, plates 1 and 2 gradually start to exhibit simple harmonic motion under the influence of the flow field, and the amplitude of the dimensionless vertical displacement steadily increases before reaching a stable state. Both plates swing periodically in a firstorder bending mode, displaying a swinging state very similar to the vortexinduced vibration of a single plate. Furthermore, the swing process of plates 1 and 2 exhibits symmetry. As the lift coefficient changes, it can be observed that it eventually settles into a stable periodic change over time, with the fluid loads on plates 1 and 2 also being symmetrical. This indicates that the motion states of the two plates have to a certain extent become independent of each other.
A spectral analysis of the dimensionless vertical displacements and lift coefficients of the two plates has been performed, and the results are presented in Fig. 20. The spectrogram in Fig. 20(a) reveals a high degree of consistency between the dimensionless vertical displacements and lift coefficients of plates 1 and 2. Both exhibit a shared dimensionless main frequency of 0.123. Under this condition, the plates vibrate in a transitional state, and the dualfrequency characteristics also vanish. The process gradually shifts back to individual plate vibration, indicating a relatively independent state.
Figure 21 displays the phase trajectories of the vibration of two flexible plates for S/D = 2. The phase trajectories in this state are depicted as an overlapping ring, which indicates that the coupling between the flow field and structure has stabilized. Additionally, the two trajectories exhibit clear symmetry characteristics, further demonstrating that the two plates are in a stable, independent motion state.
Figure 22 illustrates the distribution of vortex structures in the x–y plane at eight distinct instants (S/D = 1.5). During the time period from T = 256.5 to T = 259.065, plates 1 and 2 exhibit their maximum vertical displacements by swinging in opposite directions. During the time period from T = 259.065 to T = 269.325, both plates simultaneously swing inward toward each other and return to a horizontal position. Finally, during the time period from T = 269.325 to T = 274.455, the two plates swing outward once again. In this process, two distinct forms of vortex shedding will occur. First, symmetrical vortex clusters will gradually form between the two plates, and as they move away from the plates, they will eventually take on a mushroomshaped vortex structure, as outlined by the red circles. Second, the vortex outside the plate will gradually move downstream in tandem with the movement of the plate, eventually falling away from its position far from the plate, as outlined by the blue boxes.
To verify the decoupling state, the distance between the two parallel square cylinders is increased to 2.5D. The dimensionless vertical displacement sand lift coefficients of plates 1 and 2 are obtained, as illustrated in Fig. 23. As can be seen, when both plates are stable, periodic harmonic motion occurs. The fluid load on the plates also exhibits a periodic pattern of change. In comparison with Fig. 19, it can be observed that the vibrational state and load change pattern of the plates when S/D = 2.5 are consistent with those when S/D = 2, demonstrating that increasing the distance between the two square cylinders does not affect the motion pattern of the plates. This highlights the fact that the motions of the two flexible plates are uncoupled, demonstrating a state of decoupling. The results of a spectral analysis of the dimensionless vertical displacements and lift coefficients of the two plates are presented in Fig. 24. Additionally, Fig. 25 displays the phase trajectories of the vibration of the two flexible plates under condition S/D = 2.5. Both figures demonstrate similar trends when compared with condition S/D = 2.0.
Figure 26 illustrates the distribution of vortex structures in the x–y plane at 12 distinct instants (S/D = 2.5) to analyze the development of vortex structure in the decoupling mode. During the time period from T = 256.5 to T = 261.63, both plates 1 and 2 swing inward simultaneously. The vortices a, b, c, and d that detach from both sides remain close to each other under the influence of the swing’s inertia force. This initiates the formation of the mushroom vortex structure, referred to as vortex I. At the next time T = 264.195, plates 1 and 2 swing outward simultaneously, and vortex I moves downstream. With the progression of time, the influence of the plates’ swing on the vortex decreases, and the vortex moves downstream as a stable mushroom vortex I. During the time period from T = 261.63 to T = 264.195, plates 1 and 2 swing outward together, causing the vortices L1 and L2 to separate from the outside of the two plates under the influence of the plates’ swinging inertia force. Similarly, from T = 271.89 to T = 277.02, plates 1 and 2 swing inward, leading to the formation of a stable mushroom vortex II, which then moves downstream in a stable shape. In the interval from T = 277.02 to T = 279.585, the vortices L3 and L4 detached from the outsides of the two plates fall away from the plates owing to the inertia force of the plates swinging outward. Therefore, the creation of the mushroomshaped vortex between the two plates results from the plates’ swinging inward inertia force. When the two plates swing outward, a symmetrical vortex will form at a distance from the plates.
VI. CONCLUSIONS
In this study, a twoway FSI algorithm based on tight coupling has been employed to simulate the flapping vortex dynamics of two sidebyside flexible plates submerged in the wake of a square cylinder. The variations of flapping modes, motion frequency, amplitude, and other parameters associated with the spacing between the two square cylinders have been analyzed. Additionally, the mechanism by which the spacing influences the coupled motion has been systematically investigated. The conclusions drawn from this research are as follows:

As the distance between two flexible plates behind parallel square cylinders increases, the vibration of the plates exhibits four different flapping modes: outofphase flapping, inphase flapping, transition flapping, and decoupled flapping. The interaction between the square cylinders affects the flapping mode of the plates. At a distance of 0.5D, the interaction is strong, resulting in an outofphase flapping mode. As the distance increases to 1.0D, the interaction becomes weaker, and the flapping mode becomes inphase. At 1.5D, the interaction continues to weaken, causing the plate vibrations to transition from nonperiodic to periodic. At 2.0D, the interaction disappears, and both plates swing independently with stable periodic motion. Finally, at 2.5D, there is no interaction between the plates.

Each flapping mode shows distinct vortex patterns. Outofphase flapping struggles to maintain stable vortices, which often fall off the outer sides of the plates. Inphase flapping allows clear vortex development, fading as the plates swing. Transition flapping creates uneven vortex shedding due to plate interaction, weaker inside than outside. Decoupled flapping forms stable, alternating vortices resembling a symmetric Kármán vortex street. These modes illustrate varied vortex behaviors during plate movement, ranging from unstable formations to symmetric patterns, showing how plate interaction influences vortex intensity and symmetry during flapping.
ACKNOWLEDGMENTS
The authors are grateful to the anonymous reviewers for their critical and constructive review of the manuscript. This study was cosupported by the National Natural Science Foundation of China (Grant No. U2106225), the Open Research Subject of Key Laboratory of Fluid and Power Machinery (Xihua University), Ministry of Education (Grant No. LTDL2023003), the Outstanding Youth Foundation of Jiangsu Province of China (Grant No. BK20211547), the Excellent Scientific and Technological Innovation Team Project in Colleges and Universities of Jiangsu Province [Grant No. SKJ(2021)1], the Natural Science Foundation of Jiangsu Province (Grant No. 21KJA470002), and the Senior Talent Foundation of Jiangsu University (Grant No. 18JDG034).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Bin Xu: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Writing – original draft (equal). Hao Wang: Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal). Weibin Zhang: Formal analysis (equal); Resources (equal); Visualization (equal). Yilin Deng: Resources (equal); Visualization (equal). Xi Shen: Validation (equal). Desheng Zhang: Funding acquisition (equal); Project administration (equal); Supervision (equal). B. P. M. (Bart) van Esch: Writing – review & editing (equal).
DATA AVAILABILITY
Data available on request from the authors.
NOMENCLATURE
 a and b

squares of Frobenius norms of A and B, respectively
 CFD

computational fluid dynamics
 CSD

computational solid dynamics
 D

diameter of square cylinder
 DES

detached eddy simulation
 E

Young’s modulus
 F

blending function
 f

body force
 FIV

flowinduced vibration
 FSI

fluid–structure interaction
 I

identity tensor
 k

turbulent kinetic energy
 p

pressure
 P

production limiter
 S

gap between two square cylinders
 T

dimensionless time
 t

time
 u

velocity
 x

space variable
 y

distance to nearest wall
Greek
Subscripts
NOMENCLATURE
 a and b

squares of Frobenius norms of A and B, respectively
 CFD

computational fluid dynamics
 CSD

computational solid dynamics
 D

diameter of square cylinder
 DES

detached eddy simulation
 E

Young’s modulus
 F

blending function
 f

body force
 FIV

flowinduced vibration
 FSI

fluid–structure interaction
 I

identity tensor
 k

turbulent kinetic energy
 p

pressure
 P

production limiter
 S

gap between two square cylinders
 T

dimensionless time
 t

time
 u

velocity
 x

space variable
 y

distance to nearest wall