Merge pull request #921 from pybamm-team/issue-825-tplus

Issue 825 tplus

examples/scripts/compare_comsol/compare_comsol_DFN.py

@@ -78,18 +78,22 @@ def get_interp_fun(variable_name, domain):
     def myinterp(t):
         try:
             return interp.interp1d(
-                comsol_t, variable, fill_value="extrapolate", bounds_error=False
+                comsol_t, variable, fill_value="extrapolate", bounds_error=False,
             )(t)[:, np.newaxis]
         except ValueError as err:
             raise ValueError(
-                """Failed to interpolate '{}' with time range [{}, {}] at time {}.
-                Original error: {}""".format(
-                    variable_name, comsol_t[0], comsol_t[-1], t, err
-                )
+                (
+                    "Failed to interpolate '{}' with time range [{}, {}] at time {}."
+                    + "Original error: {}"
+                ).format(variable_name, comsol_t[0], comsol_t[-1], t, err)
             )
 
     # Make sure to use dimensional time
-    fun = pybamm.Function(myinterp, pybamm.t, name=variable_name + "_comsol")
+    fun = pybamm.Function(
+        myinterp,
+        pybamm.t * pybamm_model.timescale.evaluate(),
+        name=variable_name + "_comsol",
+    )
     fun.domain = domain
     fun.mesh = mesh.combine_submeshes(*domain)
     fun.secondary_mesh = None
@@ -109,7 +113,7 @@ def myinterp(t):
         fill_value="extrapolate",
         bounds_error=False,
     ),
-    pybamm.t,
+    pybamm.t * pybamm_model.timescale.evaluate(),
 )
 comsol_voltage.mesh = None
 comsol_voltage.secondary_mesh = None

pybamm/input/parameters/lead-acid/electrolytes/sulfuric_acid_Sulzer2019/parameters.csv

@@ -4,6 +4,7 @@ Name [units],Value,Reference,Notes
 # Electrolyte properties,,,
 Typical electrolyte concentration [mol.m-3],5650,,
 Cation transference number,0.7,,
+1 + dlnf/dlnc,1,,
 Partial molar volume of water [m3.mol-1],1.75E-05,,
 Partial molar volume of anions [m3.mol-1],3.15E-05,,
 Partial molar volume of cations [m3.mol-1],1.35E-05,,

pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Kim2011/parameters.csv

@@ -4,6 +4,7 @@ Name [units],Value,Reference,Notes
 # Electrolyte properties,,,
 Typical electrolyte concentration [mol.m-3],1200,,
 Cation transference number,0.4,Reported as a function in Kim2011 (Implement later),
+1 + dlnf/dlnc,1,,
 Electrolyte diffusivity [m2.s-1],[function]electrolyte_diffusivity_Kim2011,,
 Electrolyte conductivity [S.m-1],[function]electrolyte_conductivity_Kim2011,,
 ,,,

pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Marquis2019/parameters.csv

@@ -4,6 +4,7 @@ Name [units],Value,Reference,Notes
 # Electrolyte properties,,,
 Typical electrolyte concentration [mol.m-3],1000,Scott Moura FastDFN,
 Cation transference number,0.4,Scott Moura FastDFN,
+1 + dlnf/dlnc,1,,
 Electrolyte diffusivity [m2.s-1],[function]electrolyte_diffusivity_Capiglia1999,, 
 Electrolyte conductivity [S.m-1],[function]electrolyte_conductivity_Capiglia1999,,
 ,,,

pybamm/input/parameters/lithium-ion/electrolytes/lipf6_Nyman2008/parameters.csv

@@ -4,6 +4,7 @@ Name [units],Value,Reference,Notes
 # Electrolyte properties,,,
 Typical electrolyte concentration [mol.m-3],1000,Chen 2020,
 Cation transference number,0.2594,Chen 2020,
+1 + dlnf/dlnc,1,,
 Electrolyte diffusivity [m2.s-1],[function]electrolyte_diffusivity_Nyman2008,Nyman 2008," "
 Electrolyte conductivity [S.m-1],[function]electrolyte_conductivity_Nyman2008,Nyman 2008," "
 ,,,

pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py

