Nonlinear finite element system simulation of piezoelectric vibration-based energy harvesters

Abstract

Piezoelectric vibration-based energy harvesters consist of an electromechanical structure and an electric circuitry, influencing each other. We propose a novel approach that allows a finite element based system simulation of nonlinear electromechanical structures coupled to nonlinear electric circuitries. In the finite element simulation the influence of the electric circuit on the electromechanical structure is considered via the vector of external forces, using an implicit time integration scheme. To demonstrate the applicability of the new simulation method an active power circuit is considered. Several examples of piezoelectric vibration-based energy harvesters, connected to standard or synchronized switch harvesting on inductor (SSHI) circuits, showing linear or nonlinear mechanical behavior, are studied to validate the proposed simulation method against numerical results reported in the literature. The advocated method allows for consistent and efficient simulations of complete nonlinear energy harvesters using only one software tool.

1. Introduction

Energy harvesting is of great importance in the digital world to enable self-powered systems like wireless sensors or other small electronic components. Vibration-based energy harvesting converts ambient vibration energy to electrical energy. Compared to other energy conversion principles like electrostatic (Mitcheson et al., 2004) or electromagnetic (Glynne-Jones et al., 2004; Williams and Yates, 1996) transduction mechanisms the advantage of the piezoelectric effect (Roundy et al., 2004; Rupitsch, 2019) is the high power density and the ease of application of piezoelectric materials (Erturk and Inman, 2011).

A piezoelectric vibration-based energy harvester (PVEH) is composed of an electromechanical structure and an electric circuit to extract the energy. Various types of piezoelectric transducers are proposed in the literature, for example unimorph (a single piezoelectric layer bonded to a substrate layer) and bimorph (a single substrate layer with two symmetric piezoelectric layers bonded on each side) cantilever beams (Safaei et al., 2019). Beam-type harvesters are very efficient when excited at their resonance frequency but become less efficient when the excitation is below or above their resonance frequency. Numerical simulations are a valuable tool to develop PVEHs as they increase the understanding of the energy harvesting problem and allow to match the electromechanical structure and the electric circuit according to the excitation range.

The electromechanical structure and the electric circuit influence each other. Therefore, an accurate prediction of the behavior of a PVEH requires accurate modeling of both the electromechanical structure and the electric circuit.

Analytical solutions are possible for relatively simple geometries of the electromechanical structure in combination with simple electric circuits. Erturk and Inman (2008b) clarified some issues in the mathematical modeling of piezoelectric-based energy harvesting. They identified some problems of single degree of freedom models, mainly related to inaccuracies in the modeling of the electromechanical coupling. Furthermore, they suggested an analytical distributed parameter model based on the Euler-Bernoulli beam theory (Erturk and Inman, 2008a). These analytical models are valuable tools to increase understanding, to get fast solutions for simple PVEH and to validate numerical simulations (e.g. De Marqui et al., 2009; Rupp et al., 2009; Yang and Tang, 2009). To overcome the restrictions of the analytical methods, numerical methods as the Finite Element Method (FEM) are frequently applied to simulate more complex or arbitrarily shaped PVEHs.

While being computationally expensive, the FEM allows for optimization of the electromechanical structure (Noh and Yoon, 2012; Park et al., 2012; Wein et al., 2013), for modeling of nonlinear electro-elastic behavior (Vu et al., 2007), for detailed material parameter analysis (Daniels et al., 2013), for the consideration of multiple segmented plate structures with different electric circuit connection patterns (Lumentut and Shu, 2018, 2021) and, for example, for the simulation of crack propagation (Abdollahi and Arias, 2012). Various approaches are introduced in the literature how to couple the simulation of electric circuits with an FE analysis of electromechanical structures.

In Elvin and Elvin (2009) an explicit coupling is presented, which utilizes an available FEM package for a purely mechanical simulation and a standard circuit simulation software for the electric circuit. The coupling between the circuit simulation and the FE simulation is done via introducing equivalent piezoelectric loads in the mechanical model and applying equivalent electric voltages in the electric model. As already concluded in Elvin and Elvin (2009), the drawback of this method is that the explicit solution technique can be computationally expensive.

In Zhu et al. (2009) an electromechanical structure with an attached load resistance is solved with the commercial FE Software ANSYS. However, this system simulation is restricted to a combination of linear circuit elements. De Marqui et al. (2009) also considered only a load resistance as a passive electric circuit and developed an electromechanically coupled finite element plate model for predicting the electrical power output of piezoelectric energy harvester plates. Furthermore, commercial FE packages are utilized to validate analytical models (Schoeftner and Buchberger, 2013), to simulate electromechanical structures with simple electric circuits (Abdelkefi et al., 2014; Akbar and Curiel-Sosa, 2019) or to identify parameters for an analytical modeling approach (Xiong and Oyadiji, 2014).

In Wu and Shu (2015) an equivalent steady-state load impedance of a standard and an SSHI circuit is derived. The load impedances can be directly included in available FE software to allow for steady state simulations of PVEHs. Notably, this approach is extensible to multiple piezoelectric oscillators. Cheng et al. (2016) proposed two coupled models: The first one is an equivalent circuit model of the electromechanical structure. The equivalent circuit model parameters are estimated by FE analysis. The second one integrates the impedance of the circuit into a FE framework. Furthermore, a comparison of the two modeling techniques is carried out.

Recently, Gedeon and Rupitsch (2018) developed a system simulation method based on the FEM. The matrices of a FEM discretization are directly imported into an electric circuit simulation with Simscape (Matlab/Simulink). All capacitive and electromechanical coupling effects are included and a model order reduction technique is applied to the FEM simulation to reduce the computational costs. No nonlinearities of the electromechanical structure are considered in this work.

However, all FE based methods reported in the literature, which are not coupled to circuit simulation software, are limited to linear circuit elements and passive electric circuits. The FE based methods, which are combined with an external circuit simulation tool consider only linear electromechanical structures or the coupling between the electromechanical structure simulation and the electric circuit simulation is not very efficient.

To overcome the mentioned drawbacks of existing FE methods for PVEHs, a novel system simulation approach for nonlinear electric circuits coupled to nonlinear electromechanical structures using only the FEM is proposed. The applicability of the presented approach is demonstrated by simulating a unimorph PVEH attached to an SSHI electric circuit. Furthermore, a bimorph PVEH with nonlinear elastic behavior is considered in combination with nonlinear electric circuitry.

2. Governing equations of piezoelectricity

Within this contribution index notation in accordance to Institute of Electrical and Electronics Engineers (IEEE) (1988) is applied. The piezoelectric body $\mathrm{\Omega }$ is characterized by the help of the linear mechanical strain tensor $S_{\mathrm{ij}}$ and the electric field $E_i$

$$S_{\mathrm{ij}}=\frac{1}{2}[u_{i,j}+u_{j,i}]$$

(1)

$$E_i=-{\phi }_{,i}$$

(2)

Here, $u_i$ is the mechanical displacement and $\phi$ is the electric voltage. The mechanical equation is given by the balance of linear momentum as

$$T_{\mathrm{ij},i}=\rho {\stackrel{··}{u}}_j$$

(3)

Here, $\rho$ is the material density, $T_{\mathrm{ij}}$ are the components of the mechancial stress tensor and the double dot symbolizes the second derivative with respect to time. The electric equation is described by Gauss’ law considering that piezoelectric materials are insulating (no free volume charges) as

$$D_{i,i}=0$$

(4)

Here, $D_i$ is the dielectric displacement. The set of equations (1)–(4) has to be completed with a constitutive law that specifies the material behavior. Within this contribution both, a linear piezoelectric constitutive law and a nonlinear piezoelectric constitutive law, are applied. The linear piezoelectric constitutive law can be modeled with the help of the enthalpy density $H$ (IEEE, 1988)

$$H=\frac{1}{2}{c}_{\mathrm{ijkl}}^E{S}_{\mathrm{ij}}{S}_{\mathrm{kl}}-e_{\mathrm{kij}}{E}_k{S}_{\mathrm{ij}}-\frac{1}{2}{\epsilon }_{\mathrm{ij}}^S{E}_i{E}_j$$

(5)

