甲醇/水二元间歇精馏塔的建模与仿真
间歇精馏过程是过程工业中的重要分离过程之一，由于间歇化工过程适合小批量、多品种和高附加值的精细化学品的生产，近年来其工艺研究和开发得到过程工业界的日益高度重视。
本文以一套甲醇与水二元间歇精馏塔试验装置为研究对象，应用动态物料平衡和相平衡等机理建模法，建立了固定回流比的间歇精馏工艺过程原理模型，以Matlab为平台，编写了模拟间歇精馏塔的一次开车过程，在不同回流比下，对每层塔板、冷凝器以及再沸器中甲醇的浓度的变化，以及第2、4层塔板和冷凝器、再沸器中的温度变化进行了仿真和模拟计算。
甲醇/水二元间歇精馏塔试验装置简介
精馏实验前，原料一次性加入再沸器。再沸器夹套通水蒸汽加热，水蒸气流量用一套流量控制器控制。塔内共设有6块筛板式塔板。塔顶甲醇蒸汽流经冷水冷凝器冷凝，部分回流。工艺流程图如图1所示：
图1甲醇/水二元间歇精馏塔实验装置工艺流程图
动态数学模型的建立
为了建立该甲醇/水二元间歇精馏塔的动态数学模型，做了以下主要假设：
(1)塔内蒸汽流量V(mol/h)和回流量L(mol/h)保持恒定；无挟带和渗漏，冷凝器全冷凝。
(2)每块塔板上的滞留液量恒定(Mi,i=1,2,…,6)；蒸汽滞液量可忽略。
(3)塔板为理想塔板，塔板效率恒定，气液相平衡，气相为理想气体。
(4)绝热操作，塔板和再沸器内气液相完全混合。
基于以上假设，该间歇精馏塔的再沸器、塔顶冷凝器和塔板1到塔板6的动态物料平衡方程如下：
(1)
(2)
(3)
(4)
(5)
(6)
i=2,3,4,5
(7)
式中，Ms为再沸器的滞留液量；Md为塔顶冷凝器的滞留液量；Xi(i=1,2,…,6,d,s)是塔板i(i=1,2,…,6)、再沸器和塔顶冷凝器的液相甲醇浓度(摩尔分数)；Yi(i=1,2,…,6,s)是塔板i(i=1,2,…,6)和再沸器的气相甲醇浓度(摩尔分数)。
假定精馏塔中二元混合物的气相甲醇浓度Y和液相甲醇浓度X的关系，可以近似的以恒定的相对挥发度来表示，即：
， j=1,2,…6,s (8)
静态物料平衡方程和回流比计算式为：
V=L+D (9)
R=L/D (10)
式中，D为塔顶冷凝后收集的甲醇产品流量，mol/h；R为回流比。
经以前的实验研究，得到下式的温度-浓度经验模型，在进行间歇精馏塔仿真控制研究中，下式常用于从已知的浓度计算各塔板的温度。
j=2,4,d,s (11)
利用欧拉法将式(1)—(7)离散化，转变为差分方程依次为(12)—(18)，为采样时间：
(12)
(13)
(14)
(15)
(16)
(17)
j=2,3,4,5
(18)
根据式(12)—(18)以及(8)、(9)、(10)，利用matlab编程，模拟一次间歇精馏从开车到结束的过程，分别对各层塔板以及塔顶冷凝器和塔釜再沸器的甲醇气液相浓度的变化进行仿真。初始条件为：回流比1.8和2.5，塔釜甲醇初始浓度0.45(摩尔分数)，6个塔板的滞液量0.01kmol，冷凝器的初始滞液量为0.05kmol，再沸器的初始滞液量为1.6kmol，塔内蒸汽流量为3.2/3600kmol/s，相对挥发度为3.48，精馏时间为5000s，采样时间为0.5s。
得到各层在r=1.8时从开始到结束的动态浓度变化曲线、第2，4塔板、冷凝器、再沸器温度变化曲线以及在不同回流比(r=1.8和r=2.5)下塔顶冷凝器浓度变化曲线、第六块塔板浓度变化曲线、第一块塔板浓度变化曲线、再沸器浓度变化曲线以及再沸器的温度变化曲线，如下图2—图8所示：
从图2中可以看出在仿真过程中各层的浓度的变化，随着塔板高度增高，浓度是依次增大的，由于冷凝器的滞液量是在变化的，所以浓度比第六层塔板有一定的滞后。随着时间的变化，再沸器中的浓度越来越低，在2500s后蒸发出来的基本上是水蒸气，故冷凝器的浓度开始下降，当水蒸干时，原来的混合液全都进入冷凝器，故5000s时冷凝器浓度在初始值0.45(摩尔分数)。
从图3中可以看出2、4层塔板及冷凝器、再沸器的温度变化，随着精馏过程进行，塔板温度逐渐升高。冷凝器温度降低，在浓度较高时保持在65-70。当蒸发的基本为水时温度回升到初始值附近。
各层r=1.8从开始直至结束的动态浓度变化曲线10.90.80.70.60.50.40.30.20.10 050010001500200025003000Time(s)3500400045005000xdx6x5x4x3x2x1xs x
图2各层从开始到结束的动态浓度变化曲线
各层r=1.8从开始直至结束的动态温度变化曲线10510095908580757065 0tdt4t2ts t50010001500200025003000Time(s)3500400045005000
图3第2，4塔板、冷凝器、再沸器温度变化曲线
冷凝器产品不同回流比从开始直至结束的动态浓度变化曲线0.950.90.850.80.75r=1.8r=2.5 xd0.70.650.60.550.50.45 050010001500200025003000Time(s)3500400045005000
图4不同回流比(r=1.8和r=2.5)下塔顶冷凝器浓度变化曲线
第六块塔板不同回流比从开始直至结束的动态浓度变化曲线10.90.80.70.6r=1.8r=2.5 x60.50.40.30.20.10 050010001500200025003000Time(s)3500400045005000
图5不同回流比(r=1.8和r=2.5)下第六块塔板浓度变化曲线
第一块塔板不同回流比从开始直至结束的动态浓度变化曲线0.7r=1.8r=2.5 0.60.50.4x10.30.20.10 050010001500200025003000Time(s)3500400045005000
图6不同回流比(r=1.8和r=2.5)下第一块塔板浓度变化曲线
再沸器不同回流比从开始直至结束的动态浓度变化曲线0.450.40.350.30.25r=1.8r=2.5 xs0.20.150.10.050 050010001500200025003000Time(s)3500400045005000
图7不同回流比(r=1.8和r=2.5)下再沸器浓度变化曲线
再沸器不同回流比从开始直至结束的动态温度变化曲线105 100r=1.8r=2.59590ts85807570 050010001500200025003000Time(s)3500400045005000
图8不同回流比(r=1.8和r=2.5)下再沸器的温度变化曲线
从图4-8是比较了在两种不同回流比下浓度和温度的变化，可以看出当回流比比较大时在相同时刻各塔板的浓度都比较小回流比时高，但是精馏的过程时间比较长，达到最高浓度的时间也较晚。对于温度，也可以看出较大回流比时的在相同时刻的温度比小回流比的低。
程序如下：
% 初始化
clear all
r=1.8;r1=2.5;
z=0.45;
m=0.01;%
v=3.2/3600;%
d(1,1)=v/(r+1);d(2,1)=v/(r1+1);
l(1,1)=v-d(1,1);l(2,1)=v-d(2,1);
a=3.48;
i=1;
tt=5000;
detat=0.5;
md_0=0.05;%
ms_0=1.6;%
for t=0:detat:tt
if t==0
%初始值
for k=1:2
md(k,i)=md_0;
ms(k,i)=ms_0;
%x
xd(k,i)=z;
x1(k,i)=z;
x2(k,i)=z;
x3(k,i)=z;
x4(k,i)=z;
x5(k,i)=z;
x6(k,i)=z;
xs(k,i)=z;
%y
%yd(i)=a*xd(i)/(1+(a-1)*xd(i));
y1(k,i)=a*x1(k,i)/(1+(a-1)*x1(k,i));
y2(k,i)=a*x2(k,i)/(1+(a-1)*x2(k,i));
y3(k,i)=a*x3(k,i)/(1+(a-1)*x3(k,i));
y4(k,i)=a*x4(k,i)/(1+(a-1)*x4(k,i));
y5(k,i)=a*x5(k,i)/(1+(a-1)*x5(k,i));
y6(k,i)=a*x6(k,i)/(1+(a-1)*x6(k,i));
ys(k,i)=a*xs(k,i)/(1+(a-1)*xs(k,i));
end
else
%滞液量的计算
for k=1:2
md(k,i)=(v-l(k,1))*detat+md(k,i-1);
ms(k,i)=(l(k,1)-v)*detat+ms(k,i-1);
%x浓度的计算
xd(k,i)=(v-l(k,1))*(y6(k,i-1)-xd(k,i-1))*detat/md(k,i-1)+xd(k,i-1);
x1(k,i)=(v*(ys(k,i-1)-y1(k,i-1))+l(k,1)*(x2(k,i-1)-x1(k,i-1)))*detat/m+x1(k,i-1);
x2(k,i)=(v*(y1(k,i-1)-y2(k,i-1))+l(k,1)*(x3(k,i-1)-x2(k,i-1)))*detat/m+x2(k,i-1);
x3(k,i)=(v*(y2(k,i-1)-y3(k,i-1))+l(k,1)*(x4(k,i-1)-x3(k,i-1)))*detat/m+x3(k,i-1);
x4(k,i)=(v*(y3(k,i-1)-y4(k,i-1))+l(k,1)*(x5(k,i-1)-x4(k,i-1)))*detat/m+x4(k,i-1);
x5(k,i)=(v*(y4(k,i-1)-y5(k,i-1))+l(k,1)*(x6(k,i-1)-x5(k,i-1)))*detat/m+x5(k,i-1);
x6(k,i)=(v*(y5(k,i-1)-y6(k,i-1))+l(k,1)*(y6(k,i-1)-x6(k,i-1)))*detat/m+x6(k,i-1);
xs(k,i)=(l(k,1)*(x1(k,i-1)-xs(k,i-1))-v*(ys(k,i-1)-xs(k,i-1)))*detat/ms(k,i-1)+xs(k,i-1);
%y浓度的计算
y1(k,i)=a*x1(k,i)/(1+(a-1)*x1(k,i));
y2(k,i)=a*x2(k,i)/(1+(a-1)*x2(k,i));
y3(k,i)=a*x3(k,i)/(1+(a-1)*x3(k,i));
y4(k,i)=a*x4(k,i)/(1+(a-1)*x4(k,i));
y5(k,i)=a*x5(k,i)/(1+(a-1)*x5(k,i));
y6(k,i)=a*x6(k,i)/(1+(a-1)*x6(k,i));
ys(k,i)=a*xs(k,i)/(1+(a-1)*xs(k,i));
end
end
%温度t的计算
for k=1:2
t2(k,i)=80.1-15.44*x2(k,i)+20*exp(-x2(k,i)/0.1277);
t4(k,i)=80.1-15.44*x4(k,i)+20*exp(-x4(k,i)/0.1277);
td(k,i)=80.1-15.44*xd(k,i)+20*exp(-xd(k,i)/0.1277);
ts(k,i)=80.1-15.44*xs(k,i)+20*exp(-xs(k,i)/0.1277);
end
for k=1:2
if ms(k,i)<1e-5
break;
end
end
i=i+1;
end
%各层r=1.8浓度变化图
figure(1);
t=0:detat:tt;
plot(t,xd(1,:),'b.--',t,x6(1,:),'r-.',t,x5(1,:),'k-.',t,x4(1,:),'c-.',t,x3(1,:),'m-.',t,x2(1,:),'y-.',t,x1(1,:),'g-.',t,xs(1,:),'b-.')
xlabel('Time(s)')
ylabel( 'x')
title('各层r=1.8从开始直至结束的动态浓度变化曲线')
legend('xd', 'x6','x5','x4','x3','x2','x1','xs')
grid on
%r=1.8时2，4，s，d温度变化图
figure(2);
t=0:detat:tt;
plot(t,td(1,:),'b.--',t,t4(1,:),'r-.',t,t2(1,:),'k-.',t,ts(1,:),'g-.')
xlabel('Time(s)')
ylabel( 't')
title('各层r=1.8从开始直至结束的动态温度变化曲线')
legend('td', 't4','t2','ts')
grid on
%冷凝器浓度图
figure(3);
t=0:detat:tt;
plot(t,xd(1,:),'b.-',t,xd(2,:),'r-.')
xlabel('Time(s)')
ylabel( 'xd')
title('冷凝器产品不同回流比从开始直至结束的动态浓度变化曲线')
legend('r=1.8', 'r=2.5')
grid on
%第六块塔板浓度图
figure(4);
t=0:detat:tt;
plot(t,x6(1,:),'b.-',t,x6(2,:),'r-.')
xlabel('Time(s)')
ylabel( 'x6')
title('第六块塔板不同回流比从开始直至结束的动态浓度变化曲线')
legend('r=1.8', 'r=2.5')
grid on
%第一块塔板浓度图
figure(5);
t=0:detat:tt;
plot(t,x1(1,:),'b.-',t,x1(2,:),'r-.')
xlabel('Time(s)')
ylabel( 'x1')
title('第一块塔板不同回流比从开始直至结束的动态浓度变化曲线')
legend('r=1.8', 'r=2.5')
grid on
%再沸器浓度图
figure(6);
t=0:detat:tt;
plot(t,xs(1,:),'b.-',t,xs(2,:),'r-.')
xlabel('Time(s)')
ylabel( 'xs')
title('再沸器不同回流比从开始直至结束的动态浓度变化曲线')
legend('r=1.8', 'r=2.5')
grid on
%再沸器温度变化图
figure(7);
t=0:detat:tt;
plot(t,ts(1,:),'b.-',t,ts(2,:),'r-.')
xlabel('Time(s)')
ylabel( 'ts')
title('再沸器不同回流比从开始直至结束的动态温度变化曲线')
legend('r=1.8', 'r=2.5')
grid on
因篇幅问题不能全部显示，请点此查看更多更全内容