基于IDDES框架的γ -Reθ轉(zhuǎn)捩模型

易淼榮，趙慧勇，樂嘉陵，肖保國，鄭忠華

中國空氣動力研究與發(fā)展中心 高超聲速沖壓發(fā)動機技術(shù)重點實驗室，綿陽 621000

邊界層轉(zhuǎn)捩對高超聲速飛行器的設(shè)計有利有弊。一方面，邊界層轉(zhuǎn)捩導(dǎo)致氣動熱的急劇增加，嚴(yán)重影響熱防護系統(tǒng)的設(shè)計，同時轉(zhuǎn)捩還影響飛行器的摩阻、壓阻、飛行器的性能和氣動力控制；另一方面，在高超聲速進氣道上實現(xiàn)邊界層強制轉(zhuǎn)捩，可以減少流動分離，提高流動抗逆壓的能力，便于進氣道的起動和超燃發(fā)動機的正常工作。在長達一個多世紀(jì)的研究后，人們對于邊界層轉(zhuǎn)捩問題的研究取得了長足的進步，發(fā)展了線性穩(wěn)定性、弱非線性穩(wěn)定性、瞬態(tài)增長、三波共振、二次失穩(wěn)等理論[1]描述邊界層轉(zhuǎn)捩的過程，發(fā)展了靜音風(fēng)洞來模擬天上真實的來流擾動條件[2]。但是，針對一個給定的較復(fù)雜構(gòu)型，人們?nèi)圆荒艽_切地給出邊界層會在什么條件下、什么位置由層流轉(zhuǎn)捩為湍流，特別是在高超聲速條件下，由于受到太多因素及各種因素之間互相耦合作用的影響，轉(zhuǎn)捩預(yù)測變得更加困難。

目前能夠進行邊界層轉(zhuǎn)捩研究的數(shù)值方法大致可分為4種：第1種是通過精細的空間和時間分辨率來解析邊界層發(fā)展過程中擾動的物理發(fā)展過程，從而捕獲邊界層由層流轉(zhuǎn)捩為湍流的詳細信息，主要為DNS(Direct Numerical Simulation)和LES(Large Eddy Simulation)方法，被用來研究典型狀態(tài)下的轉(zhuǎn)捩機制和詳細的轉(zhuǎn)捩過程[3-4]，但是由于該方法對網(wǎng)格、格式和計算量的要求極大，很難針對高雷諾數(shù)、復(fù)雜構(gòu)型開展研究。第2種是通過線性和非線性穩(wěn)定性理論，通過對擾動增長率的研究從而獲得邊界層內(nèi)擾動波的增長直至失穩(wěn)過程，如LST(Linear Stability Theory)和PSE(Parabolized Stability Equation)方法[5-6]，其有明確的物理機制，能夠模擬非穩(wěn)定波線性(甚至是弱非線性階段)的發(fā)展，但無法模擬Breakdown和湍斑生成的過程，也無法模擬Bypass轉(zhuǎn)捩。第3種是通過對大量試驗結(jié)果進行總結(jié)分析，提取出關(guān)鍵的流場參數(shù)，給定一個標(biāo)準(zhǔn)，當(dāng)參數(shù)達到該標(biāo)準(zhǔn)時，則認為轉(zhuǎn)捩發(fā)生，稱之為轉(zhuǎn)捩準(zhǔn)則，如NASP(National AeroSpace Plane)的尖前緣平板轉(zhuǎn)捩準(zhǔn)則：Reθ/Mae=305[7]，其中Reθ為以動量厚度為長度尺度定義的雷諾數(shù)，Mae為邊界層外緣馬赫數(shù)。該方法形式簡單，操作方便，但是通用性不好，只能針對類似的構(gòu)型在類似的條件下開展轉(zhuǎn)捩預(yù)測，而針對其他情況的預(yù)測就會出現(xiàn)偏差。

