Uncertainty Analysis of Piezoelectric Vibration Energy Harvesters Using a Finite Element Level-Based Maximum Entropy Approach

Piezoelectric energy harvesting has received significant attention in recent years as a viable solution to self-powered wireless sensors for emerging applications including wearable electronics, structural health monitoring, internet of things (IoT) and robotics. Ambient vibrations are abundant in many applications, for example, industrial machines, moving vehicles and aircraft, building and bridges, and human motions. Piezoelectric materials have a crystalline structure inside which the atoms are not symmetrically arranged. Deforming the structure modifies the balance of the electric charges and results in a net electric charge on the crystal surface, which is called the direct piezoelectric effect. This effect has typically been used for sensor devices such as accelerometers, microphones and load cells, and also makes piezoelectrics a suitable material for vibration energy harvesting, where the piezoelectric material is strained as a result of vibration. In addition, piezoelectrics have attractive features such as high energy density and compact and simple architecture [1]. These all contributed to the growing interest in piezoelectric vibration energy harvesting for self-powered microsystems. A general overview of the research and development in piezoelectric vibration energy harvesting can be found in Refs. [2–8].

Most of the analytical models of piezoelectric energy harvesters (PEHs) have been developed based on the assumption that the vibration excitation is harmonic. Also, it is desirable to make vibration energy harvesters lightly damped to utilize their large structural response for greater power output. As a result, the effective harvesting bandwidth is usually narrow and the power output is very sensitive to the “matching” between the excitation frequency and the natural frequency of the system. Though nonlinearity has been introduced to broaden the bandwidth and reduce this sensitivity [9,10] and has shown to overperform the linear configuration by an uncertainty propagation study [11], the uncertainty in the system still has a significant influence on the power performance. Recognizing the uncertainty in most environmental vibration excitations, i.e., variation of amplitude and frequency, researchers have studied the effect of excitation randomness on harvested power and attempted to developed models in a stochastic manner. Lefeuvre et al. [12] performed theoretical and experimental studies, and compared the power performance of the standard AC–DC and synchronous electric charge extraction techniques in the case of broadband, random vibration. Halvorsen [13] developed a closed-form model of linear resonant energy harvesters driven by broadband vibrations, and obtained the Fokker–Planck equation describing the probability density function of the harvested power. Following a similar approach, Adhikari et al. [14] derived expressions for the mean normalized power of system of stack configuration subjected to Gaussian white noise base acceleration. They studied the cases when the system was connected to a resistive load and connected to a resistive load and inductor in parallel. Seuaciuc-Osório and Daqaq [15] presented a theoretical analysis of the response of energy harvesters to excitations having a time-varying frequency. Yoon and Youn [16] applied a statistical time–frequency analysis to quantify the harvested power of a piezoelectric energy harvester under nonstationary random vibrations. Based on a distributed-parameter electro-elastic formulation, Zhao and Erturk [17] presented analytical and numerical solutions, and experimental validations of piezoelectric energy harvesting from broadband random vibrations. Based on the Wiener path integral technique, a methodology was developed by Petromichelakis et al. [18] to determine and optimize stochastic response of nonlinear electromechanical energy harvesters.

On the other hand, to enhance the power or energy conversion performance of the system and provide design guidance for energy harvesters, optimization studies have been conducted on the geometrical parameters [19–22] or the electrical parameters [23–26]. Additionally, topology optimization methods have also been applied [27–29]. However, as pointed out by Franco and Varoto [30], most of these optimization strategies usually seek for a single optimal parameter at a time. A piezoelectric energy harvester is an electromechanically coupled system whose performance is simultaneously affected by multiple parameters in the materials, mechanical, and electrical domain. Moreover, during the modeling and manufacturing processes, uncertainty is inevitably introduced into the mean model, for which the optimal parameters are obtained. As a result, the actual performance will deviate from predictions.

