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

浅夏110

8 D. V2 X6 }, }- _3 g题意解析：

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

1 {' h+ [! D+ G9 I" y  Y& p( v

1 ?" I7 h" e1 d% |% }(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为

" N  B9 P( i$ X
& j8 S- b$ Y9 I" R+ ]8 s
9 {; c5 R; S1 [- o; M5 I# r  g" M" u% B而起始及边界条件同上。
) N2 F3 u7 L: o1 i
7 ]5 _, U% z! O* `' |在获得浓度分布后，即可以 Fick’s law


& p7 `; T* M" D4 X8 s$ I: C' n
1 ^) e: I) q9 B$ K! W计算流通量。
, p' t& H3 S( W* M( g& \- w
) n# Y# e/ k! BMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
$ q! D0 a- ?( a3 Q


, z) R  M! K5 j" |. @4 k利用以上的处理结果，可编写 MATLAB 参考程序如下：

function ex20_3_2
: w) m: Q+ k4 `%*****************************
* @) T6 S- \+ a% 扩散系统之浓度分布
%*****************************
9 a4 B9 I' ?3 Q3 o* f) _clear
clc
' P1 u( E9 R% z$ |! Nglobal DAB k CA0
%******************************
% 给定数据
% S% B6 r0 Y% m3 H" Q! S. n%******************************
CA0=0.01;
' O* l) D$ ~) ]9 K5 iL=0.1;
- ^) Q2 C5 V5 \DAB=2e-9;
6 F+ I4 K% D- M4 K9 T& Mk=2e-7;
h=10*24*3600;
6 A; P3 f: G, H- m! N$ o% n%*******************************
' y9 B% u( J' N2 x7 ~% 取点
7 _- `/ E% e- z8 o+ F1 q9 j%*******************************
# P0 G- K* l7 at=linspace(0,h,100);
. B0 @5 V# j7 m6 ?) @z=linspace(0,L,10);
' K$ U' K# w: ~! c%*******************************
* R* r# R9 Q, X8 @% case (a)
$ S7 I4 d+ W$ j% E) K7 r%*******************************
m=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
5 [: W) z0 T5 n" eCA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
  N* H/ n" }- H3 |/ T/ d NAz(i)=-dCAdz_i*DAB;
. O0 I) ?3 k* ^* T* T! B( Pend
figure(1)
subplot(211)
surf(z,t/(24*3600),CA)
title('case (a)')
- B9 g) W, H0 j" }+ c5 l) Txlabel('length (m)')
ylabel('time (day)')
zlabel('conc. (mol/m^3)')
! ]9 ]5 _7 @2 M3 Q6 |- tsubplot(212)
, Y( Y! V2 f& b" V7 Kplot(t/(24*3600),NAz'*24*3600)
7 Z" J9 q" W9 E- M8 dxlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%************************************
% case (b)
%************************************
! F8 z' k, E& h6 g, _' um=0;
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
; ]6 n7 a# g2 I4 Efor i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
. ]4 o: ^+ X4 R0 xend
%
% ~* x9 Q- c, y( Z2 f, B/ B# Xfigure(2)
8 }0 |2 V7 y7 L3 _subplot(211)
surf(z,t/(24*3600),CA)
title('case (b)')
" n7 K8 J4 [0 R# Z/ f) U2 N/ mxlabel('length (m)')
ylabel('time (day)')
zlabel('conc. (mol/m^3)')
! V9 }/ ?2 ^' ~subplot(212)
$ L+ o2 ]' H9 I% \/ y) @' Hplot(t/(24*3600),NAz'*24*3600)
0 V4 @! X- U9 [7 p4 R9 C* Mxlabel('time (day)')
. B6 m: |1 {* S+ G/ I( i9 k& n1 r0 w) Vylabel('flux (mol/m^2.day)')
- S( Y1 I0 _1 |+ ?# _%********************************************
% PDE 函数
%********************************************
% case (a)
5 D+ V: j4 I5 k2 f# x%********************************************
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
; h1 r0 l7 l* t9 uc=1;
$ U$ m) m% O& F  ?f=DAB*dCAdz;
s=0;
%*********************************************
% case (a)
%*********************************************
3 f6 F7 g9 o' H% N/ vfunction [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
3 ~: B* P  D2 Y5 ?# }6 m3 Eglobal DAB k CA0
! d- i2 q2 y5 E+ ?) h# d8 ]& fc=1;
f=DAB*dCAdz;
/ i  |& ]; C) W, ^% e: Y  Qs=k*CA;
%**********************************************
% 初始条件函数
%**********************************************
function CA_i=ex20_3_2ic(z)
CA_i=0;
%************************************************
6 t2 R$ F( `1 a: ~8 q* g; q% 边界条件函数
%************************************************
* h( l, b7 F4 ~1 Y2 _. ~  H, V( `function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
- l' j; }, ?& spl=CAl-CA0;
3 k& Z: |, w3 [( [- ~2 Fql=0;
pr=0;
# c6 L6 Q" r: I$ `* nqr=1/DAB;

————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
  }2 q5 l9 n% f9 o& e( p原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
7 ~) E3 D* l4 I$ f, C

2 g; m5 c% |9 X/ Q