QQ登录

只需要一步,快速开始

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

    $ z5 l5 V7 W1 B题意解析:
    , }7 |3 [3 Z5 j: H4 ]0 E7 V& H! P& M* r: m5 }
    (a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:
    & b9 K2 e+ u' W0 O  s+ s" M' ^( A  O% j& x1 J8 ?
    . g$ @5 s! l1 b7 C9 c5 Z
    ! I4 \0 x: m; W9 p4 H  n5 V) [+ N
    (b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为
    ) X  S; E& i' l. _! N9 y# E
    * |# N* `9 p4 y9 ~/ }! {0 l6 L, W4 c* _( K( I. \. g
    3 }1 h5 f& S* K# ?$ ]  ]$ |
    而起始及边界条件同上。, W/ ?  e. e$ Q7 v( Y

    8 ~8 r. c# s9 g* C, J# |在获得浓度分布后,即可以 Fick’s law3 U2 `3 ?8 \" J+ I/ \- J) L5 `
    ) n0 t" x7 W! \5 I3 t5 p
    + W3 X8 L# j& ^- j0 V2 I  y. D

    + }5 y/ `1 S5 M/ K1 [) i6 Q计算流通量。
    : I, i5 b) K' F+ d' ^0 ?5 c
    - k" z( Q, K/ PMATLAB 程序设计: 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
    + E* D. r- Q5 o& C) k' w; ]6 c8 g4 [/ n

    1 N1 D( v  B* y% Y% @
    5 r9 ]$ G' ~: l4 e6 w; |) m# R4 F) j3 {利用以上的处理结果,可编写 MATLAB 参考程序如下:
    : a& \) U  a% b" P
    4 A/ n. n  [9 f3 k4 A5 Wfunction ex20_3_20 r7 X# x/ m" N0 f* n$ ^- p7 c4 f
    %*****************************1 {# M1 Y, J" j% S$ f) R/ C: w# e" c
    % 扩散系统之浓度分布
    ' q0 r/ N9 L2 f( H/ Y" o4 G7 i%*****************************3 P: E# e3 N! G6 \1 u" {0 F, X0 d) ^
    clear4 S6 X  L/ e4 a2 U5 K8 f
    clc% ?7 L4 m0 q- a$ G# Z6 \
    global DAB k CA0
    , J( g6 q: ?3 F- p%******************************+ Z2 `& {* |( n$ g- Y: E& f
    % 给定数据
    3 b9 {/ C' L8 v; e%******************************' K$ n7 l) o" w' Q* o
    CA0=0.01;
    + }( q3 y8 S4 I5 n, w1 }L=0.1;& I4 J+ _! L  D: D6 }" o  R# D
    DAB=2e-9;6 T! n8 }8 d4 j* A
    k=2e-7;( ?' ~4 [. Y7 g, @1 C! R; u% Q) l
    h=10*24*3600;' h( L* [$ `& m$ j" g2 `
    %*******************************
    0 Z+ I3 a8 j: J3 r% 取点3 k4 g9 F  I% u% {6 s# X$ N' d
    %*******************************- V8 u' S3 X0 v+ U
    t=linspace(0,h,100);
    ! M$ w& z, y, W, r) |z=linspace(0,L,10);
    3 {' y5 c9 B" p: A* _& L  [* q%*******************************' h9 R" N2 R- d8 z/ g
    % case (a)
    ; P1 _& Z7 R& v6 ?- P  s/ s( e2 f%*******************************. n; j' I' v3 H
    m=0;
    - H0 K; a" h* g; g: v( dsol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);! N; U* y6 ^6 p
    CA=sol(:,:,1);0 ?6 {! X3 X% c6 T0 v
    for i=1:length(t); B* {) i7 _! a1 R) |0 G
    [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);/ b1 n# j4 X- E: G
    NAz(i)=-dCAdz_i*DAB;
    2 Y5 M+ ^  @- H! P8 Yend( [9 F* K6 Y  a7 S( u
    figure(1)3 l& r4 q2 M# a  @; H; [" }
    subplot(211)
    & I3 W4 {( S) }  K: [/ Tsurf(z,t/(24*3600),CA)
    ! o! J5 K0 ?" n5 a0 G- z# b/ Ttitle('case (a)') # E1 p  D/ V0 p' E. h
    xlabel('length (m)'), m' c2 t3 U5 {( `
    ylabel('time (day)')
    ' N2 n9 X4 h6 Lzlabel('conc. (mol/m^3)')$ H8 ?* F9 \. G8 a
    subplot(212)0 Y6 L" k8 h$ l  L: J$ v0 [6 n+ |0 p
    plot(t/(24*3600),NAz'*24*3600)
    4 S5 ~; E! |* y* I! uxlabel('time (day)')
    ' ?  q$ t8 U+ j* B6 Iylabel('flux (mol/m^2.day)')
    ; \' Q8 ], g2 x' B) w) P%************************************* a1 j# Q4 S% U3 e2 ^4 E9 s2 r# m; [
    % case (b)
    . p! z, L, k% n' X5 j& g& ~%************************************5 K% V, x) v; E7 s: U
    m=0;) F! f& f/ q+ Q: g5 T( A
    sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);5 A# g7 J2 r! M! b* Z8 }! _- z, M0 h3 t
    CA=sol(:,:,1);
    0 G& u, Q3 I7 W* Mfor i=1:length(t)
    # B6 o7 y+ C, ^0 B; H2 T6 c, T& p [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
    . x0 N; \3 t* n, @! ^  i9 y/ B NAz(i)=-dCAdz_i*DAB;# k6 @; }! A3 S9 A  x0 E
    end
    ' m- k8 {: p. E* o! m; E- N5 I%$ k& t4 K8 B3 R  t0 p/ `6 `
    figure(2)
    ( x5 O: T" d0 B% `/ @( wsubplot(211)
    1 T4 k9 S" R# G2 osurf(z,t/(24*3600),CA)
      y+ h$ j/ x' G# r5 gtitle('case (b)')7 t: V; r- O, j8 z
    xlabel('length (m)')
    4 G5 E( D( {) e$ X0 eylabel('time (day)')
    / W% |5 L, ]4 o' s$ D2 E$ Lzlabel('conc. (mol/m^3)')
    ! G$ g% U/ E$ S& A7 Lsubplot(212)
    + i  K$ C" {# K7 a( Aplot(t/(24*3600),NAz'*24*3600)0 E) h. m5 \! `0 i" S
    xlabel('time (day)')
    + [: q" ~7 Z, i' `- _ylabel('flux (mol/m^2.day)')7 P1 m7 ^5 f4 r$ P! ]2 n% H7 o; q7 D7 m
    %********************************************
    4 a2 h* s, B; }% C0 c1 s1 H, b% PDE 函数; r# a8 W7 V/ \
    %********************************************
    / D2 E' F7 }' q# d& U! G( Q% case (a)+ P* O9 S" L8 c* V9 z+ F
    %********************************************( L( s  f( I5 w: z0 h( X
    function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
    + M5 g' W1 `4 l& W2 lglobal DAB k CA0
    - G  r9 V) L) O/ `" s: oc=1;
    ; d: y" r/ D3 W* L% R; _: h: P) e8 Nf=DAB*dCAdz;
    % w& Y! B3 b% T6 j- u) c6 rs=0;0 t' ?+ \- p6 ]. ]. v5 O
    %*********************************************
    / f/ z& I( Z" n4 }9 x2 t& I% case (a)
    , m: H0 K5 ^; a9 d! a3 j%*********************************************' F# T0 O. I2 p/ Z( U- Q7 E5 H3 N
    function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)( i% J: j2 M3 G" g6 E. b
    global DAB k CA0
    4 Q3 t6 b7 ~" h) U4 Zc=1;/ A; j/ e  x" u7 r( H# i
    f=DAB*dCAdz;
    , o  j* r1 [' E! b5 Ps=k*CA;9 m. F  e2 n* _
    %**********************************************3 Q% w( Z4 v$ \- I
    % 初始条件函数! F+ C8 \" D  v
    %**********************************************
    5 U& D& J/ `5 t& G% ?& Nfunction CA_i=ex20_3_2ic(z)2 G: N' c' _; q! I/ K  \
    CA_i=0;
    ) M. u" B  S$ @! d7 q7 p%************************************************ ! d0 v) j' e: H: B/ s9 O
    % 边界条件函数" v/ s) Y1 @* u$ q7 s' d
    %************************************************
    ( `9 r, g* p- e& `7 }/ Zfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
    / H3 J% X" B; O2 x  M3 O0 tglobal DAB k CA0
    % E0 ~" Y) y% Epl=CAl-CA0;
      q% p# v- O7 ~ql=0;
    5 `5 ?  \, c9 l& R8 ?pr=0;
    * Z. V) ^: E" W( _) X  cqr=1/DAB;
    0 j+ k# l( V9 b* t- M. w3 m7 {6 K/ Q8 t" V  u+ i( l: r, W
    ————————————————: c$ u; P9 S( J" g' M; Q& p* \3 w
    版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。. N. _+ ~8 `7 W
    原文链接:https://blog.csdn.net/qq_29831163/article/details/89711694# g( D% \6 @: X8 b( J; m6 b
    2 M. X) p* S% G* q3 G! ]! U

    7 p% v) W6 t8 Y6 u, a
    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:38 , Processed in 0.369863 second(s), 50 queries .

    回顶部