基于格子Boltzmann方法后臺(tái)階流動(dòng)特性的研究

喬振軒，王樹青，王淑彥

（1.東北石油大學(xué)石油工程學(xué)院，黑龍江大慶163318；2.大慶油田有限責(zé)任公司第四采油廠，黑龍江大慶163000）

后臺(tái)階流動(dòng)是流體力學(xué)中一種非常典型的分離流動(dòng)，如壩面溢流、管道中的突擴(kuò)、繞飛行器和建筑物氣流流動(dòng)等。在這種流動(dòng)中，流體流過臺(tái)階后會(huì)發(fā)生分離，在臺(tái)階角處形成回流渦，待繞過回流渦后再與壁面接觸，然后繼續(xù)沿著壁面流動(dòng)。在分離流動(dòng)研究中，最常用的簡(jiǎn)化模型就是后臺(tái)階流動(dòng)，因?yàn)槠鋷缀涡螤詈?jiǎn)單，具有分離流的流動(dòng)特性，所以被廣泛研究討論，在工程設(shè)備上有很多應(yīng)用。

雷諾數(shù)是衡量流體流動(dòng)狀態(tài)的主要無因次數(shù)，在后臺(tái)階流動(dòng)中雷諾數(shù)主要用來表征臺(tái)階后產(chǎn)生的渦旋狀況。H.Nowruzi等[1]采用能量梯度法研究了不同雷諾數(shù)和膨脹比下二維層流后臺(tái)階（BFS）流動(dòng)的不穩(wěn)定性，且通過速度、渦度和能量梯度函數(shù)的分布來研究流動(dòng)的不穩(wěn)定性。V.I.Terekhov等[2]通過PIV實(shí)驗(yàn)研究了臺(tái)階上肋片對(duì)后臺(tái)階流動(dòng)的影響，獲得了單個(gè)肋片對(duì)平均流速和湍流波動(dòng)的影響規(guī)律，分析了肋片的位置和高度對(duì)再循環(huán)區(qū)域的影響。

格子Boltzmann方法（LBM）作為一種復(fù)雜流動(dòng)的數(shù)值模擬方法，具有形式簡(jiǎn)單、復(fù)雜邊界易于處理等優(yōu)點(diǎn)，受到了國(guó)內(nèi)外學(xué)者的廣泛關(guān)注[3-5]。大渦模擬(LES)方法是對(duì)湍流脈動(dòng)的一種空間平均，通過某種濾波函數(shù)將大尺度的渦和小尺度的渦分離開，大尺度的渦直接模擬，小尺度的渦用模型來封閉，因而它是求解高雷諾數(shù)湍流流動(dòng)的主要方法。將LBM與LES結(jié)合起來的LBM-LES的方法提高了計(jì)算的準(zhǔn)確度和速度，目前常被用于模擬復(fù)雜流場(chǎng)的高雷諾數(shù)流動(dòng)。O.Malaspinas等[4]提出了基于LBM的大渦模擬的原始壁模型，它依賴于第一個(gè)離壁單元內(nèi)的速度分布或湍流邊界層方程的解，該模型具有較好的穩(wěn)定性和準(zhǔn)確性。H.Si等[5]采用LBM-LES的方法模擬高雷諾數(shù)下的流體流動(dòng)，結(jié)果表明LBM結(jié)合大渦模擬模型可以有效地模擬高雷諾數(shù)流動(dòng)。A.Pradhan等[6]采用LBM-LES模擬通道中的湍流，研究結(jié)果表明，這種方法適合復(fù)雜幾何結(jié)構(gòu)的流場(chǎng)并且計(jì)算穩(wěn)定。K.N.Premnath等[7]采用LBM-LES模擬了后臺(tái)階流動(dòng)，對(duì)再循環(huán)和再附著區(qū)域內(nèi)的湍流信息進(jìn)行統(tǒng)計(jì)，與直接模擬結(jié)果和實(shí)驗(yàn)結(jié)果吻合。A.F.A.Gawad等[8]采用LBMLES研究了在二維通道內(nèi)障礙物后面的湍流特性和渦旋脫落機(jī)理，結(jié)果表明，對(duì)于一個(gè)和兩個(gè)方形障礙物的兩種情況下渦旋脫落具有不同的特征。

本文采用LBM-LES研究高雷諾數(shù)下的后臺(tái)階流動(dòng)，得到后臺(tái)階流動(dòng)的演變過程，并且與文獻(xiàn)中的數(shù)據(jù)作對(duì)比，以此來驗(yàn)證模擬結(jié)果的正確性。

