求matlab编写的wilson或newmark法程序例子

zycshensave

小弟是新手，学习中遇到很多问题，尤其是自己编程的时候，发现不会的东西太多了，而且缺少实例和程序对照学习的资料，各位老大们，能分享一些关于结构振动分析的wilson 或newmark法的主子程序以及简单的梁结构实例吗，最好是matlab编的，万分感谢了，希望自己的学习能够得到大家的帮助！！谢谢

only8810205028

1. %function [d,v,a] = Wilson( K, M, C, f, d1, v1, dt, tend )
   
2. % 利用Wilson-theta (Newmark) 解析解法计算结构的动力响应
   
3.    K=[2,-1;-1,3];
   
4.   M=[4,0;0 8];
   
5.   C=0;
   
6. dt=0.1;% dt ----- 时间步长
   
7. tend=60;% tend --- 结束时间
   
8. w=600;
   
9. [n,n] = size( K ) ;
   
10. 
    
11. %Wilson算法
    
12. theta = 1.37 ;
    
13. tao = theta*dt ;
    
14. alpha0 = 6/tao^2 ;
    
15. alpha1 = 3/tao ;
    
16. alpha2 = 2*alpha1 ;
    
17. alpha3 = tao/2 ;
    
18. alpha4 = alpha0/theta ;
    
19. alpha5 = -alpha2/theta ;
    
20. alpha6 = 1-3/theta ;
    
21. alpha7 = dt/2 ;
    
22. alpha8 = dt^2/6 ;
    
23. K1 = K + alpha0*M + alpha1*C ;
    
24. d = zeros( n, floor(tend/dt) + 1 ) ;
    
25. v = zeros( n, floor(tend/dt) + 1 ) ;
    
26. a = zeros( n, floor(tend/dt) + 1 ) ;
    
27. f = zeros( n, floor(tend/dt) + 1 ) ;
    
28. d(1,1) = 1.2 ; %初始位移
    
29. d(2,1)=0;
    
30. v(:,1) = 0 ;%初始速度
    
31. f(:,1) =0 ;%初始载荷
    
32. a(:,1) = inv(M)*(f(:,1)-K*d(:,1)-C*v(:,1)) ;  %初始加速度
    
33. t=0:dt:tend;
    
34. for i=2:1:length(t)
    
35. %t(i) = (i-1)*dt ;
    
36. f(:,i)=0*100*sin(w*t(i));
    
37. ftheta = floor(theta) ;
    
38. %fq = f(i-1+ftheta-1)+ (theta-ftheta)*( f(i+ftheta-1) - f(i+ftheta-2) ) ;
    
39. fq=f(:,i-1)+ftheta*(f(:,i)-f(:,i-1));
    
40. f1 = fq + M*(alpha0*d(:,i-1)+alpha2*v(:,i-1)+2*a(:,i-1))+ C*(alpha1*d(:,i-1)+2*v(:,i-1)+alpha3*a(:,i-1)) ;
    
41. dq= inv(K1)*f1 ;
    
42. a(:,i) = alpha4*(dq-d(:,i-1)) + alpha5*v(:,i-1) + alpha6*a(:,i-1) ;
    
43. v(:,i) = v(:,i-1) + alpha7 * ( a(:,i) + a(:,i-1) ) ;
    
44. d(:,i) = d(:,i-1) + dt*v(:,i-1) + alpha8 * ( a(:,i)+2*a(:,i-1) ) ;
    
45. end
    
46. %解析解
    
47. 
    
48. for i=1:length(t)
    
49. x(i)=0.4*cos(0.5*t(i))+0.8*cos(1.581*0.5*t(i));
    
50. end
    
51. 
    
52. %Newmark算法
    
53. gama = 0.5 ;
    
54. beta = 0.25 ;
    
55. [n,n] = size( K ) ;
    
56. Nalpha0 = 1/beta/dt^2 ;
    
57. Nalpha1 = gama/beta/dt ;
    
58. Nalpha2 = 1/beta/dt ;
    
59. Nalpha3 = 1/2/beta - 1 ;
    
60. Nalpha4 = gama/beta - 1 ;
    
61. Nalpha5 = dt/2*(gama/beta-2) ;
    
62. Nalpha6 = dt*(1-gama) ;
    
63. Nalpha7 = gama*dt ;
    
64. NK1 = K + Nalpha0*M + Nalpha1*C ;
    
65. Nd = zeros( n, floor(tend/dt) + 1 ) ;
    
66. Nv = zeros( n, floor(tend/dt) + 1 ) ;
    
67. Na = zeros( n, floor(tend/dt) + 1 ) ;
    
68. Nd(1,1) = 1.2 ; %初始位
    
69. Nd(2,1)=0;
    
70. f(:,1) =0 ;%初始载荷
    
71. Na(:,1) = inv(M)*(f(:,1)-K*Nd(:,1)-C*Nv(:,1)) ;  %初始加速度
    
72. t=0:dt:tend;
    
73. for i=2:1:length(t)
    
74. f(:,i)=0*100*sin(w*t(i));
    
75. f2 = f(:,i) + M*(Nalpha0*Nd(:,i-1)+Nalpha2*Nv(:,i-1)+Nalpha3*Na(:,i-1))+ C*(Nalpha1*Nd(:,i-1)+Nalpha4*Nv(:,i-1)+Nalpha5*Na(:,i-1)) ;
    
76. Nd(:,i) = inv(NK1)*f2 ;
    
77. Na(:,i) = Nalpha0*(Nd(:,i)-Nd(:,i-1)) - Nalpha2*Nv(:,i-1) - Nalpha3*Na(:,i-1) ;
    
78. Nv(:,i) = Nv(:,i-1) + Nalpha6*Na(:,i-1) + Nalpha7*Na(:,i) ;
    
79. end
    
80. 
    
81. plot(t,d(1,:),'-b^',t,Nd(1,:),'g+',t,x,'r')
    
82. xlabel('t');
    
83. ylabel('u(t)');
    
84. title('位移与时间的关系');
    
85. 
    
86. 两自由度的三种解法
    
87. 
    
88. 
    
89. 
    

复制代码

刘德贵

好东西~

平行线不相交

好东西，谢谢

mxlzhenzhu

http://forum.vibunion.com/thread-138097-1-1.html

24hoursonline

厉害厉害 我都不知道怎么用 直接运行吗？，，

charlesxing

mxlzhenzhu 发表于 2016-4-8 17:33
http://forum.vibunion.com/thread-138097-1-1.html

这个可以的