基于格子Boltzmann 方法的液潤(rùn)表面減阻規(guī)律

秦聲雷，侯國(guó)祥，郭文強(qiáng)，周斌斌，姜思遠(yuǎn)

華中科技大學(xué) 船舶與海洋工程學(xué)院，湖北 武漢 430074

0 引 言

經(jīng)典的流體力學(xué)假設(shè)固液邊界為無(wú)滑移邊界條件，即壁面的切向速度恒為零，固液界面不發(fā)生滑移。這種假設(shè)在宏觀上得到驗(yàn)證，并被廣泛應(yīng)用于科學(xué)分析和工程應(yīng)用中。但是有研究表明[1-3]在微尺度流動(dòng)中，界面上存在滑移現(xiàn)象。吳承偉等[4]對(duì)流動(dòng)邊界滑移問(wèn)題進(jìn)行了綜述，邊界滑移現(xiàn)象受諸多因素影響，如表面粗糙度、液體粘性、接觸角、速度梯度和溫度等[5-8]。格子Boltzmann方法（LBM）基于分子動(dòng)理學(xué)，作為一種介觀方法，適用于對(duì)微尺度流動(dòng)的模擬。近年來(lái)，黃橋高等[9]應(yīng)用LBM 的偽勢(shì)模型[10-11]對(duì)超疏水表面減阻特性進(jìn)行研究，得出表觀滑移發(fā)生在近壁區(qū)低密度層上的結(jié)論。王凱等[12]應(yīng)用LBM，對(duì)固液界面滑移現(xiàn)象的影響因素進(jìn)行了研究，表明滑移長(zhǎng)度不顯著依賴于剪切率。

傳統(tǒng)的減阻理論認(rèn)為，超疏水表面在水下可以形成氣膜層而不與外層液體直接接觸，原來(lái)的固液接觸面變?yōu)闅庖航佑|面，使得表面承受的剪切應(yīng)力大大減小，并在氣液界面產(chǎn)生滑移速度。任劉珍等[13]綜述并總結(jié)了超疏水表面減阻所面臨的問(wèn)題——在高剪切速率和高壓作用下氣液界面會(huì)發(fā)生失穩(wěn)，從而導(dǎo)致減阻效果不佳。呂鵬宇等[14]對(duì)氣液界面的穩(wěn)定性和演化規(guī)律進(jìn)行了總結(jié)。事實(shí)上，將微結(jié)構(gòu)中的氣體替換為潤(rùn)滑油，可以在保持良好減阻特性的同時(shí)，還較大地保證滑移界面的穩(wěn)定性。Wong 等[15]合成了一種液潤(rùn)多孔表面（SLIPS），發(fā)現(xiàn)這種表面能使液滴具有極小的滾動(dòng)角、可以承受高壓和沖擊（氮?dú)猸h(huán)境中達(dá)到676 atm）、具有良好的自我修復(fù)性和光學(xué)透明性。這種新型表面被稱為液潤(rùn)表面（liquidinfused surface，LIS），在抗冰、防污、減阻等研究領(lǐng)域引起了廣泛的關(guān)注。在減阻方面，Solomon等[16]通過(guò)改變油液粘性比，在動(dòng)力粘性比為260 時(shí)，得到了18 μm 的滑移長(zhǎng)度，減阻效率為16%。Rosenberg 等[17]設(shè)計(jì)了一組粘性比為0.3～1.4，溝槽尺度為100～800 μm 的液潤(rùn)表面，發(fā)現(xiàn)當(dāng)粘性比越大、溝槽寬度越寬、雷諾數(shù)越大時(shí)，減阻效果越好，但在個(gè)別極端的情況下會(huì)出現(xiàn)反常。Chang等[18]的研究也表明減阻效果會(huì)隨溝槽寬度的增加而增加，但受溝槽深度的影響較小。

