幾何不確定性區(qū)間分析及魯棒氣動(dòng)優(yōu)化設(shè)計(jì)

宋鑫， 鄭冠男,*， 楊國(guó)偉， 姜倩

(1. 中國(guó)科學(xué)院力學(xué)研究所流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室, 北京 100080；2. 中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院, 北京 100049)

航空航天工業(yè)中，飛行器設(shè)計(jì)的各參數(shù)往往被看作確定的量，但實(shí)際上，工程中是無(wú)法避免不確定性存在的。例如飛行器結(jié)構(gòu)設(shè)計(jì)由于制造水平、測(cè)量誤差等會(huì)造成參數(shù)的不確定性，飛行器飛行環(huán)境及載荷作用也會(huì)存在更強(qiáng)的不確定性等[1]。不確定性因素嚴(yán)重影響了與其相對(duì)的確定性設(shè)計(jì)的精度和可靠度，因此，想要獲得更可靠更穩(wěn)定的設(shè)計(jì)，必須考慮不確定性因素的影響。考慮不確定性的設(shè)計(jì)方法主要包括魯棒(穩(wěn)健)優(yōu)化設(shè)計(jì)方法[2]及可靠性優(yōu)化設(shè)計(jì)方法[3]。魯棒優(yōu)化設(shè)計(jì)旨在降低飛行器性能對(duì)不確定因素的靈敏度，而可靠性設(shè)計(jì)旨在降低發(fā)生故障或失效的概率，提高飛行器的可靠性。本文主要研究翼型的魯棒優(yōu)化設(shè)計(jì)方法。

與傳統(tǒng)的確定性優(yōu)化相比，魯棒優(yōu)化設(shè)計(jì)需要在優(yōu)化過(guò)程中不斷進(jìn)行不確定性分析。其主要包括概率方法和非概率方法。概率方法又包括蒙特卡羅(Monte Carlo)方法、泰勒展開、隨機(jī)展開等，如徐明等[4]采用蒙特卡羅方法對(duì)旋翼的最優(yōu)轉(zhuǎn)速進(jìn)行了不確定性分析；Salehi等[5]和鄔曉敬等[6]采用非浸入式混沌多項(xiàng)式方法對(duì)離心泵的流場(chǎng)及翼型跨聲速隨機(jī)氣動(dòng)特性等分別進(jìn)行了不確定性分析；戴玉婷和楊超[7]考慮廣義剛度的隨機(jī)不確定性，基于浸入式隨機(jī)展開方法分析了不確定性對(duì)顫振邊界的影響；Papadimitriou等[8]考慮流動(dòng)相關(guān)參數(shù)及幾何不確定性，并采用離散網(wǎng)格技術(shù)進(jìn)行不確定性分析，對(duì)阻力系數(shù)目標(biāo)進(jìn)行了魯棒優(yōu)化設(shè)計(jì)。概率方法一般要事先給定不確定性變量的概率分布函數(shù)，但是工程中很難獲得充足的樣本來(lái)估計(jì)分布函數(shù)的參數(shù)或驗(yàn)證概率分布函數(shù)的合理性，而人為假設(shè)又會(huì)導(dǎo)致較大的分析誤差[9]，相比之下獲得不確定性因素的取值范圍較容易，因此研究人員發(fā)展了“基于邊界表征”的非概率不確定性方法，其中區(qū)間分析方法應(yīng)用較多。屈小章等[10]采用區(qū)間方法描述風(fēng)機(jī)翼型的幾何不確定性變量，對(duì)其氣動(dòng)性能進(jìn)行了穩(wěn)健優(yōu)化設(shè)計(jì)；Zheng和Qiu[11]采用區(qū)間分析的方法對(duì)考慮不確定性的氣動(dòng)力和氣動(dòng)熱進(jìn)行了分析；張軍紅和韓景龍[12]提出了一種考慮區(qū)間不確定性的機(jī)翼顫振優(yōu)化方法。

