运用Excel计算纯流体饱和热力学性质

陈漓 莫小梅 孙美娟

(百色学院 a.材料科学与工程学院； b.化学与环境工程学院， 广西 百色 533000)



运用Excel计算纯流体饱和热力学性质

陈漓a莫小梅b孙美娟a

(百色学院 a.材料科学与工程学院； b.化学与环境工程学院， 广西 百色 533000)

利用Excel的数组公式和相关数据处理功能，通过简单的设计，解决立方型状态方程计算流体饱和热力学性质所遇到的非线性高次方程、多次迭代问题，有效避开手工计算难度大和计算机计算需编程问题等，而且计算过程直观、可视，利于掌握，计算结果准确、客观。该方法可为研究者、学习者和教学者提供参考。

Excel；饱和蒸汽压；改进PR状态方程；热力学性质

0 引言

计算纯流体饱和热力学性质，常常需要状态方程，由于立方型状态方程参数少且有较好的计算精度，因此得到广泛的运用。在众多立方形方程中PR状态方程应用最为普遍。PR方程在预测液体体积上较之之前方程有所改善，但要进一步外推就比较困难。DAKEWU和SHULINCHEN修正PR方程，使方程的计算精度有显著改进，在计算纯流体的饱和蒸气压、饱和液相密度方面有良好的准确度。在此我们以改进PR状态方程为模型结合Excel电子表格计算纯流体的饱和热力学性质。

纯流体处于饱和状态时，一般用4个状态参量描述，温度T、压强p、汽相摩尔体积Vν和液相摩尔体积Vl。然而系统平衡时只有一个独立变量，通常取温度T或压强p，故有两类计算过程:第一类是蒸汽压计算；第二类是沸点计算。

以第一类的蒸气压计算为例说明，并以改进PR状态方程为模型[1]

(1)

纯流体在饱和状态时需满足汽液平衡准则式:φν=φl

(2)

其中φν为汽相逸度，φl为液相逸度。

结合状态方程，其逸度系数表达式可写为：

(3)

纯流体的汽液平衡的计算，首先要根据方程输入临界参数和偏心因子，计算给定温度T下的改进PR方程常数a，b。假设p的初值，求状态方程得到汽相摩尔体积Vν和液相摩尔体积Vl，由此判别方程式(2)是否满足收敛条件，若不满足，通过调节p，直到方程式(2)收敛，此时的p、Vν和Vl就是方程组式(1)和式(2)的解。由于在计算中涉及到高次方程运算，需借助计算机进行编程计算，操作过程较为繁琐。如果能用办公软件Excel来解决，将给计算带来很大的便利。

Excel是微软办公软件的一个重要部分，它可以进行各种数据处理、统计分析等操作。对于热力学数据的处理和运算也具有良好的交互界面，可简单明了地展示计算原理与过程，有利于解决繁琐的热力学系统计算问题[2-5]。

1 纯流体汽液平衡的理论基础

1.1 计算饱和蒸汽压

饱和蒸汽压p进行首轮迭代计算时，要预先估计蒸汽压的初值。运用Ednister偏心因子估算式ω=[3Tbr/7(1-Tbr)]lgpc-1能得到临界参数和偏心因子估计蒸汽压的初值的计算式[6]：

(4)

对于饱和蒸汽压p的迭代式，可从式(2)利用Newton迭代法得到

p(k+1)=p(k)-[RT(lnφν-lnφ1)/(Vν-Vl)](k)

(5)

1.2 计算汽相、液相摩尔体积

对于立方形状态方程摩尔体积的计算，当T

Vk+1=Vk-f(Vk)/f′(Vk)

(6)

为了方便迭代，将改进PR方程用摩尔体积的三次展开式来表示，即

(7)

f ′(V)为f(V)的一阶导。

汽相摩尔体积Vν的求取通常以理想气体体积V=RT/p为初值，而液相摩尔体积以V=b为初值，代入式(6)，得到Vl值后再代入等式右边，一直迭代下去，直到满足Vl<ε。利用Excel的计算功能，可方便同时迭代出Vν和Vl，这样就极大提高计算的效率。

1.3 计算偏离焓、偏离熵等热力学性质

流体的各种热力学性质，比如焓、熵值和热容等，可以利用偏离函数并结合状态方程计算出来[7]。

(8)

(9)

(10)

偏离定压摩尔热容：

(11)

2 在Excel中进行纯流体饱和热力学性质的计算

为计算纯流体汽液平衡下的饱和热力学性质，在Excel中进行如下操作：

