基于格子Boltzmann 方法的冪律流體二維頂蓋驅動流轉捩研究*

張恒 任峰 胡海豹

(西北工業大學航海學院,西安 710072)

研究非牛頓流體轉捩問題,可為調控非牛頓流體動力特性提供理論基礎.相對于牛頓流體轉捩問題,非牛頓流體轉捩研究較少,缺乏轉捩雷諾數精細預報方法.論文以格子Boltzmann 方法為核心求解器,以典型非牛頓流體冪律模型為例,開展了冪律流體二維頂蓋驅動流轉捩模擬,給出剪切變稀和剪切增稠流體的第一轉捩雷諾數,并分析了轉捩雷諾數附近流場時頻域特性及模態分布.結果表明,剪切變稀流體和剪切增稠流體的第一轉捩雷諾數與牛頓流體差異顯著,且在轉捩臨界雷諾數附近監控點處速度分量均呈現周期性變化趨勢.通過對流場速度和渦量的本征正交分解發現,不同類型的流體在轉捩臨界雷諾數附近,前兩階模態均為流場的主模態,能量占比超過95%,且同類型流體不同雷諾數的主模態間具有相似的結構.

1 引 言

頂蓋驅動流是計算流體力學中經典物理模型之一,被廣泛應用于層流流動、流動轉捩等問題的研究[1?4].研究流動轉捩問題,可為界定流場狀態和調控流體動力特性提供理論依據,具有重要學術意義.轉捩問題涉及多種流動狀態,屬復雜流體力學問題,已有大量學者開展過相關研究.例如,Peng 等[5]采用有限差分法研究了頂蓋驅動流的轉捩雷諾數,并給出臨界雷諾數附近流場的時頻域特性.Bruneau 和Saad[6]進一步提高網格分辨率,改進對流項差分格式,佐證了Peng 等[5]、Auteri 等[7]、Sahin 和Owens[8]關于第一轉捩雷諾數在8000 附近的結論.需要指出的是,上述研究均是以牛頓流體為研究對象.

非牛頓流體的剪應力與剪切應變率之間為非線性關系,會表現出剪切增稠[9]、剪切變稀[10]等不同于牛頓流體的運動行為.非牛頓流體廣泛存在于自然界和工業生產中,如化工中的聚乙烯、聚氯乙烯,食品中的淀粉液,人體內的血液等.非牛頓流體種類眾多,其中具有代表性的是冪律流體.冪律流體按照冪律指數主要分為剪切變稀流體、剪切增稠流體.作為特例,牛頓流體可視為指數為1的冪律流體.學者們開展了大量的冪律流體實驗與仿真的研究,如Johnston 等[11]發現相比于牛頓流體,用冪律流體模擬血液流動,可獲更精準的壁面切應力;Zhao 等[12]模擬了冪律流體在微通道中的電滲流動;Hojjat 等[13]通過實驗研究了含納米顆粒流體的剪切變稀行為.但是,關于冪律流體轉捩問題的鮮有研究報道.

與有限差分和有限體積等傳統方法不同,格子Boltzmann 方 法(lattice Boltzmann method,LBM)[14]通過對分布函數進行碰撞和遷移處理,可實現宏觀流場時空演化過程的模擬,具有并行性高,邊界處理簡單等優勢,已在牛頓流體仿真領域取得大量研究成果[15?17].近年來,也有學者將LBM拓展至非牛頓流體模擬,如Boyd 等[18]在LBM 框架下實現了冪律流體中速度導數的快速求解;Wang和Ho[19]在LBM的框架下模擬了剪切變稀的血液流動;Chen 等[20]模擬了冪律流體通過多孔介質的流動.然而,將LBM 應用于冪律流體轉捩問題研究仍有待探索.

