基于Navier-Stokes 方程殘差的隱式大渦模擬有限元模型1)

陳林烽

(江蘇科技大學,江蘇鎮江 212003)

引言

作為湍流數值模擬方法的一種,大渦模擬方法(large eddy simulation,LES) 在未來的工業應用中具有極大的潛力.大渦模擬方法的主體思想是針對湍流不同尺度渦的不同特點,弓入介于大渦尺度和小渦尺度之間的濾波器對Navier-Stokes 方程進行過濾,保留濾波器寬度以上的大渦特征,對所保留的大尺度渦結構進行直接模擬,對小尺度渦結構采用不可解尺度模型模擬[1-2].相比于直接數值模擬(Direct numerical simulation,DNS),它過濾了小尺度信息,計算時網格數量需求更低,計算耗時變得更短;相比于雷諾平均方法(Reynolds average Navier-stokes,RANS),它通過調節網格尺寸,可以保留不同程度的大尺度信息,因此除了平均尺度之外,它還能計算得出不同程度的大尺度旋渦[3-4].在過去的幾十年里,大渦模擬方法的研究一直是計算流體力學科研學者們的重點關注問題[5-6].

大渦模擬方法的核心問題是構建不可解尺度(不可解尺度) 模型.自大渦模擬方法提出以來,已經有多種不可解尺度模型被提出.Smagorinsky[7]在1963 年提出了渦黏形式的不可解尺度模型,隨后Lilly[8]利用湍動能譜確定了Smagorinsky 模型系數.該模型最早應用于大氣工程中的大渦模擬計算.實際使用過程中發現該模型在近壁區耗散過大,無法捕捉層流到湍流的轉捩過程.Bardina 等[9]提出流動中可解尺度脈動向不可解尺度脈動輸運的動量由可解尺度脈動中的最小尺度部分產生,并且可解尺度脈動的最小尺度脈動速度和過濾掉的小尺度脈動速度相似,從而提出尺度相似模式.該模型在實際應用時湍動能耗散太小,往往導致計算發散.Germano 等[10]提出二次過濾并假設二次過濾后的亞格子應力等于粗、細網格上的亞格子應力差,從而得出動態不可解尺度模式.動態確定模型系數時,需要對流場信息進行統計平均.對于湍流脈動存在統計均勻方法,通常采用空間統計平均來代替系綜平均.對于沒有統計均勻方向的流動,Maneveau 等[11-12]提出沿質點軌跡平均的動態確定模式系數方法.實際計算結果表明,對于規則幾何邊界處的湍流,動態不可解尺度模式可以得到較好的計算結果.對于不規則幾何邊界處的湍流,二次過濾計算動態不可解尺度模型系數時往往會出現錯誤.隨后,Nicoud 等[13]提出了壁面自適應不可解尺度模型(WALE).

隨著大渦模擬方法的發展,Hughes 等[14]在有限元方法的框架下將多尺度方法應用于大渦模擬方法.在有限元的框架下,假設Navier-Stokes 方程的解的形函數由可解尺度和不可解尺度形函數疊加組成,弓入對應的權函數,將Navier-Stokes 方程的有限元變分形式分解為可解尺度和不可解尺度系統.Hughes 等[14]將可解尺度進一步分解為大尺度與小尺度,采用傳統的Smagorinsky 不可解尺度模型模擬不可解尺度,然后將不可解尺度模型僅作用于小尺度(或大尺度).該方法無需使用過濾函數對Navier-Stokes 方程的變量進行過濾,通過對形函數進行尺度分解實現解的尺度分解.不可解尺度模型僅作用于可解尺度中的小尺度(或大尺度),更貼近湍動能輸運規律.Hughes等[15-16]通過各向同性湍流和槽道湍流的數值實驗驗證了基于有限元多尺度概念的大渦模擬方法較傳統的大渦模擬方法更為準確,證實了該方法的未來潛力.

