Notes_Example_RichComb.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import all the necessary toolboxes...\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.optimize import least_squares\n",
    "from scipy.interpolate import interp1d\n",
    "from numpy import*\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Rich Combustion\n",
    "\n",
    "One kmol of propane $\\left(C_3H_{8}\\right)$ is burned in air under rich conditions in an engine. Assuming rich, incomplete combustion (i.e., 30% of the propane goes to CO in the exhaust), determine (i)the equivalence ratio, and (ii)the combustion temperature.\n",
    "\n",
    "$C_3H_8 + \\phi \\left(O_2 + 3.76 N_2\\right) \\rightarrow n_1CO_2 + n_2H_2O + n_3N_2 + n_4O_2 + n_5CO + n_6H_2$\n",
    "\n",
    "where, $n_1 = 0.7*3$ and $n_5 = (1-0.7)*3$\n",
    "\n",
    "(note in this example we're using propane becuase the default chemical mechanisms don't have C8 molecules in them ... also smaller mechanisms are faster)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phi is calculated to be 1.43\n"
     ]
    }
   ],
   "source": [
    "# Set up the knowns\n",
    "fuel = 1 #kmol of C3H6\n",
    "C_fuel = 3\n",
    "H_fuel = 8\n",
    "\n",
    "# Note that we were told 70% of the octane converted to CO2 while 30% was converted to CO:\n",
    "pB = 0.7 #70% burned to CO2\n",
    "\n",
    "a = pB*(C_fuel+H_fuel/4)\n",
    "a_s = (C_fuel+H_fuel/4)\n",
    "phi = a_s/a\n",
    "print(\"Phi is calculated to be %1.2f\" % phi)\n",
    "\n",
    "n1 = pB*(C_fuel)\n",
    "n5 = (1-pB)*C_fuel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fsolve approach\n",
    "This section demonstrates how to use the \"fsolve\", the next section will cover a slightly different numerical solver appraoch.\n",
    "\n",
    "\"fsolve\" is a simple and versatile optimization algorithm that is suitable for a wide range of problems. However, it may not always converge and has limited options for optimization. On the other hand, \"least_squares\" is a more powerful optimization algorithm that is specifically designed for solving least-squares problems, but it is limited to this type of problem and may be slower than \"fsolve\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The temperature is  876.40 K\n",
      "x_CO2 = 0.1042\n",
      "x_CO = 0.0446\n"
     ]
    }
   ],
   "source": [
    "# Define our equations to solve\n",
    "def eqn(vars):\n",
    "    t, T, K, d1, a1, b1, c1, n2, n6, n3 = vars #this line sets the order of the output for the lest_squares apprach\n",
    "    eq1 = t - T/1000\n",
    "    eq2 = K - exp(2.743 - 1.761/t - 1.611/t**2 + 0.2803/t**3)\n",
    "    eq3 = d1 - 2*a_s*(1-1/phi)\n",
    "    eq4 = a1 - (1-K)\n",
    "    eq5 = b1 - (H_fuel/2 + K*C_fuel - d1*(1-K))\n",
    "    eq6 = c1 - (-C_fuel*d1*K)\n",
    "    eq7 = K - n2*n5/(n1*n6)\n",
    "    eq8 = n2 - (H_fuel/2 - d1 + n5)\n",
    "    eq9 = n3 - 3.76*a_s/phi\n",
    "    eq10 = n6 - (d1 - n5)\n",
    "    return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10]\n",
    "\n",
    "# Now call the fsolve function to solve our system of equatins\n",
    "t, T, K, d1, a1, b1, c1, n2, n6, n3 = fsolve(eqn, (1, 1000, 1, 1, 1, 1, 1, 1, 1, 1))\n",
    "\n",
    "print(\"The temperature is  %1.2f K\" % T)\n",
    "#print(\"n_CO2 = %1.2f\" % n1)\n",
    "#print(\"n_CO = %1.2f\" % n5)\n",
    "\n",
    "# Calculating mole fractions\n",
    "ntot = n1+n2+n3+n5+n6 #total moles\n",
    "x_CO2 = n1/ntot\n",
    "x_CO = n5/ntot\n",
    "\n",
    "print(\"x_CO2 = %1.4f\" % x_CO2)\n",
    "print(\"x_CO = %1.4f\" % x_CO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least-squares appraoch\n",
    "Note that sometimes it is advantageous to be able to set bounds on our solver. Using the \"least_squares\" solver in scipy allows us to do this (note that \"fsolve\" has no such ability). *HOWEVER*, we're not guaranteed that we've found an actual answer with \"least_squares\", so we need to check the \"cost\" function, which is a measure of the residuals. This \"cost\" needs to be really small, like less than 10-5 or -6 for us to believe the answer with any certainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check that 4.644396126208125e-16 is VERY small, which implies a good \"zero\" was found.\n",
      "The temperature is  876.40 K\n",
      "x_CO2 = 0.1042\n",
      "x_CO = 0.0446\n"
     ]
    }
   ],
   "source": [
    "# Same soln approach, but using least_squares approach\n",
    "lb = [0, 0, 0, -100, -100, -100, -100, 0, 0, 0] # lower bounds\n",
    "ub = 10000*ones(10) #upper bounds\n",
    "x0 = ones(10) #initial guess\n",
    "soln = least_squares(eqn, x0, bounds = (lb,ub))\n",
    "t, T, K, d1, a1, b1, c1, n2, n6, n3 = soln.x\n",
    "print(f'Check that {soln.cost**.5} is VERY small, which implies a good \"zero\" was found.')\n",
    "\n",
    "print(\"The temperature is  %1.2f K\" % T)\n",
    "#print(\"n_CO2 = %1.2f\" % n1)\n",
    "#print(\"n_CO = %1.2f\" % n5)\n",
    "\n",
    "# Calculating mole fractions\n",
    "ntot = n1+n2+n3+n5+n6 #total moles\n",
    "x_CO2 = n1/ntot\n",
    "x_CO = n5/ntot\n",
    "\n",
    "print(\"x_CO2 = %1.4f\" % x_CO2)\n",
    "print(\"x_CO = %1.4f\" % x_CO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative Problem Statement: setting T\n",
    "\n",
    "Alternatively, we could have specified the representative reaction temperature *instead* of the product species. Note that in this example we've already set the $\\phi$ above.\n",
    "\n",
    "$T_{rxn}$ = 876 K\n",
    "\n",
    "...in this case, we don't need a multi equation solver:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phi is calculated to be 1.43\n",
      "The temperature is  876.00 K\n",
      "x_CO2 = 0.1041\n",
      "x_CO = 0.0448\n"
     ]
    }
   ],
   "source": [
    "# Set up the knowns\n",
    "fuel = 1 #kmol of C3H6\n",
    "C_fuel = 3\n",
    "H_fuel = 8\n",
    "\n",
    "phi = 1.43\n",
    "print(\"Phi is calculated to be %1.2f\" % phi)\n",
    "\n",
    "T = 876\n",
    "t = T/1000\n",
    "K = exp(2.743 - 1.761/t - 1.611/t**2 + 0.2803/t**3)\n",
    "a1 = (1-K)\n",
    "d1 = 2*a_s*(1-1/phi)\n",
    "b1 = (H_fuel/2 + K*C_fuel - d1*(1-K))\n",
    "c1 = (-C_fuel*d1*K)\n",
    "n3 = 3.76*a_s/phi\n",
    "n5 = ((-b1+(b1**2-4*a1*c1)**0.5)/(2*a1))\n",
    "n1 = (C_fuel-n5)\n",
    "n2 = (H_fuel/2 - d1 + n5)\n",
    "n6 = (d1 - n5)\n",
    "\n",
    "print(\"The temperature is  %1.2f K\" % T)\n",
    "#print(\"n_CO2 = %1.2f\" % n1)\n",
    "#print(\"n_CO = %1.2f\" % n5)\n",
    "\n",
    "# Calculating mole fractions\n",
    "ntot = n1+n2+n3+n5+n6 #total moles\n",
    "x_CO2 = n1/ntot\n",
    "x_CO = n5/ntot\n",
    "\n",
    "print(\"x_CO2 = %1.4f\" % x_CO2)\n",
    "print(\"x_CO = %1.4f\" % x_CO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%==================\n",
    "\n",
    "Now let's check our answers with cantera:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The equivalence ratio prior to equilb calc is: 1.43\n",
      "The mole fraction of CO2 is 0.1049 \n",
      "The mole fraction of CO is 0.0434 \n",
      "\n",
      "  gri30:\n",
      "\n",
      "       temperature   876 K\n",
      "          pressure   1.0132e+05 Pa\n",
      "           density   0.36265 kg/m^3\n",
      "  mean mol. weight   26.068 kg/kmol\n",
      "   phase of matter   gas\n",
      "\n",
      "                          1 kg             1 kmol     \n",
      "                     ---------------   ---------------\n",
      "          enthalpy       -1.9281e+06       -5.0262e+07  J\n",
      "   internal energy       -2.2075e+06       -5.7545e+07  J\n",
      "           entropy            8876.9         2.314e+05  J/K\n",
      "    Gibbs function       -9.7043e+06       -2.5297e+08  J\n",
      " heat capacity c_p              1331             34696  J/K\n",
      " heat capacity c_v              1012             26381  J/K\n",
      "\n",
      "                      mass frac. Y      mole frac. X     chem. pot. / RT\n",
      "                     ---------------   ---------------   ---------------\n",
      "                H2         0.0079364           0.10262           -19.463\n",
      "                 H        6.4331e-13        1.6637e-11           -9.7316\n",
      "                OH        7.4047e-15         1.135e-14           -50.296\n",
      "               H2O          0.065307          0.094501           -60.027\n",
      "               CH3        1.2485e-12        2.1647e-12            -32.21\n",
      "               CH4        0.00052509        0.00085321           -41.942\n",
      "                CO          0.046654           0.04342            -43.58\n",
      "               CO2           0.17708           0.10489           -84.144\n",
      "               HCO        8.5473e-14        7.6784e-14           -53.311\n",
      "              CH2O        3.1302e-09        2.7176e-09           -63.043\n",
      "             CH3OH        9.6454e-11        7.8471e-11           -82.506\n",
      "              C2H2        8.9313e-14        8.9416e-14           -25.494\n",
      "              C2H4        1.2236e-10         1.137e-10           -44.958\n",
      "              C2H6        4.1852e-10        3.6282e-10           -64.421\n",
      "             CH2CO        2.2624e-12         1.403e-12           -66.059\n",
      "               NH3        2.5739e-05        3.9397e-05           -41.668\n",
      "               HCN        4.8344e-08        4.6631e-08            -25.22\n",
      "              HNCO        2.6551e-08        1.6087e-08           -65.784\n",
      "                N2           0.70247           0.65367           -24.946\n",
      "            CH3CHO        2.1938e-12        1.2981e-12           -85.522\n",
      "     [  +33 minor]         1.756e-14        1.3396e-14  \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# First you need to import the requried toolboxes\n",
    "import cantera as ct\n",
    "\n",
    "# Let's create the gas object\n",
    "gas1 = ct.Solution('gri30.yaml')\n",
    "gas1.set_equivalence_ratio(phi=phi, fuel='C3H8', oxidizer={'O2':1.0, 'N2':3.76}) #note we pulled \"phi\" from above\n",
    "gas1.TP = T, 101325     #let's pull the temp from up above\n",
    "#check that we've set the equivalence ratio correctly\n",
    "g1eqr = gas1.equivalence_ratio()\n",
    "print(\"The equivalence ratio prior to equilb calc is: %1.2f\" % g1eqr)\n",
    "\n",
    "# Solve the const enthaply equation\n",
    "gas1.equilibrate('TP')\n",
    "\n",
    "# Let's search for specific products\n",
    "mf_co2 = gas1['CO2'].X\n",
    "print(\"The mole fraction of CO2 is %1.4f \" % mf_co2)\n",
    "mf_co = gas1['CO'].X\n",
    "print(\"The mole fraction of CO is %1.4f \" % mf_co)\n",
    "\n",
    "gas1()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## i-Octane\n",
    "Now, let's switch to a larger mechanism and try our iso-octane example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The equivalence ratio prior to equilb calc is: 1.43\n",
      "The mole fraction of CO2 is 0.1142 \n",
      "The mole fraction of CO is 0.0443 \n",
      "\n",
      "  gas:\n",
      "\n",
      "       temperature   846 K\n",
      "          pressure   1.0132e+05 Pa\n",
      "           density   0.38182 kg/m^3\n",
      "  mean mol. weight   26.506 kg/kmol\n",
      "   phase of matter   gas\n",
      "\n",
      "                          1 kg             1 kmol     \n",
      "                     ---------------   ---------------\n",
      "          enthalpy       -1.9452e+06        -5.156e+07  J\n",
      "   internal energy       -2.2106e+06       -5.8594e+07  J\n",
      "           entropy            8708.3        2.3082e+05  J/K\n",
      "    Gibbs function       -9.3124e+06       -2.4684e+08  J\n",
      " heat capacity c_p            1305.1             34593  J/K\n",
      " heat capacity c_v            991.43             26279  J/K\n",
      "\n",
      "                      mass frac. Y      mole frac. X     chem. pot. / RT\n",
      "                     ---------------   ---------------   ---------------\n",
      "                 H        2.0768e-13        5.4611e-12           -9.7162\n",
      "                H2         0.0073325          0.096407           -19.432\n",
      "               H2O           0.05413          0.079643           -61.269\n",
      "                N2           0.70036           0.66266            -24.85\n",
      "                CO          0.046863          0.044347            -44.02\n",
      "               CO2           0.18968           0.11424           -85.857\n",
      "              CH2O        3.1262e-09        2.7597e-09           -63.453\n",
      "               HCO        3.4865e-14        3.1847e-14           -53.736\n",
      "             HOCHO         9.898e-09        5.7003e-09           -105.29\n",
      "             CH3OH        1.2956e-10        1.0718e-10           -82.885\n",
      "               CH4         0.0016359         0.0027028           -41.049\n",
      "               CH3        1.5421e-12        2.7187e-12           -31.332\n",
      "              C2H6        3.1532e-09        2.7795e-09           -62.665\n",
      "              C2H4         4.592e-10        4.3386e-10           -43.232\n",
      "              C2H2        1.9716e-13         2.007e-13             -23.8\n",
      "            CH3CHO        1.1872e-11         7.143e-12           -85.069\n",
      "             CH2CO        6.3561e-12        4.0078e-12           -65.636\n",
      "              C3H8        1.8399e-14        1.1059e-14           -84.281\n",
      "              C3H6        3.2957e-14        2.0759e-14           -64.849\n",
      "     [ +855 minor]        1.2177e-14        8.9066e-15  \n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\baumgardner\\AppData\\Local\\Temp\\ipykernel_32680\\1526225502.py:3: UserWarning: NasaPoly2::validate: \n",
      "For species EC6H13O2, discontinuity in s/R detected at Tmid = 1388.0\n",
      "\tValue computed using low-temperature polynomial:  107.80575085204809\n",
      "\tValue computed using high-temperature polynomial: 107.11626904383155\n",
      "\n",
      "  gas2 = ct.Solution('ic8_ver3_mech.yaml')\n"
     ]
    }
   ],
   "source": [
    "# Let's try the same for iso-octane\n",
    "T2 = 846\n",
    "gas2 = ct.Solution('ic8_ver3_mech.yaml')\n",
    "gas2.set_equivalence_ratio(phi=phi, fuel='IC8H18', oxidizer={'O2':1.0, 'N2':3.76}) #note we pulled \"phi\" from above\n",
    "gas2.TP = T2, 101325     #let's pull the temp from up above\n",
    "#check that we've set the equivalence ratio correctly\n",
    "g2eqr = gas2.equivalence_ratio()\n",
    "print(\"The equivalence ratio prior to equilb calc is: %1.2f\" % g2eqr)\n",
    "\n",
    "# Solve the const enthaply equation\n",
    "gas2.equilibrate('TP')\n",
    "\n",
    "# Let's search for specific products\n",
    "mf_co2 = gas2['CO2'].X\n",
    "print(\"The mole fraction of CO2 is %1.4f \" % mf_co2)\n",
    "mf_co = gas2['CO'].X\n",
    "print(\"The mole fraction of CO is %1.4f \" % mf_co)\n",
    "\n",
    "gas2()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "502.667px",
    "left": "1070px",
    "right": "20px",
    "top": "120px",
    "width": "350px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}