Sensors and Actuators APhysical Thermal cycle modeling of

Document Sample
Sensors and Actuators APhysical Thermal cycle modeling of Powered By Docstoc
					                                                                    Sensors and Actuators A 152 (2009) 192–202

                                                                  Contents lists available at ScienceDirect

                                                      Sensors and Actuators A: Physical
                                                     journal homepage:

Thermal cycle modeling of electrothermal microactuators
Mohammad Mayyas a,∗ , Panos S. Shiakolas b , Woo Ho Lee a , Harry Stephanou a
    Automation & Robotics Research Institute, University of Texas at Arlington, Jack Newell Blvd. S., Fort Worth Texas 76118, USA
    Department of Mechanical and Aerospace, University of Texas at Arlington, 500 West First Street, 211 Woolf Hall, Arlington, TX 76019-0018, USA

a r t i c l e           i n f o                            a b s t r a c t

Article history:                                           A methodology for computing the transient heating (charging) and cooling (discharging) time for a serially
Received 25 May 2008                                       connected electro thermo (E-T) microbeams actuator were analytically approximated from fundamental
Received in revised form 4 February 2009                   characteristics including input voltage, material and geometry. Specifically, this paper presents closed
Accepted 5 March 2009
                                                           form solutions for folded and bend beam MEMS actuators. The E-T heating time was obtained from the
Available online 2 April 2009
                                                           simplified transient heat conduction equations using a trial solution method combined with the Laplace
                                                           domain based on polynomial shape function, while the cooling response of thermal actuator is derived
                                                           from steady state temperature distribution. The proposed approximate methods are utilized in illustrative
Thermal cycle
                                                           examples and validated by Finite Difference Approximation (FDA), and Finite Element Modeling (FEM).
Actuators                                                                                                                                Published by Elsevier B.V.

1. Introduction

    Electrothermal (E-T) micromechanical actuation based on asymmetrical thermal expansion has popularly been employed as one of the
major actuation principles in microsystems. Thermally stable E-T actuators have a dynamic cycle that consists of heating, dwell (engaging)
and cooling times. The maximum operating frequency of E-T actuators, excluding the dwell time, is measured from the thermal heating and
cooling rates. Traditionally, MEMS actuators have faster response times compared to macro scale actuators [1]. It is desired to periodically
drive thermal actuators for power input cycle time less than or equal to thermal cycle (full duty) [2] to provide enough time for the actuator
to dissipate the heat generated by the high frequency and amplitude input current/voltage. If the actuator does not completely dissipate the
heat, it will either retain a perturbed static deflection as if it were driven by a DC input or in a worst case scenario the heat will accumulate
in the structural layers causing the temperature to keep rising and eventual thermal failure. Thus, the structural dynamic response of an
E-T microactuator is of interest and is largely determined by its heating and cooling time.
    In case of E-T actuators with a narrow air gap, fast cooling can be achieved by heat conduction across a substrate; for example, in E-T
vibromotor the bandwidth of thermal–elastic response can be up to a 1.0 kHz range at a low voltage input [3]. It is also reported that the
bandwidth of a lateral thermal actuator fabricated from poly silicon is up to few dozen of kHz at resonance [4,19]. The thermal response of
MEMS devices is usually computed using numerical and lumped models [5], whereas little attention has been given to the measurement of
the transient thermal responses due to the limitation of current sensing capability [1]. A SPICE lumped model was developed to incorporate
electrical loading, transient responses, and deflections of a polysilicon thermal actuator. The model predicts insightful overall performances,
but accurate results require dense network fragmentation [2]. The reduced order dynamic model of a MUMPs thermal actuator was obtained
using experimental data and subsequent used in finite element (FE) simulation demonstrated high speed actuation up to hundreds of kHz
with input shape control [6]. Other previous efforts in modeling the thermal response of E-T actuators include the simplified partial
differential equation (PDE)/ordinary differential equation (ODE) formulation of a lineshape and U-shape microstructure [1,7,8] and the
finite element model (FEM) of thermal microgrippers and bent beams [9,23]. Nodal analysis method by fit and Fourier transformation has
been implemented to calculate the temperature at beam’s nodes and the thermo mechanical coupling by using circuit simulation software
[20,21]. Although the aforementioned techniques relate to thermal cycle analysis, the thermal transient analysis of complex microstructures
is difficult to obtain analytically and in addition, it requires time consuming computation. Furthermore, they may not clearly illustrate a
way to relate efficient high speed performance to transient thermal conditions.

    ∗ Corresponding author. Tel.: +1 8172725887.
      E-mail address: (M. Mayyas).

0924-4247/$ – see front matter. Published by Elsevier B.V.
                                                                M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202                   193

    This paper presents the development of a methodology based on a quadratic approximation function that enables the evaluation of
the transient thermal cycle behavior for an arbitrary number of serially connected microbeam E-T actuators. The developed methodology
is also evaluated for bent and folded beam actuators by comparing it with FDA, FEM and the experimentally obtained results.

2. Thermal cycle model

2.1. Approximate ramping time

   In practice, the thermal actuator should be operated below a frequency that it neither causes thermal failure nor a static deflection offset.
