3.1. Vibration modes in uncontrolled FSI system under varying stiffness

Numerical simulations of (2.3)–(2.8) are performed with a fixed Reynolds number, solid-to-fluid density ratio and Poisson’s ratio, while varying the non-dimensional Young’s modulus, such that

(3.1) \begin{equation} \textit{Re} = 150, \quad M = 1, \quad v_s = 0.35, \quad 2\times 10^4 \leqslant E_s \leqslant 6\times 10^6. \end{equation}

The selected parameter regime is chosen to explore a novel flapping regime, distinct from yet comparable to those in previous studies. Although we adopt the same Reynolds numbers as in Wu et al. (Reference Wu, Qiu, Shu and Zhao2014), our study employs a flexible splitter plate set-up instead of a filament. Furthermore, the range of Young’s modulus investigated is also novel, differing from that studied by Lee & You (Reference Lee and You2013) ( $1.61\times 10^4 \leqslant E_s\leqslant 9.66 \times 10^4$ ). The choice of $M = 1$ allows us to investigate instability mechanisms driven by FSI rather than by wake instabilities alone. Moreover, the regime investigated in this study corresponds to biological conditions observed in small aquatic organisms, such as ascidian and fish larvae (McHenry, Azizi & Strother Reference McHenry, Azizi and Strother2003; China et al. Reference China, Levy, Liberzon, Elmaliach and Holzman2017), which typically operate at low Reynolds numbers and density ratios near unity.

Figure 3.

Three types of plate responses under different Young’s moduli: (a) large-amplitude flapping (LAF) at $E_s = 4\times 10^4$ , (b) small-deflected flapping (SDF) at $E_s = 4\times 10^5$ , and (c) small-amplitude flapping (SAF) $E_s = 4\times 10^6$ . The top row (ai–ci) shows snapshot of the plate in saturated oscillation phases. The bottom row (aii–cii) presents the corresponding phase-space trajectories over time from $t = 0$ to 300.

Figure 4.

Snapshots of the vorticity field. Panel (a) is the solution before the onset of instability. Panels (b)–(d) correspond to the three representative flapping regimes in figure 3, i.e. LAF, SDF and SAF, respectively. The black lines with arrows represent the streamlines. Black lines with arrows indicate streamlines and blue arrows denote the flapping direction.

The simulations are initialised with a uniform zero flow field, and the inlet velocity is smoothly increased from zero to $U_{\infty }^*$ . As the flow develops, two symmetric (with respect to the $x$ axis) recirculating vortices form behind the cylinder, located above and below the elastic plate. These recirculating regions tend to slightly compress the elastic plate in the direction of $x \lt 0$ , leading the system to an undeflected equilibrium (UE) phase. Subsequently, self-excited instabilities emerge, giving rise to three different types of limit cycles, each corresponding to the selected $E_s$ : $4\times 10^4$ , $4\times 10^5$ and $4\times 10^6$ . To illustrate the characteristics of these three cases, figure 3 presents flow snapshots over one oscillation cycle in the saturated flapping phase (first row) along with phase portraits (second row) constructed from the plate tip displacement $\delta$ and the velocity ${\textrm d}\delta /{\textrm d}t$ over the time interval $0 \lt t \lt 300$ . For $E_s = 4 \times 10^4$ , an inflection point appears in the centreline deformation of the plate, and the maximum transverse deviation increases, causing the plate’s end to trace a characteristic figure-eight pattern as shown in figure 3(ai). This regime is referred to as large-amplitude flapping (LAF), characterised by displacement exceeding half the cylinder diameter. Moreover, the onset of instability exhibits characteristics of a Hopf bifurcation, with oscillatory behaviour emerging as the system becomes unstable, as shown in the phase trajectory in figure 3(aii). The system transitions to a limit-cycle state after $t \gt 160$ . Representative snapshots of the vorticity field during this regime are presented in figure 4(b). Two shear layers with opposite signs form on the upper and lower surfaces of the rigid cylinder, further shedding vortices in the downstream wake. The elastic plate interacts strongly with the larger vortex shedding, particularly near its trailing edge, where it appears to act as a ‘vortex cutter’, enhancing vortex detachment. This results in a ‘2S’ vortex shedding pattern that consists of two single vortices per cycle. The dominant flapping frequency in this regime is measured to be $f = 0.212$ . When the Young’s modulus increases to $E_s = 4 \times 10^5$ , small-deflected flapping (SDF) emerges, characterised by oscillations around an off-centre position deflected from the symmetry axis $y = 0$ . As illustrated in figure 3(bi), the plate deflects downward; however, a solution with upward deflection of the same magnitude can also be obtained. The selection between these two symmetric states is determined by the direction of the initial perturbation, which in a numerical simulation is implicitly introduced by factors such as the inherent micro-asymmetries of the computational mesh. The time evolution of $\delta$ shown in figure 3(bii) clearly displays the mean deflection. Between $t = 0$ and $t = 125$ , the plate tip gradually shifts downward without oscillation, a behaviour previously reported in numerical studies (Bagheri, Mazzino & Bottaro Reference Bagheri, Mazzino and Bottaro2012; Lācis et al. Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014). For $t \gt 125$ , oscillations emerge and grow exponentially, eventually saturating into a periodic limit cycle beyond $t \gt 220$ . The fundamental frequency is $f = 0.168$ and the oscillation amplitude is substantially smaller than that of the LAF mode. Consequently, the mean deflection exerts a significant stabilising effect on wake oscillations. Figure 4(c) presents snapshots of vorticity in the deflected limit cycle. The plate’s mean position consistently remains in the $y \lt 0$ region, with small oscillations superimposed. Vortices are released but with diminished intensity compared with LAF and at locations further from the plate surface. The ‘2S’ mode is absent, which suggests the overall more streamlined shape impedes vortex release in the wake. Moreover, a triangular lid-driven recirculation zone is formed between the rear of the cylinder and the plate, trapping a vortex $V^*$ that remains confined throughout the flapping cycle without escaping into the wake.

