Temperature prediction by a fractional heat conduction model for the bi-layered spherical tissue in the hyperthermia experiment

International Journal of Thermal Sciences 145 (2019) 105990

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Temperature prediction by a fractional heat conduction model for the bilayered spherical tissue in the hyperthermia experiment

T

Bo Yua,*, Xiaoyun Jiangb a b

School of Mathematics and Statistics, Shandong University, Weihai, Weihai 264209, China School of Mathematics, Shandong University, Jinan 250100, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Fractional heat conduction Bi-layered spherical tissue Finite difference Conjugate gradient method Fractional sensitivity equation

In this paper, a novel time-fractional heat conduction model is presented for the bi-layered spherical tissue in the hyperthermia experiment, and an efficient numerical technique to estimate the unknown fractional parameter in the proposed time-fractional heat conduction model is investigated. Firstly, for the solution of the proposed timefractional heat conduction model, the technique of variable transformation and the implicit finite difference method are employed. Secondly, based on the measured experimental data, an efficient numerical technique to estimate the unknown order of the fractional derivative is investigated. The fractional sensitivity equation is first obtained by means of the digamma function, and the conjugate gradient method based on the fractional derivative is constructed to obtain the optimal estimation of the unknown order of the fractional derivative in the proposed time-fractional heat conduction model. Lastly, compared with the hyperthermia experimental data, it is observed that the calculated temperature increase values agree well with the measured temperature increase values in the experiment, which implies that the proposed time-fractional heat conduction model is effective in modeling the heat transfer in the hyperthermia experiment, and the proposed numerical technique to estimate the unknown fractional parameter is efficient.

1. Introduction In the past few years, fractional calculus has received worldwide attention due to its broad range of applications in various areas of science and engineering. Many literatures have shown that the powerlaw behavior is a hallmark of many biological phenomena observed at different scales and at various levels of organization [1–4] and that a rheological behavior that conforms to the power law can be described by using methods of fractional calculus [5–7]. Fractional calculus provides an excellent instrument for the description of memory and hereditary properties of various materials and processes [8]. For example, fractional calculus provides novel mathematical tools for modeling physical and biological processes [9]. Hence, the concept of fractional order modeling in biological systems has received significant interest in the research community [10–13]. Magin et al. [14] described the formulation of the bioheat transfer in one dimension in terms of the fractional order differentiation with respect to time. Their study demonstrated that fractional calculus could provide a unified approach to examine the periodic heat transfer in peripheral tissue regions. Recently, fractional heat conduction equation in composite medium has received the scholars' interest. Ilic et al. [15] solved a one-

*

dimensional space fractional diffusion equation in a composite medium consisting of two layers in contact both analytically and numerically. Povstenko [16,17] studied fractional heat conduction in a composite medium consisting of two semi-infinite regions being in perfect thermal contact and in a composite medium with a spherical inclusion. Jiang and Chen [18] set up a time fractional thermal diffusion equation in one-dimensional composite medium, and gave analytical and numerical methods to solve the equation of fractional thermal diffusion in Mlayers composite medium. Zhuang et al. [19] studied an inverse problem of parameter estimation for time-fractional heat conduction in a composite medium using carbon-carbon experimental data by the Levenberg-Marquardt iterative method. The research of the inverse problems involving fractional differential or integral equations is a relatively novel topic. Battaglia [20] solved an inverse heat conduction problem using a non-integer identified model. Murio [21] investigated the numerical solution of the Caputo time fractional inverse heat conduction problem on a finite slab in the presence of measured (noisy) data. Qian [22] considered the determination of the boundary temperature from one measured transient data temperature at some interior point of a one-dimensional semi-infinite conductor. Murio [23] introduced a simple algorithm based on

Corresponding author. E-mail addresses: [email protected] (B. Yu), [email protected] (X. Jiang).

https://doi.org/10.1016/j.ijthermalsci.2019.105990 Received 23 January 2017; Received in revised form 31 December 2018; Accepted 19 June 2019 Available online 02 July 2019 1290-0729/ © 2019 Elsevier Masson SAS. All rights reserved.

International Journal of Thermal Sciences 145 (2019) 105990

B. Yu and X. Jiang

demonstrates that the classical bio-heat transfer equation may not completely describe the thermal behavior captured in the experiment [30]. Based on the theory of the fractional calculus [8], we here present a novel time-fractional heat conduction model for the bi-layered spherical tissue which combines the heat transfer in the tumor (0 ≤ r ≤ R ) and healthy tissues (R < r ≤ a ) with constant physiological parameters in the following form:

space marching mollification techniques and Gruünwald-Letnikov fractional derivatives for the numerical solution of the recovery of the boundary temperature and heat flux functions from one measured transient temperature data at some interior point of a one-dimensional semi-infinite conductor. Wei et al. [24] developed a coupled method for inverse source problem of spatial fractional anomalous diffusion equations. Ghazizadeh et al. [25] studied an inverse heat conduction problem in processed meat for simultaneous estimation of relaxation parameter and order of fractionality in fractional single-phase-lag heat equation. In our previous studies, we studied an inverse problem to estimate the unknown order of the Riemann-Liouville fractional derivative for a fractional Stokes' first problem for a heated generalized second grade fluid [26], studied the numerical identification of the two fractional derivatives in the two-dimensional fractional Cable equation [27], and studied the numerical algorithms to estimate relaxation parameters and Caputo fractional derivative for a fractional thermal wave model in spherical composite medium [28]. In this paper, we first present a novel time-fractional heat conduction model for the bi-layered spherical tissue in the hyperthermia experiment, and then propose an efficient numerical technique to estimate the unknown fractional parameter in the proposed time-fractional heat conduction model. From the methodological point of view, this kind of technique to deal with the fractional inverse heat conduction problem in a spherical composite medium is very novel and never be observed in the literature. The structure of this paper is organized as follows. In Section 2, we present the time-fractional heat conduction model for the bi-layered spherical tissue in the hyperthermia experiment. With change of variable, the original problem in the spherical coordinate system can be transformed into a new tractable problem, implicit finite difference method is proposed to solve the new problem in Section 3. In Section 4, for the parameter estimation of the proposed time-fractional heat conduction model, fractional sensitivity equation is first obtained by means of the digamma function, and then conjugate gradient method based on fractional derivative is constructed to obtain the optimal estimation of the unknown fractional derivative. Compared with the hyperthermia experimental data, the optimal estimation of the fractional parameter is obtained, and some detailed discussions are given in Section 5. Finally, we give some conclusions.

ρ1 c1

∂αT ρ2 c2 α2 ∂t

ρ2 c2

k2 ∂ ⎛ 2 ∂T2 ⎞ r , R < r ≤ a, r 2 ∂r ⎝ ∂r ⎠

(3)

(4)

t

(5)

Fractional derivative α is a convolution of the temperature derivative and kernel function. It implies thermal diffusion has time memory and the temperature change is related to history process. The boundary conditions are obtained as follows by the assumptions that the temperature is continuous at the surface of the sphere and the heat flux is continuous across the interface:

T1 (0, t ) is finite,

(6)

T2 (a, t ) = T0,

(7)

T1 (R, t ) = T2 (R, t ),

(8)

k1

∂T1 (R, t ) ∂T (R, t ) = k2 2 . ∂r ∂r

(9)

The initial conditions are given by

Ti (r , 0) = T0, i = 1,2.

(10)

3. Numerical treatment of the fractional heat conduction model To solve the fractional heat conduction model (3)–(4), we first introduce a new variable θ defined by

θ = r (T − T0),

(11)

to transform the original problem in the spherical coordinate system into a new tractable problem. Under the technique of variable transformation, we have

∂αθ ∂ α T ∂ 2θ 1 ∂ ⎛ 2 ∂T ⎞ ∂θ θ ∂T =r α, = r , = +r . ∂t α ∂t ∂r 2 r ∂r ⎝ ∂r ⎠ ∂r r ∂r

(12)

Hence, the original problem (3)–(4) can be rewritten in terms of θ, respectively, as

ρ1 c1

∂T2 ∂T k ∂ = 22 ⎛r 2 2 ⎞, r ≥ R. ∂t r ∂r ⎝ ∂r ⎠

0 < α ≤ 1, t > 0,

f (n) (τ ) 1 ⎧ ⎪ Γ(n − α ) ∫ (t − τ )α − n + 1 dτ , n − 1 < α < n, ∂αf (t ) = 0 ∂t α ⎨ ⎪ f (n) (t ), α = n. ⎩

