An implicitly coupled finite element-electronic circuit simulator method for efficient system simulations of piezoelectric energy harvesters

Abstract

A piezoelectric vibration-based energy harvester (PVEH) consists of an electromechanical structure and an electric circuit, influencing each other. The finite element method (FEM) allows for the simulation of arbitrarily shaped and nonlinear electromechanical structures, while an electronic circuit simulator (ECS) can be applied to simulate arbitrary electric circuits. We propose a novel implicit FEM-ECS coupling method, which combines the advantages of the FEM and the flexibility of an ECS. The Newton-Raphson method is applied to achieve rapid convergence at the interface between the FEM and ECS simulations, so that arbitrary PVEHs can be efficiently analyzed. The new implicit FEM-ECS coupling is validated against explicitly and monolithically coupled FEM-ECS simulations of PVEHs from literature. Moreover, an application example consisting of a nonlinear electromechanical structure and a self-powered SSHI circuit is considered and a shock-like base excitation is applied. The examples demonstrate the practical applicability of the novel coupling method for realistic simulations of PVEHs.

1. Introduction

A piezoelectric vibration-based energy harvester (PVEH) is a combined device consisting of an electromechanical structure in conjunction with an electric circuit. Such a device converts ambient vibration energy into electric energy to drive low-power electronics (Erturk and Inman, 2011b). PVEHs provide clean energy and are an alternative to conventional non-environmentally friendly power sources such as batteries (Safaei et al., 2019). The piezoelectric effect, which describes an interaction between mechanical and electrical quantities, is exploited as the energy conversion mechanism of PVEHs (Rupitsch, 2018).

The most common electromechanical transducer in vibration energy harvesting is a cantilever beam (Safaei et al., 2019). Linear beam-type PVEHs for example in unimorph (a single piezoelectric layer bonded to a substrate layer) or bimorph (a single substrate layer with two symmetric piezoelectric layers bonded on each side) configuration are very efficient when excited at resonance and often a tip mass is added to further improve the energy output (Safaei et al., 2019). However, if the mechanical excitation of linear PVEHs does not exactly match the resonance frequency, the power output drops dramatically (Gammaitoni et al., 2009; Liang et al., 2021). To overcome this issue of the narrow bandwidth of linear PVEHs, nonlinearities are intentionally introduced resulting in a significant extension of the bandwidth (Erturk and Inman, 2011a; Leadenham and Erturk, 2015; Zhou et al., 2016). In addition to these intentionally introduced nonlinearities the electromechanical structure itself typically shows nonlinear constitutive and dissipative behavior for moderate to large levels of mechanical excitation (Leadenham and Erturk, 2020; Stanton et al., 2010, 2012). Consequently, nonlinearities of the electromechanical structure are an important characteristic of PVEHs and must be considered in their design and modeling.

An electric circuit is used to extract the energy from the electromechanical structure. During oscillation, the electromechanical structure naturally generates alternating electric signals (Erturk and Inman, 2011b). However, charging an electric storage component requires a DC signal, meaning that the electric circuit must perform an AC-to-DC conversion and is hence forced to be nonlinear. Due to the direct and indirect piezoelectric effect the electric circuit and the electromechanical structure mutually influence each other. Accordingly, an accurate simulation of a PVEH requires the modeling of both the nonlinear electromechanical structure and the nonlinear electric circuit.

Analytical models of PVEHs are generally limited to simple geometries and circuits, and offer little flexibility with respect to changes in boundary conditions (Erturk and Inman, 2008). As a consequence, numerical methods have evolved and are generally more flexible than analytical models, but their disadvantage is the high computational cost. In Elvin and Elvin (2009b) and Yang and Tang (2009) an electronic circuit simulator (ECS), which allows for the simulation of arbitrary electric circuits, is applied for system simulations of complete PVEHs. The mechanical degrees of freedom are decoupled via modal analysis, which limits the approaches to linear electromechanical structures. Furthermore, equivalent electric circuit models of PVEHs considering dissipative and stiffness nonlinearities are introduced in Silva et al. (2018) and Elvin (2014).

Often, the finite element method (FEM) is applied for the simulation of PVEHs and allows to solve the full set of 3D piezoelectric equations (Hegendörfer et al., 2022a; Yamamoto et al., 2021). FEM simulations of electromechanical structures are coupled with linear passive electric circuit elements in Wang and Ostergaard (1999), Zhu et al. (2009), De Marqui Junior et al. (2009), Abdelkefi et al. (2014), Akbar and Curiel-Sosa (2019), and Ishihara et al. (2021). However, practical electric circuits for PVEHs necessarily require nonlinear components. In Wu and Shu (2015) the equivalent steady state impedances of the standard and synchronized switch harvesting on inductor (SSHI) circuit are derived and directly included in an FEM software using available passive electric circuit elements.

Recently, in Hegendörfer et al. (2022b) a system simulation method that allows finite element simulations of nonlinear electromechanical structures coupled to nonlinear electric circuits is introduced. The FEM solves the complete set of piezoelectric equations while the electric circuit is considered via the boundary conditions. Since no ECS is used, realistic electric circuit components cannot be considered.

In summary, the FEM offers powerful simulation capabilities for the electromechanical structure, but provides only limited capabilities for the consideration of electric circuits. Conversely, an ECS can simulate arbitrary electric circuits, but offers only limited possibilities for simulating electromechanical structures. Therefore, there are efforts to combine the two simulation methods to overcome their inherent limitations.

In Elvin and Elvin (2009a) an explicit coupling between a purely mechanical FEM simulation and an ECS simulation for the electric circuit is presented. In the approach, the electric behavior of the electromechanical structure of a PVEH is simplified to a lumped capacitance, which may not be appropriate especially for strongly coupled electromechanical structures (Gedeon and Rupitsch, 2018). An explicit coupling means, that no iterations between the FEM and ECS simulations are conducted and equilibrium between the two simulation domains is not ensured. Consequently, small time step sizes are necessary to obtain accurate results, which leads to high computational costs.

Monolithic FEM-ECS coupling approaches reported in the literature (Aftab et al., 2012; Gedeon and Rupitsch, 2018; Xiang et al., 2018) exploit a multiphysics simulation environment, which includes an ECS. The algebraic equations stemming from an FEM discretization of the electromechanical structure are reduced and integrated along with the electric circuit model into one system using a multiphysics simulation environment. Subsequently, the resulting system can be solved consistently in the multiphysics simulation environment. This modeling strategy can take into account arbitrary electric circuits but nonlinearities of the electromechanical structure or changing mechanical boundary conditions cannot be considered.

