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

浅夏110

题意解析：

. w1 \, }7 `: W+ G' F0 C(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：

& S9 D7 c$ j8 w& X3 q
; k% }& z: h9 Y) B
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
! _) A. l- v# B# J" X1 T

  {) x0 R4 r) O: X0 F
而起始及边界条件同上。
) B/ ]' X: {0 n
0 ]0 D4 K6 t) {1 R( @在获得浓度分布后，即可以 Fick’s law



计算流通量。

# ?# u5 \6 c5 FMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
6 _# n2 l- x9 A" O# D2 P

# D4 L6 }! b; f& b2 _
利用以上的处理结果，可编写 MATLAB 参考程序如下：
& C  x3 J4 S9 p6 R
4 Z6 @: P& X. M3 Mfunction ex20_3_2
) j3 Z9 P# E. w# v%*****************************
% 扩散系统之浓度分布
7 I- {1 z1 ]; s% m7 I%*****************************
9 R$ D6 H! p( w/ Z3 d. d. Pclear
  d- j  o& Z4 @7 c7 I; Dclc
global DAB k CA0
%******************************
) Q: O& X0 \% R, q0 ~% 给定数据
3 a9 B2 @3 g* ]; m1 x: V. A8 @3 k. U%******************************
! S, P7 |; e" \" V  i: mCA0=0.01;
$ G7 v9 J7 W' k+ \# x% S; K  v5 f, BL=0.1;
6 w1 |3 n! J, c; N2 a( D' {. F& c/ PDAB=2e-9;
. V9 I4 r# \9 S* {: Mk=2e-7;
1 r  s; \8 Z5 ^  u2 y) D+ lh=10*24*3600;
% x# M& H6 l* n: p. g. o* K%*******************************
' P* f# O; l) ]. z0 k% 取点
%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
%*******************************
6 K3 i- Y; i9 |( `% case (a)
" r; y9 j! X4 B0 d3 h%*******************************
m=0;
2 d5 S( U7 {; B* x: c; P, r' P6 vsol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
. W+ Y$ `9 P& A: i) }4 t2 f* ?( ~- T) s1 pfor i=1:length(t)
' z; ?; B* Q3 d) W# C  K5 l [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
7 O, Z! o( ^$ ]. H  O" e' s! ~: Z9 Gend
figure(1)
( b* m- I% J! psubplot(211)
surf(z,t/(24*3600),CA)
; M7 v0 q' s* ?6 vtitle('case (a)')
! ]8 I, V# L5 v0 l0 S' m/ hxlabel('length (m)')
/ O7 C$ B6 q9 g, d5 vylabel('time (day)')
zlabel('conc. (mol/m^3)')
' j1 F7 l* Q; m$ f) N+ Ysubplot(212)
plot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
- G5 d! e8 R5 Uylabel('flux (mol/m^2.day)')
%************************************
% case (b)
%************************************
' z7 N% M. M* }* ]m=0;
! }) W: l* S* c$ \9 s' K+ ^sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
" m: D4 C9 z) b. N" sCA=sol(:,:,1);
for i=1:length(t)
1 w. R! z% B; r" k+ t5 Y/ E4 }9 [ [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
/ f  t/ n; v3 I3 D, c  g8 L NAz(i)=-dCAdz_i*DAB;
end
9 z! A6 ^9 W  ]% J6 E%
) a2 l0 y4 ]( b0 ]. z$ O& tfigure(2)
, ~4 O  Y4 o- s4 q3 Z1 ?. Q# W3 ~) hsubplot(211)
7 `) {$ c" x9 A) H% Q* hsurf(z,t/(24*3600),CA)
title('case (b)')
. r1 k2 q+ b# E/ X6 |# p! ^xlabel('length (m)')
7 ~  w8 S  K$ wylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
, |3 z4 w5 t+ ~- y# K9 \plot(t/(24*3600),NAz'*24*3600)
* g0 S- }* e- zxlabel('time (day)')
3 L$ L2 g7 a! b. M. Qylabel('flux (mol/m^2.day)')
%********************************************
% PDE 函数
- B0 Z* G0 i% z0 ^) [, o! h+ S%********************************************
% case (a)
6 A; |( u9 f  O; l%********************************************
8 k) ~# i0 L8 `function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
c=1;
) I9 ~! w4 _* |f=DAB*dCAdz;
- E' _* v/ ~! m( A- z  Js=0;
%*********************************************
* W7 f4 }: @( M5 C6 q% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
1 [7 i, ?4 o1 X2 h0 G# Pc=1;
f=DAB*dCAdz;
s=k*CA;
0 a/ A% \. f' [, X1 Q7 Y$ ^4 G+ @%**********************************************
% ~+ F0 l& T! u# n3 \& u, i& o  J" u% 初始条件函数
%**********************************************
, w# R' W+ K1 J- S  y2 C) Nfunction CA_i=ex20_3_2ic(z)
CA_i=0;
%************************************************
% 边界条件函数
! c9 o2 M1 o$ Q9 v( ?%************************************************
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
pl=CAl-CA0;
ql=0;
pr=0;
$ A5 R" j2 W  z$ l; l6 xqr=1/DAB;
- L) J; M) v; u" ~9 D# P
————————————————
8 c1 ]& T4 f: ]9 k1 d版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
5 D: r0 h6 j: U( M

* s2 B+ T. o) s7 i3 D; j, ]9 Q