二维Caputo型时间分数阶热传导方程非齐次初边值问题的研究

Study on Nonhomogeneous Initial Boundary Value Problems of Two-Dimensional Caputo Type Time Fractional Heat Conduction Equation

摘要: 本文基于Povstenko型分数阶热弹性理论推导出Caputo型时间分数阶热传导方程，应用分离变量法和Laplace变换法，求解了二维有热源项的时间分数阶热传导方程的非齐次初边值问题，并给出详细的求解过程。最后通过具体的算例，分析了分数阶阶数 在不同的初边值条件下对于热传导过程的影响。

Abstract: In this paper, Caputo type time fractional heat conduction equation is derived based on Povstenko type fractional thermoelasticity theory. The nonhomogeneous initial boundary value problem of two-dimensional time fractional heat conduction equation with heat source term is solved by using variable separation method and Laplace transform method, and the detailed solution process is given. Finally, through a specific example, the influence of fractional order   on heat conduction process under different initial and boundary value conditions is analyzed.

1. 引言

随着分数阶微积分应用范围的增加，很多的专家学者渐渐地加入了分数阶微分方程理论与应用的研究中，无论对哪一方面，分数阶微分方程的求解都有重要意义，并且相较于数值解而言，方程的精确解能够全局地反映参数变化时系统的动力学性质的变化。因此，对于分数阶微分方程的精确解和求解方法的研究就变得十分重要。

随着社会技术的快速发展，人们提出了许多有效的关于分数阶微分方程的求解方法，例如分离变量法、格林函数法 [1]、Adomian分解法 [2]、同伦分析法 [3] 以及积分变换法等。Wyss [4] 在1986年研究了时间分数阶扩散方程，使用Mellin变换得到了含有Fox函数形式的解析解。2000年Gorenfo [5] 在求解时间分数阶扩散–波动方程的过程中，通过Laplace变换得到了含有Wright函数形式的尺度不变解。Agrawal [6] 应用分离变量方法求给出了Caputo型分数阶扩散方程在齐次Dirichlet边界条件下的Mittag-Leffler函数形式解析解。Huang和Liu [7] 将时间–空间分数阶扩散方程延伸到时间–空间分数阶对流弥散方程，再通过积分变换法得到了该方程的解析解。2012年，Atanacković [8] 基于时空Cattaneo热传导定律，推导出Caputo型空间分数阶热传导方程，并利用Laplace变换和Fourier变换得到了解的显式形式。2015年，Meilanov [9] 给出了具有扩散和对流传热机制的分数阶热传导方程的解，并研究探讨了分数阶导数对温度分布随时间和坐标变化的影响。2019年，Siedlecka [10] 对时间变量采用Laplace变换法，对空间变量使用特征函数级数展开法，求出了有限区域内Caputo时间分数阶单相滞后热传导问题的解析解。何松林 [11] 等人利用Laplace变换法求解出分数阶振子自由振动方程的含有Mittag-Leffler函数形式的解析解，并以此来研究和分析分数振子的运动规律和性质。2021年，Sylvain [12] 使用经典积分变换法求解出了稳态分数阶平流扩散方程的解析解，并据此来模拟空气污染物在有限介质中的扩散。

近些年来，人们发现使用了一些很特殊的材料如黏弹性材料、软物质材料、多孔材料等，而这些材料的热弹性行为与经典的热传导理论不同，所产生的反常传导、反常扩散需要用分数阶微分方程来描述。近些年来，学者们在经典的热弹性理论的基础上，将分数阶微分算子引入热传导、热弹性力学、黏弹性力学中进行修正，随后，不同类型的分数阶热弹性理论相继问世，如Povstenko [13] [14] 型分数阶热弹性理论、Youssef [15] 型分数阶热弹性理论以及Ezzat [16] 型分数阶热弹性理论。因此，在这些分数阶热弹性理论的基础上，如何求解分数阶热传导方程就成为一个重要的研究课题。本文基于Povstenko型分数阶热弹性理论，给出了Caputo型时间分数阶热传导方程具体的求解过程。

2. 预备知识

2.1. Caputo分数阶微分算子

函数 $f\left(t\right)$ 的 $\alpha$ 阶Caputo型分数阶导数定义为 [17] [18]