To mathematically predict the tissue temperature distributions in the hyperthermia experiment, Andrä et al. [29] modeled a breast carcinomas surrounded by heathy tissues as a small solid sphere without considering the effects of the blood perfusion and metabolism. In their work, they restricted the considerations to small tumors in relatively homogeneous tissue with low blood perfusion such as muscle tissue or fat, so that the perfusion may be taken into account by modification of the heat conduction equation. Assume the heat source of constant power density P is concentrated within a small sphere of radius R surrounded by a medium of homogeneous heat conductivity, heat transfers in the radius direction symmetrically, and the temperature distribution in the tumor and healthy tissues depends only on the distance r from the centre of the sphere and time t [29]. The mathematical model is proposed in the following form [29]:

0 ≤ r ≤ R, t > 0,

=

0 ≤ r ≤ R,

where ρi , ci, ki, Ti (i = 1,2) denote the mass density, specific heat capacity, heat conductivity and temperature for the interior and exterior of the sphere, respectively, a denotes the boundary of tissue, T0 is regarded α as the arterial temperature, P is the heat power density, ∂ f α(t ) is the ∂t Caputo time-fractional derivative defined by Ref. [8].

2. Fractional heat conduction model

∂T ∂T k ∂ ρ1 c1 1 = 12 ⎛r 2 1 ⎞ + P , ∂t r ∂r ⎝ ∂r ⎠

∂αT1 ∂T k ∂ = 12 ⎛r 2 1 ⎞ + P , ∂t α r ∂r ⎝ ∂r ⎠

(1)

∂αθ1 ∂ 2θ = k1 21 + Pr , ∂t α ∂r

∂αθ ρ2 c2 α2 ∂t

(2)

=

∂ 2θ k2 22 , ∂r

0 ≤ r ≤ R,

R < r ≤ a.

0 < α ≤ 1, t > 0,

(13) (14)

The initial and boundary conditions are also rewritten in terms of θ as follows

From Ref. [29], we can find out that there is an apparent difference between the calculated temperature increase values from the classical bio-heat transfer equation and the measured temperature increase values in the experiment indicated with the symbols. This phenomenon 2

θ1 (0, t ) = 0,

(15)

θ2 (a, t ) = 0,

(16)

International Journal of Thermal Sciences 145 (2019) 105990

B. Yu and X. Jiang

θ1 (R, t ) = θ2 (R, t ), k1 ⎛ ⎝

(17)

∂θ1 (R, t ) θ (R, t ) ⎞ ∂θ (R, t ) θ (R, t ) ⎞ − 1 = k2 ⎛ 2 − 2 , ∂r R ⎠ ∂r R ⎠ ⎝

θi (r , 0) = 0, i = 1,2.

available. The solution of the inverse problem is obtained through the minimization of the least squares norm between the measured temperature increase vector ΔT in the experiment and the calculated temperature increase vector ΔTc from the direct problem in the following form

(18) (19)

S (α ) = (ΔT − ΔTc)T (ΔT − ΔTc),

It is worthy to mention that the temperature increase ΔT is

ΔT = T − T0 =

(29)

where

θ , r

(20)

ΔT = [ΔT (t1,1), ΔT (t1,2), …, ΔT (t1, n ), ΔT (t2,1), ΔT (t2,2), …, ΔT (t2, n)]T ,

θ

(30)

and the value of ΔT = r at r = 0 is indeterminate and must be replaced by its limit as r → 0 . In order to solve the above new problem (13)–(14) numerically, we first introduce a uniform grid of mesh points: 0 = r1,0 < r1,1 < … < r1, M = R = r2,0 < r2,1 < … < r2, M = a , 0 = t0 < t1 < … < tN = T , where Δri = ri, j − ri, j − 1, i = 1,2, j = 1,2, …, M , and Δt = tk − tk − 1, k = 1,2, …, N . For 0 < α < 1, the Caputo time-fractional derivative can be discretized as

∂αθ (r , t ) t = tk + 1 (Δt )−α |r = ri, j = ∂t α Γ(2 − α ) 1)1 − α

and

ΔTc = [ΔTc (t1,1, α ), ΔTc (t1,2, α ), …, ΔTc (t1, n, α ), ΔTc (t2,1, α ), ΔTc , (t2,2, α ), …, ΔTc (t2, n, α )]T

ti, j denotes the measured time nodes, n is the total number of the measured data from the experiment. 4.1. Conjugate gradient method

k

∑ bl (θik,j−l+1 − θik,j−l), l=0

Conjugate gradient method is a straightforward and powerful iterative technique for solving linear and nonlinear inverse problems of parameter estimation [31]. In this paper, for estimating the unknown parameter α, conjugate gradient method is introduced to minimize the above least squares norm. The iterative algorithm for minimization of the above norm S (α ) is

