Thermal cycle modeling of electrothermal microactuators

Sensors and Actuators A 152 (2009) 192–202

Contents lists available at ScienceDirect

Sensors and Actuators A: Physical journal homepage: www.elsevier.com/locate/sna

Thermal cycle modeling of electrothermal microactuators Mohammad Mayyas a,∗ , Panos S. Shiakolas b , Woo Ho Lee a , Harry Stephanou a a b

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

Article history: Received 25 May 2008 Received in revised form 4 February 2009 Accepted 5 March 2009 Available online 2 April 2009 Keywords: Thermal cycle MEMS Actuators Electrothermal

a b s t r a c t A methodology for computing the transient heating (charging) and cooling (discharging) time for a serially connected electro thermo (E-T) microbeams actuator were analytically approximated from fundamental characteristics including input voltage, material and geometry. Specifically, this paper presents closed form solutions for folded and bend beam MEMS actuators. The E-T heating time was obtained from the 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 examples and validated by Finite Difference Approximation (FDA), and Finite Element Modeling (FEM). 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: [email protected] (M. Mayyas). 0924-4247/$ – see front matter. Published by Elsevier B.V. doi:10.1016/j.sna.2009.03.015

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)

 I 2

∂T d Cp dx = ∂t

wh

o







S d2 T 1 + (T − Ts ) dx − − o 2 dx + h dx

T − T  s

RT



dx + (whcu + 2hhcs )(T − Tamb )dx

(1)

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 considered. 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, 1 ∂ ∂2  = + ε. ˛p ∂t ∂x2

(2)

where ˛p =

 = (T − T∞ )

he = whcu + 2hhcs



T∞ =

Ts +

J 2 o ε o



o

d Cp ε=

S hRT

o

+

he o

−

J 2 o 

.

o

(3)

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 : ⎪ ⎪ ⎪ ⎪ ⎪ ⎨x = L +g :

Ts − T∞h s T∞c − T∞h ; Hh − Hc = s

Hh =

T∞f ⎪ ⎪ x = L + g + Lc : Hh − Hf (L + g + Lc , s) = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Ts − T∞f ⎩ x = 2L + g : Hf = s

dHh dHc − wc =0 dx dx . − T∞c dHf dHc ; wc − wf =0 s dx dx

wh

(4)

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) − dx2



s + ˛p εk ˛p



Hk (x, s) +

Ti = 0 . . . k = {h, c, f } ˛p

(5)

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 (6) ˜ k (x, s, T∞k , ˛, ε) = a0k + a1k x + a2k x2 . . . k = {h, c, f } H

(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

Rf (x, s) = Ti /˛p + a0f 1f + a1f 2f + a2f 3f

2L + g ≥ x ≥ L + g + Lc

(7)

where the coefficients of undetermined residuals are



1k =

−

 2k =

s + ˛p εk ˛p



s + ˛p εk − x ˛p



⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

⎪ ⎪ ⎪  ⎪ ⎪ s + ˛p εk 2 ⎪ ⎪ ⎪ x 3k = 2 − ⎭ ˛

· · ·k = {h, c, f }

(8)

p

The next step is to find the best set of parameters of the basis function approximation, which minimize the error between the approximated 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 coefficients 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)



Rh dx = 0 0



(2L+g)

(L+g+Lc )

Rc dx = 0 (L+g)

Rf dx = 0. (L+g+Lc )

(9)

M. Mayyas et al. / Sensors and Actuators A 152 (2009) 192–202

195

˜ into both the auxiliary conditions in Eq. (4) and The undetermined coefficients are evaluated by substituting the approximating H optimization criteria in Eq. (7)

⎡

⎤−1

Aa . [a] = ⎣ .. ⎦ Ab

[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

 b=

Ts − T∞h s

T∞c − T∞h s

0

T∞f − T∞c s

Ts − T∞f

0

−

s

Ti (L + g) ˛p

−

Ti (Lc ) ˛p

−

Ti (L − Lc ) ˛p

T

The boundary conditions and flux continuity between the microbeams yield matrix Aa as

⎡

1 0 ⎢ 1 (L + g) ⎢ 1 ⎢0 Aa = ⎢ 0 ⎢0 ⎣0 0 0 0

0 0 0 (L + g)2 −1 −(L + g) 2(L + g) 0 −(wc /wh ) 0 1 (L + g + Lc ) 0 0 1 0 0 0

0 −(L + g) −2(L + g)(wc /wh ) (L + g + Lc )2 2(L + g + Lc ) 0

0 0 0 −1 0 1

0 0 0 −(L + g + Lc ) −(wf /wc ) (2L + g)

(11)

⎤

0 ⎥ 0 ⎥ 0 ⎥ ⎥. −(L + g + Lc )2 ⎥ −2(L + g + Lc )(wf /wc ) ⎦ (2L + g)2

(12-a)

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

⎡

Ab =

⎢− ⎢ ⎢ ⎣

(L + g) (s + ˛p εh ) ˛p 0

−

(L + g)2 (s + ˛p εh ) 2˛p 0

0

2L + 2g −

(L + g)3 (s + ˛p εh ) 3˛p 0

0

0

0 Lc − (s + ˛p εc ) ˛p 0

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

0

0

0

0

0

0

−

(L − Lc ) (s + ˛p εf ) ˛p

−

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

(2L + g)2 − (L + g + Lc )2 (s + ˛p εf ) 2˛p

2(L − Lc ) −

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

⎤ ⎥ ⎥ ⎥ ⎦

(12-b)

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)] . ⎪ ⎪ ⎩ A (s) ⎭ b

