陳 浩，袁先旭，畢 林，華如豪，司芳芳，唐志共，傅亞陸，*

(1. 空氣動力學(xué)國家重點(diǎn)實(shí)驗室，綿陽 621000；2. 中國空氣動力研究與發(fā)展中心，綿陽 621000)

0 引 言

新一代高機(jī)動飛行器研發(fā)過程中，更高的戰(zhàn)技性能需求也給空氣動力學(xué)研究提出了許多新的挑戰(zhàn)。例如，國際上的第五代戰(zhàn)斗機(jī)[1-2]有許多區(qū)別于以往戰(zhàn)機(jī)的重要技術(shù)特點(diǎn)，其中最重要的有隱身(Stealth)、超聲速巡航(Supersonic Cruise)、超機(jī)動飛行性能(Supermaneuverability)和短距起降性能(Short takeoff and landing)，統(tǒng)稱為“4S”能力。“4S”性能的實(shí)現(xiàn)不僅需要高性能的推力矢量技術(shù)，還需要優(yōu)良的氣動布局以及對大迎角狀態(tài)下前體渦演化的理解和控制。

目前，先進(jìn)戰(zhàn)斗機(jī)如F-22和F-35均采用脊形前體，以降低可探測性和實(shí)現(xiàn)良好的超聲速性能。脊形前體的脊形邊緣能吸收雷達(dá)波，有效減少雷達(dá)散射截面積（RCS），同時，脊形前體高度混合外形也有利于超高速飛行。此外，由于脊形前體的背風(fēng)流場的渦結(jié)構(gòu)很強(qiáng)[3]，且較穩(wěn)定，與機(jī)翼前緣渦相互干擾，常常能夠延遲機(jī)翼的完全失速[4]，增加最大升力，并能在有側(cè)滑時提供穩(wěn)定的滾轉(zhuǎn)力矩，因此提高了飛機(jī)的橫向穩(wěn)定性。著名的Su-27飛機(jī)的“眼鏡蛇”機(jī)動就是利用上述原理提升的橫向/方向穩(wěn)定性來擴(kuò)大其大迎角性能，度過機(jī)翼失速[5]。此外，穩(wěn)定的渦流場也能提高飛機(jī)平飛和大迎角飛行時的橫向/方向穩(wěn)定性[6]。

1963年，脊形前體首次用在洛克希德A-12偵察機(jī)（Lockheed A-12）上，后續(xù)逐步推廣應(yīng)用于YF-12直升機(jī)、D-21無人偵察機(jī)和SR-71偵察機(jī)[7]。盡管國外對脊形前體的靜態(tài)氣動特性已開展了大量研究工作，并應(yīng)用于工程實(shí)際，但由于風(fēng)洞試驗和數(shù)值模擬方法等研究手段的局限，對其大迎角狀態(tài)下流動演化機(jī)理和影響的研究尚不充分，而深入了解脊形前體大迎角氣動特性，對于工程應(yīng)用和學(xué)科發(fā)展而言，均有明確的研究需求。

本文瞄準(zhǔn)上述需求，針對脊形前體大迎角湍流大分離流動特點(diǎn)，開展基于DES類RANS/LES湍流數(shù)值模擬方法的應(yīng)用研究，并以此為基礎(chǔ)，重點(diǎn)研究了橫截面形狀、迎角等參數(shù)對脊形前體的大迎角氣動特性和空間渦流場結(jié)構(gòu)演化的影響規(guī)律，為更好地了解其動態(tài)特性和進(jìn)行先進(jìn)高機(jī)動飛行器設(shè)計提供參考。

1 數(shù)值模擬方法

在脊形前體渦流場以及前體渦與翼型渦干擾的復(fù)雜流場中，存在湍流脈動、流動分離等，它們在流場中占據(jù)重要位置，直接影響到流場渦結(jié)構(gòu)的形成、演化、破碎與脫落等。

1.1 湍流方程和模型

一般曲線坐標(biāo)系下，湍流方程可以寫作統(tǒng)一的形式：

