QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 2558|回复: 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 |邮箱已经成功绑定

    4 G/ r% f' @3 P# E. o+ ]/ L1 m题意解析:! D5 A! t; O. p, w6 M
    8 u9 A& y/ F  B3 c# K
    (a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:
    # a( e7 e# R5 {( E8 W7 F/ \9 u. [( r  {

    * Y+ G. m3 U+ P! k
    5 |6 Y9 g2 g( f5 }( |(b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为
    ' F  x/ l/ G; R) R& x
    9 P6 X0 k5 K; }3 t. y% R- ]" P6 Q( O' q+ y; t3 k3 V5 O2 \
    9 T/ \2 y2 ~1 Z% ^+ v
    而起始及边界条件同上。  C7 ^; |: G* Y; K2 ^

    2 Q5 L9 |, N& ~- e# L# o在获得浓度分布后,即可以 Fick’s law
    : Z+ U' C# T$ r4 [& Y
    ' Y6 }3 ^9 C: A
    # d& I. L( y1 c$ s9 b! t+ @! H
    9 o5 o2 ?% P) H& Q计算流通量。
    5 _# ?* P- q3 ^9 h4 q
    " z+ C9 T1 K# n) ?! s8 W% `2 CMATLAB 程序设计: 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
    0 K+ g/ v* A' E4 D6 `- n) [
    . Y* u% s9 D% y) \- `# p6 w' C* p9 E! ?% k* R3 c
    7 x9 c, W. k% m" A- e
    利用以上的处理结果,可编写 MATLAB 参考程序如下:
    8 l) R* m/ w. I- H; T: y
    ; q- R  K' E0 u7 u" r; s4 o4 afunction ex20_3_23 g5 x, W9 P8 W8 b/ p! x+ q: V
    %*****************************
    0 [* ^6 u4 n5 ?% 扩散系统之浓度分布
    8 F9 }. b, W2 R, ~& B8 l3 I% S%*****************************0 L/ b* }* s- X6 F7 t; f
    clear( t1 }9 k  X0 i' c* l
    clc
    " B& O2 w9 M7 v% iglobal DAB k CA0( a" U) a, S1 N; E
    %******************************3 a& v! t  i' M1 a/ f- L
    % 给定数据& o5 _4 j5 g+ z) x5 v' Q
    %******************************
      z) J' L+ t! O: G0 q" Z6 G& V, O" h1 XCA0=0.01;
    " {! t. T7 c0 a/ n- z& \6 mL=0.1;
    / t. X9 h2 S' e; iDAB=2e-9;1 ^6 R1 k. D9 Q5 G4 n
    k=2e-7;
    / }* n0 q; l; |1 U$ M, y* rh=10*24*3600;
    5 J3 Y+ D& j5 w: }* f%*******************************7 D6 ~1 N& l  ?
    % 取点
    $ e; G" U0 i, f$ |%*******************************
    ( E/ g+ z  }) O/ Xt=linspace(0,h,100);
    ( I5 \; @- r: k2 {4 xz=linspace(0,L,10);
    1 n. a6 X5 r" @) i%*******************************
    : W1 ]& N/ I; l4 l0 O# z% case (a)8 d2 v! j$ F/ W5 Y/ x/ Z
    %*******************************
    3 q+ {/ k4 T# U- b3 |! S# um=0;
    ' i- O0 }0 |+ s- D; msol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);4 `' A' n! I& T2 N
    CA=sol(:,:,1);5 Y1 x7 {0 y& a6 y: J% q3 T
    for i=1:length(t)  ^! n( x: r% q5 [
    [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
    ) `0 q# {% r1 ~0 W3 S: A. x+ m NAz(i)=-dCAdz_i*DAB;
    7 t6 _& h/ G: T+ l( v$ bend2 o' }' [, b! E9 S3 U
    figure(1)7 x. g4 p" H+ @0 @1 L. g- U
    subplot(211)
    3 L% {9 @5 |2 c& R6 [% y: h9 msurf(z,t/(24*3600),CA)( R% m- ]5 ]! W9 o; D4 X; a
    title('case (a)') ' A  Y- o9 o4 Z" m/ |
    xlabel('length (m)')
    6 }  K8 v" ~8 O8 e8 |5 V0 y+ S! q" Rylabel('time (day)')) I1 S4 s& v4 Z. U
    zlabel('conc. (mol/m^3)')
    4 s  Y; w. i! y# U( Esubplot(212)/ L% `  \8 q  C& `6 m! U
    plot(t/(24*3600),NAz'*24*3600)( Y6 A, ^' P: \, u
    xlabel('time (day)')
    4 ^$ T4 v- a" |7 m  i) c  \ylabel('flux (mol/m^2.day)')
    6 P, Q7 m5 D3 c" A%************************************
    6 E7 ~( J3 F1 c% F4 s' @9 d% case (b)* K/ F- u1 ~3 t" O6 @
    %************************************
    " C/ B9 a) y- ?m=0;
    * a1 `! a9 ~' v! \0 F: hsol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
    8 z( J* S# C$ e+ aCA=sol(:,:,1);2 O3 B; z! z- F& c
    for i=1:length(t)
    . I0 l% m5 T- F* P# q9 P [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);" q% K  N4 i% V- L' g6 O0 k8 g7 H
    NAz(i)=-dCAdz_i*DAB;
    ! O4 Y1 x  _3 W$ C1 Jend- h  e6 M" q( Q
    %
    1 u6 ]1 y& D$ ^1 \- ]6 t# @: a6 cfigure(2)
    & g4 ~' R) c- E: K6 Rsubplot(211)
    , Y, s) a' V7 g8 M7 W$ zsurf(z,t/(24*3600),CA)( _1 r7 @- x3 }4 g& ]% \+ g5 l
    title('case (b)')3 y( Y, W; W, N$ d$ x# Z
    xlabel('length (m)')
    ' R3 n0 S( e2 Z$ G5 k% cylabel('time (day)')9 m, T# s% U  P5 Q% v5 {1 {' C+ y; x& z& q
    zlabel('conc. (mol/m^3)')
    2 t: C5 I$ a2 L4 Gsubplot(212)! H: _+ U! ?3 M
    plot(t/(24*3600),NAz'*24*3600)$ ^7 B( _1 m* ?' z! A1 t0 e" P( i
    xlabel('time (day)')
    2 S8 M, P3 `' w5 {5 l. |( Q* S. |ylabel('flux (mol/m^2.day)'), @" ?" B6 ?! }/ f- T6 x
    %********************************************& z6 V& q  L3 O% S; _% C
    % PDE 函数4 }! h4 U9 _& o8 [# n" w/ u$ [% i
    %********************************************
    , P% L" U: L" j: t4 d4 F5 W% case (a)
    ) l; L: p( J" q9 ~7 J%********************************************+ f* @) B; p9 W( z: @" e
    function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)! X) a* R8 I5 |( ]9 }
    global DAB k CA0
    ( ]* e$ J4 R" u( c0 ]( [( s: Ec=1;
    0 z& g* h0 b3 _- D. ]6 `8 N7 gf=DAB*dCAdz;7 p1 }3 H1 c& \& O# R9 \1 U
    s=0;3 e0 |2 k1 i% f5 o/ `$ E9 y
    %*********************************************5 F$ Q8 ~4 T# ^$ m  e: z
    % case (a)
    ( g' H  G9 H7 d5 Y- s. r: V%*********************************************7 A" c! J5 ]9 C- l7 ?1 L( X
    function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
    ' R& j7 {0 _, P4 Y: zglobal DAB k CA0
    ; L: W; S, t9 d# v0 yc=1;/ a7 V$ X- }0 O# ^* Y& E3 @% x
    f=DAB*dCAdz;: m9 a- B  ]" V' G/ Z8 H; J9 K
    s=k*CA;+ |, H8 v0 i& e8 ~
    %**********************************************$ p6 e# t0 j9 {1 w
    % 初始条件函数
    1 t5 q3 _4 U: W% N* m%**********************************************& \' `% H' |" Q' D2 S
    function CA_i=ex20_3_2ic(z)8 U3 t5 c& h& k
    CA_i=0;
    7 K9 T$ T* M1 n- t/ X%************************************************ 5 d) ~* A3 c$ _' z# [
    % 边界条件函数
      q" V; g6 O' O" m( f" E6 T  G%************************************************
    ! T8 o1 E3 k3 G6 S6 @function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)% k: m) F: I  {# E, z$ u  B. ~2 R2 R
    global DAB k CA0
    - l% E# O* ?$ x8 bpl=CAl-CA0;8 b2 H" r! a$ l( R* A" ~
    ql=0;& c# K. s& K8 a/ ?0 Q9 i5 {) c
    pr=0;
    % J: s9 u3 X2 c: F" G3 x- l) uqr=1/DAB;
    9 ~$ j/ L  g( ^, J4 M. o) E5 p5 c& a1 B
    ) e* c1 n) R: \  p$ r4 }————————————————
    ; M+ a5 E+ [+ \8 v版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。  f: W$ `# ^, m/ O
    原文链接:https://blog.csdn.net/qq_29831163/article/details/89711694% R( i5 e% d* w& D  ^
    " z+ t4 t" I2 M1 s& `9 S
    & P, Y  v# W6 }8 ~* G4 i
    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-7-30 16:48 , Processed in 0.331096 second(s), 50 queries .

    回顶部