Here, $c_{\mathrm{ijkl}}^E$ are the components of the elasticity tensor at constant electric field, $e_{\mathrm{kij}}$ is the piezoelectric constant tensor, and ${\epsilon }_{\mathrm{ij}}^S$ is the dielectric constant at constant strain. The following relations between $H$ , $T_{\mathrm{ij}}$ , $D_i$ , $S_{\mathrm{ij}}$ , and $E_i$ hold

$$T_{\mathrm{ij}}=\frac{\partial H}{\partial S_{\mathrm{ij}}}$$

(6)

$$D_i=-\frac{\partial H}{\partial E_i}$$

(7)

For details about the aforementioned relations please refer to the literature (IEEE, 1988; Rupitsch, 2019). The linear piezoelectric constitutive law resulting from equations (5)–(7) reads

$$T_{\mathrm{ij}}=c_{\mathrm{ijkl}}^E{S}_{\mathrm{kl}}-e_{\mathrm{kij}}{E}_k$$

(8)

$$D_i=e_{\mathrm{ikl}}{S}_{\mathrm{kl}}+{\epsilon }_{\mathrm{ij}}^S{E}_j$$

(9)

For small strains, which for PVEHs are usually related to small base accelerations, the linear piezoelectric model gives satisfactory results. However, in Stanton et al. (2012) the influences of nonlinearities in the elastic, electroelastic and damping behavior were analyzed. It was shown that a nonlinear elastic material model results in better agreement with experimental results. The model was introduced in a 1D setting and postulates a nonlinear polynomial function for the stress $T_{11}$

$$T_{11}=c_{1111}^E{S}_{11}-e_{311}{E}_3+c_4{S}_{11}^3+c_6{S}_{11}^5$$

(10)

with the additional parameters $c_4$ and $c_6$ . Here, the nonlinear elasticity model should be applied in 3D FE simulations and therefore the linear constitutive equations (5) and (8) are extended as follows

$$\begin{array}{c}\multicolumn{1}{c}{H=\frac{1}{2}{c}_{\mathrm{ijkl}}^E{S}_{\mathrm{ij}}{S}_{\mathrm{kl}}+\frac{1}{4}{c}_4{\left[S_{\mathrm{mn}}{A}_{\mathrm{mn}}\right]}^4+\frac{1}{6}{c}_6{\left[S_{\mathrm{mn}}{A}_{\mathrm{mn}}\right]}^6}\\ -e_{\mathrm{kij}}{E}_k{S}_{\mathrm{ij}}-\frac{1}{2}{\epsilon }_{\mathrm{ij}}^S{E}_i{E}_j\hfill \end{array}$$

(11)

$$\begin{array}{c}\multicolumn{1}{c}{T_{\mathrm{ij}}=c_{\mathrm{ijkl}}^E{S}_{\mathrm{kl}}-e_{\mathrm{kij}}{E}_k+c_4{\left[S_{\mathrm{mn}}{A}_{\mathrm{mn}}\right]}^3{A}_{\mathrm{ij}}}\\ +c_6{\left[S_{\mathrm{mn}}{A}_{\mathrm{mn}}\right]}^5{A}_{\mathrm{ij}}\hfill \end{array}$$

(12)

with the structural tensor $A_{\mathrm{ij}}=a_i{a}_j$ and the vector $a_i=e_{1i}$ being equal to the Cartesian basis vector and denoting the direction with the nonlinear elastic behavior. A nonlinear elastic behavior is thus only introduced in 1-direction which is a reasonable approach since this is the direction of the largest stresses and strains. Consistently with Stanton et al. (2012), the fifth order nonlinear term in equation (12) becomes later important for the bimorph PVEH, because a relatively heavy tip mass and a high level of base excitation is applied. The material behavior with respect to the electric field is still linear and the linear relation (9) for the dielectric displacement is also maintained.

Equations (1)–(4) along with the linear constitutive law (equations (8) and (9)) or the nonlinear constitutive law (equations (9) and (12)) can be solved with appropriate boundary conditions (Lerch, 1990)

$$u_i={\overline{u}}_i\mathrm{on}\partial {\mathrm{\Omega }}_{\mathrm{Du}}$$

(13)

$$\phi =\overline{\phi }\mathrm{on}\partial {\mathrm{\Omega }}_{D\phi }$$

(14)

$$T_{\mathrm{ij}}{n}_j={\overline{t}}_i\mathrm{on}\partial {\mathrm{\Omega }}_{\mathrm{Nu}}$$

(15)

$$D_i{n}_i=-\overline{Q}\mathrm{on}\partial {\mathrm{\Omega }}_{N\phi }$$

(16)

Here, the prescribed surface traction ${\overline{t}}_i$ and the free surface charge density $\overline{Q}$ are introduced. The boundary $\partial \mathrm{\Omega }$ of $\mathrm{\Omega }$ consists of subsets that do not overlap, such that $\partial {\mathrm{\Omega }}_D\cup \partial {\mathrm{\Omega }}_N=\partial \mathrm{\Omega }$ and $\partial {\mathrm{\Omega }}_D\cap \partial {\mathrm{\Omega }}_N=Ø$ .

3. Discretization of the equations

A fundamental step of the FEM is to transform the partial differential equations from their strong formulation (equations (3) and (4) along with equations (1) and (2), a constitutive law and boundary conditions) into their weak formulation. To obtain the desired weak formulation, equation (3) is multiplied by the vector-valued test function ${\eta }_j$ , vanishing on the Dirichlet boundary $\partial {\mathrm{\Omega }}_{\mathrm{Du}}$ , and integrated over $\mathrm{\Omega }$ . Similarly, equation (4) is multiplied with a scalar valued test function $\xi$ that becomes zero on the Dirichlet boundary $\partial {\mathrm{\Omega }}_{D\phi }$ , and integrated over the domain. Applying integration by parts and the divergence theorem and introducing the boundary conditions lead to the resulting weak forms

$${\int }_{\mathrm{\Omega }}\rho {\eta }_j{\stackrel{··}{u}}_j\mathrm{dV}+{\int }_{\mathrm{\Omega }}{\eta }_{j,i}{T}_{\mathrm{ij}}\mathrm{dV}={\int }_{\partial {\mathrm{\Omega }}_{\mathrm{Nu}}}{\eta }_j{\overline{t}}_j\mathrm{dA}$$

(17)

$${\int }_{\mathrm{\Omega }}{\xi }_{,i}{D}_i\mathrm{dV}=-{\int }_{\partial {\mathrm{\Omega }}_{N\phi }}\xi \overline{Q}\mathrm{dA}$$

(18)

In the FEM, the domain $\mathrm{\Omega }$ is subdivided into small discrete elements ${\mathrm{\Omega }}^e$ (“finite elements”). The unknown solutions, the displacement field $u_j$ , and the electric voltage $\phi$ , are approximated elementwise by means of polynomial ansatz functions and nodal values as $u_j^e=N_j^I{u}^I$ and ${\phi }^e={\hat{N}}^I{\phi }^I$ , whereby $N_j^I$ denotes the vector valued ansatz function (for the displacements) of degree of freedom I and ${\hat{N}}^I$ denotes the scalar valued ansatz function (for the electric potential) of degree of freedom I. The test functions ${\eta }_i$ and $\xi$ are approximated by means of the same ansatz functions. Introducing these approximations into the weak forms (17) and (18) result in two coupled vector-valued equations for the unknown nodal displacements and electric potential values

$$\underset{e=1}{\stackrel{n_{\mathrm{el}}}{\mathcal{A}}}[\underset{f^{\mathrm{dyn},u,I}}{\underbrace{{\int }_{{\mathrm{\Omega }}^e}{N}_j^I\rho {\stackrel{··}{u}}_j\mathrm{dV}}}+\underset{f^{\mathrm{int},u,I}}{\underbrace{{\int }_{{\mathrm{\Omega }}^e}{N}_{j,i}^I{T}_{\mathrm{ij}}\mathrm{dV}}}-\underset{f^{\mathrm{ext},u,I}}{\underbrace{{\int }_{\partial {\mathrm{\Omega }}_{\mathrm{Nu}}^e}{N}_i^I{\overline{t}}_i\mathrm{dA}}}]=0$$

(19)

