Nonlinear Analyses of Porous Functionally Graded Sandwich Piezoelectric Nano-Energy Harvesters under Compressive Axial Loading

: In this study, a sandwich piezoelectric nano-energy harvester model under compressive axial loading with a core layer fabricated of functionally graded (FG) porous material is presented based on the nonlocal strain gradient theory (NSGT). The von Karman type geometric nonlinearity and the axial loading were considered. The electromechanical governing equations were obtained using Hamilton’s principle. The nonlinear vibration frequencies, root mean square (RMS) voltage output and static buckling were obtained using the Galerkin method. The effects of different types of porous distribution, porosity coefﬁcients, length scale parameters, nonlocal parameters, ﬂexoelectricity, excitation frequencies, lumped mass and axial loads on the natural frequency and voltage output of nanobeams were investigated. Results show that the porous distributions, porosity coefﬁcient of porous materials, the excitation frequencies and the axial load have a large effect on the natural frequency and voltage output of the sandwiched piezoelectric nanobeams. When the NSGT is considered, the critical buckling load depends on the values of the nonlocal parameters and strain gradient constants. In addition, the electromechanical conversion efﬁciency of the post-buckling process is signiﬁcantly higher than that of the pre-buckling process. The ﬂexoelectric effect can signiﬁcantly increase the RMS voltage output of the energy harvester.


