1.   引　言

黏系数是流体力学中的重要概念, 也是化学化工领域的重要参数^[1]. 在Newtonian流体的Navier-Stokes描述中^[2,3], 应力张量与两个不同的输运系数密切相关, 即剪切黏度${\eta }_{{\mathrm{s}}}$和体积黏度${\eta }_{{\mathrm{v}}}$. 前者用于表征流体在体积不变的前提下发生剪切形变所引起的动量损失, 而后者则用于表征流体在形状不变的前提下发生体积形变而引起的动量损失. 长期以来, 人们对Newtonian流体中的剪切黏度进行了深入且细致的研究^[4–8]. 相比之下, 体积黏度却未受到广泛关注, 甚至在一些研究中被忽略不计, 即Stokes假设^[9]. 对单原子理想气体或不可压缩流体而言, 该假设并无不妥. 然而, 在诸如冲击波^[10]、超声衰减^[11]和流体振荡^[12]等过程中, 体积黏度起到不可忽视的作用. 正因如此, 流体的体积黏度正日益引起相关领域学者的研究兴趣.

超声吸附技术的发展为实验研究流体的体积黏度提供了可能的途径. 根据声学知识^[13], 体积黏度可由声波吸附系数算得. 然而, 计算中所涉及的等体和等压热容、零频声速、热传导和剪切黏度等物理量使得准确测定流体的体积黏度极为困难^[14,15]. 相比之下, 分子动力学模拟在预测体积黏度方面显得更具功效. 采用平衡分子动力学(equilibrium molecular dynamics, EMD)方法^[16–18], 体积黏度可以结合Green-Kubo关系或Einstein关系, 并通过计算应力张量的自相关函数在$t\to \infty$时的极限值来获得. 而在非平衡分子动力学(non-equilibrium molecular dynamics, NEMD)模拟^[19–21]中, 则需首先构建体系的非平衡稳态, 进而通过施加适当的外部激励来计算流体的响应. 上述两种分子动力学方法在体积黏度的理论研究中发挥着重要的作用. 除此之外, 结合流体力学和统计力学的基本原理, 人们也提出一些理论方法^[22,23]用于计算流体的体积黏度. 不同于EMD和NEMD, 此类方法仅需体系的结构信息(如径向分布函数)和粒子间两体势作为输入条件, 因而是体积黏度研究方法的有益补充.

大量研究已表明, 非均匀流体的结构和性质与均匀体相存在很大差异^[24–27]. 这种差异性主要源自于体系的非均匀结构和有限尺寸效应. 按其定义, 体积黏度是流体元体积对外部激励响应的量度. 事实上, 流体元的体积变化在很大程度上受制于其邻域内的流体结构. 据此, 不难理解非均匀流体的体积黏度也应与体系的局域结构密切相关, 并应呈非均匀分布. 然而, 目前已有的关于体积黏度的研究主要针对均匀体相流体, 而非均匀流体中体积黏度的相关研究却鲜有报道, 这也导致关于非均匀流体的体积黏度及其微观机理至今尚不明晰.

本文在Maxwell黏弹理论范畴内, 提出一种可用于计算非均匀流体体积黏度的理论方法. 该方法首先基于体相流体的体积黏度和力学模量, 并结合权重密度近似(weighted density approximation, WDA)定义了体系的局域弛豫时间. 进一步, 基于经典密度泛函理论(classical density functional theory, CDFT)所算得的平衡密度分布, 计算了体系中的力学模量分布. 最后, 结合Maxwell弛豫模型, 计算得到了体系中体积黏度的非均匀分布.

2.   均匀流体的黏弹弛豫模型

对于线性黏弹材料, 构造其本构关系通常有两种方法: 力学类比和Boltzmann叠加原理^[28]. 根据力学类比法, 线性黏弹性为可视为弹簧(spring)和活塞(dashpot)的线性组合. 二者分别描述体系中的弹性响应和黏性响应. 常用的力学类比模型主要有: Maxwell模型、Kelvin-Voigt模型、标准线性固体(standard linear solid)模型.

在图1中, 令并联的上、下两支中具有相同的体积应变, 并记为${\varepsilon}_{{\mathrm{v}}}$. 进一步假设上、下支中的压缩应力分别为${\varPi }_{{\mathrm{v}}1}$和${\varPi }_{{\mathrm{v}}2}$, 则总应力为${\varPi }_{{\mathrm{v}}}= {\varPi }_{{\mathrm{v}}1}+{\varPi }_{{\mathrm{v}}2}$. 从物理角度来看, 上支描述体系中的弹性成分, 即:

$${\varPi }_{{\mathrm{v}}1}={\kappa }_{0}{\varepsilon}_{{\mathrm{v}}},$$

(1)

式中${\kappa }_{0}$为体积模量的非弛豫成分. 事实上, 下支中的弹簧和活塞组合即为Maxwell模型的力学结构. 因二者串联, 故弹簧和活塞中的应力均为${\varPi }_{{\mathrm{v}}2}$. 与此同时, 下支中的应变率应为弹簧和活塞中的应变率之和, 即:

$$\frac{{\varPi }_{{\mathrm{v}}2}}{{\eta }_{1}}+\frac{1}{{\kappa }_{1}}\frac{{\mathrm{d}}{\varPi }_{{\mathrm{v}}2}}{{\mathrm{d}}t}=\frac{{\mathrm{d}}{\varepsilon}_{{\mathrm{v}}}}{{\mathrm{d}}t}={\dot{\varepsilon}}_{{\mathrm{v}}},$$

(2)

式中, 第一、二两项分别为黏性和弹性的贡献, ${\kappa }_{1}$为体积模量的弛豫成分, ${\eta }_{1}$为体积黏度系数. 事实上, (2)式可进一步变形为

