CO2混合工質的氣液相平衡的混合規則對比與預測研究

吳子睿，孫瑞，石凌峰，田華，王軒，舒歌群

（1 中國科學技術大學熱科學和能源工程系，安徽合肥 230027； 2 天津大學內燃機燃燒學國家重點實驗室，天津 300072）

引 言

二氧化碳(CO2)屬于自然工質，由于其高熱穩定性、安全環保并且高能量密度實現部件小型化等優點，使得CO2動力循環在新一代熱功轉化循環、內燃機余熱、中低溫熱能、火電、核電和太陽能利用等領域具有廣泛的應用前景。同時CO2作為循環介質，也是一種對CO2的直接資源化利用方式，符合國家雙碳政策，但是CO2純工質循環運行壓力高，應用中存在安全隱患，且臨界溫度較低對冷源要求很高[1]。碳氫類(HCs)工質雖能提高能源利用效率，但HCs易燃易爆；氫氟烴類工質(HFCs)的全球變暖潛能值（GWP）偏高，不符合《〈蒙特利爾議定書〉基加利修正案》的要求[2]；氫氟烯烴類工質(HFOs)汽化潛熱低，熱力性能較差。因此很難找到性能優異、熱穩定性高又對環境友好安全的純工質，而CO2混合工質動力循環通過不同的組分配比，可以實現既能提高能源利用效率又能對環境友好的要求。CO2+HCs混合工質可以消除可燃性，降低易燃易爆的風險；CO2+HFCs混合工質可以降低GWP 值；CO2+HFOs混合工質可以提高工質的熱力性能。

Bell 等[3]總結了目前能作為制冷工質的各類混合工質的氣液相平衡實驗數據，其中CO2混合工質的氣液相平衡實驗數據有21 種。混合工質氣液相平衡理論計算的重點是選擇合適的混合規則，提高計算精度，混合規則的形式各種各樣，包括vdW 混合規則、WS 混合規則、MHV1 混合規則等。各類CO2混合工質適用于何種混合規則沒有明確給出，Kim等[4]對于CO2+propane混合工質使用vdW 混合規則 計 算；對 于CO2+HFC?134a 混 合 工 質，Duran?Valencia 等[5]使用vdW 混合規則計算氣液相平衡數據，而Lim 等[6]和Silva?Oliver 等[7]選擇的是WS 混合規 則；Valtz 等[8?9]選 擇WS 混 合 規 則 計 算 了CO2+HFC?152a 和CO2+HFC?227ea 的氣液相平衡數據；對于CO2+HFO?1234yf 混合工質，Juntarachat 等[10]使用vdW 混合規則計算氣液相平衡數據；而Wang等[11]使用了WS 混合規則計算了CO2+HFO?1234ze(E)混合工質。綜上所述，目前各類CO2混合工質使用的混合規則形式不統一，系統分析不同混合規則形式對于各類CO2混合工質的適用性具有重要意義。另外，由于氣液相平衡實驗需要耗費大量的時間和精力，因此提出適合CO2混合工質氣液相平衡的預測模型，也是一項具有重大意義的工作。

本文結合CO2混合工質動力循環的應用背景，選取出11 種CO2混合工質的氣液相平衡實驗數據，包括7 種CO2+HFCs/HFOs 二元體系和4 種CO2+HCs二元體系。選用PR 方程結合三種不同形式的混合規則(vdW、WS、MHV1)計算這11 種二元體系的氣液相平衡數據，并與各體系公開發表的文獻實驗值相比較，得出計算結果與實驗數據之間的相對偏差，從而分析與討論不同混合規則對于各類CO2混合工質的適用性。最后提出一種差值模型，對CO2混合工質氣液相平衡模型進行了預測。

1 氣液相平衡的理論基礎公式

氣液相平衡性質反映的是流體氣液兩相之間的相互關系，描述的是混合物氣液平衡時，其溫度T、壓力p和氣相組分yi、液相組分xi之間的關系。氣液相平衡理論需要狀態方程結合混合規則與活度系數模型描述得出。

1.1 PR狀態方程

本文采用PR 方程作為計算混合工質氣液相平衡的基礎。PR 方程由Peng 和Robinson 在1976 年提出[12]，其表達式為：

式中，a和b是狀態方程系數，可由臨界參數和偏心因子計算得到：

α(T)函數定義為：

狀態方程被用于描述混合工質的熱力學性質時，需要引入混合規則。

1.2 混合規則

1.2.1 vdW 混合規則 描述混合工質的熱力學性質時，vdW[13]混合規則被廣泛應用，作為一種常數型混合規則，其形式如下：

