From c76262747e7863f9d36f9cd377e4136bdbfd0414 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 22 Jul 2019 22:45:02 +0100 Subject: [PATCH] #546 some progress with Newman-Tiedemann --- pybamm/expression_tree/unary_operators.py | 2 +- pybamm/models/standard_variables.py | 15 ++++++++--- .../submodels/convection/base_convection.py | 12 ++++++--- .../electrode/ohm/surface_form_ohm.py | 13 ++++++--- .../full_stefan_maxwell_conductivity.py | 4 ++- ...urface_form_stefan_maxwell_conductivity.py | 27 +++++++++++++++---- .../full_stefan_maxwell_diffusion.py | 1 - .../submodels/interface/lead_acid_oxygen.py | 2 +- .../porosity/full_reaction_driven_porosity.py | 4 ++- 9 files changed, 60 insertions(+), 20 deletions(-) diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index c81f9878bd..b644f4d9b4 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -449,7 +449,7 @@ class IndefiniteIntegral(Integral): **Extends:** :class:`Integral` """ - def __init__(self, child, integration_variable): + def __init__(self, child, integration_variable, domain=None): if isinstance(integration_variable, list): if len(integration_variable) > 1: raise NotImplementedError( diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py index d42363f23e..cd9e161fa7 100644 --- a/pybamm/models/standard_variables.py +++ b/pybamm/models/standard_variables.py @@ -122,9 +122,18 @@ ) # Electrolyte pressure -pressure_n = pybamm.Variable("Negative electrolyte pressure", ["negative electrode"]) -pressure_s = pybamm.Variable("Separator electrolyte pressure", ["separator"]) -pressure_p = pybamm.Variable("Positive electrolyte pressure", ["positive electrode"]) +pressure_n = pybamm.SecondaryBroadcast( + pybamm.Variable("Negative electrolyte pressure", ["negative electrode"]), + "current collector", +) +pressure_s = pybamm.SecondaryBroadcast( + pybamm.Variable("Separator electrolyte pressure", ["separator"]), + "current collector", +) +pressure_p = pybamm.SecondaryBroadcast( + pybamm.Variable("Positive electrolyte pressure", ["positive electrode"]), + "current collector", +) pressure = pybamm.Concatenation(pressure_n, pressure_s, pressure_p) # Temperature diff --git a/pybamm/models/submodels/convection/base_convection.py b/pybamm/models/submodels/convection/base_convection.py index e80891fe1b..166afde8b6 100644 --- a/pybamm/models/submodels/convection/base_convection.py +++ b/pybamm/models/submodels/convection/base_convection.py @@ -123,9 +123,15 @@ def _separator_velocity(self, variables): # Simple formula for velocity in the separator dVbox_dz = pybamm.Concatenation( - pybamm.Broadcast(0, "negative electrode"), - pybamm.PrimaryBroadcast(-d_vbox_s__dx, "separator"), - pybamm.Broadcast(0, "positive electrode"), + pybamm.SecondaryBroadcast( + pybamm.Broadcast(0, "negative electrode"), "current collector" + ), + pybamm.SecondaryBroadcast( + pybamm.PrimaryBroadcast(-d_vbox_s__dx, "separator"), "current collector" + ), + pybamm.SecondaryBroadcast( + pybamm.Broadcast(0, "positive electrode"), "current collector" + ), ) v_box_s = pybamm.outer(d_vbox_s__dx, (x_s - l_n)) + pybamm.PrimaryBroadcast( v_box_n_right, "separator" diff --git a/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py b/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py index 5bcaf04af6..1433aca96e 100644 --- a/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py +++ b/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py @@ -33,6 +33,9 @@ def get_coupled_variables(self, variables): i_e = variables[self.domain + " electrolyte current density"] eps = variables[self.domain + " electrode porosity"] + if isinstance(i_boundary_cc, pybamm.Broadcast): + i_boundary_cc = i_boundary_cc.orphans[0] + i_s = i_boundary_cc - i_e if self.domain == "Negative": @@ -45,10 +48,12 @@ def get_coupled_variables(self, variables): delta_phi_p = variables["Positive electrode surface potential difference"] conductivity = param.sigma_p * (1 - eps) ** param.b - phi_s = ( - -pybamm.IndefiniteIntegral(i_s / conductivity, x_p) - + pybamm.boundary_value(phi_e_s, "right") - + pybamm.boundary_value(delta_phi_p, "left") + phi_s = -pybamm.IndefiniteIntegral( + i_s / conductivity, x_p + ) + pybamm.PrimaryBroadcast( + pybamm.boundary_value(phi_e_s, "right") + + pybamm.boundary_value(delta_phi_p, "left"), + "positive electrode", ) variables.update(self._get_standard_potential_variables(phi_s)) diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py index 06959cf524..ecebfd2198 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py @@ -51,7 +51,9 @@ def set_algebraic(self, variables): sum_j = sum( pybamm.Concatenation( variables[reaction["Negative"]["aj"]], - pybamm.Broadcast(0, "separator"), + pybamm.SecondaryBroadcast( + pybamm.Broadcast(0, "separator"), "current collector" + ), variables[reaction["Positive"]["aj"]], ) for reaction in self.reactions.values() diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py index 4ca2b3fe3b..d5b7d5ffdf 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py @@ -74,6 +74,9 @@ def set_boundary_conditions(self, variables): c_e = variables[self.domain + " electrolyte concentration"] delta_phi = variables[self.domain + " electrode surface potential difference"] + if isinstance(i_boundary_cc, pybamm.Broadcast): + i_boundary_cc = i_boundary_cc.orphans[0] + if self.domain == "Negative": c_e_flux = pybamm.BoundaryFlux(c_e, "right") flux_left = -pybamm.BoundaryValue(i_boundary_cc / sigma_eff, "left") @@ -151,6 +154,9 @@ def _get_neg_pos_coupled_variables(self, variables): c_e = variables[self.domain + " electrolyte concentration"] delta_phi = variables[self.domain + " electrode surface potential difference"] + if isinstance(i_boundary_cc, pybamm.Broadcast): + i_boundary_cc = i_boundary_cc.orphans[0] + i_e = conductivity * ( (param.chi(c_e) / c_e) * pybamm.grad(c_e) + pybamm.grad(delta_phi) @@ -180,10 +186,15 @@ def _get_sep_coupled_variables(self, variables): phi_e_n = variables["Negative electrolyte potential"] eps_s = variables["Separator porosity"] + if isinstance(i_boundary_cc, pybamm.Broadcast): + i_boundary_cc = i_boundary_cc.orphans[0] + chi_e_s = param.chi(c_e_s) kappa_s_eff = param.kappa_e(c_e_s) * (eps_s ** param.b) - phi_e_s = pybamm.boundary_value(phi_e_n, "right") + pybamm.IndefiniteIntegral( + phi_e_s = pybamm.PrimaryBroadcast( + pybamm.boundary_value(phi_e_n, "right"), "separator" + ) + pybamm.IndefiniteIntegral( chi_e_s / c_e_s * pybamm.grad(c_e_s) - param.C_e * i_boundary_cc / kappa_s_eff, x_s, @@ -213,6 +224,9 @@ def nasty_hack_to_get_phi_s(self, variables): i_boundary_cc = variables["Current collector current density"] i_e = variables[self.domain + " electrolyte current density"] + if isinstance(i_boundary_cc, pybamm.Broadcast): + i_boundary_cc = i_boundary_cc.orphans[0] + i_s = i_boundary_cc - i_e if self.domain == "Negative": @@ -225,10 +239,13 @@ def nasty_hack_to_get_phi_s(self, variables): delta_phi_p = variables["Positive electrode surface potential difference"] conductivity = param.sigma_p * (1 - eps) ** param.b - phi_s = ( - -pybamm.IndefiniteIntegral(i_s / conductivity, x_p) - + pybamm.boundary_value(phi_e_s, "right") - + pybamm.boundary_value(delta_phi_p, "left") + + phi_s = -pybamm.IndefiniteIntegral( + i_s / conductivity, x_p + ) + pybamm.PrimaryBroadcast( + pybamm.boundary_value(phi_e_s, "right") + + pybamm.boundary_value(delta_phi_p, "left"), + "positive electrode", ) return phi_s diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py index b45db477f2..4992ca9c89 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py @@ -60,7 +60,6 @@ def set_rhs(self, variables): N_e = variables["Electrolyte flux"] # i_e = variables["Electrolyte current density"] - # TODO: check lead acid version in new form # 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( diff --git a/pybamm/models/submodels/interface/lead_acid_oxygen.py b/pybamm/models/submodels/interface/lead_acid_oxygen.py index 5844471616..5b3429116e 100644 --- a/pybamm/models/submodels/interface/lead_acid_oxygen.py +++ b/pybamm/models/submodels/interface/lead_acid_oxygen.py @@ -42,7 +42,7 @@ def _get_exchange_current_density(self, variables): """ c_e = variables[self.domain + " electrolyte concentration"] # If c_e was broadcast, take only the orphan - if isinstance(c_e, pybamm.Broadcast): + if isinstance(c_e, pybamm.Broadcast) and c_e.broadcast_type != "secondary": c_e = c_e.orphans[0] if self.domain == "Negative": diff --git a/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py b/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py index 2ffd28e6c1..b10b8edb2f 100644 --- a/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py +++ b/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py @@ -32,7 +32,9 @@ def get_coupled_variables(self, variables): j_p = variables["Positive electrode interfacial current density"] deps_dt_n = -self.param.beta_surf_n * j_n - deps_dt_s = pybamm.Broadcast(0, ["separator"]) + deps_dt_s = pybamm.SecondaryBroadcast( + pybamm.Broadcast(0, ["separator"]), "current collector" + ) deps_dt_p = -self.param.beta_surf_p * j_p deps_dt = pybamm.Concatenation(deps_dt_n, deps_dt_s, deps_dt_p)