$${\varPi }_{{\mathrm{v}}2}+\frac{{\eta }_{1}}{{\kappa }_{1}}\frac{{\mathrm{d}}{\varPi }_{{\mathrm{v}}2}}{{\mathrm{d}}t}={\eta }_{1}{\dot{\varepsilon}}_{{\mathrm{v}}}.$$

(3)

为方便起见, 定义算符: $\mathcal{A}=1+{\tau }_{{\mathrm{v}}}\dfrac{{\mathrm{d}}}{{\mathrm{d}}t}$, 其中${\tau }_{{\mathrm{v}}}\equiv {{\eta }_{1}}/{{\kappa }_{1}}$为Maxwell弛豫时间. 据此, (3)式可改写为 $\mathcal{A}{\varPi }_{2}= {\eta }_{1}{\dot{\varepsilon}}_{{\mathrm{v}}}$, 即:

$${\varPi }_{{\mathrm{v}}2}={\mathcal{A}}^{-1}{\eta }_{1}{\dot{\varepsilon}}_{{\mathrm{v}}}={\kappa }_{1}\left(1-{\mathcal{A}}^{-1}\right){\varepsilon}_{{\mathrm{v}}}.$$

(4)

(4)式计算中用到${\eta }_{1}\dfrac{{\mathrm{d}}}{{\mathrm{d}}t}={\kappa }_{1}\left(\mathcal{A}-1\right)$.

结合(1)式和(4)式, 体系中的总应力可以表达为

$${\varPi }_{{\mathrm{v}}}=\left[{\kappa }_{0}+{\kappa }_{1}\left(1-{\mathcal{A}}^{-1}\right)\right]{\varepsilon}_{{\mathrm{v}}}.$$

(5)

对简谐压缩应变${\varepsilon}_{{\mathrm{v}}}\sim{{\mathrm{e}}}^{{\mathrm{i}}\omega t}$而言, $\dfrac{{\mathrm{d}}}{{\mathrm{d}}t}\sim{\mathrm{i}}\omega$且$\mathcal{A}=1+ {\mathrm{i}}\omega {\tau }_{{\mathrm{v}}}$. 据此, 体系的复弹性${K}^{*}$和复黏度${\eta }^{*}$可分别表示为

$${K}^{*}\left(\omega \right)=\frac{{\varPi }_{{\mathrm{v}}}}{{\varepsilon}_{{\mathrm{v}}}}={\kappa }_{0}+{\kappa }_{1}\frac{{\mathrm{i}}\omega {\tau }_{{\mathrm{v}}}}{1+{\mathrm{i}}\omega {\tau }_{{\mathrm{v}}}},\tag{6a}$$

$${\eta }^{*}\left(\omega \right)=\frac{{\varPi }_{{\mathrm{v}}}}{{\dot{\varepsilon}}_{{\mathrm{v}}}}=\frac{{K}^{*}\left(\omega \right)}{{\mathrm{i}}\omega }.\tag{6b}$$

按定义, ${K}^{*}\left(\omega \right)$和${\eta }^{*}\left(\omega \right)$的实部应分别为体系的弹性系数${K}_{{\mathrm{v}}}\left(\omega \right)$和黏度系数${\eta }_{{\mathrm{v}}}\left(\omega \right)$, 即:

$${K}_{{\mathrm{v}}}\left(\omega \right)={\kappa }_{0}+{\kappa }_{1}\frac{{\left(\omega {\tau }_{{\mathrm{v}}}\right)}^{2}}{1+{\left(\omega {\tau }_{{\mathrm{v}}}\right)}^{2}},\tag{7a}$$

$${\eta }_{{\mathrm{v}}}\left(\omega \right)={\kappa }_{1}\frac{{\tau }_{{\mathrm{v}}}}{1+{\left(\omega {\tau }_{{\mathrm{v}}}\right)}^{2}}. \tag{7b}$$

事实上, (7a)式和(7b)式中的${\kappa }_{0}$和${\kappa }_{1}$可由其在频率空间中的边界条件得出. 当$\omega \to 0$时, ${K}_{{\mathrm{v}}}\left(0\right) = {K}_{0}$; 当$\omega \to \infty$时, ${K}_{{\mathrm{v}}}\left(\infty \right) = {K}_{\infty }$. 其中, ${K}_{0}$和${K}_{\infty }$分别为体系的零频和高频体积模量. 据此可得: ${\kappa }_{0}={K}_{0}$, ${\kappa }_{1}={K}_{2}\equiv {K}_{\infty }-{K}_{0}$. ${K}_{2}$为体积模量的弛豫贡献. 据此, 在零频极限下, 体系的体积黏度系数${\eta }_{{\mathrm{v}}}$可表示为

$${\eta }_{{\mathrm{v}}}={K}_{2}{\tau }_{{\mathrm{v}}}.$$

(8)

需指出的是, (8)式将作为本文研究非均匀流体体积黏度的出发点. 然而, 此式描述的是均匀体相流体. 因此, 下面需将其推广至非均匀流体的情况.

3.   非均匀黏度的Maxwell弛豫模型

在非均匀流体的研究中, 相应的体相性质常被用于参考甚至被用作相关计算的基本出发点. 在此类计算中, 通常的作法是直接借鉴DFT中对自由能密度泛函的处理方式, 即: 由某种局域的平均密度或权重密度来替代体相结果中的体密度, 进而得到非均匀流体的相关性质. 其中最具代表性的是由Bitsanis等^[29]提出的局域平均密度模型 (local average density model, LADM). 事实上, LADM及其改进版本^[30]已被广泛用于研究受限流体的动态性质. 而且, 研究中所采用的权重函数和近似方法也广泛涉及van der Waals模型、 Fischer-Methfessel模型、硬棒模型以及Tarazona模型. 然而, 其预测结果与模拟结果的比较表明^[31], LADM的可靠性在很大程度上依赖于计算中所选取的权重函数类型. 由此可见, 直接对体相的动态性质进行权重密度近似并不能很好地预测非均匀体系的动态性质. 为了克服这一缺陷, 本文以(8)式为出发点, 提出一种可用于计算非均匀流体体积黏度的新颖方法.

