康海貴，李承功

（大連理工大學(xué)海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，大連116023）

格子Boltzmann方程模擬高雷諾數(shù)三維方腔流

康海貴，李承功

（大連理工大學(xué)海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，大連116023）

為分析方腔流內(nèi)部流場的特性和驗(yàn)證格子Boltzmann方法模擬湍流的能力，應(yīng)用標(biāo)準(zhǔn)Smagorinsky渦粘性模型與多松弛時(shí)間格子Boltzmann方程（Multiple Relaxation Time Lattice Boltzmann Equation，MRT-LBE）組合對高雷諾數(shù)（Re=10 000）三維方腔流進(jìn)行數(shù)值研究，計(jì)算了時(shí)間平均量，如速度，均方根速度、雷諾應(yīng)力以及中心斷面（y=W/2）處的流線等高線。模擬結(jié)果與已有實(shí)驗(yàn)和數(shù)值模型結(jié)果比較可知，MRT-LBE能夠精準(zhǔn)地計(jì)算剪切驅(qū)動(dòng)方腔內(nèi)流場的變化。另外，將基于圖形處理（graphic processor unit，GPU）的計(jì)算統(tǒng)一設(shè)備架構(gòu)（Compute Unified Device Architecture，CUDA）并行技術(shù)引入到基于MRT-LBE的Smagorinsky模型以提高計(jì)算效率，計(jì)算效率提高達(dá)200倍。

多松弛時(shí)間格子Boltzmann方程；標(biāo)準(zhǔn)Smagorinsky渦粘性模型；高雷諾數(shù)三維方腔流；GPU；CUDA

格子Boltzmann方法［1］（Lattice Boltzmann Method，LBM）是克服了格子自動(dòng)機(jī)方法（Lattice Gas Cellular Automata，LGCA）一系列缺點(diǎn)并繼承其優(yōu)點(diǎn)發(fā)展而來的一種新的數(shù)值方法。與傳統(tǒng)的求解Navier-Stokes動(dòng)量方程組的數(shù)值方法相比，LBM從介觀的角度通過分布函數(shù)演化流體粒子的碰撞和遷移過程，從而描述流體的宏觀運(yùn)動(dòng)，該方法具有邊界處理簡單，易于編程和適合并行的優(yōu)點(diǎn)，因此引起了各領(lǐng)域?qū)W者的極大關(guān)注，并已經(jīng)被成功應(yīng)用于模擬眾多復(fù)雜流體流動(dòng)現(xiàn)象［2-3］。

近年來，LBM還被用來解決海洋和海岸工程領(lǐng)域的相關(guān)問題。例如Zhou等（2000）提出了淺水方程的LBM模型及包含湍流模型的形式，并利用該模型模擬了圓柱繞流，渠流和潮流等水動(dòng)力問題［4］。張金鳳和張慶河（2011）利用LBM對河口海岸細(xì)粒黏性泥沙進(jìn)行數(shù)值研究，分析了三維分形絮團(tuán)的沉降特性，對球形顆粒在靜水中沉降引起的紊動(dòng)流場進(jìn)行了數(shù)值模擬［5］。在開發(fā)利用海洋能發(fā)電上，王樹杰等（2011）應(yīng)用LBM對水輪機(jī)的引導(dǎo)槽裝置進(jìn)行了水動(dòng)力優(yōu)化［6］。Janβen等（2012）應(yīng)用基于完全非線性淺水波方程的LBM數(shù)值模擬了海洋長波的傳播和爬坡過程［7］。最近，LBM與VOF方法相結(jié)合直接求解Navier-Stokes方程已經(jīng)被用來計(jì)算自由表面流體流動(dòng)問題，如破碎波［8］等。

