合成進口湍流方法在大渦模擬中的應(yīng)用

劉俊杰，仲冬冬，葛寧

(南京航空航天大學 能源與動力學院，江蘇 南京 210016)

0 引言

與工程中常用的雷諾平均N-S方程相比，大渦模擬(large eddy simulation,LES)在理論上可以更準確地捕捉流場細節(jié)，對于分離和其他流動的模擬具有明顯的優(yōu)勢。隨著計算機硬件水平的提高，工程領(lǐng)域逐漸采用大渦模擬。與雷諾平均方法不同，對于湍流的大渦模擬，有必要在計算域的入口處提供最接近真實的瞬時湍流速度。研究表明，湍流的入口條件對平面射流、空間擴展的邊界層流和后臺階流動有重要影響[1]。為了了解進口瞬時擾動速度設(shè)定方法的特點，研究大渦模擬對促進大渦模擬的工程應(yīng)用具有重要意義。

為解決真實模擬的問題，KEATING A[1]、SAGAUT P等[2]發(fā)表的文獻提出了許多解決辦法。大部分都使用平均速度分解，指定了平均速度剖面，然后試圖在平均速度上疊加真實的湍流脈動，這些方法叫做合成渦方法(synthesis eddy method)。這一系列方法的期望輸出是一個非定常湍流流場，所呈現(xiàn)的流場信息為低階統(tǒng)計量均值和均方根速度以及與之理想的兩點時空相關(guān)性。低階統(tǒng)計量均值較易匹配，而均方根速度以及與之理想的兩點時空相關(guān)性相比較而言可能更難重建，因為它包括尺度和結(jié)構(gòu)方面的信息。然而，這些重構(gòu)是必要的，用以模仿真實物理機理存在于實際的壁面約束流動。通常來說，入流數(shù)據(jù)不具有這種性質(zhì)，流動必須沿一定距離進行重構(gòu)調(diào)整，通常用初始邊界層厚度δ0來無量綱化。在δ0以下稱為“自適應(yīng)距離”，該距離決定了方法的性能，因為它確定了用于生成真實湍流計算域的流向長度[3]。

根據(jù)參考文獻[1]，入流方法可分為3類：回收調(diào)節(jié)法、前體數(shù)據(jù)庫法和合成湍流法。基本上，回收調(diào)節(jié)法使用周期邊界條件的變形，考慮到流體各向異性。在平板邊界層的情況下，回收調(diào)節(jié)法的核心思想是根據(jù)湍流邊界層的相似理論，將邊界層下游回收平面的時均參數(shù)和脈動參數(shù)提取出來，根據(jù)邊界層相似理論對參數(shù)進行調(diào)整，然后疊加到入口平面上，使入口湍流邊界層保持合理的湍流脈動信息。這種方法的特別之處在于通過回收過程實現(xiàn)湍流脈動場建立在計算域的內(nèi)部，從而避免了高額的計算負擔。

前體數(shù)據(jù)庫這種技術(shù)已被廣泛應(yīng)用，并在參考文獻[3-4]中進行了修正，用以模擬大范圍的流動。生成流入數(shù)據(jù)的更通用的方法是從單獨的前體數(shù)據(jù)庫中提取數(shù)據(jù)。通過提供適當?shù)乃俣炔▌颖壤?，?shù)據(jù)庫可以為各種配置和雷諾數(shù)提供數(shù)據(jù)。盡管如此，這些方法仍需要很大的存儲容量。合成湍流法是收集所有合成湍流波動的方法，大多數(shù)情況下是通過隨機序列進行的。然后，目標是修改這些序列，以使它們在數(shù)值或黏性耗散下不消失，并且其統(tǒng)計屬性與實際流程的統(tǒng)計屬性接近。LUND等使用一種基本的隨機波動方法，根據(jù)該方法將目標雷諾應(yīng)力分配給白噪聲并引入作為入口。但該方法生成的數(shù)據(jù)既缺乏湍流結(jié)構(gòu)又缺乏非線性能量傳遞。

本文提出了一種應(yīng)用在大渦模擬計算中產(chǎn)生入流湍流條件的方法。這是JARRIN等人合成渦一種變形，是使用由雷諾應(yīng)力張量的Cholesky分解的2階矩來添加隨機速度信號，信號由具有規(guī)定幾何形狀、隨機符號和位置的湍流結(jié)構(gòu)疊加而成。該方法對隨機信號的定義進行了修正，使其可以分為幾種模式，具有不同的時間、長度和速度尺度，也具有不同的渦度，中心思想是更真實地再現(xiàn)在湍流邊界層流動壁面法線方向上的尺度分布。本文采用基于k-ωSST湍流模型，結(jié)合 WENO_ZQ高精度格式對平板同時模擬真實進口條件，獲得了精細的流場結(jié)構(gòu)，并利用渦識別方法和數(shù)值紋影圖等流動顯示方法識別流場。

1 數(shù)值方法

1.1 方法介紹

本文采用SANDBERG R D[5]提出的方法。這個便捷的方法是基于幾個離散波數(shù)的疊加。在入口的3個速度分量中引入擾動，采用多個離散波之和。

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ù)在時間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)

上述公式的每個參數(shù)都允許對入口擾動進行微調(diào)以獲得湍流目標狀態(tài)?？梢允褂妙l率β1(n)以及波數(shù)β2(n)和β3(n)來調(diào)整時間和長度刻度。計算域在螺距方向和跨度方向上是周期性的，這要求必須選擇波數(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)的入口方向來實現(xiàn)，最后可以使用ψi(n)指定相移。在每個方向上，使用系數(shù)Ii獨立地調(diào)整湍流強度。在詳盡的初步研究中發(fā)現(xiàn)，基于4種不同的波數(shù)組合j，總共使用16個波獲得了良好的結(jié)果。表1中n、φ2、φ3給出了每個j的波數(shù)。每個波數(shù)組合j中的每一個波移動四分之一周期，相應(yīng)的ψi由表1的其他參數(shù)給出。

表1 湍流生成參數(shù)

1.2 求解器介紹

本研究內(nèi)容是基于課題組自主開發(fā)的CFD求解器NUAA-Turbo2.0進行的，控制方程與數(shù)值解法詳見參考文獻[6]。求解器采用有限體積結(jié)構(gòu)化網(wǎng)格、顯式3階龍格庫塔(Runge-Kutta methods)和全隱式雙時間步LU-SGS時間推進方法，同時發(fā)展了Roe和AUSM-up通量分裂方法、MUSCL、WENO_JS和WENO_ZQ等高精度方法重構(gòu)界面量。

1.3 研究對象

選取一個經(jīng)典的平板算例進行驗證，參考Naguib實驗數(shù)據(jù)[7]。出口為靜壓，絕熱壁面為非滑移邊界，周向為遠場邊界，跨度方向為周期邊界條件。計算域和網(wǎng)格分布如圖1所示。網(wǎng)格總量為420萬。滿足網(wǎng)格無關(guān)性要求。以進口邊界層厚度δ作為計算參考長度。計算域為8.62δ×2.15δ×5.23δ[8]，與實驗?zāi)Ｐ鸵恢?，對?yīng)流動方向、展向方向和法向分別是網(wǎng)格在流向x和展向方向y上均勻分布；沿法線方向z在靠近壁面處進行加密處理，保證靠近墻的第一層高度z+≤1。頻率值[9]模擬采用MPI技術(shù)，沿流向分成10個塊進行并行計算。

圖1 平板算例圖

1.4 計算設(shè)置

計算采用基于SST 湍流模型的RANS方法進行模擬，然后在定常條件收斂下?lián)Q用LES方法模擬進行非定常計算。同一研究對象采用同一套結(jié)構(gòu)化網(wǎng)格，采用6階中心差分格式處理黏性項，除RANS采用WENO 3階格式外，其余都采用具有5階精度的WENO_ZQ格式重構(gòu)界面量。選取的物理時間步長可保證尾跡區(qū)CFL數(shù)≤1。

2 結(jié)果分析和討論

湍流邊界層(TBL)是自然界中廣泛存在的一種復(fù)雜的流動結(jié)構(gòu)?，F(xiàn)有的湍流模型和流動機理主要基于不可壓縮的湍流邊界層，例如邊界層分層理論，發(fā)夾渦理論和高低速條狀結(jié)構(gòu)。在湍流邊界層中，發(fā)現(xiàn)湍流邊界層沒有呈現(xiàn)無序狀態(tài)，但是相干渦結(jié)構(gòu)維持了邊界層的內(nèi)部運動。根據(jù)經(jīng)典的發(fā)夾式渦旋理論，當兩個相鄰的反向旋渦在下游發(fā)展時，將在旋渦頭之間形成沿展向連接的渦旋結(jié)構(gòu)，并最終發(fā)展成完整的發(fā)夾式渦旋結(jié)構(gòu)。WU X H[10]在零壓力梯度下對不可壓縮湍流邊界層進行了直接數(shù)值模擬，結(jié)果證實了發(fā)夾渦旋的存在，并且湍流邊界層中的大多數(shù)發(fā)夾旋渦幾乎是對稱的發(fā)夾結(jié)構(gòu)。

為了便于觀察，用圖2顯示了通過Q準則識別的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)：當動量厚度雷諾數(shù)(基于邊界層動量厚度的雷諾數(shù))<500時，總體渦旋結(jié)構(gòu)較短，呈馬蹄形渦旋或渦旋環(huán)。當動量厚度雷諾數(shù)>2 000時，渦旋結(jié)構(gòu)呈現(xiàn)細長的發(fā)夾渦。結(jié)果表明，隨著雷諾數(shù)的增加，發(fā)夾渦破裂并呈現(xiàn)出不對稱的藤條結(jié)構(gòu)。

圖2 Q準則識別的三維擬序結(jié)構(gòu)

圖3顯示了不同的雷諾數(shù)對發(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)被認為與湍流邊界層的維持和發(fā)展密切相關(guān)。近年來，大量的數(shù)值研究證實了這種結(jié)構(gòu)不僅存在于不可壓縮流動中，而且存在于可壓縮湍流邊界層中。

圖4 邊界層底部的速度帶結(jié)構(gòu)

圖5顯示了流向方向在距進口截面約為8δ處沿法線方向的平均速度剖面。黑色實線表示此LES計算的結(jié)果。MORKOVIN M V[14]提出可壓縮流邊界層中的參數(shù)分布可以通過數(shù)學關(guān)系轉(zhuǎn)換與不可壓縮流中的參數(shù)分布聯(lián)系起來。變換后的速度分布仍然滿足經(jīng)典的壁面定律分布，例如線性定律、對數(shù)定律和尾跡定律，這就是著名的Morkovin假設(shè)。

圖5 計算的結(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)

3 結(jié)語

對于湍流邊界層的數(shù)值模擬，主要結(jié)論如下：在低聲速湍流邊界層中，流場中的擬序渦總體上表現(xiàn)出較短的整體結(jié)構(gòu)，表現(xiàn)為馬蹄形渦或渦旋環(huán)，這種結(jié)構(gòu)取決于傳入流的動量厚度和雷諾數(shù)。密度加權(quán)轉(zhuǎn)換后的平均速度分布仍遵循經(jīng)典的壁面律，例如線性律、對數(shù)律和尾跡律。