若已知体相流体的体积黏度${\eta }_{{\mathrm{v}}}^{{\mathrm{b}}}$和弛豫模量${K}_{2}^{{\mathrm{b}}}$, 结合(8)式, 可定义体系的局域弛豫时间${\tau }_{{\mathrm{v}}}\left(z\right)$:

$${\tau }_{{\mathrm{v}}}\left(z\right)=\frac{{\eta }_{{\mathrm{v}}}^{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)}{\left[{K}_{\infty }^{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)-{K}_{0}^{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)\right]},$$

(9)

式中, 上角标b表示体相. 此外, $\bar{\rho }\left(z\right)$为权重密度:

$$\bar{\rho }\left(z\right)=\int \rho \left({z}'\right)\omega \left(\left|{z}'-z\right|\right){\mathrm{d}}{z}',$$

(10)

其中$\omega \left(z\right)$为权重函数, $\rho \left(z\right)$为密度分布. 众所周知, 对于某一选定的流体粒子, 其弛豫过程主要决定于其邻域内的流体结构及其所处区域内的能垒高低. 相比之下, 距离愈远, 流体的结构对选定粒子的弛豫过程的影响愈弱. 有鉴于此, 可借助高斯型权重函数来描述这一正态分布特征:

$$\omega \left(z\right)=\frac{1}{\sqrt{2\pi }\xi }{{\mathrm{e}}}^{-{{z}^{2}}/({2{\xi }^{2}})},$$

(11)

其中$\xi$是用于调控权重函数分布的参数.

对非均匀流体而言, 结构上的非均匀性必将导致其性质的非均匀性, 包括弛豫模量${K}_{2}$. 为此, 可进一步定义体系的局域弛豫模量:

$${K}_{2}\left(z\right)={K}_{\infty }\left(z\right)-{K}_{0}\left(z\right),$$

(12)

其中${K}_{\infty }\left(z\right)$和${K}_{0}\left(z\right)$非均匀流体的高频和零频体积模量, 二者的具体计算将在下节给出. 结合(8)式、(9)式和(12)式, 非均匀流体的体积黏度可最终表示为

$${\eta }_{{\mathrm{v}}}\left(z\right)={K}_{2}\left(z\right)\frac{{\eta }_{{\mathrm{v}}}^{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)}{\left[{K}_{\infty }^{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)-{K}_{0}^{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)\right]}.$$

(13)

不难发现, (13)式中${K}_{2}(z)$和${K}_{\infty }^{{\mathrm{b}}}(\bar{\rho }(z)) - {K}_{0}^{{\mathrm{b}}}(\bar{\rho }(z))$均是对体系中的局域弛豫模量的描述. 其中, 后者是将体相结果中的体密度替代为LADM权重密度$\bar{\rho }(z)$而得. 然而, 这一简单的替代将导致粒子间关联在一定程度上被低估. 相比之下, ${K}_{2}(z)$则是直接由粒子分布函数和两体势函数计算所得, 故其可更多地计入粒子间关联对弛豫模量的贡献. 为了说明该问题, 图2给出了狭缝中Lennard-Jones (LJ)流体的局域弛豫模量的空间分布. 图中结果表明, 在所选定的计算条件下, 狭缝中各点的${K}_{\infty }^{{\mathrm{b}}}(\bar{\rho }(z))- {K}_{0}^{{\mathrm{b}}}(\bar{\rho }(z))$均小于${K}_{2}(z)$的值. 当${\rho }_{{\mathrm{b}}}^{*}=0.8$时, 粒子间较强的短程关联还可使${K}_{2}(z)$曲线呈现出显著的振荡趋势. 然而, 上述振荡趋势在${K}_{\infty }^{{\mathrm{b}}}(\bar{\rho }(z)) -{K}_{0}^{{\mathrm{b}}}(\bar{\rho }(z))$曲线中却未能呈现. 以上结果均说明权重密度$\bar{\rho }(z)$的引入在一定程度上削弱了短程关联对弛豫模量的贡献.

(13)式表明, 体相流体的体积黏度是计算非均匀流体体积黏度的先决条件. 为便于后续的计算, 本文采用Heyes基于NEMD所提出的拟合结果^[32]:

$${\eta }_{{\mathrm{v}}}^{{\mathrm{b}}*}=\sum _{i=0}^{3}{\rho }^{*i}\left({c}_{i1}+{c}_{i2}{T}^{*{c}_{T}}\right),$$

(14)

其中 ${c}_{T}=0.05$, 其他拟合系数${c}_{i1}$和${c}_{i2}$在表1中给出. 需指出的是, (14)式中各量均采用无量纲形式: ${\eta }^{*}=\eta {\sigma }^{2}/\sqrt{m\varepsilon }$, ${\rho }^{*}=\rho {\sigma }^{3}$, ${T}^{*}={k}_{{\mathrm{B}}}T/\varepsilon$. $m$, $\sigma$和$\varepsilon$分别为流体粒子的质量、直径和作用强度, ${k}_{{\mathrm{B}}}$为Boltzmann常数, $T$为绝对温度.

