[Bug]: hysteresis state error

[llh929423719]

PyBaMM Version

24.5

Python Version

3.9

Describe the bug

Under some specific parameter sets, hysteresis state will exceed the range of (-1,1)，adding bounds in WyciskOpenCircuitPotential submodel seems to have no effect.

Steps to Reproduce

my hysteresis paramset:
'Secondary: Negative particle hysteresis decay rate': 0.005,
'Secondary: Negative particle hysteresis switching factor': 0.2,

Relevant log output

No response

[rtimms]

@mleot

[mleot]

@llh929423719 can you provide more details on the specific base parameter set, the experimental protocol, the model used, and the solver? I would like to try and replicate the error you are experiencing.

[llh929423719]

paramset:

import pybamm
import os
import numpy as np

y_low = 0.18
y_up = 0.819
y_mid = (y_low + y_up) / 2

x_gr_low = 0.053410
x_gr_up = 0.86455
x_gr_mid = (x_gr_low + x_gr_up) / 2
x_si_low = 0.18072
x_si_up = 0.82146
x_si_mid = (x_si_low + x_si_up) / 2


def si_lith_yk(sto):
    si_lith_sto = np.array([0.000446504811095948, 0.00563431441023561, 0.0181660574920546,
                            0.0334360461827143, 0.0460270060067983, 0.0613297753940689,
                            0.0784344256531157, 0.108160701549205, 0.135206891136893,
                            0.185718772288891, 0.250676366228057, 0.292182015346653,
                            0.340003235771988, 0.378801305094702, 0.404968231803938,
                            0.445567124556602, 0.485262433268812, 0.525858153695997,
                            0.565552404966381, 0.606152355160871, 0.651267053232150,
                            0.698186804500684, 0.744206144054243, 0.791126952764602,
                            0.834437655080453, 0.852483957285695, ])

    si_lith_ocp = np.array([0.994727595054768, 0.742820630901517, 0.627991965557005,
                            0.547139434807918, 0.497922919892966, 0.453391400988814,
                            0.405343888012681, 0.342057580879760, 0.309235644037588,
                            0.276399960451676, 0.249414045048080, 0.237673268452228,
                            0.223585499723215, 0.211846309290103, 0.204801103123313,
                            0.188374273074835, 0.170776326203885, 0.150834559525251,
                            0.132064967110916, 0.116809782605822, 0.103895245582934,
                            0.0909796511182185, 0.0804078764611865, 0.0686639275398565,
                            0.0569220935021789, 0.0522249369103773,
                            ])

    return pybamm.Interpolant(x=si_lith_sto, y=si_lith_ocp,
                              children=sto, interpolator='linear', )


# si_lith_yk.__name__ = 'silicon_ocp_lithiation_yk'
def si_delith_yk(sto, ):
    si_delith_sto = np.array([0.000443332485617468, 0.00853329117665252, 0.0220278350411900,
                              0.0364249055043548, 0.0535263834379232, 0.0724413741033554,
                              0.0940650020064959, 0.113889921363339, 0.138224301387840,
                              0.165268376091875, 0.194118561434991, 0.225679498899467,
                              0.256341082090795, 0.293325638681645, 0.328506199517066,
                              0.368198335903798, 0.406084361651449, 0.440361338446417,
                              0.474634085474080, 0.509808301658545, 0.557614717898313,
                              0.600915903237728, 0.640612269391765, 0.684821268538938,
                              0.726322687890229, 0.765114412561986, 0.816524590545307,
                              0.839051802488901, 0.849834536790252, 0.853386483884322,
                              ])

    si_delith_ocp = np.array([0.991212658424613, 0.954886888091460, 0.906841489998980, 0.858795563185586,
                              0.807233113579298, 0.765042770878264, 0.724022487557875, 0.690033134939623,
                              0.652526202086650, 0.617360974157708, 0.583366334330325, 0.552885044970358,
                              0.525919220961458, 0.504807923623092, 0.484869329269937, 0.463756445768832,
                              0.441472974166168, 0.420363262990541, 0.394566969641374, 0.367598502027909,
                              0.337107695691506, 0.314821051763364, 0.298394750435798, 0.281965805503668,
                              0.265538446734276, 0.246769383040854, 0.209246588560489, 0.169397422062572,
                              0.116667027959291, 0.0522244081894643,
                              ])

    return pybamm.Interpolant(x=si_delith_sto, y=si_delith_ocp,
                              children=sto, interpolator='linear', )


# si_delith_yk.__name__ = 'silicon_ocp_delithiation_yk'

def nca_ls_ocp(sto):
    """
    Parameters
    ----------
    sto: :class:`pybamm.Symbol`
        Electrode stochiometry

    Returns
    -------
    :class:`pybamm.Symbol`
        Open-circuit potential
    """

    u_eq = (
            4.5
            - 2.02303 * sto
            + 8.18849 * (sto ** 2)
            - 23.6282 * (sto ** 3)
            + 27.8839 * (sto ** 4)
            - 11.4928 * (sto ** 5)
            - 0.5623e-4 * pybamm.exp(109.451 * sto - 100.006)
    )
    # update 240311
    return u_eq

def nca_ls_ocp_new(sto):
    """
    Parameters
    ----------
    sto: :class:`pybamm.Symbol`
        Electrode stochiometry

    Returns
    -------
    :class:`pybamm.Symbol`
        Open-circuit potential
    """

    u_eq = (
            4.5
            - 2.9883 * sto
            + 11.6890 * (sto ** 2)
            - 29.6247 * (sto ** 3)
            + 32.7684 * (sto ** 4)
            - 13.0161 * (sto ** 5)
            - 0.5623e-4 * pybamm.exp(109.451 * sto - 100.006)
    )
    # update 240311
    return u_eq


def graphite_ls_electrolyte_exchange_current_density(
        c_e, c_s_surf, c_s_max, T
):
    """

    Parameters
    ----------
    c_e : :class:`pybamm.Symbol`
        Electrolyte concentration [mol.m-3]
    c_s_surf : :class:`pybamm.Symbol`
        Particle concentration [mol.m-3]
    c_s_max : :class:`pybamm.Symbol`
        Maximum particle concentration [mol.m-3]
    T : :class:`pybamm.Symbol`
        Temperature [K]

    Returns
    -------
    :class:`pybamm.Symbol`
        Exchange-current density [A.m-2]
    """
    m_ref = 2.65e-7*pybamm.constants.F/1000 ** 0.5  # (A/m2)(m3/mol)**1.5 - includes ref concentrations
    E_r = pybamm.Parameter("Primary: Negative electrode exchange-current density activation energy")
    arrhenius = np.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return m_ref * arrhenius * c_e ** 0.5 * c_s_surf ** 0.5 * (c_s_max - c_s_surf) ** 0.5