基于此,本文以二維頂蓋驅動流為物理模型,以格子Boltzmann 方法為核心求解器,采用非牛頓流體冪律模型,開展了非牛頓流體二維頂蓋驅動流轉捩模擬,并借助快速傅里葉變換、本征正交分解(proper orthogonal decomposition,POD)等方法,分析了轉捩點附近的流場時頻域特性及模態分布.

2 數值方法

2.1 物理模型

頂蓋驅動流模型如圖1 所示,上頂板為運動壁面,速度為U0,其余3 個壁面均為靜止無滑移壁面.模型水平和豎直方向上長度分別為Lx和Ly,且Lx=Ly=L0,L0為特征長度.圖1 中,坐標原點O位于底面與左壁面交點處.為定量表征流場狀態與雷諾數(Reynolds number,Re)之間的關系,選取監控點坐標為(0.6L0,0.6L0).

圖1 算例模型示意圖Fig.1.Schematic diagram of calculation model.

冪律流體運動黏性系數的計算公式如下:

式中,S為流場中格點處應變率張量,γ˙ 為剪切率,n為冪律指數,ν為運動黏性系數,m為冪律系數(與Re有關的常量).當n=1 時為牛頓流體;n<1為剪切變稀流體;n>1 時為剪切增稠流體.

雷諾數的定義為[21]

2.2 數值方法

采用LBM 多松弛時間模型[22],并使用基于統一計算設備構架的顯卡并行加速技術[23]來提高計算速度.其中,多松弛時間模型表達式為

式中,f(x,t) 為t時刻流場中坐標位置為x處的分布函數,ei為離散速度,δt為時間步長,M為從速度空間向矩空間進行轉換的矩陣,為平衡態分布函數.

針對文中研究的二維頂蓋驅動流,選用經典D2Q9 模型[24],離散速度ei為

LBM 平衡態分布函數采用不可壓模型[25],其表達式為

式中,ρ0為不可壓流體的密度;p為壓強;u為格點處的速度;權系數wi的取值為w0=4/9 ,w1?4=1/9,w5?8=1/36 ;cs為格子聲速.松弛參數對角矩陣選為Λf=[1.0,1.1,1.0,1.0,1.2,1.0,1.2,1/τ,1/τ],松弛時間τ與運動黏性關系為ν=(τ ?0.5).

流場中宏觀速度與壓強的求解公式為

為保證流場時空分辨率,本文最終選用的網格分辨率為512 × 512,計算時間總長2000T0(T0為參考時間),即1600 萬個格子時間步長.靜壁邊界采用半步長反彈格式[26],動壁采用動量補償格式[27].為防止運動黏性系數過大或者過小造成計算結果發散,這里設置運動黏性系數上限和下限[28]分別為 2 0ν0,ν0/20 .不失一般性,文中取n=0.75,1.25 分別代表剪切變稀流體和剪切增稠流體.具體模擬流程如圖2 所示.

圖2 LBM 與冪律模型耦合計算流程示意圖Fig.2.Diagram of coupling calculation process between LBM and power law model.

2.3 POD 方法

POD[29,30]是一種廣泛應用于流場降階模型中的方法,可用于流場特征信息的提取.以速度POD為例,已知N個時刻流場速度信息u(x,ti),i=1—N,則有

另 外,Σ=diag(σ1,σ2,···σi,···)中各σi模長的平方可用來表征對應模態的“能量”.按照模態能量的大小進行排序即可得到流場的主要模態.

2.4 第一轉捩雷諾數計算

在頂蓋驅動流中,隨著Re增加,呈現出的流動狀態往往是從層流逐步過渡到湍流,并存在Hopf 分岔現象.通過在流場的特定位置處設置監控點的方法[5,31],可以研究Hopf 分岔現象.開始產生分岔現象時的雷諾數,即第一轉捩雷諾數可以通過如下的方式計算:

其中 max[u(t)] 與 min[u(t)] 分別為流場充分發展后,監控點處水平速度u的最大值與最小值.鑒于與Re成正比關系,因此可以通過計算不同雷諾數下的,使用最小二乘法線性擬合出Re-圖.擬合所得的直線與=0的交點,即流場從穩定向周期態轉變的第一轉捩臨界雷諾數.