While some optimization methods have been coupled to stochastic response analysis tools for the response optimization of energy harvesters subjected to external random excitation, e.g., in Ref. [18], researchers have also recognized the importance of quantifying the uncertainties in the system parameters and further accounting for them in the analysis and design processes of harvesters, e.g., see Refs. [11] and [30–39]. These investigations have used parametric uncertainty, i.e., they have considered variations in some of the parameters of the system which they have modeled as random variables, independent of each other when multiple such parameters are varied. Geometrical dimensions, properties of the structural component and/or of the piezoelectric device, and characteristics of the electrical circuits have been treated as random. Moreover, the structural model of the harvester has typically been an equivalent single-degree-of-freedom (SDOF) although finite elements have also been used [36,39]. These investigations have yielded two key findings. First and foremost, they have demonstrated that a small level of uncertainty typically induces a dramatic change in power output. Relatedly, they have also shown the interest of designing under uncertainty. That is, the performance of the harvester optimally designed based on the mean model is no longer optimum once uncertainty is introduced. These findings clearly demonstrate the need to include uncertainty in the analysis and design of piezoelectric energy harvesters and the focus of the present investigation is on developing a general framework to consider in such efforts the uncertainties on the geometry/material properties of the structural and piezoelectric components.

While the above investigations have focused on introducing uncertainty on certain parameters of simple structural models of the harvesters, the real-world uncertainty is expected to affect all components and parameters of the model in a coupled manner. For example, uncertainty in a natural frequency originates from uncertainty on the structural properties and/or geometry and thus also implies uncertainty in the mode shapes, in the coupling with the piezoelectric elements, as well as in other natural frequencies. To account for all these effects, it is necessary to have a global, or nonparametric, modeling of uncertainty and such a modeling is best performed on a finite element description of the harvester, i.e., structure and piezoelectric components, especially to enable the consideration of complex geometries. Such an uncertainty modeling is proposed here based on the very recent work in Refs. [40–42], where the developed nonparametric approach is nonintrusive to the finite element software, requiring only the capability to output finite element mass, stiffness, electromechanical coupling, and piezoelectric capacitance matrices for the mean harvester design. The approach in Refs. [40–42] follows the original work of Soize in Ref. [43], in which the maximum entropy principle is applied and the uncertainty is modeled directly at the level of a reduced-order model of the structure, not within the finite element formulation as is done with the stochastic finite element approach. In this nonparametric maximum entropy approach, the uncertainty is lumped into the mass and/or stiffness matrix of the structure corresponding to a specified, deterministic basis. The joint probability density function of the elements of these matrices is then selected to: (i) have means which equal the corresponding deterministic matrices of the mean structural model, (ii) satisfy the mathematical requirements existing for these matrices (symmetry, positive definiteness), (iii) satisfy an integrability condition which guarantees the existence of the response in mean square, and (iv) maximize the entropy. This maximum entropy approach has been applied to matrices and fields with a variety of properties for which the construction of realizations has been detailed and is often very straightforward permitting a broad set of applications, see Refs. [44] and [45] for extensive review. The specific formulation in Refs. [40–42], also used in this study, is a modification of the original methodology of Soize, in which the elemental mass and stiffness matrices are recognized as reduced order model matrices of the structure within each of the finite elements with the interpolation functions serving as basis functions inside the elements. Once the mean model elemental matrices have been obtained, uncertainty is introduced element per element for each sample of the uncertain harvester, the structural-piezoelectric matrices are reassembled and the simulation, e.g., output power determination, can be carried out. Random elemental matrices are simulated in such a way that all sources of uncertainties, material properties and geometry, are accounted for as long as the components remain linear. The uncertainty modeling then becomes equivalent to the random field modeling of these elemental matrices, which is described in full in Sec. 3 of the paper. For the broadest applicability, the use of finite element software that does not include an explicit piezoelectric element is also considered and it is shown by extending the thermal analogy, that the modeling of the complete system can still be accomplished if the finite element model includes thermal effects.

In this study, the above approach is applied to a bimorph PEH, for which an analytical model and an optimization formula have been developed [24] to obtain the optimal load resistance for the specified maximum power output at a given excitation frequency. A finite element model of the energy harvester is first constructed with Nastran using the piezoelectric-thermal analogy and the uncertainty modeling is implemented to allow uncertainty in structural and/or piezoelectric properties. These capabilities are then utilized to study the effects of uncertainty on the power output of the energy harvester, both with weak and strong coupling and with different tuning scenarios of the resistive load to optimize the power output. Section 2 provides an overview of important power characteristics of piezoelectric vibration energy harvesters, along with their finite element modeling, including the thermal analogy necessary in Nastran, the finite element software chosen for this analysis. The uncertainty modeling is reviewed in Sec. 3 and the application to the bimorph PEH is detailed in Sec. 4. Section 5 concludes the paper.

