薛幫猛, 張文升, 張志雄

(中國(guó)商飛北京民用飛機(jī)技術(shù)研究中心 民用飛機(jī)設(shè)計(jì)數(shù)字仿真北京市重點(diǎn)實(shí)驗(yàn)室，北京 102211)

0 引 言

飛機(jī)/發(fā)動(dòng)機(jī)集成(下文簡(jiǎn)稱(chēng)：飛發(fā)集成)是翼吊布局飛機(jī)研制中的關(guān)鍵技術(shù)。隨著涵道比和尺寸的增大，發(fā)動(dòng)機(jī)短艙相對(duì)機(jī)翼的近距耦合安裝不可避免，發(fā)動(dòng)機(jī)和機(jī)翼間的氣動(dòng)干擾問(wèn)題加劇。研究表明[1]：發(fā)動(dòng)機(jī)短艙/吊掛對(duì)內(nèi)、外兩側(cè)機(jī)翼前緣附近的“上洗”氣流分別產(chǎn)生增強(qiáng)和削弱的效果。國(guó)內(nèi)多位學(xué)者也對(duì)飛發(fā)集成構(gòu)型的復(fù)雜流場(chǎng)特征開(kāi)展了分析研究[2-4]。實(shí)際應(yīng)用中，安裝效應(yīng)對(duì)機(jī)翼壓力分布和展向載荷分布的影響范圍超過(guò)50%翼展[1,5]。因此，在機(jī)翼設(shè)計(jì)中直接考慮短艙/吊掛的干擾作用變得很有必要。德國(guó)宇航院Hohciscl和巴西航空工業(yè)公司Guilherme等分別在文獻(xiàn)[5]和文獻(xiàn)[6]中詳盡闡述了在飛發(fā)集成的環(huán)境中開(kāi)展部件氣動(dòng)設(shè)計(jì)的重要性，Bousquet[7]總結(jié)了ONERA在發(fā)動(dòng)機(jī)集成方面的風(fēng)洞試驗(yàn)研究。

在設(shè)計(jì)工具方法方面，由于跨聲速階段阻力對(duì)幾何外形變化的響應(yīng)過(guò)于敏感，依靠人工修形的經(jīng)驗(yàn)設(shè)計(jì)方法難以最大限度地挖掘減阻潛力。搭建自動(dòng)的優(yōu)化系統(tǒng)開(kāi)展氣動(dòng)外形優(yōu)化設(shè)計(jì)是近年來(lái)研究的熱點(diǎn)領(lǐng)域，主要涉及幾何外形參數(shù)化、計(jì)算網(wǎng)格自動(dòng)生成、氣動(dòng)性能評(píng)估、優(yōu)化搜索算法等方面。東京大學(xué)Hyoung-Jin等[8-10]使用伴隨梯度優(yōu)化方法開(kāi)展了針對(duì)機(jī)翼/機(jī)身/短艙/吊掛構(gòu)型(下文簡(jiǎn)稱(chēng)：WBNP構(gòu)型)的多點(diǎn)優(yōu)化設(shè)計(jì)工作，CFD計(jì)算用的是Euler求解器。在WBNP復(fù)雜外形幾何參數(shù)化和網(wǎng)格變形方法方面，Takanaka[11]、Truong[12]和左英桃[13]等都做了大量有意義的研究。胡曉東[14]等在民機(jī)短艙安裝位置優(yōu)化中得出了安裝位置對(duì)全機(jī)氣動(dòng)特性的影響規(guī)律。清華大學(xué)張宇飛等[15]嘗試了運(yùn)用多目標(biāo)優(yōu)化算法在短艙干擾下開(kāi)展機(jī)翼一體化設(shè)計(jì)，有較高的工程應(yīng)用價(jià)值，但由于計(jì)算量巨大，設(shè)計(jì)周期是個(gè)問(wèn)題。

