A piezoelectric thin-film microactuator in the form of an asymmetrically laminated diaphragm is developed as an intracochlear hearing aid. Experimentally, natural frequencies of the microactuator bifurcate with respect to an applied bias voltage. To qualitatively explain the findings, we model the lead-zirconate-titanate (PZT) diaphragm as a doubly curved, asymmetrically laminated, piezoelectric shallow shell defined on a rectangular domain with simply supported boundary conditions. The von Karman type nonlinear strain–displacement relationship and the Donnell–Mushtari–Vlasov theory are used to calculate the electric enthalpy and elastic strain energy. Balance of virtual work between two top electrodes is also considered to incorporate an electric-induced displacement field that has discontinuity of in-plane strain components. A set of discretized equations of motion are obtained through a variational approach.

## Introduction

For the last decade, bistability behavior of micro-electrical-mechanical systems (MEMS) has attracted significant attention of researchers [1–17]. One unique phenomenon of these bistable devices is the ability to migrate from one equilibrium state to another via snap-through. The snap-through phenomenon could enable new applications for electrostatic MEMS actuators, such as switches [18–20], microvalves and microrelays [21–23], and band-pass filters [10]. The study of snap-through also facilitates better designs against pull-in instability of electrostatic MEMS actuators. In fact, many MEMS structures are fabricated via high-temperature processes. Thermal stresses resulting from these processes often cause the MEMS structures to warp, thus producing a bistable system potentially for snap-through to occur.

In structural mechanics, snap-through instability is a well-studied topic [24]. Snap-through instability could occur in one-dimensional structures (e.g., arches [25–27]) or two-dimensional structures (e.g., shells [28,29]). For MEMS devices, the study of snap-through instability is focused primarily on electrostatic actuators. For these systems, an applied voltage and a curved, elastic structure form an electrostatic load causing the snap-though. The constitutive equation of the structure follows Hooke's law. There are, however, many other modes of actuation in MEMS devices, such as piezoelectric and thermal actuations.

For a piezoelectric actuator, the applied voltage affects a curved structure via the converse piezoelectric effect. Distinct compared to electrostatic actuators, the applied voltage changes the stiffness and thus the stability of the structure, resulting in rich bistability behaviors. For example, for a bimorph buckled beam under distributed actuation, Maurini et al. [30] showed that depending on the magnitude and phase of axial actuation, the beam can have different numbers of equilibria and the stability of the equilibria change accordingly. Furthermore, the voltage–displacement curve of a bimorph curved beam of partial surface piezoelectric coverage subjected to bias voltages was observed separated into two disjoint portions before and after snap through in a finite element analysis conducted by Varelis and Saravanos; see Fig. 7 of Ref. [31]. None of the above has been observed in electrostatic actuators.

The purpose of this paper is to study snap-through phenomenon of a piezoelectric MEMS actuator. The piezoelectric actuator takes the form of a two-dimensional diaphragm with initial warping, and the snap-through is induced via an applied voltage. This paper is motivated by the development of a lead-zirconate-titanate (PZT), thin-film, acoustic microactuator intended for an intra-cochlear hearing aid [32]; see Fig. 1. The microactuator consists of four parts: a silicon diaphragm, a bulk silicon substrate, a PZT thin-film layer, and a pair of electrodes. (Note that the four parts in Fig. 1 are not drawn to scale.) The diaphragm is anchored to the silicon substrate. The diaphragm can be formed, for example, by etching the silicon substrate from the bottom side to form a cavity under the diaphragm. On top of the silicon diaphragm is a layer of PZT thin film with a pair of bottom and top electrodes. When a driving voltage is applied to the electrodes, the PZT thin film extends and contracts flexing the diaphragm like a speaker.

The diaphragm of the microactuator, however, is not perfectly flat. In fact, it is warped into a thin shell due to thermal stresses induced in the fabrication process; see Fig. 14 of Ref. [33]. As a result, the microactuator is likely to have multiple equilibrium positions to enable snap-through. The snap-through instability, if reached, may cause several concerns. First, when a PZT microactuator is used in cochlea, it may need to be poled regularly to maintain its piezoelectricity. If snap-through instability occurs during poling, it may affect the gain and bandwidth of the microactuator, as different equilibrium states imply different local stiffness and natural frequencies. There is also a risk that a standard 5 V power source could drive the micro-actuator into various domains of attraction under normal operation.

It is not a trivial task to study snap-through instability of such PZT thin-film micro-actuators. There are several major challenges that need to be overcome. First, there must be some preliminary experimental results suggesting that snap-though instability may have occurred in the operation of the PZT thin-film micro-actuator. In other words, theoretical modeling of the snap-through instability should be justified.

