Nearly (but not quite) physical results from a custom model - error in my use of FunctionParameters?

[isaacbasil]

Hi everyone,

I am programming a custom model for a Li-ion cell in Pybamm, which now solves without throwing any errors (model details below). The results are not absurd, but something is definitely wrong - I have attached some pictures of results compared to the DFN. I suspect that the problem might be coming from an error in the way I am defining one of the functions in the model. I have written a function, "j_ave" which is essentially the Butler Volmer equation. This function calls 2 other functions, one for exchange current density "i0" and one for open circuit potential "u0". I have tried:

1. Defining j_ave, i0 and u0 all as FunctionParameters.
2. Defining i0 and u0 as FunctionParameters, and j_ave just as a regular function in the model definition.
3. Defining j_ave, i0 and u0 all as regular functions.

The reason I think that j_ave is not outputting correct results is because it seems to be completely insensitive to any change. For example, I can add a factor of 10^6 to the output of j_ave (multiplying the electrochemical reaction rate by 10^6) and it hardly changes the results. In the second and third plot that I have attached, you can see that the reaction rate differs from the DFN by about a factor of 10 in both electrodes.

Any recommendations on what I may be doing wrong or on other approaches I could try for writing this j_ave function would be much appreciated :)

Model description

You can think of the model as similar to the DFN model, but instead of solving in an extra "pseudo" dimension for concentration within active material, the active material phase is also homogenised (like the other 3 governing equations). Thus, the model is completely macroscale and 1D. The electrolyte and active material phase equations are coupled by the electrochemical reaction source term, j_ave.
In the plots below, the dotted green line is the model that I am talking about here, compared to the DFN using the exact same parameter values. (please ignore the blue line called LIBAT).

Here are the most important parts of the model definition in pybamm, with the j_ave function below:

# calculate reaction rates
a_j_n = av_n * self.j_ave("Negative", c_e_n, c_s_n, c_n_max, phi_e_n, phi_s_n,
                          self.T)  # reaction rate for negative electrode

j_s = pybamm.PrimaryBroadcast(0, "separator") #reaction rate for separator

a_j_p = av_p * self.j_ave("Positive", c_e_p, c_s_p, c_p_max, phi_e_p, phi_s_p,
                          self.T)  # reaction rate for positive electrode

a_j = pybamm.concatenation(a_j_n, j_s, a_j_p)

# calculate current
i_n = - eps_s_n * (sigma_n / tau_s_n) * pybamm.grad(phi_s_n)  # negative electrode volume averaged current
i_e = (eps_e / tau_e) * (- sigma_e * pybamm.grad(phi_e) - theta * 1 / c_e * pybamm.grad(c_e))
i_p = - eps_s_p * (sigma_p / tau_s_p) * pybamm.grad(phi_s_p)  # negative electrode (macroscale) current

# calculate flux
N_s_n = - (eps_s_n * D_s_n / tau_s_n) * pybamm.grad(c_s_n)  # negative AM flux
N_e = - (eps_e * D_e / tau_e) * pybamm.grad(c_e) # electrolyte flux
N_s_p = - (eps_s_p * D_s_p / tau_s_p) * pybamm.grad(c_s_p)  # negative AM flux

# define negative electrode mass balance
dc_s_n_dt = 1 / eps_s_n * (
            - pybamm.div(N_s_n) - a_j_n)  # define the rhs equation
self.model.rhs = {c_s_n: dc_s_n_dt}  # add the equation to rhs dictionary

# define electrolyte mass balance
dc_e_dt = 1 / eps_e * (-pybamm.div(N_e) + (
            1 - transp_no) * a_j)  # define the rhs equation
self.model.rhs[c_e] = dc_e_dt  # add the equation to rhs dictionary

# define positive electrode mass balance
dc_s_p_dt = 1 / eps_s_p * (
            - pybamm.div(N_s_p) - a_j_p)  # define the rhs equation
self.model.rhs[c_s_p] = dc_s_p_dt  # add the equation to rhs dictionary

# define negative electrode charge balance
self.model.algebraic[phi_s_n] = L_x ** 2 * (
            - pybamm.div(i_n) - self.F * a_j_n)  

# define electrolyte charge balance
self.model.algebraic[phi_e] = L_x ** 2 * (- pybamm.div(i_e) + self.F * a_j)

# define positive electrode charge balance
self.model.algebraic[phi_s_p] = L_x ** 2 * (- pybamm.div(i_p) - self.F * a_j_p)

##### define initial conditions #####

self.model.initial_conditions[c_s_n] = c_s_n_0
self.model.initial_conditions[c_e] = c_e_0
self.model.initial_conditions[c_s_p] = c_s_p_0
# Initial conditions must also be provided for algebraic equations, as an
# initial guess for a root-finding algorithm which calculates consistent initial
# conditions
self.model.initial_conditions[phi_s_n] = pybamm.Scalar(0)
self.model.initial_conditions[phi_s_p] = 3.9832  # initial_u0_p - initial_u0_n
# initial guess for electrolyte phi as 0, this is what is used in basic DFN
self.model.initial_conditions[phi_e] = 0

##### define boundary conditions #####

self.model.boundary_conditions[c_s_n] = {
    "left": (pybamm.Scalar(0), "Neumann"),
    "right": (pybamm.Scalar(0), "Neumann"),
}

self.model.boundary_conditions[c_s_p] = {
    "left": (pybamm.Scalar(0), "Neumann"),
    "right": (pybamm.Scalar(0), "Neumann"),
}

self.model.boundary_conditions[c_e] = {
    "left": (pybamm.Scalar(0), "Neumann"),
    "right": (pybamm.Scalar(0), "Neumann"),
}

self.model.boundary_conditions[phi_s_n] = {
    "left": (pybamm.Scalar(0), "Dirichlet"),
    "right": (pybamm.Scalar(0), "Neumann"),
}

grad_phi_rhs = - (I / A) / sigma_p
self.model.boundary_conditions[phi_s_p] = {
    "left": (pybamm.Scalar(0), "Neumann"),
    "right": (grad_phi_rhs, "Neumann"),
}

self.model.boundary_conditions[phi_e] = {
    "left": (pybamm.Scalar(0), "Neumann"),
    "right": (pybamm.Scalar(0), "Neumann"),
}


def j_ave(self, key, c_e, c_s_surf, c_s_max, phi_e, phi_s, T):
    '''the butler volmer equation, which returns the surface averaged molar reaction rate'''

    def u0_graphite(c_s, c_s_max):
        """
        The u0 function for graphite
        """
        x = c_s / c_s_max

        u1x = (
                -0.01059423355572770 * pybamm.tanh((x - 0.01453708425609560) / (9.089868397988610e-5)) +
                0.02443615203087110 * pybamm.tanh((x - 0.5464261369950400) / (0.6270508166379020)) -
                0.01637520788053810 * pybamm.tanh((x - 0.5639025014475490) / (0.07053886409518520)) -
                0.06542365622896410 * pybamm.tanh((x - 0.5960370524233590) / (1.409966536648620)) -
                0.04173226059293490 * pybamm.tanh((x - 0.1787670587868640) / (0.07693844911793470)) -
                0.4792178163846890 * pybamm.tanh((x + 0.003845707852011820) / (0.04112633446959460)) -
                0.04364293924074990 * pybamm.tanh((x - 0.09449231893318330) / (-0.02046776012570780)) -
                0.08241166396760410 * pybamm.tanh((x - 0.07746685789572230) / (0.03593817905677970)) +
                0.6594735004847470
        )

        u2x = (
                (((((4.466520763420150 * x - 1.163251395483540)
                    * x - 597.3328922807720) * x + 2173.321763906580)
                  * x - 3005.282176279630) * x + 1857.779353226540)
                * x - 431.7897183333920
        )

        u_eq = u1x + (u2x - u1x) / (1.0 + pybamm.exp(-100.0 * (x - 0.8627886168611160)))

        return u_eq

    def u0_nmc(c_s, c_s_max):
        """
        The u0 function for nmc
        """
        x = c_s / c_s_max
        logX = pybamm.log(x)

        u_eq = (
                13.49050 - 10.96038 * x +
                8.203617 * pybamm.exp(1.358699 * logX) -
                3.10758 * 1e-6 * pybamm.exp(127.1216 * x - 114.2593) -
                7.033556 * pybamm.exp(-3.362749 * 1e-2 * logX)
        )

        return u_eq

    def i0_graphite(c_e, c_s_surf, c_s_max):
        """"
        A function to calculate the exchange current density in the graphite
        """

        k_n = 6.48e-7  # taken from chen2020, equivalent to Fk
        alpha_n = 0.5  # the anodic strength coeff
        ce_ref = 1  # the reference concentration

        i0_graph = k_n * c_s_surf ** (1 - alpha_n) * (c_e / ce_ref) ** alpha_n * (c_s_max - c_s_surf) ** (alpha_n)

        return i0_graph

    def i0_nmc(c_e, c_s_surf, c_s_max):
        """"
        A function to calculate the exchange current density in the nmc
        """
        k_p = 4.824e-6  # taken from Mohtat2020
        alpha_p = 0.5  # the anodic strength coeff
        ce_ref = 1  # the reference concentration
        F = 96485

        i0_nmc = k_p * c_s_surf ** (1 - alpha_p) * (c_e / ce_ref) ** alpha_p * (c_s_max - c_s_surf) ** (alpha_p)

        return i0_nmc

    if key == "Negative":
        u0 = u0_graphite(c_s_surf, c_s_max)
        i0 = i0_graphite(c_e, c_s_surf, c_s_max)
    else:
        u0 = u0_nmc(c_s_surf, c_s_max)
        i0 = i0_nmc(c_e, c_s_surf, c_s_max)

    F = 96485
    R = 8.314
    alpha = 0.5
    eta = phi_s - phi_e - u0
    j0 = i0 / F

    return j0 * (pybamm.exp(alpha * F / (R * T) * eta)
                 - pybamm.exp(- (1 - alpha) * F / (R * T) * eta)) 

