Defining User Defined Functions Taking Electrode SOC as Input

[MirAbbasAli2A]

Hello -

In an attempt to add the effects of electrode SOC on the solid phase diffusivity, I wanted to input the data from a curve plotting electrode SOC vs. electrode diffusivity. I wish to call on this data by using the electrode SOC as an input and use this to adjust the current electrode diffusivity definition that already contains the Arrhenius relationship (i.e. Electrode diffusivity = Arrhenius correction x SOC correction). I am aware that the O'Regan dataset already has provisions for this but it doesn't have any of the aging parameters I needed for my work. I was wondering if there is a way to perform the requested function with PyBamm using any alternative dataset (I am using OKane2022 for reference).

Thank you for your time.

[MirAbbasAli2A]

Hello,

Does anyone from the team have any suggestions regarding this request?

[brosaplanella]

If I understand correctly, you want to define a diffusion coefficient that depends on SoC? You can add this as a parameter, for example you can see how it is done in the ORegan2022 dataset and follow that route. If your values come from data, you can create an interpolant from data and pass that as your parameter.

[MirAbbasAli2A]

@brosaplanella , thank you for your reply. I understand how to create the interpolant from the data. The part I am confused about however is how to ensure that the electrode SOC from pybamm is used as the input? In the pybamm source code I have seen it referenced as "sto". Would it be something as simple as :

diff_interp = pybamm.Interpolant(data[:,0], data[:,-1], sto)? #where column one in data is electrode SOC and the last column is the diffusivity

[MirAbbasAli2A]

@brosaplanella, an update:

I implemented the same function in the ORegan2022 parameter set as below:

T = model.variables["Volume-averaged cell temperature [K]"]
sto = model.variables["Negative electrode stoichiometry"]
sto_pos = model.variables["Positive electrode stoichiometry"]

