C- -Re 高超聲速三維邊界層轉捩預測模型

向星皓，張毅鋒，袁先旭，涂國華，萬兵兵，陳堅強,2,*

1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000

2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

高超聲速三維邊界層橫流轉捩的準確預測是高超聲速飛行器氣動設計面臨的主要困難之一。由于飛行器三維曲面構型和飛行姿態(tài)角(攻角、側滑角)的影響，邊界層內往往存在與無黏流線方向垂直的流動分量，稱為橫流流動[1]。在真實飛行中，由橫流速度拐點導致的橫流失穩(wěn)機制在三維邊界層中普遍存在，并成為三維邊界層轉捩的主導因素[2]。

基于雷諾平均Navier-Stokes(RANS)方程的轉捩模型是解決轉捩實際工程問題的重要手段。在高超聲速橫流轉捩模型方面，王亮[3]提出了k-ω-γ新型湍流/轉捩模型，對高超聲速圓錐預測效果良好；周玲等[4]基于傳統(tǒng)橫流雷諾數(shù)判據(jù)對k-ω-γ轉捩模型進行了改進；張毅鋒等[5]改進的γ-Reθ t轉捩模型實現(xiàn)了對零攻角尖錐、橢錐表面轉捩陣面模擬；徐晶磊等[6]提出了KDO(Kinetic Dependent Only)單方程湍流/轉捩一體化模型，實現(xiàn)了對橢球體、尖錐表面的轉捩預測。此外，國內外進行了大量的低速橫流轉捩模型與判據(jù)研究[2,7]，取得了較好的預測效果，但在高超橫流轉捩方面的模型研究較少。

目前高超聲速橫流轉捩領域由于缺乏足夠豐富的變參數(shù)試驗結果，例如表面粗糙度、噪聲等，高超聲速橫流轉捩判據(jù)與預測模型的構建與應用遇到較大困難。本文采用線性穩(wěn)定性分析eN方法對高超聲速橫流轉捩數(shù)據(jù)進行拓展，結合橫流強度與表面粗糙度構造當?shù)鼗母叱曀贆M流轉捩判據(jù)，面向高超聲速三維邊界層對基于Chant 2.0計算平臺的高超聲速修正γ-Reθ轉捩模型進行完全當?shù)鼗臋M流拓展，建立適用于高超聲速三維邊界層轉捩預測的C-γ-Reθ橫流轉捩預測模型。

1 計算方法

1.1 控制方程

采用有限體積法求解RANS方程。無黏通量采用AUSM(Advection Upstream Splitting Method)格式、二階精度MUSCL(Monotonic Upwind Scheme for Conservation Law)格式離散，利用Van-Leer和Minmod混合限制器捕捉激波間斷。黏性通量采用二階中心差分格式離散，時間推進采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)隱式方法，采用MPI(Message Passing Interface)技術進行并行計算。

1.2 基礎轉捩模型

湍流模型采用兩方程k-ωSST(Shear Stress Transfer)模型?；A轉捩模型采用基于Chant 2.0高超計算平臺的高超聲速修正γ-Reθ轉捩模型[5]，能夠預測無攻角情況下的高超聲速平板、直錐、裙錐等的流向轉捩。

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式中：Mae為邊界層外緣馬赫數(shù)；C1～C4為修正系數(shù)[5]；y為壁面距離；v為黏性系數(shù)。

高超聲速修正包括壓力梯度參數(shù)修正與湍流普朗特數(shù)修正[5]，使之能夠預測無攻角情況下的高超聲速尖錐流向轉捩，如圖1[5]所示，圖中x為軸向位置，對比了無修正、壓力梯度參數(shù)λθ修正以及兩種(λθ+Prt)修正疊加γ-Reθ模型預測的尖錐壁面溫度回復因子r分布。

圖1 尖錐壁面溫度恢復因子r分布[5]Fig.1 Temperature recovery factor r distribution on sharp cone surface[5]

2 構建橫流轉捩模型

能體現(xiàn)流動物理機制的高保真模型是提高CFD計算結果精準度的重要基礎[7]。從現(xiàn)有模型和判據(jù)的構造出發(fā)，說明γ-Reθ轉捩模型無法預測橫流轉捩的原因，并闡釋基于傳統(tǒng)橫流雷諾數(shù)判據(jù)的橫流模型在高超聲速橫流轉捩預測中失效的機制。在此基礎上通過線性穩(wěn)定性eN方法對高超聲速橫流轉捩數(shù)據(jù)進行拓展，利用橫流強度與壁面粗糙度構造當?shù)鼗母叱曀贆M流轉捩判據(jù)，對γ-Reθ轉捩模型進行高超橫流拓展，初步建立適用于高超聲速三維邊界層的C-γ-Reθ轉捩預測模型。

2.1 γ-Reθ轉捩模型的橫流預測不適用性

動量厚度雷諾數(shù)輸運方程的生成項定義為

(8)

式中：cθ t為生成項系數(shù)，cθ t=0.03；Reθ t為經(jīng)驗關系擬合雷諾數(shù)；Fθ t為邊界層識別函數(shù)，使源項在邊界層外開啟，在邊界層內關閉。因此在邊界層外，輸運量趨于當?shù)財M合值Reθ t；在邊界層內，輸運量主要依靠擴散作用。該邊界層識別函數(shù)也被用于橫流源項構造。

圖2 橢球表面摩擦力系數(shù)分布[10]Fig.2 Skin friction coefficient distributions on prolate spheroid[10]

2.2 低速橫流判據(jù)/模型的高超聲速不適用性

傳統(tǒng)橫流雷諾數(shù)判據(jù)一般將橫流雷諾數(shù)定義為Recrossflow=Wmaxδ10%/v，其中Wmax為沿壁面法向的最大橫流速度，δ10%為橫流速度為最大橫流速度10%的位置(壁面較遠處)對應的壁面法向距離。

天津大學趙磊[9]發(fā)現(xiàn)，同一流場中采用邊界層名義厚度δ99無量綱化的橫流速度型幾乎完全重合，如圖3[9]所示，其中us和uc分別為按其最大值歸一化的勢流方向速度和橫流方向速度。這意味著不同站位下的δ10%/δ99基本上是一定的，令其比值等于Cδ，則有

(9)

圖3 超聲速后掠鈍頭平板不同站位歸一化勢流、橫流速度分布[9]Fig.3 Normalized potential flow and crossflow profiles at different locations for hypersonic swept blunt plate[9]

式中:δ99和v與來流相關。

通過超聲速、高超聲速尖錐數(shù)值計算發(fā)現(xiàn)，最大橫流速度Wmax與來流馬赫數(shù)緊密相關，并且變化幅度較大。圖4為固定雷諾數(shù)、來流馬赫數(shù)Ma=3～7情況下，尖錐側面90°子午線最大橫流速度沿軸向的分布。對于橫流速度的馬赫數(shù)相關性，趙磊在鈍頭后掠平板流動中也有類似發(fā)現(xiàn)[9]。由圖4可知，尖錐表面橫流強度大部分情況下超過12%，該強度也超過了Reed和Haynes[11]提出的橫流準則適用范圍,即最大橫流速度Ucf,max與邊界層外緣速度Ue比值不超過8%(Ucf,max/Ue≤8%)。固定來流雷諾數(shù)前提下，最大橫流速度具有較強的馬赫數(shù)相關性，導致采用單一臨界值的橫流雷諾數(shù)判據(jù)在圖4所示高超聲速情況下不再適用。

圖4 不同馬赫數(shù)下尖錐表面90°子午線最大橫流速度分布Fig.4 Max crossflow velocity distribution for 90° meridian line on sharp cone surface with different Mach numbers