現(xiàn)有的研究主要是針對(duì)液潤(rùn)表面的溝槽形式、尺寸、油液粘性比等對(duì)減阻效果的影響，研究的假設(shè)也都是油液絕對(duì)不溶。事實(shí)上，不存在絕對(duì)不溶的物質(zhì)，難溶具有相對(duì)性。本文基于LBM，主要從潤(rùn)滑油難溶性方面，對(duì)潤(rùn)滑油的減阻規(guī)律進(jìn)行探究，為設(shè)計(jì)適當(dāng)?shù)奈⒈砻娼Y(jié)構(gòu)提供參考意見(jiàn)。

1 數(shù)值方法

1.1 格子Boltzmann 方法

采用單松弛的格子Boltzmann 方法（LBGK）[19]，其不含源項(xiàng)演化方程為

1.2 Shan-Chen 多組分模型

1.3 模型有效性驗(yàn)證

組分間分子作用力的引入引起了動(dòng)量的變化，上式添加力項(xiàng)有

設(shè)置150 ×150 的計(jì)算域，中心處設(shè)置半徑不等的圓域（圖1）。圓域內(nèi)設(shè)置組分1 密度為1.0、組分2 密度為0，圓域外設(shè)置組分2 密度為1.0、組分1 密度為0，組分間作用強(qiáng)度Gc設(shè)置為1.85。左右為周期邊界，上下為非平衡態(tài)反彈格式[21]。最終得到內(nèi)外壓強(qiáng)差dp和半徑倒數(shù)1/r的關(guān)系，如圖2 所示（本文數(shù)值模擬中變量單位均為格子單位，直線為輔助視圖）。圖2 中dp和1/r之間為線性關(guān)系，且連線經(jīng)過(guò)原點(diǎn)，這表明數(shù)值模擬結(jié)果與式（26）一致。據(jù)此，可以計(jì)算出此時(shí)的表面張力為0.085。由理論分析和數(shù)值模擬結(jié)果可知，此多組分模型能夠有效反映不同組分間的作用力，可以用來(lái)模擬油液分界面上產(chǎn)生的疏水減阻現(xiàn)象。

圖1 計(jì)算模型示意圖Fig. 1 Sketch of computational model

圖2 G c=1.85時(shí)液滴dp 和1/r 的關(guān)系Fig. 2 The relationship between pressure difference dp and reciprocal radius 1/r with parameter Gc=1.85

2 數(shù)值模擬

2.1 壁面潤(rùn)濕性數(shù)值模擬

壁面的潤(rùn)濕性可由接觸角GA來(lái)描述，當(dāng)GA大于90°時(shí)表現(xiàn)為疏水，小于90°時(shí)表現(xiàn)為親水。其疏水特性可以通過(guò)調(diào)整壁面作用強(qiáng)度Gads,σ實(shí)現(xiàn)模擬。采用200×100 的計(jì)算區(qū)域，左右采取周期邊界，上下采取非平衡態(tài)反彈格式（圖3）。初始化時(shí)，在距壁面0.85R處放置一個(gè)半徑R=40的圓域。在圓內(nèi)設(shè)置組分1 的密度為1.0，組分2的密度為0。在圓外設(shè)置組分1 密度為0，組分2的密度為1.0。組分分子間作用強(qiáng)度Gc取1.85，弛豫時(shí)間均取1.0。得到結(jié)果如圖4 和圖5 所示。

圖3 計(jì)算模型示意圖Fig. 3 Sketch of computational model

圖4 不同壁面參數(shù)G ads,σ時(shí)的接觸角Fig. 4 Contact angle with different wall parameters Gads,σ

圖5 不同壁面參數(shù)G ads,1時(shí)的接觸角Fig. 5 Contact angle with different wall parameters Gads,1

從圖5 可以觀察到，隨著壁面參數(shù)Gads,1的不斷變大，其對(duì)應(yīng)的接觸角也不同程度地增大。結(jié)果表明接觸角的大小與壁面參數(shù)Gads,1有正相關(guān)關(guān)系。從物理層面上來(lái)說(shuō)，壁面參數(shù)Gads,1反映固液界面排斥力的大小，Gads,1越大，排斥力越大，壁面疏水性自然也越強(qiáng)。因此接觸角與Gads,1具有正相關(guān)關(guān)系。

2.2 潤(rùn)滑油溶解變化數(shù)值模擬

在Shan-Chen 模型[20]中Gc存在一個(gè)臨界值Gc,crit=1/(ρ1+ρ2)， 當(dāng)Gc小 于Gc,crit時(shí)，兩組分會(huì)發(fā)生混溶，并產(chǎn)生穩(wěn)定的溶解密度解。當(dāng)Gc大于Gc,crit時(shí)，組分之間會(huì)發(fā)生分離并形成穩(wěn)定界面。由于計(jì)算模型并不復(fù)雜，且橫向采用周期邊界條件，為節(jié)省計(jì)算資源，采用較小的30×100 的計(jì)算域，上下采用非平衡態(tài)反彈格式，Gads,σ取為0。初始設(shè)置為上半平面只有組分1，下半平面只有組分2（潤(rùn)滑油），密度皆為1.0（圖6）。得到結(jié)果如圖7 和圖8（a）所示。

圖6 計(jì)算模型示意圖Fig. 6 Sketch of computational model

圖7 潤(rùn)滑油密度分布圖Fig. 7 The distribution of lubricant

圖8 潤(rùn)滑油在上半平面的溶解密度 ρ2和 組分間作用強(qiáng)度 Gc 的關(guān)系圖（虛線為輔助視圖）Fig. 8 The relationship between parameter G c and lubricant's dissolved density in half upper space (the dotted line is for auxiliary view)

從圖7(a)中可以看出，當(dāng)Gc=1.0，潤(rùn)滑油完全溶于組分1 中，兩組分之間不能形成穩(wěn)定界面。當(dāng)Gc>1.0，組分之間可以形成穩(wěn)定界面，且潤(rùn)滑油在上半平面的溶解密度隨著Gc的增加迅速下降，并在Gc=1.8 附近逐漸趨于零。圖8（b）重畫(huà)了文獻(xiàn)[22] 中 ρtotal=2的 部 分 圖 像（原 文 ρtotal用 ρi表 示，本文 ρtotal=1），對(duì)比結(jié)果，本文與文獻(xiàn)[22]一致。

2.3 微結(jié)構(gòu)液潤(rùn)表面滑移流動(dòng)數(shù)值模擬

為了驗(yàn)證液潤(rùn)表面確實(shí)具有減阻效果，我們?cè)O(shè)置了計(jì)算域?yàn)?30×250，溝槽深度為20，寬度為30，溝槽內(nèi)設(shè)置組分2 為潤(rùn)滑油，外部流場(chǎng)設(shè)置為組分1。為簡(jiǎn)便起見(jiàn)，各組分密度均設(shè)置為1.0，初始溶解密度均為0.02。左右為周期邊界，壁面采用標(biāo)準(zhǔn)反彈格式。為了滿足LBGK 的低馬赫數(shù)條件，頂部（Y=250）水平速度設(shè)置為0.01（圖9），并采用非平衡態(tài)外推格式[23]。取Gc= 1.8，Gads,1=0.2，Gads,2=?0.2。

圖9 計(jì)算模型示意圖Fig. 9 Sketch of computing model

達(dá)到平衡狀態(tài)時(shí)，外部流動(dòng)區(qū)域組分1 密度為0.991，組分2 溶解密度為0.031。溝槽內(nèi)組分2密度為0.98，組分1 溶解密度約為0.032 7。將收斂后計(jì)算域中所有組分2 密度求和約為5 586.44，與初始設(shè)置的組分2 總質(zhì)量29×20×7+0.02×230×331=5 582.6 對(duì)比，相對(duì)誤差為 6.87×10?4，仍然滿足質(zhì)量守恒。模擬結(jié)果如圖10 所示，圖中微溝槽用黑色標(biāo)記。

圖10 微結(jié)構(gòu)密度分布圖Fig. 10 The density distribution of the microstructure with microflow

