Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

56, 2, pp. 339-349, Warsaw 2018
DOI: 10.15632/jtam-pl.56.2.339

FRACTIONAL HEAT CONDUCTION IN A SPHERE UNDER

MATHEMATICAL AND PHYSICAL ROBIN CONDITIONS

STANISŁAW KUKLA, Urszula Siedlecka

Czestochowa University of Technology, Institute of Mathematics, Częstochowa, Poland

e-mail: stanislaw.kukla@im.pcz.pl; urszula.siedlecka@im.pcz.pl

In this paper, the effect of a fractional order of time-derivatives occurring in fractional heat
conduction models on the temperature distribution in a composite sphere is investigated.
The researchconcernsheat conduction in a sphere consistingof a solid sphere anda spherical
layerwhichare inperfect thermal contact.The solutionof the problemwith a classicalRobin
boundary condition and continuity conditions at the interface in an analytical form has
been derived. The fractional heat conduction is governed by the heat conduction equation
with the Caputo time-derivative, a Robin boundary condition and a heat flux continuity
condition with the Riemann-Liouville derivative. The solution of the problem of non-local
heat conduction by using the Laplace transform technique has been determined, and the
temperature distribution in the sphere by using a method of numerical inversion of the
Laplace transforms has been obtained.

Keywords: heat conduction, fractional heat equation, Robin boundary condition

1. Introduction

The classical heat conduction model based on the Fourier law has a non-physical property that
the heat propagates with an infinite speed (Özişik, 1993). This property is a consequence of
the dependence between the heat flux vector and the temperature gradient which is established
by the Fourier law. This disadvantage does not appear when the non-local time dependence
between the flux vector and the temperature gradient is assumed (Povstenko, 2014; Sur and
Kanoria, 2014). This assumption leads to a differential equation and/or boundary conditions
with derivatives of a non-integer order. The properties of fractional derivatives and different
analytical methods to solve fractional differential equations are presented in (Atanacković et
al., 2014; Klimek, 2009; Leszczyński, 2011; Magin, 2006; Mainardi, 2010; Povstenko, 2015).
Approximate numerical methods were applied to solving fractional initial-boundary problems
in numerous papers, for example in (Blaszczyk and Ciesielski, 2017; Ciesielski and Błaszczyk,
2013; Dimitrov, 2014).

The heat conductionmodelled by using the fractional order derivative is the subject ofmany
papers. A mathematical model of one-dimensional heat conduction in a slab was proposed in
paper (Žecová and Terpák, 2015). The Grünwald-Letnikov derivative with respect to a time
variable was used. A solution to the problem of fractional heat conduction in a two-layered
slab with the Caputo time-derivative in the heat conduction equation was presented in (Kukla
andSiedlecka, 2015). Heat transfer for non-contacting face seals described by the time-fractional
heat conduction equation in the cylindrical coordinate systemwas considered in (Blasiak, 2016).
The fractional model of thermal energy transport in rigid bodies was derived in (Raslan, 2016).
The effect of the order of the Marchand-type derivative in the heat transfer equation on the
temperature distribution in a rigid conductor was numerically investigated. An application of
the fractional order theory to a problem of thermal stress distribution in a spherical shell was



340 S. Kukla, U. Siedlecka

studied in (Zingales, 2014). In the paper by Atangana and Bildik (2013), the time fractional
calculuswas employed in themathematicalmodel of groundwater flow.Applications of fractional
order systems to an ultracapacitor and beam heating problems were presented in (Dzieliński
et al., 2010). An application of fractional calculus in continuum mechanics to a problem of
linear elasticity under small deformation was shown in (Sumelka and Blaszczyk, 2014). Some
applications of the fractional calculus were also discussed in the papers (Abbas, 2012; Dalir and
Bashour, 2010; Rahimy, 2010).

