Lab 02-03: Torque-Motion Relationship

Objectives

Understand the relationship τ = Iα for rotational motion
Measure moment of inertia through simulation
Analyze gyroscopic effects
Study coupled rotational and translational dynamics

Materials

Type	Name	Version	Tier
software	Python	3.10+	required
software	MuJoCo	3.0+	required
software	NumPy	1.24+	required
software	Matplotlib	3.5+	required

Instructions

Step 1: Single-Axis Rotation Under Constant Torque

import mujoco
import numpy as np
import matplotlib.pyplot as plt

# Create a rotating bar model
bar_xml = """
<mujoco>
  <option gravity="0 0 0" timestep="0.001"/>

  <worldbody>
    <light diffuse="1 1 1" pos="0 0 3"/>

    <body name="bar" pos="0 0 1">
      <joint type="free"/>
      <geom type="box" size="0.5 0.05 0.05" mass="2.0" rgba="0.8 0.3 0.2 1"/>
    </body>
  </worldbody>
</mujoco>
"""

model = mujoco.MjModel.from_xml_string(bar_xml)
data = mujoco.MjData(model)

# Calculate theoretical moment of inertia
# For a rectangular bar rotating about its center (z-axis):
# I_zz = (1/12) * m * (l_x² + l_y²)
mass = model.body_mass[1]
lx, ly = 1.0, 0.1  # Full length and width (2 * half-size)
I_theoretical = (1/12) * mass * (lx**2 + ly**2)
print(f"Theoretical I_zz: {I_theoretical:.6f} kg⋅m²")

# Apply constant torque and measure angular acceleration
mujoco.mj_resetData(model, data)

torque_values = [0.5, 1.0, 2.0, 4.0]
results = {}

for tau in torque_values:
    mujoco.mj_resetData(model, data)

    times = []
    angular_velocities = []

    for i in range(1000):
        # Apply torque about z-axis (index 5 in 6D wrench)
        data.xfrc_applied[1] = [0, 0, 0, 0, 0, tau]
        mujoco.mj_step(model, data)

        times.append(data.time)
        angular_velocities.append(data.qvel[5])  # z angular velocity

    # Measure angular acceleration (slope of angular velocity)
    alpha_measured = np.polyfit(times, angular_velocities, 1)[0]
    I_measured = tau / alpha_measured

    results[tau] = {
        'times': times,
        'omega': angular_velocities,
        'alpha': alpha_measured,
        'I': I_measured
    }

    print(f"τ = {tau:.1f} Nm: α = {alpha_measured:.3f} rad/s², I = {I_measured:.6f} kg⋅m²")

# Plot angular velocity vs time for different torques
plt.figure(figsize=(10, 6))
for tau, data_dict in results.items():
    plt.plot(data_dict['times'], data_dict['omega'], label=f'τ = {tau} Nm')

plt.xlabel('Time (s)')
plt.ylabel('Angular Velocity (rad/s)')
plt.title('Angular Velocity Under Different Torques')
plt.legend()
plt.grid(True)
plt.savefig('torque_response.png', dpi=150)
plt.show()

# Average measured inertia
avg_I = np.mean([r['I'] for r in results.values()])
print(f"\nAverage measured I: {avg_I:.6f} kg⋅m²")
print(f"Error from theoretical: {abs(avg_I - I_theoretical)/I_theoretical * 100:.2f}%")

Verification: Measured I should be close to theoretical value (~0.168 kg⋅m²). Understanding: α = τ/I, so doubling torque doubles angular acceleration.

Step 2: Pendulum Dynamics

# Create a simple pendulum
pendulum_xml = """
<mujoco>
  <option gravity="0 0 -9.81" timestep="0.001"/>

  <worldbody>
    <light diffuse="1 1 1" pos="0 0 3"/>

    <body name="pivot" pos="0 0 2">
      <joint name="hinge" type="hinge" axis="0 1 0" limited="false"/>
      <geom type="sphere" size="0.05" mass="0.1"/>

      <body name="pendulum" pos="0 0 -0.5">
        <geom type="capsule" size="0.02" fromto="0 0 0.5 0 0 -0.5" mass="1.0"/>
      </body>
    </body>
  </worldbody>
</mujoco>
"""

