Monthly Archives: November 2011

by

The fall of the slinky I

The video

The video shows how a spring suspended from one of its ends reacts when its dropped. It can be observed that the lower end “doesn’t know what happened” until a wave propagates to it. In this post we will make a computer simulation of its behavior, to see if we can reproduce the phenomenon, and in the next we will apply a more analytic approach to the same problem.

Discretization

As we are studying the actions of gravity and inertia, we need to model the slinky as a massive spring. In (macroscopic) reality the mass and elasticity are continuously distributed throughout the slinky (macroscopically speaking) but, for the purpose of simulating its behavior with a computer model, we will represent the object as a series of masses connected with massless ideal springs:

Slinky discretization, showing how the discretized element properties relate to the original ones.

If we divide a slinky of total mass M in N small masses, each of them will have the value m = M/N. As the spring constants of series springs add as the resistances of parallel resistors, if we have a slinky with overall spring constant K divided into N-1 smaller springs, each of them will have a bigger spring constant, k = K(N-1) and (obviously) a smaller unstressed length, l_0 = L/(N-1).

Writing the simulation

Now it’s just a question of applying Newton’s laws of motion and the ideal spring equation to get a system of ordinary differential equations:

\displaystyle \frac{d^2x_i}{dt^2} = g + \frac{k}{m}\left[(x_{i+1}-x_i-l_0)H[i]-(x_i-x_{i-1}-l_0)H[N-1-i]\right]

where x_i are the coordinates of the masses (with i going from 0 to N-1), g is the acceleration of gravity, k is the spring constant of each massless spring, m is the value of each small mass, l_0 is the unstressed length of each massless spring and H[n] is the discrete Heaviside step function (used to avoid depending on undefined values).

This second order system can be reduced to system of first order ODEs in the usual way. Integrating it using RK4, we get the following Python code:

def get_xdot(x):
    sk = K * (NX - 1)
    sl = L / (NX - 1)
    mm = M / NX
    xdot = [x[NX + i] if i < NX else G
            for i in range(2 * NX)]
    for i in range(NX - 1):
        a = sk * (x[i + 1] - x[i] - sl) / mm
        xdot[NX + i] += a
        xdot[NX + i + 1] -= a
    return xdot

def rk4_step(x, dt):
    k1 = get_xdot(x)
    k2 = get_xdot([x[i] + dt * k1[i] / 2.0 for i in range(len(x))])
    k3 = get_xdot([x[i] + dt * k2[i] / 2.0 for i in range(len(x))])
    k4 = get_xdot([x[i] + dt * k3[i] for i in range(len(x))])
    return [x[i] + dt * (k1[i] + 2.0 * (k2[i] + k3[i]) + k4[i]) / 6.0
            for i in range(len(x))]

Now we need to define the initial conditions. If we just start with the masses separated by a distance l_0 and at rest, we won’t match the initial conditions of the slinky, because it was being stretched by the action of gravity. It’s not very difficult to compute initial conditions that leave the system at rest if it’s being held:

def initial_x():
    sl0 = L / (NX - 1)
    mm = M / NX
    sk = K * (NX - 1)
    w = M - mm
    x = [0.0 for i in range(2 * NX)]
    for i in range(1, NX):
        x[i] = x[i - 1] + sl0 + w / sk
        w -= mm
    return x

The remaining code is just plumbing and matplotlib presentation code. The whole program can be seen at the repository.

Running the simulation

If we run the simulation with the parameters

NT = 1000 # number of timesteps
NX = 40   # number of masses
T = 1.0   # simulation duration

L = 0.5   # slinky length
K = 1.0   # slinky overall spring constant
M = 1.0   # slinky mass
G = 1.0   # gravitational acceleration

we get the following results:

A simulation where the springs are too soft, giving some negative spring lengths (and consequent overlapping) near t = 1.

In this plot the gray lines represent the trajectory of the small masses and the black lines the trajectory of the slinky’s ends.

Clearly the springs are too soft and we are getting unphysical results, as we spring lengths go negative when t nears 1. To fix that, let’s run the simulations with a greater spring constant, K = 5:

A simulation where the springs have a physically reasonable constant, giving an intriguing behavior to its ends.

Now we get a more reasonable result, showing a phenomenon that is more similar to the one observed in the video: the bottom of the slinky remains in place while the top begins to fall. Now we can check if the slinky remains stationary when held by m_0:

def get_xdot(x):
    sk = K * (NX - 1)
    sl = L / (NX - 1)
    mm = M / NX
    xdot = [x[NX + i] if i < NX else G
            for i in range(2 * NX)]
    for i in range(NX - 1):
        a = sk * (x[i + 1] - x[i] - sl) / mm
        xdot[NX + i] += a
        xdot[NX + i + 1] -= a
    xdot[NX] = 0.0 # holding the slinky
    return xdot

Simulation to check if the slinky remains stationary when held from its upper end.

Conclusions

This simulation validates the main point of the original video: the lower end “doesn’t know” the upper end was released until a compression wave reaches it, at t \approx 0.45 in our simulation. But the detailed behavior differs, as the slinky only shows a compression wave once it reaches the nonlinear regime (when is no more space between the spires).

In the next post we will show an analysis of this nonlinear behavior and the analytical solution to the idealized slinky drop that was numerically simulated in this post.

Advertisement