论文复现（螺杆压缩机典型故障振动分析及诊断应用研究）-CSDN博客

本文链接：https://blog.csdn.net/Lvyou2008/article/details/132811312

声明：

本次论文复现仅为笔者本人的个人理解，正确与否欢迎斧正，并不对论文质量作出评价。
本文提供的代码仅作参考

本次论文复现的是《螺杆压缩机典型故障振动分析及诊断应用研究》，关注的是文章中第三章螺杆压缩机振动信号分析及故障特征提取方式，第27页中描述的二维全息图公式

由于三角函数和差公式转换较为复杂，将xi，和yi用复数方式表达，方便后续推导应该

$x_{i}=A_{1i}cos\left ( \omega _{i} t+\psi _{1i}\right )$

$\frac{x_{i}}{A_{1i}}=cos\left ( \omega_{i} t+\psi _{1i}\right )=\frac{e^{i\left ( \omega _{i}t+\psi _{1i}\right )}-e^{-i\left ( \omega _{i}t+\psi _{1i}\right )}}{2}$ \

$\frac{y_{i}}{A_{2i}}=cos\left ( \omega_{i} t+\psi _{2i}\right )=\frac{e^{i\left ( \omega _{i}t+\psi _{2i}\right )}-e^{-i\left ( \omega _{i}t+\psi _{2i}\right )}}{2}$

$\frac{x^{2}}{A_{1i}^{2}}+\frac{y^{2}}{A_{2i}^{2}}-\frac{2x_{i}y_{i}}{A_{1i}A_{2i}}cos\left ( \psi _{2i} - \psi_{1i} \right )=\frac{1}{2}cos\left ( 2\psi _{2i} -2\psi_{1i} \right )$

此处的推导结果与文章的结果相异，根据笔者推导的结果，与椭圆曲线相对应，得出系数如下

椭圆曲线：

$ax^{2}+bxy+cy^{2}+dx+ey+f=0$

得出：

$a=\frac{1}{A_{1i}^{2}}$

$b=\frac{-2cos\left ( \psi _{2i} -\psi _{1i}\right )}{A_{1i}A_{2i}}$

$c = \frac{1}{A_{2i}}$

$d=0$

$e=0$

$f = -\frac{1}{2}cos\left ( \psi _{2i} -\psi _{1i}\right )$

根据椭圆曲线长轴公式可知如下

$\Delta =b^{2}-4ac$

$\delta _{1} = f\Delta +\left ( ae^{2}+cd^{2} -bde\right )$

$\delta _{2}=d^{2}+e^{2}-4\left ( a+c\right )f$

$\sigma _{1}=ae^{2}+cd^{2}-bde$

$\sigma_{2}=\sqrt{\left ( a+c \right )+\Delta }$

当δ1＞0时，实半长轴为

$\sqrt{\left | \frac{2\delta _{1}}{\Delta \times \left [ \sigma _{2} -\left ( a+c \right )\right ]} \right |}$

虚半长半轴为

$\sqrt{\left | \frac{-2\delta _{1}}{\Delta \times \left [ \sigma _{2} +\left ( a+c \right )\right ]} \right |}$

当δ1＜0时，实半长轴为

$\sqrt{\left | \frac{-2\delta _{1}}{\Delta \times \left [ \sigma _{2} +\left ( a+c \right )\right ]} \right |}$

虚半长半轴为

$\sqrt{\left | \frac{2\delta _{1}}{\Delta \times \left [ \sigma _{2} -\left ( a+c \right )\right ]} \right |}$

以网上凯斯西储大学中轴承数据为算例，计算振动模态的全息图谱代码如下。

clear all;

clf;

tic

D = load('E:\电机振动模态和故障预警\凯斯西储大学数据\12k Fan End Bearing Fault Data\273.mat');

Data1 = getfield(D,'X273_DE_time');

Data2 = getfield(D,'X273_FE_time');

Data3 = getfield(D,'X273_BA_time');

Data = cat(2,Data1,Data2,Data3);

t_s = 0.008;

t_start = 0;

f_s = 1/t_s;

for i = 1:2

    t_end = length(Data(:,i))*0.002;

    Druation = t_end - t_start;

    Sampling_piont = Druation/t_s + 1;

    f_x = 0:(f_s/(Sampling_piont-1)):f_s;

    y(:,i) = Data(:,i);

    y_f(:,i) = fft(y(:,i));

    Fuzhi(:,i) = abs(y_f(:,i));

    Psi(:,i) = angle(y_f(:,i));

end

Fuzhi1 = Fuzhi((1:length(f_x)),1);

Fuzhi2 = Fuzhi((1:length(f_x)),2);

Psi1 = Psi((1:length(f_x)),1);

Psi2 = Psi((1:length(f_x)),2);

f = diag(f_x);

t = (t_start:0.002:4)';

T = repmat(t,1,length(Psi1));

G = T*f;

psi1 = repmat(Psi1',length(t),1);

psi2 = repmat(Psi2',length(t),1);

B1 = cos(G + psi1);

B2 = cos(G + psi2);

Fuzhi_x = diag(Fuzhi1);

Fuzhi_y = diag(Fuzhi2);

x = B1*Fuzhi_x;

y = B2*Fuzhi_y;

for i =1:20

    figure(i);

    plot(x(:,i+6660*0.15),y(:,i+6660*0.15));

    hold on;

    axis equal;

end

toc