chap24.ipynb

Kernel: Python 3 (ipykernel)

Printed and electronic copies of Modeling and Simulation in Python are available from No Starch Press and Bookshop.org and Amazon.

Rotation

Modeling and Simulation in Python

Copyright 2021 Allen Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

# download modsim.py if necessary

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)
    
download('https://raw.githubusercontent.com/AllenDowney/' +
         'ModSimPy/master/modsim.py')

# import functions from modsim

from modsim import *

This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. Click here to access the notebooks: https://allendowney.github.io/ModSimPy/.

In this chapter and the next we'll model systems that involve rotating objects. In general, rotation is complicated. In three dimensions, objects can rotate around three axes and many objects are easier to spin around some axes than others. If the configuration of an object changes over time, it might become easier or harder to spin, which explains the surprising dynamics of gymnasts, divers, ice skaters, etc. And when you apply a twisting force to a rotating object, the effect is often contrary to intuition. For an example, see this video on gyroscopic precession: http://modsimpy.com/precess.

We will not take on the physics of rotation in all its glory; rather, we will focus on simple scenarios where all rotation and all twisting forces are around a single axis. In that case, we can treat some vector quantities as if they were scalars; that is, simple numbers.

The fundamental ideas in these examples are angular velocity, angular acceleration, torque, and moment of inertia. If you are not already familiar with these concepts, don't worry; I will define them as we go along, and I will point to additional reading.

The Physics of Toilet Paper

As an example of a system with rotation, we'll simulate the manufacture of a roll of toilet paper, as shown in this video: https://youtu.be/Z74OfpUbeac?t=231. Starting with a cardboard tube at the center, we will roll up 47 m of paper, a typical length for a roll of toilet paper in the U.S. (see http://modsimpy.com/paper).

The following figure shows a diagram of the system: $r$ represents the radius of the roll at a point in time. Initially, $r$ is the radius of the cardboard core, $R_{min}$. When the roll is complete, $r$ is $R_{max}$.

I'll use $\theta$ to represent the total rotation of the roll in radians. In the diagram, $d\theta$ represents a small increase in $\theta$, which corresponds to a distance along the circumference of $r~d\theta$.

I'll use $y$ to represent the total length of paper that's been rolled. Initially, $\theta=0$ and $y=0$. For each small increase in $\theta$, there is a corresponding increase in $y$:

$$dy = r~d\theta$$

If we divide both sides by a small increase in time, $dt$, we get a differential equation for $y$ as a function of time.

$$\frac{dy}{dt} = r \frac{d\theta}{dt}$$

As we roll up the paper, $r$ increases. Assuming it increases by a fixed amount per revolution, we can write

$$dr = k~d\theta$$

where $k$ is an unknown constant we'll have to figure out. Again, we can divide both sides by $dt$ to get a differential equation in time:

$$\frac{dr}{dt} = k \frac{d\theta}{dt}$$

Finally, let's assume that $\theta$ increases at a constant rate of $\omega = 300$ rad/s (about 2900 revolutions per minute):

$$\frac{d\theta}{dt} = \omega$$

This rate of change is called an angular velocity. Now we have a system of differential equations we can use to simulate the system.

Setting Parameters

Here are the parameters of the system:

Rmin = 0.02      # m
Rmax = 0.055     # m
L = 47           # m
omega = 300      # rad / s

Rmin and Rmax are the initial and final values for the radius, r. L is the total length of the paper. omega is the angular velocity in radians per second.

Figuring out k is not easy, but we can estimate it by pretending that r is constant and equal to the average of Rmin and Rmax:

Ravg = (Rmax + Rmin) / 2

In that case, the circumference of the roll is also constant:

Cavg = 2 * np.pi * Ravg

And we can compute the number of revolutions to roll up length L, like this.

revs = L / Cavg

Converting rotations to radians, we can estimate the final value of theta.

theta = 2 * np.pi * revs
theta

Finally, k is the total change in r divided by the total change in theta.

k_est = (Rmax - Rmin) / theta
k_est

At the end of the chapter, we'll derive k analytically, but this estimate is enough to get started.

Simulating the System

The state variables we'll use are theta, y, and r. Here are the initial conditions:

init = State(theta=0, y=0, r=Rmin)

And here's a System object with init and t_end:

system = System(init=init, t_end=10)

Now we can use the differential equations from the previous section to write a slope function:

def slope_func(t, state, system):
    theta, y, r = state
    
    dydt = r * omega
    drdt = k_est * omega
    
    return omega, dydt, drdt

As usual, the slope function takes a time stamp, a State object, and a System object.

The job of the slope function is to compute the time derivatives of the state variables. The derivative of theta is angular velocity, omega. The derivatives of y and r are given by the differential equations we derived.

And as usual, we'll test the slope function with the initial conditions.

slope_func(0, system.init, system)

We'd like to stop the simulation when the length of paper on the roll is L. We can do that with an event function that passes through 0 when y equals L:

