[转帖]静止的单自由度Duhamel积分

happy

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. %%%%%%%本程序只适用于初始状态为静止的单自由度Duhamel积分，内部采用Simpson积分
   
3. %%%%%%%可计算无阻尼和阻尼强迫震动
   
4. %%%%%%%输入可为数组或函数荷载
   
5. %%%%%%%只能得出位移时程图
   
6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
7. 
   
8. clear all;
   
9. %w=30
   
10. w=input('输入自振频率 ω: ');
    
11. %n=10
    
12. n=input('输入荷载步数n ：');
    
13. %m=96.6/32.3
    
14. m=input('输入单元质量m ：');
    
15. %T=0.05
    
16. T=input('输入外荷载持续时间 T：');
    
17. %si=0.0
    
18. si=input('输入阻尼系数ξ：');
    
19. %k=2700
    
20. k=input('输入刚度系数k：');
    
21. deta=T/n;
    
22. wD=w*(1-si^2)^0.5;
    
23. G=deta/(m*wD)/3;
    
24. t1=[0:deta:T];
    
25. reply = input('输入荷载数组直接回车或敲击Y，输入函数荷载输入N： [Y]: ','s');
    
26. if isempty(reply)
    
27. reply = 'Y';
    
28. end
    
29. if reply=='Y'
    
30. p2=input('输入荷载数组并用[]抱住 为n＋1个元素：例如[0 19.3200 38.6400 57.9600 77.2800 96.6000 77.2800 57.9600 38.6400 19.3200 0]:\n');
    
31. 
    
32. %p2=[0:19.32:96.6,77.28:-19.32:0]
    
33. exp1=exp(-2*si*w*deta);
    
34. exp2=4*exp(-si*w*deta);
    
35. A(1)=p2(1)*cos(wD*0);
    
36. for i=2:n/2+1
    
37. A(i)=(A(i-1)+p2(2*i-3)*cos(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*cos(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*cos(wD*2*(i-1)*deta);
    
38. end
    
39. 
    
40. B(1)=p2(1)*sin(wD*0);
    
41. for i=2:n/2+1
    
42. B(i)=(B(i-1)+p2(2*i-3)*sin(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*sin(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*sin(wD*2*(i-1)*deta);
    
43. end
    
44. 
    
45. t1=[0:2*deta:T];
    
46. disp('位移随时间的变化');
    
47. V=G*(A.*sin(wD*t1)-B.*cos(wD*t1))
    
48. disp('弹性力随时间的变化');
    
49. f=k*V
    
50. for i=1:n/2+1
    
51. p22(i)=p2(2*i-1);
    
52. end
    
53. disp('加速度随时间的变化');
    
54. v11=(p22-k*V)/m
    
55. subplot(1,2,1),plot(t1,V),title('位移'),grid on,
    
56. 
    
57. 
    
58. t2=[T:deta:T*5];
    
59. sw=sin(w*t2);
    
60. cw=cos(w*t2);
    
61. disp('扩大四倍后位移随时间的变化');
    
62. y=G*A(n/2+1)*sw-G*B(n/2+1)*cw
    
63. subplot(1,2,2),
    
64. 
    
65. plot(t1,V,'r') ,hold on
    
66. plot(t2,y,'r');
    
67. title('加载时间增加四倍时的位移位移时程曲线')
    
68. grid on;
    
69. else
    
70. 
    
71. syms tao t;
    
72. y=input('函数荷载(关于未知数tao)：');
    
73. %y=-10000/8*tao+100;
    
74. v1=1/m/wD*int(y*exp(-si*w*(t-tao))*sin(wD*(t-tao)),tao,0,t);
    
75. v01=diff(v1,t);
    
76. v02=diff(v01,t)
    
77. t1=[0:T/n:T];
    
78. for i=1:n+1
    
79. v2(i)=subs(v1,t,t1(i));
    
80. v011(i)=subs(v01,t,t1(i));
    
81. v022(i)=subs(v02,t,t1(i));
    
82. end
    
83. disp('荷载作用时间内位移变化');
    
84. v2
    
85. disp('荷载作用时间内位移最大值');
    
86. max(v2)
    
87. disp('荷载作用时间内速度变化');
    
88. v011
    
89. disp('荷载作用时间内速度最大值');
    
90. max(v011)
    
91. disp('荷载作用时间内加速度变化');
    
92. v022
    
93. disp('荷载作用时间内加速度最大值');
    
94. max(v022)
    
95. subplot(2,3,1),plot(t1,v2),grid on,title('荷载作用时间内位移时程曲线');
    
96. subplot(2,3,2),plot(t1,v011),grid on,title('荷载作用时间内速度时程曲线');
    
97. subplot(2,3,3),plot(t1,v022),grid on,title('荷载作用时间内加速度时程曲线')
    
98. A=1/m/wD*int(y*exp(si*w*tao-si*w*t)*cos(wD*tao),tao,0,T);
    
99. B=1/m/wD*int(y*exp(si*w*tao-si*w*t)*sin(wD*tao),tao,0,T);
    
100. 
     
101. A1=subs(A,T);
     
102. B1=subs(B,T);
     
103. V33=A1*sin(wD*t)-B1*cos(wD*t);
     
104. V033=diff(V33,t);
     
105. V0033=diff(V033,t);
     
106. 
     
107. t2=[T:T/n:5*T];
     
108. v33=A1*sin(wD*t2)-B1*cos(wD*t2);
     
109. for i=1:4*n+1
     
110. V03(i)=subs(V033,t,t2(i));
     
111. V003(i)=subs(V0033,t,t2(i));
     
112. end
     
113. disp('扩大四倍后位移随时间的变化');
     
114. v33
     
115. disp('扩大四倍后位移最大值');
     
116. max(v33)
     
117. disp('扩大四倍后速度随时间的变化');
     
118. V03
     
119. disp('扩大四倍后速度最大值');
     
120. max(V03)
     
121. disp('扩大四倍后加速度随时间的变化');
     
122. V003
     
123. disp('扩大四倍后加速度最大值');
     
124. max(V003)
     
125. subplot(2,3,4),plot(t1,v2),hold on,title('加载时间增加四倍时的位移时程曲线'),plot(t2,v33),grid on;
     
126. subplot(2,3,5),plot(t1,v011),hold on,title('加载时间增加四倍时的速度时程曲线'),plot(t2,V03),grid on;
     
127. subplot(2,3,6),plot(t1,v022),hold on,title('加载时间增加四倍时的加速度时程曲线'),plot(t2,V003),grid on;
     
128. end

复制代码

来自钢结构

[ 本帖最后由 suffer 于 2006-10-9 19:39 编辑 ]

assist

谢谢，很好用，只是我想要个多自由度的。