INVERSE THERMAL PROBLEM VIA MODEL ORDER REDUCTION:
DETERMINING MATERIAL PROPERTIES OF A MICROHOTPLATE
Tamara Bechtold, Dennis Hohlfeld, Evgenii B. Rudnyi, Hans Zappe and Jan G. Korvink
IMTEK, University of Freiburg, 79110 Freiburg, Germany, Email: bechtold@imtek.de,
Phone: +49 761 203 7387, Fax: +49 761 203 7382
ABSTRACT
In this paper we present a new possibility to determine
unknown material thermal properties based on experimen-
tal measurements. We use a parametrized finite element
(FE) model and fit it to a measured transient curve. Using
brute force, this process requires an extensive computa-
tional time. However, it is possible to speed it up signifi-
cantly by using model order reduction (MOR).
1. INTRODUCTION
An important engineering task is to built a validated model
for each characterized novel MEMS device. It is possible
to fit an RC-ladder network to the measured results [1], but
in this way we don?t obtain a complete physical picture of
the device. As in most MEMS applications the whole tem-
perature field has to be known, a more detailed FE model
is required. Unfortunately, a common problem here is that
the material properties of the employed thin film materials,
like e. g. thermal conductivity and heat capacity ,
strongly depend on fabrication conditions and may also be
specific for the device under the test. In such case, it is pos-
sible to extract the material properties by fitting a parame-
trized FE model to a measured transient curve. However,
the conventional optimization process is highly time con-
suming, because in each iteration a time integration of a
full-scale model must be performed (see Fig. 1). As in
most MEMS applications the size of an accurate finite ele-
ment model easily exceeds 100,000 ordinary differential
equations (ODEs), we suggest an alternative approach
based on model order reduction. The right path in Fig. 1
shows that in each iteration of optimization loop, the sug-
gested approach requires only the time integration of the
reduced model (with less than 50 ODEs) and hence brings
along an enormous saving in computational time. By
defining an objective function, which characterizes the dif-
ference between simulated and measured results, data fit-
ting cycle is performed. In this paper we demonstrate that
Arnoldi-based model order reduction combined with opti-
mization can be used for the efficient extraction of material
thermal properties of a microhotplate. We limit ourselves
to a simplified two dimensional (2D) finite element model
with the goal to test the computational environment.
2. MEMS CASE STUDIE - SILICON-BASED
MICROHOTPLATE
The fabricated structure resembles a micro-hotplate since
it features a membrane for thermal isolation and integrated
resistors for heat generation (see Fig. 2). This class of
structures is employed in a variety of other micro-fabri-
cated devices such as gas sensors [3] and infrared sources
[4], [5].
The membrane features a thin-film metal resistor for tem-
perature modulation through Joule heating. A second resis-
tor, which serves as a temperature sensor, is placed
adjacent to the heater. This resistor is configured for a four-
point measurement of its resistance. In order to achieve a
preferably circular symmetric and homogenous tempera-
? c
p
Fig. 1 Parameter extraction process for transient thermal
problems via optimization and model order reduction (see
also [2]).
time integration
of the full-scale FE model
evaluate the objective
function
build a parametrized FE model
with parameters k and r*Cp
change parameter
values
convergence criterium
fulfilled?
yes
no
start
stop
model order reduction
time integration
of the reduced model
conventional
approach
suggested
approach
Belgirate, Italy, 28-30 September 2005
ture distribution at the center of the square membrane, both
resistors are arranged as shown in Fig. 3.
The characterization of the static and transient thermal
properties of the filter membrane is performed on a tem-
perature controlled mount. For characterization of
dynamic temperature changes, a constant current of
100 is passed through the outer terminals of the sens-
ing resistor while the voltage between the inner terminals
is measured using an oscilloscope. In order to use the sens-
ing resistor as a temperature sensor, the linear temperature
coefficient of the material?s resistivity has to be known.
The temperature coefficient is measured by acquiring the
sensor?s resistance at various temperatures. These are pre-
cisely set by the Peltier-mount. The electrical resistance
depends linearly on the temperature over the investigated
temperature range. This is modeled through:
(1)
with R
0
as the resistance for ?T = 0. A temperature coeffi-
cient of ? = 2.293K
-1
is obtained for a metallization of
150nm platinum with 50nm titanium. The electro-thermal
dynamics of the system is determined through the specific
heat capacity, material density and thermal conductivity of
the membrane material. Further influences on transient
thermal behavior are the geometry parameters, like mem-
brane shape and location of resistors. The transient thermal
response of the filter membrane is characterized by apply-
ing rectangular heat pulses to the heating resistor using a
function generator. The signal output is configured as a
voltage source with a fixed output impedance of 50 .
After applying the heating power, the membrane?s temper-
ature increases until a maximum value is reached. This
temperature is defined as the steady-state value. After set-
ting the power to zero, the heat stored in the membrane?s
volume is dissipated to the surrounding media by conduc-
tion and free convection. Thus, the temperature drops
down to its initial value. The thermal response over a
whole period is presented in Fig. 4.
As already mentioned, we use a 2D model whose FE mesh
is shown in Fig. 5. This simplification requires the omis-
sion of both out-of-plane thermal conduction and convec-
tion from the membrane?s surface. Model contains 4,402
nodes and Dirichlet boundary conditions are set
at the outer edges of the simulation domain which results
in 4,182 ODEs.
The layer parameters are differentiated into three sections:
areas where only the membrane layers are present, sections
with the metal thin-films on the membrane and the silicon
frame. The material properties are thickness related values,
which are calculated as:
Fig. 2 Microstructured silicon/nitride membrane.
Fig. 3 Schematic view of the thin-film resistors for heating
and temperature sensing. Heating resistor is operated at
constant voltage. The sensing resistor is configured for
four-point measurement.
tunable cavity
thin-film heater
silicon substrate
air
silicon nitride membrane
100 ?m
sense
resistor
heating
resistors R
heat1
R
heat2
I
1
I
2
V
sense1
V
sense2
R
heat1
R
heat2
I
1
I
2
V
sense1
V
sense2
?A
R ?T()R
0
1 ??T+()=
Fig. 4 Temperature modulation of a silicon/nitride mem-
brane. A square wave heating power with a frequency of
25 Hz was applied to the membrane.
?
temperature change (K)
time (ms)
heating power (mW)
010203040
0
10
20
30
40
0.0
1.5
3.0
T 0?C=
, , (2)
where d
i
is a thickness of each layer and , and are
its thermal conductivity, mass density and heat capacity
respectively. In order to achieve the consistency between
the numerical model and the measurement data, the effec-
tive thermal conductivity and volumetric heat capacity
are set as fit parameters in order to determine
their actual values via optimization of the reduced model.
3. MODEL ORDER REDUCTION
MOR allows a formal conversion of the physical model,
that is, governing partial differential equation to a low-
dimensional ordinary differential equation system (see Fig.
1). The intermediate level is a device level, that is, a high
dimensional ODE system. The first conversion of the
physical to the device model we perform via the finite ele-
ment (FE) discretisation. The heat transfer within a hot-
plate is described through the following equations:
, (3)
where is the thermal conductivity in W/mK at the
position r, is the specific heat capacity in J/kgK,
is the mass density in kg/m
3
, is the tempera-
ture distribution and is the heat generation rate per
unit volume in W/m
3
. Assuming that the heat generation is
uniformly distributed within the heater, and that the system
matrices are temperature independent around the working
point, the finite element based spatial discretization of (3)
leads to a large linear ODE system of the form:
(4)
where C and K are the global heat capacity and heat con-
ductivity matrices, F is the load vector (matrix) and E is
the output vector. (4) is a starting point for model order
reduction. As already mentioned for our test case it con-
tains 4,402 ODEs. Software tool mor4ANSYS [6] uses
Arnoldi reduction algorithm, which can be viewed as a
projection, from the full space to the reduced Krylov-sub-
space:
(5)
with and . This projection is
based on the transformation of the state vector T to the vec-
tor of generalized coordinates z, subjected to some small
error :
(6)
and subsequent left side multiplication of (4) with the V
T
.
The transformation matrix , where are the
dimensions of the reduced and the full system, respec-
tively, is a direct output of Arnoldi algorithm. The property
of the Krylov-subspace (5) is such that the transfer func-
tion of (4) is approximated through the first r coefficients
of its Taylor series around an arbitrary chosen frequency.
Figure 3 in [7] shows the transfer function of the full and
reduced order model of a similar microhotplate-based
structure, for the expansion around zero. Due to choice of
the expansion point we observe a good approximation in
the low frequency domain, as can be expected. This corre-
sponds to the good approximation of the steady-state in
transient thermal response (shown in Figure 4 in [7]).
Fig. 5 FE mesh of the 2D model with 4,402 nodes.
Fig. 6 Place of model order reduction within a conversion
process from physical to compact model.
?'
m
?
i
d
i
i
?
= ?'
m
?
i
d
i
i
?
= c'
m
1
?
i
d
i
i
?
----------------- c
i
?
i
d
i
i
?
?=
?
i
?
i
c
i
?'
m
c'
m
?'
m
?
FEM MOR
? x
r
= A
r[]
x
r
+b
r
u
y
r
=c
r
T
x
r
Reduced system
of r << n ODEs
System of
n ODEs
Physics &
Geometry
physical level device level system level
??T?()? Q ?c
t?
?T
?+ 0= Q j
2
R=
? r()
C
p
r()
? r() Trt,()
Qrt,()
CT
?
? KT?+ FI
2
t()RT()?=
yE
T
T?=
K
r
Ab,{}span b A
2
b ? A
r 1?
b,,,
??
??
??
=
AK?
1?
C= bK?
1?
F=
?
TVz?+?=
VR
nxr
? rn?
Apart from choosing an expansion point, the MOR user
must also chose an order for the reduced model which will
fulfill the desired accuracy. In [7] we have suggested an
error indicator criterion which is based on the convergence
of relative error between two reduced models with subse-
quent order. In case of the microhotplate order 10 has been
chosen.
Concerning the saving in computational times due to
MOR, it is worth of noting that the time integration of the
2D microhotplate model in ANSYS with 30 time-steps
takes 150s, whereas order reduction and the subsequent
time integration of the reduced order 10 model in mathe-
matica last only about 2s. However, the writing of files in
each time-step by ANSYS, takes up to 50% of integration
time for small model sizes. For the higher dimensional
models the ratio between ANSYS integration and MOR
process is smaller, so that the total saving in computational
effort is about 20 times. In [8] we have shown that the time
needed for the Arnoldi-based model order reduction of
thermal MEMS models is comparable with a single sta-
tionary solution of the original system. We use this enor-
mous increase in efficiency to extract material thermal
parameters via optimization.
4. OPTIMIZATION
Fig. 7 shows the flexible optimization environment cou-
pled to MOR process. Mathematica is used for scripting,
visualization and small size computations. Its function eval
takes as arguments the fitting parameters and and
calls the external programs ANSYS (for rebuilding the FE
model with changed material properties) and mor4ansys
(for creating a reduced model). It further integrates the
reduced model and evaluates the objective function, which
is defined as a quadratic error between the measured and
the computed curves. Its value is transferred back to DOT
optimizer [9] which communicates with Mathematica via
Mathlink (our implementation can be found at http://evge-
nii.rudnyi.ru/soft/dot/).
Fig. 8 shows the measured temperature response and the
simulated temperature response for the reduced order 10
model before the optimization (initial values for the mate-
rial parameters were chosen as and
). Fig. 9 shows the mea-
sured and the simulated temperature response after 35
cycles of optimization (end values for the material parame-
ters were and
).
Fig. 7 Implementation of the reduced model optimization.
Mathematica
eval
(?,c
p
)
Mathlink
DOT optimizer
ANSYS
mor4ansys
integrate reduced
model, compute
TF = (T
measured
?
reduced
)
2
?
obj. func.F
call eval
call DOT
? c ??
Fig. 8 Measured and simulated temperature response be-
fore optimization.
Fig. 9 Measured and simulated temperature response after
35 cycles of simulation.
? 6.35WmK??=
? c
p
? 525 10
3
Jm
3
K?()??=
? 5WmK??=
? c
p
? 662 10
3
Jm
3
K?()??=
010203040
0
10
20
30
time (ms)
temperature change (K)
measured
simulated
010203040
0
10
20
30
time (ms)
temperature change (K)
measured
simulated
5. CONCLUSION
We have developed a software environment which can be
used for the solution of the inverse thermal problem via
model order reduction and optimization. We have applied
it to extract material thermal parameters of a fabricated
microhotplate. Due to the fact that we have used a simpli-
fied 2D model under neglecting out-of-plane thermal con-
duction and convection and that we have performed
optimization with only 2 parameters, the achieved numeri-
cal values may not be correct. Nevertheless, we have dem-
onstrated that the tool works in principle. The next step
should be to apply it the more realistic 3D model under
consideration of convection. It is further possible to change
the parameter values at the level of the reduced model [10]
(see the right short path in Fig. 1) and to so completely
avoid building of a new FE model and subsequent model
order reduction.
11. REFERENCES
[1] V. Sz?kely, M. Renz, ?Thermal dynamics and the
time constant domain?, IEEE Trans. on Comp. Pack.
Techn., 23(3), pp. 587-594, (2000).
[2] J. S. Han, E. B. Rudnyi, J. G. Korvink, ?Efficient
optimization of transient dynamic problems in
MEMS devices using model order reduction?, J.
Micromech. Microeng., 15(4), pp. 822-832, (2005).
[3] M. Graf, D. Barrettino, S. Taschini, C. Hagleitner,
A. Hierlemann, H. Baltes, "Metal oxide-based
monolithic complementary metal oxide semicon-
ductor gas sensor microsystem", Analytical Chemis-
try, 76, pp. 4437-45, (2004).
[4] R. W. Bernstein, A. Ferber, I. Johansen, S. T. Moe,
H. Rogne, and D. T. Wang, "Optical MEMS for
infrared gas sensors", Proc. IEEE/LEOS Interna-
tional Conference on Optical MEMS, pp. 137-8,
(2000).
[5] A. Pollien, J. Baborowski, N. Ledermann, and P.
Muralt, "New material for thin film filament of
micromachined hot-plate", Proc. 11th International
Conference on Solid-State Sensors and Actuators,
Munich, Germany, pp. 848-51, (2001).
[6] E. B. Rudnyi, J. G. Korvink, ?mor4ansys (version
1.8): Compact Behavioral Models from ANSYS by
Means of Model Order Reduction?, at http://
www.imtek.uni-freiburg/simulation/mor4ansys,
(2005).
[7] T. Bechtold, E. B. Rudnyi, J. G. Korvink, ?Error
indicators for fully automatic extraction of heat-
transfer macromodels for MEMS?, J. Micromech.
Microeng., 15(3), pp. 430-440, (2005).
[8] T. Bechtold, ?Model Order Reduction of Electro-
Thermal MEMS?, PhD thesis, University of
Freiburg, (2005).
[9] DOT Users Manual version 4.20 (Colorado Springs,
CO: Vanderplaats Research and Development), at
http://www.vrand.com/DOT.html
[10] L. H. Feng, E. B. Rudnyi, J. G. Korvink, ?Preserving
the film coefficient as a parameter in the compact
thermal model for fast electro-thermal simulation?,
IEEE Transactions on Computer-Aided Design of
Integrated Circuits and Systems, to be published,
(2005).