The fact that the mechanical/structural frequency is much faster than thermal frequency can be realized by obtaining analytical solutions
which characterize both the heating and cooling times of an E-T microactuator made of a series of non-uniform cross-section microbeams.
The resistive joule heating of an element in an E-T device is equal to the heat conducted in the element minus the losses out of the element.
The resulting thermal distribution will cause expansion in the connected beams causing the actuator to deflect according to its structural
configuration. The modeling of a thermal actuator can be simplified as a one-dimensional problem when the length scale of a beam is
much greater than its width and thickness, i.e., the Biot number Bi = (tv kp )/(wkv ) less than 0.0002 [21]. tv and w are the thicknesses of
the air gap and width of microbeam. k and kv are the thermal conductivities of the microbeam and air, respectively. In case of a packaged
micro device operating at low temperature, radiation and convection heat losses are negligible compared to the conduction heat loss to the
substrate [15]. The heat lost to the surroundings by radiation found to be less than 1% of the even at high temperature [22]. Meanwhile the
thermal radiation between suspended structure and substrate is negligible for separation of 1 m and above. However, heat conduction
through airgap becomes insignificant under several tens of nanometers [21]. It is assumed that the heat conduction coefficient of suspended
structure is constant (kp = 0 ) and not a function of temperature, and that the temperature of pads and substrate Ts is constant. Thus, the
temperature profile of an actuator due to an input current, I, could be evaluated or estimated by assuming a system of linear PDE for each
microbeam segment. A series of rectangular beams, with constant height, h, and varying width, w, float on top of a substrate with an air
gap and their two ends are anchored at the pads. The input current density is j = I/wh. For an actuator with n arbitrary links (microbeams),
the rate of heat change per area is equal to heat generation per unit volume per time minus heat losses according to (1)
              ∂T              I                                                      d2 T      S   T − Ts
       d Cp      dx =                   o       1 + (T − Ts ) dx − −             o        dx +              dx + (whcu + 2hhcs )(T − Tamb )dx   (1)
              ∂t             wh                                                      dx2       h     RT

where the shape factor S = h/w(2tv /h + 1) + 1 amplifies the heat flow into the substrate [7,11]. In general, a floating microbeam will exhibit
heat loses in the form of conduction and convention and radiation at the upper and lower faces. However, the thickness of the air gap
between the microbeam lower face and the substrate dictates the mechanism of heat loss at the lower face. In this analysis, the airgap
thickness is considered to be in the range of 1–10 m indicating that there is conductive and not convective heat loss [7,21]. Also, the
temperature differences between the microbeam and the substrate and the ambient are small and heat losses due to radiation are also not
    Thus, the conductive thermal resistance between the lower face of the microbeam and the substrate Ri is given by RT = tv /kv + ts /ks ,
where ts and ks are the thicknesses of the substrate and the thermal conductivities of the air and the substrate, respectively. Convection
is considered on the upper face of a thermal actuator with the convection coefficient hcu and film temperature Tfilm . While convection
through side walls can be negligible for a thermal actuator with a thin device layer like those fabricated using PolyMUMPs, the convection
heat transfer of a thick device layer might be considered. It is assumed that the resistivity of the microbeam material is linear with respect
to temperature changes with thermal coefficient . The thermal conductivity of the beam material o is assumed to be constant. The lower
surface substrate temperature is assumed to be constant during short actuation times and equal to ambient Tamb = Ts .
    The transient parabolic heat conduction PDE of an element in Eq. (1) is expressed as,
      ∂2    1 ∂Â
           =       + εÂ.                                                                                                                        (2)
      ∂x2    ˛p ∂t
      Â = (T − T∞ )          ˛p =
                                        d Cp
                                                   S           he       J2   o
      he = whcu + 2hhcs             ε=                     +        −                .                                                          (3)
                                                 hRT   o        o            o
                        J2   o
      T∞ =      Ts +
                        ε    o

    The temperature distribution of a folded beam microactuator consisting of three serially connected beams (hot, cold, and flexure) must
be continuous starting from the right side pad (hot arm side) and ending at the left side pad (flexure arm side) as shown in Fig. 1.
    For continuous heat flux through the arms, the following Dirichlet conditions are necessary but not sufficient. The Laplace transformation
of initial and boundary conditions with H = {Â} gives
      ⎪x = 0 :                         Ts − T∞h
      ⎪                           Hh =
      ⎪                                    s
      ⎪                                      T∞c − T∞h
      ⎪x = L +g :
      ⎨                           Hh − Hc =            ;                     wh
                                                                                       − wc
                                                 s                                  dx       dx
                                                                                                                .                               (4)
      ⎪ x = L + g + Lc : Hh − Hf (L + g + Lc , s) = T∞f                          − T∞c          dHc      dHf
      ⎪                                                                                ;    wc      − wf     =0
      ⎪                                                                          s               dx       dx
      ⎪                       Ts − T∞f
      ⎩ x = 2L + g :     Hf =
194                                                       M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202

                     Fig. 1. Schematic of an attached thermal beam actuator. (a) Top view and (b) cross section with a single layer handle substrate.

   The major difficulty in solving the Laplace form of Eqs. (2) and (4) is simultaneously solving the PDEs of three systems. The changes
in conditions (ε < 0,ε = 0 and ε > 0) result in several temperature responses [13]. A common case is the exponential of low magnitude
temperature profile occurring at low excitation signal with ε > 0. At higher magnitudes, temperature changes its profile at ε = 0, and becomes
sinusoidally distributed for ε < 0. In such wavy profile, high order numerical approximations are required to avoid the suppression and
numerical oscillation of temperature distribution. Hence, an approximate solution is proposed that can optimize and fit different conditions.
   The temperature boundary conditions (BCs) for all the beams are T(x,0) = Ti . The temperature boundary conditions at the pads are not
necessarily the same as the initial conditions for all the beams. These BCs are transformed from the time domain into the Laplace domain,
in order to transform them into only functions of x. The ordinary differential equations of the three systems are rewritten in the Laplace
domain as
      d2 Hk (x, s)       s + ˛p εk                     Ti
                   −                     Hk (x, s) +      = 0 . . . k = {h, c, f }                                                                      (5)
          dx2               ˛p                         ˛p

    A trial solution in terms of finite sum of basis-functions is formulated using the independent space variable x, time phase lag, and the
indeterminate coefficients. A second-order polynomial approximation is chosen as the basis function for the current analysis as shown in

      Hk (x, s, T∞k , ˛, ε) = a0k + a1k x + a2k x2 . . . k = {h, c, f }
      ˜                                                                                                                                                 (6)

   The choice of a second-order approximating polynomial function is due to the fact that it can approximate the exponential temperature
