关于t^a型热传导方程解的定性研究及数值模拟

Qualitative Study and Numerical Simulation of the Solution of t^a Heat Conduction Equation

摘要: 本文主要研究t^a型热传导方程。从t^a型热传导方程的导出入手，运用分离变量法求解此方程，得到t^a型热传导方程的通解。其次，利用t^a型热传导方程的边界条件以及初值条件，对其解的存在性、唯一性和稳定性进行了研究。最后，利用Matlab编程和有限差分法对热传导方程的解进行数值模拟，并讨论时间的幂次对热传导方程数值解的影响。

Abstract: This paper is to study the heat conduction equation related to t^a . First, we start from the origin of the heat conduction equation. The general solution of the heat conduction equation is obtained by using the separation of variables method. Secondly, we study the existence, uniqueness and stabil-ity of the solution of  t^a heat conduction equation under three boundary conditions. Finally, Matlab programming and finite difference method are used to simulate the solution of the heat conduction equation. And the influence of the power of time and the numerical solution of the heat conduction equation is discussed.

1. 引言

人类的生存发展与热能的传递过程息息相关。人们通过观察温度的差异所引起的热能的传递，对其规律性进行探索和总结，从而有了传热学。而热传导方程 [1] 是在人们研究热传导这一类自然现象的基本规律后产生的。人们将热传导的基本规律写成偏微分方程 [2] 的形式，并且在实践过程中来验证和发展这些规律。

20世纪以来，国际上面临一系列问题，能源危机、资源短缺促使人民重新认识能量流动，如何节约传热过程中的能量损耗问题成为了热点研究方向，极大地加强了传热学与不同领域之间的结合。另外，由于20世纪计算机技术的飞速发展，运用数值模拟的方法 [3] 对传热学中的微分方程进行数值求解成为了可能。学术界进一步地在这个基础上研究传热问题和发展新兴领域，也出现了传热学上新的分支——数值传热学。至今，由于微机电系统和纳机电系统研究的迅速发展，在世界范围内都兴起了研究微米–纳米尺度范围的传热和流动问题的热潮，形成了微米–纳米传热学的研究方向。

热能在物质中的传递有三种方式：热传导、热对流与热辐射。在本文中，我们仅仅只研究热能传递的其中一种——热传导过程 [4] [5] [6] 。在研究热传导方程之前，我们首先要理解热传导是指物体保持固有的位置，不发生任何相对运动，仅仅依靠分子、原子及自由电子等微观粒子的热运动而产生的热能传递。在热传导、扩散等物理现象中，我们将会遇到与热传导方程相似的这一类偏微分方程 [7] 。

标准的热传导方程中将热传导系数看作一个常量，对热传导方程的初边值问题 [8] ，利用分离变量的方法求得方程的通解并对热传导方程解的一系列性质定理进行了论述以及证明 [9] [10] [11] [12] 。但是在实际生产生活中热传导系数与温度等因素有关，那么在热量传递的过程中热传导参数会随着时间而变化。因此，本文主要就是研究热传导方程中的热传导系数与时间的幂次有关，对此时的热传导方程解的存在性、唯一性和稳定性进行研究并对此热传导方程的解进行数值模拟。

2. $t^{\alpha }$ 型热传导方程

2.1. $t^{\alpha }$ 型热传导方程的导出

在热传导问题中，我们以函数 $u\left(x,y,z,t\right)$ 表示空间中物体G在位置 $\left(x,y,z\right)$ 时刻t的温度。热能在物质中的传递过程的傅里叶实验定律表明，物体在无穷小时段 $\text{d}t$ 内沿着法线方向 $n$ 流过一个无穷小面积 $\text{d}S$

的热量 $\text{d}Q$ 与物体温度沿曲面 $\text{d}S$ 法线方向的方向导数 $\frac{\partial u}{\partial n}$ 成正比，即

$\text{d}Q=-k\left(x,y,z\right)\frac{\partial u}{\partial n}\text{d}S\text{d}t,$ (1)

其中 $k\left(x,y,z\right)$ 称为物体在点 $\left(x,y,z\right)$ 处的热传导系数。(1)式中的负号的出现是由于热力学第一定律告诉我们热量总是从温度高的一侧流向低的一侧，因此， $\text{d}Q$ 因和 $\frac{\partial u}{\partial n}$ 异号。

在物体G内任取一闭曲面 $\Gamma$ ，所包围的区域记为 $\Omega$ ，由(1)式，从时刻 $t_1$ 到 $t_2$ 流进此闭曲面的全部热量为

$Q={\displaystyle {\int }_{t_1}^{t_2}\left[{\displaystyle \underset{\Gamma }{\iint }k\left(x,y,z\right)\frac{\partial u}{\partial n}\text{d}S}\right]\text{d}t},$