Langtry[12]基于當?shù)販u強度和表面粗糙度提出了當?shù)鼗退贆M流轉捩拓展方案，針對γ-Reθ模型進行了橫流拓展。其橫流轉捩判據(jù)實現(xiàn)了當?shù)鼗?，脫離了傳統(tǒng)橫流雷諾數(shù)的定義形式，判據(jù)主要應用于低速流動。圖5為在Chant 2.0平臺實現(xiàn)的Langtry低速橫流轉捩模型對于亞聲速條件下NFL(2)-0415后掠翼于不同來流雷諾數(shù)下的風洞狀態(tài)，圖中Rec為弦長雷諾數(shù)，(x/c)tr為以弦長c進行無量綱化的轉捩位置，模型預測的橫流轉捩位置與實驗結果[13]較為符合，模型性能顯著優(yōu)于γ-Reθ原始模型。

圖5 后掠翼表面橫流轉捩預測Fig.5 Crossflow transition prediction on swept wing surface

經(jīng)測試，Langtry橫流轉捩模型在高超聲速橫流轉捩預測中會出現(xiàn)轉捩位置滯后或者不發(fā)生橫流轉捩的情況。如圖6所示，低速橫流模型預測的尖錐迎風面轉捩陣面相比風洞試驗[14]顯著延后，在試驗窗口內未能觀測到轉捩現(xiàn)象。轉捩模型和風洞試驗的轉捩位置見表1。

圖6 攻角α=2°時迎風面低速橫流模型計算斯坦頓數(shù)St分布與風洞試驗[14]溫升ΔT分布Fig.6 Stanton number St distribution of low-speed cross-flow model calculation and temperature rise ΔT distribution of wind tunnel experiment[14] at windward side when angle of attack α=2°

表1 轉捩位置Table 1 Transition location

2.3 構建高超聲速橫流轉捩判據(jù)

參照Langtry低速橫流轉捩判據(jù)的構建形式，以穩(wěn)定性拓展的直接數(shù)值模擬數(shù)據(jù)[15]為基礎構造了適用于高超聲速橫流轉捩的預測判據(jù)。

現(xiàn)有橫流雷諾數(shù)轉捩判據(jù)不適用于高超聲速橫流轉捩，而高超聲速橫流轉捩判據(jù)的構建受限于試驗數(shù)據(jù)的缺乏，嘗試采用線性穩(wěn)定性eN方法的拓展技術[16]對變表面粗糙度的高超聲速橢錐橫流轉捩DNS(Direct Numerical Simulation)數(shù)據(jù)[15]進行拓展。表面粗糙度是指由于表面加工導致的分布式粗糙度，粗糙高度單位為μm，在實際風洞試驗中由觸針法測量。根據(jù)“臨界橫流雷諾數(shù)ReCF,crit-表面粗糙高度h”數(shù)據(jù)組構造考慮表面粗糙度的橫流轉捩判據(jù)，并對判據(jù)進行當?shù)鼗幚?。拓展后的“ReCF,crit-h”數(shù)據(jù)一共分為9組，每組3個數(shù)據(jù)點，共計27個數(shù)據(jù)點，如圖7[17]所示?？紤]攻角和表面粗糙度，研究多重因素影響下的轉捩臨界橫流雷諾數(shù)ReCF的變化規(guī)律。

高超聲速橢錐在α=0°, 1°, 2°下的高超聲速(Ma=6)橫流轉捩臨界雷諾數(shù)呈現(xiàn)出明顯的規(guī)律性，如圖7所示。判據(jù)規(guī)律具體表現(xiàn)為：① 在相同來流狀態(tài)與飛行姿態(tài)下，對數(shù)坐標下的橢錐表面臨界橫流雷諾數(shù)ReCF,crit與表面粗糙度呈現(xiàn)出較強的線性相關性；② 有攻角條件下該線性擬合關系的斜率與截距幾乎不受來流雷諾數(shù)的影響；③ 隨攻角變化，橢錐表面的橫流強度發(fā)生變化，“ReCF-h”線性關系沿縱軸整體發(fā)生平移。

圖7 橢錐轉捩數(shù)據(jù)的穩(wěn)定性拓展[17]Fig.7 Stability extension of transition data for elliptic cone[17]

基于圖7[17]的數(shù)據(jù)與規(guī)律發(fā)現(xiàn)，參考低速橫流轉捩判據(jù)[12]的形式，采用表面粗糙度與無量綱渦量強度Hcrossflow實現(xiàn)判據(jù)當?shù)鼗瑯嬙炝丝紤]表面粗糙度的當?shù)鼗叱曀贆M流轉捩判據(jù)：

(10)

式中：lμ為粗糙度參考高度，lμ=1 μm；CCF-1和CCF-2為系數(shù)，由最小二乘法獲得，CCF-1=-9.618，CCF-2=128.33；f(Hcrossflow)為考慮橫流強度的抬升函數(shù)。抬升函數(shù)以0°攻角橢錐數(shù)據(jù)為基準進行定義：

有人說“困難像彈簧，你弱它就強”，壓力也不例外。在工作壓力面前，校長如何自我調適呢？我認為多想多做，腳踏實地地朝著目標前行；與壓力為伴，多修煉自身，多讀書，讀好書。在學校管理實踐中，不斷打磨自己，不斷提高自己的工作素質和能力，才能出色完成學校的管理工作，從而減輕工作給自己帶來的壓力。

(11)

f(Hcrossflow)=6 000|0.106 6-ΔHcrossflow|+

50 000(0.106 6-ΔHcrossflow)2

(12)

式中：ΔHcrossflow和Hcrossflow分別為橫流強度抬升項[12]和當?shù)販u量強度。

2.4 基于γ-Reθ轉捩模型的橫流拓展

橫流效應以源項形式嵌入到動量厚度雷諾數(shù)Reθ t輸運方程中：

(13)

式中：DCF為源項，定義為

(14)

3 高超聲速三維邊界層轉捩算例驗證

在高超聲速橫流轉捩判據(jù)建立完成后，對高超聲速修正的γ-Reθ轉捩模型進行了橫流拓展，將改進后模型命名為C-γ-Reθ高超聲速三維邊界層轉捩預測模型?；谠撃Ｐ蛯Ω叱曀偌忮F風洞試驗[14]和某飛行試驗的多個狀態(tài)進行數(shù)值模擬，驗證了尖錐外形三維邊界層轉捩陣面形態(tài)。

3.1 帶攻角錐的常規(guī)風洞試驗

尖錐幾何參數(shù)包括頭部鈍度、半錐角和長度。風洞試驗[14]由中國空氣動力研究與發(fā)展中心完成，轉捩陣面由紅外熱圖試驗獲得，來流馬赫數(shù)Ma=6，雷諾數(shù)Re=1.0×107，攻角α=0°～10°。數(shù)值計算網(wǎng)格量為110萬，物面第1層法向網(wǎng)格間距y+<1，來流湍流度為Tu∞≈0.8%，表面粗糙度h=0.2 μm。來流總/靜溫、總/靜壓等具體設置見文獻[14]。采用粗、中、密3套網(wǎng)格(M1、M2、M3)進行網(wǎng)格無關性驗證，網(wǎng)格量分別約為30萬、60萬、80萬。如圖8所示，熱流區(qū)別主要集中在轉捩發(fā)生后的湍流區(qū)域，M1差異較大，而M2與M3之間分布差異幾乎可忽略不計。因此，采用中等網(wǎng)格M2開展進一步多狀態(tài)計算。

圖8 α=0°時不同網(wǎng)格尖錐中心線St分布Fig.8 St distributions on center line for different grids of sharp cone at α=0°

圖9為不同攻角下風洞試驗熱流分布[14]與數(shù)值預測結果對比，攻角依次為0°、2°、4°、6°、8°、10°。風洞試驗測量為溫升ΔT分布，模型預測為表征熱流的斯坦頓數(shù)St分布，二者均可用于判斷轉捩位置。風洞試驗和CFD計算中綠色區(qū)域為層流，紅色區(qū)域為湍流，明顯分界線為轉捩。