事实上, 体积黏度的非均匀分布也可通过将(14)式中的体密度$\rho$直接替换为权重密度$\bar{\rho }\left(z\right)$来得到, 即LADM方案. 相比于此方案, (13)式的优势体现于如下两方面. 一方面, (13)式中的分子项和分母项均会因权重密度$\bar{\rho }\left(z\right)$的引入而导致其短程关联贡献被部分削弱. 然而, 由于体积黏度和弛豫模量对流体密度均具有正相关的依赖关系, 这使得$\bar{\rho }\left(z\right)$对二者比值(即局域弛豫时间${\tau }_{{\mathrm{v}}}\left(z\right)$)的影响会小于其对体积黏度${\eta }_{{\mathrm{v}}}^{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)$的直接影响. 另一方面, 正如图2所示, ${K}_{2}\left(z\right)$在描述局域弛豫模量中的关联贡献时比LADM方案更具优势, 这使得(13)式可以更多地计入关联效应对体积黏度的贡献.

4.   非均匀流体中的力学模量

4.1   非均匀流体的高频体积模量

对各向同性的均匀流体而言, 体系的高频力学模量可由Zwanzig-Mountain公式^[33]来算得. 然而, 该公式并不适于非均匀流体的相关计算. 事实上, 在前期研究^[34]中, 本文作者基于胡克定律和应力张量的Irving-Kirkwood公式已计算得到可用于描述非均匀流体的高频力学模量, 且该结果在体相极限下可回归至Zwanzig-Mountain公式. 根据文献[34]结果, 高频体积模量可表示为

$${K}_{{\mathrm{\infty }}}\left(z\right)={K}_{{\mathrm{\infty }}}^{\left({\mathrm{k}}\right)}\left(z\right)+{K}_{{\mathrm{\infty }}}^{\left({\mathrm{p}}\right)}\left(z\right),$$

(15)

式中, 动能贡献${K}_{{\mathrm{\infty }}}^{\left({\mathrm{k}}\right)}\left(z\right)={5}/{3}{k}_{{\mathrm{B}}}T\rho \left(z\right)$精确可知, 而势能贡献${K}_{{\mathrm{\infty }}}^{\left({\mathrm{p}}\right)}\left(z\right)$可表示为

$${K}_{{\mathrm{\infty }}}^{\left({\mathrm{p}}\right)}\left(z\right)=-\frac{2}{9}{I}_{1}\left(z\right)+\frac{1}{9}{I}_{2}\left(z\right),$$

(16)

其中辅助函数${I}_{m}(m={\mathrm{1, 2}})$定义式为

$$\begin{split}{I}_{m}\left(z\right)=\;&\frac{1}{2}\int {\mathrm{d}}\boldsymbol{r}{r}^{m}{\varphi }^{\left(m\right)}\left(r\right)\\ &\times \int_{0}^{1}{\mathrm{d}}\zeta {\rho }^{\left(2\right)}\left(z-\zeta \boldsymbol{r},z-\zeta \boldsymbol{r}+\boldsymbol{r}\right). \end{split}$$

(17)

这里$\varphi \left(r\right)$为流体粒子两体势, ${\varphi }^{\left(m\right)}\left(r\right)$为$\varphi \left(r\right)$的$m$阶导数, ${\rho }^{\left(2\right)}\left({\boldsymbol{r}}_{1}, {\boldsymbol{r}}_{2}\right) = \rho \left({\boldsymbol{r}}_{1}\right)\rho \left({\boldsymbol{r}}_{2}\right)g\left({r}_{12}\right)$为两体密度函数, $g\left({r}_{12}\right)$为径向分布函数, 且${r}_{12}=\left|{\boldsymbol{r}}_{2}-{\boldsymbol{r}}_{1}\right|$. 对于体相LJ流体而言, (15)式可由其压强$P$和内能$E$表示为

$${K}_{\infty }^{{\mathrm{b}}}=\frac{20}{3}{k}_{{\mathrm{B}}}T{\rho }_{{\mathrm{b}}}+7P-8\rho E,$$

(18)

其中${\rho }_{{\mathrm{b}}}$为体密度; $P$和$E$可由MBWR状态方程给出^[35].

由(15)式不难发现, ${K}_{{\mathrm{\infty }}}\left(z\right)$的计算需将流体分子的密度分布和结构信息作为输入条件. 这些均可由CDFT来得出. 在CDFT的数值计算中, 本文分别采用基本度量理论^[36]和权重密度近似^[37]来处理硬球作用和长程吸引对Helmholtz自由能的贡献, 并由Picard迭代的方法来计算体系的平衡密度分布. 关于CDFT和MBWR的计算细节可见补充材料.

4.2   非均匀流体的零频体积模量

根据经典热力学理论, 均匀体相的等温体积模量${K}_{T}^{{\mathrm{b}}}$可表示为

$${K}_{T}^{{\mathrm{b}}}={\rho }_{{\mathrm{b}}}^{2}{\left(\frac{\partial {\rho }_{{\mathrm{b}}}}{\partial \mu }\right)}_{T}^{-1},$$

(19)

其中$\mu$为体系化学势. 采用统计力学求和规则的压缩方案, 可将上述结果推广至非均匀体系^[38]:

$${K}_{T}\left(z\right)={\left[\rho \left(z\right)\right]}^{2}{\left[\frac{\partial \rho \left(z\right)}{\partial \mu }\right]}_{T}^{-1}.$$

(20)

然而, 前文中的${K}_{0}^{{\mathrm{b}}}$和${K}_{0}\left(z\right)$均属绝热体积模量. 事实上, 在零频极限下, 绝热模量和等温模量可通过比热比$\gamma \equiv {C}_{{\mathrm{P}}}/{C}_{{\mathrm{V}}}$建立联系, 即:

$${K}_{0}^{{\mathrm{b}}}={\gamma }_{{\mathrm{b}}}{K}_{T}^{{\mathrm{b}}},\tag{21a}$$

