condistill.m

这是一些关于MATLAB的小程序

function ConDistill
% 连续多组分（三元）精馏塔的模拟计算 
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/14 $

clear all
clc

global F z1 z2 z3 R alpha1 alpha2 alpha3 M1 MN M Nt Nf V1 V D L L1 W
F = 40;             % 进料流量, kmol/hr
R = 5;              % 回流比 
z1 = 0.6;           % 苯的进料组成（摩尔分率） 
z2 = 0.25;          % 甲苯的进料组成（摩尔分率）
z3 = 1-z1-z2;       % 苯乙烯的进料组成（摩尔分率）
M1 = 75;            % 塔顶冷凝器中的滞液量(kmol)
M = 10;             % 塔板上的滞液量(kmol)
MN = 150;           % 塔釜中的滞液量(kmol)
q = 1;              % 饱和进料

tf = 35;
dt = 1;

% 相对挥发度 
alpha1 = 2.75;
alpha2 = 1;
alpha3 = 0.4;         

Nt = 10;    % 塔板总数
Nf = 5;     % 进料位置

V1 = 150;   % 从塔釜蒸发上来的蒸汽流量(kmol/hr)

% 精馏段
V = V1-(1-q)*F;
D = V/(R+1);
L = V-D;

% 提馏段
L1 = L+F;
W = L1-V1; 

% 初始化x1和x2---开车时塔内所有板上的x1和x2分别与进料的z1和z2相同
x1 =  z1 * ones(1,Nt);
x2 =  z2 * ones(1,Nt);

[t,y] = ode45(@DistMassBalances,[0:dt:tf],[x1 x2])

% 输出结果
x1 = y(:,1:Nt);         % 苯的液相组成（摩尔分率）
x2 = y(:,Nt+1:2*Nt);    % 甲苯的液相组成（摩尔分率）
x3 = 1-x1-x2;           % 苯乙烯的液相组成（摩尔分率）
plot(t,x1(:,1),'r-',t,x2(:,1),'k--',t,x3(:,1),'b:',t,x1(:,end),'r.-',t,x2(:,end),'k-.',t,x3(:,end),'b.--')
xlabel('Time (h)')
ylabel('x_1_1, x_1_2, x_1_3, x_N_1, x_N_2, x_N_3')
title('塔顶和塔釜产品从进料开始直至稳态的动态浓度曲线')
legend('x_1_1','x_1_2','x_1_3','x_N_1','x_N_2','x_N_3')
% 稳态图
figure
plate = 1:Nt;
plot(plate,x1(end,:),'r.-',plate,x2(end,:),'k.--',plate,x3(end,:),'b.:')
xlabel('塔板')
ylabel('稳态时苯,甲苯,苯乙烯的组成')
title('稳态时精馏塔的浓度曲线')
legend('苯','甲苯','苯乙烯')

% ------------------------------------------------------------------
function dydt = DistMassBalances(t,y)       % 物料平衡方程组 
global F z1 z2 z3 R alpha1 alpha2 alpha3 M1 MN M Nt Nf V1 V D L L1 W
x1 = y(1:Nt);       % 组分1（苯）
x2 = y(Nt+1:2*Nt);  % 组分2（甲苯）
x3 = 1-x1-x2;       % 组分3（苯乙烯）

% 气相平衡
denom = alpha1*x1+alpha2*x2+alpha3*x3;
y1 = alpha1*x1./denom;
y2 = alpha2*x2./denom;

% 对塔顶冷凝器（i = 1）  
i = 1;
dx1dt(i) = (V*y1(i+1)-(L+D)*x1(i))/M1; 
dx2dt(i) = (V*y2(i+1)-(L+D)*x2(i))/M1;

% 精馏段（i = 2~Nf-1）  
for i=2:Nf-1
dx1dt(i) = (L*(x1(i-1)-x1(i))+V*(y1(i+1)-y1(i)))/M;
dx2dt(i) = (L*(x2(i-1)-x2(i))+V*(y2(i+1)-y2(i)))/M;
end

% 进料板（i = Nf）
i = Nf;     
dx1dt(i) = (F*z1+L*x1(i-1)-L1*x1(i)+V1*y1(i+1)-V*y1(i))/M;
dx2dt(i) = (F*z2+L*x2(i-1)-L1*x2(i)+V1*y2(i+1)-V*y2(i))/M;

% 提馏段（Nf+1~Nt-1）
for i=Nf+1:Nt-1
dx1dt(i) = (L1*(x1(i-1)-x1(i))+V1*(y1(i+1)-y1(i)))/M; 
dx2dt(i) = (L1*(x2(i-1)-x2(i))+V1*(y2(i+1)-y2(i)))/M; 
end

% 塔釜（i = Nt）   
i = Nt;
dx1dt(i) = (L1*x1(i-1)-V1*y1(i)-W*x1(i))/MN;
dx2dt(i) = (L1*x2(i-1)-V1*y2(i)-W*x2(i))/MN;

dydt = [dx1dt dx2dt]';