動力干擾下寬體客機機翼多目標優化設計

薛幫猛，張文升，孫學衛，吳宇昂

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

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

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

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

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

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

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

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

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

(DFrame+DPre+DPost)

(1)

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

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

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

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

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

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

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

2 優化設計方法

本節介紹優化設計中涉及的幾何外形參數、計算網格生成、CFD計算和優化流程等方面。

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

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

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

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

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

圖7 優化設計流程Fig.7 Flow chart of optimization design

3 機翼優化設計

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

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

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

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

圖8 優化過程中3個目標函數的歷史Fig.8 History of three cost functions in optimization

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

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

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

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

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

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

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

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

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

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

構型阻力系數變化量/countMa=0.83Ma=0.85Ma=0.87WBP+PN+2.0-1.5-8.3WBP+TFN+1.3-0.8-8.2

4 結 論

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

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

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