气液平衡的计算方法

合肥学院

Hefei University 《化工热力学》过程论文

题目：气液平衡的计算方法

系别：化学与材料工程系

专业：化学工程与工艺

学号：1303021001

姓名：于晓飞

教师：高大明

气液平衡的计算方法

摘要：气液平衡计算是化学过程中一项十分重要的计算。气液平衡的计算方法有几种，活度系数法，状态方程法（EOS 法），GEMC 和GDI 方法计算流体气液相平衡。在气液平衡的计算中有三种泡点计算 、露点计算和闪蒸计算，这里我们对闪蒸计算不做研究。

关键词:气液平衡 计算方法GEMC GDI 正文：

气液平衡计算的基本公式及计算类型：相平衡的判据应用于气液平衡，即为：

=

f

V i

?f L i

?

(i=1,2,3,…,N)式中，f i

?为混合物中组分i 的逸度；上标V 指的是气相；

上标L 指的是液相。

上式既是气液平衡的准则，有事气液平衡计算的基本公式。具体应用时，需要建立混合物中组分的逸度

f

V i

?、

f L i

?

与体系的温度、压力以及气液相平衡组成

关系。

1.1活度系数法

根据溶液热力学力论，将液相中组分的逸度与组分的活度系数相联系，简称活度系数法。对液相，由活度与活度系数的定义式得出

f L i

?

=

f

x i

i

i

θ

γ

式中，f i

θ

为标准态的逸度，以取Lewis-Randall 定则为基准的标准态，即纯液体i 在体系的温度下的逸度。

f

i

θ

=

f

L i

=

dp p

p RT

s i V p L

i S i

S i

?exp φ

式中，指数项dp p

p RT

s i V L

i ?exp 称为Poynting 因子，其意义是压力对

f

i

θ

影响的校正。

对气相

将f

L i

?与f

V i

?表达式带入式中，得

φ

?V i

p y i =dp p RT

p

L i

S i

S i

i

i

s

i

V

p x ?

exp

φγ (i=1,2，…，N)

式中，y i

和x i 分别为汽、液相中组分i 的摩尔分数；φ

?V i

为气相混合物中组分i 在体系温度T ，体系压力p 下的逸度系数；γi 为液相中组分i 的活度系数；p S

i

为

纯组分i 在体系温度T 时的饱和蒸气压；φS

i

为纯组分i 在体系温度T 与其饱和

蒸气压

p

S i

时的逸度系数；V L

i 为纯组分i 在体系温度T 时液相的摩尔体积。

1.2GEMC 方法计算原理

GEMC 方法可同时在两个模拟盒子中进行蒙特卡罗(MC)模拟,二者相对独立,但保持热力学相关,即满足相平衡条件(压力、温度和化学势相等),其温度T 、总体积V 和两个盒子中的总粒子数N 保持不变.为达到相平衡,在模拟过程中需要进行3种不同的蒙特卡罗移动,按不同的接受概率接受,以满足相平衡条件:

1)两个模拟盒的粒子分别在盒内自由移动,包括粒子的平动、转动等,以达到盒内平衡,其接受概率为

P=min[1,exp(-ΔU/kT)] （1） 2)在保持总体积不变的条件下,在两个模拟盒子间进行体积的涨落,以达到两模拟盒压力相等,其接受概率为

P=min ???

?????-++++-V V 22