流入的热量使物体内部温度发生变化，在时间间隔 $\left(t_1,t_2\right)$ 中物体温度从 $u\left(x,,y,z,t_1\right)$ 变化到 $u\left(x,,y,z,t_2\right)$ ，应该要吸收的热量是

${\displaystyle \underset{\Omega }{\iiint }c\left(x,y,z\right)\rho \left(x,y,z\right)\left[u\left(x,y,z,t_1\right)-u\left(x,y,z,t_2\right)\right]\text{d}x\text{d}y\text{d}z},$

其中c为比热， $\rho$ 为密度。于是由“热量守恒”有

${\displaystyle {\int }_{t_1}^{t_2}\left[{\displaystyle \underset{\Gamma }{\iint }k\left(x,y,z\right)\frac{\partial u}{\partial n}\text{d}S}\right]\text{d}t}={\displaystyle \underset{\Omega }{\iiint }c\rho \left[u\left(x,y,z,t_1\right)-u\left(x,y,z,t_2\right)\right]\text{d}x\text{d}y\text{d}z},$ (2)

利用格林公式可以将(2)式化为

${\displaystyle {\int }_{t_1}^{t_2}{\displaystyle \underset{\Omega }{\iiint }\left[\frac{\partial }{\partial x}\left(k\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial u}{\partial y}\right)+\frac{\partial }{\partial z}\left(k\frac{\partial u}{\partial z}\right)\right]\text{d}x\text{d}y\text{d}z}}={\displaystyle \underset{\Omega }{\iiint }c\rho \left({\displaystyle {\int }_{t_1}^{t_2}\frac{\partial u}{\partial t}\text{d}t}\right)\text{d}x\text{d}y\text{d}z},$

通过交换函数的积分性质，并且由于 $t_1,t_2$ 和区域 $\Omega$ 是任意的，我们得到了

$c\rho \frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(k\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial u}{\partial y}\right)+\frac{\partial }{\partial z}\left(k\frac{\partial u}{\partial z}\right).$

在实际的生活中，热传导系数是表征材料导热性能优劣的参数，不同材料的热传导系数不同，对于同一种材料，热传导系数还会与温度等系数有关。在本文中，我们将研究

$\frac{k\left(t\right)}{c\rho }\propto t^{\alpha },\left(0\le \alpha \le 1\right)$

在这样的情形下，即得

$\frac{\partial u}{\partial t}=t^{\alpha }\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right),\left(0\le \alpha \le 1\right)$ (3)

(3)式就是本文所说的 $t^{\alpha }$ 型热传导方程。

2.2. 一维 $t^{\alpha }$ 型热传导方程的分离变量法

对于一维 $t^{\alpha }$ 型热传导方程的初边值问题，以下我们将第一类边界条件的初边值问题为例详细的讨论对于系数与时间有关的热传导方程的求解过程。我们用分离变量法解如下的初边值问题：

$\{\begin{array}{l}{u}_t-t^{\alpha }{u}_{xx}=0,(4)\\ u\left(x,0\right)=\phi \left(x\right),\left(t>0,0<x<l,0\le \alpha \le 1\right),(5)\\ u\left(0,t\right)=0,u\left(l,t\right)=0,(6)\end{array}$

令

$u\left(x,t\right)=X\left(x\right)T\left(t\right),$

这里的 $X\left(x\right)$ 和 $T\left(t\right)$ 分别表示仅与x有关和仅与t有关的函数。把它代入方程(4)，得到

$X{T}^{\prime }=t^{\alpha }{X}^{″}T,$

即

$\frac{T^{\prime }}{t^{\alpha }T}=\frac{X^{″}}{X}.$

这等式只有在两边均等于常数时才成立。令此常数为 $-\lambda$ ，则有

$T^{\prime }+\lambda t^{\alpha }T=0$ (7)

$X^{″}+\lambda X=0.$ (8)

将 $u\left(x,t\right)=X\left(x\right)T\left(t\right)$ 代入方程(6)，得到

$X\left(0\right)=0,X\left(l\right)=0.$ (9)

我们先对 $\lambda$ 的取值范围进行讨论。

1) 当 $\lambda <0$ 时，此时方程(8)的通解为

$X\left(x\right)=A{\text{e}}^{\sqrt{-\lambda }x}+B{\text{e}}^{-\sqrt{-\lambda }x},$

由边界条件(9)可以得到 $A=B=0$ ，所以 $X\left(x\right)\equiv 0$ 。

2) 当 $\lambda =0$ 时，方程(8)的解为

$X\left(x\right)=Ax+B,$

由边界条件(9)易得 $A=B=0$ ，因而 $X\left(x\right)\equiv 0$ 。

3) 当 $\lambda >0$ 时，方程(8)的通解具有如下形式

