Time-Fractional Heat Conduction in a Plane with Two External Half-Infinite Line Slits under Heat Flux Loading

Abstract

The time-fractional heat conduction equation follows from the law of conservation of energy and the corresponding time-nonlocal extension of the Fourier law with the “long-tail” power kernel. The time-fractional heat conduction equation with the Caputo derivative is solved for an infinite plane with two external half-infinite slits with the prescribed heat flux across their surfaces. The integral transform technique is used. The solution is obtained in the form of integrals with integrand being the Mittag–Leffler function. A graphical representation of numerical results is given.

1. Introduction

The conventional Fourier law states the linear dependence between the heat flux vector $\mathbf{q}$ and the gradient of temperature:

$$\mathbf{q}=-k\mathrm{grad}T,$$

(1)

with k being the thermal conductivity of a solid.

The Fourier law (1) in combination with the law of conservation of energy

$$C\frac{\partial T}{\partial t}=-\mathrm{div}\mathbf{q}$$

(2)

leads to the classical parabolic heat conduction equation

$$\frac{\partial T}{\partial t}=a\Delta T,$$

(3)

where $a=k/C$ denotes the thermal diffusivity coefficient, C is the heat capacity,

In non-classical theories, the Fourier law (1) and the heat conduction Equation (3) are substituted by alternative equations which give an account of the complex internal structure of a body and reflect processes taking place at the microscopic level. Gurtin and Pipkin [1] generalized the Fourier law to time-nonlocal dependence between the heat flux vector $\mathbf{q}$ and the gradient of temperature. Nigmatullin [2] proposed the time-nonlocal extension of the Fourier law in the form:

$$\mathbf{q}\left(t\right)=-k{\int }_0^{t}K(t-\tau )\mathrm{grad}T\left(\tau \right)\mathrm{d}\tau ,$$

(4)

with $K(t-\tau )$ being the time-nonlocality kernel. With the law of conservation of energy (2), the constitutive Equation (4) gives the heat conduction equation with memory [2]

$$\frac{\partial T}{\partial t}=a{\int }_0^{t}K(t-\tau )\Delta T\left(\tau \right)\mathrm{d}\tau .$$

(5)

It is obvious that the standard Fourier law (1) and the parabolic heat conduction Equation (3) are obtained when

$$K(t-\tau )=\delta (t-\tau ),$$

(6)

where $\delta (t-\tau )$ is the Dirac delta function.

Different extended theories of heat conduction correspond to the appropriate selection of the memory kernel $K(t-\tau )$. “Full sclerosis” (localized heat conduction) [3,4] is described by the memory kernel being the time derivative of the Dirac delta function

$$K(t-\tau )=-{\delta }^{\prime }(t-\tau )$$

(7)

and results in the Helmholtz equation for temperature

$$T=a\Delta T.$$

(8)

“Full memory” [5,6] infers that there is no fading of memory, the kernel $K(t-\tau )$ is constant

$$K(t-\tau )=1,$$

(9)

and temperature obeys the wave equation

$$\frac{{\partial }^{2}T}{\partial t^2}=a\Delta T.$$

(10)

“Short-tail” memory [7,8] is associated with the exponential time-nonlocality kernel and leads to the Cattaneo telegraph equation for temperature. “Middle-tail” memory [9,10] involves kernels being different Mittag–Leffler functions and leads to the time-fractional telegraph equations.

The time-nonlocal dependences between the heat flux vector $\mathbf{q}$ and the gradient of temperature with the “long-tail” power kernel [3,4,11,12,13,14]

$$\mathbf{q}\left(t\right)=-\frac{k}{\mathsf{\Gamma }\left(\alpha \right)}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^t{(t-\tau )}^{\alpha -1}\mathrm{grad}T\left(\tau \right)\mathrm{d}\tau ,0<\alpha \le 1,$$

(11)

