QQ登录

只需要一步,快速开始

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

    * o4 D6 @6 E' b' K8 s* b( `6 L9 w1 S% |题意解析:7 g4 [( \# W# E8 K' X- l' A4 A/ h2 A
    4 a8 z2 L8 L$ ?+ D9 y- R
    (a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:
    / s& X% n; g  h7 e# N7 A; f
    1 h& ~4 G# U, ~9 ^- V- O0 b# S
    2 [; F4 \8 p5 f) Y/ ?* D+ c3 S% U7 U4 ~
    (b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为6 r7 X  H, l4 v9 {: a; g1 ~& i
    4 w$ {' h* F- l! J

    : E4 f: A/ D- ?, ^  c/ W
    6 z: }6 ^* T+ F9 Y' [" B1 ]" q9 B而起始及边界条件同上。: s% G; ~, I" @  I& N! g
      N9 H/ U; \3 y: k3 L
    在获得浓度分布后,即可以 Fick’s law
    5 y* i* U, o. A1 \, \
    2 V* ]* C% ?, P+ ~7 t
    - g/ j- A0 M6 A( Y3 N1 M& v, R8 |" z. j
    计算流通量。
    0 c6 U" b5 s  k8 ?) L
    3 L( V  Z# a' c4 t$ ~6 MMATLAB 程序设计: 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
    - U2 J' v( d$ }
    : e6 f; |3 [6 x/ a* R  U, L' a' n3 r9 g! x: S/ \
    ; z1 F$ F9 i! h, R+ ]; z
    利用以上的处理结果,可编写 MATLAB 参考程序如下:: S: _8 R0 d! Y  E

    0 I% b: j2 e7 j' Efunction ex20_3_2, T5 G9 l: q7 @  a
    %*****************************, [. P7 B  F! R) E4 \9 r! A  `, ^3 g# m
    % 扩散系统之浓度分布
    5 O) D8 j% c% V/ O! X: k2 t& A%*****************************9 ~0 {1 E3 X, ?! k
    clear' N4 U2 P# I2 N8 q- l
    clc' Y4 x+ O; k0 @0 k* X& n3 P
    global DAB k CA01 g& O, ?" X  I. o- ]2 z
    %******************************4 z% F& ~: O" o  R9 w- L
    % 给定数据
    7 s: O! B: G( W& K: t* q3 f%******************************8 |5 }5 o* G( \5 A- C" {; O$ C( ^
    CA0=0.01;* h6 A; P" p* T) [$ U" _" Q
    L=0.1;
    * J# P, M! b! J* o( K2 {5 rDAB=2e-9;: Y- ]- y# p' L3 N% Z
    k=2e-7;9 S/ |3 ]& l& i' T0 Y9 O' y3 l2 C
    h=10*24*3600;
    4 G( Y0 I* A4 u5 i0 C3 n" {. u%*******************************! h6 U; g9 o% R
    % 取点% Z! `+ R6 `  w
    %*******************************1 n4 g4 K$ J8 N) q
    t=linspace(0,h,100);
    2 y1 p% m$ f, @+ U8 J" N( U* uz=linspace(0,L,10);1 I: t& M+ `; @  o9 J
    %*******************************0 Q& H3 o0 I/ {' }6 r" O' y0 _( d- s
    % case (a)& C; t" O0 n  `/ S; O$ t7 B
    %*******************************$ x$ J3 i5 C. J7 {, Q; R
    m=0;
    + W% D  d4 r5 ?- V7 hsol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
    ( u6 L% t: N0 J0 L9 ACA=sol(:,:,1);0 S; A3 X2 V# N! |7 j; I
    for i=1:length(t)4 w" g( C5 Q! I6 b' w  U3 B9 Y
    [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);0 q2 v" j( H/ b4 G& s$ v/ A
    NAz(i)=-dCAdz_i*DAB;
    1 E% H% C6 |: hend
    9 x% z/ W; t7 N# {0 \' v. d' D6 |6 Pfigure(1)
    # H% L: h. P2 C; j: @# ]  S3 rsubplot(211). {5 Y3 O8 n' k' Y6 r
    surf(z,t/(24*3600),CA)% Z4 `0 K$ y  J
    title('case (a)') , ?3 Q; ~  f7 h$ a2 e# L8 `
    xlabel('length (m)')
    / B3 ~. H& i6 w  h9 A; O2 {1 pylabel('time (day)')
    - |8 {9 Z$ `9 \; Jzlabel('conc. (mol/m^3)')
    8 O' w0 m- J$ i, s7 esubplot(212)
    : ]+ [) Y  s. e' Y% a6 \! _1 Tplot(t/(24*3600),NAz'*24*3600)
    7 j2 k+ h# ^7 O. X- ~# C6 vxlabel('time (day)')
    & X9 {1 n5 b4 U0 B. u6 |! o& Iylabel('flux (mol/m^2.day)')3 f" `8 C( d7 H3 T7 ~6 h
    %************************************( \$ e; ^( `' R/ O* _
    % case (b)
    , V+ x: A3 S. T/ Y( O5 U%************************************
    ! D2 V* Y: C% ]9 N. k* x4 Pm=0;9 f# y" L' W5 M. _3 q( @( ~
    sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
    9 ?0 }" u( M/ P$ z$ j* rCA=sol(:,:,1);: c! C- @0 O4 w/ c: v7 k
    for i=1:length(t)
    ( ~, A7 j& t- X6 L$ P6 j [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);) \# ?* D. {+ s1 j
    NAz(i)=-dCAdz_i*DAB;
    2 t! v5 p" q; Z  {. Nend
    % Z3 \" R6 Q! H) a) E% `%8 Q3 H" G1 u+ t
    figure(2), b7 Y9 C- y7 ~1 K8 e( B( d* w
    subplot(211)
    ; p- m6 |7 B( O2 ~' g7 Fsurf(z,t/(24*3600),CA)0 ]+ k0 l, ~6 m
    title('case (b)')) ?- ^% g9 E7 [* Y' X/ z% C
    xlabel('length (m)')
    + D  ]3 h1 C  x' t& s( f3 ^, pylabel('time (day)')
    5 x3 z( ^+ {6 Jzlabel('conc. (mol/m^3)')
    - A/ P. i* Y$ ~5 D' Fsubplot(212)
    - t* P1 e& D) }3 N0 s, Q! mplot(t/(24*3600),NAz'*24*3600)) [5 o  Z, I8 N6 g7 ?/ B, P
    xlabel('time (day)'), x  C( ~9 e8 n; p. {7 I
    ylabel('flux (mol/m^2.day)')% r5 K: J# x% y) y
    %********************************************
    1 q- h" q: J$ w% M4 ~+ L3 s% PDE 函数
    ! J) z( o) _. P% d6 B( e+ g' x%********************************************
    3 d+ |$ Z" T) e! W4 }% A( K8 |! G% case (a)
    0 q( Q+ g2 w+ j: X& _& f%********************************************" C; e) g! o6 D- z2 I, P
    function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
    * @/ o: g9 o% v' f, tglobal DAB k CA0
    $ f5 f3 U/ ], ]+ Q, A: x' F- Mc=1;8 e7 |% {! @: ^6 D- @8 R  Z4 r
    f=DAB*dCAdz;
    ( ]/ c5 D: e9 ~s=0;
    " b5 Z0 ~: W- M1 s% T%*********************************************+ w: ?: {0 s8 @4 j- m
    % case (a)4 g5 c+ a( ^* W" @; A# C6 e
    %*********************************************
    . R6 ?& ]8 q9 r1 b7 ?& Dfunction [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)% G' M! ~8 i# F# z, K
    global DAB k CA0
    9 u  _0 u* r. U& sc=1;5 |( v+ K! u) {  {
    f=DAB*dCAdz;
    , A: a  [6 g  s2 |$ o' p* N5 Es=k*CA;4 J4 y! a7 x1 G- {0 R1 x+ P7 W
    %**********************************************  h8 a5 z* H9 G5 v7 h, r6 @( X
    % 初始条件函数, t& z9 q# j/ z$ v! N% Y) H
    %**********************************************
    % J% h( N0 U) G* hfunction CA_i=ex20_3_2ic(z)# Y$ ~  F- A3 D& k4 g% f1 X" ]; s
    CA_i=0;" Z+ o/ f! S" T5 R' L7 }# F  q6 u
    %************************************************ $ @9 q6 f; C2 N- i: K
    % 边界条件函数
    % t  v* }; A3 a( Y, c3 h%************************************************
    ! r# `$ \$ Z0 C. o: F8 x8 n# ^function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)4 z' Y7 p+ u9 X
    global DAB k CA09 y1 M* S# _4 }5 n' d0 \3 w
    pl=CAl-CA0;
    ( Y9 @( U. V+ {0 V+ E5 L1 }  @ql=0;
    8 p6 G' P" j- Y6 {8 Ypr=0;
    5 }7 l& d/ K+ s# T1 O7 q+ I! |qr=1/DAB; , y6 y+ |0 y8 C. F1 e+ W
    # \/ d' I$ {8 I% C* [* K
    ————————————————
    # u1 b) W* ]/ X2 k0 W版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
    % s3 t. n" g9 L% i; b" M& E( H) a4 v% }原文链接:https://blog.csdn.net/qq_29831163/article/details/89711694
    ' h1 g& S7 @) I" H! `4 P9 ?$ y% u3 R4 t

    " j' N& b' c6 S  r6 G& }# q- C
    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 15:42 , Processed in 0.382790 second(s), 52 queries .

    回顶部