薛幫猛，張文升，孫學(xué)衛(wèi)，吳宇昂

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

大型民用飛機機翼在產(chǎn)生升力的同時，阻力貢獻也占到了全機一半，機翼的氣動性能嚴重影響飛機的經(jīng)濟性。常規(guī)翼吊氣動布局中，機翼、機身、發(fā)動機短艙、吊掛之間相互干擾，形成一個緊密聯(lián)系的部件系統(tǒng)(簡稱WBPN構(gòu)型)。20世紀90年代德國宇航院的研究顯示[1]：吊裝涵道比為8～9的渦扇發(fā)動機會使機翼超過60%翼展范圍的升力顯著損失，而吊掛也將造成內(nèi)、外翼展向升力系數(shù)分布的不連續(xù)。文獻[2]認為短艙/吊掛對機翼的干擾源于其誘導(dǎo)了額外的上、下洗作用?？梢悦鞔_的是，真實飛機的機翼需要在比均勻來流復(fù)雜得多的條件下工作。與單獨設(shè)計機翼或翼身組合體外形相比，在更接近全機外形的部件系統(tǒng)中優(yōu)化設(shè)計機翼，對提高最終整機的氣動性能更有意義。巴西航空工業(yè)公司Oliveira等[3]也指出了在飛機/發(fā)動機集成構(gòu)型中優(yōu)化設(shè)計部件外形的必要性。東京大學(xué)Nakahashi團隊[4-7]使用伴隨梯度優(yōu)化方法進行了一系列針對WBPN構(gòu)型的多點優(yōu)化設(shè)計研究。最近，密歇根大學(xué)Martins團隊[8]用梯度優(yōu)化算法實現(xiàn)了支線機全機構(gòu)型中機翼外形的優(yōu)化設(shè)計。中國在飛/發(fā)集成構(gòu)型的優(yōu)化設(shè)計方面也做出了很多卓有成效的研究工作，左英桃等[9]使用伴隨梯度優(yōu)化方法在WBPN構(gòu)型下，同時對DLR-F6的機翼和短艙實施了減阻優(yōu)化設(shè)計。張宇飛等[10]嘗試了用遺傳優(yōu)化算法在短艙干擾下實施機翼一體化設(shè)計。筆者所在團隊[11]利用超級計算機和大規(guī)模并行計算，實現(xiàn)了在可接受的時間周期內(nèi)完成NASA CRM(Common Research Model)飛/發(fā)集成構(gòu)型下的機翼多目標優(yōu)化過程。

除了機身、短艙和吊掛等有形部件的干擾，發(fā)動機噴流對機翼的影響同樣不容忽視。模擬發(fā)動機進排氣的流場中，在噴流的誘導(dǎo)作用下，其附近流場的速度大小和方向都將發(fā)生改變。翼吊氣動布局中，發(fā)動機噴流不但增加噴流區(qū)內(nèi)掛架上的刮擦阻力，對機翼下表面流動的加速作用會直接使升力損失，并間接增加巡航阻力。相比于窄體機，吊裝更大涵道比和尺寸發(fā)動機的遠程寬體客機，發(fā)動機噴流對機翼的氣動干擾問題更為嚴峻。文獻[12-16]在渦扇發(fā)動機噴流對機翼的干擾效應(yīng)和氣動性能影響方面做了大量研究。風(fēng)洞試驗中噴流干擾的研究手段主要是使用渦輪風(fēng)扇模擬器(Turbofan Powered Simulators，TPS)[14-19]，但這種試驗研究非常昂貴，不適合大量開展。另外，使用數(shù)值計算手段模擬發(fā)動機進排氣時物面不封閉，直接由物面積分得到的阻力并不是真正意義上的機體阻力。此時，要將發(fā)動機推力和機體阻力劃分清楚，需要經(jīng)過較為復(fù)雜的推力/阻力分離計算[14-16]，這對噴流干擾下的機體減阻設(shè)計帶來了挑戰(zhàn)。

本文嘗試直接在發(fā)動機噴流作用下開展寬體客機機翼氣動外形的多目標減阻優(yōu)化設(shè)計。首先，用動量定理分析動力短艙流場，選擇合適的變量作為優(yōu)化設(shè)計目標。接下來，以某寬體客機為對象，研究巡航工況下噴流對機翼氣動干擾的程度，明確在噴流作用下設(shè)計機翼的必要性。最后使用前期建立的基于雷諾平均Navier-Stokes(RANS)方程計算分析和遺傳算法尋優(yōu)的多目標優(yōu)化系統(tǒng)[11]，在CFD(Computational Fluid Dynamics)計算中以進排氣邊界條件模擬動力短艙，實施機翼外形的多目標優(yōu)化設(shè)計。

1 動力外形流場分析和噴流干擾計算

CFD計算圖1所示的動力短艙[14]繞流時，進排氣模擬是通過邊界條件來實現(xiàn)的，風(fēng)扇入口處(E12)給定靜壓，風(fēng)扇出口(E17)和核心機出口(E7)則給定總壓比和總溫比。圖中VPre、VPost分別為上、下游遠場速度大??；DPre、DPost分別為進、排氣流管內(nèi)物面積分阻力大?。沪誔re、φPost分別為進、排氣流管表面積分阻力大小；DFrame為機體表面積分阻力大小。

圖1 動力短艙流場控制體劃分[14]Fig.1 Control volume division of powered-on nacelle flow-field[14]

對進、排氣流管分別列流動方向動量方程，推導(dǎo)后可得推力區(qū)凈推力FN和阻力區(qū)總阻力DExt的表達式，推力與阻力之差為

FN-DExt=(F17+F7-F12)-

(DFrame+DPre+DPost)

(1)

式中：F17+F7-F12為排氣與進氣邊界上“毛推力”之差；DFrame+DPre+DPost則為整個流場中所有壁面積分流向力的大小，也即CFD計算直接給出的“阻力”，如前述因壁面不封閉它并不是真正的阻力。給定進排氣邊界和來流條件，僅修改機體外形時F17+F7-F12大小不會變。在機體部件減阻優(yōu)化設(shè)計中，降低的無論是阻力還是推力損失，都會使推力與阻力之差增大，這正是期望的優(yōu)化目標方向。此時FN-DExt和DFrame+DPre+DPost的變化量大小相等，因此以壁面積分流向力最小化為目標也是合理的。

本文中CFD計算使用內(nèi)部CFD程序SFlow[11]，該程序在多塊點對點結(jié)構(gòu)化網(wǎng)格上以有限體積法求解雷諾平均Navier-Stokes方程。無黏通量項用Roe平均迎風(fēng)通量差分分裂格式離散，黏性通量項用中心差分格式離散，時間推進計算采用隱式LU-SGS(Lower Upper-Symmetric Gauss Seidel)方法。SFlow有S-A(Spalart Allmaras)一方程和SST(Shear Stress Transport)兩方程湍流模型，本文的CFD計算采用SST兩方程湍流模型。MPI(Message Passing Interface)方式并行和基于貪婪算法的并行任務(wù)分配策略使SFlow有很高的并行效率。SFlow具有定升力計算功能，在迭代計算過程中，通過不斷調(diào)整來流攻角，逼近設(shè)定的升力系數(shù)值。SFlow程序經(jīng)過大量標準模型的驗證計算，其計算精度也在很多工程項目中確認。圖2展示了NASA CRM機翼-機身-水平尾翼構(gòu)型的計算結(jié)果(圖中Cp為壓力系數(shù)，x/c為無量綱弦向位置)，機翼壓力分布與NASA網(wǎng)站提供的風(fēng)洞試驗結(jié)果一致性良好。計算狀態(tài)為：馬赫數(shù)Ma=0.85，升力系數(shù)CL=0.534， 基于平均氣動弦長的雷諾數(shù)Re=3×107。計算中使用了包含1 200萬單元的網(wǎng)格，并且用網(wǎng)格變形方法使機翼計算外形的扭轉(zhuǎn)變形與風(fēng)洞試驗中測得的變形量一致。

是否有必要直接在噴流作用下設(shè)計機翼，取決于噴流干擾程度。圖3給出了某寬體客機的翼身組合體、機翼/機身/吊掛/通氣短艙、機翼/機身/吊掛/動力短艙(分別簡稱WB、WBP+TFN、WBP+PN構(gòu)型)外形，在馬赫數(shù)Ma=0.85，迎角α=1.95°和基于平均氣動弦長雷諾數(shù)Re=5×107狀態(tài)下表面等壓線的CFD計算結(jié)果。通氣短艙(TFN)和動力短艙(PN)的外形面形狀和流量系數(shù)相同。計算網(wǎng)格方面，3個外形的網(wǎng)格單元數(shù)都在1 000萬左右(表1)，物面第1層尺度足以保證y+在1左右。WBP+TFN和WBP+PN的網(wǎng)格在機翼和機身周圍的拓撲結(jié)構(gòu)和點分布規(guī)律是一致的，該對比關(guān)注升力和壓力分布，不做阻力的對比。

圖2 NASA CRM機翼截面壓力分布計算與試驗對比Fig.2 Comparison of NASA CRM wing sectional pressure distributions between CFD calculation and test

圖3 安裝和噴流干擾下的機翼表面等壓線圖Fig.3 Contours of wing surface pressure coefficient under installation and jet effects

由圖3可見，飛/發(fā)集成構(gòu)型內(nèi)翼為弱激波，去掉短艙吊掛后的WB構(gòu)型內(nèi)翼卻表現(xiàn)為位置靠后的強激波。圖4給出了內(nèi)翼26%半展長位置機翼壓力分布對比，在噴流作用下，內(nèi)翼激波位置前移，且強度略有增加。圖5對比了機翼展向升力系數(shù)分布(橫坐標η為展向相對位置)，可以看到從WB到WBP+TFN，再到WBP+PN，機翼升力是逐步降低的。從表1可以定量地看到，安裝通氣短艙和吊掛后，升力系數(shù)損失了0.018，加入噴流后升力系數(shù)又損失0.021。如果最終配平后全機巡航升力系數(shù)為0.5，那么這兩部分損失疊加高達全機升力系數(shù)的8%。因此，從干擾導(dǎo)致升力損失這個角度來說，翼吊布局寬體客機機翼氣動設(shè)計不僅應(yīng)該考慮短艙吊掛的安裝效應(yīng)，還應(yīng)該考慮噴流干擾。

表1 發(fā)動機安裝和噴流干擾下的升力系數(shù)Table 1 Lift coefficients under installation and jet effects of engine

圖4 安裝和噴流干擾下的機翼截面壓力系數(shù)分布對比Fig.4 Comparison of wing sectional pressure coefficient distributions under installation and jet effects

圖5 安裝和噴流干擾下的機翼展向升力系數(shù)分布對比Fig.5 Comparison of wing spanwise lift coefficient distributions under installation and jet effects

2 優(yōu)化設(shè)計方法

本節(jié)介紹優(yōu)化設(shè)計中涉及的幾何外形參數(shù)、計算網(wǎng)格生成、CFD計算和優(yōu)化流程等方面。

優(yōu)化設(shè)計的對象為WBP+PN構(gòu)型下的機翼外形，具體來說是有限個機翼控制剖面。使用CST(Class function/Shape function Transformation)參數(shù)化方法[20-22]表達三維機翼外形，設(shè)計變量可以是控制剖面的扭轉(zhuǎn)角或CST參數(shù)。該參數(shù)化方法具有一些設(shè)計者關(guān)心的直觀參數(shù)，可直接顯式控制翼型前緣半徑、尾錐角和后緣厚度等幾何特征。設(shè)計變量以擾動量的形式疊加到初始外形(Baseline外形)的基本參數(shù)上，基本參數(shù)需在優(yōu)化開始前對初始控制剖面經(jīng)過參數(shù)反算獲得。

