合成湍流對空腔流動RANS-LES混合模擬結果的影響

郭啟龍, 李 辰, 劉朋欣, 孫 東, 張涵信

(1. 中國空氣動力研究與發展中心 空氣動力學國家重點實驗室, 綿陽 621000;2. 中國空氣動力研究與發展中心 計算空氣動力研究所, 綿陽 621000;3. 北京航空航天大學 國家計算流體力學實驗室, 北京 100191)

0 引 言

高速的空腔流動會產生強烈的振動與噪聲，是航空航天領域中一類受到廣泛關注的流動問題[1-3]。數值模擬是研究空腔流動機理與振動/噪聲控制方法的重要手段。從流動分離的角度來看，空腔流動是一類典型的具有固定分離點(空腔的前緣點)的大范圍分離流動，其分離區的特征尺度與空腔尺寸相當。在較高雷諾數條件下，空腔流動中既包含由外形引起的大尺度三維旋渦結構，又包含隨機無序的小尺度的湍流渦旋，這使得在高雷諾數條件下很難獲得高速空腔流動的高保真數值模擬結果。

直接數值模擬(DNS)和大渦模擬(LES)[4-6]能夠較為精確地分辨空腔流動中包含大范圍分離的多尺度湍流，但需要消耗很大的計算資源(尤其是在高雷諾數情況下)，目前很難大范圍地用于高雷諾數空腔流動的預測。雷諾平均Navier-Stokes(RANS)方法通過半經驗的模型來?；械耐牧髅}動，只對最大尺度的平均流進行求解，其計算量相對較小，但現有的RANS方法在含有大范圍分離的流動中難以給出準確的預測結果。

近年來RANS-LES混合方法(Hybrid RANS-LES method, HRLM)在大范圍分離流動的模擬中獲得了廣泛的關注[7-9]。這種方法一般是僅在遠離壁面的大尺度分離區分辨大渦并采用LES進行模擬，而在其他區域(如附體的湍流邊界層內)采用RANS來模擬，從而對于高雷諾數壁面湍流的模擬兼具了RANS的高效率和LES的湍流模擬能力。根據混合方式可分為分區混合方法和非分區混合方法[10]。

非分區混合方法通過建立依賴網格尺度和流動參數的經驗關系式來決定RANS和LES的區域，其流場中沒有明確的RANS和LES的區分界面，如Speziale[11]構造的模型。而在分區混合方法中，需要基于對流動機理的認識人為確定交界面，并在兩側分別使用RANS和LES，如Embedded LES[12]。

無論是否人為分區，RANS解和LES解之間的匹配區域在HRLM中是普遍存在的，即所謂的“灰區”題。一般情況下，HRLM計算的可靠性由LES決定，因此從RANS解向LES解過度導致的“灰區”也是HRLM中關注最廣泛的問題，在這個過程中分辨的湍流脈動與模化的湍流應力是不匹配的，雖然混合方法進入了LES的分支，但是數值解并不能正確的給出可分辨的“大渦”與相應的亞格子應力，這就需要有一段適應區域以供LES逐步建立起真實的湍流脈動成分。

分區或非分區HRLM均能適用于高雷諾數空腔流動的模擬，例如文獻[13,14]中的Embedded LES和文獻[15]中的脫體渦模擬(DES)。由于來流邊界層的分離點及分離區域的位置是固定的，因此HRLM 的“灰區”通常位于分離點附近(分區模擬時人為確定的分界面位于接近空腔前緣的上游位置)。

在分區模擬的交界面上添加人工合成的湍流脈動能夠有效地減小下游適應區的長度以使LES盡快完成真實亞格子應力的恢復。在使用HRLM對附體的湍流邊界層流動或中/小范圍分離湍流進行模擬時，合成湍流能夠有效減小“灰區”的不利影響[16]，但是對于類似空腔流動的大范圍分離湍流，關于合成湍流對HRLM模擬的影響并沒有明確的結論。直觀上看，來流邊界層中添加合理的脈動有助于更真實地描述分離后的初始剪切層中固有的Kelvin-Helmholtz不穩定性，但是Sagaut等[17]認為流動中絕對主導的自持振蕩機制使得是否精細地描述來流邊界層中的真實湍流結構對模擬結果的影響十分有限。在文獻[13]中的Embedded LES中，交界面設置在空腔前緣位置導致了初始剪切層仍然以二維的旋渦結構為主。Lawson等[18]通過添加合成渦(SEM)[19-20]給出了來流脈動對Rossiter模態的影響，但其計算采用的展向對稱邊界條件可能掩蓋剪切層中的展向失穩模態。

本文對馬赫數0.6的空腔流動開展了分區和非分區RANS-LES混合模擬，具體計算采用基于S-A模型的改進的延遲DES(IDDES)方法[21]，IDDES作為DES的一種改進，具備同時模擬附體邊界層和大范圍分離流動的能力，其在高速高雷諾數空腔流動的模擬中已經得到的廣泛的應用[22-23]。當來流邊界層中含有合適的湍流脈動時，IDDES模擬等同于Wall-Modeled LES(壁面模型由RANS部分提供)，因此IDDES也能夠用于大范圍分離湍流的分區模擬(類似Embedded LES)[16]。分區模擬的交界面上添加了人工合成的湍流脈動，主要目的即在于給出有/無合成湍流(Synthetic Turbulence，ST)時所得到空腔流動的時均特性和動態特性，通過詳細地對比分析明確來流邊界層中ST對空腔流動HRLM模擬結果的影響。

1 計算設置

1.1 空腔外形與計算條件

本文模擬采用的空腔外形參考文獻[2]中的風洞實驗，空腔的長度(L)、深度(D)和寬度(W)之間的比值為L∶D∶W=6∶1∶2，在空腔上游設置有一段長度等于空腔長度(L)的平板(如圖1)。來流馬赫數Ma=0.6，基于來流參數的雷諾數ReL=ρ∞U∞L/μ∞=2.46×106。

圖1 空腔流動計算設置示意圖Fig.1 Schematic of the case-setup of the cavity flow

計算采用多塊對接網格(如圖2所示)，在腔內包括369×181×185個點(分別對應x、y、z方向)，空腔外部包括585×181×301個點，整個網格共計含有4.42×107個網格點，網格分布在各個壁面及空腔開口處的剪切層附近加密，壁面第一層網格對應壁面單位約為1.0。

圖2 計算網格示意圖Fig.2 Schematic of the computational grids

1.2 方法與模型

本文針對可壓縮Navier-Stokes方程進行求解，其在三維笛卡爾坐標系(x,y,z)下的無量綱守恒形式為：

(1)

式中，Q=[ρ,ρu,ρv,ρw,e]為守恒變量，ρ為密度、p為壓力、u=(u,v,w)為速度矢量，總能e可表示為：

(2)

E、F、G表示空間各個方向上的無黏通量，Ev、Fv、Gv表示黏性通量，各通量的具體表達式參見文獻[24]。

基于S-A模型的IDDES輸運方程為：

(3)

(4)

(5)

壁面邊界采用無滑移絕熱邊界條件；對于除固壁外的其他各個開邊界(包括來流邊界、出口邊界、遠場邊界)，為了抑制聲波虛假反射對模擬結果的“污染”，本文采用文獻[26]的特征邊界處理，并在開邊界與流動區域之間額外添加了一層邊界緩沖區，緩沖區內的控制方程添加了人工阻尼項，從而確保旋渦等大尺度流動結構順利流出邊界而不會造成反射[27-28]。

為了對比人工合成湍流對流動模擬結果的影響，本文設置了兩個算例，在第一個算例中使用文獻[29]提出的ST構造方法(記為“Case 0: ST-on”)，即在選定的交界面上添加湍流速度脈動u′，其計算方式為：

(6)

(7)

(8)

其中，上標“n”表示第n個模態的參數，q是由當地的von Karman模型能譜近似給出的模態幅值，kn=k·dn為三維的波數矢量，k表示波數矢量的模，d為均勻分布于單位半徑的球面上的隨機方向矢量，σ為隨機生成的、與d垂直的單位矢量，φ是均勻分布于[0,2π]區間的隨機相位，r′為表征擾動傳播的位置矢量。該ST構造方法具有效率高、魯棒性強的優點，更具體的細節參見文獻[29]。此時數值模擬相當于分區HRLM。ST添加的流向位置在空腔上游距離前緣0.5D處(如圖1中紅色虛線所示)，根據此處的邊界層厚度估算，0.5D的流向長度能夠保證足夠的適應區距離。在第二個算例中不添加任何人工脈動(記為“Case 1: ST-off”)，此時全場使用IDDES進行計算，相當于非分區HRLM。

圖3中給出了底面中線(z=0)上的平均壓力系數Cp和整體聲壓級OASPL分布?？梢钥闯鰞煞N情況下計算得到的底面中線上的Cp在前半部分與實驗結果符合較好，而在x/L>0.6的范圍內兩組計算的結果均略高于實驗值，其中Case 0與實驗值更接近，這將在后面結合流場結構進一步分析；兩個算例的OASPL與實驗結果具有較好的一致性。整體來說，數值模擬與實驗結果的對比證明了本文計算的可靠性。