${}_0{}^{C}D{}_t^{\alpha }f\left(t\right)=\{\begin{array}{l}\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_0^t\frac{f^{\left(n\right)}\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha -n+1}}\text{d}\tau ,n-1<\alpha \le n,}\\ \frac{{\text{d}}^{n}f\left(t\right)}{\text{d}{t}^n},\alpha =n\in N.\end{array}$ (1)

对连续的函数 $f\left(t\right)$ 且导数 ${}_0{}^{C}D{}_t^{\lambda +\beta }f\left(t\right)$ 存在，则有

${}_0{}^{C}D{}_t^{\lambda }{}_0{}^{C}D{}_t^{-\beta }f\left(t\right)={}_0{}^{C}D{}_t^{\lambda -\beta }f\left(t\right)\left(\lambda >0,\beta >0\right)$ (2)

对函数 $f\left(t\right)$ 的Caputo型分数阶导数进行Laplace变换

$\mathcal{L}\left\{{}_0{}^{C}D{}_t^{\alpha }f\left(t\right);s\right\}=s^{\alpha }\mathcal{L}\left[f\left(t\right)\right]-{\displaystyle \underset{k=0}{\overset{n-1}{\sum }}{f}^{\left(k\right)}\left(0\right)s^{\alpha -k-1}},n-1<\alpha \le n$ (3)

2.2. Mittag-Leffler函数

单参数的Mittag-Leffler函数 [19]

$E_{\alpha }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(\alpha k+1\right)}}\left(\alpha >0\right)$ (4)

双参数的Mittag-Leffler函数

$E_{\alpha ,\beta }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(\alpha k+\beta \right)}}\left(\alpha >0,\beta >0\right)$ (5)

双参数Mittag-Leffler函数的Laplace变换 [20]

(6)

其中

$E_{\alpha ,\beta }^{\left(k\right)}\left(z\right)=\frac{{\text{d}}^k}{\text{d}{z}^k}{E}_{\alpha ,\beta }\left(z\right)={\displaystyle \underset{j=0}{\overset{\infty }{\sum }}\frac{\left(j+k\right)!z^j}{j!\Gamma \left(\alpha \left(j+k\right)+\beta \right)}}$

Mittag-Leffler函数的重要运算公式 [21]

${\displaystyle {\int }_0^x{t}^{\alpha -1}}{E}_{\alpha ,\alpha }\left(-a{t}^{\alpha }\right){\left(x-t\right)}^{\beta -1}{E}_{\alpha ,\beta }\left(-b{\left(x-t\right)}^{\alpha }\right)\text{d}t=\frac{E_{\alpha ,\beta }\left(-b{x}^{\alpha }\right)-E_{\alpha ,\beta }\left(-a{x}^{\alpha }\right)}{a-b}{x}^{\beta -1},a\ne b$ (7)

3. 二维Caputo型时间分数阶热传导方程的解析解

3.1. Caputo型时间分数阶热传导方程的导出

Povstenko提出的分数阶热弹性理论

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

观察Caputo型分数阶导数定义，两式作比较，我们可以得到时间分数阶Fourier定律

$q=-k\left({}_0{}^{C}D{}_t^{1-\alpha }\mathrm{grad}\text{}\text{}\text{}\text{}\text{}\text{}\text{}T\right),0<\alpha \le 1$ (9)

根据能量守恒，局部热通量为

$-\nabla \cdot q=\rho C\frac{\partial T}{\partial t}$ (10)

其中，C是恒压下的比热容， $\rho$ 是密度。

联立(9)和(10)式

$\rho C\frac{\partial T}{\partial t}=k{}_0{}^{C}D{}_t^{\alpha }\Delta T, 0<\alpha \le 1$ (11)

再根据Caputo型分数阶导数的性质，有

${}_0{}^{C}D{}_t^{\alpha }T=c^2\Delta T, 0<\alpha \le 1$ (12)

其中 $c^2=\frac{k}{\rho C}$， $\Delta$ 是Laplace算子，称上式为Caputo型时间分数阶热传导方程。

3.2. 二维非齐次分数阶热传导方程的齐次初边值问题

对于二维齐次边界条件的齐次时间分数阶热传导方程

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }u\left(x,y,t\right)=c^2\left[\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial y^2}\right],0<x\le a,0<y\le b,t>0,0<\alpha \le 1\\ u\left(x,y,0\right)=\phi \left(x,y\right),\text{}0<x\le a,0<y\le b\\ u\left(0,y,t\right)=0,u\left(a,y,t\right)=0,0<y\le b,t\ge 0\\ u\left(x,0,t\right)=0,u\left(x,b,t\right)=0,0<x\le a,t\ge 0\end{array}$ (13)

