# An Analytical Solution for Transient Heat Conduction in a Composite Slab with Time-Dependent Heat Transfer Coefficient

Authors: Ryoichi Chiba
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

An analytical solution is derived for one-dimensional transient heat conduction in a composite slab consisting ofn layers, whose heat transfer coefficient on an external boundary is an arbitrary function of time. The composite slab, which has thermal contact resistance at n
-
1 interfaces, as well as an arbitrary initial temperature distribution and internal heat generation, convectively exchanges heat at the external boundaries with two different time-varying surroundings. To obtain the analytical solution, the shifting function method is first used, which yields new partial differential equations under conventional types of external boundary conditions. The solution for the derived differential equations is then obtained by means of an orthogonal expansion technique. Numerical calculations are performed for two composite slabs, whose heat transfer coefficient on the heated surface is either an exponential or a trigonometric function of time. The numerical results demonstrate the effects of temporal variations in the heat transfer coefficient on the transient temperature field of composite slabs.


---

## Body

## 1. Introduction

Generally, heat transfer coefficients (HTCs) vary with time, and their variations can be random or, in certain cases, periodic. For example, because of irregular fluid motion around turbine blades, the HTC of the blade surfaces fluctuates [1]. The surface HTCs of the fuel elements in a boiling water reactor and a pebble bed reactor also vary over time [2]. The same holds for the HTCs between a casting and its metal moulds [3], on diesel fuel droplets subjected to transient heating [4], and on solids enveloped by pulsating flows of liquid or gas in internal combustion engines [5]. In addition, during the thermal processing of foods, the pattern of air circulation may be altered around them, resulting in changes of the surface HTC [6]. The time variation in the surface aspects of objects is also one of the possible reasons (e.g., oxidation, dust contamination, and fissuring) [7]. Such a temporal change in the HTC leads to alterations in the object’s temperature.Transient heat conduction with boundary conditions including a time-dependent HTC has been investigated since the late 1960s. Studies using (semi)numerical methods include the combined use of the finite difference and Laplace transform methods [8] and the parameter-group transformation followed by the Runge-Kutta shooting method [9]. Early studies using approximate methods include those of Ivanov et al. [10–12], which transformed the governing equation into a nonlinear equation using the change of variables and omitted the nonlinear term from the derived equation by restricting the research object to thin bodies. Kozlov [13] reduced the original problem to finding a solution to an infinite number of simultaneous ordinary differential equations, approximating it by a finite number. In addition, the application of the Laplace-Carson integral transform [7], finite integral transform methods with time-dependent eigenvalues [14, 15], and Lie point symmetry analysis [16] were reported. Further, an analytical method based on the Laplace transform and bifrequency transfer function was presented, but the inverse transform is extremely difficult [17]. An analytical solution to the heat conduction problem with arbitrary time-dependent HTC at the boundary had not been obtained for a long period because of the difficulty posed by eigenfunctions and eigenvalues that depend on time. However, Lee and colleagues [18] resolved this situation in 2010; they successfully derived explicit-form solutions by means of the shifting function method [19, 20].On the other hand, there have been few studies with composite media. Il’chenko [21] obtained an approximate solution to the heat conduction problem of a two-layered plate with time-dependent HTC by replacing the continuous time function of the HTC with a piecewise function. Prikhod’ko [22] derived a transient temperature solution of a two-layer plate using a special series expansion of Green’s function, while assuming that the temperature is uniform across the thickness in one layer. Yener and Ozisik [23] extended an analytical method based on a finite integral transform [14] to solve the problem for finite composite media with variable boundary condition parameters. All of them are, however, approximate methods, leading to less rigorous analysis.Composite (or multiregion) media have a wide application in many industrial, environmental, and biological fields (e.g., in turbines, heat exchangers, fuel cells, and brake systems). Thus, it is industrially important to study this kind of heat conduction problems for composite media. To be specific, the need to analyse multiregion problems including contact resistances has been underlined owing to the canning of fuel elements in a reactor [24] and design of thermal insulation systems [21, 22]. In those systems, ambient temperature may vary with time as well as HTC. In studying such engineering problems, analytical solutions are highly valuable as they provide greater insight into the solution behaviour, which is typically missing in any numerical solution. They can be used to benchmark numerical solutions as well. The foregoing facts motivate the present study.In this paper, we address the one-dimensional transient heat conduction problem of an infinite composite slab with general time-dependent HTC at one external boundary. The shifting function method developed by Chen et al. [18] and an orthogonal expansion method [25, 26] are used to obtain an analytical solution for the multiregion heat conduction caused by time-varying ambient temperatures. Main differences from the approach for single layered slab [20] are the initial shifting of solution in each layer, the addition of interface conditions, and the eigenfunction expansion of functions based on the Vodicka type of orthogonality relationship, which would be the novelty of the present analytical method. Numerical calculations are performed to verify the solution for several test cases. Subsequently, two numerical examples are given to quantify the effects of the temporal variations in the HTC on the temperature profiles in the slabs.


## 2. Problem Formulation

