- Email: [email protected]
An analytical method to solve heat conduction in layered spheres with time-dependent boundary conditions
Physics Letters A 351 (2006) 274–282 www.elsevier.com/locate/pla
An analytical method to solve heat conduction in layered spheres with time-dependent boundary conditions Xiaoshu Lu a,b,∗ , Martti Viljanen a a Laboratory of Structural Engineering and Building Physics, Department of Civil and Environmental Engineering, Helsinki University of Technology,
PL 2100, FIN-02015 Hut, Espoo, Finland b Department of Physiology, Finnish Institute of Occupational Health, Topeliuksenkatu 41 a A, FIN-00250 Helsinki, Finland
Received 30 June 2005; received in revised form 28 October 2005; accepted 7 November 2005 Available online 17 November 2005 Communicated by A.P. Fordy
Abstract An analytical method is proposed to solve the equation of heat conduction in a layered sphere subject to a time-dependent boundary temperature. It is well known that for such problems in general, eigenvalue and residue computation poses a challenge, which can become too complicated to handle with many layers. In this Letter, the proposed analytical method is free of eigenvalue and residue calculations. A closed-form approximate solution is derived with high accuracy. 2005 Elsevier B.V. All rights reserved. PACS: 44.10.+i; 44.05.+e Keywords: Layered sphere; Heat conduction; Analytical method; Laplace transform
1. Introduction The study of transient heat conduction in a composite slab is important in physics and engineering areas such as thermodynamics, fuel cells, electrochemical reactors, high density microelectronics, food products, solidification processes and many others. In terms of the geometry of the slab, the layered sphere is extensively used in studying the thermal properties of composite media by assuming that spherical particles are embedded in the composite matrix. Thermal analysis of such a composite matrix has wide applications in manufacturing also. The existing studies are often based on the analytical analysis of heat conduction in a layered sphere. The purpose of this article is to establish a closed-form approximate analytical solution for transient heat conduction in a layered sphere subject to a time-dependent boundary condition. Analytical methods for the transient heat conduction problem in the spherical geometry are completely analogous to those in Cartesian coordinates. The commonly applied techniques are finite integral transforms which are often employed to a single layer slab, Green’s functions, orthogonal expansions and Laplace’s transform [1]. A traditional application of the first two techniques often leads to eigenvalue problems. In a single layer slab, for instance, the eigenfunction links the space and time variables when applying separation of variables. However, in a multilayer slab, the eigenfunctions can present boundary conditions at the contacted layers. Hence, eigenvalue problems may exist even for the steadystate heat conduction problem in a composite slab (e.g. [2]). The associated eigenvalue calculations may become much more complicated.
* Corresponding author.
E-mail address: [email protected] (X. Lu). 0375-9601/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2005.11.017
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
275
Similarly, an application of the third technique, the Laplace transform, often yields a residue computation. In a composite slab, the residue computation is found by directly and numerically searching for the roots of a hyperbolic equation, finding the derivatives of the equation, and evaluating and summing the residues. The calculation procedure is tedious if the slab has more than two layers, as numerical searching for roots has to be made with a very fine increment for the inverse Laplace transform to prevent missing roots which can lead to a wrong inverse [3]. Due to the complex nature of the associated eigenvalue and residue problems, the advantage of analytical methods over numerical methods is sometimes hard to recognise. Closed-form solutions are lacking for heat conduction problems in composite slabs though there has been extensive research on this topic. Most of the analytical solutions in literature are available formally: numerical programs for searching eigenvalues and residues are necessary. In Cartesian coordinates, a closed-form solution was given for the transient heat conduction problem in a three-layer composite slab subject to an unchanged temperature [4]. A closed-form solution in an n-layer composite slab was conjectured [4]. In cylindrical coordinates, Fredman developed a semi-analytical method for composite layer diffusion [5]. The method relied on a local analytical solution for a single material layer combined with a numerical solution scheme for the material boundary states [5]. In a spherical geometry, an analytical–numerical technique was employed to predict the temperature and heat flux histories in a bi-layer composite sphere [6]. In a multi-dimensional composite slab, similar analytical techniques are often adopted for the heat conduction problem (e.g. [7] with Green’s function). The associated eigenvalue problems are much more complicated as the eigenvalues may become imaginary. Attention must be paid when computing eigenvalues since spacing between successive eigenvalues changes between zero and the maximal value. Numerically, imaginary eigenvalues can produce instability [7]. For the steady-state heat conduction problem, HajiSheikh et al. presented an analytical solution in a three-dimensional two-layer slab [8]. In both papers, an effort has been made to develop an efficient algorithm to compute the eigenvalues. Most of the above-cited papers produce an exact solution for the heat conduction problem. There is no doubt that the exact solutions are extremely valuable whenever they appear. However, all these works need numerical schemes for eigenvalues which are tedious if the slab has more than three layers. Recently, a novel analytical method was developed to solve the one-dimensional transient heat problem in a composite slab with a time-dependent boundary temperature [9]. Unlike most of the traditional methods, the new method involves no numerical computation such as a numerical search for eigenvalues and residues. The boundary conditions are given more generally. An approximate closed-form solution was provided with high accuracy. In this Letter, the developed analytical method is extended to a spherical structure. The objective of this study is to derive a more general closed-form solution for the transient heat conduction in a layered sphere subject to a time-dependent temperature change. It is believed that the results can be beneficial to the study of thermal properties of composite media. 2. Mathematical formulation 2.1. Governing equations Consider an n-layer sphere having constant thermal conductivity, diffusivity and radius in each layer which are denoted by λj , kj and rj , j = 1, . . . , n. A schematic figure is shown in Fig. 1. The spherical shells are [rj −1 , rj ], j = 1, . . . , n, with r0 = 0. The general heat conduction problem for the composite sphere can be described by the following equation for temperatures Tj (t, r) in spherical coordinates: kj ∂ ∂Tj ∂Tj = 2 r2 , r ∈ [rj −1 , rj ], j = 1, . . . , n. (2.1a) ∂t ∂r r ∂r The boundary conditions are given as Tj (t, r) = Tj +1 (t, r),
r = rj , j = 1, . . . , n − 1,
Fig. 1. Schematic diagram of an n-layer composite sphere.
(2.1b)
276
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
∂Tj ∂Tj +1 (t, r) = −λj +1 (t, r), r = rj , j = 1, . . . , n − 1, ∂r ∂r ∂Tn (t, r) = −α T∞ (t) − Tn (t, r) , r = rn , −λn ∂r and the initial condition is −λj
Tj (0, r) = 0,
r ∈ [rj −1 , rj ], j = 1, . . . , n.
(2.1c) (2.1d)
(2.1e)
A perfect thermal contact is assumed at the contact interfaces. Note that r1 = 0 cannot be defined in Eq. (2.1a). For simplicity, it is assumed that the outer boundary temperature varies periodically T∞ (t) = cos(ωt + ϕ). Cases with more general boundary temperatures will be discussed later. The surface heat transfer coefficient is α. The initial temperature is set to zero for the sake of calculation convenience. 2.2. Dimensionless parameters Writing Uj = rTj ,
j = 1, . . . , n,
(2.2)
leads to the following constant-coefficient heat conduction equation for Eq. (2.1a): ∂Uj ∂ 2 Uj , r ∈ [rj −1 , rj ], j = 1, . . . , n. = kj ∂t ∂r 2 The boundary and initial conditions (2.1b)–(2.1e) are
(2.3a)
Uj = Uj +1 , r = rj , j = 1, . . . , n − 1, Uj 1 ∂Uj 1 ∂Uj +1 Uj +1 λj − 2 = λj +1 − 2 , r = rj , j = 1, . . . , n − 1, r ∂r r ∂r r r Un 1 ∂Un Un λn − 2 = α T∞ (t) − , r = rn , r ∂r r r Uj (0, r) = 0, r ∈ [rj −1 , rj ], j = 1, . . . , n,
(2.3b)
T∞ (t) = cos(ωt + ϕ).
(2.3f)
(2.3c) (2.3d) (2.3e)
It should be noticed that, because of the form of the boundary conditions (2.3b)–(2.3e), the calculated results for a composite sphere do not follow immediately from the corresponding results for a composite slab [1]. First, we proceed to the analytical solution by assuming a complex form of the boundary temperature T∞ (t) = exp(iωt + iϕ) = cos(ωt + ϕ) + i sin(ωt + ϕ). Clearly, the real part of the solution is the sought-after solution. We shall adopt the same notation for the corresponding solution of Eq. (2.1). Furthermore, for calculation convenience, denote U∞ (t) = rn T∞ (t). Define the following dimensionless parameters: t∗ =
tkn , rn2
r∗ =
r , rn
kj∗ =
kj , kn
λ∗j =
λj , λn
α∗ =
αrn . λn
(2.4)
Then the spherical shells are changed as [rj∗−1 = rj −1 /rn , rj∗ = rj /rn ], j = 1, . . . , n. Eq. (2.3a) becomes ∂Uj ∂ 2 Uj = kj∗ ∗ 2 , r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n. ∗ ∂t ∂r The boundary and initial conditions (2.3b)–(2.3f) are Uj = Uj +1 , r ∗ = rj∗ , j = 1, . . . , n − 1, Uj 1 ∂Uj +1 Uj +1 ∗ 1 ∂Uj ∗ λj ∗ ∗ − ∗ 2 = λj +1 ∗ − ∗ 2 , r ∗ = rj∗ , j = 1, . . . , n − 1, r ∂r r ∂r ∗ r r Un Un 1 ∂Un λ∗n ∗ ∗ − ∗ 2 = α ∗ U∞ t ∗ − ∗ , r ∗ = 1, r ∂r r r ∗ ∗ ∗ ∗ Uj 0, r = 0, r ∈ rj −1 , rj , j = 1, . . . , n, U∞ t ∗ = rn exp iω∗ t ∗ + iϕ ,
(2.5a)
(2.5b) (2.5c) (2.5d) (2.5e) (2.5f)
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
277
where ω∗ =
ωrn2 . kn
(2.5g)
3. Solution method Performing a Laplace transformation on Eq. (2.5) yields ∂ 2 U¯ j , r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n. ∂r ∗ 2 The boundary conditions are s U¯ j = kj∗
U¯ j = U¯ j +1 , ¯ ∗ 1 ∂ Uj λj ∗ ∗ − r ∂r 1 ∂ U¯ n λ∗n ∗ ∗ − r ∂r
(3.1a)
r ∗ = rj∗ , j = 1, . . . , n − 1, U¯ j 1 ∂ U¯ j +1 U¯ j +1 ∗ − ∗ 2 , r ∗ = rj∗ , j = 1, . . . , n − 1, = λj +1 ∗ r ∂r ∗ r∗ 2 r U¯ n U¯ n = α ∗ U¯ ∞ (s) − ∗ , r ∗ = 1. ∗ 2 r r
(3.1b) (3.1c) (3.1d)
The Laplace transform of f (t ∗ ) is defined as [1]: f¯(s) =
∞
exp −st ∗ f t ∗ dt ∗ .
(3.2a)
0
The Laplace transform of a convolution is given by [1]: f1 ∗ f2 (s) = f¯1 (s)f¯2 (s) ⇔ f1 t ∗ ∗ f2 t ∗ =
t ∗
f1 (τ )f2 t ∗ − τ dτ.
(3.2b)
0
The solution of Eq. (3.1a) is easily constructed as U¯ j = Aj sinh qj r ∗ − rj∗−1 + Bj cosh qj r ∗ − rj∗−1 , r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n, (3.3) where qj = s/kj∗ and coefficients Aj and Bj are determined by applying the boundary conditions (3.1b)–(3.1d) as follows: B1 = 0,
(3.4a)
Aj sinh ξj + Bj cosh ξj − Bj +1 = 0,
j = 1, . . . , n − 1,
(3.4b)
Aj (cosh ξj − ηj sinh ξj ) + Bj (sinh ξj − ηj cosh ξj ) − hj Aj +1 +
λj +1 ηj Bj +1 = 0, λj
j = 1, . . . , n − 1,
(3.4c)
An hA + Bn hB = α ∗ U¯ ∞ ,
(3.4d)
where
∗ ∗
kj λ 1 j +1 ∗ ∗ ηj = for j = 1, . . . , n, hj = ∗ ∗ for j = 1, . . . , n − 1, ξj = qj rj − rj −1 , ∗ qj rj λj kj +1 ∗ ∗ α α hA = qn cosh ξn + − ηn sinh ξn , − ηn cosh ξn . hB = qn sinh ξn + qn qn
The coefficients Aj and Bj can be solved from Eq. (3.4) by applying Cramer’s rule. Let
0 sinh ξ1 cosh ξ1 − η1 sinh ξ1 ··· ∆(s) = 0 0 0
λ2 η1 λ1
··· 0
··· ··· ··· ··· ···
0 0 0 ··· sinh ξn−1
0 0 0 ··· cosh ξn−1
0 0 0 ··· 0
0
0
···
cosh ξn−1 − ηn−1 sinh ξn−1
sinh ξn−1 − ηn−1 cosh ξn−1
−hn−1
0
0
···
λn ηn−1 λn−1
0
0
hA
hB
1 cosh ξ1 sinh ξ1 − η1 cosh ξ1 ··· 0
0 0 −h1 ··· 0
0 0
0 −1
0 0 0 ··· −1
(3.5a)
278
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
↓ 0 ··· 0 α ∗ U¯ ∞
Column 2j −1
↓ 0 ··· 0 α ∗ U¯ ∞
Column 2j
··· ··· ··· ··· ··· ··· · · · ··· , , ∆ ∆2j −1 (s) = (s) = 2j ∆(s) ∆(s) ∆(s) ∆(s) ··· ··· ··· ··· det(∆2j −1 (s)) det(∆2j (s)) Aj = Bj = = δ2j −1 (s)U¯ ∞ , = δ2j (s)U¯ ∞ , det(∆(s)) det(∆(s)) where δ2j −1 (s) and δ2j (s) are the matrix determinants obtained from Eqs. (3.5b), (3.5c). Finally, from Eq. (3.3) we get U¯ j = Aj sinh qj r ∗ − rj∗−1 + Bj cosh qj r ∗ − rj∗−1 = Gj s, r ∗ U¯ ∞ , where
Gj s, r ∗ = δ2j −1 (s) sinh qj r ∗ − rj∗−1 + δ2j (s) cosh qj r ∗ − rj∗−1 ,
r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n.
(3.5b)
(3.5c)
(3.6a)
(3.6b)
In deducing the inverse Laplace transform of U¯ j , we shall make an approximation to avoid evaluating residues. Let gj (t ∗ , r ∗ ) be the inverse Laplace transform of Gj (s, r ∗ ). Eq. (3.6a) is then expressed as (see Eq. (3.2b))
∗
Uj t , r
∗
∗
= gj t , r
∗
∗ U∞ t ∗ =
t ∗
∗
gj τ, r U∞
∞ ∞ t − τ dτ = − gj τ, r ∗ U∞ t ∗ − τ dτ
∗
0
∞ ≈
0
gj τ, r ∗ U∞ t ∗ − τ dτ =
0
∞
t∗
gj τ, r ∗ rn exp iω∗ t ∗ − τ + iϕ dτ
0
= rn exp iω∗ t ∗ + iϕ
∞
exp −iω∗ τ gj τ, r ∗ dτ = rn Gj iω∗ , r ∗ exp iω∗ t ∗ + iϕ = Gj iω∗ , r ∗ U∞ t ∗ .
0
(3.7)
↑
Laplace transform of
gj
at
iω∗ = Gj (iω∗ ,r ∗ ),
see Eq. (3.2a)
Note that the inverse Laplace transform gj (t ∗ , r ∗ ) was acting only as a symbolic function. By taking an advantage of the exponential expression of U∞ (t ∗ ), gj (t ∗ , r ∗ ) was replaced by its Laplace transform Gj (iω, r) which is already available in Eq. (3.6b). 4. Solution method for general boundary conditions For a general time-dependent boundary condition, we approximate it using Fourier series as T∞ (t) = a0 +
∞
ak cos(ωk t + ϕk ).
(4.1)
k=1
Hence, the solution can be split up into two subproblems, namely, Uj = Uj1 + Uj2 : (1) solution Uj1 : boundary condition is constant, a0 ; (2) solution Uj2 : boundary condition is time-dependent ∞ k=1 ak cos(ωk t + ϕk ). The second subproblem is equivalent to that with the boundary condition cos(ωt + ϕ) due to the linearity of the system. For the first subproblem with the constant boundary condition, U∞ = a0 , Eq. (3.6a) gives a 0 a0 U¯ j1 = Gj s ∗ , r ∗ U¯ ∞ = Gj s ∗ , r ∗ (4.2) = δ2j −1 (s) sinh qj r ∗ − rj∗−1 + δ2j (s) cosh qj r ∗ − rj∗−1 . s s The matrix determinants δ2j −1 (s) and δ2j (s) are functions of hyperbolic functions sinh ξ and cosh ξ which can be approximated by power series such as sinh ξ = ξ/1! + ξ 3 /3! + · · · = U¯ j1 ≈ Finally, Uj1 =
s/kj∗ (rj∗ − rj∗−1 ) + · · · . Hence, linearisation of Eq. (4.2) gives
const . const 1 ∗ s + const 2
(4.3)
const const 2 ∗ exp − t . const 1 const 1
(4.4)
Note that the series may ascend as powers of
√
s. An approximate inverse of U¯ j1 is then obtained in the same way as that in Eq. (4.3).
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
279
5. Calculation examples Since an approximate formula has been used in reaching the mathematical solution, this section is aimed at assessing the accuracy of the formula. A selected five-layer sphere is presented here as a calculation example. Its shells are numbered by layer one to five from the inner core. A schematic figure is shown in Fig. 2. Thermal properties and dimensions of the composite are listed in Table 1. The composition is an extension of the three-layer slab employed as an exterior wall structure in our test house to evaluate the accuracy of our model for predicting the evolution of the temperature distributions inside the wall subject to interior and exterior climatic conditions. The interior and exterior climatic temperatures were measured. In this calculation, the boundary conditions were taken from such measurements and then fitted with cosine functions as 5 2πt ai cos − ϕi , T∞ (t) = a0 + (5.1) ωi 1
where the fitting parameters are listed in Table 2. The surface heat transfer coefficient, α, is assumed to be 25 W m−2 K−1 . In order to evaluate the accuracy of the analytical solution, calculations were made with both the analytical solution presented in Eq. (3.7) and a numerical model. In the numerical model, the finite volume method was used to do the spatial discretisation procedure for Eq. (2.1). The Crank–Nicolson time marching procedure was then applied to the discretised equations. The final algebraic equations were solved using Newton’s iterative method [10]. In the calculation, the central regions of layers 1 and 3 were chosen as the calculation points which are marked in Fig. 2. Figs. 3, 4 show the comparisons of the transient temperature variations at the calculation points using both analytical and numerical methods. The temperatures were stored in files as hourly values and shown in figures as hourly and daily values. At both calculation points, the maximal discrepancy between the analytical and numerical values is about 0.37 ◦ C at small t with a relative error of 6%. The comparison results for small values of time t are displayed in Fig. 5. To demonstrate the effect of heat capacity (ability to store heat) of the materials, the thicknesses of both layers 2 and 3 were extended to 1.0 m. The boundary conditions were kept the same (Eq. (5.1)). The calculated temperature variations in the central
Fig. 2. Schematic diagram of the five-layer sphere. Table 1 Material properties of the layered sphere, see Fig. 2 Layer 1 2 3 4 5
Thermal conductivity (W m−1 K−1 )
Thermal diffusivity (m2 s−1 )
Shell thickness (mm)
0.12 0.0337 0.9 0.147 0.23
1.5 × 10−7
25 200 100 100 13
1.47 × 10−6 3.75 × 10−7 1.61 × 10−7 4.11 × 10−7
Table 2 Parameters for Eq. (5.1)
a0 17.29892
ω1 120.0
ω2 30.0
ω3 10.0
ω4 5.0
ω5 1.0
ϕ1 28.56693
ϕ2 7.790886
ϕ3 6.483743
ϕ4 3.647276
ϕ5 3.295951
a1 2.3712
a2 −1.77464
a3 −0.53307
a4 −0.05364
a5 −0.12107
280
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
Fig. 3. Comparison of analytical and numerical results in the central region of layer 1.
Fig. 4. Comparison of analytical and numerical results in the central region of layer 3.
Fig. 5. Comparison of analytical and numerical results in the central region of layer 1.
Fig. 6. Comparison of analytical and numerical results in the central region of layer 1. Thicknesses of layers 2 and 3 are extended to 1 m.
region of layer 1 with numerical and analytical models are displayed in Fig. 6. As expected, fluctuating temperature change was obtained in the material with a larger heat capacity (a smaller thermal diffusivity). In this calculation, the maximal discrepancy between the numerical and analytical values is about 0.2 ◦ C (relative error less than 4%). These two examples demonstrate a high accuracy of the developed analytical solution.
6. Discussion
The solution of transient heat conduction for an n-layer sphere is explicitly expressed. We make the following observations. • The calculation includes only simple computation of a matrix determinant which can be easily accomplished by commercial mathematical packages like M APLE, M ATLAB and M ATHEMATICA. For any layer, only five sparse matrices are involved. The calculation load is small and the computing time is short. Furthermore, an increase of a slab’s layer applies only to the dimension of the matrices. The computing load will not be affected much. • The linearisation of G(s, r ∗ ) can be tedious in Eq. (4.2) to obtain the approximately ‘transient’ change for the average temperature if the layer number is more than three. Some mathematical packages, say M APLE, may be helpful which can handle symbolic computations. In authors’ experience, a steady-state solution as an approximate one for the average temperature is quite sufficient if studies do not focus very much on the initial temperature change. • Compared with numerical methods, the developed method is easier to implement and a possible instability in the numerical solution is avoided. The accuracy is good.
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
281
• The main contribution of this Letter, which distinguishes it from all previous work, is an elegant way to perform the inverse Laplace transformation of U¯ j to obtain the solution Uj in Eq. (3.7). The usual way of doing the inverse is to find the inverse Laplace transformation gj (t ∗ , r ∗ ) of the function Gj (s, r ∗ ) by the residue theorem. Once gj (t ∗ , r ∗ ) is available, the exact analytical solution is obtained. It is easy to see that, if the composition has more than two layers, then it is difficult to find the inverse Laplace transformation of the complicated Gj (s, r ∗ ) by the residue theorem. Furthermore, there are many other questions have to be considered when trying to numerically search for roots of the inverse. That is the main reason why a closed-form solution was given for the transient heat conduction problem in a three-layer composition only to date [4], though there has been extensive research on such topic. In this Letter in the inverse calculation of U¯ j , gj (t ∗ , r ∗ ), was acting only as a symbolic function. The real involvement in the calculation is the Laplace transformation Gj (s, r ∗ ) of gj (t ∗ , r ∗ ) which is explicitly expressed in Eq. (3.6). This achievement results in a closed-formed approximate solution but the exact solution. • There are difficulties in making an evaluation of the accuracy of the proposed method mainly dueto the complicated function ∞ of Gj (s, r ∗ ) which makes it very hard to get the inverse gj (t ∗ , r ∗ ). Therefore, the neglected term t ∗ gj (τ, r ∗ )U∞ (t ∗ − τ ) dτ cannot be estimated with various functions of U∞ . However, one thing is for sure that the neglected term approaches to zero if t is big enough. Therefore, a conservative interpretation of the accuracy of the proposed method is that the closed-form analytical solution is exact for quasi-steady-state heat conduction problems. For solving such problems in general, eigenvalue problems often exist for multilayer compositions (e.g. [2]). The proposed method is free of eigenvalue problems. • As we mentioned, the proposed method relies on an approximation in Eqs. (3.7) and (4.3). For small times t , the accuracy may be low which is probably due to a poor approximate in Eq. (4.3). We have made an error estimate of Eq. (3.7) for a simple heat conduction equation in which the exact solution is available. The result suggested that the magnitude of the omitted term is small even for small values of t . Therefore, to increase the accuracy, the analytical method for the transient heat conduction problem in a composite slab with a constant boundary temperature change in [4] can be employed. Then a numerical program is needed for the eigenvalues. • It is known that any periodic function can be represented as its Fourier series. Non-periodic function can be approximated as Fourier series in the extended interval. For both periodic and non-periodic functions, transient heat conduction must be solved during a short time interval compared with its period. In fact, the steady-state may not be reached during several periods or during normal calculational cycles. Therefore, the results obtained in this Letter are transient solutions for general boundary conditions. As demonstrated in the calculation examples, the boundary temperatures are not necessarily periodic. In Fig. 4, the boundary temperature is non-periodic function, even though the approximation involves a sum of cosine functions with a maximal period of 30 days. A transient calculation was made for 20 days which is much shorter than the duration of one period. • Finally, it is easy to see that the assumption of perfect thermal contact between layers can be removed without adding any content into the calculation method. It is only for demonstrational convenience that we prescribed this limitation. 7. Conclusions In this Letter, an analytical approach to heat conduction in a layered sphere subject to a time-dependent temperature change has been presented. The boundary temperature was approximated using Fourier series. The basic technique in obtaining the analytical solution is the Laplace transform which often yields a residue computation in the exact solution. In this approximate method, a novel mathematical computation was made in deducing an approximate inverse Laplace transform. As a result, it is free of the residue evaluation and the developed solution is an approximate. The method is shown to have considerable potential for solving the equations of heat conduction. Compared to analytical solutions that are mainly possible with some assumptions, the developed method has no restrictions and is flexible with regard to the boundary condition. Moreover, the benefit of the results is that its simple and concise mathematical form of the solutions can be used to analyse physical properties in combination with material properties in the heat transfer process. Agreement with the numerical solutions is good. In a general conduction or diffusion application context, however, numerical schemes are usually necessary and may be very complex due to the special geometry of the material. Finally, it is worth mentioning that the analysis for layered spheres cannot be extended straightforwardly to layered cylinders. Acknowledgements The authors would like to thank the Academy of Finland for supporting this research, and the referees for their valuable criticisms of the early manuscript. References [1] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Oxford Univ. Press, Oxford, UK, 1959. [2] P. Snabre, A. Perez, A. Arhaliass, Eur. Phys. J: Appl. Phys. 1 (1998) 315.
282
[3] [4] [5] [6] [7] [8] [9] [10]
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
M.C.B. Gough, Modelling heat flow in buildings: an eigenfunction approach, Ph.D. Thesis, University of Cambridge Press, 1982. Y. Sun, I.S. Wichman, Int. J. Heat Mass Transfer 47 (2004) 1555. T.P. Fredman, Heat Mass Transfer 39 (2003) 285. C. Tsai, C. Hung, Int. J. Heat Mass Transfer 46 (2003) 5137. A. Haji-Sheikh, J.V. Beck, Int. J. Heat Mass Transfer 45 (2002) 1865. A. Haji-Sheikh, J.V. Beck, D. Agonafer, Int. J. Heat Mass Transfer 46 (2003) 2363. X. Lu, P. Tervola, J. Phys. A: Math. Gen. 38 (2005) 81. S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980.
An analytical method to solve heat conduction in layered spheres with time-dependent boundary conditions Xiaoshu Lu a,b,∗ , Martti Viljanen a a Laboratory of Structural Engineering and Building Physics, Department of Civil and Environmental Engineering, Helsinki University of Technology,
PL 2100, FIN-02015 Hut, Espoo, Finland b Department of Physiology, Finnish Institute of Occupational Health, Topeliuksenkatu 41 a A, FIN-00250 Helsinki, Finland
Received 30 June 2005; received in revised form 28 October 2005; accepted 7 November 2005 Available online 17 November 2005 Communicated by A.P. Fordy
Abstract An analytical method is proposed to solve the equation of heat conduction in a layered sphere subject to a time-dependent boundary temperature. It is well known that for such problems in general, eigenvalue and residue computation poses a challenge, which can become too complicated to handle with many layers. In this Letter, the proposed analytical method is free of eigenvalue and residue calculations. A closed-form approximate solution is derived with high accuracy. 2005 Elsevier B.V. All rights reserved. PACS: 44.10.+i; 44.05.+e Keywords: Layered sphere; Heat conduction; Analytical method; Laplace transform
1. Introduction The study of transient heat conduction in a composite slab is important in physics and engineering areas such as thermodynamics, fuel cells, electrochemical reactors, high density microelectronics, food products, solidification processes and many others. In terms of the geometry of the slab, the layered sphere is extensively used in studying the thermal properties of composite media by assuming that spherical particles are embedded in the composite matrix. Thermal analysis of such a composite matrix has wide applications in manufacturing also. The existing studies are often based on the analytical analysis of heat conduction in a layered sphere. The purpose of this article is to establish a closed-form approximate analytical solution for transient heat conduction in a layered sphere subject to a time-dependent boundary condition. Analytical methods for the transient heat conduction problem in the spherical geometry are completely analogous to those in Cartesian coordinates. The commonly applied techniques are finite integral transforms which are often employed to a single layer slab, Green’s functions, orthogonal expansions and Laplace’s transform [1]. A traditional application of the first two techniques often leads to eigenvalue problems. In a single layer slab, for instance, the eigenfunction links the space and time variables when applying separation of variables. However, in a multilayer slab, the eigenfunctions can present boundary conditions at the contacted layers. Hence, eigenvalue problems may exist even for the steadystate heat conduction problem in a composite slab (e.g. [2]). The associated eigenvalue calculations may become much more complicated.
* Corresponding author.
E-mail address: [email protected] (X. Lu). 0375-9601/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2005.11.017
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
275
Similarly, an application of the third technique, the Laplace transform, often yields a residue computation. In a composite slab, the residue computation is found by directly and numerically searching for the roots of a hyperbolic equation, finding the derivatives of the equation, and evaluating and summing the residues. The calculation procedure is tedious if the slab has more than two layers, as numerical searching for roots has to be made with a very fine increment for the inverse Laplace transform to prevent missing roots which can lead to a wrong inverse [3]. Due to the complex nature of the associated eigenvalue and residue problems, the advantage of analytical methods over numerical methods is sometimes hard to recognise. Closed-form solutions are lacking for heat conduction problems in composite slabs though there has been extensive research on this topic. Most of the analytical solutions in literature are available formally: numerical programs for searching eigenvalues and residues are necessary. In Cartesian coordinates, a closed-form solution was given for the transient heat conduction problem in a three-layer composite slab subject to an unchanged temperature [4]. A closed-form solution in an n-layer composite slab was conjectured [4]. In cylindrical coordinates, Fredman developed a semi-analytical method for composite layer diffusion [5]. The method relied on a local analytical solution for a single material layer combined with a numerical solution scheme for the material boundary states [5]. In a spherical geometry, an analytical–numerical technique was employed to predict the temperature and heat flux histories in a bi-layer composite sphere [6]. In a multi-dimensional composite slab, similar analytical techniques are often adopted for the heat conduction problem (e.g. [7] with Green’s function). The associated eigenvalue problems are much more complicated as the eigenvalues may become imaginary. Attention must be paid when computing eigenvalues since spacing between successive eigenvalues changes between zero and the maximal value. Numerically, imaginary eigenvalues can produce instability [7]. For the steady-state heat conduction problem, HajiSheikh et al. presented an analytical solution in a three-dimensional two-layer slab [8]. In both papers, an effort has been made to develop an efficient algorithm to compute the eigenvalues. Most of the above-cited papers produce an exact solution for the heat conduction problem. There is no doubt that the exact solutions are extremely valuable whenever they appear. However, all these works need numerical schemes for eigenvalues which are tedious if the slab has more than three layers. Recently, a novel analytical method was developed to solve the one-dimensional transient heat problem in a composite slab with a time-dependent boundary temperature [9]. Unlike most of the traditional methods, the new method involves no numerical computation such as a numerical search for eigenvalues and residues. The boundary conditions are given more generally. An approximate closed-form solution was provided with high accuracy. In this Letter, the developed analytical method is extended to a spherical structure. The objective of this study is to derive a more general closed-form solution for the transient heat conduction in a layered sphere subject to a time-dependent temperature change. It is believed that the results can be beneficial to the study of thermal properties of composite media. 2. Mathematical formulation 2.1. Governing equations Consider an n-layer sphere having constant thermal conductivity, diffusivity and radius in each layer which are denoted by λj , kj and rj , j = 1, . . . , n. A schematic figure is shown in Fig. 1. The spherical shells are [rj −1 , rj ], j = 1, . . . , n, with r0 = 0. The general heat conduction problem for the composite sphere can be described by the following equation for temperatures Tj (t, r) in spherical coordinates: kj ∂ ∂Tj ∂Tj = 2 r2 , r ∈ [rj −1 , rj ], j = 1, . . . , n. (2.1a) ∂t ∂r r ∂r The boundary conditions are given as Tj (t, r) = Tj +1 (t, r),
r = rj , j = 1, . . . , n − 1,
Fig. 1. Schematic diagram of an n-layer composite sphere.
(2.1b)
276
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
∂Tj ∂Tj +1 (t, r) = −λj +1 (t, r), r = rj , j = 1, . . . , n − 1, ∂r ∂r ∂Tn (t, r) = −α T∞ (t) − Tn (t, r) , r = rn , −λn ∂r and the initial condition is −λj
Tj (0, r) = 0,
r ∈ [rj −1 , rj ], j = 1, . . . , n.
(2.1c) (2.1d)
(2.1e)
A perfect thermal contact is assumed at the contact interfaces. Note that r1 = 0 cannot be defined in Eq. (2.1a). For simplicity, it is assumed that the outer boundary temperature varies periodically T∞ (t) = cos(ωt + ϕ). Cases with more general boundary temperatures will be discussed later. The surface heat transfer coefficient is α. The initial temperature is set to zero for the sake of calculation convenience. 2.2. Dimensionless parameters Writing Uj = rTj ,
j = 1, . . . , n,
(2.2)
leads to the following constant-coefficient heat conduction equation for Eq. (2.1a): ∂Uj ∂ 2 Uj , r ∈ [rj −1 , rj ], j = 1, . . . , n. = kj ∂t ∂r 2 The boundary and initial conditions (2.1b)–(2.1e) are
(2.3a)
Uj = Uj +1 , r = rj , j = 1, . . . , n − 1, Uj 1 ∂Uj 1 ∂Uj +1 Uj +1 λj − 2 = λj +1 − 2 , r = rj , j = 1, . . . , n − 1, r ∂r r ∂r r r Un 1 ∂Un Un λn − 2 = α T∞ (t) − , r = rn , r ∂r r r Uj (0, r) = 0, r ∈ [rj −1 , rj ], j = 1, . . . , n,
(2.3b)
T∞ (t) = cos(ωt + ϕ).
(2.3f)
(2.3c) (2.3d) (2.3e)
It should be noticed that, because of the form of the boundary conditions (2.3b)–(2.3e), the calculated results for a composite sphere do not follow immediately from the corresponding results for a composite slab [1]. First, we proceed to the analytical solution by assuming a complex form of the boundary temperature T∞ (t) = exp(iωt + iϕ) = cos(ωt + ϕ) + i sin(ωt + ϕ). Clearly, the real part of the solution is the sought-after solution. We shall adopt the same notation for the corresponding solution of Eq. (2.1). Furthermore, for calculation convenience, denote U∞ (t) = rn T∞ (t). Define the following dimensionless parameters: t∗ =
tkn , rn2
r∗ =
r , rn
kj∗ =
kj , kn
λ∗j =
λj , λn
α∗ =
αrn . λn
(2.4)
Then the spherical shells are changed as [rj∗−1 = rj −1 /rn , rj∗ = rj /rn ], j = 1, . . . , n. Eq. (2.3a) becomes ∂Uj ∂ 2 Uj = kj∗ ∗ 2 , r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n. ∗ ∂t ∂r The boundary and initial conditions (2.3b)–(2.3f) are Uj = Uj +1 , r ∗ = rj∗ , j = 1, . . . , n − 1, Uj 1 ∂Uj +1 Uj +1 ∗ 1 ∂Uj ∗ λj ∗ ∗ − ∗ 2 = λj +1 ∗ − ∗ 2 , r ∗ = rj∗ , j = 1, . . . , n − 1, r ∂r r ∂r ∗ r r Un Un 1 ∂Un λ∗n ∗ ∗ − ∗ 2 = α ∗ U∞ t ∗ − ∗ , r ∗ = 1, r ∂r r r ∗ ∗ ∗ ∗ Uj 0, r = 0, r ∈ rj −1 , rj , j = 1, . . . , n, U∞ t ∗ = rn exp iω∗ t ∗ + iϕ ,
(2.5a)
(2.5b) (2.5c) (2.5d) (2.5e) (2.5f)
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
277
where ω∗ =
ωrn2 . kn
(2.5g)
3. Solution method Performing a Laplace transformation on Eq. (2.5) yields ∂ 2 U¯ j , r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n. ∂r ∗ 2 The boundary conditions are s U¯ j = kj∗
U¯ j = U¯ j +1 , ¯ ∗ 1 ∂ Uj λj ∗ ∗ − r ∂r 1 ∂ U¯ n λ∗n ∗ ∗ − r ∂r
(3.1a)
r ∗ = rj∗ , j = 1, . . . , n − 1, U¯ j 1 ∂ U¯ j +1 U¯ j +1 ∗ − ∗ 2 , r ∗ = rj∗ , j = 1, . . . , n − 1, = λj +1 ∗ r ∂r ∗ r∗ 2 r U¯ n U¯ n = α ∗ U¯ ∞ (s) − ∗ , r ∗ = 1. ∗ 2 r r
(3.1b) (3.1c) (3.1d)
The Laplace transform of f (t ∗ ) is defined as [1]: f¯(s) =
∞
exp −st ∗ f t ∗ dt ∗ .
(3.2a)
0
The Laplace transform of a convolution is given by [1]: f1 ∗ f2 (s) = f¯1 (s)f¯2 (s) ⇔ f1 t ∗ ∗ f2 t ∗ =
t ∗
f1 (τ )f2 t ∗ − τ dτ.
(3.2b)
0
The solution of Eq. (3.1a) is easily constructed as U¯ j = Aj sinh qj r ∗ − rj∗−1 + Bj cosh qj r ∗ − rj∗−1 , r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n, (3.3) where qj = s/kj∗ and coefficients Aj and Bj are determined by applying the boundary conditions (3.1b)–(3.1d) as follows: B1 = 0,
(3.4a)
Aj sinh ξj + Bj cosh ξj − Bj +1 = 0,
j = 1, . . . , n − 1,
(3.4b)
Aj (cosh ξj − ηj sinh ξj ) + Bj (sinh ξj − ηj cosh ξj ) − hj Aj +1 +
λj +1 ηj Bj +1 = 0, λj
j = 1, . . . , n − 1,
(3.4c)
An hA + Bn hB = α ∗ U¯ ∞ ,
(3.4d)
where
∗ ∗
kj λ 1 j +1 ∗ ∗ ηj = for j = 1, . . . , n, hj = ∗ ∗ for j = 1, . . . , n − 1, ξj = qj rj − rj −1 , ∗ qj rj λj kj +1 ∗ ∗ α α hA = qn cosh ξn + − ηn sinh ξn , − ηn cosh ξn . hB = qn sinh ξn + qn qn
The coefficients Aj and Bj can be solved from Eq. (3.4) by applying Cramer’s rule. Let
0 sinh ξ1 cosh ξ1 − η1 sinh ξ1 ··· ∆(s) = 0 0 0
λ2 η1 λ1
··· 0
··· ··· ··· ··· ···
0 0 0 ··· sinh ξn−1
0 0 0 ··· cosh ξn−1
0 0 0 ··· 0
0
0
···
cosh ξn−1 − ηn−1 sinh ξn−1
sinh ξn−1 − ηn−1 cosh ξn−1
−hn−1
0
0
···
λn ηn−1 λn−1
0
0
hA
hB
1 cosh ξ1 sinh ξ1 − η1 cosh ξ1 ··· 0
0 0 −h1 ··· 0
0 0
0 −1
0 0 0 ··· −1
(3.5a)
278
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
↓ 0 ··· 0 α ∗ U¯ ∞
Column 2j −1
↓ 0 ··· 0 α ∗ U¯ ∞
Column 2j
··· ··· ··· ··· ··· ··· · · · ··· , , ∆ ∆2j −1 (s) = (s) = 2j ∆(s) ∆(s) ∆(s) ∆(s) ··· ··· ··· ··· det(∆2j −1 (s)) det(∆2j (s)) Aj = Bj = = δ2j −1 (s)U¯ ∞ , = δ2j (s)U¯ ∞ , det(∆(s)) det(∆(s)) where δ2j −1 (s) and δ2j (s) are the matrix determinants obtained from Eqs. (3.5b), (3.5c). Finally, from Eq. (3.3) we get U¯ j = Aj sinh qj r ∗ − rj∗−1 + Bj cosh qj r ∗ − rj∗−1 = Gj s, r ∗ U¯ ∞ , where
Gj s, r ∗ = δ2j −1 (s) sinh qj r ∗ − rj∗−1 + δ2j (s) cosh qj r ∗ − rj∗−1 ,
r ∗ ∈ rj∗−1 , rj∗ , j = 1, . . . , n.
(3.5b)
(3.5c)
(3.6a)
(3.6b)
In deducing the inverse Laplace transform of U¯ j , we shall make an approximation to avoid evaluating residues. Let gj (t ∗ , r ∗ ) be the inverse Laplace transform of Gj (s, r ∗ ). Eq. (3.6a) is then expressed as (see Eq. (3.2b))
∗
Uj t , r
∗
∗
= gj t , r
∗
∗ U∞ t ∗ =
t ∗
∗
gj τ, r U∞
∞ ∞ t − τ dτ = − gj τ, r ∗ U∞ t ∗ − τ dτ
∗
0
∞ ≈
0
gj τ, r ∗ U∞ t ∗ − τ dτ =
0
∞
t∗
gj τ, r ∗ rn exp iω∗ t ∗ − τ + iϕ dτ
0
= rn exp iω∗ t ∗ + iϕ
∞
exp −iω∗ τ gj τ, r ∗ dτ = rn Gj iω∗ , r ∗ exp iω∗ t ∗ + iϕ = Gj iω∗ , r ∗ U∞ t ∗ .
0
(3.7)
↑
Laplace transform of
gj
at
iω∗ = Gj (iω∗ ,r ∗ ),
see Eq. (3.2a)
Note that the inverse Laplace transform gj (t ∗ , r ∗ ) was acting only as a symbolic function. By taking an advantage of the exponential expression of U∞ (t ∗ ), gj (t ∗ , r ∗ ) was replaced by its Laplace transform Gj (iω, r) which is already available in Eq. (3.6b). 4. Solution method for general boundary conditions For a general time-dependent boundary condition, we approximate it using Fourier series as T∞ (t) = a0 +
∞
ak cos(ωk t + ϕk ).
(4.1)
k=1
Hence, the solution can be split up into two subproblems, namely, Uj = Uj1 + Uj2 : (1) solution Uj1 : boundary condition is constant, a0 ; (2) solution Uj2 : boundary condition is time-dependent ∞ k=1 ak cos(ωk t + ϕk ). The second subproblem is equivalent to that with the boundary condition cos(ωt + ϕ) due to the linearity of the system. For the first subproblem with the constant boundary condition, U∞ = a0 , Eq. (3.6a) gives a 0 a0 U¯ j1 = Gj s ∗ , r ∗ U¯ ∞ = Gj s ∗ , r ∗ (4.2) = δ2j −1 (s) sinh qj r ∗ − rj∗−1 + δ2j (s) cosh qj r ∗ − rj∗−1 . s s The matrix determinants δ2j −1 (s) and δ2j (s) are functions of hyperbolic functions sinh ξ and cosh ξ which can be approximated by power series such as sinh ξ = ξ/1! + ξ 3 /3! + · · · = U¯ j1 ≈ Finally, Uj1 =
s/kj∗ (rj∗ − rj∗−1 ) + · · · . Hence, linearisation of Eq. (4.2) gives
const . const 1 ∗ s + const 2
(4.3)
const const 2 ∗ exp − t . const 1 const 1
(4.4)
Note that the series may ascend as powers of
√
s. An approximate inverse of U¯ j1 is then obtained in the same way as that in Eq. (4.3).
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
279
5. Calculation examples Since an approximate formula has been used in reaching the mathematical solution, this section is aimed at assessing the accuracy of the formula. A selected five-layer sphere is presented here as a calculation example. Its shells are numbered by layer one to five from the inner core. A schematic figure is shown in Fig. 2. Thermal properties and dimensions of the composite are listed in Table 1. The composition is an extension of the three-layer slab employed as an exterior wall structure in our test house to evaluate the accuracy of our model for predicting the evolution of the temperature distributions inside the wall subject to interior and exterior climatic conditions. The interior and exterior climatic temperatures were measured. In this calculation, the boundary conditions were taken from such measurements and then fitted with cosine functions as 5 2πt ai cos − ϕi , T∞ (t) = a0 + (5.1) ωi 1
where the fitting parameters are listed in Table 2. The surface heat transfer coefficient, α, is assumed to be 25 W m−2 K−1 . In order to evaluate the accuracy of the analytical solution, calculations were made with both the analytical solution presented in Eq. (3.7) and a numerical model. In the numerical model, the finite volume method was used to do the spatial discretisation procedure for Eq. (2.1). The Crank–Nicolson time marching procedure was then applied to the discretised equations. The final algebraic equations were solved using Newton’s iterative method [10]. In the calculation, the central regions of layers 1 and 3 were chosen as the calculation points which are marked in Fig. 2. Figs. 3, 4 show the comparisons of the transient temperature variations at the calculation points using both analytical and numerical methods. The temperatures were stored in files as hourly values and shown in figures as hourly and daily values. At both calculation points, the maximal discrepancy between the analytical and numerical values is about 0.37 ◦ C at small t with a relative error of 6%. The comparison results for small values of time t are displayed in Fig. 5. To demonstrate the effect of heat capacity (ability to store heat) of the materials, the thicknesses of both layers 2 and 3 were extended to 1.0 m. The boundary conditions were kept the same (Eq. (5.1)). The calculated temperature variations in the central
Fig. 2. Schematic diagram of the five-layer sphere. Table 1 Material properties of the layered sphere, see Fig. 2 Layer 1 2 3 4 5
Thermal conductivity (W m−1 K−1 )
Thermal diffusivity (m2 s−1 )
Shell thickness (mm)
0.12 0.0337 0.9 0.147 0.23
1.5 × 10−7
25 200 100 100 13
1.47 × 10−6 3.75 × 10−7 1.61 × 10−7 4.11 × 10−7
Table 2 Parameters for Eq. (5.1)
a0 17.29892
ω1 120.0
ω2 30.0
ω3 10.0
ω4 5.0
ω5 1.0
ϕ1 28.56693
ϕ2 7.790886
ϕ3 6.483743
ϕ4 3.647276
ϕ5 3.295951
a1 2.3712
a2 −1.77464
a3 −0.53307
a4 −0.05364
a5 −0.12107
280
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
Fig. 3. Comparison of analytical and numerical results in the central region of layer 1.
Fig. 4. Comparison of analytical and numerical results in the central region of layer 3.
Fig. 5. Comparison of analytical and numerical results in the central region of layer 1.
Fig. 6. Comparison of analytical and numerical results in the central region of layer 1. Thicknesses of layers 2 and 3 are extended to 1 m.
region of layer 1 with numerical and analytical models are displayed in Fig. 6. As expected, fluctuating temperature change was obtained in the material with a larger heat capacity (a smaller thermal diffusivity). In this calculation, the maximal discrepancy between the numerical and analytical values is about 0.2 ◦ C (relative error less than 4%). These two examples demonstrate a high accuracy of the developed analytical solution.
6. Discussion
The solution of transient heat conduction for an n-layer sphere is explicitly expressed. We make the following observations. • The calculation includes only simple computation of a matrix determinant which can be easily accomplished by commercial mathematical packages like M APLE, M ATLAB and M ATHEMATICA. For any layer, only five sparse matrices are involved. The calculation load is small and the computing time is short. Furthermore, an increase of a slab’s layer applies only to the dimension of the matrices. The computing load will not be affected much. • The linearisation of G(s, r ∗ ) can be tedious in Eq. (4.2) to obtain the approximately ‘transient’ change for the average temperature if the layer number is more than three. Some mathematical packages, say M APLE, may be helpful which can handle symbolic computations. In authors’ experience, a steady-state solution as an approximate one for the average temperature is quite sufficient if studies do not focus very much on the initial temperature change. • Compared with numerical methods, the developed method is easier to implement and a possible instability in the numerical solution is avoided. The accuracy is good.
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
281
• The main contribution of this Letter, which distinguishes it from all previous work, is an elegant way to perform the inverse Laplace transformation of U¯ j to obtain the solution Uj in Eq. (3.7). The usual way of doing the inverse is to find the inverse Laplace transformation gj (t ∗ , r ∗ ) of the function Gj (s, r ∗ ) by the residue theorem. Once gj (t ∗ , r ∗ ) is available, the exact analytical solution is obtained. It is easy to see that, if the composition has more than two layers, then it is difficult to find the inverse Laplace transformation of the complicated Gj (s, r ∗ ) by the residue theorem. Furthermore, there are many other questions have to be considered when trying to numerically search for roots of the inverse. That is the main reason why a closed-form solution was given for the transient heat conduction problem in a three-layer composition only to date [4], though there has been extensive research on such topic. In this Letter in the inverse calculation of U¯ j , gj (t ∗ , r ∗ ), was acting only as a symbolic function. The real involvement in the calculation is the Laplace transformation Gj (s, r ∗ ) of gj (t ∗ , r ∗ ) which is explicitly expressed in Eq. (3.6). This achievement results in a closed-formed approximate solution but the exact solution. • There are difficulties in making an evaluation of the accuracy of the proposed method mainly dueto the complicated function ∞ of Gj (s, r ∗ ) which makes it very hard to get the inverse gj (t ∗ , r ∗ ). Therefore, the neglected term t ∗ gj (τ, r ∗ )U∞ (t ∗ − τ ) dτ cannot be estimated with various functions of U∞ . However, one thing is for sure that the neglected term approaches to zero if t is big enough. Therefore, a conservative interpretation of the accuracy of the proposed method is that the closed-form analytical solution is exact for quasi-steady-state heat conduction problems. For solving such problems in general, eigenvalue problems often exist for multilayer compositions (e.g. [2]). The proposed method is free of eigenvalue problems. • As we mentioned, the proposed method relies on an approximation in Eqs. (3.7) and (4.3). For small times t , the accuracy may be low which is probably due to a poor approximate in Eq. (4.3). We have made an error estimate of Eq. (3.7) for a simple heat conduction equation in which the exact solution is available. The result suggested that the magnitude of the omitted term is small even for small values of t . Therefore, to increase the accuracy, the analytical method for the transient heat conduction problem in a composite slab with a constant boundary temperature change in [4] can be employed. Then a numerical program is needed for the eigenvalues. • It is known that any periodic function can be represented as its Fourier series. Non-periodic function can be approximated as Fourier series in the extended interval. For both periodic and non-periodic functions, transient heat conduction must be solved during a short time interval compared with its period. In fact, the steady-state may not be reached during several periods or during normal calculational cycles. Therefore, the results obtained in this Letter are transient solutions for general boundary conditions. As demonstrated in the calculation examples, the boundary temperatures are not necessarily periodic. In Fig. 4, the boundary temperature is non-periodic function, even though the approximation involves a sum of cosine functions with a maximal period of 30 days. A transient calculation was made for 20 days which is much shorter than the duration of one period. • Finally, it is easy to see that the assumption of perfect thermal contact between layers can be removed without adding any content into the calculation method. It is only for demonstrational convenience that we prescribed this limitation. 7. Conclusions In this Letter, an analytical approach to heat conduction in a layered sphere subject to a time-dependent temperature change has been presented. The boundary temperature was approximated using Fourier series. The basic technique in obtaining the analytical solution is the Laplace transform which often yields a residue computation in the exact solution. In this approximate method, a novel mathematical computation was made in deducing an approximate inverse Laplace transform. As a result, it is free of the residue evaluation and the developed solution is an approximate. The method is shown to have considerable potential for solving the equations of heat conduction. Compared to analytical solutions that are mainly possible with some assumptions, the developed method has no restrictions and is flexible with regard to the boundary condition. Moreover, the benefit of the results is that its simple and concise mathematical form of the solutions can be used to analyse physical properties in combination with material properties in the heat transfer process. Agreement with the numerical solutions is good. In a general conduction or diffusion application context, however, numerical schemes are usually necessary and may be very complex due to the special geometry of the material. Finally, it is worth mentioning that the analysis for layered spheres cannot be extended straightforwardly to layered cylinders. Acknowledgements The authors would like to thank the Academy of Finland for supporting this research, and the referees for their valuable criticisms of the early manuscript. References [1] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Oxford Univ. Press, Oxford, UK, 1959. [2] P. Snabre, A. Perez, A. Arhaliass, Eur. Phys. J: Appl. Phys. 1 (1998) 315.
282
[3] [4] [5] [6] [7] [8] [9] [10]
X. Lu, M. Viljanen / Physics Letters A 351 (2006) 274–282
M.C.B. Gough, Modelling heat flow in buildings: an eigenfunction approach, Ph.D. Thesis, University of Cambridge Press, 1982. Y. Sun, I.S. Wichman, Int. J. Heat Mass Transfer 47 (2004) 1555. T.P. Fredman, Heat Mass Transfer 39 (2003) 285. C. Tsai, C. Hung, Int. J. Heat Mass Transfer 46 (2003) 5137. A. Haji-Sheikh, J.V. Beck, Int. J. Heat Mass Transfer 45 (2002) 1865. A. Haji-Sheikh, J.V. Beck, D. Agonafer, Int. J. Heat Mass Transfer 46 (2003) 2363. X. Lu, P. Tervola, J. Phys. A: Math. Gen. 38 (2005) 81. S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980.