基于COSMO-RS的單工質(zhì)氣液相平衡熱物性計算方法

(上海理工大學(xué)能源與動力工程學(xué)院 上海200093)

由于環(huán)境因素的限制，制冷系統(tǒng)中對制冷劑的要求越來越高[1]，越來越多的專家學(xué)者投入到新型制冷劑的研發(fā)中。氣液相平衡數(shù)據(jù)對于制冷劑的工程應(yīng)用和設(shè)計意義重大，是熱力學(xué)中不可或缺的組成部分。氣液相平衡數(shù)據(jù)主要可以通過狀態(tài)方程法、活度系數(shù)法計算及模擬仿真和實驗測定。吳獻忠等[2]利用UNIFAC基團貢獻法預(yù)測了混合制冷劑R32+R227ea、R227ea+R134a和R152a+R227ea及三元混合制冷劑R32+R125+R134a的氣液相平衡數(shù)據(jù)。宋有強等[3]采用PR狀態(tài)方程及其混合法則計算了混合制冷劑CO2+R600的氣液相平衡數(shù)據(jù)。董學(xué)強等[4]利用基于循環(huán)法的高精度可視化中低溫相平衡裝置測量了溫度為255、265、275、288 K時R134a/R290的氣液相平衡數(shù)據(jù)。狀態(tài)方程法計算嚴(yán)格，但參數(shù)多，應(yīng)用范圍具有局限性；活度系數(shù)法雖然計算簡單，但需要大量的實驗數(shù)據(jù)，對于計算新型制冷劑時還有困難；傳統(tǒng)意義上可以通過實驗手段獲得氣液相平衡數(shù)據(jù)，但實驗結(jié)構(gòu)設(shè)計復(fù)雜、費用大，對于易燃易爆、帶有毒性的制冷劑還具有一定的危險性[5]。

隨著計算化學(xué)理論和計算機技術(shù)的發(fā)展，模擬仿真逐漸變成計算氣液相平衡的主要工具之一[6]。趙勝喜等[7]利用吉布斯蒙特卡羅法模擬R410A的氣液相平衡并得到臨界狀態(tài)點。陳秀萍等[8]對比了分子動力學(xué)、吉布斯系綜蒙特卡羅和真實溶劑似導(dǎo)體屏蔽模型3種方法的熱物性模擬結(jié)果，前兩種計算方法由于目前分子間相互作用力理論不夠完善，缺少專用于制冷劑的力場，勢能函數(shù)的精度不夠高，例如，力場中的勢能參數(shù)具有局限性，不能描述所有物質(zhì)，且大部分力場不能很好的處理分子間的靜電作用，因此還存在一定誤差。真實溶劑似導(dǎo)體屏蔽模型(COSMO-RS)[9]通過物質(zhì)的分子結(jié)構(gòu)信息來計算物質(zhì)的氣液相平衡，避開了勢能函數(shù)精度的影響，只需要通過調(diào)節(jié)分子的尺寸因子。該方法允許快速并自洽的處理溶液中的分子，且在量子化學(xué)水平上最大程度的近似溶劑的特點：通過經(jīng)典的介質(zhì)理論描述溶劑與溶質(zhì)的靜電相互作用。陳秀萍等[10-11]采用COSMO-RS模擬了R290/R227ea和R290/R1234ze的氣液相平衡數(shù)據(jù)，誤差≤4%。Bai Zhenmin等[12]利用真實溶劑似導(dǎo)體屏蔽模型模擬了甲醛與水的氣液相平衡，并與實驗值對比，具有很好的一致性。但對于COSMO-RS在計算眾多種類單工質(zhì)制冷劑的氣液相平衡時準(zhǔn)確度是否會出現(xiàn)變化缺乏研究。因此本文針對無機物、碳?xì)浠衔锖蜌浞鸁N3種不同化學(xué)類型的單工質(zhì)制冷劑進行氣液相平衡計算，并對其精度進行研究。

1 COSMO-RS理論

1.1 COSMO-RS理論

COSMO-RS是一種連續(xù)介質(zhì)溶劑化模式，模型中將連續(xù)介質(zhì)的介電常數(shù)設(shè)為無窮大(理想導(dǎo)體)，這樣將屏蔽電荷限制在界面上，從而使分子和溶劑間沒有電場，導(dǎo)體內(nèi)也沒有電荷[13-14]。COSMO-RS只需要物質(zhì)的結(jié)構(gòu)信息就可以通過對分子進行量子化學(xué)計算來預(yù)測多元體系的相平衡及其他熱物性，是目前量子化學(xué)和工程熱力學(xué)之間最好的連接手段。作為一種描述流體和熱力學(xué)性質(zhì)的方法，由于不需要實驗數(shù)據(jù)，COSMO-RS發(fā)展迅速。