(a) 壓力系數

(b) 整體聲壓級

2 結果與討論

本節從流動的時均特性和動態特性兩方面分析添加ST對HRLM計算結果的影響。

2.1 流動時均特性對比分析

通過對非定常流動進行時間平均，可以對流場結構的主要特征進行分析。圖4中給出了展向中截面上的流線分布(背景云圖為壓力系數)，可以看出中截面上的流動由兩個環流區構成，其中較大的環流區靠近前壁面，其流向尺度占據了空腔長度的一半以上，較小的環流區靠近于后壁。當有ST時，后壁環流的尺寸明顯大于無ST的結果，且具有更強的渦量，較強的渦量會導致更低的壓力分布，這也解釋了圖3(a)中的時均壓力系數分布在后壁附近出現的差別。

(a) Case 0: ST-on

(b) Case 1: ST-off

Zhang等[30]對相近工況下(Ma=0.6，L/D=6)的空腔流動開展了實驗和LES模擬研究，其雷諾數分別為ReL,exp=1.875×106以及ReL,LES=6×104，與本文有所不同，但是一般認為雷諾數對空腔流動中主控的自持振蕩的影響較小。圖5中給出了Zhang等[30]得到的展向截面上時均流線分布，其中計算和實驗均體現了與本文結果相似的主要環流結構，本文的計算結果相比于文獻中的LES模擬更接近實驗，而未添加合成湍流的后壁環流區略微偏小。

(a) 實驗測量結果

(b) LES模擬結果

剪切層在空腔流致振蕩中扮演了重要角色，基于平均場可以得到剪切層的渦量厚度δ和動量厚度θ，用來考察剪切層的定量特性。二者的定義分別為:

(9)

(10)

式中y1、y2、U1、U2的定義見文獻[31, 32]。圖6中給出了δ和θ沿流向的變化，可以看出有無ST對剪切層平均厚度的影響很小，渦量厚度曲線大致可分為兩段線性增長的區域，這與Beresh等[32]和Larcheveque等[31]的結果類似，兩個區域的線性增長率分別為dδ/dx≈0.18和dδ/dx≈0.09，略低于Beresh等[32]得到的結果。圖6(a)中，本文結果與Zhang等的計算和實驗數據[30]取得了較好的一致性。動量厚度近似以dθ/dx≈0.043的斜率線性增長。添加ST使Case 0得到的動量厚度在4

為了估計在剪切層內三維渦結構的平均對流速度，我們沿著剪切層中心線(y=0，z=0)，計算了時空自相關函數，其定義為:

(11)

(a) 渦量厚度

(b) 動量厚度

其中〈*〉表示時間平均，v表示法向脈動速度(能較好的反映對流效應)，自相關函數具有較大值的位置代表流動結構在時空平面內的運動軌跡，該軌跡的斜率drx/drt近似代表了剪切層內流動結構的平均對流速度。圖6中給出了兩個算例中自相關函數在rx-rt平面內的分布，可以估算出Case 0中對流速度約為drx/drt=0.539，而Case 1中對流速度約為drx/drt=0.571，基本符合Rossiter公式中表征對流速度的常數取值(0.55～0.57)。二者之間的差異說明上游添加的ST使剪切層渦結構的平均對流速度略有降低，在剪切層厚度變化(圖5)非常接近的前提下，ST對對流速度的影響可能是通過改變上游的流動三維機制引入的，這一點將在后文進一步闡述。

圖7 沿剪切層中心的時空自相關函數Fig.7 Spatial-temporal auto-correlation along the center of the shear layer

2.2 流動動態特性對比分析

在圖8和圖9中分別給出了由Q準則表示的瞬時旋渦結構和展向中間截面(z=0)上的瞬時渦量分布。在圖8(a)中可以清楚的看到與上游邊界層中所添加的ST對應的三維旋渦結構。ST對結果的影響主要集中在剪切層中。在Case 1中，來流邊界層僅由RANS給出了平均剖面(圖8(b))，進入剪切層區域后，在初始階段迅速形成展向近似均勻的渦結構，經過一定的發展距離流動的三維特性逐漸增強，并開始出現明顯的流向渦結構。從圖9中也能看出，兩個算例在剪切層初始階段的流動結構尺度上存在差別，經過大約2D的流向距離，這種定性上的差別逐漸變得不明顯。

圖8 Q準則顯示的瞬時旋渦結構Fig.8 Instantaneous vortical structures displayed by the Q-criterion

(a) Case 0: ST-on

(b) Case 1: ST-off