首先，我们应用分离变量法，设 $u\left(x,y,t\right)=U\left(x,y\right)T\left(t\right)$，带入上述方程可获得以下两个微分方程

$\frac{{\partial }^{2}U\left(x,y\right)}{\partial x^2}+\frac{{\partial }^{2}U\left(x,y\right)}{\partial y^2}+\lambda U\left(x,y\right)=0$ (14)

${}_0{}^{C}D{}_t^{\alpha }T\left(t\right)+c^2\lambda T\left(t\right)=0$ (15)

再设 $U\left(x,y\right)=X\left(x\right)Y\left(y\right)$，带入(14)式，分离边界条件，整理得到关于x和y的两个线性微分方程

$\{\begin{array}{l}{X}^{″}\left(x\right)+\mu X\left(x\right)=0\\ X\left(0\right)=X\left(a\right)=0\end{array}$ (16)

$\{\begin{array}{l}{Y}^{″}\left(y\right)+\left(\lambda -{\beta }^2\right)Y\left(y\right)=0\\ Y\left(0\right)=Y\left(b\right)=0\end{array}$ (17)

其中 $\lambda ,\mu$ 为分离常数，对于二阶齐次线性微分方程(16)，该问题具有特征值 ${\beta }_m$ 和特征函数 $X_m\left(x\right)$，且在 $L^2\left[0,L\right]$ 空间中满足下列条件 [22]

${\displaystyle {\int }_0^a{X}_m^2\left(x\right)\text{d}x}={\Vert X_m\left(x\right)\Vert }^2$

当分离常数 $\mu \ge 0$ 时，方程才有解，我们记 $\mu ={\beta }^2\left(\beta >0\right)$，它的通解为 $X\left(x\right)=C\sin \beta x+D\cos \beta x$，带入边界条件可求得到特征值 ${\beta }_m$ 以及相应的特征函数 $X_m(\; x\; )$

${\beta }_m=\frac{m\pi }{a},m=0,1,2,\cdots$ (18)

$X_m\left(x\right)=\sin \frac{m\pi }{a}x,m=1,2,\cdots$ (19)

再考虑微分方程(17)，当 $\lambda -{\beta }^2\le 0$ 时，方程没有非平凡解，记 ${\gamma }^2=\lambda -{\beta }^2\left(\lambda >0\right)$，得到特征值 ${\gamma }_n$ 以及对应的特征函数 $Y_n(\; y\; )$

${\gamma }_n=\frac{n\pi }{b},n=1,2,\cdots$ (20)

$Y_n\left(y\right)=\sin \frac{n\pi }{b}y,n=1,2,\cdots$ (21)

注意 $\lambda ={\gamma }^2+{\beta }^2$，且 $Y_n\left(y\right)$ 在 $L^2\left[0,L\right]$ 空间中满足下列条件

${\displaystyle {\int }_0^b{Y}_n^2\left(y\right)\text{d}y}={\Vert Y_n\left(y\right)\Vert }^2$

对于方程(15)，我们使用Laplace变换，有

(22)

其中 ${\lambda }_{mn}={\gamma }^2+{\beta }^2=\frac{m^2{\pi }^2}{a^2}+\frac{n^2{\pi }^2}{b^2}$，经过整理可得

$\mathcal{L}\left[T_{mn}\left(t\right)\right]={\displaystyle \underset{k=0}{\overset{j-1}{\sum }}{T}_{mn}^{\left(k\right)}\left(0\right)}\frac{s^{\alpha -k-1}}{s^{\alpha }+c^2{\lambda }_{mn}}$ (23)

其中 $T_{mn}^{\left(k\right)}\left(0\right)$ 是广义傅里叶系数，由初始条件 $u\left(x,y,0\right)=\phi \left(x,y\right)$ 确定，因 $0<\alpha \le 1$，上式中 $j=1$， $k=0$，我们得到

$T_{mn}\left(0\right)=\frac{1}{{\Vert X_n\left(x\right)\Vert }^{\text{2}}{\Vert Y_m\left(y\right)\Vert }^{\text{2}}}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^b\phi }\left(x,y\right)X_m\left(x\right)Y_n\left(y\right)\text{d}x}\text{d}y,m,n=1,2,\cdots$ (24)

再由Mittag-Leffler函数的Laplace变换公式得

