基于Darcy-Stokes 耦合模型的多孔介質(zhì)顆粒懸浮液等效黏性系數(shù)計算1)

胡 洋彭 巍 李德才

* (北京交通大學(xué)機械與電子控制工程學(xué)院,北京 100044)

? (清華大學(xué)摩擦學(xué)國家重點實驗室,北京 100084)

引言

顆粒懸浮液廣泛存在于自然界及工程應(yīng)用領(lǐng)域,其黏性特征對懸浮液的流動行為有著重要的影響[1-4].早在1905 年,Einstein[5]就研究了低濃度條件下球形固體顆粒對流體黏性的影響,給出了懸浮液等效黏度系數(shù)計算的一個強有力的理論框架,并得到了低濃度球形顆粒懸浮液的著名的Einstein 黏性

公式,即 μs=μf(1+Λφ),這里分別 μs,μf和 φ 為懸浮

液等效黏性系數(shù)、液體本身黏性系數(shù)和固體顆粒體積分?jǐn)?shù),Λ=2.5 即為此時的特性黏度.Einstein 的工作引發(fā)了大量的后繼研究,很多學(xué)者對更為復(fù)雜情形下的顆粒懸浮液等效黏性系數(shù)計算進行了研究,主要包括了高濃度條件、考慮慣性效應(yīng)以及非球形顆粒3 個方向的拓展.在高顆粒濃度的研究方面,典型工作包括:Batchelor 和Green[6]首先考慮了固體顆粒間的相互作用,將等效黏性系數(shù)公式拓展至體積分?jǐn)?shù)的二次項,在顆粒大小相等情形下有μs=μf(1+2.5φ+5.2φ2).其后Batchelor[7]進一步考慮了流體中布朗運動的影響,將公式修正為μs=μf(1+2.5φ+6.2φ2).近來,Zhu 等[8]結(jié)合前人的研究成果,提出了一個簡單的等效黏性系數(shù)關(guān)系式,能夠覆蓋較寬的顆粒濃度范圍.在考慮慣性效應(yīng)的研究方面,典型工作包括:Lin 等[9]首先從理論上給出了一個修正公

式,即 μs=μf[1+(2.5+1.34Re1.5)φ],其中Re為雷諾數(shù).Kulkarni 和Morris[10],Yeo 和Maxey[11-12]分別采用介觀和宏觀尺度的直接數(shù)值模擬方法,研究了有限雷諾數(shù)條件下剪切流中顆粒懸浮液的流變特性,并給出了較高濃度條件下等效黏性系數(shù)隨雷諾數(shù)的變化規(guī)律.在非球形顆粒的研究方面,典型工作包括:Jeffery[13]首先考察了低濃度條件下橢球形顆粒懸浮液的等效黏性系數(shù),并得到了長橢球體和扁球體顆粒懸浮液在不同偏心率情形下對應(yīng)的特性黏度.Yamamoto 和Matsuoka[14]采用粒子模擬方法研究了低濃度棒狀顆粒的等效黏性系數(shù),獲得了剛體棒和柔性棒兩種情形下特性黏度隨方位角和軌道參數(shù)的變化規(guī)律.Huang 等[15]采用數(shù)值方法研究了低濃度條件下長橢球體和扁球體顆粒懸浮液的等效黏性系數(shù),結(jié)果表明低雷諾數(shù)時特性黏度和雷諾數(shù)呈線性相關(guān),而在高雷諾數(shù)時呈非線性關(guān)系.