(21)

l1 − α ,

where bl = (l + − l = 0,1,2, ⋯. In order to guarantee the stability and accuracy of the numerical method, we here adopt the implicit finite difference schemes as follows

ρ1 c1

1 k+1 k+1 θ1,k + (Δt )−α k j + 1 − 2θ1, j + θ1, j − 1 + Pr1, j, bl (θ1,k j− l + 1 − θ1,k j− l ) = k1 ∑ Γ(2 − α ) l = 0 (Δr1)2

α k + 1 = α k + β kdk ,

(22)

θ2,k +j +11 − 2θ2,k +j 1 + θ2,k +j −11 (Δt )−α k ρ2 c2 bl (θ2,k −j l + 1 − θ2,k −j l ) = k2 . ∑ Γ(2 − α ) l = 0 (Δr2)2

(31)

(32)

dk

βk

is the search step size, is the direction of descent and the where superscript k is the number of iterations. The direction of descent is given by

(23)

dk = ∇S (α k ) + γ kdk − 1.

The boundary and initial conditions can also be discretized as

(33)

∇S (α k )

k+1 θ1,0 = 0,

(24)

θ2,k +M1 = 0,

(25)

is the component of the gradient direction evaluated Here, at iteration k, the expression of the gradient direction is obtained by differentiating equation (29) with respect to the unknown parameter α, that is,

1 k+1 θ1,k + M = θ2,0 ,

(26)

∇S (α k ) = −2(Jk)T (ΔT − ΔTc (α k )),

(34)

Jk

1 k+1 k+1 k+1 k+1 1 θ1,k + θ1,k + M − θ1, M − 1 M ⎞ ⎛ θ2,1 − θ2,0 − θ2,0 ⎞, k1 ⎛⎜ − ⎟ = k2 ⎜ ⎟ Δr1 R ⎠ Δr2 R ⎠ ⎝ ⎝

(27)

θi0, j = 0, i = 1,2.

(28)

is the fractional sensitivity matrix (Jacobian matrix in the where fractional sense). The elements of the fractional sensitivity matrix are the first derivative of the estimated temperature increase with respect to the unknown parameter α. The conjugation coefficient γ k is given as

γ 0 = 0,

By means of Thomas algorithm, the numerical solution can be easily obtained by solving (22)–(28).

γk =

[∇S (α k )][∇S (α k ) − ∇S (α k − 1)] , k = 1,2, …. [∇S (α k − 1)]2

The search step size 4. Parameter estimation

βk = The above fractional heat conduction model contains parameters α, ρ1, c1, ρ2 , c2 , k1, k2 , P, R and a, all of which are assumed known in order to obtain the solution of the proposed fractional heat conduction model. In the inverse analysis, the measured temperature increase values in the experiment are already known, the fractional derivative α is regarded as being unknown while other parameters ρ1, c1, ρ2 , c2 , k1, k2 , P, R and a are already known in the hyperthermia experiment. Fractional derivative α is a very important parameter in determining the temperature distribution. It is a convolution of the temperature derivative and kernel function which implies that thermal diffusion has time memory and the temperature change is related to history process. The computed temperature distribution is forced to match the experimental temperature distribution, hence, the optimal estimation of the fractional derivative α is obtained. In Ref. [29], the measured temperature increase ΔT as a function of time for different reduced distances r / R from the composite centre are

[Jkdk ]T [ΔT

− ΔTc [Jkdk ]T [Jkdk ]

βk

(35)

is given by

(α k )]

.

(36)

The iteration begins with an appropriate initial guess for the unknown parameter α and is stopped when the following stopping criterion is satisfied

S (α k + 1) < ε,

(37)

where the value of the tolerance ε is chosen so that sufficiently stable solutions are obtained. 4.2. Computation of the fractional sensitivity matrix The fractional sensitivity equations are obtained through differentiating (13)–(14) with respect to the unknown parameter α, that is,

ρi ci 3

∂ ⎛ ∂αθi (r , t ) ⎞ ∂ 2J = ki 2i , i = 1,2, ∂α ⎝ ∂t α ⎠ ∂r

(38)

International Journal of Thermal Sciences 145 (2019) 105990

B. Yu and X. Jiang

where Ji (r , t ) =

∂θi (r , t ) , ∂α

matrix can be easily obtained.

i = 1,2 .

Concerning the discrete form of the term ∂αθi (r , t ) ∂t α

∂ ∂α

(

=

∂ ∂α

=

(Δt )−α Γ(2 − α )

(

(

∂αθi (r , t ) ∂t α

), we have

4.3. Computational algorithm

)

(Δt )−α Γ(2 − α )

(

∂ ∂α

∑l = 0k bl (θik, j− l + 1 − θik, j− l )

Assume that the measured temperature increase vector ΔT is already known, the initial guess α 0 is chosen for the unknown fractional derivative α, set k = 0 :

)

∑l = 0k ⎡bl (Jik, j− l + 1 − Jik, j− l ) + dl (θik, j− l + 1 − θik, j− l ) ⎢ ⎣

1

+ − Γ(2 − α )

∂Γ(2 − α ) ∂α

Step 1. Solve the fractional heat conduction model (22)–(28) with the available estimate α k to gain the computed temperature increase vector

)

− ln Δt bl (θik, j− l + 1 − θik, j− l ) ⎤ ⎥ ⎦

ΔTck = [ΔTc (t1,1, α k ), ΔTc (t1,2, α k ), …, ΔTc (t1, n, α k ), ΔTc ; (t2,1, α k ), ΔTc (t2,2, α k ), …, ΔTc (t2, n, α k )]T

(Δt )−α Γ(2 − α )

∑l = 0k ⎡bl (Jik, j− l + 1 − Jik, j− l ) + dl (θik, j− l + 1 − θik, j− l ) ⎢ ⎣ + (ψ (2 − α ) − ln Δt ) bl (θik, j− l + 1 − θik, j− l )] =