方腔流通常是指頂部平板以恒定速度驅(qū)動(dòng)規(guī)則區(qū)域內(nèi)封閉的不可壓流體（例如水）的流動(dòng)，由于在方腔流的流動(dòng)中可以觀察到幾乎所有可能發(fā)生在不可壓流體中的流動(dòng)現(xiàn)象（如渦，二次流，Taylor-G?rtler-like（TGL）渦，流體分岔，不穩(wěn)定，瞬變以及湍流等），因此對方腔流已經(jīng)進(jìn)行了廣泛的實(shí)驗(yàn)和數(shù)值研究［9-21］。最初的方腔流研究始于Burggraf等（1966），他們應(yīng)用解析和數(shù)值的方法計(jì)算了封閉區(qū)域內(nèi)由具有恒定速度的頂部平板驅(qū)動(dòng)的穩(wěn)定流場［9］，隨后Pan等（1967）在實(shí)驗(yàn)中應(yīng)用攝影技術(shù)觀測了矩形方腔內(nèi)的流場結(jié)構(gòu)，但是他們忽略了三維作用對流場的影響［10］。Koseff等（1984）最早應(yīng)用可視化技術(shù)對雷諾數(shù)Re≤10 000的方腔流進(jìn)行了深入的實(shí)驗(yàn)研究，詳細(xì)討論了腔內(nèi)流場的三維特征，角渦，TGL渦，流場的不穩(wěn)定狀態(tài)，邊墻摩擦對腔內(nèi)流場的影響等［11-13］。Prasad和Koseff（1989）根據(jù)Koseff和Street的實(shí)驗(yàn)結(jié)果分析了表征湍流脈動(dòng)強(qiáng)度的二階統(tǒng)計(jì)量（均方根速度，雷諾應(yīng)力），進(jìn)一步考察了邊墻，雷諾數(shù)和TGL渦對腔內(nèi)流場的影響［14］。2000年Shankar和Deshpande對方腔流給出了詳細(xì)的綜述概括［15］。最近，Guermond等（2002）對三維方腔流（寬長高比為1：1：2，Re=1 000）的初始階段流場進(jìn)行了實(shí)驗(yàn)和數(shù)值研究［16］。Leriche等（2000）［17］和Bouffanais等（2007）［18］對Re=12 000的三維方腔流分別進(jìn)行了直接數(shù)值模擬和大渦模擬研究。在應(yīng)用LBM數(shù)值計(jì)算上，d′Humières等（2002）應(yīng)用MRT-LBE對Re≤4 000三維對角驅(qū)動(dòng)方腔流進(jìn)行模擬證明了MRT模型的的穩(wěn)定性要由于LBGK模型［19］。De等（2009）應(yīng)用MRT-LBE對Re≤1 000的三維方腔流進(jìn)行數(shù)值模擬，以研究在長寬比A=1條件下不同深寬比對腔內(nèi)三維流場結(jié)構(gòu)的影響和評(píng)估MRT模型計(jì)算三維流場的能力［20］。Premnath等（2009）應(yīng)用廣義格子Boltzmann方程模擬了Re=12 000的三維方腔流［21］。

本文首先對三維方腔流初始流場進(jìn)行求解以驗(yàn)證標(biāo)準(zhǔn)的Smagorinsky渦粘性模型［22］與MRT-LBE［19］組合的MRT-SMAG模型的有效性，其次為了評(píng)估LBM模擬湍流的能力為今后建立數(shù)值波浪水槽提供技術(shù)基礎(chǔ)，利用該模型對高雷諾數(shù)三維方腔流進(jìn)行數(shù)值實(shí)驗(yàn)，并與物理模型實(shí)驗(yàn)結(jié)果和其他數(shù)值結(jié)果進(jìn)行比較驗(yàn)證，并考察了寬高比對腔內(nèi)流場的影響。為了提高LBM求解三維流場的計(jì)算效率和針對LBM易于并行的特性，采用了最新的GPU-CUDA并行計(jì)算［23］以縮減模型中粒子的碰撞和遷移時(shí)間。

1 控制方程

1.1 MRT-LBE

對于D維空間具有Q個(gè)粒子運(yùn)動(dòng)方向的多松弛時(shí)間格子Boltzmann模型的演化方程為［19］

本文采用D3Q19格子模型［19］（圖1），因此離散的粒子速度矢量為

1.2 基于Smagorinsky渦粘性模型的MRT-LBE模型

大渦模擬（Large Eddy Simulation，LES）［24］的主要原理是直接求解大尺度（或稱可解尺度）的湍流運(yùn)動(dòng)，而小尺度（或稱不可解尺度）的紊動(dòng)作用則利用相應(yīng)的渦粘性模型模擬出來以減少計(jì)算耗費(fèi)。本文主要采用的Smagorinsky（1963）［22］所提出的渦粘性模型求解殘余應(yīng)力，該模型的主要形式為

