Gradient

Our approach is to apply Newton’s second law in a tangential and normal basis (TNB basis or Serret–Frenet frame) with unit vectors u_t and u_n, as shown in figure 5.17 [@hibbeler2015, § 13.5]. The advantage of using this coordinate system is that the variable (normal) force applied to the bead by the hoop is decoupled from the force along the (tangential) trajectory of the bead.

The position vector r has its tail at the origin O of the Cartesian coordinate system and its head at the center of the bead. The easiest parameter to describe the circular curve with is the angle θ from the positive x axis. This parametric representation of r is, in Cartesian coordinates, $$ \bm{r}(\theta) = \begin{bmatrix} R \cos\theta \\ R \sin\theta \end{bmatrix}. $$ However, when working with a TNB basis, we should parameterize the position vector with the arclength s [@oreilly2019,§ 3.1.2]. In general, the process of changing parameters from parameter θ to arclength s generally proceeds as follows. Compute the arclength s from its definition [@kreyszig2011,p. 385], for some a ∈ ℝ, $$ s = \int_a^\theta \sqrt{\bm{r}'(\tilde{\theta}) \cdot \bm{r}'(\tilde{\theta})}\ d\tilde{\theta},$$ then solve for θ(s). The arclength parameter s can now replace θ by substitution, r(θ) → r(θ(s)) → r(s).

However in our case the (circular) arclength has the well-known relationship with the corresponding angle, s = Rθ. Solving for θ(s), θ(s) = s/R. So the arclength parameterization of the position vector is $$ \bm{r}(s) = \begin{bmatrix} R \cos(s/R) \\ R \sin(s/R) \end{bmatrix}. $$

From here, we proceed in Python. Begin by loading the necessary packages.

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, cumulative_trapezoid

Define a symbolic parametric representation of the position vector r(s) in Cartesian coordinates. In our code, we will use T instead of θ.

R = sp.symbols("R", positive=True)  # Radius of the hoop
s = sp.symbols("s", real=True)  # Arc length
r_s = sp.Matrix([R * sp.cos(s/R), R * sp.sin(s/R)])

Now we can compute the gravitational force in the tangential direction. The gravitational force is $\bm{F}_g = -m g \hat{\bm{j}}$, where m is the mass of the bead. The unit tangent vector u_t can be conveniently computed from the arclength parameterization of r as $$\bm{u}_t = \frac{d\bm{r}}{ds} = \begin{bmatrix}-\sin(s/R) \\ \cos(s/R)\end{bmatrix}.$$ The gravitational force in the tangential direction is then $$\bm{F}_g(s) \cdot \bm{u}_t = \bm{F}_g \cdot \frac{d\bm{r}}{ds} = \begin{bmatrix}0 \\ -m g\end{bmatrix} \cdot \begin{bmatrix}-\sin(s/R) \\ \cos(s/R)\end{bmatrix} = -m g \cos(s/R).$$ Computing this expression symbolically:

m, g = sp.symbols("m, g", positive=True) # Mass, grav. acceleration
F_g = sp.Matrix([0, -m * g])  # Grav. force (Cartesian coordinates)
u_t = r_s.diff(s)  # Unit tangent vector (Cartesian coordinates)
    # This has unit length because r_s is parameterized by arc length
F_g_t = F_g.dot(u_t)  # Grav. force in the tangential direction

Convert the symbolic expression to a Python function for numerical evaluation.

params = {R: 1, m: 1, g: 9.81}  # Parameters
F_g_t_fun = sp.lambdify((s), F_g_t.subs(params), "numpy")

Apply Newton’s second law to the mass to form a dynamic system state space model. The state vector is $\begin{bmatrix}s & v\end{bmatrix}^\top$, where s is the arc length and v is the velocity in the tangential direction. To define the state space model, we need to compute the time derivatives of the state variables. By definition, the velocity is the rate of change of arc length with respect to time, so $$\frac{ds}{dt} = v.$$ So this is the first equation in our system of ODEs comprising the state space model. The tangential acceleration is defined as the rate of change of velocity with respect to time, $$a = \frac{dv}{dt}.$$ Applying Newton’s second law to the bead in the tangential direction, we have ma = F_g(s) ⋅ u_t =  − mgcos (s/R). Substituting the definition of acceleration and the gravitational force, we have $$\frac{dv}{dt} = -g \cos(s/R).$$ This is the second equation in our state space model. Define the state space model in Python:

def system(t, X, params):
    """Bead on a hoop dynamic system, sans damping."""
    s, v = X  # Unpack the state vector
    ds_dt = v  # Velocity in the tangential direction
    dv_dt = 1/m.subs(params) * F_g_t_fun(s)
        # Acceleration in the tangential direction
    return [ds_dt, dv_dt]

Numerically solve the system of ODEs for an initial state vector [s₀,v₀].

X0 = [sp.pi/3, 0]  # Initial state vector [s0, v0]
t = np.linspace(0, 10, 201)  # Time array
X_t = solve_ivp(
    system, (t[0], t[-1]), 
    X0, t_eval=t, args=(params,)
).y  # Solve the system

Now write a function to plot the trajectory of the bead on the hoop through time.

def plot_trajectory(t, X_t):
    """Plot the trajectory of the bead on the hoop through time."""
    fig, ax = plt.subplots(2, 1, sharex=True)
    ax[0].plot(t, X_t[0])
    ax[0].set_ylabel("Arc length $s$ (m)")
    ax[1].plot(t, X_t[1])
    ax[1].set_ylabel("Velocity $v$ (m/s)")
    ax[1].set_xlabel("Time $t$ (s)")
    return fig

