Lecture22-T2-blank.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 22 Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# First import some libraries\n",
    "#%matplotlib notebook \n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import scipy\n",
    "import scipy.integrate\n",
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This cell just tries to make graphs look nicer\n",
    "try:\n",
    "    import seaborn as sns\n",
    "except ImportError:\n",
    "    # This block will be run if there's an ImportError, i.e you don't have seaborn installed.\n",
    "    sns = False\n",
    "    print (\"If you want to try different figure formatting, \"\n",
    "           \"type 'conda install seaborn' at an anaconda command prompt or terminal. \"\n",
    "           \"See https://stanford.edu/~mwaskom/software/seaborn/ for details\")\n",
    "    # If not using seaborn, we can still control the size of the figures this way\n",
    "    from pylab import rcParams\n",
    "    rcParams['figure.figsize'] = 6, 4\n",
    "else:\n",
    "    # This block will be run if there is no ImportError\n",
    "    sns.set_style(\"ticks\")\n",
    "    sns.set_context(\"poster\",rc={\"figure.figsize\": (6, 4)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In terms of number of moles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Reactor volumes\n",
    "Vliquid = 4000. # dm3\n",
    "Vgas = 5000. # dm3\n",
    "\n",
    "# Molar gas constants\n",
    "R_JmolK = 8.3145 # J/mol/K\n",
    "R_atmLmolK = 8.206e-2 # Atm.L /mol/K\n",
    "\n",
    "def dvector_dt(vector, t):\n",
    "    \"\"\"\n",
    "    Time derivative of the vector containing Na, Nb, Nc, Nd, Ns, T \n",
    "    \"\"\"\n",
    "    Na, Nb, Nc, Nd, Ns, T = vector\n",
    "\n",
    "    \"\"\"\n",
    "    Reaction 1.\n",
    "    A + B --> C + 1/2 D \n",
    "    \"\"\"\n",
    "\n",
    "    DHRx1 = -45400. # J/mol\n",
    "    A1 = 4e14 #per hour\n",
    "    E1 = 128e3 #J/mol\n",
    "    k1 = ###\n",
    "    r1 = ###\n",
    "    \n",
    "    # Rates of generation (mol/hr)\n",
    "    raV = ###\n",
    "    rbV = ###\n",
    "    rcV = ###\n",
    "    rdV = ###\n",
    "    rsV = ###\n",
    "    \n",
    "    # Pressure buildup and venting equations for gas headspace.\n",
    "    P = ###  \"D\" is the only gas, in volume Vgas\n",
    "    \n",
    "    # Normal operation: pressure control valve releases all the gas (D) being made\n",
    "    Fvent = rdV\n",
    "    # Plus a bit more or less as needed to restore pressure to 4.4atm\n",
    "    Fvent += (P - 4.4)*1e5\n",
    "    if Fvent > 11400:\n",
    "        # Pressure control valve not big enough.\n",
    "        # Flow choked through 1\" pipe\n",
    "        Fvent = (P-1.)*3360\n",
    "    if P > 28.2:\n",
    "        # Safety rupture disc breaks.\n",
    "        # Additional pressure relief through 4\" pipe\n",
    "        Fvent += (P-1.)*15360\n",
    "\n",
    "    # Mole balances\n",
    "    dNa_dt = ### \n",
    "    dNb_dt = ### \n",
    "    dNc_dt = ### \n",
    "    dNd_dt = ###\n",
    "    dNs_dt = ###\n",
    "    \n",
    "    # Total heat capacity\n",
    "    # Sum of product of cp and mass of liquid-phase components in reactor\n",
    "    CPtot = 1.26e7 # J/K\n",
    "    \n",
    "    UA = 2.77e6 # J/hr/K\n",
    "    Ta = 373.15 # K\n",
    "    dT_dt = ###\n",
    "    \n",
    "    return dNa_dt, dNb_dt, dNc_dt, dNd_dt, dNs_dt, dT_dt\n",
    "\n",
    "\n",
    "# Initial conditions\n",
    "T0 = 422 # K \n",
    "P0 = 4.4 # initial pressure in the reactor, atm absolute\n",
    "\n",
    "Na0 = ###\n",
    "Nb0 = ###\n",
    "Nc0 = ###\n",
    "Nd0 = ###\n",
    "Ns0 = ###\n",
    "initial_vector = [Na0, Nb0, Nc0, Nd0, Ns0, T0]\n",
    "times = np.linspace(0,4,1000) # hours\n",
    "result = scipy.integrate.odeint(dvector_dt, initial_vector, times)\n",
    "\n",
    "Na, Nb, Nc, Nd, Ns, T = result.T # transpose and split into columns\n",
    "P = ###\n",
    "\n",
    "plt.plot(times, Na/1e3, label='$N_A$')\n",
    "plt.plot(times, Nb/1e3, label='$N_B$')\n",
    "plt.plot(times, Nc/1e3, label='$N_C$')\n",
    "plt.plot(times, Nd/1e3, label='$N_D$')\n",
    "plt.plot(times, Ns/1e3, label='$N_S$')\n",
    "plt.ylim(0,40)\n",
    "plt.legend() \n",
    "plt.xlabel('Time (h)')\n",
    "plt.ylabel('Amount (kmol)')\n",
    "plt.show()\n",
    "plt.plot(times, T)\n",
    "plt.xlabel('Time (h)')\n",
    "plt.ylabel('Temperature (K)')\n",
    "plt.ylim(400,600)\n",
    "plt.show()\n",
    "plt.plot(times, P)\n",
    "plt.xlabel('Time (h)')\n",
    "plt.ylabel('Pressure (atm)')\n",
    "plt.ylim(0,45)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## For when you need it\n",
    "if False:\n",
    "    \"\"\"\n",
    "    Reaction 2.\n",
    "    S --> 3D + stuff\n",
    "    \"\"\"\n",
    "    DHRx2 = -3.2E5 # J/mol\n",
    "    A2 = 1e84 #per hour\n",
    "    E2 = 800e3 #J/mol\n",
    "    k2 = A2*np.exp(-E2/(R_JmolK*T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "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
}