## 2 Design and Modeling
### Modeling  
The electromechanical coupling dynamical model can be formulated while using the extended Hamilton's method. The kinetic energy of the MBBP-REH can be given by  
\[T=\frac{1}{2}\int_{0}^{l}m_{pb}\Big{(}\dot{u}(x,t)^{2}+\dot{w}(x,t)^{2}\Big{)}dx+ \frac{1}{2}m_{c}\Bigg{(}\dot{u}\Big{(}\frac{l}{2},t\Big{)}^{2}+\dot{w}\Big{(} \frac{l}{2},t\Big{)}^{2}\Bigg{)}+\frac{1}{2}m_{a}\dot{u}(l,t)^{2} \tag{1}\]  
where \(m_{pb}\) is the mass per unit length of the buckled beam, and it can be given as \(m_{pb}=r_{s}b_{s}h_{s}+r_{p}b_{p}h_{p}\big{[}H(x-l_{1})+H\big{(}x-l_{1}-l_{p} \big{)}\big{]}\), \(\rho_{s}\) and \(\rho_{p}\) are the density of the substrate and the piezoelectric patches, \(b_{s}\) and \(b_{p}\) are the width of the substrate and the piezoelectric patches, \(h_{s}\) and \(h_{p}\) are the thickness of the substrate and the piezoelectric patches, \(l_{1}\) and \(l_{p}\) are the location and the length of the piezoelectric patches, \(H(x)\) is the Heaviside function, \(m_{c}\) is the mass of the center magnet, \(m_{a}\) is the mass of the end moving mass with a magnet, \(w(x,t)\) and \(u(x,t)\) are the transverse and longitudinal displacement of the buckled beam, and the superscript dot represents the derivative with respect to the time, and \(l\) is the length of the buckled beam. The potential energy of the buckled beam is expressed as  
\[U=\frac{1}{2}\int_{0}^{l}EI(w^{\prime\prime})^{2}dx-J_{p}V(t)\int_{0}^{l_{p}}w^ {\prime\prime}dx+\frac{1}{2}C_{p}V(t)^{2} \tag{2}\]  
where (\({}^{\prime}\)) denotes the derivative with respect to location \(x,J_{p}\) is the electromechanical coupling term, which can be given as \(J_{p}=d_{31}E_{p}b_{p}h_{p}\big{(}h_{s}h_{p}+h_{p}^{2}\big{)}\big{/}\big{(}2h _{p}\big{)}\big{[}H(x-l_{1})-H\big{(}x-l_{1}-l_{p}\big{)}\big{]}\), where \(d_{31}\) is the piezoelectric strain constant. \(C_{p}=e_{33}^{S}b_{p}l_{p}/\big{(}2h_{p}\big{)}\) is the capacitance of the piezoelectric patches, where \(e_{33}^{S}\) is the permittivity at constant strain. \(V(t)\) is the generated voltage of the piezoelectric patches. The bending stiffness is given by  
\[EI=\frac{E_{s}b_{s}h_{s}^{3}}{12}\big{[}H(x)-H(x-l_{1})+H\big{(}x-l_{1}-l_{p} \big{)}\big{]}+\frac{E_{s}b_{s}h_{s}^{3}+2E_{p}b_{p}\Big{(}3h_{s}^{2}h_{p}+6h_ {s}h_{p}^{2}+4h_{p}^{3}\Big{)}}{12}\big{[}H(x-l_{1})-H\big{(}x-l_{1}-l_{p}\big{)} \big{]} \tag{3}\]  
where \(E_{s}\) and \(E_{p}\) are the Young's modulus of the elastic beam and the piezoelectric layers, respectively. The magnets that are used in the harvester can be modeled as point dipoles by following the reference [37]. The distance \(\mathbf{r}_{i}\) from the \(i\)th rotational magnet to the center magnet can be given as  
\[\mathbf{r}_{i}=(R\sin\beta_{i}-w)\mathbf{\hat{e}}_{x}-(R-R\cos\beta_{i}+d) \mathbf{\hat{e}}_{y} \tag{4}\]  
where \(R\) is the rotational radius of excitation magnets, \(\beta_{i}\) is the angular displacement of the excitation magnet \(E_{i}\), and \(d\) is the minimum center distance between the excitation magnets and the center magnet along the \(y\) direction. For the center magnet attached to the buckled beam, the magnetic moment vectors can be expressed as  
\[\mathbf{\mu}_{C}=M_{C}V_{c}\mathbf{\hat{e}}_{x} \tag{5}\]  
where \(M_{C}\) is the magnitude of the magnetization vector, and can be calculated as \(M_{C}=B_{r}/u_{0}\) where \(B_{r}\) is the residual flux density of the magnet and \(u_{0}\) is the permeability of free space. \(V_{C}\) is the volume of the magnet \(C\).