Motor Current Signature Analysis for Gear Train Fault Detection

This example shows how to identify faults in a gear train using motor current signature analysis (MCSA) of a current signal driving a hobby-grade servo. MCSA is a useful method for the diagnosis of faults that induce torque or speed fluctuations, and has been proven to be ideal for motor fault analysis. Gear fault detection using traditional vibration instruments is challenging, especially in cases where the gear train is not easily accessible for instrumentation with accelerometers or other vibration sensors, like the inner workings of a nuclear power plant. This example illustrates how to apply current signature analysis to extract spectral metrics to detect faults in drive gears of a hobby-grade servo motor. The simplified workflow to obtain the spectral metrics from the current signal is as follows:

  1. Compute nominal rpm to detect frequencies of interest.

  2. Construct frequency bands where fault signals may be present.

  3. Extract power spectral density (PSD) data.

  4. Compute spectral metrics at the frequency bands of interest.

Hardware Overview

For this example, the electrical current data was collected from a standard Futaba S3003 hobby servo, which was modified for continuous rotation. Servos convert the high speed of the internal DC motor to high torque at the output spline. To achieve this, servos consist of a DC motor, a set of nylon or metal drive gears, and the control circuit. The control circuit was removed to allow the current signal to the DC motor to be directly monitored. The tachometer signal at the output spline of the servo was collected using an infrared photointerrupter along with a 35 mm diameter, eight-holed, standard hobby servo wheel. The eight holes in the wheel were equally spaced and the IR photointerrupter was placed such that there were exactly eight pulses per rotation of the servo wheel horn. The servo and photointerrupter were held in place by custom 3-D printed mounts.

The DC motor was driven at a constant 5 volts, and with four pairs of gears providing 278:1 speed reduction, the shaft speed at the spline was about 50 rpm. The current consumption was calculated using Ohm's law by measuring the voltage drop across a 0.5 Ohm resistor. Since the change in current measurement values was too small for detection, the current signal was amplified using a single-supply sensor interface amplifier. The amplified current signal was then filtered using an anti-aliasing fifth-order elliptic low-pass filter to smooth it and to eliminate noise before sending it to an Arduino Uno through an analog-to-digital converter (ADC).

As the flowchart shows, the current signal was first amplified and filtered using the amplifier and the anti-aliasing low-pass filter, respectively. The Arduino Uno sampled the current signal through an ADC at 1.5 kHz and streamed it to the computer along with the tachometer pulses as serial data at a baud rate of 115200. A MATLAB script fetched the serial data from the Arduino Uno and wrote it to a comma-separated values (CSV) file. The CSV files were then read and processed using MATLAB to extract the spectral metrics.

Servo Gear Train

The Futaba S3003 servo consists of four pairs of nylon gears as illustrated in the above figure. The pinion P1 on the DC motor shaft meshes with the stepped gear G1. The pinion P2 is a molded part of the stepped gear G1 and meshes with the stepped gear G2. The pinion P3, which is a molded part of gear G2, meshes with the stepped gear G3. Pinion P4, which is molded with G3, meshes with the final gear G4 that is attached to the output spline. The stepped gear sets G1 and P2, G2 and P3, and G3 and P4 are free spinning gears - that is, they are not attached to their respective shafts. The set of drive gears provides a 278:1 reduction going from a motor speed of 13901 rpm to about 50 rpm at the output spline when the motor is driven at 5 volts. The following table outlines the tooth count and theoretical values of output speed, gear mesh frequencies, and cumulative gear reduction at each gear mesh.

A total of 10 healthy data sets were collected before the faults were introduced in the stepped gears G2 and G3. Since the gears were molded out of nylon, simulated cracks were introduced in both the gears by cutting slots in the tooth space with a hobby knife. The tooth space is the gap between two adjacent teeth measured along the pitch circle of the spur gear. The slot depths were about 70 percent of the gear radius. A total of 10 faulty data sets were recorded after introducing the faults in the gears G2 and G3.

Visualize Data