$X\left(x\right)=A\cos \sqrt{\lambda }x+B\sin \sqrt{\lambda }x,$

由边界条件 $X\left(0\right)=0$ 得 $A=0$ 。于是由边界条件 $X\left(l\right)=0$ 可知

$X\left(l\right)=B\sin \sqrt{\lambda }l=0,$

因此，由式(8)和式(9)所构成的本征值问题的本征值为

$\lambda ={\left(\frac{n\pi }{l}\right)}^2\left(n=1,2,\cdots \right),$ (10)

其相应的本征函数即方程(8)的解为

$X_n\left(x\right)=B_n\sin \frac{n\pi }{l}x\left(n=1,2,\cdots \right).$

将式(10)代入方程(7)，可以求得

$T_n\left(t\right)=C_n{\text{e}}^{-\frac{{\left(n\pi \right)}^2}{l^2\left(\alpha +1\right)}{t}^{\alpha +1}}\left(n=1,2,\cdots \right),$

于是我们可以得到一列可分离变量的特解

$u_n\left(x,t\right)=D_n{\text{e}}^{-\frac{{\left(n\pi \right)}^2}{l^2\left(\alpha +1\right)}{t}^{\alpha +1}}\sin \frac{n\pi }{l}x\left(n=1,2,\cdots \right).$

由于方程(6)及边界条件(7)和(8)都是线性的，故可利用叠加原理构造级数形式的解

$\begin{array}{l}u\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{D}_n}{\text{e}}^{-\frac{{\left(n\pi \right)}^2}{l^2\left(\alpha +1\right)}{t}^{\alpha +1}}\sin \frac{n\pi x}{l},\\ D_n=\frac{2}{l}{\displaystyle {\int }_0^l\phi \left(\xi \right)\sin \frac{n\pi \xi }{l}\text{d}\xi }\left(n=1,2,\cdots \right).\end{array}$

3. $t^{\alpha }$ 型热传导方程解的性质

3.1. 解的存在性

首先我们用傅里叶变换来求解 $t^{\alpha }$ 型热传导方程的柯西问题：

$\{\begin{array}{l}\frac{\partial u}{\partial t}=t^{\alpha }\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right)\\ u\left(x,0\right)=\phi \left(x\right)\end{array}\left(0\le \alpha \le l\right)$ (11)

在研究这一个非齐次的 $t^{\alpha }$ 型热传导方程之前，我们先求解齐次 $t^{\alpha }$ 型热传导方程：

$\{\begin{array}{l}\frac{\partial u}{\partial t}=t^{\alpha }\frac{{\partial }^{2}u}{\partial x^2}\\ u\left(x,0\right)=\phi \left(x\right)\end{array}\left(0\le \alpha \le l\right)$ (12)

对x进行傅里叶变换，记 $F\left[u\left(x,t\right)\right]=\tilde{u}\left(\lambda ,t\right)$ ， $F\left[\phi \left(x\right)\right]=\tilde{\phi }\left(\lambda \right)$ 。

对式(10)的两边关于x进行傅里叶变换，可以得到

$\{\begin{array}{l}\frac{\text{d}\tilde{u}}{\text{d}t}=-t^{\alpha }{\lambda }^2\tilde{u}\\ \tilde{u}\left(\lambda ,0\right)=\tilde{\phi }(\lambda )\end{array}$

上式的解为

$\tilde{u}\left(\lambda ,t\right)=\tilde{\phi }\left(\lambda \right){\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}.$

函数 ${\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}$ 的傅里叶逆变换为

$\begin{array}{c}{F}^{-1}\left({\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}\right)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}}{\text{e}}^{i\lambda x}\text{d}\lambda \\ =\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}}\left(\cos \lambda x+i\sin \lambda x\right)\text{d}\lambda \\ =\frac{1}{\pi }{\displaystyle {\int }_0^{\infty }{\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}\cos \lambda x\text{d}\lambda }.\end{array}$

利用复变函数的积分计算得

${\displaystyle {\int }_0^{\infty }{\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}\cos \lambda x\text{d}\lambda }=\sqrt{\frac{\alpha +1}{t^{\alpha +1}}}\frac{\sqrt{\pi }}{2}{\text{e}}^{-\frac{\alpha +1}{4t^{\alpha +1}}{x}^2}$

所以

$F^{-1}\left[{\text{e}}^{-\frac{{\lambda }^2}{\alpha +1}{t}^{\alpha +1}}\right]=\frac{\sqrt{\alpha +1}}{2\sqrt{\pi t^{\alpha +1}}}{\text{e}}^{-\frac{\alpha +1}{4t^{\alpha +1}}{x}^2}.$

因此，利用傅里叶变换性质，我们可以得到齐次 $t^{\alpha }$ 型热传导方程的柯西问题的解为

