基于非支配排序遗传算法的涡轮发动机转子系统装配参数优化

冯睽睽，张发平，王武宏，张文杰，张田会

(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)