3 計算結果與分析

3.1 算例驗證

為了驗證計算的準確性,這里首先給出了牛頓流體頂蓋驅動流模擬結果.圖3 為Re=5000 時,方腔中心線上的速度分布與Peng 等[5]結果的對比,其中X-v曲線顯示的是豎直方向上方腔的中心線上速度分量分布,即縱坐標y/L=0.5 上的X和v之間的對應關系;Y-u曲線顯示的是水平方向上的中心線上速度分量分布,即橫坐標x/L=0.5上的Y和u之間的關系.Peng 等[5]在牛頓流體頂蓋驅動流轉捩問題研究中,使用的是限差分法,本文的計算結果與Peng 等[5]對比良好.值得注意的是Peng 等[5]求解對流項使用的是七階迎風格式,擴散項使用的是六階中心差分格式.

圖3 頂蓋驅動流在中心線處的速度分布Fig.3.Velocity distribution at center line of lid-driven flow.

3.2 第一轉捩雷諾數

按照文中2.4 節的方法得到牛頓流體、剪切變稀流體和剪切增稠流體的第一轉捩雷諾數依次為7835,5496 和11546,結果如圖4 所示.

圖4 轉捩臨界雷諾數Fig.4.Transition critical Reynolds number.

將本文計算的牛頓流體轉捩雷諾數與文獻的結果進行比對,如表1 所列,可以看出,本文計算的牛頓流體轉捩臨界雷諾數相對于Peng 等[5]和Poliashenko 與Aidu[32]的計算結果而言較大,而相對于Cazemier 等[31]和Bruneau 與Saad[6]的計算結果而言較小,但是從總的相對誤差的值來看,誤差上限不超過5%,可以看出本文計算可靠性較高.

表1 牛頓流體轉捩臨界雷諾數計算結果的對比Table 1. Comparison of calculation results of critical Reynolds number for Newtonian fluid transition.

相比于牛頓流體,剪切變稀流體的轉捩雷諾數較低,即剪切變稀流體更容易從定常狀態向周期態轉變.而由剪切變稀流體本身的定義可知,剪切速率越大,其黏性會越低,流動狀態也越容易發生轉變,這與計算結果表現出的性質是相符合的.與牛頓流體相比,剪切增稠流體的轉捩臨界雷諾數較高,即剪切增稠流體更不易從定常狀態向周期態轉變.而根據剪切增稠流體本身的定義可知,剪切速率越大,其黏性會越大,流動狀態也越不易發生轉變.

3.3 速度監控點時頻域特性

選取三類流體在轉捩雷諾數附近的流場,進行時頻域分析.圖5 給出了三類流體在轉捩雷諾數附近,速度監控點處的時頻域特性圖.其中圖5(a)—(c)給出了三類流體在轉捩前后,水平方向速度隨時間的變化,可以看出,三類流體在轉捩發生前的定常流動狀態下,速度為一水平直線,在轉捩雷諾數附近,速度均呈現周期性變化的特點,在轉捩后湍流狀態下,速度有明顯的大幅度波動.對水平速度作快速傅里葉變換得圖5(d)—(f),在發生轉捩前,由于其直流部分即f=0 為其主要分量,在頻譜圖上顯示為一點,在轉捩雷諾數附近,水平速度均存在一個主頻f1與諧頻f2,且滿足f2=2f1,但不同類型流體的f1值并不相同,在轉捩后湍流狀態下,出現了明顯的多頻現象.圖5(g)—(i)為三類流體在轉捩雷諾數附近,水平速度u與豎直方向速度v之間形成的速度相圖.可以發現不同類型的流體的速度相圖的變化范圍并不相同,但均為單一光滑封閉曲線.