Let us consider a composite slab consisting ofn layers, as shown in Figure 1, which has finite thickness H but extends to infinity in the other two dimensions. The composite slab has an arbitrary initial temperature distribution T
i
n
(
z
) and is subjected to two different, arbitrary time-varying ambient temperatures, y
1
(
t
) and y
2
(
t
), at its external boundaries. The boundaries have different HTCs, h
1 and h
2
(
t
), one of which depends on time. At the n
-
1 interfacial boundaries, heat is transferred through different constant conductance coefficients r
i, under the condition of heat flux continuity. A time- and space-dependent internal heat generation Q
i
(
z
,
t
) is also considered in each layer. The material of each layer is isotropic and homogeneous with constant thermophysical properties.Figure 1
Schematic of an infinite composite slab with a time-dependent heat transfer coefficient on one external boundary.The temperature distribution in theith layer of n solidly joined layers is given by the diffusion equation:(1)
1
κ
i
∂
T
i
∂
t
=
∂
2
T
i
∂
z
2
+
Q
i
z
,
t
λ
i
for
i
=
1,2
,
…
,
n
.This composite slab is subject to the following initial and external boundary conditions:(2)
T
i
=
T
i
n
z
at
t
=
0
for
i
=
1,2
,
…
,
n
,
λ
1
∂
T
1
∂
z
-
h
1
T
1
-
y
1
t
=
0
at
z
=
0
,
λ
n
∂
T
n
∂
z
+
h
2
t
T
n
-
y
2
t
=
0
at
z
=
H
.At the internal boundaries, the slab is subject to the following interface conditions:(3)
λ
i
∂
T
i
∂
z
+
r
i
T
i
-
T
i
+
1
=
0
at
z
=
z
i
for
i
=
1,2
,
…
,
n
-
1
,
(4)
λ
i
∂
T
i
∂
z
=
λ
i
+
1
∂
T
i
+
1
∂
z
at
z
=
z
i
for
i
=
1,2
,
…
,
n
-
1
,where (3) implies the discontinuity of temperature at the interfaces and (4) states that the heat flux is continuous at the interfaces. When r
i
→
∞, (3) leads to T
i
-
T
i
+
1
=
0, implying the continuity of temperature or perfect thermal contact at the interfaces. For convenience, we nondimensionalise (1) as follows:(5)
1
κ
-
i
∂
θ
i
∂
τ
=
∂
2
θ
i
∂
ξ
2
+
Q
-
i
ξ
,
τ
for
i
=
1,2
,
…
,
n
,and the initial, boundary, and interface conditions may be stated as(6)
θ
i
=
θ
i
n
ξ
at
τ
=
0
for
i
=
1,2
,
…
,
n
,
(7)
∂
θ
1
∂
ξ
-
B
1
θ
1
-
Y
1
τ
=
0
at
ξ
=
0
,
(8)
∂
θ
n
∂
ξ
+
B
2
τ
θ
n
=
ψ
τ
at
ξ
=
1
,
(9)
∂
θ
i
∂
ξ
+
R
i
θ
i
-
θ
i
+
1
=
0
at
ξ
=
ξ
i
for
i
=
1,2
,
…
,
n
-
1
,
(10)
λ
-
i
∂
θ
i
∂
ξ
=
λ
-
i
+
1
∂
θ
i
+
1
∂
ξ
at
ξ
=
ξ
i
for
i
=
1,2
,
…
,
n
-
1
,where the dimensionless quantities are defined as follows:(11)
B
1
=
h
1
H
λ
1
,
B
2
τ
=
h
2
t
H
λ
n
,
Q
-
i
ξ
,
τ
=
Q
i
z
,
t
H
2
λ
i
T
0
,
R
i
=
r
i
H
λ
i
,
Y
j
τ
=
y
j
t
T
0
for
j
=
1,2
,
θ
i
ξ
,
τ
=
T
i
z
,
t
T
0
,
θ
i
n
ξ
=
T
i
n
z
T
0
,
κ
-
i
=
κ
i
κ
0
,
λ
-
i
=
λ
i
λ
0
,
ξ
=
z
H
,
ξ
i
=
z
i
H
,
τ
=
κ
0
t
H
2
,
ψ
τ
=
B
2
τ
Y
2
τ
.For very large values of R
i at the interfacial boundaries, the problem is reduced to that of temperature continuity at the interfaces.To make the following analysis manageable,B
2
(
τ
) is separated in the form of(12)
B
2
τ
=
δ
+
F
τ
,where δ is the initial value of the Biot number, that is, B
2
(
0
). Substituting (12) into (8) yields(13)
∂
θ
n
∂
ξ
+
δ
θ
n
=
-
F
τ
θ
n
+
ψ
τ
at
ξ
=
1
.


## 3. Analysis

### 3.1. Shifting Function Method

To solve the system of (5)–(10), we first introduce the substitution [20](14)
θ
i
ξ
,
τ
=
v
i
ξ
,
τ
+
f
τ
g
i
ξ
for
i
=
1,2
,
…
,
n
,where v
i
(
ξ
,
τ
) is the transformed function, g
i
(
ξ
) is a shifting function to be determined, and f
(
τ
) is the auxiliary function defined as(15)
f
τ
=
-
F
τ
θ
n
1
,
τ
+
ψ
τ
.Substituting (14) into (5)–(7), (9), (10), and (13) yields(16)
1
κ
-
i
∂
v
i
ξ
,
τ
∂
τ
+
d
f
τ
d
τ
g
i
ξ
=
∂
2
v
i
ξ
,
τ
∂
ξ
2
+
f
τ
d
2
g
i
ξ
d
ξ
2
+
Q
-
i
ξ
,
τ
for
i
=
1,2
,
…
,
n
,
(17)
v
i
ξ
,
0
+
f
0
g
i
ξ
=
θ
i
n
ξ
for
i
=
1,2
,
…
,
n
,
(18)
∂
v
1
0
,
τ
∂
ξ
+
f
τ
d
g
1
0
d
ξ
-
B
1
v
1
0
,
τ
+
f
τ
g
1
0
-
Y
1
τ
=
0
,
(19)
∂
v
n
1
,
τ
∂
ξ
+
f
τ
d
g
n
1
d
ξ
+
δ
v
n
1
,
τ
+
f
τ
g
n
1
=
f
τ
,
(20)
∂
v
i
ξ
i
,
τ
∂
ξ
+
f
τ
d
g
i
ξ
i
d
ξ
+
R
i
v
i
ξ
i
,
τ
+
f
τ
g
i
ξ
i
-
v
i
+
1
ξ
i
,
τ
-
f
τ
g
i
+
1
ξ
i
=
0
for
i
=
1,2
,
…
,
n
-
1
,
(21)
λ
-
i
∂
v
i
ξ
i
,
τ
∂
ξ
+
f
τ
d
g
i
ξ
i
d
ξ
=
λ
-
i
+
1
∂
v
i
+
1
ξ
i
,
τ
∂
ξ
+
f
τ
d
g
i
+
1
ξ
i
d
ξ
for
i
=
1,2
,
…
,
n
-
1
.Following the procedure of the shifting function method, we suitably assume the shifting functiong
i
(
ξ
) to be a function satisfying the following system of equations:(22)
d
2
g
i
ξ
d
ξ
2
=
a
for
i
=
1,2
,
…
,
n
,
(23)
d
g
1
0
d
ξ
-
B
1
g
1
0
=
0
,
(24)
d
g
n
1
d
ξ
=
1
,
(25)
g
n
1
=
0
,
(26)
d
g
i
ξ
i
d
ξ
+
R
i
g
i
ξ
i
-
g
i
+
1
ξ
i
=
0
for
i
=
1,2
,
…
,
n
-
1
,
(27)
λ
-
i
d
g
i
ξ
i
d
ξ
=
λ
-
i
+
1
d
g
i
+
1
ξ
i
d
ξ
for
i
=
1,2
,
…
,
n
-
1
,where a is an unknown constant. The general solution of (22) is given by(28)
g
i
ξ
=
a
2
ξ
2
+
b
i
ξ
+
c
i
,where b
i and c
i
(
i
=
1,2
,
…
,
n
) are additional unknown constants. A total of (
2
n
+
1
) unknown constants are readily determined from (23)–(27).Substituting (25) into (14) yields(29)
θ
n
1
,
τ
=
v
n
1
,
τ
.The substitution of (22) and (25) into (16) results in the following partial differential equation:(30)
∂
v
i
ξ
,
τ
∂
τ
-
κ
-
i
∂
2
v
i
ξ
,
τ
∂
ξ
2
-
g
i
ξ
F
τ
v
˙
n
1
,
τ
+
-
g
i
ξ
F
˙
τ
+
κ
-
i
a
F
τ
v
n
1
,
τ
=
-
g
i
ξ
ψ
˙
τ
+
κ
-
i
a
ψ
τ
+
κ
-
i
Q
-
i
ξ
,
τ
for
i
=
1,2
,
…
,
n
,where the “over-dot” denotes the time derivative. Considering (23)–(27), the associated initial, boundary, and interface conditions are obtained as follows:(31)
v
i
ξ
,
0
=
θ
i
n
ξ
-
ψ
0
g
i
ξ
for
i
=
1,2
,
…
,
n
,
(32)
∂
v
1
0
,
τ
∂
ξ
-
B
1
v
1
0
,
τ
-
Y
1
τ
=
0
,
(33)
∂
v
n
1
,
τ
∂
ξ
+
δ
v
n
1
,
τ
=
0
,
(34)
∂
v
i
ξ
i
,
τ
∂
ξ
+
R
i
v
i
ξ
i
,
τ
-
v
i
+
1
ξ
i
,
τ
=
0
for
i
=
1,2
,
…
,
n
-
1
,
(35)
λ
-
i
∂
v
i
ξ
i
,
τ
∂
ξ
=
λ
-
i
+
1
∂
v
i
+
1
ξ
i
,
τ
∂
ξ
for
i
=
1,2
,
…
,
n
-
1
.


