CompMech03-IVPs_project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initial Value Problems - Project\n",
    "\n",
    "![Initial condition of firework with FBD and sum of momentum](../images/firework.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We 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 (c), and the __momentum__ of the firework $m\\mathbf{v}$ and the propellent $dm \\mathbf{u}$ in (d). \n",
    "\n",
    "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\n",
    "\n",
    "$$\\begin{equation}\n",
    "m\\frac{dv}{dt} = u\\frac{dm}{dt} -mg - cv^2.~~~~~~~~(1)\n",
    "\\end{equation}$$\n",
    "\n",
    "If we 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](https://www.math24.net/rocket-motion/) [1,2], \n",
    "\n",
    "$$\\begin{equation}\n",
    "m\\frac{dv}{dt} = u\\frac{dm}{dt}~~~~~(2.a)\n",
    "\\end{equation}$$\n",
    "\n",
    "$$\\begin{equation}\n",
    "\\frac{m_{f}}{m_{0}}=e^{-\\Delta v / u},~~~~~(2.b) \n",
    "\\end{equation}$$\n",
    "\n",
    "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. \n",
    "\n",
    "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. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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]. \n",
    "\n",
    "$\\frac{d~state}{dt} = f(state)$\n",
    "\n",
    "$\\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]$\n",
    "\n",
    "Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) to integrate the `simplerocket` 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. \n",
    "\n",
    "Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s. \n",
    "\n",
    "_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_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams.update({'font.size': 22})\n",
    "plt.rcParams['lines.linewidth'] = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simplerocket(state,dmdt=0.05, u=250):\n",
    "    '''Computes the right-hand side of the differential equation\n",
    "    for the acceleration of a rocket, without drag or gravity, in SI units.\n",
    "    \n",
    "    Arguments\n",
    "    ----------    \n",
    "    state : array of three dependent variables [y v m]^T\n",
    "    dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
    "    u    : speed of propellent expelled (default is 250 m/s)\n",
    "    Returns\n",
    "   -------\n",
    "   derivs: array of three derivatives [v (u/m*dmdt) -dmdt]^T\n",
    "   '''\n",
    "    \n",
    "    \n",
    "    dstate = np.zeros(np.shape(state))\n",
    "    return dstate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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]. \n",
    "\n",
    "$\\frac{d~state}{dt} = f(state)$\n",
    "\n",
    "$\\left[\\begin{array}{c} v\\\\a\\\\ \\frac{dm}{dt} \\end{array}\\right] = \n",
    "\\left[\\begin{array}{c} v\\\\ \\frac{u}{m}\\frac{dm}{dt}-g-\\frac{c}{m}v^2 \\\\ \\frac{dm}{dt} \\end{array}\\right]$\n",
    "\n",
    "Use [two integration methods](../notebooks/03_Get_Oscillations.ipynb) 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. \n",
    "\n",
    "Integrate the function until mass, $m_{f}=0.05~kg$, using a mass rate change of $\\frac{dm}{dt}=0.05$ kg/s, . \n",
    "\n",
    "Compare solutions between the `simplerocket` and `rocket` integration, what is the height reached when the mass reaches $m_{f} = 0.05~kg?$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n",
    "g=9.81\n",
    "derivs = np.array([state[1], (u/state[2])dmdt-g-(c/state[2])((state[1])**2), -dmdt])\n",
    "return derivs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rocket(state,dmdt=0.05, u=250,c=0.18e-3):\n",
    "    '''Computes the right-hand side of the differential equation\n",
    "    for the acceleration of a rocket, with drag, in SI units.\n",
    "    \n",
    "    Arguments\n",
    "    ----------    \n",
    "    state : array of three dependent variables [y v m]^T\n",
    "    dmdt : mass rate change of rocket in kilograms/s default set to 0.05 kg/s\n",
    "    u    : speed of propellent expelled (default is 250 m/s)\n",
    "    c : drag constant for a rocket set to 0.18e-3 kg/m\n",
    "    Returns\n",
    "    -------\n",
    "    derivs: array of three derivatives [v (u/m*dmdt-g-c/mv^2) -dmdt]^T\n",
    "    '''\n",
    "    dstate = np.zeros(np.shape(state))\n",
    "    return dstate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 \n",
    "\n",
    "$f_{m}= f_{m}(\\frac{dm}{dt},~parameters)$\n",
    "\n",
    "where $\\frac{dm}{dt}$ is the variable we 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$, we have found the correct root. \n",
    "\n",
    "Plot the height as a function of time and use a star to denote detonation at the correct height with a `'*'`-marker\n",
    "\n",
    "Approach the solution in two steps, use the incremental search [`incsearch`](../notebooks/04_Getting_to_the_root.ipynb) with 5-10 sub-intervals _we want to limit the number of times we call the function_. Then, use the modified secant method to find the true root of the function.\n",
    "\n",
    "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.$\n",
    "\n",
    "b. Use the modified secant method to find the root of the function $f_{m}$.\n",
    "\n",
    "c. Plot your solution for the height as a function of time and indicate the detonation with a `*`-marker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_m(dmdt,m0=0.25, c=0.18e-3, u=250):\n",
    "    ''' define a function f_m(dmdt) that returns \n",
    "    height_desired-height_predicted[-1]\n",
    "    here, the time span is based upon the value of dmdt\n",
    "    \n",
    "    arguments:\n",
    "    ---------\n",
    "    dmdt: the unknown mass change rate\n",
    "    m0: the known initial mass\n",
    "    c: the known drag in kg/m\n",
    "    u: the known speed of the propellent\n",
    "    \n",
    "    returns:\n",
    "    --------\n",
    "    error: the difference between height_desired and height_predicted[-1]\n",
    "        when f_m(dmdt)= 0, the correct mass change rate was chosen\n",
    "    '''\n",
    "    \n",
    "    return error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Math 24 _Rocket Motion_. <https://www.math24.net/rocket-motion/\\>\n",
    "\n",
    "2. Kasdin and Paley. _Engineering Dynamics_. [ch 6-Linear Momentum of a Multiparticle System pp234-235](https://www.jstor.org/stable/j.ctvcm4ggj.9) Princeton University Press \n",
    "\n",
    "3. <https://en.wikipedia.org/wiki/Specific_impulse>\n",
    "\n",
    "4. <https://www.apogeerockets.com/Rocket_Motors/Estes_Motors/13mm_Motors/Estes_13mm_1_4A3-3T>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}