A Fractional Single-Phase-Lag Model of Heat Conduction for Describing Propagation of the Maximum Temperature in a Finite Medium

In this paper, an investigation of the maximum temperature propagation in a finite medium is presented. The heat conduction in the medium was modelled by using a single-phase-lag equation with fractional Caputo derivatives. The formulation and solution of the problem concern the heat conduction in a slab, a hollow cylinder, and a hollow sphere, which are subjected to a heat source represented by the Robotnov function and a harmonically varying ambient temperature. The problem with time-dependent Robin and homogenous Neumann boundary conditions has been solved by using an eigenfunction expansion method and the Laplace transform technique. The solution of the heat conduction problem was used for determination of the maximum temperature trajectories. The trajectories and propagation speeds of the temperature maxima in the medium depend on the order of fractional derivatives occurring in the heat conduction model. These dependencies for the heat conduction in the hollow cylinder have been numerically investigated.


Introduction
The classical Fourier's law of the heat conduction establishes the relationship between the heat flux vector and the gradient of the temperature [1] q(r, t) = −k∇T(r, t) where q is the heat flux vector, r is the point in the considered region, t is the time, k is the thermal conductivity of the material, ∇ is the gradient operator and T is the temperature. This relationship implies a nonphysical infinite speed of a thermal signal in the medium. To avoid this disagreement between the mathematical model and the observations, the single-phase-lag was introduced to the heat conduction model. Namely, the relationship (Equation (1)) is replaced by the following one [2] q(r, t + τ) = −k∇T(r, t) where τ is the phase lag of the heat flux. Expanding the left-hand side of Equation (2) into the Taylor series with respect to variable τ and taking into account two terms of this series (assuming that the phase lag τ is small), Equation (2) can be approximated by q(r, t) + τ ∂q ∂t (r, t) = −k∇T(r, t) Equation (3) was proposed by Cattaneo [3] and Vernotte [4] and currently it is known as the Cattaneo-Vernotte constitutive equation [5].
The heat flux vector occurring in Equation (3) can be eliminated by using the energy equation [6] − ∇·q(r, t) + g(r, t) = ρC p ∂T ∂t (4) where g(r, t) is the rate of the heat generation per unit volume, ρ is the density of the material and C p is the specific heat capacity. As a result, the single-phase-lag heat conduction equation is obtained: where κ = k/ ρC p is the thermal diffusivity. This hyperbolic type of equation describes the heat transfer in the wave form. The finite speed of the heat wave in the medium determines the square root of the ratio κ/τ. In the literature, there are numerous applications of the Equation (5) for modelling of the heat conduction. For instance, recently, works [7][8][9][10][11] have been published in which the hyperbolic model of heat conduction was applied. The description of transport processes with the use of fractional derivatives was proposed by Compte and Metzler [12]. In the presented mathematical model, a generalization of the Cattaneo-Vernotte constitutive equation was utilized. This generalization consists of replacing the time-derivative in the constitutive equation (Equation (3)) by the fractional derivative. The resulting generalized constitutive equation has the following form: q(r, t) + τ α ∂ α q ∂t α = −k∇T(r, t), 0 < α ≤ 1 (6) where ∂ α ∂t α denotes the fractional derivative with respect to variable t of order α. Scalar multiplying Equation (6) by the vector ∇ and then using Equation (4), the heat flux vector q can be eliminated. As a result, the fractional heat conduction equation is obtained: A large variety of generalizations of the constitutive equation for the heat transfer and their applications are presented in the literature [13][14][15][16][17][18][19].
A generalization of the heat conduction model can be obtained by replacing the time derivative in Equation (4) with the fractional derivative of order β. The obtained generalized fractional energy equation has the form In this equation, the coefficient ν β−1 is introduced to keep the accordance of dimensions. The generalized constitutive equation (Equation (6)) in combination with the generalized energy equation (Equation (8)), results in the single-phase-lag heat conduction equation with fractional derivatives The notation of the first term on the left-hand side of Equation (9) is dictated by the fact that in general, the fractional derivatives are noncommutative operators [20]. It can be noted that this equation for β = 1 is of the form of fractional Equation (7), and for α = β = 1, it has the form of hyperbolic Equation (5). In the literature, there are no works devoted to the propagation of the maximum temperature in a medium based on the heat conduction model in which the non-local and phase-lag properties are considered.
In applications of the fractional calculus, the Riemann-Liouville and the Caputo derivatives are often used. The definitions and properties of these fractional derivatives are presented in the books by Diethelm [21], Kilbas et al. [22], Mainardi [20], Podlubny [23] and Povstenko [24]. In this paper, the Caputo derivative and its properties will be used. The fractional Caputo derivative of order α with the lower limit zero, is defined as From Equation (10), it follows that the fractional derivatives occurring in the differential equation describing the state of the system contain information about its past state. This non-local property of the fractional derivatives is the important advantage of using fractional calculus in mathematical modeling.
The fractional differential Equation (9) is completed by initial and boundary conditions. A solution to this initial-boundary value problem is the temperature distribution as a function of time and space variables. This function for a fixed time variable can achieve a local maximum value with respect to the space variable. The point of the maximum propagates with a finite speed in the considered region. The propagation problem of the maximum point of a fundamental solution to a fractional equation was considered by Luchko et al. [25]. The presented results concern the Cauchy problem for a one-dimensional time-fractional diffusion-wave equation in an unbounded region.
In this paper, a solution to the heat conduction problem according to the time-fractional single-phase-lag model is presented. The solution in the Laplace transform domain includes the one-dimensional heat conduction in a slab, a hollow cylinder, and a hollow sphere. The obtained temperature distribution for the tracking of the propagation of the maximum temperature in the considered region was used. The presented numerical results concern the hollow cylinder with the Robin-Neumann boundary conditions, which is subjected to a variable ambient temperature or impulsive heat source.

