DPW4翼／身／平尾組合體的數值模擬

王運濤，張書俊，孟德虹

（中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000）

0 引 言

隨著大型網格生成軟件、大規模并行計算技術和后置處理軟件的發展，CFD幾乎可以模擬所有高度復雜飛行器外形的繞流流場［1］，CFD已經與風洞試驗一道，成為飛行器氣動設計工作者最重要的研究手段。飛行器設計者使用CFD工具面臨的首要問題是軟件的計算精度和效率問題。多重網格技術、低速預處理技術和并行計算技術的綜合應用，極大地提高了CFD軟件的計算效率，但CFD的驗證與確認問題依然是當前計算空氣動力學研究的熱點［2］。

為研究CFD模擬的精準度問題，國際上先后組織了許多專題研討會，其中比較具有代表性的如歐洲的計算空氣動力學研究項目 ECARP［3］（European Computational Aerodynamics Research Project）和AIAA的DPW（Drag Prediction Workshop）系列會議。為了開展典型運輸機構型的阻力預測精度研究，AIAA阻力計算的工作小組先后在2001年6月、2003年6月、2006年6月召開了三次阻力計算的工作會議，選擇的研究構型包括DLR－F4翼身組合體構型、DLR－F6翼／身／架／艙復雜組合體構型、DLR－F6翼身組合體及其修型構型等。該系列會議通過規定統一的計算狀態、統一的計算網格和精細的試驗結果，對世界范圍內各家研究機構提供的CFD軟件計算結果進行評估，明確CFD的發展現狀和下一步的努力方向。

AIAA的第四屆阻力計算的工作會議DPW 4（AIAA 4thDrag Prediction Workshop）于2009年6月召開。本次會議的基本宗旨與前三次會議相同，在總結前三次會議經驗的基礎上，在組織方式上做了一些調整。第一是重新設計了研究模型，即共同研究模型 CRM（Common Research Model）［4］，采用了現代運輸機的翼／身／平尾組合體構型，設計馬赫數為0.85；第二是事先不提供任何風洞試驗結果，即采用盲評的方式組織此次專題研討活動。

本文基于TRIP3.0（TRIsonic Platform V3.0）軟件和多塊對接網格技術，采用粗、中、細三套結構網格，選擇SST兩方程湍流模型和MUSCL型的ROE格式，通過求解任意坐標系下的雷諾平均NS方程（RANS），數值模擬了DPW－IV翼／身／平尾組合體構型的氣動特性。數值模擬內容包括：網格收斂性研究、下洗影響和阻力發散馬赫數計算。在與National Transonic Facility（NFT）的試驗結果對比基礎上，本文詳細研究了網格密度、來流馬赫數和平尾偏角對總體氣動特性和壓力分布的影響。通過研究發現TRIP軟件能夠較好地模擬該構型的氣動特性；網格密度對壓差阻力的影響較大，對摩擦阻力的影響較?。痪W格密度對壓力分布的影響不明顯。

1 TRIP3.0軟件的簡介

TRIP是中國空氣動力研究與發展中心自行研發的大型通用CFD軟件，已經成功推廣到國內19家氣動研究單位，極大地促進了國內CFD應用水平的提高。TRIP軟件從0.0版本、1.0版本、2.0版本，發展到現在的3.0版本。TRIP3.0版本在兼容2.0版本［5］所有功能的基礎上，與本項研究相關的幾個方面的改進包括：

·改進了多重網格收斂加速技術，克服了多重網格技術不能在并行計算中使用的缺點；

·對算法進行了優化，減少了冗余運算或重復運算，每一步的計算效率同比提高30%以上；

·提高了SA一方程湍流模型和SST兩方程湍流模型的求解精度，增強了軟件的魯棒性；

·增加了雙精度計算選項，進一步提高了數值模擬的精度。

2 計算網格

本文采用的計算網格由Boeing公司提供［6］，網格結構為多塊對接網格（1－to－1），網格分為粗網格（Coarse）、中等網格（Medium）和細網格（Fine）三種。這三套網格的y＋分別約為1、1／3和4／9，壁面第一層網格分別不大于0.0375mm、0.025mm和0.0166mm。全機網格的詳細信息如表1所示，網格拓撲和表面網格如圖1所示。

計算采用的平均氣動弦長MAC＝7.00532m，參考面積Sref＝383.6895552m2，力矩參考點：Xm＝33.67786m，Ym＝4.51993m，Zm＝－11.90625m。

圖1 翼／身／尾組合體網格拓撲以及表面網格（中等）Fig.1 Surface grid for DPW4wing／body／tail（Medium）

3 計算結果和分析

本文的數值模擬采用TRIP3.0軟件，選擇二階精度的通量差分（FDS）類型的MUSCL差分格式［7］，Menter’sk－ωSST 兩方程湍流模型［8］，綜合采用局部時間步長、三重網格加速收斂和大規模并行計算技術，數值模擬了兩種工況的流場，其中 WBH0、WBHm2和 WBHp2的平尾迎角分別為0°、－2°和2°，WBHNo為CRM不帶平尾的模型。

3.1 CASE 1計算結果與分析

3.1.1 網格收斂性研究（CASE1a）

網格收斂性研究的計算構型為 WBH0，平尾偏角為0°，數值模擬采用粗網格、中等網格、細網格，計算狀態如下：

圖2給出了WBH0組合體采用粗網格、中等網格和密網格的阻力系數、摩阻系數和壓差阻力系數的計算結果，其中橫坐標為網格點數的－2／3次冪。表2給出了CL＝0.500（±0.001）翼身組合體對應的迎角、阻力系數及其分量，可以看出，采用上述三套網格，在固定升力系數下，迎角、總阻力系數、壓差阻力系數、摩擦阻力系數及俯仰力矩系數均是單調變化的。從粗網格到密網格，迎角從2.305°增加到2.32°，總阻力系數從0.02722減少到0.02682，壓差阻力系數從0.01218增加到0.01232，摩擦阻力系數從0.01504減少到0.01450，俯仰力矩系數從－0.05196增加到－0.04357。圖中同時給出了文獻［9］的計算結果。文獻［9］中采用了非結構網格，粗、中、細三種網格規模分別為410萬、1170萬和3410萬，空間離散格式采用二階中心格式，湍流模型采用SST兩方程湍流模型。

根據Richardson外推公式，本文得到的固定升力系數下的迎角、阻力系數、力矩系數的網格無關性結果分別為：2.340、0.02701和－0.04025，與文獻［10］的對比如表3所示。文獻［10］提供了19家單位29個計算結果的算術平均值及標準方差，本文的計算結果均落在文獻［10］的計算誤差帶內。文獻［10］同時給出了試驗結果的平均值及誤差，與相應的試驗結果的平均值相比較，阻力的計算誤差接近0.0008，迎角相差0.687°，力矩出現了符號的差異。導致計算結果與風洞試驗結果的差異的主要原因首先應該是風洞試驗模型添加了轉捩帶，而數值模擬采用了全轉捩的模擬方式；其次，是風洞試驗中采用了垂直尾部支撐的方式，試驗模型尾部的支撐形式對氣動特性有明顯的影響，對力矩特性的影響尤其顯著。具體的原因需要通過更加精細的試驗和數值模擬進行分析。

圖2 DPW4翼／身／平尾組合體的網格收斂性研究Fig.2 Grid convergence study for DPW4wing／body／tail

表2 WBH0翼／身／平尾組合體的網格收斂性研究Table 2 Grid convergence study for WBH0wing／body／tail

表3 WBH0組合體外推結果與文獻的比較Table 3 WBH0data extrapolated to continuum

圖3給出了平尾迎角為0°的翼／身／尾組合體分別采用粗網格、中等網格和密網格，CL＝0.500±0.001時，平尾附近的壓力分布及表面流線云圖?？梢钥闯觯煌木W格密度下均存在一定的分離區，三套網格給出的流線形態基本一致，密網格得到的平尾后部的分離區要小一些。圖4給出了其在13.06%、28.28%、39.71%、50.24%、60.28%、72.68%、84.56%、95%等8個典型站位上的壓力分布比較?？梢钥闯觯跈C翼典型站位上，粗網格、中等網格和密網格的壓力分布幾乎完全重合。

圖3 壓力分布及表面流線云圖（CL＝0.500±0.001）Fig.3 DPW4wing／body／tail surface streamline

圖4 典型站位壓力分布（CL＝0.500）Fig.4 Cpdistributions for typical wing stations

3.1.2 配平影響（CASE1b）

計算狀態如下：

·M＝0.85

·α＝0°，1.0°，1.5°，2°，2.5°，3.0°，4.0°

·平尾迎角：－2°，0°，2°，無平尾

·Re＝5.0×106（基于MAC＝7.00532m）

