An Analytical–Numerical Method for Simulating the Performance of Piezoelectric Harvesters Mounted on Wing Slats

Recently, there have been significant developments in the field of vibration energy harvesters to feed sensors for structural health monitoring in aeronautics and other high technology fields. Within the framework of the EU InComEss project, new eco-friendly piezoelectric materials are under development. A foreseen application is vibration energy harvesting from the wing slats of aircraft. Semi-analytical models of the vibrating slat make it possible to estimate the maximum voltage that can be generated by a piezoelectric patch bonded to a slat surface. A more detailed analysis must consider details of the three-dimensional geometry of both the harvester and the bonding layer. This can only be carried out with multiphysics finite element software. A finite element model of a whole slat would require a large computational effort as millions of elements are typically needed to model very thin piezoelectric layers. To simplify this analysis, an integrated analytical–numerical method is proposed in this paper. A large-scale analytical model of the whole slat was used to calculate loads on the portion of the slat where a piezoelectric patch was attached. Then, a small-scale finite element model of the portion of the slat with the piezoelectric patch was used to calculate the open circuit voltage generated by the patch. The response of the harvester to random excitation, typical of aeronautic applications, was calculated. The effects of the details of the harvester design on the generated voltage were analyzed and discussed.


Introduction
Vibration energy harvesting by means of piezoelectric (PE) materials is a promising technology for feeding remote sensor nodes and other microelectronic equipment, especially when weight, space and accessibility constraints are important. In aeronautics, several applications of vibration energy harvesting have been proposed [1,2], but the most important application is feeding sensors for structural health monitoring (SHM) [3][4][5]. In the framework of the EU InComEss project [6], new lead-free piezoelectric materials are under development. A foreseen application is vibration energy harvesting from wing slats. Slats are movable aerodynamic surfaces on the leading edge of the wing that increase the lift force on the wing when they are deployed. In [7], an analytical model based on the modal superposition approach [8] was used to simulate PE harvesters mounted on a deployed slat. A comparison was made between the performance of two possible design solutions: cantilever harvesters and PE patches directly bonded to the slat surface. Numerical results showed that both solutions were able to generate relevant voltage and power. The cantilever harvester optimally tuned to the most excited mode of vibration of the slat showed the best performance, since it exploited the resonance phenomenon, but this device required more mass, volume and stress inside the brittle piezoelectric material than the PE patch. Because the mass, volume and reliability of PE devices are very important issues in slat case-studies, the PE patch was selected for further development within the framework of the InComEss project. The deployed slat in [7] showed the largest vibrations; however, most of the energy could be harvested during the flight when the slat was retracted. The analytical model presented in [7] was extended in this study to investigate this operating condition. Moreover, to analyze the effect of details such as the three-dimensional geometry of the PE patch and thin adhesive layers used to bind the PE patch to the slat, an integrated approach combining a large-scale analytical model and a small-scale finite element model is proposed.
The large-scale model, implemented in MATLAB ® , is fully analytical and is based on the modal expansion approach. It describes a whole slat excited by the broadband acceleration spectrum typical of wings, and it enables the calculation of the bending moment and shear force acting on the portion of the slat where the harvester is mounted. The contact between the slat and the wing edge is simulated by means of a distributed stiffness, and the variation in the value of the contact stiffness makes it possible to simulate both retracted and deployed slats.
The small-scale model, implemented in COMSOL Multiphysics ® , is numerical and is based on a multiphysics finite element (FE) method. It describes the sandwich structure of the harvester mounted on an equivalent portion of the slat, which is excited by the loads calculated by means of the large-scale model. This makes it possible to compute the voltage generated by the PE material considering the effects of the adhesive layer and of the other layers that compose the sandwich.
The paper is organized as follows: Section 2 describes the problem and the structure of the combined analytical-numerical approach. Section 3 provides the input data: slat dimensions, PE patch dimensions, mechanical properties, and electrical properties. Section 4 deals with the large-scale analytical model of the whole slat. Section 5 deals with the small-scale FE model. Section 6 shows the numerical results and, in particular, the effects of the adhesive thickness on the voltage generated by the harvester. Finally, conclusions are drawn.

