Development of an accurate heat conduction

Microsyst Technol (2015) 21:345–353 DOI 10.1007/s00542-014-2079-x

TECHNICAL PAPER

Development of an accurate heat conduction model for micromachined convective accelerometers Brahim Mezghani · F. Tounsi · M. Masmoudi

Received: 22 October 2013 / Accepted: 8 January 2014 / Published online: 29 January 2014 © Springer-Verlag Berlin Heidelberg 2014

Abstract This paper presents detailed development of a heat conduction model adapted for micromachined convective accelerometers. Fitting expressions obtained from FEM data are used in derived spherical model expressions to come up with an accurate analytical model of heat conduction main parameters: common mode temperature and heat transfer coefficient. Two variables have been used in FEM analysis: applied heater temperature and micromachined cavity depth. The latter parameter has a large impact on the overall conductive behavior of thermal accelerometers since it fixes the volume where produced heat bubble can expand. In addition, this depth is one of the two main reasons behind the produced isotherms deformation. Isotherm bending is also due to the high aspect ratio in heater width and height imposed by the technology. Since the produced hot bubble form is closely related to sensor design and heater temperature, then spheres deformation modeling has been used to derive conduction model equations. Two distinct equivalent radius modeling are studied and are used to express conduction behavior in analytical expressions. For comparison, Hardee’s conduction solution is given and it has been found that our solution gives a more accurate reading of common mode temperature. Therefore, Hardee’s solution has to be revised if it is to be used for convective accelerometers. It is also shown that derived expressions are still valid for various sensor designs and that conductive behavior of thermal accelerometers can be predicted in an early stage and for all possible design geometries and biasing temperatures.

B. Mezghani (*) · F. Tounsi · M. Masmoudi EMC Research Group, Department of Electrical Engineering, National Engineering School of Sfax (ENIS), University of Sfax, Route Soukra Km3.5, BP 1173, 3038 Sfax, Tunisia e-mail: [email protected]