### 3.2. Orthogonal Expansion Technique

The next step is to find a solution to the multiregion heat conduction problem expressed by (30)–(35). In the literature, various analytical methods are available for transient heat conduction in composite media, such as the orthogonal expansion technique, Laplace transform method, method of separation of variables, Green’s function method, and finite integral transform technique. In the present work, Vodicka’s method [25], which has been proven effective in analysing multiregion problems with general space- and time-dependent heat sources and nonhomogeneous boundary conditions, is applied to solve the system of (30)–(35).The solution to (30)–(35) is assumed to be(36)
v
i
ξ
,
τ
=
∑
m
=
1
∞
ϕ
m
τ
X
i
m
ξ
-
L
i
ξ
Y
1
τ
for
i
=
1,2
,
…
,
n
,with the following function:(37)
L
i
ξ
=
C
i
ξ
+
D
i
for
i
=
1,2
,
…
,
n
.The constants C
i and D
i are obtained from the following relationships:(38)
d
L
1
0
d
ξ
-
B
1
L
1
0
=
B
1
,
d
L
n
1
d
ξ
+
δ
L
n
1
=
0
,
d
L
i
ξ
i
d
ξ
+
R
i
L
i
ξ
i
-
L
i
+
1
ξ
i
=
0
,
λ
-
i
d
L
i
ξ
i
d
ξ
=
λ
-
i
+
1
d
L
i
+
1
ξ
i
d
ξ
,
for
i
=
1,2
,
…
,
n
-
1
.The functionX
i
m
(
ξ
) is the solution to the eigenvalue problem corresponding to (30)–(35) and is given as follows:(39)
X
i
m
ξ
=
A
i
m
cos
⁡
d
i
m
ξ
+
B
i
m
sin
⁡
d
i
m
ξ
,where d
i
m
=
γ
m
/
κ
-
i and γ
m is an eigenvalue. The unknown constants A
i
m and B
i
m are determined from the following conditions, which are obtained by substituting (36)–(39) into the boundary and interface conditions, (32)–(35):(40)
d
X
1
m
0
d
ξ
-
B
1
X
1
m
0
=
0
,
d
X
n
m
1
d
ξ
+
δ
X
n
m
1
=
0
,
d
X
i
m
ξ
i
d
ξ
+
R
i
X
i
m
ξ
i
-
X
i
+
1
m
ξ
i
=
0
,
λ
-
i
d
X
i
m
ξ
i
d
ξ
=
λ
-
i
+
1
d
X
i
+
1
m
ξ
i
d
ξ
,
for
i
=
1,2
,
…
,
n
-
1
.In other words, a system of 2
n linear homogeneous equations is obtained for the constants A
i
m and B
i
m. The eigenvalues γ
m
(
m
=
1,2
,
…
) are obtained from the condition under which this system of equations has a nontrivial solution and are therefore positive roots of the following transcendental equation in matrix notation:(41)
G
·
E
1
·
E
2
⋯
E
n
-
1
·
u
=
0
,where(42)
G
=
-
B
1
d
1
m
,
E
i
=
J
i
·
K
i
+
1
,
J
i
=
λ
-
i
d
i
m
cos
⁡
d
i
m
ξ
i
-
d
i
m
cos
⁡
d
i
m
ξ
i
-
R
i
sin
⁡
d
i
m
ξ
i
λ
-
i
d
i
m
sin
⁡
d
i
m
ξ
i
-
d
i
m
sin
⁡
d
i
m
ξ
i
+
R
i
cos
⁡
d
i
m
ξ
i
,
K
i
+
1
=
R
i
cos
⁡
d
i
+
1
m
ξ
i
R
i
sin
⁡
d
i
+
1
m
ξ
i
-
λ
-
i
+
1
d
i
+
1
m
sin
⁡
d
i
+
1
m
ξ
i
λ
-
i
+
1
d
i
+
1
m
cos
⁡
d
i
+
1
m
ξ
i
,
u
=
d
n
m
cos
⁡
d
n
m
+
δ
sin
⁡
d
n
m
d
n
m
sin
⁡
d
n
m
-
δ
cos
⁡
d
n
m
.By substituting (36) into (31), the following equation is obtained:(43)
W
i
ξ
=
∑
m
=
1
∞
ϕ
m
0
X
i
m
ξ
=
θ
i
n
ξ
-
ψ
0
g
i
ξ
+
L
i
ξ
Y
1
0
.The eigenfunction X
i
m
(
ξ
) given by (39) has an orthogonal relationship with discontinuous weighting function; it is expressed as follows [26]:(44)
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
X
i
l
ξ
X
i
k
ξ
d
ξ
=
c
o
n
s
t
.
l
=
k
0
l
≠
k
.Functions g
i
(
ξ
), W
i
(
ξ
), L
i
(
ξ
), and κ
-
i
Q
-
i
(
ξ
,
τ
) and a constant κ
-
i can be expanded in an infinite series of X
i
m
(
ξ
):(45)
g
i
ξ
=
∑
m
=
1
∞
α
m
X
i
m
ξ
,
W
i
ξ
=
∑
m
=
1
∞
w
m
X
i
m
ξ
,
L
i
ξ
=
∑
m
=
1
∞
l
m
X
i
m
ξ
,
κ
-
i
Q
-
i
ξ
,
τ
=
∑
m
=
1
∞
q
m
τ
X
i
m
ξ
,
κ
-
i
=
∑
m
=
1
∞
η
m
X
i
m
ξ
,where the expansion coefficients are given by(46)
α
m
=
1
N
m
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
g
i
ξ
X
i
m
ξ
d
ξ
,
w
m
=
1
N
m
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
W
i
ξ
X
i
m
ξ
d
ξ
,
l
m
=
1
N
m
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
L
i
ξ
X
i
m
ξ
d
ξ
,
q
m
τ
=
1
N
m
∑
i
=
1
n
λ
-
i
∫
ξ
i
-
1
ξ
i
Q
-
i
ξ
,
τ
X
i
m
ξ
d
ξ
,
η
m
=
1
N
m
∑
i
=
1
n
λ
-
i
∫
ξ
i
-
1
ξ
i
X
i
m
ξ
d
ξ
,in which the norm N
m is(47)
N
m
=
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
X
i
m
ξ
2
d
ξ
.The substitution of (36) and (45) into (30) yields the following equation:(48)
∑
m
=
1
∞
ϕ
˙
m
τ
X
i
m
ξ
-
l
m
Y
˙
1
τ
X
i
m
ξ
+
γ
m
2
ϕ
m
τ
X
i
m
ξ
-
g
i
ξ
F
τ
ϕ
˙
m
τ
X
n
m
1
+
-
g
i
ξ
F
˙
τ
+
a
F
τ
ϕ
m
τ
X
n
m
1
+
α
m
L
n
1
X
i
m
ξ
F
τ
Y
˙
1
τ
+
F
˙
τ
Y
1
τ
-
a
F
τ
L
n
1
Y
1
τ
η
m
X
i
m
ξ
+
α
m
ψ
˙
τ
X
i
m
ξ
-
a
ψ
τ
η
m
X
i
m
ξ
-
q
m
τ
X
i
m
ξ
=
0
.Equation (48) is satisfied if the expression in the braces is zero; that is,(49)
ϕ
˙
m
τ
X
i
m
ξ
+
γ
m
2
ϕ
m
τ
X
i
m
ξ
-
g
i
ξ
F
τ
ϕ
˙
m
τ
X
n
m
1
+
-
g
i
ξ
F
˙
τ
+
a
F
τ
ϕ
m
τ
X
n
m
1
=
l
m
Y
˙
1
τ
X
i
m
ξ
-
α
m
L
n
1
X
i
m
ξ
F
τ
Y
˙
1
τ
+
F
˙
τ
Y
1
τ
+
a
F
τ
L
n
1
Y
1
τ
η
m
X
i
m
ξ
-
α
m
ψ
˙
τ
X
i
m
ξ
+
a
ψ
τ
η
m
X
i
m
ξ
+
q
m
τ
X
i
m
ξ
.Equation (49) is multiplied by (
λ
-
i
/
κ
-
i
)
X
i
m
(
ξ
), integrated over the region [
ξ
i
-
1
,
ξ
i
], and then summed up for i
=
1,2
,
…
,
n to yield(50)
ϕ
˙
m
τ
+
γ
m
2
-
β
m
F
˙
τ
+
a
χ
m
F
τ
1
-
β
m
F
τ
ϕ
m
τ
=
ς
m
τ
,where β
m and χ
m are given by(51)
β
m
=
X
n
m
1
·
α
m
,
χ
m
=
X
n
m
1
·
η
m
,and ζ
m
(
τ
) is defined as(52)
ς
m
τ
=
Y
˙
1
τ
l
m
-
F
τ
L
n
1
α
m
+
L
n
1
Y
1
τ
a
F
τ
η
m
-
α
m
F
˙
τ
-
ψ
˙
τ
α
m
+
a
ψ
τ
η
m
+
q
m
τ
1
-
β
m
F
τ
.By solving (50) with the condition ϕ
m
(
0
) = w
m, which is obtained from the comparison between (43) and (45)2, we obtain ϕ
m
(
τ
) as(53)
ϕ
m
τ
=
e
-
∫
0
τ
γ
m
2
-
β
m
F
˙
s
+
a
χ
m
F
s
/
1
-
β
m
F
s
d
s
w
m
+
∫
0
τ
ζ
m
μ
·
e
∫
0
μ
γ
m
2
-
β
m
F
˙
s
+
a
χ
m
F
s
/
1
-
β
m
F
s
d
s
d
μ
.After substituting (28), (29), (36), (37), and (39) back to (14), the dimensionless temperature in the ith layer ξ
∈
[
ξ
i
-
1
,
ξ
i
] of the composite slab is derived as(54)
θ
i
ξ
,
τ
=
∑
m
=
1
∞
ϕ
m
τ
A
i
m
cos
⁡
d
i
m
ξ
+
B
i
m
sin
⁡
d
i
m
ξ
-
a
2
ξ
2
+
b
i
ξ
+
c
i
F
τ
A
n
m
cos
⁡
d
n
m
+
B
n
m
sin
⁡
d
n
m
-
Y
1
τ
C
i
ξ
+
D
i
-
a
2
ξ
2
+
b
i
ξ
+
c
i
F
τ
C
n
+
D
n
+
a
2
ξ
2
+
b
i
ξ
+
c
i
ψ
τ
.In the methodology for deriving the temperature solution, main differences from the counterpart for single layer [20] are the following: (i) the solution is first shifted in each layer by using not a single shifting function across the total thickness but shifting functions that differ for each layer (see (14)); (ii) since there exist interfaces, the interface conditions (see (20) and (21)) are added; and (iii) various functions are expanded by using eigenfunctions that fulfil an orthogonal relationship withdiscontinuous weighting function (see (44)).