$$\mathbf{q}\left(t\right)=-\frac{k}{\mathsf{\Gamma }(\alpha -1)}{\int }_0^t{(t-\tau )}^{\alpha -2}\mathrm{grad}T\left(\tau \right)\mathrm{d}\tau ,1<\alpha \le 2,$$

(12)

can be interpreted in terms of fractional calculus:

$$\mathbf{q}\left(t\right)=-k{D}_{RL}^{1-\alpha }\mathrm{grad}T\left(t\right),0<\alpha \le 1,$$

(13)

$$\mathbf{q}\left(t\right)=-k{I}^{\alpha -1}\mathrm{grad}T\left(t\right),1<\alpha \le 2,$$

(14)

and in combination with the law of conservation of energy (2) result in the time-fractional heat conduction equation

$$\frac{{\partial }^{\alpha }T}{\partial t^{\alpha }}=a\Delta T,0<\alpha \le 2,$$

(15)

with the Caputo fractional derivative.

Here, $I^{\alpha }f\left(t\right)$ and $D_{RL}^{\alpha }f\left(t\right)$ are the Riemann–Liouville fractional integral and the Riemann–Liouville fractional derivative [15,16]

$$I^{\alpha }f\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-\tau )}^{\alpha -1}f\left(\tau \right)\mathrm{d}\tau ,\alpha >0,$$

(16)

$$D_{RL}^{\alpha }f\left(t\right)=\frac{{\mathrm{d}}^n}{\mathrm{d}{t}^n}\left[\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^t{(t-\tau )}^{n-\alpha -1}f\left(\tau \right)\mathrm{d}\tau \right],n-1<\alpha <n.$$

(17)

$\frac{{\mathrm{d}}^{\alpha }f\left(t\right)}{\mathrm{d}{t}^{\alpha }}$ is the Caputo fractional derivative

$$\frac{{\mathrm{d}}^{\alpha }f\left(t\right)}{\mathrm{d}{t}^{\alpha }}=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^t{(t-\tau )}^{n-\alpha -1}\frac{{\mathrm{d}}^{n}f\left(\tau \right)}{\mathrm{d}{\tau }^n}\mathrm{d}\tau ,n-1<\alpha <n.$$

(18)

$\mathsf{\Gamma }\left(\alpha \right)$ denotes the Gamma function.

Equation (15) combines the whole spectrum of equations starting with the localized heat conduction (the elliptic Helmholtz equation for temperature when $\alpha \to 0$) through the standard parabolic heat conduction equation (when $\alpha =1$) and ending with the ballistic heat conduction (the hyperbolic wave equation for temperature when $\alpha =2$).

Fractional calculus (the theory of integrals and derivatives of non-integer order) has many important applications in physics, geophysics, chemistry, continuum mechanics, rheology, engineering, biology, etc. (see, for example, Refs. [17,18,19,20,21,22,23,24,25,26,27] and references therein).

Starting from the fundamental papers by Wyss [28] and Mainardi [29,30], considerable interest has been shown in solving the time-fractional heat conduction equation (diffusion-wave equation) (15). The book [31] sums up investigations in this field.

Strength of solids depends significantly on the presence in a real solid such defects as cracks, slits, holes, inclusions, etc. During deformation in the vicinity of such defects, there appears considerable stress concentration which can lead to formation of new and growth of existing crack and, hence, to local or global fracture of a solid. Thermal shock conditions which can arise during sudden cooling of a solid can result in very high thermal stresses near a crack, therefore, thermal shock is an important fracture mechanism of brittle materials. The pioneering papers [32,33] laid the groundwork for study of thermoelastic problems for solids with cracks and slits based on the classical heat conduction Equation (3). The literature on this subject is quite voluminous (see, for example, Refs. [34,35,36,37,38,39,40,41,42] and references therein). More general heat conduction equations imply formulation of associated theories of thermal stresses. Thermoelasticity associated with the time-fractional heat conduction Equation (15) was proposed in [11]; the book [26] sums up investigations in this field. The first paper on fractional heat conduction in a solid with a crack was published by the authors [43]. In the present paper, the time-fractional heat conduction equation with the Caputo derivative is solved for an infinite plate with two external half-infinite slits whose surfaces are exposed to the heat flux loading. The thermal stress aspect of the problem including investigation of the stress intensity factor will be studied in future research.

