1 前言

氯氟烃和氢氯氟烃等传统制冷剂，因具备较高全球变暖潜能值(GWP)和臭氧消耗值(ODP)逐渐被新兴制冷剂淘汰^[1]。天然制冷剂CO₂，无毒不易燃，化学性质稳定，单位容积制冷量大^[2-4]，其ODP为零且GWP极小^[5-6]。离子液体(ILs)作为由不同阴阳离子组成的绿色溶液，具有较高的化学稳定性、热稳定性、气体溶解能力和分子结构可设计性^[7-8]，研究发现多数离子液体对CO₂有较高吸收能力^[9-11]。Morgan等^[12-13]研究发现温度为293.15~363.15 K、压力为0.1~50 MPa时，CO₂在[emim][Tf₂N]中的溶解度高达0.8 mol·L⁻¹。He等^[14]研究发现温度为283.15~363.15 K、压力为0.2~5 MPa时，CO₂在[BMP][Tf₂N]中的摩尔溶解度最高可达0.532 4 mol·L⁻¹。众多研究表明ILs可作为良好的CO₂吸收剂，应用于吸收式制冷循环系统装置中^[15-16]。

很多学者对不同CO₂-ILs体系的汽液相平衡进行实验研究，得到不同温度、压力下CO₂-ILs体系的相平衡数据^[17-19]。以相平衡数据为基础，选择合适的模型可预测热力学性质^[20-21]。王卫^[22]测定了CO₂-[omim][BF₄]、CO₂-[omim][Tf₂N]的汽液相平衡数据，通过UNIFAC模型预测了相平衡时CO₂的溶解度，平均相对偏差分别为2.54%、3.87%。何丽娟等^[3]用PR+NRTL方法关联Schilderman等的实验数据，得出CO₂和[emim][Tf₂N]二元混合体系在高温高压下的超额混合焓，表明CO₂液相摩尔分数大于0.43时体系的超额混合焓均为负值，混合为放热过程。胡艺等^[23]通过实验测定离子液体的吸收焓，发现CO₂与离子液体之间的相互作用受二者之间形成的键长、键角影响。Safarov等^[24]测定了CO₂在离子液体[emim][CH₃SO₃]中的溶解度，根据亨利定律常数关联出混合物的溶解吉布斯自由能、溶解焓和溶解熵。溶解焓和溶解熵均为负值，表明CO₂溶于离子液体时存在放热过程且溶液有序性较高。众多研究表明，热力学性质对于研究CO₂及其他气体在离子液体中的溶解度具有重要的参考价值。

目前许多学者致力于研究CO₂-ILs体系的吉布斯自由能等性质，但相关的超额性质研究较少。为了进一步研究不同类别的CO₂-ILs二元体系相平衡特性和热力学性质，本研究选取阴离子相同、阳离子烷基链长度不同的4种离子液体([C₂-mim][Tf₂N]、[C₄-mim][Tf₂N]、[C₆-mim][Tf₂N]、[C₈-mim][Tf₂N])，整理了Shin等^[17]的实验数据，采用SRK状态方程+WS混合规则+UNIFAC活度系数模型，探究了温度、压力、CO₂溶解度和ILs阳离子烷基链长度对混合体系超额吉布斯自由能、超额焓、超额熵的影响。

2 热力学模型 2.1 SRK状态方程

SRK状态方程^[25]如下所示：

$ p = \frac{{RT}}{{{V_{\text{m}}} - b}} - \frac{a}{{{V_{\text{m}}}({V_{\text{m}}} + b)}} $

(1)

$ a = {a_{\text{c}}}\alpha ({T_{\text{r}}},\omega ) $

(2)

$ {a_{\text{c}}} = {0.427_{}}48\frac{{{R^2}T_{\text{c}}^2}}{{{p_{\text{c}}}}} $

(3)

$ b = {0.086_{}}64\frac{{R{T_{\text{c}}}}}{{{p_{\text{c}}}}} $

(4)

$ {\alpha ^{0.5}} = 1 + \left( {0.48 + 1.574\omega - 0.176{\omega ^2}} \right)\left( {1 - T_{\text{r}}^{0.5}} \right) $

(5)

式中：V_m为摩尔体积，m³·mol⁻¹；p为压力，Pa；$ {p_{\text{c}}} $为临界压力，MPa；T为温度，K；a、b为引力参数和斥力参数；$ {a_{\text{c}}} $为临界温度与临界压力的函数表达式；R为摩尔气体常数，J·(mol·K)⁻¹；$\alpha ({T_{\text{r}}},\omega )$为温度和偏心因子的函数；T_r为温度的参数；ω为偏心因子；T_c为临界温度，K。临界性质和偏心因子列在表 1中。