Integrated Analytical-Numerical Method
The slat is the mobile leading-edge flap of a wing of an aircraft, which is operated to increase the angle of attack of the wing during low-speed maneuvers, such as take-off and landing. The slats are retracted during flight and are forced by the aerodynamic loads to remain in contact with the leading edge of the wing. The deployment and retraction of the slat is performed by servomechanisms. In the framework of this research, it was assumed that the slat is moved by two servomechanisms, so the whole deployed slat could be schematized as a pinned beam with overhangs, as shown in Figure 1a. A distributed stiffness was added to allow for the contact force between the wing and the retracted slat. Figure 1b shows the scheme of the retracted slat. The slat vibrates since it is excited by the wing through the slat supports. The vibration levels of the supports are equal to the vibration levels of the wing's leading edge corresponding with the connecting points. The slat is not a rigid body; hence, it deforms due to inertia force. Figure 2 shows that the vibration levels at the supports are different, since vibrations along the wing of the aircraft increase from the root towards the tip of the wing [7]. In Figure 2, w r (x, t) represents the displacement of the deformed slat with respect to the undeformed configuration and a(t) is the acceleration level of one of the two supports of the slat. Actuators 2023, 12  From a physical point of view, the slat can be modelled as a pinned beam excited by a trapezoidal distribution of forces (as proposed in [7]) which represents the inertia force due to the vibration of the wing; see Figure 3. In aeronautics, acceleration levels are given in the frequency domain according to standard specifications. In this research, the standard specification RTCA-DO-160 CAT S curve E was adopted because it refers to wing components. Figure 4 shows the PSD of the acceleration of the external support of the slat. Therefore, a frequency domain analysis was needed to calculate the power spectral density (PSD) of the voltage generated by a PE patch mounted on the slat from the acceleration PSD, and then the voltage RMS value could be obtained with Parseval's Theorem.  From a physical point of view, the slat can be modelled as a pinned beam excited by a trapezoidal distribution of forces (as proposed in [7]) which represents the inertia force due to the vibration of the wing; see Figure 3. In aeronautics, acceleration levels are given in the frequency domain according to standard specifications. In this research, the standard specification RTCA-DO-160 CAT S curve E was adopted because it refers to wing components. Figure 4 shows the PSD of the acceleration of the external support of the slat. Therefore, a frequency domain analysis was needed to calculate the power spectral density (PSD) of the voltage generated by a PE patch mounted on the slat from the acceleration PSD, and then the voltage RMS value could be obtained with Parseval's Theorem. From a physical point of view, the slat can be modelled as a pinned beam excited by a trapezoidal distribution of forces (as proposed in [7]) which represents the inertia force due to the vibration of the wing; see Figure 3.  From a physical point of view, the slat can be modelled as a pinned beam excited by a trapezoidal distribution of forces (as proposed in [7]) which represents the inertia force due to the vibration of the wing; see Figure 3. In aeronautics, acceleration levels are given in the frequency domain according to standard specifications. In this research, the standard specification RTCA-DO-160 CAT S curve E was adopted because it refers to wing components. Figure 4 shows the PSD of the acceleration of the external support of the slat. Therefore, a frequency domain analysis was needed to calculate the power spectral density (PSD) of the voltage generated by a PE patch mounted on the slat from the acceleration PSD, and then the voltage RMS value could be obtained with Parseval's Theorem. In aeronautics, acceleration levels are given in the frequency domain according to standard specifications. In this research, the standard specification RTCA-DO-160 CAT S curve E was adopted because it refers to wing components. Figure 4 shows the PSD of the acceleration of the external support of the slat. Therefore, a frequency domain analysis was needed to calculate the power spectral density (PSD) of the voltage generated by a PE patch mounted on the slat from the acceleration PSD, and then the voltage RMS value could be obtained with Parseval's Theorem. The piezoelectric patch is very thin and flexible and covers a very small part of the whole slat. Therefore, the vibrations of the slat caused by the motion of the supports are not influenced by the presence of the patch. Conversely, the deformation of the portion of the slat skin where the patch is attached determines the strain inside the PE patch and the generated voltage. The adhesive and protective layers that are inserted between the slat skin and the PE patch may also influence the generated voltage. For the above-mentioned reasons, a computational method based on the synergic use of a large-scale model and a small-scale model was developed. Figure 5 shows a flow chart that describes the steps and the interactions between the large-scale and small-scale models.  The piezoelectric patch is very thin and flexible and covers a very small part of the whole slat. Therefore, the vibrations of the slat caused by the motion of the supports are not influenced by the presence of the patch. Conversely, the deformation of the portion of the slat skin where the patch is attached determines the strain inside the PE patch and the generated voltage. The adhesive and protective layers that are inserted between the slat skin and the PE patch may also influence the generated voltage. For the above-mentioned reasons, a computational method based on the synergic use of a large-scale model and a small-scale model was developed. Figure 5 shows a flow chart that describes the steps and the interactions between the large-scale and small-scale models.  The piezoelectric patch is very thin and flexible and covers a very small part of the whole slat. Therefore, the vibrations of the slat caused by the motion of the supports are not influenced by the presence of the patch. Conversely, the deformation of the portion of the slat skin where the patch is attached determines the strain inside the PE patch and the generated voltage. The adhesive and protective layers that are inserted between the slat skin and the PE patch may also influence the generated voltage. For the above-mentioned reasons, a computational method based on the synergic use of a large-scale model and a small-scale model was developed. Figure 5 shows a flow chart that describes the steps and the interactions between the large-scale and small-scale models.  The large-scale model refers to the whole slat, is based on the modal expansion approach, and is implemented in MATLAB. The aim of this model is to calculate the frequency response functions (FRFs) of the points along the slat caused by a harmonic distributed inertia force having the shape depicted in Figure 3. In particular, the FRFs between bending moment, shear force, and slat acceleration are obtained. They enable the calculation of the forces and moments acting on every portion of the slat.
The small-scale model includes the PE patch and the portion of the slat where the patch is attached. It is developed with the FE method and implemented in COMSOL. The actual cross-section of the slat, which is rather complex, is transformed into an equivalent cross-section that generates in the slat skin the same strain pattern as the actual cross-section when the same loads are present. The small-scale FE is used to calculate, for a series of frequencies, the harmonic response of the system stimulated by the inertia distributed force and by the moments and forces deriving from the large-scale model at the assigned frequency. The fitting of the results of these analyses gives the numerical FRF (FRF v ( f )) between slat acceleration and the open circuit voltage generated by the PE patch. The FRF v ( f ) relates the PSD of the OCV to the PSD of the applied acceleration: The open circuit voltage RMS value can be obtained from Parseval's Theorem as The OCV is an important merit figures since it allows for the evaluation of the electric performance of the harvester, as discussed in [9].