COSMO-RS建立在分子間片段對相互作用的概念上，即在COSMO-RS模型中將分子分解成面積大小相等的片段，分子間相互作用的概念建立在表面片段相互作用的物理觀點上。用相互接觸的片段對的凈屏蔽電荷密度σ和σ′來計量兩片段在真實體系和理想導(dǎo)體中的能量差別Emisfit，如式(1)所示。

(1)

式中：aeff為有效接觸表面積， nm2；α′為能量因子，可由靜電理論計算。若研究體系中還存在強極性物質(zhì)，還需考慮氫鍵的相互作用EHB，如式(2)所示。

(2)

A=min(0,σdon+σhb)

(3)

B=max(0,σacc-σhb)

(4)

式中：σdon=min(σ,σ′)，σacc=max(σ,σ′)；σhb為氫鍵閾值；chb為強度系數(shù)。σhb、chb必須與實驗值吻合。

式(2)、式(3)化簡后主要有兩種可能性：一是當(dāng)兩個相接觸的表面屏蔽電荷中負(fù)極部分的電荷密度σ<-σhb；二是正極部分的電荷密度σ′>σhb，氫鍵為零。此時，正電荷是接受體，負(fù)電荷是供體。氫鍵正比于剩余部分電荷(σdon+σhb)(σacc-σhb)。

除了考慮Emisfit和Ehb外，還需要考慮分子片段間的范德華作用力Evdw，如式(5)所示。

(5)

分子片段間相互作用的總能量可由式(6)得到：

E(σ,σ′)=Emisfit+EHB+Evdw

(6)

1.2 σ-profile

分子表面的屏蔽電荷密度σ是采用COSMO-RS理論最重要的數(shù)據(jù)。每個片段的表面屏蔽電荷密度能用一個合適的平均值表示，能滿足相互配對。可以采用量子化學(xué)方法得到這些片段的表面電荷密度，并經(jīng)過一個基于式(7)在半徑為rav的面積上的平均。

(7)

式中：dij為片段i和片段j之間的距離，nm；ri為根據(jù)片段面積計算的片段i的平均半徑，nm；rav為可調(diào)參數(shù)。

由分子表面所有片段的表面屏蔽電荷密度構(gòu)成的整個分子的表面屏蔽電荷密度分布直方圖，稱為分子的表面屏蔽電荷密度分布函數(shù)P(σ)，亦稱為分子的“σ-profile”，它表示在分子表面上找到具有平均屏蔽電荷密度為σ的面積的數(shù)量。

對于混合物S系綜的表面屏蔽電荷密度分布函數(shù)Ps(σ)，是體系中各種組分Xi的PXi(σ)按其摩爾分?jǐn)?shù)xi的加權(quán)平均，如式(8)所示。

PS(σ)=∑i∈SxiPXi(σ)

(8)

為方便進行統(tǒng)計熱力學(xué)計算，對Ps(σ)進行歸一化處理。由于每個分子的總表面積可以利用該分子的PXi(σ)在整個σ變化區(qū)間上積分求得，因而系綜的標(biāo)準(zhǔn)屏蔽電荷密度分布函數(shù)可定義為：

(9)

式中：AXi為組分Xi的分子的總表面積，nm2，AS為系綜S的分子的平均總表面積，nm2。

1.3 單工質(zhì)飽和蒸氣壓計算

假定一個沒有自由表面的大量液體系統(tǒng)，即分子的表面片段之間接觸，該假設(shè)僅適用于遠(yuǎn)離臨界狀態(tài)的液體。在該假設(shè)條件下，系綜的統(tǒng)計熱力學(xué)可以按照式(10)計算:

(10)

(11)

式中：μS(σ)為系綜S在溫度T下，接觸面積為aeff的片段化學(xué)勢，kJ。由于μS(σ)存在于方程的兩邊，因此該方程采用迭代求解，μS(σ)初始值為零。

該片段的活度系數(shù)可由式(12)獲得：

(12)

由式(10)和式(12)直接定義溶質(zhì)X在系綜S中的化學(xué)勢：

(13)

(14)

(15)

式中：AS為溶質(zhì)X分子的平均表面積，nm2；λ為調(diào)節(jié)參數(shù)。

在COSMO-RS理論中，除了預(yù)測液相工質(zhì)的熱力學(xué)性質(zhì)，COSMO-RS理論還提出了一個方法，用于預(yù)測純組分氣相化學(xué)勢。