2.2 WS混合规则

WS混合规则^[26]如下所示：

$ {a_{\text{m}}} = {b_{\text{m}}}RTD $

(6)

$ {b_{\text{m}}} = \frac{Q}{{1 - D}} $

(7)

$ D = \sum\limits_i {\frac{{{a_i}}}{{{b_i}RT}}} + \frac{1}{C}\frac{{{G^{\text{E}}}}}{{RT}} $

(8)

$ Q = \sum\limits_i {\sum\limits_j {{x_i}{x_j}{{(b - \frac{a}{{RT}})}_{ij}}} } $

(9)

$ {(b - \frac{a}{{RT}})_{ij}} = \frac{1}{2}({b_i} + {b_j}) - \frac{{\sqrt {{a_i}{a_j}} }}{{RT}}(1 - {k_{ij}}) $

(10)

式中：D为介电常数；$ {a_{\text{m}}} $、$ {b_{\text{m}}} $分别代表混合组分的引力和斥力参数；$ {a_i} $、a_j为组分i、j的引力参数；$ {b_i} $、b_j为组分i、j的斥力参数；$ {x_i} $、$ {x_j} $为组分i、j的摩尔分数；Q、C为混合规则常数；k_ij为交互作用参数；G^E为超额吉布斯自由能，J·mol⁻¹。

2.3 UNIFAC活度系数模型

Fredenslund等^[27]基于量子化学理论(UNIQUAC)和基团贡献概念提出UNIFAC模型。Lei等^[28]对UNIFAC模型进行了改进。该模型由组合项和剩余项组成，其中组合项$ \ln \gamma _i^{\text{c}} $主要反映分子大小和形状的差别，剩余项$ \ln \gamma _i^{\text{R}} $表示基团之间相互作用的影响，表示如下：

$ \ln {\gamma _i} = \ln \gamma _i^{\text{c}} + \ln \gamma _i^{\text{R}} $

(11)

$ 组合项：\ln \gamma _i^{\text{c}} = \ln \frac{{{\varphi _i}}}{{{x_i}}} + \frac{Z}{2}{q_i}\ln \frac{{{\theta _i}}}{{{\varphi _i}}} + {l_i} - \frac{{{\varphi _i}}}{{{x_i}}}\sum\limits_j {{x_j}} {l_j} $

(12)

$ {l_i} = \frac{Z}{2}({r_i} - {q_i}) - ({r_i} - 1){,_{}}Z = 10 $

(13)

$ {\theta _i} = \frac{{{q_i}{x_i}}}{{\sum\limits_j {{q_j}{x_j}} }} \text{ } {\varphi _i} = \frac{{{v_i}{x_i}}}{{\sum\limits_j {{v_j}{x_j}} }} $

(14)

$ {r_i} = \sum\limits_k {v_k^{(i)}} {R_k} \text{ } {q_i} = \sum\limits_k {v_k^{(i)}} {Q_k} $

(15)

式中：$ {\gamma _i} $为归一化活度系数；r_i和q_i分别为纯物质分子的van der Waals体积和表面积；θ_i和φ_i为i组分的表面积分数和体积分数；$v_k^{(i)}$为i中k基团的个数；l_i为i组分的常数；Z为晶格配位数；R_k、Q_k为基团k的体积参数和表面积参数。

$ 剩余项：\ln \gamma _i^{\text{R}} = \sum\limits_k {v_k^{(i)}} (\ln {\mathit{\Gamma } _k} - \ln \mathit{\Gamma } _k^{(i)}) $

(16)

$ \ln {\mathit{\Gamma } _k} = {Q_k}[1 - {\ln _{}}(\sum\limits_m {{\mathit{\Theta } _m}{\mathit{\Psi } _{mk}}} ) - \sum\limits_m {(\frac{{{\mathit{\Theta } _m}{\mathit{\Psi } _{mk}}}}{{\sum\limits_n {{\mathit{\Theta } _n}{\mathit{\Psi } _{nm}}} }}} )] $

(17)

$ \ln \mathit{\Gamma } _k^{(i)} = {Q_k}[1 - {\ln _{}}(\sum\limits_m {\mathit{\Theta } _m^{(i)}} {\mathit{\Psi } _{mk}}) - \sum\limits_m {\frac{{\mathit{\Theta } _m^{(i)}{\mathit{\Psi } _{km}}}}{{\sum\limits_n {\mathit{\Theta } _m^{(i)}{\mathit{\Psi } _{nm}}} }}} ] $

(18)

式中：Г_k为k基团的剩余活度系数，$ \mathit{\Theta } _m^{(i)} $为i组分中k基团的表面积分数，$ {\mathit{\Psi } _{mk}} $为基团m与基团k的相互作用项，定义为

