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 in Fourier’s equation by an adequate differential quotient. The next steps of mathematical manipulations and also the numerical algori-thm are similar to a typical boundary element approach. In the paper 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 X e=1 ∂2T (x, t) ∂x2 e = a∇ 2 T (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 , t denote 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 ∆t must be introduced 0 = t0 < t1< t2< ... < tf −1< tf < ... < tF
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∇ 2 T (x, tf) (1.3) or ∇2T (x, tf ) − 1 a∆tT (x, t f) + 1 a∆tT (x, t f −1) = 0 (1.4)
Using the weighted residual criterion one obtains (Brebbia et al., 1984; Majchrzak and Mochnacki, 1995)
Z Ω h ∇2T (x, tf ) −a∆t1 T (x, tf) + 1 a∆tT (x, t f −1)i 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πdexp −√d a∆t d = v u u t 3 X 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) −a∆t1 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 λ Z Γ U∗ (d)q(x, tf) dΓ = (1.8) = 1 λ Z Γ Q∗ (d)T (x, tf) dΓ + 1 a∆t Z Ω U∗ (d)T (x, tf −1) dΩ where q = −λ∂T/∂n, Q∗ = −λ∂U∗ /∂n.
2. Homogeneous spherical domain
According to the idea presented by Brebbia 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 discretizequ-ation in time.
The spherical co-ordinate system should be introduced
x1 = r cos ϕ sin θ x2= r sin ϕ sin θ x3 = r cos θ (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 =
q
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Γ = R2sin θ dθdϕ dΩ = r2sin θ dθdϕdr (2.3)
Thus, equation (1.8) takes the form T (ξ, tf) +R 2 λ q(R, t f) 2π Z 0 π Z 0 U∗ (d) sin θ dθdϕ = = R 2 λ T (R, t f) 2π Z 0 π Z 0 Q∗ (d) sin θ dθdϕ + (2.4) + 1 a∆t R Z 0 r2 2π Z 0 π Z 0 U∗ (d) sin θ dθdϕT (r, tf −1) dr or T (ξ, tf) +R 2 λ T ∗ (ξ, R)q(R, tf) = (2.5) = R 2 λ q ∗ (ξ, R)T (R, tf) + 1 a∆t R Z 0 r2T∗ (ξ, r)T (r, tf −1 dr where T∗ (ξ, r) = 2π Z 0 π Z 0 U∗ (d) sin θ dθdϕ = = 1 4π 2π Z 0 π Z 0 1 p r2+ ξ2− 2rξexp − s 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ξ h exp−|r − ξ| a∆t − exp−|r + ξ| a∆t i (2.7) q∗ (ξ, r) = λ rT ∗ (ξ, r) + λ 2ξr h sgn (r − ξ) exp−|r − ξ| a∆t − − sgn (r + ξ) exp−|r + ξ|a∆t i
The boundary equation (for ξ → R) is of the form (cf. equation (2.5)) T (R, tf) + R 2 λ T ∗ (R, R)q(R, tf) = (2.8) = R 2 λ q ∗ (R, R)T (R, tf) + 1 a∆t R Z 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) = R 2 λ q ∗ (ξ, R)T (R, tf) − (2.9) −R 2 λ T ∗ (ξ, R)q(R, tf) + 1 a∆t R Z 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 h r2∂Tm(r, t) ∂r i (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)
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(ξ, tf) + hr2 λmT ∗ m(ξ, R)qm(r, tf) iRm Rm−1 = (3.3) =hr 2 λm q∗ m(ξ, r)Tm(r, tf) iRm Rm−1 + 1 am∆t Rm Z Rm−1 r2T∗ m(ξ, r)Tm(r, tf −1) dr
Equation (3.3) can be written as follows
Tm(ξ, tf) + qm(ξ, Rm)qm(Rm, tf) − qm(ξ, Rm−1)qm(Rm−1, tf) = (3.4) = hm(ξ, Rm)Tm(Rm, tf) − hm(ξ, Rm−1)Tm(Rm−1, tf) = pm(ξ) where qm(ξ, r) = r 2 λmT ∗ m(ξ, r) hm(ξ, r) = r2 λmq ∗ m(ξ, r) (3.5) while pm(ξ) = 1 am∆t Rm Z Rm−1 r2T∗ m(ξ, r)Tm(r, tf −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(Rm−1+ , Rm) gm(R− m, Rm−1) gm(R − m, Rm) # " qm(rm−1, tf) qm(Rm, tf) # = = " 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, tf) Tm(Rm, tf) # + (3.7) + " pm(Rm−1) pm(Rm) # or " gm 11 g12m gm 21 g22m # " q m(Rm−1, tf) qm(Rm, tf) # = " hm 11 hm12 hm 21 hm22 # " T m(Rm−1, tf) Tm(Rm, tf) # + " pm 1 pm 2 # (3.8)
The final system for multi-layers domain results from the 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, and we assume the boundary conditions for r = R0 and r = R2 in the form
r = R0: q1(r, tf) = qb
r = R2: q2(r, tf) = α[T2(r, tf) − 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, tf) = q2(r, tf) = q(R1, tf) T1(r, tf) = T2(r, tf) = T (R1, tf) (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
" g1 11 g121 g1 21 g 1 22 # " q b q(R1, tf) # = " h1 11 h112 h1 21 h 1 22 # " T 1(R0, tf) T (R1, tf) # + " p1 1 p1 2 # " g2 11 g122 g2 21 g222 # " q(R 1, tf) α[T2(R2, tf) − T∞] # = " h2 11 h212 h2 21 h222 # " T (R 1, tf) T2(R2, tf) # + (3.11) + " p2 1 p2 2 #
" −h1 11 −h112 g121 −h1 21 −h122 g221 # T1(R0, tf) T (R1, tf) q(R1, tf) = " −g1 11qb+ p11 −g1 21qb+ p12 # (3.12) " −h2 11 g112 αg122 − h212 −h2 21 g212 αg222 − h222 # T (R1, tf) q(R1, tf) T2(R2, tf) = " αg2 12T ∞ + p2 1 αg2 22T ∞ + p2 2 #
Finally, one obtains the following system of equations
−h1 11 −h112 g121 0 −h1 21 −h122 g221 0 0 −h2 11 g112 αg122 − h212 0 −h2 21 g212 αg222 − h222 T1(R0, tf) T (R1, tf) q(R1, tf) T2(R2, tf) = −g1 11qb+ p11 −g1 21qb+ p12 αg2 12T ∞ + p2 1 αg2 22T ∞ + p2 2
at the same time
q2(R2, tf) = α[T2(R2, tf) − 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(ξ, tf) = gm(ξ, Rm−1)qm(Rm−1, tf) − gm(ξ, Rm)qm(Rm, tf) + +hm(ξ, Rm)Tm(Rm, tf) − 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, tf) = q2(r, tf) = q(R1, tf) T2(r, tf) = T1(r, tf) − Zq(R1, tf) (3.15) It should be pointed out that for Z = 0 the last condition takes form (3.10).
The resolving system for the problem discussed can be written as follows −h1 11 −h112 g121 0 −h1 21 −h122 g221 0 0 −h2 11 g112 + Zh211 αg122 − h212 0 −h2 21 g212 + Zh221 αg222 − h222 T1(R0, tf) T (R1, tf) q(R1, tf) T2(R2, tf) = −g1 11qb+ p11 −g1 21qb+ p12 αg2 12T ∞ + p2 1 αg2 22T ∞ + p2 2 (3.16) at the same time
T2(R1, tf) = T1(R1, tf) − Zq(R1, tf) (3.17)
As previously, the internal temperatures in successive layers can be found using 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.1 m, λ = 35 W/(mK), c = 4.875 · 106J/(m3K), Tb = 100◦
C, ∆t = 4.17 s (∆F0 = 0.003) and ∆t = 8.35 s (∆F0 = 0.006), while the second solution
(Fig. 4) was obtained for R = 0.2 m, λ = 1 W/(mK), c = 1.75 · 106J/(m3K),
Tb = 100◦
C, ∆t = 210 s (∆F0 = 0.003) and ∆t = 420 s (∆F0 = 0.006). In
Figures 3 and 4 the exact solution is also marked.
The very good accuracy of numerical solution was 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, λ = 35 W/(mK), c = 4.875·106J/(m3K), α = 350 W/(m2K) (Biot’s number Bi=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.05 m) which is spread within a steel
Fig. 3. Heating curves (example 1)
Fig. 5. Cooling curves
spherical shell of thickness 0.01 m. 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 from the need of assurring a good contact between the layers at the am-bient temperature. The following thermophysical parameters of sub-domains are assumed: λ1 = 53 W/(mK), c1 = 3.917 · 106J/(m3K), λ2 = 30 W/(mK),
c2 = 4.875 · 106J/(m3K). The heat transfer coefficient on the outer surface
r = R2 : α = 30 W/(m2K), 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.5 s (Fig. 6) and ∆t = 1.25 s (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.25 s.
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.001 m2K/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 : α = 30 W/(m2K), T∞= 20◦C, the
Fig. 6. Temperature field in domain considered (∆t = 2.5 s)
Fig. 8. The solution with thermal resistance (∆t = 1.25 s)
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
3. Curan D.A.S., Cross M., Lewis B.A., 1980, Solution of parabolic differen-tial equations by the BEM using discretisation in time, Appl. Math. Modelling,
4, 398-400
4. Kącki E., 1992, Partial Differential Equations in Physical and Technical Pro-blems, Warsaw WNT (in polish)
5. Majchrzak E., Mochnacki B., 1995, Application of the BEM in 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. Szopa R., 1999, Modelling of Solidification and Crystallization Using the Com-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ąpieniu wystepującej w równaniu Fouriera 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.
