数学建模社区-数学中国
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:29
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
* m2 @! @; `; ^- O m) S* N+ R0 M
题意解析：
9 u7 o1 q0 c' n
5 j) L; K" M. Y3 R) M
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
& Y4 J$ y; G# Z# p$ S
9 H) l) d" w- b7 f g, Q$ R( i; G
1 w$ z7 [ W" z$ y
0 W7 F7 [8 v9 b* Q! e6 @9 M* F5 _
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
' X3 w% b' u* O9 q9 o/ e
' Y3 a: ~& W! t+ v# T L
9 h" @' ~7 ^6 B3 n/ e
) H, ^( l" m0 r4 }
而起始及边界条件同上。
: l! n& V- c* n- C$ P/ ~
6 |' W) @ p, i" d, t5 g
在获得浓度分布后，即可以 Fick’s law
& L, d% }- ~ q7 f
5 ], _2 g# O) X0 E
0 R" ]5 ?( [8 \1 z2 H3 o8 ?
- e; p0 T+ o( @7 P9 X: h( u
计算流通量。
" r0 M( I* A2 r3 y. X
6 K* @5 ]) U2 j4 s# w
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
( O6 k) z8 [4 H7 {% k8 D8 _/ d& F
' I3 ]. K, \0 |& o. Y* V
3 L+ @& y% C% L4 D4 H+ P
( }/ O* F* J5 v+ K7 ~1 x
利用以上的处理结果，可编写 MATLAB 参考程序如下：
0 @6 I' p2 l k$ K
# d; x9 o* O7 }0 q8 r' L& Q
function ex20_3_2
; O7 R. V. Q3 s
%*****************************
$ T6 t) [$ m% M
% 扩散系统之浓度分布
3 e6 t8 O f" d* F4 z3 R0 _5 ^% Z
%*****************************
) e, Y! z& [& B
clear
/ H, M$ I" Q# T% U# ?+ |% E
clc
6 q f" O% U: c
global DAB k CA0
0 }5 ~8 ^7 Q* H; p
%******************************
% @( `+ f& C. t
% 给定数据
5 q5 D* e$ z( H, r" C$ ]$ {- ~
%******************************
E0 h4 n# S( L" L
CA0=0.01;
, V$ E" m3 T* ~4 |0 B
L=0.1;
W! p5 M5 d# H* c1 h P- L- [
DAB=2e-9;
) P* z5 h4 Q4 L2 h
k=2e-7;
, w2 q6 t/ y& s- f8 x
h=10*24*3600;
% }9 L$ K% o9 X( q
%*******************************
0 m3 e2 j( `! ?3 d
% 取点
" K0 E: R0 z6 F6 Y
%*******************************
! d% h( h+ E) i
t=linspace(0,h,100);
1 x7 F; _ Q5 }* q, i3 G/ T% m
z=linspace(0,L,10);
" v) g* J' {7 N2 [; W
%*******************************
) D) w5 a0 R3 p* j
% case (a)
7 q3 \% p4 Y2 @7 M4 j
%*******************************
$ ^' C' T3 L- o0 x5 ^7 U1 }
m=0;
! o; K; z' u! @/ m+ R% N1 G0 n6 A
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
0 Y# E* D1 _0 p( U
CA=sol(:,:,1);
4 P; k x' M5 d$ }: a9 X
for i=1:length(t)
2 h; O3 r9 J6 j% Q5 h9 V6 l. D
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
0 ]: e$ y" e, J: @7 B% j! l/ w
NAz(i)=-dCAdz_i*DAB;
. \; ?' B$ ?, p
end
5 s; L* x1 c0 Q9 I1 _" f
figure(1)
, N; Z9 @7 S# ^- T
subplot(211)
6 n! Z4 k" W7 L
surf(z,t/(24*3600),CA)
+ U2 |! C( {9 {: Y& z" L/ [, X
title('case (a)')
- r/ N% k9 J. R# c; O
xlabel('length (m)')
. n+ x, p1 m0 `2 x$ Q: N9 I" j
ylabel('time (day)')
) `$ p" G. v/ n: h6 p0 a
zlabel('conc. (mol/m^3)')
( S' y5 H. f2 h5 p2 z b2 q
subplot(212)
- ]( ?6 W+ x/ l, {
plot(t/(24*3600),NAz'*24*3600)
& ?% B' T! I' }9 Q) _
xlabel('time (day)')
' n% @) _# L1 Y1 P
ylabel('flux (mol/m^2.day)')
- s) o- V' `# K
%************************************
5 n$ G# B7 s! i* }% _* w2 g
% case (b)
( n( ?- Y1 S% [/ U" Z/ J4 a
%************************************
3 g) f; b$ h: o5 L: \1 U
m=0;
/ P+ |# ~* H/ G: y4 U" x* q J
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
; c, `" Y7 r) O X8 _+ b
CA=sol(:,:,1);
1 `5 f4 U& p/ C1 a; K ~7 o) x
for i=1:length(t)
0 F( D+ N* q" ~% \" X! R' `
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
+ ` v0 a: O" W: v
NAz(i)=-dCAdz_i*DAB;
2 P, v3 t" j5 J1 }$ ^% _: \1 q
end
7 z% C: [# j& _2 m! y
%
; |8 v3 W/ W* A& f% a1 w
figure(2)
* W- f! p- A C* c+ n$ n. f
subplot(211)
! @/ p5 l( X2 B- ]
surf(z,t/(24*3600),CA)
, S/ x# ]8 B5 G
title('case (b)')
* U1 V/ q2 F0 D: H. @$ H- X a0 Z
xlabel('length (m)')
; \+ k; G+ h8 X U, j! _
ylabel('time (day)')
8 q2 Z/ e% J: d- I, V, A4 N" h ^$ W
zlabel('conc. (mol/m^3)')
8 N# j9 [0 N+ W* O
subplot(212)
4 l) P% |8 S4 w& @4 L' B
plot(t/(24*3600),NAz'*24*3600)
% u _' z" Q' h* |% u# H! P
xlabel('time (day)')
" [0 N6 G* Q/ M. ]- C/ V Y0 t
ylabel('flux (mol/m^2.day)')
4 S U: r- j5 _
%********************************************
& W1 ^6 Z1 a+ T g/ b; \% S7 J
% PDE 函数
4 L' h$ a2 e# s5 }2 k. W- q/ s5 A4 T
%********************************************
) Q, L9 E0 \6 q6 r w) v
% case (a)
3 ?6 z0 E, ^+ M. J5 l
%********************************************
5 t9 s& C* |0 A1 V1 k: A
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
9 m! ] F1 T& J- W v1 q, k
global DAB k CA0
1 q+ G+ X3 e d J/ l
c=1;
1 N# z) M; s9 o+ }( r& }) r6 P1 Q2 j
f=DAB*dCAdz;
{6 b2 I8 o5 X. \/ R! S- D& }. F
s=0;
. [6 ~, A: q+ ?0 U) @; f/ u1 D
%*********************************************
6 M' n2 X4 R6 a6 a
% case (a)
5 K) h6 d. ~4 C) b9 k) Y; m
%*********************************************
2 E! i1 H' r& h9 ]" J6 t
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
6 u2 k1 @5 Q0 `- P( |3 a
global DAB k CA0
; C& b$ E$ m- w, N7 |& {0 U
c=1;
1 ?+ u& w! Y1 ~: p' m8 N2 z1 c( {$ y
f=DAB*dCAdz;
, N/ Y% O1 @( e* i
s=k*CA;
& s& E: W. ]4 H s# W( l7 d
%**********************************************
( g7 C! `' h7 Q( [0 ^8 n8 I
% 初始条件函数
- |" a1 H U& p9 c/ m7 f6 z
%**********************************************
" D+ g5 G8 E) X# T* |% U+ F
function CA_i=ex20_3_2ic(z)
8 Q: N( A" l6 r w3 Q5 \
CA_i=0;
2 d }) K4 o7 B. l/ @
%************************************************
% Y$ x4 t# Q4 q+ Z
% 边界条件函数
+ S0 ~3 W$ @) _/ P5 y# W: j8 Z
%************************************************
1 H* I6 w' U0 v" O- |
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
1 h. A4 |1 I( q) [7 ~7 v5 N
global DAB k CA0
s9 g& d1 u) a) M/ ]
pl=CAl-CA0;
( H8 P" W& O2 y4 C8 n
ql=0;
! y @) d: K# u2 `
pr=0;
+ R& P+ |) Y$ Q* A
qr=1/DAB;
# }7 l! I. t2 l* [; f0 O
# Z9 ^" X+ R, T: G) _0 c% Q, d
————————————————
) Z' h$ B% B- E! ]* W! [ s
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
9 y3 W( F+ B7 V5 G
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
+ ?; X/ }# v6 B* v( A; J; S/ [
+ S% v% h" d I
. p2 s& o) v: ^% L
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
* m2 @! @; `; ^- O m) S* N+ R0 M
1 w$ z7 [ W" z$ y
0 R" ]5 ?( [8 \1 z2 H3 o8 ?
3 L+ @& y% C% L4 D4 H+ P
,0);0 ]: e$ y" e, J: @7 B% j! l/ w,0);+ ` v0 a: O" W: v