单自由度体系杜哈梅积分.doc

function y=kst(t0,t1,t2,ts,m,b0,b1,w0,c) t0=input(请输入起始时间:t0= );t1=input(请输入荷载消失时间:t1= );t2=input(请输入想要的时间:t2= ); ts=input(请输入时间步长:ts= ); m=input(请输入质量:m= ) ;b0=input(请输入荷载截距:b0= );b1=input(荷载消失时的荷载:b1= );k=input(请输入刚度:k= ) ; c=input(请输入阻尼比:c= ); w0=sqrt(k/m);w1=w0*sqrt(1-c^2); t=t0:ts:t2; for i=1:(length(t)) x=linspace(t(1),t(length(t))) p=interp1([t0 t1],[b0 b1],t); p(find(isnan(p)==1)) = 0; px=linspace(p(1),p(length(t))); a=px.*exp(c*w0*x).*cos(w1*x); A=trapz(x,a); b=px.*exp(c*w0*x).*sin(w1*x); B=trapz(x,b); y=exp(-c*w0*t).*(A.*sin(w1*t)-B.*cos(w1*t))./(m*w1) v=diff(y) a0=diff(y,2) end ymax=max(y) figure plot(t,y); 此程序为复合梯形法计算冲击荷载作用下的杜哈梅积分。 以P(t)=-1250000*(t+0.08)的冲击荷载为例，质量：m=6.4；阻尼比c=0.05；刚度：k=34847.77 N/m.将参数输入程序得到以下结果： 冲击荷载作用下的体系的速度图 冲击荷载作用下的体系加速度图