homog_params input to define bulk reactions

docs/source/user_guide/homog_params.rst

@@ -0,0 +1,86 @@
+Homogeneous Reaction Parameters
+===============================
+
+The parameters of homogeneous/bulk reactions can be provided through
+:attr:`echemfem.EchemSolver.homog_params`, a list containing one dictionary for
+each reaction. Below is a list of different keys that should appear in each
+dictionary
+
+* :Key: ``"stoichiometry"``
+  :Type: :py:class:`dict`
+  :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"``
+  :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.
+* :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.
+
+Assuming an arbitrary homogeneous reaction system where the forward and
+backward rate constants for reaction :math:`i` are :math:`k_i^f` and
+:math:`k_i^b`, respectively, and the stoichiometry coefficient of species
+:math:`j` in reaction :math:`i` is :math:`s_{j,i}`, then, according to the law
+of mass action, the reaction for species :math:`j` is given by
+
+.. math::
+
+   R_j = \sum_i s_{j, i} \left( k_i^f \prod_{n, s_{n, i} < 0} c_n^{-s_{n, i}} - k_i^b \prod_{n, s_{n, i} > 0} c_n^{s_{n, i}} \right),
+
+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}`. 
+
+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
+terms of dimensionless "activities" instead of concentrations. Then, the rate
+constants have the same units (typically 1/s) and the equilibrium constants are
+dimensionless. With this definition of rate constants, we instead write the
+reaction for species :math:`j` as
+
+
+.. math::
+
+   R_j = \sum_i s_{j, i} c_\mathrm{ref} \left( k_i^f \prod_{n, s_{n, i} < 0} a_n^{-s_{n, i}} - k_i^b \prod_{n, s_{n, i} > 0} a_n^{s_{n, i}} \right),
+
+where :math:`a_n = c_n / c_\mathrm{ref}` is the activity of species :math:`n`
+and :math:`c_\mathrm{ref} = 1\text{M}` is a reference concentration, which
+needs to be provided in the correct units to ``homog_params``.
+
+Here is an example of a reaction system that can be implement via
+``homog_params``. Simplified CO2 - bicarbonate homogeneous reactions:
+
+.. math::
+
+   \mathrm{CO}_2 + \mathrm{OH}^- \xrightarrow[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-}
+
+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
+other species in this reactions, with ``"name"`` key values: ``"CO2"``,
+``"OH"``, ``"HCO3"``, and ``"CO3"``.  We get the following ``homog_params``
+list:
+
+.. code::
+
+    [{"stoichiometry": {"CO": -1,
+                        "OH": -1,
+                        "HCO3": 1},
+      "forward rate constant": k1f,
+      "backward rate constant": k1b
+     },
+     {"stoichiometry": {"HCO3": -1,
+                        "OH": -1,
+                        "CO3": 1},
+      "forward rate constant": k2f,
+      "backward rate constant": k2b
+     }]
+
+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
+``examples/carbonate_homog_params.py`` for two equivalent implementations of
+the bicarbonate bulk reactions.

docs/source/user_guide/index.rst

@@ -25,6 +25,7 @@ To create a customized solver, the user should set the following inputs:
    physical_params
    boundary_conditions
    echem_params
+   homog_params
 
 Here is a barebone example of how a user might define their own solver.
 Several concrete examples can be found in `echemfem/examples

docs/source/user_guide/physical_params.rst

@@ -103,7 +103,7 @@ Below is a list of other keys that can appear in the dictionary
 
 * :Key: ``"bulk reaction"``
   :Type: a function