Overall, the FEM-ECS coupling methods reported in the literature are limited to linear electromechanical structures or the electromechanical structure is simplified and the coupling between the FEM and the ECS simulation is not very efficient. The main advantage of the FEM, namely the capability to accurately simulate nonlinear electromechanical structures, cannot be realized. To address this problem, we present a novel implicit coupling method that allows to exploit the full capabilities of the FEM for simulating nonlinear electromechanical structures, while separately using an ECS for analyzing the electric circuit. An implicit coupling means, that equilibrium between the FEM and the ECS simulations is ensured at each time step by means of an iterative solution method.

Firstly, the piezoelectric equations along with a finite element framework are presented. Subsequently, the main contribution of this work, a novel implicit FEM-ECS coupling method, is presented. Next, the new coupling method is validated against simulations of a unimorph PVEH from literature. Finally, an application example demonstrating the practical applicability of the new implicit coupling method is considered.

2. FEM simulation of the electromechanical structure

For the simulation of the electromechanical structure a FEM framework is applied, which we already presented in Hegendörfer et al. (2022b). For completeness, the used piezoelectric equations are shortly summarized together with the FEM framework. Details can be found in Hegendörfer et al. (2022b). Later, this FEM framework is coupled via a novel implicit method to an ECS.

In this contribution electromechanical structures with a linear and a nonlinear piezoelectric constitutive law are considered. The linear piezoelectric constitutive law relates the mechanical stresses $T_{\mathrm{ij}}$ and the dielectric displacements $D_i$ with the mechanical strain $S_{\mathrm{ij}}$ and the electric field $E_i$ and is suitable for small base excitations and reads

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

(1)

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

(2)

Here, $c_{\mathrm{ijkl}}^E$ are the components of the elasticity tensor at constant electric field and ${\epsilon }_{\mathrm{ij}}^S$ is the dielectric constant at constant strain. The mechanical and electrical fields are coupled via the piezoelectric coupling tensor $e_{\mathrm{kij}}$ .

For large base excitations nonlinear elasticity becomes important and must be considered in the simulations (Stanton et al., 2012). In Stanton et al. (2012) the additional parameters $c_4$ and $c_6$ are postulated for a 1D material model for PZT-5A to include nonlinear elastic behavior via a polynomial function. In Hegendörfer et al. (2022b) this nonlinear elasticity model is extended to enable 3D FEM simulations and yields

$$\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}$$

(3)

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. The linear relation for the dielectric displacement $D_i$ (equation (2)) is maintained. These constitutive equations are introduced into the balance of linear momentum (for the mechanical problem) and Gauss’ law (for the electro-static problem) and solved with the finite element method.

Ultimately, after applying the FEM for the piezoelectric problem and the Bossak-Newmark method for direct time integration, the (nonlinear) vector-valued equation to be solved results as

