上官燕琴, 王 嫻通訊作者：王嫻（1977－）,女,漢族,吉林,副教授,博士生導(dǎo)師,主要研究方向：基于CPU/GPU體系的大規(guī)模并行計(jì)算、湍流的數(shù)值模擬、兩相流的數(shù)值模擬,E-mail：wangxian＠m(xù)ail.xjtu.edu.cn, 李躍明（西安交通大學(xué)航天航空學(xué)院機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安碑林 710049）

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

上官燕琴, 王 嫻??通訊作者：王嫻（1977－）,女,漢族,吉林,副教授,博士生導(dǎo)師,主要研究方向：基于CPU/GPU體系的大規(guī)模并行計(jì)算、湍流的數(shù)值模擬、兩相流的數(shù)值模擬,E-mail：wangxian＠m(xù)ail.xjtu.edu.cn, 李躍明
（西安交通大學(xué)航天航空學(xué)院機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安碑林 710049）

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

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

0 引言

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

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

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

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

本文的主要工作是：基于LBM-LES并結(jié)合多GPU并行技術(shù)使用上億網(wǎng)格計(jì)算三維平板射流的流動(dòng),得到高度復(fù)雜湍流流場(chǎng)的精細(xì)擬序結(jié)構(gòu),分析其摻混機(jī)理.

1 數(shù)值計(jì)算方法

1.1 格子Boltzmann方程

在LBM中,首先,將Boltzmann方程在分離的網(wǎng)格點(diǎn)中離散成速度分布[12]

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

將其代入式（1）,得到

其中,feqi是平衡分布函數(shù),λ是松弛時(shí)間.

為了得到準(zhǔn)確的計(jì)算結(jié)果,合適的平衡分布函數(shù)的選取是很重要的,本文選用的平衡分布函數(shù)

本文使用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,與此同時(shí),對(duì)于理想氣體而言,其壓強(qiáng)為

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

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

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

遷移步：

其中,fi和分別表示碰撞前與碰撞后的分?jǐn)?shù).

1.2 大渦模擬：Smagorinsky亞格子應(yīng)力模型

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

經(jīng)過(guò)濾波之后的格子Boltzmann方程[4,15－16]變成

在標(biāo)準(zhǔn)Smagorinsky模型中,渦粘性υt可以通過(guò)濾波長(zhǎng)度尺度Δx和過(guò)濾后的應(yīng)變率張量＝（＋）/2計(jì)算得到

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

2 計(jì)算模型

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

圖2 計(jì)算模型Fig.2 Flow configuration

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

3 結(jié)果與討論

3.1 多GPU的計(jì)算性能

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

3.2 湍流的擬序結(jié)構(gòu)

經(jīng)過(guò)一個(gè)多世紀(jì)的研究表明：湍流是多尺度有結(jié)構(gòu)的不規(guī)則流體運(yùn)動(dòng).也就是說(shuō),湍流并不是完全的不規(guī)則運(yùn)動(dòng),其多尺度不規(guī)則的運(yùn)動(dòng)中存在著可辨識(shí)的、有明確統(tǒng)計(jì)周期與外形的流動(dòng)結(jié)構(gòu)——擬序結(jié)構(gòu).這種擬序運(yùn)動(dòng)主宰著湍流的脈動(dòng)生成與發(fā)展,因此擬序結(jié)構(gòu)對(duì)流場(chǎng)的混合以及流體的擴(kuò)散等起著重要作用[20－21].所以,平板射流流場(chǎng)中擬序結(jié)構(gòu)的精細(xì)捕捉與分析有助于研究主流－射流摻混機(jī)理.

圖3是選用Q判據(jù)作為漩渦識(shí)別方法[22]得到的1.792×104LBM時(shí)間步長(zhǎng)的瞬態(tài)湍流擬序結(jié)構(gòu).在該圖中可以清晰地看到平板射流流場(chǎng)中的典型擬序結(jié)構(gòu)：位于射流前緣的馬蹄渦,位于射流孔之后的發(fā)卡渦以及剪切層渦.

圖3 1.792×104LBM步長(zhǎng)的瞬態(tài)湍流擬序結(jié)構(gòu)Fig.3 Instantaneous coherent structures at1.792×104LBM-steps

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

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

3.3 二次流動(dòng)特性分析

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

