Analysis of the Nonlinear Response of Piezo-Micromirrors with the Harmonic Balance Method

Abstract

In this work, we address the simulation and testing of MEMS micromirrors with hardening and softening behaviour excited with patches of piezoelectric materials. The forces exerted by the piezoelectric patches are modelled by means of the theory of ferroelectrics developed by Landau–Devonshire and are based on the experimentally measured polarisation hysteresis loops. The large rotations experienced by the mirrors also induce geometrical nonlinearities in the formulation up to cubic order. The solution of the proposed model is performed by discretising the device geometry using the Finite Element Method, and the resulting large system of coupled differential equations is solved by means of the Harmonic Balance Method. Numerical results were validated with experimental data collected on the devices.

1. Introduction

Applications of Micro-Electro-Mechanical Systems (MEMS) micromirrors are virtually countless and rapidly progressing. Among others, one can cite the Intel^® Realsense LiDAR Camera™ (Santa Clara, CA, USA), Microsoft^® Hololens2™ (Redmond, WA, USA) augmented reality glasses, consumer pico projectors [1], automotive Head-Up-Displays, as well as medical applications [2].

Their technology is rapidly progressing, associating to the early stage devices that relied on electrostatic actuation [3,4,5,6,7,8] more recent systems based on magnetic- [9] and piezoactuation [10]. In particular, piezoactuation will play a major role in the next generation of micromirrors and is the object of intensive research. Indeed, frequency stability is the key to successful operation and every source of nonlinear behaviour, if not well dominated, can be disruptive for the device. The precise control of these devices requires accurate knowledge of the frequency drift of the resonant mode when the actuation voltage increases, hence a full characterisation of the mirrors is required, together with the development of predictive tools.

Among typical piezomaterials for applications in MEMS, Aluminium Nitride (AlN) has better compatibility with semiconductor technology while Lead-Zirconate Titanate (PZT) presents a stronger electromechanical coupling and is thus preferred in actuators. Recently, Aluminium Scandium Nitride (AlScN) has been shown to preserve the positive aspects of AlN with improved coupling coefficients [11]. One should also mention the increasing importance of lead-free materials [12], which are nowadays attracting an important research effort.

In this investigation we focus on mirrors rotating about one resonant fast axis. At least three sources of potential nonlinearities can be identified: material, fluidic and geometrical. If ferroelectric materials such as PZT [13] are utilised, material nonlinearities play an important role since electric fields of the order of 10${}^7$ V/m are generated within the piezo film and the hysteretic nonlinear material behaviour is activated. This issue was tackled recently by Frangi et al. [14] who proposed and validated a modelling strategy to account for the effects of polarisation according to the Landau-Devonshire theory of ferroelectrics. The same approach is adopted in this paper. Micromirrors undergo rotations in open air in the range of tens of degrees and activate complex flow patterns, hence the equivalent quality factor is amplitude dependent. Here, the quality factor is identified from experimental data according to simplified formulas (e.g., [15]). In this work, we focus on the third type on nonlinear behaviour. Mechanical structures in large transformations are prone to geometrical nonlinearities and often display hardening or softening behaviour, mode interaction and other complex phenomena such as internal resonances. When a micromirror is actuated at low energy, its oscillation will be proportional to the linear torsional mode of the structure. At increasing energy, however, this proportionality will be lost even if the mirror, in absence of internal resonances, will keep oscillating periodically. Periodic motion of nonlinear dynamic systems falls in the framework of the Nonlinear Normal Modes (NNM) [16,17], which are an extension to the nonlinear regime of the linear modes. Their numerical computation starting from large finite element models can be performed using a combination of shooting and continuation techniques [18,19,20,21] or with the Harmonic Balance Method (HBM) [22,23], namely an extension of Fourier analysis to nonlinear problems.

In the present work, we focus on the simulation of two piezo-micromirrors using a large displacements formulation and by modelling the nonlinear piezoelectric behaviour of the piezoelectric patches using the theory of ferroelectrics developed by Landau and Devonshire. In particular, we enhance the formulation introduced in [14] by accounting for the large displacements the piezoelectric patches are subjected to during device actuation, hence providing a more comprehensive treatment of the consequence of geometric nonlinearities on the dynamic response of piezo MEMS. We prove that the modified formulation allows removing apparent linear frequency drifts observed in full order models when the average polarisation of the piezoelectric patches increases.