Figure 3(c) presents the small-amplitude flapping (SAF) regime at $E_s = 4\times 10^6$ . The oscillation initially grows exponentially, subsequently saturating into a periodic limit cycle after $t \gt 265$ . Unlike the SDF state, the plate returns to a symmetric flapping pattern, though with reduced amplitude and without centreline inflection, different from that observed under LAF conditions. The dominant flapping frequency is $f = 0.203$ . As illustrated in figure 4(d), the flow field also exhibits the ‘2S’ vortex shedding pattern, as $V^*$ is released and can enter the wake. A closer examination of the flow around the plate reveals that small secondary vortices emerge around the plate tip during motion: negative vorticity when the plate moves downward and positive vorticity when the plate moves upward. These vortices lack sufficient strength to separate and persist in the wake, but they generate enough downwash (or upwash) effects to influence the surrounding larger vortices.

Figure 5.

Characteristics of the three flapping regimes obtained from direct numerical simulation. For different values of $E_s$ , panel (a) shows the transverse tip displacement $\delta$ of the plate in the limit-cycle regime. The mean displacement is computed as $\delta _{\textit{mean}} = 1/2(\delta _{\textit{max}}+\delta _{\textit{min}})$ , while the amplitude $\delta _{\textit{amp}} = 1/2(\delta _{\textit{max}}-\delta _{\textit{min}})$ is indicated by error bars centred at the mean. Panel (b) presents the first two main flapping frequencies $f_0$ and $f_1$ . The LAF, SDF and SAF regimes are highlighted with light yellow and purple and blue. The critical transition Young’s moduli between the three regimes are $E_{s1} = 1.26 \times 10^5$ and $E_{s2} = 7.78 \times 10^5$ .

To further explore the system behaviour over a broad range of Young’s modulus values, figure 5 summarises the flapping characteristics of the elastic plate in the saturated limit-cycle regime. Figure 5(a) shows the transverse tip displacement $\delta$ and figure 5(b) presents the first two dominant flapping frequencies $f_0$ and $f_1$ obtained via a fast Fourier transform. Moreover, $f_1$ is the second harmonic of $f_0$ ( $f_1 = 2f_0$ ). As $E_s$ increases, the system transitions through three distinct dynamical regimes. For low stiffness ( $E_s \lesssim 1.26 \times 10^5$ ), the system remains symmetric about the midline, with large flapping amplitudes ( $\delta _{\textit{amp}} \gt 0.25$ ) and zero-mean displacement, corresponding to the LAF regime. Beyond this critical value is the SDF regime, where the plate oscillates around an off-centred mean position ( $\delta _{\textit{mean}} \lt 0$ ) with lower frequencies ( $f \lt 0.2$ ) and reduced amplitudes. As the stiffness further increases beyond $E_s = 7.78 \times 10^5$ , the system regains symmetry, entering the SAF regime, characterised by weaker oscillations and the absence of large-scale wake interactions. For clarity in the subsequent analysis, we select $E_s = 4 \times 10^4$ , $4 \times 10^5$ and $4 \times 10^6$ as representative parameters of the LAF, SDF and SAF regimes, respectively.

3.2. Koopman reduced-order modelling

Figure 6.

(a) Full input sequence $u_{{\textit{train}}}$ , where $D_1$ is used to construct the Koopman ROM and $D_2$ is used for validation. (b) A portion of $u_{{\textit{train}}}$ from $t= 200$ to $220$ ; the blue shaded region represents an input corresponding to $\mathcal{K}_{\bar {u}_2}$ . (c) Displacement response at probe P $_3$ under the full input shown in (a). (d) A portion of displacement response under the input $\bar {u}_2 = 0.15$ only. The abscissa $t^*$ is a redefined continuous-time coordinate.

