Energy Harvesting Performance of Plate Wing from Discrete Gust Excitation

Energy harvesting from aeroelastic response tends to have a wide application prospect, especially for small-scale unmanned aerial vehicles. Gusts encountered in flight can be treated as a potential source for sustainable energy supply. The plate model is more likely to describe a low aspect ratio, thin plate wing structure. In this paper, the Von Kármán plate theory and 3D doublet lattice method, coupled with a piezoelectric equation, are used to build a linear state-space equation. Under the load of “one-minus-cosine” discrete gust, the effects of flow speed and gust amplitude, thickness of piezoelectric ceramic transducer (PZTs) layers, and mounted load resistance are investigated. Results reveal that the PZTs layers on the wing root of the leading edge can obtain the highest electrical parameters. The flow velocity, thickness of the PZTs layers and load resistance are used to optimize energy harvesting data.


Introduction
In the last few years, numerous studies have reached a consensus that energy harvesting (EH) from aeroelastic response has a potential application foreground [1][2][3].Aeroelastic flutter characteristics can be a new approach to generate electricity from wind power, or be a direct replenishment to provide energy for an aircraft itself [4].
The expanded direction of EH is to generate wind power by using the flow-induced vibration from wing section or thin plate.It is considered as a substitute for a traditional wind turbine on miniaturization.Mckinney et al. [5] designed a device, which drives the wing section into an aeroelastic flutter under a certain flow speed and transforms plunge motion into rotation to drive an electromagnetic generator.Bryant et al. [6,7] proposed the other way to perform aeroelastic flutter EH by using piezoelectric materials, and researching the base structure character, which replaced the wing section with a rigid plate.Orrego et al. [8] investigated equipment for collecting wind power with a polyvinylidene fluoride (PVDF) membrane, a type of organic piezoelectric material.The adopted inverted-flag configuration can maintain the flutter motion in a wider range of wind speeds and the flutter boundary can be changed by tailoring the length-width ratio.An airfoil with double plunge degrees of freedom was proposed by Wu et al. [9].Wu's model can collect extra EH power from airfoil plunge motion, which can be a complement for the single pitch motion EH system in Reference [1].
With the development of small-scale unmanned aerial vehicles (UAVs), a contradiction exists between weight reduction and the penalty associated with a power source.The duration of power supply, which is limited by the energy density of carried fuel or batteries, restricts the ability to perform tasks.EH from the phenomenon of aeroelastic response was suggested as an additional energy refill method, as was solar radiation.The main aeroelastic responses of the wing structure are flutter, limit cycle oscillations (LCOs) and gust response, although they are considered as negative phenomena since in many cases these unwanted vibrations may cause a reduction in life span and structural damage.
Giacomello et al. [10] used a plate covered with ionic metal polymeric materials (IPMC) for EH of LCOs.Dunnmon et al. [11] placed a piezoelectric material at the root of a plate, and both numerical simulation and experimental results showed that the device can gain about 17% total energy from flow.Tang et al. [12] developed a piezoelectric-aeroelastic coupled system based on the composite material laminate plate and analyzed power efficiency by changing the yawed flow angle.Aeroelastic vibration EH from a cantilevered plate wing with embedded piezoceramics was demonstrated by Marqui et al. [13].The model was based on the frequency-domain and discussed the output electrical feature under the variables of different airflow speed and electrical boundary conditions.Wu et al. [14] built a coupled pitch-plunge-leadlag airfoil-based piezoaeroelastic EH model and compared the power outputs from the plunge and leadlag motions.The influence of nonlinear dynamic behaviors for the proposed EH system was discussed.
Gust response is more likely to be encountered in the practical flight [15].Many studies have been completed on gusts, such as analysis of gust response [16], multi-objective gust optimization [17], and gust alleviation control [18], although these purposes aim to alleviate the influence of gust.From a different point of view, gusts can provide a potential source for EH.Pozzi et al. [19] studied EH from a wing structure under gust loads using engineering softwares Nastran and Ansys.Wu et al. [9] built a coupled piezo-aeroelastic model for the analysis of EH from the gust response and evaluated the energy conversion efficiency and material utilization efficiency.
The aim of this paper is to investigate EH of a plate wing model under the discrete gust load based on previous work by Wu et al. [20].From the results of the beam model EH [20], an aeroelastic response converges rapidly via relevant high bending and twisting stiffness, and the time point of the highest voltage output is around the end of gust encounter phase.The different EH phenomena of the beam model can be found in Bruni et al. [21], for the obvious aeroelastic response behavior, the highest power output (as well as voltage) occurs in the phase of free response phase.EH from gust response behavior of the plate model wing structure was adopted and evaluated as a notable low bending and twisting stiffness model [13].According to the two above-mentioned beam models, plate structure seems to be more of the latter.To analyze its EH characteristics, the proportions of EH output at each phase need to be assessed under certain conditions.Furthermore, the data differences with the location variation of the chord direction can also be calculated using the 2-D plate feature.In this paper, the kinetic equation is established using linear plate theory, unsteady aerodynamics and coupled with the piezoelectric equation.A discrete "one-minus-cosine" gust is loaded as a base excitation.EH properties of the model are studied based on the time-domain state-space equations.The performance of energy output of aeroelastic response is quantified in the voltage, electrical power and energy.For further analysis, the effect of flow speed and gust intensity, piezoelectric ceramic transducer (PZTs) thickness, location and load resistance value are discussed by comparing the data of the energy output density.

