非均匀流体的体积黏度: Maxwell 弛豫模型*

孙宗利 康艳霜 张君霞

1) (华北电力大学数理系,保定 071003)

2) (华北电力大学数理系,河北省物理学与能源技术重点实验室,保定 071003)

3) (河北农业大学理学院,保定 071001)

基于黏弹力学中的Maxwell 弛豫理论,提出一种新颖的、可用于计算非均匀体积黏度的理论方法.该方法通过体相流体的黏性和弹性特征来计算体系的局域弛豫时间,进而结合体系的局域弛豫模量来计算体积黏度的非均匀分布.作为应用,计算了受限于平行狭缝中的Lennard-Jones 流体的非均匀体积黏度,系统地研究了体密度、温度、缝宽和吸附强度等因素对体积黏度的影响.结果表明,上述因素的变化均可显著地调制狭缝中流体的体积黏度.其中,体密度和吸附强度的增大均有益于体积黏度的增强,而温度的升高将则会削弱其体积黏度.此外,毛细凝聚的发生也对受限流体的体积黏度具有显著的调制作用.

1 引言

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

超声吸附技术的发展为实验研究流体的体积黏度提供了可能的途径.根据声学知识[13],体积黏度可由声波吸附系数算得.然而,计算中所涉及的等体和等压热容、零频声速、热传导和剪切黏度等物理量使得准确测定流体的体积黏度极为困难[14,15].相比之下,分子动力学模拟在预测体积黏度方面显得更具功效.采用平衡分子动力学(equilibrium molecular dynamics,EMD)方法[16–18],体积黏度可以结合Green-Kubo 关系或Einstein 关系,并通过计算应力张量的自相关函数在t→∞时的极限值来获得.而在非平衡分子动力学(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 中,令并联的上、下两支中具有相同的体积应变,并记为εv.进一步假设上、下支中的压缩应 力分别为Πv1和Πv2,则总应力为Πv=Πv1+Πv2.从物理角度来看,上支描述体系中的弹性成分,即:

图1 Maxwell 黏弹模型示意图 (a) 弹簧单元;(b) 活塞单元;(c) 标准线性固体模型Fig.1.A sketch of the Maxwell viscoelastic model: (a) Spring unit;(b) piston unit;(c) standard linear solid model.

式中κ0为体积模量的非弛豫成分.事实上,下支中的弹簧和活塞组合即为Maxwell 模型的力学结构.因二者串联,故弹簧和活塞中的应力均为Πv2.与此同时,下支中的应变率应为弹簧和活塞中的应变率之和,即:

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

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

对简谐压缩应变εv∼eiωt而言,且A=1+iωτv.据此,体系的复弹性K∗和复黏度η∗可分别表示为

按定义,K∗(ω)和η∗(ω) 的实部应分别为体系的弹性系数Kv(ω) 和黏度系数ηv(ω),即:

事实上,(7a)式和(7b)式中的κ0和κ1可由其在频率 空间中 的边界 条件得 出.当ω→0 时,Kv(0)=K0;当ω→∞时,Kv(∞)=K∞.其中,K0和K∞分别为体系的零频和高频体积模量.据此可得:κ0=K0,κ1=K2≡K∞-K0.K2为体积模量的弛豫贡献.据此,在零频极限下,体系的体积黏度系数ηv可表示为

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

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

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

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

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

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

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

图2 狭缝中LJ 流体的局域弛豫模量的分布.图中实线为 K2(z)的结果,虚线为 的 结果.计算中 的参数取为 T∗=1.5,H∗=6.0 .此 外,约化模 量Fig.2.Profiles of the local relaxation modulus of LJ fluid in slits.In the figure,the solid and dashed lines stand for the results of K2(z) and ,respectively.In the calculations,the parameters are set as T∗=1.5,H∗=6.0 .In addition,the modulus is reduced as

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

其中cT=0.05,其他拟合系数ci1和ci2在表1 中给出.需指出的是,(14)式中各量均采用无量纲形式:,ρ∗=ρσ3,T∗=kBT/ε.m,σ和ε分别为流体粒子的质量、直径和作用强度,kB为Boltzmann 常数,T为绝对温度.

表1 (14)式中的拟合系数 ci1和ci2Table 1.Fitting parameters of ci1 and ci2 in the Eq.(14).

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

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

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

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

这里φ(r) 为流体粒子两体势,φ(m)(r)为φ(r)的m阶导数,ρ(2)(r1,r2)=ρ(r1)ρ(r2)g(r12) 为两体密度函数,g(r12) 为径向分布函数,且r12=|r2-r1| .对于体相LJ 流体而言,(15)式可由其压强P和内能E表示为

其中ρb为体密度;P和E可由MBWR 状态方程给出[35].

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

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

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

式中,γb和分别为体相比热比和非均匀体系中的局域比热比.就本文所研究的LJ 流体而言,二者均可借助MBWR 状态方程予以计算.

5 结果与讨论

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

图3 狭缝中LJ 流体的剪切黏度分布.计算中的流体参数取为 =0.291,T∗=2.0Fig.3.Profiles of shear viscosity of LJ fluid in slits.In the calculations,the fluid parameters are set as =0.291 and T∗=2.0 .

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

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

图4 狭缝中LJ 流体的平均体积黏度随缝宽 H∗ 的变化.其中,实线、虚线和点分别为本文结果、LADM 结果和基于GK 的解析结果Fig.4.Dependence of the averaged volume viscosity of LJ fluid on the pore width H .In each figure,the solid lines,dashed lines and symbols stand for the results from this work,LADM and GK-based method,respectively.

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

其中,εfs和σfs分别为固-液作用的强度和直径,ρs和∆分别为固体原子面密度和层间距.此外,εfs和σfs均满足Lorentz 混合规则.据此,狭缝中的外势可表达为

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

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

图5 在 H∗=6.0 情况下,体积黏度 (z∗) 随体密度的变化 (a) T∗=2.0 ;(b) T∗=1.0Fig.5.Influence of bulk density on the volume viscosity(z∗) under the condition of H∗=6.0 : (a) T∗=2.0;(b) T∗=1.0 .

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

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

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

不难理解,受限流体发生体积形变和剪切形变的相对难易程度将会受到狭缝的影响与调制.然而,狭缝的限制效应和流体的相态变化对这一特征的调制机制仍不明晰.为此,定义.在此基础上,将在的范围内分别计算T∗=2.0和T∗=1.0 情况下受限流体的的数值,计算结果在图6 中给出.

图6 在 H∗=6.0 情况下,受限流 体中比值 随体密度的变化Fig.6.Influence of bulk density on the ratio of the confined fluids,under the condition of H∗=6.0 .

5.2 温度对体积黏度的影响

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

图8 在 H∗=6.0 情况下,受限流体中比值 随温度的变化Fig.8.Influence of temperature on the ratio of the confined fluids,under the condition of H∗=6.0 .

5.3 缝宽对体积黏度的影响

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

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

图9 在 T∗=1.5 和=0.8 条件下,(a) 体积黏度 (z∗)随缝宽的变化;(b) 比值 随缝宽的变化.图(b)中虚线为同一条件下体相液态的实验结果[40]Fig.9.Influence of pore width on (a) the volume viscosity (z∗) and (b) the ratio ,under the conditions of T∗=1.5 and =0.8 .The dashed line in panel (b) denotes the experimental result[40] under the same conditions.

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

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

狭缝两壁的吸附强度也是影响流体结构与性质的另一重要外因.为了进一步理解狭缝的自身特征对其中所限流体体积黏度的影响,本文将在T∗=1.5和H∗=6.0的条件下计算吸附强度对的影响.考虑到狭缝中可能引发的毛细凝聚及其对体积黏度的调制,计算中已选取两个不同的体系.相关的计算结果可见图10 和图11.

图10 在 H∗=6.0 情况下,体积黏度 (z∗) 随吸附势强度的变化 (a) =0.1 ;(b) =0.6Fig.10.Influence of adsorption strength on the volume viscosity (z∗) under the condition of H∗=6.0 : (a) =0.1 ;(b) =0.6 .

图11 在 H∗=6.0 情况下,受限流 体中比 值随吸附势强度的变化Fig.11.Influence of adsorption strength on the ratio of the confined fluids,under the condition of H∗=6.0 .

6 结论

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

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