From 4e8e30fe92fce9c216f0b2944f5db20acf2067db Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Fri, 5 Apr 2024 13:01:43 -0700 Subject: [PATCH] homog_params with irreversible reactions and specified reaction order --- docs/source/conf.py | 7 +++- docs/source/user_guide/homog_params.rst | 43 +++++++++++++++++++++---- echemfem/solver.py | 10 ++++-- examples/catalyst_layer_Cu_full.py | 2 +- 4 files changed, 50 insertions(+), 12 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 0b9a55b..f92d10b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -3,7 +3,7 @@ # -- Project information project = 'EchemFEM' -copyright = u'2022-2023, Lawrence Livermore National Laboratory' +copyright = u'2022-2024, Lawrence Livermore National Laboratory' author = 'Thomas Roy' release = '0.1' @@ -21,6 +21,11 @@ 'sphinx.ext.mathjax', ] +mathjax3_config = { + 'loader': {'load': ['[tex]/mhchem']}, + 'tex': {'packages': {'[+]': ['mhchem']}}, +} + intersphinx_mapping = { 'python': ('https://docs.python.org/3/', None), 'numpy': ('https://docs.scipy.org/doc/numpy/', None), diff --git a/docs/source/user_guide/homog_params.rst b/docs/source/user_guide/homog_params.rst index 0f9265a..cbdebac 100644 --- a/docs/source/user_guide/homog_params.rst +++ b/docs/source/user_guide/homog_params.rst @@ -7,17 +7,20 @@ each reaction. Below is a list of different keys that should appear in each dictionary * :Key: ``"stoichiometry"`` - :Type: :py:class:`dict` + :Type: :py:class:`dict` of :py:class:`int` :Description: This entry defines the stoichiometry of the reaction. Each key in this dictionary should be the name of a species (as defined in the ``"name"`` key values of ``conc_params``). The corresponding value is the stoichemetry coefficient for the species in this reaction. It should be an integer: negative for a reactant, positive for a product. * :Key: ``"forward rate constant"`` :Type: :py:class:`float` or Firedrake expression :Description: Rate constant for the forward reaction. -* :Key: ``"backward rate constant"`` or ``"equilibrium constant"`` +* :Key: ``"backward rate constant"`` or ``"equilibrium constant"`` (optional) :Type: :py:class:`float` or Firedrake expression - :Description: Rate constant for the backward reaction or equilibrium constant of the reaction. Only of one these is required. + :Description: Rate constant for the backward reaction or equilibrium constant of the reaction. For a reversible reaction, only of one these is required. For an irreversible reaction, leave unset. * :Key: ``"reference concentration"`` (optional) :Type: :py:class:`float` or Firedrake expression :Description: Value of 1M in the appropriate units. By setting this, the rate constants are assumed to have the same units (typically 1/s), and the equilibrium constant to be nondimensional. +* :Key: ``"reaction order"`` (optional) + :Type: :py:class:`dict` of :py:class:`int` + :Description: This entry can be used to enforce reaction order. Each key in this dictionary should be the name of a species (as defined in the ``"name"`` key values of ``conc_params``). The corresponding value (positive integer) is reaction order of the species in this reaction. If unset, the absolute value of the stoichiometry coeffcient is the reaction order. Assuming an arbitrary homogeneous reaction system where the forward and backward rate constants for reaction :math:`i` are :math:`k_i^f` and @@ -32,6 +35,9 @@ of mass action, the reaction for species :math:`j` is given by where :math:`c_n` is the concentration of species :math:`n`. If the equilibrium constant :math:`K_i` is provided instead of the backward rate constant, then it is automatically recovered via :math:`k_i^b = \dfrac{k_i^f}{K_i}`. +The above formula assumes single-step reactions so that the reaction orders are +given by the stoichiometry coefficients. If that is not the case, the exponents +can be replaced by the ``"reaction order"`` inputs. Note that here the rate constants can have different units, depending on the number of reactants or products. It is also common to write the reactions in @@ -54,9 +60,9 @@ Here is an example of a reaction system that can be implement via .. math:: - \mathrm{CO}_2 + \mathrm{OH}^- \xrightarrow[k_1^b]{k_1^f} \mathrm{HCO}_3^- + \mathrm{CO}_2 + \mathrm{OH}^- \xrightleftharpoons[k_1^b]{k_1^f} \mathrm{HCO}_3^- - \mathrm{HCO}_3^- + \mathrm{OH}^- \xrightarrow[k_2^b]{k_2^f} \mathrm{CO}_3^{2-} + \mathrm{HCO}_3^- + \mathrm{OH}^- \xrightleftharpoons[k_2^b]{k_2^f} \mathrm{CO}_3^{2-} + \mathrm{H}_2\mathrm{O} Here, we assume the dilute solution case where :math:`\mathrm{H}_2\mathrm{O}` concentration is not tracked. In ``conc_params``, we will have entries for the @@ -79,8 +85,31 @@ list: "backward rate constant": k2b }] +Now here is another bicarbonate approximation with a high pH approximation, +where the reaction system is simplified into a single irreversible reaction: + +.. math:: + + \mathrm{CO}_2 + 2\mathrm{OH}^- \xrightarrow[]{k_1^f} \mathrm{CO}_3^{2-} + \mathrm{H}_2\mathrm{O} + +Since this is not actually a single-step reaction, we need to specify the +reaction order for :math:`\mathrm{OH}^-`. We can use the following +``homog_params`` list: + +.. code:: + + [{"stoichiometry": {"CO2": -1, + "OH": -2, + "CO3": 1}, + "reaction order": {"OH": 1}, + "forward rate constant": k1f + }] + + As an alternative to the ``homog_params`` interface, the reactions can be written directly and passed as a function in ``physical_params["bulk -reaction"]`` (:doc:`physical_params`). See ``examples/carbonate.py`` and +reaction"]`` (:doc:`physical_params`). See ``examples/cylindrical_pore.py`` for +the direct implementation of the high pH approximation of the bicarbonate bulk +reactions. See ``examples/carbonate.py`` and ``examples/carbonate_homog_params.py`` for two equivalent implementations of -the bicarbonate bulk reactions. +the full bicarbonate bulk reactions. diff --git a/echemfem/solver.py b/echemfem/solver.py index 2f02965..1146d25 100644 --- a/echemfem/solver.py +++ b/echemfem/solver.py @@ -1223,7 +1223,7 @@ def bulk_reaction(u): reaction_f[idx] = reaction["forward rate constant"] if reaction.get("backward rate constant"): reaction_b[idx] = reaction["backward rate constant"] - else: + elif reaction.get("equilibrium constant"): reaction_b[idx] = reaction["forward rate constant"] / \ reaction["equilibrium constant"] if reaction.get("reference concentration"): @@ -1232,11 +1232,15 @@ def bulk_reaction(u): c_ref = 1.0 for name in reaction["stoichiometry"]: s = reaction["stoichiometry"][name] + order = abs(s) # default + if reaction.get("reaction order"): + if reaction["reaction order"].get(name): + order = reaction["reaction order"][name] a = u[self.i_c[name]] / c_ref if s < 0: - reaction_f[idx] *= a ** (-s) + reaction_f[idx] *= a ** order if s > 0: - reaction_b[idx] *= a ** s + reaction_b[idx] *= a ** order for i in self.idx_c: name = self.conc_params[i]["name"] diff --git a/examples/catalyst_layer_Cu_full.py b/examples/catalyst_layer_Cu_full.py index 85d5b2c..bcf77ae 100644 --- a/examples/catalyst_layer_Cu_full.py +++ b/examples/catalyst_layer_Cu_full.py @@ -1,4 +1,4 @@ -# import matplotlib.gridspec as gridspec +import matplotlib.gridspec as gridspec import matplotlib.pyplot as plt from firedrake import * from echemfem import EchemSolver