基于格子Boltzmann方法的平板射流大渦模擬

上官燕琴, 王 嫻通訊作者：王嫻（1977－）,女,漢族,吉林,副教授,博士生導師,主要研究方向：基于CPU/GPU體系的大規模并行計算、湍流的數值模擬、兩相流的數值模擬,E-mail：wangxian＠mail.xjtu.edu.cn, 李躍明（西安交通大學航天航空學院機械結構強度與振動國家重點實驗室,西安碑林 710049）

基于格子Boltzmann方法的平板射流大渦模擬

上官燕琴, 王 嫻??通訊作者：王嫻（1977－）,女,漢族,吉林,副教授,博士生導師,主要研究方向：基于CPU/GPU體系的大規模并行計算、湍流的數值模擬、兩相流的數值模擬,E-mail：wangxian＠mail.xjtu.edu.cn, 李躍明
（西安交通大學航天航空學院機械結構強度與振動國家重點實驗室,西安碑林 710049）

應用多GPU技術,將格子Boltzmann方法與大渦模擬相結合（LBM-LES）,使用1.12×108網格,對雷諾數Re＝4 000,傾斜角α＝30°,吹風比M＝0.5工況下的平板單孔射流進行了大規模高性能數值模擬研究.合理的定性與定量結果驗證了LBM-LES模擬平板射流的有效性與可行性.使用上億的計算網格捕捉了精細的湍流擬序結構,有利于主流與射流之間的摻混機理研究.此外,使用6個K20M GPU并行計算,模擬了71 680 LBM時間步長,僅耗時15 402秒,計算性能達到521.24MLUPS,即每秒更新5.212 4×108個網格點的數據.

格子Boltzmann方法（LBM）；大渦模擬（LES）；多GPU；三維平板射流（JICF）

0 引言

平板射流（JICF）是工程中極其常見的一種流動模型,也是一種高度復雜的湍流流動.雖然它的邊界條件很簡單,但可作為眾多具有重要工程價值的復雜流動的簡化模型,例如：煙縷擴散,燃燒器的燃油噴射以及渦輪燃氣發動機葉片的氣膜冷卻等[1].因此研究平板射流的流動并分析主流與射流之間的摻混機理,在實際工程應用中有著十分重要的意義.迄今為止,已有很多關于平板射流的實驗研究與數值研究結果,但由于其流動情況的復雜性以及實驗與數值方法的局限性,至今對于平板射流的摻混機理的認知程度還極其有限.為了得到更精確的平板射流數值結果,我們嘗試應用格子Boltzmann方法（LBM）結合多GPU并行技術,使用上億網格計算平板射流的流動.

LBM是近三十年來興起的一種新的計算流體力學數值模擬方法.由于它具有天然的并行性、適于處理復雜幾何邊界問題和容易編程等優點[2],已逐漸形成一項引人注目的數值模擬技術.目前,LBM已在計算流體力學領域取得了很大的成功,它在常規流場的模擬,多孔介質,多相流,多組分流以及電磁流等領域有著很廣泛的應用前景[3],但其在模擬湍流的可行性與正確性方面還有待進一步驗證[4－5].

平板射流是高度復雜的湍流流動,湍流模型的選取對其計算精度的影響很大.目前關于平板射流的數值研究大部分都是基于雷諾平均方程（RANS）的,但Acharya等人已證明RANS結合湍流模型無法很好地捕捉平板射流中的流動動態從而使其過小地預測了氣膜冷卻中溫度場的側向傳熱[6].隨后Tyagi等人首次嘗試使用大渦模擬（LES）計算平板射流,并得到了精確的計算結果[7].之后越來越多的研究人員開始使用LES模擬平板射流及氣膜冷卻并得到了精確的結果[8－10].本文將LES與LBM相結合用于平板射流的數值研究.

大量的計算資源耗費是當前湍流精確計算的瓶頸之一,可編程GPU（GPGPU）的出現在一定程度上解決了這個瓶頸問題.現在,GPGPU已發展成為一種高度并行化、多線程、多核、多種編程接口的處理器[11],這樣僅需在普通的個人計算機上安裝一個或多個GPU接口單元便可完成大規模的并行計算,大大方便了計算流體力學工作者進行相關工作.這也使得結合多GPU并行技術完成高精度湍流數值模擬成為一種可能和趨勢.

本文的主要工作是：基于LBM-LES并結合多GPU并行技術使用上億網格計算三維平板射流的流動,得到高度復雜湍流流場的精細擬序結構,分析其摻混機理.

1 數值計算方法

1.1 格子Boltzmann方程