## 3.1. Shifting Function Method

To solve the system of (5)–(10), we first introduce the substitution [20](14)
θ
i
ξ
,
τ
=
v
i
ξ
,
τ
+
f
τ
g
i
ξ
for
i
=
1,2
,
…
,
n
,where v
i
(
ξ
,
τ
) is the transformed function, g
i
(
ξ
) is a shifting function to be determined, and f
(
τ
) is the auxiliary function defined as(15)
f
τ
=
-
F
τ
θ
n
1
,
τ
+
ψ
τ
.Substituting (14) into (5)–(7), (9), (10), and (13) yields(16)
1
κ
-
i
∂
v
i
ξ
,
τ
∂
τ
+
d
f
τ
d
τ
g
i
ξ
=
∂
2
v
i
ξ
,
τ
∂
ξ
2
+
f
τ
d
2
g
i
ξ
d
ξ
2
+
Q
-
i
ξ
,
τ
for
i
=
1,2
,
…
,
n
,
(17)
v
i
ξ
,
0
+
f
0
g
i
ξ
=
θ
i
n
ξ
for
i
=
1,2
,
…
,
n
,
(18)
∂
v
1
0
,
τ
∂
ξ
+
f
τ
d
g
1
0
d
ξ
-
B
1
v
1
0
,
τ
+
f
τ
g
1
0
-
Y
1
τ
=
0
,
(19)
∂
v
n
1
,
τ
∂
ξ
+
f
τ
d
g
n
1
d
ξ
+
δ
v
n
1
,
τ
+
f
τ
g
n
1
=
f
τ
,
(20)
∂
v
i
ξ
i
,
τ
∂
ξ
+
f
τ
d
g
i
ξ
i
d
ξ
+
R
i
v
i
ξ
i
,
τ
+
f
τ
g
i
ξ
i
-
v
i
+
1
ξ
i
,
τ
-
f
τ
g
i
+
1
ξ
i
=
0
for
i
=
1,2
,
…
,
n
-
1
,
(21)
λ
-
i
∂
v
i
ξ
i
,
τ
∂
ξ
+
f
τ
d
g
i
ξ
i
d
ξ
=
λ
-
i
+
1
∂
v
i
+
1
ξ
i
,
τ
∂
ξ
+
f
τ
d
g
i
+
1
ξ
i
d
ξ
for
i
=
1,2
,
…
,
n
-
1
.Following the procedure of the shifting function method, we suitably assume the shifting functiong
i
(
ξ
) to be a function satisfying the following system of equations:(22)
d
2
g
i
ξ
d
ξ
2
=
a
for
i
=
1,2
,
…
,
n
,
(23)
d
g
1
0
d
ξ
-
B
1
g
1
0
=
0
,
(24)
d
g
n
1
d
ξ
=
1
,
(25)
g
n
1
=
0
,
(26)
d
g
i
ξ
i
d
ξ
+
R
i
g
i
ξ
i
-
g
i
+
1
ξ
i
=
0
for
i
=
1,2
,
…
,
n
-
1
,
(27)
λ
-
i
d
g
i
ξ
i
d
ξ
=
λ
-
i
+
1
d
g
i
+
1
ξ
i
d
ξ
for
i
=
1,2
,
…
,
n
-
1
,where a is an unknown constant. The general solution of (22) is given by(28)
g
i
ξ
=
a
2
ξ
2
+
b
i
ξ
+
c
i
,where b
i and c
i
(
i
=
1,2
,
…
,
n
) are additional unknown constants. A total of (
2
n
+
1
) unknown constants are readily determined from (23)–(27).Substituting (25) into (14) yields(29)
θ
n
1
,
τ
=
v
n
1
,
τ
.The substitution of (22) and (25) into (16) results in the following partial differential equation:(30)
∂
v
i
ξ
,
τ
∂
τ
-
κ
-
i
∂
2
v
i
ξ
,
τ
∂
ξ
2
-
g
i
ξ
F
τ
v
˙
n
1
,
τ
+
-
g
i
ξ
F
˙
τ
+
κ
-
i
a
F
τ
v
n
1
,
τ
=
-
g
i
ξ
ψ
˙
τ
+
κ
-
i
a
ψ
τ
+
κ
-
i
Q
-
i
ξ
,
τ
for
i
=
1,2
,
…
,
n
,where the “over-dot” denotes the time derivative. Considering (23)–(27), the associated initial, boundary, and interface conditions are obtained as follows:(31)
v
i
ξ
,
0
=
θ
i
n
ξ
-
ψ
0
g
i
ξ
for
i
=
1,2
,
…
,
n
,
(32)
∂
v
1
0
,
τ
∂
ξ
-
B
1
v
1
0
,
τ
-
Y
1
τ
=
0
,
(33)
∂
v
n
1
,
τ
∂
ξ
+
δ
v
n
1
,
τ
=
0
,
(34)
∂
v
i
ξ
i
,
τ
∂
ξ
+
R
i
v
i
ξ
i
,
τ
-
v
i
+
1
ξ
i
,
τ
=
0
for
i
=
1,2
,
…
,
n
-
1
,
(35)
λ
-
i
∂
v
i
ξ
i
,
τ
∂
ξ
=
λ
-
i
+
1
∂
v
i
+
1
ξ
i
,
τ
∂
ξ
for
i
=
1,2
,
…
,
n
-
1
.