Introduction
With the continuous development of advanced manufacturing technology, people have made a variety of electronic devices such as telephone watches, wireless headphones, and microelectronic systems. At present, most electronic devices need external batteries to provide power for normal operation. However, the energy supply cycle is limited and the battery needs to be replaced regularly. The extensive use and replacement of batteries have brought a series of problems such as environmental pollution. To find a way to overcome the difficulty, researchers have tried to harvest energy from the surrounding environment to replace the demand for batteries. For example, converting energy such as electromagnetic [1], piezoelectric [2], and electrostatic [3] into electrical energy. Piezoelectric energy harvesters are widely used because of their high energy density, no magnetic interference, and easy processing, miniaturization, and integration in structure. At present, many researchers have studied piezoelectric energy harvesters. For instance, Erturk and Inman [4,5] conducted theoretical and experimental studies on the electromechanical responses of cantilevered unimorph and bimorph piezoelectric energy harvesters based on the Euler-Bernoulli beam theory. An exact analytical solution for a cantilevered piezoelectric energy harvester was proposed, and the modal electromechanical coupling and the dependence of the electrical output on the locations of the electrode were discussed. Wang and Wang [6] developed an analytical model of a nanoscale unimorph piezoelectric energy harvester containing the flexoelectric effect. An analytical model of a piezoelectric cantilevered energy harvester by through-width delamination was proposed by Zeng et al. [7]. Stamatellou et al. [8] experimentally investigated the performance of a film-type piezoelectric energy harvester under simultaneous pneumatic and base excitation. Wang and Wang [9] investigated the energy harvesting performance of a piezoelectric circular nanomembrane under human blood pressure with surface effects.
In recent decades, a lot of research has been performed on piezoelectric energy harvesters. In most studies, the proposed piezoelectric energy harvesters are linear, and work only in a limited excitation frequency spectrum. In order to obtain a wider operating frequency band for the piezoelectric energy harvester, researchers try to consider the nonlinearity in the study of piezoelectric energy harvesters. At present, the research on the nonlinearity of piezoelectric energy harvesters has made great progress. For instance, Costa et al. [10] studied the nonlinear dynamical parameters of the bistable vibration of a piezoelectric energy harvester. Farokhi et al. [11] studied the nonlinear electromechanical behavior of a constrained bimorph piezoelectric energy harvester by using a coupled multimodal fully nonlinear model. It is not difficult to find that considering the nonlinearity cannot only improve the calculation accuracy but also improve the harvesting energy efficiency of energy harvester.
In addition, researchers also found that the broadband behavior of piezoelectric energy harvester can be achieved by applying axial loads to a fixed beam. For example, Masna and Daqaq [12] established an electromechanical model of a clammed single-deformation piezoelectric energy harvester subjected to transverse excitation and static axial load, and made a comprehensive discussion of the energy harvesting of the loaded beam using an axial load. Their results show that the axial load can achieve tunable characteristics and power enhancement. Chen and Yan [13] proposed a nonlinear electromechanical model of an energy harvester based on an axial pre-loaded piezoelectric beam, and combined the flexoelectric effect, which can be operated in pre-buckling (monostable) and post-buckling (bistable) configurations. Varelis and Saravanos [14] proposed a theoretical framework for analyzing the response of piezoelectric actuators and sensor composite laminates pre-and post-buckling. Based on the Euler-Bernoulli beam theory, Zhang et al. [15] investigated the free transverse vibration and buckling characteristics of a double beam system under compressive axial loads. The free buckling and vibration of a set of parallel Timoshenko beams under the axial load was studied by Foroozandeh et al. [16]. In this study, a nonlinear and axial load are considered while analyzing the vibration and bending properties of the sandwich piezoelectric nanobeam.
Piezoelectric materials have received increasing attention because of their excellent electromechanical coupling properties, and are widely used in the field of energy harvesting. When there is mechanical energy in the environment where piezoelectric materials are located, such as wind energy, noise, mechanical vibration, these mechanical energies can cause the deformation of piezoelectric materials. Because of the piezoelectric effect, the voltage potential is generated on the surface of piezoelectric nanomaterials, so the transformation of mechanical energy to electrical energy is realized. It is worth noting that at micro-and nanoscale, the electromechanical properties of piezoelectric nanomaterials and nanostructures are significantly different from those of macroscopic piezoelectric materials. It provides guidance for the development of a piezoelectric micro energy harvester.
Piezoelectric materials show a significant flexoelectric effect when the size is at nanoscale. The flexoelectricity is the electromechanical coupling effect present in a dielectric, which is similar to piezoelectricity. It not only exhibits the polarization response and strain gradient, but also exhibits the electromechanical coupling effect [17]. The flexoelectricity is the basic property of dielectric materials, which not only exists in piezoelectric materials, but also almost exists in all dielectrics. It is noteworthy that the electricity caused by the flexoelectric effect is much larger compared with the piezoelectric effect of nanomaterials. Therefore, the flexoelectric effect is usually considered when studying nanostructures [18][19][20][21][22]. Zeng et al. [18] studied the vibration of piezoelectric sandwich nanobeams based on the theory of nonlocal strain gradients and discussed the flexoelectric effect in detail. Fan et al. [23] modeled the beam as a laminated beam and investigated the effect of flexoelectric layer stiffness on the dynamic response of the structure. Moura et al. [24] investigated the combined piezoelectric and flexoelectric effects in the resonant dynamics of nanocantilever beams. Since the effect of pliability generated in piezoelectric nanostructures cannot be ignored, this study takes into account the flexoelectric effect while analyzing the vibration of functionally graded sandwich piezoelectric nanobeams.
The sandwich structure has excellent properties such as high stiffness and light weight, so it is widely used in engineering fields, such as biomedicine, architecture, automotive and aerospace. The sandwich structure is a composite structure formed by connecting two thin surface layers and one core layer. The porous material is a kind of lightweight material that widely exists in nature, such as metal foams. The porous material has the advantages of high strength and low density, and is often used as energy absorption, noise reduction, flame retardant, heat insulation [25][26][27][28], and is therefore widely used in the automotive industry, aerospace, and other fields. The composite structure with porous materials as the core layer of the sandwich structure is called the sandwich porous structure. Generally, the sandwich porous structure cannot only ensure enough stiffness, but also reduce the quality of the whole structure. In recent years, various theoretical and experimental studies have been conducted on the static and dynamic properties of sandwich porous beams and plates to meet various needs [29][30][31][32]. The functionally graded (FG) porous materials are a hot research topic for new advanced composites [33][34][35]. The internal pores of FG porous materials show gradient distribution in a certain direction. FG porous materials can be varied by local density to obtain the desired natural frequency range of the structure [25]. Previous researchers have focused on studying the dynamic analysis, static bending, and buckling of FG porous plate and beam structures. For example, Zine et al. [36] investigated the bending response of FG porous plates based on the cubic shear deformation theory. Based on two-variable shear deformation and modified coupled stress theory, Sadoughifar et al. [37] investigated the nonlinear bending of FG porous thick circular/annular nano-plates using a generalized differential quadrature method (GDQM). Based on the Timoshenko beam theory, Chen et al. [38] analyzed the static buckling and elastic buckling of shear deformable FG porous beams. In the above study, it can be found that porous beams or porous plates exhibit better mechanical properties than ordinary beams or plates of the same mass.
The structure size [39] exhibits non-negligible scale effects when it is at the micro-and nanoscale. Although researchers have performed a great deal of work using classical continuum mechanics, their application on the nanoscale is questionable. Because in classical continuum mechanics, small-scale effects such as surface effects, lattice spacing, van der Waals forces, chemical bonding, and electric forces are neglected. However, both experimental and atomic simulation results show that at the nanoscale, these small-scale effects cannot be ignored [40]. This is because at small sizes, the lattice spacing between atoms becomes important and the interior of the material is no longer homogenized into a continuum [41]. Therefore, for the scale effects of micro-and nanoscale structures, non-classical continuum theories have been proposed by previous authors to evaluate the mechanical properties of micro-and nanostructures, such as surface elasticity theory [42][43][44], strain gradient elasticity theory [45][46][47][48], and nonlocal elasticity theory [49,50]. In the last decade, numerous researchers [51][52][53] have studied a large number of size dependencies on microand nanostructures based on the nonlocal elasticity theory and strain-gradient elasticity theory. These studies have shown that nonlocal theory can only capture the stiffnesssoftening effect and the strain gradient theory can only capture the stiffness-hardening effect. Recently, researchers have combined the nonlocal elasticity theory and the strain gradient elasticity theory together to propose a new theory, namely, the nonlocal strain gradient theory (NSGT) [54]. The NSGT combines the advantages of both theories so that not only the strain-gradient stress field but also the non-local stress field can be investigated. In addition, the NSGT is able to describe both the stiffness-hardening and stiffness-softening effects. Based on the NSGT, a large number of researchers have investigated the size-dependent bending, buckling, wave propagation, and vibration of FG microand nanoscale structures [18,[55][56][57]. For example, Zeng et al. [18] investigated the free vibration of piezoelectric sandwich nanobeams with flexoelectricity based on the NSGT. Sahmani et al. [57] analyzed the size-dependent nonlinear axial postbuckling properties of functionally graded porous micro/nanoplates enhanced by graphene platelets based on the NSGT. Sharifi et al. [58] analyzed the functional gradient piezoelectric nanoplate vibration using the NSGT. Yademellat et al. [19] performed a nonlinear electromechanical analysis of micro/nanobeams based on the NSGT taking into account the flexoelectric and piezoelectric effects.
Generally, the objects attached to the beam and plate belong to the subsystem of the whole system, and the subsystem is usually simplified as a lumped (or concentrated) mass in the calculation. In engineering practice, researchers often encounter vibration problems occurring in beam-like structures with lumped masses attached to them, such as automobiles, spacecraft, etc., [59][60][61]. From previous studies, it is found that the rotational inertia of the lumped mass has a significant effect on the dynamic characteristics of mass loaded beams, especially for higher order vibration modes [62][63][64]. In addition, ignoring the effect of concentrated mass may overestimate the natural frequency and may lead to incorrect system response. Therefore, it is important to consider the effect of the lumped mass. Considering the lumped mass not only improves the calculation accuracy of the system frequency, but also increases the efficiency of the captured energy of the piezoelectric energy harvester [65]. For example, Wang [59] performed a free vibration analysis on a lumped mass uniform beam containing translational and rotational inertia. Sun et al. [66] studied the mechanism of the effect of lumped mass on the flutter characteristics of structures, and the optimal ranges for the weight and location of the lumped mass are provided.
Therefore, in this work, we propose a piezoelectric and flexoelectric energy harvester based on the NSGT for FG porous material sandwich nanobeams under axial compressive loading. According to the amplitude of axial load, the energy harvester can work in prebuckling state or postbuckling state. The work is organized as follows. In Section 2, the basic equations are derived. An energy harvester model sandwich piezoelectric beam composed of FG porous material was developed, in which the FG layer is considered with three porosities. The von Karman type nonlinear vibration behavior was considered. Then the nonlinear electromechanical coupling control equation was obtained using the Hamiltonian principle. In Section 3, the static buckling, natural frequency, and output voltage are solved using the Galerkin method. In Section 4, the model is validated and the effects of different types of distribution, porosity coefficients, length scale parameters, nonlocal parameters, excitation frequency, and axial loads on the vibration frequency and output voltage of the nanobeam are discussed. In Section 5, the main results of this paper are summarized, and some conclusions are drawn.

Mathematical Formulation
Shown in Figure 1 is the sketch of a piezoelectric energy harvester with a simply supported nanobeam settled on a base. The beam is composed of three layers, i.e., a core layer manufactured of functionally graded porous materials, and two piezoelectric layers attached perfectly on the core layer and no slip between layers. The adhesive thickness of the bonding interface is neglected. The elastic constant of the piezoelectric layer and core layer are denoted as C 11p and C 11c (z), respectively. The mass density of the piezoelectric layer and core layer are denoted as ρ p and ρ c , respectively. The height of two piezoelectric layers is h p , and of the core layer is h c . The length and the width of three layers are the same, which is denoted as L, b, respectively. The height of the entire beam is h t = 2h p + h c . There is a lumped mass m s attached on the midspan of the nanobeam to improve the energy capture efficiency of the energy harvester. A compressive load is applied axially on the nanobeam. The base is applied a displacement excitation w b (t) = W 0 cos(ωt). The two piezoelectric layers are connected in series in circuit, and connected to an external resistor R.
The two piezoelectric layers are connected in series in circuit, and connected to an external resistor R . The core layer is manufactured of functionally graded porous materials. In this study, three typical porosity distributions are considered [38,67,68]. The elastic properties and   for the open-cell model [68]. A diagrammatic sketch of those three distributions is shown in Figure 2. For the convenience of comparing discussion, the total mass of the entire beam with those three types of distributions is assumed to be equal, which provides The core layer is manufactured of functionally graded porous materials. In this study, three typical porosity distributions are considered [38,67,68]. The elastic properties and mass density are described by Equation (1) for distribution 1, Equation (2) for distribution 2, and Equation (3) for distribution 3.
where C max (ρ max ) and C min (ρ min ) are the maximum value and minimum value of elastic constant (mass density) for porosity distributions, respectively. χ 0 and χ m represent the porosity coefficients which are expressed as χ 0 = 1 − C min /C max and χ m = 1 − ρ min /ρ max , respectively. The relationship between χ 0 and χ m is expressed as for the open-cell model [68]. A diagrammatic sketch of those three distributions is shown in Figure 2. For the convenience of comparing discussion, the total mass of the entire beam with those three types of distributions is assumed to be equal, which provides  The approximate process of the governing equations and boundary conditions of the model obtained is shown in Figure 3. Based on the Euler-Bernoulli beam assumption, the displacement field of the beam are: The approximate process of the governing equations and boundary conditions of the model obtained is shown in Figure 3. The approximate process of the governing equations and boundary conditions of the model obtained is shown in Figure 3. Based on the Euler-Bernoulli beam assumption, the displacement field of the beam are: Based on the Euler-Bernoulli beam assumption, the displacement field of the beam are: Considering the von Karman geometric nonlinearity, the strain fields can be expressed as: It is assumed that the electric field only exists in the thickness direction which is expressed as E Z = −V/h p for both piezoelectric layers. Considering the NSGT with flexoelectric effect, the constitutive equations for the nanobeam are written as [18]: where e 0 a and l are the nonlocal parameter and the strain gradient constant, respectively. e 311 , µ 3113 and a 33 are the piezoelectric coefficient, the electric field-strain gradient coupling coefficient and the dielectric constant, respectively. D z is the electric displacement. The variation in the electric entropy can be calculated as [18]: By substituting Equation (6) into Equation (10), and performing the calculation, one obtains: where By combining Equation (12) with Equations (7)-(9), one has: where Neglecting the rotational inertia of the beam, the total kinetic energy of the beam is written as: where The variation in the total work performed by the compressive axial load P, damping force and the electric charge output Q(t) is expressed as: where c v is the viscous damping coefficient. The Hamilton principle is utilized as: By substituting Equations (11)-(15) into Equation (16), and performing the calculation, one obtains: where (x − x 0 ) is defined as: The two piezoelectric layers are connected in series in the circuit. Therefore, the current in the circuit is calculated as Q(t). Then, the nonlinear electromechanical coupling equations can be obtained as: By substituting Equation (13), into Equations (19)-(21), the nonlinear electromechanical coupling equations are rewritten as: For the simply supported boundary, the lateral displacements and the moments at both ends of the beam are 0, i.e., In the axial direction, the compressive load P equals to the inner force N xx at both ends, i.e.: Compared with the lateral vibration, the axial vibration of a slender beam is less important witch can be neglected, i.e., .. u = 0. As a result, from Equation (22) Combining with the boundary condition Equation (26) with Equation (27), one obtains: ∂u ∂x With consideration of Equation (28), the nonlinear electromechanical coupling equations can be deduced as:

Static Buckling Analyses
The static buckling solution can be obtained by dropping the damping terms and time-related terms. Equation (29) is reduced as: The lateral displacement is assumed in the form Equation (31) can be written as 4 , ∂w ∂x equals the left hand side of Equation (31). Then the Galerkin method is applied as: By substituting Equation (32) into Equation (33) and applying the integration, one obtains: In the prebuckling process the lateral displacement is zero. Solving the expression of P in Equation (34) with W s = 0, the Eigen critical buckling force is obtained as: In Equation (35), B 11 is related to the porosity distribution, and whose expression is defined before, which can be solved as: In the postbuckling process, the maximum vertical displacement can be obtained by solving the expression of W s in Equation (34), i.e.: where B 11 is shown in Equation (36), and the other porosity related coefficient A 11 is calculated as:

Natural Frequency Analysis
Firstly, the free vibration analysis is performed to obtain the natural frequency of the structure. Setting the base excitation w b , damping coefficient c v and voltage potential V to be zero, neglecting the nonlinear terms, and assuming harmonic excitation, the lateral displacement is assumed in the form at the prebuckling state: Appl. Sci. 2021, 11, 11787 11 of 27 By substituting Equation (39) into Equation (29), applying the Galerkin method, one obtains: −L 2 π 2 P L 2 + (e 0 a) 2 π 2 We jω n t + π 4 bB 11 L 2 + l 2 π 2 We jω n t Equation (40) can be written as: ..
where ω pre is so-called linear natural frequency at prebuckling state and: It is seen in Equation (42) that the natural frequency is related to the compressive load P.
When the beam is buckled, the lateral displacement is in the from: where W s sin( πx L ) is the static vertical displacement, and the expression of W s is provided in Equation (37) By substituting Equation (43) into Equation (29), neglecting the nonlinear terms, and applying the Galerkin method, one obtains: Equation (44) can be written as: ..
where ω post is the natural frequency at postbuckling state and: In Equation (46) the porosity related coefficient B 11 is shown in Equation (36), and the other porosity related coefficient I 0 is calculated as:

Electromechanical Responses
When the base excitation is applied (w b = W b cos(ωt)), the vibration response is assumed to be of the following form: By substituting Equation (48) in Equations (29) and (30), and applying the Galerkin method, one obtains: where ω 2 n is defined in Equation (42), and: The harmonic balance method is used to solve Equations (49) and (50). Setting: where W s is the static displacement by buckling. When the beam is in the state of prebuckling (P<P cr ), W s is zero, and when the beam is in the state of postbuckling (P≥P cr ), W s is calculated by Equation (37). Substituting Equations (52) and (53) into Equations (49) and (50), the coefficients for cos ωt and sin ωt should be balance. Then the differential equations are converted to a series of nonlinear algebraic equations such as: in which, a 1 , b 1 , c 1 , d 1 , c 2 , and d 2 can be solved numerically. The generated voltage is the alternating voltage, which is described as a function of time. In order to evaluate the energy capture efficiency conveniently, the root mean square (RMS) voltage V rms is introduced as follows: where T V is the period of the alternating voltage. The RMS voltage is the square root of the mean over one cycle of the square of the instantaneous voltage. In fact, the RMS voltage is equal to the value of the direct current that would produce the same average power dissipation in a resistive load [69].

Results and Discussion
Firstly, the validation of the present model is validated. In Table 1, the natural frequencies ( ω n L 2 /h t 12ρ max /C max ) with consideration of the nonlocal strain gradient are compared with those in Ref. [70]. For comparison purposes, the electric constants and mass densities of the core layer and piezoelectric layers are assumed to be the same, and the piezoelectricity and flexoelectricity are neglected. The lumped mass is set as zero, and the axially load is released. The length and the thickness of the beam are assumed to be L = 50 nm and h t = 1 nm, respectively. Results show good agreements with Ref. [70]. In addition, in Table 1, it can be seen that the results obtained in the present study are slightly higher than those in Ref. [70] by 0.019%. The reason why the results in the present study are slightly higher than those in Ref. [70], is that the rotational inertia is neglected in the present study, which is considered in Ref. [70]. Table 1. Comparison of the natural frequencies ((ω n L 2 /h t ) 12ρ max /C max ) with available date with consideration of nonlocal strain gradient (L = 50 nm, h t = 1 nm).  In Figure 4, the dimensionless natural frequency ((ω n L 2 /h t ) 12ρ max /C max ) with the dimensionless axial load (PL 2 /(π 2 bB 11 )) are compared with those in Ref. [13]. The material properties are E = 126GPa, e 311 = −6.5C/m 2 , µ 3113 = −1.0 × 10 −7 C/m, a 33 = 1.8 × 10 −8 C/(V · m). Results show good agreements with Ref. [13]. It seems that the natural frequency decreases with the increase in the compressive load P at the prebuckling state and reaches zero when P = P cr , and increases with the increase in the compressive load P at the postbuckling state. A similar phenomenon is also found in the experimental work [71].
In the text bellow, for the convenience of discussion, some values of material constants and geometric dimensions are provided (unless otherwise specified). The material constants of the piezoelectric material are c 11p = 126GPa, ρ p = 7800kg/m 3 , e 311 = 6.5C/m 2 , µ 3113 = 0.1µC/m, and a 33 = 1.3 × 10 −8 C/(V · m). The maximum elastic constant of the core layer is C max = 200GPa, the maximum mass density of the core layer is ρ max = 8700 kg/m 3 . The geometric dimension of the beam is fixed to be L : b : h t : h p = 100:10:1:0.1. The nonlocal parameter is set as e 0 a = 0.1L, the strain gradient constant is l = 0.2L.
In Figure 5, effects of the lumped mass on the natural frequency ((ω n L 2 /h t ) 12ρ max /C max ) for distribution 1 are displayed. The porosity is κ 0 = 0 (uniform distribution without porosity), and the nonlocal parameter and strain gradient constant are e 0 a = 0.1L, and l = 0.2L. It seems that at both the prebuckling state and postbuckling state, the lumped mass decreases the natural frequency. In this case, when m s = 0.1I 0 L, m s = 0.5I 0 L, and m s = I 0 L, the natural frequencies are reduced to 94.7%, 79.6% and 68.1% of those without lumped mass, respectively. 33 1.8 10 a =  / ( ) C V m  . Results show good agreements with Ref. [13]. It seems that natural frequency decreases with the increase in the compressive load P at the prebu ling state and reaches zero when cr P P = , and increases with the increase in the comp sive load P at the postbuckling state. A similar phenomenon is also found in the exp imental work [71].   In Figure 6, effects of the nonlocal parameter and strain gradient constant on the nat ural frequency ( ( ) ) and buckling load for distribution 1 are dis played. The porosity is 0 κ (uniform distribution without porosity). It seems that at th prebuckling state, the resonance frequency decreases because the natural frequency de π PL bB Figure 5. Effect of the lumped mass on the natural frequency ( ω n L 2 /h t 12ρ max /C max ) for distribution 1 (κ 0 = 0, e 0 a = 0.1h t , l = 0.2h t , m s = 0). In Figure 6, effects of the nonlocal parameter and strain gradient constant on the natural frequency ( ω n L 2 /h t 12ρ max /C max ) and buckling load for distribution 1 are displayed. The porosity is κ 0 (uniform distribution without porosity). It seems that at the prebuckling state, the resonance frequency decreases because the natural frequency decreases as the axial load increases; this is opposite at the postbuckling state. When the nonlocal strain gradient is considered, the critical buckling load is dependent on the values of the nonlocal parameter e 0 a and strain gradient constant l. The nonlocal parameter e 0 a diminishes the critical buckling load, while the strain gradient constant l enlarges the critical buckling load. It is because the nonlocal strain gradient theory can explain both the strengthening and weakening of the size-effect of the structure. Whether the overall effect of the nonlocal strain gradient strengthening or weakening is related to the relative quantity between the nonlocal parameter e 0 a and the strain gradient constant l. On the other hand, the nonlocal strain gradient shows the opposite effect on the natural frequency at the prebuckling and postbuckling state. The nonlocal parameter e 0 a decreases the natural frequency at the prebuckling state, while it increases the natural frequency at the postbuckling state. The strain gradient constant l increases the natural frequency at the prebuckling state, while it decreases the natural frequency at the postbuckling state. Shown in Figure 7 55.0% and 91.1%, respectively, compar with that of without axial load. The resonance amplitude also decreases as the axial loa increases. When the axial load equals the buckling load cr P , the resonance phenomeno cannot be observed because the natural frequency of the structure decreases to zero as t structure is critically buckled. In such the case, the dimensionless amplitude increas from 0 to 1 as the excitation frequency increases. π PL bB Figure 6. Effects of the nonlocal parameter and strain gradient constant on the natural frequency ( ω n L 2 /h t 12ρ max /C max ) and buckling load for distribution 1 (κ 0 = 0, m s = 0). Figure 7 are the dimensionless amplitude (W max /W 0 )-excitation frequency ( ω n L 2 /h t 12ρ max /C max ) curves for different axial loads at the prebuckling state (P ≤ P cr , m s = 0, R = 1kΩ, W 0 = 0.01h t , κ 0 = 0 for distribution 1). It seems that at the prebuckling state, the dimensionless amplitude increases to a peak value and then decreases with the increase in dimensionless excitation frequency. Moreover, when P = 0.2P cr , P = 0.4P cr , P = 0.6P cr , P = 0.8P cr and P = P cr , the maximum dimensionless amplitude decreases by 10.4%, 22.4%, 36.9%, 55.0% and 91.1%, respectively, compared with that of without axial load. The resonance amplitude also decreases as the axial load increases. When the axial load equals the buckling load P cr , the resonance phenomenon cannot be observed because the natural frequency of the structure decreases to zero as the structure is critically buckled. In such the case, the dimensionless amplitude increases from 0 to 1 as the excitation frequency increases.

Shown in
, 0 0 κ = for distribution 1). It seems that the root mean square voltage output increases to a peak and then decreases as the excitation frequency increases in the prebuckling procedure. The peak value is acquired when the excitation frequency is close to the natural frequency of the structure. Further, as the axial load increases, the peak value decreases. When the axial load is equal to the buckling load, the RMS voltage output increases as the excitation frequency increases, and the value of the RMS voltage remains at a low level.  Figure 8 are the dimensionless root mean square voltage ( a 33 e 311 V rms W max )-excitation frequency (ω n L 2 /h t 12ρ max /C max ) curves for different axial loads at the prebuckling state (P≤ P cr , m s = 0, R = 1kΩ, W 0 = 0.01h t , κ 0 = 0 for distribution 1). It seems that the root mean square voltage output increases to a peak and then decreases as the excitation frequency increases in the prebuckling procedure. The peak value is acquired when the excitation frequency is close to the natural frequency of the structure. Further, as the axial load increases, the peak value decreases. When the axial load is equal to the buckling load, the RMS voltage output increases as the excitation frequency increases, and the value of the RMS voltage remains at a low level. Shown in Figure 9 are the dimensionless vibrational amplitude ((W − W s ) max /W 0 )excitation frequency ((ω n L 2 /h t ) 12ρ max /C max ) curves for different axial loads at the postbuckling state (P ≥ P cr , m s = 0, R = 1kΩ, W 0 = 0.01h t , κ 0 = 0 for distribution 1). It is worth noting that the vibrational amplitude only calculates the vibrational component of the beam, and displacement due to the axial load (W s ) is taken off in Figure 9, because it would not produce an electric field from this part of the displacement. In contrast with the curves at the prebuckling state (shown in Figure 7), the resonance frequency increases as the axial load increases due to the decreases in the natural frequency as the axial load increases at postbuckling state. Further, the resonance amplitude increases as the axial load increases, which is different from the results at the prebuckling state. Furthermore, the value of the dimensionless resonance vibration amplitude at postbuckling state is larger than those at prebuckling state when the same excitation is applied. excitation frequency increases in the prebuckling procedure. The peak value is acquire when the excitation frequency is close to the natural frequency of the structure. Furthe as the axial load increases, the peak value decreases. When the axial load is equal to th buckling load, the RMS voltage output increases as the excitation frequency increases, an the value of the RMS voltage remains at a low level.  . 2021, 11, x FOR PEER REVIEW 18 Shown in Figure 9 are the dimensionless vibrational amplitude (

Shown in
, 0 0 κ = for distribution 1). worth noting that the vibrational amplitude only calculates the vibrational componen the beam, and displacement due to the axial load ( s W ) is taken off in Figure 9, becau would not produce an electric field from this part of the displacement. In contrast w the curves at the prebuckling state (shown in Figure 7), the resonance frequency incre as the axial load increases due to the decreases in the natural frequency as the axial l increases at postbuckling state. Further, the resonance amplitude increases as the a load increases, which is different from the results at the prebuckling state. Furtherm the value of the dimensionless resonance vibration amplitude at postbuckling sta larger than those at prebuckling state when the same excitation is applied.
, 0 0 κ = for distribution 1). As excitation frequency is close to the natural frequency of the structure, the root mean squ  Figure 10 are the dimensionless root mean square voltage ( a 33 e 311 V rms W max )-excitation frequency ((ω n L 2 /h t ) 12ρ max /C max ) curves for different axial loads at the postbuckling state (P ≥ P cr , m s = 0, R = 1kΩ, W 0 = 0.01h t , κ 0 = 0 for distribution 1). As the excitation frequency is close to the natural frequency of the structure, the root mean square voltage reaches a peak. As the axial load increases, the peak value increases, which is different from the postbuckling procedure shown in Figure 8. Further, the peak value of the root mean square voltage output for postbuckling procedure is remarkably larger than that for the prebuckling procedure. It provides a new idea to improve the mechanical electrical conversion efficiency of the structure. The mechanical electrical conversion efficiency in the postbuckling procedure is higher than that in the prebuckling procedure. Shown in Figure 11 are the dimensionless resonance amplitude ( max

Shown in
) curves for different distributions of the core layer ( 0. The resonance amplitude is acquired when the e tion frequency matches the natural frequency of the structure. It seems that when the load is zero, the resonance amplitude shows no difference for three types of distrib In the prebuckling procedure, as the axial load increases, the resonance amplitud creases, and the difference between the three types of distribution increases. Distrib 2 has the largest resonance amplitude, and distribution 1 has the smallest resonanc plitude among these three types of distribution. In the postbuckling procedure, as t ial load increases, the resonance amplitude increases. Distribution 1 has the largest nance amplitude, and distribution 2 has the smallest resonance amplitude among three types of distribution. 16 0 Figure 10. Dimensionless root mean square voltage output-excitation frequency curves for different axial load at postbuckling state (P ≥ P cr , m s = 0, R = 1kΩ, W 0 = 0.01h t , κ 0 = 0 for distribution 1). Figure 11 are the dimensionless resonance amplitude (W max /W 0 )-axial load ( 12L 2 P π 2 bC max h 3 t ) curves for different distributions of the core layer (m s = 0.5I 0 L, W 0 = 0.01h t , κ 0 = 0.5, R = 1kΩ). The resonance amplitude is acquired when the excitation frequency matches the natural frequency of the structure. It seems that when the axial load is zero, the resonance amplitude shows no difference for three types of distribution. In the prebuckling procedure, as the axial load increases, the resonance amplitude decreases, and the difference between the three types of distribution increases. Distribution 2 has the largest resonance amplitude, and distribution 1 has the smallest resonance amplitude among these three types of distribution. In the postbuckling procedure, as the axial load increases, the resonance amplitude increases. Distribution 1 has the largest resonance amplitude, and distribution 2 has the smallest resonance amplitude among these three types of distribution.

Shown in
2 has the largest resonance amplitude, and distribution 1 has the smallest resonan plitude among these three types of distribution. In the postbuckling procedure, as ial load increases, the resonance amplitude increases. Distribution 1 has the large nance amplitude, and distribution 2 has the smallest resonance amplitude amon three types of distribution.
Shown in Figure 12 are the dimensionless optimal root mean square voltage output ( a 33 e 311 V rms W max )-axial load ( 12L 2 P π 2 bC max h 3 t ) curves for different distributions of the core layer (P ≤ P cr , m s = 0.5I 0 L, W 0 = 0.01h t , κ 0 = 0.5, R = 1kΩ). In the prebuckling procedure, it can be seen that the optimal RMS voltage decreases as the axial load increases, and there is little difference between the optimal RMS voltage of distribution 1 and distribution 3, and distribution 2 has the largest optimal RMS voltage. In the postbuckling procedure, as the axial load increases, the optimal RMS voltage increases. Distribution 1 has the largest optimal RMS voltage, and distribution 2 has the smallest optimal RMS voltage among these three type distributions.
Shown in Figure 13 are the dimensionless natural frequency ((ω n L 2 /h t ) 12ρ max /C max )porosity coefficient (κ 0 ) curves for different distributions of the core layer (m s = 0.5I 0 L, W 0 = 0.01h t , R = 1kΩ, P = 0). The influence of the porosity coefficient of the sandwich layer on the vibration frequency of the composite beam is also significant. In the figure, for the same value of the porosity coefficient (κ 0 ), the mass of the three types of distribution is equal. It can be found that the vibration frequencies of distribution 1 and distribution is little difference between the optimal RMS voltage of distribution 1 and distribution and distribution 2 has the largest optimal RMS voltage. In the postbuckling procedure, the axial load increases, the optimal RMS voltage increases. Distribution 1 has the large optimal RMS voltage, and distribution 2 has the smallest optimal RMS voltage amon these three type distributions.
The influence of the porosi coefficient of the sandwich layer on the vibration frequency of the composite beam is al significant. In the figure, for the same value of the porosity coefficient ( 0 κ ), the mass the three types of distribution is equal. It can be found that the vibration frequencies distribution 1 and distribution 3 decrease rapidly with the increase in the porosity coef cient, while the decrease in the vibration frequency of distribution 2 is not obvious. addition, for the same porosity coefficient, the vibration frequency of distribution 2 is th largest, the vibration frequency of distribution 3 is the second largest, and the vibratio frequency of distribution 1 is the smallest. The dimensionless optimal root mean square voltage output ( a 33 e 311 V rms W max ) -porosity coefficients (κ 0 ) curves for different distributions of the core layer are shown in Figure 14 (m s = 0.5I 0 L, W 0 = 0.01h t , R = 1kΩ, P = 0). The influence of the porosity coefficient of the sandwich layer on the root mean square voltage output of the composite beam is also significant. In the figure, the mass of three types of distribution is assumed to be equal. It can be seen from the figure that the optimal RMS voltage output for the three type distributions decreases rapidly with the increasing porosity coefficient. Besides, for the same porosity coefficient, distribution 2 has the largest optimal RMS voltage output, distribution 3 has the second largest optimal RMS voltage output, and distribution 1 has the smallest optimal RMS voltage output among the three types of distribution.
The influence of the porosity c the sandwich layer on the root mean square voltage output of the composite significant. In the figure, the mass of three types of distribution is assumed to can be seen from the figure that the optimal RMS voltage output for the three butions decreases rapidly with the increasing porosity coefficient. Besides, porosity coefficient, distribution 2 has the largest optimal RMS voltage outp tion 3 has the second largest optimal RMS voltage output, and distribution 1 h est optimal RMS voltage output among the three types of distribution.  Figure 15 are the dimensionless amplitude (W max /W 0 ) -excitation amplitude (W 0 ) curves for linear/nonlinear vibration and different porosity coefficients (κ 0 ) (m s = 0, R = 1kΩ, P = 0 for distribution 1). It can be found that the dimensionless amplitude of nonlinear vibration decreases as the excitation amplitude increases, while the dimensionless amplitude of the linear vibration does not change with the excitation amplitude. This is because the stiffness of the structure increases when considering nonlinear vibrations. In addition, for the same excitation amplitude, the larger the porosity coefficient, the larger the dimensionless amplitude.

Shown in
Shown in Figure 16 are the dimensionless RMS voltage output ( a 33 e 311 V rms W max )excitation amplitude (W 0 ) curves for linear/non-linear vibration and different porosity coefficients (κ 0 ) (m s = 0, R = 1kΩ, P = 0 for distribution 1). It can be seen that the RMS voltage output of nonlinear vibration increases with the increase in the excitation amplitude, while the RMS voltage output of linear vibration does not vary with the excitation amplitude. This indicates that considering nonlinear vibration can enhance the mechanical electrical conversion efficiency.
plitude of nonlinear vibration decreases as the excitation amplitude increases dimensionless amplitude of the linear vibration does not change with the exci plitude. This is because the stiffness of the structure increases when considering vibrations. In addition, for the same excitation amplitude, the larger the poro cient, the larger the dimensionless amplitude.  It is assumed that the ratio of the thickness of the piezoelectric layer to ness of the beam / p t h h is assumed to remain constant. Therefore, as the of the beam increases, so does the thickness the piezoelectric layer, the RM It is assumed that the ratio of the thickness of the piezoelectric layer to the total thickness of the beam h p /h t is assumed to remain constant. Therefore, as the total thickness of the beam increases, so does the thickness the piezoelectric layer, the RMS voltage output is proportional to the total thickness of the beam. In order to study the sizedependency of the flexoelectric effect on the performance of energy harvesters. The dimensionless RMS voltage output is processed, so that when the PE is considered only, the dimensionless RMS voltage output is kept constant with the variation on the total thickness of the beam. Figure 17 shows the effect of the total beam thickness on the dimensionless root mean square voltage output with or without considering the flexoelectric effect (FE) (m s = 0, P = 0, R = 1kΩ, κ 0 = 0.5, W 0 = 0.01h t for distribution 1). It can be seen that the dimensionless RMS voltage output is independent of the total beam thickness when only piezoelectric effect (PE) is considered. While the dimensionless RMS voltage output decreases with the increase in the total beam thickness when the piezoelectric and flexoelectric effects are considered simultaneously, which indicates that the flexoelectricity has a significant effect on the performance of the small-scale energy harvester. When the structure size is smaller, the flexoelectric effect is more pronounced.  for distribution 1). It can RMS voltage output increases with the increase in the thickness ratio o when the thickness ratio increases to a certain value, the RMS voltage out little. This is because the total beam thickness is certain, and the thickne lectric layer no longer changes when the thickness of the piezoelectric la value comparable with the total beam thickness. In addition, under the the RMS voltage output considering both piezoelectricity and flexoelec than the RMS voltage output considering only piezoelectricity.

