杜一鳴, 龐 超, 舒博文, 高正紅

(西北工業(yè)大學航空學院 翼型葉柵空氣動力學重點實驗室, 西安 710072)

0 引 言

隨著計算流體力學技術在航空領域的廣泛應用，計算的驗證、確認與可信度問題成為CFD繼續(xù)深入應用時不可回避的問題之一[1-2]。學術界對CFD可信度問題的研究源自上世紀90年代計算機技術的飛速發(fā)展，并逐漸形成了以復現標模氣動力數據或流場特征為主要手段的風洞試驗能力和CFD模擬能力相互確認方法[3]，以此評估當前CFD方法在求解復雜氣動問題時的能力，其核心目標是實現準確的巡航阻力預測，這也是當前工業(yè)飛機設計領域對雷諾平均Navier-Stokes(RANS)方法的主要需求。以客機為例，參考影響力較為廣泛的歷屆AIAA阻力預測會議(Drag Prediction Workshop, DPW)[4-5]可以發(fā)現，隨著RANS方法的日趨成熟，其研究內容已不再局限于氣動力預測，還延伸出了分離流動預測、湍流模型適用性、收斂解的網格生成要求等主題，其中網格的發(fā)展和使用始終都是不可忽視的關鍵研究課題。

從第三屆阻力預測會議開始，嵌套網格(Chimera Grids)就已作為標準網格的一部分對外發(fā)布[6-7]?？梢姲殡S著結構和非結構網格的討論，融合了貼體網格和笛卡爾網格優(yōu)勢的嵌套網格已經成為降低網格生成難度，優(yōu)化網格拓撲結構并平衡網格量，同時獲得精確氣動力預測值的策略之一。本文的工作即在討論結構化嵌套網格結合傳統(tǒng)CFD方法對于跨聲速客機構型的氣動力預測能力，評估嵌套網格相比于傳統(tǒng)對接網格的優(yōu)劣，驗證相關求解器的可信度，為今后相關研究提供參考。

由于嵌套網格并非AeCW-1的標準網格，遵循組委會“網格生成指南”自主研制了規(guī)模以3.3倍增長的稀-中等-密-極密四套結構化嵌套網格系列，并采用近年開源的NASA CFL3D求解器和自研OFS3D求解器開展了網格收斂性和抖振特性研究，并與組委會對接網格結果進行了對比研究。

1 CHN-T1氣動標模

CHN-T1(China-Transport 1)[3]是由中國空氣動力研究與發(fā)展中心(CARDC)研發(fā)的用于風洞試驗和CFD可信度確認的、具有窄體機身-超臨界機翼等幾何特征的單通道客機標模。設計全機外形(BWHVNP)包含機身、機翼、平尾、垂尾、短艙、掛架及起落架整流包等部件。

為了突出流場主要特點，此次AeCW-1采用簡化的計算模型，基礎研究構型為“機身-機翼-平尾-垂尾”構型(BWHV，即Config.1)。該構型自2013年以來經過了全面的風洞試驗，氣動力和流場數據可靠。為模擬風洞干擾，帶支撐構型(Config.2)也是研究構型之一，用于支撐、靜氣彈變形和雷諾數影響的計算。表1列出了風洞試驗模型及本文計算采用的主要參數。

表1 BWHV風洞試驗構型主要參數Table 1 Primary parameters of BWHV configuration

2 計算網格生成

目前，如何進行分區(qū)界面的計算信息傳遞是CFD網格技術的核心，主要的處理方式包括對接(1-1 Blocking)、拼接(Patching)和嵌套等。其中，嵌套網格方法[8]具備靈活的分區(qū)策略，各部分通過彼此的交錯、重疊區(qū)域進行信息傳遞，從根本上降低了網格拓撲的生成難度，不但實現了網格量的合理分配，而且能夠更好地保證物面網格和空間體網格的正交性；此外在處理有活動部件的氣動外形(如直升機旋翼)時，嵌套網格具有天然優(yōu)勢。這些也使得嵌套網格方法成為越來越多的商用和自研求解器包含的基礎功能之一。

