Add "voltage as a state" option

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Adds an option "voltage as a state" that can be "false" (default) or "true". If "true" adds an explicit algebraic equation for the voltage. ([#4507](https://github.com/pybamm-team/PyBaMM/pull/4507))
 - Improved `QuickPlot` accuracy for simulations with Hermite interpolation. ([#4483](https://github.com/pybamm-team/PyBaMM/pull/4483))
 - Added Hermite interpolation to the (`IDAKLUSolver`) that improves the accuracy and performance of post-processing variables. ([#4464](https://github.com/pybamm-team/PyBaMM/pull/4464))
 - Added `BasicDFN` model for sodium-ion batteries ([#4451](https://github.com/pybamm-team/PyBaMM/pull/4451))

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -210,6 +210,9 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 solve an algebraic equation for it. Default is "false", unless "SEI film
                 resistance" is distributed in which case it is automatically set to
                 "true".
+            * "voltage as a state" : str
+                Whether to make a state for the voltage and solve an algebraic equation
+                for it. Default is "false".
             * "working electrode" : str
                 Can be "both" (default) for a standard battery or "positive" for a
                 half-cell where the negative electrode is replaced with a lithium metal
@@ -321,6 +324,7 @@ def __init__(self, extra_options):
                 "heterogeneous catalyst",
                 "cation-exchange membrane",
             ],
+            "voltage as a state": ["false", "true"],
             "working electrode": ["both", "positive"],
             "x-average side reactions": ["false", "true"],
         }

src/pybamm/models/submodels/electrode/base_electrode.py

@@ -119,7 +119,7 @@ def _get_standard_current_collector_potential_variables(
         V_cc = phi_s_cp - phi_s_cn
 
         # Voltage
-        # Note phi_s_cn is always zero at the negative tab
+        # Note phi_s_cn is always zero at the negative tab by definition
         V = pybamm.boundary_value(phi_s_cp, "positive tab")
 
         # Voltage is local current collector potential difference at the tabs, in 1D
@@ -128,10 +128,12 @@ def _get_standard_current_collector_potential_variables(
             "Negative current collector potential [V]": phi_s_cn,
             "Positive current collector potential [V]": phi_s_cp,
             "Local voltage [V]": V_cc,
+            "Voltage expression [V]": V - delta_phi_contact,
             "Terminal voltage [V]": V - delta_phi_contact,
-            "Voltage [V]": V - delta_phi_contact,
             "Contact overpotential [V]": delta_phi_contact,
         }
+        if self.options["voltage as a state"] == "false":
+            variables.update({"Voltage [V]": V - delta_phi_contact})
 
         return variables
 

src/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py

@@ -19,9 +19,23 @@ def get_fundamental_variables(self):
             "Current [A]": I,
             "C-rate": I / self.param.Q,
         }
+        if self.options.get("voltage as a state") == "true":
+            V = pybamm.Variable("Voltage [V]")
+            variables.update({"Voltage [V]": V})
 
         return variables
 
+    def set_initial_conditions(self, variables):
+        if self.options.get("voltage as a state") == "true":
+            V = variables["Voltage [V]"]
+            self.initial_conditions[V] = self.param.ocv_init
+
+    def set_algebraic(self, variables):
+        if self.options.get("voltage as a state") == "true":
+            V = variables["Voltage [V]"]
+            V_expression = variables["Voltage expression [V]"]
+            self.algebraic[V] = V - V_expression
+
 
 class ExplicitPowerControl(BaseModel):
     """External circuit with current set explicitly to hit target power."""

src/pybamm/models/submodels/external_circuit/function_control_external_circuit.py

@@ -48,13 +48,19 @@ def get_fundamental_variables(self):
             "Current [A]": I,
             "C-rate": I / self.param.Q,
         }
+        if self.options.get("voltage as a state") == "true":
+            V = pybamm.Variable("Voltage [V]")
+            variables.update({"Voltage [V]": V})
 
         return variables
 
     def set_initial_conditions(self, variables):
         # Initial condition as a guess for consistent initial conditions
         i_cell = variables["Current variable [A]"]
         self.initial_conditions[i_cell] = self.param.Q
+        if self.options.get("voltage as a state") == "true":
+            V = variables["Voltage [V]"]
+            self.initial_conditions[V] = self.param.ocv_init
 
     def set_rhs(self, variables):
         # External circuit submodels are always equations on the current
@@ -71,6 +77,10 @@ def set_algebraic(self, variables):
         if self.control == "algebraic":
             i_cell = variables["Current variable [A]"]
             self.algebraic[i_cell] = self.external_circuit_function(variables)
+        if self.options.get("voltage as a state") == "true":
+            V = variables["Voltage [V]"]
+            V_expression = variables["Voltage expression [V]"]
+            self.algebraic[V] = V - V_expression
 
 
 class VoltageFunctionControl(FunctionControl):

tests/unit/test_models/test_full_battery_models/test_base_battery_model.py

@@ -51,6 +51,7 @@
 'thermal': 'x-full' (possible: ['isothermal', 'lumped', 'x-lumped', 'x-full'])
 'total interfacial current density as a state': 'false' (possible: ['false', 'true'])
 'transport efficiency': 'Bruggeman' (possible: ['Bruggeman', 'ordered packing', 'hyperbola of revolution', 'overlapping spheres', 'tortuosity factor', 'random overlapping cylinders', 'heterogeneous catalyst', 'cation-exchange membrane'])
+'voltage as a state': 'false' (possible: ['false', 'true'])
 'working electrode': 'both' (possible: ['both', 'positive'])
 'x-average side reactions': 'false' (possible: ['false', 'true'])
 """
@@ -472,6 +473,17 @@ def test_save_load_model(self):
 
         os.remove("test_base_battery_model.json")
 
+    def test_voltage_as_state(self):
+        model = pybamm.lithium_ion.SPM({"voltage as a state": "true"})
+        assert model.options["voltage as a state"] == "true"
+        assert isinstance(model.variables["Voltage [V]"], pybamm.Variable)
+
+        model = pybamm.lithium_ion.SPM(
+            {"voltage as a state": "true", "operating mode": "voltage"}
+        )
+        assert model.options["voltage as a state"] == "true"
+        assert isinstance(model.variables["Voltage [V]"], pybamm.Variable)
+
 
 class TestOptions:
     def test_print_options(self):