劉艷，白俊強(qiáng)，華俊，2，劉南，王波

1.西北工業(yè)大學(xué)航空學(xué)院，陜西西安 710072；2.中國(guó)航空研究院，北京 100012 3.中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074

基于Stochastic Kriging的柔性機(jī)翼穩(wěn)健性優(yōu)化設(shè)計(jì)

劉艷1，白俊強(qiáng)1，華俊1，2，劉南1，王波3

1.西北工業(yè)大學(xué)航空學(xué)院，陜西西安 710072；2.中國(guó)航空研究院，北京 100012 3.中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074

采用隨機(jī)代理模型方法對(duì)柔性機(jī)翼氣動(dòng)外形進(jìn)行穩(wěn)健性優(yōu)化設(shè)計(jì)。相比確定性優(yōu)化設(shè)計(jì)，穩(wěn)健性設(shè)計(jì)能夠考慮設(shè)計(jì)變量和參數(shù)的擾動(dòng)，保持設(shè)計(jì)結(jié)果在不確定性影響下的性能穩(wěn)定。采用高精度的氣動(dòng)/結(jié)構(gòu)耦合求解器（耦合Navier-Stokes方程和結(jié)構(gòu)靜力學(xué)方程）分析柔性機(jī)翼的變形情況和氣動(dòng)效率。為了提高優(yōu)化效率，建立隨機(jī)Kriging（Stochastic Kriging，SK）代理模型，將確定性的Kriging代理模型發(fā)展到隨機(jī)空間，通過(guò)有限次輸入得到數(shù)據(jù)的固有不確定性。對(duì)柔性M6機(jī)翼的氣動(dòng)外形進(jìn)行穩(wěn)健性優(yōu)化設(shè)計(jì)，結(jié)果表明：相比確定性代理模型的穩(wěn)健性優(yōu)化結(jié)果，應(yīng)用隨機(jī)代理模型的優(yōu)化結(jié)果的設(shè)計(jì)點(diǎn)阻力系數(shù)減小2.8 counts，在可變馬赫數(shù)范圍內(nèi)阻力系數(shù)均值減小3.2 counts，優(yōu)化結(jié)果具有較高的設(shè)計(jì)點(diǎn)氣動(dòng)效率和阻力發(fā)散特性，并且優(yōu)化后構(gòu)型的翼根彎矩有明顯減小，體現(xiàn)隨機(jī)代理模型在穩(wěn)健性優(yōu)化設(shè)計(jì)系統(tǒng)中的優(yōu)勢(shì)，同時(shí)也說(shuō)明建立的SK代理模型具有較高的預(yù)測(cè)精度。

柔性機(jī)翼；穩(wěn)健性優(yōu)化；靜氣動(dòng)彈性力學(xué)；隨機(jī)Kriging代理模型；Navier-Stokes方程

目前工程中常用的優(yōu)化設(shè)計(jì)方法是確定性優(yōu)化設(shè)計(jì)，但該方法中設(shè)計(jì)變量、目標(biāo)函數(shù)和約束函數(shù)均為確定性問(wèn)題，因此不能夠反應(yīng)工程設(shè)計(jì)、產(chǎn)品制造以及服役環(huán)境對(duì)設(shè)計(jì)結(jié)果的影響，其中客觀存在大量不確定性因素會(huì)對(duì)設(shè)計(jì)方案產(chǎn)生影響，而傳統(tǒng)優(yōu)化設(shè)計(jì)的最優(yōu)解對(duì)于這些不確定性因素非常敏感［1-4］。與傳統(tǒng)的確定性優(yōu)化［5］相比，穩(wěn)健性優(yōu)化［6-7］的一個(gè)顯著特點(diǎn)就是在優(yōu)化過(guò)程中考慮不確定性因素的存在，減小設(shè)計(jì)目標(biāo)對(duì)不確定性因素的靈敏度，保持設(shè)計(jì)結(jié)果在不確定性影響下性能穩(wěn)定。穩(wěn)健性設(shè)計(jì)的求解方法可以分為確定性方法和隨機(jī)化方法兩類(lèi)。確定性方法［8］是將不確定性問(wèn)題通過(guò)一定的分析轉(zhuǎn)化為確定性問(wèn)題，然后利用相關(guān)優(yōu)化算法求解。而通過(guò)直接優(yōu)化擾動(dòng)量和約束的方法來(lái)處理不確定性問(wèn)題，因此也有時(shí)被稱為數(shù)值優(yōu)化方法。隨機(jī)化方法［9-10］通常都是基于直接搜索法開(kāi)展的研究，不需要顯式的一階或二階導(dǎo)數(shù)，而只需得到下一個(gè)位置處的目標(biāo)函數(shù)響應(yīng)，有效減少了對(duì)問(wèn)題的依賴性。

