國家數值風洞(NNW)工程在高超聲速中的應用研究進展

陳琦,陳堅強,袁先旭,郭啟龍,萬釗,邱波,李辰,張毅鋒,*

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

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

長期以來，中國CFD軟件市場基本被國際商業軟件壟斷，國產CFD軟件在體系化、工程化以及市場化等方面與歐美國家存在較大差距。為此，中國在2018年啟動了國家數值風洞(NNW)工程，通過自主發展算法模型、研發系列CFD軟件、建設高性能計算機系統，形成自主知識產權、國內開放共享的空氣動力數值模擬平臺，滿足中國航空航天等領域對自主CFD軟件產品的需求[1]。

為使國內CFD從業人員了解國家數值風洞工程，促進NNW系列軟件在全國各領域的推廣和應用，本文以高超聲速流動為例，選取了復雜流動特征高保真模擬、復雜干擾機制高清晰認知以及多學科耦合應用3個專題，從不同角度展示NNW系列軟件的應用效果和功能特點。其中，復雜流動特征專題包含微型渦流發生器在激波邊界層干擾中的應用，一種雷諾平均/大渦模擬 (Reynolds-Averaged Navier-Stockes equations/Large Eddy Simulation, RANS/LES)混合模型在圓柱底部繞流中的應用，可壓縮解析壁面函數在熱流高效預測中的應用，以及采用直接數值模擬(DNS)在圓錐邊界層Mack模態的線性和非線性演化過程中的應用等4部分內容。復雜干擾機制方面，揭示了鈍舵縫隙內部渦串干擾機制，分析了吸氣式飛行器進氣道唇口激波振蕩高低頻轉化過程，研究了噴流誘導回流激波產生底部高溫區的現象。在多學科耦合應用方面，研究了舵面操縱下導彈的動態響應特性，計算了氣動伺服彈性效應對導彈控制系統的影響，分析了化學非平衡效應對動態穩定性的影響規律。

1 復雜流動特征高保真模擬

1.1 激波-邊界層干擾精細模擬

激波-邊界層干擾是激波與邊界層相互作用的一種典型的非線性、非定常、黏性主導的復雜流動，是一種普遍存在于高超聲速飛行器的流動現象。由于它可能會帶來總壓降低、流動非定常振蕩、局部高熱流等負面效應，因此對飛行器的氣動特性和飛行安全有重要的影響，是高速飛行器設計中必須考慮的關鍵問題。在數值模擬時，由于流場中同時存在激波與多尺度湍流邊界層，故要求所采用的數值方法在結構分辨能力和計算穩定性之間達到平衡。激波-邊界層干擾問題的直接數值模擬，主要應用了NNW工程發展和集成的高精度WENN/WENO系列算法[1-3]、激波識別技術[4-5]和湍流入口生成技術[6-7]。

圖1給出了計算域內(圖中X、Y、Z分別表示3個方向空間無量綱坐標)超聲速壓縮拐角的瞬時旋渦結構，用流向速度著色，Q=0.01。可以看出，在干擾前后流場中旋渦分布與尺度存在差異：在干擾區前旋渦結構的尺度要比之后的略小，并且更靠近邊界層外層，經過激波后，湍流強度變強，邊界層變厚。圖2給出了物面壓力分布與實驗數據的比較，其中X為流向坐標，δ為入口邊界層厚度，Pw表示壁面壓力，Pinf表示以來流速度和密度計算的參考壓力，計算得到的壓力分布與實驗值整體符合較好。圖2 中，在拐角處出現了一個壓力拐點，在拐點前壓力的增長速度變小，這是分離區的位置。

圖2 物面壓力分布與實驗對比Fig.2 Pressure distribution on wall contrast to experiment data

在高超聲速的激波-邊界層干擾方面也進行了直接數值模擬，獲得了精細的旋渦結構和低頻非定常運動過程；并且還對微型渦流發生器(MVG)在激波邊界層干擾中的應用開展了研究，圖3給出了其平均及瞬時密度分布。計算采用Duan等[9]的高超聲速邊界層DNS來流條件：來流馬赫數Ma∞=5.0，雷諾數Reθin=2 300，來流溫度T∞=228 K，壁面溫度Tw=433 K。網格及邊界的詳細設置參見文獻[10,11]。圖中:MW表示馬赫波，結構“S1”和“S2”分別表示前緣激波和后緣激波，分區“I”和“II”分別表示近穩態區和間歇性大尺度渦流區。

圖3 平均及瞬時密度分布Fig.3 Average and instantaneous density contours

研究發現，微型渦流發生器能夠有效減小分離區長度，并且對于干擾區后的熱流也有顯著的調制作用。進一步研究微型渦流發生器后速度和湍流脈動剖面的自相似特性，發現與超聲速情況類似，高超聲速條件下微型渦流發生器后的速度和湍流脈動剖面存在自相似性，并且給出了自相似律，如圖4和圖5所示[11]，其中散點為不同流向位置的實驗數據，U、V表示流向和法向的平均速度，u′、v′分別表示流向和法向脈動速度的均方根，u′v′表示雷諾切應力，各個物理量的下標“scale”表示對其進行相似變換，具體的變換公式參見文獻[11]。

圖4 歸一化的流向和法向速度分布[11]Fig.4 Distribution of scaled streamwise and normal velocity[11]

圖5 歸一化的u′v′u′v′分布[11]Fig.5 Distribution of scaled turbulent quantities u′v′u′v′[11]

1.2 湍流高精度模擬

Prandtl的混合長度模型與Smagorinsky亞格子模型十分相似，通過過渡函數將兩種模型組合起來是比較好的構造脫體渦模擬(DES)的方式，已經有學者做過類似的工作，并取得了較好的效果[12]。孫東等深入分析了亞格子應力的組成以及彼此之間的關系后，結合現有的算法發展出了ML-DES算法[13-14]。袁先旭等在文獻[14]提出的混合模型的基礎上，引入壁面和網格拉伸比的影響，對亞格子尺度和過渡函數進行了修正；同時針對傳統模型采用突然變換的判斷函數，導致過渡不連續的問題，在Fujiik[15]提出的過渡函數基礎上，基于不可壓湍流邊界層特性，研究提出了一種新型過渡函數，形成了一種新的RANS-LES混合模型，具體可參考文獻[12]。

在RANS-LES混合模擬中，整體預測精度主要取決于LES部分的精度[16-18]。而LES對于計算格式、計算網格都有著較高的要求，如何正確配置計算格式和網格對于誤差控制具有重要的意義。鄧小兵[19]系統開展了LES計算中非線性多尺度系統的數值誤差分析，為數值格式的設計以及網格分布提出了指導原則。研究發現，顯式濾波能夠準確定義濾波可分辨尺度和亞濾波尺度，避免了網格尺度大小的空間離散誤差對大尺度運動的污染，能夠實現對數值誤差的有效控制；并且通過對計算過程中混淆誤差、色散誤差以及耗散誤差的分析(圖6[19])，發展了適用于RANS-LES計算的高精度格式構造方法。圖6中給出了使用多種高精度格式開展均勻各向同性湍流LES計算時所表現出的混淆差、色散誤差和耗散誤差影響，其中k表示空間波數，Δ表示網格尺度，E(k)為湍動能能譜，“CBCexp”為不同時刻實驗測得的能譜數據，“optC4”表示鄧小兵提出的一種4階中心型格式[19]，“Spectral-like”表示經過譜特性優化的中心型緊致格式(截斷波數為ωc)，“E2”和“E4”分別表示2階和4階中心差分格式，“Up3”和“Up5”分別表示3階和5階迎風偏置格式。

圖6 混淆誤差、色散誤差和耗散誤差對均勻各向同性湍流LES計算結果的影響[19]Fig.6 Influence of aliasing error, dispersion error and dissipation error on LES result of homogeneous isotropic turbulence[19]

使用超聲速圓柱底部流動問題對比了基于Spalart-Allmaras(SA)模型的模擬方法和基于RANS/LES混合模型模擬方法的計算結果。自由來流Ma∞=2.46，單位雷諾數Re/L=45×106/m，T∞=145 K，P∞=31 415 Pa。圓柱半徑R=31.75 mm，圓柱長為8R，計算網格規模約1 000萬。圖7(a)給出了SA模型和混合RANS/LES模擬得到的瞬時密度梯度分布云圖，可以看到混合RANS/LES計算結果清晰地呈現出充分發展的湍流脈動結構。圖7(b)給出了用壓力著色的底部渦量等值面云圖(壓力用來流動壓進行無量綱壓化)。可以看出，混合RANS/LES模型捕捉到了底部區域的主要流動結構，而SA模型僅能給出平均量。

圖7 超聲速圓柱底部繞流流動結構Fig.7 Flow structure around bottom of a supersonic cylinder

圖8為圓柱底部壓力系數Cp和速度U不同方法下的分布情況，其中r和x分別為尾跡區徑向坐標和流向坐標，U為流向平均速度，參考的實驗數據來自文獻[20]。相比于SA模型計算結果，混合RANS/LES模型得到的計算結果與實驗值更加接近。

圖8 RANS/LES混合模型與實驗和SA模型結果對比Fig.8 Comparison of results of RANS/LES method，experiment and SA model

1.3 熱流高效高精度模擬

高超聲速飛行器湍流數值模擬通常使用RANS方法，需要向壁面網格加密(y+≈1)，帶來計算收斂慢、計算穩定性下降等問題。因此，主流商業軟件，例如Fluent和CFD++等，對工程外形進行湍流模擬時通常使用壁面函數(通常20

文獻[21]總結了近50年壁面函數的發展，其中對于可壓縮流動，在標準壁面函數的基礎上，通過添加壓力梯度[22]、可壓縮性和壁面傳熱效應[23]等修正，擴展了壁面函數的應用范圍，但對于存在分離和再附等逆壓梯度的復雜流動，其適用性仍無法確定[24]。

文獻[25]發展的解析壁面函數通過在壁面附近解析求解簡化的湍流邊界層方程，得到速度和溫度的解析表達式，在低速分離流動中表現良好。文獻[26]基于可壓縮湍流邊界層的特征，在簡化方程中添加了對流項分布和能量方程中黏性耗散項的影響，發展了一種適用于可壓縮流動的解析壁面函數，可有效消除原始解析壁面函數在激波邊界層干擾算例數值模擬中出現的非物理振蕩和壁面熱流預測精度差等問題，接近密網格低雷諾數模型預測結果。

圖9比較了不同來流馬赫數時的可壓縮湍流邊界層[27-29]內的溫度分布，圖中Ma=5的實驗結果來自文獻[28]，Ma=8.2的實驗值來自于文獻[29]，LS model表示采用Launder-Sharma模型的計算結果，T表示溫度，Tw表示壁面溫度，T∞表示來流溫度，δ表示邊界層厚度，y表示壁面距離。通過圖9可以看出，隨著馬赫數的增加，壁面邊界層內溫度變化加快，導致流體層流黏性系數等參數變化。為此，發展了高超聲速湍流邊界層解析壁面函數方法，并集成到NNW軟件中，為工程上快速精確的湍流模擬提供了一種思路。

圖9 可壓縮湍流邊界層近壁面溫度分布Fig.9 Near-wall temperature distribution in turbulence compressible boundary layer

解析壁面函數在壁面附近通過可壓縮湍流邊界層假設簡化輸運方程，考慮對流項和壓力梯度的影響，其中對流項在壁面網格內變化較大，因此基于壁面網格點的對流項可使用無量綱壁面距離的平方項進行模擬，同時可壓縮湍流邊界層內能量方程黏性耗散項不能忽略。

通過給定湍流黏性系數分布，進行兩次積分可得到壁面速度和溫度表達式，進而得到可壓縮流動的解析壁面函數(MAWF)。結合解析壁面函數的構造過程和層流黏性系數分布規律，在黏性底層內構造一種拋物型(Parabolic)分布：

(1)

圖10比較了馬赫數Ma=8.2算例分離區前后壁面網格內層流黏性系數分布，其中密網格數值結果(LS)表明黏性系數在底層內變化較大，約30%，而本文構造的層流黏性系數相較于MAWF在壁面網格內使用壁面值，接近于密網格結果，合理地描述了強可壓縮流動中黏性底層內的變化。通過相似的積分過程，可以得到考慮層流黏性系數變化的解析壁面函數，即para-MAWF。

圖10 黏性系數分布(Ma=8.2)Fig.10 Distribution of viscosity variation (Ma=8.2)

斜激波邊界層干擾會產生激波反射、分離和再附等復雜流動現象，是考核壁面函數常用的算例，這里選取Ma=5.0和Ma=8.2兩個狀態進行比較。Ma=5.0算例的來流條件和網格等可參考文獻[26]。Ma=8.2算例[29]來流壓力為430 N/m2，來流溫度為81 K，壁面溫度為300 K，計算網格中包含激波產生器(β=10°)來模擬入射激波。計算中所有網格均經過網格無關性測試，最終低雷諾數湍流模型使用的密網格為480×105，分離區前y+≈0.6，壁面函數使用的粗網格為160×60，分離區前y+≈30。計算網格和模型尺寸如圖11所示。

圖11 Ma=8.2激波邊界層算例網格示意圖Fig.11 Sketch map of mesh for Ma=8.2 shock wave boundary layer calculation example

圖12比較了Ma=5激波產生器10°的Stanton數St分布，橫坐標中x1表示斜激波無黏反射位置，其中EXP為實驗數據[28]，AWF、SWF和MAWF分別表示原始解析壁面函數、標準壁面函數以及添加可壓縮修正的解析壁面函數的計算結果，para-MAWF為本文添加黏性底層層流黏性系數分布后的解析壁面函數數值預測結果，LS表示低雷諾數模型Launder-Sharmak-ε模型計算結果。AWF結果在分離起始位置出現非物理振蕩，主要原因是AWF對流項在壁面網格內的分布，將導致AWF預測能力降低。而MAWF認為對流項在壁面處為零，以無量綱壁面距離的平方項遞增，從物理上更符合邊界層內流動規律，從而提高預測精度。MAWF相較于AWF，壁面熱流的預測能力大幅上升，主要得益于能量方程簡化中保留了黏性耗散項，使得粗網格預測的壁面熱流精度大幅改善，接近密網格LS的結果。SWF在分離區前后均無法準確預測壁面熱流。添加層流黏性系數變化的解析壁面函數，對于本算例影響較小，接近于MAWF結果。圖13比較了不同湍流模擬方案消耗的時間對比，其中解析壁面函數消耗的時間僅為密網格LS的5%左右。

圖12 Ma=5.0算例壁面Stanton數分布Fig.12 Stanton number distribution for Ma=5.0

圖13 Ma=5.0耗費計算機時間對比Fig.13 Comparison of CPU time for Ma=5.0

圖14比較了Ma=8.2算例激波產生器10°的壁面熱流Q/Q∞結果，其中EXP-SB和EXP-TC均來自實驗數據[28]，分別表示使用塞式熱流傳感器(SB：Schmidt-Boelter gage)和熱電偶(TC：thermocouples)測量的壁面熱流，LS表示低雷諾數模型Launder-Sharmak-ε模型計算的結果，MAWF為文獻[26]發展的解析壁面函數計算結果，para-MAWF為本文添加黏性底層層流黏性系數分布后的解析壁面函數數值預測結果，從圖中可知相較于MAWF，考慮層流黏性系數分布的解析壁面函數得到壁面熱流在分離區顯著降低，更接近于實驗值，相較于密網格LS模型降低了約30%。密網格LS結果嚴重高度壁面熱流，和文獻[30]中密網格低雷諾數模型結論一致，均出現高估高超聲速算例壁面熱流現象。

圖14 Ma=8.2算例壁面熱流分布Fig.14 Heat-flux distribution on wall for Ma=8.2

1.4 轉捩過程高精度模擬

在過去的20多年里，亞利桑那大學的Fasel課題組采用直接數值模擬方法對0°攻角錐表面第二模態的自然轉捩過程開展了大量的研究[31-35]，他們發現在第二模態的非線性階段會出現流向條帶。而在橫流效應下(雖然在迎風對稱面附近的橫流速度較小)，迎風中心線附近第二模態的自然轉捩與0°攻角鈍錐表面存在何種異同是值得研究的課題。因此，本文采用直接數值模擬，研究了飛行工況下迎風中心線上最不穩定Mack模態的線性和非線性演化過程，并將其與0°攻角圓錐邊界層中Mack模態的自然轉捩過程做了對比。

如圖15所示，研究對象為一個半錐角7°、頭部半徑5 mm、長度700 mm的鈍錐。來流條件為飛行工況：馬赫數Ma∞=6、單位長度雷諾數Re∞=2.079×107/m、溫度T∞=216 K、攻角α=5°。壁溫Tw=500 K。直角坐標系(X,Y,Z)和貼體坐標系(ξ,η,θ)中的速度分量分別為(u,v,w)和(un,vn,wn)，周向角θ在背風和迎風對稱面分別為0°和180°。

圖15 坐標系和模型示意圖Fig.15 Sketch of coordinates and computational model

數值模擬主要分為4步:

第1步采用二階精度的有限體積方法獲得包含頭部和激波的基本流，網格尺寸為Nξ×Nη×Nζ=1 000(流向)×401(法向)×181(周向)。第1層網格離壁面的距離從駐點的4×10-3mm線性增長到計算域出口的8×10-3mm。

第2步為減少網格與激波不平行帶來的數值振蕩，利用第1步的流場獲得激波位置，生成與激波平行的網格，并采用和第1步相同的格式和網格量重新計算層流場。在前兩步中，考慮到流動的對稱性，計算僅使用半模。

第3步截取一個不包含激波和頭部的子計算域，采用高階有限差分獲得密網格下收斂的層流解，網格大小為Nξ×Nη×Nζ= 8 000×240×380。子計算域的入口位于x=150 mm處，在該流向位置，背風和迎風對稱面上遠場離壁面的距離分別為當地邊界層厚度的2.5倍和8倍。入口、遠場邊界條件以及初始場從上一步中插值得到。由于關心的是迎風面擾動的演化過程，為減少計算量，整個背風面為周向緩沖區，周向網格從側方對稱面向背風對稱面逐漸稀疏，周向緩沖區網格點總數為60，因此，在迎風面，網格點的周向間距為1°，同樣，流向也設置緩沖區(約300 mm)，網格量為100。在周向也采用了類似的網格分布[36]。法向網格采用指數拉伸，在每個扇面(k=1, 2,…，Nζ)，網格點離壁面的法向距離滿足：

(2)

式中：Lη(i,k)是當地遠場到壁面的法向距離;i=1, 2,…,Nξ;j=1, 2,…,Nη;參數b通過第1層網格離壁面的距離η(i, 2)=2.5×10-3mm來確定。邊界層內網格點約為120～130。

第4步在第3步的層流場中施加非定常壁面吹吸。吹吸形式為

vbs=Asin((x-x1)/(x2-x1))3rand(z,t)

(3)

式中：vbs為沿物面法向的速度；x1和x2為吹吸槽在流向的邊界，吹吸幅值A為來流速度的10-3。圖16為線性穩定性得到的不同頻率Mack模態的N值(增長率沿流向的積分)沿流向的變化。頻率f=1.5 MHz的Mack模態達到的N值最大，其中性點在x=300 mm附近，因此，將吹吸頻率設為1.5 MHz，吹吸槽位于x1=280 mm和x2=282 mm之間。

圖17為迎風面附近的時均摩阻系數Cf，很明顯在模型的下游出現了條帶，這與0°攻角錐在風洞工況下的結果是類似的，這表明，雖然有攻角時邊界層是三維的，但是第二模態的失穩形式與二維高超聲速邊界層是類似的，即都是基頻共振(Fundamental Resonance)。

圖16 線性穩定性理論得到不同頻率Mack模態N值沿流向的變化Fig.16 Variation of N-factors of Mack modes at different frequencies in streamwise direction by linear stability theory

圖17 迎風面附近的時均摩阻Fig.17 Temporal-averaged skin friction around windward plane

圖18為x=450 mm和500 mm處時均摩阻在周向的相關函數，橫坐標為周向角度的間隔。結果表明條帶之間的間距約為2°, 按弧長計算約為邊界層厚度的2～3倍。這比飛行工況下二維湍流邊界層中第二模態轉捩引起的條帶之間的間距要小，雖然在不同工況下條帶之間的絕對距離沒有可比性。

圖18 x=450和500 mm處時均摩阻周向的相關函數Fig.18 Azithumal correlation function of temporal- averaged skin friction at x=450,500 mm

圖19為轉捩前期溫度的典型均方根及其頻譜，圖19(a)給出了x=400 mm處溫度的時均值(虛線)和均方根(云圖)，其中縱坐標η表示距物面的距離，橫坐標為周向角；圖19(b)針對圖19(a) 中的3個采樣點“△□○”，給出了采樣點處的溫度脈動頻譜。結果表明，溫度均方根在近壁面和邊界層附近存在兩個極值，這表明轉捩確實是由第二模態失穩引起的。從圖19(b)還可以看出，第二模態的頻率約為1.6 MHz，這與圖16中穩定性分析得到的結果吻合。

圖19 x=400 mm處溫度的均方根及其頻譜分布Fig.19 Distribution of RMS and spectra of temperature at x=400 mm

2 復雜干擾機制高清晰認知

2.1 鈍舵縫隙局部流場數值模擬

高超聲速飛行器在穿過大氣層時會面臨嚴酷的氣動加熱[37-38]，鈍舵作為高超聲速飛行器典型的凸出部件，外部的高氣流繞過鈍舵時，會在鈍舵前產生激波。而激波打到飛行器物面上會發生激波-邊界層干擾，同時鈍舵與機體交匯處通常會發生分離與再附，進而形成一系列局部高熱流區。為了完成變舵偏等飛行動作，鈍舵底部往往會留有一定的縫隙結構，縫隙的存在會使得鈍舵結構變得更加復雜[39-40]，復雜的流動結構對于鈍舵縫隙的數值模擬提出了較大的挑戰。與此同時，由于鈍舵前緣、舵軸及安裝面常常是燒蝕的“重災區”，整個鈍舵縫隙結構是熱防護設計的重點。

高超聲速鈍舵縫隙數值模擬有兩個難點，一是縫隙尺寸和飛行器尺寸之間相差巨大，這對網格的生成提出了很高的要求;二是縫隙內的黏性流動和縫隙外的高超主流在流動速度上相差巨大，這對計算格式提出了很高要求。在計算網格方面，需要根據豐富的熱環境網格繪制經驗，繪制過渡充分均勻光滑，近壁面正交性和尺度良好，空間無“邊界層”密度的結構網格。在計算格式方面為了能在捕捉激波的同時，很好地模擬縫隙的低速黏性流動，采用了Van Leer[41]加焓守恒修正格式。Van Leer分裂具有很好的激波捕捉能力，但它的缺點是在邊界層、剪切層計算中格式耗散過大，降低了黏性流動模擬精度。因此，在Van Leer通量分裂的基礎上引入了焓守恒修正[42]，減小了邊界層、剪切層中的格式耗散。有效抑制了非物理耗散，在一定程度上解決了流場分辨率與計算穩定性之間的矛盾。

選用李艷麗和李素循[43]的無縫隙鈍舵作為計算模型，模型結構如圖20所示，詳細尺寸和來流參數見文獻[43]。

圖20 鈍舵模型[43]Fig.20 Geometry of blunt rudder[43]

圖21給出了鈍舵計算流場結構與試驗流場結構的示意圖，圖22給出了鈍舵前緣中心線平板熱流對比圖，縱坐標為無量綱化的熱流，其中參考熱流qref為相同來流條件下，沒有鈍舵干擾時對應鈍舵位置的平板熱流值，橫坐標為以舵軸前緣為零點，用鈍舵前緣直徑進行無量綱化的平板中心線位置無量綱值。從圖21、圖22可以看出計算結果無論從流場結構還是表面熱流均與文獻[43]的試驗結果較好吻合，這在一定程度上驗證了計算結果的可靠性。

圖21 鈍舵流場結構對比Fig.21 Comparison of flow structure around blunt rudder

圖22 鈍舵前緣平板熱流對比Fig.22 Comparison of heat flux at leading edge of blunt rudder

為更好符合實際飛行器鈍舵情況，基于無縫隙鈍舵，建立帶縫隙鈍舵簡化模型，舵軸直徑12 mm 位于舵正中，舵縫隙高4 mm。網格第1層高度取0.001 mm，網格法向增占比控制在1.1以內。圖23給出了鈍舵縫隙對稱面馬赫數云圖和二維流線分布。從圖中可以看到，繞鈍舵流動，舵前緣形成脫體激波。平板邊界層發生分離形成的分離激波與脫體激波碰撞形成復雜的波系結構，舵根部上游區域出現了馬蹄渦結構，如圖23(a) 所示。其次，鈍舵前緣流體填充至縫隙內部，從馬赫數云圖可以看出，縫隙內部流體速度迅速降低，變為低速的黏性流體。由于分離再附的作用，低速黏性流體在縫隙的上底面和下底面形成了多個旋渦結構，如圖23(b)所示。

圖23 對稱面馬赫數和流線分布Fig.23 Distribution of Mach number and stream line on symmetric plane

為了更好觀察發現流場結構，分析流動特性，從鈍舵縫隙模型前緣到后緣，以舵軸為對稱面截取了7個垂直截面，沿著流動方向依次為P1～P7截面，將縫隙上下物面熱流與7個縫隙截面二維流線進行對照比較，如圖24、圖25所示。對于舵軸上游截面，各截面內均出現了渦串結構。從P1截面開始出現四渦結構，到達舵軸上游縫隙中部位置附近(P2截面)，渦的數量多達6個，相應縫隙上下物面再附區域橫向有明顯的馬蹄形高熱流區，這主要是由于來流受舵縫和舵軸的影響，在上下表面多次出現分離再附引起的。從P2～P3截面可以看出，縫隙底部旋渦向舵面兩側發展，并很快發展到兩側的舵面上，并在舵側面形成了一條高熱流帶。從P4截面來看，受舵軸的影響，縫隙內旋渦結構在舵軸處基本消失。對于舵軸下游縫隙區域，P5、P6和P7截面縫隙內部流線較為簡單，物面沒有熱流高值區，舵軸后緣由于是背風區，因而出現了熱流低值區。縫隙外部馬蹄渦沿著舵面不斷向上抬升，到鈍舵的后緣附近(P7截面)，旋渦結構發展到最大，其影響范圍也最遠，同時受旋渦運動的影響，縫隙上底面側棱附近出現了一條明顯的高熱流帶。總的來看，流場結構清晰合理、壁面熱流光滑連續，流場結構基本可以解釋壁面熱流分布規律，計算結果具有較高可信度。

圖24 下表面不同站位附近熱流和流線分布Fig.24 Heat flux and stream line at different cross section on lower surface

圖25 上表面不同站位附近熱流和流線分布Fig.25 Heat flux and stream line at different cross section on upper surface

2.2 進氣道激波振蕩的物理機制

吸氣式高超聲速飛行器進氣道經常面臨激波/邊界層、激波/激波相互作用等[44]，使得進氣道唇口產生非定常的激波振蕩，而進氣道內激波串的相互作用也可能產生自激振蕩的現象[45-46]，嚴重時會使得進氣道壅塞并導致發動機不啟動現象[47-50]，需要加以控制[51-52]。目前進氣道大都采用前體/進氣道一體化設計以得到更好的發動機推力性能和飛行性能[12]。有必要對整機采用內外流一體化計算，考察不同來流工況下，進氣道激波振蕩機制以及飛行器整體氣動力特性的變化，支撐進氣道的優化設計。

為了模擬飛行器在不同工況下的進氣道激波振蕩過程，基于NNW通用軟件，開發了內外流一體化RANS-LES混合模擬方法，提高了對進氣道分離流動及復雜波系感染的分辨能力。

計算模型采用了某工程升力體外形，計算狀態為：高度20 km，來流溫度為216.65 K，攻角和側滑角均為0°；為了研究來流速度對進氣道波系結構的影響，設置4種來流馬赫數：Ma=1.5， 3.0， 5.0，7.0。計算網格在進氣道和壁面處進行了加密，網格規模2 000萬。

吸氣式飛行器進氣道的流動受到多種因素的影響，比如來流速度、飛行高度、飛行攻角等。本節通過設置不同的來流馬赫數，研究來流速度對進氣道激波振蕩過程的影響。圖26給出了來流Ma=1.5的流動結構，主圖為密度梯度分布云圖，并給出了紅色方框內的流線分布圖。可以看到，當飛行馬赫數較低時，進氣道上壓縮面邊界層發生流動分離，并產生大尺度的分離泡；分離泡周期性地變大或縮小并產生低頻振蕩；繼而與斜激波相互作用，使得斜激波在唇口來回運動并同樣產生低頻振蕩。在進氣道入口附近產生一道正激波，正激波后壓力較高，使得進氣道壅塞。飛行器的氣動力整體上也呈現出低頻振蕩特征。

圖26 Ma=1.5流動結構Fig.26 Flowfield structures at Ma=1.5

當飛行馬赫數增大至3.0時，如圖27所示，分離泡結構發生顯著變化，在進氣道上、下壓縮面均產生了分離泡，部分斜激波打入進氣道內部。在進氣道入口附近產生馬赫桿結構，來流只能從馬赫桿兩側進入進氣道，使得進氣道部分壅塞。唇口激波低頻振蕩，進氣道內激波串高頻振蕩。此時的氣動力整體上呈現出低頻振蕩與高頻振蕩共存的狀態。

圖27 Ma=3.0流動結構Fig.27 Flowfield structures at Ma=3.0

馬赫數進一步增大至5.0時(圖28)，只在進氣道上壓縮面存在分離泡，且位置更加靠后，斜激波可以直接打入進氣道內部，可以正常進氣。進氣道內激波串相互作用，形成高頻振蕩。氣動力整體上呈現出高頻振蕩的特征。馬赫數7.0時的流動結構與馬赫數5.0相似。

圖28 Ma=5.0流動結構Fig.28 Flowfield structures at Ma=5.0

總體來說，隨著飛行馬赫數不同，入口激波形狀、激波振蕩頻率及氣動力振蕩頻率均存在差異。馬赫數變化會導致激波高低頻振蕩的轉化。在飛行器確定設計工況時，需要考慮飛行馬赫數對進氣道進氣性能的影響。

2.3 噴流干擾誘導回流激波機制

高超聲速飛行器在臨近空間環境高速飛行時，由于來流大氣稀薄、來流動壓低導致控制舵面效率不足，使得采用單一的舵面姿態控制系統難以滿足飛行器機動性和彈道控制的需求，必須采用噴流控制或噴流/舵面復合控制[53]。橫向噴流由于結構設計簡單、不工作時對流場影響小、氣動響應迅速、控制效率高等特點[54-57]，被廣泛應用于高超聲速飛行器氣動控制系統。但在高超聲速來流條件下，橫向噴流與來流會形成很強的干擾流場結構，會在噴口上下游出現邊界層分離與再附、自由剪切層失穩，多激波系、多旋渦系以及激波/邊界層干擾等復雜流動現象[58-60]。

國外學者早在20世紀50年代末就開始了對側向噴流控制技術及其干擾機理的研究，國內相關工作起步較晚，始于20世紀80年代。研究采用理論分析、工程估算、風洞試驗和數值模擬的方法，分析了噴口參數(位置、形狀、數量等)、噴流參數(馬赫數、流量、噴流壓比、噴流介質和組分等)、來流參數(馬赫數、雷諾數、迎角等)對橫向噴流干擾特性的影響[61]，從流場建立過程方面分析了噴流非定常流場建立與消退過程[62]，從物理化學過程方面分析了真實氣體效應對橫向噴流干擾的影響[63]。這些研究較為系統地給出了橫向噴流流場特征和噴流干擾機理。但同時也注意到，橫向噴流多位于模型前體和中部，位于模型后體的情況研究較少，特別是扁平升力體類構型，由于發動機安裝位置限制，噴流置于飛行器底部的橫向噴流研究少見；噴流和舵面同時工作時，噴流對舵效影響的研究較少。

基于NNW軟件，針對底部橫向噴流開展數值模擬研究，分析底部橫向噴流流場特性及對舵面效應的影響規律，為相關飛行器控制系統設計提供技術支撐。研究對象為升力體構型，控制機構包括底部的控制舵和多通道噴管組合結構，噴管布局示意見圖29。

置于底部的噴流發動機采用單噴/多噴的組合實現對俯仰、偏航、滾轉3個通道的控制。1#噴管為俯仰控制，4#噴管為偏航控制，2#和3#同時開啟為滾轉控制。底部橫向噴流最顯著的特征是噴流出口壓力與當地壓力比值極大，高超聲速飛行器底部壓力接近真空，噴流開啟之后，燃氣氣體極速膨脹，對底部流場結構和附近的舵面流場結果產生較大影響，進而影響整體氣動特性和舵面控制效率。

2.3.1 底部橫向噴流干擾流場分析

數值模擬了3個控制通道的典型流場，計算網格采用結構對接網格，全場網格量約2 000萬，圖30給出了1#噴管開啟時噴管附近的流場結構。高壓氣流在高速離開噴管后迅速膨脹，受彈體底部和腹鰭物面的限制，在這些部位附近形成復雜的激波結構，導致底部及下表面形成相對的高壓分布，在下表面的物面附近產生強逆壓梯度，導致流動發生分離，在腹鰭底部附近的部分區域形成回流。下表面的高壓區起到了增大低頭力矩的作用。在氣流到達底部及腹鰭邊沿緣后，與前體激波層內的超聲速氣流發生碰撞而在各自的前進方向形成激波，同時又和激波層外的高超聲速來流相互作用形成強激波，從而整個區域呈現顯著的激波/激波干擾現象，形成范圍很大、強度很強的弓形激波結構，導致局部區域的高溫分布(最大值達9 000 K以上，如圖31所示)，在如此高溫條件下，真實氣體效應可能會變得尤為顯著。

圖30 1#噴管開啟的流場結構Fig.30 Flow field structure at the open of 1# nozzle

圖31 底部截面溫度分布Fig.31 Distribution of temperature on bottom

2.3.2 噴流對舵面效率的影響

從流場結構分析來看，噴流形成的干擾區距離升降舵較近，甚至出現噴流燃氣直接打到舵面的情況，將會對舵面壓力產生較大影響。表1分析了3個通道噴流對升降舵舵效的影響程度，噴流開啟后升降舵舵效有明顯變化，俯仰和滾轉通道方向為迎風面噴流，顯著增大舵面迎風面壓力，使得舵面效率增大約20%，而偏航方向噴流產生的吹吸效應對舵面壓力的影響比較復雜，會使得舵面效率降低約10%。噴流對舵面效率的顯著影響，需要在飛行器控制時引起關注。

3 多學科耦合應用

3.1 方形截面導彈舵面操縱動態響應特性

可操縱性是氣動設計的主要內容之一[64]，開展飛行器在舵面操縱下的動態響應過程模擬，獲取動態響應參數，對提升飛行器的動態響應品質至關重要。如研究飛行器在升降舵階躍舵偏操縱下攻角的振蕩過程，可以獲取超調量、峰值時間、調整時間以及上升時間等時域響應特性[65]；而研究諧波型舵偏操縱下飛行器運動過程，可以獲取響應帶寬、諧振峰值等頻域的響應特性。上述時域和頻域響應參數都是飛行器設計的重要指標。

模擬飛行器在舵面操縱下的動態響應過程，主要應用到NNW軟件的動網格技術和非定常流動/飛行器運動一體化耦合模擬技術。

計算選用方形截面導彈外形(圖32)，其主要特點是雷達反射面小，有利于隱身；可充分利用內埋彈艙空間，提高載彈量；平面外形能提供更高的升力等。針對該外形，分析了飛行器俯仰運動對階躍舵偏操縱的動態響應特性，半場計算網格約350萬。

圖32 方形截面導彈外形和表面網格Fig.32 Geometry and grid of square cross section misssile

舵面的偏轉按以下函數形式給出：

(4)

式中：δ0和δm分別為初始舵偏角和舵偏幅值，δ0=0°，δm根據需要選取。在實際計算中，為減少初始流場建立過程對計算的干擾，開始階段飛行器先不引入舵面運動，飛行器在自身氣動力作用下自由運動，在0.56 s時再引入舵面偏轉，在0.61 s時舵面回復到原位置。

圖33 轉動慣量對導彈動態響應特性影響Fig.33 Effect of rotational inertia on missile dynamic response characteristics

圖34比較了質心位置對動態響應特性的影響。在xc/L=0.47時，飛行器在攻角±10°范圍內是靜不穩定，0°配平攻角表現為鞍點；超出這個攻角范圍之后，導彈轉變為靜穩定的，兩個配平攻角±13°是穩定的結點。從圖34(a)模擬結果可以看到，從配平攻角釋放后，當舵偏引入的擾動較小時(δm≤8°)，導彈的運動圍繞配平點附近振蕩并逐漸收斂；但當擾動較大時(δm=12°)，俯仰運動會跨過靜穩定區域，進入靜不穩定區域，在較大攻角范圍內持續振蕩，運動特性極為復雜，不利于導彈的操縱。為了改善導彈的響應特性，可以通過調整質心，改變飛行器的靜動態氣動特性。圖34(b)展示了質心調整到xc/L=0.42時，導彈的動態響應過程。調整質心后，飛行器只有唯一的配平攻角0°，且在較大攻角范圍內都是靜穩定的，從模擬結果來看，即便是較大的擾動，飛行器在擾動消失后，其俯仰運動仍然是收斂的。

圖34 質心位置對導彈動態響應特性影響Fig.34 Effect of center of mass position on missile dynamic response characteristics

3.2 大長細導彈快速拉起過程的氣動伺服彈性特性

為追求“更快速、更靈活、更精準”等戰術指標，以空空導彈、火箭彈等為代表的新型彈箭呈現出更大的長細比、更輕的結構占比和更大的過載等主被動設計特點，從而導致其彈性變形越來越突出[66-67]。結構低階固有頻率與飛行力學短周期模態比較接近，氣動彈性和控制系統之間的耦合作用顯著增強，會不同程度影響控制系統的精度與飛行品質，甚至導致氣動伺服彈性失穩問題[68]。因此，在飛控系統設計中，對飛行器進行氣動伺服彈性分析，已經成為現代飛行器設計過程中必不可少的重要環節[69-71]。以往工程設計中，多采用定常線化假設和剛體運動/彈性振動解耦模型在頻域范圍內開展分析[72]，鑒于此類彈箭飛行姿態角等會在瞬間發生劇烈變化，導致飛行動力學問題研究中最為關鍵的氣動特性呈現出非線性和非定常特征，同時多學科之間顯著耦合效應也使得以往的氣動力簡化方法難以準確反映飛行器氣動載荷的實時動態變化[73]。

開展大長細彈箭快速機動過程中的氣動/結構/控制耦合問題的精細模擬，主要需解決3個關鍵技術：同時適應物面大尺度剛體運動和小尺度彈性變形的網格動態邊界處理技術，能夠獲取高精準氣動力/載荷特性的高效高魯棒非定常流場CFD求解技術，以及能同時包含剛體運動和結構振動耦合效應的高保真動力學模型等。基于NNW通用軟件的CFD求解模塊和動網格技術，通過二次開發發展了同時包含氣動耦合和慣性耦合的剛體運動/彈性振動耦合動力學模型(如方程(5)所示)[74]，并引入PID控制環節(如方程(6)所示)，建立了基于CFD的彈性飛行器氣動伺服彈性問題數值仿真方法。

(5)

(6)

式中：m、I0、V和ω為分別為飛行器的總質量、轉動慣量矩陣、質心速度以及繞質心的角速度;ωss表示矢量ω的斜對稱矩陣;F和M為作用在飛行器上的合外力和相對于質心的合外力矩;hij為第i階和第j階結構振型的協變量矩陣;ηi和ηj為對應第i階和第j階結構振型的廣義位移;Mij和Mii為廣義質量矩陣的元素;Qi為對應第i階結構振型的廣義氣動力;δcs為操縱面偏角值;KP、KI和KD分別為PID控制器比例、積分和微分環節的增益值;α和αcmd分別為飛行器的實時迎角和目標指令迎角;q為飛行器繞質心的俯仰角速度;θe為陀螺儀所在位置由結構變形引起的俯仰方向變形角。

計算網格包括彈身的靜止域網格以及兩片升降舵的運動域網格(圖35)，圖36給出了3套體網格經過嵌套網格方法處理前后不同位置的切片，可看出所用的嵌套網格方法生成的網格嵌套區域和嵌套邊界都比較合理。本文設計的來流馬赫數為3.0，飛行高度設定為10 km，基于彈身直徑的雷諾數為9.16×106。為了描述模型的結構彈性變形，這里選取模態分析所獲取的具有代表性的4階模態用于彈性特性計算，包含彈身的前兩階彎曲模態和升降舵的一階彎曲與一階扭轉模態，具體數據見文獻[75]。

圖35 計算網格Fig.35 Computational mesh

圖36 嵌套網格處理前后的計算網格切片Fig.36 Clips of computational grid before and after nesting operation

仿真導彈以馬赫數3.0的速度飛行，在t=0 s 時刻給定20°迎角的拉起指令，不同模型的響應歷程如圖37所示。不考慮彈性效應時，圖37(a)的仿真結果表明，該導彈在所設計的控制環節作用下飛行迎角能夠在0.5 s左右達到控制指令，且不存在超調量，可以實現穩準快的操縱要求。而同時考慮結構振動對閉環控制系統的影響時，在不加入振動抑制措施的情況下，如圖37(b)所示，彈性模型的飛行迎角也能夠短時間內快速趨向于指令要求，但之后會在目標迎角附近振蕩，且呈不斷發散特征，這表明原來按照剛體模型設計的控制律應用到彈性模型時，由于納入彈性振動效應，原來穩定的閉環飛行動力學系統會出現典型的氣動伺服彈性發散問題。

考慮到結構在受到動態載荷激勵后必然產生彈性振動，在操縱反饋過程中，傳感器所采集到的運動信號中包含的振動信號會對導彈的實時剛體運動速度、質心位置等平動特征和角速度、姿態角等轉動特征造成“誤判”。圖37(c)給出了只考慮彈性效應對平動特征的影響，即在控制環節中去掉了彈性效應對傳感器轉動信號的干擾，可以看出，飛行迎角能夠在0.6 s內快速振蕩收斂到目標，這與剛體模型達到目標指令的響應時間基本一致。

圖37 不同模型升降舵偏角和迎角的響應歷程Fig.37 Response history of control surface’s deflection angles and angles of attack for different models

從飛行力學特征模態角度分析相關機理，俯仰角變化屬于飛行器短周期模態的運動特征。對于短周期模態，由于運動歷程發展迅速，需要控制系統做出快速操縱反應，因此在飛控系統設計時除了要求其穩定外，還需保證合適的阻尼。考慮到彈體一階彎曲模態的固有頻率也與短周期模態頻率更接近，特別是引入反饋控制環節后，彈性振動對傳感器采集到的角速度等造成了干擾，從而使對穩定性有苛刻要求的短周期模態產生耦合失穩現象。相比之下，質心運動屬于導彈長周期模態的運動特征，由于發展緩慢，只要不出現非周期發散情況，也容易被控制系統來操縱。因此，即使傳感器納入了彈性振動對速度和加速度的干擾，也不容易對長周期模態的穩定性產生定性的影響，只是定量改變了響應的歷程。

3.3 考慮化學非平衡效應的動態穩定性計算

飛行器的動態穩定性是指飛行器受到擾動后保持原來飛行狀態的能力，通常用動導數進行表征，是動態氣動特性設計、彈道設計及動態穩定性分析中的關鍵參數。高超聲速飛行器再入時，其周圍空氣溫度可達上萬度，氣體在高溫作用下發生離解、電離，并發生復雜的化學反應，出現化學平衡/非平衡效應和熱力學平衡/非平衡效應，即所謂的高溫真實氣體效應，對飛行器力/熱特性、目標特性以及動態特性等產生復雜影響。

研究真實氣體效應對飛行器動態穩定性的影響，涉及熱化學反應、非定常流動以及飛行器運動等多種特征時間尺度問題的耦合模擬。鑒于數值模擬的難度，以往的研究多關注平衡氣體對動態穩定性的影響，對化學非平衡氣體則較少關注[76-78]。NNW軟件建立了比較完備的真實氣體效應模擬能力，在自主開發的軟件框架下有機集成非定常、運動等功能，具備化學非平衡效應對動導數影響的研究能力。在真實氣體效應模擬方面，NNW軟件目前具有熱力學一溫度和兩溫度模型，可以考慮5、7、11組分等不同氣體組分的化學反應過程，壁面催化特性包含完全催化、完全非催化以及有限催化等。在研究化學非平衡效應對飛行器動態穩定性影響時，化學反應模型采用Park的7組分模型，采用強迫振蕩方法計算飛行器的動導數。

圖38比較了考慮真實氣體效應的返回艙標模熱流計算結果，圖中標注為LENS的曲線是實驗結果，標注DPLR的為國外軟件數值計算結果，均來自文獻[79]，標注FCW、NCW和r=0.01的曲線為NNW軟件的計算結果，其中，FCW為完全催化條件計算結果，NCW為完全非催化條件計算結果，r=0.01為有限催化條件計算結果，計算結果與文獻值吻合很好。圖39比較了ELECTRE鈍錐外形采用完全氣體和化學非平衡模型計算的對稱面溫度分布云圖，計算馬赫數Ma=13.0，高度=53.5 km，計算網格量約為300萬。從圖39可以看出，考慮化學非平衡效應后，由于化學反應本身是一個吸熱的過程，計算得到的激波層內氣體的溫度明顯降低，同時激波位置也更貼體。

圖39 鈍錐溫度分布云圖Fig.39 Temperature contour of blunt cone

圖40 俯仰力矩遲滯圈比較Fig.40 Comparison of hysteresis loop for pitching moment

4 結 論

1) NNW軟件集成了較為完備的轉捩模型、湍流模型、RANS-LES混合模型，開發了WENO、WENN等高精度算法，基本可同時實現強激波捕捉和邊界層高精度模擬。

2) 針對不同應用需求，用戶可基于NNW軟件進行二次開發，滿足特殊場景的用戶需求。

3) NNW軟件具有良好的架構設計，可實現跨學科之間的耦合模擬，如化學非平衡流動、飛行器運動、結構振動以及控制系統等之間的耦合模擬。

致 謝

向應用團隊成員王新光、董思衛、華如豪等深表謝意，感謝NNW軟件開發團隊、測試團隊和推廣應用團隊的技術支持。