Dynamic analysis of moving mass by Newmark beta method

17 views (last 30 days)
Hi,
I am new to matlab and have couple of doubts regarding my coding for the dynamic analysis of a beam of moving mass.
I am modelling a beam of 11 nodes and moving mass with a velocity of 1m/s. i wouldlike to plot the displacement of node5 when the load is in each nodes . for example displacement of node5 when load is in node 1,displacement of node5 load is in node 2,displacement of node5 when load is in node3..........displacement of node5 when load is in ith node.
For each time step, i should get the output of displacement of one node. i have attached the code which i have done,i am getting displacement as incremental for each nodes and the displacement is not changing for the before nodes.
could anyone help me in this issue.
Thanks
  3 Comments
Shabnam Sulthana Mohamed Isaque
% Beam properties
L=10; %length in m length 10m and 11 nodes
A=0.00767; %area in m^2
E=2.1e17; %youngs modulus in N/m^2
ne1=10; %no of elements
I=0.0000305500; %moment of inertia in m
nnp=ne1+1; %no of node points
Le1=L\ne1; %length of each element
M = 60; % 60kg/m
P=10000; %Force in 10kN
p=7.860e12; % in kg/mm^3. beam density (per unit volume):kg/m^3
ks=45000000; %stiffness of the spring in 450kN/m
gam=1/2;
beta=1/4; % average acceleration method
C=0; %damping
v = 1; %initial velocity = 1m/s
% 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
numberElements=size(elementNodes,1);
% numberNodes: number of nodes
numberNodes=11;
displacement = zeros(numberNodes,1);
force = zeros(1,numberNodes);
K = zeros(numberNodes);
m=zeros(numberNodes); % mass matrix
% elementNodes: connections at elements
ii=1:numberElements;
elementNodes(:,1)=ii;
elementNodes(:,2)=ii+1;
% Time step
t = linspace(0,10,10);
dt = t(2)-t(1);
n = length(m); %no of nodes
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) ;
% Initial Conditions
depl(1) = 0 ;
vel(1) = 1 ;
P = zeros(n,1) ;
P = 10000;
accl(1) = 0;
%accl(:,1) = (P(1) - C*vel(:,1) - K*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
for i = 1:(length(t)-1)
DPbar = P + A*vel(i)+B*accl(i);
Dv = DPbar/kbar;
depl(i+1) = depl(i) + Dv
veldv = gam*Dv/(beta*dt) - gam*vel(i)/beta + dt*accl(i)*(1-0.5*gam/beta);
accldv= Dv/(beta*dt*dt) - vel(i)/(beta*dt) - accl(i)/(2*beta);
vel(i+1) = vel(i) + veldv;
accl(i+1) = accl(i) + accldv;
end
Jim Riggs
Jim Riggs on 5 Sep 2019
Edited: Jim Riggs on 5 Sep 2019
A few things need clarification:
1) Number of elements is defined 2 different times;
ne1 = 10
numberElements = size(elementNodes,1)
2) elementNodes is defined more than once:
elementNodes = [...] % hard coded
elementNodes(:,1) = ii % computed
elementNodes(:,2)=ii+1
3) P is defined at the top as a scalar, then as an array, then re-defined as a scalar:
P = 10000
....
P = zeros(n,1)
P = 10000
Is P supposed to be a vector or a scalar?
4) You define a set of six constants (a1, a2,...a6) but then don't use them where you could to simplify the equations; e.g.
kbar = ks + C*a1 + M*a2;
A = M*a3 + C*a4;
B = M*a5 + C^2*a6;
...and inside the for loop:
veldv = Dv*a1 - vel(i)*a4 - accl(i)*a6;
accldv = DC*a2 - vel(i)*a3 - accl(i)*a5;
5) The calculation of A overwrites the initial definition of A
A = 0.00767; % area in m^2
...
A = M*a3 + C*a4;
So, appart from this confusion in the code, I still do not understand the question. Perhaps a drawing would help describe the relationship of the moving mass and the beam.

Sign in to comment.

Answers (1)

Samsun Anadolu
Samsun Anadolu on 4 Dec 2020
Good Job, thanks for sharing.
I am a beginner numerical analysis problem.
I have checked similar codes; I couldn't find any clue about the boundary conditions of the system.
Do these systems have a boundary condition? Could you please explain it?

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!