龍海斌,吳裕平
(中國直升機設計研究所,江西 景德鎮 333001)
飛行器標準機身模型通常用于風洞測量系統校驗、機身氣動特性CFD計算驗證與確認等。在直升機機身風洞試驗和數值模擬計算發展過程中,相繼發展了ROBIN、ROBIN mod7、HELIFUSE和NUAA等標準機身模型。ROBIN(Rotor Body Interaction,簡稱ROBIN)機身模型的三維外形可用超橢圓方程進行描述,主要被美國NASA等機構用于旋翼/機身干擾研究、直升機機身氣動特性數值模擬驗證與確認等工作。1979年美國針對ROBIN機身模型進行了無旋翼狀態(光機身)風洞測量試驗,得到了部分表面位置的壓力系數[1]。機身的長度為旋翼半徑R的2倍。壓力系數的測量位置包括機身開始端截面、向凸臺過渡段截面、凸臺開始截面、凸臺中段、凸臺結束位置以及尾梁截面等典型位置,如圖1所示。測量過程中包含-10°、-5°、0°、5°攻角和-15°、-5°、0°、5°側滑角狀態。由于有風洞試驗結果可以對比分析,而且機身外形比較簡單,全球很多科研人員對ROBIN機身模型的繞流進行了數值模擬。本文對近年來ROBIN機身模型的數值模擬技術進行了梳理與分析,包括網格劃分之前的網格類型的選擇,網格劃分過程中的數量控制,求解過程中的計算方法選取,數值模擬結果分析等[2-32]。

圖1 ROBIN機身外形與表面測壓孔截面位置圖[33]
在機身氣動特性數值模擬之前,首先要選擇機身計算域網格劃分類型。目前常用的網格類型有結構網格、非結構網格和混合網格等。其中結構網格的質量比較高,劃分速度也比較快,在部分情況下計算結果的準確度比較高。但是,結構網格通常只適用于外形比較簡單的氣動特性數值模擬。非結構網格對復雜外形的適應能力比較強,在對復雜外形流體域劃分網格時,非結構網格的數量相對較少。
由于ROBIN機身外形相對比較簡單,因此部分研究人員在CFD計算過程中采用結構化網格,而且結構網格劃分得比較細密,如圖2所示。

圖2 結構網格示意圖[5]
在統計的29例ROBIN機身模型數值模擬過程中,有20例采用非結構網格,占總數的69%,說明大部分研究人員還是傾向于應用非結構網格進行機身氣動特性數值模擬。主要原因有兩個:一是ROBIN機身數值模擬通常作為數值模擬方法的驗證算例,后續該數值模擬方法還要用于更復雜的機身繞流、外掛物流場甚至旋翼/機身干擾流場等計算,為了保證數值模擬方法驗證的有效性,需采用一致的網格劃分方法。二是隨著非結構網格的發展與應用,非結構網格的劃分過程比較簡單,耗費的時間也比較少,而且數值模擬結果的準確度和精度越來越高,因此非結構網格的應用越來越廣泛。典型的非結構四面體網格如圖3所示[6]。

圖3 非結構四面體網格示意圖
在機身氣動特性數值模擬過程中,網格數量不僅影響數值模擬結果的準確度,而且與整個數值模擬耗費的時間密切相關。同時,網格數量與網格劃分的時間也緊密相關。在型號研制過程中還要考慮整個項目周期的要求等,在部分情況下還要考慮計算過程中的內存占用,因此需要限制流體域劃分的網格數量。
網格劃分比較稀疏時,可能造成部分機身外形失真,保形能力比較差,導致數值模擬結果有誤差,同時部分流場細節難以模擬出來。而網格劃分比較細密時,有可能導致部分網格的質量比較差,比如在非四面體網格劃分過程中產生“金字塔”形網格,計算過程中難以收斂。網格數量多時,不僅網格劃分耗費時間多,而且計算過程中耗費的時間也長,占用內存等計算資源也比較多。
在ROBIN機身數值模擬過程中,劃分的網格數量大多數在百萬量級,最少的為40.9萬四面體非結構網格,最多的達到1064萬六面體結構網格。在公開發表的19個數值模擬過程中,網格或網格點數量在各數量級的分布如表1所示。從表中可以看出網格或網格點數量在100萬~1000萬之間的占57.9%。部分網格如圖4和圖5所示。

表1 網格或網格點的數量級分布

圖4 46.9萬非結構網格示意圖[12]

圖5 1064萬結構網格示意圖[13]
在-5°攻角時,采用不同的網格數量數值模擬得到的凸臺結束位置(x/R=1.0008)的壓力系數與風洞試驗結果的對比如圖6所示。