圖5 三類流體在速度監控點處的時頻域特性 (a),(d) n=1,Re=7500,8500,16000;(g) n=1,Re=8500;(b),(e) n=0.75,Re=5500,6500,10000;(h) n=0.75,Re=6500;(c),(f) n=1.25,Re=9500,13000,20000;(i) n=1.25,Re=13000Fig.5.Time-frequency spectrum characteristics at the velocity monitoring point for three types of fluids:(a),(d) n=1,Re=7500,8500,16000;(g) n=1,Re=8500;(b),(e) n=0.75,Re=5500,6500,10000;(h) n=0.75,Re=6500;(c),(f) n=1.25,Re=9500,13000,20000;(i) n=1.25,Re=13000.

圖6 為監控點處牛頓流體、剪切變稀與剪切增稠流體在不同雷諾數情況下的速度相圖,各封閉曲線中心點用虛線進行連接.圖6(a)為牛頓流體在不同雷諾數下的速度相圖,其中Re=7500,8000時,速度波動極小,在相圖上呈現的是1 個點,而在Re=8500—9500 時,相圖為光滑封閉曲線,曲線整體向右上方移動且包圍的面積在不斷增加.圖6(b)為剪切變稀流體在不同雷諾數下的速度相圖,在Re=5500 時,速度相圖為1 個點,而Re=5700—6000 時,速度相圖為封閉曲線,隨著Re增加向上方移動,且包圍的面積在不斷增加.圖6(c)為剪切增稠流體在不同雷諾數下的速度相圖,其中Re=11500,12000,12500 時,速度相圖為單一點,而Re=13000—15000 時,速度相圖為光滑封閉曲線.速度相圖整體變化趨勢為先向右再向左上方移動,向右移動的過程中,豎直方向上速度v幾乎沒有變化,而向左上方移動的過程中,封閉曲線包圍的面積逐漸增大,有較為明顯的變化規律.

圖6 在監控點處的速度相圖 (a) 牛頓流體;(b) 剪切變稀流體;(c) 剪切增稠流體Fig.6.Velocity phase diagrams at the monitoring point:(a) Newtonian fluid;(b) shear-thinning fluid;(c) shear-thickening fluid.

3.4 流場模態分析

進一步使用POD 對牛頓流體、剪切變稀和剪切增稠流體的流場速度、渦量和黏性進行處理,結果如圖7 和圖8 所示.圖7(a)—(c)分別為牛頓流體、剪切變稀流體和剪切增稠流體的水平速度POD能量占比圖.從圖7 可以發現,不同類型的流體在第一臨界轉捩雷諾數附近且發生轉捩后,均呈現前兩階模態能量占主導地位的特點,能量占比超95%,且前兩階模態能量占比接近.在圖7(a)Re=10500 時,出現了多模態的結果,這里的流場可能是進入Hopf 第二分岔第二模式中[5].在圖7(c)Re=11500,12500 時,一階模態能量占比大幅度高于二階模態的能量占比,但結合圖6(c)可以發現,此時流場速度波動量較小幾乎為定常流場,一階模態與二階模態相對于平均流場幾乎可以忽略.

圖7 速度u的各階模態的能量占比 (a) 牛頓流體;(b) 剪切變稀流體;(c) 剪切增稠流體Fig.7.Energy share of each order of mode for velocity u:(a) Newtonian fluid;(b) shear thinning-fluid;(c) shear-thickening fluid.

圖8(a)—(c)分別為牛頓流體、剪切變稀流體和剪切增稠流體的渦量POD 能量占比圖.渦量POD 能量占比與水平速度POD 能量占比結果相比,同種流體在同一雷諾數下,渦量POD 一階模態的占比略高于水平速度POD 能量占比,渦量POD 中前兩階模仍占主導地位.

圖8 渦量的各階模態的能量占比 (a) 牛頓流體;(b) 剪切變稀流體;(c) 剪切增稠流體Fig.8.Vortex energy share of each order of mode:(a) Newtonian fluid;(b) shear-thinning fluid;(c) shear-thickening fluid.

