Problem in creating physics engine for a bouncing spring-mass cube

I am trying to create a physics engine to visualize a bouncing spring-mass cube. The cube consists of 8 masses and 12 linear springs. A fixed value is used for the masses and the spring constants. Although I am using simple Newtoniam mechanics to set all required physics element, I got stuck and don't know how to caculate the net forces on each mass at every time (T+dt). I appreciate your support with it

%Dynamic Simulation

% 1)Time increament: T = T+dt;

% 2) Interaction:

% - Springs: Caculate compression/extension using F=k(l-li)

% - Masses: Tally forces on each mass

% - Add gravity, collision, external forces (if any)

% 3) Integration step

% - Calculate a = F/m

% - Update v = v +dt*a;

% - Update P = P +v*dt

% Repeat

% if mass is fixed, set v=0;

% if dampening, reduce speed v=v*0.999

% if colliding with ground, apply restoring forces

%Collision with the ground:

% - Apply a ground reaction force Fg once a node goes beloe ground (z<0)

% - F=(0,0,-z*Kg); Kg is ground spring coeaffient: 100,000

clc

clear

close all

%initial conditions, all units are in Kg,m, and s

p = zeros(8,3,10000); % 3D vector for position of 8 elements

v = zeros(8,3,10000); % 3D vector for velocity of 8 elements

a = zeros(8,3,10000); 3D vector for acceleration of 8 elements

m = 0.1 %scalar value

F = ones(8,3,10000); %Arbitrary 1N force in 3 dimensions (x,y,z)

Fs = zeros(8,3,10000); %Spring restoring force

Fr = zeros(8,3,10000); %Reaction force during collisions or hitting ground

Ft = zeros(8,3,10000); %Total force

dt = 0.0001; %arbitrary value - 10000 time steps are needed to reach 1 sec

Li = 0.1; %arbitrary value for spring constants

L = zeros(8,3,10000); %Array to store the spring extension length at each time step

DL = L-Li; %calculating difference in spring elongation

k = 10000; %arbitrary spring constants

for T = 1:10000 %loop 10000 times to calculate the forces on each mass

for i = 1:8 %loop for 8 masses

for j = 1:3 % to calculate forces in 3 dimensions

Ft(i,j,T) = F(i,j,T)+Fs(i,j,T)+Fr(i,j,T);

a(i,j,T) = Ft(i,j,T)/m(i,j,T);

v(i,j,T) = v(i,j,T) + a(i,j,T)*dt;

p(i,j,T) = p(i,j,T) + v(i,j,T)*dt;

L(i,j,T) = abs(p(i,j,T));

DL(i,j,T) = L(i,j,T) - Li;

F(i,j,T) = m(i,j,T)*a(i,j,T);

Fs(i,j,T) = k*DL(i,j,T);

Fr(i,j,T) = -F(i,j,T);

end

end

end

[Yongjian Feng]

Sounds like a physics question, not a matlab question.

[Ahmed Ghazal]

mix of both actually

Answer 1

[Divyajyoti Nayak]

Hi @Ahmed Ghazal,

To calculate the net forces on the block at every time step, the forces need to be updated at every time step using the already calculated values in the previous/current time step. In the code provided the values at each time step are being calculated with the initialized values and not with previously calculated values, hence the output values are constant.

%F(i,j,T) is being calculated with only initial value of F(i,j,T)
Ft(i,j,T) = F(i,j,T)+Fs(i,j,T)+Fr(i,j,T);

This should be changed to something like this:

%F(i,j,T+1) is being calculated with previously calculated value of F(i,j,T)
Ft(i,j,T+1) = F(i,j,T)+Fs(i,j,T)+Fr(i,j,T); 
a(i,j,T+1) = Ft(i,j,T+1)/m;  
v(i,j,T+1) = v(i,j,T) + a(i,j,T+1)*dt;
p(i,j,T+1) = p(i,j,T) + v(i,j,T+1)*dt;
L(i,j,T+1) = abs(p(i,j,T+1));
DL(i,j,T+1) = L(i,j,T+1) - Li;
F(i,j,T+1) = m*a(i,j,T+1);
Fs(i,j,T+1) = k*DL(i,j,T+1); 
Fr(i,j,T+1) = -F(i,j,T+1);