Choosing model options, adding modules and parameter

[kabravenkatesh]

Hi everyone,
I am trying to understand few basic details of the pybamm model. So, I am trying to simulate the Li-ion battery with the help of a DFN model. I have two questions:

Within the pybamm framework there are multiple combination of options that can be chosen for a model.

1. How do I know what are the possible 'options' choices? I mean is there a place where I can find a documentation for them? Do I need to look into the source code?

model = pybamm.lithium_ion.DFN(build=False, options = {
"SEI": "reaction limited", 'surface form':'differential'
}
)

2. I want to introduce another equation and couple it with the existing set of equations. For example I wanted to change the temperature from a lumped model to a spatial one. First question is, How do I solve the energy conservation equation within the pybamm framework (or for that matter any new set of equation) within the pybamm framework. And secondly how do I couple the dependence of the electrochemical variables onto the localized temperature field.

3. How do I include an additional dependence for a parameter. For example, Ai2020 dataset has following relation for electrolyte diffusivity:

def De(c_e, T):
D_c_e = 10 ** (-8.43 - 54 / (T - 229 - 5e-3 * c_e) - 0.22e-3 * c_e)
return D_c_e

Hypothetically, let us say I want to modify the dependance of the electrolyte diffusivity to include the solid phase Li concentration as well.

def modified_De(c_e, c_s, T):
D_c_e = mult*10 ** (-8.43 - 54 / (T - 229 - 5e-3 * c_e) - 0.22e-3 * c_e)*c_s/cs_max
return D_c_e

Apart from defining above function and changing to parameter. How do I make such a change in the code?

Thanks in advance

[kabravenkatesh]

@brosaplanella @tinosulzer @DrSOKane @rtimms Any help would be greatly appreciated. TIA!

[ChaddKiggins]

To help with a few of the questions:

1. Avilable model options can be found here: https://docs.pybamm.org/en/latest/source/api/models/base_models/base_battery_model.html

2. PyBaMM has a few thermal model options. This link has a description of the available thermal submodels:
   https://docs.pybamm.org/en/stable/source/examples/notebooks/models/thermal-models.html

For spatial models, the options above work with pouch cells in certain dimensions. If you have another format, say cylindrical, you may need to use an external model. For example, to solve the spatial temperture distribution in a cylindrical cell, I used PyBaMM to calculate the total heating of the DFN model. I then took this heating term and coupled it to an external Finite Difference model to look at how the spatial temperature distribution evolved over charge/discharge. I also discretized the cylindrical cell into multiple DFN models so local DFN models would be dependant on the local temperature (cold spots would behave different than hot spots).

This link may also give some insight in how to do more in the PyBaMM environment: https://docs.pybamm.org/en/stable/source/examples/notebooks/models/unsteady-heat-equation.html

3. This link gives a good idea of what input parameters each Lithium Ion parameter expects:
   https://docs.pybamm.org/en/latest/_modules/pybamm/parameters/lithium_ion_parameters.html#LithiumIonParameters

To change what a Lithium Ion Parameter is dependent on, I think you would need to modify the source code. Below is the code for the electrolyte diffusivity.

  def D_e(self, c_e, T):
      """Dimensional diffusivity in electrolyte"""
      tol = pybamm.settings.tolerances["D_e__c_e"]
      c_e = pybamm.maximum(c_e, tol)
      inputs = {"Electrolyte concentration [mol.m-3]": c_e, "Temperature [K]": T}
      return pybamm.FunctionParameter("Electrolyte diffusivity [m2.s-1]", inputs)

[kabravenkatesh]

Hi, @ChaddKiggins This is super helpful. Thanks a lot for your response!

I wish to know how are you going about coupling the heat generation. Is it bi-directionally coupled so that the temperature field affects the electrochemical parameters or a one-way coupling that only considering heat generation from the model?

[ChaddKiggins]

No worries at all!

It's bi-directional so the thermal model will increase the temperature and this increased temperature will inform properties such as diffusivity, exchange current, etc. As an example, here is a 1D case where I looked at the temperature distribution in the radial direction of a cylindrical cell.

Here, I am roughly assuming I have 30 winds of the jelly roll from the core to the exterior. For each wind I am going to have a local P2D model to simulate the electrochemistry in that local position. So I'll have an array of 30 P2D's.

Nr = 30 #Number of nodes (or winds) for my jelly roll in radial direction
p2ds = [] #List to hold all the P2D models
for _ in range(Nr):
    model = P2dnode(parameters = get_parameter_values()) # P2D node is Internal Class
    model.parameters.update({"Current function [A]": "[input]", "Ambient temperature [K]": "[input]", "Initial temperature [K]":"[input]"})
    model.parameters.set_initial_stoichiometries(1)
    model.setup_model()
    p2ds.append(model)

Then I'll solve each of these 30 models over a given time step. Here the total simulation time is set to 600 seconds but I'm solving in 10 second steps. After this 10 second solve, I'll get the average total heating (Q) from each of the 30 P2D nodes. I'll then pass this heating and initial temperature into the external thermal model.

#Setup Time intervals for stepping the solution
ti = 0
dt_bat = 10 #P2D simulation will solve for 10 seconds
tf = 600 #Run sim for 720 seconds or until break condition is met

