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

浅夏110

  N# u1 _3 t0 O  _0 r4 c题意解析：
* g6 m, d% k# A, I9 `
# c. y3 l  s6 q1 t) a(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
) h0 w3 g& l4 b' k
: H, F: r7 G: h7 t, I

(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
) q+ j( X/ O4 ?$ w/ H9 [. v" y
7 s  I3 a$ I9 b/ m
. Z8 h" m" ]; B% _
. ]$ p0 }' z0 [  b* ?4 u1 g而起始及边界条件同上。

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



计算流通量。

MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
# \+ k- Y: N3 l7 C% _
4 e9 P# q* Y& l, @) Z

& _. v5 _5 a" `% R8 F' V利用以上的处理结果，可编写 MATLAB 参考程序如下：
* S7 {8 C9 H, j% S1 h
function ex20_3_2
%*****************************
) a: P, @0 f5 B4 _$ G. @: I, t% 扩散系统之浓度分布
" G" c: I' L- U/ O& P7 Q. E* N%*****************************
' V# \1 i& K  Cclear
clc
global DAB k CA0
%******************************
% 给定数据
; w7 I$ e, m; u7 n. x' u6 Z%******************************
9 n6 [) O! I* xCA0=0.01;
) V2 j0 J' V0 j. l, c! n" kL=0.1;
, {9 L8 w9 a1 t/ @DAB=2e-9;
k=2e-7;
h=10*24*3600;
%*******************************
% 取点
% a6 _2 n, U; M1 x  K+ {%*******************************
t=linspace(0,h,100);
9 |2 q9 Z% H4 F. v, F4 a9 yz=linspace(0,L,10);
3 Y6 W* K+ m# f6 J; S%*******************************
% case (a)
. w* \8 s1 r+ H( t( A9 R%*******************************
$ E8 k3 x3 U8 um=0;
( s% v# q* d# S/ T8 Hsol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
, ~' n& L7 w  m: O2 V7 `CA=sol(:,:,1);
for i=1:length(t)
8 U0 D9 [9 j. [5 \$ }* m, X [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
) Q' E8 x0 X- z3 eend
- V0 L* R+ R, q2 V! T4 U! Gfigure(1)
/ \5 n( Q1 E. C6 a! S1 t( Bsubplot(211)
! {, `* f* Q% _$ C6 S, Ssurf(z,t/(24*3600),CA)
title('case (a)')
xlabel('length (m)')
ylabel('time (day)')
: t) j+ A" P6 r5 Xzlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
+ h; v& \* }4 F4 exlabel('time (day)')
' r! V, V% R+ D5 i7 Lylabel('flux (mol/m^2.day)')
7 ^. D/ D' F1 _1 G%************************************
% case (b)
1 M+ i9 u+ _$ R8 ?/ p%************************************
m=0;
+ S7 _9 r7 p5 o7 n7 P" Q+ Ysol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
3 R3 j$ B6 b7 e0 X! x9 g/ @; A& ^# ZCA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
%
figure(2)
- J) i* M: x4 isubplot(211)
surf(z,t/(24*3600),CA)
8 E" S5 N5 A$ D: }4 Utitle('case (b)')
xlabel('length (m)')
+ Y- F5 ?0 y8 p" L% S, cylabel('time (day)')
zlabel('conc. (mol/m^3)')
* v# f( n. Y* [) zsubplot(212)
plot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%********************************************
% PDE 函数
%********************************************
7 e) k; v2 x  [  K, @% case (a)
%********************************************
* Y  `7 v7 V8 }function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
( Z: t2 _5 {. |6 Y$ i- ?global DAB k CA0
- Q0 o! G3 Z: I, [  ?4 zc=1;
8 t6 A( u. h5 z) Z" [/ }; Gf=DAB*dCAdz;
s=0;
& o8 `2 D! I( ?0 I! Z5 ^%*********************************************
% case (a)
%*********************************************
1 v$ X: m4 \, ]2 R/ x" wfunction [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
s=k*CA;
%**********************************************
6 I6 f3 c$ P$ W% J; j! P4 s% 初始条件函数
! C  x. j0 B/ ^6 y/ b. c%**********************************************
function CA_i=ex20_3_2ic(z)
CA_i=0;
%************************************************
% 边界条件函数
( Y8 Y3 r2 {  G3 |% M1 w% R%************************************************
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
5 g' L% e6 o+ ^% u& hglobal DAB k CA0
pl=CAl-CA0;
ql=0;
pr=0;
! G: W) T* }  v7 E9 m2 i! |qr=1/DAB;

5 M+ Y; U! E, C' |————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
' X' O& r4 \" H" g原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694


/ u) ^) U) g- `/ z2 k