384
MODELING HEAT TRANSFER IN THIN LAYERS
OF GRANULAR COMPOSITE MATERIALS BY CELLULAR
AUTOMATA
B.B.Khina and I.P.Samtsevich
Physico

Technical Institute, National Academy of Sciences, Minsk, Belarus
A new computer

oriented model for simulating heat transfer
in thin layers of granular
composite materials and calculating the effective thermal conductivity is developed
using the cellular automata method. Heat transfer in gaskets for heat removal from
heat

dissipating units used in Power Electronics is examined.
The layer thickness
accommodates only about 20 filler particles separated by the binder interlayers, and
the thermal conductivity coefficients of a filler (oxide or ceramic particles) and a
binder (epoxy or rubber resin) differ by two orders of magnitude,
thus the structure of
a composite (particle size, shape and spatial arrangement) affects substantially the
overall heat conductivity factor, especially for elongated filler particles. Numerical
simulation for AlN

filled silicon rubber is performed, and th
e effect of the spatial
arrangement of the filler particles on the effective thermal conductivity is evaluated.
INTRODUCTION
Composite materials with a polymeric binder and ceramic filler can possess a
unique combination of mutually contradicting propert
ies, e.g. high thermal
conductivity in combination with high insulation resistance. Thin layers (typically of
about 1 mm) of such materials are used as gaskets for heat removal from heat

dissipating units (transistors, microprocessors, etc.) in Power Elect
ronics. Elaboration
of new electronic devices for special applications necessitates the development of
novel composite materials with increased thermal conductivity retaining their
insulating properties. Also, there are stringent demands upon the mechanica
l
properties of such materials, viz. they must possess elasticity, plastic and adhesive
properties to fit to the roughness of the contacting surfaces (e.g., the transistor base
and radiator), along with durability and stability of thermal and mechanical pr
operties
at elevated temperatures (80

120
C in present

day devices, and up to 200
C in the
future prospect). For producing such composite materials, an electrically insulating
ceramic powder with relatively high thermal conductivity, e.g. BN, AlN, Al
2
O
3
,
ZnO
is used as a filler, and such substances as epoxy resin, synthetic rubber/silicon rubber
are employed as a matrix. It is known that the thermal conductivity,
, of the
constituents differ by two orders of magnitude: for example, for silicon rubber
=
(15.8
21.6)
10

2
W/(m
K) (1) while for Al
2
O
3
30 W/(m
K) (2). Besides, the
thermal properties of ceramic particles reported in literature vary in a wide range
depending on the production method and the crystallographic orientation of grains in
polycrys
talline samples used for measuring (e.g., for polycrystalline BN with wurtzite
lattice
varies from 30 to 60 W/(m
K) (3)).
Hence, further advance in this area demands for the development of
mathematical models capable of predicting heat conduction of such
composites
bearing in mind the ultimate goal of optimizing the structure of the material to attain
superior thermal properties. However, the specific feature of the problem is that the
385
layer thickness (H ~ 1 mm) accommodates only about 10 to 20 filler par
ticles with a
typical diameter of 20

50
m separated by interlayers of the binder (it is obvious that
the contact of several brittle ceramic particles inside the material should be avoided in
order to retain elastic properties). In this situation the analy
tical formulas widely used
for calculating thermal conduction of bulk materials (4,5 and many others), where
volume

averaging gives good agreement with experimental data, are not suitable.
Moreover, for elongated filler particles their spatial arrangement
must affect
substantially the overall thermal conductivity of the composite layer. It seems
reasonable to employ stochastic methods, which acquire wide use in modeling
disperse systems (6), for studying thermal properties of such materials.
In connection w
ith the above, the objective of this research is to develop a
computer

oriented method for modeling heat conduction in thin

layer composite
materials, where the thermal conductivity factors of the constituents differ by two
orders of magnitude, with a part
icular stress on simulating the effect of microstructure
on the target properties. To attain the goal, we use the "cellular automata" approach
(7) possessing wide potentiality for simulating heat/mass transfer in a domain with
complex internal structure.
FORMULATION OF A MODEL
Master equation
A two

