基于改进粒子群算法的非对称传动主轴多目标优化

张运涛， 李以农， 张志达， 罗法氿， 王 成

(1. 重庆大学 机械传动国家重点实验室，重庆 400030；2. 北方车辆研究所，北京 100072)

传动轴是旋转机械传递动力和运动的关键零件，其使用性能对整个系统的功能实现和可靠性起着至关重要的作用[1-2]。对于特种车辆，由于传动系统固有的非对称结构，当车辆启停、换挡和转向时，会使传动轴两端受到非对称交变载荷和冲击载荷作用，往往会造成传动轴一端因为受到过大的载荷而出现疲劳断裂，导致整个车辆无法正常工作[3-5]。可见，传动轴结构参数是影响特种车辆传动性能的关键。因此，对传动轴的结构参数优化具有重要的工程意义。

目前，对于结构优化问题，通常采用解析法和数值计算法两种[6]。解析法即利用数学分析的方法，根据目标函数导数的变化规律与函数极值的关系，求取目标函数的极值点。但当目标函数比较复杂或为非凸函数时，很难应用解析法求取由目标函数各项偏导数组成的方程组，故解析法现在已较少使用。数值计算法是以适当步长沿着目标函数下降的方向，逐步逼近目标函数最优解的方法。Huang等[7]建立了螺旋桨轴弯扭耦合振动模型，利用高阶Runge-Kutta法进行数值计算完成了螺旋桨轴的结构尺寸优化。Wang等[8]通过数值计算建立了轴距变化与滚动阻力的关系，结果表明优化轴距可以提高车辆的牵引性能，并仿真验证了优化结果的有效性。随着近代电子计算机计算能力的飞速提升，数值计算法得到了更加广泛的应用。Hong等[9]针对超精密磨床静液压主轴设计提出了一种人机集成设计框架，并对主轴进行热力学性能优化和动态优化，极大地提高了主轴的性能。杜官将等[10]利用APDL建立机床主轴的参数化有限元模型，以主轴的重量为优化目标，应用ANSYS优化设计功能对主轴的支承跨距、外径和悬伸长度等进行了优化设计，结果表明在保证机床各种性能的前提下，主轴重量得到有效的减小。

上述的优化方法在求解单目标优化问题时具有收敛速度快、搜索效率高等优点，但对于具有复杂梯度信息的问题，则无能为力；同时，面对多峰、多模态的优化问题，传统优化方法会很快陷入局部最优解，且很难跳出[11]。针对上述问题，众多智能算法应运而生，如人群搜索法、遗传算法、粒子群算法和蚁群算法等[12]。Yi等[13]根据主轴电机与传动系统之间的能量流，建立了主轴系统的能量模型，以结构体积最小和单位体积能耗最小为目标，建立了主轴电机和传动系统多目标参数优化模型，基于多目标机器学习优化算法求解Pareto最优解并通过实例验证了所提优化方法。朱成实等[14]针对传统优化设计方法在解决主轴优化设计中出现的问题，引入惯性权重值适应性递减的粒子群算法，结果表明主轴结构参数优化效果明显。陈东菊等[15]以主轴系统参数作为设计变量，利用遗传算法对液体静压主轴系统的运动误差进行优化分析，优化后主轴系统横、垂向误差运动及主轴倾角的优化效率分别为41.22%，25.21%和66.16%，优化效果显著。

现有文献对传动轴的结构优化大都是基于传动轴具有对称的结构布置，传递载荷时载荷均匀分配。然而，特种车辆由于综合传动系统结构布局限制，传动主轴两端的输出为非对称结构，系统刚度和阻尼等参数存在明显差别，导致左右轴段抗振能力存在较大差别，同时系统输入载荷也会不被左右两侧输出均分[16]。

有鉴于此，本文综合考虑传动主轴-轴承系统内外多源激励，采用集中质量法建立系统非线性振动模型，利用Runge-Kutta法进行数值求解，获得稳态工况下系统弯扭耦合振动响应以及动载荷和振动能量的分布特点。基于改进的粒子群优化 (particle swarm optimization，PSO)算法，以系统的振动载荷和振动能量为优化目标对系统结构参数进行多目标优化，并通过仿真验证了优化结果的有效性，优化结果为非对称传动系统的结构优化问题提供了重要的指导依据。

1 主轴-轴承系统动力学模型

某特种车辆的传动主轴是非对称布置的三转子盘结构，中间转子盘为复合行星排的行星架，作为主轴结构的动力输入结构；左右两边的转子盘为左、右侧行星汇流排的齿圈，作为主轴系统的输出结构分别连接左、右侧的车轮。考虑输入/输出、轴承支撑等因素的影响，基于集中质量法简化建立了如图1所示的主轴-轴承系统弯扭耦合动力学模型[17-19]。

图1 主轴-轴承系统动力学模型Fig.1 Dynamic model of spindle bearing system

图1中：mi和Ji分别为三个转子盘(包括轴承内圈)的集中质量和转动惯量；mpij为各转子盘两端轴承座(包括轴承外圈)的集中质量；kt1,2和ct1,2为主轴两端的扭转刚度和扭转阻尼；ks1,2和cs1,2为主轴两端的弯曲刚度和弯曲阻尼；kij为各转子盘两端轴承的接触刚度；cij为各转子盘两端轴承的结构阻尼；kp和cp为轴承座的支持刚度和阻尼；l1和l2分别为主轴左右端轴段长度，由于结构布局限制，主轴为非对称结构，l1明显大于l2。文中，i=1,2,3依次表示转子盘1、转子盘2、转子盘3；j=1,2依次表示转子盘左、右两端。下文出现下标i，j，定义相同。

1.1 转子盘动力学模型

对系统单个转子盘进行受力分析，考虑转子盘的偏心距，其受力分析如图2所示。

图2 转子盘受力分析图Fig.2 Stress analysis diagram of rotor disc

图2中：oi为转子盘几何中心，坐标为(xi,yi)；c为转盘质心，坐标为(xc,yc)；o为坐标原点；e为偏心距；ω为转子盘转速；θ为转子盘扭角；φ=ωt+θ为转盘转角。根据几何关系可以得到质心c的坐标为

(1)

1.2 轴承系统动力学模型

系统中轴承采用深沟球轴承，考虑轴承内外圈水平和竖直方向振动，建立4自由度动力学模型[20-21]，如图3所示。图3中：xin,o和yin,o为轴承内、外圈在轴承横截面内水平方向和竖直方向的振动位移；min为内圈质量；mo为外圈质量；kc为Hertz接触刚度；ko为外圈刚度；cin为内圈阻尼；co为外圈阻尼。

图3 深沟球轴承动力学模型Fig.3 Dynamic model of deep groove ball bearing

根据牛顿第二定律，可得轴承系统的动力学方程为

(2)

式中：fx为x向接触力；fy为y向接触力；Fr为施加在轴承内圈的y向载荷。

假设轴承外圈固定在轴承座上，内圈转动时，滚珠在滚道上作纯滚动。设滚珠个数为N0，则第n个滚珠经过时间t后转过的角度θn为

(3)

式中：ωc=ωind/(Dzc+dzc)为保持架角速度；ωin为内圈转速；dzc为轴承内径；Dzc为轴承外径。

轴承第n个滚珠在角位移θn处与滚道的径向接触变形量为

δn=(xin-xo)cosθn+(yin-yo)sinθn-s

(4)

式中，s为径向游隙。

根据Hertz接触理论可得滚珠与滚道之间的非线性接触力为

(5)

式中：Hn为接触系数，表征非线性接触力是否存在； 当δn>0时，Hn=1，否则Hn=0；m为常数，对于球轴承m=3/2。