(1)以乙烷为例，分别在B3到E3单元格输入临界温度、临界压强、气体常数和偏心因子等数值。A6到F6单元格为方程常数a，b的值以及相应一阶导数、二阶导数的数值。

(2)在B10单元格输入某一温度“220”如图1所示，即假设此时沸点温度为220 K，运用式(4)计算饱和蒸气压的初值，即在A12单元格输入“=C3*10^(7*(1+E3)/3*(1-B3/B10))”，把计算结果以粘贴数值的方式粘贴到A10单元格里。

(3)在单元格C9和D9分别输入Vν的初值(V=RT/p)“=D3*B10/(A10*10^6)”和Vl的初值(V=b)“=B6”，单元格C10-C16和D10-D16分别输入Newton迭代式(6)，如C10-C16单元格输入如下数组公式“{=C9:C15-(C9:C15^3-(D3*B10/(A10*10^6)-0.645*B6)*C9:C15^2+1/(A10*10^6)*(D6-1.645*B6*D3*B10-2.29*A10*10^6*B6^2)*C9:C15-D6*B6/(A10*10^6)+0.645*D3*B10*B6^2/(A10*10^6)+0.645*B6^3)/(3*C9:C15^2-2*(D3*B10/(A10*10^6)-0.645*B6)*C9:C15+1/(A10*10^6)*(D6-1.645*B6*D3*B10-2.29*A10*10^6*B6^2))}”。同样D10-D16单元格中式子与C10-C16内容相类似，所不同之处仅仅是把式中“C9:C15”改为“D9:D15”即可。一般来说大多数情况下迭代5-6次即可得到满意结果,如图1所示。为确保迭代的精度我们进行了7次迭代，单元格C16和D16得到的数值分别为汽相摩尔体积Vν和液相摩尔体积Vl，参见图1。单元格F9和G9分别为式(7)f(V)的Vν和Vl收敛情况，f(V)应趋于零。

(4)接下来进一步计算逸度系数。在单元格F11和G11分别输入计算汽相和液相逸度的关系式(3)，在单元格F11输入“=A10*10^6*F10/(D3*B10)-1-LN(A10*10^6*(F10-B6)/(D3*B10))-D6/(2.299135707*B6*D3*B10)*LN((F10+1.9720678535*B6)/(F10-0.3270678535*B6))”。同理在单元格G11输入液相逸度的关系式，计算后如果不满足|lnφν-lnφl|≤ε则需要调整蒸汽压p进行迭代计算，迭代关系式参见式(5)，即在单元格A14输入式(5)，即“=(A10*10^6*(1-(F11-G11)/(A10*10^6*F10/(D3*B10)-A10*10^6*G10/(D3*B10))))/10^6”。把A14计算结果以粘贴数值的方式到A10单元格里，重复粘贴数值，Excel将对整个过程进行自动迭代计算，直至满足|lnφν-lnφl|≤ε，一般迭代3-4次即能得到满意的结果，此时得到p、Vν和Vl分别为饱和蒸汽压、汽相摩尔体积和液相摩尔体积的数值。

(5)计算其它的热力学性质。在单元格F13到F16和G13到G16分别输入式(8)、式(9)、式(10)和式(11)，这样我们同时得到乙烷汽相和液相偏离焓、偏离熵、偏离定容摩尔热容和偏离定压摩尔热容等数值。例如，在单元格F13输入“=A10*10^6*F10-D3*B10-1/(2.299135707*B6)*(D6-B10*E6)*LN((F10+1.9720678535*B6)/(F10-0.3270678535*B6))”，计算得到的数值为汽相偏离焓。

3 应用分析

为了检验Excel在计算纯流体饱和性质的可靠性，我们以乙烷为研究对象，对该气体饱和蒸汽压等性质进行计算。在单元格B4、B5分别输入乙烷的临界温度Tc和临界压强pc，单元格B7输入偏心因子ω的数值参见图1，计算乙烷温度为160～305 K的饱和蒸汽压、摩尔体积和偏离性质。计算的结果参见表1。

图1 乙烷饱和性质的计算