圖6 網格數量變化時數值模擬結果對比圖
從圖中可以看出,在一定數量范圍內,網格數量的變化對數值模擬結果的準確度影響比較小。
在ROBIN機身模型數值模擬過程中,主要采用求解Eluer方程和求解N-S方程兩種方法。在早期,部分研究人員采用面元法[17,18]求解,這種方法通過布置流動的奇點來求解氣動問題。面元法的優點是計算速度非常快,網格劃分簡單,但是數值模擬結果的精度比較低,現在已經很少使用。求解Euler方程的方法在數值模擬過程中沒有考慮流體的粘性,因此該方法計算速度相對比較快,計算過程中占用的計算資源也比較少。因此,后續對更復雜的組合機身外形、加裝外掛物和旋翼/機身干擾流動等進行數值模擬時更方便和快速。由于沒有考慮空氣的粘性,因此對部分流動復雜區域的模擬能力很弱。目前大多數人員采用求解N-S方程的方法來對機身模型氣動特性進行數值模擬。在統計的31個算例中,有19個采用求解N-S方程的方法,占61.3%。目前求解N-S方程的方法主要有雷諾平均(RANS)方法、大渦模擬(LES)方法和直接數值模擬(DNS)方法等。其中雷諾平均(RANS)方法對網格數量的要求最低,消耗的計算資源也比較少,數值模擬所花費的時間也最少,因此應用最廣泛。在運用雷諾平均(RANS)方法對機身模型氣動特性進行數值模擬時,需要引入湍流模型來使得方程封閉。常用的湍流模型有B-L模型、S-A模型和k-ω模型等,其中S-A湍流模型應用最多,占到2/3以上。采用各湍流模型的數量和百分比如表2所示。采用不同湍流模型得到的數值模擬與風洞試驗結果如圖7所示。

表2 各湍流模型應用情況表

圖7 不同湍流模型模擬結果對比圖[20]
部分研究人員采用尺度自適應模擬(SAS)方法對ROBIN機身模型的氣動特性進行了數值模擬[22]。該方法通過引入第二長度尺度,可以實現雷諾平均(RANS)和大渦模擬(LES)方法的混合數值模擬。由于該方法提出的時間比較短,目前應用比較少。還有個別研究人員采用渦量輸運的方法對機身模型的氣動特性進行數值模擬[23]。
目前檢驗氣動特性數值模擬結果準確性的主要方法是將其與風洞試驗結果、理論解析值或其他公開發表的數值模擬結果等進行對比分析。ROBIN機身模型風洞試驗中給出了機身表面的壓力系數,因此研究人員將自己計算得到的機身表面壓力系數結果與風洞試驗結果進行了對比分析。大多數數值模擬過程中的機身模型長度為2m;由于直升機在大速度前飛為低頭姿態,因而計算過程中選取的攻角多為0°或-5°,側滑角為0°;選取的截面位置通常為x/L=0.0517、0.3074、0.4669、0.6003、1.0008,1.1620等。從對比分析的情況來看,大多數數值模擬結果與風洞試驗值比較接近,在靠近凸臺和腹部支撐附近的壓力系數計算值與風洞試驗結果有一定的差別。這是由于在數值模擬過程中沒有考慮腹部支撐機構,同時這些區域存在一定的流動分離,而目前數值模擬方法在模擬流動分離時存在一定的困難。
對機身模型進行數值模擬,除了得到機身表面壓力系數,還能得到機身附近的速度、流線和渦量變化等流場細節。這些數據對分析機身周圍流動的機理,直升機機身減阻增升設計等都非常有幫助。采用面元法和求解N-S方法得到機身表面流線圖。從圖8中可以看出兩種方法得到的前機身和中機身流線圖基本上一致,但是在尾梁后下端位置的流線有一定的差別。不同計算方法得到的渦量如圖9所示,從圖中可以看出兩者的主要差別在凸臺尾流區的渦量,說明目前數值方法對尾流區的復雜流動模擬還存在一定的誤差。

圖8 -10°攻角時不同方法得到的機身表面流線圖[24]

圖9 Q準則得到的瞬時渦結構等值圖[22]
早期ROBIN機身模型風洞試驗給出的是機身表面壓力系數,而在實際工程應用中更關注的是數值模擬方法計算機身氣動力和力矩結果的準確度和精度。南京航空航天大學唐韜等[33]在某開口回流式風洞中對ROBIN機身模型的力和力矩進行了測量。試驗過程中的來流速度分別為9m/s、18m/s和27m/s。后續在ROBIN機身數值模擬技術的相關研究中可參考此次風洞試驗的結果。
通過對ROBIN機身模型數值模擬過程中的網格類型選取、網格數量控制、求解方法選擇和與試驗對比分析等進行總結與分析,可得出如下結論:
1)機身模型數值模擬過程中有比較多的參數和變量,如網格類型、網格數量、湍流模型等。改變其中的一個參數或變量或許對提高數值模擬結果的準確度有一定的效果。但是如果需要大幅度提高數值模擬結果的準確度和精度,則需要綜合考慮多方面因素的影響。
2)在機身模型數值模擬過程中,網格類型、網格數量和湍流模型等的選取不僅要考慮數值模擬結果的準確度等,同時還需要綜合考慮可用的計算資源,計算時間的限制,計算方法的連續性等。在工程應用領域,對ROBIN機身模型劃分100萬左右的非結構網格可滿足相應的需求。在數值模擬過程中網格類型、網格數量和湍流模型等可考慮配套使用。
3)未來需要開發多種類型的直升機標準機身模型,以支持不同類型直升機的發展。在機身模型風洞試驗過程中,不僅測量機身表面壓力系數,還需要測量機身的阻力、升力、俯仰力矩和偏航力矩等力和力矩系數,同時還可以考慮測量機身表面附近的速度分布、尾流區域的渦系等,以方便后續的數值模擬研究工作的開展。
4)在數值模擬研究中考慮制定相應的標準數值模擬方法。在直升機型號研制過程中,由于項目周期和可用計算資源的限制,在一定程度上可以說機身氣動特性數值模擬的速度越快越好。這就需要在標準機身模型數值模擬研究中找出兼顧效率和準確度的標準方法。