MATLAB Answers

Dynamic analysis of moving mass by Newmark beta method

32 views (last 30 days)
Shabnam Sulthana Mohamed Isaque
Answered: Samsun Anadolu on 4 Dec 2020
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.


Jim Riggs
Jim Riggs on 4 Sep 2019
It will be much easier to read your code if you paste it into your question, then select all of the code and press alt-enter. This formats it as matlab code.
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
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
% numberNodes: number of nodes
displacement = zeros(numberNodes,1);
force = zeros(1,numberNodes);
K = zeros(numberNodes);
m=zeros(numberNodes); % mass matrix
% elementNodes: connections at elements
% 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;
Jim Riggs
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
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?


Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by