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

浅夏110

" G7 s: A5 o+ p& g题意解析：

(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：

- O3 o/ X% s1 e- R( ?7 w. I

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


' s& M+ I' E# `4 ^/ ~$ H. g3 d
而起始及边界条件同上。

在获得浓度分布后，即可以 Fick’s law
* b+ f( r% u% E# g8 ~


计算流通量。

MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下

8 D1 [2 B- U# ]' v) w$ Q

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

function ex20_3_2
2 R( v$ n1 U& U* D%*****************************
% 扩散系统之浓度分布
2 A2 Q" R" W' ?/ M" C%*****************************
- I# V) j: J0 ^clear
+ I8 C& q/ [( y9 O  f( Pclc
4 l( k) g( M9 o* s6 j9 Vglobal DAB k CA0
%******************************
% 给定数据
: X1 c2 \0 l5 Z4 p%******************************
CA0=0.01;
L=0.1;
" ^3 s& K5 n' x' K- \& v  u/ t. M% `DAB=2e-9;
% x- J% x0 h( |* ?k=2e-7;
2 `% Q. W" ~  zh=10*24*3600;
%*******************************
% 取点
%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
%*******************************
% case (a)
%*******************************
m=0;
3 C8 @' z! O- V) V5 x2 t+ Ksol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
/ T# \$ ~0 f" h: j; c; X3 cCA=sol(:,:,1);
( S# D3 v) X6 `7 F/ cfor i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
# A: W: T. F% b  X8 Wend
, S# Z% S! P+ a: B! T! G, d+ E7 E2 bfigure(1)
6 t6 S1 n% u& v' S8 ^+ j2 tsubplot(211)
: R2 p$ x' T: f2 s$ t: l$ Q& Osurf(z,t/(24*3600),CA)
, V% @3 f, o9 i' N) M( N) Q9 Etitle('case (a)')
xlabel('length (m)')
/ |. I/ Y) k# ]3 N& I# }- oylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
; @2 ~1 o# Q* b) p* B4 kxlabel('time (day)')
, p! D! X7 @# |. L! W/ Bylabel('flux (mol/m^2.day)')
%************************************
/ i8 a; I* w4 ^% case (b)
0 E5 |/ Y/ S8 c4 J- p( c%************************************
7 o; }1 d, f, ~m=0;
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
for i=1:length(t)
! A2 O2 A- c6 _" Z  q. P [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
0 A$ `/ H9 B* s( q%
figure(2)
7 V9 n7 S* ]- T' F7 r  r. \) Isubplot(211)
- A+ U+ {. A$ G( Tsurf(z,t/(24*3600),CA)
title('case (b)')
; b# E0 n) m% O9 t' l' S8 U$ vxlabel('length (m)')
ylabel('time (day)')
8 h. A" i& U3 g1 @; X/ l$ e8 szlabel('conc. (mol/m^3)')
  R" p7 k. {5 s; dsubplot(212)
plot(t/(24*3600),NAz'*24*3600)
, \- W/ X1 r- K1 mxlabel('time (day)')
ylabel('flux (mol/m^2.day)')
  k# O' q6 k% b/ C5 {7 j0 p) w%********************************************
; p4 I, c& m6 q- b) o, I3 c/ `: y% PDE 函数
. [" k* `6 l# y+ H% _7 e7 @%********************************************
. y, f" w$ S% F( V% case (a)
%********************************************
2 u/ K. i6 t  Ufunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
s=0;
, X! d" S) l7 {5 A# e/ @; r# x%*********************************************
% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
2 `) a8 p0 J6 G% x" O6 S7 z& Sglobal DAB k CA0
c=1;
% @! w1 X/ ]7 K' Vf=DAB*dCAdz;
8 B3 K& C4 Z9 ds=k*CA;
$ Q  t) C1 N: s8 J. R/ g: C%**********************************************
6 R8 G) M6 W5 O% C  M% 初始条件函数
%**********************************************
: o- K; f# w7 W( D" }4 ]8 Afunction CA_i=ex20_3_2ic(z)
( V, h- a) f( T1 Y+ i  X6 ]CA_i=0;
%************************************************
6 ^, g$ f- `8 s$ v% 边界条件函数
5 {& n- o' W/ R- i%************************************************
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
. c/ a0 [$ o/ f8 B" dglobal DAB k CA0
pl=CAl-CA0;
' z. f* u& E* `8 w; z2 X  cql=0;
pr=0;
qr=1/DAB;

# J% v$ I) S! k; z————————————————
, ]0 n' r* q6 M" `版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
3 i' q4 p- t2 O' v- f原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694

- @1 K+ h: s7 N) J$ r0 |