1190

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 22, NO. 5, OCTOBER 2013

Full Operational Range Dynamic Modeling of

Microcantilever Beams
Mansour Abtahi, Gholamreza Vossoughi, and Ali Meghdari

AbstractMicrocantilever beams are frequently utilized in

microelectromechanical systems. The operational range of microcantilever beams under electrostatic force can be extended
beyond pull-in in the presence of an intermediate dielectric layer,
which has a significant effect on the behavior of the system.
Three possible configurations of the beam over the operational
voltage range are floating, pinned, and flat configurations. In
this paper, a systematic method for deriving dynamic equation of
microcantilevers for all configurations is presented. First, a static
study is performed on deflection profile of the microcantilever
under electrostatic force. After that, a polynomial approximate shape function with a time-dependent variable for each
configuration is defined. Using Hamiltons principle, dynamic
equations of microcantilever in all configurations are derived. A
comparison between modeling results and previous experimental
data that has been used for validation of the model shows a good
agreement.
[2012-0218]
Index TermsDynamic modeling, electrostatic force, floating,
microcantilever beams, pinned and flat configurations.

I. Introduction

HE ELECTROSTATIC force is often regarded as the

main drive source for microsystems [1]. Electrostatically actuated beams or plates are used in a wide range of
applications such as microswitches, capacitive pressure sensors, accelerometers, micromotors, resonant sensors, and many
others [1], [2]. Modeling and simulation of electrostatically
actuated devices play an important role in the design phase
for predicting device characteristics.
The behavior of electrostatic MEMS parallel plate actuators
before pull-in is studied extensively in the literature [3][11].
Many MEMS devices operate beyond pull-in, e.g., capacitive
switches [12], [13], zipper varactors [14], [15], and tunable
coplanar waveguide (CPW) resonators [16], [17]. Little research has been performed to predict the behavior of the
electrostatically actuated microbeams beyond pull-in.
Gilbert et al. [18], [19] used CoSolve-EM, a coupled
solver for 3-D quasistatic electromechanical problems, to
simulate MEMS actuators beyond pull-in. This technique is
computationally expensive, particularly for problems involving

Manuscript received August 14, 2012; revised December 23, 2012; accepted
March 18, 2013. Date of publication May 22, 2013; date of current version
September 27, 2013. Subject Editor N. Aluru.
The authors are with the Center of Excellence in Design, Robotics
and Automation, Department of Mechanical Engineering, Sharif University
of Technology, 1458889694 Tehran, Iran (e-mail: [email protected];
[email protected]; [email protected]).
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/JMEMS.2013.2256108

mechanical contacts. To find contact length and stiction force

after pull-in, some approaches use static analysis to approximate beam deflection by iterative methods [20][23]. Gorthi
et al. [24] depicted three possible static configurations of the
beam over the total operational voltage range called floating,
pinned, and flat configurations (defined and shown in the next
section); the latter two are referred to as arc-type and S-type
in the literature. Using numerical methods, static governing
equation has been solved numerically in three configurations.
Dynamic stability analyses have been performed for the floating and pinned configurations. Thus, all previous analyses on
the behavior of microcantilevers beyond pull-in are performed
under static conditions.
This paper presents a systematic method for deriving dynamic equations of the electrostatically actuated microcantilevers which is valid over full operational range. In this
method, using static analysis, a deflection profile for each of
the floating, pinned, and flat configurations is selected. The
deflection profiles are a polynomial functions that approximate
the numerical solution of static governing equations for each
configuration. Each configuration has only one time-dependent
variable in the corresponding approximate shape function.
This results in a simplified one-DOF system for each of
configurations. Using Hamiltons principle, the dynamics of
configuration variables can be derived and given the profiles
approximate functions, dynamic behavior of microcantilevers
can be modeled. The simulation results are in good agreement
with the existing experimental data.
The proposed dynamic model is simple, precise, and computationally efficient. The proposed method can also be used for
modeling of microbeams with different boundary conditions
(i.e., simply supported or fixedfixed boundary condition) and
under different forces such as electrostatic, electromagnetic,
capillary, or van der Waals forces.
This paper is presented in six sections. In Section II, we
give the description of microcantilever and its possible configurations. The static and dynamic modeling will be presented
in Sections III and IV and model validation will be given in
Section V. Finally, conclusions are drawn in Section VI.

II. System Description

The electrostatically actuated microcantilever beam is
schematically illustrated in Fig. 1. By applying a voltage, V,
across the plates, an electrostatic attractive force is induced
which leads to a decrease of the air gap. The magnitude of