If the theoretical modeling is necessary, the second challenge that ensues is lack of existing models. Although general nonlinear piezoelectric shell lamination (multilayered) theories can be found in the literature [31,34–43], when used to derive models, various assumptions are considered. For example, a bimorph configuration is often assumed, where one or multiple elastic layers are symmetrically laminated across the thickness direction, with two identical piezoelectric layers being on the top and bottom surfaces, respectively [37–41,44,45]. (This configuration is referred to as “symmetric lamination” hereafter.) Such a bimorph configuration greatly simplifies modeling of piezoelectric actuation. One immediate example is that it allows for pure bending actuation, where the top and bottom piezoelectric layers extend and contract against each other, excluding in-plane stretching actuation. In fact, for an asymmetrically laminated bi-morph or mono-morph configuration, bending and stretching actuation are coupled with each other, e.g., see Ref. [30], and should be taken into account. Furthermore, for piezoelectric smart structures, the thickness of piezoelectric layers is often much thinner than elastic layers and assumed to be negligible, and so are their elastic properties [44,46,47]. However, as will be shown later in this paper, when the elastic properties of the piezoelectric layers are not negligible, the applied voltage affects the stiffness of a curved structure via the converse piezoelectric effect.

In reality, piezoelectric MEMS devices are often asymmetrically laminated to minimize fabrication complexity. Moreover, the piezoelectric and elastic layers are quite comparable in thickness to maximize actuation strength. For example, the intra-cochlear acoustic microactuator in Ref. [32] has multiple layers, including parylene, gold, chromium, PZT, platinum, titanium, silicon nitride, silicon oxide, and a second parylene layer. The layers above and below the PZT thin film total about 0.5 *μ*m and 1.1 *μ*m, respectively, while the PZT thin film itself is about 1 *μ*m in thickness. By all measures, the PZT actuator diaphragm is asymmetrically laminated and the thickness of the PZT thin film is not negligible.

The third challenge is that a realistic PZT thin-film microactuator employs partial electrodes to maximize actuator displacement; see Ref. [48]. The use of partial electrodes results in discontinuous electric fields when the voltage is applied. Studies of electro-mechanical coupling with continuous electric fields are available in the literature [42,43]. Studies of discontinuous electric fields, however, remain largely open.

When asymmetric lamination and partial electrode coverage or more general configurations are considered, finite element methods are predominantly used [31,35–38,42,43]. Although finite element methods facilitate studying general configurations, physical insights are often lost during the process.

From the discussions above, it is evident that existing theoretical studies of piezoelectric shells are not adequate for the design and development of piezoelectric MEMS actuators. Snap-through instability of PZT thin-film microactuators, induced by an applied voltage, has not been extensively studied despite its potential impact.

For the rest of the paper, we first present preliminary experimental results that suggest occurrence of snap-through instability in the PZT thin-film microactuator. The experimental results demonstrate that vibration modes bifurcate when the applied voltage exceeds a threshold. We then develop a mathematical model to qualitatively explain the bifurcation. Specifically, we model the PZT diaphragm as a doubly curved, asymmetrically laminated, piezoelectric shallow shell defined on a rectangular domain with simply supported boundary conditions. Assumptions and coordinate systems of the shell model are introduced and explained. Kinematic relations and constitutive equations are introduced to formulate electric enthalpy and elastic strain energy of the shell model. Virtual work done by the boundary membrane forces and bending moments are derived and simplified. A particular shape function of in-plane displacement fields that satisfy, in the sense of *weak form*, a set of continuity and equilibrium equations along the boundary of two partial electrodes is derived to incorporate the discontinuous electric fields resulting from the partial electrodes. Finally, variational forms of the potential and kinetic energy are derived, and discretized equations of motion are obtained through Hamilton's principle and the Rayleigh-Ritz method using a set of admissible shape functions and the particular shape function.

## Preliminary Experiment Study

The PZT thin-film microactuator shown in Fig. 1 has been fabricated and tested [32,49]. The diaphragm has a size of 0.8 mm by 0.8 mm and a thickness of 2–3 *μ*m. Moreover, the diaphragm has a bottom electrode and two top electrodes (i.e., one inner and one outer) as illustrated in Fig. 1. Detailed description of the test specimen can be found in Refs. [32] and [49] and is not repeated here. The experiment was done in an aqueous environment. The PZT thin-film microactuator was submerged in water. A voltage drove the microactuator and a laser Doppler vibrometer measured the velocity of the vibrating diaphragm. A spectrum analyzer processed the driving voltage and measured velocity to obtain a frequency response function (FRF) that characterizes the PZT thin-film microactuator. Detailed description of the experimental setup can also be found in Refs. [32] and [49] and is not repeated here again.

With the test specimen and the experimental setup, we conducted a preliminary experiment to investigate possible snap-through instability of the PZT thin-film microactuator. For the first step, a considerable direct current (DC) bias voltage (e.g., 2 V) was applied across the outer and bottom electrodes for a sufficiently long time (>20 min.). In the meantime, the inner electrode remained floating. The hypothesis was that the bias voltage could deflect the actuator diaphragm from one stable equilibrium configuration to another if the bias voltage was large enough. The next step was to remove the DC voltage, and a swept sine measurement was subsequently conducted. For the swept sine measurement, a small voltage (in the order of mV) was applied across the inner and bottom electrodes. At the same time, the outer electrode remained floating. The purpose of the swept sine measurement was to obtain FRFs, from which natural frequencies of the microactuator diaphragm could be extracted experimentally. The hypothesis was that natural frequencies reflected the linearized stiffness coefficients around a stable equilibrium position. If snap-through instability occurred, the actuator diaphragm would move to a different stable equilibrium position. Therefore, the linearized stiffness coefficients would change and so would the natural frequencies. These two steps were repeated by increasing or decreasing the DC bias voltage hoping to maneuver the actuator diaphragm from one stable equilibrium position to another and then back.

