陳藝夫，馬宇航，藍(lán)慶生，孫衛(wèi)平，史亞云，楊體浩，*，白俊強(qiáng)

1．西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

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

3．中航通飛華南飛機(jī)工業(yè)有限公司，珠海 519000

4．西安交通大學(xué) 航空學(xué)院，西安 710049

所有工程系統(tǒng)的性能都會(huì)受到某種程度不確定性的影響。由于受到加工制造或飛行條件及環(huán)境等諸多不確定性因素的影響，飛行器的性能容易出現(xiàn)較大變化。傳統(tǒng)的氣動(dòng)設(shè)計(jì)屬于確定性設(shè)計(jì)，有可能導(dǎo)致氣動(dòng)性能對(duì)不確定因素異常敏感，甚至?xí)?lái)一定的安全隱患。因此，在飛行器設(shè)計(jì)階段，關(guān)注不確定因素對(duì)飛行器性能的影響顯得尤為重要。

美國(guó)NASA 于1998年對(duì)工程上多學(xué)科設(shè)計(jì)領(lǐng)域的需求進(jìn)行了調(diào)查［1］，調(diào)查結(jié)果表明，需要提高飛機(jī)在多種不確定性因素下的不確定性能。2002年，NASA 又推出了有關(guān)航空航天領(lǐng)域不確定性下多學(xué)科優(yōu)化設(shè)計(jì)的規(guī)劃［2］，規(guī)劃中明確指出了飛行器不確定性下設(shè)計(jì)的難點(diǎn)和挑戰(zhàn)，為未來(lái)一段時(shí)間內(nèi)的研究奠定了基礎(chǔ)。規(guī)劃推出后，NASA 蘭利研究中心致力于對(duì)跨聲速翼型進(jìn)行不確定性優(yōu)化設(shè)計(jì)。從Huyse［3］和Walters［4］等對(duì)二維翼型考慮馬赫數(shù)不確定性的優(yōu)化設(shè)計(jì)研究中，可以看出確定性優(yōu)化和不確定性優(yōu)化之間存在較大差異，進(jìn)一步說(shuō)明了在翼型設(shè)計(jì)中考慮不確定性影響的重要性。

為開(kāi)展考慮不確定性的魯棒優(yōu)化設(shè)計(jì)，首先需要建立高效精確的不確定性分析方法。在不確定性分析中，很重要的一步是將不確定源的特征轉(zhuǎn)化為輸出變量（Quantity of Interest，QoI）不確定性估計(jì)的過(guò)程，稱之為不確定性的傳播。不確定性傳播方法可分為4 大類：解析法［5］、取樣法［6］、隨機(jī)響應(yīng)面模型法和代理模型法。

現(xiàn)代氣動(dòng)設(shè)計(jì)大量依賴于使用CFD 工具來(lái)預(yù)測(cè)所考慮的氣動(dòng)外形的性能，在減少了昂貴的風(fēng)洞試驗(yàn)的同時(shí)，效率因此得到顯著提升，設(shè)計(jì)周期進(jìn)一步縮短。Pelletier［7］和Luckring［8］等在其研究中闡述了CFD 模擬中不確定性的來(lái)源和分類。Walters 等［4］總結(jié)了針對(duì)CFD 中不確定性的分析方法。Schaefer 等［9］在其最新的綜述類文章中對(duì)當(dāng)前航空航天領(lǐng)域CFD 模擬的不確定性分析方法進(jìn)行了總結(jié)。Roelofs等［10］在其文章中總結(jié)了在不確定性條件下所進(jìn)行的優(yōu)化設(shè)計(jì)研究工作，重點(diǎn)關(guān)注的是分配和量化這些不確定性的方法和建模方法。研究發(fā)現(xiàn)，概率方法仍然是表示不確定性的最流行理論。多項(xiàng)式混沌展開(kāi)法和隨機(jī)配置方法越來(lái)越受歡迎，但蒙特卡洛模擬仍被廣泛使用。

近年來(lái)，多項(xiàng)式混沌法在CFD 氣動(dòng)設(shè)計(jì)不確定性分析中的應(yīng)用最為廣泛。Knio 等［11］和Najm［12］在其綜述中總結(jié)了多項(xiàng)式混沌法在CFD不確定性分析中的發(fā)展和應(yīng)用現(xiàn)狀。國(guó)內(nèi)外學(xué)者采用多項(xiàng)式混沌法對(duì)CFD 氣動(dòng)設(shè)計(jì)中的不確定性問(wèn)題開(kāi)展了大量研究工作?；诜乔度胧蕉囗?xiàng)式混沌法，Loeven 等［13］研究了一維活塞問(wèn)題和二維亞聲速翼型的不確定性問(wèn)題。Simon［14］和Chassaing［15］等對(duì)二維翼型跨聲速下馬赫數(shù)和迎角的不確定影響進(jìn)行研究，敏感性分析表明，不確定參數(shù)之間存在很強(qiáng)的非線性耦合。Yamazaki［16］則利用阻力分解方法分析了二維翼型的不確定性量化問(wèn)題。在Suga 等［17］的研究中，多項(xiàng)式混沌方法被用于二維超聲速雙翼飛機(jī)翼型的氣動(dòng)不確定性量化。Hosder 等［18］分別對(duì)膨脹波拐角角度的幾何不確定性和跨聲速三維機(jī)翼的馬赫數(shù)、迎角不確定性開(kāi)展了研究。West 等［19］開(kāi)發(fā)了一種多保真多項(xiàng)式混沌法，并將其成功應(yīng)用于超聲速客機(jī)的不確定性分析問(wèn)題中。

在內(nèi)流研究中該方法也已得到了廣泛應(yīng)用。Ghisu［20］、Abraham［21］和Li［22］等利用非嵌入式多項(xiàng)式混沌法研究了燃?xì)鉁u輪壓縮機(jī)在不確定工作條件下的響應(yīng)。Liu 等［23］對(duì)渦輪葉片在風(fēng)速不確定下的總體性能進(jìn)行了討論。di Stefano 等［24］研究了湍流模型閉合系數(shù)不確定性對(duì)超燃沖壓發(fā)動(dòng)機(jī)流場(chǎng)解的影響。Guo 等［25］通過(guò)任意多項(xiàng)式混沌的稀疏近似來(lái)評(píng)估制造可變性對(duì)壓縮機(jī)葉片氣動(dòng)性能的影響。Huang等［26］提出了一種結(jié)合非嵌入式多項(xiàng)式混沌法和Smolyak 稀疏網(wǎng)格的不確定性量化方法，對(duì)轉(zhuǎn)子葉尖的氣動(dòng)和傳熱性能不確定性進(jìn)行量化。

國(guó)內(nèi)方面，王曉東等采用嵌入式多項(xiàng)式混沌方法分別研究了Burgers 方程的隨機(jī)響應(yīng)［27］和隨機(jī)方腔流動(dòng)中不確定性因素的影響［28］。鄔曉敬等［29］采用非嵌入式多項(xiàng)式混沌法對(duì)翼型跨聲速隨機(jī)氣動(dòng)特性進(jìn)行不確定性及全局敏感性分析。另外，對(duì)多項(xiàng)式混沌法的改進(jìn)也在持續(xù)進(jìn)行。Zhao等［30］提出了一種高效且新穎的自適應(yīng)多保真稀疏多項(xiàng)式混沌耦合Kriging 方法，該方法具有靈活性高，非線性建模能力強(qiáng)的優(yōu)點(diǎn)。在Schaefer的最新研究［31］中提出了一種空間精確多項(xiàng)式混沌的新統(tǒng)計(jì)方法，該方法可以在整個(gè)飛行包線或設(shè)計(jì)空間內(nèi)插補(bǔ)或外推偶然和認(rèn)知不確定性。盡管已經(jīng)有大量學(xué)者利用多項(xiàng)式混沌法開(kāi)展不確定性分析研究，但不確定性分析中所需的計(jì)算成本仍然較高，尤其是當(dāng)不確定性變量維度和多項(xiàng)式階數(shù)增大時(shí)。

除了以上所述對(duì)不確定性進(jìn)行分析和量化的研究，降低氣動(dòng)特性對(duì)不確定性變量的敏感性也是氣動(dòng)優(yōu)化設(shè)計(jì)領(lǐng)域的一大挑戰(zhàn)，這需要開(kāi)展基于不確定性下的魯棒優(yōu)化設(shè)計(jì)。

