民用渦扇發動機總體多學科設計優化研究

羅 鉅 郭福水

中國航發商用航空發動機有限責任公司，上海，200241

0 引言

民用航空發動機的市場競爭非常激烈，其產品在研發時除需考慮燃油經濟性和使用維護成本之外，噪聲和污染物排放等問題也越來越多地受到關注[1]，因此，在發動機設計階段，特別是初始的概念設計階段就有必要綜合考慮除油耗外的污染物排放和噪聲等多項性能指標。

目前，國外主要航空發動機制造商和研究機構均已開發出了適用于航空發動機概念設計階段或初步設計階段的總體多學科設計優化平臺，這些平臺能夠在整個循環參數范圍內進行快速、自動化的多目標方案尋優，從而為部件進一步的詳細設計提供數據支持和優化方向導引。如美國NASA的推進系統數值仿真(numerical propulsion system simulation，NPSS)系統，可以完成發動機不同部件和不同學科的計算分析，進而通過創建數值“試驗單元”來完成發動機的仿真分析[2]。 普惠公司開發出了用于初步設計階段的航空發動機多學科設計優化工具(preliminary multi-disciplinary optimization, PMDO)，該工具集成了零維和一維分析工具，能夠實現熱力循環計算、流道設計、機械設計、質量評估和成本評估的單學科和多學科綜合尋優，也可計算得到部件的氣動效率。羅·羅公司在20世紀70年代便開始了初步階段多學科設計系統Genesis的開發。在2003年，Jones發表了關于Genesis在一臺無人飛行器發動機上的設計應用報告。整個系統主要分析計算發動機性能、結構、單臺成本、飛行性能、研發費用以及維護和處理費用[3]。MTU公司的發動機多學科初步設計工具MOPDES的開發始于20世紀90年代初期，MOPDES包含了熱力、氣動、結構、質量、噪聲和成本等主要評估模塊，以輪擋燃油為評估目標參數[4]。隨著優化計算技術和電子計算機硬件水平的不斷提升，近年來，國內外有關航空發動機各部件的氣動、傳熱和強度等的多學科優化設計研究也在不斷深入[5-7]。

在現階段，我國還沒有完全建立起一套適用于民用大涵道比渦扇發動機的總體多學科設計優化平臺，面臨的主要問題在于缺乏用于民用大涵道比渦扇發動機概念設計階段使用的排放、噪聲和尺寸/質量估算的快速分析工具。為了提高研制效率和質量，縮短研發周期，降低研制風險和成本，構建這樣的平臺很有必要且具有重要的工程實際意義。

本文以某型民用航空發動機概念設計為例，基于協同優化策略構建了發動機總體綜合設計優化平臺，在原有方案基礎上進行了總體方案多目標優化研究，分析了循環參數變化對不同學科指標的影響。

1 物理建模

某民用航空發動機是一型雙軸大涵道比直驅渦扇發動機，由1級風扇、3級增壓級、10級高壓壓氣機、1個環形燃燒室、2級高壓渦輪及6級低壓渦輪組成,見圖1。

圖1 某型渦扇發動機示意圖Fig.1 Illustration of a turbofan engine

本文主要選取該型發動機的風扇、增壓級、高壓壓氣機、燃燒室、高壓渦輪和低壓渦輪這幾個核心部件開展優化研究。其中，風扇、增壓級和低壓渦輪劃歸為低壓轉子系統，高壓壓氣機和高壓渦輪劃歸為高壓轉子系統。在下文中，參數的下標“(0)”表示巡航安裝狀態，即設計點狀態;“(1)”表示標準天起飛狀態；“(2)”表示高溫起飛狀態；“(3)”表示高溫高原起飛狀態；“(4)”表示最大爬升狀態；“(5)”表示起飛著陸循環(land and take-off，LTO)中的起飛狀態；“(6)”表示LTO中的爬升狀態；“(7)”表示LTO中的進近狀態；“(8)”表示LTO中的慢車狀態；“(9)”表示飛越狀態；“(10)”表示進近狀態；“(11)”表示考核風扇與低壓渦輪強度的紅線狀態。