## 3.2. Orthogonal Expansion Technique

The next step is to find a solution to the multiregion heat conduction problem expressed by (30)–(35). In the literature, various analytical methods are available for transient heat conduction in composite media, such as the orthogonal expansion technique, Laplace transform method, method of separation of variables, Green’s function method, and finite integral transform technique. In the present work, Vodicka’s method [25], which has been proven effective in analysing multiregion problems with general space- and time-dependent heat sources and nonhomogeneous boundary conditions, is applied to solve the system of (30)–(35).The solution to (30)–(35) is assumed to be(36)
v
i
ξ
,
τ
=
∑
m
=
1
∞
ϕ
m
τ
X
i
m
ξ
-
L
i
ξ
Y
1
τ
for
i
=
1,2
,
…
,
n
,with the following function:(37)
L
i
ξ
=
C
i
ξ
+
D
i
for
i
=
1,2
,
…
,
n
.The constants C
i and D
i are obtained from the following relationships:(38)
d
L
1
0
d
ξ
-
B
1
L
1
0
=
B
1
,
d
L
n
1
d
ξ
+
δ
L
n
1
=
0
,
d
L
i
ξ
i
d
ξ
+
R
i
L
i
ξ
i
-
L
i
+
1
ξ
i
=
0
,
λ
-
i
d
L
i
ξ
i
d
ξ
=
λ
-
i
+
1
d
L
i
+
1
ξ
i
d
ξ
,
for
i
=
1,2
,
…
,
n
-
1
.The functionX
i
m
(
ξ
) is the solution to the eigenvalue problem corresponding to (30)–(35) and is given as follows:(39)
X
i
m
ξ
=
A
i
m
cos
⁡
d
i
m
ξ
+
B
i
m
sin
⁡
d
i
m
ξ
,where d
i
m
=
γ
m
/
κ
-
i and γ
m is an eigenvalue. The unknown constants A
i
m and B
i
m are determined from the following conditions, which are obtained by substituting (36)–(39) into the boundary and interface conditions, (32)–(35):(40)
d
X
1
m
0
d
ξ
-
B
1
X
1
m
0
=
0
,
d
X
n
m
1
d
ξ
+
δ
X
n
m
1
=
0
,
d
X
i
m
ξ
i
d
ξ
+
R
i
X
i
m
ξ
i
-
X
i
+
1
m
ξ
i
=
0
,
λ
-
i
d
X
i
m
ξ
i
d
ξ
=
λ
-
i
+
1
d
X
i
+
1
m
ξ
i
d
ξ
,
for
i
=
1,2
,
…
,
n
-
1
.In other words, a system of 2
n linear homogeneous equations is obtained for the constants A
i
m and B
i
m. The eigenvalues γ
m
(
m
=
1,2
,
…
) are obtained from the condition under which this system of equations has a nontrivial solution and are therefore positive roots of the following transcendental equation in matrix notation:(41)
G
·
E
1
·
E
2
⋯
E
n
-
1
·
u
=
0
,where(42)
G
=
-
B
1
d
1
m
,
E
i
=
J
i
·
K
i
+
1
,
J
i
=
λ
-
i
d
i
m
cos
⁡
d
i
m
ξ
i
-
d
i
m
cos
⁡
d
i
m
ξ
i
-
R
i
sin
⁡
d
i
m
ξ
i
λ
-
i
d
i
m
sin
⁡
d
i
m
ξ
i
-
d
i
m
sin
⁡
d
i
m
ξ
i
+
R
i
cos
⁡
d
i
m
ξ
i
,
K
i
+
1
=
R
i
cos
⁡
d
i
+
1
m
ξ
i
R
i
sin
⁡
d
i
+
1
m
ξ
i
-
λ
-
i
+
1
d
i
+
1
m
sin
⁡
d
i
+
1
m
ξ
i
λ
-
i
+
1
d
i
+
1
m
cos
⁡
d
i
+
1
m
ξ
i
,
u
=
d
n
m
cos
⁡
d
n
m
+
δ
sin
⁡
d
n
m
d
n
m
sin
⁡
d
n
m
-
δ
cos
⁡
d
n
m
.By substituting (36) into (31), the following equation is obtained:(43)
W
i
ξ
=
∑
m
=
1
∞
ϕ
m
0
X
i
m
ξ
=
θ
i
n
ξ
-
ψ
0
g
i
ξ
+
L
i
ξ
Y
1
0
.The eigenfunction X
i
m
(
ξ
) given by (39) has an orthogonal relationship with discontinuous weighting function; it is expressed as follows [26]:(44)
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
X
i
l
ξ
X
i
k
ξ
d
ξ
=
c
o
n
s
t
.
l
=
k
0
l
≠
k
.Functions g
i
(
ξ
), W
i
(
ξ
), L
i
(
ξ
), and κ
-
i
Q
-
i
(
ξ
,
τ
) and a constant κ
-
i can be expanded in an infinite series of X
i
m
(
ξ
):(45)
g
i
ξ
=
∑
m
=
1
∞
α
m
X
i
m
ξ
,
W
i
ξ
=
∑
m
=
1
∞
w
m
X
i
m
ξ
,
L
i
ξ
=
∑
m
=
1
∞
l
m
X
i
m
ξ
,
κ
-
i
Q
-
i
ξ
,
τ
=
∑
m
=
1
∞
q
m
τ
X
i
m
ξ
,
κ
-
i
=
∑
m
=
1
∞
η
m
X
i
m
ξ
,where the expansion coefficients are given by(46)
α
m
=
1
N
m
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
g
i
ξ
X
i
m
ξ
d
ξ
,
w
m
=
1
N
m
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
W
i
ξ
X
i
m
ξ
d
ξ
,
l
m
=
1
N
m
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
L
i
ξ
X
i
m
ξ
d
ξ
,
q
m
τ
=
1
N
m
∑
i
=
1
n
λ
-
i
∫
ξ
i
-
1
ξ
i
Q
-
i
ξ
,
τ
X
i
m
ξ
d
ξ
,
η
m
=
1
N
m
∑
i
=
1
n
λ
-
i
∫
ξ
i
-
1
ξ
i
X
i
m
ξ
d
ξ
,in which the norm N
m is(47)
N
m
=
∑
i
=
1
n
λ
-
i
κ
-
i
∫
ξ
i
-
1
ξ
i
X
i
m
ξ
2
d
ξ
.The substitution of (36) and (45) into (30) yields the following equation:(48)
∑
m
=
1
∞
ϕ
˙
m
τ
X
i
m
ξ
-
l
m
Y
˙
1
τ
X
i
m
ξ
+
γ
m
2
ϕ
m
τ
X
i
m
ξ
-
g
i
ξ
F
τ
ϕ
˙
m
τ
X
n
m
1
+
-
g
i
ξ
F
˙
τ
+
a
F
τ
ϕ
m
τ
X
n
m
1
+
α
m
L
n
1
X
i
m
ξ
F
τ
Y
˙
1
τ
+
F
˙
τ
Y
1
τ
-
a
F
τ
L
n
1
Y
1
τ
η
m
X
i
m
ξ
+
α
m
ψ
˙
τ
X
i
m
ξ
-
a
ψ
τ
η
m
X
i
m
ξ
-
q
m
τ
X
i
m
ξ
=
0
.Equation (48) is satisfied if the expression in the braces is zero; that is,(49)
ϕ
˙
m
τ
X
i
m
ξ
+
γ
m
2
ϕ
m
τ
X
i
m
ξ
-
g
i
ξ
F
τ
ϕ
˙
m
τ
X
n
m
1
+
-
g
i
ξ
F
˙
τ
+
a
F
τ
ϕ
m
τ
X
n
m
1
=
l
m
Y
˙
1
τ
X
i
m
ξ
-
α
m
L
n
1
X
i
m
ξ
F
τ
Y
˙
1
τ
+
F
˙
τ
Y
1
τ
+
a
F
τ
L
n
1
Y
1
τ
η
m
X
i
m
ξ
-
α
m
ψ
˙
τ
X
i
m
ξ
+
a
ψ
τ
η
m
X
i
m
ξ
+
q
m
τ
X
i
m
ξ
.Equation (49) is multiplied by (
λ
-
i
/
κ
-
i
)
X
i
m
(
ξ
), integrated over the region [
ξ
i
-
1
,
ξ
i
], and then summed up for i
=
1,2
,
…
,
n to yield(50)
ϕ
˙
m
τ
+
γ
m
2
-
β
m
F
˙
τ
+
a
χ
m
F
τ
1
-
β
m
F
τ
ϕ
m
τ
=
ς
m
τ
,where β
m and χ
m are given by(51)
β
m
=
X
n
m
1
·
α
m
,
χ
m
=
X
n
m
1
·
η
m
,and ζ
m
(
τ
) is defined as(52)
ς
m
τ
=
Y
˙
1
τ
l
m
-
F
τ
L
n
1
α
m
+
L
n
1
Y
1
τ
a
F
τ
η
m
-
α
m
F
˙
τ
-
ψ
˙
τ
α
m
+
a
ψ
τ
η
m
+
q
m
τ
1
-
β
m
F
τ
.By solving (50) with the condition ϕ
m
(
0
) = w
m, which is obtained from the comparison between (43) and (45)2, we obtain ϕ
m
(
τ
) as(53)
ϕ
m
τ
=
e
-
∫
0
τ
γ
m
2
-
β
m
F
˙
s
+
a
χ
m
F
s
/
1
-
β
m
F
s
d
s
w
m
+
∫
0
τ
ζ
m
μ
·
e
∫
0
μ
γ
m
2
-
β
m
F
˙
s
+
a
χ
m
F
s
/
1
-
β
m
F
s
d
s
d
μ
.After substituting (28), (29), (36), (37), and (39) back to (14), the dimensionless temperature in the ith layer ξ
∈
[
ξ
i
-
1
,
ξ
i
] of the composite slab is derived as(54)
θ
i
ξ
,
τ
=
∑
m
=
1
∞
ϕ
m
τ
A
i
m
cos
⁡
d
i
m
ξ
+
B
i
m
sin
⁡
d
i
m
ξ
-
a
2
ξ
2
+
b
i
ξ
+
c
i
F
τ
A
n
m
cos
⁡
d
n
m
+
B
n
m
sin
⁡
d
n
m
-
Y
1
τ
C
i
ξ
+
D
i
-
a
2
ξ
2
+
b
i
ξ
+
c
i
F
τ
C
n
+
D
n
+
a
2
ξ
2
+
b
i
ξ
+
c
i
ψ
τ
.In the methodology for deriving the temperature solution, main differences from the counterpart for single layer [20] are the following: (i) the solution is first shifted in each layer by using not a single shifting function across the total thickness but shifting functions that differ for each layer (see (14)); (ii) since there exist interfaces, the interface conditions (see (20) and (21)) are added; and (iii) various functions are expanded by using eigenfunctions that fulfil an orthogonal relationship withdiscontinuous weighting function (see (44)).