飛行器生產(chǎn)時(shí)本身存在制造誤差，在飛行過(guò)程中也會(huì)由于載荷的作用以及結(jié)冰、燒蝕等現(xiàn)象使飛行器氣動(dòng)外形產(chǎn)生局部變化，進(jìn)而影響氣動(dòng)性能。另外，飛行器跨聲速飛行時(shí)，流場(chǎng)具有強(qiáng)非線性，幾何外形變化對(duì)氣動(dòng)性能的影響會(huì)更大、更復(fù)雜。以往研究中，對(duì)幾何不確定性變量或采用假設(shè)的概率方法描述[8,13]，容易造成較大的分析誤差；或采用基于特征幾何參數(shù)的區(qū)間描述[10-11]，無(wú)法表達(dá)翼型任意的局部變化。且大多數(shù)非概率方法研究均基于區(qū)間數(shù)運(yùn)算分析方法，該方法適用于分析函數(shù)在不確定區(qū)間內(nèi)為單調(diào)函數(shù)的情況，不適用于復(fù)雜的非線性氣動(dòng)問(wèn)題。

針對(duì)以上問(wèn)題，本文研究了幾何不確定性的通用區(qū)間描述及快速非線性區(qū)間分析方法。在此基礎(chǔ)上建立了魯棒優(yōu)化設(shè)計(jì)流程，并通過(guò)標(biāo)準(zhǔn)算例驗(yàn)證了方法的有效性。

1 區(qū)間不確定性分析

1.1 直接操作自由變形參數(shù)化方法

采用區(qū)間方法進(jìn)行幾何不確定性分析，首先要通過(guò)合適的參數(shù)化方法建立翼型幾何不確定性的區(qū)間描述。工程中現(xiàn)已發(fā)展出多種翼型參數(shù)化方法[14]，如類別形狀函數(shù)變換(Class Shape Transformation, CST)方法、Hicks-Henne型函數(shù)法、PARSEC方法、B樣條法等。上述方法或通過(guò)間接變量描述翼型，或通過(guò)如前緣半徑、最大厚度位置等特征參數(shù)描述，但均無(wú)法直接控制翼型上某點(diǎn)的變化，不能靈活表達(dá)翼型表面凹凸等局部變形，故不適合描述任意的幾何不確定性。本文采用直接操作自由變形(DFFD)方法[15]對(duì)翼型參數(shù)化，以翼型表面點(diǎn)的直接變形為參數(shù)建立幾何不確定性量區(qū)間描述。

對(duì)于一般自由變形(FFD)方法，首先在變形物體周圍布置控制體，如圖1(a)所示，圖中，x/c為弦向位置,y/c為縱向位置，x、y分別為實(shí)際弦向及縱向坐標(biāo)，c為翼型弦長(zhǎng)。然后定義變形物體上各點(diǎn)實(shí)際坐標(biāo)關(guān)于局部坐標(biāo)的函數(shù)為

圖1 FFD方法及DFFD方法控制翼型變形Fig.1 Airfoil deformations controlled by FFD and DFFD methods

(1)

在定義的局部坐標(biāo)系中，若FFD控制點(diǎn)的位移量為ΔPi,j,k，相應(yīng)可得變形物體上各點(diǎn)的位移量為

(2)

則變形物體上各點(diǎn)的最終位置為

x′(s,t,u)=x(s,t,u)+Δx(s,t,u)

(3)

