黃 勇 孟德虹* 羅振先 孫 巖 王運(yùn)濤 王曉軍 孫俊峰

（1.中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽(yáng) 621000；2.北京航空航天大學(xué)，北京 100083）

0 引言

從2010 年以來(lái)，國(guó)內(nèi)超級(jí)計(jì)算機(jī)的實(shí)測(cè)運(yùn)算速度已跨越P級(jí)(PetaFlops，即每秒千萬(wàn)億次浮點(diǎn)運(yùn)算)門(mén)檻，正向E級(jí)(ExaFlops，即每秒百億億次浮點(diǎn)運(yùn)算)進(jìn)發(fā)。“天河一號(hào)A”達(dá)到了2.507 P、“星云”達(dá)到了1.27 P?！疤旌佣?hào)”的峰值運(yùn)算速度為54.9 P，持續(xù)速度為33.9 P。2016 年，“神威·太湖之光” 的峰值運(yùn)算速度達(dá)到了125.4 P(約0.125 E)，持續(xù)速度為93.0 P。因此，可以說(shuō)我國(guó)已經(jīng)進(jìn)入基于超級(jí)計(jì)算機(jī)的高性能計(jì)算時(shí)代。

然而，就像美國(guó)懷俄明大學(xué)機(jī)械工程系教授Mavriplis等人[1]2007 年在美國(guó)看到的那樣：“計(jì)算方法正在成為除理論和實(shí)驗(yàn)之外科學(xué)探索的第三支柱，基于這樣的認(rèn)識(shí)，科學(xué)界越來(lái)越推崇高性能計(jì)算。但另一方面，計(jì)算工程卻持續(xù)地降格為一種在商品化計(jì)算機(jī)硬件上運(yùn)行的成熟軟件應(yīng)用，以為工程問(wèn)題不夠復(fù)雜，不需要部署如此大規(guī)模的先進(jìn)硬件。”Mavriplis等人[1]認(rèn)為千萬(wàn)億次計(jì)算機(jī)給航空基礎(chǔ)研究帶來(lái)了機(jī)遇，“工程應(yīng)用能從高性能計(jì)算方法的進(jìn)取計(jì)劃中得到同樣的好處，這些問(wèn)題至少同科學(xué)問(wèn)題一樣重要，尤其是那些關(guān)系到國(guó)家競(jìng)爭(zhēng)實(shí)力的問(wèn)題”。計(jì)算航空科學(xué)領(lǐng)域的“耦合多學(xué)科的數(shù)值模擬和航空航天飛行器全任務(wù)包線(xiàn)內(nèi)的設(shè)計(jì)優(yōu)化”就是這類(lèi)問(wèn)題，因?yàn)閱螌W(xué)科模擬的問(wèn)題相對(duì)更為成熟。

作為計(jì)算航空科學(xué)的重要組成部分，CFD的發(fā)展和應(yīng)用在過(guò)去五十年里極大地得益于計(jì)算機(jī)硬件的進(jìn)步，下一步從RANS 計(jì)算向LES 方向發(fā)展，也需要計(jì)算機(jī)硬件更大的進(jìn)步來(lái)支持。但是，美國(guó)斯坦福大學(xué)的Witherden 和Jameson[2]指出，目前高性能計(jì)算硬件的變化，給CFD數(shù)值算法設(shè)計(jì)帶來(lái)了新的約束。如處理器運(yùn)算能力和內(nèi)存帶寬的不平衡發(fā)展、大規(guī)模并行加速器(GPU和協(xié)處理器)帶來(lái)的異構(gòu)/異步問(wèn)題、節(jié)點(diǎn)通訊的延時(shí)問(wèn)題等等，使一種算法要適用于高效的大規(guī)模并行計(jì)算，也變得很不容易了。設(shè)計(jì)數(shù)值算法要從原來(lái)盡量減少浮點(diǎn)運(yùn)算數(shù)量，轉(zhuǎn)到考慮算法運(yùn)算強(qiáng)度上去；對(duì)傳輸?shù)拿孔止?jié)數(shù)據(jù)，都要考慮執(zhí)行了多少有用的計(jì)算。也就是說(shuō)，對(duì)一些早期開(kāi)發(fā)的大型工程應(yīng)用程序而言，要高效地使用目前的超級(jí)計(jì)算機(jī)，也是件困難的事。

國(guó)內(nèi)情況更是如此：超級(jí)計(jì)算機(jī)在航空航天領(lǐng)域的應(yīng)用還不夠廣泛深入，航空航天領(lǐng)域原先自研的大型應(yīng)用程序，對(duì)最新超級(jí)計(jì)算機(jī)的計(jì)算潛力挖掘得還不夠。為此，科技部國(guó)家重點(diǎn)研發(fā)計(jì)劃“高性能計(jì)算”重點(diǎn)專(zhuān)項(xiàng)，在2016 年啟動(dòng)了“數(shù)值飛行器原型系統(tǒng)開(kāi)發(fā)”項(xiàng)目，旨在提高航空航天領(lǐng)域自研應(yīng)用軟件適應(yīng)國(guó)內(nèi)超級(jí)計(jì)算機(jī)最新硬件的能力?？紤]大規(guī)模設(shè)計(jì)變量的氣動(dòng)-結(jié)構(gòu)綜合優(yōu)化方法開(kāi)發(fā)是其中的一個(gè)子項(xiàng)，本文討論這個(gè)領(lǐng)域的一個(gè)主要問(wèn)題：考慮靜氣彈變形影響的高速大展弦比機(jī)翼設(shè)計(jì)方法。

國(guó)內(nèi)有關(guān)文獻(xiàn)可追溯到1999 年詹浩等人[3]的改進(jìn)余量修正法工作，當(dāng)時(shí)通過(guò)全速勢(shì)方程和工程梁理論耦合計(jì)算機(jī)翼的彈性變形。2009 年詹浩等人[4]在余量修正反設(shè)計(jì)方法中，氣動(dòng)載荷及機(jī)翼結(jié)構(gòu)彈性變形通過(guò)基于非結(jié)構(gòu)網(wǎng)格的三維Euler方程耦合結(jié)構(gòu)力學(xué)平衡方程求解。2012 年王曉江等人[5]再次用余量修正法發(fā)展基于靜氣動(dòng)彈性的機(jī)翼型架外形設(shè)計(jì)方法，但采用三維NS 方程的有限體積法和結(jié)構(gòu)有限元方法弱耦合計(jì)算變形。

對(duì)線(xiàn)性小變形，迭代/修正方法的“性?xún)r(jià)比”其實(shí)挺高。2002 年梁強(qiáng)等人[6]耦合三維NS 方程與機(jī)翼彈性平衡方程，計(jì)算幾何扭轉(zhuǎn)角分布，對(duì)其迭代更新進(jìn)行機(jī)翼型架外形設(shè)計(jì)。2011 年楊國(guó)偉等人[7]基于混合網(wǎng)格的NS 方程和結(jié)構(gòu)柔度矩陣方法的耦合求解變形，不進(jìn)行迭代計(jì)算，直接采用結(jié)構(gòu)負(fù)載的卸載方法，從巡航飛行外形反推地面型架外形所需要的結(jié)構(gòu)變形量。2012 年林躍勝等人[8]對(duì)反向加載得到的初始型架外形，采用修正的偶極子格網(wǎng)法計(jì)算氣動(dòng)力，發(fā)展機(jī)翼型架外形的快速設(shè)計(jì)方法。2014 年黃江濤等人[9]也采用結(jié)構(gòu)負(fù)載的卸載方法，但在型架外形逆推、修正及靜氣動(dòng)彈性復(fù)核中，驗(yàn)證了其CFD網(wǎng)格變形和CFD/CSD數(shù)據(jù)交換技術(shù)的效率和精度。

在基于梯度優(yōu)化設(shè)計(jì)方法中，國(guó)內(nèi)目前有三類(lèi)工作。只用氣動(dòng)梯度的，如2009 年熊俊濤等人[10]在機(jī)翼的控制論優(yōu)化設(shè)計(jì)方法中，通過(guò)Euler方程和結(jié)構(gòu)靜平衡方程耦合，求解機(jī)翼在跨聲速巡航飛行狀態(tài)下達(dá)到結(jié)構(gòu)彈性平衡時(shí)的彈性變形及真實(shí)氣動(dòng)載荷。氣動(dòng)梯度和結(jié)構(gòu)梯度都用的，如2018 年徐兆可等人[11]用基于梯度的SQP算法進(jìn)行機(jī)翼氣動(dòng)結(jié)構(gòu)優(yōu)化設(shè)計(jì)，其中氣動(dòng)力通過(guò)Euler方程與結(jié)構(gòu)平衡方程的靜氣彈耦合計(jì)算得到；在此基礎(chǔ)上，升阻力系數(shù)的梯度通過(guò)求解流場(chǎng)伴隨方程獲得，翼盒結(jié)構(gòu)質(zhì)量的梯度直接計(jì)算求得，應(yīng)力梯度通過(guò)擾動(dòng)結(jié)構(gòu)設(shè)計(jì)變量采用有限差分方法獲得。還有計(jì)劃用氣動(dòng)和結(jié)構(gòu)交叉梯度的，如2018 年黃江濤等人[12]的工作，為高效地計(jì)算優(yōu)化目標(biāo)相對(duì)氣動(dòng)和結(jié)構(gòu)設(shè)計(jì)變量的梯度，針對(duì)RANS 和結(jié)構(gòu)平衡方程的耦合系統(tǒng)，研究了其對(duì)應(yīng)耦合伴隨方程的延遲數(shù)值計(jì)算方法。