$ {\mathit{\Psi } _{nm}} = {\exp _{}}(\frac{{ - {a_{mn}}}}{T}) $

(19)

式中：$ {a_{mn}} $为基团m和基团n之间的相互作用参数。

UNIFAC活度系数模型用于离子液体，目前常用的拆分方法有3种。本研究采用Lei等^[28]的拆分方式，将离子液体[C₂-mim][Tf₂N]的阴阳离子看作一个整体(即[−mim][Tf₂N]，这样可忽略阴阳离子间的静电作用)，将阳离子的烷基链(丁基)拆成1个─CH₃，1个─CH₂ ([C₄-mim][Tf₂N]、[C₆-mim][Tf₂N]、[C₈-mim][Tf₂N]拆分方法类似，拆分示意图如图 1所示)。各个基团间的相互作用参数采用Lei等^[28]回归得到的UNIFAC活度系数模型相互作用参数。

2.4 计算流程模型

对于二元混合体系汽液相平衡准则：

$ \hat f_i^{\text{v}} = \hat f_i^{\text{l}} $

(20)

汽相和液相组分逸度分别用组分逸度系数和活度系数来计算，平衡关系式如下：

$ p\hat \varphi _i^{\text{v}}{y_i} = f_i^{\text{l}}{x_i}{\gamma _i} $

(21)

式中：$ \hat \varphi _i^{\text{v}} $为组分i的气相逸度系数，$ {y_i} $为组分i的气相摩尔分数，$ {\gamma _i} $为液相活度系数，$ \hat f_i^{\text{v}} $为组分i的气相逸度，$ \hat f_i^{\text{l}} $为组分i的液相逸度，$ f_i^{\text{l}} $为纯组分i的逸度，$ f_i^{\text{l}} \approx p_i^{\text{s}}{\exp _{}}\left[ {\frac{{V_i^{\text{l}}(p - p_i^{\text{s}})}}{{RT}}} \right] = p_i^{\text{s}}{\mathit{\Phi } _i} $，其中$ {\mathit{\Phi } _i} $为Poynting因子；$ p_i^{\text{s}} $为组分i的饱和蒸气压，可由Shiflett和Yokozeki^[29]建议的安托因方程求得V_i^l为组分i的液相摩尔体积。

本研究采用状态方程+活度系数(EOS+γ)法计算了4组CO₂-ILs体系的相平衡数据。图 2是EOS+γ法计算CO₂-ILs体系相平衡数据的流程图。通过迭代计算压力值，使气相组成归一化条件成立，结束迭代，输出结果。

通过整理Shin等^[17]的汽液相平衡实验数据，计算了4组CO₂-ILs体系的相平衡数据，并关联了该体系的超额热力学性质。4组CO₂-ILs二元体系相平衡实验数据及计算数据见表 2。表中p^cal是通过迭代气泡点压力得到的计算值，p^exp是实验值，x为CO₂在液相中的摩尔分数。将G^E的定义与$ {\gamma _{\text{i}}} $关联：

$ \frac{{{G^{\text{E}}}}}{{RT}} = \sum\limits_i {{x_i}\ln {\gamma _i}} $

(22)

由G^E，可得到超额焓H^E与超额熵S^E：

$ {H^{\text{E}}} = - {T^2}\left( {\frac{{\partial \left( {\frac{{{G^{\text{E}}}}}{T}} \right)}}{{\partial T}}} \right) $

(23)

$ {S^{\text{E}}} = \frac{{{H^{\text{E}}} - {G^{\text{E}}}}}{T} = - \left( {\frac{{\partial {G^{\text{E}}}}}{{\partial T}}} \right) $

(24)

3 结果与讨论 3.1 超额热力学性质 3.1.1 超额吉布斯自由能G^E的结果分析

利用SRK状态方程、WS混合规则以及UNIFAC活度系数模型得到G^E的计算结果如图 3所示。

由图 3可知：当CO₂液相摩尔分数为定值时，CO₂-[C_n-mim][Tf₂N](n=2, 4, 6, 8)型二元混合体系的G^E分别随温度、压力和ILs阳离子烷基链长度增大而减小。其中压力对G^E影响较大。当x=0.233、ILs阳离子碳基链长度n=2时，G^E出现最大值，为−1 066.88 J·mol⁻¹。

随着CO₂液相摩尔分数的增加，G^E受温度影响减弱。ILs密度随温度的升高而减小，自由体积增大，吸附CO₂的能力增强，同时自由体积随CO₂液相摩尔分数的增大而减小，G^E的变化趋于平缓。随着ILs烷基链长度递增，自由体积增大，G^E的数值随之增大。

