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

浅夏110

题意解析：
; D" D  D7 G6 _# ~2 _3 x6 j4 y
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：

% |$ `0 u% @% a0 y/ j% M

(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
+ t3 n* |8 ?& F  j
1 q7 s+ U" l+ c: m+ j) h1 {

' P" g8 J+ l5 ^: `3 e4 _0 o( P而起始及边界条件同上。
3 V! y$ W/ c! Q- g5 I
9 \7 C) M* m, L. E在获得浓度分布后，即可以 Fick’s law
8 Q/ b( T! ?, E0 s( V. c

! i6 Y' l4 g  c8 _$ |/ w, `( C
. \& L& Z, [- i* S5 M1 v. j& n计算流通量。

* D% v  z2 |% N* NMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
/ s+ r3 M  H) g' ~3 T. @* k' B

* L/ _* v: i" w6 G% D
利用以上的处理结果，可编写 MATLAB 参考程序如下：
/ |; `6 h2 Z# f& P. W+ N: o1 M$ I
: W- O/ ^6 B& T& Q! afunction ex20_3_2
%*****************************
% 扩散系统之浓度分布
8 o7 j% S& f. C7 x%*****************************
( p9 l: g% @9 X) E7 Q& l2 [clear
clc
global DAB k CA0
%******************************
% 给定数据
1 r3 Z- h+ g) l, v4 P4 a%******************************
' s1 C9 {) Q4 G( K  dCA0=0.01;
L=0.1;
6 q; @1 I2 C" j* I7 P/ cDAB=2e-9;
k=2e-7;
* ^& \, ?; ^) ~4 D0 Yh=10*24*3600;
: ^! Q( f& w. X3 J8 |$ H. n%*******************************
% 取点
%*******************************
t=linspace(0,h,100);
8 j7 y7 o( G" H7 a* T1 Mz=linspace(0,L,10);
%*******************************
% case (a)
%*******************************
8 e# G5 I, C* sm=0;
* s" Y7 O6 I/ d* ?' g$ Msol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
2 E" j# A' J$ l$ s/ t' ECA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
. E- [. i; i1 ?. T1 Fend
figure(1)
- M' R) B; A; R2 Fsubplot(211)
) }+ @' H( V  osurf(z,t/(24*3600),CA)
title('case (a)')
* p/ J8 U4 O% ?3 lxlabel('length (m)')
5 O: }4 p' K; J4 q6 tylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
" f) h3 ^" e4 o7 qplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%************************************
" Q1 |8 x+ p% n- Q. P7 R% case (b)
%************************************
4 C" f0 h/ B7 Z: y% F# O3 em=0;
$ c) b/ _* x: ^4 A2 C2 O3 B* W9 [9 @- Rsol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
3 p0 P0 d9 V0 u) v, @$ I& q% PCA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
* y0 n8 D) V4 }; f7 i%
figure(2)
2 `- ?+ ~/ N+ k' _8 csubplot(211)
' T( z: F# b8 W: g. |surf(z,t/(24*3600),CA)
, k! N, f; M  P7 Btitle('case (b)')
' ?% B- {0 p6 j( k* Lxlabel('length (m)')
( Q7 g5 A' c$ i' v8 ^* Nylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
8 t0 ~" z, g! [! m" ~& t5 {3 xylabel('flux (mol/m^2.day)')
%********************************************
% PDE 函数
%********************************************
% case (a)
%********************************************
7 b+ g/ o9 R' P4 Hfunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
% m: F, C4 E9 S, u0 z3 x  U& Vglobal DAB k CA0
6 u  g4 N$ F  M+ `c=1;
f=DAB*dCAdz;
s=0;
%*********************************************
* u1 D6 q, f% h# x4 p6 j% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
  M7 Y" X1 G1 Z# }c=1;
9 M: o6 K2 Q- sf=DAB*dCAdz;
s=k*CA;
0 s& K' t7 S) E%**********************************************
; `1 c1 Y* C/ n  X- }. f+ B% 初始条件函数
%**********************************************
4 [7 ~" o( d; j, s/ \7 ifunction CA_i=ex20_3_2ic(z)
* [- o5 _1 W' C5 H- k8 [% cCA_i=0;
%************************************************
% 边界条件函数
%************************************************
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
pl=CAl-CA0;
ql=0;
pr=0;
qr=1/DAB;

5 p9 z/ R4 ]9 Y( r6 a1 W————————————————
6 V* N/ J: F( L1 p! W# s$ A  j版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
4 e& R2 V: U5 g: D
; `/ E9 H0 o  R' G% f% O$ X  w
$ _8 `. ^; H9 B