Flutter Phenomenon in Flow Driven Energy Harvester–A Unified Theoretical Model for “Stiff” and “Flexible” Materials

Here, we report a stable and predictable aero-elastic motion in the flow-driven energy harvester, which is different from flapping and vortex-induced-vibration (VIV). A unified theoretical frame work that describes the flutter phenomenon observed in both “stiff” and “flexible” materials for flow driven energy harvester was presented in this work. We prove flutter in both types of materials is the results of the coupled effects of torsional and bending modes. Compared to “stiff” materials, which has a flow velocity-independent flutter frequency, flexible material presents a flutter frequency that almost linearly scales with the flow velocity. Specific to “flexible” materials, pre-stress modulates the frequency range in which flutter occurs. It is experimentally observed that a double-clamped “flexible” piezoelectric P(VDF-TrFE) thin belt, when driven into the flutter state, yields a 1,000 times increase in the output voltage compared to that of the non-fluttered state. At a fixed flow velocity, increase in pre-stress level of the P(VDF-TrFE) thin belt up-shifts the flutter frequency. In addition, this work allows the rational design of flexible piezoelectric devices, including flow-driven energy harvester, triboelectric energy harvester, and self-powered wireless flow speed sensor.

. Major process integration steps. (a) The piezoelectric 0.02 µm AlN /0.2 µm Mo/1.2 µm AlN stack were consecutively deposited on a SOI wafer; (b) Pattern the top and bottom Al electrodes; (c) Pattern the piezoelectric stack; (d) 20 µm Si+ 1 µm SiO 2 + 80 µm Si etch + Backside wafer Grinding to 550 µm; (e) Backside 250 µm partial Si etch; (f) Top cap wafer fabrication; (g) bonding top cap wafer with the device wafer and (h) Inlet and outlet hole forming by DRIE on top cap wafer, after that bonding bottom wafer to device wafer and final partial dicing to expose top electrodes. Table S1. The piezoelectric and mechanical properties of the "stiff" material AlN, "flexible" material P(VDF-TrFE) and silicon.

S2.1.1 Classical beam theory
First, the vibration of belt under zero air velocity is studied. Assuming that the deformation of bending and torsion for the belt is small enough, and when it bends and twists naturally, the motion of it follows Euler-Bernoulli theory for bending and St. Venant theory for torsion S1,S2 : where x is along the flow direction, z is along the axis of beam or belt, y is the normal direction of the belt surface, and t is time. v is the deformation of belt (or beam) in the y direction, and ϕ is the twist angle of each element of belt. I x is the second axial moment of area for belt about the x-direction. E is Young's modules, G is the shear modules, and m is the mass of the belt/beam per unit length, J is the torsional constant (stiffness), K m is the radius of gyration, A is the cross section area. For the thin rectangular belts investigated in the current study, J = Bh 3 3 S3 .

S2.1.2 Beam with axial force
With an axial force inside the belt, the governing equation of bending vibration can be rewritten as S4,S5 : where T is the axial internal force, which is a constant for belts under extremely pre-stretching. Similarly, the torsional vibration of belt with a constant internal axial force can be analyzed in the following way. Assuming the torsional angle ϕ at location z, then the angle θ (see Figure S4) between the deformed string in the belt and the original axis z is Figure S4. The schematic of the torsional deformation of the belt. Left hand side is the whole picture of belt which is divided into threads; right hand side is the section between the two blue lines and the internal stress in the thread.
Assuming the axial tension is uniform within the cross section of the belt, so it can be expressed in the form where A is the cross section area, and T is the total axial stretching force. Figure S5. The schematic of the torsional deformation of a string inside the belt.

5/16
Similar to the string, the rotational force per unit length of belt (see Figure S5) due to the different deformation angle of axial string on the small cross section area dA is The torque of this force about the axis OO (or z axis) is Integrate it over the cross section where I is the polar moment of area for rectangular beam defined as I = r 2 dA = B 3 h+Bh 3 12 . Then combing with torque due to the shear stress GJ ∂ 2 ϕ ∂ z 2 ( J = Bh 3 3 for rectangular thin belt in the current study S3 .), and considering the angular acceleration, the governing equation can be obtained as where m is the mass of belt per unit length, and K m is the radius of gyration, and mK 2 m = I α is the moment of inertial.