本文中所涉及的主要符號定義見表1,這些量在后文計算時均作為量綱一量。

表1 符號定義Tab.1 Symbol definition

1.1 總體層評估

1.1.1熱力循環分析

熱力循環中的設計點和非設計點性能計算使用的是二次開發版的Gasturb V11.0程序。經修改后，該程序的輸入文件和輸出文件均使用普通文本形式進行讀寫，因此可以實現自動化計算和調用。

1.1.2尺寸/質量估算

發動機尺寸/質量估算采用的是筆者所在公司自主開發的計算程序。該程序采用了一種尺寸和氣動載荷相結合的評估方法，以部件Smith圖為基礎，將發動機尺寸估算劃分為參考點的確定和基于參考點的優化兩個部分。

1.1.3污染物排放估算

燃燒室的NOx、CO、UHC(未燃碳氫化合物)等污染氣體排放估算主要基于國際民航組織(international civil aviation organization，ICAO)數據庫的排放數據，選擇使用T3-p3法計算不同工況下的排放特性。

1.1.4整機噪聲估算

整機噪聲的估算采用基于ANNOP(aircraft noise prediction program)的發動機噪聲預測模型，可以預測發動機不同循環參數方案下的整機和部件的噪聲水平、頻譜和指向性。

1.1.5輪擋燃油估算

本文中輪擋燃油的估算部分采用的是一個基于工程經驗的關系式，具體表達式如下：

F=Δm/(0.115m0)+ΔSw/(75.2Sw0)

(1)

式中，F為輪擋燃油變化率，%；Δm為計算質量與原始方案質量的差值；m0為原始方案的質量;ΔSw為計算加權耗油率與原始方案加權耗油率差值；Sw0為原始方案的加權耗油率。

1.2 部件層評估

1.2.1風扇性能評估

根據工程經驗，風扇部件選擇了33個設計變量，分別是環量分布系數c、環量分布控制點坐標C、葉片攻角I、葉片最大厚度T等。通過NUMECA和ANSYS三維氣動軟件和結構強度軟件獲得不同設計變量組合下的風扇氣動和強度計算結果，基于這些數據生成適合風扇氣動分析的徑向基函數代理模型，由此來快速評估風扇不同設計方案下的氣動性能和靜強度。

1.2.2壓氣機氣動評估

增壓級和高壓壓氣機一維氣動設計采用UNION2程序進行計算。為更形象化地顯示計算結果，基于EXCEL的VBA語言編程，將一維氣動設計的主要輸入和輸出數據顯式化地集成在EXCEL中，完成輸入/輸出數據的自動提取與自動化運算。

1.2.3燃燒室性能評估

燃燒室性能評估模塊的主要工作是針對燃燒室總體性能設計中的初步流道設計的結果進行燃燒室性能估算，估算的性能參數包括燃燒效率、燃燒室總壓損失、出口溫度分布系數(overall temperature distribution factor, OTDF)、徑向溫度分布系數(radial temperature distribution factor, RTDF)和燃燒室污染物排放指標等參數，該一維分析模塊由筆者所在公司自行開發。

1.2.4渦輪氣動評估

高/低壓渦輪一維氣動設計程序TF1D由筆者所在公司自行開發，它根據總體層提供的熱力循環參數，指定各級匹配參數(功率分配、反力度)，通過對已知流道數據進行分析并結合氣動熱力計算和損失迭代，評估已知流量下的相應渦輪特性。

1.2.5部件強度評估

除風扇外，增壓級、高壓壓氣機、高壓渦輪和低壓渦輪的結構強度分析模塊的主要組成部分均相同，分別由參數化建模、自動化分網和強度分析三部分構成，使用UG二次開發程序、Hypermesh和ANSYS實現相應功能。

2 多學科設計優化數學建模

2.1 優化策略

圖2 協同優化策略過程Fig.2 Process of the CO strategy

結合工程設計特點，將總體層優化器作為系統級優化器，其調動總體性能計算涉及設計點計算、非設計點計算、尺寸/質量估算、污染物排放估算、噪聲估算、輪擋燃油估算相關模塊；各子系統級包括風扇、增壓級、高壓壓氣機、燃燒室和高/低壓渦輪各部件的氣動和結構強度分析模塊。系統級優化嵌套子系統級進行優化，其每進行一次仿真分析，各子系統便進行(并行)一次子系統級優化(即需進行多次仿真分析，得到子系統級最優解)。子系統級得到的最優J值向上傳遞到系統級優化器，作為一致性約束條件。

2.2 優化模型

2.2.1系統級優化模型

系統級優化模型中的優化目標、設計變量和約束條件的表述如下：

minF，m，D，O

findB(0)，πF(0)，πI(0)，πc(0)，T4(0)

s.t. πH(0)≤A1，πL(0)≤A2，T3(1)≤A3，

T3(2)≤A4，T3(3)≤A5，T4(1)≤A6，

T4(2)≤A7，T4(4)≤A8，T45(2)≤A9，

JFan≤A10，JIPC≤A11，JHPC≤A12，

JCOM≤A13，JHPT≤A14，JLPT≤A15

需要說明的是，此處的約束值Ai(i=1,2,…)在不同的子系統中對應不同的數值。

2.2.2風扇優化模型

風扇子系統的優化目標為JFan最小，其具體定義如下：

(2)

其中，帶有上標“Fan”的量符號代表風扇子系統分析模塊計算得到的相關參數，不帶上標的相同符號表示系統級中的對應參數(下同)；下標數字含義如第1節中所述；下標符號F代表風扇，下文中量符號下標I、C、COM、H、L分別代表增壓級、燃燒室、高壓渦輪和低壓渦輪各部件。

根據風扇設計經驗，選擇了33個風扇葉型設計變量，分別為：環量分布控制點位置C(i)、環量分布系數c(i)、積疊軸周向偏移控制點位置S(j)、積疊軸周向偏移量s(k)、葉片攻角展向控制點位置I(m)、葉片最大厚度分布展向控制點位置T(m)，其中i=1,2,…,8；j=1,2，…,4；k=1,2,3；m=1,2，…,5。

選擇約束條件為紅線狀態的風扇葉片葉根最大等效應力σF(11)不超過材料許用應力。

2.2.3增壓級優化模型

增壓級子系統的優化目標為JIPC最小，其具體定義如下：

(3)

其中，帶有上標“IPC”的量符號代表增壓級子系統分析模塊計算得到的相關參數。

增壓級氣動設計變量包括：第1級轉子進口處的氣流軸向速度、進口平均半徑處的氣流軸向速度、出口氣流軸向速度、第1級轉子進口的絕對氣流角、中間級反力度、首級加功量因子、出口級加功量因子7個參數；各葉排的擴壓因子和De-Haller數均滿足相應約束條件。

經過對增壓級轉子結構強度優化設計變量的篩選，增壓級轉子結構強度優化中共有10個設計變量，主要包括：增壓級1級和2級各3個尺寸參數(鼓筒內壁寬度、榫槽底厚度和榫槽高度),增壓級3級的4個尺寸參數(鼓筒內壁寬度、輪盤高度、榫槽底厚度和盤心厚度)。增壓級轉子強度計算的約束條件(周向屈服強度的安全系數、徑向屈服強度的安全系數、周向破裂裕度和徑向破裂裕度)均滿足相應的限制條件。

2.2.4高壓壓氣機優化模型

高壓壓氣機子系統的優化目標為JHPC最小，其具體定義如下：

(4)

其中，帶有上標“HPC”的量符號代表高壓壓氣機子系統分析模塊計算得到的相關參數。

高壓壓氣機子系統的氣動設計變量為首級加功量因子、末級加功量因子、首級進口氣流角和中間級反力度。各葉排的擴壓因子和De-Haller數均滿足相應約束條件。

高壓壓氣機轉子結構強度的設計變量為每級輻板厚度、盤心厚度和盤心徑向厚度，10級盤總共30個參數。高壓壓氣機轉子強度計算的約束條件(周向屈服強度的安全系數、徑向屈服強度的安全系數、周向破裂裕度和徑向破裂裕度)均滿足相應的限制條件。圖3所示為高壓壓氣機7級盤的幾何結構設計參數。其中L1、L2、L3為轉子葉片前緣、中部和尾緣的軸向長度；H1、H2、H3為輪盤、葉片前緣根部和葉片尾緣根部的徑向高度；W1、W2為葉片輪盤和輻板的厚度。

圖3 高壓壓氣機R7的尺寸參數Fig.3 Size parameters of the HPC R7

2.2.5燃燒室優化模型

燃燒室子系統的優化目標為JCOM最小，其具體定義如下：

(5)

式中，DCOM為燃燒室污染物排放指標。

燃燒室子系統的設計變量包括火焰筒參考截面系數、前置擴壓器長高比、突擴間隙與進口高度之比、頭部裝置流向長度和火焰筒長高比等14個燃燒室主要幾何和氣動相關參數；總體層與燃燒室模塊的總壓恢復系數差值、燃燒效率差值、起飛/爬升和進近工況的主燃級燃油所占比例和燃燒區空氣比例均滿足相應約束條件。

2.2.6高壓渦輪優化模型

高壓渦輪子系統的優化目標為JHPT最小，其具體定義如下：

(6)

其中，帶有上標“HPT”的量符號代表高壓渦輪子系統分析模塊計算得到的相關參數。

高壓渦輪的氣動設計變量為第1級反力度、第2級反力度和功分配比3個參數;高壓渦輪的氣動約束條件為每級出口氣流折轉角、出口馬赫數、各級軸速比和出口氣流角滿足相應約束。

高壓渦輪的結構強度設計變量為輻板的4個位置厚度、過渡圓弧角度(5個)、盤心厚度和盤心徑向厚度(4個)，共13個參數。高壓渦輪的強度約束條件為周向屈服強度的安全系數、徑向屈服強度的安全系數、周向破裂裕度和徑向破裂裕度滿足相應約束。圖4所示為高壓渦輪第1級轉子結構強度設計的尺寸參數。

圖4 高壓渦輪R1的尺寸參數Fig.4 Size parameters of the HPT R1

2.2.7低壓渦輪優化模型

低壓渦輪子系統的優化目標為JLPT最小，其具體定義如下：

(7)

其中，帶有上標“LPT”的量符號代表低壓渦輪子系統分析模塊計算得到的相關參數。

低壓渦輪的氣動設計變量為每級反力度和功分配比，總共12個參數；氣動約束變量為每級葉片的升力系數、每級轉子的折轉角和出口氣流角，均滿足相應的約束條件。

通過靈敏度分析，綜合判斷以上參數對結構應力和質量的貢獻程度，最后確定選取輻板厚度、盤心厚度、盤心高度和輪盤內徑4個尺寸參數作為輪盤結構強度優化的設計變量。因此，低壓渦輪轉子結構強度優化中共有24個設計變量，強度約束條件為周向屈服強度的安全系數、徑向屈服強度的安全系數、周向破裂裕度和徑向破裂裕度均滿足相應的約束條件。

3 MDO平臺的搭建

3.1 Isight優化平臺