FFD方法雖然變形簡(jiǎn)單，但變形物面的點(diǎn)和FFD控制點(diǎn)并沒(méi)有直觀的聯(lián)系，無(wú)法通過(guò)FFD控制點(diǎn)的變形精確控制物面某點(diǎn)的變形。DFFD方法在FFD方法的基礎(chǔ)上發(fā)展而來(lái)，其不同于一般FFD方法的是除在翼型周圍布置控制體和控制點(diǎn)外，在翼型上也布置直接操作點(diǎn)(pilot points)。通過(guò)直接操作點(diǎn)的位移反求FFD控制體上控制點(diǎn)的位移，進(jìn)而通過(guò)FFD方法控制翼型的變化。由于將直接操作點(diǎn)的位移作為控制參數(shù)，因此可以使用高階的FFD控制體而不增加變量維數(shù)。反求的過(guò)程可看作求解滿足直接操作點(diǎn)位移約束的FFD控制點(diǎn)的位移組合，假設(shè)直接操作點(diǎn)的位移相互獨(dú)立，Sf(f=1,2,…,h) 為變形前直接操作點(diǎn)的原始坐標(biāo)，Tf(f=1,2,…,h)為變形后各直接操作點(diǎn)的坐標(biāo)。由式(2)可得

(4)

式中：δi,j,k為待求的FFD控制點(diǎn)的位移。

1.2 非線性區(qū)間分析方法

區(qū)間分析方法[16]采用區(qū)間數(shù)學(xué)的概念，將每一維不確定量v表示為區(qū)間形式：

AI=[AL,AR]={v|AL≤v≤AR,v∈R}

(5)

式中：上標(biāo)I、L、R分別表示區(qū)間、區(qū)間下界、區(qū)間上界。由區(qū)間2個(gè)端點(diǎn)組成的一對(duì)有序?qū)崝?shù)稱為區(qū)間數(shù)。區(qū)間數(shù)的中值記作AC=(AL+AR)/2，區(qū)間半徑或離差記作Aw=(AR-AL)/2。

普通的線性區(qū)間函數(shù)問(wèn)題可以通過(guò)區(qū)間數(shù)運(yùn)算以及區(qū)間擴(kuò)張獲得關(guān)于自變量的區(qū)間函數(shù)值，但跨聲速條件下的氣動(dòng)問(wèn)題屬于復(fù)雜的非線性問(wèn)題，且沒(méi)有顯式的函數(shù)關(guān)系。本文通過(guò)優(yōu)化方法進(jìn)行非線性區(qū)間分析，即采用單目標(biāo)粒子群進(jìn)化算法，分別進(jìn)行最大化目標(biāo)、最小化目標(biāo)共兩次尋優(yōu)獲得目標(biāo)函數(shù)區(qū)間的上界和下界。

除采用直接優(yōu)化進(jìn)行不確定分析外，本文同時(shí)采用構(gòu)建代理模型結(jié)合優(yōu)化方法進(jìn)行不確定分析，以提高效率。本文選擇具有優(yōu)異的非線性函數(shù)近似能力的Kriging模型[17]作為代理模型，建立不確定變量與氣動(dòng)性能參數(shù)的函數(shù)關(guān)系。即通過(guò)拉丁超立方設(shè)計(jì)在不確定性變量取值空間選取一定量樣本并進(jìn)行CFD計(jì)算，通過(guò)樣本集建立Kriging模型，對(duì)模型分別進(jìn)行最大最小化目標(biāo)尋優(yōu)。為進(jìn)一步提高代理模型精度，在初始樣本的基礎(chǔ)上采用最大均方差準(zhǔn)則(MSE)進(jìn)行加點(diǎn)，即選取當(dāng)前模型在建模空間內(nèi)預(yù)測(cè)均方差最大的點(diǎn)進(jìn)行精確CFD分析并加入到樣本集內(nèi)，重新構(gòu)建代理模型。

1.3 不確定性分析算例

本文針對(duì)NACA0012翼型，進(jìn)行了考慮幾何不確定性的阻力性能分析。為簡(jiǎn)化問(wèn)題，考慮不確定性時(shí)仍保證翼型對(duì)稱。首先采用DFFD方法對(duì)上半翼型進(jìn)行參數(shù)化，建立10×8階FFD控制體，并在翼型上選取7個(gè)直接操作點(diǎn)，其弦向位置如表1所示。為保證前后緣固定，首尾2個(gè)直接操作點(diǎn)固定，改變其余5個(gè)直接操作點(diǎn)的y向位移以控制翼型變化，因此共有5個(gè)不確定性變量，其變化范圍均為[-0.001 5,0.001 5]。采用粒子群算法優(yōu)化求解時(shí)，除粒子的位置和速度約束在變化范圍內(nèi)，無(wú)其他約束。

