damped2dof.m

用线性加速度法来解双自由度体系

  clear;clc
  format long;   %数据格式
  d_t=0.002;              %步长  
  m1=806;m2=25;
  f1=8.02;f2=10.43;
  c1=4*pi*m1*f1*0.011;
  c2=4*pi*m2*f2*0.5;
  k1=4*pi^2*m1*f1^2;
  k2=4*pi^2*m2*f2^2;
  M=[m1,0;0,m2];           %质量矩阵
  C=[c1+c2,-c2;-c2,c2];    %阻尼矩阵
  K=[k1+k2,-k2;-k2,k2];    %刚度矩阵

  t=0:d_t:4;               %分划 
  y(1,1)=input('请输入楼板位移初值y(0)=');    %初始位移
  y(2,1)=input('请输入人体位移初值y(0)=');    %初始位移
  y1(1,1)=input('请输入楼板速度初值v(0)=');   %初始速度
  y1(2,1)=input('请输入人体速度初值v(0)=');   %初始速度
  for i=1:length(t)-1                        %循环
      Fi(:,i)=0-K*y(:,i)-C*y1(:,i);
      y2(:,i)=inv(M)*Fi(:,i);
      K1=K+(6/(d_t)^2)*M+(3/d_t)*C;
      d_p1(:,i)=0+M*((6/d_t)*y1(:,i)+3*y2(:,i))+C*(3*y1(:,i)+(d_t/2)*y2(:,i));
      d_y(:,i)=inv(K1)*d_p1(:,i);
      d_y1(:,i)=(3/d_t)*(d_y(:,i)-d_t*y1(:,i)-(d_t^2/6)*y2(:,i));
      y(:,i+1)=y(:,i)+d_y(:,i);
      y1(:,i+1)=y1(:,i)+d_y1(:,i);
      i=i+1;
  end
  Fi(:,length(t))=0-K*y(:,length(t))-C*y1(:,length(t));
  y2(:,length(t))=inv(M)*Fi(:,length(t));
  
  figure(1)
  plot(t,y(1,:),'r-');  %画图
  title('位移曲线')
  legend('t-y')
  xlabel('t')
  ylabel('y')
  grid on;
  pause;
  
  figure(2)
  plot(t,y1(1,:),'b-');
  title('速度曲线')
  legend('t-v')
  xlabel('t')
  ylabel('v')
  grid on;
  pause
  
  figure(3)
  plot(t,y2(1,:),'m-')
  title('加速度曲线')
  legend('t-a')
  xlabel('t')
  ylabel('a')
  grid on;