数学建模社区-数学中国
标题:
偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:27
标题:
偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布
例 4 触煤反应装置内温度及转换率的分布
$ d! v* a% q* u$ k' }2 t$ ]+ E. N
8 {/ q W d! [% h4 ~
以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下：
; \* m5 K/ B3 x* K: C
+ X" \/ S% t1 I
- k. H1 m% a) @
( ^' U. P( B, D2 K
其中T 为温度(℃)， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为
$ t1 m; }' k1 h$ f" v! y& ^
4 z$ r# h X, U
- G+ b6 f: P5 e
) j4 K! z0 {1 c! q5 P1 S; W9 L7 p3 V
此外，式中之相关数据及操作条件如下：
8 e& l* f; L/ f+ v: l
0 r5 n; o6 ]; o+ d: u2 `4 d
(i)反应速率式
0 V$ P# S0 s& K f5 C
; b- Z: }0 R, n k9 @3 Y/ u
. K# z1 l! X6 _& ]( v t8 s
+ t" K6 ^1 ^- }3 M! F* l
其中 P 表示分压(atm)，而速率参数为
& m- k7 E. v& G
7 t& Y( t; e& F
! n$ \: k2 h- ^$ _
) \$ ]0 x+ x/ D
上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。
( k/ w! A- T ~/ s4 R
; X' r8 a4 e. S9 i0 u& A r; ^
(ii)操作条件及物性数据
% n2 T9 l; u: ]. ^
3 [- E5 h. V$ v0 o
" Y6 [& O1 g3 g$ F+ c+ V+ L- V1 p
! f- n" e) E6 R, E/ M: Z: o
$ p' x: t3 @; T3 x- m* m& z
! |3 E# F) X$ B: b
题意解析：
3 A1 c% i2 n$ S3 c
; ]3 G3 m6 P$ k% L$ p; X/ c
# t8 Z) Y6 F6 ~
% ~5 Y9 p& B) n8 E G
: n' ]- n" s0 x) j& t
将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。
, D) I: d9 M/ s5 q- g. d
7 f" ^ s: `6 s4 t! A9 ~- b1 |
MATLAB 程序设计 将原方程改写成如式（35）的标准式
8 Y6 {; }, t0 p n% r, y
6 w( P/ u! ]: p) B6 i. e, I
1 I6 b! A6 @& g' z
. L( G v1 F( ]) H# K
因此
+ W* N( N+ T' s
( m! ]" K0 q! z0 R
M- T. ?3 L% C1 x: _$ L
9 o6 o+ A: h* i/ K( N5 U1 [$ m" T
根据以上的分析，可编写 MATLAB 程序求解此 PDE 问题，其参考程序如下：
]+ B& m* R& p( _
! h& a9 W$ V3 J: K7 h& r, I
function ex60_3_1
# o6 s; u7 V- ^6 i
%******************************
* u6 g/ u' P6 [* }2 b
% 触媒反应器内温度及转化率的分布
$ e/ ?8 t5 x1 c
%******************************
) F3 ~: k* d2 B9 Z6 U, U
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
2 B+ {! d: l0 e. L5 H8 j
%******************************
1 B! R7 f9 v5 A+ l; A
% 给定数据
4 r5 {( Q) u2 n+ Y, {# a5 I6 f
%******************************
4 [ V, C& n2 ^7 Z4 \) }
Pt=1.25; %总压(atm)
$ J3 m, G8 I% W) z2 a" j- M0 O! r7 m
rw=0.025; %管径(m)
9 ~& Z5 x$ }) ?0 _. [- p( e& N0 k. j
Tw=100+273; %壁温(℃)
/ {+ K1 {8 R: A: t3 b- T
G=631; %质量流率(kg/m2hr)
, U- }; x0 e e# D, _
M=30;
: G( s. c% ^1 M E) k& j2 A$ Q" ]
y0=0.0323;
) ]" {; r! w' `- q. Y8 [) S
Mav=4.47;
2 x1 C" c/ U5 G# f
rho_B=1200;
$ r& F6 ]( d; J4 b6 u7 g4 L( Y( Z% D
Cp=1.74;
$ L% B* m) {2 M- x0 h4 E1 c
dHr=-49250;
$ d0 ~8 u" m! h% K. `
h0=65.8;
- d/ O. g3 I+ R4 t& }/ K, U
T0=125+273;
8 ]) |" S7 j& Q* t) d5 D
Lw=1;
- r, W7 Y; r2 l5 R( w1 i5 y
u=8.03;
$ ~* A9 x6 u: A; R7 ?
R=1.987;
% y4 ^4 r' U7 Y4 b3 I1 [0 P6 O7 I
ke=0.65;
6 A& q$ \+ e* l! I4 x4 m
hw=112;
0 X, M) p( Y3 P% O' M
De=0.755;
- [6 o% N$ B2 ^2 o/ D2 X
%********************
6 G' u6 e9 J) Z" m( c% S3 u* l1 b
m=1;
4 @ J2 g3 g/ l: ~
%********************
, h8 e' K6 }! d7 y6 y% g
% 取点
& y7 B7 H( x. E
%********************
1 w* a' T. k" y% \% e* j
r=linspace(0,rw,10);
% {% D5 T+ N- _% B
L=linspace(0,Lw,10);
# Q$ U1 ~; f+ p @( Q* }) E7 k
%***********************
( g0 @& x* K; q# A [: }
% 利用 pdepe 求解
" Q9 V& o+ ?5 V X u
%***********************
( l6 A* i5 v& U C9 r. W
sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L);
1 o7 U& b4 t, |
T=sol(:,:,1); %温度
# P. N2 D& M* L, O8 f7 T
f=sol(:,:,2); %反应率
0 }: [; I q" i
%***********************
. O4 i# k7 F, D" k' e- i& Z/ w
% 绘图输出
5 A# C7 S- O6 ?3 s
%***********************
% g4 k$ `! Y7 M
figure(1)
' w# ?3 z) _1 `
surf(L,r,T'-273)
0 e( K! l1 l V( Z0 u) t
title('temp')
( w5 k/ l% z5 i# p" F9 O
xlabel('L')
6 o/ @2 U4 @. m$ Q8 \3 A$ o
ylabel('r')
+ l' a" a- e1 o. K
zlabel('temp (0C)')
& h' `1 C3 z! ?) h( n
%
5 O) n! h& |. D! m0 u% {5 j
figure(2)
6 H n0 c# d) |+ I0 M7 m6 R
surf(L,r,f')
: d9 k9 E" e# a8 C& z0 y7 G
title('reaction rate')
5 P; s3 w0 I7 A
xlabel('L')
7 L$ N3 s5 F' S$ _, V# y" O# |
%初始条件函数
( D. S# z' S( M/ H7 x! o
%**********************************
$ F! z ?0 f( N6 c/ T
function u0=ex20_3_1ic(x)
: p$ u" o9 l0 Q7 x
u0=[125+273 0]';
; T6 F" S' L7 ~ h" Z
%**********************************
, J) C. P! J4 H S' j$ ?
% 边界条件档
" @( [0 m; |; s5 J! U3 a
%**********************************
3 L; G; o6 G9 R1 ?% q' u6 y
function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L)
; A1 x, s( t; a# J3 c9 V2 f& v
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
^1 k% G- j' [
pl=[0 0]';
; b9 ]2 x U6 _/ s4 M: d
ql=[1 1]';
) Y( t, h/ }% _/ v7 G% U. {2 w9 ?
pr=[hw*(ur(1)-Tw) 0]';
: b6 Q# Z/ f4 D3 t3 \: F
qr=[G*Cp 1]';
4 @0 `, m; F' j& U3 u4 \4 b( l! t
ylabel('r')
/ O' f3 A/ u3 M
zlabel('reaction rate')
* [( f0 W' z; [' z1 F8 X
%*************************************************
' |$ \& ~3 A4 X) M3 o8 J
% PDE 函数
1 d5 @! ?+ ]! P% n8 J' [
%*************************************************
, S% F0 h j+ w+ R3 m9 m" F
function [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)
H! T- B+ T. l9 ~6 \2 z9 g3 e6 Y
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
" C/ D, u) C- j9 e
T=u1(1);
# H' ]. @% r# [6 c6 |6 d: u3 J/ P
f=u1(2);
. F: I* \* x0 v3 \4 U
%
$ a" `2 k* c2 R1 F
k=exp(-12100/(R*T)+32.3/R);
4 p2 p5 {/ p/ r5 D4 w
Kh=exp(15500/(R*T)-31.9/R);
* @$ W% v. R+ Y' ~3 K0 X
Kb=exp(11200/(R*T)-23.1/R);
- m: H3 ~. T* M& `
Kc=exp(8900/(R*T)-19.4/R);
. s) w+ g& J; x" U: S
%
, i# i! T1 z9 d, [2 Z3 }9 X: z: F
a=1+M-3*f;
- o6 y# D' O1 y; P) P; I
ph=Pt*(M-3*f)/a;
3 H: {" l2 b c6 D9 I# J
pb=Pt*(1-f)/a;
$ w2 @7 I( N* I: V" s
pc=Pt*f/a;
* j* y4 [5 `, U/ l6 `1 |$ |2 N
%
; j/ ]# V/ a/ \5 j" j
rA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4;
' r* s+ v8 ^* x1 U6 F, I5 C" E
%
# k2 R$ C2 G* k |& a
c1=[1 1]';
$ q' u% v$ a' T5 \
f1=[ke/(G*Cp) De/u]'.*DuDr;
}; }; L7 l' w, H6 i
%s1=[ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw)
$ X2 [5 e: o1 g$ i
s1=[-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0)];
! N* b! b# B" k/ ~3 a
%**********************************
" r! W2 Y9 ] z# g8 A! _/ x
: C' ]3 d- T& H1 a
————————————————
1 D6 t8 D J( v7 J p
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
n( [8 _! |$ i8 i, K+ L
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711536
% y; {# w( r6 i: a3 O
. X4 U: p( ^: {3 E. ^1 _: ~
8 R5 D% s; R8 A% x; I' ^
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
- k. H1 m% a) @
$ p' x: t3 @; T3 x- m* m& z
# t8 Z) Y6 F6 ~
1 I6 b! A6 @& g' z
( m! ]" K0 q! z0 R