*2.1. Design Description*

The mechanical resonator design consists of two identical 80 mm long arms (referred to as outer beams), mechanically coupled through a common end to a 60 mm long inner beam, which extends in turn towards the fixed end [34]. The tip masses are of identical weight *m* = 7.6 g. The first two resonance frequencies and mode shapes, as obtained from a harmonic analysis, are shown in Figure 1. The deformation of the piezoelectric layers results in a surface charge distribution and consequently a voltage across the patch electrodes. The maximum power output of the harvester can be derived from this voltage and the patch capacitance as explained in [34]. The phase di fference between the inner and outer patch depends on the frequency of operation. We identified the phase di fference to be 200 and 30◦ at the first, respectively the second mode. A connection of both patches shall consider and adapt the polarity of the voltage output. Realistic vibration excitation might excite both modes simultaneously so that individual power processing is mandatory.

In order to increase the frequency agility of the system, we propose magnetic frequency tuning. The approach uses permanent magnets on the resonating structure together with fixed external magnets. The interaction between the fixed and movable magnets creates an attracting or repelling force, which adds to the mechanical restoring force of the structure. The change in sti ffness of the structure leads to a

frequency up- or down-tuning. In [34] we characterized the resonator and demonstrated bidirectional frequency tuning of 15%. In order to simulate the effect of the magnetic force, the use of the full FE model is not efficient anymore, due to the very long solution time needed for such a simulation. Therefore, in [34] we proposed a compact model of the resonator and we were able to fully simulate the behavior of such a system. However, the integration of the piezoelectric transducer elements is not feasible with a lumped modelling approach. Therefore, in this paper, we derive a reduced order model of the full FE model and implement it in a system-level simulation.

**Figure 1.** (**a**) Geometry description of the harvester design together with (**b**) its simulated displacement, respectively power output. The transfer function illustrates the dual frequency operation of the structure under a base acceleration of 0.5 g. Simultaneously, it demonstrates comparable power output levels.

#### *2.2. Reduced Order Model of the Piezoelectric Energy Harvester*

The piezoelectric energy harvester model implemented in ANSYS® Mechanical (V2019, R1) was thoroughly described in [34]. Considering its high computational cost for a transient simulation, Krylov subspace-based model order reduction (MOR) methods, also known as rational interpolation [38–40] were implemented to generate a highly compact but accurate reduced order model (ROM). Furthermore, based on this ROM, a circuit-device co-simulation of the piezoelectric energy harvester became feasible in the system-level simulation.

## Model Order Reduction

The finite element method provides a spatially discretized mathematical description of the piezoelectric energy harvester. It is represented by a second order multiple-input multiple-output (MIMO) system of the form:

$$\sum\_{N} \underbrace{\left\{ \underbrace{\begin{bmatrix} M\_{11} & 0\\ 0 & 0 \end{bmatrix}}\_{M} \right\}}\_{M} + \underbrace{\left\{ \begin{array}{cccc} E\_{11} & 0\\ 0 & 0 \end{array} \right\}}\_{E} \underbrace{\begin{bmatrix} \dot{x}\_{1} \\ 0 \end{bmatrix}}\_{E} + \underbrace{\left\{ \begin{array}{cccc} K\_{11} & K\_{12} \\ K\_{21} & K\_{22} \end{array} \right\}}\_{x} \underbrace{\begin{array}{c} x\_{1} \\ x\_{2} \end{array}}\_{x} \right) = \underbrace{\begin{bmatrix} B\_{1} \\ B\_{2} \end{bmatrix}}\_{B} \mu \tag{1}$$

where *M*, *E*,*K* ∈ R*N*×*<sup>N</sup>* are the mass, damping, and stiffness matrices, respectively. *B* ∈ R*N*×*p* and *C* ∈ R*q*×*<sup>N</sup>* are the input and gathering matrices, with *u* ∈ R*p* and *y* ∈ C*q* being user defined input and output vectors. *x*1 and *x*2 present the nodal displacement and electrical potentials in the state vector *x* ∈ C*N*. Considering the large size of system (1) obtained in this work (*N* = 166,789 DoF), a highly compact but accurate ROM was generated:

$$\sum\_{r} \begin{cases} \underbrace{\boldsymbol{V}^{T} \boldsymbol{M} \boldsymbol{V} \cdot \dot{\boldsymbol{z}}}\_{\boldsymbol{M}\_{r}} + \underbrace{\boldsymbol{V}^{T} \boldsymbol{E} \boldsymbol{V} \cdot \dot{\boldsymbol{z}}}\_{\boldsymbol{E}\_{r}} + \underbrace{\boldsymbol{V}^{T} \boldsymbol{K} \boldsymbol{V} \cdot \boldsymbol{z}}\_{\boldsymbol{K}\_{r}} = \underbrace{\boldsymbol{V}^{T} \boldsymbol{B} \cdot \boldsymbol{u}}\_{\boldsymbol{B}\_{r}} \\ \boldsymbol{y} = \underbrace{\boldsymbol{\mathcal{C}} \boldsymbol{V} \cdot \boldsymbol{z}}\_{\boldsymbol{C}\_{r}} \end{cases} \tag{2}$$