Figure 2 shows the measured natural frequencies versus the DC bias voltage varied between 5 V and −5 V. Measured FRFs are also included Fig. 2 for comparison. The experiment started with a zero DC biased voltage (see section 1 in Fig. 2). In the measured FRF, there was only one significant resonance peak with a natural frequency around 24 kHz. As the DC biased voltage increased from 0 to 2.5 V, the resonance peak remained unchanged. When the DC bias voltage exceeded 2.5 V, the FRF changed abruptly (see section 2 in Fig. 2). Specifically, the 24 kHz resonance peak disappeared and two new resonance peaks emerged, one around 34 kHz and the other around 20 kHz. These two resonance peaks did not change significantly when the DC bias voltage was increased to 3.5 V (see section 3 in Fig. 2), indicating the PZT diaphragm had migrated from the initial stable equilibrium position to another stable equilibrium position.

In the second part of the experiment, the DC-biased voltage was reduced and removed. The two resonance peaks remained at around 33 kHz and 19 kHz. In other words, the PZT diaphragm stayed in the second stable equilibrium position and did not return to its initial stable equilibrium state. The measured FRF remained the same even when the biased voltage was reversed (i.e., negative DC bias). When the DC-biased voltage was reduced to −4 V, the FRF experienced another abrupt change. The two resonance peaks at 33 kHz and 19 kHz disappeared. Instead, one new resonance peak appeared at 31 kHz (see section 4 in Fig. 2). The natural frequency and FRF were not affected as the biased voltage was further decreased (see section 5 in Fig. 2) or increased (see section 6 in Fig. 2), indicating the PZT diaphragm migrated from the second stable equilibrium position to a third stable equilibrium position.

In summary, the natural frequencies of the PZT diaphragm bifurcate when the DC bias voltage exceeds a threshold. Moreover, unlike the snap-through observed in electrostatic actuators or elastic structures, the trend of bifurcation is not reversible, suggesting something unique occurred to the piezoelectric actuator. The only tangible explanation is that the DC bias voltage has created bending moments and in-plane stresses that migrate the PZT diaphragm from one stable equilibrium to another via snap-through, and the mechanism that kept it from going back to its original state remains to be investigated. These experimental results indicate that snap-through instability is highly possible for a PZT, thin-film, diaphragm-type microactuator. A detailed mathematical model to analyze snap-through of piezoelectric thin shells is critically needed.

## Assumptions

Figure 3 shows the shallow shell model to be developed in this paper. The shell model consists of three layers. The middle layer 2 is a piezoelectric layer, representing the PZT thin film. Layer 1 is an elastic layer representing all materials below the PZT thin film, such as the bottom electrode, silicon nitride, silicon oxide, and silicon. Layer 3 is an elastic layer representing all materials above the PZT thin film, such as the top electrode and parylene coating. Moreover, layer 3 consists of an inner electrode and an outer electrode neighboring to each other. (For the PZT microactuator in Fig. 1, the gap between the inner and outer electrodes is 20 *μ*m.) For the rest of the paper, quantities related to the *i*th layer are denoted by a superscript (*i*) or a subscript *i*. Furthermore, a special notation $(\xaf)$ is reserved for a reference surface, known as the modulus-weighted midsurface $S\xaf$, which will be defined later in this section. Finally, quantities without any subscript or superscript are related to the shell laminate, which consists of all three layers.

In developing the shallow shell laminate model, the following assumptions and definitions are made:

- (1)
The diaphragm of the PZT micro-actuator is warped due to thermal stresses induced in the fabrication process. Rigorously speaking, the micro-actuator should be modeled as a thermally buckled plate. However, it is a difficult task to estimate the induced thermal stresses. As a first attempt, the micro-actuator is simplified as a shell laminate model which is assumed to be initially stress free, and yet curved with constant curvature

*R*and_{x}*R*(spherical type). The constant curvatures are used to account for the warping._{y} - (2)
For the PZT micro-actuator model in Fig. 3, the actuator domain Ω consists of an inner electrode domain Ω

^{–}and an outer electrode domain Ω^{+}. Furthermore, Ω is considered as a rectangular domain defined as $\Omega ={(x,y)|0\u2264x\u2264a,0\u2264y\u2264b}$ and enclosed by an outer boundary $\u2202S\xaf$, where*a*and*b*represent the domain widths. By the same token, Ω^{–}is defined as $\Omega \u2212={(x,y)|xle\u2264x\u2264xre,yle\u2264y\u2264yre}$, and enclosed by an inner boundary $\u2202\Omega \u2212$, where $xle,xre,yle$ and $yre$ represent the four vertices of the domain. The outer electrode domain is then $\Omega +=\Omega \u2212\Omega \u2212$, and is enclosed by $\u2202S\xaf$ and an inner boundary $\u2202\Omega +$_{,}which coincides with $\u2202\Omega \u2212$. - (3)
For each layer, there exists a local coordinate system $C(i)$ with orthogonal curvilinear coordinates $(\alpha 1,\alpha 2,\alpha 3(i))$; see Fig. 3(b). Furthermore, all three layers are shallow shells satisfying the following assumptions [29].