$u\left(x,t\right)=\frac{\sqrt{\alpha +1}}{2\sqrt{\pi t^{\alpha +1}}}{\displaystyle {\int }_{-\infty }^{\infty }\phi \left(\xi \right){\text{e}}^{-\frac{{\left(x-\xi \right)}^2\left(\alpha +1\right)}{4t^{\alpha +1}}}}\text{d}\xi .$ (13)

在求解出齐次 $t^{\alpha }$ 型热传导方程的柯西问题后，由齐次化原理，非齐次方程柯西问题的解也可以写为

$u\left(x,t\right)={\displaystyle {\int }_0^{l}w\left(x,t;\tau \right)\text{d}\tau },$

其中， $w=w\left(x,t;\tau \right)$ 为下面的定解问题的解：

$\{\begin{array}{l}\frac{\partial w}{\partial t}=t^{\alpha }\frac{{\partial }^{2}w}{\partial x^2},t>\tau ,\\ w\left(x,\tau \right)=f\left(x,\tau \right).\end{array}$

这样，由叠加原理及上述求得的解，我们就求得了问题(11)的解为

$u\left(x,t\right)=\frac{\sqrt{\alpha +1}}{2\sqrt{\pi t^{\alpha +1}}}{\displaystyle {\int }_{-\infty }^{\infty }\phi \left(\xi \right){\text{e}}^{-\frac{{\left(x-\xi \right)}^2\left(\alpha +1\right)}{4t^{\alpha +1}}}}\text{d}\xi +\frac{\sqrt{\alpha +1}}{2\sqrt{\pi }}{\displaystyle {\int }_0^t{\displaystyle {\int }_{-\infty }^{\infty }\frac{f\left(\xi ,\tau \right)}{\sqrt{{\left(t-\tau \right)}^{\alpha +1}}}{\text{e}}^{-\frac{{\left(x-\xi \right)}^2\left(\alpha +1\right)}{4{\left(t-\tau \right)}^{\alpha +1}}}}\text{d}\xi }\text{d}\tau .$ (14)

定理1若 $\phi \left(x\right)$ 连续且有界，则 $t^{\alpha }$ 型热传导方程柯西问题的解是存在的。

证明：设 $\left|\phi \left(x\right)\right|\le M$ ，则注意到 ${\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{-{\zeta }^2}}\text{d}\zeta =\sqrt{\pi }$ ，由式(13)得

$\left|u\left(x,t\right)\right|\le M\frac{1}{\sqrt{\pi }}\left|{\displaystyle {\int }_{-\infty }^{\infty }\frac{\sqrt{\alpha +1}}{2\sqrt{t^{\alpha +1}}}{\text{e}}^{-\frac{{\left(x-\xi \right)}^2\left(\alpha +1\right)}{4t^{\alpha +1}}}}\text{d}\xi \right|=M\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{-\zeta }}\text{d}\zeta =M$

这就说明了式(13)是有界的。

下面证明式(13)满足方程(12)，在式(11)中，令 $\zeta =\frac{\left(\xi -x\right)\sqrt{\alpha +1}}{2\sqrt{t^{\alpha +1}}}$ ，得到

$u\left(x,t\right)=\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\phi \left(x+\frac{2\zeta \sqrt{t^{\alpha +1}}}{\sqrt{\alpha +1}}\right){\text{e}}^{-{\zeta }^2}\text{d}\zeta }.$

又有

$\phi \left(x_0\right)=\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\phi \left(x_0\right){\text{e}}^{-{\zeta }^2}\text{d}\zeta },$

因此

$u\left(x,t\right)-\phi \left(x_0\right)=\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\left[\phi \left(x+\frac{2\zeta \sqrt{t^{\alpha +1}}}{\sqrt{\alpha +1}}\right)-\phi \left(x_0\right)\right]{\text{e}}^{-{\zeta }^2}\text{d}\zeta }.$

对于所给的 $\epsilon >0$ ，取 $N>0$ 足够大，使

$\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_N^{\infty }{\text{e}}^{-{\zeta }^2}}\text{d}\zeta \le \frac{\epsilon }{6M},\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{-N}{\text{e}}^{-{\zeta }^2}}\text{d}\zeta \le \frac{\epsilon }{6M},$

固定N，由 $\phi \left(x\right)$ 的连续性，可找到 $\delta >0$ ，使当 $\left|x-x_0\right|\le \delta ,0<t\le \delta$ 时成立

$\left|\phi \left(x+\frac{2\zeta \sqrt{t^{\alpha +1}}}{\sqrt{\alpha +1}}\right)-\phi \left(x_0\right)\right|\le \frac{\epsilon }{3}\left(-N\le \zeta \le N\right)$