(14)

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 T˜k (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

 I = min



Ii =

kv ks wi 0 (2tv + h + wi ) + (kv ts + ks tv )(wi hcui + 2hhcsi ) h0 (kv ts + ks tv )

 k = {h, c, f }

(17)

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

⎡ ⎤ ⎥ wc dTc (L + g + Lc ) = wf dTf (L + g + Lc ) ⎥ ⎣ dx dx ⎦ ⎥, ⎦ w dTh (L + g) = w dTc (L + g)

⎢ T (2L + g) = Ts ⎢ f ⎢ ⎣ Th (L + g) = Tc (L + g)

h

c

dx

Tc (L + g + Lc ) = Tf (L + g + Lc )

(18)

dx

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

⎡

1

1

⎢ √ˇ / o (L+g) ⎢e h ⎢ ⎡C ⎤ ⎢ √ 1h ⎢ e ˇh / o (L+g) ⎢ C ⎢ ⎢ C2h ⎥ ⎢ 1c ⎥ = ⎢ ⎣ C2c ⎦ ⎢ ⎢ 0 C1f ⎢ C2f ⎢ ⎢ ⎢ 0 ⎣ 0

√ e− ˇh / √ −e− ˇh /

0



−

wc wh

 

 ˇc e ˇh

e

e

ˇc

ˇc

(L+g)

o

wc wh

0

o (L+g)

ˇc −√ˇc / e ˇh

o (L+g)

−e



0

−

ˇf / o (L+g+Lc )

wf

ˇf

wc



ˇc

e

 e

⎤−1

⎥ ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥  ⎥ ˇf / o (L+g+Lc ) ⎥ −e− ⎥  ⎥  ⎥ wf ˇf − ˇf / o (L+g+Lc ) ⎥ e wc ⎦ ˇc  0

0



o (L+g+Lc )

0

0

o (L+g+Lc )

√ −e− ˇc /

(L+g+Lc )

o

0



√ e− ˇc /

(L+g+Lc )

o

ˇc

0

√ −e− ˇc /

(L+g)

o

 0

0

ˇc

−e

o (L+g)

o (L+g)

0

ˇf / o (L+g+Lc )

ˇf / o (2L+g)

e−

ˇf / o (2L+g)

⎡ Ts − ˛ ⎤ h

˛ −˛

⎢ c 0 h⎥ ⎥ ×⎢ ⎣ ˛f − ˛c ⎦ .

(19)

0 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

⎡ 1 T¯ = 2L + g

⎣˛h (L + g) + ˛c (Lc ) + ˛f (Lf ) +





L+g

Th (x)dx + 0



2L+g

L+g+Lc

Tc (x)dx + L+g

L+g+Lc

⎤ Tf (x)dx⎦ .

(20)

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 T¯ = 2L + g

 o

+

ˇf

 k1 =

ˇh o

 k4 =



˛h (L + g) + ˛c (Lc ) + ˛f (Lf ) +

o

ˇh

[C1 (ek1 − 1) − C2 (e−k1 − 1)] +

o



ˇc



C4 (e−k2 − e−k3 ) − C5 (e−k2 − e−k3 )







C5 (e−k4 − e−k5 ) − C6 (e−k6 − e−k5 )

 (L + g),

k2 =

ˇc o



ˇf



o

(2L + g),

k5 =

 (L + g + Lc ),

ˇf o

k3 =

. ˇc o

 (L + g + Lc ),

k6 =

(21)

(L + g)

ˇf o

(2L + g)

The governing ODE which represents the decay (discharge) of temperature is d

+ d = 0. dt

(22)

¯ where = T − Ts , and 1/ d = [S/(R T 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 ¯ is approximated by averaging the equivalent area coefficient, respectively, and they are extracted from S and he given that average width w of arms as ¯ = [(L + g)wh + Lc wc + Lf wf ]/[2L + g] w S = h/w(2tv /h + 1) + 1

.

(23)

¯ cu + 2hhcs he = wh The solution for the temperature decays in the structure from T¯ into Ts is given by T (t) = Ts + (T¯ − Ts )e−1/ d t .

(24)

The approximate time ttk , needed for a thermal actuator to decay from its lumped temperature, T¯ , to K% of Ts is given by



ttk = log

T¯ − Ts Ts [ − 1]

d

(25)

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, Tin+1 =

k

˛p t



x2



n Ti−1 − Tin

2 + k ε x2 −



x2 k ˛ t p



n + Ti+1 + k T∞ k ε x2

k = h, c, f.

(26)

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] T (x, t) = e

−˛p εt

∞  

e

n=1

−˛p (nx/L)2 t

sin

 n x   2  L

L

0

L

[Ts − Ts.s (x)] sin

 nx  L



dx

(27-a)

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 3

Density, d (kg/m ) Thermal conductivity (W m−1 ◦ C−1 ) Thermal expansion, ˛ (◦ C−1 ) Thermal capacity, Cp (J kg−1 ◦ C−1 ) Temperature coefficient,  (◦ C−1 ) Electrical resistivity, o ( m)

Si device

Si wafer handle

2330 100 3.1 × 10−6 787 1.25 × 10−3 1.5 × 10−4

2330 30 3.1 × 10−6 787 1.25 × 10−3 2.5 × 10−2

Table 2 Dimensions of lineshape beam actuators. Parameter

SOI, folded (␮m)

SOI, single (␮m)

Air gap, tv Thickness of substrate layer, ts Width of the hot arm, wh Width of cold arm, wc Width of cold flexure arm, wf Length of hot arm, Lh Length of cold arm, Lc Length of flexure arm, Lf Structure thickness, h Gap distance between cold and hot arms, g

2 300 10 52 10 370 271 99 100 20

2 300 10 10 10 602 0 0 100 –

where Ts.s (x) is the steady state temperature distribution given by √ T∞ − (T∞ − Ts ) cosh[ ε(x − L/2)] . Ts,s(x) = √ cosh[ ε(L/2)]

(27-b)

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. (1)–(16) T˜ (x, t) = Ts + (x2 −



c x)



b

d

+

  −    a d b c c d

e−(d /c )t



(28)

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] FDA Approximate—proposed

423.6 437.2 439.3

3.2% 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 G(s) =

−5.833 × 10−5 s + 0.0589 0.000505s + 0.3043

15 ≤ ω ≤ 265

(29)

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

Approximated Experimental

Thermal

Structural

Heat 4/( t )max (s)

Cool ttk (s)

Ramp (s)

Retract (s)

0.00844 –

0.02280 –

– 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, respectively.

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 requirements. Acknowledgments 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. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

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. J.T. Butler, V.M. Bright, W.D. Cowan, Average power control and positioning of polysilicon thermal actuators, Sensors Actuators 72 (1999) 88–97. M. Pai, N.C. Tien, Low voltage electrothermal vibromotor for silicon optical, bench applications, Sensors Actuators 83 (2000) 237–243. J.H. Comtois, V.M. Bright, Application for surface micromachined polysilicon thermal actuator, Sensors Actuators A 58 (1997) 19–25. C.-C. Lee, W. Hsu, Optimization of an electro-thermally and laterally driven microactuator, Microsystem Technol. 9 (2003) 331–334. 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) 1470–1475. L. Lin, M. Chiao, Electrothermal responses of lineshape microstructures, Sensors Actuators A 55 (1996) 35–41. E.T. Enikov, K. Lazarov, PCB- integrated metallic thermal micro-actuators, Sensors Actuators 105 (2003) 76–82. C.S. Pany, W. Hsu, An electro-thermally and laterally driven polysilicon microactuator, J. Micromech. Microeng. 7 (1997) 7–13. 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. 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. Q.-A. Huang, N.K.S. Lee, Analysis and design of polysilicon thermal flexure actuator, J. Micromech. Microeng. 9 (1999) 64–70. 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. ANSYS Inc. Southpointe, 275 Technology Drive. Canonsburg, PA, 15317, http://www.ansys.com/products/cfx. N.D. Mankame, G.K. Ananthasuresh, Comprehensive thermal modelling and characterization of an electro-thermal-compliant microactuator, J. Micromech. Microeng. 11 (2001) 452–462. A. Atre, Analysis of out-of-plane thermal microactuators, J. Micromech. Microeng. 16 (2006) 205–213. Material is available at: www.matweb.com. 3D MEMS profiler, “WYKO NT1100” www.veeco.com. MATLAB 7.2, The MathWorks Inc Nov. 2005: available at: http://www.mathworks.com/. 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. 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. Q.-A. Huang, N.K.S. Lee, Analytical modeling and optimization for a laterally-driven polysilicon thermal actuator, Microsyst. Technol. 5 (1999) 133–137. 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.

Biographies 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.