由于問(wèn)題規(guī)模爆炸式增加，考慮不確定性的穩(wěn)健性設(shè)計(jì)計(jì)算量很大。為了提高效率，可以將隨機(jī)化方法與各種代理模型相結(jié)合以處理穩(wěn)健性設(shè)計(jì)問(wèn)題。主要是根據(jù)一定數(shù)量的離散樣本數(shù)據(jù)，通過(guò)構(gòu)造響應(yīng)的代理模型進(jìn)行目標(biāo)特性分析，在保證一定分析精度的基礎(chǔ)上大幅縮減計(jì)算量，進(jìn)而加快整個(gè)設(shè)計(jì)系統(tǒng)的效率。這種方法由于相比傳統(tǒng)的蒙特卡羅模擬法（Monte Carlo Simulation，MCS）計(jì)算量小，可以得到可靠的低階矩信息等優(yōu)點(diǎn)，在不確定性領(lǐng)域得到了眾多關(guān)注。

綜上，本文針對(duì)柔性機(jī)翼在飛行狀態(tài)不確定性情況，以應(yīng)用較多的Kriging代理模型為原型建立隨機(jī)Kriging（Stochastic Kriging，SK）代理模型［11］，構(gòu)建基礎(chǔ)問(wèn)題，進(jìn)行統(tǒng)計(jì)學(xué)推導(dǎo)以及確立相關(guān)參數(shù)求解方案，得到可靠的估計(jì)均值和方差，并通過(guò)基礎(chǔ)算例驗(yàn)證這一方法在不確定性問(wèn)題中較強(qiáng)的實(shí)用性。通過(guò)與基于確定性Kriging（deterministic Kriging，DK）代理模型的穩(wěn)健性優(yōu)化結(jié)果進(jìn)行研究對(duì)比，驗(yàn)證隨機(jī)代理模型在穩(wěn)健性優(yōu)化設(shè)計(jì)中的優(yōu)勢(shì)。

1　理論基礎(chǔ)

1.1氣動(dòng)/結(jié)構(gòu)耦合理論基礎(chǔ)

由于要研究考慮飛行狀態(tài)不確定性影響的柔性機(jī)翼設(shè)計(jì)，需要機(jī)翼的靜彈性分析方法能夠捕捉不同飛行狀態(tài)下對(duì)流場(chǎng)造成的影響。因此建立耦合雷諾平均Navier-Stokes（RANS）方程和結(jié)構(gòu)靜力學(xué)方程的高精度靜氣動(dòng)彈性分析方法。氣動(dòng)模型與結(jié)構(gòu)模型之間載荷傳遞、彈性變形插值采用RBF插值方法進(jìn)行，參考文獻(xiàn)［11-12］中有詳細(xì)介紹，此處僅給出靜氣彈分析流程如圖1所示，圖中CFD求解器為求解RANS方程，CSD求解器采用柔度影響系數(shù)矩陣方法，收斂條件為兩次彈性變形差小于給定誤差值，動(dòng)網(wǎng)格采用四元數(shù)方法，載荷傳遞與位移插值均采用徑向基函數(shù)插值方法。

圖1　靜氣彈分析流程圖

1.2有限樣本的Stochastic Kriging建模

考慮不確定性的重復(fù)樣本點(diǎn)的情況下，將DK代理模型從確定性空間轉(zhuǎn)化到隨機(jī)空間中去［10］。為了解決包含各處非均等的不確定性和有限采樣點(diǎn)的問(wèn)題，對(duì)空間中的隨機(jī)樣本將通過(guò)如下形式表示：

式中，ε（x）被定義為固有不確定性，用來(lái)表示擾動(dòng)項(xiàng)，并被假設(shè)為正態(tài)分布，εj（x）是固有不確定性的一個(gè)實(shí)現(xiàn)，表示的是在x處重復(fù)樣本的擾動(dòng)，并且被假設(shè)為與其他重復(fù)樣本的固有和非固有不確定性M（x）之間相互獨(dú)立，η（x）代表的是固有不確定性的方差。這里的輸入變量為若干組（xi，Yj（xi），Ni；i =1，2，…，n），其中Ni表示在位置xi處的重復(fù)樣本的個(gè)數(shù)，這里的Ni是有限值，各個(gè)位置處樣本的均值可以相應(yīng)的表示為

此方法也可采用最優(yōu)線性估計(jì)的形式進(jìn)行預(yù)測(cè)，以得到未知位置的均值：

式中，wi（t）表示的是賦予各已知位置均值的權(quán)值量，也就是通過(guò)已知均值的權(quán)重來(lái)估計(jì)未采樣位置處的均值。對(duì)于不確定性的Stochastic Kriging模型［10］均方差的構(gòu)建形式可表示為

在對(duì)均方差的構(gòu)建中，很明顯，重復(fù)樣本的固有不確定性信息被考慮進(jìn)模型的構(gòu)建過(guò)程中，以表示假設(shè)的真值和估計(jì)的形式。通過(guò)對(duì)公式（4）進(jìn)行展開(kāi)：

顯然，公式（5）中平方和展開(kāi)的最后一項(xiàng)為零。根據(jù)上述獨(dú)立性假設(shè)，公式（5）即可寫(xiě)為

對(duì)公式（6）的推導(dǎo)中，比較有代表性的一項(xiàng)可以展開(kāi)為

式中，p、q表示的是輸入變量的某一位置，1≤p、q≤n。而后利用各重復(fù)樣本之間相互獨(dú)立的假設(shè)，公式（7）可以展開(kāi)為

為了簡(jiǎn)化公式中的符號(hào)，基于二階平穩(wěn)假設(shè)［14］的假設(shè)，這里先將（8）式簡(jiǎn)寫(xiě)為如下的形式

2　基于穩(wěn)健性設(shè)計(jì)的柔性機(jī)翼設(shè)計(jì)方法

建立的SK代理模型通過(guò)對(duì)基礎(chǔ)問(wèn)題的構(gòu)建、統(tǒng)計(jì)學(xué)推導(dǎo)以及相關(guān)參數(shù)求解方案的確立，得到了可靠的估計(jì)均值和方差，下面通過(guò)算例驗(yàn)證了這一方法在不確定性問(wèn)題中較強(qiáng)的適用性。柔性機(jī)翼穩(wěn)健性優(yōu)化設(shè)計(jì)流程如圖2所示。

圖2　穩(wěn)健性優(yōu)化設(shè)計(jì)流程圖

3　算例分析

本節(jié)用設(shè)計(jì)實(shí)例論述SK代理模型的可行性和魯棒性。以M6機(jī)翼為初始構(gòu)型，設(shè)計(jì)點(diǎn)狀態(tài)為：Ma =0.84，Cl=0.271 54，Re=1.172×107，考慮來(lái)流馬赫數(shù)不確定性變化范圍Ma∈［0.82，0.83，0.84，0.85，0.86］。3.1節(jié)采用基于DK代理模型進(jìn)行穩(wěn)健性優(yōu)化得到優(yōu)化解。在3.2節(jié)應(yīng)用SK代理模型進(jìn)行柔性機(jī)翼的穩(wěn)健性設(shè)計(jì)得到優(yōu)化解。3.3節(jié)對(duì)3.1、3.2節(jié)的最優(yōu)結(jié)果進(jìn)行對(duì)比，并對(duì)結(jié)果進(jìn)行分析。所有優(yōu)化設(shè)計(jì)在滿足機(jī)翼體積不減小的約束條件下進(jìn)行。

采用改進(jìn)的FFD參數(shù)化方法［14］，該方法克服了原始FFD方法的設(shè)計(jì)變量相關(guān)性強(qiáng)的缺點(diǎn)，適用于精細(xì)化外形參數(shù)化，F(xiàn)FD控制框及控制點(diǎn)如圖3所示，共有3個(gè)控制截面，每個(gè)剖面有12個(gè)控制點(diǎn)控制z坐標(biāo)，共36個(gè)控制點(diǎn)。再加外側(cè)2個(gè)剖面的扭轉(zhuǎn)角以及機(jī)翼后掠角，因此共有39個(gè)設(shè)計(jì)變量（xi，i=1，2，…，39）。設(shè)計(jì)變量范圍見(jiàn)表1，表中ID代表設(shè)計(jì)變量序號(hào)，lb代表設(shè)計(jì)變量變化的下限，ub代表設(shè)計(jì)變量變化的上限。

表1　39個(gè)設(shè)計(jì)變量范圍

圖3　FFD控制框及控制點(diǎn)

3.1基于確定性代理模型的穩(wěn)健性優(yōu)化

應(yīng)用DK代理模型對(duì)柔性M6機(jī)翼在滿足容積約束條件下，考慮馬赫數(shù)不確定變化范圍Ma∈［0.82，0.83，0.84，0.85，0.86］時(shí)，進(jìn)行穩(wěn)健性優(yōu)化設(shè)計(jì)，提高設(shè)計(jì)點(diǎn)的氣動(dòng)特性，同時(shí)改善優(yōu)化后機(jī)翼的阻力發(fā)散特性。

優(yōu)化模型可表述為：

在設(shè)計(jì)空間內(nèi)生成869個(gè)樣本點(diǎn)，選取19個(gè)為測(cè)試樣本點(diǎn)，其余850個(gè)構(gòu)建改進(jìn)的Kriging代理模型。應(yīng)用改進(jìn)的遺傳算法每代種群數(shù)為300個(gè)，進(jìn)化代數(shù)為30代，在保證機(jī)翼體積不減小的情況下進(jìn)機(jī)翼的穩(wěn)健性優(yōu)化設(shè)計(jì)。優(yōu)化結(jié)果在3.3節(jié)中進(jìn)行詳細(xì)討論。

3.2基于隨機(jī)代理模型的穩(wěn)健性優(yōu)化

基于構(gòu)建的SK代理模型對(duì)M6機(jī)翼進(jìn)行保證容積約束的穩(wěn)健性優(yōu)化設(shè)計(jì)，利用構(gòu)建的優(yōu)化系統(tǒng)在給定的約束條件下改善阻力發(fā)散特性。設(shè)計(jì)狀態(tài)為馬赫數(shù)范圍Ma［0.82，0.83，0.84，0.85，0.86］，優(yōu)化的目標(biāo)為降低阻力系數(shù)在擾動(dòng)范圍內(nèi)的均值和方差。優(yōu)化模型可表述為：

式中

為保證與使用確定性代理模型進(jìn)行優(yōu)化設(shè)計(jì)的采樣個(gè)數(shù)相同，將選用3.1節(jié)相同樣本點(diǎn)和優(yōu)化算法。與運(yùn)用確定性代理模型不同的是，在不同設(shè)計(jì)變量處進(jìn)行馬赫數(shù)變化區(qū)間內(nèi)5次均勻分布重復(fù)采樣建立隨機(jī)代理模型，通過(guò)在不同設(shè)計(jì)變量處的均值和方差作為SK模型的輸入變量，得到設(shè)計(jì)變量在未知位置處的響應(yīng)均值和方差。

選取869個(gè)樣本點(diǎn)中19個(gè)作為測(cè)試樣本點(diǎn)，用于驗(yàn)證SK代理模型是否合理。圖4給出SK代理模型預(yù)測(cè)的均值、方差與實(shí)際計(jì)算結(jié)果的對(duì)比曲線，結(jié)果吻合得很好，證明建立的SK代理模型的精度可以滿足之后的穩(wěn)健性優(yōu)化設(shè)計(jì)的需求。

圖4　SK代理模型預(yù)測(cè)的阻力均值與方差

3.3結(jié)果分析

圖5為優(yōu)化前后幾何模型對(duì)比圖。

圖5　優(yōu)化前后幾何模型對(duì)比

可以得出優(yōu)化后構(gòu)型的后掠角均增大，而SK模型優(yōu)化結(jié)果的后掠角最大。