Hughes 等[17]在將多尺度方法運用于對流擴散方程和Stokes 方程時,借助Green 函數求出了不可解尺度系統的解析解.通過將不可解尺度的解析解代入到可解尺度系統后,可以使用較少的網格數量計算得到理想的數值模擬結果.根據不可解尺度的解析解形式,可以發現不可解尺度是可解尺度方程殘余量的函數.對于Navier-Stokes 方程,鑒于方程的復雜性,無法通過不可解尺度系統求得解析形式.在Hughes 等[18-22]的工作的啟發下,本文將根據Navier-Stokes 方程的不可解尺度系統建立基于大尺度方程殘差的不可解尺度模型,然后將其代入Navier-Stokes方程的大尺度系統并直接求解之得到Navier-Stokes方程的大尺度解.

1 數值方法

本文將采用有限元方法對流體控制方程Navier-Stokes 方程

進行離散求解[23].令W={w,q} 為對應于Navier-Stokes 變量U={u,p}的權函數,則Navier-Stokes 方程的有限元變分形式可寫為

式中(,)Ω表示前后兩個變量乘積在計算域內部的有限元單元內的積分,Ω 表示計算域內部空間.

2 多尺度方法

2.1 多尺度分解

鑒于湍流流動內存在多種尺度的旋渦,可以將Navier-Stokes 方程的解看成是大尺度解和小尺度解疊加而成(如圖1),即

相應地,可以將權函數分解為大尺度和小尺度

將式(6)和式(7)代入Navier-Stokes 方程的變分形式(3)中可以得到分別投影到大尺度權函數和小尺度權函數的兩組方程[24-27]

圖1 Navier-Stokes 方程解的尺度分解示意圖Fig.1 Sketch of scale decomposition of solution to Navier-Stokes equations

這里,將式(8)稱為大尺度系統,將式(9)稱為小尺度系統.其中

式中第二、三和四行為小尺度相關項.通常情況,?u′/?t和ν?2u′為較小值,于是式(10) 第二行的前兩項可忽略[28].為了避免在方程中出現,?u′和?p′項,對式(10)中的對應項采用分布積分[26],于是

對比于傳統大渦模擬方法的亞格子應力包含的Leonard,Corss 和Reynolds 應力項,式(11)中包含有對應項[30]:

值得注意的是,在以上大尺度方程的推導過程中還沒有使用任何模型.

對于大尺度方程式(11),如果能找到{u′,p′} 的精確表達式,便可以通過求解式(11) 得到精確的大尺度解.式(11)與式(9)相互耦合,于是可以試圖通過式(9)去找到{u′,p′}表達式.鑒于式(9)與式(11)同樣復雜,以下將通過分析式(9) 的具體形式,去建立{u′,p′}的近似模型.

2.2 不可解尺度模型

將式(9)展開后可以得到

將式(12)的大尺度和小尺度項進行整理得到

注意到式(15) 和式(16) 右邊為大尺度相關項.式中?·和N()分別為連續性和Navier-Stokes 大尺度方程殘差,即

于是,式(15)可簡化為

假設?p′為極小量,式(15)可以再次簡化為

其等效于

忽略u′·?()的影響,小尺度u′的方程最終為

式中c1,c2和c3為常數,G為單元局部坐標與全局坐標之間的轉換系數矩陣

ξ為單元局部坐標,x為全局坐標,G:G表示矩陣G的雙點積.將式(22)代入式(21)得

于是解得

由式(19)和式(16)并忽略式(19)中u′·?(′)的影響可以得到

式(32)兩邊同時點乘gτm得到

將式(33)代入式(34),在計算p′時假設rm已經滿足rm=0,于是得到p′的模型為

令τc=(gτm·g)-1,則

將式(27)和式(36)代入大尺度方程(11)可以得到

注意到,此時大尺度方程僅包含大尺度變量,即方程得到封閉.

3 時間推進方法

在數值模擬時,將直接去求解大尺度方程式(37).這里采用Jansen 等[34]提出的generalized-α 方法進行時間步推進.以及將介紹generalized-α 方法應用在式(37)的詳細細節.

