喬振軒，王樹青，王淑彥

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

后臺階流動是流體力學中一種非常典型的分離流動，如壩面溢流、管道中的突擴、繞飛行器和建筑物氣流流動等。在這種流動中，流體流過臺階后會發(fā)生分離，在臺階角處形成回流渦，待繞過回流渦后再與壁面接觸，然后繼續(xù)沿著壁面流動。在分離流動研究中，最常用的簡化模型就是后臺階流動，因為其幾何形狀簡單，具有分離流的流動特性，所以被廣泛研究討論，在工程設備上有很多應用。

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

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

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

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

1.1 格子Boltzmann方法（LBM）

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

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

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

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

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

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

圖1 D3Q19模型Fig.1 D3Q19 model

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

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

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

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

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

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

連續(xù)方程：

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

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

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

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

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

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

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

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

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

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

圖2 三維后臺階流動模型Fig.2 The three-dimensional model of back-facing step flow

表1 后臺階流動基礎模型運算參數(shù)設置Table 1 The basic parameters of back-facing step flow simulation

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

2 結(jié)果與討論

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

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

2.2 后臺階流動回流區(qū)研究

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

式中，H為臺階高度；S為回流區(qū)域長度。

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

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

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

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

2.3 雷諾數(shù)對后臺階流動的影響

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

圖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ù)下后臺階流體流動特性的研究

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

基于LBM對雷諾數(shù)對后臺階流動的影響進行了研究，所的結(jié)論如下：

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

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

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