5.0
10 -8 Figure 17. Effect of the total beam thickness on the dimensionless RMS voltage output with or without considering the flexoelectric effect (m s = 0, P = 0, R = 1kΩ, κ 0 = 0.5, W 0 = 0.01h t for distribution 1). Figure 18 shows the effect of the ratio of the piezoelectric layer thickness to the total beam thickness on the dimensionless root mean square voltage output (h t = 3 × 10 −4 m, m s = 0, P = 0, R = 1kΩ, κ 0 = 0.5, W 0 = 0.01h t for distribution 1). It can be seen that the RMS voltage output increases with the increase in the thickness ratio of the beam, and when the thickness ratio increases to a certain value, the RMS voltage output changes very little. This is because the total beam thickness is certain, and the thickness of the piezoelectric layer no longer changes when the thickness of the piezoelectric layer increases to a value comparable with the total beam thickness. In addition, under the same conditions, the RMS voltage output considering both piezoelectricity and flexoelectricity is greater than the RMS voltage output considering only piezoelectricity.
when the thickness ratio increases to a certain value, the RMS voltage output changes ver little. This is because the total beam thickness is certain, and the thickness of the piezoe lectric layer no longer changes when the thickness of the piezoelectric layer increases to value comparable with the total beam thickness. In addition, under the same condition the RMS voltage output considering both piezoelectricity and flexoelectricity is greate than the RMS voltage output considering only piezoelectricity.

