1 前言

经典成核理论(classical nucleation theory，CNT)忽略了表面张力的曲率修正(修正首项由Tolman长度表示)^[1]，成核速率的计算结果与实验值存在明显偏差^[2]。目前，Tolman长度还不能实验测定，通常采用密度泛函理论(density functional theory，DFT)或梯度平方理论(square gradient theory，SGT)进行理论计算^[3-4]。但是，DFT和SGT只适用于远离临界点的两相平衡区，在临界点附近两者都失效。Blokhuis等^[5]将Tolman长度与流体的容积性质相联系，使得可以用状态方程预测Tolman长度。但是常用的状态方程没有考虑临界密度涨落对热力学性质的影响，因此需要采用重正化群理论。但是由White建立的全局重正化群理论(global renormalization group，GRG)^[6]及其推广^[7]没有正确反映流体的临界非对称性^[8-9]，不能用来计算和非对称性有关的Tolman长度^[10-11]。最近，作者改进了GRG理论^[12]，并将改进的GRG理论与阱宽可变的统计缔合流体(statistical associating fluid theory for molecules with variable range，SAFT-VR)方程^[13]结合，计算了几种简单有机分子的临界表面张力和Tolman长度^[14]。然而，SAFT-VR方程对于强极性和缔合流体不够准确，而微扰链-SAFT方程(Perturbed-chain SAFT, PC-SAFT)^[15-16]更严格地考虑了分子间的色散相互作用，因此本研究将改进的GRG方法与PC-SAFT方程结合，用于醇、醋酸、水等工业上有着广泛应用的缔合流体，为修正经典成核理论提供理论依据。

2 理论和计算方法 2.1 重正化计算方法

原始全局重正化群虽然能正确描述流体的多数临界性质，但是未能正确描述临界非对称性，因此不能用来计算Tolman长度。场论分析^[17]和实例计算^{[12, 14]}都表明，为正确描述临界非对称性，需将原始GRG的哈密顿量改为如下形式：

$ H = \int {\left\{ {{f_{^{{\text{EOS}}}}}n({\mathit{\boldsymbol{r}}}) + [{a_{{\text{vdw}}}} + {u_1}{k_{\text{B}}}Tn({\mathit{\boldsymbol{r}}})]L_0^2{{\left| {\nabla n({\mathit{\boldsymbol{r}}})} \right|}^2}} \right\}} {\text{d}}r $

(1)

式中：H为哈密顿量；f为自由能密度，是摩尔Helmholtz自由能与摩尔体积之比；f_EOS为由状态方程计算的自由能密度；n为分子数密度，n(r)表示n是位置矢量r的函数；$ \nabla $为梯度算符；k_B为玻尔兹曼常数；T为热力学温度；a_vdw为分子间吸引作用的van der Waals参数；L₀为截断长度；u₁为待定参数。

由于临界非对称性由分子间相互作用决定，因此，如果在建立状态方程时准确考虑了分子间相互作用，则完整标度理论^[18]中表征临界非对称性的关键参数——场混合系数a₃就可由状态方程直接求出^[10]。因此本研究采用这样的方法确定u₁的值：选取合适的u₁，使得由式(1)结合重正化群理论计算得到的a₃与直接由f_EOS预测的a₃一致。u₁的这种取法使得式(1)从临界非对称的角度看是自洽的。由重正化群理论计算a₃实质上就是计算气液平衡压力和定容热容，在临界点附近二者存在如下所示的线性关系^[10]：

$ \hat P'' \approx \frac{1}{{1 + {a_3}}}\frac{{\hat C_{V, {\text{m}}}^{{\text{res}}}}}{{\hat T}} $

(2)

