可壓縮性修正湍流模型的高超聲速飛行器流場模擬研究

劉海涌，劉朝陽，劉存良

(1.第二炮兵工程大學 動力工程系，西安 710025；2. 西北工業大學 動力與能源學院，西安 710072；3.第二炮兵駐航天科技集團第四研究院軍事代表室，西安 710025)

?

可壓縮性修正湍流模型的高超聲速飛行器流場模擬研究

劉海涌1，2，劉朝陽3，劉存良2

(1.第二炮兵工程大學 動力工程系，西安 710025；2. 西北工業大學 動力與能源學院，西安 710072；3.第二炮兵駐航天科技集團第四研究院軍事代表室，西安 710025)

采用經可壓縮修正后的標準k-ε湍流模型對一類高超聲速飛行器模型的氣動特性進行了研究。計算中主要考慮了馬赫數和飛行攻角的影響。計算結果表明，通過可壓縮性修正，標準k-ε湍流模型能準確捕捉激波結構和復雜的壓力變化規律；在激波系和預壓縮面間構成的流動域內形成了高壓均勻流動，實現了飛行器設計目標；在激波結構和激波發展趨勢方面，計算結果和實驗數據符合得很好，而兩者在前體預壓縮面上的壓比值方面略有差異；經可壓縮性修正后湍流模型的計算結果要明顯優于原始湍流模型，可作為高超聲速氣動計算的有效方法。

高超聲速飛行器；湍流模型；預壓縮面；激波

0 引言

高超聲速流場具有真實氣體效應、磁流體干擾效應和力熱結構耦合效應等多重流動特性，因而相關研究需擴展到分子動力學、電磁流體力學以及流固耦合等交叉學科領域，存在巨大挑戰[1]。由于實驗研究不能提供高超聲速飛行包線全面細致的流動參數，數值模擬在其氣動特性研究中發揮了重要作用。目前，RANS方法是應用最廣，實際應用中通常對其控制方程加以修正，并對修正過程引入的新變量進行合理處理，實現方程封閉，從而減小計算量。k-ε模型是廣泛應用的一種雙方程湍流模型，能夠準確地模擬大多數流動現象，且計算效率高，在內、外流場中得到廣泛應用[2-4]。在高超聲速研究方面，k-ε湍流模型的應用也廣泛開展。Shyy等對k-ε模型的可壓縮性修正進行了研究[5]，在湍動能的生成項、耗散項和輸運項中引入了馬赫數的影響，指出速度膨脹區的額外耗散必須引入可壓縮修正。韓省思等在文獻[6-9]中對標準k-ε模型的可壓縮性修正進行了研究，這些研究考慮了Durbin和Chang實現性、Heinz湍流生成項和Sarkar可壓縮性修正，結果表明即便是復雜的高超聲速橫向射流，修正模型都在分離區內取得了與實驗結果相一致的計算結果。劉景源等對高超聲速流動計算中的兩方程湍流模型進行了對比研究[10]，研究結果表明在高超聲速橫掠流動問題中，雙方程模型的可壓縮性修正非常必要，經可壓縮性修正后的計算結果與實驗數據非常接近。耿云飛等對不同湍流模型在高超聲速流動中的適用性進行了對比評估[11]，發現對于高超聲速流動特別是大攻角流動，可壓縮性修正不可或缺。Fulco采用渦點陣方法對升力面的下洗流進行了可壓縮性修正[12]，引入可壓縮性修正因子后的計算獲得了精確結果。

綜上所述，標準k-ε湍流模型還將在保留其獨有特性(基本公式、計算穩定性和可控轉變等)的基礎上持續發展，但在高超聲速技術領域，數值模擬中有必要對k-ε湍流模型進行修正。為了獲取高超聲速飛行器詳細的流動參數，同時驗證湍流模型在高超聲速條件下的有效性，本文建立了一類高超聲速飛行器模型，采用經可壓縮性修正后的標準k-ε湍流模型開展數值模擬研究，并通過在類似模型上獲得的實驗數據對計算結果進行了驗證。

1 湍流模型、可壓縮性修正與近壁面處理

標準k-ε模型的輸運方程可表達為

Gb-ρε-YM+Sk

(1)

(2)

式中k為湍動能；ε為湍動能耗散率；Gk和Gb分別為由當地速度梯度和浮力導致的湍動能生成項；YM為可壓縮湍流中膨脹波動形成的耗散率；C1ε、C2ε、C3ε為常數，C1ε=1.44，C2ε=1.92，C3ε=0.09；σk和σε分別為與k和ε相對應的湍流普朗特數，σk=1.0，σε=1.3。

在標準k-ε模型的可壓縮性修正中考慮了膨脹耗散修正、壓力膨脹修正和激波不穩定修正三方面。其中膨脹耗散修正和壓力膨脹修正參考文獻[13]相關內容，可表達為