參照網格收斂性對網格規(guī)模的要求，首先采用自下而上的網格生成策略對BWHV構型生成中等密度網格，再按相同比例系數在空間三個方向分別進行加密或稀疏得到四套密度不同的網格系列。

2.1 中等網格

采用Pointwise軟件生成了如圖1所示物面網格。全機共分為20個網格塊(18塊貼體空間網格+2塊遠場網格)，圖中表格給出了各塊貼體空間網格的拓撲結構。考慮到流動在翼面兩端都有較強的三維效應，所以將翼根和翼梢處的物面網格均沿展向截去，分別用翼根套環(huán)(Collar)和翼尖帽(Cap)網格填充，避免在可能的復雜流動區(qū)域進行邊界插值操作。此外，與DPW會議所采用的機身網格一體化生成策略(圖2)不同，為了避免機身網格“浸入”翼面內部，降低挖洞難度，本文對機身物面網格采用分塊生成策略，這樣可以自主控制翼身連接處各塊網格的邊界，便于進行網格邊界插值。

圖1 CHN-T1 BWHV構型表面嵌套網格分布Fig.1 Overset surface grid layout for BWHV configurations of CHN-T1

圖2 DPW機身一體化嵌套網格生成策略Fig.2 Compact strategy for wing-body grid generation adopted by DPW series

對于機身以及各個翼面的物面網格，首先采用超限插值(TFI)方法生成代數網格，再利用Thomas-Middlecoff方法對其進行橢圓型方程迭代，直至最大殘值降至5×10-6量級；而機頭、機尾、翼尖和翼身連接處物面網格通過求解雙曲型方程推進得到；貼體空間網格則通過沿物面法向推進求解雙曲型方程得到。

以上方法可保證物面附近網格較好的正交性和較高的整體網格質量。遠場網格分為內層和外層兩塊，均為直角笛卡爾網格。為了更好的進行網格裝配，內層網格在機體附近進行了適當加密，加密標準是網格間距略小于貼體空間網格的最外層。圖3展示了幾處典型的局部網格形態(tài)。

2.2 網格序列

為了進行網格收斂性研究，以中等密度網格為基礎，在各個方向進行加密或加粗進而得到密網格、特密網格和稀網格，加密或加粗系數分別為1.5和0.75。圖4給出了機頭、翼身套環(huán)和主翼單塊物面網格疏密情況對比，其中i、j、k分別對應三個維度的網格點數，而“×M”(Medium)和“×F”(Fine)分別代表相應維度點數及該塊網格量與中等網格或密網格對應值的倍數關系。值得注意的是，各套網格的翼身套環(huán)(collar)存在一定差別：該塊機身網格由雙曲型方程推進求解得到，無法精確保證推出距離；而機翼網格截取自wing網格塊，由于不同網格的網格間距差別較大，為了滿足網格點倍數的要求，導致網格在機翼和機身處的推出距離不同。

表2給出了上述自研嵌套網格系列與DPW-IV會議標模CRM翼身-平尾構型嵌套網格[7]的參數對比，可見自研嵌套網格參數與DPW及AeCW-1組委會中等結構化對接網格相近。

(a) Farfield Cartesian blocks (b) Streamwise section of the (c) Grid on the symmetry plane and

overset grid around the wing O-grid of the wing and tails

圖3 典型區(qū)域網格形態(tài)Fig.3 Grid configuration at typical regions

圖4 網格序列局部對比Fig.4 Comparison of local surface grid

2.3 網格裝配

嵌套網格生成之后，要對網格進行挖洞和插值，合稱網格裝配(Grid Assembly)。挖洞使得“浸入”物面內部的網格單元從計算域中移除，常見方法有洞映射(Hole Map)方法[9]和射線(X-ray)求交法[10]等，插值使得未賦予特殊邊界條件的網格塊邊界單元獲得流場信息，常見方法有模板游走法(Stencil Walking)[11-12]、反變換法(Inverse Map)[13]和自適應叉樹法(Alternating Digital Tree，ADT)[14]等。本文利用射線求交法及反變換方法進行網格裝配。在進行網格挖洞時，由于采用了機身網格分塊生成策略，因此只需要對遠場笛卡爾網格進行挖洞。為了保證洞邊界網格單元大小匹配，使用洞面優(yōu)化技術，將機體表面一定距離之內的遠場網格點均作為洞點挖去。挖洞結果如圖5所示。

