多自由度振动程序

西风独自凉

例题见张亚辉结构动力学基础


function newmark(t_max,dt)
%t_max为持续时间，dt为时间步长
M=1*[1 0 0;0 1 0;0 0 1];
K=1*[3 -1 0;-1 2 -1;0 -1 1];
C=0.015*M+0.02*K;%c=0.015*[0 0 0;0 0 0;0 0 0];
u=[0 0 0]';
v=[0 0 0]';
a=[0 0 0]';
t(1)=0;               %时间
x(:,1)=u;             %位移
x1(:,1)=v;            %速度
x2(:,1)=a;            %加速度
%newmark参数
gama=0.5;
delta=0.25;
a0=1/(delta*dt^2);
a1=gama/(delta*dt);
a2=1/(delta*dt);
a3=1/(2*delta)-1;
a4=gama/delta-1;
a5=dt*(gama/(2*delta)-1);
a6=dt*(1-gama);
a7=gama*dt;
%等效刚度矩阵
K1=K+a0*M+a1*C;
i=1;
t(1)=0;
while t(i)<t_max
      Q=-1*[sin(t(i)) sin(t(i)) sin(t(i))]';
      q(:,i+1)=Q+M*(a0*x(:,i)+a2*x1(:,i)+a3*x2(:,i))+C*(a1*x(:,i)+a4*x1(:,i)+a5*x2(:,i));
      x(:,i+1)=inv(K1)*q(:,i+1);
      x2(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*x1(:,i)-a3*x2(:,i);
      x1(:,i+1)=a1*(x(:,i+1)-x(:,i))-a4*x1(:,i)-a5*x2(:,i);
      i=i+1;
      t(i)=t(i-1)+dt;
end  
x1=x(1,:)';
x2=x(2,:)';
x3=x(3,:)';
xx1=x1(1:t_max/dt)';
xx2=x2(1:t_max/dt)';
xx3=x3(1:t_max/dt)';
t=dt:dt:t_max;
plot(t,xx1,'-',t,xx2,'-.',t,xx3,'--')
xlabel('t[sec]')
ylabel('x(m)')
legend('x1(t)','x2(t)','x3(t)')
grid




输入：newmark(30,0.01)

yuanmeng2015

感谢     

小星星

程序很有用，谢谢楼主无私奉献！！！