1 數(shù)學(xué)模型及數(shù)值方法

1.1 格子Boltzmann方法（LBM）

玻爾茲曼方程是描述粒子速度分布函數(shù)f時(shí)空變化的守恒方程：

式中，ξ表示粒子的速度，Ω(f)是表示碰撞的影響的碰撞算子。對(duì)該方程進(jìn)行離散可以得到格子Boltzmann方程（LBE）：

由于碰撞算子設(shè)計(jì)極為復(fù)雜的非線性積分，Boltzmann方程在應(yīng)用中受到的限制很大。人們?yōu)榇颂岢隽撕芏喾N簡(jiǎn)化的碰撞算子或碰撞模型。BGK模型[9]是最簡(jiǎn)單也最常用的碰撞算子模型，該模型的碰撞算子表示為：

引入BGK模型后，格子Boltzmann方程可近似為：

式中，fi(x，t)是粒子在t時(shí)刻x處離散速度為ci的分布函數(shù)，τ為松弛時(shí)間，Δt為時(shí)間步長(zhǎng)，feqi為對(duì)應(yīng)的平衡態(tài)分布函數(shù)。

BGK模型中，由Y.H.Qian等[10]提出的DnQb模型（n為空間維數(shù)，b為離散速度數(shù)）最有代表性。本文采用D3Q19模型，在三維空間內(nèi)有19個(gè)離散速度，如圖1所示。

圖1 D3Q19模型Fig.1 D3Q19 model

在D3Q19模型中，平衡態(tài)分布函數(shù)可表示為：

式中，ωi為權(quán)系數(shù)，cs為聲速。.

可以根據(jù)式（6）來計(jì)算宏觀密度和宏觀速度：

1.2 基于格子Boltzmann方法的大渦模擬（LBM-LES）

大渦模擬的基本方法是通過直接求解N-S方程來模擬大尺度的運(yùn)動(dòng)，通過亞格子格式模擬小尺度運(yùn)動(dòng)。大渦模擬的第一步就是進(jìn)行濾波操作：

式中，w是速度、密度等流場(chǎng)宏觀量，G為空間濾波函數(shù)。假設(shè)過濾過程和求導(dǎo)過程可以交換，將濾波運(yùn)算用于N-S方程，則可得到相應(yīng)的濾波控制方程。

連續(xù)方程：

目前，較常用的湍流亞格子尺度模型是Smagorinsky 模型[11]：

式中，δij為 Krone-cker張量；νt為渦黏性系數(shù)；||為大尺度應(yīng)變張量Sij的模；C為Smagorinsky常數(shù)，C＞0；Δ為濾波寬度。

本文采用基于格子Boltzmann方法的大渦模擬（LBM-LES），仿照Smagorinsky亞格子模型的思路，保持模型中的平衡分布函數(shù)不變，以及原模型中的馳豫時(shí)間與黏性的關(guān)系，即：

式（13）中，總黏性系數(shù)νtot分為物理動(dòng)力黏性ν0和渦黏性νt兩部分：

根據(jù)過濾節(jié)點(diǎn)的分布函數(shù)的非平衡部分求得[12]：

聯(lián)立上述方程，即可求得|-S|。

至此，D3Q19的演化方程變?yōu)椋?/p>

將 τtot、||帶入式（15），即可求得 BGK 模型中的總碰撞算子，代入演化方程式（19），即可求解。

1.3 計(jì)算參數(shù)和邊界條件

后臺(tái)階流動(dòng)數(shù)值模擬的計(jì)算域幾何結(jié)構(gòu)如圖2所示。流體沿X軸方向從通道左側(cè)流入，速度由雷諾數(shù)計(jì)算得出，通過改變?nèi)肟谒俣雀淖兝字Z數(shù)，從通道右側(cè)出口流出，出口為壓力出口，壓力為101 kPa，上下壁面和臺(tái)階壁面為標(biāo)準(zhǔn)反彈格式，而通道的前后壁面為非平衡外推邊界條件。臺(tái)階位于通道左側(cè)，長(zhǎng)度為L(zhǎng)，寬度為D，高度為H。通道長(zhǎng)度為L(zhǎng)c，寬度為Dc，高度為Hc，具體幾何尺寸和流體物性參數(shù)如表1所示。

圖2 三維后臺(tái)階流動(dòng)模型Fig.2 The three-dimensional model of back-facing step flow

