# how to perform Newmark method for fixed beam?

3 views (last 30 days)
hi,
i am new to matlab. i am doing dynamic analysis of a beam using newmark method with different time steps at the midpoint point of node 6 of the beam. the outputs are not correct as i expected. Displacement is linear and not parabloic. could anyone help me out where did i do the mistake.
% Beam properties
L=10; %length in m length 10m
A=0.00767; %area in m^2
E=2.1e12; %youngs modulus in N/m^2
n=10; %no of elements
I=0.000636056; %moment of inertia in m^4
nn=n+1; %no of node points
Le=L\n; %length of each element
m = 588.42; % 60kg/m
P=10000; %Force in 10kN
p = 7800;
%ks = 4.5e+11; %stiffness
Gamma=1/2;
Beta=1/4; % average acceleration method
C=0; %damping
%global stiffness and mass matrix for fixed beam
Ka = E*I/Le^3 *[24 -12 0 0 0 0 0 0 0;
-12 24 -12 0 0 0 0 0 0;
0 -12 12 -12 0 0 0 0 0;
0 0 -12 12 -12 0 0 0 0;
0 0 0 -12 12 -12 0 0 0;
0 0 0 0 -12 12 -12 0 0;
0 0 0 0 0 -12 12 -12 0;
0 0 0 0 0 0 -12 24 -12;
0 0 0 0 0 0 0 -12 24];
Ma = p*A*Le/2 * [588.42 0 0 0 0 0 0 0 0;
0 588.42 0 0 0 0 0 0 0;
0 0 588.42 0 0 0 0 0 0;
0 0 0 588.42 0 0 0 0 0;
0 0 0 0 588.42 0 0 0 0;
0 0 0 0 0 588.42 0 0 0;
0 0 0 0 0 0 588.42 0 0;
0 0 0 0 0 0 0 588.42 0;
0 0 0 0 0 0 0 0 588.42];
%C = Dalpha * Ka + Dbeta * Ma;
% integration constant
c1 = 1/Beta/dt^2;
c2 = Gamma/Beta/dt;
c3 = 1/Beta/dt;
c4 = 1/2/Beta-1;
c5 = Gamma/Beta-1;
c6 = (Gamma/2/Beta-1)*dt;
c7 = (1-Gamma)*dt;
c8 = Gamma*dt;
% % Time step
t = linspace(1,10,11); % 10 timestep
dt = 1;
nt = length(t);
P =zeros(9);
P(:,1) = 10000;
%intial conditions
Fbar = zeros(9);
Dv = zeros(9);
veldv = zeros(9);
accldv = zeros(9);
D = zeros(9);
V = zeros(9);
Aa = zeros(9);
D(:,1) = 0;
V(:,1) = 0;
Aa(:,1) = (Ma^-1)*(P(:,1) - C*V(:,1) - Ka*D(:,1));
Kbar = c1*Ma+c2*C+Ka;
%A = Ma /(Beta*dt) + Gamma*C/beta;
%B = Ma /(2*Beta) +dt*C*((0.5*Gamma/Beta)-1)*C;
for i = 1:(nt-1)
Fbar(:,i+1) = P(:,1) + Ma * (c1 * D(:,i) + c3 * V(:,i) + c4 * Aa(:,i)) + C * (c2 * D(:,i) + c5 * V(:,i) + c6 * Aa(:,i));
Dv(:,i+1) = Kbar^-1 * Fbar(:,i+1);
D(:,i+1) = D(:,i) + Dv(:,i+1);
accldv(:,i+1) = ((c1 * D(:,i+1)) - D(:,i)) - (c3 *V(:,i)) - (c4 * Aa(:,i));
Aa(:,i+1) = Aa(:,i) + accldv(:,i+1);
veldv(:,i+1) = V(:,i) + (c7 * Aa(:,i)) + (c8 * Aa(:,i+1));
V(:,i+1) = V(:,i) + veldv(:,i+1);
end