数学建模社区-数学中国
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:29
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
" n4 b( [- a0 W$ L8 v, [
题意解析：
- d: T& Y1 P* D3 E2 g, Z' ~. O1 I1 v
% {. ]+ Q7 e3 [7 |9 o7 h2 K
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
2 H0 M9 |6 B3 l1 r! F
j3 q$ n; m4 i
4 \, b2 w& J( b( n
- A! d; M! }+ }0 D$ Q. z3 w3 W; e
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
# R/ {1 w) ^* P' m5 X" ^& r
7 K0 F9 r( J) v' W: `1 C2 m3 [
2 Q. w. e+ F! |: F0 x$ A+ L
/ l: X8 B' r* W4 d" |
而起始及边界条件同上。
! U; S6 P2 u. n+ q, K8 j
8 e0 Q0 f4 O/ l% ]3 z
在获得浓度分布后，即可以 Fick’s law
3 s) k3 s. I' G6 g3 g0 R, ^
1 J `" @0 E8 `* j$ \
+ R* R8 [( W" o" K5 ]
5 F3 I, H/ H6 x' _
计算流通量。
/ @- F5 T" J$ U
0 K8 |" b2 ?4 T5 ^
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
+ [4 h1 k, H: Y/ j- [' `
" P6 m5 Q5 ^! C7 j2 P6 u2 ~
! R0 o4 o- I5 M# v! |% V
( ]: U4 Y' d2 I; H0 V q
利用以上的处理结果，可编写 MATLAB 参考程序如下：
0 d$ ?9 D! L4 W: m. |3 \
$ K" `. Z7 I4 H' P* J4 K
function ex20_3_2
; ~ o/ A$ b0 L! t
%*****************************
9 I( q# t: V1 T3 s- ]2 u/ V4 t: Q
% 扩散系统之浓度分布
: m( f% K7 c- O p$ K2 i- j8 G7 Q
%*****************************
2 u5 p' W2 i0 R& [# j: C9 }
clear
3 j0 d3 W" t, e/ p5 X3 z( w% ^9 X
clc
2 z5 w! ~' S8 [1 V J( C
global DAB k CA0
: D! p$ l1 R! [- M' d# m* d
%******************************
6 V7 j( L" G, ~ T
% 给定数据
7 u2 X$ e' X, L4 J8 k/ d' w% N Q
%******************************
1 j: J) m4 V# D+ f, J1 F
CA0=0.01;
* }% Z& r- W& U; N
L=0.1;
9 ~1 m3 @9 ?3 r- s+ I
DAB=2e-9;
. y& u. }: W- ?* @( w1 _
k=2e-7;
# A. v2 W+ k7 q
h=10*24*3600;
4 [1 ~% t/ Y* A) s: T) Y; M
%*******************************
4 D! U) j0 E: h1 P% R7 C: l+ I
% 取点
" ?) _+ H5 H0 G
%*******************************
- `9 D2 T7 t% N/ r( m; i9 Y
t=linspace(0,h,100);
0 M# f' F) x# h* K* [% H
z=linspace(0,L,10);
$ s: v8 U J% S8 ~4 W
%*******************************
1 {6 v0 D0 f$ A: D
% case (a)
3 V+ i1 m& M- o3 h H
%*******************************
; i- i9 f" d, f0 G& E N) q" X
m=0;
W9 S& Z+ {/ I$ l
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
/ q g+ d" I4 u0 N
CA=sol(:,:,1);
6 \) p: f6 p3 r9 V
for i=1:length(t)
8 R2 @: n& J8 H
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
2 H" m" D% q* m- M6 z% d
NAz(i)=-dCAdz_i*DAB;
, v- b3 ~$ j* X! R
end
9 J$ {4 f/ V" n4 H, s) H- e5 A: b
figure(1)
! {# l: a' z% S5 ]+ P
subplot(211)
2 j! @$ A1 T g1 D: v( D
surf(z,t/(24*3600),CA)
1 a5 v/ H* z7 O& O" \: O$ N% R
title('case (a)')
2 |9 C% e" N7 q% B! s: d3 I
xlabel('length (m)')
' U; j* f* ~7 p! @
ylabel('time (day)')
: J7 n/ x" A( Q& E2 V8 p) N
zlabel('conc. (mol/m^3)')
) s+ k3 r. h; g3 V0 n+ @& w9 J" l
subplot(212)
1 z! {2 [ A6 a5 d/ n
plot(t/(24*3600),NAz'*24*3600)
0 L& o6 C$ L4 y0 ]- Q- K) z; o d
xlabel('time (day)')
9 w: C9 v$ K7 `# ]/ y- ^# K
ylabel('flux (mol/m^2.day)')
# L7 l$ P* _/ |) _ u3 c
%************************************
8 V2 n3 i3 a4 Y; E4 W$ a- O
% case (b)
/ M! V9 @8 ^ K5 K d
%************************************
9 F, r& G; L& R _6 {( {5 x) c
m=0;
# M2 m; k {5 U9 G2 T
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
/ U* C1 f- l5 x( \( K9 p
CA=sol(:,:,1);
* G- I5 ~' Q0 z1 W
for i=1:length(t)
$ y( T! y& C9 N5 s N( q3 z. f
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
0 U) q2 M1 Q: L6 |) z) h2 y) C
NAz(i)=-dCAdz_i*DAB;
5 q9 a m7 @% [: U/ ~( m, D
end
: Y8 w4 @# t2 _$ @; D: E* x
%
/ B$ |3 ]0 V" A1 Y" g9 R
figure(2)
$ d! h. Y3 Q" d. F
subplot(211)
: [" p; [+ ]$ e8 U, P- k
surf(z,t/(24*3600),CA)
0 G7 C. v# U7 d+ x- j0 F
title('case (b)')
: L8 R& y2 W, y: |
xlabel('length (m)')
1 L% z5 l+ T+ \/ j# H
ylabel('time (day)')
: U+ o" s7 O* ?5 S3 I
zlabel('conc. (mol/m^3)')
2 X n W1 W+ r: T
subplot(212)
0 \! u( w# b! }$ |" K3 k" i7 h
plot(t/(24*3600),NAz'*24*3600)
& d2 A3 `$ K- u
xlabel('time (day)')
; U% L8 T1 M7 R9 k
ylabel('flux (mol/m^2.day)')
3 v2 o& O; u( o4 ]
%********************************************
( a9 \& x. L' _7 l
% PDE 函数
& i* N0 ^. h. t3 Z
%********************************************
# a3 X. j8 x6 P3 m8 }7 p K
% case (a)
) a# c: T7 Q6 Q" E
%********************************************
" R; ]4 U7 y9 z, e2 \$ s# O' j
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
6 u- d% E8 n: p$ t% f
global DAB k CA0
' _+ |1 N Q( z/ l
c=1;
: |1 ] b7 F7 n5 j
f=DAB*dCAdz;
6 f H) d; m! w3 j1 r
s=0;
Q! Z1 @% }% ]% b- X9 c
%*********************************************
! c+ Y; x# ~* q: M1 ?- y
% case (a)
2 |* F7 g( O, A
%*********************************************
5 k" Y9 ~2 T- Q- O& x" X
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
7 l! n. W! E( }4 k% ]0 e5 `
global DAB k CA0
; q3 X, J: ^7 s8 g2 Q9 {1 A, v& k- l; F
c=1;
+ [: `; v: \* M" @; Q6 F
f=DAB*dCAdz;
, e6 j; |! A2 O
s=k*CA;
; ]5 S2 Z( S4 k, g6 h+ m' ^
%**********************************************
5 x l( o6 `- V- A; O, x0 D
% 初始条件函数
& r) P; k+ O/ x
%**********************************************
3 }! S/ y4 W5 j7 Y# B" u( `
function CA_i=ex20_3_2ic(z)
' s% H5 c: C. Q6 }2 B
CA_i=0;
! d2 u; E B' \! ^3 W9 w8 J
%************************************************
4 T/ n% K6 ]% W; M9 w9 v$ d
% 边界条件函数
' x+ Z7 a4 w' w! `8 I
%************************************************
5 L5 V- O% d' R+ B# R* V$ r
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
+ U) c/ ?3 G- c/ J
global DAB k CA0
1 m5 ~$ _' @. |; R* x; D9 ~5 W) v, N
pl=CAl-CA0;
, r" q& K$ l" L* k/ U8 T" _/ z' j U
ql=0;
4 @$ g. s, M. x
pr=0;
, {" y1 X# C y
qr=1/DAB;
" k& b+ k0 U4 k8 p5 y# M: G
6 Y. B5 X8 d. @
————————————————
# W& [" D j5 O
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
4 ?! @, W# g; n0 R& O( B
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
1 o9 i) s6 I! X! ~
/ L/ F& l6 w9 W; c' }
0 o: }$ \9 |1 V4 G. R3 X
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
" n4 b( [- a0 W$ L8 v, [
2 Q. w. e+ F! |: F0 x$ A+ L
! R0 o4 o- I5 M# v! |% V
,0);,0);