陳江濤，趙嬌，章超，劉深深，張耀冰，吳曉軍

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

氣動力和力矩的精細預測是商用運輸機氣動性能評估和外形優(yōu)化設計的關(guān)鍵因素之一。近些年來隨著網(wǎng)格生成技術(shù)、流場求解和高性能計算機等因素的持續(xù)發(fā)展，CFD在預測復雜外形氣動性能方面取得了長足的進步，評估三維真實運輸機外形氣動性能成為可能。為了評估現(xiàn)階段CFD預測商用運輸機氣動性能的能力，AIAA應用空氣動力學委員會在2001年發(fā)起了AIAA阻力預測會議，來自世界各地的研究機構(gòu)和公司使用不同的程序、網(wǎng)格、湍流模型和數(shù)值格式對翼身組合體外形進行了模擬研究[1-6]。

在復雜工程外形的數(shù)值模擬中，網(wǎng)格類型、規(guī)模和分布、湍流模型、數(shù)值格式等，都會不同程度影響模擬結(jié)果。如何評估這些模擬方法對模擬結(jié)果的影響，并識別對模擬結(jié)果影響較顯著的因素，對于后續(xù)CFD發(fā)展努力的方向有積極的借鑒指導意義。Rumsey和Allison[7]在抖振邊界條件下對某民用運輸機進行了數(shù)值模擬研究，考慮了不同軟件、湍流模型和機翼后緣處的網(wǎng)格處理等影響因素，發(fā)現(xiàn)湍流模型對模擬結(jié)果的影響最大。國內(nèi)也有很多學者從不同的角度，比如離散精度[8-9]、網(wǎng)格拓撲[10]、數(shù)值方法[11-12]等，研究了不同模擬方法對阻力會議外形模擬的影響。但是這些研究大多沒有考慮多種因素同時作用下對預測的影響，比如Rumsey和Allison[7]在研究湍流模型對模擬影響的時候，保持使用的網(wǎng)格和程序不變，這樣得到的結(jié)論有一定的局限性。因此如何在多種影響因素共同作用下評估模擬結(jié)果的變化規(guī)律，識別出對模擬影響較大的因素是本文試圖解決的問題。

為了綜合研究不同因素對商用運輸機外形阻力預測的影響，本文以AIAA第六屆阻力會議外形NASA CRM[13]為研究對象，同時考慮了網(wǎng)格、湍流模型、無黏通量格式和體心梯度求解方法等因素，使用枚舉法和正交試驗設計兩種方法研究了不同因素對阻力預測的影響大小，識別出了關(guān)鍵因素，為下一步方法的發(fā)展指出了方向。

1 計算外形、網(wǎng)格和方法

本文考慮的外形是AIAA第六屆阻力會議的NASA CRM模型，如圖1所示。該模型代表了典型的現(xiàn)代跨聲速運輸機，設計巡航狀態(tài)是馬赫數(shù)Ma=0.85，升力系數(shù)CL=0.5。該外形在NASA NTF風洞、NASA Ames 11 ft風洞[14-15]和歐洲ETW風洞[16]都進行了測試，有大量數(shù)據(jù)可供CFD確認用。但是由于不同風洞運行環(huán)境的差異等因素，試驗數(shù)據(jù)也存在著明顯的不確定性。圖2給出了NTF風洞和Ames風洞多個車次的試驗結(jié)果對比，通過插值可以得到該模型在定CL=0.5條件下在NTF風洞測得的阻力系數(shù)CD為0.024 90～0.024 97，在Ames風洞測得的阻力系數(shù)CD為0.024 14～0.024 34。本文更多關(guān)注的是模擬方法對阻力預測的影響，沒有通過試驗結(jié)果評判模擬方法的優(yōu)劣。

圖1 NASA CRM外形

圖2 CRM不同風洞試驗數(shù)據(jù)對比

這里關(guān)注的是定升力系數(shù)CL=0.5條件下風洞試驗狀態(tài)的阻力預測，來流雷諾數(shù)為5×106。研究中綜合考慮了網(wǎng)格規(guī)模、湍流模型、無黏通量格式和體心梯度求解方法對阻力預測的影響。下面分別介紹不同因素的情況。

首先使用了依次加密的3套非結(jié)構(gòu)網(wǎng)格來進行模擬，網(wǎng)格嚴格按照組委會提供的規(guī)范來生成，加密過程中保持網(wǎng)格整體加密，最粗的一套網(wǎng)格表面如圖3所示，網(wǎng)格信息列于表1。

計算使用課題組發(fā)展的MFlow程序[17], 這是基于二階有限體積法的非結(jié)構(gòu)網(wǎng)格求解器，能夠處理多種網(wǎng)格單元類型。單元內(nèi)使用線性重構(gòu)使得空間精度達到二階。研究中使用了兩種梯度求解方法來計算體心的梯度，分別為基于節(jié)點的格林高斯方法[18](GN)和基于節(jié)點的最小二乘方法[19](LS)?？紤]了3種不同的無黏通量格式，分別為Roe[20]、AUSMPW+[21](圖片標識中簡寫為AUS)、HLLC+[22](圖片標識中簡寫為HLL)。計算中假定流動為全湍流，使用SA一方程模型[23]和剪切應力輸運(SST)k-ω兩方程模型[24]。兩種模型都是使用了QCR修正[25-26]之后的版本。研究發(fā)現(xiàn)[8,27-28]，對于該外形，QCR修正能夠提高翼身結(jié)合處流動的模擬精度，計算結(jié)果與試驗結(jié)果更加吻合。

圖3 粗網(wǎng)格表面網(wǎng)格

表1 計算使用的網(wǎng)格信息

本文的計算從自由來流開始，所有的計算密度殘差都下降了兩個量級以上，氣動力的振蕩在0.2% 以內(nèi)，如圖4所示，因此計算收斂性引入的影響可以忽略不計。

圖4 密度殘差和氣動力系數(shù)收斂曲線

2 阻力預測影響分析

2.1 枚舉法結(jié)果

本文計算考慮了3種網(wǎng)格、3種無黏通量格式、2種湍流模型、2種體心梯度求解方法對阻力預測的影響。本節(jié)將所有因素可能的搭配都進行了模擬，計算安排如表2所示, 其中網(wǎng)格C代表粗網(wǎng)格，M代表中等網(wǎng)格、F代表密網(wǎng)格。

對于總阻力，使用粗網(wǎng)格、SA模型和GN方法的預測值高于上限值。對于壓差阻力，使用粗網(wǎng)格、SSTk-ω模型、AUSMPW+和GN方法的預測值高于上限值，使用密網(wǎng)格、SA模型和LS方法的預測值低于下限值。對于摩擦阻力，均在上下限之間。需要指出的是，計算結(jié)果在統(tǒng)計的上下限之外不代表該結(jié)果是有問題的。

表2 計算安排設置

圖5 阻力系數(shù)預測值的統(tǒng)計結(jié)果

圖6給出了不同模擬方法搭配下的網(wǎng)格收斂性曲線。從圖中可以看出隨著網(wǎng)格加密，總阻力和壓差阻力減小，摩擦阻力增大。摩擦阻力收斂圖上SA和SSTk-ω的結(jié)果區(qū)分的很明顯。壓差阻力圖上湍流模型和梯度求解方法影響都很大，特別是湍流模型?？傋枇D上梯度求解方法影響很大，隨著網(wǎng)格加密不同湍流模型的解之間的差別變小。

圖6 阻力預測的網(wǎng)格收斂特性曲線