dimensional domain is considered: a layer of a composite material with a
thickness H and width W. Conductive heat transfer through the material is described
by a Fourier equation written in a dimensionless for
m
(
,
)
(
,
)
(
,
)
(
,
)
x
y
c
x
y
x
x
y
x
y
x
y
y
(1)
with the initial condition
(x, y,
=0) =
0
= 0,
(2)
where x=X/L and y=Y/L are dimensionless coordinates where X and Y are the
corresponding dimensional coordinates, the 0X axis runs across the sample,
= t/t
0
is
the dim
ensionless time, L and t
0
are the characteristic size and time scale, respectively,
L = [t
0
0
/(
0
c
0
)]
1/2
;
0
,
0
and c
0
are the characteristic values of thermal conductivity,
density and specific heat, correspondingly; the dimensionless thermal parameters
are
defined as
=
/
0
,
=
/
0
,
c = c/c
0
,
= (T
T
0
)/(T*
T
0
), where T is temperature
and T* is the characteristic maximal temperature (T*=80

120
C), T
0
is the initial
temperature which is equal to the ambient temperature.
At the upper surface of th
e flat sample (x=0) we pose the 2nd order boundary
condition to equation (1):
x
q
x
x
0
0
(
)
(3)
386
where q
x=0
is the dimensionless heat flux per unit area dissipated by the electronic
device attached to the upper surface, q = QL/[
0
(T*
T
0
)], Q is
the corresponding
dimensional parameter. For power transistors used e.g., in TV sets, the typical value is
Q = 3 to 7 W/cm
2
(8). Equation (3) permits simulating heat removal not only in the
steady state thermal regime (when Q=const) but also during the ini
tial rapid heating of
an electronic device, because most failures occur during transient regimes.
The lower surface (X=H, or x=h=H/L) is considered to be kept at the ambient
temperature T
0
(1st order boundary condition):
(x=h) = 0
(4)
At the lateral s
urfaces y=0 and y=w=W/L the adiabatic boundary conditions (the
absence of heat exchange with the surroundings) are postulated:
y
y
y
y
w
0
0
(5)
Implementation of cellular automata model
The cellular automata model for numerical solution of mast
er equation (1) with
initial (2) and boundary conditions (3)

(5) is built as follows. Since the temperature
interval considered in the given problem is relatively narrow, about 100

150
C, the
temperature dependence of all the physical parameters may be ne
glected, and thus the
model is treated as a linear one with respect to temperature. A discrete equation
describing heat exchange between adjacent cells is actually an analog of an explicit
finite

difference scheme, which can be obtained by the integration

interpolation
method (9). Using simple formulas
(i
1,j) = (1/2)[
(i,j) +
(i
1,j)] and
(i,j
1) =
(1/2)[
(i,j) +
(i,j
1)] for interpolating the
value in half