温度T/K饱和汽压p/MPa摩尔体积汽相Vν/m3·mol-1液相Vl/m3·mol-1偏离焓汽相[HR]ν/J·mol-1液相[HR]l/J·mol-1偏离熵汽相[SR]ν/J·(mol·K)-1液相[SR]l/J·(mol·K)-1160002126225×10-25022×10-5-289508-156550863-01076-977709180007781875×10-25238×10-5-925438-150629562-03116-834806200021567300×10-35511×10-5-2285760-144285276-07085-717083220049043363×10-35867×10-5-4753996-137265643-13765-616090240096791728×10-36356×10-5-8817424-129190432-24176-525730260172119451×10-47081×10-5-15216885-119398562-40066-440765280282805238×10-48317×10-5-25472798-106428417-65514-354641300437522585×10-41153×10-4-45660183-84250988-118712-247348305484071903×10-41503×10-4-57951513-69336618-154363-191691

从表1计算的数据可以看出随着温度的增大，蒸发加剧密度增大饱和蒸汽压升高，当温度升高到305 K时，饱和汽压逐步接近临界压强。同时温度升高汽相摩尔体积Vν将随着饱和汽压升高而减小，液相摩尔体积Vl则随着温度升高而增大，在临界点附近它们数值趋于一致。同样偏离理想气体性质焓和熵，汽相随饱和汽压升高，越远离理想气体状态，液相随温度升高，偏离程度减小，在临界点附近它们数值趋于一致。

把表1计算的饱和汽压、摩尔体积、偏离焓和偏离熵结果与国内外计算热力学性质软件得到的数值是一致的，这也说明该方程能够较好的描述纯流体的热力学性质。

上面的饱和热力学性质的计算是以温度T为独立变量进行的，同样也可以以蒸汽压p为独立变量来进行饱和热力学性质的计算，只需在单元格B14输入如下的迭代式即可

在单元格B14输入“=B10*(1+(F11-G11)/((F13-G13)/(D3*B10)))”，多次迭代直到满足|lnφν-lnφl|≤ε，此时温度T为饱和蒸汽压沸点的数值。

4 结论

通过上面的例子可以看出，运用Excel进行纯流体汽液两相饱和性质的计算，无需编程，只需在表格里输入纯流体的临界性质和偏心因子的数值，通过简单的迭代设计，Excel即可计算各种纯流体的饱和热力学性质。整个运算过程有良好的交互界面，并能有效地提高计算结果的准确度和可靠性。

Exce1具有强大的运算功能，但以往多用于简单计算。实际上，充分利用Exce1提供的运算功能，只需输入相应的关系式，就能完成复杂的热力学工程计算，计算过程简单、直观，学生非常容易掌握。由于Exce1具有较强的数组公式和公式的复制功能，所以使用Exce1计算纯流体热力学性质具有无可比拟的优越性。

[1]DAKEWU,SHULINCHEN.AmodifiedPeng-Robinsonequationonstate[J].ChemicalEngineeringCommunications， 1997, 156(1)：215-225.

[2]陈方.用EXCEL2000解决方程求根问题[J].湖南工业职业技术学院学报，2003，4(1)：14-16.

[3]班玉凤，常圣泉，朱海峰，等.EXCEL在非理想系泡露点计算中的应用[J].计算机应用与软件，2011，28(10)：275-277.

[4]于志家，陈传棋，李香琴，等.应用Excel进行泡点与露点计算[J].化工高等教育，2012，126(4)：73-76.

[5]王智娟，胡粉娥，杨晓丽，等.Excel在化工热力学计算中的应用[J].广州化工，2013，41(3)：167-169.

[6]陈钟秀，顾飞燕，胡望明，等.化工热力学[M].3版.北京：化学工业出版社，2012：287-288.

[7]陈则韶.高等工程热力学[M].2版.合肥：中国科学技术大学出版社，2014：120-125.

[责任编辑 刘景平]

Using Excel to Calculate the Thermodynamic Properties of Pure Fluid Saturation

CHEN Lia， MO Xiaomeib， SUN MeiJuana

(a.School of Materials Science and Engineering;b.School of Chemistry and Environmental Engineering, Baise University, Baise, Guangxi 533000, China)

The study uses Excel array formulas and data processing functions, through a simple design, to solve the problems of nonlinear equation of higher degree and multiple iterations faced when calculating thermodynamic properties of fluid saturation. It can effectively avoid the problems of manual calculation being difficult and computer calculation needing programming; the calculation process is intuitive, visual, conducive to grasp, and the calculation results are accurate and objective. The method provides a reference for researchers, learners and educators.

Excel; saturated vapor pressure; modified PR equation of state; thermodynamic properties

O642

A

1672-9021(2016)05-0113-05

陈漓(1962-)，男，广西桂林人，百色学院材料科学与工程学院副教授，主要研究方向：热力学。

广西高校科学技术研究基金资助项目(2013YB244)。

2016-09-12