不同迎角下脊形前體繞流數值模擬研究

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

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

0 引 言

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

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

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

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

1 數值模擬方法

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

1.1 湍流方程和模型

一般曲線坐標系下，湍流方程可以寫作統一的形式：

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

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

1.2 數值離散

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

式中，下標i代表網格單元中心位置，i±1/2代表單元界面位置。

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

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

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

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

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

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

2 計算外形、網格和條件

2.1 計算外形

本文基于某圓形橫截面的前體，分別設計了不同脊形角和高寬比的四個脊形前體模型，表1給出了不同外形的參數。選用的脊形前體模型是參考Hall的實驗模型數據[22]，最大半寬bmax為38.1 mm，全長為L= 133.35 mm。圖1給出了實驗模型外形數據和本文根據數據點擬合出來的前體橫截面模型的對比。模型寬b、高h與體長x的關系式為：

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

圖1 脊形前體計算模型橫截面與風洞實驗模型數據的比較Fig. 1 Comparison of the cross section between models used in the computation and in the wind tunnel experiment

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

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

2.2 計算網格

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

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

圖3 B1模型的計算網格Fig. 3 Computation mesh of the B1 model

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

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

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

2.3 計算條件

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

3 數值結果

3.1 氣動力特性

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

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

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

3.2 空間流場結構

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

圖5 不同前體氣動力系數隨迎角變化(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等值面渦條變得更明顯，背風面負壓增大，主渦強度增加，其位置向上向內移動，并相繼形成了二次渦和三次渦。在迎角40°時，主渦破裂位置移到前體底部以前，升力系數達到最大值。在迎角50°時，主渦破裂位置已前移到距前體頭部前半體長以內，軸向位置x/L= 0.7以后背風展向壓力第一第二壓力極小值已消失，分別表征二次渦和三次渦已破裂。在迎角60°～70°時，主渦在前體頭部附近即發(fā)生破裂，背風面壓力分布基本呈平坦形狀，表明背風區(qū)流場基本處于完全破裂的大分離狀態(tài)。

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

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

圖9 迎角60°時不同截面物面壓力系數分布(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°后有輕微的非對稱，表明其橫側向氣動力很小。

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

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

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

從圖7～圖9軸向位置周向壓力分布和圖12軸向x/L= 1.0處截面流線圖可見，B3模型的繞流流場也保持較好的對稱性。可以看出，在相同迎角下，上下高寬比較小時，前體主渦強度更大，但較強主渦破裂的臨界迎角較小，二次渦破裂延遲。原因是當上高寬比大于下高寬比時，上體對左右渦干擾的隔離作用更強，相對地可以推遲主渦破裂。

圖6（d）對應的是圓形截面前體B5模型的流場結構演化。在迎角6°時，可見有前體主渦形成的跡象。在迎角10°時，B5模型的前體主渦最弱，繼續(xù)向底部延伸。在迎角20°時，模型的底部附近可見較細的集中主渦渦條。在迎角30°時，B5模型的前體主渦渦條變得由頭部至底部清晰可見，相對B1模型的較細，未見破裂。迎角40°時，B5模型的前體主渦渦條有旋擰現象，破裂不明顯。迎角50°時，B5模型的前體主渦開始在前體前半部發(fā)生破裂。迎角60°～70°時，B5模型的前體主渦在前體頭部附近即發(fā)生破裂。背風面壓力分布基本呈平坦形狀，其量值與脊形前體模型相比略高。從軸向位置周向壓力分布和軸向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)

通過上述不同橫截面前體在不同來流參數下氣動特性的對比，以及氣動力特性規(guī)律與空間渦結構關聯的分析，主要有以下幾點結論：

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

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

4 總 結

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

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

2）脊形角越小，前體渦越強，渦升力越大，產生前體主渦、二次渦、三次渦的臨界迎角就越小，前體渦產生和破裂的臨界迎角越小。

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

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