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

浅夏110

; p( d  A: \$ w9 G% W题意解析：
) }5 q7 \$ t. @6 a0 y  l
) u* _  u8 Q1 p1 X+ M(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：


( Y8 |, _1 E3 _5 W6 f" R" ]8 n
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为


! R$ p9 c. V0 Z$ j9 [+ a4 [  T
& H4 d; a& V! B而起始及边界条件同上。

在获得浓度分布后，即可以 Fick’s law
3 Y0 h0 M2 s" r* P3 k3 h7 i- Q; H& u
+ y5 G+ G" W( U1 f3 @0 V3 b
! m9 H+ `8 b* k$ W  h. ?
计算流通量。
6 t; t  P( B+ |9 y, X  v4 Y  X: h
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
8 ^8 p6 o; L  r; r


利用以上的处理结果，可编写 MATLAB 参考程序如下：

function ex20_3_2
$ U& w) t$ c8 s2 h1 |. }%*****************************
% 扩散系统之浓度分布
1 ]/ o. P# P& X7 _' Z5 X$ z%*****************************
clear
# B5 J+ R/ B3 x; \4 u! w4 X3 Y2 H0 r0 bclc
global DAB k CA0
%******************************
% 给定数据
9 S, T  x! A$ o, t: p" H%******************************
CA0=0.01;
! O( k4 M# e% sL=0.1;
$ X5 c+ D6 s. dDAB=2e-9;
k=2e-7;
h=10*24*3600;
. V; N' S7 D+ m* H! U8 N( a$ X# K%*******************************
) Z5 W! r  k7 B- ^0 [: j* W3 j% 取点
. U+ A7 F4 f. j6 ]# k4 l1 L%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
%*******************************
. z) }. K( l# H5 Z% case (a)
%*******************************
, r  C1 M$ {$ n' `# Hm=0;
0 k! \. q; c, Y) R$ A8 Nsol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
8 C& ?7 a9 B7 \/ Yfigure(1)
1 V' [4 N! i/ E/ d; g0 ysubplot(211)
surf(z,t/(24*3600),CA)
% b8 Z$ j; x# ]title('case (a)')
! n' r, w. H9 ?$ p" Cxlabel('length (m)')
  O0 x1 {" z" _3 J- o: eylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
; J; F$ r3 O( Q  Q3 Cplot(t/(24*3600),NAz'*24*3600)
3 E" V( f7 G5 h" _  Bxlabel('time (day)')
* n4 z' h- x7 r3 Q- [3 a" _ylabel('flux (mol/m^2.day)')
%************************************
% case (b)
%************************************
m=0;
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
% x8 f3 }4 w( u: [1 d# ~$ {for i=1:length(t)
: J4 {0 V) L7 k1 F4 D# v [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
%
figure(2)
* h* X, x6 r% }5 e' C, Y1 W, `subplot(211)
surf(z,t/(24*3600),CA)
1 E: b8 Z" G3 U" g" k: t; `title('case (b)')
* ~: o* H* j  L5 R7 F4 Wxlabel('length (m)')
ylabel('time (day)')
zlabel('conc. (mol/m^3)')
5 ^- B& w9 K3 b+ Psubplot(212)
plot(t/(24*3600),NAz'*24*3600)
) I! O: h3 V$ |: t5 oxlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%********************************************
+ V; g" k' H0 v3 J1 E/ Y4 C% PDE 函数
6 e% S4 {, p- `0 ~; M# |%********************************************
' ^9 Y9 a" @4 S: B1 @% case (a)
. O7 p8 f. z+ |%********************************************
, h6 v. m$ v. B+ C0 D$ \1 }2 tfunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
c=1;
1 r- ^! o" h* j* w# nf=DAB*dCAdz;
s=0;
% o5 @9 L% n; L# V* J+ i$ v. E: ]%*********************************************
% case (a)
%*********************************************
: m7 T3 D0 \8 c# ?9 J4 r1 J) qfunction [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
5 S: F0 Q2 E6 O9 p5 G/ }% t( Eglobal DAB k CA0
c=1;
f=DAB*dCAdz;
% @# {6 `2 u2 z% Ks=k*CA;
5 r( V) \" L8 m2 O( }8 I: M%**********************************************
% 初始条件函数
%**********************************************
function CA_i=ex20_3_2ic(z)
CA_i=0;
%************************************************
9 j  p/ f& Z1 X' r4 w% 边界条件函数
+ Q: j2 W5 e' t. J& d%************************************************
$ g3 t& S2 a. l3 @function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
pl=CAl-CA0;
& C. h" ~' u) o/ N: t( K7 Qql=0;
pr=0;
" X. ~% h$ ?' e) e! B7 ]5 bqr=1/DAB;
/ G: z$ @" j5 u4 B8 V4 \  p* n
————————————————
2 Z& h! D3 R  \& h! F版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
9 V2 _6 z* `9 R. `, J) J* Z
5 ?6 _1 C1 h$ L' a* U