(16)

根據(jù)熱力學(xué)理論，單一組分的飽和蒸氣壓可根據(jù)式(17)計算：

(17)

在COSMO-RS中，尺寸因子s是可調(diào)的，其缺省值為1.0，實際可根據(jù)實驗值擬合。本文以氣液相平衡數(shù)據(jù)為目標(biāo)值，對各類制冷劑工質(zhì)的尺寸因子進行擬合，獲得了準(zhǔn)確的氣液相平衡數(shù)據(jù)。調(diào)節(jié)s將對溶質(zhì)的液相化學(xué)勢產(chǎn)生影響，從而影響式(17)，其中的理論依據(jù)還有待進一步探討。

2 基于COSMO-RS的單工質(zhì)的氣液相平衡模擬

制冷劑按化學(xué)組成分類，主要有無機物、碳?xì)浠衔?HCs)和鹵化物(氟利昂)3類。由于具有相似化學(xué)組成的化合物在COSMO中生成屏蔽電荷的分布具有相似性，對尺寸因子s的規(guī)律尋找也具有一定的指引。基于臭氧層的保護，氟利昂制冷劑中CFCs已被淘汰使用，HCFCs也已被限定使用期限。氫氟烴類簡稱HFCs，臭氧層的破壞系數(shù)為0，但氣候變暖潛能值(GWP)很高，被視為CFCs和HCFCs類物質(zhì)的重要過渡性替代物質(zhì)。目前HFCs的大氣濃度較低，對于氣候變化的作用很小[15]。制冷系統(tǒng)的工作溫度范圍(蒸發(fā)溫度～冷凝溫度)通常為-20～60 ℃，因此本文模擬的溫度范圍為250～330 K，部分工質(zhì)的模擬溫度根據(jù)其臨界點做出調(diào)整。外延到其他溫區(qū)范圍時，根據(jù)外延溫區(qū)范圍的增大，誤差變大，約為±10%，可以通過調(diào)節(jié)s值來減小該誤差。模擬給定的迭代起點溫度大部分設(shè)定為273.15 K，部分迭代起點溫度為模擬溫度上限和下限的中間值。在計算中，s以與制冷劑特性查詢軟件NIST的最大誤差最小來調(diào)解。

s與分子化學(xué)組成和工作溫度范圍密切相關(guān)。由于制冷劑的工作環(huán)境有一定的限制，因此將工作溫度限制在一定溫度范圍內(nèi)(250～330 K)。限制了溫度對s的影響。同時對于制冷劑化學(xué)組成的分類也降低了s變化規(guī)律的不確定性。

2.1 無機物制冷劑的氣液相平衡模擬

圖1 無機物制冷劑的氣液相平衡曲線模擬值與NIST值對比Fig.1 The comparison curve of inorganic refrigerant between simulated value and NIST value of gas-liquid equilibrium

無機物是制冷劑發(fā)展的第一個階段[16]，無機物制冷劑不破壞臭氧層，產(chǎn)生溫室效應(yīng)極小，且獲取容易、價格低廉。其中CO2和NH3是最常見的兩種無機物制冷劑。利用COSMO-RS理論對這兩種無機物制冷劑進行氣液相平衡數(shù)據(jù)的模擬，模擬結(jié)果如圖1所示。圖1中，左上角曲線圖為氣液相平衡壓力模擬值與NIST值的相對誤差，右上角為該分子的表面電荷密度圖。模擬中發(fā)現(xiàn)，當(dāng)s=1.0時，計算值誤差較大。因此，逐漸改變s，發(fā)現(xiàn)通過調(diào)節(jié)分子的s，可使模擬的氣液相平衡曲線逐漸接近NIST的實驗值。圖1(a)中左上角的圖顯示了s對相對誤差的影響。圖中模擬結(jié)果與NIST值之間的相對誤差通過式(18)計算：

(18)

式中：pS為模擬壓力，MPa；pN為對應(yīng)溫度和壓力下NIST中的壓力， MPa。

由圖1(a)可知，模擬的氣液相平衡曲線與NIST提供的氣液相平衡曲線趨勢完全一致。當(dāng)CO2的s=4.4時，在低溫區(qū)和高溫區(qū)相對誤差較大，但仍在±2.5%以內(nèi)。圖1(b)中當(dāng)NH3的s=0.37時，最大相對誤差δmax在±1.0%以內(nèi)。

2.2 碳?xì)浠衔镏评鋭┑臍庖合嗥胶饽M