本文中的發動機總體多學科設計優化平臺通過商用計算機輔助優化軟件Isight進行搭建。該軟件具有柔性靈活的仿真流程集成功能，用戶可以通過其Simcode程序基礎組件實現不同程序的輸入、執行和結果提取。其自帶的流程組件Process Component可靈活地將多個應用程序集成在一起，從而形成復雜的工作流程，實現程序間的參數自動傳遞。當選取多目標優化算法后，Isight軟件的Optimization模塊能夠對計算結果自動分析，在整個設計空間內自動迭代獲得相關Pareto解[8]。

圖5所示為基于Isight軟件搭建的發動機總體多學科設計優化平臺架構示意，由系統層傳入各部件層設計優化所需要的總體參數，如進口總溫、總壓、流量、轉速、總壓比、膨脹比等氣動參數以及各部件進出口內外徑、軸向長度、各級葉片子午面軸向長度、半徑和軸向間距等幾何信息，各子系統層使用2.2節所述的相應數學模型進行氣動和強度的設計優化，其優化目標為其子系統的一致性約束值最小。最后，將優化計算得到的最小一致性約束值返回給系統層(總體層)作為該層的一個約束變量。若滿足約束，則說明此時部件設計能夠達到滿足總體性能優化時的方案指標。

3.2 優化算法

在系統層的優化器中，使用的是非支配排序遺傳算法(non-dominated sorting genetic algorithm，NSGA-Ⅱ)，這是應用最為廣泛和成功的一種多目標優化算法[9-11]。NSGA-Ⅱ 算法主要包括編碼、種群初始化、適應度評估、選擇、交叉及變異六個步驟。與此同時，NSGA-Ⅱ算法引入了快速非支配排序算法和精英策略，采用了擁擠度和擁擠度比較算子，以便確保種群的多樣性。為了加快收斂，縮短計算時間，各子系統優化器選用的是自適應模擬退火算法(adaptive simulation annealing，ASA)，它比傳統的模擬退火算法具有更優良的全局求解能力和計算效率[12]。

4 計算結果和分析

使用Isight軟件建立的發動機總體多學科設計優化平臺，對某型民用大涵道比渦扇發動機進行了方案多目標尋優。采用了DELL T3600臺式工作站執行計算，該工作站配置了4個主頻為2.8GHz的英特爾 Xeno(R) E5-1603 CPU，其安裝內存為16G。系統層計算選用上述NSGA-Ⅱ優化算法，子系統層各部件優化器選用自適應模擬退火算法，各部件子系統優化器計算步數設置為30步。

圖5 發動機MDO平臺架構Fig.5 Structure of the engine MDO platform

根據前人對NSGA-Ⅱ算法的參數化研究發現，在步數較小且相同時，適當加大種群數，優化結果相對更好[12-14]，因此設置種群數為32，代數為15。由于初始方案已經是對應單目標的最優解，且尋優范圍較小，因此暫不對交叉率和變異率等參數進行修改，保持為默認值。系統層最終計算步數為57步，計算總時間為68 h，總共得到了33個可行解，其中Pareto解為12個。

從12個Pareto解中選取對應不同優化目標分別最佳的6組方案，如表2所示。由于初始方案已經是前期優化得到的結果，此次優化得到的結果相對而言變化較小，其中不同優化目標的最優結果分別為：整機質量減小3.39%，輪擋燃油變化率下降0.18%，NOx排放指標下降1.02%，飛越工況噪聲下降0.13%，邊線工況噪聲下降0.19%，進近工況噪聲下降0.13%。此外值得注意的是，經過對比分析可以發現，最優輪擋燃油變化率和污染物排放指標在同一個方案下取得，此時涵道比接近于上限值；三個工況噪聲的變化趨勢基本上是一致的，例如飛越工況和進近工況噪聲的最小值均對應于同一個循環參數方案組合。

對于實際工程設計而言，如果污染物排放和噪聲達到適航要求的相應標準，則最為關注的是燃油消耗和發動機整機質量，因為這兩個指標直接關系著經濟性與是否滿足飛機制造方的技術指標要求。因此本文著重圍繞總體循環參數對輪擋燃油變化率和整機質量的影響展開分析。