def event_func(t, state, system):
    theta, y, r = state
    return L - y

We can test it with the initial conditions:

event_func(0, system.init, system)

Now let's run the simulation:

results, details = run_solve_ivp(system, slope_func,
                                  events=event_func)
details.message

Here are the last few time steps.

results.tail()

The time it takes to complete one roll is about 4.2 seconds, which is consistent with what we see in the video.

results.index[-1]

The final value of y is 47 meters, as expected.

final_state = results.iloc[-1] 
final_state.y

The final value of r is 0.55 m, which is Rmax.

final_state.r

The total number of rotations is close to 200, which seems plausible.

radians = final_state.theta
rotations = radians / 2 / np.pi
rotations

As an exercise, we'll see how fast the paper is moving. But first, let's take a closer look at the results.

Plotting the Results

Here's what theta looks like over time.

def plot_theta(results):
    results.theta.plot(color='C0', label='theta')
    decorate(xlabel='Time (s)',
             ylabel='Angle (rad)')
    
plot_theta(results)

theta grows linearly, as we should expect with constant angular velocity.

Here's what r looks like over time.

def plot_r(results):
    results.r.plot(color='C2', label='r')

    decorate(xlabel='Time (s)',
             ylabel='Radius (m)')
    
plot_r(results)

r also increases linearly.

Here's what y looks like over time.

def plot_y(results):
    results.y.plot(color='C1', label='y')

    decorate(xlabel='Time (s)',
             ylabel='Length (m)')
    
plot_y(results)

Since the derivative of y depends on r, and r is increasing, y grows with increasing slope.

In the next section, we'll see that we could have solved these differential equations analytically. However, it is often useful to start with simulation as a way of exploring and checking assumptions.

Analytic Solution

Since angular velocity is constant:

$$\frac{d\theta}{dt} = \omega \quad\quad (1)$$

We can find $\theta$ as a function of time by integrating both sides:

$$\theta(t) = \omega t$$

Similarly, we can solve this equation

$$\frac{dr}{dt} = k \omega$$

to find

$$r(t) = k \omega t + R_{min}$$

Then we can plug the solution for $r$ into the equation for $y$:

$$\begin{aligned} \frac{dy}{dt} & = r \omega \quad\quad (2) \\ & = \left[ k \omega t + R_{min} \right] \omega \nonumber\end{aligned}$$

Integrating both sides yields:

$$y(t) = \left[ k \omega t^2 / 2 + R_{min} t \right] \omega$$

So $y$ is a parabola, as you might have guessed.

We can also use these equations to find the relationship between $y$ and $r$, independent of time, which we can use to compute $k$. Dividing Equations 1 and 2 yields

$$\frac{dr}{dy} = \frac{k}{r}$$

Separating variables yields

$$r~dr = k~dy$$

Integrating both sides yields

$$r^2 / 2 = k y + C$$

Solving for $y$, we have

ParseError: KaTeX parse error: Undefined control sequence: \label at position 44: … \̲l̲a̲b̲e̲l̲{eqn3}

When $y=0$, $r=R_{min}$, so

$$R_{min}^2 / 2 = C$$

When $y=L$, $r=R_{max}$, so

$$L = \frac{1}{2k} (R_{max}^2 - R_{min}^2)$$

Solving for $k$ yields

ParseError: KaTeX parse error: Undefined control sequence: \label at position 53: …}^2) \̲l̲a̲b̲e̲l̲{eqn4}

Plugging in the values of the parameters yields 2.8e-5 m/rad, the same as the "estimate" we computed in Section xxx.

k = (Rmax**2 - Rmin**2) / (2 * L)
k

In this case the estimate turns out to be exact.

Summary

This chapter introduces rotation, starting with an example where angular velocity is constant. We simulated the manufacture of a roll of toilet paper, then we solved the same problem analytically.

In the next chapter, we'll see a more interesting example where angular velocity is not constant. And we'll introduce three new concepts: torque, angular acceleration, and moment of inertia.

But first, you might want to work on the following exercise.

Exercises

This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. You can access the notebooks at https://allendowney.github.io/ModSimPy/.

Exercise 1

Since we keep omega constant, the linear velocity of the paper increases with radius. We can use gradient to estimate the derivative of results.y.

dydt = gradient(results.y)

Here's what the result looks like.

dydt.plot(label='dydt')
decorate(xlabel='Time (s)',
         ylabel='Linear velocity (m/s)')

With constant angular velocity, linear velocity is increasing, reaching its maximum at the end.

max_linear_velocity = dydt.iloc[-1]
max_linear_velocity

Now suppose this peak velocity is the limiting factor; that is, we can't move the paper any faster than that. In that case, we might be able to speed up the process by keeping the linear velocity at the maximum all the time.

Write a slope function that keeps the linear velocity, dydt, constant, and computes the angular velocity, omega, accordingly. Then, run the simulation and see how much faster we could finish rolling the paper.

# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here