圖4 不同計(jì)算網(wǎng)格規(guī)模下同一瞬態(tài)時(shí)刻（1.792×104LBM時(shí)間步長(zhǎng)）湍流擬序結(jié)構(gòu)對(duì)比圖Fig.4 Instant（1.792×104LBM steps）coherent structures computed with different grids

圖5 y＝Ly/2平面的1.792×104LBM步長(zhǎng)瞬態(tài)相對(duì)壓力云圖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步長(zhǎng)瞬態(tài)流線圖和渦量云圖.由該結(jié)果可以看到,在射流孔中心位置上已經(jīng)開(kāi)始形成反對(duì)稱渦,這與圖5顯示的在射流孔中心位置后即刻出現(xiàn)的強(qiáng)低壓區(qū)相一致.而且,與其它不同位置的渦量相比,該處的渦量值是最大的,也就說(shuō)明這是的反對(duì)稱渦的強(qiáng)度是最大的.反對(duì)稱渦是由主流與射流之間的剪切效應(yīng)而產(chǎn)生的,使主流被射流夾帶到射流下方.而由圖6可以看到,距離射流孔越遠(yuǎn),反對(duì)稱渦渦強(qiáng)度越小.這也與由圖5得到的結(jié)論相一致.而且隨著與射流孔的距離越來(lái)越大,反對(duì)稱渦的尺寸也越來(lái)越大,然而流向位置為x/D＝5.0是一個(gè)轉(zhuǎn)折點(diǎn),在x/D＝5.0之后反對(duì)稱渦的強(qiáng)度與尺寸都隨著與射流孔的距離的變大而變小.這是因?yàn)樵趚/D/5.0以及之后的位置已經(jīng)位于耗散區(qū)內(nèi),湍流耗散在該區(qū)域占主導(dǎo)作用.也就意味著在耗散區(qū)主流與射流之間由于速度差而產(chǎn)生的剪切效應(yīng)逐漸因湍流耗散而消失.

3.4 平均速度與定量比較

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

圖6 不同流向位置的垂直于平板的縱剖面的1.792×104LBM步長(zhǎng)瞬態(tài)流線圖和渦量云圖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的平均流向速度以及平均垂向速度,其中實(shí)線代表的是LBM-LES的計(jì)算結(jié)果,而實(shí)心三角形代表的是實(shí)驗(yàn)結(jié)果[23－24].為了便于定量結(jié)果的比較,圖7中顯示的均是無(wú)量綱結(jié)果,其中平均速度使用主流速度u∞進(jìn)行無(wú)量綱化,而垂向距離使用孔徑D進(jìn)行無(wú)量綱化.由圖7可以清楚地看到：除了鄰近壁面區(qū)域（z/D＜2.0）計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果之間有較大偏差外,其余區(qū)域的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果均可較好地吻合.該圖從定量上說(shuō)明了使用LBMLES可以準(zhǔn)確地計(jì)算平板射流這種復(fù)雜的流動(dòng)情況.

在射流孔出口前緣位置（x/D＝－0.5）,流向速度u與垂向速度w除了z/D＜0.5的區(qū)域以外,其余區(qū)域均還保持著原本的主流速度分布情況,這是由于在射流孔前緣上主流與射流還未開(kāi)始摻混,所以射流對(duì)主流速度分布的影響區(qū)域極小.在出口后緣即x/D＝0.5位置上,流向速度u在近壁面區(qū)域?yàn)樨?fù)值,這是因?yàn)楫?dāng)射流從出口噴射出來(lái)后,在后緣位置產(chǎn)生了“回流”現(xiàn)象.

4 結(jié)論

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

多GPU技術(shù)的應(yīng)用,使計(jì)算性能得到了顯著的提高,使得上億網(wǎng)格的湍流非穩(wěn)態(tài)計(jì)算在數(shù)小時(shí)之內(nèi)即可完成,有效地解決了高精度湍流數(shù)值模擬研究發(fā)展的瓶頸問(wèn)題之一：計(jì)算消耗大.使用上億的計(jì)算網(wǎng)格量不僅能夠捕捉到精確的初級(jí)擬序渦結(jié)構(gòu),還可以得到精細(xì)的次級(jí)擬序渦結(jié)構(gòu),有助于主流與射流之間摻混機(jī)理的研究.

圖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方法的理論及應(yīng)用[M].北京：科學(xué)出版社,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] 邱翔,劉宇翔.湍流的相干結(jié)構(gòu)[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

國(guó)家973計(jì)劃（2003CB30570202）和國(guó)家自然科學(xué)基金（11302165）資助項(xiàng)目

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