2. Statement of the Problem

Consider a plane with two slits located on the half-lines $y=0$, $-\infty <x<-l$ and $y=0$, $l<x<\infty$. To study the heat conduction problem, we examine the time-fractional heat conduction equation

$$\begin{array}{c}\frac{{\partial }^{\alpha }T}{\partial t^{\alpha }}=a\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right),\hfill \\ -\infty <x<\infty ,-\infty <y<\infty ,0<t<\infty ,0<\alpha \le 2,\hfill \end{array}$$

(19)

and assume zero initial conditions

$$\begin{array}{c}t=0:T=0,0<\alpha \le 2,\hfill \end{array}$$

(20)

$$\begin{array}{c}t=0:\frac{\partial T}{\partial t}=0,1<\alpha \le 2,\hfill \end{array}$$

(21)

and heat flux loading according to the constitutive Equations (13) and (14)

$$\begin{array}{c}y=0^{+}:\phantom{-}k{D}_{RL}^{1-\alpha }\frac{\partial T}{\partial y}=q^{+}(x,t),l<\left|x\right|<\infty ,0<\alpha \le 1,\hfill \end{array}$$

(22)

$$\begin{array}{c}y=0^{+}:\phantom{-}k{I}^{\alpha -1}\frac{\partial T}{\partial y}=q^{+}(x,t),l<\left|x\right|<\infty ,1<\alpha \le 2,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & y=0^{-}:-k{D}_{RL}^{1-\alpha }\frac{\partial T}{\partial y}=q^{-}(x,t),l<\left|x\right|<\infty ,0<\alpha \le 1,\hfill \end{array}$$

(24)

$$\begin{array}{c}y=0^{-}:-k{I}^{\alpha -1}\frac{\partial T}{\partial y}=q^{-}(x,t),l<\left|x\right|<\infty ,1<\alpha \le 2.\hfill \end{array}$$

(25)

The following conditions at infinity are imposed

$$\underset{y\to \pm \infty }{lim}T(x,y,t)=0,$$

(26)

$$\underset{t\to \infty }{lim}T(x,y,t)=0.$$

(27)

In the subsequent text, we will restrict our consideration to the problem with constant heat flux loading

$$q^{+}(x,t)=q^{-}(x,t)=q_0=\mathrm{const}.$$

(28)

In this case, due to the symmetry of the problem with respect to the x-axis, we can study only the upper half-plane $y\ge 0$. The problem is reformulated as follows: solve the heat conduction equation

$$\begin{array}{c}\frac{{\partial }^{\alpha }T}{\partial t^{\alpha }}=a\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right),\hfill \\ -\infty <x<\infty ,0<y<\infty ,0<t<\infty ,0<\alpha \le 2,\hfill \end{array}$$

(29)

with zero initial conditions

$$\begin{array}{c}t=0:T=0,0<\alpha \le 2,\hfill \end{array}$$

(30)

$$\begin{array}{c}t=0:\frac{\partial T}{\partial t}=0,1<\alpha \le 2,\hfill \end{array}$$

(31)

under constant heat flux loading

$$\begin{array}{c}y=0:k{D}_{RL}^{1-\alpha }\frac{\partial T}{\partial y}=q_0,l<\left|x\right|<\infty ,0<\alpha \le 1,\hfill \end{array}$$

(32)

