高焓風洞一體化數(shù)值模擬及其對氣動特性影響

傅楊奧驍, 董維中, 丁明松, 劉慶宗, 江 濤

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

0 引 言

高超聲速飛行器在再入和滑翔過程中，如果飛行速度很高(一般認為馬赫數(shù)10以上)，會出現(xiàn)高溫氣體非平衡效應(yīng)[1]，這種效應(yīng)會對飛行器的氣動特性產(chǎn)生很大影響。目前對于高溫氣體非平衡效應(yīng)的研究大多結(jié)合試驗進行，由于飛行試驗成本較高，目前試驗研究大多是在地面高焓試驗設(shè)備中進行的。近幾十年來高焓試驗設(shè)備的研制獲得了高度的重視，世界各地發(fā)展出了各種不同的高焓試驗設(shè)備[2-4]。在眾多種類的設(shè)備中，高焓激波風洞可以模擬的總溫總壓高，運行成本低，是目前國際高焓流動試驗研究應(yīng)用的主力試驗手段[5]。

高焓激波風洞通過強激波壓縮得到高溫高壓狀態(tài)的試驗氣體，然后通過噴管的膨脹加速得到超高速試驗氣流。形成強激波一般有多種方式，常見的有加熱輕氣體、燃燒加熱、變截面、爆轟驅(qū)動等[6]。20世紀90年代，俞鴻儒[7]等提出了通過爆轟驅(qū)動獲得強激波和高焓試驗氣流的理論，在此基礎(chǔ)上設(shè)計的JF-10氫氧爆轟驅(qū)動激波風洞取得了巨大成功。風洞利用驅(qū)動段爆轟燃燒產(chǎn)生的高溫高壓氣體作為驅(qū)動氣體，為國際首創(chuàng)，具有有效試驗時間長、耗費低、運行平穩(wěn)、復現(xiàn)性好的特點。風洞可以獲得超高速高焓試驗氣流，具有模擬高溫氣體非平衡效應(yīng)的能力。

駐室中的試驗氣體可認為是處于熱化學平衡態(tài)的多組元混合物。但隨著氣體在噴管中的不斷加速，特別是經(jīng)過喉道后的迅速膨脹，流動來不及達到當?shù)販囟取毫ο碌臒峄瘜W平衡，喉道下游還會出現(xiàn)明顯的熱化學凍結(jié)現(xiàn)象[8-9]。在這種情況下，高焓風洞試驗數(shù)據(jù)的分析和使用，與常規(guī)風洞試驗相比，要困難和復雜得多，一般需要借助數(shù)值計算的方法進行研究。

在現(xiàn)有的高焓風洞數(shù)值研究中，常常將噴管流動和模型試驗解耦處理，即先通過試驗測量或者數(shù)值模擬等方式得到風洞噴管出口的參數(shù)(或者略作處理，如平均[10-11]、外延[12]等)，然后將其作為試驗?zāi)Ｐ偷木鶆騺砹鏖_展研究。通過這種解耦處理，可以相對獨立地進行噴管流動與模型試驗研究，因而應(yīng)用十分廣泛。曾明等在2006年[10-11]對JF-10激波風洞流場的參數(shù)測量進行了數(shù)值重建；2013年，Clemente等[12]對意大利SCIROCCO等離子風洞中的EXPERT太空艙試驗進行了數(shù)值模擬；Ishihara等在2013年[13]和2014年[14]對日本HIEST自由活塞激波風洞中的鈍頭體試驗進行了數(shù)值模擬；2015年，Holden等[15]綜述了在LENS系列激波風洞中進行的試驗及其數(shù)值模擬；2016年，Vasilevskii等[16]對俄羅斯VAT-104風洞的鈍錐體試驗進行了數(shù)值模擬。盡管解耦方法應(yīng)用較多，但它在試驗?zāi)Ｐ土鲃友芯恐泻雎粤藝姽艹隹诤驮囼灦瘟鲌龊诵膮^(qū)參數(shù)的非均勻性，不能較為真實地反應(yīng)高焓風洞試驗的過程。

本文以更加真實地反映風洞試驗物理過程為目的，基于一體化數(shù)值模擬思路，考慮噴管出口和試驗段流場的非均勻性，開展高焓激波風洞JF-10典型運行狀態(tài)下噴管/試驗段/模型流場的一體化數(shù)值模擬，分析噴管出口和試驗段流場非均勻性的產(chǎn)生原因及其對試驗?zāi)Ｐ土鲌鎏匦缘挠绊懀⑴c兩種噴管/試驗段模型流場解耦求解的方法進行對比，分析多個模型位置、攻角試驗條件下不同計算方法對流場參數(shù)分布、試驗?zāi)Ｐ捅砻鏌崃鞣植肌⒛Ｐ蜌鈩恿ο禂?shù)的影響。

1 計算方法及物理化學模型

1.1 控制方程及求解方法

風洞流場包含明顯的熱化學凍結(jié)現(xiàn)象，因此流場控制方程采用包含熱化學非平衡源項的三維Navier-Stokes方程，其無量綱形式為：

(1)

Q=(ρi,ρEV,ρ,ρu,ρv,ρw,ρE)T

W=(wi,wV,0,0,0,0,0)T

式中Q是守恒變量，Re為雷諾數(shù)，F(xiàn)、G、H和FV、GV、HV分別對應(yīng)x、y、z三個方向的對流項和黏性項，W為熱化學非平衡源項。計算采用結(jié)構(gòu)網(wǎng)格，控制方程離散采用有限差分方法，對流項通過AUSMPW+(Advection Upstream Splitting Method by Pressure-based Weight functions)格式離散，黏性項通過中心差分格式離散，時間項采用LU-SGS(Lower-Upper Symmetric Gauss Seidel)格式。非平衡源項、對流項和黏性項均采用全隱式處理。具體計算公式和方法詳見文獻[17-18]。

1.2 熱化學反應(yīng)模型

由于風洞試驗氣體為空氣，因此本文選用7組分(O2、N2、NO、O、N、NO+、e)Dunn-Kang化學反應(yīng)模型[18]；熱力學模型采用Park兩溫度模型[19]，考慮振動非平衡效應(yīng)。具體計算公式和方法詳見文獻[18]。

1.3 表面催化特性模型

2 計算方法校驗

針對JF-10激波風洞，開展了噴管/試驗段一體化數(shù)值模擬研究，并與試驗及文獻結(jié)果進行對比分析。JF-10激波風洞的噴管為錐形噴管，半錐角為7.1°，出口直徑為500 mm，喉道直徑為11 mm，總長為2 m。試驗段長2 m，內(nèi)徑1200 mm，計算幾何示意圖如圖1所示。駐室總壓為19.6 MPa，總溫為7920 K，試驗氣體為空氣。壁面條件采用等溫壁Tw=300 K條件、完全催化條件(FCW)，入口條件采用亞聲速入流邊界條件。

圖1 計算幾何示意圖

圖2(a)和圖2(b)給出了沿噴管軸線的溫度和部分組分質(zhì)量分數(shù)分布與文獻[10]結(jié)果的對比，橫坐標為0處為噴管喉道，文獻采用的是二維軸對稱計算，而本文為三維計算。

(a)溫度分布

從圖2中可以看出，本文計算得到的平動溫度、振動溫度、氧氣、氮氣等典型參數(shù)分布與文獻[10]符合良好。同時可以看出，在噴管出口處，溫度(平動-轉(zhuǎn)動溫度)約450 K，振動溫度在3500 K以上，存在熱力學振動溫度的凍結(jié)現(xiàn)象；溫度450 K的平衡態(tài)空氣中不會有O、N、NO等組分[10,18]，而在噴管出口處，氣流中仍存在一定濃度的O、N、NO等組分，可見，由于流速較高，氣流出現(xiàn)化學反應(yīng)的凍結(jié)現(xiàn)象。

圖3給出了激波風洞中皮托管的數(shù)值重建結(jié)果與試驗測量結(jié)果對比。皮托壓力傳感器安裝在半徑Rn=0.0095 m的球頭圓柱體頂點，通過計算風洞中三維球頭繞流流場，開展皮托壓力測量的數(shù)值重建。以噴管/試驗段一體化數(shù)值模擬中得到的皮托管安置處的自由來流參數(shù)作為來流條件，計算球頭熱化學非平衡流場，采用等溫壁、完全催化壁(FCW)條件。圖3(a)給出了d=100 mm處數(shù)值重建的壓力分布云圖，圖3(b)給出了風洞軸線上皮托壓力測量的數(shù)值重建結(jié)果與風洞試驗結(jié)果的對比，其中橫坐標d表示皮托管距噴管出口的距離。可以看出，數(shù)值重建結(jié)果與試驗結(jié)果符合良好。

(a)球頭壓力分布云圖

總的來說，計算結(jié)果與試驗及文獻結(jié)果吻合良好，計算方法和程序可信。

3 計算結(jié)果及分析

針對JF-10高焓激波風洞，開展了噴管/試驗段/試驗?zāi)Ｐ鸵惑w化數(shù)值模擬。計算域及計算網(wǎng)格如圖4所示，試驗?zāi)Ｐ蜑殁g錐體，模型總長為224.4 mm，底部半徑為75 mm，頭部半徑為50 mm，半錐角為6.9°，圖中d表示試驗?zāi)Ｐ途鄧姽艹隹诰嚯x。駐室總壓19.6 MPa，總溫7920 K，試驗氣體為空氣，試驗?zāi)Ｐ捅诿鏈囟炔捎玫葴乇赥w=300 K條件，表面催化特性采用完全催化(FCW)和非催化(NCW)條件。

(a)計算域示意圖

圖5給出了未放置模型時風洞流場典型參數(shù)分布，圖5(a)、(b)分別為沿流向溫度和速度、馬赫數(shù)分布，其中橫坐標x<2 m區(qū)域為噴管，x>2 m區(qū)域為試驗段。由圖可以看出，在試驗段，高超聲速氣流沿流向存在一定的非均勻性：沿流動方向，平動溫度有一定下降，馬赫數(shù)緩慢升高。事實上，這種沿流向的非均勻性在試驗中也得到了驗證，從圖3(b)的皮托壓力試驗測量可以看出，皮托壓力從噴管出口(d=0 mm)處的21kPa一直下降到距噴管出口d=600 mm處的12.4 kPa。

圖5(c)給出了本文計算得到的噴管馬赫數(shù)分布云圖，由圖可以看出，在噴管出口附近，存在y方向的非均勻性，噴管軸線附近區(qū)域的馬赫數(shù)略高于遠離軸線區(qū)域。圖5(d)為試驗段中距噴管出口d=200 mm處沿徑向靜壓分布，其中縱坐標y表示距風洞軸線的距離。由圖可以看出，在試驗段核心區(qū)域，氣流也存在y方向的非均勻性，中心軸線區(qū)域(y=0 m)靜壓較低，在兩端界面附近(|y|=0.22 m)出現(xiàn)局部峰值，最大相差12.6%。

試驗段流向(x方向)和徑向(y方向)的非均勻性產(chǎn)生原因，可結(jié)合圖5和圖6從兩個方面進行分析：(1)由于噴管擴張段存在一定的擴張角，氣流膨脹擴張，噴管出口處氣流存在非x方向速度，從圖6(a)給出的噴管出口處無量綱徑向速度(通過軸向速度Vx無量綱化)沿徑向分布以及圖6(b)給出的噴管徑向速度分布云圖可以看出，x-y平面，在y=0.2 m附近，氣流的y方向速度可達到500 m/s左右，約為軸向速度的10%，這種非x方向速度，會使氣流在試驗段呈擴張趨勢；(2)噴管壁面不可避免的存在黏性效應(yīng)，在噴管擴張段，沿x正方向，馬赫數(shù)逐漸增大，黏性效應(yīng)逐漸明顯，壁面附近溫度由于摩擦顯著升高，邊界層不斷增厚，形成高溫區(qū)域，使y方向的壓力、溫度、馬赫數(shù)等參數(shù)分布不均勻，這種非均勻性會隨著流動向下游擴展，影響區(qū)域逐漸增大。

(a)溫度與速度分布

圖6 噴管中垂直速度分量分布

由于這種流動的非均勻性，試驗?zāi)Ｐ吞幱谠囼灦蔚牟煌瑓^(qū)域時，其來流條件和流場流動特性將會發(fā)生變化，會對其氣動力/熱特性的預(yù)測結(jié)果產(chǎn)生一定影響。而傳統(tǒng)的解耦方法，忽略了流動的非均勻性，有可能帶來預(yù)測誤差。下面將針對這一現(xiàn)象，具體展開分析。

3.1 不同模型位置下的差別

在實際試驗中，由于試驗?zāi)Ｐ统叽缁蛘邷y量點設(shè)置等因素，模型在試驗段中的安放位置有多種可能。這里我們將試驗?zāi)Ｐ途鄧姽艹隹诰嚯xd分別取為0 mm、100 mm、200 mm和300 mm。考慮試驗段流動的非均勻性，將噴管/試驗段/試驗?zāi)Ｐ鸵惑w化數(shù)值模擬，記為Coupled。傳統(tǒng)的解耦方法，取噴管出口中心處參數(shù)，作為試驗?zāi)Ｐ偷木鶆騺砹鲄?shù)進行數(shù)值模擬，記為Case1，該方法忽略了試驗段噴流流場流向和徑向的非均勻性。另一種解耦方法，取試驗?zāi)Ｐ退幬恢玫牟ㄇ皩ΨQ軸線處參數(shù)，作為試驗?zāi)Ｐ途鶆騺砹鲄?shù)進行數(shù)值模擬，記為Case2，這種解耦方法忽略了試驗段噴流流場徑向的非均勻性。

為了排除網(wǎng)格差異的影響，Case1、Case2方法的計算網(wǎng)格與一體化模擬Coupled方法中,試驗?zāi)Ｐ筒糠值木W(wǎng)格保持一致。同時還開展了兩套網(wǎng)格的計算研究：采用試驗?zāi)Ｐ筒糠至飨颉⒎ㄏ颉⒅芟蚓W(wǎng)格點數(shù)分別為141×61×31和141×81×41，最小網(wǎng)格間距分別為5×10-6m和2×10-6m的網(wǎng)格進行一體化計算，記為Grid1和Grid2，其熱流分布如圖7所示，可以看出，兩套網(wǎng)格的熱流值變化不超過2%。計算結(jié)果受網(wǎng)格影響較小，因此本文計算中如無特殊聲明，均采用Grid2網(wǎng)格。

圖7 不同計算網(wǎng)格的模型表面熱流分布對比

圖8給出了模型距噴管距離d=200 mm時放置模型前后的馬赫數(shù)分布云圖。可以看出，置入試驗?zāi)Ｐ停谠囼災(zāi)Ｐ蜕嫌螀^(qū)域，流場變化幅度較小；而在試驗?zāi)Ｐ拖掠危渭げㄅc試驗段超聲速氣流外邊界存在強烈的相互干擾，氣流膨脹擴張，流動結(jié)構(gòu)十分復雜。

圖8 馬赫數(shù)分布云圖(d=200 mm)

圖9和圖10分別給出了d=200 mm時模型駐點線上的溫度和部分組分質(zhì)量分數(shù)分布。可以看出，在大體上，一體化計算(Coupled)得到的流場參數(shù)分布，與解耦方法計算(Case1、Case2)結(jié)果基本一致，但在局部細節(jié)上存在差異：解耦方法得到的波后溫度、分子組分質(zhì)量分數(shù)較一體化方法稍低，而原子組分則稍高；解耦方法Case2較Case1，計算結(jié)果更接近一體化計算；解耦方法Case2和Case1計算結(jié)果差異較小，可以看出，這種條件下，均勻來流的絕對值細微區(qū)別(即流向非均勻性)對于駐點線參數(shù)分布影響不大。

圖9 模型駐點線上的溫度分布(d=200 mm, FCW)

圖10 模型駐點線上的部分質(zhì)量組分分布(d=200 mm, FCW)

圖11給出了不同模型位置頭部駐點熱流的對比。由圖可以看出：采用一體化計算(Coupled)得到的頭部駐點熱流隨著距離d的增加，逐漸下降，在完全催化條件下(FCW)，其下降幅度更顯著；與一體化計算相比較，采用傳統(tǒng)的解耦方法Case1計算得到頭部駐點熱流隨著距離d的增加，偏差逐漸增大，在d=300 mm時偏差可達10%左右；采用解耦方法Case2，其計算結(jié)果也與一體化計算也存在一定差異，但同Case1方法相比，其偏差幅度整體上小一些；解耦方法Case2和一體化計算Coupled差別小于解耦方法Case1與Case2差別，這說明對于頭部駐點區(qū)域而言，噴流流向非均勻影響大于徑向非均勻性影響。

圖11 不同位置下的駐點熱流對比

圖12進一步給出了d=200 mm時模型表面熱流分布。可以看出：沿流動方向，一體化計算與解耦計算的熱流差異呈增大趨勢，在身部區(qū)域熱流值最大可相差20%；Case2較Case1更接近Coupled的結(jié)果；在靠近底部的錐身區(qū)域，解耦方法Case2和一體化Coupled的計算差別，大于解耦方法Case1與Case2的差別，造成這種現(xiàn)象的原因，是由于在身部區(qū)域，解耦方法Case1與Case2的差別主要由流向非均勻性引起。而一體化計算Coupled與解耦方法Case2的差異，既包含了徑向非均勻性的影響，同時也包含了由于模型身部在流向位置上的差異引起的影響。

圖12 模型表面熱流分布(d=200 mm)

圖13給出了不同模型位置時試驗?zāi)Ｐ偷淖枇ο禂?shù)CD。可以看出：一體化計算結(jié)果(Coupled)與解耦計算(Case1、Case2)存在一定差異，解耦計算得到的阻力系數(shù)偏大，最大可達5%左右；隨著距離d增大，一體化計算Coupled得到的阻力系數(shù)呈增大趨勢，解耦方法Case1與一體化計算結(jié)果差異減小，而解耦方法Case2與一體化的差異基本保持不變；解耦方法Case2和一體化Coupled之間的計算差別，大于解耦方法Case1與Case2之間的差別，造成這種現(xiàn)象的原因，是由于流向非均勻性導致的Case1與Case2的差異，只是這種來流參數(shù)絕對值的差異對于無量綱的氣動系數(shù)影響較小，而一體化計算中考慮了流動的擴張，擴張導致的徑向非均勻性影響了流場結(jié)構(gòu)、壓力分布，因此造成了更大的氣動力系數(shù)差異。

圖13 不同位置下的阻力系數(shù)對比

模型表面壓力的差別，是氣動力系數(shù)計算差異的主要來源之一。圖14(a)給出了本文計算的模型表面的壓力分布(d=200 mm)。圖14(b)給出了文獻[12]的計算結(jié)果(文獻的計算條件和外形，與本文不一致)，圖中紅線為解耦計算結(jié)果，藍線考慮了試驗段流場的非均勻性。可以看出，考慮試驗段流場的非均勻性時，計算結(jié)果與試驗結(jié)果符合更好。對比圖14(a)和圖14(b)，可以看出，盡管本文計算條件和外形與文獻不一致，但其規(guī)律變化相同：采用解耦方法，計算得到的模型表面壓力偏高，沿流動方向，其差距呈增大趨勢。由此可見，要準確模擬高焓風洞試驗段的非平衡流動，得到較為精確的氣動力/熱特性，必須考慮試驗段流場的非均勻性，采用一體化或類似的數(shù)值模擬方法。

(a)d=200 mm

3.2 不同模型攻角下的差別

保持試驗?zāi)Ｐ臀恢貌蛔儯Ｐ途鄧姽艹隹诰嚯xd=200 mm，試驗?zāi)Ｐ惋w行攻角分別取為0°、5°、10°、15°。

圖15給出了不同飛行攻角時模型頭部駐點熱流值和迎風面肩部頂點的熱流值對比，圖中Coupled、Case1和Case2含義與3.1節(jié)相同。可以看出，隨攻角增大，一體化計算(Coupled)的頭部駐點熱流小幅下降，在15°攻角時下降了5%左右，這是由于徑向的非均勻性造成的，而解耦計算(Case1、Case2)的熱流值則基本保持不變；三者的迎風面肩部熱流呈上升趨勢；解耦方法得到的熱流要高于一體化計算結(jié)果，二者的差距隨著攻角增大而增大，在15°攻角時Case1與Coupled差別可達23%；Case2熱流值較Case1更接近于Coupled。

(a)駐點熱流對比

圖16進一步給出了FCW條件、10°攻角下的試驗?zāi)Ｐ捅砻鏌崃鞣植迹瑘D中Windward表示迎風面，Leeward表示背風面。可以看出，模型攻角的變化使模型表面熱流分布發(fā)生了明顯變化，隨著攻角增大，迎風一側(cè)的熱流升高；在頭部區(qū)域，一體化計算熱流值在Case1與Case2之間，更接近Case2，而Case1與Case2差異明顯，這說明該區(qū)域流向非均勻性影響較大；而在靠近底部的錐身區(qū)域，解耦計算熱流結(jié)果要高于一體化計算結(jié)果，Case1與Case2的差異小于Coupled與Case2的差異，這與圖11和圖12現(xiàn)象一致。

圖16 模型表面熱流分布(α=10°)

圖17給出了不同攻角下得到的氣動力系數(shù)對比，其中CL為升力系數(shù)，CMZ為俯仰力矩，其中力矩參考點為頭部駐點。可以看出：解耦計算方法和一體化計算方法結(jié)果存在一定差異，隨攻角增大，差異呈增大趨勢；解耦方法Case1與Case2差別相對較小，這說明對于該條件下的氣動力特性而言，噴流徑向非均勻影響大于流向非均勻性影響，這與圖13結(jié)論一致。

圖17 不同攻角下的氣動力系數(shù)對比

4 結(jié) 論

本文考慮風洞試驗中的高溫空氣化學反應(yīng)，振動激發(fā)及非平衡效應(yīng)，通過數(shù)值求解三維熱化學非平衡Navier-Stokes方程，完善了高焓激波風洞流場一體化數(shù)值模擬計算方法和計算程序。考核校驗表明：數(shù)值計算結(jié)果與文獻及試驗結(jié)果符合良好，數(shù)值模擬結(jié)果可信度高。

1)通過計算發(fā)現(xiàn)，試驗段流場在流向和徑向存在一定的非均勻性。試驗?zāi)Ｐ吞幱谠囼灦蔚牟煌瑓^(qū)域時，其流場流動特性將會發(fā)生一定變化，對試驗?zāi)Ｐ蜌鈩恿μ匦浴鈩訜岘h(huán)境的預(yù)測結(jié)果產(chǎn)生影響。在進行數(shù)值模擬和試驗方案設(shè)計時有必要考慮試驗段流場的非均勻性。

2)通過不同解耦方法(Case1、Case2)與一體化數(shù)值模擬結(jié)果對比可以看出：流場非均勻性造成的氣動力熱差異，隨著模型與噴管出口距離d增大，呈增大的趨勢；對試驗?zāi)Ｐ筒煌瑓^(qū)域的影響存在差別，在駐點區(qū)域，流向非均勻性影響較大，在身部區(qū)域，徑向的非均勻性影響較大。

3)解耦計算(Case1、Case2)與一體化計算的差異在不同攻角下的差別，主要體現(xiàn)了流場徑向非均勻性的影響。隨飛行攻角增大，差異整體上呈增大趨勢；三者在身部區(qū)域差別較大，駐點附近區(qū)域相對較小，在迎風面的肩部頂點，最大差別可達20%以上。

從本文的結(jié)果來看，高焓風洞試驗段流場的非均勻性對于模型氣動力/熱的影響，與試驗?zāi)Ｐ偷陌惭b位置、攻角、外形相關(guān)，為了得到更為準確的試驗數(shù)據(jù)，需要結(jié)合具體試驗的一體化數(shù)值模擬的結(jié)果，對試驗方案以及試驗數(shù)據(jù)進行修正。基于一體化數(shù)值模擬的思想，數(shù)值計算可以更加深入地理解高焓風洞模擬特點，改善風洞試驗數(shù)值模擬精度。

本文的數(shù)值模擬方法還可以應(yīng)用于JF-12激波風洞[20]等其他高焓設(shè)備的模擬研究。