def silicon_ls_electrolyte_exchange_current_density(
        c_e, c_s_surf, c_s_max, T
):
    """

    Parameters
    ----------
    c_e : :class:`pybamm.Symbol`
        Electrolyte concentration [mol.m-3]
    c_s_surf : :class:`pybamm.Symbol`
        Particle concentration [mol.m-3]
    c_s_max : :class:`pybamm.Symbol`
        Maximum particle concentration [mol.m-3]
    T : :class:`pybamm.Symbol`
        Temperature [K]

    Returns
    -------
    :class:`pybamm.Symbol`
        Exchange-current density [A.m-2]
    """
    m_ref = 2.65e-7*pybamm.constants.F/1000 ** 0.5  # (A/m2)(m3/mol)**1.5 - includes ref concentrations
    E_r = pybamm.Parameter("Secondary: Negative electrode exchange-current density activation energy")
    arrhenius = np.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return m_ref * arrhenius * c_e ** 0.5 * c_s_surf ** 0.5 * (c_s_max - c_s_surf) ** 0.5


def ds_gr_ls(sto, T):
    D_ref = 1.12e-13
    E_D_s = pybamm.Parameter("Primary: Negative electrode diffusivity activation energy")
    arrhenius = np.exp(E_D_s / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return D_ref * arrhenius * 0.3


def ds_si_ls(sto, T):
    D_ref = 4.0e-15
    E_D_s = pybamm.Parameter("Secondary: Negative electrode diffusivity activation energy")
    arrhenius = np.exp(E_D_s / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return D_ref * arrhenius * 0.1


def ds_nac_ls(sto, T):
    D_ref = 20.0e-16
    E_D_s = pybamm.Parameter("Positive electrode diffusivity activation energy")
    arrhenius = np.exp(E_D_s / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return D_ref * arrhenius * 30 * 0.5


def silicon_ocp_lithiation_Mark2016(sto):
    """
    silicon Open-circuit Potential (OCP) as a a function of the
    stochiometry. The fit is taken from the Enertech cell [1], which is only accurate
    for 0 < sto < 1.

    References
    ----------
    .. [1] Verbrugge M, Baker D, Xiao X. Formulation for the treatment of multiple
    electrochemical reactions and associated speciation for the Lithium-Silicon
    electrode[J]. Journal of The Electrochemical Society, 2015, 163(2): A262.

    Parameters
    ----------
    sto: double
       Stochiometry of material (li-fraction)

    Returns
    -------
    :class:`pybamm.Symbol`
        OCP [V]
    """
    p1 = -96.63
    p2 = 372.6
    p3 = -587.6
    p4 = 489.9
    p5 = -232.8
    p6 = 62.99
    p7 = -9.286
    p8 = 0.8633

    U_lithiation = (
            p1 * sto ** 7
            + p2 * sto ** 6
            + p3 * sto ** 5
            + p4 * sto ** 4
            + p5 * sto ** 3
            + p6 * sto ** 2
            + p7 * sto
            + p8
    )
    return U_lithiation - 0.02


def silicon_ocp_delithiation_Mark2016(sto):
    """
    silicon Open-circuit Potential (OCP) as a a function of the
    stochiometry. The fit is taken from the Enertech cell [1], which is only accurate
    for 0 < sto < 1.

    References
    ----------
    .. [1] Verbrugge M, Baker D, Xiao X. Formulation for the treatment of multiple
    electrochemical reactions and associated speciation for the Lithium-Silicon
    electrode[J]. Journal of The Electrochemical Society, 2015, 163(2): A262.

    Parameters
    ----------
    sto: double
       Stochiometry of material (li-fraction)

    Returns
    -------
    :class:`pybamm.Symbol`
        OCP [V]
    """
    p1 = -51.02
    p2 = 161.3
    p3 = -205.7
    p4 = 140.2
    p5 = -58.76
    p6 = 16.87
    p7 = -3.792
    p8 = 0.9937

    U_delithiation = (
            p1 * sto ** 7
            + p2 * sto ** 6
            + p3 * sto ** 5
            + p4 * sto ** 4
            + p5 * sto ** 3
            + p6 * sto ** 2
            + p7 * sto
            + p8
    )
    return U_delithiation - 0.02


def silicon_ocp_avg(sto):
    p1 = (-96.63 - 51.02) / 2
    p2 = (372.6 + 161.3) / 2
    p3 = (-587.6 - 205.7) / 2
    p4 = (489.9 + 140.2) / 2
    p5 = (-232.8 - 58.76) / 2
    p6 = (62.99 + 16.87) / 2
    p7 = (-9.286 - 3.792) / 2
    p8 = (0.8633 + 0.9937) / 2

    U_avg = (
            p1 * sto ** 7
            + p2 * sto ** 6
            + p3 * sto ** 5
            + p4 * sto ** 4
            + p5 * sto ** 3
            + p6 * sto ** 2
            + p7 * sto
            + p8
    )
    return U_avg - 0.02


def silicon_ocp_delithiation_llh_fit(sto):
    p1 = -49.4254
    p2 = 154.9471
    p3 = -195.3232
    p4 = 131.4335
    p5 = -54.7137
    p6 = 15.8751
    p7 = -3.6803
    p8 = 0.9910

    U_delithiation = (
            p1 * sto ** 7
            + p2 * sto ** 6
            + p3 * sto ** 5
            + p4 * sto ** 4
            + p5 * sto ** 3
            + p6 * sto ** 2
            + p7 * sto
            + p8
    )
    return U_delithiation


def silicon_ocp_lithiation_llh_fit(sto):
    p1 = -14.8592
    p2 = 63.6160
    p3 = -118.6097
    p4 = 124.4181
    p5 = -77.8073
    p6 = 28.1588
    p7 = -5.5981
    p8 = 0.7317

    U_lithiation = (
            p1 * sto ** 7
            + p2 * sto ** 6
            + p3 * sto ** 5
            + p4 * sto ** 4
            + p5 * sto ** 3
            + p6 * sto ** 2
            + p7 * sto
            + p8
    )
    return U_lithiation


def silicon_ocp_delithiation_llh_origin_fit(sto):
    p1 = -239.6098
    p2 = 650.1135
    p3 = -701.3641
    p4 = 391.8107
    p5 = -127.1997
    p6 = 26.8382
    p7 = -4.5558
    p8 = 0.9934

    U_delithiation = (
            p1 * sto ** 7
            + p2 * sto ** 6
            + p3 * sto ** 5
            + p4 * sto ** 4
            + p5 * sto ** 3
            + p6 * sto ** 2
            + p7 * sto
            + p8
    )
    return U_delithiation


def silicon_ocp_lithiation_llh_origin_fit(sto):
    p1 = -44.8392
    p2 = 163.9680
    p3 = -261.1214
    p4 = 233.9558
    p5 = -124.9688
    p6 = 38.6309
    p7 = -6.5599
    p8 = 0.7322

    U_lithiation = (
            p1 * sto ** 7
            + p2 * sto ** 6
            + p3 * sto ** 5
            + p4 * sto ** 4
            + p5 * sto ** 3
            + p6 * sto ** 2
            + p7 * sto
            + p8
    )
    return U_lithiation

def silicon_ocp_Eeq_llh_origin_fit(sto):

    return (silicon_ocp_lithiation_llh_origin_fit(sto)+silicon_ocp_delithiation_llh_origin_fit(sto))/2

def silicon_LGM50_electrolyte_exchange_current_density_Chen2020(
        c_e, c_s_surf, c_s_max, T
):
    """
    Exchange-current density for Butler-Volmer reactions between silicon and LiPF6 in
    EC:DMC.

    References
    ----------
    .. [1] Chang-Hui Chen, Ferran Brosa Planella, Kieran O’Regan, Dominika Gastol, W.
    Dhammika Widanage, and Emma Kendrick. "Development of Experimental Techniques for
    Parameterization of Multi-scale Lithium-ion Battery Models." Journal of the
    Electrochemical Society 167 (2020): 080534.

    Parameters
    ----------
    c_e : :class:`pybamm.Symbol`
        Electrolyte concentration [mol.m-3]
    c_s_surf : :class:`pybamm.Symbol`
        Particle concentration [mol.m-3]
    c_s_max : :class:`pybamm.Symbol`
        Maximum particle concentration [mol.m-3]
    T : :class:`pybamm.Symbol`
        Temperature [K]

    Returns
    -------
    :class:`pybamm.Symbol`
        Exchange-current density [A.m-2]
    """

    m_ref = (
            6.48e-7 * 28700 / 278000
    )  # (A/m2)(m3/mol)**1.5 - includes ref concentrations
    E_r = 35000
    arrhenius = np.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return m_ref * arrhenius * c_e ** 0.5 * c_s_surf ** 0.5 * (c_s_max - c_s_surf) ** 0.5


def nmc_LGM50_ocp_Chen2020(sto):
    """
    LG M50 NMC open-circuit potential as a function of stochiometry, fit taken
    from [1].

    References
    ----------
    .. [1] Chang-Hui Chen, Ferran Brosa Planella, Kieran O’Regan, Dominika Gastol, W.
    Dhammika Widanage, and Emma Kendrick. "Development of Experimental Techniques for
    Parameterization of Multi-scale Lithium-ion Battery Models." Journal of the
    Electrochemical Society 167 (2020): 080534.

    Parameters
    ----------
    sto: :class:`pybamm.Symbol`
        Electrode stochiometry

    Returns
    -------
    :class:`pybamm.Symbol`
        Open-circuit potential
    """

    u_eq = (
            -0.8090 * sto
            + 4.4875
            - 0.0428 * np.tanh(18.5138 * (sto - 0.5542))
            - 17.7326 * np.tanh(15.7890 * (sto - 0.3117))
            + 17.5842 * np.tanh(15.9308 * (sto - 0.3120))
    )

    return u_eq


def nca_j0(c_e, c_s_surf, c_s_max, T):
    """
    Exchange-current density for Butler-Volmer reactions between NMC and LiPF6 in
    EC:DMC.

    References
    ----------
    .. [1] Chang-Hui Chen, Ferran Brosa Planella, Kieran O’Regan, Dominika Gastol, W.
    Dhammika Widanage, and Emma Kendrick. "Development of Experimental Techniques for
    Parameterization of Multi-scale Lithium-ion Battery Models." Journal of the
    Electrochemical Society 167 (2020): 080534.

    Parameters
    ----------
    c_e : :class:`pybamm.Symbol`
        Electrolyte concentration [mol.m-3]
    c_s_surf : :class:`pybamm.Symbol`
        Particle concentration [mol.m-3]
    c_s_max : :class:`pybamm.Symbol`
        Maximum particle concentration [mol.m-3]
    T : :class:`pybamm.Symbol`
        Temperature [K]

    Returns
    -------
    :class:`pybamm.Symbol`
        Exchange-current density [A.m-2]
    """
    m_ref = 1.0e-7 * pybamm.constants.F / 1000 ** 0.5  # (A/m2)(m3/mol)**1.5 - includes ref concentrations
    E_r = pybamm.Parameter("Positive electrode exchange-current density activation energy")
    arrhenius = np.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return m_ref * arrhenius * c_e ** 0.5 * c_s_surf ** 0.5 * (c_s_max - c_s_surf) ** 0.5


def electrolyte_diffusivity_ls_46(c_e, T):
    """
    Parameters
    ----------
    c_e: :class:`pybamm.Symbol`
        Dimensional electrolyte concentration
    T: :class:`pybamm.Symbol`
        Dimensional temperature

    Returns
    -------
    :class:`pybamm.Symbol`
        Solid diffusivity
    """
    p1 = 4.074e6
    p2 = -2.386
    p3 = -4405.6
    p4 = 680.36
    D_c_e = p1 * pybamm.exp(p2 * c_e / 1000) * pybamm.exp(p3 / T) * pybamm.exp(
        p4 * c_e / 1000 / T) * 1e-10  # unit in m2/s

    return D_c_e * 3


def electrolyte_conductivity_ls_46(c_e, T):
    """
    Conductivity of LiPF6 in EC:EMC (3:7) as a function of ion concentration. The data
    comes from [1], with Arrhenius temperature dependence added from [2].

    Parameters
    ----------
    c_e: :class:`pybamm.Symbol`
        Dimensional electrolyte concentration
    T: :class:`pybamm.Symbol`
        Dimensional temperature

    Returns
    -------
    :class:`pybamm.Symbol`
        Solid diffusivity
    """
    p1 = 1.581
    p2 = 233.6
    p3 = -1.84
    p4 = 0.986
    p5 = 1.2e-3
    p6 = 0.022
    c_e_m = c_e / 1000
    sigma_e = p1 * (1 + (T - p2)) * c_e_m * (
            1 + p3 * pybamm.sqrt(c_e_m) + p4 * (1 + p5 * pybamm.exp(1000 / T)) * c_e_m) / (
                      1 + c_e_m ** 4 * (p6 * pybamm.exp(1000 / T))) / 10

    return sigma_e


def electrolyte_t0plus(c_e, T):
    p1 = -38.727
    p2 = 12.940
    p3 = 0.25
    p4 = -0.292
    p5 = -0.083
    p6 = -4e-4
    p7 = -0.027
    p8 = 0.0013
    p9 = 1.29e-4
    c_e_m = c_e / 1000
    t0plus = p1 + p2 * c_e_m + p3 * T + p4 * c_e_m ** 2 + p5 * c_e_m * T + p6 * T ** 2 + p7 * c_e_m ** 3 + p8 * c_e_m ** 2 * T + p9 * c_e_m * T ** 2
    return t0plus


def electrolyte_diffusivity_Nyman2008(c_e, T):
    """
    Diffusivity of LiPF6 in EC:EMC (3:7) as a function of ion concentration. The data
    comes from [1]

    References
    ----------
    .. [1] A. Nyman, M. Behm, and G. Lindbergh, "Electrochemical characterisation and
    modelling of the mass transport phenomena in LiPF6-EC-EMC electrolyte,"
    Electrochim. Acta, vol. 53, no. 22, pp. 6356–6365, 2008.

    Parameters
    ----------
    c_e: :class:`pybamm.Symbol`
        Dimensional electrolyte concentration
    T: :class:`pybamm.Symbol`
        Dimensional temperature

    Returns
    -------
    :class:`pybamm.Symbol`
        Solid diffusivity
    """

    D_c_e = 8.794e-11 * (c_e / 1000) ** 2 - 3.972e-10 * (c_e / 1000) + 4.862e-10

    # Nyman et al. (2008) does not provide temperature dependence

    return D_c_e


def electrolyte_conductivity_Nyman2008(c_e, T):
    """
    Conductivity of LiPF6 in EC:EMC (3:7) as a function of ion concentration. The data
    comes from [1].

    References
    ----------
    .. [1] A. Nyman, M. Behm, and G. Lindbergh, "Electrochemical characterisation and
    modelling of the mass transport phenomena in LiPF6-EC-EMC electrolyte,"
    Electrochim. Acta, vol. 53, no. 22, pp. 6356–6365, 2008.

    Parameters
    ----------
    c_e: :class:`pybamm.Symbol`
        Dimensional electrolyte concentration
    T: :class:`pybamm.Symbol`
        Dimensional temperature

    Returns
    -------
    :class:`pybamm.Symbol`
        Solid diffusivity
    """

    sigma_e = (
            0.1297 * (c_e / 1000) ** 3 - 2.51 * (c_e / 1000) ** 1.5 + 3.329 * (c_e / 1000)
    )

    # Nyman et al. (2008) does not provide temperature dependence

    return sigma_e


# # Load data in the appropriate format
# path, _ = os.path.split(os.path.abspath(__file__))
# graphite_ocp_Enertech_Ai2020_data = pybamm.parameters.process_1D_data(
#     "graphite_ocp_Enertech_Ai2020.csv", path=path
# )
#
#
# def graphite_ocp_Enertech_Ai2020(sto):
#     name, (x, y) = graphite_ocp_Enertech_Ai2020_data
#     return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")


def graphite_ocp_ls(sto):
    ocp = (
            0.0646988 + 0.8 * pybamm.exp(-27.7807 * (sto - 0.0327108))
            - 0.012 * pybamm.tanh((sto + 30.7083) / 1.96366)
            - 0.0118 * pybamm.tanh((sto + 31.1913) / 3.7934)
            - 0.0035 * pybamm.tanh((sto - 0.22) / 0.02)
            - 0.0095 * pybamm.tanh((sto - 0.19) / 0.013)
            - 0.0145 * pybamm.tanh((sto - 0.49) / 0.02)
            - 0.08 * pybamm.tanh((sto - 1.03) / 0.055)

    )
    return ocp


def dudt_mcmb(sto, c_s_max):
    du = (-0.0263 * sto ** 5 + 0.0693 * sto ** 4 - 0.0708 * sto ** 3 + 0.0349 * sto ** 2 - 0.0080 * sto + 0.0006)

    return du


# Call dict via a function to avoid errors when editing in place
def get_parameter_values():
    """
    Parameters for a composite graphite/silicon negative electrode, from the paper
    :footcite:t:`Ai2022`, based on the paper :footcite:t:`Chen2020`, and references
    therein.

    SEI parameters are example parameters for composite SEI on silicon/graphite. Both
    phases use the same values, from the paper :footcite:t:`Yang2017`
    """

    return {
        "chemistry": "lithium_ion",
        # sei
        "Primary: Ratio of lithium moles to SEI moles": 2.0,
        "Primary: Inner SEI reaction proportion": 0.5,
        "Primary: Inner SEI partial molar volume [m3.mol-1]": 9.585e-05,
        "Primary: Outer SEI partial molar volume [m3.mol-1]": 9.585e-05,
        "Primary: SEI reaction exchange current density [A.m-2]": 1.5e-07,
        "Primary: SEI resistivity [Ohm.m]": 200000.0,
        "Primary: Outer SEI solvent diffusivity [m2.s-1]": 2.5000000000000002e-22,
        "Primary: Bulk solvent concentration [mol.m-3]": 2636.0,
        "Primary: Inner SEI open-circuit potential [V]": 0.1,
        "Primary: Outer SEI open-circuit potential [V]": 0.8,
        "Primary: Inner SEI electron conductivity [S.m-1]": 8.95e-14,
        "Primary: Inner SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
        "Primary: Lithium interstitial reference concentration [mol.m-3]": 15.0,
        "Primary: Initial inner SEI thickness [m]": 2.5e-09,
        "Primary: Initial outer SEI thickness [m]": 2.5e-09,
        "Primary: EC initial concentration in electrolyte [mol.m-3]": 4541.0,
        "Primary: EC diffusivity [m2.s-1]": 2e-18,
        "Primary: SEI kinetic rate constant [m.s-1]": 1e-12,
        "Primary: SEI open-circuit potential [V]": 0.4,
        "Primary: SEI growth activation energy [J.mol-1]": 0.0,
        "Secondary: Ratio of lithium moles to SEI moles": 2.0,
        "Secondary: Inner SEI reaction proportion": 0.5,
        "Secondary: Inner SEI partial molar volume [m3.mol-1]": 9.585e-05,
        "Secondary: Outer SEI partial molar volume [m3.mol-1]": 9.585e-05,
        "Secondary: SEI reaction exchange current density [A.m-2]": 1.5e-07,
        "Secondary: SEI resistivity [Ohm.m]": 200000.0,
        "Secondary: Outer SEI solvent diffusivity [m2.s-1]": 2.5000000000000002e-22,
        "Secondary: Bulk solvent concentration [mol.m-3]": 2636.0,
        "Secondary: Inner SEI open-circuit potential [V]": 0.1,
        "Secondary: Outer SEI open-circuit potential [V]": 0.8,
        "Secondary: Inner SEI electron conductivity [S.m-1]": 8.95e-14,
        "Secondary: Inner SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
        "Secondary: Lithium interstitial reference concentration [mol.m-3]": 15.0,
        "Secondary: Initial inner SEI thickness [m]": 2.5e-09,
        "Secondary: Initial outer SEI thickness [m]": 2.5e-09,
        "Secondary: EC initial concentration in electrolyte [mol.m-3]": 4541.0,
        "Secondary: EC diffusivity [m2.s-1]": 2e-18,
        "Secondary: SEI kinetic rate constant [m.s-1]": 1e-12,
        "Secondary: SEI open-circuit potential [V]": 0.4,
        "Secondary: SEI growth activation energy [J.mol-1]": 0.0,
        'Secondary: Negative particle hysteresis decay rate': 0.005,
        'Secondary: Negative particle hysteresis switching factor': 0.2,
        "Positive electrode reaction-driven LAM factor [m3.mol-1]": 0.0,
        # cell
        "Negative current collector thickness [m]": 8e-06,  # modified by lou@240313
        "Negative electrode thickness [m]": 7.85e-05,  # modified by lou@240313
        "Separator thickness [m]": 1.3e-05,  # modified by lou@240313
        "Positive electrode thickness [m]": 7.2e-05,  # modified by lou@240313
        "Positive current collector thickness [m]": 1.2e-05,  # modified by lou@240313
        "Electrode height [m]": 0.074,  # modified by lou@240313
        "Electrode width [m]": 0.057,  # modified by lou@240313
        "Cell cooling surface area [m2]": 0.0085,  # modified by lou@240313
        "Cell volume [m3]": 3.84e-05,  # modified by lou@240325
        "Cell thermal expansion coefficient [m.K-1]": 1.1e-06,
        "Negative current collector conductivity [S.m-1]": 58411000.0,
        "Positive current collector conductivity [S.m-1]": 36914000.0,
        "Negative current collector density [kg.m-3]": 8960.0,
        "Positive current collector density [kg.m-3]": 2700.0,
        "Negative current collector specific heat capacity [J.kg-1.K-1]": 1126.0,  # modified@240429
        "Positive current collector specific heat capacity [J.kg-1.K-1]": 1126.0,  # modified@240429
        "Negative current collector thermal conductivity [W.m-1.K-1]": 401.0,
        "Positive current collector thermal conductivity [W.m-1.K-1]": 237.0,
        "Nominal cell capacity [A.h]": 2.6,
        "Current function [A]": 2.6,
        "Contact resistance [Ohm]": 0,
        # negative electrode
        "Negative electrode density [kg.m-3]": 2300.0,  # modified@240429
        "Negative electrode conductivity [S.m-1]": 26.8,  # modified@240513
        "Primary: Maximum concentration in negative electrode [mol.m-3]": 31915.0,  # modified@240429
        "Primary: Initial concentration in negative electrode [mol.m-3]": 31915.0 * 0.94,
        # modified@240429
        "Primary: Negative electrode diffusivity [m2.s-1]": ds_gr_ls,  # modified@240429
        "Primary: Negative electrode OCP [V]": graphite_ocp_ls,  # modified@240429
        "Negative electrode porosity": 0.365,  # modified@240429
        "Primary: Negative electrode active material volume fraction": 0.573,  # modified@240429
        "Primary: Negative particle radius [m]": 13.309e-06 / 2,  # modified@240429
        "Negative electrode Bruggeman coefficient (electrolyte)": 2.37,  # modified@240513
        "Negative electrode Bruggeman coefficient (electrode)": 0,
        "Negative electrode charge transfer coefficient": 0.5,  # modified@240429
        "Negative electrode double-layer capacity [F.m-2]": 0.2,
        "Primary: Negative electrode exchange-current density [A.m-2]"
        "": graphite_ls_electrolyte_exchange_current_density,  # modified@240429
        "Primary: Negative electrode density [kg.m-3]": 2300.0,  # modified@240429
        "Negative electrode specific heat capacity [J.kg-1.K-1]": 1126.0,  # modified@240429
        "Negative electrode thermal conductivity [W.m-1.K-1]": 1.7,
        "Primary: Negative electrode OCP entropic change [V.K-1]": dudt_mcmb,  # modified@240429

        "Secondary: Maximum concentration in negative electrode [mol.m-3]": 3.4466e5,  # modified@240429
        "Secondary: Initial concentration in negative electrode [mol.m-3]": 3.4466e5 * 0.792,
        # modified@240429
        "Secondary: Negative electrode diffusivity [m2.s-1]": ds_si_ls,  # modified@240429
        "Secondary: Negative electrode OCP [V]": silicon_ocp_Eeq_llh_origin_fit,  # modified@240514
        "Secondary: Negative electrode lithiation OCP [V]"
        "": silicon_ocp_lithiation_llh_origin_fit,  # modified@240513
        "Secondary: Negative electrode delithiation OCP [V]"
        "": silicon_ocp_delithiation_llh_origin_fit,  # modified@240513
        "Secondary: Negative electrode active material volume fraction": 0.02498,  # modified@240429
        "Secondary: Negative particle radius [m]": 6.252e-06 / 2,  # modified@240429
        "Secondary: Negative electrode exchange-current density [A.m-2]"
        "": silicon_ls_electrolyte_exchange_current_density,
        "Secondary: Negative electrode density [kg.m-3]": 2200.0,  # modified@240513
        "Secondary: Negative electrode OCP entropic change [V.K-1]": 0.0,
        "Primary: Negative electrode diffusivity activation energy": 28000,  #@240809
        "Secondary: Negative electrode diffusivity activation energy": 48000,  # @240809
        "Primary: Negative electrode exchange-current density activation energy": 57500,
        "Secondary: Negative electrode exchange-current density activation energy": 57500,
        # positive electrode
        "Positive electrode conductivity [S.m-1]": 3.5,  # modified@240513
        "Maximum concentration in positive electrode [mol.m-3]": 49637.0,  # modified@240429
        "Positive electrode diffusivity [m2.s-1]": ds_nac_ls,  # modified@240513
        "Positive electrode OCP [V]": nca_ls_ocp_new,  # modified@240513
        "Positive electrode porosity": 0.23702,  # modified@240429
        "Positive electrode active material volume fraction": 0.689,  # modified@240429
        "Positive particle radius [m]": 12.699e-06 / 2,  # modified@240429
        "Positive electrode Bruggeman coefficient (electrolyte)": 1.78,  # modified@240429
        "Positive electrode Bruggeman coefficient (electrode)": 0,
        "Positive electrode charge transfer coefficient": 0.5,
        "Positive electrode double-layer capacity [F.m-2]": 0.2,
        "Positive electrode exchange-current density [A.m-2]"
        "": nca_j0,
        "Positive electrode density [kg.m-3]": 4800.0,  # modified@240513
        "Positive electrode specific heat capacity [J.kg-1.K-1]": 1126.0,  # modified@240429
        "Positive electrode thermal conductivity [W.m-1.K-1]": 2.1,
        "Positive electrode OCP entropic change [V.K-1]": 0.0,
        "Positive electrode diffusivity activation energy": 35780,  # add by lou@240605
        "Positive electrode exchange-current density activation energy": 31000,  # add by lou@240710
        # separator
        "Separator porosity": 0.4413,  # modified@240513
        "Separator Bruggeman coefficient (electrolyte)": 2.24,  # modified@240513
        "Separator density [kg.m-3]": 397.0,
        "Separator specific heat capacity [J.kg-1.K-1]": 1126.0,  # modified@240429
        "Separator thermal conductivity [W.m-1.K-1]": 0.16,
        # electrolyte
        "Initial concentration in electrolyte [mol.m-3]": 1200.0,  # modified@240429
        "Cation transference number": 0.34,  # modified@240429
        "Thermodynamic factor": 1.0,
        "Electrolyte diffusivity [m2.s-1]": electrolyte_diffusivity_ls_46,  # modified@240513
        "Electrolyte conductivity [S.m-1]": electrolyte_conductivity_ls_46,  # modified@240513
        "Electrolyte diffusivity activation energy": 40000,  # add by lou@240605 not use
        "Electrolyte conductivity activation energy": 40000,  # add by lou@240605 not use
        # experiment
        "Reference temperature [K]": 298.15,
        "Total heat transfer coefficient [W.m-2.K-1]": 35.0,  # modified@240429
        "Ambient temperature [K]": 298.15,
        "Number of electrodes connected in parallel to make a cell": 14.0,  # modified@240429
        "Number of cells connected in series to make a battery": 1.0,  # modified@240429
        "Lower voltage cut-off [V]": 2.5,
        "Upper voltage cut-off [V]": 4.25,  # modified@240429
        "Open-circuit voltage at 0% SOC [V]": 2.5,
        "Open-circuit voltage at 100% SOC [V]": 4.25,  # modified@240429
        "Initial concentration in negative electrode [mol.m-3]": 31915.0 * 0.05,  # modified@240429
        "Initial concentration in positive electrode [mol.m-3]": 49637.0 * 0.13,  # modified@240429
        "Initial temperature [K]": 298.15,
        # citations
        "citations": ["Chen2020", "Ai2022"],
    }

testbench

import os
import matplotlib.pyplot as plt
import numpy as np
import pybamm
import timeit
from matplotlib import style
import pybamm
import ls_46_yk
import sys
import time
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

pybamm.set_logging_level("INFO")

start = timeit.default_timer()
model = pybamm.lithium_ion.DFN(
    {
        "particle phases": ("2", "1"),
        "open-circuit potential": (("single", "Wycisk"), "single"),
        "thermal": "lumped",
        "cell geometry": "pouch",
    }
)
param_dict = ls_46_yk.get_parameter_values()
param = pybamm.ParameterValues(param_dict)

C_rate = 0.1
capacity = param["Nominal cell capacity [A.h]"]
I_load = C_rate * capacity

t_eval = np.linspace(0, 36000, 1000)
exp = pybamm.Experiment([
    (
        "Discharge at 1C until 2.5 V (3 second period)",
        "Rest for 100 second (2 second period)",
        "Rest for 500 second (10 second period)",
        "Rest for 3000 second (50 second period)",
        "Charge at 0.3C until 4.25 V (5 second period)",
        "Rest for 100 second (2 second period)",
        "Rest for 500 second (10 second period)",
        "Rest for 3000 second (50 second period)",
    )
])
param["Current function [A]"] = I_load

sim = pybamm.Simulation(

    model,
    parameter_values=param,
    experiment=exp,
    solver=pybamm.CasadiSolver(dt_max=50,atol=1e-3,rtol=1e-3),
)
solution = sim.solve(calc_esoh=False)
stop = timeit.default_timer()

plt.rcdefaults()
plot = pybamm.QuickPlot(
    [solution
     ],
    [
        # "Negative electrode porosity",
        # "Positive electrode porosity",
        "Volume-averaged cell temperature [C]",
        # "Terminal voltage [V]",
        # "Negative lithium plating thickness [m]",
        # "R-averaged negative particle concentration [mol.m-3]",
        # "Negative total SEI thickness [m]",
        "Negative primary particle concentration [mol.m-3]",
        "Negative secondary particle concentration [mol.m-3]",
        "Positive particle concentration [mol.m-3]",
        "Current [A]",
        "Electrolyte concentration [mol.m-3]",
        "Terminal voltage [V]",
        "Negative electrode secondary hysteresis state",
        "X-averaged negative electrode secondary hysteresis state",
        "Negative electrode secondary hysteresis state distribution",
        "Sum of x-averaged negative electrode primary volumetric interfacial current densities [A.m-3]",
        "Sum of x-averaged negative electrode secondary volumetric interfacial current densities [A.m-3]",
        # "Negative electrode primary exchange current density [A.m-2]",
        # "Negative electrode secondary exchange current density [A.m-2]",
        # "Negative electrode potential [V]",
        # "Negative electrolyte potential [V]",
    ],
    time_unit="seconds",
    spatial_unit="um",
)
plot.dynamic_plot()

[llh929423719]

After reading the literature, I wondered if it was due to the current problem. I modified the current calculation method and it seemed to be improved. The following is my modified wycisk_ocp(mainly modified set_rhs):

#
# from Wycisk 2022
#
import pybamm
from . import BaseOpenCircuitPotential


class WyciskOpenCircuitPotential(BaseOpenCircuitPotential):
    """
    Class for open-circuit potential with hysteresis based on the approach outlined by Wycisk :footcite:t:'Wycisk2022'.
    This approach employs a differential capacity hysteresis state variable. The decay and switching of the hysteresis state
    is tunable via two additional parameters. The hysteresis state is updated based on the current and the differential capacity.
    """

    def get_fundamental_variables(self):
        domain, Domain = self.domain_Domain
        phase_name = self.phase_name
        h = pybamm.Variable(
            f"{Domain} electrode {phase_name}hysteresis state",
            domains={
                "primary": f"{domain} electrode",
                "secondary": "current collector",
            },
            bounds=(pybamm.Scalar(-1),pybamm.Scalar(1)),
        )
        return {
            f"{Domain} electrode {phase_name}hysteresis state": h,
        }

    def get_coupled_variables(self, variables):
        domain, Domain = self.domain_Domain
        phase_name = self.phase_name
        phase = self.phase

        if self.reaction == "lithium-ion main":
            T = variables[f"{Domain} electrode temperature [K]"]
            h = variables[f"{Domain} electrode {phase_name}hysteresis state"]

            # For "particle-size distribution" models, take distribution version
            # of c_s_surf that depends on particle size.
            domain_options = getattr(self.options, domain)
            if domain_options["particle size"] == "distribution":
                sto_surf = variables[
                    f"{Domain} {phase_name}particle surface stoichiometry distribution"
                ]
                # If variable was broadcast, take only the orphan
                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
                    T, pybamm.Broadcast
                ):
                    sto_surf = sto_surf.orphans[0]
                    T = T.orphans[0]
                T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
                h = pybamm.PrimaryBroadcast(h, [f"{domain} {phase_name}particle size"])
            else:
                sto_surf = variables[
                    f"{Domain} {phase_name}particle surface stoichiometry"
                ]
                # If variable was broadcast, take only the orphan
                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
                    T, pybamm.Broadcast
                ):
                    sto_surf = sto_surf.orphans[0]
                    T = T.orphans[0]

            variables[
                f"{Domain} electrode {phase_name}hysteresis state distribution"
            ] = h

            # Bulk OCP is from the average SOC and temperature
            sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"]
            c_scale = self.phase_param.c_max
            variables[f"Total lithium in {phase} phase in {domain} electrode [mol]"] = (
                sto_bulk * c_scale
            )  # c_s_vol * L * A

            ocp_surf_eq = self.phase_param.U(sto_surf, T)
            variables[f"{Domain} electrode {phase_name}equilibrium OCP [V]"] = (
                ocp_surf_eq
            )

            T_bulk = pybamm.xyz_average(pybamm.size_average(T))
            ocp_bulk_eq = self.phase_param.U(sto_bulk, T_bulk)
            variables[f"{Domain} electrode {phase_name}bulk equilibrium OCP [V]"] = (
                ocp_bulk_eq
            )

            inputs = {f"{Domain} {phase_name}particle stoichiometry": sto_surf}
            lith_ref = pybamm.FunctionParameter(
                f"{self.phase_param.phase_prefactor}{Domain} electrode lithiation OCP [V]",
                inputs,
            )
            delith_ref = pybamm.FunctionParameter(
                f"{self.phase_param.phase_prefactor}{Domain} electrode delithiation OCP [V]",
                inputs,
            )
            H = abs(delith_ref - lith_ref) / 2
            variables[f"{Domain} electrode {phase_name}OCP hysteresis [V]"] = H

            # determine dQ/dU
            if phase_name == "":
                Q_mag = variables[f"{Domain} electrode capacity [A.h]"]
            else:
                Q_mag = variables[
                    f"{Domain} electrode {phase_name}phase capacity [A.h]"
                ]

            dU = self.phase_param.U(sto_surf, T_bulk).diff(sto_surf)
            dQdU = Q_mag / dU
            variables[
                f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
            ] = dQdU

            H_x_av = pybamm.x_average(H)
            h_x_av = pybamm.x_average(h)
            variables[f"X-averaged {domain} electrode {phase_name}hysteresis state"] = (
                h_x_av
            )

            # check if psd
            if domain_options["particle size"] == "distribution":
                # should always be true
                if f"{domain} particle size" in sto_surf.domains["primary"]:
                    # check if MPM Model
                    if "current collector" in sto_surf.domains["secondary"]:
                        ocp_surf = ocp_surf_eq + H_x_av * h_x_av
                    # must be DFN with PSD model
                    elif (
                        f"{domain} electrode" in sto_surf.domains["secondary"]
                        or f"{domain} {phase_name}particle size"
                        in sto_surf.domains["primary"]
                    ):
                        ocp_surf = ocp_surf_eq + H * h
            # must not be a psd
            else:
                ocp_surf = ocp_surf_eq + H * h

            H_s_av = pybamm.size_average(H_x_av)
            h_s_av = pybamm.size_average(h_x_av)

            ocp_bulk = ocp_bulk_eq + H_s_av * h_s_av

            dUdT = self.phase_param.dUdT(sto_surf)

        variables.update(self._get_standard_ocp_variables(ocp_surf, ocp_bulk, dUdT))
        return variables

    def set_rhs(self, variables):
        domain, Domain = self.domain_Domain
        phase_name = self.phase_name
        domain_param = self.domain_param
        L0 = domain_param.L
        param = self.param
        epss =  variables[
            f"{Domain} electrode {phase_name}active material volume fraction"
        ]
        a = variables[
            f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]"
        ]
        acc=param.A_cc
        current = variables[
            f"{Domain} electrode {phase_name}interfacial current density [A.m-2]"
        ]
        # current = variables[
        #     f"{Domain} electrode {phase_name}volumetric interfacial current density [A.m-3]"
        # ]
        # current = pybamm.x_average(current)
        # check if composite or not
        if phase_name != "":
            Q_cell = variables[f"{Domain} electrode {phase_name}phase capacity [A.h]"]
        else:
            Q_cell = variables[f"{Domain} electrode capacity [A.h]"]

        dQdU = variables[
            f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
        ]
        dQdU = dQdU.orphans[0]
        K = self.phase_param.hysteresis_decay
        K_x = self.phase_param.hysteresis_switch
        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]

        # dhdt = (
        #     K * (current / (Q_cell * (dQdU**K_x))) * (1 - pybamm.sign(current) * h)
        # )  #! current is backwards for a halfcell
        dhdt = (
            K * (current*L0*acc*a*epss / (Q_cell * (dQdU**K_x))) * (1 - pybamm.sign(current) * h)
        )
        self.rhs[h] = dhdt

    def set_initial_conditions(self, variables):
        domain, Domain = self.domain_Domain
        phase_name = self.phase_name
        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
        self.initial_conditions[h] = pybamm.Scalar(-1)