由于粗糙顆粒明顯改變了邊界層內(nèi)流場的結(jié)構(gòu)，其誘導(dǎo)的尾跡渦具有G?rtler渦的性質(zhì)，這種反向旋轉(zhuǎn)的流向渦對導(dǎo)致邊界層內(nèi)的速度剖面出現(xiàn)大的變形，使得流向速度在法向和展向均存在拐點。G?rtler渦結(jié)構(gòu)在經(jīng)歷瞬態(tài)增長后，可能會與第2模態(tài)之間發(fā)生非線性耦合作用[21]或因其條帶結(jié)構(gòu)本身的不穩(wěn)定性[22]而二次失穩(wěn)并最終導(dǎo)致轉(zhuǎn)捩的發(fā)生。如果采用足夠精細的數(shù)值模擬方法，通過對渦的形成和破碎過程進行捕捉，理論上其可以捕獲流向渦的形成、發(fā)展、失穩(wěn)到轉(zhuǎn)捩為湍流的整個過程。要捕獲這些過程，目前可以采用的方法有DNS、LES和RANS/LES混合方法，而目前最有可能大范圍工程應(yīng)用的是RANS/LES混合方法。Yoon等[23]首先將DES方法應(yīng)用于高超聲速邊界層內(nèi)單個粗糙顆粒誘導(dǎo)轉(zhuǎn)捩的算例，結(jié)果表明DES模型可以分辨大渦結(jié)構(gòu)，可以用來模擬粗糙顆粒誘導(dǎo)的強制轉(zhuǎn)捩；肖志祥等采用IDDES方法分別研究了高超聲速下粗糙顆粒誘導(dǎo)平板強制轉(zhuǎn)捩[24]和凹槽誘導(dǎo)轉(zhuǎn)捩[25]的情況，給出了轉(zhuǎn)捩渦結(jié)構(gòu)的發(fā)展過程，表明采用IDDES方法對強制轉(zhuǎn)捩開展研究是可行的，但其在進行模擬時，采用的是強制粗糙顆粒或者凹槽上游的流動為全層流的情況，因此在流過粗糙顆粒或凹槽之后，有一個從層流到湍流的轉(zhuǎn)捩過程。IDDES方法之所以能夠模擬強制轉(zhuǎn)捩，是因為其能夠解析更小尺度的渦，能夠捕獲粗糙顆粒或凹槽誘導(dǎo)出來的渦結(jié)構(gòu)，但由于IDDES方法仍然是一個針對湍流的計算方法，當(dāng)沒有間歇因子來控制其湍動能方程時，針對沒有明顯大渦結(jié)構(gòu)、由小擾動經(jīng)過線性增長和非線性增長之后破碎再轉(zhuǎn)捩為湍流的過程，其卻無法捕獲。而且純粹的IDDES方法也無法反映來流湍流度對轉(zhuǎn)捩的影響。根據(jù)普渡大學(xué)對于X-51A的試驗結(jié)果[26]，來流湍流度對粗糙顆粒誘導(dǎo)的轉(zhuǎn)捩仍然有重要的影響，因此來流湍流度也是一個不可忽略的因素。

針對γ-Reθ轉(zhuǎn)捩模型和IDDES方法各自的優(yōu)點和缺點，如果將γ-Reθ轉(zhuǎn)捩模型與IDDES方法結(jié)合起來，則理論上可以克服各自的缺點，從而同時滿足對自然轉(zhuǎn)捩和強制轉(zhuǎn)捩的模擬。實際上，早在2011年Sorensen等[27]就嘗試過將γ-Reθ轉(zhuǎn)捩模型與DES方法結(jié)合起來，并進行了大迎角機翼的計算，隨后Alam等[28]則將DDES方法與k-kL-ω轉(zhuǎn)捩模型[29]進行了結(jié)合，并將其應(yīng)用在低壓渦輪機葉片的轉(zhuǎn)捩預(yù)測中。喬磊等[30]也嘗試過將DDES方法與γ-Reθ轉(zhuǎn)捩模型相結(jié)合，對低速平板算例進行了計算。以上結(jié)果均表明轉(zhuǎn)捩模型與DES類方法結(jié)合后，能夠像原始的γ-Reθ轉(zhuǎn)捩模型一樣對自然轉(zhuǎn)捩進行模擬。但以上均為試探性研究，研究的構(gòu)型都較為簡單，都是低速領(lǐng)域的轉(zhuǎn)捩，沒有進一步發(fā)揮基于IDDES框架的轉(zhuǎn)捩模型相對于傳統(tǒng)轉(zhuǎn)捩模型能夠捕捉更精細的流場結(jié)構(gòu)從而更精確地模擬強制轉(zhuǎn)捩的這一優(yōu)勢。

