复杂精馏塔的数学模型建立和模拟计算
复杂精馏塔的数学模型建立
和模拟计算
应用化学0311
马迪卓
200349032
前言
简单蒸馏与平衡蒸馏只能将混合物进行初步的分离。为了获得较高纯度的产品，应通过多级的蒸馏过程，使混合物的气、液两相经过多次混合接触和分离。进行质量和热量的传递，使混合物中的组分达到更高程度的分离，这一目标可采用精馏的方法予以实现。
精馏塔是一圆形筒体，塔内装有多层塔板或填料，塔顶设有冷凝器，塔底装有再沸器，塔中部适宜位置设有进料板。
精馏操作是分离液体混合物最为广泛采用的重要分离单元之一。在精馏过程中所采用的塔通常是一股进料分离出两种产品，即塔顶馏出轻组分产品D,塔底获得重组分产品W,称该塔为常规精馏塔或简单精馏塔。
当分离多股组分相同而组分各异的混合物时，常采用多股进料塔。当分离多组分混合物时，其中某产品纯度要求不是很高，其挥发能力仅次于塔顶产品，故可采用抽侧线产品塔，以降低投资费。如果进料中含有少量轻组分，若要将其冷凝下来需要较高品位的冷剂。为此，可采用带部分冷凝器的精馏塔，将该少量轻组分以气象排出，然后单独处理。如果精馏过程允许带水，且要分离比水轻的组分，则可采用直接蒸汽加热方式，这样既提高了塔釜传热效率，又省出一台再沸器。如果混合物中待分离的两组分相对挥发度接近1或形成恒沸物则应选择加入质量分离剂进行分离，如恒沸精馏或萃取精馏。对于那些批量小，组分经常变化的混合物分离通常选择间歇精馏。
为了对上述精馏塔进行模拟计算，则需分别建立数学模型。本文下面就介绍用三对角线矩阵法来实现对复杂精馏塔的数学模拟。
数学建模
1．数学模型
复杂精馏塔模型如图。
设塔含N个理论板或平衡级。理论版序号从上之下编排，塔顶冷凝器序号为1，塔底再沸器序号为N。塔内理论板数为(N-2)该塔的每一个理论板上均有一个进料Fj,气相侧线Uj和一个中间冷却（或加热）器热流量为Qj。塔内的任意理论塔板如图。
平衡级的严格计算应满足4组方程式物料衡算方程（M方程），相平衡方程（E方程），摩尔分数加和方程即归一方程（S方程），热衡算方程（H方程）
A．物料衡算方程（M方程）
对塔内任意一个平衡级j上的组分i进行物料恒算可得
B．相平衡方程（E方程）
C．摩尔分数加和方程即归一方程（S方程）
D. 热衡算方程（H方程）
以上公式中，i=1,2,3,……m
j=1,2,3,……N
x(i,j)——平衡级j上的组分i的液相摩尔分数；
y(i,j)——平衡级j上的组分i的气相摩尔分数；
z(i,j)——平衡级j上的组分i的进料摩尔组成；
K(i,j)——平衡级j上的组分i的平衡常数；
L(j)——平衡级j流出进入下一级的液相流量，kmol/h;
V(j)——平衡级j流出进入上一级的气相流量，kmol/h;
hj(j)——平衡级j上液相摩尔焓,kJ/kmol;
Hj(j)——平衡级j上气相摩尔焓,kJ/kmol;
Hfj(j)——平衡级j上进料摩尔焓,kJ/kmol;
以上M、E、S、H方程分别含mN、mN、N、及N个方程，共含2N（m+1）个方程，该方程组共含变量为（4mN+9N）个。为求解M、E、S、H方程，应根据系统自由度的概念，给定其中2mN+7N个设计变量，使之具有唯一解。在需要给定的设计变量中，部分是由已知条件给定，如Fj、zfi,j、Wj、Vj、Pj等。而另一部分则由补充方程，例如物流焓与组成、温度压力的关系等，确定、、。选择热力学状态方程确定。在求解M、E、S、H中2N（m+1）个方程时由于方程多、变量多，且变量间存在非线性关系，采用解析方法求解较难，通常采用数值方法，按顺序进行求解。
2．数学模型求解
先将M方程、E方程联立，有E方程消去M方程中气相组成并整理得
并将上式表示为以下简单形式：
式中
组分i在塔内各理论级上液相组成的分布关系式可表示为以下非线性方程的形式。
也可表示为以下矩阵形式
或
上面的矩阵即三对角线矩阵方程，是一非线性方程组，其求解过程是：首先根据初值条件求得系数阵A即列矩阵D，使上面的非线性方程组线性化，而线性方程组可采用多种方法求解。
计算步骤
1．给定塔模拟计算的条件
进料液的给定条件：
压力：常压；
烃化液流量：F=150kmol/h;
烃化液温度：T=60度
烃化液组成（摩尔分数）：苯＝0.5982;
甲苯＝0.00602;
邻二甲苯＝0.3091;
乙苯＝0.08668;
2．假定各级的液相组成及温度分布初值
设第一块板的液相组成为 x11=0.98 x21=0.001 x31=0 x41=0
设第NT板的液相组成为 x1NT=0.002; x2NT=0.2035; x3NT=0.3125; x4NT=0.482
设j表示理论板的块数，j=(2:NT-1)。建立各组分在各板上的组成的线性方程：
x1j=0.98-(0.98-0.002).*(j'-1)/(NT-1)
x2j=0.001+(0.2035-0.001).*(j'-1)/(NT-1)
x3j=0.3125.*(j'-1)/(NT-1)
x4j=0.482.*(j'-1)/(NT-1)
经过循环计算得出每块板上各组分的组成。建立矩阵X=[x1j x2j x3j x4j]，便于下面计算中应用。
查表得到各组分的沸点温度：
T1b=353.3; T2b=383.8; T3b=409.3; T4b=417.6;
建立线性方程计算各板上的泡点温度：
Tj=x1j.*T1b+x2j.*T2b+x3j.*T3b+x4j.*T4b
得到矩阵Tj
3．按线性分布方法给定各级上压力分布
假定各板的总压力由以下两个线性方程表示：
j1=(1:14); Pj1=760-(14-j1').*5
j2=(15:32); Pj1=760-(15-j1').*7.6
Pj=[Pj1;Pj2]
4．计算各级泡点确定各级泡点、平衡常数及平衡气相组成
利用安托因方程计算各组分的饱和压力，查表的安托因常数：
A1=15.9008; A2=16.0137; A3=16.0195; A4=16.156;
B1=2788.5; B2=3096.52; B3=3279.47; B4=3395.57;
C1=-52.36; C2=-53.67; C3=-59.95; C4=-59.46;
安托因方程： lnp=A-B/(C+T)
则 P1o=exp(A1-B1./(C1+Tj));
P2o=exp(A2-B2./(C2+Tj));
P3o=exp(A3-B3./(C3+Tj));
P4o=exp(A4-B4./(C4+Tj));
由公式 Ｋij=Pio/Pj 得Kij
则有矩阵 K=[P1o./Pj P2o./Pj P3o./Pj P4o./Pj]
利用公式 Ｙ=K*X 分别计算各组分的气相组成值
Y1=P1o.*x1j./Pj
Y2=P2o.*x2j./Pj
Y3=P3o.*x3j./Pj
Y4=P4o.*x4j./Pj
获得矩阵 Y=[Y1, Y2 ,Y3, Y4]
计算各板上各组分气相组成的和值：
Yj=Y1+Y2+Y3+Y4;
所得矩阵即为Yij的矩阵表示。
线性回归计算abs(yij-1)
当abs(yij-1)>0.0001时，改变Tj的值，进行循环迭代，使abs(yij-1)的值小于0.0001，最后输出Tj，用于进一步的计算。
当abs(yij-1)>0.0001
Tj(j)=Tj(j)-[Yj(j)-1];
最终将Yj(j)收敛于１，计算出泡点温度。物料衡算（M），相平衡衡算（E），摩尔分数加和方程（S）及热量衡算（H）。
泡点计算框图如下：
5．计算各级气，液相及进料的函值
M1=78g/mol; M2=92g/mol; M3=106g/mol; M4=106g/mol
建立矩阵运算求得各级板的摩尔质量平均数，建立的矩阵为M。
M=X*[78 92 106 106]'
查表得到各组分在不同温度下的比热熔Cp,和气化热r。
建立如下求比热熔和汽化热的线性方程：
当 Tj(j)-273.15<100时
Cp1(j)=0.00294*(Tj(j)-273.15)+1.653;
Cp2(j)=0.00362*(Tj(j)-273.15)+1.612;
Cp3(j)=0.00378*(Tj(j)-273.15)+1.668;
Cp4(j)=0.00346*(Tj(j)-273.15)+1.657;
r1(j)=-0.617*(Tj(j)-273.15)+456.38;
r2(j)=-0.544*(Tj(j)-273.15)+434.23;
r3(j)=-0.476*(Tj(j)-273.15)+429.79;
r4(j)=-0.484*(Tj(j)-273.15)+418.75;
当Tj(j)-273.15>100时
Cp1(j)=0.00498*(Tj(j)-273.15)+1.448;
Cp2(j)=0.00400*(Tj(j)-273.15)+1.591;
Cp3(j)=0.00345*(Tj(j)-273.15)+1.690;
Cp4(j)=0.00377*(Tj(j)-273.15)+1.631;
r1(j)=-0.888*(Tj(j)-273.15)+486.70;
r2(j)=-0.694*(Tj(j)-273.15)+450.80;
r3(j)=-0.587*(Tj(j)-273.15)+442.05;
r4(j)=-0.600*(Tj(j)-273.15)+431.60;
由上述的线性方程，就可以计算出每个温度下的比热熔和汽化热。
由下列公式求出各板上的气化焓，和液化焓及进料的摩尔焓。
M= [0.5982 0.00602 0.3091 0.08668]* [78 92 106 106]'
hj=M.*Cpj.*(Tj(j)-273.15);
Hj=hj+r;
Hfj= M.* Cpm* t；
6．计算系数矩阵A及D的值
求解系数阵A中各元素：
当j=1时，冷凝器，其进料Fj及、Wj为0，
则 A1=0
D1=0
式中 U1————为塔顶液相采出流量；
V1————为塔顶气相采出流量。
当2jN-1时，塔内任意理论塔板j则可由上面关系式求得。
当j=N时，塔底再沸器，在再沸器中一般情况下其、侧线采出、均为零，也无下方理论板，故=0，则有
由于塔进料流量及组成均给定，故很容易求得列矩阵D。而系数阵A与塔内气、液流量Vj、Lj及各组分相平衡常数Kij有关，而平衡常数Kij又与各板上温度、压力及物流组成相关，则为一非线性方程。为求解系数阵A，可按照先进行各级泡点计算，然后计算各物流焓，最后确定气、液流量Vj、Lj的顺序求解。
设初值 R（回流比）＝1.65;
V(1)=80;
L(1)=R*V(1);
V(2)=(R+1)*V(1);
运用以下两个公式进行循环迭代：
V(j)=(V(j)*(Hj(j)-hj(j))+L(j-1)*(hj(j)-hj(j-1)))/(Hj(j+1)-hj(j));
L(j)=V(j)-V(1);
由于进料之后精馏塔上下的Ｌ，Ｖ有差别，所以需要进行分段的计算。
7．解三对角线矩阵，获得各板液相组成分布
由
得 对求出的各板液相组成的分布进行圆整：
8．计算收敛判据并判定各板液相组成是否满足计算精度。
如不满足精度要求进行循环计算直到满足精度要求为止，循环计算式如下：
计算结果
见附录
灵敏度分析
图（1） 图（2）
图（3） 图（4）
图（5）
图（1）：代表各组分随塔板数的增加的变化情况。
图（2）：随进料板的改变塔顶采出液的变化情况。
图（3）：表示各板的温度变化 。
图（4）：随塔板数的改变塔顶采出液的变化情况。
图（5）：随回流比的改变塔顶采出液的变化情况。
结论
苯的含量随板数的增加不断减少，甲苯的含量随板数的增加先增加后减少，邻二甲苯的含量随板数的增加不断增加，乙苯的含量随板数的增加不断增加。改变回流比塔顶的苯的组成改变,组成先随回流比的增加而增加,达到一定指之后随回流比的增加而减少。塔顶苯的组成随塔板数的增加而增加，但当塔板数增加到一定值之后组成含量不再增加。
附录
1．程序流程图
2．各组分在不同温度下的比热镕Cp,和气化热r
１——苯，２——甲苯，３——邻二甲苯，４——乙苯
各组分的泡点温度：
T1b=353.3;T2b=383.8;T3b=409.3;T4b=417.6;
3．参考文献
复杂精馏塔的数学模型建立和模拟计算.doc
