###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, G.F. Forsyth.

# Coding Assignment: Rocket

The equations of motion for a rocket in purely vertical flight are given by

\begin{align}
\frac{dh}{dt} &= v\\
(m_s+m_p) \frac{dv}{dt}& = -(m_s+m_p)g + \dot{m}_pv_e - \frac{1}{2}\rho v|v|AC_D
\end{align}

$h$ is the altitude of the rocket

$m_s = 50kg$ is the weight of the rocket shell

$g = 9.81 \frac{m}{s^2}$

$\rho = 1.091 \frac{kg}{m^3}$ is the average air density (assumed constant throughout flight)

$A = \pi r^2$ is the maximum cross sectional area of the rocket, where $r = 0.5 m$

$v_e = 325 \frac{m}{s}$ is the exhaust speed

$C_D = 0.15 $ is the drag coefficient

$m_{po} = 100 kg$ at time $t = 0$ is the initial weight of the rocket propellant

The mass of the remaining propellant is given by:

$$m_p = m_{po} - \int^t_0 \dot{m}_p d\tau$$

where $\dot{m}_p$ is the time-varying burn rate given by the following figure:

Propellant Burn Rate



Using Euler's method with a timestep of $\Delta t=0.1s$, create a Python script to calculate the altitude and velocity of the rocket from launch until crash down.  






## Assessment:

To check your answers, you can register for [MAE 6286: Practical Numerical Methods with Python](http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about).

1. At time $t=3.2s$, what is the mass (in kg) of rocket propellant remaining in the rocket?

2. What is the maximum speed of the rocket in $\frac{m}{s}$?
    At what time does this occur (in seconds)? 
    What is the altitude at this time (in meters)? 
    
3. What is the rocket's maximum altitude during flight (in meters)? At what time (in seconds) does this occur?

4. At what time (in seconds) does the rocket impact the ground? What is the velocity of the rocket (in $\frac{m}{s}$) at time of impact?

## Derivation of the rocket equations

In case you are kind of confused about the rocket equations, here we show how to get to them. 

We will start from the widely used [Reynolds' transport theorem](http://en.wikipedia.org/wiki/Reynolds_transport_theorem). For an extensive property $N$ in a volume $\Omega$, with corresponding intensive property $\eta = N/m$, the Reynolds' transport theorem reads

$$
\frac{DN}{Dt} =  \frac{\partial}{\partial t} \int_\Omega \eta \rho {\rm dV}  + \oint_{\partial \Omega} \eta \rho \mathbf{V} \cdot \mathbf{n} {\rm d} S
$$

where $m$ is mass, dV is an element of volume, $\mathbf{V}$ the velocity, $\mathbf{n}$ the unit normal vector pointing out of the control volume and $\rho$ density. We will use a control volume that encloses the rocket, and moves with it. If the velocity of the control volume is $\mathbf{V}_{\rm CV}$, we can write the velocity as $\mathbf{V} = \mathbf{V}_{\rm CV} + \mathbf{V}^\prime$, and the Reynolds' transport theorem becomes

$$
\frac{DN}{Dt} =  \frac{\partial}{\partial t} \int_\Omega \eta \rho {\rm dV}  + \oint_{\partial \Omega} \eta \rho \mathbf{V}_{\rm CV} \cdot \mathbf{n} {\rm d} S + \oint_{\partial \Omega} \eta \rho \mathbf{V}^\prime \cdot \mathbf{n} {\rm d} S.
$$

The second term in the right hand side of this equation can be rewritten using the divergence theorem

$$
\oint_{\partial \Omega} \eta \rho \mathbf{V}_{\rm CV} \cdot \mathbf{n} {\rm d} S = \int_{\Omega} \eta \rho \nabla \cdot \mathbf{V}_{\rm CV} {\rm dV},
$$

but for a non-deforming control volume, $\nabla \cdot \mathbf{V}_{\rm CV} = 0$, yielding

$$
\frac{DN}{Dt} =  \frac{\partial}{\partial t} \int_\Omega \eta \rho {\rm dV}  + \oint_{\partial \Omega} \eta \rho \mathbf{V}^\prime \cdot \mathbf{n} {\rm d} S.
$$  

Replacing $N = m \mathbf{V}^\prime$, we get the momentum conservation equation with a moving non-deforming control volume:

$$
\frac{D(m\mathbf{V}^\prime)}{Dt} =  \frac{\partial}{\partial t} \int_\Omega \mathbf{V}^\prime \rho {\rm dV}  + \oint_{\partial \Omega} \mathbf{V}^\prime \rho (\mathbf{V}^\prime \cdot \mathbf{n}) {\rm d} S.
$$  

For non-accelerating control volume, we can apply Newton's second law to compute the forces: 

$$
\frac{D(m\mathbf{V}^\prime)}{Dt} = \frac{D(m\mathbf{V})}{Dt} =  \frac{D\mathbf{p}}{Dt} = \sum \mathbf{F}.
$$ 

However, in this case the control volume is moving with the accelerating rocket, and we need to be a bit more careful. We can write $\frac{D(m\mathbf{V}^\prime)}{Dt} = \frac{D(m\mathbf{V})}{Dt} - \frac{D(m\mathbf{V}_{\rm CV})}{Dt}$, and considering the control volume is not deforming, $\frac{D(m\mathbf{V}_{\rm CV})}{Dt} = m\mathbf{a}_{\rm CV}$, where $\mathbf{a}_{\rm CV}$ is the acceleration of the control volume. On the other hand, $\mathbf{V}$ is taken from an inertial (non-accelerating) frame of reference, and Newton's second law is valid there, then

$$
\frac{D(m\mathbf{V}^\prime)}{Dt} = \frac{D(m\mathbf{V})}{Dt} - \frac{D(m\mathbf{V}_{\rm CV})}{Dt} = \sum \mathbf{F} - m \mathbf{a}_{\rm CV}.
$$

This way we get a momentum equation for an accelerating control volume:

$$
\sum \mathbf{F} - m \mathbf{a}_{\rm CV} =  \frac{\partial}{\partial t} \int_\Omega \mathbf{V}^\prime \rho {\rm dV}  + \oint_{\partial \Omega} \mathbf{V}^\prime \rho (\mathbf{V}^\prime \cdot \mathbf{n}) {\rm d} S.
$$

Let's examine this last equation to make it specific for the rocket problem. First, the rocket is subject to two forces: gravity and drag. This makes the sum of forces:

$$
\sum \mathbf{F} = (m_c + m_p) \mathbf{g} - \frac{1}{2} \rho_a \mathbf{v} |\mathbf{v}| A  C_D.
$$

where $\mathbf{v}$ is the velocity of the rocket, equal to the velocity of the control volume $\mathbf{V}_{\rm CV}$, and $\rho_a$ is the density of air.

Now, let's deal with the volume integral. Check out Figure 1: the control volume encloses the whole rocket, and the border is slightly away from the nozzle. We can make a cut on the control volume to analyze it in two parts: the top part consisting of the container only, and the bottom part from the nozzle down. The top part is a large container that moves with the control volume, and the relative velocities of the propellant inside are negligible. In the bottom part, the fluid is moving quite fast, but the amount of fluid there is constant in time. These two considerations allow us to drop the volume integral.



#### Figure 1. Control volume of the rocket.

Done with the volume integral, now the surface integral. If we consider that the burnt propellant is coming out of the control volume with a constant velocity profile, we can rewrite the surface integral as

$$
\oint_{\partial \Omega} \mathbf{V}^\prime \rho_p (\mathbf{V}^\prime \cdot \mathbf{n}) {\rm d} S = \mathbf{V}^\prime \oint_{\partial \Omega} \rho_p (\mathbf{V}^\prime \cdot \mathbf{n}) {\rm d} S = \mathbf{V}^\prime \dot{m}_p = \mathbf{v}_e \dot{m}_p,
$$

where $\mathbf{v}_e$ is the exhaust velocity of the burnt propellant coming out of the nozzle with respect to the control volume, $\rho_p$ is the density of the propellant, and $\dot{m}_p$ is the mass flow rate of propellant. By saying that $\rho = \rho_p$ we are considering that the mass of air is negligible compared to the mass of propellant coming out of the nozzle.

Finally, the momentum equation for the rocket is

$$
(m_c + m_p) \mathbf{g} - \frac{1}{2} \rho_a \mathbf{v} |\mathbf{v}| A  C_D - (m_c + m_p)\frac{{\rm d} \mathbf{v}}{{\rm d}t} =  \mathbf{v}_e \dot{m}_p.
$$

In practical terms, we're only interested in the $y$ component. Vectors $\mathbf{g}$ and $\mathbf{v}_e$ are pointing down, which gives the scalar equation:

$$
-(m_c + m_p) g - \frac{1}{2} \rho_a v^2 A  C_D - (m_c + m_p)\frac{{\rm d} v}{{\rm d}t} =  -v_e \dot{m}_p.
$$

This is the equation you will work with for your assignment!

---

###### The cell below loads the style of the notebook.


```python
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
```




<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>
<style>

@font-face {
    font-family: "Computer Modern";
    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

#notebook_panel { /* main background */
    background: rgb(245,245,245);
}

div.cell { /* set cell width */
    width: 750px;
}

div #notebook { /* centre the content */
    background: #fff; /* white background for content */
    width: 1000px;
    margin: auto;
    padding-left: 0em;
}

