孔維萱，閻 超，趙 瑞

（北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191）

0 引 言

在高超聲速飛行條件下，邊界層轉(zhuǎn)捩直接關(guān)系到飛行器摩擦阻力、熱交換及流動(dòng)的分離位置等，若能準(zhǔn)確預(yù)測(cè)轉(zhuǎn)捩位置或延遲轉(zhuǎn)捩的發(fā)生，可有效改進(jìn)飛行器氣動(dòng)性能，大大降低燃料消耗，使熱防護(hù)設(shè)計(jì)更為靈活。另一方面發(fā)動(dòng)機(jī)唇口流動(dòng)品質(zhì)是發(fā)動(dòng)機(jī)工作的重要參數(shù)之一，流動(dòng)為湍流態(tài)時(shí)燃?xì)獬煞謸交斐浞?，保證燃燒效率；正確控制流動(dòng)狀態(tài)，強(qiáng)制流動(dòng)發(fā)生轉(zhuǎn)捩已成為高超聲速一體化飛行器設(shè)計(jì)的重點(diǎn)難題。研究邊界層轉(zhuǎn)捩及其預(yù)測(cè)方法，對(duì)于飛行器氣動(dòng)設(shè)計(jì)和進(jìn)氣道一體化設(shè)計(jì)具有重要的工程意義。

風(fēng)洞實(shí)驗(yàn)一直是研究流動(dòng)轉(zhuǎn)捩的重要手段，近年來靜風(fēng)洞的研究成果［1－3］成為模擬真實(shí)飛行狀態(tài)流動(dòng)轉(zhuǎn)捩和驗(yàn)證數(shù)值方法的重要依據(jù)。隨著計(jì)算流體力學(xué)的發(fā)展，數(shù)值試驗(yàn)已逐漸成為與風(fēng)洞實(shí)驗(yàn)同等重要的流體力學(xué)研究方法。其中基于模式理論發(fā)展的轉(zhuǎn)捩模式，以其穩(wěn)定性和經(jīng)濟(jì)性成為最具有工程應(yīng)用前景的轉(zhuǎn)捩模擬和預(yù)測(cè)方法。Mayle等［4］提出非湍流脈動(dòng)能概念，?；D(zhuǎn)捩前的流動(dòng)脈動(dòng)，隨后 Warren等［5］將其耦合在湍流模式中，形成了一方程的轉(zhuǎn)捩模式。Menter等［6］提出了一種避免求解積分量的基于當(dāng)?shù)刈兞康霓D(zhuǎn)捩模式，此模式在眾多亞聲速轉(zhuǎn)捩流動(dòng)中表現(xiàn)優(yōu)秀；值得一提的是該模式特別針對(duì)分離導(dǎo)致的轉(zhuǎn)捩對(duì)間歇因子提出了修正，并在部分算例中得到良好的效果。符松和王亮［7］提出了適用于超聲速流動(dòng)的轉(zhuǎn)捩模式，該模式統(tǒng)一求解非湍流脈動(dòng)和湍動(dòng)能，并以脈動(dòng)與平均流動(dòng)的關(guān)系作為轉(zhuǎn)捩起始的判據(jù)，考慮多種不穩(wěn)定擾動(dòng)模態(tài)對(duì)了流動(dòng)的影響。

本文采用γ－Reθ模式對(duì)超聲速、高超聲速流動(dòng)的典型模型進(jìn)行計(jì)算，研究了來流馬赫數(shù)、來流雷諾數(shù)、攻角狀態(tài)和頭部鈍化半徑對(duì)流動(dòng)轉(zhuǎn)捩的影響，就γ－Reθ模式針對(duì)上述流動(dòng)參數(shù)變化的適用性進(jìn)行討論。

1 計(jì)算方法

1．1 控制方程

計(jì)算采用有限體積法求解可壓縮全N－S方程，無粘通量由Roe的FDS格式求解，粘性通量采用中心差分格式進(jìn)行離散，時(shí)間推進(jìn)采用LU－SGS隱式方法。

1．2 轉(zhuǎn)捩模式

