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

浅夏110

例 4 触煤反应装置内温度及转换率的分布
) ?6 ~) k" Q* o3 S$ z
以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下：
3 o0 Y2 Q. Q5 ]. d. F1 P7 g# z
" N$ V) ?. X/ e8 j
/ `; u  E% r; h" e% d
其中T 为温度(℃)， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为

5 F" s* y6 v/ R: J4 V7 q
, U6 D$ d# ~- u& e8 Z4 ~0 p
+ \& a, s. b, u此外，式中之相关数据及操作条件如下：

. V) c4 B* t5 W/ M6 e. ?! H(i)反应速率式



4 f: D1 c; [( |( o3 J( l其中 P 表示分压(atm)，而速率参数为
6 m0 n! ?7 Q( c' X9 q$ o* [


4 F" m6 b: @) _上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。

  p! W. r. t; {1 |(ii)操作条件及物性数据



6 w) `5 r* ~5 h. G* U, J

7 M$ I8 m9 s$ ~. Z题意解析：
- P9 s5 k4 Q5 \4 j5 L


# z$ y+ s0 c( m" O
将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。
3 A  E* [: J! o. X
MATLAB 程序设计 将原方程改写成如式（35）的标准式
6 B5 B6 o# }3 ]. B) B, y2 Q

. t! X9 [1 T* R- _: r
     因此
4 ~3 \7 q$ B; B$ I- N" O
4 \2 \% W0 j% l

* P- ~! k2 v4 g- p' r0 B根据以上的分析，可编写 MATLAB 程序求解此 PDE 问题，其参考程序如下：
$ N* P7 ^; e+ i: Y- g
+ c; O6 y* j/ S8 d8 p3 _+ a4 v& q- s' gfunction ex60_3_1
' t  i, O0 E7 u) G%******************************
% 触媒反应器内温度及转化率的分布
5 g. W& z: H. s4 w- E: E9 w3 t1 J, o%******************************
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
%******************************
" a( B# B+ N0 h+ v% 给定数据
* p4 b8 \: k8 {) h' m  p6 z; @9 D- B%******************************
0 o9 U7 ~! ?6 r- U# y: s- x3 tPt=1.25; %总压(atm)
( J. t4 q9 j0 Q! x/ A9 ]/ Q1 hrw=0.025; %管径(m)
Tw=100+273; %壁温(℃)
& a( V0 h  j) m( BG=631; %质量流率(kg/m2hr)
M=30;
y0=0.0323;
Mav=4.47;
* p/ a# `" F+ r6 s9 ~$ O3 @# Erho_B=1200;
Cp=1.74;
dHr=-49250;
4 T* a- p$ o" ~7 }- g+ jh0=65.8;
8 @2 N+ _1 B. K5 bT0=125+273;
Lw=1;
: ]& B( L7 H0 v/ b( m. G- Eu=8.03;
- I. Y% `7 I3 K; DR=1.987;
ke=0.65;
. \5 K# P5 R; k6 v( ihw=112;
De=0.755;
%********************
m=1;
%********************
% 取点
%********************
$ l. ~% M# g. f! c. Gr=linspace(0,rw,10);
7 v$ z* a% R0 `/ H, E4 d  zL=linspace(0,Lw,10);
%***********************
% 利用 pdepe 求解
%***********************
4 b3 C9 h$ y, W8 W2 ?sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L);
T=sol(:,:,1); %温度
f=sol(:,:,2); %反应率
9 {1 S' F$ o1 Z%***********************
% 绘图输出
%***********************
! s+ H0 Q, N2 ]- K& m# Qfigure(1)
9 X/ K2 v. s) T$ Q6 Jsurf(L,r,T'-273)
" y! s6 }# F1 P0 @( x8 etitle('temp')
8 S& O) _3 B- k! D% ~xlabel('L')
ylabel('r')
5 o% j( L8 a+ W1 V6 f# }4 Wzlabel('temp (0C)')
%
0 b8 ?) _/ e4 Y# c6 ?( q% [figure(2)
' }6 u' }# k0 Ksurf(L,r,f')
9 H6 L7 G! p- o* e) Xtitle('reaction rate')
xlabel('L')
%初始条件函数
/ u/ y5 {3 [, V8 ?& B%**********************************
function u0=ex20_3_1ic(x)
6 z5 k" v7 h2 H1 o) d* o: I- au0=[125+273 0]';
%**********************************
% 边界条件档
%**********************************
function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L)
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
pl=[0 0]';
) j& k4 e% X; n. u( {. Aql=[1 1]';
pr=[hw*(ur(1)-Tw) 0]';
qr=[G*Cp 1]';
# N. X9 G& P( E/ ]# Q/ K) Lylabel('r')
4 v  {' X4 ^& Y6 K' e: _" [( Pzlabel('reaction rate')
, c, D! ~/ Z5 ?- I5 r; v# w0 m%*************************************************
% PDE 函数
%*************************************************
function [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
1 Q& p/ v+ r. s' T1 kT=u1(1);
5 x) x: u2 r# V+ s( If=u1(2);
8 ^! I; \. L3 R4 l  O% i: N%
5 @9 \0 q% y% g" lk=exp(-12100/(R*T)+32.3/R);
- Q. i, V) X( D/ S, ^4 BKh=exp(15500/(R*T)-31.9/R);
Kb=exp(11200/(R*T)-23.1/R);
Kc=exp(8900/(R*T)-19.4/R);
- h& ~# A" g- a' h5 F%
a=1+M-3*f;
! w( Z- `. c- Q- Yph=Pt*(M-3*f)/a;
( k3 P( R1 O( w- p& dpb=Pt*(1-f)/a;
pc=Pt*f/a;
%
rA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4;
6 ]0 E& s7 Y; Y0 v% {1 d: A' E%
c1=[1 1]';
: v( C0 X: K8 W1 x; d0 ^8 _f1=[ke/(G*Cp) De/u]'.*DuDr;
%s1=[ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw)
  ^; B' P+ w1 W# l) gs1=[-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0)];
%**********************************
4 b- t9 [. m9 |# O* i+ S
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
+ _, O+ ~' v' E/ d原文链接：https://blog.csdn.net/qq_29831163/article/details/89711536
4 f$ R% L& y0 K4 x" ~6 t9 W0 ~
" }, ?4 W; n5 w) A
; B8 N& b; t1 I7 }