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

浅夏110

题意解析：

$ y: ?- w/ F2 p. u. T. l3 j(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
- J: O; o6 _+ p6 [) R
* u( @' [, `5 h1 Q

(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
( m3 f1 F* {# P4 f5 B% `/ T8 u

  j* T) y# @) s4 H
8 U' r) u$ T- f& q$ k1 i$ e而起始及边界条件同上。
+ W! v2 Y, y7 A% ^9 O8 R' r% _* U
( S" ]; A0 a4 y1 x在获得浓度分布后，即可以 Fick’s law


3 w0 `* Q4 W* C% D
计算流通量。
: N- @1 q* M5 f' Q. i8 E! c5 Z
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下

& Y5 j/ Q, u; }; O: i- ]
1 Z! m) A- r% r# o6 c- G6 J
+ M$ l. L8 g+ r  H6 m" p8 y; a2 y利用以上的处理结果，可编写 MATLAB 参考程序如下：

function ex20_3_2
) M) S5 N8 f' W: [# N%*****************************
% 扩散系统之浓度分布
; J; r! v+ d+ n. y%*****************************
' ]0 I1 u) Q7 }clear
clc
global DAB k CA0
8 V" b2 h! H8 L: y%******************************
% 给定数据
%******************************
CA0=0.01;
L=0.1;
4 u: Y. j+ o* e3 |3 ^- n- U: zDAB=2e-9;
k=2e-7;
h=10*24*3600;
- b2 H& \  K3 z  t9 m%*******************************
  y% r+ T4 O8 a# G0 v  {! l% 取点
%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
%*******************************
% case (a)
& _( N% d/ Y$ j/ Z0 O& V! z%*******************************
$ S8 e; [8 Z' Q* xm=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
; B* A! l$ J! C# a. z: g0 \CA=sol(:,:,1);
for i=1:length(t)
6 B, Z' e/ Z6 }; j [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
! X- u) A6 y4 E NAz(i)=-dCAdz_i*DAB;
9 `) l5 @3 b5 Z  b( \0 h- V. I. K* @end
figure(1)
subplot(211)
surf(z,t/(24*3600),CA)
title('case (a)')
xlabel('length (m)')
6 \, p! ]( w9 E4 l1 Dylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
7 S. o$ T' {' W. ?xlabel('time (day)')
* O, ^, d8 I9 g. k5 t& C8 ?0 w; dylabel('flux (mol/m^2.day)')
%************************************
+ v' F( U9 b# V6 Y- M" s4 c% case (b)
%************************************
* ~4 }  U# X; U" J, l! e2 Gm=0;
& l9 U) L1 ~! ~0 O) |3 Y' Msol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
. M; N5 `# f' Q* Y& q' U! H%
figure(2)
; m2 q8 L! U  R/ j. v+ H  q. nsubplot(211)
surf(z,t/(24*3600),CA)
title('case (b)')
* A. f8 }5 i1 Z$ }+ U3 F8 exlabel('length (m)')
ylabel('time (day)')
: L  J  f6 \& E2 {' ~7 R* gzlabel('conc. (mol/m^3)')
0 }) `1 g5 A  C. Y: Usubplot(212)
( [& G# u9 C% V6 splot(t/(24*3600),NAz'*24*3600)
3 a, }0 y: P# x# x$ ]0 u. nxlabel('time (day)')
+ n& K! {5 r: J: `2 h! Z' D7 D3 oylabel('flux (mol/m^2.day)')
%********************************************
% PDE 函数
%********************************************
% case (a)
%********************************************
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
7 V3 n, S! R2 z" K7 Uc=1;
# S# E6 {4 r* kf=DAB*dCAdz;
1 `- O3 l/ S1 b0 h' r$ j' Ps=0;
%*********************************************
% case (a)
! U6 A: T. R: `3 h%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
, _( E: [* J7 P& p0 {' ?: d7 Ec=1;
f=DAB*dCAdz;
s=k*CA;
%**********************************************
% 初始条件函数
/ q. y& q2 ~8 p; w$ P: i- Z%**********************************************
& m+ ~& {+ ?* m6 sfunction CA_i=ex20_3_2ic(z)
2 [" ?8 ?* c. N5 k" Q* s, w& zCA_i=0;
+ @7 x. M; g- f7 Z8 `0 k6 W5 h%************************************************
% 边界条件函数
%************************************************
8 q  K) E0 F( S) s: H% `# p4 {function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
5 \7 f' n9 T0 N1 }: |global DAB k CA0
pl=CAl-CA0;
ql=0;
pr=0;
' r+ w/ m- L  L( b! q" v( N5 V' Vqr=1/DAB;
7 q$ G# N8 G, C. D  Q
, y+ J  v7 [7 `8 B# ^9 I————————————————
; d. I& i% z7 s6 j6 Z版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
+ w  u# Q2 c9 g原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
; e  t' u' K) g+ S! N
8 m& b8 }: g/ q: K