The first device, depicted in Figure 1, rotates around a fixed axis given by the torsional springs directly attached to the substrate and when excited at large rotation amplitudes it shows hardening behaviour. On the contrary, in the latter device, represented in Figure 2, the central mirror is suspended to a gimbal-like structure and the position of the rotation axis is a priori unknown. A further difference between this device and the former one is that, when excited at large rotation amplitudes, the dynamic response of the latter is softening. A rather straightforward simulation strategy was applied by Frangi et al. [14] to the first device assuming that the rotation axis is fixed. Indeed, the rotational inertia with respect to the rotation axis is constant under the assumption of a pseudo-rigid body behaviour, and the nonlinear stiffness can be computed by imposing a given angle to the rigid body and computing the nonlinear couple via the finite element method. However, this strategy fails if the rotation axis of the device is unknown, as for the latter mirror analysed in this work. Therefore, we here resort to a full-order model. We apply a finite element discretisation of the linear momentum balance equation defined over the system geometry and then we solve the coupled system of nonlinear differential equations by means of the HBM to compute the steady oscillatory response of the dynamic system. Therefore, we do not apply any simplifying assumption on the motion of the devices.

The remainder of the paper is organised as follows. After a brief description of the micromirrors and of the experimental setup exploited for their characterisation in Section 2, the model describing the dynamic response of piezoelectric micromirrors and its implementation in the framework of the HBM is detailed in Section 3. In Section 4, we address the validation of the proposed model by comparing the numerical and experimental results for the two selected micromirrors. Finally, in Section 5, some conclusions and future prospectives are discussed.

2. Devices and Experimental Data

The tested micromirrors are reported in Figure 1 and Figure 2. They were fabricated by STMicroelectronics^® using monocrystalline silicon with [110] orientation [24] aligned with the ${\mathbf{e}}_1$ axis. The sol-gel PZT patches, the passivation layers and the electrodes were deposited with the P$\epsilon$tra Thin-Film Piezoelectric technology developed by STMicroelectronics^®.

Figure 1 presents the first micromirror, hereafter denoted M1. The monocrystalline silicon layer of the device has a thickness of 20 $\mathsf{\mu }$m. The circular reflective surface of the device has a diameter of 3000 $\mathsf{\mu }$m, and a circular reinforcement is added beneath it to minimise dynamic deformations. The reflective surface is suspended by means of two torsional springs connected to the substrate. The main axes of the springs pass through the center of the mirror. The anchor points of the device are visible in green in Figure 1b. Actuation is achieved by means of four trapezoidal beams with PZT patches highlighted in orange in the image. Patches are organised in two groups labelled 1 and 2, respectively, as highlighted in Figure 1b. Upon actuation, the beams apply through the folded springs a vertical force that excites the torsional mode of the micromirror, as represented in Figure 1c,d. The associated resonance frequency is 1950 Hz.

A picture of the second micromirror, hereafter labelled M2, is reported in Figure 2a and its geometry is schematically reported in Figure 2b. Material and orientation of the device are identical to mirror M1, but the thickness of the silicon layer is 150 $\mathsf{\mu }$m. The reflective surface is a disk with a diameter of 1600 $\mathsf{\mu }$m with a large reinforcement added on the bottom surface to minimise dynamic deformations. The disk is suspended to a flexible gimbal and the device is actuated at its torsional mode, that is a torsion around the ${\mathbf{e}}_1$ axis, by means of eight PZT patches organised into two sets labelled 1 and 2. The resonance frequency of the torsional mode of the structure is 25,500 Hz.

Both mirrors are actuated by means of a unipolar potential law, applied with a shift of the signal phase by $\pi$ to the two patches groups 1 and 2. The actuation laws for the two patches are given by the formulas:

$$V_1=\frac{V_0}{2}(1+cos\left(2\pi ft\right))V_2=\frac{V_0}{2}(1-cos\left(2\pi ft\right))$$

(1)

where $V_0$ denotes the maximum voltage and f is the excitation frequency. The Frequency Response Function (FRF) of the devices was experimentally determined from the opening angle of the micromirros by measuring the deflection of a laser beam incident on the device. The experiments being performed in frequency control allowed measuring only stable branches. For mirror M1, only a forward sweep of the frequency is performed. On the other hand, M2 was measured by performing both an upward and a downward sweep of the frequency. The FRF of mirror M1 is reported in Figure 3a, where the maximum tilting angle for a given frequency is reported. The frequency axis is normalised by the eigenfrequency of the torsional mode, that is 1950 Hz. As highlighted in the chart, the micromirror shows hardening behaviour with formation of an unstable branch for $V_0$ values above 15 V, which is highlighted in the experimental data by the sudden amplitude jump. The dynamic response of mirror M2 is shown in Figure 3b. Contrary to mirror M1, this device shows softening behaviour and it shows an unstable branch at any excitation voltage value. This different behaviour is mainly due to the fact that the axis of rotation is not clearly defined in M2, and the large eccentric reinforcement underneath the mirror induces a coupling with translational motions, thus reducing the apparent stiffness. Regarding the data reported in Figure 3b, the excitation frequency is normalised by the resonance frequency of the torsional mode of the device, that is 25,500 Hz.