圖7和圖8分別給出了定升力狀態(tài)下的迎角α和俯仰力矩系數(shù)Cm隨網(wǎng)格加密的變化趨勢。除了使用SA模型得到的迎角會有非單調(diào)變化,大部分結(jié)果都呈現(xiàn)單調(diào)的變化。總的來說，迎角和俯仰力矩隨著網(wǎng)格加密變化都不大，SA模型計算得到的迎角因網(wǎng)格因素導致的變化不超過0.004°，SST模型為0.02°，SA模型得到的俯仰力矩系數(shù)變化不超過0.000 3，SST模型為0.002。網(wǎng)格收斂性研究基于Taylor展開，需要指定網(wǎng)格尺度h。在三維復雜外形計算中，經(jīng)常使用N-1/3來近似h。假定粗、中、細3套網(wǎng)格計算得到的目標變量分別為f3、f2、f1, 3套網(wǎng)格的網(wǎng)格尺度分別為h3、h2、h1。通過簡單的證明可以知道：只有當0<(f1-f2)/(f2-f3)ln(h2/h1)/ln(h3/h2), 這時計算結(jié)果隨著網(wǎng)格加密呈現(xiàn)單調(diào)發(fā)散的趨勢，基于Taylor展開的網(wǎng)格收斂性分析不再適用。這里需要指出的是，用N-1/3來近似h可能會給網(wǎng)格收斂性分析帶來較大的誤差，特別是在非結(jié)構(gòu)混合網(wǎng)格的框架下，如何構(gòu)造相容的網(wǎng)格序列開展網(wǎng)格收斂性研究也是驗證與確認領(lǐng)域研究的難點。

圖8 俯仰力矩系數(shù)預測的網(wǎng)格收斂特性曲線

對于不同模擬方法對阻力預測的影響，可以定性地做出判斷。首先考慮湍流模型的影響，在固定網(wǎng)格、無黏通量格式和體心梯度求解方法的前提下，對于壓差阻力，SSTk-ω模型比SA大，隨著網(wǎng)格加密差量增大；對于摩擦阻力，SSTk-ω模型比SA小，但隨著網(wǎng)格加密差距減小，總的量級大于壓差阻力，因此總阻力SSTk-ω模型比SA小，這點在粗網(wǎng)格上更加明顯，隨著網(wǎng)格加密差距減小。

考慮梯度求解方法的影響：GN比LS壓差阻力大，摩擦阻力稍大，總阻力大，這點在粗網(wǎng)格上和使用SSTk-ω模型時更加明顯。

而無黏通量格式對阻力預測的影響幾乎可以忽略不計。需要指出的是這一結(jié)論是基于使用的這3種格式而言的，如果使用其他與這幾種格式性能有明顯差別的格式，比如Steger-Warming、Van-Leer等，很有可能得到不一致的結(jié)論。

2.2 基于聚類算法的敏感性分析

4種模擬因素同時作用下，需要對每種因素的影響大小進行量化和判斷，識別出對模擬結(jié)果影響較大的因素，為深入研究提供努力的方向。

本節(jié)使用了Kmeans聚類算法[29]來分析目標變量和輸入因素之間的潛在聯(lián)系。Kmeans聚類算法通過衡量不同點間的歐式距離來表征樣本點之間的相似度，被聚為一類的樣本點具有較高的相似度，通過分析同類樣本中輸入因素的特征，來辨別不同因素的影響大小。

分別針對總阻力、壓差阻力和摩擦阻力，將所有樣本點的目標變量分為兩類，如圖9～圖11所示。黑色實線為大體的兩類之間的分割線，黑色虛線分別為兩類樣本目標變量的中心值，不同顏色的標志代表各因素的不同方法。

以總阻力為例，圖9給出了所有樣本通過4個因素分類的情況，并用黑色實線將各個因素的結(jié)果分為兩類。在圖9(a)中第2類包括了全部細網(wǎng)格、大量中網(wǎng)格和少量粗網(wǎng)格樣本，不能有效區(qū)分網(wǎng)格因素；圖9(b)中SA樣本幾乎均等的分布在第1類和第2類，不能有效區(qū)分湍流模型的影響；圖9(c)中3種無黏通量格式樣本同樣幾乎均等分布在第1類和第2類，不能有效區(qū)分無黏通量格式的影響，在圖9(d)中，第1類大部分為GN結(jié)果，第2類大部分為LS結(jié)果，起到了有效區(qū)分，所以梯度求解方法對于總阻力的影響最大。類似地，對于壓差阻力，湍流模型發(fā)揮了最大的影響，對于摩擦阻力，湍流模型發(fā)揮了最大的影響。

圖9 總阻力系數(shù)的聚類分析結(jié)果(分2類)