在LBM中,首先,將Boltzmann方程在分離的網格點中離散成速度分布[12]

其中,fi是離散點的速度分布函數,ei是離散點第i方向的速度,與此同時,N是每個離散點所含速度方向的總數.常見模型有：D2Q9,D3Q13,D3Q15,D3Q19,D3Q27.上式中的Ωi是碰撞項,用Boltzmann-BGK近似[13]可以得到

將其代入式（1）,得到

其中,feqi是平衡分布函數,λ是松弛時間.

為了得到準確的計算結果,合適的平衡分布函數的選取是很重要的,本文選用的平衡分布函數

本文使用D3Q19模型,其速度分布為：e0（0,0,0）,e1（1,0,0）,e2（－1,0,0）,e3（0,1,0）,e4（0,－1,0）, e5（0,0,1）,e6（0,0,－1）,e7（1,1,0）,e8（－1,1,0）,e9（1,－1,0）,e10（－1,－1,0）,e11（1,0,1）,e12（－1,0, 1）,e13（1,0,－1）,e14（－1,0,－1）,e15（0,1,1）,e16（0,－1,1）,e17（0,1,－1）,e18（0,－1,－1）.

圖1 D3Q19模型Fig.1 D3Q19 LBM model

聲速為cs＝1/3,與此同時,對于理想氣體而言,其壓強為

進一步在空間x方向和時間t上離散方程（3）,可以得到它的完全離散形式

式（8）就是著名的LBGK模型.其中,τ＝λ/δt是無量綱的松弛時間.粘性υ可以從上式中推導得到

一般情況下,假設δt＝1,式（8）可以演化成以下兩個基本步驟碰撞步：

遷移步：

其中,fi和分別表示碰撞前與碰撞后的分數.

1.2 大渦模擬：Smagorinsky亞格子應力模型

大渦模擬（LES）的基本思想是通過對控制方程濾波,以實現波的大小尺度分離,再直接求解大尺度的渦旋,對于小空間尺度的渦旋則利用亞格子應力模型進行求解[14].Smagorinsky亞格子應力模型是常用的一種模型,該模型假設雷諾應力僅依賴于局部的應變率.

經過濾波之后的格子Boltzmann方程[4,15－16]變成

在標準Smagorinsky模型中,渦粘性υt可以通過濾波長度尺度Δx和過濾后的應變率張量＝（＋）/2計算得到

式（16）中的CS是Smagorinsky系數.本文參照作者之前基于LBM-LES的壁面約束流動數值模擬[17],取CS＝0.13.此外,本文選擇了有限差分方法計算應變率張量.

2 計算模型

圖2為本文計算模型.其中,x、y、z分別表示流向、展向與垂向,對應方向的計算區域長度為Lx、Ly與Lz.計算網格數Lx×Ly×Lz＝896×320×390,總網格數約為1.12×108.射流孔孔徑D為平板展向長度的1/8,即D＝Ly/8.射流孔中心即為坐標原點,且射流孔中心位于主流出口下游5D的Ly/2平面上.u∞為自由來流速度,uj為射流速度.吹風比被定義為M＝ρjuj/ρ∞u∞（假設主流密度與射流密度相同）,其中,射流出口是均勻出口速度分布,uj是射流在冷卻孔出口界面內的平均速度.射流速度方向相對于主流速度方向的傾斜角度為α.雷諾數定義為Re＝ρu∞D/υ,υ為運動粘性系數.本次計算工況設定為：Re＝4 000,α＝30°,M＝0.5.

圖2 計算模型Fig.2 Flow configuration

邊界條件的設置為：主流入口和射流入口均采用入口速度邊界條件,平板為無滑移壁面,主流出口采用對流出口邊界條件,展向前后邊界為周期性邊界條件,上邊界采用充分發展邊界條件.

3 結果與討論

3.1 多GPU的計算性能

計算網格設置為896×320×390,總計算網格量達到1.12×108.計算的工況：雷諾數Re＝4 000,射流傾斜角α＝30°,吹風比M＝0.5.本文采用了6個K20M GPU并行計算,應用CUDA-MPI進行數據傳輸[18－19],模擬了7.168×104LBM時間步長,耗時15 402秒,計算性能達到521.24 MLUPS.

3.2 湍流的擬序結構

經過一個多世紀的研究表明：湍流是多尺度有結構的不規則流體運動.也就是說,湍流并不是完全的不規則運動,其多尺度不規則的運動中存在著可辨識的、有明確統計周期與外形的流動結構——擬序結構.這種擬序運動主宰著湍流的脈動生成與發展,因此擬序結構對流場的混合以及流體的擴散等起著重要作用[20－21].所以,平板射流流場中擬序結構的精細捕捉與分析有助于研究主流－射流摻混機理.

