# Axisymmetric Solutions to Time-Fractional Heat Conduction Equation in a Half-Space under Robin Boundary Conditions

Authors: Y. Z. Povstenko
Journal: International Journal of Differential Equations
 (2012)

---

## Abstract

The time-fractional heat conduction equation with the Caputo derivative of the order0<α<2 is considered in a half-space in axisymmetric case under two types of Robin
boundary condition: the mathematical one with the prescribed linear combination of
the values of temperature and the values of its normal derivative and the physical
condition with the prescribed linear combination of the values of temperature and the
values of the heat flux at the boundary.


---

## Body

## 1.   Introduction

The generalized Fourier law, the time-nonlocal dependence between the heat flux vectorq and the temperature gradient gradT with the “long-tail memory” power kernel [1, 2] (see also [3]) as (1.1)q(t)=-kΓ(α)∂∂t∫0t(t-τ)α-1gradT(τ)dτ,0<α≤1;q(t)=-kΓ(α-1)∫0t(t-τ)α-2gradT(τ)dτ,1<α≤2,
where Γ(α) is the gamma function, can be interpreted in terms of the fractional calculus:
(1.2)q(t)=-kDRL1-αgradT(t),0<α≤1,q(t)=-kIα-1gradT(t),1<α≤2,
and in combination with the law of conservation of energy as
(1.3)ρmCdTdt=-divq
results in the time-fractional heat conduction equation with the Caputo fractional derivative
(1.4)∂αT∂tα=aΔT,0<α≤2.
Here, ρm is the mass density, C denotes the specific heat capacity, and a=k/(ρmC) is the thermal diffusivity coefficient.Recall that the Riemann-Liouville fractional integralIαf(t) and derivative DRLαf(t) are defined as follows (see [4–6]):
(1.5)Iαf(t)=1Γ(α)∫0t(t-τ)α-1f(τ)dτ,α>0,(1.6)DRLαf(t)=dmdtm[1Γ(m-α)∫0t(t-τ)m-α-1f(τ)dτ],m-1<α<m,
whereas the Caputo fractional derivative has the following form [5–7]:
(1.7)dαf(t)dtα=1Γ(m-α)∫0t(t-τ)m-α-1dmf(τ)dτmdτ,m-1<α<m.A detailed explanation of derivation of time-fractional heat conduction equation (1.4) from the constitutive equations (1.2) and the law of conservation of energy (1.3) can be found in [8]. Here we briefly present the main idea. In the case 0<α≤1, as a consequence of (1.2), (1.3), and (1.6), we have
(1.8)∂T∂t=a∂∂t[1Γ(α)∫0t(t-τ)α-1ΔT(τ)dτ],
or after integration with respect to time as
(1.9)T(t)-T(0)=aIαΔT.
Applying to both sides of (1.9) the Caputo derivative ∂α/∂tα, and taking into account that for α>0 [6],
(1.10)∂α∂tαIαT(t)=T(t),
we obtain the time-fractional heat conduction equation (1.4) for  0<α≤1.Similarly, for1<α≤2,  we get
(1.11)∂T∂t=aIα-1ΔT.
Applying  ∂α-1/∂tα-1 to  both sides of (1.11) gives (1.4)  for  1<α≤2 as the following equality fulfills [5]
(1.12)∂α-1∂tα-1∂T∂t=∂αT∂tα.If the heat conduction equation is investigated in a bounded domain, the boundary conditions should be prescribed. The mathematical Robin boundary condition is a specification of a linear combination of the values of temperature and the values of its normal derivative at the boundary of the considered domain(1.13)(c1T+c2∂T∂n)|S=F0(xS,t),
with some nonzero constants c1 and c2, while the physical Robin boundary condition specifies a linear combination of the values of temperature and the values of the heat flux at the boundary of the domain. For example, the condition of convective heat exchange between a body and the environment with the temperature Te(1.14)q·n|S=h(T|S-Te),
where h is the convective heat transfer coefficient, leads to
(1.15)(hT+kDRL1-α∂T∂n)|S=hTe(xS,t),0<α≤1,(hT+kIα-1∂T∂n)|S=hTe(xS,t),1<α≤2.The literature on mathematical aspects concerning correctness of initial-boundary-value problems for time-fractional diffusion equation and form and properties of its solutions is quite extensive (see, e.g., [9–16], among others). Geometrical explanation of fractional calculus is given in [17–19].There are only a few papers [20, 21] in which the fractional diffusion is investigated under the mathematical Robin boundary condition. In previous publications, problems for a cylinder [22] and a sphere [23] under mathematical and physical Neumann boundary conditions were considered. In the present paper, for the first time,  the solutions to time-fractional heat conduction equation in a half-space are studied under both the mathematical and physical Robin boundary conditions. The Laplace integral transform with respect to time t, the Hankel transform with respect to the spatial coordinate r, and the sin-cos-Fourier transforms with respect to spatial coordinate z are used. The solutions under the mathematical and physical Neumann boundary conditions are obtained as particular cases.


