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

浅夏110

2 o/ j- ]( Y  C8 C9 ], y题意解析：

: k# h8 U  w# J0 B9 r' L( M7 q( J(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
) l! n# ]) v: \8 L" ]1 j# R
' l& B5 Z' W# t' K$ H

(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为


& F! o1 p, D" w/ J# a6 r
$ j5 r! {+ Q0 p2 J7 P& A8 x而起始及边界条件同上。

在获得浓度分布后，即可以 Fick’s law
: p4 R2 p0 m2 m4 K& S- T# {$ ^8 ~

, \% T" i. n+ N3 p; }( R* y( W- l
计算流通量。

6 j) N$ ?! H( B+ F. B' [& eMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下



. }, L6 W  j7 u. ~4 M利用以上的处理结果，可编写 MATLAB 参考程序如下：

function ex20_3_2
" G) ]; o  ]5 S) _- f' b) m%*****************************
  I9 Y. H9 \- V- e# N$ Y% 扩散系统之浓度分布
$ r+ ~2 }- {/ J& O2 ^+ l%*****************************
clear
clc
4 P* J+ [! G3 q, w6 ]global DAB k CA0
4 P2 o. j5 b4 S) J0 z7 P%******************************
* e; ?2 N" L/ s: u0 w$ u% 给定数据
%******************************
: }1 N, F; c9 ~CA0=0.01;
* i( m: Z8 t7 L7 L7 HL=0.1;
DAB=2e-9;
k=2e-7;
1 W1 t5 Y+ |, B" W# E. mh=10*24*3600;
%*******************************
3 l3 W% Q3 Q' e( q2 k% 取点
%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
%*******************************
% \/ N) N4 T  @5 z0 U. g% case (a)
& o; m6 K% x7 s1 Z' b! c%*******************************
. J& n" Y' t- x0 j2 E+ pm=0;
) A  q2 N+ X4 W' d2 _8 I: \, jsol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
4 u" w+ Y* j, D" }! t; M, Qfor i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
NAz(i)=-dCAdz_i*DAB;
end
2 Y$ V' C2 G  Hfigure(1)
subplot(211)
8 ?+ m- V3 P) C9 U0 ssurf(z,t/(24*3600),CA)
title('case (a)')
3 t1 Y- o% e8 x; C9 l7 Gxlabel('length (m)')
4 _' @- x9 q' i6 F' U7 @/ Wylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
9 k& |! u  Z' N. I2 t- q. b' n3 pplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%************************************
% case (b)
! `0 z* v& A2 @%************************************
m=0;
9 v7 c- o6 [% F# csol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
9 @+ ], P, Q: \; _6 PCA=sol(:,:,1);
2 P4 [0 Y- t! y/ m# Q" k. C1 ?for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
4 x7 O4 V2 _4 J, o, ]$ O NAz(i)=-dCAdz_i*DAB;
& T8 ~" M; ~% Xend
%
figure(2)
subplot(211)
surf(z,t/(24*3600),CA)
title('case (b)')
xlabel('length (m)')
ylabel('time (day)')
5 |, h$ Z' k. q# g/ h: U2 w: X4 Ozlabel('conc. (mol/m^3)')
( [* n7 [0 e: V- i& esubplot(212)
plot(t/(24*3600),NAz'*24*3600)
9 c- J5 C7 Z- o( F8 m" gxlabel('time (day)')
! D: T6 Y8 v% M! z+ Yylabel('flux (mol/m^2.day)')
) x; e8 g" N, h, B0 m* `$ j%********************************************
% PDE 函数
4 Y5 i) `% ^: P, a) J9 o$ o* b9 z%********************************************
% case (a)
%********************************************
( I9 I' v2 s, r% Lfunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
1 n- b1 H/ }3 X8 X$ t: \& Oc=1;
f=DAB*dCAdz;
" O6 C. n7 G2 [& Q- is=0;
- O0 i, D1 N9 P, d4 T%*********************************************
$ S, l+ o6 ?$ X+ a3 A3 c. k% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
' F9 Y0 t$ u* b! F1 l; `s=k*CA;
4 e$ y/ i# R+ K%**********************************************
% 初始条件函数
%**********************************************
' g% C- P% G, u' g2 T3 q3 t' p4 Hfunction CA_i=ex20_3_2ic(z)
; Q0 A9 E8 {+ i6 ?CA_i=0;
%************************************************
% 边界条件函数
%************************************************
0 ?' M2 o# Q, Pfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
! C0 r: V* {) e# J0 ]. K; S# cpl=CAl-CA0;
: t6 n0 }$ K  b3 ^' ]& ?6 h! Yql=0;
2 \( K8 O4 G* T5 j4 K: T% k- Ppr=0;
& T  Z3 ~$ x) T2 \8 w: a( p( rqr=1/DAB;
7 S8 _( {) x2 u" B, V
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
2 \; k+ f& w( O1 Y原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694


0 s, M1 G- j  U4 j