2.4 帶支撐構型網格

由于支撐與后體機身夾角過小，很難利用雙曲型方程生成空間貼體網格，故采用截取對接網格的方法生成。支撐網格及局部細節(jié)如圖6所示。首先移植組委會中等結構化對接網格對應位置處的部分貼體網格至縫隙處(細節(jié)圖中藍色網格)，再生成一塊笛卡爾網格填補剩余空間(細節(jié)圖中紅色網格)。其余位置拓撲及網格均與基礎構型保持一致，該網格簡記為M2網格。

考慮靜氣彈變形的Case2c在帶支撐構型的基礎上增加了機翼變形。但本文暫未實現自動化的網格變形技術，故未針對Case2c生成嵌套網格。下文僅采用CFL3D在組委會提供的各迎角對接網格上進行了數值評估，以討論靜氣彈變形的影響。

圖5 挖洞后網格形態(tài)Fig.5 Grid after hole cutting

圖6 帶支撐構型支撐處網格Fig.6 Grid details around the support of Config.2

3 求解器與計算平臺

3.1 求解器

CFL3D結構化網格求解器[15]由NASA Langley研究中心于20世紀80年代初開始開發(fā)，并于2017年在GitHub上開源。CFL3D在有限體積框架下求解三維可壓縮非定常RANS方程，能夠處理結構化對接、拼接和嵌套網格，有多種湍流模型可供選擇。

OFS3D結構化嵌套網格求解器由課題組自主研發(fā)，具有二維/三維，定常/非定常，薄層近似/完全N-S方程求解能力，目前已植入了多種空間離散格式及時間推進格式，同樣擁有多種常用湍流模型可以調用。

為保證計算結果的可比性，兩種求解器均調用SA和SSTk-ω湍流模型進行數值模擬。RANS方程對流項和壓強項采用二階迎風Roe格式求解(不加熵修正)，并統(tǒng)一采用完全N-S形式和二階中心差分格式處理粘性項，時間推進格式選用隱式近似因式分解方法。多重網格技術用于加速收斂，CFL數因求解器和網格不同，采用2.0～5.0之間的值，事實證明穩(wěn)定性和收斂性良好。為論述簡潔，兩求解器分別簡記為“C”和“O”。

3.2 計算平臺

兩種求解器均采用Openmpi-v1.8進行并行編譯及計算，并在課題組自建高性能集群上運行。集群搭載主頻3.3 GHz 的Intel志強系列CPU，運行環(huán)境為64位Linux，計算能力約為14 Tflops。

計算效率方面，由于CFL3D不具備網格自動再分區(qū)策略，難以實現負載平衡，僅能調用20個線程，計算效率較低。而OFS3D可根據計算總線程數和網格進行自動網格再分，從而實現負載平衡及高加速比。表3給出了中等網格下調用相同計算資源(20線程)且達到相近收斂階時兩程序計算效率的對比(單位：s)，可見若增加計算資源，OFS3D計算效率還將提高，優(yōu)勢更加明顯。

表3 BWHV風洞試驗構型主要參數Table 3 Primary parameters of BWHV configuration

4 結果與討論

由于計算資源和篇幅的限制，主要對基礎構型的網格收斂性和抖振特性進行重點研究，同時對模型支撐干擾和靜氣彈影響進行初步計算分析，具體計算狀態(tài)見表4。CFL3D由于計算效率較低，暫未完成SA模型特密網格的計算，對于SST模型也僅完成了中等網格的計算。

表4 本文計算狀態(tài)Table 4 Simulation conditionsof present cases

4.1 網格收斂性

4.1.1 氣動力系數