1 Introduction In recent years, the interest in convective or thermal micromachined accelerometers has been largely increasing. Operation of this MEMS accelerometer is based on the principle of heated gas molecules convection inside a sealed cavity. Highly accurate dual temperature sensors detect the change in temperature profile resulting from very small changes in acceleration (or inclination). Simplicity of sensor structure and its possible monolithic integration with on-chip processing electronics, offering low production cost, are only few benefits which contributed to their widespread use in consumer electronics. Moreover, the absence of seismic mass added to their compact size makes them extremely robust and reliable (Garraud et al. 2011a; Bahari and Leung 2011; Kozlov and Agafonov 2011; Lin 2012; Hua et al. 2011; Chen and Shen 2008). MEMSIC is one of the first companies to integrate this sensor type along with signal processing circuitry on a standard monolithic CMOS wafer (http://www.memsic.com/index.cfm). Potential applications include automotive industry, communication and video projection systems and many varieties of other applications. These accelerometers have been widely studied in literature, but their optimization still needs to be investigated because of the complex effects of design parameters on both conduction and convection behaviors of the sensor (Garraud et al. 2011b; Usung et al. 2011; Rekik et al. 2011; Mezghani et al. 2013a; Hardee 1966; Mack and Hardee 1968; Luo et al. 2002; Courteaud et al. 2008). In addition, requirements for monolithic process compatibility impose strict constraint on sensor design parameters as input power, cavity size and package geometry (Nguyen et al. 2013; Tsang et al. 2008). Among many thermal convectionbased accelerometers, standard CMOS design ones are of

13

346

particular interest. This includes single, dual and three axis convective accelerometers. Unlike single axis thermal accelerometers, where cylindrical model is used, dual axis counterparts are generally more adapted for spherical model. However, thermal accelerometers specific design parameters can be quite different from one sensor to another. Therefore, to be able to properly use spherical model expressions, FEM simulations could be used to fit these parameters into different modeling equations. The result will be complete general analytical model expressions, which are function of different design parameters and heating temperature. Heat conduction phenomenon is the basis behind any heat flow since it sets the initial temperature distribution in the cavity. In this paper, the conductive behavior of micromachined convective accelerometers is studied through the radius deformation modeling of the obtained conduction heat bubble. Based on a newly developed 3D geometrical model for thermal accelerometers, FEM simulations using the Ansys software are done on the sensor. Obtained specific values are used in detailed analytical modeling expressions of heater heat transfer coefficient and common mode temperature. Results obtained from this new modeling are compared to those given by Hardee’s model, detailed in (Hardee 1966; Mack and Hardee 1968), and FEM simulation data. Our new modeling is found to be more accurate since its governing equations are function of both biasing temperature and sensor design geometry.

2 Convective accelerometer presentation Single and dual axis MEMS thermal accelerometers have similar operating mode concerning common mode and convection phenomena. Operating principle is thoroughly explained in (Chen and Shen 2008) and (Garraud et al. 2011b; Usung et al. 2011). Different dual axis convective accelerometer configurations were fabricated then characterized and published: • Garraud et al. (2011b) presented a dual axis CMOS micromachined thermal accelerometer with a square shaped hot plate containing a metallic meander forming the heater. Temperature detectors are placed orthogonally in the heater plane. To release the structures, an additional micromachining (anisotropic etching) step is added to the CMOS chip. Detectors are made with a set of thermocouples, which are based on the Seebeck effect and provide a potential difference proportional to the temperature difference between their hot and cold junctions. This gives a better performance of the bias shift with respect to any room temperature variation.

13

Microsyst Technol (2015) 21:345–353

• Usung et al. (2011) introduced a dual-axis micromachined convective accelerometer with a diamond-shaped heater geometry, which increases sensor sensitivity. Since temperature gradient, proportional to sensitivity, is larger when using a diamond-shaped heater, therefore, this heater form is better than the square shape. To achieve a greater concentration of generated heat, a new implementation of heater complex pattern, which is thinner and longer than conventional heaters, has been used. These dual axis accelerometers can also be used as single axis and three axis ones (Nguyen et al. 2013). Many more research groups fabricated and tested other convective accelerometer designs. Along with characterization results, they basically show numerical simulations which were developed to fit only on their specific design. There are no studies which focus on the use of analytical modeling in predicting sensor performance as a function of key geometric parameters and heating temperature. We use a newly developed 3D geometry model in FEM simulations to adapt the available spherical model equations on a general accelerometer design. This 3D model has been originally presented in (Mezghani et al. 2011) for a single axis thermal accelerometer and is practically the same for dual and three axis ones. Moreover, the same technique is used to evaluate temperature values along sensitive axes. The script describing our new 3D model and used in FEM simulations has been validated through an experimental result comparison with a device realized by (Garraud et al. 2011b). In this comparison, detailed in (Mezghani et al. 2013b, c), the design parameters used in simulations are identical to those presented in (Garraud et al. 2011b). Here, we will specifically be interested in the conductive behavior modeling. In the following sections, sensor parameters and their nominal values, used in both FEM simulations and equations evaluation, are listed in Table 1. Values chosen for rd, h2, h1 and d2 are supposed to give maximum sensitivity readings for an ro value of 300 μm. Table 1 List of accelerometer parameters and nominal values Symbol

Description

Value (μm)

rd

Distance from heater center to detector

ri

Heater half width

30

ro

Cavity half width

300

h1

Bottom cavity depth

h2

Top cover height

E

Suspended structures thickness

d2

Distance from bottom cavity edge to cover

W

Top cover total width = 2 (ro + d2)

120

300 1,000 7 600 1,800


347

Top cover Heater

h2 ro

rd

h1

2r i

To

Ti

ro

ri

Detector d2

Cavity (a)

(b)

Si substrate Fig. 1 Simplified accelerometer cross-sectional view used in simulations and showing main modeling parameters

Concerning ri, the specified value offers maximum efficiency (Mezghani et al. 2013b, c). These key geometric parameters which will be used in both numerical and analytical modeling are described in the sensor simplified cross-sectional view shown in Fig. 1. The analytical modeling which will be used is based on free convection heat transfer between two concentric spheres: from an inner small sphere to an outer larger sphere. Relating this model to our accelerometer, we state that the inner small sphere, with radius ri and temperature Ti, models the heater, whereas the outer larger sphere, with radius ro and temperature To, models the substrate. In a convective accelerometer design, the heater meander is implemented on top of the central suspended plate. Due to the high aspect ratio of this central plate thickness and width imposed by the technology, produced isotherm forms will certainly be far from being spherical or circular in a cross-sectional view. Here, we have to underline the fact that H. Hardee’s solution is based on the use of a spherical isotherm form produced by actual spheres. Consequently, this difference in isotherm forms imposes the modification of Hardee’s model equations so that they will be adapted to convective accelerometers. To give a further explanation of this isotherm form difference, we show in Fig. 2a and b, respectively, spheres ideal forms and their corresponding isotherms cross section that H. Hardee normally has used in his modeling to derive his solution. In Fig. 3a and b, we show, respectively, the x-axis and z-axis cross sectional views of the actual isotherms produced by a thermal accelerometer. Since the central plate is square, then isotherms cross section on both orthogonal horizontal axes are found to have the same shape. This figure clearly shows that sensor under test inner and outer sphere shapes are far from being ideal spheres. Therefore, to properly use spherical model expressions for convective accelerometers, both spheres deformations should be modeled and included in the governing analytical expressions. The basic idea behind the proposed solution is

Fig.  2 a Inner and outer spheres used in Hardee’s solution derivation, b isotherms produced by a cubic shaped form heater and symmetric top and bottom cavities design

(a)

(b)

Fig. 3 Produced isotherms by a convective accelerometer. Cross sections on the a x-axis and b orthogonal horizontal z-axis

to apply a mixed modeling scheme: use fitting expressions obtained from numerical FEM data in derived analytical expressions obtained from spherical model.

3 Conduction behavior modeling Harry C. Hardee studied in detail free convection heat transfer in a closed chamber (Hardee 1966). He expressed both common mode and differential temperatures as a function of parameters shown in Fig. 2a. As was previously shown in (Mezghani et al. 2013b, c, d), regarding convective accelerometers, both inner and outer spheres deformation have a large effect on temperature readings and therefore affects the overall sensor performance. Specifically, the outer sphere radius and form were found to be closely related to the sensor geometry and largely affects differential temperature readings. Among all sensor geometry design parameters, bottom cavity depth, h1, is the only distance which is very difficult to control during fabrication. In fact, this parameter is definitely fixed only at the end of the post-process etching step. This is because it can be easily affected by etching solution composition, etching time, etching solution movements and

13

348


temperature. Therefore, cavity depth value can be easily modified if ideal conditions are not met during the etching process. Based on these facts, our main modeling parameter will therefore be h1. Since we are dealing with a thermal device, therefore heating temperature values should be considered as the second important modeling parameter. 3.1 Analytical governing equations To study the effect of the h1 variable only, all other parameter values, listed in Table 1, which yield to sensitivity saturation, are used. Since we are interested in modeling the conductive behavior of the sensor, therefore we start our modeling by using the Fourier law for heat conduction, expressed as:

φ = −air A grad(T )

(1)

where λair is the air conductivity and A represents the area of a sphere with radius r. Air is used since it constitutes the gas filling the cavity of the sensor under test. Then, if we take the projection on the radial axis ur, we get the heat conduction equation in spherical coordinate:

φ = −air (T )4πr 2

dT dr

(2)

In the model, the heat conduction is studied from the inner to the outer spheres, so the radius variation should be set from ri to ro and temperature from Ti to To. In addition, expression of air conductivity from 100 K to 600 K is given by (Rekik et al. 2011): air (T ) = 0 . 1 + δ1 T + δ2 T 2 + δ3 T 3 (3) where λ0 = −3.93 × 10−4 W m−1 K−1 is the air conductivity extrapolated at 0 K (the negative sign is due to the fact that we extrapolate at 0 K an expression that is only valid from 100 to 600 K). The other three parameters are given by: δ1 = −0.259 K−1, δ2 = 1.23 × 10−4 K−2 and δ3 = −3.87 × 10−8 K−3, which represent the coefficients of thermal conductivity variation for air. Integrating, we get the expression of heat conduction distribution from ri to ro:

φ=

4π0 (Ti − To )C (T ) 1 ri

−

1 ro

(4)

where Cλ(T) is given by:

C (T ) = 1 +

δ2 Ti3 − To3 δ3 Ti4 − To4 δ1 Ti2 − To2 + + (5) 2 Ti − T o 3 Ti − T o 4 Ti − T o

In convective accelerometers, heat conduction is mainly governed by heater temperature, which is produced by an external power supply. In addition, this latter generated

13

temperature will set the initial temperature of both detectors, named also common mode temperature. Inside the micromachined cavity, heater temperature will be primarily affected by its heat transfer coefficient, which is closely related to cavity depth. Moreover, common mode temperature is set from heater temperature value. Therefore, to be able to model the conductive behavior of our sensor, we are to start modeling heater temperature through heat transfer coefficient. Then, through heat conduction in the cavity, we can analytically model common mode temperature. 3.2 Heater: Heat transfer coefficient In order for the convective accelerometer to operate, it is mandatory to generate a central hot bubble surrounding the heater. This is produced by applying a bias to the resistance forming the heater which is implemented on the central suspended plate. When this heater is powered with a dc voltage, Ui, an initial temperature distribution will be set in the cavity. The average heater temperature, Ti, is related to the electrical power, Pi, with:

Ti = To + Rthi Pi = To + Rthi

Ui2 Ri

(6)

where To is the substrate temperature, Ri is the heater electrical resistance and Rthi is its thermal resistance, defined as the ratio of temperature difference over heat conduction.

Rthi =

�T (Ti − To ) = φ φ

(7)

The electrical resistance is obviously independent of heater surroundings, but only depends on heating resistance size and material properties. In contrast, the heater thermal resistance depends not only on the dimensions of its resistance meander but also on its surroundings, and specifically h1. This is the case since all other silicon walls are quite far from the central plate. The Rthi parameter is closely related to the heat transfer coefficient (HTC), hH, which is defined as:

hH =

1 [W .m−2 .K −1 ] Rthi ∗ A

(8)

where A = 8.88 × 103 μm2 is the exchange surface area between the suspended plate and surrounding fluid (air) given by 8ri (ri + e). Applying the heat conduction expression of Eq. (4) on an inner isotherm representing the heater, with equivalent radius Ri, and an outer isotherm representing the hot bubble, with equivalent radius Ro, the HTC can be expressed as:

hH =

Ro 0 C (T ) Ri (Ro − Ri )

(9)


349

From Fig. 4, we can clearly see that the HTC is proportional to the heater temperature. This phenomenon is only due to the variation of air thermal conductivity which increases with temperature. It is also concluded that: • for high h1 values, the HTC is constant. This can be explained by the fact that silicon body walls are far enough offering no heat loss. • for low h1 values, the HTC increases. As bottom substrate silicon gets closer to the heater, more heat will have to be generated to compensate for the quantity lost in the substrate.

Fig. 4 Heat transfer coefficient as a function of cavity depth for different heater temperatures

We note that HTC dependence on heater temperature is represented in this equation through Cλ (T). From this equation, the Ro dependence on cavity depth, h1, can therefore be written as:

Ro (h1 ) =

hH (h1 )Ri2 hH (h1 )Ri − 0 C (T )

(10)

where hH (h1) denotes the hH values evaluated from FEM simulations for various technologically possible values of h1. Here, the outer sphere (outer isotherm) deformation effect is being modeled in Ro (h1) expression. Therefore, to be able to study the influence of cavity depth on the outer isotherm equivalent radius, Ro (h1), first we have to evaluate the HTC as a function of cavity depth, hH (h1). This latter study is done by performing a number of FEM simulations for various cavity depth values (50–300 μm) and also for different heater temperature values (350–600 K). Obtained simulation values are plotted in Fig. 4.

This can be explained by the shape of the hot bubble created by the heater, as illustrated in Fig. 5. In the case that cavity depth is high enough; the four lateral silicon cavity side walls will limit the hot bubble size. Therefore, hot bubble size is not affected by cavity depth value and the HTC just depends on heater temperature and cavity lateral dimensions. On the other hand, the hot bubble size reduces considerably for low cavity depth values, which becomes the main limiting factor of the produced isotherms. The transition between these two behaviors occurs around h1 = 100 μm (Fig. 4), this corresponds to a distance of (ro/3) from the heater center to the cavity border. It is very important to note that this transition value is independent of heater temperature. From HTC values evaluated from FEM simulations, Eq. (10) can now be used to calculate the values of outer isotherm equivalent radius as a function of cavity depths. Concerning the heater, we use an inner sphere equivalent radius of Ri = 28 μm. This value is a trade-off between a radius obtained from the central plate and that of the heater meander. Evaluated Ro values and fitting curve are plotted in Fig. 6. This curve fitting enables us to write an analytical model expression for the outer isotherm equivalent radius. Here, we express the outer isotherm equivalent radius Ro to reflect a linear relationship for low values of h1 and a

Fig. 5 Temperature profile cross section inside cavity for two cavity depths. a h1 = 300 μm, b h1 = 75 μm

13

350

Microsyst Technol (2015) 21:345–353 ′

where C (T ) is given by: ′

C (T ) = 1 +

3 − T3 2 − T2 4 − T4 δ2 TCM δ3 TCM δ1 TCM o o o + + 2 TCM − To 3 TCM − To 4 TCM − To

(13)

Since heat conduction is constant no matter which integration limits are used, therefore from Eq. (4), applied from Ri to Ro, and Eq. (12) we can express the dependence of TCM on Ti, given by: δ δ δ1 2 2 3 3 4 − To3 + − To4 TCM − To2 + TCM TCM 2 3 4 ′ δ δ1 2 Ri Ro (Ro − d) 2 Ti − To2 + Ti3 − To3 = ′ (Ti − To ) + 2 3 Ro d(Ro − Ri ) δ3 4 + (Ti − To4 ) 4 (14)

(TCM − To ) +

Fig. 6 Radius of equivalent sphere controlling heat transfer coefficient vs. cavity depth and its fitting

saturation at a value around Ro = 145 μm for high values of h1.

ro h1 h1 ∼ Ro = 144.27. = . 2 4 4 ro 4 4 h14 + (102.81)4 h1 + 3

(11)

In Eq. (11), the denominator root order was chosen to give the best transition fit between linear and saturation regions shown in Fig. 6. As a conclusion, the complete expression of heater HTC can be found by using the general fitting formula given by Eq. (11) in Eq. (9). 3.3 Common mode: fluid conduction 3.3.1 Spherical model and FEM data solution Detectors’ initial temperature is set by heat conduction in the air filling the cavity. This parameter is denoted common mode temperature TCM. This same concept of radial heat flow in spherical geometry will be used. The relationship between TCM and heater temperature Ti can then be obtained by solving heat conduction equation at a distance d from the heater center: detector’s location in the cavity, with Ri