$${K}_{0}\left(z\right)=\gamma \left(z\right){K}_{T}\left(z\right),\tag{21b}$$

式中, ${\gamma }_{{\mathrm{b}}}$和$\gamma \left(z\right)={\gamma }_{{\mathrm{b}}}\left(\bar{\rho }\left(z\right)\right)$分别为体相比热比和非均匀体系中的局域比热比. 就本文所研究的LJ流体而言, 二者均可借助MBWR状态方程予以计算.

5.   结果与讨论

本文旨在揭示非均匀流体中体积黏度的微观机制, 以期加深对体积黏度的认知. 为此, 本节将以受限于平行狭缝中的LJ流体为对象, 采用(13)式计算其体积黏度的空间分布及其相关的调制效应. 首先, 为了验证理论方法的可靠性, 将计算结果与模拟结果进行对比是很有必要的. 然而, 据我们所知, 关于非均匀流体体积黏度的模拟数据目前尚无报道. 因此, 首先采用同种方法对非均匀LJ流体的剪切黏度${\eta }_{{\mathrm{s}}}\left(z\right)$进行了计算, 并与相应的NEMD结果^[30]进行对比, 有关${\eta }_{{\mathrm{s}}}\left(z\right)$的计算方法可见文献[39]. 图3中的计算中已采用与NEMD相同的参数设置, 其中${z}^{*}=z/\sigma$, ${H}^{*}=H/\sigma$为约化缝宽. 图3结果表明, 此类方法对${\eta }_{{\mathrm{s}}}\left(z\right)$的预测结果与NEMD模拟结果吻合得很好.

基于上述理论方法, 还计算了狭缝中LJ流体的平均黏度. 其定义如下:

$${\bar{\eta }}_{i}=\frac{1}{H}\int _{0}^{H}{\eta }_{i}\left(z\right){\mathrm{d}}z\;\;\left(i={\mathrm{v}},{\mathrm{s}}\right).$$

(22)

为了评价(13)式的有效性, 着重计算了狭缝中LJ流体的平均体积黏度${\bar{\eta }}_{{\mathrm{v}}}^{*}$随缝宽${H}^{*}$的变化, 并在图4中将其与基于Green-Kubo (GK)公式的解析结果^[40]进行比较. 对比结果显示, 当${H}^{*} > 10.0$时, 本文结果与GK结果吻合的很好. 此外, 图4结果还表明, 当${H}^{*} < 5.0$时, 粒子间的短程关联还可导致${\bar{\eta }}_{{\mathrm{v}}}^{*}$随缝宽${H}^{*}$的振荡变化. 这与Jaeger 等^[41]对纳米通道中水的体积黏度随缝宽的变化趋势是一致的. 为了进一步评价本文的理论方法, 图4还给出了LADM的计算结果. 对比结果表明, 当${H}^{*}$较小时, LADM会低估体系的平均体积黏度.

在后续的计算中, 将基于(13)式来具体研究受限于缝宽为$H$的平行狭缝中LJ流体的体积黏度. 计算中, 狭缝单壁所提供的外势均由Steele 10-4-3势能函数来描述:

$$\begin{split} {\varphi }_{{\mathrm{w}}{\mathrm{a}}{\mathrm{l}}{\mathrm{l}}}\left(z\right)=\;&2{\mathrm{\pi }}{\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}{\rho }_{{\mathrm{s}}}{\sigma }_{{\mathrm{f}}{\mathrm{s}}}^{2}\bigg[\frac{2}{5}{\left(\frac{{\sigma }_{{\mathrm{f}}{\mathrm{s}}}}{z}\right)}^{10}-{\left(\frac{{\sigma }_{{\mathrm{f}}{\mathrm{s}}}}{z}\right)}^{4} \\ & -\frac{{\sigma }_{{\mathrm{f}}{\mathrm{s}}}^{4}}{3\varDelta {\left(z+0.61\varDelta \right)}^{3}}\bigg],\end{split}$$

(23)

其中, ${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}$和${\sigma }_{{\mathrm{f}}{\mathrm{s}}}$分别为固-液作用的强度和直径, ${\rho }_{{\mathrm{s}}}$和$\varDelta$分别为固体原子面密度和层间距. 此外, ${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}$和${\sigma }_{{\mathrm{f}}{\mathrm{s}}}$均满足Lorentz混合规则. 据此, 狭缝中的外势可表达为

$${V}_{{\mathrm{e}}{\mathrm{x}}{\mathrm{t}}}\left(z\right)={\phi }_{{\mathrm{w}}{\mathrm{a}}{\mathrm{l}}{\mathrm{l}}}\left(z\right)+{\phi }_{{\mathrm{w}}{\mathrm{a}}{\mathrm{l}}{\mathrm{l}}}\left(H-z\right).$$

(24)

