From 4e8e30fe92fce9c216f0b2944f5db20acf2067db Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Fri, 5 Apr 2024 13:01:43 -0700 Subject: [PATCH 1/4] 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 From 5a94ff314814d4ccd53ba1c2ec529e2428bc0d97 Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Fri, 12 Apr 2024 10:35:40 -0700 Subject: [PATCH 2/4] added example to be used as tutorial --- echemfem/solver.py | 5 +- examples/gupta.py | 2 +- examples/gupta_advection_migration.py | 143 ++++++++++++++++++++++++++ 3 files changed, 147 insertions(+), 3 deletions(-) create mode 100644 examples/gupta_advection_migration.py diff --git a/echemfem/solver.py b/echemfem/solver.py index 1146d25..b879928 100644 --- a/echemfem/solver.py +++ b/echemfem/solver.py @@ -368,7 +368,7 @@ def _internal_dS(subdomain_id=None, domain=None): us[self.num_liquid:self.num_liquid + self.num_gas], self.pg, gas_params) # setup applied voltage - if self.flow["poisson"] or self.flow["electroneutrality"] or self.flow["electroneutrality full"]: + if physical_params.get("U_app"): U_app = physical_params["U_app"] if isinstance(U_app, (Constant, Function)): self.U_app = U_app @@ -636,7 +636,8 @@ def setup_solver(self, initial_guess=True, initial_solve=True): elif (self.flow["electroneutrality"] or self.flow["poisson"] or self.flow["electroneutrality full"]): U0 = Function(self.Vu) - U0.assign(self.U_app / 2) + if hasattr(self, "U_app"): + U0.assign(self.U_app / 2) if initial_solve: v0 = TestFunction(self.Vu) u0 = [us[i] for i in range(self.num_mass)] diff --git a/examples/gupta.py b/examples/gupta.py index 90510c3..e7aba4c 100644 --- a/examples/gupta.py +++ b/examples/gupta.py @@ -124,7 +124,7 @@ def set_boundary_markers(self): plot(C_OH, axes=ax4) ax4.set(xlabel='distance (m)', ylabel='OH concentration (M)') -plt.plot(x_bl, C_OH_bl, axes=ax5, color='k', linewidth=2) +plt.plot(x_bl, C_OH_bl, color='k', linewidth=2) ax5.set(xlabel='distance (m)', ylabel='OH concentration (M)') diff --git a/examples/gupta_advection_migration.py b/examples/gupta_advection_migration.py new file mode 100644 index 0000000..b009c8a --- /dev/null +++ b/examples/gupta_advection_migration.py @@ -0,0 +1,143 @@ +from firedrake import * +from echemfem import EchemSolver, RectangleBoundaryLayerMesh + + +class GuptaSolver(EchemSolver): + """ + Flow past a flat plate electrode for CO2 reduction. + Using electroneutral Nernst-Planck. + The homogenous bulk reactions and the constant-rate charge-transfer + kinetics are taken 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. + """ + + def __init__(self): + + Ly = 1e-3 + Lx = 5e-3 + + mesh = RectangleBoundaryLayerMesh(100, 50, Lx, Ly, 50, 1e-6, boundary=(3,)) + x, y = SpatialCoordinate(mesh) + active = conditional(And(x >= 1e-3, x < Lx-1e-3), 1., 0.) + + conc_params = [] + + conc_params.append({"name": "CO2", + "diffusion coefficient": 19.1e-10, # m^2/s + "bulk": 34.2, # mol/m3 + "z": 0, + }) + + conc_params.append({"name": "HCO3", + "diffusion coefficient": 9.23e-10, # m^2/s + "bulk": 499., # mol/m3 + "z": -1, + }) + + conc_params.append({"name": "CO3", + "diffusion coefficient": 11.9e-10, # m^2/s + "bulk": 7.6e-1, # mol/m3 + "z": -2, + }) + + conc_params.append({"name": "OH", + "diffusion coefficient": 52.7e-10, # m^2/s + "bulk": 3.3e-4, # mol/m3 + "z": -1, + }) + + conc_params.append({"name": "K", + "diffusion coefficient": 19.6E-10, # m^2/s + "bulk": 499. + 7.6e-1 + 3.3e-4, # mol/m3 + "z": 1, + }) + + homog_params = [] + + homog_params.append({"stoichiometry": {"CO2": -1, + "OH": -1, + "HCO3": 1, + }, + "forward rate constant": 5.93, + "backward rate constant": 1.34e-4 + }) + + homog_params.append({"stoichiometry": {"HCO3": -1, + "OH": -1, + "CO3": 1, + }, + "forward rate constant": 1e5, + "backward rate constant": 2.15e4 + }) + + physical_params = {"flow": ["diffusion", "electroneutrality", "migration", "advection"], + "F": 96485.3329, # C/mol + "R": 8.3144598, # J/K/mol + "T": 273.15 + 25., # K + } + + def current(cef): + j = 50. + def curr(u): + return cef * j * active + return curr + + echem_params = [] + + echem_params.append({"reaction": current(0.1), # HCOO + "stoichiometry": {"CO2": -1, + "OH": 1 + }, + "electrons": 2, + "boundary": "electrode", + }) + + echem_params.append({"reaction": current(0.05), # CO + "stoichiometry": {"CO2": -1, + "OH": 2 + }, + "electrons": 2, + "boundary": "electrode", + }) + + echem_params.append({"reaction": current(0.25), # CH4 + "stoichiometry": {"CO2": -1, + "OH": 8 + }, + "electrons": 8, + "boundary": "electrode", + }) + + echem_params.append({"reaction": current(0.2), # C2H4 + "stoichiometry": {"CO2": -2, + "OH": 12 + }, + "electrons": 12, + "boundary": "electrode", + }) + + echem_params.append({"reaction": current(0.4), # H2 + "stoichiometry": {"OH": 2 + }, + "electrons": 2, + "boundary": "electrode", + }) + + super().__init__(conc_params, physical_params, mesh, family="CG", echem_params=echem_params, homog_params=homog_params) + + def set_boundary_markers(self): + self.boundary_markers = {"inlet": (1,), + "outlet": (2,), + "bulk": (4,), + "electrode": (3,), + } + def set_velocity(self): + _, y = SpatialCoordinate(self.mesh) + self.vel = as_vector([1.91*y, Constant(0)]) # m/s + +solver = GuptaSolver() +solver.setup_solver() +solver.solve() From d3dd3f5f1ffbeb82a92cfdd5bb32fb2fbc0dd711 Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Wed, 17 Apr 2024 15:16:23 -0700 Subject: [PATCH 3/4] Adding short description for examples --- examples/BMCSL.py | 5 +++ examples/bicarb.py | 3 +- examples/bicarb_Ag_Cu_tandem_example.py | 27 ++++++++-------- examples/bortels_threeion.py | 14 +++++++++ .../bortels_threeion_extruded_3D_nondim.py | 20 +++++++++++- examples/bortels_threeion_nondim.py | 19 ++++++++++++ examples/bortels_twoion.py | 14 +++++++++ examples/bortels_twoion_nondim.py | 19 ++++++++++++ examples/carbonate.py | 31 ++++++++++++------- examples/carbonate_homog_params.py | 31 ++++++++++++------- examples/catalyst_layer_Cu.py | 8 ++++- examples/catalyst_layer_Cu_full.py | 10 +++++- examples/cylindrical_pore.py | 8 +++-- examples/gupta.py | 6 ++++ examples/gupta_advection_migration.py | 7 ++--- examples/mms_scaling.py | 10 ++++++ examples/paper_test.py | 5 ++- examples/simple_flow_battery.py | 18 ++++++++--- examples/tworxn.py | 15 ++++----- examples/tworxn_irregular.py | 22 ++++++++----- 20 files changed, 223 insertions(+), 69 deletions(-) diff --git a/examples/BMCSL.py b/examples/BMCSL.py index 66d5e8d..0714e11 100644 --- a/examples/BMCSL.py +++ b/examples/BMCSL.py @@ -3,7 +3,12 @@ @author: Nitish Govindarajan +Simple example for the BMCSL model for finite size ion effects + Code to reproduce Figure 1 from doi:10.1016/j.jcis.2007.08.006 +Biesheuvel, P.M. and Van Soestbergen, M., 2007. Counterion volume effects in +mixed electrical double layers. Journal of Colloid and Interface Science, +316(2), pp.490-499. """ # Import required libraries diff --git a/examples/bicarb.py b/examples/bicarb.py index b2dae9d..402d798 100644 --- a/examples/bicarb.py +++ b/examples/bicarb.py @@ -5,7 +5,8 @@ class CarbonateSolver(EchemSolver): def __init__(self): """ - Bicarbonate example reproduced from: + Flow past an electrode with bicarbonate bulk reactions and linear CO2 electrolysis + Example reproduced from: Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow Reactors. Industrial & Engineering Chemistry Research, 60(31), diff --git a/examples/bicarb_Ag_Cu_tandem_example.py b/examples/bicarb_Ag_Cu_tandem_example.py index 7e618f9..212d611 100644 --- a/examples/bicarb_Ag_Cu_tandem_example.py +++ b/examples/bicarb_Ag_Cu_tandem_example.py @@ -11,19 +11,19 @@ comm = MPI.COMM_WORLD rank = comm.Get_rank() -##Description: -#This example case runs a tandem electrode system, with an Ag catalyst placed -#in front of a Cu catalyst, at a Peclet number of 1E3. -#The advection/diffusion/migration/electroneutrality equations are solved -#for the bulk electrolyte, and a Tafel description is applied for the surface -#catalyst reactions. +""" +This example case runs a tandem electrode system, with an Ag catalyst placed +in front of a Cu catalyst, at a Peclet number of 1E3. +The advection/diffusion/migration/electroneutrality equations are solved +for the bulk electrolyte, and a Tafel description is applied for the surface +catalyst reactions. -#Assumed products of Cu catalyst are C2H4, C2H6O, and H2. -#Assumed products of Ag catalysis are CO and H2. +Assumed products of Cu catalyst are C2H4, C2H6O, and H2. +Assumed products of Ag catalysis are CO and H2. + +-------------------------------------------------------------------------- +list of references -##-------------------------------------------------------------------------- -#list of references -""" Bicarbonate flow reactor setup given in: Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow @@ -53,10 +53,9 @@ Microkinetics with Continuum Transport Models to Understand Electrochemical CO2 Reduction in Flow Reactors. PRX Energy, 2(3), pp.033010. +-------------------------------------------------------------------------- +end of list of references """ -##-------------------------------------------------------------------------- -#end of list of references - #--------------------------------------------------------------------------- diff --git a/examples/bortels_threeion.py b/examples/bortels_threeion.py index 17da604..15acdad 100644 --- a/examples/bortels_threeion.py +++ b/examples/bortels_threeion.py @@ -7,6 +7,20 @@ parser.add_argument("--family", type=str, default='DG') args, _ = parser.parse_known_args() +""" +2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh +mesh with refinement close to the electrodes. GMSH is used to get the mesh file +from the .geo file as follows: + + gmsh -2 bortels_structuredquad.geo + +Three ion case taken from: +Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional +upwinding method as a new simulation tool for the analysis of multi-ion +electrolytes controlled by diffusion, convection and migration. Part 1. Steady +state analysis of a parallel plane flow channel. Journal of Electroanalytical +Chemistry, 404(1), pp.15-26. +""" class BortelsSolver(EchemSolver): def __init__(self): diff --git a/examples/bortels_threeion_extruded_3D_nondim.py b/examples/bortels_threeion_extruded_3D_nondim.py index 8d9acb6..760afce 100644 --- a/examples/bortels_threeion_extruded_3D_nondim.py +++ b/examples/bortels_threeion_extruded_3D_nondim.py @@ -8,6 +8,24 @@ set_log_level(DEBUG) +""" +3D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh +mesh with refinement close to the electrodes. The mesh is then extruded in the +z-direction. Custom block preconditioners are used. The file submitter.pl can +be adapted to launch slurm jobs using perl. + +Three ion case adapted from the 2D case in: +Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional +upwinding method as a new simulation tool for the analysis of multi-ion +electrolytes controlled by diffusion, convection and migration. Part 1. Steady +state analysis of a parallel plane flow channel. Journal of Electroanalytical +Chemistry, 404(1), pp.15-26. + +Nondimensionalization and preconditioning from: +Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the +electroneutral Nernst-Planck equations. Journal of Computational Physics, 475, +p.111859. +""" parser = argparse.ArgumentParser() parser.add_argument('-ref_levels', type=int, default=0) @@ -55,6 +73,7 @@ def coarsen_form(self, form, args): class BortelsSolver(EchemSolver): def __init__(self): # Create an initial coarse mesh + # Note: will need to use "gmsh -2 mesh_name.geo" command to convert to .msh format if args.mesh == 'structured': plane_mesh = Mesh('bortels_structuredquad_nondim_coarse1664.msh') layers = 8 @@ -197,7 +216,6 @@ def set_velocity(self): # Custom solver parameters - class CoarsenPenaltyPMGPC(P1PC): def coarsen_form(self, form, replace_dict): k = 1 diff --git a/examples/bortels_threeion_nondim.py b/examples/bortels_threeion_nondim.py index 797a661..266e8c6 100644 --- a/examples/bortels_threeion_nondim.py +++ b/examples/bortels_threeion_nondim.py @@ -1,6 +1,25 @@ from firedrake import * from echemflow import EchemSolver +""" +2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh +mesh with refinement close to the electrodes. GMSH is used to get the mesh file +from the .geo file as follows: + + gmsh -2 bortels_structuredquad_nondim.geo + +Three ion case taken from: +Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional +upwinding method as a new simulation tool for the analysis of multi-ion +electrolytes controlled by diffusion, convection and migration. Part 1. Steady +state analysis of a parallel plane flow channel. Journal of Electroanalytical +Chemistry, 404(1), pp.15-26. + +Nondimensionalization from: +Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the +electroneutral Nernst-Planck equations. Journal of Computational Physics, 475, +p.111859. +""" class BortelsSolver(EchemSolver): def __init__(self): diff --git a/examples/bortels_twoion.py b/examples/bortels_twoion.py index cb10c44..d65ae4e 100644 --- a/examples/bortels_twoion.py +++ b/examples/bortels_twoion.py @@ -5,6 +5,20 @@ parser.add_argument("--family", type=str, default='DG') args, _ = parser.parse_known_args() +""" +2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom GMSH +mesh with refinement close to the electrodes. GMSH is used to get the mesh file +from the .geo file as follows: + + gmsh -2 bortels_structuredquad.geo + +Two ion case taken from: +Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional +upwinding method as a new simulation tool for the analysis of multi-ion +electrolytes controlled by diffusion, convection and migration. Part 1. Steady +state analysis of a parallel plane flow channel. Journal of Electroanalytical +Chemistry, 404(1), pp.15-26. +""" class BortelsSolver(EchemSolver): def __init__(self): diff --git a/examples/bortels_twoion_nondim.py b/examples/bortels_twoion_nondim.py index 971dc1f..4687b4d 100644 --- a/examples/bortels_twoion_nondim.py +++ b/examples/bortels_twoion_nondim.py @@ -1,6 +1,25 @@ from firedrake import * from echemfem import EchemSolver +""" +2D Flow-plate reactor with electroneutral Nernst-Planck. Using a custom gmsh +mesh with refinement close to the electrodes. GMSH is used to get the mesh file +from the .geo file as follows: + + gmsh -2 bortels_structuredquad_nondim.geo + +Two ion case taken from: +Bortels, L., Deconinck, J. and Van Den Bossche, B., 1996. The multi-dimensional +upwinding method as a new simulation tool for the analysis of multi-ion +electrolytes controlled by diffusion, convection and migration. Part 1. Steady +state analysis of a parallel plane flow channel. Journal of Electroanalytical +Chemistry, 404(1), pp.15-26. + +Nondimensionalization from: +Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the +electroneutral Nernst-Planck equations. Journal of Computational Physics, 475, +p.111859. +""" class BortelsSolver(EchemSolver): def __init__(self): diff --git a/examples/carbonate.py b/examples/carbonate.py index 2ac3721..7f5e030 100644 --- a/examples/carbonate.py +++ b/examples/carbonate.py @@ -1,20 +1,24 @@ from firedrake import * from echemfem import EchemSolver, IntervalBoundaryLayerMesh +""" +A 1D example of diffusion-reaction for CO2 electrolysis with bicarbonate bulk +reactions. The bulk reactions and charge-transfer reactions are implemented by +hand. + +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. + +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. +""" 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 @@ -115,6 +119,9 @@ def set_boundary_markers(self): solver = CarbonateSolver() solver.setup_solver() solver.solve() + +### Plotting + C_CO2, C_OH, C_H, C_CO3, C_HCO3 = solver.u.subfunctions # OH boundary layer x = solver.mesh.coordinates diff --git a/examples/carbonate_homog_params.py b/examples/carbonate_homog_params.py index aff74f0..a08d4ee 100644 --- a/examples/carbonate_homog_params.py +++ b/examples/carbonate_homog_params.py @@ -1,20 +1,24 @@ from firedrake import * from echemfem import EchemSolver, IntervalBoundaryLayerMesh +""" +A 1D example of diffusion-reaction for CO2 electrolysis with bicarbonate bulk +reactions. The charge-transfer reactions are implemented by hand but the bulk +reactions are implemented using the homog_params interface. + +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. + +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. +""" 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 @@ -142,6 +146,9 @@ def set_boundary_markers(self): solver = CarbonateSolver() solver.setup_solver() solver.solve() + +### Plotting + C_CO2, C_OH, C_H, C_CO3, C_HCO3 = solver.u.subfunctions # OH boundary layer x = solver.mesh.coordinates diff --git a/examples/catalyst_layer_Cu.py b/examples/catalyst_layer_Cu.py index 46a75a8..33c872c 100644 --- a/examples/catalyst_layer_Cu.py +++ b/examples/catalyst_layer_Cu.py @@ -6,7 +6,10 @@ import csv import pdb """ -Model fromn +1D model for the CO2 electrolysis in a copper catalyst layer of a gas diffusion +electrode (GDE). The model uses electroneutral Nernst-Planck in a porous +medium. + Corral D, Lee DU, Ehlinger VM, Nitopi S, Acosta JE, Wang L, King AJ, Feaster JT, Lin YR, Weber AZ, Baker SE. Bridging knowledge gaps in liquid-and vapor-fed CO2 electrolysis through active electrode area. Chem Catalysis. 2022 Oct 12. @@ -240,6 +243,9 @@ def set_boundary_markers(self): solver.U_app.assign(Vs) print("V = %d mV" % (Vs * 1000)) solver.solve() + + ## Plotting + cCO2, cOH, cH, cCO3, cHCO3, phi2, phi1 = solver.u.subfunctions cK = Function(solver.V).assign(2 * cCO3 + cOH + cHCO3 - cH) diff --git a/examples/catalyst_layer_Cu_full.py b/examples/catalyst_layer_Cu_full.py index bcf77ae..a4c0e46 100644 --- a/examples/catalyst_layer_Cu_full.py +++ b/examples/catalyst_layer_Cu_full.py @@ -5,7 +5,12 @@ import numpy as np import csv """ -Model fromn +1D model for the CO2 electrolysis in a copper catalyst layer of a gas diffusion +electrode (GDE). The model uses electroneutral Nernst-Planck in a porous +medium. This version of the code does not eliminate a species by substituting +the electroneutrality equation, but rather explicitly solves for +electroneutrality. + Corral D, Lee DU, Ehlinger VM, Nitopi S, Acosta JE, Wang L, King AJ, Feaster JT, Lin YR, Weber AZ, Baker SE. Bridging knowledge gaps in liquid-and vapor-fed CO2 electrolysis through active electrode area. Chem Catalysis. 2022 Oct 12. @@ -240,6 +245,9 @@ def set_boundary_markers(self): solver.U_app.assign(Vs) print("V = %d mV" % (Vs * 1000)) solver.solve() + + ## Plotting + cCO2, cOH, cH, cCO3, cHCO3, cK, phi2, phi1 = solver.u.subfunctions x_ = Function(V).interpolate(solver.mesh.coordinates[0]).vector().dat.data diff --git a/examples/cylindrical_pore.py b/examples/cylindrical_pore.py index 742bab6..0d9504e 100644 --- a/examples/cylindrical_pore.py +++ b/examples/cylindrical_pore.py @@ -2,7 +2,9 @@ from echemfem import EchemSolver import numpy as np """ -Model from +A symmetric cylindrical pore model for CO2 electrolysis using electroneutral +Nernst-Planck and simplified bicarbonate bulk reactions. + Elucidating Mass Transport Regimes in Gas Diffusion Electrodes for CO2 Electroreduction Thomas Moore, Xiaoxing Xia, Sarah E. Baker, Eric B. Duoss, and Victor A. Beck ACS Energy Letters 2021 6 (10), 3600-3606 @@ -170,18 +172,18 @@ def set_boundary_markers(self): solver.save_solutions = False solver.solve() +## Plotting + # getting active region _, Z = SpatialCoordinate(solver.mesh) active = conditional(le(Z, 5e-6), 1., 0.) - def get_boundary_dofs(V, i): u = Function(V) bc = DirichletBC(V, active, i) bc.apply(u) return np.where(u.vector()[:] == 1) - dofs = get_boundary_dofs(solver.V, 2) Z_cat = Function(solver.V).interpolate(Z).dat.data[dofs] diff --git a/examples/gupta.py b/examples/gupta.py index e7aba4c..ca4bcdb 100644 --- a/examples/gupta.py +++ b/examples/gupta.py @@ -4,6 +4,9 @@ class GuptaSolver(EchemSolver): """ + A 1D example of diffusion-reaction for CO2 electrolysis with bicarbonate bulk + reactions. + 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 @@ -94,6 +97,9 @@ def set_boundary_markers(self): solver = GuptaSolver() solver.setup_solver() solver.solve() + +## Plotting + C_CO2, C_HCO3, C_CO3, C_OH = solver.u.subfunctions # OH boundary layer x = solver.mesh.coordinates diff --git a/examples/gupta_advection_migration.py b/examples/gupta_advection_migration.py index b009c8a..bf0a8b8 100644 --- a/examples/gupta_advection_migration.py +++ b/examples/gupta_advection_migration.py @@ -4,10 +4,9 @@ class GuptaSolver(EchemSolver): """ - Flow past a flat plate electrode for CO2 reduction. - Using electroneutral Nernst-Planck. - The homogenous bulk reactions and the constant-rate charge-transfer - kinetics are taken from: + Flow past a flat plate electrode for CO2 reduction using electroneutral + Nernst-Planck. The homogenous bulk reactions and the constant-rate + charge-transfer kinetics are taken from: Gupta, N., Gattrell, M. and MacDougall, B., 2006. Calculation for the cathode surface concentrations in the electrochemical reduction of CO2 in diff --git a/examples/mms_scaling.py b/examples/mms_scaling.py index d71c92c..cbac206 100644 --- a/examples/mms_scaling.py +++ b/examples/mms_scaling.py @@ -8,6 +8,16 @@ set_log_level(DEBUG) +""" +A manufactured example of 3D electroneutral Nernst-Planck to evaluate the +convergence of the DG scheme and scaling of custom block preconditioners. The +file submitter_mms.pl can be adapted to launch slurm jobs using perl. + +DG scheme and preconditioning from: +Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the +electroneutral Nernst-Planck equations. Journal of Computational Physics, 475, +p.111859. +""" parser = argparse.ArgumentParser() parser.add_argument('-ref_levels', type=int, default=0) diff --git a/examples/paper_test.py b/examples/paper_test.py index 5e782eb..779e2ea 100644 --- a/examples/paper_test.py +++ b/examples/paper_test.py @@ -4,7 +4,10 @@ from mpi4py import MPI from petsc4py import PETSc """ -Convergence tests for Roy et al, 2021 +Convergence tests for the DG scheme in +Roy, T., Andrej, J. and Beck, V.A., 2023. A scalable DG solver for the +electroneutral Nernst-Planck equations. Journal of Computational Physics, 475, +p.111859. """ PETSc.Sys.popErrorHandler() diff --git a/examples/simple_flow_battery.py b/examples/simple_flow_battery.py index d91b507..482bc79 100644 --- a/examples/simple_flow_battery.py +++ b/examples/simple_flow_battery.py @@ -1,6 +1,19 @@ from firedrake import * from echemfem import EchemSolver, NavierStokesBrinkmanFlowSolver import argparse + + +""" +A simple Vanadium flow battery using advection-diffusion-reaction, Poisson for +the ionic potential with a predefinied conductivity and Stokes-Brinkman for the +flow. + +Model taken from +Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2022. Topology optimization +of 3D flow fields for flow batteries. Journal of The Electrochemical Society, +169(5), p.050540. +""" + parser = argparse.ArgumentParser(add_help=False) parser.add_argument("--family", type=str, default='CG') parser.add_argument("--vel_file", type=str, default=None) @@ -17,11 +30,6 @@ mesh = RectangleMesh(50, 50, electrode_length, electrode_thickness, quadrilateral=True) class FlowBatterySolver(EchemSolver): - """ - Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2022. Topology optimization - of 3D flow fields for flow batteries. Journal of The Electrochemical Society, - 169(5), p.050540. - """ def __init__(self): conc_params = [] diff --git a/examples/tworxn.py b/examples/tworxn.py index e1c109f..ad155d7 100644 --- a/examples/tworxn.py +++ b/examples/tworxn.py @@ -1,16 +1,17 @@ from firedrake import * from echemfem import EchemSolver, RectangleBoundaryLayerMesh +""" +A 2D flow past the electrode toy model with two species and +advection-diffusion. Taken from +Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of +the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow +Reactors. Industrial & Engineering Chemistry Research, 60(31), +pp.11824-11833. +""" class CarbonateSolver(EchemSolver): def __init__(self): - """ - Two reactions example reproduced from: - Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of - the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow - Reactors. Industrial & Engineering Chemistry Research, 60(31), - pp.11824-11833. - """ Ly = 0.1 Lx = 1. diff --git a/examples/tworxn_irregular.py b/examples/tworxn_irregular.py index 7fc18a8..3a50578 100644 --- a/examples/tworxn_irregular.py +++ b/examples/tworxn_irregular.py @@ -1,6 +1,21 @@ from firedrake import * from echemfem import EchemSolver, NavierStokesFlowSolver +""" +A 2D flow past an irregular electrode toy model with two species and +advection-diffusion, and Navier-Stokes for the flow. The electrode surface +consists of square blocks. The mesh can be create in GMSH using the following +command: + + gmsh -2 squares_small.geo + +The model is adapated from +Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of +the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow +Reactors. Industrial & Engineering Chemistry Research, 60(31), +pp.11824-11833. +""" + peclet=10 damkohler=10 Ly = 0.1 @@ -11,13 +26,6 @@ class CarbonateSolver(EchemSolver): def __init__(self): - """ - Two reactions example reproduced from: - Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2021. Analysis of - the Reactive CO2 Surface Flux in Electrocatalytic Aqueous Flow - Reactors. Industrial & Engineering Chemistry Research, 60(31), - pp.11824-11833. - """ C_1_inf = 1. C_2_inf = Constant(0) From 9faf64be3b223c5d5d1aa89e2de1124f38ed4c89 Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Thu, 18 Apr 2024 13:41:21 -0700 Subject: [PATCH 4/4] Add index of examples --- docs/source/examples.rst | 37 +++++++++++++++++++++++++++++++++ docs/source/index.rst | 1 + examples/gupta.py | 4 ++-- examples/simple_flow_battery.py | 4 ++-- 4 files changed, 42 insertions(+), 4 deletions(-) create mode 100644 docs/source/examples.rst diff --git a/docs/source/examples.rst b/docs/source/examples.rst new file mode 100644 index 0000000..ee02192 --- /dev/null +++ b/docs/source/examples.rst @@ -0,0 +1,37 @@ +Examples +======== + +Here is a brief overview of the examples in `echemfem/examples `_. + +* **Planar flow reactor** with electroneutral Nernst-Planck. A Butler-Volmer expression is used for the redox reaction, and there are no homogeneous bulk reactions. A custom mesh is used with refinements close to the electrodes. + + * `2D case for binary electrolyte `_ (:math:`\ce{CuSO4}`) + * `Nondimensional version of binary case `_ + * `2D case for three-ion electrolyte `_ (:math:`\ce{CuSO4}` and :math:`\ce{H2SO4}`) + * `Nondimensional version of three ion case `_ + * `3D version of the three-ion case using custom preconditioners for a Discontinuous Galerkin (DG) scheme, and nondimensional `_ + +* Simple 1D reaction-diffusion system for :math:`\ce{CO2}` electrolysis in :math:`\ce{KHCO3}` + + * `Simplified bicarbonate bulk reactions `_ + * `Full bicarbonate bulk reactions `_ + * `Different implementation of the full bicarbonate bulk reactions `_ (using the :doc:`homog_params ` interface) + +* 2D flow-past the electrode with advection-diffusion + + * `Two species toy model with shear flow `_ + * `CO2 electrolysis in bicarbonate with linear charge-transfer kinetics and shear flow `_ + * `Two species toy model with irregular electrode with Navier-Stokes for the flow `_ + +* 1D model for the :math:`\ce{CO2}` electrolysis in a copper catalyst layer of a gas diffusion electrode (GDE). The model uses electroneutral Nernst-Planck in a porous medium. + + * `Using substitution of the electroneutrality equation to eliminate a species `_ + * `Solving the electroneutrality equation explicitly `_ + +* `A simple Vanadium flow battery using advection-diffusion-reaction, Poisson for the ionic potential with a predefinied conductivity and Navier-Stokes-Brinkman for the flow `_ + +* `A symmetric cylindrical pore model for CO2 electrolysis using electroneutral Nernst-Planck and simplified bicarbonate bulk reactions `_ + +* `A tandem flow-past the electrode system with an Ag catalyst placed in front of a Cu catalyst with electroneutral Nernst-Planck and Tafel kinetics `_ + +* `Simple example for the BMCSL model for finite size ion effects `_ diff --git a/docs/source/index.rst b/docs/source/index.rst index e3040c1..c396b5a 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -24,6 +24,7 @@ Contents .. toctree:: :titlesonly: + examples api Copyright and License diff --git a/examples/gupta.py b/examples/gupta.py index ca4bcdb..13320b5 100644 --- a/examples/gupta.py +++ b/examples/gupta.py @@ -4,8 +4,8 @@ class GuptaSolver(EchemSolver): """ - A 1D example of diffusion-reaction for CO2 electrolysis with bicarbonate bulk - reactions. + A 1D example of diffusion-reaction for CO2 electrolysis with simplified + bicarbonate bulk reactions. Steady-state version of example from Gupta, N., Gattrell, M. and MacDougall, B., 2006. Calculation for the diff --git a/examples/simple_flow_battery.py b/examples/simple_flow_battery.py index 482bc79..cd1e79e 100644 --- a/examples/simple_flow_battery.py +++ b/examples/simple_flow_battery.py @@ -5,8 +5,8 @@ """ A simple Vanadium flow battery using advection-diffusion-reaction, Poisson for -the ionic potential with a predefinied conductivity and Stokes-Brinkman for the -flow. +the ionic potential with a predefinied conductivity and Navier-Stokes-Brinkman +for the flow. Model taken from Lin, T.Y., Baker, S.E., Duoss, E.B. and Beck, V.A., 2022. Topology optimization