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

浅夏110

; K; P. ~. J* I! t2 K题意解析：
- e, U% b% L% l6 u" H  ~
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
4 ?" s- k4 |- P' h" `( y2 H1 r2 p
' Q/ \7 B# I" q1 j& s6 Q0 U% e

2 o; z; n* J( P& h& D) c1 _(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
3 O: D- A0 f6 `) e2 |, ~

' T- n: F1 H! J$ V
! A8 J7 D8 X5 ]$ U而起始及边界条件同上。

, u, X/ T! s2 `) U在获得浓度分布后，即可以 Fick’s law
, V' M" a6 V/ j# Q1 e+ X2 v
- o: O0 D( N- ~) F
% f* b2 K( c; z0 T/ A- R4 D
7 \% F* H( ?, o计算流通量。
: A, {" R2 t5 Q5 z8 h# @6 N. {% a6 y
0 a6 _! T5 c: T6 {5 YMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
5 L' W7 [' _7 B4 m8 v
) Q. N2 l$ Q" n/ k

  N2 n) x2 Z, E) @& o! F利用以上的处理结果，可编写 MATLAB 参考程序如下：

function ex20_3_2
8 i% X6 g* i  i! W3 B9 C) j7 v%*****************************
3 A# }+ A! A: }3 x2 g% 扩散系统之浓度分布
%*****************************
* P2 |! f# a, @5 }: yclear
clc
global DAB k CA0
* z0 B! t4 H/ |# C! k- Y) f%******************************
% 给定数据
2 L& ]0 ~0 R5 B9 r. c%******************************
$ N4 Z+ S7 N" X* g6 v7 R: YCA0=0.01;
$ [7 x5 ]6 ^. a  l( j. z% g) AL=0.1;
DAB=2e-9;
( E2 e. U6 {$ j5 T$ qk=2e-7;
h=10*24*3600;
& e9 A& I+ D4 l. ~' a; a3 Y) t+ k%*******************************
: }% e9 }5 D* H! L4 m$ M% 取点
) f: d& \+ E  n) L. t2 v%*******************************
# z8 V: ?$ Q2 t$ P, V2 ~t=linspace(0,h,100);
z=linspace(0,L,10);
%*******************************
% case (a)
: N* ]( ^8 Q5 ]9 Z' B0 E  z%*******************************
8 D) l  m2 [; S5 Nm=0;
1 V* q" r/ X+ u4 y+ B* V& E' J! ~) ?sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
/ c" M# I+ }" o2 N1 }+ ~: P4 yfor i=1:length(t)
! }* a0 H8 ?$ a [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
figure(1)
0 h" u4 ~" A# L& j" V- Qsubplot(211)
surf(z,t/(24*3600),CA)
title('case (a)')
xlabel('length (m)')
( V. {, [1 [, o6 g2 Fylabel('time (day)')
zlabel('conc. (mol/m^3)')
( O% {1 h4 I' y( ^" b5 Ysubplot(212)
4 v! b: G6 X5 K! h/ a) I% t. pplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
; S- S& Q! [2 D* Qylabel('flux (mol/m^2.day)')
%************************************
  L" y# {7 M, x0 R4 X% case (b)
$ G, `2 |) l. [2 J/ I%************************************
; b: s7 A0 M3 K$ Q1 Zm=0;
# {& z+ L' p& H) L+ V5 ~( Z2 h# d, I$ \sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
% f7 _, i, s7 d, [CA=sol(:,:,1);
for i=1:length(t)
1 }0 U- L  L# ~1 u2 Y( g& L+ \5 P8 _ [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
, S" C3 d6 J4 G+ D- s& K0 O0 G  i NAz(i)=-dCAdz_i*DAB;
end
%
& M  O4 o( g! R# r% I+ sfigure(2)
subplot(211)
surf(z,t/(24*3600),CA)
' y- I4 f  H# E1 ?/ o, utitle('case (b)')
5 R7 B: e5 D( H4 f' s( Ixlabel('length (m)')
* W9 F, l' l3 v" R- W% P/ O  K: ?ylabel('time (day)')
9 b- U' O! m) g( Z2 d/ O! u+ O6 {  N6 ]2 uzlabel('conc. (mol/m^3)')
# q( @$ v! |8 B) E! z% Ssubplot(212)
plot(t/(24*3600),NAz'*24*3600)
* z4 W( _3 |: f- D% m& G! mxlabel('time (day)')
$ ]9 P+ A8 }- jylabel('flux (mol/m^2.day)')
" x- d5 T" P7 n; G0 `$ E; z%********************************************
2 i) ~( M1 Q/ P; t+ j% PDE 函数
%********************************************
* v, \5 c8 `. Z6 j# [( Z2 d% case (a)
%********************************************
7 l9 E" H* o* y4 |! ?  }2 _/ J5 z2 Ofunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
1 H  ~" K" ]- C& v8 V5 j2 b" P: y, ?c=1;
f=DAB*dCAdz;
s=0;
%*********************************************
% case (a)
%*********************************************
& y+ m; r% H5 ~. W6 q0 ~/ ^function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
! K8 g2 F( n' h" E% qglobal DAB k CA0
c=1;
* u" l0 D( n( o8 P+ Y0 \/ Df=DAB*dCAdz;
s=k*CA;
%**********************************************
1 z1 v- U$ U; T7 r- L% 初始条件函数
%**********************************************
function CA_i=ex20_3_2ic(z)
CA_i=0;
' O2 {% v6 b. b9 g9 J1 `# ]- L: j%************************************************
% 边界条件函数
$ S6 o! l4 R  e8 H9 o6 d; m%************************************************
7 D, z/ U" Z0 l- _1 g3 Z6 d% \function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
+ p7 j) D( e3 s/ q: B6 t; q/ Mglobal DAB k CA0
pl=CAl-CA0;
ql=0;
+ h$ v' w8 A: u/ @6 B  l5 Zpr=0;
% K8 q0 _& E' ^: X9 }" h4 t7 a$ Cqr=1/DAB;
7 ?8 Y5 o+ ~- \/ y$ h( Y
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
  _9 N8 X3 o# o- }原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
7 K5 Z4 v8 L# i9 d6 Z

- ~8 m, l/ W+ r