為了實(shí)現(xiàn)在飛發(fā)集成構(gòu)型中直接優(yōu)化機(jī)翼，本文建立了基于直接RANS計(jì)算分析和遺傳算法尋優(yōu)的多目標(biāo)優(yōu)化系統(tǒng)，針對(duì)NASA CRM的WBNP構(gòu)型開(kāi)展了機(jī)翼多目標(biāo)優(yōu)化設(shè)計(jì)工作。在翼型CST參數(shù)化方法的基礎(chǔ)上發(fā)展了一種簡(jiǎn)捷高效的三維機(jī)翼參數(shù)化描述方法，并將幾何參數(shù)化與網(wǎng)格變形集成在一起，避免了調(diào)用CAD軟件建模和網(wǎng)格生成軟件重新生成計(jì)算網(wǎng)格的復(fù)雜CFD前處理過(guò)程，從而達(dá)到了簡(jiǎn)化優(yōu)化流程、提高魯棒性的目的。此外，實(shí)現(xiàn)了在“天河2號(hào)”超級(jí)計(jì)算機(jī)上同時(shí)對(duì)超過(guò)100個(gè)方案實(shí)施分布式并行計(jì)算分析，再加上適當(dāng)降低收斂標(biāo)準(zhǔn)和一系列CFD加速收斂措施的合理使用，較好地解決了優(yōu)化周期問(wèn)題。

1 優(yōu)化方法

1.1 幾何外形參數(shù)化

民用飛機(jī)機(jī)翼氣動(dòng)設(shè)計(jì)一般針對(duì)圖1所示的控制剖面展開(kāi)，最終的機(jī)翼幾何模型則在CATIA或其他CAD軟件中基于這些控制剖面放樣生成。本文以文獻(xiàn)[16]中二維CST參數(shù)化方法為基礎(chǔ)，發(fā)展了一種簡(jiǎn)捷高效的機(jī)翼參數(shù)化方法。

圖1 機(jī)翼外形參數(shù)化Fig.1 Wing shape parameterization

對(duì)于每個(gè)控制翼型，上表面或下表面的CST參數(shù)化表達(dá)式為：

(1)

其中弦向相對(duì)坐標(biāo)：

ψ=(x-xLE)/c

(2)

法向相對(duì)坐標(biāo)：

ζ=(y-yLE)/c

(3)

后緣開(kāi)口相對(duì)大小：

ζT=ΔyTE/c

(4)

類(lèi)型函數(shù)：

(5)

形狀函數(shù)由n階Bernstein多項(xiàng)式疊加構(gòu)成：

(6)

wr為第r項(xiàng)的權(quán)重系數(shù)，不同的權(quán)重系數(shù)組合描述不同形狀的翼型，也就是翼型的CST參數(shù)。已知翼型的m個(gè)離散點(diǎn)坐標(biāo)(m>n)，使用最小二乘逼近，求解n+1階線性代數(shù)方程組即可得到描述該翼型的n階CST參數(shù)。

這里將機(jī)翼表面視為由沿展向的無(wú)限多個(gè)剖面翼型構(gòu)成，任一展向位置z(y指向上，x指向后)處剖面參數(shù)由控制剖面參數(shù)插值得到。采用如下插值方法：機(jī)翼剖面翼型的前緣點(diǎn)y坐標(biāo)yLE、當(dāng)?shù)嘏まD(zhuǎn)角ΔαT和翼型n階CST參數(shù)(w0,w1，…，wn)由控制剖面以三次樣條插值得到；前緣點(diǎn)xLE、當(dāng)?shù)叵议L(zhǎng)c、后緣開(kāi)度ζT則由附近兩個(gè)控制剖面線性插值得到。這就是本文提出的機(jī)翼三維參數(shù)化方法。給定機(jī)翼平面形狀內(nèi)任一位置(x,z)，可以唯一地算出上、下表面y坐標(biāo)yU、yL。

本文的優(yōu)化設(shè)計(jì)中，設(shè)計(jì)變量是控制剖面的扭轉(zhuǎn)角或CST參數(shù)。以上下表面同步連動(dòng)的方式，實(shí)現(xiàn)了厚度不變、僅改變剖面中弧線的優(yōu)化設(shè)計(jì)。設(shè)計(jì)變量以擾動(dòng)量的形式疊加到描述初始外形(以下稱(chēng)Baseline外形)的基本參數(shù)上，基本參數(shù)需在優(yōu)化開(kāi)始前對(duì)所有控制剖面經(jīng)過(guò)參數(shù)反算獲得。

1.2 計(jì)算網(wǎng)格與網(wǎng)格變形方法