Input Data: Slat and PE Patch Properties
The slat considered in this paper is built by the aerospace company Sonaca Group. It is made of a composite material and has the cross-section shown in Figure 6. The properties of the material and the cross-section, provided by Sonaca, are summarized in Table 1. The large-scale model refers to the whole slat, is based on the modal expansion approach, and is implemented in MATLAB. The aim of this model is to calculate the frequency response functions (FRFs) of the points along the slat caused by a harmonic distributed inertia force having the shape depicted in Figure 3. In particular, the FRFs between bending moment, shear force, and slat acceleration are obtained. They enable the calculation of the forces and moments acting on every portion of the slat.
The small-scale model includes the PE patch and the portion of the slat where the patch is attached. It is developed with the FE method and implemented in COMSOL. The actual cross-section of the slat, which is rather complex, is transformed into an equivalent cross-section that generates in the slat skin the same strain pattern as the actual crosssection when the same loads are present. The small-scale FE is used to calculate, for a series of frequencies, the harmonic response of the system stimulated by the inertia distributed force and by the moments and forces deriving from the large-scale model at the assigned frequency. The fitting of the results of these analyses gives the numerical FRF ( ( )) between slat acceleration and the open circuit voltage generated by the PE patch. The ( ) relates the PSD of the OCV to the PSD of the applied acceleration: The open circuit voltage RMS value can be obtained from Parseval's Theorem as The OCV is an important merit figures since it allows for the evaluation of the electric performance of the harvester, as discussed in [9].

Input Data: Slat and PE Patch Properties
The slat considered in this paper is built by the aerospace company Sonaca Group. It is made of a composite material and has the cross-section shown in Figure 6. The properties of the material and the cross-section, provided by Sonaca, are summarized in Table 1.   A macro fiber composite (MFC) piezoelectric patch manufactured by Smart Material GmbH (patch M8514-P2 in the datasheet of [10]) and directly bonded to the slat skin was considered here. Figure 6 shows the slat, the PE patch, and the global reference frame. The y-axis is aligned to the neutral axis of the slat cross-section. The center of the patch is located at x pc = 1.45 m and z n = 0.025 m. Figure 7 represents the scheme of the M8514-P2 patch and shows the local reference frame and the main dimensions of the patch. The patch is designed to exploit the strain along its longitudinal direction (axis 1), and it is poled along the direction perpendicular to its middle plane (axis 3). A macro fiber composite (MFC) piezoelectric patch manufactured by Smart Material GmbH (patch M8514-P2 in the datasheet of [10]) and directly bonded to the slat skin was considered here. Figure 6 shows the slat, the PE patch, and the global reference frame. The y-axis is aligned to the neutral axis of the slat cross-section. The center of the patch is located at 1.45 m and 0.025 m. Figure 7 represents the scheme of the M8514-P2 patch and shows the local reference frame and the main dimensions of the patch. The patch is designed to exploit the strain along its longitudinal direction (axis 1), and it is poled along the direction perpendicular to its middle plane (axis 3).  Table 2 shows the values of the geometric parameters depicted in Figure 7. To exploit the largest strain values, the PE patch is oriented to align its 1-axis with the x-axis of the global reference frame of the slat. Indeed, the slat vibrates in the z-direction due to the transverse load distribution; hence, the axial strain caused by the bending moment distribution is the main contribution.
In Table 3, the electromechanical properties of the considered PE patch are presented.

Large-Scale Analytical Model of the Slat
The variable inertia force that excites the slat can be expressed as  Table 2 shows the values of the geometric parameters depicted in Figure 7. To exploit the largest strain values, the PE patch is oriented to align its 1-axis with the x-axis of the global reference frame of the slat. Indeed, the slat vibrates in the z-direction due to the transverse load distribution; hence, the axial strain caused by the bending moment distribution is the main contribution.
In Table 3, the electromechanical properties of the considered PE patch are presented. Table 3. Electromechanical properties of the PE patch.

