Problem` Description
I’m trying to implement a thermal model within the DFN base model by extending the standard DFN model to include temperature dynamics instead of using a constant temperature. The model implementation includes heat generation terms (irreversible, ohmic, ionic, reversible entropic, and contact resistance heat) and a thermal diffusion equation.
Code Implementation
I’ve extended the BaseModel class to create a ThermalDFN class that includes temperature dynamics.
Error
(well i had some problems in importing the model but tried to solve them and don’t have any error in the thermalmodel) but When trying to solve the model ( with experiment … sim.solve… ) I get " ValueError: dimension mismatch"
I tried then to check the temp variation dTdt without +Q1 with putting instead of T in the boundary : pybamm.boundary_value(T, “left”) and it worked!
but if i add the Q1 with same boundaries i have this error:
pybamm.expression_tree.exceptions.ShapeError: Cannot find shape (original error: operands could not be broadcast together with shapes (21,1) (19,1) )
Questions
- Is there an issue with how I’m defining the domains for the temperature variables?
- Am I correctly implementing the thermal diffusion term with the divergence and gradient operators?
- Is there any better way to define the boundaries?
- Are there any known limitations or specific requirements when implementing spatially-distributed temperature in PyBaMM?
Any guidance on correctly implementing thermal dynamics in a DFN model would be greatly appreciated.
@brosaplanella @rtimms
btw these are the boundary conditions of the model:
and this is the code am implimenting:
def __init__(self, name="Doyle-Fuller-Newman model"):
super().__init__(name=name)
pybamm.citations.register("Marquis2019")
# `param` is a class containing all the relevant parameters and functions for
# this model. These are purely symbolic at this stage, and will be set by the
# `ParameterValues` class when the model is processed.
param = self.param
######################
# Variables
######################
# Variables that depend on time only are created without a domain
Q = pybamm.Variable("Discharge capacity [A.h]")
# Variables that vary spatially are created with a domain
c_e_n = pybamm.Variable(
"Negative electrolyte concentration [mol.m-3]",
domain="negative electrode",
)
c_e_s = pybamm.Variable(
"Separator electrolyte concentration [mol.m-3]",
domain="separator",
)
c_e_p = pybamm.Variable(
"Positive electrolyte concentration [mol.m-3]",
domain="positive electrode",
)
# Concatenations combine several variables into a single variable, to simplify
# implementing equations that hold over several domains
c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p)
# Electrolyte potential
phi_e_n = pybamm.Variable(
"Negative electrolyte potential [V]",
domain="negative electrode",
)
phi_e_s = pybamm.Variable(
"Separator electrolyte potential [V]",
domain="separator",
)
phi_e_p = pybamm.Variable(
"Positive electrolyte potential [V]",
domain="positive electrode",
)
phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p)
# Electrode potential
phi_s_n = pybamm.Variable(
"Negative electrode potential [V]", domain="negative electrode"
)
phi_s_p = pybamm.Variable(
"Positive electrode potential [V]",
domain="positive electrode",
)
# Particle concentrations are variables on the particle domain, but also vary in
# the x-direction (electrode domain) and so must be provided with auxiliary
# domains
c_s_n = pybamm.Variable(
"Negative particle concentration [mol.m-3]",
domain="negative particle",
auxiliary_domains={"secondary": "negative electrode"},
)
c_s_p = pybamm.Variable(
"Positive particle concentration [mol.m-3]",
domain="positive particle",
auxiliary_domains={"secondary": "positive electrode"},
)
# Temperature model
T = pybamm.Variable("Temperature [K]") #
# T = pybamm.Variable("Temperature [K]", domain=["negative electrode", "separator", "positive electrode"])
T_n = pybamm.Variable("Negative electrode temperature [K]", domain="negative electrode")
T_s = pybamm.Variable("Separator temperature [K]", domain="separator")
T_p = pybamm.Variable("Positive electrode temperature [K]", domain="positive electrode")
T_s_n = pybamm.Variable(
"Negative particle temperature [K]",
domain="negative particle",
auxiliary_domains={"secondary": "negative electrode"},
)
T_s_p = pybamm.Variable(
"Positive particle temperature [K]",
domain="positive particle",
auxiliary_domains={"secondary": "positive electrode"},
)
T= pybamm.concatenation(T_n, T_s, T_p)
######################
# Other set-up
######################
# Current density
i_cell = param.current_density_with_time
# Porosity
# Primary broadcasts are used to broadcast scalar quantities across a domain
# into a vector of the right shape, for multiplying with other vectors
eps_n = pybamm.PrimaryBroadcast(
pybamm.Parameter("Negative electrode porosity"), "negative electrode"
)
eps_s = pybamm.PrimaryBroadcast(
pybamm.Parameter("Separator porosity"), "separator"
)
eps_p = pybamm.PrimaryBroadcast(
pybamm.Parameter("Positive electrode porosity"), "positive electrode"
)
eps = pybamm.concatenation(eps_n, eps_s, eps_p)
# Active material volume fraction (eps + eps_s + eps_inactive = 1)
eps_s_n = pybamm.Parameter("Negative electrode active material volume fraction")
eps_s_p = pybamm.Parameter("Positive electrode active material volume fraction")
# transport_efficiency
tor = pybamm.concatenation(
eps_n**param.n.b_e, eps_s**param.s.b_e, eps_p**param.p.b_e
)
a_n = 3 * param.n.prim.epsilon_s_av / param.n.prim.R_typ
a_p = 3 * param.p.prim.epsilon_s_av / param.p.prim.R_typ
# Interfacial reactions
# Surf takes the surface value of a variable, i.e. its boundary value on the
# right side. This is also accessible via `boundary_value(x, "right")`, with
# "left" providing the boundary value of the left side
c_s_surf_n = pybamm.surf(c_s_n)
sto_surf_n = c_s_surf_n / param.n.prim.c_max
j0_n = param.n.prim.j0(c_e_n, c_s_surf_n, T_n)
eta_n = phi_s_n - phi_e_n - param.n.prim.U(sto_surf_n, T_n)
Feta_RT_n = param.F * eta_n / (param.R * T_n)
j_n = 2 * j0_n * pybamm.sinh(param.n.prim.ne / 2 * Feta_RT_n)
c_s_surf_p = pybamm.surf(c_s_p)
sto_surf_p = c_s_surf_p / param.p.prim.c_max
j0_p = param.p.prim.j0(c_e_p, c_s_surf_p, T_p)
eta_p = phi_s_p - phi_e_p - param.p.prim.U(sto_surf_p, T_p)
Feta_RT_p = param.F * eta_p / (param.R * T_p)
j_s = pybamm.PrimaryBroadcast(0, "separator")
j_p = 2 * j0_p * pybamm.sinh(param.p.prim.ne / 2 * Feta_RT_p)
a_j_n = a_n * j_n
a_j_p = a_p * j_p
a_j = pybamm.concatenation(a_j_n, j_s, a_j_p)
######################
# State of Charge
######################
I = param.current_with_time
# The `rhs` dictionary contains differential equations, with the key being the
# variable in the d/dt
self.rhs[Q] = I / 3600
# Initial conditions must be provided for the ODEs
self.initial_conditions[Q] = pybamm.Scalar(0)
######################
# Particles
######################
# The div and grad operators will be converted to the appropriate matrix
# multiplication at the discretisation stage
N_s_n = -param.n.prim.D(c_s_n, T_s_n) * pybamm.grad(c_s_n)
N_s_p = -param.p.prim.D(c_s_p, T_s_p) * pybamm.grad(c_s_p)
self.rhs[c_s_n] = -pybamm.div(N_s_n)
self.rhs[c_s_p] = -pybamm.div(N_s_p)
# Boundary conditions must be provided for equations with spatial derivatives
self.boundary_conditions[c_s_n] = {
"left": (pybamm.Scalar(0), "Neumann"),
"right": (
-j_n / (param.F * pybamm.surf(param.n.prim.D(c_s_n, T_s_n))),
"Neumann",
),
}
self.boundary_conditions[c_s_p] = {
"left": (pybamm.Scalar(0), "Neumann"),
"right": (
-j_p / (param.F * pybamm.surf(param.p.prim.D(c_s_p, T_s_p))),
"Neumann",
),
}
self.initial_conditions[c_s_n] = param.n.prim.c_init
self.initial_conditions[c_s_p] = param.p.prim.c_init
######################
# Current in the solid
######################
sigma_eff_n = param.n.sigma(T_n) * eps_s_n**param.n.b_s
i_s_n = -sigma_eff_n * pybamm.grad(phi_s_n)
sigma_eff_p = param.p.sigma(T_p) * eps_s_p**param.p.b_s
i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
# The `algebraic` dictionary contains differential equations, with the key being
# the main scalar variable of interest in the equation
# multiply by Lx**2 to improve conditioning
self.algebraic[phi_s_n] = param.L_x**2 * (pybamm.div(i_s_n) + a_j_n)
self.algebraic[phi_s_p] = param.L_x**2 * (pybamm.div(i_s_p) + a_j_p)
self.boundary_conditions[phi_s_n] = {
"left": (pybamm.Scalar(0), "Dirichlet"),
"right": (pybamm.Scalar(0), "Neumann"),
}
self.boundary_conditions[phi_s_p] = {
"left": (pybamm.Scalar(0), "Neumann"),
"right": (i_cell / pybamm.boundary_value(-sigma_eff_p, "right"), "Neumann"),
}
# Initial conditions must also be provided for algebraic equations, as an
# initial guess for a root-finding algorithm which calculates consistent initial
# conditions
self.initial_conditions[phi_s_n] = pybamm.Scalar(0)
self.initial_conditions[phi_s_p] = param.ocv_init
######################
# Current in the electrolyte
######################
#code initial i_e
# i_e = (param.kappa_e(c_e, T) * tor) * (
# param.chiRT_over_Fc(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
# )
# Domain-specific transport efficiencies
tor_n = eps_n**param.n.b_e # negative electrode
tor_s = eps_s**param.s.b_e # separator
tor_p = eps_p**param.p.b_e # positive electrode
tor = pybamm.concatenation(tor_n, tor_s, tor_p)
# Current in each domain
i_e_n = (param.kappa_e(c_e_n, T_n) * tor_n) * (
param.chiRT_over_Fc(c_e_n, T_n) * pybamm.grad(c_e_n) - pybamm.grad(phi_e_n)
)
i_e_s = (param.kappa_e(c_e_s, T_s) * tor_s) * (
param.chiRT_over_Fc(c_e_s, T_s) * pybamm.grad(c_e_s) - pybamm.grad(phi_e_s)
)
i_e_p = (param.kappa_e(c_e_p, T_p) * tor_p) * (
param.chiRT_over_Fc(c_e_p, T_p) * pybamm.grad(c_e_p) - pybamm.grad(phi_e_p)
)
i_e = pybamm.concatenation(i_e_n, i_e_s, i_e_p)
# multiply by Lx**2 to improve conditioning
# Current divergence in each domain
div_i_e_n = pybamm.div(i_e_n)
div_i_e_s = pybamm.div(i_e_s)
div_i_e_p = pybamm.div(i_e_p)
# Concatenate the divergences
div_i_e = pybamm.concatenation(div_i_e_n, div_i_e_s, div_i_e_p)
self.algebraic[phi_e] = param.L_x**2 * (div_i_e - a_j)
self.boundary_conditions[phi_e] = {
"left": (pybamm.Scalar(0), "Neumann"),
"right": (pybamm.Scalar(0), "Neumann"),
}
self.initial_conditions[phi_e] = -param.n.prim.U_init
######################
# Electrolyte concentration
######################
N_e = -tor * param.D_e(c_e, T) * pybamm.grad(c_e)
self.rhs[c_e] = (1 / eps) * (
-pybamm.div(N_e) + (1 - param.t_plus(c_e, T)) * a_j / param.F
)
self.boundary_conditions[c_e] = {
"left": (pybamm.Scalar(0), "Neumann"),
"right": (pybamm.Scalar(0), "Neumann"),
}
self.initial_conditions[c_e] = param.c_e_init
######################
# (Some) variables
######################
voltage = pybamm.boundary_value(phi_s_p, "right")
num_cells = pybamm.Parameter(
"Number of cells connected in series to make a battery"
)
######################
# Temperature model
######################
#1. Irreversible heat
# Original equation: q_ir = j_Li(φ_s - φ_e - U)
q_ir_p = j_p * eta_p #positive
q_ir_n = j_n * eta_n #negative
q_ir_s = pybamm.PrimaryBroadcast(0, "separator")
q_ir = pybamm.concatenation(q_ir_n, q_ir_s, q_ir_p)
#2. Ohmic heat
# Original equation: q_o = σeff∇φs∇φs + κeff∇φe∇φe
#Solid phase ohmic heat
q_ohm_s_n = i_s_n * pybamm.grad(phi_s_n) # negative
q_ohm_s_p = i_s_p * pybamm.grad(phi_s_p) # positive
q_ohm_s_s = pybamm.PrimaryBroadcast(0, "separator")
# Electrolyte phase ohmic heat par domaine
# Electrolyte phase ohmic heat by domain
q_ohm_e_n = i_e_n * pybamm.grad(phi_e_n) # negative
q_ohm_e_s = i_e_s * pybamm.grad(phi_e_s) # separator
q_ohm_e_p = i_e_p * pybamm.grad(phi_e_p) # positive
# Concatenation pour chaque phase
q_ohm_s = pybamm.concatenation(q_ohm_s_n, q_ohm_s_s, q_ohm_s_p)
q_ohm_e = pybamm.concatenation(q_ohm_e_n, q_ohm_e_s, q_ohm_e_p)
# Total ohmic heat
q_o = q_ohm_s + q_ohm_e
#3. Ionic heating:
# Original equation: q_i = κD_eff∇ln(ce)∇φe
# as mentionned in i calculation : kappa_D_eff = (param.kappa_e(c_e, T) * tor) * param.chiRT_over_Fc(c_e, T)
# calculate κD_eff for each domain
kappa_D_eff_n = (param.kappa_e(c_e_n, T_n) * eps_n**param.n.b_e) * param.chiRT_over_Fc(c_e_n, T_n)
kappa_D_eff_s = (param.kappa_e(c_e_s, T_s) * eps_s**param.s.b_e) * param.chiRT_over_Fc(c_e_s, T_s)
kappa_D_eff_p = (param.kappa_e(c_e_p, T_p) * eps_p**param.p.b_e) * param.chiRT_over_Fc(c_e_p, T_p)
# Calculate ionic heat
q_i_n = kappa_D_eff_n * pybamm.grad(pybamm.log(c_e_n)) * pybamm.grad(phi_e_n)
q_i_s = kappa_D_eff_s * pybamm.grad(pybamm.log(c_e_s)) * pybamm.grad(phi_e_s)
q_i_p = kappa_D_eff_p * pybamm.grad(pybamm.log(c_e_p)) * pybamm.grad(phi_e_p)
# Concatenate for full domain
q_i = pybamm.concatenation(q_i_n, q_i_s, q_i_p)
# 4. Reversible entropic heat (qre)
# Original equation: q_re = j_Li * T * (∂U/∂T)
# Reversible heat in negative electrode
dUdT_n = param.n.prim.dUdT(sto_surf_n) # ∂U/∂T for negative
q_re_n = j_n * T_n * dUdT_n
# Reversible heat in positive electrode
dUdT_p = param.p.prim.dUdT(sto_surf_p) # ∂U/∂T for positive
q_re_p = j_p * T_p * dUdT_p
# Zero in separator (no reactions)
q_re_s = pybamm.PrimaryBroadcast(0, "separator")
# Total reversible heat across all domains
q_re = pybamm.concatenation(q_re_n, q_re_s, q_re_p)
# 5. Contact resistance heat (qc)
#original eq : qc= (I^2 R_c)/ ( Av)
R_c = pybamm.Parameter("Contact resistance [Ohm]")
h = pybamm.Parameter("Electrode height [m]")
w = pybamm.Parameter("Electrode width [m]")
v = pybamm.Parameter("Cell volume [m3]")
Av = h * w * v
q_c_scalar = (I**2 * R_c) / Av
# q_c calculated as a single scalar value so to add it to other heat terms, it needs to be broadcasted
q_c_n = pybamm.PrimaryBroadcast(q_c_scalar, "negative electrode")
q_c_s = pybamm.PrimaryBroadcast(q_c_scalar, "separator")
q_c_p = pybamm.PrimaryBroadcast(q_c_scalar, "positive electrode")
q_c = pybamm.concatenation(q_c_n, q_c_s, q_c_p)
Q1 = q_ir + q_o + q_i + q_re + q_c
#Now heat equation:
# Density
rho_n = pybamm.PrimaryBroadcast(
pybamm.Parameter("Negative electrode density [kg.m-3]"),
"negative electrode"
)
rho_s = pybamm.PrimaryBroadcast(
pybamm.Parameter("Separator density [kg.m-3]"),
"separator"
)
rho_p = pybamm.PrimaryBroadcast(
pybamm.Parameter("Positive electrode density [kg.m-3]"),
"positive electrode"
)
rho = pybamm.concatenation(rho_n, rho_s, rho_p)
# Specific heat capacity
cp_n = pybamm.PrimaryBroadcast(
pybamm.Parameter("Negative electrode specific heat capacity [J.kg-1.K-1]"),
"negative electrode"
)
cp_s = pybamm.PrimaryBroadcast(
pybamm.Parameter("Separator specific heat capacity [J.kg-1.K-1]"),
"separator"
)
cp_p = pybamm.PrimaryBroadcast(
pybamm.Parameter("Positive electrode specific heat capacity [J.kg-1.K-1]"),
"positive electrode"
)
cp = pybamm.concatenation(cp_n, cp_s, cp_p)
# Thermal conductivity
lambda_n = pybamm.PrimaryBroadcast(
pybamm.Parameter("Negative electrode thermal conductivity [W.m-1.K-1]"),
"negative electrode"
)
lambda_s = pybamm.PrimaryBroadcast(
pybamm.Parameter("Separator thermal conductivity [W.m-1.K-1]"),
"separator"
)
lambda_p = pybamm.PrimaryBroadcast(
pybamm.Parameter("Positive electrode thermal conductivity [W.m-1.K-1]"),
"positive electrode"
)
lambdac = pybamm.concatenation(lambda_n, lambda_s, lambda_p)
N = - lambdac * pybamm.grad(T)
dTdt = (1 / ( rho * cp )) * ( - pybamm.div(N) + Q1)
# self.rhs = {T : dTdt}
self.rhs[T] = dTdt
# Initial conditions
T_init = pybamm.Parameter("Initial temperature [K]")
# # self.initial_conditions[T] = T_init
T_init_n = pybamm.PrimaryBroadcast(T_init, "negative electrode")
T_init_s = pybamm.PrimaryBroadcast(T_init, "separator")
T_init_p = pybamm.PrimaryBroadcast(T_init, "positive electrode")
T_init_combined = pybamm.concatenation(T_init_n, T_init_s, T_init_p)
self.initial_conditions[T] = T_init_combined
# Boundary conditions for temperature
h = pybamm.Parameter("Total heat transfer coefficient [W.m-2.K-1]")
T_amb = pybamm.Parameter("Ambient temperature [K]")
self.boundary_conditions[T] = {
"left": (h * (T_amb - T), "Neumann"), # Convective cooling at surface
"right": (0, "Neumann") # Insulated at other end
}
# The `variables` dictionary contains all variables that might be useful for
# visualising the solution of the model
self.variables = {
"Negative particle concentration [mol.m-3]": c_s_n,
"Negative particle surface concentration [mol.m-3]": c_s_surf_n,
"Electrolyte concentration [mol.m-3]": c_e,
"Negative electrolyte concentration [mol.m-3]": c_e_n,
"Separator electrolyte concentration [mol.m-3]": c_e_s,
"Positive electrolyte concentration [mol.m-3]": c_e_p,
"Positive particle concentration [mol.m-3]": c_s_p,
"Positive particle surface concentration [mol.m-3]": c_s_surf_p,
"Current [A]": I,
"Current variable [A]": I, # for compatibility with pybamm.Experiment
"Negative electrode potential [V]": phi_s_n,
"Electrolyte potential [V]": phi_e,
"Negative electrolyte potential [V]": phi_e_n,
"Separator electrolyte potential [V]": phi_e_s,
"Positive electrolyte potential [V]": phi_e_p,
"Positive electrode potential [V]": phi_s_p,
"Voltage [V]": voltage,
"Battery voltage [V]": voltage * num_cells,
"Time [s]": pybamm.t,
"Discharge capacity [A.h]": Q,
"Temperature [K]": T,
"Temperature [°C]": T - 273.15,
}
# Events specify points at which a solution should terminate
self.events += [
pybamm.Event("Minimum voltage [V]", voltage - param.voltage_low_cut),
pybamm.Event("Maximum voltage [V]", param.voltage_high_cut - voltage),
]
@property
def default_parameter_values(self):
return pybamm.ParameterValues("Chayambuka2022")