Formulation of the Problem
Let us consider the heat conduction governed by the time-fractional differential equation (Equation (9)). This equation is valid in the region which is specified by a medium in the space. We will deal with heat conduction in a slab, a hollow cylinder, and a hollow sphere. In the each of the three cases, assuming one dimensional heat conduction, we denote the space variable by "x" where a ≤ x ≤ b. For the slab, the heat conduction in the direction of x-axis of a rectangular coordinate system is considered, for the cylinder-in a radial direction of a cylindrical coordinate system and for the sphere-in a radial direction of a spherical coordinate system. The operator ∇ 2 in Equation (9) for the slab, cylinder and sphere can be written in the form [1] where p = 0 for the slab, p = 1 for the cylinder and p = 2 for the sphere. Equation (9) is complemented by boundary and initial conditions. We assume the Robin-Neumann boundary conditions: ∂T ∂x (b, t) = 0 (13) and the following initial conditions: where h a is the convective heat transfer coefficient, T a (t) is the ambient temperature, f (x) is the initial temperature, and h(x) is the fractional time-derivative of order β of the temperature at an initial time.
The function T(x, t) is a solution of the initial-boundary problem (Equations (9) and (12)-(15)) with a non-homogenous boundary condition. In the first stage of solving this problem, we present the function T(x, t) in the form of the sum The function θ(x, t) satisfies the non-homogenous differential equation and the following homogenous boundary conditions The initial conditions for the function θ(x, t) are obtained using Equations (14)-(16) The function θ(x, t) as a solution of the initial-boundary problem (Equations (17)-(22)), will be determined in the form of an orthogonal series, and then the Laplace transform technique will be used.