Large-Scale Analytical Model of the Slat
The variable inertia force that excites the slat can be expressed as where g(x) is a shape function that correlates the acceleration of any point of the slat with the acceleration a(t) of one of its supports (the external in the present study); see Figure 2. As demonstrated in [7], the shape function is where α is a parameter related to the slope of the trapezoidal distribution and L 1 and L 2 are the length of the spans of the slat; see Figure 1a. Here, the parameter α was set to 0.05, while L 1 = 0.81 m and L 2 = 1.38 m [7]. The vibrating slat is considered to be a Euler-Bernoulli beam with constant properties: modulus of elasticity E, cross-section moment of inertia I, mass density ρ, cross-section area A, and coefficient of strain-rate damping c s (a proportional damping is assumed). The equation of the forced vibrations of each of the three spans in Figure 1a where k x is the distributed stiffness. The backward electromechanical coupling term [11] is not considered in Equation (5) since the length of the PE patch is much smaller than the overall slat length and because the PE patch is much thinner than the slat. The first steps of the modal expansion approach are the calculations of the natural frequencies and modes of vibration of the undamped system.
The free undamped vibrations for the retracted slat are described by the following homogeneous equation of motion: The relative displacement can be expressed as the product of a spatial function ψ k , which represents the deformed shape of the slat, and a time function f [8]: By inserting Equation (7) in Equation (6) and separating the variables, it becomes EI ρA The right-hand side of Equation (8) has the dimension of a squared angular frequency: After some algebraic manipulations, the partial differential equation for the spatial component ψ k of w r,k is obtained: where includes the unknown natural angular frequencies. The solution of differential Equation (10) for each of three spans of the slat (k = 1, 2, 3) in Figure 1a takes the form: where A k , B k , C k and D k are unknown coefficients. By imposing the boundary conditions on the three spans of the slat and solving the corresponding eigenvalue problem, as discussed in [7], an infinite discrete set of values γ i is found and the corresponding values (A k,i , B k,i , C k,i , and D k,i ) of the unknown coefficients are calculated. From Equation (12), the ith mode of vibration of the slat becomes The corresponding natural frequency f n,i of the ith mode of vibration of the slat can be calculated with Equation (11) as follows: Equation (13) shows that the slat has the same modes of vibrations when deployed or retracted. On the other hand, Equation (14) shows that the modes are characterized by different natural frequencies due to the dependence on k x . Therefore, the interaction between the retracted slat and the wing determines an increase in the natural frequency of each mode. Figure 8 represents the natural frequencies of the first five modes of vibration of the retracted slat as a function of the distributed stiffness k x considering the mechanical properties of the composite slat shown in Table 1. includes the unknown natural angular frequencies. The solution of differential Equation (10) for each of three spans of the slat ( 1,2,3) in Figure 1a takes the form: where , , and are unknown coefficients. By imposing the boundary conditions on the three spans of the slat and solving the corresponding eigenvalue problem, as discussed in [7], an infinite discrete set of values is found and the corresponding values ( , , , , , , anⅆ , ) of the unknown coefficients are calculated. From Equation (12), the ith mode of vibration of the slat becomes The corresponding natural frequency , of the ith mode of vibration of the slat can be calculated with Equation (11) as follows: Equation (13) shows that the slat has the same modes of vibrations when deployed or retracted. On the other hand, Equation (14) shows that the modes are characterized by different natural frequencies due to the dependence on . Therefore, the interaction between the retracted slat and the wing determines an increase in the natural frequency of each mode. Figure 8 represents the natural frequencies of the first five modes of vibration of the retracted slat as a function of the distributed stiffness considering the mechanical properties of the composite slat shown in Table 1.  Figure 8 shows that when the distributed stiffness increases, all the natural frequencies of the modes converge to the same value, which means that the dynamics of the slat are close to those of a rigid body mounted on the distributed stiffness . On the contrary, when the distributed stiffness is much smaller than the flexural stiffness of the slat, it has a negligible effect on the natural frequencies.  Figure 8 shows that when the distributed stiffness increases, all the natural frequencies of the modes converge to the same value, which means that the dynamics of the slat are close to those of a rigid body mounted on the distributed stiffness k x . On the contrary, when the distributed stiffness is much smaller than the flexural stiffness of the slat, it has a negligible effect on the natural frequencies. The relative displacement w r (x, t) of any point of the slat can be expressed by using the modal expansion approach as a linear combination of the modes of vibration of the slat [8,11,12]: where φ i (x) is the ith mass-normalized mode of vibration (that derives from the nonnormalized mode [7]) and η i (t) is the ith modal coordinate. Hence, the equation of motion of the ith modal coordinate is obtained by introducing Equation (15) in Equation (5) and applying the orthogonality conditions [9,11]: ..
where f n,i and ζ i are the natural frequency and the damping ratio of the ith mode of vibration, respectively. In the frequency domain, the modal coordinates η i (t) and the acceleration a(t) are expressed as The frequency response function (FRF) between the modal coordinate amplitude η 0,i and the amplitude of the acceleration a 0 is Therefore, only considering N modes of vibration of the slat, if Equation (18) is inserted into Equation (15), the FRF between the amplitude of the relative displacement of any point of the slat and the acceleration a 0 can be calculated as follows: The bending moment and the shear force along the slat are related to the relative displacement w r (x, t) by the following equations: With harmonic excitation, the bending moment and shear force can be described in the frequency domain using the displacement FRF: Finally, the bending moment and shear force FRFs are defined as (25) Figure 9 shows the magnitude of the FRFs related to the bending moment and shear force for different values of the distributed stiffness calculated at the center of the PE patch (x pc ). The magnitude of the bending moment FRF in Figure 9 is characterized by peaks at the natural frequencies of the odd modes, whereas the magnitude of the shear force FRF is characterized by peaks at the natural frequencies of the even modes. Indeed, the bending moment distribution corresponds to symmetric loads acting on a portion of the slat, so it only excites the odd modes, whereas the shear force corresponds to an anti-symmetric load and only excites the even modes. (25) Figure 9 shows the magnitude of the FRFs related to the bending moment and shear force for different values of the distributed stiffness calculated at the center of the PE patch ( ). The magnitude of the bending moment FRF in Figure 9 is characterized by peaks at the natural frequencies of the odd modes, whereas the magnitude of the shear force FRF is characterized by peaks at the natural frequencies of the even modes. Indeed, the bending moment distribution corresponds to symmetric loads acting on a portion of the slat, so it only excites the odd modes, whereas the shear force corresponds to an anti-symmetric load and only excites the even modes.

Equivalent Model of a Portion of the Slat
Typically, a PE patch has a longitudinal dimension that is much shorter than the length of the host structure. Moreover, the cross-section dimensions of the harvester are significantly smaller than those of the cross-section of the structure. These large differences in the geometrical dimensions represent an important issue related to the development of the FE model. Indeed, many very thin elements are required to discretize the system in the region close to the PE patch, resulting in a large computational effort. To overcome this drawback, a portion of the vibrating structure containing the PE patch is first considered. Then, the cross-section of this portion is reduced to a simpler and smaller equivalent cross-section to achieve a better aspect ratio between the dimensions of the equivalent cross-section and the ones of the patch. Figure 10a shows that in each portion of the slat, the external loads are the distributed load due to inertia force (Equation (3)), the bending moment , and shear force exerted by the rest of the slat (subscripts R and L refers to the left and right sides, respectively). The bending moment and shear force can be calculated from , and by means of the equilibrium equations. Since the bending moment and shear force can also be exerted by the clamp of a cantilever beam, Figure 10b shows that the slat portion can be converted to an equivalent cantilever beam with the same length and forced by the same external loads acting on that portion of the slat. It is worth noticing that the slat portion and the cantilever beam have the same bending moment distribution. The same bending moment distribution causes the same curvature and the same strain and stress distribution, which guarantees the same performance of the PE patch.

