高等化工热力学-2相平衡计算

ProIIAspen等，可以提供很完善的计算方法。本章简单介绍相平衡计算的基本原理，[2-6,7]，这部分内容在其他专著中介绍得相对相平衡计算——有模还应输入那些为表征系统所必需的性质。本节讨论的相平衡计算，主要采用模型来输入那些性质。kx()x()z，求分相后各相组成x()和x()。定义组分k在相和K(αβ为：kK(αβ)x(α)/

对于计算相平衡的问题，最切题的普遍热力学关系式即相平衡判据。按式(1-5.15)()()，f()

f()，k1,2,,

相采用统一的状态方程。以式(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如应用第III或第IVK(k和相均采用活度因子关联式。如均采用第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如应用第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的函数，V是T、p、y的函数。运算时，先代入y pTyk=Kkxkyk需要使用状态方程，它就是为表征系统所应输入的模型。在使用状态方程计算k时，例如使用式(1-6.45)，在积分下限和压缩因子Z中，包含有未知数体积V(V)或V(L)，也要用状态方Kf*(L)x(L)/(p(V))K

kk

kk k式中k,I是T、p、x的函数，迭代运算与式(2-1.16)相同。除了计算(V)需要状态方程外，计k(L)需要活度因子关联式，它们都是为表征系统所应输入的模型。此外还需要f*(L)，它是k 统温度压力的函数，可利用式(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*是纯物质k (2-1.10)，需要知道KHk与应改为V

2-1计输入的独立变量为y以及T或p，要求输出x以及p或T。根据x必须满足归一化条件xKxk1

1，以式(2-1.1)KKyi/Kii称为方程。具体计算和求解泡点方程时类似

图2-2压力和液相组成的计算框合物的总组成z，输出为x和y以及汽相分数，定义为n(V)/

有一点需要说明，按相律，F=K-2+2-0-0=K，但现在输入了K+1个变量，因此输出的不仅xzk(1)xkyk(1Kk

kxk

xk1KKzk/(1Kk)k

子关联式以表征系统。首先可假设x和y的初值，利用闪蒸方程解得后，代入式(2-1.22)得到x和y的新值，反复迭代直至收敛。液液平T、pz，输出分相后两液相的组成x()、x()以及相分率，n()/

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

)x(

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

x()1KKzk/(1Kk)k

得x()和x()。液固平以某固体组分1在液态溶剂(可以是纯溶剂或混合溶剂)理，有f(S)

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

f*(S) 1后者是纯固体的特性，并依赖于系统的温度和压力。对于f(L111的液体，因此通常采用第IIf(L)1

11力的函数。代入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和液态、固态的恒压热容

、

系统压力下纯组分1的为Tf,1，由于压力对影响较小，Tf,1可近似地用常压或Tt,1Tf,1下，fusGm,1(Tf,1)=0。而在系统温度压力下，fusGm,1f*(Lf*(S (T)G*(L)G*(S)RTln(f*(L)/f*(S)

由于[(GT

T]p

T2(fusGm,1/T)

fus

如设

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 性质无关。一般情况下，fusCpm,1很小，对溶解度计算的影响也不大，常常可以忽略。剂1，即fusHm,2(Tf,2) 1 fusCpm,2

lnx2

1ln

TT f,2 TT当组分1和组分2同时析出时，称为最低共，其温度和组成称为低共熔温度和低共气液相平衡计算——由T、p、x推算气液平衡数据最常规的获得方法是采用各种循环式平衡釜抽取气、液相样品进行相组成分析，并同时测定平衡时的温度T和压力p，这是一种直接测定法。一些学者曾对不同气液[2-6~2-9]Tpxy(状态方程或/和过量函数模型)，即可根据相平衡原理计算得到另两个变量，这就是前面所讨论的气液相平衡的有模型计算法。如果不是输入分子热力学模型而是通过实验测定第三个变量，则原则上第四个变量可由它们通过相平衡计算确定。由于计算过程中没有引入模型，称为气液平衡的无模型计算法。由于实际测定中气相取样和组成测定是最的，一般是恒定温度T下静态法测定一系列液相组成x对应的系统总压p(溶液的饱和蒸汽压)pT(溶液的沸点)，简单地说就是T、p、x推算所谓无模型法并不是有时仍需要采用合适的模型计算一部分非关键性的决定系统特征的性质，例如气相逸度因子等。总之，至少是大大减少了模型的使用。国内外学者已建T、p、xym(G)方程的基础上建立起来的。根据应用-D(16.1)表示的逸度的-DTpxyT、p、x推算y，根据Q与活度因子的关系(隐含了D方程)Q函数法(间接法2-2.1.1Qk汽液平衡时，按判据式(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的函数，式(2-2.3)QT、p即可求得yp*、*和V*L 关。k决定于气相组成y，可利用上次迭代的y值计算，但还需要使用合适的状态方程，从这个意义上说，T、p、x推算y并不是完全的无模型，但当压力不太高时，气相非理想性远理想性，甚至可以令k=1，也不致带入严重误差。至于HE和VE，后者很小，常可忽略， 式(2-2.3)原则可以求解，但实践上却有很大，因为导数出现在exp中，是一个NewtonQ值。这种方法不依赖于任何过量函数模型，是严格的在一个框架下。大量实例计算表明，没有收敛的，不受多元系Q函数曲面类型的限制。下面介绍计算简单的Barker法和适用于多元系的曲面样条函数法，至于其它方法，感兴2-2.1.2BarkerBarker法[2-10,11]的是选择一个合适的过量函数模型，其中包括若干待定模型参数，代入式(2-2.3)T、p、xQ函数。4Barker法的候选模型。但Barker法的准Redlich-KisterQ12jQxx12j

jAj

x)j

1212

x2

jAj

x)j2xx

jjAj

x)

2jj2jj

x2

1212jA1212j

x)j2xx

1jjA1j

x)

221j2112j21T、p、1j2112j21采用最小二乘法关联得到N+1个Aj。具体计算时，我们可以采用如下的迭代过程：k由状态方程计算气相逸度因子*和k，由式(2-2.7)计算各组分的液相活度系数k py

k

k kk

由式(2-2.8)KKQxklnkQ关联式(2-2.4)Aj(j=0,…,由式(2-2.9)

yp**x

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

kk

k

比较y1和y0，如果两者不相等，则令y0=y1，转步骤(2)，进行新一轮循环迭代，直2-2.1.3QQx1x2Q(x1,x2)，Q-x1-x2形成三维空间，对于液相完全互溶的系统，Q函数是一个定义在[x10;x20;x1+x21]三角形区域内的Q=0QNQ函数值已Qi(x1(i),x2)Qi，可将Q看成是一块无限大平板纯弯曲时的变形[2-19]，Q上的负载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是未定系数，应满足以下三个条件，NNG1x1(i)dii

NNG2x2(i)dii

NNG3dii

NQdR2ln(R2) N

ji

i

N11(j

N22(j

N R2 )2

1(j 2(j

kjj的弹性常数，cj是对曲面光顺要求所加的权重，在某点上该值取为零，则所构作的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)构成的方程组，可以得到N+3个di，通过式(2-2.11)可得任意坐标下的Q函数。这就是曲面样条函数法求解偏微分方程的基本原理。ikk注意到式(2-2.12)中的R2Qi们很容易将上述曲面样条函数法求解偏微分方程的方法推广至空间。设混合物为Kikk组分系统，这时，Q是定义在[xk0,(k=1,…,K-1);

K

1]x=(x1,…,xK-1,1)T，则(x)若该区域中N个点上的Q函数值Qj已知则按照上面的讨论，通过这N个点的曲面样条函数可以表示为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

联立求解由式(2-2.16)和式(2-2.29)N+KdiQ满足偏微分方程式(2-2.21)，将曲面样条函数式(2-2.23)NDx

联立求解式(2-2.16)和式(2-2.32)组成的方程组，同样可以得到N+K个di，从而获得偏微分方(d参数)T、p、xy提供了极大方便。1,2,…,K组成的Kn组(T,p,x)数据，连KTpN=n+KQ与活度因子的关系，

pcalc

xk)(

k

k其中，ek=(kj)kj,kj=0k=j,kj=1HE、VE和Poynting 了。对前n组(T,p,x)数据，以系统总压的计算值与实验值之差为偏差函数，Fj(pcalc,jpexp,j)/pexp,j，j=1,…, FQ函数的一个非线性偏微分方程。将曲面样条函数式(2-2.23)代入式(2-2.34)F就转化为N+K=n+2K个di的非线性方程组。对于K个纯物质，应满足如下条件，Q(xek)0，k=1,…, KFnkQ(ek)，k=1,…, ndiK个条件式(2-2.26)，N+1~N+K的偏差函数取为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改进的Newton-Raphson法(2-1)，这是一个非常有效的求解非线性方程组的方法[2-20]直接Gibbs-DuhemT、p、xy的方法[2-14,16]。对于一个含有K个组分的系统，按式(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的一个非线性常微分方程，可以用任何满足精度要求的数值积分方法求解，最常用的是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)为不定值，通常的数值积分方法难以越过此一个是，随着积分的进行误差不断累积，可能出现x=1时y1的情况。距x=M-1，xk=ix，i为格点序码。将i点处的函数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的偏导数矩阵，是(M-1)M-1)阶方阵。由式(2-2.58)可见Fi仅与yi-1、yi和yi+1有关，因此J是一个三对角线矩阵，式(2-2.62)可以方便地使用消去法，不必求逆，使工作量大为减轻。为了保证收敛，式(2-2.63)中y可乘一小于1的阻尼因子。在计算式(2-2.58)时，其中(dp/dx)或(dT/dx)T、p、x数据用样条函数或多项式拟合后求导而得，j的计算可使用合适的状态方程，y用上次迭代值。mHE一点和间接法有很大区别。具体例子和讨论参见文献[221、22、23]或[214]。对于三元系，式(22.58)可以写出两个式子仍可采用有限分法但计算效率较低收敛具体算见文献[221]。mT、p、x推算y计算实T、p、xy的计算过程中，需要有计算气相逸度因子k3章中介绍的合适的状态方程，以下介绍的计算实例均选用了改进的Peng-Robinson方程[2-25,26]。在曲面样条函数法计算中，didi=0.001，对一般系统而言，Q函数是Q2-2.155℃下氯仿－乙醇二元系的计算结果与实验值的比较。从表中结性和忽略液体体积，气相组成的推算误差从0.0016增加到0.0038。22.2101.325kPa下苯－异丙醇二元系的计算结果与实验值的比较。从表中结果可见，有限差分直接法、三次样条直接法和曲面样条间接法推算的结果具有相同的准确度，相互之间能很好的吻合，但间接法具有不需要过量焓数据的优点。以三次样条直接法为例，0.00760.01060.01050.01200.008升高到0.0181xx表2-2.2苯－异丙醇二元系推算结果比较(101.325kPa) 2-2.374.84℃124.98℃2-2.3二氯甲烷－硝基甲烷二元系推算结果比t=t=x2-2.4是苯(1)－正己烷(2)－环己烷(3)70℃表2-2.5是(1)－甲醇(2)－乙醇(3)三元系101.325kPa恒压数据的推算结果逸度系数2-2.4苯(1)－正己烷(2)－环己烷(3)三元系推算结果比较表2-2.5(1)－甲醇(2)－乙醇(3)三元系推算结果比较 图23和24画出了用曲面样条函数法得到的二氯甲烷(1)－氯仿(2)(3)三元系在45时，和(1)－氯仿(2)－乙(3)在11.325kPa下的Q函数曲面。前者比较简单，三个二元系和三元系都表现出正偏差。后者相当复杂，其中氯仿－乙醇和－乙醇两个二元系是正偏差氯仿－乙醇二元系还有最低恒点而－氯仿二元系则是一个负偏差系统并有一个最高恒沸点，三元系的Q函数曲面则呈现复杂多变的形状。对于这两个三元系，推0.014后者也只有9次。这应该说是对这一方法严峻的考验。的计算实例可以参考文献[224]。图2-3二氯甲烷(1)-氯仿(2)-四氯化碳(3) 图2-4(1)-氯仿(2)-乙醇(3)三元系三元系在45℃下的Q函数曲面 在101.325kPa下的Q函数曲面酮(4)在101.325Pa下的恒压数据，推算的气相组成误差分别为y1=0.0133、y2=0.0134y3=0.0133和y3=0.0095，迭代共36次。对乙醇(1)－氯仿(2)－(3)－正己烷(4)在55℃下迭代共25次。详细结果参见文献[2-18]和[2-24]，结果令人满意。当实验得到的T、p、xpxp*xp*xp*1 2 3

123(D1D2x1

nnCkk

x

jp*是纯组分k的饱和蒸气压，对于恒压系统，以系统压力下纯组分的沸点T*p*j Ck,ij由相应的二元系(i-j)气液平衡数据的热力学一致性检汽液平衡时有T、p、x和y四个变量，相律告诉我们其中只有两个是独立变量，如果实p、x数据推算y就是这种情况。另一方面，汽液平衡时的四个变量都是实验可以直接测定得Gibbs-Duhem方程。根据式(1-7.41)，对于二元系，活度因子的Gibbs-Duhem方程可表示为

H

VTm

Qx1dln1,Ix2dln2,I

x112xdln1,Ixdln2,I12

取不同组成下曲线的斜率，然后代入式(2-3.3)，看是否满足Gibbs-Duhem方程，这种检验方性或半定量的方法应用。例如，在给定组成下，如果dln1,I/dx1是正值，那么dln2,I/dx2必须是负值；如果dln1,I/dx1等于零，则dln2,I/dx2也必须等于零。因此，斜率法能方便地用来检2.2Tpxy的方法也可用来检验气液相平衡实验数据的热力学一致性。如vanNess等[2-27]提出，由T、p、x推算得到的y与实验值比较，如果y0.01满足，即通也可以用y的推算值与实验值的平均误差来判断数据的质量，例如y0.01。这一方法还可以应用于多元系。但以y0.01作为热力学一致性的判断标准，无法反映y的离散情况及yHerington[2-28]及Redlich和Kister[2-29]发展起来的面积检验法。Q与活度因子的关系，式(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轴所包面积应等于零。典2-5x1轴所包面积等于零，这就意味着，如果根据某组气液平衡数据计算得到各组分的活度因子，并以ln(1,I/2,I)x1x1轴上方的面积(面积A)x1轴下方的面积(面积B)时，该组数据是满足热力学一致性要求的。图x1x1

式(23.6)是在恒温、恒压条件下推导得到的。但实际系统的气液平衡都是在恒温或者恒Q(17.41)代入式(23.5)，得 HE VEln(1,I/2,I) m

RT2 RT111x1x1=0x1=1x1

x11HE VEdpx0

1dx1

x0

mmRTdx1mm对于恒温数据，HE项，通常VE很小可以忽略不计，式(2-3.8)仍可应用。如为恒压数据 VE项，而HE通常是不能忽略的。但通常情况下实验条件下的过量焓数据难以获得， 以采用下面的方法检验。先由式(2-3.8)DJ比较。JJ150/

,T1

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

*L

*L

2

1

Poynting因子和气相逸度因子对活度因子的影响很小，所以面积检验法对压力测纯组分的蒸气压之比p*/p*。对于恒温数据，面积检验法几乎只是判明蒸气压之比p*/p* x-yln(1,I/2,I)x1x-y数据的离散性也许是很敏感的，我们在2.2节已经，任何严格符合热力学一致性的一组N个结点上的T、p、x、yGibbs-Duhem方程。但是，从实验误差原理分析，由于实验设备、物料纯度以法就是在这样的思想下建立起来的。它首先由Mcdermott和Ellis提出[2-31]，Ulichson和如果忽略温度压力对Q函数的影响，由T、p、x、y子应满足式(2-3.2)[xixi+1]间对式(2-3.2)F1

x

)1(2

x)(ln

i

1,i

i

i

如果实验数据不是过分稀疏，则对于一组严格符合热力学一致性的T、p、x、y实验数据，上式积分的结果应为零。由于随机误差的存在，实际上Fi也是一个随量，其数学期望应为零，且成正态分布。根据统计学理论，其两次型服从自由度为N的2分布，即FQFT(σ2)1F~2(N)F

F其中σ2F的协方差矩阵，可以根据误差传递原理根据式(2-3.13)T、p、x、y实验数据的标准误差σ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及其协方差矩阵2F并计算Q值，然后和查得的

(N比较。如果Q

(N，则在置信度为验数据满足热力学一致性；如果Q2(N)，落入域，则在置信度为过校验。0.0595Q中既包含了测量变量的随机误差，也包含了它们的系统误差(如果存在的话)，单靠Q统计Dohnal和Fenclova[2-33]提出将Fii为单数(A)和i为双数(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实验数据不存在系统偏差，随机误差在实验者所估计的误差范围内。Dohnal和Fenclova[2-33]tw对测量变量是否存在系统误差很敏感，而对随机误差很不敏感，因此它非常适合于用来检验系统误差；而统计量Qhw统计检验法是从逸度表达的Gibbs-Duhem方程出发建立起来的。对于一组实验测定得到的二元系T、p、x、y数据(共有N个点)，如果它们是严格符合热力学一致性的，则必须符合式(2-2.56)所示的Gibbs-Duhem方程。但由于随机误差的存在，T、p、x、y应视作随量，

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的协方差矩阵σ2Fσ2Fσ2FTFσ2FTFσ2FTFσ2FTFσ2FTFσ2 TT pp xx yy

yx

TxTx

Fσ2FTFσ2FTFσ2FTFσ2pxpx HH LL VVTpxyHσ2T的协方差矩阵，FTF对TNN阶对角阵；σ2是p的协方差矩阵，Fp是F对p的偏导数，均为NN阶对角阵；σ2是x的协方差矩阵，Fx是F对x的偏导数，均为NN阶对角阵；σ2是y的协方差矩阵，Fy是F对y的偏导数，均为NN阶对角阵；σ2H的协方差矩阵，FHFHNN阶对角阵；TpxyHLσ2L的协方差矩阵，FLFLN2N12LVσ2V的协方差矩阵，FVFVN2N12Vσ2yxNN阶方阵，FyxFyxNNσσ2TxNN阶方阵，FTxFTxNNσσ2pxNN阶方阵，FpxFpxNNσ上面我们是将(dy/dx)、(dp/dx)和(dT/dx)作为独立的随量来考虑的，但实际测量的T、p、x、yT、p、xy数据计算得了三次样条函数法，效果更好，下面简要介绍这一方法。根据三次样条函数法[2-3式(C-4)]，y对x的一阶导数可表示为，dy

1h

yi1yiMi1Mi

dx

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

xx yy fσ2fT xx yy和fx是f对x的偏导数，fy是f对y的偏导数，均为NN阶矩阵。 2和由于F可看作N维随量，理论上它的数学期望是零向量。可以采用随机现象的数学应服从自由度为N的2分布，即，FQFT(σ2)1F~2(N)F

从统计学意义来说，检验一组气液平衡数据的热力学一致性，就是要在承认测量变量存在一定误差的情况下，检验实验数据(样本)即检验假设：H0:

以Q为检验统计量，则在H0成立的条件下有式(2-3.15)。所以，在置信水平下的域QQ

(N)

当实验数据存在系统误差时，显然也将引起F对零的系统偏离，这时，F的数学期望将Q有贡献，其系统偏差也有影响。当我们用上述方法确定H0:

tF/S(F)~t(N

可以用来检验F的数学期望是否为零，这里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的协方差矩阵σ2中对角线上的元素，即S2Fσ2h F的方差与根据指定的测量变量的方差由式(2-3.34)Fihi

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

Fii统计量Q、t和h，在置信水平的条件下，若下列三式同时得到满足，则认为它们是符合热FiiQQ

(N)

tt1(N

hh

(N

Qht因为F和S(F)均为平均值随机误差的大小对计算结果影响应该很小实际数据的计Qht通过了，则说明该组数据不存在系统误差，但实验的真正误差比估计的要大，可以改变给定的各个测量变量的方差(增大方差)重新Qh对大部分数据的误差都是不了解的(各种数据集中也鲜有给出原始数据误差的)，因此无法用如上所说的根据实际的标准方差进行推断。一种替代的办法是根据现有实验装置的测量精度情况，为T、p、x和y规定若干个误差等级，如表2-3.1所示，用上述方法检验实验数I级表示实验精度很高，随级数增加，误差依次增大，最后第V级误差很大，实际上是一个不能被接受的等级，T、p、x、y的标准误差、、、HELm,i( 压系统)和液体摩尔体积V(L)(恒温系统)TpymH 0.05HmH

HE的实验数据，以零代替，则可以通过过量自由焓GEmmH mH

Lm,kAntoineL对于恒温系统，液体摩尔体积V*(L)对F的影响很小，其标准误差可

mV0.1VmV

2-3.1T、p、x、y

IVxy的误差由分析仪器的精度决定，T决定于温度计的精度，pT、p、xy本身相互间不是独立的，实验中一个变量的波动会引起其它量的变化，例如在对一组实际测定的T、p、x、y数据进行热力学一致性检验时，除了要承认测量仪器T、p、xy间的不独立性而引起的相关性误差也应该是容许的。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的一阶导数已由前面的三次样条函数法得到，其余的一阶导数可p dln

dlnp* py 1(1y) 2

T y dln

dlnp* y(1y)

2

T y

p T1/p

p

T FFN维正态分布的理论分布的解析表达式，为此我们采用随机现象的数学模拟方法(MonteCarlo方法)来进行可以称为MonteCarloFN组T、p、x、y的复杂函数，而后者均服从正态分布。FF(T,p,x

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

FN维正态分布的，则它的二次型应服从参数为N的2分布，即式(2-3.38)。显y的分布中抽样得N组T、p、x、y数据(随机数)。F根据抽样得到的NT、p、x、y数据，由前述方法计算F、σ2并得到QF机数Q(T,p,xy)1,...,(T,p,xy)N重复步骤(1)和(2)MMQ1、Q2、…、QM以Q的样本分布Sn(Q)近似Q的分布如果假设是真的则SM(Q)将能足够地近Q。我们的目的是要通过检验Q的2分布特性来证明F服从N维正态分布。即检验假设：0H:F(Q)F(2,N0我们采用柯尔莫哥检验法检验[2-35]Q的定义域(0,)以Qn n

iQ中小于或等于(n+1)Q值范围内样本的经验分布函数与理论分布函数F0(2,N)的最大偏差，即，nnDmaxnnmax

F(2Q,N0(nQ)F(2nQ,N00Q 当H0成立时，统计 MDn的极限分布函数(1)Lexp(2L2u2

uQ(u

MD)

u对于给定的显著性水平，当Dn>Dn()时便否定H0。临界值Dn()可查常用数理统计表i采用上述方法，我们选择了Benzene-Acetonitrile系统在25C下的恒温数据(由T、p、x数据根据Mixon有限差分法得到气相组成y，N=19)作为无误差的标准数据。指定T、p、x、y2-3.1中方差等级的第I级，在计算机上产生方差为2(i=T、p、x、y)的随机数并分别加于T、p、x和y而得到一套含有误差的模拟实验数据，这样的模拟数据共产生100QM=100Q的随机样本，并作出样本的经验分布函数Sn(Q)，由前述方法得到的统计量Dn(n=100)为0.0824。图2-8是经验分布函数Sn(Q)和理论分布函数(2分布)的比较，可见两者是非常接近的。根据柯尔莫哥i检验法，在显著性水平＝0.05下给出的域D100,0.05=0.13403，即，D100D100,0.05，所以，可以认为原假设H0QN的2FN2-8F氯仿－乙醇在55C下的气液平衡数据是一套精度比较高的数据，用上述统计检验法对该套数据进行热力学一致性检验的结果，Q、thI级上通过检验。各实验点上的偏F2-9所示。但如果忽略气相非理想性，检验的结果发生很大的变化，Q在第II级上通过，h则在第It则不能通过检验。不尽如此，对测量变量取不同tF2-10所示。因为忽略气相非理想性实际上意味着在实验数据中引入了系统偏差，上述结果说明，统计量t图2-9氯仿-乙醇系统55C下的F函数 图2-10忽略气相非理想性时，氯仿-乙醇系55CFtQh都不能通过检验。相反，有些数据能在第I和IIQ和httQh在更低一级的水平上通过检验，表明系统误差对Q比h有更大的影响。表2- 定的四个变量T、p、x和y，它们并不是完全独立的，而受Gibbs-Duhem方程的限制。如果Gibbs-Duhem方程；如果只存在随机p、xyT、p、x、y实验数据是否符合Gibbs-Duhem方程，在多大的误差范围内符合，它们是否存在系统误差？是不正常现象。除非超过某些明显的界限例如3，我们才能说它是一个坏点而予以舍弃。2-1非线性方程组的Broyden解法[2-BroydenNewton-Raphson算法，其出发点是建立相邻两次迭代间矩阵的逆矩阵的递推关系，从而可减少甚至避免求Jacobi矩阵及其逆矩阵。非线性方程组一般可以采用Newton-Raphson迭代法求解，即Xk1XkXk[Jk]1F(Xk

其中Jk为Jacobi矩阵。这一方法在计算Jacobi矩阵Jk及其逆矩阵[Jk]-1时所花时间很大。BroydenNewton-Raphson算法比较好地解决了这一问题。根据这一方法，相邻两次迭代间Jacobi矩阵的逆矩阵存在如下关系，[Jk]1[Jk1]1

XX0F0

计算JacobiJ0H0(3)HFXk(4)wk，使F(XkwkXkF(Xk：wk=1，若满足下式则转入步骤(5)F(XkwkXk)1/2

F(Xk)1/

wk[(16)1/21]/

F(XkwkXk)/F(Xk

如果经规定次数的搜索后，仍不能使欧几里得范数减小，则返回步骤(2)Xk重新计算Jacobi矩阵中的各个偏导数。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''/