{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beispiel 3.4: Reaktionsfortschritt bei komplexen Reaktionen\n",
    "Bearbeitet von Franz Braun\n",
    "\n",
    "Dieses Beispiel befindet sich im Lehrbuch auf den Seiten 20 - 22. Das hier angewendete\n",
    "Vorgehen entspricht dem im Lehrbuch vorgestellten Lösungsweg.\n",
    "\n",
    "Zunächst werden die benötigten Pakete importiert."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Import\n",
    "import numpy as np\n",
    "from scipy.optimize import root\n",
    "from tabulate import tabulate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In diesem Beispiel wird die Gleichgewichtszusammensetzung der Methanolsynthese betrachtet.\n",
    "\n",
    "\\begin{align*}\n",
    "&\\mathrm{1.\\,\\, Reaktion} \\quad &   \\mathrm{CO + 2\\,H_2}    &       \\leftrightarrow     \\mathrm{CH_3OH} \\\\\n",
    "&\\mathrm{2.\\,\\, Reaktion} \\quad &   \\mathrm{CO_2 + 3\\,H_2}  &       \\leftrightarrow     \\mathrm{CH_3OH + H_2O} \\\\\n",
    "&\\mathrm{3.\\,\\, Reaktion} \\quad &   \\mathrm{CO + H_2O}      &       \\leftrightarrow     \\mathrm{CO_2 + H_2} \\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zwei dieser Reaktionen (z.B. R1 und R2) sind linear unabhängig (vgl. Text zu Beispiel 3.4 im Lehrbuch). Deshalb müssen nur die Massenwirkungsgesetze dieser Reaktionen berücksichtigt werden. \n",
    "\n",
    "\\begin{align*}\n",
    "K_{p,1} &= \\frac{p_\\mathrm{CH_3OH}}{p_\\mathrm{CO} \\, p_\\mathrm{H_2}^2}\\\\\n",
    "K_{p,2} &= \\frac{p_\\mathrm{CH_3OH} \\, p_\\mathrm{H_2O}}{p_\\mathrm{CO_2} \\, p_\\mathrm{H_2}^3}\n",
    "\\end{align*}\n",
    "\n",
    "Anschließend werden die Gleichgewichtskonstanten und die Reaktionsbedingungen ($T,\\,p$) parametriert."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "T       = 500       # K\n",
    "p       = 45        # bar\n",
    "K_p1    = 12.0646   # bar-2\n",
    "K_p2    = 0.1678    # bar−2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analog zum Beispiel 3.3 kann die Definition der Reaktionslaufzahl $\\xi$ mit der gegebenen Ausgangszusammensetzung der Produkte angewendet werden:\n",
    "\n",
    "\\begin{align*}\n",
    "n_\\mathrm{CO} &= 1 \\,\\mathrm{mol} - \\xi_1,\\\\\n",
    "n_\\mathrm{CO_2} &= 0.1 \\,\\mathrm{mol} - \\xi_2,\\\\\n",
    "n_\\mathrm{H_2} &= 2 \\,\\mathrm{mol} - 2\\,\\xi_1 - 3\\,\\xi_2,\\\\\n",
    "n_\\mathrm{CH_3OH} &=  \\xi_1 + \\xi_2,\\\\\n",
    "n_\\mathrm{H_2O} &= \\xi_2,\\\\\n",
    "\\sum n_i &= 3,1 \\,\\mathrm{mol} - 2\\,\\xi_1 - 2\\,\\xi_2.\\\\\n",
    "\\end{align*}\n",
    "\n",
    "Mit\n",
    "\\begin{align*}\n",
    "p_i &= x_i \\, p,\\\\\n",
    "x_i &= \\frac{n_i}{n},\n",
    "\\end{align*}\n",
    "ergeben sich folgende Bestimmungsgleichungen:\n",
    "\n",
    "\\begin{align*}\n",
    "f_1(\\xi_1,\\xi_2) &= K_{p,1} \\, p^2 \\, n_\\mathrm{CO} \\, n_\\mathrm{H_2}^2 - n_\\mathrm{CH_3OH} \\, n^2 \\qquad \\enspace = 0,\\\\\n",
    "f_2(\\xi_1,\\xi_2) &= K_{p,2} \\, p^2 \\, n_\\mathrm{CO_2} \\, n_\\mathrm{H_2}^3 - n_\\mathrm{CH_3OH} \\, n_\\mathrm{H_2O} \\, n^2 = 0.\n",
    "\\end{align*}\n",
    "\n",
    "Um die Nullstellen zu bestimmten, werden $f_1(\\xi_1,\\xi_2)$ und $f_2(\\xi_1,\\xi_2)$ als Funktion definiert. Anschließend wird mit _scipy.optimize.root_ nach den Nullstellen gesucht."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_xi(xi):\n",
    "    '''\n",
    "    Umgeformtes Massenwirkungsgesetz\n",
    "     -> Funktion zur Bestimmung der Nullstellen \n",
    "    \n",
    "    Parameter:\n",
    "    ----------\n",
    "    xi : float\n",
    "        Reaktionslaufzahlen in mol\n",
    "        xi[0] -> xi_1\n",
    "        xi[1] -> xi_2\n",
    "    \n",
    "    n_i : float\n",
    "        Stoffmengenanteile in mol\n",
    "\n",
    "    '''\n",
    "\n",
    "    xi_1        = xi[0]\n",
    "    xi_2        = xi[1]\n",
    "\n",
    "    # Stoffmengenanteile\n",
    "    n_CO        = 1         -   xi_1\n",
    "    n_CO2       = 0.1       -   xi_2\n",
    "    n_H2        = 2         -   2 * xi_1    -     3 * xi_2\n",
    "    n_CH3OH     =               xi_1        +     xi_2\n",
    "    n_H2O       =                               xi_2\n",
    "\n",
    "    n           = 3.1       - 2 * xi_1      - 2 * xi_2 \n",
    "\n",
    "    # Bestimmungsgleichungen\n",
    "    f_1 = K_p1 * p**2 * n_CO  * n_H2**2 - n_CH3OH *         n**2\n",
    "    f_2 = K_p2 * p**2 * n_CO2 * n_H2**3 - n_CH3OH * n_H2O * n**2\n",
    "\n",
    "    return np.hstack((f_1,f_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Das Lösen war erfolgreich.\n",
      "Für die Reaktionslaufzahl 1 wurde der Wert  0.974  mol ermittelt.\n",
      "Für die Reaktionslaufzahl 2 wurde der Wert  0.0023  mol ermittelt.\n"
     ]
    }
   ],
   "source": [
    "xi_guess = np.array((0, 0.))          # \"initial guess\" für xi_1 und xi_2\n",
    "sol = root(f_xi, xi_guess, tol = 1e-5)  # Lösen mit scypi.optimize.root\n",
    "\n",
    "''' \n",
    "Wenn der solver erfolgreich ist, wird xi ausgegeben.\n",
    "Sollte das Lösen fehlschlagen, wird die \"solver-message\" und der \"solver-success\" ausgegeben.\n",
    "'''\n",
    "if sol.success:\n",
    "    xi_1 = sol.x[0]\n",
    "    xi_2 = sol.x[1]\n",
    "    print('Das Lösen war erfolgreich.')\n",
    "    print('Für die Reaktionslaufzahl 1 wurde der Wert ', np.round(xi_1,4),' mol ermittelt.')\n",
    "    print('Für die Reaktionslaufzahl 2 wurde der Wert ', np.round(xi_2,4),' mol ermittelt.')\n",
    "else:\n",
    "    print(sol.success)\n",
    "    print(sol.message)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternativ ist auch hier eine Lösung mit dem NEWTONschen Verfahren, wie in Beispiel 3.3 demonstriert, möglich. Diese Vorgehensweise wird hier aber nicht gezeigt.\n",
    "\n",
    "Abschließend werden die Stoffmengen im Gleichgewicht bestimmt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_CO_eq         = 1     -   xi_1\n",
    "n_CO2_eq        = 0.1   -   xi_2\n",
    "n_H2_eq         = 2     -   2 * xi_1        -     3 * xi_2\n",
    "n_CH3OH_eq      =               xi_1        +     xi_2\n",
    "n_H2O_eq        =                                 xi_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              CO     CO2     H2    CH3OH     H2O\n",
      "---------  -----  ------  -----  -------  ------\n",
      "n_i / mol  0.026  0.0977  0.045   0.9763  0.0023\n"
     ]
    }
   ],
   "source": [
    "# Ausgabe der Ergebnisse als Tabelle\n",
    "\n",
    "names_y     = np.array((['n_i / mol']))\n",
    "n_row       = np.array((n_CO_eq, n_CO2_eq, n_H2_eq, n_CH3OH_eq, n_H2O_eq))  \n",
    "\n",
    "# Array für die Tabelle\n",
    "table       = [np.hstack((names_y, np.round(n_row, 4)))]\n",
    "\n",
    "print(tabulate(table, headers = ['CO', 'CO2', 'H2', 'CH3OH', 'H2O'] ))"
   ]
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}