Bird Design.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:83256e0cd70af6a01ba0b9512ac33e93db4fb20fb3877cd7a97313d2dc92346d"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab\n",
      "from __future__ import print_function"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Using matplotlib backend: TkAgg\n",
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Carbon Fiber Mechanical Properties\n",
      "\n",
      "* E = 70 GPa (Young's modulus)\n",
      "* $\\sigma_{ult}$ = 1.6 GPa (ultimate strength)\n",
      "* density = 1550 kg/$m^3$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Failure Modes\n",
      "\n",
      "### Bending\n",
      "\n",
      "$$ \\sigma = \\frac{My}{I} $$\n",
      "\n",
      "where\n",
      "* $\\sigma$ = strain,\n",
      "* M     = moment,\n",
      "* y = distance from neutral axis,\n",
      "* I = area moment of inertia.\n",
      "\n",
      "### Buckling\n",
      "\n",
      "$$ F = \\frac{\\pi^2 E I}{KL^2} $$\n",
      "\n",
      "where \n",
      "* F = maximum or critical force (vertical load on column),\n",
      "* E = modulus of elasticity,\n",
      "* I = area moment of inertia,\n",
      "* L = unsupported length of column,\n",
      "* K = column effective length factor, whose value depends on the conditions of end support of the column, as follows.\n",
      "    * For both ends pinned (hinged, free to rotate), K = 1.0.\n",
      "    * For both ends fixed, K = 0.50.\n",
      "    * For one end fixed and the other end pinned, K = 0.699....\n",
      "    * For one end fixed and the other end free to move laterally, K = 2.0.\n",
      "* KL is the effective length of the column."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Design Problem\n",
      "\n",
      "Find the design with the most excess power given the motors.\n",
      "\n",
      "## Assuumptions\n",
      "* Rods running length of vehicle\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Power Required\n",
      "\n",
      "The absolute minimum power for flapping flight is given by:\n",
      "\n",
      "$$ P_{mp} = \\frac{1.05 k^{3/4} m^{3/2} g^{3/2} S_b^{1/4} C_{Db}^{1/4}}{\\rho^{1/2}B^{3/2}} $$\n",
      "\n",
      "$$ V_{mp} = \\frac{0.807 k^{1/4}m^{1/2}g^{1/2}}{\\rho^{1/2}B^{1/2}S_b^{1/4}C_{Db}^{1/4}}$$\n",
      "\n",
      "$$f = m^{3/8}g^{1/2}B^{-23/24}S^{-1/3}\\rho^{-3/8}$$\n",
      "where\n",
      "* B : wing span\n",
      "* $C_{Db}$ : body drag coefficient\n",
      "* g : Acceleration due to gravity\n",
      "* k : Induced power factor\n",
      "* m : All-up mass\n",
      "* S_b : body frontal area"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Power (Min Power)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# vehicle\n",
      "B = 1.2 # wing span, m\n",
      "C_Db = 0.1 # body drag coeff, m\n",
      "k = 1.1 # induced power factor\n",
      "m = 0.3 # mass\n",
      "S_b = 0.003 # body frontal area\n",
      "S = B/6\n",
      "\n",
      "cf_E = 70e6\n",
      "cf_ult_stress = 1.6e6\n",
      "cf_dens = 1500\n",
      "\n",
      "# environment\n",
      "g = 9.8\n",
      "rho = 1.225\n",
      "\n",
      "\n",
      "P_mp = 1.05*k**(0.75)*m**(1.5)*g**(1.5)*\\\n",
      "    S_b**(0.25)*C_Db**(0.25)*rho**(-0.25)*B**(-1.5)\n",
      "V_mp = 0.807*k**0.25*m**0.5*g**0.5*\\\n",
      "    rho**-0.5*B**-0.5*S_b**-0.25*C_Db**-0.25\n",
      "f_mp = m**(3.0/8)*g**0.5*B**(-23.0/24)*S**(-1.0/3)*rho**(-3.0/8)\n",
      "P_mp, V_mp\n",
      "print(\"\"\"P_mp: {P_mp:f} Watts\n",
      "V_mp: {V_mp:f} m/s\n",
      "f_mp: {f_mp:f} Hz\"\"\".format(**locals()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "P_mp: 0.541040 Watts\n",
        "V_mp: 8.880888 m/s\n",
        "f_mp: 2.652103 Hz\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Power (Max Range)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "S_d = pi*B**2/4\n",
      "A = S_b*C_Db\n",
      "V_mr = k**0.25*m**0.5*g**0.5*rho**-0.5*A**-0.25*S_d**-0.25\n",
      "P_mr = k**0.75*m**1.5*g**1.5*A**0.25*rho**-0.25*S_d**-0.75\n",
      "V_mr, P_mr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "(11.689888938416754, 0.6176216636797301)"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Power Available\n",
      "\n",
      "$$ P_{servo} = T \\omega$$\n",
      "\n",
      "where\n",
      "* $P_{servo}$ : available servo power\n",
      "* T : servo torque\n",
      "* $\\omega$ : angular velocity of servo"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Manufacturer Servo/Battery Data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ozIn2Nm = 0.00706155183333\n",
      "T_servo = 300*ozIn2Nm # oz-in\n",
      "omega_servo = deg2rad(60/0.06) # rad/s\n",
      "eff_servo = 0.001 # total efficiency, setting to match Robo Raven performance\n",
      "V_servo = 7.4 # volts\n",
      "C_bat = 0.370 # battery capacity Ah\n",
      "d_bat = 20 # battery discharge rating"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 53
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Computed Servo/Battery Data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "P_servos = 2*T_servo*omega_servo # watts, note 2 servos\n",
      "A_servos = P_servos/ V_servo\n",
      "flight_time = C_bat/A_servos*60 # minutes\n",
      "A_bat = d_bat*C_bat\n",
      "\n",
      "flight_time_mp = 4.75\n",
      "A_mp = C_bat/flight_time_mp*60\n",
      "eff_servo = P_mp/V_servo/A_mp\n",
      "\n",
      "P_mp_eff = P_mp/eff_servo\n",
      "\n",
      "print(\"\"\"\n",
      "\n",
      "min power\n",
      "===========\n",
      "flight time\\t:{flight_time_mp:10.3f} min\n",
      "P\\t\\t:{P_mp_eff:10.3f} Watts\n",
      "A\\t\\t:{A_mp:10.3f} A\n",
      "\n",
      "full power\n",
      "===========\n",
      "flight time\\t:{flight_time:10.3f} min\n",
      "P\\t\\t:{P_servos:10.3f} Watts\n",
      "A\\t\\t:{A_servos:10.3f} A\n",
      "\n",
      "battery\n",
      "===========\n",
      "A bat\\t\\t:{A_bat:10.3f} A\n",
      "\"\"\"\\\n",
      ".format(**locals()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "min power\n",
        "===========\n",
        "flight time\t:     4.750 min\n",
        "P\t\t:    34.585 Watts\n",
        "A\t\t:     4.674 A\n",
        "\n",
        "full power\n",
        "===========\n",
        "flight time\t:     2.222 min\n",
        "P\t\t:    73.948 Watts\n",
        "A\t\t:     9.993 A\n",
        "\n",
        "battery\n",
        "===========\n",
        "A bat\t\t:     7.400 A\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 86
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We see that the battery is capable of supply the required discharge rate ($A_{bat} > A_{min\\_power}$)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculated efficiency, so endurance mataches RoboRaven."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "eff_servo*P_servos - P_mr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 27,
       "text": [
        "1.0092430895726943"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Maximize\n",
      "Excess Power\n",
      "\n",
      "## Design Vector\n",
      "* B - wing span\n",
      "* S - wing area\n",
      "* l - length from nose to start of tail\n",
      "* a_t - angle of tail\n",
      "* l_t - length of tail\n",
      "* w_t - width of tail\n",
      "* d_f - diameter of fuselage rod\n",
      "* t_f - thickness of fuselage rod\n",
      "* d_w - diameter of wing rod\n",
      "* t_w - thickness of wing rod\n",
      "* n_f - number of fuselage rods\n",
      "* r_f - radius of fuselage rods\n",
      "\n",
      "## Constraints\n",
      "* Payload can be moved to set c.g. location to desired.\n",
      "* Longitudinally stable.\n",
      "* Laterally stable.\n",
      "* Wing rod can't bend more than 5 degrees.\n",
      "* Wing rod can't twist more than 5 degrees.\n",
      "* Fuselage can't bend more than 5 degrees.\n",
      "* Fuselage can't twist more than 5 degrees.\n",
      "* Wing rod survives crash on wing tip.\n",
      "* Fuselage survives crash on end.\n",
      "* Fuselage survives crash on side."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.optimize\n",
      "\n",
      "def neg_excess_power(x):\n",
      "    # design vector\n",
      "    B = x[0] # wing span\n",
      "    r_f = x[1] # diameter of fuselage rod\n",
      "    t_f = x[2]\n",
      "    \n",
      "    # guess at structure relationship for now\n",
      "    m = 4*pi*r_f*t_f*cf_dens\n",
      "    S_b = B**2\n",
      "    \n",
      "    A = S_b*C_Db\n",
      "    P_mp = 1.05*k**(0.75)*m**(1.5)*g**(1.5)*\\\n",
      "        S_b**(0.25)*C_Db**(0.25)*rho**(-0.25)*B**(-1.5)\n",
      "    return -(P_servos - P_mp)\n",
      "\n",
      "def stress_bend(x):\n",
      "    # design vector\n",
      "    B = x[0] # wing span\n",
      "    r_f = x[1] # radius of fuselage rod\n",
      "    t_f = x[2]\n",
      "    \n",
      "    M = B**2/6 # moment scales by square of length\n",
      "    y = r_f\n",
      "    I = pi*r_f**3*t_f\n",
      "    return M*y/I\n",
      "    \n",
      "def stress_buckle(x):\n",
      "    # design vector\n",
      "    B = x[0] # wing span\n",
      "    r_f = x[1] # radius of fuselage rod\n",
      "    t_f = x[2]\n",
      "    \n",
      "    I = pi*r_f**3*t_f\n",
      "    K = 2 # buckling with one end free\n",
      "    F = pi**2*cf_E*I/(K*B**2)\n",
      "    A = 2*pi*r_f*t_f\n",
      "    return F/A\n",
      "    \n",
      "def con_bending(x):    \n",
      "    return -(cf_ult_stress - stress_bend(x))\n",
      "    \n",
      "def con_buckling(x):    \n",
      "    return -(stress_buckle(x) - stress_bend(x))\n",
      "\n",
      "def callback(x):\n",
      "    print(x)\n",
      "\n",
      "x0 = [1, 1, 1]\n",
      "\n",
      "constraints= [\n",
      "    #{'type':'ineq', 'fun':con_bending},\n",
      "    #{'type':'ineq', 'fun':con_buckling},\n",
      "]\n",
      "\n",
      "options = {'disp':True, 'iprint':2}\n",
      "bounds = [(0.1,10), (0.1,10), (0.1,10)]\n",
      "#bounds = None\n",
      "method = 'slsqp'\n",
      "res = scipy.optimize.minimize(neg_excess_power, x0, method=method, bounds=bounds,\n",
      "                              constraints=constraints, options = options, callback=callback)\n",
      "res"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  NIT    FC           OBJFUN            GNORM\n",
        "[ 1.  1.  1.]\n",
        "    1     5     4.786190E+07     1.122463E+08\n",
        "Inequality constraints incompatible    (Exit mode 4)\n",
        "            Current function value: 47861898.1125\n",
        "            Iterations: 1\n",
        "            Function evaluations: 5\n",
        "            Gradient evaluations: 1\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 91,
       "text": [
        "  status: 4\n",
        " success: False\n",
        "    njev: 1\n",
        "    nfev: 5\n",
        "     fun: 47861898.112543397\n",
        "       x: array([ 1.,  1.,  1.])\n",
        " message: 'Inequality constraints incompatible'\n",
        "     jac: array([-47861971.,  71792958.,  71792958.,         0.])\n",
        "     nit: 1"
       ]
      }
     ],
     "prompt_number": 91
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#References\n",
      "* http://en.wikipedia.org/wiki/Buckling\n",
      "* http://hitecrcd.com/products/servos/digital-brushless-servos/hsb-9370th-multi-purpose-brushless-titanium-gear-servo/product\n"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}