QQ登录

只需要一步,快速开始

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

[建模教程] 偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

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

542

主题

15

听众

1万

积分

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

    [LV.6]常住居民II

    邮箱绑定达人

    群组2019美赛冲刺课程

    群组站长地区赛培训

    群组2019考研数学 桃子老师

    群组2018教师培训(呼伦贝

    群组2019考研数学 站长系列

    跳转到指定楼层
    1#
    发表于 2020-6-10 10:29 |只看该作者 |倒序浏览
    |招呼Ta 关注Ta |邮箱已经成功绑定
    # p9 h6 W3 f6 a! i& E% p; K
    题意解析:5 i* e: q9 x3 c9 H5 p

    $ @. a6 {  r9 l5 r4 C(a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:
    * s- r3 [. z0 h/ Z
    - ]8 D( {# \* s* {
    ) q; a. F" O' q* J0 T: V) y$ }( i: ?) U  h' {' e" t; a
    (b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为" C5 B4 Q; `# }
    1 m2 R. k# m6 p9 S4 }" o0 b# Y6 N

    6 `3 d) p' d" w$ a
    ) s2 [/ e5 C0 v' h9 F. x而起始及边界条件同上。
    % A: r% Q/ {) m2 T  X! a* D  \6 R9 E- b; K* k$ X0 e
    在获得浓度分布后,即可以 Fick’s law
    % v% c* X2 N" `, p" a. q5 B7 d7 e0 {( C8 S( E. R6 V

    & O) Q+ F% p. l, I: z
      x* Y3 n+ v# ]/ Y3 `+ i计算流通量。
    & w' {5 H5 [' w% H; F# J, V3 h' S! c+ {$ _' a2 V
    MATLAB 程序设计: 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
    + b" K3 ]; z+ v7 P9 b9 A( l- i- p- z
    - ^  E8 ]! j6 E7 n) o; ]) K
    3 M$ N% N$ Q9 V8 u; U) N
    利用以上的处理结果,可编写 MATLAB 参考程序如下:
    + V$ \- {4 e2 u6 j
    / B% }1 \$ T8 q- y2 _$ Qfunction ex20_3_2* ~) l2 f+ v( \: @! |
    %*****************************4 L( m3 X8 [$ J3 o7 M
    % 扩散系统之浓度分布2 W' r. R7 e' }# R
    %*****************************
    3 H# N5 @# ^% y6 r+ w; \clear
    1 K* Y- {, j2 ?6 o' `8 m6 ~! O+ aclc! M0 z3 w. ^5 w& S
    global DAB k CA0! d7 d) ?3 E$ k% `
    %******************************
    5 o9 U2 j: F7 p9 Q9 K: ~$ v% 给定数据
    7 m$ e1 @' ^) o4 u) o%******************************+ y* E4 R9 _( w+ f. k5 K  x
    CA0=0.01;
    $ R6 l8 }0 ]* |: w% |  g, FL=0.1;
    ) [; q0 ]7 u; Y  q9 E: w6 DDAB=2e-9;
    2 F( B* l* F* j$ \( `' [+ @k=2e-7;) w# U( {/ J* T2 r
    h=10*24*3600;! F4 T; L- n* v, R6 |, f' S
    %*******************************  f5 \8 N  L8 d# A
    % 取点- l% p1 G3 k, z* Q. ?8 H
    %*******************************$ j- W7 w, j# l0 Z" r
    t=linspace(0,h,100);
    , u$ b2 [& v9 g( r& i! iz=linspace(0,L,10);
    $ ^- g3 ]) J  O7 U: J: R* b%*******************************
    1 ]  k+ a. B" z% case (a)
    / t/ {, x) _7 u2 _8 B$ e%*******************************0 f' y, I  |8 Z: S* F  k
    m=0;
    , b7 ?( |: E5 i+ k% Jsol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);& E6 n; j8 f" Q3 |$ c; L
    CA=sol(:,:,1);7 ~5 C3 k' t6 s0 m) [
    for i=1:length(t)  Z2 C: V" C% E6 l
    [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
    0 m9 f8 p3 y* v) K NAz(i)=-dCAdz_i*DAB;0 @" A& |; ?% M4 b9 a( Y8 d
    end
    9 O; |; L0 Z& m+ S9 ^3 ~% mfigure(1)5 K7 J' \( j" O
    subplot(211); g$ d6 {& r$ d- }+ H) F- O
    surf(z,t/(24*3600),CA)/ V! E' v$ N- Y2 q2 t7 @$ O$ I
    title('case (a)')   v5 |! m: L# H$ r1 F/ A5 \
    xlabel('length (m)')9 K9 ~) U+ q7 ~- ~
    ylabel('time (day)')
    $ D" A% ~: p( b9 |0 Dzlabel('conc. (mol/m^3)')
    ' G- L5 ^. S( C: J) r: Xsubplot(212)) @6 \- A6 e" a# F$ w; M5 r
    plot(t/(24*3600),NAz'*24*3600), ]( t$ g; |. h9 x5 H* B
    xlabel('time (day)')# x5 a/ Y% q; A0 w' Z- b
    ylabel('flux (mol/m^2.day)'). e! }( i3 c/ j- z
    %************************************
    0 a( j" B3 U0 R$ r% case (b)
    . F$ a& g+ t0 z6 b# H7 _%************************************
    4 f; j  _: t( q6 T6 E/ U0 O% F* Fm=0;
    $ B% r- ?% P4 C& N6 [' vsol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
    ( a! g0 R) F& A; o* r$ hCA=sol(:,:,1);7 m% t) }* t) J9 [6 r$ k
    for i=1:length(t)3 R4 v* H7 Y0 T4 X0 l
    [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
    6 v' k0 B. ?6 f9 u4 P5 u1 F' [ NAz(i)=-dCAdz_i*DAB;
    ( }5 Z* u6 E8 S. X: }1 z, I! Iend5 z# c: y9 V8 }" a1 @+ v* r
    %0 L, r6 Y3 {3 g5 X. i) v
    figure(2)" [8 m% U( m  g0 X1 ^9 R4 _
    subplot(211)
    8 P; n2 o, b2 ~/ s1 J& ]6 dsurf(z,t/(24*3600),CA)
    : |7 x+ e; i4 Gtitle('case (b)')
    * e! H, U& j9 C* Bxlabel('length (m)')
    . Q$ _+ L4 B0 b0 o' Zylabel('time (day)')
    ) {$ v2 I- I6 l6 kzlabel('conc. (mol/m^3)')
    6 j  c6 F" ]4 g: m; \& b( ksubplot(212)
    ) x! v9 h) C4 X% z4 ]plot(t/(24*3600),NAz'*24*3600)
    ) N4 i; Y8 D9 q$ j7 ]2 Xxlabel('time (day)')
    ( j3 |0 L( D. ^, Iylabel('flux (mol/m^2.day)')
    : n, [+ x7 l& W1 @/ M+ F%********************************************
    2 K2 F9 w) C: m% M" }( Y+ }% PDE 函数) L% ^# s; ?1 w
    %********************************************, C2 h' o( ?4 y0 ^7 E* `3 o
    % case (a)& ~4 j! i5 _8 M- @9 o* A! o
    %********************************************
    / \# M' C' q, M# ]5 {9 k# i& Ifunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
    7 T% E$ m7 A% I7 O( T* l3 Wglobal DAB k CA0
    ( t! @% Y; k  b* |! b  R# tc=1;; f; F3 p% z' b3 S1 V
    f=DAB*dCAdz;! S8 G+ V( F8 u4 h2 t
    s=0;$ ?  Q' \1 [: H
    %*********************************************3 e& n/ O& O) c" }, ]2 X) l0 u: U
    % case (a)5 k$ {2 b$ i7 `1 j. }: p
    %*********************************************
    - N# g7 i$ q' _7 }function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
    ! T$ H' V3 n+ m# ~4 Iglobal DAB k CA0# n% \* P/ u0 S* j
    c=1;
    ( w- S1 N2 y- d. Z# j$ `. X, [f=DAB*dCAdz;
    ! E# w+ @, J0 n: h1 S& C: bs=k*CA;& x: U) G. X, o2 Y
    %**********************************************
    & k3 j8 l9 _1 p% 初始条件函数+ ]( [- K: U8 t' m3 t
    %**********************************************! v* U# z; K$ U" C2 k2 k# d% `7 O$ n
    function CA_i=ex20_3_2ic(z)$ w4 C' _: w0 _5 L: z" R
    CA_i=0;# O+ N  S3 Q, r6 @7 B, _: d, ^( T; M
    %************************************************ 5 A3 r( R( P/ g5 U! H0 O0 F
    % 边界条件函数  I2 c+ \5 ?- x3 M& F
    %************************************************
    3 g0 S7 D: p' `  ^$ hfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
    8 x! d/ N  p" ~( A2 ^global DAB k CA05 f2 O9 B" S6 G4 s- a
    pl=CAl-CA0;1 {! f$ K" B/ [. E
    ql=0;
    : w5 N; F* @* |7 dpr=0;
    2 G: n" H* C2 e4 v9 q3 aqr=1/DAB;
      ~9 s; K4 T/ k9 W" L4 X" q5 s3 t/ ~/ v& E5 w  u0 }9 p
    ————————————————
    / u8 w% l6 S! j8 Z# o8 O5 A版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。3 o, u% [7 w: x) @9 q1 L
    原文链接:https://blog.csdn.net/qq_29831163/article/details/89711694& k# B! B! ^1 ^1 X1 t( I& l" R
    : v- r) W3 w( v; l1 u. x
    ) u; V2 A2 w9 l/ t3 j
    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-6-14 08:10 , Processed in 0.574102 second(s), 51 queries .

    回顶部