- Email: [email protected]
Numerical algorithms to estimate relaxation parameters and Caputo fractional derivative for a fractional thermal wave model in spherical composite medium
Applied Mathematics and Computation 274 (2016) 106–118
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical algorithms to estimate relaxation parameters and Caputo fractional derivative for a fractional thermal wave model in spherical composite medium Bo Yu a, Xiaoyun Jiang a,∗, Chu Wang b a b
School of Mathematics, Shandong University, Jinan 250100, China Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, USA
a r t i c l e
i n f o
Keywords: Caputo fractional derivative Fractional thermal wave model Spherical composite medium Finite difference Inverse problem
a b s t r a c t In this paper, we formulate a fractional thermal wave model for a bi-layered spherical tissue. Implicit finite difference method is employed to obtain the solution of the direct problem. The inverse analysis for simultaneously estimating the Caputo fractional derivative and the relaxation time parameters is implemented by means of the Levenberg–Marquardt method. Compared with the experimental data, we can obviously find out that the estimated temperature increase values are excellently consistent with the measured temperature increase values in the experiment. We have also discussed the effect of the fractional derivative, the relaxation time parameters, the initial guess as well as the sensitivity problem. All the results show that the proposed fractional thermal wave model is efficient and accurate in modeling the heat transfer in the hyperthermia experiment, and the proposed numerical method for simultaneously estimating multiple parameters for the fractional thermal wave model in a spherical composite medium is effective. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The classical Fourier’s heat-conduction equation leads to the solutions exhibiting an infinite velocity of thermal propagation and the thermal effect is felt instantaneously. However, many researchers have questioned the validity of the classical Fourier’s model for the emergence of many new fields involving heat conduction. In fact, heat is always found to propagate at a finite speed. Now, it is acknowledged that the classical Fourier’s heat-conduction model will be invalid in such situations as involving very large temperature gradients, very low temperature nearing to absolute zero or an ultra-high heating speed [1,2]. Since Cattaneo [3] and Verenotte [4] independently introduced a time-dependent relaxation parameter to the Fourier’s law and gave a generalized non-Fourier heat conduction equation, the non-Fourier heat conduction model has gradually become attractive in practical engineering problems. For example, in biological tissues, for the tissues are highly non-homogenous and the velocity of heat transfer in tissues is limited, the heterogeneous and non-isotropic nature yields normally a strong non-instantaneous response between heat flux and temperature gradient in non-equilibrium heat transport [5]. Many literatures have shown that the thermal behavior in non-homogenous media needs a relaxation time to accumulate energy for significant heat transfer. Unfortunately, the relaxation time of the biological tissues is still not clearly known for the
∗
Corresponding author. Tel.: +86 13011720190. E-mail addresses: [email protected] (B. Yu), [email protected] (X. Jiang), [email protected] (C. Wang).
http://dx.doi.org/10.1016/j.amc.2015.10.081 0096-3003/© 2015 Elsevier Inc. All rights reserved.
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
107
complexity of the biological components. Kaminski [6] presented that the relaxation time was in the range of 20 − 30 s. Mitra et al. [7] found the relaxation time for the processed meat was about 15.5 s. Roetzel et al. [8] investigated the non-Fourier behavior in the processed meat and showed the value of the relaxation time was about 1.77 s. Dual-phase-lag heat conduction model was used by Antaki [9] to offer a new interpretation for the experimental evidence of the non-Fourier conduction in the processed meat by introducing the phase lag times of the heat flux and the temperature gradient, the results showed that the relaxation time for the heat flux was about 16 s, 14 s and the relaxation time for the temperature gradient was about 0.043 s, 0.056 s for the first and the third experiment, respectively. In the past few decades, fractional calculus has proved its efficiency in modeling the anomalous behaviors observed in various fields of science and engineering [10–19]. Fractional calculus provides an excellent tool for the description of memory and hereditary properties of various materials and processes such as physical and biological processes [20–26]. Tan et al. [23,24] developed an anomalous sub-diffusion model for calcium sparks in cardiac myocytes. Jiang and Qi [25] derived a fractional thermal wave model of skin burns caused by spatial heating, their study demonstrated that fractional models could provide a unified approach to examine the heat transfer in biological tissue. Ghazizadeh et al. [26] studied a fractional single-phase-lag heat equation for processed meat, their results showed that fractional single-phase-lag model could give the same temperature distribution as the dual-phase-lag model. Recently, fractional heat conduction equation in a composite medium has attracted more and more researchers’ attention. Ilic et al. [27] studied a one-dimensional space fractional diffusion equation in a composite medium consisting of two layers. Jiang and Chen [28] presented analytical and numerical methods to solve the time fractional anomalous thermal diffusion equation in M-layers composite medium. Inverse problem is an ill-posed problem, various studies have been presented such as [1,29–35]. However, for the fractional inverse problems, there are relatively few literatures. Murio [36] studied the numerical solution of the time fractional inverse heat conduction problem on a finite slab. Cheng et al. [37] studied the uniqueness result in an inverse problem for a one-dimensional fractional diffusion equation. Ghazizadeh et al. [26] 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. Tian et al. [38] extended two regularization methods to solve the inverse space fractional diffusion equation which is ill-posed. Yu et al.[39] proposed a numerical method to estimate the unknown order of the Riemann–Liouville fractional derivative for a fractional Stokes’ first problem for a heated generalized second grade fluid. However, to the best of the authors’ knowledge, we have not found any literature on the inverse analysis of the fractional thermal wave equation in a spherical composite medium. This motivates us to get down to the present study. The following work is arranged as follows. In Section 2, we introduce some mathematical and physical background and formulate a fractional thermal wave model for a bi-layered spherical tissue. With change of variable, we transform the original problem into a new problem in the rectangular coordinate system. In Section 3, implicit finite difference method is employed to obtain the solution of the direct problem. Levenberg–Marquardt method based on fractional derivative is adopted to obtain the optimal estimation of the fractional derivative and the relaxation time parameters in Section 4. In Section 5, compared with the experimental data, the optimal estimations of the fractional derivative and the relaxation time parameters are obtained and some further discussions are made for the corresponding theoretical analysis. Finally, some conclusions are given in Section 6.
2. Mathematical formulation To mathematically predict the tissue temperature distributions during hyperthermia treatments, Andrä et al. [40] modeled a breast carcinomas surrounded by healthy tissues as a small solid sphere with the radius R and constant heat power density Q, and obtained the temperature distribution as function of time by the technique of Laplace transformation. In the magnetic hyperthermia, magnetic particles are injected into the center of the tumor and distribute homogeneously. The magnetic particles become heat sources in the tissue when excited by an alternating magnetic field. 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 center of the sphere and time t. From the experimental study [40], 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. From [40], we can obviously 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 demonstrates that the classical bio-heat transfer equation may not completely describe the thermal behavior captured in the experiment [41]. Liu and Lin [41] proposed a dual-phase-lag heat conduction model and investigated the corresponding inverse problem to estimate the phase lag times by a hybrid scheme based on the numerical inversions of the Laplace transform using the experimental data. The dual-phase-lag model introduces more complexities and more difficulties. Nevertheless, they provided a novel mathematical model and an efficient method for such problem. Ghazizadeh et al. [26] studied a fractional single-phase-lag heat equation for processed meat, their results showed that fractional single-phase-lag model could give the same temperature distribution as the dual-phase-lag model. In many physical situations, heat conduction has been treated according to the classical Fourier’s law
q(r, t ) = −k∇ T (r, t ),
(1)
108
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
where q is the heat flux, T is the temperature and k is the thermal conductivity. Cattaneo [3] and Verenotte [4] independently introduced a time-dependent relaxation parameter to the Fourier’s law
q(r, t + τ ) = −k∇ T (r, t ),
(2)
where τ is the thermal relaxation time parameter. Applying fractional Taylor series formula [42], we get [25,26]
q(r, t ) + τ α
∂ α q(r, t ) = −k∇ T (r, t ), 0 < α ≤ 1. ∂tα
(3)
Applying the Laplace transform and the inverse Laplace transform [10], we can obtain
k q(r, t ) = − α τ
t
0
(t − τ )α ∇ T (r, τ )dτ . (t − τ )α−1 Eα,α − τα
(4)
Considering the energy balance equation
ρc
∂ T (r, t ) = −div(q(r, t )) + Q, ∂t
(5)
where ρ , c, k are the density, specific heat and thermal conductivity, respectively, Q is the constant heat power density. After eliminating the heat flux term q, we can obtain
ρc
∂ 1+α T (r, t ) ∂ T (r, t ) + τα ∂t ∂ t 1+α
= k∇ 2 T (r, t ) + Q.
(6)
Based on the above analysis, we here formulate a fractional thermal wave 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
ρ1 c1 ρ2 c2
∂ 1+α T1 ∂ T1 + τ1α ∂t ∂ t 1+α
=
∂ 1+α T2 ∂ T2 + τ2α ∂t ∂ t 1+α
k1 ∂ ∂ T1 r2 ∂r r2 ∂ r
k2 ∂ ∂ T2 = 2 r2 ∂r r ∂r
+ Q, 0 ≤ r ≤ R, 0 < α ≤ 1,
(7)
, R < r ≤ a.
(8)
Here, ρ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. τ 1 and τ 2 are the relaxation time parameters. a denotes the boundary of tissue. T0 is regarded β as the arterial temperature. Q is the constant heat power density. ∂ β is the Caputo time-fractional operator defined as [10]
⎧ ⎨
1
∂ β f (t ) = (n − β) ∂tβ ⎩ (n) f (t ),
∂t
f (n) (τ )
t 0
(t − τ )β −n+1
dτ ,
n − 1 < β < n,
β = n.
We here consider the following boundary conditions
T1 (0, t ) is finite,
(9)
T2 (a, t ) = T0 ,
(10)
T1 (R, t ) = T2 (R, t ),
(11)
k1
∂ T1 (R, t ) ∂ T2 (R, t ) = k2 , ∂r ∂r
(12)
and the initial conditions
Ti (r, 0) = T0 ,
∂ Ti (r, 0) = 0, i = 1, 2. ∂t
(13)
In order to obtain the solution of the above problem, we introduce a new dependent variable θ defined as
θ = r(T − T0 ).
(14)
Hence, the above problem can be transformed into a new problem in the rectangular coordinate system, that is, (7) and (8) can be rewritten in terms of θ , respectively, as
ρ1 c1
∂ 1+α θ1 ∂θ1 + τ1α ∂t ∂ t 1+α
= k1
∂ 2 θ1 + Qr, 0 ≤ r ≤ R, 0 < α ≤ 1, ∂ r2
(15)
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
ρ2 c2
∂ 1+α θ2 ∂θ2 + τ2α ∂t ∂ t 1+α
= k2
109
∂ 2 θ2 , R < r ≤ a. ∂ r2
(16)
The boundary and initial conditions are rewritten as
θ1 (0, t ) = 0,
(17)
θ2 (a, t ) = 0,
(18)
θ1 (R, t ) = θ2 (R, t ), ∂θ1 (R, t ) θ1 (R, t ) ∂θ2 (R, t ) θ2 (R, t ) − − k1 = k2 , ∂r R ∂r R ∂θi (r, 0) θi (r, 0) = 0, = 0, i = 1, 2. ∂t
(19) (20) (21)
Obviously, we observe that the temperature increase T as a function of time for different reduced distances r/R from the composite center is then obtained by
T = T − T0 =
θ r
.
(22)
3. Finite difference method for the direct problem In this section, we will adopt the finite difference method to solve the direct problem. The direct problem is to solve the temperature distribution θi (r, t )(i = 1, 2) from (15)–(21) where ρ 1 , c1 , ρ 2 , c2 , k1 , k2 , P, R, a, τ 1 ,τ 2 and α are already known. Ghazizadeh et al. [43] studied explicit and implicit finite difference schemes for the fractional Cattaneo equation. In this section, we will adopt the implicit finite difference schemes to solve the bi-layered composite medium problem. Firstly, we 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 = tn − tn−1 , n = 1, 2, . . . , N. Secondly, for 1 < β = 1 + α < 2, we employ the following discrete schemes
t=tn n
∂ β θ (r, t )
( t )−β = ω θ n−k+1 − 2θi,n−k + θi,n−k−1 , j j (3 − β) k=1 k i, j ∂ t β r=r
(23)
t=tn + θi,n−2 3θi,nj − 4θi,n−1 ∂θ (r, t )
j j , = ∂ t r=r 2 t
(24)
t=tn θi,nj+1 − 2θi,nj + θi,nj−1 ∂ 2 θ (r, t )
= , 2
∂r ( ri )2 r=r
(25)
i, j
i, j
i, j
where ωk = k2−β − (k − 1)2−β , k = 1, 2, . . .. Hence, we obtain the implicit finite difference schemes for (15) and (16):
ρ1 c1
n−1 n−2 n 3θ1, − 4θ1, + θ1, j j j
2 t
ρ2 c2
3θ
n 2, j
− 4θ
n−1 2, j
+ τα 1
+θ
n−2 2, j
2 t
+ τ2α
k+1 k+1 n θ1,k+1 − 2θ1, + θ1, n−k+1
( t )−β j+1 j j−1 n−k n−k−1 ωk θ1, j − 2θ1, j + θ1, j = k1 + Qr1, j , (3 − β) k=1 ( r1 )2 ( t )−β (3 − β)
n
θ
n−k n−k−1 ωk θ2,n−k+1 − 2θ2, + θ2, = k2 j j j
k=1
(26) k+1 2, j+1
− 2θ
k+1 + 2, j r2 2
( )
θ
k+1 2, j−1
.
(27)
The boundary and initial conditions can be discretized as k+1 θ1,0 = 0,
(28) (29)
θ
k+1 2,M
= 0,
θ
k+1 1,M
k+1 = 2,0 , k+1 k+1 − 1,M−1 1,M
k1
θ
θ
θi,0j = 0,
θ r 1
−
θ
k+1 1,M
R
= k2
θ
−θ r 2
k+1 2,1
k+1 2,0
−
θ
k+1 2,0
R
,
θi,0j − θi,−1j = 0, i = 1, 2. t
By solving (26)–(32), we can easily obtain the solution of the direct problem.
(30) (31)
(32)
110
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
4. Inverse analysis For the inverse problem, the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) are investigated. The measured temperature increase T as a function of time for different reduced distances r/R from the composite center are available [40]. The solution of the inverse problem is obtained through the minimization of the least squares norm in the following form
S(P) = ( T − Tc )T ( T − Tc ),
(33)
where T is the measured temperature increase vector in the experiment
T = [ T (t1,1 ), T (t1,2 ), . . . , T (t1,n ), T (t2,1 ), T (t2,2 ), . . . , T (t2,n )]T ,
(34)
Tc is the calculated temperature increase vector from the direct problem
Tc = [ Tc (t1,1 , P), Tc (t1,2 , P), . . . , Tc (t1,n , P), Tc (t2,1 , P), Tc (t2,2 , P), . . . , Tc (t2,n , P)]T ,
(35)
ti, j denotes the measured time nodes, n is the total number of the measured data from the experiment, and P = (α , τ1 , τ2 )T is regarded as being unknown. 4.1. Levenberg–Marquardt method For estimating the unknown parameter P = (α , τ1 , τ2 )T , Levenberg–Marquardt iterative method is introduced to minimize the above least squares norm [44]. The iterative algorithm for finding the unknown parameter is
Pk+1 = Pk + [(Jk )T Jk + μk k ]−1 (Jk )T ( T − Tkc ), where
Jk
(36)
is the fractional sensitivity matrix (Jacobian matrix in the fractional sense),
k is taken as
k = diag((Jk )T Jk ),
μk
is a positive scalar (damping parameter), (37)
and k is the number of iterations. The elements of the fractional sensitivity matrix are the first derivative of the estimated temperature increases at the measured time nodes ti, j with respect to the unknown parameter Pl (l = 1, 2, 3), that is, J = [Jα , Jτ1 , Jτ2 ]. The iteration begins with an appropriate initial guess P0 = (α 0 , τ10 , τ20 )T for the unknown parameter P = (α , τ1 , τ2 )T . The stopping criterion of the iterative algorithm is
k+1
P − Pk < ,
(38)
where is a chosen tolerance. It is worthy to mention that the three parameters α , τ 1 , τ 2 should vary in their physical range in the iteration process, that is, 0 < α ≤ 1, τ 1 ≥ 0, τ 2 ≥ 0. 4.2. Calculation of the fractional sensitivity matrix 4.2.1. Calculation of Jα The Jacobian equations with respect to α are obtained through differentiating (15) and (16) with respect to the unknown parameter α , that is,
ρi ci
∂ 1+α θi ∂ ∂ 1+α θi ∂ 2 Jαi ∂ Jαi α α + τi ln τi 1+α + τi = k , i = 1, 2, i ∂t ∂α ∂ t 1+α ∂t ∂ r2
(39)
∂θ (r,t )
i where Jα i (r, t ) = ∂α , i = 1, 2. The boundary and initial conditions are
Jα 1 (0, t ) = 0,
(40)
Jα 2 (a, t ) = 0,
(41)
Jα 1 (R, t ) = Jα 2 (R, t ),
k1
∂ Jα1 (R, t ) Jα1 (R, t ) − ∂r R
Jα i (r, 0) = 0,
(42)
= k2
∂ Jα2 (R, t ) Jα2 (R, t ) − , ∂r R
∂ Jαi (r, 0) = 0, i = 1, 2. ∂t
(43)
(44)
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
111
∂ ∂ 1+α θi (r,t ) , we have Concerning the discrete form of the term ∂α ∂ t 1+α
t=tn n n−k+1
( t )−β ∂ ∂ 1+α θi (r, t )
∂ = ω θ − 2θi,n−k + θi,n−k−1 j j
∂α ∂α (3 − β) k=1 k i, j ∂ t 1+α r=ri, j n n
n−k+1
( t )−β ∂ωk n−k+1 n−k n−k−1 n−k n−k−1 = − 2θi, j + θi, j + ωk Jα,i, j − 2Jα,i, j + Jα,i, j θ (3 − β) k=1 ∂α i, j k=1 n
∂ ( t )−β + ωk θi,n−k+1 − 2θi,n−k + θi,n−k−1 , j j j ∂α (3 − β) k=1
(45)
∂ω ∂ω ∂ ( t )−β ) = ( t )−β ( − ln t + where β = 1 + α , d1 = ∂α1 = 0, dk = ∂αk = −k2−β ln k + (k − 1)2−β ln (k − 1), k = 2, 3, . . . , ( ∂α (3−β) (3−β) −∂
∂α ( (3−β)) ) (3−β)
t )−β d ln (x) = ( is the digamma function. dx (3−β) ( − ln t + ψ(3 − β)), and ψ(x) = Solving (39)–(45) by implicit finite difference method similarly as in Section 3, we can obtain Jα .
4.2.2. Calculation of Jτ1 The Jacobian equations with respect to τ 1 are obtained through differentiating (15) and (16) with respect to the unknown parameter τ 1 :
ρ1 c1 ρ2 c2 that is,
ρ1 c1 ρ2 c2
1+α ∂ 1+α Jτ1 1 ∂ Jτ1 1 θ1 α −1 ∂ + τ1α + ατ 1 ∂t ∂ t 1+α ∂ t 1+α
∂ 1+α Jτ1 2 ∂ Jτ1 2 + τ2α ∂t ∂ t 1+α ∂ 1+α Jτ1 1 ∂ Jτ1 1 + τ1α ∂t ∂ t 1+α ∂ 1+α Jτ1 2 ∂ Jτ1 2 + τ2α ∂t ∂ t 1+α
where Jτ1 i (r, t ) =
= k1
∂ 2 Jτ1 1 , 0 ≤ r ≤ R, ∂ r2
(46)
= k2
∂ 2 Jτ1 2 , R < r ≤ a, ∂ r2
(47)
= k1
∂ 2 Jτ1 1 ∂ 1+α θ1 − ατ1α −1 , 0 ≤ r ≤ R, 2 ∂r ∂ t 1+α
(48)
= k2
∂ 2 Jτ1 2 , R < r ≤ a, ∂ r2
(49)
∂θi (r,t ) ∂τ1 , i = 1, 2.
The boundary and initial conditions are
Jτ1 1 (0, t ) = 0,
(50)
Jτ1 2 (a, t ) = 0,
(51)
Jτ1 1 (R, t ) = Jτ1 2 (R, t ),
(52)
∂ Jτ1 1 (R, t ) Jτ1 1 (R, t ) ∂ Jτ1 2 (R, t ) Jτ1 2 (R, t ) = k2 , − − ∂r R ∂r R ∂ Jτ1 i (r, 0) = 0, i = 1, 2. Jτ1 i (r, 0) = 0, ∂t k1
(53) (54)
Solving (48)–(54) by implicit finite difference method similarly as in Section 3, we can obtain Jτ1 . 4.2.3. Calculation of Jτ2 The Jacobian equations with respect to τ 2 are obtained through differentiating (15) and (16) with respect to the unknown parameter τ 2 :
ρ1 c1 ρ2 c2
∂ 1+α Jτ2 1 ∂ Jτ2 1 + τ1α ∂t ∂ t 1+α
= k1
∂ 2 Jτ2 1 , 0 ≤ r ≤ R, ∂ r2
∂ 1+α Jτ2 2 ∂ Jτ2 2 ∂ 1+α θ2 + τ2α + ατ2α −1 ∂t ∂ t 1+α ∂ t 1+α
= k2
∂ 2 Jτ2 2 , R < r ≤ a, ∂ r2
(55)
(56)
112
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
that is,
ρ1 c1 ρ2 c2
∂ 1+α Jτ2 1 ∂ Jτ2 1 + τ1α ∂t ∂ t 1+α ∂ 1+α Jτ2 2 ∂ Jτ2 2 + τ2α ∂t ∂ t 1+α
where Jτ2 i (r, t ) =
= k1
∂ 2 Jτ2 1 , 0 ≤ r ≤ R, ∂ r2
(57)
= k2
∂ 2 Jτ2 2 ∂ 1+α θ2 − ατ2α −1 , R < r ≤ a, 2 ∂r ∂ t 1+α
(58)
∂θi (r,t ) ∂τ2 , i = 1, 2.
The boundary and initial conditions are
Jτ2 1 (0, t ) = 0,
(59)
Jτ2 2 (a, t ) = 0,
(60)
Jτ2 1 (R, t ) = Jτ1 2 (R, t ),
(61)
k1
∂ Jτ2 1 (R, t ) Jτ2 1 (R, t ) − ∂r R
Jτ2 i (r, 0) = 0,
= k2
∂ Jτ2 2 (R, t ) Jτ2 2 (R, t ) − , ∂r R
∂ Jτ2 i (r, 0) = 0, i = 1, 2. ∂t
(62)
(63)
Solving (57)–(63) by implicit finite difference method similarly as in Section 3, we can obtain Jτ2 . 4.3. Computational algorithms Suppose that the measured temperature increase vector T in the experiment is given, the initial guess P0 = (α 0 , τ10 , τ20 )T is available for the unknown parameter P = (α , τ1 , τ2 )T , choose a value for μ0 , for example, μ0 = 10, and set k = 0, Step 1. Solve the direct problem (26)–(32) with the available estimate Pk in order to obtain the temperature increase vector
Tkc = [ Tc (t1,1 , Pk ), Tc (t1,2 , Pk ), . . . , Tc (t1,n , Pk ), Tc (t2,1 , Pk ), Tc (t2,2 , Pk ), . . . , Tc (t2,n , Pk )]T ; Step 2. Step 3. Step 4. Step 5.
(64)
S(Pk )
Compute ; Compute the fractional sensitivity matrix Jk , and compute k ; Compute Pk = [(Jk )T Jk + μk k ]−1 (Jk )T ( T − Tkc ), and compute the new estimate Pk+1 = Pk + Pk ; , and Solve the direct problem with the new estimate Pk+1 in order to obtain the new temperature increase vector Tk+1 c compute S(Pk+1 ); Step 6. If S(Pk+1 ) ≥ S(Pk ), replace μk by 10μk , and return to Step 4; Step 7. If S(Pk+1 ) < S(Pk ), accept the new estimate Pk+1 , and replace μk by μk /10; Step 8. Check the stopping criteria. If satisfied, stop the iterative; otherwise, replace k by k + 1, and return to Step 3. 5. Results and discussions Andrä et al. [40] simulated the thermal behavior of localized magnetic hyperthermia applied to a breast carcinoma. They 3 considered a small sphere tumor of radius R = 0.00315 m with a constant power density Q = 6150, 000 W/m embedded in ex3 3 tended muscle tissue. The thermal parameters were taken to be ρ1 = 1, 660, 000 g/m , ρ2 = 1, 000, 000 g/m , c1 = 2.54 J/(g K), c2 = 3.72 J/(g K), k1 = 0.778 W/(K m), k2 = 0.642 W/(K m). 5.1. Optimal estimation of the parameter P = (α , τ1 , τ2 )T The measured temperature increase values which are used by us are extracted from the symbols at reduced distances r/R = 0.65, 1.10, the error of T and t are within the extent of the used symbols. For parameter a, we here choose a = 11R. Chosen the tolerance = 10−4 and an initial guess P0 = (0.8, 10, 10)T , the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) are obtained by the Levenberg–Marquardt method. The computational results are listed in Table 1. From Table 1, we can see that the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) are obtained after 12 iterations. Fig. 1 demonstrates the comparison between the calculated temperature increase values with the optimal estimations of the fractional derivative and the relaxation time parameters and the measured temperature increase values in the experiment at
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
113
Table 1 Optimal estimation for parameter P = (α , τ1 , τ2 )T . Initial guess P0
Number of iterations
Estimated value P
(0.8, 10, 10)T
12
(0.4381, 0.5717, 1.7914)T
35 Experimental data [40] Present result
Temperature Increase
T
30
=0.4381 =0.5717 1
25
r/R=0.65
=1.7914
2
20 r/R=1.10 15 r/R=1.39 10
r/R=1.68 r/R=1.98
5
0
0
50
100
150 t
200
250
300
Fig. 1. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10, 1.39, 1.68, 1.98.
Table 2 Optimal estimation for different initial guess. Initial guess P0
Number of iterations
Estimated value P
(0.9, 10, 10)T (0.4, 10, 10)T (0.8, 1, 10)T (0.8, 20, 10)T (0.8, 10, 1)T (0.8, 10, 20)T
12 10 13 10 11 12
(0.438139, 0.571676, 1.791362)T (0.438141, 0.571681, 1.791409)T (0.438137, 0.571670, 1.791315)T (0.438141, 0.571680, 1.791403)T (0.438137, 0.571670, 1.791316)T (0.438140, 0.571677, 1.791370)T
the reduced distances r/R = 0.65, 1.10, 1.39, 1.68, 1.98. From Fig. 1, we can obviously observe that the calculated temperature increase values are in excellently consistent with the measured temperature increase values in the experiment when α , τ 1 , τ 2 are chosen as the optimal estimated values in Table 1, which is much more better than the result of [40]. This fact demonstrates that the proposed fractional thermal wave model is efficient and accurate in modeling the heat transfer in the hyperthermia experiment, and the proposed numerical method is effective for solving the corresponding the fractional inverse thermal wave problem in a composite medium.
5.2. Effect of the initial guess Table 2 demonstrates the optimal estimations of the fractional derivative α and the relaxation time parameters τ 1 and τ 2 as well as the number of iterations for different initial guess P0 = (α 0 , τ10 , τ20 )T . From Table 2, we can obviously find out that different initial guess P0 has little impact on the final optimal estimation. This also demonstrates that Levenberg–Marquardt method based on fractional derivative is efficient for solving the fractional inverse thermal wave problem in a composite medium.
114
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
35 r/R=0.65 =0.5717, 2=1.7914
30
1
Temperature Increase
T
=0.4381 =0.1
25
=0.9
20 r/R=1.10 15
10 Experimental data [40] Present result( =0.4381) Present result( =0.1) Present result( =0.9)
5
0
0
50
100
150 t
200
250
300
Fig. 2. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10 for different α .
5.3. Effect of the fractional derivative α Fractional derivative α is a very important parameter in determining the temperature distribution. Fig. 2 demonstrates the comparison between the measured temperature increase values in the experiment and the calculated temperature increase values at r/R = 0.65, 1.10 when α is chosen as the optimal estimated value 0.4381 in Table 1 and another values 0.1 and 0.9, τ 1 , τ 2 are chosen as the optimal estimated values in Table 1. From Fig. 2, we can obviously find out that the curves are higher when α is larger, and lower when α is smaller, and the curves protrude to upper left slower when α is becoming larger. Additionally, we can observe that the calculated temperature increase values agree much better with the measured temperature increase values in the experiment when α is chosen as the optimal estimated value 0.4381 than α is chosen as 0.1 and 0.9. 5.4. Effect of the relaxation time parameters τ 1 and τ 2 Relaxation time parameters τ 1 and τ 2 are important parameters in determining the temperature distribution, too. In this section, we will discuss the effect of the relaxation time parameters τ 1 and τ 2 for the temperature increase T. Fig. 3 demonstrates the comparison between the measured temperature increase values in the experiment and the calculated temperature increase values at r/R = 0.65, 1.10 when α , τ 1 , τ 2 are chosen as the optimal estimated values in Table 1 and another value τ1 = 10. From Fig. 3, we can obviously find out that the curves on the lower left are higher when τ 1 is smaller, and lower when τ 1 is larger. But it is just contrary on the top right. Additionally, we can observe that the calculated temperature increase values agree much better with the measured temperature increase values in the experiment when τ 1 is chosen as the optimal estimated value 0.5717 than τ 1 is chosen as 10. Fig. 4 demonstrates the comparison between the measured temperature increase values in the experiment and the calculated temperature increase values at r/R = 0.65, 1.10 when α , τ 1 , τ 2 are chosen as the optimal estimated values in Table 1 and another values τ2 = 10. From Fig. 4, we can obviously find out that the curves are higher when τ 2 is smaller, and lower when τ 2 is larger. Additionally, we can observe that the calculated temperature increase values agree much better with the measured temperature increase values in the experiment when τ 2 is chosen as the optimal estimated value 1.7914 than τ 2 is chosen as 10. 5.5. Sensitivity problem The sensitivity coefficients JPi (i = 1, 2, 3) are the measurements of the sensitivity of the estimated temperature increase T with respect to the changes in the parameter Pi at the reduced distances r/R = 0.65, 1.10. Here we introduce the relative
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
115
35 =0.4381, 2=1.7914
30
r/R=0.65
Temperature Increase
T
=0.5717
1
25
20 r/R=1.10 15
10 Experimental data [40] Present result( 1=0.5717)
5
=10
Present result( 1=10)
1
0
0
50
100
150 t
200
250
300
Fig. 3. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10 for different τ 1 .
35 =0.4381, 1=0.5717
30
r/R=0.65
=1.7914
Temperature Increase
T
2
25
20 r/R=1.10 15
10 Experimental data [40] =1.7914)
=10
2
5
2
Present result( 2=10) 0
0
50
100
150 t
200
250
300
Fig. 4. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10 for different τ 2 .
116
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
0.06 J J1 J
r/R=0.65
0.05
2
Sensitivity Coefficients
0.04 J
0.03 0.02 J
0.01
2
J
1
0 −0.01 −0.02
0
50
100
150 t
200
250
300
Fig. 5. Sensitivity coefficients for α , τ 1 , τ 2 at r/R = 0.65.
0.05 J r/R=1.10
J1 J2
Sensitivity Coefficients
0.04
0.03 J 0.02
0.01 J2
J1
0
−0.01
0
50
100
150 t
200
250
300
Fig. 6. Sensitivity coefficients for α , τ 1 , τ 2 at r/R = 1.10.
sensitivity coefficients
JPi =
P i ∂ T , i = 1, 2, 3. T0 ∂ Pi
(65)
Figs. 5 and 6 show the sensitivity coefficients Jα , Jτ1 and Jτ2 at r/R = 0.65, 1.10, respectively. The sensitivity coefficients are linear independency, it is possible to have a correct estimation for the unknown parameters.
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
117
6. Conclusions In this paper, a fractional thermal wave model is formulated for a bi-layered spherical tissue. The solution of the direct problem is obtained by implicit finite difference method. For the inverse problem, we obtain the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) simultaneously by means of the Levenberg–Marquardt method. Compared with the experimental data, we can obviously find out that the calculated temperature increase values are in excellently consistent with the measured temperature increase values in the experiment, which demonstrates that the proposed fractional thermal wave model is efficient and accurate in modeling the heat transfer in the hyperthermia experiment, and Levenberg– Marquardt method for solving the fractional inverse thermal wave problem in a spherical composite medium is effective and valid. From the viewpoint of methodology, this paper provides an efficient numerical method for simultaneously estimating multiple parameters in fractional thermal wave model in a spherical composite medium. Acknowledgment The project was supported by the National Natural Science Foundation of China (11472161, 11102102 and 91130017). References [1] P.T. Hsu, Y.H. Chu, An inverse non-Fourier heat conduction problem approach for estimating the boundary condition in electronic device, Appl. Math. Model. 28 (2004) 639–652. [2] V. Vishwakarma, A.K. Das, P.K. Das, Analysis of non-Fourier heat conduction using smoothed particle hydrodynamics, Appl. Therm. Eng. 31 (2011) 2963– 2970. [3] C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, C.R. Acad. Sci. 247 (1958) 431–433. [4] P. Vernotte, Paradoxes in the continuous theory of the heat equation , C.R. Acad. Sci. 246 (1958) 154-3. [5] L. Wang, F. Fan, Modeling bioheat transport at macroscale, J. Heat Transf. 133 (2011) 011010. [6] W. Kaminski, Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure, J. Heat Transf. 112 (1990) 555–560. [7] K. Mitra, S. Kumar, A. Vedevarz, M.K. Moallemi, Experimental evidence of hyperbolic heat conduction in processed meat, J. Heat Transf. 117 (1995) 568–573. [8] W. Roetzel, N. Putra, S.K. Das, Experiment and analysis for non-Fourier conduction in materials with non-homogeneous inner structure, Int. J. Therm. Sci. 42 (2003) 541–552. [9] P.J. Antaki, New interpretation of non-Fourier heat conduction in processed meat, J. Heat Transf. 127 (2005) 189–193. [10] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [11] F. Liu, C. Yang, K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math. 231 (2009) 160–176. [12] X.Y. Jiang, M.Y. Xu, H.T. Qi, The fractional diffusion model with an absorption term and modified Fick’s law for non-local transport processes , Nonlinear Anal. Real World Appl. 11 (2010) 262–269. [13] C. Li, Z. Zhao, Y. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62 (2011) 855–875. [14] L. Zheng, F. Zhao, X. Zhang, Exact solutions for generalized Maxwell fluid flow due to oscillatory and constantly accelerating plate, Nonlinear Anal. Real World Appl. 11 (2010) 3744–3751. [15] B. Jin, R. Lazarov, Z. Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations, SIAM J. Numer. Anal. 51 (2013) 445–466. [16] B. Jin, R. Lazarov, J. Pasciak, Z. Zhou, Error analysis of a finite element method for the space-fractional parabolic equation, SIAM J. Numer. Anal. 52 (2014) 2272–2294. [17] F. Zeng, C. Li, F. Liu, I. Turner, The use of finite difference/element approaches for solving the time-fractional subdiffusion equation, SIAM J. Sci. Comput. 35 (2013) A2976–A3000. [18] F. Zeng, F. Liu, C. Li, K. Burrage, I. Turner, V. Anh, A Crank–Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction– diffusion equation, SIAM J. Numer. Anal. 52 (2014) 2599–2622. [19] B. Yu, X.Y. Jiang, H. Xu, A novel compact numerical method for solving the two-dimensional non-linear fractional reaction–subdiffusion equation 68 (2015) 923–950. [20] A. Eke, P. Herman, L. Kocsis, L.R. Kozak, Fractal characterization of complexity in temporal physiological signals, Physiol. Meas. 23 (2002) R1. [21] 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. Phys. A Math. Gen. 37 (2004) R161. [22] R.L. Magin, Fractional calculus models of complex dynamics in biological tissues, Comput. Math. Appl. 59 (2010) 1586–1593. [23] W. Tan, C. Fu, C. Fu, W. Xie, H. Cheng, An anomalous subdiffusion model for calcium spark in cardiac myocytes, Appl. Phys. Lett. 91 (2007) 3901. [24] K. Li, C. Fu, H. Cheng, W. Tan, Anomalous subdiffusion of calcium spark in cardiac myocytes, Cell. Mol. Bioeng. 4 (2011) 457–465. [25] X.Y. Jiang, H.T. Qi, Thermal wave model of bioheat transfer with modified Riemann–Liouville fractional derivative, J. Phys. A. Math. Theor. 45 (2012) 485101. [26] 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. [27] 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. [28] X.Y. Jiang, S.Z. Chen, Analytical and numerical solutions of time fractional anomalous thermal diffusion equation in composite medium, ZAMM Z. Angew. Math. Mech. 95 (2015) 156–164. [29] G. Fourestey, M. Moubachir, Solving inverse problems involving the Navier–Stokes equations discretized by a Lagrange–Galerkin method, Comput. Methods Appl. Mech. Eng. 194 (2005) 877–906. [30] A.G. Fatullayev, E. Can, N. Gasilov, Comparing numerical methods for inverse coefficient problem in parabolic equation, Appl. Math. Comput. 179 (2006) 567–571. [31] M. Ebrahimian, R. Pourgholi, M. Emamjome, P. Reihani, A numerical solution of an inverse parabolic problem with unknown boundary conditions, Appl. Math. Comput. 189 (2007) 228–234. [32] A. Shidfar, J. Damirchi, P. Reihani, An stable numerical algorithm for identifying the solution of an inverse problem, Appl. Math. Comput. 190 (2007) 231–236. [33] U. Albocher, A.A. Oberai, P.E. Barbone, I. Harari, Adjoint-weighted equation for inverse problems of incompressible plane-stress elasticity, Comput. Methods Appl. Mech. Eng. 198 (2009) 2412–2420. [34] G. Deolmi, F. Marcuzzi, A parabolic inverse convection–diffusion–reaction problem solved using space-time localization and adaptivity, Appl. Math. Comput. 219 (2013) 8435–8454. [35] Y.F. Wang, K.G. Reyes, K.A. Brown, C.A. Mirkin, W.B. Powell, Nested-batch-mode learning and stochastic optimization with an application to sequential multistage testing in materials science, SIAM J. Sci. Comput. 37 (3) (2015) B361–B381.
118
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
[36] D.A. Murio, Time fractional IHCP with Caputo fractional derivatives, Comput. Math. Appl. 56 (2008) 2371–2381. [37] J. Cheng, J. Nakagawa, M. Yamamoto, T. Yamazaki, Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation, Inverse Probl. 25 (2009) 115002. [38] W. Tian, C. Li, W. Deng, Y. Wu, Regularization methods for unknown source in space fractional diffusion equation, Math. Comput. Simul. 85 (2012) 45–56. [39] X. B.Yu, 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 (2) (2015) 153–161. [40] 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. [41] K.C. Liu, C.T. Lin, Solution of an inverse heat conduction problem in a bi-layered spherical tissue, Numer. Heat Transf. Part A Appl. 58 (2010) 802–818. [42] Z.M. Odibat, N.T. Shawagfeh, Generalized Taylor’s formula, Appl. Math. Comput. 186 (2007) 286–293. [43] H.R. Ghazizadeh, M. Maerefat, A. Azimi, Explicit and implicit finite difference schemes for fractional Cattaneo equation, J. Comput. Phys. 229 (2010) 7042– 7057. [44] M.N. Özisik, Inverse Heat Transfer: Fundamentals and Applications, CRC Press, 2000.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical algorithms to estimate relaxation parameters and Caputo fractional derivative for a fractional thermal wave model in spherical composite medium Bo Yu a, Xiaoyun Jiang a,∗, Chu Wang b a b
School of Mathematics, Shandong University, Jinan 250100, China Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, USA
a r t i c l e
i n f o
Keywords: Caputo fractional derivative Fractional thermal wave model Spherical composite medium Finite difference Inverse problem
a b s t r a c t In this paper, we formulate a fractional thermal wave model for a bi-layered spherical tissue. Implicit finite difference method is employed to obtain the solution of the direct problem. The inverse analysis for simultaneously estimating the Caputo fractional derivative and the relaxation time parameters is implemented by means of the Levenberg–Marquardt method. Compared with the experimental data, we can obviously find out that the estimated temperature increase values are excellently consistent with the measured temperature increase values in the experiment. We have also discussed the effect of the fractional derivative, the relaxation time parameters, the initial guess as well as the sensitivity problem. All the results show that the proposed fractional thermal wave model is efficient and accurate in modeling the heat transfer in the hyperthermia experiment, and the proposed numerical method for simultaneously estimating multiple parameters for the fractional thermal wave model in a spherical composite medium is effective. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The classical Fourier’s heat-conduction equation leads to the solutions exhibiting an infinite velocity of thermal propagation and the thermal effect is felt instantaneously. However, many researchers have questioned the validity of the classical Fourier’s model for the emergence of many new fields involving heat conduction. In fact, heat is always found to propagate at a finite speed. Now, it is acknowledged that the classical Fourier’s heat-conduction model will be invalid in such situations as involving very large temperature gradients, very low temperature nearing to absolute zero or an ultra-high heating speed [1,2]. Since Cattaneo [3] and Verenotte [4] independently introduced a time-dependent relaxation parameter to the Fourier’s law and gave a generalized non-Fourier heat conduction equation, the non-Fourier heat conduction model has gradually become attractive in practical engineering problems. For example, in biological tissues, for the tissues are highly non-homogenous and the velocity of heat transfer in tissues is limited, the heterogeneous and non-isotropic nature yields normally a strong non-instantaneous response between heat flux and temperature gradient in non-equilibrium heat transport [5]. Many literatures have shown that the thermal behavior in non-homogenous media needs a relaxation time to accumulate energy for significant heat transfer. Unfortunately, the relaxation time of the biological tissues is still not clearly known for the
∗
Corresponding author. Tel.: +86 13011720190. E-mail addresses: [email protected] (B. Yu), [email protected] (X. Jiang), [email protected] (C. Wang).
http://dx.doi.org/10.1016/j.amc.2015.10.081 0096-3003/© 2015 Elsevier Inc. All rights reserved.
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
107
complexity of the biological components. Kaminski [6] presented that the relaxation time was in the range of 20 − 30 s. Mitra et al. [7] found the relaxation time for the processed meat was about 15.5 s. Roetzel et al. [8] investigated the non-Fourier behavior in the processed meat and showed the value of the relaxation time was about 1.77 s. Dual-phase-lag heat conduction model was used by Antaki [9] to offer a new interpretation for the experimental evidence of the non-Fourier conduction in the processed meat by introducing the phase lag times of the heat flux and the temperature gradient, the results showed that the relaxation time for the heat flux was about 16 s, 14 s and the relaxation time for the temperature gradient was about 0.043 s, 0.056 s for the first and the third experiment, respectively. In the past few decades, fractional calculus has proved its efficiency in modeling the anomalous behaviors observed in various fields of science and engineering [10–19]. Fractional calculus provides an excellent tool for the description of memory and hereditary properties of various materials and processes such as physical and biological processes [20–26]. Tan et al. [23,24] developed an anomalous sub-diffusion model for calcium sparks in cardiac myocytes. Jiang and Qi [25] derived a fractional thermal wave model of skin burns caused by spatial heating, their study demonstrated that fractional models could provide a unified approach to examine the heat transfer in biological tissue. Ghazizadeh et al. [26] studied a fractional single-phase-lag heat equation for processed meat, their results showed that fractional single-phase-lag model could give the same temperature distribution as the dual-phase-lag model. Recently, fractional heat conduction equation in a composite medium has attracted more and more researchers’ attention. Ilic et al. [27] studied a one-dimensional space fractional diffusion equation in a composite medium consisting of two layers. Jiang and Chen [28] presented analytical and numerical methods to solve the time fractional anomalous thermal diffusion equation in M-layers composite medium. Inverse problem is an ill-posed problem, various studies have been presented such as [1,29–35]. However, for the fractional inverse problems, there are relatively few literatures. Murio [36] studied the numerical solution of the time fractional inverse heat conduction problem on a finite slab. Cheng et al. [37] studied the uniqueness result in an inverse problem for a one-dimensional fractional diffusion equation. Ghazizadeh et al. [26] 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. Tian et al. [38] extended two regularization methods to solve the inverse space fractional diffusion equation which is ill-posed. Yu et al.[39] proposed a numerical method to estimate the unknown order of the Riemann–Liouville fractional derivative for a fractional Stokes’ first problem for a heated generalized second grade fluid. However, to the best of the authors’ knowledge, we have not found any literature on the inverse analysis of the fractional thermal wave equation in a spherical composite medium. This motivates us to get down to the present study. The following work is arranged as follows. In Section 2, we introduce some mathematical and physical background and formulate a fractional thermal wave model for a bi-layered spherical tissue. With change of variable, we transform the original problem into a new problem in the rectangular coordinate system. In Section 3, implicit finite difference method is employed to obtain the solution of the direct problem. Levenberg–Marquardt method based on fractional derivative is adopted to obtain the optimal estimation of the fractional derivative and the relaxation time parameters in Section 4. In Section 5, compared with the experimental data, the optimal estimations of the fractional derivative and the relaxation time parameters are obtained and some further discussions are made for the corresponding theoretical analysis. Finally, some conclusions are given in Section 6.
2. Mathematical formulation To mathematically predict the tissue temperature distributions during hyperthermia treatments, Andrä et al. [40] modeled a breast carcinomas surrounded by healthy tissues as a small solid sphere with the radius R and constant heat power density Q, and obtained the temperature distribution as function of time by the technique of Laplace transformation. In the magnetic hyperthermia, magnetic particles are injected into the center of the tumor and distribute homogeneously. The magnetic particles become heat sources in the tissue when excited by an alternating magnetic field. 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 center of the sphere and time t. From the experimental study [40], 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. From [40], we can obviously 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 demonstrates that the classical bio-heat transfer equation may not completely describe the thermal behavior captured in the experiment [41]. Liu and Lin [41] proposed a dual-phase-lag heat conduction model and investigated the corresponding inverse problem to estimate the phase lag times by a hybrid scheme based on the numerical inversions of the Laplace transform using the experimental data. The dual-phase-lag model introduces more complexities and more difficulties. Nevertheless, they provided a novel mathematical model and an efficient method for such problem. Ghazizadeh et al. [26] studied a fractional single-phase-lag heat equation for processed meat, their results showed that fractional single-phase-lag model could give the same temperature distribution as the dual-phase-lag model. In many physical situations, heat conduction has been treated according to the classical Fourier’s law
q(r, t ) = −k∇ T (r, t ),
(1)
108
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
where q is the heat flux, T is the temperature and k is the thermal conductivity. Cattaneo [3] and Verenotte [4] independently introduced a time-dependent relaxation parameter to the Fourier’s law
q(r, t + τ ) = −k∇ T (r, t ),
(2)
where τ is the thermal relaxation time parameter. Applying fractional Taylor series formula [42], we get [25,26]
q(r, t ) + τ α
∂ α q(r, t ) = −k∇ T (r, t ), 0 < α ≤ 1. ∂tα
(3)
Applying the Laplace transform and the inverse Laplace transform [10], we can obtain
k q(r, t ) = − α τ
t
0
(t − τ )α ∇ T (r, τ )dτ . (t − τ )α−1 Eα,α − τα
(4)
Considering the energy balance equation
ρc
∂ T (r, t ) = −div(q(r, t )) + Q, ∂t
(5)
where ρ , c, k are the density, specific heat and thermal conductivity, respectively, Q is the constant heat power density. After eliminating the heat flux term q, we can obtain
ρc
∂ 1+α T (r, t ) ∂ T (r, t ) + τα ∂t ∂ t 1+α
= k∇ 2 T (r, t ) + Q.
(6)
Based on the above analysis, we here formulate a fractional thermal wave 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
ρ1 c1 ρ2 c2
∂ 1+α T1 ∂ T1 + τ1α ∂t ∂ t 1+α
=
∂ 1+α T2 ∂ T2 + τ2α ∂t ∂ t 1+α
k1 ∂ ∂ T1 r2 ∂r r2 ∂ r
k2 ∂ ∂ T2 = 2 r2 ∂r r ∂r
+ Q, 0 ≤ r ≤ R, 0 < α ≤ 1,
(7)
, R < r ≤ a.
(8)
Here, ρ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. τ 1 and τ 2 are the relaxation time parameters. a denotes the boundary of tissue. T0 is regarded β as the arterial temperature. Q is the constant heat power density. ∂ β is the Caputo time-fractional operator defined as [10]
⎧ ⎨
1
∂ β f (t ) = (n − β) ∂tβ ⎩ (n) f (t ),
∂t
f (n) (τ )
t 0
(t − τ )β −n+1
dτ ,
n − 1 < β < n,
β = n.
We here consider the following boundary conditions
T1 (0, t ) is finite,
(9)
T2 (a, t ) = T0 ,
(10)
T1 (R, t ) = T2 (R, t ),
(11)
k1
∂ T1 (R, t ) ∂ T2 (R, t ) = k2 , ∂r ∂r
(12)
and the initial conditions
Ti (r, 0) = T0 ,
∂ Ti (r, 0) = 0, i = 1, 2. ∂t
(13)
In order to obtain the solution of the above problem, we introduce a new dependent variable θ defined as
θ = r(T − T0 ).
(14)
Hence, the above problem can be transformed into a new problem in the rectangular coordinate system, that is, (7) and (8) can be rewritten in terms of θ , respectively, as
ρ1 c1
∂ 1+α θ1 ∂θ1 + τ1α ∂t ∂ t 1+α
= k1
∂ 2 θ1 + Qr, 0 ≤ r ≤ R, 0 < α ≤ 1, ∂ r2
(15)
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
ρ2 c2
∂ 1+α θ2 ∂θ2 + τ2α ∂t ∂ t 1+α
= k2
109
∂ 2 θ2 , R < r ≤ a. ∂ r2
(16)
The boundary and initial conditions are rewritten as
θ1 (0, t ) = 0,
(17)
θ2 (a, t ) = 0,
(18)
θ1 (R, t ) = θ2 (R, t ), ∂θ1 (R, t ) θ1 (R, t ) ∂θ2 (R, t ) θ2 (R, t ) − − k1 = k2 , ∂r R ∂r R ∂θi (r, 0) θi (r, 0) = 0, = 0, i = 1, 2. ∂t
(19) (20) (21)
Obviously, we observe that the temperature increase T as a function of time for different reduced distances r/R from the composite center is then obtained by
T = T − T0 =
θ r
.
(22)
3. Finite difference method for the direct problem In this section, we will adopt the finite difference method to solve the direct problem. The direct problem is to solve the temperature distribution θi (r, t )(i = 1, 2) from (15)–(21) where ρ 1 , c1 , ρ 2 , c2 , k1 , k2 , P, R, a, τ 1 ,τ 2 and α are already known. Ghazizadeh et al. [43] studied explicit and implicit finite difference schemes for the fractional Cattaneo equation. In this section, we will adopt the implicit finite difference schemes to solve the bi-layered composite medium problem. Firstly, we 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 = tn − tn−1 , n = 1, 2, . . . , N. Secondly, for 1 < β = 1 + α < 2, we employ the following discrete schemes
t=tn n
∂ β θ (r, t )
( t )−β = ω θ n−k+1 − 2θi,n−k + θi,n−k−1 , j j (3 − β) k=1 k i, j ∂ t β r=r
(23)
t=tn + θi,n−2 3θi,nj − 4θi,n−1 ∂θ (r, t )
j j , = ∂ t r=r 2 t
(24)
t=tn θi,nj+1 − 2θi,nj + θi,nj−1 ∂ 2 θ (r, t )
= , 2
∂r ( ri )2 r=r
(25)
i, j
i, j
i, j
where ωk = k2−β − (k − 1)2−β , k = 1, 2, . . .. Hence, we obtain the implicit finite difference schemes for (15) and (16):
ρ1 c1
n−1 n−2 n 3θ1, − 4θ1, + θ1, j j j
2 t
ρ2 c2
3θ
n 2, j
− 4θ
n−1 2, j
+ τα 1
+θ
n−2 2, j
2 t
+ τ2α
k+1 k+1 n θ1,k+1 − 2θ1, + θ1, n−k+1
( t )−β j+1 j j−1 n−k n−k−1 ωk θ1, j − 2θ1, j + θ1, j = k1 + Qr1, j , (3 − β) k=1 ( r1 )2 ( t )−β (3 − β)
n
θ
n−k n−k−1 ωk θ2,n−k+1 − 2θ2, + θ2, = k2 j j j
k=1
(26) k+1 2, j+1
− 2θ
k+1 + 2, j r2 2
( )
θ
k+1 2, j−1
.
(27)
The boundary and initial conditions can be discretized as k+1 θ1,0 = 0,
(28) (29)
θ
k+1 2,M
= 0,
θ
k+1 1,M
k+1 = 2,0 , k+1 k+1 − 1,M−1 1,M
k1
θ
θ
θi,0j = 0,
θ r 1
−
θ
k+1 1,M
R
= k2
θ
−θ r 2
k+1 2,1
k+1 2,0
−
θ
k+1 2,0
R
,
θi,0j − θi,−1j = 0, i = 1, 2. t
By solving (26)–(32), we can easily obtain the solution of the direct problem.
(30) (31)
(32)
110
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
4. Inverse analysis For the inverse problem, the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) are investigated. The measured temperature increase T as a function of time for different reduced distances r/R from the composite center are available [40]. The solution of the inverse problem is obtained through the minimization of the least squares norm in the following form
S(P) = ( T − Tc )T ( T − Tc ),
(33)
where T is the measured temperature increase vector in the experiment
T = [ T (t1,1 ), T (t1,2 ), . . . , T (t1,n ), T (t2,1 ), T (t2,2 ), . . . , T (t2,n )]T ,
(34)
Tc is the calculated temperature increase vector from the direct problem
Tc = [ Tc (t1,1 , P), Tc (t1,2 , P), . . . , Tc (t1,n , P), Tc (t2,1 , P), Tc (t2,2 , P), . . . , Tc (t2,n , P)]T ,
(35)
ti, j denotes the measured time nodes, n is the total number of the measured data from the experiment, and P = (α , τ1 , τ2 )T is regarded as being unknown. 4.1. Levenberg–Marquardt method For estimating the unknown parameter P = (α , τ1 , τ2 )T , Levenberg–Marquardt iterative method is introduced to minimize the above least squares norm [44]. The iterative algorithm for finding the unknown parameter is
Pk+1 = Pk + [(Jk )T Jk + μk k ]−1 (Jk )T ( T − Tkc ), where
Jk
(36)
is the fractional sensitivity matrix (Jacobian matrix in the fractional sense),
k is taken as
k = diag((Jk )T Jk ),
μk
is a positive scalar (damping parameter), (37)
and k is the number of iterations. The elements of the fractional sensitivity matrix are the first derivative of the estimated temperature increases at the measured time nodes ti, j with respect to the unknown parameter Pl (l = 1, 2, 3), that is, J = [Jα , Jτ1 , Jτ2 ]. The iteration begins with an appropriate initial guess P0 = (α 0 , τ10 , τ20 )T for the unknown parameter P = (α , τ1 , τ2 )T . The stopping criterion of the iterative algorithm is
k+1
P − Pk < ,
(38)
where is a chosen tolerance. It is worthy to mention that the three parameters α , τ 1 , τ 2 should vary in their physical range in the iteration process, that is, 0 < α ≤ 1, τ 1 ≥ 0, τ 2 ≥ 0. 4.2. Calculation of the fractional sensitivity matrix 4.2.1. Calculation of Jα The Jacobian equations with respect to α are obtained through differentiating (15) and (16) with respect to the unknown parameter α , that is,
ρi ci
∂ 1+α θi ∂ ∂ 1+α θi ∂ 2 Jαi ∂ Jαi α α + τi ln τi 1+α + τi = k , i = 1, 2, i ∂t ∂α ∂ t 1+α ∂t ∂ r2
(39)
∂θ (r,t )
i where Jα i (r, t ) = ∂α , i = 1, 2. The boundary and initial conditions are
Jα 1 (0, t ) = 0,
(40)
Jα 2 (a, t ) = 0,
(41)
Jα 1 (R, t ) = Jα 2 (R, t ),
k1
∂ Jα1 (R, t ) Jα1 (R, t ) − ∂r R
Jα i (r, 0) = 0,
(42)
= k2
∂ Jα2 (R, t ) Jα2 (R, t ) − , ∂r R
∂ Jαi (r, 0) = 0, i = 1, 2. ∂t
(43)
(44)
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
111
∂ ∂ 1+α θi (r,t ) , we have Concerning the discrete form of the term ∂α ∂ t 1+α
t=tn n n−k+1
( t )−β ∂ ∂ 1+α θi (r, t )
∂ = ω θ − 2θi,n−k + θi,n−k−1 j j
∂α ∂α (3 − β) k=1 k i, j ∂ t 1+α r=ri, j n n
n−k+1
( t )−β ∂ωk n−k+1 n−k n−k−1 n−k n−k−1 = − 2θi, j + θi, j + ωk Jα,i, j − 2Jα,i, j + Jα,i, j θ (3 − β) k=1 ∂α i, j k=1 n
∂ ( t )−β + ωk θi,n−k+1 − 2θi,n−k + θi,n−k−1 , j j j ∂α (3 − β) k=1
(45)
∂ω ∂ω ∂ ( t )−β ) = ( t )−β ( − ln t + where β = 1 + α , d1 = ∂α1 = 0, dk = ∂αk = −k2−β ln k + (k − 1)2−β ln (k − 1), k = 2, 3, . . . , ( ∂α (3−β) (3−β) −∂
∂α ( (3−β)) ) (3−β)
t )−β d ln (x) = ( is the digamma function. dx (3−β) ( − ln t + ψ(3 − β)), and ψ(x) = Solving (39)–(45) by implicit finite difference method similarly as in Section 3, we can obtain Jα .
4.2.2. Calculation of Jτ1 The Jacobian equations with respect to τ 1 are obtained through differentiating (15) and (16) with respect to the unknown parameter τ 1 :
ρ1 c1 ρ2 c2 that is,
ρ1 c1 ρ2 c2
1+α ∂ 1+α Jτ1 1 ∂ Jτ1 1 θ1 α −1 ∂ + τ1α + ατ 1 ∂t ∂ t 1+α ∂ t 1+α
∂ 1+α Jτ1 2 ∂ Jτ1 2 + τ2α ∂t ∂ t 1+α ∂ 1+α Jτ1 1 ∂ Jτ1 1 + τ1α ∂t ∂ t 1+α ∂ 1+α Jτ1 2 ∂ Jτ1 2 + τ2α ∂t ∂ t 1+α
where Jτ1 i (r, t ) =
= k1
∂ 2 Jτ1 1 , 0 ≤ r ≤ R, ∂ r2
(46)
= k2
∂ 2 Jτ1 2 , R < r ≤ a, ∂ r2
(47)
= k1
∂ 2 Jτ1 1 ∂ 1+α θ1 − ατ1α −1 , 0 ≤ r ≤ R, 2 ∂r ∂ t 1+α
(48)
= k2
∂ 2 Jτ1 2 , R < r ≤ a, ∂ r2
(49)
∂θi (r,t ) ∂τ1 , i = 1, 2.
The boundary and initial conditions are
Jτ1 1 (0, t ) = 0,
(50)
Jτ1 2 (a, t ) = 0,
(51)
Jτ1 1 (R, t ) = Jτ1 2 (R, t ),
(52)
∂ Jτ1 1 (R, t ) Jτ1 1 (R, t ) ∂ Jτ1 2 (R, t ) Jτ1 2 (R, t ) = k2 , − − ∂r R ∂r R ∂ Jτ1 i (r, 0) = 0, i = 1, 2. Jτ1 i (r, 0) = 0, ∂t k1
(53) (54)
Solving (48)–(54) by implicit finite difference method similarly as in Section 3, we can obtain Jτ1 . 4.2.3. Calculation of Jτ2 The Jacobian equations with respect to τ 2 are obtained through differentiating (15) and (16) with respect to the unknown parameter τ 2 :
ρ1 c1 ρ2 c2
∂ 1+α Jτ2 1 ∂ Jτ2 1 + τ1α ∂t ∂ t 1+α
= k1
∂ 2 Jτ2 1 , 0 ≤ r ≤ R, ∂ r2
∂ 1+α Jτ2 2 ∂ Jτ2 2 ∂ 1+α θ2 + τ2α + ατ2α −1 ∂t ∂ t 1+α ∂ t 1+α
= k2
∂ 2 Jτ2 2 , R < r ≤ a, ∂ r2
(55)
(56)
112
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
that is,
ρ1 c1 ρ2 c2
∂ 1+α Jτ2 1 ∂ Jτ2 1 + τ1α ∂t ∂ t 1+α ∂ 1+α Jτ2 2 ∂ Jτ2 2 + τ2α ∂t ∂ t 1+α
where Jτ2 i (r, t ) =
= k1
∂ 2 Jτ2 1 , 0 ≤ r ≤ R, ∂ r2
(57)
= k2
∂ 2 Jτ2 2 ∂ 1+α θ2 − ατ2α −1 , R < r ≤ a, 2 ∂r ∂ t 1+α
(58)
∂θi (r,t ) ∂τ2 , i = 1, 2.
The boundary and initial conditions are
Jτ2 1 (0, t ) = 0,
(59)
Jτ2 2 (a, t ) = 0,
(60)
Jτ2 1 (R, t ) = Jτ1 2 (R, t ),
(61)
k1
∂ Jτ2 1 (R, t ) Jτ2 1 (R, t ) − ∂r R
Jτ2 i (r, 0) = 0,
= k2
∂ Jτ2 2 (R, t ) Jτ2 2 (R, t ) − , ∂r R
∂ Jτ2 i (r, 0) = 0, i = 1, 2. ∂t
(62)
(63)
Solving (57)–(63) by implicit finite difference method similarly as in Section 3, we can obtain Jτ2 . 4.3. Computational algorithms Suppose that the measured temperature increase vector T in the experiment is given, the initial guess P0 = (α 0 , τ10 , τ20 )T is available for the unknown parameter P = (α , τ1 , τ2 )T , choose a value for μ0 , for example, μ0 = 10, and set k = 0, Step 1. Solve the direct problem (26)–(32) with the available estimate Pk in order to obtain the temperature increase vector
Tkc = [ Tc (t1,1 , Pk ), Tc (t1,2 , Pk ), . . . , Tc (t1,n , Pk ), Tc (t2,1 , Pk ), Tc (t2,2 , Pk ), . . . , Tc (t2,n , Pk )]T ; Step 2. Step 3. Step 4. Step 5.
(64)
S(Pk )
Compute ; Compute the fractional sensitivity matrix Jk , and compute k ; Compute Pk = [(Jk )T Jk + μk k ]−1 (Jk )T ( T − Tkc ), and compute the new estimate Pk+1 = Pk + Pk ; , and Solve the direct problem with the new estimate Pk+1 in order to obtain the new temperature increase vector Tk+1 c compute S(Pk+1 ); Step 6. If S(Pk+1 ) ≥ S(Pk ), replace μk by 10μk , and return to Step 4; Step 7. If S(Pk+1 ) < S(Pk ), accept the new estimate Pk+1 , and replace μk by μk /10; Step 8. Check the stopping criteria. If satisfied, stop the iterative; otherwise, replace k by k + 1, and return to Step 3. 5. Results and discussions Andrä et al. [40] simulated the thermal behavior of localized magnetic hyperthermia applied to a breast carcinoma. They 3 considered a small sphere tumor of radius R = 0.00315 m with a constant power density Q = 6150, 000 W/m embedded in ex3 3 tended muscle tissue. The thermal parameters were taken to be ρ1 = 1, 660, 000 g/m , ρ2 = 1, 000, 000 g/m , c1 = 2.54 J/(g K), c2 = 3.72 J/(g K), k1 = 0.778 W/(K m), k2 = 0.642 W/(K m). 5.1. Optimal estimation of the parameter P = (α , τ1 , τ2 )T The measured temperature increase values which are used by us are extracted from the symbols at reduced distances r/R = 0.65, 1.10, the error of T and t are within the extent of the used symbols. For parameter a, we here choose a = 11R. Chosen the tolerance = 10−4 and an initial guess P0 = (0.8, 10, 10)T , the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) are obtained by the Levenberg–Marquardt method. The computational results are listed in Table 1. From Table 1, we can see that the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) are obtained after 12 iterations. Fig. 1 demonstrates the comparison between the calculated temperature increase values with the optimal estimations of the fractional derivative and the relaxation time parameters and the measured temperature increase values in the experiment at
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
113
Table 1 Optimal estimation for parameter P = (α , τ1 , τ2 )T . Initial guess P0
Number of iterations
Estimated value P
(0.8, 10, 10)T
12
(0.4381, 0.5717, 1.7914)T
35 Experimental data [40] Present result
Temperature Increase
T
30
=0.4381 =0.5717 1
25
r/R=0.65
=1.7914
2
20 r/R=1.10 15 r/R=1.39 10
r/R=1.68 r/R=1.98
5
0
0
50
100
150 t
200
250
300
Fig. 1. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10, 1.39, 1.68, 1.98.
Table 2 Optimal estimation for different initial guess. Initial guess P0
Number of iterations
Estimated value P
(0.9, 10, 10)T (0.4, 10, 10)T (0.8, 1, 10)T (0.8, 20, 10)T (0.8, 10, 1)T (0.8, 10, 20)T
12 10 13 10 11 12
(0.438139, 0.571676, 1.791362)T (0.438141, 0.571681, 1.791409)T (0.438137, 0.571670, 1.791315)T (0.438141, 0.571680, 1.791403)T (0.438137, 0.571670, 1.791316)T (0.438140, 0.571677, 1.791370)T
the reduced distances r/R = 0.65, 1.10, 1.39, 1.68, 1.98. From Fig. 1, we can obviously observe that the calculated temperature increase values are in excellently consistent with the measured temperature increase values in the experiment when α , τ 1 , τ 2 are chosen as the optimal estimated values in Table 1, which is much more better than the result of [40]. This fact demonstrates that the proposed fractional thermal wave model is efficient and accurate in modeling the heat transfer in the hyperthermia experiment, and the proposed numerical method is effective for solving the corresponding the fractional inverse thermal wave problem in a composite medium.
5.2. Effect of the initial guess Table 2 demonstrates the optimal estimations of the fractional derivative α and the relaxation time parameters τ 1 and τ 2 as well as the number of iterations for different initial guess P0 = (α 0 , τ10 , τ20 )T . From Table 2, we can obviously find out that different initial guess P0 has little impact on the final optimal estimation. This also demonstrates that Levenberg–Marquardt method based on fractional derivative is efficient for solving the fractional inverse thermal wave problem in a composite medium.
114
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
35 r/R=0.65 =0.5717, 2=1.7914
30
1
Temperature Increase
T
=0.4381 =0.1
25
=0.9
20 r/R=1.10 15
10 Experimental data [40] Present result( =0.4381) Present result( =0.1) Present result( =0.9)
5
0
0
50
100
150 t
200
250
300
Fig. 2. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10 for different α .
5.3. Effect of the fractional derivative α Fractional derivative α is a very important parameter in determining the temperature distribution. Fig. 2 demonstrates the comparison between the measured temperature increase values in the experiment and the calculated temperature increase values at r/R = 0.65, 1.10 when α is chosen as the optimal estimated value 0.4381 in Table 1 and another values 0.1 and 0.9, τ 1 , τ 2 are chosen as the optimal estimated values in Table 1. From Fig. 2, we can obviously find out that the curves are higher when α is larger, and lower when α is smaller, and the curves protrude to upper left slower when α is becoming larger. Additionally, we can observe that the calculated temperature increase values agree much better with the measured temperature increase values in the experiment when α is chosen as the optimal estimated value 0.4381 than α is chosen as 0.1 and 0.9. 5.4. Effect of the relaxation time parameters τ 1 and τ 2 Relaxation time parameters τ 1 and τ 2 are important parameters in determining the temperature distribution, too. In this section, we will discuss the effect of the relaxation time parameters τ 1 and τ 2 for the temperature increase T. Fig. 3 demonstrates the comparison between the measured temperature increase values in the experiment and the calculated temperature increase values at r/R = 0.65, 1.10 when α , τ 1 , τ 2 are chosen as the optimal estimated values in Table 1 and another value τ1 = 10. From Fig. 3, we can obviously find out that the curves on the lower left are higher when τ 1 is smaller, and lower when τ 1 is larger. But it is just contrary on the top right. Additionally, we can observe that the calculated temperature increase values agree much better with the measured temperature increase values in the experiment when τ 1 is chosen as the optimal estimated value 0.5717 than τ 1 is chosen as 10. Fig. 4 demonstrates the comparison between the measured temperature increase values in the experiment and the calculated temperature increase values at r/R = 0.65, 1.10 when α , τ 1 , τ 2 are chosen as the optimal estimated values in Table 1 and another values τ2 = 10. From Fig. 4, we can obviously find out that the curves are higher when τ 2 is smaller, and lower when τ 2 is larger. Additionally, we can observe that the calculated temperature increase values agree much better with the measured temperature increase values in the experiment when τ 2 is chosen as the optimal estimated value 1.7914 than τ 2 is chosen as 10. 5.5. Sensitivity problem The sensitivity coefficients JPi (i = 1, 2, 3) are the measurements of the sensitivity of the estimated temperature increase T with respect to the changes in the parameter Pi at the reduced distances r/R = 0.65, 1.10. Here we introduce the relative
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
115
35 =0.4381, 2=1.7914
30
r/R=0.65
Temperature Increase
T
=0.5717
1
25
20 r/R=1.10 15
10 Experimental data [40] Present result( 1=0.5717)
5
=10
Present result( 1=10)
1
0
0
50
100
150 t
200
250
300
Fig. 3. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10 for different τ 1 .
35 =0.4381, 1=0.5717
30
r/R=0.65
=1.7914
Temperature Increase
T
2
25
20 r/R=1.10 15
10 Experimental data [40] =1.7914)
=10
2
5
2
Present result( 2=10) 0
0
50
100
150 t
200
250
300
Fig. 4. Comparison between the experimental data and the calculated temperature increase values at r/R = 0.65, 1.10 for different τ 2 .
116
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
0.06 J J1 J
r/R=0.65
0.05
2
Sensitivity Coefficients
0.04 J
0.03 0.02 J
0.01
2
J
1
0 −0.01 −0.02
0
50
100
150 t
200
250
300
Fig. 5. Sensitivity coefficients for α , τ 1 , τ 2 at r/R = 0.65.
0.05 J r/R=1.10
J1 J2
Sensitivity Coefficients
0.04
0.03 J 0.02
0.01 J2
J1
0
−0.01
0
50
100
150 t
200
250
300
Fig. 6. Sensitivity coefficients for α , τ 1 , τ 2 at r/R = 1.10.
sensitivity coefficients
JPi =
P i ∂ T , i = 1, 2, 3. T0 ∂ Pi
(65)
Figs. 5 and 6 show the sensitivity coefficients Jα , Jτ1 and Jτ2 at r/R = 0.65, 1.10, respectively. The sensitivity coefficients are linear independency, it is possible to have a correct estimation for the unknown parameters.
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
117
6. Conclusions In this paper, a fractional thermal wave model is formulated for a bi-layered spherical tissue. The solution of the direct problem is obtained by implicit finite difference method. For the inverse problem, we obtain the optimal estimations of the fractional derivative α and the relaxation time parameters τi (i = 1, 2) simultaneously by means of the Levenberg–Marquardt method. Compared with the experimental data, we can obviously find out that the calculated temperature increase values are in excellently consistent with the measured temperature increase values in the experiment, which demonstrates that the proposed fractional thermal wave model is efficient and accurate in modeling the heat transfer in the hyperthermia experiment, and Levenberg– Marquardt method for solving the fractional inverse thermal wave problem in a spherical composite medium is effective and valid. From the viewpoint of methodology, this paper provides an efficient numerical method for simultaneously estimating multiple parameters in fractional thermal wave model in a spherical composite medium. Acknowledgment The project was supported by the National Natural Science Foundation of China (11472161, 11102102 and 91130017). References [1] P.T. Hsu, Y.H. Chu, An inverse non-Fourier heat conduction problem approach for estimating the boundary condition in electronic device, Appl. Math. Model. 28 (2004) 639–652. [2] V. Vishwakarma, A.K. Das, P.K. Das, Analysis of non-Fourier heat conduction using smoothed particle hydrodynamics, Appl. Therm. Eng. 31 (2011) 2963– 2970. [3] C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, C.R. Acad. Sci. 247 (1958) 431–433. [4] P. Vernotte, Paradoxes in the continuous theory of the heat equation , C.R. Acad. Sci. 246 (1958) 154-3. [5] L. Wang, F. Fan, Modeling bioheat transport at macroscale, J. Heat Transf. 133 (2011) 011010. [6] W. Kaminski, Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure, J. Heat Transf. 112 (1990) 555–560. [7] K. Mitra, S. Kumar, A. Vedevarz, M.K. Moallemi, Experimental evidence of hyperbolic heat conduction in processed meat, J. Heat Transf. 117 (1995) 568–573. [8] W. Roetzel, N. Putra, S.K. Das, Experiment and analysis for non-Fourier conduction in materials with non-homogeneous inner structure, Int. J. Therm. Sci. 42 (2003) 541–552. [9] P.J. Antaki, New interpretation of non-Fourier heat conduction in processed meat, J. Heat Transf. 127 (2005) 189–193. [10] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [11] F. Liu, C. Yang, K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math. 231 (2009) 160–176. [12] X.Y. Jiang, M.Y. Xu, H.T. Qi, The fractional diffusion model with an absorption term and modified Fick’s law for non-local transport processes , Nonlinear Anal. Real World Appl. 11 (2010) 262–269. [13] C. Li, Z. Zhao, Y. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62 (2011) 855–875. [14] L. Zheng, F. Zhao, X. Zhang, Exact solutions for generalized Maxwell fluid flow due to oscillatory and constantly accelerating plate, Nonlinear Anal. Real World Appl. 11 (2010) 3744–3751. [15] B. Jin, R. Lazarov, Z. Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations, SIAM J. Numer. Anal. 51 (2013) 445–466. [16] B. Jin, R. Lazarov, J. Pasciak, Z. Zhou, Error analysis of a finite element method for the space-fractional parabolic equation, SIAM J. Numer. Anal. 52 (2014) 2272–2294. [17] F. Zeng, C. Li, F. Liu, I. Turner, The use of finite difference/element approaches for solving the time-fractional subdiffusion equation, SIAM J. Sci. Comput. 35 (2013) A2976–A3000. [18] F. Zeng, F. Liu, C. Li, K. Burrage, I. Turner, V. Anh, A Crank–Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction– diffusion equation, SIAM J. Numer. Anal. 52 (2014) 2599–2622. [19] B. Yu, X.Y. Jiang, H. Xu, A novel compact numerical method for solving the two-dimensional non-linear fractional reaction–subdiffusion equation 68 (2015) 923–950. [20] A. Eke, P. Herman, L. Kocsis, L.R. Kozak, Fractal characterization of complexity in temporal physiological signals, Physiol. Meas. 23 (2002) R1. [21] 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. Phys. A Math. Gen. 37 (2004) R161. [22] R.L. Magin, Fractional calculus models of complex dynamics in biological tissues, Comput. Math. Appl. 59 (2010) 1586–1593. [23] W. Tan, C. Fu, C. Fu, W. Xie, H. Cheng, An anomalous subdiffusion model for calcium spark in cardiac myocytes, Appl. Phys. Lett. 91 (2007) 3901. [24] K. Li, C. Fu, H. Cheng, W. Tan, Anomalous subdiffusion of calcium spark in cardiac myocytes, Cell. Mol. Bioeng. 4 (2011) 457–465. [25] X.Y. Jiang, H.T. Qi, Thermal wave model of bioheat transfer with modified Riemann–Liouville fractional derivative, J. Phys. A. Math. Theor. 45 (2012) 485101. [26] 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. [27] 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. [28] X.Y. Jiang, S.Z. Chen, Analytical and numerical solutions of time fractional anomalous thermal diffusion equation in composite medium, ZAMM Z. Angew. Math. Mech. 95 (2015) 156–164. [29] G. Fourestey, M. Moubachir, Solving inverse problems involving the Navier–Stokes equations discretized by a Lagrange–Galerkin method, Comput. Methods Appl. Mech. Eng. 194 (2005) 877–906. [30] A.G. Fatullayev, E. Can, N. Gasilov, Comparing numerical methods for inverse coefficient problem in parabolic equation, Appl. Math. Comput. 179 (2006) 567–571. [31] M. Ebrahimian, R. Pourgholi, M. Emamjome, P. Reihani, A numerical solution of an inverse parabolic problem with unknown boundary conditions, Appl. Math. Comput. 189 (2007) 228–234. [32] A. Shidfar, J. Damirchi, P. Reihani, An stable numerical algorithm for identifying the solution of an inverse problem, Appl. Math. Comput. 190 (2007) 231–236. [33] U. Albocher, A.A. Oberai, P.E. Barbone, I. Harari, Adjoint-weighted equation for inverse problems of incompressible plane-stress elasticity, Comput. Methods Appl. Mech. Eng. 198 (2009) 2412–2420. [34] G. Deolmi, F. Marcuzzi, A parabolic inverse convection–diffusion–reaction problem solved using space-time localization and adaptivity, Appl. Math. Comput. 219 (2013) 8435–8454. [35] Y.F. Wang, K.G. Reyes, K.A. Brown, C.A. Mirkin, W.B. Powell, Nested-batch-mode learning and stochastic optimization with an application to sequential multistage testing in materials science, SIAM J. Sci. Comput. 37 (3) (2015) B361–B381.
118
B. Yu et al. / Applied Mathematics and Computation 274 (2016) 106–118
[36] D.A. Murio, Time fractional IHCP with Caputo fractional derivatives, Comput. Math. Appl. 56 (2008) 2371–2381. [37] J. Cheng, J. Nakagawa, M. Yamamoto, T. Yamazaki, Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation, Inverse Probl. 25 (2009) 115002. [38] W. Tian, C. Li, W. Deng, Y. Wu, Regularization methods for unknown source in space fractional diffusion equation, Math. Comput. Simul. 85 (2012) 45–56. [39] X. B.Yu, 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 (2) (2015) 153–161. [40] 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. [41] K.C. Liu, C.T. Lin, Solution of an inverse heat conduction problem in a bi-layered spherical tissue, Numer. Heat Transf. Part A Appl. 58 (2010) 802–818. [42] Z.M. Odibat, N.T. Shawagfeh, Generalized Taylor’s formula, Appl. Math. Comput. 186 (2007) 286–293. [43] H.R. Ghazizadeh, M. Maerefat, A. Azimi, Explicit and implicit finite difference schemes for fractional Cattaneo equation, J. Comput. Phys. 229 (2010) 7042– 7057. [44] M.N. Özisik, Inverse Heat Transfer: Fundamentals and Applications, CRC Press, 2000.