MATLAB Examples

Time History Analysis of Multi-story Building by Mode Superposition Method

In this function to solve Equation of Motions for different modes, Newmark Linear Method has been used. By modifying inputs F & ft, you can analyze the structute for either earthquake (base-seismic-excitation) or time-dependent load applied to different stories. That's your responsibility to organize phi & omega2 such f1<f2<f3<... Preferred Unit is SI e.g. K(N/mm),M(N.s^2/mm),F(N),g=9810 mm/s^2 & Ag(unitless to g). However, you are free to choose your Unit.

Contents

@ Mostafa Tazarv, Carleton University, March 2011

NM Function Help:

% INPUTS:
% n:      Number of Stories or generally, Number fo DOFs
% nom:    Number of Modes that you want to consider in analysis; nom <= n
% dt:     Time Interval (Sampling Time)
% F:      Spatial Load Distribution Vector; in the case of base seismic
%         exitation F=-g*M*[1 1 1 ... 1]'; (A vector)
% ft:     Time-Variant Component of Load (A Vector); in the
%         case of base seismic exitation ft=Ag (Ground Acceleration in g)
% M:      Mass Matrix of Stories (n x n)
% K:      Stiffness Matrix of Stories (n x n)
% zet:    Damping Ratio in percent (%); e.g. 5
% omega2: Eigen-Values Matrix (nom x nom); e.g. [phi,omega2]=eig(K,M);
% phi:    Eigen-Vectors (nom x nom);  e.g. [phi,omega2]=eig(K,M);

% OUTPUTS (Story Response):
% disp1:  Displacement Time History of Stories
% EF:     Elastic Force Time History of Stories

Newmark Function; NM

function [disp1,EF]=NM(n,nom,dt,F,ft,M,K,zet,omega2,phi)

By Newmark Linear Method

Mn=phi'*M*phi;
Kn=phi'*K*phi;
Ln=zeros(nom);
 for i=1:nom
     Ln(i,i)=(phi(:,i))'*F/Mn(i,i);  % Load Contribution Facor
 end
u=zeros(length(ft),nom);
v=zeros(length(ft),nom);
ac=zeros(length(ft),nom);
ft(end+1)=0;
for j=1:nom                          % equation of motion(Newmark linear method)
    m=Mn(j,j);
    k=Kn(j,j);
    c=2*m*(omega2(j,j))^.5*zet/100;
    K=k+3*c/dt+6*m/(dt)^2;
    a=6*m/dt+3*c;
    b=3*m+dt*c/2;
  for i=1:length(u)-1
     u(1,j)=0;                      %initial conditions
     v(1,j)=0;
     ac(1,j)=0;
     df=Ln(j,j)*Mn(j,j)*(ft(i+1)-ft(i))+a*v(i,j)+b*ac(i,j);  % delta Force
     du=df/K;
     dv=3*du/dt-3*v(i,j)-dt*ac(i,j)/2;
     dac=6*(du-dt*v(i,j))/(dt)^2-3*ac(i,j);
     u(i+1,j)=u(i,j)+du;
     v(i+1,j)=v(i,j)+dv;
     ac(i+1,j)=ac(i,j)+dac;
  end
end
ft(end)=[];
Input argument "phi" is undefined.

Error in ==> NM at 34
Mn=phi'*M*phi;

Combine Modes Then Develop Story Responses

 disp1=zeros(size(u,1),n);     % Story displacement
 acc=zeros(size(ac,1),n);      % story acceleration
 EF=zeros(size(u,1),n);        % Story force
for st=1:n                     % st is story
    for j=1:nom                % j is mode
        disp2=zeros(size(u,1),n);
        accel=zeros(size(ac,1),n);
        EF1=zeros(size(u,1),n);
        disp2(:,st)=phi(st,j)*u(:,j);                      % Chopra, pp 438; Disp=phi*U
        accel(:,st)=phi(st,j)*ac(:,j);
        EF1(:,st)=omega2(j,j)*M(st,st)*phi(st,j)*u(:,j);   % Chopra, pp 438
        disp1=disp2+disp1;                                 % Displacement of Stories
        acc=accel+acc;                                     % Acceleration of Stories
        EF=EF1+EF;                                         % Seismic Force of Stories
    end
end