数学建模社区-数学中国
标题:
偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:27
标题:
偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布
例 4 触煤反应装置内温度及转换率的分布
9 j1 s' I, v' I( b
2 M6 E0 L3 B4 s; O5 T
以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下：
9 N7 _; A8 a0 g5 k1 f8 v3 q
* {8 N; }' J+ B
" S; O: s. ^8 W, n) u4 l# N
! K) M! Z) E- d5 C: z
其中T 为温度(℃)， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为
2 ~9 Z2 e# G8 W" }( m- U9 C8 M
! Z' c+ P5 n- K5 ?) S. m# T# h
) y L {9 ]9 G4 n% Y9 V
+ Q' Z1 u A( H7 w: T2 {; G
此外，式中之相关数据及操作条件如下：
: o) S0 K3 ^: a4 R
9 H" o, c. k. I5 l4 E+ W
(i)反应速率式
" W+ c! w% m; W$ q* e: G" i
; Q y/ i# p0 A0 ~- u6 Y3 m
1 W9 A3 a, ?) H9 W
" `0 l3 X" a. q* K2 y% g4 I7 v
其中 P 表示分压(atm)，而速率参数为
# y( u- T" i5 u9 m
8 u% q+ L+ Y6 U' d/ ^
( M! N/ J9 t& J% t! R* M& }3 O b
9 @! R' J6 f: B
上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。
4 q" w; E/ Z- d% \
0 i' I A; w. y, }4 K0 t8 I
(ii)操作条件及物性数据
$ a) K: [2 O5 }! J, {
4 Z2 K/ A' [2 R
# S X2 N! t) q8 p" _1 ~
$ l: |. ]0 x0 g& Q+ o6 V0 x& M7 G
' y* p; {" m* C: b. y2 K. }" L: \% @
' y O5 ?1 o9 `2 a
题意解析：
( x! c3 c( v1 O7 [9 {/ o
9 l0 f G/ d& Q$ X" M& o! r
" e1 b8 X+ }2 [/ m5 f' }
0 q& t, q* P! d' c
4 Z/ d; B7 q/ V& m1 s, ]. d
将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。
{& n5 `; Z# Y; J) s
/ B) y' A. j0 o# O. Y. D
MATLAB 程序设计 将原方程改写成如式（35）的标准式
2 y7 f1 D. S9 `1 q* k, y. {
) G, Y* x9 U2 T2 E. q# P
! q6 o9 |. n$ t5 m o5 e6 x) S8 H. n
* P" C# T6 ~4 b4 ]& W; |
因此
4 x! _' o( \& d% q' X( u% H, \7 `
. ]. E! H% N; J$ e/ ^
# F# W! S5 ^. r2 M0 P
9 e$ L' E/ a- s: `2 f: ]
根据以上的分析，可编写 MATLAB 程序求解此 PDE 问题，其参考程序如下：
3 l* t6 _% ^4 x# ? E
" A$ a4 C8 t( [/ X1 p$ A, r1 V" u
function ex60_3_1
+ _# _2 N( \: s# }0 y3 l4 o0 O
%******************************
W3 M" g6 P, E& I7 e* i3 c9 e& o; d
% 触媒反应器内温度及转化率的分布
$ N4 b, n! h# F, ^
%******************************
H% _0 `3 E8 k2 y9 c7 t9 p5 P4 J |
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
" O( a/ p) g# b8 {7 u
%******************************
& m6 P) K; Q1 q, V% u. f
% 给定数据
6 R# g- g7 G& Y3 H$ v- k0 v
%******************************
4 a' ]3 Y& ]- |( j
Pt=1.25; %总压(atm)
9 T6 l' `" [& V [+ n
rw=0.025; %管径(m)
9 t, j0 e/ N' U/ b0 z$ W/ l
Tw=100+273; %壁温(℃)
L9 Y6 [1 J2 A/ g* S2 i3 P9 Y) P* C
G=631; %质量流率(kg/m2hr)
3 g* Z' r+ E2 Y3 X! F) s
M=30;
, o# Y b, `* F( @( n
y0=0.0323;
5 M9 p. g0 J; _8 B2 `
Mav=4.47;
$ ]- E# v4 `+ J) T' f
rho_B=1200;
9 Q8 C0 A5 P0 Z4 _2 w
Cp=1.74;
' u1 V- s+ }5 {( E. S1 q$ d
dHr=-49250;
" R1 H7 m+ f) j/ |- _5 j4 V6 W
h0=65.8;
/ d M9 I) ^/ q. ?* D
T0=125+273;
- |) v6 H! S# R' o
Lw=1;
r& I! t, Q5 w: B; x7 Y& y, d8 R
u=8.03;
; m% M3 f7 Y* H: L
R=1.987;
0 ]# b$ O% \& q4 r$ ^& O5 a+ y0 e
ke=0.65;
2 J. `* _, L7 Z# b
hw=112;
' h3 `& c$ R* n; A, v" ~, _ `
De=0.755;
! ]! W+ h6 N2 B' Y- w
%********************
' K& ? R' R* l0 }4 H& X7 q9 E
m=1;
4 _: S7 w" s2 h, u# W! z$ V
%********************
% W! M) V0 ]1 r: v6 v/ n. u/ k
% 取点
# y4 R: ] J( s- u2 o A2 ~ V
%********************
x' A, N9 F+ u% e% |7 Y4 t
r=linspace(0,rw,10);
% y4 v% w6 L+ R
L=linspace(0,Lw,10);
5 P9 {" i: i' x! c
%***********************
& L- a) t* d+ _* g
% 利用 pdepe 求解
4 t$ \: y S/ U: d6 X) _
%***********************
. n! W! n( j7 W7 i+ h; ?, ]$ p
sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L);
: @9 l I# p; e, U2 L3 d0 y9 C' P
T=sol(:,:,1); %温度
, F( {( ?( u% |8 D `1 J& c% i4 H! Z
f=sol(:,:,2); %反应率
1 ? P! p9 ^6 ?& k! `) l* k7 I8 Q* k
%***********************
8 C- j: t, b. ^, l( L; P6 K
% 绘图输出
: w! X2 F/ s2 R; U9 L' ^
%***********************
$ ^9 x; r1 x! @+ @7 Q3 [+ J9 X
figure(1)
- h# h) q1 f1 v* {7 J
surf(L,r,T'-273)
& d1 G% f3 _6 Q7 M$ x' e
title('temp')
8 f- U# G/ s ]' Y/ `( ~" M
xlabel('L')
" _+ _# P+ a# D; c8 G7 m) _
ylabel('r')
* J: ?- i0 g* g- {; n
zlabel('temp (0C)')
* H: p5 u3 Y9 t( J% {- ^
%
. ]9 L9 O5 H! i
figure(2)
* ?+ C" j6 h$ ^+ J; t
surf(L,r,f')
2 q; `; h2 N1 V) M/ m% U
title('reaction rate')
/ D X$ x* j/ ~. E4 j8 D
xlabel('L')
E: ], U# j1 x& ^; N0 B; u5 Y
%初始条件函数
5 h$ d+ m* H8 O4 h
%**********************************
/ O. [1 F3 _4 z+ A$ E$ z# f- E; A8 ^! h
function u0=ex20_3_1ic(x)
1 _* I7 y$ t) v) C
u0=[125+273 0]';
0 |0 a5 u) K3 H! a; d9 h. r
%**********************************
/ l2 I9 }! A1 a) I, M
% 边界条件档
' J+ ~- p% V! }7 C; W; `( r5 m
%**********************************
& R. m2 k" f1 \! }0 E
function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L)
+ E" ~# |/ d( f) V4 d9 U
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
& {) j9 [) x' X2 l9 A9 d
pl=[0 0]';
+ U+ d e3 f7 W& f
ql=[1 1]';
$ S7 c$ y4 G' D; R
pr=[hw*(ur(1)-Tw) 0]';
: l& p* V9 F1 R; n, x0 a( H9 y
qr=[G*Cp 1]';
7 v$ w& w$ E; ]
ylabel('r')
% C3 w) Z, ^) ~8 q- }( W
zlabel('reaction rate')
# M4 t2 P' I9 i6 `9 F" B; ?
%*************************************************
U! v/ |4 t' c3 R+ U1 v3 |1 Z4 H6 w
% PDE 函数
3 y( H7 H1 Y8 j
%*************************************************
' n$ w8 p& D0 \, k5 m
function [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)
9 x! f9 s- Q8 ]
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
* R" F# _ G1 K4 T, U7 z
T=u1(1);
" I- K% Q9 l% z5 ~+ M
f=u1(2);
; j$ j2 k# r8 ?0 q
%
6 E7 o4 s7 f2 _
k=exp(-12100/(R*T)+32.3/R);
# v6 n9 n& e/ R# T$ S# W
Kh=exp(15500/(R*T)-31.9/R);
" K" i( b1 e7 v$ w
Kb=exp(11200/(R*T)-23.1/R);
! q; m; u1 K! K. @9 e% L
Kc=exp(8900/(R*T)-19.4/R);
3 Q2 t) K" k, M
%
& A% C! L z& {1 U6 J
a=1+M-3*f;
" K" r' q" x# H- b# I6 }8 l1 D9 @ D
ph=Pt*(M-3*f)/a;
" H" h5 p( w: E8 c* Z
pb=Pt*(1-f)/a;
8 y* _- F" F1 a* A" u V1 K
pc=Pt*f/a;
% ]+ A+ _. p. q$ i5 V. A$ ~$ V3 d
%
% u9 C. p/ M: V; i! d4 M8 d( _9 o* X
rA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4;
4 k& V$ M& l% R8 m: ~, g& a% q
%
9 l. N. Q$ e. D" }/ ?4 J
c1=[1 1]';
/ ~. J! S( F. r" s% I( N% J4 M
f1=[ke/(G*Cp) De/u]'.*DuDr;
* K8 k* J; r% e; x# d
%s1=[ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw)
. [) E" P" I; z3 i0 v
s1=[-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0)];
# R( g- [# w1 W9 q) c
%**********************************
$ ~. C% J8 u+ `. T
) F: J6 W8 ?7 V1 D8 H
————————————————
/ t9 c ]- F3 _2 q4 y! z5 z7 ?
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
, H1 G& f4 B; a3 o. L' o# `3 C
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711536
5 K( C# L4 D: J; L0 S7 R0 w& z) I
& }" G" n; q: P, P5 W, d
8 M0 d& j/ s+ p" ~! v' M6 V
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
) y L {9 ]9 G4 n% Y9 V
1 W9 A3 a, ?) H9 W
( M! N/ J9 t& J% t! R* M& }3 O b
# S X2 N! t) q8 p" _1 ~
" e1 b8 X+ }2 [/ m5 f' }
! q6 o9 |. n$ t5 m o5 e6 x) S8 H. n