圖3是選用Q判據作為漩渦識別方法[22]得到的1.792×104LBM時間步長的瞬態湍流擬序結構.在該圖中可以清晰地看到平板射流流場中的典型擬序結構：位于射流前緣的馬蹄渦,位于射流孔之后的發卡渦以及剪切層渦.

圖3 1.792×104LBM步長的瞬態湍流擬序結構Fig.3 Instantaneous coherent structures at1.792×104LBM-steps

從擬序結構圖上可以看到,射流從斜圓孔射出后立即就形成一些初級結構,然后開始向下游移動并形成一系列的發卡渦.這與Tyagi等在文獻[7]中對平板射流的近壁處邊界層的研究討論相一致.而且在射流上風面附近區域有幾段較短的渦管.射流與主流橫向速度之間的相互作用是形成渦管的原因.這些渦管從射流孔周圍產生后沿著主流流向彎曲,并附著于發卡渦上方.由于此算例中射流與主流的橫向速度之間的相互作用比較弱,所以渦管產生之后,斷斷續續地附著在發卡渦上方.圖3給出的合理的湍流擬序結構圖說明了LBM-LES模擬平板射流的有效性與可行性.

圖4顯示的是在相同計算工況下由不同網格量計算得到的同一時刻（1.792×104LBM時間步長）的湍流擬序結構.圖4（a）顯示的是使用本文的計算網格設置計算得到的結果,總網格量為1.12×108.圖4（b）則是由512×200×256的計算網格設置計算得到的結果,總網格量為2.62×107,前者使用的網格量約是后者的4.27倍.由對比結果可以明顯看到：當使用的計算網格規模達到108,可以精細地捕捉到二級渦結構甚至是三級渦結構；而使用的計算網格規模為107時,只能夠得到清晰的一級渦結構.

3.3 二次流動特性分析

圖5顯示在y＝Ly/2平面的1.792×104LBM步長瞬態相對壓力云圖,其中定義相對壓力為：P＝（p－p∞）/ρu∞.從該圖上我們可以清楚地看到射流射出之后就在射流孔后緣出現了一個較強的低壓區,之后也出現了幾個低壓區,隨著位置的后移,其強度變弱.這些低壓區其實是反對稱渦生成的區域,其強度也對應代表了反對稱渦的強度.反對稱渦是在垂直于流動方向的平面生成的,屬于二次流動,其對射流與主流之間的摻混起主要作用.由該結果可以看到,反對稱渦在射流孔中心位置開始產生,然后沿著流向移動,但強度逐漸減弱.

圖4 不同計算網格規模下同一瞬態時刻（1.792×104LBM時間步長）湍流擬序結構對比圖Fig.4 Instant（1.792×104LBM steps）coherent structures computed with different grids

圖5 y＝Ly/2平面的1.792×104LBM步長瞬態相對壓力云圖Fig.5 Instantaneous pressure contours above the bottom plate in plane of y＝Ly/2 at 1.792×104LBM-steps

在圖6中可以看到分別位于斜圓孔下游0倍直徑即圓孔中心位置（圖6（a））、1倍直徑（圖6（b））、3倍直徑（圖6（c））、5倍直徑（圖6（d））、7倍直徑（圖6（e））以及9倍直徑（圖6（f））處的垂直于平板的縱剖面的1.792×104LBM步長瞬態流線圖和渦量云圖.由該結果可以看到,在射流孔中心位置上已經開始形成反對稱渦,這與圖5顯示的在射流孔中心位置后即刻出現的強低壓區相一致.而且,與其它不同位置的渦量相比,該處的渦量值是最大的,也就說明這是的反對稱渦的強度是最大的.反對稱渦是由主流與射流之間的剪切效應而產生的,使主流被射流夾帶到射流下方.而由圖6可以看到,距離射流孔越遠,反對稱渦渦強度越小.這也與由圖5得到的結論相一致.而且隨著與射流孔的距離越來越大,反對稱渦的尺寸也越來越大,然而流向位置為x/D＝5.0是一個轉折點,在x/D＝5.0之后反對稱渦的強度與尺寸都隨著與射流孔的距離的變大而變小.這是因為在x/D/5.0以及之后的位置已經位于耗散區內,湍流耗散在該區域占主導作用.也就意味著在耗散區主流與射流之間由于速度差而產生的剪切效應逐漸因湍流耗散而消失.

3.4 平均速度與定量比較