Conclusions
In this study, the nonlinear electromechanical behavior of the piezoelectric energy harvester model of FG porous sandwich nanobeams under compressive axial loading is investigated based on the NSGT. The porosity distribution of three different types of FG layers for the Euler-Bernoulli beam theory is discussed. The effects of different types of porosity distributions, porosity coefficients, length scale parameters, nonlocal parameters, flexoelectricity, excitation frequencies, lumped mass and axial loads on the vibration frequency, and voltage output of nanobeams are investigated. The results show that: The porous distribution, the porosity coefficient of porous materials, the excitation frequencies and the axial load have a large effect on the vibration frequency and voltage output of the sandwiched piezoelectric nanobeams. The axial load and excitation frequency also have an effect on the vibration amplitude of the beam. In the critical buckling state, the vibration amplitude of the beam decreases by nearly 91.1% compared with one without the axial load. The lumped mass reduces the natural frequency of the structure. When the lumped mass equals the mass of the beam, the natural frequency decreases by 32.8%. When the NSGT are considered, the critical buckling load depends on the values of the nonlocal parameters and the strain gradient constant. In addition, the peak value of the RMS voltage output of the postbuckling process is significantly larger than that of the prebuckling process, which provides a new idea to improve the mechanical electrical conversion efficiency of the structure. Furthermore, the size-dependence of the flexoelectric effect on the performance of the energy harvester is demonstrated, indicating that the flexoelectric effect can significantly increase the RMS voltage output of the energy harvester. The results have some reference value for the design and application of piezoelectric energy harvesters.