batchdistill.m

Matlab在化学工程中的应用

function BatchDistill
% 恒定滞液量的双组分间歇精馏的回流比控制
%
% 假设：
% 1. 沸腾上升的蒸汽流量（V）恒定；
% 2. 下降的液体摩尔流量（L）恒定； 
% 3. 塔板上的滞液量（M）恒定；
% 4. 再沸器和塔顶冷凝器中的滞液量恒定。
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/01/24 $

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;
tf = 300;

% 初始值
MB = charge;
dist = 0;                       % 总的馏出液量
x0 = [zeros(1,8) xB MB dist];   % x0 x1 x2 x3 x4 x5 x6 x7 xB MB dist
tspan = [0 tf];
[t,x] = ode45(@DistilEqs,tspan,x0)
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),'b-',t,x(:,8),'k--')
xlabel('时间 (h)')
ylabel('液相摩尔分率')
legend('x_0','x_1','x_2','x_3','x_4','x_5','x_6','x_7')
x0 = x(:,1);
Flag = ones(size(t,1),1);
Flag(find(t<tc)) = 0;    % t>=tc时，开始进入控制状态（flag = 1）
Rc = KC*(x0set-x0).*Flag;
figure
plot(t,Rc)
xlabel('time (h)')
ylabel('Rc')

function dxdt = DistilEqs(t,x)
global alpha M0 M V Rinit x0set KC Flag tc

% x(1...11)表示[x0 x1 x2 x3 x4 x5 x6 x7 xB MB dist]
xB = x(9);
MB = x(10);
dist = x(11);

% 当塔顶馏分的组成超过设定值x0set时，进入控制状态，记下此时的时刻tc
if x(1)>=x0set
    Flag = 1;
    tc = t;
end

%  控制方程
Rc = KC*(x0set-x(1))*Flag;
if Rc<0    % 限制Rc> = 0.0
    Rc = 0;
end

R = Rc+Rinit*(1-Flag);
D = V/(R+1);
L = V-D;

% 汽液平衡关系
i=2:9
y(i) = alpha*x(i)./(1+(alpha-1)*x(i));

% 塔顶冷凝器
dxdt(1) = (V*y(2)-(L+D)*x(1))/M0;      % d/dt(x0)

% 第1~7块塔板的易挥发组分质量平衡方程
for i=2:8
    dxdt(i) = (L*(x(i-1)-x(i))+V*(y(i+1)-y(i)))/M; 
end

% 塔釜
dxdt(9) = (L*x(8)-V*y(9))/MB;       % d/dt(xB)   
dxdt(10) = L - V;                   % d/dt(MB)   

% 馏出液总量
dxdt(11) = D;                       % d/dt(dist)

dxdt = dxdt';