基于非支配排序遺傳算法的渦輪發動機轉子系統裝配參數優化

馮睽睽，張發平，王武宏，張文杰，張田會

(1.北京理工大學 機械與車輛學院,北京 100081；2.北京機電工程研究所 軍事需求與體系設計研究室，北京 100074)

0 引言

轉子系統是發動機的重要組成部分，發動機工作時，轉子系統高速旋轉，在流道內流通高溫高壓空氣和燃氣，使轉子受熱變形。轉子長時間處在該極端環境下，會導致其結構、壽命、可靠性等質量特性指標發生變化。當質量特性動態量值超過指標邊界時，將導致以靜態質量特性裝配的轉子系統難以滿足動態使用要求，引起轉子系統碰磨、松動、不平衡量增大等問題，導致發動機結構出現損傷故障，降低系統工作效率，增加能耗，縮短使用壽命。因此，如何控制轉子系統工作條件下不平衡量的動態變化，是提升轉子系統裝配質量的重要技術工作。

目前，對轉子系統動平衡工藝的研究主要針對汽輪機、普通發動機等零部件。這些零部件質量輕，結構尺寸小，動平衡標準嚴格，在寬徑比≥5且工作轉速大于1 000 r/min時，為抑制轉子工作時質量特性下降，都必須進行動平衡試驗[1]，轉子系統的平衡精度等級一般都在G6.3以下。由于此類轉子一般處于低速或低溫下運轉，高溫與高速同時作用的情況較少，動平衡變化慢，通過動渦旋盤的質心調整，扇區厚度差異設計就能在一定程度上改善動平衡精度和品級[2]。相比于一般的航天發動機，用于導彈發射的彈用渦輪發動機是熱、力、電、磁等惡劣環境條件下長期反復使用的熱力機械，考慮其工作時需要的較大動力，額定工作轉速達到50 000 r/min，工作溫度保持在300～1 000 ℃，實際工作中承受著巨大的熱載荷，導致零件大幅度受熱膨脹變形，材料、結構的不協調引起轉子系統質心的二次偏離，進而影響了轉子系統工作時的動態平衡。此時需要減小熱載下轉子動平衡的變化量，才可提高發動機使用的安全性、可靠性、壽命和效率[3-4]。

本文從某彈用渦輪發動機高壓轉子系統的裝配工藝出發，結合高溫狀態的變形量對系統動平衡的影響，建立發動機轉子熱- 結構理論模型，分析高溫下各零件裝配角對彈用渦輪發動機轉子系統的動平衡變化影響，采用非支配排序遺傳算法(NSGA-Ⅱ)求解使動平衡變化量達到最小的裝配角組合，通過裝配工藝實現高壓轉子系統的動平衡控制。

1 轉子動平衡控制技術

彈用渦輪發動機轉子為剛性轉子，根據轉子動平衡理論，轉子系統工作時其不平衡量主要來自三部分：轉子系統原始不平衡量；高速旋轉下離心力導致的系統變形；溫度升高引起的零件膨脹。轉子系統動平衡試驗的目的是平衡不平衡力和力矩，使轉子能夠穩定運轉。程英輝等[5]通過平衡機的多次增重分析試重在不同質量和位置上引起的振動幅值和相位的變化，得出加重效應與設備平穩性的關系；張軍[6]采用先增重，再利用回轉件在試驗中的微小運動，打緊壓塊螺絲使質心微調，實現動平衡；何振[7]減去了多次試重的繁瑣工序，采取無試重動平衡方法確定風機不平衡質量的相位和滯后角，再旋轉180°相位增加配重；蒲芃成等[8]以低轉速在線動平衡技術，在力自由控制模式下，經過一次試重校正和兩次試轉實現高精度動平衡；張仕海[9]以平衡盤式結構，利用平衡盤與機床主軸的差速運動實現高速旋轉下的動平衡調整。

以上方法可實現動衡量的校正，但都是采用增材或減材的方式，導致轉子原有結構與質量的變化，而轉子系統上有多個零件，各零件裝配角的不同也會校正整體的不平衡量。因此，可通過研究不同裝配角的轉子工作狀態下動平衡演變機理來控制轉子系統的動平衡變化。本文將先試驗測定各零件的初始不平衡量，然后推導出轉子系統的熱- 結構理論模型，再根據NSGA-Ⅱ求解使不平衡變化量和偏移扭矩變化量達到最值時對應的裝配角矩陣，測定試驗中轉子系統在加熱到特定溫度下的變化量范圍，從而對轉配工藝的影響量進行評價。研究框架如圖1所示。圖1中，ΔUt和ΔTt為不平衡變化量和偏移扭矩變化量的理論值，αmax、αmin為使ΔUt和ΔTt達到最值的裝配角矩陣，ΔUm和ΔTm為不平衡變化量和偏移扭矩變化量的試驗模型值。

圖1 動平衡研究框架Fig.1 Framework of dynamic balance research

2 轉子動平衡理論模型

2.1 初始不平衡量

轉子系統在設計上屬于完全軸對稱，轉子材質不均勻、聯軸器不平衡、鍵槽不對稱以及加工時偏差[10]導致轉子系統產生一定的不平衡量。轉子系統剖視圖如圖2所示。圖2中，左起第i個零件的不平衡量記為Ui，則轉子系統n個零件的原始不平衡量集U={U1,U2，…,Ui，…,Un}，n為零件數量。

圖2 轉子系統剖視圖Fig.2 Sectional view of rotor system

轉子系統動平衡試驗中的不平衡量是由各零件工作時的離心力導致的，各零件離心力Fi與不平衡量Ui滿足(1)式：

(1)

式中：Ω為轉子系統轉速；ω為轉子系統角速度。

2.2 系統預裝配工藝

根據圖2可知，各零件裝配時不平衡量的方向不同，即各零件間都有裝配方向角，以軸端螺母為基準，零件i的不平衡量與其方向夾角為αi-1，根據各零件方向夾角推出裝配角矩陣αT=[α1，α2，…，αi，…，αn-1]T。

根據轉子動平衡理論，考慮到零件裝配時裝配方向角的不同，各零件原始不平衡量的中心主慣性軸可能與軸線處于既不平衡又不相交的不平衡狀態，需要選取兩個平面來校正，分別計算兩校正面的不平衡量。在轉子系統中，渦輪盤固定，渦輪軸與渦輪盤以螺釘法蘭連接，與斜流輪熱裝連接；軸斜轉接座與軸流輪和斜流輪過盈配合，并通過軸端螺母壓緊，可選取軸斜轉接座和渦輪盤為校正面，如圖3所示。圖3中，L為兩校正面間距，li為各零件質心與校正面的軸向距離，a和b為兩校正面。將(1)式中計算的各零件的離心力Fi平移到兩個校正面上，得到分力Fi1和Fi2，再計算其合力，從而求解等效在校正面上的不平衡量，Fa、Fb為校正面a、b的合力。

圖3 轉子系統離心力等效圖Fig.3 Centrifugal force equivalent diagram of rotor system

根據力平移理論，離心力Fi與校正面上分力滿足(2)式。若零件i在兩校正面之外，則li為負值；若零件i在兩校正面之間，則li為正值。

(2)

從圖3中可知，平移后兩校正面的各離心力之間的夾角依然為各零件的裝配方向角，根據平行四邊形規則，即可求出Fa和Fb，如(3)式所示。以兩校正面水平方向為參考方向，假設離心力與參考方向的夾角為θi1和θi2，求解合力標量值Fa和Fb如(4)式所示：

(3)

(4)

式中：Ua和Ub為校正面a和b的不平衡量；j、k為1，2，…，n的自然數。

結合(1)式～(4)式，根據試驗測定的各零件質心間距以及各零件的裝配角，計算Ua和Ub，而轉子系統的不平衡量U為兩校正面不平衡量的模之和：

U=|Ua|+|Ub|，

(5)

此時轉子系統所受的偏移扭矩為

T=|Fa·la+Fb·lb|，

(6)

式中：la和lb為校正面a和b與系統質心的距離。

2.3 轉子系統動平衡的熱- 結構理論模型

根據文獻[11]中的輪盤試驗，變形量與回轉件角速度二次方呈正比，與彈性模量呈反比，而且回轉體的尺寸最大，變形量越大。則轉子高速旋轉時各零件的變形量ΔDd,i為

(7)

式中：ρi為第i個零件的密度；μi為第i個零件的泊松比；Ri為第i個零件回轉半徑；Ei為第i個零件的彈性模量。