(3)

(4)

激波不穩定修正可表達為下述形式：

(5)

(6)

式中Pk為湍動能生成項，可壓縮性修正在YM和UDF項中進行。

k-ε模型主要適用于湍流核心區，該區域距壁面較遠。而在近壁面區域，流動參數變化梯度很大。因此，k-ε模型無法準確描述整個流動區域，需要在近壁面區域設定壁面函數。本文使用了增強型壁面函數，其主要方程可表達為

(7)

其中，

(8)

(9)

(10)

經過上述處理，整個流體域被劃分為粘性底層和充分湍流區，可提高數值模擬精度，但也對近壁面區域的網格提出了更高要求。

2 計算模型與網格劃分

設計并建立了類X-43A高超聲速飛行器計算模型，其三維結構如圖1所示，前體/進氣口、內流通道和后體的二維結構如圖2所示。飛行器前體由3段預壓縮面構成。激波生成于預壓縮面轉角處，并在設計條件下向發動機入口匯集。δ1、δ2和δ3為流動轉折角，β1、β2和β3為激波角。H0為進氣道入口高度，L0為前緣至進氣道唇口距離，l為外護罩延伸長度，Lforebody和Laftbody分別為前體和后體長度。

本文研究中未考慮飛行側滑角的影響，因此流動沿機體長度方向呈對稱形態。選用二維模型劃分網格，由于模型的法向梯度遠大于切向梯度，且壁面粘度效應要求必須在近壁面建立精細網格以捕捉粘性影響，因此在近壁面布置了15層矩形網格。為滿足第一層網格距壁面無量綱距離y+≈1的要求，同時確保采用增強型壁面函數的k-ε模型能準確捕捉粘性底層信息，模型表面第一層網格距壁面距離不超過1×10-6m。在矩形網格外圍，流動區域采用三角形網格并通過4條控制線對網格進行了局部加密。

圖1 高超聲速飛行器三維模型Fig.1 Three-dimensional model of the hypersonic vehicle

圖2 高超聲速飛行器二維結構Fig.2 Configuration of the two-dimensional model of the hypersonic vehicle

模型的邊界條件如圖3、圖4所示。計算區外邊界為壓力遠場條件，飛行器壁面為絕熱無滑移壁面。飛行器預設巡航高度25 km，所對應的計算條件：Ma=5、6、7、8，飛行攻角α=1.5°、3.5°，環境壓力和溫度分別為p0=2 600 Pa、T0=216.65 K。

圖3 計算模型的網格與邊界條件Fig.3 Grid and boundary conditions of the model

圖4 飛行器表面附近網格設置Fig.4 Grid arrangement near the surface

網格質量對數值模擬計算具有重要影響，粗糙網格會帶來較大的求解誤差，而過于精細的網格則會耗費大量的計算資源與時間，因此在數值模擬中應盡量選用精確、穩定、經濟的網格。本文選取了4種網格對模型進行了試算，其網格數量分別為約5×104、1×105、2×105和4×105。不同網格在給定位置的壓比計算結果如圖5所示。計算結果表明，當網格數量達到約2×105時，相關壓力的計算結果已達到網格無關性要求，因此本文采用總數目約220 000的較為經濟的網格開展計算。

圖5 計算結果的網格無關性檢驗Fig.5 Computational results for mesh- independent examination

3 結果與分析

3.1 計算結果的實驗驗證

張孫等建立了一種類X-43A高超聲速飛行器模型，并對其氣動性能和內部流動特性開展了實驗研究[14]。與本文所研究的飛行器模型相比，兩者的前體結構極為相似，因此文獻[14]所獲得的實驗數據可作為本文計算結果的驗證依據。圖6即為文獻[14]中的飛行器前體和發動機入口壓力紋影圖的實驗結果。

由圖6可看到，在前體預壓縮面和自由流表面形成了一系列激波系，前體預壓縮面的激波系在發動機入口附近交匯。與圖8中的計算結果對比可看出，沿飛行器前體，無論是激波系的形態，還是壓力變化趨勢都極為一致。

圖7中將沿流動通道壓力變化的實驗結果與計算結果進行了對比，其中“EXP”表示實驗數據，“CAL”表示計算結果。在前體區域，沿程壓比的實驗結果略低于計算結果，而在后體區域兩者符合較好。在實驗結果中，發動機入口位置存在明顯的壓力躍升，其波動范圍要略高于計算結果。靜壓孔壓力測量法在高超聲速條件下會造成較大的壓力損失，且在實驗過程中前體預壓縮面上存在壓力泄露，使得前體預壓縮面上壓比的實驗測量結果低于計算結果?？傮w上，計算結果與實驗結果符合較好，經可壓縮修正后的k-ε模型可為該類高超聲速飛行器結構的氣動計算提供可信數據。