- (a)
The orthogonal curvilinear surface coordinates $(\alpha 1,\alpha 2,\alpha 3(i))$ can be replaced with the rectangular coordinates $(x,y,z(i))$; see Fig. 3(b).

- (b)
The Lamé parameters of a shallow shell equal unity, i.e., $Ax(i)=Ay(i)=1$, where

*i*= 1, 2, 3. Furthermore, the radii of curvature of each layer are equal, i.e., $R\alpha (1)=R\alpha (2)=R\alpha (3)$, where $\alpha =x,y$. As a result, the differential volume of each layer can be approximated by $dV(i)=dxdydz(i)$.

- (a)
- (4)A global coordinate system $C\xaf$ with coordinates $(x,y,z\xaf)$ is defined, where the
*x*and*y*coordinates are identical to those in the local coordinate systems $C(i)$, while the $z\xaf$ coordinate is obtained by shifting the transverse coordinate $z(i)$; see Fig. 3(b). $C\xaf$ is used to determine the modulus-weighted midsurface $S\xaf$, which serves as a reference to account for the variation of modulus and Poisson's ratio in each layer across the thickness direction. The position of $S\xaf$ can be determined by a standard moment analysis. Let $z\xaf(i)$ be the distance from the midsurface of the*i*th layer to $C\xaf$. Thenwhere $z\xaf(0)$ is the vertical distance from $S\xaf$ to $C\xaf,\u2009Ki=(Y(i)hi/(1\u2212\nu i2))$ is the conventional membrane stiffness, and $Y(i)$,$z\xaf(0)=\u2211i3Kiz\xaf(i)\u2211i3Ki$(1)*h*and_{i}*ν*are the Young's modulus, thickness, and Poisson's ratio of the_{i}*i*th layer, respectively. It is convenient to set the origin of the global coordinate system $C\xaf$ to coincide with a reference point on $S\xaf$ such that the distance $z\xaf(0)$ becomes zero. Consequently,$\u2211i3Kiz\xaf(i)=0$(2)

- (5)
All three layers adopt a geometric nonlinear strain–displacement relationship used in von Karman plates. Moreover, the Donnell–Mushtari–Vlasov theory is used to simplify the model, which is well justified by the shell laminate being shallow [29]. The interested reader is referred to Ref. [46] for the details of the Donnell-Mushtari-Vlasov (DMV) theory.

- (6)
Love's assumptions are made in the kinematic relations for the shell laminate. Consequently, the following kinematic relations can be obtained [34].

- (a)Let $u=(u,v,w)T$ be the displacement vector of the modulus-weighted midsurface, where
*u*,*v*, and*w*are the displacement in the*x*,*y,*and*z*directions, respectively, and*T*indicates transpose. Also, let $\beta =(\beta x,\beta y,\beta z)T$ be the angular rotation vector, where*β*and_{x}*β*are the angles of rotation of the laminate, and_{y}*β*are the normal deformations of the laminate across the thickness. Consider an arbitrary point_{z}*P*in the laminate, located at $(x,y,z\xaf)$ in the global coordinate system. Let $U$ be the displacement vector of*P*. Then$U(x,y,z\xaf)=u(x,y)+z\xaf\beta (x,y)$(4) - (b)Love's assumptions imply that $\beta z=0$. Furthermore, the laminate's strains in the transverse direction vanish, i.e., $Szz=Sxz=Syz=0$ [50]. As a result, the laminate's rotational angles
*β*and curvatures*κ*can be simplified by the DMV theory, i.e.,where $\alpha ,\beta =x,y$, and$\beta \alpha =\u2212\u2202w\u2202\alpha ,\kappa \alpha \beta =\u2212c\u22022w\u2202\alpha \beta $(5)*c*= 2 if $\alpha \u2260\beta $ and*c*= 1 otherwise. - (c)Following the procedure in Ref. [34], the in-plane strains of the laminate can be written in the global coordinate system aswhere $\alpha ,\beta =x,y$,$S\alpha \beta =S\alpha \beta 0+z\xaf\kappa \alpha \beta $(6)
*S*are the laminate's in-plane strains and*S*^{0}are the membrane strains of the modulus-weighted mid-surface. To transform the laminate's in-plane strains to the*i*th layer's local coordinate system, Eq. (3) is substituted into Eq. (6), leading towhere $\alpha ,\beta =x,y$ and $S(i)$ and $S0(i)$ are the$S\alpha \beta (i)=S\alpha \beta 0(i)+z(i)\kappa \alpha \beta $(7)*i*th layer's in-plane strains and membrane strains of the midsurface, respectively. Finally, the normal stresses*σ*_{33}are negligible in all three layers.

