数学建模社区-数学中国
标题:
偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:27
标题:
偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布
例 4 触煤反应装置内温度及转换率的分布
P/ Y, U4 ?* Q* b. P- e/ b7 l5 N
, L% t* W& f2 g/ V4 B% I
以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下：
6 L% F; }; l7 R! c
7 j. Y5 y( Z& g' U' X6 T
. ?% k& ~- m8 I' J/ h3 ?# _1 o& l
' O3 |2 Q! P# @* B$ e" |
其中T 为温度(℃)， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为
* h @% S; ~* Q0 h% A% Q
+ {+ D+ x3 U6 ~+ R% ]5 e! j! _# O
$ b8 R9 U; H# T8 P# G0 ^
2 E7 e& ~( ^) s
此外，式中之相关数据及操作条件如下：
2 [* {% Y5 j9 w
1 S! h! g% l- \) L! r
(i)反应速率式
& |4 R7 i2 m' j# E
$ N8 m5 i3 m9 O7 j0 Q
4 |: [, S& t- a
9 d2 N- i- E r. t$ Q
其中 P 表示分压(atm)，而速率参数为
( I3 M0 r' g# r& K X5 t
5 L. u8 k7 C$ O2 }7 X6 Z2 Q5 h% ^
, A1 ^% R5 H* {* E7 ?- P3 x
_) V/ c- y' c9 o
上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。
* c4 K1 n# W! ]0 K5 {8 R: Z) q# O, `
+ i! R1 D6 x& D( U
(ii)操作条件及物性数据
, V; f3 U7 v e$ S! f9 ]2 X2 {( j* f
0 N6 r" K5 s! }) i$ ~
' [- W& z5 j/ D& n" h `% t
* u3 N5 D- C1 X. e
3 A$ Y7 x7 k/ v. H* C }
0 u- q0 f5 E a" D$ w1 |
题意解析：
3 M, S9 w7 \" q" f% A, g
. b6 j$ I3 K8 e7 {5 O% R
0 e$ K* }" {' _* J. J9 N4 H
; V, g6 h" m( \& D: ^# m
& \, \6 }5 I |7 l) x
将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。
7 F, n" G9 }) P3 f/ m
, ?+ ^- ~ L2 f/ O3 H
MATLAB 程序设计 将原方程改写成如式（35）的标准式
" m" E7 G, ]) G2 I
+ c& A$ _8 l) ]& Z7 d! ]* H$ P
& n% j2 o$ x8 `0 A& |
" Y2 ^+ x) T$ N; @; z0 v+ Q+ O
因此
: I# n1 O' o" X R1 }+ ~
" K6 U" x! \) b$ c
! R* W. p( ~1 V4 B$ s1 K
% ]5 p5 c2 [& j5 W* n
根据以上的分析，可编写 MATLAB 程序求解此 PDE 问题，其参考程序如下：
5 {- F6 K& B; s. g# L3 s
% _! M4 X4 O* M
function ex60_3_1
1 L- d, v5 N1 @# ^
%******************************
3 V5 X& y' `) l% w
% 触媒反应器内温度及转化率的分布
2 e7 z$ N/ o' {% h$ p; K) i
%******************************
) e/ `0 e3 R K. P7 K
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
8 a+ R2 ?$ a3 r
%******************************
3 w% \/ v4 \ ^
% 给定数据
# w) ?; @6 X% w4 @% r4 n
%******************************
. v2 b% A# t3 W* R4 \+ }4 G
Pt=1.25; %总压(atm)
! y4 s. U; l1 L4 D2 r {1 S
rw=0.025; %管径(m)
, M* Y# ~3 N7 o; r: B: U# s( E- `
Tw=100+273; %壁温(℃)
9 W0 T/ w# r2 P' f
G=631; %质量流率(kg/m2hr)
* O) y% C7 k6 Z3 S
M=30;
# ~$ H- g, B- |/ r) W! |
y0=0.0323;
/ H5 ~ w/ z. T+ t( g( w
Mav=4.47;
* _ U+ V* G" i# C% k- b
rho_B=1200;
$ o7 H- C/ B3 t/ M0 c0 S! o
Cp=1.74;
3 \& j) T+ M0 b9 K
dHr=-49250;
" M& l5 d q& y# ~: P8 D3 `
h0=65.8;
, F5 u8 i v. I) w* V9 ^% ~3 p
T0=125+273;
* @) W# p" L5 f+ f7 {
Lw=1;
2 m [* L3 {; w$ t3 c
u=8.03;
. ]# O$ z7 \3 f0 t
R=1.987;
& l9 T- \% t# q3 l
ke=0.65;
% U: O, H! Z4 S. S) h# ^
hw=112;
3 b! j! a( ` T, V8 x$ k9 a1 Q
De=0.755;
/ @4 `% G4 ~- z& ~; ^" l
%********************
- r7 x# e6 i7 b6 r" x
m=1;
; [1 C. A2 w9 O, O; y* K7 C: |
%********************
6 d; h6 h, O7 Z h7 U' h
% 取点
# ~6 H7 x3 L2 z- @& ^
%********************
+ g2 U; Q* W- d# Z1 G8 G2 y: e
r=linspace(0,rw,10);
4 a/ U; [0 A3 M1 ?( V. L# ~
L=linspace(0,Lw,10);
. H6 R- h9 L* s1 c& ]/ f8 r
%***********************
5 w" s4 N5 C$ z6 t
% 利用 pdepe 求解
5 ~; _. l& M9 I! k, c
%***********************
+ A1 D5 u2 Y) S; W5 K
sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L);
8 q7 P9 j6 J# U
T=sol(:,:,1); %温度
# T, u6 a! L" R6 n( \9 L
f=sol(:,:,2); %反应率
; K' ~6 h/ J8 N7 @' F# p
%***********************
; r! p& Y) w: T# O, V
% 绘图输出
% H: `; F3 S& Y& r( D
%***********************
' t- Y) J4 t: D( i! a o
figure(1)
% h! A1 B) f! a4 a4 W
surf(L,r,T'-273)
* {; u9 ?1 s4 M
title('temp')
5 @. H4 x; _. i0 }6 F8 a4 Y+ g
xlabel('L')
6 O; ~- B. x \$ d
ylabel('r')
1 @! t' J( p p; `. B. V
zlabel('temp (0C)')
! j' w3 H) K4 R! r7 y, {* Q I
%
1 o( w# d+ a; \( Z) e% j
figure(2)
( l) a( q/ n- H7 O' s/ k4 l
surf(L,r,f')
5 ~9 d) V6 X; d
title('reaction rate')
4 X* N# }3 U2 o% M$ j) ?0 }4 c
xlabel('L')
) ]. u' P A9 L! ?- Z
%初始条件函数
) T3 a4 {3 s. s
%**********************************
& ?$ z6 a( m* J5 H- ?+ e
function u0=ex20_3_1ic(x)
5 u2 _& l6 X3 h+ j7 x
u0=[125+273 0]';
, Q6 h8 G4 s4 s& G
%**********************************
( y1 d G# y! x8 j* [$ M
% 边界条件档
" D. o2 _' B( z
%**********************************
- A% N' W+ t" N9 q) i. v
function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L)
3 ]4 \ Q2 M0 w- U" T
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
2 G. k. P: n* s& ]1 e
pl=[0 0]';
/ i: v2 L. M; i# U5 Z1 \" _
ql=[1 1]';
7 z: ~+ U, d z `5 `* W e5 Z
pr=[hw*(ur(1)-Tw) 0]';
0 a0 p2 v) p# ?% M0 A# a
qr=[G*Cp 1]';
" T) ]" b* M4 N+ Z
ylabel('r')
2 p& k( Q3 s" K/ K A' |' B& e8 X
zlabel('reaction rate')
; P- T8 F; u) e( z' L0 ?3 E
%*************************************************
) N$ E) e6 R2 V8 f
% PDE 函数
* t" R; }; E1 k# F2 O/ ~! ]; | \
%*************************************************
6 E8 k+ j/ P" |, b$ ^
function [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)
) g+ y7 l7 a7 N# ?0 { c1 q/ L
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
?' }# e# ?6 A
T=u1(1);
' r2 G8 a0 `* E
f=u1(2);
4 N- [# X6 p. @( s9 k; e* ~% L2 m
%
& N1 I- q5 H8 ^/ I" z
k=exp(-12100/(R*T)+32.3/R);
. B! O# e+ i$ z. K( F/ X
Kh=exp(15500/(R*T)-31.9/R);
?& X4 r# X, N/ s
Kb=exp(11200/(R*T)-23.1/R);
9 i+ R# ^5 h" |: W" V4 p" l
Kc=exp(8900/(R*T)-19.4/R);
7 s- n H2 w% o# O7 a% |
%
" V; A6 R) N: j
a=1+M-3*f;
8 @) ?( f* q0 z* G" ]
ph=Pt*(M-3*f)/a;
5 L& E4 s; O$ D
pb=Pt*(1-f)/a;
" T+ i' N4 h, A3 M6 y' c
pc=Pt*f/a;
2 f# s* M! V! |
%
& {9 G& D- q; Y' J
rA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4;
4 |3 F8 L' D3 S0 N0 C
%
! T8 D2 g: y: r+ B) z3 E* l% A
c1=[1 1]';
2 | C2 ?8 R/ x- O S
f1=[ke/(G*Cp) De/u]'.*DuDr;
. e# T8 e- t1 ~9 X: y" ~" t
%s1=[ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw)
" B8 E# ]) N3 N- Y3 B1 m' K. B3 R
s1=[-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0)];
( v. t, v( c( }. z0 V4 ~ \8 f: I
%**********************************
' A6 Y/ X/ ]4 { G3 `
4 B8 `0 _% R, E [" Z3 l3 v0 F
————————————————
0 V3 ]' u' k* D7 R
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
5 y* g u" j8 A
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711536
/ M6 i! Z% v! [& n. W
( {0 T4 B6 X. n6 l' h {
) \3 i7 g4 ~. C1 Q7 Z( f! w8 W, F3 c
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
. ?% k& ~- m8 I' J/ h3 ?# _1 o& l
4 |: [, S& t- a
& n% j2 o$ x8 `0 A& |