function in the x-domain as it will be shown later in the manuscript. Also, note that higher order approximating polynomials were considered
for actuators at operating condition without any substantial improvement in the results. The exponential to combined (exponential and
sinusoidal) temperature profiles behavior which normally encompasses the majority of the operation E-T actuators.
   The residuals are obtained by substituting Eq. (6) into Eq. (5) for hot, cold, and flexure arms

      Rh (x, s) = Ti /˛p + a0h «1h + a1h «2h + a2h «3h               L+g ≥x ≥0
      Rc (x, s) = Ti /˛p + a0c «1c + a1c «2c + a2c «3c               L + g + Lc ≥ x ≥ L                                                                 (7)
      Rf (x, s) = Ti /˛p + a0f «1f + a1f «2f + a2f «3f               2L + g ≥ x ≥ L + g + Lc

where the coefficients of undetermined residuals are
                     s + ˛p εk           ⎪
       «1k =     −                       ⎪
                        ˛p               ⎪
                     s + ˛p εk
       «2k =     −             x             · · ·k = {h, c, f }                                                                                        (8)
                        ˛p               ⎪
                       s + ˛p εk 2       ⎪
       «3k =     2−
                                x        ⎭

   The next step is to find the best set of parameters of the basis function approximation, which minimize the error between the approxi-
mated solution of the transformed independent variable and the exact solution of the transformed dependent variable. This is performed
using weighted residual method which includes the collocation, subdomain, least square and Galerkin criteria [10].
   The subdomain optimization method is utilized to generate the remaining equations that are needed to solve the undetermined coeffi-
cients of an approximated function. This method is based on setting the residual integrals to zero over each of the subdomains. A reasonable
choice is to divide the domain into intervals based on the length of the hot, cold, and flexure arms
      (L+g)               (L+g+Lc )                  (2L+g)

            Rh dx = 0                 Rc dx = 0                 Rf dx = 0.                                                                              (9)
        0                   (L+g)                   (L+g+Lc )
                                                         M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202                                           195

   The undetermined coefficients are evaluated by substituting the approximating H into both the auxiliary conditions in Eq. (4) and
optimization criteria in Eq. (7)
             ⎡     ⎤−1
      [a] = ⎣ . ⎦
              .           [b] = A−1 [b]                                                                                                                         (10)

where the undetermined coefficients of approximating functions for the three microbeams are a = [a0h a1h a2h a0c a1c a2c a0f a1f a2f ]T and the
constant matrix b in spatial and Laplace domain is given by
             Ts − T∞h          T∞c − T∞h           T∞f − T∞c            Ts − T∞f        Ti (L + g)        Ti (Lc )       Ti (L − Lc )
      b=                                       0                  0                 −                 −              −                                          (11)
                 s                 s                    s                  s                ˛p              ˛p                ˛p

   The boundary conditions and flux continuity between the microbeams yield matrix Aa as
             ⎡                                                                                                                                        ⎤
             1    0                0      0     0                             0                  0         0                            0
           ⎢ 1 (L + g)          (L + g)2 −1  −(L + g)                     −(L + g)               0         0                            0              ⎥
           ⎢                                                                                                                                           ⎥
           ⎢0     1             2(L + g) 0  −(wc /wh )                −2(L + g)(wc /wh )         0         0                            0              ⎥
      Aa = ⎢                                                                                                                                           ⎥.     (12-a)
           ⎢0     0                0     1 (L + g + Lc )                (L + g + Lc )2          −1    −(L + g + Lc )             −(L + g + Lc )2       ⎥
           ⎣0     0                0      0     1                       2(L + g + Lc )           0     −(wf /wc )             −2(L + g + Lc )(wf /wc ) ⎦
             0    0                0      0     0                             0                 1       (2L + g)                    (2L + g)2

  The subdomain optimization method in Eq. (9) completes the indeterminism of the approximate function coefficients. Consequently,
matrix Ab is expressed both in spatial and Laplace domains and given by
                  (L + g)                       (L + g)2                            (L + g)3
             ⎢−     ˛p
                          (s + ˛p εh )     −
                                                         (s + ˛p εh )   2L + 2g −
                                                                                             (s + ˛p εh )
      Ab =   ⎢            0                             0                              0
                            0                           0                                 0

                 0                                     0                                                0
              Lc                      (L + g + Lc )2 − (L + g)2                          (L + g + Lc )3 − (L + g)3
             − (s + ˛p εc )         −                           (s + ˛p εc )       2Lc −                           (s + ˛p εc )
              ˛p                                2˛p                                                3˛p
                 0                                     0                                                0

                      0                                     0                                                             0

                      0                                     0                                                             0
                 (L − Lc )                     (2L + g)2 − (L + g + Lc )2                                    (2L + g)3 − (L + g + Lc )3                   ⎥
             −             (s + ˛p εf )    −                              (s + ˛p εf )        2(L − Lc ) −                              (s + ˛p εf )      ⎥   (12-b)
                    ˛p                                   2˛p                                                           3˛p                                ⎦

    The [Aa |Ab ]T matrix is the coefficient matrix for the system of ODEs that describe the thermal behavior of an E-T actuator. The eigenvalues
of the coefficient matrix yield the time constants of the thermal actuator

        (s) = |sI − a| = |sI − [Aa |Ab ]T | = 0.                                                                                                                (13)

   Each of the approximate polynomial function describes the general solution of the corresponding microbeams in both space and Laplace
domains. Thus, the characteristic Eq. in (13) is of k order and the roots give a representation of the dynamical responses of each beam. In
a stable thermo-dynamic system of folded beam actuator, the real part of the roots of the characteristic equation is the inverse of time
constants 1/ t1 , 1/ t2 and 1/ t3 , these three fundamental time constants determine the decay time in the exponential components of the
time varying coefficient, where it is expected that the slowest time constant occurs as a result of the time response at the hottest arm.
   The general solution for the temperature distribution of the E-T actuator is obtained by computing matrix a in the time domain and