$$\underset{e=1}{\stackrel{n_{\mathrm{el}}}{\mathcal{A}}}[\underset{f^{\mathrm{int},\phi ,I}}{\underbrace{{\int }_{{\mathrm{\Omega }}^e}{\hat{N}}_{,j}^I{D}_j\mathrm{dV}}}+\underset{f^{\mathrm{ext},\phi ,I}}{\underbrace{{\int }_{\partial {\mathrm{\Omega }}_{N\phi }^e}{\hat{N}}^I\overline{Q}\mathrm{dA}}}]=0$$

(20)

whereby $\underset{e=1}{\stackrel{n_{\mathrm{el}}}{\mathcal{A}}}$ denotes the assembly of all element contributions to the global vector-valued equations

$${\mathit{f}}^{\mathrm{dyn},u}+{\mathit{f}}^{\mathrm{int},u}-{\mathit{f}}^{\mathrm{ext},u}=0$$

(21)

$${\mathit{f}}^{\mathrm{int},\phi }+{\mathit{f}}^{\mathrm{ext},\phi }=0$$

(22)

Here, ${\mathit{f}}^{\mathrm{dyn},u}$ is the mechanical inertial force vector, ${\mathit{f}}^{\mathrm{int},u}$ is the mechanical vector of internal forces, ${\mathit{f}}^{\mathrm{ext},u}$ is the mechanical vector of external forces, ${\mathit{f}}^{\mathrm{int},\phi }$ is the electrical vector of internal forces, and ${\mathit{f}}^{\mathrm{ext},\phi }$ is the electrical vector of external forces that contains the charge density $\overline{Q}$ .

The coupling to an electric circuit is easier if the electric surface current $-\stackrel{·}{\overline{Q}}$ appears in the equation. Therefore, the time derivative of equation (22) is used for the system simulations. Furthermore, a mechanical damping force ${\mathit{f}}^{\mathrm{damp},u}$ is introduced. The resulting system of equations thus becomes

$${\mathit{f}}^{\mathrm{dyn},u}+{\mathit{f}}^{\mathrm{damp},u}+{\mathit{f}}^{\mathrm{int},u}-{\mathit{f}}^{\mathrm{ext},u}=0$$

(23)

$${\stackrel{·}{\mathit{f}}}_{\mathrm{int},\phi }-{\stackrel{·}{\mathit{f}}}_{\mathrm{ext},\phi }=0$$

(24)

These equations are quite general as they allow to introduce nonlinearities in the material model (via ${\mathit{f}}^{\mathrm{int},u}$ or ${\stackrel{·}{\mathit{f}}}^{\mathrm{int},\phi }$ ), in the damping behavior (via ${\mathit{f}}^{\mathrm{damp},u}$ ) or in the electric circuit (via ${\stackrel{·}{\mathit{f}}}^{\mathrm{ext},\phi }$ ).

4. Linear system simulation

When applying the linear piezoelectric constitutive law (equations (8) and (9)) the forces in equations (23) and (24) depend linearly on the unknown mechanical displacements and electrical voltages. Utilizing $\mathit{u}$ and $\mathit{\phi }$ as global vectors that contain the nodal values of the displacements respectively of the electric voltages, the problem can be represented in a matrix formulation as

$$\begin{array}{c}\underset{\mathit{M}}{\underbrace{\left[\begin{array}{cc}{\mathit{M}}^{\mathrm{uu}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]}}\underset{\stackrel{··}{\mathit{x}}}{\underbrace{\left[\begin{array}{c}\stackrel{··}{\mathit{u}}\hfill \\ \multicolumn{1}{c}{\stackrel{··}{\mathit{\phi }}}\end{array}\right]}}+\underset{\mathit{D}}{\underbrace{\left[\begin{array}{cc}{\mathit{D}}^{\mathrm{uu}}\hfill & 0\hfill \\ \multicolumn{1}{c}{{\mathit{K}}^{\phi u}}& {\mathit{K}}^{\phi \phi }\hfill \end{array}\right]}}\underset{\stackrel{·}{\mathit{x}}}{\underbrace{\left[\begin{array}{c}\stackrel{·}{\mathit{u}}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\mathit{\phi }}}\end{array}\right]}}\hfill \\ +\underset{\mathit{K}}{\underbrace{\left[\begin{array}{cc}{\mathit{K}}^{\mathrm{uu}}\hfill & {\mathit{K}}^{u\phi }\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]}}\underset{\mathit{x}}{\underbrace{\left[\begin{array}{c}\mathit{u}\hfill \\ \multicolumn{1}{c}{\mathit{\phi }}\end{array}\right]}}=\underset{{\mathit{f}}^{\mathrm{ext}}}{\underbrace{\left[\begin{array}{c}{\mathit{f}}^{\mathrm{ext},u}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{f}}}^{\mathrm{ext},\phi }}\end{array}\right]}}\hfill \end{array}$$

(25)

whereby ${\mathit{D}}^{\mathrm{uu}}$ is the mechanical damping matrix and the stiffness submatrices ( ${\mathit{K}}^{\mathrm{uu}}$ : mechanical stiffness matrix, ${\mathit{K}}^{u\phi }$ and ${\mathit{K}}^{\phi u}$ : coupling matrices, ${\mathit{K}}^{\phi \phi }$ : capacitance matrix, and ${\mathit{M}}^{\mathrm{uu}}$ : mechanical mass matrix) are specified in the Appendix. The submatrices are assembled in the global mass matrix $\mathit{M}$ , the global damping matrix $\mathit{D}$ , and the global stiffness matrix $\mathit{K}$ . The global vector of external forces ${\mathit{f}}^{\mathrm{ext}}$ contains ${\mathit{f}}^{\mathrm{ext},u}$ and ${\stackrel{·}{\mathit{f}}}^{\mathrm{ext},\phi }$ . The global vector of the unknowns $\mathit{x}$ consists of $\mathit{u}$ and $\mathit{\phi }$ .

Within this contribution it is assumed that ${\mathit{f}}^{\mathrm{damp},u}$ in equation (23) represents a phenomenological velocity-proportional Rayleigh-type damping force. It can be constructed with ${\mathit{K}}^{\mathrm{uu}}$ , ${\mathit{M}}^{\mathrm{uu}}$ , and the velocity $\stackrel{·}{\mathit{u}}$ as

$${\mathit{f}}^{\mathrm{damp},u}=\underset{{\mathit{D}}^{\mathrm{uu}}}{\underbrace{[{\alpha }_R{\mathit{M}}^{\mathrm{uu}}+{\beta }_R{\mathit{K}}^{\mathrm{uu}}]}}\stackrel{·}{\mathit{u}}$$

(26)

with the Rayleigh-coefficients ${\alpha }_R$ and ${\beta }_R$ .

4.1. Coupling to electric circuits

Within a FE framework the existence of an electrode on a surface reduces the number of free voltage degrees of freedom to exactly one on this surface. All voltage degrees of freedom on an electrode must have the same value ${\phi }_{\mathrm{el}}$ . Therefore, one reference degree of freedom $\mathcal{F}$ is introduced for the electrode. Moreover, it is assumed that the electromechanical structure has at least a pair of electrodes and one of them is grounded. Hence, the voltage difference between the electrodes in this case is just the electrical voltage ${\phi }_{\mathrm{el}}$ of the non-grounded electrode with the corresponding degree of freedom $\mathcal{F}$ . The electric boundary conditions for the grounded electrode is therefore a homogeneous Dirichlet boundary condition with $\overline{\phi }=0$ . The boundary condition for the non-grounded electrode is used to couple the electromechanical structure to the electric circuit. The influence of the electric circuit on the electromechanical structure is modeled by prescribing the surface current ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ leaving the electrode. The current that is flowing from the electromechanical structure appears on the right hand side of equation (25). It can be constructed via an inhomogeneous Neumann boundary condition:

$${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}=\underset{\mathrm{electrode}}{\mathcal{A}}{\int }_{\partial {\mathrm{\Omega }}_{N\phi }^e}{\hat{N}}^I{\stackrel{·}{D}}_i{n}_i\mathrm{dA}$$

(27)

Both, ${\phi }_{\mathrm{el}}$ and ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ correspond to the same reference degree of freedom $\mathcal{F}$ of the non grounded electrode. Hence, all entries of ${\stackrel{·}{\mathit{f}}}^{\mathrm{ext},\phi }$ on the right hand side vanish except for the entry corresponding to the degree of freedom $\mathcal{F}$