值得注意的是,上述研究涉及的顆粒均為不可滲透性的固體顆粒.由于多孔介質(zhì)材料在眾多領(lǐng)域中得廣泛應(yīng)用[16-17],近年來一些國內(nèi)外學(xué)者也開始關(guān)注多孔介質(zhì)顆粒懸浮液的黏性特征,主要包括理論分析和數(shù)值模擬兩方面的工作.理論分析方面的典型工作包括:Natraj 和Chen[18]研究了帶電多孔介質(zhì)顆粒懸浮液中的電黏性效應(yīng),首次給出了特性黏度和達西數(shù)的關(guān)系式.其后Ohshima[19-22]進一步拓展了Natraj 和Chen[18]的工作,分別研究了較高濃度顆粒懸浮液以及內(nèi)含實心核的多孔介質(zhì)顆粒懸浮液的等效黏性特征.上述工作均采用Darcy-Brinkman模型[23-24]描述多孔介質(zhì)區(qū)域的流場,自由流區(qū)域采用Stokes 方程,界面上則采用速度連續(xù)和應(yīng)力連續(xù)條件[25].Prakash 和Sekhar[26]則引入剪切應(yīng)力跳躍條件,給出了多孔介質(zhì)懸浮液特性黏度與達西數(shù)、剪切應(yīng)力跳躍系數(shù)的關(guān)系.數(shù)值模擬方面的主要工作包括:Xu 等[27]采用介觀尺度的格子Boltzmann方法模擬了含多孔介質(zhì)球的Couette 剪切流問題,并通過直接計算邊界剪切應(yīng)力的方法確定了等效黏性系數(shù)和達西數(shù)、雷諾數(shù)的關(guān)系.Huang 等[28]則采用格子Boltzmann 方法研究了較高濃度、考慮流體慣性效應(yīng)以及邊界限制因素影響下的多孔介質(zhì)懸浮液黏性特征,結(jié)果表明在達西數(shù)較小時,流體慣性和邊界限制對等效黏性系數(shù)有著顯著的影響,而達西數(shù)較大時,流體慣性和邊界限制的影響可忽略不計.Liu 等[29]進一步采用研究了低濃度低雷諾數(shù)條件下橢圓形顆粒懸浮液的等效黏性特征,發(fā)現(xiàn)特性黏度隨達西數(shù)增加而單調(diào)遞減,并給出了特性黏度一個簡單的經(jīng)驗公式.上述數(shù)值方面的工作均是采用介觀尺度數(shù)值方法求解描述多孔介質(zhì)流動的廣義Navier-Stokes 方程[30-31],自由流區(qū)域和多孔介質(zhì)區(qū)域統(tǒng)一處理,界面上不做特殊處理,也就是界面上速度和應(yīng)力仍保持連續(xù).

綜上所述,已有的多孔介質(zhì)顆粒懸浮液等效黏性的研究工作均基于Darcy-Brinkman 模型或更為復(fù)雜的廣義Navier-Stokes 模型.根據(jù)Nield 的分析[32],當(dāng)采用Darcy-Brinkman 模型或廣義Navier-Stokes模型在處理自由流/多孔介質(zhì)流耦合問題時,模型自由度過大,一是模型本身包含了一個未定參數(shù),即所謂的有效黏性系數(shù),同時界面上還需考慮剪切應(yīng)力的跳躍.因此,本工作基于Darcy-Stokes 耦合模型及Beavers-Joseph(-Saffman)界面條件[33-39],研究了低雷諾數(shù)條件下多孔介質(zhì)球?qū)儜?yīng)變流動的流場的影響,給出了自由流和多孔介質(zhì)區(qū)域的流場解析解,根據(jù)流場解析解計算了由于多孔介質(zhì)球的存在額外產(chǎn)生的黏性熱耗散率,獲得了多孔介質(zhì)顆粒懸浮液的一個新的等效黏性系數(shù)計算公式,并與基于Darcy-Brinkman 模型的結(jié)果進行了比較.

1 問題模型和求解方法

為了獲得低濃度多孔介質(zhì)顆粒懸浮液的等效黏度的解析表達式,首先考慮多孔介質(zhì)球?qū)儜?yīng)變流動的影響.如圖1 所示,不可壓縮流場內(nèi)包含球心位于坐標(biāo)軸原點、直徑為 2a的各向同性多孔介質(zhì)球.多孔介質(zhì)區(qū)域 Ω0內(nèi)的流體運動滿足Darcy 方程,即

圖1 包含多孔介質(zhì)球形顆粒的流場區(qū)域的示意圖Fig.1 Schematic diagram of the flow field containing a spherical porous particle

式中upm,i和ppm分別為多孔介質(zhì)區(qū)速度的第i個分量和壓強,μf和K分別為流體黏性系數(shù)和多孔介質(zhì)滲透率.自由流區(qū)域 Ω1的流場控制方程為

其中uf,i,uf,j和pf分別為自由流區(qū)速度的第i,j個分量和壓強.

自由流-多孔介質(zhì)界面上A0的條件為