[ejfdickinson]

@llh929423719 Could you explain what your changes mean in mathematical terms? I see that you corrected some of the unit peculiarities for the implied units of the decay rate. What are the other implications?

@mleot I was working with this submodel recently (thank you for the very helpful addition!). The Plett model and the Wycisk variation (edit: Wycisk makes a meal of the units, see below) defines the decay rate $\gamma$ as (effectively) $\partial \ln h / \partial x_\mathrm{Li}$. I noticed that your $K$ differs as:

$K = \frac{\gamma}{3600} a_\mathrm{vol}V_\mathrm{el}$

where $a_\mathrm{vol}$ is volumetric surface area and $V_\mathrm{el}$ is electrode volume. Hence the decay rate is not unitless but rather in its present form has implied units m^2 (for K_x == 0 - if K_x != 0 then normalising dQdU to a reference value would be tidier unitwise). I think that revising the expression to use the standard decay rate definition would be better but now raises backwards-compatibility issues.

@llh929423719 I notice that you also changed the initial condition for h in set_initial_conditions(). @mleot It would be useful to be able to set the initial hysteresis state via the parameter values, like other variable initialisation.

[llh929423719]

@ejfdickinson The formula in the Wycisk paper is $\frac{\mathrm{d}h(z,t)}{\mathrm{d}t}=\left(\frac{k(z)\cdot I(t)}{Q_{\mathrm{Cell}}}\right)\left(1-\mathrm{sgn}\left(\frac{\mathrm{d}z(t)}{\mathrm{d}t}\right)h(z,t)\right)$ where $Q_{cell}$ is the capacity of the cell and $𝐼(𝑡)$ is the current. As in pybamm, the current $𝐼(𝑡)$ uses the interfacial current density [A.m-2] distribution(which is relatively large), but $Q_{cell}$ uses electrode phase capacity [A.h](which is relatively small), so I multiply the current density by the active material volume fraction $\varepsilon_{s}$ and the total electrode volume

