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

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

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

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

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

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

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

圖1 NASA CRM外形

圖2 CRM不同風(fēng)洞試驗(yàn)數(shù)據(jù)對(duì)比

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

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

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

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

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

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

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

2 阻力預(yù)測(cè)影響分析

2.1 枚舉法結(jié)果

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

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

表2 計(jì)算安排設(shè)置

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

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

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

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

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

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

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

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

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

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

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

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

以總阻力為例，圖9給出了所有樣本通過4個(gè)因素分類的情況，并用黑色實(shí)線將各個(gè)因素的結(jié)果分為兩類。在圖9(a)中第2類包括了全部細(xì)網(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類大部分為L(zhǎng)S結(jié)果，起到了有效區(qū)分，所以梯度求解方法對(duì)于總阻力的影響最大。類似地，對(duì)于壓差阻力，湍流模型發(fā)揮了最大的影響，對(duì)于摩擦阻力，湍流模型發(fā)揮了最大的影響。

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

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

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

2.3 基于McKay的敏感性分析

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

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

(1)

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

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

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

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

雖然在壓差阻力和摩擦阻力預(yù)測(cè)中，湍流模型都是影響最大的因素，但是在總阻力預(yù)測(cè)中，湍流模型并不是影響最大的因素。這里考慮是因?yàn)橥牧髂Ｐ蛯?duì)壓差阻力和摩擦阻力的影響趨勢(shì)相反。具體來說，SA模型預(yù)測(cè)的壓差阻力較SSTk-ω模型預(yù)測(cè)偏小，但摩擦阻力預(yù)測(cè)的偏大，這樣導(dǎo)致湍流模型導(dǎo)致的總阻力變化較小。而梯度求解方法對(duì)于壓差阻力和摩擦阻力的影響是一致的，GN預(yù)測(cè)的壓阻和摩阻都比LS偏大，所以總阻力預(yù)測(cè)中，這一趨勢(shì)被放大。

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

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

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

2.4 正交試驗(yàn)設(shè)計(jì)分析

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

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

(2)

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

表4 正交試驗(yàn)設(shè)計(jì)的因素與水平因子

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

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

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

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

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

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

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

表6 考慮兩因素交互作用的擬水平法正交試驗(yàn)設(shè)計(jì)

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

對(duì)此可以通過衡量各因素的貢獻(xiàn)率來評(píng)估因素的重要性，貢獻(xiàn)率的詳細(xì)推導(dǎo)過程可以查閱文獻(xiàn)[31]，各因素對(duì)阻力預(yù)測(cè)的貢獻(xiàn)率分析結(jié)果如表7所示。

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

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

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

表7 阻力預(yù)測(cè)的貢獻(xiàn)率分析

3 結(jié) 論

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

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

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