式中ni,nj分別為沿界面指向自由流區(qū)域單位法向量的第i,j個分量,τs,i為沿界面的切平面一個向量正交系的第i個分量.式(5)和式(6)分別表示界面上質(zhì)量守恒和法向力的平衡.式(7) 即為Beavers-Joseph 界面條件,其中 αBJ即所謂的Beavers-Joseph 系數(shù).式(8)為Saffman 修改的界面條件,即文獻中提到的Beavers-Joseph-Saffman 條件.兩類條件都被用于接下來的計算.無量綱數(shù) αBJ可用于表征自由流?多孔介質(zhì)界面上速度滑移的程度,αBJ越大,兩側(cè)速度滑移則越小.αBJ的值依賴于多孔介質(zhì)界面附近的幾何和結(jié)構(gòu)特征.Beavers 和Joseph 使用泡沫金屬、鋁土制作了平均孔隙尺寸在0.033 cm～ 0.116 cm 之間的多孔介質(zhì)材料,并通過實驗方法確定αBJ的范圍為[0.1,4],一般的數(shù)值研究中采用的 αBJ的值也在該范圍之內(nèi).

考慮未受多孔介質(zhì)圓球擾動的純應(yīng)變流場具有如下線性速度分布

式中eij為常對稱張量.根據(jù)流場的不可壓縮條件,對角元eii應(yīng)滿足

由于流場的對稱性,多孔介質(zhì)小球?qū)⒈３朱o止,同時流場出現(xiàn)了一個擾動 δui,這時速度場為

自由流區(qū)域的速度擾動有如下解[40]

式中A和B為常數(shù).自由流區(qū)域的壓力場的表達式為

圓球內(nèi)多孔介質(zhì)區(qū)域的速度場和壓力場的解為

式中C為常數(shù).為了確定常數(shù)A,B和C的值,將自由流區(qū)域和多孔介質(zhì)區(qū)域速度和壓力的解析表達式(11)～式(16)代入界面條件式(5)～ 式(8),可得

這里以推導(dǎo)式(19)或式(20)為例,給出一些計算細節(jié).應(yīng)變率張量的表達式為

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

同理可得

將f(r)和g(r) 及其導(dǎo)數(shù)的表達式代入式(22)和式(23),即可得到式(19) 或式(20).式(17)～ 式(20) 簡化可得

上述三元一次方程組的解為

2 多孔介質(zhì)懸浮液等效黏性系數(shù)的計算

由上所述,由于多孔介質(zhì)球的存在,流場在線性分布的基礎(chǔ)上增加了一個擾動,從而增加了所考慮的流場邊界A1上的黏性耗散率.另一方面,從等效的觀點看,可以認為多孔介質(zhì)球的存在使得流體的黏性由 μf等效變?yōu)?μs,但流場仍保持原來的線性分布,同時兩種情形下邊界A1上的黏性耗散率相等,即有如下關(guān)系式

式中應(yīng)力張量 σs,i j和 σij的表達式為

應(yīng)用高斯散度定理,流場邊界A1上的積分可化為多孔介質(zhì)球面A0上的積分.根據(jù)Batchelor 的推導(dǎo)[40],式(36)可化為

式中V1為整個流場區(qū)域的體積.將流場擾動 δui的解析表達式代入式(39)右端的被積函數(shù),可得

應(yīng)用如下兩個恒等式

式(39)化為

其中特性黏度 ΛD為

當(dāng)Da→0,即多孔介質(zhì)顆粒變?yōu)槌Ｒ?guī)的實心顆粒,流體無法穿透壁面,ΛD→2.5 成立,式(44)退化為經(jīng)典的Einstein 黏度公式.此外,可以看到兩類界面條件式(7)或式(8)給出的結(jié)果非常接近,差別在O(Da)量級.如無特殊說明,接下來的結(jié)果部分僅采用式(45)進行討論和分析.