$${\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$$

(4)

whereby ${\mathit{r}}_{n+1}$ is the residuum at time $t_{n+1}$ , $n$ corresponds to the current time step, while $n+1$ denotes the next time step to be computed and $\alpha$ is the Bossak-Newmark parameter. The force vectors contain the mechanical and electrical fields, whereby ${\mathit{f}}^{\mathrm{dyn}}$ is the inertial force vector, ${\mathit{f}}^{\mathrm{int}}$ is the vector of internal forces, ${\mathit{f}}^{\mathrm{damp}}$ is the vector of damping forces, and ${\mathit{f}}^{\mathrm{ext}}$ is the vector of external forces. In each time step the system of equation (4) is solved with the Newton-Raphson method.

The coupling between the FEM and ECS simulations is realized via the boundary conditions on the electrode of the piezoelectric structure, the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ and the electrode voltage ${\phi }_{\mathrm{el}}$ . In the FEM simulation, two types of boundary conditions are possible:

• The electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ , that is flowing from the electromechanical structure into the electric circuit, can be prescribed in the FEM. Depending on the prescribed electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ , an electric voltage ${\phi }_{\mathrm{el}}$ is generated on the electrode (which is computed by solving equation (4)). Therefore, this approach is denoted as current controlled coupling.

• The electric voltage ${\phi }_{\mathrm{el}}$ can be prescribed in the FEM simulation. Depending on the prescribed electric voltage ${\phi }_{\mathrm{el}}$ an electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ , is generated (again computed by means of equation (4)). Because the electric voltage ${\phi }_{\mathrm{el}}$ is prescribed, this approach is denoted as voltage controlled coupling.

The current and voltage controlled approaches are used to couple the FEM simulations with an ECS.

3. Simulation of the electric circuit (ECS)

An ECS is a powerful numerical tool for analyzing general electric circuits. In the following, the principle procedure of an ECS for the transient simulation of an electric circuit is presented.

First of all, the ECS automatically formulates the circuit equations for example by means of a modified nodal approach (Ho et al., 1975), which is based on Kirchhoff’s circuit laws. The voltages in the electric circuit nodes and the currents in the branches with voltage sources are considered as unknowns to be determined. Independently of which technique is used in the ECS to obtain the circuit equations, the result is in general a system of nonlinear first order differential algebraic equations of the form

$$\mathcal{R}(\mathit{x}\left(t\right),\stackrel{·}{\mathit{x}}\left(t\right),t)=0$$

(5)

representing the behavior of the coupled elements of an electric circuit (Nichols et al., 1994). Here, $\mathit{x}\left(t\right)$ contains the electric circuit variables for example node voltages and $\mathcal{R}$ is a nonlinear vector function. Subsequently, equation (5) is discretized in time and numerical integration is applied (Nichols et al., 1994). As a result, for each time step a time independent system of nonlinear equations is obtained and solved for example with the Newton-Raphson method (Nichols et al., 1994). After obtaining the solution at one time step the ECS advances to the next time step to be computed.

4. Implicit FEM-ECS coupling method

The electromechanical structure is coupled to the electric circuit via the electrode. A real PVEH has a unique electrode voltage ${\phi }_{\mathrm{el}}$ and a unique electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ leaving the electrode. In the presented approach, the FEM simulation of the electromechanical structure and the ECS analysis of the electric circuit are separated. To correctly model the behavior of a “real” PVEH, the electrode voltage ${\phi }_{\mathrm{el}}$ and the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ must agree in the two simulation domains at each time step. This will be realized here, by means of an implicit coupling method. In general, the main difference between an implicit coupling and an explicit coupling is, that an implicit coupling ensures equilibrium between the simulation domains at each time step. This is not the case in explicit coupling. In the following, an implicit coupling method between a 3D coupled electro-mechanical FEM simulation and an ECS analysis is presented, which is either current or voltage controlled. Both implicit coupling approaches render the same results and, thus, can be used equivalently.

The relevant variables for the implicit FEM-ECS coupling are the electrode voltage ${\phi }_{\mathrm{el}}$ and the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ . Since these quantities are scalar values, the coupling numerically corresponds to a coupling at one point.

In the presented approach, the FEM and ECS simulations run independently up to their equilibrium states in each time step with their respective prescribed boundary conditions/input variables. Next, the agreement of the electrode voltage ${\phi }_{\mathrm{el}}$ and of the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ is checked between the simulation domains. If necessary, the prescribed boundary conditions/input variables of the simulations are adjusted according to the implicit coupling algorithm to achieve equilibrium between the FEM and ECS simulations. Since the FEM and ECS simulations run independently at each time step, the algorithm can thus be applied to implicitly couple any FEM and ECS software and is minimally invasive.

Owing to the implicit character of the coupling method, equilibrium is established between the FEM and ECS simulations at time $t_{n+1}$ . For the sake of simplicity, the index $n+1$ , which specifies the next time step to be computed, is omitted in the following. Firstly, in the description of the coupling method, it is assumed that the time steps in both simulation domains are of equal size, so that the simulations of the FEM and the ECS are synchronized at all time points. Then, sub-cycling is discussed, where different time step sizes can be used in the FEM and the ECS simulation.

4.1. Current controlled coupling without sub-cycling

The current controlled coupling considers the electrode voltage ${\phi }_{\mathrm{el}}$ as a function of the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ at the next time step to be computed $t_{n+1}$ , such that ${\hat{\phi }}_{\mathrm{el}}={\hat{\phi }}_{\mathrm{el}}\left({\stackrel{·}{Q}}_{\mathrm{el}}\right)$ . The behavior of the electromechanical structure is simulated with the FEM, resulting in a nonlinear function ${\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}$ for the electric electrode voltage. This function depends on the electric current leaving the electrode ${\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}$ , such that

$${\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}={\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}\left({\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}\right)$$

(6)

The electric circuit is simulated with an ECS, resulting in a relation for the electric voltage at the same electrode ${\hat{\phi }}_{\mathrm{el}}^{\mathrm{ECS}}$ as a nonlinear function of the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$ , so that

$${\hat{\phi }}_{\mathrm{el}}^{\mathrm{ECS}}={\hat{\phi }}_{\mathrm{el}}^{\mathrm{ECS}}\left({\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}\right)$$

(7)

Equilibrium between the FEM simulation and the ECS simulation is achieved when

$${\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}\left({\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}\right)-{\hat{\phi }}_{\mathrm{el}}^{\mathrm{ECS}}\left({\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}\right)=0$$

(8)

The coupling method ensures that ${\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}={\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$ holds for the converged solution (since there is also a unique electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ in a real PVEH).

As an example, Figure 1 illustrates functions (6) and (7) for a linear electromechanical structure ${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{FEM}}$ and an electric circuit consisting only of a load resistance ${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{ECS}}$ . The open circuit voltage OC is reached for the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}=0$ , while closed circuit behavior CC occurs for the electric voltage ${\phi }_{\mathrm{el}}=0$ . Since the behavior of the electromechanical structure is linear, the function ${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{FEM}}$ is linear with a constant slope $K_{C,\mathrm{lin}}$ . The relation for the electrode voltage ${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{FEM}}$ of the electromechanical structure is therefore given as

$${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{FEM}}\left({\stackrel{·}{Q}}_{\mathrm{el}}\right)=\mathrm{OC}+K_{C,\mathrm{lin}}{\stackrel{·}{Q}}_{\mathrm{el}}$$

(9)

Figure 1.

Coupling of a linear electromechanical structure to a load resistor $R$ .

The relation for the electrode voltage ${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{ECS}}$ resulting from the electric circuit follows Ohm’s law, which states that

$${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{ECS}}\left({\stackrel{·}{Q}}_{\mathrm{el}}\right)=R{\stackrel{·}{Q}}_{\mathrm{el}}$$

(10)

and also leads to a straight line in Figure 1.

The searched solution ${\mathcal{S}}_{\mathrm{lin}}$ is represented by the intersection between the two relations ${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{FEM}}$ and ${\hat{\phi }}_{\mathrm{el},\mathrm{lin}}^{\mathrm{ECS}}$ . Both relations share the same electrode voltage ${\phi }_{\mathrm{el}}$ and electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ values at their intersection. Accordingly, the coordinates of the solution ${\mathcal{S}}_{\mathrm{lin}}$ are given as

$${\stackrel{·}{Q}}_{\mathrm{el},\mathrm{lin}}^{*}=\frac{\mathrm{OC}}{R-K_{C,\mathrm{lin}}}\mathrm{and}{\phi }_{\mathrm{el},\mathrm{lin}}^{*}=\frac{R\mathrm{OC}}{R-K_{C,\mathrm{lin}}}$$

(11)

The larger the load resistor, the lower the electric current ${\stackrel{·}{Q}}_{\mathrm{el},\mathrm{lin}}^{*}$ and the higher the electrode voltage ${\phi }_{\mathrm{el},\mathrm{lin}}^{*}$ .

In general, the relations for the electromechanical structure and the electric circuit can be nonlinear. Therefore, a numerical approach is necessary to find the solution $\mathcal{S}$ for general PVEHs.

The Newton-Raphson method is used to derive the electric voltage on the electrode ${\phi }_{\mathrm{el}}$ by solving equation (8). This procedure is schematically shown in Figure 2. The orange curve represents an exemplary nonlinear electric behavior of the electromechanical structure ${\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}$ , while the blue curve represents exemplarily the behavior of the electric circuit. For the particular time step $n+1$ , the intersection point of the two functions ${\widehat{\phi }}_{el}^{FEM}$ and ${\hat{\phi }}_{\mathrm{el}}^{\mathrm{ECS}}$ is the searched solution $\mathcal{S}$ . In the point $\mathcal{S}$ , the electrode voltage ${\phi }_{\mathrm{el}}^{\mathrm{FEM}}$ and the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}$ in the FEM simulation are in equilibrium with the electrode voltage ${\phi }_{\mathrm{el}}^{\mathrm{ECS}}$ and the electric current ${\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$ in the ECS analysis and, thus, equation (8) is fulfilled.

Figure 2.

Newton-Raphson convergence algorithm between the FEM and ECS simulations for the current controlled coupling.

Starting point of the Newton-Raphson method in Figure 2 is the input for the FEM simulation, which is the electric current value ${}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}$ , whereby $(k)$ is the iteration index. This means, that the electric current ${}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}$ is prescribed in the FEM simulation as a boundary condition. As output of the FEM simulation the corresponding electric voltage value ${}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}$ is obtained. For the first iteration the electric current ${}^{(1)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}$ is chosen to be equal to the converged solution from the last time step ${\stackrel{·}{Q}}_{\mathrm{el},n}$ . The linearization of ${\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}$ at ${}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}$ renders the linear approximation of the electric behavior of the electromechanical structure

