# Exact Analytical Solution for 3D Time-Dependent Heat Conduction in a Multilayer Sphere with Heat Sources Using Eigenfunction Expansion Method

Authors: Nemat Dalir
Journal: International Scholarly Research Notices
 (2014)

---

## Abstract

An exact analytical solution is obtained for the problem of three-dimensional transient heat conduction in the multilayered sphere. The sphere has multiple layers in the radial direction and, in each layer, time-dependent and spatially nonuniform volumetric internal heat sources are considered. To obtain the temperature distribution, the eigenfunction expansion method is used. An arbitrary combination of homogenous boundary condition of the first or second kind can be applied in the angular and azimuthal directions. Nevertheless, solution is valid for nonhomogeneous boundary conditions of the third kind (convection) in the radial direction. A case study problem for the three-layer quarter-spherical region is solved and the results are discussed.


---

## Body

## 1. Introduction

Multilayer materials are composite media composed of several layers. Because of the additional benefit of combining various mechanical, physical, and thermal properties of different substances, a construction using multilayer elements is of interest. Multilayer materials are used in semicircular fiber insulated heaters, multilayer insulation materials, and nuclear fuel rods. Multilayer transient heat conduction finds applications in thermodynamics, fuel cells, and electrochemical reactors. The layered sphere is utilized to investigate the thermal properties of composite media by assuming embedded spherical particles in the composite matrix. For solving the problems of multilayer transient heat conduction, the same methods which are used in solving problems of single layer transient heat conduction are applied. These methods can be classified into two groups: analytical methods and numerical methods. Analytical methods are advantageous over numerical methods in two ways:(
1
) analytical solutions can be used as benchmark to examine and actually confirm numerical algorithms; (
2
) compared to a discrete numerical solution, the mathematical form of an analytical solution can provide better insight. It should also be mentioned that the analytical methods applied to multilayer transient conduction are analogous to those used in the single-layer transient heat conduction. These analytical methods include Green’s function method, the Laplace transform, separation of variables, and eigenfunction expansion method.Many researchers have solved the transient heat conduction problem in a composite medium. For instance, Salt [1] solved the transient heat conduction problem in a two-dimensional composite slab using an orthogonal eigenfunction expansion technique. Mikhailov and Özişik [2], using the orthogonal expansion approach, solved the problem of transient three-dimensional heat conduction in a composite Cartesian medium. Haji-Sheikh and Beck [3] used Green’s function method to obtain temperature distribution in a three-dimensional two-layer orthotropic slab. de Monte [4, 5] applied the eigenfunction expansion method to obtain the transient temperature distribution for the heat conduction in a two-dimensional two-layer isotropic slab with homogenous boundary conditions. Lu et al. [6] and Lu and Viljanen [7] combined separation of variables and Laplace transforms to solve the transient conduction in the two-dimensional cylindrical and spherical media. Singh et al. [8, 9] and Jain et al. [10, 11] used the combination of separation of variables and eigenfunction expansion methods to solve the two-dimensional multilayer transient heat conduction in spherical coordinates. Recently, Dalir and Nourazar [12] used the eigenfunction expansion method to solve the problem of three-dimensional transient heat conduction in a multilayer cylinder.Singh et al. [8, 9] and Jain et al. [10, 11] have studied 2D multilayer transient conduction problems in spherical and cylindrical coordinates. They have obtained analytical solutions for 2D multilayer transient heat conduction in spherical coordinates, in polar coordinates with multiple layers in the radial direction, and in a multilayer annulus. They have used the method of partial solutions to obtain the temperature distributions. In the method of partial solutions, the nonhomogeneous transient problem is split into two subproblems: a nonhomogeneous steady-state subproblem and a homogeneous transient subproblem. Then, the eigenfunction expansion method is used to solve the nonhomogeneous steady-state subproblem and the method of separation of variables is used to solve the homogeneous transient subproblem.The literature survey for the exact analytical solution for 3D transient heat conduction in multilayered sphere demonstrates that such a solution has not, so far, been developed. Thus, in the present paper, using the eigenfunction expansion method, an analytical triple-series solution for transient heat conduction in the 3D spherical coordinates for radial multilayer domain with spatially nonuniform and time-dependent internal heat sources is obtained. Homogenous boundary conditions of the first or second kind can be applied on surfaces ofθ = constant and ϕ = constant. However, nonhomogeneous boundary conditions of the third kind (convection) [11] are used in the r-direction.Some assumptions are made for the 3D multilayer spherical transient conduction problem. First, the problem is a boundary-value problem of conduction in spherical (r-θ-ϕ coordinates) or part-spherical multilayer geometries. Second, volumetric internal heat sources of nonuniform and time-dependent (r, θ, ϕ and t-dependent) types are present. Third, on the inner and outer radial boundaries, nonhomogeneous boundary conditions of any kind can be used but, on the boundary surfaces in the θ and ϕ-directions, only the first or second kind of homogeneous boundary condition can be applied.


## 2. Mathematical Formulation

An-layer composite spherical slab (
r
0
≤
r
≤
r
n, 0
≤
θ
≤
ψ,  and  0
≤
ϕ
≤
ω
) is considered. All the layers have perfect thermal contact and are presumed to be isotropic in thermal properties. α
i and k
i are the temperature independent thermal diffusivity and thermal conductivity of the ith layer. At t
=
0, the ith layer is at a specified temperature f
i
(
r
,
θ
,
ϕ
) and time dependent heat sources g
i
(
r
,
θ
,
ϕ
,
t
) are switched on in each radial layer. For t
>
0, homogenous boundary conditions of first or second kind are applied to the angular surfaces of θ
=
0 and θ
=
ψ and azimuthal surfaces of ϕ
=
0 and ϕ
=
ω. For the inner (
i
=
1
,
r
=
r
0
) and the outer (
i
=
n
,
r
=
r
n
) radial surfaces, all three kinds of boundary conditions are applicable.The governing differential equation of the 3D transient conduction in a multilayered sphere with heat sources is as follows:(1)
∂
2
T
i
∂
r
2
+
2
r
∂
T
i
∂
r
+
1
r
2
∂
2
T
i
∂
θ
2
+
cot
⁡
θ
r
2
∂
T
i
∂
θ
+
1
r
2
sin
2
⁡
θ
∂
2
T
i
∂
ϕ
2
+
g
i
r
,
θ
,
ϕ
,
t
k
i
=
1
α
i
∂
T
i
∂
t
,
T
i
=
T
i
r
,
θ
,
ϕ
,
t
,
r
0
≤
r
≤
r
n
,
r
i
-
1
≤
r
≤
r
i
,
1
≤
i
≤
n
,
0
≤
θ
≤
ψ
,
ψ
<
π
,
0
≤
ϕ
≤
ω
,
ω
<
2
π
,
t
≥
0
.The boundary conditions are as follows.(i)
Inner surface of 1st layer(
i
=
1
):
(2)
A
in
∂
T
1
∂
r
r
0
,
θ
,
ϕ
,
t
+
B
in
T
1
r
0
,
θ
,
ϕ
,
t
=
C
in
.
(ii)
Outer surface ofnth layer (
i
=
n
):
(3)
A
out
∂
T
n
∂
r
r
n
,
θ
,
ϕ
,
t
+
B
out
T
n
r
n
,
θ
,
ϕ
,
t
=
C
out
.
(iii)
θ
=
ψ surface (
i
=
1,2
,
…
,
n
):
(4)
T
i
r
,
θ
=
ψ
,
ϕ
,
t
=
0
or
∂
T
i
∂
θ
r
,
θ
=
ψ
,
ϕ
,
t
=
0
,
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
i
=
1
,
…
,
n
.
(iv)
ϕ
=
0 surface (
i
=
1,2
,
…
,
n
):
(5)
T
i
r
,
θ
,
ϕ
=
0
,
t
=
0
or
∂
T
i
∂
θ
(
r
,
θ
,
ϕ
=
0
,
t
)
=
0
,
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
i
=
1
,
…
,
n
.
(v)
ϕ
=
ω surface (
i
=
1,2
,
…
,
n
):
(6)
T
i
(
r
,
θ
,
ϕ
=
ω
,
t
)
=
0
or
∂
T
i
∂
θ
(
r
,
θ
,
ϕ
=
ω
,
t
)
=
0
,
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
L
i
=
1
,
…
,
n
.
(vi)
Inner interface of theith layer (
i
=
2
,
…
,
n
):
(7)
T
i
r
i
-
1
,
θ
,
ϕ
,
t
=
T
i
-
1
r
i
-
1
,
θ
,
ϕ
,
t
i
=
2
,
…
,
n
,
k
i
∂
T
i
∂
r
r
i
-
1
,
θ
,
ϕ
,
t
=
k
i
-
1
∂
T
i
-
1
∂
r
r
i
-
1
,
θ
,
ϕ
,
t
i
=
2
,
…
,
n
.
(vii)
Outer interface of theith layer (
i
=
1
,
…
,
n
-
1
):
(8)
T
i
r
i
,
θ
,
ϕ
,
t
=
T
i
+
1
r
i
,
θ
,
ϕ
,
t
i
=
1
,
…
,
n
-
1
,
(9)
k
i
∂
T
i
∂
r
(
r
i
,
θ
,
ϕ
,
t
)
=
k
i
+
1
∂
T
i
+
1
∂
r
(
r
i
,
θ
,
ϕ
,
t
)
i
=
1
,
…
,
n
-
1
.The initial condition is as follows:(10)
T
i
r
,
θ
,
ϕ
,
t
=
0
=
f
i
r
,
θ
,
ϕ
,
1
≤
i
≤
n
.
It is worth mentioning that, at r
=
r
0 and r
=
r
n, boundary conditions of first, second, or third kind are applied by appropriate selection of the coefficients in (2) and (3). It should also be mentioned that zero inner radius (
r
0
=
0
) for multilayered sphere is modeled by assigning zero values to B
in and C
in in (2) [8].