[ejfdickinson]

Thanks @llh929423719 - I think your approach makes the most sense so that $k$ has natural units (unitless) - take the ratio of current to capacity on the same scale, as you would for the full cell.

Ultimately your change just alters the numerical value of the decay rate though, rather than changing the structure of the equations. So you could achieve the same just be scaling the decay rate, right? Do you see the same improvement if you keep the 24.5 release code and scale the decay rate parameter magnitude?

[llh929423719]

Thanks @ejfdickinson , keeping the 24.5 release code and scale the decay rate parameter magnitude has the same effect.

[mleot]

Hi @ejfdickinson, thanks for these comments. You make a fair point - as I was writing the code for this, I realized that realistically the current which is impacting the hysteresis should be measured locally at the particle interface, instead of a lumped electrode-level current. I believe I made an oversight by not reflecting and adjusting the units accordingly. I think there could be an argument for having K be a united quantity, although I have not thought about it enough to make any arguments there.

Further, I do like the suggestion of being able to specify the initial hysteresis state condition in the parameter set.

@llh929423719 I think the modification that you made above is not a bad idea. My only criticism would be that I don't believe the electrode thickness needs to be included in the calculation. I believe the alternate form below would make sure K is unitless.

dhdt = (
            K * (current / ((Q_cell / (acc * epss * a)) * (dQdU**K_x))) * (1 - pybamm.sign(current) * h)
        )
        self.rhs[h] = dhdt

