diff --git a/docs/source/user_guide/homog_params.rst b/docs/source/user_guide/homog_params.rst new file mode 100644 index 0000000..0f9265a --- /dev/null +++ b/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. diff --git a/docs/source/user_guide/index.rst b/docs/source/user_guide/index.rst index 545dc3e..68f3a8a 100644 --- a/docs/source/user_guide/index.rst +++ b/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 diff --git a/docs/source/user_guide/physical_params.rst b/docs/source/user_guide/physical_params.rst index 8a723ff..d7684b8 100644 --- a/docs/source/user_guide/physical_params.rst +++ b/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. diff --git a/echemfem/solver.py b/echemfem/solver.py index b6d12f1..ac35ddd 100644 --- a/echemfem/solver.py +++ b/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 """ diff --git a/examples/carbonate_homog_params.py b/examples/carbonate_homog_params.py new file mode 100644 index 0000000..aff74f0 --- /dev/null +++ b/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)