结构动力学-Duhamel积分法-matlab编程

%SDOF-简谐力-Duhamel clear;clc; %结构参数输入（质量、刚度、阻尼比） m=1;k=100;kesi=0.05; %外荷载输入（峰值、周期、时间间隔、时间步数、积分步数） Pmax=1;T=1;dt=0.02;n1=500; for i=1:n1 t=(i-1)*dt; T1(1,i)=t; Pt=Pmax*sin(2*pi*t/T); P(1,i)=Pt; end %求解无阻尼体系自振频率 wn=sqrt(k/m); %求解有阻尼体系自振频率 wD=wn*sqrt(1-kesi^2); %指定初始运动条件 u0=0; v0=0; aa0=0; %输入积分步数 n2=1000; %求解结构反应 u(1,1)=u0; v(1,1)=v0; aa(1,1)=aa0; T2(1,1)=0; for i=2:n2 T2(1,i)=(i-1)*dt; if i>n1 P(1,i)=0; end for j=2:i ut(1,j)=P(1,j)/(m*wD)*dt*exp(-kesi*wn*(T2(1,i)-T2(1,j)))*sin(wD*(T2(1,i)-T2(1,j))); end u(1,i)=sum(ut); end for i=1:n2-1 T3(1,i)=(i-1)*dt; v(1,i)=(u(1,i+1)-u(1,i))/dt; end for i=1:n2-2 T4(1,i)=(i-1)*dt; aa(1,i)=(v(1,i+1)-v(1,i))/dt; end %求解动力反应最大值与最小值 umax=max(u);umin=min(u); vmax=max(v);vmin=min(v); aamax=max(aa);aamin=min(aa); %绘制相对位移时程图 figure(1) plot(T2,u) xlabel('Time(s)'),ylabel('Displacement(m)') title('Displacement vs Time') saveas(gcf,'Displacement vs Time,png'); %绘制相对速度时程图 figure(2) plot(T3,v) xlabel('Time(s)'),ylabel('Velocity(m/s)') title('Velocity vs Time') saveas(gcf,'Velocityt vs Time,png'); %绘制相对加速度时程图 figure(3) plot(T4,aa) xlabel('Time(s)'),ylabel('Acceleration(m/s2)') title('Acceleration vs Time') saveas(gcf,'Acceleration vs Time,png');