BUCKLING AND NON-LINEAR VIBRATIONS OF A MEMS
BIOSENSOR

Olivier Thomas

CNAM – Structural Mechanics
and Coupled System Lab.
Paris, France
olivier.thomas@cnam.fr

Liviu Nicu, Cédric Ayela

CNRS – Laboratoire d'Analyse et
Architecture des Systèmes,
Toulouse, France
ayela@laas.fr, nicu@laas.fr

Cyril Touzé

ENSTA Unité de Mécanique,
Palaiseau, France
cyril.touze@ensta.fr

Abstract

The large amplitude non-linear vibratory behavior of a stratified circular plate with a piezoelectric patch is addressed in this study. The goal is to build an efficient model of a resonating MEMS (Micro Electro Mechanical System) bio-sensor whose function is to detect in real-time the presence of a given molecule (the so called target molecule) in a liquid medium. This study introduces the plate equations of motion, that includes the effects thermal-like prestresses, the piezoelectric material, its stratified structure with two angular regions and the geometrical non-linearities stemming from large deflections. Their resolution by an expansion onto the eigenmodes of the unstressed system gives the first insights in the buckling phenomena and the nonlinear vibrations.

Key words

Non-linear vibrations, buckling, stratified plate,
piezoelectric materials

# 1 Introduction

<figure><img src="image_1.png" /><figcaption>Figure 1. A 2 × 2 resonant membranes matrix biosensor</figcaption></figure>

The large amplitude non-linear vibratory behavior of
a stratified circular plate with a piezoelectric patch is
addressed in this study. The goal is to build an effi-
cient model of a resonating MEMS (Micro Electro Me-
chanical System) bio-sensor whose purpose is to detect

in real-time the presence of a given molecule (the so
called target molecule) in a liquid medium.

The bio-sensor under study is constituted of a matrix
of micromachined stratified circular membranes (Ayela
and Nicu, 2006) (see Fig. 1). Each membrane can
be individually actuated on its fundamental vibration
mode by a piezoelectric thin patch. A biorecognition
molecule, capable of immobilizing the target molecule,
is glued on the top face of the membrane. When the
membrane is in contact with a liquid medium that con-
tains the target molecule, the latter is trapped, thus
adding mass and lowering the fundamental resonance
frequency. Detecting and measuring this frequency
shift lead to estimate the concentration of the target
molecule in the aqueous solution. At present, whereas
the manufacturing of the biosensor is well overcomed,
its fine dynamical behavior still needs to be modeled
and simulated (Ayela and Nicu, 2006) in order to de-
sign such a biosensor.

<figure><img src="image_2.png" /><figcaption>Figure 2. Static non-planar geometry of a membrane, from experi-
ments</figcaption></figure>

Two main issues have to be modeled. First, the
manufacturing processes of the membranes create pre-
stresses, the intensity and sign of which differ from one
layer to another. A consequence is that the shape of
the membrane at rest is not plane and that buckling has
to be taken into account (see Fig. 2). Secondly, non-
linear large amplitude vibrations have been experimen-
<figure><img src="image_3.png" /><figcaption>Figure 3. Dynamic non-linear behaviour of a membrane, from experiments</figcaption></figure>

tally observed in normal use of the sensor, in the form of Duffing-like resonance curves (see Fig. 3). Similar phenomena are also encountered in beam piezoelectric MEMS resonators (Li and Balachandran, 2006).

The membrane is modeled as a stratified circular plate, clamped at its boundary, with one piezoelectric layer and the other composed of linear-elastic-homogeneous-isotropic materials. Von-Kármán-like equations with a kirchhoff-love kinematics across all layers (Reddy, 1997), are used. Those PDEs are expanded onto the normal mode basis of the plate without prestress. A set of second order differential equations, coupled by linear, quadratic and cubic terms, is obtained. The static solutions lead to predict the geometry of the buckled equilibrium position. Then, the modes of the membrane (frequencies and shapes) in post-buckling vibrations can be obtained. Finally, the dynamic non-linear behaviour of the system, in the vicinity of the fundamental mode of the prestressed membrane, is obtained with the framework of non-linear normal modes (Thomas et al., 2005; Touzé and Ama-bili, 2006). In this article, only the general problem formulation and the buckled solution is exposed. Results on the non-linear vibratory behavior of the membranes together with some experimental validations will be proposed at the colloquium.