若无特别说明, 下文计算中均采用如下参数: ${\sigma }_{{\mathrm{f}}{\mathrm{s}}}/\sigma =\varDelta /\sigma =1$, ${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon =0.483$, ${\rho }_{{\mathrm{s}}}{\sigma }_{{\mathrm{f}}{\mathrm{s}}}^{2}=4.429$, $\xi /\sigma = 0.8$.

5.1   体密度对体积黏度的影响

对受限流体而言, 体密度是影响体系中关联效应的重要因素. 众所周知, 粒子间的短程排斥和长程吸引分别对体系的结构和性质产生决定性的影响. 因此, 掌握体积黏度与体密度之间的关联是理解其微观机制的重要内容. 为此, 将在${H}^{*}=6.0$的条件下研究体密度对受限流体体积黏度的影响. 通常, 纳米孔隙中的低温气体在其压强低于饱和气压时可引发毛细凝聚. 为了考察这一结构相变对体积黏度的影响, 计算将分别在${T}^{*}=2.0$和${T}^{*}=1.0$的情况下进行, 相应的计算结果在图5中给出.

图5(a)给出了${T}^{*}=2.0$时${\eta }_{{\mathrm{v}}}^{*}\left({z}^{*}\right)$随体密度${\rho }_{{\mathrm{b}}}^{*}$的变化. 图中结果表明, 当${\rho }_{{\mathrm{b}}}^{*}$较小时, ${\eta }_{{\mathrm{v}}}^{*}$在两壁附近呈现明显的峰值, 但在狭缝的中部区域却呈微弱的气体黏度. 随着${\rho }_{{\mathrm{b}}}^{*}$的变大, 狭缝中各处的${\eta }_{{\mathrm{v}}}^{*}$均呈单调增强的变化趋势, 且在中部区域出现振荡结构. 事实上, 该变化是由受限流体的结构变化和体积黏度与局域密度间的正相关性所导致的, 即${\eta }_{{\mathrm{v}}}^{*}\left({z}^{*}\right)$曲线的振荡源自${\rho }_{{\mathrm{b}}}^{*}$增大所引发的粒子层状排列. 从微观角度来看, 在密度较大的区域, 粒子间较强的排斥体积使得流体发生体积形变更为困难, 因而体积黏度更大.

图5(b)中的结果与图5(a)存在显著的区别. 图5(b)结果显示, 随着${\rho }_{{\mathrm{b}}}^{*}$的变大, 狭缝中各处的${\eta }_{{\mathrm{v}}}^{*}$均呈先降后增的变化趋势. 具体来讲, 当${\rho }_{{\mathrm{b}}}^{*}$较小时, 狭缝中部区域的${\eta }_{{\mathrm{v}}}^{*}$具有典型的液体黏度, 并随着空间位置呈振荡分布. 然而, 随着${\rho }_{{\mathrm{b}}}^{*}$的变大, ${\eta }_{{\mathrm{v}}}^{*}$的振荡减弱且幅值也随之变小. 这一结果可归因于${\rho }_{{\mathrm{b}}}^{*}$对狭缝中毛细凝聚的调制. 通常, 凝聚液相具有远大于其体相的密度和体积黏度. 然而, 随着压强逐渐接近饱和气压, 狭缝中的毛细凝聚将被削弱. 图5(b)结果还表明, 当${\rho }_{{\mathrm{b}}}^{*} > 0.6$时, ${\rho }_{{\mathrm{b}}}^{*}$的进一步增大将导致${\eta }_{{\mathrm{v}}}^{*}$的增强. 这与图5(a)变化趋势是一致的.

对比图3和图5结果不难发现, 受限流体的体积黏度和剪切黏度大小可比. 在均匀流体中, 通常可由二者的比值${\eta }_{{\mathrm{v}}}/{\eta }_{{\mathrm{s}}}$来量度体系发生体积形变和剪切形变的相对难易程度. 事实上, Cowan与Leech^[42]已通过超声衰减实验研究了液态Ar, Kr和Xe的${\eta }_{{\mathrm{v}}}/{\eta }_{{\mathrm{s}}}$. 研究结果表明, 在${\rho }_{{\mathrm{b}}}^{*}\in \left(0.6, 0.85\right)$范围内体相液体的${\eta }_{{\mathrm{v}}}/{\eta }_{{\mathrm{s}}}$随密度的增加而减小, 但随温度的升高而增大. 由此可见, 密度的增大既不利于体积形变也不利于剪切形变, 但随之增强的排斥体积对剪切形变的抑制更为显著. 相反, 温度的升高使粒子的平均动能得以增大, 这既有利于体积形变也有利于剪切形变, 但其对剪切形变的增强效果更为显著.

不难理解, 受限流体发生体积形变和剪切形变的相对难易程度将会受到狭缝的影响与调制. 然而, 狭缝的限制效应和流体的相态变化对这一特征的调制机制仍不明晰. 为此, 定义${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}\equiv {\bar{\eta }}_{{\mathrm{v}}}/{\bar{\eta }}_{{\mathrm{s}}}$. 在此基础上, 将在${\rho }_{{\mathrm{b}}}^{*}\in \left({\mathrm{0, 1}}\right)$的范围内分别计算${T}^{*}= 2.0$和${T}^{*}=1.0$情况下受限流体的${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$的数值, 计算结果在图6中给出.

计算结果显示, 在选定的${\rho }_{{\mathrm{b}}}^{*}$范围内均有${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}} < 1$. 这说明预引起相同大小的体积应变率${\dot{\varepsilon}}_{{\mathrm{v}}}$和剪切应变率${\dot{\varepsilon}}_{{\mathrm{s}}}$, 所需压缩应力${{\mathrm{\varPi }}}_{{\mathrm{v}}}$应小于剪切应力${{\mathrm{\varPi }}}_{{\mathrm{s}}}$. 此外, 图6中${T}^{*}=2.0$的结果还表明, 在${\rho }_{{\mathrm{b}}}^{*} > 0.6$的高密度(液态)区内, ${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$表现出与Cowan实验结果^[42]相同的密度依赖关系, 即随${\rho }_{{\mathrm{b}}}^{*}$的增大而减小. 然而, 当${\rho }_{{\mathrm{b}}}^{*} < 0.6$时, 情况则完全相反. 这说明在低密度(气态)区内密度的增大对体积形变的抑制效果更为显著. 对比图6中两组数据还可发现, ${T}^{*}=1.0$时${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$对${\rho }_{{\mathrm{b}}}^{*}$的依赖关系与${T}^{*}=2.0$时几乎相同. 但是, 在${\rho }_{{\mathrm{b}}}^{*} < 0.1$低密度区内, ${T}^{*}=1.0$时${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$的数值明显大于${T}^{*}=2.0$时的对应数值. 其主要原因在于狭缝中所引发的毛细凝聚以及随之增强的排斥体积.

5.2   温度对体积黏度的影响

除体密度外, 温度也是影响体系结构与性质的重要参数, 正如图5和图6结果所示. 为了更全面地理解温度对体积黏度的影响, 本节将在${H}^{*}=6.0$的条件下, 分别对${\rho }_{{\mathrm{b}}}^{*}=0.01$和${\rho }_{{\mathrm{b}}}^{*}=0.6$的受限流体进行计算, 计算结果在图7中给出. 需要说明的是, 选取以上两个不同的体密度的目的在于考察在受限空间内, 毛细凝聚的发生会以何种方式影响并调制体系的体积黏度.

图7(a)给出了${\rho }_{{\mathrm{b}}}^{*}=0.01$时${\eta }_{{\mathrm{v}}}^{*}\left({z}^{*}\right)$的计算结果. 在此种情况下, 气体压强低于饱和气压, 故体系可在低温下引发毛细凝聚. 正如图7(a)中所显示, 随着${T}^{*}$的减小, 狭缝中各处的${\eta }_{{\mathrm{v}}}^{*}$均呈单调增加的变化趋势. 当${T}^{*} < 1.2$时, 毛细凝聚的发生使得受限流体呈现典型的液态黏度, 且在狭缝的中部区域出现振荡. 随着${T}^{*}$的进一步减小, 振荡趋势将变得愈发显著. 事实上, 以上结果均是由温度对流体相态的调制效应所导致的. 通常, 温度的升高可增强粒子的平均动能. 这不利于凝聚相的形成, 并最终削弱体系中的体积黏度. 与图7(a)中不同, 图7(b)中的体系始终保持为液相. 换言之, 此情况下${T}^{*}$的变化并未引发体系的结构相变. 然而, 温度的升高依然能够增加粒子的平均动能. 这有助于体系发生体积形变, 进而降低其体积黏度. 这一点与图7(a)中${\eta }_{{\mathrm{v}}}^{*}$随${T}^{*}$的变化趋势是一致的.

图8分别在${\rho }_{{\mathrm{b}}}^{*}=0.01$和${\rho }_{{\mathrm{b}}}^{*}=0.6$的情况下计算了 ${T}^{*}$对狭缝中${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$的影响. 图8中两组结果均表明, 不同温区内${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$随${T}^{*}$可呈现截然不同的变化趋势. 在${\rho }_{{\mathrm{b}}}^{*}=0.01$的曲线中, 低温区内${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$随${T}^{*}$而增大, 该趋势与Cowan的实验结果^[42]一致. 事实上, 这一变化过程与狭缝中的毛细凝聚有关, 即温度的升高使粒子的层状排列变得疏散, 进而降低了凝聚液态的平均密度, 最终使得${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$变大. 与此相反, 在高温区${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$却呈现随${T}^{*}$的变大而减小的趋势. 这是由于粒子平均动能的增加将有益于体积形变的发生. 此外, 对比图8中两条曲线还可发现, ${\rho }_{{\mathrm{b}}}^{*}=0.01$曲线的变化趋势比${\rho }_{{\mathrm{b}}}^{*}=0.6$曲线更为剧烈. 实际上, 在低温区, 前者中凝聚液态的平均密度比后者的平均密度更大. 而且, 相比于后者中${T}^{*}$对受限液体的调制, 前者中${T}^{*}$对毛细凝聚的调制具有更高的灵敏度. 然而在高温区, 前者中气态的平均密度比后者的平均密度更小.

5.3   缝宽对体积黏度的影响

就受限流体而言, 孔隙的限制效应是导致其结构与性质区别于体相的直接原因. 其中, 孔隙尺寸和吸附势强度是两个关键因素, 二者可通过对体系结构的调制来实现对体系性质的间接影响. 为了深入理解非均匀流体的体积黏度, 将分别研究缝宽和吸附势强度对体积黏度的影响, 以期获得更多关于非均匀体积黏度的微观信息.

图9在${T}^{*}=1.5$和${\rho }_{{\mathrm{b}}}^{*}=0.8$的条件下计算了${\eta }_{{\mathrm{v}}}^{*}\left({z}^{*}\right)$和${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$随缝宽${H}^{*}$的变化. 图9(a)结果显示, 狭缝中${\eta }_{{\mathrm{v}}}^{*}$的空间分布模式与${H}^{*}$存在密切的关联. 意料之中的是, 当${H}^{*}$取不同值时, 两壁的吸附势总可在其附近区域引起粒子的聚集, 进而使该区域呈现${\eta }_{{\mathrm{v}}}^{*}$的峰值. 然而有趣的是, 缝宽${H}^{*}$每增加一个粒子直径均会导致${\eta }_{{\mathrm{v}}}^{*}$曲线新增一个峰值, 直至缝宽增至${H}^{*}\sim7$. 事实上, 这可归因于缝宽变化对狭缝中流体层状结构的调制. 当缝宽恰好可以容纳整数层粒子时, 层间较强的排斥体积不利于体系发生体积形变, 因此层状结构的形成有益于体积黏度的增强.

不仅如此, 图9(b)结果表明缝宽${H}^{*}$对${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$也具有类似的调制效果. 一方面, 当${H}^{*}$不是很大时, ${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}\sim{H}^{*}$曲线呈现明显的振荡, 且周期约为一个粒子直径. 这是受限流体的整体性质(如剩余吸附、溶剂化力等)所具有的普遍特征, 同时也说明了粒子间的排斥体积所起到的关键作用. 另一方面, 随着${H}^{*}$的增大, ${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$将逐渐减小至某一渐近值. 此外, 由图9(b)结果不难发现, 受限流体的${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$始终大于其体相值. 事实上, 狭缝对流体粒子的吸附作用使体系的平均密度得以增加, 这对于体系发生体积形变和剪切形变均是不利的. 然而, 狭缝中沿两壁法向的限制效应导致流体发生体积形变要比剪切形变更加困难.

5.4   吸附强度对体积黏度的影响

狭缝两壁的吸附强度也是影响流体结构与性质的另一重要外因. 为了进一步理解狭缝的自身特征对其中所限流体体积黏度的影响, 本文将在${T}^{*}=1.5$和${H}^{*}=6.0$的条件下计算吸附强度对${\eta }_{{\mathrm{v}}}^{*}\left({z}^{*}\right)$和${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$的影响. 考虑到狭缝中可能引发的毛细凝聚及其对体积黏度的调制, 计算中已选取${\rho }_{{\mathrm{b}}}^{*}=0.1$和${\rho }_{{\mathrm{b}}}^{*}=0.6$两个不同的体系. 相关的计算结果可见图10和图11.

图10(a)给出了${\rho }_{{\mathrm{b}}}^{*}=0.1$时吸附强度对${\eta }_{{\mathrm{v}}}^{*}$的影响. 图10(a)结果显示, 当${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$较小时, 除在两壁附近${\eta }_{{\mathrm{v}}}^{*}$取极大值外, 其余区域的${\eta }_{{\mathrm{v}}}^{*}$的分布几乎均匀, 且具有气体黏度的典型数值. 然而, 随着${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$的增强大, 狭缝中各处的${\eta }_{{\mathrm{v}}}^{*}$均呈单调增加趋势. 尤其是在狭缝的中部区域, 在${\eta }_{{\mathrm{v}}}^{*}$随${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$增强而增大的同时, 曲线还呈现出新的振荡结构. 这是由增强的吸附势所引起的毛细凝聚所导致的. 由此可见, 吸附势的增强可提高发生毛细凝聚的临界压强(密度), 因为当${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon =0.483$时, 具有${T}^{*}=1.5$和${\rho }_{{\mathrm{b}}}^{*}=0.1$的气相在${H}^{*}=6.0$的狭缝中是不能引发毛细凝聚的. 图10(b)结果显示, 当${\rho }_{{\mathrm{b}}}^{*}=0.6$时, ${\eta }_{{\mathrm{v}}}^{*}$随${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$呈单调增加趋势. 不难理解, 此时较强的吸附势有利于体积黏度的增强.

图11分别在${\rho }_{{\mathrm{b}}}^{*}=0.1$和${\rho }_{{\mathrm{b}}}^{*}=0.6$的条件下计算了吸附强度对${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$的影响. 图11结果显示, 当${\rho }_{{\mathrm{b}}}^{*}=0.1$时, ${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$随着${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$的增大呈先增后降的变化趋势. 事实上, 当${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon \ll 1$时, 狭缝的吸附能力较为微弱. 此时, 两壁附近的粒子因受排空作用而脱离狭缝的限制, 从而导致狭缝中的平均密度低于体相密度. 因此, ${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$越小, 受限流体越容易发生体积形变. 相反, 当${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$足够大时, 狭缝具有较强的吸附能力, 并可导致受限流体的平均密度高于体相密度. 与此同时, ${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$的进一步增大将增强狭缝中的毛细凝聚. 以上两个因素协同导致${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$的减小. 当${\rho }_{{\mathrm{b}}}^{*}=0.6$时, ${R}_{\eta }^{{\mathrm{p}}{\mathrm{o}}{\mathrm{r}}{\mathrm{e}}}$随着${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$的增大呈单调减小的变化趋势. 这是由于该情况下受限流体始终为液态, 且其平均密度随${\varepsilon }_{{\mathrm{f}}{\mathrm{s}}}/\varepsilon$的增大而逐渐变大.

6.   结　论

基于Maxwell黏弹弛豫理论, 本文提出一种计算非均匀流体中体积黏度的新方法. 该方法首先基于体相流体的体积黏度和弛豫模量计算体系的局域弛豫时间, 进一步结合非均匀流体中的零频和高频力学模量得出非均匀流体的体积黏度. 采用此方法, 本文计算了受限于平行狭缝中的LJ流体的体积黏度, 系统地研究了体密度、温度、缝宽和吸附强度等因素对其空间分布的调制规律, 并从微观层面阐释了其物理机制. 研究表明, 非均匀流体的体积黏度与平衡结构存在密切关联, 且体系的结构相变可显著调制体积黏度的大小. 此外, 研究还表明, 体密度和吸附强度的增大均可增强体系的体积黏度, 而温度的升高则会削弱非均匀流体的体积黏度.

需说明的是, 本文着重研究了狭缝中LJ流体的体积黏度. 因此, 所涉及的弛豫模量和弛豫时间均是针对LJ流体来计算的. 事实上, 该方法同样可应用于其他简单的模型流体(如Yukawa流体、电解液等)及其混合物. 此外, 由于该方法基于普遍的Maxwell弛豫理论, 故将其推广至具有一定分子结构的分子流体也是可行的. 譬如在高分子流体中, 其体相弛豫时间与密度的依赖关系可借助理论方法^[13,43]得出, 而局域弛豫模量可结合体系的压力张量和线性弹性理论予以计算. 当然, 与原子流体相比, 分子流体中复杂的分子结构和势能函数使该计算变得更为复杂. 从理论角度来看, 上述尝试必将有益于加深对受限复杂流体中体积黏度的理解, 并可为研究流体力学中的相关问题提供可靠的理论支持.