為了更好地說明本文使用的計算程序的準確性,我們參照了?zcan和Larsen等人完成的平板射流實驗[23－24],定量地比較了計算結果.該實驗測量的工況是：雷諾數Re＝2 400,傾斜角α＝90°,吹風比為M＝3.31.我們計算得到y＝Ly/2平面上不同流向位置的平均速度結果并與試驗結果比較.

圖6 不同流向位置的垂直于平板的縱剖面的1.792×104LBM步長瞬態流線圖和渦量云圖Fig.6 Contours of instantaneous vorticity component in streamwise directionωxand streamlines onvarious cross sections along streamwise direction at1.792×104LBM-steps

圖7展示了y＝Ly/2平面上不同流向位置（a）和（c）x/D＝－0.5和（b）和（d）x/D＝0.5的平均流向速度以及平均垂向速度,其中實線代表的是LBM-LES的計算結果,而實心三角形代表的是實驗結果[23－24].為了便于定量結果的比較,圖7中顯示的均是無量綱結果,其中平均速度使用主流速度u∞進行無量綱化,而垂向距離使用孔徑D進行無量綱化.由圖7可以清楚地看到：除了鄰近壁面區域（z/D＜2.0）計算結果與實驗結果之間有較大偏差外,其余區域的計算結果與實驗結果均可較好地吻合.該圖從定量上說明了使用LBMLES可以準確地計算平板射流這種復雜的流動情況.

在射流孔出口前緣位置（x/D＝－0.5）,流向速度u與垂向速度w除了z/D＜0.5的區域以外,其余區域均還保持著原本的主流速度分布情況,這是由于在射流孔前緣上主流與射流還未開始摻混,所以射流對主流速度分布的影響區域極小.在出口后緣即x/D＝0.5位置上,流向速度u在近壁面區域為負值,這是因為當射流從出口噴射出來后,在后緣位置產生了“回流”現象.

4 結論

應用多GPU技術,將格子Boltzmann方法與大渦模擬相結合（LBM-LES）,采用D3Q19單松弛時間模型,使用了1.12×108計算網格,得到并分析了雷諾數Re＝4 000,傾斜角α＝30°,吹風比為M＝0.5情況下三維平板單孔射流的精細的湍流擬序結構和不同流向位置的二次流動結果,定性地說明了LBM-LES模擬平板射流這種復雜流動情況的合理性.與此同時,參照?zcan和Larsen等人完成的橫向射流實驗[23－24],其實驗工況為：雷諾數Re＝2 400,傾斜角α＝90°,吹風比為M＝3.31.比較了y＝Ly/2平面上不同流向位置的平均速度,合理的計算結果定量地說明了LBM-LES可以準確地計算平板射流.

多GPU技術的應用,使計算性能得到了顯著的提高,使得上億網格的湍流非穩態計算在數小時之內即可完成,有效地解決了高精度湍流數值模擬研究發展的瓶頸問題之一：計算消耗大.使用上億的計算網格量不僅能夠捕捉到精確的初級擬序渦結構,還可以得到精細的次級擬序渦結構,有助于主流與射流之間摻混機理的研究.

圖7 y＝Ly/2平面上不同流向位置Fig.7 Profiles ofmean streamwise（a）and（b）and wall-normal（c）and（d）velocity components in plane of y＝Ly/2 at locations near jet exit x/D＝－0.5 and x/D＝0.5

[1] Bogard D G,Thole K A.Gas turbine film cooling[J].Journal of Propulsion and Power,2006,22（2）：249－270.

[2] 何雅玲,王勇,李慶.格子Boltzmann方法的理論及應用[M].北京：科學出版社,2009.

[3] Chen SY,Doolen G D.Lattice Boltzmannmethod for fluid flows[J].Annual Review of Fluid Mechanics,1998,30：329－364.

[4] Yu H,Girimaji S S,Luo L S.DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method[J].Journal of Computational Physics,2005,209：599－616.

[5] Yu H,Luo L S,Girimaji S S.LES of turbulent square jet flow using an MRT lattice Boltzmann model[J].Computers＆Fluids,2006,35：957－965.

[6] Acharya S,Tyagi M,Hoda A.Flow and heat transfer predictions for film cooling[J].Annals of the New York Academy of Sciences,2001,934（1）：110－125.

[7] TyagiM,Acharya S.Large eddy simulation of film cooling flow from an inclined cylindrical jet[J].Journal of Turbomachinery,2003,125（4）：734－742.

[8] Guo X,Schr?derW,Meinke M.Large-eddy simulations of film cooling flows[J].Computers＆Fluids,2006,35（6）：587 －606.