Plot the trajectory of the bead on the hoop through time.

fig = plot_trajectory(t, X_t)
plt.draw()

Now write a function to plot the trajectory of the bead on the hoop through space and time.

def plot_trajectory_space_time(t, X_t, r_s):
    """Plot the trajectory of the bead through space and time"""
    # First define a function to compute the position vector r(s) 
    # from the state variable arclength s 
    r_s_fun = sp.lambdify((s), r_s.subs(params), "numpy")
    r_s_t = r_s_fun(X_t[0])  # Position vector r(s) for each s in X_t
    fig = plt.figure()
    ax = fig.add_subplot(projection="3d")
    t_normalized = (t - np.min(t)) / (1.2*np.max(t) - np.min(t))
    for i in range(len(t) - 1): # Each segment with a different color
        ax.plot(
            r_s_t[0][0][i:i+2], 
            r_s_t[1][0][i:i+2], 
            [t[i], t[i+1]], 
            color=plt.cm.viridis(t_normalized[i])
        )
    ticks = np.array([-1, 0, 1])*params[R]
    ax.set_xticks(ticks)
    ax.set_yticks(ticks)
    ax.set_xlabel("$x$ (m)")
    ax.set_ylabel("$y$ (m)")
    ax.set_zlabel("Time $t$ (s)")
    ax.view_init(elev=150, azim=0, roll=90)
    return fig

Plot the trajectory of the bead on the hoop through space and time.

fig = plot_trajectory_space_time(t, X_t, r_s)
plt.draw()

Now add a viscous damping force in the anti-tangential direction.

b = sp.symbols("b", positive=True)  # Damping coefficient
v = sp.symbols("v", real=True)  # Velocity
F_d_t = -b * v  # Damping force in the tangential direction
params[b] = 0.1  # Damping coefficient
F_d_t_fun = sp.lambdify((v), F_d_t.subs(params), "numpy")  # Numerical

Apply Newton’s second law to the mass to form a dynamic system state space model.

def system_damped(t, X, params):
    """Bead on a hoop dynamic system, with damping."""
    s, v = X  # Unpack the state vector
    ds_dt = v  # Velocity in the tangential direction
    dv_dt = 1/m.subs(params) * (F_g_t_fun(s) + F_d_t_fun(v))
        # Acceleration in the tangential direction
    return [ds_dt, dv_dt]

Numerically solve the system of ODEs for the same initial state.

X_t_damped = solve_ivp(
    system_damped, (t[0], t[-1]), 
    X0, t_eval=t, args=(params,)
).y  # Solve the system

Plot the trajectory of the bead on the hoop through time with damping.

fig = plot_trajectory(t, X_t_damped)
plt.draw()

Plot the trajectory of the bead on the hoop through space and time with damping.

fig = plot_trajectory_space_time(t, X_t_damped, r_s)
plt.draw()

Now consider the work done by the gravitational and damping forces on the bead as t → ∞. From the state space model, in steady state, a = v = 0. This corresponds to the bead at rest at the bottom of the hoop, with s =  − Rπ/2.

As we have already established, the gravitational force is conservative, so the work done by it is independent of the path taken. The work done by the damping force is given by W_d = ∫_s0^s₁F_d ⋅ u_t ds. Substituting the expression for the damping force and the unit tangent vector, we have W_d = ∫_s0^s₁ − bv ds =  − b∫_s0^s₁v ds. Or, changing the variable of integration to time via the substitution ds = v dt, W_d =  − b∫_t0^t₁v² dt. We can evaluate this integral numerically by simulating the system for a relatively long time to get v(t).

t_long = np.linspace(0, 100, 3001)  # Long time array
X_t_long = solve_ivp(
    system_damped, (t_long[0], t_long[-1]), 
    X0, t_eval=t_long, args=(params,)
).y  # Solve the system

Compute the work done by the damping force.

W_d = -params[b] * \
    cumulative_trapezoid(X_t_long[1]**2, t_long, initial=0)
print(f"The work done by the damping force is approximately "
      f"{W_d[-1]:.2f} J.")

The work done by the damping force is approximately -17.91 J.

Plot the work done by the damping force through time.

fig, ax = plt.subplots()
ax.plot(t_long, W_d)
ax.set_xlabel("Time $t$ (s)")
ax.set_ylabel("Damping Force Work $W_d$ (J)")
plt.show()

Now compute the work done by the gravitational force. Begin by computing the height of the bead at the initial and final states as follows:

h0 = R * sp.sin(X0[0]/R)  # Height, initial state
h1 = R * sp.sin(-sp.pi/2)  # Height, final state

If the damping force was not present and the bead started from rest (T₀ = 0), in the downward position the kinetic energy would be equal to the difference in potential energy from the initial position; that is, T₁ = U₀ − U₁. The difference in kinetic energy from the initial position would be equal to the work done by the gravitational force, W_g = T₁ − T₀ = U₀ − U₁ = mgh₀ − mgh₁ = mg(h₀−h₁). Regardless of the presence of the damping force or the path taken, by the work-energy theorem the work done by the gravitational force is equal to this change in kinetic energy of the bead. Let’s compute this work done by the gravitational force.

W_g = (m * g * (h0 - h1)).evalf(subs=params)

18.3057092111253

We expect that over time the work done by the damping force will be equal in magnitude to the work done by the gravitational force in moving the bead from its initial position to the bottom of the hoop.