model_p = mujoco.MjModel.from_xml_string(pendulum_xml)
data_p = mujoco.MjData(model_p)

# Set initial angle
initial_angle = np.radians(30)  # 30 degrees
mujoco.mj_resetData(model_p, data_p)
data_p.qpos[0] = initial_angle

# Simulate and record motion
times = []
angles = []
angular_velocities = []

for _ in range(5000):  # 5 seconds
    mujoco.mj_step(model_p, data_p)
    times.append(data_p.time)
    angles.append(data_p.qpos[0])
    angular_velocities.append(data_p.qvel[0])

# Theoretical period for small oscillations
# T = 2π√(L/g) for simple pendulum
L = 1.0  # pendulum length
g = 9.81
T_theoretical = 2 * np.pi * np.sqrt(L / g)

# Measure period from simulation
from scipy.signal import find_peaks
peaks, _ = find_peaks(angles, height=0)
```python
```python
if len(peaks) > 1:

    periods = np.diff([times[p] for p in peaks])
    T_measured = np.mean(periods)
else:
    T_measured = 0

print(f"Theoretical period: {T_theoretical:.4f} s")
print(f"Measured period: {T_measured:.4f} s")
print(f"Error: {abs(T_measured - T_theoretical) / T_theoretical * 100:.2f}%")

# Phase space plot
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].plot(times, np.degrees(angles))
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Angle (degrees)')
axes[0].set_title('Pendulum Motion')
axes[0].grid(True)

axes[1].plot(np.degrees(angles), angular_velocities)
axes[1].set_xlabel('Angle (degrees)')
axes[1].set_ylabel('Angular Velocity (rad/s)')
axes[1].set_title('Phase Space')
axes[1].grid(True)

plt.tight_layout()
plt.savefig('pendulum_dynamics.png', dpi=150)
plt.show()

Expected: Period should be approximately 2.0 seconds for L=1m. Analysis: Phase space should show elliptical trajectories for undamped motion.

Step 3: Coupled Rotation - Spinning Top

# Create a spinning top (gyroscope)
top_xml = """
<mujoco>
  <option gravity="0 0 -9.81" timestep="0.0005"/>

  <worldbody>
    <light diffuse="1 1 1" pos="0 0 3"/>
    <geom type="plane" size="2 2 0.1"/>

    <body name="top" pos="0 0 0.5">
      <joint type="free"/>
      <geom type="cylinder" size="0.1 0.02" mass="0.5" pos="0 0 0"/>
      <geom type="cylinder" size="0.01 0.3" mass="0.1" pos="0 0 0.15"/>
    </body>
  </worldbody>
</mujoco>
"""

model_t = mujoco.MjModel.from_xml_string(top_xml)
data_t = mujoco.MjData(model_t)

# Initialize with tilt and spin
mujoco.mj_resetData(model_t, data_t)

# Tilt the top (modify quaternion for ~10 degree tilt)
tilt_angle = np.radians(10)
data_t.qpos[3:7] = [np.cos(tilt_angle/2), np.sin(tilt_angle/2), 0, 0]

# Give it spin about its axis
spin_rate = 50.0  # rad/s
data_t.qvel[5] = spin_rate  # Spin about z-axis (before tilt)

# Record motion
times = []
x_coords = []
y_coords = []
spin_rates = []

for i in range(10000):  # 5 seconds
    mujoco.mj_step(model_t, data_t)

```python
```python
    if i % 10 == 0:

        times.append(data_t.time)
        x_coords.append(data_t.qpos[0])
        y_coords.append(data_t.qpos[1])
        spin_rates.append(np.linalg.norm(data_t.qvel[3:6]))

# Plot precession (top view of center of mass trajectory)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].plot(x_coords, y_coords)
axes[0].set_xlabel('X Position (m)')
axes[0].set_ylabel('Y Position (m)')
axes[0].set_title('Top View - Precession Pattern')
axes[0].axis('equal')
axes[0].grid(True)

