%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % CE 573 Structural Dynamics % % % % SIMMUL ASSMT 1 Shear Building Part-I % % % % Written by: A.A. AlNaseem, 11/21/11 % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all % DECLARE ASSUMPTIONS FOR MODEL disp('We begin with our assumptions:') disp('1.Floors are rigid.') disp('2.Columns do not change in length.') disp('3.Column properties of floors with splices will be weight averaged.') disp('4.This is a nonlinear system and masses are lumped at nodes.') % INPUT BUILDING PROPERTIES E = 14830; %Ksi I = [8210 7701.4 7190 6317.4 5440 4642.2 3840 3620.6 3400]';%inches^4 hg = 18*12; %inches h = 12*12; %inches zeta = 0.02*[1 1 1 1 1 1 1 1 1]';%damping ratio for all modes num_floors = 9; % BUILD MASS, STIFFNESS, DAMPING MATRICES , AND OTHER CALCULATIONS M = diag([69/12 67.7/12 67.7/12 67.7/12 67.7/12 67.7/12 67.7/12 67.7/12 73.2/12]);% Kips-sec^2/inches K = E*[(72*I(1, 1)/hg^3)+(72*I(2, 1)/h^3) -(72*I(2, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(2, 1)/h^3) (72*I(2, 1)/h^3)+(72*I(3, 1)/h^3) -(72*I(3, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(3, 1)/h^3) (72*I(3, 1)/h^3)+(72*I(4, 1)/h^3) -(72*I(4, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(4, 1)/h^3) (72*I(4, 1)/h^3)+(72*I(5, 1)/h^3) -(72*I(5, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(5, 1)/h^3) (72*I(5, 1)/h^3)+(72*I(6, 1)/h^3) -(72*I(6, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(6, 1)/h^3) (72*I(6, 1)/h^3)+(72*I(7, 1)/h^3) -(72*I(7, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(7, 1)/h^3) (72*I(7, 1)/h^3)+(72*I(8, 1)/h^3) -(72*I(8, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(8, 1)/h^3) (72*I(8, 1)/h^3)+(72*I(9, 1)/h^3) -(72*I(9, 1)/h^3) 0 0 0 0 0 0 0 -(72*I(9, 1)/h^3) (72*I(9, 1)/h^3)]; % BUILD C MATRIX if (size(M,2)==1)&&(size(M,1)==1), MM = M*eye(ndof); elseif (size(M,2)==1)||(size(M,1)==1), MM = diag(M); else MM = M; end if (size(K,2)==1)&&(size(K,1)==1), KK = diag(K*ones(1,num_floors)) + diag([K*ones(1,num_floors-1) 0]) + diag(-K*ones(1,num_floors-1),-1) + diag(-K*ones(1,num_floors-1),1); elseif (size(K,2)==1)||(size(K,1)==1), KK = diag(K) + diag([K(2:end) 0]) - diag(K(2:end),-1) - diag(K(2:end),1); else KK = K; end %FINDING EIGENVALUES AND EIGENVECTORS [EIGVEC,EIGVAL] = eig(KK/MM); EIGVAL = sqrt(diag(EIGVAL)); [temp,ind] = sort(abs(EIGVAL)); EIGVAL = EIGVAL(ind)/2/pi; EIGVEC = EIGVEC(:,ind); %CLACULATE DAMPING MATRIX TwoZW = 2*diag(zeta)*diag(2*pi*EIGVAL); CC = MM*EIGVEC*TwoZW*EIGVEC'; %CALCULATING LOAD FACTOR GG = -MM * [1;1;1;1;1;1;1;1;1]; %DSIPLAYING RESULTS disp('Part(a):') disp('1.Mass Matrix:') disp(MM) disp('2.Stiffness Matrix:') disp(KK) disp('3.Damping Matrix:') disp(CC) disp('4.Load Factor Vector:') disp(GG) disp('Part(b):') disp('1.Natural Frequencies:') disp(EIGVAL) disp('1.Natural Mode Matrix:') a = 1; PHI = [9,9]; Demo = [1:9; 1:9; 1:9; 1:9; 1:9; 1:9; 1:9; 1:9; 1:9]'; while a<10 PHI (1:9, a) = ((EIGVEC(1:9,a)/((max(abs(EIGVEC(1:9,a))))))); a = a+1; end disp (PHI) % PLOT MODE SHAPES fig1 = figure(1); n = 1; set(fig1,'position',[10 60 600 600],'name','Mode Shapes: 1-3'); while n<=3 subplot(1,3,n);plot([0 (PHI(1,n))],[0 Demo(1,n)],'r',[0 (PHI(1,n))],[0 Demo(1,n)],'o'), hold on plot(PHI(1:9,n),Demo(1:9,n),'o'),hold on plot((PHI(1:9,n)),Demo(1:9,n),'r'); set(subplot(1,3,n),'XTick',[-1:0.5:1]); axis([-1 1 0 10]); title(['Mode',int2str(n)]);xlabel(['Freq = ',num2str(real(sqrt(EIGVAL(n))))]); if n == 1 ylabel(['Floor']); end n = n + 1; end fig2 = figure(2); n = 4; i = 1; set(fig2,'position',[630 60 600 600],'name','Mode Shapes: 4-6'); while n<=6 subplot(1,3,i);plot([0 (PHI(1,n))],[0 Demo(1,n)],'r',[0 (PHI(1,n))],[0 Demo(1,n)],'o'), hold on plot(PHI(1:9,n),Demo(1:9,n),'o'),hold on plot((PHI(1:9,n)),Demo(1:9,n),'r'); set(subplot(1,3,i),'XTick',[-1:0.5:1]); axis([-1 1 0 10]); title(['Mode',int2str(n)]);xlabel(['Freq = ',num2str(real(sqrt(EIGVAL(n))))]); if n == 1 ylabel(['Floor']); end n = n + 1; i = i + 1; end fig3 = figure(3); n = 7; i = 1; set(fig3,'position',[1250 60 600 600],'name','Mode Shapes: 7-9'); while n<=9 subplot(1,3,i);plot([0 (PHI(1,n))],[0 Demo(1,n)],'r',[0 (PHI(1,n))],[0 Demo(1,n)],'o'), hold on plot(PHI(1:9,n),Demo(1:9,n),'o'),hold on plot((PHI(1:9,n)),Demo(1:9,n),'r'); set(subplot(1,3,i),'XTick',[-1:0.5:1]); axis([-1 1 0 10]); title(['Mode',int2str(n)]);xlabel(['Freq = ',num2str(real(sqrt(EIGVAL(n))))]); if n == 1 ylabel(['Floor']); end n = n + 1; i = i + 1; end disp('SEE FIGURES 1-3 FOR MODE SHAPES - HIT RETURN TO CONTINUE'); disp('Part(c):') disp('The EOM for each mode are in the following format:') disp('Mn'); disp('q``'); disp('Cn'); disp('q`'); disp('Kn'); disp('q'); disp('='); disp('gamma*r'); disp('zg``'); l = 1; while l<=9 disp(l); disp(PHI(1:9,l)'*MM*PHI(1:9,l)); disp('q``'); disp(PHI(1:9,l)'*CC*PHI(1:9,l)); disp('q`'); disp(PHI(1:9,l)'*KK*PHI(1:9,l)); disp('q'); disp('='); disp(PHI(1:9,l)'*GG); disp('zg``'); l = l + 1; end disp('**The most significant mode is mode 1 since it gives the') disp('least naturl frequency ehich is the fundamental freaquecny**'); disp('Part(d): See simulink window'); LAM = [eye(num_floors)-diag([ones(1,num_floors-1)],1)]; GAM = []; t=0.1; f=1; while f<=9 GAM(f,1) = 2*pi*71.17*sin(2*pi*t)*GG(f,1); % ground acceleration f=f+1 t= t+0.1 end % BUILD STATE SPACE MODEL invMM = inv(MM); AA = [zeros(num_floors) eye(num_floors);-invMM*KK -invMM*CC]; BB = [[zeros(num_floors,1);GAM] [zeros(num_floors);invMM*LAM]]; CC = [eye(num_floors) zeros(num_floors); eye(num_floors)-diag([ones(1,num_floors-1)],-1) zeros(num_floors) zeros(num_floors) eye(num_floors); -invMM*KK -invMM*CC]; DD = [[zeros(num_floors*4,1)] [zeros(num_floors*3,num_floors);invMM*LAM]]; % create force input dt = 0.1; Tmax = 10; time = [0:dt:Tmax]'; %force = [ones(size(time))]; % step input force = [time]; % shape the earthquake force = 2*pi*71.17*sin(2*pi*force); %return sim('Sim_PatI',[0 Tmax]); % plot responses figure(1) subplot(311) plot(tout,yout(:,1)) ylabel('displ (m)') subplot(312) plot(tout,yout(:,2)) ylabel('vel (m/sec)') subplot(313) plot(tout,yout(:,3)); xlabel('time (sec)') ylabel('accel (m/sec^2)')