added example to be used as tutorial

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)]

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)')
 

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()