數值模擬方法對NASA CRM模型阻力預測的影響

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

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

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

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

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

1 計算外形、網格和方法

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

圖1 NASA CRM外形

圖2 CRM不同風洞試驗數據對比

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

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

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

圖3 粗網格表面網格

表1 計算使用的網格信息

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

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

2 阻力預測影響分析

2.1 枚舉法結果

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

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

表2 計算安排設置

圖5 阻力系數預測值的統計結果

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

圖6 阻力預測的網格收斂特性曲線

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

圖8 俯仰力矩系數預測的網格收斂特性曲線

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

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

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

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

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

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

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

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

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

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

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

2.3 基于McKay的敏感性分析

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

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

(1)

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

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

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

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

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

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

表3 McKay主影響指標的參數敏感性

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

2.4 正交試驗設計分析

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

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

(2)

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

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

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

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

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

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

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

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

圖16 湍流模型和網格因素的交互作用

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

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

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

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

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

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

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

3 結 論

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

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

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