$${\stackrel{·}{\mathit{f}}}^{\mathrm{ext},\phi }=\{\begin{array}{c}{\stackrel{·}{\overline{Q}}}_{\mathrm{el}}\mathrm{for}\mathrm{degree}\mathrm{of}\mathrm{freedom}=\mathcal{F}\hfill \\ \multicolumn{1}{c}{0\mathrm{for}\mathrm{degree}\mathrm{of}\mathrm{freedom}\ne \mathcal{F}}\end{array}$$

(28)

To consider the behavior of different electric circuits in the FE simulation, various functional relations between the current ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ and the electric voltage ${\phi }_{\mathrm{el}}$ at the electrode have to be considered. In general, the current can depend on the electric voltage and its time derivatives

$${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}={\stackrel{·}{\overline{Q}}}_{\mathrm{el}}({\phi }_{\mathrm{el}},{\stackrel{·}{\phi }}_{\mathrm{el}},{\stackrel{··}{\phi }}_{\mathrm{el}})$$

(29)

This functional relation will be specified for two different electric circuits in the example section. The outflowing current affects the mechanical displacement and the electric voltage as energy is extracted from the electromechanical structure. This modeling approach is extendable to multiple electrode configurations of a PVEH.

4.2. Time discretization

Energy harvesting is a dynamic application since its objective is converting vibrational (dynamic) energy to the electric domain. Modeling a PVEH and ignoring for example the resonance phenomenon or using static deflections patterns lead to fundamentally incorrect results (Erturk and Inman, 2008b). The importance of dynamic effects in energy harvesting necessitates the consideration of inertia and damping effects. The Bossak-Newmark scheme is applied to directly integrate equation (25) in time. This direct time integration scheme allows to capture all of the above mentioned effects. For linear problems the implicit Bossak-Newmark scheme is defined as (Wood et al., 1980)

$$\mathit{M}[[1+\alpha ]{\stackrel{··}{\mathit{x}}}_{n+1}-\alpha {\stackrel{··}{\mathit{x}}}_n]+\mathit{D}{\stackrel{·}{\mathit{x}}}_{n+1}+\mathit{K}{\mathit{x}}_{n+1}={\mathit{f}}_{n+1}^{\mathrm{ext}}$$

(30)

with

$${\stackrel{··}{\mathit{x}}}_{n+1}=\frac{1}{\beta \mathrm{\Delta }{t}^2}[{\mathit{x}}_{n+1}-{\mathit{x}}_n]-\frac{1}{\beta \mathrm{\Delta }t}{\stackrel{·}{\mathit{x}}}_n-[\frac{1}{2\beta }-1]{\stackrel{··}{\mathit{x}}}_n$$

(31)

$$\begin{array}{c}{\stackrel{·}{\mathit{x}}}_{n+1}=\frac{\gamma }{\beta \mathrm{\Delta }t}[{\mathit{x}}_{n+1}-{\mathit{x}}_n]-[\frac{\gamma }{\beta }-1]\hfill \\ {\stackrel{·}{\mathit{x}}}_n-[\frac{\gamma }{2\beta }-1]\mathrm{\Delta }t{\stackrel{··}{\mathit{x}}}_n\hfill \end{array}$$

(32)

Here, $n$ corresponds to the current time step, while $n+1$ denotes the next time step to be computed and $\mathrm{\Delta }t$ is the time step size. In order to make the time-integration scheme unconditionally stable and second order accurate, the parameters are specified as $\gamma =1/2+\alpha$ , $\beta =1/4[1+\alpha {]}^2$ . The choice of $\alpha \in ]0,1]$ introduces a slight numerical damping of higher frequencies. For more details about the Bossak-Newmark, e.g. (Hughes, 1987) scheme please refer to literature (e.g. R16). Introducing the approximations of the acceleration (31) and the velocity (32) into equation (30) results in a linear system of equations for the unknown displacements and electric voltages

$$\begin{array}{c}[\frac{1+\alpha }{\beta \mathrm{\Delta }{t}^2}\mathit{M}+\frac{\gamma }{\beta \mathrm{\Delta }t}\mathit{D}+\mathit{K}]{\mathit{x}}_{n+1}={\mathit{f}}_{n+1}^{\mathrm{ext}}+\alpha \mathit{M}{\stackrel{··}{\mathit{x}}}_n\hfill \\ +\mathit{M}[1+\alpha ][[\frac{1}{2\beta }-1]{\stackrel{··}{\mathit{x}}}_n+\frac{1}{\beta \mathrm{\Delta }t}{\stackrel{·}{\mathit{x}}}_n+\frac{1}{\beta \mathrm{\Delta }{t}^2}{\mathit{x}}_n]\hfill \\ +\mathit{D}[[\frac{\gamma }{2\beta }-1]\mathrm{\Delta }t{\stackrel{··}{\mathit{x}}}_n+[\frac{\gamma }{\beta }-1]{\stackrel{·}{\mathit{x}}}_n+\frac{\gamma }{\beta \mathrm{\Delta }t}{\mathit{x}}_n]\hfill \end{array}$$

(33)

Equation (33) allows for implicit transient simulations of linear electromechanical structures along with linear electric circuits. The presented system simulation method is implemented based on the open source C++ software library deal.ii targeted at the computational solution of partial differential equations (Arndt et al., 2019).

4.3. Linear electromechanical structure with standard electric circuit

A harvesting structure is oscillating during operation and naturally produces alternating electrical signals. Charging a storage component such as a battery or a capacitor requires a stable DC signal (Erturk and Inman, 2011). This necessitates the use of an AC-DC converter in the electric circuit. The simplest way of doing this is to use a diode bridge to supply the rectified voltage to an energy storage element (Gedeon and Rupitsch, 2018). This standard electric circuit is shown in Figure 1 and consists of a full-wave diode rectifier and a battery for energy storage. The electromechanical structure is open circuit and ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ vanishes if (Lefeuvre et al., 2009):

$$\mathrm{Case}1:|{\phi }_{\mathrm{el}}|<V_C\Rightarrow {\stackrel{·}{\overline{Q}}}_{\mathrm{el}}=0$$

(34)

Figure 1.

Unimorph PVEH with a standard electric circuit.

When $|{\phi }_{\mathrm{el}}|$ reaches the conductive voltage $V_C$ the diode bridge conducts and ${\phi }_{\mathrm{el}}$ is prescribed to be equal to the rectified voltage $V_C$ (Shu et al., 2007)

$$\mathrm{Case}2:|{\phi }_{\mathrm{el}}|=V_C\Rightarrow {\stackrel{·}{\overline{Q}}}_{\mathrm{el}}\ne 0$$

(35)

The conductive voltage $V_C$ results from the drop voltage of the diodes $V_{\mathrm{drop}}$ and the constant voltage $V_{\mathrm{DC}}$ as

$$|V_C|=2V_{\mathrm{drop}}+V_{\mathrm{DC}}$$

(36)

The typical waveforms of ${\phi }_{\mathrm{el}}$ and ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ when applying a standard circuit to a harmonically excited electromechanical structure are illustrated in Figure 2. ${\phi }_{\mathrm{el}}$ varies in time when the diode bridge is blocking and is kept equal to $V_C$ respectively $-V_C$ when the bridge conducts. ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ vanishes except when the diode bridge conducts.

Figure 2.

Typical waveforms of ${\phi }_{\mathrm{el}}$ and ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ for the standard circuit. The figure is based on a source in Shu et al. (2007).

Within the FE based system simulation, the difference between the standard circuit and the open circuit mode is the temporary limitation of ${\phi }_{\mathrm{el}}$ to $\pm V_C$ . This is accomplished in the following way:

• Case 1: This is the open circuit case, which can be simply modeled via a homogenous Neumann boundary condition in a FE simulation.

$${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}=0$$

(37)

• Case 2: In this case ${\phi }_{\mathrm{el}}$ is kept equal to $V_C$ respectively $-V_C$ . It is implemented with an inhomogeneous Dirichlet boundary condition by prescribing ${\phi }_{\mathrm{el}}$ .

$${\overline{\phi }}_{\mathrm{el}}=\pm V_C$$

(38)

In case 2 the current can be computed in a postprocessing step using equation (27). Both cases can be simulated with a standard implementation, however, the switching between them has to be provided:

Each new time step $n+1$ is computed with the homogenous Neumann boundary condition (Case 1). If $|{\phi }_{\mathrm{el},n+1}|\le V_C$ this assumption is correct and the next time step is computed. If $|{\phi }_{\mathrm{el},n+1}|>V_C$ the time step has to be recomputed with the prescribed electric voltage ${\overline{\phi }}_{\mathrm{el},n+1}=\pm V_C$ (Case 2). This logic for the standard circuit is provided in Figure 3.

Figure 3.

Logic for the definition of boundary conditions for the FE simulation of an electromechanical structure with a standard electric circuit.

As an application example the unimorph electromechanical structure from Erturk and Inman (2008a) is used and a standard electric circuit is added. Figure 1 shows the considered PVEH. A constant voltage $V_{\mathrm{DC}}$ is assumed to emulate an ideal energy storage device and a drop voltage of the diodes $V_{\mathrm{drop}}$ is defined. Therefore, non-ideal behavior of the diodes is modeled and the dissipation in the diodes is considered. In Leadenham and Erturk (2020) and Rupp et al. (2010) the diodes are modeled via the Shockley diode equation, which allows to simulate a more realistic diode behavior. This model could also be included in the present framework, but for simplicity a constant drop voltage is considered here. Table 1 presents the material properties of the substructure and geometrical and electrical parameters. The material data for PZT-5A is given in the Appendix.

Table 1.

Material properties of the substructure and geometrical parameters of the PVEH from Erturk and Inman (2008a). Furthermore, the parameters of the electric circuit are listed.

Width of the beam (mm)	20
Length of the beam (mm)	100
Young’s modulus substructure (GPa)	100
Thickness of the substructure (mm)	0.5
Thickness of the PZT (mm)	0.4
Density of the substructure (kg/ ${\mathrm{m}}^3$ )	7165
${\alpha }_R$ (rad/s)	4.886
${\beta }_R$ (s/rad)	1.2433 ×  ${10}^{-5}$
$V_{\mathrm{drop}}$ (V)	0.6
$V_{\mathrm{DC}}$ (V)	1.8

The PVEH is excited at its base with a harmonic base acceleration $a(t)$ with a magnitude of 1 m/s ² and a frequency of 48.7 Hz. The excitation frequency is close to the first eigenfrequency of the piezoelectric structure. The unimorph PVEH is discretized with $120$ quadratic hexahedral elements, the time step size $\mathrm{\Delta }t$ is set to $0.1$  ms and the parameter $\alpha$ of the Bossak-Newmark scheme is set to $0.5$ , introducing a slight numerical damping.

Figure 4 shows the time response of ${\phi }_{\mathrm{el}}$ and the time signal of the harvested energy $\mathcal{E}$ of the considered PVEH. The harvested energy can be computed as

$$\mathcal{E}={\int }_0^t{V}_{\mathrm{DC}}|{\stackrel{·}{\overline{Q}}}_{\mathrm{el}}|d\tau$$

(39)

Figure 4.

Time signal of ${\phi }_{\mathrm{el}}$ and the harvested energy of the unimorph PVEH with a standard electric circuit.

At time $t=0$ the system is at rest and the harmonic base acceleration starts. Due to the instationary settling process ${\phi }_{\mathrm{el}}$ flutters for $t<20$  ms. For the first time $|{\phi }_{\mathrm{el}}|$ reaches $|V_C|$ at around 100 ms. At this time the PVEH actually starts to harvest energy. At $t=200$  ms around $0.0023$  mJ are harvested by the PVEH.

Furthermore, the unimorph electromechanical structure was simulated in open circuit mode and for different electric load resistances in both time and frequency domain. The FE computations agree well with the results of the analytical theory of Erturk and Inman (2008a), which verifies that the presented numerical approach gives correct results.

5. Nonlinear system simulation

A linear description of the energy harvesting problem is no longer possible when more sophisticated electric circuits are considered, or nonlinearities of the electromechanical structure have to be taken into account.

In these cases, the particular forces in equations (23) and (24) do not depend linearly on the displacements and the electric voltage anymore. Therefore, the following nonlinear residuum equation has to be solved

$$\mathit{r}=\underset{{\mathit{f}}^{\mathrm{dyn}}}{\underbrace{\left[\begin{array}{c}{\mathit{f}}^{\mathrm{dyn},u}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]}}+\underset{{\mathit{f}}^{\mathrm{damp}}}{\underbrace{\left[\begin{array}{c}{\mathit{f}}^{\mathrm{damp},u}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{f}}}^{\mathrm{int},\phi }}\end{array}\right]}}+\underset{{\mathit{f}}^{\mathrm{int}}}{\underbrace{\left[\begin{array}{c}{\mathit{f}}^{\mathrm{int},u}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]}}-\underset{{\mathit{f}}^{\mathrm{ext}}}{\underbrace{\left[\begin{array}{c}{\mathit{f}}^{\mathrm{ext},u}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{f}}}^{\mathrm{ext},\phi }}\end{array}\right]}}=0$$

(40)

The coupling to nonlinear electric circuits is included in the external force ${\mathit{f}}^{\mathrm{ext}}$ and the nonlinear elastic behavior of the electromechanical structure enters the internal force ${\mathit{f}}^{\mathrm{int},u}$ . Since no further nonlinearities are taken into account, the inertial and the damping forces are similar as in the linear case, compare equation (25).

5.1. Time discretization and solution

The Bossak-Newmark method applied to the nonlinear equation (40) results in

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{r}}_{n+1}=[1+\alpha ]{\mathit{f}}_{n+1}^{\mathrm{dyn}}-\alpha {\mathit{f}}_n^{\mathrm{dyn}}+{\mathit{f}}_{n+1}^{\mathrm{damp}}+{\mathit{f}}_{n+1}^{\mathrm{int}}-{\mathit{f}}_{n+1}^{\mathrm{ext}}=0}\end{array}$$

(41)

Here, ${\mathit{r}}_{n+1}$ is the residual at time $t_{n+1}$ . For $\stackrel{·}{\mathit{x}}$ and $\stackrel{··}{\mathit{x}}$ the same approximations like in the linear case (equations (31) and (32)) are used. Due to its nonlinearity the equation is solved iteratively by the Newton-Raphson method. The solution increment is derived from

$${}^{(k+1)}\mathit{r}{=}^{(k)}\mathit{r}+\underset{{}^{(k)}\mathcal{T}}{\underbrace{\frac{{\partial }^{(k)}\mathit{r}}{\partial \mathit{x}}}}·\mathrm{\Delta }\mathit{x}=0{\Rightarrow }^{(k)}\mathcal{T}·\mathrm{\Delta }\mathit{x}={-}^{(k)}\mathit{r}$$

(42)

Here, ${}^{(k)}\mathcal{T}$ is the tangent stiffness matrix and $(k)$ is the iteration index. The update for $\mathit{x}$ is then obtained as

$${}^{(k+1)}\mathit{x}{=}^{(k)}\mathit{x}+\mathrm{\Delta }\mathit{x}$$

(43)

Subsequently, $\mathit{r}$ is computed again with the updated $\mathit{x}$ and the procedure is repeated until the norm of $\mathit{r}$ approximately vanishes. The tangent stiffness matrix ${}^{(k)}\mathcal{T}$ can be computed with the chain rule

$${}^{(k)}\mathcal{T}=\underset{{{}^{(k)}\mathcal{T}}_{\mathcal{A}}}{\underbrace{\frac{{\partial }^{(k)}\mathit{r}}{\partial \stackrel{··}{\mathit{x}}}}}\frac{\partial \stackrel{··}{\mathit{x}}}{\partial \mathit{x}}+\underset{{{}^{(k)}\mathcal{T}}_{\mathcal{B}}}{\underbrace{\frac{{\partial }^{(k)}\mathit{r}}{\partial \stackrel{·}{\mathit{x}}}}}\frac{\partial \stackrel{·}{\mathit{x}}}{\partial \mathit{x}}+\underset{{{}^{(k)}\mathcal{T}}_{\mathcal{C}}}{\underbrace{\frac{{\partial }^{(k)}\mathit{r}}{\partial \mathit{x}}}}$$

(44)

with the tangent submatrices ${{}^{(k)}\mathcal{T}}_{\mathcal{A}}$ , ${{}^{(k)}\mathcal{T}}_{\mathcal{B}}$ , and ${{}^{(k)}\mathcal{T}}_{\mathcal{C}}$ . Introducing equations (31) and (32) results in

