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

浅夏110

$ z5 l5 V7 W1 B题意解析：
, }7 |3 [3 Z5 j: H4 ]0 E7 V
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
& b9 K2 e+ u' W0 O  s+ s


(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
) X  S; E& i' l. _! N9 y# E
* |# N* `9 p4 y9 ~/ }! {0 l

而起始及边界条件同上。

8 ~8 r. c# s9 g* C, J# |在获得浓度分布后，即可以 Fick’s law



+ }5 y/ `1 S5 M/ K1 [) i6 Q计算流通量。
: I, i5 b) K' F+ d' ^0 ?5 c
- k" z( Q, K/ PMATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
+ E* D. r- Q5 o& C) k

1 N1 D( v  B* y% Y% @
5 r9 ]$ G' ~: l4 e6 w; |) m# R4 F) j3 {利用以上的处理结果，可编写 MATLAB 参考程序如下：
: a& \) U  a% b" P
4 A/ n. n  [9 f3 k4 A5 Wfunction ex20_3_2
%*****************************
% 扩散系统之浓度分布
' q0 r/ N9 L2 f( H/ Y" o4 G7 i%*****************************
clear
clc
global DAB k CA0
, J( g6 q: ?3 F- p%******************************
% 给定数据
3 b9 {/ C' L8 v; e%******************************
CA0=0.01;
+ }( q3 y8 S4 I5 n, w1 }L=0.1;
DAB=2e-9;
k=2e-7;
h=10*24*3600;
%*******************************
0 Z+ I3 a8 j: J3 r% 取点
%*******************************
t=linspace(0,h,100);
! M$ w& z, y, W, r) |z=linspace(0,L,10);
3 {' y5 c9 B" p: A* _& L  [* q%*******************************
% case (a)
; P1 _& Z7 R& v6 ?- P  s/ s( e2 f%*******************************
m=0;
- H0 K; a" h* g; g: v( dsol=pdepe(m,@ex20_3_2pdefuna,@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;
2 Y5 M+ ^  @- H! P8 Yend
figure(1)
subplot(211)
& I3 W4 {( S) }  K: [/ Tsurf(z,t/(24*3600),CA)
! o! J5 K0 ?" n5 a0 G- z# b/ Ttitle('case (a)')
xlabel('length (m)')
ylabel('time (day)')
' N2 n9 X4 h6 Lzlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
4 S5 ~; E! |* y* I! uxlabel('time (day)')
' ?  q$ t8 U+ j* B6 Iylabel('flux (mol/m^2.day)')
; \' Q8 ], g2 x' B) w) P%************************************
% case (b)
. p! z, L, k% n' X5 j& g& ~%************************************
m=0;
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
0 G& u, Q3 I7 W* Mfor i=1:length(t)
# B6 o7 y+ C, ^0 B; H2 T6 c, T& p [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
. x0 N; \3 t* n, @! ^  i9 y/ B NAz(i)=-dCAdz_i*DAB;
end
' m- k8 {: p. E* o! m; E- N5 I%
figure(2)
( x5 O: T" d0 B% `/ @( wsubplot(211)
1 T4 k9 S" R# G2 osurf(z,t/(24*3600),CA)
  y+ h$ j/ x' G# r5 gtitle('case (b)')
xlabel('length (m)')
4 G5 E( D( {) e$ X0 eylabel('time (day)')
/ W% |5 L, ]4 o' s$ D2 E$ Lzlabel('conc. (mol/m^3)')
! G$ g% U/ E$ S& A7 Lsubplot(212)
+ i  K$ C" {# K7 a( Aplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
+ [: q" ~7 Z, i' `- _ylabel('flux (mol/m^2.day)')
%********************************************
4 a2 h* s, B; }% C0 c1 s1 H, b% PDE 函数
%********************************************
/ D2 E' F7 }' q# d& U! G( Q% case (a)
%********************************************
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
+ M5 g' W1 `4 l& W2 lglobal DAB k CA0
- G  r9 V) L) O/ `" s: oc=1;
; d: y" r/ D3 W* L% R; _: h: P) e8 Nf=DAB*dCAdz;
% w& Y! B3 b% T6 j- u) c6 rs=0;
%*********************************************
/ f/ z& I( Z" n4 }9 x2 t& I% case (a)
, m: H0 K5 ^; a9 d! a3 j%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
4 Q3 t6 b7 ~" h) U4 Zc=1;
f=DAB*dCAdz;
, o  j* r1 [' E! b5 Ps=k*CA;
%**********************************************
% 初始条件函数
%**********************************************
5 U& D& J/ `5 t& G% ?& Nfunction CA_i=ex20_3_2ic(z)
CA_i=0;
) M. u" B  S$ @! d7 q7 p%************************************************
% 边界条件函数
%************************************************
( `9 r, g* p- e& `7 }/ Zfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
/ H3 J% X" B; O2 x  M3 O0 tglobal DAB k CA0
% E0 ~" Y) y% Epl=CAl-CA0;
  q% p# v- O7 ~ql=0;
5 `5 ?  \, c9 l& R8 ?pr=0;
* Z. V) ^: E" W( _) X  cqr=1/DAB;
0 j+ k# l( V9 b* t- M. w3 m
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694


7 p% v) W6 t8 Y6 u, a