桑 萌， 黄 康， 熊杨寿

(1.合肥工业大学 机械工程学院，合肥 230009；2.安徽省数字设计与制造重点实验室，合肥 230009)

行星齿轮传动因其结构紧凑、传动比大、传动效率高等特点，适用于高速大功率以及低速大扭矩的传动场景，在运行过程中常承受较为复杂的动态载荷，相比于一般机械设备更易发生故障。因此，针对行星齿轮箱的故障诊断研究具有重要意义[1-2]。

3K行星齿轮传动作为一种基本行星传动类型，与2K-H行星轮系相比，结构更为紧凑，能提供更大传动比范围和更高的承载能力，在航空航天、工业机器人、微型传动等领域有着广泛的应用前景[3-4]。目前关于行星齿轮传动的故障诊断研究主要集中在2K-H行星轮系上，对于3K行星轮系的研究还较少。

关于行星齿轮传动的故障诊断方法，其中较为常用和成熟的为基于振动传感器的箱体振动信号分析[5-6]。但由于行星轮随行星架绕中心轴的公转，行星齿轮箱中故障激励点到测量点传递路径时变，振动信号会产生调制现象，给故障诊断带来困难[7]。此后部分学者，针对箱体振动信号的调制问题，提出利用齿轴扭转振动信号的故障诊断方法[8-9]。但该方法同样依赖角度编码器等机械检测设备，存在易受外部噪声干扰和安装限制较大等缺点。

为克服上述振动信号缺陷，近年来，电机电流信号分析(motor current signature analysis，MCSA) 逐渐被用于齿轮系统故障诊断[10]。如文献[11]通过监测风力发电机定子电流信号来检测发电机组传动系统故障；文献[12]提出异步电动机电流信号的调制信号双谱分析方法，实现对齿轮传动系统的故障监测；文献[13]通过定子电流信号分析了时变工况下行星齿轮箱故障。MCSA以电机定子电流为切入点，依据电机定子电流的波动反映扭矩的波动，从而对设备进行故障诊断。在实际操作时，使电机定子电源线 (三相中任一相) 通过霍尔电流传感器或电流钳，利用霍尔磁平衡原理，将电流信号经采集卡转化为可以识别的信号，属于无创检测，具有容易获得、抗干扰能力强等特点[14]。但目前关于MCSA的理论分析大多基于数学模型的推导，将转矩波动设定为具有固定幅值和故障频率的简谐形式，并不能反映如时变啮合刚度和阻尼等齿轮系统内部参数的影响[15-16]。

鉴于此，本文以电驱动3K-II型行星传动系统为研究对象，综合考虑齿轮传动部分时变啮合刚度、阻尼、齿侧间隙、综合啮合误差等非线性因素和电机-减速器系统内部机-电-磁耦合效应，建立了系统的机电耦合动力学模型。并以此为基础，通过去工频边带频谱和包络解调谱对不同齿轮存在裂纹故障时的定子电流信号进行了分析，探讨了齿轮裂纹故障对电流信号的调制作用。分析结果表明本文提出的机电耦合模型的合理性，能够准确模拟电机-3K减速器系统电流信号中的齿轮裂纹故障特征，可为研究各类齿轮故障情况下的电机定子电流信号和相关故障诊断技术提供理论指导。

1 3K型行星传动系统机电耦合动力学建模

1.1 齿轮系统动力学模型

以3K-II型行星轮系为研究对象，建立如图1所示，包含电动机 (转子)及负载的平移-扭转耦合集中参数模型。其中，下标c,b,e,s,pn,M,L(n=1,2,…,N,N为行星轮个数)分别为行星架、齿圈b、齿圈e、太阳轮、行星轮、电动机(转子)和负载。下标spn，bpn，epn分别为太阳轮-行星轮、齿圈b-行星轮、齿圈e-行星轮齿轮副。定义j=c，b，e，s，pn，M，L；g=spn，bpn，epn。本文如无特殊说明，均沿用此下标。kjx，kjy分别为各构件沿x方向、y方向的支撑刚度；cjx，cjy分别为各构件沿x方向、y方向的支撑阻尼系数；kju为各构件扭转支撑刚度；cju为各构件扭转支撑阻尼系数；kcMu,cLu分别为与电动机(转子) 和负载相连的联轴器扭转刚度；ccMu,cLu为与电动机(转子) 和负载相连的联轴器扭转阻尼系数；kg为各齿轮副时变啮合刚度，cg为各齿轮副啮合阻尼系数，eg为各齿轮副综合啮合误差。令齿圈b固定，太阳轮为输入构件，齿圈e为输出构件，即电动机输出转子与太阳轮连接，外部负载与齿圈e连接。采用牛顿运动定律建立系统运动微分方程如下。

图1 3K-II型行星轮系平移-扭转耦合动力学模型Fig.1 Translational-torsional coupling dynamics model of the 3K-II planetary gear system

电机转子

(1)

太阳轮

(2)

ψsn=π/2+ψn-αs,ψn=2(n-1)π/n,

n=1,2,…,N,

δspn=(xs-xpn)cos(ψsn)+(ys-ypn)sin(ψsn)+

rsbθs-(-rpbθpn)-rcbθccosαs+espn(t)

齿圈b

(3)

ψbn=π/2-(ψn+αb),

δbpn=-(xb-xpn)cos(ψbn)+(yb-ypn)sin(ψbn)+

rbbθb-rpbθpn-rcbθccosαb+ebpn(t)

齿圈e

(4)

ψen=π/2-(ψn-αe)

δepn=(xe-xpn)cos(ψen)-(ye-ypn)sin(ψen)-

rebθe+rpbθpn+rcbθccosαe+eepn(t)

行星架

(5)

行星轮

(6)

负载

(7)

(8)

考虑到动态啮合过程中，各构件转速往往围绕某一值波动，而非保持不变，通过构件转速来表示时间和误差项的对应关系不够恰当。本文使用各构件相对行星架的转角位移来表示各误差激励项。

eg(θ)=Egsin(Z1θ1c+φg)

(9)

式中:Eg为各齿轮副综合啮合误差幅值；Z1和θ1c分别为各齿轮副主动轮的齿数及其相对于行星架转角；φg为各齿轮副综合啮合误差初始相位。

同理，通过势能法[17-20]计算的各齿轮副时变啮合刚度也可表示为随转角变化的形式

(10)

(11)

式中：上标i为第i个啮合齿对；Ntooth为同时参与啮合的齿对数；下标1，2分别为主、从齿轮；假定主动轮裂纹故障出现在第一个啮合齿对初始位置，取w=mod(θ1c-2π(i-1)/Z1,Z1θT)，s=mod(w,ceil(ε)θT)；其中θT=2π/Z1为主动轮相对于行星架的啮合角周期；mod为取余函数；ceil为向上取整函数；ε为齿轮副重合度。

将式(1)～式(7)进行整理，可得系统运动方程的矩阵形式如下

T(t)+F(t)

(12)

式中：M，G，Cb，Cm，Kb，Km，KΩ分别为广义质量矩阵、陀螺矩阵、支撑阻尼矩阵、啮合阻尼矩阵、支撑刚度矩阵、啮合刚度矩阵和向心刚度矩阵；T(t)和F(t)分别为外部激励力向量和内部激励力向量；q为位移向量。

1.2 永磁同步电机数学模型

永磁同步电机因其效率高、功率密度高等特点，在现代工业生产中被广泛应用。本文拟以永磁同步电机驱动的电机-减速器系统为对象建立3K型行星轮系机电耦合动力学模型。

永磁同步电机物理模型示意图，如图2所示。包括固定的三相绕组磁轴AM，BM，CM(定子) 和以电角速度ωMr转动的永磁体S-N(转子)。当定子绕组中接入三相对称交流电后，将产生旋转磁场，带动磁体转子同步旋转，用以实现电机输出轴机械转动。

图2 永磁同步电机的物理模型示意图Fig.2 Diagram of the physical model of a permanent magnet synchronous motor

为降低模型复杂程度，便于求解，本文采用Park变换，将三相绕组用相互独立的两项正交对称绕组等效替代，在dM，qM坐标系(图2中浅色磁轴) 中建立永磁同步电机数学模型[21-22]。相应的三相-两相坐标变换公式、电压方程、磁链方程、电磁扭矩计算公式分别如式(13)～式(16)所示。

(13)

(14)

(15)

TM=1.5pM(ΨMfIqM+(LdM-LqM)IdMIqM)

(16)

1.3 机电耦合动力学模型

图3 3K型行星传动系统机电耦合关系示意图Fig.3 Diagram of electromechanical coupling relationship of 3K planetary transmission system

(17)

式中:第一个方程为3K行星传动系统机械部分受迫振动模型;第二个方程为电气系统等效电路模型;U=[VdM,VqM]T为电压向量;I=[IdM,IqM]T为电流向量;Ψ=[ΨdM,ΨqM]T为定子磁链向量;ΨMf=[0,ωMrΨMf]为转子磁链向量。

2 齿轮系统故障对定子电流的影响

2.1 定子电流信号AM-FM数学模型

本节主要从数学公式推导角度讨论由故障引起的调制效应对定子电流信号的影响。

假设齿轮故障引起的转矩波动以故障特征频率ffault的正弦形式引入，则电机总体负载转矩可以表示如下

Ts(t)=Ts0+Tsfsin(ωfaultt+φfault)

(18)

式中：Ts0为瞬时平均负载转矩；Tsf为故障引入转矩波动幅值；φfault为故障初始相位；ωfault=2πffault为故障角频率，其中ffault与各齿轮存在裂纹时的故障特征频率fscf,fpcf,fbcf,fecf对应。

假定电机-减速器系统稳态条件下，电磁转矩TM=Ts0，则电机转子转速可表示为

(19)

式中，ωMr0为电机转子平均转速。

对转子转速ωMr(t)积分，可得转子转角位移

(20)

对于永磁同步电机，在定子坐标系中的转子磁动势可表示为

FMr,(Ms)(θMs,t)=FMrsin{pM[θMs-θMr(t)]}=

(21)

式中：FMr为转子磁动势幅值；ωps=pMωMr0为供电频率；pM为磁极对数；θMs为定子坐标系中的转子角位置。

而定子坐标系中的定子磁动势与齿轮故障引起的转矩波动无关[23]，可以写为

FMs,(Ms)(θMs,t)=FMssin(pMθMs-ωpst)

(22)

式中，FMs为定子磁动势幅值。

气隙磁通密度可表示为总体磁动势和气隙磁导率之积

B(θMs,t)=Λ[FMs,(Ms)(θMs,t)+FMr,(Ms)(θMs,t)]=

BMssin(pMθMs-ωpst)+

(23)

式中，气隙磁导率Λ被视为常数。通过对B(θMs,t)关于转子转角θMs积分可得定子任意项磁通量

ΨMssin(ωpst)

(24)

式中，BMs, Mr和ΨMs,Mr分别为由定子和转子造成的气隙磁通密度和磁通量幅值。

又，对于定子任意相电压有

V(t)=RMsI(t)+dΨ(t)/dt

(25)

结合式(24)、式(25)，定子电流I(t)可表示为

Is0sin(ωpst)-IMsωpscos(ωpst)-

(26)

式中：Is0为流过定子电阻的电流幅值；IMs=ΨMs/RMs；IMr=ΨMr/RMs。

由上述数学推导过程可知，齿轮故障对定子电流信号同时有调频和调幅效果，载波频率为电源供电频率ωps，调制频率为齿轮故障特征频率ωfault。这为通过定子电流信号进行电机-减速器系统故障诊断提供了依据。

2.2 定子电流信号故障特征分析

依据建立的定子电流AF-FM模型，分别从调频和调幅角度，分析电机-减速器系统齿轮故障在频谱和包络谱图中的表现特征。

2.2.1 频谱分析

根据雅可比-安格尔恒等式

(27)

可将调频函数写成

(28)

式中，Jm(z)为m阶第一类Bessel函数。则式(26)中的定子电流可表示为

(29)

对式(29)进行傅里叶变换可得

exp[-jarctan(IMsωps/Is0)]-

δ[ω-ωps]exp[jarctan(IMsωps/Is0)]}-

exp(-jmφfault)+δ[ω-(ωps+mωfault)]·

{δ{ω-[(ωps+(m+1)ωfault)]}exp[j(m+1)φfault]+

δ{ω+[(ωps+(m+1)ωfault)]}exp[-j(m+1)φfault]}-

[(ωps+(m-1)ωfault)]}exp[j(m-1)φfault]+

δ{ω+[(ωps+(m-1)ωfault)]}·

exp[-j(m-1)φfault]}

(30)

式中，δ(·)为单位脉冲函数。

由式(30)可知，定子电流信号将会在频率±(ωps±mωfault)处，即正负供电频率±ωps附近，出现边带，边带间隔为故障频率ωfault。如考虑供电频率的高阶谐波，边带可能会在±(lωps±mωfault)处出现，其中m，l为正整数。

2.2.2 包络谱分析

根据式(26)中定子电流表达式的第二项可知，其振幅受故障频率调制作用。对第二项幅值进行傅里叶变换可得包络谱

δ(ω-ωfault)exp(jφfault)]+IMrωpsδ(ω)

(31)

由式 (31) 可知，在频率±ωfault处包络谱会出现峰值。如考虑高阶谐波，峰值也会在±mωfault处出现，其中m为正整数。

3 基于3K行星轮系机电耦合模型的定子电流仿真信号分析

3.1 基于机电耦合模型的定子电流信号仿真

本文选定的3K行星传动机电耦合系统中减速器部分和电机部分主要参数,分别如表1和表2所示。

表1 机电耦合系统减速器部分主要参数Tab.1 Main parameters of the reducer part of the electromechanical coupling system

表2 机电耦合系统电机部分主要参数Tab.2 Main parameters of the motor part of the electromechanical coupling system

考虑到太阳轮与齿圈e在行星轮同齿侧啮合，齿圈b在行星轮另一侧啮合的情况，为模拟局部齿根裂纹故障影响，本文针对5种齿轮裂纹故障情况进行分析，即：①太阳轮齿根裂纹；②太阳轮-齿圈e侧行星轮齿根裂纹；③齿圈b侧行星轮齿根裂纹；④齿圈b齿根裂纹；⑤齿圈e齿根裂纹。本文主要为探讨基于电驱动3K行星传动机电耦合模型，采用电机定子电流信号进行齿轮裂纹故障定位的可行性。基于这一研究目的，假定各齿轮裂纹均为直线且长度(1.32 mm) 和角度(65°)一致，上述 5种裂纹故障情况下受故障影响的各齿轮副啮合刚度,如图4所示。

(a) 太阳轮齿根裂纹

(b) 太阳轮-齿圈e侧行星轮齿根裂纹

(c) 齿圈b侧行星轮齿根裂纹

(d) 齿圈b齿根裂纹

(e) 齿圈e齿根裂纹图4 受裂纹故障影响的啮合刚度曲线Fig.4 Meshing stiffness curves affected by crack faults

(a) 正常情况

(b) 太阳轮故障

(c) s-e侧行星轮故障

(d) b侧行星轮故障

(e) 齿圈b故障

(f) 齿圈e故障图5 定子电流信号时域波形Fig.5 Time domain waveforms of stator current signals

表3 各特征频率Tab.3 Each characteristic frequency Hz

从图5所示的6种齿轮健康状况下的定子电流信号时域波形中可以看出，正常情况和各类齿轮裂纹故障情况下的波形区别并不明显，供电频率占绝对优势，表现为简谐形式。

为提取定子电流信号中的故障特征，本节依据前文所述，分别对信号频域边带和包络谱进行分析。

3.2 频域边带分析

由图5可知，定子电流信号主要表现为以供电频率为主频率的简谐形式，直接对电流信号进行傅里叶变换，频谱主要为供电频率处的脉冲信号，使得其余特征频率被淹没在其中。因此，要分析电流信号频域边带特征需对其进行去工频处理。本文采用时移叠加的方式来去除工频，该方法易操作，便于工业应用。

3.2.1 时移叠加

时移叠加即将电流信号与其平移半个周期后的自身相加，以消除基频信号的方法，根据上文建立的电流信号AM-FM(amplitude modulation-frequency modulation)模型，取式(29)所示电流信号，时延半个周期可得

(32)

(33)

式中，Tps=2π/ωps为供电周期。

由式(32)、式(33)可以看出，时移叠加后的电流信号能有效去除电流中的基频，而保留由故障引起的调制成分。

3.2.2 去工频后电流频谱

去工频后的电流残余信号傅里叶频谱，如图6所示。从图6可以看出，去工频后的电流傅里叶谱图中消去了主导的供电频率 (fps=50 Hz)，而保留了边带。当齿轮箱存在故障时，会产生新的边带成分，这些边带与故障频率相关。

(a) 太阳轮故障

(b) s-e侧行星轮故障

(c) b侧行星轮故障

(d) 齿圈b故障

(e) 齿圈e故障图6 电流信号去除工频后的傅里叶谱Fig.6 Fourier spectrum of current signal after removing power supply frequency

由图6(a)可知，新出现的边带成分对应正负供电频率加减太阳轮故障频率的倍频，即±fps±mfscf，其中m为正整数，如170 Hz=fps+fscf，70 Hz=-fps+fscf，这与2.2节中关于定子电流信号的频谱分析结果对应。这种现象意味着太阳轮存在故障。

在图6(b)和图6(c)中，新出现的边带成分与正负供电频率加减行星轮故障频率的倍频相对应，即±fps±mfpcf，如207 Hz=fps+6fpcf，107 Hz=-fps+6fpcf，这一特征表示行星轮存在故障。

在图6(d)所示的齿圈b故障对应的频谱中，在正负供电频率加减齿圈b故障频率的倍频处出现新的边带成分，即±fps±mfbcf，如200.5 Hz=fps+5fpcf，70.5 Hz=-fps+4fpcf。

