计算机应用基础-偏微分方程求解

cAm,n

Derz 2u r 2
1

2m
c c Am1,n1 Am1,n1

c Am1, n 1

1 2m
c Am1,n
cAm1,n1 2cAm,n1 cAm1,n cAm1,n 2cAm,n




c
Am1,n

B rA m,n1/ 2 z

c11u1 c21u1
n n
c12u2 q11u1 q12u2 g1 c22u2 q21u1 q22u2 g2
其中 n为边界外法向单位向量，g, q, h, r为在边界 上定义的函数
（3）混合边界条件
n n

c11u1 c21u1
Tm1,n1 2Tm,n1 Tm1,n1 Tm1,n 2Tm,n Tm1,n
将上三式代入二维模型的物料衡算方程
6.4 偏微分方程在化工中的应用
Tm,n1
Tm,n

erz 2ucp r
2
1

2m
T T m1,n1 m1,n1
6.4 偏微分方程在化工中的应用
固定床二维反应器模型
T z

a1
2T r 2
1 r
T r

b1rc

0
x z

a 2

2x r 2

1 r
x r


b2rc

0
Z为反应管轴向距离，r为径向距离，方程系数如下：
a1

e Gcp
a2

De u
– 在 MATLAB 下键入 pdetool
• 该界面分为四个部分
– 菜单系统 – 工具栏 – 集合编辑 – 求解区域
6.3 偏微分方程求解工具箱
菜单栏
工具栏
6.3 偏微分方程求解工具箱
5.3 偏微分方程求解工具箱
工具箱求解步骤： 1. 用options设置应用模式（可选） 2. 用Draw建立几何模型 3. 用Boundary菜单设定边界条件 4. 用PDE定义偏微分方程的类型和系数 5. 用Mesh菜单进行三角形网格划分及细化 6. 用Slove进行偏微分方程求解 7. 用Plot以图形方式显示结果
6.4 偏微分方程在化工中的应用
进入反应器，相当于总质量速率为G=2500kg.h-1.m2。反应管
外用速率为F 130kg h-1烟道气与反应混合物
逆流加热反应管，烟道气出口温度为620C。其
它数据：催化剂的堆积密度＝1440kg / m3，操作
压力P 1.2bar，乙苯的反应热H＝140000kJ / mol，
y,u, u x
, u x
当A，B，C为常数时，称为拟线性偏微分方程，可 分为三类：
B2 4AC 0 椭圆型方程
B2 4AC 0 抛物型方程
B2 4AC 0 双曲型方程
6.1 偏微分方程组求解
二 偏微分方程边界条件： （1）Dirichlet 边界条件 hu=r 也称为第一类边界条件，对于偏微分方程组，
u=assempde(b,p,e,t,c,a,f,u0) U0:初始条件，用于非线性方程求解
6.2 二阶偏微分方程的求解
例6-2，利用adaptmesh函数求解拉普拉斯方程，其 在弧上满足Dirichlet条件：
u=sin(2/3*atan2(y,x))
ME_6_3
6.2 二阶偏微分方程的求解
二 抛物线型偏微分方程
Dirichlet边界条件为 h11u1 h12u2 r1 h12u1 h22u2 r2
（2）Neumann边界条件
n cu qu g
6.1 偏微分方程组求解
也称为第三类边界条件；当q=0时，则变为第二类边界 条件。对于偏微分方程组，Neumann边界条件为：
n n
Tm1,n1
1 2m
Tm1,n
Tm1,n1 2Tm,n1 Tm1,n Tm1,n 2Tm,n




Tm1,n

B H rA m,n1/2 z
uc p
类似的可求得浓度的差分方程
cAm,n1
H rA m,n
uc p
2Tm,n

Tm1,n





cAm,n1/2

cAm,n

z 2
Der u r
2
1 2m

cAm1,n cAm1,n
B rA m,n
u

c Am1, n

2cAm,n
采用差分代替 微分，离散化 进行迭代。
如图，任一节 点的径向与轴 向位置可表示 为：
r mr
z nz
M-1, n+1 M, n+1 M+1, n+1
M, n+1/2
M-1, n
M, n M+1, n
r
z
6.4 偏微分方程在化工中的应用
对于点（M,n+1/2)的各项导数都可以用该处周围六个节点的 函数值计算的差商表示，如温度可表示为：
6.2 二阶偏微分方程的求解
parabolic函数用于求解抛物型偏微分方程的解，调用格 式如下：
u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d) b: 边界条件 u0: 初始条件 tlist；时间列表 u1:对应于tlist的解向量 p,e,t :网格数据
6.2 二阶偏微分方程的求解
一 椭圆型偏微分方程
6.2 二阶偏微分方程的求解
6.2 二阶偏微分方程的求解
adaptmesh 和assempde函数用于求解椭圆型偏微分方 程的解，调用格式如下：
[u, p, e, t]=adapmesh(g,b,c,a,f) g: 求解几何区域； b: 边界条件 u:解向量 p,e,t :网格数据
n n
c12u2 q11u1 q12u2 g1 h11 c22u2 q21u1 q22u2 g2 h21
6.1 偏微分方程组求解
三 偏微分方程数值解法
1. 有限差分法 2. 正交配置法 3. MOL法 4. 有限元法
6.1 偏微分方程组求解
四 采用pdepe( )函数求解一维偏微分方程
b1