將a定義為{}的節點系數,˙a定義為a的時間倒數.于是,式(37)離散后可整理為

式中A和B分別為式(37)離散系統整理后和a的系數矩陣.弓入generalized-α 得到

式中

其中ζ=0.5 以確保時間推進精度為二階精度和無條件收斂

在每一個時間步上采用迭代方法求解式(38).假設

將式(45)和式(46)代回式(41)和式(42),再將得到的和an+αf代入式(39),并對其使用牛頓迭代

4 數值驗證

4.1 算例設置

本文針對Reτ=180 的槽道湍流問題,對基于Navier-Stokes 大尺度方程殘差的大渦數值模擬方法及其自編程序[35-37]做出驗證.其中Reτ=uτδ/ν,uτ為壁面摩擦速度,δ 為槽道的半寬度,ν 為流體運動黏性系數.流場的計算域為Lx=2δπ,Ly=2δ,Lz=4δπ/3(如圖2 所示).

圖2 槽道流示意圖Fig.2 Sketch of channel flow

計算域的網格劃分采用六面體網格(如圖3 所示),在流向(x)與展向(z)方向上采用均勻網格,在垂向(y)方向采用雙向正切函數

對壁面附近網格進行加密.本文分別使用32×32×32和48×48×48 個單元的兩種網格對流場進行數值模擬.網格的詳細參數見表1,表示壁面單元的無量綱尺寸.

圖3 48×48×48 個六面體單元組成的網格示意圖Fig.3 Mesh sketch of 48×48×48 hexahedral elements

表1 計算參數Table 1 Computation parameters

三維單元上的形函數以及有限元權函數均使用雙線性(bilinear) 函數,即形函數由每一個方向上的線性函數相乘得到,如(1-ξ)(1-ψ)(1-ζ).圖4 為垂向上由雙向正切函數得到的非均勻網格單元上的線性形函數示意圖.

圖4 垂向線性形函數示意圖Fig.4 Sketch of linear shape functions in the normal direction

數值模擬中,流向和展向采用周期性邊界條件,上、下固壁處使用無滑移邊界條件(uwall=0).流向進出口之間給定一個恒定的壓力差來維持流體流動,并在入口處給定一個壓力參考值作為壓力的邊界條件.對離散的式(37)使用gernerzlized-α 后,可以得到式(47)的線性方程系統.這里,采用基于ILU 的GMRES 迭代算法并行求解式(47)的線性方程系統.圖5為使用48×48×48 個單元計算得到的瞬時流場流向速度云圖.從入口的截面可以看到,在壁面附近有很明顯的展向旋渦.

圖5 48×48×48 個單元瞬時流場流向速度云圖Fig.5 Instantaneous contours of streamwise velocity using 48×48×48 elements

4.2 結果與分析

在流場穩定以后,對20 000 個時間步的流場數據進行了統計計算,得到流場的平均速度和脈動速度均方根值.

圖6 為流向平均速度在垂向方向(y+)的分布圖,〈U〉 表示在x和z方向上的空間及時間綜合平均統計量,y+=uτy/ν=yReτ.圖中分別給出了Jimenez的直接數值模擬和32×32×32,48×48×48 個單元的Smagorinsky 模型大渦模擬有限差分[3-4]和本文的計算結果.可以看出,在壁面附近(y+＜ 15),本文所用方法的計算結果與DNS 數據匹配的很好,而Smagorinsky 模型過大的耗散導致該區域的計算結果低于DNS 數據; 在對數率區(30 ＜y+＜50) 和黏性外層(y+＞50),本文的計算結果略大于DNS 計算結果.該區域的平均速度隨著網格數量增多變得更接近DNS 數據.根據這一結果,可以判斷該模型在32×32×32,48×48×48 個單元的計算中所提供的湍動能耗散略偏小.

圖6 流向平均速度在垂向方向上的分布Fig.6 Mean streamwise velocity(〈U〉)profiles against y+

圖7 給出了脈動速度均方根值(urms) 在垂向方向的分布圖,