Equivalent Model of a Portion of the Slat
Typically, a PE patch has a longitudinal dimension that is much shorter than the length of the host structure. Moreover, the cross-section dimensions of the harvester are significantly smaller than those of the cross-section of the structure. These large differences in the geometrical dimensions represent an important issue related to the development of the FE model. Indeed, many very thin elements are required to discretize the system in the region close to the PE patch, resulting in a large computational effort. To overcome this drawback, a portion of the vibrating structure containing the PE patch is first considered. Then, the cross-section of this portion is reduced to a simpler and smaller equivalent cross-section to achieve a better aspect ratio between the dimensions of the equivalent cross-section and the ones of the patch. Figure 10a shows that in each portion of the slat, the external loads are the distributed load q due to inertia force (Equation (3)), the bending moment M, and shear force T exerted by the rest of the slat (subscripts R and L refers to the left and right sides, respectively). The bending moment M L and shear force T L can be calculated from M R , T R and q by means of the equilibrium equations. Since the bending moment M L and shear force T L can also be exerted by the clamp of a cantilever beam, Figure 10b shows that the slat portion can be converted to an equivalent cantilever beam with the same length and forced by the same external loads acting on that portion of the slat. It is worth noticing that the slat portion and the cantilever beam have the same bending moment distribution. The same bending moment distribution causes the same curvature and the same strain and stress distribution, which guarantees the same performance of the PE patch. The OCV generated by the PE patch depends on the strain distribution within the patch. Since the PE patch is very thin compared with the slat thickness, it is assumed that the strain distribution within the patch coincides with the strain distribution of the slat just below the patch. In this scenario, the OCV generated by the PE patch bonded to the slat can be reproduced by simply bonding the patch to an underlayer, which reproduces the same strain distribution of the slat just below the patch. By assuming a 1D approximation, the strain distribution ( ) along the longitudinal direction of the slat subject to a bending moment distribution ( ) can be approximated as where / is the module of resistance of the slat cross-section, is the distance from the neutral axis of the slat, is the Young's modulus, and is the moment of inertia of the cross-section of the slat. Equation (26) shows that two structures loaded by the same bending moment distribution with different cross-section and Young's modulus values may provide the same strain distribution if they have the same module of resistance. Therefore, it is possible to assume that the equivalent beam has a rectangular cross-section and that the PE patch is bonded to the upper surface of the beam. Figure 11a shows the actual cross-section of the slat, and Figure 11b schematizes the cross-section of the equivalent beam. The OCV generated by the PE patch depends on the strain distribution within the patch. Since the PE patch is very thin compared with the slat thickness, it is assumed that the strain distribution within the patch coincides with the strain distribution of the slat just below the patch. In this scenario, the OCV generated by the PE patch bonded to the slat can be reproduced by simply bonding the patch to an underlayer, which reproduces the same strain distribution of the slat just below the patch. By assuming a 1D approximation, the strain distribution S(x) along the longitudinal direction of the slat subject to a bending moment distribution M f (x) can be approximated as where W R = EI/z n is the module of resistance of the slat cross-section, z n is the distance from the neutral axis of the slat, E is the Young's modulus, and I is the moment of inertia of the cross-section of the slat. Equation (26) shows that two structures loaded by the same bending moment distribution with different cross-section and Young's modulus values may provide the same strain distribution if they have the same module of resistance. Therefore, it is possible to assume that the equivalent beam has a rectangular cross-section and that the PE patch is bonded to the upper surface of the beam. Figure 11a shows the actual cross-section of the slat, and Figure 11b schematizes the cross-section of the equivalent beam. The OCV generated by the PE patch depends on the strain distribution within the patch. Since the PE patch is very thin compared with the slat thickness, it is assumed that the strain distribution within the patch coincides with the strain distribution of the slat just below the patch. In this scenario, the OCV generated by the PE patch bonded to the slat can be reproduced by simply bonding the patch to an underlayer, which reproduces the same strain distribution of the slat just below the patch. By assuming a 1D approximation, the strain distribution ( ) along the longitudinal direction of the slat subject to a bending moment distribution ( ) can be approximated as where / is the module of resistance of the slat cross-section, is the distance from the neutral axis of the slat, is the Young's modulus, and is the moment of inertia of the cross-section of the slat. Equation (26) shows that two structures loaded by the same bending moment distribution with different cross-section and Young's modulus values may provide the same strain distribution if they have the same module of resistance. Therefore, it is possible to assume that the equivalent beam has a rectangular cross-section and that the PE patch is bonded to the upper surface of the beam. Figure 11a shows the actual cross-section of the slat, and Figure 11b schematizes the cross-section of the equivalent beam. The module of resistance W R,eq of the equivalent beam with rectangular cross-section is W R,eq = E eq I eq z n,eq = E eq B eq H 2 where B eq and H eq are the width and thickness of the equivalent beam, respectively, and z n,eq = H eq /2. Width and thickness B eq and H eq are defined to obtain a suitable aspect ratio between the cross-sections of the beam and the PE patch. Consequently, the module of resistance W R,eq only equals the module of resistance W R of the cross-section of the slat if the Young's modulus of the equivalent beam takes this value: (28)