采用Richardson外插方法[6]進行氣動力系數收斂性研究。該方法要求在網格系列參數等效的情況下進行定升力計算，且各套網格的殘差收斂到相同的量級。但由于計算資源有限，不同密度網格很難收斂到相同量級，故本文通過力系數波動大小來判斷是否達到可接受的收斂狀態(tài)，具體的即升力系數≤±10-3，阻力系數≤±10-4，俯仰力矩系數≤±10-3。圖7給出了不同求解器進行網格收斂性計算時密度殘差平均值的收斂曲線。如果求解器具備二階精度，當以網格單元數的-2/3次冪(-1/3次冪代表網格尺度，平方對應二階精度)為橫坐標，預測值為縱坐標時，若得到一條直線，則可定性說明求解器和求解方法整體滿足精度且具備較好的收斂性，且可通過外插得到網格收斂時的結果。參考AIAA阻力預測會議的處理方式[6-7]，采用最密的兩套網格計算結果進行外插。

圖7 兩種求解器密度殘差收斂歷史Fig.7 Convergence history of density residual by two different solvers

圖8分別給出了迎角、總阻力、阻力分量系數和俯仰力矩系數的收斂性和插值結果?？梢?，在升力系數滿足收斂要求的情況下：(1) OFS3D求解器對各個氣動系數的預測均表現出了較好的線性特征，除了迎角和俯仰力矩系數在較密網格上的微小拐折，且采用兩種湍流模型時的總阻力插值點基本重合；(2) CFL3D采用SA模型預測的迎角和俯仰力矩系數收斂性較好，但阻力線性度較差，其稀網格的明顯誤差主要來源于壓差阻力，而采用SST模型時，總阻力和阻力分量與OFS3D-SST結果基本一致；(3) 對于OFS3D迎角和俯仰力矩系數的轉折，這可能與遠離力矩中心的微小分離流動對密網格較敏感有關；(4) 基于本文的嵌套網格序列，SST湍流模型的收斂性整體上不如SA湍流模型。

(a) Angle of attack and lift

(b) Drag

(c) Pressure and skin-friction drag

(d) Pitching moment

表5給出了各個氣動系數的Richardson外插結果。為與對接網格對比，采用CFL3D求解器進行了組委會結構化對接網格的計算及外插(表中“C(1-1)”列)。采用CARDC TRIP 3.0求解器104億結構對接網格(N-2/3×105值已達0.02)70 000步計算結果[16]進行對比可以發(fā)現：(1) 整體上看，OFS3D求解器兩種湍流模型的外插結果更接近百億規(guī)模網格的計算結果，表現出了較好的一致性和收斂性；(2) 基于嵌套網格的CFL3D-SA外插迎角、阻力和俯仰力矩系數與百億網格解差別稍大，阻力系數誤差約3%，但更接近試驗阻力系數0.031 97；(3) 基于對接網格的CFL3D-SST也獲得了較好的收斂狀態(tài)，各氣動系數與百億網格解較接近；(4) 除OFS3D-SA，各結果俯仰力矩系數較百億網格解普遍偏小。整體上看，本文生成的網格序列對于收斂性研究是合適的，采用的求解器也具有較高的可靠性。

由于本文網格序列完全由手動方式分別生成，考慮到局部網格質量和網格裝配精度等因素，各規(guī)模網格的參數并不完全等效，加之收斂量級存在差別，以上因素均會對網格收斂性曲線的單調性產生一定影響。

表5 CHN-T1 BWHV構型Richardson外插結果Table 5 Richardson extrapolation of BWHV configuration

4.1.2 翼面壓力分布

分別截取機翼和平尾各兩個站位(圖9)的壓力分布，其隨網格密度增加的變化趨勢如圖10所示。對于同一求解器、同一湍流模型來說，不同規(guī)模網格的差別僅當翼面存在較強激波時(機翼50%)表現得較為明顯，隨著網格密度的增大，激波厚度逐漸減小。這一趨勢在CFL3D求解器上表現得較為明顯，而OFS3D變化較小，且密網格和極密網格的計算結果幾乎沒有變化。