碳?xì)浠衔锏闹评鋭┌ㄍ闊N類和烯烴類，屬于天然工質(zhì)，因此對大氣無污染、對臭氧層勿破壞和溫室效應(yīng)幾乎為零。但碳?xì)浠衔镒鳛橹评鋭r，可燃性、安全性需要特別注意。利用COSMO-RS理論對部分碳?xì)浠衔镞M行氣液相平衡模擬，烷類碳?xì)浠衔锏挠嬎憬Y(jié)果、模擬結(jié)果與NIST的相對誤差如圖2所示。

除了甲烷的模擬溫度在130～190 K,其余3種烷類碳?xì)浠衔锏哪M溫度范圍均在250～330 K之間。CH4的s=1時，δmax在±2%以內(nèi)；乙烷、丙烷和丁烷的s分別為0.004、0.002和0.0001時，乙烷和丙烷的δmax在±1%以內(nèi)，丁烷的δmax在±3%以內(nèi)，與NIST值均具有很好的一致性。圖2(e)和圖2(f)分別為乙烯和丙烯的氣液相平衡曲線計算結(jié)果。由圖2(e)可知，當(dāng)乙烯的s=1.25時，δmax在±3%以內(nèi)；由圖2(f)可知，丙烯的s=4.1時，δmax在±1.5%以內(nèi)。

2.3 氫氟烴類制冷劑的氣液相平衡模擬

相對于碳?xì)浠衔?，氫氟烴類制冷劑的化學(xué)結(jié)構(gòu)不對稱。主要針對14種不同的氫氟烴類制冷劑進行氣液相平衡的模擬。模擬結(jié)果及與NIST的相對誤差如圖3所示。

碳?xì)浠衔镏械臍湓颖环犹娲螅叽缫蜃影l(fā)生變化。選擇合適的尺寸因子后，12種氫氟烷(HFCs)和2種氫氟烯(HFOs)類工質(zhì)在模擬溫度范圍內(nèi)，相平衡曲線模擬值與NIST值均具有很好的一致性，δ<±2%，僅C4F10的δ<±6%。

2.4 尺寸因子變化規(guī)律分析

表1所示為22種單工質(zhì)的尺寸因子。

圖4所示為烷烴化合物的s隨碳原子的變化，s隨著碳原子的增加而減小，縱坐標(biāo)取對數(shù)坐標(biāo)時，近似呈線性變化。

圖5所示為全氟碳化合物(CnF2n+2)的s隨碳原子的變化，s也隨碳原子數(shù)的增加而減小，縱坐標(biāo)取對數(shù)坐標(biāo)時，呈較好的線性分布。

甲烷的氟化物(CHnF(4-n))、乙烷的氟化物(C2HnF(6-n))、丙烷的氟化物(C3HnF(8-n))的s隨氟原子數(shù)的變化如圖6所示。由圖6可知，s隨氟原子數(shù)的增大而增大，但全氟替代物CF4、C2F6、C3F8的s突然大幅減小，甚至小于對應(yīng)的烷類。這是由于氟原子間的相互吸引力較大所致。圖6中的變化趨勢完全一致，即先緩慢增大，最后突然增大，又急劇減小。

甲烷、乙烷和丙烷的氟化物都可以根據(jù)圖表中的規(guī)律選取相應(yīng)的s范圍。對于烷類化合物及烷類氟化物，在模擬溫度范圍內(nèi)可以預(yù)測出s的數(shù)值范圍，在數(shù)值范圍內(nèi)，預(yù)測出的模擬數(shù)據(jù)對實驗具有一定的指導(dǎo)作用。因此，可以根據(jù)s的變化規(guī)律，在已知化合物分子結(jié)構(gòu)的情況下，根據(jù)一個相平衡點和s，對沒有實驗數(shù)據(jù)的未知化合物的氣液相平衡數(shù)據(jù)進行預(yù)測，其準(zhǔn)確度約為±2%。

2.5 飽和蒸氣壓方程擬合

飽和蒸氣壓的計算在工程中非常重要，可以由狀態(tài)方程導(dǎo)出，但較為復(fù)雜，故工程上通常采用直接回歸的飽和蒸氣壓方程。飽和蒸氣壓的形式很多，常用的有Antoine方程和Wagner方程。三參數(shù)的Antoine方程是工程上最常用的飽和蒸氣壓方程，特點是簡單、系數(shù)少、較為精確[17]。其形式為：

(19)