本文計算中，有的模擬因素比如網(wǎng)格和無黏通量格式，都存在3種選項，是否應該將所有樣本歸為3類，結(jié)果如圖12～圖14所示。對于總阻力，梯度求解方法發(fā)揮了最大的影響，第1類全部為GN的結(jié)果，第3類全部為LS的結(jié)果，而其他因素在3個 類里都存在不可忽視的數(shù)量標識。對于壓差阻力，湍流模型發(fā)揮了最大的影響，第1類全部為SSTk-ω結(jié)果，第3類全部為SA結(jié)果，而其他因素在3個類里都存在不可忽視的數(shù)量標識。對于摩擦阻力，同樣是湍流模型發(fā)揮了最大的影響，第1類 全部為SA結(jié)果，第2類和第3類全部為SSTk-ω結(jié)果，而其他因素在3個類里都存在不可忽視的數(shù)量標識。

在本算例中，將所有樣本點不管是分為2類還是3類，都可以得到一致的敏感性結(jié)論。即對于總阻力，梯度求解方法發(fā)揮了最大的影響。對于壓差阻力和摩擦阻力，湍流模型發(fā)揮了最大的影響。

2.3 基于McKay的敏感性分析

聚類算法只是定性地給出了因素的敏感性分析，下面使用McKay主影響分析來量化4種因素各自的敏感性，定義為

圖10 壓差阻力系數(shù)的聚類分析結(jié)果(分2類)

(1)

式中：V(Y)表示目標變量的方差估計；V[E(Y|Xi)]表示目標變量在因素Xi作用下條件期望的方差估計。該指標量化了因素Xi導致的目標方差的變化。表3分別給出了壓差阻力、摩擦阻力和總阻力預測中各因素的敏感性指標。

從表3可以分析出，對于壓差阻力，湍流模型對預測的影響最大，其次是網(wǎng)格和梯度求解方法，無黏通量格式的影響很??；對于摩擦阻力，湍流模型對預測的影響最大，網(wǎng)格、無黏通量格式和梯度求解方法對結(jié)果的影響很小，均可以忽略不計；對于總阻力，梯度求解方法的影響最大，其次是網(wǎng)格和湍流模型，無黏通量格式的影響可以忽略不計。

圖12 總阻力系數(shù)的聚類分析結(jié)果(分3類)

圖13 壓差阻力系數(shù)的聚類分析結(jié)果(分3類)

雖然在壓差阻力和摩擦阻力預測中，湍流模型都是影響最大的因素，但是在總阻力預測中，湍流模型并不是影響最大的因素。這里考慮是因為湍流模型對壓差阻力和摩擦阻力的影響趨勢相反。具體來說，SA模型預測的壓差阻力較SSTk-ω模型預測偏小，但摩擦阻力預測的偏大，這樣導致湍流模型導致的總阻力變化較小。而梯度求解方法對于壓差阻力和摩擦阻力的影響是一致的，GN預測的壓阻和摩阻都比LS偏大，所以總阻力預測中，這一趨勢被放大。

圖14 摩擦阻力系數(shù)的聚類分析結(jié)果(分3類)

表3 McKay主影響指標的參數(shù)敏感性

需要指出的是對模擬因素敏感性的分析都是基于使用的這幾種網(wǎng)格、模型和格式，對這一特定算例分析得到的。如果使用不同的模擬方法，比如其他的網(wǎng)格、模型和格式，或者針對新的算例，可能會得到不一致的結(jié)論。但本文給出了一般的因素敏感性分析方法，針對具體問題，可以仿照這一流程進行分析。

2.4 正交試驗設計分析

在上文的分析中，考慮了4種影響因素對阻力預測的影響，每種因素考慮了2～3種不同的選項，所有可能選項遍歷下來計算量不小，因此在本節(jié)中嘗試使用更有效率的正交試驗設計[30]來完成上述分析。在正交試驗設計中，試驗因素為網(wǎng)格、湍流模型、無黏通量格式和梯度求解方法，各因素劃分為表4所示的不同水平因子。

首先不考慮各個因素之間的交互作用，使用極差來辨別因素的主次，選擇正交表L9(34)，采用擬水平法進行正交試驗設計，試驗設計方案如表5所示。

(2)