$T_{mn}\left(t\right)=T_{mn}\left(0\right)E_{\alpha ,1}\left(-c^2{\lambda }_{mn}{t}^{\alpha }\right)$ (25)

因此，我们得到二维齐次边界条件下的齐次时间分数阶热传导方程(13)的解

$u\left(x,y,t\right)=X\left(x\right)Y\left(y\right)T\left(t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{T}_{mn}\left(0\right)E_{\alpha ,1}\left(-c^2{\lambda }_{mn}{t}^{\alpha }\right)}}\sin \frac{m\pi }{a}x\sin \frac{n\pi }{b}y$ (26)

其中 ${\lambda }_{mn}=\frac{m^2{\pi }^2}{a^2}+\frac{n^2{\pi }^2}{b^2}$， $T_{mn}\left(0\right)$ 由(24)式确定。

对于二维齐次边界条件的非齐次时间分数阶热传导方程

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }u\left(x,y,t\right)=c^2\Delta u\left(x,y,t\right)+F\left(x,y,t\right),0<x\le a,0<y\le b,t>0,0<\alpha \le 1\\ u\left(x,y,0\right)=\phi \left(x,y\right),\text{}0<x\le a,0<y\le b\\ u\left(0,y,t\right)=0,u\left(a,y,t\right)=0,0<y\le b,t\ge 0\\ u\left(x,0,t\right)=0,u\left(x,b,t\right)=0,0<x\le a,t\ge 0\end{array}$ (27)

它的解可以通过特征函数 $X_m\left(x\right)Y_n\left(y\right)$ 来表示，我们设为

$u\left(x,y,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{v}_{mn}\left(t\right)X_m}}\left(x\right)Y_n\left(y\right),m,n=1,2,\cdots$ (28)

其中是 $v_{mn}\left(t\right)$ 是关于t的待定函数，为了求出 $v_{mn}\left(t\right)$，我们以下面形式展开 $F\left(x,y,t\right)$，

$F\left(x,y,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{F}_{mn}\left(t\right)X_m}}\left(x\right)Y_n\left(y\right)$ (29)

其中 $F_{mn}\left(t\right)=\frac{4}{ab}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^{b}F\left(x,y,t\right)}\sin \frac{m\pi }{a}x\sin \frac{n\pi }{b}y\text{d}x}\text{d}y$，将(28)、(29)式代入方程(27)，比较两边的系数得

${}_0{}^{C}D{}_t^{\alpha }{v}_{mn}\left(t\right)+c^2{\lambda }_{mn}{v}_{mn}\left(t\right)=F_{mn}\left(t\right)$ (30)

其中 ${\lambda }_{mn}=\frac{m^2{\pi }^2}{a^2}+\frac{n^2{\pi }^2}{b^2}$，对上式两边使用Laplace变换，我们得到

(31)

其中 $v_{mn}^{\left(k\right)}\left(0\right)$ 是广义傅里叶系数，由初始条件 $u\left(x,y,0\right)=\phi \left(x,y\right)$ 确定，因 $0<\alpha \le 1$，上式中 $j=1$， $k=0$，我们有

$v_{mn}\left(0\right)=\frac{1}{{\Vert X_n\left(x\right)\Vert }^{\text{2}}{\Vert Y_m\left(y\right)\Vert }^{\text{2}}}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^b\phi \left(x,y\right)X_m}\left(x\right)Y_n\left(y\right)\text{d}x}\text{d}y$ (32)

则式(31)可写为

(33)

再由M-L函数的Laplace变换公式得

$v_{mn}\left(t\right)={\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-c^2{\lambda }_{mn}{\left(t-\tau \right)}^{\alpha }\right)F_{mn}\left(\tau \right)\text{d}\tau }+v_{mn}\left(0\right)E_{\alpha ,1}\left(-c^2{\lambda }_{mn}{t}^{\alpha }\right)$ (34)

于是，我们得到

$\begin{array}{c}u\left(x,y,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{v}_{mn}\left(t\right)X_m}}\left(x\right)Y_n\left(y\right)\\ ={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}[{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-c^2{\lambda }_{mn}{\left(t-\tau \right)}^{\alpha }\right)F_{mn}\left(\tau \right)\text{d}\tau }}}\\ +v_{mn}\left(0\right)E_{\alpha ,1}\left(-c^2{\lambda }_{mn}{t}^{\alpha }\right)]\sin \frac{m\pi }{a}x\sin \frac{n\pi }{b}y\end{array}$ (35)