$${}^{(k)}{\hat{p}}_{\mathrm{el}}{=}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}+\underset{{}^{(k)}{K}_C}{\underbrace{\frac{{}^{(k)}\partial {\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}}{\partial {\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}}}}[{\stackrel{·}{Q}}_{\mathrm{el}}-{}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}]$$

(12)

The slope ${}^{(k)}{K}_C$ is numerically computed via the difference quotient, which requires one additional FEM simulation of the current time increment. Equation (12) can be rewritten as

$${}^{(k)}{\hat{p}}_{\mathrm{el}}=\underset{{}^{(k)}{V}_{\mathrm{DC}}}{\underbrace{{}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}{-}^{(k)}{K_C}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}}}{+}^{(k)}{K}_C{\stackrel{·}{Q}}_{\mathrm{el}}$$

(13)

whereby ${}^{(k)}{V}_{\mathrm{DC}}$ is the y-axis intercept.

The linear approximation of the electric behavior of the electromechanical structure ${}^{(k)}{\hat{p}}_{\mathrm{el}}$ has to be modeled in the ECS in each iteration of the implicit coupling algorithm. To this end, a controlled voltage source is used, whereby its output voltage ${\hat{\phi }}_{\mathrm{el}}^{\mathrm{ECS}}(t)$ is a superposition of the time dependent function blocks ${}^{(k)}{v}_{\mathrm{DC}}(t)$ and ${}^{(k)}{k}_C(t)$ . Figure 3 presents the implementation of the controlled voltage source in the ECS. The electric current ${}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}(t)$ is multiplied by ${}^{(k)}{k}_C(t)$ and added to ${}^{(k)}{v}_{\mathrm{DC}}(t)$ , which yields the output voltage ${}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{ECS}}(t)$ .

Figure 3.

Implementation of the controlled voltage source for the current controlled coupling in the ECS.

For consistency between the time steps, the ECS has to use correct initial values at time $t_n$ to compute the current values at time $t_{n+1}$ . These initial values have to be derived from the converged values of the previously computed time step. The y-axis intercept $V_{\mathrm{DC},n}$ and the slope $K_{C,n}$ of the previous time step $n$ must be provided as

$${}^{(k)}{k}_C(t_n)=K_{C,n}$$

(14)

$${}^{(k)}{v}_{\mathrm{DC}}(t_n)=V_{\mathrm{DC},n}$$

(15)

Equations (14) and (15) define the initial conditions for each iteration to ensure consistency. These are used in the current iteration to compute the values at time $t_{n+1}$ , using the y-axis intercept and slope derived from the FEM simulation

$${}^{(k)}{k}_C(t_{n+1}){=}^{(k)}{K}_C$$

(16)

$${}^{(k)}{v}_{\mathrm{DC}}(t_{n+1}){=}^{(k)}{V}_{\mathrm{DC}}$$

(17)

The ECS results in the solution ${}^{(k)}\mathcal{S}$ for the controlled voltage source with coordinates ( ${}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$ , ${}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{ECS}}$ ), as illustrated in Figure 2. Convergence at the interface between the FEM and ECS simulations is thus obtained if

$$\underset{{}^{(k+1)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}}{\underbrace{{\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}\left({}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}\right)}}{-}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{ECS}}\approx 0$$

(18)

In equation (18), the evaluation of the function ${\hat{\phi }}_{\mathrm{el}}^{\mathrm{FEM}}$ at the electric current ${}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$ is required. This evaluation necessitates one FEM simulation, but at the same time represents the starting point for the next iteration, namely ${}^{(k+1)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}$ . If there is no convergence at the interface, the updated FEM electric current is given as

$${}^{(k+1)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{FEM}}{=}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$$

(19)

and the algorithm is repeated until convergence is obtained. In case of convergence, the electric current at the FEM-ECS interface is given as

$${\stackrel{·}{Q}}_{\mathrm{el}}{=}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$$

(20)

Finally, the complete states of the FEM (e.g. nodal displacements) and the ECS (e.g. voltage values of the circuit elements) simulations are saved. Accordingly, these states are loaded as initial conditions for the computation of the next timestep.

Figure 4 summarizes the algorithm for the current controlled coupling. Starting with the electric current of the previous time step ${\stackrel{·}{Q}}_{\mathrm{el},n}$ , the corresponding voltage value ${}^{(1)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}$ is computed, which requires an FEM simulation. Using this result, an additional FEM simulation is then required within the iteration loop to obtain the linearized electric behavior of the electromechanical structure. Subsequently, the solution of the current iteration ${}^{(k)}\mathcal{S}=\left({}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}{,}^{(k)}{\phi }_{\mathrm{el}}^{\mathrm{ECS}}\right)$ is computed by means of the ECS. An additional FEM simulation, using the ECS result ${}^{(k)}{\stackrel{·}{Q}}_{\mathrm{el}}^{\mathrm{ECS}}$ as a boundary condition, computes the new electrode voltage ${}^{(k+1)}{\phi }_{\mathrm{el}}^{\mathrm{FEM}}$ , which is required for the convergence check and also represents the starting point for the next iteration. The procedure is repeated until convergence is achieved at the FEM-ECS interface. Accordingly, for $k$ iterations of the coupling method, $2k+1$ FEM simulations are necessary while only $k$ ECS simulations are required. The Newton-Raphson method which is applied for the implicit coupling approach shows quadratic convergence behavior. If a linear electromechanical structure is coupled to an arbitrary nonlinear electric circuit, convergence is achieved in the first iteration.

Figure 4.

Summary of the algorithm for the current controlled coupling.

4.2. Voltage controlled coupling without sub-cycling

The voltage controlled coupling scheme is quite similar to the current controlled coupling. The electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ is considered as a function of the electric electrode voltage ${\phi }_{\mathrm{el}}$ , such that ${\stackrel{·}{\hat{Q}}}_{\mathrm{el}}={\stackrel{·}{\hat{Q}}}_{\mathrm{el}}\left({\phi }_{\mathrm{el}}\right)$ . The electric behavior of the electromechanical structure ${\stackrel{·}{\hat{Q}}}_{\mathrm{el}}^{\mathrm{FEM}}\left({\phi }_{\mathrm{el}}^{\mathrm{FEM}}\right)$ is modeled with the FEM and the behavior of the electric circuit ${\stackrel{·}{\hat{Q}}}_{\mathrm{el}}^{\mathrm{ECS}}({\phi }_{\mathrm{el}}^{\mathrm{ECS}})$ is modeled with an ECS. Convergence between the simulations is obtained when

$${\stackrel{·}{\hat{Q}}}_{\mathrm{el}}^{\mathrm{FEM}}\left({\phi }_{\mathrm{el}}^{\mathrm{FEM}}\right)-{\stackrel{·}{\hat{Q}}}_{\mathrm{el}}^{\mathrm{ECS}}\left({\phi }_{\mathrm{el}}^{\mathrm{ECS}}\right)=0$$