圖11給出了中等嵌套和對接網格上，不同求解器、不同湍流模型得到的截面壓力分布對比。仍以TRIP 3.0求解器104億結構網格計算結果為參照?？梢钥闯?，同一求解器在相同網格類型下兩種湍流模型得到的壓力分布相近(SST相較于SA預測的激波位置稍靠前、前緣吸力峰稍高)，也即主要差別來源于求解器和網格類型。對于嵌套網格，CFL3D相較于OFS3D預測的激波位置更靠前、厚度更大，同時前緣吸力峰也更高，OFS3D預測結果與百億網格解吻合較好。而對接網格的激波分辨率要明顯低于嵌套網格，85%機翼截面處的弱激波抹平嚴重，這也正是嵌套網格相較于對接網格的優(yōu)勢之一。如圖12所示，由于嵌套網格不存在近場密網格向遠場輻射的問題，在網格量相當時，嵌套網格可將更多的單元分布于機頭、機翼等關鍵位置和流動復雜處，大大提高了局部流場的分辨率。由壓力分布對比可見，中等規(guī)模嵌套網格能夠滿足求解器后續(xù)研究需求，但中等規(guī)模對接網格在激波處的網格密度稍顯不足。

4.1.3 流動分離

跨聲速客機構型的一個典型流動特征就是容易在機翼翼根后緣和機身之間出現分離泡，準確再現這個區(qū)域的流動特征也成為了歷屆DPW會議的主題之一[6-7]。CHN-T1標模在設計時考慮了翼根整流[3]，因此在風洞試驗中未見明顯的分離流動，分離區(qū)明顯小于DLR-F6翼身組合體構型[6]。圖13給出了該區(qū)域的計算表面流線圖，從左至右分別為嵌套網格的C/M/F/E-F網格計算結果。可見相對于模型整體尺寸，分離泡尺度很小(約2～9倍后緣厚度)。從稀網格到中等網格，分離在展向和流向有明顯擴大，隨后在更密的網格上基本趨于穩(wěn)定。特別地，OFS3D-SST的分離區(qū)在密網格上仍有明顯擴展，并在特密網格上有所收縮，這也導致了其俯仰力矩在密網格上的拐折(圖8(d))。對于同一規(guī)模的網格，CFL3D-SA的分離區(qū)在流向大于OFS3D-SA；OFS3D-SST的分離區(qū)在流向和展向均明顯大于OFS3D-SA的計算結果。

圖9 CHN-T1 BWHV構型機翼尾翼站位示意Fig.9 Section stations of CHN-T1 BWHV wing and tail

圖10 CHN-T1 BWHV構型網格收斂性研究：截面壓力分布Fig.10 CHN-T1 BWHV grid-convergence study: sectional pressure distribution

圖11 中等網格CHN-T1 BWHV構型截面壓力分布對比Fig.11 Comparison of CHN-T1 BWHV sectional pressure distribution in medium grid

圖12 中等嵌套和對接網格局部分布對比Fig.12 Comparison of local grid distribution in medium overset and 1-1 blocking grid

4.2 抖振特性

開展各構型氣動特性隨迎角變化的計算，不但易與風洞試驗結果直接對比，掌握支撐、氣彈變形等對氣動特性的影響程度，進一步評估求解器的可信度，還可以檢驗設計抖振邊界，反向校核風洞試驗的準確性。經過網格收斂性研究，本節(jié)計算結果均采用中等規(guī)模網格獲得。試驗數據來自文獻[17]。

4.2.1 基礎構型氣動特性

圖14給出了升阻力系數和理想阻力系數與試驗數據的對比?？梢园l(fā)現，各求解器得到的升力線斜率與試驗值吻合較好，并且較理想地復現了非線性段起始區(qū)域；OFS3D雖然在線性段與試驗值最為吻合，但最大升力系數明顯偏低。從試驗數據和數值模擬結果可以判斷，CHN-T1標模的抖振邊界大于1.5倍巡航升力系數(即0.65)，達到設計狀態(tài)[3]。