從圖10（c）可以看到，在微結(jié)構(gòu)及潤(rùn)滑油的表面上明顯存在一層低密度層，而表觀滑移則是發(fā)生在低密度層之上[9]。對(duì)低密度層外的速度進(jìn)行線性擬合并向?qū)嶋H流場(chǎng)區(qū)域外擴(kuò)展，速度等于0 的位置與壁面的距離即為滑移長(zhǎng)度。最終得到中段處滑移長(zhǎng)度為1.06。數(shù)值模擬結(jié)果表明液潤(rùn)表面確實(shí)可以產(chǎn)生滑移現(xiàn)象，具有減阻的作用。為了進(jìn)一步探討潤(rùn)滑油溶解性質(zhì)對(duì)滑移長(zhǎng)度的影響，我們簡(jiǎn)化了模型——Couette 流，將在后文介紹。

2.4 溶解密度對(duì)滑移長(zhǎng)度的影響

在Shan-Chen 模型中，溶解密度和組分分子間作用強(qiáng)度Gc有 關(guān)，Gc越大組分分子間作用強(qiáng)度越強(qiáng)，溶解密度越小。因此可以用Gc來(lái)表征組分的溶解性質(zhì)。溶解密度和Gc的關(guān)系可以由圖8（a）得到。為了減小壁面影響，僅考慮潤(rùn)滑油溶解密度對(duì)滑移長(zhǎng)度的影響，設(shè)置壁面參數(shù)Gads,1和Gads,2為0。此時(shí)壁面表現(xiàn)為既不親水也不疏水。算例為30×100 的Couette 流，頂部剪切速度恒為0.01，上半平面只有組分1，下半平面只有組分2（潤(rùn)滑油），密度皆為1.0，左右為周期邊界，上下為非平衡態(tài)反彈邊界（圖11）。通過(guò)改變組分間作用強(qiáng)度Gc的大小，可得圖12 所示的關(guān)系。值得注意的是，曲線在Gc=1.0時(shí)出現(xiàn)了斜率不連續(xù)點(diǎn)，這是因?yàn)閮煞N液體完全混溶，分界面突然消失（圖7（a））。

圖11 計(jì)算模型示意圖Fig. 11 Sketch of computational model

圖12 中曲線第1 段具有較好的線性規(guī)律，此時(shí)2 種液體完全混溶，不存在分界面。當(dāng)Gc=0時(shí)，算例等效于單組分的Couette 流，滑移長(zhǎng)度b為0。當(dāng)Gc≠0時(shí)，由于壁面附近粒子受上方其他組分粒子的排斥力，組分1 和2 的粒子會(huì)聚集在壁面，造成壁面總密度較大（圖13（a）～圖13（c）和圖14），形成一層高密度區(qū)，表現(xiàn)為增阻。組分間作用強(qiáng)度越大，壁面總密度越高（圖13（a）～圖13（c）），滑移長(zhǎng)度越小。

圖12 G ads,1=?Gads,2=0 時(shí)滑移長(zhǎng)度b 和組分間作用強(qiáng)度G c 的關(guān)系（直線為輔助視圖）Fig. 12 The relationship between slip length b and cohesion force strength Gc with wall parameter Gads,1=?Gads,2=0 (the straight line is for auxiliary view)

圖13 不同組分間作用強(qiáng)度G c時(shí)的總密度分布圖Fig. 13 Total density distribution with different cohesive force strength Gc

圖14 G c=0.1， Gads,1=?Gads,2=0時(shí) 密 度 沿y 方 向 分 布 圖（ ρ1和 ρ2分別代表組分1、組分2 密度）Fig. 14 Density value along y axis with parameters Gc=0，Gads,1=?Gads,2=0

曲線的第2 段約從Gc=1.0 到Gc=1.8，為非線性。導(dǎo)致非線性的原因主要是2 個(gè)方面。一方面是因?yàn)榇穗A段溶解密度的下降表現(xiàn)為非線性（圖8），另一方面壁面附近的密度層逐漸從高密度層轉(zhuǎn)為低密度層（后文將說(shuō)明）。此階段組分在彼此之間的溶解密度迅速下降（圖8），分界面上的排斥力快速增長(zhǎng)，導(dǎo)致中央?yún)^(qū)域總密度降低（圖13（d）～圖13（f）），表現(xiàn)為滑移長(zhǎng)度的增加。