將標(biāo)準(zhǔn)的Smagorinsky渦粘性模型與MRT-LBE組合的關(guān)鍵在于建立粒子分布函數(shù)與應(yīng)變率張量之間的關(guān)系，即通過非平衡分布函數(shù)的二階矩計(jì)算應(yīng)變率張量［25］

1.3 邊界條件

本文采用動(dòng)量改進(jìn)的半步反彈格式處理無滑移的墻邊界［19］

2 數(shù)值模擬和結(jié)果討論

2.1 數(shù)值驗(yàn)證

三維不可壓方腔流流動(dòng)的計(jì)算區(qū)域如圖2所示，在長寬高分別為L，W和H的規(guī)則區(qū)域內(nèi)充滿了不可壓縮流體（水），區(qū)域頂部平板以恒定的速度Up從沿x軸方向左向右移動(dòng)，從而驅(qū)動(dòng)區(qū)域內(nèi)流體的流動(dòng)。區(qū)域頂部為速度邊界，其余為靜止的墻邊界，因此采用動(dòng)量改進(jìn)的半步反彈格式進(jìn)行處理，具體形式為（14）式。應(yīng)用LBM進(jìn)行模擬計(jì)算通常需要將實(shí)際的物理量轉(zhuǎn)化為無量綱的格子量［27］，下面定義r，t，u和υ分別表示長度，時(shí)間，速度和運(yùn)動(dòng)粘度系數(shù)，相應(yīng)的大寫字母L，T，U表示對應(yīng)的特征量，角標(biāo)為p的量是有量綱的物理量，而角標(biāo)為lb的量表示無量綱的格子量。

本節(jié)首先對三維方腔流初始階段（即從0時(shí)刻開始到第12周期時(shí)間內(nèi)，并采用進(jìn)行歸一化≤12）流場進(jìn)行模擬以驗(yàn)證D3Q19 MRT-SMAG模型的有效性，其中方腔的長寬高分別為=6.2 cm，=12.4 cm，=6.2 cm，頂部驅(qū)動(dòng)速度為1.8 cm/s，雷諾數(shù)Re=1 000，根據(jù)雷諾數(shù)的定義可知水的運(yùn)動(dòng)粘度系數(shù)=1.116×10-6m2/s。對計(jì)算區(qū)域采用127×255×127的規(guī)則格子網(wǎng)格進(jìn)行離散，模型中頂部驅(qū)動(dòng)速度=0.1以確保（馬赫數(shù)Ma=＜0.3），因此根據(jù)單位轉(zhuǎn)換公式可以推導(dǎo)出其他格子量。Re=1 000的三維方腔流流動(dòng)表現(xiàn)為層流狀態(tài)，因此Smagorinsky常數(shù)取零，碰撞對角矩陣中的松弛率分別取為

由圖3可知，方腔對稱截面（y=W/2）上分別沿x=L/2和z=H/2的不同時(shí)刻（t=4，6，8，10，12）的瞬時(shí)水平和垂直速度的計(jì)算結(jié)果與實(shí)驗(yàn)和其他數(shù)值方法的結(jié)果符合較好，這表明了MRT-SMAG模型可以對流場的初始階段進(jìn)行求解。但是圖3所示的數(shù)值結(jié)果（本文模型的和文獻(xiàn)［16］中基于求解Navier-Stokes方程的）與實(shí)驗(yàn)結(jié)果都略有不同，這可能是由實(shí)驗(yàn)中頂部平板速度的微小變化所導(dǎo)致的。

2.2 Re=10 000三維方腔流數(shù)值結(jié)果及結(jié)果討論

