当前位置:首页 >> 理学 >>

2相平衡计算


2 相平衡计算
迄今已有很多专著介绍相平衡的计算方法,见文献[2-1~2~5]。一些大型过程模拟软件, 如 Pro II 和 Aspen 等,可以提供很完善的计算方法。本章简单介绍相平衡计算的基本原理, 至于具体的编程技巧等方面的细节可以参看上述专著。本章花比较多的篇幅介绍相平衡计算 的无模型法及其在气液平衡数据测定中的应用,以及气液平衡实验数据的热力学一致性检验
[2-6,7]

,这部分内容在其他专著中介绍得相对较少。

2.1 相平衡计算——有模型法
为解决一个具体的相平衡问题,在有了切题的普遍热力学关系式,并确定了独立变量后, 还应输入那些为表征系统所必需的性质。本节讨论的相平衡计算,主要采用模型来输入那些 性质。 相平衡问题往往表现为:已知一个相的组成 x(?),求另一相的组成 x(?);或已知总组成 z, 求分相后各相组成 x(?)和 x(?)。定义组分 k 在?相和?相中分配的相平衡常数 K k( αβ ) 为:
(α) (β ) Kk(αβ ) ? xk / xk

(2-1.1)

相平衡问题的中心,就是要计算每一组分的相平衡常数。 对于计算相平衡的问题,最切题的普遍热力学关系式即相平衡判据。按式 (1-5.15) 和 (1-6.13),

?k(? ) ? ?k( ? ) , fk(? ) ? fk( ? ) , k ? 1, 2, ?, K

(2-1.2)

在第一章中, 已介绍了两种计算逸度的方法, 即状态方程法和活度因子法, 具体应用于式(2-1.2) 时有三种选择: (1) ? 和 ? 相采用统一的状态方程。以式(1-6.21)代入式(2-1.2),
(? ) (? ) (? ) (? ) pxk ?k ? pxk ?k ( αβ ) (β ) (α) Kk ? ?k / ?k

(2-1.3) (2-1.4)

这种选择可用于气液和液液平衡的计算,特别是高压气液平衡的计算。 (2) ?相 (例如气相) 采用状态方程, ?相(例如液、 固相)采用活度因子关联式。 以式(1-6.21) 和(1-7.18)代入式(2-1.2),对于不同类型的活度或活度因子,可分别写出:
(? ) (? ) (? ) (? ) pxk ?k ? fk*( ? ) xk ? k ,I

(2-1.5)

(? ) (? ) (? ) ? KHk xk ? k ,II (? ) (? ) o (? ) ? KHk ( m) (mk / m )? k ,III (? ) (? ) o (? ) ? KHk ( c ) (ck / c )? k ,IV (?? ) (? ) (? ) Kk ? fk*( ? )? k ,I /( p?k ) (?? ) (? ) (? ) Kk ? KHk ? k ,II /( p?k(? ) )

(2-1.6) (2-1.7) (2-1.8) (2-1.9) (2-1.10)

如应用第 III 或第 IV 种活度, K k(?? ) 需重新定义或进行浓度单位换算。这种选择常用于气液平 衡和气固平衡。 (3) ?和?相均采用活度因子关联式。如均采用第 I 种活度,式(2-1.2)变为
(? ) (? ) (? ) (? ) fk* xk ? k ,I ? fk* xk ? k ,I (? ) (? ) Kk(?? ) ? ? k ,I / ? k ,I

(2-1.11) (2-1.12)

