格子Boltzmann方程模擬高雷諾數三維方腔流

康海貴，李承功

（大連理工大學海岸和近海工程國家重點實驗室，大連116023）

格子Boltzmann方程模擬高雷諾數三維方腔流

康海貴，李承功

（大連理工大學海岸和近海工程國家重點實驗室，大連116023）

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

多松弛時間格子Boltzmann方程；標準Smagorinsky渦粘性模型；高雷諾數三維方腔流；GPU；CUDA

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

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

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

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

1 控制方程

1.1 MRT-LBE

對于D維空間具有Q個粒子運動方向的多松弛時間格子Boltzmann模型的演化方程為［19］

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

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

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

將標準的Smagorinsky渦粘性模型與MRT-LBE組合的關鍵在于建立粒子分布函數與應變率張量之間的關系，即通過非平衡分布函數的二階矩計算應變率張量［25］

1.3 邊界條件

本文采用動量改進的半步反彈格式處理無滑移的墻邊界［19］

2 數值模擬和結果討論

2.1 數值驗證

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

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

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

2.2 Re=10 000三維方腔流數值結果及結果討論

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

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

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

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

2.3 計算效率

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

3 結論

本文利用D3Q19 MRT-SMAG模型對三維方腔內的流場進行了數值模擬研究，并采用最新的基于GPU的CUDA并行技術以提高MRT-SMAG模型在模擬三維問題時的計算效率，將模型計算結果與實驗和其他數值結果進行比較驗證。可以得到如下結論：

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

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

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

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

［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.），Fields 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

國家自然科學基金（50679008）

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

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

TV 131；O 242.1

A

1005-8443（2013）05-0453-08