高雷諾數流動模擬的LBM方法

陳彥曉,李孝偉,丁 玨,翁培奮

(上海大學上海市應用數學和力學研究所,上海200072)

格子Boltzmann方法(lattice Boltzmann method,LBM)是一種從介觀角度出發對流動進行數值模擬的方法,在過去的幾十年里已經發展得相對成熟[1].傳統的LBM方法建立在均勻的網格上,故也稱為標準的LBM方法,該方法要求在每一個時間步長上,微粒從一個格點嚴格運動到其相鄰的另一個格點[2].格子Boltzmann方法使用標準均勻網格,具有原理簡單、實施方便且并行性能好等優點[3],但也存在計算精度不可調、復雜邊界處理困難等缺點.

為了避免標準LBM方法中均勻網格的缺陷,國內外學者提出了許多改進方法[4].這些方法主要分為兩種:第一種是采用分區網格或自適應直角網格[5],在需要高精度求解的局部區域直接建立或通過自適應技術構造細化的直角網格,該辦法能夠有效地解決計算精度不可調的問題,主要優勢在于每個網格的運算法則與標準的LBM完全相同,只是在非均勻網格交界面上通過插值進行必要的數值傳遞;第二種是采用貼體網格.將貼體網格用于LBM方法較好地解決了曲面邊界和網格加密的問題,由此也產生了一種補充插值LBM(interpolation supplemented LBM,ISLBM)方法[6].Imamura等[7]在此基礎上又發展出了一種廣義形式補充插值LBM(generalized interpolation LBM,GILBM),可以適用于更一般的幾何外形.1997年He等[8]首次使用GILBM方法對圓柱繞流進行了模擬計算,結果證明:該方法在Re<104的亞臨界區內具有較好的數值穩定性;而當雷諾數進一步提高時,數值計算容易失穩.為了改善LBM方法處理實際問題時受高雷諾數限制的缺陷,有研究人員提出采用多松弛時間模式的LBM方法[9].研究表明,多松弛時間模式對算法的計算穩定性有一定的改進,但所能模擬的最大雷諾數僅擴大4倍.

基于此,本工作發展了一種在貼體網格GILBM框架下的當地時間步改進方法,該方法是基于傳統當地時間步法[10]發展起來的,通過結合分子運動論[11]的相關原理,使計算域內各點都可以采用滿足該點CFL條件所允許的最大時間步長.本工作對高雷諾數下NACA0012翼型繞流的數值模擬表明,采用當地時間步改進方法可以改善高雷諾數流動模擬的計算穩定性,并有效地提高了計算效率.

1 標準LBM

格子Boltzmann采用了微觀尺度上的Boltzmann方程,在BGK近似中引入了碰撞算子(Maxwell-Boltzmann分布函數),

式中,?t為時間步長.格子Boltzmann方法控制方程可以分為兩部分:碰撞和遷移.碰撞是指不同速度的粒子在到達同一格點時將依據BGK碰撞模型的法則,改變自己的分布函數.遷移是指粒子碰撞后將按固定的速度運動到最近的節點,其粒子分布函數也會隨之運動.因此,格子Boltzmann方法控制方程的演化方程分為兩部分:碰撞項和遷移項.

式中,w0=4/9,w1=w2=w3=w4=1/9,w5=w6=w7=w8=1/36,

宏觀變量可以通過分布函數fi描述,其中密度ρ和速度u分別為

通過式(7)可求得壓強為

式中,運動黏度υ和雷諾數相關,即

其中U0為遠場速度,L為特征速度.

2 GILBM

GILBM是基于ISLBM思想的廣義坐標下的LBM模型.借鑒傳統計算流體力學中的方法,引入計算區域的概念.物理區域上的坐標點(x,y)對應于計算區域上的點(ξ,η).物理平面上是貼體網格,而計算平面上的網格是均勻直角網格.網格間隔?ξ,?η一般均取1.

在物理平面上,離散速度是常矢量,但在計算平面上格點的離散速度隨位置變化而變化,為變矢量.將計算平面上格點的離散速度記為ei,?(?=0,1),根據計算平面與物理平面的關系,可以得出離散速度在物理域和計算域的轉換關系.

