Thermal DFN Model Implementation Error

[Maguixx]

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 a ShapeError: pybamm.expression_tree.exceptions.ShapeError: Cannot find shape (original error: inconsistent shapes)​

Questions

1. Is there an issue with how I'm defining the domains for the temperature variables?
2. Am I correctly implementing the thermal diffusion term with the divergence and gradient operators?
3. 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 @tinosulzer

This is the code i am implimenting:
import pybamm class ThermalDFN(pybamm.lithium_ion.BaseModel): """Doyle-Fuller-Newman (DFN) model of a sodium-ion battery, from :footcite:t:Marquis2019`.

Parameters
----------
name : str, optional
    The name of the model.

"""

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_n] = {
        "left": (h * (T_amb - T_n), "Neumann"),
        "right": (0, "Neumann")
    }
    # 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")

`

[github-actions[bot]]

Hi! We have now moved our discussions to Discourse. Please post your question there.

[valentinsulzer]

This looks like an issue with boundary conditions. You need boundary conditions for T, not T_n.

[Maguixx]

Thank you for your response. Upon implementing your suggestion to replace T_n with T, I encountered a dimensional inconsistency error, since i have these as boundary conditions

. Specifically, a ValueError: dimension mismatch occurs during the temperature gradient calculation (grad(Temperature [K])). This suggests that when the finite volume method attempts to apply boundary conditions, it encounters matrices with incompatible dimensions.
Could you advise on how to properly define the temperature variable, boundary conditions to ensure dimensional consistency in the gradient calculations?