def graphite_LGM50_diffusivity_ORegan2022(sto, T):

    a0 = 11.17
    a1 = -1.553
    a2 = -6.136
    a3 = -9.725
    a4 = 1.85
    b1 = 0.2031
    b2 = 0.5375
    b3 = 0.9144
    b4 = 0.5953
    c0 = -15.11
    c1 = 0.0006091
    c2 = 0.06438
    c3 = 0.0578
    c4 = 0.001356
    d = 2092

    D_ref = (
        10
        ** (
            a0 * sto
            + c0
            + a1 * exp(-((sto - b1) ** 2) / c1)
            + a2 * exp(-((sto - b2) ** 2) / c2)
            + a3 * exp(-((sto - b3) ** 2) / c3)
            + a4 * exp(-((sto - b4) ** 2) / c4)
        )
        * 3.0321  # correcting factor (see O'Regan et al 2021)
    )

    E_D_s = d * pybamm.constants.R
    arrhenius_v3 = exp(E_D_s / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return D_ref * arrhenius_v3

def nmc_LGM50_diffusivity_ORegan2022(sto_pos, T):


    a1 = -0.9231
    a2 = -0.4066
    a3 = -0.993
    b1 = 0.3216
    b2 = 0.4532
    b3 = 0.8098
    c0 = -13.96
    c1 = 0.002534
    c2 = 0.003926
    c3 = 0.09924
    d = 1449

    D_ref_v2 = (
        10
        ** (
            c0
            + a1 * exp(-((sto_pos - b1) ** 2) / c1)
            + a2 * exp(-((sto_pos - b2) ** 2) / c2)
            + a3 * exp(-((sto_pos - b3) ** 2) / c3)
        )
        * 2.7  # correcting factor (see O'Regan et al 2021)
    )

    E_D_s_v2 = d * pybamm.constants.R
    arrhenius_v4 = exp(E_D_s_v2 / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return D_ref_v2 * arrhenius_v4
    
param.update({"Negative electrode diffusivity [m2.s-1]": graphite_LGM50_diffusivity_ORegan2022(sto, T)})
param.update({"Positive electrode diffusivity [m2.s-1]": nmc_LGM50_diffusivity_ORegan2022(sto_pos, T)})

However, I get a wide array of key/type errors as follows:

KeyError: Multiplication(0x6df559e3795d0557, *, children=['broadcast(-1.0 / ((Negative particle radius [m] / Negative particle radius [m]) ** 2.0))', "-div(Negative electrode diffusivity [m2.s-1] * (1.0 + ((Negative electrode partial molar volume [m3.mol-1] / (Ideal gas constant [J.K-1.mol-1] * Volume-averaged cell temperature [K])) * 2.0 * Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / (9.0 - (9.0 * Negative electrode Poisson's ratio))) * (Negative particle concentration [mol.m-3] - Negative electrode reference concentration for free of deformation [mol.m-3])) * grad(Negative particle concentration [mol.m-3]))"], domains={'primary': ['negative particle'], 'secondary': ['negative electrode'], 'tertiary': ['current collector']})


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:601 in process_symbol
    return self._processed_symbols[symbol]

KeyError: Negate(-0x5744fb8bb4986d27, -, children=["div(Negative electrode diffusivity [m2.s-1] * (1.0 + ((Negative electrode partial molar volume [m3.mol-1] / (Ideal gas constant [J.K-1.mol-1] * Volume-averaged cell temperature [K])) * 2.0 * Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / (9.0 - (9.0 * Negative electrode Poisson's ratio))) * (Negative particle concentration [mol.m-3] - Negative electrode reference concentration for free of deformation [mol.m-3])) * grad(Negative particle concentration [mol.m-3]))"], domains={'primary': ['negative particle'], 'secondary': ['negative electrode'], 'tertiary': ['current collector']})


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:601 in process_symbol
    return self._processed_symbols[symbol]

KeyError: Divergence(-0x4308714659fd6741, div, children=["Negative electrode diffusivity [m2.s-1] * (1.0 + ((Negative electrode partial molar volume [m3.mol-1] / (Ideal gas constant [J.K-1.mol-1] * Volume-averaged cell temperature [K])) * 2.0 * Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / (9.0 - (9.0 * Negative electrode Poisson's ratio))) * (Negative particle concentration [mol.m-3] - Negative electrode reference concentration for free of deformation [mol.m-3])) * grad(Negative particle concentration [mol.m-3])"], domains={'primary': ['negative particle'], 'secondary': ['negative electrode'], 'tertiary': ['current collector']})


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:601 in process_symbol
    return self._processed_symbols[symbol]

KeyError: Multiplication(-0x6d08da98c5656f45, *, children=["Negative electrode diffusivity [m2.s-1] * (1.0 + ((Negative electrode partial molar volume [m3.mol-1] / (Ideal gas constant [J.K-1.mol-1] * Volume-averaged cell temperature [K])) * 2.0 * Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / (9.0 - (9.0 * Negative electrode Poisson's ratio))) * (Negative particle concentration [mol.m-3] - Negative electrode reference concentration for free of deformation [mol.m-3]))", 'grad(Negative particle concentration [mol.m-3])'], domains={'primary': ['negative particle'], 'secondary': ['negative electrode'], 'tertiary': ['current collector']})


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:601 in process_symbol
    return self._processed_symbols[symbol]

KeyError: Multiplication(-0x38493770cb765e6e, *, children=['Negative electrode diffusivity [m2.s-1]', "1.0 + ((Negative electrode partial molar volume [m3.mol-1] / (Ideal gas constant [J.K-1.mol-1] * Volume-averaged cell temperature [K])) * 2.0 * Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / (9.0 - (9.0 * Negative electrode Poisson's ratio))) * (Negative particle concentration [mol.m-3] - Negative electrode reference concentration for free of deformation [mol.m-3])"], domains={'primary': ['negative particle'], 'secondary': ['negative electrode'], 'tertiary': ['current collector']})


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:601 in process_symbol
    return self._processed_symbols[symbol]

KeyError: FunctionParameter(-0x50d8baccdbd7eb18, Negative electrode diffusivity [m2.s-1], children=['maximum(minimum(Negative particle concentration [mol.m-3] / Maximum concentration in negative electrode [mol.m-3], 0.9999999999), 1e-10)', 'broadcast(broadcast(broadcast(Volume-averaged cell temperature [K])))'], domains={'primary': ['negative particle'], 'secondary': ['negative electrode'], 'tertiary': ['current collector']})


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File ~\anaconda3\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\users\abbas\documents\phd docs\pybamm\agingmodel\brosaplanella-tec-reduced-model-978f8c3\agingmodelwithmechanismsanddynamiccycle.py:233
    sol = sim.solve(save_at_cycles=10)

  File ~\anaconda3\lib\site-packages\pybamm\simulation.py:647 in solve
    self.build_for_experiment(check_model=check_model, initial_soc=initial_soc)

  File ~\anaconda3\lib\site-packages\pybamm\simulation.py:482 in build_for_experiment
    self.set_up_and_parameterise_experiment()

  File ~\anaconda3\lib\site-packages\pybamm\simulation.py:216 in set_up_and_parameterise_experiment
    self.set_up_and_parameterise_model_for_experiment()

  File ~\anaconda3\lib\site-packages\pybamm\simulation.py:289 in set_up_and_parameterise_model_for_experiment
    parameterised_model = new_parameter_values.process_model(

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:436 in process_model
    new_rhs[new_variable] = self.process_symbol(equation)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:603 in process_symbol
    processed_symbol = self._process_symbol(symbol)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:715 in _process_symbol
    new_right = self.process_symbol(symbol.right)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:603 in process_symbol
    processed_symbol = self._process_symbol(symbol)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:723 in _process_symbol
    new_child = self.process_symbol(symbol.child)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:603 in process_symbol
    processed_symbol = self._process_symbol(symbol)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:723 in _process_symbol
    new_child = self.process_symbol(symbol.child)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:603 in process_symbol
    processed_symbol = self._process_symbol(symbol)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:714 in _process_symbol
    new_left = self.process_symbol(symbol.left)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:603 in process_symbol
    processed_symbol = self._process_symbol(symbol)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:714 in _process_symbol
    new_left = self.process_symbol(symbol.left)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:603 in process_symbol
    processed_symbol = self._process_symbol(symbol)

  File ~\anaconda3\lib\site-packages\pybamm\parameters\parameter_values.py:697 in _process_symbol
    raise TypeError(

TypeError: Parameter provided for 'Negative electrode diffusivity [m2.s-1]' is of the wrong type (should either be scalar-like or callable).

FYI, I am using the OKane2022 dataset and updating the diffusivity function with the ORegan2022 function. Why is this happening and how do I fix this?

Thank you for your continued support.

[brosaplanella]

You can do it similarly to the OCPs in Ai2020.py:

def lico2_ocp_Ai2020(sto):
    name, (x, y) = lico2_ocp_Ai2020_data
    return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")

where lico2_ocp_Ai2020_data contains the data.

Regarding your second post, you need to pass the parameters as

param.update({"Negative electrode diffusivity [m2.s-1]": graphite_LGM50_diffusivity_ORegan2022})
param.update({"Positive electrode diffusivity [m2.s-1]": nmc_LGM50_diffusivity_ORegan2022})

[MirAbbasAli2A]

Thank you @brosaplanella , this worked.

If you don't mind me asking you just 1 more question on this topic. Adding an electrode stoichiometry based diffusivity really slows down the simulation time. Do you have any suggestions to help increase the sim speed while still incorporating his feature? I kept my mesh with the following points since I am modelling all the coupled degradation mechanisms:

var_pts = {
"x_n" : 20, # x-direction, length, negative electrode
"x_s" : 20, # x-direction, length, separator
"x_p" : 20, # x-direction, length, positive electrode
"r_n" : 30, # number of volumes in the radial direction, negative particle
"r_p" : 30 # number of volumes in the radial direction, positive particle
}

model = pybamm.lithium_ion.DFN(options = {
"thermal":"lumped",
"dimensionality":0,
"cell geometry": "arbitrary",
"SEI": "solvent-diffusion limited",
"SEI porosity change": "true",
"particle mechanics": ("swelling and cracking", "swelling only"),
"SEI on cracks": "true",
"loss of active material": "stress-driven",
"calculate discharge energy": "true",
"lithium plating": "partially reversible",
"lithium plating porosity change": "true"
},
name="TDFN",
)

Any advice would be greatly appreciated. Thank you so much for your continued support.

[brosaplanella]

If your measurements are noisy, using an interpolant makes it really slow as the solver struggles to find the right time step. In those cases, it is better to fit a smooth function through the data and pass that as the parameter.

[MirAbbasAli2A]

@brosaplanella , I didn't actually end up using the interpolant method. I used the functions for diffusivity implemented in the ORegan2022 dataset which i believe fit a smooth function through the data.

[brosaplanella]

See this notebook on some tricks to improve the convergence for this parameter set.