Effects of van der Waals Force and Thermal Stresses on Pull-in Instability of Clamped Rectangular Microplates

Abstract

We study the influence of von Kármán nonlinearity, van der Waals force, and thermal stresses on pull-in instability and small vibrations of electrostatically actuated microplates. We use the Galerkin method to develop a tractable reduced-order model for electrostatically actuated clamped rectangular microplates in the presence of van der Waals forces and thermal stresses. More specifically, we reduce the governing two-dimensional nonlinear transient boundary-value problem to a single nonlinear ordinary differential equation. For the static problem, the pull-in voltage and the pull-in displacement are determined by solving a pair of nonlinear algebraic equations. The fundamental vibration frequency corresponding to a deflected configuration of the microplate is determined by solving a linear algebraic equation. The proposed reduced-order model allows for accurately estimating the combined effects of van der Waals force and thermal stresses on the pull-in voltage and the pull-in deflection profile with an extremely limited computational effort.

1 Introduction

Microelectromechanical sensors are currently used in automotive electronics, medical equipment, smart portable electronics, hard disk drives, and computer peripherals [1]. Electrostatically actuated Mi-croElectroMechanical Systems (MEMS) are a special class of MEMS that use Coulomb's force as the driving mechanism in sensing and actuation. Electrostatically actuated MEMS are used as accelerome-ters [2], pressure sensors [3], and mass flow and chemical sensors [4]. Electrostatically actuated MEMS also find applications in the characterization of the material properties of thin films, such as elastic constants and residual stresses, see for example [5]. In addition to sensing applications, electrostatically actuated MEMS are employed in high-frequency signal filtering [6] and relay switches [7]. A comprehensive review of modeling of electrostatically actuated MEMS can be found in [8].

An electrostatically actuated MEMS is generally comprised of a conductive deformable electrode suspended above a rigid grounded electrode [9]. A direct current (DC) voltage applied across the two bodies results in the deflection of the deformable electrode and the consequent change in the system capacitance. When an alternating current (AC) is superimposed on the DC voltage, harmonic motions of the system are obtained that can be used to realize resonant devices. The applied DC voltage has an upper limit beyond which the electrostatic force is not balanced by the elastic restoring force in the deformable conductor. Beyond this critical voltage, the deformable conductor snaps and touches the lower rigid plate. This phenomenon is known as the pull-in instability, and it was first experimentally observed in the late sixties [10, 11]. The critical displacement and the critical voltage associated with this instability are called the pull-in displacement and the pull-in voltage, respectively. Their accurate evaluation is crucial in the design of electrostatically actuated MEMS. In micro-mirrors [12] and micro-resonators [13], the designer avoids this instability in order to achieve stable motions; while in switching applications [7], the designer exploits this effect to optimize the performance of the device.

For a wide class of electrostatic MEMS, the deformable electrode is a flat body whose thickness h is much smaller than its characteristic in-plane dimension ℓ [14]. Such electrodes can be regarded as two-dimensional (2-D) plate-like bodies. Since h/ℓ ≪ 1, an approximate distributed model can be employed, where the system kinematics is described only through the displacement of points on the movable electrode mid-surface, see for example [15]. Linear and nonlinear problems for microplates have been studied in [16–23]. When the bending stiffness of the deformable electrode is negligible as compared to its in-plane stretching and g₀/ℓ ≪ 1, where g₀ is the initial gap, the electrode can be regarded as a linear elastic membrane. The membrane approximation is valid for ℓ/h ≥ 400, see for example [24]. Linear micromembranes have been studied in [25–29]. As discussed in [25], the membrane approximation yields accurate results for many MEMS devices such as micro-pumps made of thin glassy polymers and grating light valves comprised of stretched thin ribbons.

With the decrease in electrostatic MEMS dimensions from the micro to the nanoscale additional nanoscale surface forces, such as the Casimir force and the van der Waals force [30–33], should be considered, see for example [34]. At small scales, the nanoscale surface forces may overcome elastic restoring actions in the device and lead to the plates' sticking during the fabrication process. van der Waals force and Casimir force can both be connected with the existence of zero-point vacuum oscillations of the electromagnetic field [31–33]. The microscopic approach to the modeling of both van der Waals and Casimir forces can be formulated in a unified way using Quantum Field Theory, see for example [30–33]. It is found that the Casimir force is generally effective at larger separation distances between the bodies than the van der Waals force. Whereas the Casimir force between semi-infinite parallel plates is inversely proportional to the fourth power of the gap, van der Waals force is inversely proportional to the third power of the gap. The dependence of these forces on the dielectric properties of the plates and the filling medium is studied in details in [31], Section 4.1.1. It is important to note that van der Waals and Casimir forces cannot in general be considered to simultaneously act in MEMS, since they describe the same physical phenomenon at two different length scales.

MEMS can be subjected to considerable temperature variations during sensing operation, such as in monitoring aircraft condition and distributed satellites communication, as well as during device packaging [35]. Modeling of the effect of temperature changes on pull-in instability of electrostatically actuated microelectromechanical switches is addressed in [36–38]. In these studies, it is observed that even a moderate temperature increase may cause premature failure of the device due to a significant decrease in the pull-in voltage. On the other hand, temperature drops may lead to considerable increase in the pull-in voltage even outside the operational range of the switch [39, 40].

We consider rectangular microplates undergoing large displacements under the combined effect of electrostatic and nanoscale forces. We use the Galerkin method to develop a tractable reduced-order model for electrostatically actuated microplates. The reduced-order model is derived by taking a family of linearly independent kinematically admissible functions as basis functions for the transverse displacement and by decomposing the in-plane displacement vector as the sum of displacements for irrotational and isochoric waves in a 2-D medium. Basis functions for the transverse and the in-plane displacements are related through the nonlinear equation governing the plate in-plane motion. The governing equations of the reduced-order model are derived from the equation governing the transverse motion of the microplate. The model is specialized to the case where a single basis function for the transverse displacement is used to yield manageable solutions. In the static analysis, pull-in parameters are found by solving a system of two nonlinear algebraic equations for the transverse displacement amplitude and the load parameter. The eigenvalue problem corresponding to linear vibrations of the system about its statically deflected position is solved for the fundamental frequency. We show that the fundamental frequency goes to zero as pull-in conditions are approached. The pull-in parameters found from the eigenvalue analysis agree well with those derived from the static analysis. We investigate effects of nanoscale forces and thermal stresses on pull-in parameters and small vibrations of electrostatically actuated microplates.