方腔內(nèi)的流動(dòng)形式主要取決于方腔幾何尺寸（如寬高比）和雷諾數(shù)（在方腔尺寸和封閉流體粘度系數(shù)不變的情況下反映頂部驅(qū)動(dòng)速度的大小）。通常認(rèn)為雷諾數(shù)為2 000時(shí)，腔內(nèi)流場為穩(wěn)定的層流，當(dāng)雷諾數(shù)在2 000～3 000時(shí)，方腔底部的角渦表現(xiàn)出不穩(wěn)定狀態(tài)，隨著雷諾數(shù)的增大，腔內(nèi)流體的不穩(wěn)定狀態(tài)更加劇烈，當(dāng)雷諾數(shù)在6 000～8 000時(shí)，腔內(nèi)流場首次表現(xiàn)出湍流特征（如擴(kuò)散現(xiàn)象），當(dāng)雷諾數(shù)大于等于10 000時(shí)，方腔底部的角渦完全發(fā)展為湍流［13-14］。因此，為了評(píng)估MRT-SMAG模型模擬湍流的能力為今后建立基于LBM的數(shù)值波浪水槽提供技術(shù)基礎(chǔ)，本文利用D3Q19 MRT-SMAG模型對Re=10 000的三維方腔流進(jìn)行數(shù)值計(jì)算，對比了方腔寬高比K=W/H=1和3兩種情況下的時(shí)間平均速度以分析邊墻對內(nèi)部流場的影響，給出了K=1方腔流的時(shí)間平均均方根速度以及雷諾應(yīng)力曲線，并與實(shí)驗(yàn)和其他數(shù)值結(jié)果進(jìn)行比較。

腔的長和高分別為L=15 cm，H=15 cm，根據(jù)寬高比可知方腔寬分別取為W=15 cm和45 cm，腔內(nèi)封閉的流體為不可壓縮的水，其運(yùn)動(dòng)粘度系數(shù)取，由雷諾數(shù)定義可以推導(dǎo)出方腔頂部驅(qū)動(dòng)速度≈0.067 m/s。對計(jì)算區(qū)域分別采用127×127×127（K=W/H=1）和100×300×100（K=W/H=3）的規(guī)則格子網(wǎng)格進(jìn)行離散，并定義格子系統(tǒng)中頂部驅(qū)動(dòng)速度0.1，再將實(shí)際物理量轉(zhuǎn)換為無量綱的格子量。Smagorinsky常數(shù)=0.12，碰撞對角矩陣中的松弛率分別取為。為了獲得精確的時(shí)間平均結(jié)果，首先選取作為初始計(jì)算階段以確保腔內(nèi)流場完全發(fā)展，再對150的瞬時(shí)結(jié)果進(jìn)行時(shí)間平均以獲得穩(wěn)定的時(shí)均結(jié)果。選取不同的初始計(jì)算時(shí)間（最大為）和時(shí)間平均周期（最多平均1 000的瞬時(shí)結(jié)果）以分析其對時(shí)間平均量的影響，分析結(jié)果表明選取更長的初始計(jì)算時(shí)間和平均時(shí)間對結(jié)果的影響可以忽略不計(jì)。

圖4給出了利用D3Q19 MRT-SMAG模型模擬Re=10 000三維方腔流得到的時(shí)間平均水平速度＜u＞和垂直速度＜w＞，MRT-SMAG模型結(jié)果與實(shí)驗(yàn)［13-14］和直接數(shù)值模擬結(jié)果［17］符合較好。另外，無論是實(shí)驗(yàn)和直接數(shù)值模擬結(jié)果還是MRT-SMAG模型計(jì)算結(jié)果都表明了中心對稱斷面處寬高比K=3的時(shí)間平均速度要略大于K=1的結(jié)果，這是由邊墻對腔內(nèi)流體的摩擦作用造成，即寬高比越大，邊墻摩擦對中心對稱斷面的影響逐越低，對湍流脈動(dòng)的摻混能力越弱。