翼型的阻力預(yù)測(cè)采用課題組自主研發(fā)的基于Navier-Stokes方程的并行CFD求解器[18]，湍流模型采用k-ω兩方程模型，空間離散采用二階格式，時(shí)間推進(jìn)采用LU-SGS方法，網(wǎng)格采用美國(guó)國(guó)家航空航天局(NASA)提供的449×129標(biāo)準(zhǔn)網(wǎng)格。圖2為馬赫數(shù)0.7、迎角4°、雷諾數(shù)4.06×106條件下，計(jì)算壓力系數(shù)Cp分布與試驗(yàn)數(shù)據(jù)[6]的比較，計(jì)算精度可滿足要求。

不確定分析算例采用工況為迎角0°，馬赫數(shù)0.85，該條件下計(jì)算得阻力系數(shù)為0.0538。首先采用20個(gè)初始樣本建立了阻力系數(shù)與不確定變量的Kriging模型(Kriging模型1)。在此基礎(chǔ)上又加點(diǎn)7次，重新構(gòu)建模型(Kriging模型2)。在不確定性變化區(qū)間內(nèi)隨機(jī)選取80個(gè)樣本使用模型進(jìn)行預(yù)測(cè)并采用精確CFD分析進(jìn)行驗(yàn)證，以檢驗(yàn)代理模型精度，樣本阻力系數(shù)的相對(duì)誤差如圖3所示，采用20個(gè)初始樣本建立代理模型，相對(duì)誤差最大為5.47%，加點(diǎn)7次后，最大相對(duì)誤差減小為2.24%，該精度可以滿足工程應(yīng)用。

表1 直接操作點(diǎn)x方向位置Table 1 Position of pilot points in x direction

圖2 NACA0012翼型計(jì)算與試驗(yàn)壓力分布比較Fig.2 Comparison of pressure distribution between computation and experiment for NACA0012 airfoil

為驗(yàn)證不確定性分析的準(zhǔn)確性，采用蒙特卡羅方法進(jìn)行不確定性分析，即在設(shè)計(jì)空間內(nèi)選取200個(gè)隨機(jī)樣本進(jìn)行精確CFD分析并統(tǒng)計(jì)，并與上述3種結(jié)果對(duì)比，直接優(yōu)化、Kriging模型優(yōu)化與蒙特卡羅方法結(jié)果比較如圖4所示。結(jié)果表明直接優(yōu)化方法具備較高的精度，只通過(guò)20個(gè)初始樣本構(gòu)建Kriging模型，得到的結(jié)果與直接優(yōu)化結(jié)果存在一定誤差，通過(guò)加點(diǎn)重建Kriging模型，誤差進(jìn)一步減小。直接優(yōu)化中每一步所有CFD分析可并行進(jìn)行，建立代理模型的20個(gè)初始樣本也可并行進(jìn)行CFD分析，但加點(diǎn)過(guò)程無(wú)法并行進(jìn)行，故3種方法的分析結(jié)果、相對(duì)誤差、CFD計(jì)算次數(shù)、并行計(jì)算時(shí)間比較如表2所示。其中，采用Kringing模型1最大相對(duì)誤差小于5%，而分析效率可提高95%。3種方法得到的阻力系數(shù)變化區(qū)間上下界對(duì)應(yīng)的翼型及壓力分布如圖5所示，由圖中可以看出，各方法得到的翼型相似。

圖3 Kriging模型阻力系數(shù)相對(duì)誤差Fig.3 Drag relative coefficient errors of Kriging models

圖4 直接優(yōu)化、Kriging模型優(yōu)化與蒙特卡羅方法結(jié)果比較Fig.4 Comparison of results among direct optimization, Kriging model optimization and Monte Carle method

表2 3種方法的分析結(jié)果、誤差、CFD計(jì)算次數(shù)、計(jì)算時(shí)間比較Table 2 Comparison of analysis results, errors, CFD calculation times and computing time among three methods

圖5 阻力系數(shù)變化區(qū)間上下界對(duì)應(yīng)的翼型及壓力分布Fig.5 Airfoils and pressure distribution corresponding to upper and lower bounds of drag coefficient variation interval

2 魯棒優(yōu)化設(shè)計(jì)

2.1 優(yōu)化算例

本文采用RAE2822翼型優(yōu)化設(shè)計(jì)標(biāo)準(zhǔn)算例，流場(chǎng)求解采用前述的CFD求解器。優(yōu)化問(wèn)題描述為

(6)

式中：CD為阻力系數(shù)目標(biāo)；X為描述翼型幾何的n維設(shè)計(jì)變量；V為設(shè)計(jì)空間；考慮升力系數(shù)CL、俯仰力矩系數(shù)CM以及面積Aa約束，Ainitial為RAE2822翼型的面積，大小為0.077 873。計(jì)算工況為馬赫數(shù)0.734，雷諾數(shù)6.5×106。圖6為迎角2.79°條件下該翼型的計(jì)算與試驗(yàn)壓力分布比較結(jié)果。

首先對(duì)該問(wèn)題進(jìn)行傳統(tǒng)的確定性優(yōu)化。通過(guò)試驗(yàn)，采用DFFD方法對(duì)翼型進(jìn)行參數(shù)化時(shí)，若參數(shù)的變化范圍較大，極易使翼型產(chǎn)生大的波動(dòng)，從而不具備翼型的基本幾何特征，類似問(wèn)題也出現(xiàn)在其他直接以物面點(diǎn)位移參數(shù)化的方法中[19]。因此，DFFD方法雖可以直接控制翼型表面變形，但并不適合較大空間內(nèi)的翼型參數(shù)化建模，故本文采用CST方法對(duì)翼型進(jìn)行參數(shù)化，作為確定性的設(shè)計(jì)變量，翼型上下部分分別采用5階CST參數(shù)化。由于進(jìn)行CFD計(jì)算時(shí)直接改變迎角進(jìn)行定升力計(jì)算，因此假設(shè)升力系數(shù)等式約束直接成立，實(shí)現(xiàn)時(shí)，給定閾值[0.823,0.825]，檢測(cè)升力系數(shù)穩(wěn)定收斂到閾值內(nèi)結(jié)束計(jì)算。不等式約束處理采用罰函數(shù)法，約束函數(shù)值計(jì)算式為

(7)

優(yōu)化算法采用單目標(biāo)粒子群算法，優(yōu)化結(jié)果以及RAE2822翼型計(jì)算結(jié)果如表3所示，α為迎角。對(duì)本文的優(yōu)化結(jié)果與其他文獻(xiàn)[19-21]進(jìn)行了對(duì)比，數(shù)值結(jié)果如表4所示，圖7為優(yōu)化后的翼型對(duì)比，圖8為優(yōu)化后壓力分布對(duì)比，本文優(yōu)化后阻力系數(shù)下降35.9%，且滿足各項(xiàng)約束。與其他文獻(xiàn)結(jié)果對(duì)比，在基本相同的設(shè)計(jì)變量個(gè)數(shù)下，本文取得了較好的結(jié)果。

圖6 RAE2822翼型計(jì)算與試驗(yàn)壓力分布比較Fig.6 Comparison of pressure distribution between computation and experiment for RAE2822 airfoil

表3 RAE2822翼型及優(yōu)化翼型的計(jì)算結(jié)果Table 3 Computing results of RAE2822 and optimized airfoils

表4 優(yōu)化翼型與其他文獻(xiàn)結(jié)果對(duì)比

圖7 優(yōu)化后翼型對(duì)比Fig.7 Comparsion of optimized airfoils

圖8 優(yōu)化后翼型壓力分布對(duì)比Fig.8 Pressure distribution comparison of optimized airfoils

2.2 魯棒優(yōu)化設(shè)計(jì)方法

在2.1節(jié)算例的基礎(chǔ)上考慮翼型幾何的區(qū)間不確定性，問(wèn)題(6)可描述為

(8)

式中：U為q維不確定向量，用q維區(qū)間向量UI表示。目標(biāo)函數(shù)和約束函數(shù)均為不確定變量U的連續(xù)函數(shù)，故其可能取值也構(gòu)成區(qū)間形式，而不是確定值，因此，上述問(wèn)題無(wú)法通過(guò)傳統(tǒng)的確定性優(yōu)化方法求解。

魯棒優(yōu)化設(shè)計(jì)中，將不確定性變量導(dǎo)致的目標(biāo)性能區(qū)間作為設(shè)計(jì)目標(biāo)，因此在優(yōu)化過(guò)程中，需比較不同樣本的目標(biāo)性能區(qū)間的大小。本文中，利用≤Cw區(qū)間序關(guān)系[16]比較2個(gè)區(qū)間AI和BI的大小，對(duì)于最小化問(wèn)題其定義如下：

(9)

P(AI≤BI)=

(10)

(11)

如2.1節(jié)所述，升力系數(shù)等式約束自動(dòng)滿足，其他約束采用罰函數(shù)法處理，則問(wèn)題(11)進(jìn)一步轉(zhuǎn)換為無(wú)約束優(yōu)化問(wèn)題：

(12)

式中：

本文基于區(qū)間不確定性分析方法及自適應(yīng)多目標(biāo)粒子群算法，建立了考慮幾何區(qū)間不確定性的魯棒優(yōu)化設(shè)計(jì)方法，求解上述優(yōu)化問(wèn)題，其流程如圖9所示。

2.3 自適應(yīng)多目標(biāo)粒子群算法

本文采用粒子群算法具有形式簡(jiǎn)潔、收斂快速以及參數(shù)調(diào)節(jié)機(jī)制靈活等優(yōu)點(diǎn)。在粒子群算法中，一個(gè)粒子i可由位置向量xi=[xi,1,xi,2,…,xi,d]T∈Rd和速度向量vi=[vi,1,vi,2,…,vi,d]T∈Rd表示，其中i=1,2,…,N，N為種群中的粒子個(gè)數(shù)，d為設(shè)計(jì)變量的維數(shù)。粒子在進(jìn)化過(guò)程中根據(jù)式(13)更新速度和位置：

(13)

式中：t為迭代次數(shù)；ω≥0為慣性權(quán)重；c1、c2為加速系數(shù)；r1、r2為在[0,1]上均勻分布的隨機(jī)數(shù)；pbesti為第i個(gè)粒子的歷史個(gè)體最優(yōu)解;gbest為群體的全局最優(yōu)解。為避免陷入局部最優(yōu)，可采用精英學(xué)習(xí)策略進(jìn)行局部極值擾動(dòng)，精英學(xué)習(xí)率為L(zhǎng)r。本文自適應(yīng)多目標(biāo)粒子群算法基于Pareto占優(yōu)原則建立并維護(hù)外部最優(yōu)解集，通過(guò)引入Pareto熵表示Pareto前緣的分布均勻性，間接表示種群的多樣性。優(yōu)化迭代時(shí)，Pareto前緣發(fā)生變化，相應(yīng)的Pareto熵也發(fā)生變化，故可以采用差熵表示Pareto前緣重新分布范圍的大小，通過(guò)差熵可以推測(cè)當(dāng)前種群發(fā)現(xiàn)新解的能力，估計(jì)種群所處的進(jìn)化狀態(tài)包括收斂狀態(tài)、多樣化狀態(tài)及停滯狀態(tài)，進(jìn)而調(diào)整參數(shù)ω、c1、c2、Lr以設(shè)計(jì)適應(yīng)當(dāng)前狀態(tài)的勘探和開采策略。

