数学建模社区-数学中国
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:29
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
/ i* q. ?$ ~$ X9 {
题意解析：
7 Z/ ?8 b" D+ v
& v$ \& n4 T1 J8 M. U
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
" `$ [ u& j: W" q& j
0 {: P0 ?4 `" u: k- T% x/ p- {
; {4 r7 d8 o4 A% X
; A( B- r# n* J
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
" {, q5 B6 D2 r" T9 A; E% K; s
- i& z: F' F3 l+ G# @
8 U l1 E: ]6 J: m: @1 T
8 t ^7 I$ U+ {8 |8 q$ |
而起始及边界条件同上。
- ?. D' K/ Y0 H" ]
: U' ?, Q8 Z3 G
在获得浓度分布后，即可以 Fick’s law
, o! C; D+ ^% s) s, `
1 l; X) o: D7 L
% d: M+ w: g6 [8 j1 ~
+ L5 i. G6 l$ B! D
计算流通量。
1 F1 G, X) p- ]1 j$ R: m3 d4 K. x
& Z/ p+ h% ]9 n- ?7 {; G- V
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
7 c; P$ m2 |3 X, h3 S& g1 `
4 F0 v$ P: V3 K5 h
- B9 U1 C6 B% ^4 \/ T
# t5 E# b. X) K! ]4 p; b, ^ \
利用以上的处理结果，可编写 MATLAB 参考程序如下：
( D( s& S7 {; d L9 n$ ^0 P% ~% v0 \' b
5 N$ |2 M# n$ k% E
function ex20_3_2
9 e. Y5 K r7 N( T( f2 G( I/ |
%*****************************
+ d8 W; Y1 p/ r% O) n3 r( x L* I
% 扩散系统之浓度分布
2 H! _8 B! S. t4 W; c9 ]
%*****************************
2 h* h4 D% }; m
clear
[! ~! F+ R" O
clc
" N! [( e% y; m* D
global DAB k CA0
1 {, B4 Y3 ^! ~# T" E
%******************************
8 Z5 Z4 f: U4 F" L
% 给定数据
% |) G9 L. a) |. } Z: w9 E; L% H
%******************************
/ h, \6 A5 @- I8 X8 s5 v! I
CA0=0.01;
q& U: x4 |3 F5 i) b
L=0.1;
" T1 x' e( N; v) r6 ~
DAB=2e-9;
$ L. ~: U* B7 k
k=2e-7;
" q7 K$ L- ?" c2 N
h=10*24*3600;
" g" M/ D3 C3 I0 v9 N
%*******************************
9 t, C' B1 F% n) i- `6 }" v1 X4 Y, _
% 取点
3 o# o+ @9 ]& U! x4 Y1 K/ U* O
%*******************************
* I/ i6 _) v" M
t=linspace(0,h,100);
0 o" h9 C9 ]4 L* B( y' j
z=linspace(0,L,10);
5 \8 }. V o& C3 H- J" r
%*******************************
9 P$ L5 A# t9 V- [' P. Y7 S
% case (a)
* F; c0 ?5 p- |. G- \4 E
%*******************************
3 }. @, k& P' m6 |
m=0;
! F9 m1 B0 w4 R& F% R \/ n
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
9 m8 X" K) x7 b# g
CA=sol(:,:,1);
: C( ?0 x3 \ L2 f0 b' Z
for i=1:length(t)
* L% o" Q( D- ^- e# D- C5 X+ S
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
9 K# j9 {2 v [7 h: Z! X, X
NAz(i)=-dCAdz_i*DAB;
: \8 v- Q' ^, p, p: J+ O
end
# O: N1 f& h4 F7 B6 E6 ~
figure(1)
+ W- y! v8 l3 x& t: I
subplot(211)
& n; F7 R, H6 b- Z* |/ T
surf(z,t/(24*3600),CA)
0 |( ]5 q& p. v* I" D
title('case (a)')
5 R& O+ y7 W" X4 v
xlabel('length (m)')
) m3 O `: P! p
ylabel('time (day)')
1 k O! |+ b8 I/ r& e! B% ?* F
zlabel('conc. (mol/m^3)')
( _1 ^# D7 X ?) Y
subplot(212)
2 R/ H3 a: @$ g: S+ Z
plot(t/(24*3600),NAz'*24*3600)
- e, ~4 ?6 y# p9 \1 r( ?
xlabel('time (day)')
! s, S4 \( w- f5 z
ylabel('flux (mol/m^2.day)')
; N! y2 j) O" O
%************************************
4 L8 }0 p4 a0 k$ ^. E
% case (b)
h/ M. ^9 Y8 E, E' h. i
%************************************
) B3 Q4 l$ Y. I+ s' ]- E f" c
m=0;
- R9 ^5 B5 H8 ?- Z- E8 \7 w% V: y
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
2 l+ U; T, b& {; q4 l% m
CA=sol(:,:,1);
, v7 E( Y% Z; Z( U
for i=1:length(t)
2 o* l0 U- r/ S/ B( t- ? `
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
$ [* M( n: J4 {% U0 Y* L
NAz(i)=-dCAdz_i*DAB;
* e6 c0 J5 Y0 \" N
end
& f" j! Y( A5 [+ r d, R; Y
%
! e+ n% N$ G' E+ P1 ?" `5 {
figure(2)
7 y0 x0 S6 r& Q- T
subplot(211)
1 p% O+ D6 {: o/ a" I9 O2 P
surf(z,t/(24*3600),CA)
7 O8 _: W! {2 s8 Y8 _+ J& F) |
title('case (b)')
5 }9 y: J6 v; r2 d
xlabel('length (m)')
/ p# r8 B" u, }6 O$ l3 e
ylabel('time (day)')
4 L* Q2 \) h5 T7 @; V
zlabel('conc. (mol/m^3)')
; f+ N. h: \" }' J5 B% r) u
subplot(212)
A7 O8 i$ K) D$ N1 ], A. W
plot(t/(24*3600),NAz'*24*3600)
+ _: n* {. u2 ^3 N% Q- P. w& s
xlabel('time (day)')
; [# E7 [ ? N1 T# K3 Q
ylabel('flux (mol/m^2.day)')
7 a% b' J# n5 ^ t; ?
%********************************************
2 r7 ~8 c! Y Y @
% PDE 函数
4 a; ^4 N* e8 c) x
%********************************************
; `) B4 t" s5 Q n/ U
% case (a)
. Y7 ]& V+ ]4 c' N
%********************************************
/ }+ [& C) O8 M" T" p
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
2 B" L( m5 q) W" W8 }$ F: N
global DAB k CA0
/ x1 S& Z1 }$ p n& u- f, @
c=1;
8 |: T/ ]% x; a/ X' q$ x* J" F& h; Q" T
f=DAB*dCAdz;
+ C! l% V6 K( m1 h: A& ?. c! I+ X
s=0;
" \+ D; l+ ]* g2 J: Z, J6 ?$ h
%*********************************************
: G5 k3 a" B$ h
% case (a)
7 \7 I; q+ v* y4 g* ~, E
%*********************************************
t) H" j1 R0 _ W
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
) z* [: ^4 s6 K" g2 t: W5 v
global DAB k CA0
$ Z. {( ?5 S5 {6 @* D! O/ G
c=1;
( s$ G1 ~. y+ U! f+ A% k6 Y
f=DAB*dCAdz;
3 T4 J1 a/ }5 N) |
s=k*CA;
, c) t! o( Y; H5 f% @' w
%**********************************************
A( P7 N' S" i% D! N4 t
% 初始条件函数
7 a1 b+ m4 ~/ ]9 }4 H# n7 D; n
%**********************************************
. ]0 T. |0 I0 h: G W
function CA_i=ex20_3_2ic(z)
' Q* [8 B- P2 ~$ B, g( m H
CA_i=0;
9 n& M( a4 _/ m1 C* F
%************************************************
" p& _9 M! d, H; `& ^( X
% 边界条件函数
6 H8 H. H8 S# ~8 }
%************************************************
9 F% w% S6 x9 B1 }
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
. v% s* ]) N. ^& m% o) W
global DAB k CA0
7 N% I$ x+ @4 w1 K
pl=CAl-CA0;
# _9 v. s8 Q$ F$ j S' X
ql=0;
# @, w2 w* l& r2 t0 F9 ~) v
pr=0;
8 E/ i7 S+ }( K- e+ I
qr=1/DAB;
; h( L( `- l2 j% t9 T1 P4 h' D% u
^4 c8 N$ g- a
————————————————
& G$ M: `1 l# w+ A9 i, Y0 {
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
+ r$ L- V8 R/ S8 ?, |; ^
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
: l: ~2 P0 F) D) x- {
/ `3 q0 M0 L$ ^$ m
: J. J1 ^ S, W
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
; {4 r7 d8 o4 A% X
,0);,0);