- (a)
- (7)
The electric potential $\varphi (x,y,z(2),t)$ inside the piezoelectric layer $S(2)$ is assumed to vary linearly across the thickness direction. In addition, the in-plane electric fields

*E*and_{x}*E*are neglected and only the transverse electric field_{y}*E*is considered._{z}

The shell laminate model is developed using a” local” to” global” procedure. Formulas of individual layers are derived in each local coordinate system $C(i)$ and in terms of quantities associated with each layer $S(i)$, such as the membrane strains $S0(i)$ of the *i*th layer (layer level). The shell laminate is then obtained by assembling individual layers and transforming all the quantities to the modulus-weighted midsurface $S\xaf$ in the global coordinate system $C\xaf$ (laminate level).

## Kinematics and Electric Fields

In this section, the kinematic relations in the layer level are transformed to the global coordinate system $C\xaf$, and then added together to derive those in the laminate level.

### Kinematic Relations.

*i*th layer and the laminate, respectively, and

**T**

*is the coordinate transformation matrix from $C\xaf$ to $C(i)$, defined as*

_{i}where $I$ and $0$ are a 3 × 3 identity and zero matrix, respectively.

*w*. These nonlinear strain–displacement relations can be simplified by the DMV theory [29]. As a result, the nonlinear membrane strain–displacement relations of the laminate are written as

### Electric Fields.

where $\Delta \varphi $ is the electric bias voltage between the top and bottom electrodes.

## Constitutive Equations

In this section, a generalized form of constitutive equations of the piezoelectric and elastic layers will be obtained in terms of stress and moment resultants as follows.

### Piezoelectric and Elastic Layers.

*ν*

_{2}, and

*e*

_{31}are the Young's modulus, Poisson's ratio, and piezoelectric constant of the piezoelectric layer, respectively, and

where *D _{z}* is the electric displacement in the thickness direction,

*e*

_{31}is the piezoelectric constant, and

*ε*

_{33}is the dielectric constant.

*i*th layer, let us define an in-plane stiffness matrix $K\u0302i$ and bending rigidity matrix $D\u0302i$ as

*i*th layer respectively. With Eq. (18), let us also define $K\u0303i$ as a layer stiffness matrix for the

*i*th layer, consisting of

*i*th layer defined in Eq. (8). For the piezoelectric layer

*e*and

*m*indicate electrical and mechanical components, respectively. Finally, $fe=fe\u2212\u222afe+$ is the union of the electrical stress resultants of the inner and outer electrode domains, which individually can be expressed in terms of $\Delta \varphi $ as follows:

To obtain $fe+$ for the outer electrode, one simply needs to replace the superscript − in Eq. (22) with +.

### Laminate Boundary Forces and Moments.

*V*

^{(}

^{i}^{)}of each layer can be completely determined by the in-pane stress and moment resultants

*N*

^{(}

^{i}^{)}and

*M*

^{(}

^{i}^{)}defined in Eqs. (20) and (21) because the transverse strains

*S*,

_{xz}*S*, and

_{yz}*S*are assumed to be zero [34]. In other words

_{zz}*T*defined in Eq. (9) is used to transform individual stress and moment resultants

_{i}*N*

^{(}

^{i}^{)}and

*M*

^{(}

^{i}^{)}, and shear forces

*V*

^{(}

^{i}^{)}to the global coordinate system. Finally, the shear forces and stress and moment resultants of each layer are summed over all three layers with respect to the modulus-weighted midsurface $S\xaf$ to obtain the total boundary forces

*N*, moments

*M*, and shear forces

*V*at the laminate level. As a result,

where $\alpha ,\beta =x,y$ and $\alpha \u2260\beta $.

## Variation of Potential and Kinetic Energies and Virtual Work

### Variation of Electric Enthalpy and Elastic Strain Energy.

*δH*is [51]

*α*and

*β*. (Note that the superscript (2) is dropped in Eq. (29), because only the piezoelectric layer has electric enthalpy density.) For the elastic layers $S(1)$ and $S(3)$, the variation of the elastic strain energy density

*δU*is

*δH*, $\delta U(1)$, and $\delta U(3)$ are summed together to obtain the potential energy of the laminated shell. In other words

**T**

*in Eq. (9) to each layer stiffness matrix $K\u0303i$ defined in Eq. (19), and summing over all three layers, i.e.,*

_{i}After a variational process, Eq. (31) will lead to a set of electromechanical equations of mechanical motion that governs the displacement vector $u$ and a set of charge equations of electrostatics that governs the electric potential $\Delta \varphi $. Since the shell laminate is an actuator, the electric field is prescribed. As a result, the electric potential $\Delta \varphi $ is no longer a degree-of-freedom (DOF) but a control variable, rendering $\delta \Delta \varphi =0$. Therefore, only the first two integrals in Eq. (31) will be considered throughout the rest of the paper.

### Variation of Kinetic Energy.

*i*th layer can be written as

*ρ*is the mass density of the