Solution to the Problem
We can search for a solution to the problem (Equations (17)- (22)) in the form of a series: where the functions Φ i (x) are solutions of the following eigenproblem The general solution to Equation (24) can be written in the form where A, B are constants and the functions ϕ(x) and ψ(x) are independent, particular solutions to this equation. Taking into account the boundary conditions (Equations (25) and (26)) and using the standard procedure, we obtain the eigenvalue equation which is solved with respect to λ. The obtained roots create an infinite sequence of eigenvalues: λ i , i = 1, 2, . . .. In turn, assuming A i = ψ i (b) and using Equations (26) and (27) for The functions ϕ i (x) and ψ i (x) for the three cases of the operator ∇ 2 given by Equation (11) are presented in Table 1. Table 1. The functions ϕ i (x) and ψ i (x) for the slab (p = 0), hollow cylinder (p = 1) and hollow sphere (p = 2).
where the normalization integrals N i are determined according to the formula The eigenfunctions Φ i (x), eigenvalue equations, and normalization integrals N i for the eigenvalue problem (Equations (24)-(26)) for a slab, a hollow cylinder, and a hollow sphere, are presented in Table 2. The presented approach can be applied to the fractional single-phase-lag heat conduction problem obtained by replacing the Neumann boundary condition (Equation (13)) with the homogeneous Dirichlet boundary condition: T(b, t) = 0.
In order to derive an equation that will be used to determine the functions Λ i (t) occurring in Equation (23), we substitute the series (Equation (23)) into Equation (17), then we multiply the resulting equation by the function x p Φ j (x) and we integrate it over the interval [a, b]. Using the condition from Equation (30), we obtain the equation in the form Similarly, multiplying both sides of Equations (21) and (22) by x p Φ j (x) and integrating over the interval [a, b], the following initial conditions are obtained: We find a solution to the initial problem (Equations (32)-(34)) by using the Laplace transform where s is a complex parameter. We utilized the linearity property of the Laplace transform and the following rule [20]: Using the rule (Equation (36)), the Laplace transform of the solution to the problem (Equations (32)-(34)), after some transformation, can be presented in the form and is the Laplace transform of the function Q(t) given by Equation (20b). The complete solution to the problem in the Laplace transform domain will be determined for an established function describing the rate of the heat generation g(x, t) and functions occurring in the initial and boundary conditions: and T a (t).
We find the Laplace transform G(x, s) assuming the function g(x, t) in Equation (20a) in the form where g 0 is the strength of the heat source per unit length of the surface [26], δ(·) is the Dirac delta function, F α (λ, z) is the Robotnov function [27]: and E α,β (z) is the two-parameter Mittag-Leffler function defined by the power series The Robotnov function is called the "impulse response" of the fundamental fractional order differential equation because it satisfies the differential equation [28] Using Equations (38), (39) and (43) and the Laplace transform pair [24] L the Laplace transform of the function G(x, t) defined by Equation (20a), can be written as The functions T a (t), f (x) and h(x) occurring in the boundary and initial conditions (12) and in Equations (14) and (15), we assume that T a (t) = P 1 + P 2 sin ωt, t ≥ 0 (46a) Considering Equations (39) and (45), in Equation (37) we obtain the Laplace transform Λ i (s) in the form where For the purpose of deriving the inverse Laplace transform L −1 Λ i (s) , we use the convolution theorem [20] L Introducing the function U and using the properties in Equation (49), we can write the inverse Laplace transform Λ i (t) in the form Finally, the temperature distribution T(x, t) is given by Equations (16), (23), (27) and (51). The function U α,β,γ i (t), as an inverse Laplace transform, can be determined in an analytical form but only for some values of the fractional orders α and β. For example, if α = 0.5, β = 1.0 and γ = 0; 0.5, this function can be written in the form For γ = −0.5, the sum in Equation (52) should be complemented with a term τ α κλ 2 i √ π .

Numerical Analysis and Discussion
The temperature distribution in the medium is given by the formula which contains the inverse Laplace transform (Equation (50)). This inverse Laplace transform for established values of α, β, and γ can be determined numerically. In the literature, many different algorithms are available for numerical Laplace inversion [29][30][31]. In order to find an effective algorithm for precise numerical inversion of the Laplace transforms appearing in the presented solution, several algorithms were tested. On the basis of these numerical tests, the Fixed-Talbot algorithm for further computations was chosen. Applying this algorithm, the values of a function f (t) = L −1 f (s) are computed using the formula [29] where µ(θ) = pθ(ctgθ + i), σ(θ) = θ + (θctgθ − 1)ctgθ, p = 2M/(5t), θ k = kπ/M, i = √ −1 and M is the number of precision decimal digits.
Numerical results computed using the Fixed-Talbot procedure were compared with those obtained by using of the analytical form of the inverse Laplace transform (Equation (52)) for the function U α,β,γ 1 (t) with α = 0.5, β = 1.0 and γ = −0.5; 0.0; 0.5. In Table 3, absolute values of the relative errors   (a, b). The maximum temperature location moves in the medium with a finite speed in the direction appointed by the decreasing temperature gradient. The numerical analysis presented in this section concerns the problem of propagation of the maximum temperature in a finite hollow cylinder that is heated by a heat source or through an operation of variable ambient temperature. Numerical calculations were performed to obtain the following data: The inner and outer radii of the cylinder are a = 0.4 m and b = 0.6 m, respectively, the thermal diffusivity is κ = 8.418×10 −5 m 2 /s, the thermal conductivity is k = 204 W/(m·K), the heat transfer coefficient at the inner surface of the cylinder is h a = 800 W/ m 2 K , and the ambient temperature at the initial time t = 0 is T a (0) = 100 • C. In obtaining the numerical results, the following non-dimensional quantities were used: . The computations were carried out using the Mathematica package [32].
Let us consider the hollow cylinder heated at the inner surface by the heat source described by the Robotnov function specified by Equation (40)    Assuming that an operation of the Robotnov heat source defined by Equation (40), the non-dimensional temperatureT x,t for a fixed valuet takes a maximal value atx ∈ (0, 1) if the following condition is fulfilled: Solving this equation with respect tox fort > 0, we obtain a curve of locations of the maxima temperatures in the plane Otx. These curves, for different orders of derivatives α and β, are presented in Figure 2. The results indicate an important significance of the orders of fractional derivatives occurring in the heat conduction equation for the time of the propagation of the maximal temperature in the cylinder. The time of the propagation of the maximal temperature is significantly Shorter for higher values of the derivative orders in the heat conduction model. The curves presented in Figure 2 show that the replacement of the Caputo derivative order β = 0.9 by β = 1.0 leads to a shortening by half of the transition time of the maximum temperature in the cylinder fromx = 0 tox = 0.5, i.e., the change of the derivative order β results in a change of speed of the propagation of the maximum temperature in the medium.
The non-dimensional speed of propagation of the maximum temperaturev is given by the formula v(t) = dx(t)/dt, wherex(t) is an implicit function defined by Equation (54). Differentiating both sides of Equation (54) with respect tot and using Equation (16), we find the derivative of the functionx(t) in the formv  The solution presented in the previous section includes a case of the fractional heat conduction in a hollow cylinder when the temperature inside the cylinder changes harmonically according to Equation (46a). The numerical computation of the temperature distribution in the hollow cylinder was performed assuming that no other heat sources occur. In Figure 4, the 3D graphs and contour  The solution presented in the previous section includes a case of the fractional heat conduction in a hollow cylinder when the temperature inside the cylinder changes harmonically according to Equation (46a). The numerical computation of the temperature distribution in the hollow cylinder was performed assuming that no other heat sources occur. In Figure 4, the 3D graphs and contour  The solution presented in the previous section includes a case of the fractional heat conduction in a hollow cylinder when the temperature inside the cylinder changes harmonically according to Equation (46a). The numerical computation of the temperature distribution in the hollow cylinder was performed assuming that no other heat sources occur. In Figure 4, the 3D graphs and contour plots of the functionT x,t for P 2 /P 1 = 0.5, ω = 0.005s −1 , β = 0.9 and different values of α are presented. The maxima and minima of temperatures propagate in the hollow cylinder from the inner to outer boundary. The amplitude of the temperature decreases with the space variable in all cases of the values of α. The higher amplitude of the temperature changes at the outer boundary of the cylinder occur for the heat conduction model with the higher order fractional differential equation.