需要指出的是，表2中的初始方案為主要考慮耗油率最低的單目標優化結果。

圖6所示分別為包含Pareto解在內的可行解中加權耗油率和整機質量隨涵道比的變化。其中圓點代表可行解，正方形的點代表Pareto解，三角形的點代表原型方案。可以看到，加權耗油率隨涵道比的增大整體呈下降趨勢，整機質量則隨涵道比的增大呈上升趨勢。這主要是由于涵道比增大使得外涵流量增大，發動機的推進效率得以提高，從而使得經濟性得到改善；而外涵流量增大意味著低壓流道流量增大，部件尺寸增大，導致整機質量增大。

圖7所示分別為包含Pareto解在內的可行解中加權耗油率和整機質量隨渦輪前溫度的變化。圖7a顯示出加權耗油率隨渦輪前溫度的變化總體呈下降趨勢，這說明在原型方案附近的這段渦輪前溫度區間內，耗油率是隨渦輪前溫度單調遞減的；相比加權耗油率而言，整機質量隨著渦輪前溫度的升高也呈現出下降趨勢，但較涵道比的影響相對較弱，Pareto解前沿未出現明顯的變化規律。

表2 優化得到的Pareto解Tab.2 Optimized pareto solutions

(a)涵道比和總耗油率關系

(b)涵道比和整機質量關系圖6 加權耗油率/整機質量與涵道比的變化關系Fig.6 Changing relation between S and B, m

圖8所示為輪擋燃油變化率和整機質量的關系。從圖8可以看到，輪擋燃油變化率和整機質量之間呈現出近似二次曲線的分布規律。當整機質量在1.005附近時，存在輪擋燃油最小的情況；若維持輪擋燃油變化率不變，存在使得整機質量減小的循環參數組合(更小的涵道比和渦輪前溫度)，但噪聲和排放指標會有所上升，這是受涵道比與渦輪前溫度影響的綜合結果。

(a)渦輪前溫度和總耗油率的關系

(b)渦輪前溫度和整機質量的關系圖7 加權耗油率/整機質量與T4的變化關系Fig.7 Changing relation between T4 and S, m

圖8 整機質量和輪擋燃油變化率的關系Fig.8 Changing relation between m and F

5 結語

本文使用Isight多學科設計優化軟件搭建了基于協同優化策略的民用大涵道比渦扇發動機總體MDO平臺。通過使用NSGA-Ⅱ優化算法，開展了以輪擋燃油、整機質量、噪聲和NOx排放指標為優化目標的發動機總體方案設計優化。相對于耗時為10 d左右的常規手動迭代方法，該平臺能夠在較短時間內快速、自動化地在整個熱力循環參數設定范圍內進行多目標尋優，對應于總體層的每一個方案，風扇、增壓級、高壓壓氣機、燃燒室、高壓渦輪和低壓渦輪各部件優化器以其一致性約束最小為優化目標，分別進行氣動和結構強度的部件自動優化設計，尋求到滿足氣動和強度約束條件下的部件最優解。該平臺能夠基本滿足概念設計階段發動機總體方案設計需求。通過對優化結果的分析研究，得出以下結論：

(1)基于NSGA-Ⅱ多目標遺傳優化算法的計算具有“先優化后決策”的優點，針對預先設定的多個優化目標，可以得到一組Pareto最優解，供設計者最終決策。

(2)優化結果表明，在原型方案總體循環參數附近的范圍內，對發動機輪擋燃油和整機質量影響最大的是涵道比和渦輪前溫度，而它們對噪聲的影響程度相對較小。涵道比和渦輪前溫度的綜合作用使得Pareto解中的輪擋燃油變化率和整機質量呈二次曲線的變化規律，當整機質量在區間內時，對應有最小的輪擋燃油變化率。