To construct the Koopman ROM, we first extract training data from full-order FSI simulations. For this purpose, we employ the optimal sensor-actuator configuration detailed in Appendix B. As shown in figure 6(a), a long input sequence $u_{{\textit{train}}}$ is applied to the FOM and the resulting plate displacements at probe locations $\boldsymbol{\delta } = \{\delta _1,\delta _2,\delta _3\}$ during the post-transient period ( $t \gt 200$ ) are recorded. The collected data are then divided into a training set ( $D_1$ ) for model construction and a validation set ( $D_2$ ) for accuracy assessment. The temporal sampling interval between neighbouring data points is $\Delta t = 0.1$ , approximately $1/50$ of the flow’s dominant period. The training dataset contains $6500$ sampling points, ranging from $t = 200$ to $850$ , while the validation dataset contains $1500$ sampling points over $t \in [850,1000]$ . A detailed view of the training input $u_{{\textit{train}}}$ within $t = 200{-}220$ is shown in figure 6(b). The input is a piecewise-constant function taking values $\bar {u}_i \in \{0.3,0.15,0,-0.15,-0.3\}$ , with variable segment durations. The selection of the magnitude range for the control input $\bar {u}_i$ is governed by a critical criterion: the energy injected by the control input must be sufficient to overcome the energy of the self-sustained flapping. If this range is too narrow, the control action will have insufficient authority over the system’s dynamics. Consequently, instead of achieving the objective of complete flapping suppression, the system will transition to a new, control-induced limit-cycle state after the control is applied. In such cases, it is necessary to expand this range and retrain the control-oriented ROM. Furthermore, we recommend that the discrete values within the range of $\bar {u}_i$ be equally spaced to ensure that the resulting ROM can comprehensively capture the system’s nonlinear dynamics across the full control space.

The training data are then categorised into five separate autonomous Koopman systems according to the value of $\bar {u}_i$ . Figure 6(c) shows the total training data collected at probe P $_3$ for the LAF case. The imposed input disrupts the original periodic dynamics, forming an irregular sequence within the range $\delta \in [-0.5, 0.5]$ with significant amplitude variations. Figure 6(d) presents a portion of displacement data used to train the autonomous Koopman system $\mathcal{K}_{\bar {u}_2}$ corresponding to $\bar {u}_2 = 0.15$ , where we concatenate the non-continuous displacement data segments corresponding to the blue shaded region in figure 6(b) in a continuous-time coordinate $t^*$ . Due to the piecewise-constant nature of the input, the system output exhibits discontinuous transitions between distinct dynamic states. This design introduces various initial conditions for each segment, thereby enhancing the generalisation capability of the trained Koopman ROM. Additionally, a delay embedding method (Yuan et al. Reference Yuan, Zhou, Zhou, Wen and Liu2021) is employed to enhance the accuracy of DMD-type algorithms, i.e.

(3.2) \begin{equation} \boldsymbol{\delta }_n = [\delta (n), \delta (n-1), \ldots , \delta (n-d+1)]^T \in \mathbb{R}^d, \end{equation}

where $d$ is the dimension of time-delay embedding. The value of $d$ directly affects the accuracy of the Koopman ROM, which will be discussed in the following section.

Figure 7.

(a) Eigenvalues of the finite Galerkin matrix $\mathbb{K}$ . The colour scale indicates the residual associated with each eigenpair. (b) Residual $ \textit{res}({\lambda _{\!j}},{{v}_{\!j}})$ of the computed spectral information for the first 100 Koopman modes.

Koopman ROMs are then constructed for each autonomous system $\mathcal{K}_{\bar {u}i}$ using the training data. Figure 7(a) shows the eigenvalue distribution of the finite-dimensional $N \times N$ Galerkin matrix $\mathbb{K}$ for the autonomous system $\mathcal{K}_{\bar {u}_2}$ , obtained via the KDMD algorithm. A total of 256 eigenvalues are computed with a delay embedding dimension of $d=5$ . To assess the accuracy of these eigenvalues, we calculate the residual value for each eigenvalue according to (2.29). As shown, all eigenvalues lie within $(\lambda _r,\lambda _i) \in [-2.5,2.5]\times [-0.4,0.1]$ . Eigenvalues distant from the origin $(\lambda _r,\lambda _i) = (0,0)$ exhibit larger residuals, while those closer have smaller residuals. These high-residual modes correspond to large oscillation frequencies and growth/decay rates that are non-physical in the current FSI system and can therefore be classified as pseudo-modes. Figure 7(b) plots the relative residuals $ \textit{res} ({\lambda _{\!j}},{{v}_{\!j}} )$ for the first 100 eigenvalues. Notably, the residuals of these Koopman modes gradually increase with the modal order $k$ . The residuals of the first 32 modes are all less than 0.1, indicating higher accuracy in capturing the dominant system dynamics. Accordingly, we reconstruct the displacement using (2.27), retaining only low-residual modes to enhance the accuracy of the ROM. The relative reconstruction error is defined as

