1. Introduction
The flapping behaviours of elastic slender structures in fluid can be observed in natural contexts like flag flapping (Connell & Yue Reference Connell and Yue2007), leaf waving (Sader et al. Reference Sader, Cossé, Kim, Fan and Gharib2016) and paper shaking (Gurugubelli & Jaiman Reference Gurugubelli and Jaiman2015), to name a few. A slender structure in flow is commonly modelled as elastic beams subjected to viscous flow. The interaction of a beam with unsteady viscous flow is a fundamental fluid-structure interaction (FSI) problem, which contains complex nonlinear flow and structural dynamics, including flow separations, vorticity generations and vortex-induced vibrations. The development of an efficient method to control the flapping behaviour of an elastic beam in viscous flow is of great importance for engineering applications.
The prediction of the dynamics of elastic beam-like structures subjected to unsteady flow relies on accurate FSI models. Numerical solutions to FSI problems of a beam in uniform flow can be obtained using either a monolithic or a partitioned approach. The monolithic approach (Heil, Hazel & Boyle Reference Heil, Hazel and Boyle2008; Ryzhakov et al. Reference Ryzhakov, Rossi, Idelsohn and Oñate2010; Pfister, Marquet & Carini Reference Pfister, Marquet and Carini2019) simultaneously solves the governing equations of the fluid and the structure as a fully coupled system. However, this method requires high computational costs when large and complex systems equations are involved. The partitioned approach integrates a fluid solver and a structure solver by data transfer at the common fluid-structure interface. For a loosely coupled partitioned method (Degroote et al. Reference Degroote, Bruggeman, Haelterman and Vierendeels2008), only unidirectional data transfer of the structural and fluid domains is considered. In contrast, a strongly coupled partitioned method (Liu et al. Reference Liu, Qu, Xie and Meng2022, Reference Liu, Qu, Xie and Meng2023) accounts for both the flow excitation acting on the structure and the influence of the deformation and motion of the structure on the surrounding flow field. Multiple sub-iterations per time step in a strongly coupled partitioned method enable the added mass effect and can accurately analyse the FSIs. These characteristics make the strongly coupled partitioned method particularly suitable for simulating systems with bidirectional electromechanical coupling, which is a critical requirement for our fluid-structure-electro model. The arbitrary Lagrangian–Eulerian (ALE) framework is further selected over conventional Lagrangian approaches because it allows simultaneous tracking of large beam deformations (via mesh motion) and vortex shedding patterns (via Eulerian flow description). This method effectively balances computational accuracy and efficiency in tracking the fluid-structure interface motions (Souli & Benson Reference Souli and Benson2013).
The flapping dynamics of a slender structure have been extensively studied based on strongly coupled FSI models. Earlier studies (Paidoussis Reference Paidoussis1966; Coene Reference Coene1992) focused on linear vibration responses of elastic beams in an inviscid fluid. Later, the effects of the flow vortices, viscosity and the structural nonlinearity on the flapping dynamics of elastic beams were investigated (Moretti Reference Moretti2003; Argentina & Mahadevan Reference Argentina and Mahadevan2005). Zhang et al. (Reference Zhang, Childress, Libchaber and Shelley2000) identified two flapping modes of a thin flag, i.e. the stretched-straight state and the flapping state. Connell & Yue (Reference Connell and Yue2007) found more dynamic modes of flapping flags, including period-one limit-cycle flapping and chaotic flapping mode. Shelley & Zhang (Reference Shelley and Zhang2011) presented a comprehensive review of the interactions between fluid flows and flapping bodies. The addition of a circular cylinder upstream of a flapping beam introduces complexities of vortex-induced dynamics. Shukla, Govardhan & Arakeri (Reference Shukla, Govardhan and Arakeri2013) found two dynamic modes of a flexible beam for different Reynolds numbers. Lee & You (Reference Lee and You2013) examined the effects of the beam length and stiffness on the flow-induced vibration modes and the associated wake patterns. Recently, Pfister & Marquet (Reference Pfister and Marquet2020) investigated typical dynamic modes of a hyperelastic beam, including the steady symmetric mode, the symmetric and periodic mode, the deviated and periodic mode and the symmetric and quasi-periodic mode. Furquan & Mittal (Reference Furquan and Mittal2021) analysed the dynamic responses of flexible beam-like structures with different flexibility and inertia, and four typical vibration modes of the structure and corresponding wake patterns were identified.
Fluid-structure interaction of a slender structure not only induces flapping motions but also leads to spontaneously broken symmetry of the structure (Alben & Shelley Reference Alben and Shelley2005). Phenomena related to symmetry-breaking (S-B) motions of elastic beams subject to unsteady flow were studied by Cimbala & Chen (Reference Cimbala and Chen1994), Cimbala, Garg & Park (Reference Cimbala, Garg and Park1988), Cimbala & Garg (Reference Cimbala and Garg1991). They found that a rigid beam attached to a rotating circular cylinder in uniform flow deflects to one side in some circumstances. The symmetry breaking of the cylinder-beam system was found to be related to the length of the beam. Furthermore, the Reynolds number is also found to be the critical parameter that determines the S-B bifurcation (Xu, Sen & Gad-el-Hak Reference Xu, Sen and Gad-el-Hak1990; Xu, Sen & Gad-el-Hak Reference Xu, Sen and Gad-el-Hak1993). In recent years, Lācis et al. (Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014), Bagheri, Mazzino & Bottaro (Reference Bagheri, Mazzino and Bottaro2012) reviewed previous studies and demonstrated that the S-B instability of rigid and flexible beams attached to circular cylinders resembles the instability of an inverted pendulum, which relates to the length of the beam and the area of the back-flow region. Pfister & Marquet (Reference Pfister and Marquet2020) employed a linear fully coupled fluid-structure analysis method (Pfister et al. Reference Pfister, Marquet and Carini2019) for analysing the symmetry breaking of a hyperelastic beam attached to a fixed circular cylinder and established a relationship between the instability and the eigenvalue of the system.
The control of flapping motions of slender structures subjected to viscous flow is of great importance in engineering. Traditional methods encompass both passive control and active control approaches. For the passive control methods, the tuned mass damper (TMD) and its variants (Igusa & Xu Reference Igusa and Xu1994; Zuo Reference Zuo2009; El-Khoury & Adeli Reference El-Khoury and Adeli2013; Elias & Matsagar Reference Elias and Matsagar2017) were widely used due to their simplicity. To achieve broadband vibration suppression, nonlinear energy sinks were introduced (Ahmadabadi & Khadem Reference Ahmadabadi and Khadem2012; Kani et al. Reference Kani, Khadem, Pashaei and Dardel2016; Ding & Chen Reference Ding and Chen2020). Passive vibration control methods can not suppress the flapping motions of elastic structures adaptively. Active vibration control methods were applied to suppress vibrations for complex excitations. The widely used ‘sensor-actuator’ configuration, consisting of piezoelectric patches positioned on the structure, forms a closed-loop control system. Based on the sensed strain, the closed-loop control system achieves vibration control by introducing various control algorithms (Aktas & Esen Reference Aktas and Esen2020; Vladimír et al. Reference Vladimír, Juraj, Gálik and Justín2021). Recently, metamaterials based on periodic embedded structures have been widely applied in vibration control. Piezoelectric metamaterials can effectively achieve specific vibration bandgaps (Liu et al. Reference Liu, Zhang, Mao, Zhu, Yang, Chan and Sheng2000; Zhu et al. Reference Zhu, Chen, Zhu, Garcia-Vidal, Yin, Zhang and Zhang2013; Sugino et al. Reference Sugino, Leadenham, Ruzzene and Erturk2017; Yuksel & Yilmaz Reference Yuksel and Yilmaz2020; Prasad & Banerjee Reference Prasad and Banerjee2021). While existing studies related to piezoelectric metamaterials on vibration control of FSI systems are for structures in a steady fluid environment (Dai et al. Reference Dai, Zheng, Mao and Qu2023; Guo et al. Reference Guo, Sheng, Zeng, Wang and Li2023; Zheng et al. Reference Zheng, Qu, Dai, Chen and Mao2024), their fixed bandgap properties fail to adapt to time-varying vortex shedding frequencies. This limitation motivates our development of adaptive metamaterials for the challenge.
While previous investigations have thoroughly characterised the dynamic behaviours and bifurcation mechanisms, two critical gaps remain: (i) existing vibration control strategies predominantly target linear conditions and lack adaptability to strongly nonlinear flow environments, (ii) no study has integrated piezoelectric metamaterials for effective bifurcation control and symmetry modulation. To bridge these gaps, this paper presents a piezoelectric metamaterial design method for controlling the flapping and S-B bifurcation motions of an elastic beam behind a circular cylinder in viscous flow. The beam is equipped with a series of smart piezoelectric meta-cells connected by unidirectional circuits. The Euler–Bernoulli theory combined with the von Kármán strains is employed to describe the geometrical nonlinearity of the beam. An implicit and strongly coupled fluid-structure-electro-control numerical model of the piezoelectric metamaterial beam is developed in an arbitrary ALE framework. Beyond suppressing bifurcations, the proposed strategy enables adaptive modulation of flapping modes and wake patterns through self-sensing and adaptive-actuation perturbations. Floquet theory and flow field analyses are employed to reveal the stability boundaries and physical mechanisms underlying the observed transitions. In particular, the system reorganises fluid-structure energy exchange pathways, as revealed by wake recovery and energy harvesting metrics. This provides a fluid-mechanical basis for tunable energy conversion and propulsion under unsteady flow conditions. The mechanisms of flapping control, S-B modulation and their potential engineering applications are discussed in detail.
2. Theory and numerical model
2.1. Design model
The schematic of a piezoelectric metamaterial beam attached to a stationary circular cylinder and subjected to uniform viscous flow is given in figure 1(a). Here
$U_{\infty }$
is the magnitude of the incoming flow velocity. Flow separation occurs after the cylinder and leads to vorticity generation on the upper and lower sides of the beam, which induces nonlinear flapping of the beam. To modulate the flow-induced vibration, the surface of the beam is uniformly covered by a total number of
$N_{p}$
piezoelectric patch pairs, which are connected by unidirectional amplification circuits. Using the piezoelectric effect, the piezoelectric patches obtain the sensing voltage according to the strain of the beam. Subsequently, based on the inverse piezoelectric effect and an amplification circuit, the sensed voltage signal is converted into an amplified actuation signal. A pair of two patches and the corresponding base beam part is recognised as a meta-cell, which is the fundamental unit of the metamaterial, as shown in figure 1(b).
Schematic of the metamaterial beam clamped at a circular cylinder subject to a uniform flow.
The diameter of the circular cylinder is
$D$
, the length of the beam is
$L_{b}$
, and the thickness of the beam is
$H_{b}$
. The surface of the beam is covered by
$N_p$
piezoelectric patch pairs. The patches are distributed uniformly along the beam. Every two pairs of patches and the corresponding base beam part are recognised as a meta-cell, which is the fundamental unit of the metamaterial (the number of the meta-cells is
$N_{p}/2$
). As shown in figure 1(c), the first piezoelectric patch is the sensor unit, while the second patch is the actuator unit. Two patches have the same length and thickness represented by
$L_{p}$
and
$H_{p}$
, respectively. A unidirectional voltage amplification circuit is established with an inverting amplifier between the sensor circuit and the actuator circuit. All patches are grounded. Here
$\phi$
is defined as the voltage of every pair of piezoelectric patches. The sensing unit feels the local strain and adaptively operates through the gain amplification circuit.
As shown in figure 1(c), an amplification circuit with a gain factor
$G$
is provided for voltage amplification from the sensor unit to the actuator unit. Figure 2 shows two amplification configurations with a phase-reversal type (
$G\lt 0$
) and in-phase type (
$G\gt 0$
). In the amplification circuit, also known as a voltage-controlled voltage source circuit, the resistance
$R_{G}$
exceeds the input resistance
$R_{I}$
significantly and the
$R_{F}$
is designed to minimise current flow across itself for stability.
Conceptual illustration of amplification circuit configuration: (a) phase-reversal amplification and (b) in-phase amplification.
The voltage relationship between the sensor unit and the actuator unit is given by
\begin{equation} \phi _{{act}}=G\phi _{\textit{sen}}, \end{equation}
where the amplification feedback factor
$G$
is defined by
$G=-R_{F}/R_{I}$
for the phase-reversal amplification circuit and
$G=1+R_{F}/R_{I}$
for the in-phase amplification circuits.
2.2. Fluid model
The incompressible Navier–Stokes equations under an ALE framework are employed to describe the unsteady viscous flow by the motion and deformation of the piezoelectric metamaterial cantilever beam, given by
\begin{align} \frac{\partial v_{i}}{\partial x_{i}} &= 0,\nonumber\\[4pt] \frac{\partial v_{i}}{\partial t}+ (v_{j}-v_{m,j})\frac{\partial v_{i}}{\partial x_{j}} &= -\frac{1}{\rho _{\!f}}\frac{\partial p_{\!f}}{\partial x_{i}}+\frac{{\rm \mu} _{\!f}}{\rho _{\!f}}\frac{\partial ^{2}v_{i}}{\partial x_{j}\partial x_{i}}, \end{align}
where
$v_{i}$
is the velocity component of the fluid,
$v_{m,i}$
is the mesh velocity component in the ALE framework,
$p_{\!f}$
is the fluid pressure,
$\rho _{\!f}$
is the density of the fluid and
${\rm \mu} _{\!f}$
is the dynamic viscosity of the fluid. The ALE formulation naturally incorporates the mesh velocity
$v_{m,i}$
into the convective terms of the Navier–Stokes equations, ensuring Galilean invariance and consistent momentum transport across the moving interface. The ALE method is employed to handle the strong interaction between the moving and deforming structure and the surrounding viscous fluid. The ALE framework combines the advantages of both Lagrangian and Eulerian descriptions by allowing the computational mesh to move independently of both the fluid and the structure, enabling the accurate capture of large deformations without excessive mesh distortion. The motion of the computational mesh is governed by the Laplace smoothing equation, which ensures smooth mesh deformation, especially near the beam surface, and preserves mesh quality in highly dynamic regions. The finite-volume method is employed to discretise (2.2) on an unstructured mesh. Equation (2.2) can be rewritten with a general scalar variable
$\psi$
, i.e.
\begin{equation} \frac{\partial \psi }{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left[\left(v_{i}-v_{m.i}\right)\psi \right]=\boldsymbol{\nabla }\boldsymbol{\cdot }\left({\Gamma} \boldsymbol{\nabla }\psi \right)+Q, \end{equation}
where
${\Gamma}$
is the diffusion coefficient of the general scalar variable
$\psi$
and
$Q$
is the source term of the incompressible Navier–Stokes equations. The fluid domain discretised on the element
$C$
from
$t$
to
$t+{\unicode[Arial]{x0394}} t$
can be described by the discretised integrated form of (2.3),
\begin{align} \int _{t}^{t+{\unicode[Arial]{x0394}} t}\int _{{V_{C}}}\frac{\partial \psi }{\partial t}{\textrm{d}}V\!{\textrm{d}}t &+\, \int _{t}^{t+{\unicode[Arial]{x0394}} t}\left[\sum _{f^{*}\sim nb(C)}\left(\, \int \limits_{f^{*}}(v_{i}-v_{m,i})\psi \boldsymbol{\cdot }{\textrm{d}}S\right)\right]{\textrm{d}}t\nonumber\\[4pt] & -\,\int _{t}^{t+{\unicode[Arial]{x0394}} t}\left[\sum _{f^{*}\sim nb(C)}\left(\, \int \limits_{f^{*}}{\Gamma} \boldsymbol{\nabla }\psi \boldsymbol{\cdot }{\textrm{d}}S\right)\right]{\textrm{d}}t\nonumber\\[4pt] & = \int _{t}^{t+{\unicode[Arial]{x0394}} t}\left[\sum _{f^{*}\sim nb(C)}\left(\int \limits_{V_{C}}Q{\textrm{d}}V\right)\right]{\textrm{d}}t, \end{align}
where
$f^{*} , V_{C}$
and
$S$
are the faces of the discretised element
$C$
, the element volume and the surface vector, respectively. The subscript
$f^{*}\sim nb(C)$
denotes all the other faces of the element
$C$
next to the neighbouring elements. To compute the flow response, a transient solver for incompressible flow called pimpleFoam in OpenFOAM (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998) is employed for solving (2.4). The transient term is executed using the second-order implicit Euler method, the convective term is treated by the linear upwind difference method and the diffusion term is treated by the central difference method. As a pressure-based solver, pimpleFoam utilises a PIMPLE (pressure implicit with the splitting of the operator) algorithm, the combination of SIMPLE (semi-implicit method for pressure-linked equations) and PISO (pressure-implicit with splitting of operators) algorithms for the velocity field correction and the pressure field determination. The realisation of the velocity-pressure correction algorithm depends on constructing an inner corrector and an outer corrector. The inner corrector determines the number of pressure corrections within each iteration, while the outer corrector defines the number of iterations used to solve the governing equations per time step. The corrected pressure values obtained from the last PIMPLE loop are used.
2.3. Structural model
The deformation, motion and the electromechanical coupling effect of the piezoelectric metamaterial beam with a unidirectional amplification circuit can be described by
\begin{align} \left[m_{\textit{eff}}^{b}+m_{\textit{eff}}^{p} \unicode{x1D6F7} (x)\right]\frac{\partial ^{2}w}{\partial t^{2}}+\left[D_{\textit{eff}}^{b}+D_{\textit{eff}}^{p}\unicode{x1D6F7} (x)\right]\frac{\partial ^{4}w}{\partial x^{4}}+e_{\textit{eff}}\frac{\phi }{H_{p}}\unicode{x1D6F7} (x) &= q_{\!f}(x,t),\nonumber\\[3pt] e_{\textit{eff}}\frac{\partial ^{2}w}{\partial x^{2}}+g_{p}\frac{A_{p}}{H_{p}}\phi _{\textit{sen}} &= 0, \end{align}
where
$w(x,t)$
is the deflection of the piezoelectric beam representing the vibration response and
$\phi$
is the electric potential of the piezoelectric patches (sensors and actuators), while
$\phi _{s{en}}$
specifically refers to the electric potential on the sensor unit. The relationship between
$\phi _{s{en}}$
and
$\phi _{{act}}$
is presented in (2.1). Here
$q_{\!f}(x,t)$
is the fluid forces acting on the beam,
$m_{\textit{eff}}^{b}$
and
$m_{\textit{eff}}^{p}$
are the effective density of the beam and the piezoelectric part,
$D_{\textit{eff}}^{b}$
and
$D_{\textit{eff}}^{p}$
are the effective bending stiffness of the beam and the piezoelectric part,
$e_{\textit{eff}}$
denote the effective piezoelectric material constant,
$g_{p} , A_{p}$
and
$H_{p}$
represent the dielectric material constant, area, and the height of each piezoelectric patch, respectively, and
$\unicode{x1D6F7} (x)$
is the distribution function that describes the arrangement of the periodic piezoelectric patches on the beam.
A two-node finite element method (FEM) is employed for predicting the numerical solution of the (2.5). The dynamic equation in finite element form can be written as
\begin{equation} \boldsymbol{M}\ddot{\boldsymbol{\xi }}+\boldsymbol{C}\dot{\boldsymbol{\xi }} + (\boldsymbol{K}_{L}+\boldsymbol{K}_{\textit{NL}}+\boldsymbol{K}_{G})\boldsymbol{\xi }=\boldsymbol{F}_{\!\!f}, \end{equation}
where
$\boldsymbol{\xi }$
is the generalised displacement vector,
$\boldsymbol{M}$
and
$\boldsymbol{C}$
are the mass matrix and the damping matrix of the metamaterial beam,
$\boldsymbol{K}_{L} , \boldsymbol{K}_{\textit{NL}}$
and
$\boldsymbol{K}_{G}$
are the linear and nonlinear stiffness matrix and the additional electromechanical coupling stiffness matrix of the metamaterial beam, respectively, and
$\boldsymbol{F}_{\!\!f}$
is the fluid force vector acting on the beam. More details about the derivation of the dynamic equation and the finite element discretisation process of the piezoelectric beam are described in Appendix A. Equation (2.6) is a standard nonlinear vibration equation for FEM. The implicit Newmark-beta method combined with the Newton–Raphson iterative method is used to obtain the fluid-induced nonlinear vibration response of the piezoelectric metamaterial beam.
2.4. Fluid-structure interaction model
An implicit strongly coupled partitioned algorithm is employed for the analysis of the FSI problem. The fluid and structural solvers iterate sequentially until convergence at the fluid-structure coupling interface every time step. The boundary conditions on the interface are given by
\begin{align} -p_{\!f}\boldsymbol{n}+\boldsymbol{\tau }_{\!f}\boldsymbol{\cdot} \boldsymbol{n}&=\sigma _{s}\boldsymbol{n}, \nonumber\\ \boldsymbol{v}_{\!f}&=\dot{\boldsymbol{\xi }}, \end{align}
where
$\boldsymbol{n}$
is the unit normal vector on the fluid-structure interface,
$\boldsymbol{\tau }_{\!f}$
and
$\boldsymbol{v}_{\!f}$
are the shear stress and the flow velocity at the interface, respectively, and
$\sigma _{s}$
is the stress of the metamaterial beam on the interface. During each iteration, the fluid solver is first conducted to compute the velocity and pressure of the fluid. The fluid force is then obtained and interpolated to the finite element mesh of the beam for determining the nonlinear structural response by the structural solver. The vibration response of the beam is interpolated to the finite-volume mesh of the fluid. The fluid mesh is responsible for regulating the internal motion of the mesh nodes, adapting to both the movement of the fluid-structure interface and the deformation of the internal mesh.
The interaction between the fluid and the structural solvers demands a precise data-mapping strategy for the non-matching meshes of the fluid and the beam. The radial basis functions method (Joshi, Jaiman & Ollivier-Gooch Reference Joshi, Jaiman and Ollivier-Gooch2020) with a thin plate splines basic function is used for data mapping between the beam and the fluid. Based on the accurate data-mapping method, a predictor-corrector method is utilised to realise the strongly coupled partitioned algorithm for FSI. An inner iteration with under-relaxation is used in the corrector step to comply with the interface coupling conditions. The partitioned approach allows the fluid and structural solvers to run alternately until convergence conditions are satisfied. The details about the algorithm are shown in table 1. The interaction between the fluid and the structure is mainly governed by three non-dimensional coefficients, the Reynolds number
${Re}=\rho _{\!f}U_{\infty }D/{\rm \mu} _{\!f}$
, the mass ratio
$M_{r}=\rho _{s}/\rho _{\!f}$
and the reduced velocity
$U_{r}=U_{\infty }/f_{n}D$
, where
$f_{n}$
is the first natural frequency of the piezoelectric cantilevered beam. The numerical validation of the proposed method and the mesh convergence study are presented in Appendix B.
Computational algorithm for FSI of piezoelectric metamaterial beam.
3. Results and discussion
3.1. Flapping behaviours of the beam
Four flapping modes of the piezoelectric metamaterial beam are identified within the parameters examined in this investigation, which are known as the limit-cycle oscillation (LCO) mode, the quasi-periodic oscillation mode, the subharmonic oscillation mode and the S-B mode. Influenced by Reynolds number and reduced velocity, the piezoelectric metamaterial beam switches spontaneously between the above modes due to lock-in phenomena (Connell & Yue Reference Connell and Yue2007) and fluid-structure instability (Lācis et al. Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014).
For
$M_{r}=10.0$
and
${Re}=400$
, when
$U_{r}$
is low, the shedding vortex frequency is away from the natural frequency of the beam. The vibration response of the beam settles to a period-one LCO mode. The phase plot of the beam turns into a limit-cycle trajectory as shown in figure 3(a). Besides, the shedding vortex of the beam follows the wake pattern of a von Kármán street, which consists of a continuous series of alternately signed vortices as shown in figure 3(c). With the increase of
$U_{r}$
, the shedding vortex frequency matches the natural frequency of the beam. The resonance overlap induces large-amplitude vibration of the beam and a transition to complexity and nonlinearity from the LCO mode (Sarkpaya Reference Sarkpaya1979; YS Reference Ys1991; Williamson & Govardhan Reference Williamson and Govardhan2004; Connell & Yue Reference Connell and Yue2007). Besides the dominant frequency component
$f_{d}$
, a 1/2 of the dominant frequency component
$f_{d}/2$
is also observed. Figure 3(b) demonstrates a periodic-doubling bifurcation from the LCO mode to the subharmonic mode of the piezoelectric metamaterial beam. Furthermore, the phase portrait of the subharmonic mode exhibits closed-loop trajectories, which indicate that the system’s state repeats after a specific period that is a multiple of the driving period. Two loops are observed in figure 3(a) for the subharmonic oscillation of order 1/2. Due to the lack of symmetry, the vortex wake of the subharmonic mode of the beam becomes irregular as shown in figure 3(c). When the natural frequency of the beam moves away from the vortex-shedding frequency, the LCO mode of the beam is observed again at higher
$U_{r}$
. The LCO mode is dominated by the second-order bending mode rather than the first-order bending mode of the cantilever beam.
Flapping dynamics of the metamaterials beam attached to a fixed circular cylinder in uniform flow with
${Re}=400$
. (a) Vibration response of the free tip, (b) dominated frequency of the beam, (c) shedding vortex of different modes.
Figure 4 demonstrates the vibration response, frequency response and shedding vortex patterns of the piezoelectric metamaterial beam in viscous flow with
${Re}=200$
. When the shedding vortex frequency approaches the natural frequency of the beam around
$U_{r}=5.0$
, two vibration modes compete with each other and turn out into a quasi-periodic oscillation mode. The phase portrait in figure 4(a) indicates that trajectories densely fill a torus in phase space. The characteristics of quasi-periodic lead to a modulation of the shedding vortex as shown in figure 4(c). The shedding vortex no longer propagates in a straight line downstream but gradually deviates from its equilibrium position after separation from the beam. With the increase of
$U_{r}$
, the interaction between the elasticity of the beam and the detached flow behind the circular cylinder leads the system to a defecting oscillation mode through a S-B bifurcation. Lācis et al. (Reference Lācis, Brosse, Ingremeau, Mazzino, Lundell, Kellay and Bagheri2014) revealed that the torque inequality of the beam inside and outside the backflow region generates an unstable equilibrium point and induces a S-B oscillation of the beam with a deflecting beam toward the bottom or the top spontaneously. As shown in figure 4(a), the beam vibrates for a LCO away from the equilibrium position. The phase plot is asymmetrical, and the shedding vortex maintains a smaller intensity than other flapping modes.
Flapping dynamics of the metamaterials beam attached to a fixed circular cylinder in uniform flow with
${Re}=200$
. (a) Vibration response of the free tip, (b) dominated frequency of the beam, (c) shedding vortex of different modes.
3.2. Flapping bifurcation suppression and stability modulation
3.2.1. Bifurcation suppression and mode transition
The flow-induced vibration responses of the piezoelectric metamaterial beam equipped with unidirectional circuits are examined. The following parameters are considered for the flow:
${Re}=400$
and
$M_{r}=10.0$
. The configuration parameters of the piezoelectric metamaterial beam follow the same setting as the case in the mesh convergence test in Appendix B and figure 24. The displacement and frequency response of the beam with
$G=120$
is shown in figure 5. With small
$U_{r}$
, the displacement amplitude and the frequency response remain the same with and without control. As
$U_{r}$
increases, the circuit significantly changes the flapping modes of the beam. The beam does not undergo a sudden change to the subharmonic mode at the original critical bifurcation reduced velocity around
$U_{r}=8.0$
. In the case of
$U_{r}=10.0$
(A) and
$U_{r}=13.0$
(B), the amplitude of the beam is suppressed
$57.83\,\%$
and
$25.61 \,\%$
after turning on the control circuit. The critical reduced velocity is delayed from
$U_{r}=8.0$
to
$U_{r}=12.0$
. The circuit brings a wider range of the little-amplitude LCO mode. As shown in figure 5(b), the flapping frequency of the beam does not jump to 1/2 of the dominant limit-cycle frequency component
$f_{d}/2$
anymore when the natural frequency of the beam approaches the vortex shedding frequency. The sudden change of the beam occurs around
$U_{r}=12.0$
. Unlike the periodic-doubling bifurcation, the beam does not enter a subharmonic mode of vibration with large amplitude, but rather a single-frequency dominated mode of vibration with high-frequency components after the sudden change, which is shown in figure 5(c) when
$U_{r}=18.0$
(C). An increase in the gain factor leads to a shift in the dominant flapping frequency, and the LCO frequency spectrum is broadened. This transition indicates that the piezoelectric metamaterial effectively disrupts the natural resonant frequency of the FSI system, redistributing the vibrational energy across multiple frequency components and thereby mitigating large-amplitude oscillations.
Flapping control effect of the piezoelectric metamaterial with gain factor
$G=120$
: (a) variation of displacement against reduced velocity, (b) variation of frequency against reduced velocity, (c) vibration response of beam with ‘control off’ and ‘control on’ at
$U_{r}=10.0$
(A),
$U_{r}=13.0$
(B) and
$U_{r}=18.0$
(C).
(a) Bifurcation diagrams of the piezoelectric metamaterial beam (
$G=120$
). (b) Tip-velocity displacement phase portraits of the beam for
$U_{r}=7.0$
(III),
$U_{r}=12.0$
(IV),
$U_{r}=15.0$
(V) and
$U_{r}=18.0$
(VI).
Figure 6 illustrates the vibration mode evolution of the beam with
$G=120$
. As
$U_{r}$
increases, the beam first transforms from a very small-amplitude LCO mode (I) to a medium amplitude LCO mode (III) through a quasi-periodic transition region (II). Different from the bifurcation of the uncontrolled beam, during the sudden change, the controlled beam undergoes another quasi-periodic transition (IV) for mode switching between the stable LCO mode (III) and the new flapping modes (V) with cloud-like (CL mode) phase plots as shown in figure 6(b). The transition occurs as a result of a decrease in beam stiffness and the additional stiffness effect of the metamaterials. In the case of a further increase of the reduced velocity, high-frequency components are introduced into the modes of the beam. The CL mode (V) dominated by a single frequency, gradually evolves into the quasi-CL mode (VI). The essence of the sudden change comes from the mode competition between the flow-induced vibration and the piezoelectric metamaterial-induced mode. Periodic piezoelectric patches with ‘sensor-actuator’ meta-cells bring additional frequency components through adaptive actuation locally. The competition between multiple frequency components leads to the quasi-periodic transition. With the increase of
$U_{r}$
, the beam becomes more flexible and the piezoelectric metamaterial-induced mode dominates the vibration. The vibration of the beam shifts from a flow-induced-vibration-dominated LCO mode (III) to a metamaterial-dominated CL mode (V) through the quasi-periodic transition. The CL mode in region V is dominated by the actuation of the piezoelectric metamaterial. Multiple loops are observed in the non-dimensional displacement-velocity phase plot (
$\xi ^{*}=w/h , \dot{\xi }^{*}=\dot{w}/h$
), which tends to occur at the maximum vibration amplitude due to the locally adaptive actuation of the beam, as shown in figure 6(b). When the beam undergoes continuous bending to larger amplitudes, the sensor units adaptively sense the stresses, prompting the actuator units to generate corresponding deformation at moments through the circuit. This causes rapid inverse changes in beam displacement and velocity. The high-frequency fluctuations at the maximum of the displacements in the vibration response illustrated in figures 6(b) and 6(c) can be seen to correspond to the multiple loops in the phase plot figure 6(b).
A deeper investigation into the period-doubling bifurcation behaviour and the evolution of its dynamic modes following modulation hinges on the stability analysis of the system. The Floquet theory (Barkley & Henderson Reference Barkley and Henderson1996; Pernot & Lamarque Reference Pernot and Lamarque2001) is employed to analyse the stability of the nonlinear vibration response of the piezoelectric beam. According to the Floquet theory, the state transition relationship of a periodic system that satisfies
$\dot{x}=\boldsymbol{A}(t)x$
and
$\boldsymbol{A}(t+T)=\boldsymbol{A}(t)$
can be described by
$\dot{\boldsymbol{\unicode{x1D6F7} }}(t)=\boldsymbol{A}(t)\boldsymbol{\unicode{x1D6F7} }(t)$
. The matrix
$\boldsymbol{\unicode{x1D6F7} }(t)$
is defined as the monodromy matrix. The Floquet multipliers are obtained by solving the eigenvalue problem of the monodromy matrix, expressed as
$\boldsymbol{\unicode{x1D6F7} }(T)\boldsymbol{v}_{i}={\rm \mu} _{i}\boldsymbol{v}_{i}$
. The absolute value of the Floquet multiplier characterises the stability of the system: when
$| {\rm \mu} _{i}| \gt 1$
, the system exhibits unstable disturbance growth. When
$| {\rm \mu} _{i}| \lt 1$
, it indicates stable disturbance attenuation. Additionally, the Floquet exponent is defined as
$\lambda _{i}=1/T\ln {\rm \mu} _{i}$
, where the real part of
$\lambda _{i}$
represents the exponential growth rate. It serves as an indicator of system stability as well. Specifically, when
${\textrm{Re}}(\lambda _{i})\gt 0$
, the system is unstable; when
${\textrm{Re}}(\lambda _{i})\lt 0$
, the system is in a stable state. In this study, Floquet stability analysis is conducted to examine the nonlinear dynamic response at the cantilever beam’s tip. The phase space is constructed by combining displacement-time and velocity-acceleration state vectors. The dynamic matrix of the system is derived via least squares estimation, while the principal period is determined through autocorrelation function analysis. Utilising the dynamic matrix and identified period, the monodromy matrix can be formulated. Subsequently, eigenvalue analysis is performed to assess the Floquet stability of the nonlinear vibration response of the beam.
Figure 7 shows the Floquet multipliers evolution and the Floquet exponents distribution in the complex space of the nonlinear vibration response of the piezoelectric beam against the reduced velocities. When the system operates in the LCO mode without control, its dynamic response exhibits a stable characteristic. However, due to the strong fluid-structure coupling effect, the system transitions into the subharmonic mode near
$U_{r}=7.0{\sim}15.0$
. Consequently, the Floquet multiplier immediately exceeds 1, leading the system into an unstable state. The period-doubling bifurcation, a mechanism that induces chaos, further elevates the Floquet multiplier, thereby destabilising the system. Nevertheless, at higher reduced velocities, the system enters the second-order LCO mode and regains stability. Figure 7(b) demonstrates that, in the absence of control, the stable modes of the system are predominantly distributed at both ends of the studied reduced velocity interval. In contrast, the system response controlled by the piezoelectric metamaterial exhibits enhanced stability. Within the range of the investigated parameters, the Floquet multipliers of most states remain below 1 (I, III and V), whereas the Floquet multipliers exceed 1 in intervals associated with mode transitions II, IV and VI. This instability corresponds to quasi-periodic modes during mode transitions. Instability caused by mode competition diminishes more rapidly as the dominance of a particular mode increases. Therefore, compared with the uncontrolled state, the controlled system demonstrates significantly improved stability. This indicates that within a broader range of reduced velocities, the controlled system avoids period-doubling bifurcation and sustains robust stability.
Stability and bifurcation analysis of the piezoelectric metamaterial beam flapping based on the Floquet theory: (a,c) the Floquet multipliers, and (b,d) the Floquet exponents distribution in complex space against the reduced velocities of the beam nonlinear vibration response when ‘control off’ and ‘control on’.
3.2.2. Wake-beam interaction and energy transfer pathway
To understand the physical mechanisms underlying the observed dynamic regimes, the analysis concerning the flow field structures, vortex dynamics and energy exchange characteristics across different flapping states is presented. Figure 8 shows the deformation envelope and tip displacement of the piezoelectric metamaterial beam under three reduced velocities:
$U_{r}=10.0 , U_{r}=13.0$
and
$U_{r}=18.0$
, with a gain factor
$G=120$
. The unidirectional circuit suppresses the subharmonic vibration and turns the beam into quasi-periodic vibration when
$U_{r}=10.0$
. This yields a stable, regular 2S wake (see figure 8
a) rather than irregular wakes in figure 8(c). The beam exhibits a CL flapping mode with a large amplitude and localised high-frequency perturbations near displacement extrema at higher reduced velocities (
$U_{r}=13.0$
). These perturbations introduce rapid inversions in displacement and velocity, enhancing vortex stretching near the beam tip. Simultaneously, local mid-span deformation disrupts wake coherence, causing vortex splitting and resulting in a symmetric 2P shedding pattern (see figure 8
b). When
$U_{r}=18.0$
, the rapid inverse changes of displacement and velocity due to the adaptive actuation mechanism with self-sensing capability are significant. More high-frequency fluctuation components and a more symmetrical 2P mode is observed. Paired vortices with opposite directions gradually move away from the axis of symmetry. The observed transition from the 2S to 2P mode reflects a significant restructuring of the FSI dynamics. As
$U_{r}$
increases, stronger beam deformation and spatially distributed actuation enhance local flow acceleration and mid-span perturbations. These effects, compounded by the rebound effect, lead to multiple vortex shedding sites and a higher-order 2P topology. The emergence of the 2P mode thus results from the coupled influence of structural nonlinearity, adaptive electromechanical forcing and vortex pairing dynamics.
Flapping control of the metamaterial beam (
$G=120$
): the deformation envelope of the beam (first row), the transversal displacement at the beam tip (second row), and the vorticity contour around the beam (third row) when (a)
$U_{r}=10.0$
, (b)
$U_{r}=13.0$
and (c)
$U_{r}=18.0$
.
The dynamic mode decomposition (DMD) method, a data-driven technique for extracting coherent structures with specific frequency components (see Appendix C for details) is conducted to analyse the flow structures behind the beam. Here
$N=800$
snapshots of the wake behind the cylinder-beam system are selected from
$t=0.2\, {\textrm{s}}$
to
$t=0.6\, {\textrm{s}}$
with
${\unicode[Arial]{x0394}} t=0.0005\, {\textrm{s}}$
for the DMD data matrix assemble. The first m-order (
$m=25$
) modes are truncated for singular value decomposition (SVD) operation. The frequency spectra of the beam displacement responses and the corresponding DMD modes of the wake are shown in figure 9. For
$U_{r}=10.0$
, the vibration of the beam is mainly dominated by the frequency components
$f_{1} , f_{2}$
and
$3f_{1}$
. The frequency component
$f_{1}$
is mainly due to the FSI of the beam, and the frequency component
$3f_{1}$
is produced by the geometric nonlinearity of the beam. The DMD modes of the flow corresponding to these two frequency components are associated with periodical and symmetrical vortices shed from the tip of the beam. The frequency component
$f_{2}$
introduces a perturbation to the original system and presents a mode competition with the FSI mode (
$f_{1}$
and
$3f_{1}$
). It makes it difficult to enter the periodic-doubling bifurcation and turns the system into a quasi-periodic mode with low amplitude. Figure 9(b) shows a single-frequency dominated vibration and wake pattern at
$U_{r}=13.0$
. The vibration is mainly dominated by
$f_{1}$
and odd-order frequency components
$3f_{1}$
and
$5f_{1}$
brought by geometric nonlinearity. The distribution of vortices with a clear discontinuous region in DMD modes between the centre region and the upper and lower regions is observed. This discontinuity comes from the rapid inverse changes of displacement and velocity at the beam tip. Figure 9(c) illustrates a symmetrical single-frequency dominated vibration and wake mode at
$U_{r}=18.0$
. The vibration response of the beam contains the dominant frequency
$f_{1}$
, odd-order frequency component
$3f_{1}$
and component near the dominant frequency
$f_{2}$
. A three-part flow structure can be observed as well due to the rapid rebound at the upper and lower regions.
Spectra of displacement responses of the piezoelectric metamaterial beam (first column) and DMD of wake when (a)
$U_{r}=10.0$
, (b)
$U_{r}=13.0$
and (c)
$U_{r}=18.0$
.
Figure 10 exhibits the contours of the instantaneous fluid pressure of the piezoelectric metamaterial beam. Here
$t_{1} , t_{2} , t_{3}$
and
$t_{4}$
denote instantaneous moments uniformly distributed over a period
$T$
. In the case of
$U_{r}=10.0$
, a distinct low-pressure centre is generated at the mid-rear position of the beam, where two oppositely directed vortices are induced and move downstream. This contributes to a simple 2S wake mode. The beam moves exactly towards the other side as the intensity of the low-pressure centre on one side increases and moves downstream.
The instantaneous pressure distribution in the near field of the cylinder-beam system with
$G=120$
when (a)
$U_{r}=10.0$
and (b)
$U_{r}=13.0$
.
For
$U_{r}=13.0$
, more pressure centres are found in the middle and the end of the beam, suggesting that more complex small-scale vortices are induced by rapid changes in the velocity and displacement of the beam. The beam comes into contact with the centre of low pressure (A, B, D, E) and undergoes rapid changes in displacement and velocity. Figure 11 illustrates the situation where the beam is in contact with the centre of low pressure in both the upper (A, B, C) and lower (D, E, F) rebound zones. This implies that the vibration rebound effect of a beam is related to its contact with the centre of low pressure. This also indicates a strongly coupled fluid-structure feedback mechanism: flow separation near the beam enhances localised suction forces, which in turn accelerate the beam segment locally, causing sharp reversals in velocity and displacement. These flow-induced ‘snap-through’ events generate high-frequency velocity fluctuations at the beam tip, consistent with the phase portraits observed in figure 6. The emergence of multiple rebound zones can be interpreted as the result of modal competition between the first bending mode and higher frequency components introduced by the adaptive actuation. This competition alters the timing and spatial distribution of vortex shedding, effectively splitting coherent structures and triggering the 2P transition. This interaction is central to the metamaterial’s ability to suppress subharmonic modes and reconfigure wake topology.
The relationship between the vibration response of the beam and the pressure distribution with
$G=120$
when
$U_{r}=18.0$
.
The emergence of high-frequency fluctuations is also observed, particularly in high-gain-factor conditions. Piezoelectric meta-cells create localised stress concentrations. Their self-sensing and adaptive feedback enable non-uniform stiffness and precise actuation timing along the beam. These non-uniformities generate localised high-wavenumber structural waves during flapping events. The elastic waves propagate at different phase velocities, causing frequency separation and temporal spreading of the response. This structural dispersion is likely amplified by the feedback mechanism, leading to the observed rapid local changes in displacement and the appearance of high-frequency components in the vibration response (see figure 11).
Flapping control effect of the metamaterial beam: (a) space-time contour showing transversal fluid force acting on the beam, (b) the frequency spectra of the fluid force, and (c) the energy transfer rate across the length of the beam with
$G=120$
when
$U_{r}=10.0$
. The first row represents the ‘control off’ situation and the second row represents the ‘control on’ situation.
Flapping control effect of the metamaterial beam: (a) space-time contour showing transversal fluid force acting on the beam, (b) the frequency spectra of the fluid force, and (c) the energy transfer rate across the length of the beam with
$G=120$
when
$U_{r}=13.0$
. The first row represents the ‘control off’ situation and the second row represents the ‘control on’ situation.
Figures 12 and 13 show the space-time distribution of the transversal fluid force, the frequency spectra of the fluid force and the energy transfer rate across the length of the uncontrolled and controlled beam. The energy transfer rate (per unit length) of the flexible beam immersed in unsteady flow is defined (Kumar, Navrose & Mittal Reference Kumar, Navrose and Mittal2016; Furquan & Mittal Reference Furquan and Mittal2021) as
$P(x,t)=v_{n}(x,t) [p(x,t)]$
, where
$v_{n}(x,t)$
and
$[ p(x,t)]$
are the normal component of the beam velocity and the difference between normal pressure on the two sides. The positive value
$P(x,t)$
relates to the gain in energy of the beam. For
$U_{r}=10.0$
, the fluid force acting on the uncontrolled beam is mainly composed of low-frequency subharmonic components, which are irregularly distributed in space and time. The metamaterial modulates the surrounding unsteady flow by suppressing the subharmonic components. Besides, modulation of the FSI system is accomplished by changing the position of the fluid energy transfer in addition to suppressing vibrations. As shown in figure 12(c), the energy transfer rate of the controlled beam exhibits a concentration of positive values at the free end of the beam in contrast to the irregular distribution of the uncontrolled beam. The concentration indicates that the introduction of the metamaterial localises the vortex shedding behaviour of the beam at the tip rather than in the centre of the beam, effectively preventing the generation of complex fluid excitation. This limits structural vibration to first-order mode-dominated LCO or quasi-periodic vibration with a small amplitude, rather than subharmonic vibration with a large amplitude. However, as for
$U_{r}=13.0$
, strong localised actuation at the mid-span of the beam leads to a new energy input site emerging there. This correlates with the pressure-driven rebound phenomena observed in figure 11, where low-pressure centres at the mid-span enhance local beam motion and serve as energy injection points. These dual input zones (tip and mid-span) correspond to a mixed-mode regime, the CL-2P mode, where the primary structure mode and high-frequency components coexist. The distributed energy transfer at multiple locations promotes vortex splitting.
3.2.3. Parametric sensitivity of actuation on flapping response
Figure 14 illustrates the impact of the gain factor on the flapping control effect. The increased gain factor enables the metamaterial to be more responsive. The flapping control of the cylinder-beam system is realised during the original subharmonic region. The amplitude of the controlled beam decreases as the gain factor increases. A stronger gain factor also results in a more delayed jump point for the beam, implying a stronger actuation capability to suppress the large-amplitude bifurcation behaviour. The increased gain factor also leads to a stronger vibration rebound effect. The rapid reversal in displacement, velocity and acceleration at the free end of the beam with maximum displacement is transferred to the vortices. The paired vortex moves up and down, respectively, to form a larger wake region and a 2P mode rather than a 2S mode, which is shown in figure 14(c). Although increasing the gain factor leads to reduced vibration amplitude and effective suppression of period-doubling bifurcations, figure 14(c) reveals a noticeable broadening of the wake. With a higher gain factor, the beam introduces stronger and more localised velocity gradients near its tip, particularly during the deceleration and reversal phases of its motion. These gradients amplify the shear layer interactions and promote the entrainment of surrounding lower-momentum fluid into the wake. As a result, the wake width increases despite the reduced structural displacement. This suggests that the control-induced modulation does not merely suppress flapping but actively alters the vortex formation process and increases the wake’s momentum mixing capacity.
The effect of the gain factor: (a) the vibration response, (b) the frequency response of the metamaterial beam with various gain factors, and (c) the wake modes of the cylinder-beam system when
$U_{r}=18.0$
with
$G=80$
(A),
$G=100$
(B) and
$G=120$
(C).
Figure 15 shows the effect of
$N_{p}$
on the flapping control effect of the metamaterial beam. In the case of
$U_{r}=13.0$
with
$G=120$
, the investigation of
$N_{p}$
is accomplished through the activation of varying numbers of piezoelectric patches (see figure 15
a). As
$N_{p}$
increases, the vibration amplitude of the beam decreases. The dominant frequency of the beam vibration shifts away from the subharmonic frequency. The enhancement of the patch number allows for enhanced control of the wake pattern (see figure 15
c). The presence of the characteristic vortex pairs with a subharmonic response off the main axis is observed when the patch number is relatively small (A), indicating that the period-doubling bifurcation is not fully suppressed. In contrast, at larger patch numbers (B and C), regular 2S and 2P modes are observed. The period-doubling bifurcation is suppressed. The wake mode is changed from an irregular vortex street to a stable vortex street. In the context of flapping control of a beam under the uniformly distributed excitation of the flow, the benefits of the self-sensing and adaptive actuation piezoelectric metamaterial control effect can be more effectively leveraged.
The effect of the number of the piezoelectric patches: (a) the vibration response, (b) the frequency response of the metamaterial beam against
$N_{p}$
, and (c) the wake modes of the cylinder-beam system when
$U_{r}=13.0$
with
$G=120$
and
$N_{p}=4$
(A),
$N_{p}=10$
(B) and
$N_{p}=14$
(C).
The parametric investigation across a range of gain factors and piezoelectric patch numbers reveals that the effectiveness of bifurcation suppression is dependent on the configuration of the piezoelectric metamaterial and the FSI. Specifically, there exists a moderate gain-factor range that allows the piezoelectric beam to block the onset of period-doubling bifurcations by eliminating frequency drift. When the gain factor is too low, the beam remains passive to wake instabilities; when the gain factor is too high, excessive local actuation introduces high-frequency perturbations, which may destabilise the vortex street and introduce quasi-periodicity or even chaos. The optimal gain window is thus system dependent. Furthermore, short beams may need lower gains for stability, while longer beams may benefit from higher gains to counter flow-induced loads.
3.3. Mechanisms of S-B and recovery transitions
3.3.1. Bifurcation dynamics and mode topology evolution
As illustrated in figure 16, the beam deflects away from the centreline and maintains a stable vibration mode with low amplitude at a low
${Re}$
and high
$U_{r}$
. Due to the unstable equilibrium point of the fluid moment acting on the beam immersed in the back-flow region and non-back-flow region, the beam undergoes a S-B bifurcation behaviour similar to that of an inverted pendulum in the parameter space when the Reynolds number and the length of the beam are appropriate. The S-B bifurcation is sensitive to the variation of the parameters. Although the changes in the system symmetry are intriguing in physics and promising in engineering applications, very few researchers focus on this phenomenon and provide feasible methods for symmetry modulation. Self-sensing and adaptive actuation piezoelectric metamaterial is able to adaptively sense the changes in local strains of the beam and generate corresponding deformation. The adaptive actuation mechanism with self-sensing capability can modulate the symmetry of the system adaptively and precisely for the parameter-sensitive system by adjusting the configuration of the proposed metamaterial and the gain factor.
Symmetry modulation of the piezoelectric metamaterial beam when
$U_{r}=15.0$
and
${Re}=200$
: (a) bifurcation diagrams, (b) the deflecting angle of the beam, (c) the amplitude of the beam, (d) the dominant frequency of the beam against different gain factors.
Figure 16(a) shows the bifurcation diagram of the controlled beam when
$U_{r}=15.0$
and
${Re}=200$
with different gain factors. Figures 16(b), 16(c) and 16(d) denote the deflecting angle, the vibration amplitude and the dominant frequency of the beam with the corresponding gain factors. Perturbation caused by the additional stiffness of the metamaterial leads to symmetry modulation for the beam. The S-B mode and symmetry-recovery (S-R) mode can be identified as two modes of the deflecting beam. As the gain factor increases, the beam maintains a S-B mode with a decreasing deflecting angle and an increasing vibration amplitude (stage I). Then, a quasi-periodic transition (stage II) leads to a symmetry reversal for turning the S-B mode into the S-R mode. The beam vibrates symmetrically about the central axis with a large amplitude (stage III). Stronger gain factors also bring more high-frequency components due to stronger actuation and lead to a quasi-periodic S-R mode (stage IV) eventually.
Such S-B and S-R mode transitions also manifest in the fluid response of the cylinder–beam system. Figure 17 illustrates the modulation of the lift and drag force coefficients on the cylinder as a function of the gain factor. The results also reveal a distinct symmetry transition bifurcation, in which two fundamentally different regimes (S-B mode and S-R mode) exert markedly different influences on the fluid force coefficients. In the low-gain regime (S-B mode,
$G\lesssim 10$
), relatively small lift force amplitudes and a persistent non-zero mean lift force on the cylinder are observed, as shown in figures 17(a) and 17(b). The directional aerodynamic bias may be exploited in micro flying robot applications for passive flow deflection or steering control. As the gain factor increases, the transition from small-amplitude, net-biased lift to large-amplitude, zero mean lift highlights a non-trivial trade-off between force magnitude and directional bias. This mechanism offers potential utility for energy harvesting systems. The drag response also exhibits a systematic trend with increasing the gain factor. As shown in figures 17(c) and 17(d), both the root-mean-square (r.m.s.) and mean values of the drag coefficient increase gradually with the gain factor. This growth suggests that the transition toward symmetric, large-amplitude vibration leads to enhanced wake unsteadiness and increased momentum entrainment from the mean flow. The intensified vibrations of the beam in the symmetric regime generate stronger periodic perturbations in the near wake, promoting vortex roll-up, thereby increasing unsteady lift and drag loading. The rising drag is thus not merely a by-product of gain factors, but a consequence of increased energy extraction from the mean flow to support the motion of the coherent vortex flow structure.
Fluid force coefficients modulation effect of the piezoelectric metamaterial beam when
$U_{r}=15.0$
against gain factors: (a) the r.m.s. value, (b) the mean value of the lift force coefficient, (c) the r.m.s. value and (d) the mean value of the drag force coefficient on the cylinder.
3.3.2. Stability landscape and transition bridging regime
While the previous section outlined the observable mode transitions, it remains critical to characterise the stability properties of each regime. Figure 18(a) shows the evolution of the largest Floquet multiplier from the S-B to the S-R mode. It is evident that the system remains linearly stable (
$\max (\| {\rm \mu} \| )\lt 1$
) in both the S-B and S-R regimes. However, within the intermediate range of
$G=5.5\!-\!10$
, the maximum multiplier slightly exceeds 1, indicating the loss of linear stability of the corresponding periodic orbit. The instability range aligns well with the parameter region where a sharp transition in force coefficients and mode shapes is observed, suggesting a dynamical route connecting two distinct stable periodic regimes through a quasi-periodic transitional state. Figure 18(b) provides a complementary view by mapping the Floquet exponents in the complex plane. In the low-gain regime (dark red points,
$G\lesssim 5$
), the exponents are clustered near the imaginary axis, with small real parts, indicating a stable periodic solution. As the gain factor increases towards the instability regime, a subset of Floquet exponents drifts rightward (toward positive real parts) near
${\textrm{Re}}(\lambda )=0$
. For larger
$G\gtrsim 10$
, the Floquet exponents shift back toward the imaginary axis and real parts remain negative. This marks the restoration of periodic stability, indicating that the high-gain S-R mode is a stable, large-amplitude, symmetric limit cycle. The evolution of Floquet exponents in this transition reveals not only the stability boundary but also the nature of the symmetry transitions: from a stable asymmetric periodic orbit, through a quasi-periodic torus regime, and eventually toward a stable symmetric periodic orbit. The quasi-periodic regime serves as a bridge between modal topologies and highlights the delicate sensitivity of the FSI system to the gain factor. These results provide a linearised characterisation of the nonlinear response, affirming that the observed mode switch is not a smooth parameter continuation but is instead punctuated by local loss of stability and dynamic complexity.
Stability analysis of the piezoelectric metamaterial beam flapping from S-B mode to S-R mode based on the Floquet theory: (a) the Floquet multipliers, and (b) the Floquet exponents distribution in complex space against the gain factors of the beam nonlinear vibration response.
3.3.3. Energy transfer restructuring and efficiency enhancement
The wake recovery rate
$R(\tilde{x})$
quantifies the degree to which the velocity in the wake recovers to the free-stream velocity
$U_{\infty }$
, reflecting the redistribution of kinetic energy within the wake. Here
$R(\tilde{x})$
is defined as
\begin{equation} R(\tilde{x})=\frac{U(\tilde{x})-U_{\textit{min}}}{U_{\infty }-U_{\textit{min}}}, \end{equation}
where
$U(\tilde{x})$
represents the streamwise flow velocity at a distance
$\tilde{x}$
downstream of the beam tip, while
$U_{\textit{min}}$
corresponds to the minimum flow velocity recorded along the post-tip sampling line (as shown in figure 24 in Appendix B). When fluid flows past the cylinder-beam system, a low-velocity region forms downstream, accompanied by vortex shedding and momentum loss. Evaluating the wake recovery rate of the beam allows for an assessment of the momentum redistribution effects of the controlled system compared with the uncontrolled system. Figure 19(a) shows the wake recovery rate of the system at
$\tilde{x}=5D , \tilde{x}=10D$
and
$\tilde{x}=30D$
under typical conditions (
$G=4,\, 8,\, 10,\, 20\text{ and }100$
, corresponding to stage I, II, III and IV, respectively). Although the system transitions from the S-B mode to the S-R mode, the relative variation in wake recovery rate remains relatively small. Within the studied parameter range,
$R(\tilde{x}=5D)$
consistently exceeds 0.6, and its growth gradually decelerates with increasing downstream distance
$\tilde{x}$
. A comparison of
$R(\tilde{x})$
across different stages reveals that a higher recovery rate is observed in stage III (
$G=10$
), suggesting that the flow velocity recovers to the free-stream condition. This indicates a faster wake recovery, where the symmetric vibration of the piezoelectric beam exerts a shorter-lived disturbance on the vortex shedding. In contrast, a lower wake recovery rate is observed in stage I (
$G=4$
), with the streamwise velocity failing to fully restore to the free-stream condition even at far-downstream positions. This implies that the asymmetric, small-amplitude vibration of the beam has a weaker influence on the flow field, leading to a prolonged, slow-recovery wake structure. In the S-B mode (small
$G$
), the piezoelectric beam maintains a small-amplitude LCO with symmetry breaking. The wake is dominated by cylinder-induced vortices. The vibration system exerts minimal influence on the flow field. Conversely, in the S-R mode (large
$G$
), the system undergoes symmetry recovery, with strong vibration energy input into the flow field. The resulting disturbances attenuate rapidly in the wake, accompanied by enhanced energy dissipation. The system exhibits superior momentum re-injection characteristics, wherein flow disturbances (particularly in the wake region) recover streamwise momentum and re-establish near-free-stream velocity profiles.
(a) Wake recovery rate of the system
$R(\tilde{x})$
at
$\tilde{x}=5D , \tilde{x}=10D$
and
$\tilde{x}=30D$
with different gain factors. (b) Energy harvesting efficiency
$\eta$
of the modulated system considering the response of the beam tip and the whole beam.
Energy harvesting efficiency
$\eta$
represents the capacity of a piezoelectric beam to extract kinetic energy from the fluid, which is defined as
\begin{equation} \eta =\frac{C_{\eta }\left\langle \dot{\xi }^{2}\right\rangle }{\frac{1}{2}\rho U_{\infty }^{3}DH_{b}}, \end{equation}
where
$\langle \dot{\xi }^{2}\rangle$
is the mean square value of the beam velocity and
$C_{\eta }$
is the efficiency coefficient (set
$C_{\eta }=1$
). Figure 19(b) demonstrates the energy harvesting efficiency
$\eta$
of the modulated system considering the response of the beam tip and the whole beam. For cantilever beams dominated by first-mode oscillations, the tip exhibits both maximal vibration response and peak energy harvesting efficiency. The global efficiency is necessarily lower due to the constrained energy extraction near the beam root, where flow separation and low-momentum regions persist. This spatial heterogeneity underscores the importance of modal dominance in optimising energy harvesting systems, as averaging over the beam length inherently attenuates the net efficiency. The contrast between the S-B and S-R modes is most pronounced in their energy harvesting efficiency. For low-gain-factor values, the system exhibits markedly poor efficiency (
$\eta \lt 0.04$
), indicative of a bottleneck in fluid-structure energy transfer. In this regime the energy exchange remains inefficient, with minimal momentum coupling between the beam and the flow. Conversely, at higher gain-factor values, the beam transitions to a symmetric large-amplitude LCO, characterised by enhanced velocity fluctuations and a pronounced improvement in energy harvesting efficiency (
$\eta \gg \eta _{{S}\hbox{-}{B}}$
). This shift signifies the establishment of an efficient energy transfer pathway, where the kinetic energy extraction by the piezoelectric beam is amplified. The underlying mechanism relies on a S-B bifurcation: by modulating the gain factor
$G$
, the system can be driven into a high-efficiency regime (
$\eta \gt 0.6$
), demonstrating precise controllability over energy harvesting performance.
While § 3.2.3 provides a parameter-based rationale for selecting an optimal gain factor, its ideal value is not universal but must be adjusted depending on the specific control objective. In flapping suppression applications the gain factor should be tuned to minimise flapping amplitude and suppress bifurcation occurrence. A conservative gain-factor value that enforces a symmetric, low-amplitude LCO is optimal. This is particularly useful for stealth, noise suppression or passive flow alignment applications. In contrast, for energy harvesting purposes, it is often beneficial to operate the system in a high-amplitude, symmetric regime, where the flapping energy extraction is maximised. This requires a higher gain factor to transition the beam into the S-R mode, where energy harvesting efficiency increases significantly. However, exceeding the threshold gain may lead to efficiency saturation or even degradation. Therefore, the gain-factor parameter acts as a control dial between two regimes: stability-focused control and energy-maximising control. The trade-off between stability and energy extraction must be considered carefully when designing application-specific control strategies.
Symmetry modulation effect of the piezoelectric metamaterial beam when
$U_{r}=15.0$
: the deformation envelope of the beam (first row), the transversal displacement at the beam tip (second row) and the vorticity contour around the beam (third row) with (a)
$G=4$
, (b)
$G=8$
and (c)
$G=10$
.
Figure 20 demonstrates the vibration response and the corresponding wake patterns of the beam in the S-B mode, the transition mode and the S-R mode when
$U_{r}=15.0$
with gain factor
$G=4 , G=8$
and
$G=10$
, respectively. With the increase of the gain factor, the envelope of the beam gradually approaches the centreline with an increasing amplitude. The metamaterial introduces more frequency components to suppress the S-B bifurcation, which promotes the transition from the deviated S-B mode to the S-R mode. Though 2S wake patterns are observed in both the S-B mode and S-R mode, the vibration of the S-R mode with large amplitude stretches the shedding vortex. The asymmetric vibration response (S-B mode) implies a net lift force acting on the cylinder, while the symmetric vibration response (S-R mode) indicates a lift force with zero mean value on the cylinder (see figure 17). The corresponding lift transition may enable the development of flexible propulsion systems for small-scale aircraft, underwater vehicles and micro flying robot platforms. By adjusting the gain factor, the system can be set to the S-B mode, generating a net lift force with an average value greater than zero, ensuring stable navigation. In contrast, when the system operates in the S-R mode, the average lift coefficient approaches zero, prompting the system to enter a stationary state. Besides, vibration amplitude variations between figures 20(a) and 20(c) demonstrate the feasibility of implementing symmetry modulation in underwater energy harvesting systems. Through the transition from S-B to S-R mode under specific flow conditions, the system achieves amplified vibration amplitudes, facilitating efficient ocean-to-mechanical energy conversion. This mechanical energy is subsequently transformed into electrical energy via the piezoelectric effect, providing sustainable power for unmanned underwater systems.
Dynamic mode decomposition of the wake of the deflecting beam when
$U_{r}=15.0$
with (a)
$G=4$
and (b)
$G=10$
. The corresponding DMD mode energy distribution for (c)
$G=4$
and (d)
$G=10$
.
To identify the difference between 2S wake patterns of the beam with the S-B mode and S-R mode, the frequency response of the beam and the corresponding DMD analysis are presented in figure 21. For the S-B mode with
$G=4$
, the actuation of the metamaterial perturbs the beam for slight changes of deflecting angle and vibration amplitude, but the beam still maintains the S-B mode. The vibration with the S-B mode is dominated by the fundamental frequency
$f_{1}$
and the second-order harmonic frequency
$2f_{1}$
. However, the influence of the additional stiffness on the potential energy of the beam makes the beam with the S-R mode (
$G=10$
) vibrate symmetrically at a lower frequency and a larger amplitude. The vibration with the S-R mode is dominated by the fundamental frequency
$f_{1}$
and the third-order harmonic frequency
$3f_{1}$
. As shown in figures 21(a) and 21(b), although the dominant DMD modes are comparable in the S-B mode and S-R mode, the dominant mode in the S-R mode has significantly low-frequency characteristics and wider spatial distribution, which is distinct evidence that the vortex is being stretched. Besides, the noticeable second-order frequency component in the vibration response of the S-B mode indicates that the vibration is affected by the periodic vortex shedding from the symmetrical circular cylinder at the front. The third-order component in the vibration response of the S-R mode means a strong interaction of the nonlinear vibration of the beam and the unsteady detached flow. As shown in figures 21(c) and 21(d), the energy contribution of the S-B dominant DMD mode is higher than that of the S-R mode, implying that the former is dominated by stable flow with periodic vortex shedding, while the latter has more harmonic components due to shedding vortex stretching caused by large-amplitude vibrations of the beam.
Bagheri et al. (Reference Bagheri, Mazzino and Bottaro2012) have mentioned that the broken symmetry of the cylinder-beam system is induced by the vorticity brought by the deformation and the motion of the beam through changing the pressure distribution near the cylinder wake. The length of the beam is proven to be an effective parameter that determines the vortex generation and pressure distribution near the cylinder wake. However, the proposed metamaterial can change the pressure distribution by only turning on and off the circuit rather than modifying the length of the beam. As shown in figure 22, for the S-B mode, the negative vortices at the free end of the beam are less intense. The shear layer around the beam is dominated by the positive vortices. In contrast to the S-B mode, with the increase of the gain factor (symmetry enhancement of the system), the free end of the beam in the S-R mode induces a stronger negative vorticity. The vorticity tends to merge with the negative vortex sheet in the lower region of the beam, promoting vortex shedding and stretching. The difference in the vorticity near the cylinder leads to a difference in pressure distribution as well.
Instantaneous vorticity (first row) and time-averaged pressure distribution (second row) of the metamaterial beam when
$U_{r}=15.0$
with (a,d)
$G=4$
, (b,e)
$G=8$
and (c,f)
$G=10$
.
Figure 22(d–f) illustrates the time-averaged pressure distribution near the beam and the cylinder when
$U_{r}=15.0$
with
$G=4 , G=8$
and
$G=10$
. Two symmetrical low-pressure centres are identified at the upper and lower positions of the circular cylinder for the S-R beam. In contrast, only one low-pressure centre is discovered at the lower position of the circular cylinder for the S-B beam. The difference in pressure distribution is also exhibited at the upper and lower sides of the beam. For the S-B beam, a much sparser pressure distribution is observed along the upper side of the beam while the lower side preserves a denser pressure distribution. In distinct contrast to the S-R mode, the pressure distribution is uniform on both sides of the beam. Meanwhile, the other significant difference in pressure distribution between the S-B mode and the S-R mode of the beam is located at the free end of the beam. As shown in figure 22(f), a distinct low-pressure centre is identified at the free end of the beam, which indicates that the vortex generation is related to the motion of the free end of the beam. The generation of vortices in the S-B mode, however, is more likely to be caused by flow separation from the stationary circular cylinder and the beam with low-amplitude vibration. As the gain factor is increased, the beam gradually induces and stretches the negative vortex at the free end, modifying the pressure distribution as well as the vortex shedding pattern. In contrast to the influence of the length of the beam on the surrounding vortices shedding, as Bagheri et al. (Reference Bagheri, Mazzino and Bottaro2012) have discussed, the transition brought by the proposed metamaterial from S-B to S-R is essentially from a flow-separation-dominated FSI mode to a vibration-dominated FSI mode.
The transition from S-B to S-R modes represents a fundamental reorganisation of energy exchange mechanisms within the FSI system. In the S-B regime the beam undergoes weak, deflected oscillations primarily driven by asymmetric vortex shedding induced by flow separation around the cylinder. The structural feedback to the flow is limited, leading to low energy harvesting efficiency, slow wake recovery and a persistent net lift force. As the gain increases, the system enters a S-R regime characterised by large-amplitude vibrations of the beam. These oscillations actively modulate the wake, enhance unsteady vortex shedding and inject energy into the surrounding flow, resulting in a shift towards a vibration-dominated wake structure. This transformation is corroborated by the emergence of high-frequency components in the beam response, increased Floquet stability and improved energy extraction from the flow.
This transition is further clarified by examining the spatial redistribution of energy transfer along the beam. In the S-B regime, energy input from the fluid is predominantly concentrated near the cylinder-beam junction, where pressure asymmetry persists. This localised interaction leads to weak excitation of structural modes and limited momentum feedback. In contrast, the S-R regime exhibits a tip-dominated energy input pattern, where the large deflections of the beam promote coherent vortex generation and sustain stronger fluid-structure coupling. The wake structure becomes more organised, the momentum deficit in the near wake recovers more rapidly and the energy harvesting efficiency increases substantially. Together, these findings indicate that the symmetry transition involves not only changes in vibration topology but also a reorganisation of how energy is transferred, distributed and dissipated within the coupled system.
4. Conclusions
This paper primarily provides a piezoelectric metamaterial design method for flapping control and symmetry modulation of a large deformable beam attached to a circular cylinder and subjected to a viscous flow. A strongly coupled fluid-structure-electro-control model of the piezoelectric metamaterial beam is established. The flapping dynamics of the system exhibit four typical flapping modes: the LCO mode, the subharmonic mode, the quasi-periodic mode and the S-B mode. Transitions among these modes are governed by nonlinear bifurcations, including periodic-doubling bifurcations and spontaneous S-B bifurcations. The piezoelectric meta-cells connected with unidirectional circuits can suppress periodic-doubling bifurcations of the beam through an adaptive actuation mechanism with self-sensing, leading to vibration amplitude reduction and delayed critical reduced velocity. The beam initially exhibiting the LCO mode can be transformed into a CL mode characterised by a vibration rebound effect, where rapid displacement-velocity inversion introduces multiple frequency components and localised perturbations. The perturbations alter the FSI mechanism, leading to significant changes in the downstream wake. Symmetry modulation is achieved through perturbations induced by the piezoelectric meta-cells, enabling transitions from S-B to S-R modes. The mode switching stems from the change of vortex generation patterns, i.e. the transition of the flow-separation-dominated vortex pattern to the vibration-dominated vortex pattern. This mode switching reorganises the energy exchange pathways between the fluid and the structure. Quantitative evaluation using wake recovery rate and energy harvesting efficiency demonstrates that the control system not only modifies flapping kinematics but also enhances the system’s ability to extract and regulate energy from the flow. The demonstrated symmetry modulation and vibration control capabilities offer significant potential for underwater energy harvesting systems, where mode switching optimises energy conversion efficiency across varying flow conditions.
Funding
This work was supported by the National Natural Science Foundation of China (Grant Nos. 12425202, U2141244, 12393781, 12121002, 12202260, 12402021), the Oceanic Interdisciplinary Program of Shanghai Jiao Tong University (Grant No. SL2021ZD104), the Project supported by Shanghai Municipal Science and Technology Major Project, and National Key R&D Program of China (Grant No. 2022YFB3902702, SQ2022YFB3900050). P.G. acknowledges the support from the National Natural Science Foundation of China (Grant No. 12202267), Shanghai Pujiang Program (Grant No. 22PJ1405300).
Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Appendix A. Modelling details of the piezoelectric metamaterial beam
The dynamic equation describing the deformation and motion of the piezoelectric beam is presented in § 2.3. This section will elaborate on the detailed derivation of the control equation for the structural model and its finite element discretisation. The deformations of the piezoelectric metamaterial beam in the
$x-z$
plane (see figure 2) can be described by the displacement field
$\xi =[ \tilde{u} \tilde{w} ]^{T}$
, given by
\begin{align} \tilde{u}(x,z,t)&=u(x,t)-z\frac{\partial w}{\partial x},\nonumber\\[2pt] \tilde{w}(x,z,t)&=w(x,t), \end{align}
where
$\tilde{u}$
and
$\tilde{w}$
are the generalised displacement components of the beam at the centreline in the
$x$
and
$z$
directions, respectively. When the frequency of the vortex shedding of the flow matches the natural frequency of the beam, nonlinear vibrations of the beam with large deflections may occur. Thus, the von Kármán strain-displacement assumption is employed here for the beam by considering the geometrical nonlinearity, given as
\begin{align} \varepsilon _{xx}&=\frac{\partial \tilde{u}}{\partial x}+\frac{1}{2}\left(\frac{\partial \tilde{w}}{\partial x}\right)^{2},\\[2pt] \gamma _{xz}&=\frac{\partial \tilde{u}}{\partial z}+\frac{\partial \tilde{w}}{\partial x}, \end{align}
where
$\varepsilon _{xx}$
and
$\gamma _{xz}$
are the in-plane strain and shear strain of the beam, respectively.
The piezoelectric metamaterial beam consists of the beam and the piezoelectric patches. The material properties of the beam follow Hooke’s law, and piezoelectric patches are characterised by linear piezoelectric constitutive equations:
\begin{equation} \begin{array}{ll} & \sigma _{ij}=C_{ijkl}\varepsilon _{kl}-e_{kij}E_{k},\\[4pt] & D_{i}=e_{ikl}\varepsilon _{kl}+g_{ik}E_{k}. \end{array} \end{equation}
Here
$\sigma _{ij}$
and
$\varepsilon _{kl}$
are the stress and strains of the piezoelectric patches,
$E_{k}$
and
$D_{i}$
denote the electric field vector and the electric displacement vector of the piezoelectric patches and
$C_{ijkl} , e_{kij}$
and
$g_{ik}$
are elastic, piezoelectric and dielectric material constants, respectively. As the thickness of the piezoelectric patches is much smaller than their length, a linear relationship between the electric field and the electric potential field is assumed:
\begin{equation} E_{k}=-\frac{\phi _{k}}{H_{p}}. \end{equation}
The dynamic equation of the piezoelectric metamaterial beam in fluid can be obtained by
\begin{equation} \int _{t_{1}}^{t_{2}}(\delta T-\delta U+\delta W){\textrm{d}}t=0, \end{equation}
where
$\delta T$
and
$\delta U$
are the virtual kinetic energy and potential energy of the piezoelectric metamaterial beam,
$\delta W$
is the virtual work of the beam contributed by the fluid forces. The kinetic energies of the piezoelectric metamaterial beam consist of the kinetic energy of the beam and the piezoelectric patches on the top and the bottom of the beam, given as
\begin{align} T &= \frac{1}{2}\int _{0}^{L_{b}}\rho _{s}\left[\left(\dot{u}-z\frac{\partial \dot{w}}{\partial x}\right)^{2} \!+\! \dot{w}^{2}\right]B_{s}{\textrm{d}}x{\textrm{d}}z+\frac{1}{2}\sum _{k=1}^{N_{p}^{t}}\int _{0}^{L_{p}}\\ \nonumber &\quad\times \rho _{p}\left[\!\left(\dot{u} \!-\! z\frac{\partial \dot{w}}{\partial x}\right)^{2}+\dot{w}^{2}\right]B_{p}{\textrm{d}}x{\textrm{d}}z + \frac{1}{2}\sum _{k=1}^{N_{p}^{b}}\int _{0}^{L_{p}}\rho _{p}\left[\!\left(\dot{u}-z\frac{\partial \dot{w}}{\partial x}\right)^{2}+\dot{w}^{2}\right]B_{p}{\textrm{d}}x{\textrm{d}}z, \end{align}
where
$N_{p}^{t}$
and
$N_{p}^{b}$
is the number of piezoelectric patches on the top and bottom of the beam, respectively,
$\rho _{s}$
and
$\rho _{p}$
are the density of the beam and the piezoelectric patches, respectively, and
$B_{s}$
and
$B_{p}$
are the width of the beam and the piezoelectric patches. The dot above the displacement variable represents differentiation to time.
The potential energy
$U$
is also contributed by both the beam and the piezoelectric patches, given as
\begin{align} U&=\frac{1}{2}\int _{0}^{L_{b}}\sigma _{ij}\varepsilon _{ij}B_{s}{\textrm{d}}x{\textrm{d}}z+\frac{1}{2}\sum _{k=1}^{N_{p}^{t}}\int _{0}^{L_{p}}\left(\sigma _{ij}\varepsilon _{ij}-D_{i}E_{i}\right)B_{p}{\textrm{d}}x{\textrm{d}}z \nonumber\\ &\quad +\,\frac{1}{2}\sum _{k=1}^{N_{p}^{b}}\int _{0}^{L_{p}}\left(\sigma _{ij}\varepsilon _{ij}-D_{i}E_{i}\right)B_{p}{\textrm{d}}x{\textrm{d}}z. \end{align}
The virtual work contributed by the fluid forces is defined as
\begin{equation} \delta W=\int _{0}^{L_{b}} \left[f_{u}(x)\delta u+f_{w}(x)\delta w \right]B_{s}{\textrm{d}}x, \end{equation}
where
$f_{u}$
and
$f_{w}$
is the shear stress and normal stress components of the fluid forces acting on the neutral axis of the metamaterial beam. As the shear stress of the fluid is much smaller than the normal stress,
$f_{u}$
is neglected in the computation. The dynamic equations of the piezoelectric beam influenced by the surrounding flow can be given by
\begin{align} \left[m_{\textit{eff}}^{b}+m_{\textit{eff}}^{p}\unicode{x1D6F7} (x)\right]\frac{\partial ^{2}w}{\partial t^{2}}+\left[D_{\textit{eff}}^{b}+D_{\textit{eff}}^{p}\unicode{x1D6F7} (x)\right]\frac{\partial ^{4}w}{\partial x^{4}}+e_{\textit{eff}}\frac{\phi }{H_{p}}\unicode{x1D6F7} (x) & = p_{\!f}(x,t),\nonumber\\ e_{\textit{eff}}\frac{\partial ^{2}w}{\partial x^{2}}+g_{p}\frac{A_{p}}{H_{p}}\phi _{\textit{sen}} & = 0, \end{align}
which are also presented in § 2.3. The effective density of the beam and the piezoelectric part are defined as
$m_{\textit{eff}}^{b}=\rho _{b}L_{b}$
and
$m_{\textit{eff}}^{p}=\rho _{p}(L_{p}^{t}+L_{p}^{d})$
, respectively. The effective bending stiffness of the beam and the piezoelectric part is given by
$D_{\textit{eff}}^{b}=E_{s}I_{s}$
and
$D_{\textit{eff}}^{p}=E_{p}(I_{p}^{t}+I_{p}^{d})$
. In the superscript and subscript of the symbol,
$(\boldsymbol{\cdot })^{b}$
(or
$(\boldsymbol{\cdot })_{b}$
) and
$(\boldsymbol{\cdot })^{p}$
(or
$(\boldsymbol{\cdot })_{p}$
) denote the beam element and the piezoelectric element, respectively, whereas
$(\boldsymbol{\cdot })^{t}$
and
$(\boldsymbol{\cdot })^{d}$
indicate the upper and lower piezoelectric patches, respectively. The effective piezoelectric material constant is presented by
$e_{\textit{eff}}=e_{p}(I_{p}^{t}+I_{p}^{d})$
. The piezoelectric distribution function is defined as
$\unicode{x1D6F7} (x)=\sum _{i=1}^{N_{p}}\varphi _{i}(x)$
, where
$\varphi _{i}(x)$
is defined as
$\varphi _{i}(x)=1$
, when
$x\in [x_{i},x_{i}+L_{p}]$
. Otherwise
$\varphi _{i}(x)=0$
.
The FEM is employed to establish the numerical model of the piezoelectric metamaterial beam. The beam is discretised using two-node elements. The axial displacement
$u$
and the transversal displacement
$w$
of the beam in each element are interpolated by linear and Hermite functions, respectively. The discretised equations of motion of the piezoelectric metamaterial beam can be written as
\begin{align} \boldsymbol{M}\ddot{\boldsymbol{\xi }}+\boldsymbol{C}\dot{\boldsymbol{\xi }}+(\boldsymbol{K}_{L}+\boldsymbol{K}_{\textit{NL}})\boldsymbol{\xi }-\boldsymbol{\varTheta }_{G} \boldsymbol{\phi }&=\boldsymbol{F}_{\!\!f},\nonumber\\ {\boldsymbol{\varTheta }_{S}}^{\textrm{T}}\boldsymbol{\xi }+\boldsymbol{C}_{p} \boldsymbol{\phi }&=\boldsymbol{0}, \end{align}
where
$\boldsymbol{\varTheta }_{G}$
and
$\boldsymbol{\varTheta }_{S}$
are reconstructed coupling matrices for the metamaterial beam with amplification circuit,
$\boldsymbol{C}_{p}$
is the capacitance matrix. Equation (A.10) can be further rewritten as
\begin{equation} \boldsymbol{M}\ddot{\boldsymbol{\xi }}+\boldsymbol{C}\dot{\boldsymbol{\xi }}+(\boldsymbol{K}_{L}+\boldsymbol{K}_{\textit{NL}}+\boldsymbol{K}_{G}) \boldsymbol{\xi }=\boldsymbol{F}_{\!\!f}, \end{equation}
where the effective stiffness matrix
$\boldsymbol{K}_{G}$
is defined as
$\boldsymbol{K}_{G}=\boldsymbol{\varTheta }_{G}\boldsymbol{C}_{p}^{-1}\boldsymbol{\varTheta }_{S}^{\textrm{T}}$
.
The natural frequency of a cantilever beam covered by a single piezoelectric patch with a short circuit and an open circuit.
Appendix B. Numerical validation and mesh convergence
The computational models of the piezoelectric metamaterial beam in fluid are validated. For validating the structural model of the beam, table 2 provides the natural frequencies of a cantilever beam covered by a single piezoelectric patch with a short and open circuit. The beam is assumed to be aluminum with elastic modulus
$E=70\, {\textrm{GPa}}$
, Poisson ratio
$\nu _{s}=0.33$
and density
$\rho _{s}=2700\, {\textrm{kg}}\,{\textrm{m}}^{-3}$
. The piezoelectric patches on the beam are assumed to be PZT-5A with the material properties as
$C_{11}=61\, {\textrm{GPa}} , C_{55}=21.1\, {\textrm{GPa}}, g_{33}=13.33\times 10^{-9}\, {\textrm{F}}\,{\textrm{m}}^{-1} , e_{31}=0.33\, {\textrm{C}}\,{\textrm{m}}^{-2}$
and
$\rho _{p}=7750\, {\textrm{kg}}\,{\textrm{m}}^{-3}$
. The results of the proposed method have a good agreement with the results of Erturk (Reference Erturk2011), which verifies the accuracy of the proposed electro-mechanical FEM model.
Table 3 provides the natural frequencies of a clamped-clamped beam covered by multiple piezoelectric patches (
$N_{p}=16$
) with a short (
$G=-100$
) and open circuit with an aluminum beam and PZA-5A piezoelectric patches. The length of the piezoelectric patch
$L_{p}=15\, {\textrm{mm}}$
, the height of the piezoelectric patch
$H_{p}=0.2\, {\textrm{mm}}$
and the length of the base beam
$L_{b}^{unit}=2L_{p}$
. The results obtained by the present method agree well with the converged solutions obtained by commercial software COMSOL. The piezoelectric metamaterial beam equipped with unidirectional amplification circuits possesses a non-Hermitian skin effect. When the excitation point is on the actuator side and the response point is on the sensor side of the beam, the elastic wave is suppressed during propagation, whereas it is enhanced when the excitation and response positions are reversed. Figure 23 shows the displacement responses of the piezoelectric metamaterial beam excited by external forces acting at different positions. Two points are selected at each end of the beam at a distance
$L_{p}/2$
from the boundary. These two points are used alternately as the excitation and response points. The external forces, set as
$F(t)=F_{0}\sin (2\pi f_{e}t)$
, are employed at the excitation point to obtain the frequency response function of the beam, where
$f_{e}$
is the excitation frequency of the external forces. In figure 23, ‘off’ and ‘on’ indicate that the beam is with a closed and open circuit, respectively. It is observed from figure 3 that the computed vibrational responses of the piezoelectric metamaterial beam are in good agreement with the solutions determined by COMSOL.
The natural frequency of a clamped-clamped beam covered by multiple piezoelectric patches (
$N_{p}=16$
) with a short circuit (
$G=-100$
) and an open circuit.
Validation of the electro-mechanical model of the piezoelectric metamaterial beam: (a) frequency responses of displacement of the beam (external force acting at the right side of the beam), (b) frequency responses of displacement of the beam (external force acting at the left side of the beam), (c) time-history responses of displacement of the beam (external force acting at the right side of the beam), (d) time-history responses of displacement of the beam (external force acting at the left side of the beam).
To verify the accuracy of the proposed FSI method, the present results are compared with those of Turek et al. (Reference Turek, Hron, Razzaq, Wobker and Schäfer2010) for a benchmark problem of an elastic beam attached to a fixed circular cylinder and subjected to uniform flow. The beam length is
$L=0.35\, {\textrm{m}}$
, the beam thickness is
$H=0.02\, {\textrm{m}}$
and the radius of the cylinder is
$r=0.05\, {\textrm{m}}$
. The material properties of the beam follow those of Turek et al. (Reference Turek, Hron, Razzaq, Wobker and Schäfer2010) with Young’s modulus
$E=5.6\, {\textrm{MPa}}$
, shear modulus
$G_{s}=2\, {\textrm{MPa}}$
, Poisson ratio
$\nu _{s}=0.4$
and density
$\rho _{s}=1000\, {\textrm{kg}}\,{\textrm{m}}^{-3}$
. A parabolic velocity distribution is assumed at the inlet boundary, which is defined as
\begin{equation} \begin{array}{l} u_{\!f}\!\left(t,y\right)=1.5\overline{U}_{\infty} \displaystyle\frac{y\left(y-y_{0}\right)}{\left(y_{0}/2\right)^{2}}\displaystyle\frac{1-\textrm{cos}(\pi t/2)}{2},\quad t\lt 2.0,\\[12pt] u_{\!f}\!\left(t,y\right)=1.5\overline{U}_{\infty }\displaystyle\frac{y\left(y-y_{0}\right)}{\left(y_{0}/2\right)^{2}},\qquad\quad \qquad \qquad\, t\geq 2.0, \end{array} \end{equation}
where
$\overline{U}_{\infty }$
is the maximum magnitude of the velocity at the inlet boundary and
$\overline{U}_{\infty }$
is set as
$\overline{U}_{\infty }=2\, {\textrm{m}}\,{\textrm{s}}^{-1}$
in this case. Table 4 illustrates the computation results of the proposed method and three other FSI coupling algorithms from the work of Turek et al. (Reference Turek, Hron, Razzaq, Wobker and Schäfer2010). Method 1 is the Rannacher method, which uses a monolithic variational formulation FSI solver based on a unified Eulerian framework to predict the flow of the fluid and the deformation of the structure. Method 2 is the Turek and Hron method, which uses a fully implicit, monolithic ALE-FEM approach for the FSI response of the beam. Method 3 is the Bletzinger method, which uses a partitioned approach that combines their in-house codes and OpenFOAM for the fluid and structural response. The velocity and corresponding frequency response in the
$x$
axis and
$y$
axis of the proposed method are consistent with the reference results.
Flow-induced vibration response of the flexible beam.
The beam attached to a fixed circular cylinder is investigated in this study. As shown in figure 24, the centre of the circular cylinder is fixed at the position with
$x=0,\, y=0$
. The inlet and outlet boundaries are placed at a distance of
$25D$
and
$50D$
from the centre of the cylinder, respectively. The top and bottom boundaries with no-slip boundary conditions are separated by a distance of
$50D$
. Uniform flow described by
$U_{\infty }$
is defined at the inlet boundary. Meanwhile, the no-slip condition is satisfied at the surface of the cylinder and the fluid-structure interface is defined at the surface of the beam. In this study, the length of the beam is
$L=3.5D$
, where the diameter of the cylinder is given by
$D=0.02\, {\textrm{m}}$
. The thickness and the width of the beam are defined as
$H_{b}=0.1D$
and
$W_{b}=1.0\, {\textrm{m}}$
. The density of the beam is
$\rho _{s}=2700{\textrm{kg}}\,{\textrm{m}}^{-3}$
, the density of the piezoelectric material is
$\rho _{p}=7780 {\textrm{kg}}\,{\textrm{m}}^{-3}$
and the elastic modulus and the piezoelectric coupling constant of the piezoelectric material are given by
$C_{11}=60\, {\textrm{GPa}}$
and
$e_{31}=-10.2\, {\textrm{C}}\,{\textrm{m}}^{-2}$
. The numerical investigation is carried out by considering Reynolds number
${Re}=200$
and
${Re}=400$
, mass ratio
$M_{r}=10$
and various values of the reduced velocity
$U_{r}$
.
Schematic of the computation domain for a cantilevered metamaterial beam attached to a circular cylinder in uniform viscous flow.
The flow domain, for a cantilevered beam attached to a circular cylinder, is discretised by a body-fitted unstructured two-dimensional mesh. The mesh, referred to as mesh 1, consists of 103 468 cells. Besides, a finer mesh, referred to as mesh 2, consisting of twice the number of cells in mesh 1, is provided for the mesh convergence validation. Table 5 lists the r.m.s. value of the normalised amplitude of vibration of the tip of the beam (defined as
$A_{{rms}}^{*}=A_{{rms}}/H_{b}$
), the r.m.s. value of the lift coefficient of the circular cylinder
$C_{{L},{rms}}$
and the mean value of the drag coefficient
$\overline{C}_{\textrm{D}}$
of the circular cylinder for flow past a fixed circular cylinder with a flexible metamaterial beam for
${Re}=400 , M_{r}=10.0 , U_{r}=5.0$
and
$N_{p}=16$
. The numerical results with a time-step size of
${\unicode[Arial]{x0394}} t=2\times 10^{-5}\, {\textrm{s}}$
from mesh 1 and mesh 2 are in excellent agreement. Only an error of 0.2 % exists between the computation results from the two meshes for the r.m.s. value of the normalised amplitude of vibration at the tip of the beam. Thus, mesh 1 is sufficient for all computations in this study. Besides, as shown in rows 1 and 2 of table 5, the results for
${\unicode[Arial]{x0394}} t=2\times 10^{-5}\, {\textrm{s}}$
and
${\unicode[Arial]{x0394}} t=1\times 10^{-5}\, {\textrm{s}}$
with mesh 1 are in excellent agreement as well. The largest deviation between the two numerical results is only 0.3 %. Thus,
${\unicode[Arial]{x0394}} t=2\times 10^{-5}\, {\textrm{s}}$
is used for all computations in this study.
Flow past a circular cylinder with a metamaterial beam for
${Re}=400 , M_{r}=10.0 , U_{r}=5.0$
and
$N_{p}=16$
: the effect of mesh resolution and time step size.
To mitigate the influence of the boundary effect in the flow field computation domain on the vibration of the FSI response, particularly for large-amplitude cases, in this study we performed a truncation convergence analysis of the computational domain associated with the flow field. The findings are summarised in table 6. This analysis was conducted under the condition where the system exhibits subharmonic vibrations with large amplitudes, considering the parameters
$M_{r}=10.0 , U_{r}=10.0 , N_{p}=16$
and
${Re}=400$
. Based on the baseline computational domain (as depicted in figure 24), the distances between the cylinder-beam region and the flow field boundaries were systematically varied. These variations included the lateral distance
$D_{{in}}$
from the inlet to the fixed cylinder, the lateral distance
$D_{{out}}$
from the outlet to the free end of the beam, and the longitudinal distance
$D_{{t}/{b}}$
from the top and bottom boundaries to the centre of the cylinder. Consequently, the vibration responses (r.m.s. amplitude
$A_{{rms}}^{*}$
, dominated harmonic frequency
$f_{\xi ,d}$
, subharmonic frequency
$f_{\xi ,s}$
) and flow field responses (r.m.s. lift coefficient amplitude
$C_{{L},{rms}}$
, average drag coefficient
$\overline{C}_{D}$
) of the FSI system with varying computational domain sizes were investigated. As shown in table 6, altering
$D_{{in}} , D_{{out}}$
and
$D_{{t}/{b}}$
yielded results that closely align with the baseline, demonstrating that the selected baseline computational domain ensures both convergence and consistency in the computations. Specifically, for systems undergoing large-amplitude subharmonic vibrations, adopting the computational domain illustrated in figure 20, along with the mesh configuration described earlier, guarantees the reliability of the numerical results presented in this study.
Flow past a circular cylinder with a metamaterial beam for
${Re}=400 , M_{r}=10.0 , U_{r}=10.0$
and
$N_{p}=16$
: the effect of computational domain size on key results.
Appendix C. Algorithm details of DMD
The DMD method can extract the dominant flow structures together with their corresponding frequency components. Based on these, the unsteady vorticity field
$\omega (x,y,t)$
can be reconstructed by the following equation,
\begin{equation} \omega (x,y,t)=\sum _{i=1}^{N}a_{i}\boldsymbol{\cdot }e^{{\lambda _{i}}t}\boldsymbol{\cdot }b_{i}(x,y), \end{equation}
where
$N$
is the number of DMD modes,
$a_{i}$
and
$\lambda _{i}$
are the corresponding initial value and eigenvalue of ith DMD mode, respectively, and
$b_{i}(x,y)$
is the ith DMD mode function. As a data-driven mode decomposition method for capturing coherent structures with specific frequency components, the DMD method can be realised through a snapshot-based approach. In DMD analysis, two contiguous vorticity fields can be described by a linear mapping,
\begin{equation} \boldsymbol{\omega }_{i+1}=\boldsymbol{A}\boldsymbol{\omega }_{i}, \end{equation}
where
$\boldsymbol{A}$
is the mapping matrix. The eigenfunctions and the corresponding eigenvalues of the mapping matrix represent the main fluid structures and frequencies. To obtain the mapping matrix, snapshots of the vorticity field
$\boldsymbol{\omega }_{i}$
with
$(n+1)$
time steps with
${\unicode[Arial]{x0394}} t$
between consecutive snapshots are collected for matrices
$\boldsymbol{\varOmega }_{x}=[\boldsymbol{\omega }_{0},\boldsymbol{\omega }_{1},\ldots ,\boldsymbol{\omega }_{n-1}]$
and
$\boldsymbol{\varOmega }_{y}=[\boldsymbol{\omega }_{1},\boldsymbol{\omega }_{2},\ldots ,\boldsymbol{\omega }_{n}]$
. Next, SVD is performed for matrix
$\boldsymbol{\varOmega }_{x}$
,
\begin{equation} \boldsymbol{\varOmega }_{x}=\boldsymbol{U} \! \boldsymbol{\varSigma }\boldsymbol{V}^{T}, \end{equation}
where
$\boldsymbol{U}$
and
$\boldsymbol{V}$
are unitary matrices and
$\boldsymbol{\varSigma }$
is a matrix containing the singular values of the data matrix
$\boldsymbol{\varOmega }_{x}$
. Next, a matrix truncation should be performed for matrix
$\boldsymbol{U} , \boldsymbol{V}$
and
$\boldsymbol{\varSigma }$
for capturing main modes of the system: Here
$\boldsymbol{U}_{m} , \boldsymbol{V}_{m}$
and
$\boldsymbol{\varSigma }_{m}$
.
$\boldsymbol{U}_{m}$
and
$\boldsymbol{V}_{m}$
are the results of the preservation of the first m columns of the matrix
$\boldsymbol{U}$
and
$\boldsymbol{V}$
while
$\boldsymbol{\varSigma }_{m}$
is obtained by preserving the first m diagonal elements of the matrix
$\boldsymbol{\varSigma }$
. The linear mapping matrix
$\tilde{\boldsymbol{\kern-2pt A}}$
can be given as
\begin{equation} \kern2pt \tilde{\boldsymbol{\kern-2pt A}}=\boldsymbol{\varOmega }_{y}{\boldsymbol{\varOmega }_{x,m}}^{\dagger }, \end{equation}
where
$\kern2pt \tilde{\boldsymbol{\kern-2pt A}}$
is used to represent the linear mapping matrix. The eigenvalues
$\lambda _{i}$
and eigenfunctions
$b_{i}$
can be obtained by solving the eigenvalue problem for the matrix
$\kern2pt \tilde{\boldsymbol{\kern-2pt A}}$
. Figure 25 demonstrates the details of the DMD algorithm. Here DC indicates data collection from the flow fields to the data matrix, MS indicates matrix assembling from the data matrix to matrices
$\boldsymbol{\varOmega }_{x}$
and
$\boldsymbol{\varOmega }_{y}$
, MT indicates matrix truncation for
$\boldsymbol{U} , \boldsymbol{V}$
and
$\boldsymbol{\varSigma }$
and MM indicates the matrix multiplication operation for
$\kern2pt \tilde{\boldsymbol{\kern-2pt A}}$
.
Schematic of the DMD algorithm details.
The eigenfunctions are spatially expanded to give the fluid structure (vorticity field) at the corresponding frequency, as shown in figures 9 and 21. The total energy
$E_{i}$
of any mode over the entire time period is defined by squaring each element of the modal function and then summing them together.