axes[1].plot(times, spin_rates)
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Spin Rate (rad/s)')
axes[1].set_title('Spin Rate Over Time')
axes[1].grid(True)

plt.tight_layout()
plt.savefig('spinning_top.png', dpi=150)
plt.show()

print(f"Initial spin rate: {spin_rate:.1f} rad/s")
print(f"Final spin rate: {spin_rates[-1]:.1f} rad/s")

Observation: The top should exhibit precession - circular motion of its axis. Analysis: Precession rate depends on spin rate, mass distribution, and tilt.

Step 4: Rotational Inertia Tensor

# Compare rotation about different axes
shapes_xml = """
<mujoco>
  <option gravity="0 0 0" timestep="0.001"/>

  <worldbody>
    <light diffuse="1 1 1" pos="0 0 5"/>

    <body name="long_bar" pos="-2 0 1">
      <joint type="free"/>
      <geom type="box" size="0.5 0.05 0.05" mass="1.0" rgba="1 0 0 1"/>
    </body>

    <body name="disc" pos="0 0 1">
      <joint type="free"/>
      <geom type="cylinder" size="0.2 0.02" mass="1.0" rgba="0 1 0 1"/>
    </body>

    <body name="cube" pos="2 0 1">
      <joint type="free"/>
      <geom type="box" size="0.15 0.15 0.15" mass="1.0" rgba="0 0 1 1"/>
    </body>
  </worldbody>
</mujoco>
"""

model_s = mujoco.MjModel.from_xml_string(shapes_xml)
data_s = mujoco.MjData(model_s)

def measure_inertia_by_axis(body_idx, axis, torque_mag=1.0):
    """Measure moment of inertia about specified axis."""
    mujoco.mj_resetData(model_s, data_s)

    # Construct torque vector
    torque = [0, 0, 0, 0, 0, 0]
    torque[3 + axis] = torque_mag  # axis: 0=x, 1=y, 2=z

    times = []
    angular_vels = []

    for i in range(500):
        data_s.xfrc_applied[body_idx] = torque
        mujoco.mj_step(model_s, data_s)

        times.append(data_s.time)
        # Get angular velocity for this body
        qvel_start = (body_idx - 1) * 6 + 3  # Skip linear velocities
        angular_vels.append(data_s.qvel[qvel_start + axis])

    alpha = np.polyfit(times, angular_vels, 1)[0]
    return torque_mag / alpha

# Measure inertia for each shape about each axis
shapes = ['Long Bar', 'Disc', 'Cube']
body_indices = [1, 2, 3]
axes = ['X', 'Y', 'Z']

print("Measured Moments of Inertia (kg⋅m²):")
print("-" * 50)
print(f"{'Shape':<12} {'I_x':>10} {'I_y':>10} {'I_z':>10}")
print("-" * 50)

for name, body_idx in zip(shapes, body_indices):
    inertias = [measure_inertia_by_axis(body_idx, ax) for ax in range(3)]
    print(f"{name:<12} {inertias[0]:>10.6f} {inertias[1]:>10.6f} {inertias[2]:>10.6f}")

Analysis:

Long bar: I_x ≈ I_y << I_z (easy to spin about long axis)
Disc: I_x ≈ I_y > I_z (like a frisbee)
Cube: I_x ≈ I_y ≈ I_z (symmetric)

Expected Outcomes

Verified τ = Iα relationship
Measured pendulum period matching theory
Observed gyroscopic precession
Compared inertia tensors for different shapes

Rubric

Criterion	Points	Description
Torque-Acceleration	25	Correct measurement and analysis
Pendulum Dynamics	25	Period measurement and phase space
Gyroscopic Effects	25	Precession observation and analysis
Inertia Tensor	15	Comparison across shapes and axes
Documentation	10	Clear explanations

Total: 100 points

Objectives​

Materials​

Instructions​

Step 1: Single-Axis Rotation Under Constant Torque​

Step 2: Pendulum Dynamics​

Step 3: Coupled Rotation - Spinning Top​

Step 4: Rotational Inertia Tensor​

Expected Outcomes​

Rubric​