Aeroelectroelastic Model
The piezoaeroelastic model is designed as a thin cantilevered plate wing and covered with piezoelectric plates, shown in Figure 1.The piezoelectric materials, covered by continuous and perfectly conductive electrodes with certain thickness, are placed along the elastic axis symmetrically.The transducers and the substructure are assumed to be perfectly bonded together.The wing is in subsonic flow and subject to a vertical discrete gust load.The load resistance is mounted for the charging circuit.Every blocked symmetric PZTs on both sides can be calculated separately.

Structural Model
The structure model is based on the bending motion equation of Von Karman large deflection equations, coupled with piezoelectric layers.The equation of motion of the plate wing is ( , , ) ( , , ) ( , , )

E Mw x y t Cw x y t Kw x y t L
where w is the displacement at plate location ) x, y ( at time t ; M , C and K are the mass, damping and stiffness matrix; E L is the aerodynamic force matrix.
The e-form piezoelectric equations are [22]: where , , T D E and S are the stress vector, electric displacement vector, electric field vector, and strain vector; , c e and ε are the elastic stiffness matrix, piezoelectric stress matrix, and matrix of permittivity components.The superscripts E and S indicate that corresponding parameters are measured at a constant electric field and a constant strain, respectively; the superscript T indicates transposition.
The Rayleigh-Ritz approach is used to discretize the system.The deflection is assumed to be: where φ is the th i order modal function of the piezoelectric cantilevered wing; u is the th i order modal time coordinate.
Combining Equations ( 1)-( 3), a generalized form of the coupled aeroelectroelastic flow-induced vibration based EH model [13] can be expressed as: where U is the vector describing the structural states; Θ is the piezoelectric coefficient matrix; o V is the piezoelectric voltage output vector; p C is the equivalent piezoelectric capacitance matrix; e R is the applied load resistance.The formula derivation can also be found in Reference [20].

Aerodynamic Model
The linear unsteady aerodynamic model is the doublet lattice method (DLM) [23], which is a finite element method for the oscillatory, subsonic, inviscid lifting surface theory.The thin wing can

Structural Model
The structure model is based on the bending motion equation of Von Karman large deflection equations, coupled with piezoelectric layers.The equation of motion of the plate wing is where w is the displacement at plate location (x, y) at time t; M, C and K are the mass, damping and stiffness matrix; L E is the aerodynamic force matrix.The e-form piezoelectric equations are [22]: where T, D, E and S are the stress vector, electric displacement vector, electric field vector, and strain vector; c, e and ε are the elastic stiffness matrix, piezoelectric stress matrix, and matrix of permittivity components.The superscripts E and S indicate that corresponding parameters are measured at a constant electric field and a constant strain, respectively; the superscript T indicates transposition.The Rayleigh-Ritz approach is used to discretize the system.The deflection is assumed to be: where φ is the ith order modal function of the piezoelectric cantilevered wing; u is the ith order modal time coordinate.
Combining Equations ( 1)-( 3), a generalized form of the coupled aeroelectroelastic flow-induced vibration based EH model [13] can be expressed as: where U is the vector describing the structural states; Θ is the piezoelectric coefficient matrix; V o is the piezoelectric voltage output vector; C p is the equivalent piezoelectric capacitance matrix; R e is the applied load resistance.The formula derivation can also be found in Reference [20].

Aerodynamic Model
The linear unsteady aerodynamic model is the doublet lattice method (DLM) [23], which is a finite element method for the oscillatory, subsonic, inviscid lifting surface theory.The thin wing can be assumed as a lift surface and divided into several "box" elements with doublet singularities of constant strength in chordwise and parabolic strength in spanwise direction.These segments are arranged in columns parallel to the freestream, surface edges and fold lines lie on the boundaries.The 1/4 chord line of each box is taken to contain a distribution of acceleration potential doublets of uniform but unknown strength.The control point, arranged at the cross of 3/4 chord and half span lines, is verified boundary condition.
The integral pressure-wash equation for multiple interfering surfaces is: Where w a is the induced downwash at the surface; K is the kernel function for nonplanar surface; N is the number of divided lifting surfaces; ∆p is the dimensionless lifting pressure coefficient; ε and σ are dummy variables of integration over the area S of the wing in spanwise direction y and chordwise direction x.
In matrix notation, Equation ( 5) becomes where D a is the downwash factor matrix related to the kernel function; ρ ∞ is the density of air; v is the airflow speed; w g is the downwash caused by the discrete gust.In general, the gust is represented in discrete or continuous forms.A relative realistic discrete gust model that allows for the spatial scale of the gust zone is the "one-minus-cosine" model.The vertical velocity of the gust v g is [23]: where v g0 is the gust velocity amplitude; H is the generalized wind discrete scale; s is the distance that the wing travels in the gust.Note that this model is based on the frequency-domain procedures, which is not conveniently adaptable to the consideration of aeroelastic gust response.The time-domain equation of motion is more suitable for establishing a liner state-space model for the discrete gust analysis.The minimum state approach can describe the time-domain aerodynamic force and the discrete gust by using a rational function of the Laplace variable s [24], (8) where F(s) is the frequency-domain aerodynamic force matrix; D w , R w , E f , E g are approximation matrix; b h is the half chord length.
Combining Equations ( 5)-( 8), the time-domain aerodynamic lift is approximated to where x a is augmenting aerodynamic state vector.

State-Space Model
Combining Equation ( 4) with the aerodynamic model of Equation ( 9), the aeroelastic model in the time-domain state-space form is obtained, .
where M, C, K are the combined mass, damping and stiffness matrix obtained from Equation (4);

Performance Analysis
The power output is where V o1 and V o2 are the voltage from upper and lower PZTs.
The electric energy harvested is where T is the time span for the energy output.
The energy output density is introduced to evaluate the performance of the harvester, namely, the energy output per unit volume of the PZTs, where V c is the volume of PZTs layers.

Results and Discussions
The initial geometric and material properties of the plate wing and PZTs layers are given in Tables 1 and 2. The wing structure is made of aluminum Al 7075, and the PZTs material is PZT-5H.The operating conditions for the subsequent results are as follows; cruise speed of 40 m/s, 1500 m above sea level, the gust shape was for a discrete "one-minus-cosine" function with peak intensity of 5 m/s, a gust wavelength of 100 m, and the gust interacts with the wing at the initial time.The displacement is measured at the leading edge of the wing tip, where the EH data was collected from a pair of PZTs layers at the leading edge of the wing root for the highest energy output density.
The range of gust intensity was chosen from ESDU 04024 [25].For 40 m/s cruise speed, 1500 m altitude and two hours flight, the largest gust intensity occurred more than once is 5 m/s.The comparison of gust amplitude is from a higher probability condition 3 m/s (more than twice during a flight) to an extreme condition 10.50 m/s (2% probability during a flight).
The flutter boundary should be determined initially as an important aeroelastic behavior.The linear flutter speed can be obtained from the state-space matrices (Equation ( 10)) using linear eigenvalue analysis.The stability and frequency of the system are determined by the real and imaginary parts of the eigenvalues.The flutter speed can be calculated when the real part of eigenvalues of the aeroelectroelastic system become positive.
The flutter characteristics of the wing model adopted is shown in terms of damping and frequency at an increasing airflow speed in Figure 2. The load resistance of EH system is 5.00 × 10 4 Ω.Results show that the real part of the eigenvalues first become positive at an airflow speed of 154 m/s and that the critical mode is a coupled 1 bending-1 torsional mode, as shown in Figure 2a. Figure 2b shows the frequency behaviors of the first three modes.The first six modes are listed in Table 3.It is noticed that all the modes present bending-torsional coupling for the plate model.Results show that the real part of the eigenvalues first become positive at an airflow speed of 154 m/s and that the critical mode is a coupled 1 bending-1 torsional mode, as shown in Figure 2a. Figure 2b shows the frequency behaviors of the first three modes.The first six modes are listed in Table 3.It is noticed that all the modes present bending-torsional coupling for the plate model.A preliminary study was carried out to investigate the number of mode shapes needed for converged results.The energy output density is chosen as a calibration criterion, which can be calculated from the combination of Equations ( 3), ( 10)-( 13).Figure 3 shows the energy output density of the PZTs plates on both sides of the leading edge of the wing root.

Time-History Analysis
The time-history analysis of the gust response can be calculated from the state-space matrices (Equation ( 10)).The Fourth Order Runge-Kutta method is used to solve the motion equation of the system.
Results are presented in Figure 4. Figure 4a shows the tip displacement at the leading edge of the wing.Figure 4b shows the voltage output at the leading edge symmetric PZTs layers of the wing root, which has the highest EH output data.Figure 4c demonstrates the power output at the leading edge symmetric PZTs layers of the wing root, obtained from Equation (11).
The gust load lasts for 0.05 s in this condition, the gust duration is 2.50 s and the first 2.50-second response is selected of the voltage and power output.The detail EH energy data of every PZTs layers is shown in Figure 5, with the gust encounter presented in Figure 5a and free response power output in Figure 5b.The displacement of wing tip continues to increase until the peak 4.40 × 10 −2 m at 2.50 s during the end of gust wavelength.After the gust moves away from the wing the voltage and power output reach their peaks, −6.44 V and 8.29 × 10 −4 W when the tip displacement springs back to −2.43 × 10 −2 m at 2.58 s.For all PZTs layers, EH power during the first second of gust response is 14.52 × 10 −4 J, 95.68% of the total EH (15.18 × 10 −4 J).During the gust interaction, the data are 1.30 × 10 −5 J, 0.86%.
Results demonstrate that the voltage and power output have phase differences with the tip displacement, yet it is consistent with the stress changing.Comparing with the percent of EH data, EH is mainly from the first few seconds (from 2.50 s in Figure 4) when the wing disengaged from the Results show that the relative error between EH data of 5-bending coupled 4-twisting modes and the data of 5-bending coupled 3-twisting modes is 0.73%.Therefore, the results presented include the lowest 5 bending and 4 twisting modes.

Time-History Analysis
The time-history analysis of the gust response can be calculated from the state-space matrices (Equation ( 10)).The Fourth Order Runge-Kutta method is used to solve the motion equation of the system.
Results are presented in Figure 4. Figure 4a shows the tip displacement at the leading edge of the wing.Figure 4b shows the voltage output at the leading edge symmetric PZTs layers of the wing root, which has the highest EH output data.Figure 4c demonstrates the power output at the leading edge symmetric PZTs layers of the wing root, obtained from Equation (11).The gust load lasts for 0.05 s in this condition, the gust duration is 2.50 s and the first 2.50-second response is selected of the voltage and power output.The detail EH energy data of every PZTs layers is shown in Figure 5, with the gust encounter presented in Figure 5a and free response power output in Figure 5b.The displacement of wing tip continues to increase until the peak 4.40 × 10 −2 m at 2.50 s during the end of gust wavelength.After the gust moves away from the wing the voltage and power output reach their peaks, −6.44 V and 8.29 × 10 −4 W when the tip displacement springs back to −2.43 × 10 −2 m at 2.58 s.For all PZTs layers, EH power during the first second of gust response is 14.52 × 10 −4 J, 95.68% of the total EH (15.18 × 10 −4 J).During the gust interaction, the data are 1.30 × 10 −5 J, 0.86%.Results demonstrate that the voltage and power output have phase differences with the tip displacement, yet it is consistent with the stress changing.Comparing with the percent of EH data, EH is mainly from the first few seconds (from 2.50 s in Figure 4) when the wing disengaged from the gust zone.It can be considered that the main influence of the gust load is the modal displacement, modal velocity, augmenting aerodynamic state vector and voltage of structure and PZTs layers' circuit system at the end of gust zone, as shown in Equation (10), which is also the initial condition at the beginning of free response phase.
Energy output differences of wing-like EH system in the chordwise direction PZTs layers need to be concerned.In the beam structure, it is neglected for the model limitation.As a 2-D model, the plate structure can demonstrate the energy output differences in the chordwise direction.The reason for the EH inconsistency in Figure 5 is that the differences are balanced by cycles of the torsional vibration (Figure 5b), unlike the disparity caused by a single one-direction displacement in gust encounter phase (Figure 5a).
The main differences of EH characteristic between the plate wing and the beam structure results of Reference [19] is the structure behavior.The torsion and bending stiffness of the plate model wing is weaker than the beam structure wing with the same wing area.Thus, the proportion of EH data of plate wing free response in whole phase is more than the beam wing at the same scale of cruise speed and gust amplitude.

Effect of Gust Velocity Amplitude and Airflow Velocity of the Energy Output Density
Under the circumstances of various flow velocity and gust amplitude, the changing of energy output density is presented in Figure 6, obtained from Equations (10) and (13).
Results show that both flow velocity and gust amplitudes exert an important part in EH.The slope of curves reveals the influence of gust intensities, and the EH density increases significantly which is impacted by the rising of both flow velocity and gust amplitude.For encountering gust, more severe conditions, discussed above, will bring a better performance of EH efficiency.The increase of gust intensity is more efficiency than the flow velocity's, the latter is limited by the safety margin of flutter speed.

Effect of Thickness and Location of the Energy Output Density
The energy output value of every PZTs layers at different locations can be obtained from the solution of Equations ( 10) and ( 13).The thickness range of PZTs layers is from 1.00 × 10 −5 m to 1.00 × 10 −3 mm.The relationship between the energy output density and the location of the PZTs is given in Figure 7.A similar variation of the output density with the location is observed.The first four pairs of layers from the wing root at leading edge are chosen for the higher energy output density.With the rising of thickness, the energy output density peaks when the abscissa value reaches 6.00 × 10 −5 m, followed by a rapid decrease.The tendency above happens to all listed PZTs layers at different locations, which also has a great influence on EH in a certain thickness range, as indicated in Figure 7.
It should be emphasized that the analysis above aims at one specific "one-minus-cosine" discrete gust load.Once the gust intensity or wind scale changes, the dynamic response of the wing would change.Therefore, the position which is beneficial to EH may be different accordingly.Generally, the amplitude and scale of the gust loads are not constant, which makes it difficult to find the optimum location.However, the parametric analysis could still be used to identify better designs.
The results indicate that the impact of PZTs thickness need to be considered even at the certain part of maximum stress.A reasonable PZTs laminate structure will improve EH efficiency.Results show that both flow velocity and gust amplitudes exert an important part in EH.The slope of curves reveals the influence of gust intensities, and the EH density increases significantly which is impacted by the rising of both flow velocity and gust amplitude.For encountering gust, more severe conditions, discussed above, will bring a better performance of EH efficiency.The increase of gust intensity is more efficiency than the flow velocity's, the latter is limited by the safety margin of flutter speed.

Effect of Thickness and Location of the Energy Output Density
The energy output value of every PZTs layers at different locations can be obtained from the solution of Equations ( 10) and (13).The thickness range of PZTs layers is from 1.00 × 10 −5 m to 1.00 × 10 −3 mm.The relationship between the energy output density and the location of the PZTs is given in Figure 7.A similar variation of the output density with the location is observed.The first four pairs of layers from the wing root at leading edge are chosen for the higher energy output density.With the rising of thickness, the energy output density peaks when the abscissa value reaches 6.00 × 10 −5 m, followed by a rapid decrease.The tendency above happens to all listed PZTs layers at different locations, which also has a great influence on EH in a certain thickness range, as indicated in Figure 7.

Effect of Load Resistance of the Energy Output Density
The effect of the external load resistance is shown in Figure 8, obtained from Equations ( 10) and ( 13).The energy output density increases with the load resistance at the beginning and reaches a peak value of 246.83 J m −3 at 5.00 × 10 4 Ω.
A steady fall of the conversion efficiency follows with the further increase of the resistance.As It should be emphasized that the analysis above aims at one specific "one-minus-cosine" discrete gust load.Once the gust intensity or wind scale changes, the dynamic response of the wing would change.Therefore, the position which is beneficial to EH may be different accordingly.Generally, the amplitude and scale of the gust loads are not constant, which makes it difficult to find the optimum location.However, the parametric analysis could still be used to identify better designs.
The results indicate that the impact of PZTs thickness need to be considered even at the certain part of maximum stress.A reasonable PZTs laminate structure will improve EH efficiency.

Effect of Load Resistance of the Energy Output Density
The effect of the external load resistance is shown in Figure 8, obtained from Equations ( 10) and ( 13).The energy output density increases with the load resistance at the beginning and reaches a peak value of 246.83 J m −3 at 5.00 × 10 4 Ω.

Effect of Load Resistance of the Energy Output Density
The effect of the external load resistance is shown in Figure 8, obtained from Equations ( 10) and (13).The energy output density increases with the load resistance at the beginning and reaches a peak value of 246.83 J m −3 at 5.00 × 10 4 Ω.
A steady fall of the conversion efficiency follows with the further increase of the resistance.As the resistance approaches infinity, the voltage output tends to be a constant value, i.e., the opencircuit voltage.Thus, the electric power output tends to be zero.Thus, the energy output density has similar trends with the variation of the load resistance.
The determination of optimal resistance is significant for the design of more complex external circuit, such as one capable of energy collection and storage.

Conclusions
A time-domain state-space aeroelastic model coupled with piezoelectric equations is calculated for gust response analysis.The base excitation of the EH system is "one-minus-cosine" discrete gusts.Initial working conditions and the density of energy as a non-dimensional parameter are defined to evaluate EH performance.Results show that output data are influenced by various factors; PZTs at A steady fall of the conversion efficiency follows with the further increase of the resistance.As the resistance approaches infinity, the voltage output tends to be a constant value, i.e., the open-circuit voltage.Thus, the electric power output tends to be zero.Thus, the energy output density has similar trends with the variation of the load resistance.
The determination of optimal resistance is significant for the design of more complex external circuit, such as one capable of energy collection and storage.

Conclusions
A time-domain state-space aeroelastic model coupled with piezoelectric equations is calculated for gust response analysis.The base excitation of the EH system is "one-minus-cosine" discrete gusts.Initial working conditions and the density of energy as a non-dimensional parameter are defined to evaluate EH performance.Results show that output data are influenced by various factors; PZTs at the wing root of leading edge can generate maximum electrical related parameters under the same condition, the phase difference exists between tip displacement and electrical output, EH efficiency is more significant for a flow speed close to the flutter boundary, and a certain thickness range of PZTs can provide the highest energy output density, as well as load resistance.The electrical output of EH from a single gust reaches grades of milliwatt and millijoule, which can be a potential source for powering micro-electronic systems mounted on small-scale UAVs.

Aerospace 2018, 5 , 13 Figure 2 .
Figure 2. First three bending-twisting modes of the wing covered by PZTs layers with 50,000 Ω load resistance.(a) Damping evolution with increasing flow velocity; (b) frequency evolution with increasing flow velocity.

Figure 2 .
Figure 2. First three bending-twisting modes of the wing covered by PZTs layers with 50,000 Ω load resistance.(a) Damping evolution with increasing flow velocity; (b) frequency evolution with increasing flow velocity.

13 Figure 3 .
Figure 3. Energy output with bending-twisting modes in the Rayleigh-Ritz method.

Figure 3 .
Figure 3. Energy output with bending-twisting modes in the Rayleigh-Ritz method.

Figure 4 .
Figure 4. Time histories of (a) displacement at leading edge of wing tip; (b) voltage output at the leading edge symmetric PZTs layers of wing root; (c) power output at the leading edge symmetric PZTs layers of wing root.Flight working condition: 40 m/s chordwise flow speed; 5 m/s "one-minuscosine" discrete gust intensity; 100 m gust wavelength.

Figure 5 .
Figure 5. EH power output during (a) gust encounter phase; (b) EH power output during free response phase.Flight working condition: 40 m/s chordwise flow speed and 5 m/s "one-minuscosine" discrete gust intensity; 100 m gust wavelength; Wing span length from the root (0 m) to the tip (0.80 m); chord length from leading edge (0 m) to trailing edge (0.20 m).

Figure 4 .
Figure 4. Time histories of (a) displacement at leading edge of wing tip; (b) voltage output at the leading edge symmetric PZTs layers of wing root; (c) power output at the leading edge symmetric PZTs layers of wing root.Flight working condition: 40 m/s chordwise flow speed; 5 m/s "one-minus-cosine" discrete gust intensity; 100 m gust wavelength.

Figure 4 .
Figure 4. Time histories of (a) displacement at leading edge of wing tip; (b) voltage output at the leading edge symmetric PZTs layers of wing root; (c) power output at the leading edge symmetric PZTs layers of wing root.Flight working condition: 40 m/s chordwise flow speed; 5 m/s "one-minuscosine" discrete gust intensity; 100 m gust wavelength.

Figure 5 .
Figure 5. EH power output during (a) gust encounter phase; (b) EH power output during free response phase.Flight working condition: 40 m/s chordwise flow speed and 5 m/s "one-minuscosine" discrete gust intensity; 100 m gust wavelength; Wing span length from the root (0 m) to the tip (0.80 m); chord length from leading edge (0 m) to trailing edge (0.20 m).

Figure 5 .
Figure 5. EH power output during (a) gust encounter phase; (b) EH power output during free response phase.Flight working condition: 40 m/s chordwise flow speed and 5 m/s "one-minus-cosine" discrete gust intensity; 100 m gust wavelength; Wing span length from the root (0 m) to the tip (0.80 m); chord length from leading edge (0 m) to trailing edge (0.20 m).

Figure 6 .
Figure 6.Energy output density with gust velocity amplitude under different airflow velocity (v) at the leading edge symmetric PZTs layers of wing root, 100 m gust wavelength.

Figure 6 .
Figure 6.Energy output density with gust velocity amplitude under different airflow velocity (v) at the leading edge symmetric PZTs layers of wing root, 100 m gust wavelength.

13 Figure 7 .
Figure 7. Energy output density with thickness under different locations of wing span (L) at the leading edge symmetric PZTs layers from wing root.

Figure 7 .
Figure 7. Energy output density with thickness under different locations of wing span (L) at the leading edge symmetric PZTs layers from wing root.

Figure 7 .
Figure 7. Energy output density with thickness under different locations of wing span (L) at the leading edge symmetric PZTs layers from wing root.

Figure 8 .
Figure 8. Energy output density with load resistance at the leading edge symmetric PZTs layers of wing root.

Figure 8 .
Figure 8. Energy output density with load resistance at the leading edge symmetric PZTs layers of wing root.

Table 1 .
The geometric and material parameters of the wing.

Table 2 .
Geometric and material parameters of the PZT-5H layers.

Table 3 .
First six bending-twisting mode shapes and natural frequencies of the PZTs covered plate

Table 3 .
First six bending-twisting mode shapes and natural frequencies of the PZTs covered plate wing.