分离工程II-07多级分离计算一三对角矩阵法

1第七章多组分多级分离计算

——三对角矩阵法

前章介绍的逐级计算法时是早期提出的，并在实践中不断得到改进的一种方法，特别是收敛法的引入，大大拓宽了适用领域，对于平衡级的物料平衡关系式基本相同的，结构简单的常规塔，用逐级计算法计算常能获得满意的结果。但是，对于功能较多，结构复杂的复述各平衡级的关系式差异较大给逐级计算带来困难，对于这类问题经常采用矩阵法求解。逐级法是逐个方程一2一顺次求解，矩阵法是将全部平衡级的关系式组成联立方程组同时求解的方法，这种方法一次求解的变量多，计算量大，手续复杂，用手算的办法是很难办到的，然而，由于它可以具体描述每一平衡级的不同情况，建立相应的关系式，可以满足多种复杂塔计算的要求，对于复杂冗长的计算，借助电子计算机大容量、快速运算的功能已不成其为问题。因而，矩阵法特别是被称为托马斯算法(Thomasalgorithm)的三对角矩阵法近年来得到广泛的研究于应用。37-1三对角线矩阵算法7-1-1通用模型塔的建立为了设计一个通用的模型塔，设模型塔共有N个理想平衡级，冷凝器为第一级，再沸器为第N级，除第1级和第N级外，每个平衡级都有进料Fj和侧线出料Sj(j=2,3,…,N-1)，并有加热或冷却设备(即有±Qj)。根据具体条件可将塔简化成任何一个实际塔，不需要的量可定为零。图7-1是一个通用模型塔，它也可以简化成简单塔或任何塔。

4图7-1多组分复杂精流塔模型D=V1+S1

V1S1(F1=G1=0)V2L123

LNQC=Q1GNQR=QNN(GN=SN=0)S2S3Sj-1SjSj+1SN-3SN-2SN-1j-1jj+1N-2N-1G2F2G3F3Gj-1Fj-1Gj+1Fj+1GjFN-2GN-2GN-1FN-1

V2Q2

V3Q3

VjQj

VN-1QN-1

LN-1

LN-2

Lj

L2

VN-2QN-2冷凝器分离罐再沸器5图7-2通用平衡级模型图第j级Lyj,Hj,Tj,PjGjFjLj-1Vj+1LjSjQj(+)(-)位头V阀F阀ZjHFjTFjPFjYj+1,Hj+1,Tj+1,Pj+1xj,hj,Tj,Pjxj-1,hj-1,Tj-1,Pj-1Vj6

设想一个功能齐全的通用的平衡级模型如上图所示。用它表示任一种类级联分离设备的任一平衡级，当表示汽液平衡关系时，图中Vj及Lj代表汽、液相流率，下标j表示级序，级序自上而下排列，下标i为组分标号。Qj为传递热量，自本级向外传递以“+”号，向本级传入标“-”号。Hj及hj分别为第j级汽相及液相物流摩尔焓。7

进入本级的原料Fj可以是单相的或是两相的，随组成Zji，摩尔焓HFj，温度TFj及在压力PFj而定，原料压力可以等于本级压力Pj，也可以高于Pj，假如高于Pj则压差(PFj-Pj)将通过F阀绝热减压降为零。 由上部j-1级流来的液流Lj-1，组成为xj-1,i，焓值为hj-1，温度为Tj-1，压力为Pj-1，此压力可能等于或小于本级压力Pj，若Pj-1<Pj，此压力将因经静压头而得以绝热地提高至Pj。8

同样，由下部j+1级流来的汽流Vj+1,i，焓值为Hj+1，温度为Tj+1，压力为Pj+1，其压差(Pj+1-Pj)通过V阀绝热减压而降为零。 具有强度性质yji，Hj，Tj及Pj的蒸汽流Vj离开本级可以分成两股，一股为侧线采出流，流率为Gj，另一股以Vj的流率送往j-1级。而最上一级(j=1)的上升蒸汽即作为产品送出设备。同时，具有强度性质xji，hj，Tj及Pj的液流(Lj+Sj)与气流(Vj+Gj)处于相平衡状态，它们离开本级亦可分成两股，一股为侧线采出流Sj，另一股Lj流向j+1级，对于最下一级(j=N)的下流液体既作为塔底产品离开设备。9

热量Qj或为正或为负值代表段间冷却器或段间加热器的热负荷，而在两端的Q1及QN则为冷凝器及再沸器的热负荷。7-1-2通用数学模型的建立 对于一个平衡级可以建立如下的四类方程，称为MESH方程。组分物料平衡方程—M方程，当物系有C个组分时，每级可有C个方程：(7-1)10