圖2 碳?xì)浠衔镏评鋭┑臍庖合嗥胶馇€模擬值與NIST值對比Fig.2 The comparison curve of HCs refrigerant between simulated value and NIST value of gas-liquid equilibrium

分類工質(zhì)尺寸因子s工質(zhì)尺寸因子s無機物CO2(R744)4.4NH3(R717)0.37烷烴CH4(R50)1C2H6(R170)0.02C3H8(R290)0.002C4H10(R600)0.000 1烯烴C2H4(R1150)1.25C3H6(R1270)4.1HFCsCH2F2(R32)1.8CHF3 (R23)4.9CF4 (R14)0.001 2C2H4F2 (R152a)1.5C2H3F3 (R143a)2.5C2H2F4 (R134a)3.2C2HF5 (R125)8.5C2F6 (R116)0.000 3C3H3F5 (R245fa)2.7C3HF7 (R227ea)15C3F8 (R218)0.000 06C4F100.000 02HFOsC3H2F4(R1234ze)3.7C3H2F4(R1234yf)6.1

圖3 氫氟烴類制冷劑的氣液相平衡曲線模擬值與NIST值對比Fig.3 The comparison curve of HFCs refrigerant between simulated value and NIST value of gas-liquid equilibrium

圖4 烷烴類(CnH(2n+2))工質(zhì)s值碳原子數(shù)的變化Fig.4 The variation trend of scale factor of CnH(2n+2)with the number of carbon atoms

圖5 全氟化烷類(CnF(2n+2))工質(zhì)s隨碳原子數(shù)的變化Fig.5 The variation trend of scale factor of CnF(2n+2) with the number of carbon atoms

圖6 甲烷氟化物、乙烷氟化物、丙烷氟化物工質(zhì)s隨氟原子數(shù)的變化Fig.6 The variation trend of scale factor of CHnF(4-n)，C2HnF(6-n) and C3HnF(8-n) with the number of carbon atoms

式中：p為壓力，MPa；T為溫度，℃，a1、a2、a3均為常數(shù)由模擬數(shù)據(jù)擬合得出。表2所示為Antoine方程中各工質(zhì)的三參數(shù)及擬合出的Antoine方程中飽和壓力與模擬出的飽和壓力之間的平均誤差，平均偏差定義為：

(20)

式中：pA為Antoine狀態(tài)方程計算的壓力值，MPa；pS為模擬出的飽和壓力，MPa。

由表2可知，擬合曲線與模擬數(shù)據(jù)之間的平均偏差很小，擬合出的Antoine方程在模擬的溫度范圍內(nèi)與模擬數(shù)據(jù)具有較高的一致性。

表2 Antoine方程中各工質(zhì)的參數(shù)擬合及擬合方程與模擬值偏差Tab.2 The parameters in Antoine equation and the deviation between the value of fitting equation and simulated value

3 結(jié)論

本文采用COSMO-RS模型研究了純制冷劑工質(zhì)的氣液相平衡的熱物性,總結(jié)了不同種類單工質(zhì)尺寸因子的變化規(guī)律。對比模擬結(jié)果與NIST提供的數(shù)據(jù)，無機物CO2和NH3的相對誤差分別在±3%和±1%以內(nèi)；烷類和烯類(HCs)的相對誤差均在±3%以內(nèi)；氫氟烴類(HFCs)的相對誤差均在±2%以內(nèi)，除了C4F10的相對誤差在±6%以內(nèi)。因此，COSMO-RS氣液相平衡模擬結(jié)果與NIST提供的數(shù)據(jù)具有很好的一致性。得到尺寸因子的變化規(guī)律如下：

1)烷烴類工質(zhì)CnH2n+2的尺寸因子隨碳原子數(shù)的增加而減小。

2)全氟化烷類工質(zhì)CnF2n+2的尺寸因子隨碳原子數(shù)的增加而減小。

3)氫氟烷類工質(zhì)CnHF2n+1(HFCs)的尺寸因子隨氟原子數(shù)的增加而增大；但當(dāng)氫原子全部被氟原子替代后，尺寸因子又大幅減小，甚至小于相對應(yīng)的烷烴類。

4)氫氟烯類工質(zhì)CnHF2n-1(HFOs)的尺寸因子隨氟原子數(shù)的增加而增大。

5)根據(jù)模擬數(shù)據(jù)擬合出Antoine方程的系數(shù)。

根據(jù)以上尺寸因子的變化規(guī)律，結(jié)合分子結(jié)構(gòu)，以及至少一點的氣液相平衡數(shù)據(jù)，可以預(yù)測新型單工質(zhì)制冷劑的氣液相平衡曲線。