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

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

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

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

1 COSMO-RS理論

1.1 COSMO-RS理論

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

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

(1)

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

(2)

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

(3)

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

(4)

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

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

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

(5)

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

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

(6)

1.2 σ-profile

分子表面的屏蔽電荷密度σ是采用COSMO-RS理論最重要的數據。每個片段的表面屏蔽電荷密度能用一個合適的平均值表示，能滿足相互配對?？梢圆捎昧孔踊瘜W方法得到這些片段的表面電荷密度，并經過一個基于式(7)在半徑為rav的面積上的平均。

(7)

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

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

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

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

(8)

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

(9)

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

1.3 單工質飽和蒸氣壓計算

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

(10)

(11)

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

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

(12)

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

(13)

(14)

(15)

式中：AS為溶質X分子的平均表面積，nm2；λ為調節參數。

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

(16)

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

(17)

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

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

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

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

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

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

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

(18)

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

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

2.2 碳氫化合物制冷劑的氣液相平衡模擬

碳氫化合物的制冷劑包括烷烴類和烯烴類，屬于天然工質，因此對大氣無污染、對臭氧層勿破壞和溫室效應幾乎為零。但碳氫化合物作為制冷劑時，可燃性、安全性需要特別注意。利用COSMO-RS理論對部分碳氫化合物進行氣液相平衡模擬，烷類碳氫化合物的計算結果、模擬結果與NIST的相對誤差如圖2所示。

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

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

相對于碳氫化合物，氫氟烴類制冷劑的化學結構不對稱。主要針對14種不同的氫氟烴類制冷劑進行氣液相平衡的模擬。模擬結果及與NIST的相對誤差如圖3所示。

碳氫化合物中的氫原子被氟原子替代后，尺寸因子發生變化。選擇合適的尺寸因子后，12種氫氟烷(HFCs)和2種氫氟烯(HFOs)類工質在模擬溫度范圍內，相平衡曲線模擬值與NIST值均具有很好的一致性，δ<±2%，僅C4F10的δ<±6%。

2.4 尺寸因子變化規律分析

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

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

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

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

甲烷、乙烷和丙烷的氟化物都可以根據圖表中的規律選取相應的s范圍。對于烷類化合物及烷類氟化物，在模擬溫度范圍內可以預測出s的數值范圍，在數值范圍內，預測出的模擬數據對實驗具有一定的指導作用。因此，可以根據s的變化規律，在已知化合物分子結構的情況下，根據一個相平衡點和s，對沒有實驗數據的未知化合物的氣液相平衡數據進行預測，其準確度約為±2%。

2.5 飽和蒸氣壓方程擬合

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

(19)

圖2 碳氫化合物制冷劑的氣液相平衡曲線模擬值與NIST值對比Fig.2 The comparison curve of HCs refrigerant between simulated value and NIST value of gas-liquid equilibrium

分類工質尺寸因子s工質尺寸因子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))工質s值碳原子數的變化Fig.4 The variation trend of scale factor of CnH(2n+2)with the number of carbon atoms

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

圖6 甲烷氟化物、乙烷氟化物、丙烷氟化物工質s隨氟原子數的變化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均為常數由模擬數據擬合得出。表2所示為Antoine方程中各工質的三參數及擬合出的Antoine方程中飽和壓力與模擬出的飽和壓力之間的平均誤差，平均偏差定義為：

(20)

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

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

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

3 結論

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

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

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

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

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

5)根據模擬數據擬合出Antoine方程的系數。

根據以上尺寸因子的變化規律，結合分子結構，以及至少一點的氣液相平衡數據，可以預測新型單工質制冷劑的氣液相平衡曲線。