Merge pull request #3330 from pybamm-team/general-thermal-bcs

General thermal bcs

CHANGELOG.md

@@ -8,6 +8,8 @@
 - Implement the MSMR model ([#3116](https://github.com/pybamm-team/PyBaMM/pull/3116))
 
 ## Bug fixes
+
+- Fixed a bug where there was a missing thermal conductivity in the thermal pouch cell models ([#3330](https://github.com/pybamm-team/PyBaMM/pull/3330))
 - Fixed a bug that caused incorrect results of “{Domain} electrode thickness change [m]” due to the absence of dimension for the variable `electrode_thickness_change`([#3329](https://github.com/pybamm-team/PyBaMM/pull/3329)).
 - Fixed a bug that occured in `check_ys_are_not_too_large` when trying to reference `y-slice` where the referenced variable was not a `pybamm.StateVector` ([#3313](https://github.com/pybamm-team/PyBaMM/pull/3313)
 - Fixed a bug with `_Heaviside._evaluate_for_shape` which meant some expressions involving heaviside function and subtractions did not work ([#3306](https://github.com/pybamm-team/PyBaMM/pull/3306))

pybamm/models/submodels/current_collector/potential_pair.py

@@ -80,12 +80,13 @@ def set_boundary_conditions(self, variables):
 
         param = self.param
         applied_current_density = variables["Total current density [A.m-2]"]
-        cc_area = self._get_effective_current_collector_area()
+        total_current = applied_current_density * param.A_cc
 
-        # cc_area appears here due to choice of non-dimensionalisation
-        pos_tab_bc = (
-            -applied_current_density * cc_area / (param.p.sigma_cc * param.p.L_cc)
-        )
+        # In the 1+1D model, the behaviour is averaged over the y-direction, so the
+        # effective tab area is the cell width multiplied by the current collector
+        # thickness
+        positive_tab_area = param.L_y * param.p.L_cc
+        pos_tab_bc = -total_current / (param.p.sigma_cc * positive_tab_area)
 
         # Boundary condition needs to be on the variables that go into the Laplacian,
         # even though phi_s_cp isn't a pybamm.Variable object
@@ -100,10 +101,6 @@ def set_boundary_conditions(self, variables):
             },
         }
 
-    def _get_effective_current_collector_area(self):
-        """In the 1+1D models the current collector effectively has surface area l_z"""
-        return self.param.L_z
-
 
 class PotentialPair2plus1D(BasePotentialPair):
     """Base class for a 2+1D potential pair model"""
@@ -117,21 +114,18 @@ def set_boundary_conditions(self, variables):
 
         param = self.param
         applied_current_density = variables["Total current density [A.m-2]"]
-        cc_area = self._get_effective_current_collector_area()
+        total_current = applied_current_density * param.A_cc
 
         # Note: we divide by the *numerical* tab area so that the correct total
         # current is applied. That is, numerically integrating the current density
         # around the boundary gives the applied current exactly.
-
         positive_tab_area = pybamm.BoundaryIntegral(
             pybamm.PrimaryBroadcast(param.p.L_cc, "current collector"),
             region="positive tab",
         )
 
         # cc_area appears here due to choice of non-dimensionalisation
-        pos_tab_bc = (
-            -applied_current_density * cc_area / (param.p.sigma_cc * positive_tab_area)
-        )
+        pos_tab_bc = -total_current / (param.p.sigma_cc * positive_tab_area)
 
         # Boundary condition needs to be on the variables that go into the Laplacian,
         # even though phi_s_cp isn't a pybamm.Variable object
@@ -160,7 +154,3 @@ def set_boundary_conditions(self, variables):
                 "positive tab": (pos_tab_bc, "Neumann"),
             },
         }
-
-    def _get_effective_current_collector_area(self):
-        """Return the area of the current collector."""
-        return self.param.L_y * self.param.L_z

pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py

@@ -82,7 +82,7 @@ def set_rhs(self, variables):
 
         self.rhs = {
             T_av: (
-                pybamm.laplacian(T_av)
+                pybamm.div(self.param.lambda_eff(T_av) * pybamm.grad(T_av))
                 + Q_av
                 + total_cooling_coefficient * (T_av - T_amb)
             )
@@ -131,25 +131,28 @@ def set_boundary_conditions(self, variables):
         bottom_cooling_coefficient = (
             param.n.h_tab * neg_tab_area * neg_tab_bottom_bool
             + param.p.h_tab * pos_tab_area * pos_tab_bottom_bool
-            + 10 * non_tab_bottom_area
+            + param.h_edge(L_y / 2, 0) * non_tab_bottom_area
         ) / total_area
 
         # just use left and right for clarity
         # left = bottom of cell (z=0)
         # right = top of cell (z=L_z)
+        lambda_eff = param.lambda_eff(T_av)
         self.boundary_conditions = {
             T_av: {
                 "left": (
-                    bottom_cooling_coefficient
-                    * pybamm.boundary_value(
-                        T_av - T_amb,
+                    pybamm.boundary_value(
+                        bottom_cooling_coefficient * (T_av - T_amb),
                         "left",
-                    ),
+                    )
+                    / pybamm.boundary_value(lambda_eff, "left"),
                     "Neumann",
                 ),
                 "right": (
-                    -top_cooling_coefficient
-                    * pybamm.boundary_value(T_av - T_amb, "right"),
+                    pybamm.boundary_value(
+                        -top_cooling_coefficient * (T_av - T_amb), "right"
+                    )
+                    / pybamm.boundary_value(lambda_eff, "right"),
                     "Neumann",
                 ),
             }

pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py

@@ -82,9 +82,11 @@ def set_rhs(self, variables):
         # correct mass matrix when discretised. The first argument is the source term
         # and the second argument is the variable governed by the equation that the
         # source term appears in.
+        # Note: not correct if lambda_eff is a function of T_av - need to implement div
+        # in 2D rather than doing laplacian directly
         self.rhs = {
             T_av: (
-                pybamm.laplacian(T_av)
+                self.param.lambda_eff(T_av) * pybamm.laplacian(T_av)
                 + pybamm.source(Q_av, T_av)
                 + pybamm.source(yz_surface_cooling_coefficient * (T_av - T_amb), T_av)
                 + pybamm.source(
@@ -112,10 +114,12 @@ def set_boundary_conditions(self, variables):
         )
 
         negative_tab_bc = pybamm.boundary_value(
-            -h_tab_n_corrected * (T_av - T_amb), "negative tab"
+            -h_tab_n_corrected * (T_av - T_amb) / self.param.n.lambda_cc(T_av),
+            "negative tab",
         )
         positive_tab_bc = pybamm.boundary_value(
-            -h_tab_p_corrected * (T_av - T_amb), "positive tab"
+            -h_tab_p_corrected * (T_av - T_amb) / self.param.p.lambda_cc(T_av),
+            "positive tab",
         )
 
         self.boundary_conditions = {

pybamm/parameters/lithium_ion_parameters.py

@@ -53,6 +53,7 @@ def _set_parameters(self):
         self.h_edge = self.therm.h_edge
         self.h_total = self.therm.h_total
         self.rho_c_p_eff = self.therm.rho_c_p_eff
+        self.lambda_eff = self.therm.lambda_eff
 
         # Macroscale geometry
         self.L_x = self.geo.L_x

pybamm/parameters/thermal_parameters.py

@@ -72,6 +72,16 @@ def rho_c_p_eff(self, T):
             + self.p.rho_c_p_cc(T) * self.geo.p.L_cc
         ) / self.geo.L
 
+    def lambda_eff(self, T):
+        """Effective thermal conductivity [W.m-1.K-1]"""
+        return (
+            self.n.lambda_cc(T) * self.geo.n.L_cc
+            + self.n.lambda_(T) * self.geo.n.L
+            + self.s.lambda_(T) * self.geo.s.L
+            + self.p.lambda_(T) * self.geo.p.L
+            + self.p.lambda_cc(T) * self.geo.p.L_cc
+        ) / self.geo.L
+
 
 class DomainThermalParameters(BaseParameters):
     def __init__(self, domain, main_param):