Lab 02-02: Simulating Free Fall and Impact

Objectives

Analyze energy conservation during free fall
Study impact dynamics and coefficient of restitution
Understand momentum transfer during collisions
Compare simulation results with analytical predictions

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: Energy Conservation Analysis

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

# Create bouncing ball model
ball_xml = """
<mujoco>
  <option gravity="0 0 -9.81" timestep="0.0005"/>

  <default>
    <geom condim="1" solimp="0.9 0.95 0.001" solref="0.02 1"/>
  </default>

  <worldbody>
    <light diffuse="1 1 1" pos="0 0 5"/>
    <geom type="plane" size="5 5 0.1" rgba="0.8 0.8 0.8 1"/>

    <body name="ball" pos="0 0 3">
      <joint type="free"/>
      <geom type="sphere" size="0.1" mass="1.0" rgba="1 0.2 0.2 1"/>
    </body>
  </worldbody>
</mujoco>
"""

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

def compute_energy(model, data):
    """Compute kinetic and potential energy."""
    mass = model.body_mass[1]
    height = data.qpos[2]
    vel = np.linalg.norm(data.qvel[:3])

    kinetic = 0.5 * mass * vel**2
    potential = mass * 9.81 * height

    return kinetic, potential

# Simulate and record energy
mujoco.mj_resetData(model, data)

times = []
kinetic_energy = []
potential_energy = []
total_energy = []
heights = []

for _ in range(4000):  # 2 seconds
    mujoco.mj_step(model, data)

    ke, pe = compute_energy(model, data)
    times.append(data.time)
    kinetic_energy.append(ke)
    potential_energy.append(pe)
    total_energy.append(ke + pe)
    heights.append(data.qpos[2])

# Plot energy over time
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

axes[0].plot(times, heights)
axes[0].set_ylabel('Height (m)')
axes[0].set_title('Ball Trajectory')
axes[0].grid(True)

axes[1].plot(times, kinetic_energy, label='Kinetic', alpha=0.7)
axes[1].plot(times, potential_energy, label='Potential', alpha=0.7)
axes[1].plot(times, total_energy, label='Total', linewidth=2)
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Energy (J)')
axes[1].set_title('Energy vs Time')
axes[1].legend()
axes[1].grid(True)

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

print(f"Initial total energy: {total_energy[0]:.3f} J")
print(f"Final total energy: {total_energy[-1]:.3f} J")
print(f"Energy loss: {(1 - total_energy[-1]/total_energy[0])*100:.1f}%")

Expected: Total energy should decrease at each bounce due to inelastic collision. Analysis: Note how energy is lost during contact but conserved during flight.

Step 2: Coefficient of Restitution

def measure_restitution(model, data, drop_height):
    """Measure coefficient of restitution from drop test."""
    mujoco.mj_resetData(model, data)
    data.qpos[2] = drop_height

    bounce_heights = []
    prev_height = drop_height
    going_down = True
    min_height_seen = drop_height

    for _ in range(10000):
        mujoco.mj_step(model, data)

        height = data.qpos[2]
        vel_z = data.qvel[2]

        # Detect bounce (transition from going down to going up)
```python
```python
        if going_down and vel_z > 0.1 and height < prev_height + 0.01:

            going_down = False
            min_height_seen = height

        # Detect peak (transition from going up to going down)
```python
```python
        if not going_down and vel_z < -0.1:

            bounce_heights.append(prev_height)
            going_down = True

```python
```python
            if len(bounce_heights) >= 5:

                break

        prev_height = height

    # Compute coefficient of restitution
    # e = sqrt(h_after / h_before) for each bounce
    restitutions = []
    heights_list = [drop_height] + bounce_heights

    for i in range(len(heights_list) - 1):
```python
```python
        if heights_list[i] > 0.01:  # Avoid division issues

            e = np.sqrt(heights_list[i+1] / heights_list[i])
            restitutions.append(e)

    return heights_list, restitutions