As reported in Section 3, modelling of piezoelectric materials with the technique proposed in [14] requires the knowledge of the time history of the average polarization field within the piezoelectric material used to actuated the device. Therefore, the average polarization P was measured on one PZT patch for each device using a standard Sawyer–Tower circuit [14,25] for each $V_0$ value in increasing order. It is worth stressing that the shape of the hysteresis loops is rather thin and somehow similar to that of relaxor ferroelectrics. This behaviour can be justified considering that the sol-gel deposited PZT considered herein presents rather small ferroelectric domains which have a progressive switching, as recently discussed in [26].

Results regarding the fifth cycle collected during both unipolar and bipolar cycles for both devices are reported in Figure 4. The data show the classical hysteretic behaviour of ferroelectric materials, highlighting the nonlinear response of the material with respect to the applied voltage.

To avoid fatigue effects due to switching between negative and positive values of the polarisation, the PZT patches are normally operated with unipolar cycles. This implies that only unipolar polarisation cycles will be used to model these devices. Elaborating the hysteresis cycles in Figure 4, one can obtain the polarisation histories $P_1(V_0,t)$ and $P_2(V_0,t)$ imposed on the two patches, respectively, assuming that the polarisation is not affected by the frequency of the voltage law, which is reasonable for the working frequencies of these devices.

3. Model for Material and Geometrical Nonlinearities

According to Landau–Devonshire theory [13,27], the polarisation induces inelastic strains

$${\mathbf{S}}^P\left[\mathit{P}\right]=\mathbb{Q}:\left(\mathit{P}\otimes \mathit{P}\right)$$

(2)

where $\mathbb{Q}$ denotes the electrostrictive tensor [13]. When a voltage bias is applied to the electrodes, the electric field in the very thin PZT layer can be expressed as $\mathbf{E}=E{\mathbf{e}}_3=-(V/h){\mathbf{e}}_3$, with h thickness of the PZT layer. It is also reasonable to postulate that, on average, the polarisation of the film $\mathit{P}$ also has the single component $P{\mathbf{e}}_3$, leading to a simplified expression of the polarisation strain components in Equation (2):

$$\begin{array}{c}{S}_1^P=S_2^P=Q_{12}{P}^2,S_3^P=Q_{11}{P}^2,S_4^P=S_5^P=S_6^P=0\end{array}$$

(3)

where a Voigt notation is adopted.

According to Saint Venant-Kirchhoff constitutive law, the second Piola–Kirchhoff stress tensor $\mathit{T}$ in the piezo film can be expressed as:

$$\begin{array}{c}\mathit{T}=\mathbb{C}:(\mathit{S}-{\mathit{S}}^P)={\mathit{T}}^S-{\mathit{T}}^P\mathrm{with}{\mathit{T}}^S=\mathbb{C}:\mathit{S},{\mathit{T}}^P=\mathbb{C}:{\mathit{S}}^P\end{array}$$

(4)

where $\mathbb{C}$ denotes the fourth-order elasticity tensor [24]. The difference between the total Green–Lagrange strain tensor $\mathit{S}$ and ${\mathit{S}}^P$ denotes elastic strains. Total strains are defined in terms of the displacement gradient $\nabla \mathit{u}$ as:

$$\mathit{S}=\frac{1}{2}\left(\nabla \mathit{u}+{\nabla }^T\mathit{u}+{\nabla }^T\mathit{u}·\nabla \mathit{u}\right)$$

According to Equation (3), under the assumption of transverse isotropic behaviour, the non-vanishing components of ${\mathit{T}}^P$ are

$$\begin{array}{cc}\hfill & T_1^P=T_2^P=(C_{11}{Q}_{12}+C_{12}{Q}_{12}+C_{13}{Q}_{11})P^2={\alpha }_1{P}^2\hfill \\ \hfill & T_3^P=(C_{33}{Q}_{11}+2C_{13}{Q}_{12})P^2={\alpha }_3{P}^2\hfill \end{array}$$