c 2013 IEEE
1057-7157 

ABTAHI et al.: FULL OPERATIONAL RANGE DYNAMIC MODELING OF MICROCANTILEVER BEAMS

1191

TABLE I
Description of the Variables and Parameters in (1)

Fig. 1.

Schematic of a cantilever beam under electrostatic actuation.

Variable or Parameter
y
x
E
I
0
r
W
L
ti
tp
h

Description
Lateral deflection of the beam
Position along the length
Youngs modulus of the beam material
Moment of inertia of the beam section
Permittivity of free space (8.851012 F/m).
Dielectric constant
Beam width
Beam length
Insulator layer thickness
Plate or beam thickness
Zero bias height
Density of the beam material
Applied voltage

3) Flat configuration: as indicated in Fig. 2(c), in this

configuration, there is a flat part that is in contact with
the dielectric layer. The length of the flat part, l, is the
variable of this configuration. The bending moment at
the separating point of the flat part is zero.

III. Static Study

The static governing equation for the mechanical deformation of an EulerBernoulli beam is as follows [24]:

Fig. 2. Possible configurations of the cantilever beam. (a) Floating configuration. (b) Pinned configuration. (c) Flat configuration.

the electrostatic force is inversely proportional to the square

of the distance. Thus, increasing the applied voltage results
in a dramatic increase in electrostatic force. Since elastic
forces increase linearly with respect to the air gap, a point
occurs where the elastic restoring force can no longer balance
the attractive electrostatic force and the structure becomes
mechanically unstable and the pull-in occurs.
During the pull-in instability, three different phases can be
considered. Each phase is associated with a specific configuration in the boundary conditions at the free end of the cantilever
beam. For each configuration, a variable is considered that
changes by varying the applied voltage. These configurations
are as follows [24].
1) Floating configuration: as indicated in Fig. 2(a), in this
configuration, the cantilever beam has no contact with
the dielectric layer. The deflection of the beam at the free
end, , is the variable of this configuration. The bending
moment and shear force at the endpoint are zero.
2) Pinned configuration: as indicated in Fig. 2(b), in
this configuration, the free end of the cantilever beam
touches the dielectric layer but is free to rotate about
the contact point. The angle between the beam and the
dielectric layer at the contact point, , is the variable of
this configuration. The bending moment at the endpoint
is zero.



4 y
h + y(x)
0 W V 2
1
EI 4 =
1 + 0.65
(1)
x
2 (h + ti + y(x))2
W
r

3
where I = W tp 12 and the right part of the equation is the
electrostatic force per unit length with fringing field correction
factor. The variables and parameters in the aforementioned
equation are described in Table I. The boundary conditions at
the clamped end of the beam for the configurations are
y(0) = 0,

y(0)
= 0.
x

(2)

At the other end of the beam, the boundary conditions for

each configuration are as follows [24]:


2 y1 (L)
3 y1 (L)
=
0,
=0
(3-1)
2
x3
 2 x
y2 (L)
= 0, y2 (L) = h
(3-2)
2
x
2 y3 (L l)
y3 (L l)
Flat
= 0,
= 0.
2
(3-3)
x
x
Config.
y3 (L l) = h
(3)
For the right side of the separation point [see Fig. 2(c)], we
have 2 y3 x2 = 0, so the bending moment at the right side
of the separation point is equal to zero. Thus, the right side
cannot apply or tolerate any bending moment. So the bending
moment at the separation point is equal to zero. Slope at the
separation point is zero because of the vertical upward force
applied to the beam at the separation point and the distributed
electrostatic force applied to the beam. It can be noted that, in
this paper, the dielectric layer is assumed to be rigid. These
Floating
Config.
Pinned
Config.

1192

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 22, NO. 5, OCTOBER 2013

boundary conditions for the flat configuration will change by

assuming the dielectric layer to be an elastic foundation.
The differential equation (1) and boundary conditions (2)
and (3) compose three boundary value problems for three
configurations. For approximating these boundary value problems, at least fourth degree polynomials are necessary. To
increase the precision of the approximation, three fifth degree
polynomials have been considered as the approximate shape
functions in the following:
y1 (x) = a0 + a1 x + + a5 x5 , for 0 x L
y2 (x) = b0 + b1 x + + b5 x5 , for 0 x L
y3 (x) = c0 + c1 x + + c5 x5 , for 0 x L l

(4)

where y1 (x), y2 (x), and y3 (x) denote lateral deflection of the