(3.3) \begin{equation} \varepsilon =\int _0^{T_0} \frac { |\delta _r - \delta _t| }{|\delta _t| } {\textrm d}t, \end{equation}

where $\delta _r$ and $\delta _t$ represent the reconstructed and true displacement, respectively; and $T_0 = 1.5$ is the final time we select.

Figure 8.

(a) Relative error $\varepsilon$ as a function of time-delay embedding order $d$ and the ROM rank $k$ . The yellow star at $(d,k) = (5,140)$ represents the parameters used in this study.

The rank $k$ of the Koopman ROM and the delay embedding dimension $d$ are the two most influential parameters governing ROM accuracy. Figure 8 presents a contour map illustrating how the ROM error $\varepsilon$ varies across delay embedding dimensions $d \in [0,9]$ and ROM ranks $k\in [20,200]$ . We first establish the necessity and effectiveness of delay embedding when using sparse measurements by fixing $k = 120$ . Without delay embedding ( $d = 0$ ), the ROM error remains substantially high (i.e. $\varepsilon = 0.32$ ). As delay embedding is introduced, ROM accuracy improves markedly, with error decreasing rapidly as $d$ increases for $0\lt d \leq 4$ . For $d\gt 4$ , the error $\varepsilon$ shows little further reduction across all $k$ , confirming the effectiveness of moderate delay embedding under sparse measurements.

Similarly, we fix $d=5$ and examine how the ROM error varies with model rank $k$ . Interestingly, the results reveal a non-monotonic trend: increasing ranks does not necessarily yield better accuracy. This is attributed to the inclusion of spurious modes at higher ranks (as illustrated in figure 7), which may degrade the performance of the ROM. In this study we construct the Koopman ROM using $k = 120$ and $d = 5$ (marked by an asterisk), yielding a reconstruction error $\varepsilon$ below 0.03. This choice offers a balanced trade-off between model accuracy and computational cost, as further increases in embedding dimension lead to considerable overhead without significant improvement in performance.

Figure 9.

Comparison between the true and reconstructed displacement responses at three probe locations. The top row (ai–ci) shows the results for (a) the LAF case with $\bar {u}_1 = 0.3$ , (b) the SDF case with $\bar {u}_3 = 0$ , and (c) the SAF case with $\bar {u}_5 = -0.3$ . The bottom row (aii–cii) presents the corresponding relative error for each probe.

Based on the selected modelling parameters, we construct Koopman ROMs for each combination of flapping mode (LAF, SDF, SAF) and control input level ( $\bar {u}_1$ to $\bar {u}_5$ ), resulting in a total of 15 autonomous systems. Figure 9 presents representative comparisons between the reconstructed displacement $\delta (t)$ with the full-order results within the range $t \in [0, 1.5]$ for three selected cases: the LAF case with $\mathcal{K}_{\bar {u}_1}$ , the SDF case with $\mathcal{K}_{\bar {u}_3}$ and the SAF case with $\mathcal{K}_{\bar {u}_5}$ . As shown in figure 9(ai–ci), the Koopman ROMs closely match the reference displacements at all three probe locations, especially in capturing the local peaks. Figure 9(aii–cii) displays the real-time relative error $(\delta _t - \delta _r)/\delta _m$ , where $\delta _m$ is the maximum displacement across the three probes in $t \in [0,1.5]$ . As expected for Koopman-based models, the reconstruction error increases with time. Therefore, we set the prediction horizon to $t_p = 1.0$ to ensure prediction accuracy within MPC, which is approximately one-fifth of the oscillation period. Within this range, the ROM maintains an average error below 5 %, which proves sufficient for effective closed-loop control. In our framework, each control cycle begins with the probes measuring the current system state $\boldsymbol{\delta }$ . This measurement serves as the initial condition for the Koopman ROMs, which then predicts the system’s future dynamics over a prediction horizon $t_p$ . An optimal control sequence is computed based on this prediction, but only the initial control action is executed on the full-order FSI system over a single control time step $t_c$ . This iterative loop continues throughout the control duration.

Figure 10.

Time evolution of the displacement response $\delta$ at the trailing edge of the plate (ai–ci) and the corresponding control input $u$ (aii–cii) for the (a) LAF, (b) SDF and (c) SAF cases. Control is activated at $t = 300$ (red dashed line), after the system reaches a fully developed state (grey shaded region). The five time instants $t_1 = 302$ , $t_2 = 310$ , $t_3 = 320$ , $t_4 = 330$ and $t_5 = 360$ (marked by yellow diamonds) are selected for detailed visualisation in figure 11.

3.3. Suppressing the plate’s flapping by ROM-MPC controller