Solutions to time-fractional heat conduction problems in a spherical coordinate system are
presented inmany papers. In the paper byNing and Jiang (2011), for the problem of fractional
heat conduction in a sphere, the method of the Laplace transform and the variable separation
were used. An analytical solution to the problem of the time-fractional radial heat conduction
in a multilayered sphere under the Robin boundary condition was presented by Kukla and
Siedlecka (2017). Fundamental solutions to the Cauchy problem and to the source problem of
the heat conduction fractional equation in a spherical coordinate system in an analytical form
were derived by article Povstenko and Klekot (2017).

The fractional heat conduction equation is complemented by initial andboundaryconditions.
Mathematical and physical formulations of the initial and boundary conditions can be conside-
red in fractional heat conduction models (Povstenko, 2013). The mathematical formulations of
Dirichlet, Neumann andRobin boundary conditions are the same as these in the classical theory
of heat conduction. Also, the physical Dirichlet condition has the same form as the classical bo-
undary condition of thefirst kind,while the physicalNeumannand the physicalRobinboundary
conditions contain the fractional time-derivative. If two solids are in perfect thermal contact,
the physical formulation of the condition of heat flux equality through the contact surface also
contain the fractional time-derivative (Povstenko, 2013).

The solution to the problem of linear fractional heat conduction in a sphere under ma-
thematical boundary conditions can be determined in an analytical form. However, in solving
such problems of heat conduction under physical Neumann or Robin boundary conditions, an
approximate methods must be used. Application of the Laplace transform method to a linear
problem allows one to obtain a solution in the Laplace domain. For the fractional heat conduc-
tion problems under physical Neumann or Robin boundary conditions and physical continuity
conditions, the inverse Laplace transform in an analytical form can not be determined. The so-
lution to the problem is obtained by applying numerical inversion of the Laplace transform.The
methods for numerical inversion of the Laplace transform, which are used in the classical pro-
blems, can be also applied to the Laplace transform obtained in fractional analysis. A review of
themethods to numerical inversion of the Laplace transformwas presented byKuhlman (2013).
An application of selected methods to determine the inverse Laplace transform in fractional
calculus were presented in (Brzeziński andOstalczyk, 2016; Sheng et al., 2011). Modification of
the method introduced by Gaver (1966) was presented in (Abate and Valkó, 2004; Valkó and
Abate, 2004).

In this paper, the fractional heat conduction problem in a solid sphere under the mathe-
matical and physical boundary condition is studied. The considered sphere consists of an inner
sphere and a spherical layer. We assume perfect thermal contact of the inner sphere and the
spherical layer which is modeled by mathematical or physical conditions. The exact solution of
the problem for the mathematical boundary condition and the solution in the Laplace domain
for the physical formulation of boundary and continuity conditions are presented. The effect of
the order of theRiemman-Liouville derivative in theRobin physical condition and in the contact
condition at the interface on the temperature distribution in the sphere has been numerically
investigated.



Fractional heat conduction in a sphere under mathematical... 341

2. Formulation of the problem

Weconsider theproblemofheat conduction ina spherewhich consists of a solid sphereoccupying
the region 0¬ r ¬ r1 and a spherical layer defined by r1 ¬ r ¬ b, in the spherical coordinates
system. The time-fractional heat conduction in the inner sphere (i = 1) and in the spherical
layer (i=2) is governed by the following equation

1

r2
∂

∂r

(
r2
∂Ti
∂r

)
=
1

ai

∂αiTi
∂tαi

i=1,2 (2.1)

where ai is the thermal diffusivity, λi is the thermal conductivity and αi denotes the fractional
order of the left Caputo derivative with respect to time t. The Caputo derivative is defined by
(Podlubny, 1999)

C
aD
α
t f(t)=

dαf(t)

dtα
=

1

Γ(m−α)

t∫

a

(t− τ)m−α−1
dmf(τ)

dτm
dτ m−1<α<m (2.2)

We consider the case of a=0andα∈ (0,1]. Note, that the thermal diffusivity coefficient can be
interpreted as ameasure of the distance on which the thermal front propagates in amedium at
the given time. The thermal conductivity is a measure of the ability of the medium to transfer
heat.