The rest of the paper is organized as follows. In section 2, we present the governing nonlinear equations of motion for a von Kármán microplate under the simultaneous effects of thermal loading, electrostatic force, and nanoscale forces. We present expressions of the distributed loads due to either the Casimir or the van der Waals forces. In section 3, we introduce a reduced-order model for the considered device that is capable of accurately predicting its dynamics. The derivation of the reduced-order model follows a procedure typically used for studying deformations of thin two-dimensional structures. That is, in-plane inertial effects are neglected, and the resulting equation is solved for in-plane displacements in terms of transverse deflections which are then substituted in the equation governing the evolution of transverse deflection. Once transverse deflections have been computed, in-plane displacements can be found. In section 3, we also briefly outline the technique used to solve equations for the reduced order model. In section 4, we present our results, that include the pull-in parameters and the fundamental frequencies for rectangular microplates. We specialize our results to the case where nanoscale effects can be subsumed into the van der Waals force and we investigate the effect of the van der Waals force on the pull-in instability and the lowest frequency of the predeformed plate. Conclusions are summarized in section 5.

2 Problem statement

We consider a rectangular plate-like body of longer side ℓ and shorter side φℓ, φ ∈]0,1], occupying the three-dimensional region Ω × (−h/2, h/2) as depicted in Figure 1. The mid-surface Ω is described by rectangular Cartesian coordinates (x₁, x₂) aligned with the longer and with the shorter sides, respectively. We assume that the initial gap g₀ between the two conductors and the thickness h of the deformable plate are much smaller than the characteristic length ℓ. In the proposed model, we assume that g₀ and h can be of the same order of magnitude. Therefore, the maximum displacement that the device can undergo is of the order of the plate thickness h, but it is much smaller than the characteristic length ℓ, since h/ℓ ≪ 1. This implies that strains in the deformable electrode are small. Under these assumptions, we use the von Kármán plate theory to account for large deflections and small strains, see for example [41]. Neglecting the effect of the rotatory inertia, the von Kármán plate equations in Cartesian coordinates are [41]

$$\varrho h\ddot{\omega }+D\frac{{\partial }^2}{\partial x_k\partial x_k}\left(\frac{{\partial }^2\omega }{\partial x_j\partial x_j}\right)-h\frac{\partial }{\partial x_j}\left({\sigma }_{jk}\frac{\partial \omega }{\partial x_k}\right)-F_e-F_s=0,$$

(1a)

$$\begin{array}{cc}\varrho {\ddot{u}}_i-\frac{\partial {\sigma }_{ij}}{\partial x_j}=0,& i=1,2.\end{array}$$

(1b)

We use the Einstein summation convention, meaning that when an index variable appears twice in a single term we are summing over the range {1,2}; free index also span the range {1,2}. Furthermore, D = Eh³/(12(1 − ν²)) is the bending stiffness of the plate, E, ν, and ϱ are Young's modulus, Poisson's ratio and mass density of the plate material, that is assumed to be homogeneous and isotropic, u_i and w are the in-plane and the out-of-plane displacements of a point on the mid-surface, F_e is the Coulomb force, and F_s is the nanoscale surface force. From equation (la), we note that Coulomb and nanoscale forces act along the normal to the microplate mid-surface. This is due to the assumption of perfectly conducting material. Expressions for the Coulomb and the nanoscale surface forces are discussed below. We note that in the linear theory of a plate the transverse displacement and in-plane displacements are decoupled. Therefore, in the linear theory, the Coulomb force and nanoscale forces do not induce in-plane deformations of the microplate.

In the von Kármán plate theory, the components ε_ij of the in-plane strain tensor are given by

$${\varepsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}+\frac{\partial \omega }{\partial x_i}\frac{\partial \omega }{\partial x_j}\right).$$

(2)

We assume that the strain tensor admits the additive decomposition ${\varepsilon }_{ij}={\varepsilon }_{ij}^{\text{e}}+{\varepsilon }_{ij}^{\mathrm{\Theta }}$, where ${\varepsilon }_{ij}^{\text{e}}$ is the elastic strain and ${\varepsilon }_{ij}^{\mathrm{\Theta }}$ the thermal one. We set ${\varepsilon }_{ij}^{\mathrm{\Theta }}=\chi \mathrm{\Theta }{\delta }_{ij}$, where Θ is the uniform temperature rise, χ is the coefficient of thermal expansion assumed to be constant, and