S2.1.3 large deformation
However, when the bending deformation/deflection of the belt is large, the axial force will be affected by the associated axial force then. Assuming this axial force is still a constant along the axis, so the additional strain of the belt along the axis due to the bending deformation is Thus, the total axial force considering bending deformation of the belt becomes where T 0 is the pre-existing axial force due to pre-stretching. Replacing the axial force in Eqs. (S2.3) and (S2.4) by the total axial force in Eq. (S2.5) yields When T 0 L 2 EI x , the bending stiffness is dominant, and then Eq. (S2.6) reduces to the theory of beam under large deformation given by Conway S6,S7 .
On the other hand, when T 0 L 2 EI x , the pre-stretching force T 0 is dominant, and then Eq. (S2.6) reduces to the Kirchhoff equation for nonlinear vibration of strings S8-S12 .

S2.2.1 Small deformation
When the deformations of bending is very small EA 3) and (S2.4) can be used to solve the problem. They are general linear partial differential equations (PDEs) with constant parameters, applying the variable separation approach as following: Let h (t) and α (t) respectively to be the bending deformation and torsional angle of belt/beam at a typical position along the axis (normally the maximum deformation along the axis), representing the features of vibration in time space. On the other hand, V (z) and Φ (z) are the deformation distribution along the axis, representing vibration mode (as shown in Figs. S3-??). Then Thus, the ordinary differential equations (ODEs) can be obtained instead of the original partial differential equations (PDEs): where ω h0 and ω α0 are the natural frequency of vibration for each mode of bending and torsion with and for torsion Thus, the vibration modes and respective frequency can be solved based on the boundary conditions of belts at its two ends. In the current study, this process is conducted in ABAQUS software.

S2.2.2 Large deformation
When the deformation of bending for the belt is large ( EA dz is close or even higher than T 0 ), the nonlinear vibration equations Eqs. (S2.6) and (S2.7) need to be used to solve the problem. Although these two PDEs are nonlinear, the leading order of the vibration can be linearized and approximately described bÿ where h (t) and α (t) are also respectively the maximum deformation along the axis similar to the situation with small deformation, representing the features of vibration in time space. Speciafically, Eqs. (S2.6) and (S2.7) were solved numerically first, and then the key features of the leading order of the vibrations were extracted using Fast Fourier transform (FFT) in order to obtain ω h0 and ω α0 .