## 3. Solution Methodology

The eigenfunction expansion method is used to solve the problem. In the eigenfunction expansion method, first, by using the associated eigenvalue problem(
∇
2
φ
=
-
λ
2
φ
), the eigenfunctions are attained at every spatial direction of the problem. The associated eigenvalue problem is solved by the use of separation of variables. Afterward, the dependent variable and the available nonhomogeneity in the governing differential equation of the problem are separately written as series expansions of the eigenfunctions. In heat conduction problems, the dependent variable is temperature and the available nonhomogeneity is the volumetric heat source. The series expansions are then substituted into the differential equation. By performing some mathematical manipulations, an ordinary differential equation (ODE) is finally obtained for the independent variable. The solution of the problem is completed by solving this ODE, which is a first order ODE in the case of heat conduction problems.As stated before, the method of partial solutions was used by Jain and Singh [11] for solving 2D transient heat conduction problems, the reason being that the heat source is independent of time. However, the method of partial solutions cannot be used for solving the present 3D transient heat conduction problem because the heat source depends on time. Thus, due to time dependence of the heat source, the partial solution of the steady-state subproblem cannot include the heat source term and the partial solutions method cannot be used. Therefore, to the best knowledge of the authors, the most efficient tool for solving the 3D heat conduction problem of the present paper is the eigenfunction expansion method.For the transient problem of present paper, the associated eigenvalue problem is written as follows:(11)
∇
2
φ
i
=
-
λ
2
φ
i
⟹
∂
2
φ
i
∂
r
2
+
2
r
∂
φ
i
∂
r
+
1
r
2
∂
2
φ
i
∂
θ
2
llllllllllllllilll
+
cot
⁡
θ
r
2
∂
φ
i
∂
θ
+
1
r
2
sin
2
⁡
θ
∂
2
φ
i
∂
ϕ
2
llllllllllllll
=
-
λ
2
φ
i
.
Using the method of separation of variables, (11) is solved as follows:
(12)
φ
i
r
,
θ
,
ϕ
=
R
i
r
Θ
i
θ
Φ
i
ϕ
,
(13)
R
i
′′
Θ
i
Φ
i
R
i
Θ
i
Φ
i
+
2
r
R
i
′
Θ
i
Φ
i
R
i
Θ
i
Φ
i
+
1
r
2
R
i
Θ
i
′′
Φ
i
R
i
Θ
i
Φ
i
+
cot
⁡
θ
r
2
R
i
Θ
i
′
Φ
i
R
i
Θ
i
Φ
i
+
1
r
2
sin
2
⁡
θ
R
i
Θ
i
Φ
i
′′
R
i
Θ
i
Φ
i
=
-
λ
2
R
i
Θ
i
Φ
i
R
i
Θ
i
Φ
i
,
(14)
R
i
′′
R
i
+
2
r
R
i
′
R
i
+
1
r
2
Θ
i
′′
Θ
i
+
cot
⁡
θ
r
2
Θ
i
′
Θ
i
+
1
r
2
sin
2
⁡
θ
Φ
i
′′
Φ
i
=
-
λ
2
,
(15)
R
i
′′
R
i
+
2
r
R
i
′
R
i
+
1
r
2
Θ
i
′′
Θ
i
+
cot
⁡
θ
r
2
Θ
i
′
Θ
i
=
-
1
r
2
sin
2
⁡
θ
Φ
i
′′
Φ
i
-
λ
2
,
(16)
sin
2
⁡
θ
r
2
R
i
′′
R
i
+
2
r
R
i
′
R
i
+
Θ
i
′′
Θ
i
+
cot
⁡
θ
Θ
i
′
Θ
i
+
λ
2
r
2
=
-
Φ
i
′′
Φ
i
-
=
+
ν
2
,
(17)
Φ
i
′′
+
ν
i
m
2
Φ
i
=
0
⟶
Φ
i
m
ϕ
=
c
1
sin
ν
i
m
ϕ
+
c
2
cos
⁡
ν
i
m
ϕ
,
r
2
R
i
′′
R
i
+
2
r
R
i
′
R
i
+
Θ
i
′′
Θ
i
+
cot
⁡
θ
Θ
i
′
Θ
i
+
λ
2
r
2
=
ν
2
sin
2
⁡
θ
,
(18)
⟶
r
2
R
i
′′
R
i
+
2
r
R
i
′
R
i
+
λ
2
r
2
=
-
Θ
i
′′
Θ
i
-
cot
⁡
θ
Θ
i
′
Θ
i
+
ν
2
sin
2
⁡
θ
=
+
β
2
,
(19)
r
2
R
i
′′
+
2
r
R
i
′
+
(
λ
i
p
n
2
r
2
-
β
n
2
)
R
i
=
0
⟶
R
i
p
n
(
r
)
=
1
r
c
3
J
β
n
+
0.5
λ
i
p
n
r
llllllllllllllllllllllllllllllllllllll
+
c
4
Y
β
n
+
0.5
λ
i
p
n
r
,
Θ
i
′′
+
(
cot
⁡
θ
)
Θ
i
′
+
β
2
-
ν
2
sin
2
⁡
θ
Θ
i
=
0
.
By using the following change of variable:
(20)
μ
=
cos
⁡
⁡
θ
⟶
1
-
μ
2
=
sin
2
⁡
θ
the first and second derivatives of Θ with respect to θ are obtained as follows:
(21)
Θ
′
=
d
Θ
i
d
θ
=
d
Θ
i
d
μ
d
μ
d
θ
=
-
sin
θ
d
Θ
i
d
μ
,
Θ
i
′′
=
d
2
Θ
i
d
θ
2
=
d
d
θ
d
Θ
i
d
θ
=
d
d
θ
-
sin
θ
d
Θ
i
d
μ
=
-
cos
⁡
θ
d
Θ
i
d
μ
-
sin
θ
d
d
θ
d
Θ
i
d
μ
=
sin
2
⁡
θ
d
2
Θ
i
d
μ
2
-
cos
⁡
⁡
θ
d
Θ
i
d
μ
.
Substituting (21) in (19) results in the following:
(22)
Θ
i
′′
+
(
cot
⁡
θ
)
Θ
i
′
+
β
2
-
ν
2
sin
2
⁡
θ
Θ
i
=
0
⟶
sin
2
⁡
θ
d
2
Θ
i
d
μ
2
-
cos
⁡
θ
d
Θ
i
d
μ
+
cos
⁡
θ
sin
θ
-
sin
θ
d
Θ
i
d
μ
+
β
2
-
ν
2
sin
2
⁡
θ
Θ
i
=
0
⟶
(
1
-
μ
2
)
d
2
Θ
i
d
μ
2
-
2
μ
d
Θ
i
d
μ
+
β
n
2
-
ν
i
m
2
1
-
μ
2
Θ
i
=
0
.
If β
n
2
=
n
(
n
+
1
), (22) is the associated Legendre equation. Its solution is written as follows:
(23)
Θ
i
n
μ
=
c
5
P
i
n
ν
i
m
(
μ
)
+
c
6
Q
i
n
ν
i
m
(
μ
)
⟶
Θ
i
n
(
θ
)
=
c
5
P
i
n
ν
i
m
(
cos
⁡
θ
)
+
c
6
Q
i
n
ν
i
m
cos
⁡
θ
→
Q
i
n
ν
i
m
cos
⁡
0
=
Q
i
n
ν
i
m
1
=
∞
→
c
6
=
0
Θ
i
n
θ
=
P
i
n
ν
i
m
cos
⁡
θ
.
The problem eigenvalues are ζ
i
m
n
p
2
=
λ
i
p
n
2
+
ν
i
m
2. It should be stated that the heat fluxes continuity at the interfaces of the radial layers gives the following:
(24)
Φ
i
m
=
Φ
m
,
ν
i
m
=
ν
m
,
λ
i
p
n
=
λ
1
p
n
α
1
α
i
.
The eigenfunctions R
i
p
n
(
r
), Θ
i
n
(
θ
), and Φ
i
m
(
ϕ
) in r-, θ-, and ϕ-directions are derived by applying the boundary conditions in each direction in the following equation:
(25)
R
i
p
n
r
=
1
r
c
3
J
n
n
+
1
+
0.5
λ
i
p
n
r
+
c
4
Y
n
n
+
1
+
0.5
λ
i
p
n
r
,
Θ
i
n
θ
=
P
i
n
ν
m
cos
⁡
θ
,
Φ
m
ϕ
=
c
1
sin
ν
m
ϕ
+
c
2
cos
⁡
ν
m
ϕ
.
It is assumed that the solution of the problem is in the form of a triple-series expansion of the derived eigenfunctions as follows:
(26)
T
i
r
,
θ
,
ϕ
,
t
=
∑
m
=
1
∞
∑
n
=
1
∞
∑
p
=
1
∞
T
i
m
n
p
t
R
i
p
n
r
Θ
i
n
θ
Φ
m
ϕ
.
The heat source term is also written as a triple-series expansion of the eigenfunctions such that
(27)
g
i
r
,
θ
,
ϕ
,
t
=
∑
m
=
1
∞
∑
n
=
1
∞
∑
p
=
1
∞
g
i
m
n
p
t
R
i
p
n
r
Θ
i
n
θ
Φ
m
ϕ
,
where the coefficient g
i
m
n
p
(
t
) is obtained by the use of the orthogonality property as follows:
(28)
g
i
m
n
p
t
=
∫
0
ω
∫
0
ψ
∫
r
i
-
1
r
i
g
i
(
r
,
θ
,
ϕ
,
t
)
r
2
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
×
R
i
p
n
r
Θ
i
n
θ
Φ
m
ϕ
d
r
d
θ
d
ϕ
∫
0
ω
×
∫
0
ω
∫
0
ψ
∫
r
i
-
1
r
i
r
2
R
i
p
n
2
r
Θ
i
n
2
θ
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
×
Φ
m
2
ϕ
d
r
d
θ
d
ϕ
∫
0
ω
-
1
.
Substitution of (26) and (27) in (1) results in the following:(29)
T
i
m
n
p
(
t
)
R
i
p
n
′′
(
r
)
Θ
i
n
(
θ
)
Φ
i
m
(
ϕ
)
+
2
r
T
i
m
n
p
(
t
)
R
i
p
n
′
(
r
)
Θ
i
n
(
θ
)
Φ
i
m
(
ϕ
)
+
1
r
2
T
i
m
n
p
(
t
)
R
i
p
n
(
r
)
Θ
i
n
′′
(
θ
)
Φ
i
m
(
ϕ
)
+
cot
⁡
θ
r
2
T
i
m
n
p
(
t
)
R
i
p
n
(
r
)
Θ
i
n
′
(
θ
)
Φ
i
m
(
ϕ
)
+
1
r
2
sin
2
⁡
θ
T
i
m
n
p
t
R
i
p
n
r
Θ
i
n
θ
Φ
i
m
′′
ϕ
+
1
k
i
g
i
m
n
p
t
R
i
p
n
r
Θ
i
n
θ
Φ
i
m
ϕ
=
1
α
i
T
i
m
n
p
′
t
R
i
p
n
r
Θ
i
n
θ
Φ
i
m
ϕ
,
(30)
⟶
d
T
i
m
n
p
t
d
t
+
-
α
i
R
i
p
n
′′
r
R
i
p
n
r
+
2
r
R
i
p
n
′
r
R
i
p
n
r
+
1
r
2
Θ
i
n
′′
θ
Θ
i
n
θ
+
cot
⁡
θ
r
2
Θ
i
n
′
θ
Θ
i
n
θ
+
1
r
2
sin
2
⁡
θ
Φ
i
m
′′
ϕ
Φ
i
m
ϕ
︸
=
Γ
i
m
n
p
T
i
m
n
p
t
=
α
i
k
i
g
i
m
n
p
t
.Equation (30) is a first-order nonhomogeneous ODE and has the following solution:
(31)
T
i
m
n
p
t
=
α
i
k
i
e
-
Γ
i
m
n
p
t
∫
τ
=
0
τ
=
t
g
i
m
n
p
τ
e
Γ
i
m
n
p
τ
d
τ
+
a
i
1
e
-
Γ
i
m
n
p
t
.
Application of the initial condition (10) on (26) gives the following:
(32)
T
i
r
,
θ
,
ϕ
,
t
=
0
=
f
i
(
r
,
θ
,
ϕ
)
=
∑
m
=
1
∞
∑
n
=
1
∞
∑
p
=
1
∞
T
i
m
n
p
0
R
i
p
n
r
Θ
i
n
θ
Φ
i
m
ϕ
.
The coefficient a
i
1 of (31) is found by the use of the orthogonality property for obtained eigenfunctions as follows:
(33)
a
i
1
=
T
i
m
n
p
(
0
)
=
∫
0
ω
∫
0
ψ
∫
r
i
-
1
r
i
f
i
r
,
θ
,
ϕ
r
2
R
i
p
n
r
Θ
i
n
θ
Φ
i
m
ϕ
d
r
d
θ
d
ϕ
∫
0
ω
∫
0
ψ
∫
r
i
-
1
r
i
r
2
R
i
p
n
2
r
Θ
i
n
2
θ
Φ
i
m
2
ϕ
d
r
d
θ
d
ϕ
.
The solution of differential equation (1) having (2) to (9) as boundary conditions and (10) as initial condition is (26) with (31) as the coefficients.