基于中等規模的網格，進一步研究比較了M＝0.85來流狀態下，平尾迎角分別為－2°、0°、2°和無平尾翼身組合體的氣動特性。根據平尾迎角為－2°、0°和2°時，Cm～CL曲線，可以插值得到不同升力系數下平尾的配平迎角，如圖5所示。升力系數小于0.34時，配平迎角為正值，絕對值隨著升力系數的增加而逐漸減小；升力系數大于0.34時，配平迎角為負值，絕對值隨著升力系數的增加而增加，在升力系數0.6時，配平迎角絕對值達到最大；升力系數進一步增加，配平迎角絕對值逐漸減小。

圖5 下洗研究Fig.5 Downwash study

根據配平迎角隨升力系數的變化曲線和阻力系數隨升力系數的變化曲線，進一步得到配平后阻力系數隨升力系數的變化曲線，翼身組合體與配平后的阻力差量變化趨勢見圖6，升力系數小于0.42時，配平后引起的阻力增量基本不變，量值在20個阻力單位（1阻力單位＝0.0001）；升力系數大于0.42后，配平后引起的阻力增量隨升力系數的增加而基本線性增加，在升力系數0.64左右達到最大，量值為55個阻力單位。

圖6 不同升力系數下平尾配平引起的阻力增量Fig.6 Delta effects of drag and and corresponding lift values

3.2 CASE 2計算結果及分析

計算狀態如下：

·M＝0.70，0.75，0.80，0.83，0.85，0.86，0.87

·CL＝0.4，0.45，0.5（±0.001）

·中等網格

·Re＝5.0×106（基于MAC＝7.00532m）

基于中等規模的網格，在不同升力系數下，WBH0的阻力特性隨馬赫數變化趨勢如圖7所示?？梢钥闯?，馬赫數小于0.83，阻力系數變化平穩；馬赫數繼續增大，阻力特性則迅速惡化。

圖7 阻力系數隨來流馬赫數變化曲線（ih＝0）Fig.7 Drag rise curves for various Mach numbers（ih＝0）

4 結 論

本文采用TRIP（TRIsonic Platform）軟件和結構網格技術，數值模擬了DPW4（AIAA 4thDrag Prediction Workshop）翼／身／平尾組合體構型，詳細研究了網格密度、馬赫數、平尾迎角對總體氣動特性和壓力分布的影響，通過與相應的試驗結果和計算結果的比較，得到以下一些基本結論：

（1）本文的計算結果與DPW4計算結果的平均值取得了較好的一致，風洞試驗結果與計算結果在迎角、力矩特性方面的差異有待于進一步分析；

（2）網格密度對總阻力系數有一定的影響，相比較而言對壓差阻力的影響大一些，對摩擦阻力系數影響較小；

（3）在合理的網格分布前提下，網格密度對壓力分布的影響不明顯，粗網格、中等網格和密網格的壓力分布幾乎完全重合。

致謝：感謝TRIP軟件的其他開發人員，張玉倫、洪俊武和王光學等同志在多重網格技術和湍流模型方面的研究工作。

［1］TINOCO E N，BOGUE D R.Progress toward CFD for full flight envelope［J］.AeronauticalJournal，2005：451－460.

［2］WILLIAM L O，TIMOTHY G T.Verification and validation in computational fluid dynamics［J］.Progressin AerospaceSciences，2002，38：209－272.

［3］HAASE W，CHAPUT E，LESCHZINER M A.The European high lift project EUROLIFT Ⅱ－objectives，approach，and structure［R］.AIAA 2007－4296.

［4］VASSBERG J C，DEHAAN M A，RIVERS S M，et al.Development of a common research model for applied CFD validation studies［R］.AIAA 2008－6919.

［5］王運濤，王光學，張玉倫.采用TRIP2.0軟件計算DLR＿F6構型的阻力［J］.空氣動力學學報，2009，27（1）：108－113.

［6］ftp：／／cmb24.larc.nasa.gov／outgoing／DPW4／multiblock＿boeing

［7］Van LEER B.Towards the ultimate conservation differencescheme II，monoticity and conservation combined in a second order scheme［J］.J.Comp.Phya.，1974，14：361－370.

［8］MENTER F R.Improved two－equationk－ωturbulence modelsfor aerodynamic－flows［R］.NASA TM－103975，1992.

［9］ELIASSON P，PENG S H.Influence of turbulence modelling and grid resolution in computations of the DPW－4CRM configuration［R］.AIAA 2010－1416.

［10］VASSBERG J C，TINOCO E N，et al.Summary of the fourth AIAA CFD drag prediction workshop［R］.AIAA 2010－4547.