Adding a new field equation to DFN: coupling mechanics with electrochemistry

PyBaMM Q&A
1

Hi all,

I’m trying to add a custom equation to the existing DFN model in PyBaMM, but I’m struggling to figure out how to properly integrate it into the framework so that it’s solved simultaneously with the rest of the model equations.

My use case involves coupling electrochemical and mechanical effects, specifically the two-way interaction between the particle domain and the electrode domain. To start, I want to add an equation describing displacement in the electrode due to particle swelling, and later feed it back into the model via hydrostatic pressure, which influences particle dynamics. For now, I’m working in one dimension only—along the electrode thickness—to keep things manageable.

The initial equation I’d like to implement is:

div(grad(u)) = (Ω/3) * grad(c_s_avg)

Here:
- u is the displacement field I want to solve for,
- Ω is a constant parameter,
- c_s_avg is the average particle concentration in the electrode control volume
(i.e., “R-averaged negative particle concentration [mol.m-3]” in PyBaMM notation),
- grad and div are the standard gradient and divergence operators in one dimension.

I’ve successfully calculated the displacement field in post-processing (which was straightforward), but I now plan to compute strain and stress from this displacement, and use the stress to influence the particle surface flux. This creates a feedback loop that makes the system implicitly coupled, and I’d like all equations to be solved together.

However, I haven’t figured out how to inject this additional equation into the existing DFN model. I’ve reviewed the following tutorials from https://docs.pybamm.org/en/latest/source/examples/notebooks/creating_models/:

  • Half-cell model example
  • Simple SEI model
  • Creating a submodel

These examples are helpful, but they all build models from scratch using the BaseModel class. What I’d really like to do is extend the DFN model—either by directly adding my algebraic equation or by wrapping it into a custom submodel and plugging it into the DFN structure.

Is there a recommended way to do this? Any tips on how to proceed, or pointers to similar examples involving the DFN or another built-in model, or at least an explanation of what the error ‘negative particle’ means, would be much appreciated.

Here’s a snippet of the code I’ve been working with together with the error that I am getting (just let me know if you’d like the full block):

model = pybamm.lithium_ion.DFN()

#defining a submodel class with my equation
class MySubmodel(pybamm.BaseSubModel):
    def __init__(self, param, domain, options=None):
        super().__init__(param, domain, options=options)

    # def get_fundamental_variables(self):
    #     pass 

    def get_coupled_variables(self, variables):
        # create displacement variable
        u = pybamm.Variable("Displacement [m]", domain="negative electrode")

        variables = {
            "Displacement [m]": u,
        }

        return variables
    
    def set_rhs(self, variables):
        pass

    def set_algebraic(self, variables):
        # extract the variables we need
        u = variables["Displacement [m]"]
        c_s_rav = variables["R-averaged negative particle concentration [mol.m-3]"]
        
        Omega = self.param["Negative electrode partial molar volume [m3.mol-1]"]

        # Define the algebraic equation
        ddu_dxx = pybamm.div(pybamm.grad(u)) - Omega / 3 * pybamm.grad(c_s_rav)

        # add it to the submodel dictionary
        self.algebraic = {u:ddu_dxx}

    def set_boundary_conditions(self, variables):
        # extract the variables we need
        u = variables["Displacement [m]"]
        e0_right = pybamm.Scalar(-0.024)

        # add the boundary conditions to the submodel dictionary
        self.boundary_conditions = {
            u: {"left": (pybamm.Scalar(0), "Dirichlet"),  
                "right": (e0_right, "Neumann")}  
        }

    def set_initial_conditions(self, variables):
        # extract the variable we need
        u = variables["Displacement [m]"]

        # define the initial displacement parameter
        u0 = pybamm.Scalar(0)

        self.initial_conditions = {u: u0}
        
model.submodels = {
    "MyMechanics": MySubmodel(None, "Negative"),
}

model.update()


#model parameters, geometry, mesh ...
param = pybamm.ParameterValues("Ai2020")

x_i = pybamm.SpatialVariable("xi", domain="negative electrode", coord_sys="cartesian")
L_n = param["Negative electrode thickness [m]"]

geometry = {
    "negative electrode": {x_i: {"min": 0, "max": L_n}},
}

spatial_methods = {"negative electrode": pybamm.FiniteVolume()}

submesh_types = {"negative electrode": pybamm.Uniform1DSubMesh}

var_pts = {x_i: 20} 

mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

solver = pybamm.IDAKLUSolver()

sim = pybamm.Simulation(
    model,
    geometry=geometry,
    parameter_values=param,
    submesh_types=submesh_types,
    var_pts=var_pts,
    spatial_methods=spatial_methods,
    solver=solver,
)

#solve the equations
t = [0, 100] 
solution = sim.solve(model, t)

Output:

KeyError: 'negative particle'