因此

$\begin{array}{c}\left|u\left(x,t\right)-\phi \left(x_0\right)\right|=\left|\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{\infty }\left[\phi \left(x+\frac{2\zeta \sqrt{t^{\alpha +1}}}{\sqrt{\alpha +1}}\right)-\phi \left(x_0\right)\right]{\text{e}}^{-{\zeta }^2}\text{d}\zeta }\right|\\ \le \frac{2M}{\sqrt{\pi }}{\displaystyle {\int }_{-\infty }^{-N}{\text{e}}^{-{\zeta }^2}}\text{d}\zeta +\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_{-N}^N\left|\phi \left(x+\frac{2\zeta \sqrt{t^{\alpha +1}}}{\sqrt{\alpha +1}}\right)-\phi \left(x_0\right)\right|{\text{e}}^{-{\zeta }^2}\text{d}\zeta }+\frac{2M}{\sqrt{\pi }}{\displaystyle {\int }_N^{\infty }{\text{e}}^{-{\zeta }^2}}\text{d}\zeta \\ \le \frac{1}{\sqrt{\pi }}\frac{\epsilon }{3}{\displaystyle {\int }_{-N}^N{\text{e}}^{-{\zeta }^2}}\text{d}\zeta +4M\frac{\epsilon }{6M}\le \frac{\epsilon }{3}+\frac{2\epsilon }{3}=\epsilon \end{array}$

这样，我们就证明了式(13)是柯西问题(12)的有界解。

3.2. 解的唯一性和稳定性

极值原理是描述扩散、热传导等现象的热传导方程的重要特性。现在我们对齐次 $t^{\alpha }$ 型热传导方程

$u_t=t^{\alpha }{u}_{xx}$

的极值原理进行阐述。

极值原理 设 $u\left(x,t\right)$ 在矩形 $R_T\left\{0\le x\le l,0\le t\le T\right\}$ 上连续，并且在矩形内部满足齐次 $t^{\alpha }$ 型热传导方程，

则它在矩形的两个侧边( $x=0$ 及 $x=l,0\le t\le T$ )及底边( $t=0,0\le x\le l$ )上取到其最大值和最小值。换而言之，如果以 ${\Gamma }_T$ 表示 $R_T$ 的两侧边及底边所组成的边界曲线，那么成立着

$\underset{R_T}{\max }u\left(x,t\right)=\underset{{\Gamma }_T}{\max }u\left(x,t\right),\underset{R_T}{\max }u\left(x,t\right)=\underset{{\Gamma }_T}{\max }u\left(x,t\right).$

定理2 由极值原理，可以得到下述的 $t^{\alpha }$ 型热传导方程的初边值问题

