fix calculation of bulk ocp in hysteresis models (#5169)

* fix calculation of bulk ocp

* remove print statements

Author: rtimms

CHANGELOG.md

@@ -9,6 +9,7 @@
 
 ## Bug fixes
 
+- Fix a bug in the calculation of "Bulk" OCP terms in hysteresis models ([#5169](https://github.com/pybamm-team/PyBaMM/pull/5169))
 - Fixed a bug where the final duration of a drive cycle would not be inferred correctly. ([#5153](https://github.com/pybamm-team/PyBaMM/pull/5153))
 
 # [v25.8.0](https://github.com/pybamm-team/PyBaMM/tree/v25.8.0) - 2025-08-04

src/pybamm/models/submodels/interface/open_circuit_potential/base_hysteresis_ocp.py

@@ -78,13 +78,15 @@ def _get_coupled_variables(self, variables):
 
         if self.reaction == "lithium-ion main":
             # Get equilibrium, delithiation and lithiation OCPs
-            sto_surf, T = self._get_stoichiometry_and_temperature(variables)
+            sto_surf, sto_bulk, T, T_bulk = self._get_stoichiometry_and_temperature(
+                variables
+            )
             U_eq = self.phase_param.U(sto_surf, T)
             U_eq_x_av = self.phase_param.U(sto_surf, T)
             U_lith = self.phase_param.U(sto_surf, T, "lithiation")
-            U_lith_x_av = pybamm.x_average(U_lith)
+            U_lith_bulk = self.phase_param.U(sto_bulk, T_bulk)
             U_delith = self.phase_param.U(sto_surf, T, "delithiation")
-            U_delith_x_av = pybamm.x_average(U_delith)
+            U_delith_bulk = self.phase_param.U(sto_bulk, T_bulk, "delithiation")
 
             H = U_lith - U_delith
             H_x_av = pybamm.x_average(H)
@@ -122,11 +124,8 @@ def _get_coupled_variables(self, variables):
                 ocp_surf = (1 + h) / 2 * U_delith + (1 - h) / 2 * U_lith
 
             # Size average
-            U_delith_s_av = pybamm.size_average(U_delith_x_av)
-            U_lith_s_av = pybamm.size_average(U_lith_x_av)
             h_s_av = pybamm.size_average(h_x_av)
-
-            ocp_bulk = (1 + h_s_av) / 2 * U_delith_s_av + (1 - h_s_av) / 2 * U_lith_s_av
+            ocp_bulk = (1 + h_s_av) / 2 * U_delith_bulk + (1 - h_s_av) / 2 * U_lith_bulk
             dUdT = self.phase_param.dUdT(sto_surf)
 
         variables.update(self._get_standard_ocp_variables(ocp_surf, ocp_bulk, dUdT))

src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py

@@ -132,7 +132,9 @@ def _get_stoichiometry_and_temperature(self, variables):
         domain_options = getattr(self.options, domain)
         phase_name = self.phase_name
 
+        sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"]
         T = variables[f"{Domain} electrode temperature [K]"]
+        T_bulk = pybamm.xyzs_average(T)
 
         # For "particle-size distribution" models, take distribution version
         # of sto_surf that depends on particle size.
@@ -156,4 +158,4 @@ def _get_stoichiometry_and_temperature(self, variables):
                 sto_surf = sto_surf.orphans[0]
                 T = T.orphans[0]
 
-        return sto_surf, T
+        return sto_surf, sto_bulk, T, T_bulk

src/pybamm/models/submodels/interface/open_circuit_potential/one_state_differential_capacity_hysteresis_ocp.py

@@ -47,9 +47,8 @@ def get_coupled_variables(self, variables):
 
         if self.reaction == "lithium-ion main":
             # determine dQ/dU
-            sto_surf, T = self._get_stoichiometry_and_temperature(variables)
+            sto_surf, _, T, T_bulk = self._get_stoichiometry_and_temperature(variables)
             sto_surf = pybamm.size_average(sto_surf)
-            T_bulk = pybamm.xyz_average(pybamm.size_average(T))
             if phase_name == "":
                 Q_mag = variables[f"{Domain} electrode capacity [A.h]"]
             else:

src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py

@@ -12,15 +12,15 @@ def get_coupled_variables(self, variables):
         phase_name = self.phase_name
 
         if self.reaction == "lithium-ion main":
-            sto_surf, T = self._get_stoichiometry_and_temperature(variables)
+            sto_surf, sto_bulk, T, T_bulk = self._get_stoichiometry_and_temperature(
+                variables
+            )
 
             ocp_surf = self.phase_param.U(sto_surf, T)
-            dUdT = self.phase_param.dUdT(sto_surf)
-
-            sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"]
-            T_bulk = pybamm.xyzs_average(T)
             ocp_bulk = self.phase_param.U(sto_bulk, T_bulk)
 
+            dUdT = self.phase_param.dUdT(sto_surf)
+
             variables.update(
                 {
                     f"{Domain} electrode {phase_name}equilibrium open-circuit potential [V]": ocp_surf,