Finite Element Model
The numerical FE model of the patch can be implemented in COMSOL Multiphysics ® . As previously mentioned, the equivalent beam has a rectangular cross-section, and the PE patch is bonded to the top side of the beam. Table 2 shows that the PE patch has an overall length and width of 0.1 m and 0.018 m, respectively. However, the active portion of the patch, containing the piezoelectric material, has smaller dimensions. Indeed, the active portion is enclosed by a polymeric material (Kapton). The length of the equivalent beam was chosen to be slightly longer than the active portion of the PE patch. The width of the beam is equal to the width of the patch. The determination of the thickness of the equivalent beam must consider that the FE model will be used to perform a frequency domain analysis. The equivalent cantilever beam has its own natural frequency ( f n,eq ) that can be excited by the applied loads but does not correspond to an actual natural frequency of the slat: f n,eq = 1.875 2 2πL 2 eq E eq B eq H 3 eq 12ρ eq A eq . (29) Hence, the thickness H eq was defined to have a natural frequency of the equivalent beam far from the range of frequency of interest (above 2000 Hz). Table 4 shows the geometrical features of the FE model. The cross-section area and moment of inertia of the equivalent beam can be easily calculated using the formulas of the rectangular cross-section. Finally, the equivalent Young's modulus is obtained from Equation (28). The density of the beam is determined by assuming the mass per unit length µ eq of the beam equal to that of the slat µ, yielding Table 5 presents the mechanical parameters of the equivalent beam.

Parameter Unit Value
Young's modulus (E eq ) GPa 5.00 × 10 14 Density (ρ eq ) kg m −3 1.51 × 10 5 Cross-section area (A eq ) m 2 2.80 × 10 −5 Cross-section moment of inertia (I eq ) m 4 9.33 × 10 −12 Natural frequency ( f n,eq ) Hz 2297 Since the geometry of the equivalent beam in Figure 11b is symmetric with respect to the x-z plane, the model only considers one half of the system to reduce the computational effort; see Figure 12a.  Hz 2297 Since the geometry of the equivalent beam in Figure 11b is symmetric with respect to the x-z plane, the model only considers one half of the system to reduce the computational effort; see Figure 12a. The constitutive equations of a PE material allow for electromechanical coupling in the multiphysics problem. The strain-charge and stress-forms are available in COMSOL Multiphysics ® . In this analysis, the strain-charge form was adopted, since it leads to a drastic reduction in the number of material constants when some assumptions are made. The constitutive equations in the strain-charge form are as follows [12]: where are the mechanical strain components (6 × 1), represents the electric displacement components (3 × 1), represents the mechanical stress components (6 × 1), represents the electric field components (3 × 1), represents the compliance constants in a constant electric field (6 × 6), represents the piezoelectric strain constants (3 × 6), and represents the permittivity constants at constant stress (3 × 3). Therefore, the full characterization of a piezoelectric material requires 63 material constants. The following assumptions that are usually made when thin piezoelectric layers are considered [13] can be introduced: 1. Electrodes of the PE material acting along the 3-axis of the local reference frame of the material, i.e., 0, ≠ 0. 2. Thin beam, i.e., The constitutive equations of a PE material allow for electromechanical coupling in the multiphysics problem. The strain-charge and stress-forms are available in COMSOL Multiphysics ® . In this analysis, the strain-charge form was adopted, since it leads to a drastic reduction in the number of material constants when some assumptions are made. The constitutive equations in the strain-charge form are as follows [12]: where S i are the mechanical strain components (6 × 1), D i represents the electric displacement components (3 × 1), T j represents the mechanical stress components (6 × 1), E j represents the electric field components (3 × 1), s E ij represents the compliance constants in a constant electric field (6 × 6), d ij represents the piezoelectric strain constants (3 × 6), and ε T ij represents the permittivity constants at constant stress (3 × 3). Therefore, the full characterization of a piezoelectric material requires 63 material constants. The following assumptions that are usually made when thin piezoelectric layers are considered [13] can be introduced:

1.
Electrodes of the PE material acting along the 3-axis of the local reference frame of the material, i.e., E 1 = E 2 = 0, E 3 = 0.
In this case, the 9 × 9 matrix of material constants in Equation (31) can be reduced to a 2 × 2 symmetric matrix with only three material constants. The reduced constitutive equation becomes [6] S The compliance constant s E 11 and the piezoelectric strain constant d 31 are provided in the datasheet of the considered PE patch [10]. The permittivity constant ε T 33 is not provided in the datasheet; however, it can be calculated as [4,6] ε T 33 = ε 0 where ε 0 is the vacuum permittivity. Tables 3 and 6 represent the electromechanical parameters of the PE patch implemented in COMSOL. Two mechanical boundary conditions are applied to the model (Figure 12a). Firstly, a fixed constraint to the clamped surface of the cantilever beam is imposed. Secondly, a symmetry condition is imposed on the surfaces of the beam and the patch belonging to the plane of symmetry. Finally, both the beam and the piezo are loaded by a distributed load per unit length corresponding to the inertia force obtained from Equations (3) and (4).
Two electrostatic boundary conditions are applied (Figure 12b). A floating potential is imposed on the top side of the PE patch, whereas its bottom side is set to ground potential. In this way, the value of the floating potential, which is computed from the solution of the multiphysics problem, corresponds to the OCV at the patch terminals.
The geometry is discretized by using a mapped mesh (second-order FE elements), which is obtained by extruding a quadrilateral mesh (element size range 0.5-1 mm). To obtain accurate results, four elements are used along the thickness of the beam and six elements are used along the thickness of the patch. Figure 13 shows the mapped mesh used to discretize the geometry.