本文采用γ－Reθ模式模擬流動(dòng)的轉(zhuǎn)捩。γ－Reθ模式是Menter及合作者［6］提出的轉(zhuǎn)捩模式，引入渦雷諾數(shù)概念將邊界層相關(guān)量的求解當(dāng)?shù)鼗?，建立輸運(yùn)方程求解間歇因子和臨界動(dòng)量厚度雷諾數(shù)，在SST湍流模式的基礎(chǔ)上實(shí)現(xiàn)了對(duì)于流動(dòng)轉(zhuǎn)捩的計(jì)算模擬和預(yù)測(cè)。

1．2．1 當(dāng)?shù)乩字Z數(shù)輸運(yùn)方程

轉(zhuǎn)捩起始位置的求解對(duì)于轉(zhuǎn)捩計(jì)算和預(yù)測(cè)至關(guān)重要，γ－Reθ模式中求解“臨界動(dòng)量厚度雷諾數(shù)~Reθt”輸運(yùn)方程從而給出轉(zhuǎn)捩起始位置的相關(guān)信息，與由經(jīng)驗(yàn)公式或轉(zhuǎn)捩判據(jù)確定轉(zhuǎn)捩起始位置的方法相比，此方法能夠更多地反映當(dāng)?shù)亓鲃?dòng)特性。當(dāng)?shù)乩字Z數(shù)輸運(yùn)方程表示為：

其中生成項(xiàng)為：

轉(zhuǎn)捩動(dòng)量厚度雷諾數(shù)Reθt由來流湍流度和壓力梯度等參數(shù)擬合的經(jīng)驗(yàn)公式，表示為：

模式常數(shù)為：

1．2．2 間歇因子輸運(yùn)方程

將間歇因子作為權(quán)值耦合層流與湍流區(qū)域的計(jì)算，采用輸運(yùn)方程求解間歇因子，間歇因子輸運(yùn)方程為：

間歇因子的生成項(xiàng)和耗散項(xiàng)分別表示為：

決定間歇因子增長的關(guān)鍵參數(shù)是Fonset，由下式確定：

其中Reθc表示間歇因子開始增長處的動(dòng)量厚度雷諾數(shù)，可表示為的函數(shù)：

確定間歇因子生成項(xiàng)的另一關(guān)鍵參數(shù)為Flength，表示轉(zhuǎn)捩區(qū)的長度，同樣以的函數(shù)形式給出如下：

其它模式常數(shù)為：

γ壁面法向通量為0；為抑制入口處湍動(dòng)能的耗散現(xiàn)象，γ設(shè)置為1。

1．2．3 模式最終形式

間歇因子體現(xiàn)在原SST模式中湍動(dòng)能生成項(xiàng)和耗散項(xiàng)中，γ－Reθ轉(zhuǎn)捩模式的最終形式如下：

2 數(shù)值實(shí)驗(yàn)

2．1 數(shù)值格式對(duì)轉(zhuǎn)捩模式的影響

數(shù)值格式本身的特性是影響數(shù)值模式的關(guān)鍵因素之一，這里選取超聲速平板流動(dòng)算例，采用二階格式包括Roe的FDS格式（分別組合min－mod和superbee限制器）、AUSM＋格式和五階 WENO格式進(jìn)行計(jì)算，研究不同格式對(duì)轉(zhuǎn)捩模擬的影響。超聲速平板流動(dòng)條件為：來流馬赫數(shù)4．5，來流雷諾數(shù)6．433×106／m，來流溫度61．1K，來流湍流度0．1%。計(jì)算域設(shè)置為1m×0．33m，保證法向第一層網(wǎng)格y＋＜1。入口給定速度、壓力、k和ω。