substituting it into the approximate shape function. The coefficient vector a is obtained in time domain by taking the inverse Laplace of the
right-hand side of Eq. (10)
               ⎧⎡       ⎤−1  ⎫
               ⎪ Aa (s)
               ⎨             ⎪
      [a(t)] =  ⎣ .. ⎦ [b(s)] .                                                                                                                                 (14)
               ⎪   .
               ⎩ A (s)       ⎪
196                                                                M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202

   The general solution consists of both the transient and steady state temperature distributions which are obtained by substituting Eq.
(14) in the temperature model

         Tk (x; t) = T∞k + aok (t) + a1k (t)x + a2k (t)x2 . . . k = {h, c, f }
         ˜                                                                                                                                                                                                                                (15)

   Although the presented analysis considered only three sequentially connected folded beams, this approach can be easily extended to
any finite number of connected microbeams. In such case, the matrices in Eqs. (12-a) and (12-b) are obtained by introducing and expanding
additional B.C.s between the microbeams and relaxing the residuals for each one. The matrices Aa and Ab can be constructed for any k
serially connected microbeams with the temperature approximated to a polynomial function of order n, where the generalized matrix a is
square positive definite of size (k × n) × (k × n).

2.2. Approximate cooling time

   The steady state heat conduction equation (S.S.H.C.E.) is obtained by dropping the time partial derivative in Eq. (1). Such simple linear
ordinary differential equations (ODEs) are found in many literatures [12,13]. The temperature distribution of steady state condition in hot,
cold, and flexure arms are easily derived for ε > 0
                          √                 √
       Tk (x) = ˛k + C1k e ˇk / o x + C2k e− ˇk / o x . . . k = {h, c, f }                                                              (16)

