HW4-Q4-question.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 4 Question 4\n",
    "\n",
    "The reaction between Hydrogen and Bromine \n",
    "$$ H_2 + Br_2 \\rightarrow 2HBr$$\n",
    "has a complex reaction rate expression\n",
    "$$r_{rxn} = -r_{Br_2} = \\frac{k_1 C_{H_2} C_{Br_2}^{1/2}}{k_2 + C_{HBr}/C_{Br_2}}$$\n",
    "\n",
    "First we import some useful things"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "from scipy.optimize import fsolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## a)\n",
    "Define a function `rate(C_H2, C_Br2, C_HBr, k1, k2)` that takes five variables (the concentrations of $H_2$, $Br_2$, and $HBr$, and the rate constants $k_1$ and $k_2$) and returns the rate of forward reaction, $r_{rxn}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell will test your function at a few values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert rate(1.,1.,1.,1.,1.)==0.5\n",
    "assert rate(1.,1.,1.,k1=10.,k2=1.) == 5\n",
    "assert rate(10.,1.,10.,k1=10.,k2=10.) == 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the following plots will show the complex shape of the rate expression, with apparent orders with respect to one species depending on the relative concentration of another species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N=20\n",
    "h2 = np.linspace(1e-6, 1.0, N)\n",
    "br2 = np.linspace(1e-6, 1.0, N)\n",
    "X, Y = np.meshgrid(h2, br2)\n",
    "Z = rate(X, Y, 1e-8, 1., 1.)\n",
    "cs = pyplot.contourf(X, Y, Z)\n",
    "pyplot.title(\"Low $HBr$\")\n",
    "pyplot.xlabel(\"$H_2$\")\n",
    "pyplot.ylabel(\"$Br_2$\")\n",
    "cb = pyplot.colorbar()\n",
    "pyplot.show()\n",
    "br2 = np.linspace(1e-10,1e-5)\n",
    "pyplot.plot(br2,rate(1.,br2,1e-8,1.,1.))\n",
    "pyplot.title(\"Low $HBr$, High $H_2$\")\n",
    "pyplot.xlabel(\"$Br_2$\")\n",
    "pyplot.ylabel(\"Reaction rate\")\n",
    "pyplot.show()\n",
    "br2 = np.linspace(1e-10,1e-5)\n",
    "pyplot.plot(br2,rate(1.,br2,1.,1.,1.))\n",
    "pyplot.title(\"High $HBr$, High $H_2$\")\n",
    "pyplot.xlabel(\"$Br_2$\")\n",
    "pyplot.ylabel(\"Reaction rate\")\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## b) \n",
    "Define a function `rateX(X)` that returns the rate of reaction as a function of conversion $X$ of $H_2$ (i.e. $H_2$ is your limiting reactant A), assuming temperature and pressure are constant and:\n",
    "$$C_{H_2,0} = 20\\ mol/m^3$$\n",
    "$$C_{Br_2,0} = 20\\ mol/m^3$$\n",
    "$$C_{HBr,0} = 0\\ mol/m^3$$\n",
    "$$k_1 = 10.\\ mol^{-0.5} m^{1.5} s^{-1}$$\n",
    "$$k_2 = 1.0$$\n",
    "\n",
    "The function `rateX` should evaluate the variables needed and then call the function `rate`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell will test your function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert round(rateX(0.),5) == round(894.4271909999159,5)\n",
    "assert round(rateX(0.5),5) == round(105.40925533894598,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## c)\n",
    "\n",
    "For an isothermal isobaric plug flow reactor with feed volumetric flow rate $v_0=1\\ m^3/s$ and feed concentrations as given above, define a function `X_final(Vpfr)` which returns the conversion at the exit of a PFR  with the given volume `Vpfr` in m<sup>3</sup>.  It should use your function `rateX`, and `odeint` from the `scipy.integrate` module.\n",
    "\n",
    "This is essentially solving\n",
    "$$X_{final} = \\int_0^{V_{PFR}} \\frac{-r_A}{F_{A0}} dV $$\n",
    "but noting that $-r_A$ is given as a function of $X$.\n",
    "\n",
    "To avoid problems trying to integrate with limits $\\int_0^0$ you might start your function with \n",
    "```python\n",
    "if Vpfr == 0:\n",
    "    return 0\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell will test your `X_final` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert round(X_final(0.5),5) == round(0.86591397042924045,5)\n",
    "assert X_final(0) == 0\n",
    "assert X_final(1e5) > 0.999\n",
    "assert round(X_final(0.01),5) == round(0.26405829302912914,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example shows how you can specify default values for variables in a function definition, and combine them with `if` statements to control the behavior of your function:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "4\n",
      "This was called with variable2 = True!\n",
      "6\n",
      "This was called with variable2 = True!\n",
      "8\n"
     ]
    }
   ],
   "source": [
    "# variable2 is given a default value of False \n",
    "# when we define the function:\n",
    "def simple_function(variable1, variable2=False):\n",
    "    if variable2:\n",
    "        # by default this won't execute:\n",
    "        print(\"This was called with variable2 = {0}!\".format(variable2))\n",
    "    return 2*variable1\n",
    "# so if we call the function without specifying variable2,\n",
    "# it is assumed False, and simple_function just returns double the input:\n",
    "print(simple_function(1))\n",
    "print(simple_function(2))\n",
    "# but if we call it with variable2 set to True, it behaves differently:\n",
    "print(simple_function(3,True))\n",
    "# or we can set it by name:\n",
    "print(simple_function(4,variable2=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use this technique to modify your `X_final` function so that:\n",
    " * if you call it without specifying a second variable, i.e. as `X_final(V)` it will behave as before and just return the final conversion (useful if you are going to call `X_final` as part of an optimization or `fsolve` call to find V, for example).\n",
    " * but if you call this as `X_final(V, plot=True)` then it will not only return the final conversion but will also draw a plot of conversion as a function of Volume (useful for one-off uses once you know what V you want to use).\n",
    "\n",
    "Clue:\n",
    "```python\n",
    "def X_final(Vprf, plot=False):\n",
    "    ### some stuff\n",
    "    if plot:\n",
    "        pyplot.plot()\n",
    "        ### etc...\n",
    "        pyplot.show()\n",
    "    return result\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should still work as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert round(X_final(0.5),5) == round(0.86591397042924045,5)\n",
    "assert X_final(0) == 0\n",
    "assert X_final(1e5) > 0.999\n",
    "assert round(X_final(0.01),5) == round(0.26405829302912914,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## d) \n",
    "Use `fsolve` and your function `X_final` to find the volume of PFR that would give a final conversion of 0.95, and store this in a variable called `volume`.  Then use `X_final(volume, plot=True)` to plot the graph of Conversion vs Volume along the length of this reactor (i.e. from $V=0$ to whatever volume gives $X=0.95$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell will check your value of `volume`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert round(X_final(volume),5) == 0.95\n",
    "assert round(volume,1) == 2.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}