作者：趙兵
日期：2020-02-17
目錄
單自由度系統中質量、阻尼和剛度變化對頻率響應函數(FRF)影響圖的繪制
1. 背景
2. VISIO繪制
3. Matlab繪制
(1) M變化時
(2) K變化時
(3) C變化時
4. 參考文章
1. 背景
寫文章時需要用到幾張圖,下面是從PDF上截圖截出來的，用來表示單自由度系統在沖擊激勵下的頻率響應曲線，當K（剛度）,C（阻尼）,M（質量）變化時，頻率響應曲線的變化情況。
 |
 |
 |
| 圖 1 單自由度系統剛度，阻尼，質量影響曲線 |
||
用圖1放在文章中，不太美觀，一看就是影印的，要是一般的文章還好，如果用來發表的就拉下檔次了。所以就嘗試了下重繪，重繪有兩個方案，
（1）用PPT、visio或者Adobe Illustrator等繪圖工具自己繪制
（2）用matlab先得到曲線的函數，進一步把函數顯示出來
2. Visio繪制
嘗試了一下，質量變化影響曲線圖相對還好畫一些，但阻尼變化的曲線就不是那么容易畫了，畫出來總覺得差點意思。很難保證各曲線之間的間隔均勻；當然肯定是能做到，但自己非此方面的熟手，所以有兩個解決方案，一個是嘗試用Adobe Illustrator進行一定時間的專門學習，但時間成本太高；或者是淘寶上尋求供外包解決，也是能解決的；但最終我還是決定先嘗試下第二種方案，先得到曲線的函數，然后用軟件把函數繪出來。
 |
| 圖 2 Visio繪制曲線圖 |
3. Matlab繪制
單自由度系統的物理模型如圖 3 所示，
 |
| 圖 3單自由度系統模型 |
它的動力學方程為
其中
$M$:質量; $C$:阻尼; $K$:剛度; $\ddot{x}$:加速度;$\dot{x}$:速度; $x$:位移; $f$:外力; $t$:時間。
這里的f(t)為脈沖激勵：
使用漢寧窗生成：$ f(t<T_c)= \frac{A}{2} * cos(\frac{2*\pi* t(t<T_c)}{T_c}) $
% % 激勵采用脈沖激勵，脈沖激勵為Hanning函數
function f = hanning_imp(t, Tc, A)
f = zeros(size(t));
f(t < Tc) = A / 2 * (1 - cos(2*pi * t(t < Tc) / Tc));
end
可以根據這個函數得到一個脈沖激勵
畫圖此圖：
dt = 0.00001;
t = 0:dt:200;
Tc = 0.001;
A = 10;
u = hanning_imp(t, Tc, A);
plot(t,u,'LineWidth',1.5 );
axis([-5 200 -0.1 10.5 ]);
xlabel('t/s')
ylabel('Amp/N')
text(75,8,'Impact Force');
|
（a）整體圖 |
(b)局部放大圖 |
|
| 圖 4 沖擊力
|
|
|
建立系統方程，求解頻率響應函數（FRF）
function [freq_x , amp_y]=frf_bing(m,k,c) % m 質量 % k 剛度 % c 阻尼 num = 1; den = [m c k]; sys = tf(num, den); %采樣頻率(Hz) 100Hz 實際並不需要這么高的采樣頻率，但是如果采樣時間太小，hanning脈沖不完整 % 為了得到准確的響應dt一定要小，否則做出的相位可能不對 dt = 0.00001; fs = 1/dt; t = 0:dt:200; Tc = 0.001; A = 10; u = hanning_imp(t, Tc, A); y = lsim(sys, u, t); y = y'; N = length(u); fy = fft(y); fu = fft(u); ft = fy ./ fu; f = (0:N-1) * fs ./ N; ft_r = real(ft); ft_i = imag(ft); part = (f < 30); freq_x=f(part); amp_y= abs(ft(part)); End
(1) M變化時
取m=100,120,140,160,180,200分別繪制FRF的響應曲線。
clc;
clear;
close all;
%%
k = 1000; %初始化k
c=100; %初始化c
M=100:20:200; %初始化m, 取m=100,120,140,160,180,200分別繪圖
f1= figure(1);
hold on
for i= 1:length(M)
[a(:,i) , b(:,i)]=frf_bing(M(i),k,c);
end
plot(a,log(b),'b');
axis([0 1 -8.2 -5.3]);
title('Log(FRF)');
xlabel('Frequency')
ylabel('Amplitude')
set(gca,'XTick',[],'YTick',[]);
text(0.3,-7,'M↑');
annotation('arrow',[0.7 0.3],[0.6 0.8]) ;
f1.Position=([ 0 0 400 300])
% 保存為emf 矢量格式
set(gcf,'unit','centimeters','position',[10 5 6.5 4.8]);
print(f1,'-dmeta','M.emf')
可以生成一幅想要的曲線，並且可以保存為矢量格式，無論放大多少倍，圖片還是很清晰。
 |
| 圖 5 質量變化的影響 |
(2) K變化時
取K=1000,1200,1400,1600,1800,2000分別繪制FRF的響應曲線。
clc;
clear;
close all;
%%
k = 1000:200:2000;
c=100;
M=100;
f1= figure(1);
hold on
for i= 1:length(k)
[a(:,i) , b(:,i)]=frf_bing(M,k(i),c);
end
plot(a,log(b),'b');
axis([0 1 -8.2 -5.3]);
title('Log(FRF)');
xlabel('Frequency')
ylabel('Amplitude')
set(gca,'XTick',[],'YTick',[]);
text(0.6,-7,'K↑');
annotation('arrow',[0.4 0.8],[0.63 0.6]) ;
f1.Position=([ 0 0 400 300])
% 保存為emf 矢量格式
set(gcf,'unit','centimeters','position',[10 5 6.5 4.8]);
print(f1,'-dmeta','K.emf')
|
|
| 圖 6剛度變化的影響 |
(3) C變化時
取K=1000,1200,1400,1600,1800,2000分別繪制FRF的響應曲線。
clc; clear; close all; %%
fontsizevalue=18; c=100:20:200; f1= figure(1); hold on for i= 1:length(c) [a(:,i) , b(:,i)]=frf_bing(c(i)); end plot(a,log(b),'b'); axis([0 1 -8 -5.5]); t1=title('Log(FRF)'); xl=xlabel('Frequency') y1=ylabel('Amplitude') t1.FontSize =fontsizevalue; y1.FontSize =fontsizevalue; xl.FontSize =fontsizevalue; set(gca,'XTick',[],'YTick',[]); text(0.5,-7,'C↑'); annotation('arrow',[0.5 0.5],[0.9 0.5]) ; f1.Position=([ 0 0 400 300]) % 保存為emf 矢量格式 set(gcf,'unit','centimeters','position',[10 5 6.5 4.8]); print(f1,'-dmeta','K.emf')
 |
| 圖 7阻尼變化的影響 |
4. 參考文章
1. CSDN博主「whoispo」文 https://blog.csdn.net/WhoisPo/article/details/46865401