注意到等效黏性系數(shù) μs和Beavers-Joseph 系數(shù)αBJ相關(guān),圖2 給出了 αBJ取不同的值時特性黏度ΛD隨達西數(shù)Da的變化.可以看到 αBJ增大時,特性黏度ΛD也隨之增大.此外,αBJ越大,ΛD增加的幅度也越小,αBJ=1.25和 αBJ=1.5 時的兩條曲線非常接近.如前所述,αBJ越大表明自由流?多孔介質(zhì)界面上速度滑移越小,即自由流區(qū)域的速度有較大的擾動,從而導(dǎo)致特性黏度 ΛD增加.另一方面,Da和 ΛD之間呈現(xiàn)非線性變化.當(dāng) 10?6≤Da≤10?4時,多孔介質(zhì)顆粒的滲透率很小,接近于固體實心顆粒,因而特性黏度ΛD接近于2.5.當(dāng)10?4≤Da≤10?1時,此時多孔介質(zhì)球內(nèi)部的流體阻礙明顯減弱,特性黏度 ΛD快速下降,因而等效黏性系數(shù)更加接近于流體本身的黏性,這與已有的理論和數(shù)值方法給出的結(jié)果是一致的.

圖2 αBJ=0.5,0.75,1,1.25和 1.5時特性黏度 ΛD隨達西數(shù) Da 的變化Fig.2 Intrinsic viscosity ΛDversus the Darcy number Da when αBJ=0.5,0.75,1,1.25 and1.5

3 與Darcy-Brinkman 模型的結(jié)果的比較

如前所述,除了Darcy 模型,一些學(xué)者也采用Darcy-Brinkman 模型對多孔介質(zhì)懸浮液的等效黏性系數(shù)進行了計算.Darcy-Brinkman 模型的控制方程為

式中 μeff為有效黏性系數(shù).自由流區(qū)域和多孔介質(zhì)區(qū)域界面上采用如下剪切應(yīng)力跳躍模型[26]

其中 ζ 為界面應(yīng)力跳躍系數(shù).當(dāng) μeff=μf時,Prakash 和Sekhar[26]得出了參數(shù) Λ 如下形式的表達式

將式(43)和式(47)做展開,即

很顯然,如果Beavers-Joseph 系數(shù)和界面應(yīng)力跳躍系數(shù)滿足條件 αBJ=1?ζ,Darcy 模型和Darcy-Brinkman 模型在低達西數(shù)條件下將給出較為一致的結(jié)果.圖3 給出了兩類模型在 10?4≤Da≤10?1范圍內(nèi)特性黏度 Λ 的變化,考慮的3 組參數(shù)( αBJ,ζ) 如下:( 1,0),(0.8,0.2)和 (1.2,?0.2).可以看到達西數(shù)較小時,兩類模型的確給出了幾乎相同的結(jié)果.

圖3 兩類模型在 10?4≤Da≤10?1 范圍內(nèi)特性黏度 Λ 的變化(續(xù))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 范圍內(nèi)特性黏度 Λ 的變化Fig.3 Relationship between the intrinsic viscosity Λ and the Darcy number based on the Darcy model and Darcy-Brinkman model.

4 結(jié)論

本文首先研究了線性分布的流場中多孔介質(zhì)球引起的擾動問題.在低雷諾數(shù)條件下,自由流區(qū)域滿足Stokes 方程,而多孔介質(zhì)球內(nèi)部的流場采用Darcy 模型描述.多孔介質(zhì)球表面滿足質(zhì)量守恒律、法向應(yīng)力平衡條件以及描述速度滑移的Beavers-Joseph 條件或Beavers-Joseph-Saffman 條件.通過求解上述Darcy-Stokes 耦合模型,獲得了自由流/多孔介質(zhì)區(qū)域速度場和壓力場的解析表達式.其后根據(jù)Batchelor 提出的定義計算了低濃度多孔介質(zhì)顆粒懸浮液的等效黏性系數(shù),得到了其與達西數(shù)、Beavers-Joseph 系數(shù)和顆粒體積分?jǐn)?shù)的定量關(guān)系.分析結(jié)果表明等效黏性系數(shù)隨著Beavers-Joseph 系數(shù)增加而增加,在低達西數(shù)條件下多√孔介質(zhì)懸浮液特性黏度有如下漸進公式.此外,所得到的等效黏性計算公式與基于Darcy-Brinkman 模型的結(jié)果進行了比較,結(jié)果顯示當(dāng)Beavers-Joseph 系數(shù)和界面應(yīng)力跳躍系數(shù)之和為1 時,兩類計算公式在低達西數(shù)條件下給出幾乎一致的結(jié)果.