利用文獻[12]中的熱裝試驗，轉子系統加熱后各零件的膨脹量ΔDe，i采用熱膨脹理論計算，如(8)式所示：

ΔDe，i=λti·Ri·(t-t0)，

(8)

式中：λti為第i個零件在溫度t時的熱膨脹系數；t為溫度點；t0為常溫點。

(7)式和(8)式相加，可以獲得轉子各零件尺寸在高速旋轉和受熱條件下的總變化量，則不平衡量的變化情況如(9)式所示：

(9)

故工作狀態下的各零件不平衡量為

(10)

由以上各零件的不平衡量，結合(1)式計算各零件的離心力；再根據圖3中離心力等效方法求解兩校正面的離心力合力，計算出轉子系統工況下的不平衡量U′和偏移扭矩T′，如(11)式、(12)式所示：

(11)

T′=|F′ala+F′blb|.

(12)

因此，轉子系統動平衡變化量的理論模型為

(13)

3 裝配工藝因素的多目標優化

模型中需要通過求解ΔU和ΔT的最優值來確定轉子系統的裝配角矩陣，屬于多目標問題的優化[13-15]。相比于傳統優化算法從單個初始值迭代求最優解易造成局部最優的局面，遺傳算法從串集開始搜索，覆蓋面大，利于全局擇優，而NSGA-Ⅱ方法降低了非劣排序遺傳算法的復雜性，具有運行速度快、解集收斂性好的優點。因此，本文采用NSGA-Ⅱ[16]求解ΔU和ΔT的最優值及其對應的裝配角矩陣。算法流程如圖4所示。

圖4 算法流程Fig.4 Algorithm flow chart

3.1 初始條件

轉子系統動平衡的熱- 結構理論模型中，各零件的密度不隨工作環境的改變而變化，如表1所示。

表1 轉子系統材質Tab.1 Materials of rotor system

根據2.2節的裝配角矩陣集αT可知，每個裝配角的范圍為αi∈[0°,360°)(以基準離心力順時針方向為正)，算法中將每個裝配角等分為360份，因此共有P=360n-1種組合，初始種群需要從P中抽取。由于傳統的隨機抽取不能均勻提取裝配角矩陣的空間特征，降低了初始種群的覆蓋面，可選用最優拉丁超立方抽樣方法[17]對多維矩陣空間進行分層抽樣，使抽取的樣本盡可能均勻分布，具有較好的均衡性和空間填充性。初始種群P0中α1-α2的抽樣分布如圖5所示。

圖5 初始種群分布圖Fig.5 Initial population distribution

熱- 結構理論模型中轉子系統的ΔU和ΔT是由溫度決定的。表2和表3表明，溫度與零件材料的彈性模量呈反比，與熱膨脹系數呈正比，故ΔU和ΔT在升溫過程中會呈現非線性變化，溫度點t不同，相應的αT最優值也會不同。

表2 轉子各零件材料不同溫度下的彈性模量Tab.2 Elastic moduli of rotor parts at different temperatures

表3 轉子各零件材料不同溫度下的熱膨脹系數Tab.3 Thermal expansion coefficients of rotor parts at different temperatures

3.2 確定多目標函數

在基于動平衡的多目標問題優化中，需要求解最小化不平衡變化量和偏移矩陣變化量對應的裝配角矩陣，則多目標適應度函數為

(14)

式中：Ti為各零件偏移扭矩；ΔU(αT)和ΔT(αT)為裝配角αT對應的不平衡量和偏移扭矩的變化量；gm(αT)為函數F(αT)的約束條件，總數為M；所有設計變量均在定義域內取值，始終滿足F(αT)>0.

3.3 種群優化

3.3.1 種群選擇

由圖4可知：相比于傳統遺傳算法中直接通過遺傳算法來得到新種群，NSGA-Ⅱ在進行遺傳運算前還需先對種群作快速非支配排序[18]，將種群進行分級處理；再計算種群個體的擁擠度，保證種群的多樣性；最后根據排序和擁擠度選擇合適的新種群。

快速非支配排序流程如圖6所示。在種群中，若個體1在所有目標函數的適應度上都不劣于個體2，并且至少在一個目標上優于個體2，則說明個體1可支配個體2. 圖6中，N(i)為支配個體i的個體數，S(i)為被個體i支配的個體數，K為迭代數，K=1，2，…，H.最終，每個個體都有一個非支配排序Rank(K)。