## 2 Governing equations

<figure><img src="image_4.png" /><figcaption>Figure 4. Cross section of a layered membrane.</figcaption></figure>

A circular stratified plate with a circular piezoelectric (PZT) film at centre is considered (Fig. 4). The

silicon (Si), silicon oxyde (SiO₂) and electrode (Ti/Pt) layers are constituted of isotropic homogeneous materials. The electrodes, plugged to an external voltage generator, create in the PZT layer an electric field in the transverse-z direction. This enables both to polarize the PZT layer (in the z direction) and to deform it in the radial-e_r direction to drive the system in vibration. The PZT layer material is homogeneous and transversely-isotropic in the longitudinal plane of the plate. The classical laminated plate theory (Reddy, 1997) is used, with von-Kármán-like non-linear strain-displacement relationships. The electric field, directed in the transverse-z direction, is assumed uniform through the thickness of the plate. An axisymmetric geometry is considered and a polar (r, θ, z) coordinate system is used. The plate is composed of two parts: a centre circular region, for r ∈ [0, R₁] and an annular outer region for r ∈ [R₁, R₂] (Fig. 4). The plate is clamped at its circular edge, for r = R₂.

We consider here only axisymmetric deformations of the plate. The strains in a point of coordinates (r, θ, z) write $\epsilon = \{\epsilon_r \epsilon_\theta\}^T = \epsilon + z\kappa$, with $\epsilon = \{\epsilon_r \epsilon_\theta\}^T$, $\kappa = \{\kappa_r \kappa_\theta\}^T$ and

$$ \epsilon_r = u_{r,r} + 1/2 w_{r,r}^2, \quad \epsilon_\theta = u_r/r, \quad (1) $$

$$ \kappa_r = -w_{rr}, \quad \kappa_\theta = -w_r/r. \quad (2) $$

In the above equations, $u_r$ and $w$ are respectively radial and transverse displacements and (·),r is the first derivative of (·) with respect to r. The constitutive relation for the linear piezoelectric material writes:

$$ \sigma = Q(\epsilon + z\kappa) - \sigma_{0P}1, \quad (3) $$

with

$$ \sigma_{0P} = \sigma_0 + \frac{e_{31}}{h_p} V(t), \quad (4) $$

where $\sigma = \{\sigma_r \sigma_\theta\}^T$ are the stresses, $Q$ is the plate stiffness operator and $1 = \{1 \ 1\}^T$. $\sigma_0$ represents the manufacturing prestresses, $V(t)$ is the voltage applied between the two Ti/Pt electrodes, $e_{31}$ is the piezoelectric constant and $h_p$ is the PZT layer thickness. By integrating Eq. (3) over the thickness, the following constitutive relations are obtained, with respect to the membrane forces $N = \{N_r N_\theta\}^T$, the moments $M = \{M_r M_\theta\}^T$ and the stiffness matrices $A$ (extensional), $B$ (bending-extensional) and $D$ (bending) (Reddy, 1997):

$$ \begin{pmatrix} N \\ M \end{pmatrix} = \begin{pmatrix} A & B \\ B & D \end{pmatrix} \begin{pmatrix} \epsilon \\ \kappa \end{pmatrix} - \begin{pmatrix} N_{0P}1 \\ M_{0P}1 \end{pmatrix} \quad (5) $$

where $N_{0P}$ and $M_{0P}$ are tensions and moments created by both the prestresses and the piezoelectric converse effect, stemming from the $\sigma_{0P}$ term in Eq. (3).
By inverting Eq. (5), one obtains (Chia, 1980):

