Matlab中的模態分析

來自專欄 勒博聲振檢測試驗室

勒博好久沒寫東西了，主要是忙著寫大論文，最近有網友問勒博，頻響函數FRF是如何得到的。這個問題很難嗎？ 勒博回憶了一下，其實對新手來說確實蠻難的，雖然課本上將了一堆模態分析，但是我們都是用較為成熟的商業軟體進行,只需要對軟體「點點點」即可。勒博其實為此困擾過很久。

今天勒博就簡單講講我們一般說的模態分析。

(1)

勒博理解的模態分析本身包含兩部分：

• 模態分解
• 模態測試

模態分解就是已知式（1）的M、C、K（中括弧[]表示方陣），而模態測試就是已知振動位移x（或振動速度，振動加速度），還有外部激勵F（大括弧{}表示它是一個列向量）

1.模態分解

現在用過簡單模型來解釋，這是一個無阻尼的三自由度模型，很簡單的吧，我們設

很簡單，我們可以得到方陣M與K

k1=1;k2=1;k3=1;k4=1;M=[2 0 0;0 1 0;0 0 3];K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3+k4];[col,row]=size(M);[D,w2]=eig(K,M);mr=D*M*D;%第r階模態質量、廣義質量kr=D*K*D;%第r階模態剛度、廣義剛度w0=w2^0.5;%無阻尼固有頻率syms H w hfor jj=1:col for kk=1:col for rr=1:col h=D(rr,jj)*D(rr,kk)/(mr(rr,rr)*(w2(rr,rr)-w^2)); if rr==1 H(jj,kk)=h; else H(jj,kk)=H(jj,kk)+h; end end endendw=0:0.001:2;for kk=1:9subplot(3,3,kk) plot(w,eval(H(kk)))end

代碼中的H 就是我們得到的頻響函數FRF。可以看出，我們是先得到模態參數（固有頻率、阻尼比、模態振型），然後根據模態參數擬合得到FRF。我們可以看見關於FRF的對稱性

同理，我們在有限元軟體中，畫好網格，輸入合適的密度、泊松比、楊氏模量、阻尼也是完成上述過程。代碼中具體理論推導大部分課本中都有詳細敘述。

我們可以看見，特徵分解是模態的核心，也就是說，eig解決了絕大部分問題。這個命令如果是阻尼系統中一樣適用，我們只需要將擴展成廣義的方程即可。

A=[C,M;M,zeros(col,row)];B=[K,zeros(col,row);zeros(col,row),-M]; [Theta,LAMD]=eig(B,-A);%λ LAMD 蘭木達 θ Theta 西塔 Ar=Theta*A*Theta;%第r階 Br=Theta*B*Theta;%第r階 r=1,2N [SS1,LAMD_index]= sort(LAMD); %排列從小到大，共軛出現 for rr=1:N SR(rr)=real(SS1(rr*2-1));% 衰減係數 SI(rr)=imag(SS1(rr*2-1));%阻尼固有頻率 w0(rr)=sqrt(SR(rr).^2+SI(rr).^2);%無阻尼固有頻率end Theta_sort=Theta(:,LAMD_index); P11=Theta_sort(1:N,1:2:2*N);%ψ Sr特徵向量 P1 P22=Theta_sort(1:N,2:2:2*N);%ψ* Sr*共軛特徵向量 P2% 模態歸一化 AR %% 求模態係數 mr=P22*M*P11;%第r階模態質量、廣義質量 kr=P22*K*P11;%第r階模態剛度、廣義剛度 cr=P22*C*P11;%第r階模態阻尼、廣義阻尼

然後勒博使個壞，關於阻尼的FRF 可以讀者可以嘗試自己編一下，提示一下無量綱阻尼係數 臨界阻尼比 。

而我們的頻響函數是什麼頻響函數呢？ 位移頻響還是加速度頻響？需要深入研究的朋友多多思考，勒博點到為止了。

2. 模態實驗

一般科研單位都用成熟的商業軟體，比如LMS TEST.LAB，如下圖所示，我們輸入振動數據 （可以加速度、速度、位移甚至雜訊,）,下標表示一個測試點,測試數據是何時間有關的列向量，對於EMA實驗，我們還要輸入激振力 （力錘或者激振器輸輸出的力信號）,力錘為例，激振下表n表示測點（下圖LMStestlab界面，如有侵權請私信告知）

而頻響函數FRF 就是根據 和 （下圖數據與照片摘自為聖小珍教授授課ppt，如有侵權請私信告知）

(2)

time=xx(:,1);%圖中時間序列acce=xx(:,2);%圖中的振動數據forces=xx(:,3);%圖中激振力數據subplot(2,1,1)plot(time,acce);gridsubplot(2,1,2)plot(time,forces);gridxlim([0 0.0005])fs=1/xx(2,1);N=length(time);omega=linspace(0,fs/2,N/2);fftacce=fftshift(fft(acce));fftforces=fftshift(fft(forces));H= fftacce(N/2+1:end)./fftforces(N/2+1:end); figure(Name, 直接FRF ); semilogy(omega,abs(H));grid xlim([4500 6500])

頻響函數就是這麼得到的。

但是實際中，我們的測試實驗由於延遲、干擾等因素，都是用功率譜估計來求

GFF=fftforces(N/2+1:end).*conj(fftforces(N/2+1:end));%F自功率譜 GXX=fftacce(N/2+1:end).*conj(fftacce(N/2+1:end));%X自功率譜 GXF=fftacce(N/2+1:end).*conj(fftforces(N/2+1:end));%xf互功率譜 GFX=fftforces(N/2+1:end).*conj(fftacce(N/2+1:end));%fx互功率譜

得到功率譜後，我們開始進行H1、H2、Hv估計，就是圖3中右下角的那個選項

H1=GXF./GFF;H2=GXX./GFX;

勒博繼續挖坑，Hv估計這裡留個讀者自己研究。

與模態分解不同，模態測試時先得到頻響函數，然後再根據頻響函數用一些方法去識別模態參數（如實部的曲線擬合）.

好了，勒博又要去寫自己的論文了， 希望這篇文章對您有所啟發。

（寫的有點亂，如要轉發，請私信）