其中 $v_{mn}\left(0\right)$ 由式(32)确定， ${\lambda }_{mn}=\frac{m^2{\pi }^2}{a^2}+\frac{n^2{\pi }^2}{b^2}$，这就是二维齐次边界条件的非齐次时间分数阶热传导(27)的解。

3.3. 二维非齐次分数阶热传导方程的非齐次初边值问题

对于二维非齐次边界条件的非齐次时间分数阶热传导方程

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }u\left(x,y,t\right)=c^2\Delta u\left(x,y,t\right)+F\left(x,y,t\right),0<x\le a,0<y\le b,t>0,0<\alpha \le 1\\ u\left(x,y,0\right)=\phi \left(x,y\right),\text{}0<x\le a,0<y\le b\\ u\left(0,y,t\right)=g_1\left(y,t\right),u\left(a,y,t\right)=g_2\left(y,t\right),0<y\le b,t\ge 0\\ u\left(x,0,t\right)=h_1\left(x,t\right),u\left(x,b,t\right)=h_2\left(x,t\right),0<x\le a,t\ge 0\end{array}$ (36)

首先，作变量替换 $u\left(x,y,t\right)=w\left(x,y,t\right)+v\left(x,y,t\right)$，并选择 $w\left(x,y,t\right)$ 满足齐次边界条件，我们可得到两个分数阶微分方程

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }w\left(x,y,t\right)=c^2\Delta w\left(x,y,t\right)+F\left(x,y,t\right)\\ w\left(x,y,0\right)=\phi \left(x,y\right)\\ w\left(0,y,t\right)=0,w\left(a,y,t\right)=0\\ w\left(x,0,t\right)=0,w\left(x,b,t\right)=0\end{array}$ (37)

及

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }v\left(x,y,t\right)=c^2\Delta v\left(x,y,t\right)\\ v\left(x,y,0\right)=0\\ v\left(0,y,t\right)=g_1\left(y,t\right),v\left(a,y,t\right)=g_2\left(y,t\right)\\ v\left(x,0,t\right)=h_1\left(x,t\right),v\left(x,b,t\right)=h_2\left(x,t\right)\end{array}$ (38)

对于方程(37)，其求解过程上一节已经给出，接下来求解方程(38)，设 $v\left(x,y,t\right)=K\left(x,y,t\right)+M\left(x,y,t\right)$，通过变量替换将非齐次边界条件齐次化，并选择 $M\left(x,y,t\right)$ 在x方向满足边界条件

$M\left(0,y,t\right)=g_1\left(y,t\right),M\left(a,y,t\right)=g_2\left(y,t\right)$

可得到

$M\left(x,y,t\right)=\frac{a-x}{a}{g}_1\left(y,t\right)+\frac{x}{a}{g}_2\left(y,t\right)$ (39)

将 $v\left(x,y,t\right)$ 带入方程(38)，整理得

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }K\left(x,y,t\right)=c^2\Delta K\left(x,y,t\right)+\tilde{F}\left(x,y,t\right)\\ K\left(x,y,0\right)=\tilde{\phi }\left(x,y\right)\\ K\left(0,y,t\right)=0,K\left(a,y,t\right)=0\\ K\left(x,0,t\right)={\tilde{h}}_1\left(x,t\right),K\left(x,b,t\right)={\tilde{h}}_2\left(x,t\right)\end{array}$ (40)

其中

