In PH121 you will learn(or have already learned) how the position, velocity, and acceleration of an object are related to one another. One of the more challenging cases occurs when the acceleration of the object is not constant. The acceleration could vary with position, velocity, or time. In these cases, you may not be able to use calculus to find the position and/or velocity of an object. Instead, you’ll have to repeatedly use these equations \[ v_{i+1} = v_i + a_i \Delta t \] \[ x_{i+1} = x_i + v_i \Delta t \]
with a very small \(\Delta t\). If \(\Delta t\) is sufficiently small, assuming that the acceleration is constant is a reasonable approximation. Below we will explain the physics of a few important situations where Euler’s method is a good choice to solve the problem.
The drag force (due to air resistance) is given by:
\[ F_\text{drag} = {1\over 2} C \cdot \rho_\mathrm{air} \cdot A \cdot v^2 \]
where \(C\) is the drag coefficient, \(\rho_\mathrm{air}\) is the density of air, and \(A\) is the cross-sectional area of the projectile. For an object that is dropped (or thrown) vertically, the acceleration is given by the following
\[ \begin{align*} -m g - {1\over 2} \rho A C v^2 &= m a\\ -g - {1\over 2 m} \rho A C v^2 &= a\\ \end{align*} \]
and we can use the standard kinematic equations with a very small \(\Delta t\) to determine the position and velocity at future times. \[ y_\mathrm{n+1} = y_n + v_n \Delta t \] \[ v_{n+1} = v_n - a \Delta t \]
In the code cell below you will find an example of Euler’s method for a small styrofoam ball that is dropped from the top of a \(100\)-m-tall building. The styrofoam ball has a mass \(m = 100\) g, drag constant \(C = 0.02\), radius \(r = 0.25\) m and the density of air is \(1.29\) kg/m\(^3\).
from matplotlib import pyplot as plt
from numpy import pi
dt = 0.01
v = [0]
y =[100]
t = [0]
g = 9.8
C = 0.02
rho = 1.29
r = .25
A = pi * r**2
m = .1
while y[-1] > 0:
a = 1/2/m * rho * A * C * v[-1]**2 - g
v.append(v[-1]+ a * dt)
y.append(y[-1] + v[-1] * dt)
t.append(t[-1] + dt)
fig = plt.figure()
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
ax1.plot(t,y,'r.')
ax2.plot(t,v,'r.')
plt.show()
For two dimensional motion, we need to decompose the drag force into x and y components.
\[ \begin{align*} F_x &= F_D \cos(\theta)\\ &= F_D {v_x\over v}\\ &= -{1\over 2} \rho A C v^2 {v_x\over v}\\ &= -{1\over 2} \rho A C v v_x\\ \end{align*} \]
A similar expression for the y component of the drag force is \[ F_y = -{1\over 2} \rho A C v v_y\\ \]
The liftoff of a rocket provides some interesting physics to consider. The mass of the rocket + fuel is decreasing as the fuel burns off. Hence, the acceleration of the rocket is very much not a constant. However, in your code you don’t want to decrease the mass of the rocket-fuel system forever or the mass will become lower than the mass of the rocket alone, and eventually will go negative (yikes!). Watch out for that.
The force of the thrust on the rocket is given by: \[ F = v_e {dM\over dt}\] where \({dM\over dt}\) is the rate of mass loss due to the fuel burning and \(v_e\) is the speed with which the exhaust exits the rocket. Can you do a unit analysis on this function and convince yourself that it has the right units?
The force of air drag on the rocket is quadratic. In other words, it is given by: \[F_\text{drag} = {1\over 2} \rho A C v^2\] where \(\rho\) is the density of air, \(A\) is the cross-sectional area of the rocket, and \(C\) is the drag coefficient. Also notice that the direction of the drag force changes at the peak of the trajectory. For the upward journey, the drag force points downward, and for the downward part of the journey, the drag force points upward. To handle this efficiently when writing your code, I recommend calculating this force like this: \[ F_D = {1\over 2} \rho A C v |v|\] so that the sign of this calculation changes naturally at the peak of the trajectory.
Since the rocket may travel a considerable distance into the atmosphere, we must take into account that the density of the atmosphere varies with altitude according to: \[\rho(h) = (1.09 - {0.0065 h\over 300})^{2.5}\]
\[ {1\over 2 m} \rho A C = {B \over m} = .0039 + {.0058\over 1 + e^{v - 35\over 5}} \]