∂αJi (r , t ) ∂t α

+

Check the stopping criteria. If satisfied, stop the Step 2. Compute iteration; otherwise, continue. Step 3. Compute the fractional sensitivity matrix Jk by (47)–(52);

(Δt )−α Γ(2 − α )

∑l = 0k ⎡dl (θik, j− l + 1 − θik, j− l ) ⎢ ⎣ k−l+1 + (ψ (2 − α ) − ln Δt ) bl (θi, j − θik, j− l )], =

Step 4. Compute the gradient direction ∇S (α k ) and the conjugation coefficient γ k ;

(39)

Step 5. Compute the direction of descent dk ;

where (40)

Step 6. Compute the search step size β k ;

is the digamma function [26]. and ψ (x ) = Hence, the fractional sensitivity equations have been transformed into the following forms:

Step 7. Compute the new estimate α k + 1;

d 0 = 0,

dl = l1 − α ln l − (l + 1)1 − α ln(1 + l),

l = 1,2, ⋯,

ln Γ(x ) dx

ρi ci

Step 8. Replace k by k + 1 and return to Step 1. The graphical representation of the parameter estimation process is depicted in Fig. 1.

∂αJi (r , t ) ∂t α

(Δt )−α − ρi ci Γ(2 − α ) ∑l = 0k ⎡dl (θik, j− l + 1 − θik, j− l ) ⎢ ⎣ + (ψ (2 − α ) − ln Δt ) bl (θik, j− l + 1 − θik, j− l )], i = 1,2.

= ki

∂2Ji ∂r 2

5. Results and discussions (41)

Andrä et al. [29] simulated the thermal behavior of localized magnetic hyperthermia applied to a breast carcinoma. They considered

The boundary and initial conditions are

J1 (0, t ) = 0,

(42)

J2 (a, t ) = 0,

(43)

J1 (R, t ) = J2 (R, t ),

(44)

k1 ⎛ ⎝

(53)

S (α k ) .

J (R, t ) ⎞ J (R, t ) ⎞ ∂J1 (R, t ) ∂J (R, t ) − 1 = k2 ⎛ 2 − 2 , ∂r R ⎠ ∂ r R ⎠ ⎝

Ji (r , 0) = 0, i = 1,2.

(45) (46)

We here employ the implicit finite difference scheme similarly as in Section 3 (Δt )−α

ρi ci Γ(2 − α ) ∑l = 0k bl (Jik, j− l + 1 − Jik, j− l ) = ki

Jik, j++11 − 2Jik, j+ 1 + Jik, j+−11 (Δri)2

(Δt )−α − ρi ci Γ(2 − α ) ∑l = 0k ⎡dl (θik, j− l + 1 − θik, j− l ) ⎢ ⎣ + (ψ (2 − α ) − ln Δt ) bl (θik, j− l + 1 − θik, j− l )], i = 1,2.

(47)

The boundary and initial conditions can be discretized as k+1 J1,0

= 0,

(48)

J2,k +M1 = 0,