式中：pj為第j個因素的水平數(shù)目。極差越大，則說明該因素對結(jié)果的影響越大，即該因素越重要。

表4 正交試驗設計的因素與水平因子

表5 擬水平法下的L9(34)正交表

通過極差分析可得不同因素對結(jié)果影響的主次順序：① 對于壓差阻力，湍流模型的影響>網(wǎng)格的影響>梯度求解方法的影響>無黏通量格式的影響；② 對摩擦阻力，湍流模型的影響>網(wǎng)格的影響>梯度求解方法的影響>無黏通量格式的影響；③ 對總阻力，梯度求解方法的影響>網(wǎng)格的影響>湍流模型的影響>無黏通量格式的影響。這里得到的結(jié)論與枚舉法的敏感性分析結(jié)果一致，但是計算量只是枚舉法的1/4。

除此之外，還可以看出在不同水平下網(wǎng)格、湍流模型對壓差阻力和摩擦阻力的影響趨勢是相反的，而梯度求解方法對壓差阻力和摩擦阻力的影響趨勢是相同的，這也解釋了在總阻力預測中，梯度求解方法影響最大的原因。

圖15 不同因素不同水平下阻力均值的變化曲線

上面的分析中沒有考慮各個因素之間的交互作用。從枚舉法的結(jié)果中可以看到，兩個因素之間存在一定的交互作用，如圖16所示，可以明顯看出網(wǎng)格和湍流模型之間存在一定的交互作用，這在統(tǒng)計學上稱為正交互。

為了進一步考查各因素之間的交互作用貢獻，根據(jù)混合水平正交表L18(2×37)安排試驗，如表6所示，第3、6、7、8列作為誤差列。需要說明的是，混合水平正交表L18(2×37)可以考察置于第1列的二水平因子與置于第2列的三水平因子的交互作用，但不能考察其他列間的交互作用，故

圖16 湍流模型和網(wǎng)格因素的交互作用

表6 考慮兩因素交互作用的擬水平法正交試驗設計

此處只求解湍流模型和網(wǎng)格的交互作用，梯度求解方法與網(wǎng)格、湍流模型間的交互作用可用相同方法進行分析。

對此可以通過衡量各因素的貢獻率來評估因素的重要性，貢獻率的詳細推導過程可以查閱文獻[31]，各因素對阻力預測的貢獻率分析結(jié)果如表7所示。

通過正交表結(jié)果分析，可以得到以下結(jié)論：

1) 貢獻率分析的結(jié)果與極差分析保持一致：對壓差阻力，湍流模型對指標影響最大，其次是網(wǎng)格；對摩擦阻力，湍流模型對指標影響最大；對總阻力，梯度求解法對指標影響最大，這也與枚舉法的敏感性分析結(jié)果一致。

2) 梯度求解方法與網(wǎng)格、梯度求解方法與湍流模型的交互作用對壓差阻力、摩擦阻力以及總阻力的貢獻率均可忽略；湍流模型與網(wǎng)格的交互作用引起的數(shù)據(jù)波動只對總阻力有較顯著的貢獻，但也與誤差引起的數(shù)據(jù)波動的貢獻率差不多，故也可忽略。

表7 阻力預測的貢獻率分析

3 結(jié) 論

使用枚舉法來考察所有因素的綜合作用，通過基于聚類分析的定性敏感性分析方法和基于Mckay主影響分析的定量方法，識別出了對總阻力預測影響較大的因素為梯度求解方法。而湍流模型是影響壓差阻力和摩擦阻力預測的最關(guān)鍵因素。

為了減輕計算負擔，使用正交試驗設計方法完成上述研究，敏感性分析結(jié)果與枚舉法基本一致，但計算量大大減少，這位解決類似的因素影響分析提供了切實可行的思路。

需要指出的是，本文得到的模擬因素敏感性結(jié)論都是基于使用的這幾種網(wǎng)格、模型和格式，對這一特定算例分析得到的。如果使用不同的模擬方法，比如其他的網(wǎng)格、模型和格式，或者針對新的算例，可能會得到不一致的結(jié)論。但本文給出了一般的因素敏感性分析方法。針對具體問題，可以仿照這一流程進行分析。