Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

2, 40, 2002

NUMERICAL MODELLING OF HEAT TRANSFER

IN SPHERICAL DOMAINS BY MEANS OF THE BEM

USING DISCRETISATION IN TIME

Bohdan Mochnacki

Romuald Szopa

Institute of Mathematics and Computer Sciences, Technical University of Częstochowa

e-mail: moch@matinf.pcz.czest.pl

A combined variant of the BEM called in literature the BEM using di-
scretization in time consists in an approximation of the time derivative
appearing inFourier’s equationbyanadequate differential quotient.The
next steps ofmathematicalmanipulations and also the numerical algori-
thmare similar to a typical boundaryelement approach. In thepaper the
method is applied to numerical computations concerning a non-steady
heat diffusion in homogeneous and non-homogeneous spherical domains.
In the final part of the paper the results of computations are presented.

Key words: heat transfer, boundary element method

1. Introduction

At first, the well known linear Fourier equation for 3D domain oriented in
Cartesian co-ordinate system is considered

x∈Ω :
∂T(x, t)

∂t
= a

3
∑

e=1

∂2T(x, t)

∂x2e
= a∇2T(x, t) (1.1)

where x = (x1,x2,x3), a = λ/c is the heat diffusion coefficient (λ is the
thermal conductivity, while c is the specific heat per unit volume), T , tdenote
temperature and time, respectively. On the outer surface Γ of the system
boundary conditions are given, the initial condition is also known.
In this place a time grid with a constant step ∆tmust be introduced

0= t0 <t1 <t2 < ... < tf−1 <tf < ... < tF <∞ (1.2)



358 B.Mochnacki, R.Szopa

Considering the transition tf−1 → tf one transforms equation (1.1) to the
form (Curan et al., 1980; Sichert, 1989)

T(x, tf)−T(x, tf−1)
∆t

= a∇2T(x, tf) (1.3)

or

∇2T(x, tf)−
1

a∆t
T(x, tf)+

1

a∆t
T(x, tf−1)= 0 (1.4)

Using the weighted residual criterion one obtains (Brebbia et al., 1984;
Majchrzak andMochnacki, 1995)

∫

Ω

[

∇2T(x, tf)− 1
a∆t
T(x, tf)+

1

a∆t
T(x, tf−1)

]

U∗(d) dΩ =0 (1.5)

where U∗ is a fundamental solution, in particular to the considered task it is
the following function (Brebbia et al., 1984)

U∗(d) =
1

4πd
exp
(

− d√
a∆t

)

d=

√

√

√

√

3
∑

e=1

(xe− ξe)2 (1.6)

where ξ = (ξ1,ξ2,ξ3) is a point at which the concentrated heat source is
applied (Brebbia et al., 1984).

The function U∗ fulfills the equation

∇2U∗(d)−
1

a∆t
U∗(d) =−δ(ξ,x) (1.7)

where δ(ξ,x) is the Dirac function.

Using the 2nd Green formula and property (1.7) one obtains the WRM
criterion in the form

T(ξ, tf)−
1

λ

∫

Γ

U∗(d)q(x, tf) dΓ =

(1.8)

=
1

λ

∫

Γ

Q∗(d)T(x, tf) dΓ +
1

a∆t

∫

Ω

U∗(d)T(x, tf−1) dΩ

where q=−λ∂T/∂n, Q∗ =−λ∂U∗/∂n.



Numerical modelling of heat transfer... 359

2. Homogeneous spherical domain

According to the idea presented byBrebbia et al. (1984) (the basic variant
of the BEM and cylindrical object has been considered), equation (1.8) must
be integrated assuming that Ω corresponds to the interior of the cylinder,
while Γ to its surface. The similar approach has been proposed by Bokota
(1989) in relation to the spherical domain.

In the paper this concept is applied in order to obtain the boundary equ-
ation in the case of BEM using discretization in time.

The spherical co-ordinate system should be introduced

x1 = rcosϕsinθ x2 = r sinϕsinθ x3 = rcosθ (2.1)

Additionally, we assume the position of Cartesian system for which ξ =
(ξ1,ξ2,ξ3) = (0,0,ξ) and then that the distance between this point and the
point considered x (see Fig.1) is equal to

d=
√

