连续多组分三元精馏塔的模拟计算 下载本文

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=tc时,开始进入控制状态(flag = 1) i=2:9

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