# Nonlinear Dynamics of a Piezoelectric Bistable Energy Harvester Using the Finite Element Method
## 2 Piezoelectric Energy Harvester
### Finite Element Model  
The formulation of a nonlinear finite element method can be expressed by the following generic equation by considering the displacement \(u\), and the mass matrix, \(M\), the damping matrix, \(C\), and a forcing term, \(f\),  
\[M\ddot{u}_{n+1}+C\dot{u}_{n+1}+\dot{f}_{n+1}(u_{n+1})=f_{n+1} \tag{12}\]  
where \(\dot{f}\) represents a nonlinear stiffness that depends on the displacement. The Newmark method can be employed to solve this equation, with the residual vector \(\boldsymbol{\tau}_{n+1}(u_{n+1})\) being expressed by  
\[\boldsymbol{\tau}_{n+1}(u_{n+1})=f_{n+1}-\dot{f}_{n+1}-M\ddot{u}_{n+1}-C\dot {u}_{n+1} \tag{13}\]  
\begin{table}
\begin{tabular}{c c} \hline \hline
**Parameters** & **Piezoceramic** \\ \hline Beam-magnet horizontal spacing (mm) & \(d_{m}=4.7\) \\ Spacing between magnets (mm) & \(d_{g}=14\) \\ Magnets length (mm) & \(l_{m}=4\) \\ Magnets width (mm) & \(b_{m}=4\) \\ Magnetization (A/m) & \(12\times 10^{5}\) \\ Vacuum permeability (N/A\({}^{2}\)) & \(1.26\times 10^{-6}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Properties of magnets.  
The Newton-Raphson method is employed in the solution,  
\[\mathbf{r}_{n+1}\Big{(}\mathbf{u}_{n+1}^{k}\Big{)}+\frac{\partial\mathbf{\tau}\Big{(}\mathbf{u}_{n +1}^{k}\Big{)}}{\partial\mathbf{u}_{n+1}}\Delta\mathbf{u}_{n+1}^{k}=0 \tag{14}\]  
where \(\Delta\mathbf{u}_{n+1}^{k}\) is an increment of \(\mathbf{u}_{n+1}\) as described in detail in Ref. [47].  
The energy harvesting system is modeled and analyzed using ANSYS Academic Research 17. Figure 3 shows the FE model of the harvester presented in Figure 1. The substrate beam is modeled with 2D structural element PLANE82, which has 8 nodes and 2 DOFs per node, presenting translation motion in \(x\) and \(y\) directions. The PZT material is modeled using 2D 8 nodes element, PLANE223, having two translational DOFs and one additional DOF of voltage for each node. In addition, a resistor element (CIRC94) is connected to the PZT elements. A dissipation mechanism is introduced assuming damping proportional to system's mass and stiffness (Rayleigh method) defined by parameters \(\alpha\) and \(\beta\), respectively.  
Base excitation is assumed to be harmonic, being imposed by coupling the nodes located at \(x=0\) in a component and applying displacement boundary condition. In addition, the voltage DOF on the upper and lower surfaces are coupled to provide uniform electrical potentials and thus emulate the electrodes of each piezoelectric layer.  
The magnetic force exerted at the beam free-end is the essential nonlinear effect involved in the energy harvesting device. This magnetic force is modeled using a nonlinear spring element as proposed by [PERSON] and [PERSON] [43]. ANSYS COMBIN39 element is employed to simulate the magnetic interaction by approximating the nonlinear force with piecewise linear segments. Figure 4 shows the variation in magnetic force (\(f_{Mz}\)) with lateral displacement using Equation (1) and the linear fit supplied to COMBIN39 element. This force is applied at a node in the free-end-neutral-axis of the beam. Parameter values used to calculate \(f_{Mz}\) are presented in Table 2.  
Figure 4: Magnetic force variation with lateral displacement.  
Figure 3: Schematic of FE model of energy harvester with a magnetic force \(f_{M}\), a resistive load \(R\), boundary displacement conditions \(U\) and coupled voltage DOF to emulate the electrodes.