[valentinsulzer]

I think your Butler-Volmer equations are fine. Try checking carefully for where you put the factors of a (as in 3*eps/R) and F as different models have different standards. Your interfacial current is off by a factor of around 10 which is what's causing the discrepancy in all the other variables (your model works fine but it thinks the current is 1/10th of what it actually is).

By the way, is the "LIBAT data" public?

[isaacbasil]

Thanks for getting back :) Yeah, the interfacial current density is off by a factor of 10. But what I find strange is that if I just add a factor of 10 (why not?) so that the BV equation becomes
10 * j0 * (pybamm.exp(alpha * F / (R * T) * eta) - pybamm.exp(- (1 - alpha) * F / (R * T) * eta))
the results seem to be more or less unchanged. Even when I make it:
10**6 * j0 * (pybamm.exp(alpha * F / (R * T) * eta) - pybamm.exp(- (1 - alpha) * F / (R * T) * eta))
There's hardly any difference in the results. That's what's making me think there's an error in the way I've included the BV in my governing system of equations.
Unfortunately I'm not allowed to share the LIBAT data, sorry!

[valentinsulzer]

The interfacial current density is governed by the Ohm's law equation (and boundary condition) for the potentials, and then the Butler-Volmer equation determines what the overpotential is given a fixed interfacial current density. Of course it's a bit more coupled than that but that's the general idea. So the error is in Ohm's law or its boundary condition. For example why is the source terms a_j * F not a_j / F? It can't be just that because then you'd be off by F**2 (which is not 10) but that's the kind of thing I'd check rather than BV

[valentinsulzer]

In fact, a / F = 3 * eps / R / F is roughly 10 for typical values of eps and R. So you could be off by that factor somewhere

[isaacbasil]

Ok cheers, I will go back through and have another look for factors like that which might be missing in my governing eqs and BCs.

[isaacbasil]

The error was indeed in the BC on Ohm's law, thanks for that pointer! I was using an intrinsic conductivity rather than an effective conductivity, which threw things off by a factor of ~8. Works perfectly now, cheers.

[valentinsulzer]

Also, you're welcome to model the battery however you want, but homogenizing the particle equations so that they appear as macroscale equations is pretty questionable, as explained in this paper https://iopscience.iop.org/article/10.1088/2516-1083/ac7d31/pdf

However, the homogenised battery model (unlike the porous medium flow equation, see [136]) retains a
microscale variable in order to model lithium transport within the electrode particles. This is caused by two
main reasons. First, lithium-ion motion within the active materials that form the particles is extremely slow,
so that significant lithium-ion concentration gradients form on the microscale in order to drive the requisite
fluxes. Secondly, lithium-ions cannot flow directly between the discrete microscopic electrode particles.
Therefore, when the microscale model is homogenised the particle equations do not upscale and only the
electrolyte equations (6) and the current flow equations in the matrix (1b) and (1d) do (see [52]). These
macroscopic equations couple to a series of microscale problems, parameterised by the macroscopic variable
X, describing lithium-ion transport within the particles at representative points in the electrodes.

[isaacbasil]

Yeah that's a fair point! I'm not expecting it to be a good model, it's mainly a stepping stone so that I can get my own model working in Pybamm before I start increasing the complexity.