## 4. Verification of Analytical Solution

To illustrate the applicability of the solution derived, first, a single-layered heat conduction problem, whose solution was previously known [20], is modelled as a multilayer problem, with adjacent layers having identical physical properties. This problem (hereinafter called Case A) involves a time-dependent HTC (or Biot number) on one boundary surface of the slab. The following numerical parameters are adopted [20]: θ
i
n = 0.336, B
1 = 0, B
2
(
τ
) = 1.2
-
exp
⁡
(
-
2
τ
), Y
2(τ) = 1.5
-
0.5
exp
⁡
(
-
τ
), and Q
-
i
=
0. For the numerical computation, we use the analytical solution, (see (54)), taking n
=
5 and ξ
i
(
i
=
1,2
,
3,4
) evenly spaced between 0 and 1, assuming perfect contact at the interfaces and uniform thermal conductivity and diffusivity within the entire region; that is, R
i
=
∞, λ
-
i
=
1, and κ
-
i
=
1.Second, we consider another test case with a constant Biot number to verify our analytical solution for multiregion heat conduction. This test case (called conveniently Case B here) represents transient diffusion in a composite slab consisting ofn
=
10 layers with evenly spaced ξ
i. The material properties are taken as λ
-
i
=
κ
-
i
=
1 for i
=
1,3
,
5,7
,
9 and λ
-
i
=
κ
-
i
=
0.1 for i
=
2,4
,
6,8
,
10 [27]. The constants appearing in the initial and boundary conditions are θ
i
n = 0, B
1
=
∞, B
2 = 0, and Y
1 = 1. Moreover, imperfect thermal contact is applied at the interfaces (R
i
=
0.5
/
κ
-
i
,
i
=
1,2
,
…
,
9).Figures2 and 3 plot the analytical solution (see (54)) and numerical results previously published for Cases A and B, respectively. The derived analytical solution is computed using the first 20 terms in the infinite series. In both cases, the analytical solution is in good agreement with the solutions of the existing literature. Tables 1 and 2 show the relative error defined by(55)
e
r
r
o
r
=
max
i
=
1,2
,
…
,
n
⁡
θ
i
ξ
i
,
τ
-
Θ
ξ
i
,
τ
max
i
=
1,2
,
…
,
n
⁡
θ
i
ξ
i
,
τ
,where θ
i is the derived solution and Θ is the solution obtained using the analytical method [20] for Case A and one of either the semianalytical method [27] or Fokas method [28] for Case B. As shown in the tables, the relative error is within approximately 0.2%.Table 1
Relative error in Case A.
Fourier number(
τ
)
1
4
Error
2.08 × 10−3
1.39 × 10−4Table 2
Relative error in Case B.
Fourier number(
τ
)
0.007
2
10
Error (semianalytical)
1.56 × 10−3
5.64 × 10−4
8.62 × 10−4
Error (Fokas method)
5.32 × 10−5
5.04 × 10−5
5.01 × 10−5Figure 2
Solution verification for Case A: the derived analytical solution is plotted as a solid line and the solution presented in [20] is plotted as circles.Figure 3
Solution verification for Case B: the derived analytical solution is plotted as a solid line and the solution presented in [27] is plotted as circles.