如?相用第 I 种活度,?相用第 II 种活度,式(2-1.2)变为
(? ) (? ) (? ) (? ) (? ) fk*(? ) xk ? k ,I ? KHk xk ? k ,II (?? ) (? ) (? ) ) Kk ? KHk ? k ,II /( f k*(? )? k(? ,I )

(2-1.13) (2-1.14)

如应用第 III 或第 IV 种活度, K k(?? ) 需重新定义或进行浓度单位换算。这种选择常用于液液平 衡和液固平衡。

2.1.1 气液平衡
(L) (V) 气液平衡的计算有三种, 即泡点计算、 露点计算和闪蒸计算。 下面以 xk ? xk ,yk ? xk 或 (g) ( VL ) , Kk ? Kk 分别代表液相组成、气相组成和气液平衡常数。 xk

2.1.1.1 泡点计算 目的是求液相在一定压力下的沸点即泡点,或一定温度下的蒸气压,以及平衡的气相组 成。输入的独立变量为 x 以及 T 或 p,要求输出 y 以及 p 或 T。最切题的普遍热力学关系式为

fk(V) ? fk(L) ,这是一个 K 维方程组。实践中,可简化为一元方程求根,中心是建立泡点方程。
根据 y 必需满足归一化条件 ?k ?1 yk ? 1 ,以式(2-1.1)代入,
K

?K x
k ?1

K

k k

?1

(2-1.15)

这就是泡点方程,式中 Kk 可采用上面介绍过的任一个式子计算。例如采用式(2-1.4),可得

? (?
k ?1

K

(L) k k

x / ?k( V ) ) ? 1

(2-1.16)

式中 ?i( L) 是 T、p、x 的函数, ?i( V) 是 T、p、y 的函数。运算时,先代入 y 的初值,方程中就 只含一个未知数,即 p 或 T,解得后根据 yk=Kkxk 计算 y 的新值,反复迭代直至收敛。为计算

?k 需要使用状态方程,它就是为表征系统所应输入的模型。在使用状态方程计算?k 时,例如
使用式(1-6.45),在积分下限和压缩因子 Z 中,包含有未知数体积 V(V)或 V(L),也要用状态方 程求取。又如采用式(2-1.9),

?f
k ?1

K

*( L ) (L) k k k ,I

x?

/( p?k( V ) ) ? 1

(2-1.17)

(V) 式中?k,I 是 T、p、x 的函数,迭代运算与式(2-1.16)相同。除了计算 ?k 需要状态方程外,计算

) *( L ) ,它是系 ? k(L , I 需要活度因子关联式,它们都是为表征系统所应输入的模型。此外还需要 f k

统温度压力的函数,可利用式(1-6.6)积分,
p * *( L ) ? ? f k*( L ) (T , p) ? f k* (T , pk ) exp?? Vm ? , k dp ? / RT ? * ? ? ? ?? p k ?

(2-1.18)

*( L ) * * * 式中 pk 和 Vm , k 是纯液体 k 的饱和蒸气压和摩尔体积, exp 项常称为 Poynting 因子。 f k (T , pk )

可按下式求得,
* * * * * fk*(L) (T , pk ) ? fk*(V) (T , pk ) ? pk ?k (T , pk )

(2-1.19)

* * 式中 ?k (T , pk ) 是纯物质 k 饱和蒸气的逸度因子,可用纯物质的状态方程计算。如果采用式

(2-1.10),需要知道 KHk,可直接采用实验数据,或采用估算用的半经验关联式,后者可参阅
* 文献[2-2]和[2-3]。 应该注意 KHk 与 f k* 类似, 也与压力有关, 换算式类似于式(2-1.18), 只是 Vm ,k ? 应改为 Vm , k ,即无限稀释下的偏摩尔体积。

图 2-1 泡点温度和气相组成的计算框图

2.1.1.2 露点计算 目的是求气相在一定压力下的露点,或一定温度即露点下的压力,以及平衡的液相组成。 输入的独立变量为 y 以及 T 或 p,要求输出 x 以及 p 或 T。根据 x 必须满足归一化条件

?
K i ?1

K k ?1 k

x ? 1 ,以式(2-1.1)代入,得

?y /K
i

i

?1

(2-1.20)

称为露点方程。具体计算和求解泡点方程时类似。

图 2-2 露点压力和液相组成的计算框图

2.1.1.3 闪蒸计算 目的是求一定温度压力下,混合物分相后的汽液两相组成。输入的独立变量为 T、p 和混 合物的总组成 z,输出为 x 和 y 以及汽相分数?,定义为

? ? n( V ) / n

(2-1.21)

有一点需要说明,按相律,F=K-2+2-0-0=K,但现在输入了 K+1 个变量,因此输出的不仅 x 和 y 等强度性质,还输出了一个另外的性质即?,经验表明,这样的计算效率较高。根据物料 衡算,并以式(2-1.1)代入,得

zk ? (1 ? ? ) xk ? ?yk ? (1 ? ? ? ?Kk ) xk
由于 x 必需归一, ?k ?1 xk ? 1 ,以上式代入得
K

(2-1.22)

?z
k ?1

K

k

/(1 ? ? ? ?K k ) ? 1

(2-1.23)

称为闪蒸方程。式中 Kk 可用本节中的任一个表达式计算,计算时需要输入状态方程或活度因 子关联式以表征系统。首先可假设 x 和 y 的初值,利用闪蒸方程解得?后,代入式(2-1.22)得 到 x 和 y 的新值,反复迭代直至收敛。

2.1.2 液液平衡
液液平衡的计算和闪蒸计算非常类似,输入 T、p 和混合物总组成 z,输出分相后两液相 的组成 x(?)、x(?)以及?相分率?,

? ? n( ? ) / n
(?? ) (? ) (? ) 将 Kk 简记为 Kk,根据物料衡算并以式(2-1.1)代入,得 ? xk / xk (? ) (? ) (? ) zk ? (1 ??) xk ? ?xk ? (1 ?? ? ?Kk ) xk
(? ) ? 1,以上式代入得 由于 x(?)必需归一, ? k ?1 xk K

(2-1.24)

(2-1.25)

?z
k ?1

K

k

/(1 ? ? ? ?K k ) ? 1

(2-1.26)

这就是计算液液平衡的方程,在输入状态方程或活度系数关联式以表征系统后,即可迭代求 得 x(?)和 x(?)。

2.1.3 液固平衡
以某固体组分 1 在液态溶剂(可以是纯溶剂或混合溶剂)中的溶解度为例,按照相平衡原 理,有 f1(S) ? f1(L) 。固相中的溶剂可略去,因此不必列出 f2 的等式,并且可写出 f1(S) ? f1*(S) ,

后者是纯固体的特性,并依赖于系统的温度和压力。对于 f1(L) ,由于系统温度压力下一般不 存在纯组分 1 的液体,因此通常采用第 II 种活度来表示, f1(L) ? KH1x1? 1,II ,KH1 也是温度和压 力的函数。代入 f1(S) ? f1(L) ,得

f1*(S) ? KH1 x1?1,II

(2-1.27)

式中?1,II 的计算需要活度因子关联式,连同 f1*(S ) 和 KH1,是必需输入以表征系统的性质。独立 变量为 T 和 p,输入后即可输出溶解度,即 x1。 输入 f1*(S ) 和 KH1 并不是唯一的选择,另一种更实用的方法是:对 f1(L) ,采用第一种活度,

f1(L) ? f1*(L) x1?1,I , f1*( L ) 可看作是系统温度压力下过冷液体 1 的逸度。代入 f1(S) ? f1(L) ,得 f1*(S) ? f1*(L) x1? 1,I
(2-1.28)

(L) (S ) 比值 f1*(L) / f1*(S) 可由组分 1 的熔化热?fusHm,1 和液态、固态的恒压热容 C pm ,1 、 C pm ,1 求得。设在

系统压力下纯组分 1 的熔点为 Tf,1,由于压力对熔点影响较小,Tf,1 可近似地用常压熔点或三 相点 Tt,1 代替。在 Tf,1 下,?fusGm,1(Tf,1)=0。而在系统温度压力下,?fusGm,1 与 f1*(L) / f1*(S) 通过下 式相联系,
*(L) *(S) *(L) ?fus Gm,1(T ) ? Gm / f1*(S) ) ,1 ? Gm,1 ? RT ln( f1

(2-1.29)

由于 [?(G T ) ?T ] p ? ? H T 2 ,因此

? H ? ?(?fus Gm,1 / T ) ? ? ? fus 2 m,1 ? ? ?T T ? ?p
( L) (S ) 如设 ?fus Cpm,1 ? Cpm ,1 ? C pm,1 不随温度而变,可得

(2-1.30)

?fus Hm,1 (T ) ? ?fus Hm,1 (Tf,1) ? ?fus Cpm,1 (T ? Tf,1)
代入式(2-1.30),并在 Tf,1 至 T 间积分得

(2-1.31)

?fus Gm,1 (T ) ? ?fus Hm,1 (Tf,1) ? (1 ? T / Tf,1) ? ?fus Cpm,1T ln(Tf ,1 / T ) ? ?fus Cpm,1 (T ? Tf,1)
以式(2-1.29)代入式(2-1.28),

(2-1.32)

ln x1 ? ??fus Gm,1 (T ) RT ? ln ? 1, I ?? ?fus H m,1 (Tf,1) ? 1 1 ? ?fus C pm ,1 ? ? Tf,1 ?? ? ? ?? ln( T / T ) ? ? ? f,1 ? T ? 1? ?? ? ln ? 1, I ?T T ? R R ? ?? f,1 ? ? ?
(2-1.33)

(L) (S ) *( S ) 以式(2-1.32)代入即可计算溶解度 x1。Tf,1、?fusHm,1(Tf,1)、 C pm 和 KH1 作为 ,1 、 C pm ,1 可代替 f1

表征系统的性质,前者比后者更容易得到。 如果溶液为理想溶液,则式(2-1.33)变为,
ln x1 ? ? ? fus H m,1 (Tf,1) ? 1 1 ? ? fus C pm ,1 ? ? Tf,1 ?? ? ? ?? ?ln(T / Tf,1) ? ? ? T ? 1? ?? ? ? R R ? ?? ? ? T Tf,1 ?

(2-1.34)

称为理想溶解度。由式(2-1.34)可见,理想溶解度只与溶质 2 的纯物质性质有关,而与溶剂的 性质无关。一般情况下,?fusCpm,1 很小,对溶解度计算的影响也不大,常常可以忽略。 如果在所关注的温度区间,溶剂 2 也可以以固体的形式析出,则式(2-1.33)同样适用于溶 剂 1,即
ln x2 ? ? ? ? C ? fus H m, 2 (Tf,2 ) ? 1 ? ? 1 ? ? fus pm , 2 ? ? R R ? T Tf,2 ? ? ? Tf,2 ?? ?ln(T / Tf,2 ) ? ? ? T ? 1? ?? ? ln ? 2, I ? ?? ?

(2-1.35)

当组分 1 和组分 2 同时析出时,称为最低共熔点,其温度和组成称为低共熔温度和低共 熔组成,它们可以通过联立求解式(2-1.33)和(2-1.35)得到。

2.2 气液相平衡计算——由 T、p、x 推算 y
气液平衡数据最常规的获得方法是采用各种循环式平衡釜抽取气、液相样品进行相组成 分析,并同时测定平衡时的温度 T 和压力 p,这是一种直接测定法。一些学者曾对不同气液 平衡釜的操作原理、结构、优缺点等作过详细的论述[2-6~2-9]。按照相律,在气液平衡的 T、p、 x 和 y 四个变量中,只有两个是独立变量,当给定了四个变量中的任意两个,并输入合适的能 表征该系统特征的分子热力学模型(状态方程或/和过量函数模型),即可根据相平衡原理计算 得到另两个变量,这就是前面所讨论的气液相平衡的有模型计算法。如果不是输入分子热力 学模型而是通过实验测定第三个变量,则原则上第四个变量可由它们通过相平衡计算确定。 由于计算过程中没有引入模型,称为气液平衡的无模型计算法。由于实际测定中气相取样和 组成测定是最困难的,一般是恒定温度 T 下由静态法测定一系列液相组成 x 对应的系统总压 p(溶液的饱和蒸汽压),或恒定压力 p 下由沸点仪测定一系列液相组成对应的系统温度 T(溶液 的沸点),然后由热力学原理计算与液相组成对应的气相组成 y,简单地说就是 T、p、x 推算 y。所谓无模型法并不是绝对的,有时仍需要采用合适的模型计算一部分非关键性的决定系统 特征的性质,例如气相逸度因子等。总之,至少是大大减少了模型的使用。国内外学者已建 立了多种 T、p、x 推算 y 的方法,它们都是在 Gibbs-Duhem(G-D)方程的基础上建立起来的。 根据应用 G-D 方程方式上的不同,可以归结为两大类:其一是直接法,它是将式(1-6.11)表示 的逸度的 G-D 方程同时应用于气液两相而得到联系 T、p、x 和 y 的共存方程,解此共存方程

即可实现由 T、p、x 推算 y 的目的;另一种是间接法,它首先计算过量吉氏函数 Q,根据 Q 与活度因子的关系(隐含了 G-D 方程)计算液相活度因子, 从而实现间接计算气相组成的目的。

2.2.1 Q 函数法(间接法)
2-2.1.1 Q 函数法原理 汽液平衡时,按判据式(1-6.13), fk(V) ? fk(L) (k=1, …, K),如气相采用逸度因子、液相采用 第 I 种活度因子分别计算气液相的非理想性,得
* * *L * pyk?k ? pk ?k xk? k ,I exp[ Vm , k ( p ? pk ) / RT ] , k ? 1, 2, ?, K

(2-2.1)

整理上式可得系统总压 p,
* * *L * p ? ? pyk ? ? pk ?k xk ? k , I exp[ Vm , k ( p ? pk ) / RT ] / ? k k ?1 k ?1 K K

(2-2.2)

式中 ? k , I 用式(1-7.46、7.47)代入,得
? K ?1 ? ?Q ? ? ?Q ? ? ? ? ? ? p?? exp?Q ? ? ? x ? j ? ? ?x ? ?k ? x k ?1 j ? 1 ? k ? x[ k , K ] ? ? j ? x[ j , K ] ? E ? ? V E ?? ?p ? ?? K ?1 K ?1 ? ?T ? ? ?p ? ? ?T ? Hm m ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? xj ? ? ? xj ? 2 ? ? ? ? ? ? ? ? RT ?? ?xk ? x[ k , K ] j ?1 ? ?x j ? RT ?? ?xk ? x[ k , K ] j ?1 ? ?x j ? ? ? x[ j , K ] ? x[ j , K ] ? ? ? ? ?
K * * *L * pk ? k xk exp[ Vm , k ( p ? pk ) / RT ]

(2-2.3)

*L * * 注意当 i=K, 式中对 xK 的偏导数全为零。 式(2-2.3)的意义在于: 如果暂时不考虑 pk 、?k 、 Vm ,k 、
E E 和 Vm ,则式中除了 Q 以外,其它的变量就是已输入的 T、p、x。而 Q 函数正是 T、 ?k 、 H m

p、x 的函数,式(2-2.3)实质上是一个 Q 函数的偏微分方程,只要有足够数量的一系列 T、p、 x 的实验数据, 原则上可以解得 Q=Q(T, p, x)。 有了 Q, 可用式(1-7.46、 47)计算?k,I, 代入式(2-2.1)
*L * * 即可求得 y。至于那些暂时放在一边的变量:其中 pk 、 ?k 和 Vm , k 是纯组分性质,与混合物无

关。?k 决定于气相组成 y,可利用上次迭代的 y 值计算,但还需要使用合适的状态方程,从 这个意义上说,T、p、x 推算 y 并不是完全的无模型,但当压力不太高时,气相非理想性远 没有液相的那样强烈,在压力较低时,采用截止到第二维里系数的维里方程足以估算这种非
E E 理想性,甚至可以令?k=1,也不致带入严重误差。至于 H m 和 Vm ,后者很小,常可忽略,前

者对于恒温数据不起作用,对于恒压数据,实践证明略去后影响不大。总之,这一方法基本 上不使用模型,或者严格地说,不使用液相活度因子模型,而它是气液平衡计算中最关键的 模型。 式(2-2.3)原则可以求解,但实践上却有很大困难,因为导数出现在 exp 中,是一个超越

型的偏微分方程,没有解析解,只能通过数值方法求解。国内外学者已发展了多种方法,根 据所采用数值方法的不同,可以分为几种类型:最早是 Barker[2-10,11]提出的方法,其核心是选 择一个过量函数模型代入式(2-2.3),利用一系列 T、p、x 的实验数据,拟合得到模型参数和 Q 函数。第二种是 Mixon 等[2-12]发展的有限差分法,它以差分来逼近式(2-2.3)中的导数,然后利 用 Newton 法迭代求得离散格点上的 Q 值。这种方法不依赖于任何过量函数模型,是严格的 无模型法[2-13,14]。它对二元系的计算非常成功,得到广泛的应用。但用于三元系时,收敛速度 极慢,且求解过程不稳定。第三种是样条函数法,包括适用于二元系的三次样条函数法[2-15,16] 和适用于任意组分数的曲面样条函数法[2-16, 17,18]。特别是曲面样条函数法,它不仅能方便地用 于二元系和三元系,也能成功地应用于多元系,更重要的是不同组分数的计算方法可以统一 在一个框架下。大量实例计算表明,没有收敛的困难,不受多元系 Q 函数曲面类型的限制。 下面介绍计算简单的 Barker 法和适用于多元系的曲面样条函数法,至于其它方法,感兴 趣的读者可以参考相关的文献和著作。 2-2.1.2 Barker 法 Barker 法[2-10,11]的核心是选择一个合适的过量函数模型,其中包括若干待定模型参数, 代入式(2-2.3)后,利用一系列 T、p、x 的实验数据,拟合得到模型参数,这就得到 Q 函数。 原则上,第 4 章中介绍的各种过量函数或活度因子模型均可以作为 Barker 法的候选模型。但 这些模型都是针对特定的对象而建立起来的,都有一定的适用范围,这就使得 Barker 法的准 确度受到所选模型可靠性和适用性的限制,对多元系问题会更突出。解决的办法是尽可能选 用灵活性大的经验模型, 例如对于二元系, 可以采用如下的 Redlich-Kister 型的经验 Q 函数模 型,

Q ? x1 x2 ? j ? 0 Aj ( x1 ? x2 ) j
N

?

?
N

(2-2.4)

相对应的活度因子为

2 ln ? 1, I ? x2 ? j ?0 Aj ( x1 ? x2 ) j ? 2 x1x2 ? j ?0 jAj ( x1 ? x2 ) j ?1 N

ln ? 2, I

? ? x ??
2 1

N j ?0

? ? A ( x ? x ) ?? 2 x x ??
j j 2 1 1 2

N j ?0

jAj ( x2 ? x1 ) j ?1

? ?

(2-2.5) (2-2.6)

根据 T、p、x 实验数据的多少和计算精度的要求,可以选择不同的 N 值。 式(2-2.4)是一个关于 Aj 的线性方程,如果已知不同组成下的 Q 函数值,可以非常方便地 采用最小二乘法关联得到 N+1 个 Aj。具体计算时,我们可以采用如下的迭代过程: (1) 假设气相为理想气体、液相为理想溶液,计算气相组成的初值 y0;
* (2) 由状态方程计算气相逸度因子 ?k 和?k,由式(2-2.7)计算各组分的液相活度系数;

* * *L * ? k ,I ? pyk?k /?pk ?k xk exp[ Vm ,k ( p ? pk ) / RT ]?

(2-2.7)

(3) 由式(2-2.8)计算各实验点的过量吉氏函数 Q;
Q ? ? xk ln ? k ,I
k ?1 K

(2-2.8)

(4) 由计算得到的过量吉氏函数 Q 关联式(2-2.4)中的未知参数 Aj(j=0, …, N); (5) 由式(2-2.5、2.6)计算各组分的活度因子?k,I; (6) 由式(2-2.9)计算各组分新的气相组成 y1;
* * *L * yk ? pk ?k xk? k ,I exp[ Vm ,k ( p ? pk ) / RT ] / p?k , k ? 1, ?, K

(2-2.9)

(7) 比较 y1 和 y0,如果两者不相等,则令 y0=y1,转步骤(2),进行新一轮循环迭代,直 至达到规定的计算进度。 2-2.1.3 曲面样条函数法 Q 函数法 T、p、x 推算 y 的核心是构造 Q 函数曲面,并由 T、p、x 实验数据通过式(2-2.3) 确定该 Q 函数曲面。先看三元系,这时 Q 函数是 x1 和 x2 的函数 Q(x1, x2),Q-x1- x2 形成三维 空间,对于液相完全互溶的系统,Q 函数是一个定义在[x1?0; x2 ?0; x1+x2?1]三角形区域内的 连续曲面, 且在三角形的三个顶点处 Q=0。 设在 Q 函数定义域内的 N 个结点上的 Q 函数值已 知, Qi(x1(i), x2(i))(i=1, …, N)。 为构造一个曲面样条函数 Q(x1, x2), 要求在每个结点上 Q(x1(i), x2(i))= Qi,可将 Q 看成是一块无限大平板纯弯曲时的变形[2-19],Q 即为板的挠度,它和作用于该板 上的负载 q 之间存在如下微分方程,

D?4Q ? q

(2-2.10)

式中 D 为板的刚度。给定 N 个结点上的挠度,并考虑到远离所加负载处曲面变得平坦所引出 的条件,以及将分布负载代替点负载,然后应用以弹性力作用于平板代替强制平板通过 N 个 结点的假定,曲面样条函数即可表达为:
Q( x1 , x2 ) ? ?[di Ri2 ln(Ri2 ? ? )] ? d N ?1 x1 ? d N ? 2 x2 ? d N ? 3 ? c
i ?1 N

(2-2.11)

式中

Ri2 ? ( x1 ? x1(i ) )2 ? ( x2 ? x2(i ) )2
相应地可写出两个一阶偏导数为
N ? ? Ri2 ?? ? ?Q ? ? 2 ? ? ? 2 d ( x ? x ) ? ln(Ri2 ? ? ) ? ? ? i 1 1( i ) ? ? ?x ? ?? ? d N ?1 i ?1 ? ? 1 ? x2 ? Ri ? ? ??

(2-2.12)

(2-2.13)

N ? ? Ri2 ?? ? ?Q ? ? 2 ? ? ? 2 d ( x ? x ) ? ln(Ri2 ? ? ) ? ? i 2 ? 2(i ) ? ? ?x ? ?? ? d N ? 2 i ?1 ? ? 2 ? x1 ? Ri ? ? ??

(2-2.14)

式中?是个小量, 可根据曲面曲率大小取值, 曲率变化平缓时, ?=1~10-1, 变化大时?=10-4~10-6, 一般可取?=1~10-2,奇性曲面?=10-5~10-6。di 是未定系数,应满足以下三个条件,
G1 ? ? x1(i ) di ? 0
i ?1 N N

(2-2.15)

G2 ? ? x2(i ) di ? 0
i ?1 N

(2-2.16)

G3 ? ? di ? 0
i ?1

(2-2.17)

式(2-2.11)的含义在于:某点的挠度是 N 个结点上所施加负载所引起的总后果。在 N 个结点上 则有:
2 2 Q j ? ? di Rij ln(Rij ? ? ) ? d N ?1 x1( j ) ? d N ? 2 x2( j ) ? d N ? 3 ? c j i ?1 N

(2-2.18)

式中
2 Rij ? ( x1( j ) ? x1(i ) )2 ? ( x2( j ) ? x2(i ) )2

(2-2.19) (2-2.20)

c j ? 16?D / k j

kj 是结点 j 的弹性常数,cj 是对曲面光顺要求所加的权重,在某点上该值取为零,则所构作的 样条曲面通过该点的 Q 坐标处, 即通过该已知结点(x1(i), x2(i), Qi)。 求解由式(2-2.15、 2.16、 2.17) 和(2-2.18)组成的 N+3 个线性方程组,可以解得 N+3 个 di 。从而得到曲面样条函数 Q(x1, x2)。 当各个独立点上的 Q 函数是未知的,但它满足某个偏微分方程,即
D(Q) ? 0

(2-2.21)

将式(2-2.11)代入上式,对每个独立点有,

Di ( x1(1) , x2(1) , ..., x1( N ) , x2( N ) , d1, ..., dN ?1, d N ? 2 , dN ?3 ) ? 0, i ? 1, ..., N

(2-2.22)

联系式(2-2.18)中 di 必须满足的三个条件, 即式(2-2.15~2.17), 可得到联系 N+3 个 di 的方程组。 如果 D(Q)是个线性偏微分方程,得到的是线性方程组,反之则得到非线性方程组。求解由式 (2-2.22)和式(2-2.15~2.17)构成的方程组,可以得到 N+3 个 di,通过式(2-2.11)可得任意坐标下 的 Q 函数。这就是曲面样条函数法求解偏微分方程的基本原理。 注意到式(2-2.12)中的 Ri2 实际上是 Q 函数定义域中的任意一点与结点 i 的距离的平方, 我

们很容易将上述曲面样条函数法求解偏微分方程的方法推广至多维空间。设混合物为 K 元多 组分系统,这时,Q 是定义在[xk?0, (k=1, …, K-1) ;

?

K ?1 k ?1 k

x ? 1 ]区域内的连续曲面。定义

x=(x1, …, xK-1, 1)T, 则 Q=Q(x)。 若该区域中 N 个点上的 Q 函数值 Qj 已知, 则按照上面的讨论, 通过这 N 个点的曲面样条函数可以表示为
Q(x) ? ?[di Ri2 ln(Ri2 ? ? )] ? gT x
i ?1 N

(2-2.23)

这里 Ri2 的意义仍然是 Q 函数定义域中的任意一点与结点 i 的距离的平方,只不过现在是在多 维(K 维)空间中,因此有
Ri2 ? ? ( xk ? xk (i ) )2 ? (x ? xi )T (x ? xi )
k ?1 K ?1

(2-2.24) (2-2.25)

g ? (d N ?1, ..., d N ? K )T
这样式(2-2.23)的含义与式(2-2.11)的含义是一致的。N 个 di 之间满足如下关系

?

N

i ?1

d i xi ? 0

(2-2.26)

0 为有 K 个元素的零向量。曲面样条函数的一阶导数为,
N ? ? Ri2 ?? ? dQ ? ? ln(Ri2 ? ? ) ? ? ? ? 2? ?di (x ? xi )? 2 ? ?? ? g ? dx ? i ?1 ? ? Ri ? ? ??

(2-2.27)

因为 x 的第 K 个分量为常数 1,所以

? dQ ? ? ? ?0 ? dx ? K
在 N 个结点上则有,
2 2 Q j (x j ) ? ?[di Rij ln(Rij ? ? )] ? gT x j i ?1 N

(2-2.28)

(2-2.29)

式中
2 Rij ? ? ( xk ( j ) ? xk (i ) )2 ? (x j ? xi )T (x j ? xi ) k ?1 K ?1

(2-2.30)

相应地,有
2 N ? ? Rij ?? ? dQ ? 2 ?? ? g ? ln( R ? ? ) ? ? ? 2? ?di (x j ? xi )? ij ? R2 ? ? ? ? dx ? j i ?1 ? ij ? ?? ? ?

(2-2.31)

联立求解由式(2-2.16)和式(2-2.29)组成的方程组,可以得到 N+K 个 di。

若 Q 满足偏微分方程式(2-2.21),将曲面样条函数式(2-2.23)代入,在 N 个结点上则有,

Dj (x j , x1, ..., x( N ) , d1, ..., d N ? K ) ? 0,

j ? 1, ..., N

(2-2.32)

联立求解式(2-2.16)和式(2-2.32)组成的方程组,同样可以得到 N+K 个 di,从而获得偏微分方 程式(2-2.21)的近似解。值得注意的是,求解多维空间的偏微分方程与求解三维空间的偏微分 方程的方法是相同的,只要结点数不变,需要求解的曲面样条函数中的参数的数目变化很小 (每增加一个维度,需要增加一个 d 参数),这就为多元系由 T、p、x 推算 y 提供了极大方便。 对于由组分 1, 2, …, K 组成的 K 元多组分系统,假定已由实验测得 n 组(T, p, x)数据,连 同 K 个纯物质的 T 和 p 数据,总共有 N=n+K 个实验点。根据过量函数 Q 与活度因子的关系, 我们可以得到系统总压的计算式为,
pcalc ? ?
k ?1 K * * pk ?k xk

?k

T ? ? ? ? K ? dQ ? exp?Q ? ? ? (x ? ek )? ? ? pk ? ? ? dx ? ? ? k ?1

(2-2.33)

E E 其中,ek=(?kj),对 k?j, ?kj =0;对 k=j, ?kj =1。上式中已经忽略了 H m 、 Vm 和 Poynting 因子的

贡献,因为它们对间接法推算结果的影响很小,可能的原因是由此引起的误差被 Q 函数吸收 了。对前 n 组(T, p, x)数据,以系统总压的计算值与实验值之差为偏差函数,

Fj ? ( pcalc,j ? pexp, j ) / pexp, j ,j=1, …, n

(2-2.34)

则 F 是 Q 函数的一个非线性偏微分方程。将曲面样条函数式(2-2.23)代入式(2-2.34),则 F 就 转化为 N+K=n+2K 个 di 的非线性方程组。对于 K 个纯物质,应满足如下条件,

Q(x ? ek ) ? 0 ,k=1, …, K
所以,后 K 点数据的偏差函数可简单取为

(2-2.35)

Fn ? k ? Q(ek ) , k=1, …, K
外加 n 个 di 之间必须满足的 K 个条件式(2-2.26),N+1~N+K 的偏差函数取为,
FN ? k ? ?i ?1 di xi (k ) , k=1, …, K
N

(2-2.36)

(2-2.37)

xi(k)为 xi 的第 k 个元素。 这样,由式(2-2.34 、 2.36 、 2.37)组成了有 n+2K 个关系的非线性方程组,从中可解得 N+K=n+2K 个 di。再由

yk ? pk / pcalc , k=1, …, K
可得各组(T, p, x)数据对应的气相组成 y=(y1, …, yK)T。令,

(2-2.38)

F ? ( F1, ..., FN ? K )T

(2-2.39)

X ? (d1, ..., d N ? K )T
则上述 N+K 个偏差函数组成的非线性方程组可用向量函数的形式表示为
F(X) ? 0

(2-2.40)

(2-2.41)

这 里 的 0 为 有 N+K 个 元 素 的 零 向 量 。 对 式 (2-2.41) 的 求 解 可 采 用 Broyden 改 进 的 Newton-Raphson 法(见附录 2-1),这是一个非常有效的求解非线性方程组的方法[2-20]。

2.2.2 直接法
2.2.2.1 直接法原理 直接法是从逸度的 Gibbs-Duhem 方程出发建立起来的 T、p、x 推算 y 的方法[2-14,16]。对于 一个含有 K 个组分的系统,按式(1-6.11),其液相逸度的 Gibbs-Duhem 方程为,
f xk d ln ? p k ?1
K (L) k o

?

?x H
k ?1 k

K

o m, k

(L) ? Hm

RT 2

(L) Vm dT ? dp RT

(2-2.42)

当处理一系列 T、p、x 实验数据时, T、p 和 f k(L) 均可形式上表达为液相组成 x1、x2、…、xK-1 的函数,上式变为
? ? ln( f p )? ? xk ? ? dx j ? ? ? ? ?x j k ?1 j ?1 ? ? x[ j , K ]
K K ?1 (L) k θ

?x H
k ?1 k

K

o m, k

(L) ? Hm 2

RT

?? ? ?x

K ?1

? ?T ? ? dx j ? j ?1 ? j ? x[ j , K ]

(2-2.43)

?

V RT

( L ) K ?1 m

?? ? ?x

? ?p ? ? dx j ? 0 ? j ?1 ? j ? x[ j , K ]

变换求和次序,
K ? o (L) xk H m ? ,i ? H m (L) ? ? K ?1 ? K ? ? ln( f k( L ) / p o ) ? ? ?T ? Vm ?p ? i ?1 ? ? ? ? ? ? ? dx j ? 0 ? x ? ? ? ? ? k ? ?x ? ? ?x ? ? ? RT 2 j ?1 ? k ?1 j ? ? x[ j , K ] ? j ? x[ j , K ] RT ? ?x j ? x[ j , K ] ? ? ? ?

(2-2.44)

由于 K-1 个 dxj 可以独立变化,其系数必需等于零,

? ? ? ln( f k( L ) p o ) ? k ?1 ? ? x ? ? k? ? ? x k ?1 j ? ? x[ j , K ]
K

K

o (L) xk H m ,k ? H m

RT 2

? ?T ? V ( L ) ? ?p ? ? ? ? ? m ? ?0 ? ?x ? ? ?x ? RT ? j ? x[ j , K ] ? j ? x[ j , K ]

(2-2.45)

气液平衡时, fk(L) ? f k( V) ? pyk?k ,代入上式,得

K ? ? ln( p / p o ) ? ? ? ln yk ? ? x ? xk ? ? ? k? ? ? ?x j k ?1 ? ? x[ j , K ] k ?1 ? ?x j K

K ? ? ? ln ?k ? ? ? xk ? ? ? ? x[ j , K ] k ?1 ? ?x j (L) m

? ? ? ? x[ j , K ]

?
其中

?x H
i ?1 k

K

o m, k

(L) ? Hm 2

(2-2.46)

RT

? ?T ? V ? ? ? ? ?x ? ? j ? x[ j , K ] RT

? ?p ? ? ? ?0 ? ?x ? j ? ? x[ j , K ]

? ? ln yk xk ? ? ? ?x k ?1 j ?
K K

K ?1 ? ? xk xK ?? ?yk ? ? ? ?? ?y ?y ? ?? ?x ? k ? 1 k K ? ? ? x[ j , K ] ? j

? ? ? ? x[ j , K ]

(2-2.47)

? ? ln( p / p o ) ? 1? ?p ? ? ? ? ? x ? ? k? ? ? ? ? x p ? x k ?1 j ? ? x[ j , K ] ? j ? x[ j , K ]
(L) *( L ) E Hm ? ? xk H m ,k ? H m k ?1 K

(2-2.48)

(2-2.49) (2-2.50) (2-2.51)

o *( L) Hm , k ? H m, k ? Lk
K

( L) *( L ) E Vm ? ? xkVm , k ? Vm k ?1

*( L ) E E Lk 和 Vm , k 分别是纯液体 k 的蒸发焓和摩尔体积; H m 和 Vm 分别是液体混合物的过量摩尔焓和

过量摩尔体积,后者很小,一般情况下可以忽略不计。代入式(2-2.46),并用 Fj 表示,
K ?1 ? xk xK ?? ?yk ? Fj ? ? ? ?y ?y ? ?? ?x k ?1 ? k K ?? j K ? ? ? ln ? k ? ? ? ? ? xk ? ? ? ? ? x j ? x[ j , K ] k ?1 ? ? x[ j , K ] K K E ? 1 ? x V *( L ) Hm ? ?k xk Lk ? ?T ? ? ? k k m, k ? ? ? ? 2 ? ? ? RT RT ? ?x j ? x[ j , K ] ? p

?? ?p ? ?? ? ?0 ?? ?x j ? ? x[ j , K ] ??

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

(2-2.52)

因为?k 是 T、p 和 y 的函数,可以写出

? ? ln ?k ? ? ?x j ?

K ?1 ? ? ? ln ?k ? ? ?? ? ? ? x[ j ,K ] i?1 ? ?yi

? ? ?yi ? ? ? ? ln ?k ? ? ln ?k ? ? ? ?T ? ? ? ? ? ?? ? ? ? ? ? ? ? ?x ? ? y[i ,K ] ? j ? x[ j ,K ] ? ?T ? y ? ?x j ? x[ j ,K ] ? ?p

? ? ? ? ?p ? ? (2-2.53) ? ? ?x ? ? y ? j ? x[ j , K ]

代入式(2-2.52),得

K ?1 K K ?1 ? ?yi ? ? xk xK ?? ?yk ? ? ? ln ? k ? ? ? ? ? ? ? ? Fj ? ? ? ? ? x ? ? k ?y ? ? ? ? ? ?x ? yK ?? ? x ? y k ?1 ? k k ? 1 i ? 1 j i j ? ? ? ? x[ j , K ] ? x[ j , K ] y[ i , K ] ? K E K ? H m ? ? xk Lk ? ?T ? ? ? ln ? k ? ? k ? ? ? ?? ? x ? ? ? k ? ? RT 2 ? x ? T ? ? ? ? k ?1 j ? x[ j , K ] y ? ? ? *( L ) K ? 1 ? K xkVm ? ? ln ? k ,k ?? ? k ? ? xk ? ? ?p p RT ? k ?1 ? ? ?0

(2-2.54)

? ? ?? ? ?p ? ? ? ? ? ?x ? ?y ? ? ? j ? x[ j , K ]

上式第二项交换 i 和 k,并整理后得
K ?1 ? K ?? ?y ? ? ln ?i ? x x ?? k ? F j ? ? ? k ? K ? ? xi ? ? ? ? yK i ?1 ? ?yk ? y[ k , K ] ?? ?x k ?1 yk ? ?? j

? ? ? ? x[ j , K ]

K E K ? Hm ? ?T ? ? ? k xk Lk ? ? ln ? k ? ? ? ? ? ?? ? x ? ? ? k 2 ? ? RT ? x ? T ? ?y ? ? k ?1 ? ? ? j ? x[ j , K ] *( L ) K ? 1 ? K xkVm ? ? ? ln ? k ? ? ? ,k ? ?p ? ? ?? ? k ? ? ? xk ? ? ?p ? ? ?x ? p RT ? k ?1 ? ?y ? ? ? ? j ? x[ j , K ] ?0

(2-2.55)

此式共 K-1 个,每个式子有三项,分别含有 yk、T、p 随液相组成变化的偏导数。由于一系列
E * (L) T、p、x 已由实验提供,如进一步输入 Lk、Vm , k 和 H m ,?k 则由合适的状态方程计算,所需 y

使用上次迭代值,用适当方法求解式(2-2.55),原则上可以直接求得 y。 对于对于二元系,式(2-2.55)变为:
? x 1? x ? ln ?1 ? ln ? 2 ? dy F ?? ? y ? 1 ? y ? x ?y ? (1 ? x) ?y ? ? dx ? ?
E ? Hm ? xL1 ? (1 ? x) L2 ? ln ?1 ? ln ? 2 ? d T ?? ?x ? (1 ? x) 2 RT ?T ?T ? ? ? dx *( L ) *( L ) ? 1 xVm p ? ln ?1 ? ln ? 2 ? ,1 ? (1 ? x)Vm , 2 ?d ?? ? ? x ? ( 1 ? x ) ?p RT ?p ?p ? ? ? dx ?0

(2-2.56)

2.2.2.2 数值积分法 式(2-2.56)可进一步改写成,

E ? xL1 ? (1 ? x) L2 dy ?? H m ? ln ?1 ? ln ? 2 ? d T ? ?? ?x ? (1 ? x) 2 RT dx ?? ?T ?T ? ? dx *( L ) *( L ) ? 1 xVm ? ln ?1 ? ln ?2 ,1 ? (1 ? x)Vm , 2 ?? ? ?x ? (1 ? x) ?p RT ?p ?p ?

? dp? ? ? ? dx ? ? ? ?

(2-2.57)

? x 1? x ? ln ?1 ? ln ? 2 ? /? ? y ? 1 ? y ? x ?y ? (1 ? x) ?y ? ? ? ?
如果实验测定了恒温下整个浓度范围内的 p-x 数据或恒压下的 T-x 数据, 上式就是 y 与 x 的一个非线性常微分方程,可以用任何满足精度要求的数值积分方法求解,最常用的是 Runge-Kutta 法,参见附录 2-2[2-14]。 积分求解过程中,初始值可以取(x=0, y=0),然后从 x=0 积分到 x=1。也可以取初始值为 (x=1, y=1),然后从 x=1 反向积分至 x=0。 计算时的主要困难在于具有恒沸点的系统。从式(2-2.57)可以看到,如果忽略气相非理想 性,则恒沸点时有(dp/dx)=0 或(dT/dx)=0,(dy/dx)为不定值,通常的数值积分方法难以越过此 点。克服这一困难的办法是分别从(x=0, y=0)和(x=1, y=1)两端向恒沸点积分。数值积分法的另 一个困难是,随着积分的进行误差不断累积,可能出现 x=1 时 y?1 的情况。 2.2.2.3 有限差分法 将整个液相组成范围从 x=0 到 x=1 划分为等距离的 M+1 个格点,M 可取 100~1000,间 距?x=M-1,xk=i?x,i 为格点序码。将 i 点处的函数 Fi 写成差分形式,
? i?x 1 ? i?x ? ? ln ?1 ? ? ? ln ? 2 ? ? yi ?1 ? yi ?1 ? ? ? ?? Fi ? ? ? y ? 1 ? y ? i?x? ?y ? ? (1 ? i?x)? ?y ? ? 2?x ? ?i ? ?i ? i ? i ? H E ? i?x L1,i ? (1 ? i?x) L2,i ? ? ln ?1 ? ? ? ln ? 2 ? ?? d T ? ? ? m,i ? i?x? ? ? (1 ? i?x)? ? ?? ? 2 RTi ? ?T ?i ? ?T ?i ?? d x ?i ? *( L ) *( L ) ? 1 i?xVm ? ? ln ?1 ? ? ? ln ? 2 ? ? dp ? ,1 ? (1 ? i?x )Vm , 2 ?? ?? ? ? ? ? ? ? i ? x ? ( 1 ? i ? x ) ? ? ? ? ? ? ?p ? RTi ? ?p ?i ? ?p ?i ?? d x ?i ? i ?0

(2-2.58)

此方程除两端点(x=0, 1)外, 共可列出 M-1 个, 未知数是 M-1 个格点上的 yi, 原则上可以求解。 令

F ? ( F1, F2 ,?, FM ?1 )T y ? ( y1, y2 ,?, yM ?1 )T
可将 M-1 个式(2-2.58)写成函数向量式,
F ( y) ? 0

(2-2.59) (2-2.60)

(2-2.61)

利用 Newton-Raphson 迭代法求解,

?y( k ) ? ?(J ( k ) )?1 F ( y( k ) ) y ( k ?1) ? y( k ) ? ?y( k )

(2-2.62) (2-2.63)

式中 k 为迭代序号,J 为 Jacobi 矩阵,即 F 对 y 的偏导数矩阵,是(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 用上次迭代值。
E 直接法可以成功地用于二元系气相组成的推算,当应用于恒压系统时, H m 不可忽略,这

一点和间接法有很大区别。具体例子和讨论参见文献[2-21、22、23]或[2-14]。对于三元系, 式(2-2.58)可以写出两个式子,仍可采用有限差分法,但计算效率较低,收敛困难,具体算例 见文献[2-21]。 对式 2-2.56)的求解,也可以借助三次样条函数的方法[2-24]。但总体上来说,直接法的计算 效率没有间接法高,尤其是直接法推广至三组分以上的多元系非常困难,按式(2-2.55),这时 需要求解偏微分方程组,而间接法则只需要求解一个偏微分方程式(2-2.3)。

2.2.3 T、p、x 推算 y 计算实例
由 T、p、x 推算 y 的计算过程中,需要有计算气相逸度因子?k 的状态方程,可以选用第 3 章中介绍的合适的状态方程,以下介绍的计算实例均选用了改进的 Peng-Robinson 方程
[2-25,26]

。在曲面样条函数法计算中,di 的初值全部取为 di=0.001,对一般系统而言,Q 函数是

不存在奇性的,曲率变化比较平缓,所以取 ? =0.001,cj 均取为零以使得到的曲面样条函数通 过各结点的 Q 坐标值处。 表 2-2.1 是不同方法对 55℃ 下氯仿-乙醇二元系的计算结果与实验值的比较。从表中结 果可见,有限差分直接法、有限差分间接法、三次样条直接法和曲面样条间接法推算的结果 具有相同的准确度,相互之间能很好的吻合。计算结果表明,是否考虑气相非理想性和是否 忽略液体体积,对计算结果的影响甚微。以三次样条直接法为例,如果同时忽略气相非理想 性和忽略液体体积,气相组成的推算误差从 0.0016 增加到 0.0038。 表 2-2.2 是对 101.325kPa 下苯-异丙醇二元系的计算结果与实验值的比较。从表中结果 可见,有限差分直接法、三次样条直接法和曲面样条间接法推算的结果具有相同的准确度, 相互之间能很好的吻合,但间接法具有不需要过量焓数据的优点。以三次样条直接法为例, 如果忽略过量焓的影响,气相组成的推算误差将从 0.0076 增大到 0.0106;如果忽略气相非理

想性,则误差增大到 0.0105;如果同时忽略气相非理想性和液相过量焓,误差增大到 0.0120。 如果采用有限差分直接法,忽略过量焓的结果将导致气相组成的推算误差从 0.0058 升高到 0.0181。可见过量焓对直接法推算结果的影响还是不小的。

表 2-2.1 氯仿-乙醇二元系推算结果比较(55℃ ) x 0.00 0.01 0.02 0.05 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 0.95 0.98 0.99 1.00 p/kPa 37.312 38.377 39.443 42.748 48.572 60.262 69.480 76.076 80.691 83.829 85.813 86.957 86.772 85.622 83.985 83.196 82.372 ?y y(exp) 0.0000 0.0372 0.0726 0.1694 0.3054 0.4988 0.6108 0.6774 0.7252 0.7598 0.7914 0.8264 0.8774 0.9226 0.9634 0.9800 1.0000 y(a) 0.0000 0.0371 0.0722 0.1697 0.3079 0.5027 0.6127 0.6803 0.7272 0.7632 0.7938 0.8276 0.8789 0.9208 0.9620 0.9805 1.0000 0.0020 y(b) 0.0000 0.0366 0.0717 0.1696 0.3077 0.5024 0.6123 0.6799 0.7268 0.7628 0.7935 0.8274 0.8790 0.9209 0.9630 0.9806 1.0000 0.0018 y(c) 0.0000 0.0361 0.0718 0.1697 0.3088 0.5008 0.6133 0.6784 0.7278 0.7604 0.7964 0.8263 0.8786 0.9207 0.9620 0.9805 1.0000 0.0016 y(d) 0.0000 0.0372 0.0725 0.1702 0.3089 0.5032 0.6130 0.6804 0.7274 0.7631 0.7941 0.8269 0.8783 0.9205 0.9614 0.9801 1.0000 0.0016

注:exp-实验值;a-有限差分直接法;b-有限差分间接法;c-三次样条直接法;d-曲面样条间接法

表 2-2.2 苯-异丙醇二元系推算结果比较(101.325kPa) x 0.000 0.084 0.126 0.153 0.199 0.240 0.291 0.357 0.440 0.492 0.556 0.624 0.762 0.836 0.887 0.936 0.972 1.000 t/℃ 82.2 78.5 77.1 76.3 75.3 74.4 73.6 73.0 72.4 72.3 72.2 72.0 72.4 73.0 73.8 75.5 77.5 80.1 HE/J.mol-1 0.0 1054.37 1305.41 1456.03 1615.02 1721.72 1824.22 1899.54 1914.18 1889.08 1809.58 1686.15 1317.96 979.06 732.20 447.69 209.20 0.0 ?y y(exp) 0.000 0.208 0.276 0.316 0.371 0.410 0.451 0.493 0.535 0.558 0.583 0.612 0.673 0.717 0.760 0.825 0.901 1.000 y(a) 0.0000 0.2052 0.2754 0.3151 0.3718 0.4139 0.4572 0.4954 0.5358 0.5578 0.5884 0.6272 0.6793 0.7148 0.7573 0.8218 0.9046 1.0000 0.0058 y(b) 0.0000 0.2066 0.2810 0.3214 0.3738 0.4195 0.4630 0.4956 0.5432 0.5508 0.5835 0.6339 0.6959 0.7065 0.7581 0.8189 0.9037 1.0000 0.0076 y(c) 0.0000 0.2099 0.2798 0.3196 0.3779 0.4205 0.4626 0.5013 0.5383 0.5587 0.5931 0.6251 0.6688 0.7114 0.7528 0.8182 0.8986 1.0000 0.0053

注:exp-实验值;a-有限差分直接法;b-三次样条直接法;c-曲面样条间接法

表 2-2.3 列出了对 74.84℃ 和 124.98℃ 下二氯甲烷-硝基甲烷二元系的推算结果。从表中 可见,有限差分间接法、三次样条直接法、三次样条间接法和曲面样条间接法四种方法所得 气相组成相互之间能很好的吻合。对其它系统的计算结果也取得令人满意的结果[2-24]。

表 2-2.3 二氯甲烷-硝基甲烷二元系推算结果比较 x 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 p/kPa 62.789 81.883 99.470 115.758 130.942 145.216 158.757 171.667 184.027 195.916 207.414 218.600 229.555 240.364 251.117 261.901 272.805 283.915 295.297 t = 74.84℃ y(a) y(b) 0.3603 0.3606 0.5310 0.5323 0.6313 0.6326 0.6979 0.6989 0.7458 0.7468 0.7824 0.7834 0.8116 0.8124 0.8357 0.8363 0.8562 0.8565 0.8740 0.8742 0.8898 0.8900 0.9040 0.9045 0.9172 0.9177 0.9295 0.9301 0.9414 0.9418 0.9529 0.9531 0.9643 0.9643 0.9757 0.9758 0.9875 0.9874 y(c) 0.3581 0.5301 0.6313 0.6984 0.7465 0.7831 0.8122 0.8363 0.8566 0.8744 0.8900 0.9042 0.9173 0.9297 0.9414 0.9530 0.9643 0.9758 0.9875 y(d) 0.3616 0.5324 0.6323 0.6989 0.7465 0.7832 0.8122 0.8365 0.8568 0.8747 0.8903 0.9047 0.9177 0.9302 0.9418 0.9534 0.9762 0.9645 0.9875 p/kPa 263.895 320.260 372.645 421.478 467.450 511.294 553.636 594.605 634.176 672.326 709.088 744.707 779.470 813.667 847.563 881.380 915.334 949.778 985.729 y(a) 0.2580 0.4117 0.5143 0.5883 0.6449 0.6904 0.7285 0.7610 0.7892 0.8140 0.8362 0.8586 0.8753 0.8933 0.9106 0.9277 0.9448 0.9623 0.9809 t = 124.98℃ y(b) 0.2580 0.4116 0.5145 0.5890 0.6461 0.6917 0.7293 0.7613 0.7891 0.8137 0.8360 0.8565 0.8756 0.8936 0.9110 0.9279 0.9450 0.9625 0.9811 y(c) 0.2557 0.4103 0.5140 0.5887 0.6458 0.6914 0.7295 0.7619 0.7901 0.8147 0.8369 0.8570 0.8759 0.8937 0.9111 0.9280 0.9452 0.9625 0.9811 y(d) 0.2599 0.4140 0.5167 0.5906 0.6473 0.6928 0.7309 0.7634 0.7915 0.8164 0.8385 0.8587 0.8775 0.8954 0.9125 0.9295 0.9463 0.9635 0.9815

注:a-有限差分间接法;b-三次样条直接法;c-三次样条间接法;d-曲面样条间接法

表 2-2.4 是苯(1)-正己烷(2)-环己烷(3)三元系 70℃ 恒温数据的推算结果,计算时已略去 逸度校正。从表中可见,曲面样条函数法、有限差分直接法和有限差分间接法三种方法推算 的结果不仅相互之间很吻合, 与实验结果比较的误差也只有?y1=0.0068 和?y2=0.0055, 效果是 令人满意的。 表 2-2.5 是丙酮(1)-甲醇(2)-乙醇(3)三元系 101.325kPa 恒压数据的推算结果, 逸度系数 和过量焓均略去。由表可见,曲面样条函数法和有限差分直接法两种方法的推算值与实验值 符合的情况都令人满意,曲面样条函数法的推算结果更好。可见对于恒压系统,间接法对过 量焓确是不敏感,但直接法则有明显影响。
表 2-2.4 苯(1)-正己烷(2)-环己烷(3)三元系推算结果比较(70℃ ) x1 0.125 0.250 0.375 0.500 0.625 0.750 0.875 0.125 0.250 0.375 0.500 0.625 0.750 0.875 0.125 0.250 0.375 0.500 0.625 0.750 0.875 x2 0.2188 0.1875 0.1563 0.1250 0.0937 0.0625 0.0312 0.4375 0.3750 0.3125 0.2500 0.1875 0.1250 0.0625 0.6492 0.5565 0.4638 0.3710 0.2782 0.1855 0.0928 实验值 p/kPa 88.62 84.70 85.42 84.93 83.87 81.73 78.13 91.06 91.94 91.33 89.99 87.49 84.79 80.06 98.54 98.45 96.90 94.66 91.77 88.06 81.31 ?y y1 0.1523 0.2830 0.3775 0.4825 0.5740 0.6820 0.8225 0.1460 0.2640 0.3670 0.4540 0.5600 0.6575 0.8080 0.1290 0.2350 0.3420 0.4270 0.5325 0.6440 0.7890 y2 0.2860 0.2470 0.2125 0.1780 0.1490 0.1125 0.0625 0.5040 0.4360 0.3760 0.3290 0.2640 0.2030 0.1175 0.7020 0.6170 0.5375 0.4575 0.3810 0.2910 0.1730 曲面样条函数法 y1 y2 0.1365 0.2881 0.2678 0.2466 0.3781 0.2100 0.4763 0.1762 0.5781 0.1424 0.6820 0.1115 0.8222 0.0686 0.1327 0.5134 0.2500 0.4479 0.3503 0.3883 0.4479 0.3318 0.5593 0.2648 0.6622 0.2082 0.8081 0.1215 0.1214 0.7140 0.2373 0.6195 0.3288 0.5466 0.4337 0.4616 0.5341 0.3825 0.6420 0.2963 0.8013 0.1634 0.0068 0.0053 有限差分直接法 y1 y2 0.1415 0.2872 0.2664 0.2485 0.3770 0.2127 0.4785 0.1784 0.5792 0.1433 0.6899 0.1048 0.8245 0.0588 0.1327 0.5178 0.2496 0.4509 0.3538 0.3900 0.4514 0.3315 0.5515 0.2707 0.6659 0.2009 0.8093 0.1143 0.1245 0.7071 0.2354 0.6185 0.3354 0.5382 0.4304 0.4616 0.5295 0.3813 0.6458 0.2872 0.7957 0.1659 0.0068 0.0055 有限差分间接法 y1 y2 0.1414 0.2872 0.2664 0.2485 0.3769 0.2128 0.4785 0.1784 0.5793 0.1433 0.6900 0.1048 0.8246 0.0588 0.1326 0.5178 0.2495 0.4510 0.3538 0.3900 0.4515 0.3315 0.5516 0.2706 0.6659 0.2009 0.8093 0.1143 0.1244 0.7071 0.2354 0.6185 0.3355 0.5382 0.4347 0.4615 0.5296 0.3813 0.6458 0.2872 0.7956 0.1660 0.0068 0.0055

表 2-2.5 丙酮(1)-甲醇(2)-乙醇(3)三元系推算结果比较(101.325kPa) x1 0.017 0.019 0.011 0.060 0.060 0.059 0.113 0.113 0.115 0.190 0.186 0.176 0.283 0.255 0.379 0.357 0.491 0.469 0.638 0.619 0.768 0.899 x2 0.330 0.535 0.850 0.168 0.531 0.793 0.188 0.562 0.838 0.201 0.465 0.709 0.279 0.563 0.074 0.375 0.074 0.331 0.075 0.231 0.082 0.068 实验值 p/kPa 72.3 69.3 65.6 72.4 67.7 64.3 69.9 65.1 61.8 67.0 63.7 61.2 63.2 60.5 63.3 60.2 61.2 58.6 58.9 57.2 57.2 56.0 ?y y1 0.048 0.045 0.034 0.157 0.143 0.133 0.257 0.235 0.220 0.378 0.343 0.315 0.461 0.408 0.569 0.513 0.652 0.591 0.739 0.681 0.815 0.893 y2 0.434 0.645 0.882 0.233 0.585 0.785 0.236 0.572 0.761 0.229 0.457 0.630 0.298 0.505 0.093 0.350 0.095 0.316 0.097 0.246 0.103 0.080 曲面样条函数法 y1 y2 0.0454 0.4314 0.0467 0.6532 0.0257 0.8820 0.1597 0.2266 0.1408 0.6040 0.1248 0.7879 0.2756 0.2288 0.2290 0.5881 0.2266 0.7466 0.3910 0.2220 0.3447 0.4642 0.3020 0.6395 0.4786 0.2698 0.3937 0.5171 0.6046 0.0812 0.5192 0.3414 0.6762 0.0789 0.6051 0.2934 0.7667 0.0778 0.7198 0.2041 0.8178 0.0950 0.8895 0.0899 0.0122 0.0127 有限差分直接法 y1 y2 0.0484 0.4474 0.0483 0.6513 0.0244 0.8954 0.1686 0.2214 0.1393 0.6044 0.1239 0.7953 0.2793 0.2261 0.2338 0.5821 0.2191 0.7585 0.3987 0.2187 0.3508 0.4605 0.3132 0.6323 0.4873 0.2738 0.4162 0.4981 0.6113 0.0728 0.5308 0.3392 0.6883 0.0712 0.6105 0.2960 0.7715 0.0721 0.7107 0.2159 0.8326 0.0810 0.9067 0.0704 0.0177 0.0138

图 2-3 和 2-4 画出了用曲面样条函数法得到的二氯甲烷(1)-氯仿(2)-四氯化碳(3)三元系 在 45℃ 时,和丙酮(1)-氯仿(2)-乙醇(3)在 101.325kPa 下的 Q 函数曲面。前者比较简单,三 个二元系和三元系都表现出正偏差。后者相当复杂,其中氯仿-乙醇和丙酮-乙醇两个二元 系是正偏差,氯仿-乙醇二元系还有最低恒沸点,而丙酮-氯仿二元系则是一个负偏差系统, 并有一个最高恒沸点,三元系的 Q 函数曲面则呈现复杂多变的形状。对于这两个三元系,推 算都取得很好效果, 计算的气相组成与实验值的均方误差均在 0.01 左右, 迭代次数前者 4 次, 后者也只有 9 次。这应该说是对这一方法严峻的考验。更多的计算实例可以参考文献[2-24]。

图 2-3 二氯甲烷(1)-氯仿(2)-四氯化碳(3)

图 2-4 丙酮(1)-氯仿(2)-乙醇(3)三元系

三元系在 45℃ 下的 Q 函数曲面

在 101.325kPa 下的 Q 函数曲面

曲面样条函数法也已成功地应用于四元系的推算。对环己烷(1)-苯(2)-异丙醇(3)-甲乙 酮 (4) 在 101.325Pa 下的恒压数据,推算的气相组成误差分别为 ?y1=0.0133 、 ?y2=0.0134 、 ?y3=0.0133 和?y3=0.0095,迭代共 36 次。对乙醇(1)-氯仿(2)-丙酮(3)-正己烷(4)在 55℃ 下 的恒压数据,推算的气相组成误差分别为?y1=0.0144、?y2=0.0145、?y3=0.0288 和?y3=0.0201, 迭代共 25 次。详细结果参见文献[2-18]和[2-24],结果令人满意。 当实验得到的 T、p、x 数据有较大的误差时,会对推算的结果产生一定影响,特别是三 元或多元系统,其实验测定的误差更大一些,可以对原始实验数据进行适当的光滑化处理。 例如,对于三元恒温系统,可以对实验数据用下列经验式进行关联:
* * * p ? x1 p1 ? x2 p2 ? x3 p3 ? x1 x2 S12 ? x1 x3S13 ? x2 x3S23

? x1 x2 x3 ( D1 ? D2 x1 ? D3 x2 ? D4 x1 x2 ? D5 x1 x3 ? D6 x2 x3 ? D7 x1 x2 x3 )
Sij ? ? Ck ,ij ( xi ? x j )k
k ?0 n

(2-2.64)

(2-2.65)

* * 式中 pk 是纯组分 k 的饱和蒸气压,对于恒压系统,以系统压力下纯组分的沸点 Tk* 代替 pk 。

Ck,ij 由相应的二元系(i-j)的蒸气压或沸点数据关联而得。

2.3 气液平衡数据的热力学一致性检验
汽液平衡时有 T、p、x 和 y 四个变量,相律告诉我们其中只有两个是独立变量,如果实 验测定了第三个变量,则理论上可以根据热力学关系计算得到第四个变量,前面介绍的由 T、 p、x 数据推算 y 就是这种情况。另一方面,汽液平衡时的四个变量都是实验可以直接测定得 到的。则实验测得的四个变量必须符合热力学关系,这就是热力学一致性的要求。考察实验 得到的气液平衡数据是否符合热力学一致性要求,称为热力学一致性检验。 热力学一致性校验的基础是 Gibbs-Duhem 方程。根据式(1-7.41),对于二元系,活度因子 的 Gibbs-Duhem 方程可表示为
x1d ln ? 1,I ? x2d ln ? 2,I ? ?
E E Hm Vm T d ? dp RT 2 RT

(2-3.1)

如果忽略温度压力对 Q 函数的影响,则上式变为,

x1d ln ? 1, I ? x2d ln ? 2, I ? 0
上式对 x1 求导,得

(2-3.2)

x1

d ln ? 1,I d ln ? 2,I ? x2 ?0 dx1 dx1

(2-3.3)

理论上,我们可以直接用式(2-3.3)检验气液平衡数据是否符合热力学一致性,亦即根据 气液平衡实验数据,由式(2-2.7)计算各组分的液相活度因子,以 ln?1,I 和 ln?2,I 对 x1 作图,并获 取不同组成下曲线的斜率,然后代入式(2-3.3),看是否满足 Gibbs-Duhem 方程,这种检验方 法称为斜率法。看起来斜率法既简单又严格,但却不太有实用价值,因为要准确获取曲线的 斜率是比较困难的。因此,斜率法只能提供一种粗略的热力学一致性检验方法,只能作为定 性或半定量的方法应用。例如,在给定组成下,如果 dln?1,I/dx1 是正值,那么 dln?2,I/dx2 必须 是负值;如果 dln?1,I/dx1 等于零,则 dln?2,I/dx2 也必须等于零。因此,斜率法能方便地用来检 验实验数据中的严重误差。 2.2 节介绍的由 T、 p、 x 推算 y 的方法也可用来检验气液相平衡实验数据的热力学一致性。 如 van Ness 等[2-27]提出,由 T、p、x 推算得到的 y 与实验值比较,如果 ?y ? 0.01 满足,即通 过校验,哪一点不满足即那一点有问题,如有相当多的数据不满足则整个数据系列有疑问。 也可以用 y 的推算值与实验值的平均误差来判断数据的质量,例如 ?y ? 0.01。这一方法还可 以应用于多元系。但以 ?y ? 0.01作为热力学一致性的判断标准,无法反映?y 的离散情况及 y 是否存在系统偏差。 下面介绍几种常用的热力学一致性检验方法。 2.3.1 面积检验法 一种比较简单而有效的定量检验气液平衡数据热力学一致性的方法是 Herington[2-28]及 Redlich 和 Kister[2-29]发展起来的面积检验法。 根据过量吉氏函数 Q 与活度因子的关系,式(2-2.8),对于二元系有

Q ? x1 ln ? 1,I ? x2 ln ? 2,I
在恒温、恒压下对 x1 求导,得

(2-3.4)

d ln ? 1, I d ln ? 2, I dQ ? x1 ? ln ? 1, I ? x2 ? ln ? 2, I dx1 dx1 dx1
将式(2-3.3)代入上式,得
dQ ? ln(? 1, I / ? 2, I ) dx1

(2-3.5)

(2-3.6)

对 x1 从 x1=0 到 x1=1 积分,并注意到 Q(x1=0)=0 和 Q(x1=1)=0,得

?

x1 ?1 dQ ? dx1 ? ? ln 1 dx1 ? 0 x1 ? 0 dx x1 ? 0 ?2 1 x1 ?1

(2-3.7)

上式表明,如以 ln(?1,I/?2,I)对 x1 作图,在 x1=0~1 的区间内曲线与 x1 轴所包面积应等于零。典 型曲线如图 2-5 所示,曲线与 x1 轴所包面积等于零,这就意味着,如果根据某组气液平衡数 据计算得到各组分的活度因子,并以 ln(?1,I/?2,I)对 x1 作图,当 x1 轴上方的面积(面积 A)等于 x1 轴下方的面积(面积 B)时,该组数据是满足热力学一致性要求的。

图 2-5 30~35℃ 范围内,乙醇(1)/甲基环己烷(2)二元系的活度系数比值随组成的变化关系

由于实验数据总有一定的误差, 通常 x1 轴上方的面积不会严格等于 x1 轴下方的面积。 存 在一定大小的误差应该是允许的。一般情况下,如果

D?

(面积A) ? (面积B) ? 100 ? 2 ~ 5 (面积A) ? (面积B)

(2-3.8)

即可认为该组气液平衡实验数据是符合热力学一致性的。 式(2-3.6)是在恒温、恒压条件下推导得到的。但实际系统的气液平衡都是在恒温或者恒 压条件下获得的, 这时需要考虑温度或压力变化对 Q 函数的影响。 将式(1-7.41)代入式(2-3.5), 得
E dQ H E dT Vm dp ? ln(? 1, I / ? 2, I ) ? m2 ? dx1 RT dx1 RT dx1

(2-3.9)

对 x1 从 x1=0 到 x1=1 积分,得
E E x ?1? H ?1 Vm dp ? m dT ?x ?0 ln ? 2 dx1 ? ?x ?0 ? ? RT 2 dx ? RT dx ? ?dx1 1 1? ? x1 ?1
1 1 1

(2-3.10)

E E 对于恒温数据, H m 项消失,通常 Vm 很小可以忽略不计,式(2-3.8)仍可应用。如为恒压数据, E E 项消失,而 H m 通常是不能忽略的。但通常情况下实验条件下的过量焓数据难以获得,可 Vm

以采用下面的方法检验。先由式(2-3.8)计算 D 值,然后与另一数量 J 比较。J 由下式计算:

J ? 150 ? / Tm , ? ? T1 ? T2

(2-3.11)

Tm 为整个组成范围内的最低沸点,T1 和 T2 分别为组分 1 和 2 的沸点温度。如果有恒沸点,则 可按下图计算 ?。

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

如果(D-J)<10,则一般可以认为实验数据是符合热力学一致性的。 面积校验法简单易行,缺点是缺乏严密的误差分析,使结果带有一定的任意性。例如, 面积检验法由于采用活度因子的比值,即
* * *L * * * *L * ? 1,I / ? 2,I ? y1?1 p2 ?2 x2 exp[ Vm Vm , 2 ( p ? p2 ) / RT ] /?y2?2 p1?1 x1 exp[ ,1 ( p ? p1 ) / RT ]?

(2-3.12)

因为低压下 Poynting 因子和气相逸度因子对活度因子的影响很小,所以面积检验法对压力测 量误差很不敏感。van Ness 指出[2-30],面积检验法作图所需的数据仅仅是 x 和 y 数据以及两个
* * * * 纯组分的蒸气压之比 p2 。对于恒温数据,面积检验法几乎只是判明蒸气压之比 p2 是 / p1 / p1

否适合于这组测定的 x-y 值,以 ln(?1,I/?2,I)对 x1 作图对于 x-y 数据的离散性也许是很敏感的, 但是它没有告诉我们有关这些数据间的内在一致性。 2.3.2 局部面积检验法 我们在 2.2 节已经指出,任何严格符合热力学一致性的一组 N 个结点上的 T、p、x、y 数据必然满足 Gibbs-Duhem 方程。但是,从实验误差原理分析,由于实验设备、物料纯度以 及操作熟练程度的限制,使 T、p、x、y 的实验测定不可避免地都带有一定误差,包括随机误 差和系统误差,因此这种遵守不是绝对的。由于存在误差,因此应将 T、p、x、y 的实验数据 看成是随机变量,通常可以假设它们符合正态分布,我们可以严格使用统计误差分析理论检 验 T、p、x、y 实验数据的热力学一致性,这样才能使一致性校验避免任意性。局部面积检验 法就是在这样的思想下建立起来的。它首先由 Mcdermott 和 Ellis 提出 [2-31] , Ulichson 和 Stevenson[2-32]进一步发展了这种方法,Dohnal 和 Fenclova[2-33]对这一方法作了深入研究,使其 更趋完善。 如果忽略温度压力对 Q 函数的影响,由 T、p、x、y 实验数据计算得到的二元系活度因

子应满足式(2-3.2)。 在相邻的两实验点[xi, xi+1]间对式(2-3.2)进行积分, 并使用梯形积分公式得,
Fi ? 1 1 ( xi ?1 ? xi )(ln ? 1,i ?1 ? ln ? 1,i ) ? (2 ? xi ?1 ? xi )(ln ? 2,i ?1 ? ln ? 2,i ) 2 2

(2-3.13)

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

(2-3.14)

2 其中 σ F 是 F 的协方差矩阵,可以根据误差传递原理根据式(2-3.13)由 T、p、x、y 实验数据的

2 2 2 标准误差 σT 、 σ2 p 、 σ x 和 σ y (由实验者提供)计算得到。

图 2-7 自由度为 6 的?2 分布概率密度

图 2-7 画出了自由度为 6 的?2 分布概率密度曲线,?2 的期待值即自由度,对图 2-7 的?2 分布,

?2 的期待值即为 6。 图中纵坐标为概率密度, 由图可见,? 2 ? ?12?? (6) 的概率为?,? 2 ? ?12?? (6)
的概率为(1-?)。一般统计学手册中均有列出在各种自由度下,概率 Pr [ ? 2 ? ?12?? ( N )] ? 1 ? ? 的 由于 F 的二次型 Q 服从自由度为 N 的?2 分布, Q 有可能具有 0 ~ ? 的各种数值, ?12?? ( N ) 值。 但是由图 2-7 可见,Q(即?2)很大时的概率密度是很小的。如果我们取概率为 0.95,?=0.05,

?称为臵信度(或显著性水平),由统计学手册可查出 1-?=0.95 时的 ?12?? ( N ) ,下式成立,
Pr [Q ? ?12?? ( N )] ? ?
(2-3.15)

即 Q 值比 ?12?? ( N ) 还要大的概率为?即 5%,而 Q 比 ?12?? ( N ) 小的概率为 1-?即 95%。一致性
2 校验即可根据这一特点进行。 首先根据实验得出的 T、 p、 x、 y 数据计算 F 及其协方差矩阵 ? F ,

并计算 Q 值,然后和查得的 ?12?? ( N ) 比较。如果 Q ? ?12?? ( N ) ,则在臵信度为?的条件下,实 验数据满足热力学一致性;如果 Q ? ?12?? ( N ) ,落入拒绝域,则在臵信度为?下,数据不能通

过校验。?通常即取 0.05,表示校验考虑了 95%的可能性。 上述方法可以从总体上判断一组实验数据的热力学一致性或实验误差的大小。但统计量 Q 中既包含了测量变量的随机误差,也包含了它们的系统误差(如果存在的话),单靠 Q 统计 量是无法分辨出这两种误差的。 为了检验测量变量是否存在系统误差, Dohnal 和 Fenclova[2-33] 提出将 Fi 分成 i 为单数(A)和 i 为双数(B)的两组。这样,每组中的 Fi 是相互独立的,则 Fi 的 加权平均值 Fw 是 F 的无偏估计。
Fw ?

??

M i ?1

Fi / S 2 ( Fi ) / ?i ?11 / S 2 ( Fi )
M

??

?

, w ? A, B

(2-3.16)

而 Fw 的方差的无偏估计则为
S 2 ( Fw ) ?
M M 1 ( Fi ? Fw ) 2 / S 2 ( Fi ) / ?i ?11 / S 2 ( Fi ) ? i ?1 M ?1

?

??

?

, w ? A, B

(2-3.17)

统计量 tw 可以检验 Fi 的均值是否为零,

tw ? Fw / S (Fw ) ~ t (M ? 1) , w ? A, B
而统计量 hw 则可以用于检验 Fi/S(Fi)的方差是否为 1,
hw ? ?i ?1 ( Fi ? Fw ) 2 / S 2 ( Fi ) ~ ? 2 ( M ? 1) , w ? A, B
M

(2-3.18)

(2-3.19)

在一定的臵信水平?下,如果能通过式(2-3.15)、(2-3.18)和(2-3.19)五种检验,则可以认为所测 定的该组 T、p、x、y 实验数据不存在系统偏差,随机误差在实验者所估计的误差范围内。 Dohnal 和 Fenclova[2-33]通过对参数敏感性的检验发现,统计量 tw 对测量变量是否存在系 统误差很敏感,而对随机误差很不敏感,因此它非常适合于用来检验系统误差;而统计量 Q 和 hw 则对随机误差比较敏感, 因此可以用于检验测量变量的随机误差是否与实验者给出的结 果一致。 虽然上述方法比其它方法都优越,但用梯形近似有很大缺点,因为当两相邻点的位臵各 向相反方向移动时,梯形面积不变,此外梯形近似以直线代替曲线,在实验数据不多时,将 引入较大失真。 2.3.3 统计检验法 2.3.3.1 方法的建立 统计检验法是从逸度表达的 Gibbs-Duhem 方程出发建立起来的。对于一组实验测定得到 的二元系 T、p、x、y 数据(共有 N 个点),如果它们是严格符合热力学一致性的,则必须符合 式(2-2.56)所示的 Gibbs-Duhem 方程。但由于随机误差的存在,T、p、x、y 应视作随机变量, 而作为它们的函数的 F 也必然是一个随机变量,理论上其数学期望应为零。对于第 i 个结点, 式(2-2.56)可表示为

? x 1 ? xi ? ? ln ?1 ? ? ? ln ? 2 ? ?? dy ? ? ? Fi ? ? i ? ? xi ? ? ( 1 ? x ) i ? ?y ? ? ?y ? ? ? ? dx ? ? ? ?i ? ?i ? ? yi 1 ? yi ?? ?i
E ? Hm ? xi L1,i ? (1 ? xi ) L2,i ? ? ln ?1 ? ? ? ln ? 2 ? ?? d T ? ? ? ,i ? x ? ( 1 ? x ) ? ? ? ? ?? ? i i RTi 2 ? ?T ?i ? ?T ?i ?? d x ?i ? *( L ) *( L ) ? 1 xiVm ? ? ln ?1 ? ? ? ln ? 2 ? ? ? d p ? ,1, i ? (1 ? xi )Vm , 2 , i ?? ? ? ? xi ? ? ( 1 ? x ) i ? ? ?p ? ? ?p ? ? ?? dx ? RTi ? pi ? ? ?i ? ?i ? ?i ? ?

(2-3.20)

?0

上式可简单表示为

? ? dy ? ? d T ? ? d p ? E *( L ) *( L ) Fi ? F ? ? , H m,i , L1,i , L2,i , Vm,1,i , Vm, 2,i ? Ti , pi , xi , yi , ? dx ? , ? d x ? , ? ? ?i ? ?i ? d x ?i ?
为简单计,定义

? ? ??0 ?

(2-3.21)

F ? ?F1, F2 , ..., FN ? T ? ?T1, T2 , ..., TN ?
T

T

(2-3.22) (2-3.23)

p ? ? p1, p2 , ..., pN ? x ? ?x1, x2 , ..., xN ?
T

T

(2-3.24) (2-3.25) (2-3.26)
T

y ? ? y1, y2 , ..., yN ?

T

? ? dT ? ? dT ? ? dT ? ? Tx ? ? ? ? dx ? , ? dx ? , ..., ? dx ? ? ? ?1 ? ?2 ? ?N ? ?? ? ? dp ? ? dp ? ? dp ? ? px ? ? ? ? dx ? , ? dx ? , ..., ? dx ? ? ? ? ?N ? ? ? ?1 ? ? 2 ? ? dy ? ? dy ? ? dy ? ? yx ? ? , , ..., ? ? ? ? ? ? ? ? dx ? ? dx ? N ? ? ? ?1 ? dx ? 2
E E E H ? Hm ,1 , H m , 2 , ..., H m , N

(2-3.27)

T

(2-3.28)

T

(2-3.29) (2-3.30)

?

?

T

L ? ?L1,1, L2,1, L1, 2 , L2, 2 , ..., L1, N , L2, N ,

?

T

(2-3.31)

*( L ) *( L ) *( L ) *( L ) *( L ) *( L ) V ? Vm ,1,1 , Vm , 2,1 , Vm ,1, 2 , Vm , 2 , 2 , ..., Vm ,1, N , Vm , 2 , N

?

?

T

(2-3.32)

则式(2-3.21)可表示为

F ? F?T, p, x, y, yx, Tx, px, H, L, V ? ? 0
2 F 的协方差矩阵 σ F 可以利用误差传递规则计算

(2-3.33)

2 2 T 2 T T 2 T 2 T 2 T σF ? FTσ T FT ? Fpσ p Fp ? Fxσ 2 x Fx ? Fy σ y Fy ? Fyx σ yx Fyx ? FTx σ Tx FTx 2 T 2 T 2 T 2 T ? Fpxσ p x Fpx ? FH σ H FH ? FLσ L FL ? FV σ V FV

(2-3.34)

上式中,
2 是 T 的协方差矩阵,FT 是 F 对 T 的偏导数,均为 N?N 阶对角阵; σT

2 是 p 的协方差矩阵,Fp 是 F 对 p 的偏导数,均为 N?N 阶对角阵; σp

σ2 x 是 x 的协方差矩阵,Fx 是 F 对 x 的偏导数,均为 N?N 阶对角阵;

σ2 y 是 y 的协方差矩阵,Fy 是 F 对 y 的偏导数,均为 N?N 阶对角阵;
σ2 H 是 H 的协方差矩阵,FH 是 F 对 H 的偏导数,均为 N?N 阶对角阵;
2 是 L 的协方差矩阵,FL 是 F 对 L 的偏导数,均为 N?2N 阶对角阵,其子块为 1?2 的矩阵; σL

FV 是 F 对 V 的偏导数,均为 N?2N 阶对角阵,其子块为 1?2 的矩阵; σ2 V 是 V 的协方差矩阵,

σ2 yx 是 yx 的协方差矩阵,为 N?N 阶方阵,Fyx 是 F 对 yx 的偏导数,为 N?N 阶对角阵;
2 σT x 是 Tx 的协方差矩阵,为 N?N 阶方阵,FTx 是 F 对 Tx 的偏导数,为 N?N 阶对角阵;

2 σp x 是 px 的协方差矩阵,为 N?N 阶方阵,Fpx 是 F 对 px 的偏导数,为 N?N 阶对角阵;

上面我们是将(dy/dx)、(dp/dx)和(dT/dx)作为独立的随机变量来考虑的,但实际测量的是 T、p、x、y 数据,各结点上的这三个一阶导数需采用数值方法根据 T、p、x 和 y 数据计算得 到。在胡英等[2-14,16,34]发展的整体统计校验法中,采用的是三点抛物线法,刘洪来[2-24]则采用 了三次样条函数法,效果更好,下面简要介绍这一方法。根据三次样条函数法[见附录 2-3 式 (C-4)],y 对 x 的一阶导数可表示为,

1 y ? yi M i ?1 ? M i ? dy ? ? hi ? ? ? ? hi M i ? i ?1 2 hi 6 ? dx ?i
其中 hi=xi+1-xi,Mi 和 Mi+1 由式(C-6)计算得到。用矩阵形式可表示为,
yx ? f (y, x)

(2-3.35)

(2-3.36)

根据误差传递原理,有
2 T 2 T σ2 yx ? f xσ x f x ? f yσ y f y

(2-3.37)

2 2 fx 是 f 对 x 的偏导数,fy 是 f 对 y 的偏导数,均为 N?N 阶矩阵。 σ T x 和 σ px 可以采用同样方法

求得。

由于 F 可看作 N 维随机变量,理论上它的数学期望是零向量。可以采用随机现象的数学 模拟(Monte Carlo)方法证明,F 服从均值为零的 N 维正态分布,按统计学原理,F 的二次型 Q 应服从自由度为 N 的?2 分布,即,
2 ?1 Q ? FT (σF ) F ~ ? 2 (N )

(2-3.38)

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

H0 : ? ? ? F

(2-3.39)

以 Q 为检验统计量,则在 H0 成立的条件下有式(2-3.15)。所以,在臵信水平?下的拒绝域为

Q ? ?12?? ( N )

(2-3.40)

当实验数据存在系统误差时,显然也将引起 F 对零的系统偏离,这时,F 的数学期望将 不再为零,不仅测量的随机误差对 Q 有贡献,其系统偏差也有影响。当我们用上述方法确定 了某组 T、p、x、y 数据不符合热力学一致性时,我们是无法判断到底是测量变量的系统误差 还市其随机误差所致的。对实验数据的使用者来说,系统误差和随机误差的意义是不同的, 若证明实验数据存在系统误差,则在未搞清是哪个测量变量存在系统误差及其大小前,使用 这组数据需要非常小心。若实验数据无系统误差而仅存在随机误差,尽管这种误差比较大, 通过合适的模型对其进行光滑化后仍可使用。因此我们有必要检验实验数据是否存在系统偏 差,即检验,

H 0 : ? ? ?F
根据统计学原理,统计量

(2-3.41)

t ? F / S ( F ) ~ t ( N ? 1)
可以用来检验 F 的数学期望是否为零,这里 F 是 F 的加权平均值,
F?

(2-3.42)

??

N

i ?1

Fi / S 2 ( Fi ) / ?i?11 / S 2 ( Fi )
N

??

?
?? ?

(2-3.43) (2-3.44)

S 2 (F ) ?

N N 1 ( Fi ? F ) 2 / S 2 ( Fi ) / ?i ?11 / S 2 ( Fi ) ? i ?1 N ?1

?

2 2 其中 S2(Fi)是 F 的协方差矩阵 σ F 中对角线上的元素,即 S 2 ( Fi ) ? σF ii 。而统计量 h 则可以用于

检验 F 的方差与根据指定的测量变量的方差由式(2-3.34)计算的 F 的方差是否一致。
h ? ?i ?1 ( Fi ? F ) 2 / S 2 ( Fi ) ~ ? 2 ( N ? 1)
N

(2-3.45)

根据以上结论, 气液平衡数据热力学一致性检验的计算过程是:对于一组实验得到的 T、
2 p、 x、 y 数据, 根据实验条件估计各个测量变量的方差, 由误差传递原理计算 F 的协方差 σ F 及

统计量 Q、t 和 h,在臵信水平?的条件下,若下列三式同时得到满足,则认为它们是符合热 力学一致性的。

Q ? ?12?? ( N )

(2-3.46) (2-3.47) (2-3.48)

t ? t1?? ( N ? 1)
h ? ?12?? ( N ? 1)

反之则认为不符合热力学一致性。若 Q 和 h 满足而 t 不满足,则说明该组数据存在系统误差。 因为 F 和 S ( F ) 均为加权平均值,随机误差的大小对其计算结果影响应该很小,实际数据的计 算结果也证明了这一点。但若 Q 和 h 未通过检验而 t 通过了,则说明该组数据不存在系统误 差,但实验的真正误差比估计的要大,可以改变给定的各个测量变量的方差 (增大方差)重新 计算,只要在可忍受的误差范围内能通过 Q 和 h 检验,我们仍可以认为它们是符合热力学一 致性的,但这只是测量变量存在较大随机误差的意义下的热力学一致性。 2.3.3.2 随机误差等级的确定 对于实验数据的使用者来说,除了少数数据可以在原始文献上找到各测量变量的误差外, 对大部分数据的误差都是不了解的 (各种数据集中也鲜有给出原始数据误差的),因此无法用 如上所说的根据实际的标准方差进行推断。一种替代的办法是根据现有实验装臵的测量精度 情况,为 ? T 、 ? p 、 ? x 和 ? y 规定若干个误差等级,如表 2-3.1 所示,用上述方法检验实验数 据能在哪一个误差等级下通过热力学一致性检验, 以此来判断它的质量。I 级表示实验精度很 高,随级数增加,误差依次增大,最后第 V 级误差很大,实际上是一个不能被接受的等级, 在这一级被通过或不被通过即可认为不满足热力学一致性。
E 除了 T、p、x、y 的标准误差 ? T 、 ? p 、 ? x 、 ? y 外,还需要为过量焓 H m 、汽化热 Lm,i(恒

*( L ) 压系统)和液体摩尔体积 Vm 由于它们并非与 T、 p、 x、 y 同时测量, , i (恒温系统)规定标准误差。

而是单独测量或是根据各种经验或半经验方法计算得到的,其误差情况比较复杂。可以采用 如下标准: 对于过量焓,可取
E ? H ? 0.05H m E E 如果没有 H m 的实验数据,以零代替, ? H 则可以通过过量自由焓 Gm 计算,

(2-3.49)

E ? H ? 1.5Gm

(2-3.50)

对于气化热 Lm,k,通常可以由饱和蒸气压的 Antoine 方程计算得到,标准误差可取为

? L ? 0.005L
*( L ) 对于恒温系统,液体摩尔体积 Vm , k 对 F 的影响很小,其标准误差可取为
*( L ) ?V ? 0.1Vm

(2-3.51)

(2-3.52) 表 2-3.1 T、p、x、y 的误差等级
等级 I II III IV V σx 0.00625 0.00125 0.0025 0.005 0.01 σy 0.00125 0.0025 0.005 0.01 0.02 σT/K 0.025 0.05 0.1 0.2 0.4 σp/Pa 33.375 66.75 133.5 267 534

2.3.3.3 误差相关性的考虑 以上考虑的标准误差纯粹是由测量上的随机误差引起的,例如 x 和 y 的误差由分析仪器 的精度决定,?T 决定于温度计的精度,?p 则决定于压力表的精度,它们之间是相互独立的。 但 T、p、x 和 y 本身相互间不是独立的,实验中一个变量的波动会引起其它量的变化,例如 恒温条件下,当 x 波动时,p 和 y 都将发生变化,温度的波动同样也会引起 p 和 y 的改变,相 互间的这种影响取决于其导数的大小,与具体的系统有关。 在对一组实际测定的 T、p、x、y 数据进行热力学一致性检验时,除了要承认测量仪器 误差的客观存在外,由于 T、p、x 和 y 间的不独立性而引起的相关性误差也应该是容许的。 基于上述考虑,它们的标准误差应作如下修正[2-34]:
2 02 ?x ??x

(2-3.53)

对于恒温系统

? ??
2 y

02 y

? ?y ? 02 ? ?y ? 02 ? ? ? ?T ? ? ? ?x ? ?T ? ? ?x ?

2

2

(2-3.54) (2-3.55)

2 02 ?T ? ?T

? ??
2 p

02 p

? ?p ? 02 ? ?p ? 02 ? ? ? ?T ? ? ? ?x ? ?T ? ? ?x ?

2

2

(2-3.56)

对于恒压系统

? ??
2 y

02 y

? ?y ? 02 ? ?y ? 02 ?? ? ?p ? ? ? p ? ? ?x ? ? x ? ? ? ?

2

2

(2-3.57)

? ??
2 T

02 T

? ?T ? 02 ? ?T ? 02 ?? ? ?p ? ? ? p ? ? ?x ? ? x ? ? ? ?

2

2

(2-3.58) (2-3.59)

2 ?p ? ? 02 p

上述各式中的 T、p、y 对 x 的一阶导数已由前面的三次样条函数法得到,其余的一阶导数可 以由下列各式计算,为简单计,推导时作了气相为理想气体的假设。

? ?p ? ? ?? ? ?T ?

* * ? d ln p1 ? d ln p2 ? ? p? y ? (1 ? y) ? d T d T ? ?

(2-3.60)

* * ? d ln p1 ? d ln p2 ? ?y ? ? ? ? y ( 1 ? y ) ? ? ? ? dT ? ? ?T ? ? dT ?

(2-3.61)

? ?y ? y ? ? ?p ? ??? p ? ? ? ?T ? ? ?p ? ? ? ?p ? ? ? 1 /? ?T ? ? ? ? ?
2.3.3.4 F 函数正态分布特性的检验

(2-3.62)

(2-3.63)

上述热力学一致性检验的统计学方法是建立在共存方程的偏差函数 F 服从 N 维正态分布 这一基本假设基础上的。这一假设的正确与否直接关系到上述方法检验的结果的可信程度。 下面我们检验这一假设的正确性。 由于 F 是 T、p、x、y 的复杂函数,无法从理论上根据 T、p、x、y 的正态分布特性导出 F 的理论分布的解析表达式,为此我们采用随机现象的数学模拟方法(Monte Carlo 方法)来进行 研究。 广义地说,凡是利用随机数对模型系统进行模拟以产生数值形式的几率分布的方法,都 可以称为 Monte Carlo 方法。 作为计算数学的一个分支, 它通常被用来解决一些用常规方法不 能解决的因素复杂的问题。这里我们要研究的是共存方程的偏差函数 F 的分布函数,它是 N 组 T、p、x、y 的复杂函数,而后者均服从正态分布。

F ? F?(T , p, x y)1, ..., (T , p, x y) N

?

(2-3.64)

假如 F 是服从 N 维正态分布的,则它的二次型应服从参数为 N 的?2 分布,即式(2-3.38)。显 然如果我们能证明式(2-3.38)成立,则 F 服从 N 维正态分布的假设也成立。检验的计算过程可 表示为: (1) 根据 T、p、x、y 的正态分布特性,在计算机上用随机抽样的方法,分别从 T、p、x、 y 的分布中抽样得 N 组 T、p、x、y 数据(随机数)。

2 (2) 根据抽样得到的 N 组 T、p、x、y 数据,由前述方法计算 F、 σ F 并得到 Q 的一个随

机数 Q?(T , p, x y)1, ..., (T , p, x y) N ? 。 (3) 重复步骤(1)和(2)M 次,得 M 个随机数 Q1、Q2、…、QM。 (4) 以 Q 的样本分布 Sn(Q)近似 Q 的分布。 如果假设是真的, 则 SM(Q)将能足够地逼近 Q。 我们的目的是要通过检验 Q 的?2 分布特性来证明 F 服从 N 维正态分布。即检验假设:

H 0 : F (Q) ? F ( ? 2 , N )
我们采用柯尔莫哥洛夫检验法检验[2-35]。 将统计量 Q 的定义域(0, ?)以?Q 为间隔分割,并定义样本的经验分布函数为

?0 S n (Q ) ? ? ?i / M

Q?0 n?Q ? Q ? (n ? 1)?Q n ? 1, 2, ..., ?

这里 i 是样本 Q 中小于或等于(n+1)?Q 的个数。 当样本数 M 较多时,Sn(Q)将是?2(N)的一个近似。取检验统计量 Dn,表示在 Q 的整个取 值范围内样本的经验分布函数与理论分布函数 F0(?2, N)的最大偏差,即,

Dn ? maxSn (Q) ? F0 ( ? 2 ? Q, N )
0?Q??

? maxSn ( n?Q) ? F0 ( ? 2 ? n?Q, N )
0?Q??

当 H0 成立时,统计量 M Dn 的极限分布函数为

? ? L ?2L2u 2 ) ? ? ( ?1) e xp( Q(u ? M Dn ) ? ?L??? ?0 ?

u?0 u?0

对于给定的显著性水平?,当 Dn>Dn(?)时便否定 H0。临界值 Dn(?)可查常用数理统计表。 采用上述方法,我们选择了 Benzene-Acetonitrile 系统在 25?C 下的恒温数据(由 T、p、x 数据根据 Mixon 有限差分法得到气相组成 y,N=19)作为无误差的标准数据。指定 T、p、x、 y 的方差为表 2-3.1 中方差等级的第 I 级,在计算机上产生方差为 ? i2 (i= T、p、x、y)的随机数 并分别加于 T、p、x 和 y 而得到一套含有误差的模拟实验数据,这样的模拟数据共产生 100 套。对于每套数据用前面的方法计算统计量 Q,得到容量为 M=100 的统计量 Q 的随机样本, 并作出样本的经验分布函数 Sn(Q),由前述方法得到的统计量 Dn(n=100)为 0.0824。图 2-8 是 经验分布函数 Sn(Q)和理论分布函数(?2 分布)的比较,可见两者是非常接近的。根据柯尔莫哥 洛夫检验法,在显著性水平?=0.05 下给出的拒绝域为

Pr(D100 ? D100,0.05 ) ? 0.05

由数理统计表查得 D100,0.05=0.13403,即,D100< D100,0.05,所以,可以认为原假设 H0 成立,即 Q 服从自由度为 N 的?2 分布,而 F 服从 N 维正态分布。

图 2-8 F 函数的理论分布与经验分布比较

氯仿-乙醇在 55?C 下的气液平衡数据是一套精度比较高的数据,用上述统计检验法对 该套数据进行热力学一致性检验的结果,Q、t 和 h 均在第 I 级上通过检验。各实验点上的偏 差函数 F 值如图 2-9 所示。但如果忽略气相非理想性,检验的结果发生很大的变化,Q 在第 II 级上通过,h 则在第 I 级上通过检验,而 t 则不能通过检验。不尽如此,对测量变量取不同 的随机误差,对统计量 t 几乎没有影响,此时各实验点上的偏差函数 F 值如图 2-10 所示。因 为忽略气相非理想性实际上意味着在实验数据中引入了系统偏差,上述结果说明,统计量 t 对系统误差非常敏感,而对随机误差不敏感,用它来检验一组气液平衡实验数据是否存在系 统误差是比较合适的。

图 2-9 氯仿-乙醇系统 55?C 下的 F 函数值

图 2-10 忽略气相非理想性时,氯仿-乙醇系统 55?C 下的 F 函数值

表 2-3.2 是对部分系统检验的结果,大部分数据取自文献[2-36]。为了比较,表中同时列 出了用面积检验法和 y-检验法检验的结果。从表中可以看到,当一组 T、p、x、y 数据不能通 过 t 检验时,并不一定 Q 和 h 都不能通过检验。相反,有些数据能在第 I 和 II 级上通过 Q 和

h 检验,却不能通过 t 检验。通常的情况是,如果不能通过 t 检验,则 Q 将比 h 在更低一级的 水平上通过检验,表明系统误差对 Q 比 h 有更大的影响。
表 2-3.2 二元气液平衡实验数据的热力学一致性检验

对大量实验数据检验的结果表明,统计检验法与 y-检验法的结果比较接近,而与面积检 验法的结果差别比较大。 任何实验不可能不存在误差。但不同的实验者由于其实验技能、操作的熟练程度不同, 所得实验数据的误差大小也有差异,如果所采用的仪器不存在系统误差,则实验误差就是一 种随机误差。如果不是到了难以忍受的程度,随机误差的存在即是容许的。在气液平衡中测 定的四个变量 T、p、x 和 y,它们并不是完全独立的,而受 Gibbs-Duhem 方程的限制。如果 这四个测量变量都不存在任何误差,则它们应严格遵循 Gibbs-Duhem 方程;如果只存在随机 误差,则它们也应该是随机地偏离这一限制方程;如果实验中发生了系统偏差,则得到的 T、 p、x 和 y 数据必然也系统地偏离该方程。热力学一致性检验的目的就是检验 T、p、x、y 实验 数据是否符合 Gibbs-Duhem 方程,在多大的误差范围内符合,它们是否存在系统误差? 总之,任何包含随机误差的数据是不可能完全地满足热力学一致性的。我们所要检验的 是在一定的显著性水平以及一定的误差等级下,数据是否满足热力学一致性。从统计学的观 点,测量变量是在它们的真值的周围随机地分布的,少数数据距真值较远,不能简单地认为 是不正常现象。除非超过某些明显的界限例如 3 ? ,我们才能说它是一个坏点而予以舍弃。 谈论个别实验点的一致性一般来说没有太大的意义;考虑到数据的随机性,而对整体作出判

断,才比较符合统计学的原理。这正是这一方法的要点。上面的检验是在一定的显著性水平 上的检验,例如 ? ? 0.05 ,表明即使作出判断,我们仍有一定的概率即 5%可能犯错误。统计 学的原理告诉我们,绝对可靠的东西是没有的,在一定的臵信水平上作出判断更合乎实际情 况。

附录 2-1 非线性方程组的 Broyden 解法[2-20] Broyden 法是一种改进的 Newton-Raphson 算法,其出发点是建立相邻两次迭代间 Jacobi 矩阵的逆矩阵的递推关系,从而可减少甚至避免求 Jacobi 矩阵及其逆矩阵。 非线性方程组一般可以采用 Newton-Raphson 迭代法求解,即
Xk ?1 ? Xk ? ?Xk

(A-1) (A-2)

?Xk ? [?J k ]?1 F(Xk )

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

[J k ]?1 ? [J k ?1 ]?1 ? ?[J k ?1 ]?1
其方法如下: (1) 假定变量 X 的初值为 X0,并计算,

(A-3)

F0 ? F(X0 )
(2) 计算 Jacobi 矩阵 J0 及其逆矩阵

(A-4)

H0 ? ?[J 0 ]?1
(3) 用 H 和 F 的最新值计算
?Xk ? H k F k

(A-5)

(A-6)

(4) 求取阻尼因子 wk,使得 F( X k ? wk ?X k ) ? F( X k ) : 首先令 wk=1,若满足下式则转入步骤(5)。
F( Xk ? wk ?Xk )
1/ 2

? F( X k )

1/ 2

(A-7)

否则,用下式计算 wk,

wk ? [(1 ? 6? )1/ 2 ? 1] / 3?

(A-8) (A-9)

? ? F ( X k ? wk ?X k ) / F ( X k )

如果经规定次数的搜索后,仍不能使欧几里得范数减小,则返回步骤(2),用 Xk 重 新计算 Jacobi 矩阵中的各个偏导数。 (5) 在步骤(4)的计算中,已求得
Xk ?1 ? Xk ? wk ?Xk

(A-10) (A-11)

Fk ?1 ? F(Xk ?1 )
检验 Fk+1 是否满足收敛条件,若不满足,则计算
Zk ? F k ?1 ? F k

(A-12) (A-13)

H k ?1 ? H k ?

(H k Zk ? wk ?Xk )(?Xk )T H k (?Xk )T H k Zk

然后返回步骤(2),直至满足收敛要求。 以上计算过程中,步骤(4)是可以灵活改变的,可以采用更为有效的一维搜索方法。

附录 2-2 非线性微分方程的 Runge-Kutta 法积分求解[2-37] 设微分方程为
dy ? f ( x, y ) dx

(B-1)

从某初始值(x0, y0)出发,为求下一点 x1=x0+?x 时,y 值的增量?y,可按下列步骤计算:

?y' ? ?xf ( x0 , y0 ) ?y' ' ? ?xf ( x0 ? ?x / 2, y0 ? ?y' / 2) ?y' ' ' ? ?xf ( x0 ? ?x / 2, y0 ? ?y' ' / 2) ?y' ' ' ' ? ?xf ( x0 ? ?x, y0 ? ?y' ' ' )
?y ? (?y'??y' ' ' ' ) / 6 ? (?y' '??y' ' ' ) / 3