優(yōu)化中使用網(wǎng)格變形技術(shù)實現(xiàn)新外形計算網(wǎng)格的快速生成。基于必要的網(wǎng)格收斂性研究，計算精度和時間的權(quán)衡，本文優(yōu)化設(shè)計中使用包含445個網(wǎng)格塊和1 065萬單元的結(jié)構(gòu)化網(wǎng)格。初始網(wǎng)格在ANSYS ICEM-CFD軟件中生成，物面第1層尺度保證y+≈1，法向增長率為1.2。 這里的網(wǎng)格變形是與機翼參數(shù)化方法集成在一起的。給定一組機翼參數(shù)后，對于Baseline機翼表面的每一個網(wǎng)格點，根據(jù)其在平面形狀中的展向、弦向位置，用參數(shù)化方法可以算出其在新機翼外形中的坐標。算出所有表面網(wǎng)格點的新坐標后也就生成了新的表面網(wǎng)格。接下來以指數(shù)衰減規(guī)律將表面網(wǎng)格角點的位移，傳遞到每個網(wǎng)格塊的角點。最后用無限插值(TFI)方法插值逐步得到內(nèi)部網(wǎng)格線、面、體新的空間坐標。同時，機翼厚度、容積等幾何信息也都計算出來，這些量用來判斷方案是否滿足幾何約束。

為了實現(xiàn)優(yōu)化過程中的快速CFD計算評估，計算將以收斂好的Baseline外形流場解為初場，迭代計算擾動后外形的流場，同時使用多重網(wǎng)格和當?shù)貢r間步長等加速收斂措施。經(jīng)大量測試證明，用96核并行，10 min左右即可完成600次迭代，阻力系數(shù)和力矩系數(shù)的收斂程度可以滿足要求。圖6為有代表性的新外形阻力系數(shù)(CD)計算收斂歷史(圖中1 count表示阻力系數(shù)為0.000 1)。

本文以遺傳算法為驅(qū)動，調(diào)動外形參數(shù)化、網(wǎng)格變形和CFD計算形成設(shè)計循環(huán)，實施多目標優(yōu)化設(shè)計。優(yōu)化流程如圖7所示,每個新個體都是在Baseline外形的基礎(chǔ)上，疊加由尋優(yōu)算法給出的擾動量得到，表面網(wǎng)格也會隨擾動相應(yīng)變形。同時，新外形的幾何特征也會被分析出來，并判斷是否滿足幾何約束。如果不滿足幾何約束，該個體將被淘汰，不再進行后續(xù)的CFD計算分析。對于滿足幾何約束的個體，由計算結(jié)果判斷是否滿足氣動約束，根據(jù)目標函數(shù)值判定個體優(yōu)劣。完成一代種群的分析后，優(yōu)化算法會生成新一代個體。本文的優(yōu)化案例采用具備精英策略的非支配排序遺傳算法NSGA-II[23]實現(xiàn)多目標尋優(yōu)。

圖6 阻力系數(shù)收斂歷史Fig.6 History of drag coefficient convergence

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

3 機翼優(yōu)化設(shè)計

Baseline外形為前面計算的寬體客機WBP+PN構(gòu)型，該方案在前期已經(jīng)過初步人工修形設(shè)計，巡航馬赫數(shù)為0.85時已有較為理想的升阻比和壓力分布形態(tài)，但阻力發(fā)散性能不佳。下面將在此基礎(chǔ)上優(yōu)化機翼的扭轉(zhuǎn)分布和剖面形狀，在滿足幾何和氣動約束下減小阻力。除了巡航馬赫數(shù)為0.85的狀態(tài)，還考慮馬赫數(shù)為0.83、0.87兩個狀態(tài)。對本次優(yōu)化設(shè)計的預(yù)期是：在保證馬赫數(shù)為0.85時阻力不增的前提下，盡可能降低馬赫數(shù)為0.87時的阻力，同時馬赫數(shù)為0.83的阻力不高于馬赫數(shù)為0.85的。

優(yōu)化問題定義為：進行3點3目標優(yōu)化，3個設(shè)計點分別為飛行馬赫數(shù)為0.83、0.85和0.87，將每個設(shè)計點的壁面積分“阻力”作為獨立的目標函數(shù)。設(shè)計變量包括9個控制剖面的扭轉(zhuǎn)角和8階CST參數(shù)共171個。設(shè)計約束有升力系數(shù)、低頭力矩系數(shù)和迎角，其中升力系數(shù)的等式約束通過CFD程序的定升力計算實現(xiàn)。幾何約束包括：剖面最大相對厚度為給定值；弦向15%、72%兩處厚度不小于初始值的95%；油箱容積不小于初始值。

