数学建模社区-数学中国
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:29
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
# Q* n( v; Q. R6 e1 N* p J
题意解析：
" t: b$ q( D+ \
! Y# p. X& ?! \2 p4 O7 T! t7 Y
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
8 |3 P! o, @" W: _5 O
! h( G6 C$ ?: v( A) t& b
A9 g% c- w9 s2 i1 I4 ]6 Z
6 A7 k. S5 G: E8 G8 U/ p
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
2 H8 _) Q1 {% Y. h2 O- h# A
$ d. P% O9 O- A2 t- d+ j |, v
& F1 ^, J0 T7 Y3 Z4 X; w
J6 ?7 K' m2 B
而起始及边界条件同上。
) c: X" T e6 y% r- X- H! V
$ r) V' O' t8 Q" g8 X6 T& A
在获得浓度分布后，即可以 Fick’s law
" W. u1 a4 ]) D! x
% Y0 o' p6 i1 Y( z+ S( P5 ?% m8 D
( l3 S' B* b Q# {# R6 I' A$ V/ G0 D
9 d9 H* f: G _$ C* P4 Y8 t6 @
计算流通量。
6 b2 A' ?) X! {
% D1 N M+ G' t, @) K; j8 G
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
$ u) f5 e5 N/ q2 D% p" R" k
9 Q t% K d1 \7 R( B# u
9 A! O( t N, b1 I
: \9 r# G8 h: i5 G' A' F5 C" C
利用以上的处理结果，可编写 MATLAB 参考程序如下：
( A' n. O: I. j7 y
7 K% }5 V! _* y. A
function ex20_3_2
) [4 R3 w% X" L Y, n& ?3 u
%*****************************
. b' a* g" t0 ^6 F* e; k
% 扩散系统之浓度分布
0 _- _9 J. G, \, x
%*****************************
' B$ B5 ]% G8 V! v. r0 _
clear
4 J4 {9 C( a; q% K N
clc
, Y. `: R* V" [ n1 }2 w. g2 A
global DAB k CA0
% f c. a' ?) c" U! w4 d
%******************************
: s7 d, X% Y& Y3 _4 h1 I& Z
% 给定数据
/ s, A" {" W2 n, [
%******************************
) l, I! u5 y! w
CA0=0.01;
% o% P q( w# K9 L; b
L=0.1;
- V1 Q2 u" J7 t. X
DAB=2e-9;
* ~" `: T- [8 A) P# l" y! X
k=2e-7;
) |% I( D5 ~3 j0 f: ^6 h
h=10*24*3600;
# M; ~2 x d& j3 E) m
%*******************************
5 g/ ^/ i2 k$ v8 G, I1 T3 J, C- R
% 取点
1 A+ [& p! J6 z) \& h
%*******************************
" M- @, z% S. m2 U+ K! s7 `
t=linspace(0,h,100);
5 e I" q. Q$ Z3 a7 I
z=linspace(0,L,10);
5 `2 G* [! q }( } N- ~
%*******************************
; D) B/ o# W/ `4 t
% case (a)
1 ~% ?) g: f0 I$ C( g
%*******************************
3 a: h- O: i1 F
m=0;
* P3 O$ ^; U6 E) N! J, z Q2 C
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
`2 e, R. ~/ a" K0 ]5 n5 ^
CA=sol(:,:,1);
; C: L. k; H$ E I
for i=1:length(t)
$ j/ ~+ j" l2 y4 E+ A' @
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
, u4 a5 j6 y$ T9 W
NAz(i)=-dCAdz_i*DAB;
) ? u0 P' {4 m% h7 L( L
end
3 A/ {0 s& U3 F
figure(1)
4 e# h4 x4 Z, ~' I" c
subplot(211)
6 a/ ~" B* _$ n0 S7 \
surf(z,t/(24*3600),CA)
' l) V& `) U% k
title('case (a)')
: d9 t# Z2 X8 ?" B/ \
xlabel('length (m)')
0 F) P# I: ]( W1 i- h, Y0 J* T, _
ylabel('time (day)')
. ^5 m; Z& Y9 [1 G+ ~! ]
zlabel('conc. (mol/m^3)')
" {- f4 s/ N3 g" E
subplot(212)
- J9 u$ {7 p# X+ n9 P* U! H" Q
plot(t/(24*3600),NAz'*24*3600)
+ M) @4 g5 |, r; h8 O) U T
xlabel('time (day)')
! I3 z/ s6 @7 A. k
ylabel('flux (mol/m^2.day)')
3 @6 }4 J' o" N0 v) U: e8 R6 O
%************************************
) T+ G6 e U7 l' j8 m% o: @/ y* V
% case (b)
. W, W% w6 u/ [+ m
%************************************
. [( L# O; C K( V
m=0;
3 W! F" O7 Q+ @% f' m" l9 i
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
h4 f+ }; ]) ^. g6 q0 D1 \3 e
CA=sol(:,:,1);
1 e* @ {3 R% q
for i=1:length(t)
) A: X! k# m! y$ U. H* V0 M7 }% L
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
; R) n+ x( f7 N6 O6 }2 n
NAz(i)=-dCAdz_i*DAB;
- C. j! E* c" K4 ?& E2 ?% w
end
) n8 X3 \3 {9 ]' @8 C3 n
%
: b% A4 Y% d! w8 `; |9 A2 k" T
figure(2)
3 G' x8 Z$ M4 `! E% y3 }
subplot(211)
. c/ m" \$ t' e, C, q- |
surf(z,t/(24*3600),CA)
' y- X0 @3 m! U. t: v* r
title('case (b)')
. [3 ?8 S- @9 F0 o, @' K7 s
xlabel('length (m)')
/ B. K M# h0 W" [ J8 P
ylabel('time (day)')
: j* h1 i" ]4 O8 H
zlabel('conc. (mol/m^3)')
5 J( g9 L: d P+ ]( Y$ q/ C: N8 L
subplot(212)
3 k+ O+ A; z) e3 h$ ^! x
plot(t/(24*3600),NAz'*24*3600)
3 Z) P! k, O* K3 h- `
xlabel('time (day)')
* N3 S0 n" G% P* e B
ylabel('flux (mol/m^2.day)')
9 B( Z4 k! N7 T# h$ Z) P6 D
%********************************************
& A* t# w+ b8 @6 c$ k
% PDE 函数
5 f) L& d w( e5 ?
%********************************************
( Q, G! _6 D2 e1 c
% case (a)
( _7 x$ z# r5 j3 i+ u
%********************************************
2 t4 O" V! Q5 P( Q' w J3 B
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
# x1 L9 A* b: A! B2 z k5 D
global DAB k CA0
8 j- p3 b' _+ K8 h- X, c6 f
c=1;
( [ X# R5 @: L1 @; D
f=DAB*dCAdz;
% G/ h9 p+ t# A; i2 V
s=0;
( `# p/ D2 N$ G5 S1 E) o
%*********************************************
- S" x5 g" h6 E x, t' L
% case (a)
' L% h2 \+ r, j4 P/ b# C1 {+ l
%*********************************************
8 y4 Y' L9 _) [- g9 E
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
$ [- P1 M8 a+ U9 J7 ~
global DAB k CA0
1 W- B! f2 W, D9 {
c=1;
' s% L1 z% Q% j0 a3 O) \6 u, t
f=DAB*dCAdz;
$ N* _8 ~% ~9 R! M6 }$ K
s=k*CA;
& [0 _! g% S& @; Q- [9 Z
%**********************************************
3 J7 ?8 l H: Z3 {- f
% 初始条件函数
1 B& f( r$ b0 J' I$ _5 T
%**********************************************
, H9 o/ M I! \0 i) D" U9 ?% U
function CA_i=ex20_3_2ic(z)
1 Y! T, M, |+ O. [2 w
CA_i=0;
8 V4 o l8 ~; |9 ?* d3 F9 @
%************************************************
! l7 k- q5 b1 d/ f5 s* Z
% 边界条件函数
; N, ?5 l/ Q" ?. _/ ^ I
%************************************************
0 N# x8 o* o c( _/ k& {# [
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
: Y7 u0 @. D( }5 C4 o
global DAB k CA0
' _ k8 D7 z n9 ~! N" [
pl=CAl-CA0;
$ `4 v% z; F, G, ~
ql=0;
1 I5 ?5 G) h" c: P
pr=0;
' S0 d6 S6 y# S X, g6 Z
qr=1/DAB;
3 q% y# V4 S( ?2 ~! \1 b
& H! ^# C; Y2 b' q5 {5 X
————————————————
7 I2 u- v; {4 N
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
8 S5 F5 [# ]/ k% q$ W/ R
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
) j. q" }6 u. P3 j
& N w- ?6 O4 d8 J
6 j* x+ e* y8 b) E3 l
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
A9 g% c- w9 s2 i1 I4 ]6 Z
( l3 S' B* b Q# {# R6 I' A$ V/ G0 D
9 A! O( t N, b1 I
,0);,0);