杭州锐达数字技术有限公司
查看: 872|回复: 4

[综合讨论] 振动方程,数值计算,特征频率

[复制链接]
发表于 2017-6-3 22:17 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
对三自由度振动微分方程组采用线性加速度法进行数值计算,可以得到三个自由度的位移、速度及加速度的时程曲线,对三个自由度中某一自由度的加速度进行FFT变换时,从加速度的频响曲线中应该可以看到三个特征频率,却知道看两个,少了一个特征频率,请各位大侠指教。振动微分方程组的质量矩阵m=[1 0 0;0 1 0;0 0 1]刚度矩阵为k=50*[2 -1 0;-1 2 -2;0 -1 2],理论计算结果得到的三个特征频率为0.5825Hz,1.5915Hz,2.1741Hz。数值计算程序如下。
close all;
clc;
clear all;
m=[1 0 0;0 1 0;0 0 1];                                           %质量矩阵
% c=[2 -1 0;-1 2 -1;0 -1 2];                                       %阻尼矩阵
c=[0 0 0;0 0 0;0 0 0];      
% k=50*[2 -1 0;-1 2 -2;0 -1 2];
k=50*[2 -1 0;-1 2 -2;0 -1 2];                                    %刚度矩阵
x0=[1 1 1]';                                                     %初位移
v0=[1 1 1]';                                                     %初速度
delt=1/(12);                                                     %时间步长
fs=1/delt;
time=25;                                                         %仿真时间
n=time/delt;                                                     %循环次数
disp=zeros(n,3);                                   %设定n行3列存储加速度矩阵
minv=inv(m+delt*c/2+delt^2*k/6);
i=1;
for t=0:delt:time
%     f=[6*sin(3.5*t) -2*cos(2*t) 1.5*sin(1.5*t)]';
    f=[0 0 0]';
    if t==0
        a0=inv(m)*(f-k*x0-c*v0);                            %计算初始加速度
    else
        a=minv*(f-c*(v0+delt/2*a0)-k*(x0+delt*v0+delt^2*a0/3));%计算加速度
        v=v0+delt*(a0+a)/2;                                    %计算速度
        x=x0+delt*v0+delt^2/3*a0+delt^2/6*a;                     %计算位移
        a0=a;v0=v;x0=x;i=i+1;
    end
    disp(i,:)=a0;                                   %3个自由度的加速度矩阵

end
t=0:delt:time;
figure
% plot(t,disp(:,1),t,disp(:,2),t,disp(:,3));
plot(t,disp(:,1));    %三个自由度方向上加速度仿真曲线
grid on ;
xlabel('时间(s)'),title('3自由度时程曲线');



A=disp(:,1);                                   %第一自由度方向的加速度向量
M=length(A);
n=0:M-1;
figure
plot((0:M-1)*fs/M,20*log10(abs(fft(disp(:,1)))));
grid;



仿真的图为
QQ图片20170603221536.png
回复
分享到:

使用道具 举报

 楼主| 发表于 2017-6-4 22:07 | 显示全部楼层
有哪位大侠帮我看一看问题出来哪里?
发表于 2017-6-14 10:13 | 显示全部楼层
我看到有人帮你改出来了  是把初始激励改为零吧  
发表于 2017-10-8 16:30 | 显示全部楼层
楼主解决了吗
发表于 2019-4-19 21:20 | 显示全部楼层
学习下
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表