基于Darcy-Stokes 耦合模型的多孔介质颗粒悬浮液等效黏性系数计算1)

胡 洋彭 巍 李德才

* (北京交通大学机械与电子控制工程学院,北京 100044)

† (清华大学摩擦学国家重点实验室,北京 100084)

引言

颗粒悬浮液广泛存在于自然界及工程应用领域,其黏性特征对悬浮液的流动行为有着重要的影响[1-4].早在1905 年,Einstein[5]就研究了低浓度条件下球形固体颗粒对流体黏性的影响,给出了悬浮液等效黏度系数计算的一个强有力的理论框架,并得到了低浓度球形颗粒悬浮液的著名的Einstein 黏性

公式,即 µs=µf(1+Λφ),这里分别 µs,µf和 φ 为悬浮

液等效黏性系数、液体本身黏性系数和固体颗粒体积分数,Λ=2.5 即为此时的特性黏度.Einstein 的工作引发了大量的后继研究,很多学者对更为复杂情形下的颗粒悬浮液等效黏性系数计算进行了研究,主要包括了高浓度条件、考虑惯性效应以及非球形颗粒3 个方向的拓展.在高颗粒浓度的研究方面,典型工作包括:Batchelor 和Green[6]首先考虑了固体颗粒间的相互作用,将等效黏性系数公式拓展至体积分数的二次项,在颗粒大小相等情形下有µs=µf(1+2.5φ+5.2φ2).其后Batchelor[7]进一步考虑了流体中布朗运动的影响,将公式修正为µs=µf(1+2.5φ+6.2φ2).近来,Zhu 等[8]结合前人的研究成果,提出了一个简单的等效黏性系数关系式,能够覆盖较宽的颗粒浓度范围.在考虑惯性效应的研究方面,典型工作包括:Lin 等[9]首先从理论上给出了一个修正公

式,即 µs=µf[1+(2.5+1.34Re1.5)φ],其中Re为雷诺数.Kulkarni 和Morris[10],Yeo 和Maxey[11-12]分别采用介观和宏观尺度的直接数值模拟方法,研究了有限雷诺数条件下剪切流中颗粒悬浮液的流变特性,并给出了较高浓度条件下等效黏性系数随雷诺数的变化规律.在非球形颗粒的研究方面,典型工作包括:Jeffery[13]首先考察了低浓度条件下椭球形颗粒悬浮液的等效黏性系数,并得到了长椭球体和扁球体颗粒悬浮液在不同偏心率情形下对应的特性黏度.Yamamoto 和Matsuoka[14]采用粒子模拟方法研究了低浓度棒状颗粒的等效黏性系数,获得了刚体棒和柔性棒两种情形下特性黏度随方位角和轨道参数的变化规律.Huang 等[15]采用数值方法研究了低浓度条件下长椭球体和扁球体颗粒悬浮液的等效黏性系数,结果表明低雷诺数时特性黏度和雷诺数呈线性相关,而在高雷诺数时呈非线性关系.

值得注意的是,上述研究涉及的颗粒均为不可渗透性的固体颗粒.由于多孔介质材料在众多领域中得广泛应用[16-17],近年来一些国内外学者也开始关注多孔介质颗粒悬浮液的黏性特征,主要包括理论分析和数值模拟两方面的工作.理论分析方面的典型工作包括:Natraj 和Chen[18]研究了带电多孔介质颗粒悬浮液中的电黏性效应,首次给出了特性黏度和达西数的关系式.其后Ohshima[19-22]进一步拓展了Natraj 和Chen[18]的工作,分别研究了较高浓度颗粒悬浮液以及内含实心核的多孔介质颗粒悬浮液的等效黏性特征.上述工作均采用Darcy-Brinkman模型[23-24]描述多孔介质区域的流场,自由流区域采用Stokes 方程,界面上则采用速度连续和应力连续条件[25].Prakash 和Sekhar[26]则引入剪切应力跳跃条件,给出了多孔介质悬浮液特性黏度与达西数、剪切应力跳跃系数的关系.数值模拟方面的主要工作包括:Xu 等[27]采用介观尺度的格子Boltzmann方法模拟了含多孔介质球的Couette 剪切流问题,并通过直接计算边界剪切应力的方法确定了等效黏性系数和达西数、雷诺数的关系.Huang 等[28]则采用格子Boltzmann 方法研究了较高浓度、考虑流体惯性效应以及边界限制因素影响下的多孔介质悬浮液黏性特征,结果表明在达西数较小时,流体惯性和边界限制对等效黏性系数有着显著的影响,而达西数较大时,流体惯性和边界限制的影响可忽略不计.Liu 等[29]进一步采用研究了低浓度低雷诺数条件下椭圆形颗粒悬浮液的等效黏性特征,发现特性黏度随达西数增加而单调递减,并给出了特性黏度一个简单的经验公式.上述数值方面的工作均是采用介观尺度数值方法求解描述多孔介质流动的广义Navier-Stokes 方程[30-31],自由流区域和多孔介质区域统一处理,界面上不做特殊处理,也就是界面上速度和应力仍保持连续.

