偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

浅夏110

5 i5 P2 l; b/ Z, @# X2 i! h- |题意解析：

(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
; e3 K* }: U1 R5 s+ B* M5 w
) q& D, _; @' z3 |0 i
1 x; n+ S. v5 p7 F" U
2 w. ~5 P" @6 t& o9 S(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
: D* u1 W' H: J. n

  [+ T' U* N' D
而起始及边界条件同上。

/ I6 G, `# K9 Z* N  I5 y; U在获得浓度分布后，即可以 Fick’s law



, m1 h; H( T$ C7 \计算流通量。

& D9 `3 f# F0 mMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下

: Q# e' ]* V+ u
/ [! [) t* T2 l1 ?" V$ {0 |" f
利用以上的处理结果，可编写 MATLAB 参考程序如下：
9 `1 _; V- u% s" F9 F7 e* K
- f- m, y2 s/ pfunction ex20_3_2
%*****************************
% 扩散系统之浓度分布
%*****************************
clear
4 _$ d; Y: j& A* j, e! W+ _3 _clc
/ q# D% [. w7 y6 _- M9 d7 Bglobal DAB k CA0
%******************************
% 给定数据
%******************************
CA0=0.01;
L=0.1;
8 M7 n# _; r0 ?6 I; c0 z, QDAB=2e-9;
+ U" c$ k  ~; O8 P7 d# Dk=2e-7;
: y# L2 N5 Z# ]3 qh=10*24*3600;
, f9 ~& \2 j- F+ ~& M& B4 t, d- O%*******************************
- {4 J; c- ?- e) n# T  g+ t# `! B9 _% 取点
* F; Q! `% r' b%*******************************
t=linspace(0,h,100);
% C! j' C7 a1 V  d6 z2 Zz=linspace(0,L,10);
%*******************************
& Y4 n% k3 F; V% case (a)
%*******************************
m=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
4 C$ a% u8 C' Z- ?) K. SCA=sol(:,:,1);
- {) g3 v. ?) ?1 O# sfor i=1:length(t)
2 x" c: H/ i' P$ r, q3 Q- {8 L [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
3 Y. h5 [. @4 m- P# Bend
! I" _1 t- i. Z1 |figure(1)
subplot(211)
' I9 o; y' u2 O) B. osurf(z,t/(24*3600),CA)
title('case (a)')
, k) K( h7 o: _8 n- f" G, u3 pxlabel('length (m)')
! ^" ?  U5 N+ b9 p# eylabel('time (day)')
zlabel('conc. (mol/m^3)')
! E/ W( B9 M5 y+ o! B+ ]+ Wsubplot(212)
plot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%************************************
0 h; r; T$ x3 ?0 |% case (b)
, x# \& z1 F: ?" ]: R8 {' L%************************************
. n( q8 U& s9 W; R8 R! x- @# Km=0;
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
$ @5 ~$ c7 c: T4 e( OCA=sol(:,:,1);
for i=1:length(t)
8 q( e8 c6 E0 X+ t [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
%
figure(2)
0 G  \, x( f/ j& W8 n( Jsubplot(211)
$ ]" l! N2 j9 ?/ Xsurf(z,t/(24*3600),CA)
title('case (b)')
, E0 G5 o* T6 V# G" q6 V9 |xlabel('length (m)')
ylabel('time (day)')
" K, A* \, t( e& n" W; qzlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
8 _) G) r' ]7 Qxlabel('time (day)')
5 K( M! Z, r. Z; N/ e1 Qylabel('flux (mol/m^2.day)')
%********************************************
' V+ r' f! B" o$ B% PDE 函数
+ B, a+ X, ]' o& U$ a: i%********************************************
% case (a)
%********************************************
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
8 g! }) \1 \% h. |) t8 X3 y4 H4 Uglobal DAB k CA0
9 t# H5 C% Y% E' F# j, sc=1;
f=DAB*dCAdz;
# W( o0 t: K$ Js=0;
+ m0 c- g( `6 `; e4 J$ G# w% Q( c%*********************************************
: L  F/ [: h0 F% I# W: K. c( n/ L# w% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
2 o  w4 k5 p3 B* ^5 }9 Lglobal DAB k CA0
0 T! V) s# A# U. p; u' Vc=1;
f=DAB*dCAdz;
/ R4 q0 d  F) A4 L4 fs=k*CA;
%**********************************************
8 @6 f$ o9 ^3 J3 x. X% 初始条件函数
2 ~  h9 k9 d! P, [%**********************************************
function CA_i=ex20_3_2ic(z)
4 m) @; T5 L# C/ C( v2 b3 b; {5 F' iCA_i=0;
%************************************************
% 边界条件函数
%************************************************
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
2 o2 i6 E/ X0 W  m) T% gpl=CAl-CA0;
ql=0;
/ W8 U7 Z5 e# {, ~6 d/ ppr=0;
qr=1/DAB;
( y0 v  J! t3 V8 ]2 k3 R4 {; t
————————————————
9 M( z- L" v2 u$ k4 L  A2 m* ?版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
% V( P/ D! O! U- _4 i原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
" u7 l8 l" s, N
0 n% q6 }8 r* Y# r: s
* H, Z" R1 J8 y( r0 V' {9 L: Z* j1 N