引言

旋转机械是国防及工业等领域非常重要的动力装置，但是转子的振动异常会频繁造成机器报警或停车故障，严重影响生产率。转子振动的主要原因是质量不平衡^[1]，所以为解决这一问题，需要了解不平衡的来源以及具体的参数。

目前进行不平衡识别的方法主要分为两类：一是基于模型的识别，如模态扩展^[2]、模态分析和傅里叶变换^[3-4]，这种方法对转子有限元模型的精确性要求较高，而且需要准确获取转子模态信息，工作量很大；二是基于算法的识别，即利用启发式^[5]、最小角回归^[6]、伪随机优化^[7]、BP神经网络^[8]、稳健回归分析^[9]等算法进行寻优，这些算法都可用于识别不平衡参数，但存在设置参数多、容易陷入局部最优等局限。

果蝇算法是潘文超^[10]提出的一种全局寻优算法，具有程序简单、设置参数少的优点，已被广泛应用于多个领域，但目前还未见将果蝇算法用于转子不平衡参数识别的文献。本文提出一种基于反问题的等效不平衡力识别理论，结合改进的果蝇优化算法来识别转子的不平衡参数，以期为旋转机械的故障诊断提供一种更为有效的方法。

1 转子不平衡识别与优化方法 1.1 转子动力学及不平衡识别理论

采用有限元方法建立转子运动微分方程为

$ \mathit{\boldsymbol{M}}{{\mathit{\boldsymbol{\ddot r}}}_0}\left( t \right) + \mathit{\boldsymbol{D}}{{\mathit{\boldsymbol{\dot r}}}_0}\left( t \right) + \mathit{\boldsymbol{K}}{\mathit{\boldsymbol{r}}_0}\left( t \right) = {\mathit{\boldsymbol{F}}_0}\left( t \right) $

(1)

其中，M、D、K分别是系统的质量矩阵、包含陀螺效应的阻尼矩阵和刚度矩阵，r₀(t)为转子的位移向量，F₀(t)为广义力向量。在式(1)中加上不平衡量，转子的运动方程可由式(2)描述

$ \mathit{\boldsymbol{M\ddot r}}\left( t \right) + \mathit{\boldsymbol{D\dot r}}\left( t \right) + \mathit{\boldsymbol{Kr}}\left( t \right) = {\mathit{\boldsymbol{F}}_0}\left( t \right) + \Delta \mathit{\boldsymbol{F}}\left( t \right) $

(2)

其中ΔF(t)为转子不平衡故障产生的等效不平衡力向量。在式(1)、(2)中，滑动轴承油膜力采用线性化表达式，则该转子系统可视为线性系统；同时假定同一转速下不平衡激励前后的轴承刚度和阻尼近似不变，只考虑不平衡力的作用，则将式(2)减去式(1)可得

$ \mathit{\boldsymbol{M}}\Delta \mathit{\boldsymbol{\ddot r}}\left( t \right) + \mathit{\boldsymbol{D}}\Delta \mathit{\boldsymbol{\dot r}}\left( t \right) + \mathit{\boldsymbol{K}}\Delta \mathit{\boldsymbol{r}}\left( t \right) = \Delta \mathit{\boldsymbol{F}}\left( t \right) $

(3)

其中，Δr(t)=r(t)－r₀(t)，Δ${\mathit{\boldsymbol{\dot r}}} $(t)= ${\mathit{\boldsymbol{\dot r}}} $ (t)－${\mathit{\boldsymbol{\dot r}}} $₀(t)，Δ $ {\mathit{\boldsymbol{\ddot r}}}$ (t)=$ {\mathit{\boldsymbol{\ddot r}}}$(t)－$ {\mathit{\boldsymbol{\ddot r}}}$₀(t)。

转子的位移Δr(t)和等效不平衡力ΔF(t)表示成复数形式：

$ \left\{ \begin{array}{l} \Delta \mathit{\boldsymbol{r}}\left( t \right) = {\mathit{\boldsymbol{R}}_0}{{\rm{e}}^{{\rm{i}}\left( {\mathit{\Omega }t + \alpha } \right)}}\\ \Delta \mathit{\boldsymbol{F}}\left( t \right) = \mathit{\boldsymbol{U}}{\mathit{\Omega }^2}{{\rm{e}}^{{\rm{i}}\left( {\mathit{\Omega }t + \alpha } \right)}} \end{array} \right. $

(4)

式中R₀、U分别为位移幅值向量和不平衡幅值向量，Ω为转速。将式(4)代入式(3)，则有

$ {\mathit{\boldsymbol{R}}_0} = {\left( { - \mathit{\boldsymbol{M}}{\mathit{\Omega }^2} + {\rm{i}}\mathit{\Omega }\mathit{\boldsymbol{D}} + \mathit{\boldsymbol{K}}} \right)^{ - 1}}\mathit{\boldsymbol{U}}{\mathit{\Omega }^2} $

(5)

记H=Ω²(－MΩ²+iΩD+K)^－1，则有

$ {\mathit{\boldsymbol{R}}_0} = \mathit{\boldsymbol{HU}} $

(6)

提取H中对应测点所在节点的行和不平衡所在节点的列构成矩阵B，提取U中对应的不平衡量所在的行并构成向量Q，则有

$ {\mathit{\boldsymbol{R}}_0} = \mathit{\boldsymbol{BQ}} $

(7)

采集转子的不平衡响应代入式(7)，反求Q，则有

$ \mathit{\boldsymbol{Q}} = \mathit{\boldsymbol{B}}{\left( \mathit{\Omega } \right)^{ - 1}}{\mathit{\boldsymbol{R}}_0} $

(8)

进而有ΔF(t)=QΩ²e^i(Ωt+α)，即实测等效不平衡力可被求解出来变成已知量。含有未知不平衡量的理论不平衡力如式(9)

$ \Delta \mathit{\boldsymbol{F}}\left( {\beta ,t} \right) = u{\mathit{\Omega }^2}\left[ \begin{array}{l} \cos \left( {\mathit{\Omega t} + \varphi } \right)\\ \sin \left( {\mathit{\Omega t} + \varphi } \right) \end{array} \right] $

(9)

其中β包含不平衡量的幅值u和相位角φ，通过理论不平衡力与实测等效不平衡力构建优化目标函数，再采用优化策略对目标函数进行优化，即可识别出不平衡参数β。

1.2 目标函数的建立

目标函数是不平衡的理论值(待估计)与估计的等效不平衡力之间差值的最小绝对值，即

$ \min \left| {\Delta \mathit{\boldsymbol{F}}\left( {\beta ,t} \right) - \Delta \mathit{\boldsymbol{F}}\left( t \right)} \right| $

(10)

理论不平衡力ΔF(β, t)将其表示成复数形式，则

$ \Delta \mathit{\boldsymbol{F}}\left( {\beta ,t} \right) = u{\mathit{\Omega }^2}{{\rm{e}}^{{\rm{i}}\left( {\mathit{\Omega }t + \alpha } \right)}} $

(11)

同理，求解出Δr(t)后将其代入式(3)，则可将ΔF(t)作为已知量用相同的方式表示出来。假设实测不平衡力幅值大小为|F|，相位角为α，则

$ \Delta \mathit{\boldsymbol{F}}\left( t \right) = \left| \mathit{\boldsymbol{F}} \right|{{\rm{e}}^{{\rm{i}}\left( {\mathit{\Omega }t + \alpha } \right)}} $

(12)

为消去时间这个无关变量，只保留待识别变量，将式(11)除以式(12)并取自然对数，则式(9)变为

$ \min \left| {\ln \left( {\frac{{u{\mathit{\Omega }^2}}}{{\left| \mathit{\boldsymbol{F}} \right|}}} \right) + {\rm{i}}\left( {\varphi - \alpha } \right)} \right| $

(13)

目标函数是未知不平衡参量u和φ的二元函数，在可行域内运用多种优化策略对目标函数进行优化，当目标函数收敛到最小值时即可识别出不平衡参数。本文主要通过果蝇算法(FOA)优化目标函数，并与模拟退火算法(SA)、遗传算法(GA)的识别结果进行对比。

1.3 果蝇优化算法

果蝇优化算法是一种基于果蝇觅食行为推演出的寻找全局最优解的算法。果蝇利用自身敏锐的嗅觉搜集空气中的各种气味，利用视觉发现食物或同伴聚集的位置，并向该方向飞去搜寻食物。本文针对原始果蝇算法步长固定的局限性，提出分阶段变步长的改进措施，分别设置果蝇移动范围参数a、b, 增长因子c和缩减因子d。改进的果蝇算法寻优步骤如下(主要改进在第8、9步)：

1) 随机初始化果蝇群体位置；

2) 赋予果蝇个体利用嗅觉搜寻食物的随机距离与方向；

3) 计算味道浓度判定值，此值为距离的倒数；

4) 将味道浓度判定值代入味道浓度判定函数以求出该果蝇个体位置的味道浓度；

5) 找出该果蝇群体中味道浓度最高的果蝇(求极小值)；

6) 保留最佳味道浓度值与对应的果蝇的坐标值，此时果蝇群体利用视觉往该位置飞去；

7) 进入迭代寻优，重复执行步骤2)~5)，并判断味道浓度是否优于前一迭代味道浓度，若是则实行步骤6)；

8) 进入第一阶段跌代寻优，设置增长因子c>1，改变果蝇个体的初始步长，更新个体位置重复步骤2)~5)，判断味道浓度是否优于前一迭代味道浓度，若是则实行步骤6)，直至达到迭代次数的一半；

9) 进入第二阶段寻优，设置缩减因子d < 1，改变果蝇个体的初始步长，更新个体位置重复步骤2)~5)，判断味道浓度是否优于前一迭代味道浓度，若是则实行步骤6)，直至达到最大迭代次数，结束。

2 仿真分析 2.1 实验台动力学建模

建立转子实验台的有限元模型如图 1所示，转轴共划分为13个节点、12个单元，位移传感器测点分别位于第3、10节点，轮盘位于第7节点。通过仿真计算转子一阶临界转速为2800r/min，实测一阶临界转速为2820r/min，表明所建有限元模型是有效的，接近实际转子系统。

2.2 模拟及结果分析

不平衡参数具体识别过程为：首先计算转子未加不平衡量时6000r/min工作转速下的初始位移，然后在轮盘上施加如表 1所示的不平衡量，计算出转子在相同转速下的位移后代入式(3)求解等效不平衡力，再由等效不平衡力与理论不平衡力建立目标函数，分别运用FOA、SA、和GA 3种算法进行目标函数优化。迭代寻优过程如图 2所示，识别结果如表 2所示。

图 2的迭代过程表明，遗传算法的总迭代次数最少，模拟退火算法最多；而在使函数值下降的有效迭代次数中，果蝇算法的有效迭代次数最少，目标函数值下降的速度最快，因此在实际运行中其迭代时间最短，效率最高。

表 2的识别结果表明，果蝇算法对不平衡幅值的识别误差最大为0.23%，相位识别最大误差为0.24%，均小于其他两种算法，证明了果蝇算法识别性能的优越性。

3 实验验证 3.1 实验台

本文搭建滑动轴承支撑的单盘转子实验台如图 3所示，实验台主要参数见表 3。转轴通过波纹管联轴器与额定功率为0.3kW的直流电机相联，最高转速达到12000r/min，3只电涡流位移传感器距离两端轴承75mm，速度和加速度传感器各一只测量基座的振动。实验台配套有数据采集系统及相应的分析软件，设置采样频率为1024Hz，采集时间为2s，对采集到的时域响应数据进行正弦曲线拟合，得到不平衡响应的幅值和初相位，用来估算不平衡力。

3.2 实验结果分析

利用2.2节的步骤进行不平衡识别和分析，识别结果如表 4所示，3种优化算法对目标函数的迭代过程如图 4所示。

表 4的识别结果表明，果蝇算法对不平衡量的幅值、相位识别的最大误差分别为-4.66%、11.07%，均小于其余两种算法。图 4的迭代过程表明果蝇算法的有效迭代次数最少，实际迭代时间也最短。综合比较，果蝇算法的识别精度和识别效率最优。

4 结束语

本文通过建立转子-轴承系统有限元模型，构建了关于不平衡参数寻优的目标函数，并基于改进的果蝇算法对目标函数进行优化，克服了基本果蝇算法迭代步长固定的缺陷。仿真和实验结果表明果蝇算法能准确识别出不平衡参数，识别效果最佳。果蝇算法的缺点是算法的稳定性差，容易陷入局部最优，后续的研究内容将考虑对果蝇算法的进一步改进。