# NDOF Linear Shear Building - Forcing Function problem set up

5 views (last 30 days)
Tariq Azizi on 2 Mar 2017
I have developed this example where a simple NDOF shear building is analyzed under Step Load [P], for damped/undamped condition using newmark's method. However, I am having difficulties setting up the Forcing Matrix as F(x)=F1+F2*Sin(Omega*t); Right now it is just a step load matrix, with pre-defined forces in each floor.
if true
% code
end
% Code by : Azizi |
% CEE 536 - Structural Dynamics |
% |
%--------------------------------------------------------------------------
%
% This MATLAB program solves the NDOF linear shear building model using
% newmark's numerical integration method. In order to use the program;
% first describe the initial conditions and bounder's as matrixes. Values
% to input are ndof,M,K,C,P,and type of newmark acceleration.
%
%--------------------------------------------------------------------------
clear ;clc ;
EI = 5*10^9 ;
h = 100 ;
K = (EI/h^3)*[6 -4 0 0 0 ; % K - Stiffness Matrix (INPUT)
-4 10 -6 0 0 ;
0 -6 11 -4 0 ;
0 0 -4 6 -2 ;
0 0 0 -2 2];
ndof = length(K) ;
m = 200 ;
M = m*diag([2 1 1 1 1]) ; % M - Mass Matrix (INPUT)
%%%%%%%Modal Analysis
[V,D]=eig(K,M);
[W,k]=sort(diag(D));
V=V(:,k);
Factor=diag(V'*M*V);
Phi=V*inv(sqrt(diag(Factor)));
Omega=diag(sqrt(Phi'*K*Phi));
w = diag(Omega(1:ndof)) ;
%-------------------------------------------------------%
%%%%%%%%%%%%%Provide damping option &&&&&&&
sys = 'damped' ; %Select if system Damped
switch sys
case 'undamped' % C - If System undamped
C = zeros(size(K)) ;
case 'damped'
C= 2*w*0.2; % C - Damping Matrix (Damped)
end
%-------------------------------------------------------
%%%%%%%%%%%%%NEWMARK'S METHOD : LINEAR SYSTEM &&&&&&&
acc = 'Linear' ; % acc - Type of Newmark's Method to be used
switch acc
case 'Average'
gama = 1/2 ;beta = 1/4 ;
case 'Linear'
gama = 1/2 ;beta = 1/6 ;
end
% Time step
ti = 0. ;
tf = 20. ;
dt = 0.1 ;
t = ti:dt:tf ;
nt = fix((tf-ti)/dt) ;
n = length(M) ;
% Constants used in Newmark's integration
a1 = gama/(beta*dt) ; a2 = 1/(beta*dt^2) ;
a3 = 1/(beta*dt) ; a4 = gama/beta ;
a5 = 1/(2*beta) ; a6 = (gama/(2*beta)-1)*dt ;
P = 100*[10 20 30 40 50]' ; % P - Force MAtrix
disp = zeros(n,nt) ;
vel = zeros(n,nt) ;
accl = zeros(n,nt) ;
U = zeros(ndof,nt) ;
% Initial Conditions for displacement and acceleration
% Individual input for initial position of each building can be made:
disp(:,1) = zeros ; %Building is assumed at rest
vel(:,1) = zeros ;
U(:,1) = Phi*disp(:,1) ;
P0 = zeros(n,1) ;
accl(:,1) = M\(P-C*vel(:,1)-K*disp(:,1)) ;
Kcap = K+a1*C+a2*M ;
a = a3*M+a4*C ;
b = a5*M+a6*C ;
% Tme step starts
for i = 1:nt
delP = P0+a*vel(:,i)+b*accl(:,i) ;
delu = Kcap\delP ;
delv = a1*delu-a4*vel(:,i)-a6*accl(:,i) ;
dela = a2*delu-a3*vel(:,i)-a5*accl(:,i);
disp(:,i+1) = disp(:,i)+delu ; % disp - modal displacement's
vel(:,i+1) = vel(:,i)+delv ; % vel - modal velocities
accl(:,i+1) = accl(:,i)+dela ; % accl - modal accelerations
U(:,i+1) = Phi*disp(:,i+1) ; % U - system's displacement
end
%%%%%%%%GRAPHICAL PLOT%%%%%%
figure(1);
subplot(2,3,1);
plot(t,disp(1,:),'-','Color','blue') ;
xlabel('Time in sec') ; ylabel('displacement');
title('First Floor Displacement');
subplot(2,3,2);
plot(t,disp(2,:),'-','Color','red');
xlabel('Time in sec') ; ylabel('displacement');
title('Second Floor Displacement');
subplot(2,3,3)
plot(t,disp(3,:),'-','Color','yellow');
xlabel('Time in sec') ; ylabel('displacement');
title('Third Floor Displacement');
subplot(2,3,4)
plot(t,disp(4,:),'-','Color','green');
xlabel('Time in sec') ; ylabel('displacement');
title('Fourth Floor Displacement');
subplot(2,3,5)
plot(t,U(5,:),'-','Color','black');
xlabel('Time in sec') ; ylabel('displacement');
title('Building Displacement');