圖6 快速非支配排序流程Fig.6 Flowchart of fast non-dominant sorting

經過快速非支配排序后的種群，已知每個個體的非支配排序，在選擇新父種群時，非支配排序Rank(K)越小，越優先選取。而對于非支配排序相同的個體，則需要計算其擁擠度，如圖7所示。確定目標函數的對應值，畫出包含個體i-1、i和i+1的最小長方形，其長寬和id即為個體i的擁擠度，id越大，越優先選擇。

圖7 擁擠度排序計算Fig.7 Calculation of crowding degree sorting

3.3.2 編碼與運算

按照3.3.1節選擇新的種群Pn，對Pn中的所有個體采用二進制編碼，如新種群中某一裝配角矩陣為αT=[55°,103°,161°,227°,274°]，則αT對應的二進制編碼B(αT)為

(15)

各裝配角取值范圍αi∈[0 rad,2π rad)，故種群中所有個體的元素編碼定義域D∈[000000000,101100111]。編碼完成后，對種群做選擇、交叉和變異的遺傳運算。

經過快速非支配排序和擁擠度的計算后，種群的優劣已經做了區分。選擇Rank(K)排名靠前，id值大的個體直接遺傳到下一代，以保證優良基因的延續。

對于選擇運算剩余的個體，按照自身的交叉概率Pc和變異概率Pm參與交叉和變異運算。每個個體的交叉概率和變異概率隨對應的目標適應度函數值自適應[19]變化，如(16)式、(17)式所示：

(16)

(17)

無論是交叉還是變異，個體所有元素的編碼值都要在D∈[000000000,101100111]定義域內。最后作種群的迭代，直至達到最大的遺傳代數，從而輸出目標函數的最優矩陣集，如(18)式所示：

(18)

4 實例驗證

4.1 算法求解

由圖2可知，整個轉子系統有多個動平衡零件，包括軸端螺母、軸流輪、軸斜轉接座、斜流輪、渦輪軸、渦輪盤等零件，其中高壓渦輪盤與渦輪軸由螺栓法蘭連接，裝配組成高壓轉子。各零件的不平衡量由平衡機測定，平衡機示意圖如圖8所示。

圖8 平衡機示意圖Fig.8 Schematic diagram of balancing machine

采用圖5的值作為模型初始種群；轉子系統在高溫下的材料特性也會發生變化，選取某彈用渦輪發動機穩定工作時的溫度t=600 ℃，擬定常溫點t0=20 ℃，則模型的溫度差t-t0=580 ℃；查表1、表2和表3，選取600 ℃工況時轉子各零件的彈性模量和熱膨脹系數；轉子系統的轉速則選用動平衡試驗的平衡轉速Ω=1 200 r/min. 以上參數值作為模型的初始條件輸入。

根據多目標適應度函數，采用3.3節的方法開展種群優化，在遵循定義域D內變換的原則進行遺傳運算，遺傳代數N選取50代。圖9和圖10分別表示算法迭代過程的三維圖和二維圖。

圖9 不平衡量和偏移扭矩的變化量最小化迭代Fig.9 Minimization iterations of amount of unbalance and variable quantity of offset torque

圖10 不平衡量和偏移扭矩的變化量收斂曲線Fig.10 Convergence curves of amount of unbalance and variable quantity of offset torque

圖11 不平衡量和偏移扭矩的變化量最小化迭代對應的裝配角Fig.11 Minimization iterations of amount of balance and variable quantity of offset torque corresponding to assembly angle

4.2 試驗設計

為了驗證求解的矩陣集的有效性，對轉子系統進行動平衡試驗。試驗步驟如下：

步驟2將轉子系統裝配在平衡機上，啟動平衡機，將轉速調至平衡轉速(1 200 r/min±50 r/min)，由于動平衡機無法達到工作轉速，從而用較低的平衡轉速去模擬工作轉速，測定轉子系統常溫下兩校正面的不平衡量Ua,m和Ub,m，并計算系統初始不平衡量Um.

步驟3根據步驟2中的Ua,m和Ub,m計算兩校正面的離心力，求解系統初始偏移扭矩Tm.

