Topology Optimization of Piezoelectric Energy Harvesters for Enhanced Open-Circuit Voltage Subjected to Harmonic Excitations

Energy harvesting devices made of piezoelectric material are highly anticipated energy sources for power wireless sensors. Tremendous efforts have been made to improve the performance of piezoelectric energy harvesters (PEHs). Noticeably, topology optimization has shown an attractive potential to design PEHs with enhanced energy conversion efficiency. In this work, an alternative yet more practical design objective was considered, where the open-circuit voltage of PEHs is enhanced by topologically optimizing the through-thickness piezoelectric material distribution of plate-type PEHs subjected to harmonic excitations. Compared to the conventional efficiency-enhanced designs, the open-circuit voltage of PEHs can be evidently enhanced by the proposed method while with negligible sacrifice on the energy conversion efficiency. Numerical investigations show that the voltage cancellation effect due to inconsistent voltage phases can be effectively ameliorated by optimally distributed piezoelectric materials.


Introduction
Energy harvesting devices made of piezoelectric material are highly anticipated energy sources for power wireless sensors [1,2]. Tremendous efforts have been made to improve the performance of piezoelectric energy harvesters (PEHs). According to the literature, PEHs were mainly designed for energy conversion efficiency [3,4], output power [5,6] and bandwidth [7,8]. From a practical point of view, the electric power generated by PEHs is of the most interest. Therefore, numerous researches focused on designing PEHs to achieve higher output power with broader operating bandwidth. To this end, various shapes were proposed for the piezoelectric patch layer of cantilever-type PEHs, such as rectangles [9], triangles [10] and trapezoids [11]. It was shown that the voltage frequency response and harvested power of PEHs could be enhanced by modifying the in-plane geometry or introducing through-thickness geometrical features within the piezoelectric layer, such as cross-section [12], cavities [13] and gaps [14] of varying sizes. Though the voltage or output power of PEHs can be increased to a certain extent by these designs, the overall efficiency is maintained at a low level for the confined design space.
Noticeably, topology optimization [15] has shown an attractive potential to design PEHs with enhanced energy conversion efficiency [16,17]. To our best knowledge, Zheng et al. [18] are the first to introduce topology optimization for PEHs, considering the energy conversion efficiency as the design objective. Chen et al. [19] proposed to topologically design PEHs under harmonic loads with the adoption of level-set methods [20,21]. Kim et al. [22] demonstrated that the intrinsic and objective-dependent conditions should be satisfied simultaneously to achieve stable convergence of the topological evolution of PEHs when considering three penalty exponents. Vatanabe et al. [23] proposed to design functionally graded piezocomposite materials using topology optimization for enhanced electromechanical coupling coefficient. Almeida and Pavanello [24] topologically designed the through-thickness geometrical features of a bimorph harvester using the bi-directional evolutionary structural optimization method [25]. He et al. [26] proposed a multi-material topology optimization method to enhance the energy conversion efficiency of PEHs with simultaneous distribution of the elastic, piezoelectric and void materials. Previous investigations showed that the amount of piezoelectric material embedded in structures is positively correlated to the energy conversion efficiency [22,27]. The PEH energy conversion efficiency enhancement by means of topology optimization was comprehensively reviewed recently by Amlashi et al. [28].
It should be noted that enhancing the energy conversion efficiency of PEHs does not necessarily result in a high output voltage or power, which is, however, of more importance in practical uses of miniature sensors such as ambient light, temperature and pressure sensors. It was shown that the output voltage or power can be enhanced through the designs of geometry, tip mass position, piezoelectric material layout and electrodes [29]. Rupp et al. [30] were the first to pave the way for topology optimization of PEHs with respect to the power considering harmonic excitations. It was shown that in order to enhance the output power, the piezoelectric patch needs to be topologically optimized to alter the structural modes such that the structure is simultaneously tuned to the driving frequency and prevents charge cancellation. Following this framework, Lee and Youn [31] designed a piezoelectric skin patched to an outdoor condensing unit to harvest vibration energy and experimentally validated the performance; Noh and Yoon [27] demonstrated the effects of penalization of the piezoelectric material interpolation model in topology optimization on the harvesting efficiency. Notice that the aforementioned designs for output power mainly focused on the in-plane material distribution of the piezoelectric patch, while the designs of the through-thickness geometrical features of plate-type PEHs were limited except for the works reported by Vatanabe et al. [32] and Wein et al. [33].
This work followed the same framework developed by Rupp et al. [30] and focused on the design of the through-thickness piezoelectric material distribution of plate-type PEHs subjected to harmonic excitations. The open-circuit voltage was considered as the design objective for simplicity consideration, which is equivalent to a power design assuming an infinitely large resistor connected. Meanwhile, the energy conversion efficiency was also considered for design. The results were compared to those obtained from opencircuit voltage designs. The innovation of this paper was to achieve the performance optimization of piezoelectric materials along the thickness direction at both low and high frequencies. The remainder of this paper is organized as follows. Section 2 briefly reviews the finite element formulations of piezoelectric structures. Section 3 establishes the topology optimization model and considers, respectively, the open-circuit voltage and the energy conversion efficiency as the design objective. The analytical sensitivities derived from the two design objectives were also demonstrated. Section 4 validates the method through a series of designs of a benchmark clamped-clamped piezoelectric plate. The conclusion is drawn in Section 5.