圖6為阻力發(fā)散特性曲線對(duì)比可得，2種優(yōu)化結(jié)果在設(shè)計(jì)點(diǎn)處的阻力特性均得到改善，并且優(yōu)化后阻力發(fā)散特性均得到明顯提高。通過(guò)DK代理模型和SK代理模型的優(yōu)化結(jié)果對(duì)比曲線可得，SK代理模型的結(jié)果明顯優(yōu)于DK代理模型，SK模型對(duì)隨機(jī)問(wèn)題表現(xiàn)出了更強(qiáng)的預(yù)測(cè)能力。

圖6　阻力發(fā)散特性曲線對(duì)比

表2給出DK代理模型優(yōu)化解在設(shè)計(jì)點(diǎn)位置阻力降低21.56%。在5個(gè)馬赫數(shù)阻力系數(shù)均值降低21.33%，方差降至2.149×10-8。SK模型優(yōu)化解在設(shè)計(jì)點(diǎn)位置阻力降低23.30%，在可變馬赫數(shù)范圍為阻力系數(shù)均值降低23.23%，方差降至5.543 ×10-9。

表2　DK、SK代理模型穩(wěn)健性優(yōu)化結(jié)果對(duì)比

圖7為0.84馬赫數(shù)下運(yùn)用DK代理模型和SK代理模型進(jìn)行的穩(wěn)健性優(yōu)化設(shè)計(jì)結(jié)果在2個(gè)典型剖面的壓力分布對(duì)比，優(yōu)化構(gòu)型消除了機(jī)翼上表面激波，壓力分布恢復(fù)平緩。優(yōu)化后機(jī)翼前緣峰值有所降低，有利于阻力發(fā)散特性改善有利，應(yīng)用SK代理模型優(yōu)化結(jié)果的峰值比DK代理模型優(yōu)化的結(jié)果要低。

圖8為優(yōu)化前后機(jī)翼考慮柔性變形后無(wú)量綱翼根彎矩的對(duì)比曲線。由圖中可得：采用DK代理模型優(yōu)化解的翼根彎矩較原始構(gòu)型有所增加，而采用SK代理模型優(yōu)化解則明顯減小。可以得出，SK代理模型優(yōu)化后構(gòu)型沒(méi)有增加結(jié)構(gòu)的要求下，明顯降低了翼根剖面的彎矩。

圖7　典型剖面壓力系數(shù)曲線對(duì)比圖

圖8　無(wú)量綱翼根彎矩對(duì)比曲線

4　結(jié) 論

通過(guò)研究柔性機(jī)翼的基于DK與SK穩(wěn)健性優(yōu)化設(shè)計(jì)結(jié)果的對(duì)比，可看出基于SK代理模型的穩(wěn)健性優(yōu)化設(shè)計(jì)更能有效提高柔性機(jī)翼的阻力特性及阻力發(fā)散特性。應(yīng)用SK代理模型在優(yōu)化過(guò)程中考慮不確定因素的存在，可以有效解決穩(wěn)健性優(yōu)化設(shè)計(jì)問(wèn)題。SK代理模型優(yōu)化結(jié)果較DK代理模型優(yōu)化結(jié)果在設(shè)計(jì)點(diǎn)處阻力系數(shù)減小2.8 counts，在可變馬赫數(shù)范圍內(nèi)阻力系數(shù)均值減小3.2 counts。同時(shí)SK代理模型穩(wěn)健型優(yōu)化后機(jī)翼對(duì)翼根彎矩有明顯減小。顯示出所構(gòu)建SK代理模型在不確定性問(wèn)題的研究中較DK代理模型有優(yōu)勢(shì)。

［1］ Zang T，Hemsch M J，Hilburger M W，et al.Needs and Opportunities for Uncertainty-Based Multi Disciplinary Design Methods for Aerrospace Vehicles［R］.NASA TM-2002-211462

［2］ Zeeshan Q，Yunfeng D，Rafque A F，et al.Multidisciplinary Robust Design and Optimization of Multistage Boost Phase Interceptor［R］.AIAA-2010-2920

［3］ Yu X，Du X.Reliability-Based Multidisciplinary Optimization for Aircraft Wing Design［J］.Structure and Infrastructure Engineering 2006，2（3/4）：277-289

［4］ Sues R H，Oakley D R，Rhodes G S.MDO of Aeropropulsion Components Considering Uncertainty［R］.AIAA-1996-4062