I suspect the issue you encountered may have been from modifying the initial hysteresis state to a value very close to it's limit (I don't believe the hysteresis state variable should every actually reach 1 or -1). In the original bug, did you modify the initial condition at all?

[llh929423719]

@mleot thanks for reply, I modify the initial condition just for fit the experiment data (cell is fully charged before discharge), but setting initial condition to zero does not make improvement. I take the electrode thickness into calculation just because the $Q_{cell}$ is calaulated cross the whole electrode, and $I(t) [A.m^{-2}] \cdot L_{0} [m] \cdot A_{cc} [m^2] \cdot a [m^{-1}]$ has the same unit $A$ as in Wycisk's paper Equ[9]
$\frac{\mathrm{d}h(z,t)}{\mathrm{d}t}=K\cdot\left(\frac{I(t)}{Q_{\mathrm{Cell}}\cdot C_{\mathrm{diff}}(z)}\right)\left(1-\mathrm{sgn}\left(\frac{\mathrm{d}z(t)}{\mathrm{d}t}\right)h(z,t)\right)$.
But maybe I have misunderstood some of the variables in the pybamm model.

[ejfdickinson]

@mleot @llh929423719

To keep electrode thickness out and to keep the focus on particle-scale properties, would it be better to normalise to the volumetric charge capacity of the pure material instead? That is:

$\frac{\partial h}{\partial t} \propto K\cdot\left(\frac{i_\mathrm{surf}a_\mathrm{vol}}{\phi_\mathrm{act}Fc_\mathrm{max}}\right)$

This means that the formulation is considering the particle-scale ratio of current to charge capacity. This is directly compatible with the need to normalise only by the capacity of the particle phase under consideration.

At present you are getting a whole electrode property (product of volumetric charge capacity with electrode dimensions) and then having to divide back out by the electrode dimensions to get back to the volumetric charge capacity. So just grab the voilumetric charge capacity directly?

edit: removed discussion of single-phase vs total capacity which is already handled in the current code

[rtimms]

I’m in favour of updating this formulation. It will be helpful to have reusable and easily interpretable parameters, especially as we look to integrate more hysteresis models.

[ejfdickinson]

@rtimms Would it be worth opening a discussion for this? I've used the WyciskOpenCircuitPotential feature a fair amount in the last couple of months and have a number of suggestions, but better that these are aired before committing too unilaterally.

[rtimms]

I’m happy to use this issue to discuss the changes. Could you post your suggestions and we can figure out a plan to move forward?

[ejfdickinson]

@rtimms @mleot @llh929423719

Suggestions for improvements:

• Solve a general one-state equation:

$\frac{\partial h}{\partial t} = \gamma\cdot\left(\frac{i_\mathrm{surf}a_\mathrm{vol}}{\phi_\mathrm{act}Fc_\mathrm{max}}\right)\cdot \left(1 - \mathrm{sgn}(i_\mathrm{surf})h\right)$

