Simulating Dynamics - 2

Just a different flavor of the solution posted by @Karan Chatrath . I have a $2 \times 2$ linear system to solve for $T$ and $\ddot{x}$ , so my expressions are correspondingly more complex than in the $4 \times 4$ linear system. Let the pulley height be $h$

The string length constraint gives (after a great deal of chain-ruling):

$\frac{m_2 g - T}{m_2} = x^2 \dot{x}^2 (x^2 + h^2)^{-3/2} - (x \ddot{x} + \dot{x}^2)(x^2 + h^2)^{-1/2}$

And the dynamics of the block with friction give the other equation ( $\theta$ is the angle between the string and the horizontal):

$T \cos \theta - \mu(m_1 g - T \sin \theta) = -m_1 \ddot{x}$

Solve the system for $T$ and $\ddot{x}$ on each time step, and use numerical integration. The block reaches $x = 0$ at approximately $t = 3.87$

 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

import math
import numpy as np

h = 8.0
x0 = math.sqrt(14.0**2.0 - h**2.0)

m1 = 3.0
m2 = 2.0
u = 0.5
g = 10.0

dt = 10.0**(-4.0)

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

t = 0.0

x = x0
xd = 0.0
xdd = 0.0

while x >= 0.0:

    x = x + xd*dt
    xd = xd + xdd*dt

    cos = x/math.sqrt(x**2.0 + h**2.0)
    sin = h/math.sqrt(x**2.0 + h**2.0)

    q1 = (x**2.0 + h**2.0)**(-1.0/2.0)
    q2 = (x**2.0 + h**2.0)**(-3.0/2.0)

    J11 = 1.0/m2
    J12 = -x*q1

    J21 = cos + u*sin
    J22 = m1

    J = np.array([[J11,J12],[J21,J22]])

    vec1 = g - (x**2.0)*(xd**2.0)*q2 + (xd**2.0)*q1
    vec2 = u*m1*g

    vec = np.array([vec1,vec2])

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

    T = Sol[0]
    xdd = Sol[1]

    t = t + dt

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

print dt
print t

#>>> 
#0.001
#3.868
#>>> ================================ RESTART ================================
#>>> 
#0.0001
#3.8658
#>>> 

The equations of motion have been re-arranged in a matrix form in the solution. I have not necessarily used the same sign convention as you might have.

$m_2g-T = m_2\ddot{y}$ $T\cos{\theta} - \mu N = m_1 \ddot{x}$ $T\sin{\theta} + N = m_1 g$

Finally, the constraint equation is:

$x^2+64=(14-y)^2$

Double differentiation gives the fourth equation using which $T$ , $N$ , $\ddot{x}$ and $\ddot{y}$ can be solved for. I did that by re-arranging the equations in a matrix form.

$\theta = \arccos\left(\frac{-x}{14-y}\right)$

$\left[\begin{matrix} m_1&0&\mu& -\cos{\theta} \\ 0&m_2&0&1 \\ 0&0&1&\sin{\theta}\\2x&28-2y&0&0\end{matrix}\right]\left[\begin{matrix} \ddot{x}\\ \ddot{y} \\ N \\ T\end{matrix}\right]=\left[\begin{matrix} 0\\m_2g \\ m_1g \\ 2\dot{x}^2-2\dot{y}^2\end{matrix}\right]$

$\implies A \left[\begin{matrix} \ddot{x}\\ \ddot{y} \\ N \\ T\end{matrix}\right]=b$ $\implies \left[\begin{matrix} \ddot{x}\\ \ddot{y} \\ N \\ T\end{matrix}\right] = A^{-1}b$

Uncommented simulation code attached below. Used the modified Euler step for numerical integration:

 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

clear all
clc

m1 = 3;
m2 = 2;

g  = 10;
mu = 0;

x(1) = -14*cos(asin(8/14));
y(1) = 0;

dx(1) = 0;
dy(1) = 0;

dt    = 1e-5;
tf    = 3.8655;

t     = 0:dt:tf;

for k = 1:length(t)-1


    A = [m1 0 mu -x(k)/(y(k)-14);0 m2 0 1;0 0 1 (1 - x(k)^2/(y(k) - 14)^2)^(1/2);2*x(k) 28-2*y(k) 0 0];
    b = [0;m2*g;m1*g;2*dy(k)^2-2*dx(k)^2];

    S = inv(A)*b;

    ddx = S(1);
    ddy = S(2);

    dx(k+1) = dx(k) + ddx*dt;
    dy(k+1) = dy(k) + ddy*dt;

    x(k+1)  = x(k) + dx(k+1)*dt;
    y(k+1)  = y(k) + dy(k+1)*dt;
end