李 偉, 王運濤, 洪俊武,*, 孟德虹, 李 樺

(1. 國防科技大學 空天科學學院, 長沙 410073; 2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 綿陽 621000)

0 引 言

基于Navier-Stokes方程的數值模擬技術已經被廣泛地用于多種飛行器復雜流場模擬，而對數值模擬技術進行驗證確認以明確下一步發(fā)展方向的研究，一直是近年來CFD研究的熱點[1]。國際上已經有很多組織舉辦了CFD驗證確認的研討會，如歐洲計算動力研討會[2](European Computational Aerodynamics Research Project, ECARP)、美國的阻力預測會議[3-7](Drag Prediction Workshop, DPW)、高升力預測會議[8](High Lift Prediction Workshop, HiLiftPW)等。國內CFD工作者也進行了許多關于CFD驗證確認的研究工作[9-14]。

為了進一步促進國內CFD驗證和確認工作的穩(wěn)步開展，評估CFD當前技術狀態(tài)，探索CFD的發(fā)展方向。2018年4月，中國空氣動力研究與發(fā)展中心(CARDC)計算空氣動力研究所聯合國內相關單位組織了第一屆航空CFD可信度研討會(1st Aeronautic CFD Credibility Workshop，AeCW-1)。

大會選擇的運輸機標模構型(CHiNa-Transport, CHN-T1)，是由CARDC自行設計的一種翼身平立尾 (Wing/Body/Tails, WBT)單通道運輸機巡航構型，其外形如圖1所示[15]。

圖1 CHN-T1翼身平立尾構型Fig.1 CHN-T1 wing/body/tails configuration

該構型于2016年9月在CARDC跨聲速風洞(FL-26)完成了風洞試驗[16]，該風洞是一座試驗段尺寸為2.4m×2.4m的引射式、半回流、暫沖式跨聲速增壓風洞。CHN-T1模型的風洞試驗照片如圖2所示[15]。另外，該模型還在歐洲跨聲速風洞(European Transonic Windturnnel, ETW)進行了高雷諾數風洞試驗，Re=15.0×106。

圖2 CHN-T1標模在2.4米風洞中試驗照片Fig.2 Photo of CHN-T1 standard model tested in 2.4 m wind tunnel

大會提供的計算模型與風洞試驗模型尺寸相同：展長為1.548 2 m, 平均氣動弦長為0.193 7 m, 機身長度為1.574 4 m。另外，組委會還提供了帶支撐裝置的計算模型(Wing/Body/Tails/Support, WBTS)以及通過風洞測量機翼的變形數據構造得到的靜氣彈計算模型(Wing/Body/Tails/Support/Deforming, WBTSD)。

另外，會議組委會確定了本次會議的研究內容：針對標模WBT構型進行網格收斂性研究(case1)；針對標模WBTS構型進行有/無支撐影響研究(case2a、case2b)；針對標模WBTSD構型進行靜氣動彈性研究(case2c)；針對標模WBTS構型進行雷諾數影響研究(case3)。

本文基于由CARDC自主研發(fā)的TRIP3.0軟件平臺，針對大會提出的計算狀態(tài)，開展了CHN-T1模型的數值模擬，主要目的是評估支撐裝置、機翼靜氣動彈性變形和雷諾數效應對CHN-T1模型氣動特性數值模擬結果的影響。

1 計算方法與計算網格

1.1 計算方法

本文的研究工作基于CARDC研發(fā)的亞跨超CFD軟件平臺(TRIP 3.0)開展。TRIP軟件經過多年的應用發(fā)展[17-20]，已經成為一個非常成熟的高精度數值模擬軟件平臺。

本文的研究采用有限體積法和結構網格技術求解曲線坐標系下的雷諾平均Navier-Stokes方程(RANS)，控制方程的無粘項采用二階MUSCL[22]差分格式進行離散，并利用Roe格式進行通量分裂；粘性項采用二階中心格式進行離散；離散方程組的求解采用LU-SGS[23]方法。計算中采用SST兩方程湍流模型進行全湍流模擬。另外，采用多重網格技術和大規(guī)模并行計算技術加速收斂。

1.2 計算網格

本文研究中所使用的網格均為采用商業(yè)軟件生成的多塊結構網格，且提交給會議組委會作為大會的結構網格標準網格。進行網格生成時，單個網格塊上三個維度各維度網格點數均為4N+1，以確保能夠進行三重網格計算加速收斂。為確保網格能有效的捕捉流場細節(jié)，表面網格分布按照網格生成規(guī)范[24]進行：翼型前/后緣網格流向尺度約為當地弦長的0.1%，機翼翼根和翼尖約為半展長的0.1%，附面層內網格增長率以及流線的拉伸率應小于1.25，附面層第一層網格高度滿足y+<1。空間網格拓撲采用H型，附面層網格拓撲采用O型，而在翼身結合部和平立尾與機身的結合部位附面層網格采用H型，以防止出現數值模擬流動分離區(qū)域偏小的現象[25]。

針對會議組委會確定的研究內容，分別生成了對應的網格：網格收斂性研究內容(case1)，首先對CHN-T1標模WBT構型生成了一套中等規(guī)模(Medium)多塊對接網格作為基準網格，并以此為基礎進行網格重構，從而分別獲得粗網格(Coarse)和細網格(Fine)。(另外，本文還額外生成了一套百億級的極細(Extremely fine)計算網格進行模擬，來進一步驗證網格收斂性)。其網格拓撲及網格表面分布如圖3所示；有/無支撐影響研究內容(case2a、case2b)，選擇WBTS構型在WBT構型中等網格基礎上，生成了中等規(guī)模的帶支撐網格，其網格拓撲及網格表面分布如圖4所示；靜氣動彈性影響研究內容(case2c)，選擇WBTSD構型在case2b網格基礎上，重新貼體生成了靜氣動彈性變形后的網格case2c，圖5給出了靜氣彈網格與未變形網格在機翼上的差異；雷諾數影響研究內容(case3)，在WBTS構型網格基礎上，更改附面層內第一層距離與網格點數，得到高雷諾數下的網格case3。不同網格的具體參數見表1。

表1 CHN-T1構型不同網格參數Table1 Grid parameters for CHN-T1 configuration

圖3 CHN-T1構型case1網格拓撲及表面網格Fig.3 Surface mesh and topology of CHN-T1 configuration for case1

圖4 CHN-T1構型case2b網格拓撲及表面網格Fig.4 Surface mesh and topology of CHN-T1 configuration for case2b

圖5 CHN-T1構型case2c表面網格Fig.5 Surface mesh of CHN-T1 configuration for case2c

2 CHN-T1標模網格收斂性研究

為了研究定升力系數下網格收斂性，采用標模WBT構型的粗、中、細3套網格以及百億極細網格，為了便于表述分別記為case1-c、case1-m、case1-f，case1-e。四套網格的具體參數上一節(jié)已經給出。來流條件為：Ma=0.78，CL= 0.500(±0.001)，基于平均氣動力弦長(MAC=0.193 7 m)的雷諾數為Re=3.3×106。

圖6給出采用粗、中、細三套網格以及百億極細網格獲得的阻力系數收斂歷程，其橫坐標為迭代步數(Iter)。從圖中可以看出，每套網格均可得到收斂的氣動力結果。其中三套網格在5000步左右阻力系數已經接近收斂，收斂性良好(至10 000步時阻力系數無變化)，且隨著網格的加密，阻力系數單調減??；百億極細網格在12 000步左右阻力系數接近收斂，收斂性良好(至70 000步阻力系數無變化)，且阻力系數收斂結果幾乎與細網格的值重合。可以看出，阻力系數滿足網格收斂性。

圖6 不同網格密度下的阻力系數收斂歷程Fig.6 Convergence history of drag coefficient

表2給出了不同網格密度下得到的迎角α、阻力系數CD、升力系數CL、壓差阻力系數CD,p、摩擦阻力系數CD,f、理想阻力系數CD,idea=CD-CL2/(ARπ)、和俯仰力矩系數Cm。表中同時給出了數值結果的插值解, 其值采用Richardson插值(Richardson Extrapolation, R.E.)方法[23,25]基于中等網格和細網格的結果得到。圖7給出了不同網格密度下獲得的阻力系數、定升迎角、俯仰力矩系數隨網格密度的變化曲線，其中橫坐標為網格單元數的-2/3次冪。

表2 CHN-T1構型case1氣動特性Table 2 Aerodynamic characters of CHN-T1 configuration for case1

從圖表中結果可以看出，對于粗中細三套網格，隨著網格加密，阻力系數及其摩阻分量與壓阻分量、定升迎角均單調變化(定升迎角增大，阻力系數減小，摩阻增大，壓差阻力減小，力矩系數減小)，結果體現出良好的收斂性。因此，后續(xù)的研究內容均采用中等網格規(guī)模進行模擬。

百億極細網格的結果除壓差阻力分量外，其余結果均較插值結果略偏大，不過其阻力系數及定升迎角等結果均滿足上述收斂性，只有力矩系數結果偏離了前面的單調性。由于百億極細網格力矩系數的偏離幅度量值極小(較細網格力矩系數結果偏大0.000 42)，而影響力矩系數變化的因素較多也十分復雜，而且此類超大規(guī)模網格的數值模擬也缺乏足夠的數據對照，因此百億極細網格力矩系數未保持單調變化的原因有待進一步研究。

(a) Grid convergence of CD

(b) Grid convergence of α

(c) Grid convergence of Cm

3 CHN-T1標模支撐裝置影響研究

CHN-T1標模在進行風洞試驗時，需在機身后體安裝支撐裝置以支持測量天平的數據測量。其外形如圖2所示。支撐裝置的存在，會影響平尾附近的流動。直接利用標模WBT構型進行模擬與風洞試驗的結果會存在一定程度的差異。需要進行標模WBTS構型的模擬以減小與風洞試驗中的模型差異。

WBTS構型中的支撐裝置與風洞試驗中的支撐裝置一致。為了研究支撐裝置對模型氣動特性的影響，選擇網格case1-m與case2b進行模擬，來流條件為：Ma=0.78，α=-2°～4.5°，基于平均氣動力弦長(MAC=0.1937 m)的雷諾數為Re=3.3×106。

圖8給出了CHN-T1標模有/無支撐構型氣動特性隨迎角變化曲線。從圖中可以看出，有/無支撐構型的氣動力特性差異不大，同迎角下升力系數有支撐構型較無支撐構型增大約0.015，而無支撐構型的升力與試驗結果更接近；在零度迎角后隨迎角增大，有支撐構型的阻力特性較無支撐構型增大幅度逐漸加大，在當前的計算范圍內，阻力系數的變化約在2%，且有支撐構型的阻力與試驗結果更接近。有支撐構型力矩系數較無支撐構型明顯降低，同迎角下力矩系數大約減小0.05，且有支撐構型力矩系數更接近試驗結果。

由于支撐安裝在后機身平尾下方，直觀感受可以知道支撐的存在對平尾流場具有很大的影響。為了更細致地研究支撐裝置對全機的具體影響部件，本文進行了分部件積分。選擇α=3°的結果進行比較分析。

首先給出了各個部件有/無支撐下的氣動特性結果，以及各個部件的氣動特性變化量在整體氣動特性中的比例，如表3。表中結果可以看出，平尾部件所占的比重最大：阻力增大4.03%，升力增大1.96%，俯仰力矩減小126.8%。即添加支撐以后，平尾的氣動特性受到的影響最大。另外，從表中的數值也可以知道，添加支撐以后總體阻力增大1.91%，升力增大2.23%，俯仰力矩減小131.7%。即添加支撐對力矩特性的影響更大。

另外，圖9給出了有/無支撐條件下，平尾展向截面上/下表面壓力系數分布曲線。從圖中可以看出，添加支撐以后，平尾各截面上的壓力曲線包圍的面積均減小，即平尾上的向下的壓力減小(也可認為是升力增大)。由于平尾距離飛機重心的力臂較大，平尾上的壓力變化引起的力矩變化更明顯，因此造成抬頭力矩減小幅度較大。圖10給出了有/無支撐條件下，平尾上表面壓力系數云圖。從圖中可以看出，添加支撐以后平尾上表面壓力減小，平尾下表面壓力增大，因此平尾上的升力增大。與圖9所描述的狀態(tài)相同，平尾上產生一個較大的低頭力矩，使得飛機的俯仰力矩減小。

(a) Drag coefficient

(b) Lift coefficient

(c) Moment coefficient

圖8 case2a和case2b氣動力系數Fig.8 Aerodynamic characters for case2a and case2b

(a) Cp distribution at section 1(η=0.28)

(b) Cp distribution at section 2(η=0.5)

(c) Cp distribution at section 3(η=0.95)

圖9 case2a和case2b平尾上各截面表面壓力系數分布曲線Cp分布
Fig.9Cpdistribution of the horizonal tail at each section for case2a and case2b

(a) Top view

(b) Bottom view

4 CHN-T1標模靜氣動彈性影響研究

CHN-T1標模在進行風洞試驗時，受氣動壓力的影響，隨著迎角的增加，在機翼上會發(fā)生越來越大的彈性變形。機翼形變以后，其氣動特性也會發(fā)生相應的變化。需要進行標模WBTSD構型的模擬，來研究靜氣動彈性變形的影響。