beam for floating, pinned, and flat configurations, respectively,
and ai , bi , and ci (i = 0, . . . , 5) are unknown constants. By
applying boundary conditions (2) and (3) in (4), we have
y1 (x) = (6a4 L2 + 20a5 L3 )x2
+(4a4 L 10a5 L2 )x3
+a4 x4 + a5 x5
(3b4 L4 + 7b5 L5 3h) 2
y2 (x) =
x
2L2 5
4
(5b4 L + 9b5 L h) 3

x
3
2L
4
5
+b4 x + b5 x
(c5 (L l)5 + 6h)x2
y3 (x) =
(L l)2
(3c5 (L l)5 + 8h)x3
+
(L l)3
3(c5 (L l)5 + h)x4

+ c5 x 5
(L l)4

3(c5 (L l)5 + h)x4

+ c5 x 5 .
(7)
(L l)4
There is one constant for each configuration, a5 , b5 , and c5 ,
that can be used to minimize the error of the approximate
shape function with respect to the accurate solution of the
boundary value problem. The boundary value problem (1) with
boundary conditions (2) and (3) can be solved numerically using boundary value problems solver (i.e., bvp4c or bvp5c
in MATLAB). Now, a5 , b5 , and c5 can be selected to minimize
error between approximate shape functions (7) and numerical
solutions. By trial and error, a5 , b5 , and c5 are selected as

5h + 4L
5h
a5 =
, b5 =
, c5 =
.
(8)
7L5
L5
(L l)5

One can use finite element method (FEM) simulation as

the accurate solution and try to minimize error between
approximate shape functions (7) and FEM simulation results.
By substituting (8) into (7), the approximate shape functions
are then obtain as follows:

for 0 x L

for 0 x L

16x2 2x3
6x4 x5
+ 3
+
2
7L
L
7L4 7L5
(h + L)x2 (7h + 7L)x3
y2 (x, ) =
+
L2
L3
4
(12h + 10L)x
(5h + 4L)x5

+
L4
L5

2
3
hx
7hx

(L l)2
(L l)3

,0 x L l
y3 (x, l) =
12hx4
5hx5

(L l)4
(L l)5

h
,L l < x L
y1 (x, ) =

for 0 x L l.

(5)
The variable of the floating and pinned configuration, and
, respectively, can be defined as follows:
y1 (L) = 3a4 L4 + 11a5 L5 =
y2 (L)
b4 L4 + 3b5 L5 + 3h
=
= .
x
2L

(6)

The approximate shape functions must be a function of the

corresponding configuration variable. In (5), the approximate
shape function for the flat configuration y3 (x) is a function
of its variable l. For the floating and pinned configurations,
(6) can be used to substitute the unknown constants a4
and b4 by unknown variables and , respectively. The
configuration variables are unknown and change by varying
the applied voltage. By substituting (6) in (5), approximate
shape functions become a function of the variables d, a, and
l, so we have
2(a5 L5 + )x2 (14a5 L5 + 4)x3
+
L2
3L3
5
4
(11a5 L + )x

+ a5 x 5
3L4

y1 (x, ) =

(b5 L5 + 6h 3L)x2 (3b5 L5 + 8h 5L)x3

+
L2
L3
(3b5 L5 + 3h 2L)x4

+ b5 x 5
L4
(c5 (L l)5 + 6h)x2 (3c5 (L l)5 + 8h)x3
y3 (x, l) =
+
(L l)2
(L l)3
y2 (x, ) =

(9)
By considering a5 , b5 , and c5 to be zero, fourth degree
approximate shape functions are obtained. Difference between fourth and fifth degree approximate shape functions
and numerical solutions are shown in Fig. 3. The maximum
error of fourth and fifth degree approximate shape functions
with respect to numerical solution is about 15% and 7%,
respectively. In fact, we tried to obtain three approximate
shape functions that are function of x and the corresponding configuration variable while having the minimum error
compared to the numerical solutions. One can use the sixth
or a greater order polynomials to increase the precision but
finding suitable constants (i.e., a6 , b6 , and c6 ) will be more
complicated.
In (9), the beam end deflection, , can vary up to h. At the
boundary between the floating and pinned configurations, by
substituting =h in (9), we have
y1 (x, )|=h =

16hx2 2hx3
6hx4
hx5
+

+
7L2
L3
7L4 7L5 .

(10)

ABTAHI et al.: FULL OPERATIONAL RANGE DYNAMIC MODELING OF MICROCANTILEVER BEAMS

1193

been considered as a one-degree-of-freedom (1-DOF) system