圖1給出壁面摩阻系數(shù)沿流向分布情況，比較各種格式計(jì)算結(jié)果可知：同為二階格式時(shí)，限制器的選擇對(duì)計(jì)算結(jié)果的影響明顯，其中壓縮性較強(qiáng)的superbee限制器計(jì)算得到的轉(zhuǎn)捩位置明顯提前于耗散性較大的min－mod限制器結(jié)果；不同格式計(jì)算結(jié)果差異不大；相對(duì)于二階格式，高階格式對(duì)流動(dòng)結(jié)構(gòu)的刻畫更為精確，數(shù)值耗散小，但計(jì)算穩(wěn)定性差，并對(duì)存儲(chǔ)量和計(jì)算量提出了更高的要求。對(duì)于使用γ－Reθ模式的轉(zhuǎn)捩預(yù)測(cè)，計(jì)算格式的影響主要體現(xiàn)在格式的耗散性上，耗散較小的計(jì)算格式得到的轉(zhuǎn)捩位置較耗散性大的格式提前，是由于γ－Reθ轉(zhuǎn)捩模式中控制間歇因子γ生成的函數(shù)Fonset與粘性負(fù)相關(guān)，即體現(xiàn)為數(shù)值粘性越小轉(zhuǎn)捩發(fā)生位置越提前。在原有模式參數(shù)設(shè)置的基礎(chǔ)上，耗散性較大的格式計(jì)算結(jié)果與DNS數(shù)據(jù)吻合較好，綜合考慮計(jì)算穩(wěn)定性等因素，下文各算例均采用Roe格式組合min－mod限制器的空間離散方法進(jìn)行計(jì)算。

圖1 馬赫4．5平板壁面摩阻系數(shù)分布Fig．1 Skin－friction coefficient distribution for M4．5plate

2．2 算例描述

應(yīng)用γ－Reθ模式計(jì)算以下算例：（1）馬赫3．5、半錐角為5°尖錐；（2）馬赫5．91、半錐角5°尖錐；（3）馬赫5．91、半錐角5°裙錐和（4）馬赫7．16、半錐角7°尖錐。其中馬赫3．5尖錐計(jì)算了6種雷諾數(shù)條件，討論來流雷諾數(shù)對(duì)流動(dòng)轉(zhuǎn)捩的影響。針對(duì)馬赫5．91尖錐和裙錐的計(jì)算，討論不同攻角狀態(tài)對(duì)轉(zhuǎn)捩的影響。7°半錐角尖錐算例中給出了三種頭部鈍化半徑的計(jì)算結(jié)果，討論頭部鈍化半徑對(duì)流動(dòng)轉(zhuǎn)捩的影響。

各算例流動(dòng)參數(shù)和計(jì)算設(shè)置列于表1，需要說明的是，針對(duì)各算例分別設(shè)置湍動(dòng)能及比耗散率的入口條件，以盡量保證計(jì)算來流湍流度與實(shí)驗(yàn)或參考文獻(xiàn)相同。

表1 各算例流動(dòng)參數(shù)及網(wǎng)格設(shè)置Table 1 Flow condition and grid refinement

以下給出各算例計(jì)算結(jié)果和討論。

2．3 雷諾數(shù)對(duì)轉(zhuǎn)捩的影響

計(jì)算長為550mm、半錐角為5°的尖錐，討論不同來流雷諾數(shù)條件下流動(dòng)轉(zhuǎn)捩情況。流動(dòng)條件參照實(shí)驗(yàn)［8］設(shè)置，來流馬赫數(shù)為3．5，來流湍流度為0．1%，來流雷諾數(shù)分別為：2．69×107、3．70×107、4．69×107、5．66×107、6．59×107和7．49×107。計(jì)算網(wǎng)格設(shè)置（流向×法向×周向，下同）為121×121×81，法向第一層網(wǎng)格高度為1．0E－6m；流動(dòng)入口給定速度、溫度、壓力、k和ω，壁面提無滑移、絕熱條件，出口提遠(yuǎn)場(chǎng)條件。

圖2 馬赫3．5尖錐計(jì)算結(jié)果Fig．2 Skin－friction coefficient distribution for M3．5cone

圖3 馬赫3．5尖錐計(jì)算結(jié)果Fig．3 ReTvs．Re／m from M3．5cone results

圖2 給出馬赫3．5圓錐計(jì)算結(jié)果，圖中實(shí)線（由左側(cè)坐標(biāo)標(biāo)示）表示壁面摩阻系數(shù)，離散點(diǎn)（由右側(cè)坐標(biāo)標(biāo)示）表示壁面恢復(fù)溫度風(fēng)洞實(shí)驗(yàn)值，恢復(fù)溫度定義為r＝（Taw－Te）／（T0－Te）。隨來流雷諾數(shù)的增大，γ－Reθ模式計(jì)算結(jié)果呈現(xiàn)出：轉(zhuǎn)捩起始位置逐漸向上游移動(dòng)、轉(zhuǎn)捩后摩阻系數(shù)的峰值逐漸增大的變化趨勢(shì)。在較小的來流雷諾數(shù)條件下，模式計(jì)算的轉(zhuǎn)捩起始位置明顯比實(shí)驗(yàn)提前，隨著來流雷諾數(shù)的增大，這一差異逐漸減小。

圖3比較了飛行實(shí)驗(yàn)、風(fēng)洞實(shí)驗(yàn)和γ－Reθ模式計(jì)算的轉(zhuǎn)捩雷諾數(shù)隨來流雷諾數(shù)的變化趨勢(shì)。γ－Reθ模式計(jì)算結(jié)果與飛行實(shí)驗(yàn)變化趨勢(shì)相同，可以證實(shí)γ－Reθ模式能夠正確反映流動(dòng)轉(zhuǎn)捩隨來流雷諾數(shù)的變化；但模式計(jì)算的轉(zhuǎn)捩雷諾數(shù)整體小于飛行實(shí)驗(yàn)結(jié)果，說明γ－Reθ模式未考慮可壓縮性對(duì)流動(dòng)的穩(wěn)定作用，在超聲速流動(dòng)（Ma～3．5左右）條件下預(yù)測(cè)的轉(zhuǎn)捩位置提前，轉(zhuǎn)捩雷諾數(shù)較小。而風(fēng)洞實(shí)驗(yàn)結(jié)果只在來流雷諾數(shù)較小時(shí)與飛行實(shí)驗(yàn)結(jié)果吻合良好，而隨著來流雷諾數(shù)的增大，風(fēng)洞自身的噪聲很大程度上影響了流動(dòng)轉(zhuǎn)捩的結(jié)果。

2．4 攻角對(duì)轉(zhuǎn)捩的影響

計(jì)算裙錐和圓錐0°～4°攻角流動(dòng)，討論小攻角對(duì)流動(dòng)轉(zhuǎn)捩的影響。圓錐長635mm、半錐角為5°；裙錐前半部為長254mm、半錐角為5°尖錐，后半部為曲率半徑2364mm的尾裙，特別設(shè)計(jì)的曲率半徑使流動(dòng)在尾裙區(qū)域具有逆壓梯度，并保持邊界層厚度不變。頭部鈍化半徑均為0．8mm，裙錐模型幾何外形如圖4所示。

圖4 尖錐裙幾何外形示意圖Fig．4 Geometry for flared－cone

流動(dòng)參數(shù)參照風(fēng)洞實(shí)驗(yàn)［9－10］設(shè)置，來流馬赫數(shù)為5．91，來流雷諾數(shù)為 9．348E6／m，來流湍流度為0．1%，來流溫度為56．2K。網(wǎng)格設(shè)置為121×121×41，法向第一層網(wǎng)格高度1．0E－6m；邊界條件設(shè)置與上一算例相同。

圖5為馬赫5．91尖錐裙流動(dòng)計(jì)算結(jié)果。圖5（a）給出0°攻角時(shí)對(duì)稱面壓力等值線分布，除尖錐前緣產(chǎn)生的斜激波外，尾裙區(qū)域存在一簇較弱的壓縮波。圖5（b）給出0°攻角母線、2°攻角迎風(fēng)、背風(fēng)子午線壓力系數(shù)分布。γ－Reθ模式計(jì)算的壓力系數(shù)與實(shí)驗(yàn)值吻合良好，最大差異在裙錐尾部；圖5（b）中可以看到，實(shí)驗(yàn)結(jié)果2°攻角背風(fēng)面尾裙中部位置出現(xiàn)壓力系數(shù)的局部下降，這一流動(dòng)未能由γ－Reθ模式計(jì)算再現(xiàn)。實(shí)驗(yàn)測(cè)得0°攻角轉(zhuǎn)捩發(fā)生在x＝371mm位置，2°攻角背風(fēng)側(cè)轉(zhuǎn)捩位置提前至x＝310mm，迎風(fēng)側(cè)基本保持層流狀態(tài)不發(fā)生轉(zhuǎn)捩。計(jì)算結(jié)果與實(shí)驗(yàn)相同，圖5（c）給出0°和2°攻角壁面法向第一層網(wǎng)格間歇因子分布，0°攻角間歇因子在x＝370mm開始增長；2°攻角時(shí)背風(fēng)子午線間歇因子開始增長位置提前至x＝300mm，迎風(fēng)子午線間歇因子保持在很小數(shù)值。

圖5 M5．91裙錐流動(dòng)計(jì)算結(jié)果Fig．5 Results for M5．91flared－cone

圖6 給出馬赫5．91直圓錐計(jì)算結(jié)果。圖6（a）為0°攻角壁面熱流系數(shù)與實(shí)驗(yàn)對(duì)比，計(jì)算結(jié)果與實(shí)驗(yàn)值基本吻合，但轉(zhuǎn)捩起始位置較為提前。圖6（b）給出1°、2°和4°攻角壁面法向第一層網(wǎng)格間歇因子分布。同裙錐計(jì)算結(jié)果相似，當(dāng)攻角由0°變?yōu)?°時(shí)，背風(fēng)側(cè)轉(zhuǎn)捩位置明顯提前，迎風(fēng)側(cè)轉(zhuǎn)捩位置推后；隨著攻角的增加，迎、背風(fēng)側(cè)轉(zhuǎn)捩位置逐漸向下游移動(dòng)，且周向轉(zhuǎn)捩線傾斜程度逐漸增大。

分析本節(jié)兩算例的計(jì)算結(jié)果，γ－Reθ模式通過求解當(dāng)?shù)乩字Z數(shù)輸運(yùn)方程得到轉(zhuǎn)捩起始位置的相關(guān)信息，并結(jié)合當(dāng)?shù)乩字Z數(shù)確定間歇因子，計(jì)算能夠一定程度上反應(yīng)邊界層的三維特性，在未引入橫向不穩(wěn)定擾動(dòng)?；那闆r下，能夠正確預(yù)測(cè)小攻角對(duì)流動(dòng)轉(zhuǎn)捩的影響。

圖6 馬赫5．91直圓錐計(jì)算結(jié)果Fig．6 Results for M5．91cone

2．5 頭部鈍化對(duì)轉(zhuǎn)捩的影響

最后討論頭部鈍化半徑對(duì)流動(dòng)轉(zhuǎn)捩的影響，選取長2300mm、半錐角為7°的圓錐進(jìn)行研究，頭部鈍化半徑（Rn／mm）分別取為2．5、5．0和6．35，以頭部鈍化半徑為特征長度的雷諾數(shù)（ReRn×10－4）分別為2．55、5．1和6．477。流動(dòng)條件參照實(shí)驗(yàn)［11］給定，來流馬赫數(shù)為7．16，來流湍流度為0．675%，來流溫度為231．7K，壁面提絕熱條件。計(jì)算網(wǎng)格設(shè)置為121×121×41，法向第一層網(wǎng)格高度為1．0E－6m；邊界條件設(shè)置與上一算例相同。

圖7為馬赫7．16圓錐流動(dòng)計(jì)算結(jié)果，圖7（a）和（b）給出頭部鈍化半徑分別為2．5mm和5．0mm時(shí)壁面熱流值與實(shí)驗(yàn)對(duì)比，熱流峰值的位置和大小等計(jì)算結(jié)果與實(shí)驗(yàn)值基本相同，但同樣存在預(yù)測(cè)轉(zhuǎn)捩起始位置較早的問題。頭部鈍化半徑的變化對(duì)于轉(zhuǎn)捩的影響，主要源于壓力梯度的變化；由于球頭與錐體相接處存在曲率半徑的突變，在該位置產(chǎn)生較大逆壓梯度，從而影響流動(dòng)轉(zhuǎn)捩的發(fā)生。圖7（c）給出圓錐前緣局部壓力分布及壁面壓力梯度（?p／?s）分布，可見隨著頭部鈍化半徑的增加，逆亞梯度峰值位置逐漸向下游移動(dòng)；壓力梯度的變化對(duì)流動(dòng)轉(zhuǎn)捩的影響表現(xiàn)為，隨著頭部鈍化半徑的增加，轉(zhuǎn)捩位置向下游移動(dòng)，轉(zhuǎn)捩后的熱流峰值有一定程度的降低。模式能夠正確模擬頭部半徑對(duì)流動(dòng)轉(zhuǎn)捩的影響，針對(duì)壓力梯度進(jìn)行了修正提高了對(duì)于此類問題的適用性。

