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 [117]. 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 [1820], microvalves and microrelays [2123], 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 [2527]) 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,3443], 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 [3741,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,3538,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 ith layer are denoted by a superscript (i) or a subscript i. Furthermore, a special notation (¯) is reserved for a reference surface, known as the modulus-weighted midsurface S¯, 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. (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 Rx and Ry (spherical type). The constant curvatures are used to account for the warping.

  2. (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 Ω={(x,y)|0xa,0yb} and enclosed by an outer boundary S¯, where a and b represent the domain widths. By the same token, Ω is defined as Ω={(x,y)|xlexxre,yleyyre}, and enclosed by an inner boundary Ω, where xle,xre,yle and yre represent the four vertices of the domain. The outer electrode domain is then Ω+=ΩΩ, and is enclosed by S¯ and an inner boundary Ω+, which coincides with Ω.

  3. (3)

    For each layer, there exists a local coordinate system C(i) with orthogonal curvilinear coordinates (α1,α2,α3(i)); see Fig. 3(b). Furthermore, all three layers are shallow shells satisfying the following assumptions [29].

    • (a)

      The orthogonal curvilinear surface coordinates (α1,α2,α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α(1)=Rα(2)=Rα(3), where α=x,y. As a result, the differential volume of each layer can be approximated by dV(i)=dxdydz(i).

  4. (4)
    A global coordinate system C¯ with coordinates (x,y,z¯) is defined, where the x and y coordinates are identical to those in the local coordinate systems C(i), while the z¯ coordinate is obtained by shifting the transverse coordinate z(i); see Fig. 3(b). C¯ is used to determine the modulus-weighted midsurface S¯, 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¯ can be determined by a standard moment analysis. Let z¯(i) be the distance from the midsurface of the ith layer to C¯. Then 
    z¯(0)=i3Kiz¯(i)i3Ki
    (1)
    where z¯(0) is the vertical distance from S¯ to C¯,Ki=(Y(i)hi/(1νi2)) is the conventional membrane stiffness, and Y(i), hi and νi are the Young's modulus, thickness, and Poisson's ratio of the ith layer, respectively. It is convenient to set the origin of the global coordinate system C¯ to coincide with a reference point on S¯ such that the distance z¯(0) becomes zero. Consequently, 
    i3Kiz¯(i)=0
    (2)
Moreover, the coordinate transformation between C(i) and C¯ can be written as 
z(i)=z¯z¯(i)
(3)
  1. (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.

  2. (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 β=(βx,βy,βz)T be the angular rotation vector, where βx and βy are the angles of rotation of the laminate, and βz are the normal deformations of the laminate across the thickness. Consider an arbitrary point P in the laminate, located at (x,y,z¯) in the global coordinate system. Let U be the displacement vector of P. Then 
      U(x,y,z¯)=u(x,y)+z¯β(x,y)
      (4)
    • (b)
      Love's assumptions imply that β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., 
      βα=wα,καβ=c2wαβ
      (5)
      where α,β=x,y, and c = 2 if αβ 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 as 
      Sαβ=Sαβ0+z¯καβ
      (6)
      where α,β=x,y, S are the laminate's in-plane strains and S0 are the membrane strains of the modulus-weighted mid-surface. To transform the laminate's in-plane strains to the ith layer's local coordinate system, Eq. (3) is substituted into Eq. (6), leading to 
      Sαβ(i)=Sαβ0(i)+z(i)καβ
      (7)
      where α,β=x,y and S(i) and S0(i) are the ith layer's in-plane strains and membrane strains of the midsurface, respectively. Finally, the normal stresses σ33 are negligible in all three layers.
  3. (7)

    The electric potential ϕ(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 Ex and Ey are neglected and only the transverse electric field Ez is considered.

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 ith layer (layer level). The shell laminate is then obtained by assembling individual layers and transforming all the quantities to the modulus-weighted midsurface S¯ in the global coordinate system C¯ (laminate level).

Kinematics and Electric Fields

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

Kinematic Relations.

After substituting the coordinate transformation (3) and z(i)=0 into Eq. (6), the membrane strains and curvatures of the ith layer are written in the global coordinate system C¯ as 
s(i)=Tis
(8)
where s(i)=(Sxx0(i),Syy0(i),Sxy0(i),κxx,κyy,κxy)T and s=(Sxx0,Syy0,Sxy0,κxx,κyy,κxy)T are the strain vectors of the ith layer and the laminate, respectively, and Ti is the coordinate transformation matrix from C¯ to C(i), defined as 
Ti=[Iz¯(i)I0I]
(9)

where I and 0 are a 3 × 3 identity and zero matrix, respectively.

Tzou and Bao [34] derived von Karman-type nonlinear strain–displacement relations for thin piezoelectric shells with large transverse deflections 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 
s=[Bl+Bnl]u
(10)
In Eq. (10) 
Bl=[()x0()y0000()y()x000()Rx()Ry02()x22()y22()xy]T
(11)
and 
Bnl=[00000000000012()x()x12()y()y()x()y000]T
(12)
are linear and nonlinear differential operators, respectively.

Electric Fields.

Since the top and bottom layers include electrodes, the transverse electric field appears only inside the piezoelectric layer S(2). Furthermore, the top layer S(3) consists of an inner electrode and an outer electrode (cf. Fig. 1 or Fig. 3(a)). Therefore, the electric potential ϕ(x,y,z(2),t) in layer S(2) is a union of the electric potentials established by the inner and outer electrodes, i.e., 
ϕ(x,y,z(2),t)=ϕ(x,y,z(2),t)ϕ+(x,y,z(2),t)
(13)
where the superscripts − and + refer to the inner and outer electrodes, respectively. Based on Assumption 7, the electric potential varies linearly across the thickness direction. Therefore, the transverse electric field is 
Ez(2)=[Δϕh2Δϕ+h2]
(14)

where Δϕ 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.

According to Ref. [51] and Assumption 6(b) and Assumption 7, the constitutive equations for the piezoelectric layer with z(2) being the poling direction take the form of 
(σxx(2)σyy(2)σxy(2))=[Y(2)1ν22ν2Y(2)1ν220ν2Y(2)1ν22Y(2)1ν22000Y(2)2(1+ν2)](Sxx(2)Syy(2)Sxy(2))(e31Ez(2)e31Ez(2)0)
(15)
where Y(2), ν2, and e31 are the Young's modulus, Poisson's ratio, and piezoelectric constant of the piezoelectric layer, respectively, and 
Dz=e31(Sxx(2)+Syy(2))+ϵ33Ez(2)
(16)

where Dz is the electric displacement in the thickness direction, e31 is the piezoelectric constant, and ε33 is the dielectric constant.

The elastic layers S(1) and S(3) are assumed to be isotropic. Consequently, the constitutive equations take the form of 
(σxx(i)σyy(i)σxy(i))=[Y(i)1νi2νiY(i)1νi20νiY(i)1νi2Y(i)1νi2000Y(i)2(1+νi)](Sxx(i)Syy(i)Sxy(i)),i=1,3
(17)
For the ith layer, let us define an in-plane stiffness matrix K̂i and bending rigidity matrix D̂i as 
K̂i=[KiνiKi0νiKiKi000(1νi)Ki2],D̂i=[DiνiDi0νiDiDi000(1νi)Di2]
(18)
where Ki=(Y(i)hi/(1νi2)) and Di=(Y(i)hi3/12(1νi2)) are the conventional membrane stiffness and bending rigidity of the ith layer respectively. With Eq. (18), let us also define K̃i as a layer stiffness matrix for the ith layer, consisting of 
K̃i=[K̂i00D̂i]
(19)
After substituting Eq. (7) into the constitutive equations (15) and (17) and integrating them from hi/2 to hi/2 in the local coordinate system (x,y,z(i)), the constitutive equations can be expressed in terms of stress and moment resultants, and the layer stiffness matrix (19). For the elastic layers 
fm(i)=K̃is(i),i=1,3
(20)
where fm(i)=(Nxxm(i),Nyym(i),Nxym(i),Mxxm(i),Myym(i),Mxym(i))T is a force (per unit length) vector consisting of the mechanical stress and moment resultants of the elastic layers, and s(i) is the strain vector of the ith layer defined in Eq. (8). For the piezoelectric layer 
f(2)=K̃2s(2)fe(2)=fm(2)fe(2)
(21)
where f(2)=(Nxx(2),Nyy(2),Nxy(2),Mxx(2),Myy(2),Mxy(2))T is a force (per unit length) vector consisting of the stress and moment resultants of the piezoelectric layer. Note that f(2) is composed of fm(2), the mechanical stress and moment resultants obtained under short-circuit condition, and fe=(Nxxe,Nyye,0,0,0,0)T, the electrical stress resultants. Note that the electric moment resultants Mxxe=Myye=0 because the electric field (14) is symmetric with respect to the midsurface of the piezoelectric layer and the integration was carried out from h2/2 to h2/2. Also, note that the superscripts e and m indicate electrical and mechanical components, respectively. Finally, fe=fefe+ is the union of the electrical stress resultants of the inner and outer electrode domains, which individually can be expressed in terms of Δϕ as follows: 
fe=e31(Δϕ,Δϕ,0,0,0,0)T
(22)

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

Laminate Boundary Forces and Moments.

The boundary forces and moments for the entire shell laminate are derived as follows: First, in the current shell formulation, transverse shear forces 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 Sxz, Syz, and Szz are assumed to be zero [34]. In other words 
Vαz(i)=Vαzm(i)Vαze(i)=[Qαzm(i)+Mαβm(i)βNααm(i)βαNαβm(i)ββ][Nααe(i)βα]
(23)
where α,β=x,y and αβ, and Qm(i) are defined as 
Qαzm(i)=Mααm(i)α+Mαβm(i)β
(24)
where α,β=x,y and αβ. Second, the coordinate transformation matrix Ti defined in Eq. (9) is used to transform individual stress and moment resultants 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¯ to obtain the total boundary forces N, moments M, and shear forces V at the laminate level. As a result, 
Nαβ=i=13Nαβm(i)Nαβe(2)
(25)
where α,β=x,y. Note that Nαβe(2)=0 when αβ. Furthermore, 
Mαα=i=13[Mααm(i)+Nααm(i)z¯(i)]Nααe(2)z¯(2)
(26)
where α=x,y and Nz¯ is the moment created by the in-plane stress resultants due to the moment arm z¯. Moreover 
Vαz=i=13Vαzm(i)Vαze(2)
(27)
where α=x,y. Finally, 
Vαβm=i=13[Nαβm(i)+Mαβm(i)Rβ]
(28)

where α,β=x,y and αβ.

Variation of Potential and Kinetic Energies and Virtual Work

Variation of Electric Enthalpy and Elastic Strain Energy.

For the piezoelectric layer S(2), the variation of the electric enthalpy δH is [51] 
δH=V[σαβδSαβDzδEz]dV,α,β=x,y
(29)
where Einstein's summation notation is used for α 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 
δU(i)=V(i)σαβ(i)δSαβ(i)dV(i),i=1,3
(30)
The variation of the total potential energy of the shell laminate is derived as follows: First, the constitutive equations (15)(17) are substituted into Eqs. (29) and (30) and integrated in the thickness direction from hi/2 to hi/2 in the individual local coordinate systems. Second, the constitutive equations in the layer level (cf. Eqs. (20) and (21)) are used to express the potential energy in terms of the layer strain vector s(i). Third, the potential energy is transformed to the global coordinate system and thus expressed in terms of s, the laminate strain vector of the shell laminate, via Eq. (8). Fourth, the nonlinear strain–displacement relations (10) are used to express the potential energy in terms of the displacement vector u of the modulus-weighted midsurface. Finally, δH, δU(1), and δU(3) are summed together to obtain the potential energy of the laminated shell. In other words 
δV[u,Δϕ]=ΩδuT(BlT+BnlT)K(Bl+Bnl)udAΩδuT(BlT+BnlT)T2TfedA+δΔϕe31Ω[1TT2(Bl+Bnl)u]dAδΔϕCPΔϕ
(31)
where 1=(1,1,0,0,0,0)T,Cp=(ε33V(2)/h22) is the equivalent capacitance of the piezoelectric layer, and K is the laminate stiffness matrix, which is obtained by applying coordinate transformation Ti in Eq. (9) to each layer stiffness matrix K̃i defined in Eq. (19), and summing over all three layers, i.e., 
K=i=13TiTK̃iTi=i=13[K̂iz¯(i)K̂iz¯(i)K̂i(z¯(i))2K̂i+D̂i]
(32)

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 Δϕ. Since the shell laminate is an actuator, the electric field is prescribed. As a result, the electric potential Δϕ is no longer a degree-of-freedom (DOF) but a control variable, rendering δΔϕ=0. Therefore, only the first two integrals in Eq. (31) will be considered throughout the rest of the paper.

Variation of Kinetic Energy.

When formulating shallow shells, in-plane and rotary inertia effects are often neglected [28,45,52]. The negligible rotary inertia can be accounted for by the DMV theory [29]. Furthermore, Tzuo and Bao [34] asserted that the negligible in-plane inertia is due to the DMV theory. As a result, the kinetic energies of the ith layer can be written as 
δT(i)=ρihiΩδw˙(i)w˙(i)dA
(33)
where ρi is the mass density of the ith layer. In Eq. (33), the operator (˙) 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., 
δT[w]=i=13δT(i)
(34)

Virtual Work Along Boundaries.

Let us define ub=(un,ut,w,βn)T as a vector consisting of the displacement and slope components in a local n − t coordinate system at the outer boundary S¯ (see Fig. 4), where n and t denote the outward normal and tangential direction at a boundary. For the rectangular domain Ω={(x,y)|0xa,0yb} considered in this paper, ub can be written as [51] 
x=0,a:ub=(u,v,w,βx)T,(u,v,w,βx)Ty=0,b:ub=(v,u,w,βy)T,(v,u,w,βy)T
(35)

By the same token, let us define ub 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 S¯, and thus ub shares the same sign convention as ub. Therefore, to obtain ub 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+=ub except for w+=w since the transverse direction is the same.

Similarly, let us define fb=fbmfbe=(Nnnm,Nntm,Vnm,Mnnm)T(Nnne,Nnte,Vne,Mnne)T as a vector consisting of the mechanical and electrical force and moment components in the local n − t coordinate system along S¯. Then, according to Ref. [51] 
x=0,a:fb=(Nxx,Vxym,Vxz,Mxx)T,(Nxx,Vxym,Vxz,Mxx)Ty=0,b:fb=(Nyy,Vyxm,Vyz,Myy)T,(Nyy,Vyxm,Vyz,Myy)T
(36)

By the same token, let us define fb 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 except for Vxz+=Vxz and Vyz+=Vyz.

Finally, the virtual work along S¯ is written as 
δW=S¯δubTfb
(37)

The virtual work along Ω or Ω+ can be obtained by simply replacing S¯, ub, and fb in Eq. (37) with Ω or Ω+,ub or ub+, and fb or fb+.

Variational Formulation and Discretization

Assuming that no external forces are present and the electric potential Δϕ is prescribed, a variational formulation (with respect to time t) for the piezoelectric shell laminate leads to 
δTδV=ΩδuTL[u,w˙,Δϕ]dAΩ(δub)T(fbfb+)dSS¯δubTfbdS
(38)
In Eq. (38), δ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˙Δϕ] governs the motion of the shell laminate. Ideally, when the Hamilton principle t0t1(δTδV)dt=0 is applied, Eq. (38) leads to a set of equations of motion with boundary conditions 
L[u,w˙,Δϕ]=0onΩ,ub=0orfb=0atS¯
(39)
Furthermore, since Ω coincide with Ω+, a set of continuity equations (in displacements and slopes) and equilibrium equations (in forces and moments) need to be additionally satisfied 
ub=ub+,expectforw+=w,fbfb+=0,expectforVn=Vn+atΩ
(40)
To discretize (39), for example, using Galerkin's method, one could assume 
u=k=1Nψk(Ω)qk(t),onΩ;u+=k=1Nψk+(Ω+)qk+(t),onΩ+
(41)

where N is the number of terms retained, ψk(Ω) and ψk+(Ω+) are two sets of shape functions defined on Ω and Ω, respectively, and qk(t) and qk+(t) are the generalized coordinates.

In the discretization process, a major challenge is to choose proper shape functions that satisfy the continuity and equilibrium equations defined in Eq. (40). Finding such ψk and ψk+ that exactly satisfy Eq. (40) can be formidable because it involves 16 nonlinear algebraic equations (fb+=fb has four equations on each side of Ω). As an alternative, one may satisfy (40) in a weak form, i.e., 
Ω(δub)T(fbfb+)dS=0
(42)
and thus resort to the Rayleigh–Ritz method (lieu of Galerkin's method) by seeking a set of admissible functions and particular functions 
u=k=1Nψk(Ω)qk(t)+ψe(Ω)qe(t)
(43)
to discretize (38) over the entire laminate domain Ω. The admissible functions ψk serve to satisfy the essential boundary conditions ub at S¯ and the particular function ψe 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).
When using admissible functions to perform discretization, the virtual work associated with natural boundary conditions has to be included in the process since they are not satisfied by the admissible functions [53]. As a result, the following variational formulation will be used to derive the equation of motion for the shell laminate 
δTδV+Ω(δub)T(fbfb+)dS+S¯δubTfbdS=0
(44)

Shape Function Expansion

To find a particular function that satisfies the weak form in Eq. (42), the displacement field of the shell laminate is assumed to be of u=um+ue, where 
um=i=1Ij=1Jϕi(x)ψj(y)(pij(t)qij(t)rij(t))ue=i=1Ij=1J(μ(x)ψj(y)pje(t)θ(y)ϕi(x)qie(t)0)
(45)
In Eq. (45), pij(t),qij(t), and rij(t) are the generalized coordinates, and ϕi and ψj are orthogonal, admissible functions that satisfy the essential boundary conditions and orthogonality 0aϕiϕj=δij0aϕiϕj and 0bψiψj=δij0bψiψj, where δij is the Kronecker delta. Furthermore, ϕi and ψj 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, μ(x) and θ(y) are assumed to be stepwise functions that satisfy 
μxμ+x=1,θyθ+y=1
(46)

and pje and qie are their generalized coordinates.

After substituting Eq. (45) into the weak form in Eq. (42), pje and qie can be determined as follows: First, let us focus on the side x=xre. After substituting Eq. (45) into Eqs. (42) through (10), (20), (21), and (25)(28), the virtual work at x=xre can be written in terms of the particular displacement function ue and electric potential Δϕ as follows: 
Ω(δub)T(fbfb+)dS=yleyreδu[K(uexue+x)+e31(ΔϕΔϕ+)]dyyleyreδwβx[K(uexue+x)+e31(ΔϕΔϕ+)]dy+yleyreδβxz¯(2)e31(ΔϕΔϕ+)dy
(47)

where K=i=13Ki. Note that the integrals are evaluated at x=xre and that um, vm, w, βy, κxx, κyy, and κxy 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.

When ΔϕΔϕ+,(ue/x) needs to be a stepwise function that has a discontinuity at x=xre for the first two terms in Eq. (47) to vanish. This is the reason why μ(x) and θ(y) are assumed to be the stepwise functions defined in Eq. (46). Incorporating similar stepwise functions to satisfy the continuity and compatibility equations at the boundaries of partial and concatenated piezoelectric beam laminates can also be found in Refs. [56] and [57], respectively. Substitution of Eq. (46) into Eq. (47) leads to 
Ω(δub)T(fbfb+)dS=i=1Ij,k=1Jδpijϕiyleyre[Kψjψkpke+e31(ΔϕΔϕ+)]dy+i=1Ij,k=1Jδrijϕiyleyreβx[Kψjψkpkee31(ΔϕΔϕ+)]dy+yleyreδβxz¯(2)e31(ΔϕΔϕ+)dy
(48)
Note that the integrals are evaluated at x=xre. Finally, equating the first or second integral in Eq. (48) to zero and making use of orthogonality, one can derive pje as 
pje(t)=e31K[Δϕ(t)Δϕ(t)+]yleyreψj(y)dy0b[ψj(y)]2dy
(49)
The same analysis can be done to y=yre or y=yle to derive qie. In other words 
qie(t)=e31K[Δϕ(t)Δϕ(t)+]xlexreϕi(x)dx0a[ϕi(x)]2dx
(50)

Determination of the final forms of μ(x) and θ(y) depends on the boundary condition at S¯, which will be demonstrated in Sec. 9.

A Reference System

To demonstrate how to discretize Eq. (44), we will consider, for the rest of the paper, a reference system where the shell laminate is subject to a set of simply supported boundary conditions. We will use the reference system to study its snap-through phenomena for the rest of the paper. For the reference system, the simply supported boundary conditions at S¯ are 
u=0,Mnn=0atS¯
(51)
where u=(u,v,w)T is the displacement vector of the modulus-weighted midsurface S¯, and Mnn is the normal moment defined in Eq. (36). Based on these boundary conditions, the shape functions ϕi(x)=sin(iπx/a) and ψj(y)=sin(jπy/b) are assumed. Furthermore, with the boundary conditions u=v=0 at x=0,a and y=0,b,μ(x) and θ(y) can be determined from Eq. (46) as 
μ(x)={(xrexle)xa,0xxle(axre+xle)xaxle,xle<x<xre(xrexle)(xa1),xrexa
(52)
and 
θ(y)={(yreyle)yb,0yyle(byre+yle)ybyle,yle<y<yre(yreyle)(yb1),yreyb
(53)
Provided that the shape functions μ(x) and θ(y) nullify the first two terms in Eq. (47), the weak form of equilibrium equations in Eq. (44) can be reduced to 
Ω(δub)T(fbfb+)dS=[Nxxe(2)Nxxe(2)+]z¯(2)yleyre(δβx|x=xreδβx|x=xle)dy[Nyye(2)Nyye(2)+]z¯(2)xlexre(δβy|y=yreδβy|y=yle)dx
(54)
Also, the simply supported boundary conditions require the virtual work at S¯ in Eq. (38) to be 
S¯δubTfbdS=Nxxe+z¯(2)0b(δβx(2)|x=aδβx(2)|x=0)dyNyye+z¯(2)0a(δβy(2)|y=bδβy(2)|y=0)dx
(55)

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 ψk(Ω) to discretize the equation of motion are easier to select.

Discretized Equation of Motion

After substituting the shape function expansion in Eq. (45) into the variational formulation (44) with the virtual work ((54) and (55)) and carrying out the variation, the following discretized equations of motion can be derived: 
i=1Ij=1JLpij[u,v,w,Δϕ]=0,i=1Ij=1JLqij[u,v,w,Δϕ]=0i=1Ij=1JMijr¨ij+L̃rij[u,v,w,Δϕ]=0
(56)

where Mij=mk=1Il=1JΩϕiψjϕkψldA and m is the mass per unit area of the shell laminate. Definition of the nonlinear operators Lpij[u,v,w,Δϕ],Lqij[u,v,w,Δϕ], and L̃rij[u,v,w,Δϕ] is provided in Appendix.

Note that pkl and qkl in Lpij[u,v,w,Δϕ]=0 and Lqij[u,v,w,Δϕ]=0 linearly depend on each other; see Eqs. (A1), (A2) and (A15) in the Appendix. Thus, Lpij[u,v,w,Δϕ]=0 and Lqij[u,v,w,Δϕ]=0 can be solved for pkl and qkl. Substitution of the solutions into L̃rij[u,v,w,Δϕ]=0 can eliminate pkl and qkl from it, leading to a set of equations consisting of only rs. As a result, the definitive equations of motion are rewritten as 
i=1Ij=1JMijr¨ij+Lrij[u,v,w,Δϕ]=0
(57)
In Eq. (57) 
Lrij[u,v,w,Δϕ]=[k=1Il=1JKijkl+k,m=1Il,n=1JSijklmnnle(ΔϕΔϕ+)+k=1Il=1JSijklnle+Δϕ++k=1Il=1JSijklnleΔϕ]rkl+k,m=1Il,n=1JQijklmnnlrklrmn+k,m,o=1Il,n,p=1JΓijklmnopnlrklrmnrop+k=1Il=1JFijkle(ΔϕΔϕ+)+Fije+Δϕ++FijeΔϕBije(ΔϕΔϕ+)Bije+Δϕ+
(58)

where K is the linear stiffness, Q is the quadratic stiffness, Γ is the cubic stiffness, Snle is the nonlinear softening or stiffening coefficient (depending on the signs of Δϕ and Δϕ+) 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 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 and Snle+, as manifested by the coefficients Snle. As seen in Eq. (58), if the electric potentials on the inner and outer electrodes are equal, or Δϕ=Δϕ+, i.e., the shell laminate is uniformly actuated, the Snle terms will then vanish. In other words, the coefficient Snle 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. (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. (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. (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. (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

In Eq. (56) 
Lpij[u,v,w,Δϕ]=k=1Il=1JAijklpkl+k=1Il=1JDijklqkl+k=1Il=1JEijklrkl+k,m=1Il,n=1JA1ijklmnnlrklrmn+k=1Il=1JA1ijkle(ΔϕΔϕ+)A1ije+Δϕ+A1ijeΔϕ
(A1)
and 
Lqij[u,v,w,Δϕ]=k=1Il=1JDklijpkl+k=1Il=1JBijklqkl+k=1Il=1JFijklrkl+k,m=1Il,n=1JA2ijklmnnlrklrmn+k=1Il=1JA2ijkle(ΔϕΔϕ+)A2ije+Δϕ+A2ijeΔϕ
(A2)
and 
L̃rij[u,v,w,Δϕ]=[k=1Il=1JCijkl+k,m=1Il.n=1JA1ijklmnnle(ΔϕΔϕ+)k=1Il=1JA2ijklnle+Δϕ+k=1Il=1JA2ijklnleΔϕ]rkl+k,m=1Il,n=1JA3ijklmnnlrklrmn+k,m=1Il,n=1JB1ijklmnnlrklpmn+k,m=1Il,n=1JB2ijklmnnlrklqmn+k,m,o=1Il,n,p=1JCijklmnopnlrklrmnrop+k=1Il=1JA3ijkle(ΔϕΔϕ+)A3ije+Δϕ+A3ijeΔϕBije(ΔϕΔϕ+)Bije+Δϕ++k=1Il=1JEklijpkl+k=1Il=1JFklijqkl
(A3)

where A to F are linear stiffness coefficients, A1nl to A3nl and B1nl to B2nl are nonlinear quadratic stiffness coefficients, Cnl 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+/ to A3e+/ are linear electric force coefficients and A2nle+/ nonlinear electric force coefficients in the outer and inner electrode domains Ω+ and Ω due to the converse piezoelectric effect, and finally Be and Be+ are linear electric force coefficients due to the virtual work at Ω+ and S¯ defined in Eqs. (54) and (55).

In deriving Eqs. (A1)(A3), several definitions of vectors are introduced as follows: From Eq. (22) 
fe=e31Δϕ1,fe+=e31Δϕ+1
(A4)
where 1=(1,1,0,0,0,0)T. From Eqs. (11) and (43) 
Blue=e31(ΔϕΔϕ+)Nijecij
(A5)
where 
Nije=[μ,xψj0μψj,y000θ,yϕi0θϕi,x000]T,cij=1K(yleyreψjdy0b[ψj]2dy,xlexreϕidx0a[ϕi]2dx)T
(A6)
From Eqs. (12) and (43) 
Blum=[npijl,nqijl,nrijl](pij,qij,rij)T
(A7)
where 
[npijl,nqijl,nrijl]=[ϕi,xψj0ϕiψj,y0000ϕiψj,yϕi,xψj000ϕiψjRxϕiψjRy0ϕi,xxψjϕiψj,yy2ϕi,xψj,y]T
(A8)
and 
Bnlum=nijklnl(pijpkl,qijqkl,rijrkl)T
(A9)
where 
nijklnl=[12ϕi,xϕk,xψjψl12ϕiϕkψj,yψl,yϕi,xϕkψjψl,y000]T
(A10)
The coefficients in Eqs. (A1)(A3) can be found as follows: 
Aijkl=Ω(npijl)TKnpkll,Bijkl=Ω(nqijl)TKnqkllCijkl=Ω(nrijl)TKnrkll,Dijkl=Ω(npijl)TKnqkllEijkl=Ω(npijl)TKnrkll,Fijkl=Ω(nqijl)TKnrkll
(A11)
To obtain quantities ()klij, one simply needs to interchange ij and kl. Note that Aijkl = Aklij, Bijkl = Bklij, and Cijkl = Cklij. Furthermore, Dijkl = Dklij, Eijkl = Eklij, and Fijkl = Fklij only when the shell laminate is symmetric in the xy plane. Also, 
A1ijklmnnl=Ω(npijl)TKnklmnnlA3ijklmnnl=Ω(nrijl)TKnklmnnl+2Ω(nijklnl)TKnrmnlB1ijklmnnl=2Ω(nijklnl)TKnpmnlCijklmnopnl=2Ω(nijklnl)TKnmnopnl
(A12)
and 
A1ijkle=e31Ω(npijl)TKNklecklA1ijklmnnle=2e31Ω(nijklnl)TKNmnecmnA2ijklnle+=2e31Ω(nijklnl)TT2T1A1ije+=e31Ω+(npijl)TT2T1
(A13)
Furthermore 
Bije=e31z¯(2)(ϕi,x|x=xleϕi,x|x=xre)yleyreψjdy+e31z¯(2)(ψj,y|y=yleψj,y|y=yre)xlexreϕidxBije+=e31z¯(2)(ϕi,x|x=0ϕi,x|x=a)0bψjdy+e31z¯(2)(ψj,y|y=0ψj,y|y=b)0aϕidx
(A14)

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,A2ije,A3ije and A2ijklnle.

Note that pkl and qkl in Eqs. (A1) and (A2) linearly depend on each other, and thus can be rearranged into a matrix form 
[AijklDijklDklijBijkl](pklqkl)=(fpkl(r,Δϕ)fqkl(r,Δϕ))
(A15)
In Eq. (A15) 
fpkl(r,Δϕ)=k=1Il=1JEijklrklk,m=1Il,n=1JA1ijklmnnlrklrmnk=1Il=1JA1ijkle(ΔϕΔϕ+)+A1ije+Δϕ++A1ijeΔϕ
(A16)
and fqkl(r,Δϕ) can be obtained by substituting Fijkl, A2ijklmnnl,A2ijkle,A2ije+ and A2ije for Eijkl, A1ijklmnnl,A1ijkle,A1ije+ and A1ije, respectively. Assume there exists an inverse matrix (typically obtained by numerical methods) such that 
[ÃijklD̃ijklD̃klijB̃ijkl](fpkl(r,Δϕ)fqkl(r,Δϕ))=(pklqkl)
(A17)
Substitution of pkl and qkl in Eq. (A17) into Eq. (A3), pkl and qkl can be eliminated from Eq. (A3), leading to Eq. (58). Finally, coefficients in Eq. (58) can be found as follows: 
Fijkle=A3ijkleEklijm=1In=1J(ÃijmnA1ijmne+D̃ijmnA2ijmne)Fklijm=1In=1J(D̃mnijA1ijmne+B̃ijmnA2ijmne)Fije+=A3ije++k=1Il=1JEklij(ÃijklA1kle++D̃ijklA2kle+)+k=1Il=1JFklij(D̃klijA1kle++B̃ijklA2kle+)+m=1In=1JB2ijklmnnl(D̃mnijA1mne+B̃ijmnA2mne)Sijklnle+=m=1In=1JB1ijklmnnl(ÃijmnA1mne++D̃ijmnA2mne+)+m=1In=1JB2ijklmnnl(D̃mnijA1mne++B̃ijmnA2mne+)A2ijklnle+Sijklmnnle=B1ijklmnnlo=1Ip=1J(ÃijopA1ijope+D̃ijopA2ijope)B2ijklmnnlo=1Ip=1J(D̃opijA1ijope+B̃ijopA2ijope)+A1ijklmnnleKijkl=CijklEklijm=1In=1J(ÃijmnEijmn+D̃ijmnFijmn)Fklijm=1In=1J(D̃mnijEijmn+B̃ijmnFijmn)Qijklmnnl=A3ijklmnnlB1ijklmnnlo=1Ip=1J(ÃijopEijop+D̃ijopFijop)B2ijklmnnlo=1Ip=1J(D̃opijEijop+B̃ijopFijop)Eklijo=1Ip=1J(ÃijopA1ijmnopnl+D̃ijopA2ijmnopnl)Fklijo=1Ip=1J(D̃opijA1ijmnopnl+B̃ijopA2ijmnopnl)Γijklmnopnl=CijklmnopnlB1ijklmnnlq=1Ir=1J(ÃijqrA1ijopqrnl+D̃ijqrA2ijopqrnl)B2ijklmnnlq=1Ir=1J(D̃qrijA1ijopqrnl+B̃ijqrA2ijopqrnl)
(A18)

Note that Fije and Sijklnle can be obtained by substituting Ω for Ω+.

References

References
1.
Krylov
,
S.
, and
Seretensky
,
S.
,
2006
, “
Higher Order Correction of Electrostatic Pressure and Its Influence on the Pull-In Behavior of Microstructures
,”
J. Micromech. Microeng.
,
16
(
7
), pp.
1382
1396
.
2.
Zhang
,
Y.
,
Wang
,
Y.
,
Li
,
Z.
,
Huang
,
Y.
, and
Li
,
D.
,
2007
, “
Snap-Through and Pull-In Instabilities of an Arch-Shaped Beam Under an Electrostatic Loading
,”
J. Microelectromech. Syst.
,
16
(
3
), pp.
684
693
.
3.
Nayfeh
,
A. H.
,
Younis
,
M. I.
, and
Abdel-Rahman
,
E. M.
,
2007
, “
Dynamic Pull-in Phenomenon in Mems Resonators
,”
Nonlinear Dyn.
,
48
(
1–2
), pp.
153
163
.
4.
Krylov
,
S.
,
Ilic
,
B.
,
Schreiber
,
D.
,
Seretensky
,
S.
, and
Craighead
,
H.
,
2008
, “
The Pull-In Behavior of Electrostatically Actuated Bistable Microstructures
,”
J. Micromech. Microeng.
,
18
(
5
), p. 055026.
5.
Batra
,
R. C.
,
Porfiri
,
M.
, and
Spinello
,
D.
,
2008
, “
Effects of Van Der Waals Force and Thermal Stresses on Pull-In Instability of Clamped Rectangular Microplates
,”
Sensors
,
8
(
12
), pp.
1048
1069
.
6.
Das
,
K.
, and
Batra
,
R.
,
2009
, “
Symmetry Breaking, Snap-Through and Pull-In Instabilities Under Dynamic Loading of Microelectromechanical Shallow Arches
,”
Smart Mater. Struct.
,
18
(
11
), p.
115008
.
7.
Das
,
K.
, and
Batra
,
R.
,
2009
, “
Pull-In and Snap-Through Instabilities in Transient Deformations of Microelectromechanical Systems
,”
J. Micromech. Microeng.
,
19
(
3
), p. 035008.
8.
Ouakad
,
H. M.
, and
Younis
,
M. I.
,
2010
, “
The Dynamic Behavior of Mems Arch Resonators Actuated Electrically
,”
Int. J. Non-Linear Mech.
,
45
(
7
), pp.
704
713
.
9.
Medina
,
L.
,
Gilat
,
R.
, and
Krylov
,
S.
,
2014
, “
Symmetry Breaking in an Initially Curved Pre-Stressed Micro Beam Loaded by a Distributed Electrostatic Force
,”
Int. J. Solids Struct.
,
51
(
11–12
), pp.
2047
2061
.
10.
Ouakad
,
H. M.
, and
Younis
,
M. I.
,
2014
, “
On Using the Dynamic Snap-Through Motion of Mems Initially Curved Microbeams for Filtering Applications
,”
J. Sound Vib.
,
333
(
2
), pp.
555
568
.
11.
Medina
,
L.
,
Gilat
,
R.
,
Robert Ilic
,
B.
, and
Krylov
,
S.
,
2016
, “
Experimental Dynamic Trapping of Electrostatically Actuated Bistable Micro-Beams
,”
Appl. Phys. Lett.
,
108
(
7
), p. 073503.
12.
Saghir
,
S.
, and
Younis
,
M. I.
,
2016
, “
An Investigation of the Static and Dynamic Behavior of Electrically Actuated Rectangular Microplates
,”
Int. J. Non-Linear Mech.
,
85
, pp.
81
93
.
13.
Saghir
,
S.
,
Bellaredj
,
M.
,
Ramini
,
A.
, and
Younis
,
M. I.
,
2016
, “
Initially Curved Microplates Under Electrostatic Actuation: Theory and Experiment
,”
J. Micromech. Microeng.
,
26
(
9
), p.
095004
.
14.
Medina
,
L.
,
Gilat
,
R.
, and
Krylov
,
S.
,
2016
, “
Bistable Behavior of Electrostatically Actuated Initially Curved Micro Plate
,”
Sens. Actuators A: Phys.
,
248
, pp.
193
198
.
15.
Jallouli
,
A.
,
Kacem
,
N.
,
Bourbon
,
G.
,
Le Moal
,
P.
,
Walter
,
V.
, and
Lardies
,
J.
,
2016
, “
Pull-In Instability Tuning in Imperfect Nonlinear Circular Microplates Under Electrostatic Actuation
,”
Phys. Lett. A
,
380
(
46
), pp.
3886
3890
.
16.
Medina
,
L.
,
Gilat
,
R.
, and
Krylov
,
S.
,
2017
, “
Modeling Strategies of Electrostatically Actuated Initially Curved Bistable Micro Plates
,”
Int. J. Solids Struct.
,
118–119
, pp.
1
13
.
17.
Li
,
L.
, and
Zhang
,
Q.-C.
,
2017
, “
Nonlinear Dynamic Analysis of Electrically Actuated Viscoelastic Bistable Microbeam System
,”
Nonlinear Dyn.
,
87
(
1
), pp.
587
604
.
18.
Saif
,
M. T. A.
,
2000
, “
On a Tunable Bistable Mems-Theory and Experiment
,”
J. Microelectromech. Syst.
,
9
(
2
), pp.
157
170
.
19.
Han
,
J. S.
,
Ko
,
J. S.
,
Kim
,
Y. T.
, and
Kwak
,
B. M.
,
2002
, “
Parametric Study and Optimization of a Micro-Optical Switch With a Laterally Driven Electromagnetic Microactuator
,”
J. Micromech. Microeng.
,
12
(
6
), p.
939
.
20.
Receveur
,
R. A.
,
Marxer
,
C. R.
,
Woering
,
R.
,
Larik
,
V. C.
, and
de Rooij
,
N.-F.
,
2005
, “
Laterally Moving Bistable Mems DC Switch for Biomedical Applications
,”
J. Microelectromech. Syst.
,
14
(
5
), pp.
1089
1098
.
21.
Wagner
,
B.
,
Quenzer
,
H.
,
Hoerschelmann
,
S.
,
Lisec
,
T.
, and
Juerss
,
M.
,
1996
, “
Bistable Microvalve With Pneumatically Coupled Membranes
,”
Ninth Annual International Workshop on Micro Electro Mechanical Systems, Investigation of Micro Structures, Sensors, Actuators, Machines and Systems
, (
MEMS'96
), San Diego, CA, Feb. 11–15, pp.
384
388
.
22.
Wilcox
,
D. L.
, and
Howell
,
L. L.
,
2005
, “
Fully Compliant Tensural Bistable Micromechanisms (Ftbm)
,”
J. Microelectromech. Syst.
,
14
(
6
), pp.
1223
1235
.
23.
Ko
,
J. S.
,
Lee
,
M. G.
,
Han
,
J. S.
,
Go
,
J. S.
,
Shin
,
B.
, and
Lee
,
D.-S.
,
2006
, “
A Laterally-Driven Bistable Electromagnetic Microrelay
,”
ETRI J.
,
28
(
3
), pp.
389
392
.
24.
Simitses
,
G. J.
,
1990
,
Dynamic Stability of Suddenly Loaded Structures
,
Springer
,
New York
.
25.
Hsu
,
C.
,
1968
, “
Stability of Shallow Arches against Snap-Through Under Timewise Step Loads
,”
ASME J. Appl. Mech.
,
35
(
1
), pp.
31
39
.
26.
Lin
,
J.-S.
, and
Chen
,
J.-S.
,
2003
, “
Dynamic Snap-Through of a Laterally Loaded Arch Under Prescribed End Motion
,”
Int. J. Solids Struct.
,
40
(
18
), pp.
4769
4787
.
27.
Chen
,
J.-S.
, and
Lin
,
J.-S.
,
2004
, “
Dynamic Snap-Through of a Shallow Arch Under a Moving Point Load
,”
ASME J. Vib. Acoust.
,
126
(
4
), pp.
514
519
.
28.
Ye
,
Z.
,
1997
, “
The Non-Linear Vibration and Dynamics Instability of Thin Shallow Shells
,”
J. Sound Vib.
,
202
(
3
), pp.
303
311
.
29.
Ventsel
,
E.
, and
Krauthammer
,
T.
,
2001
,
Thin Plates and Shells: Theory: Analysis, and Applications
,
CRC Press
, Boca Raton, FL.
30.
Maurini
,
C.
,
Pouget
,
J.
, and
Vidoli
,
S.
,
2007
, “
Distributed Piezoelectric Actuation of a Bistable Buckled Beam
,”
Eur. J. Mech.-A/Solids
,
26
(
5
), pp.
837
853
.
31.
Varelis
,
D.
, and
Saravanos
,
D. A.
,
2006
, “
Coupled Mechanics and Finite Element for Non-Linear Laminated Piezoelectric Shallow Shells Undergoing Large Displacements and Rotations
,”
Int. J. Numer. Methods Eng.
,
66
(
8
), pp.
1211
1233
.
32.
Luo
,
C.
,
Cao
,
G.
, and
Shen
,
I.
,
2013
, “
Development of a Lead-Zirconate-Titanate (PZT) Thin-Film Microactuator Probe for Intracochlear Applications
,”
Sens. Actuators A: Phys.
,
201
, pp.
1
9
.
33.
Luo
,
C.
,
Cao
,
G.
, and
Shen
,
I.
,
2012
, “
Enhancing Displacement of Lead-Zirconate-Titanate (PZT) Thin-Film Membrane Microactuators Via a Dual Electrode Design
,”
Sens. Actuators A: Phys.
,
173
(
1
), pp.
190
196
.
34.
Tzou
,
H.
, and
Bao
,
Y.
,
1997
, “
Nonlinear Piezothermoelasticity and Multi-Field Actuations—Part 1: Nonlinear Anisotropic Piezothermoelastic Shell Laminates
,”
ASME J. Vib. Acoust.
,
119
(
3
), pp.
374
381
.
35.
Varelis
,
D.
, and
Saravanos
,
D. A.
,
2004
, “
Coupled Finite Element for Non-Linear Laminated Piezoelectric Composite Shells With Applications to Buckling and Postbuckling Behavior
,”
AIAA
Paper No. AIAA 2004-1715.
36.
Kundu
,
C.
,
Maiti
,
D.
, and
Sinha
,
P.
,
2007
, “
Post Buckling Analysis of Smart Laminated Doubly Curved Shells
,”
Compos. Struct.
,
81
(
3
), pp.
314
322
.
37.
Lentzen
,
S.
,
Kłosowski
,
P.
, and
Schmidt
,
R.
,
2007
, “
Geometrically Nonlinear Finite Element Simulation of Smart Piezolaminated Plates and Shells
,”
Smart Mater. Struct.
,
16
(
6
), p.
2265
.
38.
Nanda
,
N.
,
2010
, “
Non-Linear Free and Forced Vibrations of Piezoelectric Laminated Shells in Thermal Environments
,”
IES J. Part A: Civ. Struct. Eng.
,
3
(
3
), pp.
147
160
.
39.
Pradyumna
,
S.
, and
Gupta
,
A.
,
2011
, “
Nonlinear Dynamic Stability of Laminated Composite Shells Integrated With Piezoelectric Layers in Thermal Environment
,”
Acta Mech.
,
218
(
3–4
), pp.
295
308
.
40.
Boroujerdy
,
M. S.
, and
Eslami
,
M.
,
2013
, “
Nonlinear Axisymmetric Thermomechanical Response of Piezo-Fgm Shallow Spherical Shells
,”
Arch. Appl. Mech.
,
83
(
12
), pp.
1681
1693
.
41.
Boroujerdy
,
M. S.
, and
Eslami
,
M.
,
2014
, “
Axisymmetric Snap-Through Behavior of Piezo-Fgm Shallow Clamped Spherical Shells Under Thermo-Electro-Mechanical Loading
,”
Int. J. Pressure Vessels Piping
,
120–121
, pp. 19–26.
42.
Pasquali
,
M.
, and
Gaudenzi
,
P.
,
2015
, “
A Nonlinear Piezoelectric Shell Model: Theoretical and Numerical Considerations
,”
J. Intell. Mater. Syst. Struct.
,
27
(
6
), pp.
724
742
.
43.
Pasquali
,
M.
, and
Gaudenzi
,
P.
,
2015
, “
A Nonlinear Formulation of Piezoelectric Shells With Complete Electro-Mechanical Coupling
,”
Meccanica
,
50
(
10
), pp.
2471
2486
.
44.
Tzou
,
H.
, and
Zhou
,
Y.
,
1997
, “
Nonlinear Piezothermoelasticity and Multi-Field Actuations—Part 2: Control of Nonlinear Deflection, Buckling and Dynamics
,”
ASME J. Vib. Acoust.
,
119
(
3
), pp.
382
389
.
45.
Zhou
,
Y.-H.
, and
Tzou
,
H.
,
2000
, “
Active Control of Nonlinear Piezoelectric Circular Shallow Spherical Shells
,”
Int. J. Solids Struct.
,
37
(
12
), pp.
1663
1677
.
46.
Tzou
,
H.
, and
Wang
,
D.
,
2003
, “
Micro-Sensing Characteristics and Modal Voltages of Linear/Non-Linear Toroidal Shells
,”
J. Sound Vib.
,
254
(
2
), pp.
203
218
.
47.
Giannopoulos
,
G.
,
Monreal
,
J.
, and
Vantomme
,
J.
,
2007
, “
Snap-Through Buckling Behavior of Piezoelectric Bimorph Beams—I: Analytical and Numerical Modeling
,”
Smart Mater. Struct.
,
16
(
4
), p.
1148
.
48.
Lee
,
C.-C.
,
Cao
,
G.
, and
Shen
,
I.
,
2010
, “
Effects of Residual Stresses on Lead-Zirconate-Titanate (PZT) Thin-Film Membrane Microactuators
,”
Sens. Actuators, A: Phys.
,
159
(
1
), pp.
88
95
.
49.
Luo
,
C.
,
Tai
,
W.
,
Yang
,
C.-W.
,
Cao
,
G.
, and
Shen
,
I.
,
2016
, “
Effects of Added Mass on Lead-Zirconate-Titanate Thin-Film Microactuators in Aqueous Environments
,”
ASME J. Vib. Acoust.
,
138
(
6
), p.
061015
.
50.
Tzou
,
H.
,
1993
, “
Piezoelectric Shells: Distributed Sensing and Control of Continua
,”
Solid Mechanics and Its Applications
,
Kluwer Academic, Dordrecht, The Netherlands
.
51.
Shen
,
I.
,
1997
, “
Active Constrained Layer Damping Treatments for Shell Structures: A Deep-Shell Theory, Some Intuitive Results, and an Energy Analysis
,”
Smart Mater. Struct.
,
6
(
1
), pp.
89
101
.
52.
Nath
,
Y.
,
Mahrenholtz
,
O.
, and
Varma
,
K.
,
1987
, “
Non-Linear Dynamic Response of a Doubly Curved Shallow Shell on an Elastic Foundation
,”
J. Sound Vib.
,
112
(
1
), pp.
53
61
.
53.
Reddy
,
J. N.
,
2006
,
Theory and Analysis of Elastic Plates and Shells
,
CRC Press
,
Boca Raton, FL
.
54.
Amabili
,
M.
,
2005
, “
Non-Linear Vibrations of Doubly Curved Shallow Shells
,”
Int. J. Non-Linear Mech.
,
40
(
5
), pp.
683
710
.
55.
Dickinson
,
S.
, and
Di Blasio
,
A.
,
1986
, “
On the Use of Orthogonal Polynomials in the Rayleigh-Ritz Method for the Study of the Flexural Vibration and Buckling of Isotropic and Orthotropic Rectangular Plates
,”
J. Sound Vib.
,
108
(
1
), pp.
51
62
.
56.
Gibbs
,
G. P.
, and
Fuller
,
C. R.
,
1992
, “
Excitation of Thin Beams Using Asymmetric Piezoelectric Actuators
,”
J. Acoust. Soc. Am.
,
92
(
6
), pp.
3221
3227
.
57.
Cottone
,
F.
,
Gammaitoni
,
L.
,
Vocca
,
H.
,
Ferrari
,
M.
, and
Ferrari
,
V.
,
2012
, “
Piezoelectric Buckled Beams for Random Vibration Energy Harvesting
,”
Smart Mater. Struct.
,
21
(
3
), p.
035021
.