表1 后臺(tái)階流動(dòng)基礎(chǔ)模型運(yùn)算參數(shù)設(shè)置Table 1 The basic parameters of back-facing step flow simulation

由于網(wǎng)格的分辨率同時(shí)影響計(jì)算的精度和效率，所以需要驗(yàn)證網(wǎng)格的無關(guān)性。分別采用100×20×20、150×30×30、200×40×40和 250×50×50四種網(wǎng)格劃分情況模擬后臺(tái)階流動(dòng)，得出四種計(jì)算結(jié)果的漩渦大小差異很小，但較密的網(wǎng)格可以預(yù)測(cè)出詳細(xì)的漩渦內(nèi)部結(jié)構(gòu)，綜合考慮計(jì)算精度與成本，選擇200×40×40這種網(wǎng)格進(jìn)行下列的模擬計(jì)算。

2 結(jié)果與討論

2.1 瞬時(shí)流場(chǎng)渦結(jié)構(gòu)的演變過程

圖3表示不同時(shí)刻臺(tái)階后部流體的流線分布，即描述流場(chǎng)的渦結(jié)構(gòu)的瞬時(shí)變化。從圖3中可以看出，隨著時(shí)間逐漸推移，流體流過臺(tái)階后首先經(jīng)歷了分離層的分離，在臺(tái)階后形成了與臺(tái)階尺寸相近的沿著流動(dòng)方向的旋渦。旋渦內(nèi)的流體進(jìn)行的是有旋流動(dòng)，旋轉(zhuǎn)方向?yàn)轫槙r(shí)針。隨著時(shí)間的推移，旋渦狀結(jié)構(gòu)的尺寸逐漸增大，緩慢進(jìn)入到充分發(fā)展階段，在臺(tái)階下游形成穩(wěn)定的回流區(qū)域，這是后臺(tái)階流動(dòng)過程中表現(xiàn)出的主要特性。此外，從圖3中還可以看出，有少量流體沿展向方向流動(dòng)，說明對(duì)后臺(tái)階流動(dòng)進(jìn)行三維的數(shù)值模擬是必要的。

2.2 后臺(tái)階流動(dòng)回流區(qū)研究

在研究后臺(tái)階流體流動(dòng)時(shí)，通常用再附著點(diǎn)與分離點(diǎn)之間的距離表示回流區(qū)域的長(zhǎng)度，而回流區(qū)域長(zhǎng)度與臺(tái)階高度的比值K是對(duì)后臺(tái)階流體流動(dòng)特性進(jìn)行定量分析的重要參數(shù)。K值計(jì)算公式如下：

式中，H為臺(tái)階高度；S為回流區(qū)域長(zhǎng)度。

回流區(qū)域的長(zhǎng)度可以通過對(duì)流體的軸向速度進(jìn)行分析而取得。首先，分離點(diǎn)的位置是固定的，位于臺(tái)階后壁面上邊緣，其次，再附著點(diǎn)位于通道下壁面，假設(shè)當(dāng)通道下壁面附近流體軸向速度不再為負(fù)時(shí)的位置為（x1，y1，z1），那么就可以根據(jù)該點(diǎn)來確定再附著點(diǎn)的位置。

圖4表示流動(dòng)穩(wěn)定后通道中沿流動(dòng)方向，不同位置處流體的軸向速度沿徑向的分布曲線。從圖4中可以看出，臺(tái)階后壁面附近（X/l=22）流體的軸向速度均為負(fù)值，但是數(shù)值比較小，整體上接近于0，此處軸向速度沿徑向的分布沒有明顯的規(guī)律；靠近漩渦區(qū)域內(nèi)部（X/l=50，X/l=60，X/l=80）流體的軸向速度沿徑向的分布規(guī)律性較為明顯，沿著流動(dòng)方向，回流區(qū)減弱，在旋渦中心位置出現(xiàn)最大值，而后逐漸減小；當(dāng)X/l=94時(shí)，流體的速度值全部為正。通過計(jì)算得出中軸面上流動(dòng)穩(wěn)定后K值約為3.38。

圖3 瞬時(shí)流場(chǎng)渦結(jié)構(gòu)的演變Fig.3 Evolution of vortex structure in transient flow field

圖5 表示通道內(nèi)不同軸向截面旋渦的K值變化。通道中心區(qū)域臺(tái)階后回流區(qū)域長(zhǎng)度與K值較為接近。隨著位置偏離，臺(tái)階后回流區(qū)域長(zhǎng)度逐漸變短，K值隨之減小。鑒于中軸面上漩渦的K值最大，因此在進(jìn)行對(duì)比時(shí)均選取中軸面上的K值代表漩渦的K值。