r2+ ξ2−2rξcosθ (2.2)

Fig. 1. Spherical co-ordinate system

The surface and volume elements of the sphere can be expressed as follows

dΓ =R2 sinθ dθdϕ dΩ= r2 sinθ dθdϕdr (2.3)

where R is the radius of the sphere.



360 B.Mochnacki, R.Szopa

Thus, equation (1.8) takes the form

T(ξ,tf)+
R2

λ
q(R,tf)

2π
∫

0

π
∫

0

U∗(d)sinθ dθdϕ=

=
R2

λ
T(R,tf)

2π
∫

0

π
∫

0

Q∗(d)sinθ dθdϕ+ (2.4)

+
1

a∆t

R
∫

0

r2
2π
∫

0

π
∫

0

U∗(d)sinθ dθdϕT(r,tf−1) dr

or

T(ξ,tf)+
R2

λ
T∗(ξ,R)q(R,tf)=

(2.5)

=
R2

λ
q∗(ξ,R)T(R,tf)+

1

a∆t

R
∫

0

r2T∗(ξ,r)T(r,tf−1 dr

where

T∗(ξ,r)=

2π
∫

0

π
∫

0

U∗(d)sinθ dθdϕ=

=
1

4π

2π
∫

0

π
∫

0

1
√

r2+ ξ2−2rξ
exp
(

−
√

r2+ ξ2−2rξ
a∆t

)

sinθ dθdϕ

(2.6)

q∗(ξ,r) = −λ
∂T∗(ξ,r)

∂r

After the integration one obtains

T∗(ξ,r) =

√
a∆t

2rξ

[

exp
(

−|r− ξ|
a∆t

)

− exp
(

−|r+ ξ|
a∆t

)]

(2.7)

q∗(ξ,r) =
λ

r
T∗(ξ,r)+

λ

2ξr

[

sgn(r−ξ)exp
(

−|r− ξ|
a∆t

)

−

− sgn(r+ξ)exp
(

−
|r+ ξ|
a∆t

)]



Numerical modelling of heat transfer... 361

The boundary equation (for ξ→R) is of the form (cf. equation (2.5))

T(R,tf)+
R2

λ
T∗(R,R)q(R,tf)=

(2.8)

=
R2

λ
q∗(R,R)T(R,tf)+

1

a∆t

R
∫

0

r2T∗(R,r)T(r,tf−1) dr

The first stage of numerical computations consists in determination of the
boundary heat flux (if the boundary temperature is given) or boundary tem-
perature (if the boundary heat flux is given) – equation (2.9). In the second
stage the temperatures at the set of internal points ξ∈ (0,R) for time the tf
can be found on the basis of equation

T(ξ,tf)=
R2

λ
q∗(ξ,R)T(R,tf)−

(2.9)

−R
2

λ
T∗(ξ,R)q(R,tf)+

1

a∆t

R
∫

0

r2T∗(ξ,r)T(r,tf−1) dr

The solution obtained constitutes a pseudo-initial condition for the next lo-
op of computations. It should be pointed out that the integral appearing in
equations (2.8) and (2.9) can be found using numerical methods e.g. Gaussian
quadratures. The integral over the first internal cell (r∈ [0,∆r]) is a singular
one, but in a numerical realisation it does not cause essential difficulties.

3. Non-homogeneous spherical domain

The non-steady temperature field in the domain considered is described
by a system of equations

Rm−1 <r<Rm :
∂Tm(r,t)

∂t
=
am
r2
∂

∂r

[

r2
∂Tm(r,t)

∂r

]

(3.1)

where m=1,2, ...,M.
For r=Rm,m=1,2, ...,M −1 the continuity conditions in the form

r=Rm :











−λm
∂Tm(r,t)

∂r
=−λm+1

∂Tm+1(r,t)

∂r

Tm(r,t) =Tm+1(r,t)

(3.2)



362 B.Mochnacki, R.Szopa

are given. For r =R0 and r=RM the boundary temperatures or boundary
heat fluxes are known. For the time t = 0 the initial temperatures are also
given.
Equation (2.5) for the spherical shell r∈ (Rm−1,Rm) takes the form

Tm(ξ,t
f)+

[ r2

λm
T∗m(ξ,R)qm(r,t

f)
]Rm

Rm−1
=

(3.3)

=
[ r2

λm
q∗m(ξ,r)Tm(r,t

f)
]Rm

Rm−1
+
1

am∆t

Rm
∫

Rm−1

r2T∗m(ξ,r)Tm(r,t
f−1) dr

Equation (3.3) can be written as follows

Tm(ξ,t
f)+qm(ξ,Rm)qm(Rm, t

f)−qm(ξ,Rm−1)qm(Rm−1, tf)=
(3.4)

=hm(ξ,Rm)Tm(Rm, t
f)−hm(ξ,Rm−1)Tm(Rm−1, tf)= pm(ξ)

where

qm(ξ,r)=
r2

λm
T∗m(ξ,r) hm(ξ,r)=

r2

λm
q∗m(ξ,r) (3.5)

while

pm(ξ)=
1

am∆t

Rm
∫

Rm−1

r2T∗m(ξ,r)Tm(r,t
f−1) dr (3.6)

For ξ→R+m−1 and ξ→R
−

m−1 one obtains a system of equations

[

gm(R
+
m−1,Rm−1) gm(R

+
m−1,Rm)

gm(R
−

m,Rm−1) gm(R
−

m,Rm)

][

qm(rm−1, t
f)

qm(Rm, t
f)

]

=

=

[

hm(R
+
m−1,Rm−1)−1 hm(R

+
m−1,Rm)

hm(R
−

m,Rm−1) hm(R
−

m,Rm)−1

][

Tm(Rm−1, t
f)

Tm(Rm, t
f)

]

+(3.7)

+

[

pm(Rm−1)

pm(Rm)

]

or
[

gm11 g
m
12

gm21 g
m
22

][

qm(Rm−1, t
f)

qm(Rm, t
f)

]

=

[

hm11 h
m
12

hm21 h
m
22

][

Tm(Rm−1, t
f)

Tm(Rm, t
f)

]

+

[

pm1

pm2

]

(3.8)



Numerical modelling of heat transfer... 363

Thefinal system formulti-layers domain results fromthe coupling of equations
(3.8) by the continuity conditions given for r=Rm,m=1, ...,M −1.
For example, we consider the two-layer system (M =2), Figure 2, andwe

assume the boundary conditions for r=R0 and r=R2 in the form

r=R0 : q1(r,t
f)= qb

r=R2 : q2(r,t
f)=α[T2(r,t

f)−T∞]
(3.9)

where α is the heat transfer coefficient and T∞ is the ambient temperature.
The continuity condition for r=R1 can be written in the form (cf. equation
(3.2))

r=R1 :

{

q1(r,t
f)= q2(r,t

f)= q(R1, t
f)

T1(r,t
f)=T2(r,t

f)=T(R1, t
f)

(3.10)

Fig. 2. Two-layer spherical domain

We put equations (3.10) to (3.8) for m = 1,2. Additionally, taking into
account the boundary conditions for r=R0 and r=R2 we have

[

g111 g
1
12

g121 g
1
22

][

qb

q(R1, t
f)

]

=

[

h111 h
1
12

h121 h
1
22

][

T1(R0, t
f)

T(R1, t
f)

]

+

[

p11

p12

]

[

g211 g
2
12

g221 g
2
22

][

q(R1, t
f)

α[T2(R2, t
f)−T∞]

]

=

[

h211 h
2
12

h221 h
2
22

][

T(R1, t
f)

T2(R2, t
f)

]

+

(3.11)

+

[

p21

p22

]

Well-ordered systems (3.11) can be written in the form



364 B.Mochnacki, R.Szopa

[

−h111 −h112 g112
−h121 −h122 g122

]









T1(R0, t
f)

T(R1, t
f)

q(R1, t
f)









=

[

−g111qb+p11
−g121qb+p12

]

(3.12)

[

−h211 g211 αg212−h212
−h221 g221 αg222−h222

]









T(R1, t
f)

q(R1, t
f)

T2(R2, t
f)









=

[

αg212T
∞+p21

αg222T
∞+p22

]

Finally, one obtains the following system of equations

















−h111 −h112 g112 0
−h121 −h122 g122 0
0 −h211 g211 αg212−h212
0 −h221 g221 αg222−h222

































T1(R0, t
f)

T(R1, t
f)

q(R1, t
f)

T2(R2, t
f)

















=

















−g111qb+p11
−g121qb+p12
αg212T

∞+p21

αg222T
∞+p22

















at the same time

q2(R2, t
f)=α[T2(R2, t

f)−T∞] (3.13)
The knowledge of boundary values for r = R0, r = R1 and r = R2 allows
one to find the internal temperatures at the time tf using the equation (cf.
formula (3.4))

Tm(ξ,t
f)= gm(ξ,Rm−1)qm(Rm−1, t

f)−gm(ξ,Rm)qm(Rm, tf)+

+hm(ξ,Rm)Tm(Rm, t
f)−hm(ξ,Rm−1)Tm(Rm−1, tf)+pm(ξ)

The similar algorithm can be used in the case of non-zero thermal resi-
stance Z between sub-domains considered. Then the continuity condition can
be written in the form

r=R1 : −λ1
∂T1(r,t)

∂r
=
T1(r,t)−T2(r,t)

Z
=−λ2

∂T2(r,t)

∂r
(3.14)

or

r=R1 :

{

q1(r,t
f)= q2(r,t

f)= q(R1, t
f)

T2(r,t
f)=T1(r,t

f)−Zq(R1, tf)
(3.15)

It should be pointed out that for Z =0 the last condition takes form (3.10).



Numerical modelling of heat transfer... 365

The resolving system for the problem discussed can be written as follows
















−h111 −h112 g112 0
−h121 −h122 g122 0
0 −h211 g211+Zh211 αg212−h212
0 −h221 g221+Zh221 αg222−h222

































T1(R0, t
f)

T(R1, t
f)

q(R1, t
f)

T2(R2, t
f)

















=

















−g111qb+p11
−g121qb+p12
αg212T

∞+p21

αg222T
∞+p22

















(3.16)
at the same time

T2(R1, t
f)=T1(R1, t

f)−Zq(R1, tf) (3.17)

Aspreviously, the internal temperatures in successive layers canbe foundusing
equation (3.15).

4. Examplary computations

The first example (Szopa, 1999) concerns a boundary initial problem for
which a constant boundary temperature is known T(R,t) = Tb, while for
t=0 : T(r,0)= 0. The problem can be solved in the exact way (Kącki, 1992).
Test computations show, that a very good accuracy of numerical solutions can
be obtained in a wide range of time steps, and they are practically the same
as the exact result for the dimensionless time interval ∆F0 ∈ [0.001,0.009].
In Figures 3 and 4 the heating curves at selected points (r = 0.05R,

0.5R, 0.75R, 0.95R) are shown. The first solution (Fig.3) was obtained for
R=0.1m, λ=35W/(mK), c=4.875 ·106J/(m3K), Tb =100◦C,∆t=4.17s
(∆F0 = 0.003) and ∆t = 8.35s (∆F0 = 0.006), while the second solution
(Fig.4) was obtained for R= 0.2m, λ= 1W/(mK), c= 1.75 ·106J/(m3K),
Tb = 100

◦C, ∆t = 210s (∆F0 = 0.003) and ∆t = 420s (∆F0 = 0.006). In
Figures 3 and 4 the exact solution is also marked.
The very good accuracy of numerical solutionwas also obtained in the case

of the Robin boundary condition: q(R,t) = α[T(R,t)−T∞], where α is the
heat transfer coefficient and T∞ is the ambient temperature.
In Figure 5 the numerical and exact solution (symbols) for R= 0.2, λ=

35W/(mK), c=4.875·106J/(m3K),α=350W/(m2K)(Biot’s numberBi=1),
T∞ =0, initial temperature T(r,0)= 100◦C,∆F0 =0.002, 0.003, 0.006, 0.01.
The cooling curves correspond to points r=0.05R, 0.5R, 0.75R, 0.95R.
The next examples concern non-homogeneous domains. We consider a

sphere made from cast iron (R1 = 0.05m) which is spread within a steel



366 B.Mochnacki, R.Szopa

Fig. 3. Heating curves (example 1)

Fig. 4. Heating curves (example 2)



Numerical modelling of heat transfer... 367

Fig. 5. Cooling curves

spherical shell of thickness 0.01m. The initial temperature of the sphere (the
internal radius R0 = 10

−5m) is equal T10 = 20
◦C, while the initial tempe-

rature of the shell T20 = 300
◦C. The differentiation of initial temperatures

results fromtheneed of assurring a good contact between the layers at the am-
bient temperature. The following thermophysical parameters of sub-domains
are assumed: λ1 = 53W/(mK), c1 = 3.917 ·106J/(m3K), λ2 = 30W/(mK),
c2 = 4.875 · 106J/(m3K). The heat transfer coefficient on the outer surface
r = R2 : α = 30W/(m

2K), ambient temperature T∞ = 20◦C. For r = R0
the adiabatic condition is assumed. The interior is divided into 25 linear ele-
ments, time step ∆t=2.5s (Fig.6) and ∆t=1.25s (Fig.7).

The results are compared with the numerical solution obtained using a
repeatedly verified FDM program (symbols in Figures 6 and 7). It should be
pointed out that the solutions are similar – a somewhat better agreement one
obtains for the time step ∆t=1.25s.

The last example concerns the problem of heat conduction in the domain
considered for the case of non-zero thermal resistance between the sphere and
shell. It is assumed that Z = const = 0.001m2K/W. The geometry of the
domain and the values of thermophysical parameters are the same as previo-
usly. For r=R0 : q(r,t) = 0, for r=R2 :α=30W/(m

2K), T∞ =20◦C, the



368 B.Mochnacki, R.Szopa

Fig. 6. Temperature field in domain considered (∆t=2.5s)

Fig. 7. Temperature field in domain considered (∆t=1.25s)



Numerical modelling of heat transfer... 369

Fig. 8. The solution with thermal resistance (∆t=1.25s)

initial temperatures T10 = 20
◦C, T20 = 300

◦C. The results of computations
are shown in Figure 8. The continuous lines illustrate the solution obtained,
while the symbols the FDM solution. The agreement of these results is quite
satisfactory.

Summing up, the algorithms presented in this paper can be used in the
numerical modelling of heat conduction proceeding in spherical domains.

Acknowledgement

This paper is a part of research project BS-05-301 realised at the Institute of

Mathematics and Computer Sciences.

References

1. Bokota A., 1989, CPBP 02.09, Report 02.06, Częstochowa, (not published)

2. Brebbia C.A., Telles J.C.F., Wrobel L.C., 1984, Boundary Element
Techniques, Springer-Verlag, Berlin, NewYork



370 B.Mochnacki, R.Szopa

3. Curan D.A.S., Cross M., Lewis B.A., 1980, Solution of parabolic differen-
tial equations by the BEMusing discretisation in time,Appl. Math. Modelling,
4, 398-400

4. Kącki E., 1992,Partial Differential Equations in Physical and Technical Pro-
blems, WarsawWNT (in polish)

5. MajchrzakE.,MochnackiB., 1995,Application of theBEMin the thermal
theory of foundry,Engineering Analysis with Boundary Elements, 16, 99-121

6. Sichert W., 1989,Berechnung von instationaren thermishen Problemen mit-
tels der Randelement-methode, Erlangen

7. SzopaR., 1999,Modelling of Solidification andCrystallizationUsing theCom-
bined Variant of the BEM, Publ. of the Silesian Technical University, Metal-
lurgy, 54, Gliwice (in polish)

Model numeryczny przepływu ciepła w obszarach sferycznych

z wykorzystaniem kombinowanej metody elementów brzegowych

Streszczenie

Kombinowany wariant metody elementów skończonych, nazywany w literaturze
MEB, z dyskretyzacją czasu polega na zastąpieniuwystepującej w równaniuFouriera
pochodnej temperatury po czasie odpowiednim ilorazem różnicowym. Dalsze etapy
przekształceń matematycznych i konstrukcji algorytmu numerycznego nie odbiegają
od typowego podejścia charakteryzującego klasycznąmetodę elementów brzegowych.
W pracy metodę kombinowaną wykorzystano do modelowania nieustalonej dyfuzji
ciepła w obszarach sferycznych jednorodnych i niejednorodnych. W końcowej części
przedstawiono przykłady obliczeń numerycznych.

Manuscript received February 27, 2001; accepted for print September 25, 2001