大量研究工作已經(jīng)基于多項(xiàng)式混沌法展開(kāi)。根據(jù)魯棒優(yōu)化設(shè)計(jì)中是否使用梯度，可分為 PCE非梯度類魯棒優(yōu)化設(shè)計(jì)、PCE 梯度類魯棒優(yōu)化設(shè)計(jì)?；诜翘荻阮悆?yōu)化方法，Zein［32］結(jié)合多項(xiàng)式混沌法提出了一種信任區(qū)域優(yōu)化方法，在信任區(qū)域上建立代理模型以完成魯棒優(yōu)化設(shè)計(jì)。Wang等［33］耦合了非嵌入式多項(xiàng)式混沌法與多目標(biāo)遺傳算法，并將其應(yīng)用于壓縮機(jī)轉(zhuǎn)子葉片的魯棒優(yōu)化設(shè)計(jì)中。Palar 等［34］將進(jìn)化算法和非嵌入式多項(xiàng)式混沌法結(jié)合，開(kāi)展了跨聲速翼型的魯棒優(yōu)化設(shè)計(jì)，目標(biāo)為降低不確定性并最大化升阻比。國(guó)內(nèi)方面，鄔曉敬［35］耦合非嵌入式多項(xiàng)式混沌法和Kriging 代理模型發(fā)展了自適應(yīng)魯棒優(yōu)化設(shè)計(jì)方法，對(duì)二維翼型進(jìn)行了考慮不確定性的優(yōu)化設(shè)計(jì)。

使用梯度優(yōu)化算法仍然是解決大規(guī)模設(shè)計(jì)變量氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題最有效的途徑之一，對(duì)于考慮不確定性的優(yōu)化設(shè)計(jì)其收益將更為明顯。Padulo 等［36］提出了一種簡(jiǎn)化正交技術(shù)，可在不顯著增加計(jì)算成本的條件下用于考慮較少設(shè)計(jì)變量的翼型氣動(dòng)優(yōu)化設(shè)計(jì)中。Shankaran 等［37］提出了一個(gè)魯棒優(yōu)化框架，結(jié)合多項(xiàng)式混沌法和伴隨理論來(lái)有效地獲得梯度，以實(shí)現(xiàn)考慮馬赫數(shù)不確定性的魯棒優(yōu)化設(shè)計(jì)。Schillings 等［38］基于多項(xiàng)式混沌和梯度算法對(duì)二維歐拉流動(dòng)不確定性下的氣動(dòng)外形開(kāi)展魯棒優(yōu)化設(shè)計(jì)研究，目標(biāo)函數(shù)定義為均值和方差的組合。Rallabhandi 等［39］使用PCE方法將伴隨導(dǎo)出的梯度與不確定度量化相耦合來(lái)同時(shí)對(duì)音爆進(jìn)行魯棒優(yōu)化，同時(shí)降低音爆響度的均值和標(biāo)準(zhǔn)差。Boopathy 等［40］提出了一種半嵌入式隨機(jī)Galerkin 方法，該方法可以將確定性導(dǎo)數(shù)用于伴隨方法，從而產(chǎn)生隨機(jī)的Galerkin 伴隨解。Sabater 等［41］通過(guò)耦合非嵌入式高斯過(guò)程模型和伴隨方程法搭建了魯棒優(yōu)化設(shè)計(jì)框架，為了適應(yīng)更多的不確定性輸入使用了梯度增強(qiáng)型的Kriging 模型，并成功應(yīng)用于考慮了馬赫數(shù)和迎角不確定性的跨聲速二維翼型的優(yōu)化設(shè)計(jì)，所提出的方法可以同時(shí)減少阻力的平均值和標(biāo)準(zhǔn)差。

結(jié)合以上總結(jié)和討論，目前基于多項(xiàng)式混沌法的不確定性分析和梯度優(yōu)化設(shè)計(jì)相關(guān)研究存在的不足為不確定性分析過(guò)程中計(jì)算成本較高，優(yōu)化設(shè)計(jì)并未包含多種設(shè)計(jì)狀態(tài)或輸入不確定性，且對(duì)優(yōu)化結(jié)果的分析不足。因此，本文的研究?jī)?nèi)容為首先利用伴隨方程法求解目標(biāo)函數(shù)對(duì)不確定性變量的導(dǎo)數(shù)，發(fā)展一種梯度增強(qiáng)型多項(xiàng)式混沌法，以降低優(yōu)化設(shè)計(jì)中不確定性分析模塊的計(jì)算成本，采用亞聲速和跨聲速二維翼型算例對(duì)該方法進(jìn)行驗(yàn)證，并量化不確定性變量的全局敏感性。其次，建立不確定性分析模塊統(tǒng)計(jì)矩梯度求解方法并搭建優(yōu)化設(shè)計(jì)系統(tǒng)。最后，對(duì)低亞聲速和跨聲速二維翼型進(jìn)行考慮馬赫數(shù)和迎角不確定性的梯度優(yōu)化設(shè)計(jì)，并對(duì)優(yōu)化結(jié)果進(jìn)行分析討論。

1 數(shù)值模擬方法和伴隨理論

1.1 數(shù)值模擬方法

采用有限體積法求解三維可壓縮雷諾平均Navier-Stokes（RANS）方程，其積分形式如式（1）所示。湍流模型使用Spalart-Allmaras（SA）方程模型。

式中：V為控制體的體積邊界；Q為流場(chǎng)守恒變量；FI和FV分別是無(wú)黏通量和有黏通量，在直角坐標(biāo)系下可分解為3 個(gè)方向的分量。

空間離散采用基于有限體積法的JST 空間離散方法。時(shí)間推進(jìn)上采用廣義最小殘差（Generalized Minimal Residual，GMRES）方法，該方法作為一種Newton-Krylov 子空間方法的分支，是目前較為先進(jìn)的線性系統(tǒng)迭代求解方法之一。

1.2 求解器驗(yàn)證

驗(yàn)證模型為ONERA-M6 機(jī)翼。計(jì)算狀態(tài)為Ma=0.84，Re=11.72×106，α=3.06°，該狀態(tài)具有較為完備的試驗(yàn)數(shù)據(jù)［42］。利用商業(yè)軟件ICEM 劃分多塊結(jié)構(gòu)O-H 型網(wǎng)格，網(wǎng)格量為158 萬(wàn)，邊界層第1 層高度為1×10-6。物面和對(duì)稱面網(wǎng)格如圖1 所示。

圖1 M6 機(jī)翼網(wǎng)格Fig.1 M6 wing mesh

沿機(jī)翼展向分別截取20%展長(zhǎng)、44%半展長(zhǎng)和65%半展長(zhǎng)3 個(gè)剖面，將求解器的計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，如圖2 所示。圖中，圓形、方形和菱形數(shù)據(jù)點(diǎn)為實(shí)驗(yàn)數(shù)據(jù)（Exp），實(shí)線為計(jì)算得到的截面壓力分布，可以看出計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)在各個(gè)機(jī)翼展向截面吻合均較好，說(shuō)明采用的CFD 求解器滿足優(yōu)化設(shè)計(jì)的精度要求。

圖2 M6 機(jī)翼計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig.2 Comparison of M6 wing calculation results with experimental data

1.3 伴隨方程

在梯度優(yōu)化設(shè)計(jì)中，高效精確地獲得目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)成為關(guān)鍵一環(huán)。伴隨方程法求解梯度的效率幾乎不受設(shè)計(jì)變量個(gè)數(shù)的影響，對(duì)于具有高維特性的氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題具有顯著優(yōu)勢(shì)。

給定目標(biāo)函數(shù)J和設(shè)計(jì)變量X。計(jì)算網(wǎng)格是設(shè)計(jì)變量的函數(shù)G(X)。對(duì)于RANS 方程，流場(chǎng)解是計(jì)算網(wǎng)格和設(shè)計(jì)變量的函數(shù)Q(G(X)，X)。

目標(biāo)函數(shù)通過(guò)對(duì)收斂的流場(chǎng)解Q進(jìn)行積分得到，有

對(duì)于收斂的流場(chǎng)解殘差R為0，即

當(dāng)給定流場(chǎng)解Q時(shí)，根據(jù)鏈?zhǔn)角髮?dǎo)法則得到目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的Jacobian矩陣：

為了避免求解計(jì)算最為耗時(shí)的dQ/dX一項(xiàng)，引入伴隨方程：

式中：ψ為伴隨向量。可表示為