2.3 雷諾數(shù)對(duì)后臺(tái)階流動(dòng)的影響

流體的流動(dòng)參數(shù)對(duì)后臺(tái)階流體流動(dòng)特性有不可忽視的影響，其中雷諾數(shù)是十分重要的影響因素之一。本節(jié)將分別針對(duì)低雷諾數(shù)和高雷諾數(shù)兩種情況來研究不同雷諾數(shù)下后臺(tái)階流動(dòng)特性。

圖4 流體軸向速度沿徑向分布Fig.4 The radial distribution of the fluid axial velocity

圖5 不同軸向截面處K值Fig.5 Variation of K value versus axial sections

2.3.1 低雷諾數(shù)下后臺(tái)階流體流動(dòng)特性的研究

圖6表示低雷諾數(shù)情況下，不同雷諾數(shù)時(shí)通道內(nèi)流體的流線圖。從圖6中可以看出，在低雷諾數(shù)情況下，雷諾數(shù)對(duì)通道內(nèi)再附著點(diǎn)的位置和臺(tái)階后回流區(qū)域的長(zhǎng)度具有較為重要的影響。隨著雷諾數(shù)逐漸增大，通道中臺(tái)階后的旋渦狀結(jié)構(gòu)更加明顯，旋渦狀結(jié)構(gòu)內(nèi)流體的運(yùn)動(dòng)形態(tài)也更加清晰。但是隨著雷諾數(shù)增加，旋渦狀結(jié)構(gòu)內(nèi)部的流體受黏滯力影響下降，與通道上部流體之間的速度梯度隨之變小。臺(tái)階后再附著點(diǎn)與分離點(diǎn)的距離隨著雷諾數(shù)的增加而變大。

圖6 不同低雷諾數(shù)下旋渦區(qū)域流線分布Fig.6 Streamlines distribution of vortex versus low Reynolds numbers

圖7 表示回流區(qū)長(zhǎng)度隨雷諾數(shù)的變化曲線。從圖7中可以看出，隨著雷諾數(shù)的增大，慣性力與黏性力的比值逐漸增加，慣性力對(duì)流體的流動(dòng)的控制相對(duì)變大，黏性力對(duì)流體的影響相對(duì)變小，由黏性力控制的回流過程逐漸變得困難，導(dǎo)致回流區(qū)域的長(zhǎng)度變大。

圖7 回流區(qū)長(zhǎng)度隨低雷諾數(shù)的變化曲線Fig.7 Variation of the length of the recirculation zone versus low Reynolds numbers

圖8 表示K值隨Re的變化曲線。從圖8中可以看出，K值隨雷諾數(shù)的增加而增加，但增長(zhǎng)幅度逐漸變小。說明了在低雷諾數(shù)時(shí)慣性力的增加對(duì)回流的影響更大，隨著雷諾數(shù)的增加慣性力的影響逐漸減弱。S.N.Sinha等[13]在風(fēng)洞中布置了一個(gè)后臺(tái)階裝置，模擬三維后臺(tái)階流動(dòng)，B.F.Armaly等[14]的實(shí)驗(yàn)則是在二維的裝置中完成的。將計(jì)算結(jié)果與Sinha、Armaly的模擬結(jié)果進(jìn)行對(duì)比，可以發(fā)現(xiàn)計(jì)算結(jié)果與Sinha的結(jié)果基本吻合，與Armaly的結(jié)果有一定差異，這是因?yàn)锳rmaly的模擬結(jié)果是通過二維計(jì)算得到的，而本文的計(jì)算為三維模擬，同時(shí)也說明了雷諾數(shù)的增加，三維展向流動(dòng)增強(qiáng)。

圖8 不同雷諾數(shù)時(shí)K值變化Fig.8 Variation of K value versus low Reynolds numbers

