||
function [X,Xd,Xdd,t] = cendiff(M,C,K,f,dt,Xi,Xdi) %================================================================ %function [X,t] = newmark(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 central difference method. The central difference % 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 % 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(:,2)=Xi; Xd(:,1)=Xdi; % Check Input Excitation % ====================== if(fdimc==n) f=f'; end m=size(f,2); % Coefficients % ============ c0 = 1/(dt*dt) ; c1 = 1/(2*dt) ; c2 = 2*c0 ; %Initialize Mass Matricies % ============================== Meff = c0*M + c1*C ; Minv = inv(Meff) ; % Initial Acceleration % ==================== R0=f(:,1); Xdd(:,1)= inv(M)*(f - K*Xi - C*Xdi); % Perform First Step % ================== X(:,1) = X(:,2) - 1/(2*c1)*Xd(:,1) + 1/(2*c0)*Xdd(:,1) ; f(:,1) = f(:,1) - (K - c2*M)*X(:,2) - (c0*M - c1*C)*X(:,1); X(:,3) = Minv*f(:,1); Xd(:,2) = c1*(X(:,3) - X(:,1)); Xdd(:,2) = c0*(X(:,3) - 2*X(:,2) + X(:,1)); %Perform Subsequent Steps % ======================== for i=1:size(f,2)-1; X(:,i+1) = X(:,i+2) - 1/(2*c1)*Xd(:,i+1) + 1/(2*c0)*Xdd(i+1) ; f(:,i+1) = f(:,i+1) - (K - c2*M)*X(:,i+2) - (c0*M - c1*C)*X(:,i+1); X(:,i+3) = Minv*f(:,i+1); Xd(:,i+2) = c1*(X(:,i+3) - X(:,i+1)); Xdd(:,i+2) = c0*(X(:,i+3) - 2*X(:,i+2) + X(:,i+1)); 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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 13:20
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社