本次優(yōu)化種群規(guī)模為512，初始種群在設(shè)計空間內(nèi)隨機產(chǎn)生。每批次同時對128個個體計算評估，每個個體串行計算3個設(shè)計點的氣動性能。計算中同時使用共12 288核CPU的計算資源，CPU型號是英特爾Xeon E5-2692，時鐘頻率為2.2 GHz。在80 h內(nèi)，完成接近20 000個個體的計算分析，遺傳優(yōu)化近40代。圖8展示了優(yōu)化過程中3個目標函數(shù)的演化歷史，從下沿輪廓來看，馬赫數(shù)為0.83、0.85下的兩個設(shè)計目標已經(jīng)趨于收斂，但馬赫數(shù)為0.87目標函數(shù)值的下降趨勢未盡。圖9為馬赫數(shù)為0.85和0.87兩個目標函數(shù)值的散點分布圖(圖中Candidates代表已分析方案)，Pareto前緣上的方案數(shù)量還不夠密集。圖中藍色方形標志為Baseline方案所在位置，綜合考量3個設(shè)計點阻力系數(shù)、力矩系數(shù)、壓力分布和幾何展向過渡等因素后，在Pareto前緣上選定紅色菱形標志的個體(ID19257)為此次優(yōu)化的最優(yōu)結(jié)果。相比于Baseline方案，最優(yōu)解壁面積分流向力系數(shù)下降情況為：馬赫數(shù)為0.83時增加2.0 counts， 馬赫數(shù)為0.85時降低1.5 counts，馬赫數(shù)為0.87時降低8.3 counts。可以選擇馬赫數(shù)為0.83阻力略有增加的方案是因為Baseline方案馬赫數(shù)為0.83的阻力系數(shù)比馬赫數(shù)為0.85的小3 counts左右。

圖10為剖面優(yōu)化后馬赫數(shù)為0.85時機翼表面等壓線分布與Baseline方案對比。盡管優(yōu)化后阻力降低了，但機翼表面等壓線沿展向分布較亂。事實上，由于自動優(yōu)化中各個控制剖面擁有各自獨立的自由度，加上進化類優(yōu)化算法的隨機性，優(yōu)化后機翼幾何展向過渡、等壓線等值后掠、良好的壓力分布形態(tài)很難同時保證。這里使用人工修形設(shè)計加以改進，修形原則是：在優(yōu)化結(jié)果的基礎(chǔ)上微調(diào)設(shè)計變量，保持激波強度和機翼展向升力分布基本不變，達到從翼根到翼梢剖面形狀和壓力分布漸進變化的效果。當然，這對設(shè)計者的經(jīng)驗和能力有較高要求。圖11給出了人工修形設(shè)計后(Refined)的機翼表面等壓線圖，修形后，中外翼部分等壓線基本上是等值后掠的，并且方案在氣動性能上沒有損失。

圖8 優(yōu)化過程中3個目標函數(shù)的歷史Fig.8 History of three cost functions in optimization

圖9 優(yōu)化過程中兩個目標函數(shù)的散點圖Fig.9 Scatter diagrams of two cost function candidates in shape optimization

圖10 馬赫數(shù)為0.85時優(yōu)化前后機翼表面等壓線圖Fig.10 Contours of wing surface pressure coefficient before and after optimization at Ma=0.85

圖11 人工修形后馬赫數(shù)為0.85時機翼表面等壓線圖Fig.11 Contours of wing surface pressure coefficient after refinement at Ma=0.85

圖12給出了機翼弦向15%、72%處相對厚度(t/c)和最大相對厚度沿展向分布的對比情況，優(yōu)化和人工修形后機翼厚度基本無損失。圖13對比了馬赫數(shù)為0.85時優(yōu)化前后4處剖面壓力分布，其中前兩處位于掛架內(nèi)側(cè)的內(nèi)翼部分。優(yōu)化后，激波在內(nèi)翼略有前移和減弱，中段和外翼部分則顯著前移和減弱，這有利于降低激波阻力。人工修形后，激波后流動的二次加速現(xiàn)象完全消失，壓力分布形態(tài)較為理想。圖14對比了馬赫數(shù)為0.85 時機翼展向升力(升力系數(shù)CL與當?shù)叵议Lc的乘積)分布，優(yōu)化前后無明顯變化。