T = np.ones(Nr) * T_in #Initial temperature array for the sim
while ti < tf:
    #Step 1: Solve P2D models
    Q = []
    for i, p2d in enumerate(p2ds):
        solver = pybamm.CasadiSolver(mode="safe", atol=1e-06, rtol=1e-06, dt_max=10)
        inputs = {"Current function [A]": 20, "Ambient temperature [K]": T[i], "Initial temperature [K]": T[i]}
        p2d.step_solution = solver.step(model=p2d.model, old_solution=p2d.step_solution, npts=5, dt=dt_bat, inputs=inputs, save=True)
        print(p2d.step_solution["Ambient temperature [K]"].data[-1])
        q = p2d.step_solution["Total heating [W.m-3]"].data
        q = np.average(q)
        Q.append(q)
    Q = np.array(Q)
    print(Q)
    print(T)

    solver = OneDThermSolver_Cylindral(Nr, R, t_final, dt, rho, cp, k, h, T_amb, Q, T) #External thermal model
    T = solver.backward_euler()

    print(T)

    if p2d.step_solution["Terminal voltage [V]"].data[-1] < 2.51:
        break

    ti += dt_bat

print(T)

solver.plot_1D()

The thermal model is using Backwards Euler for time stepping and Centered Differences for spatial discretization. It is currently setup with a Neumann boundary at the "jelly roll" core for a no flux condition, and a robin boundary at the jelly-roll exterior for convection.

class OneDThermSolver_Cylindral:
    def __init__(self, Nr, R, t_final, dt, rho, cp, k, h, T_amb, Q, T_in, can = False):
        .....
        
    def backward_euler(self):

        #Setup 
        #Representation of sparse matrix and right hand side
        diagonal = np.zeros(self.Nr)
        lower = np.zeros(self.Nr-1)
        upper = np.zeros(self.Nr-1)
        b = np.zeros(self.Nr)

        #Precompute the sparse matrix (2nd order accurate)
        diagonal[:] = 1 + 2*self.F
        lower[:] =  self.F *self.dr/(self.r[1:]*2) - self.F
        upper[:] = -self.F * self.dr/(self.r[:-1]*2) - self.F
    
        A = scipy.sparse.diags(
        diagonals=[diagonal, lower, upper],
        offsets = [0, -1, 1],
        shape=(self.Nr,self.Nr),
        format="lil") #Set to "lil" so we can insert boundary conditions

        #Insert left axisymmetric boundary condition representing no flux 
        A[0, 1] = -8*self.F
        A[0, 2] = self.F

        #Insert right convective robin boundary condition
        A[-1, -3] = self.F/2
        A[-1, -2] = -4*self.F
        A[-1, -1] = 1 + 3*self.dr*self.F*self.h/self.k + self.dt*self.h*self.alpha/(self.R*self.k) +7/2*self.F

        #Set sparse matrix to "csr" for solving
        A = A.tocsr()

        #Setup initial temperature array
        self.T = np.ones(self.Nr) * self.T_in
        #T = np.zeros(self.Nr)

        for t in range(0, self.time_steps - 1):
            #Create b matrix interior nodes
            b[1:-1] = self.T[1:-1] + self.S[1:-1]*self.dt #Add the source term to the initial temps
            #Insert left boundary condition b term
            b[0] = self.T[0] + self.S[0]*self.dt
            #Insert right boundary condition b term
            b[-1] = self.T[-1] + 3*self.dr*self.h*self.F*self.T_amb/self.k + self.dt *self.h*self.alpha*self.T_amb/(self.R*self.k) + self.S[-1] * self.dt
            #Solve
            self.T[:] = scipy.sparse.linalg.spsolve(A, b)

        return self.T

After each time step, the output temperature of the thermal model will be the input temperature to the next round of sims for the 30 P2Ds. As the temperature starts to diverge (due to the boundary conditions), each of the 30 P2D models will have a different initial temperature, therefore their local electrochemistry will be different due to the Arrhenius dependence on many of the parameters, for example the exchange current densities.

I'm sure there are more elegant ways to do this but this works well to get a spatial temperature distribution as well as see how the local electrochemistry diverges as the temperature distribution diverges. It's also easy to extend to 2D/3D.

Let me know and I'd be happy to share the full code. Didn't want to write it all here as I didn't want to completely bog down the conversation lol

[TomTranter]

To add to @ChaddKiggins great answer I worked on a similar problem a few years ago now and published some papers on coupling pybamm models with thermal models for jellyrolls and extended the work to packs which became liionpack .

https://iopscience.iop.org/article/10.1149/1945-7111/aba44b/meta
https://iopscience.iop.org/article/10.1149/1945-7111/abd44f/meta
https://joss.theoj.org/papers/10.21105/joss.04051

My method was similar to Chad's where I step and solve the electrochemistry then update the temperature distribution then feed that back and calculate new heat sources. I got the solving parallelized so that I was running a CC 1C discharge with about 1300 models around a 2D jellyroll on my desktop with 12 cores in about 10 minutes. I was taking steps of about 30 s I think. The code hasn't been updated for a while but it's getting more interest now and I will be updating it soon. Currently I haven't particularly advertised that it's out there and open source as I don't have much time to support it.

Also I think the solving speed could be greatly improved by doing away with the step solve and the long term plan is to actually migrate all the solving including the thermal to PyBaMM proper and just use this package to generate the meshes.

https://github.com/TomTranter/pybamm_pnm/tree/main