_{i}*i*th layer. In Eq. (33), the operator $(\u02d9)$ denotes the time derivative. It follows that the total kinetic energy of the shell model is simply the summation of the kinetic energy of all three layers, i.e.,

### Virtual Work Along Boundaries.

*n*−

*t*coordinate system at the outer boundary $\u2202S\xaf$ (see Fig. 4), where

*n*and

*t*denote the outward normal and tangential direction at a boundary. For the rectangular domain $\Omega ={(x,y)|0\u2264x\u2264a,0\u2264y\u2264b}$ considered in this paper,

**u**

*can be written as [51]*

_{b}By the same token, let us define $ub\u2212$ and $ub+$ as the boundary displacement and slope components in the individual *n* − *t* coordinate systems at *∂*Ω^{−} and *∂*Ω^{+}, respectively. Comparing Figs. 3 and 4, it can be seen that the inner electrode boundary *∂*Ω^{−} has the same normal and tangential direction as $\u2202S\xaf$, and thus $ub\u2212$ shares the same sign convention as **u*** _{b}*. Therefore, to obtain $ub\u2212$ at $x=xle,xre$ and $y=yle,yre$, one simply needs to replace the components in Eq. (35) with equivalent counterparts of superscripts −. The normal and tangential directions along

*∂*Ω

^{+}, on the other hand, have an opposite sign. Therefore, $ub+=\u2212ub\u2212$ except for $w+=w\u2212$ since the transverse direction is the same.

*n*−

*t*coordinate system along $\u2202S\xaf$. Then, according to Ref. [51]

By the same token, let us define $fb\u2212$ and $fb+$ as the boundary force and moment components in the individual *n* − *t* coordinate systems at *∂*Ω^{−} and *∂*Ω^{+}, respectively. Again referred to Fig. 4, the boundary forces and moments components along *∂*Ω^{−} and *∂*Ω^{+} are the same in their individual local *n* − *t* coordinate systems. Therefore, $fb+=fb\u2212$ except for $Vxz+=\u2212Vxz\u2212$ and $Vyz+=\u2212Vyz\u2212$.

The virtual work along *∂*Ω^{−} or *∂*Ω^{+} can be obtained by simply replacing $\u2202S\xaf$, **u*** _{b}*, and

**f**

*in Eq. (37) with*

_{b}*∂*Ω

^{−}or $\u2202\Omega +,\u2009ub\u2212$ or $ub+$, and $fb\u2212$ or $fb+$.

## Variational Formulation and Discretization

*t*) for the piezoelectric shell laminate leads to

*δT*and

*δV*are the variation of the kinetic and potential energies of the laminated shell (cf. Eqs. (31) and (34)), and superscripts + and − indicate quantities at the boundary

*∂*Ω

^{+}and

*∂*Ω

^{−}Moreover, the nonlinear operator $L[u,w\u02d9\Delta \varphi ]$ governs the motion of the shell laminate. Ideally, when the Hamilton principle $\u222bt0t1(\delta T\u2212\delta V)dt=0$ is applied, Eq. (38) leads to a set of equations of motion with boundary conditions

*∂*Ω

^{−}coincide with

*∂*Ω

^{+}, a set of continuity equations (in displacements and slopes) and equilibrium equations (in forces and moments) need to be additionally satisfied

where *N* is the number of terms retained, $\psi k\u2212(\Omega \u2212)$ and $\psi k+(\Omega +)$ are two sets of shape functions defined on Ω^{−} and Ω^{−}, respectively, and $qk\u2212(t)$ and $qk+(t)$ are the generalized coordinates.

*∂*Ω

^{−}). As an alternative, one may satisfy (40) in a

*weak form*, i.e.,

*ψ*serve to satisfy the essential boundary conditions

_{k}**u**

*at $\u2202S\xaf$ and the particular function*

_{b}*ψ*to satisfy the weak form in Eq. (42). The particular function

^{e}*ψ*is similar to a particular solution that is used to satisfy inhomogeneous essential boundary conditions in the Rayleigh-Ritz method (see Ref. [53]). For Donnell's and Novozhilov's nonlinear elastic shell models, such a function can be found in Amabili's work [54], which was added to admissible displacement fields in order to exactly satisfy the classical simply supported boundary conditions. However, to the best of the authors' knowledge, such a function to satisfy Eq. (42) for piezoelectric shell laminates is still missing. Later, in Sec. 8, we will present such a particular function that partially satisfies Eq. (42).

^{e}## Shape Function Expansion

*ψ*are orthogonal, admissible functions that satisfy the essential boundary conditions and orthogonality $\u222b0a\varphi i\varphi j=\delta ij\u222b0a\varphi i\varphi j$ and $\u222b0b\psi i\psi j=\delta ij\u222b0b\psi i\psi j$, where

_{j}*δ*is the Kronecker delta. Furthermore, $\varphi i$ and