for $i_\mathrm{surf}$ the interfacial phase-specific current density (anodic current density positive).

• Allow a scalar input or arbitrary functional input for $\gamma$, or allow $\gamma$ to be computed using the specific Wycisk formula with parameters $K$, $n$ and the OCP:

$\gamma = K \cdot \left(\frac{\frac{\partial q_\mathrm{vol}}{\partial U}}{\left.\frac{\partial q_\mathrm{vol}}{\partial U}\right|_\mathrm{ref}}\right)^{-n}$

where the evaluation at a reference point (implicitly at $x_\mathrm{Li}=x_\mathrm{Li,ref}$ which in turn defines $K$ as the decay rate at the reference stoichiometry) is to ensure that we take an arbitrary power only of a unitless quantity.

• Allow general (particle-wise) initialisation of the hysteresis state, not necessarily = 0.

[mleot]

@ejfdickinson Thanks for the suggestions.

Regarding your first suggestion: I think is a neat implementation of this portion of the equation.

Regarding the second suggestion: Is this specific equation referenced in the Wycisk paper? This does seem like a better way to implement the differential capacity with the power term. How do you suggest the reference point be implemented? An additional parameter for a specific stoichiometry from which the reference value is extracted?

I also like the idea of being able to specify a custom function, but the implementation for that seems potentially challenging. It might be nice in an alternate implementation to have a $$\gamma$$ which accounts for rate invariance (mentioned in the Baker 2017 paper)link.