而当齿圈e出现故障时，相应的在正负供电频率加减齿圈e故障频率的倍频，即±fps±mfecf处，出现新的边带，如222 Hz=fps+6fecf，93 Hz=-fps+4fecf。

通过对比5种齿轮故障情况可知，由太阳轮故障引起的边带成分虽然也能识别，但与其他几种故障相比不够明显。这与3K型行星轮系传动比较大而导致的输入端传递扭矩较小，且3K-II型的一体式双联行星排结构使得齿轮对s-p(n)的有效啮合齿宽相比于其他齿轮对较大相关。此外，值得注意的是，在各齿轮故障引起的边带成分中包括供电频率fps的负值，这一现象与Chen等研究中的结果存在差别，这可能与故障频率和供电频率相近甚至超过供电频率有关。

3.3 包络谱分析

通过2.2节中的分析可知，频谱边带中包含故障频率和供电频率的组合成分，并不能直观显示出故障特征。而包络谱中峰值与故障特征直接相关，本节通过采用Hilbert变换幅值解调的方式获得包络谱。

各齿轮裂纹故障下的包络谱，如图7所示。与正常情况相比，在齿轮故障情况下，包络谱中会出现新的峰值，这些峰值位置与故障频率对应。

(a) 太阳轮故障

(b) s-e侧行星轮故障

(c) b侧行星轮故障

(d) 齿圈b故障

(e) 齿圈e故障图7 电流信号包络谱Fig.7 Envelope spectrum of current signal

对于太阳轮故障，从图7(a)中可以看出，在频率120 Hz处(fscf)包络谱幅值有所增加，这与太阳轮故障频率一致。与其他齿轮故障相比，图7(a)中太阳轮故障特征不够明显，这一现象同上文频率边带分析对应。

在图7(b)和图7(c)中，可以看出在行星轮故障频率及其倍频处相比于正常情况幅值明显增加，如131 Hz=5fpcf，157 Hz=6fpcf，这表明行星轮存在故障。

当齿圈b存在故障时，可以从图7(d)中看出在齿圈b故障频率及其倍频处出现明显的峰值，如60.5 Hz=2fbcf，180.5 Hz=6fbcf。这一现象说明了齿圈b存在故障。

而对于齿圈e故障，与正常情况相比，图7(e)中出现了新的峰值成分，这些峰值位置与齿圈e故障频率及其倍频对应(28.5 Hz=fecf，172 Hz=6fecf)。这同样说明了齿圈e存在故障。

通过上述分析，说明电机定子电流信号能够准确反映3K行星齿轮箱的故障特征信息，电流信号同时受故障特征的幅值调制和频率调制作用。也说明了本文提出的机电耦合模型的合理性，能够准确模拟电机-3K减速器系统电流信号中的齿轮故障特征，为研究各类齿轮故障情况下的电机定子电流信号提供理论指导。

4 结 论

本文通过综合考虑齿轮传动部分时变啮合刚度、齿侧间隙、综合啮合误差等非线性因素和电机-减速器系统内部机-电-磁耦合效应，建立了电驱动3K行星传动系统的机电耦合动力学模型。并基于该模型，计算了5种裂纹故障情况下的时变啮合刚度，得到了不同齿轮裂纹故障情况下的定子电流仿真信号。通过时移叠加去工频频谱和包络解调谱对电流信号进行分析，得到下列结论：

(1) 定子电流信号受故障特征的频域调制作用。在通过时移叠加去工频后的频谱中，齿轮故障表现为围绕供电频率及供电频率负值的变频带，边带间隔为故障频率。且由于3K行星轮系大速比特点，各齿轮故障频率与供电频率接近或超过供电频率，在频率正轴出现的边带成分中包含供电频率负值。

(2) 定子电流信号受故障特征的幅值调制作用。在包络解调谱中，齿轮故障表现为在故障频率及其倍频处出现峰值。

(3) 由于3K行星轮系传动比较大，电机输入端扭矩较小，且3K-II型特殊的一体式双联行星排结构令齿轮对s-p(n)的有效啮合齿宽相比于其他齿轮副较大。太阳轮故障在定子电流频谱边带和包络谱中虽能识别，但相较于其他齿轮故障，表现不够明显。

(4) 通过机电耦合模型得到的电机定子电流信号分析结果与通过AM-FM数学模型分析结果吻合，表明本文提出的机电耦合模型的合理性，能够准确模拟电机-3K减速器系统电流信号中的齿轮故障特征，为研究各类齿轮故障情况下的电机定子电流信号提供理论指导。