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

浅夏110

题意解析：

' ?7 k2 H: ?/ N9 \4 K% O(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
' h- O5 U; s7 y3 \

9 K3 q4 a5 }2 s; A
2 E6 ^# `+ b! ^( U' t: O7 b4 S8 s(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
6 @' v3 L9 e; \% @


而起始及边界条件同上。

在获得浓度分布后，即可以 Fick’s law
3 n& K( @5 f# P$ d% J2 T5 B


: I& h& I# `- e& n; e计算流通量。

6 b3 c$ s5 L# V, QMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下

/ _# r3 }2 x; c. |3 `" S

利用以上的处理结果，可编写 MATLAB 参考程序如下：
2 q+ W. G" W7 L, I4 q# s
. j# i1 Q. i% v7 _- y* f* E2 u( lfunction ex20_3_2
, x5 @) r+ @, `. N9 F$ A. a' u%*****************************
% 扩散系统之浓度分布
%*****************************
2 Q* _! O0 n$ P0 kclear
clc
0 a$ m" q5 o: g9 N% c* L1 Wglobal DAB k CA0
. u, \6 U1 o# c  v+ l%******************************
9 M, w6 Y3 Y6 v# W! H' x2 y% 给定数据
%******************************
CA0=0.01;
L=0.1;
DAB=2e-9;
k=2e-7;
h=10*24*3600;
, J& @& v0 `$ L7 p- T! d%*******************************
% 取点
%*******************************
t=linspace(0,h,100);
: \" N0 z  K* y! M" Gz=linspace(0,L,10);
%*******************************
% case (a)
2 C# _: l  y/ ?- A! D/ y; `%*******************************
m=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
) L8 l% M" R7 J8 s1 s: lCA=sol(:,:,1);
! E0 B0 Z9 ^. b6 ufor i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
- \4 u5 H" O  @' {. g) W NAz(i)=-dCAdz_i*DAB;
5 F) _. o4 k0 T4 pend
figure(1)
: m. |; l: K% P2 gsubplot(211)
surf(z,t/(24*3600),CA)
0 s+ h1 h6 G+ V" I3 }3 W% otitle('case (a)')
xlabel('length (m)')
ylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
8 T4 z+ a+ m' c! p+ p5 xplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
2 t9 `1 v& K, o- d8 rylabel('flux (mol/m^2.day)')
) `% v# X8 e8 b3 F8 j%************************************
% case (b)
( Z2 j; X  h6 `6 F; N9 ]%************************************
m=0;
# r" \- R7 l/ U# w: C8 Y3 Bsol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
+ C" ~5 F/ O: Q1 @+ ~. H; _) [9 t5 s' Wfor i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
%
& j' y7 `( T3 c( Z; T7 w1 W- G% ffigure(2)
subplot(211)
surf(z,t/(24*3600),CA)
title('case (b)')
xlabel('length (m)')
2 Z: N8 ~7 c# w$ z5 wylabel('time (day)')
6 c) G. M: d% e3 ^8 h" Pzlabel('conc. (mol/m^3)')
3 w) K/ J! a3 r! {+ U( P' m* Msubplot(212)
* I* Q+ l. D2 R3 {6 f' E; Iplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%********************************************
% PDE 函数
' e1 y+ M  Q+ @+ Q) l8 s%********************************************
% case (a)
  E! }9 u$ I% ?%********************************************
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
6 D9 X+ `. s& s, e! Z4 pglobal DAB k CA0
c=1;
f=DAB*dCAdz;
s=0;
4 e: L( t! \& |  u) q%*********************************************
% case (a)
1 S7 c: m8 s3 ~9 j/ G%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
0 c* {2 U2 s( Q! z, i% J9 W# M; Gglobal DAB k CA0
c=1;
f=DAB*dCAdz;
s=k*CA;
%**********************************************
% 初始条件函数
%**********************************************
( L& A9 }! r7 A+ ]  e( Cfunction CA_i=ex20_3_2ic(z)
CA_i=0;
%************************************************
% 边界条件函数
0 O, H+ z' w5 c( u%************************************************
! q& v) G6 @" V1 _* V" tfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
) |0 F& ]! h# `8 S. b' `! dglobal DAB k CA0
pl=CAl-CA0;
ql=0;
2 H' T  e9 x: ?- V- P3 }0 ^6 Jpr=0;
qr=1/DAB;
* H) p* O6 u- w  E  P
————————————————
) r( z8 O# o. L! D6 |版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
6 e8 r  W, T! {1 J0 P9 J" _原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
) L; p( ]; V, a; q4 v$ Z) g
1 V0 P( k  y! |( x4 y