## 2. Mathematical Preliminaries

### 2.1. Laplace Transform

The Laplace transform is defined as(2.1)ℒ{f(t)}=f*(s)=∫0∞‍f(t)e-stdt,
where s is the transform variable.The inverse Laplace transfrom is carried out according to the Fourier-Mellin formula:(2.2)ℒ-1{f*(s)}=f(t)=12πi∫c-i∞c+i∞f*(s)estds,t>0,
where c is a positive fixed number.The Laplace transform of the Riemann-Liouville fractional integral of the orderα is carried out according to the formula similar to the Laplace transform of m-fold primitive of a function as
(2.3)ℒ{Iαf(t)}=1sαf*(s).The Caputo derivative for its Laplace transform requires the knowledge of the initial values of the functionf(0+) and its integer derivatives f(k)(0+) of the order k=1,2,…,m-1(2.4)ℒ{dαf(t)dtα}=sαf*(s)-∑k=0m-1f(k)(0+)sα-1-k,m-1<α<m,
whereas the Riemann-Liouville derivative for its Laplace transform rule requires the knowledge of the initial values of the fractional integral Im-αf(0+) and its derivatives of the order k=1,2,…,m-1(2.5)ℒ{DRLαf(t)}=sαf*(s)-∑k=0m-1DkIm-αf(0+)sm-1-k,m-1<α<m.
The reader interested in applications of integral transforms in fractional calculus is referred to [4–7, 24].


### 2.2. Hankel Transform

The Hankel transform is used to solve problem in cylindrical coordinates in the domain0≤r<∞ and is defined as
(2.6)ℋ{f(r)}=f*(η)=∫0∞f(r)Jν(ηr)rdr,ℋ-1{f*(η)}=f(r)=∫0∞f*(η)Jν(ηr)ηdη,
where Jν(r) is the Bessel function of the order ν.The following formula is fulfilled:(2.7)ℋ{d2f(r)dr2+1rdf(r)dr-ν2r2f(r)}=-η2f*(η).


### 2.3. Sin-Cos-Fourier Transform

In the case of boundary condition of the third kind with the prescribed boundary value of linear combination of a function and its normal derivative(2.8)z=0࢙:-dfdz+Hf=φ0,
the  following sin-cos-Fourier transform [25] is employed:
(2.9)ℱ{f(z)}=f*(ξ)=∫0∞K(z,ξ)f(z)dz,ℱ-1{f*(ξ)}=f(z)=2π∫0∞K(z,ξ)f*(ξ)dξ,
with the kernel
(2.10)K(z,ξ)=ξcos(zξ)+Hsin(zξ)ξ2+H2.Application of the sin-cos-Fourier transform to the second derivative of a function gives(2.11)ℱ{d2f(z)dz2}=-ξ2f*(ξ)-ξξ2+H2[df(z)dz-Hf(z)]|z=0.


## 2.1. Laplace Transform

