高等化工熱力學(xué)-2相平衡計算

ProIIAspen等，可以提供很完善的計算方法。本章簡單介紹相平衡計算的基本原理，[2-6,7]，這部分內(nèi)容在其他專著中介紹得相對相平衡計算——有模還應(yīng)輸入那些為表征系統(tǒng)所必需的性質(zhì)。本節(jié)討論的相平衡計算，主要采用模型來輸入那些性質(zhì)。kx()x()z，求分相后各相組成x()和x()。定義組分k在相和K(αβ為：kK(αβ)x(α)/

對于計算相平衡的問題，最切題的普遍熱力學(xué)關(guān)系式即相平衡判據(jù)。按式(1-5.15)()()，f()

f()，k1,2,,

相采用統(tǒng)一的狀態(tài)方程。以式(1-6.21)代入式(2-px()()px()(

K(αβ)(β)/(α)

px()()f*()x()(

kK()x()(

Hk kK()(m()/mo)(

kK()(c()/co)(

kK()f*()()/(p()

k K()K()()/(p()

Hkk k如應(yīng)用第III或第IVK(k和相均采用活度因子關(guān)聯(lián)式。如均采用第I種活度，式(2-1.2)f*x()()f*x()(

k k k kK()()/(

k k如相用第I種活度，相用第II種活度，式(2-1.2)f*()x()()K()x()(

k

Hk

kK()K()()/(f*()()

k

kk如應(yīng)用第III或第IVK(k氣液平、氣液平衡的計算有三種即泡點計算計算和閃蒸計算下面以xx(L)，yx(V)、k x(g)，k

K(VL)分別代表液相組成、氣相組成和氣液平衡常

kkf(V)kk

kyk

yk1，以式(2-1.1)KKKkxkk

Kk可采用上面介紹過的任一個式子計算。例如采用式(2-1.4)K((L)K

/(V))

k式中L是T、p、x的函數(shù)，V是T、p、y的函數(shù)。運(yùn)算時，先代入y pTyk=Kkxkyk需要使用狀態(tài)方程，它就是為表征系統(tǒng)所應(yīng)輸入的模型。在使用狀態(tài)方程計算k時，例如使用式(1-6.45)，在積分下限和壓縮因子Z中，包含有未知數(shù)體積V(V)或V(L)，也要用狀態(tài)方Kf*(L)x(L)/(p(V))K

kk

kk k式中k,I是T、p、x的函數(shù)，迭代運(yùn)算與式(2-1.16)相同。除了計算(V)需要狀態(tài)方程外，計k(L)需要活度因子關(guān)聯(lián)式，它們都是為表征系統(tǒng)所應(yīng)輸入的模型。此外還需要f*(L)，它是k 統(tǒng)溫度壓力的函數(shù)，可利用式(1-6.6)f*(L)(T,p)f*(T,

k k

式中p*和V*(L)k的飽和蒸氣壓和摩爾體積，exp項常稱為Poynting因子。f*(T f*(L)(T,p*)f*(V)(T,p*)p**(T,

k 式中*(T,p*是純物質(zhì)k (2-1.10)，需要知道KHk與應(yīng)改為V

2-1計輸入的獨立變量為y以及T或p，要求輸出x以及p或T。根據(jù)x必須滿足歸一化條件xKxk1

1，以式(2-1.1)KKyi/Kii稱為方程。具體計算和求解泡點方程時類似

圖2-2壓力和液相組成的計算框合物的總組成z，輸出為x和y以及汽相分?jǐn)?shù)，定義為n(V)/

有一點需要說明，按相律，F(xiàn)=K-2+2-0-0=K，但現(xiàn)在輸入了K+1個變量，因此輸出的不僅xzk(1)xkyk(1Kk

kxk

xk1KKzk/(1Kk)k

子關(guān)聯(lián)式以表征系統(tǒng)。首先可假設(shè)x和y的初值，利用閃蒸方程解得后，代入式(2-1.22)得到x和y的新值，反復(fù)迭代直至收斂。液液平T、pz，輸出分相后兩液相的組成x()、x()以及相分率，n()/

K(x(x(Kk，根據(jù)物料衡算并以式(2-1.1) z(1)x()x()(1

)x(

kk由于x()必需歸一，kk

x()1KKzk/(1Kk)k

得x()和x()。液固平以某固體組分1在液態(tài)溶劑(可以是純?nèi)軇┗蚧旌先軇?理，有f(S)

f(Lf2的等式，并且可寫出f(S)

f*(S) 1后者是純固體的特性，并依賴于系統(tǒng)的溫度和壓力。對于f(L111的液體，因此通常采用第IIf(L)1

11力的函數(shù)。代入f(S)11

f(L)1f*(S)1

1變量為T和p，輸入后即可輸出溶解度，即x1。1f*(SKH1并不是唯一的選擇，另一種更實用的方法是：對f(Lf(L)

f*(L)x

，f*(L1的逸度。代入f(S)

f(L) 1 f*(S)f*(L)x

1比值f*(L)/f*(S1的熔化熱fusHm,1和液態(tài)、固態(tài)的恒壓熱容

、

系統(tǒng)壓力下純組分1的為Tf,1，由于壓力對影響較小，Tf,1可近似地用常壓或Tt,1Tf,1下，fusGm,1(Tf,1)=0。而在系統(tǒng)溫度壓力下，fusGm,1f*(Lf*(S (T)G*(L)G*(S)RTln(f*(L)/f*(S)

由于[(GT

T]p

T2(fusGm,1/T)

fus

如設(shè)

T

代入式(2-1.30)Tf,1T

lnx1fusGm,1(T)RTlnfusHm,1(Tf,1) 1 fusCpm,1

TT TT

1ln f,1 以式(2-1.32)代入即可計算溶解度x1。Tf,1、fusHm,1(Tf,1)、

、、

可代替f*(S和KH111如果溶液為理想溶液，則式(2-1.33)lnx1

fusHm

TTTT

fusCpm,1

TT

f,1 性質(zhì)無關(guān)。一般情況下，fusCpm,1很小，對溶解度計算的影響也不大，常?？梢院雎?。劑1，即fusHm,2(Tf,2) 1 fusCpm,2

lnx2

1ln

TT f,2 TT當(dāng)組分1和組分2同時析出時，稱為最低共，其溫度和組成稱為低共熔溫度和低共氣液相平衡計算——由T、p、x推算氣液平衡數(shù)據(jù)最常規(guī)的獲得方法是采用各種循環(huán)式平衡釜抽取氣、液相樣品進(jìn)行相組成分析，并同時測定平衡時的溫度T和壓力p，這是一種直接測定法。一些學(xué)者曾對不同氣液[2-6~2-9]Tpxy(狀態(tài)方程或/和過量函數(shù)模型)，即可根據(jù)相平衡原理計算得到另兩個變量，這就是前面所討論的氣液相平衡的有模型計算法。如果不是輸入分子熱力學(xué)模型而是通過實驗測定第三個變量，則原則上第四個變量可由它們通過相平衡計算確定。由于計算過程中沒有引入模型，稱為氣液平衡的無模型計算法。由于實際測定中氣相取樣和組成測定是最的，一般是恒定溫度T下靜態(tài)法測定一系列液相組成x對應(yīng)的系統(tǒng)總壓p(溶液的飽和蒸汽壓)pT(溶液的沸點)，簡單地說就是T、p、x推算所謂無模型法并不是有時仍需要采用合適的模型計算一部分非關(guān)鍵性的決定系統(tǒng)特征的性質(zhì)，例如氣相逸度因子等?？傊辽偈谴蟠鬁p少了模型的使用。國內(nèi)外學(xué)者已建T、p、xym(G)方程的基礎(chǔ)上建立起來的。根據(jù)應(yīng)用-D(16.1)表示的逸度的-DTpxyT、p、x推算y，根據(jù)Q與活度因子的關(guān)系(隱含了D方程)Q函數(shù)法(間接法2-2.1.1Qk汽液平衡時，按判據(jù)式(1-6.13)f(V)k

kk第Ipyp**x

exp[V*L(pp*)/RT]，k1,2,,

k

kkk

KpK

p**x

exp[V*L(pp*)/RT]/

KkKk

kkk

k

式中k,I用式(1-7.46、7.47)代入，Kp*

exp[V*L(pp*)/RT

K

Qpkk exp

xj k

k

j

xHET

K T

x[k,K VEp

x[j,KjK p j

m

xj

m

xj RT2x

j

x

RTx

j

x kx[k,K

jx[j,K]

kx[k,K

jx[j,K] HE和VEQT、p、xQT p、x的函數(shù)，式(2-2.3)QT、p即可求得yp*、*和V*L 關(guān)。k決定于氣相組成y，可利用上次迭代的y值計算，但還需要使用合適的狀態(tài)方程，從這個意義上說，T、p、x推算y并不是完全的無模型，但當(dāng)壓力不太高時，氣相非理想性遠(yuǎn)理想性，甚至可以令k=1，也不致帶入嚴(yán)重誤差。至于HE和VE，后者很小，?？珊雎裕?式(2-2.3)原則可以求解，但實踐上卻有很大，因為導(dǎo)數(shù)出現(xiàn)在exp中，是一個NewtonQ值。這種方法不依賴于任何過量函數(shù)模型，是嚴(yán)格的在一個框架下。大量實例計算表明，沒有收斂的，不受多元系Q函數(shù)曲面類型的限制。下面介紹計算簡單的Barker法和適用于多元系的曲面樣條函數(shù)法，至于其它方法，感興2-2.1.2BarkerBarker法[2-10,11]的是選擇一個合適的過量函數(shù)模型，其中包括若干待定模型參數(shù)，代入式(2-2.3)T、p、xQ函數(shù)。4Barker法的候選模型。但Barker法的準(zhǔn)Redlich-KisterQ12jQxx12j

jAj

x)j

1212

x2

jAj

x)j2xx

jjAj

x)

2jj2jj

x2

1212jA1212j

x)j2xx

1jjA1j

x)

221j2112j21T、p、1j2112j21采用最小二乘法關(guān)聯(lián)得到N+1個Aj。具體計算時，我們可以采用如下的迭代過程：k由狀態(tài)方程計算氣相逸度因子*和k，由式(2-2.7)計算各組分的液相活度系數(shù)k py

k

k kk

由式(2-2.8)KKQxklnkQ關(guān)聯(lián)式(2-2.4)Aj(j=0,…,由式(2-2.9)

yp**x

exp[V*L(pp*)/RT]/

kk

k

比較y1和y0，如果兩者不相等，則令y0=y1，轉(zhuǎn)步驟(2)，進(jìn)行新一輪循環(huán)迭代，直2-2.1.3QQx1x2Q(x1,x2)，Q-x1-x2形成三維空間，對于液相完全互溶的系統(tǒng)，Q函數(shù)是一個定義在[x10;x20;x1+x21]三角形區(qū)域內(nèi)的Q=0QNQ函數(shù)值已Qi(x1(i),x2)Qi，可將Q看成是一塊無限大平板純彎曲時的變形[2-19]，Q上的負(fù)載q 2.10)N個NQ(x,x)[dR2ln(R2)] N

x

i

i

N1

N2

NR2(x )2(x

Q N

2di(x1x1(i))

)dN

x

i

R2 2 1 2Q N

2di(x2x2(i))

)dN

i

R2 1 2 1式中1~10-110-4~1-6，一般可取1~1-2，奇性曲面10-5~1-6。di是未定系數(shù)，應(yīng)滿足以下三個條件，NNG1x1(i)dii

NNG2x2(i)dii

NNG3dii

NQdR2ln(R2) N

ji

i

N11(j

N22(j

N R2 )2

1(j 2(j

kjj的彈性常數(shù)，cj是對曲面光順要求所加的權(quán)重，在某點上該值取為零，則所構(gòu)作的Q(x1(i),x2(i),Qi)(2-7)和(2-2.18)N+3N+3diQ(x1,x2)QD(Q)

如果D(Q)是個線性偏微分方程，得到的是線性方程組，反之則得到非線性方程組。求解由式(2-2.22)和式(2-2.15~2.17)構(gòu)成的方程組，可以得到N+3個di，通過式(2-2.11)可得任意坐標(biāo)下的Q函數(shù)。這就是曲面樣條函數(shù)法求解偏微分方程的基本原理。ikk注意到式(2-2.12)中的R2Qi們很容易將上述曲面樣條函數(shù)法求解偏微分方程的方法推廣至空間。設(shè)混合物為Kikk組分系統(tǒng)，這時，Q是定義在[xk0,(k=1,…,K-1);

K

1]x=(x1,…,xK-1,1)T，則(x)若該區(qū)域中N個點上的Q函數(shù)值Qj已知則按照上面的討論，通過這N個點的曲面樣條函數(shù)可以表示為NQ(x)[dR2ln(R2)]gTN

i ii iKR2 (x

)2(xx)T(xx

ik

k N NgN N這樣式(2-2.23)的含義與式(2-2.11)的含義是一致的。Ndi

iiiNdxiii

0KdQ N

2di(xxi)

)

dx

i

R2 xK1dQ

NNQ(x)[dR2ln(R2)]gTN

i

i KR2

)2

x)T(xx

k(jk

k

dQ N

2d

x) ln(R2)g

dx

R2

i1

聯(lián)立求解由式(2-2.16)和式(2-2.29)N+KdiQ滿足偏微分方程式(2-2.21)，將曲面樣條函數(shù)式(2-2.23)NDx

聯(lián)立求解式(2-2.16)和式(2-2.32)組成的方程組，同樣可以得到N+K個di，從而獲得偏微分方(d參數(shù))T、p、xy提供了極大方便。1,2,…,K組成的Kn組(T,p,x)數(shù)據(jù)，連KTpN=n+KQ與活度因子的關(guān)系，

pcalc

xk)(

k

k其中，ek=(kj)kj,kj=0k=j,kj=1HE、VE和Poynting 了。對前n組(T,p,x)數(shù)據(jù)，以系統(tǒng)總壓的計算值與實驗值之差為偏差函數(shù)，F(xiàn)j(pcalc,jpexp,j)/pexp,j，j=1,…, FQ函數(shù)的一個非線性偏微分方程。將曲面樣條函數(shù)式(2-2.23)代入式(2-2.34)F就轉(zhuǎn)化為N+K=n+2K個di的非線性方程組。對于K個純物質(zhì)，應(yīng)滿足如下條件，Q(xek)0，k=1,…, KFnkQ(ek)，k=1,…, ndiK個條件式(2-2.26)，N+1~N+K的偏差函數(shù)取為NiNii

dx(k)，k=1,…, ixi(k)xiki這樣，由式(2-2.34、2.36、2.37)組成了有n+2KN+K=n+2Kdiykpk/pcalc，k=1,…, 可得各組(Tp,x)y=(y1,…,yK)TF(F,...,

NX(d,...,

NN+KF(X)

這里的0為有N+K個元素的零向量。對式(2-2.41的求解可采用Broyden改進(jìn)的Newton-Raphson法(2-1)，這是一個非常有效的求解非線性方程組的方法[2-20]直接Gibbs-DuhemT、p、xy的方法[2-14,16]。對于一個含有K個組分的系統(tǒng)，按式(1-6.11)，其液相逸度的Gibbs-Duhem方程為，KxHK

H f

Vxkdln k1 dTmdp

k

RT kkKK xHKK

H 1ln(f

pθ)

m1T

dxjk

k

V(L)K1p

x[j,K

2 j1xjx[j,K2

m

dxjRTj1xj

x[j,K

KxHK

HK1

ln(f(L)/po)

T

V(L)p xk

i

m

dxj

j1k

RT

x

RTx

x[j,K

jx[j,K

jx[j,K]K-1個dxjK xHK

H(L) K ln(fK

po)

mT

V(L)pxk

k

m

k

x[j,K

xjx[j,K

RTxjx[j,K氣液平衡時，f(L)f(V)py k ln(p/po)

lny

lnxk

xk k

xk kk

x[j,K

k

x[j,K

k

x[j,KKxHK

H

mT

V(L)pi

m RT

jx[j,K

RT

jx[j,K lny

K1

xyxk k

K k

k

k1

x x[j,K x[j,KK ln(p/po) 1K xk

k

pxj x[j,K x[j,KjKH(L)xH*(L)HK

k1H H*(L)

KV(L)xV*(L)VK

km,k k1Lk和V*L)分別是純液體kHE和VE 過量摩爾體積，后者很小，一般情況下可以忽略不計。代入式(2-2.46)FjK1

xKyk

lnk Fj

xk k1

TK jx[jT

k

KxV*(L)KxV*(L)p

，j=1,2,…,K- HEKxHEKx

kk

km,k

RT

jx[j,K

jx[j,Kp因為kT、pyplnk

yi

lnk

T

lnk

p

i1

j j

x

yxj

x x[j,K

y[i,K]

x[j,K

x[j,K

x[j,K代入式(2-2.52)K1

xKyk

K1lnk

yi Fj

xk

k1

K

k

i1

j j

xx[j,K

y[i,K]

x[j,K

kkk x

k HEKxKHEKxKT

x k

y

jx[j,K

KxV*(

p kkm,kx k

k

x k

y

jx[j,KikK1

ln

yjFjkj

i

kk1

i

yk

xHE

kkk

y[k,K]KTKT

x[j,KKxKx

k

x

k

y

jx[j,K KxV*(

p kkm,kx k

kk

k

y

jx[j,KK-1yk、T、pT、p、xLk、V*(L)HE，k F

1

ln2 1 HExL(1x)

2x 1(1x) 2 RT

dp m,2x 1(1x) 2

式(2-2.56) HE

ln 2x

x) 2

RTxV

dp

m

x

x) 2

1

ln2

dx 1

p-xT-xyx的一個非線性常微分方程，可以用任何滿足精度要求的數(shù)值積分方法求解，最常用的是Runge-Kutta法，參見附錄2-2[2-14]。積分求解過程中，初始值可以取(x=0,y=0)x=0x=1(x=1,y=1)x=1x=0性，則恒沸點時有(dp/dx)=0或(dT/dx)=0，(dy/dx)為不定值，通常的數(shù)值積分方法難以越過此一個是，隨著積分的進(jìn)行誤差不斷累積，可能出現(xiàn)x=1時y1的情況。距x=M-1，xk=ix，i為格點序碼。將i點處的函數(shù)Fi寫成差分形式，

1

ln2yi1yiFi

(1

1

i

H

ln

dT 2,iix 1

2 i i

i

dx

ixV*(L)(1ix)V

ln

dp m,2ix 1(1ix) 2 p p

dxF(F,F,,

My(y,y,,

MM-1個式(2-2.58)F(y)

Newton-Raphsony(k)(J(k))1F(y(k)

y(k1)y(k)y(k

k為迭代序號，JJacobiFy的偏導(dǎo)數(shù)矩陣，是(M-1)M-1)階方陣。由式(2-2.58)可見Fi僅與yi-1、yi和yi+1有關(guān)，因此J是一個三對角線矩陣，式(2-2.62)可以方便地使用消去法，不必求逆，使工作量大為減輕。為了保證收斂，式(2-2.63)中y可乘一小于1的阻尼因子。在計算式(2-2.58)時，其中(dp/dx)或(dT/dx)T、p、x數(shù)據(jù)用樣條函數(shù)或多項式擬合后求導(dǎo)而得，j的計算可使用合適的狀態(tài)方程，y用上次迭代值。mHE一點和間接法有很大區(qū)別。具體例子和討論參見文獻(xiàn)[221、22、23]或[214]。對于三元系，式(22.58)可以寫出兩個式子仍可采用有限分法但計算效率較低收斂具體算見文獻(xiàn)[221]。mT、p、x推算y計算實T、p、xy的計算過程中，需要有計算氣相逸度因子k3章中介紹的合適的狀態(tài)方程，以下介紹的計算實例均選用了改進(jìn)的Peng-Robinson方程[2-25,26]。在曲面樣條函數(shù)法計算中，didi=0.001，對一般系統(tǒng)而言，Q函數(shù)是Q2-2.155℃下氯仿－乙醇二元系的計算結(jié)果與實驗值的比較。從表中結(jié)性和忽略液體體積，氣相組成的推算誤差從0.0016增加到0.0038。22.2101.325kPa下苯－異丙醇二元系的計算結(jié)果與實驗值的比較。從表中結(jié)果可見，有限差分直接法、三次樣條直接法和曲面樣條間接法推算的結(jié)果具有相同的準(zhǔn)確度，相互之間能很好的吻合，但間接法具有不需要過量焓數(shù)據(jù)的優(yōu)點。以三次樣條直接法為例，0.00760.01060.01050.01200.008升高到0.0181xx表2-2.2苯－異丙醇二元系推算結(jié)果比較(101.325kPa) 2-2.374.84℃124.98℃2-2.3二氯甲烷－硝基甲烷二元系推算結(jié)果比t=t=x2-2.4是苯(1)－正己烷(2)－環(huán)己烷(3)70℃表2-2.5是(1)－甲醇(2)－乙醇(3)三元系101.325kPa恒壓數(shù)據(jù)的推算結(jié)果逸度系數(shù)2-2.4苯(1)－正己烷(2)－環(huán)己烷(3)三元系推算結(jié)果比較表2-2.5(1)－甲醇(2)－乙醇(3)三元系推算結(jié)果比較 圖23和24畫出了用曲面樣條函數(shù)法得到的二氯甲烷(1)－氯仿(2)(3)三元系在45時，和(1)－氯仿(2)－乙(3)在11.325kPa下的Q函數(shù)曲面。前者比較簡單，三個二元系和三元系都表現(xiàn)出正偏差。后者相當(dāng)復(fù)雜，其中氯仿－乙醇和－乙醇兩個二元系是正偏差氯仿－乙醇二元系還有最低恒點而－氯仿二元系則是一個負(fù)偏差系統(tǒng)并有一個最高恒沸點，三元系的Q函數(shù)曲面則呈現(xiàn)復(fù)雜多變的形狀。對于這兩個三元系，推0.014后者也只有9次。這應(yīng)該說是對這一方法嚴(yán)峻的考驗。的計算實例可以參考文獻(xiàn)[224]。圖2-3二氯甲烷(1)-氯仿(2)-四氯化碳(3) 圖2-4(1)-氯仿(2)-乙醇(3)三元系三元系在45℃下的Q函數(shù)曲面 在101.325kPa下的Q函數(shù)曲面酮(4)在101.325Pa下的恒壓數(shù)據(jù)，推算的氣相組成誤差分別為y1=0.0133、y2=0.0134y3=0.0133和y3=0.0095，迭代共36次。對乙醇(1)－氯仿(2)－(3)－正己烷(4)在55℃下迭代共25次。詳細(xì)結(jié)果參見文獻(xiàn)[2-18]和[2-24]，結(jié)果令人滿意。當(dāng)實驗得到的T、p、xpxp*xp*xp*1 2 3

123(D1D2x1

nnCkk

x

jp*是純組分k的飽和蒸氣壓，對于恒壓系統(tǒng)，以系統(tǒng)壓力下純組分的沸點T*p*j Ck,ij由相應(yīng)的二元系(i-j)氣液平衡數(shù)據(jù)的熱力學(xué)一致性檢汽液平衡時有T、p、x和y四個變量，相律告訴我們其中只有兩個是獨立變量，如果實p、x數(shù)據(jù)推算y就是這種情況。另一方面，汽液平衡時的四個變量都是實驗可以直接測定得Gibbs-Duhem方程。根據(jù)式(1-7.41)，對于二元系，活度因子的Gibbs-Duhem方程可表示為

H

VTm

Qx1dln1,Ix2dln2,I

x112xdln1,Ixdln2,I12

取不同組成下曲線的斜率，然后代入式(2-3.3)，看是否滿足Gibbs-Duhem方程，這種檢驗方性或半定量的方法應(yīng)用。例如，在給定組成下，如果dln1,I/dx1是正值，那么dln2,I/dx2必須是負(fù)值；如果dln1,I/dx1等于零，則dln2,I/dx2也必須等于零。因此，斜率法能方便地用來檢2.2Tpxy的方法也可用來檢驗氣液相平衡實驗數(shù)據(jù)的熱力學(xué)一致性。如vanNess等[2-27]提出，由T、p、x推算得到的y與實驗值比較，如果y0.01滿足，即通也可以用y的推算值與實驗值的平均誤差來判斷數(shù)據(jù)的質(zhì)量，例如y0.01。這一方法還可以應(yīng)用于多元系。但以y0.01作為熱力學(xué)一致性的判斷標(biāo)準(zhǔn)，無法反映y的離散情況及yHerington[2-28]及Redlich和Kister[2-29]發(fā)展起來的面積檢驗法。Q與活度因子的關(guān)系，式(2-2.8)Qx1ln1,Ix2ln

x1dQ

xdln2,Iln

11

將式(2-3.3)dQln(

1,I/2,I

x1x1=0x1=1Q(x1=0)=0Q(x1=1)=0x11dQdx

x11ln1dx

x x

ln(1,I/2,I)x1x1=0~1x1軸所包面積應(yīng)等于零。典2-5x1軸所包面積等于零，這就意味著，如果根據(jù)某組氣液平衡數(shù)據(jù)計算得到各組分的活度因子，并以ln(1,I/2,I)x1x1軸上方的面積(面積A)x1軸下方的面積(面積B)時，該組數(shù)據(jù)是滿足熱力學(xué)一致性要求的。圖x1x1

式(23.6)是在恒溫、恒壓條件下推導(dǎo)得到的。但實際系統(tǒng)的氣液平衡都是在恒溫或者恒Q(17.41)代入式(23.5)，得 HE VEln(1,I/2,I) m

RT2 RT111x1x1=0x1=1x1

x11HE VEdpx0

1dx1

x0

mmRTdx1mm對于恒溫數(shù)據(jù)，HE項，通常VE很小可以忽略不計，式(2-3.8)仍可應(yīng)用。如為恒壓數(shù)據(jù) VE項，而HE通常是不能忽略的。但通常情況下實驗條件下的過量焓數(shù)據(jù)難以獲得， 以采用下面的方法檢驗。先由式(2-3.8)DJ比較。JJ150/

,T1

(a)最高恒沸混合物 (b)最低高恒沸混合物圖2-6式(2-3.11)中的計算

*L

*L

2

1

Poynting因子和氣相逸度因子對活度因子的影響很小，所以面積檢驗法對壓力測純組分的蒸氣壓之比p*/p*。對于恒溫數(shù)據(jù)，面積檢驗法幾乎只是判明蒸氣壓之比p*/p* x-yln(1,I/2,I)x1x-y數(shù)據(jù)的離散性也許是很敏感的，我們在2.2節(jié)已經(jīng)，任何嚴(yán)格符合熱力學(xué)一致性的一組N個結(jié)點上的T、p、x、yGibbs-Duhem方程。但是，從實驗誤差原理分析，由于實驗設(shè)備、物料純度以法就是在這樣的思想下建立起來的。它首先由Mcdermott和Ellis提出[2-31]，Ulichson和如果忽略溫度壓力對Q函數(shù)的影響，由T、p、x、y子應(yīng)滿足式(2-3.2)[xixi+1]間對式(2-3.2)F1

x

)1(2

x)(ln

i

1,i

i

i

如果實驗數(shù)據(jù)不是過分稀疏，則對于一組嚴(yán)格符合熱力學(xué)一致性的T、p、x、y實驗數(shù)據(jù)，上式積分的結(jié)果應(yīng)為零。由于隨機(jī)誤差的存在，實際上Fi也是一個隨量，其數(shù)學(xué)期望應(yīng)為零，且成正態(tài)分布。根據(jù)統(tǒng)計學(xué)理論，其兩次型服從自由度為N的2分布，即FQFT(σ2)1F~2(N)F

F其中σ2F的協(xié)方差矩陣，可以根據(jù)誤差傳遞原理根據(jù)式(2-3.13)T、p、x、y實驗數(shù)據(jù)的標(biāo)準(zhǔn)誤差σ2、σ2、σ2和σ2(由實驗者提供)計算得到。F 2-76的22-76的2分布概率密度曲線，22-7的2262

(6的概率為，2

的概率為(1-)P[22(N1 2

(NFQN的2分布，Q0~2-7可見，Q(即2)稱為置信度(或顯著性水平)1-=0.952(NP[Q

(N)]

Q

(N還要大的概率為5Q

(N小的概率為1-即95FTpxyF及其協(xié)方差矩陣2F并計算Q值，然后和查得的

(N比較。如果Q

(N，則在置信度為驗數(shù)據(jù)滿足熱力學(xué)一致性；如果Q2(N)，落入域，則在置信度為過校驗。0.0595Q中既包含了測量變量的隨機(jī)誤差，也包含了它們的系統(tǒng)誤差(如果存在的話)，單靠Q統(tǒng)計Dohnal和Fenclova[2-33]提出將Fii為單數(shù)(A)和i為雙數(shù)(B)的兩組。這樣，每組中的FiFiFwFwiiiiiFMF/S2(F)/M1/S2(wiiiii

S2(F) M(FF

)2/S2(F)/M1/S2(F)

wA,

M

i

i twFitwFw/S(Fw)~t(M

hwFi/S(Fi)wi ihM(FF)2/S2(F)~2(Mwi i

定的該組T、p、x、y實驗數(shù)據(jù)不存在系統(tǒng)偏差，隨機(jī)誤差在實驗者所估計的誤差范圍內(nèi)。Dohnal和Fenclova[2-33]tw對測量變量是否存在系統(tǒng)誤差很敏感，而對隨機(jī)誤差很不敏感，因此它非常適合于用來檢驗系統(tǒng)誤差；而統(tǒng)計量Qhw統(tǒng)計檢驗法是從逸度表達(dá)的Gibbs-Duhem方程出發(fā)建立起來的。對于一組實驗測定得到的二元系T、p、x、y數(shù)據(jù)(共有N個點)，如果它們是嚴(yán)格符合熱力學(xué)一致性的，則必須符合式(2-2.56)所示的Gibbs-Duhem方程。但由于隨機(jī)誤差的存在，T、p、x、y應(yīng)視作隨量，

1

ln2dyFi

1

xi

(1xi

dx i H

x

(1x)

ln

dTi i 2,ix 1i

(1x) 2 i i

i

i

i

dx

xV*(L)(1x)V ln

lndpii i m,2,ix 1ii

(1x) 2

i

dx

dy

dT

dp

*(

*(L)FiFTi,pi,xi,yi,dx,dx,dx,Hm,i,L1,i,L2,i,Vm,1,i,

FF,F

i

i TT,T,...,T pp,p,...,p xx,x,...,x yy,y,...,y dTdT dTTx , ,..., 1 NTdpdp dpTpx , ,...,

1 Ndydy

dyTyx , ,..., T

dx1dx

dxNHH

,,

,...,HE

m,Lm,

,

,

,

,...,L1,

,L2,

,

VV*(L),V*(L),V*(L),V*(L),...,V*(L),V*(L)

m,2,則式(2-3.21)

FF的協(xié)方差矩陣σ2Fσ2Fσ2FTFσ2FTFσ2FTFσ2FTFσ2FTFσ2 TT pp xx yy

yx

TxTx

Fσ2FTFσ2FTFσ2FTFσ2pxpx HH LL VVTpxyHσ2T的協(xié)方差矩陣，F(xiàn)TF對TNN階對角陣；σ2是p的協(xié)方差矩陣，F(xiàn)p是F對p的偏導(dǎo)數(shù)，均為NN階對角陣；σ2是x的協(xié)方差矩陣，F(xiàn)x是F對x的偏導(dǎo)數(shù)，均為NN階對角陣；σ2是y的協(xié)方差矩陣，F(xiàn)y是F對y的偏導(dǎo)數(shù)，均為NN階對角陣；σ2H的協(xié)方差矩陣，F(xiàn)HFHNN階對角陣；TpxyHLσ2L的協(xié)方差矩陣，F(xiàn)LFLN2N12LVσ2V的協(xié)方差矩陣，F(xiàn)VFVN2N12Vσ2yxNN階方陣，F(xiàn)yxFyxNNσσ2TxNN階方陣，F(xiàn)TxFTxNNσσ2pxNN階方陣，F(xiàn)pxFpxNNσ上面我們是將(dy/dx)、(dp/dx)和(dT/dx)作為獨立的隨量來考慮的，但實際測量的T、p、x、yT、p、xy數(shù)據(jù)計算得了三次樣條函數(shù)法，效果更好，下面簡要介紹這一方法。根據(jù)三次樣條函數(shù)法[2-3式(C-4)]，y對x的一階導(dǎo)數(shù)可表示為，dy

1h

yi1yiMi1Mi

dx

2 hi=xi+1-xi，MiMi+1由式(C-6)yxf(y,

xx yy fσ2fT xx yy和fx是f對x的偏導(dǎo)數(shù)，fy是f對y的偏導(dǎo)數(shù)，均為NN階矩陣。 2和由于F可看作N維隨量，理論上它的數(shù)學(xué)期望是零向量?？梢圆捎秒S機(jī)現(xiàn)象的數(shù)學(xué)應(yīng)服從自由度為N的2分布，即，F(xiàn)QFT(σ2)1F~2(N)F

從統(tǒng)計學(xué)意義來說，檢驗一組氣液平衡數(shù)據(jù)的熱力學(xué)一致性，就是要在承認(rèn)測量變量存在一定誤差的情況下，檢驗實驗數(shù)據(jù)(樣本)即檢驗假設(shè)：H0:

以Q為檢驗統(tǒng)計量，則在H0成立的條件下有式(2-3.15)。所以，在置信水平下的域QQ

(N)

當(dāng)實驗數(shù)據(jù)存在系統(tǒng)誤差時，顯然也將引起F對零的系統(tǒng)偏離，這時，F(xiàn)的數(shù)學(xué)期望將Q有貢獻(xiàn)，其系統(tǒng)偏差也有影響。當(dāng)我們用上述方法確定H0:

tF/S(F)~t(N

可以用來檢驗F的數(shù)學(xué)期望是否為零，這里F是F的平均值F

F/S2(F)/N1/S2(F

iiiS2(F) iii

(FF)2/S2(F)/N1/S2(F

N

i

i S2(Fi)F的協(xié)方差矩陣σ2中對角線上的元素，即S2Fσ2h F的方差與根據(jù)指定的測量變量的方差由式(2-3.34)Fihi

(FF)2/S2(F)~2(N

Fii統(tǒng)計量Q、t和h，在置信水平的條件下，若下列三式同時得到滿足，則認(rèn)為它們是符合熱FiiQQ

(N)

tt1(N

hh

(N

Qht因為F和S(F)均為平均值隨機(jī)誤差的大小對計算結(jié)果影響應(yīng)該很小實際數(shù)據(jù)的計Qht通過了，則說明該組數(shù)據(jù)不存在系統(tǒng)誤差，但實驗的真正誤差比估計的要大，可以改變給定的各個測量變量的方差(增大方差)重新Qh對大部分?jǐn)?shù)據(jù)的誤差都是不了解的(各種數(shù)據(jù)集中也鮮有給出原始數(shù)據(jù)誤差的)，因此無法用如上所說的根據(jù)實際的標(biāo)準(zhǔn)方差進(jìn)行推斷。一種替代的辦法是根據(jù)現(xiàn)有實驗裝置的測量精度情況，為T、p、x和y規(guī)定若干個誤差等級，如表2-3.1所示，用上述方法檢驗實驗數(shù)I級表示實驗精度很高，隨級數(shù)增加，誤差依次增大，最后第V級誤差很大，實際上是一個不能被接受的等級，T、p、x、y的標(biāo)準(zhǔn)誤差、、、HELm,i( 壓系統(tǒng))和液體摩爾體積V(L)(恒溫系統(tǒng))TpymH 0.05HmH

HE的實驗數(shù)據(jù)，以零代替，則可以通過過量自由焓GEmmH mH

Lm,kAntoineL對于恒溫系統(tǒng)，液體摩爾體積V*(L)對F的影響很小，其標(biāo)準(zhǔn)誤差可

mV0.1VmV

2-3.1T、p、x、y

IVxy的誤差由分析儀器的精度決定，T決定于溫度計的精度，pT、p、xy本身相互間不是獨立的，實驗中一個變量的波動會引起其它量的變化，例如在對一組實際測定的T、p、x、y數(shù)據(jù)進(jìn)行熱力學(xué)一致性檢驗時，除了要承認(rèn)測量儀器T、p、xy間的不獨立性而引起的相關(guān)性誤差也應(yīng)該是容許的。2

y

202

02

02

T

2

p

202

02

02

T

x y y202 02

T

T202 02

p

x 2

T、p、yx的一階導(dǎo)數(shù)已由前面的三次樣條函數(shù)法得到，其余的一階導(dǎo)數(shù)可p dln

dlnp* py 1(1y) 2

T y dln

dlnp* y(1y)

2

T y

p T1/p

p

T FFN維正態(tài)分布的理論分布的解析表達(dá)式，為此我們采用隨機(jī)現(xiàn)象的數(shù)學(xué)模擬方法(MonteCarlo方法)來進(jìn)行可以稱為MonteCarloFN組T、p、x、y的復(fù)雜函數(shù)，而后者均服從正態(tài)分布。FF(T,p,x

...,(T,p,xy)N

FN維正態(tài)分布的，則它的二次型應(yīng)服從參數(shù)為N的2分布，即式(2-3.38)。顯y的分布中抽樣得N組T、p、x、y數(shù)據(jù)(隨機(jī)數(shù))。F根據(jù)抽樣得到的NT、p、x、y數(shù)據(jù)，由前述方法計算F、σ2并得到QF機(jī)數(shù)Q(T,p,xy)1,...,(T,p,xy)N重復(fù)步驟(1)和(2)MMQ1、Q2、…、QM以Q的樣本分布Sn(Q)近似Q的分布如果假設(shè)是真的則SM(Q)將能足夠地近Q。我們的目的是要通過檢驗Q的2分布特性來證明F服從N維正態(tài)分布。即檢驗假設(shè)：0H:F(Q)F(2,N0我們采用柯爾莫哥檢驗法檢驗[2-35]Q的定義域(0,)以Qn n

iQ中小于或等于(n+1)Q值范圍內(nèi)樣本的經(jīng)驗分布函數(shù)與理論分布函數(shù)F0(2,N)的最大偏差，即，nnDmaxnnmax

F(2Q,N0(nQ)F(2nQ,N00Q 當(dāng)H0成立時，統(tǒng)計 MDn的極限分布函數(shù)(1)Lexp(2L2u2

uQ(u

MD)

u對于給定的顯著性水平，當(dāng)Dn>Dn()時便否定H0。臨界值Dn()可查常用數(shù)理統(tǒng)計表i采用上述方法，我們選擇了Benzene-Acetonitrile系統(tǒng)在25C下的恒溫數(shù)據(jù)(由T、p、x數(shù)據(jù)根據(jù)Mixon有限差分法得到氣相組成y，N=19)作為無誤差的標(biāo)準(zhǔn)數(shù)據(jù)。指定T、p、x、y2-3.1中方差等級的第I級，在計算機(jī)上產(chǎn)生方差為2(i=T、p、x、y)的隨機(jī)數(shù)并分別加于T、p、x和y而得到一套含有誤差的模擬實驗數(shù)據(jù)，這樣的模擬數(shù)據(jù)共產(chǎn)生100QM=100Q的隨機(jī)樣本，并作出樣本的經(jīng)驗分布函數(shù)Sn(Q)，由前述方法得到的統(tǒng)計量Dn(n=100)為0.0824。圖2-8是經(jīng)驗分布函數(shù)Sn(Q)和理論分布函數(shù)(2分布)的比較，可見兩者是非常接近的。根據(jù)柯爾莫哥i檢驗法，在顯著性水平＝0.05下給出的域D100,0.05=0.13403，即，D100D100,0.05，所以，可以認(rèn)為原假設(shè)H0QN的2FN2-8F氯仿－乙醇在55C下的氣液平衡數(shù)據(jù)是一套精度比較高的數(shù)據(jù)，用上述統(tǒng)計檢驗法對該套數(shù)據(jù)進(jìn)行熱力學(xué)一致性檢驗的結(jié)果，Q、thI級上通過檢驗。各實驗點上的偏F2-9所示。但如果忽略氣相非理想性，檢驗的結(jié)果發(fā)生很大的變化，Q在第II級上通過，h則在第It則不能通過檢驗。不盡如此，對測量變量取不同tF2-10所示。因為忽略氣相非理想性實際上意味著在實驗數(shù)據(jù)中引入了系統(tǒng)偏差，上述結(jié)果說明，統(tǒng)計量t圖2-9氯仿-乙醇系統(tǒng)55C下的F函數(shù) 圖2-10忽略氣相非理想性時，氯仿-乙醇系55CFtQh都不能通過檢驗。相反，有些數(shù)據(jù)能在第I和IIQ和httQh在更低一級的水平上通過檢驗，表明系統(tǒng)誤差對Q比h有更大的影響。表2- 定的四個變量T、p、x和y，它們并不是完全獨立的，而受Gibbs-Duhem方程的限制。如果Gibbs-Duhem方程；如果只存在隨機(jī)p、xyT、p、x、y實驗數(shù)據(jù)是否符合Gibbs-Duhem方程，在多大的誤差范圍內(nèi)符合，它們是否存在系統(tǒng)誤差？是不正常現(xiàn)象。除非超過某些明顯的界限例如3，我們才能說它是一個壞點而予以舍棄。2-1非線性方程組的Broyden解法[2-BroydenNewton-Raphson算法，其出發(fā)點是建立相鄰兩次迭代間矩陣的逆矩陣的遞推關(guān)系，從而可減少甚至避免求Jacobi矩陣及其逆矩陣。非線性方程組一般可以采用Newton-Raphson迭代法求解，即Xk1XkXk[Jk]1F(Xk

其中Jk為Jacobi矩陣。這一方法在計算Jacobi矩陣Jk及其逆矩陣[Jk]-1時所花時間很大。BroydenNewton-Raphson算法比較好地解決了這一問題。根據(jù)這一方法，相鄰兩次迭代間Jacobi矩陣的逆矩陣存在如下關(guān)系，[Jk]1[Jk1]1

XX0F0

計算JacobiJ0H0(3)HFXk(4)wk，使F(XkwkXkF(Xk：wk=1，若滿足下式則轉(zhuǎn)入步驟(5)F(XkwkXk)1/2

F(Xk)1/

wk[(16)1/21]/

F(XkwkXk)/F(Xk

如果經(jīng)規(guī)定次數(shù)的搜索后，仍不能使歐幾里得范數(shù)減小，則返回步驟(2)Xk重新計算Jacobi矩陣中的各個偏導(dǎo)數(shù)。Xk1Xkwk

Fk1

檢驗Fk+1是否滿足收斂條件，若不滿足，則ZkFk1

Hk

(HkZkwkXk)(Xk)T(Xk)THk

然后返回步驟(2)2-2Runge-Kutta法積分求解[2-dy

f(x,

從某初始值(x0,y0)x1=x0+x時，y值的增量yy'xf(x0,y''xf(x0x/2,y0y'/y'''xf(x0x/2,y0y''/