## 5. Illustrative Examples

In this section, we apply the analytical solution described in Section3 to two industrial applications. Because we focus on the effect of time-dependent HTC at the heated surface on the transient temperature distribution of composite media, only the cases with constant ambient temperature and no internal heat generation are considered here.


### 5.1. Bilayered Plate with Exponentially Varying HTC

Consider the case of a bilayered plate of aluminium and lead, whose initial temperature is zero. Equal thicknesses for both layers are used, which leads toξ
1 = 1/2, with perfect thermal contact. For τ
>
0, the surrounding medium on the side of lead is raised to a temperature of Y
2 = 1, while the other surface of the plate is insulated (B
1 = 0). For the Biot number at ξ = 1, we choose B
2
τ
=
p
1
-
p
2
exp
⁡
(
-
s
τ
), where p
1, p
2, and s are arbitrary constants. The following material properties have been used in the calculations [26]: λ
-
1
=
5.850, λ
-
2
=
1, κ
-
1
=
4.096, and κ
-
2
=
1.The resulting temperature distributions are given in Figures4 and 5 for various values of location ξ and Fourier number τ. In Figure 4, the temporal variation in B
2 is also shown in the top part, which ranges from 0.2 to 1.2. In addition, Table 3 lists the temperatures at both the insulated and heated surfaces, for different number of terms in the series, to show the convergence behaviour of the present analytical solution. From Table 3, it can be seen that the convergence of the present solution is very fast. When the Fourier number is larger than 0.2, the error for the one-term approximation solution is less than 0.3%. This fast convergence of the solution is one of the notable features that the analytical solutions obtained through the shifting function method have [18, 20]. Note that all the graphical representations shown in Section 5 are based on the results obtained from the four-term approximation.Table 3
Dimensionless temperatureθ for different number of terms in the series of (54), for s = 1, p
1 = 1.2, and p
2 = 1.
τ
θ
(
0
,
τ
)
θ
(
1
,
τ
)
1 term
2 terms
4 terms
10 terms
1 term
2 terms
4 terms
10 terms
0.1
−0.0064
−0.0024
−0.0024
−0.0024
0.1116
0.1006
0.1006
0.1006
0.2
0.0110
0.0119
0.0119
0.0119
0.1590
0.1563
0.1563
0.1563
0.3
0.0344
0.0345
0.0345
0.0345
0.2051
0.2045
0.2045
0.2045
0.4
0.0622
0.0622
0.0622
0.0622
0.2498
0.2497
0.2497
0.2497
0.5
0.0934
0.0935
0.0935
0.0935
0.2928
0.2928
0.2928
0.2928
1
0.2723
0.2723
0.2723
0.2723
0.4807
0.4807
0.4807
0.4807
1.5
0.4461
0.4461
0.4461
0.4461
0.6241
0.6241
0.6241
0.6241
2
0.5903
0.5903
0.5903
0.5903
0.7300
0.7300
0.7300
0.7300
3
0.7853
0.7853
0.7853
0.7853
0.8624
0.8624
0.8624
0.8624Figure 4
Influence of parameters on the temperature variation of the bilayered plate at two locations: (a) ξ
=
0 and (b) ξ
=
1 for p
1 = 1.2 and p
2 = 1.
(a)
(b)Figure 5
Temperature distributions along the thickness of the bilayered plate atτ = 1 and 4 with different s values for p
1 = 1.2 and p
2 = 1.As shown in Figure4, the temporal variations in the temperatures correspond to the variation in B
2, being the most remarkable for s = 2, for which the rate of increase in B
2 is the highest. The time-dependent Biot number B
2 greatly affects the temperature of the heated surface for a small value of τ, having a lagged influence on the temperature at the insulated boundary. It is worth noting that the influence on the insulated boundary is not negligible. Figure 5 shows that different s values make a measurable difference in the profiles of θ throughout the thickness.


### 5.2. Three-Layered Composite with Periodically Varying HTC

Consider the case of a three-layered composite made of tin, aluminium, and lead with interfaces atξ
1 = 1/3 and ξ
2 = 2/3 [26]. Again, we assume perfect thermal contact. The initial and boundary conditions are the same as those in Section 5.1, except that B
2
τ
=
p
1
-
p
2
cos
⁡
(
ω
τ
), where ω is a constant. The following thermal conductivities and diffusivities have been used [26]: λ
-
1
=
1.900, λ
-
2
=
5.850, λ
-
3
=
1, κ
-
1
=
1.633, κ
-
2
=
4.096, and κ
-
3
=
1.Figure6 illustrates the temperature variations at two external boundaries for different ω values. The periodically varying B
2 ranges from 0.2 to 2.2, as shown in the top part of the figures. Visible oscillations are observed in both temperatures for ω = 2 and 4, for which the change cycle of B
2 is short. Moreover, unlike the case in Section 5.1, there are time slots during which the magnitude relation of the Biot numbers is not necessarily consistent with that of the temperatures for different ω values. The temperature profiles of the composite at τ = 1 and 4 are plotted in Figure 7, which are all monotonic functions of ξ with discontinuous gradient at the interfaces. It can be observed that the effects of the temporal variations in B
2 propagate through the whole thickness of the three-layered composite.Figure 6
Influence of parameterω on the temperature variation of the three-layered composite at two locations: (a) ξ = 0 and (b) ξ = 1 for p
1 = 1.2 and p
2 = 1.
(a)
(b)Figure 7
Temperature distributions along the thickness of the three-layered composite atτ = 1 and 4 with different ω values for p
1 = 1.2 and p
2 = 1.


