QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 2901|回复: 0
打印 上一主题 下一主题

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

[复制链接]
字体大小: 正常 放大
浅夏110 实名认证       

542

主题

15

听众

1万

积分

  • TA的每日心情
    开心
    2020-11-14 17:15
  • 签到天数: 74 天

    [LV.6]常住居民II

    邮箱绑定达人

    群组2019美赛冲刺课程

    群组站长地区赛培训

    群组2019考研数学 桃子老师

    群组2018教师培训(呼伦贝

    群组2019考研数学 站长系列

    跳转到指定楼层
    1#
    发表于 2020-6-10 10:27 |只看该作者 |正序浏览
    |招呼Ta 关注Ta |邮箱已经成功绑定
    例 4 触煤反应装置内温度及转换率的分布* {! P5 [$ a5 Q3 w$ b; d/ S" _; E

    . D3 |( p% ^4 G( D以外部热交换式的管形固定层触煤反应装置,进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下:& [  D' |* r' ^3 {2 E- t, b. X7 Q+ U

    , D' C% C  M( a4 O4 j- }2 {5 ~" |9 B1 ~

    * q* m. j$ l* r9 \, Q3 w8 i- [5 ~ 其中T 为温度(℃), f 为反应率,L 为轴向距离,r 为径向距离。此系统的边界条件为; I7 K/ G9 d: T& B9 q$ S, `& i

    . F. U; X9 p1 X5 [" P! d
    * a  E& }6 ^& Q) F& z
    , ^, r/ ~* D; G0 F此外,式中之相关数据及操作条件如下:: X: {4 e! \& x; b3 {/ A+ x! {

    , k$ {: O& x* u; P4 v; K(i)反应速率式
    $ n$ f. x& O4 }( V: v+ E5 O; c7 ^3 V

    : a" ~! O$ V, H
    ; w& Q, L; B1 `6 q+ H8 @8 N2 t1 m其中 P 表示分压(atm),而速率参数为
    3 Q; `) ^5 G; |& h2 c& N! T
    9 ]& Y# P! T& [% R# y/ a1 z- T
    / ^* n% W' K& w2 u0 g
    % M* |% r4 ?" h1 r: e1 p$ l( n上式中,下标 B,H 及 C 分别代表苯,氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。
    . X) X0 s" ~  Z8 `  j: \: K) T4 i7 Z  Q* V. _' ^
    (ii)操作条件及物性数据
    : N2 h- |- O2 }8 e  G/ Q4 b' m
    4 `6 s: _+ b& }* z1 t
      @" K$ r& \' U$ v! t2 o. W: h' P9 T2 g6 A! }& ]

    7 G6 W( `& `# y4 o, W1 u8 p9 k" Z# ^
    题意解析:
    . M1 }" H: a% b; U8 D# a
    2 U0 Q( p) z; I. Y* Y8 z3 _2 n( x) C- I, J) F
    / ?2 b0 z$ L5 Y8 U+ f* K
    , a# L7 V6 r* F! ?2 a4 m7 r, s0 N4 C
    将上式,连同反应速率式,带入平衡方程式中,配合边界条件,可利用 pdepe 求解。
    $ I( U7 M8 L9 {; q  w4 |6 O7 b/ Y9 u. F+ S2 i7 ~, Y% K
    MATLAB 程序设计 将原方程改写成如式(35)的标准式
    # \. G5 k0 I" S0 X, A
    % m" f) x: V: \* k$ Y6 ]8 I$ J) v, F$ L; x/ ]% J8 |
    0 G' i% g) P& l- R9 k1 q+ F
         因此
    * p* d& W; X3 e: `( h1 p/ D; ]/ N
    3 D9 ~& ~/ q, M! `* q! c# s- Y7 D, s+ W) T3 |: ~) Q% ^3 p
    0 `/ d1 P  o8 v' O7 Z+ Q
    根据以上的分析,可编写 MATLAB 程序求解此 PDE 问题,其参考程序如下:
    . f# [0 C8 s% K9 q2 p) A. E5 e$ D0 x
    function ex60_3_1/ Z9 t6 A. I/ m) ?  K
    %******************************' i7 s* b% L+ B/ {) ?" F
    % 触媒反应器内温度及转化率的分布; |7 ?9 Q2 B$ _- B; X* ~
    %******************************
    & [7 O# n+ F" Q  a/ h$ Cglobal Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
    : X4 @! a3 V3 m8 h%******************************3 o4 p+ }9 D- i3 G+ U
    % 给定数据9 b/ @, m: k) v% X. B' d
    %******************************
    ! U  Z/ l! g, S# H: p* UPt=1.25; %总压(atm)8 t* a0 A" I4 A' P1 d& s1 |! Q
    rw=0.025; %管径(m)
    ; ^- W  X5 p/ K1 n+ fTw=100+273; %壁温(℃)
    & H: `, P7 V3 U' jG=631; %质量流率(kg/m2hr)  N% o0 F0 U- O; D  m
    M=30;
    " b0 O( u$ M& O" Ny0=0.0323;, Y+ N: o. V1 n( F7 W8 Q
    Mav=4.47;. i0 T# f: \3 ^( x  \
    rho_B=1200;
    4 @) q' X) O. r- T6 [! E( {Cp=1.74;' o$ f- H+ W( t6 ~+ J/ y
    dHr=-49250;
    " H* g9 n5 T6 _) e$ Jh0=65.8;
      V$ U( u5 y2 @; QT0=125+273;/ X& [+ @& U' E, s
    Lw=1;+ m- w& M- l# i3 U* V6 T, G4 R3 e! z. ^
    u=8.03;
      z8 z& c2 `* T$ H- |R=1.987;9 C  ]$ z- D+ I7 @+ F7 C1 t& `
    ke=0.65;
    & `8 N# A8 t6 V' }! j' thw=112;& H$ H$ r9 u2 C6 H* u* U+ O/ B
    De=0.755;
    0 p2 H; [8 V+ d- i0 T% [%******************** - P* x. d; f3 w- J! F5 i
    m=1;
      p: S) ?5 z. w1 N%********************% I! w5 r& N1 }0 ^$ V
    % 取点
    ) M0 @: Y" _  b( g%********************
    0 l# ]5 V: g+ g; k& c% A$ |9 C/ K+ Pr=linspace(0,rw,10);
    1 P) M* G; i; y+ b. W- ~/ O3 ]L=linspace(0,Lw,10);
    - o4 h! F9 o' ]# u7 V" i%************************ C) \; s0 U, Y$ T$ L+ @" c7 y
    % 利用 pdepe 求解5 t9 [  @& N! L6 q$ R* `
    %***********************
    $ g2 j: B0 z1 ^sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L);
    + W+ o8 s# R: n3 @T=sol(:,:,1); %温度
    1 [! `  x& U+ ~f=sol(:,:,2); %反应率
    0 }* U; f$ m+ w/ G0 I) [%***********************" G  k, q3 n% S& U
    % 绘图输出
    + j# G; q% _+ B/ Q%***********************& [. D! d2 }0 ~% }
    figure(1)4 w" \5 i. Z. |" I
    surf(L,r,T'-273)8 e# l8 `* A) g3 X+ f
    title('temp'), `1 Q4 X2 q! e( B
    xlabel('L')0 Z( y7 O8 B! [6 m( U
    ylabel('r')
    7 l0 ?' X7 [: p9 H" M' M6 zzlabel('temp (0C)')
    6 W, g* |: X; e0 h%
    ! W  {0 H3 F2 p+ ?6 k+ Nfigure(2)
    7 C0 Z" B8 V" [surf(L,r,f')+ l/ m  y) b7 C, s! c( l( V
    title('reaction rate')
    7 a% b% p( S9 c7 |  ixlabel('L')5 `; ?8 s8 s3 z: V: s" S, y
    %初始条件函数* R5 D* u. N9 P; r+ o6 A1 K
    %**********************************
    + F) ?3 i  t% n/ a8 ^function u0=ex20_3_1ic(x)
    6 O; r0 e' i4 ?' V: ~u0=[125+273 0]';; o# u. H5 \( [* R3 ~* s& l# L
    %**********************************
      f( O( J, j1 f# e; F2 [% 边界条件档
    4 L& ]* [, f2 Q9 `- A3 A. a%**********************************4 Q6 p$ y4 e, s/ L& E6 S
    function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L)
    5 u+ r. n4 v) p. v, r/ v7 W8 dglobal Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
    1 |  O$ k0 X, j; Y2 L5 ^  i+ }pl=[0 0]';
      {. w) Q& n( G8 y* Jql=[1 1]';5 g  q' M2 }% x  y
    pr=[hw*(ur(1)-Tw) 0]';
    2 ~# N' n. V8 X* Rqr=[G*Cp 1]';
    # [( r, I7 b- P% M( g, c  v- Rylabel('r')/ }! c* `' F1 W1 ^
    zlabel('reaction rate'): Z( O+ C# `- t  y
    %*************************************************, S4 o# n3 M5 F3 B, I
    % PDE 函数% h. }( |# r( y+ v) Z
    %*************************************************
    1 X0 [8 r+ ^/ w3 p' n* O: ifunction [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)9 w) |" J9 e7 Q$ @$ Z; ]$ x
    global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
    + n- [8 j& w$ r# e' `. o! b5 B$ j6 _T=u1(1);/ [) A+ E, ^5 R6 }( f
    f=u1(2);
    / {5 N  `9 l  s, a3 }0 v%/ L9 x  W& d) |
    k=exp(-12100/(R*T)+32.3/R);% n# y( }: l, M# w% t# l  n
    Kh=exp(15500/(R*T)-31.9/R);/ U8 {; q, d6 p8 b6 O- Q
    Kb=exp(11200/(R*T)-23.1/R);# ~! N. k: Y/ }
    Kc=exp(8900/(R*T)-19.4/R);
    $ z& X( t4 V+ J  l# l0 r%$ c% o( Q% r  k/ O
    a=1+M-3*f;
    0 r5 B3 v. ]' r. Z1 ~- ^7 C/ sph=Pt*(M-3*f)/a;
    / E, v% p, F! P3 {! \' Vpb=Pt*(1-f)/a;
    0 i% [9 N2 w( Y3 q- z3 e9 _pc=Pt*f/a;7 a$ v$ K0 P- D: W2 i$ M& s
    %6 f7 O: [$ C# [5 L9 I1 H
    rA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4;
    ) y( m: B: v7 x3 O0 o%3 E" z# A( x$ S3 m# @5 v3 v( \
    c1=[1 1]';% H9 D. O: f3 a4 T# A
    f1=[ke/(G*Cp) De/u]'.*DuDr;
    $ Y) d& `! e" ^* h4 f: B4 j%s1=[ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw)2 y9 p9 }) n- ]! o5 J
    s1=[-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0)];
    2 b$ ^4 \. k- Y7 T7 o. D; v%**********************************
    " A; Y- ^& u4 s7 Q% D3 r/ s9 L5 R7 Q1 I8 p+ v  i
    ————————————————& L) M0 p2 M9 ^) M% }! Q4 O
    版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
    5 `( @( B. y! B% a原文链接:https://blog.csdn.net/qq_29831163/article/details/89711536
    * ^1 @* o9 |$ [  Y3 U' b8 W2 I/ U# p  ~/ @0 K

    ; Z0 \/ p' W" _& O! U4 A7 G
    zan
    转播转播0 分享淘帖0 分享分享0 收藏收藏0 支持支持0 反对反对0 微信微信
    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2026-9-13 20:28 , Processed in 0.535447 second(s), 51 queries .

    回顶部