Simulate a semi-realistic baseball trajectory, involving air resistance and spin.
import numpy as np
import ode
%matplotlib inline
import matplotlib.pyplot as plt
plt.matplotlib.style.use('ggplot')
Occurs at high Reynolds numbers, i.e., turbulent flow. Only approximate:
$$ \mathbf{F}_2 = -b_2 v^2 \mathbf{\hat v} = -b_2 v \mathbf{v} $$Relate $b_2$ to the medium that is displaced:
$$ b_2 = \frac{1}{2} C_D \rho A $$where $C_D$ is called the drag coefficient.
Magnus effect: airflow is changed around a spinning object. The Magnus force is
$$ \mathbf{F}_M = \alpha \boldsymbol{\omega} \times \mathbf{v} $$where $\boldsymbol{\omega}$ is the ball's angular velocity in rad/s (e.g., 200/s for a baseball).
For a sphere the proportionality constant $\alpha$ can be written
$$ \mathbf{F}_M = \frac{1}{2} C_L \rho A \frac{v}{\omega} \boldsymbol{\omega} \times \mathbf{v} $$where $C_L$ is the lift coefficient, $\rho$ the air density, $A$ the ball's cross section. (Advantage of defining $C_L$ this way: when spin and velocity are perpendicular, the Magnus force is simply $F_M = \frac{1}{2} C_L \rho A v^2$.)
$C_L$ is mainly a function of the spin parameter
$$ S = \frac{r\omega}{v} $$with the radius $r$ of the ball. In general we write
$$ \mathbf{F}_M = \frac{1}{2} C_L \frac{\rho A r}{S} \boldsymbol{\omega} \times \mathbf{v} $$For a baseball, experimental data show approximately a power law dependence of $C_L$ on $S$
$$ C_L = 0.62 \times S^{0.7} $$All together:
\begin{align} \mathbf{F}_M &= \alpha\ \boldsymbol{\omega} \times \mathbf{v}\\ v &= \sqrt{\mathbf{v}\cdot\mathbf{v}}\\ S &= \frac{r\omega}{v}\\ C_L &= 0.62 \times S^{0.7}\\ \alpha &= \frac{1}{2} C_L \frac{\rho A r}{S} \end{align}(quadratic drag $-\frac{b_2}{m} v \mathbf{v}$ included.)
Implement the full baseball equations of motions:
For the cross product you can look at numpy.cross().
We will live-code the baseball simulation in class (i.e., build it from scratch), but if you want to work on this problem on your own and need some starter code, see below.
Full solution will be posted as baseball_solution.ipynb
.
Simulate baseball throw for initial velocity $\mathbf{v} = (30\,\text{m/s}, 0)$.
Plot x vs y and x vs z (to see curving).
Try out different spins; a good value is $\boldsymbol{\omega} = 200\,\text{rad/s} \times (0, 1, 1)$.
Simulate the baseball throw with
Plot the three scenarios in 2D planes: x-y (side view) and x-z (top view).
Use simple matplotlib
3D plot. (BONUS: Make it work with vpython)
If we use the ipympl
backend for matplotlib then we will be able to interactively rotate our matplotlib 3D graphics. (Note: If this does not seem to work, disable adblockers and allow javascript on the page.)
%matplotlib ipympl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(r[:,1], r[:,3], r[:,2], 'o', label="no spin")
# ...
# hand of the catcher, 0.2m above homeplate
ax.plot([18.4, 18.4], [0, 0], [0, 0.2], color="black", lw=6)
ax.set_xlabel("$x$ (m)")
ax.set_ylabel("$z$ (m)")
ax.set_zlabel("$y$ (m)")
ax.legend(loc="upper left", numpoints=1)
ax.figure.tight_layout()
(from Wolfram Alpha)
rho_air = 1.275 # kg/m^3
mu_air = 1.845e-5 # Pa s
L = 0.05 # m
v = 200 # m/s
def ReynoldsNumber(v, L, rho=rho_air, mu=mu_air):
return rho*v*L/mu
ReynoldsNumber(v, L)
691056.9105691056
This means that we really should have been using quadratic air resistance for the projectile that we simulated with linear air resistance... but we did this because (1) it was easy and (2) we will come back to the problem later when we look at root finding, and we can actually solve the trajectory analytically.