Thecondition at thecentre of the sphere, the continuity conditions at the interface, theRobin
boundary condition on the outer surface and the initial condition are (Povstenko, 2013a,b)

|T(0, t)|<∞ T1(r1, t)=T2(r1, t) (2.3)

λ1D
1−β1
RL

∂T1
∂r
(r1, t)=λ2D

1−β2
RL

∂T2
∂r
(r1, t)

λ2D
1−β2
RL

∂T2
∂r
(b,t)= a∞

(
T∞(t)−T2(b,t)

) (2.4)

T(r,0)=Fi(r) (2.5)

where a∞ is the outer heat transfer coefficient and T∞ is the ambient temperature. The left
Riemann-Liouville fractional derivativeD

1−β
RL which occurs in equations (2.4) is defined by (Die-

thelm, 2010)

D
β
RLf(t)=

d

dt

( 1
Γ(1−β)

t∫

0

f(τ)

(t− τ)β
dτ
)

0<β¬ 1 (2.6)

Conditions (2.4) forβ1 =α1 andβ2 =α2 forα1,α2 ∈ (0,1) are called the physical conditions
(Rahimt, 2010; Raslan, 2016). If β1 = β2 = 1, the conditions are called the mathematical
conditions. In this case, the D0RL means an identity operator and can be omitted in equations
(2.4).

3. Solution to the problem

The fractional heat conduction problem defined by equations (2.1) and (2.3)-(2.5) can be trans-
formed to a new problem for functionsUi(r,t) by using the formula

Ui(r,t)= r (Ti(r,t)−T∞(t)) i=1,2 (3.1)



342 S. Kukla, U. Siedlecka

Taking into account relationship (3.1) in equation (2.1) and conditions (2.3)-(2.5), we obtain a
formulation of the initial-boundary problem in the form

ai
∂2Ui(r,t)

∂r2
=
∂αiUi(r,t)

∂tαi
+r
dαiT∞(t)

dtαi
i=1,2 (3.2)

U1(0, t)= 0 U1(r1, t)=U2(r1, t) (3.3)

λ1D
1−β1
RL

(∂U1(r1, t)
∂r

−
1

r1
U1(r1, t)

)
=λ2D

1−β2
RL

(∂U2(r1, t)
∂r

−
1

r1
U2(r1, t)

)

λ2D
1−β2
RL

(∂U2
∂r
(b,t)− 1

b
U2(b,t)

)
=−a∞U2(b,t)

(3.4)

Ui(r,0)= r (Fi(r)−T∞(0)) i=1,2 (3.5)
The solution to initial-boundaryproblem (3.2)-(3.5) forβ1 =β2 =1 (mathematical formulation)
and for β1 =α1, β2 =α2 (physical formulation) will be presented below.

3.1. Mathematical conditions

An analytical solution to time-fractional heat conduction problem (3.2)-(3.5) under mathe-
matical conditions (3.4) for α1 =α2 =αwill be determined by using the method of separation
of variables. As a result, we find a solution to the problem in the form of a series

Ui(r,t)=
∞∑

k=1

Λk(t)Φi,k(r) i=1,2 (3.6)

The functions Φ1,k(r) and Φ2,k(r) for k = 1,2, . . . are obtained as a solution to the corre-
sponding eigenvalue problem

d2Φi,k(r)

dr2
+
γ2k
ai
Φi,k(r)= 0 i=1,2 (3.7)

Φ1,k(0)= 0 Φ1,k(r1)=Φ2,k(r1) (3.8)

λ1
dΦ1(r1)

dr
+
1

r1
(λ2−λ1)Φ1(r1)=λ2

dΦ2(r1)

dr

dΦ2(b)

dr
=
(1
b
−
a∞
λ2

)
Φ2(b) (3.9)

The eigenfunctions Φi,k(r) are given by

