由于自己科研水平较低，记录的各种体会更多的是给自己做个小结，错误之处，欢迎大家指正。

使用MATLAB求解偏微分方程或者方程组，大致有三类方法。第一种是使用MATLAB中的PDE Toolbox，PDE Toolbox既可以使用图形界面，也可以使用命令行进行求解。最近版本的MATLAB中将工具箱称为APPS，有点像手机系统上的各种应用，Toolbox也确实有点像是基于MATLAB平台的各种应用。PDE Toolbox只能求解二维问题，另外不能求解刚性问题，PDE Toolbox是基于有限元方法的，其中装配矩阵函数assempde只适用于非刚性矩阵，所以类似反应过程的这类刚性问题对流扩散方程是无法使用PDE Toolbox求解的。第二种就是使用MATLAB中的m语言进行编程计算，相比与Fortran和C等语言，MATLAB中编程计算有许多库函数可以使用，对于大型矩阵的运算也要方便的多。当然使用m语言编程计算也有劣势，如循环计算的效率很低，所以在MATLAB中尽量使用向量运算代替循环运算。第三种就是使用pdepe函数，MATLAB中pdepe函数的作用是求解一维抛物型和椭圆形方程。

通过pdepe函数的m文件代码或者查阅help文档中所附的参考文献，可以得知pdepe函数是基于线上法(method of lines，MOL)求解一维的抛物型和椭圆型方程的。通过Find Files可以在MATLAB中找到pdepe函数的m文件代码，除去说明型文字，整个pdepe函数仅有250行左右的代码。线上法的核心在于使用近似方法对偏微分方程的空间微分项进行离散，空间微分项离散后，偏微分方程将转化为不再显含空间独立变量的常微分方程组。MATLAB中pdepe函数主要基于的文献[1]即针对线上法的空间项离散提出一种新的方法。实际上MOL在求解偏微分方程的诸多方法中，只能算是非主流。

pdepe函数求解的偏微分方程形式为：

$c(x,t,u,{{\partial u} \over {\partial x}}){{\partial u} \over {\partial t}} = {x^{ - m}}{\partial \over {\partial x}}({x^m}f(x,t,u,{{\partial u} \over {\partial x}})) + s(x,t,{{\partial u} \over {\partial x}})$

可以将其中的f理解为通量，s为源项。方程初始条件的形式为：

$u(x,{t_0}) = {u_0}(x)$

控制方程边界条件的形式为：

$p(x,t,u) + q(x,t)f(x,t,u,{{\partial u} \over {\partial x}}) = 0$$" style="font-family:宋体;line-height:29px;text-align:center;$

pdepe函数的调用形式为：sol =pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)。通过pdefun函数来指定c，f，s的表达式，通过icfun函数来指定初始条件u0的表达式，通过bcfun函数来指定p，q在左右边界上pl，pr，ql，qr四个量的表达式。m的值可以为0、1、2，是为了方便求解圆柱和球坐标下的一维问题，一般直线坐标问题m取0。求解偏微分方程组时，只需要将各个量表示为向量的形式即可。

pdepe函数无法求解较为复杂的问题，所以选用一个比较简单的气固对流换热过程作为示例。这个简单问题是环境气流流进入给定温度的颗粒填充床中，填充颗粒物性参数给定，对流换热系数也给定，求解气固相的温度变化。气固相的控制方程如下：

$\varepsilon {\partial \over {\partial t}}({\rho _{\rm{g}}}{c_{\rm{g}}}{T_{\rm{g}}}) = - {\partial \over {\partial x}}\left( {{\rho _{\rm{g}}}{c_{\rm{g}}}{T_{\rm{g}}}{u_{\rm{m}}}} \right) + {{{\partial ^2}} \over {\partial {x^2}}}(\varepsilon {k_{\rm{g}}}{T_{\rm{g}}}) + h{a_{\rm{v}}}({T_{\rm{s}}} - {T_{\rm{g}}})$

$\left( {1 - \varepsilon } \right){\partial \over {\partial t}}({\rho _{\rm{s}}}{c_{\rm{s}}}{T_{\rm{s}}}) = (1 - \varepsilon ){k_{\rm{s}}}{{{\partial ^2}{T_{\rm{s}}}} \over {\partial {x^2}}} - h{a_{\rm{v}}}({T_{\rm{s}}} - {T_{\rm{g}}})$