The file mcsaData.mat contains servoData, a 10-by-2 table of timetables where each timetable corresponds to one data set. The first column of servoData contains 10 timetables of healthy data while the second column contains 10 timetables of faulty data. Each data set contains about 11 seconds of data sampled at 1500 Hz.

Load the data.

load('mcsaData.mat','servoData')
servoData
servoData=10×2 table
        healthyData            faultyData     
    ___________________    ___________________

    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}
    {16384x2 timetable}    {16384x2 timetable}

head(servoData.healthyData{1,1})
ans=8×3 timetable
         Time         Pulse    Signal
    ______________    _____    ______

    0 sec               0      66.523
    0.00066667 sec      0      62.798
    0.0013333 sec       0      63.596
    0.002 sec           0      64.128
    0.0026667 sec       0      60.669
    0.0033333 sec       0      62.798
    0.004 sec           0      65.459
    0.0046667 sec       0      56.678

Each timetable in servoData contains one data set with the tachometer signal in the first column and the current signal in the second column.

For this example, consider one set each of healthy and faulty data, and visualize the current signal and tachometer pulses.

healthyData = servoData.healthyData{1,1};
faultyData = servoData.faultyData{1,1};
healthyCurrent = healthyData.Signal;
faultyCurrent = faultyData.Signal;
healthyTacho = healthyData.Pulse;
faultyTacho = faultyData.Pulse;
tHealthy = healthyData.Time;
tFaulty = faultyData.Time;

figure
ax1 = subplot(221);
plot(tHealthy,healthyCurrent)
title('Current Signal - Healthy Gears')
ylabel('Current (mA)')
ax2 = subplot(222);
plot(tFaulty,faultyCurrent)
title('Current Signal - Faulty Gears')
ylabel('Current (mA)')
ax3 = subplot(223);
plot(tHealthy,healthyTacho)
title('Tachometer Pulse - Healthy Gears')
ylabel('Pulses, 8/rev')
ax4 = subplot(224);
plot(tFaulty,faultyTacho)
title('Tachometer Pulse - Faulty Gears')
ylabel('Pulses, 8/rev')
linkaxes([ax1,ax2,ax3,ax4],'x');
ax1.XLim = seconds([0 2]);

The output spline of the servo has eight pulses per rotation, corresponding to the eight holes on the servo wheel.

Compute Nominal Rpm

Compute the nominal speed to detect frequencies of interest in the gear system and match them correctly with the frequencies on the power spectra. Using the sampling frequency value of 1500 Hz, visualize the tachometer signal and compute the nominal rpm using tachorpm.

fs = 1500;
figure
tachorpm(healthyTacho,fs,'PulsesPerRev',8)

figure
tachorpm(faultyTacho,fs,'PulsesPerRev',8)

rpmHealthy = mean(tachorpm(healthyTacho,fs,'PulsesPerRev',8))
rpmHealthy = 49.8550
rpmFaulty = mean(tachorpm(faultyTacho,fs,'PulsesPerRev',8))
rpmFaulty = 49.5267

Observe that there is a very small difference in the output shaft speed between the healthy and faulty data sets. The actual nominal rpm values are also close to the theoretical value of 50 rpm. Hence, consider the same value of 49.9 rpm for both the healthy and faulty signal analysis.

Construct Frequency Bands

Constructing frequency bands is an important prerequisite for computing the spectral metrics. Using the tooth count of the drive gears in the gear train and the nominal rpm, first compute the frequencies of interest. The frequencies of interest are actual output speed values in hertz whose values are close to the theoretical values listed in the table below.

G4 = 41; G3 = 35; G2 = 50; G1 = 62;
P4 = 16; P3 = 10; P2 = 10; P1 = 10;
rpm = rpmHealthy;

FS5 = mean(rpm)/60;
FS4 = G4/P4*FS5;
FS3 = G3/P3*FS4;
FS2 =  G2/P2*FS3;
FS1 =  G1/P1*FS2;
FS = [FS1,FS2,FS3,FS4,FS5]
FS = 1×5

  231.0207   37.2614    7.4523    2.1292    0.8309