_{ij}*ψ*are assumed to be sufficiently smooth such that their derivatives are continuous functions (so are the corresponding strains). The above assumptions are well justified by ample research results of elastic shell studies. For example, Dickinson and Di Blasio [55] proposed a recursive formula based on the Gram–Schmidt process to generate successive orthogonal polynomials that satisfy classical plate boundary conditions, such as SS-SS (simply supported on all sides) and C-C (clamped on all sides). Also, $\mu (x)$ and $\theta (y)$ are assumed to be stepwise functions that satisfy

_{j}and $pje$ and $qie$ are their generalized coordinates.

*u*and electric potential $\Delta \varphi $ as follows:

^{e}where $K=\u2211i=13Ki$. Note that the integrals are evaluated at $x=xre$ and that *u _{m}*,

*v*,

_{m}*w*,

*β*,

_{y}*κ*,

_{xx}*κ*, and

_{yy}*κ*do not appear in Eq. (47) because they are assumed to be continuous across $x=xre$. Also, in deriving Eq. (47), Eq. (2) has been used to eliminate the membrane strains $Sxx0$ from the last integral.

_{xy}Determination of the final forms of $\mu (x)$ and $\theta (y)$ depends on the boundary condition at $\u2202S\xaf$, which will be demonstrated in Sec. 9.

## A Reference System

*M*is the normal moment defined in Eq. (36). Based on these boundary conditions, the shape functions $\varphi i(x)=sin(i\pi x/a)$ and $\psi j(y)=sin(j\pi y/b)$ are assumed. Furthermore, with the boundary conditions $u=v=0$ at $x=0,a$ and $y=0,b,\u2009\mu (x)$ and $\theta (y)$ can be determined from Eq. (46) as

_{nn}As a final remark, the reference system is not chosen to model the microactuator that yields the result in Fig. 2. Due to uneven etch rates, exact diaphragm geometry of the microactuator is rather complex. It has a uniform circular membrane area at the center but a slightly tapered residual silicon region between the circular membrane and the rectangular anchor; see Fig. 6 in Ref. [49]. The piezoelectric thin film is also not uniform. The inner electrode is not rectangular because it needs to extend to outside for electric connection; the outer electrode is not a rectangular ring, because it needs to give way to the inner electrode connection; see Fig. 9 in Ref. [33]. All things considered, it is impractical to model the real microactuator tested. Instead, the reference system is chosen to study snap through for two reasons. First, elastic shells on rectangular domains with simply supported boundary conditions are better studied for snap-through phenomena. Second, shape functions $\psi k(\Omega )$ to discretize the equation of motion are easier to select.

## Discretized Equation of Motion

where $Mij=m\u2211k=1I\u2211l=1J\u222b\Omega \varphi i\psi j\varphi k\psi ldA$ and *m* is the mass per unit area of the shell laminate. Definition of the nonlinear operators $Lpij[u,v,w,\Delta \varphi ],\u2009Lqij[u,v,w,\Delta \varphi ]$_{,} and $L\u0303rij[u,v,w,\Delta \varphi ]$ is provided in Appendix.

*p*and

_{kl}*q*in $Lpij[u,v,w,\Delta \varphi ]=0$ and $Lqij[u,v,w,\Delta \varphi ]=0$ linearly depend on each other; see Eqs. (A1), (A2) and (A15) in the Appendix. Thus, $Lpij[u,v,w,\Delta \varphi ]=0$ and $Lqij[u,v,w,\Delta \varphi ]=0$ can be solved for

_{kl}*p*and

_{kl}*q*. Substitution of the solutions into $L\u0303rij[u,v,w,\Delta \varphi ]=0$ can eliminate

_{kl}*p*and

_{kl}*q*from it, leading to a set of equations consisting of only $r\u2032s$. As a result, the definitive equations of motion are rewritten as

_{kl}where *K* is the linear stiffness, *Q* is the quadratic stiffness, Γ is the cubic stiffness, *S*^{nle} is the nonlinear softening or stiffening coefficient (depending on the signs of $\Delta \varphi \u2212$ and $\Delta \varphi +$) due to the electric potential, *F* is the electric force coefficient due to the converse piezoelectric effect, and *B* is the electric force coefficient due to boundary conditions. The solution procedure to derive Eq. (58) and definition of all coefficients in Eq. (58) can be found in Appendix.

There are two things worth noting in Eq. (58). First, due to the von Karman-type strain displacement relation, the converse piezoelectric effect leads to nonlinear softening or stiffening effect as manifested by the coefficients $Snle\u2212$ and $Snle+$. In other words, an applied voltage can change the stiffness of the shell laminate. Furthermore, these coefficients would vanish if the elastic properties of the piezoelectric layer had been ignored, which was assumed in many studies for piezoelectric smart structures, e.g., in Refs. [44], [46], and [47]. Second, it is also worth noting that the presence of the inner and outer electrodes has a similar effect with $Snle\u2212$ and $Snle+$, as manifested by the coefficients *S*^{nle}. As seen in Eq. (58), if the electric potentials on the inner and outer electrodes are equal, or $\Delta \varphi \u2212=\Delta \varphi +$, i.e., the shell laminate is uniformly actuated, the *S*^{nle} terms will then vanish. In other words, the coefficient *S*^{nle} indicates a specialty of nonuniform or distributed actuation of piezoelectric shells or plates, which largely remains unexplored.