where ˇ = (S/hRT + he − J2 o ), ˛d = Ts + J2 o /ˇ, and (C1h , . . ., C2f ) are constants in the temperature solution found by solving steady state heat
conduction equation. A microactuator reaches stable steady state in Eq. (16) if {∀ˇk > 0}. For the three beams (normally occur at low input
current [13], i.e., the maximum allowed current that produces exponential steady state response is

                                 kv ks wi   0 (2tv   + h + wi ) + (kv ts + ks tv )(wi hcui + 2hhcsi )
         I = min      Ii =                                                                                                                k = {h, c, f }                                                                                  (17)
                                                          h 0 (kv ts + ks tv )

   The temperature distribution of a three beam microstructure is continuous and it starts from pad temperature at hot arm and ends at
pad temperature at the flexure arm. The heat flux is continuous at the arms (microbeam interface) and the following Dirichlet conditions
are used to solve for the six unknown constants
         ⎡                                            ⎤
             Th (0) = Ts                                 ⎡                                          ⎤
         ⎢ T (2L + g) = Ts                            ⎥ wc dTc (L + g + Lc ) = wf dTf (L + g + Lc )
         ⎢ f                                          ⎥ ⎣         dx                     dx         ⎦
         ⎢                                            ⎥,                                                                                                                                                                                  (18)
         ⎣ Th (L + g) = Tc (L + g)                    ⎦ w dTh (L + g) = w dTc (L + g)
                                                                   h                                   c
                                                                                 dx                              dx
             Tc (L + g + Lc ) = Tf (L + g + Lc )

      Applying boundary conditions into S.S.H.C.E gives unknown constants C

                      ⎡                                                                                                                                                                                                             ⎤−1
                               1               1                                  0                                        0                                             0                                       0
                  ⎢ √ˇ / o (L+g)          √                                           ˇc                             √                                                                                                                ⎥
                  ⎢e h                  e− ˇh /    o (L+g)              −e                 (L+g)                  −e− ˇc /      o (L+g)                                  0                                       0                    ⎥
                  ⎢                                                                    o                                                                                                                                              ⎥
         ⎡C ⎤ ⎢ √                                                                                                                                                                                                                     ⎥
            1h    ⎢ e ˇh / o (L+g)        √                        wc        ˇc              ˇc             wc        ˇc −√ˇc /                                                                                                       ⎥
           C      ⎢                    −e− ˇh /      o (L+g)   −                e                  (L+g)
                                                                                                                         e             o (L+g)                           0                                       0                    ⎥
         ⎢ C2h ⎥ ⎢                                                 wh        ˇh                o            wh        ˇh                                                                                                              ⎥
         ⎢ 1c ⎥ = ⎢                                                                                                                                                                                                                   ⎥
         ⎣ C2c ⎦ ⎢                                                               ˇc                                √                                                                                                                  ⎥
                  ⎢       0                    0                       e              (L+g+Lc )                  e− ˇc /    o (L+g+Lc )                    −e        ˇf / o (L+g+Lc )
                                                                                                                                                                                                     −e−     ˇf /    o (L+g+Lc )      ⎥
                  ⎢                                                               o                                                                                                                                                   ⎥
           C2f    ⎢                                                                                                                                                                                                                   ⎥
                  ⎢                                                              ˇc                                 √                                                                                                                 ⎥
                  ⎢       0                    0                       e              (L+g+Lc )                  −e− ˇc /      o (L+g+Lc )        −
                                                                                                                                                      wf        ˇf
                                                                                                                                                                     e       ˇf / o (L+g+Lc )
                                                                                                                                                                                                wf     ˇf
                                                                                                                                                                                                            e−       ˇf / o (L+g+Lc ) ⎥
                  ⎣                                                               o                                                                   wc        ˇc                              wc     ˇc                             ⎦
                               0               0                                  0                                        0                                e        ˇf / o (2L+g)
                                                                                                                                                                                                      e−     ˇf / o (2L+g)

                          ⎡ Ts − ˛ ⎤
                             ˛ −˛
                       ⎢ c 0 h⎥
                      ×⎢         ⎥
                       ⎣ ˛f − ˛c ⎦ .                                                                                                                                                                                                      (19)

                             Ts − ˛f

   As the thermal actuator reaches steady state, the heat is stored in the structure and kept constant unless the input current is removed.
The structure dissipates heat through conduction from the lower device surface to the substrate surface. The folded beam structure is
lumped by a single block element under average temperature
                        ⎡                                                  L+g                        L+g+Lc                          2L+g
         T=             ⎣˛h (L + g) + ˛c (Lc ) + ˛f (Lf ) +                      Th (x)dx +                       Tc (x)dx +                     Tf (x)dx⎦ .                                                                              (20)
            2L + g
                                                                             0                             L+g                      L+g+Lc
                                                                       M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202                                                         197

   Specifically, when ˇ > 0 the average temperature is computed from Eqs. (16), (19) and (20) as

             1                                                                     o                                                    o
      T=                     ˛h (L + g) + ˛c (Lc ) + ˛f (Lf ) +                        [C1 (ek1 − 1) − C2 (e−k1 − 1)] +                       C4 (e−k2 − e−k3 ) − C5 (e−k2 − e−k3 )
           2L + g                                                               ˇh                                                    ˇc

            +                 C5 (e−k4 − e−k5 ) − C6 (e−k6 − e−k5 )
                                                                                                                                                                                      .     (21)
                    ˇh                                ˇc                                            ˇc
      k1 =               (L + g),       k2 =                  (L + g + Lc ),       k3 =                   (L + g)
                     o                                    o                                           o

                    ˇf                                    ˇf                                          ˇf
      k4 =               (2L + g),        k5 =                    (L + g + Lc ),       k6 =                   (2L + g)
                     o                                        o                                           o

   The governing ODE which represents the decay (discharge) of temperature is
         +      d    = 0.                                                                                                                                                                  (22)
where = T − Ts , and 1/ d = [S/(RT h) + he ]/(Cp d ). d is the decay time constant. The final value needed for the thermal actuator to return
to its initial position after being deflected is approximately 5 d . S and he are the equivalent total shape factor and the heat convection
coefficient, respectively, and they are extracted from S and he given that average width w is approximated by averaging the equivalent area
of arms as
      w = [(L + g)wh + Lc wc + Lf wf ]/[2L + g]

      S = h/w(2tv /h + 1) + 1                                           .                                                                                                                  (23)

      he = whcu + 2hhcs
   The solution for the temperature decays in the structure from T into Ts is given by
      T (t) = Ts + (T − Ts )e−1/
                    ¯                        dt   .                                                                                                                                        (24)
   The approximate time ttk , needed for a thermal actuator to decay from its lumped temperature, T , to K% of Ts is given by
                      T − Ts
      ttk = log                                                                                                                                                                            (25)
                     Ts [Ä − 1]

3. Simulation

3.1. Verification

    The verification of the proposed technique for evaluating the approximate rising time is performed comparing the results with those
obtained through finite difference approximation method and the exact solution of a special case for a microbeam already reported in the
literature [7]. The steady state temperature profile of a thermal folded beam actuator was previously derived and experimentally confirmed
[12], accordingly the analytic verification was not considered on the derived cooling time of the lumped folded beam actuator. Finite element
modeling is utilized here to validate the lumped model of cooling time responses.
    The general solution of a rising temperature distribution may be obtained by Forward-Time Centered-Space method (FTCS). The finite
difference equation of Eqs. (2) and (3) can be expressed as,
                      p           t                                                    x2
      Tin+1 =                            n
                                        Ti−1 − Tin        2 + k ε x2 −                            + Ti+1 + k T∞ k ε x2
                                                                                                                               k = h, c, f.                                                (26)
                             x2                                                k˛
                                                                                 p          t
where n and i are integers that refer to the time and space mesh, and x and t are space and time grid resolutions, respectively.
    The time and space grid resolutions must be chosen to satisfy Crank–Nicholson convergence condition. The non-homogeneity in the
structure is accounted for by the change in constant properties along the spatial x-direction. For a folded beam actuator, k ˛p , k ε and k T∞
are functions of temperature and they change as the temperature profile changes across flexure, cold and hot arms.
    The continuity of flux and temperature at intermediate beam joints is automatically satisfied through space and time marching. Notice
that the derived equation in (26) is applicable to the analysis of thermal actuators with any number of serially connected microbeams {k ∈ I:
1, 2, . . .}. An actuator with non-homogeneity in material and structure configuration can be iterated by simply changing the k ˛p , k ε and
k T values while marching in the x-domain.
    A special case study simplifies the solution of the general temperature response for a folded beam actuator or V-shaped actuator with
uniform widths which is thermally equivalent to a single attached thermal microbeam of uniform width and total length L. A folded beam
actuator or a V-shaped actuator with arms of uniform width is a special case that simplifies the solution for the temperature profile.
Actuators of uniform width are equivalent to a single thermal microbeam of total length L. The exact general temperature response of such
microbeam has been investigated by Liwei [7]
                                  ∞                                                        L
                                                      2              n x       2                                         nx
      T (x, t) = e−˛p εt                e−˛p (nx/L) t sin                                      [Ts − Ts.s (x)] sin            dx                                                          (27-a)
                                                                      L        L       0
198                                                    M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202

                            Table 1
                            Material properties of folded beam actuators [11,12,16,17].

                            Parameter                                                     Si device                           Si wafer handle
                            Density, d (kg/m )                                            2330                                2330
                            Thermal conductivity (W m−1 ◦ C−1 )                           100                                 30
                            Thermal expansion, ˛ (◦ C−1 )                                 3.1 × 10−6                          3.1 × 10−6
                            Thermal capacity, Cp (J kg−1 ◦ C−1 )                          787                                 787
                            Temperature coefficient, (◦ C−1 )                              1.25 × 10−3                         1.25 × 10−3
                            Electrical resistivity, o ( m)                                1.5 × 10−4                          2.5 × 10−2

                      Table 2
                      Dimensions of lineshape beam actuators.

                      Parameter                                                           SOI, folded ( m)                          SOI, single ( m)

                      Air gap, tv                                                           2                                         2
                      Thickness of substrate layer, ts                                    300                                       300
                      Width of the hot arm, wh                                             10                                        10
                      Width of cold arm, wc                                                52                                        10
                      Width of cold flexure arm, wf                                         10                                        10
                      Length of hot arm, Lh                                               370                                       602
                      Length of cold arm, Lc                                              271                                         0
                      Length of flexure arm, Lf                                             99                                         0
                      Structure thickness, h                                              100                                       100
                      Gap distance between cold and hot arms, g                            20                                         –

where Ts.s (x) is the steady state temperature distribution given by
                 T∞ − (T∞ − Ts ) cosh[ ε(x − L/2)]
     Ts,s(x) =                  √                  .                                                                                                                (27-b)
                           cosh[ ε(L/2)]
    The exact transient time constant is evaluated from the transient solution in Eq. (27-a) and given by
 n = [(S o /(hRT ) − J2 o )/( d Cp ) + (nx/L)2 /( d Cp )]−1 . A set of finite and sequentially connected microbeams of uniform width are equivalent
to a folded beam actuator of uniform width wf = wc = wh .
    Thus, the approximate general temperature response of a lineshape microbeam can be obtained from the procedures described in Eqs.

                                                     a d   −     b c
       T (x, t) = Ts + (x2 −
       ˜                            c x)
                                                +                      e−(   d / c )t                                                                                 (28)
                                            d              c d

where a = (Ts − T∞ − Ti ), b = ˛p ε(Ts − T∞ ), c = − 5L2 /6, and d = ˛p (2 − 5εL2 /6).
   The parameters used in the simulations including geometry, material properties and heat transfer coefficients are shown in Tables 1 and 2.
The natural convection coefficient of a microbeam is obtained from [15], and it is a function of geometric size, surface roughness, device
surface temperature and convective face orientation. Initial temperatures of substrate lower face and contact pads are defined at Ti = 20 ◦ C
and Ts ∼ 20 ◦ C.
   The general temperature response of a single lineshape SOI microactuator evaluated using FDA simulation is presented in Fig. 2. The
fastest temperature response occurs next to the pad boundary where the temperature is fixed. However, the time required to reach steady
state increases as the location moves towards the microbeam center which has the slowest time response. The slowest time required to
settle within 2% of the final value is depicted from the temperature response at the midsection point as shown in Fig. 2b. Also, the exact
slowest rising time response can be obtained at (n = 1, x = L/2). The approximate rising time constant can either be calculated directly from
Eq. (28) or by choosing the largest time constant evaluated by Eq. (13).

Fig. 2. A single SOI microbeam simulated using FDA for input of (35 mA or 5 V), and mesh grid of ( t = 0.5 s,   t = 10 s): (a) temperature response; T(x,t), and (b) slowest
temperature response at Tmax (t) = T(x = L/2;t).
                                                       M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202                                                 199

                              Table 3
                              Ramping time of a beam with uniform widths using 4/( t ).

                              Method                                            Heating time ( s)                        Error (Liwei base)

                              Liwei result [7]                                  423.6
                              FDA                                               437.2                                    3.2%
                              Approximate—proposed                              439.3                                    3.7%

Fig. 3. Attached folded beam actuated at 4.4 V square input (a) ramping up to steady state temperature distribution using FEM and (b) comparison of the average transient
cooling response using approximate (proposed) and FEM methods [14].

    The proposed approach is verified by employing it to calculate the heating time for a single lineshape microbeam of uniform width as
presented by Liwei [7] and through FDA. The results of the heating time constants evaluated by the three methods are presented in Table 3
and indicate that the excellent agreement between them verifies the proposed approach for evaluating the heating time.
    The SOI folded beam actuator with initial body temperature is tested for 4.4 V step voltage input and then switched off. Results in Fig. 3
shows that the lumped method corresponds well with the FEM. in the next section, the method will be used to compare the thermal cycle
of folded beam actuators fabricated on different material.

3.2. Examples on thermal cycle of folded beam actuator

    The thermal cycle responses for folded beam actuator is evaluated for actuators fabricated using SOI using the geometries and material
properties (Tables 1 and 2). The performance of a device largely depends on the material and fabrication processes. A folded beam actuator
is fabricated on an SOI wafer composed of four sandwiched layers; 300 m thick silicon handle wafer, 2 m thick SiO2 with air gap, 100 m
thick silicon device structure, and metal coated pads.
    A commercial 3D MEMS dynamic profiler, Wyko/Veeco NT1100 DMEMS profiler [18], was used to measure the dynamical response of a
SOI folded beam actuator. The lateral deflection of a folded beam tip has been recorded by an interferometric position sensor that probes
the full output cycle with ı = 5◦ increment.
    The frequency response transfer function, G(s), of a thermally stable actuator is obtained for first order model using sinusoidal frequency
input at 12 V amplitude. A fine search for the locations of poles and zeros are subjected to minimize the magnitude of a norm between
experimental magnitudes and the magnitude of the iterated locations

                −5.833 × 10−5 s + 0.0589
       G(s) =                                     15 ≤ ω ≤ 265                                                                                                     (29)
                  0.000505s + 0.3043
   For frequencies higher than 215 Hz, the structural dynamic response of square input signal input is obtained for 50% duty cycle and 24 V
amplitude. Fig. 4 displays response for frequencies ranges from 215 to 1015 Hz. At ω ≥ 265 Hz, the ramping region interferes with retracting
region without reaching a steady state. The main reason is that the excitation power frequency generates an effective thermal cycle that is
slower than mechanical cycle, i.e., the actuator is excited repeatedly before its structure completely responds.
   The input signal that stimulate no dwell time is obtained from Eq. (29) and plot in Fig. 5. The significance of this region is that both
structural and thermal cycles are close, but with structural response always faster as depicted from Table 4, where the computed thermal
time of heating and cooling responses are based on the settling to the 98% of S.S temperature and 1.02% of initial temperature Ts .

Table 4
Thermal and structural time comparison for SOI folded beam actuators at 24 V.

SOI folded beam actuator                         Thermal                                                                     Structural

                                                 Heat 4/( t )max (s)                        Cool ttk (s)                     Ramp (s)                         Retract (s)

Approximated                                     0.00844                                    0.02280                          –                                –
Experimental                                     –                                          –                                0.00878                          0.01042
200                                               M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202

                                     Fig. 4. Dynamic response of the SOI thermal folded actuator attached on 1.0 cm2 die.

                       Fig. 5. Mechanical response of the SOI folded beam actuator simulated for square input of 50% duty cycle and 24 V.

   The analytical thermal cycle analysis of the SOI folded beam suggested that the actuator could operate safely a full cycle of charging
and discharging without monitoring any heating on the MEMS device layer at the voltage input frequency ω ≤ 13.6 Hz. At higher frequency
input, the average temperature stored on the thermal actuator increases due to the successive overlap. However, the accumulation of heat
did not result in an observable mechanical offset as observed in the range of 15 ≤ ω ≤ 265 Hz. This can be explained by the dynamic coupling
effects in thermo-elastic actuators. Purely elastic MEMS devices governed by high order transfer function have frequency response that
undergoes attenuation on the deflection. Meanwhile, additive heat storage on such device causes additional deflection which is appreciated
when the total MEMS package overheat.
   Moreover, it is observed that the folded beam actuator tends to saturate at frequency higher than 1015 Hz, i.e., the increase in frequency
results in a more static deflection and a less ripple deflection as noticed in Fig. 4. This trend continues until structural dynamics is no longer
observed and the folded beam actuator only generates static deflection at ω → ∞.

4. Discussion

   The proposed approximate method for computing the heating time constant of three serially connected homogeneous beams of different
widths resulted in three fundamental time constants that correspond to the approximate transient effect of each arm of the folded beam
actuator. For example, the slowest heating time of the SOI folded beam actuator based on the approximate method was calculated to take
16 ms for a 4.4 V input. Although the derived method was in excellent agreement with FDA in the case of uniform widths, the FDA shows a
small difference in the simulated heating time of folded beam, which is 13 ms. Also, a small difference in the results also obtained by the
proposed method was observed in the steady state temperature profiles of folded beam when the folded beam started to discharge heat
at an average temperature of 82 ◦ C for an exact based solution. Meanwhile the FEM simulation revealed average temperature of 100 ◦ C.
Thus, the cooling time of the lumped model took a little longer than the simulated FEM with heating and cooling times of 15.6 and ∼13 ms,
                                                          M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202                                                    201

   The importance of mechanical response experiment is to reveal the practical limitation of the method under the model’s assumptions. The
approximate models have close physical insight to real part for low temperatures, where the material properties conditions are temperature
independent, the boundary conditions has fixed temperatures and the radiation is insignificant.
   The transient analysis using numerical approaches (FDA and FEM) was primitively slow compared to proposed method, often unstable
and inaccurate attributed to the large model size resulting from the fine spatial and time discretization, solution round off errors and ill
conditioned matrices, limitations of algorithm’s stability and convergence. The proposed procedure provides an analytical methodology
that can easily and quickly compute the thermal cycle for folded beam actuators. It is observed that the thermal cycle of E-T folded beam
actuator explain the structural response and also indicate its stability.
   Additional simulations revealed that an attached actuator conducts more heat to the substrate than a suspended actuator, and thus it
can generate faster responses. Fast thermal and structural cycles can generally be achieved by scaling down the folded beam dimension
providing a thermally conductive substrate at low temperature, reducing the air gap between the device layer and the substrate, and last
but not least, reducing input power by either reducing input voltage amplitude or reducing the effective input duty cycle.

5. Conclusions

   A novel approximate method for rapidly evaluating the thermal cycle of serially connected microbeams of varying width acting as E-T
actuators was presented. The thermal cycle results obtained by the proposed method were in good agreement with results obtained from
time consuming numerical approaches such as FDA and FEM. The discharging time of E-T actuator is generally higher than the charging
time at high input voltage. The computational performance of the obtained analytical expressions could be employed in future to reverse
engineer E-T MEMS devices based on their scale, boundary conditions and physical properties in order to meet specified performance


   This research was partially supported from Grant # N00014-05-1-0587 from the Office of Naval Research. The authors would like to
thank Prof. Dan Popa, and Dr. Ping Zhang and Dr. Nitin Uppal for their assistance with fabrication and insightful discussions on E-T.


 [1]   R. Hickey, D. Sameoto, T. Hubbard, et al., Time and frequency response of two-arm micromachined thermal actuators, J. Micromech. Microeng. 13 (2003) 40–46.
 [2]   J.T. Butler, V.M. Bright, W.D. Cowan, Average power control and positioning of polysilicon thermal actuators, Sensors Actuators 72 (1999) 88–97.
 [3]   M. Pai, N.C. Tien, Low voltage electrothermal vibromotor for silicon optical, bench applications, Sensors Actuators 83 (2000) 237–243.
 [4]   J.H. Comtois, V.M. Bright, Application for surface micromachined polysilicon thermal actuator, Sensors Actuators A 58 (1997) 19–25.
 [5]   C.-C. Lee, W. Hsu, Optimization of an electro-thermally and laterally driven microactuator, Microsystem Technol. 9 (2003) 331–334.
 [6]   D.O. Popa, B.H. Kang, J.T. Wen, H.E. Stephanou, G. Skidmore, A. Geisberger, Dynamic modeling and input shaping of thermal bimorph MEMS actuators, ICRA (2003)
 [7]   L. Lin, M. Chiao, Electrothermal responses of lineshape microstructures, Sensors Actuators A 55 (1996) 35–41.
 [8]   E.T. Enikov, K. Lazarov, PCB- integrated metallic thermal micro-actuators, Sensors Actuators 105 (2003) 76–82.
 [9]   C.S. Pany, W. Hsu, An electro-thermally and laterally driven polysilicon microactuator, J. Micromech. Microeng. 7 (1997) 7–13.
[10]   S. Kiwan, M. A1-Nimr, M. A1-Sharo’a, Trial solution methods to solve the hyperbolic heat conduction equation, Int. Comm. Heat Mass Transfer 27 (2000) 865–876.
[11]   C.D. Lott, T.W. McLain, J.N. Harb, L.L. Howell, Modeling the thermal behavior of a surface-micromachined linear-displacement thermomechanical microactuator, Sensors
       Actuators A 101 (2002) 239.
[12]   Q.-A. Huang, N.K.S. Lee, Analysis and design of polysilicon thermal flexure actuator, J. Micromech. Microeng. 9 (1999) 64–70.
[13]   M. Mayyas, P.S. Shiakolas, A study on the thermal behavior of electrothermal microactuators due to various voltage inputs, Proceedings of IMECE 2006, Paper No.
       IMECE2006-15321, Chicago IL, November 2006.
[14]   ANSYS Inc. Southpointe, 275 Technology Drive. Canonsburg, PA, 15317,
[15]   N.D. Mankame, G.K. Ananthasuresh, Comprehensive thermal modelling and characterization of an electro-thermal-compliant microactuator, J. Micromech. Microeng.
       11 (2001) 452–462.
[16]   A. Atre, Analysis of out-of-plane thermal microactuators, J. Micromech. Microeng. 16 (2006) 205–213.
[17]   Material is available at:
[18]   3D MEMS profiler, “WYKO NT1100”
[19]   MATLAB 7.2, The MathWorks Inc Nov. 2005: available at:
[20]   R.G. Li, Q.A. Huang, W.H. Li, A nodal analysis method for simulating the behavior of electrothermal microactuators, Microsyst. Technol. 14 (2008) 119–129.
[21]   R. Hickey, D. Sameoto, T. Hubbard, M. Kujath, Time and frequency response of two-arm micromachined thermal actuators, J. Micromech. Microeng. 13 (2003) 40–46.
[22]   Q.-A. Huang, N.K.S. Lee, Analytical modeling and optimization for a laterally-driven polysilicon thermal actuator, Microsyst. Technol. 5 (1999) 133–137.
[23]   B. Borovic, F.L. Lewis, D. Agonafer, E.S. Kolerser, M.M. Hossain, D.O. Popa, Method for determining a dynamic state-space model for control of the thermal MEMS devices,
       J. Microelectromech. Syst. 14 (2005) 961–970.


Mohammad A. Mayyas received PhD in mechanical engineering and MSME all from University of Texas at Arlington, 2003–2007, BSME from Jordan University of Science
and Technology, 2001. He is currently a research scientist with the Automation & Robotics Research Institute/UT-Arlington. His research interests include micro-surfaces
reconstruction, synthesis and analysis of microsystems, microassembly and micro-actuations techniques, mechatronics and micro-robotics.

Panos S. Shiakolas received PhD from University of Texas at Arlington in 1992, MSME from University of Texas at Austin in 1988 and BSME from University of Texas at Austin
in 1986. He joined the Mechanical and Aerospace Engineering Department of the University of Texas at Arlington as an assistant professor in Fall 1996. Before joining UTA, he
has performed extensive research in the areas of rapidly configured robotic and manufacturing systems. Since joining UTA, he has performed extensive research in the areas
of applying robotics-based automation in aerospace manufacturing applications. He is currently focusing his research efforts in rapid microfabrication, micromechatronics,
microrobotics, microassembly and MEMS.

Woo Ho Lee holds a BS degree in mechanical design and production engineering from Hanyang University, Seoul Korea and a MS degree in production engineering from the
Korea Advanced Institute of Science and Technology in Taejon, Korea. He worked at the Living System Research Laboratory of Gold Start Electronics Co. from 1990 to 1993,
and at the Robotics and Fluid Power Control Laboratory of the Korea Institute of Science and Technology from 1993 to 1995. In 2001, he received a PhD degree in mechanical
engineering from Rensselaer Polytechnic Institute, focusing on reconfigurable robotic systems. Upon graduation, he joined the Center for Automation Technologies as a
research associate. He is presently a faculty associate with the Automation & Robotics Research Institute at The University of Texas at Arlington. He has been involved in the
development of dynamic modeling and tolerance analysis algorithms for microrobotic assembly, the design of MEMS sensors and actuators, and the packaging of biosensors.
202                                                      M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202

Harry Stephanou is a professor of electrical engineering and the Director of the Automation & Robotics Research Institute at The University of Texas at Arlington. He is also the
Founder of BMC and Chairman of its Board of Directors. He has received a PhD in electrical engineering from Purdue University in 1976, focusing on robotics and intelligent
systems. From 1976 to 1985, he was with the Exxon Production Research Company in Houston, TX. In 1980, he was appointed Head of the Systems Research Section in the
Long Range Research Division. From 1985 to 1990, he was on the faculty of the School of Information Technology and Engineering at George Mason University, in Fairfax, VA.
From 1987 to 1988, he also served as part-time Program Director for Robotics and Machine Intelligence at the National Science Foundation. Between 1990 and 2004, he was
with Rensselaer Polytechnic Institute in Troy, NY as Director of the Center for Automation Technologies. His current research interests are in modular microrobotic systems
and distributed micromachines. He has held several leadership positions including Vice President, Conference Program Chairman and Associate Editor at the IEEE Society for
Robotics and Automation.

Shared By: