Basic lead acid model

CHANGELOG.md

@@ -1,8 +1,12 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
-## Optimizations
+## Features
 
 -   Made `QuickPlot` compatible with Google Colab ([#935](https://github.com/pybamm-team/PyBaMM/pull/935))
+-   Added `BasicFull` model for lead-acid ([#932](https://github.com/pybamm-team/PyBaMM/pull/932))
+
+## Optimizations
+
 -   Sped up model building ([#927](https://github.com/pybamm-team/PyBaMM/pull/927))
 -   Changed default solver for lead-acid to `CasadiSolver` ([#927](https://github.com/pybamm-team/PyBaMM/pull/927))
 

pybamm/models/full_battery_models/lead_acid/__init__.py

@@ -5,3 +5,4 @@
 from .loqs import LOQS
 from .higher_order import BaseHigherOrderModel, FOQS, Composite, CompositeExtended
 from .full import Full
+from .basic_full import BasicFull

pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py

@@ -80,16 +80,3 @@ def set_reactions(self):
             }
             self.reactions["main"]["Negative"]["s_ox"] = 0
             self.reactions["main"]["Positive"]["s_ox"] = 0
-
-    def set_soc_variables(self):
-        "Set variables relating to the state of charge."
-        # State of Charge defined as function of dimensionless electrolyte concentration
-        soc = self.variables["X-averaged electrolyte concentration"] * 100
-        self.variables.update({"State of Charge": soc, "Depth of Discharge": 100 - soc})
-
-        # Fractional charge input
-        if "Fractional Charge Input" not in self.variables:
-            fci = pybamm.Variable("Fractional Charge Input", domain="current collector")
-            self.variables["Fractional Charge Input"] = fci
-            self.rhs[fci] = -self.variables["Total current density"] * 100
-            self.initial_conditions[fci] = self.param.q_init * 100

pybamm/models/full_battery_models/lead_acid/basic_full.py

@@ -0,0 +1,286 @@
+#
+# Basic lead-acid model
+#
+import pybamm
+from .base_lead_acid_model import BaseModel
+
+
+class BasicFull(BaseModel):
+    """Porous electrode model for lead-acid, from [2]_.
+
+    This class differs from the :class:`pybamm.lead_acid.Full` model class in that it
+    shows the whole model in a single class. This comes at the cost of flexibility in
+    comparing different physical effects, and in general the main DFN class should be
+    used instead.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+
+    References
+    ----------
+    .. [2] V Sulzer, SJ Chapman, CP Please, DA Howey, and CW Monroe. Faster lead-acid
+           battery simulations from porous-electrode theory: Part II. Asymptotic
+           analysis. Journal of The Electrochemical Society 166.12 (2019), A2372–A2382..
+
+
+    **Extends:** :class:`pybamm.lead_acid.BaseModel`
+    """
+
+    def __init__(self, name="Full model"):
+        super().__init__({}, name)
+        # `param` is a class containing all the relevant parameters and functions for
+        # this model. These are purely symbolic at this stage, and will be set by the
+        # `ParameterValues` class when the model is processed.
+        param = self.param
+
+        ######################
+        # Variables
+        ######################
+        # Variables that depend on time only are created without a domain
+        Q = pybamm.Variable("Discharge capacity [A.h]")
+        # Variables that vary spatially are created with a domain
+        c_e_n = pybamm.Variable(
+            "Negative electrolyte concentration", domain="negative electrode",
+        )
+        c_e_s = pybamm.Variable(
+            "Separator electrolyte concentration", domain="separator",
+        )
+        c_e_p = pybamm.Variable(
+            "Positive electrolyte concentration", domain="positive electrode",
+        )
+        # Concatenations combine several variables into a single variable, to simplify
+        # implementing equations that hold over several domains
+        c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
+
+        # Electrolyte potential
+        phi_e_n = pybamm.Variable(
+            "Negative electrolyte potential", domain="negative electrode",
+        )
+        phi_e_s = pybamm.Variable(
+            "Separator electrolyte potential", domain="separator",
+        )
+        phi_e_p = pybamm.Variable(
+            "Positive electrolyte potential", domain="positive electrode",
+        )
+        phi_e = pybamm.Concatenation(phi_e_n, phi_e_s, phi_e_p)
+
+        # Electrode potential
+        phi_s_n = pybamm.Variable(
+            "Negative electrode potential", domain="negative electrode",
+        )
+        phi_s_p = pybamm.Variable(
+            "Positive electrode potential", domain="positive electrode",
+        )
+
+        # Porosity
+        eps_n = pybamm.Variable(
+            "Negative electrode porosity", domain="negative electrode",
+        )
+        eps_s = pybamm.Variable("Separator porosity", domain="separator")
+        eps_p = pybamm.Variable(
+            "Positive electrode porosity", domain="positive electrode",
+        )
+        eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
+
+        # Pressure (for convection)
+        pressure_n = pybamm.Variable(
+            "Negative electrolyte pressure", domain="negative electrode",
+        )
+        pressure_s = pybamm.Variable(
+            "Separator electrolyte pressure", domain="separator",
+        )
+        pressure_p = pybamm.Variable(
+            "Positive electrolyte pressure", domain="positive electrode",
+        )
+        pressure = pybamm.Concatenation(pressure_n, pressure_s, pressure_p)
+
+        # Constant temperature
+        T = param.T_init
+
+        ######################
+        # Other set-up
+        ######################
+
+        # Current density
+        i_cell = param.current_with_time
+
+        # Tortuosity
+        tor = pybamm.Concatenation(
+            eps_n ** param.b_e_n, eps_s ** param.b_e_s, eps_p ** param.b_e_p
+        )
+
+        # Interfacial reactions
+        j0_n = param.j0_n_S_ref * c_e_n
+        j_n = (
+            2
+            * j0_n
+            * pybamm.sinh(param.ne_n / 2 * (phi_s_n - phi_e_n - param.U_n(c_e_n, T)))
+        )
+        j0_p = param.j0_p_S_ref * c_e_p ** 2 * param.c_w(c_e_p)
+        j_s = pybamm.PrimaryBroadcast(0, "separator")
+        j_p = (
+            2
+            * j0_p
+            * pybamm.sinh(param.ne_p / 2 * (phi_s_p - phi_e_p - param.U_p(c_e_p, T)))
+        )
+        j = pybamm.Concatenation(j_n, j_s, j_p)
+
+        ######################
+        # State of Charge
+        ######################
+        I = param.dimensional_current_with_time
+        # The `rhs` dictionary contains differential equations, with the key being the
+        # variable in the d/dt
+        self.rhs[Q] = I * param.timescale / 3600
+        # Initial conditions must be provided for the ODEs
+        self.initial_conditions[Q] = pybamm.Scalar(0)
+
+        ######################
+        # Convection
+        ######################
+        v = -pybamm.grad(pressure)
+        l_s = pybamm.geometric_parameters.l_s
+
+        # Difference in negative and positive electrode velocities determines the
+        # velocity in the separator
+        v_box_n_right = param.beta_n * i_cell
+        v_box_p_left = param.beta_p * i_cell
+        d_vbox_s__dx = (v_box_p_left - v_box_n_right) / l_s
+
+        # Simple formula for velocity in the separator
+        dVbox_dz = pybamm.Concatenation(
+            pybamm.PrimaryBroadcast(0, "negative electrode"),
+            pybamm.PrimaryBroadcast(-d_vbox_s__dx, "separator"),
+            pybamm.PrimaryBroadcast(0, "positive electrode"),
+        )
+        beta = pybamm.Concatenation(
+            pybamm.PrimaryBroadcast(param.beta_n, "negative electrode"),
+            pybamm.PrimaryBroadcast(0, "separator"),
+            pybamm.PrimaryBroadcast(param.beta_p, "positive electrode"),
+        )
+        self.algebraic[pressure] = pybamm.div(v) + dVbox_dz - beta * j
+        self.boundary_conditions[pressure] = {
+            "left": (pybamm.Scalar(0), "Dirichlet"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.initial_conditions[pressure] = pybamm.Scalar(0)
+
+        ######################
+        # Current in the electrolyte
+        ######################
+        i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
+            param.chi(c_e) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)
+        )
+        self.algebraic[phi_e] = pybamm.div(i_e) - j
+        self.boundary_conditions[phi_e] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.initial_conditions[phi_e] = -param.U_n(param.c_e_init, param.T_init)
+
+        ######################
+        # Current in the solid
+        ######################
+        i_s_n = -param.sigma_n * (1 - eps_n) ** param.b_s_n * pybamm.grad(phi_s_n)
+        sigma_eff_p = param.sigma_p * (1 - eps_p) ** param.b_s_p
+        i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
+        # The `algebraic` dictionary contains differential equations, with the key being
+        # the main scalar variable of interest in the equation
+        self.algebraic[phi_s_n] = pybamm.div(i_s_n) + j_n
+        self.algebraic[phi_s_p] = pybamm.div(i_s_p) + j_p
+        self.boundary_conditions[phi_s_n] = {
+            "left": (pybamm.Scalar(0), "Dirichlet"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.boundary_conditions[phi_s_p] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (i_cell / pybamm.boundary_value(-sigma_eff_p, "right"), "Neumann"),
+        }
+        # Initial conditions must also be provided for algebraic equations, as an
+        # initial guess for a root-finding algorithm which calculates consistent initial
+        # conditions
+        self.initial_conditions[phi_s_n] = pybamm.Scalar(0)
+        self.initial_conditions[phi_s_p] = param.U_p(
+            param.c_e_init, param.T_init
+        ) - param.U_n(param.c_e_init, param.T_init)
+
+        ######################
+        # Porosity
+        ######################
+        beta_surf = pybamm.Concatenation(
+            pybamm.PrimaryBroadcast(param.beta_surf_n, "negative electrode"),
+            pybamm.PrimaryBroadcast(0, "separator"),
+            pybamm.PrimaryBroadcast(param.beta_surf_p, "positive electrode"),
+        )
+        deps_dt = -beta_surf * j
+        self.rhs[eps] = deps_dt
+        self.initial_conditions[eps] = param.epsilon_init
+        self.events.extend(
+            [
+                pybamm.Event(
+                    "Zero negative electrode porosity cut-off", pybamm.min(eps_n)
+                ),
+                pybamm.Event(
+                    "Max negative electrode porosity cut-off", pybamm.max(eps_n) - 1
+                ),
+                pybamm.Event(
+                    "Zero positive electrode porosity cut-off", pybamm.min(eps_p)
+                ),
+                pybamm.Event(
+                    "Max positive electrode porosity cut-off", pybamm.max(eps_p) - 1
+                ),
+            ]
+        )
+
+        ######################
+        # Electrolyte concentration
+        ######################
+        N_e = -tor * param.D_e(c_e, T) * pybamm.grad(c_e) + c_e * v
+        s = pybamm.Concatenation(
+            -pybamm.PrimaryBroadcast(param.s_plus_n_S, "negative electrode"),
+            pybamm.PrimaryBroadcast(0, "separator"),
+            -pybamm.PrimaryBroadcast(param.s_plus_p_S, "positive electrode"),
+        )
+        self.rhs[c_e] = (1 / eps) * (
+            -pybamm.div(N_e) / param.C_e
+            + (s - param.t_plus(c_e)) * j / param.gamma_e
+            - c_e * deps_dt
+        )
+        self.boundary_conditions[c_e] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.initial_conditions[c_e] = param.c_e_init
+        self.events.append(
+            pybamm.Event(
+                "Zero electrolyte concentration cut-off", pybamm.min(c_e) - 0.002
+            )
+        )
+
+        ######################
+        # (Some) variables
+        ######################
+        voltage = pybamm.boundary_value(phi_s_p, "right")
+        # The `variables` dictionary contains all variables that might be useful for
+        # visualising the solution of the model
+        pot = param.potential_scale
+
+        self.variables = {
+            "Electrolyte concentration": c_e,
+            "Current [A]": I,
+            "Negative electrode potential [V]": pot * phi_s_n,
+            "Electrolyte potential [V]": -param.U_n_ref + pot * phi_e,
+            "Positive electrode potential [V]": param.U_p_ref
+            - param.U_n_ref
+            + pot * phi_s_p,
+            "Terminal voltage [V]": param.U_p_ref - param.U_n_ref + pot * voltage,
+        }
+        self.events.extend(
+            [
+                pybamm.Event("Minimum voltage", voltage - param.voltage_low_cut),
+                pybamm.Event("Maximum voltage", voltage - param.voltage_high_cut),
+            ]
+        )
+

pybamm/models/full_battery_models/lead_acid/full.py

@@ -6,7 +6,7 @@
 
 
 class Full(BaseModel):
-    """Porous electrode model for lead-acid, from [1]_, based on the Full
+    """Porous electrode model for lead-acid, from [1]_, based on the Newman-Tiedemann
     model.
 
     Parameters

tests/integration/test_models/test_full_battery_models/test_lead_acid/test_compare_basic_models.py

@@ -0,0 +1,52 @@
+#
+# Compare basic models with full models
+#
+import pybamm
+
+import numpy as np
+import unittest
+
+
+class TestCompareBasicModels(unittest.TestCase):
+    def test_compare_full(self):
+        basic_full = pybamm.lead_acid.BasicFull()
+        full = pybamm.lead_acid.Full({"convection": True})
+
+        parameter_values = pybamm.ParameterValues(
+            chemistry=pybamm.parameter_sets.Sulzer2019
+        )
+        parameter_values["Current function [A]"] = 10
+
+        # Solve basic Full mode
+        basic_sim = pybamm.Simulation(
+            basic_full, solver=pybamm.CasadiSolver(), parameter_values=parameter_values,
+        )
+        t_eval = np.linspace(0, 400)
+        basic_sim.solve(t_eval)
+        basic_sol = basic_sim.solution
+
+        # Solve main Full model
+        sim = pybamm.Simulation(
+            full, solver=pybamm.CasadiSolver(), parameter_values=parameter_values
+        )
+        t_eval = np.linspace(0, 400)
+        sim.solve(t_eval)
+        sol = sim.solution
+
+        # Compare solution data
+        np.testing.assert_array_almost_equal(basic_sol.y, sol.y, decimal=4)
+        np.testing.assert_array_almost_equal(basic_sol.t, sol.t, decimal=4)
+        # Compare variables
+        for name in basic_full.variables:
+            np.testing.assert_array_almost_equal(
+                basic_sol[name].entries, sol[name].entries, decimal=4
+            )
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()