Simulating Dynamics - 3

Fun problem. Simulation code is below. Aside from computing the basic dynamics, the key thing about this problem is to keep track of two different forces:

1) The combined force from gravity and the spring parallel to the ramp
2) The maximum available friction force (equal to the friction coefficient multiplied by the normal force)

Every time the block's velocity becomes zero, you have to check force (1) to see if it is bigger/smaller than force (2). If it is bigger, the motion continues. If it is smaller, the block stays where it is in perpetuity. As it turns out, the block comes to rest at a displacement of $\approx 2.438$ at time $t \approx 4.6$ . It stays there from then on.

  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
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110

import math
import numpy as np

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

# Constants

xh = 7.0   # hinge coordinates
yh = 9.0

x0 = 8.0   # coordinates of start point
y0 = 8.0

theta = math.pi/4.0  # ramp angle

uk = 0.5   # friction coefficient
g = 10.0  # gravity

tf = 5.0  # final time

m = 8.0   # mass

k = 50.0  # spring constant
L0 = math.sqrt(2.0)  # spring natural length

dt = 10.0**(-5.0)  # time step

u1x = -1.0/math.sqrt(2.0)  # unit vector down ramp
u1y = -1.0/math.sqrt(2.0)

u2x = -1.0/math.sqrt(2.0)  # unit vector normal to ramp
u2y = 1.0/math.sqrt(2.0)

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

# Initialize simulation

t = 0.0
count = 0

D = 0.0  # distance from start point along ramp
Dd = 0.0
Ddd = 0.0

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

while t <= tf:

    D = D + Dd*dt     # numerical integration
    Dd = Dd + Ddd*dt

    x = x0 + D*u1x   # xy coordinates
    y = y0 + D*u1y

    Sx = xh - x     # vector from block to hinge
    Sy = yh - y

    S = math.hypot(Sx,Sy)   # spring length

    sx = Sx/S   # unit vector from block to hinge
    sy = Sy/S

    Fs = k*(S-L0) # spring force magnitude

    Fsx = Fs * sx  # spring force vector
    Fsy = Fs * sy

    # Resolve spring force into components down ramp and perpendicular to ramp

    # Fsx = F1*u1x + F2*u2x   
    # Fsy = F1*u1y + F2*u2y

    M11 = u1x
    M12 = u2x

    M21 = u1y
    M22 = u2y

    M =  np.array([[M11,M12],[M21,M22]])
    vec = np.array([Fsx,Fsy])

    Sol = np.linalg.solve(M, vec)

    F1 = Sol[0]
    F2 = Sol[1]

    N = m*g*math.cos(theta) - F2   # Normal force

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

    Fgs = m*g*math.sin(theta) + F1    # Force parallel to ramp due to gravity and spring
    Facc = Fgs + Ff                   # Friction force

    Ddd = Facc/m  # acceleration of block along ramp

    t = t + dt
    count = count + 1

    if count % 1000 == 0:
        print t,D,Dd,(Fgs/Ff)

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

print ""
print ""  
print dt
print D

My approach is completely numerical. The main approach with this kind of problem is to compute the net force and thereby determine the acceleration of the block. The forces acting on the block are:

Gravitational force:

This force makes the block move down the ramp, as well as creates a normal force between the ramp and the block. Let's assume a coordinate system where $y$ is normal to the ramp and $x$ is orthogonal to $y$ (in the slope of the ramp).

$\displaystyle \bold{F}_{gx} = mg \sin (\theta)$

$\displaystyle \bold{F}_{gy} = -mg \cos(\theta)$

Spring force:

The spring force components are the most difficult to work out with the new set of bases directions (stated above).

Let the angle the spring makes with the ramp be $\phi$ .

$\displaystyle \sin(\phi) = \frac{\sqrt{2}}{d^2+2}$ (where $d$ is the displacement)

$\displaystyle \cos(\phi) = \frac{d}{\sqrt{d^2 + 2}}$

The magnitude of spring force:

$\displaystyle ks = k(\sqrt{d^2+2} - \sqrt{2})$

Components:

$\displaystyle \bold{F}_{s x} = -k(\sqrt{d^2+2} - \sqrt{2}) \cos(\phi)$

$\displaystyle \bold{F}_{s y} = k(\sqrt{d^2+2} - \sqrt{2}) \sin(\phi)$

Frictional force:

The normal force is the net $y$ force, given by:

$\bold{N} = \bold{F}_{gy} - \bold{F}_{s y}$

So the frictional force is $\pm \mu \bold{N}$ (depending on velocity of block).

Total force in $x$ -direction:

$\bold{F}_x = \bold{F}_{gx} + \bold{F}_{s x} \pm \mu \bold{N}$

Now, we can code the forces in and simulate the motion of the block. Here's a graph of the displacement of the block:

Here's my code:

 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

time = 0;
deltaT = 10^-4;

velocity = 0;
displacement = 0.1;

%constants

L = sqrt(2);
k = 50;
m = 8;
g = 10;
mu = 0.5;

%ramp angle
sin_theta = 1/sqrt(2);
cos_theta = 1/sqrt(2);

displacementValue = [displacement];
timeValue = [time];
velocityValue = [velocity];

while time <= 5

  displacementValue = [displacementValue, displacement];
  timeValue = [timeValue, time];
  velocityValue = [velocityValue, velocity];



  %component multipliers for spring force:    
  sin_phi = L/(hypot(displacement,L));
  cos_phi = displacement/hypot(displacement,L);

  %spring force magnitude:
  stretch = hypot(L,displacement)-L;
  springForce = k*stretch;

  %Gravitational force components:
  Fgx = m*g*sin_theta;
  Fgy = -m*g*cos_theta;

  %Spring force components:
  Fsx = -springForce*cos_phi;
  Fsy = springForce*sin_phi;

  %normal force computation:
  normalForce = Fsy + Fgy;


  %spring force (dependent on direction of velocity):
  if velocity < 0
    frictionForce = -mu*normalForce;
  else
    frictionForce = mu*normalForce;
  end

  %total force:
  totalForceX = Fgx + Fsx + frictionForce;

  %Euler numerical integration:
  acceleration = totalForceX/m;
  velocity = velocity + acceleration*deltaT;
  displacement = displacement + velocity*deltaT;

  time = time + deltaT;

end


plot(timeValue,displacementValue)

hold on

plot(timeValue, velocityValue)

hold off