综上所述,已有的多孔介质颗粒悬浮液等效黏性的研究工作均基于Darcy-Brinkman 模型或更为复杂的广义Navier-Stokes 模型.根据Nield 的分析[32],当采用Darcy-Brinkman 模型或广义Navier-Stokes模型在处理自由流/多孔介质流耦合问题时,模型自由度过大,一是模型本身包含了一个未定参数,即所谓的有效黏性系数,同时界面上还需考虑剪切应力的跳跃.因此,本工作基于Darcy-Stokes 耦合模型及Beavers-Joseph(-Saffman)界面条件[33-39],研究了低雷诺数条件下多孔介质球对纯应变流动的流场的影响,给出了自由流和多孔介质区域的流场解析解,根据流场解析解计算了由于多孔介质球的存在额外产生的黏性热耗散率,获得了多孔介质颗粒悬浮液的一个新的等效黏性系数计算公式,并与基于Darcy-Brinkman 模型的结果进行了比较.

1 问题模型和求解方法

为了获得低浓度多孔介质颗粒悬浮液的等效黏度的解析表达式,首先考虑多孔介质球对纯应变流动的影响.如图1 所示,不可压缩流场内包含球心位于坐标轴原点、直径为 2a的各向同性多孔介质球.多孔介质区域 Ω0内的流体运动满足Darcy 方程,即

图1 包含多孔介质球形颗粒的流场区域的示意图Fig.1 Schematic diagram of the flow field containing a spherical porous particle

式中upm,i和ppm分别为多孔介质区速度的第i个分量和压强,µf和K分别为流体黏性系数和多孔介质渗透率.自由流区域 Ω1的流场控制方程为

其中uf,i,uf,j和pf分别为自由流区速度的第i,j个分量和压强.

自由流-多孔介质界面上A0的条件为

式中ni,nj分别为沿界面指向自由流区域单位法向量的第i,j个分量,τs,i为沿界面的切平面一个向量正交系的第i个分量.式(5)和式(6)分别表示界面上质量守恒和法向力的平衡.式(7) 即为Beavers-Joseph 界面条件,其中 αBJ即所谓的Beavers-Joseph 系数.式(8)为Saffman 修改的界面条件,即文献中提到的Beavers-Joseph-Saffman 条件.两类条件都被用于接下来的计算.无量纲数 αBJ可用于表征自由流−多孔介质界面上速度滑移的程度,αBJ越大,两侧速度滑移则越小.αBJ的值依赖于多孔介质界面附近的几何和结构特征.Beavers 和Joseph 使用泡沫金属、铝土制作了平均孔隙尺寸在0.033 cm～ 0.116 cm 之间的多孔介质材料,并通过实验方法确定αBJ的范围为[0.1,4],一般的数值研究中采用的 αBJ的值也在该范围之内.

考虑未受多孔介质圆球扰动的纯应变流场具有如下线性速度分布

式中eij为常对称张量.根据流场的不可压缩条件,对角元eii应满足

由于流场的对称性,多孔介质小球将保持静止,同时流场出现了一个扰动 δui,这时速度场为

自由流区域的速度扰动有如下解[40]

式中A和B为常数.自由流区域的压力场的表达式为

圆球内多孔介质区域的速度场和压力场的解为

式中C为常数.为了确定常数A,B和C的值,将自由流区域和多孔介质区域速度和压力的解析表达式(11)～式(16)代入界面条件式(5)～ 式(8),可得

这里以推导式(19)或式(20)为例,给出一些计算细节.应变率张量的表达式为

利用恒等式nkxk=r和 τs,kxk=0,则可得

同理可得

将f(r)和g(r) 及其导数的表达式代入式(22)和式(23),即可得到式(19) 或式(20).式(17)～ 式(20) 简化可得

上述三元一次方程组的解为

2 多孔介质悬浮液等效黏性系数的计算

由上所述,由于多孔介质球的存在,流场在线性分布的基础上增加了一个扰动,从而增加了所考虑的流场边界A1上的黏性耗散率.另一方面,从等效的观点看,可以认为多孔介质球的存在使得流体的黏性由 µf等效变为 µs,但流场仍保持原来的线性分布,同时两种情形下边界A1上的黏性耗散率相等,即有如下关系式

式中应力张量 σs,i j和 σij的表达式为