points, where i and j denote
the cell coordinate along the 0x and 0y axes, correspondingly, we
derive the equation
c
i
j
i
j
i
j
n
i
j
n
a
i
j
i
j
i
j
n
i
j
n
a
i
j
i
j
i
j
n
i
j
n
a
i
j
i
j
i
j
n
i
j
n
a
(
,
)
(
,
)
(
,
,
)
(
,
,
)
[
(
,
)
(
,
)]
(
,
,
)
(
,
,
)
[
(
,
)
(
,
)]
(
,
,
)
(
,
,
)
[
(
,
)
(
,
)]
(
,
,
)
(
,
,
)
1
1
2
1
1
1
2
1
1
1
2
1
1
1
2
[
(
,
)
(
,
)]
(
,
,
)
(
,
,
)
i
j
i
j
i
j
n
i
j
n
a
1
1
(6)
where
is a temporal step and n denotes the step number, a is the cell size.
For explicit schemes, the stability criterion exists,
=
/[(a
2
min(
c
)]
0.25,
which imposes limitation on the
value.
A two

color picture representing the structure of a composite material (a model
image or a real microstructure entered into computer from an optical microscope) is
superimposed on a field composed of square cells (not less then 100
100 cells), and
the ther
mal parameters prescribed to a particular cell, viz.
c(i,j),
(i,j) and
(i,j) are
387
associated with the color of the corresponding image fragment (e.g., black means
matrix and white denotes filler). After that, a numerical procedure is invoked which
simulat
es the evolution of the temperature field according to Eq.(6) with
corresponding boundary conditions. Calculations are continued until the steady

state
regime of heat transfer is attained.
The effective thermal conductivity factor of the composite material
is calculated
for the steady

state regime when the heat flux, Q, becomes constant along the 0X axis:
eff
Q
T(
X
T
H
(
)
)
/
0
0
(7)
Local thermal conductivity in a longitudinal section (parallel to the 0Y axis),
s
(X),
which must depend on the volume fracti
on, V, of the filler phase in the corresponding
section, was determined in the course of simulation by the formula
s
W
W
X
X
Y
T(
X
Y
X
dY
T(
X
Y
X
dY
(
)
(
,
)
,
)
,
)
0
0
(8)
RESULTS OF COMPUTER SIMULATION AND DISCUSSION
Calculations were performed for AlN

filled silicon rubber whose t
hermal
properties are listed in Table
1
. For boundary condition (3) we took the dimensional
heat flux Q=5 W/cm
2
(see Ref.(8)).
Table
1
. Thermal parameters of the constituents
Substance
Ⱐ,⽣/
3
c, J/(g
䬩
Ⱐ,⼨/
䬩
Re晥牥湣es
䅬A
㌮ㄲ
〮㜳0
㌰⸰
⠳(
s
楬楣潮畢扥
〮㤷0
ㄮ1
〮ㄷ
⠱(
F楧畲猠
1
through
3
present the results of computer simulation of heat transfer in
thin layers of composite materials having a close volume fraction of the filler (about
40%) and similar particle shape and size but differ
ent microstructure, i.e. spatial
arrangement of the filler particles: ordered structure with elongated particles placed
upright (Fig.
1
,a) and aflat (Fig.
2
,a), and disordered structure (Fig.
3
,a). The
corresponding temperature maps (in grayscale) after the a
ttainment of a steady

state
regime are shown in Figs.
1
,b through
3
,b, the temperature

scaling palette being placed
below each figure (darker color denotes higher temperature). The lines separated
differently colored areas in Figs.
1
,b
3
,b are isotherms.
Fro
m the comparison of the composite structures and the corresponding
temperature maps it is seen that hot (i.e., superheated) areas in the upper part of the
specimen appear where the heat

conductive particles are absent (Fig.
2
,b).
388
Fig.
1
. Results of model
ing heat
transfer in a composite layer with the
filler particles placed upright: (a)
microstructure, (b) temperature map
Fig.
2
. Results of modeling heat
transfer in a composite layer with the
filler particles placed aflat: (a)
microstructure, (b) t
emperature map
The average temperature at the upper surface of the composite layer after the
attainment of the steady

state regime and the effective thermal conductivity calculated
by Eq.(7) were the following: T(X=0) = 70
C and
eff
= 0.69 W/(m
K) f
or upright
arrangement of particles (Fig.
1
,a), T(X=0) = 115
C and
eff
= 0.44 W/(m
K) for
structure with particles placed aflat (Fig.
2
,a), and T(X=0) = 85
C and
eff
= 0.59
389
W/(m
K) for disordered microstructure (Fig.
3
,a). Thus, the calculated
eff
value
varies
in a relatively wide range [

25%, +17%] with respect to
eff
for a sample with chaotic
distribution of the filler particles. So, the effect of spatial arrangement of the filler
particles in thin layers of composite materials on the thermal conductio
n is substantial
for the conditions considered in this work, viz. heat removal from power electronic
devices. The optimal microstructure is that presented in Fig.
1
,a. Low thermal
conductivity of the structure shown in Fig.
2
,a is connected with the presence
of
continuous horizontal interlayers of the low

conducting binder.
Fig.
3
. Results of modeling heat
transfer in a composite layer with
chaotically oriented filler particles:
(a) microstructure, (b) temperature
map
For comparison, a volume

averaged
formula (Ref.(4)) for a bulk material with
cubical inclusions
eff
=
m
m
V/[(1
f
/
m
)
1
(1
V)/3]
(9)
where
m
and
f
are thermal conductivity factors of the matrix and filler,
correspondingly, yields
eff
= 0.5 W/(m
K) for the same volume fr
action of the filler
(V=0.4). This differs by 15% from the calculated value for the disordered
microstructure (Fig.
3
,a).
CONCLUSION
A new approach for studying conductive heat transfer through thin

layer
composite materials with greatly differing thermal
conductivity factors of the
constituents is developed. It is demonstrated that spatial arrangement of high

conducting filler particles in a thin

layered composite material used as a heat

removing gasket has a substantial influence on the overall thermal c
onduction: the
390
calculated effective thermal conductivity factor varies in a relatively wide range with
changing the microstructure of the material. This is important for power electronic
devices, especially for those used in special applications.
The resu
lts of cellular automata modeling lay a basis for optimizing the structure
of a composite material in order to improve its thermal properties.
In this work, we implied perfect adhesion at the filler particle/binder interface,
i.e. an absence of an interfac
ial thermal barrier (the effect of such barrier has been
evaluated analytically in Ref.(5) for a volume

averaged approach, which is applicable
for bulk composite materials). However, the proposed model can be generalized to
take into account this effect ju
st by including chaotically distributed thin interlayers of
an additional low

conducting phase at the particle/matrix interfaces.
REFERENCES
1.
Shatz M: Silicon rubber. Prague, SNTL, 1971 (in Czech; Russian translation:
Moscow, Khimiya, 1975).
2.
Chudnovskiy
A F: Thermophysical properties of disperse materials. Moscow,
Fizmatgiz, 1962 (in Russian).
3.
Kosolapova T Ya, Andreeva T V, Bartnitskaya T S, Gnesin G G, Makarenko G
N, Osipova I I, Prilutskiy E V: Nonmetallic refractory compounds. Moscow,
Metallurgia, 198
5 (in Russian).
4.
Dul'nev G N, Zarichnyak Yu P: Thermal conductivity of mixtures and
composite materials: a handbook. Leningrad, Energiya, 1974 (in Russian).
5.
Lipton R: 'Design of particle reinforced heat conducting composites with
interfacial barriers'. Jour
nal of Composite Materials 1998 32 (14) 1322

31.
6.
Mavrin S V, Stengach A V, Potanin A A: 'Stochastic model of disperse systems'.
Inzhenerno

Fizicheskiy Zhurnal 1999 72 (2) 245

50 (in Russian; English
translation: Journal of Engineering Physics (USA)).
7.
Wolfr
am S (ed.): Theory and applications of cellular automata. World
Scientific: Singapore, 1986.
8.
Borodin B A, Lomakin V M, Mokryakov V V, Petukhov V M, Khrulev A K:
Power semiconductor devices: transistors. Moscow, Radio i Svyaz, 1985 (in
Russian).
9.
Samarskiy A
A: Introduction to the theory of finite

difference schemes.
Moscow, Nauka, 1971 (in Russian).
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο