Dynamic analysis of moving mass by Newmark beta method

Asked by Shabnam Sulthana Mohamed Isaque on 4 Sep 2019

Latest activity Edited by Jim Riggs on 5 Sep 2019

11 views (last 30 days)

tenmetre.m

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

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 on 5 Sep 2019

% 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 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.

You are now following this question

Shabnam Sulthana Mohamed Isaque (view profile)