Matlab在求解扩散系统之浓度分布中的应用

读书报告
Matlab在求解扩散系统之浓度分布
中的应用
池雨
一、问题的提出
管中储放静止液体 B ，高度为L=10 ㎝，放置于充满A 气体的环境中。

假设与B 液体接触面之浓度为C A0=0.01mol/m 3，且此浓度不随时间改变而改变，即在操作时间内(h=10天)维持定值。

气体A 在液体B 中之扩散系数为D AB =2×10−9m 2/s 。

试决定A 与B 不发生反应；情况下，气体A 溶于液体B 中之流通量(flux)。

参考如图所示的装置。

二、知识背景
Fick 第一定律：实验表明，在稳态扩散的条件下，单位时间内通过垂直于扩散方向的单位截面面积的扩散物质的通量与浓度梯度成
正比。

数学表达式为： Fick 第二定律：根据质量平衡关系即在微小体积中积存的物质
（留入的物质量）J 1—（留出的物质量）J 2得出 因此Fick 第二定律的数学表达为：
三、问题求解
根据题意不同时间t 和距界面厚度不同处x 的浓度C=f(z,t); 因气体A 与液体 B 不发生反应，故其扩散现象的质量平衡方程根据Fick
dx
dc
D
J -=dx J J dt dc 2
1-=2
2
x
c
D t c ∂∂=∂∂
第二定律。

依题意，其初始及边界条件为： I.C. C A0(Z,0)=0, Z>0 B.C. C A (0,t)=C A0, t ≥0 ;
在获得浓度分布后即可应用Fick 第一定律求得流通量即 四、matlab 程序设计
偏微分方程（Partial Differential Equation ，简称PDE ）就是涉及到两个自变量以上的微分方程。

在化学工程领域，为了更好的进行过程设计、优化和控制，经常需要了解化工设备（如反应器）中的温度、浓度和速度在不同空间上的分布以及随时间的动态变化规律，因此涉及到许多偏微分方程的问题。

Matlab 函数pdepe
（）可用于求解偏微分方程，模型为：
用以解含上述初始值及边界值条件的偏微分方程MATLAB 命令 pdepe 的用法如下：
若要获得特定位置及时间下的解，可配合以 pdeval 命令。

使用格式如下：
,0≥=∂∂=t t
C L
Z A 0
)(z
=∂∂-=z A AB
t AZ C D N
第一步：Fick 化为标准形式后即
第二步：编写偏微分方程的系数向量函数： function [c,f,s]=ex20_3_2pdefun(z,t,CA,dCAdt) c=1/D AB ; f=dCAdt; s=0;
第三步：编写起始条件：function CA_i=ex20_3_2ic(z) CA_i=0;
第四步：编写边界条件：
function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t) global DAB k CA0 pl=CAl-CA0; ql=0; pr=0; qr=1/DAB;
第五步：取点。

例如 t=linspace(0,h,100); z=linspace(0,L,10);
第六步 利用 pdepe 求解。

m=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
),,,(z ),,,(1
),
,,(c ==∂∂∂∂=
∂∂=∂∂m z
c c t z s c z c c t z f D z c c t z A
A A
A A AB
A A
Sol(j,k,i)第一维代表时间t，第二维代表空间位置z，第三维代表解向量u的第i个元素，即ui=sol（：，：，i）
第七步显示结果。

u=sol(:,:,1);
surf(x,t,u)
title('标题')
xlabel('位置')
ylabel(''conc. (mol/m^3)'' )
zlabel('u')
若要获得特定位置及时间下的解，可配合以pdeval 命令。

使用格式如下
本题的程序编写为：
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0);
NAz(i)=-dCAdz_i*DAB;
End
附录（matlab程序）：
function ex20_3_2
% 扩散系统之浓度分布
clear
clc
global DAB k CA0
% 给定数据
a=input('请输入要求通量处厚度a= ')
CA0=0.01;
L=0.1;
DAB=2e-9;
h=10*24*3600;
% 取点
t=linspace(0,h,100);
z=linspace(0,L,10);
m=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t); CA=sol(:,:,1);
for i=1:length(t)
[CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),a);
NAz(i)=-dCAdz_i*DAB;
end
subplot(211)
surf(z,t/(24*3600),CA)
title('case (a)')
xlabel('length (m)')
ylabel('time (day)')
zlabel('conc. (mol/m^3)')
subplot(212)
plot(t/(24*3600),NAz'*24*3600)
xlabel('time (day)')
ylabel('flux (mol/m^2.day)')
%************************************************ % PDE 函数
%************************************************ function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
c=1/DAB;
f=dCAdz;
s=0;
%************************************************ % 初始条件函数
%************************************************ function CA_i=ex20_3_2ic(z)
CA_i=0;
%************************************************ % 边界条件函数
%************************************************ function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
pl=CAl-CA0;
ql=0;
pr=0;
qr=1/DAB;。