-  :Description: Homogeneous/bulk reactions to be added to the right-hand side of mass conservation equations.
+  :Description: Homogeneous/bulk reactions to be added to the right-hand side of mass conservation equations. These can instead be set using ``homog_params`` (:doc:`homog_params`).
 
     Args:
         u: solution state. The value of the different concentrations can be recovered through ``u([self.i_c["species name"])`` within a ``echemfem.EchemSolver`` object.

echemfem/solver.py

@@ -60,6 +60,7 @@ def __init__(
             physical_params,
             mesh,
             echem_params=[],
+            homog_params=[],
             gas_params=[],
             stats_file="",
             overwrite_stats_file=False,
@@ -72,6 +73,7 @@ def __init__(
         self.physical_params = physical_params
         self.conc_params = conc_params
         self.echem_params = echem_params
+        self.homog_params = homog_params
         self.gas_params = gas_params
         self.num_c = len(conc_params)
         self.num_g = len(gas_params)
@@ -380,6 +382,9 @@ def _internal_dS(subdomain_id=None, domain=None):
         elif self.flow["advection"]:
             self.set_velocity()
 
+        if self.homog_params:
+            self.setup_bulk_reactions()
+
         self.setup_forms(us, v)
 
         self.u = u
@@ -1205,6 +1210,51 @@ def init_solver_parameters(
         solver_parameters.update(pc)
         self.solver_parameters = solver_parameters
 
+    def setup_bulk_reactions(self):
+        """ Setting up the homogeneous bulk reactions using homog_params
+
+        Setup using the law of mass action
+        """
+        def bulk_reaction(u):
+            reactions = [0] * self.num_liquid
+            reaction_f = [0] * len(self.homog_params)
+            reaction_b = [0] * len(self.homog_params)
+            for idx, reaction in enumerate(self.homog_params):
+                reaction_f[idx] = reaction["forward rate constant"]
+                if reaction.get("backward rate constant"):
+                    reaction_b[idx] = reaction["backward rate constant"]
+                else:
+                    reaction_b[idx] = reaction["forward rate constant"] / \
+                                        reaction["equilibrium constant"]
+                if reaction.get("reference concentration"):
+                    c_ref = reaction["reference concentration"]
+                else:
+                    c_ref = 1.0
+                for name in reaction["stoichiometry"]:
+                    s = reaction["stoichiometry"][name]
+                    a = u[self.i_c[name]] / c_ref
+                    if s < 0:
+                        reaction_f[idx] *= a ** (-s)
+                    if s > 0:
+                        reaction_b[idx] *= a ** s
+
+            for i in self.idx_c:
+                name = self.conc_params[i]["name"]
+                for idx, reaction in enumerate(self.homog_params):
+                    if name in reaction["stoichiometry"]:
+                        s = reaction["stoichiometry"][name]
+                        if reaction.get("reference concentration"):
+                            c_ref = reaction["reference concentration"]
+                        else:
+                            c_ref = 1.0
+                        reactions[i] += s * c_ref * (reaction_f[idx] - reaction_b[idx])
+
+            return reactions
+        if self.physical_params.get("bulk reaction"):
+            print = PETSc.Sys.Print
+            print("WARNING: Overwriting bulk reaction from physical_params with the one in homog_params")
+        self.physical_params["bulk reaction"] = bulk_reaction
+
     def effective_diffusion(self, D, phase="liquid"):
         """ Bruggeman correlation for effective diffusion in a porous medium
         """

examples/carbonate_homog_params.py

@@ -0,0 +1,183 @@
+from firedrake import *
+from echemfem import EchemSolver, IntervalBoundaryLayerMesh
+
+
+class CarbonateSolver(EchemSolver):
+    """
+    Steady-state version of example from
+    Gupta, N., Gattrell, M. and MacDougall, B., 2006. Calculation for the
+    cathode surface concentrations in the electrochemical reduction of CO2 in
+    KHCO3 solutions. Journal of applied electrochemistry, 36(2), pp.161-172.
+    but using the bicarbonate bulk reactions from
+    Schulz, K.G., Riebesell, U., Rost, B., Thoms, S. and Zeebe, R.E., 2006.
+    Determination of the rate constants for the carbon dioxide to bicarbonate
+    inter-conversion in pH-buffered seawater systems. Marine chemistry,
+    100(1-2), pp.53-65.
+    """
+
+    def __init__(self):
+
+        delta = 0.0001
+        mesh = IntervalBoundaryLayerMesh(100, delta, 100, 1e-7)
+
+        # Reaction Rates in SI units
+        k1f = 2.23  # m3/mol.s
+        k1r = 9.71e-5  # 1/s
+        k2f = 3.71e-2  # 1/s
+        k2r = 2.67e1  # m3/mol.s
+        k3f = 3.06e5  # 1/s
+        k3r = 6.0e6  # m3/mol.s
+        k4f = 5.0e7  # m3/mol.s
+        k4r = 59.44  # 1/s
+        k5f = 2.31e7  # m3/mols. (Note - typo in Schulz)
+        k5r = 1.4  # mol/m3.s
+        # Bulk Conditions
+        C_HCO3_bulk = 500.  # mol/m3
+        C_CO32_bulk = 6.5  # mol/m3
+        C_OH_bulk = k3f / k3r * C_CO32_bulk / C_HCO3_bulk  # mol/m3
+        C_H_bulk = (k5r / k5f) / C_OH_bulk  # mol/m3
+        C_CO2_bulk = (k1r / k1f) * C_HCO3_bulk / C_OH_bulk  # mol/m3
+
+        homog_params = []
+
+        homog_params.append({"stoichiometry": {"CO2": -1,
+                                               "OH": -1,
+                                               "HCO3": 1},
+                             "forward rate constant": k1f,
+                             "backward rate constant": k1r
+                             })
+
+        homog_params.append({"stoichiometry": {"CO2": -1,
+                                               "H": 1,
+                                               "HCO3": 1},
+                             "forward rate constant": k2f,
+                             "backward rate constant": k2r
+                             })
+
+        homog_params.append({"stoichiometry": {"CO3": -1,
+                                               "OH": 1,
+                                               "HCO3": 1},
+                             "forward rate constant": k3f,
+                             "backward rate constant": k3r
+                             })
+
+        homog_params.append({"stoichiometry": {"CO3": -1,
+                                               "H": -1,
+                                               "HCO3": 1},
+                             "forward rate constant": k4f,
+                             "backward rate constant": k4r
+                             })
+
+        homog_params.append({"stoichiometry": {"OH": -1,
+                                               "H": -1},
+                             "forward rate constant": k5f,
+                             "backward rate constant": k5r
+                             })
+
+        conc_params = []
+
+        conc_params.append({"name": "CO2",
+                            "diffusion coefficient": 19.1e-10,  # m^2/s
+                            "bulk": C_CO2_bulk,  # mol/m3
+                            })
+
+        conc_params.append({"name": "OH",
+                            "diffusion coefficient": 54.0e-10,  # m^2/s
+                            "bulk": C_OH_bulk,  # mol/m3
+                            })
+
+        conc_params.append({"name": "H",
+                            "diffusion coefficient": 96.0e-10,  # m^2/s
+                            "bulk": C_H_bulk,  # mol/m3
+                            })
+
+        conc_params.append({"name": "CO3",
+                            "diffusion coefficient": 7.0e-10,  # m^2/s
+                            "bulk": C_CO32_bulk,  # mol/m3
+                            })
+
+        conc_params.append({"name": "HCO3",
+                            "diffusion coefficient": 9.4e-10,  # m^2/s
+                            "bulk": C_HCO3_bulk,  # mol/m3
+                            })
+
+        physical_params = {"flow": ["diffusion"],
+                           "F": 96485.,  # C/mol
+                           }
+
+        super().__init__(conc_params, physical_params, mesh, family="CG", homog_params=homog_params)
+
+    def neumann(self, C, conc_params, u):
+        name = conc_params["name"]
+        if name in ["HCO3", "CO3", "H"]:
+            return Constant(0)
+        j = 50.  # A/m^2
+        F = self.physical_params["F"]
+        zeffCH4 = 8.
+        zeffC2H4 = 12.
+        zeffCO = 2.
+        zeffHCOO = 2.
+        zeffH2 = 2.
+        cefCH4 = 0.25
+        cefC2H4 = 0.20
+        cefCO = 0.05
+        cefHCOO = 0.10
+        cefH2 = 0.40
+        if name == "CO2":
+            print(-(j / F) * (cefHCOO / zeffHCOO + cefCO / zeffCO
+                              + cefCH4 / zeffCH4 + 2 * cefC2H4 / zeffC2H4))
+            return -(j / F) * (cefHCOO / zeffHCOO + cefCO / zeffCO
+                               + cefCH4 / zeffCH4 + 2 * cefC2H4 / zeffC2H4)
+        if name == "OH":
+            return (j / F) * (cefHCOO / zeffHCOO + 2 * cefCO / zeffCO
+                              + 8 * cefCH4 / zeffCH4 + 12 * cefC2H4 / zeffC2H4
+                              + 2 * cefH2 / zeffH2)  # note error in Gupta et al.
+
+    def set_boundary_markers(self):
+        self.boundary_markers = {"bulk dirichlet": (1,),  # U_solid = U_app, C = C_0
+                                 "neumann": (2,),
+                                 }
+
+
+solver = CarbonateSolver()
+solver.setup_solver()
+solver.solve()
+C_CO2, C_OH, C_H, C_CO3, C_HCO3 = solver.u.subfunctions
+# OH boundary layer
+x = solver.mesh.coordinates
+C_OH_bl = Function(solver.V).assign(C_OH).dat.data[100:]
+x_bl = x.dat.data[100:]
+
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib.ticker import FormatStrFormatter
+filename = "carbonate.png"
+fig = plt.figure(constrained_layout=True, figsize=(16, 8))
+spec = gridspec.GridSpec(ncols=3, nrows=2, figure=fig)
+ax1 = fig.add_subplot(spec[0, 0])
+ax2 = fig.add_subplot(spec[0, 1])
+ax3 = fig.add_subplot(spec[1, 0])
+ax4 = fig.add_subplot(spec[1, 1])
+ax5 = fig.add_subplot(spec[0, 2])
+ax6 = fig.add_subplot(spec[1, 2])
+
+plot(C_CO2, axes=ax1)
+ax1.set(xlabel='distance (m)',
+        ylabel='CO$_2$ concentration (M)')
+plot(C_HCO3, axes=ax2)
+ax2.set(xlabel='distance (m)',
+        ylabel='HCO$_3$ concentration (M)')
+plot(C_CO3, axes=ax3)
+ax3.set(xlabel='distance (m)',
+        ylabel='CO$_3$ concentration (M)')
+plot(C_OH, axes=ax4)
+ax4.set(xlabel='distance (m)',
+        ylabel='OH concentration (M)')
+plot(C_H, axes=ax5)
+ax5.set(xlabel='distance (m)',
+        ylabel='H concentration (M)')
+plt.plot(x_bl, C_OH_bl, axes=ax6, color='k', linewidth=2)
+ax6.set(xlabel='distance (m)',
+        ylabel='OH concentration (M)')
+
+plt.savefig(filename)