Fractional Order Dual-Phase-Lag Model of Heat Conduction in a Composite Spherical Medium

In the paper, a solution of the fractional dual-phase-lag heat conduction problem is presented. The considerations are related to the heat conduction in a multi-layered spherical medium with azimuthal symmetry. The final form of the analytical solution is given in a form of the double series of spherical Bessel functions and Legendre functions. Numerical calculations concern the study of the effect of the order of the Caputo derivative on the temperature distribution in a composite solid sphere, hemisphere and spherical cone.


Introduction
Mathematical modeling of the heat conduction is an important stage in the design of systems subjected to thermal loads, because it allows the appropriate selection of materials to avoid adverse phenomena related to the occurrence of thermal stresses. For example, the temperature distribution in the tested system may lead to the formation of thermal stresses, which may cause micro-cracks in the components of this system [1][2][3], and cyclical changes in the thermal load of the device may cause undesirable vibrations of the components of this device [4][5][6]. On the other hand, the heat source causing high temperature of biological tissue may destroy diseased tissue, but it may also damage healthy tissue [7,8]. Different mathematical models are used to accurately describe the heat conduction in the considered bodies.
The classical mathematical model of heat conduction is derived from Fourier's law of the heat conduction. The Fourier law establishes a proportionality between the heat flux vector and the temperature gradient [9] q(r, t) = −λ ∇T(r, t), where q is the heat flux vector, λ is the thermal conductivity of the material, r is the point in the considered region, t is the time, ∇ is the gradient operator and T is the temperature. Although the Fourier law quite accurately describes the heat conduction in most practical macroscopic problems, the relationship (1) implies an unrealistic infinite speed of the heat propagation; this means that the sudden temperature change at some point in the domain will be felt everywhere and instantaneously at distant points in the domain-hence, Fourier's law can be treated as having the unphysical property [10]. To eliminate this drawback of the mathematical model, the Fourier law (1) is replaced by the following relationship [11] q r, t + τ q = −λ∇T(r, t + τ T ), (2) − ∇ · q(r, t) + g(r, t) = ρC p ∂T(r, t) ∂t , where g(r, t) is the volumetric rate of the heat generation, ρ is the density of the material and C p is the specific heat capacity. Eliminating the heat flux from Equations (3) and (4), we obtain the heat conduction equation in the form where ∇ 2 is the Laplace operator and κ = λ ρC p is the thermal diffusivity. The dual-phase-lag heat shown in Equation (5) has been applied in mathematical modeling of the heat transfer in functionally graded materials [13][14][15], ultrafast pulselaser heating problems [16,17], porous media [18][19][20], nanocomposites [21,22], and living tissue [23][24][25]. If τ q = τ T = 0 in Equation (5), then the classical parabolic heat conduction equation is obtained. For τ T = 0 and τ q > 0, one obtains the single-phase-lag equation of hyperbolic type. The wave character of this heat conduction equation was used in the investigations of propagations of heat waves in papers [26,27].
A generalization of the dual-phase-lag model of the heat conduction is obtained by replacement of the derivatives in the constitutive equation and energy equation by the derivatives of non-integer order. In the present paper, the Caputo derivative of non-integer order α is used and is defined by [28] The properties of the derivative can be found in many books devoted to fractional calculus, for instance in the books [28][29][30][31].
Replacing the time derivatives in the constitutive Equation (3) and energy Equation (4) by the Caputo time fractional derivatives, one obtains − ∇ · q(r, t) + g(r, t) = ρC p ν α−1 ∂ α T(r, t) where the coefficient ν is introduced to keep the accordance of dimensions. Combining Equations (7) and (8), one obtains the heat conduction equation with the Caputo fractional derivatives in the form We use this equation for modeling the heat conduction in a multi-layered body. According to the authors' knowledge, this approach to the problem of two-dimensional time-fractional heat conduction in a multi-layered medium was not used in the literature. Due to the definition of the Caputo derivative, we can note that the differential Equation (9) describes the temperature distribution in a body, taking into account the history of changes of the temperature.
The fractional Equation (9) describes the heat conduction in a region which is specified by the considered medium. In this paper, we consider the fractional heat conduction in a spherical body. The Laplace operator which occurs in Equation (9) in the spherical co-ordinates is defined as follows where r is the radial coordinate, ϕ is the polar and θ is the azimuthal angle coordinate.
Assuming that the temperature distribution in the body is azimuthally invariant (azimuthal symmetry), the last term in the bracket in Equation (10) can be dropped. A simplified form of the Laplace operator can be obtained by introducing a new variable µ which is related to the polar angle ϕ by the relationship Taking into account the azimuthal symmetry, the Laplace operator in the new coordinates can be written in the form In this paper, we present a dual-phase-lag fractional heat conduction model for a multilayered spherical body with azimuthal symmetry. An analytical solution of the problem is derived in the form of the double series of spherical Bessel functions and Legendre functions. Numerical computations of the temperature distribution in the body include a composite solid sphere, hemisphere and spherical cone. Effects of the fractional order of time-derivatives and phase lagging occurring in the heat conduction model on the temperature distribution in the considered bodies are investigated numerically.
The problem considered in this paper is a continuation of the research that has been presented in papers [32][33][34][35][36]. A solution of the fractional heat conduction problem in the solid sphere was presented in [32]. The authors studied a fractional single-phase lag heat conduction problem in the whole-space 1D domain [34], in the slab [35] and in the hollow cylinder [33]. A fractional dual-phase-lagging heat conduction in the whole-space 1D domain was presented in [36].