in each configuration. In this section, the dynamic governing
equation for each configuration is derived using Hamiltons
principle. The time derivative of the beam deflection for the
floating configuration can be expressed as
dy1 (x, (t))
=
dt

16x2 2x3
6x4
x5
+ 3
+
2
4
7L
L
7L
7L5

(t).

(15)

It must be noted that the only time-dependent variable in

the floating configuration is the beam end deflection . The
kinetic energy of the beam can be calculated as

1 L
dy1 (x, (t)) 2
T1 =
Wtp (
) dx
2 0
dt
(16)
4492
=
Wtp L 2 (t).
33957
The mechanical bending strain energy of the beam is given
by
Fig. 3. Difference between fourth and fifth degree approximate shape functions and numerical solutions.

For angle of the beam at the endpoint, we have

y1

9h
=
= .
x
= h
7L
x=L
After substituting =

(11)

9h
in (9), we have
7L

hx2
7hx3 12hx4
5hx5

L2
L3
L4
L5
(13)
hx2
7hx3 12hx4
5hx5
y3 (x, l)|l=0 = 2 3 +
5 .
L
L
L4
L
The right-hand side in (13) is again the same. It can be
easily deduced that variables , , and l can take values in the
following range:
y2 (x, )|=0 =

Floating Config. : < h

0

Flat Config. :

0 < l.

9h
7L

1
2

EI(
0

2 y1 2
2904 EI2 (t)
)
dx
=
.
x2
1715 L3

(17)

Hamiltons principle requires specification of the virtual

work. The electrostatic and damping forces can be considered
as distributed external forces. Hence, the work done by the
external forces is
 L
W1 =
f (x, t)y1 (x, (t))dx
(18)
0

6hx4 hx5
16hx2 2hx3
+

+
. (12)
y2 (x, )| 9h =
7L2
L3
7L4 7L5
=
7L
The right-hand sides in (10) and (12) are the same. Similarly, at the boundary between pinned and flat configurations,
we have

Pinned Config. :

U1 =

(14)

IV. Dynamic Modeling

With regards to the beam deflection profile for each configuration as (9), , , and l can be considered as a timedependent variable for floating, pinned, and flat configuration,
respectively. In fact, we assume that the applied voltage to the
microcantilever only provokes the first vibration mode of the
microcantilever in each configuration, so microcantilever has

where f(x, t) is the electrostatic and damping force per unit

length as a function of x and t and can be formulated as
1
f (x, t) = 
2

0 W V 2 (t)

2
ti
h + + y (x, (t))
r


h + y (x, (t))
1 + 0.65
cy (x, (t))
W

(19)

where c indicates the damping coefficient and y(x,(t)) can be

expressed by y1 (x,(t)), y2 (x,(t)), or y3 (x,l(t)) corresponding
to each configuration. Damping coefficient is a simple substitute for damping sources such as squeeze film damping,
structural damping, and so on. However, squeeze film and
adhesion forces can be simply considered in the modeling
process. The Hamilton principle states that
 t2
(T1 V1 + W1 )dt = 0
(20)
t1

where  in the aforementioned equation indicates the variation, operationally equivalent to a total differential. After substituting (16)(18) in (20) and some mathematical operations,
the dynamic equation of motion for the first configuration of
the system can be derived as follows:
3
= 99 5808EI(t) 1715L Q (t)
(t)
44920
L4 Wtp
L
16x2 2x3
6x4
x5
Q (t) = 0 f (x, t)( 2 + 3
+
)dx.
7L
L
7L4 7L5

(21)

1194

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 22, NO. 5, OCTOBER 2013

Similarly, the dynamic equation of motion for the pinned

configuration can be obtained as
110 EI(44L(t) 62h) 7L2 Q (t)
9
L5 Wtp
L
x2 7x3
10x4 4x5
Q (t) = 0 f (x, t)( + 2 3 + 4 )dx.
L
L
L
L
And for the flat configuration, we have
=
(t)

1
(109296EIh2
1646(L l)3 Wtp h2
+823Wtp h2l(t)2 (L l(t))2 4620Ql (L l(t))4 )
 Ll(t)
2hx2
21hx3
f (x, t)(

Ql (t) = 0
3
(L l(t))
(L l(t))4
4
48hx
25hx5
+

)dx.
5
(L l(t))
(L l(t))6

(22)

l(t) =

if( = h and > 0)

(23)

(24)

After this, by increasing the applied voltage, the contact