(21)

In the FEM simulation, the voltage controlled coupling scheme results in Dirichlet boundary conditions for the electric voltage ${\phi }_{\mathrm{el}}$ at the electrode. The electric current ${\stackrel{·}{Q}}_{\mathrm{el}}$ can then be computed in a post-processing step.

In the ECS, the linearized behavior of the electromechanical structure is modeled as a controlled current source, consisting of two time dependent function blocks.

The detailed equations and illustrations of the voltage controlled coupling are given in Appendix A.

4.3. Implicit coupling with sub-cycling

The presented implicit coupling method allows to choose different time step sizes for the FEM and ECS simulations, that is sub-cycling of one method is possible. Then, the coupling takes place at the synchronization points and the time steps for FEM and ECS between these points can be selected independently. Typically, the electrical frequencies in energy harvesting applications are orders of magnitude higher than the mechanical frequencies. Thus, smaller time step sizes are required in the ECS than in the FEM simulations, and sub-cycling of the ECS is desired. The FEM time steps are considered as synchronization points between the two simulation domains and sub-cycling is only applied in the ECS simulation.

In this case intermediate values for the tangents and y-axis intercepts in the interval $[t_n,t_{n+1}]$ must be specified. Their values at the discrete time points $t_n$ and $t_{n+1}$ are specified in equations (14)–(17) for the current controlled coupling and in equations (25)–(28) for the voltage controlled coupling. For the sub-cycles, linear interpolation between these values is applied as

$$X(t)=\frac{X_n[t-t_{n+1}]}{t_n-t_{n+1}}+\frac{X_{n+1}[t-t_n]}{t_{n+1}-t_n}$$

(22)

whereby $X$ denotes the tangents or y-axis intercepts. Equation (22) defines continuous values for the tangents and y-axis intercepts within the time interval $[t_n,t_{n+1}]$ , which are integrated in the ECS simulations. This enables the ECS to sub-cycle and to use powerful solvers with automatically adapted time step sizes that solve system (5) and ensure a certain accuracy and stability.

5. Validation of the implicit coupling

The FEM simulation of the electromechanical structure is implemented based on the open source C++ Finite Element library deal.ii targeted at the efficient solution of partial differential equations (Arndt et al., 2019). This FEM library can perform all steps of electromechanical coupled FEM simulations, including meshing, assembling the system matrix, solving linear systems of equations and post-processing the results, and is used to solve equation (4) for the electromechanical structure. For the simulation of the electric circuit the multiphysics simulator MATLAB/Simulink Simscape is applied, which allows for the simulation of arbitrary electric circuits, for the development of control systems, and, for example, for design optimization. To ensure computational efficiency, the iterative characteristic of the presented implicit coupling method requires that the FEM and ECS software can store and load their internal states to update the initial conditions after the computation of each timestep (Busch, 2012). The current state of the electric circuit simulation can be conveniently saved and loaded in MATLAB/Simulink Simscape and the in-house FEM code also offers this option.

To validate the presented current and voltage controlled FEM-ECS coupling and to demonstrate its capabilities, a unimorph PVEH with different nonlinear electric circuits is simulated. The results of the implicit coupling method are compared against results reported in the literature obtained by explicitly and monolithically coupled FEM-ECS simulations. Moreover, the computational efficiency of the different coupling methods is discussed.

5.1. Unimorph PVEH with half diode bridge circuit

Figure 5 shows a unimorph electromechanical structure introduced in Erturk and Inman (2008) consisting of a layer of substructure and a layer of PZT-5A. The PVEH is excited at its base with a harmonic base acceleration $a(t)$ with a magnitude of $1\mathrm{m}/{\mathrm{s}}^2$ and a frequency of $48\mathrm{Hz}$ .

Figure 5.

Unimorph PVEH with half diode bridge circuit.

The electric circuit used is a half diode bridge circuit consisting of two diodes $D_1$ and $D_2$ , a capacitor $C_s$ , and a load resistor $R_d$ . Table 1 presents the material properties of the substructure and geometrical and electrical parameters, while the material data for PZT-5A is given in Appendix B. This application example was originally developed in Elvin and Elvin (2009a).

Table 1.

Material properties of the substructure and geometrical parameters of the unimorph PVEH from Erturk and Inman (2008). 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 $[\mathrm{kg}/{\mathrm{m}}^3]$	7165
Diodes $D_1$ and $D_2$	1N914
$R_d$ [ $\mathrm{M}\Omega$ ]	$10$
$C_s$ [ $\mu \mathrm{F}$ ]	$3.3$

For the FEM simulation the electromechanical structure is discretized with 120 quadratic hexahedral elements. The same quadratic ansatz functions are used to approximate the displacements and the electric voltage. For the non-piezoelectric substructure, elements with only mechanical degrees of freedom but without electrical degrees of freedom are used. A constant time step size of $0.1\mathrm{ms}$ is used in the FEM simulation, while sub-cycling of the ECS is enabled. In the ECS analysis, the SPICE diode model available in Matlab/Simulink Simscape is used to model the behavior of the diodes $D_1$ and $D_2$ . The FEM time steps are defined as synchronization points, so that equilibrium is established between the FEM and ECS simulations in all FEM time steps.

Figure 6 compares the results of the current and voltage controlled coupling approaches against the results of the explicitly coupled FEM-ECS simulation from Elvin and Elvin (2009a). The differences are mainly due to oscillations of the electric voltage ${\phi }_{\mathrm{el}}$ in Elvin and Elvin (2009a), which could be caused by the explicit coupling. Furthermore, in Elvin and Elvin (2009a) thin shell elements are used for the simulation of the mechanical behavior of the structure and the electric behavior is simplified to a lumped capacitance. In this contribution, 3D continuum elements are applied for the simulation of the structure, precisely taking into account the coupled electromechanical behavior. Overall, the implicit current and voltage controlled coupling approaches render the same results and a good agreement with the simulation results reported in Elvin and Elvin (2009a) is observed.

Figure 6.

Comparison of the piezoelectric voltage ${\phi }_{\mathrm{el}}$ of the current and voltage controlled coupling methods against the results in Elvin and Elvin (2009a) for the unimorph PVEH with half diode bridge circuit.

Furthermore, in Elvin and Elvin (2009a) the accuracy of the explicit coupling method is studied for different numbers of time steps per oscillation period $N$ . To this end, the unimorph PVEH with a simple load resistance $(1\mathrm{M}\Omega )$ as electric circuit is considered and a harmonic base acceleration with a frequency of 48.8 Hz is applied. After $1\mathrm{s}$ the relative error in the maximum voltage between a monolithic simulation and the explicit coupling is evaluated and is shown in Figure 7. For $N=10$ the relative error is around $3$ 3% and decreases for an increasing number of time steps per oscillation period $N$ .

