龍海斌，吳裕平

(中國直升機設(shè)計研究所，江西 景德鎮(zhèn) 333001)

0 引言

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

圖1 ROBIN機身外形與表面測壓孔截面位置圖[33]

1 網(wǎng)格類型

在機身氣動特性數(shù)值模擬之前，首先要選擇機身計算域網(wǎng)格劃分類型。目前常用的網(wǎng)格類型有結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格和混合網(wǎng)格等。其中結(jié)構(gòu)網(wǎng)格的質(zhì)量比較高，劃分速度也比較快，在部分情況下計算結(jié)果的準確度比較高。但是，結(jié)構(gòu)網(wǎng)格通常只適用于外形比較簡單的氣動特性數(shù)值模擬。非結(jié)構(gòu)網(wǎng)格對復(fù)雜外形的適應(yīng)能力比較強，在對復(fù)雜外形流體域劃分網(wǎng)格時，非結(jié)構(gòu)網(wǎng)格的數(shù)量相對較少。

由于ROBIN機身外形相對比較簡單，因此部分研究人員在CFD計算過程中采用結(jié)構(gòu)化網(wǎng)格，而且結(jié)構(gòu)網(wǎng)格劃分得比較細密，如圖2所示。

圖2 結(jié)構(gòu)網(wǎng)格示意圖[5]

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

圖3 非結(jié)構(gòu)四面體網(wǎng)格示意圖

2 網(wǎng)格數(shù)量

在機身氣動特性數(shù)值模擬過程中，網(wǎng)格數(shù)量不僅影響數(shù)值模擬結(jié)果的準確度，而且與整個數(shù)值模擬耗費的時間密切相關(guān)。同時，網(wǎng)格數(shù)量與網(wǎng)格劃分的時間也緊密相關(guān)。在型號研制過程中還要考慮整個項目周期的要求等，在部分情況下還要考慮計算過程中的內(nèi)存占用，因此需要限制流體域劃分的網(wǎng)格數(shù)量。

網(wǎng)格劃分比較稀疏時，可能造成部分機身外形失真，保形能力比較差，導(dǎo)致數(shù)值模擬結(jié)果有誤差，同時部分流場細節(jié)難以模擬出來。而網(wǎng)格劃分比較細密時，有可能導(dǎo)致部分網(wǎng)格的質(zhì)量比較差，比如在非四面體網(wǎng)格劃分過程中產(chǎn)生“金字塔”形網(wǎng)格，計算過程中難以收斂。網(wǎng)格數(shù)量多時，不僅網(wǎng)格劃分耗費時間多，而且計算過程中耗費的時間也長，占用內(nèi)存等計算資源也比較多。

在ROBIN機身數(shù)值模擬過程中，劃分的網(wǎng)格數(shù)量大多數(shù)在百萬量級，最少的為40.9萬四面體非結(jié)構(gòu)網(wǎng)格，最多的達到1064萬六面體結(jié)構(gòu)網(wǎng)格。在公開發(fā)表的19個數(shù)值模擬過程中，網(wǎng)格或網(wǎng)格點數(shù)量在各數(shù)量級的分布如表1所示。從表中可以看出網(wǎng)格或網(wǎng)格點數(shù)量在100萬～1000萬之間的占57.9%。部分網(wǎng)格如圖4和圖5所示。

表1 網(wǎng)格或網(wǎng)格點的數(shù)量級分布

圖4 46.9萬非結(jié)構(gòu)網(wǎng)格示意圖[12]

圖5 1064萬結(jié)構(gòu)網(wǎng)格示意圖[13]

在-5°攻角時，采用不同的網(wǎng)格數(shù)量數(shù)值模擬得到的凸臺結(jié)束位置(x/R=1.0008)的壓力系數(shù)與風洞試驗結(jié)果的對比如圖6所示。

圖6 網(wǎng)格數(shù)量變化時數(shù)值模擬結(jié)果對比圖

從圖中可以看出，在一定數(shù)量范圍內(nèi)，網(wǎng)格數(shù)量的變化對數(shù)值模擬結(jié)果的準確度影響比較小。

3 求解方法

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

表2 各湍流模型應(yīng)用情況表

圖7 不同湍流模型模擬結(jié)果對比圖[20]

部分研究人員采用尺度自適應(yīng)模擬(SAS)方法對ROBIN機身模型的氣動特性進行了數(shù)值模擬[22]。該方法通過引入第二長度尺度，可以實現(xiàn)雷諾平均(RANS)和大渦模擬(LES)方法的混合數(shù)值模擬。由于該方法提出的時間比較短，目前應(yīng)用比較少。還有個別研究人員采用渦量輸運的方法對機身模型的氣動特性進行數(shù)值模擬[23]。

4 結(jié)果分析