曲線的第3 段具有較好的線性規(guī)律，此時(shí)組分間的溶解密度達(dá)到極小，且趨于0，溶解密度不再隨Gc發(fā)生顯著變化（圖8），界面上的排斥力主要由Gc決 定。顯然排斥力和Gc為 線性關(guān)系，Gc越大，排斥力越大，中心區(qū)域界面上的總密度則越?。▓D13（f）～圖13（h）），滑移長(zhǎng)度越大。同時(shí)值得注意的是，壁面附近原來(lái)的高密度層已經(jīng)轉(zhuǎn)變?yōu)榈兔芏葘?。主要原因是壁面附近組分1 和2 密度分布不均勻。當(dāng)Gc較大時(shí)，在壁面附近組分2的粒子數(shù)量遠(yuǎn)遠(yuǎn)少于組分1。微量的粒子在排斥力作用下更容易聚集在壁面，而聚集在壁面的粒子又會(huì)反過(guò)來(lái)排斥組分2 的粒子，使得壁面附近組分2 密度下降。由于壁面組分2 的密度高于組分1，所以組分2 密度變化影響高于組分1，因此，總密度在壁面附近表現(xiàn)為低密度層（圖15）。

圖15 G c=1.8 ， G ads,1=?Gads,2=0時(shí)密度沿y 方向分布圖Fig. 15 Density value along y axis with parameters Gc=1.8，Gads,1=?Gads,2=0

2.5 剪切率對(duì)滑移長(zhǎng)度的影響

沿用2.4 節(jié)中的計(jì)算模型，固定Gc=1.8，并改變模型頂部的水平速度，結(jié)果如圖16 所示。對(duì)于潤(rùn)滑油的減阻特性來(lái)說(shuō)，其產(chǎn)生的滑移長(zhǎng)度并不顯著依賴于剪切率，這一點(diǎn)與傳統(tǒng)疏水壁面減阻特性[12]相似。這一相似性，主要是因?yàn)楫?dāng)潤(rùn)滑油難溶時(shí)，油液界面可以產(chǎn)生較穩(wěn)定的排斥力，而這與疏水壁面產(chǎn)生的排斥力是相似的。在穩(wěn)定的排斥力作用下，滑移長(zhǎng)度不顯著依賴于剪切率。值得注意的是在文獻(xiàn)[12] 中，水平速度保持不變，但改變豎向距離。與本文相比，雖然方式略有不同，但實(shí)質(zhì)上變化的都是剪切率?u/?y。

圖16 剪切率對(duì)滑移長(zhǎng)度的影響Fig. 16 Influence of shearing rate on slip length

3 結(jié) 論

本文利用格子Boltzmann 方法對(duì)液潤(rùn)表面的滑移現(xiàn)象進(jìn)行了數(shù)值模擬，并研究了潤(rùn)滑油溶解密度、剪切率對(duì)滑移長(zhǎng)度的影響，得到如下結(jié)論：

1） 通過(guò)理論分析Shan-Chen 多組分模型壓力項(xiàng)的表達(dá)式，以及對(duì)拉普拉斯定律的數(shù)值驗(yàn)證，表明Shan-Chen 多組分模型可以有效地反映不同組分間分子作用力，進(jìn)而可以用來(lái)模擬油液分界面上產(chǎn)生的疏水減阻現(xiàn)象。

2） 在其他物理性質(zhì)相似的情況下，潤(rùn)滑油難溶程度越高，油液界面上排斥力越大，誘導(dǎo)的低密度區(qū)密度越低，產(chǎn)生的滑移長(zhǎng)度越大。

3） 在不考慮潤(rùn)滑油流出的情況下，液潤(rùn)表面產(chǎn)生的滑移長(zhǎng)度不顯著依賴于剪切率。