Figure 7.

Relative error in the maximum voltage between the explicit coupling in Elvin and Elvin (2009a) and a monolithic simulation and between the proposed implicit voltage and current controlled couplings and a monolithic simulation.

Similarly to Elvin and Elvin (2009a), the accuracy of the proposed implicit coupling method is evaluated by comparing the implicitly coupled FEM-ECS simulation results against monolithic FEM simulations. The load resistance is considered in the monolithic FEM simulations via the vector of external forces, similar to the method presented in Hegendörfer et al. (2022b). Different numbers of time steps per oscillation period $N$ are tested and for each time step size chosen, the relative error in the maximum voltage between the implicitly and monolithically coupled simulations vanishes, as shown in Figure 7. Consequently, consistently with Matthies et al. (2006), the accuracy of the implicit coupling method is equivalent to a monolithical coupling. As a result, the implicit coupling allows for larger time step sizes than an explicit coupling.

In the implementation of the explicit coupling method in Elvin and Elvin (2009a), the entire time history must be recalculated at each time step, since it is not possible to restart the simulations from a previous time step. This fact leads to a significant unnecessary additional computational effort. In contrast, our implementation using the implicit coupling method can store and load the states of the FEM and ECS simulations for the calculation of each time step to overcome this problem.

It should be noted, that convergence between the FEM and ECS simulations of the unimorph PVEH is achieved in the first Newton-Raphson iteration for both electric circuits, the half diode bridge and the load resistor. Although the half diode bridge circuit is strongly nonlinear, only one iteration is necessary because the electromechanical structure is linear.

5.2. Unimorph PVEH with SSHI circuit

Figure 8 shows the linear unimorph electromechanical structure from the previous example, which is here equipped with an SSHI circuit. A constant voltage $V_{\mathrm{const}}=1.8\mathrm{V}$ is assumed to emulate an ideal energy storage device and a drop voltage of the diodes $V_{\mathrm{drop}}=0.6\mathrm{V}$ is specified. The PVEH is harmonically excited at its base with a magnitude of $1\mathrm{m}/\mathrm{s}$ and a frequency of $48.7\mathrm{Hz}$ , which is close to the first eigenfrequency of the piezoelectric structure. This application example was originally developed in Gedeon and Rupitsch (2018).

Figure 8.

Unimorph PVEH with SSHI circuit.

For details about the SSHI circuit please refer to the literature (e.g. Guyomar et al., 2005; Hegendörfer et al., 2022b; Shu et al., 2007). A self-adaptive switching control for the SSHI circuit is used, which triggers the switch when the piezoelectric voltage $|{\phi }_{\mathrm{el}}|$ decreases with the following restrictions: Due to the instationary settling process, when the base excitation starts, triggering the switch is prohibited during the first $t_{\mathrm{switch}}=10.8\mathrm{ms}$ . Moreover, 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}}=4\mathrm{ms}$ is introduced after the switch has been opened, during which triggering the switch is prohibited.

In order to enable efficient and accurate simulations, variable time step sizes are introduced in the FEM simulation: For the process of voltage inversion, during which the switch is closed, a time step size of $2·{10}^{-4}\mathrm{ms}$ is applied since it happens almost instantaneously. Subsequently, after the switch is opened again, the time step size is set to $5·{10}^{-3}\mathrm{ms}$ for a time span $t_{\mathrm{osc}}=0.25\mathrm{ms}$ to capture the higher frequency oscillations. For the remaining time period a time step size of ${10}^{-1}\mathrm{ms}$ is used. Sub-cycling of the ECS is enabled and all FEM time steps are selected as synchronization points. For the simulation of the electric circuit in Matlab/Simulink Simscape the implicit variable step solver daessc is chosen. The maximum and minimum step sizes of the daessc solver are set to auto, which means that Simscape determines these parameters automatically. Normally Simscape uses about 60 time steps during the sub-cycling between synchronization points. However, if the time period between two synchronization points changes significantly, which is the case when the switch is triggered, the ECS needs hundreds of time steps. If no sub-cycling of the ECS is allowed, this automatic adaption of the ECS time steps is not possible, which leads to numerical instabilities of the electric circuit simulation and oscillations of the electric voltage ${\phi }_{\mathrm{el}}$ . In contrast, when sub-cycling is enabled, no stability problems occur.

Figure 9 compares the results for the electrode voltage ${\phi }_{\mathrm{el}}$ of the current and voltage controlled implicitly coupled simulations against the result of a reduced and monolithically coupled simulation using Matlab/Simulink in Gedeon and Rupitsch (2018). Because of the abrupt inversion of the piezoelectric voltage when the switch of the SSHI circuit is triggered, higher frequency modes are excited, which are shown in the zoom-in in Figure 9. Within the considered time period of 80 ms the largest differences between the results of the presented implicit coupling and the approach from Gedeon and Rupitsch (2018) occur during the instationary settling process. After a time span of about $45\mathrm{ms}$ a very good agreement of the two models can be observed.

Figure 9.

Comparison of the piezoelectric voltage ${\phi }_{\mathrm{el}}$ of the current and voltage controlled coupling methods against the results in Gedeon and Rupitsch (2018) for the unimorph PVEH with SSHI circuit.

The results of the current and voltage controlled couplings are almost identical, but there are very small differences due to tolerances of the solver and the numerical evaluation of the slopes $K_V$ and $K_C$ . The results do not change when different maximum and minimum time step sizes are manually specified in the ECS.

The monolithic FEM-ECS couplings reported in the literature (Aftab et al., 2012; Gedeon and Rupitsch, 2018; Xiang et al., 2018) apply a model order reduction technique to the FEM model. Therefore, the monolithic coupling methods are computationally more efficient than the implicit coupling method using the full FEM model. However, the implicit coupling method can exploit the full capabilities of the FEM for the consideration of nonlinear electromechanical structures. The reduced order models in contrast are limited to linear piezoelectric structures.

Considering the different PVEH simulations, their results with the current and voltage controlled coupling methods agree with the results reported in the literature. It is confirmed that the presented system simulation method with an implicit coupling of FEM and ECS gives accurate results. The presented implicit coupling method can be applied to all FEM and ECS tools that can read and export time step results.

6. Application example