式中：Φ為湍流變量（如k、ω、v等），F(xiàn)、G、H為對流通量，F(xiàn)v、Gv、Hv為黏性通量，S p、S D、D分 別為生成項、破壞項和擴(kuò)散項。

由于本文研究的飛機(jī)前體大迎角分離流動的復(fù)雜性，綜合目前各種湍流方法的優(yōu)缺點(diǎn)，確保模擬結(jié)果能準(zhǔn)確地反映氣動力非線性特性，本文采用混合RANS/LES方法中的IDDES方法[8-10]，該方法結(jié)合了DDES和WMLES方法[8]，可有效解決DES方法中存在的對數(shù)層不匹配[8]的問題，同時可節(jié)省計算量，并且在包含壁面的復(fù)雜流動問題上，能得到比DDES更好的結(jié)果。相對于通常的LES和（D）DES，IDDES采用了新的亞格子尺度[10]，其與網(wǎng)格大小和壁面距離皆有關(guān)，可實(shí)現(xiàn)湍流求解模型的轉(zhuǎn)換。該模型引入混合函數(shù)[10]來修正由于RANS和LES交界面相互作用而損耗過多的雷諾應(yīng)力，提高模型保真度。近年來，很多學(xué)者對于該方法進(jìn)行研究、應(yīng)用和發(fā)展，并取得了大量成果[11-13]。

1.2 數(shù)值離散

迎風(fēng)型NND格式以其形式簡單，計算量小，且數(shù)值耗散較小，穩(wěn)定性較好等優(yōu)勢，得到了廣泛的應(yīng)用。本文數(shù)值模擬中對無黏項的空間離散都采用迎風(fēng)型NND格式[14]。

式中，下標(biāo)i代表網(wǎng)格單元中心位置，i±1/2代表單元界面位置。

對于方程（2），界面無黏通量為：

其他方向的無黏通量也是類似的表述。

可以看出，本文采用的空間離散格式為二階精度。對于具有旋渦、分離等特征的復(fù)雜流動問題，采用混合算法或高精度格式可以改善格式耗散特性[15]，并在放寬網(wǎng)格尺度要求的情況下獲得高精度的數(shù)值解。NND格式本身屬于耗散格式相對較小的TVD格式，本文在保證網(wǎng)格量的前提下開展研究，且網(wǎng)格尺度滿足該格式對于流動拓?fù)湫螒B(tài)的捕捉[16]。張涵信院士團(tuán)隊曾基于該格式開展了復(fù)雜流動問題的模擬，取得了較好的效果[16]。本文將在后續(xù)工作中進(jìn)行高精度格式應(yīng)用的對比研究。

采用迎風(fēng)型NND能獲得較高的間斷分辨率，但也會由于其采用的通量格式而可能帶來非物理解。例如，當(dāng)其通量Jacobi矩陣的特征值很小時，會違反熵條件，產(chǎn)生非物理解，這時需要引入熵修正。本文采用的是Harten熵修正，其表達(dá)式為：

一般0.05＜δ＜0.25，在Ma≤ 0.8時，不采用熵修正，即δ= 0；在Ma＞ 0.8時采用熵修正，δ= 0.1。本文計算中Ma= 0.4，未采用熵修正。

為了保證非定常的時間計算精度，同時又具有較高的計算效率，時間推進(jìn)采用Jameson提出的雙時間步[17-19]方法。對于本文發(fā)展的DES類混合湍流模型和匹配的非定常算法，課題組成員已經(jīng)通過典型算例進(jìn)行了考核和驗證[20-21]，確保本文數(shù)值方法的有效性和可靠性。

2 計算外形、網(wǎng)格和條件

2.1 計算外形

本文基于某圓形橫截面的前體，分別設(shè)計了不同脊形角和高寬比的四個脊形前體模型，表1給出了不同外形的參數(shù)。選用的脊形前體模型是參考Hall的實(shí)驗?zāi)Ｐ蛿?shù)據(jù)[22]，最大半寬bmax為38.1 mm，全長為L= 133.35 mm。圖1給出了實(shí)驗?zāi)Ｐ屯庑螖?shù)據(jù)和本文根據(jù)數(shù)據(jù)點(diǎn)擬合出來的前體橫截面模型的對比。模型寬b、高h(yuǎn)與體長x的關(guān)系式為：

