跨聲速風洞高速段一體化數值模擬研究

楊小川，王運濤,洪俊武,黃知龍，孫運強，江雄

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

(2.中國空氣動力研究與發展中心 設備設計及測試技術研究所，綿陽 621000)

0 引 言

隨著航空航天飛行器逐漸步入精細化設計的時代，特別是追求高亞聲速經濟巡航的民航客機以及跨聲速高機動特性的戰斗機，高性能跨聲速風洞的需求日趨緊迫[1]。

目前，全世界有兩座高性能跨聲速風洞，分別是美國NTF(1983年第一次運行)和歐洲ETW(1993年第一次運行)，投資分別為同時期的8 500萬美元和5.62億馬克，由于ETW較NTF晚十年建設，在某些指標上更為先進。關于建設ETW跨聲速風洞的起因，其推動力來自于1968年9月AGARD流體力學專家座談會，會上報道了發生在19世紀60年代出現的眾多跨聲速空氣動力學問題，如跨聲速飛行面臨的阻力劇增及焦點后移等氣動問題[2]。同時，跨聲速風洞能模擬真實飛行雷諾數，得到更為精準的試驗數據供設計參考，如C-141飛機在超臨界狀態下，超臨界機翼升阻和力矩特性對雷諾數特別敏感，真實飛行與風洞試驗結果在機翼中部站位的激波位置相差近20%弦長[3]，而超臨界翼型和機翼設計技術又是現代跨聲速運輸類飛機氣動力設計中最重要的關鍵技術之一[4]，因此高性能跨聲速風洞對現代飛行器設計意義重大。

對于跨聲速狀態而言，風洞試驗段模型主要面臨三個問題：①激波/膨脹波壁面反射；②氣流擁塞；③相對自由來流的氣流變化。因此跨聲速風洞設計的關鍵在于消除/抑制上述物理問題，進而發展了槽壁或孔壁等技術，以及駐室和抽吸氣等輔助手段[5]。關于跨聲速狀態風洞試驗，直到1945年還沒有可行的方案使試驗段馬赫數大于0.9。隨后，Wright發展了一種縱向槽壁理論模型，用于消除試驗段壁板和開口區域邊界層干擾，并于1946年在安裝有槽道試驗段的風洞中獲得成功，試驗段馬赫數達到跨聲速范圍，該成果很快應用到蘭利的8和16 ft(1 ft=0.304 8 m)高速風洞中[2]。

在風洞數值模擬研究方面，國外，M.Kohzai等[6]采用滲透邊界對跨聲速風洞試驗段進行支撐系統干擾模擬，并與試驗結果進行對比；B.Goffert等[5]采用三維Euler方程對跨聲速風洞中的翼型進行壁板有無開槽影響研究，并與試驗結果中翼型壓力系數分布、PSP分布及三維效應進行對比，其中開槽計算結果與試驗吻合較好；S.A.Glazkov等[7]采用98萬網格對簡化的單槽道試驗段模型進行三維粘性模擬，模型包括擴張角0.55°的槽壁、7°外開角的再導入調節片(或稱指片Finger Flap)以及駐室，并與T125風洞試驗進行對比分析，發現試驗中存在明顯的橫向流動，而采用單槽道試驗段計算結果無橫向流動；A.R.Gorbushin等[8]首次采用具有低溫求解模塊的EWT-TsAGI軟件對ETW風洞進行了帶CRM模型的高速段一體化模擬(計算網格包括風洞試驗段、槽壁、再導入調節片、駐室、彎刀尾撐機構以及CRM模型等部件，半模網格量達到1.7億)，該軟件考慮了低溫增壓狀態下氮氣的真實氣體效應，計算得到的試驗段下壁面壓力分布與試驗值吻合較好。

國內，范潔川[9]對低速風洞洞壁干擾進行了大量研究，并得到合理的開閉比和開槽形式等參數；叢成華等[10]采用三維數值模擬方法對不同抽氣量、抽氣口位置以及壁板厚度、駐室容積、開閉比及槽型等參數進行對比分析，得到不同條件對跨聲速試驗段流動的影響情況；Qu K等[11]采用雷諾平均N-S方程對二維翼型進行有無開孔情況下的壁面干擾分析，并對不同孔壁高度和開孔率進行研究；劉光遠等[12]對不同孔壁開孔率、引射縫、擴開角等狀態下試驗段核心流均勻性進行研究，得到了大型跨聲速風洞試驗段壁板參數影響的主要規律；江雄等[13]采用基于ARK方法模擬氮氣低溫高壓真實氣體效應的RANS求解軟件，對典型運輸機構型DLR-F6進行低溫真實氣體效應影響研究，結果顯示真實氣體效應引起的氣動力差異相對雷諾數效應很小。上述研究對于提高我國跨聲速風洞試驗段設計能力起到推動作用，但關于跨聲速風洞高速段一體化的復雜問題的研究仍相對較少。

本文采用中國空氣動力研究與發展中心自主研發的“亞跨超CFD軟件平臺”(TRIP3.0)[14-15]，對某跨聲速風洞回路高速段進行一體化數值模擬初步研究，控制方程為雷諾平均N-S方程，采用有限體積法離散控制方程，空間方向無粘項離散采用MUSCL-Roe格式，湍流模型運用SST湍流模型。初步得到數值計算收斂評判方法以及不同初始條件和槽道擴張角對試驗段流場品質的影響，并對槽道流動進行分析。

1 計算方法

1.1 控制方程

假設采用慣性笛卡兒坐標系，忽略徹體力，則Euler/Navier-Stokes方程為

(1)

其中，

式中：ρ,u,v，w,p,e和h分別為密度、x，y，z方向的絕對速度分量、壓力、內能和總焓。

當NVIS=0時，方程形式為Euler方程；當NVIS=1時，方程形式為N-S方程。

1.2 拼接網格方法

拼接網格也稱滑移網格(Sliding Mesh)，多用于“剪刀縫”或構型復雜的計算外形以及存在多體相對運動的物體，與重疊網格功能類似。其基本特點是：存在拼接的網格塊之間相互獨立，在交接面上網格可不對應，即不用考慮相鄰塊的拓撲關系，僅需要網格塊之間存在交接面。為了保證計算的收斂性，交界面兩側網格尺度及分布應盡可能接近。

拼接采用面積加權的線性插值方式，保證兩側通量守恒。該計算策略降低了網格的生成難度，特別是結構網格生成。另外，拼接面兩側網格塊不需要進行對應，對單側網格塊進行加密，另一側網格塊不受影響，使計算網格規模適當降低[16]。

1.3 邊界條件

在邊界條件方面，物面上采用無滑移條件，入口邊界給定總壓P0，取穩定段總壓；出口邊界給定背壓Pb，且存在

σ=Pb/PIN

(2)

式中：PIN為入口靜壓或方程無量綱化靜壓；σ為出口壓比。

根據試驗馬赫數，不斷更改擴壓段出口背壓，使試驗段實際馬赫數與試驗馬赫數相近。

2 算例驗證

采用非對稱平板擴壓器(Plane Asymmetric Diffuser，簡稱PAD)算例，進行內流場中分離與再附等粘性流問題驗證，主要對計算湍流模型進行考核。關于非對稱平板擴壓器有較多試驗，而擴張角由10°減小到8.5°，能減小分離區域，獲得更典型的流動二維性以及穩定的分離和再附點，即非定常效應影響更小[17]，故選取8.5°擴張角非對稱平板擴壓器作為研究對象。

算例中擴壓器共分為三段：入口段、擴壓段以及出口段，擴壓段擴張角為8.5°。該擴壓器入口來流速度約為19.48 m/s，雷諾數為40 000。入口段、擴壓段以及出口段流向的網格分布分別為121、165和137個點，展向和法向網格分布為69×81，采用全對接結構網格，網格量約為230萬。

擴壓段截面速度分布與試驗結果對比如圖1所示，可以看出：計算值與文獻[17]中試驗值吻合較好，在x/H為15～25截面靠近傾斜壁板一側區域存在差異，主要是該區域存在較大范圍分離所致。

擴壓段截面速度云圖及流線如圖2所示，分離區域如圖3所示，可以看出：計算得到的分離點在x/H=9附近，再附點在x/H=31附近，這與試驗結果相吻合，初步驗證了計算方法的可行性；同時計算得到的分離區域較試驗值均偏小，這可能和網格尺度、湍流模型等因素有關。

圖1 擴壓段截面速度分布

圖2 擴壓段截面速度云圖及流線

圖3 擴壓段分離區域示意圖

3 計算對象及網格

3.1 計算對象

計算對象為某跨聲速風洞回路高速段外形，包括穩定段、收縮段、噴管、試驗段、二喉道、擴張段以及包圍試驗段的駐室。試驗段采用槽壁構型，上下各6根條型槽壁，左右為實壁，包含彎刀、尾撐以及再導入調節片。計算中再導入調節片偏角0°、二喉道開角為0°，且駐室及擴壓段無抽吸氣處理，同時尾撐頭部為圓錐體。

3.2 計算網格

在實際計算中需對真實風洞進行適當簡化，主要在彎刀和駐室兩個部分進行簡化。將伸入駐室的部分彎刀適當扣除，同時將駐室簡化為軸對稱外形，且駐室前段空間適當減小。計算主要針對跨聲速狀態，并對不同時間步長、初始化條件和槽壁擴張角進行分析。

計算網格為結構網格，由于模型可視為上下和左右對稱，因此選取四分之一模型網格作為計算域。為了保證計算精度，特別是計算槽道粘性作用，需在壁面和槽道生成附面層網格，而風洞流道內存在圓方相互轉換以及槽壁臺階等因素影響，采用全對接結構網格生成難度大，故選取拼接網格方法進行處理。

計算網格由三個獨立的全對接結構網格組成，網格拼接處理示意圖如圖4所示，其中三個網格區域共進行兩次拼接處理，分別位于噴管與試驗段結合處和試驗段與二喉道結合處。圖中Cross Field#2為試驗段及對應駐室部分網格，網格量約為5 200萬；Cross Field#1為穩定段、收縮段、噴管段及對應駐室部分網格，網格量約為470萬；Cross Field#3為二喉道、擴壓段及對應駐室部分網格，網格量約為390萬。整個計算域網格量約為6 060萬，由于槽道流動是跨聲速風洞氣動設計的關鍵之處，試驗段網格在槽道拐折等變化較大處均進行適當加密，試驗段網格示意圖如圖5所示。

圖5 跨聲速風洞構型網格示意圖

4 結果分析

對管道流動而言，收斂依據較難判斷，通常采用出口流量或流場參數變化等作為收斂參考，但對于風洞而言，試驗段流場均勻性是評判風洞品質的核心參數，因此本文選取試驗段監測點作為收斂判斷依據。風洞對稱面示意圖如圖6所示，在圖中模型區域軸向位置近似選取一前一后兩個固定點，即監測點Point 1和監測點Point 2，計算中實時輸出兩個監測點的馬赫數，通過兩點馬赫數變化量以及兩者馬赫數差量變化作為判斷依據，X1～X3為流向x方向截面站位位置分布。

圖6 對稱面計算監測點及截面站位位置示意圖

風洞試驗段X站位截面及槽道位置示意圖如圖7所示，其中虛線區域為模型區域，中心點為中軸線處點，60%點為模型區域對角線上60%處點，用于表示中心線或60%站位處軸向馬赫數分布情況。

圖7 試驗段橫截面參考點位置及槽道位置示意圖

4.1 收斂步數

為了保證計算結果接近收斂解，選取兩個典型狀態作為參考，分別為亞聲速狀態和跨聲速狀態。通過分析模型區前后兩個監測點馬赫數的變化以及兩者馬赫數差量變化，來判斷計算收斂情況。

4.1.1 亞聲速狀態

計算狀態：試驗段Ma=0.58，再入調節片偏角0°，上下槽壁，槽壁擴張角0.0°，計算初始化條件為Ma=0.0，為了盡可能保證計算的收斂性，定常計算迭代步數達到200萬步。

監測點馬赫數隨時間步變化曲線如圖8所示，模型區馬赫數從速度為0.0開始逐漸加速到Ma=0.58附近，從圖8(a)可以看出：在10萬步左右監測點馬赫數基本達到穩定值，在10萬～200萬步的計算歷程中，馬赫數均在小幅波動；監測點馬赫數差量在10萬前逐漸縮小到0.005 5附近，且在10萬～200萬步的計算歷程中微幅波動；模型區域監測點1較監測點2馬赫數大，僅為馬赫數差量約0.005 5。從圖8(b)可以看出：監測點1穩定在Ma=0.583 3附近波動，波動幅度達±0.36%；監測點2穩定在Ma=0.578 0附近波動，波動幅度達±0.35%；兩個監測點馬赫數差量穩定在0.005 5附近，波動基本維持在0.005 0～0.006 0之間，整個模型區核心流馬赫數分布均方根偏差為0.001 84。

(a) 計算歷程

(b) 馬赫數波動

4.1.2 跨聲速狀態

計算狀態：試驗段Ma=0.94，再入調節片偏角0°，上下槽壁，槽壁擴張角0.3°，計算初始化條件為Ma=0.9，為了保證計算收斂性，定常計算迭代步數達到60萬步。

監測點馬赫數隨時間步變化曲線如圖9所示，模型區馬赫數從速度Ma=0.9開始逐漸加速超聲速后又回到Ma=0.94附近，從圖9(a)可以看出：在12萬步左右監測點馬赫數差量基本達到穩定值，監測點馬赫數差量逐漸縮小并穩定在0.008 5附近，且12萬～60萬步的計算歷程中微幅波動；模型區域監測點1較監測點2馬赫數大，馬赫數差量約為0.008 5。從圖9(b)可以看出：兩個監測點在12萬步后馬赫數基本穩定；監測點1穩定在Ma=0.944 1附近波動，波動幅度達±0.17%，監測點2穩定在Ma=0.935 6附近波動，波動幅度達±0.17%；兩個監測點馬赫數差量穩定在0.008 5附近，波動基本維持在0.006～0.011之間，整個模型區核心流馬赫數分布均方根偏差為0.002 1。

(a) 計算歷程

(b) 馬赫數波動

通過對上述兩種狀態的收斂性分析發現：對于試驗段馬赫數為低亞聲速狀態，計算在10萬步左右基本收斂；對于試驗段馬赫數為跨聲速狀態，計算在12萬步左右基本收斂。為了保證計算收斂性，后續計算迭代步數均選取20～40萬步。

4.2 初始條件

對比狀態：計算域初始化流場中馬赫數賦值分別為Ma=0.0、Ma=0.9，其他條件均相同，即再導入調節片偏角0°，上下槽壁擴張角0.3°，出口背壓相同，計算步數均為23萬步。

在上述兩種不同初始化流場馬赫數條件下，對稱面及Slot 2槽道截面馬赫數云圖、X站位馬赫數云圖以及槽道流動細節對比情況分別如圖10～圖12所示，其中圖12在x軸方向坐標尺度適當壓縮，以便于觀察流線方向。

(a) Ma=0.0

(b) Ma=0.9

(a) Ma=0.0

(b) Ma=0.9

(a) Ma=0.0

(b) Ma=0.9

從圖10～圖12可以看出：①兩種初值條件得到的馬赫數分布，在對稱面、槽道截面以及X站位處基本相同；②氣流從穩定段開始經過聲速噴管逐漸加速到跨聲速，并在試驗段形成范圍較大較穩定的跨聲速區域，經彎刀尾撐形成局部加速區，同時在二喉道附近減速到亞聲速，并繼續在擴壓段減速；③氣流高速區基本維持在試驗段內，駐室和槽道未出現明顯的高速氣流；④試驗段槽道附近出現明顯的橫向流動；⑤槽道內低速氣流進入試驗段，同時氣流從駐室流向槽道和試驗段，然后在尾撐彎刀附近流回槽道，最后經過再入調節片流回二喉道或駐室，該現象可能是駐室無抽吸氣處理，導致駐室壓力過大，將槽道氣流“壓入”試驗段所致。

兩種不同初始化流場馬赫數條件下，兩個監測點馬赫數收斂曲線以及試驗段軸向馬赫數分布情況如圖13所示，其中黑線和紅線分別代表模型區中心線前后兩監測點馬赫數，藍線代表兩者之差，可以看出：當初值為Ma=0.0時，模型區速度由0.0開始逐漸增大并穩定在0.9附近，兩點差量在5萬步基本收斂，且收斂前抖動劇烈，收斂后抖動大幅減小且穩定在0.01內；當初值為Ma=0.9時，模型區馬赫數先由0.9劇增到1.53，然后減小并穩定在0.9附近，兩點差量基本穩定在0.01內且始終保持小幅抖動；從收斂速度來看，初值為Ma=0.0時，計算到5萬步接近收斂，而初值為Ma=0.9時，計算在3萬步左右收斂，即初值Ma=0.9較初值Ma=0.0計算收斂速度更快。

在圖13(c)試驗段軸向馬赫數分布上，兩種初值條件得到的收斂后結果基本相近，初值Ma=0.0得到的模型區馬赫數較初值Ma=0.9變化趨勢一致，且初值Ma=0.0條件得到的馬赫數略高，這可能是因為管道流動，特別是跨聲速槽道流動復雜，分離區域較多，氣流流動對計算因素變化反應敏感所致。模型區的中心線馬赫數差量即為前后兩個監測點差量，60%站位點(如圖7所示)軸向馬赫數分布較中心線分布更為平緩，馬赫數差量較中心線更小，這主要是模型區后方存在尾撐，而尾撐對中心線方向影響較60%站位點更明顯，且越靠近尾撐，干擾越明顯。

(a) Ma=0.0

(b) Ma=0.9

(c) 試驗段軸向馬赫數分布

通過對比分析，得到兩種初值條件對收斂結果總體影響較小(試驗段馬赫數存在一定差異)，特別是各截面流場分布和槽道流動方向上，兩者結果基本相同，同時兩種初值條件均在5萬步前收斂，5萬步到23萬步變化很小，且兩者前后監測點馬赫數差量均在0.008 5附近。差異之處在于：①收斂速度。初值Ma=0.9狀態在3萬步收斂，而初值Ma=0.0狀態在5萬步收斂；②收斂馬赫數。計算收斂后，初值Ma=0.0狀態得到的前方監測點馬赫數約為0.94，而初值Ma=0.0狀態約為0.93；③收斂過程。初值Ma=0.0狀態監測點馬赫數單調增加到0.9附近，而初值Ma=0.9狀態先增大到1.53后減小到0.90附近，且出現明顯超聲速氣流。

4.3 槽道擴張角

對比狀態：槽壁擴張角分別為0.0°和0.3°，其他條件均相同，即再入調節片偏角0°，上下槽壁，出口背壓相同，計算步數均為39萬步。

槽壁兩種擴張角狀態下，對稱面及Slot 2槽道截面馬赫數云圖、X站位馬赫數云圖以及Slot 2槽道流動細節對比情況分別如圖14～圖16所示。從圖14可以看出：兩種擴張角得到的流場結構基本相同，其中擴張角0.3°狀態得到的馬赫數較擴張角0.0°略高。從圖15可以看出：兩種擴張角得到的流場結構存在一定差異，主要在槽道附近氣流速度上，擴張角0.3°狀態中槽道低速氣流進入試驗段，而擴張角0.0°狀態則是試驗段高速氣流進入槽道。從圖16可以看出：當槽壁擴張角為0.0°時，氣流從試驗段和駐室流入槽道內，然后經過再入調節片流回二喉道或駐室，而當槽壁擴張角為0.3°時，氣流從駐室流向槽道和試驗段，然后再流回槽道內，最后經過再入調節片流回二喉道或駐室。這可能是因為槽壁擴張角0.0°狀態下，槽壁未擴張引起試驗段截面尺寸不變，而附面層逐漸增厚，氣流實際通道面積減小，壓力增加，將試驗段氣流“頂進”槽道。

(a) 擴張角0.0°

(b) 擴張角0.3°

(a) 擴張角0.0°

(b) 擴張角0.3°

(a) 擴張角0.0°

(b) 擴張角0.3°

不同槽壁擴張角狀態下兩個監測點馬赫數收斂曲線以及試驗段軸向馬赫數分布情況如圖17所示，可以看出：兩種槽壁擴張角狀態下，前后兩點馬赫數均在5萬步接近收斂，從5萬步到39萬步變化較小。同時，槽壁擴張角0.3°得到的兩點馬赫數均較槽壁擴張角0.0°略大，這與圖14和圖15中各站位馬赫數分布結論一致。從圖17(a)和(b)可以看出：槽壁擴張角0.0°狀態下前后兩點馬赫數差量保持在0.018 5附近，而槽壁擴張角0.3°狀態下則保持在0.008 5附近，兩者馬赫數差量相差0.01，因此槽壁擴張角0.3°較槽壁擴張角0.0°更能獲到均勻的試驗段流場。圖17(a)中槽壁擴張角0.3°狀態中，模型區中心線前后兩點馬赫數差量抖動較為強烈。從圖17(c)可以看出：槽壁擴張角0.3°狀態下的試驗段馬赫數較槽壁擴張角0.0°大，且在模型區域內馬赫數分布更加平緩，特別是在接近尾撐附近，槽壁擴張角0.3°得到的馬赫數較槽壁擴張角0.0°受尾撐干擾更小。

(a) 擴張角0.0°

(b) 擴張角0.3°

(c) 試驗段軸向馬赫數分布

通過對比分析，得到兩種槽壁擴張角對結果影響明顯，特別是試驗段軸向馬赫數分布上，槽壁擴張角0.3°得到的模型區馬赫數差量為0.008 5附近，而槽壁擴張角0.0°狀態到達0.018 5，即槽壁擴張角0.3°狀態得到的試驗段模型區域流場品質更高，且同樣條件下得到的試驗段馬赫更大。風洞高速段槽道處馬赫數云圖如圖18所示，可以看出：氣流經穩定段到收縮段噴管速度逐漸增加，在試驗段區域保持在跨聲速附近，且試驗段對稱面馬赫數分布均勻。

圖18 槽壁擴張角0.3°槽道馬赫數云圖

5 結 論

(1) 采用模型區前后兩個監測點馬赫數變化以及兩者馬赫數差量變化作為數值計算收斂評判依據，方法可行且能得到較穩定的模型區流場。

(2) 不同計算初始化條件對收斂結果總體影響較小(試驗段馬赫數存在一定差異)，特別是各截面流場分布和槽道流動方向上，兩者結果基本相同。

(3) 在跨聲速狀態適當增加槽壁擴張角，可提高試驗段模型區域流場品質。

(4) 槽道內氣流流動方向受槽壁擴張角影響明顯，且與試驗段和駐室壓力關系密切。

由于管道流動較外流復雜，且跨聲速風洞高速段一體化數值模擬可參考的研究工作較少，特別是試驗段槽壁流動情況，關于試驗段槽道附近流動等工作還需深入分析。