步驟4將轉子系統放在溫控箱中做溫升試驗，溫控箱按照預先設定好的加熱速率加熱至設定的工作溫度600 ℃，溫控箱可以保證溫差在1 ℃以內，加熱完成后做平衡轉速試驗，測定兩校正面在工況下的不平衡量U′a,m和U′b,m，得到系統工況下不平衡量U′m.

步驟5根據步驟4中的U′a,m和U′b,m計算兩校正面的離心力，求解系統工況下的偏移扭矩T′m.

步驟6計算兩種裝配工藝下不平衡變化量以及偏移扭矩變化量max {ΔUt}、max {ΔTt}、min {ΔUt}和min {ΔTt}。

由于試驗中存在誤差，可按照相同輸入條件反復進行多組試驗求證。試驗值與理論值對比如表4所示。

表4 理論與試驗數據對比Tab.4 Comparison of theoretical and experimental data

由表4可知，試驗值接近于理論值，誤差小于9%，驗證了算法的有效性。無論是動平衡模型算法解析還是試驗驗證，轉子各零件裝配角的不同會導致動平衡參數變化量的不同。試驗中，ΔU和ΔT的最大值與最小值分別相差了12.7倍和4.9倍，表明裝配工藝對高溫動平衡變化具有較大影響，合理的裝配工藝(主要是各零件的相對安裝角度)會減小因溫度升高而導致動平衡量的改變，從而改善轉子系統在高溫工況下的質量特性。

5 結論

本文以轉子系統動平衡量的高溫影響為研究對象，建立轉子系統在高溫狀態下的熱- 結構動平衡理論模型，采用NSGA-Ⅱ開展多目標模型優化，并設計動平衡試驗加以驗證，結合試驗測試和模型優化分析了轉子各零件裝配工藝對動平衡質量特性的影響。所得主要結論如下：

1)基于高溫工況下的熱- 結構動平衡理論模型的計算值和試驗測量值基本吻合，表明NSGA-Ⅱ對多目標模型優化的適用性，該算法可以進一步推廣到渦輪發動機其他質量特性指標的優化。

2)基于NSGA-Ⅱ的多目標模型優化方法，揭示了不同零件裝配角對動平衡的影響，實現了大幅減小轉子系統工作狀態下動平衡指標的目的，有利于改善轉子系統的質量特性，可用于改進轉子系統的裝配工藝。

參考文獻(References)

[1] 陳立芳,王維民,高金吉.航空發動機自動平衡技術發展綜述[J].航空動力學報,2019,34(7):1530-1541.

CHEN L F,WANG W M,GAO J J.Summary on the development of auto-balancing technology of aero-engine[J].Journal of Aerospace Power,2019,34(7):1530-1541.(in Chinese)

[2] 叢培田,曹松松,張振厚,等.干式渦旋真空泵的動平衡工藝與振動品質研究[J].真空科學與技術學報,2016,36(6):613-617.

CONG P T,CAO S S,ZHANG Z H,et al.Dynamic balance and vibration suppression of dry scroll vacuum pump [J].Chinese Journal of Vacuum Science and Technology,2016,36(6):613-617.(in Chinese)

[3] 丁鋒,栗祥,韓帥.EEMD與NRS在渦槳發動機轉子故障診斷中的應用[J].航空動力學報,2018,33(6):1423-1431.

DING F,LI X,HAN S.Application of EEMD and NRS in turboprop engine rotor fault diagnosis[J].Journal of Aerospace Power,2018,33(6):1423-1431.(in Chinses)

[4] 游令非,張建國,周霜,等.航空發動機限壽件疲勞可靠度計算新方法[J].航空學報,2019,40(12):109-120.

YOU L F,ZHANG J G,ZHOU S,et al.A new method for calculating fatigue reliability of life limiting parts of aero-engine [J].Acta Aeronautica et Astronautica Sinica,2019,40(12):109-120.(in Chinese)

[5] 程英輝,崔靜波,項青春,等.轉子動平衡中試加重量和位置與振幅和相位關系的研究[J].冶金動力,2019(9):69-70.

CHENG Y H,CUI J B,XIANG Q C,et al.A study on the relationship between trial weight/position and vibration amplitude/phase in dynamic balance of rotor [J].Metallurgical Power,2019(9):69-70.(in Chinese)