$\{\begin{array}{l}{u}_t=t^{\alpha }{u}_{xx}+f\left(x,t\right),\\ u\left(x,0\right)=\phi \left(x\right),\\ u\left(0,t\right)=u_1\left(t\right),u\left(l,t\right)=u_2\left(t\right),\end{array}\left(0\le \alpha \le 1\right)$ (15)

在区域 $R_T$ 上的解是唯一的，而且连续地依赖于边界 ${\Gamma }_T$ 上所给定的初始条件及边界条件。

证明：假设问题(15)有两个解 $u_1$ 及 $u_2$ ，则其差为 $u=u_1-u_2$ 在区域 $R_T$ 内满足齐次方程

$u_t=t^{\alpha }{u}_{xx},$

而在 ${\Gamma }_T$ 上取零值，于是由极值原理就能得到在 $R_T$ 上 $u\equiv 0$ ，即这两个解在 $R_T$ 全同。

其次，如果初边值问题的两个解 $u_1$ 和 $u_2$ 在 ${\Gamma }_T$ 上满足

$\left|u_1-u_2\right|\le \epsilon ,$

则由极值原理就能得到在 $R_T$ 内也成立

$\left|u_1-u_2\right|\le \epsilon .$

这就证明了初边值问题(15)的解的稳定性。

4. $t^{\alpha }$ 型热传导方程解得数值模拟

从前面我们了解了初边值条件对于求解热传导方程的必要性，求解一个导热问题的热传导方程关键是要获得物体中的温度分布，有了温度分布导热量不难利用傅里叶定律得出所需要的热传导方程的解。对于某些工程类的问题，主要的目的是获得通过导热所传递的热量。

4.1. $t^{\alpha }$ 型热传导方程的差分法

4.1.1. $t^{\alpha }$ 型热传导方程的显示差分格式

考察一维 $t^{\alpha }$ 型热传导方程的初边值问题

$\{\begin{array}{l}{u}_t-t^{\alpha }{u}_{xx}=0\left(0<x<l,t>0,0\le \alpha \le 1\right),\\ {u|}_{t=0}=\phi \left(x\right)\left(0<x<l\right),\\ {u|}_{x=0}=u_1\left(t\right),{u|}_{x=l}=u_2\left(t\right)\left(t>0\right),\end{array}$ (16)

解的数值模拟。

为了建立这一类初边值问题(16)的差分格式，首先我们在它的求解区域内 $R\left(0\le x\le l,t\ge 0\right)$ 上作矩形网络，在x轴上以步长 $\Delta x=\frac{l}{J}$ (J是一个正整数)，把区间 $\left[0,l\right]$ 进行J等分，并于各分点作平行于t轴的网格线。在t轴一步长 $\Delta t$ 作平行于x轴的网格线。如上所作的两族直线的交点 $\left(x_j,t_n\right)=\left(j\Delta x,n\Delta t\right)$ 称为网格节点。

假设热传导方程的初边值问题的解在区域内部适当光滑时，对任意的区域内部的节点 $\left(x_j,t_n\right)$ 利用泰勒公式可以得到

$\frac{u_j^{n+1}-u_j^n}{\Delta t}-{\left(\frac{\partial u}{\partial t}\right)}_j^n=\frac{\Delta t}{2}\frac{{\partial }^{2}u\left(x_j,\tilde{t}\right)}{\partial t^2}\left(t_n\le \tilde{t}\le t_{n+1}\right),$

$\frac{u_j^{n+1}-2u_j^n+u_{j-1}^n}{{\left(\Delta x\right)}^2}-{\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_j^n=\frac{{\left(\Delta x\right)}^2}{12}\frac{{\partial }^{4}u\left(\tilde{x},t_n\right)}{\partial x^4}\left(\left|\tilde{x}-x_j\right|\le \Delta x\right).$

由上面的两式，我们可以得到

$\begin{array}{l}\frac{u_j^{n+1}-u_j^n}{\Delta t}-{\left(\Delta t\right)}^{\alpha }\frac{u_j^{n+1}-2u_j^n+u_{j-1}^n}{{\left(\Delta x\right)}^2}\\ =\frac{\Delta t}{2}\frac{{\partial }^{2}u\left(x_j,\tilde{t}\right)}{\partial t^2}-{\left(\Delta t\right)}^{\alpha }\frac{{\left(\Delta x\right)}^2}{12}\frac{{\partial }^{4}u\left(\tilde{x},t_n\right)}{\partial x^4}.\end{array}$

基于此网格节点，式中的初始条件和边界条件可以化为

$\{\begin{array}{l}{u}_j^0=\phi \left(j\Delta x\right)\left(j=1,\cdots ,J-1\right),\\ u_0^n=u_1\left(n\Delta t\right),u_J^n=u_2\left(n\Delta t\right)\left(n=0,1,2,\cdots \right).\end{array}$

就这样，我们将 $t^{\alpha }$ 型热传导方程的初边值问题数值解用差分法表示了出来，得到显式差分格式。

可以很明显的看到一个问题，上述所作的差分方程的解与原来初边值的热传导方程(16)一般是不相同

的。其中，我们将所忽略掉的项，通常叫做截断误差。这里的误差是 $O\left(\Delta t\right)+O\left({\left(\Delta x\right)}^2\right)$ 。

令

$\lambda =\frac{{\left(\Delta t\right)}^{\alpha +1}}{{\left(\Delta x\right)}^2},$

那么我们的显式差分格式可以写成

$\{\begin{array}{l}{U}_j^{n+1}=\lambda U_{j+1}^n+\left(1-2\lambda \right)U_j^n+\lambda U_{j-1}^n\\ U_j^0=\varphi \left(j\Delta x\right)\\ U_o^n=u_1\left(n\Delta t\right),U_J^n=u_2\left(n\Delta t\right)\end{array}\left(j=1,\cdots ,J-1;n=0,1,2,\cdots \right).$ (17)

它表示未知函数在第 $n+1$ 排任一节点上之值U依赖于它在第n排上的三个节点上的值。

4.1.2. 差分格式的稳定性

设计算初始条件时在各节点上产生的舍入误差是

${\epsilon }_j\left(j=1,2,\cdots ,J-1\right),$

此时，我们求得的差分方程的近似解 ${\tilde{U}}_j^n$ 满足

$\{\begin{array}{l}{\tilde{U}}_j^{n+1}=\lambda {\tilde{U}}_{j+1}^n+\left(1-2\lambda \right){\tilde{U}}_j^n+\lambda {\tilde{U}}_{j-1}^n\\ {\tilde{U}}_j^0=\varphi \left(j\Delta x\right)+{\epsilon }_j\\ {\tilde{U}}_0^n=u_1\left(n\Delta t\right),{\tilde{U}}_j^n=u_2\left(n\Delta t\right)\end{array}\left(j=1,2,\cdots ,J-1;n=0,1,2,\cdots \right).$ (18)

我们用 $V_j^{n+1}={\tilde{U}}_j^n-U_j^n$ 表示差分方程近似解与差分方程精确解的偏差，得到 $V_j^n$ 满足

$\{\begin{array}{l}{V}_j^{n+1}=\lambda V_{j+1}^n+\left(1-2\lambda \right)V_j^n+\lambda V_{j-1}^n\\ V_j^0={\epsilon }_j\\ V_0^n=V_J^n=0\end{array}\left(j=1,2,\cdots ,J-1;n=0,1,2,\cdots \right).$ (19)

这样，我们能够给出热传导方程的差分格式的稳定性定义。如果对于任意给定的 $\epsilon >0$ ，存在与 $\Delta t$ 和

$\Delta x$ 无关的 $\delta$ ，使得当 $\underset{j}{\max }\left|V_j^0\right|\le \delta$ 时，对于一切 $n=0,1,\cdots ,N$ ，成立

$\underset{j}{\text{max}}\left|V_j^n\right|\le \epsilon .$

我们称方程(17)是稳定的。

定理3当 $\lambda \le \frac{1}{2}$ 时， $t^{\alpha }$ 型热传导方程的初边值问题(16)显式差分格式是稳定的。

证明：当 $\lambda \le \frac{1}{2}$ 时，由(19)可知，对任意的j，有

$\left|V_j^{n+1}\right|\le \lambda \left|V_{j+1}^n\right|+\left(1-2\lambda \right)\left|V_j^n\right|+\lambda \left|V_{j-1}^n\right|\le \underset{j}{\max }\left|V_j^n\right|.$

所以

$\underset{j}{\text{max}}\left|V_{j+1}^n\right|\le \underset{j}{\max }\left|V_j^n\right|,$

对于任意的n都成立，因此对一切n成立

$\underset{j}{\text{max}}\left|V_j^n\right|\le \underset{j}{\max }\left|V_j^0\right|=\underset{j}{\max }\left|{\epsilon }_j\right|.$

由此证得差分格式的稳定性。

4.1.3. 隐式格式以其稳定性

从上面的讨论可以看到，用显式格式进行数值求解 $t^{\alpha }$ 型热传导方程的定界问题，优点是计算比较简

便，但是由于必须满足稳定性条件 $\lambda =\frac{{\left(\Delta t\right)}^{\alpha +1}}{{\left(\Delta x\right)}^2}\le \frac{1}{2}$ ，因此，t方向要满足 $\Delta t\le \sqrt[\alpha +1]{\Delta x^2/2}$ 。

我们为了提高数值解的精确度，必须缩小步长 $\Delta x$ ，此时 $\Delta x$ 的变化还要相应变小，考虑到 $\alpha$ 的取值，变小的程度会与 $\alpha$ 有关。下面我们为了避免显式格式所造成的计算时长大大加长的缺点，提出一种隐式差分格式。

若 $u\left(x,t\right)$ 是初边值问题的精确解，利用泰勒展开公式可以得到

$\begin{array}{l}\frac{u_j^{n+1}-u_j^n}{\Delta t}-{\left(\Delta t\right)}^{\alpha }\frac{u_{j+1}^{n+1}-2u_j^{n+!}+u_{j-1}^{n+1}}{{\left(\Delta x\right)}^2}\\ =-\frac{\Delta t}{2}\frac{{\partial }^{2}u\left(x_j,\tilde{t}\right)}{\partial t^2}-{\left(\Delta t\right)}^{\alpha }\frac{{\left(\Delta x\right)}^2}{12}\frac{{\partial }^{4}u\left(\tilde{x},t_{n+1}\right)}{\partial x^4}.\end{array}$

略去上式右端部分，就可以得到热传导方程的另一差分近似。这样，我们将初边值条件像前面一样处理，就得带了我们的隐式差分格式：

$\{\begin{array}{l}-\lambda U_{j-1}^{n+1}+\left(1+2\lambda \right)U_j^{n+1}-\lambda U_{j+1}^{n+1}=U_j^n\\ U_0^{n+1}=u_1\left(\left(n+1\right)\Delta t\right)\\ U_J^{n+1}=u_2\left(\left(n+1\right)\Delta t\right)\end{array}\left(j=1,\cdots ,J-1;n=0,1,2,\cdots \right).$

由上面的推导，容易知道我们的隐式差分格式的截断误差同样是 $O\left(\Delta t\right)+O\left({\left(\Delta x\right)}^2\right)$ 。这样的差分格

式无需对歩长加任何条件，所以是无条件稳定的。

4.2. 数值模拟

前面我们介绍了一维 $t^{\alpha }$ 型热传导方程的初边值问题的分析解。在实际工程计算中我们会遇到这样一类热传导方程，这些热传导问题具有复杂几何形状或者边界条件，因为数学上的困难没有得到很好的解决，无法获得解析解。另一方面，最近三十多年来，随着计算机这门技术的高速发展，对一般性的物理问题进行离散化求出数值解的方法发展非常迅速，这样的技术发展促进了数学物理理论的发展壮大，使得理论指导实践更加贴近。这些数值方法包括了有限差分法，有限元法及有限边界法等。

下面我们要举例应用Matlab建模构建图形，更加形象直观的理解热量分布与时间、空间的分布关系。对如下一维 $t^{\alpha }$ 型热传导方程的初边值问题解进行数值模拟

$\{\begin{array}{l}{u}_t=t^{\alpha }{u}_{xx}\\ u\left(x,0\right)=200\\ u\left(0,t\right)=u\left(l,t\right)=0\end{array}\left(0\le x\le \pi ,0\le \alpha \le 1,t>0\right).$

我们可以看到当 $\alpha =0$ 时，此方程与我们过去的一维热传导方程是一样的，是热传导系数与时间无关的

情形，当 $0<\alpha \le 1$ 此方程的数值解会随着 $\alpha$ 的变化而变化，下面我们比较地来看 $\alpha =0,\alpha =\frac{1}{4},\alpha =\frac{1}{2},\alpha =1$ 的情形。

Figure 1. Finite difference method solution (The coefficient is 0)

图1. 有限差分法解(系数为 $\alpha =0$ )

Figure 2. Finite difference method solution (The coefficient is 1/4)

图2. 有限差分法解(系数为 $\alpha =1/4$ )

Figure 3. Finite difference method solution (The coefficient is 1/2)

图3. 有限差分法解(系数为 $\alpha =1/2$ )

Figure 4. Finite difference method solution (The coefficient is 1)

图4. 有限差分法解(系数为 $\alpha =1$ )

图1数值模拟的稳定性系数为0.4523.所以，此一维热传导方程的差分格式的解收敛于原初边值问题并且是稳定的。图2数值模拟的稳定性系数为0.1315，图3数值模拟的稳定性系数为0.0382，图4数值模拟的稳定性系数为0.0032。我们发现，在我们x的步长与t的步长保持不变的情况下，我们的热传导系数 $\alpha$ 的变化会导致我们的有限差分格式的稳定性受到影响。我们通过Matlab实验可以得到如下结论：

当热传导系数 $\alpha =0$ 时，我们要研究的差分格式的位移步长和时间步长适当的话，利用计算机可以很好的模拟出一维 $t^{\alpha }$ 型热传导方程的数值解。当 $0<\alpha \le 1$ 时，x的步长与t的步长比例会影响到对一维 $t^{\alpha }$ 型热传导方程的数值解，随着 $\alpha$ 的增大，传热速率增加，热传导系数的增大导致传热加快，描述出了在实际生活中容易受到温度影响的材料的热传递方程解的数值模拟。

5. 小结

在本文中，我们对 $t^{\alpha }$ 型热传导方程这一类偏微分方程有了更加深刻的认识，在数学物理方程中的

$\frac{\partial u}{\partial t}=a^2\frac{{\partial }^{2}u}{\partial x^2}$ 热传导方程存在性、唯一性以及稳定性研究的基础上，从实际材料中热传导系数会随着时间和温度的变化而导致传热系数变化的问题出发，建立了 $\frac{\partial u}{\partial t}=t^{\alpha }\frac{{\partial }^{2}u}{\partial x^2}$ 的一维 $t^{\alpha }$ 型热传导方程，这是一个系

数与温度的幂次有关的偏微分方程。这样，我们利用所学的分离变量及常微分方程的知识，求出它的通解。在这个基础上，我们由此推出这一类偏微分方程解的存在性与唯一性，并进行论证。然后，本文利用了对 $t^{\alpha }$ 型热传导方程解进行数值模拟的基本方法——有限差分法，这样我们按照时间t增加的方向对每一个具有代表性的点进行逐点求解，在对新建立的与时间t的 $\alpha$ 次有关的热传导方程的解进行数值模

拟后，我们发现 $\alpha$ 在 $\left[0,1\right]$ 区间，随着 $\alpha$ 的增加， $t^{\alpha }$ 也会逐渐增大，热能的扩散会越来越快，这一结论也

与我们的生活实践相符合，当我们的材料受热的影响，温度升高而导致热传导系数增大时，热的扩散就会加快。而这一结论也同样可以应用于关于扩散的其他方面。

基金项目

本论文的完成受到湖南省自然科学基金(2019JJ50659)，长沙理工大学双一流国际合作拓展项目(2019IC39)的资助。