[9] Renze P,Schr?der W,Meinke M.Large-eddy simulation of film cooling flows with variable density jets[J].der Flow, Turbulence and Combustion,2008,80（1）：119－132.

[10] Sakai E,Takahashi T,Watanabe H.Large-eddy simulation of an inclined round jet issuing into a crossflow[J].International Journal of Heat and Mass Transfer,2014,69：300－311.

[11] Nvidia Corporation.NVIDIA CUDA compute unified device architecture programming guide（Version 0.8.2）[EB/OL].[2015－01－18].http：//moss.csc.ncsu.edu/～mueller/cluster/nvidia/0.8/NVIDIA_CUDA_Programming_Guide_0.8.2, pdf,2007－04－24.

[12] Cercignani C.Theory and application of the boltzmann equation[M].Scottish Academic Press,London,UK,1975.

[13] Bhatnagar P L,Gross E P,Krook M.A model for collision processes in gases.I.Small amplitude processes in charged and neutral one-component systems[J].Physical Review,1954,94（3）：511－525.

[14] Galperin B,Orszag S A.Large eddy simulation of complex engineering and geophysical flows[M].New York：Cambridge University Press,1993.

[15] Hou S,Sterling J,Chen S,Doolen G D.A lattice Boltzmann subgrid model for high Reynolds number flows[J].Pattern Formation and Lattice Gas Automata,1996,6：151－166.

[16] Somers JA.Direct simulation of fluid flow with cellular automata and the lattice-Boltzmann equation[J].Applied Scientiic Research,1993,51（1/2）：127－133.

[17] Wang X,Shangguan Y,Ondera N,Kobayashi H,Aoki T.Direct numerical simulation and large eddy simulation on a turbulentwall-bounded flow using lattice Boltzmann method and multiple GPUs[J].Mathematical Problems in Engineering, 2014：742432.

[18] Ogawa S,Aoki T.GPU Computing for 2-dimensional incompressible-flow simulation based on multi-grid method[J].Transactions of JSCES,2009：20090021.

[19] Wang X,Aoki T,Multi-GPU performance of incompressible flow computation by lattice Boltzmann method on GPU cluster [J].Parallel Computing,2011,37：521－535.

[20] Pope SB.Turbulent flows[M].Cambridge：Cambridge University Press,2000.

[21] 邱翔,劉宇翔.湍流的相干結構[J].自然雜志,2004,26（4）：187－193.

[22] Dubief Y,Delcayre F.On coherent-vortex identification in turbulence[J].Journal of Turbulence,2000,1（1）：1－22.

[23] ?zcan O,Larsen P S.An experimental study of a turbulent jet in cross-flow by using LDA[M].Lyngby：DTU Mechanical Engineering,2001.

[24] ?zcan O,Larsen P S.Laser Doppler anemometry study of a turbulent jet in crossflow[J].AIAA Journal,2003,41（8）：1614－1616.

High-Performance Numerical Simulation of Jet in Cross-Flow Based on Lattice Boltzmann M ethod

SHANGGUAN Yanqin,WANG Xian,LIYueming
（State Key Laboratory for Strength and Vibration ofMechanical Structures,School of Aerospace, Xi′an Jiaotong University,28 Xianning West Road,Xi′an 710049,China）

Large eddy simulation（LES）with 1.12×108mesheswas performed on three-dimensional jet in cross-flow（JICF）using lattice Boltzmann method（LBM）and multiple Graphic Processing Units（Multi-GPUs）.Reynolds number based on free-stream velocity and jet diameter is 4 000,streamwise inclination angle of jet isα＝30°,and jet-to-cross-flow velocity ratio is set as 0.5.Validity and feasibility of LBM-LESon simulating jet in cross-flow were verified by reasonable qualitative and quantitative results.Fine coherent structures were captured by large-scaled simulation which benefits study on mixing mechanism of jet-in-cross-flow（JICF）.Moreover,6 K20M GPUs were adopted in simulation and it took 15 402 seconds for LBM-GPU solver to simulate 71 680 LBM steps resulting in calculated performance of 521.24 MLUPS.

lattice Boltzmannmethod（LBM）；large eddy simulation（LES）；multiple Graphic Processing Units（Multi-GPUs）；jet in cross-flow（JICF）

1001-246X（2015）06-0669-08

O358

A

2015－03－11；

2015－03－14

國家973計劃（2003CB30570202）和國家自然科學基金（11302165）資助項目

上官燕琴（1991－）,女,博士研究生,主要研究方向：湍流的數值模擬,E-mail：sgyq.6510336＠stu.xjtu.edu.cn