## 4. Case Study Problem

We consider a three-layer quarter-sphere(
0
≤
r
≤
r
3
,
0
≤
θ
≤
π
/
2,  and  0
≤
ϕ
≤
π
) which is initially (
t
=
0
) at uniform unit temperature. For time t
>
0, thermal convection occurs, from the outer radial surface at r
=
r
3 at (zero) ambient temperature. However, the θ
=
0, θ
=
π
/
2, ϕ
=
0, and ϕ
=
π surfaces are at uniform and constant zero temperatures. These boundary conditions lead to the following: A
in
=
1, A
out
=
k
3, B
in
=
0, B
out
=
h, C
in
=
0, and C
out
=
0. Additionally, the uniformly distributed heat source g
i, i
=
1
,
…
,
3, is turned on in each layer at t
=
0. The governing differential equation for the 3D transient heat conduction with heat sources in this three-layer quarter-spherical region is as follows:
(34)
∂
2
T
i
∂
r
2
+
2
r
∂
T
i
∂
r
+
1
r
2
∂
2
T
i
∂
θ
2
+
cot
⁡
θ
r
2
∂
T
i
∂
θ
+
1
r
2
sin
2
⁡
θ
∂
2
T
i
∂
ϕ
2
+
g
i
k
i
=
1
α
i
∂
T
i
∂
t
,
T
i
=
T
i
r
,
θ
,
ϕ
,
t
,
0
≤
r
≤
r
3
,
r
i
-
1
≤
r
≤
r
i
,
1
≤
i
≤
3
,
0
≤
θ
≤
π
2
,
0
≤
ϕ
≤
π
,
t
≥
0
.
The boundary conditions have the following forms:
(35)
∂
T
1
∂
r
0
,
θ
,
ϕ
,
t
=
0
,
(36)
k
3
∂
T
3
∂
r
r
3
,
θ
,
ϕ
,
t
+
h
T
3
r
3
,
θ
,
ϕ
,
t
=
0
,
(37)
T
i
r
,
π
2
,
ϕ
,
t
=
0
,
(38)
T
i
r
,
θ
,
0
,
t
=
0
,
(39)
T
i
(
r
,
θ
,
π
,
t
)
=
0
.
(i)
Inner interface of theith layer (
i
=
2,3
):
(40)
T
i
r
i
-
1
,
θ
,
ϕ
,
t
=
T
i
-
1
r
i
-
1
,
θ
,
ϕ
,
t
,
k
i
∂
T
i
∂
r
r
i
-
1
,
θ
,
ϕ
,
t
=
k
i
-
1
∂
T
i
-
1
∂
r
r
i
-
1
,
θ
,
ϕ
,
t
.
(ii)
Outer interface of theith layer (
i
=
1,2
):
(41)
T
i
r
i
,
θ
,
ϕ
,
t
=
T
i
+
1
r
i
,
θ
,
ϕ
,
t
,
(42)
k
i
∂
T
i
∂
r
r
i
,
θ
,
ϕ
,
t
=
k
i
+
1
∂
T
i
+
1
∂
r
r
i
,
θ
,
ϕ
,
t
.The initial condition is as follows:
(43)
T
i
(
r
,
θ
,
ϕ
,
0
)
=
1
,
1
≤
i
≤
3
.
According to (25), by the use of the eigenfunction expansion method, the following solutions in z-, r-, and θ-directions are obtained from the associated eigenvalue problem:
(44)
Φ
m
ϕ
=
c
1
sin
ν
m
ϕ
+
c
2
cos
⁡
ν
m
ϕ
,
Θ
i
n
θ
=
P
i
n
ν
m
cos
⁡
θ
,
R
i
p
n
r
=
1
r
c
3
J
n
n
+
1
+
0.5
λ
i
p
n
r
+
c
4
Y
n
n
+
1
+
0.5
λ
i
p
n
r
.
The eigenvalues are ζ
i
m
n
p
2
=
λ
i
p
n
2
+
ν
i
m
2. The heat flux continuity conditions at the interfaces imply the following:
(45)
λ
i
p
n
=
λ
1
p
n
α
1
α
i
.
The eigenfunctions R
i
p
n
(
r
), Θ
i
n
(
θ
), and Φ
m
(
ϕ
) in the r, θ- and ϕ-directions are obtained by applying the relevant boundary conditions in each direction. Application of the boundary conditions in the ϕ-direction due to (44) gives the following:
(46)
Φ
m
ϕ
=
c
1
sin
ν
m
ϕ
+
c
2
cos
⁡
ν
m
ϕ
⟶
Φ
m
(
0
)
=
0
⟶
c
2
=
0
⟶
Φ
m
(
θ
)
=
c
1
sin
ν
m
ϕ
Φ
m
(
π
)
=
0
⟶
c
1
sin
ν
m
π
=
0
→
c
1
≠
0
sin
ν
m
π
=
0
=
sin
m
π
.
Then the ϕ-direction eigenvalues and eigenfunction are obtained as follows:
(47)
ν
m
=
m
,
m
=
1
,
2
,
…
Φ
m
ϕ
=
sin
m
ϕ
.
Using the θ-direction boundary condition, (37), on Θ in (44) results in the following relation:
(48)
Θ
i
n
θ
=
P
i
n
m
cos
⁡
θ
→
Θ
i
n
π
/
2
=
0
P
i
n
m
cos
⁡
π
2
=
0
⟶
P
i
n
m
(
0
)
=
0
,
where P
i
n
m
(
0
)
=
0 is only satisfied when n are odd integers; that is, n
=
1
,
3
,
5
,
…. Thus the θ-direction eigenvalues and the eigenfunction are as follows:
(49)
n
=
1
,
3
,
5
,
…
Θ
i
n
θ
=
P
i
n
m
cos
⁡
θ
.
Applying the r-direction boundary conditions, that is, (35) and (36), gives the following:
(50)
R
i
p
n
r
=
1
r
c
3
J
n
n
+
1
+
0.5
λ
i
p
n
r
+
c
4
Y
n
n
+
1
+
0.5
λ
i
p
n
r
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
i
=
1
,
2
,
3
⟶
R
′
(
0
)
=
0
⟶
R
(
0
)
=
finite
→
Y
n
(
n
+
1
)
+
0.5
(
0
)
=
∞
c
4
=
0
⟶
R
i
p
n
(
r
)
=
c
3
1
r
J
n
(
n
+
1
)
+
0.5
(
λ
i
p
n
r
)
k
R
′
(
r
3
)
+
h
R
(
r
3
)
=
0
⟶
k
c
3
1
r
3
J
n
(
n
+
1
)
+
0.5
′
(
λ
i
p
n
r
3
)
-
k
c
3
1
2
r
3
r
3
J
n
n
+
1
+
0.5
(

λ
i
p
n
r
3

)
+
h
c
3
1
r
3
J
n
(
n
+
1
)
+
0.5
(
λ
i
p
n
r
3
)
=
0
→
c
3
≠
0
1
r
3
h
-
k
2
r
3
k
J
n
(
n
+
1
)
+
0.5
′
(
λ
i
p
n
r
3
)
lllllllllllllllll
+
h
-
k
2
r
3
J
n
(
n
+
1
)
+
0.5
(
λ
i
p
n
r
3
)
=
0
⟶
k
J
n
(
n
+
1
)
+
0.5
′
(
λ
i
p
n
r
3
)
+
h
-
k
2
r
3
J
n
n
+
1
+
0.5
(
λ
i
p
n
r
3
)
=
0
.
Thus the r-direction eigencondition and eigenfunction are derived as (
ν
m
=
m
):
(51)
k
J
n
(
n
+
1
)
+
0.5
′
λ
i
p
n
r
3
+
h
-
k
2
r
3
J
n
(
n
+
1
)
+
0.5
λ
i
p
n
r
3
=
0
,
R
i
p
n
r
=
1
r
J
n
n
+
1
+
0.5
λ
i
p
n
r
.
According to (28) to (33), the coefficients g
i
m
n
p
(
t
), a
i
1, and Γ
i
m
n
p are obtained as follows:
(52)
g
i
m
n
p
t
=
∫
0
π
∫
0
π
∫
r
i
-
1
r
i
g
i
r
2
1
r
J
β
n
+
0.5
λ
i
p
n
r
P
i
n
m
θ
llllllllllllllllllllll
×
sin
⁡
m
ϕ
d
r
d
θ
d
ϕ
∫
0
π
×
∫
0
π
∫
0
π
∫
r
i
-
1
r
i
r
J
β
n
+
0.5
2
λ
i
p
n
r
P
i
n
m
θ
2
lllllllllllllllllllllllll
×
sin
2
⁡
m
ϕ
d
r
d
θ
d
ϕ
∫
0
π
-
1
=