式中：$ \hat p $、$ \hat T $、$ \hat C_{V, {\text{m}}}^{{\text{res}}} $为约化的压力、温度和扣除了理想气体贡献的摩尔定容热容，由$ \hat p = {p \mathord{\left/ {\vphantom {p {({n_{\text{c}}}{k_{\text{B}}}{T_{\text{c}}})}}} \right. } {({n_{\text{c}}}{k_{\text{B}}}{T_{\text{c}}})}} $，$ \hat T = {T \mathord{\left/ {\vphantom {T {{T_{\text{c}}}}}} \right. } {{T_{\text{c}}}}} $，$ \hat C_{V, {\text{m}}}^{{\text{res}}} = C_{V, {\text{m}}}^{{\text{res}}}/R $计算，其中p、T、$ C_{V, {\text{m}}}^{{\text{res}}} $为对应的未约化物理量；n_c、T_c为临界分子数密度和临界热力学温度；R为摩尔气体常数；$ \hat p'' $表示在两相平衡条件下对温度求二阶导数。利用式(2)对临界区p″和C_{V, m}的计算结果线性回归，可求得a₃。由于式(2)仅在临界点附近成立，本研究将温度区间$ {10^{ - 4}} < \left| {{\Delta }\hat T} \right| < {10^{ - 2}} $($ {\Delta }\hat T = \hat T - 1 $)定义为临界区。计算p″和C_{V, m}时需要首先计算自由能对温度的二阶偏导数，本研究用重正化方法先计算自由能对温度的一阶偏导数(即内能密度$ e = {{\partial ({f \mathord{\left/ {\vphantom {f {({k_{\text{B}}}T)}}} \right. } {({k_{\text{B}}}T)}})} \mathord{\left/ {\vphantom {{\partial ({f \mathord{\left/ {\vphantom {f {({k_{\text{B}}}T)}}} \right. } {({k_{\text{B}}}T)}})} {\partial {{({k_{\text{B}}}T)}^{ - 1}}}}} \right. } {\partial {{({k_{\text{B}}}T)}^{ - 1}}}} $)，再对其进行一次数值差分求出二阶偏导数^[12]。

实际计算中将密度区间[0, n_max] 50 000等分(n_max为体系允许的最大数密度)，计算的是每个等分点上的f和e。详细过程如下：

$ {\tilde f_i}(n) = {\tilde f_{i - 1}}(n) + {\text{δ }}{\tilde f_i}(n){, _{}}{\tilde f_0}(n) = {\tilde f_{^{_{{\text{EOS}}}}}}(n){;_{}}i = 1, 2, \cdots {;_{}}\tilde f(n) = \mathop {\lim }\limits_{i \to \infty } {\tilde f_i}(n), $

(3a)

$ {\text{δ }}{\tilde f_i}(n) = - \frac{1}{{L_i^3}}\ln \frac{{\int_0^{{{\min }_{}}(n, {n_{\max }} - n)} {\exp \left\{ { - L_i^3{{\tilde G}_{i, s}}(n, x)} \right\}{\text{d}}x} }}{{\int_0^{{{\min }_{}}(n, {n_{\max }} - n)} {\exp \left\{ { - L_i^3{{\tilde G}_{i, l}}(n, x)} \right\}{\text{d}}x} }}{, _{}}{L_i} = {2^i}{L_0}, $

(3b)

$ {\tilde G_{i, s}}(n, x) = \frac{{{{\tilde f}_{i - 1}}(n + x) + {{\tilde f}_{i - 1}}(n - x)}}{2} - {\tilde f_{i - 1}}(n) + \left( {\frac{{{a_{{\text{vdw}}}}}}{{{k_{\text{B}}}T}} + {u_1}n} \right)\frac{{{x^2}\mathit{\Delta} }}{{{2^{2i}}}}, $

(3c)

$ {\tilde G_{i, l}}(n, x) = \frac{{{{\tilde f}_{i - 1}}(n + x) + {{\tilde f}_{i - 1}}(n - x)}}{2} - {\tilde f_{i - 1}}(n) + \frac{{{a_{{\text{vdw}}}}{x^2}}}{{{k_{\text{B}}}T}}, $

(3d)