$\{\begin{array}{l}\tilde{F}\left(x,y,t\right)=c^2\frac{{\partial }^{2}M\left(x,y,t\right)}{\partial y^2}-{}_0{}^{C}D{}_t^{\alpha }M\left(x,y,t\right)\\ \tilde{\phi }\left(x,y\right)=-\frac{a-x}{a}{g}_1\left(y,0\right)-\frac{x}{a}{g}_2\left(y,0\right)\\ {\tilde{h}}_1\left(x,t\right)=h_1\left(x,t\right)-\frac{a-x}{a}{g}_1\left(0,t\right)-\frac{x}{a}{g}_2\left(0,t\right)\\ {\tilde{h}}_2\left(x,t\right)=h_2\left(x,t\right)-\frac{a-x}{a}{g}_1\left(b,t\right)-\frac{x}{a}{g}_2\left(b,t\right)\end{array}$

因 $K\left(x,y,t\right)$ 在x方向满足齐次边界条件，它的解可以通过特征函数 $X_m\left(x\right)=\sin \frac{m\pi }{a}$ 来表示，设为

$K\left(x,y,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{K}_m\left(y,t\right)X_m\left(x\right)}$ (41)

其中 $K_m\left(y,t\right)=\frac{2}{a}{\displaystyle {\int }_0^{a}K\left(x,y,t\right)\sin \frac{m\pi }{a}x\text{d}x}$，同样也可以通过特征函数 $X_m\left(x\right)$ 来表示 $\tilde{F}\left(x,y,t\right)$，即

$\tilde{F}\left(x,y,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{F}_m\left(y,t\right)X_m\left(x\right)}$ (42)

其中 $F_m\left(y,t\right)=\frac{2}{a}{\displaystyle {\int }_0^a\tilde{F}\left(x,y,t\right)\sin \frac{m\pi }{a}x\text{d}x}$，将(41)式和(42)式带入方程(40)整理得到

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }{K}_m\left(y,t\right)=c^2\left[-{\left(\frac{m\pi }{a}\right)}^2{K}_m\left(y,t\right)+\frac{{\partial }^2{K}_m\left(y,t\right)}{\partial y^2}\right]+F_m\left(y,t\right)\\ {K_m\left(y,t\right)|}_{t=0}={\phi }_m\left(y\right)\\ {K_m\left(y,t\right)|}_{y=0}=h_{m1}\left(t\right),{K_m\left(y,t\right)|}_{y=b}=h_{m2}\left(t\right)\end{array}$ (43)

其中

$\{\begin{array}{l}{\phi }_m\left(y\right)=\frac{2}{a}{\displaystyle {\int }_0^a\tilde{\phi }\left(x,y\right)\sin \frac{m\pi }{a}x\text{d}x}\\ h_{m1}\left(t\right)=\frac{2}{a}{\displaystyle {\int }_0^a{\tilde{h}}_1\left(x,t\right)\sin \frac{m\pi }{a}x\text{d}x},h_{m2}\left(t\right)=\frac{2}{a}{\displaystyle {\int }_0^a{\tilde{h}}_2\left(x,t\right)\sin \frac{m\pi }{a}x\text{d}x}\end{array}$

可以看出方程(43)是关于 $y,t$ 非齐次时间分数阶热传导方程非齐次初边值问题，设

$K_m\left(y,t\right)=P_m\left(y,t\right)+Q_m\left(y,t\right)$ (44)

并选择 $Q_m\left(y,t\right)$ 满足边界条件

${Q_m\left(y,t\right)|}_{y=0}=h_{m1}\left(t\right),{Q_m\left(y,t\right)|}_{y=b}=h_{m2}(\; t\; )$

容易得到

$Q_m\left(y,t\right)=\frac{b-y}{b}{h}_{m1}\left(t\right)+\frac{y}{b}{h}_{m2}\left(t\right)$ (45)

将 $K_m\left(y,t\right)=P_m\left(y,t\right)+Q_m\left(y,t\right)$ 带入方程(43)，得到

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }{P}_m\left(y,t\right)=c^2\frac{{\partial }^2{P}_m\left(y,t\right)}{\partial y^2}-c^2{\left(\frac{m\pi }{a}\right)}^2{P}_m\left(y,t\right)+{\tilde{F}}_m\left(y,t\right)\\ {P_m\left(y,t\right)|}_{t=0}={\tilde{\phi }}_m\left(y\right)\\ {P_m\left(y,t\right)|}_{y=0}=0,{P_m\left(y,t\right)|}_{y=b}=0\end{array}$ (46)

其中

$\{\begin{array}{l}{\tilde{F}}_m\left(y,t\right)=F_m\left(y,t\right)-c^2{\left(\frac{m\pi }{a}\right)}^2{Q}_m\left(y,t\right)-{}_0{}^{C}D{}_t^{\alpha }{Q}_m\left(y,t\right)\\ {\tilde{\phi }}_m\left(y\right)={\phi }_m\left(y\right)-\frac{b-y}{b}{h}_{m1}\left(0\right)+\frac{y}{b}{h}_{m2}(\; 0\; )\end{array}$

方程(46)的解可以通过特征函数 $Y_n\left(x\right)=\sin \frac{n\pi }{b}$ 来表示，设为

$P_m\left(y,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{T}_{mn}\left(t\right)Y_n\left(y\right)}$ (47)

${\tilde{F}}_m\left(y,t\right)$ 同样可以用下面形式展开，即

${\tilde{F}}_m\left(y,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\tilde{F}}_{mn}\left(t\right)Y_n\left(y\right)}$ (48)