Finite Element Analysis of Piezoelectric Structures
The piezoelectric material can convert mechanical into electrical energy and vice versa. The linearly coupled mechanical and electrical constitutive behavior of piezoelectric material can be written in the following stress-charge form as [34]: where σ and ε are the mechanical stress and strain tensors, respectively, and D is the electric displacement vector. C E and κ S are the elasticity and dielectric permittivity tensors at constant electric field E and constant mechanical strain ε, correspondingly. e is the piezoelectric coefficient tensor, and the superscript T denotes matrix transposition. The generalized matrix dimension of C E and κ S are 6 by 6 and 3 by 3 symmetric matrices, respectively, and e is a 6 by 3 matrix. As the plane strain assumption is adopted, the dimensions of the matrices C E and κ S are thus 3 by 3 and 2 by 2, respectively, and the e is a 3 by 2 matrix. Due to the brittleness of piezoceramic, the piezoelectric energy harvesting structure consists of a piezoelectric domain and an elastic (aluminum or copper) support layer. Figure 1 shows the schematic representation of the modeled piezoelectric energy harvesting structure under the plane strain assumption. The solid arrows with point represent a dynamic pressure load acting on the bottom surface. The direction of PZT material is polarized along the z(3) direction, which is indicated by the solid double arrow.The mechanical properties (stiffness) of these physical electrodes can be neglected, while for the electric property, the top and bottom electrodes represent two equipotential surfaces in the model. where and are the mechanical stress and strain tensors, respectively, and is th electric displacement vector. and are the elasticity and dielectric permittivity ten sors at constant electric field and constant mechanical strain , correspondingly. i the piezoelectric coefficient tensor, and the superscript T denotes matrix transposition The generalized matrix dimension of and are 6 by 6 and 3 by 3 symmetric matrice respectively, and is a 6 by 3 matrix. As the plane strain assumption is adopted, th dimensions of the matrices and are thus 3 by 3 and 2 by 2, respectively, and the is a 3 by 2 matrix.
Due to the brittleness of piezoceramic, the piezoelectric energy harvesting structur consists of a piezoelectric domain and an elastic (aluminum or copper) support layer. Fig  ure 1 shows the schematic representation of the modeled piezoelectric energy harvesting structure under the plane strain assumption. The solid arrows with point represent a dy namic pressure load acting on the bottom surface. The direction of PZT material is polar ized along the z(3) direction , which is indicated by the solid double arrow.The mechanica properties (stiffness) of these physical electrodes can be neglected, while for the electri property, the top and bottom electrodes represent two equipotential surfaces in the mode As shown in Figure 1, the design domain is assumed to be discretized into N squar finite elements, and each element is associated with the continuously defined topology variable , , denoting the existence of piezoelectric material: 0 , 1 with e = 1, 2, …, N.
(2 Following the piezoelectric material with penalization (PEMAP) model [35], which is an extension of the well-known solid isotropic material with penalization (SIMP) mode [15], the topology variable is linked to material constitutive parameters in the following form: As shown in Figure 1, the design domain is assumed to be discretized into N square finite elements, and each element is associated with the continuously defined topology variable ρ pzt,e , denoting the existence of piezoelectric material: 0 ≤ ρ pzt,e ≤ 1 with e = 1, 2, . . . , N. ( Following the piezoelectric material with penalization (PEMAP) model [35], which is an extension of the well-known solid isotropic material with penalization (SIMP) model [15], the topology variable is linked to material constitutive parameters in the following form: pzt,e κ pzt (6) where {m pzt ,C pzt } and {m nonpzt ,C nonpzt } are, respectively, the mass and elastic tensors of piezoelectric and non-piezoelectric materials; e pzt and κ pzt are the piezoelectric coupling and dielectric coefficients, respectively; p uu , p uφ and p φφ are penalization coefficients for the elastic, piezoelectric coupling and dielectric properties, respectively. According to the above material interpolation model, the piezoelectric material is obtained with ρ pzt,e = 1; elastic material (aluminum) is obtained with ρ pzt,e = 0. Note that intermediate values (0 < ρ pzt,e < 1) represent composite material elements without physical meanings. These artificial densities should be penalized in a manner analogous to other continuous optimization approximations of 0-1 problems [15]. Details about the setting of the penalization coefficients are provided in Appendix A. By using Hamilton's variational principle and neglecting the damping effect, the detailed finite element derivation and calculation of the piezoelectric element matrix are introduced in reference [34]. It can be seen in Figure 1 that each square element has four nodes, and each node has three degrees of freedom (DOFs) (two mechanical DOFs regarding the displacements in x and z directions and one electric potential DOF). The finite element formulation of piezoelectric structures subjected to harmonic excitations can be written as follows: where F t and Q t are the time-harmonic excitations of mechanical force and electrical charge and Φ(n × 1) are the amplitudes of harmonic force, electrical charge, displacement and potential, respectively.
are, respectively, the global mass, stiffness, piezoelectric coupling and dielectric matrices. The external charge Q is set to zero for energy harvesting problems. The above harmonic analysis formulation can be further derived to the following simplified form [34]: where ω is the imposed excitation frequencies ( f = ω/2π). The piezoelectric plate for energy harvesting is commonly sandwiched between two silver electrode layers. In practice, the electrode layers are assumed to be zero, and their effect is equivalently modeled by imposing equipotential boundary conditions [28]. The potentials of the nodes on the electrode layers are gathered from the global potential vector by using a Boolean matrix, and their values are enforced to be the same. This study focused on the application of topology optimization to optimize PEHs by developing an in-house FE code under dynamic loads.

Optimization Model and Sensitivity Analysis
The design objectives were considered, respectively, in this work, namely, the energy conversion efficiency and the open-circuit voltage. The energy conversion efficiency is defined as the ratio of the generated electric energy to the external mechanical work [26]: Here, Π S = (1/4)U T K uu U and Π E = (1/4)Φ T K φφ Φ are the elastic and dielectric energies under harmonic analysis. Note that harmonic stress and strain components are generally not in phase with each other, so the two energies are estimated as the real parts of their cycle averaged values.
The open-circuit voltage is defined as the electric potential difference between the top and bottom electrode layers of the piezoelectric patch, as shown in Figure 1. Since equipo-tential boundary conditions are imposed at the two layers, the potential of an arbitrarily prescribed node on the top layer can be adopted for denoting the open-circuit voltage: where Φ oc is the global electric potential vector, L dummy is a selecting vector in which the component corresponding to the prescribed node is set to one while the others are zero.
In the following analysis, the open-circuit voltage is measured without connecting to an external harvesting circuit, i.e., an infinite external resistance is assumed. By considering the above two properties as design objectives, the mathematical model for topology optimization of piezoelectric materials subjected to harmonic excitations is formulated as follows: where V ρ pzt and V 0 are, respectively, the volumes of piezoelectric material and design domain, and V is the prescribed volume fraction.
In order to favor gradient-based mathematical programming algorithms [36,37], sensitivities of the two design objectives with respect to topology variables are derived using the adjoint method. The sensitivity of energy conversion efficiency w.r.t. ρ pzt,e is derived as: The derivatives of the elastic and dielectric energies w.r.t. ρ pzt,e can be obtained as [26]: where λ 1 , µ 1 , λ 2 and µ 2 are solutions of the following two adjoint problems: The derivatives of the constitutive matrices w.r.t. ρ pzt in Equation (13) can be straightforwardly derived according to the defined PEMAP material interpolation model.
In order to simplify the mathematical minimization problem, the design objective of open-circuit voltage can be equivalently reformulated in the following form: Similarly, the sensitivity of open-circuit voltage w.r.t. ρ pzt,e based on the adjoint variable method is derived as: Considering that (∂F/∂ρ pzt,e ) = 0 and Q = 0, the sensitivity equation can be further reduced by algebraic operations: In which the derivatives of the two-state variables {∂U/∂ρ pzt,e , ∂Φ/∂ρ pzt,e } can be eliminated with the solutions of the following adjoint problem: (18) Additionally, the final sensitivity estimation form is derived as: To subsequently calculate Equation (19), all the derivative terms can be obtained from the interpolation scheme in Equations (5)-(8). Herein, we did not develop it.

Numerical Examples
As is shown in Figure 2a In order to simulate the effect of the two electrode layers in Figure 2a, equipotential boundary conditions were imposed, and their influence on the energy harvesting performance was demonstrated in Figure 3 and Table 1. It can be observed that the potential distribution within the piezoelectric layer was changed as a result of the equipotential boundary conditions. Notice that a piezoelectric patch without electrodes is practically useless as the voltage varies along the length (x) direction, and thus there exists no specific output open-circuit voltage or power. Therefore, the attachment of electrode layers is necessary, and it inevitably brings the issue of charge or voltage cancellation [38,39], resulting in the sacrifice of the energy harvesting performance. In the following, the PZT-5A material distribution was topologically optimized for energy conversion efficiency and open-circuit voltage in Sections 4.1 and 4.2, respectively. The filter radius was set as two times the element length. The volume fraction of PZT-5A was constrained below 30%. The method of moving asymptotes (MMA, [40]) was adopted as the optimizer. Additionally, the post-processing scheme with a threshold was adopted to eliminate intermediate densities [41].  In order to simulate the effect of the two electrode layers in Figure 2a, equipotential boundary conditions were imposed, and their influence on the energy harvesting performance was demonstrated in Figure 3 and Table 1. It can be observed that the potential distribution within the piezoelectric layer was changed as a result of the equipotential boundary conditions. Notice that a piezoelectric patch without electrodes is practically useless as the voltage varies along the length (x) direction, and thus there exists no specific output open-circuit voltage or power. Therefore, the attachment of electrode layers is necessary, and it inevitably brings the issue of charge or voltage cancellation [38,39], resulting in the sacrifice of the energy harvesting performance. In the following, the PZT-5A material distribution was topologically optimized for energy conversion efficiency and open-circuit voltage in Sections 4.1 and 4.2, respectively. The filter radius was set as two times the element length. The volume fraction of PZT-5A was constrained below 30%. The method of moving asymptotes (MMA, [40]) was adopted as the optimizer. Additionally, the post-processing scheme with a threshold was adopted to eliminate intermediate densities [41]. In order to simulate the effect of the two electrode layers in Figure 2a, equipotential boundary conditions were imposed, and their influence on the energy harvesting performance was demonstrated in Figure 3 and Table 1. It can be observed that the potential distribution within the piezoelectric layer was changed as a result of the equipotential boundary conditions. Notice that a piezoelectric patch without electrodes is practically useless as the voltage varies along the length (x) direction, and thus there exists no specific output open-circuit voltage or power. Therefore, the attachment of electrode layers is necessary, and it inevitably brings the issue of charge or voltage cancellation [38,39], resulting in the sacrifice of the energy harvesting performance. In the following, the PZT-5A material distribution was topologically optimized for energy conversion efficiency and open-circuit voltage in Sections 4.1 and 4.2, respectively. The filter radius was set as two times the element length. The volume fraction of PZT-5A was constrained below 30%. The method of moving asymptotes (MMA, [40]) was adopted as the optimizer. Additionally, the post-processing scheme with a threshold was adopted to eliminate intermediate densities [41].

Design for Enhanced Energy Conversion Efficiency
The energy conversion efficiency enhanced designs obtained by topology optimization for variant excitation frequencies are shown in Figure 4 and Table 2. For each of the excitation frequencies, the designed PZT-5A material distribution enhances an average 10 times the efficiency than that of the original design with top-layered PZT-5A material. Meanwhile, it can also be observed that the designed PZT-5A material distribution is highly dependent on the excitation frequency. Generally speaking, PZT-5A material is preferentially distributed in highly strained regions so as to generate more dielectric energy. This explains why PZT-5A materials are distributed around the two clamped ends for low-frequency excitations while they are distributed towards the center for high-frequency excitations. Frequency responses of the energy conversion efficiencies of the enhanced designs are shown in Figure 5. The frequency responses of the designs of Figure 4a,b are almost identical, which can actually be anticipated from their similar PZT-5A distribution. The designs of Figure 4c,d evidently outperforms the designs of Figure 4a,b in terms of the efficiency for high excitation frequencies, while the opposite for low excitation frequencies.

Design for Enhanced Energy Conversion Efficiency
The energy conversion efficiency enhanced designs obtained by topol optimization for variant excitation frequencies are shown in Figure 4 and Table 2. For e of the excitation frequencies, the designed PZT-5A material distribution enhances average 10 times the efficiency than that of the original design with top-layered PZT material. Meanwhile, it can also be observed that the designed PZT-5A mate distribution is highly dependent on the excitation frequency. Generally speaking, PZT material is preferentially distributed in highly strained regions so as to generate m dielectric energy. This explains why PZT-5A materials are distributed around the clamped ends for low-frequency excitations while they are distributed towards the cen for high-frequency excitations. Frequency responses of the energy conversion efficien of the enhanced designs are shown in Figure 5. The frequency responses of the design Figure 4a,b are almost identical, which can actually be anticipated from their similar P 5A distribution. The designs of Figure 4c,d evidently outperforms the designs of Fig  4a,b in terms of the efficiency for high excitation frequencies, while the opposite for excitation frequencies.      The iterative history and topological evolution at 4 kHz and 16 kHz dur optimization steps are shown in Figure 6. As the iteration progresses, the top configuration becomes clearer. Meanwhile, the energy conversion efficiency gr increases and finally tends to be stable. For further demonstration, the Von Mises stress and voltage contours efficiency-enhanced designs of Figure 4a,d are shown in Figure 7. From the v mode shape of the optimized designs in Figure 7a,b, it can be seen that the material is preferentially distributed at high stressed regions for efficiency enhan purposes. On the contrary, the voltage contours in Figure 7c,d show that posit negative voltages are canceled along the length direction of the top electrod The iterative history and topological evolution at 4 kHz and 16 kHz during the optimization steps are shown in Figure 6. As the iteration progresses, the topological configuration becomes clearer. Meanwhile, the energy conversion efficiency gradually increases and finally tends to be stable. The iterative history and topological evolution at 4 kHz and 16 kHz during the optimization steps are shown in Figure 6. As the iteration progresses, the topological configuration becomes clearer. Meanwhile, the energy conversion efficiency gradually increases and finally tends to be stable. For further demonstration, the Von Mises stress and voltage contours of the efficiency-enhanced designs of Figure 4a,d are shown in Figure 7. From the vibration mode shape of the optimized designs in Figure 7a,b, it can be seen that the PZT-5A material is preferentially distributed at high stressed regions for efficiency enhancement purposes. On the contrary, the voltage contours in Figure 7c,d show that positive and negative voltages are canceled along the length direction of the top electrode layer, resulting in open-circuit voltages of low magnitude. In order to be consistent with the settings of the original design in Figure 2, equipotential boundary conditions are imposed on the top layer to simulate the electrode layer, and the bottom layer is grounded, as shown in Figure 7c,d. The same electrical boundary conditions are imposed in the For further demonstration, the Von Mises stress and voltage contours of the efficiencyenhanced designs of Figure 4a,d are shown in Figure 7. From the vibration mode shape of the optimized designs in Figure 7a,b, it can be seen that the PZT-5A material is preferentially distributed at high stressed regions for efficiency enhancement purposes. On the contrary, the voltage contours in Figure 7c

Design for Enhanced Open-Circuit Voltage
Topology optimization is further conducted considering the open-circuit voltage the design objective to maximize, and the results are shown in Figure 8. Unlike previous efficiency enhanced designs, open-circuit voltage enhanced designs do exhibit an obvious evolution tendency of PZT-5A material distribution on the excitat frequency. Table 3 shows the optimization results of the open-circuit voltage enhan designs at variant excitation frequencies. It was noticed that in order to enhance the op circuit voltage, PZT-5A material is distributed such that the first or second structu eigenfrequency is tuned to approach the excitation frequency for this design ca Convergence challenges encountered near the resonance for the excitation freque ranged from 6.5 kHz to 9.5 kHz because the material distribution considerably chan the structural eigenfrequencies and vice versa.

Design for Enhanced Open-Circuit Voltage
Topology optimization is further conducted considering the open-circuit voltage as the design objective to maximize, and the results are shown in Figure 8. Unlike the previous efficiency enhanced designs, open-circuit voltage enhanced designs do not exhibit an obvious evolution tendency of PZT-5A material distribution on the excitation frequency. Table 3 shows the optimization results of the open-circuit voltage enhanced designs at variant excitation frequencies. It was noticed that in order to enhance the open-circuit voltage, PZT-5A material is distributed such that the first or second structural eigenfrequency is tuned to approach the excitation frequency for this design case. Convergence challenges encountered near the resonance for the excitation frequency ranged from 6.5 kHz to 9.5 kHz because the material distribution considerably changes the structural eigenfrequencies and vice versa.

Design for Enhanced Open-Circuit Voltage
Topology optimization is further conducted considering the open-circuit voltage the design objective to maximize, and the results are shown in Figure 8. Unlike previous efficiency enhanced designs, open-circuit voltage enhanced designs do exhibit an obvious evolution tendency of PZT-5A material distribution on the excitat frequency. Table 3 shows the optimization results of the open-circuit voltage enhan designs at variant excitation frequencies. It was noticed that in order to enhance the op circuit voltage, PZT-5A material is distributed such that the first or second structu eigenfrequency is tuned to approach the excitation frequency for this design ca Convergence challenges encountered near the resonance for the excitation freque ranged from 6.5 kHz to 9.5 kHz because the material distribution considerably chan the structural eigenfrequencies and vice versa.       Figure 10a. It can also be seen that some fluctuations appear in the conv gence curve because a part of the piezoelectric materials distributes near the top ed which can produce voltage oscillations. By the post-processing scheme with thresho the final open-circuit voltage is 18.78 V, and the energy conversion efficiency is 12.0 Figure 10b shows that the piezoelectric material tends to distribute towards the cente the structure at an excitation frequency of 16 kHz.   Figure 10a. It can also be seen that some fluctuations appear in the convergence curve because a part of the piezoelectric materials distributes near the top edge, which can produce voltage oscillations. By the post-processing scheme with threshold, the final open-circuit voltage is 18.78 V, and the energy conversion efficiency is 12.02%. Figure 10b shows that the piezoelectric material tends to distribute towards the center of the structure at an excitation frequency of 16 kHz.
The frequency response results of the open-circuit voltage and energy conversion efficiency for the designs obtained under excitation frequencies of 4 kHz and 16 kHz are compared in Figure 11. The first two eigenfrequencies of the original design in Figure 2 are 7088 Hz and 16,121 Hz. As PZT-5A material has a comparable modulus to the substrate, the newly designed harvesters dynamically behave similarly to the original design with adjacent eigenfrequencies. Abrupt jumps in the frequency response curves can be observed at the resonance frequencies. For both excitation frequencies (4 kHz and 16 kHz), voltageenhanced designs outperform the efficiency-enhanced designs in terms of the open-circuit voltages over almost the complete frequency range. Meanwhile, the comparison results show also that the enhanced designs for specific operation frequencies can, in fact, benefit the energy harvesting performance over a much broader working bandwidth. iterations in Figure 10a. It can also be seen that some fluctuations appear in the convergence curve because a part of the piezoelectric materials distributes near the top edge, which can produce voltage oscillations. By the post-processing scheme with threshold, the final open-circuit voltage is 18.78 V, and the energy conversion efficiency is 12.02%. Figure 10b shows that the piezoelectric material tends to distribute towards the center of the structure at an excitation frequency of 16 kHz.  Figure 11. The first two eigenfrequencies of the original design in Figure 2 are 7088 Hz and 16121 Hz. As PZT-5A material has a comparable modulus to the substrate, the newly designed harvesters dynamically behave similarly to the original

Conclusions
This work proposed a topology optimization approach to achieve piezoelectric energy harvesters with enhanced energy conversion efficiency or open-circuit voltage. The piezoelectric material layout along the thickness direction was topologically optimized for the two considered design objectives. It can be concluded from the numerical results that: (i) for realizing a high energy conversion efficiency, the piezoelectric materials should be placed at higher stress/strain regions to generate more electricity energy; (ii) for realizing a high open-circuit voltage, the distribution of piezoelectric materials should be tuned to match the external excitation frequency to obtain larger electric power under constant electric resistance. Meanwhile, it can be inferred that the voltage cancellation issue can be evidently ameliorated by means of open-circuit enhancement designs. In the future, the design of PEHs circuits with varying resistance needs to be integrated into the piezoelectric dynamic system. Furthermore, the polarization directions and the manufacturability of piezoelectric materials will be considered to achieve high-efficient PEHs with a broader operational bandwidth.

Conclusions
This work proposed a topology optimization approach to achieve piezoelectric energy harvesters with enhanced energy conversion efficiency or open-circuit voltage. The piezoelectric material layout along the thickness direction was topologically optimized for the two considered design objectives. It can be concluded from the numerical results that: (i) for realizing a high energy conversion efficiency, the piezoelectric materials should be placed at higher stress/strain regions to generate more electricity energy; (ii) for realizing a high open-circuit voltage, the distribution of piezoelectric materials should be tuned to match the external excitation frequency to obtain larger electric power under constant electric resistance. Meanwhile, it can be inferred that the voltage cancellation issue can be evidently ameliorated by means of open-circuit enhancement designs. In the future, the design of PEHs circuits with varying resistance needs to be integrated into the piezoelectric dynamic system. Furthermore, the polarization directions and the manufacturability of piezoelectric materials will be considered to achieve high-efficient PEHs with a broader operational bandwidth.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.