基于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 時,兩類計算公式在低達西數條件下給出幾乎一致的結果.