$$\begin{array}{c}y=0:k{I}^{\alpha -1}\frac{\partial T}{\partial y}=q_0,l<\left|x\right|<\infty ,1<\alpha \le 2,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & y=0:\frac{\partial T}{\partial y}=0,0\le \left|x\right|<l,0<\alpha \le 2,\hfill \end{array}$$

(34)

and zero conditions at infinity

$$\underset{y\to \infty }{lim}T(x,y,t)=0,$$

(35)

$$\underset{t\to \infty }{lim}T(x,y,t)=0.$$

(36)

3. Solution of the Problem

To solve the problem, we will employ the Laplace transform with respect to time t, the exponential Fourier transform with respect to the spatial coordinate x, and the cos-Fourier transform with respect to the spatial coordinate y.

Recall the Laplace transform rules for fractional integrals and derivatives [15,16]:

$$\mathcal{L}\left\{I^{\alpha }f\left(t\right)\right\}=\frac{1}{s^{\alpha }}{f}^{*}\left(s\right),$$

(37)

$$\mathcal{L}\left\{{\mathrm{D}}_{RL}^{\alpha }f\left(t\right)\right\}=s^{\alpha }{f}^{*}\left(s\right)-\sum _{k=0}^{n-1}{D}^k{I}^{n-\alpha }f\left(0^{+}\right)s^{n-1-k},n-1<\alpha <n,$$

(38)

$$\mathcal{L}\left\{\frac{{\mathrm{d}}^{\alpha }f\left(t\right)}{\mathrm{d}{t}^{\alpha }}\right\}=s^{\alpha }{f}^{*}\left(s\right)-\sum _{k=0}^{n-1}{f}^{\left(k\right)}\left(0^{+}\right)s^{\alpha -1-k},n-1<\alpha <n.$$

(39)

Here, s denotes the Laplace transform variable; the asterisk marks the transform.

The integral transform technique results in the solution in the transform domain

$$\tilde{\tilde{T}}{}^{*}(\xi ,\eta ,s)=-\frac{\sqrt{2}a{q}_0}{k\sqrt{\pi }}\left[\pi \delta \left(\xi \right)-\frac{sin\left(l\xi \right)}{\xi }\right]\frac{s^{\alpha -2}}{s^{\alpha }+a\left({\xi }^2+{\eta }^2\right)}.$$

(40)

Here, both Fourier transforms are marked by the tilde, $\xi$ is the transform variable of the exponential Fourier transform, $\eta$ is the transform variable of the cos-Fourier transform, $\delta \left(\xi \right)$ is the Dirac delta function appearing in the integral

$${\int }_{-\infty }^{\infty }cos\left(x\xi \right)\mathrm{d}x=2\pi \delta \left(\xi \right).$$

(41)

Inverting the integral transforms with taking into account the formula for the inverse Laplace transform [15,16]

$${\mathcal{L}}^{-1}\left\{\frac{s^{\alpha -\beta }}{s^{\alpha }+b}\right\}=t^{\beta -1}{E}_{\alpha ,\beta }(-b{t}^{\alpha }),$$

(42)

where $E_{\alpha ,\beta }\left(z\right)$ is the Mittag–Leffler function in two parameters $\alpha$ and $\beta$

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathsf{\Gamma }(\alpha k+\beta )},\alpha >0,\beta >0,z\in C,$$

(43)

we get the solution

$$\begin{array}{c}T(x,y,t)=-\frac{2a{q}_{0}t}{\pi k}{\int }_0^{\infty }{E}_{\alpha ,2}\left(-a{\eta }^2{t}^{\alpha }\right)cos\left(y\eta \right)\mathrm{d}\eta \hfill \\ +\frac{2a{q}_{0}t}{{\pi }^{2}k}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }{E}_{\alpha ,2}\left[-a\left({\xi }^2+{\eta }^2\right)t^{\alpha }\right]\frac{sin\left(l\xi \right)}{\xi }cos\left(x\xi \right)cos\left(y\eta \right)\mathrm{d}\xi \mathrm{d}\eta .\hfill \end{array}$$