基于代理模型的優(yōu)化設(shè)計(jì)方法，考慮計(jì)算資源限制，僅在樣本分析中采用高保真耦合模型。2006年劉金輝等人[13]基于響應(yīng)面的多目標(biāo)多約束優(yōu)化設(shè)計(jì)方法中，通過(guò)三維Euler方程的有限體積法和機(jī)翼彈性平衡方程的有限元方法的耦合迭代計(jì)算靜氣彈變形。2014 年左英桃等人[14]基于Kriging模型利用遺傳算法優(yōu)化設(shè)計(jì)，其中結(jié)構(gòu)分析用有限元方法計(jì)算，氣動(dòng)載荷用RANS 方程計(jì)算，流固耦合界面數(shù)據(jù)插值采用RBF插值方法。2017 年聶雪媛等人[15]采用Kriging模型結(jié)合改進(jìn)的粒子群優(yōu)化方法，基于楊國(guó)偉的CFD/CSD耦合分析，對(duì)大飛機(jī)的結(jié)構(gòu)剛度進(jìn)行了優(yōu)化設(shè)計(jì)研究。

不用代理模型的智能優(yōu)化設(shè)計(jì)方法，考慮計(jì)算資源限制，需要快速的耦合計(jì)算，采用的模型保真度低些。2013 年楊佑緒等人[16]基于遺傳算法提出一種氣動(dòng)彈性綜合優(yōu)化設(shè)計(jì)方法，采用等效板多板模型計(jì)算翼面結(jié)構(gòu)動(dòng)力學(xué)特性，利用面元法計(jì)算氣動(dòng)力，在考慮顫振速度、靜氣動(dòng)彈性變形約束的情況下，以結(jié)構(gòu)質(zhì)量最小為目標(biāo)開(kāi)展了優(yōu)化設(shè)計(jì)方法研究。2017 年楊體浩等人[17]使用改進(jìn)的微分進(jìn)化算法優(yōu)化，在靜氣彈耦合計(jì)算上氣動(dòng)力用快速的全速勢(shì)方程加附面層修正計(jì)算，而結(jié)構(gòu)減重優(yōu)化用商業(yè)軟件實(shí)現(xiàn)。

回顧以上設(shè)計(jì)方法研究(不包括MDO方法和全商業(yè)軟件集成的工作)可以看出，過(guò)去二十年里，靜氣彈耦合計(jì)算用的氣動(dòng)和結(jié)構(gòu)模型，在不可逆轉(zhuǎn)地向計(jì)算資源消耗型的高保真模型發(fā)展；以此為基礎(chǔ)的設(shè)計(jì)方法，如迭代/修正法、余量修正法、伴隨方程優(yōu)化法、耦合伴隨方程優(yōu)化法、代理模型優(yōu)化法、智能優(yōu)化法等，都盡量避免在計(jì)算資源消耗上進(jìn)一步雪上加霜。而事實(shí)上，如前所述，目前國(guó)內(nèi)已經(jīng)具有很好的高性能計(jì)算資源，所以，有必要根據(jù)型號(hào)的實(shí)際需求，利用好這些資源，發(fā)展更有力的優(yōu)化設(shè)計(jì)工具。

2018 年薛幫猛等人[18]在這方面邁出了一步。在“天河二號(hào)”上探索民機(jī)飛發(fā)集成復(fù)雜構(gòu)型的精細(xì)優(yōu)化設(shè)計(jì)，用遺傳算法直接調(diào)用RANS 分析開(kāi)展全局尋優(yōu)的多目標(biāo)優(yōu)化，以期得到可靠的減阻效果。本文則進(jìn)一步考慮機(jī)翼靜氣彈變形的影響，基于高保真氣彈分析，以在超算上建立大規(guī)模設(shè)計(jì)變量(十萬(wàn)量級(jí))的氣動(dòng)—結(jié)構(gòu)綜合優(yōu)化方法為目標(biāo)，開(kāi)展相關(guān)工作。

1 優(yōu)化方法

1.1 機(jī)翼外形參數(shù)化方法

基于描述幾何外形的某種數(shù)學(xué)模型，采用有限個(gè)參數(shù)來(lái)控制機(jī)翼外形變化的方法，稱(chēng)為機(jī)翼外形參數(shù)化方法。在優(yōu)化設(shè)計(jì)中，可認(rèn)為是對(duì)設(shè)定連續(xù)設(shè)計(jì)空間的一種離散近似，它對(duì)尋優(yōu)的質(zhì)量和效率有著重要影響。作為它的基礎(chǔ)，Guenov等人[19]對(duì)翼型參數(shù)化提出了完整性、正交性、無(wú)暇性、簡(jiǎn)約性和直覺(jué)性的評(píng)價(jià)要求。完整性指是否能按指定容差描述任何翼型；正交性指確保每個(gè)翼型與一套參數(shù)一一對(duì)應(yīng)；無(wú)暇性要求參數(shù)化方法不會(huì)生成不良外形，諸如交叉的或退化的翼型；簡(jiǎn)約性指翼型主要幾何特征產(chǎn)生顯著變化只需少量參數(shù)；直覺(jué)性要求與參數(shù)物理意義相關(guān)，可簡(jiǎn)化輸入界限選擇或設(shè)計(jì)判斷。

盡管目前已有很多參數(shù)化方法，但似乎還沒(méi)有一種按這五條標(biāo)準(zhǔn)來(lái)看都是完美的。僅就完整性而言，Rendall等人[20]用2 174 個(gè)翼型測(cè)試了七種參數(shù)化方法；其中構(gòu)造類(lèi)方法純粹根據(jù)指定參數(shù)描述外形，測(cè)試了CST、B樣條、PARSEC三種；變形類(lèi)方法基于已有外形變形生成新外形，測(cè)試了Hicks-Henne鼓包函數(shù)方法、RBF域元法、Bezier面的FFD方法和SVD模態(tài)提取法四種。該研究發(fā)現(xiàn)，這七種方法按Kulfan 的容差要求描述所有翼型需要20 個(gè)～25個(gè)參數(shù)，用13 個(gè)～18 個(gè)參數(shù)僅能覆蓋其中80%的翼型。也就是說(shuō)，這些參數(shù)化方法的完整性與其簡(jiǎn)約性就不能兩全其美。由此推測(cè)，增加參數(shù)數(shù)量有利于增大參數(shù)化方法對(duì)設(shè)計(jì)空間的覆蓋能力，可能促使優(yōu)化要向大規(guī)模變量方向發(fā)展才更有潛力。

國(guó)內(nèi)目前翼型優(yōu)化大多用二維CST方法參數(shù)化，而機(jī)翼優(yōu)化大多用FFD方法。考慮到不僅要用機(jī)翼平面參數(shù)和盡量多的型面參數(shù)，而且要有較好的無(wú)暇性、直覺(jué)性。本文應(yīng)用三維類(lèi)別函數(shù)/形狀函數(shù)變換方法[21]CST(Class-function/Shape-function Transformation)定義機(jī)翼。

設(shè)(x，y，z)為機(jī)翼表面流向、展向和縱向的右手系笛卡爾坐標(biāo)，用下標(biāo)U，N，T分別表示上翼面、前緣點(diǎn)和后緣點(diǎn)。再設(shè)b為機(jī)翼展長(zhǎng)，用c、xN、αTW表示某展向位置的弦長(zhǎng)、前緣點(diǎn)流向坐標(biāo)和幾何扭轉(zhuǎn)角，則翼面無(wú)量綱的流向和展向坐標(biāo)(Ψ，η)可定義為：

這時(shí)機(jī)翼上表面的縱向坐標(biāo)可表示為(下翼面類(lèi)似)：

其中，機(jī)翼當(dāng)?shù)責(zé)o量綱切向位移(計(jì)當(dāng)?shù)厍熬夵c(diǎn)，描述機(jī)翼上反)：

機(jī)翼當(dāng)?shù)責(zé)o量綱后緣半厚度：

類(lèi)別函數(shù)定義為：

形狀函數(shù)定義為：

這里系數(shù)Aij(i=1，…，n，j=1，…，m)是待求的參數(shù)，以后用作控制翼面形狀的設(shè)計(jì)參數(shù)。Kulfan的CST方法里Si(Ψ)和Sj(η)是Bernstein 多項(xiàng)式，可寫(xiě)成(二者類(lèi)似寫(xiě)一個(gè)，下同)：