which varies as the excitation frequency changes. If the load resistance is optimally tuned for all frequencies, this results in the power envelope of the system shown in Fig. 1(a), representing the maximum possible power through the tuning of load resistance. Graphically, this envelope is essentially the outer profile of the power curves if we plot them for a continuous and infinite range of load resistance instead of the five resistance values.

which is a function of the mechanical damping ratio. A system of coupling higher than the critical coupling is defined as strongly coupled, and it is weakly coupled if the coupling is lower than the critical coupling. For damping ratio 0.02, the critical coupling coefficient is 0.0816, and the curves for k² = 0.18 and 0.25 are of strongly coupled systems. There are two power limit peaks in the power envelope of a strongly coupled resistive energy harvester: one near the short-circuit natural frequency and the other near the open-circuit natural frequency. Note that the critical coupling changes with the type of energy harvesting circuit interface [47], which has been utilized as a method for enhanced power performance through innovative circuit designs. However, the system is still subjected to the same power limit regardless of the interface circuit type [49].

Even though once the coupling reaches its critical value, the power saturates and further increasing of the coupling does not lead to enhanced power; a strongly coupled system still offers few benefits. For example, the higher coupling induces more damping and further reduces the structural response and stress, which helps extend the fatigue life of the system. In addition, it can be seen from Fig. 2 that the frequency bandwidth, over which the harvested power is relatively large, is much wider at high coupling. This allows the system to be more robust as to the change in the excitation source frequency through a “correction” tuning of the electrical load to match the actual frequency.

where C_p is the total piezoelectric capacitance of the PZT layer, obtained by the summation of all the elements of the capacitance matrix, $Cp(FE)$⁠, and q is the total electric charge induced on the layer, obtained by summing the columns of the charge vector, $q¯(FE)$⁠.

Substituting this relationship into Eqs. (12a) and (12b) yields the governing equations

which serves as the governing FE equations to be solved.

To perform the uncertainty analysis at the finite element level, the FE matrices in Eqs. (11a) and (11b) need to be obtained. A convenient way is to extract them from a commercial fe software. However, not all commercial finite element software allows directly modeling of piezoelectric element, which prevents piezoelectric matrices from direct extraction. This is the case, for example, of MSC Nastran, which is a widely used multidisciplinary finite element package that is capable of performing static, dynamic, and thermal analysis of structures in both linear and nonlinear domains. Nevertheless, it is possible to resort to some analogy between piezoelectric and structural or thermal properties to extract these FE matrices. This is implemented in this study, described as follows.

In general, the FE matrices to extract can be classified into three groups: (a) structural matrices; (b) piezoelectric-structural coupling matrix; and (c) piezoelectric matrix.

• The structural mass and stiffness matrices, $M(FE)$ and $K(FE)$⁠, are output by using the Nastran DMAP alters, the high-level Nastran language commands. The damping matrix $C(FE)$ is then computed by using the Rayleigh damping model, i.e.,

  $C(FE)=αM(FE)+βK(FE)$

  (16)

  where $α$ and $β$ are the Rayleigh damping coefficients, determined according to the damping ratio of the mode of interest. The Rayleigh damping model is a widely used viscous damping model where the damping is considered to be associated with mass and stiffness and expressed as a linear combination of the mass and the linear stiffness matrices, which is convenient to use with a finite element model or in the modal expansion approach.