(49)

1 k+1 J1,k + M = J2,0 ,

(50)

k+1 k+1 k+1 k+1 1 1 J1,k + J1,k + M − J1, M − 1 M ⎞ ⎛ J2,1 − J2,0 − J2,0 ⎞, k1 ⎛⎜ − ⎟ = k2 ⎜ ⎟ Δ r R Δ r R ⎠ 1 2 ⎝ ⎠ ⎝

(51)

Ji0, j = 0, i = 1,2.

(52) Fig. 1. Graphical representation of the parameter estimation process.

Solving (47)–(52) with Thomas algorithm, the fractional sensitivity 4

International Journal of Thermal Sciences 145 (2019) 105990

B. Yu and X. Jiang

a small sphere tumor of radius R = 0.00315m with a constant power density P = 6150000W / m3 embedded in extended muscle tissue. The thermal parameters were taken to be ρ1 = 1660000g / m3 , c2 = 3.72Ws /(gK ) , c1 = 2.54Ws /(gK ) , k1 = ρ2 = 1000000g / m3 , 0.778W /(Km) , k2 = 0.642W /(Km) . From the experimental study, the distribution of the temperature increase as a function of exposure time was measured with thermocouples. Measured values of temperature increase at various reduced distances r / R = 0.65, 1.10, 1.39, 1.68, 1.98 were presented with symbols.

5.1. Parameter estimation Experimental data are extracted from the symbols at various reduced distances r / R = 0.65, 1.10 , where 11 measurement values are at r / R = 0.65 and another 11 measurement values are at r / R = 1.10 , the error of ΔT and t are within the extent of the used symbols. The thermal parameters ρ1, ρ2 , c1, c2 , k1, k2 are taken as above. For the parameter a, we here choose a = 11R to ensure the experimental data points just fall on the numerical grid point. Chosen a initial guess α 0 , the estimation of the order of the fractional derivative is implemented by the proposed numerical method. The computational results are listed in Table 1, where α inv is the optimal estimation of the order of the fractional derivative. From Table 1, we can obviously find out that, for different initial guess α 0 = 0.9, 0.8, 0.7 , the estimated values α inv are almost the same when the iteration number k ≥ 3, which demonstrates that different initial guess α 0 has little impact on the final estimation of the order of the fractional derivative. Fig. 2 depicts the estimation process for different initial guess α 0 = 0.9, 0.8, 0.7 . All these facts demonstrate that the proposed numerical technique is effective and efficient for the estimation of the order of the fractional derivative in the proposed time-fractional heat conduction model for the bi-layered spherical tissue in the hyperthermia experiment.

Fig. 2. Estimation results of the fractional derivative for different initial guess α0.

5.2. Simulation of the time-fractional heat conduction model To qualify the accuracy of our present results, the comparison of the present results with those given in Ref. [29] is made. In Fig. 3, we first depict the comparison between the measured temperature increase values in the experiment [29], the estimated temperature increase values when α is chosen as the estimated value α inv = 0.98211 at the reduced distances r / R = 0.65, 1.10 , respectively. Then we also depict the comparison of the results at various reduced distances r / R = 1.39, 1.68, 1.98 in the same Fig. 3. From Fig. 3, we can obviously find out that the estimated temperature increase values agree well with the measured temperature increase values in the experiment when α is chosen as the estimated value α inv = 0.98211. This fact demonstrates that the proposed timefractional heat conduction model is effective and accurate in modeling the heat transfer in the hyperthermia experiment. In addition, the present result seems a little better than the result in the classical model [29].

Fig. 3. Comparison between the experimental data and the estimated temperature increase values at various reduced distances.

5.3. Influence of the fractional derivative for the time-fractional heat conduction model Fractional derivative α is a significant parameter in predicting the temperature distribution in the hyperthermia experiment. It is a convolution of the temperature derivative and kernel function which implies that thermal diffusion has time memory and the temperature change is related to history process. In this subsection, we will discuss the influence of the fractional derivative α for the temperature increase ΔT . Fig. 4 demonstrates the comparison between the measured temperature increase values in the experiment [29], the calculated temperature increase values when α is chosen as the estimated value α = 0.98211 and another value α = 0.9 at the reduced distances r / R = 0.65, 1.10 . From Fig. 4, we can obviously find out that the temperature increase ΔT are higher when α is larger, and lower when α is smaller. This fact demonstrates that the heat conduction process in this experiment is an anomalous diffusion process.