Figure 10 on the left shows a bimorph piezoelectric structure introduced in Erturk and Inman (2009) with a tip mass $m_t$ . The cantilevered bimorph electromechanical structure consists of two layers of PZT-5A (material parameters in Appendix B) bracketing a passive substructure layer. The piezoceramic layers are poled in opposite directions. Later, a large base excitation is applied and therefore the nonlinear elastic behavior of PZT-5A has to be taken into account (Stanton et al., 2012). Neglecting the nonlinear elastic behavior of PZT-5A for large base excitations leads to a significant overestimation of the harvested energy in the simulations (Hegendörfer et al., 2022b). To consider nonlinear elasticity, the nonlinear constitutive law for PZT-5A, equations (2) and (3), is applied. The parameters of the bimorph electromechanical structure are provided in Table 2.

Figure 10.

Application example consisting of a nonlinear bimorph electromechanical structure and a self-powered SSHI circuit.

Table 2.

Parameters of the bimorph electromechanical structure from Erturk and Inman (2009).

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 $[\mathrm{kg}/{\mathrm{m}}^3]$	9000
Tip mass $m_t$ [kg]	0.012

The nonlinear bimorph electromechanical structure is coupled with an efficient self-powered SSHI circuit presented in Eltamaly and Addoweesh (2017). This SSHI circuit consists of diodes and transistors and is shown in Figure 10 together with the nonlinear bimorph electromechanical structure. One advantage of this SSHI circuit is that only the two capacitors $C_2$ and $C_3$ are used to detect maxima and minima of the electrode voltage ${\phi }_{\mathrm{el}}$ . The SSHI circuit is used to charge the storage capacitor $C_r$ and its working principle is explained in detail in Eltamaly and Addoweesh (2017).

In the MATLAB/Simulink Simscape simulations the same models are available for the semiconductor components as used in the circuit simulator SPICE. For the simulation of the transistors $Q_1$ , $Q_2$ , $Q_3$ , and $Q_4$ the Gummel-Poon model, which is presented in Gummel and Poon (1970), is solved and the SPICE model for the diodes $D_1$ , $D_2$ , $D_3$ , and $D_4$ of the full-bridge rectifier is used. The parameters and specifications of the electrical components used for the self-powered SSHI circuit are shown in Table 3, while different values for the storage capacitor $C_r$ and inductor $L_1$ with parasitic resistance $R_1$ are applied in the simulations.

Table 3.

Electric circuit parameters and specifications of the semiconductor components used for the self-powered SSHI circuit.

Capacitor $C_1$ [nF]	2
Capacitor $C_2$ [pF]	1
Diodes $D_1$ , $D_2$ , $D_3$ , and $D_4$	1N914
PNP transistors $Q_1$ and $Q_3$	2N2907
NPN transistors $Q_2$ and $Q_4$	2N2222

Since the simulations are executed in the time domain, discrete time step sizes must be introduced. FEM and ECS simulations must be conducted at each time step, so the time step size significantly affects the accuracy and the computational costs. A smaller time step size means higher computational costs, but the simulations become more accurate. If discontinuities occur in the simulations very small time step sizes are required to resolve them. Conversely, if the states of the simulation do not change rapidly, large time step sizes are desired for computational efficiency. Therefore, efficient simulations in the time domain require an adaptation of the time step sizes to the current conditions in the simulation. Consequently, an algorithm is necessary which adapts the time step sizes automatically to ensure efficient and accurate simulations.

The commercial software Matlab/Simulink Simscape used to simulate the electric circuit provides powerful solvers with automatically adjusted time step sizes. The ECS can run sub-cycles and all solvers of Matlab/Simulink Simscape are available. As a result, the electric circuit can be solved efficiently.

For the adaptive control of the applied time step sizes in the FEM simulation, a purely empirical method, presented in Hibbitt and Karlsson (1979), which is based on limiting a residual calculated between two FEM time steps, is used. Furthermore, to accurately capture the interaction between the FEM and ECS simulations the maximum change of the absolute value of the electrode voltage $|{\phi }_{\mathrm{el}}|$ between two consecutive FEM time steps is limited to a value ${\phi }_{\mathrm{el},\max }$ . Maximum $t_{\max }$ and minimum $t_{\min }$ time step sizes are introduced and the time step size can be adjusted between these limits.

In the next step, the ECS simulation results for the SSHI circuit are compared against results reported in the literature to validate the ECS analysis. Finally, a shock-like base excitation is applied to the piezoelectric structure and the results of the coupled FEM-ECS simulation are discussed.

6.1. Validation of the ECS simulation

A similar self-powered SSHI circuit as presented in Figure 10 is also considered in Zhang et al. (2021), coupled to a linear electromechanical structure. In Zhang et al. (2021), the response of the linear piezoelectric structure for the first vibration mode is represented by means of an equivalent circuit model, to allow the simulation of the PVEH in an ECS. To validate our Matlab/Simulink Simscape model of the SSHI circuit the same equivalent circuit model of the electromechancal structure is simulated together with the SSHI circuit. The same storage capacitance $C_r=100\mu F$ , inductance $L_1=50\mathrm{mH}$ , and resistance $R_1=60\Omega$ values as in Zhang et al. (2021) are applied for the validation, however, different diode models are used. In Zhang et al. (2021) only the forward drop voltage value of the diodes is available, but in this contribution the SPICE model for the diodes is applied.

Figure 11 compares the results for the voltage across the storage capacitor $C_r$ of our Matlab/Simulink Simscape simulation against the results in Zhang et al. (2021). After $60\mathrm{s}$ , in Zhang et al. (2021) a voltage across the storage capacitor $C_r$ of about $8.1\mathrm{V}$ is reached, while the own simulation yields a voltage of $7.8\mathrm{V}$ . This deviation of only around 3.7% means a very good agreement between the two simulations and validates the electric circuit simulation in Matlab/Simulink Simscape.

Figure 11.

Comparison of the present Matlab/Simulink Simscape simulation results for the self-powered SSHI circuit against the results reported in Zhang et al. (2021).

6.2. Shock-like excitation

A shock-like base excitation $a(t)$ with a magnitude of 0.2 mm and a duration of 10 ms is applied, which is shown in Figure 12. The electric circuit values for the capacitor $C_r=1\mu F$ , inductor $L_1=1\mathrm{mH}$ , and resistor $R_1=1\Omega$ are defined. Sub-cycling of the ECS is enabled and the solver ode15s is chosen in Matlab/Simulink Simscape. The bimorph electromechanical structure is discretized with 90 quadratic hexahedral piezoelectric continuum elements. Using an Intel Core i7-9700 CPU the simulation takes about 8 h and the adaptive time stepping algorithm requires around 4300 time steps.

Figure 12.

Shock-like base excitation $a(t)$ with a magnitude of 0.2 mm and a duration of 10 ms.