The Laplace transform is defined as(2.1)ℒ{f(t)}=f*(s)=∫0∞‍f(t)e-stdt,
where s is the transform variable.The inverse Laplace transfrom is carried out according to the Fourier-Mellin formula:(2.2)ℒ-1{f*(s)}=f(t)=12πi∫c-i∞c+i∞f*(s)estds,t>0,
where c is a positive fixed number.The Laplace transform of the Riemann-Liouville fractional integral of the orderα is carried out according to the formula similar to the Laplace transform of m-fold primitive of a function as
(2.3)ℒ{Iαf(t)}=1sαf*(s).The Caputo derivative for its Laplace transform requires the knowledge of the initial values of the functionf(0+) and its integer derivatives f(k)(0+) of the order k=1,2,…,m-1(2.4)ℒ{dαf(t)dtα}=sαf*(s)-∑k=0m-1f(k)(0+)sα-1-k,m-1<α<m,
whereas the Riemann-Liouville derivative for its Laplace transform rule requires the knowledge of the initial values of the fractional integral Im-αf(0+) and its derivatives of the order k=1,2,…,m-1(2.5)ℒ{DRLαf(t)}=sαf*(s)-∑k=0m-1DkIm-αf(0+)sm-1-k,m-1<α<m.
The reader interested in applications of integral transforms in fractional calculus is referred to [4–7, 24].


## 2.2. Hankel Transform

The Hankel transform is used to solve problem in cylindrical coordinates in the domain0≤r<∞ and is defined as
(2.6)ℋ{f(r)}=f*(η)=∫0∞f(r)Jν(ηr)rdr,ℋ-1{f*(η)}=f(r)=∫0∞f*(η)Jν(ηr)ηdη,
where Jν(r) is the Bessel function of the order ν.The following formula is fulfilled:(2.7)ℋ{d2f(r)dr2+1rdf(r)dr-ν2r2f(r)}=-η2f*(η).


## 2.3. Sin-Cos-Fourier Transform

In the case of boundary condition of the third kind with the prescribed boundary value of linear combination of a function and its normal derivative(2.8)z=0࢙:-dfdz+Hf=φ0,
the  following sin-cos-Fourier transform [25] is employed:
(2.9)ℱ{f(z)}=f*(ξ)=∫0∞K(z,ξ)f(z)dz,ℱ-1{f*(ξ)}=f(z)=2π∫0∞K(z,ξ)f*(ξ)dξ,
with the kernel
(2.10)K(z,ξ)=ξcos(zξ)+Hsin(zξ)ξ2+H2.Application of the sin-cos-Fourier transform to the second derivative of a function gives(2.11)ℱ{d2f(z)dz2}=-ξ2f*(ξ)-ξξ2+H2[df(z)dz-Hf(z)]|z=0.


## 3. Solution to the Problem under Mathematical Robin Boundary Condition