Formulation of the Problem
Let us consider the fractional heat conduction in a spherical body (with coordinates r ≡ {r, ϕ} consisting of n concentric spherical layers which are defined by: r i−1 r r i (i = 1, . . . , n), 0 ϕ φ, where 0 <φ π. The body is a full solid sphere whenφ tends to π, forφ = π/2, the body is a hemisphere, and for 0 <φ < π/2, the body is a spherical cone (Figure 1). The heat conduction in the i-th spherical layer is governed by the following timefractional differential equation where T i ≡ T i (r, µ, t) is the temperature, λ i is the thermal conductivity, κ i is the thermal diffusivity, g i ≡ g i (r, µ, t) is the volumetric heat source in the i-th layer and α ∈ (0, 1] denotes the fractional order of the Caputo derivative, µ ∈ µ ϕ , 1 and µ ϕ = cos(φ).
The fractional differential Equation (13) for i = 1, . . . , n is complemented by boundary and initial conditions and the conditions providing the perfect thermal contact of the neighboring layers. The conditions are as follows: ∂T i (r, µ, t) ∂µ where a ∞ and T ∞ are the outer heat transfer coefficient and ambient temperature, respectively. It should be noted that in general, the boundary conditions (17) and (18) should be formulated using the dual-phase-lag model [37,38]. However, if values of the relaxation times are identical in two neighboring layers, than we can assume the simplification in the above boundary conditions. The initial conditions are assumed in the form

Solution of the Problem
An analytical solution to the initial-boundary problem (13)- (20) can be presented in the form of a sum where the function φ i satisfies the fractional differential Equation (13) and homogenous boundary conditions (transient problem) and ψ i is a solution of a steady-state problem.
Substituting the function T i in the form (21) into Equation (13), we receive differential equations for the functions φ i and ψ i . For the function φ i , one obtains an equation with fractional time derivative and for the function ψ i , one obtains the Laplace equation Suppose that φ i is a function of the form Taking into account Equation (24) in the homogenous differential equation obtained by omitting the last term in Equation (22) and separating the space and time variables, we receive the Helmholtz equation for the function Λ i where Substituting the functions ψ i and Λ i into Equations (23) and (25), separating the variables and assuming the separation constant as β(β + 1) where β is a real number, one obtains the three homogenous differential equations: Lagrange equation, Euler equation and spherical Bessel equation: The functions M, R Λ i , and R ψ i satisfy boundary conditions which are obtained by taking the functions (21), (24) and (26) in conditions (14)- (18). The function M satisfies the conditions: The functions R Λ i (r), R ψ i (r) satisfy the same conditions at r = 0 and at interfaces r = r i , i = 1, 2, . . . , n − 1 (the superscripts are omitted): where ϑ i = λ i+1 λ i . Moreover, the function R Λ n (r) satisfies the homogenous condition at and the function R ψ n (r), as the radial part of the function ψ n (r, The solutions of Equations (27)- (29) with appropriate conditions among (30)- (35) are presented in Section 3.1.

Solution of the Lagrange Equation
The solution to Equation (27), which satisfies the condition (30), is the function where P β (µ) is the Legendre function of the first kind and A 1 is a constant. Using the derivative of the Legendre function [39] dP β (µ) and the boundary condition (31), one obtains the following equation The roots of this equation for µ ϕ = −1 (solid sphere) are β m = m, m = 0, 1, 2, . . . and for µ ϕ = 0 (hemisphere) are β m = 2m, m = 0, 1, 2, . . .. In these cases, the eigenfunctions P m (µ) and P 2m (µ) where m is a positive integer number are the Legendre polynomials. The roots to Equation (38) for µ ϕ ∈ (−1, 0) ∪ (0, 1) are numerically determined. The functions P β m (µ), m = 0, 1, 2, . . ., create an orthogonal set of functions. The orthogonality condition of the functions can be derived using an indefinite integral of the product of functions P β (µ) and P β (µ) where β = β . Utilizing Equation (27), this integral can be expressed as follows where C is an arbitrary constant. Taking into account the antiderivative given by Equation (39) and the boundary condition (31), one obtains In order to find the antiderivative of the square of the Legendre function, the integral occurring in Equation (39), employing Equation (37), is rewritten in the form Evaluating limits of the expressions occurring on the left and right-hand sides of Equation (41) as β tends to β, the square of the norm N µ m of the Legendre function for m = 1, 2, . . ., one obtains in the form where β m are the roots of the Equation (38). For m = 0, we have N µ 0 = 1 − µ ϕ . Finally, the orthogonality condition of the functions P β m (µ), m = 0, 1, 2, . . ., can be written as where δ m n is the Kronecker delta.
The system of Equations (45) and (46) is complemented by an equation which is obtained on the basis of the boundary condition (35) for the function ψ n . The functions ψ i are given by Equation (26) as a product of two functions. To fulfil the condition (35), we assume that Substituting the function ψ n into the condition (35), multiplying the received equation by P β m (µ), integrating with respect to µ in the interval µ ϕ , 1 and using the orthogonality condition (43), one obtains the condition for the function R Employing Equation (44) for i = n in Equation (48), the following condition is obtained where the integral I m may be expressed in an analytical form as T a for m = 1, 2, . . ..

Solution of the Spherical Bessel Equation in the Multi-Layered Region
The general solution to Equation (29) can be written as follows where i are constants, and j β m and y β m are spherical Bessel functions of the first and second kind, respectively. The spherical Bessel functions are defined by (see Ref. [9]) where J β and Y β are the Bessel functions of the first and second kind, respectively. The functions R Λ i,m (r) are defined for r ∈ [r i−1 , r i ], i = 1, 2, . . . , n, wherein r 0 = 0. As the function y β (z) tends to infinity when z tends to zero, i.e., lim z→0 y β (z) = −∞, then taking into account the condition (14), we assume C 2, 1 = 0 in Equation (51) for i = 1.
Substituting the function R Λ i,m into the continuity conditions (33) and the boundary condition (34), we obtain a system of 2n − 1 homogenous equations with unknown C 1, 1 , C 1,2 , C 2, 2 , . . ., C n, 1 , C n, 2 . This system of equations in the matrix form can be written as where C = [C 1,1 C 1,2 C 2,2 . . . C 1,n C 2,n ] T is the column matrix of the unknowns and A = a ij 1 i,j 2n−1 is the coefficients matrix of the equations system. The non-trivial solution of Equation (54) exists if the characteristic equation is satisfied Next, this equation is solved numerically with respect to ω for β = β m , m = 0, 1, 2 . . .. Next, for the roots ω j, m , j = 1, 2, . . . , m = 0, 1, 2, . . . of Equation (55), the coefficients C 1,2 , C 2,2 , . . . , C 1,n , C 2,n (index m is omitted) are successfully determined as a solution of Equation (54) with C 1,1 = 1. Taking into account the eigenvalues ω j, m and the coefficients C 1,i , C 2,i in Equation (51), we obtain the functions R Λ i, j, m . These functions satisfy the orthogonality condition which can be derived using the differential Equation (29), the conditions at the interfaces (33) and the boundary condition (34). The orthogonality condition has the form where The functions Λ i , defined by Equation (26), can be now written as

Solution of the Time-Fractional Differential Equation
The functions Λ i, j, m are used to express the functions φ i in the double series form where the functions θ j, m satisfy the nonhomogeneous Equation (22). Substituting the functions (59) into Equation (22) and using the orthogonality conditions (43) and (56), we receive the fractional differential equation for the functions θ j, m in the form This differential equation is complemented by initial conditions which are derived using Equations (19), (20) and (59) and the orthogonality conditions (43) and (56). The initial conditions for the function θ j, m have the form In order to determine a solution of the initial value problem (60) and (61), we apply the Laplace transform technique. The Laplace transformf (s) of the function f (t) is defined as The Laplace transform of the Caputo derivative of order α ∈ (0, 1] is given by The application of the Laplace transform to Equation (60) allows us to express the transform θ j,m (s) in the form For the purpose of determining the inverse Laplace transform L −1 θ j, m (s) , we find firstly four inverse transforms: j, m (s) . Introducing notation: where z ± j, m = p j, m ± ϑ j, m . Employing the Laplace transform pair given by where E γ α, β is a three parameter Mittag-Leffler function, which is also known as the Prabhakar function [28], and we obtain the inverse transform K j,m (t) = L −1 K j, m (s) for ϑ j, m = 0 in the form If ϑ j, m = 0, then on the basis of (67), we have Similarly, we receive Using the property of the Laplace transforms, we can find the inverse Laplace transform L −1 s −1K j, m (s) as an integral of the function K j,m (t). Utilizing the formula for integration of the Mittag-Leffler and Prabhakar functions in Equations (69) and (70), we obtain The inverse Laplace transform L −1 s α−1K j, m (s) can be obtained using Formula (68) Employing the functions K j,m , K * j,m , K * * j,m and the Laplace transformθ j, m given by Equation (65), we can write the functions θ j, m as If the functions g i are given by where Q 1 is a constant, then the functions θ j, m can be written in the form sin Ω 1, j, 0 r 1 − Ω 1, j, 0 r 1 cos Ω 1, j, 0 r 1 K * * j,0 (t).
Taking the function θ j, m into account in Equation (59), we obtain the function φ i . This function and the function ψ i given by Equation (47) are used in Equation (21), which determines the temperature T i in the i-th layer of the body under consideration.

Numerical Examples and Discussion
The presented analytical solution of the fractional heat conduction problem was used to compute the temperature distributions in a layered spherical cone to investigate the effect of the phase-lags and the order of the fractional time-derivative on the heat conduction in the medium. The analysis concerns the heat conduction in the cone of constant initial temperature with an inner or outer heat source. We introduce the non-dimensional timē t = t κ n /b 2 , the non-dimensional radiir = r/b andr i = r i /b where r i determines the boundary of the i-th spherical layer, i = 1, . . . , n. Temperature T in the cone is defined as T(r, µ,t) = T i (r, µ, t) where r =r b, t =t b 2 /κ n and i is the layer number; i.e., the condition r i−1 < r ≤ r i is fulfilled. The spherical cone under consideration consists of five concentric layers whereinr i = 0.2 i, i = 1, . . . , 5. The physical data assumed in the computation are given in Table 1.  The role of the fractional order occurring in the mathematical model of the fractional heat conduction and its effect on the temperature distribution in the spherical cone was investigated for α ∈ [0.4; 0.6]. Computer simulation data were generated by the Mathematica software.
The first example concerns the heat conduction in a spherical cone without inner and outer heat sources. The cone is specified by the radius b = 0.5 m and the half-angle at the apex of the coneφ = π/3. The initial temperature, the ambient temperature and the heat transfer coefficient were assumed as: F i (r, µ) = T 0 = 20 • C, T a = 0 • C and a ∞ = 1200 W/(m 2 K), respectively. The computations were performed for the phase lags  In the case of the heat conduction in the spherical cone without an inner heat source, we predict that the function φ i occurring in Equation (21) satisfies the condition for i = 1, . . . , n and all r ∈ [0, b] and µ ∈ µ ϕ , 1 , i.e., the steady state of temperature distribution: T i (r, µ, t) = ψ i (r, µ) occurs in the cone. For the purpose of the formal proof of this statement, the asymptotic formula may be used for the Mittag-Leffler function E α,β given by Employing this formula in Equation (72), we obtain that lim t→∞ K * * j,m (t) = 1/q j,m . Taking into account the first term on the right-hand side of Equation (76) for consideration of the heat conduction without an inner source, we have lim t→∞ θ j, m (t) = 0 for j = 1, 2, . . . and m = 0, 1, 2, . . . Hence, on the basis of Equation (59), the condition (78) is fulfilled. Thus, the steady-state temperature distribution in the i-th layer of the solid sphere is given by Equation (21) in which the first term on the right-hand side can be omitted. In this case, the temperature in the considered region is independent from the phase-lags and the order of fractional derivatives which occur in the heat conduction model. In Figure 3, the 3D graph and contour plot of the function T i = ψ i (r, µ) illustrating the change of temperature in a hemisphere when its surface is heated by an outer heat source defined by Equation (50) with µ 0 = 1/2 and T a = 40 • C shown. The physical data that characterized the hemisphere, which were used in computation, are the same as in the first example.  The Robin condition (18) on the spherical surface of the cone describes the heat exchange with the environment. If the ambient temperature is higher than the temperature of the spherical surface of the cone, then the spherical surface is heated, and the heat flows in the cone in the radial direction. In Figure 4, curves are presented that illustrate the changes of temperatures in the cone when the ambient temperature is specified by Equation (50) with T 0 = 0 • C and T a = 20 • C. The calculations were performed assuming that there is no internal heat source. The next example concerns the spherical cone heated by the internal heat source described by Equation (76). The curves in Figure 5 present temperature distributions in the spherical cone for volumetric heat source Q 1 = 2 × 10 6 W/m 3 and four values of dimensionless time. Decreasing the temperature while increasing the radial distance from the heat source follows as a result of the heat exchange through the spherical surface of the cone with the environment of zero temperature. It should be noted that the temperatures in a part of the cone close to the heat source are higher for the higher fractional orders of the heat conduction models, and the temperatures in a part of the cone close to the spherical surface are higher for the lower orders of the conduction models.

Conclusions
The analytical solution of the heat conduction problem for a spherical cone with an internal and external heat source is presented. The problem formulation and its solution include also the special cases for the hemisphere and solid sphere. The fractional heat conduction model with a Caputo time-derivative includes phase-lags of the temperature and heat flux. The solution in the form of a double series of spherical Bessel functions and Legendre functions was derived. Numerical investigation of the effect of the fractional order of the time derivative occurring in the heat conduction equation is presented. According to the mathematical model, the fractional order of the differential equation has significant importance for temperature distribution in the spherical cone. It was stated that the higher order of fractional derivative in the heat conduction model leads to the higher temperature in the heated body and lower temperature in the cooled body. As was expected, for each value of the derivative order α, the temperature in all points of the cone tends with time to a steady state. However, as can be seen in presented examples, the temperature decreases/increases faster for larger values of the order occurring in the mathematical model of the fractional heat conduction.
In the future, we would like to develop a numerical model for solving the considered problem and make a comparative analysis of the obtained results. We think that in the case of a numerical solution, the constant values of parameters in the considered model can be replaced by functional dependencies, e.g., the temperature-dependent thermophysical parameters of materials in each layer. In addition, in the future research, we would like to include, among others, an external moving heat source in the heat conduction model presented in this paper and find analytical and/or numerical solutions for them. We would like to find the practical applications of the considered fractional heat conduction model, e.g., in modeling the influence of solar energy on the heating of the spherical cone [40].
Author Contributions: The authors (S.K., U.S., M.C.) claim to have contributed equally and significantly in this paper. All authors have read and agreed to the published version of the manuscript.