圖9 為牛頓流體在不同雷諾數下各階模態的流場模態的速度場云圖,圖9(a)和圖9(d)分別為Re=8500 與Re=9000 時平均水平速度場云圖,圖9(b)和圖9(e) 分別為Re=8500 與Re=9000時水平速度一階模態云圖,圖9(c)和圖9(f) 分別為Re=8500 與Re=9000 時水平速度二階模態云圖.從速度平均場來看,Re=8500 與Re=9000的結果相近;從一階模態的云圖來看,峰值區域均靠近壁面,且峰值區域的形狀相似;從二階模態的云圖來看,二階模態的結果相近,且峰值主要分布在壁面處.

圖9 牛頓流體的水平速度的各階模態圖 (a),(b),(c) Re=8500 時,平均場、一階模態與二階模態;(d),(e),(f) Re=9000 時,平均場、一階模態與二階模態Fig.9.Modal diagrams of horizontal velocity of Newtonian fluid:(a),(b) and(c) The mean field,the first and second modes when Re=8500;(d),(e),(f) mean field,first-order mode and second-order mode when Re=9000.

圖10 為剪切變稀流體在不同Re下各階模態的流場模態的速度場云圖,圖10(a)和圖10(d)分別為Re=5700 與Re=6000 時平均速度場云圖,圖10(b)和圖10(e)分別為Re=5700 與Re=6000時速度一階模態云圖,圖10(c)和圖10(f)分別為Re=5700 與Re=6000 時速度二階模態云圖.從平均場的結果來看,Re=5700 與Re=6000 結果相近;從一階和二階模態的結果來看,模態峰值分布區域集中在壁面附近,峰值分布區域的形狀相似,但峰值分布區的模態值相反.

圖11 為剪切增稠流體在不同Re下各階模態的流場模態的速度場云圖,圖11(a)和圖11(d)分別為Re=13000 與Re=13500 時平均速度場云圖,圖11(b)和圖11(e)分別為Re=13000 與Re=13500 時速度一階模態云圖,圖11(c)和圖11(f)分別為Re=13000 與Re=13500 時速度二階模態云圖.從平均場的結果來看,Re=13000 與Re=13500 結果相近;從一階和二階模態的結果來看,峰值分布區域集中在壁面附近,且峰值分布區域的形狀相似,這與牛頓流體和剪切變稀流體的結果是類似的.剪切增稠流體一階和二階模態峰值分布區的模態值相反.

圖11 剪切增稠流體的水平速度的各階模態圖 (a),(b),(c) Re=13000 時,平均場、一階模態與二階模態;(d),(e),(f) Re=13500 時,平均場、一階模態與二階模態Fig.11.Modal diagrams of horizontal velocity of shear-thickening fluid:(a),(b),(c) The mean field,the first and second modes when Re=13000;(d),(e),(f) mean field,first-order mode and second-order mode when Re=13500.

4 結 論

論文針對典型的非牛頓流體—冪律流體,利用格子Boltzmann的計算框架,以頂蓋驅動流為物理模型,研究得到冪律指數在0.75 和1.25 時的第一臨界轉捩Re分別為5496 和11546.在轉捩Re附近且發生轉捩后,流場監控點的速度隨時間呈周期性變化,存在一個主頻與諧頻.與轉捩前速度相圖為單一點不同,轉捩后速度相圖為光滑封閉曲線,且隨著Re增加,曲線中心點朝一定的方向移動且包絡面積不斷增加,但不同流體中心點移動的方向不相同.流場的速度、渦量POD 分解的結果均表明,流場在轉捩臨界雷諾數附近,其前兩階模態為其主要模態,能量占比超95%.從模態云圖的結果看,一階模態與二階模態速度峰值分布區域主要在壁面附近,且同種流體不同Re下的平均場結果相近,一階模態與二階模態的峰值分布區域形狀相似.