圖9 風洞試驗[14]與C-γ-Reθ模型預測迎風面轉捩位置對比Fig.9 Transition location comparison of wind tunnel experiment[14] and C-γ-Reθ model predictions on windward surface

由于模型預測的部分狀態(tài)(α=2°, 4°)中心線轉捩位置相對靠后，為完整觀測轉捩陣面形態(tài)，設置圖8模型預測結果橫坐標范圍為[0, 900] mm，大于風洞試驗觀測窗口范圍[0, 800] mm。由圖9數(shù)值計算結果與風洞試驗結果對比可知，構建的轉捩模型能夠較為準確地預測尖錐迎風面中心線Mack模態(tài)以及側面橫流不穩(wěn)定性導致的轉捩，轉捩位置與三維轉捩陣面整體形態(tài)都得到了較好的模擬，并且能夠反映三維轉捩陣面隨攻角的變化趨勢。大攻角(α=8°, 10°)下模型預測的轉捩位置比試驗結果更為靠前，可能的原因是大攻角下迎風面邊界層較薄，表面粗糙度更易引起轉捩。

3.2 帶攻角錐的飛行試驗

某圓錐飛行試驗轉捩段的飛行馬赫數(shù)約為6，其轉捩陣面由試驗模型表面的分布式薄壁測溫傳感器測量得到。在該飛行試驗的一定彈道范圍內測得了清晰的轉捩陣面隨攻角和雷諾數(shù)的變化情況，且經(jīng)過穩(wěn)定性分析和相關測量數(shù)據(jù)分析，證明了飛行器表面橫流轉捩的存在。該橫流流動由攻角引起的周向壓力梯度造成。飛行試驗測量的轉捩數(shù)據(jù)和對應飛行狀態(tài)可以用于驗證C-γ-Reθ高超聲速三維邊界層轉捩預測模型。

如圖10所示，C-γ-Reθ轉捩模型較好地模擬了帶攻角圓錐飛行試驗的橫流轉捩陣面。在不同飛行高度對應不同飛行狀態(tài)下，圓錐±90°子午線附近區(qū)域轉捩陣面呈近似對稱狀態(tài)，轉捩模型計算的橫流轉捩陣面與飛行試驗測量結果較接近。

圖10 飛行試驗與C-γ-Reθ轉捩模型表面熱流分布Fig.10 Heat flux distributions on surface of flight test and C-γ-Reθ transition model

4 結 論

基于Chant 2.0高超計算平臺的γ-Reθ轉捩計算模塊構造了高超聲速橫流轉捩判據(jù)，進行了橫流轉捩模型拓展和驗證，改進后的高超聲速三維邊界層轉捩預測模型命名為C-γ-Reθ轉捩模型。研究工作得到如下結論：

1) 利用已有的DNS計算數(shù)據(jù)，經(jīng)過穩(wěn)定性分析的eN方法拓展，獲得了邊界層轉捩臨界橫流雷諾數(shù)ReCF,crit與表面粗糙度h和當?shù)販u量強度Hcrossflow的函數(shù)關系，構造出當?shù)鼗母叱曀贆M流轉捩判據(jù)。橫流判據(jù)的建立對于高超聲速橫流轉捩模型的建立以及橫流拓展具有重要意義。

2)C-γ-Reθ轉捩模型在多個橫流轉捩測試算例中表現(xiàn)良好，計算結果與測量結果符合較好，模型基本能夠模擬多狀態(tài)下的圓錐風洞試驗和飛行試驗的三維邊界層轉捩陣面形態(tài)。

建立的高超聲速三維邊界層轉捩模型是基于RANS方程的工程預測方法，對于橫流不穩(wěn)定性誘導轉捩的高超聲速邊界層轉捩預測具有一定的意義與價值。目前該模型僅在少數(shù)典型算例上進行了數(shù)值驗證，下一步將在更多工況、更加復雜外形的高超橫流轉捩算例中進行測試與改進，進一步考核與提升該模型在工程應用中的適用性。