Orthotropic material.
In this case, the 9 × 9 matrix of material constants in Equation (31) can be reduced to a 2 × 2 symmetric matrix with only three material constants. The reduced constitutive equation becomes [6] The compliance constant and the piezoelectric strain constant are provided in the datasheet of the considered PE patch [10]. The permittivity constant is not provided in the datasheet; however, it can be calculated as [4,6] + where is the vacuum permittivity. Tables 3 and 6 represent the electromechanical parameters of the PE patch implemented in COMSOL. Two mechanical boundary conditions are applied to the model (Figure 12a). Firstly, a fixed constraint to the clamped surface of the cantilever beam is imposed. Secondly, a symmetry condition is imposed on the surfaces of the beam and the patch belonging to the plane of symmetry. Finally, both the beam and the piezo are loaded by a distributed load per unit length corresponding to the inertia force obtained from Equations (3) and (4).
Two electrostatic boundary conditions are applied (Figure 12b). A floating potential is imposed on the top side of the PE patch, whereas its bottom side is set to ground potential. In this way, the value of the floating potential, which is computed from the solution of the multiphysics problem, corresponds to the OCV at the patch terminals.
The geometry is discretized by using a mapped mesh (second-order FE elements), which is obtained by extruding a quadrilateral mesh (element size range 0.5 1 mm). To obtain accurate results, four elements are used along the thickness of the beam and six elements are used along the thickness of the patch. Figure 13 shows the mapped mesh used to discretize the geometry.

Validation of the Integrated Analytical-Numerical Method
If the PE layer is directly bonded to the slat surface, the OCV FRFs can be calculated both by means of the integrated model proposed in this paper and by means of the analytical model presented in [7] modified to consider contact stiffness k x . Figure 14 compares the magnitude of the calculated OCV FRF considering k x = 0 Nm −2 (deployed slat, Figure 14a) and k x = 10 6 Nm −2 (retracted slat, Figure 14b). Figure 13. Mapped mesh used to discretize the geometry.

Validation of the Integrated Analytical-Numerical Method
If the PE layer is directly bonded to the slat surface, the OCV FRFs can be calculated both by means of the integrated model proposed in this paper and by means of the analytical model presented in [7] modified to consider contact stiffness . Figure 14 compares the magnitude of the calculated OCV FRF considering 0 − (deployed slat, Figure 14a) and 10 6 − (retracted slat, Figure 14b).  Figure 14a,b shows that the integrated and analytical models, respectively, were in good agreement, even if different values of were adopted.

Effect of Contact Stiffness on Generated Voltage
The deployment and retraction of the slat modify the contact stiffness between the slat and the wing and influence the generated voltage. To analyze this effect, a parametric analysis was carried out by means of the integrated model. The scheme in Figure 5 shows that the numerical OCV FRF, obtained using COMSOL, could be imported in MATLAB to calculate the using Equations (1) and (2). Figure 15 represents the generated by the PE patch as a function of the distributed stiffness . It can be noted that the effect of on voltage output was only valuable above 10 6 Nm −2 . For large values of the generated voltage drastically decreased because the slat began to behave as a rigid body connected to the wing edge. These results agreed with the trend of the natural frequencies against contact stiffness, which is presented in Figure 8.  Figure 14a,b shows that the integrated and analytical models, respectively, were in good agreement, even if different values of k x were adopted.

Effect of Contact Stiffness on Generated Voltage
The deployment and retraction of the slat modify the contact stiffness between the slat and the wing and influence the generated voltage. To analyze this effect, a parametric analysis was carried out by means of the integrated model. The scheme in Figure 5 shows that the numerical OCV FRF, obtained using COMSOL, could be imported in MATLAB to calculate the V RMS using Equations (1) and (2). Figure 15 represents the V RMS generated by the PE patch as a function of the distributed stiffness k x . It can be noted that the effect of k x on voltage output was only valuable above 10 6 Nm −2 . For large values of k x the generated voltage drastically decreased because the slat began to behave as a rigid body connected to the wing edge. These results agreed with the trend of the natural frequencies against contact stiffness, which is presented in Figure 8.

Effect of Interposed Layers on Performance
The integrated model was found to be suitable for analyzing the effect of the interposed layers between the PE patch and the slat skin. Typical layers are the ones of the adhesive and Kapton, which are used to seal the PE material. The PE layer can be embedded into a more complex electric generator, such as the one presented in [6,14], which was a hybrid Thermo-Piezoelectric Generator (TPEG). Hence, not only adhesive and Kapton layers but also thermoelectric and conductive layers may be interposed between the PE layer and the surface of the structure. Finally, the FE model also allows for the evaluation of the performance of several PE generators piled one onto the other to increase the total generated electric power. Figure 15.
vs. . The values of were derived from the numerical OCV FRFs obt in COMSOL.

Effect of Interposed Layers on Performance
The integrated model was found to be suitable for analyzing the effect of the i posed layers between the PE patch and the slat skin. Typical layers are the ones o adhesive and Kapton, which are used to seal the PE material. The PE layer can be em ded into a more complex electric generator, such as the one presented in [6] and which was a hybrid Thermo-Piezoelectric Generator (TPEG). Hence, not only adh and Kapton layers but also thermoelectric and conductive layers may be interposed tween the PE layer and the surface of the structure. Finally, the FE model also allow the evaluation of the performance of several PE generators piled one onto the oth increase the total generated electric power.
In the framework of this research, a layer of an isotropic polymeric material added to the FE model. Figure 16 shows the cross-section of the model.  In the framework of this research, a layer of an isotropic polymeric material was added to the FE model. Figure 16 shows the cross-section of the model. Figure 15.
vs. . The values of were derived from the numerical OCV FRFs obtained in COMSOL.