(B-2) (B-3) (B-4) (B-5) (B-6) (B-7)

y1 ? y0 ? ?y
得到(x1, y1)后,继续用上述方法计算(x2, y2),…,直至求解结束。

附录 2-3 三次样条函数[2-38] 对于区间[a, b]的给定划分?={xi}N: a=x0 < x1 < … < xN < xN+1=b, 可以唯一地得到满足如下 条件的分段三次多项式 S(x): (1) S(x)在每一子区间(xi , xi+1)上是三次多项式,(i=0, 1, …, N)

(2) S(x)? C2[a, b] (3) 在各结点上满足

S ( xi ) ? yi , i ? 0, 1, ..., N , N ? 1
'' '' S ' ' (a) ? y0 , S ' ' (b) ? yN ?1

(C-1) (C-2)

这样的分段三次多项式 S(x)成为三次样条函数,这里满足上述条件的 S(x)可写成

S ( x) ? M i

( xi ?1 ? x)3 ( x ? xi )3 M h ( x ? x) M h ( x ? xi ) ? M i ?1 ? ( yi ? i i ) i ?1 ? ( yi ?1 ? i ?1 i ) 6hi 6hi 6 hi 6 hi

(C-3)

( xi ? x ? xi ?1 )
相应的一阶和二阶导数为

