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

浅夏110

' o$ ?8 P$ B9 @  J题意解析：

6 t$ O: O  t# ?9 N(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
3 D6 M2 F6 I) ^/ M3 k" P* x6 b

# y. l: w2 K# n$ f
7 Z1 J4 c: O# k' `% ?(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为



, x% Z6 R2 Y' X0 l而起始及边界条件同上。
  i. O, p7 K; f( C: Z! Z  M8 [
在获得浓度分布后，即可以 Fick’s law
/ O3 X$ _+ c+ U; h3 ^
: x% K7 z' y6 u

2 ~' r: h  ]  G& [计算流通量。

MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
) S# e5 q4 A9 _
- }# {5 [& ]1 I$ Y( y. w5 S
0 O9 Z) p+ R2 t+ ?2 o
# i0 n" S4 ?- ?5 e' k利用以上的处理结果，可编写 MATLAB 参考程序如下：
  n% c6 u+ q$ U. u9 L1 M
function ex20_3_2
%*****************************
& `$ y1 [0 L, ^& ^1 F% 扩散系统之浓度分布
%*****************************
clear
1 d( w; g5 ]- y9 h, `6 Uclc
4 O* s# d! O) S; d6 X# x' Qglobal DAB k CA0
%******************************
$ S% C( r( S  o) J% 给定数据
% O' }) G- [& G4 |%******************************
( g4 v% i3 O1 V- a* R: b. vCA0=0.01;
L=0.1;
DAB=2e-9;
$ c9 u4 P! x& V- c* g6 Xk=2e-7;
h=10*24*3600;
%*******************************
% 取点
%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
4 d4 W2 Z, M; j6 o% B3 K# h%*******************************
8 `. I8 u. C  U) i, s$ G, [  x% case (a)
* ?. l0 x# Y! g3 C  s$ F%*******************************
! v( C% w/ `3 R; S' {" z8 ?& im=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
for i=1:length(t)
4 P# L# H% a1 P' W8 h8 U8 q# x [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
6 U7 K" J) z0 }2 m NAz(i)=-dCAdz_i*DAB;
' r/ i; _$ S4 d" iend
2 D( V8 s# k7 i* |% jfigure(1)
subplot(211)
3 Y9 P4 f& J0 _, v* ysurf(z,t/(24*3600),CA)
: q  G+ ]9 g* d0 I" h" l0 Atitle('case (a)')
xlabel('length (m)')
* y- Q+ H4 [1 ^: z" kylabel('time (day)')
+ |& \5 B1 K: Q9 b8 @$ J- gzlabel('conc. (mol/m^3)')
2 _# e9 h" b! o% usubplot(212)
& [2 N6 B3 f% g# n# w3 zplot(t/(24*3600),NAz'*24*3600)
5 K: m0 `$ r& yxlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%************************************
  F5 J8 q. W9 B2 e% case (b)
%************************************
m=0;
6 I' F/ ^' N# V$ o, W. z5 r% }2 b. K8 wsol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
1 P+ c( d- O* O" S# ~$ WCA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
1 |+ [. U7 ?+ W) send
4 T/ _* L4 G. [* x- Y%
figure(2)
subplot(211)
surf(z,t/(24*3600),CA)
; i$ U! H5 F: Q. s8 z, ftitle('case (b)')
4 M/ {7 Q4 k6 ~) ~0 R7 |& r* Qxlabel('length (m)')
ylabel('time (day)')
. H9 \- \! F$ b' \zlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
! W/ C$ c" Z* v5 ?6 }% o%********************************************
% PDE 函数
%********************************************
) L1 Q: m# Y# O. W4 ^% case (a)
%********************************************
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
+ t  V0 z' T7 k  Xglobal DAB k CA0
* K3 r# J6 Y6 a' t* V) z! h/ Ic=1;
1 \. j' v1 v1 ^4 Q" df=DAB*dCAdz;
s=0;
%*********************************************
; \) ?+ Y' {; u% m0 x9 P% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
( ?: V! P" |9 ?3 A1 R! Dglobal DAB k CA0
5 S- ^! g: b1 gc=1;
# _+ H& w# Y& ]: }' Wf=DAB*dCAdz;
s=k*CA;
%**********************************************
% 初始条件函数
%**********************************************
function CA_i=ex20_3_2ic(z)
$ }6 M6 {- G1 Y* m$ q" {( I( `CA_i=0;
- k+ e" L2 T/ W/ ]/ T2 [%************************************************
) `8 t% e0 p/ ~) K% 边界条件函数
%************************************************
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
. \, k. G0 k- Zpl=CAl-CA0;
1 Q. l5 U9 R1 @' s, F  _ql=0;
pr=0;
9 M( R* T3 R4 l" V; f0 Y4 Rqr=1/DAB;

6 M6 b# _5 f( E, |4 G4 y————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694

9 y, S. T/ M9 x
* ^7 A, m3 D" M6 O