其中 ${\tilde{F}}_{mn}\left(t\right)=\frac{2}{b}{\displaystyle {\int }_0^b{\tilde{F}}_m\left(y,t\right)Y_n\left(y\right)\text{d}y}$，将(47)式和(48)式带入方程(46)整理可得

${}_0{}^{C}D{}_t^{\alpha }{T}_{mn}\left(t\right)+c^2\left[{\left(\frac{m\pi }{a}\right)}^2+{\left(\frac{n\pi }{b}\right)}^2\right]T_{mn}\left(t\right)={\tilde{F}}_{mn}\left(t\right)$ (49)

对上式两边使用Laplace变换，

(50)

记 ${\lambda }_{mn}=\frac{m^2{\pi }^2}{a^2}+\frac{n^2{\pi }^2}{b^2}$， $T_{mn}\left(0\right)$ 由初始条件 ${P_m\left(y,t\right)|}_{t=0}={\tilde{\phi }}_m\left(y\right)$ 推出，因 $0<\alpha \le 1$，上式中 $j=1$， $k=0$，我们得到

$T_{mn}\left(0\right)=\frac{2}{b}{\displaystyle {\int }_0^b{\tilde{\phi }}_m\left(y\right)\sin \frac{n\pi y}{b}\text{d}y}$ (51)

则式(50)可写为

(52)

再由Mittag-Leffler函数的Laplace变换公式得

$T_{mn}\left(t\right)={\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-c^2{\lambda }_{mn}{\left(t-\tau \right)}^{\alpha }\right){\tilde{F}}_{mn}\left(\tau \right)\text{d}\tau }+T_{mn}\left(0\right)E_{\alpha ,1}\left(-c^2{\lambda }_{mn}{t}^{\alpha }\right)$ (53)

从而得到方程(46)的解

$\begin{array}{c}{P}_m\left(y,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{T}_{mn}\left(t\right)Y_n\left(y\right)}\\ ={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}[{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-c^2{\lambda }_{mn}{\left(t-\tau \right)}^{\alpha }\right){\tilde{F}}_{mn}\left(\tau \right)\text{d}\tau }}\\ +T_{mn}\left(0\right)E_{\alpha ,1}\left(-c^2{\lambda }_{mn}{t}^{\alpha }\right)]\sin \frac{n\pi y}{b}\end{array}$ (54)

综上，我们得到方程(36)的解

$\begin{array}{c}u\left(x,y,t\right)=w\left(x,y,t\right)+v\left(x,y,t\right)\\ =w\left(x,y,t\right)+K\left(x,y,t\right)+M\left(x,y,t\right)\\ =w\left(x,y,t\right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{K}_m\left(y,t\right)X_m\left(x\right)}+M\left(x,y,t\right)\\ =w\left(x,y,t\right)+{\displaystyle \underset{m=1}{\overset{\infty }{\sum }}\left[P_m\left(y,t\right)+Q_m\left(y,t\right)\right]X_m\left(x\right)}+M\left(x,y,t\right)\end{array}$ (55)

其中 $w\left(x,y,t\right)$ 由(35)式给出， $M\left(x,y,t\right)$ 由(39)式给出， $P_m\left(y,t\right)$ 由(54)式给出， $Q_m\left(y,t\right)$ 由(45)式给出。

4. 数值算例

在二维空间中，考虑一块板状物体，侧面绝缘，且温度自由地扩散到温度为零的介质中去，并假定初始温度分布为 $\phi \left(x,y\right)$，在没有热源，且考虑前一温度历史影响的情况下，观察杆内的温度分布以及变化规律。它可用以下模型描述

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }u\left(x,y,t\right)=c^2\left[\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial y^2}\right],0<x\le a,0<y\le b,t>0,0<\alpha \le 1\\ u\left(x,y,0\right)=\phi \left(x,y\right)=xy,\text{}0<x\le a,0<y\le b\\ u\left(0,y,t\right)=0,u\left(a,y,t\right)=0,0<y\le b,t\ge 0\\ u\left(x,0,t\right)=0,u\left(x,b,t\right)=0,0<x\le a,t\ge 0\end{array}$ (56)

我们取 $\phi \left(x,y\right)=xy$，得到的解析解为