## Conclusions

In this paper, we reported the experimental findings of a PZT thin-film microactuator that was submerged in water. We have also developed a mathematical model hoping to qualitatively explain the experimental findings. With the presentations above, we reach the following conclusions.

- (1)
Experimental results indicate that natural frequencies of the microactuator bifurcate when the applied DC bias voltage reaches specific thresholds. The bifurcation indicates a change of local stiffness and thus a migration between stable equilibrium states of the microactuator. Furthermore, the bifurcation is irreversible, unlike those seen in electrostatic actuators or elastic curve structures. A snap-through phenomenon specific to nonlinear piezoelectric actuators is hypothesized to facilitate the migration when the threshold DC bias voltage is reached.

- (2)
To explore the possibility of snap-through induced by a DC bias voltage, we develop a multilayered, asymmetrically laminated, doubly curved, geometrically nonlinear, shallow shell model. The model is applied to a rectangular domain with simply boundary conditions and two partial electrodes, specifically, an inner electrode and an outer electrode. Discretized equations of motion are derived via a variational formulation and assumed shapes of displacement fields.

- (3)
The geometrical nonlinearity couples in-plane displacements with transverse displacement in the derived equations of motion. In-plane equation of motion is algebraic, implying that in-plane DOFs can be represented in terms of transverse DOFs. Transverse equations of motion are differential equations involving nonlinear terms from both in-plane and transverse DOFs. After elimination of the in-plane DOFs, the resulting equations of motion are coupled second-order differential equations with linear, quadratic, and cubic nonlinear restoring forces.

- (4)
The presence of partial electrodes will lead to a softening or stiffening effect that is similar to what can be observed in nonlinear piezoelectric actuators involved with large deformation. Also, it will lead to a discontinuous electric field, if the applied voltage on the two partial electrodes is not the same. The discontinuity will lead to a discontinuous in-plane strain field across the boundary of the electrodes. The strain discontinuity must be compensated via specific functions $ue(x,y,t)$ in the displacement expansion to ensure proper convergence of the potential energy. The discontinuity will also lead to residual virtual work along the electrode boundaries that needs to be incorporated in the derivation of the discretized equation of motion.

## Acknowledgment

This material is based upon work supported by the National Science Foundation under Grant No. CBET-1159623. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

## Funding Data

National Science Foundation (Grant No. CBET-1159623).

### Appendix

where *A* to *F* are linear stiffness coefficients, $A1nl$ to $A3nl$ and $B1nl$ to $B2nl$ are nonlinear quadratic stiffness coefficients, *C*^{nl} are nonlinear cubic stiffness coefficients, $A1e$ to $A3e$ are linear electric force coefficients, $A1nle$ are nonlinear electric force coefficients due to the particular shape functions defined in Eq. (45), $A1e+/\u2212$ to $A3e+/\u2212$ are linear electric force coefficients and $A2nle+/\u2212$ nonlinear electric force coefficients in the outer and inner electrode domains Ω^{+} and Ω^{–} due to the converse piezoelectric effect, and finally *B ^{e}* and $Be+$ are linear electric force coefficients due to the virtual work at

*∂*Ω

^{+}and $\u2202S\xaf$ defined in Eqs. (54) and (55).

*ij*and

*kl*. Note that

*A*=

_{ijkl}*A*,

_{klij}*B*=

_{ijkl}*B*, and

_{klij}*C*=

_{ijkl}*C*. Furthermore,

_{klij}*D*=

_{ijkl}*D*,

_{klij}*E*=

_{ijkl}*E*, and

_{klij}*F*=

_{ijkl}*F*only when the shell laminate is symmetric in the

_{klij}*xy*plane. Also,

To obtain $A2ijklmnnl$ and $B2ijklmnnl$, one simply needs to substitute $nqijl$ for $npijl$. Similarly, one simply needs to substitute $nqijl$ for $npijl$ to obtain $A2ijkle$ and $A2ije+$, and substitute $nrijl$ to obtain $A3ijkle$ and $A3ije+$. Furthermore, one needs to substitute Ω^{−} for Ω^{+} to obtain $A1ije\u2212,\u2009A2ije\u2212,\u2009A3ije\u2212$ and $A2ijklnle\u2212$.

*F*, $A2ijklmnnl,\u2009A2ijkle,\u2009A2ije+$ and $A2ije\u2212$ for

_{ijkl}*E*, $A1ijklmnnl,\u2009A1ijkle,\u2009A1ije+$ and $A1ije\u2212$, respectively. Assume there exists an inverse matrix (typically obtained by numerical methods) such that

_{ijkl}*p*and

_{kl}*q*in Eq. (A17) into Eq. (A3),

_{kl}*p*and

_{kl}*q*can be eliminated from Eq. (A3), leading to Eq. (58). Finally, coefficients in Eq. (58) can be found as follows:

_{kl}Note that $Fije\u2212$ and $Sijklnle\u2212$ can be obtained by substituting Ω^{−} for Ω^{+}.