$ {e_i} = \frac{{\partial {{\tilde f}_i}}}{{\partial {{({k_{\text{B}}}T)}^{ - 1}}}} = \frac{{\partial ({{\tilde f}_{i - 1}} + {\text{δ }}{{\tilde f}_i})}}{{\partial {{({k_{\text{B}}}T)}^{ - 1}}}} = {e_{i - 1}} + {\text{δ }}{e_i}{, _{}}{e_0} = {e^{{\text{EOS}}}} = \frac{{\partial {{\tilde f}_{^{{\text{EOS}}}}}}}{{\partial {{({k_{\text{B}}}T)}^{ - 1}}}}{;_{}}i = 1, 2, \cdots {;_{}}e = \mathop {\lim }\limits_{i \to \infty } {e_i}, $

(4a)

$ {\text{δ}}{e_i} = \frac{{\partial ({\text{δ }}{{\tilde f}_i})}}{{\partial {{({k_{\text{B}}}T)}^{ - 1}}}} = \frac{{\int_0^{{{\min }_{}}(n, {n_{\max }} - n)} {{g_{i, s}}{{\exp }_{}}\left\{ { - L_i^3{{\tilde G}_{i, s}}} \right\}{\text{d}}x} }}{{\int_0^{{{\min }_{}}(n, {n_{\max }} - n)} {{{\exp }_{}}\left\{ { - L_i^3{{\tilde G}_{i, s}}} \right\}{\text{d}}x} }} - \frac{{\int_0^{{{\min }_{}}(n, {n_{\max }} - n)} {{g_{i, l}}{{\exp }_{}}\left\{ { - L_i^3{{\tilde G}_{i, l}}} \right\}{\text{d}}x} }}{{\int_0^{{{\min }_{}}(n, {n_{\max }} - n)} {{{\exp }_{}}\left\{ { - L_i^3{{\tilde G}_{i, l}}} \right\}{\text{d}}x} }}, $

(4b)

$ {g_{i, s}} = \frac{{\partial {{\tilde G}_{i, s}}}}{{\partial {{({k_{\text{B}}}T)}^{ - 1}}}} = \frac{{{e_{i - 1}}(n + x) + {e_{i - 1}}(n - x)}}{2} - {e_{i - 1}}(n) + \frac{{{a_{{\text{vdw}}}}{x^2}\mathit{\Delta} }}{{{2^{2i}}}}, $

(4c)

$ {g_{i, l}} = \frac{{\partial {{\tilde G}_{i, l}}}}{{\partial {{({k_{\text{B}}}T)}^{ - 1}}}} = \frac{{{e_{i - 1}}(n + x) + {e_{i - 1}}(n - x)}}{2} - {e_{i - 1}}(n) + {a_{{\text{vdw}}}}{x^2}, $

(4d)

式中：字符顶上的“∼”表示除以k_BT，例如$ {\tilde f_i}(n) = {{{f_i}(n)} \mathord{\left/ {\vphantom {{{f_i}(n)} {({k_{\text{B}}}T)}}} \right. } {({k_{\text{B}}}T)}} $；i表示第i次迭代，f_i、L_i和e_i分别为第i次迭代的自由能密度、截断长度和内能密度；δf表示自由能密度的增量；x为积分变量，n和x均位于等分点上，积分用梯形公式近似计算；Δ是用来描述密度涨落波包平均尺寸的参数。a_vdw和Δ如下计算：

$ {a_{{\text{vdw}}}} = - \frac{{{m^2}}}{2}\int {{u_{{\text{att}}}}(r){\text{d}}r} , \;\;\mathit{\Delta} = \frac{{5\int {{u_{{\text{att}}}}(r){r^2}{\text{d}}r} }}{{3L_0^2\int {{u_{{\text{att}}}}(r){\text{d}}r} }} $

(5)

式中：m为分子的等效链节数；u_att(r)为等效链节间的吸引势，r为位置矢量r的模，是标量。重正化迭代的收敛判据为两次迭代的自由能之间相对误差小于10⁻¹²。

2.2 表面张力和Tolman长度

Tolman^[1]在Gibbs表面热力学理论的框架下推导出了球形液滴的表面张力与曲率半径的关系：