angle decreases to zero and transition from pinned to flat
configuration occurs. This transition is physically similar to
rolling motion. Thus, in this transition, the kinetic energy
does not change. By decreasing the applied voltage, the beam
returns to pinned and floating configurations. Transition from
flat to pinned configuration and transition from pinned to
floating configuration do not involve impact and thus cannot
change the kinetic energy. Based on these explanations, initial
conditions for the rate of change of the configuration variable
at each transition can be extracted. For the kinetic energy in
each configuration, similar to (16), we have
4492
T1 =
Wtp L 2 (t)
33957
9
Wtp L3 2 (t)
T2 =
(25)
1540
823 Wtp h2l2 (t)
.
T3 |l=0 =
4620
L
For transition from pinned to flat configuration, we have

3
2469 L2
.
(26)
T3 |l=0 = T2 l =
823
h

9h
and > 0)
7L
Pinned Floating

if( = 0 and < 0)

9h
7L
= 0

Floating Pinned
if( =

Considering a beam in the rest position with zero applied

voltage, by increasing the applied voltage, the beam bends
downward. If the applied voltage excels the pull-in voltage, the
tip of the beam will touch the dielectric layer. At this moment,
there is an impact between the beam end and dielectric layer
surface. Since investigation on the impact in this system is
beyond the scope of this paper, the impact has been considered
totally inelastic. This means after the collision, the beam end
sticks to the dielectric surface and will not bounce back. Since
the beam end velocity becomes zero after the collision, (16)
implies that the kinetic energy of the beam after impact will
be zero. Here, 1-DOF modeling causes the kinetic energy to
vanish after the totally inelastic collision. If one models the
microcantilever as a system with more DOFs, some of the
kinetic energy would be conserved after the collision. Thus,
after transition from floating to pinned configuration, the time
derivative of the pinned configuration variable will also be
zero
T2 = T1 = 0 = 0.

For transition from flat to pinned configuration, we have

2469 h l
|
.
(27)
T2 = T3 l=0 =
9
L2
For transition from pinned to floating configuration, we have

63 5615

T1 = T2 =
L .

(28)
22460
Thus for transition between different configurations, the
following conditions and initial values are used:

=h
63
5615
=
L
22460

l=0

3
2469 L2
l =
823
h
if(l = 0 and l < 0)

= 0

Flat Pinned
2469 h l
=
.
9
L2
Pinned Flat

(29)
Using dynamic equations of motion (21)(23), the variation
of the system variables d, a, and l can be calculated in the
range specified by (14). Also, the conditions in (29) can be
used for switching between the configurations. By using the
approximate shape functions in (9), the behavior of the system
can be simulated.

V. Results
It this section, simulation results of the model is compared
with published data to validate the model. The RungeKutta
method is adopted for solving the resulting differential equation of motion.
The first example is a polysilicon cantilever beam over
a substrate without dielectric layer that has been analyzed
dynamically via full-Lagrangian schemes in [10]. The beam
is 80 m long, 0.5 m thick, and 10 m wide. The initial
gap between the beam and the substrate is 0.7 m. Youngs
modulus of 169 GPa and a mass density of 2231 kg/m3 are
employed in the simulations. The damping factor c is taken to
be zero.
Fig. 4 shows the cantilever beam end deflection for an
applied step voltage of 0.5 V that has been compared with
modeling results reported in [10]. As seen from this graph, a
very good agreement is achieved about the resonant frequency
and maximum deflection. Fig. 5 shows the dynamic pull-in of
the cantilever beam at an applied step voltage of 2.15 V. A
completely similar pull-in behavior is reported in [10] but at
the applied voltage of 2.12 V. There is less than 2% difference
between the obtained pull-in voltages.

ABTAHI et al.: FULL OPERATIONAL RANGE DYNAMIC MODELING OF MICROCANTILEVER BEAMS

Fig. 4.

Cantilever beam end deflection by a step voltage of 0.5 V.

Fig. 5.

Pull-in behavior of a beam by a step voltage of 2.15 V.

The second example is a polysilicon cantilever beam with

experimental results reported in [11]. The beam is 20 mm long,
57 m thick, and 5 mm wide. The initial gap between the beam
and the substrate is 92 m. Youngs modulus of 155.8 GPa and
a mass density of 2330 kg/m3 are employed in the simulations.
The damping factor c is 1.57 kg/ms that is about the critical
damping and eliminates oscillations in the response. Fig. 6
shows the final cantilever beam end deflection for various
applied voltages that is compared with experimental results
reported in [11]. The results obtained by the present model
completely match the experimental data (less than 3% error).
The third example is a series of the polysilicon cantilever
beams shown in Fig. 7 over the substrate with silicon nitride
dielectric layer that have been tested in [20]. The beams width
and thickness are 10 m and 2 m, respectively, and lengths
range from 200 to 550 m. The initial gap between the beams
and the substrate is 2 m. Youngs modulus of 154 GPa and a
mass density of 2330 kg/m3 are employed in the simulations.
The damping factor c is taken to be zero. The dielectric layer
thickness is 20 nm and its constant is 5.7. In Fig. 8, the model
predictions for pull-in voltages are compared to the measured
pull-in voltages for cantilevers of different lengths reported
in [20]. There is a good agreement (less than 12%) between
model predictions and experimental results.

1195

Fig. 6.

Beam end deflection by applying various voltages.

Fig. 7.
[20].

Array of cantilevers of various lengths used in the measurements

Fig. 8.

Pull-in voltages for different cantilever lengths.

After pull-in has occurred, by increasing the voltage, more

of the beam touches the dielectric layer. The percentage of the
contact length relative to the beam length was measured for
various voltages in [20]. In this case, the cantilever length
is 450 m and the damping factor c is considered to be
0.01 kg/ms that is about critical damping and eliminates

1196

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 22, NO. 5, OCTOBER 2013

Fig. 9. Percentage of the contact length relative to the beam length for
various applied voltages.
TABLE II
Parameters of the Modeled Cantilever Beam
L = 450 m
W = 10 m
E = 154 GPa
ti = 20 nm
0 = 1

Fig. 10.

tp = 2 m
h = 2 m
= 2330 kg/m3
r = 5.7
c = 0.01 kg/ms

oscillations in the response. The parameters of this cantilever

beam are summarized in Table II. Fig. 9 shows a quite good
agreement (less than 7% error) between the model results and
experimental data.
The good agreement between the model results and previous
experimental or analytical data shows that the model is reliable
and can be used to extract the dynamic behavior of the system.
Now, considering the cantilever beam described in Table II,
Fig. 10 shows the variation of three configuration variables for
a predefined voltage profile with different maximum voltage
values. As seen from the graph, there are different pull-in times
for different voltage values. For maximum voltage value equal
to 3.3 and 3.5 V, after entering the flat configuration, contact
length oscillations damp out rapidly. For better demonstration,
Fig. 11 shows the zoomed variation of and l. As seen
from the graph with the maximum voltages equal to 3.103
V, after entering the pinned configuration, the contact angle
will decrease to near zero and after removing the voltage,
microbeam will return to the floating configuration.
Gorthi et al. [24] defined a nondimensional quantity
ti
h =
.
r h

(30)

They showed that voltage limits of different configurations

depend on the h value. Considering the cantilever beam
described in Table II, voltage limits of different configurations
for different insulator layer thicknesses, ti , has been obtained
from the simulation results of the present model. Table III
shows this voltage limits.
These results (obtained from simulations) are in complete
agreement with the results in [24]. The current microcantilever

(ti =0.02 m, h=0.0018)

exhibits bistability between the floating

Applied voltage graph and variation of configuration variables.

Fig. 11. Applied voltage graph and variation of configuration variables

(zoomed view).

TABLE III
Voltage Limits of Different Configurations for Different
Insulator Layer Thicknesses Obtained From the Simulation
Results
ti (m)
h
max
Vfloat
min
Vpin
max
Vpin
min
Vflat

[24]

0.02
0.0018
2.71
0.20
1.94

0.5
0.0439
2.87
0.99
3.47

0.8
0.0702
2.99
1.30
4.19

0.45
0 < h < 0.03
max
max
Vfloat
> Vpin

2.49
0.03 < h < 0.07
max
max
Vfloat
< Vpin

3.25
0.07 < h < 0.40
max
max
Vfloat
< Vpin

min < V max

Vflat
float

min < V max

Vflat
float

min > V max

Vflat
float

ABTAHI et al.: FULL OPERATIONAL RANGE DYNAMIC MODELING OF MICROCANTILEVER BEAMS

1197

voltage). Given the small damping coefficient, the damped natural frequency is used as an estimate of the natural frequency.
This routine can be simulated to obtain the natural frequency
of the beam for various applied voltages in different configurations. Fig. 13 shows the natural frequency of the beam
with respect to the natural frequency of the unactuated beam
for various applied voltages in different configurations. In the
floating configuration, as the voltage is increased and the beam
deflects toward the substrate, the natural frequency of the beam
around the deflected position starts to decrease. In the pinned
configuration, there is a decrease in natural frequency when
the applied voltage is close to the maximum voltage applicable
max
. In the flat configuration, as
in the pinned configuration, Vpin
the voltage is increased, the contact length and subsequently
natural frequency of the beam increase steadily.