1.2.2 MHV1 混合規則 MHV1 混合規則是Michelsen[14]對HV[15]混合規則進行改進得到的，其是一種經典的GE?EOS 混合規則,以零壓為參考態，具體形式如下：

1.2.3 WS 混 合 規 則 Wong 和Sandler 于1992 年 提出Wong?Sandler(WS)混合規則[16]，該法則基于超額Helmholtz 自由能，與密度無關且滿足第二維里系數邊界條件,其形式如下：

三種混合規則各有優勢，WS和MHV1混合規則對高度非理性體系和強極性體系有著非常好的描述能力；而vdW 混合規則只有一個交互參數，形式簡單，具有較強的物理意義。

1.3 活度系數模型

1968 年Renon 和Prausnitz[17]修正了溶液局部組成表達式, 在雙流體理論的基礎上提出了NRTL 活度系數模型,其能很好地描述二元體系的相平衡性質,對于常見的二元系統其表達式為:

組分1、2的活度系數方程為:

2 混合物氣液相平衡計算

三種混合規則中的相互作用參數需要結合實驗數據，使用目標函數優化計算得出。本文的目標函數綜合考慮了壓力計算相對偏差與氣相組分濃度絕對偏差，如下所示：

本文分別在CO2/HFCs、CO2/HFOs 和CO2/HCs 二元體系中選取了CO2(1)+HFC?152a(2)[9]、CO2+HFO?1234ze(E)[11]、CO2+propane[4]三類混合工質，給出了不同混合規則的比較。

2.1 CO2+HFC-152a混合工質各種混合規則對比

CO2(1)+ HFC?152a(2)混合工質的壓力相對偏差、氣相摩爾分數絕對偏差以及p-x-y曲線如圖1、圖2 所示，其中圖2 中的散點表示不同溫度下的氣液相平衡的實驗點，從圖中可以看出對于CO2(1)+HFC?152a(2)混合工質，WS+NRTL 混合規則壓力相對偏差的絕對值在1%以內，氣相摩爾分數絕對偏差的絕對值在0.01 以內，其計算精度最高；在低于CO2臨界溫度304.13 K 的溫區(亞臨界區域)，三種混合規則計算精度相差不大，但在高于CO2臨界溫度的溫區(超臨界區域)，WS+NRTL 混合規則的優勢就體現了出來，這是因為WS 混合規則對高度非理性體系和強極性體系有著非常好的描述能力。

圖1 CO2(1)+HFC?152a(2)的壓力相對偏差與氣相摩爾分數絕對偏差Fig.1 The relative deviation of the pressure and the absolute deviation of the component mole fraction of CO2(1)+HFC?152a(2)

圖2 CO2(1)+HFC?152a(2)的p?x?y曲線Fig.2 p?x?y diagram of CO2(1)+HFC?152a(2)

2.2 CO2+HFO-1234ze（E）混合工質各種混合規則對比

從圖3、圖4 中可以看出對于CO2(1)+HFO?1234ze(E)(2)混合工質，WS+NRTL 混合規則效果最好，MHV1+NRTL 次之，vdW 最差；在低于CO2臨界溫度區域(亞臨界區域)三種混合規則計算精度相差不大，但在高于CO2臨界溫度區域(超臨界區域)，WS+NRTL 混合規則的優勢較高，MHV1+NRTL 的計算精度也較高，誤差較大的地方其實在液相區，這是因為對于MHV1 混合規則，液相摩爾體積是一個恒定的常數導致的。

圖3 CO2(1)+HFO?1234ze(E)(2)的壓力相對偏差與氣相摩爾分數絕對偏差Fig.3 The relative deviation of the pressure and the absolute deviation of the component mole fraction of CO2(1)+HFO?1234ze(E)(2)

圖4 CO2(1)+HFO?1234ze(E)(2)的p?x?y曲線Fig.4 p?x?y diagram of CO2(1)+HFO?1234ze(E)(2)

2.3 CO2+propane混合工質各種混合規則對比

從圖5、圖6 中可以看出對于CO2(1)+propane(2)混合工質，vdW 混合規則的壓力偏差和氣相摩爾分數絕對偏差與WS 和MHV1 混合規則相差不大，甚至在高溫區(高于CO2臨界溫度304.13 K)更具優勢，此時vdW 混合規則形式簡單的優勢突顯出來，這是因為該二元體系非理想性不強導致的。本文后續會繼續介紹CO2(1)+烷烴(2)混合工質使用vdW 混合規則的優勢，并通過已知溫區下的二元交互作用參數推算未知溫區下的二元交互作用參數，從而推算未知溫區下的氣液相平衡曲線。

