薛幫猛, 張文升, 張志雄

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

0 引 言

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

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

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

1 優(yōu)化方法

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

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

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

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

(1)

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

ψ=(x-xLE)/c

(2)

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

ζ=(y-yLE)/c

(3)

后緣開口相對大小：

ζT=ΔyTE/c

(4)

類型函數(shù)：

(5)

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

(6)

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

這里將機(jī)翼表面視為由沿展向的無限多個剖面翼型構(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ù)叵议Lc、后緣開度ζT則由附近兩個控制剖面線性插值得到。這就是本文提出的機(jī)翼三維參數(shù)化方法。給定機(jī)翼平面形狀內(nèi)任一位置(x,z)，可以唯一地算出上、下表面y坐標(biāo)yU、yL。

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

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

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

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

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

1.3 CFD計算分析和自動后處理

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

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

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

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

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

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

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

圖5 CFD自動后處理Fig.5 Automatic post-processing of CFD results

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

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

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

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

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

研究發(fā)現(xiàn)，高通量計算應(yīng)用的主要瓶頸在于登錄節(jié)點(diǎn)的線程處理能力和整個超算平臺的I/O系統(tǒng)。若要提高計算規(guī)模，必須設(shè)法降低單個計算作業(yè)消耗的系統(tǒng)資源及通訊量。采用Python語言編寫腳本，在胖節(jié)點(diǎn)上實(shí)現(xiàn)遠(yuǎn)程ssh連接、計算作業(yè)提交及相應(yīng)文件上傳。計算任務(wù)提交完畢之后，關(guān)閉對應(yīng)ssh連接。計算完成后，再由Python腳本調(diào)用登錄節(jié)點(diǎn)本地的shell腳本將計算結(jié)果下載回胖節(jié)點(diǎn)。在緩解I/O系統(tǒng)壓力方面，通過代碼優(yōu)化、流程改進(jìn)等措施盡可能減少文件創(chuàng)建以及讀寫操作，采用這種“登錄節(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è)計的標(biāo)準(zhǔn)模型，設(shè)計過程是先在翼身組合體構(gòu)型(下文簡稱：WB構(gòu)型)中設(shè)計一副性能優(yōu)良并且足夠魯棒的機(jī)翼，以保證集成短艙、吊掛后性能損失不大。為了降低網(wǎng)格生成難度，動力裝置用一個單圈通氣短艙模擬，其流量系數(shù)與典型商用飛機(jī)巡航時相當(dāng)。WB構(gòu)型巡航設(shè)計點(diǎn)為馬赫數(shù)Ma=0.85，基于平均氣動弦長的雷諾數(shù)Re=4×107，升力系數(shù)CL=0.5。此工況下WB構(gòu)型和相同迎角(1.875°)下WBNP構(gòu)型表面等壓線、機(jī)翼剖面壓力分布的CFD計算結(jié)果對比分別如圖7、圖8所示。圖8中前兩處剖面位于吊掛內(nèi)側(cè)，后兩處位于外側(cè)。集成短艙吊掛后，內(nèi)翼上表面在附加上洗作用下吸力峰抬高，激波前移增強(qiáng)，內(nèi)翼下表面在機(jī)翼、吊掛和短艙構(gòu)成的通道內(nèi)顯著加速。吊掛外側(cè)機(jī)翼上表面在下洗作用下吸力峰降低。根據(jù)設(shè)計經(jīng)驗(yàn)，內(nèi)翼激波增強(qiáng)會對阻力發(fā)散特性不利。

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

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

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

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

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

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

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

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

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

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

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

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

圖12 優(yōu)化前后表面等壓線圖對比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)分布對比Fig.13 Comparison of wing twist distribution

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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