[6] 張軍.硬盤動不平衡測量和校正方法研究[D].哈爾濱：哈爾濱工業大學,2016.

ZHANG J.Research of dynamic unbalance detection and adjustment method in HDDs [D].Harbin：Harbin Institute of Techno-logy,2016.(in Chinese)

[7] 何振.風機動平衡無試重校正方法研究[D].哈爾濱：哈爾濱工業大學,2018.

HE Z.Method of dynamic balance without trial weights of draught fan [D].Harbin：Harbin Institute of Technology,2018.(in Chinese)

[8] 蒲芃成,張剴,于金鵬,等.力自由不平衡控制下的高速磁懸浮飛輪系統在線動平衡[J].光學精密工程,2017,25(7):1796-1806.

PU F C,ZHANG K,YU J P,et al.On-line balancing of high- speed magnetic suspended fly wheel system under force free control of imbalance [J].Optics and Precision Engineering,2017,25(7):1796-1806.(in Chinese)

[9] 張仕海.高速機床主軸內置式雙面在線動平衡裝置及關鍵技術研究[D].北京：北京工業大學,2012.

ZHANG S H.Research on built-in,double-face online dynamic balance device for high-speed spindle and its key technologies [D].Beijing：Beijing University of Technology,2012.(in Chinese)

[10] SONG L K,BAI G C,FEI C W.Dynamic surrogate modeling approach for probabilistic creep-fatigue life evaluation of turbine disks [J].Aerospace Science and Technology，2019,95:105439.

[11] 房友龍,劉永葆,賀星.離心力作用下渦輪轉子的徑向變形[J].艦船電子工程,2011,31(2):177-180.

FANG Y L,LIU Y B,HE X.Study of the radial elongations of the turbine rotor with the centrifugal effect [J].Ship Electronic Engineering,2011,31(2):177-180.(in Chinese)

[12] 王水霞,張璞,張金娜.熱載荷下葉輪變形仿真分析[J].風機技術,2015,57(6):54-59.

WANG S X,ZHANG P,ZHANG J N.Simulation analysis on distortion of centrifugal impellers under thermal load [J].Chinese Journal of Turbomachinery,2015,57(6):54-59.(in Chinese)

[13] DECAIX J,BALARAC G,DREYER M,et al.RANS and LES computations of the tip-leakage vortex for different gap widths [J].Journal of Turbulence,2015,16(4):309-341.

[14] 駱清國,趙耀,桂勇,等.基于多種群協同進化免疫多目標優化算法的百葉窗優化研究[J].兵工學報,2019,40(4):689-696.

LUO Q G,ZHAO Y,GUI Y,et al.Research on optimization of louvered fin with hybrid multi-objective optimizations based on EDA and AIS [J].Acta Armamentarii,2019,40(4):689-696.(in Chinese)

[15] 呂國俊,曹建軍,鄭奇斌,等.基于多目標蟻群優化的單類支持向量機相似重復記錄檢測[J].兵工學報,2020,41(2):324-331.

Lü G J,CAO J J,ZHENG Q B,et al.Detection of similar duplicate records based on OCSVM and multi-objective ant colony optimization[J].Acta Armamentarii,2020,41(2):324-331.(in Chinese)

[16] WANG S,ALI S,YUE T,et al.Integrating weight assignment strategies with NSGA-II for supporting user preference multiobjective optimization[J].IEEE Transactions on Evolutionary Computation，2018,22(3):378-393.

[17] ZHOU W P,YANG J F,LIU M Q.Optimal maximin L2-distance Latin hypercube designs [J].Journal of Statistical Planning and Inference，2020,207:113-122.

[18] Hossein Nourianfar,Hamdi Abdi.Solving the multi-objective economic emission dispatch problems using fast non-dominated sorting TVAC-PSO combined with EMA [J].Applied Soft Computing Journal，2019,85:105770.

[19] 張浩為,謝軍偉,張昭建,等.基于混合自適應遺傳算法的相控陣雷達任務調度[J].兵工學報,2017,38(9):1761-1770.

ZHANG H W,XIE J W,ZHANG Z J,et al.Task scheduling of phased array radar based on hybrid adaptive genetic algorithm[J].Acta Armamentarii,2017,38(9):1761-1770.(in Chinese)