表1 脊形前體計算外形Table 1 Shape parameters of ridged precursors studied

圖1 脊形前體計算模型橫截面與風(fēng)洞實(shí)驗?zāi)Ｐ蛿?shù)據(jù)的比較Fig. 1 Comparison of the cross section between models used in the computation and in the wind tunnel experiment

圖2給出了三個前體的整體外形圖，因為B3和B4外形分別是B1和B2外形上下翻轉(zhuǎn)得到的，所以不再給出整體外形圖以及計算網(wǎng)格示意圖。

圖2 不同前體計算模型Fig. 2 Different models for the forebody

2.2 計算網(wǎng)格

首先考察網(wǎng)格疏密度對計算結(jié)果的影響。選取B1前體，分別生成三套不同規(guī)模的網(wǎng)格，如表2所示，計算網(wǎng)格均為C-H型，因此僅給出G2計算網(wǎng)格，如圖3所示。

表2 不同網(wǎng)格規(guī)模Table 2 Different grid sizes

圖3 B1模型的計算網(wǎng)格Fig. 3 Computation mesh of the B1 model

取來流馬赫數(shù)Ma= 0.4，參考面積Sref= 5 086.35 mm2，參考長度Lref= 38.1 mm，單位雷諾數(shù)Re/Lref= 9.3×106/m，T∞= 300 K，迎角為30°，湍流模型為SA模型。

從圖4中可以看出，三種網(wǎng)格的模擬結(jié)果基本一致，因此，在一定網(wǎng)格規(guī)模以上，網(wǎng)格大小對氣動力系數(shù)的影響很小，但稀網(wǎng)格耗時最短，且收斂快。由于IDDES對網(wǎng)格質(zhì)量要求較高，為了后續(xù)計算中能給出較為準(zhǔn)確的空間渦結(jié)構(gòu)圖，滿足IDDES模擬需求，且綜合考慮計算效率，以下計算均選用G2網(wǎng)格。其他脊形前體模型的網(wǎng)格生成也參照G2網(wǎng)格的規(guī)模。圓形前體的計算網(wǎng)格為O型，網(wǎng)格規(guī)模為800萬。

圖4 氣動力系數(shù)收斂曲線Fig. 4 Convergence curves for aerodynamic coefficients

2.3 計算條件

本文主要比較靜態(tài)不同迎角下脊形前體（B1～B5模型）的氣動特性和空間旋渦流場。取海拔0 km的大氣參數(shù)，參考長度Lref= 38.1 mm，單位雷諾數(shù)Re/Lref= 9.3×106/m，俯仰力矩中心點(diǎn)Xm= 80.01 mm。取來流馬赫數(shù)Ma= 0.4，迎角分別為α= 6°、10°、20°、30°、40°、50°、60°、70°，側(cè)滑角為0°。

3 數(shù)值結(jié)果

3.1 氣動力特性

從圖5中馬赫數(shù)為0.4時不同截面前體氣動力系數(shù)隨迎角的變化可以看出，在所考察的迎角范圍內(nèi)，脊形前體的升力系數(shù)都明顯高于圓形橫截面前體B5的升力系數(shù)。脊形角為7.5°的B1模型的升力系數(shù)極大值出現(xiàn)在40°迎角附近，而脊形角為90°的B2模型和脊形角為180°的B5模型的升力系數(shù)極大值出現(xiàn)在50°迎角附近，這表明在脊形角較小時，產(chǎn)生的渦升力較大。α＞ 20°時，隨著迎角增加，B1模型CD最大，B5模型CD最小，這表明脊形角越小，相應(yīng)的阻力系數(shù)越大。

B1、B2和B5模型的升阻比對比結(jié)果表明，在α＜30°時，脊形角越小，升阻比越大。隨著迎角增大，脊形角對模型升阻比影響減小。對比B1、B3和其分別上下翻轉(zhuǎn)后的B2、B4模型可以看出，B1和B3模型K差別很小，B2和B4模型K差別很小，這表明，在考察的迎角范圍內(nèi)，h/b變化對模型的升阻比影響很小。

對比圖5（d）中俯仰力矩系數(shù)B1、B2和B5可以看出，中到大迎角時，隨著脊形角減小，Cm值增大，脊形前體縱向穩(wěn)定性降低。對比B1、B3和B2、B4可以看出，中小迎角時，上半截面h/b小的脊形前體縱向穩(wěn)定性較差，隨著迎角增加到40° ＜α＜ 60°之間，上半截面h/b較大的脊形前體縱向穩(wěn)定性較差。

3.2 空間流場結(jié)構(gòu)

圖6為不同迎角時用壓力著色的Q等值面。圖中壓力值為無量綱量（當(dāng)?shù)貕簭?qiáng)與來流壓強(qiáng)的比值）。圖7～圖9為大迎角時各軸向位置截面物面壓力系數(shù)周向分布，其采用的是時均處理。圖7為各前體在大迎角時在x/L= 1.0截面處截面流線分布。

圖5 不同前體氣動力系數(shù)隨迎角變化(Ma = 0.4)Fig. 5 Aerodynamic coefficients variation with the angle of attack for different forebodies (Ma = 0.4)

圖6 不同迎角和模型用壓力著色的渦量圖(Q = 0.003，Ma = 0.4；B1～B5：不同計算模型)Fig. 6 Vorticity iso-surfaces coloured by the pressure at different attack angles(Q = 0.003, Ma = 0.4; B1～B5: Different computation models)

從圖6（a）可以看出，具有尖銳脊形前緣（脊形角7.5°）的B1模型，在迎角6°時，自由來流在脊形前緣分離，形成自由剪切層，自由剪切層向上卷起，在前緣附近形成集中的主渦。在迎角10°～30°時，前體主渦Q等值面渦條變得更明顯，背風(fēng)面負(fù)壓增大，主渦強(qiáng)度增加，其位置向上向內(nèi)移動，并相繼形成了二次渦和三次渦。在迎角40°時，主渦破裂位置移到前體底部以前，升力系數(shù)達(dá)到最大值。在迎角50°時，主渦破裂位置已前移到距前體頭部前半體長以內(nèi)，軸向位置x/L= 0.7以后背風(fēng)展向壓力第一第二壓力極小值已消失，分別表征二次渦和三次渦已破裂。在迎角60°～70°時，主渦在前體頭部附近即發(fā)生破裂，背風(fēng)面壓力分布基本呈平坦形狀，表明背風(fēng)區(qū)流場基本處于完全破裂的大分離狀態(tài)。

圖7 迎角40°時不同截面物面壓力系數(shù)分布(Ma = 0.4)Fig. 7 Pressure coefficient distributions of different sections at 40° angle of attack (Ma = 0.4)

圖8 迎角50°時不同截面物面壓力系數(shù)分布(Ma = 0.4)Fig. 8 Pressure coefficient distributions of different sections at 50° angle of attack (Ma = 0.4)

圖9 迎角60°時不同截面物面壓力系數(shù)分布(Ma = 0.4)Fig. 9 Pressure coefficient distributions of different sections at 60° angle of attack (Ma = 0.4)

從圖7～圖9幾個軸向位置周向壓力分布和圖10中軸向x/L= 1.0處截面流線圖可見，B1模型的繞流流場保持較好的對稱性，僅在迎角大于50°后有輕微的非對稱，表明其橫側(cè)向氣動力很小。

從圖6（b）可以看出，對于脊形角為90°的B2模型，在迎角6°時，在較鈍脊形前緣的附近形成集中的主渦，與B1模型相似。在迎角10°～20°時，相對于B1主渦，其Q等值面渦條較細(xì)，緊貼脊形前緣。在迎角30°時，B2模型的主渦渦條有旋擰現(xiàn)象，破裂不明顯，在空間流線可見二次渦和三次渦的位置。在迎角40°時，B2模型的主渦開始在前體底部附近發(fā)生破裂。在迎角50°時，B2模型的主渦破裂位置已前移到距前體頭部，二次渦和三次渦出現(xiàn)破裂，其升力系數(shù)在此迎角達(dá)到最大值，相比B1模型要高。在迎角60°～70°時，B2模型的主渦在前體頭部附近發(fā)生破裂。背風(fēng)面壓力分布基本呈平坦趨勢，表明背風(fēng)區(qū)流場基本處于完全破裂的大分離狀態(tài)，其量值與B1模型的大致相當(dāng)。

從圖7～圖9周向壓力分布和圖11軸向x/L= 1.0處截面流線圖可見，B2模型的繞流流場在迎角小于50°前保持較好的對稱性，在迎角大于60°后有非對稱出現(xiàn)，相比B1模型較明顯，進(jìn)而產(chǎn)生一定橫側(cè)向氣動力。對比B1模型，相同迎角下，脊形角小時前體主渦強(qiáng)度更大，但較強(qiáng)主渦破裂的臨界迎角較小。

圖6（c）對應(yīng)的是B3模型，其與B1具有相同尖銳脊形前緣（脊形角7.5°），但其上下高寬比較大（h/b= 1）。在迎角10°～30°時，B3模型的前體主渦演化與B1模型的類似，三次渦相對不明顯。在迎角40°時，B3模型的主渦破裂位置更靠前，二次渦破裂位置更靠后，其升力系數(shù)也在此迎角達(dá)到最大值。在迎角50°時，B3模型的前體主渦破裂位置和B1模型相近。在迎角60°～70°時，B3模型的前體主渦在前體頭部附近即發(fā)生破裂。背風(fēng)面壓力分布基本呈平坦形狀，其量值與B1模型相比略高。

從圖7～圖9軸向位置周向壓力分布和圖12軸向x/L= 1.0處截面流線圖可見，B3模型的繞流流場也保持較好的對稱性?？梢钥闯觯谙嗤窍?，上下高寬比較小時，前體主渦強(qiáng)度更大，但較強(qiáng)主渦破裂的臨界迎角較小，二次渦破裂延遲。原因是當(dāng)上高寬比大于下高寬比時，上體對左右渦干擾的隔離作用更強(qiáng)，相對地可以推遲主渦破裂。

圖6（d）對應(yīng)的是圓形截面前體B5模型的流場結(jié)構(gòu)演化。在迎角6°時，可見有前體主渦形成的跡象。在迎角10°時，B5模型的前體主渦最弱，繼續(xù)向底部延伸。在迎角20°時，模型的底部附近可見較細(xì)的集中主渦渦條。在迎角30°時，B5模型的前體主渦渦條變得由頭部至底部清晰可見，相對B1模型的較細(xì)，未見破裂。迎角40°時，B5模型的前體主渦渦條有旋擰現(xiàn)象，破裂不明顯。迎角50°時，B5模型的前體主渦開始在前體前半部發(fā)生破裂。迎角60°～70°時，B5模型的前體主渦在前體頭部附近即發(fā)生破裂。背風(fēng)面壓力分布基本呈平坦形狀，其量值與脊形前體模型相比略高。從軸向位置周向壓力分布和軸向x/L= 1.0處截面流線圖13可見，B5模型的繞流流場也保持較好的對稱性。

圖10 前體x/L = 1.0截面處流線分布(B1, Ma = 0.4)Fig. 10 Streamline distributions at the section x/L = 1.0 (B1，Ma = 0.4)

圖11 前體x/L = 1.0截面處流線分布(B2, Ma = 0.4)Fig. 11 Streamline distributions at the section x/L = 1.0 (B2, Ma = 0.4)

圖12 前體x/L=1.0截面處流線分布(B3，Ma = 0.4)Fig. 12 Streamline distributions at the section x/L = 1.0 (B3, Ma = 0.4)

圖13 前體x/L=1.0截面處流線分布(B5，Ma = 0.4)Fig. 13 Streamline distributions at the section x/L = 1.0 (B5, Ma = 0.4)

通過上述不同橫截面前體在不同來流參數(shù)下氣動特性的對比，以及氣動力特性規(guī)律與空間渦結(jié)構(gòu)關(guān)聯(lián)的分析，主要有以下幾點(diǎn)結(jié)論：

1）迎角變化時，較小脊形角外形的前體主渦強(qiáng)度較大，前體渦產(chǎn)生和破裂的臨界迎角較小。在迎角較小時，脊形角越小，較強(qiáng)的前體渦產(chǎn)生的升力越大，升阻比越高；隨著迎角進(jìn)一步增加，不同脊形角的前體渦均破裂，升阻特性趨于一致。在小迎角時，脊形角越小，前體渦升力越大，誘導(dǎo)抬頭力矩，因而縱向穩(wěn)定性變差；在大迎角時，前體渦完全破裂，縱向穩(wěn)定性改善。與圓形橫截面前體相比，脊形前體主渦盡管破裂提前，但相繼產(chǎn)生較強(qiáng)的二次渦和三次渦，使得主渦破裂后升力減小平緩，失速特性明顯優(yōu)于圓形截面前體。

2）迎角變化時，對于相同脊形角的前體，上半截面高寬比（h/b）越小，前體主渦強(qiáng)度越大，前體渦破裂臨界迎角減小。在迎角較小時，上半截面高寬比小的脊形前體的前體渦雖較強(qiáng)，但二次渦較弱，因而渦升力僅略大，而迎風(fēng)阻力面積也增大，阻力略增，故升阻比相當(dāng)。隨著迎角進(jìn)一步增加，上半截面高寬比小的脊形前體，由于左右渦的相互干擾更強(qiáng)，主渦提前破裂，但二次渦破裂延遲，因而升力系數(shù)減小，但渦誘導(dǎo)阻力也相應(yīng)減小，故升阻比仍然相當(dāng)。此外，迎角較小時，上半截面h/b小的脊形前體，由于渦升力更大，誘導(dǎo)的抬頭力矩也更大，因而縱向穩(wěn)定性變差。隨著迎角增加，上半截面h/b小的脊形前體渦提前破裂，渦升力下降，誘導(dǎo)的抬頭力矩也減小，故縱向穩(wěn)定性有所改善。

4 總 結(jié)

本文針對脊形前體飛行器大迎角湍流大分離流動計算的困難，采用DES類混合湍流模型，以及與之匹配的非定常算法，開展脊形前體在不同迎角下的靜態(tài)繞流數(shù)值模擬研究，對五種不同橫截面形狀的前體，數(shù)值研究了不同迎角下背風(fēng)渦流場結(jié)構(gòu)和氣動力特性，分析了不同脊形角和上下半截面h/b不同時對前體氣動力系數(shù)和空間渦結(jié)構(gòu)的影響，主要得到如下結(jié)論：

1）隨著迎角的增大，前體主渦經(jīng)歷了從形成到增強(qiáng)再到破裂的演化過程，氣動升力系數(shù)也隨之先升高后減小。渦破裂位置和最大升力系數(shù)對應(yīng)的迎角與脊形前體形狀有關(guān)。

2）脊形角越小，前體渦越強(qiáng)，渦升力越大，產(chǎn)生前體主渦、二次渦、三次渦的臨界迎角就越小，前體渦產(chǎn)生和破裂的臨界迎角越小。

3）脊形前體不同橫截面形狀對渦演化的影響規(guī)律表現(xiàn)為：相同脊形角的前體，上下高寬比越小，前體主渦強(qiáng)度越大，前體渦破裂臨界迎角越小。

4）與圓形橫截面前體相比，脊形前體有著較平緩的升力減小趨勢，即具備更好的失速特性。并且，脊形前體渦形態(tài)較為穩(wěn)定，對陣風(fēng)擾動的抵抗能力強(qiáng)，相應(yīng)地，橫航向氣動特性會更優(yōu)良。