目前檢驗氣動特性數(shù)值模擬結(jié)果準確性的主要方法是將其與風洞試驗結(jié)果、理論解析值或其他公開發(fā)表的數(shù)值模擬結(jié)果等進行對比分析。ROBIN機身模型風洞試驗中給出了機身表面的壓力系數(shù)，因此研究人員將自己計算得到的機身表面壓力系數(shù)結(jié)果與風洞試驗結(jié)果進行了對比分析。大多數(shù)數(shù)值模擬過程中的機身模型長度為2m；由于直升機在大速度前飛為低頭姿態(tài)，因而計算過程中選取的攻角多為0°或-5°，側(cè)滑角為0°；選取的截面位置通常為x/L=0.0517、0.3074、0.4669、0.6003、1.0008,1.1620等。從對比分析的情況來看，大多數(shù)數(shù)值模擬結(jié)果與風洞試驗值比較接近，在靠近凸臺和腹部支撐附近的壓力系數(shù)計算值與風洞試驗結(jié)果有一定的差別。這是由于在數(shù)值模擬過程中沒有考慮腹部支撐機構(gòu)，同時這些區(qū)域存在一定的流動分離，而目前數(shù)值模擬方法在模擬流動分離時存在一定的困難。

對機身模型進行數(shù)值模擬，除了得到機身表面壓力系數(shù)，還能得到機身附近的速度、流線和渦量變化等流場細節(jié)。這些數(shù)據(jù)對分析機身周圍流動的機理，直升機機身減阻增升設(shè)計等都非常有幫助。采用面元法和求解N-S方法得到機身表面流線圖。從圖8中可以看出兩種方法得到的前機身和中機身流線圖基本上一致，但是在尾梁后下端位置的流線有一定的差別。不同計算方法得到的渦量如圖9所示，從圖中可以看出兩者的主要差別在凸臺尾流區(qū)的渦量，說明目前數(shù)值方法對尾流區(qū)的復(fù)雜流動模擬還存在一定的誤差。

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

圖9 Q準則得到的瞬時渦結(jié)構(gòu)等值圖[22]

早期ROBIN機身模型風洞試驗給出的是機身表面壓力系數(shù)，而在實際工程應(yīng)用中更關(guān)注的是數(shù)值模擬方法計算機身氣動力和力矩結(jié)果的準確度和精度。南京航空航天大學唐韜等[33]在某開口回流式風洞中對ROBIN機身模型的力和力矩進行了測量。試驗過程中的來流速度分別為9m/s、18m/s和27m/s。后續(xù)在ROBIN機身數(shù)值模擬技術(shù)的相關(guān)研究中可參考此次風洞試驗的結(jié)果。

5 總結(jié)與討論

通過對ROBIN機身模型數(shù)值模擬過程中的網(wǎng)格類型選取、網(wǎng)格數(shù)量控制、求解方法選擇和與試驗對比分析等進行總結(jié)與分析，可得出如下結(jié)論：

1)機身模型數(shù)值模擬過程中有比較多的參數(shù)和變量，如網(wǎng)格類型、網(wǎng)格數(shù)量、湍流模型等。改變其中的一個參數(shù)或變量或許對提高數(shù)值模擬結(jié)果的準確度有一定的效果。但是如果需要大幅度提高數(shù)值模擬結(jié)果的準確度和精度，則需要綜合考慮多方面因素的影響。

2)在機身模型數(shù)值模擬過程中，網(wǎng)格類型、網(wǎng)格數(shù)量和湍流模型等的選取不僅要考慮數(shù)值模擬結(jié)果的準確度等，同時還需要綜合考慮可用的計算資源，計算時間的限制,計算方法的連續(xù)性等。在工程應(yīng)用領(lǐng)域，對ROBIN機身模型劃分100萬左右的非結(jié)構(gòu)網(wǎng)格可滿足相應(yīng)的需求。在數(shù)值模擬過程中網(wǎng)格類型、網(wǎng)格數(shù)量和湍流模型等可考慮配套使用。

3)未來需要開發(fā)多種類型的直升機標準機身模型，以支持不同類型直升機的發(fā)展。在機身模型風洞試驗過程中，不僅測量機身表面壓力系數(shù)，還需要測量機身的阻力、升力、俯仰力矩和偏航力矩等力和力矩系數(shù)，同時還可以考慮測量機身表面附近的速度分布、尾流區(qū)域的渦系等，以方便后續(xù)的數(shù)值模擬研究工作的開展。

4)在數(shù)值模擬研究中考慮制定相應(yīng)的標準數(shù)值模擬方法。在直升機型號研制過程中，由于項目周期和可用計算資源的限制，在一定程度上可以說機身氣動特性數(shù)值模擬的速度越快越好。這就需要在標準機身模型數(shù)值模擬研究中找出兼顧效率和準確度的標準方法。