求得遷移運動速度在計算平面上的分量后,通過對ei,?在時間步長的積分

得到物理平面上碰撞后粒子在計算平面相應的位置.由于該位置一般并不位于網格格點上,因此可以通過距離該點最近的網格節點上的分布函數插值得到該點上的分布函數數值.

雖然可以選擇不同的插值格式進行計算,但為了避免虛假的線性速度,線性插值的數值精度是遠遠不夠的.為了在保證計算精度的前提下力求計算穩定,在流場內部區域節點采用二階迎風插值,因此碰撞后粒子的分布函數為

式中,ak,bl,sk,dl(k,l=?1,0,1,2)為插值系數,md=sgn(dξ),nd=sgn(dη)用來確定插值過程中插值點的取值方向.對于與邊界相鄰的網格節點,由于網格節點數目的限制,無法使用二階迎風插值,本工作采用二階中心插值求解,碰撞后粒子的分布函數為

3 邊界條件

與傳統的計算流體力學方法相同,在使用格子Boltzmann方法進行流場模擬分析時,邊界的處理對計算結果的精度、計算的數值穩定性以及計算效率都有很大影響.研究表明,非平衡態外推格式[12]在時間和空間上都具有二階精度,而且該格式操作簡單、數值穩定且適用范圍廣.因此,本工作采用非平衡態外推格式作為計算過程中的邊界處理方法.實施過程為

4 當地時間步改進技術

LBM作為一種介觀模擬方法,與采用相同數量網格數的傳統NS方程求解器相比,其計算時間要長的多,這是因為其時間推進方法為顯示格式,會受到穩定性條件的限制.只有當時間步取整個流場允許的最小值時,才能滿足時間精確度.如果被用來處理宏觀問題,即較大雷諾數時,根據公式Re=U0L/υ可知,要么增大特征長度L,要么減小運動黏度υ.當L固定,只能使υ取非常小的值.由式(9)可知,LBM控制方程碰撞項線性化處理后,使得松弛時間τ取定值,如果υ值變小,則時間步也應變得更小,這就引起了數值計算的不穩定問題.

為保證數值迭代的穩定性,時間步長?t需滿足柯朗-弗里德里希斯-列維(Courant-Friedrichs-Lewy,CFL)條件,目的首先是為了保證數值精度,其次在計算區域內粒子每次的遷移距離不會超越一個格子.參照傳統計算流體力學中的處理方法,全局性時間步為

式中:λ為CFL數,0<λ 6 1;min{1/c(i,j)}遍取計算平面所有的格點上離散速度的值.本工作引入當地時間步法,各點都可以采用該點穩定條件允許的最大時間步長,這樣流場處處以穩定極限向前發展,計算速度明顯提高.一般認為,引入當地時間步法可使計算時間減少一個量級.引入當地時間步法的時間步為

通過算例驗證可知,在低雷諾數下引入當地時間步法可以顯著提高計算效率,但是依然不能改善高雷諾數流動下計算不穩定的問題.為了解決這一問題,本工作結合分子運動論的相關原理,引入3個不同的時間尺度:碰撞時間?t1=d/ξ、平均自由時間?t2=l/ξ和宏觀時間?t3=L/ξ,其中d,l,L分別是氣體分子的大小、平均自由程和宏觀量變化尺度,ξ是氣體分子平均運動速率.假設為考慮平均自由程的碰撞時間,則有

式中,c1=0.5.理論上,平均自由時間應等于分子平均自由程與分子運動特征速度之比.但在實際計算中,由于網格分辨率有限,遷移時間?t2可以按式(19)來計算,式(21)即為只考慮兩粒子接觸時的碰撞時間.通過氣體分子平均運動速率ξ,可得出碰撞時間?t1(i,j)和平均自由時間?t2(i,j)的關系為

式中,c2=l/d>1.由于分布函數fi在?t1(i,j)的時間間隔中變化很小,fi變化的時間尺度為?t01(i,j).因此,為改善計算的收斂情況,可令時間步長隨網格點不同自動調節,定義局部時間步長為