(a) Ma=7,α=1.5° (b) Ma=8,α=1.5° (c) Ma=7,α=3.5° (d) Ma=8,α=3.5°

(a) Ma=7

(b) Ma=8

3.2 計算結果分析

圖8給出了計算條件為Ma=7、8和α=1.5°、3.5°時飛行器周圍壓力的計算結果。求解使用的湍流模型為經可壓縮修正后的標準k-ε模型，壓力計算結果用遠場壓力p0進行了無量綱化處理。從圖8可看到，在3處預壓縮面設計拐角處，形成了3條激波帶。沿預壓縮面長度方向，激波強度逐漸增強。在發動機內部也有激波存在，并使當地壓力急劇升高。激波在由機體和發動機壁形成的通道內多次反射并在反射點附近形成局部高壓區。氣流在激波系和預壓縮面間形成的流動空間內逐次流經激波系漸次增壓，最后形成了高壓的均勻流場。在發動機唇口位置也形成了一系列強激波系，該激波系將為飛行器帶來較大阻力。在同一飛行攻角下，隨著馬赫數由7增加至8，壓力的最大壓縮比增加了約30%，而且激波傾角略有減小。在相同馬赫數下，隨著飛行攻角由1.5°增加至3.5°，預壓縮面上的激波系強度增加，自由流表面激波系減弱，預壓縮面上的激波系交匯點向發動機入口靠近。

(a) Ma=7，α=1.5°

(b) Ma=8,α=1.5°

(c) Ma=7,α=3.5°

(d) Ma=8,α=3.5°

圖9給出了不同馬赫數和攻角下前體預壓縮面、發動機內上表面和后體表面上的壓力數據計算結果。所使用的湍流模型為經可壓縮修正后的標準k-ε模型。由圖9可看到，在前體預壓縮面上存在3處清晰的壓力躍升，躍升位置與預壓縮面轉角位置符合良好。在發動機入口位置，也存在一次壓力躍升，這表明在該位置形成了強激波系。發動機內部的壓力一直保持在較高水平，而在發動機尾部壓力突然降低。發動機內部的激波反射使其內部壓力形成雙峰分布形態。隨著馬赫數的增加，前體預壓縮面、發動機內上表面和后體表面上的壓力提高，且發動機內部壓力峰值形態發生變化，較高的壓力峰值隨馬赫數的增加而不斷向發動機入口靠近。這表明在發動機入口位置，激波強度不斷提升。與馬赫數相比，飛行攻角的增加同樣會在總體上提升表面壓力，但對壓力峰值形態未產生明顯影響。

(a) Ma=5～8,α=1.5°

(b) Ma=5～8,α=3.5°

3.3 與未進行可壓縮性修正標準k-ε模型計算結果的對比

圖10給出了普通標準k-ε模型在Ma=7、8和α=1.5°、3.5°下飛行器周圍的壓力分布的計算結果。與經過可壓縮性修正后的標準k-ε模型的計算結果相比，在未進行可壓縮性修正的情況下，最大壓比的計算結果要明顯低于前者的計算值。在前體預壓縮區域，普通標準k-ε模型未能捕捉到轉折角處的激波系。圖11給出了普通標準k-ε模型在不同飛行馬赫數和攻角下，前體預壓縮面、發動機內上表面和后體表面上的壓力變化的計算結果。與圖9中經過可壓縮修正后的湍流模型的計算結果相比，普通模型未能捕捉到前緣處的壓力變化；而在前體預壓縮面上，壓力呈線性升高趨勢，未能捕捉到轉折角處的壓力躍升；在發動機出口位置，壓力出現劇烈波動，未實現平滑過渡；在整個流動通道中，沿程壓比值整體水平較低。

上述結果表明，在高超聲速條件下，經過可壓縮性修正后的標準k-ε模型能更為準確的反映壓力變化，捕捉激波能力增強，計算精度和穩定性提高，能夠提供更為準確的計算結果。

(a) Ma=7，α=1.5°

(b) Ma=8,α=1.5°

(c) Ma=8,α=3.5°

(d) Ma=8,α=3.5°

(a) Ma=5～8,α=1.5°

(b) Ma=5～8,α=3.5°

4 結論

(1)計算結果表明飛行器模型實現了初始設計目標。在激波系和預壓縮面間構成的流動域內形成了高壓均勻流動，為飛行器提供了很高升力。在發動機內部，流動減速并形成高壓環境，有利于發動機點火和穩定燃燒。經可壓縮性修正后，標準k-ε模型能準確捕捉激波結構和壓力變化特性，在前體預壓縮面上壓力表現出階梯躍升，在發動機內部壓力躍升并保持較高壓力水平。

