Design Scalability Study of the Γ -Shaped Piezoelectric Harvester Based on Generalized Classical Ritz Method and Optimization

: This paper studies the design scalability of a Γ -shaped piezoelectric energy harvester ( Γ EH) using the generalized classical Ritz method (GCRM) and differential evolution algorithm. The generalized classical Ritz method (GCRM) is the advanced version of the classical Ritz method (CRM) that can handle a multibody system by assembling its equations of motion interconnected by the constraint equations. In this study, the GCRM is extended for analysis of the piezoelectric energy harvesters with material and/or orientation discontinuity between members. The electromechanical equations of motion are derived for the PE harvester using GCRM, and the accuracy of the numerical simulation is experimentally validated by comparing frequency response functions for voltage and power output. Then the GCRM is used in the power maximization design study that considers four different total masses—15 g, 30 g, 45 g, 60 g—to understand design scalability. The optimized Γ EH has the maximum normalized power density of 23.1 × 10 3 kg·s·m − 3 which is the highest among the reviewed PE harvesters. We discuss how the design parameters need to be de-termined at different harvester scales.


Introduction
The high demand for self-powered electronic devices, such as sensors and small electronics, has increased and energy harvesting is expected to gain more attraction. Piezoelectricity is one of the energy conversion principles used in energy harvesting systems, and it has shown advantages such as a high voltage output and ease of fabrication both in macro-and microscales [1][2][3]. These advantages have led to the use of piezoelectricity in various energy harvesting applications, and their design studies-structural modifications and electric circuit designs-have been actively conducted during the last two decades [4,5]. Previous piezoelectric energy (PE) harvesting studies considered various dynamic energy sources, such as HVAC (heating, ventilation, and air conditioning) systems (outdoor condensing unit, duct, compressor, pipe, fan belt cage) [6][7][8], human motions (vibration [9,10] and shoe compression [11][12][13]), pavements (or tiles) [14][15][16][17], vibrating bridges [18], wind-induced motions [19][20][21], raindrops [22,23], tires [24], and even inside human bodies [25][26][27][28]. The aforementioned works developed numerical models for PE harvesters and performed design studies to reduce the development cost while satisfying the power requirement.
There are three typical methods to develop numerical models for PE harvesters: the finite element method (FEM) [29,30], analytical methods [1,[31][32][33], and the classical Ritz method (CRM) [1,[34][35][36][37][38][39]. Each method has its advantages and disadvantages. The FEM is advantageous in modeling complicated features such as multiple materials, complex geometric profiles, and various boundary conditions. However, the FEM inherently needs a large number of degrees of freedom (DOFs) and is computationally inefficient [39][40][41][42]. The analytical method shows high computational efficiency, but it is limited to a system having simple structural configurations and a limited number of boundary conditions. The CRM is a preferred method for a model with a moderate geometry complexity and multimaterial configuration, e.g., a two-segment self-charging piezoelectric cantilever [34]. Besides, the CRM does not need many DOFs and has excellent computational efficiency compared to the FEM. Meirovitch and Kwak [39] showed that the convergence speed of the CRM is faster than that of the FEM in the modal analysis for a nonuniform rod in axial vibration. Morales [38] demonstrated how the CRM can be used for N-story framed structures and showed a faster convergence speed of the CRM than FEM. However, there are two disadvantages of the CRM: (1) admissible trial functions satisfying all the geometric boundary conditions of the system must be used to approximate the deformation field, and (2) for systems consisting of multiple bodies as in [38], the so-called kinematical procedure, which represents kinematical quantities (i.e., displacement, velocity, and acceleration) of a certain body based on the adjoining body, has to be included in a derivation process of equations of motion to satisfy the compatibility condition. Given that the aforementioned disadvantages prevent the CRM being used to model general systems having multiple structures and many PE harvesters that consist of multiple structures (e.g., L-shape [32,33], M-shape [43,44], zigzag [45][46][47], E-shape [48], multiple cantilevers [49], and Γ-shape in this study) have better performance than a cantilever PE harvester in terms of power output or broadband characteristic, the CRM needs to be upgraded.
Jeong and Yoo [41,42] developed a modified version of the CRM, so-called generalized CRM (GCRM), that always uses the same set of trial functions regardless of the system. Their studies focused on modeling and analyzing structures with geometrical nonlinearity. This study extends the use of GCRM for electromechanical linear PE harvesters with general geometry, by including the piezoelectric constitutive relations. For simplicity, the proposed modified GCRM for PE harvesters is referred to as GCRM-P hereafter. The efficiency of the proposed GCRM-P method enables its smooth adaptation to a shape optimization process for PE harvesters to enhance design compactness and durability while meeting power requirements. This paper performs a design optimization study for a specific harvester configuration named Γ-shaped PE harvester (ΓEH) that assembles two uniform beams at a right angle (90°) (Figure 1a). This configuration has uniform strain distribution along its vertical structure, whereas a cantilevered PE harvester (CEH) (Figure 1b) has nonuniform strain distribution. The ΓEH is analyzed using the GCRM-P model and experimentally validated by comparing the voltage frequency response. In the design optimization study, four different mass constraints (15 g, 30 g, 45 g, 60 g) are considered, to understand design scalability. Furthermore, the energy harvesting performance of its optimal design is evaluated based on the normalized power density (NPD) [4,50]. The goals of this study are summarized as follows: • Development of a GCRM-P model to predict linear electromechanical behaviors of PE harvesters having multiple structural members, • Experimental validation for the GCRM-P model in terms of energy harvesting performance, • Study of design scalability-design optimization of ΓEH under different mass scales, and comparison of the power output performance with the other recent PE harvester studies.

Generalized CRM for Piezoelectric Harvester (GCRM-P)
This section derives electromechanical equations of motion for a piezoelectric harvester system using GCRM-P. The system is assumed to exhibit 2D linear stress−strain behavior and the material properties are assumed to be constant. We also assume that the beam is so slender that the rotary inertia effect of the cross-section is negligible. In the piezoelectric beam, the perfect bonding condition is assumed between the piezoelectric and substrate layers. Figure 2 shows the configuration of a unit piezoelectric beam element (before and after deformation) as well as its electrical connection to a load resistor. The beam element is a composite having both piezoelectric layers (top and bottom layers) and substrate layer

Electromechanical Energy Formulations
According to [1], the potential energy term is the width of the beam, and where ( ) k y is the location of the generic point along ( ) The kinetic energy term , sin cos where • denotes time differentiation and ( ) k η is the mass per unit length given as The internal electrical energy in the piezoelectric layers is given as is the equivalent internal capacitance of the piezoelectric layer given as and ( )

, 33
S k ε is the permittivity of the piezoelectric material at constant strain.

Constraint Equations
The GCRM uses constraint equations and the Lagrange multiplier method to meet compatibility and boundary conditions. The constraint equations for the unit element in Figure 1a are given as The constraint equations impose the continuity of DOFs across the neighboring elements. The boundary conditions can be applied by assigning corresponding values to Equations (15)- (20). For example, a fixed-free boundary condition can be implemented by (15), (17) and (19).

Spatial Discretization and Electromechanical Equations of Motion
The deformation variables are trial functions. The choice of the trial function affects the rate of convergence, computational time, and numerical stability [51]. Different from the CRM, the GCRM can use nonadmissible trial functions. Among others, this study selects the Legendre polynomials ( Figure 3) that were successfully used for both GCRM-based static [41] and dynamic [42] analyses. A mathematical expression for the Legendre polynomials is given as follows: where ! denotes the factorial of the given integer. According to the electromechanical Lagrange equations of motion with Lagrange multipliers, we have where e N and n N are the total numbers of elements and nodes. In Equation (24), ξ represents a coordinate in the system, i.e., γ represents either a mechanical forcing function (when ξ is a coordinate for displacement) or an electrical charge output ( ( ) ( ) . It is noted that the mechanical forcing effect by ground movement is already included in the kinetic energy formulations.
Between the two types of usual damping mechanisms-viscous air damping and strain-rate damping-only the latter one is considered because the system is assumed not to be operated in viscous fluids or microscale devices. The strain-rate damping effect is modeled by adding the modal damping matrix to the derived equations of motion [1,52]. Following previous studies [31,36,53], the modal damping ratio of 0.025 is employed for the first mode response.

Uniform Strain Distribution in -Shaped Structure
The CEH and ΓEH consist of piezoelectric and substrate layers and tip mass. In this section, however, the simplified structures without the piezoelectric layers and tip mass are analyzed for the simplicity of the analysis (detailed analysis and design are conducted in Section 4). Figure 4a shows the simplified CEH whose length, Young's modulus, and the second moment of the cross-section area are C L , C E , and C I . C x is the coordinate of the arbitrary point from the origin point Figure 4b shows the simplified ΓEH whose vertical length, horizontal length, Young's modulus, and the second moments of the vertical and horizontal beams' cross-section areas are Γ,V  This section performs a simply static analysis to explain how a uniform strain is generated in ΓEH. In both structures in Figure 4, the tip load F is applied so that it causes the tip displacements, C δ and Γ δ , respectively. The internal moments from each structure are given as Based on the Euler−Bernoulli beam theory, the curvatures are represented as where w is the transverse displacement and the superscript ( ) '' means that the variable is differentiated with respect to the spatial variable ( C x or Γ,V x ) twice. Substituting Equations (25) and (26) into Equations (27) and (28) gives To conclude, the bending curvature of the cantilever ( C '' w ) linearly decreases, whereas that of the vertical beam in Γ-shaped structure ( Γ,V '' w ) is constant and uniform.
The uniform curvature results in a uniform strain distribution and higher power generation. It is noted that one previous study-the L-shaped PE harvesters [32,33]-had a similar angled geometry, but had a different contribution to increase the bandwidth. It is necessary to validate the uniform strain distribution of the Γ-shaped structure by finite element analysis (FEA). For the FEA analysis, a commercial software ANSYS v19.2 (Ansys, Inc., Canonsburg, PA, USA) was employed to conduct a static analysis. In the FEA model, 1500 and 3764 quadratic hexahedron solid elements were used for the simplified CEH and ΓEH, respectively. The specifications in terms of geometry and material are not given here because the point of this FEA is to show the uniform strain distribution of the proposed structure rather than values. Figure 5 shows the strain distributions, boundary conditions, and loading conditions of the FEA models. The color in Figure  5 shows the magnitude of the von Mises strain induced by the applied static load F (the red and blue colors mean high and low strains). It is noted that the simplified CEH had the nonuniform strain distribution in its longitudinal direction as predicted in the analytical solution (Equation (29)), i.e., its strain varied from high (red) to low (blue) along the longitudinal direction (Figure 5a), whereas the simplified ΓEH had the uniform strain distribution (Figure 5b) in its vertical beam.

Configuration of -Shaped Harvester
The configuration of the ΓEH is given in Figure 6. The PIC151 ceramic (PI-Ceramic GmbH, Lederhose, Germany), polylactic acid (PLA), and lead were used for the piezoelectric layers, substrate layers, and tip mass, and their material properties (Table 1) were assigned to the numerical model. The piezoelectric layers were placed only in the vertical part where a uniform strain distribution was observed.
The GCRM-P model does not need the serious mesh generation as required in an FEM model, but only a few nodes where there is a change in material property, crosssectional property, or orientation of a structural member. Considering that the vertical beam is not fully covered by piezoelectric material, there needs to be a dummy element (which is rigid and massless) to place a gap between vertical and horizontal beams. The tip mass is modeled as a point mass at the center of the tip mass. Therefore, a node has to be defined there and is connected to the tip of the horizontal beam with a dummy element. Accordingly, five nodes (boxed numbers) and four elements (circled numbers) ware assigned as shown in Figure 6. In the design problem of this study, there are nine ( where A and f are 0.16 g and 40 Hz, respectively, and ( ) X V t  is assumed to be zero.

Shape Optimization
This section introduces the experimental validation of the GCRM-P model developed in Section 2 and the design formulation for the ΓEH. In the experimental validation, the accuracy of the GCRM-P model was verified by comparing the voltage output FRFs from the experiment and simulation. With the experimentally validated model, the shape optimization was performed and the design change trend was discussed when the total mass of the harvester changed.

Vibration source
(a)

Experimental Validation of the GCRM-P Model for EH
In this section, the simulation analysis result is validated experimentally, using the ΓEH prototype as specified in Table 2. For prototyping, the substrates of the ΓEH were fabricated part by part using a 3D printer (Prusa I3 MK3, Prusa Research, Prague, Czech Republic), with a 100% infill density, and assembled afterward. Then, the piezoelectric layers with the electrode and tip mass made of the lead were glued on the substrate using the cyanoacrylate adhesive.
The experimental setup in Figure 8 was prepared to measure the frequency response function (FRF) for the voltage and power output of the fabricated prototype. An electromagnetic shaker (Type 4809, Brüel & Kjaer, Naerum, Denmark), a function generator (33250A, Agilent, Santa Clara CA, USA), and a power amplifier (Type 2718, Brüel & Kjaer, Naerum, Denmark) were used to excite the energy harvesting beam, and the accelerometer in Section 4 was used to measure the acceleration of the base. The energy harvester was connected to a variable load resistance from which the voltage output was obtained using an oscilloscope (Tektronix DPO4054B, Tektronix, Beaverton, OR, USA). Then, the power output could be calculated by applying Ohm's law. The experimental and numerical FRFs for the voltage and power outputs are shown in Figure 9. The symbols and lines denote the experimental and numerical results, respectively. It is noted that the results are normalized by gravitational acceleration (g). The voltage outputs were measured at two load resistance conditions: one was the open-circuit condition ( oc R ) and the other was the optimized load resistance condition ( opt R ≈ 1/(ω ( ) 1 , P eq C )) given in Table 2, and the power output was measured only at the optimized load resistance condition. In both voltage and power results, the experimental FRFs were slightly shifted to the right-hand side. This stiffening effect could result from the thickness of electrodes in between the layers, unevenly distributed adhesive, substrates more thickly printed than intended, or a combined effect of some of the aforementioned possible causes. Nonetheless, we conclude that both results showed a reasonable agreementthe errors on the maximum open-circuit voltage output, voltage output at the optimized load resistance, and power output were about 11.4%, 2.6%, and 5.5%, respectively.

Design Formulation
The energy harvester design problem is formulated to maximize the peak value of power output. The design constraints consider cost-effectiveness and compactness as well as durability, and are formulated to limit: total mass ( total M ), horizontal length ( hor L ), ver- where PIC151 ρ , PLA ρ , and Lead ρ are density values of the materials given in Table 1. (38) is the amplitude of the downward nodal displacement at node 5.
The subscript " st dy + " indicates that the value is evaluated considering static (gravity effect) and dynamic (base excitation) loads. It is noted that the gravity effect can be removed or changed based on the installation orientation of the harvester. The value of 30 MPa in Equation (40) We select to use the differential evolution (DE) algorithm for design optimization [55,56], which evolves a pool of design candidate populations generation to generation. , ..., , 1, ..., . After the mutation operation, a crossover operation is employed to each pair of the target and mutant vectors , ..., , 1, ..., The binomial crossover is the most commonly used crossover operation given as if rand (0, 1) or , 1, ..., and 1, ..., otherwise where CR is a user-defined crossover rate (0 < CR < 1) that controls the fraction of vector components copied from the mutant vector, rand j (0, 1) is j th random real number in the range of [0, 1], and rand j is a randomly selected integer within the range [1, D].

After the trial vectors
where the function ( ) f ⋅ is a fitness function to be maximized. These three operationsmutation, crossover, and selection-are repeated until the termination criterion is met in terms of the number of maximum function evaluations.
In this study, we chose one of the DE variants, L-SHADE-EIG [57,58], combined with the oracle penalty method [59] to handle a constrained problem. This DE variant gained great attention as it won the first prize in a real-parameter single objective optimization competition at 2015 IEEE congress on evolutionary computation. In the algorithm implementation, the numbers of initial populations and the maximum function evaluation for the algorithm termination are set depending on the number of design variables (D=9) to be 10 D × and 1000 D × (or 90 and 9000) respectively, and other settings are adopted from the original L-SHADE-EIG algorithm. The fitness function is chosen as the peak value of power output obtained with the given input acceleration conditions.

Shape Optimization Results and Discussions
In this section, four different system mass values are attempted ( allowed M = 15 g, 30 g, 45 g, and 60 g) for the optimization and the corresponding design change trend is studied. The optimal design parameters at different allowed M values are found as shown in Table 3, Figure 10 (in the normalized scale), as well as Figure  11. The normalized values in Figure 10 were calculated by the following formulation.

Unnormalized value Lower bound Normalized value=
Upper bound Lower bound − − (45)   In all four cases, the mass constraint (Equation (32)) is found active, obviously indicating that a larger power is generated from a heavier harvester. More interesting observations include: (1)  h also increase so that the beam stiffness is increased and the natural frequency of the system stays close to the excitation frequency (40 Hz). It is noted from the optimization results ( Figure 10) that both the length and the width of vertical piezoelectric beam ( P L and b) increase but their increasing trends are different. as the stress violation is concerned-the bending moment of inertia is a cubic function of the thickness, which is more effective to decrease stress than the width (linear function). The tip mass parameters are increased correspondingly to maintain the harvester's resonant frequency. Several performance measures are compared from four different designs as shown in Figure 12. The stress values increase as a larger mass is allowed as shown in Figure 12a.  (37)) and deteriorate the durability of the energy harvester. Figure 12b shows the change of the piezoelectric material volume that increases gradually. However, the piezoelectric material volumes are found not active (100 mm 3 ) in all of the cases. This implies that using the maximum allowable piezoelectric material does not necessarily guarantee the higher power output, because of increased stiffness and decreased strain generation in the piezoelectric material. The harvester clearance from the ground ( ver d , Figure 12e) is found positive from all of the cases by satisfying Equation (35).
The power output and the normalized power density (NPD) is compared in Figure  12f. The NPD is defined as with the R 2 (coefficient of determination) of 0.999. Based on this formula, energy harvester designers can approximately estimate a total mass that will be needed to meet the power requirement. An additional optimization study at allowed M = 37.5 g found the power output 1.857 mW (the optimized design variables can be found in Table 4) and confirms this linear relation-only 0.6% different from a predicted value (1.845 mW).    Table 5 shows the NPDs and modified NPDs of various PE harvesters and their specification in terms of piezoelectric material, excitation amplitude, piezoelectric material volume, system volume, power output, and total mass. The modified NPD is a new index defined similarly as the NPD. However, the power output is divided by the volume of the system rather than that of the piezoelectric material. Therefore, the modified NPD can be referred to when the compactness is more important than the cost. All four optimized ΓEHs in this study show higher NPDs (15.49 × 10 3~2 3.10 × 10 3 kg·s·m −3 ) than those of any other PE harvesters reviewed in this table (0.016 × 10 3~1 1.55 × 10 3 kg·s·m −3 ). It is noted that the highest NPD among the referenced works, i.e., 11.55 × 10 3 kg·s·m −3 by Yang, et al. [63], was enabled by its large mass (100 g) and the high piezoelectric coupling constant of 275 pC/N. However, more than doubled NPD can be obtained from our design using a smaller mass ( allowed M = 60 g) and a smaller piezoelectric coupling constant (210 pC/N). In terms of the modified NPD, the proposed designs with allowed M = 45 g and 60 g have higher values than the design in [63]. This shows the effectiveness of the optimized ΓEH. By the way, the NPD increases as the mass becomes larger, but shows a saturated trend (around 24 × 10 3 kg·s·m −3 ), because the constraints of vertical dimension and stress (Equations (34) and (36)) become active.

Conclusions
This paper performed design optimization of a new ΓEH using the GCRM-P and DE algorithm. Some conclusive remarks are as follows: (1) The accuracy of the proposed GCRM-P model used for the frequency response analysis of the ΓEH was experimentally validated with the error 5.5% for the peak power frequency. (2) The proposed DE-based approach successfully provided the optimized solutions with the high NPDs, while satisfying the six design constraints. Specifically, we could obtain higher harvester NPDs than the multiple mesoscale PE harvesters from recent studies. (3) The linear relation between the harvester mass (Mallowed) and power performance does not necessarily mean that all the design variables are linearly scaled-they need to be carefully chosen to maximize the power output while satisfying all of the constraints, especially the stress and the natural frequency measures.
This study focused on the power and NPD at a fixed frequency as a performance index. A more practical design study considering broadband characteristics, durability, manufacturing cost, and compactness of the ΓEH remains as future work.