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

胡 洋彭 巍 李德才

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

? (清華大學(xué)摩擦學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100084)

引言

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

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

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

式,即 μ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ù),并得到了長(zhǎng)橢球體和扁球體顆粒懸浮液在不同偏心率情形下對(duì)應(yīng)的特性黏度.Yamamoto 和Matsuoka[14]采用粒子模擬方法研究了低濃度棒狀顆粒的等效黏性系數(shù),獲得了剛體棒和柔性棒兩種情形下特性黏度隨方位角和軌道參數(shù)的變化規(guī)律.Huang 等[15]采用數(shù)值方法研究了低濃度條件下長(zhǎng)橢球體和扁球體顆粒懸浮液的等效黏性系數(shù),結(jié)果表明低雷諾數(shù)時(shí)特性黏度和雷諾數(shù)呈線(xiàn)性相關(guān),而在高雷諾數(shù)時(shí)呈非線(xiàn)性關(guān)系.

值得注意的是,上述研究涉及的顆粒均為不可滲透性的固體顆粒.由于多孔介質(zhì)材料在眾多領(lǐng)域中得廣泛應(yīng)用[16-17],近年來(lái)一些國(guó)內(nèi)外學(xué)者也開(kāi)始關(guān)注多孔介質(zhì)顆粒懸浮液的黏性特征,主要包括理論分析和數(shù)值模擬兩方面的工作.理論分析方面的典型工作包括:Natraj 和Chen[18]研究了帶電多孔介質(zhì)顆粒懸浮液中的電黏性效應(yīng),首次給出了特性黏度和達(dá)西數(shù)的關(guān)系式.其后Ohshima[19-22]進(jìn)一步拓展了Natraj 和Chen[18]的工作,分別研究了較高濃度顆粒懸浮液以及內(nèi)含實(shí)心核的多孔介質(zhì)顆粒懸浮液的等效黏性特征.上述工作均采用Darcy-Brinkman模型[23-24]描述多孔介質(zhì)區(qū)域的流場(chǎng),自由流區(qū)域采用Stokes 方程,界面上則采用速度連續(xù)和應(yīng)力連續(xù)條件[25].Prakash 和Sekhar[26]則引入剪切應(yīng)力跳躍條件,給出了多孔介質(zhì)懸浮液特性黏度與達(dá)西數(shù)、剪切應(yīng)力跳躍系數(shù)的關(guān)系.數(shù)值模擬方面的主要工作包括:Xu 等[27]采用介觀尺度的格子Boltzmann方法模擬了含多孔介質(zhì)球的Couette 剪切流問(wèn)題,并通過(guò)直接計(jì)算邊界剪切應(yīng)力的方法確定了等效黏性系數(shù)和達(dá)西數(shù)、雷諾數(shù)的關(guān)系.Huang 等[28]則采用格子Boltzmann 方法研究了較高濃度、考慮流體慣性效應(yīng)以及邊界限制因素影響下的多孔介質(zhì)懸浮液黏性特征,結(jié)果表明在達(dá)西數(shù)較小時(shí),流體慣性和邊界限制對(duì)等效黏性系數(shù)有著顯著的影響,而達(dá)西數(shù)較大時(shí),流體慣性和邊界限制的影響可忽略不計(jì).Liu 等[29]進(jìn)一步采用研究了低濃度低雷諾數(shù)條件下橢圓形顆粒懸浮液的等效黏性特征,發(fā)現(xiàn)特性黏度隨達(dá)西數(shù)增加而單調(diào)遞減,并給出了特性黏度一個(gè)簡(jiǎn)單的經(jīng)驗(yàn)公式.上述數(shù)值方面的工作均是采用介觀尺度數(shù)值方法求解描述多孔介質(zhì)流動(dòng)的廣義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ì)流耦合問(wèn)題時(shí),模型自由度過(guò)大,一是模型本身包含了一個(gè)未定參數(shù),即所謂的有效黏性系數(shù),同時(shí)界面上還需考慮剪切應(yīng)力的跳躍.因此,本工作基于Darcy-Stokes 耦合模型及Beavers-Joseph(-Saffman)界面條件[33-39],研究了低雷諾數(shù)條件下多孔介質(zhì)球?qū)儜?yīng)變流動(dòng)的流場(chǎng)的影響,給出了自由流和多孔介質(zhì)區(qū)域的流場(chǎng)解析解,根據(jù)流場(chǎng)解析解計(jì)算了由于多孔介質(zhì)球的存在額外產(chǎn)生的黏性熱耗散率,獲得了多孔介質(zhì)顆粒懸浮液的一個(gè)新的等效黏性系數(shù)計(jì)算公式,并與基于Darcy-Brinkman 模型的結(jié)果進(jìn)行了比較.

1 問(wèn)題模型和求解方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

同理可得

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

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

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

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

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

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

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

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

式(39)化為

其中特性黏度 ΛD為

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

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

圖2 αBJ=0.5,0.75,1,1.25和 1.5時(shí)特性黏度 ΛD隨達(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 模型對(duì)多孔介質(zhì)懸浮液的等效黏性系數(shù)進(jìn)行了計(jì)算.Darcy-Brinkman 模型的控制方程為

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

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

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

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

圖3 兩類(lèi)模型在 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 兩類(lèi)模型在 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é)論

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