∫
0
π
g
i
∫
r
i
-
1
r
i
r
3
/
2
J
β
n
+
0.5
λ
i
p
n
r
d
r
∫
0
π
P
i
n
m
θ
d
θ
l
l
l
l
l
l
×
∫
0
π
sin
⁡
m
ϕ
d
ϕ
×
∫
r
i
-
1
r
i
r
J
β
n
+
0.5
2
λ
i
p
n
r
d
r
∫
0
π
P
i
n
m
θ
2
d
θ
l
l
l
l
l
l
l
l
l
×
∫
0
π
sin
2
⁡
m
ϕ
d
ϕ
∫
r
i
-
1
r
i

-
1
=

-
1
m
cos
⁡
⁡
m
θ
θ
=
0
θ
=
π
g
i
μ
i
n
m
r
J
m
+
1
μ
i
n
m
r
r
=
r
i
-
1
r
=
r
i
l
l
l
l
l
l
×
-
1
m
cos
⁡
⁡
m
θ
θ
=
0
θ
=
π
-
1
m
cos
⁡
⁡
m
ϕ
ϕ
=
0
ϕ
=
π
×


μ
i
n
m
r
2
2
J
m
2
μ
i
n
m
r
+
J
m
+
1
2
μ
i
n
m
r

r
=
r
i
-
1
r
=
r
i
l
l
l
l
l
l
l
l
l
l
×
π
2
×
L
2
μ
i
n
m
r
2
2
-
1
=