综上得轴承在水平和竖直方向总的接触力为

(6)

1.3 主轴-轴承系统弯扭耦合动力学模型

考虑三个转子盘集中质量点在自身平面内水平方向和竖直方向上的自由度x和y以及随主轴轴线扭转的自由度θ，根据质心运动定理和动量矩定理，建立主轴-轴承系统的动力学微分方程如下：

转子1的动力学微分方程为

(7)

转子2的动力学微分方程为

(8)

转子3的动力学微分方程为

(9)

转子处轴承座的动力学微分方程为

(10)

式中：T2为转子盘2的外部激励扭矩，即系统输入扭矩；T1和T3分别为转子盘1、转子盘3的负载扭矩，即系统输出扭矩；Tjl4和Tjl8分别为4倍和8倍发动机转频正弦激励扭矩；fxij/fyij，fcxij/fcyij分别为轴承处的x/y向接触力和x/y向结构阻尼力，结构阻尼力表达式为

(11)

式中，cij为各轴承的结构阻尼。

2 主轴-轴承系统振动特性分析

本文研究的某特种车辆传动主轴-轴承系统在一挡起步工况下出现主轴右端断裂故障，如图4所示。为探究主轴断裂机理，开展一挡工况下的主轴振动特性分析，采用Runge-Kutta法对系统动力学方程进行求解，取计算时间步长2×10-6s。一挡工况参数如表1所示。给定主轴结构几何参数如表2所示。系统支撑轴承依据文献[22]计算可得轴承系统各仿真参数。

图4 主轴断裂故障件外观Fig.4 Appearance of failure parts of spindle fracture

表1 系统一挡工况参数表Tab.1 System first gear working condition parameter table

表2 主轴结构几何参数表Tab.2 Table of geometric parameters of spindle structure

2.1 主轴-轴承系统扭转振动特性分析

(12)

(13)

式中：Tb1，Tb3分别为左右轴段的波动扭矩；WT1，WT3分别为左右轴段的扭转振动能量。

一挡工况下，主轴的扭转角位移、扭转角加速度、波动扭矩和扭转振动能量，如图5所示。该工况下，主轴左右轴段扭转振动响应有较大差异，左右轴段的扭转角位移幅值和扭转角加速度几乎相同，但左端轴段由于扭转刚度较小在偏离平衡位置更大的位置来回扭转；右端轴段的波动扭矩峰值达到900 N·m，远大于左端轴段的400 N·m，这是导致主轴右端断裂的一个主要原因；同时，主轴右段的扭转振动能量峰值达到200 W，也大于左端轴段的扭转振动能量峰值。由此可见，一挡工况下主轴右端轴段的扭转振动剧烈程度要远大于主轴左端。

图5 主轴扭转振动特性Fig.5 Torsional vibration characteristics of spindle

2.2 主轴-轴承系统弯曲振动特性分析

与扭转振动特性分析相同，将相邻转子盘质量点x，y向位移和加速度相减便可得到主轴轴段x，y向位移和加速度，同时引入弯曲振动动载荷和弯曲振动能量，定义为

(14)

(15)

(16)

式中：f1，f3分别为左右轴段弯曲动载荷；W1，W3分别为左右轴段弯曲振动能量。

一挡工况下，主轴的弯曲振动位移、加速度、弯曲动载荷和弯曲振动能量，如图6所示。在该工况下，主轴左右轴段的弯曲振动亦有明显差异。由于左端轴段的弯曲刚度较小，左端轴段x，y向位移都要大于右端轴段位移，且左端轴段的横向弯曲加速度要略大于右端轴段横向弯曲加速度；右端轴段的弯曲动载荷峰值达到150 N，而左端轴段的弯曲动载荷峰值仅为15 N，虽然右端的横向弯曲位移较小，但右端轴段的弯曲刚度较大，可见轴段弯曲刚度是影响弯曲动载荷的主要因素；同时，主轴右段的弯曲振动能量也显著大于左端轴段的扭转振动能量。

图6 主轴弯曲振动特性Fig.6 Bending vibration characteristics of spindle

考虑转子盘质量偏心，转子盘弯曲振动与扭转振动存在弯扭耦合效应，且随着偏心距的增加耦合效应愈明显。综合系统弯扭耦合振动可以得出结论：系统在一挡起步工况下需要传递大扭矩，扭转振动响应要远远剧烈于弯曲振动响应；右端波动扭矩峰值达到了900 N·m，而右端弯曲动载荷峰值仅为150 N，扭转振动能量的数值也远远大于横向振动能量数值，由此可见主轴是“扭断”的。对比Zhao等对同一研究对象开展的有限元仿真研究，如图7所示，当主轴的输入扭矩为26 000 N·m时，左右轴段上的应力分布值分别为285.9 MPa和502.6 MPa，右侧轴段的应力值是左侧的2倍左右，与本文数值计算得到的结论基本一致。为进一步验证本文仿真结果的正确性，对主轴进行断口检查、金相观察和拉伸冲击测试等试验，如图8所示，得出主轴的断裂模式为扭转疲劳断裂，裂纹起源于距端面约50 mm处的花键底部，周向多源萌生向心部拓展，并最终断裂，亦验证了本文数值仿真所得结论。

图7 主轴应力分布图Fig.7 Stress distribution of main shaft

图8 断裂萌生与拓展示意图Fig.8 Sketch map of fracture initiation and development

3 改进PSO算法

3.1 标准PSO算法

(17)

(18)

式中：ω为惯性因子；t为迭代次数；vij为粒子速度，xij为粒子位置；c1，c2为学习因子；r1，r2为介于[0,1]的随机数。

3.2 改进PSO算法

PSO算法针对不同问题优化时，无法保证每一次的结果都收敛到全局最优解，易陷入局部最优解。为了提高PSO算法的全局搜索能力和搜索效率，对标准粒子群算法进行如下改进：

3.2.1 添加递减惯性因子

较大的惯性因子有利于提高算法的全局寻优能力，反之则可以保证粒子在最优解附近精细搜索。为了平衡算法的全局和局部的搜索能力，在迭代过程中线性地减小ω的值，定义为[24]

(19)

式中：t为粒子群当前迭代次数；tmax为粒子群总迭代次数；ωmax和ωmin分别为最大、最小惯性因子。

研究表明，对于多数优化问题，在ωmax和ωmin相同情况下，凹函数递减策略优于线性策略[25]，定义为

(20)

3.2.2 改进学习因子

(21)

3.2.3 变异操作

基于标准PSO算法引入一种新颖的变异操作，在粒子进行前一半迭代的过程中，对上一步得到的粒子群进行变异操作，随着迭代次数的叠加，变异操作对粒子群的影响随之减小，如图9所示，用以提高粒子的全局搜索能力。变异操作的伪代码如表3所示。

图9 受变异操作影响粒子所占百分比Fig.9 Percentage of particles affected by mutation operation

表3 变异操作伪代码Tab.3 Variation operation pseudo code

4 主轴结构参数多目标优化

基于所提出的改进PSO算法对传动主轴结构参数进行优化设计。将主轴结构参数优化看作为一个具有n个优化目标和u个约束条件的多目标非线性优化问题，则多目标参数优化问题描述为

Ffitness(x)=best{f1(x),f2(x),…,fn(x)},x∈T

(22)

T={x∈Rm,G∶gi min≤gi(x)≤gi max,i=1,…,u}

(23)

式中：f1(x),f2(x),…,fn(x)分别为第1～第n个优化设计目标；T为约束条件；Rm为设计变量域；G为变量约束；n为设计目标维数；m为设计变量维数；u为约束个数。

4.1 设计优化目标

主轴断裂的主要原因为：主轴左右轴段动态扭矩和扭转振动能量由于结构不对称相差较大，主轴右端扭转振动响应较左端轴段剧烈。因此，选取主轴左右轴段波动扭矩差值和振动能量差值为优化目标，同时为防止主轴扭转振动更加剧烈，将主轴左右轴段波动扭矩总和及振动能量总和作为优化目标：

目标函数1

f1(x)=min(T3-T1)

(24)

目标函数2

f2(x)=min(WT3-WT1)

(25)

目标函数3

f3(x)=min(T1+T3)

(26)

目标函数4

f4(x)=min(WT3+WT1)

(27)

4.2 设计优化变量

通过灵敏度分析发现，为了降低波动扭矩和扭转振动能量，应该先优化扭转刚度kt2和转动惯量J2，其次为扭转刚度kt1和轴承游隙。但转动惯量J2为累加的整车惯量无法改变，轴承游隙为装配工艺所决定，因此应优先考虑优化主轴的扭转刚度降低波动扭矩和扭转振动能量。

扭转刚度计算公式为[26]

(28)

(29)

由上式可知，扭转刚度K与材料的剪切弹性模量G、轴段横截面对圆心的极惯性矩Ip和轴段长度l有关，但主轴常用材料的剪切弹性模量大致相同，同时系统中主轴轴段长度受到结构限制不可更改，因此只能够通过优化主轴内径和外径来改变Ip来实现主轴的扭转刚度优化。当主轴的外径D保持不变，内径d即使取到0.5D这样大的数值时，Ip依然为原来的15/16，表明轴段内径对极惯性矩的影响很小，应优先考虑优化主轴外径D来寻得系统的最佳结构参数。同时，对于主轴左右轴段的扭转刚度来说，当外径相同时，扭转刚度的比值趋近于一个定值，那么结构不对称带来的抗振能力差异和载荷分配不均就会一直存在，因此应考虑将主轴的左右轴段外径取不同的数值，以实现最优的主轴扭转振动响应。

综合考虑，设计优化变量为

x=[dldrDlDr]

(30)

式中，dl,r和Dl,r分别为主轴左右轴段的内径和外径。

4.3 设计变量约束

优化变量的基准值为原始主轴结构参数，由减少主轴扭转振动波动扭矩和振动能量为要求，确定变量dl，dr，Dl，Dr的取值范围，约束条件为

dl,dr∈[0,20]，Dl,Dr∈[50,70]

(31)

传动主轴结构在稳态工况下，由于内外激励共同作用，承受多种交变载荷，在进行优化设计时必须考虑主轴危险部位的应力和应变，确保其有足够的强度和刚度，因此需要约束主轴危险截面的应力和应变小于许用值。为此引入惩罚函数，表示为

(32)

式中：Fi(x)为新的目标函数；M为罚因子，是一个正常数；m为约束条件个数；gj(x)为约束条件。当M充分大时，Fi(x)的最优解能逼近原始约束问题的最优解。

4.4 多目标优化框架

利用MATLAB软件编写主轴模型和改进PSO算法，在改进PSO算法中嵌入主轴模型进行求解，计算主轴波动扭矩和扭转振动能量作为优化算法的目标函数。多目标优化的设置为：种群大小N为100，最大迭代次数为tmax为100，最大惯性因子ωmax为0.9，最小惯性因子ωmin为0.4，变异率Pm为0.5。本文所提出的多目标优化流程如图10所示。

图10 基于改进PSO的多目标优化流程Fig.10 Multi objective optimization process based on improved PSO

4.5 Pareto解集选优

(33)

定义支配函数φk，第k个解的支配值为

(34)

式中：l为Pareto解集中解的数目，根据模型的优化结果得到；n为主轴优化设计目标个数。

由支配函数φk计算公式，可得Pareto解集中每个非劣解的支配值，支配值反应了该解的综合性能，选择具有较大支配值的解为最优解。