圖12 優(yōu)化前后機翼弦向3處厚度分布對比Fig.12 Comparison of distributions of wing thickness at three chordwise positions before and after optimization

為進一步驗證優(yōu)化取得的阻力下降效果，在優(yōu)化加人工修形后的機翼上安裝第2節(jié)中的通氣短艙進行CFD驗證計算。圖15展示了馬赫數(shù)為0.85時Baseline和人工修形后方案的表面等壓線分布。圖16對比了通氣構(gòu)型馬赫數(shù)為0.87時優(yōu)化前后4處剖面壓力分布，優(yōu)化加人工修形后的主要變化是吸力峰略有提高，激波前移并減弱。表2列出了安裝PN和TFN條件下阻力的變化量，PN構(gòu)型下為壁面積分流向力系數(shù)。馬赫數(shù)為0.87時，兩種條件下的阻力變化量基本一致，但馬赫數(shù)為0.83和0.85時有差異，這樣的差異源于優(yōu)化設(shè)計中是否考慮噴流干擾。如前所述，優(yōu)化前Baseline方案已在馬赫數(shù)為0.85時具備理想的升阻比性能，因此阻力的下降不顯著?？偟膩碚f，本次優(yōu)化設(shè)計達到了預(yù)期效果，特別是在馬赫數(shù)為0.87時，減阻高達全機阻力的3%。本次 優(yōu)化過程持續(xù)80 h， 共消耗計算機時96 萬核時。在精細設(shè)計階段，以如此的時間和計算成本換取這樣的性能收益是很有現(xiàn)實意義的。

圖13 馬赫數(shù)為0.85時優(yōu)化前后機翼剖面壓力分布對比Fig.13 Comparison of distributions of wing sectional pressure at Ma=0.85 before and after optimization

圖14 馬赫數(shù)為0.85時優(yōu)化前后機翼展向升力分布對比Fig.14 Comparison of distributions of wing spanwise lift at Ma=0.85 befor and after optimization

圖15 通氣短艙下優(yōu)化前后機翼表面等壓線圖對比Fig.15 Comparison of contours of wing surface pressure coefficient before and after optimization with TFN integrated

圖16 馬赫數(shù)為0.87時優(yōu)化前后機翼剖面壓力分布對比Fig.16 Comparison of distributions of wing sectional pressure at Ma=0.87 before and after optimization表2 機翼優(yōu)化后阻力系數(shù)變化量Table 2 Variation of drag coefficient after wing shape optimization

構(gòu)型阻力系數(shù)變化量/countMa=0.83Ma=0.85Ma=0.87WBP+PN+2.0-1.5-8.3WBP+TFN+1.3-0.8-8.2

4 結(jié) 論

1) 應(yīng)用動量定理分析動力短艙流場的進排氣流管和外部流動，明確了機翼外形優(yōu)化設(shè)計中以推力與阻力之差最大化或壁面積分“阻力”最小化為目標的合理性。

2) 以寬體客機Baseline機翼方案為對象，研究了短艙/吊掛以及噴流的干擾作用。CFD計算結(jié)果顯示，從WB到WBP+TFN，再到WBP+PN，機翼升力因受干擾而逐步損失，且噴流引起的升力損失大于短艙/吊掛。從升力損失角度看，翼吊布局寬體客機機翼氣動設(shè)計應(yīng)該同時考慮短艙吊掛的安裝效應(yīng)和噴流干擾。

3) 運行前期搭建的優(yōu)化系統(tǒng)，80 h內(nèi)完成了近20 000個設(shè)計方案的計算評估，遺傳優(yōu)化近40代。在滿足所有約束條件的前提下，選擇的最優(yōu)方案在取得馬赫數(shù)為0.87時阻力系數(shù)下降超過8 counts的同時，也使馬赫數(shù)為0.85時的阻力略有降低，阻力發(fā)散性能明顯提高。人工修形設(shè)計后，機翼幾何展向過渡和壓力分布形態(tài)更為理想，且氣動性能無損失。動力構(gòu)型下取得的減阻成果，在通氣短艙構(gòu)型下得到了進一步驗證和確認。