(44)

To simplify the second term in Equation (44), we will use the polar coordinate system in the $(\eta ,\xi )$-plane:

$$\xi =\rho cos\theta ,\eta =\rho sin\theta$$

(45)

and the substitution $u=cos\theta$. Having regard to the following integral [31]

$${\int }_0^1\frac{sin\left(qw\right)cos\left(p\sqrt{1-w^2}\right)}{w\sqrt{1-w^2}}\mathrm{d}w=\frac{\pi }{2}{\int }_0^q{J}_0\left(\sqrt{p^2+z^2}\right)\mathrm{d}z,$$

(46)

where $J_n\left(z\right)$ is the Bessel function of the first kind; after standard mathematical treatment, we finally obtain

$$\begin{array}{ccc}\hfill T(x,y,t)& =& -\frac{2a{q}_{0}t}{\pi k}{\int }_0^{\infty }{E}_{\alpha ,2}\left(-a{\eta }^2{t}^{\alpha }\right)cos\left(y\eta \right)\mathrm{d}\eta \hfill \\ & +& \frac{a{q}_{0}t}{\pi k}{\int }_0^{\infty }\rho E_{\alpha ,2}\left(-a{\rho }^2{t}^{\alpha }\right)\hfill \\ & \times & {\int }_0^1\left\{(l+x)J_0\left[\rho \sqrt{y^2+{(l+x)}^2{u}^2}\right]+(l-x)J_0\left[\rho \sqrt{y^2+{(l-x)}^2{u}^2}\right]\right\}\mathrm{d}u\mathrm{d}\rho .\hfill \end{array}$$

(47)

Now, we present several special cases of the solution (48) associated with integer values of $\alpha$.

3.1. Localized Heat Conduction ($\alpha \to 0$)

In this case,

$$E_{0,2}(-x)=\frac{1}{1+x}.$$

(48)

Taking into account that [44,45]

$${\int }_0^{\infty }\frac{1}{x^2+c^2}cos\left(bx\right)\mathrm{d}x=\frac{\pi }{2c}{\mathrm{e}}^{-bc},b>0,c>0,$$

(49)

$${\int }_0^{\infty }\frac{x}{x^2+c^2}{J}_0\left(bx\right)\mathrm{d}x=K_0\left(bc\right),b>0,c>0,$$

(50)

where $K_n\left(z\right)$ is the modified Bessel function, we get

$$\begin{array}{ccc}\hfill T(x,y,t)& =& -\frac{q_0\sqrt{a}t}{k}exp\left(-\frac{y}{\sqrt{a}}\right)+\frac{q_{0}t}{k\pi }{\int }_0^1\{(l+x)K_0\left[\sqrt{\frac{y^2+{(l+x)}^2{u}^2}{a}}\right]\hfill \\ & +& (l-x)K_0\left[\sqrt{\frac{y^2+{(l-x)}^2{u}^2}{a}}\right]\}\mathrm{d}u.\hfill \end{array}$$

(51)

3.2. Classical Heat Conduction ($\alpha =1$)

For the standard heat conduction equation,

$$E_{1,2}(-x)=\frac{1-{\mathrm{e}}^{-x}}{x}.$$

(52)

Taking into consideration the following integrals evaluated by the authors,

$${\int }_0^{\infty }\frac{1-e^{-c{x}^2}}{x^2}cos\left(bx\right)\mathrm{d}x=\sqrt{\pi c}\left[exp\left(-\frac{b^2}{4c}\right)-\frac{\sqrt{\pi }b}{2\sqrt{c}}\mathrm{erfc}\left(\frac{b}{2\sqrt{c}}\right)\right],b>0,c>0,$$

(53)

$${\int }_0^{\infty }\frac{1-e^{-c{x}^2}}{x}{J}_0\left(bx\right)\mathrm{d}x=\frac{1}{2}{\mathbf{E}}_1\left(\frac{b^2}{4c}\right),b>0,c>0,$$

