分離工程II-07多級分離計算一三對角矩陣法

1第七章多組分多級分離計算

——三對角矩陣法

前章介紹的逐級計算法時是早期提出的，并在實踐中不斷得到改進的一種方法，特別是收斂法的引入，大大拓寬了適用領(lǐng)域，對于平衡級的物料平衡關(guān)系式基本相同的，結(jié)構(gòu)簡單的常規(guī)塔，用逐級計算法計算常能獲得滿意的結(jié)果。但是，對于功能較多，結(jié)構(gòu)復(fù)雜的復(fù)述各平衡級的關(guān)系式差異較大給逐級計算帶來困難，對于這類問題經(jīng)常采用矩陣法求解。逐級法是逐個方程一2一順次求解，矩陣法是將全部平衡級的關(guān)系式組成聯(lián)立方程組同時求解的方法，這種方法一次求解的變量多，計算量大，手續(xù)復(fù)雜，用手算的辦法是很難辦到的，然而，由于它可以具體描述每一平衡級的不同情況，建立相應(yīng)的關(guān)系式，可以滿足多種復(fù)雜塔計算的要求，對于復(fù)雜冗長的計算，借助電子計算機大容量、快速運算的功能已不成其為問題。因而，矩陣法特別是被稱為托馬斯算法(Thomasalgorithm)的三對角矩陣法近年來得到廣泛的研究于應(yīng)用。37-1三對角線矩陣算法7-1-1通用模型塔的建立為了設(shè)計一個通用的模型塔，設(shè)模型塔共有N個理想平衡級，冷凝器為第一級，再沸器為第N級，除第1級和第N級外，每個平衡級都有進料Fj和側(cè)線出料Sj(j=2,3,…,N-1)，并有加熱或冷卻設(shè)備(即有±Qj)。根據(jù)具體條件可將塔簡化成任何一個實際塔，不需要的量可定為零。圖7-1是一個通用模型塔，它也可以簡化成簡單塔或任何塔。

4圖7-1多組分復(fù)雜精流塔模型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

設(shè)想一個功能齊全的通用的平衡級模型如上圖所示。用它表示任一種類級聯(lián)分離設(shè)備的任一平衡級，當(dāng)表示汽液平衡關(guān)系時，圖中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，此壓力將因經(jīng)靜壓頭而得以絕熱地提高至Pj。8

同樣，由下部j+1級流來的汽流Vj+1,i，焓值為Hj+1，溫度為Tj+1，壓力為Pj+1，其壓差(Pj+1-Pj)通過V閥絕熱減壓而降為零。 具有強度性質(zhì)yji，Hj，Tj及Pj的蒸汽流Vj離開本級可以分成兩股，一股為側(cè)線采出流，流率為Gj，另一股以Vj的流率送往j-1級。而最上一級(j=1)的上升蒸汽即作為產(chǎn)品送出設(shè)備。同時，具有強度性質(zhì)xji，hj，Tj及Pj的液流(Lj+Sj)與氣流(Vj+Gj)處于相平衡狀態(tài)，它們離開本級亦可分成兩股，一股為側(cè)線采出流Sj，另一股Lj流向j+1級，對于最下一級(j=N)的下流液體既作為塔底產(chǎn)品離開設(shè)備。9

熱量Qj或為正或為負值代表段間冷卻器或段間加熱器的熱負荷，而在兩端的Q1及QN則為冷凝器及再沸器的熱負荷。7-1-2通用數(shù)學(xué)模型的建立 對于一個平衡級可以建立如下的四類方程，稱為MESH方程。組分物料平衡方程—M方程，當(dāng)物系有C個組分時，每級可有C個方程：(7-1)10

(2)組分相平衡方程—E方程，每級有C個方程：