$${}^{(k)}\mathcal{T}={\frac{1+\alpha }{\beta \mathrm{\Delta }{t}^2}}^{(k)}{\mathcal{T}}_{\mathcal{A}}+{\frac{\gamma }{\beta \mathrm{\Delta }t}}^{(k)}{\mathcal{T}}_{\mathcal{B}}{+}^{(k)}{\mathcal{T}}_{\mathcal{C}}$$

(45)

${}^{(k)}\mathcal{T}$ must be computed and equation (42) has to be solved in every Newton-Raphson iteration. The tangent submatrices ${{}^{(k)}\mathcal{T}}_{\mathcal{A}}$ , ${{}^{(k)}\mathcal{T}}_{\mathcal{B}}$ , and ${{}^{(k)}\mathcal{T}}_{\mathcal{C}}$ are specified in the Appendix.

5.2. Linear electromechanical structure with SSHI circuit

The so called synchronized switch damping technique was introduced by Richard et al. (1999) to address the problem of vibration damping for mechanical structures. Based on this technique Guyomar et al. (2005) derived the SSHI circuit for energy harvesting applications. It was reported that the SSHI circuit can increase the harvested power by up to 900% as compared to a standard circuit depending on both, the electromechanical coupling of the system as well as the excitation signal (Badel et al., 2006; Guyomar et al., 2005).

Figure 5 presents the unimorph electromechanical structure from the previous example coupled to an SSHI circuit. The SSHI circuit adds a switch and an inductance ( $L$ ) to the standard circuit. Furthermore, a parasitic resistance in the inductance ( $R_L$ ) is considered because it can strongly influence the behavior and efficiency of a PVEH with an SSHI circuit (Guyomar and Lallart, 2011).

Figure 5.

Unimorph PVEH with an SSHI circuit.

Both, the inductance and the switch are connected in parallel to the piezoelectric element. A self-adaptive control of the switch is applied that triggers on the maxima and minima of ${\phi }_{\mathrm{el}}$ .

The switch is closed when $|{\phi }_{\mathrm{el}}|$ starts to decrease. When the switch is closed the following condition is fulfilled

$$\mathrm{Case}3:|{\phi }_{\mathrm{el}}|<V_C\Rightarrow {\stackrel{··}{\overline{Q}}}_{\mathrm{el}}=\frac{{\phi }_{\mathrm{el}}-R_L{\stackrel{·}{\overline{Q}}}_{\mathrm{el}}}{L}$$

(46)

When the switch is closed the rectifier bridge is in open circuit mode and the inductance is connected to the electromechanical structure. The piezoelectric element acts like a capacitance and thus results in an electric resonant circuit consisting of the piezoelectric capacitance and the inductance. Hence, the electric voltage and the electric current start to oscillate as soon as the switch is closed. The switch is opened again when the current in the inductor is zero-crossing, which results in an inversion of ${\phi }_{\mathrm{el}}$ . Because the eigenfrequency of the electric resonant circuit is usually much higher than the mechanical frequencies, the inversion of ${\phi }_{\mathrm{el}}$ is almost instantaneous.

Figure 6 presents the typical waveform of ${\phi }_{\mathrm{el}}$ and ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ arising from an SSHI circuit. At $t=0$ the switch is opened and because $|{\phi }_{\mathrm{el}}|<V_C$ the electromechanical structure is in open circuit mode and ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ vanishes. When $|{\phi }_{\mathrm{el}}|$ reaches $V_C$ the diode bridge conducts and ${\phi }_{\mathrm{el}}$ is kept equal to $\pm V_C$ . Hence, there is an electric current ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}\ne 0$ . When $|{\phi }_{\mathrm{el}}|$ starts to decrease the switch is closed, which leads to an inversion of ${\phi }_{\mathrm{el}}$ . When the current in the inductance crosses zero the switch is opened again. Through the almost instantaneous inversion of ${\phi }_{\mathrm{el}}$ , the SSHI circuit aims to extend the time fraction during which the diodes are conducting.

Figure 6.

Typical waveforms of ${\phi }_{\mathrm{el}}$ and ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ for an electric circuit of the SSHI circuit. The figure is based on a source in Shu et al. (2007).

The situation when the switch is opened is simulated similarly as in the previous example (Cases 1 and 2). When the switch is closed the influence of the inductance on the electromechanical structure is considered via an inhomogeneous Neumann boundary condition, integrating equation (46) in time (Case 3).

By applying the implicit second-order accurate trapezoidal rule, ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ can be computed as

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\overline{Q}}}_{\mathrm{el},n+1}=\frac{2L-\mathrm{\Delta }t{R}_L}{2L+\mathrm{\Delta }t{R}_L}{\stackrel{·}{\overline{Q}}}_{\mathrm{el},n}}\\ +\frac{\mathrm{\Delta }t}{2L+\mathrm{\Delta }t{R}_L}[{\phi }_{\mathrm{el},n+1}+{\phi }_{\mathrm{el},n}]\hfill \end{array}$$

(47)

This relation for ${\stackrel{·}{\overline{Q}}}_{\mathrm{el},n+1}$ is included in the external force vector ${\stackrel{·}{\mathit{f}}}_{n+1}^{\mathrm{ext},\phi }$ in equation (40). The tangent submatrix ${\mathcal{T}}_{\mathcal{C}}$ , compare equation (A11), which is needed for the iterative solution of the problem, can now be specified by means of the derivative

$$\frac{\partial {\stackrel{·}{\overline{Q}}}_{\mathrm{el},n+1}}{\partial {\phi }_{\mathrm{el},n+1}}=\frac{\mathrm{\Delta }t}{2L+\mathrm{\Delta }t{R}_L}$$

(48)

which is placed on the main diagonal of the tangent submatrix for the corresponding degree of freedom $\mathcal{F}$ of the electrode, compare equation (28).

In the following, the logic of the self adaptive SSHI circuit is described: Assume that the switch is opened at time $n$ . Then, timestep $n+1$ is also computed with opened switch and in open circuit mode (Case 1) with three possible results:

• If $|{\phi }_{\mathrm{el},n+1}|<V_C$ but $|{\phi }_{\mathrm{el},n+1}|>|{\phi }_{\mathrm{el},n}|$ the solution is accepted and the next time step is simulated (Case 1).

• If $|{\phi }_{\mathrm{el},n+1}|>V_C$ timestep $n+1$ is recomputed with prescribed ${\phi }_{\mathrm{el},n+1}$ to $\pm V_C$ (Case 2).

• If $|{\phi }_{\mathrm{el},n+1}|$ decreases, such that $|{\phi }_{\mathrm{el},n+1}|<|{\phi }_{\mathrm{el},n}|$ , timestep $n+1$ is recomputed with closed switch (Case 3).

If the switch is closed at timestep $n$ the next time step is computed with closed switch (Case 3) with two possible results:

• If the direction of ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ does not change, that is $\mathrm{sgn}({\stackrel{·}{\overline{Q}}}_{\mathrm{el},n+1})=\mathrm{sgn}({\stackrel{·}{\overline{Q}}}_{\mathrm{el},n})$ (sgn is the signum function), the switch remains closed and the next timestep is simulated (Case 3).

• If the current in the inductance is zero-crossing, that is $\mathrm{sgn}({\stackrel{·}{\overline{Q}}}_{\mathrm{el},n+1})\ne \mathrm{sgn}({\stackrel{·}{\overline{Q}}}_{\mathrm{el},n})$ timestep $n+1$ is recomputed with opened switch in open circuit mode (Case 1).

In Figure 7 the logic of the SSHI circuit is summarized. Furthermore, it has to be prevented that the switch is triggered at inappropriate times. After closing the switch, ${\phi }_{\mathrm{el}}$ is inverted nearly instantaneously because the electrical oscillation frequency is by orders higher than the mechanical excitation frequency. The instantaneous inversion of ${\phi }_{\mathrm{el}}$ acts like an actuation on the electromechanical structure, which leads to an oscillation with a higher order frequency of both, ${\phi }_{\mathrm{el}}$ as well as the displacement of the structure. To prevent that the switch is triggered because of oscillations of ${\phi }_{\mathrm{el}}$ , caused by the process of voltage inversion, a time span $t_{\mathrm{ns}}$ after the switch has been opened is defined during which closing of the switch is prohibited. Furthermore, closing the switch during the instationary settling process is similarly prohibited during a time span $t_{\mathrm{switch}}$ in the beginning of the simulation.