圖9 考慮幾何不確定性的魯棒優(yōu)化設(shè)計(jì)流程Fig.9 Robust optimization design process considering geometric uncertainties

優(yōu)化算法的具體流程如圖10所示，文獻(xiàn)[22]采用此優(yōu)化算法對(duì)不同測(cè)試算例進(jìn)行了驗(yàn)證，并與其他常用多目標(biāo)優(yōu)化算法進(jìn)行了對(duì)比，結(jié)果顯示該算法總體上具有更好的多樣性和收斂性。

2.4 結(jié)果分析

采用2.2節(jié)魯棒優(yōu)化設(shè)計(jì)方法對(duì)問(wèn)題(12)優(yōu)化。幾何不確定性變量仍采用DFFD方法參數(shù)化產(chǎn)生，對(duì)超臨界翼型的上半部分和下半部分分別進(jìn)行參數(shù)化，各布置7個(gè)直接操作點(diǎn)，采用16×12階控制體，直接操作點(diǎn)及控制體建立如圖11所示。翼型前緣和后緣的直接操作點(diǎn)仍保持固定，故共10個(gè)不確定性變量，變化區(qū)間仍為[-0.0015 ,0.001 5]。為提高計(jì)算效率，對(duì)計(jì)算精度與計(jì)算時(shí)間進(jìn)行折中，選取通過(guò)25個(gè)初始樣本構(gòu)建Kriging模型的方式進(jìn)行阻力系數(shù)及力矩系數(shù)的不確定性分析，直接優(yōu)化的方式進(jìn)行翼型面積的不確定分析。代理模型精度校驗(yàn)得阻力系數(shù)預(yù)測(cè)的最大誤差為7.4%，力矩系數(shù)預(yù)測(cè)的最大誤差為0.63%。多目標(biāo)優(yōu)化所得Pareto解集如圖12所示。

由于多目標(biāo)優(yōu)化中粒子群容量及進(jìn)化步數(shù)設(shè)置較小，且跨聲速問(wèn)題機(jī)理復(fù)雜，所得的Pareto前緣不甚明顯，但工程中只需選取一個(gè)或幾個(gè)最優(yōu)解即可，故選取圖中Pareto前緣的一點(diǎn)作為本次魯棒優(yōu)化的最優(yōu)解，并對(duì)最優(yōu)翼型進(jìn)行CFD分析，得其阻力系數(shù)CD；同時(shí)采用相同的參數(shù)化及代理模型建立方法對(duì)確定性優(yōu)化得到的最優(yōu)翼型及初始RAE2822翼型進(jìn)行一次不確定性分析，其相關(guān)數(shù)值如表5所示。2種優(yōu)化方法得到的最優(yōu)翼型及對(duì)應(yīng)壓力分布如圖13所示。對(duì)RAE2822翼型及最優(yōu)翼型進(jìn)行不確定性分析，得到阻力目標(biāo)、力矩與面積約束區(qū)間上下界對(duì)應(yīng)的翼型及壓力分布，如圖14所示。

圖11 直接操作點(diǎn)及FFD控制體Fig.11 Pilot points and FFD control body

圖12 多目標(biāo)魯棒優(yōu)化的Pareto解集Fig.12 Pareto set of multi-objective robust optimization

