newmarkb.rar_K._结构动力响应程序_荷载识别

%%modified by Tang 2007-05-14 at Tongji University function [X,Xd,Xdd,t] = newmarkb(M,C,K,f,dt,gamma,beta,Xi,Xdi) %================================================================ %function [X,t] = newmark(M,C,K,f,dt,gamma,beta,Xi,Xdi) % % (c) Jerome Peter Lynch, (all rights reserved) % Stanford University % February 14, 2001 % % 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 newmark alpha-beta method.okok.org The newmark alpha-beta % 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= Newmark Const (1/6 or 1/4 usually) % gam = Newmark Const (1/2) % 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); % Check Input Excitation % ====================== if (fdimc == n) f=f'; end; m = size(f,2); % Coefficients % ============ c0 = 1/(beta*dt*dt) ; c1 = gamma/(beta*dt) ; c2 = 1/(beta*dt) ; c3 = 1/(beta*2) - 1 ; c4 = gamma/beta - 1 ; c5 = 0.5*dt*(gamma/beta - 2 ) ; c6 = dt*(1 - gamma ) ; c7 = dt* gamma; % Initialize Stiffness Matricies %============================== Keff = c0*M + c1*C + K ; Kinv = inv(Keff) ; %Initial Acceleration %==================== R0 = f(:,1); Xddi= inv(M)*(R0 - K*Xi - C*Xdi); %Perform First Step %================== f(:,1) = f(:,1) + M*(c0*Xi+c2*Xdi+c3*Xddi) ... +C*(c1*Xi+c4*Xdi+c5*Xddi) ; X(:,1) = Kinv*f(:,1) ; Xdd(:,1)= c0*(X(:,1)-Xi) - c2*Xdi - c3*Xddi ; Xd(:,1) = Xdi + c6*Xddi + c7*Xdd(:,1) ; %Perform Subsequent Steps % ======================== for i=1:size(f,2)-1; f(:,i+1) = f(:,i+1) + M * ... (c0*X(:,i)+c2*Xd(:,i)+c3*Xdd(:,i)) ... + C*(c1*X(:,i)+c4*Xd(:,i)+c5*Xdd(:,i)) ; X(:,i+1) = Kinv*f(:,i+1) ; Xdd(:,i+1)= c0*(X(:,i+1)-X(:,i))-c2*Xd(:,i)-c3*Xdd(:,i) ; Xd(:,i+1) = Xd(:,i)+c6*Xdd(:,i)+c7*Xdd(:,i+1) ; end; % Strip Off Padded Response Zeros % =============================== % X = X'; % Generate the Time Vector % ======================== for i=0:1:size(X,2)-1 t(i+1,1) = i*dt; end;