Table 1 Estimation results of the fractional derivative. Initial guess α 0

Iteration number k

Estimated value α inv

0.9 0.8 0.7

k≥3 k≥3 k≥3

0.98211 0.98211 0.98211

5.4. Sensitivity analysis of the fractional derivative α The sensitivity coefficients Jα are measurements of sensitivity of the 5

International Journal of Thermal Sciences 145 (2019) 105990

B. Yu and X. Jiang

spherical tissue in the hyperthermia experiment. With change of variable, we transform the original problem in the spherical coordinate system into a new tractable problem. Implicit finite difference method is employed to solve the problem. (2). For the inverse problem, we investigate an efficient numerical method to estimate the unknown order of the fractional derivative. Fractional sensitivity equation is obtained by means of the digamma function, the conjugate gradient method based on the fractional derivative is constructed to obtain the optimal estimation of the unknown order of the fractional derivative. Compared with the experimental data, it is observed that the estimated temperature increase values agree well with the measured temperature increase values in the experiment. (3). The proposed time-fractional heat conduction model is effective and accurate in modeling the heat transfer in the hyperthermia experiment. Fractional calculus provides a novel mathematical tool to model the heat transfer phenomenon. (4). The proposed numerical inversion method is effective and efficient for the estimation of fractional derivative in the proposed timefractional heat conduction model in the bi-layered spherical tissue in the hyperthermia experiment. From the methodological point of view, this kind of numerical technique to deal with the fractional inverse heat conduction problem in a spherical composite medium is very novel, and it can be applied to other similar parameter estimation problems involving fractional differentiating or integral equations.

Fig. 4. Influence of the fractional derivative.

Acknowledgements The authors thank the editors and the referees for their valuable comments and helpful suggestions. The project was supported by Shandong Provincial Natural Science Foundation, China (Grant number ZR2017MA030), China Postdoctoral Science Foundation (Grant numbers 2018T110679 and 2016M602127) and National Natural Science Foundation of China (Grant numbers 11771254 and 11672163). References [1] B. Suki, A.L. Barabasi, K.R. Lutchen, Lung tissue viscoelasticity: a mathematical framework and its molecular basis, J. Appl. Phys. 76 (1994) 2749–2759. [2] B. Suki, A.L. Barabasi, Z. Hantos, F. Peták, H.E. Stanley, Avalanches and power-law behaviour in lung inflation, Nature 368 (1994) 615–618. [3] R.G. Larson, The Structure and Rheology of Complex Fluids, New York Oxford university press, 1999. [4] B. Fabry, G.N. Maksym, J.P. Butler, M. Glogauer, D. Navajas, J.J. Fredberg, Scaling the microrheology of living cells, Phys. Rev. Lett. 87 (2001) 148102. [5] R.L. Magin, Fractional Calculus in Bioengineering, Begell House Publishers, 2006. [6] W.C. Tan, C.Q. Fu, C.J. Fu, W.J. Xie, H.P. Cheng, An anomalous subdiffusion model for calcium spark in cardiac myocytes, Appl. Phys. Lett. 91 (2007) 183901. [7] R.L. Magin, Fractional calculus models of complex dynamics in biological tissues, Comput. Math. Appl. 59 (2010) 1586–1593. [8] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [9] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Physiol. Anthropol. 37 (2004) R161–R208. [10] M. Kohandel, S. Sivaloganathan, G. Tenti, K. Darvish, Frequency dependence of complex moduli of brain tissue using a fractional Zener model, Phys. Med. Biol. 50 (2005) 2799–2805. [11] D. Craiem, F.J. Rojo, J.M. Atienza, R.L. Armentano, G.V. Guinea, Fractional-order viscoelasticity applied to describe uniaxial stress relaxation of human arteries, Phys. Med. Biol. 53 (2008) 4543–4554. [12] J. Singh, P.K. Gupta, K.N. Rai, Solution of fractional bioheat equations by finite difference method and HPM, Math. Comput. Model. 54 (2011) 2316–2325. [13] X.Y. Jiang, H.T. Qi, Thermal wave model of bioheat transfer with modified Riemann-Liouville fractional derivative, J. Physiol. Anthropol. 45 (2012) 485101. [14] R.L. Magin, S.C. Boregowda, C. Deodhar, Modeling of pulsating peripheral bio-heat transfer using fractional calculus and constructal theory, Int. J. Des. Nat. Ecodyn. 1 (2007) 18–33. [15] M. Ilic, I. Turner, F. Liu, V. Anh, Analytical and numerical solutions of a one-dimensional fractional-in-space diffusion equation in a composite medium, Appl. Math. Comput. 216 (2010) 2248–2262. [16] Y.Z. Povstenko, Fractional heat conduction in infinite one-dimensional composite medium, J. Therm. Stress. 36 (2013) 351–363.