Next, construct the frequency bands for the all the output speeds which include the following frequencies of interest using faultBands.

  • FS1 at 231 Hz, its first two harmonics, and 0:1 sidebands of FS2

  • FS2 at 37.3 Hz, its first two harmonics, and 0:1 sidebands of FS3

  • FS3 at 7.5 Hz and its first two harmonics

  • FS4 at 2.1 Hz and its first two harmonics

[FB1,info1] = faultBands(FS1,1:2,FS2,0:1);
[FB2,info2] = faultBands(FS2,1:2,FS3,0:1);
[FB3,info3] = faultBands(FS3,1:2);
[FB4,info4] = faultBands(FS4,1:2);
FB = [FB1;FB2;FB3;FB4]
FB = 16×2

  187.9838  199.5348
  225.2452  236.7962
  262.5066  274.0577
  419.0045  430.5556
  456.2659  467.8170
  493.5273  505.0784
   28.8776   30.7407
   36.3299   38.1929
   43.7822   45.6452
   66.1390   68.0021
      ⋮

The output FB is a 16-by-2 array of frequency ranges for these frequencies of interest.

Extract Power Spectral Density (PSD) Data

Compute and visualize the power spectrum of the healthy and faulty data. Also plot the frequencies of interest on the spectrum plot by using the information in the structure info.

figure
pspectrum(healthyCurrent,fs,'FrequencyResolution',0.5)
hold on
pspectrum(faultyCurrent,fs,'FrequencyResolution',0.5)
hold on
info1.Labels = regexprep(info1.Labels,'F0','FS1');
info1.Labels = regexprep(info1.Labels,'F1','FS2');
helperPlotXLines(info1,[0.5 0.5 0.5])
info2.Labels = regexprep(info2.Labels,'F0','FS2');
info2.Labels = regexprep(info2.Labels,'F1','FS3');
helperPlotXLines(info2,[0.5 0.5 0.5])
info3.Labels = regexprep(info3.Labels,'F0','FS3');
helperPlotXLines(info3,[0.5 0.5 0.5])
info4.Labels = regexprep(info4.Labels,'F0','FS4');
helperPlotXLines(info4,[0.5 0.5 0.5])
legend('Healthy Data','Faulty Data','Location','South')
hold off

The plot in blue indicates the healthy spectrum while the plot in red indicates the spectrum of the faulty data. From the plot, observe the increase in amplitudes of fault frequencies:

  • 1FS1 at 231 Hz, its second harmonic 2FS1 at 462 Hz, and their respective sidebands

Zoom in on the plot to observe the increase in amplitudes of the following remaining frequencies:

  • 1FS2 at 37.2 Hz and its sidebands

  • 1FS3 at 7.5 Hz and its second harmonic 2FS3 at 15 Hz

  • 1FS4 at 2.1 Hz and its second harmonic 2FS4 at 4.2 Hz

Use pspectrum to compute and store the PSD and the frequency grid data for the healthy and faulty signals respectively.

[psdHealthy,fHealthy] = pspectrum(healthyCurrent,fs,'FrequencyResolution',0.5);
[psdFaulty,fFaulty] = pspectrum(faultyCurrent,fs,'FrequencyResolution',0.5);

Compute Spectral Metrics

Using the frequency bands and the PSD data, compute the spectral metrics for the healthy and faulty data using faultBandMetrics.

