wilson.m

动力学教程中的数值方法

function [X,Xd,Xdd,t] = wilson(M,C,K,f,dt,beta,Xi,Xdi) 
%================================================================ 
%function [X,t] = wilson(M,C,K,f,dt,gamma,beta,Xi,Xdi) 
% 
%  (c) Ajax, (all rights reserved) 
%  Central South University 
%  October 13, 2006 
% 
%  This function is intended to perform the numerical 
%  integration of a structural system subjected to an 
%  external dynamic excitation such as a wind or earthquake. 
%  The structural model is assumed to be a lumped mass shear 
%  model.The integration scheme utilized for this analysis 
%  is the wilson-θ method. The wilson-θ 
%  method is an implicit time steping scheme so stability of 
%  the system need not be considered.
% 
%  Input Variables: 
%  [M] = Mass Matrix (nxn) 
%  [C] = Damping Matrix (nxn) 
%  [K] = Stiffness Matrix (nxn) 
%  {f} = Excitation Vector (mx1) 
%  dt = Time Stepping Increment 
%  beta= Wilson Const (1.4 usually)  
%  Xi = Initial Displacement Vector (nx1) 
%  Xdi = Initial Velocity Vector (nx1) 
% 
%  Output Variables: 
%  {t} = Time Vector (mx1) 
%  [X] = Response Matrix (mxn) 
%================================================================ 

  n = size(M,1); 
  fdimc = size(f,2); 
  fdimr = size(f,1); 
  X(:,1) = Xi;
  Xd(:,1) = Xdi;
  
%  Check Input Excitation 
%  ====================== 
   if(fdimc==n)
       f=f';
   end
   m=size(f,2);
   
%  Coefficients 
%  ============ 
        c0 = 6/(beta*dt*dt) ; 
        c1 = 3/(beta*dt) ; 
        c2 = 2*c1 ; 
        c3 = beta*dt/2; 
        c4 = c0/beta ; 
        c5 = -c2/beta ; 
        c6 = 1 - 3/beta ; 
        c7 = dt/2 ; 
        c8= dt*dt/6 ;

%Initialize Stiffness Matricies 
% ============================== 
        Keff = c0*M + c1*C + K ; 
        Kinv = inv(Keff) ;

         
% Initial Acceleration
% ==================== 
        R0=f(:,1);
        Xdd(:,1) = inv(M)*(f - K*Xi - C*Xdi)

% Perform First Step 
% ==================
        R = f(:,1) + beta*(f(:,1)-f(:,1)) + M*(c0*X(:,1)+c2*Xd(:,1)+2*Xdd(:,1)) ... 
                 +C*(c1*X(:,1)+2*Xd(:,1)+c3*Xdd(:,1)) ; 
        Xb  = Kinv*R ; 
        Xdd(:,2)= c4*(Xb-X(:,1)) + c5*Xd(:,1) + c6*Xdd(:,1) ; 
        Xd(:,2) = Xd(:,1) + c7*(Xb+Xd(:,1)) ; 
        X(:,2) = X(:,1) + dt*Xd(:,1) + c8*(Xb+2*Xdd(:,1));
        f(:,2) = M*Xdd(:,2) + C*Xd(:,2) + K*X(:,2);

%Perform Subsequent Steps 
% ======================== 
        for i=1:size(f,2)-1; 
           R = f(:,1) + beta*(f(:,2)-f(:,1)) + M*(c0*X(:,1)+c2*Xd(:,1)+2*Xdd(:,1)) ... 
                 +C*(c1*X(:,1)+2*Xd(:,1)+c3*Xdd(:,1)) ; 
           Xb  = Kinv*R ; 
           Xdd(:,2)= c4*(Xb-X(:,1)) + c5*Xd(:,1) + c6*Xdd(:,1) ; 
           Xd(:,2) = Xd(:,1) + c7*(Xb+Xd(:,1)) ; 
           X(:,2) = X(:,1) + dt*Xd(:,1) + c8*(Xb+2*Xdd(:,1));
           f(:,2) = M*Xdd(:,2) + C*Xd(:,2) + K*X(:,2); 
           X(:,1)=X(:,2);
           Xd(:,1)=Xd(:,2);
           Xdd(:,1)=Xdd(:,2);
        end; 

         
%  Strip Off Padded Response Zeros 
%  =============================== 
  X = X'
  Xd=Xd'
  Xdd=Xdd'
  
%  Generate the Time Vector 
%  ======================== 
  for i=0:1:size(X,1)-1 
    t(i+1,1) = i*dt; 
  end; 
  t