S2.3 Analysis of 2DOF fluttering system
Thus, the ordinary differential equations (ODEs) for the deformations of bending and torsion at the typical position along the axis of belt (normally the maximum deformation along the belt's axis) are not only accurate for the vibrating belt with small deformation, but is also applicable for the leading order of belt's vibration with large deformation. Generally, there is friction and dissipation in the vibration system, bringing the damping effect into the governing equations. In addition, the external lift and torque due to the aerodynamic force are also added in, and then the governing equations are obtained for this two degree of freedom (2DOF) system S13-S17 S2.16) and where ζ h0 and ζ α0 are the original damping ratio for bending and torsion without incoming flow, and I α = mK 2 m . L and M are the aerodynamic lift force and torque per unit length of belt respectively.
According to Scanlan and Theodorsen's theory S13-S16 , the lift force and torque is determined by the local motion of belt/beam as where K = Bβ /U is the non-dimensional frequency, β is the final fluttering frequency of belt, B is the belt width in the flow direction, U is the velocity of the coming air flow, and ρ a is the density of air. Besides, H * i and A * i are the generalized flutter derivatives related to the cross section feature of the belt and the non-dimensional frequency K.
With substituting Eqs. (S2.18) and (S2.19) back into Eqs. (S2.16) and (S2.17) and rearranging the same order terms, it is obtained in matrix form that M eẌ + C eẊ + K e X = 0, (S2.20) where X = h α . It can be further simplified into the following form that Assuming its solution is y = Ψ e λt , then we got (λ R + S)Ψ = 0 whose Eigen vector is and the Eigen-function is where j is complex number √ −1, the subscript h is related to bending, and α is related to torsion. Solving such nonlinear eigenvalue problem needs an iterative process between Eqs. S2.18, S2.19 and S2.23 until β in them coincides. Normally, the mechanical dynamic matrix R −1 S of a 2DOF system has four Eigen-values, which are two pairs of conjugates and show the features of the whole vibrating system. ω i is the angular frequency of the new system, and χ i represents the damping feature of the coupled vibrating system. When χ i > 0, the system is stable, the corresponding amplitude will damp to zero after long period; however, when χ i < 0, it is an unstable system, the corresponding amplitude will magnify by itself until it breaks down. Since the flutter derivatives are also dependent on the resultant vibration frequency, so solving the fluttering equations needs iterations between time-domain and frequency domain in order to obtain the correct flutter frequency and the corresponding damping ratio S17 .
Then the typical bending deformation and the torsion angle in the time-domain can be written as h = A hh e −χ h t cos (β h t + θ hh ) + A hα e −χ α t cos (β α t + θ hα ) (S2.24) and α = A αh e −χ h t cos (β h t + θ αh ) + A αα e −χ α t cos (β α t + θ αα ) (S2.25) where χ h and β h corresponds to the damping factor and angular frequency due to bending, χ α and β α correspond to the damping factor and angular frequency due to torsion. When χ i > 0, the system with the corresponding vibration is stable at the specific frequency; however, when χ i < 0, it is an unstable system at the corresponding frequency. It can also be found that the both bending and torsional deformation can vibrate at two frequencies, of which one will finally dominate another (depending on the damping factor of each component).

The ratio between the bending and torsion can be determined by terms of Eigen-vector
where subscript 1 and 2 respectively means the first and second component of Eigen-vector, and i = h, α respectively means the Eigen-vector correlated to bending or torsion.

9/16
In addition, the final flutter frequency f i and damping ratio ζ i at different air incoming velocity can be derived as:

S2.4 Results for "stiff" material based on theoretical prediction
The above theory and the natural frequency obtained by ABAQUS (bending: 10128Hz, torsional: 13410Hz) are employed to calculate the vibration frequencies of belts based on "stiff" material and corresponding damping ratio as shown in Tab. S2 and Figure S6 in order to evaluate the onset of fluttering and flutter frequency. The vibration frequencies and damping ratios are grouped into two series: the torsional frequency, and the other is bending frequency. It shall be noted that bending and torsion motions contain two components which are associated to bending and torsional frequencies respectively (see Eqs. S2.24 and S2.25). It is observed that when the air flow velocity exceeds the critical value (75m/s herein), the torsion damping ratio ζ α becomes negative, so the motion components associated to torsional frequency becomes divergent for both bending and torsional motions. In other words, the bending and torsional amplitude become larger and larger when the air velocity exceeds onset velocity, and so the flutter occurs.

S3.1 "Stiff" AlN/Si belt
The Euler-Bernoulli model s shown in Figure S7. The AlN and Si substrate both bend about a common neutral axis which is no longer the neutral axis of the belt. Perfect bonding is assumed, and the AlN is considered to be a layer of the belt. This neutral axis is calculated by a modulus-weighted algorithm.

11/16
The equation for the distance to the neutral axis can be written as where E a , E b are the Young's Modulus of the AlN and Si belt; h a , h b are the thickness of the AlN and Si belt.
To simplify the calculations, the average strain in the AlN is determined and used to find the voltage. The average strain in AlN ε a is S3.4) and M is the bending moment of the belt. So .

(S3.5)
The voltage on the AlN surfaces is related to the stress by where g 31 is the AlN voltage constant and σ a is stress in AlN. Finally, the voltage output voltage can be obtained as

