数学建模社区-数学中国
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
[打印本页]
作者:
浅夏110
时间:
2020-6-10 10:29
标题:
偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
) W3 y/ A% A/ [+ V* ]- k- O2 q
题意解析：
& u( v, ~7 T9 G6 ?" I% O, q: A
0 k5 @" Y7 b0 r: D2 J+ T# c7 U; Q
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
' a- V9 J, r8 t& Z! v/ E7 i7 v/ L
& G" V" U) r9 L" D
+ Z" p7 h* F+ |1 ]/ S
1 ?) m) z* m) t: F3 n
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
+ [% U$ N) S; _* d! @7 ?4 ?% u
' t* G- M( n3 }
5 Z# Z2 j1 H) }1 E+ p* x) t/ u
% @4 P7 k: X# L4 V
而起始及边界条件同上。
+ C5 H) b: t; ]! B
5 w p& a1 x) U5 u3 a3 S
在获得浓度分布后，即可以 Fick’s law
; C& ~4 C# o" A' x! ^ h
$ P" `5 {: R5 j2 x$ y
4 k) f7 Q* a4 T: \+ S
' m- y, E$ Z$ K6 Z% V$ Q$ I8 s
计算流通量。
1 }0 ^1 @2 p9 ^9 f* h) Y7 w9 O
. w& w9 }: c7 T3 u; v* g0 `! y
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
3 {% `6 O4 p+ D- [ h% A
$ C# B! I2 r; j `1 F
8 z, ?+ `7 ^! f3 f; j. \! r. l
/ Q5 w0 H1 I* \% t' ^
利用以上的处理结果，可编写 MATLAB 参考程序如下：
# p( \ ~1 F7 i. C, N
9 E' S& C* [; A# h( b8 _) M
function ex20_3_2
, V4 c! ]# r# V" c
%*****************************
) E7 c$ f3 a; u2 X
% 扩散系统之浓度分布
, ?$ g8 \ d! g( }& |: {% }
%*****************************
, G# X5 q `' P3 Y5 z
clear
0 s9 f. a/ s3 {) E' ~; P* M; C( q2 C
clc
" {% z, V c7 H' m
global DAB k CA0
5 r$ {# f' R" d; h) I
%******************************
. K9 x2 |8 }# r) F+ f+ @' ^
% 给定数据
% i1 W _6 e' }, [
%******************************
8 t4 s6 A/ z- q2 [5 r9 T
CA0=0.01;
! q) [1 d2 v% N/ e$ u
L=0.1;
; n9 L4 ^' ?% J4 E
DAB=2e-9;
! g8 z6 a% ^7 } F7 b, l. M
k=2e-7;
' z/ n8 O5 r- R7 R. J3 `4 S
h=10*24*3600;
$ [3 c; R! [) O8 {
%*******************************
% s& v5 O6 m* A
% 取点
/ k4 E: P. G# [
%*******************************
$ o, y2 j$ z6 E" z. B% \) z
t=linspace(0,h,100);
4 f# o0 x/ U/ n+ X5 I3 s% A( b8 B
z=linspace(0,L,10);
0 t% v& J6 [, x8 Z
%*******************************
3 I) `% [1 x1 @7 f% a1 O
% case (a)
( H `% }$ B" Z7 l6 e
%*******************************
& d) L+ n4 n( q3 P: {
m=0;
4 N) V F# t& B2 U& l; M9 R
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
) M. ^4 [" M7 H8 w! k8 E: r
CA=sol(:,:,1);
- s' H8 n& F# k1 r
for i=1:length(t)
1 k( o/ s) [7 n/ {+ V b
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
$ W# a7 ]* H8 V
NAz(i)=-dCAdz_i*DAB;
. V: G: y7 ]" T! u2 x* D2 ^+ ^. s" y& A
end
4 U4 x. g5 n5 E/ |7 ]- z: [/ y7 O
figure(1)
5 G$ `- G$ I d; O ~1 g
subplot(211)
0 K5 q5 U) J; K: Y& U, I8 V
surf(z,t/(24*3600),CA)
+ |. {" w8 ?' V) k& V5 K
title('case (a)')
, Q1 X- [: {+ i) D4 \: i
xlabel('length (m)')
0 T3 |" i2 J/ r0 P* G0 I
ylabel('time (day)')
. u6 T9 n+ c; A. k7 N
zlabel('conc. (mol/m^3)')
1 f6 i, x9 ?' t s1 d' R: G% _- l
subplot(212)
# c! @$ p- i+ E: k5 y
plot(t/(24*3600),NAz'*24*3600)
" d6 m5 s6 h, y4 q/ m
xlabel('time (day)')
' G- Q7 B' q+ P
ylabel('flux (mol/m^2.day)')
9 B/ E, H3 G5 j( t! X
%************************************
( u' v. n: x" M& o5 W% G% K
% case (b)
, l9 E# Y; ?6 j' @1 X7 L6 @
%************************************
! g" ^( f8 h# g3 ~/ u2 p- J) `
m=0;
: U* t4 S) E# j5 t3 m( h3 W% @6 u
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
0 ?' ^' i) c+ c+ t, Y- @; B; O
CA=sol(:,:,1);
6 q6 h# R O- u# u# ^& s, I
for i=1:length(t)
' E" L7 d" ` i, ?; c+ [$ y8 ^
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,
,0);
% b/ ~6 V2 q% W5 t9 ?: z; f
NAz(i)=-dCAdz_i*DAB;
' N! c6 m! P& G4 c) `) U% c
end
. Q, @* R/ @) H! X8 C7 \
%
, C8 @9 Q! V: L( k6 R& E: i1 _5 L
figure(2)
) W+ j; e8 e% S6 B9 P! q$ m7 N
subplot(211)
4 p& Q8 X4 M, S* i! p1 S
surf(z,t/(24*3600),CA)
# r# g; z5 j3 K8 F* o
title('case (b)')
9 Z8 e" ?3 e/ L+ v! ?2 S
xlabel('length (m)')
, ?; z# A/ n; d3 K
ylabel('time (day)')
9 ], g- A; k) d* U
zlabel('conc. (mol/m^3)')
; H. q* \1 ?0 ~) z8 T! o q. `; [
subplot(212)
8 @3 y i p8 m5 V
plot(t/(24*3600),NAz'*24*3600)
' `5 x* d/ z Q2 Y) t
xlabel('time (day)')
8 N- C' y/ ]1 ?
ylabel('flux (mol/m^2.day)')
. h% b# A" N% \6 J& C
%********************************************
2 ^) d7 S0 w5 v& t5 X+ _# W9 k" B
% PDE 函数
+ [7 v9 l$ D/ h
%********************************************
9 W3 i' _9 \# {3 N
% case (a)
/ R r& i ?/ T, A) a
%********************************************
! b9 u: y- M8 [8 U. {
function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
/ L+ x$ D; m! l! R6 Z
global DAB k CA0
: N! J. \& d+ y4 @# _) h) m: A% Z
c=1;
6 F" }' Z0 y1 J: R
f=DAB*dCAdz;
$ i0 M X! u' U9 u9 F
s=0;
% ^* g$ U1 p ^# `# G7 z* i; q# I
%*********************************************
- n$ z/ G. M$ B( O
% case (a)
) V# K: D) s: v8 E
%*********************************************
, C; v6 i; v/ ], C
function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz)
" @" v& i9 q, {% Y1 [
global DAB k CA0
' {2 B' a* @7 I
c=1;
/ z6 C+ j: U) T
f=DAB*dCAdz;
( p: l" i" O9 W* D
s=k*CA;
6 f) t7 V I5 Z- g% l/ g( ?5 N/ x
%**********************************************
: x* V( Z# ?$ A5 ^
% 初始条件函数
& b- C8 V2 P5 r2 I7 a4 z. N* m
%**********************************************
5 X# p4 G8 [, i# G9 u7 G' i9 @) t8 l
function CA_i=ex20_3_2ic(z)
: O. M+ Z! I; x: P
CA_i=0;
4 e: p' {! s3 O- l/ h" h7 Y
%************************************************
' n: W% H& L. p+ S1 Z* e6 Y
% 边界条件函数
" L: t' x0 B& }# i. H
%************************************************
; Z+ {5 J1 \7 @* I7 f9 [
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
% B0 L1 _9 s5 }1 m" \& U- M% t
global DAB k CA0
/ c ?! }2 `6 ~
pl=CAl-CA0;
, _! V& m$ s" L# f* ~) B
ql=0;
% y2 v0 N# E. B! i2 g( @
pr=0;
5 Q B1 y, c) |) t5 S0 n
qr=1/DAB;
9 }1 U( A9 r' w! ]% u5 t0 m
$ D! a* Q+ ?7 B% `8 u, n' E
————————————————
e! z1 \+ X9 `- R( S( M
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
8 A3 V+ |" a1 v& F+ x) S M
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
/ Y$ t8 q7 D" ] V
0 `! s, n$ c; {3 K q
& ]* @9 q+ z$ g" s @
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
+ Z" p7 h* F+ |1 ]/ S
5 Z# Z2 j1 H) }1 E+ p* x) t/ u
,0);$ W# a7 ]* H8 V,0);% b/ ~6 V2 q% W5 t9 ?: z; f