Conclusions
The solution to the heat conduction problem based on the fractional single-phase-lag model was derived. The formulation and solution of the problem concern the heat conduction in a slab, a hollow cylinder, and a hollow sphere. The considered 1D problem with time-dependent Robin and homogenous Neumann boundary conditions were solved by the use of the eigenfunction expansion method with respect to the space variable and the Laplace transform technique with respect to time. It was assumed that the medium is exposed to a heat source represented by the Robotnov function and the sinusoidaly changing ambient heat. The derived temperature distribution in a hollow cylinder helped in the numerical investigation of propagation of the maximum temperature. It was stated that fractional orders of the Caputo derivatives occurring in the differential equation governing the heat conduction have a significant influence on the propagation of the maximum temperature in the cylinder subjected to the Robotnov heat source. The trajectories of the maximum temperatures show that a slight increase of the fractional derivative orders can cause a considerable decrease of the occurrence time of the maximum temperature in the cylinder. It was observed that the propagation speed of the maximum temperature in the cylinder is higher for a higher fractional order of the differential heat conduction equation. The maximum and minimum temperature propagates in the medium also when the ambient temperature changes harmonically. The amplitude

Conclusions
The solution to the heat conduction problem based on the fractional single-phase-lag model was derived. The formulation and solution of the problem concern the heat conduction in a slab, a hollow cylinder, and a hollow sphere. The considered 1D problem with time-dependent Robin and homogenous Neumann boundary conditions were solved by the use of the eigenfunction expansion method with respect to the space variable and the Laplace transform technique with respect to time. It was assumed that the medium is exposed to a heat source represented by the Robotnov function and the sinusoidaly changing ambient heat. The derived temperature distribution in a hollow cylinder helped in the numerical investigation of propagation of the maximum temperature. It was stated that fractional orders of the Caputo derivatives occurring in the differential equation governing the heat conduction have a significant influence on the propagation of the maximum temperature in the cylinder subjected to the Robotnov heat source. The trajectories of the maximum temperatures show that a slight increase of the fractional derivative orders can cause a considerable decrease of the occurrence time of the maximum temperature in the cylinder. It was observed that the propagation speed of the maximum temperature in the cylinder is higher for a higher fractional order of the differential heat conduction equation. The maximum and minimum temperature propagates in the medium also when the ambient temperature changes harmonically. The amplitude of these changes in the medium decrease with increasing distance from the heated boundary. This decreasing of the amplitudes of the temperature oscillations is greater for smaller orders of the fractional derivative in the heat conduction model. This observation leads to a physical interpretation of the parameter α as a thermal damping coefficient in the fractional heat conduction model.