(5)

In tensor form:

$${\mathit{T}}^P=P^2\mathsf{\Lambda },\mathrm{with}\mathsf{\Lambda }={\alpha }_1({\mathbf{e}}_1\otimes {\mathbf{e}}_1)+{\alpha }_1({\mathbf{e}}_2\otimes {\mathbf{e}}_2)+{\alpha }_3({\mathbf{e}}_3\otimes {\mathbf{e}}_3)$$

(6)

Equation (5) holds pointwise in the PZT thin film, but experiments only provide $P_1(V_0,t)$ and $P_2(V_0,t)$ for each piezo patch, i.e. the time history of an average polarisation. Secondly, even if the value of P should also depend on the local elastic strain $\mathit{S}$, we introduce the assumption that, in presence of extremely large electric fields, the influence of $\mathit{S}$ on P can be neglected. As a consequence $P_1(V_0,t)$ and $P_2(V_0,t)$ acquire the meaning of pointwise polarisation in the film.

We choose to enforce equilibrium conditions for the micromirrors in large rotations applying the principle of virtual power in reference (undeformed) configuration:

$$\begin{array}{c}\hfill {\int }_{{\mathsf{\Omega }}_0}{\rho }_0\ddot{\mathit{u}}·\mathit{w}\mathrm{d}{\mathsf{\Omega }}_0+{\int }_{{\mathsf{\Omega }}_0}\mathit{S}:\mathbb{C}:\mathrm{sym}({\nabla }^T\mathit{w}·\mathit{F})\mathrm{d}{\mathsf{\Omega }}_0-\\ \hfill {\displaystyle {\displaystyle \sum }_{i=1}^{N_P}}{P}_i^2{\int }_{{\mathsf{\Omega }}_P^i}\mathsf{\Lambda }:({\nabla }^T\mathit{w}·\mathit{F})\mathrm{d}{\mathsf{\Omega }}_0=0\forall \mathit{w}\in \mathcal{C}\left(0\right)\end{array}$$

(7)

where ${\mathrm{\Omega }}_0$ is the domain in material configuration, ${\rho }_0$ is the density in material configuration, $\mathit{w}$ is a test function defined over the space of functions $\mathcal{C}\left(0\right)$ that vanishes on the boundary where Dirichlet boundary conditions are prescribed on the unknown displacements $\mathit{u}$, $\mathit{F}=\mathit{I}+\nabla \mathit{u}$ is the transformation gradient, ${\mathsf{\Omega }}_P^i$ denotes the ith piezoelectric material domain and $N_P$ is the number of piezoelectric domains. Equation (7) highlights that the different actuation patches of the system can be subjected to different actuation voltage histories. $\mathit{S}:\mathbb{C}:\mathrm{sym}({\nabla }^T\mathit{w}·\mathit{F})$ generates nonlinearities up to cubic order in terms of the unknown displacement field $\mathit{u}$ and it is the main source of geometrical nonlinearities. The last integral in Equation (7) represents the forcing term due to the piezoactuation, with $P_i$ provided by experiments as commented above. It is worth mentioning that, according to the unipolar polarisation histories in Figure 4, the polarisation oscillates around a non-zero average value which has the same effect as a large pre-stress on the structure. As a consequence, $P_i^2\parallel \nabla \mathit{u}\parallel$ has the same order of magnitude as the terms of internal forces in the left hand side and $P_i^2\mathsf{\Lambda }:({\nabla }^T\mathit{w}·\nabla \mathit{u})$ cannot be neglected. In a reduced order model, such as the one formulated in [14], the average value of the polarisation gets automatically filtrated by the skewsymmetric test function. However, in a full-order simulation such as the one discussed herein, the additional term is required to avoid spurious frequency shifts.

Equation (7) is discretised with a standard finite element approach leading the following system of time-dependent coupled differential equations:

$$\left[\mathbb{M}\right]\left\{\ddot{\mathbb{U}}\right\}+\left[\mathbb{C}\right]\left\{\dot{\mathbb{U}}\right\}+\left\{\mathbb{F}\left(\mathbb{U}\right)\right\}-\sum _{i=1}^{{\mathrm{N}}_p}{\mathrm{P}}_i^2\left[V_i\right]\left({\left\{{\mathbb{F}}^p\right\}}_i+{\left[{\mathbb{K}}^p\right]}_i\left\{\mathbb{U}\right\}\right)=\mathbb{G}(\ddot{\mathbb{U}},\dot{\mathbb{U}},\mathbb{U},P_i^2)=0$$