Baseline外形的計(jì)算網(wǎng)格在ANSYS ICEM-CFD軟件中生成(如圖2所示)，由239塊結(jié)構(gòu)化網(wǎng)格構(gòu)成，總網(wǎng)格單元數(shù)約800萬(wàn)。相應(yīng)的機(jī)翼表面網(wǎng)格由110片構(gòu)成。Baseline外形的網(wǎng)格將作為網(wǎng)格變形的輸入。本文的網(wǎng)格變形方法是與1.1節(jié)的機(jī)翼參數(shù)化方法集成在一起的。在給定一組機(jī)翼參數(shù)后，對(duì)于Baseline機(jī)翼表面的每一個(gè)網(wǎng)格點(diǎn)，根據(jù)其在平面形狀中的展向、弦向位置，應(yīng)用本文的參數(shù)化方法可以算出其在新機(jī)翼外形中的坐標(biāo)。這樣就可以把Baseline機(jī)翼表面網(wǎng)格“投影”到新機(jī)翼表面上。接下來(lái)以指數(shù)衰減規(guī)律將表面網(wǎng)格的位移，傳遞到每個(gè)網(wǎng)格塊的角點(diǎn)。最后用無(wú)限插值(TFI)方法插值得到內(nèi)部網(wǎng)格點(diǎn)新的空間坐標(biāo)。與此同時(shí)，機(jī)翼厚度、容積等幾何信息也都計(jì)算出來(lái)，這些量用來(lái)判斷是否滿足幾何約束條件。

圖2 計(jì)算網(wǎng)格Fig.2 Computational mesh

在變形過(guò)程中，網(wǎng)格拓?fù)浣Y(jié)構(gòu)不變，空間網(wǎng)格點(diǎn)的分布規(guī)律也變化不大。實(shí)際應(yīng)用中，該網(wǎng)格變形方法對(duì)機(jī)翼扭轉(zhuǎn)分布和剖面形狀的改變具有足夠的魯棒性，最重要的是，整個(gè)變形過(guò)程可在10 s左右完成。此外，表面和空間網(wǎng)格的變形可以分開(kāi)進(jìn)行，變形后的表面網(wǎng)格只相當(dāng)于空間網(wǎng)格變形的“邊界條件”。在優(yōu)化過(guò)程中，只需向“天河2號(hào)”傳輸表面網(wǎng)格，文件大小不到空間網(wǎng)格的1%。這樣可以顯著降低存儲(chǔ)系統(tǒng)讀寫(xiě)操作強(qiáng)度，利于提高并行計(jì)算規(guī)模。

1.3 CFD計(jì)算分析和自動(dòng)后處理

使用內(nèi)部CFD程序SFlow對(duì)外形方案計(jì)算分析，該CFD程序在多塊點(diǎn)對(duì)點(diǎn)結(jié)構(gòu)化網(wǎng)格上，以有限體積法求解雷諾平均N-S方程。無(wú)黏通量項(xiàng)用Roe平均迎風(fēng)通量差分分裂格式(FDS)離散，黏性通量項(xiàng)用中心差分格式離散，時(shí)間推進(jìn)計(jì)算采用隱式LU-SGS方法。該程序有SA一方程和SST兩方程湍流模型可選用。本文的CFD計(jì)算用SST兩方程湍流模型。采用MPI方式并行和基于貪婪算法的并行任務(wù)分配策略具有相當(dāng)高的并行效率。該CFD程序具有定升力計(jì)算功能，在迭代計(jì)算過(guò)程中，通過(guò)不斷調(diào)整來(lái)流迎角，逼近設(shè)定的升力系數(shù)值。

網(wǎng)格分辨率和迭代計(jì)算收斂標(biāo)準(zhǔn)嚴(yán)重影響氣動(dòng)力計(jì)算時(shí)間和可信度。使用多大網(wǎng)格量，如何設(shè)定收斂標(biāo)準(zhǔn)，是優(yōu)化設(shè)計(jì)前必須解決的問(wèn)題。圖3為針對(duì)WBNP構(gòu)型，網(wǎng)格單元數(shù)分別為200萬(wàn)、800萬(wàn)、1600萬(wàn)和6400萬(wàn)的4套網(wǎng)格在馬赫數(shù)0.85、升力系數(shù)0.5下的網(wǎng)格收斂性曲線。單元數(shù)1600萬(wàn)的網(wǎng)格是作者在氣動(dòng)性能計(jì)算工作中經(jīng)常使用的，在這里與6400萬(wàn)網(wǎng)格計(jì)算的阻力相差不到3 counts。單元數(shù)200萬(wàn)的網(wǎng)格由1600萬(wàn)的網(wǎng)格經(jīng)過(guò)半粗化得到。單元數(shù)800萬(wàn)網(wǎng)格是在1600萬(wàn)的網(wǎng)格基礎(chǔ)上，適當(dāng)降低機(jī)身、短艙/吊掛以及遠(yuǎn)場(chǎng)網(wǎng)格分辨率后得到，機(jī)翼附近網(wǎng)格分辨率基本不變。

圖3 網(wǎng)格收斂性曲線Fig.3 Grid convergence study

圖4為使用單元數(shù)800萬(wàn)的網(wǎng)格，以收斂好的Baseline流場(chǎng)解為初場(chǎng)，擾動(dòng)后外形CFD計(jì)算的收斂歷史。在“天河2號(hào)”上使用48核并行，同時(shí)使用多重網(wǎng)格和當(dāng)?shù)貢r(shí)間步長(zhǎng)等加速收斂措施，13 min左右即可完成600次迭代，阻力系數(shù)和力矩系數(shù)的收斂程度可以滿足要求。經(jīng)過(guò)計(jì)算精度和時(shí)間的權(quán)衡，下文的優(yōu)化設(shè)計(jì)中使用單元數(shù)800萬(wàn)網(wǎng)格，而驗(yàn)證計(jì)算則使用1600萬(wàn)的網(wǎng)格。

胡冬林（1955-2017），長(zhǎng)春人，吉林省作家協(xié)會(huì)副主席。出版有《原始森林手記》《鷹屯烏拉田野札記》《狐貍的微笑》、長(zhǎng)篇?jiǎng)游镄≌f(shuō)《野豬王》、兒童科幻長(zhǎng)篇小說(shuō)《巨蟲(chóng)公園》以及單篇散文《青羊消息》《拍濺》《金角鹿》等，是我國(guó)生態(tài)文學(xué)寫(xiě)作的領(lǐng)先作家。

圖4 力系數(shù)收斂歷史Fig.4 Force coefficient convergence history

CFD計(jì)算結(jié)束后，為優(yōu)化設(shè)計(jì)而發(fā)展的后處理程序從流場(chǎng)計(jì)算結(jié)果中自動(dòng)提取流場(chǎng)特征量和氣動(dòng)力系數(shù)。這里主要實(shí)現(xiàn)了指定截面壓力分布特征量(吸力峰、激波位置等)和展向載荷分布、升力系數(shù)分布(圖5)的提取分析。這些特征量和力系數(shù)都可作為優(yōu)化設(shè)計(jì)的目標(biāo)或約束。

圖5 CFD自動(dòng)后處理Fig.5 Automatic post-processing of CFD results

1.4 超級(jí)計(jì)算機(jī)上的優(yōu)化設(shè)計(jì)流程

本文的優(yōu)化設(shè)計(jì)系統(tǒng)以Windows/Linux跨平臺(tái)的方式運(yùn)行于“天河2號(hào)”超級(jí)計(jì)算機(jī)(下文簡(jiǎn)稱(chēng)“超算”)上。具體來(lái)說(shuō)是以操作系統(tǒng)為Windows Server的胖節(jié)點(diǎn)(32核心，64線程，256GB內(nèi)存)為整個(gè)優(yōu)化設(shè)計(jì)系統(tǒng)的控制臺(tái)，完成計(jì)算任務(wù)預(yù)處理，連接登錄節(jié)點(diǎn)提交計(jì)算任務(wù)，以及優(yōu)化搜索計(jì)算。登陸節(jié)點(diǎn)和計(jì)算節(jié)點(diǎn)操作系統(tǒng)為L(zhǎng)inux，分別負(fù)責(zé)計(jì)算任務(wù)排隊(duì)提交與計(jì)算結(jié)果回傳、CFD計(jì)算與后處理。經(jīng)過(guò)反復(fù)測(cè)試，這樣的系統(tǒng)架構(gòu)運(yùn)行穩(wěn)定，超算系統(tǒng)本身承受的壓力較小。

具體優(yōu)化流程如圖6所示。每個(gè)新個(gè)體都是在Baseline外形的基礎(chǔ)上，疊加由尋優(yōu)算法給出的擾動(dòng)量得到，相應(yīng)的表面網(wǎng)格也會(huì)做相應(yīng)的擾動(dòng)變形。與此同時(shí)，新外形的幾何特征也會(huì)被分析出來(lái)，并判斷是否滿足幾何約束。如果不滿足幾何約束，該個(gè)體將被淘汰，不再進(jìn)行后續(xù)的CFD計(jì)算分析。這樣做的好處在于盡早剔除不滿足約束的個(gè)體，避免無(wú)效的CFD計(jì)算并加快優(yōu)化進(jìn)程。對(duì)于滿足幾何約束的個(gè)體，其表面網(wǎng)格將被傳輸給“天河二號(hào)”，通過(guò)登錄節(jié)點(diǎn)提交計(jì)算作業(yè)，在計(jì)算節(jié)點(diǎn)上依次完成空間網(wǎng)格變形、CFD計(jì)算和自動(dòng)后處理過(guò)程。最后，計(jì)算結(jié)果傳回胖節(jié)點(diǎn)，判斷是否滿足氣動(dòng)約束，根據(jù)目標(biāo)函數(shù)值判定個(gè)體優(yōu)劣。完成一代種群的分析后，優(yōu)化算法會(huì)生成新一代個(gè)體。本文的優(yōu)化案例采用具備精英策略的非支配排序遺傳算法 NSGA-II實(shí)現(xiàn)多目標(biāo)尋優(yōu)。

圖6 優(yōu)化設(shè)計(jì)流程Fig.6 Flow chart of optimization

本文開(kāi)展的基于大規(guī)模并行計(jì)算的民機(jī)復(fù)雜外形優(yōu)化設(shè)計(jì)在超算應(yīng)用領(lǐng)域?qū)儆诟咄款?lèi)型計(jì)算，即單個(gè)計(jì)算作業(yè)占用資源較少(本文中使用48核持續(xù)不到15 min)，但同時(shí)進(jìn)行的計(jì)算作業(yè)數(shù)量龐大(數(shù)百甚至數(shù)千個(gè))。每個(gè)計(jì)算作業(yè)的提交和監(jiān)控都必須通過(guò)登錄節(jié)點(diǎn)完成，每個(gè)作業(yè)都需要拷貝文件、新建文件和讀寫(xiě)文件。高通量應(yīng)用最為考驗(yàn)超級(jí)計(jì)算機(jī)的綜合性能，對(duì)計(jì)算資源、I/O系統(tǒng)、網(wǎng)絡(luò)系統(tǒng)、文件存儲(chǔ)和文件管理等系統(tǒng)資源消耗巨大。超算的登錄節(jié)點(diǎn)是一個(gè)公共平臺(tái)，主要用于管理員對(duì)超算環(huán)境維護(hù)管理，及所有用戶提交計(jì)算作業(yè)和上傳下載數(shù)據(jù)，有時(shí)還作為各種自編程序的編譯環(huán)境。因此登錄節(jié)點(diǎn)資源不允許單個(gè)用戶長(zhǎng)時(shí)間或大量占用。

研究發(fā)現(xiàn)，高通量計(jì)算應(yīng)用的主要瓶頸在于登錄節(jié)點(diǎn)的線程處理能力和整個(gè)超算平臺(tái)的I/O系統(tǒng)。若要提高計(jì)算規(guī)模，必須設(shè)法降低單個(gè)計(jì)算作業(yè)消耗的系統(tǒng)資源及通訊量。采用Python語(yǔ)言編寫(xiě)腳本，在胖節(jié)點(diǎn)上實(shí)現(xiàn)遠(yuǎn)程ssh連接、計(jì)算作業(yè)提交及相應(yīng)文件上傳。計(jì)算任務(wù)提交完畢之后，關(guān)閉對(duì)應(yīng)ssh連接。計(jì)算完成后，再由Python腳本調(diào)用登錄節(jié)點(diǎn)本地的shell腳本將計(jì)算結(jié)果下載回胖節(jié)點(diǎn)。在緩解I/O系統(tǒng)壓力方面，通過(guò)代碼優(yōu)化、流程改進(jìn)等措施盡可能減少文件創(chuàng)建以及讀寫(xiě)操作，采用這種“登錄節(jié)點(diǎn)和I/O系統(tǒng)減負(fù)”方案后，登錄節(jié)點(diǎn)CPU占用率和內(nèi)存消耗明顯減少，優(yōu)化系統(tǒng)運(yùn)行更加穩(wěn)定。

2 優(yōu)化案例

2.1 CRM安裝效應(yīng)分析

NASA CRM是以巡航馬赫數(shù)0.85的寬體客機(jī)為背景設(shè)計(jì)的標(biāo)準(zhǔn)模型，設(shè)計(jì)過(guò)程是先在翼身組合體構(gòu)型(下文簡(jiǎn)稱(chēng)：WB構(gòu)型)中設(shè)計(jì)一副性能優(yōu)良并且足夠魯棒的機(jī)翼，以保證集成短艙、吊掛后性能損失不大。為了降低網(wǎng)格生成難度，動(dòng)力裝置用一個(gè)單圈通氣短艙模擬，其流量系數(shù)與典型商用飛機(jī)巡航時(shí)相當(dāng)。WB構(gòu)型巡航設(shè)計(jì)點(diǎn)為馬赫數(shù)Ma=0.85，基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)Re=4×107，升力系數(shù)CL=0.5。此工況下WB構(gòu)型和相同迎角(1.875°)下WBNP構(gòu)型表面等壓線、機(jī)翼剖面壓力分布的CFD計(jì)算結(jié)果對(duì)比分別如圖7、圖8所示。圖8中前兩處剖面位于吊掛內(nèi)側(cè)，后兩處位于外側(cè)。集成短艙吊掛后，內(nèi)翼上表面在附加上洗作用下吸力峰抬高，激波前移增強(qiáng)，內(nèi)翼下表面在機(jī)翼、吊掛和短艙構(gòu)成的通道內(nèi)顯著加速。吊掛外側(cè)機(jī)翼上表面在下洗作用下吸力峰降低。根據(jù)設(shè)計(jì)經(jīng)驗(yàn)，內(nèi)翼激波增強(qiáng)會(huì)對(duì)阻力發(fā)散特性不利。

(a) WB構(gòu)型 (b) WBNP構(gòu)型

圖8 短艙吊掛安裝前后機(jī)翼剖面壓力分布對(duì)比Fig.8 Wing sectional pressure distribution before and after engine installation

圖9對(duì)比了短艙、吊掛集成前后機(jī)翼展向載荷分布，圖中橫坐標(biāo)Eta為機(jī)翼展向相對(duì)位置，縱坐標(biāo)Load為當(dāng)?shù)厣ο禂?shù)與弦長(zhǎng)的乘積。在迎角不變的情況下，WBNP構(gòu)型機(jī)翼升力顯著損失的范圍超過(guò)60%翼展。而且在吊掛的阻斷作用下，機(jī)翼載荷分布在吊掛位置存在間斷，這是在飛發(fā)集成構(gòu)型設(shè)計(jì)中需要考慮吊掛的一個(gè)重要原因。CRM的設(shè)計(jì)初衷除了用于驗(yàn)證數(shù)值計(jì)算工具和標(biāo)定風(fēng)洞外，還將用作考察優(yōu)化設(shè)計(jì)工具的標(biāo)準(zhǔn)模型。WB構(gòu)型的機(jī)翼載荷分布被設(shè)計(jì)成十分接近橢圓，因此就單獨(dú)WB構(gòu)型來(lái)說(shuō)，很難在機(jī)翼厚度不減的情況下同時(shí)在多個(gè)狀態(tài)點(diǎn)上取得顯著減阻效果。而集成了短艙吊掛后的CRM在激波阻力和誘導(dǎo)阻力方面都應(yīng)該有可降低的潛力。

圖9 短艙吊掛安裝前后機(jī)翼展向載荷分布對(duì)比Fig.9 Wing load distribution before and after engine installation

2.2 機(jī)翼彎扭分布優(yōu)化設(shè)計(jì)

平面形狀確定后，機(jī)翼高速氣動(dòng)性能在很大程度上取決于厚度分布和彎扭分布的組合。厚度主要影響激波阻力，在內(nèi)部容積、結(jié)構(gòu)布置等條件的約束下，厚度分布往往沒(méi)有太多縮減的余地。通過(guò)改變機(jī)翼彎扭分布可以調(diào)整展向升力分布和載荷分布，進(jìn)而影響誘導(dǎo)阻力、俯仰力矩以及低速失速特性。使用本文的優(yōu)化系統(tǒng)，在厚度不減的前提下對(duì)CRM WBNP構(gòu)型的機(jī)翼扭轉(zhuǎn)分布和控制剖面中弧線形狀進(jìn)行優(yōu)化設(shè)計(jì)。優(yōu)化問(wèn)題定義如下：

minCD_Ma=0.83&CD_Ma=0.85&CD_Ma=0.87

CL=0.5

s.t.Cm_Ma=0.85≥-0.145

α_Ma=0.85≤2.2°

(7)