Figure 7.

Logic for the definition of boundary conditions for the FE simulation of an electromechanical structure with an SSHI circuit.

To validate the proposed simulation method a PVEH with SSHI circuit is simulated and the results are compared with results from literature. To this end, the same unimorph electromechanical structure from the previous chapter along with an SSHI circuit ( $L$  = 0.1 mH and $R_L=0$   $\mathrm{\Omega }$ ) is considered. The model considered here was originally utilized by Gedeon and Rupitsch (2018) to show the applicability of their respective approach.

In order to efficiently and accurately simulate the system, variable time step sizes are introduced. The voltage inversion (during which the switch is closed) requires the smallest time step of ${10}^{-4}$  ms since it happens almost instantaneously. After the switch is opened the time step size is set to ${10}^{-3}$  ms within a time span $t_{\mathrm{osc}}=0.25$  ms to capture the higher order oscillation of the electromechanical structure with sufficient precision. For the remaining time period a time step size of ${10}^{-1}$  ms is used. In addition, the time spans $t_{\mathrm{ns}}$ and $t_{\mathrm{switch}}$ are chosen to $4$ and $10.8$  ms, respectively. The same discretization of the unimorph PVEH as in the previous chapter with $120$ quadratic hexahedral elements is used. The parameter $\alpha$ of the Bossak-Newmark scheme is set to $0.5$ .

Figure 8 compares the resulting electric voltage computed with the presented FE simulation against the results obtained from a simulink simulation of a reduced model of Gedeon and Rupitsch (2018) for the SSHI circuit. The largest differences between the results of the presented method and this reduced model occur during the instationary settling process. After a time span of about 45 ms a very good agreement of the two models can be observed. This overall good agreement verifies that the presented FE system simulation approach is correct.

Figure 8.

Comparison of ${\phi }_{\mathrm{el}}$ for the SSHI circuit of the reduced model of Gedeon and Rupitsch (2018) and the presented method.

Furthermore, Figure 9 shows ${\phi }_{\mathrm{el}}$ and the harvested energy of the PVEH with the SSHI circuit up to 200 ms. The actual energy harvesting starts at around 60 ms when $|{\phi }_{\mathrm{el}}|$ reaches $V_C$ for the very first time. The harvested energy of the SSHI circuit after $200$  ms is around $0.0033$  mJ which is approximately 44% more energy harvested than with the standard circuit.

Figure 9.

Time signals of ${\phi }_{\mathrm{el}}$ and the harvested energy of the unimorph PVEH with an SSHI circuit till $200$  ms.

To investigate the influence of the parasitic resistance of the inductance $R_L$ on the behavior of the unimorph PVEH with SSHI circuit, four different parasitic resistances $R_L$ are introduced: $0$ , $0.5$ , $1.5$ , and $10$   $\mathrm{\Omega }$ . Figure 10 presents ${\phi }_{\mathrm{el}}$ for the time period 0–200 ms for the different values of $R_L$ and a constant inductance $L=0.1$  mH. Nearly the same evolution of ${\phi }_{\mathrm{el}}$ is obtained for $R_L=0$ , $R_L=0.5$ , and $R_L=1.5$ $\mathrm{\Omega }$ . Only the evolution of ${\phi }_{\mathrm{el}}$ for $R_L=10$   $\mathrm{\Omega }$ differs significantly: After voltage inversion, $|{\phi }_{\mathrm{el}}|$ for $R_L=10$   $\mathrm{\Omega }$ is significantly smaller than $|{\phi }_{\mathrm{el}}|$ for $R_L=0$   $\mathrm{\Omega }$ , as can be observed in the zoom-in of Figure 10. Moreover, after a time period of $200$  ms, the configurations with $R_L=0$ , $R_L=0.5$ , $R_L=1.5$   $\mathrm{\Omega }$ harvest around 0.0033 mJ, while the configuration with $R_L=10$   $\mathrm{\Omega }$ harvests around 0.0031 mJ, resulting in a reduction of harvested energy of about 9%. The parasitic resistance of the inductance $R_L$ damps the electric oscillation of ${\phi }_{\mathrm{el}}$ , thus, with increasing $R_L$ the value of ${\phi }_{\mathrm{el}}$ after voltage inversion decreases. The lower $R_L$ is, the lower are the electric losses during inversion of ${\phi }_{\mathrm{el}}$ . At low values of $R_L$ , the process of voltage inversion is efficient and leads to a high efficiency of the SSHI circuit.

Figure 10.

Time signals of ${\phi }_{\mathrm{el}}$ of the unimorph PVEH with an SSHI circuit till $200$  ms for different values of $R_L$ .

For the considered inductance ( $L=0.1$  mH), the parasitic resistance $R_L$ is usually smaller than $0.5$   $\mathrm{\Omega }$ . Figure 10 illustrates, that the results for such small values of $R_L$ do not significantly differ from the results when $R_L$ is not considered. Therefore, an ideal inductance without a parasitic resistance $R_L$ is considered in the following.

In the open circuit mode, the SSHI circuit works as a Synchronized Switch Damping and thus reduces the displacements (Richard et al., 2000). The presented simulation framework is fully coupled and accounts for this damping phenomena, in addition to the damping which results from the harvested energy. To illustrate this functionality of the simulation framework, a unimorph electromechanical structure with an electric circuit of the synchronized switch damping on inductor (SSDI) technique is simulated. The SSDI circuit consists of a switch and an inductance, compare to Figure 11. Similar to the SSHI technique, the switch is closed when $|{\phi }_{\mathrm{el}}|$ reaches a maximum value and subsequently opened when ${\phi }_{\mathrm{el}}$ is inverted. When operating in open circuit mode the SSHI circuit is equivalent to the SSDI technique.

Figure 11.

Unimorph electromechanical structure with an SSDI circuit.

In the simulations, the same unimorph electromechanical structure with the same mechanical, electrical, and excitation parameters as in the previous example is used. Figure 12 presents the piezoelectric voltage ${\phi }_{\mathrm{el}}$ of the open circuit configuration and the SSDI configuration. Within the considered time, the SSDI circuit generates higher voltages than the open circuit configuration because of the rapid inversion of the piezoelectric voltage ${\phi }_{\mathrm{el}}$ . Figure 13 compares the tip displacement of the unimorph electromechanical structure in open circuit mode against the SSDI configuration. The tip displacement of the SSDI configuration is significantly smaller than the tip displacement of the open circuit configuration. Because of the rapid inversion and magnification of the piezoelectric voltage ${\phi }_{\mathrm{el}}$ the SSDI damps the vibration of the structure (Richard et al., 1999, 2000). Moreover, the influence of the SSDI technique changes the resonance frequency. Consequently, the tip displacement is additionally reduced, because the excitation frequency is not close to the eigenfrequency of the SSDI configuration. These coupling effects between the mechanical and electrical fields are completely captured by the presented simulation framework, as illustrated in Figures 12 and 13.

Figure 12.

Comparison of the piezoelectric voltage ${\phi }_{\mathrm{el}}$ of the open circuit configuration and the SSDI configuration.

Figure 13.

Comparison of the tip displacement of the open circuit configuration and the SSDI configuration.

5.3. Nonlinear electromechanical structure with SSHI circuit

In the following the applicability of the presented system simulation approach for nonlinear electromechanical structures combined with nonlinear electric circuitry is demonstrated. An application example is considered that is restricted to nonlinear elasticity since this is the primary source of nonlinearity in energy harvesting applications (Stanton et al., 2012). Other possible sources of nonlinearities like geometric nonlinearities (Behjat and Khoshravan, 2012), nonlinear damping (Stanton et al., 2012; Yang and Upadrashta, 2016) as well as strong electric fields (Beige, 1983), and others are neglected.

The SSHI circuit from the previous chapter is added to the bimorph electromechanical structure introduced in Erturk and Inman (2009). The cantilevered bimorph electromechanical structure consists of two layers of PZT-5A (material parameters in the Appendix) bracketing a passive substructure layer. The piezoceramic layers are poled in opposite directions. To consider nonlinear elasticity, the nonlinear constitutive law for PZT-5A is applied (equations (9) and (12)). The substructure is assumed to be linear elastic. Figure 14 presents the considered bimorph PVEH, which has a mass mounted on its tip, and Table 2 provides its parameters. The same frequency $47.8$  Hz of the harmonic base acceleration $a(t)$ as for the unimorph harvester is applied, since this excitation frequency is close to the first open circuit resonance frequency of the bimorph harvester.