Next, we apply the established Koopman ROM to control the full-order fluid–structure coupled system across all three vibration modes (LAF, SDF and SAF). The control objective is to suppress the flapping motion and stabilise the flap to an undeformed equilibrium state, i.e. $\boldsymbol{\delta }_{\textit{ref}} = 0$ . The control input $u$ is obtained by solving the optimisation problem (2.37), and it is updated every $t_c = 0.2$ control horizon, approximately one twenty-fifth of the oscillation period. Figure 10 presents the time evolution of response $\delta$ and control input $u$ before and after control activation for the LAF, SDF and SAF cases. The control is activated at $t = 300$ , by which time the flow has entered a fully developed state, i.e. exhibiting limit-cycle periodic motion. Notably, the control algorithm exhibits strong robustness to varying initial conditions, consistently achieving successful suppression. For the LAF case (figure 10 ai), the uncontrolled displacement amplitude is $\delta _{\textit{amp}} = 0.41$ . After activating control, $\delta$ undergoes a slight amplification cycle, reaching a maximum value of 0.45 before beginning to decay. This transient increase is attributed to the strong actuation injecting additional disturbance energy into the system, reshaping the mean flow field structure and thereby stimulating non-modal energy amplification. Similar responses have been reported in previous studies on jet control strategies (Li & Zhang Reference Li and Zhang2022; Jiang & Cao Reference Jiang and Cao2024). Figure 10(aii) shows the control input $u$ , which reaches its maximum value during the initial few cycles and gradually decreases as the plate’s oscillation decreases. Furthermore, we find that the dominant flapping frequency after control closely matches the frequency identified in the resolvent analysis, consistent with findings in Gross et al. (Reference Gross, Marks and Sondergaard2024). This suggests that resolvent analysis can also offer valuable guidance for controlling FSI systems. To visualise the control effect more intuitively, figure 11 presents flow snapshots at five representative time instants ( $t_1$ – $t_5$ ) during the control process. As shown in figure 11(ai–aiii), the blowing–suction actuators continuously modify the vortex structures above and below the plate, progressively altering the vorticity field and adjusting the plate’s deformation accordingly. The ‘2S’ vortex shedding mode is gradually suppressed from $t_1$ to $t_3$ , with vortices becoming narrower and more aligned with the flow direction. Simultaneously, the plate undergoes different deformation states under fluid loading. By $t_4$ and $t_5$ , the original ‘2S’ vortex shedding mode has been effectively suppressed, and the flow field evolves toward a more streamlined configuration with reduced and evenly distributed vorticity. The entire system exhibits an undeformed equilibrium state. For the SDF case, the uncontrolled plate exhibits biased flapping, with its displacement significantly deviating from the target equilibrium position (figure 10 bi). When the control is applied, the asymmetric blowing-suction jets (figure 10 bii) disrupt the original oscillation mode (figure 11 bi). After $t \gt 305$ , the plate crosses into the $\delta \gt 0$ region and undergoes another amplification cycle. Figure 11(bii–biii) shows that during this phase, the flow field becomes more unstable and begins to develop the ‘2S’ separated vortices. After $t \gt 320$ , the flow field enters an oscillation decay phase (figure 11 biv), with both plate displacement and control input $u$ returning to a symmetric oscillatory state. By $t \gt 330$ , the plate reaches an unbiased equilibrium, successfully achieving the control objective. The control behaviour in the SAF case resembles that of the LAF case. As shown in figure 10(c), the system exhibits symmetric decay of plate oscillations under control. The corresponding flow field snapshots in figure 11(c) further confirm this trend, showing a monotonic suppression of vortex shedding and flow stabilisation.

Figure 11.

Five snapshots of vorticity contours during control: (a) LAF case, (b) SDF case and (c) SAF case corresponding to figure 10.

To illustrate the evolution of coherent flow structures and stability characteristics before and after control, we perform Koopman modal analysis for the three flapping cases. To extract Koopman modes, we collect 500 flow field snapshots at a time interval of $\Delta t = 0.1$ over $t = 300$ – $350$ . Reduced-order modelling using the full state of the flow field can achieve high accuracy without the need for time-delay embedding $(d=0)$ (Colbrook et al. Reference Colbrook, Ayton and Szőke2023). A comparison between the reconstructed flow field and the ground truth reveals that the average error is less than 1 %. Figure 12(a–c) presents the first 13 eigenvalues (frequency $f = \lambda _r /2\pi$ and growth rate $\lambda _i$ , sorted by residuals) before and after control. For all three cases, the leading four modes exhibit growth rates of the order of $\mathcal{O}(10^{-3})$ or smaller, indicating that these modes are neutrally stable. Here $M_0$ is the mean mode with $f = 0$ and the $M_1$ mode corresponds to the dominant oscillation frequency. It can be seen that, after control is applied, the eigenvalue distribution changes significantly: except the mean mode $M_0$ , all other modes exhibit increased decay rates and changing frequencies. As shown in figure 12(a), for the LAF case, the eigenvalue of $M_1$ shifts from $(\lambda _r/2\pi , \lambda _i) = (0.212, 1.83 \times 10^{-5})$ to $ (0.238, -1.58 \times 10^{-3})$ . Similar trends are observed for the SDF and SAF cases.