CHN-T1標模WBTSD構型中，不同迎角下機翼形變不同。為了研究靜氣動彈性變形的影響，選擇網格case2b與case2c進行模擬，來流條件為：Ma=0.78，α=-2°～4.5°，基于平均氣動力弦長(MAC=0.193 7 m)的雷諾數為Re=3.3×106。

圖11給出了CHN-T1標模有/無靜氣彈變形氣動特性隨迎角變化曲線。從圖中可以看出，機翼發(fā)生變形后，升阻力特性與力矩特性的變化幅度均不大。且隨著迎角的增加，阻力系數與升力系數減小的幅度增大，力矩系數逐漸增大。

(a) Drag coefficient of case2b and case2c

(b) Lift coefficient of case2b and case2c

(c) Moment coefficient of case2b and case2c

5 CHN-T1標模雷諾數影響研究

CHN-T1標模在進行風洞試驗時，試驗雷諾數(Re=3.3×106)遠小于真實飛行雷諾數。為了研究雷諾數的影響，組委會還在ETW風洞進行了高雷諾數風洞試驗(Re=15×106)。選擇網格case2b與case3進行模擬，來流條件為：Ma=0.78，CL= 0.500(±0.001)，基于平均氣動力弦長(MAC=0.1937 m)的雷諾數分別為Re=3.3×106、Re=15.0×106。

表4給出了不同雷諾數下，定升力的氣動特性結果與風洞試驗結果，以及數值結果與試驗結果的差量ΔS-E和不同雷諾數的數值結果差量ΔRe。從表中結果可以看出，不同雷諾數下，計算結果與試驗結果均符合良好，且高雷諾數下的計算結果與試驗結果的一致性更好；雷諾數增大以后，阻力系數及其摩阻分量和壓阻分量均減小，定升迎角減小，力矩系數增大，其中摩阻分量相比壓阻分量減少量更大。

圖12給出了不同雷諾數下，CL=0.500(±0.001)，機翼與平尾展向截面上/下表面壓力系數分布曲線。從圖中可以看出，雷諾數增大以后，機翼上表面壓力減小，特別是激波后壓力減小明顯，機翼下表面壓力略有增加，因此機翼上升力增大；對于平尾，上表面壓力增大，下表面壓力減小，因此平尾上產生的向下壓力增大，即升力減小，抬頭力矩增大。

表4 case2b和case3，氣動特性Table 4 Aerodynamic characters for case2b and case3

(a) Cp distribution of the wing at section 1(η=0.28)

(b) Cp distribution of the wing at section 2(η=0.5)

(c) Cp distribution of the wing at section 3(η=0.95)

(d) Cp distribution of the tail at section 1(η=0.28)

(e) Cp distribution of the tail at section 2(η=0.5)

(f) Cp distribution of the tail at section 3(η=0.95)

圖12 Case2b和case3機翼與平尾上各截面表面壓力系數分布曲線Cp分布
Fig.12Cpdistribution of the horizonal tail and wing at each section for case2b and case3

6 結 論

本文基于TRIP3.0軟件平臺，采用多塊結構網格技術，通過求解曲線坐標系下的RANS方程，模擬了第一屆航空CFD可信度研討會提供的運輸機標模構型。通過理論分析和試驗結果對比，基本結論如下：

(1) 采用不同網格密度得到的結果，體現出良好的網格收斂性，百億級的極細網格得到的結果依然符合網格收斂性；

(2) 支撐裝置主要影響平尾上的氣動特性，使得升力/阻力系數增大，俯仰力矩系數減小，其中對俯仰力矩系數影響更大；機翼的靜彈性變形使得升力/阻力系數減小，俯仰力矩系數增大，然而對氣動力的影響量均比較小；雷諾數增大使得阻力系數/俯仰力矩系數減小，主要影響氣動力；

(3) 數值模擬結果與試驗結果對比顯示出良好的一致性，表明對于此類運輸機標模構型，TRIP軟件平臺的模擬結果具有良好的適用性。

(4) 考慮了支撐裝置、機翼靜彈性變形后，數值模擬得到的結果與試驗結果仍有差異，表明在對于運輸機標模的模擬上，仍有值得深入研究的問題。

致謝:本文的相關內容是在TRIP軟件開發(fā)小組的共同努力下完成的，在此對課題組成員張書俊副研究員、楊小川工程師、孫巖助理研究員、岳皓研究實習員等同志表示衷心的感謝！