QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 2552|回复: 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 L8 }0 w0 `) o& _% [) D题意解析:
    " s5 ^; g( a! _, h" o/ Z- h6 _0 W
    ! N4 M  W$ N7 a( Q% v(a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:
    # f4 n% T! _& s8 C+ m% k" ?* r, v1 o2 K$ h' f5 I9 y

    5 b0 V) Q( j# [6 c; w' O) E- Q  i8 X5 A4 S+ {3 H/ V/ E* p
    (b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为& a: r# `5 K: z  E) m
    7 A2 i4 Z+ s2 E8 g, r! }6 r
    * k5 B5 b( B; x- X9 V5 Y0 i

    8 t6 ~1 [( J1 H) V) V6 W: ?1 m/ L0 @而起始及边界条件同上。; P) }: j2 A* ?; C

    & s8 A; v+ m: ?在获得浓度分布后,即可以 Fick’s law
    . z: m" e& V! [4 i' ^4 t1 ~: I5 W) k- j3 t" Q- v' [( \0 d$ H

    / g4 ~* g/ ?: }3 f4 [- E' o
    " g' u+ o8 h8 v; B; b5 K+ S) Q计算流通量。$ _9 x# B9 g: I$ u

    , V/ @) I* t' u4 f8 B- [; ~" VMATLAB 程序设计: 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下# z4 k) r' N+ G' g. p4 u; B) J
    ' D$ @: d3 i4 L0 t! ?6 o
    ) q. p. ]" c: A/ s

    . k$ {: R8 T3 r+ y: D利用以上的处理结果,可编写 MATLAB 参考程序如下:* w. T; h) Y) N7 @. L

    + s5 f. x5 ^2 e+ _4 _3 o9 H5 x1 Gfunction ex20_3_2
    $ L4 h2 g* x' Q; t" T5 s%*****************************: l1 [6 {. ?. L) m! X2 g5 S4 N& ~
    % 扩散系统之浓度分布5 r2 y1 _2 n9 ]: j* }
    %*****************************
    ; `: J& e. C1 jclear1 P7 T3 E! B+ F, R8 b$ S4 Q+ Z8 O
    clc2 f- }% I: `6 f9 G+ f1 R8 ~
    global DAB k CA02 n. w9 u; e' }  z7 p
    %******************************* e1 V* L- r. s# F
    % 给定数据
    ( ^: d, R: Q: {% P# Z# S; D%******************************! D2 }2 r' s7 I
    CA0=0.01;
    . y/ S  |) p* ?9 g0 UL=0.1;! g9 C% x! o" l
    DAB=2e-9;
    ! }1 g0 t! _2 ~" {, h* q. o7 @( Yk=2e-7;
    # T7 T( x  V) }1 I, Xh=10*24*3600;
    5 {1 V1 n$ Q# A8 t' t8 a%*******************************" X1 p0 A. w0 @2 W
    % 取点8 ~) a, s' i4 K
    %*******************************9 F; V5 Y1 C  S, I
    t=linspace(0,h,100);
    % A( W9 R# |6 f7 Rz=linspace(0,L,10);
    3 q& A& U* n$ f- u# d%*******************************
    - L% w) N; a2 b! [% case (a)
    $ p& U& {0 C' `5 W, Z! w! o! r0 C& r%*******************************
    + |% A5 \7 D) v8 dm=0;  J( i  C9 L0 Y/ {7 D
    sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
    + {& f# A  u5 jCA=sol(:,:,1);
    1 w/ g  d3 v3 D' b3 p2 A) t; Sfor i=1:length(t)
    ; B; G, S7 Y0 e# N" C, Z0 q% R [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
    7 |' s) l! i% m# f, X# { NAz(i)=-dCAdz_i*DAB;
    0 W% F8 H: K# y! N; yend6 @* n, |8 x& r5 c' \# a5 U* c" C
    figure(1)
    7 \, R" I. a% w- zsubplot(211)
    # X" @& B' C5 ~& ^. O( Wsurf(z,t/(24*3600),CA)
    - {' {' j; y% [. I) y- jtitle('case (a)')
    + u, Q* b( {! L4 @2 }6 H0 A: ?xlabel('length (m)')
    $ x9 n/ ?  i- H; [ylabel('time (day)')
    * t- \- |. v" A; Ezlabel('conc. (mol/m^3)')0 O9 X8 |, J) t# f
    subplot(212)
    * _. U8 P' E! i) F: \! ?9 Mplot(t/(24*3600),NAz'*24*3600)# r: m& X" {* U
    xlabel('time (day)')
    . `4 K* R" `6 `, P5 Lylabel('flux (mol/m^2.day)'). r- `' Q0 Z: ^( A: C6 K, H
    %************************************
    2 n% k- I% }/ Y2 d8 r" s9 y% case (b)8 B. [) K0 X  x
    %************************************* A- P! @/ i1 ^3 E  X( n  S
    m=0;
    / f  S- r0 V) u0 Usol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
    / X9 }! k6 @9 i: Q; u  }CA=sol(:,:,1);
    # I' `5 y/ {: D3 Y7 o! Cfor i=1:length(t)
    5 a$ `7 t. g! v, s7 W6 f [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);7 \% a; r& v% `/ p  k
    NAz(i)=-dCAdz_i*DAB;
    0 [. O0 ^9 `  D9 c  U5 u) v+ Uend
    " J) ?& p- Z2 ]# o& c4 O%& @- a8 O) L& l% q
    figure(2)
    " Z1 W4 Z0 u$ |! m% Ssubplot(211)3 b' W8 R$ n7 C
    surf(z,t/(24*3600),CA)
    : E9 E" B; m; ~% d- Mtitle('case (b)')
    7 H0 @7 _- A2 @* e. I7 r! m1 Jxlabel('length (m)'), E4 [0 O3 b: T# A1 l
    ylabel('time (day)')
    0 ^& B& ^+ h' S3 s- `$ P& mzlabel('conc. (mol/m^3)')
    * p6 Y8 b6 @7 K- ^subplot(212)
    - y! v# E% o5 V3 d1 o( a+ v+ gplot(t/(24*3600),NAz'*24*3600)
    " F5 M0 ]) v3 q  Hxlabel('time (day)')
    % ]- H! i1 W7 M& @, q, wylabel('flux (mol/m^2.day)')( y1 Q1 ]0 v% b, Y( p
    %********************************************( Y$ M- O! F# ~" A1 ?9 I" o
    % PDE 函数
    ! F, d4 L- o' G4 A7 j  T# M8 ^%********************************************
    % Q) r5 }. C# G0 T# w5 A% case (a)7 X. j' |" Y/ g% R8 v  d
    %********************************************% `, _3 k9 ~- b/ ~
    function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
    / G( o1 c! i  f0 q" U* fglobal DAB k CA0+ w. A! w" o; h9 d* k2 v* q
    c=1;
    " q: x0 A& _# J( u# a2 e, jf=DAB*dCAdz;
    ; g0 g; v6 D4 O" \% cs=0;
    6 r2 H  @1 j3 N9 i$ m2 R% q# Q, z%*********************************************; ^! w' H% k4 N$ q6 S3 r
    % case (a)( ^5 L8 W, N, R! K
    %*********************************************8 D) o8 \( G2 x; c; E
    function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)% P: e# C8 b: f( F: [4 D- M
    global DAB k CA0
    " x- V  ^9 P& u/ Z% bc=1;
    8 H, R$ ^+ X; F( Vf=DAB*dCAdz;
    " ]: h- E) P9 F/ @. n* ~# {s=k*CA;+ L1 l9 k9 \; ^6 F" T4 B* N
    %**********************************************
    / Q& `8 j# t: C. k  q$ t$ P% 初始条件函数; L4 j; Z2 H6 a4 A; V* u* V
    %**********************************************: E7 J0 K! s( P$ A
    function CA_i=ex20_3_2ic(z)# o& ^2 u. x6 E9 p( Y3 _% ]3 B3 H
    CA_i=0;
    ) a3 P, l- f2 e6 \$ E%************************************************
    ; B; T6 h; Z! N0 \% 边界条件函数' I5 c9 R  m, X) B) n
    %************************************************# W$ ]6 y3 o% x. k$ r( k
    function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
    / C& p& ~, I! T8 }global DAB k CA0
    2 G- g  w" y4 v7 h( d  {pl=CAl-CA0;
    ( b8 z+ l  m( |9 ?6 s; O# R: l# Mql=0;
    / O  X, R; a4 p& [- `pr=0;8 N. B2 d! |" x4 T( J5 {( l/ W
    qr=1/DAB; 9 `2 c. V6 r- W) _  r6 @7 D( _

    ( A7 Y. X; V; }) d% W————————————————; A- y6 l1 o# G1 u8 r# X
    版权声明:本文为CSDN博主「wamg潇潇」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。' b3 \! O6 J! B" |+ z
    原文链接:https://blog.csdn.net/qq_29831163/article/details/89711694
    0 A1 `! c1 ?/ x5 Q# P
    0 q0 G+ x9 u) Z- V7 M8 {/ x* w! P2 z- }! I+ p, V
    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-28 23:05 , Processed in 0.499301 second(s), 50 queries .

    回顶部