• The electromechanical coupling matrix $Θ(FE)$ is obtained by considering the following static actuation equation

  $f¯(FE)=Θ(FE)v¯(FE)$

  (17)

  which shows that when a static distributed voltage $v¯(FE)$ (a $NE$-by-1 vector) is applied to the piezoelectric layer, it induces a static distributed force $f¯(FE)$ (a $NS$-by-1 vector) on the structure. From Eq. (17), one can imagine when the static voltage is such that its i_th element is one and the other elements are zeroes, the static structural force will be equivalent to the i_th column of the matrix $Θ(FE)$⁠. In this way, varying i from 1 to $NE$ one by one and collecting the corresponding structural force vectors, $Θ(FE)$ can be obtained. With Nastran linear static solution (SOL101), the structural force vector can be obtained as the reaction force by applying the static voltage and fixing all the structural degrees-of-freedom, where the static voltage is applied by the equivalent distributed temperature using the piezoelectric-thermal analogy [50,51].To clarify this analogy, consider a piezoelectric beam harvester with direction 1 defined along the length direction and direction 3 defined along the thickness direction (also the polarization direction). Due to the converse piezoelectric effect, the induced normal or bending strain in the one-direction, i.e., ε₁, due to an applied electrical field in the three-direction, E₃, is given as

  $ε1=d31E3$

  (18)

  where d₃₁ is the piezoelectric coefficient. Usually, it is a negative number, meaning that application of a positive electric field in the polarization direction, i.e., direction 3, will generate a compressive strain in direction 1. Assume the electrical field within the PZT element is uniform and given as E₃ = v/t_p, where v is the voltage across the element (between the electrodes on the top and bottom surfaces) and t is the thickness of the PZT. Equation (18) can be rewritten as

  $ε1=d31tpv$

  (19)

  On the other hand, the thermal strain induced by a temperature change is given as

  $ε1=α1(T−Tref)$

  (20)

  where α is the coefficient of thermal expansion, T the current temperature, and T_ref is the reference temperature.Comparing Eqs. (19) and (20) and setting T_ref = 0 leads to the piezoelectric-thermal analogy

  $v=T,α1=d31tp$

  (21)

  which means that the converse piezoelectric effect can be equivalently modeled by the thermoelastic effect of structures with the temperature corresponding to the applied voltage and the equivalent thermal expansion coefficient related to the piezoelectric properties as shown in Eq. (21).