5 算例分析

NACA0012可認為是最簡單的對稱低速模型,近幾十年來人們對它進行了大量的數值模擬[13]和風洞實驗[14].本工作選擇翼型繞流作為數值模擬的對象,主要目的有:①將采用當地時間步法的計算結果與試驗值進行對比,檢驗算法在高雷諾數下的可靠性;②在不同雷諾數下對翼型繞流進行數值模擬,驗證本工作中算法保持數值模擬穩定性的能力,同時檢驗算法的計算效率.

5.1 可靠性分析

為了檢驗算法的可靠性,本工作采用了Harris[15]的實驗設置和實驗數據,該實驗是1981年在美國NASA Langley Research Center實驗室3.0 m×9.0 m的跨聲速風洞中進行的,實驗來流馬赫數Ma=0.3,模型弦長為0.635 m,雷諾數Re=3×106.Harris公布的實驗數據含有不同攻角下翼型的升力系數值和翼型表面的壓力系數,因此,可以將本工作計算的結果與Harris提供的實驗數據進行對比.

計算模型采用304×90的C型網格,翼型上有128個點,最小網格(附面層第一層網格)尺度為8.6×10?5,翼型弦長為1.NACA0012流場計算網格如圖1所示,其中x,y等相關物理量均為無量綱量,進口邊界設在離翼型前緣10倍弦長處,出口邊界為15倍弦長處.

表1是NACA0012翼型在不同攻角下的升力系數值.將本工作的計算結果和實驗結果進行比較,可以看出:升力系數的變化規律是一致的,均隨著攻角的增加而變大;計算值要比實驗值偏高,這與本工作在計算中沒有考慮湍流模型有一定關系.

圖2為NACA0012翼型在?0.14?,3.86?和6.86?攻角下吸力面和壓力面壓力系數分布曲線.由翼型表面壓力系數Cp的分布可知,翼型上表面流動并沒有發生流動分離,翼型下表面的壓力系數與實驗值吻合良好,而上表面的壓力系數隨著攻角增大會出現一定的誤差.在未采用湍流模型的情況下,這樣的計算精度是可以接受的.

圖1 NACA0012翼型流場的計算網格Fig.1 Calculation grid of flow around NACA0012

表1 NACA0012翼型不同攻角下的升力系數Table 1 Lift coefficient of NACA0012 airfoil at different angles of attack

圖2 NACA0012翼型不同攻角下壓力系數分布Fig.2 Pressure coefficient distribution of NACA0012 airfoil at different angles of attack

圖3為8.86?攻角時NACA0012翼型的流線圖.由圖可知,翼型繞流為附著流，在翼型上表面附近出現明顯的附著渦結構,但未出現分離渦結構,但可以預見,隨著翼型攻角的增大,附著渦會逐漸地變大,最后形成穩定的脫落渦結構.

圖3 NACA0012翼型8.86?攻角時的流線圖Fig.3 Streamline diagram of NACA0012 airfoil at 8.86?angle of attack

表2 雷諾數和CFL數不同時NACA翼型的計算狀態參數Table 2 Calculation of state parameters for different Reynolds numbers and CFL numbers

5.2 計算效率和穩定性分析

為了進一步檢驗算法在高雷諾數下的適用性、計算效率和可靠性,本工作從算法改進前后計算狀態均良好的中雷諾數出發,對不同中高雷諾數下的NACA0012翼型繞流進行了數值模擬,旨在發現雷諾數在增大至高雷諾數過程中計算狀態的變化規律.取來流馬赫數Ma=0.1,迎角α=7?,雷諾數Re分別取5.0×104,5.0×105,2.5×106和5.0×106,通過與采用當地時間步法的計算結果進行對比,分析檢驗算法在高雷諾數下的計算效率和穩定性.另外,將雷諾數為5.0×105時的計算結果與CFL3D[4]值進行對比,檢驗算法在高雷諾數下的可靠性.