Φ1,k(r)=B1,k sinµ1,kr Φ2,k(r)=A2,k cosµ2,k(r−r1)+B2,k sinµ2,k(r−r1) (3.10)

where µi,k = γk/
√
ai and γk are the roots of the eigenvalue equation

M1λ2µ1 sinµ1r1+M2M3 =0 (3.11)

where

M1 =
(a∞
λ2
− 1
b

)
cosµ2(b−r1)−µ2 sinµ2(b−r1)

M2 =
(a∞
λ2
−
1

b

)
sinµ2(b−r1)+µ2cosµ2(b−r1)

M3 =λ1µ1cosµ1r1+
λ2−λ1
r1

sinµ1r1

The coefficients B1,k, A2,k and B2,k are determined by using continuity and boundary con-
ditions (3.8) and (3.9). AssumingB1,k =1, we obtain A2,k =sinµ1,kr1 andB2,k =M3/λ2µ2,k.



Fractional heat conduction in a sphere under mathematical... 343

The functionΛk(t), occurring in equation (3.6), is a solution to the fractional initial problem
which is obtained by using the orthogonality condition in the form

λ1
a1

r1∫

0

Φ1,k(r)Φ1,k′(r) dr+
λ2
a2

b∫

r1

Φ2,k(r)Φ2,k′(r) dr=

{
0 for k′ 6= k
Nk for k

′ = k
(3.12)

Assuming Fi(r) = Tinit = const for i = 1,2 and condition (3.12) in equation (3.2) and (3.5),
the initial problem is obtained

dαΛk(t)

dtα
+γ2kΛk(t)= 0

Λk(0)=
Tinit−T∞
Nrk

(λ1
a1

r1∫

0

rΦ1,k(r) dr+
λ2
a2

b∫

r1

rΦ2,k(r) dr
) (3.13)

A solution to problem (3.13) is given by (Diethelm, 2010)

Λk(t)=
Tinit−T∞
Nrk

Eα(−γ2ktα)
(λ1
a1

r1∫

0

rΦ1,k(r) dr+
λ2
a2

b∫

r1

rΦ2,k(r) dr
)

(3.14)

whereEα(z) is theMittag-Leffler function (Kilbas et al., 2006)

Eα(z)=
∞∑

k=0

zk

Γ(αk+1)
(3.15)

Finally, the functionsTi(r,t) are given by equations (3.1), (3.6), (3.10) and (3.14). Assuming
that the following conditions are fulfilled: a1 = a2 = a, λ1 = λ2 = λ, α1 = α2 = α and
β1 =β2 =1, we obtain the temperature T(r,t) in the homogeneous sphere

T(r,t) =T∞+
4(Tinit−T∞)

r

∞∑

k=1

bµk cosbµk− sinbµk
µk(sin2bµk−2bµk)

Eα(−γ2ktα)sinµkr (3.16)

In this case, µk = γk/
√
a and γk are the roots of the equation

(
1−
ba∞
λ

)
sinbµ− bµcosbµ=0 (3.17)

3.2. Physical conditions

We obtain a solution to problem (3.2)-(3.5) under physical boundary and continuity condi-
tions (β1 = α1, β2 = α2 in equations (3.4) and (3.5)) by using the Laplace transform method.
The Laplace transform f(s) of a function f(t) is defined by

f(s)=

∞∫

0

f(t)e−st dt (3.18)

where s is a complex parameter. Using the properties of the Laplace transform, equations
(3.2)-(3.4) can be rewritten in the Laplace domain as

d2Ui
dr2
−
sαi

ai
Ui(r,s)=

rsαi

ai

(
T∞(s)−

Fi(r)

s

)
(3.19)

U1(0,s)= 0 U1(r1,s)=U2(r1,s) (3.20)



344 S. Kukla, U. Siedlecka

λ1s
1−α1

(dU1(r1,s)
dr

−
1

r1
U1(r1,s)

)
=λ2s

1−α2
(dU2(r1,s)

dr
−
1

r1
U2(r1,s)

)

λ2s
1−α2

(dU2(b,s)
dr

− 1
b
U2(b,s)

)
=−a∞U2(b,s)

(3.21)

The general solution to equation (3.19) for i=1,2 has the form

U1(r,s)=B1 sinhS1r+
1

S1

r∫

0

P1(u)sinhS1(r−u) du

U2(r,s)=A2coshS2(r−r1)+B2 sinhS2(r−r1)+
1

S2

r∫

r1

P2(u)sinhS2(r−u) du
(3.22)

where

Si =
sαi/2
√
ai

Pi(r)=
rsαi

ai

(
T∞(s)−

Fi(r)

s

)

The constants B1, A2 and B2 are determined by using conditions (3.20)2 and (3.21). After
some transformations, the functionsU1(r,s) andU2(r,s) can be written as

U1(r,s)= B̃1 sinhS1r U2(r,s) = Ã2coshS2(r−r1)+ B̃2 sinhS2(r−r1) (3.23)

where

B̃1 =−sα1−2
a∞
λ1d
S2b
2 Ã2 =−sα1−2

a∞b

λ1d
S2bsinhS1r1

B̃2 = s
α2−2
b

r1

a∞b

λ2d

[(
1−sα1−α2λ2

λ1

)
sinhS1r1−S1r1coshS1r1

]

d= sα1−α2
λ2
λ1
S2r1 sinhS1r1[wcoshS2(b−r1)+S2bsinhS2(b−r1)]

+
]
S1r1coshS1r1−

