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

陳林烽

(江蘇科技大學,江蘇鎮(zhèn)江 212003)

引言

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

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

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

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

1 數(shù)值方法

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

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

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

2 多尺度方法

2.1 多尺度分解

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

相應地,可以將權(quán)函數(shù)分解為大尺度和小尺度

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

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

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

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

對比于傳統(tǒng)大渦模擬方法的亞格子應力包含的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) 右邊為大尺度相關(guān)項.式中?·和N()分別為連續(xù)性和Navier-Stokes 大尺度方程殘差,即

于是,式(15)可簡化為

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

其等效于

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

式中c1,c2和c3為常數(shù),G為單元局部坐標與全局坐標之間的轉(zhuǎn)換系數(shù)矩陣

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

于是解得

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

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

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

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

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

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

3 時間推進方法

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

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

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

式中

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

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

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

4 數(shù)值驗證

4.1 算例設(shè)置

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

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

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

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

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

表1 計算參數(shù)Table 1 Computation parameters

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

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

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

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

4.2 結(jié)果與分析

在流場穩(wěn)定以后,對20 000 個時間步的流場數(shù)據(jù)進行了統(tǒng)計計算,得到流場的平均速度和脈動速度均方根值.

圖6 為流向平均速度在垂向方向(y+)的分布圖,〈U〉 表示在x和z方向上的空間及時間綜合平均統(tǒng)計量,y+=uτy/ν=yReτ.圖中分別給出了Jimenez的直接數(shù)值模擬和32×32×32,48×48×48 個單元的Smagorinsky 模型大渦模擬有限差分[3-4]和本文的計算結(jié)果.可以看出,在壁面附近(y+＜ 15),本文所用方法的計算結(jié)果與DNS 數(shù)據(jù)匹配的很好,而Smagorinsky 模型過大的耗散導致該區(qū)域的計算結(jié)果低于DNS 數(shù)據(jù); 在對數(shù)率區(qū)(30 ＜y+＜50) 和黏性外層(y+＞50),本文的計算結(jié)果略大于DNS 計算結(jié)果.該區(qū)域的平均速度隨著網(wǎng)格數(shù)量增多變得更接近DNS 數(shù)據(jù).根據(jù)這一結(jié)果,可以判斷該模型在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截面上的空間點數(shù),N為時間采樣數(shù).與DNS 結(jié)果相比發(fā)現(xiàn)本文所用大渦數(shù)值模擬的流向脈動速度rms 在壁面附近(y+＜15)DNS 結(jié)果匹配良好,然而在對數(shù)率區(qū)(30 ＜y+＜50)和黏性外層(y+＞50)則略大于DNS的結(jié)果; Smagorinsky 模型的計算結(jié)果顯示流向脈動速度均方根的峰值跟DNS 結(jié)果有較大偏差.在垂向與展向上,兩種大渦模擬的計算結(jié)果均小于DNS 結(jié)果; 相比于Smagorinsky 模型的計算結(jié)果,本文所得到的結(jié)果要更接近DNS 數(shù)據(jù).流向脈動速度均方根為湍動能的主要來源,因此它的增加意味著流場整體的湍動能變大,盡管其他兩個方向上的脈動速度均方根略有減小,即表明該模型所提供的湍動能耗散略偏小,從而導致圖6 中的平均流向速度偏大.

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

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

圖9 為DNS 數(shù)據(jù)和48×48×48 網(wǎng)格計算得到的8 ＜U＜10(U為流向速度)的等值面分布圖.如圖所示,發(fā)現(xiàn)DNS 數(shù)據(jù)得到的等值面圖上可以看到很多細小的流場結(jié)構(gòu),而大渦數(shù)值模擬計算得到的流場僅存留了較大尺度的旋渦結(jié)構(gòu).

圖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 網(wǎng)格計算數(shù)據(jù)在壁面附近y+=5.4 處xz截面的流向速度云圖.從圖中可以看出,在近壁面的流場,出現(xiàn)了明顯的低速條帶結(jié)構(gòu).近壁面處出現(xiàn)低速條帶結(jié)構(gòu)為觸發(fā)湍流擬序結(jié)構(gòu)的第一號信息.這說明,盡管本文在模擬中使用的網(wǎng)格數(shù)較小,該模型還是能有效觸發(fā)近壁湍流的擬序結(jié)構(gòu).

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

5 結(jié)論

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

文中使用新提出的大渦模擬方法的自編程序?qū)崿F(xiàn)了槽道湍流的并行數(shù)值計算.數(shù)值計算結(jié)果驗證了該方法的合理性以及自編程序的收斂性.并得出,該方法可以在使用較少的網(wǎng)格數(shù)可以得到與DNS 結(jié)果接近的數(shù)值解,并在近壁區(qū)可以觀察到明顯的低速條帶結(jié)構(gòu).另外,數(shù)值結(jié)果證實該方法在網(wǎng)格數(shù)較少的情況下流向向另外兩個方向的湍動能輸運略偏低.