Figure 14.

Bimorph PVEH with SSHI circuit.

Table 2.

Parameters of the bimorph electromechanical structure from Erturk and Inman (2009). Furthermore, the parameters of the SSHI circuit are listed.

Width of the beam (mm)	31.8
Length of the beam (mm)	50.8
Young’s modulus substructure (GPa)	105
Thickness of the substructure (mm)	0.14
Thickness PZT (mm)	0.26 (each)
Density of the substructure (kg/ ${\mathrm{m}}^3$ )	9000
Rayleigh coefficient ${\alpha }_R$ (rad/s)	14.65
Rayleigh coefficient ${\beta }_R$ (s/rad)	${10}^{-5}$
Tip mass $m_t$ (kg)	0.012
Inductivity (mH)	0.1
$V_{\mathrm{drop}}$ (V)	0.6
$V_{\mathrm{DC}}[V]$	1.8

Like in the previous chapter ${\stackrel{·}{\overline{Q}}}_{\mathrm{el}}$ is computed via equation (47) when the switch of the SSHI circuit is closed. Furthermore, the nonlinear elastic material model (12) is used here, which has to be considered in the internal forces and the tangent submatrix ${\mathcal{T}}_{\mathcal{C}}$ , compare equation (A11) in Appendix. In each Newton-Raphson iteration, within the solution for time step $n+1$ , the following system of equations has to be solved

$$[\frac{1+\alpha }{\beta \mathrm{\Delta }{t}^2}\mathit{M}+\frac{\gamma }{\beta \mathrm{\Delta }t}\mathit{D}+\left[\begin{array}{cc}\frac{\partial {\mathit{f}}_{n+1}^{\mathrm{int},u}}{\partial \mathit{u}}\hfill & {\mathit{K}}^{u\phi }\hfill \\ \multicolumn{1}{c}{0}& {[-\frac{\mathrm{\Delta }t}{2L}]}^{\mathcal{F}}\hfill \end{array}\right]]\mathrm{\Delta }\mathit{x}=-{\mathit{r}}_{n+1}$$

(49)

In the simulation of the bimorph PVEH $t_{\mathrm{osc}}=$  3 ms is applied, because due to the tip mass the higher order vibration mode induced by the process of voltage inversion decays slower than in the example with the unimorph electromechanical structure without a tip mass. To prohibit triggering the switch at inappropriate times $t_{\mathrm{ns}}=4$  ms and $t_{\mathrm{switch}}=10.8$  ms are chosen, similarly to the unimorph PVEH application example. The bimorph PVEH is discretized with $90$ quadratic hexahedral elements. The parameter $\alpha$ of the Bossak-Newmark time integration scheme is set to $0.1$ , which is for this example sufficient to damp high frequencies.

Firstly, a harmonic base excitation of $1$  m/s ² is applied. Figure 15 presents the results for the harvested energy and for the electric voltage ${\phi }_{\mathrm{el}}$ for the bimorph PVEH. Within 200 ms the bimorph PVEH generates around $0.0065$  mJ which is approximately 97% more than the amount of energy generated by the unimorph PVEH with an SSHI circuit.

Figure 15.

Harvested energy and ${\phi }_{\mathrm{el}}$ of the bimorph PVEH with SSHI circuit and nonlinear electromechanical structure under a harmonic base excitation with an acceleration of $1$  m/s ².

To analyze the influence of nonlinear elasticity on the behavior of the PVEH, simulations with the linear constitutive law are performed and compared to the results with the nonlinear constitutive law. Since the magnitude of the harmonic base excitation of $1$  m/s ² is relatively small, the results of the bimorph PVEH with the linear constitutive law agree well with the results shown in Figure 15, and are not shown here. The simulation with the linear piezoelectric constitutive law results in $0.5$ % higher energy generation than the simulation with the nonlinear constitutive law.

To further analyze the nonlinearities induced by the nonlinear elasticity model, a higher harmonic base acceleration of $9.81$  m/s ² is applied to the bimorph PVEH. Firstly, the behavior in open circuit mode is studied. Figure 16 compares the open circuit voltage of the simulations with linear and nonlinear elastic PZT-5A behavior. The influence of nonlinear elasticity leads not only to a reduction of the magnitudes of ${\phi }_{\mathrm{el}}$ but also induces a phase shift compared to the linear solution.

Figure 16.

Comparison of ${\phi }_{\mathrm{el}}$ of the linear and nonlinear bimorph electromechanical structure in open circuit mode, accelerated with a harmonic base excitation of 9.81 m/s ².

Figure 17 presents the comparison of the electric voltage ${\phi }_{\mathrm{el}}$ of the PVEH coupled to an SSHI circuit with the nonlinear respective linear constitutive law. Because of the high mechanical excitation level, the advantage of the SSHI circuit, namely extending the conduction time of the diode bridge, does not significantly improve the efficiency of the PVEH. However, like in the open circuit configuration, a phase difference of ${\phi }_{\mathrm{el}}$ between the linear and nonlinear material model arises.

Figure 17.

Comparison of ${\phi }_{\mathrm{el}}$ of the linear and nonlinear bimorph electromechanical structure with SSHI circuit excited with a harmonic base acceleration of $9.81$  m/s ².

Figure 18 compares the harvested energy of the PVEH with nonlinear constitutive law with the harvested energy of the PVEH with linear constitutive law. A significantly reduced amount of energy is harvested by the nonlinear electromechanical structure. In this application example the nonlinear configuration harvests within $200$  ms around $0.095$  mJ while the simulation with the linear constitutive law renders $0.140$  mJ. Neglecting nonlinear material behavior for high levels of base excitation therefore leads to a significant overestimation in the computation of the harvested energy.

Figure 18.

Harvested energy of the bimorph PVEH with nonlinear electromechanical structure respective linear electromechanical structure under a harmonic base acceleration of $9.81$  m/s ².

In this contribution a DC battery, which is held at constant 1.8 V, is charged by the standard and SSHI circuit. Another important aspect would be to use different electric circuits to charge a storage capacitor instead of a battery (Zhang et al., 2021). We will apply the presented simulation framework in a future work to simulate the influence of a storage capacitor on the efficiency of PVEHs.

6. Conclusion

A novel finite element based system simulation approach for nonlinear electromechanical structures coupled to nonlinear electric circuits was introduced. The proposed solution scheme allows for arbitrary nonlinear effects, for example due to the material model, due to damping or due to the electric circuit. Using the Newton-Raphson method these nonlinearities are taken into account in a consistent manner. However, here nonlinearities are restricted to the material model and the electric circuit.

The presented system simulation approach is verified against an unimorph PVEH with an SSHI circuit from literature. The FE-based system simulation approach allows to couple the same unimorph PVEH to a standard circuit, confirming the higher efficiency of the SSHI circuit to harvest energy. Moreover, the tip displacement and electric voltage of the unimorph structure with open circuit configuration and an SSDI configuration is compared. The simulation results demonstrate that the proposed simulation method takes all coupling effects between the mechanical and the electrical domain into account. A further application example of a bimorph PVEH reveals the advantage of the presented system simulation method to easily account for nonlinearities in the electromechanical structure. The bimorph PVEH with SSHI circuit is simulated at a high excitation level. Nonlinear elasticity is considered in the material model of the piezoceramic representing the main source of nonlinearity in energy harvesting applications. Comparisons of simulations with linear and nonlinear material behavior of the piezoceramic demonstrate that significantly less energy is harvested when the nonlinear behavior of the piezoceramic is taken into account. Thus, accounting for nonlinearities of the electromechanical structure when the generated energy has to be computed is paramount for high levels of mechanical excitation.

The presented holistic simulation approach necessitates no external circuit simulation software and allows for efficient time integration schemes of the underlying differential equations. While exploiting the structural simulation capabilities of FEM without restrictions, the proposed method is highly flexible regarding changes of the electromechanical structure or of the electric circuit. The presented approach allows to evaluate both mechanical and electrical quantities in detail and facilitates the rapid development of PVEHs.