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.

1. 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] 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] 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

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] where is the rate of the heat generation per unit volume, is the density of the material and is the specific heat capacity. As a result, the single-phase-lag heat conduction equation is obtained: where 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: where 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 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 is of the form of fractional Equation (7), and for , 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.

2. 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 . 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 in Equation (9) for the slab, cylinder and sphere can be written in the form [1] where for the slab, for the cylinder and for the sphere.

Equation (9) is complemented by boundary and initial conditions. We assume the Robin–Neumann boundary conditions: and the following initial conditions: where is the convective heat transfer coefficient, is the ambient temperature, is the initial temperature, and is the fractional time-derivative of order β of the temperature at an initial time.

The function 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 in the form of the sum

The function satisfies the non-homogenous differential equation and the following homogenous boundary conditions where

The initial conditions for the function are obtained using Equations (14)–(16)

The function 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.

3. 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 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 and 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: , . In turn, assuming and using Equations (26) and (27) for , , we find . Hence, the eigenfunctions corresponding to the eigenvalues can be rewritten as

The functions and for the three cases of the operator given by Equation (11) are presented in Table 1.

Table 1

The functions and for the slab (), hollow cylinder () and hollow sphere ().

p	φi (x)	ψi (x)
0	cos(λ_ix)	sin(λ_ix)
1	J_0(λ_ix)	Y_0(λ_ix)
2	cos(λ_ix)x	sin(λ_ix)x

The functions satisfy the orthogonality condition in the form where the normalization integrals are determined according to the formula

The eigenfunctions , eigenvalue equations, and normalization integrals 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: .

Table 2

The eigenfunctions , eigenvalue equations and normalization integrals for the slab (), hollow cylinder () and hollow sphere ().

p	Eigenfunction	Eigenequation	Normalization Integral
0	Φ_i(x)=cos(b−x)λ_i	sin(b−a)λ−h_akλcos(b−a)λ=0	N_i=b−a2(1+sin2λ_i(b−a)2λ_i(b−a))
1	Φ_i(x)=Y_1(bλ_i)J_0(λ_ix) −J_1(bλ_i)Y_0(λ_ix)	h_akλ(J_1(bλ)Y_0(aλ)−J_0(aλ)Y_1(bλ))+J_1(bλ)Y_1(aλ)−J_1(aλ)Y_1(bλ)=0	N_i=2π^2λi2−a^22(1+(h_akλ_i)^2)·(J_1(bλ_i)Y_0(aλ_i)−J_0(aλ_i)Y_1(bλ_i))^2
2	Φ_i(x)=1x(cosλ_i(b−x) −1bλ_isinλ_i(b−x))	(1+kh_a(1a+1b(bλ)^2))sin(b−a)λbλ−(1+kh_a(1a−1b))cos(b−a)λ=0	N_i=b^2λi2−14b^2λi3sin2λ_i(b−a)+12b^2λi2·((b−a)(bλ_i)^2−a+bcos2λ_i(b−a))

In order to derive an equation that will be used to determine the functions occurring in Equation (23), we substitute the series (Equation (23)) into Equation (17), then we multiply the resulting equation by the function and we integrate it over the interval . Using the condition from Equation (30), we obtain the equation in the form

Similarly, multiplying both sides of Equations (21) and (22) by and integrating over the interval , the following initial conditions are obtained:

We find a solution to the initial problem (Equations (32)–(34)) by using the Laplace transform technique. The Laplace transform of a function is defined as where 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 where and is the Laplace transform of the function 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 and functions occurring in the initial and boundary conditions: , and .

We find the Laplace transform assuming the function in Equation (20a) in the form where is the strength of the heat source per unit length of the surface [26], is the Dirac delta function, is the Robotnov function [27]: and 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] the Laplace transform of the function defined by Equation (20a), can be written as

The functions , and occurring in the boundary and initial conditions (12) and in Equations (14) and (15), we assume that

Considering Equations (39) and (45), in Equation (37) we obtain the Laplace transform in the form where

For the purpose of deriving the inverse Laplace transform , we use the convolution theorem [20]

Introducing the function defined as and using the properties in Equation (49), we can write the inverse Laplace transform in the form

Finally, the temperature distribution is given by Equations (16), (23), (27) and (51).

The function , 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 and , this function can be written in the form where and are roots of equation: . For , the sum in Equation (52) should be complemented with a term .

4. 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 are computed using the formula [29] where , , , , and 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 with and . In Table 3, absolute values of the relative errors are presented where are obtained by using the numerical inversion of the Laplace transform and are values of the function in Equation (52). The small relative errors justify the use of the Fixed-Talbot procedure for numerical inversion of the Laplace transform given by Equation (50).

Table 3

The relative errors of the results obtained by using the Fixed-Talbot procedure and exact values of the function for , and .

t^=κt/(b−a)^2	γ=−0.5	γ=0	γ=0.5
0.5	3.36748×10^−7	2.22232×10^−6	2.52736×10^−5
1	4.85927×10^−7	2.84182×10^−6	5.47403×10^−5
1.5	2.51907×10^−7	1.83325×10^−5	5.94386×10^−4
2	5.96577×10^−7	5.66522×10^−6	1.12375×10^−3
2.5	6.85496×10^−7	7.71945×10^−6	1.77174×10^−4

The temperature for fixed time t as a function of the space variable x, defined on the interval , can achieve a maximum value in the open interval . 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 and , respectively, the thermal diffusivity is , the thermal conductivity is , the heat transfer coefficient at the inner surface of the cylinder is , and the ambient temperature at the initial time t = 0 is . 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) with . The graphs of the temperature distributions as functions of the space variable for different moments of time are shown in Figure 1. For a fixed dimensionless time , the function of variable assumes a maximum at a point . The thick red line in Figure 1 is created by the points which are observed at different times . The points of the maximum temperatures propagate with time in the direction of the region of lower temperature. The temperature of the cylinder decreases with time, therefore the maxima, after some time, are small. For this reason, the observation of the temperature maxima is limited to the interval of non-dimensional space variable . The presented curves were obtained for , , , and .

Figure 1

Temperature distributions in the hollow cylinder as functions of the space variable for , and different dimensionless time .

Assuming that an operation of the Robotnov heat source defined by Equation (40), the non-dimensional temperature for a fixed value takes a maximal value at if the following condition is fulfilled:

Solving this equation with respect to for , we obtain a curve of locations of the maxima temperatures in the plane . 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

Figure 2

The curves of locations of maxima temperatures in the hollow cylinder for different values of parameter ; (a) ; (b) ; (c) .

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 by leads to a shortening by half of the transition time of the maximum temperature in the cylinder from to , 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 temperature is given by the formula (t) = dx(t)/dt, where is an implicit function defined by Equation (54). Differentiating both sides of Equation (54) with respect to and using Equation (16), we find the derivative of the function in the form whereas , and . The curves of propagation speed of the maximum temperature in the hollow cylinder subjected to the Robotnov heat source for and different α values, are presented in Figure 3. It can be noticed that the propagation speed of the maximum temperature strongly depends on the derivative order α, occurring in the fractional heat conduction equation.

Figure 3

The curves of the propagation speed of maxima temperatures in the hollow cylinder for and different values of parameter .

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 function for , , 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.

Figure 4

Temperature distributions in the hollow cylinder as functions of non-dimensional space variable and time when the temperature inside the hollow cylinder changes sinusoidaly for the fractional derivative order ; (a) ; (b) ; (c) ; (d) .

5. 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.