Fig. 5. Sensitivity coefficients at the reduced distances r / R = 0.65, 1.10 .

estimated temperature increase ΔT with respect to the changes in the parameter α at the reduced distances r / R = 0.65, 1.10 . Here, we introduce the relative sensitivity coefficients [25,31].

Jα =

α ∂ (ΔT) . T0 ∂α

(54)

Fig. 5 depicts the sensitivity coefficients Jα at the reduced distances r / R = 0.65, 1.10 , respectively. From Fig. 5, we can obviously observe that the sensitivity coefficients are not equal to zero for t > 0 , it is possible to obtain the correct estimation of the unknown fractional derivative α. 6. Conclusions In this paper, we present a novel time-fractional heat conduction model for the bi-layered spherical tissue in the hyperthermia experiment, and propose an efficient numerical technique to estimate the unknown fractional parameter in the proposed time-fractional heat conduction model. (1). Based on the theory of the fractional calculus, we first present a novel time-fractional heat conduction model for the bi-layered

6

International Journal of Thermal Sciences 145 (2019) 105990

B. Yu and X. Jiang

[25] H.R. Ghazizadeh, A. Azimi, M. Maerefat, An inverse problem to estimate relaxation parameter and order of fractionality in fractional single-phase-lag heat equation, Int. J. Heat Mass Transf. 55 (2012) 2095–2101. [26] B. Yu, X.Y. Jiang, H.T. Qi, An inverse problem to estimate an unknown order of a Riemann-Liouville fractional derivative for a fractional Stokes' first problem for a heated generalized second grade fluid, Acta Mech. Sin. 31 (2015) 153–161. [27] B. Yu, X.Y. Jiang, Numerical identification of the fractional derivatives in the twodimensional fractional Cable equation, J. Sci. Comput. 68 (1) (2016) 252–272. [28] B. Yu, X.Y. Jiang, C. Wang, Numerical algorithms to estimate relaxation parameters and Caputo fractional derivative for a fractional thermal wave model in spherical composite medium, Appl. Math. Comput. 274 (2016) 106–118. [29] W. Andrä, C.G. d'Ambly, R. Hergt, I. Hergt, W.A. Kaiser, Temperature distribution as function of time around a small spherical heat source of local magnetic hyperthermia, J. Magn. Magn. Mater. 194 (1999) 197–203. [30] K.C. Liu, C.T. Lin, Solution of an inverse heat conduction problem in a bi-layered spherical tissue, Numer. Heat Tran. 58 (2010) 802–818. [31] M.N. Özisik, Inverse Heat Transfer: Fundamentals and Applications, Taylor Francis, NewYork, 2000.

[17] Y.Z. Povstenko, Fractional heat conduction in an infinite medium with a spherical inclusion, Entropy 15 (2013) 4122–4133. [18] X.Y. Jiang, S. Chen, Analytical and numerical solutions of time fractional anomalous thermal diffusion equation in composite medium, ZAMM 95 (2015) 156–164. [19] Q. Zhuang, B. Yu, X.Y. Jiang, An inverse problem of parameter estimation for timefractional heat conduction in a composite medium using carbon-carbon experimental data, Physica B 456 (2015) 9–15. [20] J.L. Battaglia, O. Cois, L. Puigsegur, A. Oustaloup, Solving an inverse heat conduction problem using a non-integer identified model, Int. J. Heat Mass Transf. 44 (2001) 2671–2680. [21] D.A. Murio, Time fractional IHCP with Caputo fractional derivatives, Comput. Math. Appl. 56 (2008) 2371–2381. [22] Z. Qian, Optimal modified method for a fractional-diffusion inverse heat conduction problem, Inverse Probl. Sci. Eng. 18 (2010) 521–533. [23] D.A. Murio, Stable numerical evaluation of Gruünwald-Letnikov fractional derivatives applied to a fractional IHCP, Inverse Probl. Sci. Eng. 17 (2009) 229–243. [24] H. Wei, W. Chen, H.G. Sun, X.C. Li, A coupled method for inverse source problem of spatial fractional anomalous diffusion equations, Inverse Probl. Sci. Eng. 18 (2010) 945–956.

7