表2給出了NACA0012翼型在不同雷諾數和CFL數時的計算狀態和計算耗時情況,以及相同雷諾數下LBM與采用當地時間步法后計算狀態的比較.本工作測得在CFL數為1.0時,滿足計算穩定所需的最大雷諾數為5.0×104.因此,在此基礎上進行比較分析.由式(9)可知,在松弛時間τ取定值時,如果υ值隨著雷諾數Re變小,則時間步也應同比例的變小,因此在雷諾數增大時,通過控制CFL數改變時間步的大小,可使松弛時間τ保持不變.觀察表2可知,隨著雷諾數Re的增大,未采用當地時間步法的LBM計算狀態發生了變化.當雷諾數逐漸增大,計算狀態變得不穩定;當雷諾數繼續增大時,數值計算趨于發散.圖4是雷諾數為2.5×106時計算的殘差收斂曲線,其中R為殘差,T為迭代次數.由圖可以看出:未采用當地時間步法的計算在最初階段會出現震蕩,最后經過一段時間后出現發散;而采用當地時間步法后,計算快速的趨于穩定.另外,采用當地時間步法后,高雷諾數下計算狀態不穩定的現象得到了改善,CFL數作為局部時間步的調控參數,已不再受雷諾數的約束.在計算耗時方面,隨著雷諾數的增大,計算耗時變化規律一致,但與未采用當地時間步法相比,采用當地時間步法可以縮短計算趨于穩定的時間,明顯提高了計算效率.

圖4 計算殘差收斂曲線(Re=2.5×106,α=7?)Fig.4 Calculated error convergence curve(Re=2.5× 106,α =7?)

圖5 為翼型在雷諾數為5.0×104和5.0×105時表面的壓力系數分布曲線.由圖5(a)可以看出,兩種計算結果有稍許不同,但基本吻合,表明是否采用當地時間步法對當前雷諾數下的結果影響不大.由圖5(b)可以看出,二者的計算結果卻發生了明顯的變化.采用當地時間步法后計算的結果與CFL3D值吻合較好,而未采用該方法計算的結果偏差較大,特別是翼型尾緣處的壓力系數分布.這是因為隨著雷諾數增大,時間步變得更小,造成了計算的不穩定現象.

圖5 NACA0012翼型表面壓力系數分布曲線(α=7?)Fig.5 Pressure coefficient distribution of NACA0012 airfoil surface(α =7?)

未采用當地時間步法時,NACA0012翼型的流線如圖6所示.由圖可知,翼型上表面靠近尾緣處形成了明顯的附著渦結構,這使得該位置的壓力系數明顯偏大.

圖 6 NACA0012翼型表面的流線圖(Re=5.0×105,α=7?)Fig.6 Airfoil flow chart of NACA0012 airfoil surface(Re=5.0×105,α =7?)

圖7 為雷諾數為2.5×106和5.0×106時翼型表面的壓力系數曲線.由圖7可知,曲線均沒有出現大的變形,計算結果良好.

圖7 NACA0012翼型表面壓力系數分布曲線(α=7?)Fig.7 Streamline diagram of NACA0012 airfoil(α =7?)

6 結束語

高雷諾數流動下基于LBM方法的數值模擬,在理論和應用上都是十分重要的,同時它也是一項極具挑戰性的研究課題.本工作在未使用湍流模型情況下研究了NACA0012低速繞流問題,并分別對算法可靠性、計算效率和穩定性進行了分析,通過與實驗值進行比較,發現結果雖有所差異,但仍可接受.

研究結果表明,本工作結合分子運動論的相關原理,采用了在貼體網格GILBM框架下發展的當地時間步改進技術,對NACA0012翼型不同工況進行了數值模擬,該方法使得計算區域內的各點都可以采用滿足CFL條件所允許的最大時間步長.因此,流場處處以穩定極限向前發展,這不僅改善了算法的計算穩定性,又提高了計算效率,從而增強了LBM方法在處理高雷諾數流動問題方面的適應性.本工作的算例計算結果與實驗結果十分吻合,證明了該方法的準確性.