Simulating Dynamics - 5.0

I got a bit different answer than what was expected (about $5.288$ ). Let $(\alpha, \beta) = (\frac{1}{a^2}, \frac{1}{b^2})$ . Derive an acceleration constraint equation.

$\alpha x^2 + \beta y^2 = 1 \\ 2 \alpha x \dot{x} + 2 \beta y \dot{y} = 0 \\ 2 \alpha x \ddot{x} + 2 \alpha \dot{x}^2 + 2 \beta y \ddot{y} + 2 \beta \dot{y}^2 = 0$

For convenience, I started the mass at $x = +0.1$ and let it go until $+ 3.5$ . Let $\vec{u}_1$ be the unit normal vector, and let $\vec{u}_2$ be a unit tangent vector opposite to the particle motion. For the sake of brevity, I will not derive these here. Let $F_{gy}$ and $F_{sy}$ be the vertical components of the gravity and spring forces. Write the Newton's 2nd law equations.

$m \ddot{x} = N u_{1x} + \mu N u_{2x} \\ m \ddot{y} = N u_{1y} + \mu N u_{2y} + F_{gy} + F_{sy}$

Plugging these into the constraint equation and solving for the normal force yields:

$N = \frac{-2 \beta y (F_{gy} + F_{sy}) - 2 \alpha m \dot{x}^2 - 2 \beta m \dot{y}^2}{2 \alpha x u_{1x} + 2 \alpha x \mu u_{2x} + 2 \beta y u_{1y} + 2 \beta y \mu u_{2y}}$

Having solved for $N$ , plug these back into the 2nd law equations to solve for the accelerations. Numerical integration takes care of the rest.

Solution code is attached.

  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

import math

dt = 10.0**(-6.0)

m = 30.0
g = 10.0
mu = 0.05
a = 4.0
b = 8.0
k = 5.0
L0 = 2.0

alpha = 1.0/(a**2.0)
beta = 1.0/(b**2.0)

gamma = 1.0/math.sqrt(beta)

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

t = 0.0

x = 0.1
y = gamma*math.sqrt(1.0 - alpha*(x**2.0))

xd = 0.0
yd = 0.0

xdd = 0.0
ydd = 0.0

maxres = 0.0

while x <= 3.5:

    x = x + xd*dt
    y = y + yd*dt

    xd = xd + xdd*dt
    yd = yd + ydd*dt

    u2x = -1.0
    u2y = gamma*alpha*x/math.sqrt(1.0-alpha*(x**2.0))

    u2 = math.hypot(u2x,u2y)

    u2x = u2x/u2
    u2y = u2y/u2

    u1x = u2y
    u1y = -u2x

    Fgy = -m*g

    s = 2.0*y - L0

    Fsy = -k*s

    right = -2.0*beta*y*(Fgy + Fsy) - 2.0*alpha*m*(xd**2.0) - 2.0*beta*m*(yd**2.0)
    left = 2.0*alpha*x*u1x + 2.0*alpha*x*mu*u2x + 2.0*beta*y*u1y + 2.0*beta*y*mu*u2y

    N = right/left

    xdd = (N*u1x + mu*N*u2x)/m
    ydd = (N*u1y + mu*N*u2y + Fgy + Fsy)/m

    res = math.fabs(alpha*(x**2.0) + beta*(y**2.0) - 1.0)

    if res > maxres:
        maxres = res

    t = t + dt

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

print dt
print t
print ""
print x
print y
print ""
print maxres

#>>> 
#1e-05
#5.28854999997

#3.50000446073
#3.87311446683

#1.78214652926e-05
#>>> ================================ RESTART ================================
#>>> 
#1e-06
#5.28849700039

#3.50000029177
#3.87299701653

#1.78218324365e-06
#>>> 

Simulating the dynamics by running a force computation routine is a nice way to attempt this problem.

In this kind of problem, there are two relevant components to each of these forces; tangential and normal. For that, I have calculated the unit tangent and normal vector of an ellipse:

$\displaystyle \vec{\bold{t}} = \frac{-ay}{b} \bold{\hat{i}} + \frac{bx}{a} \bold{\hat{j}}$

$\displaystyle \bold{\hat{t}} = \frac{\vec{\bold{t}}}{|\vec{\bold{t}}|}$

$\displaystyle \vec{\bold{n}} = \frac{bx}{a} \bold{\hat{i}} + \frac{ay}{b} \bold{\hat{j}}$

$\displaystyle \bold{\hat{n}} = \frac{\vec{\bold{n}}}{|\vec{\bold{n}}|}$

Normal forces

Normal reaction to gravity:

$\vec{\bold{F}}_{g \bold{\hat{n}}} = mg \cos(\theta) \bold{\hat{n}}$

Normal reaction to spring force:

$\vec{\bold{F}}_{s \bold{\hat{n}}} = k(2y-l) \cos(\theta) \bold{\hat{n}}$

Total normal force:

$\vec{\bold{N}} = \vec{\bold{F}}_{g \bold{\hat{n}}} + \vec{\bold{F}}_{s \bold{\hat{n}}}$

Tangential forces

Ramp angle:

$\displaystyle \theta = \tan^{-1} \left( \frac{dy}{dx} \right)$

Tangential gravitational force:

$\displaystyle \vec{\bold{F}}_{g \bold{\hat{t}}} = mg \sin(\theta) \bold{\hat{t}}$

Tangential spring force:

$\displaystyle \vec{\bold{F}}_{s\bold{\hat{t}}} = k(2y-l) \sin(\theta) \bold{\hat{t}}$

Friction force (occurs tangentially):

$\displaystyle \vec{\bold{F}}_{\mu\bold{\hat{t}}} = \pm \mu |\vec{\bold{N}}| \bold{\hat{t}}$

Total tangential force

$\displaystyle \vec{\bold{F}}_{t} = \vec{\bold{F}}_{g \bold{\hat{t}}} + \vec{\bold{F}}_{s\bold{\hat{t}}} \pm \vec{\bold{F}}_{\mu\bold{\hat{t}}}$

Simulation 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

time = 0;
deltaT=  10^-7;

x = -0.1;
xDot = 0;

y = sqrt(64*(1-x^2/16));
yDot = 0;


u = 0.05;
m = 30;
g = 10;
a = 4;
b = 8;
k = 5;
l = 2;

xValues = [];
times = [];

while x >= -3.5
    %xValues = [xValues; x];
    %times = [times; time];

    y = sqrt(64*(1 - (x^2/16)));
    theta = atan(-x/(2*sqrt(1-x^2/16))); 


    %unit tangent vector and unit normal to the ellipse 
    tangentVector = [-a*y/b b*x/a];
    unitTangent = tangentVector/norm(tangentVector); %unit tangent vector

    normalVector = [b*x/a a*y/b];
    unitNormal = normalVector/norm(normalVector); %unit normal vector

    %Forces
    Fg = m*g; %gravitational force magnitude
    springForce = k*(2*y - l); %spring force magnitude

    Ftgravity = Fg*sin(theta)*unitTangent;%tangential component (gravity)
    Ftspring = springForce*sin(theta)*unitTangent; %tangential component(spring)

    Fngravity = Fg*cos(theta)*unitNormal; %normal component (gravity)
    Fnspring = springForce*cos(theta)*unitNormal; %normal component (spring)

    %total normal force
    totalNormalForce = Fngravity + Fnspring;


    %friction force
    if xDot < 0 
        frictionForce = -1*u*norm(totalNormalForce)*unitTangent;    
    else
        frictionForce = u*norm(totalNormalForce)*unitTangent;
    end

    %total tangential force
    totalTangentialForce = Ftgravity + Ftspring + frictionForce;

    %numerical integration and tangential acceleration
    acceleration = totalTangentialForce/m;

    xDotDot = acceleration(1);
    xDot = xDot + xDotDot*deltaT;
    x = x + xDot*deltaT;

    time = time + deltaT;

end

%plot(times,xValues);

disp(time);