Due to the piezoelectric effect, the mechanical and electrical behavior of the bimorph electromechanical structure influence each other. Taking into account nonlinear elasticity, the electric voltage ${\phi }_{\mathrm{el}}$ depends nonlinearly on the electric current and vice versa. Therefore, the tangents ${}^{(k)}{K}_C$ in equation (12) and ${}^{(k)}{K}_V$ in equation (23) are not constant anymore, as in the case of linear electromechanical structures. Consequently, some Newton-Raphson iterations are required to achieve convergence between the FEM and ECS simulations. However, for this application example, a maximum of only two iterations for all time steps is required. This confirms the rapid convergence behavior of the presented coupling method.

Figure 13 presents the time signal of the electrode voltage ${\phi }_{\mathrm{el}}$ . At time $t=0$ , the electromechanical structure is at rest and the voltage across the storage capacitor $C_r$ vanishes. Then the shock-like base excitation $a(t)$ starts, and as can be seen in Figure 13, the oscillation of the electromechanical structure consists of the superposition of many vibration modes. As the zoom-in in Figure 13 in the middle shows, these higher frequency modes trigger the action of the SSHI circuit and cause inversions of the electrode voltage ${\phi }_{\mathrm{el}}$ at about $10.6\mathrm{ms}$ and $12.1\mathrm{ms}$ . After about $30\mathrm{ms}$ , the higher frequency modes have largely decayed due to damping. Then, the vibration of the electromechanical structure mainly follows its first vibration mode. The zoom-in in Figure 13 on the right shows the process of voltage inversion at about $96.5\mathrm{ms}$ . Because of the non-ideal behavior of the SSHI circuit, the electrode voltage ${\phi }_{\mathrm{el}}$ must drop by a certain amount of its maximum value to trigger the inversion process of the SSHI circuit. Subsequently, the rapid inversion of the electrode voltage ${\phi }_{\mathrm{el}}$ at about $96.5\mathrm{ms}$ excites higher frequency modes, which are strong enough that they are detected as maximum values by the SSHI circuit. Therefore the electrode voltage ${\phi }_{\mathrm{el}}$ is inappropriately inverted a second time immediately after the first inversion. Since the electrode voltage ${\phi }_{\mathrm{el}}$ continues to evolve in the opposite direction as assumed by the SSHI circuit, a third inversion is triggered at about $96.9\mathrm{ms}$ .

Figure 13.

Time signal of the electrode voltage ${\phi }_{\mathrm{el}}$ for the self-powered SSHI circuit coupled to the nonlinear bimorph electromechanical structure under shock-like base excitation.

During the first around $100\mathrm{ms}$ of the simulation, the maximum values of the electrode voltage ${\phi }_{\mathrm{el}}$ increase as the storage capacitor $C_r$ is charged. Then, because of the damping, the electrode voltage ${\phi }_{\mathrm{el}}$ does not reach sufficiently high values to further charge the storage capacitor $C_r$ and continues to oscillate until it vanishes. The storage capacitor $C_r$ reaches after $100\mathrm{ms}$ a voltage of around $4.5\mathrm{V}$ , which results in a harvested energy of approximately $0.01\mathrm{mJ}$ .

For practical applications, an accurate evaluation of the mechanical stresses in the electromechanical structure is paramount. It depends on them whether the electromechanical structure fails during operation. Considering that PZT-5A is a brittle ceramic material, Rankine’s theory is applied, which states that failure occurs when the maximum principal stress at any point reaches a value equal to the maximum dynamic tensile strength of PZT-5A. An advantage of using the FEM to simulate the electromechanical structure is, that the mechanical stresses in the electromechanical structure can be precisely evaluated. Figure 14 shows the time evolution of the maximum principal stresses in the PZT-5A layers. For times $t<30$  ms, the maximum principal stresses fluctuate due to the higher frequency modes. The maximum value of $16.5\mathrm{MPa}$ for the maximum principal stresses occurs at about $11.3\mathrm{ms}$ . In Berlincourt et al. (2000), the maximum dynamic tensile strength of PZT-5A is reported as $27.6\mathrm{MPa}$ . Using a safety factor of 1.5, a typical value in aerospace engineering, the maximum allowable stress is $18.4\mathrm{MPa}$ . Therefore, it is concluded that the electromechanical structure does not fail during operation.

Figure 14.

Time evolution of the maximum principal stresses in the PZT-5A layers of the electromechanical structure.

The precise consideration of the behavior of the electric circuit and the nonlinear electromechanical structure leads to a realistic simulation of a PVEH. All mechanical and electrical quantities are available in the simulation. As a result, the practicality of the presented implicit coupling method for the simulation of PVEHs is demonstrated.

7. Conclusion

This contribution addresses the problem of the coupled simulation of a PVEH consisting of an electromechanical structure and an electric circuit which influence each other. The system simulation method presented allows the FEM to be applied to the simulation of the electromechanical structure and an ECS to be used for the electric circuit by means of an implicit coupling. The full capabilities of the FEM can be exploited without restrictions and the complete set of 3D coupled piezoelectric equations is solved. Varying mechanical boundary conditions and nonlinearities of the electromechanical structure, which are essential characteristics of PVEHs, can be taken into account. Moreover, using an ECS allows for the consideration of arbitrary and nonlinear electric circuits. In summary, no simplifications, which could limit the applicability of the method, are introduced in the FEM or ECS simulations. Consequently, the proposed system simulation method allows for the simulation of arbitrary PVEHs and combines the advantages of the FEM and of an ECS.

The implicit coupling at the interface between the FEM and ECS simulations is realized with the Newton-Raphson method. In the case of a linear electromechanical structure, the implicit coupling method achieves convergence in the first Newton-Raphson iteration, regardless of the electric circuit. Moreover, in the case of the nonlinear bimorph electromechanical structure, the presented FEM-ECS coupling method requires only a few iterations between the nonlinear FEM and ECS simulations, which confirms the rapid convergence behavior of the presented coupling method. Furthermore, sub-cycling of the ECS ensures numerical stability of the implicit coupling method even if the time period between two synchronization points changes dramatically. In addition, the FEM-ECS coupling is implicit, meaning that significantly larger time steps can be used compared to an explicit coupling.

The presented implicit FEM-ECS coupling method is validated against a unimorph PVEH equipped with different nonlinear electric circuits from literature. Moreover, an application example consisting of a nonlinear electromechanical structure and a self-powered SSHI circuit is considered and a shock-like base excitation is applied. These examples demonstrate the applicability of the presented implicit coupling method for efficient simulations of general nonlinear PVEHs.

The main limitation of the presented implicit coupling method is the high computational cost, which is generally a problem in transient FEM simulations. Taking into account general nonlinear electromechanical behavior requires transient simulations and the associated effort is not always avoidable. The consideration of both complicated electromechanical structures and long time histories leads to long simulation times, which can make the implicit coupling method impractical in some cases.

It should be noted that the presented FEM-ECS coupling method is not limited to energy harvesting applications, but can also be used for simulations of general smart structures and devices.