where *V* ∈ R*N*×*r* is the orthonormal basis of the second order Krylov subspace <sup>K</sup>*r*−*K*−1*E*,<sup>−</sup>*K*−1*M*,<sup>−</sup>*K*−1*B* obtained through the Second Order Arnoldi Reduction (SOAR) method [41,42]. The full-scale state vector *x* is approximated by *x* ≈ *V*·*z*, where *z* ∈ C*r* is the reduced state vector,*r* = 35 *N*.

In the previous research [43–45], the stability of the reduced system (2) could not be guaranteed by the conventional MOR methods. Several stabilization approaches have been introduced and mathematically proven in [46–48]. In this work, we apply the efficient approach 'Schur after MOR' to generate a stable ROM of the piezoelectric energy harvester model. "Schur after MOR" projects the reduced system (2) onto a sorted orthonormal eigenbasis of matrix *Mr*:

$$\widetilde{\sum}\_{r} \begin{cases} \underbrace{T^T \mathcal{M}\_r T \cdot \ddot{q}}\_r + \underbrace{T^T E\_r T \cdot \dot{q}}\_r + \underbrace{T^T \mathcal{K}\_r T \cdot q}\_{\widetilde{E}\_r} = \underbrace{T^T \mathcal{B}\_r \cdot u}\_{\widetilde{B}\_r} \\ y = \underbrace{\mathcal{C}\_r T}\_{\widetilde{C}\_r} \cdot q \\ \end{cases} \tag{3}$$

where *T* ∈ R*r*×*r* and *M r* are the eigenvector matrix and sorted diagonal eigenvalue matrix of *Mr*, obtained through eigen decomposition of *Mr* = *TM <sup>r</sup>TT*. According to the indexes of the relative small eigenvalues in *M r*, all the system matrices in (3) can be partitioned:

$$
\begin{aligned}
\widetilde{M}\_{\boldsymbol{r}} &= \begin{bmatrix} \widetilde{M}\_{1} & 0 \\ 0 & \widetilde{M}\_{2} \end{bmatrix} \approx \begin{bmatrix} \widetilde{M}\_{1} & 0 \\ 0 & 0 \end{bmatrix}; \widetilde{E}\_{\boldsymbol{r}} = \begin{bmatrix} \overline{E}\_{11} & \overline{E}\_{12} \\ \overline{E}\_{21} & \overline{E}\_{22} \end{bmatrix} \approx \begin{bmatrix} \widetilde{E}\_{11} & 0 \\ 0 & 0 \end{bmatrix}; \\\ \widetilde{K}\_{\boldsymbol{r}} &= \begin{bmatrix} \widetilde{K}\_{11} & \overline{K}\_{12} \\ \overline{K}\_{21} & \overline{K}\_{22} \end{bmatrix}; \widetilde{B}\_{\boldsymbol{r}} = \begin{bmatrix} \overline{B}\_{1} \\ \overline{B}\_{2} \end{bmatrix}; \widetilde{C}\_{\boldsymbol{r}} = \begin{bmatrix} \widetilde{C}\_{1} & \widetilde{C}\_{2} \end{bmatrix}
\end{aligned} \tag{4}
$$

where *M* 1, *E* 11,*K* 11 ∈ <sup>R</sup>*I*×*I*, *M* 2, *E* 22,*K* 22 ∈ <sup>R</sup>(*<sup>r</sup>*−*<sup>I</sup>*)×(*<sup>r</sup>*−*<sup>I</sup>*), *<sup>E</sup>*12,*<sup>K</sup>*12 ∈ <sup>R</sup>*I*×(*<sup>r</sup>*−*<sup>I</sup>*), *<sup>E</sup>*21,*<sup>K</sup>*21 ∈ <sup>R</sup>(*<sup>r</sup>*−*<sup>I</sup>*)×*I*, *B*1 ∈ <sup>R</sup>*I*×*p*, *B*2 ∈ <sup>R</sup>(*<sup>r</sup>*−*<sup>I</sup>*)<sup>×</sup>*p*, *C*1 ∈ <sup>R</sup>*q*×*I*, *C*2 ∈ <sup>R</sup>*q*<sup>×</sup>(*<sup>r</sup>*−*<sup>I</sup>*), and *I* ∈ [1,*r*]. In order to reconstruct the reduced system in an analogous manner as system (1), the relatively small part *M* 2 and parts *E* 12, *E* 21, *E* 22 in (4) are all neglected. In this way, the Schur complement transformation can be performed and the stable piezoelectric energy harvester ROM is obtained:

$$\sum\_{r\_-\text{Schur}} \begin{cases} \widetilde{M}\_1 \cdot \ddot{q}\_1 + \widetilde{E}\_{11} \cdot \dot{q}\_1 + \left( \widetilde{K}\_{11} - \widetilde{K}\_{12} \widetilde{K}\_{22}^{-1} \widetilde{K}\_{21} \right) \cdot q\_1 = \left( \widetilde{B}\_1 - \widetilde{K}\_{12} \widetilde{K}\_{22}^{-1} \widetilde{B}\_2 \right) \cdot u \\ y = \left( \widetilde{C}\_1 - \widetilde{C}\_2 \widetilde{K}\_{22}^{-1} \widetilde{K}\_{21} \right) \cdot q\_1 + \left( \widetilde{C}\_2 \widetilde{K}\_{22}^{-1} \widetilde{B}\_2 \right) \cdot u \end{cases} \tag{5}$$