參數(shù)化描述氣動(dòng)外形，考慮到優(yōu)化中的計(jì)算開(kāi)銷(xiāo)，一般希望參數(shù)Aij的個(gè)數(shù)盡量少。但有時(shí)擔(dān)心如此約束了設(shè)計(jì)空間，妨礙尋優(yōu)，又希望參數(shù)個(gè)數(shù)多一點(diǎn)，即n，m的取值大一些。這時(shí)采用Bernstein 多項(xiàng)式的次數(shù)高，描述的外形會(huì)存在比較大的波動(dòng)。因?yàn)樾螤詈瘮?shù)的定義就是張量積Bezier曲面，所以本文推廣一下，采用遞推定義的B樣條基函數(shù)[22]Ni，p(Ψ)來(lái)描述Si(Ψ)和Sj(η)，把形狀函數(shù)定義成B樣條曲面，這樣能大幅改善參數(shù)數(shù)量比較多時(shí)帶來(lái)的波動(dòng)問(wèn)題。即：

其中，p是流向Ψ方向曲線(xiàn)可指定的次數(shù)，另設(shè)q為展向η方向曲線(xiàn)可指定的次數(shù)。Ψi是根據(jù)原剖面點(diǎn)定義的節(jié)點(diǎn)，定義同前，但首末要取重節(jié)點(diǎn)以保持原CST方法控制前緣半徑和后緣船尾角的能力。Aij可用作設(shè)計(jì)參數(shù)，在原CST方法里來(lái)源于Bernstein 基函數(shù)的最大值特性，但它的影響仍是全局的。改用具有局部支承性的B樣條基函數(shù)，Aij的影響本質(zhì)上是局部的。與本文的替換方式不同，van Tooren 等人[23]用B樣條基函數(shù)與Bernstein 基函數(shù)相乘來(lái)拓展CST。

用以上方法進(jìn)行機(jī)翼參數(shù)化描述，需四個(gè)平面設(shè)計(jì)參數(shù)，包括：機(jī)翼面積、展弦比、根梢比、1/4弦線(xiàn)后掠角。上反、扭轉(zhuǎn)分布至少要在翼根和翼梢指定兩個(gè)參數(shù)以實(shí)現(xiàn)線(xiàn)性分布，對(duì)一般的非線(xiàn)性分布可用簡(jiǎn)單的曲線(xiàn)參數(shù)化來(lái)描述，數(shù)量可定制。在指定較低的流向和展向的曲線(xiàn)次數(shù)(即p，q)后，給定若干展向站位的翼剖面數(shù)據(jù)，可擬合出型面設(shè)計(jì)參數(shù)Aij，擬合精度隨n，m增大迅速改善。

圖1 顯示，針對(duì)CRM機(jī)翼采用以上方法在p=3，q=1 時(shí)，當(dāng)m＞18，n ＞11 后表征擬合精度的最大殘值Rmax才不大幅下降，以英寸為單位的全尺寸外形坐標(biāo)可精確到千分位。圖2 給出某個(gè)Aij取值變化對(duì)機(jī)翼上表面外形的影響，具有B樣條曲面很規(guī)范的局部特性。圖3 是在固定機(jī)翼面積下變化展弦比和1/4 弦線(xiàn)后掠角的效果，結(jié)合圖2，展示了三維CST參數(shù)化方法既能修改機(jī)翼平面參數(shù)也能修改型面參數(shù)的能力。

圖1 型面參數(shù)數(shù)量對(duì)擬合精度的影響

圖2 型面參數(shù)變化的影響

圖3 平面參數(shù)變化的影響(固定面積)

1.2 靜氣動(dòng)彈性的CFD/CSD耦合分析方法

靜氣動(dòng)彈性的CFD/CSD耦合方法可分成強(qiáng)耦合(聯(lián)立求解CFD和CSD方程)與弱耦合(分開(kāi)求解CFD和CSD方程，在界面交換信息)兩類(lèi)，國(guó)內(nèi)似乎沒(méi)人發(fā)展強(qiáng)耦合方法，盡管比較困難，但對(duì)柔性大變形問(wèn)題這種方法可能更有利。弱耦合又可分成緊耦合(CFD和CSD方程迭代求解過(guò)程中交換信息)與松耦合(CFD和CSD方程迭代求解收斂后交換信息)兩種，國(guó)內(nèi)以緊耦合居多，不僅耦合時(shí)間精度高還可用于動(dòng)氣彈計(jì)算。

2006 年徐敏等人[24]建立的那種CFD/CSD耦合分析方法，很像是界面信息傳輸采用守恒體積轉(zhuǎn)換(CVT)方法的松耦合。2009 年楊國(guó)偉[25]介紹的進(jìn)展是結(jié)構(gòu)運(yùn)動(dòng)方程和氣動(dòng)主宰方程的緊耦合，同時(shí)還介紹了基于虛功原理的兩種氣動(dòng)—結(jié)構(gòu)數(shù)據(jù)插值技術(shù)IPS 和RBF，以及三種動(dòng)網(wǎng)格方法：適合多塊結(jié)構(gòu)網(wǎng)格中小變形的TFI方法、需要迭代求解的非結(jié)構(gòu)網(wǎng)格彈簧網(wǎng)格平衡法、適用大規(guī)模網(wǎng)格基于背景網(wǎng)格的高效變形方法。2013 年張偉偉等人[26]從變形能力、變形質(zhì)量和變形效率角度，重點(diǎn)總結(jié)了物理模型和數(shù)學(xué)插值兩類(lèi)的六種網(wǎng)格變形方法，并介紹了其在使用RBF進(jìn)行網(wǎng)格變形和流固耦合界面數(shù)據(jù)傳遞上的收益，建議關(guān)注混合方法、網(wǎng)格變形算法的并行化等。2018 年楊超等人[27]對(duì)氣動(dòng)彈性分析中使用的氣動(dòng)力方法進(jìn)行了綜述，指出“滿(mǎn)足工程氣動(dòng)彈性設(shè)計(jì)的需求，特別是對(duì)于大柔性飛行器的全機(jī)氣動(dòng)彈性問(wèn)題”，CFD/CSD耦合要解決與飛行力學(xué)、結(jié)構(gòu)動(dòng)力學(xué)及飛行控制等的多學(xué)科耦合計(jì)算方法和效率問(wèn)題。這對(duì)應(yīng)用超算提出了需求。

本文采用的靜氣動(dòng)彈性CFD/CSD耦合分析方法中，多塊結(jié)構(gòu)網(wǎng)格變形采用孫巖等人[28-29]發(fā)展的RBF和TFI混合方法，耦合數(shù)據(jù)傳遞采用TPS 方法，自主開(kāi)發(fā)的CFD軟件TRIP與CSD的緊耦合測(cè)試驗(yàn)證，詳見(jiàn)孫巖等人2017 年的介紹[30]。對(duì)所用動(dòng)網(wǎng)格技術(shù)，王運(yùn)濤等人[31]在“天河二號(hào)”集群上開(kāi)展了十億量級(jí)網(wǎng)格變形的并行效率測(cè)試，對(duì)1.8 億網(wǎng)格的CRM-WBH構(gòu)型及其風(fēng)洞試驗(yàn)?zāi)Ｐ偷挠邢拊Ｐ停_(kāi)展了靜氣動(dòng)彈性模塊并行效率及一致性測(cè)試，都取得了良好結(jié)果。值得再提的是，應(yīng)用這套方法及軟件，王運(yùn)濤等人[32]確認(rèn)包含支撐裝置和靜氣動(dòng)彈性變形的CRM-WB構(gòu)型氣動(dòng)特性數(shù)值模擬結(jié)果，可以更加接近試驗(yàn)結(jié)果。

需要指出的是，不同于前述文獻(xiàn)用柔度矩陣方法或結(jié)構(gòu)模態(tài)分析方法去計(jì)算機(jī)翼靜氣彈變形，本文CSD分析直接采用SiPESC.FEMS 軟件進(jìn)行有限元計(jì)算。SiPESC是大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室自主開(kāi)發(fā)的“工程與科學(xué)計(jì)算軟件集成平臺(tái)”，其中的開(kāi)放式結(jié)構(gòu)有限元分析系統(tǒng)FEMS 是一款可與NASTRAN對(duì)標(biāo)的自主軟件。直接采用有限元計(jì)算可能更有利于柔性機(jī)翼的靜氣彈大變形計(jì)算，與CFD軟件TRIP耦合，對(duì)非線(xiàn)性氣動(dòng)力的情況也有較好的適用性。

1.3 優(yōu)化算法

針對(duì)十萬(wàn)變量的超大規(guī)模優(yōu)化問(wèn)題，智能尋優(yōu)算法存在“維數(shù)詛咒”，難以應(yīng)用。所謂大規(guī)模的智能尋優(yōu)算法，變量規(guī)模大多在一千或二千以下[33]，超過(guò)五千變量的極少。文獻(xiàn)[34]的合作協(xié)同進(jìn)化算法對(duì)一些五千變量的測(cè)試函數(shù)進(jìn)行優(yōu)化，函數(shù)計(jì)算次數(shù)高達(dá)百萬(wàn)量級(jí)。西安電子科大公茂果等人2005 年左右發(fā)展的免疫遺忘克隆規(guī)劃算法[35]，對(duì)一些超高維(一萬(wàn))測(cè)試函數(shù)的優(yōu)化，函數(shù)調(diào)用次數(shù)降到了十二萬(wàn)量級(jí)；而且估計(jì)若更高維使用其算法，函數(shù)維數(shù)每增加一維，函數(shù)計(jì)算次數(shù)平均增加接近25 次。由此估計(jì)十萬(wàn)變量的優(yōu)化也需要兩百多萬(wàn)次函數(shù)計(jì)算。若這些函數(shù)調(diào)用換做費(fèi)時(shí)的CFD計(jì)算，即便使用超算恐怕也難以勝任，所以針對(duì)這種問(wèn)題，使用智能尋優(yōu)算法必須想辦法降維。

優(yōu)化過(guò)程中不直接調(diào)用CFD軟件計(jì)算，而通過(guò)模型來(lái)預(yù)測(cè)，可減少計(jì)算量，即所謂基于代理模型的智能尋優(yōu)算法。其主要計(jì)算開(kāi)銷(xiāo)集中在建模需要的樣本計(jì)算上，增大設(shè)計(jì)變量的數(shù)量規(guī)模對(duì)構(gòu)建代理模型也是一項(xiàng)挑戰(zhàn)。一方面，隨設(shè)計(jì)變量增多，建模需要的樣本規(guī)模迅速增大；廣泛采用的拉丁超立方采樣(LHS)，為保證均勻性一般也要求樣本規(guī)模至少大于變量維數(shù)，樣本規(guī)模越大均勻性越好，需要開(kāi)展的CFD計(jì)算也越多，這對(duì)應(yīng)用超算進(jìn)行大規(guī)模并行計(jì)算提出了需求[36]。王超等人[37]在70 個(gè)設(shè)計(jì)變量的翼身組合體等氣動(dòng)優(yōu)化問(wèn)題中構(gòu)建支持向量機(jī)模型，發(fā)現(xiàn)在同為1 000 的樣本規(guī)模下，采用高維模型表示方法(HDMR)比采用LHS 得到的模型預(yù)測(cè)精度更好。似乎可以推斷同等模型精度要求下，HDMR方法的樣本規(guī)?？梢愿?，即使如此，估計(jì)樣本規(guī)模也要比變量維數(shù)大一個(gè)數(shù)量級(jí)以上。

另一方面，就算能完成這種大樣本計(jì)算，高維數(shù)據(jù)的建模也是一種挑戰(zhàn)。當(dāng)然這主要指構(gòu)建數(shù)據(jù)驅(qū)動(dòng)的代理模型，包括氣動(dòng)優(yōu)化領(lǐng)域廣泛采用的參數(shù)模型(如響應(yīng)面模型、徑向基模型)和非參數(shù)概率模型(如Kriging模型、高斯過(guò)程回歸模型)等各種數(shù)據(jù)擬合模型及降階模型(如基于正規(guī)正交分解的POD模型)，都存在維數(shù)災(zāi)難。近年來(lái)，各種高維統(tǒng)計(jì)降維方法的應(yīng)用，BP神經(jīng)網(wǎng)絡(luò)和支持向量機(jī)(SVM)的各種改進(jìn)，甚至深度神經(jīng)網(wǎng)絡(luò)(DNN)等大數(shù)據(jù)學(xué)習(xí)方法的流行，使得高維數(shù)據(jù)建模技術(shù)迅速發(fā)展。但所謂的高維建模，就作者所看到的而言，建模的變量數(shù)很少有超過(guò)二、三百的，而一個(gè)典型的氣動(dòng)優(yōu)化問(wèn)題設(shè)計(jì)變量大致在20 個(gè)～1 000 個(gè)左右[36]，差距是有的。所以，再考慮到樣本計(jì)算的成本，基于物理的代理模型[38]可能更有前途，雖然適用范圍不如數(shù)據(jù)驅(qū)動(dòng)模型寬，但可能更適合氣動(dòng)相關(guān)領(lǐng)域的超大規(guī)模優(yōu)化需要。這類(lèi)模型指高效的低階預(yù)測(cè)加修正/增強(qiáng)技術(shù)，即用簡(jiǎn)化物理模型或稀化計(jì)算網(wǎng)格預(yù)測(cè)，其結(jié)果與高保真預(yù)測(cè)的差量用與維數(shù)無(wú)關(guān)的修正/增強(qiáng)技術(shù)改善；故有人又稱(chēng)之為啟發(fā)式模型或多可信度、變可信度或變復(fù)雜度模型。

智能尋優(yōu)算法在超大規(guī)模優(yōu)化中存在的困難，好比“優(yōu)化十個(gè)設(shè)計(jì)變量需要五千次函數(shù)計(jì)算”，為此本文選用梯度尋優(yōu)算法，以期“優(yōu)化五千個(gè)設(shè)計(jì)變量?jī)H用十次函數(shù)計(jì)算”。但傳統(tǒng)梯度尋優(yōu)算法的計(jì)算開(kāi)銷(xiāo)與設(shè)計(jì)變量的數(shù)目仍然相關(guān)，雖然相對(duì)弱些，但還存在所謂的“擴(kuò)展問(wèn)題”。假設(shè)設(shè)計(jì)變量數(shù)目為N，Jameson 等基于最速捷線(xiàn)問(wèn)題的優(yōu)化研究，估計(jì)[39]最速下降法(具有整體收斂性；線(xiàn)性收斂速度)的計(jì)算開(kāi)銷(xiāo)為O(N2)，偽牛頓法(對(duì)凸函數(shù)具有整體收斂性；二次終止性，超線(xiàn)性收斂速度)的為O(N)，N很大的話(huà)，優(yōu)化過(guò)程將非常漫長(zhǎng)。Jameson提出的伴隨優(yōu)化方法之所以高效，不僅因?yàn)橄颂荻扔?jì)算量與設(shè)計(jì)變量數(shù)的相關(guān)性，而且還有一項(xiàng)關(guān)鍵，即使用Soblev梯度消除了尋優(yōu)算法對(duì)設(shè)計(jì)變量數(shù)的相關(guān)性[39]。所以，不能認(rèn)為解決了梯度計(jì)算量問(wèn)題，用任意多設(shè)計(jì)變量做大規(guī)模優(yōu)化就暢通無(wú)阻了。

在工程大規(guī)模優(yōu)化研究領(lǐng)域，Vanderplaats應(yīng)用其軟件產(chǎn)品BIGDOT在2002 年完成了5 萬(wàn)變量和5萬(wàn)約束的優(yōu)化驗(yàn)證[40]，實(shí)現(xiàn)了超過(guò)13.5 萬(wàn)個(gè)變量的結(jié)構(gòu)無(wú)約束拓?fù)鋬?yōu)化，可能是航空工程領(lǐng)域提倡十萬(wàn)變量?jī)?yōu)化的源頭，國(guó)內(nèi)相同規(guī)模的結(jié)構(gòu)優(yōu)化見(jiàn)羅利龍等人2016 年工作[41]。Vanderplaats從可靠性、效率和內(nèi)存需求角度評(píng)估了一系列新老算法之后，在BIGDOT中采用了受專(zhuān)利保護(hù)的改進(jìn)SUMT外點(diǎn)法[40]處理約束、FR共軛方向用作搜索方向、多項(xiàng)式插值進(jìn)行一維線(xiàn)搜索，這套算法的函數(shù)和梯度計(jì)算次數(shù)不隨設(shè)計(jì)變量數(shù)量增加而顯著變化，具有很好的擴(kuò)展性。其經(jīng)驗(yàn)表明，進(jìn)行超大規(guī)模的約束非線(xiàn)性工程優(yōu)化，算法要在保證可靠性的基礎(chǔ)上，削弱計(jì)算開(kāi)銷(xiāo)與設(shè)計(jì)變量的相關(guān)性。即減少內(nèi)存需求和避免耗時(shí)的子優(yōu)化或優(yōu)化方向計(jì)算。

因此，本文貼近工程應(yīng)用的靜氣彈優(yōu)化，考慮其計(jì)算開(kāi)銷(xiāo)較大，不追求優(yōu)化解的精度，選用整體收斂性較好的一種改進(jìn)最速下降法[42](SD with BFGSpreconditioned)，它對(duì)初值的依賴(lài)性較弱，用不多的內(nèi)存開(kāi)銷(xiāo)，對(duì)尋優(yōu)方向進(jìn)行一種預(yù)條件處理，圖4給出了該算法對(duì)某測(cè)試函數(shù)在十萬(wàn)變量下，與設(shè)計(jì)變量數(shù)目無(wú)關(guān)的優(yōu)良收斂特性。結(jié)構(gòu)尺寸優(yōu)化的特點(diǎn)是變量和約束多，計(jì)算開(kāi)銷(xiāo)相對(duì)較小，選用該領(lǐng)域發(fā)展出來(lái)并流行的移動(dòng)漸近線(xiàn)方法[43](MMA)。該方法構(gòu)造并求解一系列可分離的子優(yōu)化問(wèn)題，這些子優(yōu)化問(wèn)題是原目標(biāo)和約束函數(shù)的顯式單調(diào)嚴(yán)格凸近似，與設(shè)計(jì)變量數(shù)目無(wú)關(guān)、收斂性可調(diào)是MMA最為顯著的優(yōu)勢(shì)。在此基礎(chǔ)上發(fā)展的序列凸規(guī)劃法(SCP)，用于超過(guò)百萬(wàn)維的結(jié)構(gòu)拓?fù)鋬?yōu)化已得到驗(yàn)證[44]。

圖4 測(cè)試函數(shù)最小化收斂歷程

1.4 設(shè)計(jì)敏度計(jì)算方法

敏度分析研究模型的輸出如何隨輸入的變化而變化，主要用于基于梯度的優(yōu)化、不確定度量化、誤差分析等[45]。它有全局敏度和局部敏度之分；全局敏度分析在更為寬廣的取值范圍，量化相對(duì)輸入的輸出響應(yīng)，適合于有較大不確定度的模型；而局部敏度分析針對(duì)固定的一套輸入，量化輸出響應(yīng)，用于高確定度的物理模型。這里將CFD和CSD計(jì)算程序視作確定性模型討論局部敏度計(jì)算，偏導(dǎo)數(shù)計(jì)算是敏度分析的主要內(nèi)容，應(yīng)用中對(duì)二者大多不加嚴(yán)格區(qū)分。

對(duì)由計(jì)算機(jī)程序定義的函數(shù)，其導(dǎo)數(shù)計(jì)算目前有三類(lèi)方法：數(shù)值微分、自動(dòng)微分和解析方法。數(shù)值微分原理上基于泰勒級(jí)數(shù)展開(kāi)做近似計(jì)算，存在截?cái)嗾`差，用函數(shù)值或函數(shù)其它信息估算導(dǎo)數(shù)。主要包括有限差分法和復(fù)數(shù)步方法兩種，計(jì)算量都與設(shè)計(jì)變量數(shù)目相關(guān)。有限差分法估算導(dǎo)數(shù)的精度依賴(lài)擾動(dòng)步長(zhǎng)的選取，減小截?cái)嗾`差需要減小擾動(dòng)步長(zhǎng)，但一直減小擾動(dòng)步長(zhǎng)會(huì)增大減法抵消誤差(相近數(shù)相減帶來(lái)的有效數(shù)字損失)，對(duì)每個(gè)設(shè)計(jì)變量存在一個(gè)平衡這兩種誤差的最優(yōu)擾動(dòng)步長(zhǎng)。這種方法大多用在把計(jì)算模型當(dāng)做“黑箱”的情況，是梯度類(lèi)優(yōu)化軟件的缺省設(shè)置(如采用SQP方法的SNOPT)，當(dāng)函數(shù)計(jì)算存在噪聲時(shí)(來(lái)源于糟糕的收斂特性或網(wǎng)格變形)，足夠大的擾動(dòng)步長(zhǎng)有平滑作用，可抓住函數(shù)變化的正確趨勢(shì)[45]。復(fù)數(shù)步方法用復(fù)變量計(jì)算實(shí)函數(shù)的導(dǎo)數(shù)，需要修改計(jì)算模型的源程序，可用于任何算法以及高階導(dǎo)數(shù)計(jì)算，國(guó)外的應(yīng)用日益增多。因?yàn)闆](méi)有減法抵消誤差，擾動(dòng)步長(zhǎng)可以小到機(jī)器零，幾乎消除了截?cái)嗾`差，使得導(dǎo)數(shù)的計(jì)算精度與函數(shù)本身的計(jì)算精度相當(dāng)，自然比有限差分方法高。Martins等人[45]認(rèn)為復(fù)數(shù)步方法本質(zhì)上可以看成用運(yùn)算符重載的自動(dòng)微分，在極限上即擾動(dòng)步長(zhǎng)為零時(shí)，可得到與自動(dòng)微分方法相同的導(dǎo)數(shù)；而在實(shí)際有限精度的算法中，只要擾動(dòng)步長(zhǎng)足夠小，就可得到相同的結(jié)果。

自動(dòng)微分(AD)也叫算法微分或計(jì)算微分，是一種基于微分鏈?zhǔn)椒▌t，對(duì)計(jì)算機(jī)程序逐行進(jìn)行求導(dǎo)的技術(shù)；與解析方法對(duì)數(shù)學(xué)表達(dá)式求導(dǎo)類(lèi)似，但要用專(zhuān)門(mén)的軟件工具(如ADIFOR和ADIC等)改造源程序。因?yàn)闆](méi)有截?cái)嗾`差，能得到與解析方法同樣準(zhǔn)確的結(jié)果，但可自動(dòng)執(zhí)行更為簡(jiǎn)單。自動(dòng)微分有正向、反向兩種計(jì)算模式，計(jì)算量分別與設(shè)計(jì)變量數(shù)目、目標(biāo)函數(shù)數(shù)目相關(guān)。但像復(fù)數(shù)步方法一樣需要更多的內(nèi)存和計(jì)算的開(kāi)銷(xiāo)，比有限差分方法更為耗時(shí)。自動(dòng)微分方法計(jì)算的導(dǎo)數(shù)最為準(zhǔn)確，可用于解析方法中近似計(jì)算導(dǎo)數(shù)的精度校核。

這里所說(shuō)的解析方法，需要把計(jì)算機(jī)程序采用的算法返回到數(shù)學(xué)表達(dá)式上，再用鏈?zhǔn)椒▌t等方法推導(dǎo)，最后編成程序進(jìn)行計(jì)算。期間可采用符號(hào)微分軟件工具如Maple、Mathematica等，但更多地需要手工操作。這種推導(dǎo)要么在連續(xù)的微積分表達(dá)式上進(jìn)行，要么在離散的代數(shù)表達(dá)式上進(jìn)行，都帶來(lái)兩類(lèi)導(dǎo)數(shù)計(jì)算方法：直接法和伴隨法。前者導(dǎo)數(shù)計(jì)算量也與設(shè)計(jì)變量數(shù)目相關(guān)(但解的是線(xiàn)性方程，不同前述解非線(xiàn)性方程)，在目標(biāo)(或約束)函數(shù)比較多時(shí)效率高。后者相反，導(dǎo)數(shù)計(jì)算量?jī)H依賴(lài)目標(biāo)(或約束)函數(shù)的數(shù)目，在設(shè)計(jì)變量比較多時(shí)效率高。典型氣動(dòng)外形優(yōu)化問(wèn)題的變量多目標(biāo)少，選用伴隨方法計(jì)算導(dǎo)數(shù)效率最高。在某些結(jié)構(gòu)減重的尺寸優(yōu)化中，設(shè)計(jì)變量和目標(biāo)(或約束)函數(shù)的數(shù)目可能相當(dāng)，而且數(shù)目都很大，這兩種方法的優(yōu)勢(shì)都不突出，常用聚合約束減少其數(shù)量的方法放緩了優(yōu)化收斂歷程，需要發(fā)展新方法[45]。解析方法求導(dǎo)數(shù)比較復(fù)雜，需要詳細(xì)掌握計(jì)算模型采用的算法，編程調(diào)試開(kāi)發(fā)時(shí)間也比較長(zhǎng)。但因?yàn)橛?jì)算機(jī)條件的限制，開(kāi)發(fā)這種高效率計(jì)算方法的人還是很多。

理論上，解析方法求得的導(dǎo)數(shù)，其數(shù)值精度與原算法相同，而實(shí)際上沒(méi)有如此樂(lè)觀。主要原因是總導(dǎo)數(shù)即設(shè)計(jì)敏度的數(shù)學(xué)表達(dá)式可以解析推導(dǎo)出來(lái)，但是其中具體的各項(xiàng)偏導(dǎo)數(shù)并不是都能解析推導(dǎo)出來(lái)的，如通過(guò)交互操作生成的飛機(jī)復(fù)雜外形CFD計(jì)算網(wǎng)格相對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)等等。絕大多數(shù)的解析方法，用有限差分法計(jì)算這些偏導(dǎo)數(shù)，Dwight等人[46]認(rèn)為這是敏度誤差的主要來(lái)源，要提高計(jì)算精度，需要對(duì)每個(gè)設(shè)計(jì)變量逐個(gè)尋找有限差分法計(jì)算用的最優(yōu)擾動(dòng)步長(zhǎng)。這是件非常繁重的工作，通常統(tǒng)一采用一個(gè)步長(zhǎng)，由此帶來(lái)誤差。此外，在連續(xù)伴隨方法中，不是對(duì)任何目標(biāo)函數(shù)，都能得到連續(xù)伴隨方程的適定邊界條件：1)有可壓縮間斷時(shí)，推導(dǎo)連續(xù)伴隨方程存在連續(xù)可微的理論假設(shè)問(wèn)題；2)RANS 的連續(xù)伴隨方程推導(dǎo)，忽略了湍流模型的影響；3)數(shù)值求解連續(xù)伴隨方程，存在離散一致性問(wèn)題；4)伴隨敏度的網(wǎng)格收斂性問(wèn)題等等[46]。在離散伴隨方法中，控制方程左端項(xiàng)雅克比矩陣影響敏度計(jì)算精度，其完整、準(zhǔn)確的推導(dǎo)是件極其繁重的工作，引入簡(jiǎn)化近似假設(shè)會(huì)給敏度計(jì)算帶來(lái)誤差[46]。如大多數(shù)工作凍結(jié)渦團(tuán)黏性系數(shù)，忽略對(duì)湍流模型方程的線(xiàn)化，帶來(lái)了誤差，有些工作對(duì)各種湍流模型輸運(yùn)方程進(jìn)行了線(xiàn)化處理，但可能使最后迭代求解的線(xiàn)性方程條件數(shù)變差，帶來(lái)收斂問(wèn)題。Dwight等人還列舉了離散伴隨方程推導(dǎo)中，帶來(lái)誤差的其它一些近似方法[46]，包括：忽略對(duì)流項(xiàng)中人工耗散系數(shù)的導(dǎo)數(shù)、線(xiàn)化時(shí)僅使用一階對(duì)流通量、線(xiàn)化簡(jiǎn)化的黏性通量、完全忽略黏性項(xiàng)以及最極端的做法，敏度計(jì)算中徹底忽略氣動(dòng)的貢獻(xiàn)，只有網(wǎng)格相對(duì)設(shè)計(jì)變量導(dǎo)數(shù)的影響，等等。用基于這些近似推導(dǎo)的雅克比計(jì)算的敏度與準(zhǔn)確結(jié)果相比，有些相對(duì)誤差可高達(dá)90%，以致于是否可用受到質(zhì)疑。

