twostory1.m

it is about matlam m file useing in matlab

%   Example MATLAB file for dynamic analysis of a frame
%   impacted by a dynamic load.  Newmark's method is used to 
%   carry out direct integration with respect to time.  
clear all
close('all', 'hidden')
%
%   Set up the mass matrix
%
density = 7850.0;
massbeam = 0.1 * 0.2 * 10 * density;
masscolumn = 0.1 * 0.1 * 10 * density;
m1floor = massbeam + 2 * masscolumn;
m2floor = massbeam + masscolumn;

mass = [m1floor 0; 0 m2floor]

%   Set up the stiffness matrix

inercol = 0.1 * 0.1^3 / 12;
inerbeam = 0.2 * 0.1^3 / 12;

emodulus = 2.0e11;

k1story = 2*(12 * emodulus * inercol / 10^3);
k2story = k1story;

stiff = [k1story+k2story -k1story; -k1story k2story]

%    Find the natural frequencies and periods of the system
[eigenvecs eigenvals] = eig(stiff, mass);
freqs = sqrt(eigenvals)

period1 = 2*pi/freqs(1,1)
period2 = 2*pi/freqs(2,2)

%    Use the constant average acceleration method
alpha = 0.5;
beta = 0.25;

%    Specify the time increment to be used for integration
dt = 0.2;
dt = 1
dt = 0.02;
total_simulation_time = 100 ;
%    Define values for Newmark-Beta integration
a0 = 1/(beta*dt^2);
a2 = 1/(beta*dt);
a3 = 1/(2*beta) - 1;
a6 = dt*(1-alpha);
a7 = alpha*dt;

%   Initialize the displacement, velocity, and acceleration
disp = [0; 0];
vel = [0; 0];
acc = [0; 0];
%
keff = stiff + a0*mass;
keff2 = inv(keff);
% 
time = dt:dt:total_simulation_time;
num_time_steps = length(time);

%   Loop through each time step
for i=1:1:num_time_steps

	if time(i)<=1
		pt = [100*time(i); 200*time(i)];
	elseif time(i)>1 & time(i)<=2
		pt = [100; 200];
	elseif time(i)>2 & time(i)<=3
		pt1 = -200*time(i) + 500;
		pt2 = 2 * pt1;
		pt = [pt1; pt2];
    elseif time(i)>3 & time(i)<=5
        pt1 = 50*time(i) - 250;
        pt2 = 2 * pt1;
        pt = [pt1; pt2];
	elseif time(i)>5
		pt = [0; 0];
	end
%
	fefft = pt + mass*(a0*disp + a2*vel + a3*acc);
	dispt = keff2 * fefft;
%
%    Save data for a graphs
	force(i,:) = [time(i), pt(1), pt(2)];

	x(i,:) = [time(i), dispt'];

	acct = a0*(dispt-disp) - a2*vel - a3*acc;

	velt = vel + a6*acc + a7*acct;

	disp = dispt;
	vel = velt;
	acc = acct;
end;
%
%   Plot the force applied to each floor on the same graph.
hold on
plottime = force(:,1);
plotforce1 = force(:,2);
plotforce2 = force(:,3);
plot(plottime, plotforce1, 'co', plottime, plotforce1, 'c-');
plot(plottime, plotforce2, 'r+', plottime, plotforce2, 'r-.');
grid on
xlabel('Time (sec)')
ylabel('Force (N)')
title('Forces Applied to a Two-Story Building')
hold off
%
x1 = x(:,1);
x2 = x(:,2);
x3 = x(:,3);
%
%   Plot the response of both floors on the same graph.
%     Also, start a new figure (new window).
figure
hold on
grid on
plot(x1, x2, 'c+', x1, x2, 'c-')
plot(x1, x3, 'rx', x1, x3, 'r-.')
xlabel('Time (sec)')
ylabel('Displacement (m)')
title('Displacement at Each Floor of a Two-Story Building')