偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布
https://img-blog.csdnimg.cn/20190430161546862.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70题意解析：
(a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下：
https://img-blog.csdnimg.cn/2019043016165643.png
(b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为
https://img-blog.csdnimg.cn/20190430161717296.png
而起始及边界条件同上。
在获得浓度分布后，即可以 Fick’s law
https://img-blog.csdnimg.cn/20190430161810270.png
计算流通量。
MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下
https://img-blog.csdnimg.cn/20190430161834231.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
利用以上的处理结果，可编写 MATLAB 参考程序如下：
function ex20_3_2
%*****************************
% 扩散系统之浓度分布
%*****************************
clear
clc
global DAB k CA0
%******************************
% 给定数据
%******************************
CA0=0.01;
L=0.1;
DAB=2e-9;
k=2e-7;
h=10*24*3600;
%*******************************
% 取点
%*******************************
t=linspace(0,h,100);
z=linspace(0,L,10);
%*******************************
% case (a)
%*******************************
m=0;
sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
for i=1:length(t)
=pdeval(m,z,CA(i,:),0);
NAz(i)=-dCAdz_i*DAB;
end
figure(1)
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)')
%************************************
% case (b)
%************************************
m=0;
sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t);
CA=sol(:,:,1);
for i=1:length(t)
=pdeval(m,z,CA(i,:),0);
NAz(i)=-dCAdz_i*DAB;
end
%
figure(2)
subplot(211)
surf(z,t/(24*3600),CA)
title('case (b)')
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 函数
%********************************************
% case (a)
%********************************************
function =ex20_3_2pdefuna(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
s=0;
%*********************************************
% case (a)
%*********************************************
function =ex20_3_2pdefunb(z,t,CA,dCAdz)
global DAB k CA0
c=1;
f=DAB*dCAdz;
s=k*CA;
%**********************************************
% 初始条件函数
%**********************************************
function CA_i=ex20_3_2ic(z)
CA_i=0;
%************************************************
% 边界条件函数
%************************************************
function =ex20_3_2bc(zl,CAl,zr,CAr,t)
global DAB k CA0
pl=CAl-CA0;
ql=0;
pr=0;
qr=1/DAB;
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711694
页:
[1]