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

浅夏110

7 M' `- \9 Q, s/ {- e0 O4 U  G题意解析：
0 r) S+ u/ h/ n4 X
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
1 ^8 `, w0 P, a
3 Z: b, M5 K1 C

, i3 M( `- Y3 n(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为

( |6 g1 j9 S' D. w* [$ o

/ ~- N* z2 M: S5 A% E8 n而起始及边界条件同上。
: o6 }5 y8 m" G
4 p7 o" E& G( `  z, H. @在获得浓度分布后，即可以 Fick’s law

7 b! s8 g6 _: T2 w5 k
+ S3 h* F$ }. p# N( t
1 k- r( Y2 Y! l0 K) f5 Z计算流通量。

) \$ e0 z& n: C$ O$ R+ p% L  pMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下

( Q: K; l+ E& G: ~. X! N' \6 M% l
+ o6 P$ _/ M& ]. F4 L
8 s5 G& Z0 z6 y' i利用以上的处理结果，可编写 MATLAB 参考程序如下：
: z( }0 ?: W0 F' P
function ex20_3_2
' e3 e9 g6 d$ D( g%*****************************
% 扩散系统之浓度分布
- f/ @" O( z7 S2 V" k$ }+ A3 \( e%*****************************
+ i+ v2 J& v1 S$ D) I7 b% S1 H$ Kclear
clc
global DAB k CA0
%******************************
1 T! H0 k3 t0 D% 给定数据
4 u3 u3 ]8 A: J%******************************
CA0=0.01;
: a) u, Y$ Z% J, J4 _+ b5 nL=0.1;
DAB=2e-9;
. e: X7 R* p) |: X$ _, x, u* gk=2e-7;
7 }8 P4 q6 _! m  r+ j  I/ Xh=10*24*3600;
! }/ r$ e4 k' a% h/ n+ P%*******************************
. }) V7 z) x1 U2 f5 T* Y% 取点
%*******************************
4 p2 v+ \+ C" \. l4 ?9 `& i6 ]t=linspace(0,h,100);
% Q2 e; [1 |2 Q% t+ A! ^z=linspace(0,L,10);
%*******************************
+ [' x5 Q- f3 E6 K% case (a)
* n2 W: w* V( m( s" O6 H3 D' G%*******************************
7 ]0 W. X2 ?8 ~  y9 |! c2 {m=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
. q' V% n/ b" q& x* R" [' I0 m; OCA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
$ P8 ~) x$ D; x* _end
figure(1)
subplot(211)
surf(z,t/(24*3600),CA)
title('case (a)')
xlabel('length (m)')
ylabel('time (day)')
zlabel('conc. (mol/m^3)')
! m( e0 Z$ z5 |& [' k) p- Bsubplot(212)
" K0 o+ `+ f8 u9 }, l% oplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%************************************
% case (b)
0 q/ `' U5 c9 ~%************************************
m=0;
2 W; U9 Q, V5 Msol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
; m: D* E, w" w" c( r3 Q2 XCA=sol(:,:,1);
: U) s8 @7 ^  Pfor i=1:length(t)
( j! }  O, r8 h& A- ?' C [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
9 [* F0 y& P+ p4 ~5 dend
% |( w; n/ n+ T! h%
figure(2)
) Q: w# l5 V. _; gsubplot(211)
  Y- L, m! `3 A# z; x' D- Vsurf(z,t/(24*3600),CA)
title('case (b)')
, r" H* y) W8 @+ R# o. L, Yxlabel('length (m)')
+ k7 A0 C6 u  B: j4 `ylabel('time (day)')
, D$ Y7 p( ^$ Czlabel('conc. (mol/m^3)')
3 Q, ]2 o" g# u4 o7 @1 csubplot(212)
: u4 O0 a- C- Z/ R: J9 v! E: U5 ?plot(t/(24*3600),NAz'*24*3600)
; ]$ K3 i6 [" ]/ g8 x* n( wxlabel('time (day)')
ylabel('flux (mol/m^2.day)')
% J1 D& U5 @3 A( f" K2 p) j%********************************************
% PDE 函数
%********************************************
6 {! O4 [5 l4 e0 T* [7 \6 r  q% case (a)
) P$ ^$ B" `! I" `( z%********************************************
$ ]/ ~) [$ b/ n' ofunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
5 ?0 ~# x4 W) Ls=0;
%*********************************************
& V9 s4 b4 C5 N" {: K7 ~% w& v- y% case (a)
%*********************************************
, W. N2 Z" a' p% j  ^! {; Wfunction [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
  `7 g7 }7 Z; \5 v- {s=k*CA;
%**********************************************
% 初始条件函数
%**********************************************
function CA_i=ex20_3_2ic(z)
+ i* H. l' w* y. l1 G. E# E( pCA_i=0;
%************************************************
* C; e% ~/ P# i" n. _% R% 边界条件函数
%************************************************
% E5 v( ?9 }; F5 l  [; Ofunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
pl=CAl-CA0;
$ x4 b. e0 x1 D" Z, z: ]' F: V$ xql=0;
pr=0;
qr=1/DAB;

————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
+ b! A0 m) V% ?  I3 d4 c4 G原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
6 R; X/ {$ s% q3 e