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

浅夏110

* O' F. d7 q6 P1 \& g  x题意解析：

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


9 [/ z: |$ [5 R& U: ]
# |+ W" M6 m% q(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
8 j1 O  j* h! \7 C
# }: h3 @4 b1 S* |- D
2 W/ G4 v) O2 n3 j" y
而起始及边界条件同上。
0 P0 d" ^5 Z4 N0 `4 O6 m; _% ]' R7 d
在获得浓度分布后，即可以 Fick’s law
/ L2 b6 M( @$ |
+ {* \7 ^. a- y1 X, z$ o# k% n

计算流通量。

7 s' @2 {; R  o+ g6 K9 tMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
  `0 o' s' U: h+ a; A4 C1 R1 E1 R
( K7 t; Q/ ~0 _5 L* B
& L- A1 T) F% B( ?* E
利用以上的处理结果，可编写 MATLAB 参考程序如下：
8 g- J6 b# b- r( V; w8 G
function ex20_3_2
+ F3 [: }% X# T8 ~7 T9 O& A: _%*****************************
, ^# `. h9 w. f; ~+ p& S% 扩散系统之浓度分布
%*****************************
7 q# }# `" U. c# t2 g5 t5 d# v4 |3 [clear
7 R3 `0 \# ]) [. O' ?! m  iclc
global DAB k CA0
4 i* ]! e( `- B8 z4 T%******************************
; |: X) w/ a" ?! B( J6 W5 a3 d: A. A. x% 给定数据
6 r; p4 S) N, n6 h' B1 T# c%******************************
CA0=0.01;
5 h0 R3 P( n$ U, v# l8 dL=0.1;
DAB=2e-9;
k=2e-7;
) s, M$ {/ T+ mh=10*24*3600;
%*******************************
% 取点
%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
1 ]$ o& X& w  \+ F5 u' s2 s/ F%*******************************
1 \: e) L. T/ P' h/ \% case (a)
. \) t7 B; \4 k# m9 I3 O6 R%*******************************
! I9 C' q( E* i( w+ c3 Km=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
, A& N; q$ _4 G; |8 \for i=1:length(t)
5 W, }9 G* k5 l( |% s  n& Y [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
$ ~9 |0 h0 z1 pend
figure(1)
subplot(211)
& o) I4 T  X. Z: G& ysurf(z,t/(24*3600),CA)
/ {; i9 F5 _% e+ r: M! {4 e5 A5 D, ~& Ftitle('case (a)')
; Z* M1 t8 R% l7 B1 Fxlabel('length (m)')
ylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
3 j+ H6 F. K- X1 `3 R3 nxlabel('time (day)')
/ c. t; H: @9 @8 u/ p  Wylabel('flux (mol/m^2.day)')
%************************************
% case (b)
( c! h2 N' V! b( b: C6 v: i%************************************
m=0;
5 M( [1 ]/ c) \7 W/ _, U1 ^sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
4 g- @; `" |7 Z  @: HCA=sol(:,:,1);
: r$ S6 L3 P/ d& N- ?: [( d8 sfor i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
* J9 v3 M) r3 X! m2 d! h5 b0 T1 m NAz(i)=-dCAdz_i*DAB;
end
, ?. A0 H6 T0 v& W- R%
figure(2)
subplot(211)
0 ~  N' M2 V2 ?! q* [8 R" N( ksurf(z,t/(24*3600),CA)
% l2 X. W  L- K) ]title('case (b)')
xlabel('length (m)')
, n. d- ?% a! O7 W1 D9 L6 Cylabel('time (day)')
) x# k8 K* a& M7 Izlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
6 `8 u6 m7 b$ q$ T2 _ylabel('flux (mol/m^2.day)')
%********************************************
- p  j: q" L7 H8 `1 f% PDE 函数
%********************************************
8 _4 r* y4 z( a  n' m8 X! ~% case (a)
%********************************************
* I  Y6 F2 l$ y5 yfunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
/ y# G6 B# k. s; Jc=1;
f=DAB*dCAdz;
s=0;
%*********************************************
% case (a)
%*********************************************
2 f+ `6 H+ K) E; x% tfunction [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
: x0 y: s- b7 n; V3 P$ Iglobal DAB k CA0
9 P0 H1 ~1 r) R4 l! i2 Y, oc=1;
) C" ^$ ~0 G9 z! Q- u8 uf=DAB*dCAdz;
s=k*CA;
%**********************************************
4 {# C7 M7 W! O: Z/ J7 Q% 初始条件函数
  P/ ?/ r! d1 {# {2 I3 p. W& V%**********************************************
. o& |& T! g2 ]& r( wfunction CA_i=ex20_3_2ic(z)
CA_i=0;
% q6 s+ \7 L2 i7 R& y" [/ \%************************************************
% 边界条件函数
%************************************************
+ j7 A$ U5 R* o# cfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
% v; w% }1 n$ r" ]  F7 dglobal DAB k CA0
* L) m) x8 y8 P5 g; Cpl=CAl-CA0;
ql=0;
8 ^3 J. y, Q7 x# @& Jpr=0;
2 X  z! ~! H, m) U5 V4 \. K2 tqr=1/DAB;

5 m; b8 k4 T6 Q9 v8 }: Y% `————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
$ l7 }4 H! g. `原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
! n5 X1 t* I) I
: j/ c8 v6 D9 g
& ]1 g; X8 y, `$ I$ \0 K