## 5.1. Bilayered Plate with Exponentially Varying HTC

Consider the case of a bilayered plate of aluminium and lead, whose initial temperature is zero. Equal thicknesses for both layers are used, which leads toξ
1 = 1/2, with perfect thermal contact. For τ
>
0, the surrounding medium on the side of lead is raised to a temperature of Y
2 = 1, while the other surface of the plate is insulated (B
1 = 0). For the Biot number at ξ = 1, we choose B
2
τ
=
p
1
-
p
2
exp
⁡
(
-
s
τ
), where p
1, p
2, and s are arbitrary constants. The following material properties have been used in the calculations [26]: λ
-
1
=
5.850, λ
-
2
=
1, κ
-
1
=
4.096, and κ
-
2
=
1.The resulting temperature distributions are given in Figures4 and 5 for various values of location ξ and Fourier number τ. In Figure 4, the temporal variation in B
2 is also shown in the top part, which ranges from 0.2 to 1.2. In addition, Table 3 lists the temperatures at both the insulated and heated surfaces, for different number of terms in the series, to show the convergence behaviour of the present analytical solution. From Table 3, it can be seen that the convergence of the present solution is very fast. When the Fourier number is larger than 0.2, the error for the one-term approximation solution is less than 0.3%. This fast convergence of the solution is one of the notable features that the analytical solutions obtained through the shifting function method have [18, 20]. Note that all the graphical representations shown in Section 5 are based on the results obtained from the four-term approximation.Table 3
Dimensionless temperatureθ for different number of terms in the series of (54), for s = 1, p
1 = 1.2, and p
2 = 1.
τ
θ
(
0
,
τ
)
θ
(
1
,
τ
)
1 term
2 terms
4 terms
10 terms
1 term
2 terms
4 terms
10 terms
0.1
−0.0064
−0.0024
−0.0024
−0.0024
0.1116
0.1006
0.1006
0.1006
0.2
0.0110
0.0119
0.0119
0.0119
0.1590
0.1563
0.1563
0.1563
0.3
0.0344
0.0345
0.0345
0.0345
0.2051
0.2045
0.2045
0.2045
0.4
0.0622
0.0622
0.0622
0.0622
0.2498
0.2497
0.2497
0.2497
0.5
0.0934
0.0935
0.0935
0.0935
0.2928
0.2928
0.2928
0.2928
1
0.2723
0.2723
0.2723
0.2723
0.4807
0.4807
0.4807
0.4807
1.5
0.4461
0.4461
0.4461
0.4461
0.6241
0.6241
0.6241
0.6241
2
0.5903
0.5903
0.5903
0.5903
0.7300
0.7300
0.7300
0.7300
3
0.7853
0.7853
0.7853
0.7853
0.8624
0.8624
0.8624
0.8624Figure 4
Influence of parameters on the temperature variation of the bilayered plate at two locations: (a) ξ
=
0 and (b) ξ
=
1 for p
1 = 1.2 and p
2 = 1.
(a)
(b)Figure 5
Temperature distributions along the thickness of the bilayered plate atτ = 1 and 4 with different s values for p
1 = 1.2 and p
2 = 1.As shown in Figure4, the temporal variations in the temperatures correspond to the variation in B
2, being the most remarkable for s = 2, for which the rate of increase in B
2 is the highest. The time-dependent Biot number B
2 greatly affects the temperature of the heated surface for a small value of τ, having a lagged influence on the temperature at the insulated boundary. It is worth noting that the influence on the insulated boundary is not negligible. Figure 5 shows that different s values make a measurable difference in the profiles of θ throughout the thickness.


## 5.2. Three-Layered Composite with Periodically Varying HTC

Consider the case of a three-layered composite made of tin, aluminium, and lead with interfaces atξ
1 = 1/3 and ξ
2 = 2/3 [26]. Again, we assume perfect thermal contact. The initial and boundary conditions are the same as those in Section 5.1, except that B
2
τ
=
p
1
-
p
2
cos
⁡
(
ω
τ
), where ω is a constant. The following thermal conductivities and diffusivities have been used [26]: λ
-
1
=
1.900, λ
-
2
=
5.850, λ
-
3
=
1, κ
-
1
=
1.633, κ
-
2
=
4.096, and κ
-
3
=
1.Figure6 illustrates the temperature variations at two external boundaries for different ω values. The periodically varying B
2 ranges from 0.2 to 2.2, as shown in the top part of the figures. Visible oscillations are observed in both temperatures for ω = 2 and 4, for which the change cycle of B
2 is short. Moreover, unlike the case in Section 5.1, there are time slots during which the magnitude relation of the Biot numbers is not necessarily consistent with that of the temperatures for different ω values. The temperature profiles of the composite at τ = 1 and 4 are plotted in Figure 7, which are all monotonic functions of ξ with discontinuous gradient at the interfaces. It can be observed that the effects of the temporal variations in B
2 propagate through the whole thickness of the three-layered composite.Figure 6
Influence of parameterω on the temperature variation of the three-layered composite at two locations: (a) ξ = 0 and (b) ξ = 1 for p
1 = 1.2 and p
2 = 1.
(a)
(b)Figure 7
Temperature distributions along the thickness of the three-layered composite atτ = 1 and 4 with different ω values for p
1 = 1.2 and p
2 = 1.


## 6. Conclusions

An analytical method based on the shifting function method has been presented for solving the one-dimensional transient heat conduction problem with time-dependent HTC for a composite slab consisting of an arbitrary number of layers. Applying the shifting function method yields new partial differential equations under conventional types of initial, boundary, and interfacial conditions. In this paper, the derived equations were solved by using an orthogonal expansion based on the Vodicka type of orthogonality relationship. Numerical calculations were performed for two representative composite slabs. The numerical results demonstrated the effects of temporal variations in the HTC on the transient temperature field, in particular, on the temperatures at two external boundaries.This paper considers the temporal variation in the HTC at only one external boundary. The solution needs to be extended so that the time-dependency of the HTCs at both boundaries can be considered. This might possibly be achieved by introducing two different shifting functions [19, 29]. The present analytical method can easily be extended to other simple geometries; it changes nothing beyond the forms of the eigenfunction, shifting function, function L
i
(
ξ
), and the orthogonal relationship of (44). The analytical solution derived is applicable to nonhomogeneous materials, such as functionally graded materials, under the laminate approximation theory [30].The presented analytical method is available for transient heat conduction problems with arbitrary time-dependent HTC. However, currently, it cannot be applied to those in which the HTC is a function of position or of temperature.



---

*Source: 4707860-2018-04-19.xml*