H Gcp
b2

uc0
6.4 偏微分方程在化工中的应用
通过反应计量关系获得各组分的分压：
pE

1.21 x
11 x
pS

pH

1.2x 11 x
rc

15100 exp 1100 /
T111xx

1.2x 2
K11 x2

温度初始边界条件：
z 0, r 0.
分压（bar）; rc的单位为kmol/(h kg催化剂）
反应速率常数：k 12600exp11000/ T
平衡常熟： K 0.027exp[0.021T 773]
T表示温度，K
已知反应管的内径2a 10cm，乙苯和水蒸汽在6000C
下分别以0.069kmol h1和0.69kmol h1的加料速率
T
z m,n1/ 2
1 z
Tm,n1 Tm,n
T
r m,n1/ 2
1 4r
Tm1,n1 Tm1,n1 Tm1,n Tm1,n

2T r 2
m,n1/ 2

1 2(r)2
1. 问题定义及参数初始化 2. 网格化 3. 求解 4. 显示结果
6.3 偏微分方程求解工具箱
例6-5 求解二维动态热传导方程
T t

2T x2

2T y 2
,
0 x 15,0 y 20,t 0
初始条件： u t0 0, 0 x 15,0 y 20
u
6.4 偏微分方程在化工中的应用
上两式中(rA)m,n+1/2是指温度为Tm,n+1/2与浓度cAm,n+1/2 条件下的反应速率，有：
Tm,n1/2
Tm,n

z 2
er

uc
p
r



1 2 2m
B
Tm1,n Tm1,n Tm1,n
第六章 偏微分方程求解
6.1 偏微分方程组求解 6.2 二阶偏微分方程的数学描述 6.3 偏微分方程的求解界面应用举例 6.4 偏微分方程在化工中的应用
6.1 偏微分方程组求解
一 偏微分方程的分类
2u A x2
B
2u xy
C
2u y 2
D u x
E
u y

Fu

f
x,

c Am1, n




6.4 偏微分方程在化工中的应用
管中心和管壁处的边值条件也可用差分方程表示。管中心：
T0,n1
T0,n

2erz ucp r
2
Байду номын сангаасT1,n1 T0,n1 T0,n T1`,n
B H rA 0,n1/2 z
T T0 T 0 r
6.4 偏微分方程在化工中的应用
径向边界条件
ra
T z

C1
T r
C1

2ae Fc'p
质量初始边界条件
z0
x0
r 0 x 0; r a x 0
r
r
6.4 偏微分方程在化工中的应用
微分方程组的差分格式离散化
一 隐式差分（Crank-Nicholson)
6.3 偏微分方程求解工具箱
【例6-3】 求解椭圆型方程
2u

x
2

2u y2

2
u 100
5 x 5, 5 y 5 x 5 y 5
采用工具箱求解
6.3 偏微分方程求解工具箱
6.3 偏微分方程求解工具箱
利用PDE工具箱命令行求解偏微分方程:
例6-3：求解热传导方程：
u

t


2u x2

2u y 2

2u z 2


0 ,
u 0
G x, y, z0 x, y, z 1
ME_6_2
6.2 二阶偏微分方程的求解
三 双曲型偏微分方程
6.3 偏微分方程求解工具箱
• 启动偏微分方程求解界面
6.1 偏微分方程组求解
边界条件的函数描述：
【例6-1】
6.1 偏微分方程组求解
6.1 偏微分方程组求解
6.1 偏微分方程组求解
偏微分方程求解程序 “c7mpde”
function [c, f, s]=c7mpde(x, t, u, du) c=[1 ; 1] ; y=u(1)-u(2) ; F=exp(5.73*y)-exp(-11.46*y) ; s=F*[-1; 1] f=[0.024*du(1); 0.17*du(2)];
边界条件： u x0 u x15 u y0 u y20 100, (t 0)
ME_6_6.m
6.4 偏微分方程在化工中的应用
在一管式催化反应器中进行乙苯的催化脱氢反应，所用原料为 乙苯和水蒸汽的气体混合物。动力学方程为：
rc

k
pE

pS pH K

式中，pE , pS , pH 分别是乙苯、苯乙烯和氢气的
6.1 偏微分方程组求解
边界条件程序”c7mbc.m” function [pa, qa, pb, qb]=c7mpbc(xa, ua, xb, ub, t) pa=[0; ua(2)]; qa=[1; 0]; pb=[ub(1)-1; 0]; qb=[0; 1];
function u0=c7mpic(x) u0=[1; 0];
>>x=0: 0.05: 1; t=0: 0.05: 2; m=0; sol=pdepe(m, @c7mpde, @c7mpic, @c7mpbc, x, t); surf(x, t, sol( :, :, 1))
6.1 偏微分方程组求解
求解函数：ME_5_1.m, 结果如下
6.2 二阶偏微分方程的求解
床层有效导热系数e 0.45w.m1.k 1,有效扩散系数
与气流流线速之比De / u 0.000427m，反应混合物比热
容C p

2.18k J
/(k
g

k
)，烟
道气比热容C
' p
1.0kJ
/(k g.k )
试计算反应管轴向和径向温度分布和转化率分布，
并计算需要多长的反应器管道才能使乙苯的转化率达到45%。