$u\left(x,y,t\right)={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\frac{4ab{\left(-1\right)}^{m+n+2}}{mn{\pi }^2}{E}_{\alpha ,1}\left[-c^2\left(\frac{m^2{\pi }^2}{a^2}+\frac{n^2{\pi }^2}{b^2}\right)t^{\alpha }\right]}}\sin \frac{m\pi }{a}x\sin \frac{n\pi }{b}y$ (57)

注意到 $E_{1,1}\left(z\right)={\text{e}}^z$，当 $\alpha =1$ 时，上式还可以写成下形式

$u\left(x,y,t\right)={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\frac{4ab{\left(-1\right)}^{m+n+2}}{mn{\pi }^2}{\text{e}}^{-c^2\left(\frac{m^2{\pi }^2}{a^2}+\frac{n^2{\pi }^2}{b^2}\right)t}}}\sin \frac{m\pi }{a}x\sin \frac{n\pi }{b}y$ (58)

这正是标准二维整数阶热传导方程的解析解，且当 $\alpha =1$ 时，选取初始条件与参考文献 [23] 中一致时，解析解的表达式也是一致的。令 $c=0.5$， $a=1$， $b=1$，在(57)式中取 $\alpha =1$，在同一坐标轴中画出(57)和(58)的图像，如下

Figure 1. When $t=0.07$, the temperature distribution in the rod

图1. 在 $t=0.07$ 时，杆内温度的分布情况

在函数值处，分数阶解析解(57)是用直线连接的，整数阶解析解(58)是用圆圈表示，从图1可以看出两图形完全重合，这说明本章的求解过程与求出的解析解是正确的。

图2给出了当 $\alpha =1$ 时，杆内温度的变化情况，可以看出两个端点位置温度恒定，在靠近 $x=1$ 端点处的温度下降的较快，在靠近 $x=0$ 端点处有一小段下降的较为平缓，但总体是下降的。

图3给出了当 $x=0.5,y=0.7$ 时，杆内温度随时间t的变化情况；当 $\alpha \in \left(0,1\right]$ 时，在 $t\le 0.65$ 时刻之前，杆内温度随着 $\alpha$ 的增加下降的越慢，在 $t>0.65$ 时刻之后，温度随着 $\alpha$ 的增加下降的越快；随着时间的推移， $\alpha$ 越大，前一时刻温度对后一时刻温度变化情况的影响越小。

Figure 2. When $y=0.2$, the temperature distribution in the rod

图2. 在 $y=0.2$ 时，杆内温度的分布情况

Figure 3. When $x=0.5,y=0.7$, effect of $\alpha$ on the temperature in the rod at different times

图3. 在 $x=0.5,y=0.7$ 时， $\alpha$ 对不同时刻杆内温度的影响

图4给出了温度随位移x的变化情况。可以看出，随着x的增加，各点温度都是缓慢增大后快速减小，且随着 $\alpha$ 的增大，杆内温度变化幅度的越大；这也说明， $\alpha$ 越小，前一时刻的记忆性越强。

5. 结论

本文在对分数阶热传导方程求解方法的研究中，首先基于Povstenko型分数阶热弹性理论，推导出Caputo型时间分数阶热传导方程，应用分离变量法求解了有热源项的时间分数阶热传导方程的各种初边值问题。其求解思路是：

1) 先构造合适的辅助函数，对方程中的两个方向的非齐次边界条件逐步进行齐次化处理，使原方程转化为齐次边界条件下的定解问题；

Figure 4.When $y=0.5,t=0.05$, effect of $\alpha$ on the temperature in the rod

图4. 在 $y=0.5,t=0.05$ 时， $\alpha$ 对杆内温度的影响

2) 再根据分离变量法将方程分解成关于t的分数阶微分方程和关于x和y的整数阶微分方程；对于整数阶微分方程，依据相应的齐次边界条件求出特征值与特征函数；

3) 对于分数阶微分方程，若方程中存在热源项，将其解假设为变系数的含特征函数形式，再依据初始条件由傅里叶变换求出变系数，最后通过Laplace变换求出含有M-L函数形式的解析解。

最后通过具体的算例分析了不同 $\alpha$ 值对热传导过程的影响，得出结论：随着时间的推移以及 $\alpha$ 的增大，前一时刻温度对后一时刻温度变化情况的影响越小，这也说明， $\alpha$ 越小，前一时刻的记忆性越强。

基金项目

国家自然科学基金(11862014)。