diff --git a/CHANGELOG.md b/CHANGELOG.md index c37300d743..842e200466 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Added three-electrode EIS support to `pybamm.EISSimulation`, including automatic default reference-electrode insertion, named positive/negative electrode impedance outputs, and component-aware Nyquist plotting. ([#5648](https://github.com/pybamm-team/PyBaMM/pull/5648)) - Added unstructured mesh support (`UnstructuredSubMesh`, generators, and interface coupling) for arbitrary 2D/3D domains. Hexahedra must have planar faces (warped hexes raise a `GeometryError`), and `UserSuppliedUnstructuredMesh` accepts tetrahedral, triangular, and quadrilateral cells only. ([#5687](https://github.com/pybamm-team/PyBaMM/pull/5687)) - Generalised `VectorField` to N components and added `Component`/`Norm` operators for multi-dimensional vector fields. ([#5686](https://github.com/pybamm-team/PyBaMM/pull/5686)) - Removed the left sidebar from the documentation home page for a cleaner landing experience. ([#5699](https://github.com/pybamm-team/PyBaMM/pull/5699)) diff --git a/examples/scripts/run_eis_three_electrode_simulation.py b/examples/scripts/run_eis_three_electrode_simulation.py new file mode 100644 index 0000000000..9efd7ac4b7 --- /dev/null +++ b/examples/scripts/run_eis_three_electrode_simulation.py @@ -0,0 +1,14 @@ +import numpy as np + +import pybamm + +model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}) + +frequencies = np.logspace(-4, 4, 50) + +eis_sim = pybamm.EISSimulation( + model, + three_electrodes=True, +) +eis_sim.solve(frequencies) +eis_sim.nyquist_plot() diff --git a/packages/pybamm/src/pybamm/simulation/eis_simulation.py b/packages/pybamm/src/pybamm/simulation/eis_simulation.py index 0c7ebdfdef..4c8f676a4f 100644 --- a/packages/pybamm/src/pybamm/simulation/eis_simulation.py +++ b/packages/pybamm/src/pybamm/simulation/eis_simulation.py @@ -29,6 +29,10 @@ class EISSimulation(BaseSimulation): A dictionary of the types of spatial method to use on each domain. skip_surface_form_check : bool, optional If True, skip the 'surface form' model option validation. Defaults to False. + three_electrodes: bool, optional + Whether to use a three-electrode model. If True, the model must have a + reference electrode inserted. Default is False. If True, returns the impedance + of both the working and counter electrodes, as well as the cell impedance. """ def __init__( @@ -40,17 +44,28 @@ def __init__( var_pts=None, spatial_methods=None, skip_surface_form_check=False, + three_electrodes=False, ): timer = pybamm.Timer() model_name = model.name parameter_values = parameter_values or model.default_parameter_values + self.three_electrodes = three_electrodes + if three_electrodes and ( + "Positive electrode 3E potential [V]" not in model.variables + or "Negative electrode 3E potential [V]" not in model.variables + ): + reference_position = ( + parameter_values["Negative electrode thickness [m]"] + + parameter_values["Separator thickness [m]"] / 2 + ) + model.insert_reference_electrode(reference_position) # Validate required variables and surface form before any processing - self._validate_model_for_eis(model, skip_surface_form_check) + self._validate_model_for_eis(model, skip_surface_form_check, three_electrodes) pybamm.logger.info(f"Setting up {model_name} for EIS") - eis_model = self._set_up_model_for_eis(model) + eis_model = self._set_up_model_for_eis(model, three_electrodes) # Compute impedance scale factor after model transformation: the # external-circuit FunctionControl swap rescales "Current [A]" by @@ -78,7 +93,9 @@ def __init__( pybamm.citations.register("Hallemans2025") @staticmethod - def _validate_model_for_eis(model, skip_surface_form_check=False): + def _validate_model_for_eis( + model, skip_surface_form_check=False, three_electrodes=False + ): """Validate that a model is suitable for frequency-domain EIS. Raises @@ -87,6 +104,14 @@ def _validate_model_for_eis(model, skip_surface_form_check=False): If the model is missing required variables or options. """ required_vars = ["Voltage [V]", "Current [A]"] + if three_electrodes: + required_vars.extend( + [ + "Positive electrode 3E potential [V]", + "Negative electrode 3E potential [V]", + ] + ) + for var in required_vars: if var not in model.variables: raise ValueError( @@ -105,7 +130,7 @@ def _validate_model_for_eis(model, skip_surface_form_check=False): ) @staticmethod - def _set_up_model_for_eis(model): + def _set_up_model_for_eis(model, three_electrodes=False): """Prepare a model for frequency-domain EIS. Creates a copy of the model with voltage and current as algebraic @@ -132,6 +157,27 @@ def _set_up_model_for_eis(model): new_model.algebraic[V_cell] = V_cell - V new_model.initial_conditions[V_cell] = new_model.param.ocv_init + if three_electrodes: + V_pos_3e = new_model.variables["Positive electrode 3E potential [V]"] + V_pos_3e_var = pybamm.Variable( + "Positive electrode 3E potential variable [V]" + ) + new_model.variables["Positive electrode 3E potential variable [V]"] = ( + V_pos_3e_var + ) + new_model.algebraic[V_pos_3e_var] = V_pos_3e_var - V_pos_3e + new_model.initial_conditions[V_pos_3e_var] = new_model.param.p.prim.U_init + + V_neg_3e = new_model.variables["Negative electrode 3E potential [V]"] + V_neg_3e_var = pybamm.Variable( + "Negative electrode 3E potential variable [V]" + ) + new_model.variables["Negative electrode 3E potential variable [V]"] = ( + V_neg_3e_var + ) + new_model.algebraic[V_neg_3e_var] = V_neg_3e_var - V_neg_3e + new_model.initial_conditions[V_neg_3e_var] = new_model.param.n.prim.U_init + # Replace current with a FunctionControl variable external_circuit_variables = pybamm.external_circuit.FunctionControl( model.param, None, model.options, control="algebraic" @@ -212,7 +258,7 @@ def _build_matrix_problem(self, inputs_dict=None): return self._cached_M, neg_J, self._cached_b @staticmethod - def _calculate_impedance(frequency, M, neg_J, b): + def _calculate_impedance(frequency, M, neg_J, b, three_electrodes=False): """Calculate impedance at a single frequency. Parameters @@ -233,9 +279,15 @@ def _calculate_impedance(frequency, M, neg_J, b): """ A = 1.0j * 2 * np.pi * frequency * M + neg_J x = spsolve(A, b) - # Voltage is penultimate, current is last (by construction in - # _set_up_model_for_eis) - return -x[-2] / x[-1] + if three_electrodes: + z_cell = -x[-4] / x[-1] + z_pos = -x[-3] / x[-1] + z_neg = -x[-2] / x[-1] + return z_cell, z_pos, z_neg + else: + # Voltage is penultimate, current is last (by construction in + # _set_up_model_for_eis) + return -x[-2] / x[-1] def solve(self, frequencies, inputs=None, initial_soc=None): """Compute impedance at the given frequencies. @@ -269,9 +321,20 @@ def solve(self, frequencies, inputs=None, initial_soc=None): M, neg_J, b = self._build_matrix_problem(inputs_dict=inputs) - zs = [self._calculate_impedance(f, M, neg_J, b) for f in frequencies] + zs = [ + self._calculate_impedance(f, M, neg_J, b, self.three_electrodes) + for f in frequencies + ] impedance = np.array(zs) * self._z_scale - self._solution = pybamm.EISSolution(frequencies, impedance) + if self.three_electrodes: + self._solution = pybamm.EISSolution(frequencies, impedance[:, 0]) + self._solution._data["Cell impedance [Ohm]"] = impedance[:, 0] + self._solution._data["Positive electrode impedance [Ohm]"] = impedance[:, 1] + self._solution._data["Negative electrode impedance [Ohm]"] = -impedance[ + :, 2 + ] + else: + self._solution = pybamm.EISSolution(frequencies, impedance) self._solution.set_up_time = self.set_up_time self.solve_time = timer.time() diff --git a/packages/pybamm/src/pybamm/solvers/solution.py b/packages/pybamm/src/pybamm/solvers/solution.py index b67ce33c7e..d78dbe3b48 100644 --- a/packages/pybamm/src/pybamm/solvers/solution.py +++ b/packages/pybamm/src/pybamm/solvers/solution.py @@ -173,7 +173,35 @@ def nyquist_plot(self, **kwargs): """ from pybamm.plotting.nyquist_plot import nyquist_plot - return nyquist_plot(self.impedance, **kwargs) + component_keys = { + "Cell": "Cell impedance [Ohm]", + "Positive electrode": "Positive electrode impedance [Ohm]", + "Negative electrode": "Negative electrode impedance [Ohm]", + } + if not all(key in self._data for key in component_keys.values()): + return nyquist_plot(self.impedance, **kwargs) + + plot_kwargs = dict(kwargs) + ax = plot_kwargs.pop("ax", None) + show_plot = plot_kwargs.pop("show_plot", True) + plot_kwargs.pop("label", None) + fig = None + for label, key in component_keys.items(): + fig_i, ax = nyquist_plot( + self._data[key], + ax=ax, + show_plot=False, + label=label, + **plot_kwargs, + ) + fig = fig or fig_i + ax.legend() + + if show_plot: # pragma: no cover + plt = pybamm.import_optional_dependency("matplotlib.pyplot") + plt.show() + + return fig, ax _DEFAULT_SOLUTION_OPTIONS = { diff --git a/packages/pybamm/tests/unit/test_eis_simulation.py b/packages/pybamm/tests/unit/test_eis_simulation.py index 9d30746f94..c6e4d1484f 100644 --- a/packages/pybamm/tests/unit/test_eis_simulation.py +++ b/packages/pybamm/tests/unit/test_eis_simulation.py @@ -29,6 +29,12 @@ def test_skip_surface_form_check(self): model, skip_surface_form_check=True ) + def test_three_electrode_inserts_reference_electrode(self): + model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}) + pybamm.EISSimulation(model, three_electrodes=True) + assert "Positive electrode 3E potential [V]" in model.variables + assert "Negative electrode 3E potential [V]" in model.variables + class TestEISSolution: """Tests for the EISSolution class.""" @@ -251,6 +257,23 @@ def test_high_freq_intercept_matches_contact_resistance(self): assert high_freq_z.real == pytest.approx(1.0, abs=1e-3) assert abs(high_freq_z.imag) < 1e-3 + @pytest.mark.parametrize( + "model_class", + [pybamm.lithium_ion.SPM, pybamm.lithium_ion.SPMe, pybamm.lithium_ion.DFN], + ) + def test_three_electrode_impedances_sum_to_cell_impedance(self, model_class): + model = model_class(options={"surface form": "differential"}) + eis_sim = pybamm.EISSimulation(model, three_electrodes=True) + frequencies = np.logspace(-2, 2, 5) + + result = eis_sim.solve(frequencies) + z_cell = result["Cell impedance [Ohm]"] + z_pos = result["Positive electrode impedance [Ohm]"] + z_neg = result["Negative electrode impedance [Ohm]"] + + np.testing.assert_allclose(z_pos + z_neg, z_cell, rtol=1e-5, atol=1e-5) + np.testing.assert_allclose(result.impedance, z_cell) + class TestNyquistPlot: """Tests for Nyquist plotting.""" @@ -290,6 +313,26 @@ def test_eis_nyquist_plot(self): assert fig is not None assert ax is not None + def test_eis_nyquist_plot_components(self): + import matplotlib + + matplotlib.use("Agg") + + impedance = np.array([1 + 0.5j, 2 + 1j, 3 + 1.5j]) + solution = pybamm.EISSolution(np.array([1.0, 10.0, 100.0]), impedance) + solution._data["Cell impedance [Ohm]"] = impedance + solution._data["Positive electrode impedance [Ohm]"] = 0.4 * impedance + solution._data["Negative electrode impedance [Ohm]"] = 0.6 * impedance + + fig, ax = solution.nyquist_plot(show_plot=False) + + assert fig is not None + assert [line.get_label() for line in ax.get_lines()] == [ + "Cell", + "Positive electrode", + "Negative electrode", + ] + def test_nyquist_plot_before_solve_raises(self): model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}) eis_sim = pybamm.EISSimulation(model)