单自由度振动动画

jqche

太牛了

牛小贱

虽然该帖子有点“久远”了，我还是顶一下，顺便上传一个本人在学习《汽车振动学》时，做的一个Matlab程序，和楼主的差不多，但是相对来说比较简单，对于初学matlab（编程）的人来说，比较容易理解，对学习matlab编程有用！以下是程序：
clear all;
clear;
clc;
rectangle('position',[12,1.2,2,0.2],'FaceColor',[0.5,0.3,0.4]);
axis([0,15,-5,15]);
%画顶板
hold on
%画直线
y=1.5:.35:5.5;
M=length(y);
x=12+mod(1:M,2)*2;
x(1)=13;
x(end-3:end)=13;
D=plot(x,y);
%弹簧 D
C=(1/8:1/4:7/8)*2*pi;r=0.35;
t1=2*r*sin(C);
F1=fill(13+3*r*cos(C),5.5+t1,'r');
% 球
% set(gca,'ytick',[0:2:9]);
% set(gca,'yticklabels',num2str([-1:3]'));
plot([0,15],[5.9,5.9],'black');
H1=plot([0,13],[5.9,5.9],'y');
% 句柄[黄线]
Q=plot(0,5.5,'color','r');
% 运动曲线;
td=[];yd=[];
text(2,14,'阻尼振动','fontsize',16);
set(gcf,'doublebuffer','on');
%有阻尼自由衰减振动分析计算程序
k=80;   %N/m
m=10;   %Kg
c=0.219*sqrt(k*m)*2;%N.s/m
x0=-2;
v0=-10;
X1=[x0,v0]';
p=sqrt(k/m);
n=c/(2*m);
f=p/(2*pi);
ksai=n/p;
T0=1/f;
T1=T0*(1/sqrt(1-ksai^2));
T1_T=T1/T0;
A2_A1=1/exp(n*T1);
s1=-n+sqrt(n^2-p^2);
s2=-n-sqrt(n^2-p^2);
s=[1,1;s1,s2];
Cd=inv(s)*X1;
T=0;
while T<12;
pause(0.2);
Dy=real(Cd(1)*exp(s1*T)+Cd(2)*exp(s2*T))/4+1.15;
Y=(y-1.5)*Dy+1.5;
Yf=Y(end)+t1;
td=[td,T];yd=[yd,Y(end)];
set(D,'ydata',Y);
set(F1,'ydata',Yf,'facecolor',rand(1,3));
set(H1,'xdata',[T,13],'ydata',[Y(end),Y(end)]);
set(Q,'xdata',td,'ydata',yd) ;
T=T+0.1;
end

沧海一声笑

牛小贱 发表于 2013-12-11 22:52

                               
登录/注册后可看大图

虽然该帖子有点“久远”了，我还是顶一下，顺便上传一个本人在学习《汽车振动学》时，做的一个Matl ...

谢谢分享！！

牛小贱

为方便大家复制代码，我重新整理一下程序！

1. clear all;
   
2. clear;
   
3. clc;
   
4. rectangle('position',[12,1.2,2,0.2],'FaceColor',[0.5,0.3,0.4]);
   
5. axis([0,15,-5,15]);
   
6. %画顶板
   
7. hold on
   
8. %画直线
   
9. y=1.5:.35:5.5;
   
10. M=length(y);
    
11. x=12+mod(1:M,2)*2;
    
12. x(1)=13;
    
13. x(end-3:end)=13;
    
14. D=plot(x,y);
    
15. %弹簧 D
    
16. C=(1/8:1/4:7/8)*2*pi;r=0.35;
    
17. t1=2*r*sin(C);
    
18. F1=fill(13+3*r*cos(C),5.5+t1,'r');
    
19. % 球
    
20. % set(gca,'ytick',[0:2:9]);
    
21. % set(gca,'yticklabels',num2str([-1:3]'));
    
22. plot([0,15],[5.9,5.9],'black');
    
23. H1=plot([0,13],[5.9,5.9],'y');
    
24. % 句柄[黄线]
    
25. Q=plot(0,5.5,'color','r');
    
26. % 运动曲线;
    
27. td=[];yd=[];
    
28. text(2,14,'阻尼振动','fontsize',16);
    
29. set(gcf,'doublebuffer','on');
    
30. %有阻尼自由衰减振动分析计算程序
    
31. k=80;   %N/m
    
32. m=10;   %Kg
    
33. c=0.219*sqrt(k*m)*2;%N.s/m
    
34. x0=-2;
    
35. v0=-10;
    
36. X1=[x0,v0]';
    
37. p=sqrt(k/m);
    
38. n=c/(2*m);
    
39. f=p/(2*pi);
    
40. ksai=n/p;
    
41. T0=1/f;
    
42. T1=T0*(1/sqrt(1-ksai^2));
    
43. T1_T=T1/T0;
    
44. A2_A1=1/exp(n*T1);
    
45. s1=-n+sqrt(n^2-p^2);
    
46. s2=-n-sqrt(n^2-p^2);
    
47. s=[1,1;s1,s2];
    
48. Cd=inv(s)*X1;
    
49. T=0;
    
50. while T<12;
    
51. pause(0.2);
    
52. Dy=real(Cd(1)*exp(s1*T)+Cd(2)*exp(s2*T))/4+1.15;
    
53. Y=(y-1.5)*Dy+1.5;
    
54. Yf=Y(end)+t1;
    
55. td=[td,T];yd=[yd,Y(end)];
    
56. set(D,'ydata',Y);
    
57. set(F1,'ydata',Yf,'facecolor',rand(1,3));
    
58. set(H1,'xdata',[T,13],'ydata',[Y(end),Y(end)]);
    
59. set(Q,'xdata',td,'ydata',yd) ;
    
60. T=T+0.1;
    
61. end
    

复制代码

补充一个多自由度的振动程序（参考张亚辉结构动力学基础）

1. function newmark(t_max,dt)
   
2. %t_max为持续时间，dt为时间步长
   
3. M=1*[1 0 0;0 1 0;0 0 1];
   
4. K=1*[3 -1 0;-1 2 -1;0 -1 1];
   
5. C=0.015*M+0.02*K;%c=0.015*[0 0 0;0 0 0;0 0 0];
   
6. u=[0 0 0]';
   
7. v=[0 0 0]';
   
8. a=[0 0 0]';
   
9. t(1)=0;               %时间
   
10. x(:,1)=u;             %位移
    
11. x1(:,1)=v;            %速度
    
12. x2(:,1)=a;            %加速度
    
13. %newmark参数
    
14. gama=0.5;
    
15. delta=0.25;
    
16. a0=1/(delta*dt^2);
    
17. a1=gama/(delta*dt);
    
18. a2=1/(delta*dt);
    
19. a3=1/(2*delta)-1;
    
20. a4=gama/delta-1;
    
21. a5=dt*(gama/(2*delta)-1);
    
22. a6=dt*(1-gama);
    
23. a7=gama*dt;
    
24. %等效刚度矩阵
    
25. K1=K+a0*M+a1*C;
    
26. i=1;
    
27. t(1)=0;
    
28. while t(i)<t_max
    
29.       Q=-1*[sin(t(i)) sin(t(i)) sin(t(i))]';
    
30.       q(:,i+1)=Q+M*(a0*x(:,i)+a2*x1(:,i)+a3*x2(:,i))+C*(a1*x(:,i)+a4*x1(:,i)+a5*x2(:,i));
    
31.       x(:,i+1)=inv(K1)*q(:,i+1);
    
32.       x2(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*x1(:,i)-a3*x2(:,i);
    
33.       x1(:,i+1)=a1*(x(:,i+1)-x(:,i))-a4*x1(:,i)-a5*x2(:,i);
    
34.       i=i+1;
    
35.       t(i)=t(i-1)+dt;
    
36. end  
    
37. x1=x(1,:)';
    
38. x2=x(2,:)';
    
39. x3=x(3,:)';
    
40. xx1=x1(1:t_max/dt)';
    
41. xx2=x2(1:t_max/dt)';
    
42. xx3=x3(1:t_max/dt)';
    
43. t=dt:dt:t_max;
    
44. plot(t,xx1,'-',t,xx2,'-.',t,xx3,'--')
    
45. xlabel('t[sec]')
    
46. ylabel('x(m)')
    
47. legend('x1(t)','x2(t)','x3(t)')
    
48. grid
    

复制代码

在Matlab命令窗口输入：输入：newmark(30,0.01)
另外，大家还可以参考一下这个帖子：http://forum.vibunion.com/thread-30794-1-1.html（关于多自由度 自由振动 MATLAB 仿真 的，里面有程序，又仿真结果，很详细！！）

牛小贱

沧海一声笑 发表于 2013-12-12 14:53

                               
登录/注册后可看大图

谢谢分享！！

客气了 ，应该的！！！楼主，是Matlab编程初学者吗？我这还有些资料，有时间给你发过去！！

Nicky_ONE

牛小贱 发表于 2013-12-13 12:44

                               
登录/注册后可看大图

客气了 ，应该的！！！楼主，是Matlab编程初学者吗？我这还有些资料，有时间给你发过去！！

奋斗小夏

牛小贱 发表于 2013-12-13 12:44
客气了 ，应该的！！！楼主，是Matlab编程初学者吗？我这还有些资料，有时间给你发过去！！

楼主，话说能不能给我发点资料  小弟最近也在编一个类似的程序，但是没有思路，不知从何下手！

牛小贱

奋斗小夏 发表于 2014-2-25 15:50
楼主，话说能不能给我发点资料  小弟最近也在编一个类似的程序，但是没有思路，不知从何下手！

我在本贴的49楼已经留下这个比较简单的单自由度振动动画的程序了：http://forum.vibunion.com/forum. ... 6579&fromuid=10
另外，我已经把【MATLAB编程(最好中文教程)】资料上传到新浪爱问资料共享上去了，你可以看一下，链接为http://ishare.iask.sina.com.cn/f/67354723.html（PS：你可以看一下Matlab在振动分析应用这方面的资料，对你会有帮助！但是，既然你也要编一个类似的程序，那我建议你还是把这个小程序吃透，关键是要理解里面每个语句的含义……）希望对你有所帮助！

奋斗小夏

牛小贱 发表于 2014-2-25 16:00
我在本贴的49楼已经留下这个比较简单的单自由度振动动画的程序了：http://forum.vibunion.com/forum.php? ...

好的 你辛苦了   非常地感谢你！

Zeng_0902

厉害

hanqzh

忒厉害了也

兵临神下