Initial Value Problems - Project#
You are going to end this module with a bang by looking at the flight path of a firework. Shown above is the initial condition of a firework, the Freedom Flyer in (a), its final height where it detonates in (b), the applied forces in the Free Body Diagram (FBD) in ©, and the momentum of the firework \(m\mathbf{v}\) and the propellent \(dm \mathbf{u}\) in (d).
The resulting equation of motion is that the acceleration is proportional to the speed of the propellent and the mass rate change \(\frac{dm}{dt}\) as such
If you assume that the acceleration and the propellent momentum are much greater than the forces of gravity and drag, then the equation is simplified to the conservation of momentum. A further simplification is that the speed of the propellant is constant, \(u=constant\), then the equation can be integrated to obtain an analytical rocket equation solution of Tsiolkovsky [1,2],
where \(m_f\) and \(m_0\) are the mass at beginning and end of flight, \(u\) is the speed of the propellent, and \(\Delta v=v_{final}-v_{initial}\) is the change in speed of the rocket from beginning to end of flight. Equation 2.b only relates the final velocity to the change in mass and propellent speed. When you integrate Eqn 2.a, you will have to compare the velocity as a function of mass loss.
Your first objective is to integrate a numerical model that converges to equation (2.b), the Tsiolkovsky equation. Next, you will add drag and gravity and compare the results between equations (1) and (2). Finally, you will vary the mass change rate to achieve the desired detonation height.
1. Create a simplerocket function that returns the velocity, \(v\),
the acceleration, \(a\), and the mass rate change \(\frac{dm}{dt}\), as a
function of the \(state = [position,~velocity,~mass] = [y,~v,~m]\) using
eqn (2.a). Where the mass rate change \(\frac{dm}{dt}\) and the propellent
speed \(u\) are constants. The average velocity of gun powder propellent
used in firework rockets is \(u=250\) m/s [3,4].
\(\frac{d~state}{dt} = f(state)\)
\(\left[\begin{array}{c} v\\a\\ \frac{dm}{dt} \end{array}\right] = \left[\begin{array}{c} v\\ \frac{u}{m}\frac{dm}{dt} \\ \frac{dm}{dt} \end{array}\right]\)
Use an integration method to
integrate the simplerocket function. Demonstrate that your solution
converges to equation (2.b) the Tsiolkovsky equation. Use an initial
state of y=0 m, v=0 m/s, and m=0.25 kg.
Integrate the function until mass, \(m_{f}=0.05~kg\), using a mass rate change of \(\frac{dm}{dt}=0.05\) kg/s.
Hint: your integrated solution will have a current mass that you can use to create \(\frac{m_{f}}{m_{0}}\) by dividing state[2]/(initial mass), then your plot of velocity(t) vs mass(t)/mass(0) should match Tsiolkovsky’s
\(\log\left(\frac{m_{f}}{m_{0}}\right) = \log\left(\frac{state[2]}{0.25~kg}\right) = \frac{state[1]}{250~m/s} = \frac{-\Delta v+error}{u}\) where \(error\) is the difference between your integrated state variable and the Tsiolkovsky analytical value.
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
def simplerocket(state,dmdt=0.05, u=250):
'''Computes the right-hand side of the differential equation
for the acceleration of a rocket, without drag or gravity, in SI units.
Arguments
----------
state : array of three dependent variables [y v m]^T
dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s
u : speed of propellent expelled (default is 250 m/s)
Returns
-------
derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T
'''
dstate = np.zeros(np.shape(state))
# your work
return dstate
m0=0.25
mf=0.05
dm=0.05
t = np.linspace(0,(m0-mf)/dm,500)
dt=t[1]-t[0]
2. You should have a converged solution for integrating simplerocket. Now, create a more relastic function, rocket that incorporates gravity and drag and returns the velocity, \(v\), the acceleration, \(a\), and the mass rate change \(\frac{dm}{dt}\), as a function of the \(state = [position,~velocity,~mass] = [y,~v,~m]\) using eqn (1). Where the mass rate change \(\frac{dm}{dt}\) and the propellent speed \(u\) are constants. The average velocity of gun powder propellent used in firework rockets is \(u=250\) m/s [3,4].
\(\frac{d~state}{dt} = f(state)\)
\(\left[\begin{array}{c} v\\a\\ \frac{dm}{dt} \end{array}\right] = \left[\begin{array}{c} v\\ \frac{u}{m}\frac{dm}{dt}-g-\frac{c}{m}v^2 \\ \frac{dm}{dt} \end{array}\right]\)
Use two integration methods to integrate the rocket function, one explicit method and one implicit method. Demonstrate that the solutions converge to equation (2.b) the Tsiolkovsky equation. Use an initial state of y=0 m, v=0 m/s, and m=0.25 kg.
Integrate the function until mass, \(m_{f}=0.05~kg\), using a mass rate change of \(\frac{dm}{dt}=0.05\) kg/s, .
Compare solutions between the simplerocket and rocket integration, what is the height reached when the mass reaches \(m_{f} = 0.05~kg?\)
def rocket(state,dmdt=0.05, u=250,c=0.18e-3):
'''Computes the right-hand side of the differential equation
for the acceleration of a rocket, with drag, in SI units.
Arguments
----------
state : array of three dependent variables [y v m]^T
dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s
u : speed of propellent expelled (default is 250 m/s)
c : drag constant for a rocket set to 0.18e-3 kg/m
Returns
-------
derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T
'''
g=9.81
dstate = np.zeros(np.shape(state))
# your work
return dstate
3. Solve for the mass change rate that results in detonation at a height of 300 meters. Create a function f_dm that returns the final height of the firework when it reaches \(m_{f}=0.05~kg\). The inputs should be
\(f_{m}= f_{m}(\frac{dm}{dt},~parameters)\)
where \(\frac{dm}{dt}\) is the variable you are using to find a root and \(parameters\) are the known values, m0=0.25, c=0.18e-3, u=250. When \(f_{m}(\frac{dm}{dt}) = 0\), you have found the correct root.
Plot the height as a function of time and use a star to denote detonation at the correct height with a '*'-marker
Approach the solution in two steps, use the incremental search
incsearch with 5-10
sub-intervals limit the number of times you call the
function. Then, use the modified secant method to find the true root of
the function.
a. Use the incremental search to find the two closest mass change rates within the interval \(\frac{dm}{dt}=0.05-0.4~kg/s.\)
b. Use the modified secant method to find the root of the function \(f_{m}\).
c. Plot your solution for the height as a function of time and indicate the detonation with a *-marker.
def f_dm(dmdt, m0 = 0.25, c = 0.18e-3, u = 250):
''' define a function f_dm(dmdt) that returns
height_desired-height_predicted[-1]
here, the time span is based upon the value of dmdt
arguments:
---------
dmdt: the unknown mass change rate
m0: the known initial mass
c: the known drag in kg/m
u: the known speed of the propellent
returns:
--------
error: the difference between height_desired and height_predicted[-1]
when f_dm(dmdt) = 0, the correct mass change rate was chosen
'''
# your work
return error
def mod_secant(func,dx,x0,es=0.0001,maxit=50):
'''mod_secant: Modified secant root location zeroes
root,[fx,ea,iter]=mod_secant(func,dfunc,xr,es,maxit,p1,p2,...):
uses modified secant method to find the root of func
arguments:
----------
func = name of function
dx = perturbation fraction
xr = initial guess
es = desired relative error (default = 0.0001 )
maxit = maximum allowable iterations (default = 50)
p1,p2,... = additional parameters used by function
returns:
--------
root = real root
fx = func evaluated at root
ea = approximate relative error ( )
iter = number of iterations'''
iter = 0;
xr=x0
for iter in range(0,maxit):
xrold = xr;
dfunc=(func(xr+dx)-func(xr))/dx;
xr = xr - func(xr)/dfunc;
if xr != 0:
ea = abs((xr - xrold)/xr) * 100;
else:
ea = abs((xr - xrold)/1) * 100;
if ea <= es:
break
return xr,[func(xr),ea,iter]
References#
Math 24 Rocket Motion. https://www.math24.net/rocket-motion/\
Kasdin and Paley. Engineering Dynamics. ch 6-Linear Momentum of a Multiparticle System pp234-235 Princeton University Press
https://www.apogeerockets.com/Rocket_Motors/Estes_Motors/13mm_Motors/Estes_13mm_1_4A3-3T