從伴隨方程的組成上看，求解伴隨方程并不依賴于設(shè)計(jì)變量。盡管氣動(dòng)優(yōu)化問(wèn)題具有明顯的高維特性，然而求解伴隨方程需要的時(shí)間與單次CFD計(jì)算所需的時(shí)間大體相當(dāng)，大大提升了計(jì)算目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量梯度的效率。通過(guò)求解伴隨方程，可以得到伴隨向量，將式（6）代入式（4）中得到：

式中：?J/?G的計(jì)算使用收斂的流場(chǎng)解，對(duì)網(wǎng)格變量求一次偏導(dǎo)數(shù)即可；dG/dX可根據(jù)參數(shù)化方法和網(wǎng)格變形算法計(jì)算得到。

1.4 梯度驗(yàn)證

使用有限差分法對(duì)伴隨方程求解所得梯度進(jìn)行驗(yàn)證。為保證梯度驗(yàn)證的普適性，隨機(jī)挑選某翼型進(jìn)行梯度驗(yàn)證，如圖3 所示，在翼型四周布置FFD 控制框。FFD 控制框弦向共12 個(gè)控制點(diǎn)，控制點(diǎn)移動(dòng)方向?yàn)閥向。在翼型前后緣施加FFD 控制點(diǎn)的移動(dòng)約束，上下表面FFD 控制點(diǎn)在y方向的移動(dòng)大小相等，方向相反，以保證前后緣位置不變。由此總的設(shè)計(jì)變量個(gè)數(shù)為10×2+2（翼型前后緣設(shè)計(jì)變量）=22。

圖3 梯度驗(yàn)證翼型幾何和FFD 控制框Fig.3 Airfoil geometry and FFD control frame for gradient validation