(8)

where $\left[\mathbb{M}\right]$ is the mass matrix, $\left\{\ddot{\mathbb{U}}\right\}$ is the nodal acceleration vector, $\left\{\mathbb{F}\right(\mathbb{U}\left)\right\}$ is the nonlinear internal force vector, ${\left\{{\mathbb{F}}^p\right\}}_i$ is a vector of piezoelectric forces and ${\left[{\mathbb{K}}^p\right]}_i$ is a piezoelectric stiffness. We add damping in the form $\left[\mathbb{C}\right]\left\{\dot{\mathbb{U}}\right\}$ and select a Rayleigh model:

$$\left[\mathbb{C}\right]=\frac{{\omega }_0}{Q}\left[\mathbb{M}\right]$$

to reproduce the experimentally measured quality factor Q for the resonating mode. ${\omega }_0$ is the torsional mode eigenfrequency.

We stress once more that the piezoelectric force is made by two terms: the first term is constant with respect to the displacement, while the second is linear with respect to the displacement itself. This last term implies that, since the polarisation is a function of time, the stiffness of the system is modulated in time by the piezoelectric force.

Numerical Solution

The system reported in Equation (8) is solved with the HBM [20,21,28,29,30,31,32]. We assume that the solution in time can be represented by a truncated Fourier series of order N:

$$\left\{\mathbb{U}\right\}=\left\{{\mathbb{U}}_0\right\}+\sum _{n=1}^N\left\{{\mathbb{U}}_c^n\right\}cos\left(n\omega t\right)+\left\{{\mathbb{U}}_s^n\right\}sin\left(n\omega t\right)$$

(9)

where N is the order of the harmonic expansion, $\left\{{\mathbb{U}}_0\right\}$, $\left\{{\mathbb{U}}_c^n\right\}$, $\left\{{\mathbb{U}}_s^n\right\}$ are vectors collecting the unknown coefficients of the expansion, and $\omega =2\pi f$. Their values are computed by substituting Equation (9) into Equation (8) and by projecting the residual $\mathbb{G}$ on the same Fourier basis used for the displacement field:

$$\begin{array}{c}{\int }_0^T\mathbb{G}(\ddot{\mathbb{U}},\dot{\mathbb{U}},\mathbb{U},P_i^2)\mathrm{d}t=\left\{0\right\},\\ {\int }_0^{T}cos\left(n\omega t\right)\mathbb{G}(\ddot{\mathbb{U}},\dot{\mathbb{U}},\mathbb{U},P_i^2)\mathrm{d}t=\left\{0\right\},\\ {\int }_0^{T}sin\left(n\omega t\right)\mathbb{G}(\ddot{\mathbb{U}},\dot{\mathbb{U}},\mathbb{U},P_i^2)\mathrm{d}t=\left\{0\right\}\end{array}$$

(10)

for n = 1:N. T is the period of oscillation. The large-scale nonlinear system is solved with a Newton–Raphson method [20]. In particular, $P_i^2$ is itself expressed in the same Fourier basis as in (9) and integrals over the time period of the cubic nonlinear terms and of their linearisation are performed exactly with the trapezoidal rule since the latter method is exponentially convergent for periodic functions.

Regarding the expansion order of the Fourier series, a sufficiently high N is required to analyse micromirrors that experience large rotations, of the order of $\pm 0.3$ rad. This limitation can be intuitively explained by simply considering a rigid mirror rotating about a fixed axis z, with angle $\theta =\mathsf{\Theta }sin\omega t$, where $\mathsf{\Theta }$ is the maximum rotation angle. In Cartesian coordinates, the displacement components of a material point with radial distance r from the axis of rotation can be expressed with a truncated series as:

$$u_x=r(cos\theta -1)\equiv r\sum _{n=1}^N{(-1)}^n\frac{{\theta }^{2n}}{\left(2n\right)!}{u}_y=rsin\theta \equiv r\sum _{n=0}^N{(-1)}^n\frac{{\theta }^{2n+1}}{(2n+1)!}$$

where N should be selected to guarantee the desired accuracy. Since $\theta =\mathsf{\Theta }sin\left(\omega t\right)$ this amounts to considering terms up to $sin\left(N\omega t\right)$, $cos\left(N\omega t\right)$ in the trigonometric series (9). All the simulations reported in the following sections are performed with $N=5$.