圖5 CO2(1)+propane(2)的壓力相對偏差與氣相摩爾分數絕對偏差Fig.5 The relative deviation of the pressure and the absolute deviation of the component mole fraction of CO2(1)+propane(E)(2)

圖6 CO2(1)+propane(2)的p?x?y曲線Fig.6 p?x?y diagram of CO2(1)+propane(2)

2.4 CO2混合工質各種混合規則對比

采用PR+WS+NRTL 模型、PR+MHV1+NRTL 模型和PR+vdW 模型對7種CO2+HFCs/HFOs 二元混合工質氣液相平衡性質進行了計算。計算結果見表1，包括了實驗點數和實驗溫區，采用PR+WS+NRTL模型，AARD(p)為0.61%，AAD(y)為0.0055；采用PR+MHV1+NRTL 模 型，AARD(p)為1.14%，AAD(y)為0.0091；采用PR+vdW 模型，AARD(p)為1.41%，AAD(y)為0.0110。3 種模型對比可知，對于CO2+HFCs/HFOs 二元體系，相比于vdW 混合規則，MHV1 混合規則對計算精度提升有限，WS 混合規則對計算精度提升明顯。

表1 三種混合規則對CO2+HFCs/HFOs混合工質的氣液相平衡計算偏差Table 1 Calculation deviations in the vapor-liquid equilibrium of the CO2+HFCs/HFOs mixtures by three mixing rules

采用三種模型對4種CO2+HCs二元混合工質氣液相平衡性質進行了計算。計算結果見表2，采用PR+WS+NRTL 模型，AARD(p)為0.57%，AAD(y)為0.0060；采 用PR+MHV1+NRTL 模 型，AARD(p)為0.73%，AAD(y)為0.0068；采用PR+vdW 模型，AARD(p)為1.01%，AAD(y)為0.0066。3 種模型對比可知，對于CO2+HCs 二元混合工質，相比于vdW 混合規則，MHV1 與WS 混合規則對計算精度提升有限，特別是對氣相組分，計算精度基本沒有提升。

表2 三種混合規則對CO2+HCs二元混合工質氣液相平衡計算偏差Table 2 Calculation deviations in the vapor-liquid equilibrium of the CO2+HCs mixtures by three mixing rules

對于CO2+烷烴類二元混合工質，在vdW 混合規則下，利用實驗數據回歸計算各個溫區下的k12如表3 所 示，從 表 中 可 以 看 出，CO2+propane[4]、CO2+n?butane[20?21]、CO2+isobutane[22?23]三 種 二 元 混 合 工 質 的二元交互作用參數k12都在0.12～0.13 附近。將這些溫區下的二元交互作用參數取平均值與先前優化出來的二元交互作用參數做對比，分別計算各個溫區下的AARD(p)和AAD(y)的值，計算結果如表4所示。

表3 CO2+烷烴類混合工質的二元交互作用參數的最優值（vdW混合規則）Table 3 Optimal value of binary interaction parameters for CO2+alkane mixtures（vdW mixing rules）

2.5 CO2+小分子烷烴混合工質的氣液相平衡的預測

從表4 中可以看出，對于CO2+小分子烷烴體系，平均二元交互作用參數與優化的二元交互作用參數相比，計算精度相差不大。因此對于這類混合工質可以用平均二元交互作用參數代替使用實驗數據優化的二元交互作用參數，并將這一結論擴展到其他沒有計算的溫區，利用平均二元交互作用參數得到的氣液相平衡數據與實驗值比較，如圖7所示。

表4 CO2+小分子烷烴二元混合物的優化結果與預測結果的比較Table 4 Comparison of the optimization results and prediction results of CO2+small molecule alkane binary mixture

如圖7所示，使用平均二元交互作用參數，計算330 K 下 的CO2+propane[26]、387 K 下 的CO2+n?butane[20]和320 K 下的CO2+isobutane[27]，計算結果與實驗點擬合較好，進一步驗證了使用平均二元交互作用參數的有效性。這一發現具有重要的意義，對于CO2+小分子烷烴二元體系，通過少量的實驗，得出二元交互作用參數的規律，預測沒有實驗數據的溫區下的氣液相平衡，可以節省大量實驗所需要的時間與精力。

圖7 vdW混合規則預測的p?x?y曲線Fig.7 p-x-y diagram predicted by vdW mixing rules

2.6 CO2+HFCs/HFOs 混合工質的氣液相平衡的預測