@@ -63,19 +63,19 @@ def set_reactions(self):
         icd = " interfacial current density"
         self.reactions = {
             "main": {
-                "Negative": {"s": param.s_n, "aj": "Negative electrode" + icd},
-                "Positive": {"s": param.s_p, "aj": "Positive electrode" + icd},
+                "Negative": {"s": -param.s_plus_n_S, "aj": "Negative electrode" + icd},
+                "Positive": {"s": -param.s_plus_p_S, "aj": "Positive electrode" + icd},
             }
         }
         if "oxygen" in self.options["side reactions"]:
             self.reactions["oxygen"] = {
                 "Negative": {
-                    "s": -(param.s_plus_Ox + param.t_plus),
+                    "s": -param.s_plus_Ox,
                     "s_ox": -param.s_ox_Ox,
                     "aj": "Negative electrode oxygen" + icd,
                 },
                 "Positive": {
-                    "s": -(param.s_plus_Ox + param.t_plus),
+                    "s": -param.s_plus_Ox,
                     "s_ox": -param.s_ox_Ox,
                     "aj": "Positive electrode oxygen" + icd,
                 },

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -43,13 +43,7 @@ def set_reactions(self):
         icd = " interfacial current density"
         self.reactions = {
             "main": {
-                "Negative": {
-                    "s": 1 - self.param.t_plus,
-                    "aj": "Negative electrode" + icd,
-                },
-                "Positive": {
-                    "s": 1 - self.param.t_plus,
-                    "aj": "Positive electrode" + icd,
-                },
+                "Negative": {"s": 1, "aj": "Negative electrode" + icd},
+                "Positive": {"s": 1, "aj": "Positive electrode" + icd},
             }
         }

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -264,7 +264,7 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         ######################
         N_e = -tor * param.D_e(c_e, T) * pybamm.grad(c_e)
         self.rhs[c_e] = (1 / eps) * (
-            -pybamm.div(N_e) / param.C_e + (1 - param.t_plus) * j / param.gamma_e
+            -pybamm.div(N_e) / param.C_e + (1 - param.t_plus(c_e)) * j / param.gamma_e
         )
         self.boundary_conditions[c_e] = {
             "left": (pybamm.Scalar(0), "Neumann"),

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py

@@ -62,14 +62,9 @@ def get_coupled_variables(self, variables):
         kappa_p_av = param.kappa_e(c_e_av, T_av) * tor_p_av
 
         chi_av = param.chi(c_e_av)
-        if chi_av.domain == ["current collector"]:
-            chi_av_n = pybamm.PrimaryBroadcast(chi_av, "negative electrode")
-            chi_av_s = pybamm.PrimaryBroadcast(chi_av, "separator")
-            chi_av_p = pybamm.PrimaryBroadcast(chi_av, "positive electrode")
-        else:
-            chi_av_n = chi_av
-            chi_av_s = chi_av
-            chi_av_p = chi_av
+        chi_av_n = pybamm.PrimaryBroadcast(chi_av, "negative electrode")
+        chi_av_s = pybamm.PrimaryBroadcast(chi_av, "separator")
+        chi_av_p = pybamm.PrimaryBroadcast(chi_av, "positive electrode")
 
         # electrolyte current
         i_e_n = i_boundary_cc_0 * x_n / l_n

pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/composite_stefan_maxwell_diffusion.py

@@ -73,12 +73,16 @@ def set_rhs(self, variables):
         }
 
     def _get_source_terms_leading_order(self, variables):
+        param = self.param
+        c_e_n = variables["Negative electrolyte concentration"]
+        c_e_p = variables["Positive electrolyte concentration"]
+
         return sum(
             pybamm.Concatenation(
-                reaction["Negative"]["s"]
+                (reaction["Negative"]["s"] - param.t_plus(c_e_n))
                 * variables["Leading-order " + reaction["Negative"]["aj"].lower()],
                 pybamm.FullBroadcast(0, "separator", "current collector"),
-                reaction["Positive"]["s"]
+                (reaction["Positive"]["s"] - param.t_plus(c_e_p))
                 * variables["Leading-order " + reaction["Positive"]["aj"].lower()],
             )
             / self.param.gamma_e

pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/first_order_stefan_maxwell_diffusion.py

@@ -57,15 +57,15 @@ def get_coupled_variables(self, variables):
 
         # Right-hand sides
         rhs_n = d_epsc_n_0_dt - sum(
-            reaction["Negative"]["s"]
+            (reaction["Negative"]["s"] - param.t_plus(c_e_0))
             * variables[
                 "Leading-order x-averaged " + reaction["Negative"]["aj"].lower()
             ]
             for reaction in self.reactions.values()
         )
         rhs_s = d_epsc_s_0_dt
         rhs_p = d_epsc_p_0_dt - sum(
-            reaction["Positive"]["s"]
+            (reaction["Positive"]["s"] - param.t_plus(c_e_0))
             * variables[
                 "Leading-order x-averaged " + reaction["Positive"]["aj"].lower()
             ]

pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py

@@ -59,15 +59,16 @@ def set_rhs(self, variables):
         deps_dt = variables["Porosity change"]
         c_e = variables["Electrolyte concentration"]
         N_e = variables["Electrolyte flux"]
-        # i_e = variables["Electrolyte current density"]
+        c_e_n = variables["Negative electrolyte concentration"]
+        c_e_p = variables["Positive electrolyte concentration"]
 
-        # source_term = ((param.s - param.t_plus) / param.gamma_e) * pybamm.div(i_e)
-        # source_term = pybamm.div(i_e) / param.gamma_e  # lithium-ion
         source_terms = sum(
             pybamm.Concatenation(
-                reaction["Negative"]["s"] * variables[reaction["Negative"]["aj"]],
+                (reaction["Negative"]["s"] - param.t_plus(c_e_n))
+                * variables[reaction["Negative"]["aj"]],
                 pybamm.FullBroadcast(0, "separator", "current collector"),
-                reaction["Positive"]["s"] * variables[reaction["Positive"]["aj"]],
+                (reaction["Positive"]["s"] - param.t_plus(c_e_p))
+                * variables[reaction["Positive"]["aj"]],
             )
             / param.gamma_e
             for reaction in self.reactions.values()

pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py

@@ -60,10 +60,10 @@ def set_rhs(self, variables):
 
         source_terms = sum(
             param.l_n
-            * rxn["Negative"]["s"]
+            * (rxn["Negative"]["s"] - param.t_plus(c_e_av))
             * variables["X-averaged " + rxn["Negative"]["aj"].lower()]
             + param.l_p
-            * rxn["Positive"]["s"]
+            * (rxn["Positive"]["s"] - param.t_plus(c_e_av))
             * variables["X-averaged " + rxn["Positive"]["aj"].lower()]
             for rxn in self.reactions.values()
         )

pybamm/parameters/standard_parameters_lead_acid.py

@@ -49,7 +49,6 @@
 
 # Electrolyte properties
 c_e_typ = pybamm.Parameter("Typical electrolyte concentration [mol.m-3]")
-t_plus = pybamm.Parameter("Cation transference number")
 V_w = pybamm.Parameter("Partial molar volume of water [m3.mol-1]")
 V_plus = pybamm.Parameter("Partial molar volume of cations [m3.mol-1]")
 V_minus = pybamm.Parameter("Partial molar volume of anions [m3.mol-1]")
@@ -171,6 +170,11 @@
 "2. Dimensional Functions"
 
 
+def t_plus(c_e):
+    "Dimensionless transference number (i.e. c_e is dimensional)"
+    return pybamm.FunctionParameter("Cation transference number", c_e * c_e_typ)
+
+
 def D_e_dimensional(c_e, T):
     "Dimensional diffusivity in electrolyte"
     return pybamm.FunctionParameter("Electrolyte diffusivity [m2.s-1]", c_e)
@@ -310,7 +314,7 @@ def U_p_dimensional(c_e, T):
 centre_z_tab_p = pybamm.geometric_parameters.centre_z_tab_p
 
 # Diffusive kinematic relationship coefficient
-omega_i = c_e_typ * M_e / rho_typ * (t_plus + M_minus / M_e)
+omega_i = c_e_typ * M_e / rho_typ * (t_plus(1) + M_minus / M_e)
 # Migrative kinematic relationship coefficient (electrolyte)
 omega_c_e = c_e_typ * M_e / rho_typ * (1 - M_w * V_e / V_w * M_e)
 C_e = tau_diffusion_e / tau_discharge
@@ -347,12 +351,10 @@ def U_p_dimensional(c_e, T):
 # Main
 s_plus_n_S = s_plus_n_S_dim / ne_n_S
 s_plus_p_S = s_plus_p_S_dim / ne_p_S
-s_n = -(s_plus_n_S + t_plus)  # Dimensionless rection rate (neg)
-s_p = -(s_plus_p_S + t_plus)  # Dimensionless rection rate (pos)
-s = pybamm.Concatenation(
-    pybamm.FullBroadcast(s_n, ["negative electrode"], "current collector"),
+s_plus_S = pybamm.Concatenation(
+    pybamm.FullBroadcast(s_plus_n_S, ["negative electrode"], "current collector"),
     pybamm.FullBroadcast(0, ["separator"], "current collector"),
-    pybamm.FullBroadcast(s_p, ["positive electrode"], "current collector"),
+    pybamm.FullBroadcast(s_plus_p_S, ["positive electrode"], "current collector"),
 )
 j0_n_S_ref = j0_n_S_ref_dimensional / interfacial_current_scale_n
 j0_p_S_ref = j0_p_S_ref_dimensional / interfacial_current_scale_p
@@ -407,7 +409,9 @@ def U_p_dimensional(c_e, T):
 ) / potential_scale
 
 # Electrolyte volumetric capacity
-Q_e_max = (l_n * eps_n_max + l_s * eps_s_max + l_p * eps_p_max) / (s_p - s_n)
+Q_e_max = (l_n * eps_n_max + l_s * eps_s_max + l_p * eps_p_max) / (
+    s_plus_n_S - s_plus_p_S
+)
 Q_e_max_dimensional = Q_e_max * c_e_typ * F
 capacity = Q_e_max_dimensional * n_electrodes_parallel * A_cs * L_x
 
@@ -490,7 +494,7 @@ def kappa_e(c_e, T):
 def chi(c_e, c_ox=0, c_hy=0):
     return (
         chi_dimensional(c_e_typ * c_e)
-        * (2 * (1 - t_plus))
+        * (2 * (1 - t_plus(c_e)))
         / (V_w * c_T(c_e_typ * c_e, c_e_typ * c_ox, c_e_typ * c_hy))
     )
 

pybamm/parameters/standard_parameters_lithium_ion.py

@@ -326,17 +326,24 @@ def U_p_dimensional(sto, T):
 alpha_prime = alpha / delta
 
 # Electrolyte Properties
-t_plus = pybamm.Parameter("Cation transference number")
+
+
+def t_plus(c_e):
+    return pybamm.FunctionParameter("Cation transference number", c_e)
+
+
+def one_plus_dlnf_dlnc(c_e):
+    return pybamm.FunctionParameter("1 + dlnf/dlnc", c_e)
+
+
 beta_surf = pybamm.Scalar(0)
-s = 1 - t_plus
 
 
 # (1-2*t_plus) is for Nernst-Planck
 # 2*(1-t_plus) for Stefan-Maxwell
 # Bizeray et al (2016) "Resolving a discrepancy ..."
-# note: this is a function for consistancy with lead-acid
 def chi(c_e):
-    return 2 * (1 - t_plus)
+    return (2 * (1 - t_plus(c_e))) * (one_plus_dlnf_dlnc(c_e))
 
 
 # Electrochemical Reactions

tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_leading_stefan_maxwell_diffusion.py

@@ -13,11 +13,11 @@ def test_public_functions(self):
         reactions = {
             "main": {
                 "Negative": {
-                    "s": param.s_n,
+                    "s": -param.s_plus_n_S,
                     "aj": "Negative electrode interfacial current density",
                 },
                 "Positive": {
-                    "s": param.s_p,
+                    "s": -param.s_plus_p_S,
                     "aj": "Positive electrode interfacial current density",
                 },
             }

tests/unit/test_parameters/test_standard_parameters_lead_acid.py

@@ -50,7 +50,7 @@ def test_parameters_defaults_lead_acid(self):
 
     def test_concatenated_parameters(self):
         # create
-        s_param = pybamm.standard_parameters_lead_acid.s
+        s_param = pybamm.standard_parameters_lead_acid.s_plus_S
         self.assertIsInstance(s_param, pybamm.Concatenation)
         self.assertEqual(
             s_param.domain, ["negative electrode", "separator", "positive electrode"]