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

读书报告

Matlab在求解扩散系统之浓度分布

中的应用

池雨

一、问题的提出

管中储放静止液体 B，高度为L=10 ㎝，放置于充满A气体的环境中。假设与B液体接触面之浓度为CA0=0.01mol/m3，且此浓度不随时间改变而改变，即在操作时间内(h=10天)维持定值。气体A在液体B中之扩散系数为DAB=2×10−9m2/s。试决定A与B不发生反应；情况下，气体A溶于液体B 中之流通量(flux)。 参考如图所示的装置。

二、知识背景

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

dcdx Fick第二定律：根据质量平衡关系即在微小体积中积存的物质

dcJ1J2（留入的物质量）J1—（留出的物质量）J2得出 dtdx2cc因此Fick第二定律的数学表达为： D2tx

正比。数学表达式为：J D三、问题求解

根据题意不同时间t和距界面厚度不同处x的浓度C=f(z,t); 因气体A与液体 B 不发生反应，故其扩散现象的质量平衡方程根据Fick

第二定律。

依题意，其初始及边界条件为： I.C. CA0(Z,0)=0, Z>0 B.C. CA(0,t)=CA0, t≥0 ;

CAtZL0,t0CAzz0在获得浓度分布后即可应用Fick第一定律求得流通量即 NAZ(t)DAB四、matlab程序设计

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

Matlab函数pdepe（）可用于求解偏微分方程，模型为：

用以解含上述初始值及边界值条件的偏微分方程MATLAB 命令 pdepe 的用法如下：

若要获得特定位置及时间下的解，可配合以 pdeval 命令。使用格式如下：

第一步：Fick化为标准形式后即

c(z,t,cA,cA1)zDABcAcAf(z,t,cA,)zzcAs(z,t,cA,)0zm0第二步：编写偏微分方程的系数向量函数： function [c,f,s]=ex20_3_2pdefun(z,t,CA,dCAdt) c=1/DAB ; 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);

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;