S ' ( x) ? ? M i
S ' ' ( x) ? M i

( xi ?1 ? x)2 ( x ? xi )2 yi ?1 ? yi M i ?1 ? M i ? M i ?1 ? ? hi 2hi 2hi hi 6
xi ?1 ? x x ? xi ? M i ?1 hi hi ( xi ? x ? xi ?1 )

( xi ? x ? xi ?1 )

(C-4)

(C-5)

其中 hi=xi+1-xi, M i ? S ' ' ( xi ) ,Mi 和 yi 之间存在如下关系

? y ? y y ? yi ?1 ? hi ?1M i ?1 ? 2(hi ?1 ? hi )M i ? hi M i ?1 ? 6? i ?1 i ? i ? ? 0 , i ? 1, 2, ..., N ? 1, N hi ?1 ? ? hi

(C-6)

参考文献 [1] 郭天民,多元汽液平衡和精馏,化学工业出版社,1983 [2] 朱自强,姚善泾,金彰礼,流体相平衡原理及其应用,浙江大学出版社,1990 [3] Prausnitz J. M., Lichthenthaler R. N., de Azevedo E. G., Molecular thermodynamics of fluid-phase equilibria, 3rd ed., Prentice Hall PTR, 1999 [4] Assael M. J., Trusler J. P. M., Tsolakis T. F., Thermophysical properties of fluids, Imperial College Press, London, 1996 [5] Prausnitz J. M., Anderson T. F., Computer calculations for the multicomponent vapor-liquid and liquid-liquid equilibria, Prentice-Hall Inc., Englewood Cliffs. 1980 [6] Malanowski S., Fluid Phase Equilibria, 8, 197(1982) [7] Malanowski S., Fluid Phase Equilibria, 9, 311(1982) [6] Weir R. D., de Loos Th W., Measurement of the thermodynamic properties of multiple phases, Elsevier, Amsterdam, 2005