本文還計(jì)算了腔內(nèi)流場的二階統(tǒng)計(jì)量以分析湍流脈動(dòng)強(qiáng)度分布，這些量均以和行歸一化。圖5-a表示在中心對稱斷面（y=W/2）上與頂部驅(qū)動(dòng)方向一致的時(shí)間平均均方根速度，即在x= L/2處沿垂向方向（z軸）的分布，而圖5-b表示在中心對稱斷面（y= W/2）上與頂部驅(qū)動(dòng)方向垂直的時(shí)間平均均方根速度，即在z=H/2處沿水平方向（x軸）的分布。另外，圖6-a和6-b分別給出了在中心對稱斷面（y=W/2）上時(shí)間平均雷諾應(yīng)力在x=L/2處沿垂直方向（z軸）的分布和在z=H/2處沿水平方向（x軸）的分布。這些結(jié)果表明了腔內(nèi)湍流是沿墻邊界產(chǎn)生，下游邊界（x=L）附近的湍流脈動(dòng)要比上游邊界（x=0）附近的強(qiáng)烈，而底部強(qiáng)邊界（z=0）附近的湍流脈動(dòng)強(qiáng)度是最大的。這與Leriche等（2000）對三維方腔流的討論是一致的。因此腔內(nèi)湍流強(qiáng)度的分布可由腔內(nèi)動(dòng)量傳遞過程解釋，即速度邊界首先將動(dòng)量傳遞到頂部邊界附近的粘性層當(dāng)中使得頂部邊界的下游處壓力增大，增大的壓力抑制了粘性層的運(yùn)動(dòng)并驅(qū)動(dòng)流體沿著下游墻邊界（x=L）向下運(yùn)動(dòng)，類似于沿平直壁面的射流，該射流沿著下游墻邊界的中心線分解為兩個(gè)近似橢圓的自由射流，隨后以一個(gè)非常小的角度沖擊到底部墻邊界導(dǎo)致湍流脈動(dòng)的產(chǎn)生［17］。通過與已有數(shù)據(jù)進(jìn)行對比可知，應(yīng)用MRT-SMAG模型計(jì)算的結(jié)果與實(shí)驗(yàn)和直接數(shù)值模擬以及大渦模擬的結(jié)果在峰值上略有偏差但是總體趨勢符合較好，表明了該模型有能力對湍流場進(jìn)行求解，其中模擬結(jié)果與已有數(shù)據(jù)的偏差可能是由頂部速度邊界的不同造成的，由于在方腔流的實(shí)驗(yàn)中運(yùn)動(dòng)平板與腔內(nèi)流體接觸會(huì)帶走少量流體，造成頂部邊界的速度分布無法確定，但是也可假設(shè)為常數(shù)，但是文獻(xiàn)［17］和［18］在對方腔流的數(shù)值計(jì)算中采用高階的關(guān)于坐標(biāo)的多項(xiàng)式模擬頂部的速度分布，而本文采用的是恒定的速度邊界。

2.3 計(jì)算效率

求解格子Boltzmann模型主要包括初始化變量，粒子的碰撞，遷移和邊界處理，以及每隔1 000個(gè)迭代步判斷結(jié)果是否收斂三個(gè)步驟。為了改進(jìn)LBM求解三維問題的計(jì)算效率，采用GPU-CUDA并行計(jì)算技術(shù)加速D3Q19 MRT-SMAG模型，具體的實(shí)現(xiàn)過程參見Li等的工作［28］。本文使用的支持CUDA計(jì)算機(jī)的具體配置為：AMD Phenom II 1100T 3.3GHz處理器，16GB內(nèi)存，以及具有1.5GB顯存的GTX580顯卡。軟件環(huán)境為Windows XP操作系統(tǒng)，CUDA Toolkit版本為3.2，GPU驅(qū)動(dòng)263.06?；贕PU的MRT-SMAG模型與僅用CPU進(jìn)行計(jì)算的效率進(jìn)行比較可知相對于僅用CPU處理，應(yīng)用GPU對模型進(jìn)行加速可以得到較高的計(jì)算比（見表1），這對應(yīng)用LBM處理實(shí)際問題具有非常重要的意義。

3 結(jié)論

本文利用D3Q19 MRT-SMAG模型對三維方腔內(nèi)的流場進(jìn)行了數(shù)值模擬研究，并采用最新的基于GPU的CUDA并行技術(shù)以提高M(jìn)RT-SMAG模型在模擬三維問題時(shí)的計(jì)算效率，將模型計(jì)算結(jié)果與實(shí)驗(yàn)和其他數(shù)值結(jié)果進(jìn)行比較驗(yàn)證?？梢缘玫饺缦陆Y(jié)論：

（1）為了驗(yàn)證模型的有效性，采用本文所建立的MRT-SMAG模型對初始階段的三維方腔流進(jìn)行模擬，計(jì)算結(jié)果很好的符合實(shí)驗(yàn)數(shù)據(jù)表明MRT-SMAG可以對三維初始流場進(jìn)行精確的預(yù)測。

