1. 连续多组分（三元）精馏塔的模拟计算 function ConDistill 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('苯','甲苯','苯乙烯')
% ------------------------------------------------------------------
1
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]'; 2.恒定滞液量的双组分间歇精馏的回流比控制 function BatchDistill % 假设：
% 1. 沸腾上升的蒸汽流量（V）恒定； % 2. 下降的液体摩尔流量（L）恒定； % 3. 塔板上的滞液量（M）恒定；
% 4. 再沸器和塔顶冷凝器中的滞液量恒定。 clear all
global alpha M0 M V Rinit x0set KC Flag tc
alpha = 3.5; % 相对挥发度
M0 = 100; % 塔顶冷凝器中的摩尔滞液量, kmol M = 5; % 塔板上的摩尔滞液量, kmol V = 10; % 沸腾上升的蒸汽流量, kmol/h Rinit = 1E10; % 全回流条件（全回流-刚开始时） charge = 2500; % 蒸馏塔负荷（处理量）, kmol xB = 0.8; % 蒸馏塔釜进料组成, 摩尔分率
% 全回流时, Flag = 0 Flag = 0;
% 控制参数
x0set = 0.9; % 馏出液组成的设定值 KC = 100; % 控制增益
dt = 0.01;
2
tf = 300; Flag = 1;
tc = t; % 初始值 end MB = charge;
dist = 0; % 总的馏出液量
% 控制方程
x0 = [zeros(1,8) xB MB dist]; % x0 x1 x2 x3 x4 x5 x6 x7 xB MB dist Rc = KC*(x0set-x(1))*Flag; tspan = [0 tf];
if Rc<0 % 限制Rc> = 0.0 [t,x] = ode45(@DistilEqs,tspan,x0)
Rc = 0; plot(t,x(:,1),'r-',t,x(:,2),'k-.',t,x(:,3),'b--',t,x(:,4),'k-',t,x(:,5),'b-.',t,x(:,6),'r--',t,x(:,7),end 'b-',t,x(:,8),'k--')
xlabel('时间 (h)')
R = Rc+Rinit*(1-Flag); ylabel('液相摩尔分率')
D = V/(R+1); legend('x_0','x_1','x_2','x_3','x_4','x_5','x_6','x_7') L = V-D; x0 = x(:,1);
Flag = ones(size(t,1),1);
% 汽液平衡关系 Flag(find(t
Rc = KC*(x0set-x0).*Flag; y(i) = alpha*x(i)./(1+(alpha-1)*x(i)); figure
plot(t,Rc)
% 塔顶冷凝器
xlabel('time (h)') dxdt(1) = (V*y(2)-(L+D)*x(1))/M0; % d/dt(x0) ylabel('Rc')
% 第1~7块塔板的易挥发组分质量平衡方程 function dxdt = DistilEqs(t,x)
for i=2:8
global alpha M0 M V Rinit x0set KC Flag tc dxdt(i) = (L*(x(i-1)-x(i))+V*(y(i+1)-y(i)))/M;
end % x(1...11)表示[x0 x1 x2 x3 x4 x5 x6 x7 xB MB dist] % 塔釜
xB = x(9); dxdt(9) = (L*x(8)-V*y(9))/MB; % d/dt(xB) MB = x(10); dxdt(10) = L - V; % d/dt(MB) dist = x(11); % 馏出液总量
dxdt(11) = D; % d/dt(dist) % 当塔顶馏分的组成超过设定值x0set时，进入控制状态，记下此时的时刻tc dxdt = dxdt'; if x(1)>=x0set
3