spectralMetrics = faultBandMetrics({[psdHealthy,fHealthy],[psdFaulty,fFaulty]},FB)
spectralMetrics=2×49 table
    PeakAmplitude1    PeakFrequency1    BandPower1    PeakAmplitude2    PeakFrequency2    BandPower2    PeakAmplitude3    PeakFrequency3    BandPower3    PeakAmplitude4    PeakFrequency4    BandPower4    PeakAmplitude5    PeakFrequency5    BandPower5    PeakAmplitude6    PeakFrequency6    BandPower6    PeakAmplitude7    PeakFrequency7    BandPower7    PeakAmplitude8    PeakFrequency8    BandPower8    PeakAmplitude9    PeakFrequency9    BandPower9    PeakAmplitude10    PeakFrequency10    BandPower10    PeakAmplitude11    PeakFrequency11    BandPower11    PeakAmplitude12    PeakFrequency12    BandPower12    PeakAmplitude13    PeakFrequency13    BandPower13    PeakAmplitude14    PeakFrequency14    BandPower14    PeakAmplitude15    PeakFrequency15    BandPower15    PeakAmplitude16    PeakFrequency16    BandPower16    TotalBandPower
    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    ______________

       0.002071           193.75         0.010883        0.50813            231.06         0.46652        0.0019579            272.5          0.01184       0.0020489           424.06         0.011224        0.54985               462          0.8955        0.0024297           493.69        0.0091051       0.0029648           29.812        0.0035477        0.015581           37.25          0.011131       0.0028858           44.688        0.0041316        0.011896            67.062          0.0072763        0.059124             74.5           0.033566         0.013218                82           0.007988         5.7905             7.4375           2.3115          0.068452            14.938          0.027653          0.79002             2.125           0.14381         0.098478             4.25           0.018058          3.9737    
      0.0098034           192.44         0.017915         3.6923            229.44          2.9976        0.0035208           266.44         0.015639        0.005705           421.75         0.019292         1.2974            459.69          3.2184        0.0053275           495.88         0.016296       0.0031679           28.938         0.004428        0.023981              37          0.014445          0.0136           44.438        0.0089112         0.01142            66.625          0.0077036        0.068403               74           0.037017         0.012242            81.438          0.0075801         7.7931              7.375           3.0192           0.15693            14.812          0.058259           2.4212             2.125           0.44071          0.55162             4.25            0.10029          9.9837    

The output is a 2-by-49 table of metrics for the frequency ranges in FB. The first row contains the metrics for healthy data while the second row contains the faulty data metrics. Observe that the following metrics have considerably higher values for the faulty data than for the healthy data:

  • PeakAmplitude2 for the frequency at 231 Hz with a difference of 3.1842 units

  • TotalBandPower with a difference of 6.01 units

Hence, use these two metrics to create a scatter plot to group the faulty and healthy data sets separately.

Create Scatter Plot

Create a scatter plot to classify healthy and faulty data using the two spectral metrics PeakAmplitude2 and TotalBandPower for all 20 data sets in the table servoData.

plotData = zeros(10,4);
for n = 1:max(size(servoData))
    hC = servoData.healthyData{n,1}.Signal;
    fC = servoData.faultyData{n,1}.Signal;
    
    [psdH,freqH] = pspectrum(hC,fs,'FrequencyResolution',0.5);
    [psdF,freqF] = pspectrum(fC,fs,'FrequencyResolution',0.5);
    
    sMetrics = faultBandMetrics({[psdH,freqH],[psdF,freqF]},FB);
    plotData(n,:) = [sMetrics{:,4}',sMetrics{:,49}'];
end
figure
scatter(plotData(:,1),plotData(:,3),[],'blue')
hold on;
scatter(plotData(:,2),plotData(:,4),[],'red')
legend('Healthy Data','Faulty Data','Location','best')
xlabel('Peak Amplitude 2')
ylabel('Total Band Power')
hold off

Observe that the healthy data sets and the faulty data sets are grouped in different areas of the scatter plot.

Hence, you can classify faulty and healthy data by analyzing the current signature of the servo motor.

To use all the available spectral metrics for classification, use Classification Learner.

Helper Function

The helper function helperPlotXLines uses the information in the structure info to plot the frequency band center lines on the power spectrum plot.

function helperPlotXLines(I,c)
for k = 1:length(I.Centers)
    xline(I.Centers(k),'Label',I.Labels(k),'LineStyle','-.','Color',c);
end
end

References

[1] Moster, P.C. "Gear Fault Detection and Classification Using Learning Machines." Sound & vibration. 38. 22-27. 2004

See Also

| | |