multicompdiff.m

《实用化工计算机模拟：MATLAB在化学工程中的应用 》这本书光盘里的程序~

function MultiCompDiff
% Simultaneous Multicomponent Diffusion of Gases
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/05/13 $

clear all 
clc
global CT NA NB NC DAB DBC DAC
NA0 = 2.396e-5;
NB0 = -3.363e-4;
NC = 0;
DAB = 1.47e-4;
DBC = 1.245e-4;
DAC = 1.075e-4;
CT = 0.2/(82.057e-3*328);

z0 = 0;
zf = 0.001;
CA0 = 2.229e-4;
CB0 = 0;
CC0 = 7.208e-3;
CAf = 0;
CBf = 2.701e-3;
CCf = 4.730e-3;
zspan = [z0 zf];
y0 = [CA0 CB0 CC0];
flux = fminsearch(@ObjFunc,[NA0 NB0],[],zspan,y0,CAf,CBf)
NA = flux(1);
NB = flux(2);
[z,y] = ode45(@ODEs,zspan,y0,[],NA,NB);
xA = y(:,1)/CT;
xB = y(:,2)/CT;
xC = y(:,3)/CT;

% 输出结果
fprintf('组分A和B从位置1到位置2的摩尔扩散通量分别是：NA = %.2e, NB = %.2e',NA,NB)
plot(z,xA,'-',z,xB,'.-',z,xC,'--')
xlabel('Distance (z),m')
ylabel('Mole fraction')
legend('x_A','x_B','x_C')

% ------------------------------------------------------------------
function f = ObjFunc(flux,zspan,y0,CAf,CBf)
NA = flux(1);
NB = flux(2);
[z,y] = ode45(@ODEs,zspan,y0,[],NA,NB);
f = abs(CAf - y(end,1)) + abs(CBf - y(end,2))

% ------------------------------------------------------------------
function dydz = ODEs(z,y,NA,NB)
global CT NC DAB DBC DAC
CA = y(1);
CB = y(2);
CC = y(3);
xA = CA/CT;
xB = CB/CT;
xC = CC/CT;
dCAdz = (xA*NB-xB*NA)/DAB+(xA*NC-xC*NA)/DAC;
dCBdz = (xB*NA-xA*NB)/DAB+(xB*NC-xC*NB)/DBC;
dCCdz = (xC*NA-xA*NC)/DAC+(xC*NB-xB*NC)/DBC;
dydz = [dCAdz; dCBdz; dCCdz];