本文希望通過結(jié)合γ-Reθ轉(zhuǎn)捩模型對自然轉(zhuǎn)捩模擬的優(yōu)點和IDDES方法在模擬強制轉(zhuǎn)捩時能分辨更精細物理結(jié)構(gòu)從而捕獲轉(zhuǎn)捩過程的優(yōu)點，將γ-Reθ轉(zhuǎn)捩模型添加到IDDES方法的框架中去，從而實現(xiàn)對自然轉(zhuǎn)捩和強制轉(zhuǎn)捩更精確的模擬。在沒有粗糙顆粒、凹槽等強制轉(zhuǎn)捩因素時，網(wǎng)格只需滿足RANS的網(wǎng)格要求即可，此時在邊界層內(nèi)RANS方法會被激活，并對轉(zhuǎn)捩進行模擬。在有粗糙顆粒、凹槽等強制轉(zhuǎn)捩因素時，采用更精細的網(wǎng)格，邊界層內(nèi)大部分區(qū)域LES方法被激活，可對渦結(jié)構(gòu)的出現(xiàn)及發(fā)展過程進行捕捉，從而實現(xiàn)對強制轉(zhuǎn)捩的模擬。

1 計算方法

本文搭建的基于IDDES框架的的γ-Reθ轉(zhuǎn)捩模型是在課題組開發(fā)的大規(guī)模并行CFD平臺AHL3D[31]上添加模塊建立的，AHL3D 采用基于網(wǎng)格平均的有限體積方法離散求解Navier-Stokes方程組，定常時間推進采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法，非定常計算則包括顯示的龍格-庫塔法和隱式的雙時間步法。無黏通量采用重構(gòu)-推進方法，重構(gòu)采用3階MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)方法，本文在此基礎(chǔ)上又添加了5階WENO(Weighted Essentially Non-Oscillatory)重構(gòu)方法。無黏通量的推進采用Steger-Warming通量分裂、AUSM+(improved Advection Upstream Splitting Method)、AUSMPW+(AUSM+ by Pressure-based Weighted function)和LDFSS(Low Diffusion Flux-Splitting Scheme)4種計算格式。黏性通量采用改進的Gauss定理計算。湍流模型中包含了本文需要用到的剪切應(yīng)力輸運(SST)雙方程模型。作者團隊在文獻[20]中已經(jīng)實現(xiàn)了γ-Reθ轉(zhuǎn)捩模型的搭建和應(yīng)用，因此本文以該轉(zhuǎn)捩模型為基礎(chǔ)，先實現(xiàn)IDDES方法，然后將IDDES方法與γ-Reθ轉(zhuǎn)捩模型相結(jié)合，重新對轉(zhuǎn)捩模型進行標(biāo)定，從而實現(xiàn)本文計算方法的建立。

1.1 γ-Reθ轉(zhuǎn)捩模型

γ-Reθ模型詳細的計算公式可參考文獻[9]。其通過間歇因子γ對SST模型湍動能k方程的產(chǎn)生項和破壞項進行控制，進而控制整個轉(zhuǎn)捩過程。間歇因子輸運方程為

(1)

式中：ρ為密度；Uj為xj方向的雷諾平均速度；μ為分子黏性系數(shù)；μt為湍流黏性系數(shù)；σf=1.0為常數(shù)；Pγ和Eγ分別為γ方程的產(chǎn)生項和破壞項，Pγ的形式為

Pγ=Flengthca1ρS(γFonset)0.5(1-ce1γ)

(2)

其中：Fonset為決定轉(zhuǎn)捩開始的閾函數(shù)；Flength為轉(zhuǎn)捩區(qū)域長度；其余相關(guān)參數(shù)的含義詳見文獻[9]。

(3)

式中：Reν為當(dāng)?shù)販u雷諾數(shù)；Reθc為臨界動量厚度雷諾數(shù)。

(4)

(5)

(6)

(7)

Reθt=f3(Tu)f4(λθ)

(8)

此外，γ-Reθ模型還考慮了分離誘導(dǎo)轉(zhuǎn)捩的作用。

(9)

γeff=max(γ,γsep)

(10)

式中：γsep為分離區(qū)的間歇因子；γeff為最終的間歇因子；Freattach為與再附區(qū)有關(guān)的參數(shù)，具體形式參見文獻[9]。其可通過改變經(jīng)驗參數(shù)s1來調(diào)節(jié)對分離誘導(dǎo)轉(zhuǎn)捩的敏感度。最后，轉(zhuǎn)捩模型對SST模型k方程的作用方式為

(11)

式中：Pk和Dk分別為SST模型中k方程的產(chǎn)生項和破壞項。

1.2 γ-Reθ轉(zhuǎn)捩模型與IDDES方法的結(jié)合

γ-Reθ模型通過添加兩個輸運方程求解間歇因子γ，且通過γ作用于k方程的產(chǎn)生項和破壞項從而與SST湍流模型結(jié)合。在形式上，添加基于SST模型的IDDES方法則更為簡單，可以通過修改湍動能輸運方程破壞項中的長度尺度函數(shù)而實現(xiàn)[32]。SST模型中湍動能方程為