# Test with different drop heights
drop_heights = [1.0, 2.0, 3.0]
results = {}

for h in drop_heights:
    heights, e_values = measure_restitution(model, data, h)
    results[h] = (heights, e_values)
    print(f"\nDrop height {h}m:")
    print(f"  Bounce heights: {[f'{hh:.3f}' for hh in heights]}")
    print(f"  Restitution coefficients: {[f'{e:.3f}' for e in e_values]}")
    print(f"  Average e: {np.mean(e_values):.3f}")

# Plot bounce decay
fig, ax = plt.subplots(figsize=(10, 6))

for h, (heights, _) in results.items():
    bounces = range(len(heights))
    ax.plot(bounces, heights, 'o-', label=f'Drop from {h}m')

ax.set_xlabel('Bounce Number')
ax.set_ylabel('Height (m)')
ax.set_title('Bounce Height Decay')
ax.legend()
ax.grid(True)
plt.savefig('bounce_decay.png', dpi=150)
plt.show()

Analysis Questions:

Is the coefficient of restitution constant across bounces?
Does it depend on drop height (impact velocity)?

Step 3: Two-Body Collision

# Create two-ball collision model
collision_xml = """
<mujoco>
  <option gravity="0 0 0" timestep="0.0001"/>

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

    <body name="ball1" pos="-1 0 1">
      <joint type="free"/>
      <geom type="sphere" size="0.1" mass="1.0" rgba="1 0.2 0.2 1"/>
    </body>

    <body name="ball2" pos="1 0 1">
      <joint type="free"/>
      <geom type="sphere" size="0.1" mass="1.0" rgba="0.2 0.2 1 1"/>
    </body>
  </worldbody>
</mujoco>
"""

model2 = mujoco.MjModel.from_xml_string(collision_xml)
data2 = mujoco.MjData(model2)

# Set initial velocities (head-on collision)
mujoco.mj_resetData(model2, data2)
data2.qvel[0] = 2.0   # ball1 moving right
data2.qvel[6] = -1.0  # ball2 moving left

# Record before collision
m1 = model2.body_mass[1]
m2 = model2.body_mass[2]
v1_before = data2.qvel[0]
v2_before = data2.qvel[6]

p_before = m1 * v1_before + m2 * v2_before
ke_before = 0.5 * m1 * v1_before**2 + 0.5 * m2 * v2_before**2

print("Before Collision:")
print(f"  Ball 1: v = {v1_before:.2f} m/s")
print(f"  Ball 2: v = {v2_before:.2f} m/s")
print(f"  Total momentum: {p_before:.2f} kg⋅m/s")
print(f"  Total KE: {ke_before:.2f} J")

# Record trajectory
times = []
x1_positions = []
x2_positions = []
v1_velocities = []
v2_velocities = []
momenta = []

for _ in range(5000):
    mujoco.mj_step(model2, data2)

    times.append(data2.time)
    x1_positions.append(data2.qpos[0])
    x2_positions.append(data2.qpos[7])
    v1_velocities.append(data2.qvel[0])
    v2_velocities.append(data2.qvel[6])
    momenta.append(m1 * data2.qvel[0] + m2 * data2.qvel[6])

# After collision
v1_after = v1_velocities[-1]
v2_after = v2_velocities[-1]
p_after = momenta[-1]
ke_after = 0.5 * m1 * v1_after**2 + 0.5 * m2 * v2_after**2

print("\nAfter Collision:")
print(f"  Ball 1: v = {v1_after:.2f} m/s")
print(f"  Ball 2: v = {v2_after:.2f} m/s")
print(f"  Total momentum: {p_after:.2f} kg⋅m/s")
print(f"  Total KE: {ke_after:.2f} J")
print(f"\nMomentum conservation: {abs(p_after - p_before) < 0.01}")
print(f"Energy loss: {(ke_before - ke_after):.3f} J ({(1-ke_after/ke_before)*100:.1f}%)")

# Plot collision
fig, axes = plt.subplots(3, 1, figsize=(10, 10))

axes[0].plot(times, x1_positions, label='Ball 1')
axes[0].plot(times, x2_positions, label='Ball 2')
axes[0].set_ylabel('X Position (m)')
axes[0].legend()
axes[0].grid(True)

axes[1].plot(times, v1_velocities, label='Ball 1')
axes[1].plot(times, v2_velocities, label='Ball 2')
axes[1].set_ylabel('X Velocity (m/s)')
axes[1].legend()
axes[1].grid(True)

axes[2].plot(times, momenta)
axes[2].axhline(y=p_before, color='r', linestyle='--', label='Initial')
axes[2].set_xlabel('Time (s)')
axes[2].set_ylabel('Total Momentum (kg⋅m/s)')
axes[2].legend()
axes[2].grid(True)

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

Verification: Momentum should be conserved. Energy may be lost (inelastic collision). Analysis: Compare with theoretical predictions for elastic/inelastic collisions.

Step 4: Vary Collision Parameters

def simulate_collision(e_coef, v1_init, v2_init):
    """Simulate collision with given restitution coefficient."""
    # Modify contact parameters
    xml = f"""
    <mujoco>
      <option gravity="0 0 0" timestep="0.0001"/>
      <default>
        <geom condim="1" solimp="{e_coef} 0.95 0.001" solref="0.001 1"/>
      </default>
      <worldbody>
        <body name="ball1" pos="-0.5 0 0">
          <joint type="free"/>
          <geom type="sphere" size="0.1" mass="1.0"/>
        </body>
        <body name="ball2" pos="0.5 0 0">
          <joint type="free"/>
          <geom type="sphere" size="0.1" mass="2.0"/>
        </body>
      </worldbody>
    </mujoco>
    """

    m = mujoco.MjModel.from_xml_string(xml)
    d = mujoco.MjData(m)

    d.qvel[0] = v1_init
    d.qvel[6] = v2_init

    # Run until well after collision
    for _ in range(5000):
        mujoco.mj_step(m, d)

    return d.qvel[0], d.qvel[6]

# Compare with analytical solution for 1D elastic collision
# v1' = ((m1-m2)*v1 + 2*m2*v2) / (m1+m2)
# v2' = ((m2-m1)*v2 + 2*m1*v1) / (m1+m2)

m1, m2 = 1.0, 2.0
v1, v2 = 3.0, -1.0

# Analytical (perfectly elastic)
v1_elastic = ((m1-m2)*v1 + 2*m2*v2) / (m1+m2)
v2_elastic = ((m2-m1)*v2 + 2*m1*v1) / (m1+m2)

print("Analytical Solution (elastic collision):")
print(f"  v1' = {v1_elastic:.3f} m/s")
print(f"  v2' = {v2_elastic:.3f} m/s")

# Simulate with different restitution values
print("\nSimulated Results:")
for e in [0.5, 0.7, 0.9, 0.99]:
    v1_sim, v2_sim = simulate_collision(e, v1, v2)
    print(f"  e={e}: v1'={v1_sim:.3f}, v2'={v2_sim:.3f}")

Analysis: As restitution approaches 1.0, results should approach elastic collision.

Expected Outcomes

Energy analysis showing conservation during flight, loss during contact
Coefficient of restitution measurements
Momentum conservation verification
Comparison with analytical collision predictions

Rubric

Criterion	Points	Description
Energy Analysis	25	Correctly tracks and analyzes energy
Restitution Measurement	25	Accurately measures coefficient
Collision Dynamics	25	Momentum/energy analysis
Parameter Study	15	Systematic variation of parameters
Visualization	10	Clear, informative plots

Total: 100 points

Objectives​

Materials​

Instructions​

Step 1: Energy Conservation Analysis​

Step 2: Coefficient of Restitution​

Step 3: Two-Body Collision​

Step 4: Vary Collision Parameters​

Expected Outcomes​

Rubric​