首先從優(yōu)化翼型的CFD計(jì)算結(jié)果可看出，無(wú)論是確定性優(yōu)化還是魯棒優(yōu)化最優(yōu)翼型的阻力系數(shù)較初始RAE2822翼型均有較大程度的下降。從翼型幾何和對(duì)應(yīng)壓力分布情況來(lái)看，優(yōu)化使得翼型上表面前緣下降后緣抬升，造成前緣吸力峰值更高，但減弱了激波強(qiáng)度，從而減小阻力。確定性優(yōu)化結(jié)果較魯棒優(yōu)化結(jié)果激波更平緩，阻力更小。從考慮阻力及其他約束的綜合不確定性分析的結(jié)果來(lái)看，原始RAE2822翼型以及確定性優(yōu)化的區(qū)間半徑目標(biāo)fw較魯棒優(yōu)化結(jié)果差距很大，但從表5也可看出,區(qū)間中值目標(biāo)fC與CFD計(jì)算結(jié)果CD相差較大，說(shuō)明不確定分析結(jié)果fw和fC較大主要是由約束函數(shù)造成的。

圖15為不確定分析時(shí)各采樣樣本的目標(biāo)函數(shù)和約束函數(shù)值比較，該結(jié)果一定程度上可顯示優(yōu)化目標(biāo)和各約束對(duì)不確定分析結(jié)果的影響。其中阻力系數(shù)結(jié)果顯示，確定性優(yōu)化結(jié)果的整體阻力系數(shù)較低且變化區(qū)間較小，但從力矩和面積約束的結(jié)果可看出，確定性優(yōu)化的結(jié)果均在一定程度上違反了約束，而魯棒優(yōu)化結(jié)果則基本滿足約束。這是因?yàn)榇_定性優(yōu)化時(shí)，由于目標(biāo)和約束的沖突，所得最優(yōu)結(jié)果一般分布在約束邊緣，當(dāng)考慮不確定性因素時(shí)，極易違反約束。由此可見(jiàn)，本文魯棒優(yōu)化流程考慮了不確定性對(duì)優(yōu)化目標(biāo)及約束的影響，相對(duì)確定性優(yōu)化，得到了綜合魯棒性更強(qiáng)的結(jié)果。實(shí)際操作中，也可以根據(jù)決策者對(duì)各優(yōu)化目標(biāo)和約束的偏好，靈活設(shè)定約束罰函數(shù)因子，得到更理想的結(jié)果。

表5 優(yōu)化結(jié)果比較Table 5 Optimization result comparison

圖13 優(yōu)化所得最優(yōu)翼型及壓力分布比較Fig.13 Optimized airfoils and pressure distribution comparison

圖14 目標(biāo)及約束區(qū)間上下界對(duì)應(yīng)的翼型及壓力分布Fig.14 Airfoils and pressure distribution corresponding to upper and lower bounds of objective and constraint intervals

圖15 不確定分析采樣樣本的阻力系數(shù)、力矩系數(shù)及面積對(duì)比Fig.15 Drag coefficient, moment coefficient and area comparison of samples for uncertainty analysis

3 結(jié) 論

1) 對(duì)NACA0012翼型的阻力系數(shù)性能進(jìn)行了跨聲速條件下考慮幾何不確定性的快速非線性區(qū)間分析。通過(guò)建模方法與直接分析方法對(duì)比，最大誤差在5%以內(nèi)，而分析效率可提高95%。

2) 以跨聲速條件下RAE2822翼型的阻力性能為目標(biāo)，并考慮升力、力矩及面積約束，分別進(jìn)行了翼型的確定性優(yōu)化設(shè)計(jì)及魯棒氣動(dòng)優(yōu)化設(shè)計(jì)。確定性優(yōu)化所得最優(yōu)翼型較初始翼型阻力減小35.9%，且嚴(yán)格滿足各項(xiàng)約束，同等精度下優(yōu)于其他文獻(xiàn)結(jié)果。魯棒優(yōu)化所得最優(yōu)翼型較初始翼型和確定優(yōu)化最優(yōu)翼型，其魯棒性得到了大幅提高。分析總體優(yōu)化目標(biāo)各項(xiàng)可知，本次優(yōu)化結(jié)果主要通過(guò)降低違反約束的概率提高魯棒性，通過(guò)調(diào)節(jié)阻力目標(biāo)與其他約束的罰函數(shù)因子，可得到偏好不同的結(jié)果。