本次實(shí)施的是3點(diǎn)3目標(biāo)優(yōu)化，3個(gè)設(shè)計(jì)點(diǎn)分別為馬赫數(shù)0.83、0.85和0.87，每個(gè)設(shè)計(jì)點(diǎn)的阻力系數(shù)作為獨(dú)立的目標(biāo)函數(shù)。設(shè)計(jì)變量包括9個(gè)控制剖面的扭轉(zhuǎn)角和表達(dá)中弧線的8階CST參數(shù)，共90個(gè)。CST參數(shù)的變化范圍設(shè)定為在Baseline值基礎(chǔ)上擾動(dòng)±0.05。翼根剖面的扭轉(zhuǎn)角變化范圍±0.5°，為保證扭轉(zhuǎn)曲線單調(diào)下降，其余控制剖面的扭轉(zhuǎn)角定義為在內(nèi)側(cè)相鄰剖面扭角絕對(duì)值基礎(chǔ)上的下降量，變化范圍設(shè)定為在Baseline基礎(chǔ)上擾動(dòng)±0.5°。這些量值的設(shè)定依賴(lài)于大量應(yīng)用后的經(jīng)驗(yàn)總結(jié)，而且在下文的優(yōu)化過(guò)程中，沒(méi)有發(fā)現(xiàn)變量向邊界收斂的現(xiàn)象。

設(shè)計(jì)約束包括：升力系數(shù)CL=0.5；設(shè)計(jì)點(diǎn)Ma=0.85的俯仰力矩系數(shù)和迎角。在CFD程序定升力計(jì)算功能的幫助下，可直接實(shí)現(xiàn)升力系數(shù)的等式約束。計(jì)算評(píng)估后，不滿足力矩系數(shù)和迎角約束的個(gè)體將會(huì)被列為不可行方案，在排序中不可行方案的等級(jí)低于所有可行方案。隨著種群逐漸進(jìn)化成熟，這種不可行個(gè)體所占比例也會(huì)越來(lái)越少。

本次優(yōu)化種群規(guī)模為256，每批次同時(shí)對(duì)128個(gè)個(gè)體進(jìn)行CFD計(jì)算，每個(gè)個(gè)體串行完成3個(gè)設(shè)計(jì)點(diǎn)的計(jì)算。在60 h內(nèi)，完成超過(guò)10 000個(gè)個(gè)體的計(jì)算分析，遺傳優(yōu)化40代。圖10為優(yōu)化過(guò)程中3個(gè)目標(biāo)函數(shù)的演化歷史，從下沿輪廓來(lái)看，優(yōu)化已經(jīng)趨于收斂。

圖10 優(yōu)化過(guò)程中目標(biāo)函數(shù)演化歷史Fig.10 History of cost functions in optimization

多目標(biāo)優(yōu)化問(wèn)題以找到Pareto前緣作為結(jié)束。最終在多個(gè)目標(biāo)間權(quán)衡選擇還需要由設(shè)計(jì)者根據(jù)應(yīng)用需要完成。分析優(yōu)化結(jié)果發(fā)現(xiàn)，馬赫數(shù)0.83和0.85兩個(gè)目標(biāo)間基本不構(gòu)成競(jìng)爭(zhēng)關(guān)系。因此最優(yōu)個(gè)體的遴選主要在馬赫數(shù)0.85和0.87兩個(gè)目標(biāo)間的Pareto前緣上進(jìn)行。圖11展示了優(yōu)化過(guò)程中成功完成計(jì)算的個(gè)體的馬赫數(shù)0.85和0.87兩個(gè)目標(biāo)函數(shù)值分布，藍(lán)色方形標(biāo)志為Baseline方案所在位置，在Pareto前緣上選定三個(gè)紅色菱形標(biāo)志的個(gè)體為此次優(yōu)化結(jié)果的備選方案。圖12對(duì)比了他們和Baseline方案的表面等壓線圖，三個(gè)備選方案的內(nèi)翼激波強(qiáng)度都明顯降低。

表1列出了他們與Baseline方案的阻力系數(shù)對(duì)比。ID號(hào)9740的個(gè)體在馬赫數(shù)0.85時(shí)阻力最低，比Baseline方案降低3.1 counts，但在馬赫數(shù)0.87時(shí)阻力僅降低3.3 counts。9415號(hào)個(gè)體在馬赫數(shù)0.87時(shí)阻力下降最多，達(dá)到14.2 counts，但在馬赫數(shù)0.85時(shí)阻力僅降低0.6 counts。7777號(hào)個(gè)體在馬赫數(shù)0.85時(shí)阻力降低2.3 counts，在馬赫數(shù)0.87時(shí)阻力降低10.4 counts，各設(shè)計(jì)點(diǎn)阻力下降比較均衡，下面重點(diǎn)分析該個(gè)體的細(xì)節(jié)特征。圖13對(duì)比了優(yōu)化前后機(jī)翼扭轉(zhuǎn)分布的差異，除了翼梢扭轉(zhuǎn)角有所增加，大部分區(qū)域扭轉(zhuǎn)角基本保持不變。圖14則對(duì)比了優(yōu)化前后9個(gè)控制剖面的形狀變化。