S3.2 "Flexible" P(VDF-TrFE) belt
The total charge generated from the deformation induced by the vibration of the P(VDF-TrFE) (see Figure S9) can be obtained by the equation: S3.9) where C = ε 0 ε r A p h p is the capacitance of the P(VDF-TrFE) belt, in Farads; A p = B p × l p is the surface area of the P(VDF-TrFE) belt on the top side; ε r is the relative static permittivity of the P(VDF-TrFE) (it is 9 for the P(VDF-TrFE) material); ε 0 is the electric constant (ε 0 ≈ 8.854 × 10 −12 F/m); h p is the thickness of the P(VDF-TrFE) belt; l p is the length of the P(VDF-TrFE) belt; E p is the Young's Modulus of the P(VDF-TrFE) belt.
The force applied on the belt is where d 31 is the piezoelectric coefficient (6 × 10 −12 C/N) of the P(VDF-TrFE) thin film.
Since the strain generated is hence the maximum bending deformation of the P(VDF-TrFE) belt is  where n stands for the number of mode, together with boundary condition, and C v and C ϕ are the maximum bending deformation and torsional angle respectively. The boundary condition is v = 0, v = 0 at the belt ends, ϕ = 0 at the belt ends.
The numerical results shown in Tab. S3 demonstrate the pair of peak frequencies which are related to the fluttering phenomenon observed in the experiment. It is found that the nonlinear effect of deformation under 0.01mm (observed in the experiment) is ignorable: the torsional frequency is almost a constant within the different deformation range; the bending frequency slightly increases a little. Thus, it is natural that the final fluttering frequency is almost a constant at different incoming air velocity, with considering the fact found that the final fluttering frequency is slightly lower than the torsional frequency. Thus it can be concluded that linear vibration equations Eqs. (S2.3) and (S2.4) are still applicable for the vibration of "stiff" AlN/Si belts vibrating at small deformation in the current experiment.

S4.2 "Flexible" P(VDF-TrFE) belt
In contrast of "stiff" AlN/Si belt, the vibration of "flexible" belt is dominated by the stretching force. So Eqs. (S2.6) and (S2.7) under T 0 = constant = 0.07N is solved with the initial condition of second order perturbation are given as together with the boundary condition v = 0, v = 0 at the belt ends, ϕ = 0 at the belt ends, where n stands for the number of mode, together with boundary condition, and C v and C ϕ are the maximum bending deformation and torsional angle respectively. The numerical results are shown in Tab. S4, and they only demonstrate the one pair of peak frequencies which are related to the fluttering phenomenon observed in the experiment. It is found that the bending deformation of belt about 0.1-0.4 mm (recorded in the experiment) has highly nonlinear effect on both bending and torsional frequencies. Thus, with considering the fact that the final fluttering frequency is determined by the torisional frequency, it is unsurprising that the final fluttering frequency dramatically increases with the increasing incoming air velocity, which is consistent with the experimental observations. It is concluded that for bending and torsional vibration for belts with such large deformation, nonlinear vibration equations Eqs. (S2.6) and (S2.7) shall be used to involve nonlinearity in order to obtain the accurate peak bending and torsional frequencies, then they can be used to calculate the final fluttering frequency based the 2DOF Theodorsen theory (see §S2). Table S4. The mathematic model retrieved natural frequencies of the bending and torsional modes and the corresponding flutter frequency f F for less pre-stressed "flexible" P(VDF-TrFE) micro-belt. If the pre-existing stretching force is increased up to T 0 = constant = 0.21N, then the relative nonlinearity of large deformation becomes much lower than that in the case with smaller stretching force T 0 = constant = 0.07N. Similar numerical simulation was conducted, and the results are listed in Tab. S5. Table S5. The mathematic model retrieved natural frequencies of the bending and torsional modes and the corresponding flutter frequency f F for more pre-stressed "flexible" P(VDF-TrFE) micro-belt.