該翼型計(jì)算狀態(tài)為Ma=0.69，α=-0.39°，Re=11.7×106，對(duì)應(yīng)于Honda 飛機(jī)巡航狀態(tài)［43］。表1 和表2 分別展示了針對(duì)目標(biāo)函數(shù)CL和CD利用有限差分方法和耦合伴隨方程計(jì)算得到的導(dǎo)數(shù)值。從對(duì)比結(jié)果上看，伴隨方程法求解得到的梯度與有限差分法吻合較好，相對(duì)誤差在10-3以下，符合有限差分截?cái)嗾`差精度。

表1 升力系數(shù)關(guān)于設(shè)計(jì)變量的梯度驗(yàn)證Table 1 Gradient verification of lift coefficients with respect to design variables

表2 阻力系數(shù)關(guān)于設(shè)計(jì)變量的梯度驗(yàn)證Table 2 Gradient verification of drag coefficients with respect to design variables

2 不確定性分析方法

2.1 多項(xiàng)式混沌法原理

考慮一個(gè)可由確定性映射表示的物理模型系統(tǒng)Q=M(X)，這里X=[X1，X2，…，XN]T，N≥1 是輸入變量向量（包括幾何設(shè)計(jì)變量和來(lái)流條件變量等）。Q=[Q1，Q2，…，QM]T，M≥1 是系統(tǒng)提供的一組輸出變量向量。

由于假設(shè)輸入向量X中部分n個(gè)變量受不確定性影響，因此這一部分的變量可由具有指定概率密度函數(shù)（Probability Density Function， PDF）的隨機(jī)向量x表示。另一部分的輸入變量為確定性變量，用D表示。這里x為真實(shí)隨機(jī)變量，其可看作標(biāo)準(zhǔn)隨機(jī)變量ξ的函數(shù)：

式中：μ和σ分別是該隨機(jī)變量的均值和方差。通過(guò)不確定性傳播，導(dǎo)致該系統(tǒng)的輸出響應(yīng)Q也是隨機(jī)變量。

多項(xiàng)式混沌法（Polynomial Chaos Expansion， PCE）的核心思想是通過(guò)含有獨(dú)立隨機(jī)變量的正交多項(xiàng)式基尋找隨機(jī)變量ξ（輸入?yún)?shù)）與感興趣量Q（QoI）之間的函數(shù)依賴關(guān)系。將隨機(jī)函數(shù)Q分解成確定和隨機(jī)兩部分，得到式（9）所示的無(wú)窮級(jí)數(shù)表達(dá)形式：

在實(shí)際使用中，取前P階有限模態(tài)截?cái)?，如式?0）所示。在氣動(dòng)優(yōu)化設(shè)計(jì)中D為確定性設(shè)計(jì)變量向量，ξ為隨機(jī)設(shè)計(jì)變量向量（如迎角、馬赫數(shù)、來(lái)流湍流度等）。那么αj(D)為該多項(xiàng)式混沌展開(kāi)第j階的確定性部分，而Ψj(ξ)則為該多項(xiàng)式混沌展開(kāi)第j階的隨機(jī)部分。

采樣總數(shù)Ns為P+1，由3 個(gè)參數(shù)確定：隨機(jī)變量的數(shù)量p，獨(dú)立基函數(shù)多項(xiàng)式的階數(shù)p和過(guò)采樣率np，如式（11）所示。

經(jīng)過(guò)文獻(xiàn)中的研究，過(guò)采樣率一般為2，即可在滿足多項(xiàng)式混沌法精度的情況下，最大限度地減少采樣點(diǎn)的個(gè)數(shù)。

多項(xiàng)式混沌法的關(guān)鍵在于求解正交多項(xiàng)式的權(quán)重系數(shù)αj(D)，采用點(diǎn)配置回歸法進(jìn)行求解。回歸方法即采用最小二乘方法求解通過(guò)確定性抽樣構(gòu)建的線性方程組。

在隨機(jī)空間中選擇Ns個(gè)樣本，并在這些點(diǎn)評(píng)估確定性問(wèn)題。最后，針對(duì)隨機(jī)變量的頻譜模式，建立了一個(gè)由Ns個(gè)方程組成的線性系統(tǒng)并求解。該線性系統(tǒng)由式（12）給出：

式中：矩陣中的ξ0代表該標(biāo)準(zhǔn)隨機(jī)變量矢量ξ在第一個(gè)采樣點(diǎn)處的值，同理ξNs-1為標(biāo)準(zhǔn)隨機(jī)變量矢量ξ在第Ns個(gè)采樣點(diǎn)處的值。最終待求系數(shù)矩陣寫成緊致格式為

多項(xiàng)式混沌法的最大優(yōu)點(diǎn)是可解析表達(dá)輸出隨機(jī)變量的均值和方差。根據(jù)對(duì)均值（數(shù)學(xué)期望）和方差的統(tǒng)計(jì)定義，由多項(xiàng)式混沌法得到的輸出隨機(jī)變量的均值和方差可表示為式（14）和式（15）：

2.2 梯度增強(qiáng)型多項(xiàng)式混沌法

多項(xiàng)式混沌法的精度和基函數(shù)的階數(shù)有關(guān)，為提高精度需要增加基函數(shù)的階數(shù)。由式（11）可知階數(shù)的增加會(huì)導(dǎo)致采樣點(diǎn)數(shù)的成倍增加，即需要調(diào)用成倍數(shù)量的CFD 數(shù)值模擬，需要大量的計(jì)算資源。

考慮到多項(xiàng)式混沌法精度和采樣點(diǎn)數(shù)的矛盾，采用伴隨方程法獲取梯度信息，進(jìn)而對(duì)多項(xiàng)式混沌法進(jìn)行增強(qiáng)。梯度增強(qiáng)型多項(xiàng)式混沌法（Gradient-enhanced Polynomial Chaos Expansion， GPCE）可以在保證不增加基函數(shù)階數(shù)的情況下顯著提高多項(xiàng)式混沌法的精度。簡(jiǎn)而言之，梯度增強(qiáng)型多項(xiàng)式混沌法可在保證相同的模擬精度下，大幅減少采樣點(diǎn)的個(gè)數(shù)。

使用這種梯度增強(qiáng)型多項(xiàng)式混沌法，可以構(gòu)建基于二階回歸的全局敏感性分析，其成本僅隨著維度的增加而線性增加。因此，這種方法甚至可以應(yīng)用于大維數(shù)問(wèn)題。

梯度增強(qiáng)型多項(xiàng)式混沌法可通過(guò)對(duì)式（10）所示的PCE 核心公式微分得到。等式的左右兩邊同時(shí)對(duì)歸一化的標(biāo)準(zhǔn)隨機(jī)變量求導(dǎo)數(shù)可以得到：

式（16）展開(kāi)寫成矩陣形式，如式（17）所示：

對(duì)于Ns個(gè)樣本ξ0，ξ1，…，ξNs-1，n個(gè)標(biāo)準(zhǔn)隨機(jī)變量ξi，ξi+1，…，ξn，n個(gè)真實(shí)隨機(jī)變量xi，xi+1，…，xn，左端項(xiàng)P+1 個(gè)多元聯(lián)合基函數(shù)Ψj分別對(duì)n個(gè)歸一化標(biāo)準(zhǔn)隨機(jī)變量ξi求導(dǎo)數(shù)?Ψj/?ξi。

式（17）右端項(xiàng)中?Q/?x為在每個(gè)確定性采樣點(diǎn)處感興趣的輸出變量Q對(duì)n個(gè)真實(shí)隨機(jī)輸入變量的確定性梯度，其可以由伴隨方程法直接得到。

由于真實(shí)隨機(jī)變量x又是標(biāo)準(zhǔn)隨機(jī)變量ξ的函數(shù)，那么?xj/?ξj可以基于不確定性變量的已知分布獲得。例如，一個(gè)均值為μ和標(biāo)準(zhǔn)差為σ的正態(tài)分布，單位化形式為

該導(dǎo)數(shù)為

如圖4 所示，通過(guò)在式（10）左右兩端增加梯度信息對(duì)線性系統(tǒng)進(jìn)行增廣處理。最終得到的基函數(shù)梯度增廣矩陣（17）大小為Ns(n+1)×Nb，右端項(xiàng)輸出矩陣大小為Ns(n+1)×Nq，PCE 待求系數(shù)矩陣的大小不變，仍然為Nb×Nq。采用梯度增強(qiáng)型PCE，需要的采樣點(diǎn)數(shù)仍為Ns，相比于同等精度的PCE 可減少nNs個(gè)采樣點(diǎn)。換言之，利用Ns個(gè)樣本點(diǎn)可構(gòu)建與Ns(n+1)個(gè)樣本點(diǎn)具有相同精度的多項(xiàng)式混沌響應(yīng)面。

圖4 梯度增強(qiáng)型多項(xiàng)式混沌法Fig.4 Gradient-enhanced polynomial chaos expansion method

2.3 基于方差分解的全局敏感性分析方法

通過(guò)不確定性分析，可以得到不確定性因素對(duì)輸出的影響，但是要確認(rèn)每一個(gè)不確定性變量對(duì)這一影響的貢獻(xiàn)程度，還需要進(jìn)行敏感性分析。敏感性分析分為全局敏感性分析和局部敏感性分析。全局敏感性分析可以克服局部敏感性分析的局限性，不僅考慮了各因素概率密度函數(shù)的分布和形態(tài)的影響，而且在分析時(shí)所有因素可同時(shí)變化，還可以不受模型限制并考察輸入?yún)?shù)在整個(gè)參數(shù)變化范圍內(nèi)對(duì)輸出的貢獻(xiàn)程度。本文采用基于方差分解的全局敏感性分析方法。

由方差降維分析，總方差可以表示為

式中：

方差分解顯示了如何將模型輸出的方差分解為可歸因于每個(gè)輸入的項(xiàng)，以及它們之間的交互作用。所有項(xiàng)加在一起就是模型輸出的總方差。根據(jù)以上的數(shù)學(xué)推導(dǎo)，定義一種基于直接方差分解的敏感性度量Sobol 指數(shù)，如式（22）所示。

式中：Si稱為“一階敏感性指標(biāo)”或“主效應(yīng)指標(biāo)”，其反映的是Xi對(duì)輸出方差的貢獻(xiàn)，由總方差進(jìn)行歸一化。同理，可以通過(guò)將方差分解中的其他項(xiàng)除以Var(Q) 來(lái)形成高階交互的Sobol 指數(shù)Sij，Sijk等：

基于多項(xiàng)式混沌法的Sobol 指數(shù)可以方便地表示為式（24）和式（25），詳細(xì)推導(dǎo)見(jiàn)文獻(xiàn)［44］。

2.4 統(tǒng)計(jì)矩梯度求解

基于梯度的不確定性優(yōu)化設(shè)計(jì)需要求解不確定性分析模塊的梯度信息，即包括隨機(jī)輸出變量均值對(duì)設(shè)計(jì)變量的梯度dμQ/dD和標(biāo)準(zhǔn)差對(duì)設(shè)計(jì)變量的梯度dσQ/dD。

由PCE 方法可知，隨機(jī)輸出變量Q的均值和方差可以解析地表達(dá)為式（14）和式（15）。將均值對(duì)設(shè)計(jì)變量進(jìn)行微分，得到式（26）。

式中：dQ/dD是每個(gè)確定性采樣點(diǎn)處輸出變量Q對(duì)設(shè)計(jì)變量的梯度，該項(xiàng)由伴隨方程法直接獲得。那么該式可以簡(jiǎn)單地理解為輸出變量均值對(duì)設(shè)計(jì)變量的梯度等于輸出變量對(duì)設(shè)計(jì)變量梯度的均值。

輸出變量方差對(duì)設(shè)計(jì)變量的梯度可以表示為

為了求式（27）的梯度值，需要確定多項(xiàng)式系數(shù)對(duì)設(shè)計(jì)變量的梯度dαj/dD。對(duì)于第個(gè)設(shè)計(jì)變量，可以通過(guò)對(duì)PCE 核心公式微分來(lái)構(gòu)造第二個(gè)多項(xiàng)式混沌展開(kāi)。式（10）的多項(xiàng)式混沌法核心公式左右兩端同時(shí)對(duì)設(shè)計(jì)變量求導(dǎo)可得

式中：左端項(xiàng)?Y(D，ξ)/?Dj可由伴隨方法求得；右端項(xiàng)中Ψi(ξ)為在構(gòu)建PCE 時(shí)選定的多項(xiàng)式基函數(shù)，那么待求項(xiàng)只剩這個(gè)新構(gòu)造的多項(xiàng)式展開(kāi)的新系數(shù)dαi(D)/dDj，同樣可利用回歸方法求解。

對(duì)于nDV個(gè)設(shè)計(jì)變量，需要構(gòu)建nDV個(gè)線性系統(tǒng)，如式（29）所示：

通過(guò)求解這nDV個(gè)線性系統(tǒng)，可得到輸出變量Q對(duì)nDV個(gè)設(shè)計(jì)變量的導(dǎo)數(shù)dQ/dD。

3 優(yōu)化設(shè)計(jì)框架

3.1 優(yōu)化設(shè)計(jì)輔助技術(shù)

采用FFD 參數(shù)化方法。FFD 參數(shù)化方法基于彈性體變形思想來(lái)表示幾何變形問(wèn)題。其核心公式為

式中：Xs為設(shè)計(jì)外形上任意一點(diǎn)的全局坐標(biāo)向量，(u，v，w)為其參數(shù)坐標(biāo)向量，u，v，w的取值范圍為(0，1)；Pi，j，k是控制框中每個(gè)控制點(diǎn)的全局坐標(biāo)向量。

選擇逆距離權(quán)重插值算法（Inverse Distance Weighting， IDW）作為本文的網(wǎng)格變形算法，以建立網(wǎng)格變形方法的實(shí)質(zhì)是建立空間網(wǎng)格與表面網(wǎng)格的映射關(guān)系。IDW 算法的核心公式為

式中：N是選為插值基點(diǎn)的個(gè)數(shù)；u(xv)為空間網(wǎng)格待插值點(diǎn)；uj為表面網(wǎng)格插值基點(diǎn)；wj為第j個(gè)位置的權(quán)重系數(shù)，與空間網(wǎng)格插值點(diǎn)到表面網(wǎng)格插值基點(diǎn)的距離成反比。

經(jīng)大量實(shí)踐發(fā)現(xiàn)，序列二次規(guī)劃算法是目前針對(duì)高維度大規(guī)模設(shè)計(jì)變量光滑非線性規(guī)劃問(wèn)題表現(xiàn)最優(yōu)的算法之一。本文使用基于序列二次規(guī)劃算法的SNOPT 作為優(yōu)化算法。

3.2 確定性優(yōu)化設(shè)計(jì)框架

結(jié)合3.1 節(jié)介紹的優(yōu)化設(shè)計(jì)輔助技術(shù)，可最終搭建如圖5 所示的確定性優(yōu)化設(shè)計(jì)系統(tǒng)。具體優(yōu)化步驟如下：

圖5 確定性優(yōu)化設(shè)計(jì)框架Fig.5 Deterministic optimization design frame

1）建立初始模型和網(wǎng)格，給定設(shè)計(jì)變量的初值X0。

2）通過(guò)FFD 方法對(duì)表面網(wǎng)格xs進(jìn)行參數(shù)化變形。

3）基于第2）步變形后的表面網(wǎng)格，通過(guò)IDW 動(dòng)網(wǎng)格技術(shù)得到變形后的空間網(wǎng)格xv。

4）求解器求解RANS 方程，計(jì)算流場(chǎng)結(jié)果Q。

5）基于流場(chǎng)結(jié)果，計(jì)算伴隨方程并求解目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度信息dJ/dX。

6）結(jié)合第4）步中流場(chǎng)計(jì)算結(jié)果和第5）步中的梯度信息，使用SNOPT 優(yōu)化算法對(duì)設(shè)計(jì)變量進(jìn)行更新。

7）重復(fù)第2）步～第6）步，直至優(yōu)化設(shè)計(jì)收斂。

3.3 不確定性優(yōu)化設(shè)計(jì)框架

考慮不確定性后，修改目標(biāo)函數(shù)為輸出變量的 均 值 和 方 差 的 線 性 組 合J=K1×μq+K2×σq，優(yōu)化設(shè)計(jì)系統(tǒng)中需引入多項(xiàng)式混沌法分析模塊來(lái)計(jì)算不確定性的傳播過(guò)程，并得到輸出變量的均值和方差及兩者對(duì)設(shè)計(jì)變量的梯度信息。最終可搭建如圖6 所示的不確定性優(yōu)化設(shè)計(jì)系統(tǒng)。具體優(yōu)化步驟如下：

圖6 不確定性優(yōu)化設(shè)計(jì)框架Fig.6 Uncertain optimization design frame

1）建立初始模型和網(wǎng)格，給定設(shè)計(jì)變量的初值X0。

2）通過(guò)FFD 方法對(duì)表面網(wǎng)格xs進(jìn)行參數(shù)化變形。

3）基于第2）步變形后的表面網(wǎng)格，通過(guò)IDW 動(dòng)網(wǎng)格技術(shù)得到變形后的空間網(wǎng)格xv。

4）設(shè)定不確定性變量，例如馬赫數(shù)M或迎角α等。根據(jù)不確定性變量在隨機(jī)空間進(jìn)行采樣，采樣數(shù)為Ns。

5）求解器在每個(gè)采樣點(diǎn)進(jìn)行確定性分析，求解RANS 方程，計(jì)算流場(chǎng)結(jié)果Q。

6）基于流場(chǎng)結(jié)果，在每個(gè)采樣點(diǎn)進(jìn)行確定性分析計(jì)算伴隨方程，得到不確定性輸出變量q對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)dq/dX。

7）基于確定性分析結(jié)果，構(gòu)建梯度增強(qiáng)型多項(xiàng)式混沌線性方程組ΨC=Q，求解系數(shù)矩陣C，得到輸出變量的均值μq和方差σq。

8）利用確定性分析所得梯度耦合多項(xiàng)式混沌系統(tǒng)的梯度，求解輸出變量均值及方差對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)dμq/dX，dσq/dX。那么目標(biāo)函數(shù)的梯度可由兩者進(jìn)行線性表示：

9）根據(jù)第8）步中的梯度信息，使用SNOPT優(yōu)化算法對(duì)設(shè)計(jì)變量進(jìn)行更新。

10）重復(fù)第2）步～第9）步，直至優(yōu)化設(shè)計(jì)收斂。

4 考慮飛行條件不確定性的翼型不確定性分析

本節(jié)將以亞聲速及跨聲速翼型算例進(jìn)行不確定性分析，同時(shí)可驗(yàn)證梯度增強(qiáng)型多項(xiàng)式混沌法的計(jì)算精度。

4.1 亞聲速翼型

選擇NACA0012 翼型進(jìn)行亞聲速狀態(tài)全湍流下多項(xiàng)式混沌法不確定性分析的精度測(cè)試。假設(shè)馬赫數(shù)滿足Ma～N(0.5，(0.02)2)的正態(tài)分布，迎角滿足α～N(1.5，(0.2)2)的正態(tài)分布。

計(jì)算高度為10 km，計(jì)算網(wǎng)格如圖7 所示，采用全湍流計(jì)算。為驗(yàn)證多項(xiàng)式混沌法的精度，分別進(jìn)行5 000、10 000 次蒙特卡洛法作為驗(yàn)證標(biāo)準(zhǔn)。蒙特卡洛法在采樣點(diǎn)總數(shù)足夠多時(shí)，其所得結(jié)果趨近于真實(shí)值。計(jì)算結(jié)果如表3 所示。2 次采樣的結(jié)果在均值上基本一致，方差上有較小的不同，數(shù)值上基本收斂，因此可作為驗(yàn)證多項(xiàng)式混沌法的標(biāo)準(zhǔn)。

表3 亞聲速翼型：蒙特卡洛法和PCE 法計(jì)算結(jié)果Table 3 Subsonic airfoil： Monte Carlo and PCE calculation results

圖7 NACA0012 網(wǎng)格Fig.7 NACA0012 mesh

為檢驗(yàn)PCE 方法的精度，定義PCE 方法的相對(duì)誤差為

式中：J=mean(CL)，mean(CD)，mean(Cmz)，std(CL)，std(CD)，std(Cmz)。

使用1～4 階的獨(dú)立基函數(shù)多項(xiàng)式，采樣點(diǎn)數(shù)從10～70 個(gè)不等，分別采用傳統(tǒng)PCE 方法和梯度增強(qiáng)型PCE 方法得到誤差在階數(shù)和采樣點(diǎn)數(shù)空間內(nèi)的分布，如圖8 所示。圖中顯示的是3 個(gè)力系數(shù)均值和方差的相對(duì)誤差總和，即

圖8 亞聲速翼型：傳統(tǒng)PCE 和梯度增強(qiáng)型PCE 的精度對(duì)比Fig.8 Subsonic airfoils： Accuracy comparison of conventional PCE and gradient-enhanced PCE

從表3 中可以看出和蒙特卡洛法的結(jié)果對(duì)比，GPCE 方法在預(yù)測(cè)力系數(shù)標(biāo)準(zhǔn)差的精度上具有明顯優(yōu)勢(shì)。同時(shí)，梯度信息的引入對(duì)低階多項(xiàng)式混沌法的精度增強(qiáng)效果較弱，但能在一定程度上減輕高階多項(xiàng)式混沌法的過(guò)擬合現(xiàn)象，從而提高精度。這一點(diǎn)從圖8 中也能看出，梯度增強(qiáng)型的PCE 方法大幅增加了采用3 階或4 階基函數(shù)時(shí)多項(xiàng)式混沌法的擬合精度。對(duì)于亞聲速全湍流狀態(tài)，采用1 階或2 階的多項(xiàng)式混沌法已經(jīng)可以滿足預(yù)測(cè)精度的要求，更高階的基函數(shù)不僅會(huì)增加采樣點(diǎn)的個(gè)數(shù)和計(jì)算成本，還會(huì)導(dǎo)致精度降低。

分別使用蒙特卡洛法和梯度增強(qiáng)型PCE 法計(jì)算得到的亞聲速翼型表面平均壓力分布系數(shù)沿弦向的分布如圖9 所示。圖中實(shí)線表示在不確定輸入的影響下壓力系數(shù)的均值，陰影部分表示壓力系數(shù)均值上下正負(fù)3 倍標(biāo)準(zhǔn)差區(qū)間。從圖中可以看出，受輸入不確定性影響最大的是翼型頭部，梯度增強(qiáng)型PCE 法模擬的均值和標(biāo)準(zhǔn)差與蒙特卡洛法幾乎完全重合，具有較高的精度。

圖9 亞聲速：表面壓力系數(shù)不確定性響應(yīng)對(duì)比Fig.9 Subsonic： Comparison of surface pressure coefficient uncertainty response

圖10 和圖11 分別展示了蒙特卡洛法和梯度增強(qiáng)型多項(xiàng)式混沌法計(jì)算得到的流場(chǎng)平均壓力系數(shù)分布和壓力系數(shù)在流場(chǎng)中的標(biāo)準(zhǔn)差云圖。從標(biāo)準(zhǔn)差云圖上看，NACA0012 翼型在給定的馬赫數(shù)和迎角輸入不確定性條件下，空間流場(chǎng)的不確定性主要集中在翼型的前50%弦長(zhǎng)，尤其是上表面的負(fù)壓峰處，受不確定性影響最為顯著。對(duì)比蒙特卡洛法得到的流場(chǎng)結(jié)果和多項(xiàng)式混沌法得到的流場(chǎng)結(jié)果，兩者高度相似。

圖10 亞聲速翼型：2 種方法流場(chǎng)壓力系數(shù)均值對(duì)比Fig.10 Subsonic airfoils： comparison of mean values of flow field pressure coefficients between two methods Subsonic airfoils

圖11 亞聲速翼型：2 種方法流場(chǎng)壓力系數(shù)標(biāo)準(zhǔn)差對(duì)比Fig.11 Subsonic airfoils： Comparison of standard deviation of flow field pressure coefficient between two methods

圖12 所示為翼型表面歸一化Sobol 指數(shù)沿弦向的分布。圖中［-1，0］為下表面，［0，1］為上表面。Sobol指數(shù)在上下表面近似呈對(duì)稱分布，在上下表面前90%弦長(zhǎng)范圍內(nèi)都由迎角的不確定性主導(dǎo)，而在前后緣的10%范圍內(nèi)受馬赫數(shù)不確定性影響較大。2個(gè)不確定性參數(shù)之間基本不存在耦合效應(yīng)。

圖12 亞聲速表面 Sobol 指數(shù)弦向變化Fig.12 Subsonic surface Sobol exponent chordwise variation

受馬赫數(shù)和迎角2 種不確定性影響的翼型空間流場(chǎng)歸一化Sobol 指數(shù)云圖如圖13 和圖14 所示。圖中顯示出在亞聲速全湍流狀態(tài)下，馬赫數(shù)不確定性主要影響翼型頭部和尾部的流場(chǎng)方差，而其他空間流域內(nèi)的不確定性主要由迎角的不確定性支配。兩種不確定性因素之間的耦合作用較小。

圖13 亞聲速Sobol 指數(shù)：迎角Fig.13 Subsonic sobol index： Angle of attack

圖14 亞聲速Sobol 指數(shù)：馬赫數(shù)Fig.14 Subsonic sobol index： Mach number

4.2 跨聲速翼型

選擇RAE2822 翼型進(jìn)行跨聲速狀態(tài)全湍流下多項(xiàng)式混沌法不確定性的精度測(cè)試，并對(duì)該狀態(tài)進(jìn)行不確定性分析和全局敏感性分析。

計(jì)算網(wǎng)格如圖15 所示，弦向布置281 個(gè)網(wǎng)格點(diǎn)，法向121 個(gè)網(wǎng)格點(diǎn)，附面層第1 層高度1×10-6。計(jì)算狀態(tài)為Re= 6.5×106，馬赫數(shù)滿足Ma～N(0.734，(0.02)2)的正態(tài)分布，迎角滿足α～N(2.79，(0.2)2)。選取的狀態(tài)流場(chǎng)中具有強(qiáng)激波，呈現(xiàn)高度非線性效應(yīng)，對(duì)驗(yàn)證工作有益［45］。分別進(jìn)行了5 000、10 000 次蒙特卡洛分析，得到力系數(shù)的計(jì)算結(jié)果如表4 所示。表中顯示5 000、10 000 次蒙特卡洛法的結(jié)果，2 次結(jié)果較為接近，可認(rèn)為蒙特卡洛法收斂，將其作為檢驗(yàn)多項(xiàng)式混沌法精度的標(biāo)準(zhǔn)。

表4 跨聲速翼型：蒙特卡洛法和PCE 法計(jì)算結(jié)果Table 4 Transonic airfoils： Monte Carlo and PCE calculation results

圖15 RAE2822 翼型計(jì)算網(wǎng)格Fig.15 RAE2822 airfoil mesh

分別使用1～4 階傳統(tǒng)PCE 方法和梯度增強(qiáng)型PCE 方法進(jìn)行該問(wèn)題的不確定性傳播，結(jié)果如表4 所示。在跨聲速條件下，流場(chǎng)非線性效應(yīng)大幅增加，從而增加了PCE 方法的預(yù)測(cè)難度。從表中可以看出，傳統(tǒng)PCE 方法和蒙特卡洛法在各個(gè)力系數(shù)均值的結(jié)果上差異較小，而在力系數(shù)標(biāo)準(zhǔn)差的預(yù)測(cè)上兩者誤差有所增大。梯度增強(qiáng)型PCE 方法的模擬精度較傳統(tǒng)PCE 方法具有一定提升，其主要體現(xiàn)在對(duì)力系數(shù)標(biāo)準(zhǔn)差的預(yù)測(cè)上，例如對(duì)比PCE（p=4）和GPCE（p=4）。

同樣，利用亞聲速下的精度驗(yàn)證方法，即式（33），選取采樣點(diǎn)數(shù)從6～50 個(gè)不等，分別采用傳統(tǒng)PCE 方法和GPCE 方法得到誤差在階數(shù)和采樣點(diǎn)數(shù)空間內(nèi)的分布，如圖16 所示。圖中顯示，在跨聲速狀態(tài)下GPCE 方法仍能在一定程度上減輕高階多項(xiàng)式混沌法的過(guò)擬合現(xiàn)象，例如4階和5 階。盡管GPCE 方法相比于傳統(tǒng)PCE 方法的優(yōu)勢(shì)在跨聲速狀態(tài)下整體上有所減弱，然而需要注意到其重點(diǎn)改善了低階、小樣本數(shù)情況下的擬合精度，這對(duì)于減少計(jì)算成本是有利的。

圖16 跨聲速翼型：傳統(tǒng)PCE 和梯度增強(qiáng)型PCE 的精度對(duì)比Fig.16 Transonic airfoils： Accuracy comparison of conventional PCE and gradient-enhanced PCE

使用蒙特卡洛法和GPCE 法計(jì)算得到的流場(chǎng)中壓力系數(shù)均值和標(biāo)準(zhǔn)差云圖如圖17 和圖18所示。對(duì)比2 種方法的結(jié)果可以看到，GPCE 法能夠較為精確地模擬不確定性在跨聲速流場(chǎng)中的傳播，包括流場(chǎng)中具有較強(qiáng)非線性的激波區(qū)域。從壓力分布均值上看，在該狀態(tài)下翼型上表面 60%弦長(zhǎng)處存在一道強(qiáng)激波。相應(yīng)的，標(biāo)準(zhǔn)差云圖中顯示出該位置受不確定性影響最明顯，標(biāo)準(zhǔn)差最大。

圖17 跨聲速翼型：2 種方法壓力系數(shù)均值對(duì)比Fig.17 Transonic airfoils： Comparison of mean values of flow field pressure coefficients between two methods

圖19 所示為翼型表面 Sobol 指數(shù)沿弦向的分布，在跨聲速狀態(tài)下，占主導(dǎo)地位的是馬赫數(shù)，在激波位置附近馬赫數(shù)、迎角及其耦合作用存在競(jìng)爭(zhēng)。

圖19 跨聲速表面 Sobol 指數(shù)弦向變化Fig.19 Transonic surface Sobol exponent chordwise variation

通過(guò)梯度增強(qiáng)型多項(xiàng)式混沌法求解的歸一化Sobol 指數(shù)云圖如圖20 所示。從圖中可以看出，在跨聲速全湍流流場(chǎng)中，馬赫數(shù)不確定性占主導(dǎo)地位，其影響流場(chǎng)中的大部分區(qū)域，尤其是激波位置處，但對(duì)下表面前緣區(qū)域影響較小，由于本例為正迎角工況，因此該處主要受迎角不確定性的影響。馬赫數(shù)和迎角的交互作用主要集中在翼型上表面激波前后區(qū)域。

圖20 跨聲速全局敏感性分析結(jié)果Fig.20 Transonic global sensitivity analysis results

5 考慮飛行條件不確定性的翼型確定性及不確定性優(yōu)化設(shè)計(jì)

5.1 低亞聲速翼型優(yōu)化設(shè)計(jì)

初始翼型選擇RAE2822，參考Cessna 狀態(tài)，Ma= 0.19，Re= 5×106，CL= 0.3。FFD 控制框如圖21 所示?？刂瓶蛟谝硇皖^部加密以更好地控制翼型前緣的型面。定義式（35）所示的確定性和不確定性優(yōu)化問(wèn)題。

圖21 翼型優(yōu)化FFD 控制框Fig.21 Airfoil optimization FFD control frame

優(yōu)化過(guò)程中設(shè)計(jì)變量總數(shù)為16 個(gè)，均為幾何設(shè)計(jì)變量。不確定性優(yōu)化中假設(shè)馬赫數(shù)和迎角分別滿足Ma～N(0.19，(0.05)2)的正態(tài)分布和α～N(0.45，(1.0)2)的正態(tài)分布，此時(shí)均值狀態(tài)下翼型升力系數(shù)為0.3。優(yōu)化過(guò)程中使用3 階梯度增強(qiáng)型多項(xiàng)式混沌法，預(yù)測(cè)誤差可保持在0.1%左右。

表5 給出了低亞聲速翼型優(yōu)化的確定性和不確定性力系數(shù)結(jié)果，其中各確定性的力系數(shù)為基準(zhǔn)流場(chǎng)結(jié)果，即Ma=μ(Ma)，α=μ(α)。確定性優(yōu)化設(shè)計(jì)結(jié)果命名為DeOpt （Deterministic Optimization）；不確定性優(yōu)化設(shè)計(jì)結(jié)果命名為UMOpt （Uncertainty Multiple factor Optimization）。從表5 中的確定性力系數(shù)可以看到，確定性優(yōu)化最終減阻10.03 counts，約為10.15%。其中主要減少的是壓差阻力，而摩擦阻力幾乎保持不變。不確定性優(yōu)化總減阻量為9.33 counts，略小于確定性優(yōu)化。另一方面，從阻力系數(shù)的不確定性統(tǒng)計(jì)響應(yīng)結(jié)果上看，不確定性優(yōu)化無(wú)論是均值還是標(biāo)準(zhǔn)差均要小于初始構(gòu)型和確定性優(yōu)化結(jié)果。相比于初始構(gòu)型，不確定性優(yōu)化結(jié)果阻力系數(shù)均值減少約17%，阻力系數(shù)標(biāo)準(zhǔn)差減少約80%。綜合以上兩個(gè)方面，可以得出結(jié)論：不確定性優(yōu)化在確定性基準(zhǔn)場(chǎng)的性能上有所犧牲，以換取在擾動(dòng)范圍內(nèi)性能的魯棒性。

表5 低亞聲速翼型優(yōu)化結(jié)果（阻力系數(shù)單位為counts）Table 5 Low subsonic airfoil optimization results （Drag coefficient in counts）

對(duì)比確定性和不確定性的結(jié)果，如圖22 和圖23 所示。UMOpt 翼型上表面前緣厚度較DeOpt 翼型有所增大，下表面前緣的前加載現(xiàn)象消失，厚度增大在壓力分布上體現(xiàn)為上下表面的順壓梯度相比于DeOpt 翼型增大，吸力峰位置前移。

圖22 低亞聲速確定性及不確定性優(yōu)化翼型對(duì)比Fig.22 Comparison of low subsonic deterministic and uncertain optimized airfoils

圖23 低亞聲速確定性及不確定性優(yōu)化壓力分布對(duì)比Fig.23 Comparison of low subsonic deterministic and uncertain optimal pressure distributions

采用蒙特卡洛法在隨機(jī)空間內(nèi)采樣并進(jìn)行確定性分析，得到的隨機(jī)樣本點(diǎn)阻力系數(shù)分布如圖24 所示。小提琴圖代表樣本點(diǎn)阻力系數(shù)的分布特性，圖中越飽滿的部分其概率密度越大，樣本點(diǎn)落在該區(qū)間的個(gè)數(shù)越多。每個(gè)小提琴圖中心的黑色粗實(shí)線代表上下四分位數(shù)的區(qū)間，中間點(diǎn)代表整體數(shù)據(jù)的中位數(shù)。因此直觀上看，小提琴圖的上下位置可近似代表阻力系數(shù)平均值的高低，越靠下的代表平均阻力系數(shù)更低；小提琴圖的形狀代表了阻力性能的魯棒性，越細(xì)長(zhǎng)的小提琴圖說(shuō)明阻力系數(shù)標(biāo)準(zhǔn)差更大魯棒性更差，越短粗的說(shuō)明阻力系數(shù)標(biāo)準(zhǔn)差更小，性能更加集中，魯棒性更好。對(duì)比3 種翼型的小提琴圖結(jié)果發(fā)現(xiàn)，UMOpt 翼型具有明顯優(yōu)勢(shì)。初始構(gòu)型無(wú)論是平均性能還是性能的魯棒性都較差，DeOpt翼型和UMOpt 翼型均大幅降低了阻力系數(shù)的平均值，但相比來(lái)講UMOpt 的性能更加向平均值集中。

圖24 低亞聲速確定性及不確定性優(yōu)化阻力系數(shù)的隨機(jī)分布Fig.24 Random distribution of drag coefficients for low subsonic deterministic and uncertain optimization

將初始構(gòu)型、DeOpt 翼型、UMOpt 翼型受馬赫數(shù)和迎角不確定性影響的空間流場(chǎng)壓力系數(shù)標(biāo)準(zhǔn)差云圖進(jìn)行對(duì)比。從圖25 中可以看出在上表面前緣位置處流場(chǎng)壓力系數(shù)標(biāo)準(zhǔn)差大小排序?yàn)閁MOpt＜DeOpt＜Initial。因此，全湍流下考慮不確定性的優(yōu)化設(shè)計(jì)相比于確定性優(yōu)化設(shè)計(jì)更能有效降低空間流場(chǎng)內(nèi)壓力系數(shù)的不確定性，體現(xiàn)了不確定性優(yōu)化設(shè)計(jì)的優(yōu)勢(shì)。

圖25 考慮馬赫數(shù)和迎角不確定性下初始、確定性及不確定性優(yōu)化結(jié)果空間流場(chǎng)壓力系數(shù)標(biāo)準(zhǔn)差云圖Fig.25 Standard deviation contour of pressure coefficient of initial， deterministic and uncertain optimization results considering uncertainty of Mach number and angle of attack

確定性和不確定性優(yōu)化歷程如圖26 所示。各個(gè)圖中確定性優(yōu)化歷程以三角形在線中標(biāo)示，不確定性優(yōu)化歷程以圓形在線中標(biāo)示，不確定性優(yōu)化過(guò)程中確定性基準(zhǔn)場(chǎng)阻力系數(shù)以灰色線表示。由于不確定性優(yōu)化中的目標(biāo)函數(shù)為J=μ(CD)+σ(CD)，因此灰色圓形標(biāo)注線和黑色圓形標(biāo)準(zhǔn)線之間的差可近似表示標(biāo)準(zhǔn)差的大小。圖中分別用不同顏色標(biāo)注了不確定性優(yōu)化過(guò)程中的關(guān)鍵節(jié)點(diǎn)位置，并在其旁邊展示了壓力分布的均值和標(biāo)準(zhǔn)差響應(yīng)結(jié)果以及翼型的變化。壓力分布圖中每一個(gè)關(guān)鍵迭代步除了可以和上一關(guān)鍵步對(duì)比壓力分布均值外，還可以始終與初始翼型對(duì)比壓力分布均值和標(biāo)準(zhǔn)差響應(yīng)范圍。

圖26 低亞聲速確定性及不確定性優(yōu)化歷程（黑色圓形線：DeOpt，確定性優(yōu)化目標(biāo)函數(shù)變化；黑色三角線：UMOpt，不確定性優(yōu)化目標(biāo)函數(shù)變化；灰色圓形線：UMOpt，不確定性優(yōu)化中基準(zhǔn)流場(chǎng)的阻力系數(shù)變化）Fig.26 Low subsonic deterministic and uncertain optimization history（Black circle line： DeOpt，change of objective function in deterministic optimization； black triangle line： UMOpt，change of objective function in uncertainty optimization； gray circle line： UMOpt， change of drag coefficient of mean flow field in uncertainty optimization）

不確定性優(yōu)化歷經(jīng)主迭代8 次，最終滿足約束和目標(biāo)函數(shù)收斂條件。優(yōu)化過(guò)程中前幾步主迭代目標(biāo)函數(shù)下降明顯，而最后幾步目標(biāo)函數(shù)變化較小，主要完成的是對(duì)壓力分布的調(diào)整和光順，優(yōu)化過(guò)程中阻力系數(shù)標(biāo)準(zhǔn)差大幅減小。不確定性優(yōu)化設(shè)計(jì)主要的優(yōu)化方向?yàn)樾》岣呱媳砻骖^部吸力峰峰值，并不斷推遲吸力峰出現(xiàn)的位置；降低下表面的順壓梯度。優(yōu)化的最終結(jié)果UMOpt 翼型基準(zhǔn)場(chǎng)阻力系數(shù)略大于DeOpt 翼型，但在翼型頭部附近的標(biāo)準(zhǔn)差響應(yīng)明顯小于DeOpt 翼型。

5.2 跨聲速翼型優(yōu)化設(shè)計(jì)

初始翼型選擇RAE2822，設(shè)計(jì)狀態(tài)參考Honda Jet，Ma= 0.70，Re= 7.93×106，CL=0.38。FFD 控制框與低亞聲速狀態(tài)下保持一致，即圖21 所示。對(duì)于該翼型分別進(jìn)行單點(diǎn)確定性優(yōu)化設(shè)計(jì)和考慮馬赫數(shù)和迎角不確定性的優(yōu)化設(shè)計(jì)。定義如式（36）所示的優(yōu)化問(wèn)題。相比于低亞聲速優(yōu)化增加力矩約束以控制低頭力矩。

在考慮馬赫數(shù)和迎角不確定性的優(yōu)化中假設(shè)馬赫數(shù)和迎角分別滿足Ma～N(0.70，(0.05)2)和α～N(0.387，(1.0)2)的正態(tài)分布，此時(shí)均值狀態(tài)下翼型升力系數(shù)為0.38。根據(jù)PCE 驗(yàn)證結(jié)果，在跨聲速優(yōu)化設(shè)計(jì)中選擇使用3 階梯度增強(qiáng)型多項(xiàng)式混沌法，可將誤差控制在1%以下。

表6 給出了跨聲速翼型優(yōu)化的確定性和不確定性力系數(shù)結(jié)果，其中各確定性的力系數(shù)為基準(zhǔn)流場(chǎng)結(jié)果，即Ma=μ(Ma)，α=μ(α)。UMOpt翼型確定性基準(zhǔn)場(chǎng)阻力略大于DeOpt 翼型，但從不確定性結(jié)果上看UMOpt 翼型的阻力系數(shù)均值（89.19 counts）相比于初始翼型減少約6%，略小于DeOpt 翼型的阻力系數(shù)均值（91.44 counts）。在阻力系數(shù)標(biāo)準(zhǔn)差響應(yīng)上，差異更加明顯，UMOpt 翼型僅為1.6 counts，與初始翼型相當(dāng)。而DeOpt 翼型的阻力系數(shù)標(biāo)準(zhǔn)差為5.11 counts，甚至高于初始翼型，表明僅考慮確定性性能的優(yōu)化結(jié)果有可能導(dǎo)致性能魯棒性變差。

表6 跨聲速翼型優(yōu)化結(jié)果（阻力系數(shù)單位為counts）Table 6 Transonic airfoil optimization results （Drag coefficient in counts）

圖27 和圖28 分別為確定性優(yōu)化和不確定性優(yōu)化的翼型結(jié)果對(duì)比和壓力分布對(duì)比。DeOpt翼型頭部半徑變大，上表面最大厚度位置前移，下表面前緣厚度變小。相應(yīng)的壓力分布上表面吸力峰增大。為了滿足力矩約束，后緣被卸載，原本的后加載被消除，升力載荷向頭部移動(dòng)。

圖27 跨聲速確定性及不確定性優(yōu)化翼型對(duì)比Fig.27 Comparison of transonic deterministic and uncertain optimal airfoils

圖28 跨聲速確定性及不確定性優(yōu)化壓力分布對(duì)比Fig.28 Comparison of transonic deterministic and uncertain optimal pressure distribution

從圖29 所示的小提琴圖中可以看出，UMOpt 翼型具有明顯優(yōu)勢(shì)。DeOpt 翼型在馬赫數(shù)和迎角的不確定性擾動(dòng)下具有比初始翼型更差的阻力性能魯棒性。相比之下，UMOpt 翼型在保持初始構(gòu)型魯棒性基本不變的情況下進(jìn)一步降低了阻力系數(shù)的平均值。

圖29 跨聲速確定性及不確定性優(yōu)化Fig.29 Random distribution of drag coefficients for transonic deterministic and uncertain optimization

跨聲速確定性和不確定性優(yōu)化歷程如圖30所示，不確定性優(yōu)化歷經(jīng)主迭代8 次，最終滿足約束和目標(biāo)函數(shù)收斂條件。優(yōu)化歷程中前半程目標(biāo)函數(shù)下降較快，后半程目標(biāo)函數(shù)變化較小，主要完成的是對(duì)壓力分布的調(diào)整和光順，優(yōu)化過(guò)程中阻力系數(shù)標(biāo)準(zhǔn)差變化較小，說(shuō)明全湍流考慮馬赫數(shù)和迎角不確定性的優(yōu)化僅降低了平均阻力性能，未明顯優(yōu)化魯棒性。不確定性優(yōu)化設(shè)計(jì)主要的優(yōu)化方向是不斷提高頭部吸力峰值，減小后部的升力載荷，將載荷向前緣移動(dòng)。

圖30 跨聲速確定性及不確定性優(yōu)化歷程（黑色圓形線：DeOpt，確定性優(yōu)化目標(biāo)函數(shù)變化；黑色三角線：UMOpt，不確定性優(yōu)化目標(biāo)函數(shù)變化；灰色圓形線：UMOpt，不確定性優(yōu)化中基準(zhǔn)流場(chǎng)的阻力系數(shù)變化）Fig.30 Transonic deterministic and uncertain optimization history（Black circle line： DeOpt， change of objective function in deterministic optimization；black triangle line： UMOpt， change of objective function in uncertainty optimization； gray circle line： UMOpt， change of drag coefficient of mean flow field in uncertainty optimization）

6 結(jié) 論

基于非嵌入式多項(xiàng)式混沌法并結(jié)合伴隨方程法建立了高效、可靠的基于梯度算法的不確定性優(yōu)化設(shè)計(jì)方法，相關(guān)研究總結(jié)如下：

1）利用伴隨方程法求得的導(dǎo)數(shù)信息發(fā)展了一種梯度增強(qiáng)型多項(xiàng)式混沌法，實(shí)現(xiàn)了高效精確的不確定性分析。通過(guò)亞聲速和跨聲速二維翼型的不確定性分析算例驗(yàn)證了所使用的梯度增強(qiáng)型多項(xiàng)式混沌法的精度。其可對(duì)力系數(shù)、空間流場(chǎng)壓力分布的不確定性響應(yīng)進(jìn)行較為全面而準(zhǔn)確的預(yù)測(cè)。

2）利用基于方差分解的全局敏感性分析方法對(duì)亞聲速和跨聲速二維翼型不確定性分析算例中不確定性變量的敏感性進(jìn)行量化。結(jié)果表明，前后緣的條狀區(qū)域常受到馬赫數(shù)不確定性的影響，駐點(diǎn)附近的流場(chǎng)壓力系數(shù)常受到迎角不確定性的影響。亞聲速狀態(tài)下對(duì)流場(chǎng)不確定性響應(yīng)產(chǎn)生主要影響的是迎角；對(duì)于存在激波的跨聲速流場(chǎng)，對(duì)流場(chǎng)壓力系數(shù)標(biāo)準(zhǔn)差貢獻(xiàn)最大的是馬赫數(shù)。

3）基于梯度增強(qiáng)型多項(xiàng)式混沌法、FFD 參數(shù)化方法、IDW 動(dòng)網(wǎng)格方法和SNOPT 優(yōu)化算法搭建了不確定性梯度優(yōu)化設(shè)計(jì)框架。利用該框架分別對(duì)低亞聲速二維翼型、跨聲速二維翼型進(jìn)行考慮馬赫數(shù)和迎角不確定性的優(yōu)化設(shè)計(jì)。

4）優(yōu)化結(jié)果表明，只考慮確定性性能的優(yōu)化設(shè)計(jì)有可能導(dǎo)致不確定性條件下性能的降低，而考慮不確定性的優(yōu)化設(shè)計(jì)相比于確定性優(yōu)化設(shè)計(jì)能夠顯著提高翼型抵抗馬赫數(shù)和迎角擾動(dòng)的能力，同時(shí)減小阻力均值和標(biāo)準(zhǔn)差，分別可達(dá)到約17%和80%。

研究結(jié)論表明，所建立的梯度增強(qiáng)型多項(xiàng)式混沌法以及與伴隨理論耦合的不確定性梯度優(yōu)化設(shè)計(jì)系統(tǒng)可以應(yīng)用于具有大規(guī)模確定性設(shè)計(jì)變量和少量不確定性變量的氣動(dòng)魯棒優(yōu)化設(shè)計(jì)。