(
1−sα1−α2

λ2
λ1

)
sinhS1r1

]
[S2bcoshS2(b−r1)+wsinhS2(b−r1)]

w=
a∞b

λ2s1−α2
−1

Assuming Fi(r) = Tinit = const for i = 1,2, the temperature distribution in the sphere is
given by

Ti(r,t)=T∞+(Tinit−T∞)
r1
r
L−1[Ui(r,s)] (3.24)

For the homogeneous sphere, the following conditions are fulfilled a1 = a2 = a, λ1 =λ2 =λ,
α1 =α2 =α, β1 = β2 =β and S1 =S2 =S. In this case, the function Ti(r,t) =T(r,t) has the
form

T(r,t) =T∞+(Tinit−T∞)
b

r
L−1[U(r,s)] (3.25)

where

U(r,s) =− 1+w
s(SbcoshSb+wsinhSb)

sinhSr

The inverse of the Laplace transform of the functions U1(r,s) and U2(r,s) are numerically
determined. The calculation has been performed by the Gaver method using the sequence of



Fractional heat conduction in a sphere under mathematical... 345

functionals presented in Gaver (1966) and Valkó and Abate (2004). Applying this method, the
approximate values of the original function Ui(r,t) = L

−1[Ui(r,s)] are determined using the
formula

Ui(r,t)≃nτ
(
2n

n

)
n∑

i=0

(−1)i
(
n

i

)
Ui(r,(n+ i)τ) (3.26)

where τ =(ln2)/t and n is a fixed positive integer number.

The functions Ti(r,t) and T(r,t) obtained for the mathematical and physical conditions
will serve for investigation of the influence of the orders of the Caputo and Riemann-Liouville
derivatives occurring in the heat conduction models on the temperature distribution in the
sphere.

4. Results of numerical calculations

The effect of the order of the fractional derivative in the heat conduction equation on the tempe-
rature distribution in the sphere has been numerically investigated. The results for the mathe-
matical boundary condition obtained by usingnumerical inversion of theLaplace transformshas
been comparedwith the exact solution. The computations were performed for the homogeneous
sphere (SphereA) and for the sphere consisting of a solid sphereanda spherical layer (SphereB).
The radius of both Spheres was b=1.0m and the interface in Sphere Bwas at r̂1 = r1/b=0.9.
The thermal diffusivitya=3.352·10−6m2/sα and the thermal conductivityλ=16W/(m·K)we-
re assumed for SphereA.The thermal diffusivities a1 =2.3·10−5m2/sα, a2 =3.352·10−6m2/sα
and the thermal conductivitiesλ1 =80W/(m·K),λ2 =16W/(m·K)were assumed for SphereB.
Subscript 1 was used for the inner sphere and subscript 2 – for the spherical layer of Sphere B.
For both Spheres, the outer heat transfer coefficient was a∞ = 200W/(m

2·K), the ambient
temperature was T∞ =100