圖7 脈動速度均方根值在垂向方向上的分布Fig.7 Root-mean-square of velocity fluctuations against y+

其中u表示任一方向的流場速度,M為xz截面上的空間點數,N為時間采樣數.與DNS 結果相比發現本文所用大渦數值模擬的流向脈動速度rms 在壁面附近(y+＜15)DNS 結果匹配良好,然而在對數率區(30 ＜y+＜50)和黏性外層(y+＞50)則略大于DNS的結果; Smagorinsky 模型的計算結果顯示流向脈動速度均方根的峰值跟DNS 結果有較大偏差.在垂向與展向上,兩種大渦模擬的計算結果均小于DNS 結果; 相比于Smagorinsky 模型的計算結果,本文所得到的結果要更接近DNS 數據.流向脈動速度均方根為湍動能的主要來源,因此它的增加意味著流場整體的湍動能變大,盡管其他兩個方向上的脈動速度均方根略有減小,即表明該模型所提供的湍動能耗散略偏小,從而導致圖6 中的平均流向速度偏大.

以上結果說明要提高該模型在低網格數條件下的計算精確度,需要改善該模型預測流向脈動速度rms 的能力.另外,根據圖7 可以判斷,在較少網格情況下使用該模型時,從流向向垂向及展向方向上的湍動能輸運略偏低.

圖8 為統計平均得到的雷諾應力曲線分布圖.如圖所示,對于Smagorinsky 模型,由于它提供過大的湍動能耗散,導致各個方向上的速度脈動均減小,因此雷諾應力均減小.與Smagorinsky 模型相比,本文計算得到的雷諾應力略大于Smagorinsky 模型的計算結果,盡管隨著網格數的增加,雷諾應力分布有所改進,但結果還是偏小.而雷諾應力值的偏小削弱了流向方向上的湍動能向垂向方向輸運.因此,改善該模型對雷諾應力的預測可以有效提高脈動速度均方根的計算精確.

圖9 為DNS 數據和48×48×48 網格計算得到的8 ＜U＜10(U為流向速度)的等值面分布圖.如圖所示,發現DNS 數據得到的等值面圖上可以看到很多細小的流場結構,而大渦數值模擬計算得到的流場僅存留了較大尺度的旋渦結構.

圖8 平均雷諾應力項〈uv〉在垂向方向上的分布Fig.8 Reynolds stress against y+

圖9 流向速度等值面分布圖:(a)DNS;(b)LES(48×48×48)Fig.9 Isosurfaces of streamwise velocity:(a)DNS;(b)present method with 48×48×48 elements

圖10 給出了48×48×48 網格計算數據在壁面附近y+=5.4 處xz截面的流向速度云圖.從圖中可以看出,在近壁面的流場,出現了明顯的低速條帶結構.近壁面處出現低速條帶結構為觸發湍流擬序結構的第一號信息.這說明,盡管本文在模擬中使用的網格數較小,該模型還是能有效觸發近壁湍流的擬序結構.

圖10 近壁面y+=5.4 處xz 截面的流向速度云圖Fig.10 Contours of streamwise velocity in the xz plane at y+=5.4

5 結論

本文在有限元的框架下,弓入多尺度的概念,提出了基于多尺度模式的大渦數值模擬方法.論文里采用有限元形函數和權函數空間的分解來實現湍流中大尺度和不可解尺度信息的區分,并推導得到大尺度系統的有限元變分方程.本文通過對不可解尺度系統中各項的特點分析,提出了不可解尺度的近似模式,從而達到封閉大尺度方法的目的.

文中使用新提出的大渦模擬方法的自編程序實現了槽道湍流的并行數值計算.數值計算結果驗證了該方法的合理性以及自編程序的收斂性.并得出,該方法可以在使用較少的網格數可以得到與DNS 結果接近的數值解,并在近壁區可以觀察到明顯的低速條帶結構.另外,數值結果證實該方法在網格數較少的情況下流向向另外兩個方向的湍動能輸運略偏低.