Hu 等[28]指出，WS 與MHV1 這類多參數狀態方程很難預測混合工質氣液相平衡性質，然而如果狀態方程中二元交互作用參數可以預測，則大大降低實驗耗費的精力。表5 列出了CO2+HFCs/HFOs 二元體系，使用PR+vdW 模型回歸計算的各個溫區下的k12的值。目前關于混合工質氣液相平衡的預測主要集中在半經驗公式的預測、基團貢獻法的預測以及理論模型的預測。結合Hu 等[28]、Chen 等[29]、Zhang 等[30]以及Duarte 等[31]對二元交互作用參數k12的處理，由于CO2與含氟制冷劑的特殊性質，提出一種CO2與含氟制冷劑的差值計算模型[式(21)～式(23)]，結合CO2和制冷劑的物性參數以及氟原子數量，具有一定的物理意義。

表5 CO2+HFCs/HFOs混合工質的二元交互作用參數的最優值（vdW混合規則）Table 5 Optimal value of binary interaction parameters for CO2+HFCs/HFOs（vdW mixing rules）

這個差值模型，只與CO2的臨界溫度、臨界壓力和偏心因子以及含氟制冷劑的偏心因子和氟原子數量有關，是一個完全預測性的模型。以該模型計算得到的二元交互作用參數與前面優化結果進行對比，計算結果如表6所示。

由表6 可知，優化計算得到的AARD(p)值為1.17%，AAD(y)值為0.0078，而預測的AARD(p)值為2.03%，AAD(y)值為0.0120，所以可以看出該預測模型的精度很高。由于CO2+含氟制冷劑的氣液相平衡實驗數據較少，這一差值模型的提出對于處理其他沒有實驗數據的混合工質具有極其重要的意義。

表6 CO2+HFCs/HFOs混合工質的優化結果與預測結果的比較Table 6 Comparison of the optimization results and prediction results of CO2+HFCs/HFOs

3 結 論

本文將PR 方程與vdW、MHV1、WS 三種混合規則相結合，建立了CO2混合工質氣液相平衡計算模型，計算了11 種CO2混合工質的相平衡性質，包括7種CO2+HFCs/HFOs二元體系和4種CO2/HCs二元體系。主要結論如下。

（1）對于CO2+HFCs/HFOs 二元體系，相比于vdW 混合規則，MHV1 混合規則對計算精度提升有限，WS 混合規則對計算精度提升明顯，因此需要通過形式較為復雜的MHV1 或WS 混合規則對CO2+HFCs/HFOs 的氣液相平衡性質進行關聯和預測。

（2）對于CO2+HCs 二元體系，因HC 分子極性弱于HFC 和HFO，其二元混合物非理想性不強，因此相比于vdW 混合規則，MHV1 與WS 混合規則對氣液相平衡性質的計算精度提升有限，特別是對氣相組分，計算精度基本沒有提升。

（3）對于CO2+小分子烷烴二元體系，可以通過PR方程結合vdW混合規則，通過部分等溫線下的實驗數據獲得平均二元交互作用參數，對于CO2+烷烴混合工質平均二元交互作用參數在0.12～0.13附近，可以通過平均二元交互作用參數對其他溫區下的氣液相平衡進行預測。

（4）對于CO2+含氟制冷劑二元體系，提出了一種差值模型預測其氣液相平衡性質，預測得到的AARD(p)值為2.03%，AAD(y)值為0.0120，預測效果較好，下一階段的研究目標是擴展這一差值模型的適用工質范圍。

符 號 說 明

AARD(p)——壓力相對偏差

AAD(y)—— 氣相摩爾分數絕對偏差

ai,bi——純物質的狀態參數

aij——混合工質組分i和組分j的相互作用項

am,bm——分別代表引力項和協體積項，是混合工質的狀態參數

C——狀態方程有關的常數，對于PR方程，C=?0.62323

GE——無窮壓力下的超額Helmholtz自由能

kij,k12——二元交互作用參數，且kij=kji,k12=k21

k1——CO2的混合因子

k2——含氟制冷劑的混合因子

pc——臨界壓力，Pa

pexp,pcal——分別為實驗壓力值和計算壓力值，Pa

q1——與狀態方程有關的常數，對于PR方程，q1=?0.528

R——通用氣體常數,J/(mol·K)

Tc——臨界溫度，K

Tr——對比溫度,K

v——體積，m3/mol

x1,y1——分別為組分1 的液相摩爾分數和氣相摩爾分數

yexp,ycal——分別為實驗氣相摩爾分數和計算氣相摩爾分數

τ12,τ21,α12,α21——可調參數, 一般由氣液相平衡實驗數據擬合得到。α12與溶液組成和溫度無關, 只取決于溶液的類型,本文取0.30,且α12=α21

ω——偏心因子