4.6 优化结果分析

优化所得Pareto解集采用支配值计算公式计算得到的支配值，如图11所示。第18号粒子具有最大支配值，对应主轴设计变量及其优化前后数值，如表4所示。

图11 Pareto最优解集的支配函数值Fig.11 The dominating function value of Pareto optimal solution

表4 设计变量优化结果Tab.4 Optimization results of design variables

根据表4所得主轴优化前后结构参数，对优化前后主轴波动扭矩和扭转振动能量进行对比分析，分别如图12和13所示。

由图12和图13分析可知，主轴经过参数优化后，系统的波动扭矩和扭转振动能量整体上来看得到明显减小，同时主轴两端轴段振动不平衡得到明显改善。右端的波动扭矩峰值从900.60 N·m降低到646.60 N·m，左端的波动扭矩峰值有所增加，从378.90 N·m增加到644.80 N·m；右端的扭转振动能量峰值从207.00 W下降到178.30 W，左端的扭转振动能量峰值几乎保持不变。从对扭转振动的优化效果来看，使主轴整体性能达到最优的多目标优化目的已达到，但是无法避免左端轴段的局部扭转动载荷有所增加。

图12 主轴波动扭矩优化前后对比Fig.12 Comparison of spindle fluctuating torque before and after optimization

图13 主轴扭转振动能量优化前后对比Fig.13 Comparison of torsional vibration energy before and after optimization

虽然主轴是在一挡起步工况下出现的主轴右端断裂故障，但仅仅针对一挡工况下的系统扭转振动进行优化可能会加剧主轴系统在其余工况下的振动响应，因此考虑主轴常用工作挡位的多工况优化设计。系统常用挡位工况参数，如表5所示，经支配函数计算公式计算得到的各工况下优化结果，如表6所示。

表5 主轴-轴承系统各挡位工况参数Tab.5 Working condition parameters of spindle bearing system in each gear

表6 设计变量优化结果对比Tab.6 Comparison of optimization results of design variables

图14 波动扭矩随左右轴端扭转刚度变化规律Fig.14 Variation of fluctuating torque with torsional stiffness of left and right shaft ends

5 结 论

综合考虑传动主轴-轴承系统内外多源非线性激励，采用集中质量法建立了系统非线性振动模型，通过数值计算得到了系统在一挡工况下的弯扭耦合振动响应，并基于此分析了主轴断裂的主要原因。针对传统PSO算法易陷入局部最优问题，对PSO的参数进行自适应调整并引入变异操作，提出一种改进的PSO算法。基于所提出改进的PSO算法对非对称主轴结构进行了多目标优化设计，得到了以下结论：

(1) 传动主轴-轴承系统虽然存在由于转子质量偏心引起的弯扭耦合效应，但主轴承受的载荷主要是扭矩，系统扭转振动响应程度要远远剧烈于弯曲振动响应；而由于布局限制，主轴两侧为长短不同的非对称非等强度结构，右侧轴段在工作时要承担更大的输出扭矩，其波动扭矩幅值为左侧的2倍左右，这使主轴右端轴段更早的达到扭转疲劳极限进而发生断裂。

(2) 针对传动主轴右侧断裂原因，利用改进的PSO算法对非对称主轴结构参数进行多工况优化设计。优化结果表明，当Dl=68.5 mm，Dr=55.1 mm时，主轴左右两侧处于等强度设计，此时系统整体的扭转振动响应达到一个最优的平衡位置，右侧波动扭矩幅值下降25%左右，右侧扭转振动能量下降15%左右，优化后的主轴系统载荷分配更加合理，优化效果明显。

(3) 本文研究工作揭示了某特种车辆传动主轴在一挡工况下发生右端断裂失效的根本原因，针对断裂原因基于改进的PSO算法对主轴结构进行了优化，所得优化结果对于类似非对称非等强度设计结构的结构优化具有指导意义。