Figure 12.

Distribution of the leading 13 Koopman eigenvalues before (orange triangles) and after (blue rhombi) control for the (a) LAF, (b) SDF and (c) SAF cases. The hollow circles represent the mean mode $M_0$ and the filled markers denote the dominant oscillation mode $M_1$ .

Figure 13.

Comparison of key Koopman modes for the uncontrolled and controlled systems. The columns correspond to the (a) LAF, (b) SDF and (c) SAF cases. The top two rows visualise the mean flow mode ( $M_0$ ), while the bottom two rows visualise the dominant oscillatory mode ( $M_1$ ). For each mode, the figure compares the uncontrolled case (rows 1 and 3) with the corresponding controlled case (rows 2 and 4). The dashed boxes in (ai) indicate three streamwise regions: R1 ( $0.5 \lt x \leqslant 2.5$ ), R2 ( $2.5 \lt x \leqslant 6$ ) and R3 ( $6\lt x \leqslant 13$ ).

We focus on analysing the characteristics of the mean mode $M_0$ and the dominant oscillation mode $M_1$ . The first two rows of figure 13 illustrate the spatial distribution of $M_0$ for the LAF, SDF and SAF cases under both uncontrolled and controlled conditions. The downstream region is partitioned into three distinct zones: R1 (the immediate wake encompassing the plate, $0.5 \lt x \leqslant 2.5$ ), R2 (the intermediate wake, $2.5 \lt x \leqslant 6$ ) and R3 (the far wake, $6 \lt x \leqslant 13$ ). Vortices are generated in R1, amplified in R2 and subsequently advected into the far-field region R3. In the LAF case (figures 13 ai and 13 aii), the mean mode undergoes significant alterations across all three regions after control is applied. Post-control, the vorticity becomes more concentrated near the flexible plate in the latter part of region R1 ( $1.5 \lt x \leqslant 2.5$ ). Furthermore, a symmetric and elongated shear layer emerges in regions R2 and R3. This phenomenon occurs because the control suppresses both the flapping motion and vortex shedding. Conversely, for the SDF case, region R1 shows minimal change between the uncontrolled and controlled states (figures 13 bi and 13 bii). In contrast, regions R2 and R3 exhibit an opposite trend compared with the LAF case: a shear layer prominent in the uncontrolled flow vanishes after control is activated. This is because the control action transiently disrupts the deflected state, compelling the system to transition towards a symmetric configuration. The SAF case (figures 13 ci and 13 cii) is analogous to the LAF case, but the post-control shear layer in regions R2 and R3 is wider than that observed in the controlled LAF case (figure 13 aii). Collectively, this comparison demonstrates that the application of control induces significant structural changes in the $M_0$ mode, confirming that the proposed strategy fundamentally alters the time-averaged characteristics of the flow field.

Next, we analyse the dominant oscillation mode $M_1$ . For the LAF case, figures 13(a iii) and 13(aiv) show that before and after control, the morphological characteristics of the vorticity in the three regions remain almost unchanged, but there are some changes in amplitude. The mode amplitude in the latter half of region R2 and in region R3 is lower than before control, indicating that the shear and separation effects of the plate are weakened during the control process. For the SAF case, as shown in figures 13(ciii) and 13(civ), the dominant vortex pattern exhibits noticeable changes mainly within region R1, reflecting the effect of jet actuation on near-wall travelling waves. We observe that the modifications of $M_1$ in these two cases are relatively minor compared with the SDF case. Figures 13(b iii) and (biv) show the SDF case. Before the control, the dominant vorticity within $M_1$ mode is primarily concentrated in R2 and R3, exhibiting a typical Kármán vortex street resembling that of a pure fluid system. After control is applied, the $M_1$ mode undergoes significant changes in all three regions. Specifically, the vorticity increases markedly in R1, while the vortex patterns in R2 and R3 become broader and more irregular compared with the uncontrolled case. In our previous work on using balanced proper orthogonal decomposition (BPOD)-ROM for cylinder VIV control (Jiang & Cao Reference Jiang and Cao2024), we found that control either modified existing modes or created new ones, and additional energy is required to sustain the newly induced modes. A similar pattern appears here: for LAF and SAF, the Koopman-based control alters the original mode $M_1$ , while for SDF, it generates a new dominant mode.

We further analyse the control mechanism from an energy transfer perspective. For the present FSI system, we calculate the energy variations of the plate and the flow field kinetic energy relative to the undeformed equilibrium state during the control process, as well as the energy injected into the flow field by actuators. First, for the Saint-Venant–Kirchhoff material, the elastic potential energy and kinetic energy densities can be calculated as