圖7 馬赫7．16尖錐計(jì)算結(jié)果Fig．7 Results for M7．16cone

3 結(jié) 論

（1）將γ－Reθ模式應(yīng)用于超聲速、高超聲速邊界層轉(zhuǎn)捩問題的計(jì)算，討論了來流雷諾數(shù)、攻角狀態(tài)和頭部鈍化半徑等關(guān)鍵參數(shù)變化時(shí)的模式性能。數(shù)值實(shí)驗(yàn)表明γ－Reθ模式具有一定的預(yù)測(cè)超聲速、高超聲速轉(zhuǎn)捩的能力。

（2）馬赫數(shù)3至7范圍內(nèi)的圓錐類流動(dòng)中，γ－Reθ模式對(duì)轉(zhuǎn)捩位置的預(yù)測(cè)普遍提前，來流雷諾數(shù)較小時(shí)尤為明顯。隨頭部鈍化半徑的增大，γ－Reθ模式計(jì)算得到的轉(zhuǎn)捩起始位置后移，對(duì)于壓力梯度的變化敏感。

（3）γ－Reθ模式未針對(duì)橫流不穩(wěn)定性進(jìn)行特別模化，但能夠正確模擬馬赫6圓錐或裙錐小攻角狀態(tài)下的邊界層轉(zhuǎn)捩。相比于0°攻角，小攻角時(shí)背風(fēng)側(cè)轉(zhuǎn)捩位置提前、迎風(fēng)側(cè)轉(zhuǎn)捩延遲；隨攻角的增加，周向轉(zhuǎn)捩線變陡峭并向下游移動(dòng)。

［1］ SCHNEIDER S．A Review of hypersonic boundary layer stability experiments in a quiet Mach 6wind tunnel［R］．AIAA 1997－1819，1997．

［2］ LACHOWICZ J，CHOKANI N，et al．Boundary－layer stability measurements in a hypersonic quiet tunnel［J］．AIAA J．，34（12）：2496－2500，1996．

［3］ CHEN F，MALIK M，et al．Boundary－layer transition on a cone and flat plat at Mach 3．5［J］．AIAA J．，27（6）：687－693，1989．

［4］ MAYLE R，SCHULZ A．The path to predicting bypass transition［J］．J．Turbomach，1997，119：405－411．

［5］ WARREN E，HARRIS J，HASSAN H．Transition model for high－speed flow［J］．AIAA J．，1995，33（8）：1391－1397．

［6］ LANGTRY R，MENTER F．Correlation－based transition modeling for unstructured parallelized computational fluid dy－namics codes［J］．AIAA J．，2009，47（12）：2894－2906．

［7］ 王亮，符松．一種適用于超音速邊界層的湍流轉(zhuǎn)捩模式［J］．力學(xué)學(xué)報(bào)，2009，41（2）：162－168．

［8］ CHEN F，MALIK M，et al．Boundary－layer transition on a cone and flat plat at Mach 3．5［J］．AIAA J，1989，27（6）：687－693．

［9］ DOGGETT G，CHOKANI N，WILKINSON S．Effects of angle of attack on hypersonic boundary layer stability in a quiet wind tunnel［J］．AIAA J，1997，35（3）：464－470．

［10］HORVATH T，BERRY S，et al．Boundary layer transition on slender cones in conventional and low disturbance Mach 6wind tunnels［R］．AIAA 2002－2742．

［11］JOHNSON H，ALBA C，et al．Boundary layer stability analysis to support the HiFiRE transition experiment［R］．AIAA 2007－0311．