# how to do Newmark method for 2 different timesteps and to calculate mipoint deflection with respect to loads in all nodes?

2 views (last 30 days)
Hi everyone,
i'm new to matlab.For my research, i am calculating dynamic analysis of beam for 2 different timesteps where i coulnt get the correct results. i have to calulate midpoint displacement of the beam when the load at each point (drawing attached). could anyone help me to solve this.
% 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
ne1=10; %no of elements
I=0.000636056; %moment of inertia in m^4
nnp=ne1+1; %no of node points
Le1=L\ne1; %length of each element
M = 588.42; % 60kg/m
P=10000; %Force in 10kN
ks = 4.5e+11; %stiffness
gam=1/2;
beta=1/4; % average acceleration method
C=0; %damping
% elementNodes: connections at elements
elementNodes=[1 2;2 3;3 4;4 5;5 6;6 7;7 8;8 9;9 10;10 11];
% numberElements: number of Elements
numberEements=size(elementNodes,1);
% numberNodes: number of nodes
numberNodes=11;
% elementNodes: connections at elements
ii=1:numberElements;
elementNodes(:,1)=ii;
elementNodes(:,2)=ii+1;
% Time step
t = linspace(0,11,11); % 11 timestep and 110 timesteps (total simulation time 11sec)
dt = t(2)-t(1);
nt = length(t);
% Constants used in Newmark's integration
a1 = gam/(beta*dt) ; a2 = 1/(beta*dt^2) ;
a3 = 1/(beta*dt) ; a4 = gam/beta ;
a5 = 1/(2*beta) ; a6 = (gam/(2*beta)-1)*dt ;
depl = zeros(nt,1) ;
vel = zeros(nt,1) ;
accl = zeros(nt,1) ;
P = zeros(nt,1) ;
% Initial Conditions
depl(1) = 0 ;
depl(11) = 0;
vel(1) = 0 ;
P(1) = 10000;
accl(1) = (P(1) - C*vel(1) - ks*depl(1))/ M;
kbar = ks + gam*C/(beta*dt) + M /(beta*dt*dt);
A = M /(beta*dt) + gam*C/beta;
B = M /(2*beta) +dt*C*((0.5*gam/beta)-1)*C;
% Time step starts
DPbar = zeros(nt,1);
Dv = zeros (nt,1);
veldv = zeros(nt,1);
accldv = zeros(nt,1);
%calculation of timestep
for i = 1:(length(t)-1)
DPbar(i) = P(1) + A*vel(i)+B*accl(i);
Dv(i) = DPbar(i)/kbar;
depl(i+1) = depl(i) + Dv(i);
veldv(i) = gam*Dv(i)/(beta*dt) - gam*vel(i)/beta + dt*accl(i)*(1-0.5*gam/beta);
accldv(i)= Dv(i)/(beta*dt*dt) - vel(i)/(beta*dt) - accl(i)/(2*beta);
vel(i+1) = vel(i) + veldv(i);
accl(i+1) = accl(i) + accldv(i);
end 