2.3.2 高雷諾數(shù)下后臺(tái)階流體流動(dòng)特性 圖9為不同雷諾數(shù)時(shí)后臺(tái)階流動(dòng)的流線圖。從圖9中可以看出，隨著雷諾數(shù)的增加，流體的非均勻流動(dòng)特性增強(qiáng)。在高雷諾數(shù)下，首先，對(duì)于處在通道上部的流體而言，由于通道內(nèi)流體處于充分發(fā)展的湍流狀態(tài)，這部分流體的流動(dòng)也變得不再穩(wěn)定，表現(xiàn)出了一定的無序性。其次，對(duì)于處在臺(tái)階后部的流體而言，臺(tái)階后依然會(huì)形成一個(gè)充分發(fā)展的漩渦狀結(jié)構(gòu)。漩渦狀結(jié)構(gòu)隨雷諾數(shù)的增加而逐漸變得明顯，之后又變得紊亂，這是由于超大的慣性力，造成了流動(dòng)的無序性。再附著點(diǎn)與分離點(diǎn)之間的距離隨著雷諾數(shù)的增大而增長(zhǎng)，增長(zhǎng)幅度逐漸降低，在Re=8 000時(shí)達(dá)到峰值，之后逐漸減小，最后趨于一個(gè)定值，說明了慣性力對(duì)通道內(nèi)流動(dòng)的影響具有一定的局限性。

Armaly等實(shí)驗(yàn)研究了雷諾數(shù)為70～8 000的二維后臺(tái)階流動(dòng)，采用激光多普勒獲得測(cè)量結(jié)果。P.M.Nadge等[15]在兩種不同ER比的風(fēng)洞中進(jìn)行了二維的后臺(tái)階流動(dòng)實(shí)驗(yàn)，并對(duì)不同雷諾數(shù)下的實(shí)驗(yàn)結(jié)果進(jìn)行了分析。將計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，結(jié)果見圖10。由圖10可知，計(jì)算結(jié)果與實(shí)驗(yàn)值略有差異，模擬結(jié)果與Armaly等的實(shí)驗(yàn)數(shù)據(jù)基本一致，與Nadge等的實(shí)驗(yàn)數(shù)據(jù)的誤差為7.69%。分析其原因主要是：首先，Nadge等的實(shí)驗(yàn)是在二維的情況下進(jìn)行的，本文的數(shù)據(jù)是在三維的模擬中獲得的，在高雷諾數(shù)下的三維湍流模擬中，存在著很強(qiáng)的展向流動(dòng)，會(huì)對(duì)回流區(qū)域的形態(tài)產(chǎn)生影響。

圖9 不同雷諾數(shù)時(shí)漩渦區(qū)域流線分析Fig.9 Streamlines distribution of vortex versus high Reynolds numbers

其次，ER值和臺(tái)階高度均對(duì)流場(chǎng)有一定影響，Nadge等的實(shí)驗(yàn)在這兩點(diǎn)都有一定差異，這一點(diǎn)也會(huì)在一定程度上導(dǎo)致結(jié)果的差異。最后，本文采用LBM-LES的方法模擬后臺(tái)階流動(dòng)在高雷諾數(shù)下的運(yùn)動(dòng)規(guī)律，模擬過程中使用的Smagorinsky常數(shù)C的取值，對(duì)湍流的模擬會(huì)產(chǎn)生一定程度的影響，但在本文的計(jì)算中，取其均值為0.1，因而會(huì)與實(shí)驗(yàn)數(shù)據(jù)產(chǎn)生誤差。通過對(duì)比，可以說明LBM-LES模型可以用來模擬高雷諾數(shù)下后臺(tái)階流動(dòng)。

圖10 不同高雷諾數(shù)時(shí)K值變化Fig.10 Variation of K value versus high Reynolds numbers

3 結(jié) 論

基于LBM對(duì)雷諾數(shù)對(duì)后臺(tái)階流動(dòng)的影響進(jìn)行了研究，所的結(jié)論如下：

（1）在低雷諾數(shù)情況下，臺(tái)階后再附著點(diǎn)與分離點(diǎn)的距離隨著雷諾數(shù)的增加而變大；K值隨雷諾數(shù)的增加而增加，但增長(zhǎng)幅度逐漸變小。

（2）在高雷諾數(shù)情況下，臺(tái)階后再附著點(diǎn)與分離點(diǎn)的距離隨著雷諾數(shù)的增加而增加，增長(zhǎng)幅度逐漸降低，當(dāng)Re=8 000時(shí)，再附著點(diǎn)與分離點(diǎn)之間的距離趨于穩(wěn)定，并達(dá)到峰值；隨著雷諾數(shù)的繼續(xù)增大，這個(gè)距離開始減小，最終趨于一個(gè)定值。

（3）對(duì)于高雷諾數(shù)下的后臺(tái)階流動(dòng)，本文采用了LES-LBM的方法進(jìn)行模擬，與前人的實(shí)驗(yàn)數(shù)據(jù)吻合得很好，表明利用LBM-LES模型對(duì)高雷諾數(shù)后臺(tái)階流動(dòng)進(jìn)行模擬具有較高的精確性。