劉俊杰,仲冬冬,葛寧
(南京航空航天大學(xué) 能源與動力學(xué)院,江蘇 南京 210016)
與工程中常用的雷諾平均N-S方程相比,大渦模擬(large eddy simulation,LES)在理論上可以更準(zhǔn)確地捕捉流場細(xì)節(jié),對于分離和其他流動的模擬具有明顯的優(yōu)勢。隨著計(jì)算機(jī)硬件水平的提高,工程領(lǐng)域逐漸采用大渦模擬。與雷諾平均方法不同,對于湍流的大渦模擬,有必要在計(jì)算域的入口處提供最接近真實(shí)的瞬時(shí)湍流速度。研究表明,湍流的入口條件對平面射流、空間擴(kuò)展的邊界層流和后臺階流動有重要影響[1]。為了了解進(jìn)口瞬時(shí)擾動速度設(shè)定方法的特點(diǎn),研究大渦模擬對促進(jìn)大渦模擬的工程應(yīng)用具有重要意義。
為解決真實(shí)模擬的問題,KEATING A[1]、SAGAUT P等[2]發(fā)表的文獻(xiàn)提出了許多解決辦法。大部分都使用平均速度分解,指定了平均速度剖面,然后試圖在平均速度上疊加真實(shí)的湍流脈動,這些方法叫做合成渦方法(synthesis eddy method)。這一系列方法的期望輸出是一個(gè)非定常湍流流場,所呈現(xiàn)的流場信息為低階統(tǒng)計(jì)量均值和均方根速度以及與之理想的兩點(diǎn)時(shí)空相關(guān)性。低階統(tǒng)計(jì)量均值較易匹配,而均方根速度以及與之理想的兩點(diǎn)時(shí)空相關(guān)性相比較而言可能更難重建,因?yàn)樗ǔ叨群徒Y(jié)構(gòu)方面的信息。然而,這些重構(gòu)是必要的,用以模仿真實(shí)物理機(jī)理存在于實(shí)際的壁面約束流動。通常來說,入流數(shù)據(jù)不具有這種性質(zhì),流動必須沿一定距離進(jìn)行重構(gòu)調(diào)整,通常用初始邊界層厚度δ0來無量綱化。在δ0以下稱為“自適應(yīng)距離”,該距離決定了方法的性能,因?yàn)樗_定了用于生成真實(shí)湍流計(jì)算域的流向長度[3]。
根據(jù)參考文獻(xiàn)[1],入流方法可分為3類:回收調(diào)節(jié)法、前體數(shù)據(jù)庫法和合成湍流法?;旧?,回收調(diào)節(jié)法使用周期邊界條件的變形,考慮到流體各向異性。在平板邊界層的情況下,回收調(diào)節(jié)法的核心思想是根據(jù)湍流邊界層的相似理論,將邊界層下游回收平面的時(shí)均參數(shù)和脈動參數(shù)提取出來,根據(jù)邊界層相似理論對參數(shù)進(jìn)行調(diào)整,然后疊加到入口平面上,使入口湍流邊界層保持合理的湍流脈動信息。這種方法的特別之處在于通過回收過程實(shí)現(xiàn)湍流脈動場建立在計(jì)算域的內(nèi)部,從而避免了高額的計(jì)算負(fù)擔(dān)。
前體數(shù)據(jù)庫這種技術(shù)已被廣泛應(yīng)用,并在參考文獻(xiàn)[3-4]中進(jìn)行了修正,用以模擬大范圍的流動。生成流入數(shù)據(jù)的更通用的方法是從單獨(dú)的前體數(shù)據(jù)庫中提取數(shù)據(jù)。通過提供適當(dāng)?shù)乃俣炔▌颖壤瑪?shù)據(jù)庫可以為各種配置和雷諾數(shù)提供數(shù)據(jù)。盡管如此,這些方法仍需要很大的存儲容量。合成湍流法是收集所有合成湍流波動的方法,大多數(shù)情況下是通過隨機(jī)序列進(jìn)行的。然后,目標(biāo)是修改這些序列,以使它們在數(shù)值或黏性耗散下不消失,并且其統(tǒng)計(jì)屬性與實(shí)際流程的統(tǒng)計(jì)屬性接近。LUND等使用一種基本的隨機(jī)波動方法,根據(jù)該方法將目標(biāo)雷諾應(yīng)力分配給白噪聲并引入作為入口。但該方法生成的數(shù)據(jù)既缺乏湍流結(jié)構(gòu)又缺乏非線性能量傳遞。
本文提出了一種應(yīng)用在大渦模擬計(jì)算中產(chǎn)生入流湍流條件的方法。這是JARRIN等人合成渦一種變形,是使用由雷諾應(yīng)力張量的Cholesky分解的2階矩來添加隨機(jī)速度信號,信號由具有規(guī)定幾何形狀、隨機(jī)符號和位置的湍流結(jié)構(gòu)疊加而成。該方法對隨機(jī)信號的定義進(jìn)行了修正,使其可以分為幾種模式,具有不同的時(shí)間、長度和速度尺度,也具有不同的渦度,中心思想是更真實(shí)地再現(xiàn)在湍流邊界層流動壁面法線方向上的尺度分布。本文采用基于k-ωSST湍流模型,結(jié)合 WENO_ZQ高精度格式對平板同時(shí)模擬真實(shí)進(jìn)口條件,獲得了精細(xì)的流場結(jié)構(gòu),并利用渦識別方法和數(shù)值紋影圖等流動顯示方法識別流場。
本文采用SANDBERG R D[5]提出的方法。這個(gè)便捷的方法是基于幾個(gè)離散波數(shù)的疊加。在入口的3個(gè)速度分量中引入擾動,采用多個(gè)離散波之和。
u′(n,t,y,z)=I1(n)sin[f(n,t)]cos[g(t,z,n)]cos[h(t,y,n)]
v′(n,t,y,z)=I2(n)cos[f(n,t)]sin[g(t,z,n)]cos[h(t,y,n)]
w′(n,t,y,z)=I3(n)cos[f(n,t)]cos[g(t,z,n)]sin[h(t,y,n)]
(1)
諧波函數(shù)在時(shí)間t以及通過沿流向x的對流,俯仰方向z和跨度y方向的變量由下式給出:
f(t,n)=β1(n)t+φ1(n)
g(t,z,n)=β2(n)(z-φ2(n)t)+φ2(n)
(2)
h(t,y,n)=β3(n)(y-φ3(n)t)+φ3(n)
上述公式的每個(gè)參數(shù)都允許對入口擾動進(jìn)行微調(diào)以獲得湍流目標(biāo)狀態(tài)。可以使用頻率β1(n)以及波數(shù)β2(n)和β3(n)來調(diào)整時(shí)間和長度刻度。計(jì)算域在螺距方向和跨度方向上是周期性的,這要求必須選擇波數(shù)以使擾動滿足相同的約束。β2(n)和β3(n)由下式給出:
β2,3(n)=2πk2,3(n)/p2,3
(3)
其中:p2和p3分別是俯仰方向和翼展方向的長度;k2,3是整數(shù)。翼展方向和俯仰方向的波數(shù)通常不一樣。此限制不適用于β1(n)??赏ㄟ^調(diào)整不同結(jié)構(gòu)的入口方向來實(shí)現(xiàn),最后可以使用ψi(n)指定相移。在每個(gè)方向上,使用系數(shù)Ii獨(dú)立地調(diào)整湍流強(qiáng)度。在詳盡的初步研究中發(fā)現(xiàn),基于4種不同的波數(shù)組合j,總共使用16個(gè)波獲得了良好的結(jié)果。表1中n、φ2、φ3給出了每個(gè)j的波數(shù)。每個(gè)波數(shù)組合j中的每一個(gè)波移動四分之一周期,相應(yīng)的ψi由表1的其他參數(shù)給出。
表1 湍流生成參數(shù)
本研究內(nèi)容是基于課題組自主開發(fā)的CFD求解器NUAA-Turbo2.0進(jìn)行的,控制方程與數(shù)值解法詳見參考文獻(xiàn)[6]。求解器采用有限體積結(jié)構(gòu)化網(wǎng)格、顯式3階龍格庫塔(Runge-Kutta methods)和全隱式雙時(shí)間步LU-SGS時(shí)間推進(jìn)方法,同時(shí)發(fā)展了Roe和AUSM-up通量分裂方法、MUSCL、WENO_JS和WENO_ZQ等高精度方法重構(gòu)界面量。
選取一個(gè)經(jīng)典的平板算例進(jìn)行驗(yàn)證,參考Naguib實(shí)驗(yàn)數(shù)據(jù)[7]。出口為靜壓,絕熱壁面為非滑移邊界,周向?yàn)檫h(yuǎn)場邊界,跨度方向?yàn)橹芷谶吔鐥l件。計(jì)算域和網(wǎng)格分布如圖1所示。網(wǎng)格總量為420萬。滿足網(wǎng)格無關(guān)性要求。以進(jìn)口邊界層厚度δ作為計(jì)算參考長度。計(jì)算域?yàn)?.62δ×2.15δ×5.23δ[8],與實(shí)驗(yàn)?zāi)P鸵恢?,對?yīng)流動方向、展向方向和法向分別是網(wǎng)格在流向x和展向方向y上均勻分布;沿法線方向z在靠近壁面處進(jìn)行加密處理,保證靠近墻的第一層高度z+≤1。頻率值[9]模擬采用MPI技術(shù),沿流向分成10個(gè)塊進(jìn)行并行計(jì)算。
圖1 平板算例圖
計(jì)算采用基于SST 湍流模型的RANS方法進(jìn)行模擬,然后在定常條件收斂下?lián)Q用LES方法模擬進(jìn)行非定常計(jì)算。同一研究對象采用同一套結(jié)構(gòu)化網(wǎng)格,采用6階中心差分格式處理黏性項(xiàng),除RANS采用WENO 3階格式外,其余都采用具有5階精度的WENO_ZQ格式重構(gòu)界面量。選取的物理時(shí)間步長可保證尾跡區(qū)CFL數(shù)≤1。
湍流邊界層(TBL)是自然界中廣泛存在的一種復(fù)雜的流動結(jié)構(gòu)。現(xiàn)有的湍流模型和流動機(jī)理主要基于不可壓縮的湍流邊界層,例如邊界層分層理論,發(fā)夾渦理論和高低速條狀結(jié)構(gòu)。在湍流邊界層中,發(fā)現(xiàn)湍流邊界層沒有呈現(xiàn)無序狀態(tài),但是相干渦結(jié)構(gòu)維持了邊界層的內(nèi)部運(yùn)動。根據(jù)經(jīng)典的發(fā)夾式渦旋理論,當(dāng)兩個(gè)相鄰的反向旋渦在下游發(fā)展時(shí),將在旋渦頭之間形成沿展向連接的渦旋結(jié)構(gòu),并最終發(fā)展成完整的發(fā)夾式渦旋結(jié)構(gòu)。WU X H[10]在零壓力梯度下對不可壓縮湍流邊界層進(jìn)行了直接數(shù)值模擬,結(jié)果證實(shí)了發(fā)夾渦旋的存在,并且湍流邊界層中的大多數(shù)發(fā)夾旋渦幾乎是對稱的發(fā)夾結(jié)構(gòu)。
為了便于觀察,用圖2顯示了通過Q準(zhǔn)則識別的3D相干結(jié)構(gòu)。不難發(fā)現(xiàn),在湍流邊界層的外層有許多發(fā)夾狀渦流,并且大多數(shù)以不對稱結(jié)構(gòu)的形式存在。ROBINSON S K[11]指出,雷諾數(shù)對湍流邊界層中渦旋的結(jié)構(gòu)和狀態(tài)有很大的影響。HEAD M R和BANDYOPADHYAY P[12]觀察并總結(jié)了不同雷諾數(shù)下零壓力梯度湍流邊界層的相干渦結(jié)構(gòu):當(dāng)動量厚度雷諾數(shù)(基于邊界層動量厚度的雷諾數(shù))<500時(shí),總體渦旋結(jié)構(gòu)較短,呈馬蹄形渦旋或渦旋環(huán)。當(dāng)動量厚度雷諾數(shù)>2 000時(shí),渦旋結(jié)構(gòu)呈現(xiàn)細(xì)長的發(fā)夾渦。結(jié)果表明,隨著雷諾數(shù)的增加,發(fā)夾渦破裂并呈現(xiàn)出不對稱的藤條結(jié)構(gòu)。
圖2 Q準(zhǔn)則識別的三維擬序結(jié)構(gòu)
圖3顯示了不同的雷諾數(shù)對發(fā)夾狀渦旋拓?fù)涞挠绊慬13]。本文中,湍流邊界層入口動量厚度雷諾數(shù)為400,屬于低雷諾數(shù)。從圖3可以看出,相干渦旋整體結(jié)構(gòu)基本較短,呈現(xiàn)出馬蹄形渦旋或渦旋環(huán),與已有結(jié)論相符。
圖3 雷諾數(shù)對相干渦旋結(jié)構(gòu)的影響
圖4為湍流邊界層底部x+≈15處的流場云圖。從圖中可以清楚地看到沿翼展方向有高、低速交替帶。這種結(jié)構(gòu)被認(rèn)為與湍流邊界層的維持和發(fā)展密切相關(guān)。近年來,大量的數(shù)值研究證實(shí)了這種結(jié)構(gòu)不僅存在于不可壓縮流動中,而且存在于可壓縮湍流邊界層中。
圖4 邊界層底部的速度帶結(jié)構(gòu)
圖5顯示了流向方向在距進(jìn)口截面約為8δ處沿法線方向的平均速度剖面。黑色實(shí)線表示此LES計(jì)算的結(jié)果。MORKOVIN M V[14]提出可壓縮流邊界層中的參數(shù)分布可以通過數(shù)學(xué)關(guān)系轉(zhuǎn)換與不可壓縮流中的參數(shù)分布聯(lián)系起來。變換后的速度分布仍然滿足經(jīng)典的壁面定律分布,例如線性定律、對數(shù)定律和尾跡定律,這就是著名的Morkovin假設(shè)。
圖5 計(jì)算的結(jié)果與先前發(fā)表的結(jié)果之間的比較
可以看出,在黏性底層和對數(shù)區(qū)域,LES得到的曲線與經(jīng)典壁面律吻合較好。
線性律:
(4)
對數(shù)律:
(5)
以上公式中的參數(shù)由下面公式給出:
(6)
(7)
(8)
(9)
在經(jīng)典的對數(shù)律中,k為Von Karman常數(shù)[15],其值為0.4~0.41,這里取0.41;C取5.25。動力黏度μw是溫度的函數(shù),可以通過Sutherland公式確定:
(10)
對于湍流邊界層的數(shù)值模擬,主要結(jié)論如下:在低聲速湍流邊界層中,流場中的擬序渦總體上表現(xiàn)出較短的整體結(jié)構(gòu),表現(xiàn)為馬蹄形渦或渦旋環(huán),這種結(jié)構(gòu)取決于傳入流的動量厚度和雷諾數(shù)。密度加權(quán)轉(zhuǎn)換后的平均速度分布仍遵循經(jīng)典的壁面律,例如線性律、對數(shù)律和尾跡律。