[9] Hala E., Pick J., Fried V., Vilim O., Vapour-liquid equilibria, 2nd Ed., Pergamon, Oxford, 1967 [10] Barker J. A., Austr. J. Chem., 6, 207(1953) [11] Abbott M. M., van Ness H. C., Fluid Phase Equilibria, 1, 3(1977) [12] Mixon F. O., Gomowski B., Carpenter B. H., Ind. Eng. Chem. Fundam., 4, 455(1965) [13] 胡英,英徐根,张鸿哲,化工学报,(2), 153(1979) [14] 胡英,流体的分子热力学,高等教育出版社,上海,1983 [15] 刘洪来, 英徐根, 胡英,化工学报,42, 393(1991) [16] 胡英,近代化工热力学,上海科学技术文献出版社,1994 [17] 刘洪来, 英徐根, 胡英, 化工学报,42, 400(1991) [18] Hu Y., Liu H.L., Prausnitz J. M., Fluid Phase Equilibria,95, 73(1994) [19] Harder R. L., Desmarais R. N., J. Aircraft, 9, 189(1972) [20] Holland C. D., Fundamentals of multicomponent distillation, McGraw-Hill, New York, pp.147-149, 1981 [21] 胡英,英徐根,化工学报,(1),27(1980) [22] Ljunglin J. J., van Ness H. C., Chem. Eng. Sci., 17,531(1962) [23] 胡英,王琨,吕瑞东,化工学报,(4),341(1980) [24] 刘洪来,华东化工学院博士学位论文,1991 [25] 刘洪来,冯国祥,胡英,天然气化工,(1), 47(1987) [26] 刘洪来,冯国祥,胡英,天然气化工,(2), 48(1987) [27] van Ness H. C., Byer S. M., Gibbs R. E., AIChE J., 19,238(1973) [28] Herington E. F. G., Nature, 160,610(1947); J. Inst. Petrol., 37,457(1957) [29] Redlich O., Kister A. T., Ind. Eng. Chem., 40, 345(1948) [30] van Ness H. C., Pure & Appl. Chem., 67, 859(1995) [31] Mcdermott C., Ellis J. R. M., Chem. Eng. Sci., 20, 293(1965) [32] Ulichson D. L., Stevenson F. D., Ind. Eng. Chem. Fundam., 11, 287(1972) [33] Dohnal V., Fenclova D., Fluid Phase Equilibia, 21, 211(1985) [34] 胡英,英徐根,唐小琪,华东化工学报,(2),205(1984) [35] 戚国安,数理统计,见《现代工程数学手册》第 IV 卷第六十篇,华中工学院出版社, 武汉,pp.107-240,1987 [36] Gmehling J., et al., Vapor-liquid equilibrium data collection, part 1-6, DECHEMA, 1977-1986 [37] 南京大学数学系计算数学专业编,常微分方程数值解法,科学出版社,北京,1979

