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

浅夏110

4 l1 \' O* z9 b题意解析：

(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
! @! Q5 n/ h% T* o6 k0 w

' ~8 l2 i/ p4 X7 B2 \
$ i) ~/ C# N7 N7 q(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为

7 _, ~5 H% H5 i; G# ^$ B

而起始及边界条件同上。

在获得浓度分布后，即可以 Fick’s law

! ]; I2 D1 ~% C6 ~: b1 B# x3 S
- ]1 A+ U( l6 Y. [' b; H$ L
0 d9 M' _1 x: \) \: Y' e计算流通量。
* r3 N6 I: D+ @2 v9 v
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下

2 s6 ]9 P' m8 I+ g) t* f
/ n3 {# ^" G% l0 _
0 {. `' S( q# O  Z; Y利用以上的处理结果，可编写 MATLAB 参考程序如下：

function ex20_3_2
%*****************************
: C, S& C: b  x- J* W  f% 扩散系统之浓度分布
%*****************************
% U0 @% x" I# x4 y7 vclear
3 B% x% D! n  `/ }4 n4 ?clc
global DAB k CA0
* S  Z+ `/ ^4 S%******************************
$ {0 j* n1 X, Y% 给定数据
5 G; n- }2 n6 L%******************************
CA0=0.01;
  p6 ^2 A$ \1 f/ q5 }) r# iL=0.1;
DAB=2e-9;
+ _0 ~3 p- a8 a( e# hk=2e-7;
h=10*24*3600;
- W' ]4 C) J% Q% v+ ?%*******************************
% 取点
%*******************************
t=linspace(0,h,100);
  Q; b$ F' S0 l- B$ Uz=linspace(0,L,10);
: p5 c- o  R' ]* y3 _) z%*******************************
% case (a)
* j6 M, o' R+ l! h%*******************************
2 K4 K4 J( a3 t4 G$ um=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
1 M- S0 f- I$ P- o. S! I; ZCA=sol(:,:,1);
for i=1:length(t)
0 S1 s9 p+ Q& ?8 f1 |4 A7 V* |% K; m [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
' M& Z/ b" q' z8 V NAz(i)=-dCAdz_i*DAB;
end
figure(1)
8 t* k, F% L/ c% |+ `subplot(211)
- f; u; g$ z: N' z6 dsurf(z,t/(24*3600),CA)
title('case (a)')
xlabel('length (m)')
9 z% J1 b5 ~1 a7 p- pylabel('time (day)')
zlabel('conc. (mol/m^3)')
/ X& \: k0 G, r  Q" Nsubplot(212)
plot(t/(24*3600),NAz'*24*3600)
" h5 X0 g8 x# h$ ]xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
% F" C$ E8 N$ H: z9 Q1 V8 L: Z%************************************
% case (b)
%************************************
m=0;
; T5 r& s+ e4 u+ msol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
+ D' k, C( _2 _5 oCA=sol(:,:,1);
2 {4 }& z8 h: i* [9 p! E% X8 ]. c& dfor i=1:length(t)
0 v( F' ^1 Y/ h" W9 T' ] [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
$ w2 {7 Y; v9 | NAz(i)=-dCAdz_i*DAB;
end
* ?% h0 p4 v4 {- t%
) d4 S% C$ n1 p* \$ }figure(2)
subplot(211)
surf(z,t/(24*3600),CA)
title('case (b)')
) i% k+ i" _; }5 W# ?; A) M/ exlabel('length (m)')
ylabel('time (day)')
8 S$ s5 l" Y+ \' M6 fzlabel('conc. (mol/m^3)')
subplot(212)
- Z9 g) ~0 k+ L; c& i" ~0 ~& S4 splot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%********************************************
% PDE 函数
%********************************************
% case (a)
%********************************************
+ I$ G" N; ?& E2 Q- ]function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
) r( N+ w6 G% ?c=1;
  A: K' \  P5 C& s: x! d& cf=DAB*dCAdz;
s=0;
6 M/ a9 x+ t9 p%*********************************************
% r8 I0 h( C5 e8 c7 y) I+ `1 C% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
s=k*CA;
%**********************************************
/ e4 H6 {- R- v2 D& l" E% 初始条件函数
: ~  o8 z, g- t% S4 \# L4 m%**********************************************
! m; P- r; c+ _# t6 [1 Gfunction CA_i=ex20_3_2ic(z)
: x- U1 l. f* h8 ]CA_i=0;
9 w9 `# M: n' U* z8 j%************************************************
% 边界条件函数
$ m$ G! @: W5 m; w* U6 Q%************************************************
/ a( Y3 ]" ~3 Afunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
' \7 p" ~! ~( D& H! Upl=CAl-CA0;
ql=0;
pr=0;
: S' l: l( T; X8 Iqr=1/DAB;

& S+ Y# i+ F9 C& }. R————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
8 {: I  B9 s$ X$ t8 B+ X5 O! z

2 x, [7 d; {7 w6 T9 Z; e- F0 }# b