$ \sigma ({R_e}) = {\sigma _\infty }\left( {1 - \frac{{2{\delta _T}}}{{{R_e}}} + \cdots } \right) $

(6)

式中：σ(R_e)为计入了曲率修正的表面张力；R_e为液滴半径；σ_∞为平坦表面张力；δ_T为Tolman长度。Blokhuis和Kuipers给出了球形液滴Tolman长度近似表达式^[5]：

$ - {\hat \delta _T} \approx \frac{{{{\hat \sigma }_\infty }}}{{{{({{\hat n}_{\text{l}}} - {{\hat n}_{\text{g}}})}^2}}}\left( {\hat n_{\text{l}}^2{{\hat \kappa }_{\text{l}}} - \hat n_{\text{g}}^2{{\hat \kappa }_{\text{g}}}} \right) $

(7)

式中：$ {\hat \delta _T} $为约化Tolman长度，$ {\hat \delta _{\text{T}}} = n_{\text{c}}^{{\text{1/3}}}{\delta _{\text{T}}} $；$ {\hat \sigma _\infty } $为约化表面张力，$ {\hat \sigma _\infty } = {{{\sigma _\infty }} \mathord{\left/ {\vphantom {{{\sigma _\infty }} {({k_{\text{B}}}{T_{\text{c}}}n_{\text{c}}^{1/3})}}} \right. } {({k_{\text{B}}}{T_{\text{c}}}n_{\text{c}}^{1/3})}} $；$ \hat n $为约化分子数密度，$ \hat n = {n \mathord{\left/ {\vphantom {n {{n_{\text{c}}}}}} \right. } {{n_{\text{c}}}}} $；κ和$ \hat \kappa $分别为等温压缩系数及其约化值，$ \hat \kappa = {k_{\text{B}}}{T_{\text{c}}}{n_{\text{c}}}\kappa $；下标l、g分别表示液相和气相。由式(7)结合完整标度理论，还可推得近临界表面张力的渐近表达式：

$ {\hat \sigma _\infty } = 0.274\;7{\left( {{{B_0^2} \mathord{\left/ {\vphantom {{B_0^2} {\mathit{\Gamma} _0^ - }}} \right. } {\mathit{\Gamma} _0^ - }}} \right)^{2/3}}{\left| {{\Delta}\hat T} \right|^{2\nu }} $

(8)

式中：v为临界指数，ν ≈0.6；临界振幅B₀和$ \mathit{\Gamma} _0^ - $由标度性质$ {\hat n_{\text{l}}} - {\hat n_{\text{g}}} \approx 2{B_0}{\left| {{\Delta}\hat T} \right|^\beta } $和$ \hat n_{\text{l}}^2{\hat \kappa _{\text{l}}} + \hat n_{\text{g}}^2{\hat \kappa _{\text{g}}} \approx 2\mathit{\Gamma} _0^ - {\left| {{\Delta}\hat T} \right|^{ - \gamma }} $取对数，线性回归得到；β、γ为临界指数。式(7)和(8)涉及的所有热力学性质均可以由自由能密度算得。

进一步可计算Tolman长度对成核功ΔW的修正，如式(9)所示；ΔW关系到成核速率J的大小^[19]，如式(10)所示，其中J₀是与过饱和度或过冷度有关的系数：

$ \Delta W = \frac{4}{3}\pi {\sigma _\infty }R_{\text{e}}^2\left( {1 - \frac{{4{\delta _{\text{T}}}}}{{{R_{\text{e}}}}}} \right) $

(9)

$ J = {J_0}\exp \;( - {{\Delta W} \mathord{\left/ {\vphantom {{\Delta W} {({k_{\text{B}}}T)}}} \right. } {({k_{\text{B}}}T)}}) $

(10)

2.3 状态方程的选取