The last suggestion makes sense and I have no arguments against it.

[ejfdickinson]

@mleot I don't think this specific equation is anywhere in Wycisk but it's just a factorisation of a scalar into the product of two parts with the effect of (a) tidying up the units and (b) correspondingly making more explicit and comprehensible the physical meaning of the parameters. Wycisk et al. would be far from the first authors to fail to introduce the tidiest mathematical expression of their own good idea, particularly when a general-purpose modelling approach is envisaged.

Setting the specific stoichiometry at which the reference value is extracted is sensible to me. That just changes the definition of K from "decay rate" to "decay rate at xLi = ???", where the Wycisk theory then uses (dq/dU)^(-n) as the scaling to get decay rate at other xLi.

[ejfdickinson]

I don't completely know enough about PyBaMM's architecture to know the difficulties of arbitrary funtional input for $\gamma$ once this parameter is exposed in the form I wrote it.

At least a scalar input option would be very useful as an 'opt-out' from the specific Wycisk relation. I have some models at the moment where I achieve scalar input by setting n == 0, but the existing implementation still always tries to evaluate $dq/dU$ so some shenanigans are needed to work around cases where _function_diff is undefined for the open-circuit potential expression (e.g. interpolants).

[mleot]

Referring to @rtimms for advice on difficulty of implementation

[rtimms]

If you make something a FunctionParameter, it can be a function or a parameter. The only thing that is fixed is the arguments of the function. What variables could $\gamma$ reasonably depend on?

An alternative approach would be to have a general model that allows for scalar $\gamma$ (or covers some sensible functional forms), that defines a _get_rate_constant(variables) method that returns the FunctionParameter. Then any extension can simply subclass and implement their own _get_rate_constant(variables) that can depend on the variables in any way.

In that case, we could have a general OneStateHysteresis model and the Wycisk model could subclass and compute $\gamma$ from dq/dU.

[ejfdickinson]

@rtimms @mleot I guess in the first instance we would expect this to be parameterised depending on lithiation extent (cs / csmax) and temperature, like diffusivity is.

Anything more would be in the realm of future custom submodels depending on academic developments in understanding of the physical phenomena that dictate the magnitude of $\gamma$. In reality I think it's an inter-particle phenomenon for LFP and probably fairly mechanically controlled for Si, so there's loads of extra physics if you want to get beyond semi-empiricism. That's definitely for the future!

[rtimms]

Sounds good! Would either of you like to take this on? IMO we should just do it as a breaking change and make sure we document how to update parameters that users have for the current implementation to match the new one.

[ejfdickinson]

I would recommend that @mleot continues actual development of this class (if he is willing!), but I would be happy to assist with code review and testing.

[rtimms]

If helpful, I’m happy to help with setting up the new class structure. Then I’ll leave the details to @mleot (again, if you’re willing!).

[mleot]

I would be happy to help! @rtimms If you want to set up the new class structure that would be very helpful!