$$ \begin{pmatrix} \epsilon \\ M \end{pmatrix} = \begin{pmatrix} A^* & B^* \\ -B^* & D^* \end{pmatrix} \begin{pmatrix} N \\ \kappa \end{pmatrix} + \begin{pmatrix} A^+ N_{0P} & 1 \\ -\tilde{M}_{0P} & 1 \end{pmatrix} \quad (6) $$

with $A^* = A^{-1}$, $B^* = -A^{-1}B$, $D = D + BB^*$, $A^+ = A_{11}^* + A_{12}^*$, $\tilde{M}_{0P} = B^+ N_{0P} + M_{0P}$ and $B^+ = B_{11}^* + B_{12}^*$. Substituting Eqs. (6) in the equilibrium equations of the plate and using dimensionless variables leads to obtaining the following equations of motions, in term of transverse displacement $w$ and Airy stress function $F$:

$$ L_D(w) + \varepsilon L_B(F) + \Delta M_{0P} + \varepsilon \Delta(B^+ N_{0P}) + \bar{m} \ddot{w} = \varepsilon L(w, F) + p, \quad (7) $$

$$ L_A(F) - L_B(w) + \Delta(A^+ N_{0P}) = -\frac{1}{2} L(w, w). \quad (8) $$

$L_A$, $L_B$ and $L_D$ are differential operators that include the discontinuity of matrices $A$, $B$ and $D$ at $r = R_1$. At any point out of this discontinuity, the operators write $L_A(\circ) = A_{11}^* \Delta \Delta \circ$, $L_B(\circ) = B_{12}^* \Delta \Delta \circ$ and $L_D(\circ) = D_{11}^* \Delta \Delta \circ$ with $\Delta$ the laplacian operator. $L(\circ, \circ)$ is the classical bilinear operator found in von Kármán-like equations of motion (Touzé et al., 2002; Thomas and Bilbao, 2008). $\bar{m}$ is the mass per unit surface of the plate, $\ddot{w}$ is the second derivative of $w$ with respect to time and $p$ is an external pressure. $F$ is related to $N$ by $N_r = F_r/r$ and $N_\theta = F_\theta$.

The clamped boundary conditions impose that $w$ and $F$ be finite in $r = 0$ and that at the edge, in $r = R_2$: $w = 0$, $w_r = 0$, $\epsilon_\theta = 0$ and $\epsilon_r - (r\epsilon_\theta)_r = 0$.

### 3 Modal expansion

The transverse unknown displacement $w$ and the stress function $F$ is expanded onto ad-hoc functions:

$$ w(r,t) = \sum_{s=1}^{+\infty} \Phi_s(r) q_s(t), \quad (9a) $$

$$ F(r,t) = \sum_{s=1}^{+\infty} \Psi_s(r) q_s(t) + \sum_{s=1}^{+\infty} \Upsilon_s(r) \eta_s(t). \quad (9b) $$

The $(\Phi_s, \Psi_s, \omega_s)$ are the eigenmodes of the short-circuited plate with zero prestresses ($V(t) \equiv 0$, $\sigma_0 = 0$), solutions of:

$$ L_D(\Phi_s) + \varepsilon L_B(\Psi_s) - \bar{m} \omega_s^2 \Phi_s = 0, \quad (10) $$

$$ L_A(\Psi_s) - L_B(\Phi_s) = 0. \quad (11) $$

together with appropriate clamped boundary conditions. The $(\Upsilon_s, \zeta_s)$ are solutions of:

$$ L_A(\Psi_s) - \zeta_s^4 \Upsilon_s = 0, \quad (12) $$

together with appropriate clamped boundary conditions. By inserting Eq. (9a,b) into the equations of motions (7) and (8) and using the orthogonality properties of the $(\Phi_s, \Psi_s, \omega_s)$ and $(\Upsilon_s, \zeta_s)$, one obtains the equivalent set of equations verified by unknown modal coordinates $q_s(t)$, $\eta_s(t)$, $\forall s \in \mathbb{N}^*$:

$$ \eta_s = -n_0^s - n_P^s V - \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} G_{pq}^s q_p q_q, \quad (13) $$

$$ \ddot{q}_s + 2\xi_s\omega_s\dot{q}_s + \omega_s^2 q_s = -A_0^s - A_P^s V - \sum_{p=1}^{+\infty} C_p^s \eta_p + \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} (D_{pq}^s q_p q_q + E_{pq}^s q_p \eta_q) + Q_s. \quad (14) $$

A modal damping term has been added in the above equations $(\xi_s$ is the damping factor of the $s$-th. mode). All coefficients appearing in the above equations $(n_0^s, n_P^s, C_{pq}^s, A_0^s, A_P^s, C_p^s, D_{pq}^s$ and $E_{pq}^s)$ are functions of $(\Phi_s, \Psi_s, \omega_s, \Upsilon_s, \zeta_s)$, with known analytical expressions, not reported here for a sake of brevity.

The initial partial differential equations (7,8) have been replaced by the equivalent discretized problem (13,14) of coupled ordinary differential equations with quadratic non-linearities. In the above formulation, the $\{\eta_s\}_{s \in \mathbb{N}^*}$ can be eliminated to obtain the following set of coupled ordinary differential equations with quadratic and cubic non-linearities, for all $s \in \mathbb{N}^*$:

$$ \ddot{q}_s + 2\xi_s\omega_s\dot{q}_s + \omega_s^2 q_s = k_0^s + k_P^s V - \sum_{q=1}^{+\infty} [\alpha_{0q}^s + \alpha_{Pq}^s V] q_q + \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} \beta_{pq}^s q_p q_q - \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} \sum_{r=1}^{+\infty} \Gamma_{pqr}^s q_p q_q q_r + Q_s. \quad (15) $$

In the above equations, the coupling terms stem from different sources. Geometrical non-linearities create the cubic terms $(\Gamma_{pqr}^s)$, in the same manner as for a homogeneous plate (Touzé et al., 2002; Thomas and Bilbao, 2008). The layered structure of the plate, if it is non-symmetrical in the transverse-z direction, adds quadratic non-linearities $(\beta_{pqr}^s)$. The voltage applied to the piezoelectric layer can be written as $V(t) = \hat{V} + \tilde{V}(t)$, with a constant part $\hat{V}$, that creates the polarization, superimposed to an oscillating part $\tilde{V}(t)$. Both $\hat{V}$ and the prestresses $(N_{0P}, M_{0P})$ add a constant forcing $(k_0^s, k_P^s \hat{V})$ as well as linear terms $(\alpha_{0q}^s, \alpha_{Pq}^s \hat{V})$. The constant terms are responsible of the non-plane static position of the plate. The linear terms modifies the static position, can create buckling and shifts the natural frequencies $\omega_s$ of the plate. The dynamic piezoelectric driving stemming from voltage $\tilde{V}(t)$ creates a direct forcing $(k_P^s \hat{V}(t))$ as well as a parametric excitation $(\alpha_{Pq}^s \hat{V}(t))$ of the plate. $Q_s$ stems from the external pressure $p$.
## 4 Discretized equations

A solution to the problem defined in the previous section is obtained by writing:

$$w(r,t) = \hat{w}(r) + \tilde{w}(r,t) = \sum_{s=1}^{+\infty} \Phi_s(r) [\hat{q}_s + \tilde{q}_s(t)], \quad (16a)$$

$$F(r,t) = \hat{F}(r) + \tilde{F}(r,t) = \sum_{s=1}^{+\infty} (\Psi_s(r) [\hat{q}_s + \tilde{q}_s(t)] + \Upsilon_s(r) [\hat{\eta}_s + \tilde{\eta}_s(t)]), \quad (16b)$$

