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

浅夏110

题意解析：
- @. `" Q5 ~3 t  J6 E. E# n: O
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
/ T) P8 Y' p- I' K2 y0 |8 b$ q


5 m1 h2 k9 K: ]1 P$ x6 Y(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
' m# G" P# @, H3 j+ A3 I

' q" ~* _; ^+ t
而起始及边界条件同上。

' D1 Y- D- L: Z' x! T% \: R4 i在获得浓度分布后，即可以 Fick’s law
% f, g8 s' `! {! D# i


$ E7 v* Y% W" V3 _6 ^! C" V4 {计算流通量。

MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
, q  M, t! z. e7 M


7 j2 |$ S( }- q/ d0 M( u! p/ K利用以上的处理结果，可编写 MATLAB 参考程序如下：

6 K3 A6 D$ s( x4 _# afunction ex20_3_2
; ?! g2 U, x9 M%*****************************
7 W1 e. s8 T7 J# R- {3 X# E% 扩散系统之浓度分布
, g- m% ^$ c8 e9 \7 K. ~5 q' g9 K9 e4 m%*****************************
clear
clc
global DAB k CA0
" p5 ~' ^6 B! g" D%******************************
2 n9 W: ^$ t/ a+ c" ~& d2 s% 给定数据
& f  A$ a4 K6 P, Q/ B%******************************
CA0=0.01;
+ f9 q) |% h& j+ I& @9 K( q3 cL=0.1;
; n7 _8 y: l" C# q' PDAB=2e-9;
9 e, Z) o8 L$ uk=2e-7;
5 y% S. e' |8 G% D0 |h=10*24*3600;
%*******************************
5 w  j/ ?; F7 C* D) m% 取点
2 k9 R* R+ U7 X$ I& R: @%*******************************
; J4 h8 i; _" D: z4 k- Tt=linspace(0,h,100);
% M3 p2 r) m8 }; \. c. v9 \z=linspace(0,L,10);
$ P4 P0 ~2 C/ \& `% [%*******************************
9 K  N2 T2 X  i: v& R- _% case (a)
2 n% Z  e/ w: s* e- X%*******************************
m=0;
sol=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;
+ |: F% f' P9 b7 b# tend
figure(1)
subplot(211)
surf(z,t/(24*3600),CA)
/ _( v+ D; y8 l, ~# qtitle('case (a)')
" r( N3 }5 T( G8 w+ M8 Vxlabel('length (m)')
ylabel('time (day)')
5 ]. t7 F( z1 i- o2 y4 \8 ezlabel('conc. (mol/m^3)')
! @0 Q% r; T1 D3 \; r0 o/ L: g' G0 ksubplot(212)
3 V# m) h' ?9 tplot(t/(24*3600),NAz'*24*3600)
& W4 q& e2 @% \% i& Fxlabel('time (day)')
% N; a) e! ?& d0 d/ u0 R$ y$ h  _2 P% Iylabel('flux (mol/m^2.day)')
%************************************
% case (b)
2 E8 g" q7 z, K8 ^%************************************
m=0;
7 ^! J/ e% v, e/ i- W, ]/ J( D' usol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
for i=1:length(t)
4 S+ T* q4 T5 R! N9 Y [CA_i,dCAdz_i]=pdeval(m,z,CA(i,,0);
5 w0 m: S) s8 C# y- x! w NAz(i)=-dCAdz_i*DAB;
end
%
figure(2)
subplot(211)
surf(z,t/(24*3600),CA)
& P7 R1 J4 M6 O# \3 Z$ Ptitle('case (b)')
xlabel('length (m)')
( G+ W  c- \& Tylabel('time (day)')
4 X. u7 h6 W8 ]6 P3 szlabel('conc. (mol/m^3)')
! ?9 R# A. ?3 l8 [! psubplot(212)
1 g9 J1 H; N' k; O5 d5 wplot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%********************************************
5 r( T8 O& F( _5 H( k) T% PDE 函数
% ~- m- Z9 @; Z* \# F5 ?%********************************************
% case (a)
%********************************************
( z2 X! E1 r7 t3 |* mfunction [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
& y$ c, p* y) c3 q6 bc=1;
f=DAB*dCAdz;
s=0;
+ j; X2 ?' ?9 `6 I3 l  a3 {%*********************************************
% g8 C% x1 t0 r% case (a)
%*********************************************
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
c=1;
6 Z# }9 _9 G" u" ~1 [& d: Bf=DAB*dCAdz;
s=k*CA;
/ j$ U: K" A: T/ E) ]" E%**********************************************
9 g+ B- t( b9 O% 初始条件函数
%**********************************************
; b9 Y) Q" W$ ]4 A# P+ H( qfunction CA_i=ex20_3_2ic(z)
! K& x3 J4 c& n+ ACA_i=0;
%************************************************
% 边界条件函数
%************************************************
4 H8 }% s9 b- g3 T* O* @1 a8 Yfunction [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
: f& O1 [; }& ~6 v/ f+ c0 q! q: d, wglobal DAB k CA0
+ u& a6 C+ U: |8 x  q: t6 Fpl=CAl-CA0;
ql=0;
6 T8 k9 u9 ^% \1 K' Bpr=0;
$ D- n& s; R8 \& Z6 N4 Dqr=1/DAB;
  W& G) P: R# O/ g
# \/ \0 V2 A- x7 Z' m0 H* q% `————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694


, T' F' i: H& S* U( ]- [