［5］ Queipo N V，Haftka R T，Shyy W，et al.Tucker Surrogate-Based Analysis and Optimization［J］.Progress in Aerospace Sciences，2005，41：1-28

［6］ Schuelera G I，Jensen H A.Computational Methods in Optimization Considering Uncertainties—An Overview［J］.Computer Methods in Applied Mechanics and Engineering，2008，198（1）：2-13

［7］ Andradottir S.A Review of Simulation Optimization Techniques［C］∥Proceedings of the 1998 Winter Simulation Conference，IEEE，Piscataway，NJ，1998：151-158

［8］ Jin R，Du X，Chen W.The Use of Metamodeling Techniques for Optimization Under Uncertainty［J］.J Struct Multidisciplinary Optim，2003，25（2）：99-116

［9］ Hans Georg Beyer，Bernhard Sendhoff，Robust Optimization—A Comprehensive Survey［J］.Computer Methods in Applied Mechanics and Engineering，2007，196（33/34）：3190-3218

［10］Wang B，Bai Junqiang，Hae Chang Gea.Stochastic Kriging for Random Simulation Metamodeling with Finite Sampling［C］∥39th Design Automation Conference（DAC），Portland，Oregon，2013

［11］Liu Y，Bai J Q，Hua J.Static Aeroelasticity Analysis of High Aspect Ratio wing and the Jig-Shape Design Based on Multi-Block Structural Grid［J］.Advanced Engineering and Materials，2013，302：377-383

［12］劉艷，白俊強(qiáng)，華俊，等.非線性流固耦合分析方法研究［C］∥中國(guó)力學(xué)大會(huì)，2013

Liu Y，Bai J Q，Hua J，et al.An Approach to the Nonlinear Fluid-Structure Coupling Technical［C］∥CCTAM2013，2013（in Chinese）

［13］Journel A G，Huijbregts CJ.Mining Geostatistics［M］.London，Academic Proess，UK，1978

［14］Lamousin H J，Waggenspack W N.NURBS-Based Free-Form Deformations［J］.IEEE Computer Graphics and Applications，1994，11：59-65

Robust Optimization of Flexible Wing Using Stochastic Kriging Surrogate Model

Liu Yan1，Bai Junqiang1，Hua Jun1，2，Liu Nan1，Wang Bo3

1.College of Aeronautics，Northwestern Polytechnical University，Xi′an 710072，China 2.Chinese Aeronautics Establishment，Beijing 100012，China 3.China Academy of Aerospace Aerodynamics，Beijing 100074，China

The stochastic surrogate model method is applied in the robust optimization design of flexible wing.Comparing with deterministic optimization design，robust design can take into consideration the disturbances of design variables and parameters.Therefore，the performance of design results can be kept stable under uncertainties.A High-fidelity fluid-structure coupling solver（coupled Navier-Stokes equations and static structure equations）is used to analyze the deformation and aerodynamic characteristics of flexible wing.In order to enhance optimization efficiency，Stochastic Kriging（SK）surrogate model，which extends deterministic Kriging（DK）surrogate model into stochastic space，is built.The intrinsic uncertainties of data is acquired by finite number of inputs.The robust optimization of flexible M6 wing illustrates that comparied with robust optimization results of DK，the drag coefficient of optimal result of SK has reduced 2.8 counts and the mean values in variable MachNumberrange have reduced 3.2 counts.The optimal result of SK has higher aerodynamic efficiency in design point and better drag divergence characteristics.It is manifested that stochastic surrogate model is favorable in robust optimizationand the SK surrogate model possesses high prediction accuracy.

aeroelasticity，computational fluid dynamics，deformation，design，drag coefficient，elastic deformation，finite element method，flexible wings，flow fields，flow charting，forecasting，genetic algorithms，geometry，Mach number，mean square error，Navier Stokes equations，optimization，pressure distributions，Reynolds numbers，statistics，stochastic models，structural dynamics，robust optimization static aeroelasticity；stochastic Kriging（SK）surrogate model

V211.3 V214.3+2

A

1000-2758（2015）06-0906-07

2015-03-31

國(guó)家“973”計(jì)劃（2014CB744804）資助

劉艷（1988—），西北工業(yè)大學(xué)博士研究生，主要從事氣動(dòng)彈性分析與多學(xué)科優(yōu)化設(shè)計(jì)研究。