用newmark法算多自由度系统动力响应有问题

martin2013

给出的方程式

这个是Fmr的有关equation

所模拟的建筑物

根据上面的方程式和模拟图，怎么通过newmark算其动力响应？
求帮助，我看了很多有关的例子，但是还是不能解决，这个方程式的右边，f(t)并没有给出，但是文章说的gaussian white-noise  ground acceleration with spectral density S0=100 cm^2 .我不知道怎么确定这个excitiation 和fmr(t).求高手解答。

我在赋予一个例子，这个是书上的。

% 利用Newmark 法计算结构的动力响应 function [d,v,a] = Newmark( K, M, C, f, d1, v1, dt, tend )

% 输入参数

  %    K  -----  刚度矩阵

  %    M  -----  质量矩阵

  %    C  -----  阻尼矩阵

  %    d1 -----  初始位移

  %    v1 -----  初始速度

  %    dt -----  时间步长

  %    tend ---  结束时间

  % 返回值

  %    d  -----  位移

  %    v  -----  速度

  %    a  -----  加速度

    gama = 0.5 ; beta = 0.25 ;

    [n,n] = size( K ) ;

alpha0 = 1/beta/dt^2 ; alpha1 = gama/beta/dt ; alpha2 = 1/beta/dt ; alpha3 = 1/2/beta - 1 ;

alpha4 = gama/beta - 1 ;  alpha5 = dt/2*(gama/beta-2) ;  alpha6 = dt*(1-gama) ;  alpha7 = gama*dt ;

      K1 = K + alpha0*M + alpha1*C ;

      d = zeros( n, floor(tend/dt) + 1 ) ;

      v = zeros( n, floor(tend/dt) + 1 ) ;

      a = zeros( n, floor(tend/dt) + 1 ) ;

      d(:,1) = d1 ;

      v(:,1) = v1 ;

      a(:,1) = M\(f(:,1)-K*d1-C*v1) ;

      for i=2:1:floor(tend/dt) + 1;

          t = (i-1)*dt ;

  f1 = f(:,i) + M*(alpha0*d(:,i-1)+alpha2*v(:,i-1)+alpha3*a(:,i-1)+ C*(alpha1*d(:,i-1)+alpha4*v(:,i-1)+alpha5*a(:,i-1)) ;

          d(:,i) = K1\f1 ;

          a(:,i) = alpha0*(d(:,i)-d(:,i-1)) - alpha2*v(:,i-1) - alpha3*a(:,i-1) ;

          v(:,i) = v(:,i-1) + alpha6*a(:,i-1) + alpha7*a(:,i) ;

      end

那我该怎么仿写这个呢？我自己写的我plot不出来。。求大神解救。。。。。。。。。。。。

martin2013

我再附上我自己的编辑的code，其中那个f荷载的力没弄清楚，还期待人帮我修改。出主意，还有fmr的也没考虑进去，因为不知道怎么附加上去


k=883.65;
a=[-k,-k,-k,-k,-k,-k,-k,-k,-k];
b=[2*k,2*k,2*k,2*k,2*k,2*k,2*k,2*k,2*k,k];
K=diag(b)+diag(a,1)+diag(a,-1);
M=diag([50,50,50,50,50,50,50,50,50,50]);
C=[0.094*M+0.0008*K];
d=ones(10,1);
T1=M*d;

dt=0.05;
T=1000;
beta=0.25;
gamma=0.5;
a0=1/(beta*dt*dt);
a1=gamma/(beta*dt);
a2=1/(beta*dt);
a3=1/(beta*2)-1;
a4=gamma/beta-1;
a5=0.5*dt*(gamma/beta-2);
a6=dt*(1-gamma);
a7=dt*gamma;
K1=K+a1*C+a0*M;

Y=randn([1,T/dt-1]);
x0=zeros(10,1+T/dt);
x1=zeros(10,1+T/dt);
x2=zeros(10,1+T/dt);
x0(:,1)=x0;
x1(:,1)=x1;
x2(:,1)=inv(M)*(f-K*x0-C*x1);
for i=2:1:1+T/dt;
    t=(i-1)*dt;
    f(:,i)=T1*Y(1,i);
    f1=f(:,1)+M*(a0*x0(:,i-1)+a2*x1(:,i-1)+a3*x2(:,i-1))+C*(a1*x0(:,i-1)+a4*x1(:,i-1s)+a5*x2(:,i-1));
    x0(:,i)=inv(K1)*f1;
    x2(:,i)=a0*(x(:,i)-x(:,i-1)-a2*x1(:,i-1)-a3*x2(:,i-1));
    x1(:,i)=x1(:i-1)+a6*x2(:,i-1)+a7*x2(:,i);
end

martin2013

求帮忙求帮忙求帮忙

hnczf739892267

楼主，你的问题解决了没有啊，我也遇到这个问题

luotao1014