# Chapter 24¶

Click here to run this chapter on Colab

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, in the same way that we sometimes treat velocity as a scalar with an implicit direction.

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, 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\):

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

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

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:

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

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

## 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 be 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
```

```
1253.3333333333335
```

Finally, `k`

is the total change in `r`

divided by the total change in `theta`

.

```
k_est = (Rmax - Rmin) / theta
k_est
```

```
2.7925531914893616e-05
```

At the end of the chapter, we’ll derive `k`

analytically, but this estimate is enough to get started.

## Simulation¶

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)
```

```
(300, 6.0, 0.008377659574468085)
```

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)
```

```
47.0
```

Now let’s run the simulation:

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

```
'A termination event occurred.'
```

Here are the last few time steps.

```
results.tail()
```

theta | y | r | |
---|---|---|---|

4.010667 | 1203.200000 | 44.277760 | 0.05360 |

4.052444 | 1215.733333 | 44.951740 | 0.05395 |

4.094222 | 1228.266667 | 45.630107 | 0.05430 |

4.136000 | 1240.800000 | 46.312860 | 0.05465 |

4.177778 | 1253.333333 | 47.000000 | 0.05500 |

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]
```

```
4.177777777777779
```

The final value of `y`

is 47 meters, as expected.

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

```
47.00000000000001
```

The final value of radius is 0.55 m, which is `Rmax`

.

```
final_state.r
```

```
0.05500000000000001
```

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

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

```
199.47419534184218
```

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

## Plotting¶

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.

But since the derivative of `y`

depends on `r`

, and `r`

is increasing, `y`

grows with increasing slope.

```
def plot_y(results):
results.y.plot(color='C1', label='y')
decorate(xlabel='Time (s)',
ylabel='Length (m)')
plot_y(results)
```

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

## Analysis¶

Since angular velocity is constant:

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

Similarly, we can solve this equation

to find

Then we can plug the solution for \(r\) into the equation for \(y\):

Integrating both sides yields:

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

Separating variables yields

Integrating both sides yields

Solving for \(y\), we have

When \(y=0\), \(r=R_{min}\), so

When \(y=L\), \(r=R_{max}\), so

Solving for \(k\) yields

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
```

```
2.7925531914893616e-05
```

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¶

### Exercise¶

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)
```

```
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
```

```
16.447499999999934
```

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.

Run the simulation and see how much faster we could finish rolling the paper.

```
# Solution
def slope_func(t, state, system):
theta, y, r = state
dydt = max_linear_velocity
omega = dydt / r
drdt = k * omega
return omega, dydt, drdt
```

```
# Solution
slope_func(0, system.init, system)
```

```
(822.3749999999967, 16.447499999999934, 0.022965259308510545)
```

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

```
'A termination event occurred.'
```

```
# Solution
t_final = results.index[-1]
t_final
```

```
2.8575771393828964
```

```
# Solution
plot_theta(results)
```

```
# Solution
plot_r(results)
```

```
# Solution
plot_y(results)
```