(12)

(13)

長度尺度函數(shù)從DES到IDDES發(fā)展歷程的具體情況可參見文獻[33]，本文僅給出lIDDES的計算形式。lIDDES與網(wǎng)格尺度Δ有關(guān)。

Δ=min{max[Cwdw,Cwhmax,hmin],hmax}

(14)

式中：hmax=max(Δx,Δy,Δz)，hmin=min(Δx,Δy,Δz),Δx、Δy、Δz分別為x、y、z這3個方向上的網(wǎng)格尺度；dw為距壁面距離；Cw為與亞格子模型無關(guān)的經(jīng)驗常數(shù)，一般取Cw=0.15。

lIDDES的形式為

(15)

(16)

fdt=1-tanh[(Cdtrdt)3]

(17)

(18)

fB=min{2exp(-9α2),1.0}

(19)

其中:α=0.25-dw/hmax。至此，lIDDES求出。

將RANS框架換成IDDES框架后，會對轉(zhuǎn)捩模型產(chǎn)生2個方面的影響：

2 IDDES方法的驗證

超聲速混合層流動是典型的以大尺度渦結(jié)構(gòu)為主導(dǎo)的可壓縮自由剪切層湍流流動。本節(jié)通過Goebel的混合層試驗[34]來檢驗本文的IDDES方法對大尺度渦結(jié)構(gòu)的捕捉能力以及對平均量和脈動量的預(yù)測能力。

Goebel的試驗構(gòu)型[34]如圖1所示，混合區(qū)總高48 mm，總長500 mm，隔板下部傾角為2.5 °，尖劈厚度為0.5 mm。混合區(qū)入口參數(shù)如表1所示。計算中為滿足混合區(qū)入口處的附面層厚度，參考文獻[33]中的做法，將隔板上部入口選在混合區(qū)入口前67 mm處，隔板下部入口選在混合區(qū)入口前42 mm處。

計算求解二維非定常Navier-Stokes方程組，采用基于SST模型的IDDES方法進行模擬。無黏通量重構(gòu)采用五階WENO方法，分裂采用AUSMPW+格式，黏性通量采用改進的高斯定理，時間推進采用隱式雙時間步法，物理時間步長為5×10-8s，子迭代最大步數(shù)為15，子迭代CFL數(shù)為0.8。上下壁面采用對稱邊界條件，隔板采用絕熱壁條件，出口采用外推。混合區(qū)網(wǎng)格量為1 001×261，垂直壁面第1層網(wǎng)格法向長度取0.01 mm。隔板上下的初始條件分別設(shè)為上下入口的來流。計算經(jīng)過1 ms后流場充分發(fā)展并穩(wěn)定，再計算1 ms,統(tǒng)計1～2 ms時的平均量。

計算得到的t=1 ms時的流場結(jié)構(gòu)如圖2所示。從圖中可以看出，本文IDDES方法能夠分辨混合層中大尺度渦結(jié)構(gòu)的形成、發(fā)展和脫落過程。t=1～2 ms期間的3個不同流向位置的時間統(tǒng)計平均流向速度U無量綱化之后的分布如圖3所示，圖中h=48 mm為混合區(qū)總高度，U2=399 m/s為混合層下層入口速度，ΔU=301 m/s為上下層入口速度之差。x=50 mm時，計算得到的混合層比試驗稍薄，這可能是由于來流條件中沒有給定脈動值，造成剪切層發(fā)展比試驗值要稍慢導(dǎo)致的。而到x=100 mm和150 mm時，計算與試驗吻合得非常好，此時旋渦開始變大并脫落，混合層厚度明顯增加。

圖1 Goebel的混合層試驗示意圖[34]Fig.1 Sketch of Goebel’s mixing layer test[34]

表1 混合區(qū)入口來流條件Table 1 Free stream conditions of mixing region inlet

圖2 二維超聲速混合層流場結(jié)構(gòu)Fig.2 Flow field structures of two-dimensional supersonic mixing layer

圖4和圖5則分別為流向脈動速度均方根urms和法向脈動速度均方根vrms的分布，流向脈動速度均方根在x=50 mm時仍然比試驗值偏小，這與平均速度顯示混合層偏薄相一致。而在x=100 mm 和150 mm位置的流向速度均方根則與試驗吻合良好，表明本文計算方法能夠有效捕捉到流場大渦結(jié)構(gòu)。而法向速度均方根則在x=50 mm時與試驗吻合較好，在x=100 mm和150 mm時較試驗值要明顯偏大，計算得到的流向、法向速度脈動在同一位置基本是同等量級的，但試驗測量到的流向脈動要大于法向脈動，計算與試驗的差別可能是由于計算是采用二維模擬而導(dǎo)致的。總體來說，計算結(jié)果與試驗結(jié)果的對比表明本文計算方法足以分辨超聲速混合層中的大渦結(jié)構(gòu)，可應(yīng)用于可壓縮湍流的模擬。