圖12 優(yōu)化前后表面等壓線圖對(duì)比Fig.12 Cp contours before and after optimization

IDCD_Ma=0.83/countCD_Ma=0.85/countCD_Ma=0.87/countBaseline243.8244.7265.87777241.9242.4255.49415242.3244.1251.69740240.0241.6262.5

圖13 機(jī)翼扭轉(zhuǎn)分布對(duì)比Fig.13 Comparison of wing twist distribution

圖14 控制剖面形狀對(duì)比Fig.14 Comparison of wing sectional shapes

圖15和圖16分別對(duì)比了馬赫數(shù)0.85和0.87優(yōu)化前后四處剖面壓力分布的變化。優(yōu)化后，在馬赫數(shù)0.85時(shí)，內(nèi)翼部分激波強(qiáng)度明顯減弱，吊掛外側(cè)則表現(xiàn)為吸力峰提高。馬赫數(shù)0.87時(shí)，各剖面激波后移減弱。

圖15 優(yōu)化前后馬赫數(shù)0.85機(jī)翼剖面壓力分布對(duì)比Fig.15 Comparison of wing sectional pressure distribution at Ma=0.85

圖16 優(yōu)化前后馬赫數(shù)0.87機(jī)翼剖面壓力分布對(duì)比Fig.16 Comparison of wing sectional pressure distribution at Ma=0.87

2.3 優(yōu)化結(jié)果驗(yàn)證

本文的機(jī)翼三維參數(shù)化方法與CATIA中的放樣成型不可能完全相同，同時(shí)優(yōu)化過(guò)程中網(wǎng)格分辨率和收斂判據(jù)也有所降低，因此優(yōu)化取得阻力下降的效果還需要進(jìn)一步驗(yàn)證。將優(yōu)化后的控制剖面在CATIA軟件中放樣成型，得到優(yōu)化后機(jī)翼的CAD模型。在ICEM-CFD軟件中重新生成單元數(shù)1600萬(wàn)的計(jì)算網(wǎng)格開(kāi)展驗(yàn)證計(jì)算。計(jì)算從自由來(lái)流初場(chǎng)開(kāi)始迭代，殘值下降4個(gè)量級(jí)后提取計(jì)算結(jié)果。圖17為優(yōu)化前后阻力發(fā)散曲線的對(duì)比，可見(jiàn)優(yōu)化過(guò)程中所取得的減阻成果基本得到驗(yàn)證。

圖17 阻力發(fā)散特性優(yōu)化效果驗(yàn)證Fig.17 Validation of drag reduction

3 結(jié) 論

探索了應(yīng)用超級(jí)計(jì)算機(jī)開(kāi)展民機(jī)飛發(fā)集成構(gòu)型精細(xì)優(yōu)化設(shè)計(jì)。直接使用RANS計(jì)算評(píng)估，同時(shí)在近距耦合部件間的干擾作用下開(kāi)展全局尋優(yōu)的多目標(biāo)優(yōu)化設(shè)計(jì)，得到的減阻效果可靠。主要工作簡(jiǎn)述如下：

1) 搭建了可用于民機(jī)WBNP構(gòu)型的優(yōu)化設(shè)計(jì)系統(tǒng)，采用遺傳算法進(jìn)行全局尋優(yōu)，包括CFD前處理、RANS計(jì)算的CFD分析和自動(dòng)后處理等功能模塊。該系統(tǒng)在超級(jí)計(jì)算機(jī)上可同時(shí)對(duì)數(shù)百個(gè)設(shè)計(jì)案例開(kāi)展計(jì)算分析。并行計(jì)算、加速收斂和適當(dāng)降低收斂標(biāo)準(zhǔn)等技術(shù)的綜合運(yùn)用使優(yōu)化過(guò)程能在可接受的時(shí)間周期內(nèi)完成。

2) 在厚度不變的情況下，對(duì)CRM WBNP構(gòu)型中的機(jī)翼彎扭分布開(kāi)展了多目標(biāo)優(yōu)化設(shè)計(jì)。60 h內(nèi)完成了40代遺傳優(yōu)化，各設(shè)計(jì)點(diǎn)阻力系數(shù)有2～10 counts的下降，阻力發(fā)散特性顯著改善。驗(yàn)證計(jì)算結(jié)果表明，優(yōu)化過(guò)程中取得的阻力下降量是可靠的。

致謝:感謝廣州超算中心對(duì)本文工作的大力支持。