• The piezoelectric-thermal analogy only models the converse piezoelectric effect only in the mechanical domain, i.e., Eq. (11a) or (14a), thus can only be used to extract the electromechanical coupling matrix $Θ(FE)$⁠. To model the direct effect on the circuit dynamics in the electrical domain and set up the other governing Eq. (11b) or (14b), the piezoelectric capacitance matrix needs to be extracted as well. This can be achieved by an analogy between the structural mass, M, and the total piezoelectric capacitance, C_p, expressed as, respectively,

  $M=∫VSρSdvS$

  (22a)

  $Cp=∫VSψT(z)εSψ(z)dvS$

  (22b)

  where $ρS$ is the mass density, $εS$ is the dielectric constant at constant strain, and $ψ(z)$ defines the distribution of the piezoelectric material across the thickness of the device domain, $VS$⁠. From Eqs. (22a) and (22b), the capacitance matrix $Cp(FE)$ can be obtained as the mass matrix of the structure, when the mass density is set to be

  $ρsE=ψT(z)εSψ(z)={ε0(ε33−d312ε33Y11)⋅1tp2 piezoelectric layers0 substrate$

  (23)

  where $ε0$ is the vacuum permittivity, $ε33$ is the relative permittivity in the three-direction, $Y11$ is the Young's modulus in the one-direction, and $tp$ is the thickness of the piezoelectric layer.

With the finite element modeling framework constructed above, fully coupled piezoelectric-structural simulations of piezoelectric energy harvesters can be performed. For validation, a finite element model of a bimorph harvester of properties given in Table 2 was constructed using Nastran, as shown in Fig. 3. The harvester is basically a cantilever beam type structure, and CQUAD4 shell elements were used for the structure FE model with 1760 elements and 1887 nodes. The composite property card (PCOMP) was used to define the substrate and the piezolayers structural properties.

The solution is implemented outside of the Nastran environment. At each excitation frequency, the amplitude of harvested power output is computed from the solution as $P=(ωq0)2R$⁠, where $q0$ denotes the amplitude of the charge.

The harvested power versus excitation frequency is plotted in Fig. 4(a) at different load resistances. It is the same configuration used by Liao and Sodano [24] in their experimental studies and the results are shown in Fig. 4(b). For further comparisons, its ansys and analytical results are shown in Figs. 4(c) and 4(d), respectively. The analytical results are obtained by using Eq. (5) with the effective system parameters evaluated by using the expression in the Appendix. Overall, it can be seen that the NASTRAN derived results are in excellent agreement with other results, which confirms the adequacy of the thermal analogy.

In this study, the finite element-level nonparametric maximum entropy approach is used for uncertainty analysis of the bimorph piezoelectric energy harvester. Based on the piezoelectric-thermal analogy, the approach is implemented in a similar way to that in the uncertainty analysis of heated structures [42], where the details of the theoretical derivation can be found. The mass, stiffness, and the electromechanical coupling matrices of the FE model are randomized in the present uncertainty analysis to reflect the uncertainties of the parameters of the coupled system globally. Note that the ordering of the degrees-of-freedom in the matrices discussed below is degree-of-freedom 1 for all nodes, degree-of-freedom 2 for all nodes, etc.

where $H(e)$ is a m-by-m random matrix where m is the number of degrees-of-freedom per node and its structure is shown in Fig. 5, $Ir$ is the r-by-r identity matrix where r is the number of nodes per element, and ⊗ denotes the Kronecker product operation.

where F is the cumulative distribution function of the standard Gaussian random variable, $FHii−1$ is the inverse of the cumulative distribution function of the Gamma random variable. As shown in Fig. 5, $σ=1/2μ$ where $μ$ is the parameter controlling the uncertainty level. Usually, an alternative dispersion parameter $δ$ is used (as in this study), which is related to $μ$ as $μ=(n+1)/(2δ2)$⁠, where n is the dimension of the matrix. Once random samples of all elemental matrices are obtained, they are assembled to construct the random sample of the mass matrix of the whole FE model.

where $K(e)$ and $Θ(e)$ are the $nS$-by-$nS$ structural stiffness and the $nS$-by-$nΘ$ electromechanical coupling matrices of an element, respectively. Here, $nS$ and $nΘ$ are the numbers of degrees-of-freedom of structural displacement and electrical voltage per element, respectively. In the present application, $nS=6×r$ and $nΘ=1×r$ where r = 4 as the CQUAD4 elements have four nodes, with six structural degrees-of-freedom, and one electrical degree-of-freedom each. Finally, $KEE(e)$ denotes a $nΘ$-by-$nΘ$ matrix relating to the strain energy induced by the electromechanical coupling, see Ref. [40] for discussion in the thermal case, and which will be discussed later in this section. The randomization of the elemental stiffness and the elemental electromechanical coupling matrices is carried out with the mean $KB(e)$ matrix, denoted as $K¯B(e)$⁠.

that is, $Φ¯D(e)$ and $Λ¯D(e)$ are the matrices of eigenvectors and eigenvalues (diagonal elements) of the matrix $Θ¯(e)$⁠.

The third matrix, $L¯EE(e)$⁠, should be obtained from the matrix $K¯EE(e)$⁠. The thermal counterpart of $K¯EE(e)$ has been shown [42,53] to be difficult to accurately estimate from a commercial finite element code nonintrusively. Hence, considering that it does not appear in the governing equation, it was proposed [42,53] that it be selected to maximize the entropy of the random samples, achieved by setting $L¯EE(e)$ to be the identity matrix.

where $HEE(e)$ only affects the $KEE(e)$ block. Since the $KEE(e)$ matrix does not appear in the governing equation, its value is irrelevant and thus $HEE(e)$ neither needs to be computed nor discussed. It is symbolically replaced by a * in the following.

where the notation $A,i$ denotes the ith row of the matrix A and diag is the operation taking a vector and creating the diagonal matrix having these elements along the diagonal. Moreover, in Eq. (43), $hEK(e)$ is a matrix of $nΘ$-by-6 components defined as independent Gaussian random fields with a specified autocorrelation function, e.g., Eq. (29), and $Jr$ denotes the row vector of dimension r with all elements equal to one.

Once random samples of all elemental matrices are obtained, they are assembled to construct the random sample of the stiffness and the coupling matrices of the whole FE model. For each random sample, the damping matrix is computed by the Rayleigh damping model using the randomized mass and stiffness matrices and the same Rayleigh damping coefficients as in the nominal FE model. Therefore, the variation of the damping ratio is introduced through the variations of mass and stiffness.

where $Lcorr$ is the corresponding correlation length. While each random field can be described by its own correlation length, they were all taken equal here to 1/10 of the beam length to exemplify the methodology.

To illustrate the uncertainty analysis process and investigate the effect of coupling and load tuning schemes on power uncertainty, two numerical beam configurations shown in Table 3 are studied. System 1 is weakly coupled with coupling coefficient k² = 0.03, and system 2 is strongly coupled with k² = 0.15. There are two sets of uncertainty studies performed for both systems: (1) uncertainty in the harvested power at a given excitation frequency and (2) uncertainty in the power envelope. We also compare the situation when the electrical load remains fixed to its pretuned or pre-optimized value to the situation when the load is “adaptively” retuned to the modified system parameters due to the introduced uncertainty. As a comparison baseline, the two systems harvest the same amount of power at a target excitation frequency, i.e., 3 mW at 135 Hz. The analytical model discussed in Sec. 2.1 has been used to guide the selection of these system properties in Table 3. The associated power behavior of the systems is shown in Fig. 6, which shows that the analytical model predicts that the two systems will harvest 3 mW at 135 Hz with an electrical load of 6001 Ω and 2806 Ω, respectively. To provide a “nominal” model for the uncertainty analysis, the two systems are also analyzed by NASTRAN using the approach outlined in Sec. 2.2 and the results are shown in Fig. 6 as well. Overall, the NASTRAN results match the analytical results quite well, except for a slightly frequency shift (about 2 Hz or 1.5% with respect to 135 Hz) between the results.

For these two configurations, the mass, stiffness, and the electromechanical coupling matrices of the FE model are randomized to globally reflect the uncertainties of the parameters of the coupled system. Since no detailed uncertainty information is available for the system parameters, the uncertainty levels of the three matrices, denoted as $δM$⁠, $δK$⁠, and $δΘ$⁠, respectively, are selected to yield specified values of the coefficient of variation of three quantities, selected here to be the total mass, the natural frequency of the first mode, and the modal coupling coefficient of the first mode, i.e., $θ1=ϕ¯1(FE)T⋅θ¯(FE)T$ where $ϕ¯1(FE)$ is the first mode. A population of 100 random samples was generated with coefficients of variation of the above three quantities set to 1.0%, 1.0%, and 1.5%, respectively, following the discussion of Sec. 3.

If some detailed information about the uncertainty in the parameters of the energy harvester is known, it can be used to determine the above uncertainty levels. For example, if the variation of the mass density is known, a number of samples of the total structural mass can be generated by a traditional stochastic approach, e.g., the Monte Carlo method, from which the true variation of the total mass is computed. Then, for each of a set of uncertainty levels $δM$⁠, a number of random samples are generated by the nonparametric approach, and the corresponding variation of total mass is computed. The uncertainty level at which the variation of total mass matches the true variation is taken. Alternatively, a more rigorous approach is to evaluate the likelihood function using the sample data for each uncertainty level, with respect to the true sample data; then the uncertainty level is determined as the one corresponding to the maximum likelihood [54,55]. The above procedure can be applied to $δK$ if the variation of the Young's modulus is known; thus, the variation of the natural frequency of the first mode can be computed, and to $δΘ$ if the variation of the piezoelectric coefficient is known; thus, the variation of the modal coupling coefficient of the first mode can be computed.

The uncertainty analysis results are presented in Fig. 7 in the form of uncertainty bands spanning the 5th to 95th percentile interval. For both beam configurations (weakly and strongly coupled), there are 100 random samples at each nominal point, i.e., for each particular frequency of excitation in the band 115 Hz to 160 Hz. For each of the 100 samples, the response to the base excitation of 1 g acceleration at the nominal point is computed and the harvested power is determined. These computations are performed for two situations: (1) the load resistance is fixed to its pre-optimized (pretuned) value of the nominal model for all 100 samples and (2) the load resistance is “adaptively” optimized (adaptively retuned) for each sample individually to account for the system change due to uncertainty using based on Eq. (8). The first scenario corresponds to a harvester produced as a black box ready to be installed for a particular application. The second scenario would correspond to a critical use of a harvester where availability of maximum power is more important than the additional cost associated with the tuning of the load.

In the case of preoptimized resistance, as seen from Figs. 7(a) and 7(b), the uncertainty effect on the power envelope is significant. For both configurations, the uncertainty band is quite broad with the relatively small uncertainty level, and the effect on the weakly coupled configuration appears to be larger than that of the strongly coupled one. The uncertainty does not appear to affect the overall maximum value of the harvested power as much, which is consistent with the piezoelectric energy harvesting theory that the overall maximum power of an energy harvester is “capped” by its system properties, for example, Eq. (9). The change in system properties due to uncertainty is not significant enough to result in a large change in the overall maximum power. On the other hand, the uncertainty analysis shows that there could be a large probability that an energy harvester sample does not yield this maximum power. This means that a factor of safety on power needs to be considered in the design process, and the uncertainty analysis serves as an effective tool to quantify it.

In the case of adaptively reoptimized resistance, as seen from Figs. 7(c) and 7(d), the uncertainty band of the weakly coupling configuration does not change as much; while that of the strongly coupling configuration is reduced by almost half in the frequency range between the two power limit frequencies, over which the power level is high. This is consistent with the piezoelectric energy harvesting theory that higher coupling leads to a broader harvesting bandwidth. The uncertainty in the system changes the system parameters and graphically moves the power envelope left and right (with small change in power). For a weakly coupled system, this movement could lead to a significant power drop at a particular frequency because of its narrow power envelope; while the effect is not as significant for a strongly coupled system because of the wider frequency range of high power. In all, the uncertainty results and observations suggest that the configuration of strong coupling be preferred in the design, since it is less sensitive to the uncertainty as shown in Figs. 7(a) and 7(b); and when the uncertainty is present, its performance can be improved, e.g., by the adaptive optimization, as shown in Figs. 7(c) and 7(d).

In the above study, the uncertainties in mass, stiffness, and coupling are all included. To examine their effects separately, the configuration in Fig. 7(b) was used as an example and analyzed in two scenarios: presence of uncertainty in mass and stiffness only, and presence of uncertainty in coupling only. The results for the two scenarios are presented in Figs. 8(a) and 8(b), respectively. For the system and under the uncertainty settings, it can be seen that the overall uncertainty of power shown in Fig. 7(b) is mostly due to the uncertainty in mass and stiffness. In addition, the effect from the mass and stiffness uncertainty appears to be “centered” about the nominal power envelope through most of the frequency range; while the effect from the coupling uncertainty alone is almost unnoticeable on the left side of the power envelope and becomes much more visible on the right side. These observations are consistent with the analytical PEH power characteristics discussed in Sec. 2.1. For Fig. 8(a), the uncertainty in mass and stiffness changes the natural frequency of the system, effectively moving the power envelope left and right about the original power envelope. In the middle frequency range between the two power limit peaks, the maximum power is still subjected to the power limit, and the pretuned load resistance becomes “mistuned” to the “new” system in the presence of uncertainty. As a result, almost all the power points fall below the nominal power envelope curve over the middle frequency range. On the other hand, the change in the coupling of a strongly coupled system does not change its short-circuit natural frequency but open-circuit natural frequency. As illustrated in Fig. 2, when the coupling increases, the left side of the power envelope (near the short-circuit natural frequency) remains almost “fixed” while the right side extends to the right as a result of increased open-circuit natural frequency. Therefore, the effect of uncertainty in coupling causes the power envelope to extends and contracts mostly on the right side. This is in agreement with the results shown in Fig. 8(b), where the uncertainty band is very noticeable on the right side but not on the left side of the power envelope.

To investigate further, the role of adaptive optimization in handling the uncertainty effect, and the distribution of the harvested power (samples) at a few selected excitation frequencies by the systems shown in Fig. 7 are presented in Fig. 9. The adaptive optimization does not offer much power improvement as to uncertainty for the weakly coupled system throughout the frequencies, e.g., Figs. 9(c), and 9(e), except some slight enhancement near the peak frequency, e.g., Fig. 9(a). This can be also stated about the strongly coupled system for frequencies outside the frequency band between the two power limit frequencies, e.g., Figs. 9(d) and 9(f). There are only a very small number of samples whose harvested power is increased by the adaptive optimization. However, as shown in Fig. 9(b), for frequencies that fall into the band, the adaptive optimization is able to not only reduce the size of the uncertainty band (or scattering) but also boost the harvested power for most of the samples. In addition, the distribution of the optimized values of the load resistance is shown in Fig. 10. It can be seen that the enhanced power in Figs. 9(a) and 9(b) is connected to a more centered distribution of the adaptively retuned optimized resistance about its associated pretuned optimal value.

The adaptive retuning of the load resistance has been shown to be an effective method for power enhancement and uncertainty reduction of strongly coupled systems. However, compared to the pretuned method, it requires an additional tuning step and its implementation increases the complexity of the system because the load needs to be adjustable. As another effort to find a pretuned and fixed load resistance for power enhancement, a new tuning scheme was investigated. Instead of using the analytical optimization Eq. (8) to determine the pretuned optimal resistance, uncertainty analysis was performed to find an optimal resistance that maximizes the mean power of the samples [24]. As an example, the harvested power of the strongly coupled system at excitation frequency 140 Hz was studied. Uncertainty was introduced into 100 samples. For a fixed pretuned resistance, the harvested power was obtained for all the samples. The results are presented in Fig. 11, where the harvested power is plotted against the pretuned resistance. The plot clearly demonstrates the load dependence of power, a characteristic of PEHs. The pretuned R_opt given analytically by Eq. (8) based on the nominal model is 7863 Ω, and the numerical study finds the load resistance maximizing the mean power (of the 100 samples) is 8200 Ω. Table 4 summarizes and compares the power characteristics of the three load-resistance tuning schemes. The pretuned R_opt scheme determines a fixed load resistance analytically from Eq. (8) based on the nominal model; the adaptive retuned R_opt scheme adjusts the load resistance still analytically from Eq. (8) but based on the updated model due to the presence of uncertainty; and the pretuned-for-the-mean method finds the optimal load resistance through the uncertainty analysis of the nominal model. It can be seen from the table that the two pretuned methods perform about the same. However, the pretuned-for-the-mean method could be a more practical approach as the pretuned R_opt scheme relies on the analytical theory, i.e., Eq. (8), which might not be available or accurate due to the modeling complexity of energy harvesters. The adaptive tuning scheme outperforms the other two. However, as mentioned above, the load needs be adjustable and adjusted individually to account for different uncertainty effects in the samples.

A general nonparametric uncertainty analysis methodology is presented for uncertainty quantifications of piezoelectric energy harvesters at the finite element level. First, a finite element modeling framework is developed, which utilizes the matrices relating to the structural properties (mass, stiffness), the piezoelectric capacitance matrix, as well as the structural-piezoelectric coupling terms of the mean harvester. For finite element software that do not have piezoelectric elements but include thermal effects on structures, the thermal analogy linking piezoelectric and temperature effects can be applied. Fully coupled piezoelectric-structural simulations can be performed using the constructed finite element modeling framework. This framework was applied to a bimorph energy harvester in Nastran and validated by the excellent agreement between its results and those obtained experimentally, analytically and from ansys.

To better account for the coupled manner that the real world uncertainty is expected to affect all components and parameters of the model, a global, nonparametric uncertainty modeling approach based on the maximum entropy principle is adopted with the constructed finite element model, and applied to two configurations of the energy harvester, one of weak and the other of strong piezoelectric-structural coupling. For a fair comparison, they are designed for the same excitation frequency at the same power level. Random samples of mass, stiffness, and electromechanical coupling matrices are generated to represent the global effects of the uncertainties of the system parameters. The uncertainty analysis is performed to study the uncertainty effect on the power output of the energy harvester with three load resistance tuning schemes: First scheme is when a fixed load resistance is given by an optimization method applied to the nominal model without considering the system changes due to uncertainty; second scheme is when the effect of uncertainty is considered and the load resistance of each random sample is adaptively optimized to the changed system parameters; third scheme is when a fixed load resistance is designed to maximum the mean power of the entire sample set. The uncertainty effect is shown to be significant even at a relatively low uncertainty level, and more pronounced for the weakly coupled configuration than the strongly coupled one. The uncertainty analysis demonstrates the advantages of a strongly coupled harvester over a weakly coupled harvester. First, its harvested power is less sensitive to uncertainty with either pre-optimized or adaptively optimized load resistance. Second, with the application of the adaptive optimization, its uncertainty effect can be significantly reduced and the harvested power can even be boosted if the target excitation frequency falls within the “harvesting bandwidth” between the two power limit frequencies. On the other hand, the application of the adaptive optimization method does not affect the uncertainty characteristic of the weakly coupled harvester as much.

The financial support of this work to the first (XQW) and third (MPM) authors by the Air Force Multi-University Research Initiative contract FA9550-15-1-0038 with Dr. Jean-Luc Cambier as Technical Monitor and the Air Force Office of Scientific Research Grant No. FA9550-16-1-0021 with Dr. Jaimie Tiley as Technical Monitor are gratefully acknowledged.

• Air Force Multi University Research Initiative (FA9550-15-1-0038).

• Air Force Office of Scientific Research (FA9550-16-1-0021; Funder ID: 10.13039/100000181).

where the subscripts s and p denote the substrate and piezoelectric materials, respectively. The length, width, and thickness are denoted by L, b, and t, respectively. In addition, ρ is the density, Y is the Young's modulus in the one-direction, $K3S$ is the relative dielectric constant measured at constant strain, and d₃₁ is the piezoelectric coefficient.