C- -Re 高超聲速三維邊界層轉(zhuǎn)捩預(yù)測模型

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

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

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

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

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

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

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)技術(shù)進行并行計算。

1.2 基礎(chǔ)轉(zhuǎn)捩模型

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

(1)

(2)

(3)

(4)

(5)

(6)

(7)

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

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

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

2 構(gòu)建橫流轉(zhuǎn)捩模型

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

2.1 γ-Reθ轉(zhuǎn)捩模型的橫流預(yù)測不適用性

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

(8)

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

圖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%的位置(壁面較遠(yuǎn)處)對應(yīng)的壁面法向距離。

天津大學(xué)趙磊[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與來流相關(guān)。

通過超聲速、高超聲速尖錐數(shù)值計算發(fā)現(xiàn)，最大橫流速度Wmax與來流馬赫數(shù)緊密相關(guān)，并且變化幅度較大。圖4為固定雷諾數(shù)、來流馬赫數(shù)Ma=3～7情況下，尖錐側(cè)面90°子午線最大橫流速度沿軸向的分布。對于橫流速度的馬赫數(shù)相關(guān)性，趙磊在鈍頭后掠平板流動中也有類似發(fā)現(xiàn)[9]。由圖4可知，尖錐表面橫流強度大部分情況下超過12%，該強度也超過了Reed和Haynes[11]提出的橫流準(zhǔn)則適用范圍,即最大橫流速度Ucf,max與邊界層外緣速度Ue比值不超過8%(Ucf,max/Ue≤8%)。固定來流雷諾數(shù)前提下，最大橫流速度具有較強的馬赫數(shù)相關(guān)性，導(dǎo)致采用單一臨界值的橫流雷諾數(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]基于當(dāng)?shù)販u強度和表面粗糙度提出了當(dāng)?shù)鼗退贆M流轉(zhuǎn)捩拓展方案，針對γ-Reθ模型進行了橫流拓展。其橫流轉(zhuǎn)捩判據(jù)實現(xiàn)了當(dāng)?shù)鼗?，脫離了傳統(tǒng)橫流雷諾數(shù)的定義形式，判據(jù)主要應(yīng)用于低速流動。圖5為在Chant 2.0平臺實現(xiàn)的Langtry低速橫流轉(zhuǎn)捩模型對于亞聲速條件下NFL(2)-0415后掠翼于不同來流雷諾數(shù)下的風(fēng)洞狀態(tài)，圖中Rec為弦長雷諾數(shù)，(x/c)tr為以弦長c進行無量綱化的轉(zhuǎn)捩位置，模型預(yù)測的橫流轉(zhuǎn)捩位置與實驗結(jié)果[13]較為符合，模型性能顯著優(yōu)于γ-Reθ原始模型。

圖5 后掠翼表面橫流轉(zhuǎn)捩預(yù)測Fig.5 Crossflow transition prediction on swept wing surface

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

圖6 攻角α=2°時迎風(fēng)面低速橫流模型計算斯坦頓數(shù)St分布與風(fēng)洞試驗[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 轉(zhuǎn)捩位置Table 1 Transition location

2.3 構(gòu)建高超聲速橫流轉(zhuǎn)捩判據(jù)

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

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

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

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

基于圖7[17]的數(shù)據(jù)與規(guī)律發(fā)現(xiàn)，參考低速橫流轉(zhuǎn)捩判據(jù)[12]的形式，采用表面粗糙度與無量綱渦量強度Hcrossflow實現(xiàn)判據(jù)當(dāng)?shù)鼗?，?gòu)造了考慮表面粗糙度的當(dāng)?shù)鼗叱曀贆M流轉(zhuǎn)捩判據(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ù)為基準(zhǔn)進行定義：

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

(11)

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

50 000(0.106 6-ΔHcrossflow)2

(12)

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

2.4 基于γ-Reθ轉(zhuǎn)捩模型的橫流拓展

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

(13)

式中：DCF為源項，定義為

(14)

3 高超聲速三維邊界層轉(zhuǎn)捩算例驗證

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

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

尖錐幾何參數(shù)包括頭部鈍度、半錐角和長度。風(fēng)洞試驗[14]由中國空氣動力研究與發(fā)展中心完成，轉(zhuǎn)捩陣面由紅外熱圖試驗獲得，來流馬赫數(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。來流總/靜溫、總/靜壓等具體設(shè)置見文獻(xiàn)[14]。采用粗、中、密3套網(wǎng)格(M1、M2、M3)進行網(wǎng)格無關(guān)性驗證，網(wǎng)格量分別約為30萬、60萬、80萬。如圖8所示，熱流區(qū)別主要集中在轉(zhuǎn)捩發(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為不同攻角下風(fēng)洞試驗熱流分布[14]與數(shù)值預(yù)測結(jié)果對比，攻角依次為0°、2°、4°、6°、8°、10°。風(fēng)洞試驗測量為溫升ΔT分布，模型預(yù)測為表征熱流的斯坦頓數(shù)St分布，二者均可用于判斷轉(zhuǎn)捩位置。風(fēng)洞試驗和CFD計算中綠色區(qū)域為層流，紅色區(qū)域為湍流，明顯分界線為轉(zhuǎn)捩。

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

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

3.2 帶攻角錐的飛行試驗

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

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

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

4 結(jié) 論

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

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

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

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