3.1.2 超额焓H^E的结果分析

利用SRK状态方程、WS混合规则以及UNIFAC活度系数模型得到H^E的计算结果如图 4所示。

由图 4可知，当CO₂摩尔分数为定值时，CO₂-[C_n-mim][Tf₂N](n=2, 4, 6, 8)型二元混合体系的超额焓H^E受温度和压力的影响较小。H^E随ILs阳离子烷基链长度、CO₂摩尔分数的增大而减小，其中CO₂摩尔分数对混合体系超额焓影响较大。当x=0.233、ILs阳离子烷基链长度n=2时，H^E出现最大值，为407.24 J·mol⁻¹。

H^E受黏度影响较大，混合体系的黏度由分子间氢键和范德华力决定，随着ILs阳离子烷基链长度和CO₂液相摩尔分数增大，烷基链与CO₂之间形成新的分子间作用力，H^E随之减小。

3.1.3 超额熵S^E的结果分析

利用SRK状态方程、WS混合规则以及UNIFAC活度系数模型得到S^E的计算结果如图 5所示。

由图 5可知，当CO₂摩尔分数为定值时，S^E随ILs阳离子烷基链长度的增大而增大，随压力的增大而减小。随着CO₂摩尔分数增大，温度对超S^E的影响逐渐减小。当x=0.301 9、ILs阳离子烷基链长度n=8时，S^E出现最大值，7.88 J·mol⁻¹·K⁻¹。

S^E代表真实溶液与理想溶液的熵的偏差。随着ILs阳离子烷基链长度的增大，混合体系的有序性降低。随着CO₂摩尔分数增大，温度和压力对溶液有序性的影响逐渐减小。

3.2 相互作用参数k_ij回归及模型精确度计算

k_ij影响相平衡模型的精确度，k_ij不同，迭代计算输出压力值也不同。为使模型预测效果更加精准，通过目标函数OF(式(25))，对CO₂-[C_n-mim][Tf₂N](n=2, 4, 6, 8)型二元混合体系的k_ij进行回归。当OF计算结果取最小值时，输出此时的k_ij。计算结果显示，ILs阳离子烷基链长度n=2, 4, 6, 8时，二元混合体系的k_ij分别为0.42、0.51、0.60和0.62时，模型的计算误差最小，此时压力p的平均相对偏差ARD(式(26))如表 3所示。随着温度升高，CO₂-ILs体系压力p的平均相对偏差ARD逐渐增大，模型预测精度逐渐降低；当温度为定值时，随着ILs阳离子烷基链的增加，压力p的ARD逐渐减小，模型预测精度逐渐升高。同时随着ILs阳离子烷基链的增加，CO₂与[C_n-mim][Tf₂N]的相互作用参数逐渐增大。

$ \text{OF} = \sum\limits_{q = 1}^{{N_{\text{E}}}} {{{\left( {\frac{{p_q^{\exp } - p_q^{\text{cal}}}}{{p_q^{\exp }}}} \right)}^{\text{2}}}} $

(25)

$ \text{ARD} = \frac{1}{N}\sum\limits_{q = {\text{1}}}^N {\frac{{\left| {p_q^{\exp } - p_q^{\text{cal}}} \right|}}{{p_q^{\exp }}}} $

(26)

式中：N_E为系统实验数据组数。

由表 3知，SRK+WS+UNIFAC模型计算CO₂-[C_n-mim][Tf₂N](n=2, 4, 6, 8)型二元混合体系具有较高精度。ILs阳离子烷基链长度越长，SRK+WS+UNIFAC模型的计算精度越高，同时该模型的精度受压力影响较大。4种离子液体的极性都随着阳离子上烷基链长度的增加而减小^[30]，因此SRK+WS+UNIFAC模型更适合计算温度压力较低、极性较弱的CO₂-ILs体系。

4 结论

采用选定状态方程+活度系数法计算了CO₂-[C_n-mim][Tf₂N](n=2, 4, 6, 8) 4种混合工质对的相平衡，得出以下结论：

(1) 超额吉布斯自由能G^E随温度、压力、ILs烷基链长度增大呈下降趋势，受压力影响变化较大。

(2) 超额焓H^E受温度、压力的影响较小，受CO₂摩尔分数影响较大，随ILs阳离子烷基链长度的增加而减小。

(3) 超额熵S^E随ILs阳离子烷基链长度的增大而增大，随着压力的升高而减小。

(4) SRK+WS+UNIFAC模型对CO₂-[C_n-mim][Tf₂N](n=2, 4, 6, 8)型二元混合体系具有较高计算精度。