where $\hat{w}(r)$ is the static position of the plate and $\tilde{w}(r,t)$ is related to its oscillations with respect to its position at rest $\hat{w}(r)$. $\hat{F}$ and $\tilde{F}$ are the static and dynamic part of stress function $F$, linked respectively to $\hat{w}$ and $\tilde{w}$.

The static solution is obtained by inserting Eq. (16a,b) into (13,14) with $\tilde{w} = \tilde{F} = 0$. One obtains the following set of algebraic equations with quadratic non-linearities, for all $s \in \mathbb{N}^*$:

$$\hat{\eta}_s = -n_0^s - n_P^s \hat{V} - \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} G_{pq}^s \hat{q}_p \hat{q}_q, \quad (17)$$

$$\omega_s^2 \hat{q}_s = -A_0^s - A_P^s \hat{V} - \sum_{p=1}^{+\infty} C_p^s \hat{\eta}_p + \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} (D_{pq}^s \hat{q}_p \hat{q}_q + E_{pq}^s \hat{q}_p \hat{\eta}_q). \quad (18)$$

Again, the $\hat{\eta}_s$ can be eliminated between Eqs. (17) and (18), to obtain the following set of algebraic equations, with quadratic and cubic non-linearities, for all $s \in \mathbb{N}^*$:

$$\omega_s^2 \hat{q}_s = k_0^s + k_P^s \hat{V} - \sum_{q=1}^{+\infty} [\alpha_{0q}^s + \alpha_{Pq}^s \hat{V}] \hat{q}_q + \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} \beta_{pq}^s \hat{q}_p \hat{q}_q - \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} \sum_{r=1}^{+\infty} \Gamma_{pqr}^s \hat{q}_p \hat{q}_q \hat{q}_r. \quad (19)$$

Solving the two above sets of equations (either Eqs. (17,18) or Eqs. (19)) leads to obtain the $\hat{q}_s$ and consequently the system static deflection $\hat{w}$. If $\hat{w} \neq 0$, the non-planar static position of the membrane depicted on Fig. 2 can be obtained.

Then, inserting Eq. (16a) in (15) with $\tilde{w} \neq 0$ and using Eqs. (19) leads to obtain the dynamic part $\hat{w}$ of the system deflection. For this, the following set of equa-

tions is obtained, for all $s \in \mathbb{N}^*$:

$$\ddot{\tilde{q}}_s + 2\xi_s\omega_s\dot{\tilde{q}}_s + \omega_s^2\tilde{q}_s = k_P^s \tilde{V} - \sum_{q=1}^{+\infty} [\alpha_{Pq}^s \tilde{V} + \tilde{\alpha}_q^s] \tilde{q}_q + \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} \tilde{\beta}_{pq}^s \tilde{q}_p \tilde{q}_q - \sum_{p=1}^{+\infty} \sum_{q=1}^{+\infty} \sum_{r=1}^{+\infty} \Gamma_{pqr}^s \tilde{q}_p \tilde{q}_q \tilde{q}_r + Q_s, \quad (20)$$

where

$$\tilde{\alpha}_p^s = \alpha_{0p}^s + \alpha_{Pp}^s \hat{V} - \sum_{i=1}^{+\infty} 2\beta_{ip}^s \hat{q}_i + \sum_{i=1}^{+\infty} \sum_{j=1}^{+\infty} (2\Gamma_{ijp}^s + \Gamma_{ipj}^s) \hat{q}_i \hat{q}_j, \quad (21)$$

$$\tilde{\beta}_{pq}^s = \beta_{pq}^s - \sum_{i=1}^{+\infty} (\Gamma_{pqi}^s + 2\Gamma_{ipq}^s) \hat{q}_i. \quad (22)$$

We recall that the $\{\hat{q}_i\}_{i \in \mathbb{N}^*}$ in the above equations are the solutions of either the set (17,18) or the set (19).