$${\delta }_{ij}=\{\begin{array}{ccc}1& \text{if}& i=j\\ 0& \text{if}& i\ne j\end{array},$$

(3)

is the Kronecker delta. Assuming the response of the material to be linear elastic, the constitutive relation under the Kirchhoff assumption is, see for example [42],

$${\sigma }_{ij}=\frac{E}{1+v}\left({\varepsilon }_{ij}^e+\frac{v}{1-v}{\varepsilon }_{kk}^e{\delta }_{ij}\right)=\frac{E}{1+v}\left({\varepsilon }_{ij}+\frac{v}{1-v}{\varepsilon }_{kk}{\delta }_{ij}\right)-\frac{E\chi \mathrm{\Theta }}{1-v}{\delta }_{ij}.$$

(4)

Substituting for σ_ij from Equation (4) into equations (1) and using Equation (2), we obtain the following partial differential equations for the displacement fields u_t and w

$$\begin{array}{l}\varrho h\ddot{\omega }+D\frac{{\partial }^2}{\partial x_k\partial x_k}\left(\frac{{\partial }^2\omega }{\partial x_j\partial x_j}\right)-B\frac{\partial }{\partial x_l}\left(\left(\frac{1}{2}(\frac{\partial u_l}{\partial x_j}+\frac{\partial u_j}{\partial x_l}\right)+{\delta }_{lj}\frac{\partial u_k}{\partial x_k}\right)\frac{\partial \omega }{\partial x_j})\\ -\frac{B}{2}\frac{\partial }{\partial x_k}\left(\frac{\partial \omega }{\partial x_i}\frac{\partial \omega }{\partial x_j}\frac{\partial \omega }{\partial x_k}\right)+\frac{E\chi \mathrm{\Theta }}{1-v}\frac{{\partial }^2\omega }{\partial x_j\partial x_j}-F_e-F_s=0,\end{array}$$

(5a)

$$\begin{array}{l}\varrho h{\ddot{u}}_i-\frac{B}{2}\left((1-v)\frac{{\partial }^2{u}_i}{\partial x_k\partial x_k}+(1+v)\frac{\partial }{\partial x_i}\left(\frac{\partial u_k}{\partial x_k}\right)\right)\\ =\frac{B}{2}\left((1-v)\frac{\partial }{\partial x_k}\left(\frac{\partial \omega }{\partial x_i}\frac{\partial \omega }{\partial x_k}\right)+v\frac{\partial }{\partial x_i}\left(\frac{\partial \omega }{\partial x_k}\frac{\partial \omega }{\partial x_k}\right)\right),i=1,2,\end{array}$$

(5b)

where B = Eh/(1 − v²).

From an electrical point of view, the system depicted in Figure 1 behaves as a variable gap capacitor. By assuming that g₀/ℓ ≪ 1 and by neglecting fringing fields, the magnitude F_e of the Coulomb or the magnitude of the electrostatic force acting on the deformable electrode along its normal is given by the parallel plate approximation [9]

$$F_e=-\frac{{ϵ}_0{V}^2}{2g_0^2{(1+\widehat{w})}^2},$$

(6)

where ŵ = w/g₀ is the non-dimensional transverse displacement, ϵ₀ is the dielectric constant in vacuum, and V is the applied direct current voltage. Therefore, within the parallel plate approximation the expression for the electrostatic force depends only on the gap g. In [43], perturbation theory is used to develop expressions for the electrostatic force with corrections to the parallel plate approximations. It is shown that the parallel plate approximation is adequate when g₀ ≪ ℓ, since the first correction multiplies the ratio $g_0^2/{\ell }^2$. Consistent with the parallel plate approximation for the electrostatic force, and the small deformations assumption in the mechanical model, the distributed load due to the Casimir force per unit surface area is given by the proximity force approximation [44, 45], which for perfect conductors yields:

$$F_c=-\frac{\hslash c{\pi }^2}{240g_0^4{(1+\widehat{w})}^4},$$

(7)

where ћ is the Plank's constant and c is the speed of light in vacuum. Corrections to Equation (7) for geometries with known and fixed departures from the parallel configurations are reported in [44, 45]. For gaps smaller than the retardation length, that is, for gaps smaller than the wavelength of the virtual transitions responsible for the quantum dipole fluctuations [31, 33], the Casimir force reduces to the van der Waals distributed force per unit surface area [46–48]

$$F_{\text{vdW}}=-\frac{A}{6\pi g_0^3{(1+\widehat{w})}^3}.$$

(8)

In Equation (8), A is the Hamaker constant, see for example [31], with values in the range [0.4,4] × 10⁻¹⁹J. *

van der Waals and Casimir forces between parallel layered metallic surfaces have been extensively studied in the literature, see for example [31, 33]. For gold-coated aluminum surfaces, it is found that van der Waals force, see Equation (8), is effective in the gap range 0.5 nm – 4 nm. For gaps in the range 4 nm − 1 μm, there is a transition between the force-distance dependence g⁻³ (van der Waals force) to the force distance dependence g⁻⁴ (Casimir force). For gaps larger than 1 μm, the interaction between the plates is described by the Casimir force, see Equation (7). Therefore, for perfect conductors, at large separation distances the interaction force is independent of the material properties of plates, whereas, as the gap decreases, the interaction force is affected by the material properties of the system. For dielectric bodies, the van der Waals force is effective at larger distances as shown in [49].

The effect of the Casimir force on pull-in parameters of von Kármán microplates has been studied in [20, 21]. In the present work, we focus only on the van der Waals force. Therefore, we set F_s = F_vdW in equation (5a).

We introduce the non-dimensional time t̂ = t/τ, where

$${\tau }^2=\frac{12\varrho {\ell }^4}{E{h}^2}(1-{\nu }^2),$$

(9)

and the non-dimensional in-plane displacement ${\widehat{u}}_i=u_i\ell /g_0^2$. Coordinates x₁ and x₂ are nondimension-alized as x̂₁ = x₁/ℓ, x̂₂ = x₂/ℓ. Henceforth, we use a superimposed dot to denote time derivative with respect to t̂. Also, we drop the superimposed hat on non-dimensional variables. Thus, equations (5) become

$$\begin{array}{l}\ddot{w}+\frac{{\partial }^2}{\partial x_k\partial x_k}\left(\frac{{\partial }^{2}w}{\partial x_j\partial x_j}\right)-12\alpha \frac{\partial }{\partial x_l}\left(\left(\frac{1}{2}\left(\frac{\partial u_l}{\partial x_j}+\frac{\partial u_j}{\partial x_l}\right)+{\delta }_{\mathit{\text{ij}}}\frac{\partial u_k}{\partial x_k}\right)\frac{\partial w}{\partial x_j}\right)\\ -6\alpha \frac{\partial }{\partial x_k}\left(\frac{\partial w}{\partial x_j}\frac{\partial w}{\partial x_j}\frac{\partial w}{\partial x_k}\right)+\beta \frac{{\partial }^{2}w}{\partial x_j\partial x_j}+\frac{\lambda }{{(1+w)}^2}+\frac{\mu }{{(1+w)}^3}=0,\end{array}$$

(10a)

$$\gamma {\ddot{u}}_i-(1-\nu )\frac{{\partial }^2{u}_i}{\partial x_k\partial x_k}-(1+\nu )\frac{\partial }{\partial x_i}\left(\frac{\partial u_k}{\partial x_k}\right)=(1-\nu )\frac{\partial }{\partial x_k}\left(\frac{\partial w}{\partial x_i}\frac{\partial w}{\partial x_k}\right)+\nu \frac{\partial }{\partial x_i}\left(\frac{\partial w}{\partial x_k}\frac{\partial w}{\partial x_k}\right),i=1,2,$$

(10b)

where

$$\begin{array}{l}\alpha =\frac{g_0^2}{h^2},\beta =12\frac{\chi \mathrm{\Theta }{\ell }^2}{h^2}(1+v),\gamma =\frac{h^2}{6{\ell }^2},\\ \lambda =\frac{6{ϵ}_0{V}^2{\ell }^4}{E{h}^3{g}_0^3}(1-v^2),\mu =\frac{2A{\ell }^4}{\pi E{h}^3{g}_0^4}(1-v^2).\end{array}$$

(11)

The non-dimensional parameters β, λ, and μ are indicators of the effects of thermal stress, the Coulomb force, and the van der Waals force, respectively, on the MEMS dynamics. We quantify the order of magnitude of different terms in equations (10) through the order of magnitude of the corresponding non-dimensional parameters. Since the von Kármán approximation holds for h/ℓ ≪ 1, we neglect the inertial term in equation (10b) to obtain

$$(1-v)\frac{{\partial }^2{u}_i}{\partial x_k\partial x_k}+(1+v)\frac{\partial }{\partial x_i}\left(\frac{\partial u_k}{\partial x_k}\right)=-(1-v)\frac{\partial }{\partial x_k}\left(\frac{\partial w}{\partial x_i}\frac{\partial w}{\partial x_k}\right)-v\frac{\partial }{\partial x_i}\left(\frac{\partial w}{\partial x_k}\frac{\partial w}{\partial x_k}\right),i=1,2.$$

(12)

We note that Equation (12) implies a time-independent link between the in-plane and the transverse displacements.

We consider the boundary Γ of Ω to be clamped. The kinematic boundary conditions for a clamped edge are [50]

$$\begin{array}{ccccc}w=0& \text{and}& \frac{\partial w}{\partial x_1}=0& \text{on}& x_2=\left\{0,\phi \right\}\end{array},\begin{array}{ccc}\frac{\partial w}{\partial x_2}=0& \text{on}& x_1=\left\{0,1\right\},\end{array}$$

(13a)

$$\begin{array}{cc}{u}_i=0,& i=1,2.\end{array}$$

(13b)

Initial conditions are not needed since we either study static deformations of the MEMS, or analyze frequencies of small vibrations around an electrostatically deformed configuration.

3 Reduced-order model

A closed-form solution of the initial-boundary-value problem defined in equation set (10) and boundary conditions in equation set (13) cannot be found. We construct an approximate solution by expressing the displacement fields u_i, with i = 1, 2, and w as

$$w(x_1,x_2,t)=\sum _{n=1}^N{\overline{w}}_{(n)}(x_1,x_2){\zeta }_{(n)}(t)={\text{W}}^{\text{T}}(x_1,x_2)\zeta (t),$$

(14a)

$$\begin{array}{cc}{u}_i(x_1,x_2,t)=\sum _{p=1}^P{\overline{u}}_{(p)i}(x_1,x_2){\xi }_{(p)}(t)={\text{U}}_i^{\text{T}}(x_1,x_2)\xi (t),& i=1,2,\end{array}$$

(14b)

where w̄_(n) and ū_(p)i are basis functions for the transverse and the in-plane displacements and ζ_(n) and ξ_(p) are the corresponding amplitude parameters or equivalently the mode participation factors. Basis functions are collected into the N- column vector W and into the P- column vector U_i, respectively, and amplitudes are collected into the N- column vector ζ and into the P- column vector ξ, respectively. Each basis function satisfies the corresponding kinematic boundary conditions prescribed through equation set (13).

In [20], it is shown that accurate solutions for pull-in parameters of von Kármán microplates can be obtained by using a properly selected single basis function for the transverse displacement in equation (14a), that is, by setting N = 1 in equation (14a). Therefore, we derive the reduced-order model by approximating the transverse displacement field as

$$w(x_1,x_2,t)=\overline{w}(x_1,x_2)\zeta (t).$$

(15)

3.1 Basis functions for in-plane displacement

Basis functions for the in-plane displacement can be determined by solving the following linear eigenvalue problem associated with equation (10b), see for example [41],

$$\begin{array}{cc}{k}^2{\overline{u}}_i+(1-v)\frac{{\partial }^2{\overline{u}}_i}{\partial x_k\partial x_k}+(1+v)\frac{\partial }{\partial x_i}\left(\frac{\partial {\overline{u}}_k}{\partial x_k}\right)=0,& i=1,2,\end{array}$$

(16)

where κ is the wave number. Following [41], we decompose the in-plane displacement as $u_i=u_i^{\text{n}}+u_i^t$, with i = l, 2, where $u_i^{\text{n}}$ and $u_i^t$ are displacements associated, respectively, with the longitudinal and the transverse waves, and satisfy

$$\frac{\partial {\overline{u}}_1^n}{\partial x_2}-\frac{\partial {\overline{u}}_2^n}{\partial x_1}=0,$$

(17a)

$$\frac{\partial {\overline{u}}_1^t}{\partial x_1}-\frac{\partial {\overline{u}}_2^t}{\partial x_2}=0.$$

(17b)

Therefore, Equation (16) is equivalent to the following pair of partial differential equations [41]

$$\begin{array}{cc}\frac{{\partial }^2{\overline{u}}_i^n}{\partial x_k\partial x_k}+{\eta }_n^2{\overline{u}}_i^n=0,& i=1,2,\end{array}$$

(18a)

$$\begin{array}{cc}\frac{{\partial }^2{\overline{u}}_i^t}{\partial x_k\partial x_k}+{\eta }_t^2{\overline{u}}_i^t=0,& i=1,2,\end{array}$$

(18b)

where ${\eta }_{\text{n}}=\kappa /\sqrt{2}$ and ${\eta }_{\text{t}}=\kappa /\sqrt{1-\nu }$ are wave numbers of the longitudinal and the transverse waves, respectively.

By integrating over the domain Ω and by applying Green's formulas to transform surface integrals into line integrals, equation set (17) implies the following set of boundary conditions for the displacements associated with the longitudinal and the transverse waves

$$\begin{array}{cccc}{u}_1^{\text{n}}(x_1,0)=0,& u_2^{\text{n}}(1,x_2)=0,& u_1^{\text{n}}(x_1,\phi )=0,& u_2^{\text{n}}(0,x_2)=0,\end{array}$$

(19a)

$$\begin{array}{cccc}{u}_2^{\text{t}}(x_1,0)=0,& u_1^{\text{t}}(1,x_2)=0,& u_2^{\text{t}}(x_1,\phi )=0,& u_1^{\text{t}}(0,x_2)=0.\end{array}$$

(19b)

Additional boundary conditions are provided by equation (13b) through

$$\begin{array}{cccc}{u}_2^{\text{n}}(x_1,0)=0,& u_1^{\text{n}}(1,x_2)=0,& u_2^{\text{n}}(x_1,\phi )=0,& u_1^{\text{n}}(0,x_2)=0,\end{array}$$

(20a)

$$\begin{array}{cccc}{u}_1^{\text{t}}(x_1,0)=0,& u_2^{\text{t}}(1,x_2)=0,& u_1^{\text{t}}(x_1,\phi )=0,& u_2^{\text{t}}(0,x_2)=0.\end{array}$$

(20b)

Therefore, basis functions for the in-plane displacement are given by

$$\begin{array}{ccc}{\overline{u}}_{(p)i}=A_{mn}sin(m\pi x_1)sin\left(\frac{n\pi x_2}{\phi }\right),& m=1,\cdots ,\overline{p},& n=1,\cdots ,\overline{q},\end{array}$$

(21)

where A_mn are constants used for normalization, and p = (m − 1)q̄ + n for m = 1, …, p̄ and n = 1, …, q̄.

3.2 Relation between ξ, and ζ

In order to express the coefficients ξ, in terms of the single coefficient ζ, we substitute expressions for the displacement fields given in equations (14b) and (15) into Equation (12), premultiply both sides of the resulting equation with the transpose of the in-plane basis functions vector U, and integrate over the domain Ω. By applying the divergence theorem and by imposing the boundary conditions in equation (13b) on the boundary integrals, we finally obtain

$${\xi }_{(p)}={\mathrm{\Xi }}_{(p)}{\zeta }^2,$$

(22)

where the scalar quantity Ξ_(p) is defined through

$$\begin{array}{ll}{\mathrm{\Xi }}_{(p)}& =-{\mathit{\int }}_{\Omega }\left((1-v)\frac{\partial {\overline{u}}_{(p)j}}{\partial x_k}\frac{\partial \overline{w}}{\partial x_k}\frac{\partial \overline{w}}{\partial x_j}+v\left(\frac{\partial {\overline{u}}_{(p)k}}{\partial x_k}\right)\frac{\partial \overline{w}}{\partial x_j}\frac{\partial \overline{w}}{\partial x_j}\right)\text{d}\Omega \\ & \times {\left({\mathit{\int }}_{\Omega }\left((1-v)\frac{\partial {\overline{u}}_{(p)j}}{\partial x_k}\frac{\partial {\overline{u}}_{(p)j}}{\partial x_k}+(1+v){\left(\frac{\partial {\overline{u}}_{(p)k}}{\partial x_k}\right)}^2\right)\text{d}\Omega \right)}^{-1}.\end{array}$$

(23)

Equation (22) implies that the full set of amplitude parameters describing the in-plane motion can be obtained from the single transverse amplitude parameter ζ through the scalar parameter Ξ_(p) defined in Equation (23). In addition, we note that the relation between the in-plane motion and the transverse motion is nonlinear and that the in-plane motion vanishes in the linear theory. Because of the P internal constraints (22), the reduced-order model for the microplate has 1 degree of freedom.

3.3 Governing equations of the reduced-order model

The reduced-order model is obtained by premultiplying both sides of equation (10a) with w̄, integrating the resulting equation over Ω, and substituting into the resulting equation the approximations given by equations (14b), (15), and (22). We thus obtain

$$\begin{array}{l}{\mathrm{\int }}_{\Omega }{\overline{w}}^2\ddot{\zeta }\text{d}\Omega +{\mathrm{\int }}_{\Omega }\overline{w}\frac{{\partial }^2}{\partial x_k\partial x_k}\left(\frac{{\partial }^2\overline{w}}{\partial x_j\partial x_j}\right)\zeta \text{d}\Omega \\ -12\alpha {\mathrm{\int }}_{\Omega }\overline{w}\frac{\partial }{\partial x_j}\left({\overline{\varepsilon }}_{jk}\frac{\partial \overline{w}}{\partial x_k}+\frac{1}{2}\frac{\partial \overline{w}}{\partial x_k}\frac{\partial \overline{w}}{\partial x_k}\frac{\partial \overline{w}}{\partial x_j}\right){\zeta }^3\text{d}\Omega \\ +\beta {\mathrm{\int }}_{\Omega }\overline{w}\frac{{\partial }^2\overline{w}}{\partial x_k\partial x_k}\zeta \text{d}\Omega +\lambda {\mathrm{\int }}_{\Omega }\frac{\overline{w}}{{(1+\overline{w}\zeta )}^2}\text{d}\Omega +\mu {\mathrm{\int }}_{\Omega }\frac{\overline{w}}{{(1+\overline{w}\zeta )}^3}\text{d}\Omega =0,\end{array}$$

(24)

Where

$${\overline{\varepsilon }}_{ij}=\sum _{p=1}^P{\mathrm{\Xi }}_{(p)}\left(\frac{1-v}{2}\left(\frac{\partial {\overline{u}}_{(p)i}}{\partial x_j}+\frac{\partial {\overline{u}}_{(p)j}}{\partial x_i}\right)+v\left(\frac{\partial {\overline{u}}_{(p)k}}{\partial x_k}\right){\delta }_{ij}\right).$$

(25)

From equations (23) and (25), we note that the normalization constants A_mn introduced in Equation (21) do not affect ε̄_ij. By using the divergence theorem and by imposing the boundary conditions in equation (13a), we obtain the following equation for the reduced-order system

$$m\ddot{\zeta }+(k_1+\beta {\kappa }_2)\zeta +\alpha k_3{\zeta }^3+\lambda f_e(\zeta )+\mu f_{\text{vdW}}(\zeta )=0,$$

(26)

where

$$m$$

(27a)

$$\begin{array}{cc}{k}_1={{\mathrm{\int }}_{\Omega }\left(\frac{{\partial }^2\overline{w}}{\partial x_k\partial x_k}\right)}^2\text{d}\Omega ,& k_2=-{{\mathrm{\int }}_{\Omega }\left(\frac{\partial \overline{w}}{\partial x_k}\frac{\partial \overline{w}}{\partial x_k}\right)}^2\text{d}\Omega ,\end{array}$$

(27b)

$$k_3=12{\mathrm{\int }}_{\Omega }({\overline{\varepsilon }}_{jk}\frac{\partial \overline{w}}{\partial x_j}\frac{\partial \overline{w}}{\partial x_k}+\frac{1}{2}{\left(\frac{\partial \overline{w}}{\partial x_k}\frac{\partial \overline{w}}{\partial x_k}\right)}^2)\text{d}\Omega ,$$

(27c)

$$\begin{array}{cc}{f}_e(\zeta )={\mathrm{\int }}_{\Omega }\frac{\overline{w}}{{(1+\overline{w}\zeta )}^2}\text{d}\Omega ,& f_{\text{VdW}}(\zeta )={\mathrm{\int }}_{\Omega }\frac{\overline{w}}{{(1+\overline{w}\zeta )}^3}\text{d}\Omega .\end{array}$$

(27d)

Equation (26) describes a nonlinear mass-spring single degree-of-freedom system, where the spring force is the combination of a number of concurring phenomena: linear bending (k₁ζ), thermal stress (k₂ζ), nonlinear membrane stretching (k₃ζ³), the Coulomb force (λf_e(ζ)), and the van der Waals force (μf_vdW(ζ)). From Equation (26), we note that an increase in the device temperature, that is, decreasing the parameter β, results in mechanical softening effect, whereas a decrease in the device temperature yields a mechanical stiffening effect. As the parameter β decreases mechanical buckling due to thermal stress can be reached.

In what follows, we discuss two equivalent methods to extract static pull-in parameters. The first method is based on the solution of the static problem, while the second one is based on the study of small vibrations of the system around its static equilibrium configurations. The method based on the vibration analysis can be used to experimentally determine the pull-in voltage without potentially damaging the MEMS [51].

3.4 Extraction of pull-in parameters from the static problem

At the onset of instability the system's tangent stiffness

$$K(\beta ,\lambda ,\mu ,\zeta )=k_1+\beta k_2+3\alpha k_3{\zeta }^2+{\lambda \text{g}}_e(\zeta )+\mu g_{\text{vdW}}(\zeta )$$

(28)

vanishes. The quantities in Equation (28) are defined by

$$g_e(\zeta )=\frac{\text{d}{f}_e}{\text{d}\zeta }(\zeta )=-{\mathrm{\int }}_{\Omega }\frac{2{\overline{w}}^2}{{(1+\overline{w}\zeta )}^3}\text{d}\Omega ,$$

(29a)

$$g_{\text{vdW}}(\zeta )=\frac{d{f}_{\text{vdW}}}{\text{d}\zeta }(\zeta )=-{\mathrm{\int }}_{\Omega }\frac{3{\overline{w}}^2}{{(1+\overline{w}\zeta )}^4}\text{d}\Omega .$$

(29b)

Therefore, at the pull-in instability, the reduced-order model satisfies Equation (26) with ζ̈ = 0 and the additional condition

$$K(\beta ,\lambda ,\mu ,\zeta )=0.$$

(30)

We solve the problem for λ = 0 to compute the critical value, say μ_cr, of the van der Waals force parameter. When μ = μ_cr the system collapses spontaneously without applying any voltage, that is, for V = 0. We note that the critical value of the van der Waals force parameter depends on the parameter β that directly affects the mechanical stiffness. By eliminating μ from equations (26) and (30) with ζ̈ = 0 and λ = 0, we obtain the nonlinear equation

$$((k_1+\beta k_2){\zeta }_{\text{cr}}+\alpha k_3{\zeta }_{\text{cr}}^3){\text{g}}_{\text{vdW}}({\zeta }_{\text{cr}})=((k_1+\beta k_2+3\alpha k_3{\zeta }_{\text{cr}}^2)f_{\text{vdW}}({\zeta }_{\text{cr}}),$$

(31)

whose positive first root gives ζ_cr. The corresponding value of μ_cr is given by

$${\mu }_{\text{cr}}=-\frac{(k_1+\beta k_2){\zeta }_{\text{cr}}+\alpha k_3{\zeta }_{\text{cr}}^3}{f_{\text{vdW}}({\zeta }_{\text{cr}})}=-\frac{k_1+\beta k_2+3\alpha k_3{\zeta }_{\text{cr}}^2}{g_{\text{vdW}}({\zeta }_{\text{cr}})}.$$

(32)

The effect of the van der Waals force on pull-in parameters λ_pi and ║w_PI║_∞ is investigated by solving equations (26) and (30) with ζ̈ = 0 and variable λ for different values of μ in the range [0, μ_cr]. Here, ║ · ║_∞ is defined as max_(x1,x2)∈Ω | · |. By solving Equation (30) for λ and by substituting in Equation (26), we obtain the following nonlinear equation in ζ_PI:

$$((k_1+\beta k_2){\zeta }_{\text{PI}}+\alpha k_3{\zeta }_{\text{PI}}^3+\mu f_{\text{vdW}}({\zeta }_{\text{PI}}))g_{\text{vdW}}({\zeta }_{\text{PI}})=(k_1+\beta k_2+3\alpha k_3{\zeta }_{\text{PI}}^2+\mu g_{\text{vdW}}({\zeta }_{\text{PI}}))f_{\text{vdW}}({\zeta }_{\text{PI}}).$$

(33)

The lowest positive root of Equation (33) gives the non-dimensional pull-in deflection ζ_PI. The corresponding non-dimensional pull-in voltage is thereby determined from Equation (26) with ζ̈ = 0 as

$${\lambda }_{\text{PI}}=-\frac{(k_1+\beta k_2){\zeta }_{\text{PI}}+\alpha k_3{\zeta }_{\text{PI}}^3+\mu f_{\text{vdW}}({\zeta }_{\text{PI}})}{f_{\text{e}}({\zeta }_{\text{PI}})}=-\frac{k_1+\beta k_2+3\alpha k_3{\zeta }_{\text{PI}}^2+\mu g_{\text{vdW}}({\zeta }_{\text{PI}})}{g_{\text{e}}({\zeta }_{\text{PI}})}.$$

(34)

3.5 Buckling thermal stress parameter from the static problem

The buckling thermal stress parameter β_B for a given λ and a given μ is determined by eliminating β from equations (26) and (28) with ζ̈ = 0. The lowest positive root of the resulting nonlinear equation

$$\lambda (f_{\text{e}}\left({\zeta }_{\text{B}}\right)-{\zeta }_{\text{B}g\text{e}}({\zeta }_{\text{B}}))+\mu (f_{\text{vdW}}({\zeta }_{\text{B}}))-{\zeta }_{\text{BgvdW}}({\zeta }_{\text{B}}))-2\alpha k_3{\zeta }_{\text{B}}^3=0,$$

(35)

gives the value ζ _B of the displacement parameter corresponding to buckling instability. The parameter β_B is obtained as

$${\beta }_{\text{B}}=-\frac{k_1+3\alpha k_3{\zeta }_{\text{B}}^2+{\lambda }_{g\text{e}}({\zeta }_{\text{B}})+\mu g_{\text{vdW}}({\zeta }_{\text{B}})}{k_2}.$$

(36)

3.6 Extraction of pull-in parameters from the linear eigenvalue problem

We determine the lowest frequency of the deflected plate at a given solution (ζ, λ, μ) through the following procedure. We perturb the equilibrium state ζ with a harmonic term exp (iωt) as ζ(t) + ζ̲ exp (iωt), where ζ̲ is a constant such that |ζ̲| ≪|ζ| and $\text{i}=\sqrt{-1}$. Substituting into Equation (26) and retaining terms linear in ζ̲, we obtain the following relation for the fundamental frequency w₀

$$K(\beta ,\lambda ,\mu ,\zeta )-w_0^{2}m=0.$$

(37)

Since the tangent stiffness defined in Equation (28) is zero at pull-in, it follows that at pull-in the lowest natural frequency of the system equals zero. This can be viewed as an alternative way of defining the static pull-in [18,19,51,52].

4 Results and discussion

Integrals appearing in the governing equations of the reduced-order model, including equations (23) and (27), are evaluated using the 21 x 21 Gauss quadrature rule. We established convergence of the adopted Gauss quadrature rule by checking that the values of the above integrals do not change by considering more quadrature points for numerical integration. Results presented below are for ν = 0.25.

For the transverse displacement, we use the following kinematically admissible function in Equation (15)

$$\begin{array}{ll}\overline{w}(x_1,x_2)& =A_0(\frac{cosh({\vartheta }_0{x}_1)-cosh({\vartheta }_0{x}_1)}{cosh{\vartheta }_0-cosh{\vartheta }_0}-\left(\frac{sinh({\vartheta }_0{x}_1)-sinh({\vartheta }_0{x}_1)}{sinh{\vartheta }_0-sinh{\vartheta }_0}\right)\\ & \times \left(\frac{cosh\left(\frac{{\vartheta }_0{x}_2}{\phi }\right)-cosh\left(\frac{{\vartheta }_0{x}_2}{\phi }\right)}{cosh{\vartheta }_0-cosh{\vartheta }_0}-\frac{sinh\left(\frac{{\vartheta }_0{x}_2}{\phi }\right)-sinh\left(\frac{{\vartheta }_0{x}_2}{\phi }\right)}{sinh{\vartheta }_0-sinh{\vartheta }_0}\right),\end{array}$$

(38)

where ϑ₀ = 4.73004 is the lowest nonzero root of the transcendental equation cosh ϑ cos ϑ = 1, and the constant A₀ is chosen by normalizing w̅ with respect to its maximum value. The basis function in Equation (38) equals the product of the first transverse modes of vibration of clamped-clamped Euler beams of lengths 1 and φ, respectively.

Since w̅ in Equation (38) is symmetric about the axes x₁= 1/2 and x₂ = φ/2, the basis functions u_(p)i in Equation (21) with mod (mod(m, 2) + mod(n, 2), 2) = 0 imply Ξ_(p) = 0 in Equation (23). Here, mod(m, n) gives the remainder in the division of m by n. Therefore, the number of degrees of freedom in equation (14b) is given by

$$P=\frac{1}{2}(\overline{p}\overline{q}-mod(\overline{p}\overline{q},2)),$$

(39)

where only basis functions u_(p)i such that mod (mod(m, 2) + mod(n, 2), 2) = 1 are considered. Results below are computed by selecting p̅ = q̅.

In [20], it is shown that convergent solutions for pull-in parameters of rectangular von Kármán mi-croplates are obtained for p̅ = q̅ = 4, corresponding to P = 8. Thus, results presented herein are computed with 9 basis functions in the reduced-order model.

4.1 Critical value of van der Waals force parameter

Figure 2 shows the variation of the critical van der Waals force parameter μ_cr with respect to the aspect ratio φ for two values of $\alpha =g_0^2/h^2$ and for λ = 0 and β = 0. Numerical data in Figure 2 are interpolated with

$${\mu }_{\text{cr}}=f(\alpha )\left(1+\sum _{k=2}^4{\phi }^{-k}\right),$$

(40)

where the function f(α), plotted in Figure 3, is given by

$$f(\alpha )=3.91\alpha +36.7.$$

(41)

We note that the critical parameter μ_cr increases rapidly as the plate aspect ratio φ decreases. We further remark that if two opposite edges are clamped and the remaining two edges are left traction free, the MEMS can be modeled as a plate undergoing cylindrical deformations [52]. In the case of slender electrodes, effect of fringing fields that has been discarded in this work becomes relevant as discussed in [53-55]. It is highly unlikely that clamped rectangular plates with φ ≪ 1 will be used in practice due to their relatively high stiffness. For φ ≪ 1 and the longer edges traction free, the MEMS can be modeled as a beam for which the analysis in [53, 54] applies.

For the Casimir force, the function f(α) in Equation (41) equals 2.31α + 35.8; see [20]. As pointed out in [20], Equation (40) gives approximate values of μ_cr

4.2 Effect of van der Waals force on pull-in parameters

In Figures 4 and 5, we plot the pull-in parameters versus μ in the range [0, μ_cr] for α = 1 and β = 0 for a square and a rectangular plate. As μ increases the pull-in parameter λ_PI decreases monotonically from its maximum value ${\lambda }_{\text{PI}}^{max}$ corresponding to μ = 0 to its minimum value 0 corresponding to μ = μ_cr.

The value μ = μ_cr. represents the intersection of the curves with the horizontal axis. With an increase in μ, the non-dimensional maximum transverse displacement decreases monotonically from its maximum value $\parallel w_{\text{PI}}{\parallel }_{\infty }^{max}$ attained for μ = 0 to its minimum value corresponding to μ = μ_cr. For a MEMS made of a specified material, Equation (11)₅ implies that ${\ell }^4/(h^3{g}_0^4)$ is proportional to μ_cr. Thus, μ_cr provides a quantitative indication of the devices' size that can be safely fabricated. This means that reduced deflection ranges are allowable for devices having a large value of μ_cr since μ_cr is inversely proportional to $g_0^4$. By comparing results depicted in Figures 4 and 5, we conclude that a change in the aspect ratio of a plate from 1 to 1/2 significantly increases ${\lambda }_{\text{PI}}^{max}$ and it does not noticeably affect the difference $\parallel {\overline{w}}_{\text{PI}}{\parallel }_{\infty }^{max}-\parallel w_{\text{cr}}{\parallel }_{\infty }.$

For μ = 0.3μ_cr and α = 4, we depict in Figure 6 deformed shapes of a square plate when λ ⋍ 200 and ‖w‖_∞ ⋍ 0.59, and of a rectangular plate with φ = 1/2 when λ ⋍ 1500 and ‖w‖_∞ ⋍ 0.59. In both cases thermal effects are discarded by selecting β = 0. Fringe plots of the van der Waals pressure (cf. Equation (8)) are also exhibited in the same plots. For the same value of ‖w‖_∞ the maximum magnitude of the van der Waals pressure for the rectangular plate is nearly an order of magnitude higher as compared to that for the square plate; note that F_vdW in Equation (8) does not depend upon the applied voltage. Its value is negative because it acts along the negative z-axis.

4.3 Effect of thermal stress on pull-in and critical van der Waals force parameters

In Figure 7, we plot the variation with β of the critical van der Waals force parameter μ_cr for a square plate and a rectangular plate with φ = 1/2, and also two values of α. We recall that the parameter β in Equation (11) is a measure of the prestress due to a uniform change Θ in the microplate temperature. For each case examined, we observe that μ_cr increases with decreasing temperature and decreases with increasing temperature. Indeed, from equations (11) and (28), a positive increment of temperature decreases the overall stiffness, whereas a negative increment increases the overall stiffness. Thus, the temperature affects the minimum size of the device that can be safely fabricated.

In the absence of van der Waals force, that is, for μ = 0, Figure 8 exhibits the variation with the prestress parameter β of the pull-in voltage λ_PI for a square plate, and a rectangular plate with p = 1/2, and also two values of α. We note that magnitudes of the compressive and tensile thermal stresses are limited, respectively, by the buckling instability of the MEMS and the tensile strength of the material of the MEMS.

4.4 Pull-in parameters from the analysis of frequencies of a deformed plate

In Figure 9, we report, for two different values of μ, the fundamental frequency w₀ of the deflected microplate versus A for α = 4 and β = 0 and for a square and a rectangular plate. The natural frequency is normalized with respect to the value w̅₀ corresponding to λ = 0; values of non-dimensional w̅₀ for Figures 9 and 10 are reported in

Table 1. Values of w̅₀ in Figures 9, 10, 11 and 12.

α	φ = 1		φ = 1/2
	μ = 0	μ ⋍ 0.3μcr	μ = 0	μ ⋍ 0.3μcr
1	36.1	34.2	98.6	93.4
4	36.1	34.2	98.6	93.5

. The corresponding dimensional frequencies equal ${\overline{w}}_0\sqrt{D/\left(\varrho {\ell }^4\right)}$. As reported in [56], for μ = 0 and λ = 0, the values of the non-dimensional frequency w̅₀ equal 36.108 and 98.592 for the square and the rectangular plate with φ = 1/2, respectively.

The variation ofw₀/w̅₀ with A is non-mono tonic due to the combined effect of the strain hardening represented by k₃ζ³ and the softening effect introduced by the Coulomb and the van der Waals forces. Indeed, from Equation (26), it is clear that the overall behavior of the Coulomb and the van der Waals forces is equivalent to a nonlinear spring with negative spring constant. When the strain-hardening effect is negligible the fundamental frequency monotonically decreases to zero, as it is typically predicted by the linear plate theory.

Results reported in Figure 10 show that for α = 1 and β = 0 the fundamental frequency monoton-ically decreases, implying that in this case the softening effect related to the Coulomb and the van der Waals forces overwhelms the strain hardening effect due to the mechanical nonlinearity. We note that the value of A corresponding to w₀ = 0 agrees with the λ_PIfound from Figures 4 and 5 which of course should be true since in both cases K (ζ, λ, μ) = 0. This provides an alternative way to find the pull-in parameters.

4.5 Effect of thermal stress on the fundamental frequency of predeformed plates

In Figures 11 and 12, we plot for μ = 0 the non-dimensional fundamental frequency of the deformed plate versus β for two values of λ. The natural frequency is normalized with respect to the value w̅₀ corresponding to λ = 0, μ = 0 and β = 0 reported in Table 1. The value of β for which w₀ = 0 corresponds to the buckling instability of the MEMS plate, and is found with the method explained in section 3.5. For β ≤ β_B, the fundamental frequency is computed from Equation (37) with μ = 0, the constant value of λ, and the corresponding ζ satisfying Equation (26).

For α = 4 (g₀/h = 2) and λ ≠ 0 the trend is non-monotonic due to the concurrent effects of the strain hardening (αk₃ζ³) and the thermal softening (βk₂ζ). For the case studied, when α = 1 the thermal softening is dominant with respect to the strain hardening, and the fundamental frequency monotonically decreases. For λ = μ = 0, the buckling is due to equal biaxial thermal stresses of β ⋍ 53 for both α = 1 and 4. In the absence of the Coulomb and the van der Waals forces the thermal stress at the onset of buckling is unaffected by the initial gap between the two plates.

5 Conclusions

We have studied effects of thermal stress and nanoscale forces on pull-in instability and resonant behavior of electrostatically actuated microplates. Mechanical nonlinearities are accounted for by mod-eling the deformable microplate using the von Kármán plate theory. The thermal stress is modeled as a homogenous residual stress depending on the microplate temperature. Nanoscale surface forces are described using the van der Waals force. We have derived a simple and tractable mass spring single degree of freedom model for analyzing the behavior of the considered device. The reduced-order model is derived by using a single basis function for the transverse displacement and eight basis functions to describe strain hardening due to membrane stretching.

Results show that the pull-in voltage and the pull-in deflection are strongly affected by thermal stresses and the van der Waals force. More specifically, as the temperature of the MEMS increases the pull-in voltage decreases. Moreover, the van der Waals force becomes more relevant as the MEMS size is reduced and can potentially lead to the spontaneous collapse of the system in absence of applied voltage.