应用高斯散度定理,流场边界A1上的积分可化为多孔介质球面A0上的积分.根据Batchelor 的推导[40],式(36)可化为

式中V1为整个流场区域的体积.将流场扰动 δui的解析表达式代入式(39)右端的被积函数,可得

应用如下两个恒等式

式(39)化为

其中特性黏度 ΛD为

当Da→0,即多孔介质颗粒变为常规的实心颗粒,流体无法穿透壁面,ΛD→2.5 成立,式(44)退化为经典的Einstein 黏度公式.此外,可以看到两类界面条件式(7)或式(8)给出的结果非常接近,差别在O(Da)量级.如无特殊说明,接下来的结果部分仅采用式(45)进行讨论和分析.

注意到等效黏性系数 µs和Beavers-Joseph 系数αBJ相关,图2 给出了 αBJ取不同的值时特性黏度ΛD随达西数Da的变化.可以看到 αBJ增大时,特性黏度ΛD也随之增大.此外,αBJ越大,ΛD增加的幅度也越小,αBJ=1.25和 αBJ=1.5 时的两条曲线非常接近.如前所述,αBJ越大表明自由流−多孔介质界面上速度滑移越小,即自由流区域的速度有较大的扰动,从而导致特性黏度 ΛD增加.另一方面,Da和 ΛD之间呈现非线性变化.当 10−6≤Da≤10−4时,多孔介质颗粒的渗透率很小,接近于固体实心颗粒,因而特性黏度ΛD接近于2.5.当10−4≤Da≤10−1时,此时多孔介质球内部的流体阻碍明显减弱,特性黏度 ΛD快速下降,因而等效黏性系数更加接近于流体本身的黏性,这与已有的理论和数值方法给出的结果是一致的.

图2 αBJ=0.5,0.75,1,1.25和 1.5时特性黏度 ΛD随达西数 Da 的变化Fig.2 Intrinsic viscosity ΛDversus the Darcy number Da when αBJ=0.5,0.75,1,1.25 and1.5

3 与Darcy-Brinkman 模型的结果的比较

如前所述,除了Darcy 模型,一些学者也采用Darcy-Brinkman 模型对多孔介质悬浮液的等效黏性系数进行了计算.Darcy-Brinkman 模型的控制方程为

式中 µeff为有效黏性系数.自由流区域和多孔介质区域界面上采用如下剪切应力跳跃模型[26]

其中 ζ 为界面应力跳跃系数.当 µeff=µf时,Prakash 和Sekhar[26]得出了参数 Λ 如下形式的表达式

将式(43)和式(47)做展开,即

很显然,如果Beavers-Joseph 系数和界面应力跳跃系数满足条件 αBJ=1−ζ,Darcy 模型和Darcy-Brinkman 模型在低达西数条件下将给出较为一致的结果.图3 给出了两类模型在 10−4≤Da≤10−1范围内特性黏度 Λ 的变化,考虑的3 组参数( αBJ,ζ) 如下:( 1,0),(0.8,0.2)和 (1.2,−0.2).可以看到达西数较小时,两类模型的确给出了几乎相同的结果.

图3 两类模型在 10−4≤Da≤10−1 范围内特性黏度 Λ 的变化(续)Fig.3 Relationship between the intrinsic viscosity Λ and the Darcy number based on the Darcy model and Darcy-Brinkman model (continued)

图3 两类模型在 10−4≤Da≤10−1 范围内特性黏度 Λ 的变化Fig.3 Relationship between the intrinsic viscosity Λ and the Darcy number based on the Darcy model and Darcy-Brinkman model.

4 结论

本文首先研究了线性分布的流场中多孔介质球引起的扰动问题.在低雷诺数条件下,自由流区域满足Stokes 方程,而多孔介质球内部的流场采用Darcy 模型描述.多孔介质球表面满足质量守恒律、法向应力平衡条件以及描述速度滑移的Beavers-Joseph 条件或Beavers-Joseph-Saffman 条件.通过求解上述Darcy-Stokes 耦合模型,获得了自由流/多孔介质区域速度场和压力场的解析表达式.其后根据Batchelor 提出的定义计算了低浓度多孔介质颗粒悬浮液的等效黏性系数,得到了其与达西数、Beavers-Joseph 系数和颗粒体积分数的定量关系.分析结果表明等效黏性系数随着Beavers-Joseph 系数增加而增加,在低达西数条件下多√孔介质悬浮液特性黏度有如下渐进公式.此外,所得到的等效黏性计算公式与基于Darcy-Brinkman 模型的结果进行了比较,结果显示当Beavers-Joseph 系数和界面应力跳跃系数之和为1 时,两类计算公式在低达西数条件下给出几乎一致的结果.