(2)组分相平衡方程—E方程，每级有C个方程：

yji-Kjixji=0式中Kji为相平衡常数(3)摩尔分率加合方程—S方程，每级有C个方程：(4)热量平衡方程—H方程，每级有1个方程：(7-2)(7-3a)(7-3b)(7-4)11将E方程带入M方程，得ME方程：为了简化MESH方程组，将L表示成V的函数，以减少未知量。为此，从冷凝器到第j级做物料恒算，如图7-3所示。从冷凝器至第j级总进料量以FS表示，总气相采出量以GS表示，总液相采出量以SS表示。则总物料恒算得：(7-6)(7-5)12式中：代入式(7-6)得可导出相邻两级汽液相流率的关系式：D=V1+S1V1S1SSGSFSLjVj+1}j图7-3冷凝器到j级物料衡算图(7-10)(7-7)(7-8)(7-9)13以上诸式中的变量不是一成不变的，可以用另外的变量加以取代。例如，可以用组分流率取代摩尔分率vji=Vjyji，lji=Ljxji。7-1-3三对角矩阵算法 三对角矩阵也称Thomas算法，它具有简便、灵活及通用性好的特点，当被引入平衡级计算，便极大地推动了多组分多级计算的发展。 我们从上面的MESH方程可知，任一组分i任一平衡级j(j1且jN)的M方程只相关于相互比14邻的三个平衡级，就是说，第j级只与第j-1级、第j+1极之间存在物料平衡关系，如果将M方程与其他方程分割并以级温度Tj及级流率Vj作迭代变量赋以确定的值，则M方程便是一组以组成xji为变量的线性方程组，由于任一式只涉及三个平衡级，当按级排成矩阵式时，就是具有三对角矩阵的结构，用Gauss消去法不难解出全部。 当式(7-10)代入ME方程，用汽相组成vji代替液相组成lji，消lji得到：15令：(7-11)(7-12)(7-13)(7-14)(7-15)16得到新的ME方程

Ajxj-1,i+Bjxji+Cjxj+1,i=Dj

列成矩阵方程：(7-16)(7-17)17系数矩阵有三对角矩阵结构，当组分数为c时，便可有c个独立的三对角矩阵方程。只要迭代变量Tj及Vj赋以定值，可以算出Kji值，则系数矩阵各元素Aj，Bj，Cj均为已知量，式(7-17)便可用高斯消去法求解。高斯消去法所用公式由变形的M方程式(7-16)导出。j=1B1x1i+C1x2i=D1令q1=D1/B1，p1=C1/B1则上式变为x1i=q1-p1x2i(7-18)18 x2i=q2-p2x3i其中：依此类推，对于任一级(2jN-1)j=jxji=qj-pjxj+1,i其中：j=NxNi=qN其中：(7-18)(7-19)(7-20)1920

计算时，由第1级开始算出q1，p1，接着按顺序递推地算出q2，p2，q3，p3，…一直算出qN，由式(7-20)便求出xNi值，随后，将xNi代入式(7-19)算出xN-1,i，然后反复使用式(7-19)算出xN-2,i，xN-3,i，…，x2,i，x1,i各值。 这种算法运算中通常不涉及两相近数相减，从而避免了舍入误差，而且算出的xji很少出现负值,

加之此法效率高，需要的计算机存储容量小，因此比起矩阵求逆法，此法的优越性是明显的，21虽然如此，此法的应用范围也是有限制的，当所需平衡级很多且各组分的吸收因子在某一段中小于1，而在另一段中大于1的情况，计算将遇到麻烦，必须改用其它更适合的方法。常用的MESH方程组的解法有以下几种：1).顺序收敛法，包括BP法，DP法，SR法2).同时收敛法，包括有多元牛顿法(或2N牛顿法)3).松弛法，主要用于不稳定工况的精流计算。227-2泡点法(BP法)

泡点法(BubblePointMethod)是三对角矩阵多组分多级计算法中的一种计算方案，由于它是在迭代过程中利用三对角矩阵解出各级组成后，通过求泡点的方法确定新一轮各级温度，因而得名。通常，为使迭代计算有稳定的收敛，泡点温度对组成不宜过分敏感，平衡常数K值范围窄的分离物系，由于其组成变化不会引起泡点温度大幅度变化，泡点法对它的计算特别有效，精馏过程多涉及此种物系。23