(3.4) \begin{equation} {{w}_{p}}=\frac {{{\lambda }_{s}}}{2}{{\left \{ {tr}\left [ \boldsymbol{E}({{\eta }_{s}}) \right ] \right \}}^{2}}+{{\mu }_{s}}{\textit{tr}}\big [\boldsymbol{E}{{({{\eta }_{s}})}^{2}} \big ], \qquad w_k = \frac {M}{2}\left (\frac {{\textrm d} \eta _{s}}{{\textrm d}t} \right )^2\!, \end{equation}

where $w_p$ is the elastic potential energy density and $w_k$ is the kinetic energy density. Total energy variations for structure and flow field are

(3.5) \begin{equation} W_s = \int _{\varOmega _s} (w_p + w_k) {\textrm d}\varOmega - \int _{\varOmega _s} \bar {w}_p {\textrm d}\varOmega , \qquad W_{\kern-1pt f} = \int _{\varOmega _f} \frac {\boldsymbol{u}^2}{2} {\textrm d}\varOmega -\int _{\varOmega _f} \frac {{U}^2}{2} {\textrm d}\varOmega , \end{equation}

where $\bar {w}_p$ and $U$ represent the UE state. Furthermore, energy injected into the flow field by the jet for $t_0$ – $t_e$ are

(3.6) \begin{equation} W_i = \int _{t_0}^{t_e} \left \{\int _{\varGamma _{\!j}} \frac {\left [{Q_{t} h(\boldsymbol{u}_{\!j})}\right ]^2}{2} {\textrm d}x\right \} {\textrm d}t, \end{equation}

where $\varGamma _{\!j}$ is the surface of the jet on the cylinder, $h(\boldsymbol{u}_{\!j})$ represents the jet velocity profile and $Q_t$ is the time-varying magnitude of the jet. We define two energy metrics: the cumulative injected energy, $W_i^t$ , which is the total kinetic energy injected from a fixed start time $t_0 = 300$ to the current time $t_e = t$ ; and the single-injected energy, $W_i^p$ , representing the kinetic energy from a single jetting action over a sliding time window of fixed duration $t_c$ (where $t_e - t_0 = t_c$ ).

Figure 14.

(a) Evolution of cumulative jet energy $W^t_i$ and total flow field energy $(W_{\kern-1pt f} + W_s)$ from $t = 300$ to $400$ for the LAF, SDF and SAF cases, with colour indicating temporal progression. Panels (b), (c) and (d) show time evolution of the actuation jet energy $W^p_i$ for the LAF, SDF and SAF cases, respectively.

Figure 14(a) illustrates the relationship between the accumulated jet energy $W_i^t$ from actuators and the flow field energy ( $W_s + W_{\kern-1pt f}$ ) over time. There exists a time delay between the initial application of control and its effectiveness. For the LAF case, during the initial control phase ( $300 \lt t \lt 330$ ), the flow field energy maintains a magnitude comparable to the uncontrolled state. When $330 \lt t \lt 360$ , the control effect begins to manifest as $W_s + W_{\kern-1pt f}$ decreases linearly with an approximate slope of $-1510$ , and $W_i^t$ falls within the range $[1.4, 1.8]$ . When $t \gt 360$ , the flow field energy continues to decline, while the jet energy approaches a limiting value with minimal temporal variation, $W_i^t \rightarrow 1.81$ . Similar behaviour is observed in the SAF case, with the key distinction being that a higher final jet energy expenditure of $W_i^t \rightarrow 2.2$ is required to achieve significant attenuation. For the SDF case, $W_s + W_{\kern-1pt f}$ initially maintains the smallest magnitude compared with the other two cases. When $2 \lt W_i^t \lt 2.3$ , the flow field energy begins to increase rapidly before experiencing substantial linear decay. We find that the SDF case requires the highest total jet energy of $W_i^t \rightarrow 2.75$ to completely suppress instability. Figures 14(b), 14(c) and 14(d) illustrate the time evolution of $W_i^p$ for each jetting action after time $t \gt 350$ . It can be observed that the jetting energy maintains a small magnitude during this period. This is because all three cases have reached a stable equilibrium state accompanied by very weak flapping, as shown in figure 10. Consequently, $W_i^p$ also gradually decays as the flapping motion attenuates. The role of the controller at this stage is to maintain this stable state and prevent the re-emergence of LAF triggered by the inherent instability of the FSI system. Based on the modal decomposition results presented above, flow control in the SDF case operates by generating new modes that consume more energy, while the LAF and SAF cases modify the original inherent modes, consuming less energy. This conclusion aligns with the findings of Jiang & Cao (Reference Jiang and Cao2024).

Figure 15.

Time evolution of elastic potential energy $W_p$ and kinetic energy $W_k$ of the plate in a range of $t \in [290,330]$ for the (a) LAF case, (b) SDF case and (c) SAF case.

Figure 16.

(a) Time evolution of $W_p$ , $W_k$ and $W_s$ from $t = 296$ to 310 for the SDF case. (b) Leading instability eigenmode calculated by global linear instability analysis. (c)-(h) Snapshots of the flow field at six representative time instances ( $t_1$ – $t_6$ ). Here $V^+$ and $V^-$ denote vortices rotating in the counterclockwise and clockwise direction, respectively.

