劉俊杰，仲冬冬，葛寧

(南京航空航天大學(xué) 能源與動(dòng)力學(xué)院，江蘇 南京 210016)

0 引言

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

為解決真實(shí)模擬的問(wèn)題，KEATING A[1]、SAGAUT P等[2]發(fā)表的文獻(xiàn)提出了許多解決辦法。大部分都使用平均速度分解，指定了平均速度剖面，然后試圖在平均速度上疊加真實(shí)的湍流脈動(dòng)，這些方法叫做合成渦方法(synthesis eddy method)。這一系列方法的期望輸出是一個(gè)非定常湍流流場(chǎng)，所呈現(xiàn)的流場(chǎng)信息為低階統(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í)際的壁面約束流動(dòng)。通常來(lái)說(shuō)，入流數(shù)據(jù)不具有這種性質(zhì)，流動(dòng)必須沿一定距離進(jìn)行重構(gòu)調(diào)整，通常用初始邊界層厚度δ0來(lái)無(wú)量綱化。在δ0以下稱為“自適應(yīng)距離”，該距離決定了方法的性能，因?yàn)樗_定了用于生成真實(shí)湍流計(jì)算域的流向長(zhǎng)度[3]。

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

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

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

1 數(shù)值方法

1.1 方法介紹

本文采用SANDBERG R D[5]提出的方法。這個(gè)便捷的方法是基于幾個(gè)離散波數(shù)的疊加。在入口的3個(gè)速度分量中引入擾動(dòng)，采用多個(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以及通過(guò)沿流向x的對(duì)流，俯仰方向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ù)都允許對(duì)入口擾動(dòng)進(jìn)行微調(diào)以獲得湍流目標(biāo)狀態(tài)。可以使用頻率β1(n)以及波數(shù)β2(n)和β3(n)來(lái)調(diào)整時(shí)間和長(zhǎng)度刻度。計(jì)算域在螺距方向和跨度方向上是周期性的，這要求必須選擇波數(shù)以使擾動(dòng)滿足相同的約束。β2(n)和β3(n)由下式給出：

β2,3(n)=2πk2,3(n)/p2,3

(3)

其中：p2和p3分別是俯仰方向和翼展方向的長(zhǎng)度；k2,3是整數(shù)。翼展方向和俯仰方向的波數(shù)通常不一樣。此限制不適用于β1(n)?？赏ㄟ^(guò)調(diào)整不同結(jié)構(gòu)的入口方向來(lái)實(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è)波移動(dòng)四分之一周期，相應(yīng)的ψi由表1的其他參數(shù)給出。

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

1.2 求解器介紹

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

1.3 研究對(duì)象

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

圖1 平板算例圖

1.4 計(jì)算設(shè)置

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

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

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

為了便于觀察，用圖2顯示了通過(guò)Q準(zhǔn)則識(shí)別的3D相干結(jié)構(gòu)。不難發(fā)現(xiàn)，在湍流邊界層的外層有許多發(fā)夾狀渦流，并且大多數(shù)以不對(duì)稱結(jié)構(gòu)的形式存在。ROBINSON S K[11]指出，雷諾數(shù)對(duì)湍流邊界層中渦旋的結(jié)構(gòu)和狀態(tài)有很大的影響。HEAD M R和BANDYOPADHYAY P[12]觀察并總結(jié)了不同雷諾數(shù)下零壓力梯度湍流邊界層的相干渦結(jié)構(gòu)：當(dāng)動(dòng)量厚度雷諾數(shù)(基于邊界層動(dòng)量厚度的雷諾數(shù))<500時(shí)，總體渦旋結(jié)構(gòu)較短，呈馬蹄形渦旋或渦旋環(huán)。當(dāng)動(dòng)量厚度雷諾數(shù)>2 000時(shí)，渦旋結(jié)構(gòu)呈現(xiàn)細(xì)長(zhǎng)的發(fā)夾渦。結(jié)果表明，隨著雷諾數(shù)的增加，發(fā)夾渦破裂并呈現(xiàn)出不對(duì)稱的藤條結(jié)構(gòu)。

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

圖3顯示了不同的雷諾數(shù)對(duì)發(fā)夾狀渦旋拓?fù)涞挠绊慬13]。本文中，湍流邊界層入口動(dòng)量厚度雷諾數(shù)為400，屬于低雷諾數(shù)。從圖3可以看出，相干渦旋整體結(jié)構(gòu)基本較短，呈現(xiàn)出馬蹄形渦旋或渦旋環(huán)，與已有結(jié)論相符。

圖3 雷諾數(shù)對(duì)相干渦旋結(jié)構(gòu)的影響

圖4為湍流邊界層底部x+≈15處的流場(chǎng)云圖。從圖中可以清楚地看到沿翼展方向有高、低速交替帶。這種結(jié)構(gòu)被認(rèn)為與湍流邊界層的維持和發(fā)展密切相關(guān)。近年來(lái)，大量的數(shù)值研究證實(shí)了這種結(jié)構(gòu)不僅存在于不可壓縮流動(dòng)中，而且存在于可壓縮湍流邊界層中。

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

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

圖5 計(jì)算的結(jié)果與先前發(fā)表的結(jié)果之間的比較

可以看出，在黏性底層和對(duì)數(shù)區(qū)域，LES得到的曲線與經(jīng)典壁面律吻合較好。

線性律：

(4)

對(duì)數(shù)律：

(5)

以上公式中的參數(shù)由下面公式給出：

(6)

(7)

(8)

(9)

在經(jīng)典的對(duì)數(shù)律中，k為Von Karman常數(shù)[15]，其值為0.4～0.41，這里取0.41；C取5.25。動(dòng)力黏度μw是溫度的函數(shù)，可以通過(guò)Sutherland公式確定：

(10)

3 結(jié)語(yǔ)

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