偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

浅夏110

例 4 触煤反应装置内温度及转换率的分布

以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下：
% }  Q% |7 d# ]( C; g5 |9 x
) M% l5 ]& e6 ^8 u) p
& g) j- j% M5 M4 p6 v9 Z) @3 j; F
* [" ^: [; H( C 其中T 为温度(℃)， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为



此外，式中之相关数据及操作条件如下：

$ m' S- R6 `! E(i)反应速率式
' g# ]% E" z1 S


6 L) Q* s; U  w  M, s+ ~其中 P 表示分压(atm)，而速率参数为

, |$ ?  W, Z7 k1 s

4 W* \5 Z+ E) r. b# U$ p上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。

+ w0 q3 c: Q3 z" j  A+ c(ii)操作条件及物性数据
+ Q; @$ S! N* a# X/ e8 N
. a, a0 V& J! I. g
) ]9 S3 R1 d4 G8 i. E8 I
4 l+ Q  P' e0 q4 e2 X
* _% G1 z! J3 |5 a4 o; x
; \: m3 s* d4 \1 B( W  h# b: _* p5 g题意解析：




将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。

& |8 ~" G7 w% B8 h$ UMATLAB 程序设计 将原方程改写成如式（35）的标准式
' v% U( X5 d. _


     因此
1 @3 N6 S, m3 V( h: o


  J) _2 R0 R. D根据以上的分析，可编写 MATLAB 程序求解此 PDE 问题，其参考程序如下：
0 b. Y; @% |) z4 {4 r
function ex60_3_1
9 J3 r. t- U1 g; P3 i5 T%******************************
5 A5 L$ o  P$ i7 M+ e% 触媒反应器内温度及转化率的分布
7 C  Q* R* M: [5 ^! K  v%******************************
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
( U% X3 t3 n& @$ F) |%******************************
' f3 C+ x0 |. K8 R# H" t1 `$ |0 h% 给定数据
%******************************
% T2 l9 w6 m/ t( l: M& oPt=1.25; %总压(atm)
rw=0.025; %管径(m)
$ R1 u+ p! j. L% K1 qTw=100+273; %壁温(℃)
G=631; %质量流率(kg/m2hr)
9 j- m4 G1 {' B; c5 d$ XM=30;
+ X& }2 f/ }  Z0 _y0=0.0323;
Mav=4.47;
& {$ [1 W& t+ o2 `# Crho_B=1200;
! \, ]2 r) Z& L! P8 zCp=1.74;
dHr=-49250;
+ X( m, ?3 h4 o$ Gh0=65.8;
/ c2 J! }. A5 _/ e+ ]; |0 nT0=125+273;
Lw=1;
' g* Q) H$ I6 [& Eu=8.03;
R=1.987;
5 T) C) ]: A6 @ke=0.65;
9 w* ]; F4 O2 t4 a% Hhw=112;
De=0.755;
%********************
m=1;
%********************
1 ~* [8 c7 v  ]% 取点
+ H+ V. H/ S: O* T$ u( B%********************
r=linspace(0,rw,10);
L=linspace(0,Lw,10);
%***********************
+ h1 u4 ^+ i6 X% 利用 pdepe 求解
%***********************
sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L);
T=sol(:,:,1); %温度
f=sol(:,:,2); %反应率
%***********************
( k! D6 b$ y8 p# }  Z% Z% 绘图输出
$ v" c+ X- o5 T7 Z%***********************
figure(1)
surf(L,r,T'-273)
title('temp')
# E8 f6 j2 N8 f1 x. C4 ?xlabel('L')
ylabel('r')
zlabel('temp (0C)')
%
2 C4 p2 _6 i1 ^figure(2)
" |, `3 t& S) f* v2 }8 i0 osurf(L,r,f')
* @* w0 m2 _4 g2 i9 @' _  Dtitle('reaction rate')
xlabel('L')
%初始条件函数
%**********************************
function u0=ex20_3_1ic(x)
u0=[125+273 0]';
%**********************************
% 边界条件档
%**********************************
( {5 g" B$ Z2 p7 V" K7 ?% \function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L)
5 W) a9 ?0 \, Jglobal Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
pl=[0 0]';
ql=[1 1]';
9 N7 l! h% |# \4 g6 J) Cpr=[hw*(ur(1)-Tw) 0]';
qr=[G*Cp 1]';
9 k4 K9 e* s3 {6 O+ F+ s, @% [9 Qylabel('r')
" `: W5 i, ^+ t9 p+ h2 r3 d9 Ezlabel('reaction rate')
; |5 A/ m" }* T! O( E%*************************************************
' P: d% r, A8 j/ d* i% PDE 函数
%*************************************************
function [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)
; s8 F4 Y. F3 A$ sglobal Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
  K- p. g+ R! S8 n. l, OT=u1(1);
$ z) _1 T8 K3 X" e% Wf=u1(2);
7 c7 S# G7 b# i& O+ e+ M. U%
  r1 A' i$ \3 I) y1 B5 D5 zk=exp(-12100/(R*T)+32.3/R);
Kh=exp(15500/(R*T)-31.9/R);
- Y1 [" T$ }8 X( S: G* ]Kb=exp(11200/(R*T)-23.1/R);
Kc=exp(8900/(R*T)-19.4/R);
%
a=1+M-3*f;
ph=Pt*(M-3*f)/a;
pb=Pt*(1-f)/a;
- L. w2 M& P+ |8 A9 l1 ipc=Pt*f/a;
%
5 M- l7 J+ Z; ?/ YrA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4;
- w9 ?, u+ t- a7 @) n. _" E* i2 F%
$ T: q3 T9 e4 G$ _6 Bc1=[1 1]';
& C. i$ }: P3 e9 L6 Uf1=[ke/(G*Cp) De/u]'.*DuDr;
%s1=[ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw)
s1=[-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0)];
%**********************************
! V9 p* L- O( f' E3 W' S. f  X
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
3 Q, P# N, W* l/ L原文链接：https://blog.csdn.net/qq_29831163/article/details/89711536
) [- Q1 |( j" W" Z/ q$ G; O

2 ]5 t! N" D" [& G% S