圖3 不同流向位置的時間統(tǒng)計平均流向速度分布Fig.3 Distribution of time-average streamwise velocity in different streamwise positions

圖4 流向脈動速度均方根分布Fig.4 Root-mean-square distribution of streamwise velocity fluctuations

圖5 法向脈動速度均方根分布Fig.5 Root-mean-square distribution of vertical velocity fluctuations

3 基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型的建立及驗證

γ-Reθ轉(zhuǎn)捩模型可以反映雷諾數(shù)對轉(zhuǎn)捩的影響，也可以通過式(8)反映壓力梯度、湍流度對轉(zhuǎn)捩的影響，但由于其對轉(zhuǎn)捩位置的預(yù)測受到計算方法、格式等多方面因素的影響，因此每一個框架上搭建的γ-Reθ轉(zhuǎn)捩模型都可能需要進行重新標(biāo)定。作者團隊基于在文獻[19-20]中采用γ-Reθ模型進行計算的經(jīng)驗，得到以下幾點結(jié)論：

1)在AHL3D框架上搭建的基于RANS的γ-Reθ轉(zhuǎn)捩模型如果直接采用原始的經(jīng)驗關(guān)系式進行高超聲速轉(zhuǎn)捩模擬，得到的轉(zhuǎn)捩位置整體會比實際轉(zhuǎn)捩位置大大提前。

2)通過重新標(biāo)定將轉(zhuǎn)捩位置整體推遲后，添加Forsythe[35]的可壓縮性修正可有效反映可壓縮性對轉(zhuǎn)捩的影響，使得該方法能在較寬馬赫數(shù)范圍內(nèi)對轉(zhuǎn)捩進行模擬。

基于以上經(jīng)驗，本文擬對基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型進行如下修改：

1)加入Forsythe的可壓縮性修正。

2)經(jīng)驗耦合關(guān)系式仍然采用Langtry的原始關(guān)系式，參考文獻[15]中的做法，將式(3)中的2.193調(diào)整為3.5。

3)式(9)中，常數(shù)2會致使γ在分離區(qū)直接變?yōu)?，導(dǎo)致k值的增加，從而致使本應(yīng)由LES方法處理的區(qū)域被程序自動判斷為RANS處理區(qū)域，這樣會導(dǎo)致分離區(qū)的精細流動因采用RANS方法而被抹去，實際上，當(dāng)將RANS方法改為IDDES方法后，對分離流的捕捉要更加精細，因此，為了適當(dāng)弱化模型對分離誘導(dǎo)轉(zhuǎn)捩模化的功能，當(dāng)需要用IDDES方法來分辨流場內(nèi)的精細流動時，將式(9)中的常數(shù)2改為1。

3.1 馬赫數(shù)為3.5的尖錐轉(zhuǎn)捩

Chen和Malik[36]在NASA蘭利研究中心的馬赫數(shù)Ma=3.5低擾動風(fēng)洞中進行了半錐角為5°的尖錐轉(zhuǎn)捩測量，其來流單位雷諾數(shù)Re由2.81 ×107/m變化到7.8×107/m，得到了相同馬赫數(shù)條件下雷諾數(shù)對轉(zhuǎn)捩位置的定量影響規(guī)律。采用該算例不僅可以驗證模型對可壓縮流轉(zhuǎn)捩的模擬能力，也可以驗證模型是否能正確反應(yīng)雷諾數(shù)對轉(zhuǎn)捩的影響規(guī)律。

本文選取了其試驗中的6個不同來流雷諾數(shù)進行計算，并與試驗得到的轉(zhuǎn)捩位置進行對比。計算設(shè)置的來流條件如表2所示，表中T0和P0分別為來流總溫和總壓。尖錐迎角為0°，壁面溫度統(tǒng)一設(shè)為絕熱壁，來流湍流度設(shè)為0.4%，來流黏性比RT=1.0。計算采用軸對稱模型進行，方法為基于IDDES框架的γ-Reθ轉(zhuǎn)捩模型，無黏通量采用WENO五階重構(gòu)和AUSMPW+格式。模型全長為381 mm。尖錐流向和法向的網(wǎng)格量為341×201，壁面第1層網(wǎng)格距離壁面0.001 mm。

由于試驗是通過測量溫度恢復(fù)系數(shù)來判斷流動是層流還是湍流，而真實的溫度壁面條件在計算中并不好給定，因此本文采用摩阻系數(shù)來判斷轉(zhuǎn)捩位置，且只與試驗對比轉(zhuǎn)捩位置，計算結(jié)果如圖6所示，其中Cf為摩阻系數(shù)，r為溫度恢復(fù)系數(shù)。從圖中可知，隨著來流雷諾數(shù)的增加，轉(zhuǎn)捩位置由約0.28 m逐漸前移至0.1 m左右。模型可以反映可壓縮條件下雷諾數(shù)變化對轉(zhuǎn)捩的影響。

表2 馬赫數(shù)為3.5的尖錐轉(zhuǎn)捩來流條件Table 2 Free stream condition of sharp cone transition tests with Mach number of 3.5

圖6 不同Re下的尖錐轉(zhuǎn)捩Fig.6 Transition of sharp cone with different Re

3.2 馬赫數(shù)為8.1的雙鍥轉(zhuǎn)捩

Neuenhahn和Olivier[37]在德國TH2高超聲速激波風(fēng)洞中進行了來流馬赫數(shù)為8.1的雙鍥中邊界層與激波干擾的試驗。試驗研究了壁溫和熵層對激波邊界層干擾的影響。在研究者后續(xù)針對該問題的CFD研究中發(fā)現(xiàn)，該試驗中包含了邊界層轉(zhuǎn)捩現(xiàn)象，因為單純的層流或湍流均無法正確得到壓力分布和分離區(qū)大小。其中Krause[16]和Zhang[17]等均采用該算例對可壓縮條件下的γ-Reθ轉(zhuǎn)捩模型進行了驗證。本文采用二維計算，模擬了尖前緣、前緣半經(jīng)R為0.5 mm和1.0 mm這3個不同模型在同一來流條件下的試驗情況。構(gòu)型如圖7所示，第1楔長度L1=180 mm，傾斜角為9°，第2楔長度L2=255 mm，拐角為11.5°。來流馬赫數(shù)為8.1，總壓為5.7 MPa，總溫為1 430 K，單位雷諾數(shù)為3.8×106/m。壁溫設(shè)為300 K，來流湍流度設(shè)為0.5%，來流黏性比RT設(shè)為1。計算方法同3.1節(jié)，網(wǎng)格總量約為3.4萬，第1、2楔面的流向網(wǎng)格數(shù)分別為121和101，法向網(wǎng)格數(shù)為131，壁面第1層網(wǎng)格距離壁面小于0.002 mm。

計算得到的壁面壓力(p)分布如圖8所示，從圖中可知，層流與湍流計算出來的拐角位置的分離泡大小有很大區(qū)別，而試驗值均介于層流和湍流計算結(jié)果之間，確實存在明顯的邊界層轉(zhuǎn)捩現(xiàn)象。同時，分離泡的大小隨R由0 mm增加到1.0 mm而逐漸增大，表明轉(zhuǎn)捩位置隨著前緣鈍度的增大而逐漸后移。計算結(jié)果反應(yīng)出了邊界層轉(zhuǎn)捩給分離泡大小帶來的影響，也反應(yīng)出了前緣鈍度(0～1.0 mm范圍內(nèi))對轉(zhuǎn)捩的影響規(guī)律。

圖7 雙楔試驗構(gòu)型示意圖[37]Fig.7 Sketch of double wedge geometry[37]

圖8 不同前緣半徑下雙楔的壁面壓力分布Fig.8 Wall pressure distribution of double wedge with different leading edge radiuses

3.3 平板上單個粗糙顆粒誘導(dǎo)的轉(zhuǎn)捩

對基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型進行修改并經(jīng)過自然轉(zhuǎn)捩算例驗證后，對單個粗糙顆粒誘導(dǎo)高超聲速邊界層轉(zhuǎn)捩的情況進行模擬，分析轉(zhuǎn)捩過程。本文選取了Tirtey等在高超聲速VKI H3風(fēng)洞中進行的強制轉(zhuǎn)捩試驗[38-39]進行模擬。Tirtey等研究了圓柱形、鉆石型、斜坡型的單個粗糙顆粒誘導(dǎo)轉(zhuǎn)捩的情況，本文選取其圓柱形粗糙顆粒在高雷諾數(shù)下的情況進行模擬。模型如圖9所示[38]，一塊長 290 mm、寬100 mm、前緣鈍度為0.5 mm的平板，在距前緣60 mm處安裝了一個粗糙顆粒，其參數(shù)為：a=4 mm，b=1.02 mm。本文只模擬了其長170 mm、寬40 mm的區(qū)域，流場區(qū)域高20 mm，如圖10所示。來流馬赫數(shù)為6，來流單位雷諾數(shù)為2.7×107/m，來流總壓為3.1 MPa，總溫為500 K，壁溫在試驗開始時均為292±2 K，試驗結(jié)束后層流區(qū)最大溫升為4 K，湍流區(qū)最大溫升為12 K，由于這里壁溫變化都不大，因此選取300 K為壁面溫度。整個流場的網(wǎng)格量約為3 500萬，平板第1層網(wǎng)格距離壁面0.002 mm，對應(yīng)的y+小于0.1。粗糙顆粒下游的流場流向、法向、展向的網(wǎng)格分別為601×151×207，對粗糙顆粒周圍的流場進行了加密。粗糙顆粒后緣的流向和展向網(wǎng)格間距分別為0.005 mm和0.03 mm，對應(yīng)無量綱網(wǎng)格尺度Δx+<0.1，Δz+<1，尾跡區(qū)的網(wǎng)格間距隨著流向距離的增大而逐漸增加，到x=170 mm位置時，Δx=0.6 mm，Δz=0.2 mm，對應(yīng)的Δx+≈18，Δz+≈6。

計算得到的以Q等值面表示的流場渦結(jié)構(gòu)如圖11所示(圖中T為溫度)，可以看出粗糙顆粒誘導(dǎo)出的馬蹄渦和剪切層除了向下游發(fā)展外，還有明顯的往展向拓寬趨勢，這符合湍流的性質(zhì)，也與試驗結(jié)果相符。圖12給出了壁面Stanton數(shù)(St)分布。可知，采用基于RANS的γ-Reθ轉(zhuǎn)捩模型得到的粗糙顆粒尾跡在x=70 mm(粗糙顆粒后緣下游10 mm處)位置時寬度約為12 mm，在x=170 mm位置的寬度約為14 mm，試驗測量得到的尾跡寬度在x=70 mm時約為9 mm，在x=170 mm時約為22 mm，而采用基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型得到的尾跡寬度在x=70 mm時約為9 mm，在x=170 mm時約為20 mm。顯然本文模擬結(jié)果與試驗情況更為接近。

圖9 平板模型示意圖[38]Fig.9 Sketch of panel model geometry[38]

圖10 流場及網(wǎng)格拓撲示意圖Fig.10 Sketch of topology of flow field and grid

圖11 以溫度染色的Q等值面圖Fig.11 Iso-surfaces of Q colored by temperature

圖12 壁面St分布Fig.12 St distribution on wall

圖13為本文計算得到的St在尾跡區(qū)(對稱面上)及非尾跡區(qū)的分布與試驗對比，可知，計算與試驗基本吻合，St在粗糙顆粒后緣均先劇烈增大(由粗糙顆粒后緣誘導(dǎo)出的尾跡渦導(dǎo)致)并迅速衰減，再經(jīng)歷一個相對緩慢的增長過程，達到極大值后開始隨著流向距離的增加逐漸減小。若定義St達到全湍流值之前的極大值點為轉(zhuǎn)捩結(jié)束位置，則計算得到的轉(zhuǎn)捩結(jié)束位置為x≈120 mm，對應(yīng)的轉(zhuǎn)捩雷諾數(shù)約為3.24×106，而試驗得到的轉(zhuǎn)捩結(jié)束位置為x≈108 mm，對應(yīng)的轉(zhuǎn)捩雷諾數(shù)約為2.92×106。二者在轉(zhuǎn)捩位置和尾跡展向發(fā)展的差別的原因可能是本文的計算方法雖然一方面捕捉到了粗糙顆粒尾渦的形成、發(fā)展及失穩(wěn)過程，另一方面也模化了自然轉(zhuǎn)捩時擾動發(fā)展導(dǎo)致轉(zhuǎn)捩的過程，由于來流擾動是以湍動能k值的形式給定的，而非直接給定具體形式的脈動，因此無法反映來流擾動與粗糙顆粒之間的耦合作用，而在該試驗中，風(fēng)洞為常規(guī)噪聲風(fēng)洞，因此邊界層內(nèi)粗糙顆粒之前的擾動可能已經(jīng)足夠大并與粗糙顆粒產(chǎn)生了耦合作用，而目前本文的方法并不能反映出這種耦合作用。

圖13 尾跡區(qū)與非尾跡區(qū)的St分布Fig.13 St distribution of wake region and out of wake region.

