{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Equilibrium by Gibbs minimization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Equilibrium criteria is\n", "\n", "$$(dG)_{T,P}=0$$\n", "\n", "$$G = G(T,P,n_i)$$\n", "\n", "$$dG = \\frac{\\partial G}{\\partial T}dT + \\frac{\\partial G}{\\partial P}dP + \\sum_i\\frac{\\partial G}{\\partial n_i}dn_i$$\n", "\n", "$$(dG)_{T,P}= \\sum_i\\frac{\\partial G}{\\partial n_i}dn_i$$\n", "\n", "* Here, we have $N_s$ species, with a given species denoted by index $i$.\n", "\n", "* Assume from now on that we have constant $T$ and $P$ and drop the $T,P$ subscripts that explicitly indicate constant $T$ and $P$.\n", "\n", "* Note the following equivalences\n", "\n", "$$dG = \\sum_i\\frac{\\partial G}{\\partial n_i}dn_i = \\sum_i\\mu_idn_i = \\nabla G\\cdot\\mathbf{dn}$$\n", "* $\\nabla$ is the gradient in the $n_i$ space.\n", "\n", "* Now, $\\mathbf{dn}$ is arbitrary, so, for $dG=0$, we have $\\nabla G=0$. That is, we are minimizing $G$ in the $n_i$ space.\n", "\n", "* At constant $T$, we can divide by $RT$ and bring inside the $\\nabla$.\n", "\n", "$$\\nabla\\frac{G}{RT} = 0$$\n", "\n", "$$\\sum_i\\frac{\\mu_i}{RT} = 0$$\n", "\n", "* Define $\\mu_i/RT$ in terms of species moles $n_i$ and total moles $n_t = \\sum_i n_i$\n", "\n", "$$\n", "\\begin{align*}\n", "\\frac{\\mu_i}{RT}&= \\frac{g_{i,f}^o}{RT} + \\ln\\left(\\frac{P_i}{P^o}\\right), \\\\\n", " &= \\underbrace{\\frac{g_{i,f}^o}{RT} + \\ln\\left(\\frac{P}{P^o}\\right)}_{\\hat{g}_i} + \\ln\\left(\\frac{n_i}{n_t}\\right).\n", "\\end{align*}\n", "$$\n", "\n", "$$\\frac{\\mu_i}{RT} = \\hat{g}_i + \\ln\\left(\\frac{n_i}{n_t}\\right)$$\n", "\n", "* **Summary: for equilibrium, solve the following for species moles $n_i$**\n", "\n", "\n", "$$\\nabla\\frac{G}{RT} = \\sum_i\\left[\\hat{g}_i + \\ln\\left(\\frac{n_i}{n_t}\\right)\\right]\\mathbf{e_i} = 0,\\phantom{xxx} i=1,\\,N_s$$\n", "\n", "Here, $\\mathbf{e_i}$ is a unit vector in the species $i$ direction.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Elemental constraints\n", "\n", "* We can't have just any $n_i$.\n", "* We need $n_i$ that conserve elements: $C$, $H$, $O$, $N$, etc.\n", " * denote elements with index $k$.\n", " * For some initial mixture (reactants), we find the equilibrium composition, subject to the constraints that elements are conserved.\n", "* $A_k$ total moles of element $k$\n", "* $a_{i,k}$ is the moles of element $k$ is species $i$\n", "* $N_e$ is the number of elements.\n", "\n", "$$A_k = \\sum_i n_ia_{i,k},\\phantom{xxx} k=1,\\,N_e$$\n", "\n", "* Constraint functions:\n", "\n", "$$c_k(n_i) = \\sum_in_ia_{i,k} -A_k= 0$$\n", "\n", "\n", "* **How do we minimize $G$, AND satisfy the elemental constraints?**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Lagrange multipliers\n", "\n", "\n", "\n", "* Minimize a function $f(x,y)$ subject to constraint $c(x,y)=0$\n", "* At the mimimum $\\nabla f = 0$\n", "* In the picture shown, $f(x,y)$ has a minimim in the center of the contours.\n", "* But subject to the constraint of being on curve $c(x,y)=0$, we want the miminimum of $f(x,y)$ along that curve.\n", "* This occurs when the isocontours of $f(x,y)$ are **tangent** to $c(x,y)=0$. \n", " * That is, the *point* where contours of $f(x,y)$ are tangent to $c(x,y)=0$.\n", "* The tangent criteria is when $\\nabla f$ and $\\nabla c$ are *in a line*\n", " * Note that $c(x,y)=0$ is our constraint *line*, but we could have a surface $z=c(x,y)$ and we consider the corresponding $\\nabla c$. Our constraint puts us on the particular contour of $z=0$ of the $c(x,y)$ surface.\n", " $$\\nabla f - \\lambda \\nabla c = 0$$\n", " * Or, $\\nabla f = \\lambda \\nabla c$, where $\\lambda$ is some scalar variable.\n", " * Vector $\\nabla f$ minus a multiple $\\lambda$ of vector $\\nabla c$ cancel. This only happens when the two vectors are parallel.\n", "* Our point is then defined by \n", "$\\nabla f - \\lambda \\nabla c = 0$, and $c(x,y)=0$, or\n", "$$\\frac{\\partial f}{\\partial x} - \\lambda \\frac{\\partial c}{\\partial x} = 0,$$\n", "$$\\frac{\\partial f}{\\partial y} - \\lambda \\frac{\\partial c}{\\partial y} = 0,$$\n", "$$c(x,y) = 0$$\n", " * This is **3 equations** in **3 unknowns: $x$, $y$, and $\\lambda$**\n", " \n", "#### More dimensions and constraints\n", "* Now, the above was 2D with 1 constraint\n", "* For 3D with 1 constraint\n", " * The $f(x,y,z)$ isocontour is tangent to the $c(x,y,z)=0$ line.\n", "* For 3D with 2 constraints\n", " * $c_1(x,y,z)=0$, $c_2(x,y,z)=0$\n", " * $f(x,y,z)$ is not necessarily tangent to $c_1$ or $c_2$, but it is tangent to the line of intersection of surfaces $c_1$ and $c_2$\n", " * In 3D $c_1=0$ and $c_2=0$ will be surfaces. When both are satisfied, we will be on the line of intersection of these two surfaces. $f$ is then minimum along this line when contours of $f$ are tangent to this line of intersection.\n", " $$\\nabla f - \\lambda_1\\nabla c_1 - \\lambda_2\\nabla c_2 = 0$$\n", "* In general, we will have at least one fewer constraint than the number of dimensions.\n", "* If we have $k$ constraints, then \n", "\n", "$$\\nabla f - \\sum_k\\lambda_k\\nabla c_k=0$$\n", "\n", "\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Equilibrium\n", "\n", "For equilibrium $f\\equiv G/RT$, and $c\\equiv c_k$, given above. Hence,\n", "$$\\nabla\\frac{G}{RT} - \\sum_k\\lambda_k\\nabla c_k= 0.$$\n", "* The gradient operators each have $N_s$ additive terms, one for each species. For the whole equation to be 0, each term will be zero.\n", "This gives\n", "$$\\hat{g}_i + \\ln\\left(\\frac{n_i}{n_t}\\right) - \\sum_k\\lambda_ka_{i,k}=0,\\phantom{xxxx}i=1,\\,N_s$$\n", "Along with the constraint equations \n", "$$\\sum_in_ia_{i,k}=A_k, \\phantom{xxx}k=1,N_e,$$\n", "and the normalization condition\n", "$$n_t=\\sum_in_i.$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Divide the second and third equations through by $n_t$, and solve the first equation for $x_i=n_i/n_t$:\n", "\n", "\n", "$$x_i = \\exp\\left(-\\hat{g}_i + \\sum_ka_{i,k}\\lambda_k\\right), \\phantom{xxx} i=1,N_s\\tag{1}$$\n", "\n", "$$\\sum_ix_ia_{i,k} = A_k/n_t, \\phantom{xxx} k=1,N_e\\tag{2}$$\n", "\n", "$$\\sum_ix_i = 1.\\tag{3}$$\n", "\n", "* **We have $N_s+N_e+1$ equations in unknowns $x_i$, $\\lambda_k$ and $n_t$**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reduce equations\n", "\n", "* We can easily have thousands of species. This would require solving thousands of coupled nonlinear equations.\n", "* Simplify by inserting Eq. (1) above into Eqs. (2) and (3):\n", "\n", "\n", "$$\\sum_ia_{i,k}\\exp\\left(-\\hat{g}_i + \\sum_{j}\\lambda_j a_{i,j}\\right) = A_k/n_t,\\phantom{xxxx} k=1,\\,N_e$$\n", "$$\\sum_i\\exp\\left(-\\hat{g}_i + \\sum_k\\lambda_ka_{i,k}\\right) = 1$$\n", "\n", "\n", "* note in the equation above we have a sum over $j$ since index $k$ is already used.\n", "* These are a set of $N_e+1$ equations in unknowns $\\lambda_k$, $n_t$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matrix form\n", "\n", "* The above equations can be written in a matrix form.\n", "\n", "\n", "$$\\mathbf{[a]}^{T}\\exp(-\\mathbf{\\hat{g}}+\\mathbf{[a]}\\boldsymbol{\\lambda}) = \\mathbf{A}/n_t,$$\n", "$$\\sum_i\\exp(-\\mathbf{\\hat{g}}+\\mathbf{[a]}\\boldsymbol{\\lambda}) = 1.$$\n", "\n", "\n", "\n", "$$\\mathbf{x} = \\exp(-\\mathbf{\\hat{g}} + \\mathbf{[a]}\\boldsymbol{\\lambda})$$\n", "\n", "\n", "* $\\mathbf{A}$, and $\\boldsymbol{\\lambda}$ are vectors of length $N_e$, and $\\mathbf{\\hat{g}}$, $\\mathbf{x}$ are vectors of length $N_s$; $\\mathbf{[a]}$ is the $N_s\\times N_e$ matrix of elemental compositions of the species.\n", "* For the species $CH_4, CO_2, CO, H_2O, N_2, H_2$, we have $\\mathbf{[a]} = $\n", " \n", "| | C | H | O | N |\n", "|--------|---|---|---|---|\n", "| $CH_4$ | 1 | 4 | 0 | 0 |\n", "| $CO_2$ | 1 | 0 | 2 | 0 |\n", "| $CO$ | 1 | 0 | 1 | 0 |\n", "| $H_2O$ | 0 | 2 | 1 | 0 |\n", "| $N_2$ | 0 | 0 | 0 | 2 |\n", "| $H_2$ | 0 | 2 | 0 | 0 |\n", "\n", "* We solve the first two blue equations for $\\boldsymbol{\\lambda}$ and $n_t$. Then we solve the green equation for the species vector $\\mathbf{x}$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Element potentials\n", "\n", "Consider the green equation\n", "\n", "$$\\mathbf{x} = \\exp(-\\mathbf{\\hat{g}} + \\mathbf{[a]}\\boldsymbol{\\lambda})$$\n", "\n", "\n", "We can rearrange this to \n", "$$\\mathbf{[a]}\\boldsymbol{\\lambda} = \\mathbf{\\hat{g}} + \\ln(x) = \\boldsymbol{\\mu}/RT $$\n", "\n", "The first and last parts are\n", "$$\\mathbf{[a]}\\boldsymbol{\\lambda} = \\boldsymbol{\\mu}/RT $$\n", "\n", "Chemical potentials $\\mu_i$ are linear combinations of element potentials $\\lambda_k$.\n", "$$\\sum_k\\lambda_ka_{i,k} = \\mu_i/RT$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution approach\n", "\n", "* Solve the above two blue equations for $\\lambda_k$ and $n_t$.\n", " * Use Newton's method.\n", "* Then, solve for species composition $x_i$ using the green equation above.\n", "* This approach can be used at a given $T$ and $P$. If enthalpy is known instead of $T$, then we can either add an enthalpy equation like $H = H(T,n_i)$, or we can guess temperature solve equilibrium, and then do an iteration on this process to find T that satisfies $H=H(T,x_i)$.\n", "\n", "* We need an initial guess for the $\\lambda_k$.\n", " * If we guess $\\mathbf{x}$, then we can invert the green equation to find $\\lambda$. But there are more $x_i$ than $\\lambda_k$. That is, there are more constraints than degrees of freedom, so we use linear least squares.\n", " * For generic $Ax=b$, where the length of $b$ is greater than the length of $x$ (and $A$ is rectangular), the linear least squares solution is $A^TAx = A^Tb$. We solve this for $x$ given $b$.\n", " * For our problem, we have $[a]^T[a](\\lambda) = [a]^T(\\hat{g}+\\ln(x))$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solver\n", "* Use Cantera for thermochemical properties\n", "* You will need [streams.py](http://ignite.byu.edu/che633/lectures/streams.py).\n", " * This has greek letters in the text, so right click and save as, don't copy and paste text from the browser.\n", "* You also need to have cantera installed. [See this link](http://www.cantera.org/docs/sphinx/html/install.html).\n", "\n", "### Analytic Jacobian\n", "Nonlinear solvers based on Newton's method use the Jacobian matrix. This can be computed numerically, but when an analytical Jacobian is available it is best to use it.\n", "\n", "Write our equations as $f$ and $h$:\n", "\n", "\n", "$$f_k = \\sum_ia_{i,k}\\exp\\left(-\\hat{g}_i + \\sum_{j}\\lambda_j a_{i,j}\\right) -A_k/n_t= 0,\\phantom{xxxx} k=1,\\,N_e$$\n", "$$h = \\sum_i\\exp\\left(-\\hat{g}_i + \\sum_k\\lambda_ka_{i,k}\\right) -1 = 0$$\n", "\n", "\n", "The Jacobian matrix is then\n", "$$J = \\begin{pmatrix}\n", "\\frac{\\partial f_1}{\\partial \\lambda_1} & \n", "\\ldots &\n", "\\frac{\\partial f_1}{\\partial \\lambda_{N_e}} & \n", "\\frac{\\partial f_1}{\\partial n_t} & \n", "\\\\\n", "\\vdots & \\ldots & \\vdots & \\vdots \\\\\n", "\\\\\n", "\\frac{\\partial f_{N_e}}{\\partial \\lambda_1} & \n", "\\ldots &\n", "\\frac{\\partial f_{N_e}}{\\partial \\lambda_{N_e}} & \n", "\\frac{\\partial f_{N_e}}{\\partial n_t} & \n", "\\\\\n", "\\frac{\\partial h}{\\partial \\lambda_1} &\n", "\\ldots &\n", "\\frac{\\partial h}{\\partial \\lambda_{N_e}} &\n", "\\frac{\\partial h}{\\partial n_t}\n", "\\end{pmatrix}$$\n", "\n", "The four main components are given by\n", "\n", "$$\\begin{align}\n", "&\\frac{\\partial f_k}{\\partial\\lambda_m} = \\sum_ia_{i,k}\\exp\\left(-\\hat{g}_i+\\sum_ja_{i,j}\\lambda_j\\right)a_{i,m}, \\\\\n", "&\\frac{\\partial f_k}{\\partial n_t } = A_k/n_t^2, \\\\\n", "&\\frac{\\partial h}{\\partial \\lambda_m} = \\sum_i\\exp\\left(-\\hat{g}_i+\\sum_ka_{i,k}\\lambda_k\\right)a_{i,m}, \\\\\n", "&\\frac{\\partial h}{\\partial n_t} = 0.\n", "\\end{align}$$\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [], "source": [ "import numpy as np\n", "import cantera as ct\n", "from scipy.optimize import fsolve\n", "import matplotlib.pyplot as plt\n", "from streams import streams" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Species XCantera XSolver\n", "-----------------------------------\n", "CH4 5.137512e-09 5.137512e-09\n", "O2 2.846952e-11 2.846952e-11\n", "N2 5.685436e-01 5.685436e-01\n", "CO2 3.037884e-02 3.037884e-02\n", "H2O 1.282186e-01 1.282186e-01\n", "CO 1.134398e-01 1.134398e-01\n", "H2 1.594184e-01 1.594184e-01\n", "OH 6.834862e-07 6.834862e-07\n", "O 7.735590e-11 7.735590e-11\n" ] } ], "source": [ "ξ = 0.1\n", "T = 1600.0\n", "P = 101325\n", "\n", "strm = streams({\"O2\":1,\"N2\":3.76}, {\"CH4\":1}, 300, 300, 101325, \"./simple.yaml\")\n", "#strm = streams({\"O2\":1,\"N2\":3.76}, {\"CH4\":1}, 300, 300, 101325, \"gri30.yaml\")\n", "\n", "spNames = [\"CH4\", \"O2\", \"N2\", \"CO2\", \"H2O\", \"CO\", \"H2\", \"OH\", \"O\"]\n", "elNames = ['C', 'H', 'O', 'N']\n", "\n", "isp = [strm.gas.species_index(i) for i in spNames]\n", "Nsp = len(spNames)\n", "Nel = len(elNames)\n", "\n", "xmix = strm.get_x_from_y(strm.get_mixing_state(ξ)[0])\n", "strm.gas.TPX = T,P,xmix\n", "\n", "#------------ get a\n", "\n", "a = np.array([[strm.gas.n_atoms(isp[i],elNames[k]) for k in range(Nel)] \n", " for i in range(Nsp)])\n", "\n", "#------------ get A\n", "\n", "A = np.dot(a.T,xmix[isp])\n", "\n", "#------------ get ghat\n", "\n", "ghat = (strm.gas.standard_gibbs_RT + np.log(P/101325))[isp]\n", "\n", "#------------ solver function\n", "\n", "def FeqTP(λnt):\n", " λ = λnt[:-1]\n", " nt = λnt[-1]\n", " \n", " F = np.zeros(len(λnt))\n", " F[:-1] = np.dot(a.T, np.exp(-ghat+np.dot(a,λ))) - A/nt\n", " F[-1] = np.sum(np.exp(-ghat+np.dot(a,λ))) - 1.0\n", " \n", " return F\n", "\n", "#------------ Jacobian function\n", "\n", "def getJac(λnt):\n", " λ = λnt[:-1]\n", " nt = λnt[-1]\n", " \n", " n = len(λnt)\n", " J = np.zeros((n,n))\n", " \n", " egaλ = np.exp(-ghat+np.dot(a,λ))\n", " \n", " for k in range(n-1):\n", " for m in range(n-1):\n", " J[k,m] = np.sum(a[:,k]*a[:,m]*egaλ)\n", " m = n-1\n", " J[k,m] = A[k]/(nt*nt)\n", " k = n-1\n", " for m in range(n-1):\n", " J[k,m] = np.sum(a[:,m]*egaλ)\n", " m = n-1\n", " J[k,m] = 0.0\n", " \n", " return J\n", "\n", "#------------ initial guesses for nt, and λ (make sure consistent with A above), then solve\n", "\n", "x = strm.get_pCC(ξ, getYorX='x')[isp] + 1E-15 # products of complete combustion as a guess.\n", "λ = np.linalg.solve(np.dot(a.T,a), np.dot(a.T, ghat+np.log(x))) # linear least squares\n", "λnt_guess = np.append(λ, 1)\n", "\n", "λnt = fsolve(FeqTP, λnt_guess, factor=0.001, maxfev=10000, fprime=getJac)\n", "\n", "#------------ recover composition\n", "\n", "x = np.exp(-ghat + np.dot(a,λnt[:-1]))\n", "\n", "#------------ compare to Cantera\n", "\n", "xmix = strm.get_x_from_y(strm.get_mixing_state(ξ)[0])\n", "strm.gas.TPX = T,P,xmix\n", "strm.gas.equilibrate(\"TP\")\n", "xCantera = strm.gas.X\n", "\n", "#------------ output results\n", "\n", "print(\"Species XCantera XSolver\")\n", "print(\"-----------------------------------\")\n", "for i in range(Nsp):\n", " print(f\"{spNames[i]:8s} {xCantera[isp][i]:.6e} {x[i]:.6e}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "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.9.7" } }, "nbformat": 4, "nbformat_minor": 4 }