formoving的个人博客分享 http://blog.sciencenet.cn/u/formoving

博文

使用MATLAB中pdepe函数求解一维偏微分方程

已有 32252 次阅读 2014-5-4 13:46 |个人分类:社会感叹|系统分类:科研笔记

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

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

   通过pdepe函数的m文件代码或者查阅help文档中所附的参考文献,可以得知pdepe函数是基于线上法(method of linesMOL)求解一维的抛物型和椭圆型方程的。通过Find Files可以在MATLAB中找到pdepe函数的m文件代码,除去说明型文字,整个pdepe函数仅有250行左右的代码。线上法的核心在于使用近似方法对偏微分方程的空间微分项进行离散,空间微分项离散后,偏微分方程将转化为不再显含空间独立变量的常微分方程组。MATLABpdepe函数主要基于的文献[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函数来指定cfs的表达式,通过icfun函数来指定初始条件u0的表达式,通过bcfun函数来指定pq在左右边界上plprqlqr四个量的表达式。m的值可以为012,是为了方便求解圆柱和球坐标下的一维问题,一般直线坐标问题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}}})$

主要参数可以简单的设定为:

ρg

1.293  kg/m3

ρs

2200  kg/m3

cg

1075  J/kg·K

cs

1200  J/kg·K

kg

0.035  W/(m·K)

ks

3.4  W/(m·K)

h

65  W/(m2·K)

av

1.69×103  m-1

ε

0.56

um

1.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}}})$

分别用u1u2代替TgTspdepe函数的求解偏微分方程形式中的cfs分别可以表示为

${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})]$

cfs向量在控制方程和边界条件的乘积运算中均为点乘。通量形式的左边界条件为:

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

左边界条件为:dTs/dx(x=0m)=0;Tg(x=0m)=300Ku1u2代替TgTs。左边界条件可以写成:

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

上式可以用于指plql右边界条件为:dTs/dx(x=1m)=0;dTg/dx(x=1m)=0。u1u2代替TgTs。右边界条件可以写成:

$[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给出输入的pqpdepe函数的详细介绍除了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).







https://blog.sciencenet.cn/blog-1305714-791201.html


下一篇:真的有600所本科院校转向职业教育?
收藏 IP: 159.226.48.*| 热度|

3 王振亭 田灿荣 曲杨

该博文允许注册用户评论 请点击登录 评论 (7 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-28 21:44

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部