## 5 Resolution

The method of resolution is here described. In a first step, all the coefficients $k_0^s, k_P^s, \alpha_{0q}^s, \alpha_{Pq}^s, \beta_{pq}^s$ and $\Gamma_{pqr}^s$ have to be numerically computed. To do this, the modes $(\Phi_s, \omega_s)$ of the short-circuited plate with no prestresses, with the geometrical and material discontinuity in $r = R_1$, have been calculated in a semi-analytical manner (with Bessel functions), by connecting two continuous annular plates in $r = R_1$. An example of expansion functions is given in Fig. 11.

In a second step, the set (19) have to be solved, to obtain the system static shape $\hat{w}(r)$ as a function of the prestress. This can be done with the Asymptotic Numerical Method (ANM) (Cochelin et al., 1994), a powerful continuation method that as recently been implemented in the Matlab software environment, under the name MANLab (Arquier, 2005; Arquier et al., 2006). As this method mandatory needs a set of equations with non-linearities up to the second order, the formulation of Eqs. (17,18) has to be used.

In a third step, a diagonalization of the linear part of Eqs. (20) lead to calculate the natural frequencies of the plate as a function of the prestresses. Finally, the fourth step is the resolution of the full non-linear set (20), in the particular case of a resonating piezoelectric excitation in the vicinity of the fundamental frequency of the plate, can be conducted with the formalism of the non-linear normal modes (Touzé and Amabili, 2006), in order to obtain the curves of Fig. 3. This final step is presently under progress and will be presented at the colloquium. In the following section, two simple systems have been selected to test the above described method of resolution.
# 6 Results

## 6.1 A perfect circular plate

<figure><img src="image_5.png" /><figcaption>Figure 5. Sketch of the prestressed circular plate</figcaption></figure>

As a preliminary test of the method of resolution, a homogeneous circular plate prestressed by a uniform compression $N_0$ and a uniform external pressure $p_0$ is considered (Fig. 5). In this case, $V = 0$, $M_{0P} = 0$, $N_{OP} = N_0$ so that $n_0^s = -N_0n_s$ and $A_0^s = 0$. Moreover, because of the system transverse symmetry (the symmetry in direction $e_z$), the bending-extensional matrix $B = 0$ which lead to $C_p^s = D_{pq}^s = 0$. The static solution thus verifies the following set, for all $i \in \{1, \dots, N_F\}$ and $j \in \{1, \dots, N_w\}$:

$$ \hat{\eta}_i = N_0 n_0^i - \sum_{p=1}^{N_w} \sum_{q=1}^{N_w} G_{pq}^i \hat{q}_p \hat{q}_q, \quad (23) $$

$$ \omega_j^2 \hat{q}_j = \sum_{p=1}^{N_w} \sum_{q=1}^{N_F} E_{pq}^j \hat{q}_p \hat{\eta}_q + p_0 \hat{Q}_j. \quad (24) $$

The above set is derived from Eqs. (17,18), truncated to the first $N_w$ transverse functions $\Phi_i$ and the first $N_F$ membrane functions $\Upsilon_j$. The $\{\hat{q}_j\}_{j \le N_w}$ and $\{\hat{\eta}_i\}_{i \le N_F}$ solutions are obtained with the above cited ANM, with $\lambda = N_0$ as the continuation parameter. In the following, $N_w = 15$ and $N_F = 30$ have been used.

The classical buckling paths of Fig. 6 are first obtained, together with the corresponding static deflection shapes. If $p_0 = 0$, the dimensionless critical loads $\lambda = N_0 = 14.68, 49.22$ and $103.5$ are obtained for the first three buckling modes. Those values correspond to those classically known in the literature (Bloom and Coffin, 2001).

By following the paths corresponding to the first buckling mode with $p_0 = 0$, one obtains the evolution of the plate static deflection shape as a function of the prestress $\lambda = N_0$, by combining the expansion functions $\Phi_j$ with Eq. (16a). Three deflections shapes for three values of $\lambda = N_0$ above the first critical load are depicted on Fig. 7. The diagrams on the right of this figure show the magnitude of each $\hat{q}_j$ for $j \le 15$ (transverse coordinates) and each $\hat{\eta}_i$ for $i \le 30$. One can observe that for small values of $\lambda$, the deflection shape is very close to the first expansion function $\Phi_1$. When $\lambda$ increases, the other mode shapes gain amplitude and modify the whole static deflection shape. Those diagrams also show the convergence of the method: to recover a static deflection shape with a maximum amplitude of the order of ten times the plate thickness $h$,

<figure><img src="image_6.png" /><figcaption>Figure 6. Perfect plate of Fig. 5. Buckling paths as a function of load parameter $\lambda$ and corresponding buckled system deflection profiles for $\lambda = 21.5$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7. Perfect plate of Fig. 5. First buckling mode deflection shapes for several values of loading parameter $\lambda$ and corresponding convergence diagrams</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8. Perfect plate of Fig. 5. Evolution of the first and second natural frequencies as functions of the loading parameter, on the various buckling paths</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9. Perfect plate of Fig. 5. First vibration mode shape for three prestress values together with the corresponding convergence diagrams</figcaption></figure>
$N_w$ = 6 transverse modes and $N_F$ = 15 membrane functions seems enough.

Finally, for each value of $\lambda$ and the corresponding $\{\hat{q}_j\}_{j \le N_w}$, the natural frequencies and corresponding mode shapes of the prestress system can be obtained by diagonalizing the linear part of Eqs. (20). Figure 8 shows the evolutions of the first two natural frequencies as a function of $\lambda$. Figure 9 shows the first system mode shape for three values of $\lambda$ above the critical load as well as the corresponding convergence diagrams.

## 6.2 A circular plate with two piezoelectric patches

<figure><img src="image_10.png" /><figcaption>Figure 10. Sketch of the circular plate with two piezoelectric patches</figcaption></figure>

As a second test, we consider a homogeneous circular plate (made in Silicium) with two identical circular piezoelectric patches (each patch is composed of a PZT layer with two Ti/Pt electrodes) glued at its center on its both sides (Fig. 10). The central region is composed of seven layers (1 central silicium layer, 4 electrodes and 2 PZT layers). The annular outer zone is the continuation of the central silicium layer. The system is considered without prestress and the same static potential difference $\hat{V}$ is applied on the two patches, so that $M_{0P} = 0$ everywhere, $N_{0P} = 0$ in the outer region and $N_{0P}$ is proportional to $\hat{V}$ in the central region. Since the system is transversely symmetric, the bending-extensional matrix $B$ vanishes again and $C_p^s = D_{pq}^s = 0$. The static solution thus verifies the following set, for all $i \in \{1, \dots, N_F\}$ and $j \in \{1, \dots, N_w\}$:

$$ \hat{\eta}_i = \hat{V} n_0^i - \sum_{p=1}^{N_w} \sum_{q=1}^{N_w} G_{pq}^i \hat{q}_p \hat{q}_q, \quad (25) $$

$$ \omega_j^2 \hat{q}_j = \sum_{p=1}^{N_w} \sum_{q=1}^{N_F} E_{pq}^j \hat{q}_p \hat{\eta}_q + p_0 \hat{Q}_j. \quad (26) $$

Figure 11 shows the first three mode shapes, for each of the three mode families ($\Phi_s$, $\Psi_s$ and $\Upsilon_s$) used in the mode expansion of Eq. 9. Figure 12 shows the buckling paths and the static deflection shapes for three values of the continuation parameter $\lambda = \hat{V}$. One can observe that the effect of the localization of the prestresses on the central region also localizes the deflection.

## 6.3 Conclusion

This paper relates the first results of a work in progress. The main characteristics of this model are (i)

<figure><img src="image_11.png" /><figcaption>Figure 11. First three mode shapes of the system of Fig. 10, for each of the three mode families used in the expansion of Eq. 9, with $R_1/R_2 = 0.3$ ('-': plate with piezoelectric patches; '-'-': prestressed homogeneous plate without patches).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12. Circular plate with piezoelectric patches of Fig. 10. Buckling paths and corresponding static deflection profiles ('-': plate with piezoelectric patches; '-'-': prestressed homogeneous plate without patches)</figcaption></figure>

a continuous model that includes the geometrical nonlinearities, the stratified structure, the geometrical and material discontinuity at the connection radius of the two plate annular parts, the manufacturing prestresses and the piezoelectric actuation; (ii) the use of expansion functions that include the discontinuity, in order to increase the convergence; (iii) the simulation of both
buckling and non-linear vibrations.

In the present article, the main milestones of the model have been included. Then, only the first results on buckling have been presented: first of all, a validation study on a perfect circular plate; then, the buckling behavior of a symmetric stratified MEMS with two piezoelectric patches and seven layers of different materials.

Results on non-linear periodic forced vibrations, in the vicinity of the first natural frequencies of the system, are under progress and will be presented at the colloquium. Along with the awaited physical insights in the non-linear behavior of the MEMS membrane, more computational issues will also be addressed. Since the Asymptotic Numerical Method will be used, coupled to a harmonic balance method, in order to transform the non-linear partial differential equations into a non-linear algebraic set of equations, two issues will be investigated. (i) Some convergence studies will be performed: how many modes in the truncated expansions of Eq. 9 and how many harmonics have to be retained to obtain converged results ? (ii) Moreover, the computational time associated to the use of the Asymptotic Numerical Method coupled with the harmonic balance method will be tested, in order to assess the efficiency of our resolution strategy.

## References

Arquier, R. (2005). Manlab : logiciel de continuation interactif (manuel utilisateur). Technical report. Laboratoire de Mécanique et d'Acoustique, CNRS. http://www.lma.cnrs-mrs.fr/~manlab. In french.

Arquier, R., S. Bellizzi, R. Bouc and B. Cochelin (2006). Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes. *Computers and Structures* **84**, 1565–1576.

Ayela, C. and L. Nicu (2006). Micromachined piezoelectric membranes with high nominal quality factors in newtonian liquid media: a lamb's model validation at the microscale. *Sens. Act. B Chem.* **123**(2), 860–868.

Bloom, F. and D. Coffin (2001). *Handbook of thin plate buckling and postbuckling*. Chapman & Hall / CRC.

Chia, C. Y. (1980). *Nonlinear analysis of plates*. Mc Graw Hill.

Cochelin, B., N. Damil and M. Potier-Ferry (1994). Asymptotic numerical methods and Padé approximants for non-linear elastic structures. *International Journal of Numerical Methods in Engineering* **37**, 1187–1213.

Li, H. and B. Balachandran (2006). Buckling and free oscillations of composite microresonators. *J. Microelectromech. Syst.* **15**(1), 42–51.

Reddy, J. N. (1997). *Mechanics of laminated composite plates*. CRC press.

Thomas, O. and S. Bilbao (2008). Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties. *Journal of Sound and Vibration*. submitted.

Thomas, O., C. Touzé and A. Chaigne (2005). Non-linear vibrations of free-edge thin spherical shells: modal interaction rules and 1:1:2 internal resonance. *Int. J. Sol. Struct.* **42**(11-12), 3339–3373.

Touzé, C. and M. Amabili (2006). Non-linear normal modes

for damped geometrically non-linear systems: application to reduced-order modeling of harmonically forced structures. *J. Sound Vib.* **298**(4-5), 958–981.

Touzé, C., O. Thomas and A. Chaigne (2002). Asymmetric non-linear forced vibrations of free-edge circular plates, part 1: theory. *J. Sound Vib.* **258**(4), 649–676.