[38] 孙家永日,样条函数与计算几何,科学出版社,北京,1982


赞助商链接
相关文章:
气液平衡的计算
气液平衡的计算方法 气液平衡的计算方法摘要 本文综合分析了多组分相平衡理论...- 8 - -2- 气液平衡的计算方法 1 前言 精馏计算是空分流程计算的关键环节,...
高一必修二化学平衡的简单计算
1 = 1 n 2 P2 n 2 V2 计算的技巧 1、平衡状态的判断例 1.可逆反应 N2+3H2 A.3υ 正(N2)=υ正(H2) C.2υ 正(H2)=3υ 逆(NH3) 2NH3 的...
2相平衡
相平衡练习题 一、是非题,下列各题的叙述是否正确,对的画√错的画× 1、 ...( 10、在相图中总可以利用杠杆规则计算两相平衡时两相的相对量。( 二、选择题...
第5章 相平衡
? 上式为以化学位表示的相平衡判据式。化学位在计算上不方便,故在解决相平衡问题 时用组分逸度代替化学位,可得: ?? ? f ?? f i i (i=1,2,3,…,N...
相平衡二
相平衡 4页 免费如要投诉违规内容,请到百度文库投诉中心;如要提出功能问题或意见...相平衡二 相平衡相平衡隐藏>> 物理化学试卷班级 姓名 分数 一、计算题 ( 共...
化学平衡计算(带答案)
计算: (1)平衡时混合气体中 H2 的体积分数是多少? (2) N2 的转化率为多少? 解析:(1)66.7%(2) 20% 13.在一定条件下已测得反应 2CO2===2CO+O2 ...
气液相平衡方面的基础知识
气液相平衡 相是指系统的某一部分具有相的物理和化学性质,具有相同的组成,并且...闪蒸计算前需要确定的是: 1>气相的摩尔分数 n v ; 2>液相的摩尔分数 n l...
经典相平衡习题
经典相平衡习题 - 一、选择题 ( 共 11 题 20 分 ) 1. 2 分 (2738) 2738 二级相变符合的爱伦菲斯方程式是: (A) dp / dT ? ?CV /(TV??) (B)....
物理化学习题答案(2)-相平衡
物理化学习题答案(2)-相平衡 - 物理化学习题答案(2) 三. 计算题 2、固态氨的饱和蒸气压为 ln( p / kPa) ? 21.01 ? 3754 , T /K 液态氨的饱和蒸气...
张7.2.(1、2)相平衡关系拉乌尔定律
张7.2.(1、2)相平衡关系拉乌尔定律_化学_自然科学_专业资料。班级: 姓名: ...2. 【小组合作】完成例题其他 8 组数据的计算——各温度下平衡汽液相中苯的...
更多相关标签: