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

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

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

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

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

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

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

1 計算方法

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

1.1 γ-Reθ轉捩模型

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

(1)

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

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

(2)

其中：Fonset為決定轉捩開始的閾函數；Flength為轉捩區域長度；其余相關參數的含義詳見文獻[9]。

(3)

式中：Reν為當地渦雷諾數；Reθc為臨界動量厚度雷諾數。

(4)

(5)

(6)

(7)

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

(8)

此外，γ-Reθ模型還考慮了分離誘導轉捩的作用。

(9)

γeff=max(γ,γsep)

(10)

式中：γsep為分離區的間歇因子；γeff為最終的間歇因子；Freattach為與再附區有關的參數，具體形式參見文獻[9]。其可通過改變經驗參數s1來調節對分離誘導轉捩的敏感度。最后，轉捩模型對SST模型k方程的作用方式為

(11)

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

1.2 γ-Reθ轉捩模型與IDDES方法的結合

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

(12)

(13)

長度尺度函數從DES到IDDES發展歷程的具體情況可參見文獻[33]，本文僅給出lIDDES的計算形式。lIDDES與網格尺度Δ有關。

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

(14)

式中：hmax=max(Δx,Δy,Δz)，hmin=min(Δx,Δy,Δz),Δx、Δy、Δz分別為x、y、z這3個方向上的網格尺度；dw為距壁面距離；Cw為與亞格子模型無關的經驗常數，一般取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框架后，會對轉捩模型產生2個方面的影響：

2 IDDES方法的驗證

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

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

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

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

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

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

圖2 二維超聲速混合層流場結構Fig.2 Flow field structures of two-dimensional supersonic mixing layer

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

圖3 不同流向位置的時間統計平均流向速度分布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θ轉捩模型的建立及驗證

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

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

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

基于以上經驗，本文擬對基于IDDES方法的γ-Reθ轉捩模型進行如下修改：

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

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

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

3.1 馬赫數為3.5的尖錐轉捩

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

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

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

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

圖6 不同Re下的尖錐轉捩Fig.6 Transition of sharp cone with different Re

3.2 馬赫數為8.1的雙鍥轉捩

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

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

圖7 雙楔試驗構型示意圖[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 平板上單個粗糙顆粒誘導的轉捩

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

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

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

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

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

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

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

圖14 粗糙顆粒附近流場結構Fig.14 Flow field structure around roughness element

圖15 對稱面流場結構Fig.15 Flow field structures on symmetry plane

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

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

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

4 結 論

1)將IDDES方法與γ-Reθ轉捩模型結合起來并發揮各自的優勢，從而建立形成了新的基于IDDES框架的轉捩模型。

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

3)計算結果表明新建立的方法在強制轉捩模擬上要優于傳統的γ-Reθ轉捩模型。