主要参数可以简单的设定为：ρg1.293  kg/m3ρs2200  kg/m3

cg1075  J/kg·Kcs1200  J/kg·K

kg0.035  W/(m·K)ks3.4  W/(m·K)

h65  W/(m2·K)av1.69×103m-1

ε0.56um1.5  m/s

代入参数后方程简化为：

$778.4{{\partial {T_{\rm{g}}}} \over {\partial t}} = - 1389.98u{{\partial {T_{\rm{g}}}} \over {\partial x}} + 0.02{{{\partial ^2}{T_{\rm{g}}}} \over {\partial {x^2}}} + 110 \times {10^3}({T_{\rm{s}}} - {T_{\rm{g}}})$

$1.16 \times {10^6}{{\partial {T_{\rm{s}}}} \over {\partial t}} = 1.5{{{\partial ^2}{T_{\rm{s}}}} \over {\partial {x^2}}} - 110 \times {10^3}({T_{\rm{s}}} - {T_{\rm{g}}})$

分别用u1，u2代替Tg，Ts。pdepe函数的求解偏微分方程形式中的c，f，s分别可以表示为：

${f_1}(x,t,{u_1},{{\partial {u_1}} \over {\partial x}}) = - 2084.97{u_1} + 0.02{{\partial {u_1}} \over {\partial x}}$

${f_2}(x,t,{u_2},{{\partial {u_2}} \over {\partial x}}) = 1.5{{\partial {u_2}} \over {\partial x}}$

$f = [{f_1};{f_2}] = [ - 2084.97{u_1} + 0.02{{\partial {u_1}} \over {\partial x}};1.5{{\partial {u_2}} \over {\partial x}}]$

$c = [{c_1};{c_2}] = [778.4;1.16{\rm{e}}6]$

$s = [{s_1};{s_2}] = [110{\rm{e}}3({u_2} - {u_1}); - 110{\rm{e}}3({u_2} - {u_1})]$

c，f，s向量在控制方程和边界条件的乘积运算中均为点乘。通量形式的左边界条件为：

$\vec p(x,t,u) + \vec q(x,t).*[{f_1};{f_2}] = 0$

左边界条件为：dTs/dx(x=0m)=0；Tg(x=0m)=300K。用u1，u2代替Tg，Ts。左边界条件可以写成：

$[6.25{\rm{e}}5;0] + [1;1].*[ - 2084.97{u_1};1.5{{\partial {u_2}} \over {\partial x}}] = 0$

上式可以用于指定pl，ql。右边界条件为：dTs/dx(x=1m)=0；dTg/dx(x=1m)=0。用u1，u2代替Tg，Ts。右边界条件可以写成：

$[50*2084.97{u_1};0] + [50;1].*[ - 2084.97{u_1} + 0.02{{\partial {u_1}} \over {\partial x}};1.5{{\partial {u_2}} \over {\partial x}}] = 0$

初始条件为：Ts(0,x)= 1000K；Tg(0,x)= 300K。即u0=[1000,1000]。

在x方向0~1m上取100个计算节点，时间从0~50s取100额计算节点。u1和u2(也就是气流和固体颗粒的温度TgTs)的计算结果如下：

这两个温度变化的计算结果是基本合理的。应该说使用pdepe求解一维的简单对流扩散方程还是比较便捷的，需要注意的是控制方程的边界条件设置。在上面求解的这个小问题中，气体的导热项是本是可以忽略的，但在使用pdepe求解控制方程的设定中则不能省掉，否则边界条件设置过程中，将无法基于f给出输入的p，q。pdepe函数的详细介绍除了MATLAB中的help文档，还可以参考[2][3]这两个介绍文档。

[1]   Skeel R D, BerzinsM. A method for the spatial discretization of parabolic equations in one spacevariable[J]. SIAM journal on scientific and statistical computing, 1990, 11(1):1-32.

[2]  Howard, P."Partial Differential Equations in MATLAB 7.0."University of Maryland(2005).

[3]  Chung, Jaywan. "Differential Equations inMATLAB 7." (2012).

转载本文请联系原作者获取授权，同时请注明本文来自邓浩鑫科学网博客。

收藏

分享

分享到: