Simulating Dynamics - 1

My solution is completely numerical, like Steven Chase's. One cannot solve this problem without implementing a numerical solution. The forces involved are too complicated and ever-changing, so it's very hard to do analytically.

My method was to compute the acceleration at each point in time by running a series of force computations. For those who don't understand the theory behind the computer simulation, the forces of the system include gravitational force, spring force, and the frictional force.

The spring force is given by the stretch of the spring and the spring constant:

$\displaystyle |\vec{\bold{F}}_{spring}| = k( \sqrt{1+x^2}-1)$

The horizontal and vertical components can be given by multiplying by $\displaystyle \cos(\theta) = \frac{x}{\sqrt{1+x^2}}$ , $\displaystyle \sin(\theta) = \frac{1}{\sqrt{1+x^2}}$

The opposing force to the upward spring force is gravitational force, given by $mg$ .

Now, we need to consider the normal force to the floor, which is the resultant force in the $y$ -axis, given by:

$\displaystyle \bold{F}_y = \bold{F}_{spring} - \bold{F}_{g} = \frac{k( \sqrt{1+x^2}-1) }{\sqrt{1+x^2}} - mg$ (normal force)

Finally, we can formulate an equation that gives us the force in the $x-axis$ :

$\displaystyle \bold{F}_x = \frac{|\vec{\bold{F}}_{spring}|x}{\sqrt{1+x^2}} - \mu \bold{F}_y$

With this clear in the mind, we can program these forces and compute the total force along the $x$ -axis, thereby determining the acceleration. We can then numerically integrate the acceleration to find the distance travelled given a small time-step. That's the logic behind the computer program.

Graphs produced with matlab:

Here's the python code (using Explicit Euler instead due to slow speed with RK 4):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

time = 0;
deltaT = 10^-4; %timestep


x = 10; %initial displacement
xDot = 0; %initial velocity


xValues = [];
timeValues = [];

%constants

k = 20;
m = 5;
g = 10;
mu = 0.7;
initialTotalEnergy = (k*(hypot(1,x)-1)^2)/2;

kineticEnergies = []; 
potentialEnergies = [];
dissipatedEnergies = [];

while time <= 5 %evolve till time = 5
  xValues = [xValues; x];
  timeValues = [timeValues; time];

  %Energy variables:  

  Kinetic = (m*xDot^2)/2;
  Potential = (k*(hypot(1,x)-1)^2)/2;
  totalEnergy = Kinetic + Potential;
  dissipatedEnergy = initialTotalEnergy - totalEnergy;

  kineticEnergies = [kineticEnergies; Kinetic];
  potentialEnergies = [potentialEnergies; Potential];
  dissipatedEnergies = [dissipatedEnergies; dissipatedEnergy];



  %spring force magnitude
  springForce = k*(hypot(1,x)-1);

  %spring force components
  Fsx = -springForce*x/(hypot(1,x));
  Fsy = springForce/(hypot(1,x));


  %gravitational force 
  gravityForceY = -m*g;


  %total forces in y = spring + gravity
  totalForceY = Fsy + gravityForceY;



  %frictional force = mu*N
  if xDot < 0 %negative or positive dependent on direction of velocity
    frictionForceX = -mu*totalForceY;
  else
    frictionForceX = mu*totalForceY;
  end


  %total force in x = spring force x + friction
  totalForceX = Fsx + frictionForceX;


  %Explicit Euler numerical integration
  xDotDot = totalForceX/m;
  xDot = xDot + xDotDot*deltaT;
  x = x + xDot*deltaT;

  %update time value
  time = time + deltaT;
end

plot(timeValues, xValues  );

Fun problem. This sort of problem benefits greatly from a code solution

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

import math

# Constants

m = 5.0
k = 20.0
L0 = 1.0
h = 1.0
g = 10.0
u = 0.7

dt = 10.0**(-6.0)

##############################

t = 0.0

x = 10.0            # Initialize simulation
xd = 0.0
xdd = 0.0

while t <= 5.0:

    x = x + xd*dt     # numerical integration
    xd = xd + xdd*dt

    L = math.hypot(x,h)  # spring length

    s = L - L0  # spring stretch

    Fs = k*s  # spring force magnitude

    Fsx = -Fs*(x/L)     # horizontal spring force
    Fsy = Fs*(h/L)      # vertical spring force

    # m*g = Fsy + N    # vertical force balance

    N = m*g - Fsy     # normal force

    if xd > 0.0:      # friction force
        Ff = -u*N     # depends on velocity sign   
    else:
        Ff = u*N

    Fx = Fsx + Ff     # net x force

    xdd = Fx/m        # x acceleration

    t = t + dt

##############################

print dt
print x

#>>> ================================ RESTART ================================
#>>> 
#0.0001
#0.668524353854
#>>> ================================ RESTART ================================
#>>> 
#1e-05
#0.673516137027
#>>> ================================ RESTART ================================
#>>> 
#1e-06
#0.67401590979
#>>>