劉祖斌,趙鵬
(1浙江大學海洋工程系,杭州310058;2英偉達半導體科技(上海)有限公司,上海201210)
結合大渦模擬的格子玻爾茲曼方法模擬高雷諾數(shù)流動
劉祖斌1,趙鵬2
(1浙江大學海洋工程系,杭州310058;2英偉達半導體科技(上海)有限公司,上海201210)
為了將格子玻爾茲曼方法(Lattice Boltzmann method)推廣到高雷諾數(shù)流動的數(shù)值模擬上,文章研究了將大渦模擬的Smagorinsky模型及其改進的IR模型應用到LBM方法中并發(fā)展了LBM-SM模型和LBM-IR模型,其后以頂蓋驅動方腔流動(Re=1 000;10 000;100 000)和后臺階流動(Re=389;1 000;10 000;100 000)為標準進行數(shù)值模擬與比較。數(shù)值結果說明該方法可以應用于高雷諾數(shù)流動的數(shù)值模擬,并克服了傳統(tǒng)LBM方法在模擬高雷諾數(shù)流場時會出現(xiàn)非物理振蕩的情況。通過分析對比,不僅得到了兩種模型能夠穩(wěn)定計算且不失真的參數(shù)范圍,并且認為LBM-IR模型有效參數(shù)的選取范圍更為廣泛。
Lattice Boltzmann method;大渦模擬;高雷諾數(shù);LBM-SM;LBM-IR
格子玻爾茲曼方法(LBM,Lattice Boltzmann Method)是一種新興的計算流體力學數(shù)值模擬方法。它從分子運動論和統(tǒng)計力學的觀點和理論出發(fā),建立微觀粒子的運動模型,用流動和碰撞兩個簡單的過程模擬流動的宏觀特性[1-3]。該方法具有計算過程簡單,復雜邊界問題處理十分靈活方便,計算局部性強,具有天然的并行性,非常適合在大規(guī)模并行計算機上運行[4]。因此,自從上世紀末期該方法發(fā)展以來,對其的研究與應用也成為了當今計算流體力學熱點。然而問題同樣也存在,LBM方法是一種直接數(shù)值模擬方法(DNS,Direct Numerical Simulation),在足夠細的計算網(wǎng)格下可以模擬高雷諾數(shù)的流動問題,但是受到計算條件以及計算成本的限制,計算往往很難做到非常大的網(wǎng)格規(guī)模,同時由于在導出Navier-Stocks方程的過程中有低馬赫數(shù)(Ma<0.3)的要求。綜合以上兩點,在使用LBM方法進行數(shù)值模擬時,流場雷諾數(shù)會受到一定限制。
大渦模擬是當今流體力學研究的熱門課題之一,作為研究湍流的主要手段之一在未來將有更為廣泛的應用[5]。大渦模擬的基本思想是采用各種濾波方法將那些控制動量輸運等的大尺度渦與小尺度渦分離。大尺度渦應用數(shù)值模擬方法直接求解,可以獲得大尺度流動的動態(tài)特征,而小尺度渦則采用亞網(wǎng)格模式模擬。使用大渦模型所需要的網(wǎng)格數(shù)目很大,若要求能分辨出最小尺度也就是Kolmogorov尺度的渦,則網(wǎng)格數(shù)目將正比于Re9/4,因此在高雷諾數(shù)流動情況下的計算效率非常低下。
將大渦模擬的模型應用于LBM方法中是解決此問題的一個有效方法。這個思想,可追溯到LBM方法建立初期,陳十一等[6]曾進行過高雷諾數(shù)方腔流動的模擬,后來國內外有不少從事大渦模擬研究的學者[7-11]也結合起來研究,但所應用的模型太少,沒有引起足夠的重視。本文將大渦模擬的Smagorinsky模型及其改進的IR模型應用到LBM方法中,分別對二維的頂蓋驅動方腔流動和后臺階流動進行了數(shù)值模擬。方腔流動的雷諾數(shù)的取值范圍從1 000到100 000,并且與已有的數(shù)值標準解進行了對比,其結果符合得較好;作為該方法的擴展,本文給出了高雷諾數(shù)的后臺階流動,雷諾數(shù)分別取了10 000,100 000,并給出了相應的分析。
結合大渦模擬的LBM方法(LBM-LES)可以在很大程度上提高數(shù)值模擬中的雷諾數(shù),又具有直接數(shù)值模擬的精細分辨率和LBM的良好并行性能,將會推動直接數(shù)值模擬(DNS)的更廣泛應用。
LBM方法是一種模擬流體流動的數(shù)值方法,它不再基于連續(xù)介質假設,而是把流體看成由許多只有質量沒有體積的微粒組成,這些微粒可以向空間的若干方向任意流動。它的主要思想就是以簡單規(guī)則的微觀粒子運動代替復雜多變的宏觀現(xiàn)象。粒子在每個時間步的運動由流動和碰撞兩步構成:碰撞,當多個粒子到達同一網(wǎng)格結點時,它們按碰撞規(guī)則相互作用并重新分布;流動,每個粒子按其速度方向移動到最近的網(wǎng)格結點。
1992年,Qian(錢躍竑)等人[2-3]提出LBM方法的單松弛模型,即:LBGK(Lattice Bhatnagar-Gross-Krook)模型。目前,LBGK模型是格子玻爾茲曼方法中應用最廣泛的模型。LBGK模型來源于LGA(Lattice Gas Automator),繼承了空間離散速度離散的原則和流動碰撞的簡單運動,拋棄了Fermi-Dirac分布及粒子碰撞理論,采用Maxwell平衡態(tài)分布和BGK碰撞理論(由Bhatnager,Gross,Krook提出),使得各向同性、Galilean不變性和壓力獨立于流速等條件都得到滿足。再用單松弛時間逼近碰撞項,其對應的玻爾茲曼演化方程為:
其中:ω為松弛因子,ω在0到2之間,小于1為亞松弛,大于1為超松弛。fi為單粒子分布函數(shù),表示在t時間上,x格點上速度為ei的有1個粒子的概率;ei稱為粒子運動速度或離散速度;為平衡態(tài)分布函數(shù),是流體動力學量(如密度,動量等)的函數(shù)。方程左側為流動項,右側即為碰撞項。由此方程可以看出,碰撞過程只與本節(jié)點有關,為局部計算過程,流動步也僅與相鄰的網(wǎng)格點之間進行數(shù)據(jù)交換。因此,該方法具有天然的并行性能。LBGK模型中的DnQb系列模型的分布函數(shù)為:
其中:cs為聲速,ti為對應粒子速度的權系數(shù),u為宏觀速度,ρ為宏觀密度。為了保證得到正確的宏觀Navier-Stocks方程,需要通過使?jié)M足質量守恒、動量守恒等約束條件來得到正確的權系數(shù)。
DnQb系列模型的粘性系數(shù)為:
對于本文所選用的D2Q9模型,權系數(shù)和離散速度的選取如下:
2.1 大渦模擬及其模型
大渦模擬的基本思想是直接數(shù)值模擬大尺度結構,而小尺度結構(其尺度小于計算網(wǎng)格尺度)對流場的大尺度結構不敏感依賴,通過建立亞網(wǎng)格模型(Subgrid Scale,SGS)來模擬。
其中:τij為湍流亞網(wǎng)格應力
其次,大渦模擬最重要的就是給出正確的亞網(wǎng)格模型。最早提出且應用至今的是Smagorinsky模型,這是由氣象學家Smagorinsky于1963年提出的。該模型認為亞格子湍流具有混合長度型的渦粘系數(shù),混合長度與過濾尺度同一量級,屬于唯象渦粘模型。其優(yōu)點是簡單,并可以通過調整參數(shù)獲得相當好的結果;缺點是完全屬于耗散型的,沒有考慮從小尺度到大尺度的能量逆轉。后人在此基礎發(fā)展了諸多模型,如Smagorinsky-Lilly模式,及動力學模式等[7]。
Smagorinsky模型中亞網(wǎng)格應力中的非各向同性部分正比于湍流渦粘性系數(shù)與大尺度變形率張量的乘積,
其中:δij是Kronedcker δ函數(shù),vt是渦粘系數(shù),C是Smagorinsky常數(shù),是大尺度應變率張量的模:
由于Smagorinsky常數(shù)的取值對計算的結果影響較大,而通常的模型如Smagorinsky-Lilly模式對取值缺乏判定的依據(jù),所以針對該常數(shù)的機理性探討非常有價值。2005年Mayers和Saguat[11]提出一個改進的Inetial-Range模型,認為模型系數(shù)C強烈依賴于兩個比值,即:大渦模擬的濾波尺度與Kolmogorov尺度的比值Δ/η和宏觀尺度與濾波尺度的比值L/Δ,具體的關系參照文獻。在實際應用中,他們給出了在無窮大雷諾數(shù)條件下的極限值C∞=0.18,在以下計算中,我們將對此進行驗證。
2.2 格子Boltzmann大渦模擬方法(LBM-LES)
在LBM方法中引入大渦模擬模型的方法,首先給出LBM方法中的濾波形式,對粒子的分布函數(shù)fi(平衡態(tài)函數(shù))定義濾波后的形式為:
然后對Boltzmann方程的微分形式進行濾波,得到濾波Boltzmann方程
在LBGK模型中,所有格子的間距均為1,這樣就使得物理空間的濾波對格點上的分布函數(shù)保持不變,且松弛時間與平均自由程等價,所以松弛時間的改變就等價于改變局部平均自由程,從而改變了局部粘性??偟恼承韵禂?shù)vt與松弛時間τt的關系仍然保持通常方法中的形式:
物理運動粘性在流動過程中保持不變,而流動結構的變化會改變渦粘性系數(shù)??偟恼承韵禂?shù)由物理運動粘性系數(shù)v0和渦粘性系數(shù)vv組成:
Smagorinsky模型[7-10],
Inetial-Range模型[11-14],
其中:C和C*為模型常數(shù),其取值范圍在0-2之間,Δ取網(wǎng)格長度為1。
從以上計算公式中可知,局部渦粘性系數(shù)是由局部應變率張量來求得,而相比于其他差分法的復雜,在格子玻爾茲曼方法中可通過局部非平衡態(tài)分布函數(shù)簡單求得,這也是相比于單純大渦模擬的一個優(yōu)勢。結果如下:
濾波非平衡動量流張量,
濾波局部應變率張量,
因此得到加入兩個模型的格子波爾茲曼方法的最終松弛時間的表達式:
LBM-Smagorinsky模型(LBM-SM),
LBM-Inetial Range模型(LBM-IR),
3.1 標準頂蓋方腔驅動流
標準頂蓋方腔驅動流是經典的粘性不可壓縮流動數(shù)值模擬的算例,從上個世紀六十年代初被提出后已經研究了四十余年,該模型問題通過假定流體四周均為固體壁面而避免了物理上復雜的進口和出口。由于它整個流場形狀規(guī)則,邊界條件簡單,因此便于計算、分析和比較[15],成為檢驗計算流體力學方法的最廣泛算例。
該問題的數(shù)學模型是假設一個寬度為1的正方形腔內充滿流體,左右和底面均固定,在頂部有一平板以固定速度由左向右移動。本文使用的計算分辨率為256×256,平板移動速度為0.1。在本文以下的結果中,將速度與網(wǎng)格均做歸一化處理,以便于與標準解[16]對比。
驅動方腔流動過程中存在轉捩過程,在高雷諾數(shù)下由層流轉變?yōu)橥牧鞯倪^程中存在一個雷諾數(shù)區(qū)間[Re1,Re 2],5 000<Re1<=7 500,Re2>10 000。當Re在[Re1,Re 2]時,驅動方腔的流動存在非定常周期解;當Re<Re1時存在穩(wěn)定的層流解,而當Re>Re2時驅動方腔的流動完全轉變?yōu)橥牧鳌?/p>
下面將給出雷諾數(shù)為1 000,10 000,100 000的結果。作為對基本算法的檢驗,對于雷諾數(shù)為1 000和10 000的數(shù)值模擬,我們比較了DNS和兩種不同的組合模型的結果,并與基準解(Benchmark)進行了比照,驗證了模型的準確性。而對于大雷諾數(shù),由于DNS已無法模擬,直接給出數(shù)值模擬的結果。本文在計算中專門對模型系數(shù)C進行了敏感性測試,以0.02為步長對0-2的區(qū)域均進行了模擬計算。
Re=1 000時,流動較穩(wěn)定,屬于層流階段。本文采用了三種方法進行了數(shù)值模擬,如圖1,其中方格、三角、虛線、實線分別代表數(shù)值基準解、LBM方法直接數(shù)值模擬、LBM-SM模型和LBM-IR模型結果,并且給出了兩條中心線上的速度分布值。由圖1可以看出標準的LBM方法與LBM-SM模型和LBM-IR模型所得的結果都對標準解吻合得很好。在模型系數(shù)選取方面,LBM-SM模型在全部取值都能完成計算,其中C=0.1能得到最好的結果;LBM-IR模型測試結果是同樣的情況。
當Re=10 000時,圖2給出了其幾何中心線上的速度圖,并與標準解進行對比,符合依然很好。但在細節(jié)上速度產生了周期性現(xiàn)象,圖4給出了方腔左下方點(0.78,0.78)速度周期變化的圖示。LBMSM模型在模型系數(shù)小于0.1區(qū)域能得到正確的流場,在0.04附近能得到最好的效果,而大于0.1會逐漸抹平邊角區(qū)的小渦結構。LBM-IR模型在大于0.1依然能得到較好的效果,在C=0.18附近效果最好。計算結果說明,模型系數(shù)C的合理取值范圍較小是標準Smagorinsky模型的一個缺點,而且系數(shù)的選取沒有合理的依據(jù),而IR模型在取值的合理性和有效性方面有較大提高。
圖1 Re=1 000時,幾何中心垂直線上的速度u分布和速度v分布Fig.1 Re=1 000,u-velocity and v-velocity distribution in geometry center vertical line
圖2 Re=10 000時,幾何中心垂直線上的速度Ux分布和速度Uy分布Fig.2 Re=10 000,Ux-velocity and Uy-velocity distribution in geometry center vertical line
圖3 Re=10 000時,左下方點(0.78,0.78)X和Y方向速度隨時間的演化圖Fig.3 Re=10 000,evolution of Uxand Uywith time at bottom left dot(0.78,0.78)
Re=100 000時,流動已完全進入了湍流階段。由于此時LBM方法會由于松弛因子太接近于2而出現(xiàn)很強的非物理數(shù)值現(xiàn)象,LBM-SM模型能保證不出現(xiàn)數(shù)值發(fā)散的區(qū)域也相對較小,在C=0.1附近;LBM-IR模型相對區(qū)域較大,在C=0.28附近。
3.2 后臺階流動
分離再附現(xiàn)象是實際工程中經常遇到的流動問題,如機翼分離渦的再附、沖壓發(fā)動機燃燒室內的流動、流體經過截面突擴的管道后的流動等等。為改善流動狀態(tài)、加強燃燒、強化換熱以及有利于混合物的摻混等,對這一流動進行研究是十分必要和有意義的。后臺階流動由于其分離點穩(wěn)定、邊界條件簡單,實驗方便,成為研究分離再附運動中最簡單也是最典型的流動。
圖4 Re=10 000時,左下方點(0.78,0.78)X方向速度隨時間的演化圖(局部放大圖)Fig.4 Re=10 000,enlarged partial view of evolution of Ux-velocity with time at bottom left dot(0.78,0.78)
在后臺階流動中,臺階后緣會產生類似平板剪切層的大渦序列,上下壁面處會出現(xiàn)一系列附著渦。從而,臺階后緣的剪切層,擴張段的分離區(qū)以及下游充分發(fā)展的槽道流,成為后臺階流動的三個主要流動組成。從20世紀70年代開始,許多研究者都對后臺階流動進行了實驗研究和數(shù)值模擬[17-19]。
后臺階流動的粘性不可壓流體的控制方程為:
這里u=(u,v,)w為速度場,p為壓力,雷諾數(shù)Re定義為:
其中:Um為入口最大速度,D為特征尺度,此處為臺階高度,即D=h,h為臺階高度,ν為流體的運動粘性系數(shù)。
圖5 后臺階流動的不同位置的速度剖面圖,(a)Re=389,(b)Re=1 000Fig5 Velocity profiles of backward facing step flow in variety of places,(a)Re=389,(b)Re=1 000
本文應用LBM-IR模型對二維高雷諾數(shù)的流動做了數(shù)值模擬。使用傳統(tǒng)的DNS進行數(shù)值模擬時,模擬高雷諾數(shù)流動需要加密計算網(wǎng)格,使得計算時間、計算成本都很高,但將LES與LBM結合的模型能夠較好地解決這一問題。本文在4 000×100的網(wǎng)格上模擬了Re為389,1 000,10 000,100 000的流動。如圖5所示,在雷諾數(shù)為389和1000的結果中,我們比較了不同位置速度剖面與實驗值[15]的比照,驗證了該模型的準確性。
如圖6所示,當雷諾數(shù)達到10 000,100 000的時候,在臺階后方形成了一系列旋渦,與理論的預測一致。圖(a),(b)分別給出了兩個雷諾數(shù)時臺階之后的流動瞬時渦量云圖。此時的流動已經發(fā)展為湍流,傳統(tǒng)的LBM方法不能直接進行數(shù)值模擬。而通過本文提出的LBM-LES方法能夠在不增加分辨率的情況下捕捉到流動的細節(jié)信息。由圖可以看出,在臺階之后先出現(xiàn)了一些小渦,隨著時間的推進,小渦逐漸合并形成了大結構的渦。并且我們的算法也能夠獲取大渦中更為復雜的結構。
圖6 后臺階流動的渦量分布圖,(a)Re=10 000,(b)Re=100 000Fig.6 Voticity contours of backward facing step flow,(a)Re=10 000,(b)Re=100 000
本文將LES湍流模型加入到LBM方法中組合成為LBM-LES模型不僅在理論上是完全可行的,以上算例也很好地驗證了這類模型的準確性。通過對高雷諾數(shù)的標準頂蓋驅動方腔流動和后臺階流動的模擬,并與標準解進行對比以及理論分析可以得到以下結論:
(1)使用LBM-SM模型和LBM-IR模型進行模擬過程中,分別得到了不同雷諾數(shù)的最佳模型系數(shù)取值范圍,而這個結果并不完全是一個確定的值,是與雷諾數(shù)相關的。
(2)在兩個組合模型的對比中可以發(fā)現(xiàn)LBM-IR模型不僅具有更完整的系數(shù)確定依據(jù),而且在系數(shù)取值范圍上有較大的提高。
本文研究的兩個模型具有普遍意義,對于大渦模擬的渦粘模型可以依據(jù)本文理論推導加入到LBM方法中,因此LSM-LES模型可以得到可持續(xù)的發(fā)展。而對于該類模型的詳細研究還是很有必要,尤其在模型系數(shù)敏感性上,尋求最佳的系數(shù)將對模擬結果的精確度產生關鍵的作用,本文作者將深入研究。
[1]郭照立,鄭楚光,李青.流體動力學的格子Boltzmann方法[M].武漢:湖北科學技術出版社,2002.
[2]Qian Y H,d’Humieres D,Lallemand P.Lattice BGK models for the Navier-Stokes equation[J].Euro phys Lett.,1992(17): 479-484.
[3]Qian Y,Succi S,Orszag S.Recent advances in lattice Boltzmann computing[J].Ann.Rev.Comp.Phys.,1995,3:195-242.
[4]趙鵬,張丹丹,汪魯兵等.格子Boltzmann并行程序的優(yōu)化與性能分析[J].微電子學與計算機,2008,25(10):185-188.
[5]崔桂香,許春曉,張兆順.湍流大渦數(shù)值模擬進展[J].空氣動力學報,2004,6,22(2):121-129. Cui Guixiang,Xu Chunxiao,Zhang Zhaoshun.Progress in large eddy simulation of turbulent flows[J].ACTA Aerodynamica Sinica,2004,6,22(2):121-129.
[6]Hou S,Sterling J,Chen S,Doolen G D.A lattice Boltzmann subgrid model for high Reynolds number flows[M].In: Lawniczak AT,Kapra lR,editors.Pattern formation and lattice gas automata.Fields Institute Communications 6,Providence,RI:AMS,1996:66-151.
[7]Guan H,Wu C J.Large-eddy simulations of turbulent cavity flows with the dynamic SGS model and LBM algorithm[J]. Transactions of Nanjing Univ.of Aeronautics&Astronautics,2001,18(Suppl.):60-63.
[8]楊帆,劉樹紅,唐學林,吳玉林.格子Boltzmann亞格子模型的研究[J].工程熱物理學報,2004,S1. Yang Fan,Liu Shuhong,Tang Xuelin,Wu Yulin.Study of subgrid turbulence model for lattice Boltzmann method[J]. Journal of Engineering Thermophysics,2004,S1.
[9]Yu H,Girimaji S S,Luo L S.DNS and LES of decaying isotropic turbulence with and without frame rotation using Lattice Boltzmann methods[J].J Comput.Phys.,2005,209:599-616.
[10]Burattini P,Lavoie P,Agrawal A,Djenidi L,Antonia R A.Power law of decaying homogeneous isotropic turbulence at low Reynolds number[J].Physical Review E,2006,73,066301.
[11]Meyers J,Sagaut P.On the model coefficients for the standard and the variational multi-scale Smagorinsky model[J].J Fluid Mech,2006,569:287-319.
[12]Meryers J,Sagaut P.Evaluation of Smagorinsky variants in large-eddy simulations of wall-resolved plane channel flows [J].Physics of Fluids,2007,19(9):1-12.
[13]Dong Y H,Sagaut P,Simon M.Inertial consistent subgrid model for large-eddy simulation based on lattice Boltzmann method[J]Physics of Fluids,2008,20(3):1-12.
[14]董宇紅,周亦航,鄧義求.基于格子Boltzmann方程的大渦模擬對湍流時空關聯(lián)性的研究[J].北京理工大學學報, 2012,32(8):876-880. Dong Yuhong,Zhou Yihang,Deng Yiqiu.Time-space corrections of isotropic turbulence by lattice Boltzmann based large eddy simulation[J].Transactions of Beijing Institute of Technology,2012,32(8):876-880.
[15]Ghia U,Ghia K N,Shin C T.High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]J Comput.Phys.,1982,48:387-411.
[16]Armaly B F,Durst F,Rereira J C F,Schonung B.Experimental and theoretical investigation of backward facing step[J].J Fluid Mech.,1983,127:473-496.
[17]Kim J,Moin P.Application of a fractional-step method to incompressible Navier-Stokes equations[J].J Comput.Phys., 1985,59:308-323.
[18]王小華,樊洪明,何鐘怡.后臺階流動的數(shù)值模擬[J].計算力學學報,2003,20(3):361-365. Wang Xiaohua,Fan Hongming,He Zhongyi.Numerical simulation of backward facing step flow[J].Chinese Journal of Computational Mechanics,2003,20(3):361-365
[19]王廣超,廖國勇.后臺階流動的非均勻格子Boltzmann方法模擬[J].計算機與現(xiàn)代化,2007,1:6-9. Wang Guangchao,Liao Guoyong.Non-uniform lattice Boltzmann method simulation of backward-facing step flow[J].Computer&Modernization,2007,1:6-9.
High-Re flow simulation with Lattice Boltzmann method combined with LES
LIU Zu-bin1,ZHAO Peng2
(1 Departement of Ocean Engineering,Zhejiang University,Hangzhou 310058,China; 2 NVIDIA Semiconductor Technology(Shanghai)Co.,Ltd,Shanghai 201210,China)
A kind of lattice Boltzmann method(LBM)combined with LES models was proposed.Smagorinsky(SM)model and improved model Inetial Range(IR)model were used to extend LBM computation scales for high Reynolds number.Simulations are performed for typical cases of high Reynolds number,cavity flow (Re=1 000;10 000;100 000)and back facing step flow(Re=389;1 000;10 000;100 000),and results are compared with benchmark results.The comparison shows that both the LBM-SM and the LBM-IR models behaved well to match the accurate solutions and could conquer the problem of oscillation caused by LBM. The analysis and comparison of the two models illustrate that the LBM-IR model has an advantage of wider area of valid model coefficient.
Lattice Boltzmann method;Large Eddy Simulation;High Reynolds number;LBM-SM;LBM-IR
O35
A
10.3969/j.issn.1007-7294.2015.05.002
1007-7294(2015)05-0484-09
2014-12-05
國家自然科學基金(10625210);上海市教委項目基金(08ZZ43)
劉祖斌(1983-),男,博士后,E-mail:liuzb@zju.edu.cn;
趙鵬(1980-),男,高級GPU架構工程,E-mail:patricz@nvidia.com。