◦C and the initial temperature was assumed as Tinit =25
◦C.

In Table 1, the non-dimensional temperature T̂ =T/Tinit in SphereA for different orders of
the Caputo derivative α at the reference time t̂= tb2/a=1.0 is presented. The calculation has
been performed for the mathematical Robin boundary condition, i.e. for β = 1.0. The results
were obtained by using the exact solution, Eq. (3.16), and by the Gaver method of numerical
inversionof theLaplace transforms,Eq. (3.26), andusingrelationship (3.1).Asimilar comparison
of numerically obtained non-dimensional temperatures have been performed for Sphere B. The
results are presented in Table 2. The relative error evaluated on the basis of the results given in
Tables 1 and 2 fulfils the condition: |Exact −NILT|/Exact < 3.6 ·10−5. The good accordance
of the results obtained for mathematical formulation of the boundary and continuity condition
allows one to use the NILTmethod to the heat conduction problem under physical formulation
of the boundary and continuity condition.

Table 1.Non-dimensional temperature T̂(r̂, t̂) for t̂=1.0, computed byusing the exact solution
and by using numerical inversion of the Laplace transform (NILT) for Sphere A

r̂
α=0.8 α=0.9 α=1.0

Exact NILT Exact NILT Exact NILT

0 1.12416 1.12412 2.41882 2.41882 3.91149 3.91145

0.25 1.19635 1.19634 2.52818 2.52818 3.91902 3.91900

0.50 1.49058 1.49058 2.84104 2.84103 3.93932 3.93932

0.75 2.21186 2.21186 3.30678 3.30678 3.96637 3.96629

1.00 3.50777 3.50776 3.83398 3.83398 3.99254 3.99251



346 S. Kukla, U. Siedlecka

Table 2.Non-dimensional temperature T̂(r̂, t̂) for t̂=1.0, computed byusing the exact solution
and by using numerical inversion of the Laplace transform (NILT) for Sphere B

r̂
α=0.8 α=0.9 α=1.0

Exact NILT Exact NILT Exact NILT

0 2.35412 2.35407 3.66988 3.66988 3.99982 3.99969

0.25 2.39527 2.39526 3.67879 3.67880 3.99983 3.99969

0.50 2.51791 2.51791 3.70498 3.70498 3.99984 3.99971

0.75 2.71870 2.71870 3.74679 3.74679 3.99987 3.99974

1.00 3.48174 3.48172 3.89885 3.89885 3.99995 3.99988

Fig. 1. Non-dimensional temperature T̂(r̂, t̂) as a function of time t̂ in Sphere A for various values of
fractional derivativesα and β: (a) r̂=0, (b) r̂=0.6, (c) r̂=0.8, (d) r̂=1.0



Fractional heat conduction in a sphere under mathematical... 347

Thenon-dimensional temperatures T̂ as functions of the time t̂ for various radial coordinates
are presented inFig. 1. The pairs of curves obtained formathematical and physical formulations
of theRobin condition show the effect of the order of theRiemann-Liouville derivative occurring
in the physical boundary condition on the temperature in the sphere. Significant differences can
be observed in the temperatures obtained for the classical heat conduction model (α= β =1)
and fractional models (α=0.8 and α=0.9), particularly in the inner points of the sphere.