復(fù)雜解析方法的敏度計(jì)算精度評(píng)估，大多根據(jù)其計(jì)算結(jié)果與有限差分方法結(jié)果的對(duì)比情況進(jìn)行判斷，但是后者作為標(biāo)準(zhǔn)也受質(zhì)疑，所以有人在稍小的問(wèn)題中，采用復(fù)數(shù)步方法或AD方法的計(jì)算結(jié)果作為評(píng)估標(biāo)準(zhǔn)。敏度在優(yōu)化算法中當(dāng)作梯度用于構(gòu)造尋優(yōu)方向，其精度對(duì)優(yōu)化的速度和精度有重要影響，越高越好。然而，Dwight等人[46]指出，對(duì)一些健壯的優(yōu)化算法，如共軛梯度法和最速下降法，用一些精度較差的敏度(各種近似的離散伴隨方法的結(jié)果)，在數(shù)量相當(dāng)?shù)牡綌?shù)內(nèi)，也能得到與用高精度敏度相同的優(yōu)化結(jié)果，但對(duì)BFGS 偽牛頓算法則不行。也就是說(shuō)，離散伴隨方法中雅克比的大多數(shù)近似，是否可用與問(wèn)題相關(guān)，在對(duì)敏度計(jì)算精度沒(méi)把握時(shí)，采用簡(jiǎn)單健壯的優(yōu)化算法是有利的，這同前述解決大規(guī)模優(yōu)化問(wèn)題的指向一致。

在設(shè)計(jì)敏度的計(jì)算上，為避免有限差分法帶來(lái)的計(jì)算量問(wèn)題，過(guò)去三十年里，先是以連續(xù)伴隨方法的研究為主，后來(lái)離散伴隨方法的研究增多。二者都能給出滿(mǎn)足需要的敏度精度。因?yàn)榈湫桶殡S敏度的主要誤差來(lái)源是計(jì)算網(wǎng)格相對(duì)設(shè)計(jì)變量的導(dǎo)數(shù)，所以連續(xù)和離散哪種方法更好，取決于具體推導(dǎo)和編程實(shí)施[46]。目前AD技術(shù)在離散伴隨方法偏導(dǎo)數(shù)的計(jì)算中取得了特別好的效果[45]，二者的深度結(jié)合可能是一個(gè)發(fā)展方向[46]。

然而，高性能計(jì)算技術(shù)的提高，尤其在E級(jí)計(jì)算環(huán)境下，設(shè)計(jì)敏度有限差分法計(jì)算量的問(wèn)題，即使在工程應(yīng)用上，也可以預(yù)期將大為緩解。首先，CFD分析已廣泛采用大規(guī)模并行計(jì)算；其次，設(shè)計(jì)敏度用有限差分法計(jì)算，相互間沒(méi)有信息交換，完全可以通過(guò)并行計(jì)算縮短時(shí)間[36，47]。雖然計(jì)算效率較低，沒(méi)有在理論上解決敏度計(jì)算量相對(duì)于設(shè)計(jì)變量數(shù)目的擴(kuò)展性問(wèn)題；但是可把各種軟件當(dāng)做“黑箱”，只有一個(gè)簡(jiǎn)單地步長(zhǎng)優(yōu)化問(wèn)題，上手快，應(yīng)用簡(jiǎn)便靈活，適用范圍廣，特別適合工程多學(xué)科問(wèn)題。原來(lái)在串行計(jì)算環(huán)境里，采用差分敏度方法，計(jì)算時(shí)間過(guò)于冗長(zhǎng)，但在目前及以后的大規(guī)模并行計(jì)算環(huán)境里，資源不是大問(wèn)題，敏度計(jì)算時(shí)間完全可以縮短到一遍流場(chǎng)計(jì)算時(shí)間的量級(jí)，情況與原來(lái)有了很大不同。相比發(fā)展復(fù)雜的氣動(dòng)-結(jié)構(gòu)耦合伴隨敏度計(jì)算方法[48-49]，差分敏度更易在短期內(nèi)得到成功實(shí)施。為此，本文在氣彈并行計(jì)算軟件的基礎(chǔ)上，采用并行的有限差分方法計(jì)算設(shè)計(jì)敏度。

2 優(yōu)化問(wèn)題描述

高速大展弦比機(jī)翼的典型優(yōu)化設(shè)計(jì)問(wèn)題是最大化航程/航時(shí)。除發(fā)動(dòng)機(jī)耗油率外，主要涉及氣動(dòng)和結(jié)構(gòu)兩個(gè)學(xué)科，需要?dú)鈩?dòng)效率最大化，結(jié)構(gòu)重量最小化。二者對(duì)機(jī)翼展向載荷分布、展弦比、相對(duì)厚度等參數(shù)的要求迥異，需采用優(yōu)化技術(shù)來(lái)定量協(xié)調(diào)。目前這兩個(gè)學(xué)科各自的優(yōu)化設(shè)計(jì)技術(shù)都很成熟，希望進(jìn)一步提高設(shè)計(jì)技術(shù)的能力，必須在設(shè)計(jì)初期進(jìn)行學(xué)科綜合，以期獲得整體最優(yōu)的協(xié)同效應(yīng)。然而，各學(xué)科自身的優(yōu)化設(shè)計(jì)問(wèn)題，已是一個(gè)復(fù)雜的多目標(biāo)多約束的非線(xiàn)性?xún)?yōu)化問(wèn)題，如何綜合存在挑戰(zhàn)。例如，氣動(dòng)外形優(yōu)化要在來(lái)自結(jié)構(gòu)設(shè)計(jì)的厚度、容積等幾何約束下，最大化巡航設(shè)計(jì)點(diǎn)的升阻比，最大化起降設(shè)計(jì)點(diǎn)的最大升力系數(shù)，還要在巡航馬赫數(shù)下滿(mǎn)足對(duì)抖振升力系數(shù)、俯仰力矩系數(shù)的約束，以及巡航升力系數(shù)下對(duì)阻力發(fā)散馬赫數(shù)的約束，等等。結(jié)構(gòu)尺寸優(yōu)化，就機(jī)翼減重而言，要考慮各種機(jī)動(dòng)、陣風(fēng)、滑行等多種載荷下的強(qiáng)度、剛度、疲勞壽命和可靠性等要求，還要考慮來(lái)自靜/動(dòng)氣彈分析的穩(wěn)定性約束。把兩個(gè)學(xué)科實(shí)際設(shè)計(jì)目標(biāo)和約束過(guò)度簡(jiǎn)化/忽略的綜合優(yōu)化，實(shí)用價(jià)值可能受限。

為充分發(fā)揮各學(xué)科的專(zhuān)業(yè)優(yōu)勢(shì)，本文采用學(xué)科分治的思路進(jìn)行綜合優(yōu)化，同時(shí)也貼近工業(yè)部門(mén)的組織分工和設(shè)計(jì)流程。技術(shù)實(shí)現(xiàn)上做了三件事：1)氣動(dòng)外形的參數(shù)化模型和結(jié)構(gòu)有限元的參數(shù)化模型，適時(shí)匹配(見(jiàn)圖5)，方便載荷、變形信息的準(zhǔn)確交換；2)基于該匹配模型，用TRIP+SiPESC耦合軟件計(jì)算靜氣彈平衡的氣動(dòng)力和載荷，分別用于氣動(dòng)優(yōu)化和結(jié)構(gòu)優(yōu)化；3)氣動(dòng)優(yōu)化輸出的匹配模型及其對(duì)應(yīng)的機(jī)動(dòng)載荷，是結(jié)構(gòu)優(yōu)化的輸入，后者更新有限元模型的尺寸屬性再返回氣動(dòng)優(yōu)化，之后循環(huán)直至收斂。形成的綜合優(yōu)化框架如圖6 所示，其中左右兩邊分別為氣動(dòng)優(yōu)化和結(jié)構(gòu)優(yōu)化，中間是有限元模型與氣動(dòng)外形的匹配合成，通過(guò)它及耦合靜氣彈計(jì)算實(shí)現(xiàn)兩學(xué)科優(yōu)化的綜合。

圖5 有限元模型與氣動(dòng)外形的匹配

圖6 氣動(dòng)結(jié)構(gòu)綜合優(yōu)化框架

這樣進(jìn)行氣動(dòng)-結(jié)構(gòu)優(yōu)化的綜合，不僅可以發(fā)揮各學(xué)科做細(xì)做實(shí)的專(zhuān)業(yè)優(yōu)勢(shì)，而且有利于各學(xué)科根據(jù)自己?jiǎn)栴}的特點(diǎn)，選用合適的優(yōu)化算法，更有效地解決問(wèn)題。而不是籠統(tǒng)地采用一種優(yōu)化算法，給困難的非線(xiàn)性約束優(yōu)化算法平添更多麻煩。比較而言，氣動(dòng)優(yōu)化的特點(diǎn)是函數(shù)計(jì)算開(kāi)銷(xiāo)大，結(jié)構(gòu)優(yōu)化的特點(diǎn)是變量和約束多。所以，這里氣動(dòng)優(yōu)化用無(wú)約束優(yōu)化算法，其中幾何約束在參數(shù)化建模時(shí)實(shí)施，不對(duì)違反約束的外形進(jìn)行CFD計(jì)算，減少計(jì)算開(kāi)銷(xiāo)，來(lái)自各設(shè)計(jì)點(diǎn)的設(shè)計(jì)目標(biāo)，采用線(xiàn)性加權(quán)的評(píng)價(jià)函數(shù)方法，轉(zhuǎn)化成單目標(biāo)優(yōu)化問(wèn)題減少計(jì)算量。而結(jié)構(gòu)優(yōu)化則用大規(guī)模變量的非線(xiàn)性約束優(yōu)化算法，可對(duì)有限元模型每個(gè)單元提出設(shè)計(jì)變量和約束，采用靜氣彈平衡的氣動(dòng)載荷，而不是其它大多數(shù)工作采用的剛體外形氣動(dòng)載荷，理論上可使結(jié)構(gòu)優(yōu)化更為精細(xì)可靠。如此，也實(shí)現(xiàn)了設(shè)計(jì)變量的分組，緩解了大規(guī)模設(shè)計(jì)變量帶來(lái)的困難。

2.1 氣動(dòng)優(yōu)化問(wèn)題

選用NASA的CRM 標(biāo)模[50]機(jī)翼作為研究對(duì)象。這款公開(kāi)的標(biāo)模具有跨聲速遠(yuǎn)程寬體民機(jī)的典型氣動(dòng)設(shè)計(jì)特征，不僅提供了機(jī)翼巡航時(shí)的名義1-g外形，還提供了相應(yīng)的結(jié)構(gòu)有限元模型。其巡航設(shè)計(jì)點(diǎn)為：

原始機(jī)翼的參考面積(Wimpress面積[50])為4 130 ft2(1 ft=0.304 8 m)，參考弦長(zhǎng)是不考慮Yehudi的梯形翼平均氣動(dòng)弦長(zhǎng)275.8 in(1 in=0.025 4 m)，力矩參考點(diǎn)取該平均氣動(dòng)弦長(zhǎng)的前1/4 弦長(zhǎng)點(diǎn)。該梯形翼四分之一弦線(xiàn)后掠角為35°，展弦比為9.0，梢根比0.275。原始機(jī)翼展長(zhǎng)b 為192.8 ft，翼身交接在10%*b/2 附近，37%*b/2 站位為后緣Yehudi拐折。展向有21 個(gè)定義截面，除拐折外其余按5%*b/2 的間隔均勻分布；由縮比模型后緣厚度不能小于0.014 in 的制造約束，導(dǎo)出翼尖后緣相對(duì)厚度為0.48%c(c是當(dāng)?shù)叵议L(zhǎng))。CRM外露翼的截面最大厚度，翼根10%*b/2 處約65 in，拐折37%*b/2 處約30 in，翼尖約10 in，最大相對(duì)厚度分別為13.8%、10.5%和9.5%。CRM外露翼的截面最大相對(duì)彎度，翼根10%*b/2處約為0.1%，拐折37%*b/2 處約為0.95%，并開(kāi)始采用超臨界翼型；大的彎度集中在60%*b/2 到90%*b/2 之間，其最大相對(duì)彎度在1.5%到1.6%之間，之后往翼尖降到0.09%。

CRM機(jī)翼展向幾何扭轉(zhuǎn)角，總的外洗超過(guò)8°。有兩段大的扭轉(zhuǎn)，估計(jì)用于改善失速特性，見(jiàn)圖7 所示，一段從翼根的4.44°降到拐折附近的0.76°，另一段在翼梢，從90%*b/2 附近的-2.57°降到翼尖的-3.75°，兩段之間布置的超臨界翼型外洗在2°以?xún)?nèi)。根據(jù)這個(gè)特征，定義了4 個(gè)設(shè)計(jì)變量來(lái)控制幾何扭轉(zhuǎn)分布，即圖7 中的DV1～DV4。對(duì)CRM模型名義1-g外形來(lái)說(shuō)，以上幾何扭轉(zhuǎn)包含了靜氣彈變形。機(jī)翼上反方面也包含了靜氣彈彎曲，相對(duì)翼根前緣點(diǎn)，翼尖前緣點(diǎn)非線(xiàn)性地向上抬起了約88 in，相當(dāng)于7.6%的半展長(zhǎng)高度，估計(jì)在線(xiàn)性靜氣彈的變形范圍。

圖7 機(jī)翼幾何扭轉(zhuǎn)分布的設(shè)計(jì)參數(shù)定義

圖8 給出了CRM機(jī)翼無(wú)彎扭的平面機(jī)翼。本文假設(shè)它是初始型架外形，期望通過(guò)優(yōu)化得到其展向幾何扭轉(zhuǎn)分布，疊加上氣彈扭轉(zhuǎn)變形，滿(mǎn)足巡航氣動(dòng)效率最大化。在前述設(shè)計(jì)點(diǎn)的固定升力系數(shù)下，最小化選用的目標(biāo)函數(shù)為：

圖8 CRM機(jī)翼參數(shù)化無(wú)彎扭外形與1-g外形對(duì)比

其中ω1，ω2為權(quán)因子，CD，Cm分別為阻力系數(shù)和俯仰力矩系數(shù)，CD0，Cm0是這兩系數(shù)的初始值，用于規(guī)范數(shù)值大小。這里引入俯仰力矩系數(shù)旨在減小配平對(duì)氣動(dòng)效率的影響。

包括阻力系數(shù)和俯仰力矩系數(shù)在內(nèi)的氣動(dòng)力系數(shù)，通過(guò)氣彈耦合程序TRIP+SiPESC計(jì)算得到。其耦合方法前文已介紹，在CFD方面，TRIP采用SA湍流模型計(jì)算全RANS 方程，對(duì)流項(xiàng)離散采用MUSCL-Roe方法，時(shí)間推進(jìn)用隱式LU-SGS 迭代。氣動(dòng)系數(shù)的無(wú)量綱參考量和力矩參考點(diǎn)與CRM機(jī)翼一致，不再贅述。

2.2 結(jié)構(gòu)優(yōu)化問(wèn)題

以NASA的CRM-V15-wingbox-1.bdf全尺寸半展長(zhǎng)翼盒有限元模型(FEM)為基礎(chǔ)，其結(jié)構(gòu)布局是標(biāo)準(zhǔn)設(shè)計(jì)手冊(cè)提供的傳統(tǒng)方案，包括兩主梁和蒙皮，以及50 多翼肋，以及桁條、翼肋緣條和加強(qiáng)筋等部件。首先，考慮氣動(dòng)載荷的傳遞需要，在前后緣操縱面區(qū)域添加了模擬蒙皮和翼肋的面單元，替換原來(lái)的剛體單元RBE3 和集中質(zhì)量單元CONM2，如圖9 所示，然后，根據(jù)前述氣動(dòng)模型除去彎扭建立的表面外形，調(diào)整節(jié)點(diǎn)坐標(biāo)與氣動(dòng)翼面一致。整個(gè)FEM模型有20 451 個(gè)節(jié)點(diǎn)，42 477 個(gè)單元，包括12 980個(gè)簡(jiǎn)單梁?jiǎn)卧狢BAR、29 470 個(gè)四邊形板單元CQUAD4 和少數(shù)三角形板單元CTRIA3 三類(lèi)。原CRM有限元模型材料采用各向同性的鋁合金，楊氏模量E=10.3 ×106psi，密度ρ=0.101 lb/in3。翼根有位移約束的節(jié)點(diǎn)分布如圖10 所示，是有限元計(jì)算的固支邊界約束。氣動(dòng)載荷加載位置的分布如圖11 所示，做靜氣彈的線(xiàn)性分析，作用力的方向假設(shè)不變。加載的2.5-g氣動(dòng)載荷通過(guò)氣彈耦合程序TRIP+SiPESC計(jì)算得到，計(jì)算條件是：