VI. Conclusion

Fig. 12. Applied voltage graph and variation of configuration variables with
repeated 1.7 V.

This paper presented a routine method for deriving full

operational range dynamic equations of the microcantilevers
under electrostatic force. Using Hamiltons principle, dynamics of configuration variables are derived. Using dynamics of
variables and approximate shape functions, the dynamic behavior of microcantilevers is simulated. Comparison between
the modeling results and previous experimental data showed
a good agreement. The proposed dynamic model is simple,
precise, and computationally efficient. The proposed method
can also be used for modeling of microbeams with different
boundary conditions (i.e., simply supported or fixedfixed
boundary condition) and under different forces such as electrostatic, electromagnetic, capillary, or van der Waals forces.

References

Fig. 13. Natural frequency of the beam for various applied voltages in
different configurations.

and pinned configurations in the voltage range 0.20.45 V,

tristability among all three configurations in the voltage range
0.451.94 V, and bistability between the pinned and flat
configurations in the voltage range 1.942.71 V. Fig. 12 shows
the variation of three configuration variables for another predefined voltage profile in which voltage equals to 1.7 V at three
different time periods. This figure shows that the microcantilever is stable in three different configurations for voltage values in the tristable range when applied voltage is well-defined.
To obtain the natural frequency, a constant voltage is applied, and the beam is allowed to achieve a stable equilibrium
configuration. After stabilization, beam will vibrate at its
damped natural frequency if the applied voltage suddenly
changes by a small amount (less than 0.1% of the applied

[1] V. Ostasevicius and R. Dauksevicius, Microsystems Dynamics (International Series on Intelligent Systems, Control, and Automation: Science
and Engineering), vol. 44, S. G. Tzafestas, Ed. The Netherlands:
Springer, 2011, pp. 34.
[2] M. Lishchynska, N. Cordero, O. Slattery, and C. OMahony, Modelling
electrostatic behaviour of microcantilevers incorporating residual stress
gradient and non-ideal anchors, J. Micromech. Microeng., vol. 15, pp.
1014, Jul. 2005.
[3] L. A. Rocha, E. Cretu, and R. F. Wolffenbuttel, Full characterisation
of pull-in in single-sided clamped beams, Sens. Actuators A, vol. 110,
pp. 301309, Feb. 2004.
[4] N. Agarwal and N. R. Aluru, Stochastic analysis of electrostatic MEMS
subjected to parameter variations, J. Microelectromech. Syst., vol. 18,
pp. 14541468, Dec. 2009.
[5] S. Chaterjee and G. Pohit, A large deflection model for the pull-in
analysis of electrostatically actuated microcantilever beams, J. Sound
Vibrat., vol. 322, pp. 969986, May 2009.
[6] G. M. Rebeiz, RF MEMS: Theory, Design and Technology. Hoboken,
NJ, USA: Wiley, 2003.
[7] A. H. Nayfeh, M. I. Younis, and E. M. Abdel-Rahman, Dynamic pullin phenomenon in MEMS resonators, Nonlinear Dynam., vol. 48, nos.
12, pp. 153163, 2007.
[8] X. Zhao, E. M. Abdel-Rahman, A. H. Nayfeh, A reduced-order model
for electrically actuated microplates, J. Micromech. Microeng., vol. 14,
pp. 900906, Jul. 2004.
[9] J. A. Pelesko, Mathematical modeling of electrostatic MEMS with
tailored dielectric properties, SIAM J. Appl. Math., vol. 62, pp.
888908, Feb. 2002.
[10] S. K. De and N. R. Aluru, Full-Lagrangian schemes for dynamic
analysis of electrostatic MEMS, J. Microelectromech. Syst., vol.13, pp.
737758, Oct. 2004.

1198

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 22, NO. 5, OCTOBER 2013

[11] Y. C. Hu, C. M. Chang, and S. C. Huang, Some design considerations

on the electrostatically actuated microstructures, Sens. Actuators A, vol.
112, pp. 155161, Apr. 2004.
[12] J. B. Muldavin and G. M. Rebeiz, High-isolation CPW MEMS shunt
switches. 1. Modeling, IEEE Trans. Microw. Theory Tech., vol. 48, no.
6, pp. 10451052, Jun. 2000.
[13] H. S. Newman, RF MEMS switches and applications, in Proc. 40th
Annu. Int. Reliab. Phys. Symp., Apr. 2002, pp. 111115.
[14] R. Mahameed and G. M. Rebeiz, Electrostatic RF MEMS tunable
capacitors with analog tunability and low temperature sensitivity, in
IEEE MTT-S Int. Microw Symp. Dig., May 2010, pp. 12541257.
[15] G. V. Ionis, A. Dec, and K. Suyama, A zipper-action differential micromechanical tunable capacitor, in Proc. Microelectromech. Syst. Conf.,
Aug. 2001, pp. 2932.
[16] E. Fourn, A. Pothier, C. Champeaux, P. Tristant, A. Catherinot, P.
Blondy, G. Tann, E. Rius, C. Person, and F. Huret, MEMS switchable
interdigital coplanar filter, IEEE Trans. Microw. Theory Tech., vol. 51,
no. 1, pp. 320324, Jan. 2003.
[17] T. Ketterl, T. Weller, and D. Fries, A micromachined tunable CPW
resonator, in IEEE MTT-S Int. Microw. Symp. Dig., vol. 1. May 2001,
pp. 345348.
[18] J. R. Gilbert, R. Legtenberg, and S. D. Senturia, 3D coupled electromechanics for MEMS: Applications of CoSolve-EM, in Proc. Int. Conf.
Microelectromech. Syst., Jan.Feb. 1995, pp. 122127.
[19] J. R. Gilbert, G. K. Ananthasuresh, and S. D. Senturia, 3D modeling
of contact problems and hysteresis in coupled electro-mechanics, in
Proc. 9th Int. Workshop Microelectromech. Syst. Conf., Feb. 1996, pp.
127132.
[20] T. Lam and R. B. Darling, Physical modeling of MEMS cantilever
beams and the measurement of stiction force, in Proc. Model. Simul.
Microsyst. 2001, pp. 418421.
[21] J. A. Knapp and M. P. de Boer, Mechanics of microcantilever beams
subject to combined electrostatic and adhesive forces, J. Microelectromech. Syst., vol. 11, pp. 754764, Dec. 2002.
[22] Z. Yin and Z. Ya-pu, Static study of cantilever beam stiction under
electrostatic force influence, Acta Mech. Sol. Sinica, vol. 17, pp.
104112, Jun. 2004.
[23] S. Basu, A. Prabhakar, and E. Bhattacharya, Estimation of stiction
force from electrical and optical measurements on cantilever beams,
J. Microelectromech. Syst., vol. 16, pp. 12541262, Oct. 2007.

[24] S. Gorthi, A. Mohanty, and A. Chatterjee, Cantilever beam electrostatic

MEMS actuators beyond pull-in, J. Micromech. Microeng., vol. 16, pp.
18001810, Sep. 2006.
Mansour Abtahi received the B.S. degree in mechanical engineering from
Isfahan University of Technology, Isfahan, Iran, in 2006, and the M.S. degree
in mechanical engineering from Sharif University of Technology, Tehran,
Iran, in 2008. He has been working toward the Ph.D. degree in mechanical
engineering at Sharif University of Technology since 2009.
His current research interests are focused on MEMS and dynamic modeling.
Gholamreza Vossoughi received the B.S., M.S., and Ph.D. degrees from the
Mechanical Engineering Department, University of Minnesota, Minneapolis,
MN, USA.
He joined Sharif University of Technology, Tehran, Iran, in 1992 where
he is presently a Professor of mechanical engineering. His research interests
are bio inspired robotics, manmachine interface and haptic systems, mechatronics, micro-mechatronics, and control systems.
Ali Meghdari received the Ph.D. degree in mechanical engineering from the
University of New Mexico, Albuquerque, NM, USA, in 1987.
He then joined the Robotics Group of Los Alamos National Laboratory
(LANL) as a Research Collaborator. In 1988, he transferred to Sharif
University of Technology (SUT), Tehran, Iran. In 1997, he was granted the
rank of Full Professor of mechanical engineering at SUT, and was recognized
by the Iranian Society of Mechanical Engineers (ISME) as the youngest (at age
37) ME Full Professor in Iran. From 2001 to 2010, he was the Vice President
of Academic Affairs, and the Director of the Center of Excellence in Design,
Robotics, and Automation (CEDRA) at Sharif University of Technology. He
has performed extensive research in the areas of robotics kinematics and
dynamics, nano-robotics and nano-manipulation, and dynamic modeling of
biomechanical systems.
Dr. Meghdari was nominated and elected as a Fellow of the American
Society of Mechanical Engineers (ASME) in 2001.