The curves presented in Fig. 2 represent the non-dimensional temperatures T̂ as functions
of the reference time t̂ for SphereB. In numerical calculations with themathematical conditions
(MC) the following values have been assumed α1 = α2 = α= 0.8, 0.9, 1.0 and β1 = β2 = 1.0.
The numerical calculations to the problemwith the physical conditions (PC) have been carried
out for:α1 =β1 =0.9,α2 =α=0.8, 0.9, 1.0 andβ2 =β=α2.Ahigher temperature is observed
for the heat conductionwith the physical boundary and continuity conditions than in themodel
with the mathematical formulation of these conditions. A significant effect on the temperature
distribution in the sphere has the order of the Caputo derivative in the heat conductionmodel.

Fig. 2. Non-dimensional temperature T̂(r̂, t̂) as a function of the time t̂ in Sphere B for various values of
fractional derivativesα and β: (a) r̂=0, (b) r̂=0.6, (c) r̂=0.8, (d) r̂=1.0



348 S. Kukla, U. Siedlecka

5. Conclusions

Asolution to theproblemof fractional heat conduction inahomogeneous sphereandacomposite
sphere consistingwith a solid sphere and a spherical layer has beenpresented.Themathematical
and physical formulations of theRobin boundary condition and the continuity conditions at the
interface have been considered. The temperature distribution in the sphere, under the physical
boundary and continuity conditions, has been obtained by using the Laplace transform techni-
que. Numerical results show a significant effect of the order of the Caputo derivative occurring
in the heat conduction equation on the time-history of temperature in the sphere. The order
of the Riemann-Liouville derivative occurring in the boundary and continuity conditions of the
fractional model of heat conduction has a smaller effect on the time-history of temperature in
the sphere than the order of the fractional Caputo derivative in the heat conduction equation.

References

1. Abate J., Valkó P.P., 2004,Multi-precision Laplace transform inversion, International Journal
for Numerical Methods in Engineering, 60, 979-993

2. Abbas I.A., 2012, On a thermoelastic fractional order model, Journal of Physics, 1, 2, 24-30

3. Atanacković T.M., Pilipović S., Stanković B., Zorica D., 2014, Fractional Calculus with
Applications in Mechanics, JohnWiley & Sons, NewYork

4. AtanganaA., BildikN., 2013,The use of fractional order derivative to predict the groundwater
flow,Mathematical Problems in Engineering, Article ID 543026

5. Blasiak S., 2016, Time-fractional heat transfer equations in modeling of the non-contacting face
seals, International Journal of Heat and Mass Transfer, 100, 79-88

6. BlaszczykT., CiesielskiM., 2017,Numerical solution ofEuler-Lagrange equationwithCaputo
derivatives,Advances in Applied Mathematics and Mechanics, 9, 173-185

7. Brzeziński D.W., Ostalczyk P., 2016, Numerical calculations accuracy comparison of the in-
verse Laplace transformalgorithms for solution of fractional order differential equations,Nonlinear
Dynamics, 84, 1, 65-77

8. CiesielskiM., BłaszczykT., 2013,An approximationof the analytical solution of the fractional
Euler-Lagrange equation, Journal of Applied Mathematics and Computational Mechanics, 12, 4,
23-30

9. DalirM.,BashourM., 2010,Applications of fractional calculus,AppliedMathematical Sciences,
4, 21, 1021-1032

10. Diethelm K., 2010, The Analysis of Fractional Differential Equations, Springer-Verlag Berlin
Heidelberg

11. Dimitrov Y., 2014, Numerical approximations for fractional differential equations, Journal of
Fractional Calculus and Applications, 5, 38, 1-45

12. Dzieliński A., Sierociuk D., Sarwas G., 2010, Some applications of fractional order calculus,
Bulletin of the Polish Academy of Sciences – Technical Sciences, 58, 4, 583-592

13. GaverD.P. Jr., 1966,Observing stochastic processes and approximate transform inversion,Ope-
rational Research, 14, 444-459

14. KilbasA., SrivastavaH.,Trujillo J., 2006,Theory andApplications of Fractional Differential
Equations, Elsevier, Amsterdam

