簡介這份資源面向機械工程、設(shè)備維護與信號處理方向的學(xué)習(xí)者和研究人員圍繞凱斯西儲大學(xué)軸承數(shù)據(jù)提供軸承故障特征頻率計算的MATLAB實現(xiàn)幫助讀者從振動信號中識別滾動體、內(nèi)圈、外圈等故障對應(yīng)的特征頻率解決故障診斷中理論公式與代碼落地脫節(jié)的問題。壓縮包共3個文件包含1個m腳本和2個png圖片整體約91KB腳本用于特征頻率計算與頻譜分析圖片則給出深溝球軸承結(jié)構(gòu)及計算截圖便于對照理解參數(shù)含義。目前已有4135人學(xué)習(xí)下載說明其在軸承故障診斷入門與課程實踐中具有一定參考價值。讀者可借助代碼復(fù)現(xiàn)基本旋轉(zhuǎn)頻率、滾動體通過頻率等計算流程結(jié)合步驟分析講解掌握數(shù)據(jù)預(yù)處理、頻譜繪制與故障判斷思路適合作為機械故障診斷、信號處理課程設(shè)計或課題研究的實操素材。1. 軸承故障特征頻率為什么你的包絡(luò)譜總是對不上號做過旋轉(zhuǎn)機械診斷的人多半有過這種經(jīng)歷拿一段凱斯西儲大學(xué)CWRU的驅(qū)動端加速度數(shù)據(jù)FFT 一看轉(zhuǎn)頻及其倍頻清清楚楚可一旦想找內(nèi)圈、外圈的故障特征頻率譜線要么淹沒在噪聲里要么算出來的位置和理論值差了一截。問題往往不在數(shù)據(jù)而在“特征頻率到底怎么算、算出來之后怎么在譜上認”。這個標題指向的正是這件事用 MATLAB 把軸承的幾何參數(shù)轉(zhuǎn)成理論故障頻率再把它和實測包絡(luò)譜對齊。它解決的是“知道有故障但說不清是哪個部件、哪一階”的問題適合手里有振動數(shù)據(jù)、想從頻譜里讀出物理含義的設(shè)備診斷工程師和做狀態(tài)監(jiān)測的學(xué)生。核心不是寫多復(fù)雜的代碼而是把幾何—頻率—譜線這條鏈路走通讓每一次峰值都能落到具體部件上。2. 從幾何參數(shù)到特征頻率四個公式和它們的物理含義2.1 軸承特征頻率的推導(dǎo)邏輯滾動軸承的故障特征頻率來自滾動體與滾道的相對運動。設(shè)節(jié)徑為 $D$滾動體直徑為 $d$接觸角為 $\alpha$滾動體個數(shù)為 $Z$軸轉(zhuǎn)頻為 $f_r$。外圈固定、內(nèi)圈隨軸旋轉(zhuǎn)時四個最常用的頻率是外圈故障頻率 BPFO$\frac{Z}{2} f_r \left(1 - \fracdrlblxrh1{D}\cos\alpha\right)$內(nèi)圈故障頻率 BPFI$\frac{Z}{2} f_r \left(1 \fracdrlblxrh1{D}\cos\alpha\right)$滾動體故障頻率 BSF$\frac{D}{2d} f_r \left[1 - \left(\fracdrlblxrh1{D}\cos\alpha\right)^2\right]$保持架頻率 FTF$\frac{1}{2} f_r \left(1 - \fracdrlblxrh1{D}\cos\alpha\right)$這幾個式子不是背下來的而是從“滾動體經(jīng)過一個缺陷點的時間間隔”推出來的。外圈有缺陷時每個滾動體經(jīng)過缺陷都會產(chǎn)生一次沖擊所以頻率是滾動體公轉(zhuǎn)頻率乘以個數(shù)內(nèi)圈有缺陷時缺陷隨內(nèi)圈轉(zhuǎn)動沖擊頻率還要疊加內(nèi)圈自身的旋轉(zhuǎn)。理解這一點后面看譜線時才知道為什么 BPFI 附近常出現(xiàn)轉(zhuǎn)頻邊帶而 BPFO 的邊帶通常不明顯。2.2 用 MATLAB 把公式寫成可復(fù)用函數(shù)直接在主腳本里寫死數(shù)字換一個軸承就要改一遍容易出錯。我一般會寫一個獨立函數(shù)把幾何參數(shù)和轉(zhuǎn)頻作為輸入返回四個頻率。下面這段代碼可以直接存成bearing_freq.mfunction [BPFO, BPFI, BSF, FTF] bearing_freq(D, d, Z, alpha_deg, fr) % D 節(jié)徑 (mm) % d 滾動體直徑 (mm) % Z 滾動體個數(shù) % alpha_deg 接觸角 (度) % fr 軸轉(zhuǎn)頻 (Hz) alpha deg2rad(alpha_deg); % 角度轉(zhuǎn)弧度 ratio (d / D) * cos(alpha); % 中間量反復(fù)用到 BPFO (Z / 2) * fr * (1 - ratio); BPFI (Z / 2) * fr * (1 ratio); BSF (D / (2 * d)) * fr * (1 - ratio^2); FTF 0.5 * fr * (1 - ratio); end邏輯說明ratio是四個公式共有的無量綱量單獨提出來既減少重復(fù)計算也方便檢查是否超過 1。接觸角用角度輸入更符合工程習(xí)慣函數(shù)內(nèi)部轉(zhuǎn)弧度。參數(shù)說明節(jié)徑和滾動體直徑單位要一致用毫米即可因為比值無量綱轉(zhuǎn)頻由轉(zhuǎn)速除以 60 得到例如 1750 rpm 對應(yīng)約 29.17 Hz。如果軸承資料給的是接觸角弧度把deg2rad去掉即可。2.3 一個可復(fù)現(xiàn)的算例與結(jié)果核對以某型深溝球軸承為例節(jié)徑 39.04 mm滾動體直徑 7.94 mm滾動體個數(shù) 9接觸角 0 度軸轉(zhuǎn)頻 29.17 Hz。調(diào)用上面的函數(shù)D 39.04; d 7.94; Z 9; alpha 0; fr 1750 / 60; [BPFO, BPFI, BSF, FTF] bearing_freq(D, d, Z, alpha, fr); fprintf(BPFO%.2f Hz\nBPFI%.2f Hz\nBSF%.2f Hz\nFTF%.2f Hz\n, ... BPFO, BPFI, BSF, FTF);運行后大致得到 BPFO 約 104.6 HzBPFI 約 157.9 HzBSF 約 68.9 HzFTF 約 11.6 Hz。這幾個數(shù)字要記下來后面做包絡(luò)譜時就是拿它們?nèi)Ψ?。注?BSF 的公式里有一個平方項當(dāng) $d/D$ 較大時結(jié)果對參數(shù)很敏感所以滾動體直徑的測量誤差會明顯影響 BSF 的定位這也是實際診斷中 BSF 最難對上的原因之一。3. 用 MATLAB 做包絡(luò)譜把故障頻率從噪聲里撈出來3.1 為什么直接 FFT 不夠要先做包絡(luò)軸承故障產(chǎn)生的沖擊是高頻衰減振蕩直接對原始信號做 FFT能量集中在幾千赫茲的共振頻帶故障特征頻率本身反而很弱。包絡(luò)分析的思路是先帶通濾波保留共振頻帶再取包絡(luò)最后對包絡(luò)做低通或直接 FFT這樣故障頻率就搬到低頻段容易識別。常見做法是 Hilbert 變換取包絡(luò)也可以用平方解調(diào)。我一般先用譜峭度或經(jīng)驗選一個帶通范圍再進包絡(luò)。3.2 完整包絡(luò)譜腳本與參數(shù)設(shè)置下面這段腳本假設(shè)你已經(jīng)把 CWRU 的驅(qū)動端數(shù)據(jù)讀成變量x采樣率fs已知常見為 12000 Hz 或 48000 Hz。代碼里帶通范圍需要根據(jù)你的數(shù)據(jù)調(diào)整這里給一個示例值fs 12000; % 采樣率按實際數(shù)據(jù)改 x x - mean(x); % 去直流 bp [2000 4000]; % 帶通范圍圍繞共振頻帶選 % 帶通濾波 [b, a] butter(4, bp / (fs / 2), bandpass); xf filtfilt(b, a, x); % 零相位濾波避免相位失真 % Hilbert 包絡(luò) env abs(hilbert(xf)); env env - mean(env); % 去直流否則零頻會很高 % 包絡(luò)譜 N length(env); Y abs(fft(env)) / N; f (0:N-1) * (fs / N); half 1:floor(N/2); plot(f(half), Y(half)); xlabel(頻率 (Hz)); ylabel(幅值); xlim([0 500]); % 故障頻率一般在低頻段 grid on;邏輯說明filtfilt做雙向濾波相位不失真對包絡(luò)提取很關(guān)鍵hilbert返回解析信號取模得到包絡(luò)去均值是為了壓掉零頻分量否則譜圖第一個點會蓋住一切。參數(shù)說明帶通范圍bp不能拍腦袋可以先對原始信號做 FFT 看共振峰在哪再把范圍卡在峰兩側(cè)濾波器階數(shù) 4 是折中太高會振鈴太低阻帶衰減不夠。xlim先看 0 到 500 Hz因為常見軸承故障頻率都在這個區(qū)間確認有峰后再放大。3.3 把理論頻率疊到包絡(luò)譜上做比對算出 BPFO、BPFI 之后最直觀的驗證方式是在包絡(luò)譜上畫豎線。下面幾行接在上面的腳本后面hold on; for k 1:5 xline(k * BPFO, r--); % 外圈及其諧波 xline(k * BPFI, b--); % 內(nèi)圈及其諧波 end legend(包絡(luò)譜, BPFO 諧波, BPFI 諧波); hold off;如果外圈故障你會看到 BPFO 及其 2 倍、3 倍處有峰內(nèi)圈故障則在 BPFI 兩側(cè)出現(xiàn)間隔為轉(zhuǎn)頻的邊帶。這一步是“認譜”的核心理論值只是參考實際峰值可能因為滑移偏幾個百分點允許 1% 到 2% 的偏差。若偏差超過 5%先檢查轉(zhuǎn)頻是否算錯再檢查軸承型號參數(shù)是否對得上。4. 避坑與排查五個讓特征頻率對不上的常見原因4.1 轉(zhuǎn)頻算錯后面全錯現(xiàn)象所有特征頻率整體偏移比例大致相同。原因把轉(zhuǎn)速 rpm 直接當(dāng) Hz 用或者用了電機同步轉(zhuǎn)速而不是實際轉(zhuǎn)速。解決轉(zhuǎn)頻 實測轉(zhuǎn)速 / 60沒有轉(zhuǎn)速計時可從譜圖上的轉(zhuǎn)頻峰反推確認后再代入公式。4.2 帶通范圍選在噪聲區(qū)現(xiàn)象包絡(luò)譜一片模糊看不到明顯峰值。原因帶通范圍沒有覆蓋共振頻帶或者范圍太寬把噪聲也放進來。解決先看原始信號 FFT找能量最高的共振峰把帶通卡在峰兩側(cè)約 ±500 Hz如果共振峰不明顯用譜峭度圖輔助選帶。4.3 濾波器相位失真導(dǎo)致包絡(luò)變形現(xiàn)象包絡(luò)譜出現(xiàn)很多不明諧波和理論頻率對不上。原因用了filter而不是filtfilt相位延遲讓包絡(luò)形狀改變。解決包絡(luò)提取前一律用零相位濾波filtfilt是標配。4.4 軸承幾何參數(shù)抄錯或單位不統(tǒng)一現(xiàn)象BPFI 和 BPFO 的比值異常比如內(nèi)圈頻率反而比外圈低。原因節(jié)徑和滾動體直徑一個用毫米一個用英寸或者滾動體個數(shù)抄錯。解決所有長度單位統(tǒng)一接觸角確認是度還是弧度算完先檢查 BPFI 是否大于 BPFO正常情況內(nèi)圈頻率更高。4.5 滑移導(dǎo)致峰值偏移誤判為無故障現(xiàn)象理論頻率處沒有峰但附近幾十赫茲有峰。原因滾動體存在滑移實際沖擊間隔比理論值略長頻率偏低。解決在理論值 ±2% 范圍內(nèi)搜索峰值不要死盯一個點如果偏移穩(wěn)定可反推實際滑移率但不必強行修正公式。5. 進階技巧用包絡(luò)譜的諧波結(jié)構(gòu)區(qū)分內(nèi)外圈故障5.1 看邊帶而不是只看主峰外圈故障時沖擊位置相對傳感器固定包絡(luò)譜上 BPFO 及其諧波干凈利落邊帶很少。內(nèi)圈故障時缺陷隨內(nèi)圈轉(zhuǎn)動與傳感器距離周期性變化BPFI 兩側(cè)會出現(xiàn)間隔為轉(zhuǎn)頻的邊帶。這個差別比主峰幅值更可靠。實操中我會把包絡(luò)譜在 BPFI 附近放大看有沒有BPFI ± fr的峰有就是內(nèi)圈沒有再看 BPFO。5.2 用倒譜做二次確認倒譜能把諧波族壓縮成一根線適合確認故障頻率的基頻。下面這段代碼接在包絡(luò)譜之后cep real(ifft(log(abs(fft(env)) eps))); q (0:N-1) / fs; plot(q(1:floor(N/2)), cep(1:floor(N/2))); xlabel(倒頻率 (s)); ylabel(幅值); xlim([0 0.1]); grid on;邏輯說明log前加eps防止零值取對數(shù)出錯倒譜橫軸是時間故障頻率的倒數(shù)處會出現(xiàn)峰。參數(shù)說明xlim取 0 到 0.1 秒對應(yīng) 10 Hz 以上頻率覆蓋常見故障頻率。如果倒譜在 1/BPFO 處有峰基本可以確認外圈。5.3 一個我常犯的錯誤早期我總想把所有峰值都解釋成故障頻率結(jié)果把轉(zhuǎn)頻諧波、齒輪嚙合頻率也硬套進去越看越亂。后來養(yǎng)成習(xí)慣先確認轉(zhuǎn)頻和它的諧波把它們從候選列表里劃掉再在剩下的峰里找 BPFO、BPFI。這個順序能省掉大量返工。另外CWRU 數(shù)據(jù)里不同負載和轉(zhuǎn)速對應(yīng)不同文件用之前一定核對轉(zhuǎn)速別拿 1730 rpm 的參數(shù)去算 1750 rpm 的數(shù)據(jù)。希望幫到你。本文還有配套的精品資源點擊獲取