改進k-ω-γ 轉捩模式對不同雷諾數下HIAD 的轉捩預測

劉宏康，陳堅強，向星皓，趙雅甜，

1．中南大學 交通運輸工程學院，長沙 410075

2．中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，綿陽 621000

3．中國空氣動力研究與發展中心 計算空氣動力研究所，綿陽 621000

隨著載人航天和行星探測任務對更大有效載荷需求的日益迫切，傳統減速器受運載火箭整流罩大小限制的弊端越來越明顯。高超聲速充氣式柔性減速器（Hypersonic Inflatable Aerodynamic Decelerator，HIAD）由于具有質量輕、可折疊、展開后阻力面積大等優勢成為了一種可行的替代方案，逐漸展示出了巨大的應用潛力。早在20 世紀末，美國國家航空航天局（NASA）便以未來大質量火星著陸探測任務為潛在工程背景，啟動了HIAD 項目，并于2009 年和2012 年先后發射充氣再入航天器IRVE-Ⅱ[1］和IRVE-3[2］以驗證該技術。飛行試驗表明，HIAD 柔性材料在氣動力作用下會變形為波紋狀（貝殼狀），從而促進流動轉捩為湍流，使得壁面熱流急劇增大，給熱防護系統的設計帶來挑戰。因此，準確預測變形后HIAD 外形的轉捩起始位置和壁面熱流至關重要。遺憾的是，之前的轉捩研究工作或針對光滑返回艙外形[3］，或基于孤立粗糙元，對這類波紋狀返回艙外形的轉捩研究十分有限。為此，2018年NASA 蘭利實驗室[4］開展了一系列風洞試驗來探究HIAD 壁面變形對邊界層轉捩和壁面熱流分布的影響規律。這些豐富而寶貴的試驗數據為轉捩數值方法的驗證和改進提供了基礎。

在當前所有針對邊界層轉捩的數值模擬方法中，基于雷諾平均N-S 方程（Reynolds Averaged Navier-Stokes Equations）的轉捩模式方法由于其經濟性和優越的魯棒性，被認為是最有可能應用于工程實踐的預測手段。目前已有數十種轉捩模式被提出，其中應用比較廣泛的有Menter 等[5-6］提 出 的γ-Reθ轉 捩 模 式 以 及 在 此 基礎上發展而來的C-γ-Reθ轉捩模式[7］，王亮和符松[8-10］提出的k-ω-γ 轉捩模式，Cho 和Chung[11］提出的k-ε-γ 轉捩模式和Walters 等[12］提出的kT-kLω 轉捩模式等。然而受限于轉捩現象本身的復雜性，目前發展的轉捩模式往往具有較為苛刻的適用條件，一般只適用于特定機制的轉捩現象，難以準確模擬復雜三維邊界層轉捩[13］。例如，基于低速試驗數據經驗擬合的轉捩模式往往在高速轉捩預測中效果不理想[14］。基于流向行波不穩定轉捩機理構造的模式，對橫流、分離等誘導轉捩的預測就會受到限制。而實際飛行器繞流往往存在多種轉捩機制相互干擾耦合[15］。陳堅強等[16］強調，目前國內外已經開始重視轉捩研究與真實飛行之間的匹配，所涉及的流動和轉捩機理也越來越復雜。因此，如何將多種轉捩失穩機制融入同一模式中是當前亟需解決的關鍵問題。

近年來，Zhou 等[17-18］拓展了k-ω-γ 轉捩模式的橫流轉捩預測功能，從而將其推廣應用至三維復雜外形的轉捩預測。Zhao[19］將考慮了橫流的k-ω-γ 轉捩模式應用于變形后HIAD 的轉捩預測，并 與γ-Reθ和kT-kL-ω 模 式 進 行 了 對 比。發現，kT-kL-ω 和γ-Reθ轉捩模式預測的結果均明顯偏離試驗值，只有k-ω-γ 模式得到了較為滿意的結果。然而遺憾的是，進一步研究發現k-ω-γ 模式對不同來流條件的適用性有待提高。當來流攻角或雷諾數發生改變，k-ω-γ 轉捩模式的預測性能開始惡化。這可能是其對一些轉捩機理考慮不充分導致的。考慮到HIAD 壁面變形后發生的流動分離現象，Zhao 等[20］進一步拓展了k-ωγ 模式的分離誘導轉捩預測功能。為便于區分，將其稱為改進k-ω-γ 轉捩模式，將未考慮分離誘導轉捩的稱為原始k-ω-γ 轉捩模式。

本文將改進后的k-ω-γ 模式應用于不同來流雷諾數下的HIAD 邊界層轉捩預測，從而：①通過與原始模式預測結果進行對比，評估和驗證改進k-ω-γ 轉捩模式對不同來流雷諾數下HIAD 轉捩現象的預測能力；②從模式構造的角度，揭示主導轉捩預測的不穩定性擾動及其相應作用范圍，為轉捩模式的發展和改進提供指導。

1 轉捩模式

1.1 原始k-ω-γ 轉捩模式

王 亮 和 符 松[8-10］受Langtry 等[5］當 地 變 量 構造方式的啟發，基于SST 湍流模式框架，同時借鑒Warren 等[21］在有效黏性系數中考慮非湍流脈動影響的思路，提出并發展了同時適用于亞聲速、超聲速和高超聲速轉捩預測的k-ω-γ 轉捩模式。Zhou 等[17-18］進一步拓展了k-ω-γ 轉捩模式的橫流轉捩預測功能。該模式包含湍動能k、湍動能比耗散率ω，以及間歇因子γ 這3 個輸運方程：

在具體形式上，k 方程和ω 方程與SST 湍流模式基本一致，不同之處在于采用有效黏性系數μeff替代了SST 模式中的湍流黏性系數μt。μeff的具體表達式為

式中：μnt為非湍流脈動黏性系數，

其中：τnt代表非湍流脈動時間尺度，它的模化考慮了第1 模態τnt1、第2 模態τnt2和橫流模態τcross的貢獻。根據穩定性理論結果，定義相對馬赫數Marel=（U-cr）/c 來區分第2 模態的作用區域：

式中：Eμ=0.5×（U-Uw）i2為相對壁面的平均流動動能；U（ys）代表廣義拐點ys處的速度；f（w）和f（Recf）分別為橫流速度和橫流雷諾數限制函數。更多模式構造細節可參見文獻[10，17］。

1.2 改進k-ω-γ 轉捩模式

k-ω-γ 轉捩模式不具備預測流動分離失穩的能力，這在一定程度上限制了其應用范圍。為此，文獻[20］拓展了k-ω-γ 轉捩模式的分離誘導轉捩預測性能。

該方法首先通過比較分子擴散時間尺度和瞬間壓力應變時間尺度之間的關系，構造了一個阻尼函數fss，以對有效長度尺度ζeff進行修正，從而避免分離泡前緣ζeff快速增長導致的有效渦黏性系數分布異常。fss的構造形式為

修正后的有效長度尺度ζeff_mod表示為

fss能很好地識別出剪切層高渦量區域。具體而言，在分離泡前緣附近強剪切層以及邊界層黏性底層很薄的區域fss=0，而在自由來流和邊界層大部分區域fss很快恢復到1。

其次，考慮到壓力梯度在分離誘導轉捩中起到重要作用，引入壓力梯度因子λζ作為指示器來直接反映壓力梯度對分離剪切層中轉捩的影響，其定義為

式中：V 是法向速度；y 為笛卡爾坐標系法向。該壓力梯度因子不僅嚴格基于當地變量和Galilean不變量，同時顯式的含有渦量Ω，以直接反映K-H不穩定中擾動在強剪切區域快速增長的特性。

最后，借由λζ構造分離間歇因子γsep：

式中：Freattach為控制函數，用以防止流動再附，表達式為

Fθt為邊界層保護函數：

最終的有效間歇因子γeff取γsep和輸運方程得到的間歇因子的最大值：

以上原始和改進的k-ω-γ 轉捩模式均已在自主開發的軟件平臺中實現，并通過高超聲速平板、尖錐、HIFiRE-5、X-51 前體等多類算例進行了驗證[22-25］。

2 物理模型及計算條件

2.1 物理模型

研究對象是NASA 蘭利實驗室2018 年風洞試驗[4］采用的模型，是IRVE-3 的0.050 8 比例縮比模型，縮比后直徑D=6 in（1 in=2.54 cm）。IRVE-3 是一種堆疊圓環型充氣減速器，主體由中心體和柔性充氣展開結構2 部分構成。中心體及其內部的復雜結構不在研究范圍內，因此將其簡化。充氣展開結構安裝在中心體的頭錐上，由6 個大的充氣式堆疊超環面和1 個位于肩部的小超環面構成，外部由F-TPS 鏈接，充氣展開后呈半錐角為60°的倒錐形。飛行過程中，在氣動載荷的作用下，F-TPS 受力變形嵌入到圓環體之間的間隙，進而形成波紋狀（貝殼狀）曲面。其變形后的外形輪廓如圖1 所示[19］。需用說明的是，為便于規律研究同時排除流固耦合效應的影響，蘭利的風洞試驗采用了一系列參數化變形的靜態模型來代替真實飛行條件下的壁面變形情況。本文選取其中一種變形程度（文獻[4］中記為Scallop-10）進行研究。

2.2 計算條件

參考蘭利實驗室針對HIAD 所做的風洞試驗結果[4］，結合對該轉捩流動已有認識[19］，選取大攻角下4 種不同來流雷諾數條件進行了數值模擬。具體來流參數設置見表1，表中：1 ft=0.304 8 m。

表1 來流條件Table 1 Freestream conditions

表1 中，Re、α、Ma∞、T∞和ρ∞分別為來流雷諾數、攻角、馬赫數、溫度和密度；Tw為壁面溫度；Tu∞為自由來湍流度，根據式（17）確定：

式中：γg為比熱比，等于1.4；PL代表風洞中的壓力脈動，根據文獻[26］通過式（18）計算：

計算采用半模多塊對接結構網格。為了能夠合理捕捉關鍵流動特征，對激波附近、壁面邊界層以及壁面曲率變化較大區域的網格進行了適當加密和優化布置，并保證邊界層內網格法向增長率不大于1.03，半模網格總量約為198 萬。針對HIAD 的轉捩預測，文獻[19，27］中開展了詳細的算法驗證和網格無關性驗證，此處不再贅述。

3 結果分析

3.1 典型流場結構

首先以來流雷諾數Re=2.16×106ft-1狀態為例，對HIAD 的典型流場結構進行分析。圖2為對稱面無量綱壓力云圖（P/P∞）與壁面熱流分布云圖h/hFR。由于頭部鈍度較大，HIAD 頭部前緣形成了一道明顯的脫體弓形激波，并且在大攻角（α=18°）來流條件下，迎風面激波較背風面更貼體。因此，波后迎風面壓強大于背風面，形成較大壓力差，驅使流動由迎風側流向背風側。其中變形后HIAD 壁面的“溝槽”結構迫使流動沿周向運動，從而導致了強橫流效應。

圖2 Re=2.16×106 ft-1時對稱面壓力與壁面熱流云圖Fig.2 Symmetry cut plane of pressure and surface of heat flux contour at Re=2.16×106 ft-1

圖3 給出了對稱面密度梯度（Density Gradient Magnitude， DGM）云圖和馬赫數等值線圖。在HIAD 的背風面由于凹凸起伏的外形存在復雜的波系結構，氣流流經波谷斜后方時被壓縮，經過波峰之后又膨脹，伴隨著氣流的減速與加速運動，形成了交替出現的壓縮波和膨脹波，這些激波/膨脹波串同時與頭部弓形激波發生相互作用，形成激波/激波干擾。并且由于HIAD的凹凸起伏構型，相鄰波峰之間的波谷處會形成類空腔流場，“空腔”內包含激波/邊界層干擾、剪切不穩定性、邊界層分離以及渦旋等復雜現象。激波和膨脹波間存在強的局部逆壓梯度。黏性和逆壓梯度的共同作用使得流動在波峰處發生分離，然后在波谷處再附，從而形成了明顯的分離線和再附線并在近壁區形成一個鏈狀低速區。

圖3 Re=2.16×106 ft-1時對稱面密度梯度及馬赫數云圖Fig.3 Symmetry cut plane of density gradient and Mach number contour at Re=2.16×106 ft-1

總體而言，壁面波紋變形使得流動結構變得更加復雜，最為顯著的便是橫流效應和流動分離現象的出現，二者都可能對邊界層轉捩產生顯著影響，從而對邊界層轉捩的精確預測帶來挑戰。

3.2 邊界層轉捩

圖4 給 出Re為 2.16×106、3.03×106、3.88×106和6.63×106ft-1時，分別基于原始和改進k-ω-γ轉捩模式預測得到的HIAD 外形壁面間歇因子（左）和壁面熱流分布（右）云圖，并給出風洞試驗的壁面熱流結果（記為EXP.）作為對比。間歇因子開始偏離0 或者壁面熱流顯著增大的位置判斷為轉捩的起始位置。其中，熱流h/hFR的計算公式可參見文獻[19］。

圖4 壁面間歇因子（左）和壁面熱流（右）云圖Fig.4 Surface intermittency （left） and heat flux （right） contour

由壁面熱流分布可以看出，曲率半徑較小的鼻端區域受熱嚴重。之后，伴隨著壁面的高低起伏，熱流在靠近波峰位置開始爬升，過波峰后又逐漸下降，從而形成了6 個連續上下的熱流波動。在波谷處由于流動分離的影響出現低壁面熱流區域。從間歇因子和熱流分布可觀察到一個由鼻端展開的扇形轉捩區，離圓心越遠的位置轉捩位置越靠前，并且熱流峰值往往出現在轉捩區。隨著雷諾數增大，轉捩位置逐漸前移。Re=2.16×106ft-1時，轉捩位置僅分布在背風面，而當Re為3.03×106～6.63×106ft-1時 迎 風 面 也開始出現了轉捩，特別是在靠近肩部的區域。同時注意到，隨著雷諾數的增大，過波峰后的低熱流區逐漸變小，這是由于隨著雷諾數的增大，轉捩位置提前，湍流抵抗逆壓梯度的能力更強，流動發生分離之后很快再附，分離泡變小。此外，雷諾數越大，轉捩區和湍流區的熱流隨之增大，高熱流的范圍逐漸擴大。

改進k-ω-γ轉捩模式準確地捕捉到了上述規律，預測得到的不同雷諾下的壁面熱流分布均與試驗值吻合較好，證明了該模式在這類波紋壁面轉捩模擬中的優越性。相比之下，原始k-ω-γ轉捩模式預測的結果與試驗存在偏差，主要體現為3 點：其一，轉捩起始位置延遲，特別是在靠近肩部的區域體現更加明顯；其二，轉捩發展過慢，轉捩區過長，使得該區域邊界層未充分發展，從而導致轉捩區熱流異常；其三，轉捩位置對雷諾數變化的響應較慢，這導致雷諾數越大，轉捩位置與實驗結果相差越大。

因此，接下來基于改進k-ω-γ模式預測結果，對不同雷諾數下的HIAD 數值模擬結果進行定量對比。圖5 首先給出了不同雷諾數下壁面中心子午線熱流分布。在迎風子午線及鼻端位置，流動為層流，不同雷諾數下的中心子午線熱流相差不大。在背風子午線，流動繞過鼻端后立刻轉捩為全湍流。因此，隨著雷諾數增大，背風子午線波峰和波谷處的熱流均不同程度的增大。由此可知，在中心子午線上，來流雷諾數對背風面熱流的影響顯著大于迎風面。

圖5 壁面中心子午線熱流Fig.5 Surface centerline heat flux

圖6 進一步比較了中心子午線的邊界層厚度變化，圖中同時給出了壁面變形高度ksc=21.87 mile（1 mile=1 609.344 m）作為參考。需要說明的是，由于無量綱總焓比不受激波和邊界層內拐點的影響，在工程外形應用中具有較高的魯棒性，因此對邊界層外緣的識別采用總焓比準則，即定義邊界層內當地總焓比達到自由來流總焓0.995 處的位置為邊界層外緣。

圖6 中心子午線邊界層厚度Fig.6 Surface boundary layer heights

可以發現，伴隨著壁面起伏，邊界層厚度在波谷處增厚，在波峰處變薄，且迎風面的邊界層厚度基本低于壁面變形高度ksc。之后隨著流動由迎風面流向背風面，邊界層沿流向逐漸發展并發生轉捩，其厚度急劇增加，并遠大于ksc。對比發現，轉捩后背風面邊界層厚度可達迎風面厚度的2～6 倍。其次，對比不同雷諾數的邊界層厚度可知，隨著來流雷諾數的增大，邊界層厚度逐漸減小，這與理論結果一致。

3.3 轉捩預測機制

由于HIAD 繞流的高度復雜性和強非線性，其邊界層轉捩是多種不穩定擾動相互耦合的結果。為探究改進k-ω-γ轉捩模式對HIAD 邊界層的預測機制，明確模式構造中各擾動的主導區域。圖7 給出了改進模式預測的流場中各不穩定擾動的影響區域。其中，γsep為所改進模式所構造的分離間歇因子，μ∞代表自由來流分子動力黏性系數，μnt1、μnt2和μcross分別為第1 模態、第2 模態和橫流模態對應的非湍流脈動黏性系數，具體定義如下：

由于波峰處和波谷處的流動特征存在顯著區別，圖7 同時提供了位于波峰（x=0.64 in）和波谷（x=0.77 in）的橫截面作為對比。顯然，流動特征不同，波峰和波谷處轉捩模式預測的主導模態也存在區別。需要特殊說明的是，結果顯示無論在波峰處還是在波谷處，不同雷諾數下μnt2在的值均為0。因此，圖中未給出μnt2的分布。這與Chang 和Hollis 等[28］得出的結論一致，他們通過試驗和穩定性分析研究了MSL 返回艙的邊界層轉捩現象，發現對于這類大鈍頭體返回艙，其前方形成的強弓形激波使過激波后邊界層外緣相對馬赫數低于第2 模態起作用的臨界值。因此，對于本文研究的HIAD 外形來說，盡管來流馬赫數高達6，但轉捩依然不受第2 模態不穩定的影響。

圖7 各不穩定模態黏性系數及分離間歇因子在波峰x=0.64 in（左）和波谷x=0.77 in（右）截面近壁區分布Fig.7 Different mode disturbance viscosity and separation intermittency distribution on crest x=0.64 in （left） and in valley x=0.57 in （right）

在波峰位置，HIAD 轉捩預測主要由改進模式構造的分離誘導間歇因子主導，第1 模態對轉捩起始位置的影響不大，起到加速邊界層失穩的效果，而橫流模態不發揮作用。對比不同雷諾數下的μnt1可以發現，在波峰位置處，第1 模態的作用區域主要集中在對稱面兩側附近，隨著雷諾數的增大，μnt1開始增大的位置雖然提前，但是其對雷諾數的變化并不敏感，作用區域始終出現在對稱面兩側附近。相比之下，γsep的作用區域更廣，其從轉捩的起始位置一直到對稱面附近，且γsep開始增大的位置對雷諾數的變化非常敏感。雷諾數越大，γsep增長起始點越提前，與前文顯示的轉捩起始位置相一致。證明了本文構造的γsep能有效反映轉捩位置隨來流雷諾數變化的規律。

在波谷位置，預測結果顯示邊界層轉捩主要由第1 模態、橫流模態以及流動分離失穩引起。波谷與波峰處的顯著差異在于第1 模態和橫流模態，由于波谷處強橫流效應的存在，其邊界層轉捩受到橫流模態的影響。隨著雷諾數增大，μcross量值顯著增大。橫流模態的作用區域由轉捩起始位置附近延伸到對稱面附近，并且波谷位置處的第1 模態作用范圍比波峰處的作用范圍更大。產生這一變化的原因可以從兩方面進行分析，從物理機理角度來看，波谷處第1 模態作用范圍的增大可能是由于橫流渦的影響。由于橫流渦的“上拋下掃”導致更多的不穩定性進入到邊界層內，進而激發邊界層內第1 模態的增大，導致其第1 模態影響區域變大；從模式構造的角度來看，橫流效應誘導邊界層失穩之后，流動發生轉捩，邊界層增厚，使得模式中的有效長度尺度ζeff增大，從而使得第1 模態對應的非湍流脈動黏性系數的增長，進而導致了第1 模態的影響區域變大。對比不同雷諾數下的μnt1、μcross和γsep分布可發現，隨著雷諾數的增大，3 者增長起始點均顯著提前，且γsep對雷諾數的變化尤為敏感，進一步揭示了引入分離誘導轉捩預測功能的重要性。

由上可知，HIAD 復雜流動失穩是多重不穩定擾動模態共同作用的結果，對其邊界層轉捩的準確預測離不開對相應轉捩機理的合理模化。

4 結 論

將拓展了分離誘導轉捩預測功能的改進k-ωγ 轉捩模式，應用于高超聲速充氣式柔性減速器（Hypersonic Inflatable Aerodynamic Decelerator，HIAD）的轉捩預測研究。評估了不同來流雷諾數下，改進k-ω-γ 模式對HIAD 邊界層轉捩的預測性能，詳細分析了其預測機制，并與原始k-ω-γ模式進行了對比。所得主要結論如下：

1）壁面變形使得HIAD 流場結構變得更加復雜，其背風側存在交替分布的膨脹波和壓縮波串，伴隨的逆壓梯度形成了相應的流動分離和再附。同時大攻角條件導致的周向壓力差會誘導橫向流動。橫流和流動分離，均會對邊界層轉捩產生重要影響，結合可能存在的第1 模態和第2模態失穩，多重不穩定擾動的共存對轉捩模式預測性能提出了挑戰。

2）改進后的k-ω-γ 轉捩模式同時具備對第1模態、第2 模態、橫流模態以及分離失穩的預測能力。該模式可準確預測HIAD 的扇形轉捩區陣面。且不同來流雷諾數下，其預測的轉捩起始點，轉捩區形狀以及壁面熱流等均與風洞試驗結果吻合良好，體現出其對該類波紋壁面復雜外形轉捩預測的優越性。相比之下，原始k-ω-γ 模式存在轉捩區發展過慢，轉捩起始點隨雷諾數變化不明顯等缺陷。

3）波紋壁面波峰與波谷處流動結構不同，使得轉捩預測機制存在顯著差異。在波峰位置，轉捩預測主要由改進模式構造的流動分離失穩主導。而在波谷位置，由第1 模態、橫流模態以及流動分離共同誘導。無論波峰或波谷處，由于過強弓形激波后，邊界層外緣相對馬赫數較小，第2模態均不發揮作用。隨雷諾數增大，轉捩預測主導模態的起始增長點前移，其中分離間歇因子對雷諾數尤為敏感，最終導致高雷諾數下轉捩位置提前。因此，對于多重不穩定擾動模態耦合誘導下的邊界層轉捩預測，需充分考慮相應流動中可能存在的失穩機制，并對其進行合理模化。