（2）為了評(píng)估LBM模擬湍流的能力，采用MRT-SMAG模型求解了Re=10 000的三維方腔流，計(jì)算結(jié)果與已有數(shù)據(jù)符合較好，通過對比寬高比K=1和K=3的時(shí)間平均速度表明寬高比越大，邊墻的摩擦作用對腔內(nèi)湍流脈動(dòng)的摻混能力越弱。另外通過分析表征湍流脈動(dòng)強(qiáng)度的二階統(tǒng)計(jì)量可知腔內(nèi)湍流是沿墻邊界產(chǎn)生，下游邊界（x=L）附近的湍流脈動(dòng)要比上游邊界（x=0）附近的強(qiáng)烈，而底部強(qiáng)邊界（z=0）附近的湍流脈動(dòng)強(qiáng)度是最大的。

（3）利用GPU-CUDA并行技術(shù)對MRT-SMAG模型進(jìn)行計(jì)算加速可以得到較高的加速比，約為200倍。

綜上所述，MRT-SMAG模型能夠較好的對湍流進(jìn)行求解，而GPU-CUDA技術(shù)又能很好的加速模型的計(jì)算效率，因此將兩者結(jié)合并引入可處理自由表面流動(dòng)的VOF技術(shù)以建立數(shù)值波浪水槽具有非常好的應(yīng)用前景和科研價(jià)值。

［1］Succi S.The lattice Boltzmann equation for fluid dynamics and beyond［M］.New York：Oxford University Press，2001.

［2］Kang S K，Hassan Y A.The effect of lattice models within the lattice Boltzmann method in the simulation of wall-bounded turbulent flows［J］.Journal of Computational Physics，2013，232：100-117.

［3］Thpmmes G，Becker J，Junk M，et al.A lattice Boltzmann method for immiscible multiphase flow simulations using the level set method［J］.Journal of Computational Physics，2009，228：1 139-1 156.

［4］Zhou J G.A lattice Boltzmann model for the shallow water equations［J］.Computer Methods in Applied Mechanics and Engineering，2002，191：3 527-3 539.

［5］Zhang JF，Zhang Q H.Lattice Boltzmann simulation of the flocculation process of cohesive sediment due to differential settling［J］.Continental Shelf Research，2011，31：S94-S105.

［6］Wang S，Xu C，Yuan P，et al.Hydrodynamic optimization of channeling device for hydro turbine based on lattice Boltzmann method［J］.Computers and Mathematics with Applications，2011，61：3 722-3 729.

［7］JANBEN C F，GRILLIST，KRAFCZYK M.Efficient simulations of long wave propagation and runup using a LBM approach on GPGPU hardware［C］//Chung J S.Proceedings of the Twenty-second International Offshore and Polar Engineering Conference. Rhodes：ISOPE，2012：17-22.

［8］JANSSEN C F，GRILLI S T，KRAFCZYK M.Modeling of wave breaking and wave-structure interactions by coupling of fully nonlinear potential flow and lattice-Boltzmann models［C］//Chung JS.Proceedings of the Twentieth International Offshore and Polar Engineering Conference.Beijing：ISOPE，2010：686-693.

［9］Burggraf O R.Analytical and numerical studies of the structure of steady separated flows［J］.Journal of Fluid Mechanics，1966，24：113-151.

［10］Pan F，Acrivos A.Steady flows in rectangular cavities［J］.Journal of Fluid Mechanics，1967，28（4）：643-655.

［11］Koseff JR，Street R L.Visualization studies of a shear driven three-dimensional recirculating flow［J］.ASME Journal of Fluids Engineering，1984，106：21-27.

［12］Koseff JR，Street R L.On end wall effects in a lid driven cavity flow［J］.ASME Journal of Fluids Engineering，1984，106：385-389.

［13］Koseff JR，Street R L.The lid-driven cavity flow：a synthesis of qualitative and quantitative observations［J］.ASME Journal of Fluids Engineering，1984，106：390-398.

［14］Prasad A K，Koseff JR.Reynolds number and end-wall effects on a lid driven cavity flow［J］.Physics of Fluids A，1989，1：208-218.