15. Klimek M., 2009,On Solutions of Linear Fractional Differential Equations of a Variational Type,
The Publishing Office of CzestochowaUniversity of Technology, Czestochowa



Fractional heat conduction in a sphere under mathematical... 349

16. KuhlmanK.L., 2013,Reviewof inverseLaplace transformalgorithms forLaplace-spacenumerical
approaches,Numerical Algorithms, 63, 2, 339-355

17. Kukla S., SiedleckaU., 2015, Laplace transform solution of the problemof time-fractional heat
conduction in a two-layered slab, Journal of Applied Mathematics and Computational Mechanics,
14, 4, 105-113

18. Kukla S., Siedlecka U., 2017, An analytical solution to the problem of time-fractional heat
conduction in a composite sphere,Bulletin of the Polish Academy of Sciences – Technical Sciences,
65, 2, 179-186

19. Leszczyński J.S., 2011,An Introduction to Fractional Mechanics, The Publishing Office of Cze-
stochowaUniversity of Technology, Czestochowa

20. Magin R.L., 2006,Fractional Calculus in Bioengineering, Begell House Publishers, Connecticut

21. Mainardi F., 2010, Fractional Calculus and Waves in Linear Viscoelasticity: an Introduction to
Mathematical Models, Imperial College Press, London

22. Ning T. H., Jiang X. Y., 2011, Analytical solution for the time-fractional heat conduction
equation in spherical coordinate system by the method of variable separation, Acta Mechanica
Sinica, 27, 6, 994-1000

23. Özişik M.N., 1993,Heat Conduction, Wiley, NewYork

24. Podlubny I., 1999,Fractional Differential Equations, Academic Press, San Diego

25. PovstenkoY., 2013a,Fractionalheat conduction inan infinitemediumwitha spherical inclusion,
Entropy, 15, 4122-4133

26. Povstenko Y., 2013b, Time-fractional heat conduction in an infinite medium with a spherical
hole under Robin boundary condition, Fractional Calculus and Applied Analysis, 16, 2, 354-369

27. Povstenko Y., 2014, Fractional heat conduction in a semi-infinite composite body,Communica-
tions in Applied and Industrial Mathematics, 6, 1, e-482

28. Povstenko Y., 2015, Linear Fractional Diffusion-Wave Equation for Scientists and Engineers,
Birkhäuser, NewYork

29. Povstenko Y., Klekot J., 2017, The fundamental solutions to the central symmetric time-
fractional heat conduction equation with heat absorption, Journal of Applied Mathematics and
Computational Mechanics, 16, 2, 101-112

30. RahimyM., 2010,Applications of fractional differential equations,AppliedMathematical Sciences,
4, 50, 2453-2461

31. Raslan W.E., 2016, Application of fractional order theory of thermoelasticity to a 1D problem
for a spherical shell, Journal of Theoretical and Applied Mechanics, 54, 1, 295-304

32. Sheng H., Li Y., Chen Y., 2011,Application of numerical inverse Laplace transform algorithms
in fractional calculus, Journal of the Franklin Institute, 384, 315-330

33. Sumelka W., Blaszczyk T., 2014, Fractional continua for linear elasticity, Archives of Mecha-
nics, 66, 3, 147-172

34. Sur A., Kanoria M., 2014, Fractional heat conductionwith finite wave speed in a thermo-visco-
elastic spherical shell, Latin American Journal of Solids and Structures, 11, 1132-1162

35. Valkó P.P., Abate J., 2004, Comparison of sequence accelerators for the Gaver method
of numerical Laplace transform inversion, Computers and Mathematics with Applications, 48,
629-636

36. Žecová M., Terpák J., 2015, Heat conduction modeling by using fractional-order derivatives,
Applied Mathematics and Computation, 257, 365-373

37. Zingales M., 2014, Fractional-order theory of heat transport in rigid bodies,Communications in
Nonlinear Science and Numerical Simulation, 19, 3938-3953

Manuscript received February 14, 2017; accepted for print April 18, 2017