QQ登录

只需要一步,快速开始

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

    5 ?1 f- S0 D0 x) l* k# n/ K) s3 X题意解析:! p) l, a# F& N( H  g( ^
    " T% T6 {2 F. `0 \& L1 U7 b. r
    (a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:! G, G. S: F- m/ r( f+ C  N

    8 o" X( j6 U: U6 r: ?4 `8 ^* O- ^0 w
    8 K! G# |! {: r: `) I0 a: o. G5 s
    (b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为
    / |9 M) r" P) w$ R. x$ A7 L( s4 E/ V: P  s

    4 j6 ^) C7 T" R8 Z# p' A7 r6 r: [$ N3 f  O: {3 O+ T
    而起始及边界条件同上。/ X  H0 ~" [5 N1 u# U5 A3 ]! y3 a0 j
    1 e  W7 |& a6 t& c9 }% g/ e
    在获得浓度分布后,即可以 Fick’s law
    5 c9 Q0 g, a5 n* j; [
    ; |. z0 Z$ ?: h$ ]0 j
    . J2 I# u8 d, x" U3 g8 i5 Y6 ]  I0 m0 n! O" H* @9 {" B
    计算流通量。
    6 I0 ^* C% m( v/ g- D, v7 O& G9 ?: o0 A8 ]4 O( A9 ~
    MATLAB 程序设计: 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
    6 A5 V! y; N, m5 W* \3 a* w0 K0 Q2 n! \* e/ O( J5 @) e! a

    ; h: l5 G9 v: r8 K
    " A7 S5 {: G' ^1 W, M# _+ i" f0 ], ^利用以上的处理结果,可编写 MATLAB 参考程序如下:
    % T. B: V: q: R" X+ `' N, h3 z1 G) ~" X: O% H' i# E
    function ex20_3_28 W' U1 g4 w, p% n8 {* G
    %*****************************: r4 q0 h8 `9 j
    % 扩散系统之浓度分布* N1 F3 O& a0 a8 B# ?
    %*****************************9 Q. w3 O8 j& [# K6 E) M
    clear; {6 e+ i" @; q7 @+ C5 T' f* V
    clc
    8 n/ k8 \; \. Hglobal DAB k CA0
    2 G( a) Y, O) J5 O& x%******************************2 l3 F: @1 d" J" ]# W4 c% S
    % 给定数据
    ; D3 |& G; W) O; d  b0 K%******************************
    8 O) K" K0 B/ l0 RCA0=0.01;
    1 v/ `9 H( n9 Q" m4 S/ h: qL=0.1;
    / ~) J8 A! x* ^# LDAB=2e-9;
    $ ~3 u- q5 k$ z& y2 Y" i6 W9 Nk=2e-7;
    2 r5 V. @1 E+ `2 R: K# g; H" e4 |h=10*24*3600;1 @  r1 e, h! G" d+ e# O
    %*******************************( J5 j8 w0 ~7 L& C
    % 取点
    $ |/ A6 m$ ?% y0 Y6 r) B%*******************************0 Q9 z* Q0 {' `
    t=linspace(0,h,100);
    $ ~& K: y) @* v# r: iz=linspace(0,L,10);
    + p7 V! b7 z1 x5 d3 @1 Z' I! K%*******************************0 m( T/ a8 n( y% y- b' [2 M# {
    % case (a)
    / ]0 f  G2 [+ u; `%*******************************4 w9 v9 \* ?: }" m
    m=0;( ?& n1 T; Q8 C2 L% `  V0 T
    sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
    . ~5 E. c' K( K' ACA=sol(:,:,1);
    + g1 t0 b, q, [; n8 Q9 N) lfor i=1:length(t)3 i- H6 u& E8 G3 W, `  O1 m# q
    [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);" ~/ L: U) \5 M  M/ B0 i* k  U
    NAz(i)=-dCAdz_i*DAB;* S% c3 d* M/ H- V6 j( j
    end
      X0 {4 M" I4 Bfigure(1)
    , [; C& Z3 k7 Jsubplot(211); V; F0 E6 P# I3 n# y. u! |5 R% t
    surf(z,t/(24*3600),CA)6 y0 }! h) c) i1 ~1 l* M8 {8 B
    title('case (a)')
    . X9 J  ?/ I+ t) s. zxlabel('length (m)')4 x  q. g, X' B, \& c! ?% F% D
    ylabel('time (day)')$ ~; s& X" z& f* I" C
    zlabel('conc. (mol/m^3)')
      T. q4 }: P8 }6 I6 ?8 }, esubplot(212): I0 T: C) d9 `& ~4 b" p! ^# b
    plot(t/(24*3600),NAz'*24*3600)& [% o; R( y2 n$ ?& ?' v: [3 R
    xlabel('time (day)')6 Y4 k8 W- ]& l- ~
    ylabel('flux (mol/m^2.day)')' f, p/ [5 Z8 A& u2 |. i
    %************************************
    $ E( ^, Q2 B' J/ N, @% case (b)
    , l+ [6 r' g' Z. u/ e%************************************  f3 F8 \2 i0 F: g
    m=0;% ^- {! C& t3 T, b8 }
    sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
    $ l1 y: ~- j! d/ \1 M/ \9 OCA=sol(:,:,1);" L  m6 U, g1 \
    for i=1:length(t)$ M8 N9 D6 k5 }, B# e
    [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
    # s2 J7 c" q# b, n NAz(i)=-dCAdz_i*DAB;
    7 ?3 B9 {6 S& Tend* [) M( k  t- f( I3 ^. }2 b! L
    %
    % |+ X' ]. S; i" w/ {figure(2)
    3 a# y: _# w7 J9 p3 usubplot(211)
    + t+ r; C+ a8 [surf(z,t/(24*3600),CA)
    4 p% T% j3 ~0 w, otitle('case (b)')
    " \# M7 o8 X9 `5 _xlabel('length (m)')# y/ V3 @; K4 p" i4 D4 K8 i
    ylabel('time (day)')
    1 s* q& |2 D3 W) xzlabel('conc. (mol/m^3)')0 P0 h0 c4 z, G9 t& d( w0 C5 V
    subplot(212), b. F5 v: k6 r( u6 J
    plot(t/(24*3600),NAz'*24*3600)
    - F  P; \/ Q5 K/ @- h2 |7 Ixlabel('time (day)')
    8 A  P/ H: Y! b3 fylabel('flux (mol/m^2.day)')* C" ?( j- {; t1 [' W
    %********************************************
    % b$ d$ K# S& Z: o& s1 z% PDE 函数; f/ i  }. b4 y: N
    %********************************************
    ! i+ o. f. Y& s& ?' V* c0 i% case (a)% n, F2 t5 [$ F$ k8 |
    %********************************************% b+ Z; n* e7 ~, W9 X/ P" [# n
    function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
    6 a7 M' U$ ]2 O9 s" bglobal DAB k CA05 R- O, N  V% M$ Z, Q* I* l) b
    c=1;
    / |% E/ z) Q0 e) e0 r7 Nf=DAB*dCAdz;. X- S6 U$ J2 P0 k$ d+ H$ E! D! c
    s=0;. _9 C7 u5 w" S3 b1 x! a% w
    %*********************************************
    " E5 [( L3 v# c& t/ F% case (a)
    " U/ ]0 A. a; Y7 Q8 h%*********************************************
    " A, _6 m& F, n$ Pfunction [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
    2 Q! e1 w, }$ J. w6 qglobal DAB k CA0
    ( v: k- @& L7 S- j) h  C# ~c=1;
    ! i4 M1 @, k7 I) h7 gf=DAB*dCAdz;
    ; k) }$ p+ `/ W( J. {9 Z+ ps=k*CA;! ]+ C1 M3 u. J6 Y5 \
    %**********************************************  Y' i$ L# Z; I% }, w
    % 初始条件函数
    7 b9 f8 O0 G3 g/ x+ w2 a%**********************************************
    . L- F% ]7 ^& I! Kfunction CA_i=ex20_3_2ic(z)) V4 Z* E( a8 o( ~
    CA_i=0;
    . L- P0 t7 f$ W( j%************************************************ 6 K8 K- J( L) A0 c, K  b
    % 边界条件函数
    2 v8 U- Q7 N0 O8 t, \1 H0 I- T%************************************************
    : X8 f# v, {- x% ^8 H- @& Cfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)6 J" ]1 ?4 G( U8 u
    global DAB k CA0
    / K& T$ g; A1 a4 vpl=CAl-CA0;& M1 H; t, L- x& x/ _5 g' K
    ql=0;! i4 g6 ~/ Y" c: C9 t" v
    pr=0;. A  @3 c( a# E. f) w
    qr=1/DAB;
    ; I( u9 m* g6 |' T; Y1 t1 K  O! B% c9 u, n
    ————————————————
    ; B9 j) R0 t2 T; n. A9 Z  s2 T: ?版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。) x3 v, o# ]2 Q" Z
    原文链接:https://blog.csdn.net/qq_29831163/article/details/89711694( h# K9 V" m2 o2 Q  M2 |. I
    ' ^' ^: S8 n3 N8 x4 _: g
    " g0 P, q0 x3 ]& D/ m5 g0 \
    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-4-12 17:36 , Processed in 0.394464 second(s), 51 queries .

    回顶部