对于普通的弱极性、非缔合流体，SAFT-VR方程已经足够准确，但是对于强极性、缔合流体，以水为例，SAFT-VR方程预测水的a₃为−0.030 9^[20]，但由两相平衡数据^[21]可以推算出其a₃应大于零，若用SAFT-VR+GRG计算，则水的Tolman长度在临界点将趋于正无穷，与完整标度理论预言(趋于负无穷)完全相反。而PC-SAFT方程预测的a₃与相平衡实验结果吻合，因此，本研究选用PC-SAFT方程，u_att(r)取为Lennard-Jones势的吸引部分，d为Lennard-Jones势表征分子大小的特征参数；ε为能量特征参数。

$ {u_{{\text{att}}}}(r) = \left\{ \begin{array}{l} 0, \;\;\;r < d \hfill \\ 4\varepsilon \left[ {{{(d/r)}^{12}} - {{(d/r)}^6}} \right], \;\;r \geqslant d \hfill \\ \end{array} \right. $

(11)

3 结果与讨论 3.1 模型参数

本研究考察的流体包括正构烷醇、醋酸和水。表 1给出计算所需的参数^[16]，包括2个和缔合有关的参数：缔合能ε_HB、有效缔合体积κ_HB。L₀^*为L₀的约化，L₀^*=L₀/d，a_{3, EOS}为状态方程预测的a₃，u₁为式(1)中的待定参数。流体的缔合类型都取2B模型^[22]。

临界指数和振幅^[12]、另一个场混合系数b₂的计算结果列于表 2。$ A_0^ - $为与热容有关的参数，b₂的计算结果与状态方程的预测值b_{2, EOS}差别不大，说明对哈密顿量的修正是自洽的，标度律α + 2β + γ = 2以及普适的振幅比^[23]$ {{\alpha A_0^ - \mathit{\Gamma} _0^ - } \mathord{\left/ {\vphantom {{\alpha A_0^ - \mathit{\Gamma} _0^ - } {B_0^2}}} \right. } {B_0^2}} = 0.022 $也得到满足。临界指数α的计算结果偏大，原因在于本研究使用的重正化方法来自Wilson的最初设想，和Wilson的原始结果^[24]以及单圈重正化结果^[25]接近。

3.2 表面张力和Tolman长度

图 1、2显示了各种缔合流体表面张力σ_∞的计算结果(实线)和较准确的经验公式^[26](数据点)的对比。严格说式(8)只有在临界点附近才准确，但由图可以看出，PC-SAFT+GRG计算的σ_∞在外推到远离临界点时依然和经验公式保持较高的一致。

图 3、4为PC-SAFT+GRG对从甲醇到己醇的正构烷醇、醋酸和水的Tolman长度的计算结果(实线)。由于缺乏Tolman长度的实验数据，本研究将取自NIST数据库^[21]的气液平衡数据以及由经验公式^[26]计算的σ_∞代入式(7)，估算其值(记作Eq+NIST，以数据点表示)，可以看出PC-SAFT+GRG与Eq+NIST符合良好。图中还用点虚线给出了水的DFT计算结果^[27]，可以看出DFT不能反映Tolman长度的临界发散行为，而改进的重正化群方法的结果与完整标度理论的预测一致。由于缺乏准确的近临界相平衡实验数据，图 4只画出了PC-SAFT+GRG计算结果。

将计算所得σ_∞和δ_T代入式(9)便能计算成核功。以R_e = 0.5 nm的微小水滴为例，图 5显示了本研究的计算结果(实线)和Eq+NIST(数据点)、DFT(点线)以及经典成核理论(即完全忽略Tolman修正，虚线)的对比，其中Eq+NIST、DFT、CNT用到的σ_∞均由经验公式^[26]计算。从图中可以看出考虑了Tolman修正的3种方法计算的成核功比较接近，它们都与经典成核理论的计算结果有显著差异，这将显著影响成核速率的预测。

4 结论

本研究将改进的重正化群理论与PC-SAFT方程结合，预测了包括醇、醋酸、水在内的各种缔合流体的表面张力和Tolman长度。PC-SAFT+GRG对缔合流体的表面张力与Tolman长度的描述不仅在临界点符合完整标度理论的预测，也能外推至非临界区，依然保持较高的准确性。经过Tolman修正的表面张力能对经典成核理论预测的成核速率作出改进。