Simple examples of jump modeling with power and work breakdowns

simple_jump.py

# This is a simple point mass "jumping" simulation. Given a point mass of a
# realistic human mass that starts some height above the ground with a force
# that acts vertically between the ground and mass that represents the force
# that legs can generate, we simulate from standstill to maximum height in the
# air and calculate the work done by all forces on the mass as well as the
# potential and kinetic energy at all times in the simulation. There is also a
# Coulomb friction term that can represent some dissapative resistance to the
# motion.
#
#   O   |g
#   ^   v
#   |
#   Fl
#   |
#   v
# ------
# //////
#
#
# License: MIT 2022 Jason K. Moore

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np

m = 75.0  # mass of human, kg
g = 9.81  # acceleration due to gravity, m/s**2
Fl_max = 1200.0  # maximum force produced by the legs, N
Ff_mag = 200.0  # magnitude of the Coulomb friction, N


# friction always opposes the motion
@np.vectorize
def calc_friction_force(t, v):
    return -np.sign(v)*Ff_mag


@np.vectorize
def calc_gravity_force(t):
    return -m*g


@np.vectorize
def calc_leg_force(t):
    if t < 0.2:
        Fl = 0.0
    elif t > 1.0:
        Fl = 0.0
    else:
        Fl = Fl_max
    return Fl


def rhs(t, s):

    v = s[1]

    Fg = calc_gravity_force(t)
    Fl = calc_leg_force(t)
    Ff = calc_friction_force(t, v)

    Pg = Fg*v  # gravitational power
    Pl = Fl*v  # leg power
    Pf = Ff*v  # friction power
    P = Pg + Pl + Pf  # total power
    Pl_abs = np.abs(Pl)  # sort of represents metabolic power needed

    ydot = v
    vdot = (Fl + Fg + Ff)/m  # a = F/m

    return ydot, vdot, Pg, Pl, P, Pl_abs, Pf


# stop the simulation when the jumper reaches maximum height
def at_max_height(t, s):
    return s[1] - 1e-12


at_max_height.terminal = True
at_max_height.direction = -1.0  # don't stop at the squat position

t0 = 0.0
tf = 2.0
y0 = 1.0  # starting height of the center of mass

# simulate to find the max height
res = solve_ivp(rhs,
                (t0, tf),
                [y0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                events=at_max_height,
                atol=1e-12,
                rtol=1e-12)
# simulate again to have a fine mesh over the full duration
tf = res.t[-1]
res = solve_ivp(rhs,
                (t0, tf),
                [y0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                t_eval=np.linspace(t0, tf, num=1001),
                atol=1e-12,
                rtol=1e-12)

# extract all the simulation results
y = res.y[0]  # mass height
v = res.y[1]  # mass speed
Wg = res.y[2]  # work done by the gravity force on the mass
Wl = res.y[3]  # work done by the leg force on the mass
W = res.y[4]  # work done by all forces acting on the mass
Wl_abs = res.y[5]  # work done by leg, sum of + and -
Wf = res.y[6]  # work done by friction

# calculate the various power terms (same as in rhs())
leg_power = calc_leg_force(res.t)*v
gravity_power = calc_gravity_force(res.t)*v
friction_power = calc_friction_force(res.t, v)*v

# calculate the energy at each time
potential_energy = np.abs(calc_gravity_force(res.t))*y
kinetic_energy = m*v**2/2
total_energy = potential_energy + kinetic_energy

E_p0 = potential_energy[0]  # potential energy when standing still (above ground)
E_pmin = np.min(potential_energy)  # potential energy at lowest squat
E_pmax = np.max(potential_energy)  # potential energy at high point in the jump

print('Change in potential energy from initial to max height: ', E_pmax - E_p0)
print('Total positive work done by by all forces', W[-1])
print('Total positive work done by leg', Wl[-1])
print('Total positive work done by gravity', Wg[-1])
print('Total positive work done by friction', Wf[-1])
print('Total positive work done by gravity + friction', Wg[-1] + Wf[-1])
print('Work done by leg to both decel and accel the mass: ', Wl_abs[-1])

fig, axes = plt.subplots(7, 1, sharex=True)

axes[0].plot(res.t, calc_gravity_force(res.t))
axes[0].plot(res.t, calc_leg_force(res.t))
axes[0].plot(res.t, calc_friction_force(res.t, v))
axes[0].legend([r'$F_g$', r'$F_l$', r'$F_f$'])
axes[0].set_ylabel('Force')

axes[1].plot(res.t, y)
axes[1].set_ylabel(r'$y$')

axes[2].plot(res.t, v)
axes[2].set_ylabel(r'$v$')

axes[3].plot(res.t, gravity_power)
axes[3].plot(res.t, leg_power)
axes[3].plot(res.t, friction_power)
axes[3].plot(res.t, leg_power + gravity_power + friction_power)
axes[3].set_ylabel('Power')
axes[3].legend([r'$P_g$', r'$P_l$', r'$P_f$', r'$P$'])

axes[4].plot(res.t, Wg)
axes[4].plot(res.t, Wl)
axes[4].plot(res.t, Wf)
axes[4].plot(res.t, W)
axes[4].plot(res.t, Wl_abs)
axes[4].legend([r'$W_g$', r'$W_l$', r'$W_f$', r'$W$', r'$W_l^*$'])
axes[4].set_ylabel('Work')

axes[5].plot(res.t, potential_energy)
axes[5].plot(res.t, kinetic_energy)
axes[5].plot(res.t, total_energy)
axes[5].legend([r'$E_p$', r'$E_k$', r'$E$'])
axes[5].set_ylabel('Energy')

# this graph shows that the total work done to the mass equals the kinetic
# energy at any time
axes[6].plot(res.t, W, res.t, kinetic_energy)
axes[6].legend([r'$W$', r'$E_k$'])

axes[-1].set_xlabel('Time')

for ax in axes:
    ax.grid()

plt.show()