［15］Shankar PN，Deshpande M D.Fluid mechanics in the driven cavity［J］.Annual Review of Fluid Mechanics，2000，32：93-136.

［16］Guermond JL，Migeon C，Pineau G，etal.Start-up flows in a three-dimensional rectangular driven cavity of aspect ratio 1：1：2 at Re=1 000［J］.Journal of Fluid Mechanics，2002，450：169-199.

［17］Leriche E，Gavrilakis S.Direct numerical simulation of the flow in a lid-driven cubical cavity［J］.Physics of Fluids，2000，12（6）：1 363-1 376.

［18］Bouffanais R，Deville M O.Large eddy simulation of the flow in a lid-driven cubical cavity［J］.Physics of Fluids，2007，19：055108.

［19］D′Humières D，Ginzburg I，Krafczyk M，et al.Multiple-relaxation-time lattice Boltzmann models in three dimensions［J］. Philosophical Transactions of the Royal Society A，2002，360：437-451.

［20］De S，Nagendra K，Lakshmisha K N.Simulation of laminar flow in a three dimensional lid driven cavity by lattice Boltzmann method［J］.International Journal of NumericalMethods for Heat&Fluid Flow，2009，19（6）：790-815.

［21］Premnath K N，Pattison M J，Banerjee S.Generalized lattice Boltzmann equation with forcing term for computation of wallbounded turbulent flows［J］.Physical Review E，2009，79：026703.

［22］Smagorinsky J.General circulation experiments with the primitive equations：I.The basic equations［J］.Monthly Weather Review，1963，91：99-164.

［23］Kirk D B，Hwu WW.Programmingmassively parallel processors：a hands-on approach［M］.San Francisco：Morgan Kaufmann，2010.

［24］Pope S.Turbulent flows［M］.New York：Cambridge University Press，2000.

［25］Hou S L，Sterling J，Chen S Y，et al.A lattice Boltzmann subgrid model for high Reynolds number flows［J］.In：Pattern Formation and Lattice Gas Automata，in：A.T.Lawniczak，R.Kapral（Eds.），F(xiàn)ields Inst Comm，vol.6，AMS，Providence，1996，6：151-166.

［26］Krafczyk M，T?lke J，Luo L S.Large-eddy simulations with a multiple-relaxation-time LBE model［J］.International Journal of Modern Physics B，2003，17：33-39.

［27］L?tt J.Hydrodynamic limitof lattice Boltzmann equation［D］.Switzerland：University of Geneva，2007.

［28］LiCG，Maa JPY，Kang H G.Solving generalized lattice Boltzmannmodel for 3-D cavity flows using CUDA-GPU［J］.Science China Physics，Mechanics&Astronomy，2012，10：1 894-1 904.

Simulation of three-dimensional cavity flow with high Reynolds num ber using lattice Boltzmann equation

KANG Hai-gui,LICheng-gong
(State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116023,China)

In order to analyze the features of the flow field in the cavity and access the capability of simulating turbulent flow of the lattice Boltzmann method(LBM),the multiple relaxation time lattice Boltzmann equation(MRT-LBE)with the standard Smagorinsky eddy viscosity model was used for the simulation of three dimensional cavity flow with high Reynolds number(Re=10 000).The time averaged variables,e.g.,the velocity, root mean square velocity,Reynolds stress,and stream line contours on the symmetry plane(y=W/2)were computed.By the comparisons of the present results with the experimental and numerical results,it shows that the changing of the turbulent flow field in the cavity driven by the shear plane could be performed accurately by the MRT-LBE with the Smagorinsky model.In additional,the compute unified device architecture(CUDA)parallel technique on graphic processor unit(GPU)was introduced into MRT-LBE with the Smagorinsky model for improving the computational efficiency,up to 200 times.

multiple relaxation time lattice Boltzmann equation;standard Smagorinsky eddy viscosity model;high Reynolds number three dimensional cavity flow;GPU;CUDA

2012-04-07；

2013-05-12

國家自然科學(xué)基金（50679008）

康海貴（1945-），男，遼寧省大連人，教授，主要從事海岸工程方面的研究。

Biography：KANG Hai-gui（1945-），male，professor.

TV 131；O 242.1

A

1005-8443（2013）05-0453-08