对于温度对组成非常敏感的情况，例如，甲烷与正丁烷的双组分物系，在28.12kg/mm2压力下，其组成与泡点温度的关系如下：甲烷在溶液中的摩尔分率 泡点温度OC

0.000135.00.018121.10.05493.30.09365.6可见，甲烷溶解于正丁烷中的浓度只要发生少许变化，对泡点温度就带来显著的影响。这种情况用泡点法作迭代计算是很难收敛的。24泡点法的计算步骤1.解三对角矩阵方程(ME方程)求组分分布根据ME方程组中等式左边系数矩阵为三对角矩阵，除主对角线及相邻两条对角线上的元素不全为零外，其余元素为零的特点，用上节课所讲的方法求解出组分浓度，注意，解出的xij(1)需进行圆整。2.用S方程求各板Tj

由各板液相组成xij(1)，利用泡点方程可求出各板温度Tj(1)。对烃类混合物，可用牛顿迭代公式求Tj(1)

，若Tj(1)与原设值一致，则温度计算结束，否则，以Tj(1)作为假设值Tj(0)

，再重复计算，直至Tj(1)与原假设值Tj(0)一致为止。25

Wang和Henke建议用如下较简单的准则作为迭代计算的收敛判据。求出各级温度后，可计算出相平衡常数，再根据相平衡关系计算出yij值。3.用H方程求各板Vj

在计算新一组流率之前，先要用已求出的Tj及xji,yji与求焓关系式hi(Tj)、Hj(Tj)求出每股物流的摩尔焓hj，Hj。将Lj+Sj=Lj-1+Vj+1+Fj-(Vj+Gj)代入H方程，整理得：(7-21)26Vj+1=[(Hj-hj)(Vj+Gj)+(hj-hj-1)Lj-1-(HFj-hj)Fj+Qj]/ [Hj+1-hj] (2jN-1)V1=(S1+V1)-S1=D-S1V2=V1+S1+L1将V1，V2代入H方程，可求出塔顶冷凝器负荷Q1，而塔釜再沸器负荷QN，可由全塔热量平衡求之，将第1级至第N级的H方程相加便导出QN公式如下：

用式(7-22)可计算出各级的汽相流率。若计算出的(7-22)(7-23)(7-24)(7-25)27Vj与假设Vj(0)不一致，以新的Vj做假设值重新计算，直至相符为止。 解出Vj便可求出Lj：

L1=RD Lj=Vj+1+(Fk-Gk-Sk)泡点法的框图如下：(7-27)(7-26)28图7-3BP法计算框图输入物性数据及N，N1(进料板数)假设Tj(0),Vj(0),xjj(0),由Kij=(Tj(0),P,xjj(0))计算Kij解三对角矩阵计算xjj(1)圆整由H方程计算Vj(1)输出结果Tj(0)<=Tj(1)Vj(0)<=Vj(1)xjj(0)<=xjj(1)NYNYT-循环V-循环297-3流率加和法(SR)

SR法的基本要点是利用组分流率概念，将ME方程化成关于组分流率的方程形式。解ME方程求出组分流率lij或vij，组分流率加和得到Lj和Vj，即由ME方程校正总流量Lj和Vj，然后再用热平衡方程(H方程)校正温度Tj。 对于BP法，ME方程求出xij，由S方程校正Tj,最后由H方程校正Lj和Vj；而SR法是由ME方程校正Lj和Vj,而用H方程校正Tj。一般来说，板上温度主要取决于组成，流量取决于热量的情况适宜用BP法。如沸点相近混合物的精30馏计算适用BP法。若流量主要取决于板上组成，而温度主要收热量的影响，如宽沸点范围混合物的精馏计算和吸收计算适宜用SR法。一般用进料的泡点和露点温度之差tDB作为选择BP法或SR法的判据，如果tDB<55oC，选BP法较好，若tDB>55oC，且tDB值越大，用SR法效果越好。至于中等tDB区域，或用2N-牛顿法，或用SR或BP法中的一种加收敛法。31

SR法计算步骤确定原始数据，包括进料变量(流率Fj，组成zji，压力PFj，温度TFj或焓HFj)，各级压力Pj，每一侧线流率Gj及Sj，出入每级传热速率Qj以及总平衡级数。确定迭代变量Vj及Tj的初值，按恒流比及温升与溶解量成比例的原则进行赋值，也可以按恒摩尔流率及线性温度分布的假设来赋值，后一方法虽更粗略，但多数情况下尚能满足需要。323.用托马斯算法求解出液相摩尔分率xji，不经圆整，直接代入下式（7-27），计算出液相流率Lj。

相应的汽相流率可由总物料衡算式(7-28)求得：

然后直接迭代至收敛。4.由H方程校正Tj。具体方法是：令j级热平衡的偏差函数为：(7-27)(7-28)(7-29)33

按牛顿近似法将Ej展开成泰勒级数且只保留一阶偏导数得：对j=1：对2jN-1(7-30)(7-30a)(7-30b)34对j=N对于2jN-1，对照式(7-30)与式(7-30b)并略去后一式中的上标，则有：(7-30c)(7-31a)(7-31b)(7-31c)35式中，偏导数和实际上表示汽、液相在温度T时的比热cV，cL，对于理想混合物：式中，ci为纯组分的比热；角标V，L分别代表汽、液相。

所求的Tj应满足Ej(k+1)=0。由式(7-31)和(7-32)可写出类似ME方程的形式：(7-32a)(7-32b)36同三对角矩阵法的泡点法相类似，用追赶法可解(7-33)37出Tj；第k+1次同第k次Tj的关系为：流率加和法液可以用类似前面所介绍的加速迭代法，对第一级组分流率都用校正因子，其定义是：三对角矩阵解出lij，vij后用j校正，得出新的lij，vij再加和求Lj，Vj。(7-34)(7-35)38输入物性数据及N，进料位置假设Tj(0),Vj(0),xjj(0),计算组分流率ljj(0)及吸收因子Ajj(0)解三对角矩阵算ljj(1),ljj(1),并计算xjj(1)由H方程计算新的Tj(1)输出结果Tj(0)<=Tj(1)Vj(0)<=Vj(1)xjj(0)<=xjj(1)NYNYV-循环T-循环流量加和求总流量，判断计算值与初始值是否一致图7-4SR法计算框图397-42N-牛顿-拉夫森法(简称N-R法)

当溶液为理想溶液或正规溶液时，吸收过程还可以用2N型牛顿-拉夫森法(简称N-R法)来计算。它与三对角矩阵法的根本区别在于后者的Tj，Vj是分别用ME方程或H方程进行校正的，而N-R法中Tj和Vj是同时进行校正的。下面我们介绍这种方法。7-4-1N-R法

N-R法需要上节课所讲的偏差函数Ej：(7-29)40此外，还需另设偏差函数Sj：Sj和Ej可视为Tj及Vj的2N个独立变量的函数，即：按牛顿近似法，将函数Sj，Ej做泰勒级数展开并取一阶偏导数，得：(7-36)(7-37)(7-38)41

(1jN)对各级校正Tj，Vj使：(7-39)(7-40)(7-41)(7-42)42于是，式(7-39,40)可用下面的矩阵形式来表示：(7-43)43或可缩写成： Jkxk=-fk式中：

xk=[T1(k),T2(k),···,TN(k),V1(k),V2(k),···,VN(k)]Tfk=[S1(k),S2(k),···,SN(k),E1(k),E2(k),···,EN(k)]TJk为雅可比矩阵，(7-44)44(7-45)45各级温度和流量的校正值Tj和Vj可通过对Jk求逆而解出。 xk=-Jk-1f或 x(k+1)=x(k)-Jk-1f向量x(k+1)=[T1(k+1),T2(k+1),···,V1(k+1),V2(k+1),···,VN(k+1)]T，解得向量就是求得了新的Tj，Vj，再由ME方程和H方程重新计算Sj和Ej，重复以上计算，直到各级Sj及Ej均近似等于零为止。(7-46)(7-47)46

N-R法的优点是当初值接近于解时收敛速度很快，适用性广，对精馏、吸收等过程均适用；缺点是计算工作量大。因为要计算(2N)2个偏导数，为此至少需对2N+1组独立变量值计算Sj和Ej。另外，一般需经多次迭代方能求解出{Tj}，{Vj}，而每迭代一次，都需要重新计算雅可比矩阵，计算量相当大。为此，C.G.布伊登提出了改进的N-R法。477-4-2改进的N-R法

C.G.布伊登对N-R法的改进主要是在迭代过程中勿需每次计算雅可比矩阵，只要计算一次J中的偏导数即可。 为了便于讨论改进的N-R法，设有任一非线性方程组f(x1,x2,,xN)=0，或写成：

Jkxk=-fk式中： xk=x(k+1)-x(k) x(k)=[x1(k),x2(k),,xN(k),]T(7-48)48

f(k)=[f1(k),f2(k),,fN(k),]T改进的N-R法的计算步骤为：(1).假设一组初值为x0，计算:f0=f(x0)(7-49)(7-50)49(2).计算近似的雅可比矩阵中各元素的值，初始J0中各元素的值用差分法计算如下： 式中，hi=0.0010xi

令H0=-J0-1(3).计算Hk和fk的值并计算

xk=Hkfk(7-51)(7-52)(7-53)50(4).搜寻Sk的值，Sk的选取应满足以下不等式：

Sk的值可在+1和-1之间搜寻，初始值可取S1(k)=1。若不满足上式，可取：

S2(k)=[(1+6)1/2-1]/3

式中：(7-55)(7-54)(7-56)51(5).计算xk+1，fk+1的值

xk+1=xk+Skxk fk+1=f(xk+1)

检验fk+1是否达到收敛要求，若未收敛计算

yk=fk+1-fk(6).计算迭代用的矩阵Hk+1的值

(7).返回步骤(3)重复以上计算。(7-59)(7-56)(7-57)(7-58)527-5松弛法

上面讨论的精馏计算法均属描述稳定状态工况下基本方程的求解法，松弛法是以不稳定状态下的物料平衡为基础进行精馏计算的方法之一。7-5-1基本关系式 图7-5表示塔第j级，Uj表示级上的存液量。当第j级处于不稳定状态时，在t时间内，由于物料进出第j平衡级的I组分物料不等，引起i组分在级内积累，导致其浓度由原来的xji(k)变化为xji(k+1)，可以用下式描述其基本关系式。53图7-5非稳态操作物料平衡图Vj·yjiLj-1·xj-1,iLj·xj,iVj+1·yj+1,i第j级Uj·xj,i54

[进料]-[出料]=[积累]其中：

则，j级i组分的物料平衡如下：(7-60)(7-61)(7-62)(7-63)55积分中值定理：当函数F(t)在区间(t,t+t)内连续，则在此区间必有一点(t+t)(01)使下式成立：微分中值定理：当函数f(t)在区间[t,t+t]内连续，则在此区间必有一点(t+

t)(01)使下式成立：(7-65)(7-64)56应用中值定理，则式(7-)可变为：取Uj为常量，且当t0时，(t+t)t，各量均趋近于t处的瞬时值，则上式可化为：(7-67)(7-66)57将导数写成差分形式，将经期间浓度的变化比作一次迭代计算后的浓度值变化，将表示时间的下标t及换成标明迭代次数的符号(k)及(k+1)，则上式化为：其中，j=t/Uj，是松弛系数。简单记为：xji(k+1)=xji(k)+jGji(7-68)(7-69)58其中，为因进出j级的i组分物料量不等而引起的积累量。通过调整j值的大小，来控制xji值经一次迭代计算后的变化幅度。运用此式，逐一组分逐级地反复进行计算，每做一轮计算，校正一次浓度，直到级内浓度不再变化，xji(k+1)=xji(k)，表明进出该级组分物料相等，Gji=0，达到最终收敛。(7-70)59

整个运算过程，正如一个平衡级分离设备从开车投料的非稳态操作到稳态操作过程。这种方法的最大优点，在于收敛稳定性好，即使用于高度非理想物系的多级分离计算，一般也能稳定收敛。607-5-2计算公式及计算步骤

松弛法具体计算与各平衡级实际的物料出入情况直接有关，以一个有N个平衡级(包括全凝器及再沸器)，多股进料的精馏塔为例，其计算公式，根据各级情况，列之如下：第1级(全凝器)：(7-71a)61

LNxN,i

VNyN,i

LN-1xN-1,iFj,zjiHFiN

VNyN,i

LN-1xN-1,i

Vj+1

yj+1,i

Lj-1xj-1,i

Ljxj,i

Vjyj,i

V2y2,i

L1x1,i

Dx1,i1j图7-6松弛法模型塔62第j级:(2jN-1)第N级(再沸器)：(7-71c)(7-71b)63计算时，逐个成分，自塔顶开始，逐级向下计算，直至塔釜。为加速收敛，每算一级得出i的新的组成xj-1,i(k+1)的值可以直接代入下一级的计算公式，求取下一级新的xj,i(k+1)，这样，式(7-71b)及式(7-71c)则可变成如下形式：(7-72b)64式(7-71a)则保持原来各项之值，可写为下列形式：(7-72a)(7-72c)65由公式可知，松弛法是一种操作型算法，要预先规定平衡级数N及各进料位置，此外，还必须预先设定汽液相流率(Vj及Lj)，温度(Tj)及组成(xji)的初值。流率初值按恒摩尔流假设进行估算，温度初值在塔顶温度及塔釜温