To better understand the plate’s energy evolution during control, figure 15 presents the time history of elastic potential energy $W_p$ and kinetic energy $W_k$ . For the LAF case, the plate’s low Young’s modulus leads to relatively small elastic potential energy $W_p$ , approximately one-tenth of $W_k$ . Consequently, the plate’s energy is predominantly governed by its kinetic component as shown in figure 15(a). After control is applied, both $W_p$ and $W_k$ increase with time until $t \lt 308$ , with the increase in kinetic energy significantly exceeding that of potential energy. For the SDF case, the plate initially exhibits a biased flapping state with increased stiffness, storing unreleased elastic potential energy around $W_p = 4$ (figure 15 b). Due to the limited flapping amplitude, the kinetic energy remains low at $W_k = 0.2$ . Once control is applied, however, $W_k$ increases significantly, peaking at 20. For the SAF case, we observe that $W_p$ and $W_k$ exhibit nearly identical amplitudes but opposite phases, as shown in figure 15(c). This indicates a cyclical process where kinetic energy converts to potential energy, which subsequently releases and transforms back into kinetic energy. Moreover, we define energy efficiency as the ratio of the dissipated mechanical energy of the structure (the sum of its kinetic and potential energy) to the cumulative energy $W_i^t$ injected by the jets. It is found that our control strategy has a high energy efficiency.

Figure 16(a) shows the time evolution of elastic potential and kinetic energy from $t = 296$ – $310$ for the SDF case. At $t = t_1$ , the jet actuation induces vortex $V_1^+$ , which is initially locked at the junction between the cylinder’s lower right side and the plate, enhancing its intensity and advects it downstream. Figure 16(b) illustrates the plate’s transition from equilibrium to deflected state through a ‘symmetry-breaking’ mechanism, which explains the presence of the locked vortex $V_1^+$ as evidenced in figure 5(c). By $t = t_2$ , vortex $V_1^+$ has reached the plate’s trailing edge and escaped the previously confining region shown in figure 16(d). During this process, both $W_p$ and $W_k$ increase simultaneously as the plate oscillates downward. At $t = t_3$ , the plate reaches maximum deflection with $W_p = 6$ and a new vortex $V_2^+$ forms (figure 16 e). As the plate reverses direction, upper-side vortices are shed into the flow, completing one cycle (figure 16 f). At $t = t_4$ , the plate’s elastic potential energy begins to release while its kinetic energy increases, establishing a phase difference between $W_p$ and $W_k$ . With the original biased oscillation being disrupted, the plate enters the $y \gt 0$ region, and the shedding dynamics begin to resemble those of the SAF case. Figures 16(g) and 16(h) illustrate the alternating vortex shedding process corresponding to times $t_4$ and $t_5$ in figure 16(a). This sequence reveals the ‘breaking’ control mechanism that the actuation exerts on the SDF solution. Specifically, the control effectively disrupts the locked vortex structure that maintains the deflected state, thereby initiating a transition toward a more symmetric flapping pattern characterised by alternating vortex shedding.

Figure 17.

(a) An illustration of the RL framework. (b) Training process for the LAF case. Testing process: (c) time evolution of the displacement response $\delta$ and (d) the corresponding control input $u$ .

Having established the effectiveness of the Koopman-MPC, we now turn to its computational efficiency, particularly in comparison to RL approaches. We compare the computational efficiency of our method with RL approaches. In both frameworks, the most time-consuming component lies in the use of full-order FSI simulations: in our method, full-order simulation results are used to train Koopman ROMs, while in RL, they provide the simulation environment for agent training. Therefore, we compare the two methods by evaluating the amount of full-order simulation time required to derive effective control strategies. For a fair comparison, all computations in the present study were performed on an Intel(R) Xeon(R) Silver 4216 processor (32 cores @ 2.1 GHz). Further details regarding the RL framework implementation and Koopman ROM training can be found in Appendix C. For the LAF case, the training results under the RL framework of figure 17(a) are shown in figure 17(b). The proximal policy optimisation (PPO) agent typically learns the control strategy after approximately 120 episodes, requiring about 125 hours, significantly exceeding the training time of the Koopman ROM (10.5 h). Regarding the control performance, the RL controller begins by applying a stronger actuation, as shown in figure 17(d). This aggressive initial action causes a transient increase in displacement (figure 17 c) but achieves a more rapid initial amplitude decay compared with the ROM-based method. However, the decay rate slows in the later stages, resulting in a total control time comparable to the ROM method for a stricter objective, such as $\delta \lt 0.05$ . These results highlight the computational efficiency of the proposed Koopman-based approach. Nevertheless, this does not suggest that it can fully replace RL methods, especially as recent advances have greatly improved the data efficiency of RL algorithms (Liu & MacArt Reference Liu and MacArt2024; Jia & Xu Reference Jia and Xu2024).