除嵌套網格CFL3D-SA外，其余求解器和湍流模型得到的阻力系數基本全部落在了試驗誤差帶內，

(a) SST model of OFS3D solver

(b) SA model of OFS3D solver

(c) SA model of CFL3D solver

圖14 CHN-T1 BWHV構型氣動特性：升阻力Fig.14 Aerodynamic characteristics of CHN-T1 BWHV: lift and drag

而各預測值與理想阻力系數試驗值均吻合較好。整體來看，SA模型得到的最大升力系數和阻力系數均較SST模型大。

圖15給出了升阻極曲線與試驗數據的對比。各預測值基本聚集在試驗結果附近，僅在迎角較大時差別稍大。從設計點附近的結果可以看出，各預測值均與試驗值吻合較好，相互之間的最大差別在8counts以內，其中OFS3D和對接網格CFL3D的結果與試驗值最為接近。

4.2.2 風洞模型支撐與靜氣彈變形的影響

雖然基礎構型升阻特性預測結果較為理想，但由于沒有考慮模型支撐及靜氣彈變形的影響，俯仰力矩特性與試驗數據相差較大。為此，分別采用OFS3D在嵌套網格上對帶支撐構型，以及CFL3D在對接網格上對考慮靜氣彈變形的帶支撐構型進行評估。圖16給出了三種構型的各組預測結果?？梢?，對于BWHV構型，各組結果重合度較好，但相比試驗值明顯偏大；采用考慮靜氣彈變形的帶支撐構型在小迎角范圍內能夠得到更準確的俯仰力矩特性，但迎角大于3°后與試驗值仍存在較大差距，而僅考慮支撐干擾的構型在此迎角范圍的預測值與試驗值則較為吻合。有關模型支撐和靜氣彈變形的影響仍需開展進一步的計算研究。

圖15 CHN-T1 BWHV構型氣動特性：升阻極曲線Fig.15 Aerodynamic characteristics of CHN-T1 BWHV: drag polar

圖16 CHN-T1氣動特性：俯仰力矩Fig.16 Aerodynamic characteristics of CHN-T1: pitching moment

5 結 論

嵌套網格相比于對接網格簡化了拓撲劃分難度，不存在近場密網格的遠場輻射問題。對于跨聲速運輸機構型，能夠將更多的網格布置于機頭、機翼等關鍵部位，提高激波等復雜流場結構的分辨率，同時保持計算可靠性。本文自研的嵌套網格系列能夠滿足網格收斂性研究的需求，但網格系列的構造方式仍有待改進，如何自動化地實現網格布點及密度調整是今后研究的重點，這樣不但能夠確保網格參數等效，還可以提高網格生成質量和效率。經過計算分析，對于超臨界運輸機構型的數值模擬，初步得到以下結論：

(1) 氣動力和力矩系數在網格系列上存在斜率明顯的收斂曲線，較大規(guī)模的網格不但能夠得到更接近收斂的氣動特性預測值，翼面等關鍵部位較小的網格尺度還能有效提高前緣吸力峰、激波和分離區(qū)的分辨率。

(2) 自研OFS3D求解器前后處理和計算效率較高，具備較好的嵌套網格收斂性和可靠性。而CFL3D嵌套網格功能和工具相對缺乏，計算效率較低。同時，嵌套網格CFL3D-SA計算結果并不理想，其嵌套網格計算能力有待進一步驗證。

(3) 初步分析表明，湍流模型對壓力分布、分離區(qū)大小和最大升力系數等的預測有一定影響，對分離區(qū)大小影響較大。更深入的湍流模型精度和可信度比較還需要進一步的計算分析。

(4) 對于采用常規(guī)測力方法的風洞試驗，數值模擬時有必要考慮模型支撐和機翼靜氣彈變形的影響，這對于獲得準確的力矩特性十分必要。

致謝:感謝清華大學航空航天學院李浩然同學提供的組委會結構密網格CFL3D求解器計算結果。西安理工大學土木建筑工程學院聶勝陽老師在標模計算過程中與作者進行了積極討論，在此一并致謝！