(54)

we arrive at

$$\begin{array}{ccc}\hfill T(x,y,t)& =& -\frac{2q_0\sqrt{at}}{k}\left[\frac{1}{\sqrt{\pi }}exp\left(-\frac{y^2}{4at}\right)-\frac{y}{2\sqrt{at}}\mathrm{erfc}\left(\frac{y}{2\sqrt{at}}\right)\right]\hfill \\ & +& \frac{q_0}{2\pi k}{\int }_0^1\left\{(l+x){\mathbf{E}}_1\left[\frac{y^2+{(l+x)}^2{u}^2}{4at}\right]+(l-x){\mathbf{E}}_1\left[\frac{y^2+{(l-x)}^2{u}^2}{4at}\right]\right\}\mathrm{d}u.\hfill \end{array}$$

(55)

Here, $\mathrm{erfc}\left(z\right)$ is the complementary error function, and ${\mathbf{E}}_n\left(z\right)$ is the exponential integral [46]

$${\mathbf{E}}_n\left(z\right)={\int }_1^{\infty }\frac{{\mathrm{e}}^{-zt}}{t^n}\mathrm{d}t.$$

(56)

We have changed the standard notation of the exponential integral ${\mathbf{E}}_n\left(z\right)$ to prevent confusion with the Mittag–Leffler function $E_{\alpha }\left(z\right)$.

3.3. Ballistic Heat Conduction ($\alpha =2$)

In the case of the wave equation for temperature,

$$E_{2,2}(-x)=\frac{sin\sqrt{x}}{\sqrt{x}},$$

(57)

the required integrals read

$${\int }_0^{\infty }\frac{sin\left(cx\right)cos\left(bx\right)}{x}\mathrm{d}x=\left\{\begin{array}{cc}{\displaystyle \pi /2,}\hfill & 0<b<c,\hfill \\ 0,\hfill & 0<c<b,\hfill \end{array}\right.$$

(58)

$${\int }_0^{\infty }sin\left(cx\right)J_0\left(bx\right)\mathrm{d}x=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\sqrt{c^2-b^2}},}\hfill & 0<b<c,\hfill \\ 0,\hfill & 0<c<b,\hfill \end{array}\right.$$

(59)

and using the Heaviside step function the solution takes the form

$$\begin{array}{ccc}\hfill T(x,y,t)& =& -\frac{\sqrt{a}{q}_0}{k}H\left(\sqrt{a}t-y\right)\hfill \\ & +& \frac{\sqrt{a}{q}_0}{\pi k}\underset{x-l}{\overset{x+l}{\int }}\left\{\begin{array}{cc}{\displaystyle \frac{1}{\sqrt{a{t}^2-y^2-u^2}},}\hfill & {\displaystyle u^2<a{t}^2-y^2}\hfill \\ 0,\hfill & {\displaystyle u^2>a{t}^2-y^2}\hfill \end{array}\right\}\mathrm{d}u.\hfill \end{array}$$

(60)

The dependence of the solution on the spatial coordinate x is shown in Figure 1 and Figure 2; the dependence of the solution on the spatial coordinate y is depicted in Figure 3 and Figure 4. In calculations, we have used the non-dimensional quantities

$$\overline{x}=\frac{x}{l},\overline{y}=\frac{y}{l},\kappa =\frac{\sqrt{a}{t}^{\alpha /2}}{l},\overline{T}=\frac{kl}{a{q}_{0}t}T.$$

(61)

4. Conclusions

We have examined the time-fractional heat conduction equation with the Caputo fractional derivative of the order $0<\alpha \le 2$ in an infinite plane with two external half-infinite cracks whose surfaces are exposed to the heat flux loading. The solution has been obtained using the integral transform technique. Several special cases corresponding to the integer values of the order of the time-derivative have been considered.