4. Numerical Results

The meshes adopted for the numerical simulation are illustrated in Figure 5. They are made by quadratic wedge elements with a total of 35,969 and 29,280 nodes, respectively, corresponding to 107,907 and 87,840 unknowns for each Fourier component. If $N=5$ is selected, this gives a total of 1,186,977 and 966,240 unknowns in the nonlinear system (10) for M1 and M2, respectively.

Numerical simulations were run with $V_0$ values equal to 10, 15 and 20 V for mirror M1 and 20, 25 and 30 V for mirror M2. Mechanical properties for silicon were taken from Hopcroft [24]. On the other hand, since very limited data are available for the material properties of PZT, the Young modulus $E=65$ GPa and Poisson coefficient $\nu =0.34$ were selected to guarantee the best possible agreement with experiments. The piezoelectric coefficients values ${\alpha }_1=-1800\mathsf{\mu }\mathrm{N}\mathsf{\mu }{\mathrm{m}}^2/{\mathrm{pC}}^2$ and ${\alpha }_3=4500\mathsf{\mu }\mathrm{N}\mathsf{\mu }{\mathrm{m}}^2/{\mathrm{pC}}^2$ correspond to the electrostrictive coefficients taken from [13] for a mole fraction of PbTiO${}_3$ equal to $x=0.5$.

An important parameter of the model is the damping, which is quantified here by the quality factor Q associated to the torsional mode of the structure. In the present work Q was calibrated on the basis of experimental data using the simple formula provided by Davis [15] under the assumption of linear damping. The estimated quality factors are summarised in

Table 1. Quality factors obtained from Equation (14) of [15] under the assumption of linear damping.

	M1			M2
V${}_0$ [V]	10	15	20	20	25	30
Q [-]	186	143	105	1386	1200	1092

.

The comparison between numerical and experimental results is reported for mirror M1 in Figure 6. The results obtained for mirror M2 are collected in Figure 7.

Several remarks are worth stressing. Both experiments and numerical tests were run in frequency control, hence unstable branches could not be detected. A continuation procedure, as proposed in [20], could be implemented to follow the frequency response function beyond the folds.

The numerical method predicts the correct amplitude also far from resonance, a regime which is only minimally affected by the quality factors. This represents a strong validation of the piezoactuation model. The simulations also correctly capture the frequency shift of the peaks and the insurgence of instabilities (i.e., folds in the frequency response function), which is a primary feature of the device.

Furthermore, the correct bandwidth of the numerical curves implies that the assumption of linear damping is reasonable at least as a first approximation to model dissipation for the present device, even though a dependence of Q on the maximum opening would be more adequate.

Finally, we focus on the effect of the $P_i^2\Lambda :({\nabla }^T\mathbf{w}·\nabla \mathbf{u})$ term on the predicted response of the structure. As mentioned in Section 3, this term is required to avoid spurious frequency shifts that would otherwise be observed in the numerical curves, an effect which is often present in pre-stressed structures. The results are shown in Figure 8, highlighting the remarkable impact that the large displacements formulation has on the computed FRFs. In particular, we observe that nonlinearity and maximum amplitude are preserved in both models and that the only difference between the two curves is the shift in resonant frequency of the structure.

5. Conclusions

In this work, we address the full-order simulation of piezoelectric micromirrors oscillating at large rotations around a fast axis. Material nonlinearities, due to the hysteretic behaviour of the ferroelectric material, are modelled with the formulation recently proposed in [14]. Geometrical nonlinearities, on the contrary, are modelled with the classical Lagrangian equilibrium condition, which allows for large displacements (and small strains) and generates a formulation with nonlinearities up to cubic order. The resulting model is discretised by means of the finite element method and the associated system of nonlinear equations is solved through the HBM. An experimental campaign was run on two different micromirrors to extract both the frequency response function and the polarisation history in the piezo-patches. The excellent agreement between numerical and experimental data highlights that the proposed approach efficiently models the nonlinear dynamic response of piezo microsystems. In particular, the formulation correctly predicts the softening and hardening behaviour of the mirrors both qualitatively and quantitatively along with the nonlinear shift of the resonating frequency. Only the quality factor of the device was extracted from experimental data according to a simplified assumption of linear damping. Even though this assumption turns out to be sufficiently accurate to predict the overall mirror response, a numerical and experimental campaign is under way using ring-down measurements of the devices and dedicated Arbitrary Lagrangian-Eulerian Navier-Stokes Simulations.