(2)實驗模型中獲取的壓力紋影圖提供了有價值的激波結構和壓力變化趨勢信息。在激波結構方面，實驗數據與計算結果良好符合。但在前體預壓縮面上的壓力變化方面，沿程壓比的計算結果高于實驗結果，實驗中測量方法造成的壓力損失是造成上述差別的主要原因?？傮w上，計算結果與實驗結果較好符合，經可壓縮修正后的k-ε模型可作為該類高超聲速飛行器結構氣動計算的有效方法。

[1] 王江峰,伍貽兆,季衛棟,等.高超聲速復雜氣動問題數值方法研究進展[J].航空學報,2015,36(1)：159-175.

[2] 謝春艷,張敬秋,王洋.改進型k-ε湍流模型在三維側噴流數值仿真中的應用[J].導彈與航天運載技術,2014,330(1):56-60.

[3] 劉薇,宋國萍,褚雙磊,韓博.基于ANSYS FLUENT的近場翼尖渦數值模擬與分析[J].飛行力學,2015,33(2):16-21.

[4] 王樂,宋衛東,楊曉霖.基于可壓縮修正k-ε模型的彈丸氣動計算[J].彈箭與制導學報,2013(1):149-151.

[5] Shyy W,Krishnamurty V S.Compressibility effects in modeling complex turbulent flows[J].Aerospace Science,1997,33:587-645.

[6] 韓省思,葉桃紅,朱雯明,等.一個新的可壓縮性修正的k-ε模型[J].空氣動力學學報,2009,27(6):677-682.

[7] 韓省思,葉桃紅,朱旻明,等.k-ε湍流模型可壓縮性修正在超音速混合層中的應用研究[J].工程熱物理學報,2007,28(6):1053-1055.

[8] 韓省思,劉亮,葉桃紅,等.三種湍流對Scramjet燃燒室冷態流場的數值模擬[J].中國科學技術大學學報,2009,39(3):321-325.

[9] 韓省思,葉桃紅,朱旻明,等.激波不穩定性效應的k-ε可壓縮修正湍流模型[J].科學通報,2008,53(22);2722-2729.

[10] 劉景源,李椿萱.高超聲速二方程湍流模型的數值模擬對比[J].北京航空航天大學學報,2007,33(10):1131-1135.

[11] 耿云飛,閻超,徐晶磊,等.高超聲速流動湍流模式評估[J].北京航空航天大學學報,2011,37(8):907-911.

[12] Fulco V,Roelof V.A new compressibility correction method to predict aerodynamic interaction between lifting surfaces[R].AIAA 2013-4299.

[13] Sarkar S,Erlebacher G,Hussaini M Y,et al.The analysis and modeling of dilatational terms in compressible turbulence[J].Journal of Fluid Mechanics,1991,227:473-493.

[14] 張孫.類X-43A高超聲速飛行器氣動力特性及其全流道流動特征的研究[D].南京：南京航空航天大學,2007.

(編輯：呂耀輝)

Research on the application of a turbulence models with compressibility effects correction for a hypersonic vehicle

LIU Hai-yong1，2，LIU Chao-yang3，LIU Cun-liang2

(1.School of Engine and Energy，Second Artillery Engineering Univ.，Xi'an 710025，China；2.School of Power and Energy，Northwestern Polytechnical Univ.，Xi'an 710072，China；3.Military Representative Department of the 2nd Artillery Force Stationed in the 4th Research Academy of CASC，Xi'an 710025，China)

A series of numerical analyses have been performed to investigate the flow characteristics of a hypersonic vehicle model. Flow fields corresponding to the Mach number and attack angle were analyzed by solving the Reynolds averaged Navier-stokes equations with the standk-εturbulence model in which the compressibility effects were considered. With the correction of compressibility effects，the standardk-εturbulence model can capture delicate shock wave structures and complex pressure variation characteristics,and their calculation results are in highly agreement. The calculated shock wave structures and pressure variation tendency agree well with experimental data，except for some differences in the variation of the static pressure on the pre-compressed surface. Compared with the originalk-εturbulence model,the improved model with compressibility effects correction obtains more accurate data,implying the validity and effectiveness of the numerical approach for analyzing relevant flow field.

hypersonic vehicle；turbulence model；pre-compressed surface；shock wave

2015-01-03；

：2015-05-17。

國家自然科學基金(51206180)；陜西省自然科學基金(2014JQ7276)。

劉海涌(1981—)，男，博士后，主要從事飛行器及發動機高溫部件熱防護技術研究。E-mail:helian_xicheng@163.com

V430

A

1006-2793(2015)06-0770-06

10.7673/j.issn.1006-2793.2015.06.004