Consider the axisymmetric time-fractional heat conduction equation in cylindrical coordinates(3.1)∂αT∂tα=a(∂2T∂r2+1r∂T∂r+∂2T∂z2),0≤r<∞,0<z<∞,0<t<∞,0<α≤2,
with zero initial conditions
(3.2)t=0࢙:T=0,0<α≤2,t=0࢙:∂T∂t=0,1<α≤2,
and the mathematical Robin boundary condition
(3.3)z=0࢙:HT-∂T∂z=f(r,t).The zero conditions at infinity are also assumed(3.4)limr→∞T(r,z,t)=0,limz→∞T(r,z,t)=0.The solution to the initial-boundary-value problem (3.1)–(3.4) can be written as
(3.5)T=∫0t∫0∞f(ρ,τ)𝒢m(r,z,ρ,t-τ)ρdρdτ,
where 𝒢m(r,z,ρ,t) is the fundamental solution being the solution of the following problem:
(3.6)∂α𝒢m∂tα=a(∂2𝒢m∂r2+1r∂𝒢m∂r+∂2𝒢m∂z2),0≤r<∞,0<z<∞,0<t<∞,0<α≤2,t=0࢙:𝒢m=0,0<α≤2,t=0࢙:∂𝒢m∂t=0,1<α≤2,z=0࢙:H𝒢m-∂𝒢m∂z=1rδ(r-ρ)δ(t).
Here,  δ(r) is the Dirac delta function.The Laplace transform with respect to timet, the Hankel transform with respect to the spatial coordinate r, and the sin-cos-Fourier transform with respect to the spatial coordinate z result in
(3.7)𝒢m***=aξJ0(ρη)ξ2+H21sα+a(ξ2+η2),
where each of the integral transforms is denoted by the asterisk.Invertion of the integral transforms leads to(3.8)𝒢m(r,z,ρ,t)=2atα-1π∬0∞ξ2cos(zξ)+ξHsin(zξ)ξ2+H2×Eα,α[-a(ξ2+η2)tα]J0(rη)J0(ρη)ηdξdη.
Here, Eα,β is the Mittag-Leffler function in two parameters α and β defined by the series representation:
(3.9)Eα,β(z)=∑m=0∞zmΓ(αm+β),α>0,β>0,z∈C.
The essential role of the Mittag-Leffler function in fractional calculus results from the following formula for the inverse Laplace transform
(3.10)ℒ-1{sα-βsα+b}=tβ-1Eα,β(-btα).
Consider several particular cases of  solution (3.8).The caseH=0 corresponds to the mathematical Neumann boundary condition with the prescribed boundary value of the normal derivative of temperature, and the solution reads
(3.11)𝒢m(r,z,ρ,t)=2atα-1π∬0∞Eα,α[-a(ξ2+η2)tα]×cos(zξ)J0(rη)J0(ρη)ηdξdη.In the case of classical heat conduction (α=1), the Mittag-Leffler function:
(3.12)E1,1(-x2)=e-x2,
and taking into account the following integrals [26, 27]:
(3.13)∫0∞xe-ax2J0(bx)J0(cx)dx=12aexp(-b2+c24a)I0(bc2a),a>0,∫0∞e-ax2cos(cx)dx=π2aexp(-c24a),a>0,∫0∞1x2+b2e-ax2cos(cx)dx=π4beab2[e-bcerfc(ab-c2a)+ebcerfc(ab+c2a)],a>0,b>0,∫0∞xx2+b2e-ax2sin(cx)dx=π4eab2[e-bcerfc(ab-c2a)-ebcerfc(ab+c2a)],a>0,b>0,
where I0(x) is the modified Bessel function, erfcx is the complementary error function, we get
(3.14)𝒢m(r,z,ρ,t)=12texp(-r2+ρ2+z24at)I0(rρ2at)×{1πat-Hexp[(atH+z2at)2]erfc(atH+z2at)}.For the wave equation (α=2),
(3.15)E2,2(-x2)=sinxx,
and after evaluation of the necessary integrals appearing in solution (3.8) (under assumption ρ=0) [26, 27]
(3.16)∫0∞xx2+p2sin(bx2+p2)J0(cx)dx={1b2-c2cos(pb2-c2),0<c<b,0,0<b<c.∫0∞cos(px)cos(qx)x2+b2dx={π2be-bpcosh(bq),0<q<p,π2be-bqcosh(bp),0<p<q,b>0,∫0∞xcos(px)sin(qx)x2+b2dx={-π2e-bpsinh(bq),0<q<p,π2e-bqcosh(bp),0<p<q,b>0,
we obtain the solution
(3.17)𝒢m(r,z,0,t)={aat2-r2[δ(at2-r2-z)-He-H(at2-r2-z)]0<r<at,0<z<at2-r2,0at<r<∞,at2-r2<z<∞Of particular interest is also the caseα=1/2 for which
(3.18)E1/2,1/2(-x)=1π-xex2erfcx=2π∫0∞e-u2-2xuudu,x≥0,𝒢m(r,z,ρ,t)=12πt∫0∞exp(-u2-r2+ρ2+z28atu)I0(rρ4atu)×{12πaut1/4-Hexp[(2auHt1/4+z22aut1/4)2]×erfc(2auHt1/4+z22aut1/4)}du.


## 4.  Solution to the Problem under Physical Robin Boundary Condition

Consider the following  axisymmetric time-fractional heat conduction equation:(4.1)∂αT∂tα=a(∂2T∂r2+1r∂T∂r+∂2T∂z2),0≤r<∞,0<z<∞,0<t<∞,0<α≤2,
under zero initial conditions
(4.2)t=0࢙:T=0,0<α≤2,t=0࢙:∂T∂t=0,1<α≤2,
and the physical Robin boundary condition
(4.3)z=0࢙:HT-DRL1-α∂T∂z=f(r,t),0<α≤1,z=0࢙:HT-Iα-1∂T∂z=f(r,t),1<α≤2,
where H=h/k.The solution to the initial-boundary-value problem (4.1)–(4.3) has the form
(4.4)T=∫0t∫0∞f(ρ,τ)𝒢p(r,z,ρ,t-τ)ρdρdτ,
where the fundamental solution 𝒢p(r,z,ρ,t) fulfills the following equation
(4.5)∂α𝒢p∂tα=a(∂2𝒢p∂r2+1r∂𝒢p∂r+∂2𝒢p∂z2),0≤r<∞,0<z<∞,0<t<∞,0<α≤2,
under the  following conditions:
(4.6)t=0࢙:𝒢p=0,0<α≤2,t=0࢙:∂𝒢p∂t=0,1<α≤2,z=0:H𝒢p-DRL1-α∂𝒢p∂z=1rδ(r-ρ)δ(t),0<α≤1,z=0:H𝒢p-Iα-1∂𝒢p∂z=1rδ(r-ρ)δ(t),1<α≤2.The Laplace transform with respect to timet leads to the following equation:
(4.7)sα𝒢p*=a(∂2𝒢p*∂r2+1r∂𝒢p*∂r+∂2𝒢p*∂z2),0≤r<∞,0<z<∞,0<α≤2,
and the  following  boundary condition:
(4.8)z=0:sα-1H𝒢p*-∂𝒢p*∂z=1rδ(r-ρ)sα-1.In this case, the kernel of the sin-cos-Fourier transform with respect to the spatial coordinatez depends on the Laplace transform variable s(4.9)K(z,ξ,s)=ξcos(zξ)+sα-1Hsin(zξ)ξ2+(sα-1H)2,
and in the transform domain we obtain
(4.10)𝒢p***=aξJ0(ρη)ξ2+(sα-1H)2sα-1sα+a(ξ2+η2).Inversion of the Laplace transform in (4.10) depends on the value of α. For 0<α≤1, we have
(4.11)𝒢p(r,z,ρ,t)=2aπ∬0∞∫0tηJ0(rη)J0(ρη)τα-1Eα,α[-a(ξ2+η2)τα]×{(t-τ)-αE2-2α,1-α[-H2ξ2(t-τ)2-2α]cos(zξ)+Hξ(t-τ)1-2αE2-2α,2-2α[-H2ξ2(t-τ)2-2α]sin(zξ)}dτdξdη,
whereas for  1<α≤2,  we get
(4.12)𝒢p(r,z,ρ,t)=2aπ∬0∞∫0tJ0(rη)J0(ρη)Eα[-a(ξ2+η2)τα]×{ξ2H2(t-τ)2α-3E2α-2,2α-2[-ξ2H2(t-τ)2α-2]cos(zξ)+ξH(t-τ)α-2E2α-2,α-1[-H2ξ2(t-τ)2-2α]sin(zξ)}dτdξdη.
The fundamental solution under physical Neumann boundary condition is obtained forH=0 as
(4.13)𝒢p(r,z,ρ,t)=2aπ∬0∞Eα[-a(ξ2+η2)tα]J0(rη)J0(ρη)ηcos(zξ)dξdη.Formulae (4.11) and (4.12) simplify for α=1/2 and α=3/2, respectively, taking into account (3.12) and that
(4.14)E1,1/2(-x2)=1π[1-2xD(x)],
where
(4.15)D(x)=e-x2∫0xeu2du
is the Dawson integral.


## 5. Conclusion

We have derived the analytical solutions to time-fractional heat conduction equation in a half-space under mathematical and physical Robin boundary conditions. The integral transform technique has been used. It should be emphasized that in the case of physical Robin boundary condition, the order of integral transforms is important as the kernel of the sin-cos-Fourier transform depends on the Laplace transform variable. The limiting caseH=0 corresponds to solutions of problems under mathematical and physical Neumann boundary conditions with the prescribed boundary value of the normal derivative and with the prescribed boundary value of the heat flux, respectively. The difference between mathematical and physical boundary conditions (as well as the difference between the solutions) disappears in the case of standard heat conduction equation (α=1).



---

*Source: 154085-2012-08-29.xml*