圖14為計算得到的粗糙顆粒附近對稱面流線和壁面摩擦力線圖與Tirtey給出的流場結(jié)構(gòu)對比。從圖中可以看出，粗糙顆粒前緣誘導(dǎo)出了微弱的斜激波，而在顆粒前緣上游則存在較大的分離區(qū)，較大的分離渦又誘導(dǎo)出二次反向分離渦，致使粗糙顆粒周圍的摩擦力線存在3條分離線和3條再附著線，分離渦向下游發(fā)展，形成了馬蹄渦結(jié)構(gòu)，數(shù)值模擬對粗糙顆粒附近的流場結(jié)構(gòu)的模擬與試驗基本一致。

圖15給出了對稱面上的流場結(jié)構(gòu)，分別為x方向的密度梯度、流向渦大小及速度脈動的均方根分布，其中u∞表示來流速度。從圖中可知，在對稱面上，粗糙顆粒上部除誘導(dǎo)出了斜激波外，還誘導(dǎo)出了較強的剪切層，該剪切層厚度很薄，擾動量級很大，但基本維持在邊界層外緣，實際上，在粗糙顆粒下游，距離壁面較近的底層也出現(xiàn)了剪切層，底層部分的剪切層在向下游發(fā)展的同時也在向上層發(fā)展，同時厚度逐漸增厚，渦結(jié)構(gòu)最終破碎，擾動擴散至整個邊界層，擾動在x=120 mm位置附近達到一個極大值之后，擾動幅值不再繼續(xù)增加，此時轉(zhuǎn)捩已經(jīng)完成。

圖14 粗糙顆粒附近流場結(jié)構(gòu)Fig.14 Flow field structure around roughness element

圖15 對稱面流場結(jié)構(gòu)Fig.15 Flow field structures on symmetry plane

圖16、圖17則分別為瞬態(tài)流場中不同法向位置的溫度分布和不同流向位置流向渦的渦量分布(圖中Δx′為不同流向截面與粗糙顆粒中心的距離)。從圖中可知，粗糙單元后存在兩層剪切層(與對稱面上的結(jié)果相一致)，上部剪切層很薄，在向下游的傳播過程中逐漸擴散到外流并最終因耗散作用而徹底消失。而下部分的剪切層則與粗糙單元兩側(cè)的馬蹄渦一開始分別各自向下游發(fā)展并緩慢擴散，直到三簇渦在展向上越來越靠近并最終開始相互耦合，之后較大的渦結(jié)構(gòu)開始破碎成較小的渦結(jié)構(gòu)，流動變得不穩(wěn)定，開始變成湍流。在x=120 mm處，即Δx′=60b時，渦結(jié)構(gòu)已經(jīng)完全破碎，渦量在展向分布趨于均勻，表明此時轉(zhuǎn)捩已經(jīng)完成，與前面壁面熱流系數(shù)和對稱面速度脈動均方根分布體現(xiàn)的結(jié)果一致。

圖16 壁面不同法向位置的瞬態(tài)溫度分布Fig.16 Instantaneous temperature distribution of different vertical distances from wall

圖17 不同流向位置的流向渦渦量云圖Fig.17 Streamwise vorticity magnitude contours of different streamwise positions

4 結(jié) 論

1)將IDDES方法與γ-Reθ轉(zhuǎn)捩模型結(jié)合起來并發(fā)揮各自的優(yōu)勢，從而建立形成了新的基于IDDES框架的轉(zhuǎn)捩模型。

2)通過超聲速混合層、高超聲速自然轉(zhuǎn)捩及粗糙顆粒誘導(dǎo)的強制轉(zhuǎn)捩等一系列算例的計算，表明新建立的方法在模擬自然轉(zhuǎn)捩時，能夠像γ-Reθ模型一樣通過基于轉(zhuǎn)捩動量厚度的經(jīng)驗關(guān)系式得到轉(zhuǎn)捩位置，而在模擬粗糙顆粒誘導(dǎo)的強制轉(zhuǎn)捩時，又能發(fā)揮IDDES方法能夠捕捉到粗糙顆粒在邊界層中誘導(dǎo)出的擾動的特點，能夠?qū)Υ植陬w粒誘導(dǎo)出的剪切層和馬蹄渦結(jié)構(gòu)的形成、耦合和失穩(wěn)過程進行較為精細地刻畫和模擬，從而達到模擬強制轉(zhuǎn)捩的目的。

3)計算結(jié)果表明新建立的方法在強制轉(zhuǎn)捩模擬上要優(yōu)于傳統(tǒng)的γ-Reθ轉(zhuǎn)捩模型。