圖9 前后緣替換剛體和集中質(zhì)量單元(紅)的面單元

圖10 翼根位移約束節(jié)點(diǎn)的分布

圖11 氣動(dòng)載荷加載位置的分布

結(jié)構(gòu)尺寸優(yōu)化設(shè)計(jì)參數(shù)，計(jì)劃取殼單元屬性PSHELL的厚度參數(shù)，以及梁的截面尺寸參數(shù)，后者影響梁?jiǎn)卧獙傩訮BAR的面積、慣性矩等等。本文取了10 000 個(gè)四邊形板單元CQUAD4 的殼單元屬性PSHELL中的厚度參數(shù)，作尺寸設(shè)計(jì)變量。優(yōu)化約束是最大應(yīng)力和最大位移，通過(guò)有限元計(jì)算軟件SiPESC.FEM 得到，它兼容商業(yè)軟件NASTRAN的bdf輸入文件格式。最小化結(jié)構(gòu)重量為優(yōu)化目標(biāo)，密度不變情況下，通過(guò)體積計(jì)算值表征。

3 初步優(yōu)化結(jié)果

雖然基于梯度的優(yōu)化算法比無(wú)需梯度的算法難用，但對(duì)設(shè)計(jì)變量在O(102)量級(jí)或更多的優(yōu)化問(wèn)題，它們是目前唯一的途徑[45]。然而，對(duì)氣動(dòng)—結(jié)構(gòu)兩學(xué)科、甚至更多學(xué)科的復(fù)雜耦合系統(tǒng)來(lái)講，設(shè)計(jì)敏度的高效計(jì)算也是一件困難的事。除發(fā)展先進(jìn)的算法外，利用超級(jí)計(jì)算機(jī)提供的高性能并行計(jì)算能力，也是一種不應(yīng)該被忽視的解決方案，為此，本文提出采用差分敏度并行計(jì)算策略代替早年的串行計(jì)算方案。為證明其可行性，這里采用小規(guī)模的優(yōu)化設(shè)計(jì)算例驗(yàn)證。初步以CRM單個(gè)機(jī)翼而不是翼身組合體為對(duì)象，減小CFD網(wǎng)格規(guī)模(單獨(dú)機(jī)翼500萬(wàn)左右)和計(jì)算時(shí)間，便于調(diào)試。初步以鋼作為結(jié)構(gòu)材料，其楊氏模量為E=30.5 ×106psi，避免鋁合金在2.5 g載荷下變形大帶來(lái)的幾何非線(xiàn)性問(wèn)題；因?yàn)楸疚牟捎渺o氣彈計(jì)算的可靠性，目前還局限在線(xiàn)性變形范圍內(nèi)，希望最大位移控制在10%*b/2內(nèi)，對(duì)CRM機(jī)翼來(lái)說(shuō)為115 in，最大應(yīng)力不大于材料的極限，估計(jì)為3.5E5 psi。初步以前述4 個(gè)展向幾何扭轉(zhuǎn)控制參數(shù)作設(shè)計(jì)變量，沒(méi)有用更多的型面控制參數(shù)作設(shè)計(jì)變量去減小波阻，旨在減少計(jì)算量便于程序調(diào)試。初步以指定范圍單元的最大應(yīng)力、最大位移為約束，不對(duì)所有單元都提約束，降低對(duì)一萬(wàn)變量約束優(yōu)化算法的挑戰(zhàn)。

在氣動(dòng)中心計(jì)算所，采用六百個(gè)核耗時(shí)三天完成上述算例計(jì)算。圖12 的氣動(dòng)結(jié)構(gòu)兩輪綜合優(yōu)化歷程顯示，在不降低氣動(dòng)效率，滿(mǎn)足最大應(yīng)力和位移約束下，實(shí)現(xiàn)減重9%，俯仰力矩系數(shù)改善近40%，驗(yàn)證了方法及軟件的有效性。圖12(a)說(shuō)明1 g載荷靜氣彈優(yōu)化過(guò)程中固定了升力系數(shù)。注意在一輪結(jié)構(gòu)優(yōu)化后，做第二輪氣彈優(yōu)化時(shí)升阻比不升反降。對(duì)比圖12(b)可發(fā)現(xiàn)，第二輪氣彈優(yōu)化過(guò)程中，俯仰力矩系數(shù)絕對(duì)值在迅速減小。由于這里優(yōu)化的評(píng)價(jià)函數(shù)是阻力系數(shù)和俯仰力矩系數(shù)的加權(quán)組合(ω1=0.8，ω2=0.2)，說(shuō)明這時(shí)俯仰力矩系數(shù)絕對(duì)值的減小，主導(dǎo)了評(píng)價(jià)函數(shù)的減小。圖12(c)給出了兩輪2.5 g載荷下的結(jié)構(gòu)優(yōu)化歷程?？梢钥闯?，主要的改善來(lái)自第一輪優(yōu)化的前兩步。圖中也給出了最大應(yīng)力、位移滿(mǎn)足約束的數(shù)值。最大應(yīng)力2.5E5 psi小于約束，最大位移83.9 in 也小于我們線(xiàn)性變形的位移約束。

圖12 氣動(dòng)-結(jié)構(gòu)兩輪綜合優(yōu)化歷程

圖13 介紹型架外形幾何扭轉(zhuǎn)優(yōu)化后，得到的1-g靜氣彈變形情況。圖13(a)展示了巡航設(shè)計(jì)點(diǎn)的表面壓力分布，沒(méi)有優(yōu)化型面參數(shù)還存在激波。從圖13(b)可看出這時(shí)翼尖有約40 in 的位移，是加載2.5-g載荷最大位移的一半左右，定性合理。圖13(c)中紅色曲線(xiàn)是1-g構(gòu)型的幾何扭轉(zhuǎn)分布，它是型架外形幾何扭轉(zhuǎn)分布(藍(lán)色曲線(xiàn))和靜氣彈扭轉(zhuǎn)變形(綠色曲線(xiàn))的疊加結(jié)果。85%*b/2 處有近4°的外洗，與原CRM機(jī)翼相當(dāng)。但向翼尖方向扭轉(zhuǎn)角減小，與一般的扭轉(zhuǎn)分布有些不同，這是由于本文定義的扭轉(zhuǎn)設(shè)計(jì)變量DV3 和DV4(見(jiàn)圖7)的敏度不同所致，似乎有些翼梢小翼的作用，用一個(gè)設(shè)計(jì)變量控制整個(gè)外段扭轉(zhuǎn)可以避免這種情況。這里的扭轉(zhuǎn)分布，在翼梢、翼根處與CRM機(jī)翼扭轉(zhuǎn)分布差別較大，僅照顧了1-g飛行氣動(dòng)效率的需要，下一步需要在翼身組合體上，同時(shí)兼顧改善失速特性的需要進(jìn)行優(yōu)化。

圖13 巡航外形的靜氣彈變形

4 結(jié)論

面向工程應(yīng)用的高保真靜氣彈優(yōu)化設(shè)計(jì)，可采用先進(jìn)算法實(shí)現(xiàn)，有高性能計(jì)算技術(shù)的支持，也可采用簡(jiǎn)單健壯的方法實(shí)現(xiàn)。通過(guò)本文基于CFD的優(yōu)化方法介紹，圍繞CRM機(jī)翼開(kāi)展的稍小規(guī)模算例測(cè)試分析，可以得到以下兩條結(jié)論：

1)差分敏度采用大規(guī)模并行技術(shù)計(jì)算，可避免原來(lái)差分敏度串行計(jì)算時(shí)間開(kāi)銷(xiāo)過(guò)大的問(wèn)題，而與耦合伴隨方法相當(dāng)，進(jìn)而快速解決氣動(dòng)-結(jié)構(gòu)高保真耦合分析的敏度計(jì)算難題；

2)采用健壯的預(yù)條件最速下降法和移動(dòng)漸近線(xiàn)法優(yōu)化，不僅整體收斂性好，而且可使優(yōu)化循環(huán)步數(shù)與設(shè)計(jì)變量數(shù)目無(wú)關(guān)，能避免大規(guī)模設(shè)計(jì)變量帶來(lái)的計(jì)算開(kāi)銷(xiāo)過(guò)大問(wèn)題。

本文基于TRIP和SiPESC軟件建立的方法具備向更大規(guī)模氣動(dòng)結(jié)構(gòu)綜合優(yōu)化發(fā)展的潛力。下一步，將應(yīng)用已獲得初步調(diào)試結(jié)果的幾何非線(xiàn)性氣彈計(jì)算代碼，針對(duì)鋁合金和復(fù)合材料機(jī)翼的翼身組合體構(gòu)型，進(jìn)行多點(diǎn)優(yōu)化設(shè)計(jì)研究。