#notebook li { /* More space between bullet points */
    margin-top:0.8em;
}

/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running { 
    border: 1px solid #111;
}

/* Put a solid color box around each cell and its output, visually linking them*/
div.cell.code_cell {
    background-color: rgb(256,256,256); 
    border-radius: 0px; 
    padding: 0.5em;
    margin-left:1em;
    margin-top: 1em;
}

div.text_cell_render{
    font-family: 'Alegreya Sans' sans-serif;
    line-height: 140%;
    font-size: 125%;
    font-weight: 400;
    width:600px;
    margin-left:auto;
    margin-right:auto;
}


/* Formatting for header cells */
.text_cell_render h1 {
    font-family: 'Nixie One', serif;
    font-style:regular;
    font-weight: 400;    
    font-size: 45pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.5em;
    margin-top: 0.5em;
    display: block;
}

.text_cell_render h2 {
    font-family: 'Nixie One', serif;
    font-weight: 400;
    font-size: 30pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    display: block;
}	

.text_cell_render h3 {
    font-family: 'Nixie One', serif;
    margin-top:16px;
    font-size: 22pt;
    font-weight: 600;
    margin-bottom: 3px;
    font-style: regular;
    color: rgb(102,102,0);
}

.text_cell_render h4 {    /*Use this for captions*/
    font-family: 'Nixie One', serif;
    font-size: 14pt;
    text-align: center;
    margin-top: 0em;
    margin-bottom: 2em;
    font-style: regular;
}

.text_cell_render h5 {  /*Use this for small titles*/
    font-family: 'Nixie One', sans-serif;
    font-weight: 400;
    font-size: 16pt;
    color: rgb(163,0,0);
    font-style: italic;
    margin-bottom: .1em;
    margin-top: 0.8em;
    display: block;
}

.text_cell_render h6 { /*use this for copyright note*/
    font-family: 'PT Mono', sans-serif;
    font-weight: 300;
    font-size: 9pt;
    line-height: 100%;
    color: grey;
    margin-bottom: 1px;
    margin-top: 1px;
}

.CodeMirror{
    font-family: "PT Mono";
    font-size: 90%;
}

</style>





370e8360d51ad64830756b43f604a5791b117fbc
ipynb
Jupyter Notebook
lessons/01_phugoid/Rocket_Assignment.ipynb
rolando-contribute/numerical-mooc
5f2115666006bf6e6367320fff46ddc1e0e32044
CC-BY-3.0
lessons/01_phugoid/Rocket_Assignment.ipynb
rolando-contribute/numerical-mooc
5f2115666006bf6e6367320fff46ddc1e0e32044
CC-BY-3.0
lessons/01_phugoid/Rocket_Assignment.ipynb
rolando-contribute/numerical-mooc
5f2115666006bf6e6367320fff46ddc1e0e32044
CC-BY-3.0
Qwen/Qwen-72B
1. YES
2. YES