Effect of Interposed Layers on Performance
The integrated model was found to be suitable for analyzing the effect of the interposed layers between the PE patch and the slat skin. Typical layers are the ones of the adhesive and Kapton, which are used to seal the PE material. The PE layer can be embedded into a more complex electric generator, such as the one presented in [6] and [14], which was a hybrid Thermo-Piezoelectric Generator (TPEG). Hence, not only adhesive and Kapton layers but also thermoelectric and conductive layers may be interposed between the PE layer and the surface of the structure. Finally, the FE model also allows for the evaluation of the performance of several PE generators piled one onto the other to increase the total generated electric power.
In the framework of this research, a layer of an isotropic polymeric material was added to the FE model. Figure 16 shows the cross-section of the model. The layer can be made of two different materials: Material 1 has the mechanical properties typical of polymers such as epoxy and Kapton; Material 2 has a Young's modulus about one order of magnitude smaller than that of Material 1. The geometrical and The layer can be made of two different materials: Material 1 has the mechanical properties typical of polymers such as epoxy and Kapton; Material 2 has a Young's modulus about one order of magnitude smaller than that of Material 1. The geometrical and mechanical properties of the added layers are shown in Table 7. It is worth noticing that the Young's modulus of the added layers was about one/two orders of magnitude smaller than the ones of the PE and slat materials. A series of parametric simulations was carried out to evaluate the generated voltage as a function of the thickness H a of the polymeric layer. It was assumed that k x = 0 Nm −2 (deployed slat). The model was discretized using the same mesh as the previous simulations. The polymeric layer was discretized using five elements along its thickness. Figure 17a shows the RMS values of the generated OC voltage as a function of the thickness of the interposed layer and the Young's modulus of the polymeric layer. Figure 17b shows the distance of the center of the PE patch from the neutral axis of the equivalent beam. mechanical properties of the added layers are shown in Table 7. It is worth noticing that the Young's modulus of the added layers was about one/two orders of magnitude smaller than the ones of the PE and slat materials. A series of parametric simulations was carried out to evaluate the generated voltage as a function of the thickness of the polymeric layer. It was assumed that 0 − (deployed slat). The model was discretized using the same mesh as the previous simulations. The polymeric layer was discretized using five elements along its thickness. Figure  17a shows the RMS values of the generated OC voltage as a function of the thickness of the interposed layer and the Young's modulus of the polymeric layer. Figure 17b shows the distance of the center of the PE patch from the neutral axis of the equivalent beam.  Figure 17 highlights that moving away the center of the PE patch from the neutral axis increased its performance, since, as shown in Equation (26), the strain within the patch increased. However, this advantageous effect gradually decreased, since the strain was not correctly transmitted from the slat to the patch when the thickness of the interposed polymeric layer was very large. Indeed, for large thicknesses, the distance from the neutral axis linearly increased, whereas the performance showed a decrease in the growth rate.
The comparison between the two curves of Figure 17 shows that if the polymeric layer had a small Young's modulus, the positive effect due to the increased distance from the neutral axis vanished for smaller values of thickness.
It is worth noticing that if the PE patch was located inside the slat structure, the increase in the polymeric layer thickness led to a small decrease in the distance from the neutral axis, with a reduction in the generated voltage.  Figure 17 highlights that moving away the center of the PE patch from the neutral axis increased its performance, since, as shown in Equation (26), the strain within the patch increased. However, this advantageous effect gradually decreased, since the strain was not correctly transmitted from the slat to the patch when the thickness of the interposed polymeric layer was very large. Indeed, for large thicknesses, the distance from the neutral axis linearly increased, whereas the performance showed a decrease in the growth rate.
The comparison between the two curves of Figure 17 shows that if the polymeric layer had a small Young's modulus, the positive effect due to the increased distance from the neutral axis vanished for smaller values of thickness.
It is worth noticing that if the PE patch was located inside the slat structure, the increase in the polymeric layer thickness led to a small decrease in the distance from the neutral axis, with a reduction in the generated voltage.

Conclusions
An integrated analytical-numerical model for the prediction of the voltage generated by a small PE patch mounted on a long slender body that vibrates due to random excitation has been presented and applied to a wing slat.
The introduction of a variable contact stiffness between the slat and the wing edge makes it possible to simulate the behavior of the PE patch, both when the slat is retracted and when the slat is deployed. Numerical results showed that the PE patch generated much less voltage when the slat was retracted, because the large contact stiffness reduced the deformability of the slat.
The transmission of the strain from the slat surface to the active PE layer through the intermediate layers (adhesive and protective) is an important issue. To this end, the integrated model was used to simulate the effect of a passive layer of increasing thickness located between the slat and the active PE layer. The results showed that a small increase in the thickness of the intermediate layer had a beneficial effect. This effect was the result of the increase in the distance between the PE layer and the neutral axis of the cross-section caused by the thicker layer. It vanished when the intermediate layer became very thick (some mm with the considered materials).
A further application of the integrated model will be the simulation of more complex sandwich structures such as generators that include both PE harvesters and thermo-electric harvesters based on the Seebeck effect. Funding: This research was funded by H2020 InComEss project-Innovative polymer-based composite systems for high-efficiency energy scavenging and storage (Project ID: 862597).