1
-
-
1
l
g
i
μ
i
n
m
r
J
m
+
1
μ
i
n
m
r
r
=
r
i
-
1
r
=
r
i
×
1
m
1
-
-
1
m
l
l
l
l
l
l
×
L
l
π
1
-
-
1
l
×


μ
i
n
m
r
2
2
J
m
2
μ
i
n
m
r
+
J
m
+
1
2
μ
i
n
m
r

r
=
r
i
-
1
r
=
r
i
l
l
l
l
l
l
l
l
l
l
×
π
L
4
μ
i
n
m
r
2
2
-
1
⟶
g
i
m
n
p
(
t
)
=
g
i
m
n
p
=
a
i
1
g
i
,
Γ
i
m
n
p
=
-
α
i
×

R
i
p
n
′′
(
r
)
R
i
p
n
(
r
)
+
1
r
R
i
p
n
′
(
r
)
R
i
p
n
(
r
)
+
1
r
2
Θ
i
n
′′
(
θ
)
Θ
i
n
(
θ
)
+
Φ
i
m
′′
(
ϕ
)
Φ
i
m
(
ϕ
)

=
-
α
i
J
m
′′
λ
i
p
n
r
J
m
λ
i
p
n
r
+
1
r
J
m
′
λ
i
p
n
r
J
m
λ
i
p
n
r
-
m
2
r
2
-
m
2
.
Using the form of a triple-series expansion of the obtained eigenfunctions, that is, (26), the solution is as follows:
(53)
T
i
r
,
θ
,
ϕ
,
t
=
∑
m
=
1
∞
∑
n
=
1
∞
∑
p
=
1
∞
T
i
m
n
p
t
1
r
J
β
n
+
0.5
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
×
λ
i
p
n
r
P
i
n
m
θ
sin
⁡
m
ϕ
,
where the coefficient T
i
m
n
p
(
t
) is attained as follows:
(54)
T
i
m
n
p
t
=
α
i
k
i
e
-
Γ
i
m
n
p
t
g
i
m
n
p
1
Γ
i
m
n
p
(
e
Γ
i
m
n
p
t
-
1
)
+
a
i
1
e
-
Γ
i
m
n
p
t
=
α
i
g
i
m
n
p
k
i
Γ
i
m
n
p
+
1
g
i
-
α
i
k
i
Γ
i
m
n
p
g
i
m
n
p
e
-
Γ
i
m
n
p
t
.
Therefore, the solution of (34) having (35) to (42) as boundary conditions and (43) as initial condition is (53) with (54) as the coefficients.


## 5. Conclusions

The exact analytical solution, that is, transient temperature distribution, is derived for the 3D transient heat conduction problem in a multilayered sphere using eigenfunction expansion method. Time-dependent and nonuniform volumetric heat generation is considered in each radial layer. Third kind nonhomogeneous boundary conditions are applied in the radial direction but the first or second kind homogenous boundary conditions are used in the angular and azimuthal directions. The heat conduction in a three-layer quarter-sphere is solved as a case study problem and the temperature distribution is found.



---

*Source: 708723-2014-11-18.xml*