(a) 流向平均速度

(b) 湍動能

圖11中給出了展向中截面的雷諾剪切應力τxy的分布，并與Zhang等[30]的計算和實驗結果進行了對比。在剪切層靠近前緣的位置，Case 0添加合成湍流的計算結果更符合實驗，而Case 1以及文獻[30]中的LES結果在此處分辨的τxy剪切應力的絕對值偏小，體現了來自“灰區”的耗損，這也與湍動能的比較體現了類似的差異(見圖10(b))。在剪切層的后部，本文計算結果與實驗測量的結果取得了很好的一致。

(a) Case 0: ST-on

(b) Case 1: ST-off

(c) Zhang等[30]的實驗測量結果

(d) Zhang等[30]的LES模擬結果

Zhang等所開展的LES模擬在入口添加了隨機Fourier模態擾動，但是從模擬結果可以看出剪切層初始階段的流動結構仍然以大尺度展向渦為主(見文獻[30]中圖9)，說明邊界層/剪切層中的小尺度湍流結構并未得到充分的激發，入口附近仍存在“灰區”。可以認為文獻[30]中的LES模擬結果與實驗測量結果之間的差異也部分地來源于剪切層初始階段的小尺度湍流結構。

圖12中給出了空腔底面中線上x/D=1.32、3.0位置上的壓力脈動的功率譜密度，通過對比可以看出添加ST對Rossiter模態(圖12中較低Strouhal數的峰值)的影響很小。在x/D=1.32位置的高頻區域(圖12(a))，Case 0的結果相比于Case 1出現了額外的高頻峰值(Strouhal數約為3.645處)，這也對應著來流邊界層中的ST對腔內流動的誘導作用；而在更下游的x/D=3.0位置，高頻區域Case 0中的峰值消失。

結合2.1節中關于平均特性的對比分析結果，可以發現ST對模擬結果的直接影響主要表現為增強剪切層初始階段的脈動量，而隨著流動的發展，有無ST對中下游位置的剪切層的影響非常微小，且圖12中的對比也說明二維平面內的自持振蕩機制并沒有因添加ST而改變，說明兩個算例在后壁附近的環流區域出現的差異并不是由二維機制導致的。

為了定量考察流動結構的三維特性，可定義剪切層湍流脈動的展向波數譜為:

(12)

(a) x/D=1.32

(b) x/D=3.0

其中，w(x,y,z,t)表示非定常展向速度，ξ為展向波數，積分限取為l=0.9D，以排除側壁對展向均勻性的影響。

圖13中給出了x-y平面內六個不同位置上時間平均的〈Sz〉。在剪切層中心(y/D=0)，Case 0在最上游的x/D=0.5位置上(圖13(a))得到的波數譜上明顯出現了在ξ≈31.6處的峰值，而Case 1的結果中沒有，說明該波數的產生與上游添加的ST有關(在圖中記為“ST-induced spanwise modes”)，而隨著流動向下游發展，在x/D=3.0和x/D=5.2位置上，兩個算例得到的展向波數譜變得非常接近。

在接近底面的不同流向位置上(圖13(b、d、f))，Case 0的展向波數譜也表現出了上述展向模態，說明ST改變了對腔內流動三維機制的模擬結果，而這種變化并未出現在對下游位置剪切層的模擬結果中。綜合以上分析可以認為有無ST所導致的HRLM模擬結果的差異主要來自于三維流動的區別，而主導的自持振蕩(本質上為二維機制)并未受到明顯影響。根據與實驗數據的對比，可以認為添加ST使空腔流動HRLM模擬結果更符合實際流動。

圖13 不同位置上的時間平均展向波數譜Fig.13 Time-averaged spanwise wavenumber spectra at different locations

3 結 論

本文對馬赫數0.6的空腔流動開展了分區及非分區的RANS-LES混合模擬，對比了在分離點上游添加合成湍流對剪切層以及腔內的分離區模擬產生的影響：

1) 在剪切層初始階段，合成湍流顯著增加了湍動能的大小，旋渦結構具有一定的展向波數，而未添加合成湍流時，剪切層初始階段的結果主要呈現出展向近似均勻的渦結構。

2) 對于腔內的流動，添加合成湍流增大了后壁附近環流區的尺寸，進而降低了靠近后壁區域的剪切層動量厚度和平均壓力分布。在添加合成湍流情況下，底面靠近后壁區域的壓力系數更加接近實驗數據。

3) 合成湍流的添加未對主導流動的自持振蕩機制產生明顯的影響，而腔內流動模擬結果差別主要來源于合成湍流改變了腔內流動的三維機制。