2111ΔV ln ΔV ln ΔU)/KT ΔU (exp(,1V N V N （2）

3)保持粒子总数不变的条件下,在两个模拟盒

子间进行粒子交换,以达到两模拟盒化学势相等,接受概率为

P=min ???

?????+-+)/KT )ΔΔ(exp()1(,1U U 21211

2V N V N （3）

式中,U 、V 和N 分别为两盒子的能量、体积和粒子

数,T 为体系温度,k 为玻尔兹曼常数. 1.3 GDI 方法计算原理

流体的蒸发焓与饱和蒸气压可通过Clapeyron 方程关联起来:

V H V V

ΔΔT dT dp =

（4） 式中,p 、T 分别为体系压力和温度,ΔHv 为蒸发焓,ΔVv 为两相的体积差.对含气体的相平衡,(4)式一般变形为

V H V

V ΔΔln pT dT p

d = （5）

(5)式右边气液两相的蒸发焓及体积差均可由分子模拟求得,进而式转变为压

力对温度的常微分方程,可通过数值法求解.因此在已知体系气液共存线上一个点作为积分起点的情况下,整个气液共存线可通过分子模拟并结合热力学积分计算得到.本工作采用三种分子的沸点作为积分参考点,液体和气体分别采用分子动力学和蒙特卡罗方法模拟,通过预测-矫正法计算积分,最终计算气液共存线。 1.4Q 函数法(间接法)

汽液平衡时，按判据式(1-6.13)，)()(L V k k f f =(k =1, …, K )，如气相采用逸度因子、液相采用第I 种活度因子分别计算气液相的非理想性，得

]/)(exp[*

*,,**RT p p V x p py k L k m k k k k k k -=I γ??，K k ,,2,1 =

(2-2.1)

整理上式可得系统总压p ，

∑∑==-==K

k k k L k m k k k k K

k k RT p p V x p py p 1

*

*,,**1

/]/)(exp[?γ?I

(2-2.2)

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

??

??????????????? ????-???? ????-???????????? ????-???? ????+????????? ????-???? ????+-=∑∑∑∑-=-==-=11],[],[11],[],[2111],[],[*

*,**exp ]

/)(exp[K j K j x j j K k x k E m K j K j x j j K k x k E m K

k K j K j x j j K k x k k k L k m k k k x p x x p RT V x T x x T RT H x Q x x Q Q RT p p V x p p ??

(2-2.3)

注意当i =K ，式中对x K 的偏导数全为零。式(2-2.3)的意义在于：如果暂时不考

虑*k p 、*k ?、L k m V *,、k ?、E m H 和E

m V ，则式中除了Q 以外，其它的变量就是已输入

的T 、p 、x 。而Q 函数正是T 、p 、x 的函数，式(2-2.3)实质上是一个Q 函数的偏微分方程，只要有足够数量的一系列T 、p 、x 的实验数据，原则上可以解得

Q =Q (T , p , x )。有了Q ，可用式(1-7.46、47)计算γk ,I ，代入式(2-2.1)即可求得

y 。至于那些暂时放在一边的变量：其中*k p 、*k ?和L k m V *,是纯组分性质，与混合物

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

差。至于E m H 和E m V ，后者很小，常可忽略，前者对于恒温数据不起作用，对于恒

压数据，实践证明略去后影响不大。总之，这一方法基本上不使用模型，或者严格地说，不使用液相活度因子模型，而它是气液平衡计算中最关键的模型。式(2-2.3)原则可以求解，但实践上却有很大困难，因为导数出现在exp 中，是一个超越型的偏微分方程，没有解析解，只能通过数值方法求解。国内外学者已发展了多种方法，根据所采用数值方法的不同，可以分为几种类型：最早是Barker ]提出的方法，其核心是选择一个过量函数模型代入式(2-2.3)，利用一系列T 、p 、

x 的实验数据，拟合得到模型参数和Q 函数。第二种是Mixon 等]发展的有限差分法，它以差分来逼近式(2-2.3)中的导数，然后利用Newton 法迭代求得离散格点

上的Q 值。这种方法不依赖于任何过量函数模型，是严格的无模型法。它对二元系的计算非常成功，得到广泛的应用。但用于三元系时，收敛速度极慢，且求解过程不稳定。第三种是样条函数法，包括适用于二元系的三次样条函数法和适用于任意组分数的曲面样条函数法。特别是曲面样条函数法，它不仅能方便地用于二元系和三元系，也能成功地应用于多元系，更重要的是不同组分数的计算方法可以统一在一个框架下。大量实例计算表明，没有收敛的困难，不受多元系Q 函数曲面类型的限制。

下面介绍计算简单的Barker 法和适用于多元系的曲面样条函数法，至于其它方法，感兴趣的读者可以参考相关的文献和著作。 1.5Barker 法

Barker 法的核心是选择一个合适的过量函数模型，其中包括若干待定模型参数，代入式(2-2.3)后，利用一系列T 、p 、x 的实验数据，拟合得到模型参数，这就得到Q 函数。原则上，第4章中介绍的各种过量函数或活度因子模型均可以作为Barker 法的候选模型。但这些模型都是针对特定的对象而建立起来的，都有一定的适用范围，这就使得Barker 法的准确度受到所选模型可靠性和适用性的限制，对多元系问题会更突出。解决的办法是尽可能选用灵活性大的经验模型，例如对于二元系，可以采用如下的Redlich-Kister 型的经验Q 函数模型，

[

]

∑=-=N

j j j x x A x x Q 02121)(

(2-2.4)

相对应的活度因子为

[

][

]

∑∑=-=-+-=N

j j j N

j j j x x jA x x x x A x 01212102122I ,1)(2)(ln γ

(2-2.5)

[

][

]

∑∑=-=-+-=N

j j j N

j j j x x jA x x x x A x 01122101221I ,2)(2)(ln γ

(2-2.6)

根据T 、p 、x 实验数据的多少和计算精度的要求，可以选择不同的N 值。

式(2-2.4)是一个关于A j 的线性方程，如果已知不同组成下的Q 函数值，可以非常方便地采用最小二乘法关联得到N +1个A j 。具体计算时，我们可以采用如下的迭代过程：

(1) 假设气相为理想气体、液相为理想溶液，计算气相组成的初值y 0;

(2) 由状态方程计算气相逸度因子*

k ?和?k ，由式(2-2.7)计算各组分的液相

活度系数；

{}]/)(exp[/*

*,**,RT p p V x p py k L k m k k k k k k -=??γI

(2-2.7)

(3) 由式(2-2.8)计算各实验点的过量吉氏函数Q ;

∑==K

k k k x Q 1I ,ln γ

(2-2.8)

(4) 由计算得到的过量吉氏函数Q 关联式(2-2.4)中的未知参数A j (j =0, …,

N );

(5) 由式(2-2.5、2.6)计算各组分的活度因子γk ,I ; (6) 由式(2-2.9)计算各组分新的气相组成y 1；

k k L k m k k k k k p RT p p V x p y ?γ?/]/)(exp[*

*,,**-=I ，K k ,,1 =

(2-2.9)

(7) 比较y 1和y 0，如果两者不相等，则令y 0=y 1，转步骤(2)，进行新一轮

循环迭代，直至达到规定的计算进度。

参考文献

[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)