yji-Kjixji=0式中Kji為相平衡常數(shù)(3)摩爾分率加合方程—S方程，每級有C個方程：(4)熱量平衡方程—H方程，每級有1個方程：(7-2)(7-3a)(7-3b)(7-4)11將E方程帶入M方程，得ME方程：為了簡化MESH方程組，將L表示成V的函數(shù)，以減少未知量。為此，從冷凝器到第j級做物料恒算，如圖7-3所示。從冷凝器至第j級總進料量以FS表示，總氣相采出量以GS表示，總液相采出量以SS表示。則總物料恒算得：(7-6)(7-5)12式中：代入式(7-6)得可導(dǎo)出相鄰兩級汽液相流率的關(guān)系式：D=V1+S1V1S1SSGSFSLjVj+1}j圖7-3冷凝器到j(luò)級物料衡算圖(7-10)(7-7)(7-8)(7-9)13以上諸式中的變量不是一成不變的，可以用另外的變量加以取代。例如，可以用組分流率取代摩爾分率vji=Vjyji，lji=Ljxji。7-1-3三對角矩陣算法 三對角矩陣也稱Thomas算法，它具有簡便、靈活及通用性好的特點，當(dāng)被引入平衡級計算，便極大地推動了多組分多級計算的發(fā)展。 我們從上面的MESH方程可知，任一組分i任一平衡級j(j1且jN)的M方程只相關(guān)于相互比14鄰的三個平衡級，就是說，第j級只與第j-1級、第j+1極之間存在物料平衡關(guān)系，如果將M方程與其他方程分割并以級溫度Tj及級流率Vj作迭代變量賦以確定的值，則M方程便是一組以組成xji為變量的線性方程組，由于任一式只涉及三個平衡級，當(dāng)按級排成矩陣式時，就是具有三對角矩陣的結(jié)構(gòu)，用Gauss消去法不難解出全部。 當(dāng)式(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系數(shù)矩陣有三對角矩陣結(jié)構(gòu)，當(dāng)組分數(shù)為c時，便可有c個獨立的三對角矩陣方程。只要迭代變量Tj及Vj賦以定值，可以算出Kji值，則系數(shù)矩陣各元素Aj，Bj，Cj均為已知量，式(7-17)便可用高斯消去法求解。高斯消去法所用公式由變形的M方程式(7-16)導(dǎo)出。j=1B1x1i+C1x2i=D1令q1=D1/B1，p1=C1/B1則上式變?yōu)閤1i=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，然后反復(fù)使用式(7-19)算出xN-2,i，xN-3,i，…，x2,i，x1,i各值。 這種算法運算中通常不涉及兩相近數(shù)相減，從而避免了舍入誤差，而且算出的xji很少出現(xiàn)負值,

加之此法效率高，需要的計算機存儲容量小，因此比起矩陣求逆法，此法的優(yōu)越性是明顯的，21雖然如此，此法的應(yīng)用范圍也是有限制的，當(dāng)所需平衡級很多且各組分的吸收因子在某一段中小于1，而在另一段中大于1的情況，計算將遇到麻煩，必須改用其它更適合的方法。常用的MESH方程組的解法有以下幾種：1).順序收斂法，包括BP法，DP法，SR法2).同時收斂法，包括有多元牛頓法(或2N牛頓法)3).松弛法，主要用于不穩(wěn)定工況的精流計算。227-2泡點法(BP法)

泡點法(BubblePointMethod)是三對角矩陣多組分多級計算法中的一種計算方案，由于它是在迭代過程中利用三對角矩陣解出各級組成后，通過求泡點的方法確定新一輪各級溫度，因而得名。通常，為使迭代計算有穩(wěn)定的收斂，泡點溫度對組成不宜過分敏感，平衡常數(shù)K值范圍窄的分離物系，由于其組成變化不會引起泡點溫度大幅度變化，泡點法對它的計算特別有效，精餾過程多涉及此種物系。23

對于溫度對組成非常敏感的情況，例如，甲烷與正丁烷的雙組分物系，在28.12kg/mm2壓力下，其組成與泡點溫度的關(guān)系如下：甲烷在溶液中的摩爾分率 泡點溫度OC

0.000135.00.018121.10.05493.30.09365.6可見，甲烷溶解于正丁烷中的濃度只要發(fā)生少許變化，對泡點溫度就帶來顯著的影響。這種情況用泡點法作迭代計算是很難收斂的。24泡點法的計算步驟1.解三對角矩陣方程(ME方程)求組分分布根據(jù)ME方程組中等式左邊系數(shù)矩陣為三對角矩陣，除主對角線及相鄰兩條對角線上的元素不全為零外，其余元素為零的特點，用上節(jié)課所講的方法求解出組分濃度，注意，解出的xij(1)需進行圓整。2.用S方程求各板Tj

由各板液相組成xij(1)，利用泡點方程可求出各板溫度Tj(1)。對烴類混合物，可用牛頓迭代公式求Tj(1)

，若Tj(1)與原設(shè)值一致，則溫度計算結(jié)束，否則，以Tj(1)作為假設(shè)值Tj(0)

，再重復(fù)計算，直至Tj(1)與原假設(shè)值Tj(0)一致為止。25

Wang和Henke建議用如下較簡單的準則作為迭代計算的收斂判據(jù)。求出各級溫度后，可計算出相平衡常數(shù)，再根據(jù)相平衡關(guān)系計算出yij值。3.用H方程求各板Vj

在計算新一組流率之前，先要用已求出的Tj及xji,yji與求焓關(guān)系式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方程相加便導(dǎo)出QN公式如下：

用式(7-22)可計算出各級的汽相流率。若計算出的(7-22)(7-23)(7-24)(7-25)27Vj與假設(shè)Vj(0)不一致，以新的Vj做假設(shè)值重新計算，直至相符為止。 解出Vj便可求出Lj：

L1=RD Lj=Vj+1+(Fk-Gk-Sk)泡點法的框圖如下：(7-27)(7-26)28圖7-3BP法計算框圖輸入物性數(shù)據(jù)及N，N1(進料板數(shù))假設(shè)Tj(0),Vj(0),xjj(0),由Kij=(Tj(0),P,xjj(0))計算Kij解三對角矩陣計算xjj(1)圓整由H方程計算Vj(1)輸出結(jié)果Tj(0)<=Tj(1)Vj(0)<=Vj(1)xjj(0)<=xjj(1)NYNYT-循環(huán)V-循環(huán)297-3流率加和法(SR)

SR法的基本要點是利用組分流率概念，將ME方程化成關(guān)于組分流率的方程形式。解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法的判據(jù)，如果tDB<55oC，選BP法較好，若tDB>55oC，且tDB值越大，用SR法效果越好。至于中等tDB區(qū)域，或用2N-牛頓法，或用SR或BP法中的一種加收斂法。31

SR法計算步驟確定原始數(shù)據(jù)，包括進料變量(流率Fj，組成zji，壓力PFj，溫度TFj或焓HFj)，各級壓力Pj，每一側(cè)線流率Gj及Sj，出入每級傳熱速率Qj以及總平衡級數(shù)。確定迭代變量Vj及Tj的初值，按恒流比及溫升與溶解量成比例的原則進行賦值，也可以按恒摩爾流率及線性溫度分布的假設(shè)來賦值，后一方法雖更粗略，但多數(shù)情況下尚能滿足需要。323.用托馬斯算法求解出液相摩爾分率xji，不經(jīng)圓整，直接代入下式（7-27），計算出液相流率Lj。

相應(yīng)的汽相流率可由總物料衡算式(7-28)求得：

然后直接迭代至收斂。4.由H方程校正Tj。具體方法是：令j級熱平衡的偏差函數(shù)為：(7-27)(7-28)(7-29)33

按牛頓近似法將Ej展開成泰勒級數(shù)且只保留一階偏導(dǎo)數(shù)得：對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式中，偏導(dǎo)數(shù)和實際上表示汽、液相在溫度T時的比熱cV，cL，對于理想混合物：式中，ci為純組分的比熱；角標V，L分別代表汽、液相。

所求的Tj應(yīng)滿足Ej(k+1)=0。由式(7-31)和(7-32)可寫出類似ME方程的形式：(7-32a)(7-32b)36同三對角矩陣法的泡點法相類似，用追趕法可解(7-33)37出Tj；第k+1次同第k次Tj的關(guān)系為：流率加和法液可以用類似前面所介紹的加速迭代法，對第一級組分流率都用校正因子，其定義是：三對角矩陣解出lij，vij后用j校正，得出新的lij，vij再加和求Lj，Vj。(7-34)(7-35)38輸入物性數(shù)據(jù)及N，進料位置假設(shè)Tj(0),Vj(0),xjj(0),計算組分流率ljj(0)及吸收因子Ajj(0)解三對角矩陣算ljj(1),ljj(1),并計算xjj(1)由H方程計算新的Tj(1)輸出結(jié)果Tj(0)<=Tj(1)Vj(0)<=Vj(1)xjj(0)<=xjj(1)NYNYV-循環(huán)T-循環(huán)流量加和求總流量，判斷計算值與初始值是否一致圖7-4SR法計算框圖397-42N-牛頓-拉夫森法(簡稱N-R法)

當(dāng)溶液為理想溶液或正規(guī)溶液時，吸收過程還可以用2N型牛頓-拉夫森法(簡稱N-R法)來計算。它與三對角矩陣法的根本區(qū)別在于后者的Tj，Vj是分別用ME方程或H方程進行校正的，而N-R法中Tj和Vj是同時進行校正的。下面我們介紹這種方法。7-4-1N-R法

N-R法需要上節(jié)課所講的偏差函數(shù)Ej：(7-29)40此外，還需另設(shè)偏差函數(shù)Sj：Sj和Ej可視為Tj及Vj的2N個獨立變量的函數(shù)，即：按牛頓近似法，將函數(shù)Sj，Ej做泰勒級數(shù)展開并取一階偏導(dǎo)數(shù)，得：(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，重復(fù)以上計算，直到各級Sj及Ej均近似等于零為止。(7-46)(7-47)46

N-R法的優(yōu)點是當(dāng)初值接近于解時收斂速度很快，適用性廣，對精餾、吸收等過程均適用；缺點是計算工作量大。因為要計算(2N)2個偏導(dǎo)數(shù)，為此至少需對2N+1組獨立變量值計算Sj和Ej。另外，一般需經(jīng)多次迭代方能求解出{Tj}，{Vj}，而每迭代一次，都需要重新計算雅可比矩陣，計算量相當(dāng)大。為此，C.G.布伊登提出了改進的N-R法。477-4-2改進的N-R法

C.G.布伊登對N-R法的改進主要是在迭代過程中勿需每次計算雅可比矩陣，只要計算一次J中的偏導(dǎo)數(shù)即可。 為了便于討論改進的N-R法，設(shè)有任一非線性方程組f(x1,x2,,xN)=0，或?qū)懗桑?/p>

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).假設(shè)一組初值為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的選取應(yīng)滿足以下不等式：

Sk的值可在+1和-1之間搜尋，初始值可取S1(k)=1。若不滿足上式，可?。?/p>

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)重復(fù)以上計算。(7-59)(7-56)(7-57)(7-58)527-5松弛法

上面討論的精餾計算法均屬描述穩(wěn)定狀態(tài)工況下基本方程的求解法，松弛法是以不穩(wěn)定狀態(tài)下的物料平衡為基礎(chǔ)進行精餾計算的方法之一。7-5-1基本關(guān)系式 圖7-5表示塔第j級，Uj表示級上的存液量。當(dāng)?shù)趈級處于不穩(wěn)定狀態(tài)時，在t時間內(nèi)，由于物料進出第j平衡級的I組分物料不等，引起i組分在級內(nèi)積累，導(dǎo)致其濃度由原來的xji(k)變化為xji(k+1)，可以用下式描述其基本關(guān)系式。53圖7-5非穩(wěn)態(tài)操作物料平衡圖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積分中值定理：當(dāng)函數(shù)F(t)在區(qū)間(t,t+t)內(nèi)連續(xù)，則在此區(qū)間必有一點(t+t)(01)使下式成立：微分中值定理：當(dāng)函數(shù)f(t)在區(qū)間[t,t+t]內(nèi)連續(xù)，則在此區(qū)間必有一點(t+

t)(01)使下式成立：(7-65)(7-64)56應(yīng)用中值定理，則式(7-)可變?yōu)椋喝j為常量，且當(dāng)t0時，(t+t)t，各量均趨近于t處的瞬時值，則上式可化為：(7-67)(7-66)57將導(dǎo)數(shù)寫成差分形式，將經(jīng)期間濃度的變化比作一次迭代計算后的濃度值變化，將表示時間的下標t及換成標明迭代次數(shù)的符號(k)及(k+1)，則上式化為：其中，j=t/Uj，是松弛系數(shù)。簡單記為：xji(k+1)=xji(k)+jGji(7-68)(7-69)58其中，為因進出j級的i組分物料量不等而引起的積累量。通過調(diào)整j值的大小，來控制xji值經(jīng)一次迭代計算后的變化幅度。運用此式，逐一組分逐級地反復(fù)進行計算，每做一輪計算，校正一次濃度，直到級內(nèi)濃度不再變化，xji(k+1)=xji(k)，表明進出該級組分物料相等，Gji=0，達到最終收斂。(7-70)59

整個運算過程，正如一個平衡級分離設(shè)備從開車投料的非穩(wěn)態(tài)操作到穩(wěn)態(tài)操作過程。這種方法的最大優(yōu)點，在于收斂穩(wěn)定性好，即使用于高度非理想物系的多級分離計算，一般也能穩(wěn)定收斂。607-5-2計算公式及計算步驟

松弛法具體計算與各平衡級實際的物料出入情況直接有關(guān)，以一個有N個平衡級(包括全凝器及再沸器)，多股進料的精餾塔為例，其計算公式，根據(jù)各級情況，列之如下：第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由公式可知，松弛法是一種操作型算法，要預(yù)先規(guī)定平衡級數(shù)N及各進料位置，此外，還必須預(yù)先設(shè)定汽液相流率(Vj及Lj)，溫度(Tj)及組成(xji)的初值。流率初值按恒摩爾流假設(shè)進行估算，溫度初值在塔頂溫度及塔釜溫