How to change a parameter inside the function graphite_cracking_rate_Ai2020

[Akila1993]

Hi Everyone & @DrSOKane

I am using the degradation models with the OKane2022 parameter set, and for my analysis, I need to change the k_cr =Paris’ law cracking rate to different rates for example 1x,10x, and 30x. The value is inside OKane2022.py as the function graphite_cracking_rate_Ai2020.

Is there a simple way of doing this, could someone guide me, please? Is the following code can be used.

Thank you.

import pybamm
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft

# Define the function for graphite cracking rate as in the original code
def graphite_cracking_rate_Ai2020(T_dim, k_cr=10.0):
    """
    Graphite particle cracking rate as a function of temperature and scaling factor for k_cr.
    
    Parameters
    ----------
    T_dim: :class:`pybamm.Symbol`
        temperature in Kelvin.
    k_cr: float, optional
        Cracking rate scaling factor. Default is 1.0.
    
    Returns
    -------
    k_cr: :class:`pybamm.Symbol`
        Cracking rate [m/(Pa.m^0.5)^m_cr].
    """
    # Apply Arrhenius law for temperature dependence
    Eac_cr = 0  # Activation energy for cracking (can be updated if available)
    arrhenius = pybamm.exp(Eac_cr / pybamm.constants.R * (1 / T_dim - 1 / 298.15))
    return k_cr * arrhenius  # 'k_cr' will be scaled dynamically

# Original function for degradation parameters
pp = pybamm.ParameterValues("OKane2022")

def get_chen2020_parameters():
    parameter_values = pybamm.ParameterValues("Chen2020")
    parameters = parameter_values
    return parameters

def get_okane2022_degradation_parameters():
    degradation_parameters = {
        'Dead lithium decay constant [s-1]': 1e-06,
        'Dead lithium decay rate [s-1]': pp["Dead lithium decay rate [s-1]"],
        'Initial plated lithium concentration [mol.m-3]': 0.0,
        'Lithium metal partial molar volume [m3.mol-1]': 1.3e-05,
        'Lithium plating kinetic rate constant [m.s-1]': 1e-09,
        'Lithium plating transfer coefficient': 0.65,
        'Negative electrode LAM constant exponential term': 2.0,
        'Negative electrode LAM constant proportional term [s-1]': 2.7778e-07,
        "Negative electrode Paris' law constant b": 1.12,
        "Negative electrode Paris' law constant m": 2.2,
        "Negative electrode Poisson's ratio": 0.3,
        "Negative electrode Young's modulus [Pa]": 15000000000.0,
        'Negative electrode cracking rate': pp["Negative electrode cracking rate"],  # Function
        'Negative electrode critical stress [Pa]': 60000000.0,
        'Negative electrode initial crack length [m]': 2e-08,
        'Negative electrode initial crack width [m]': 1.5e-08,
        'Negative electrode number of cracks per unit area [m-2]': 3180000000000000.0,
        'Negative electrode partial molar volume [m3.mol-1]': 3.1e-06,
        'Negative electrode reference concentration for free of deformation [mol.m-3]': 0.0,
        'Negative electrode volume change': pp["Negative electrode volume change"],
        'Positive electrode LAM constant exponential term': 2.0,
        'Positive electrode LAM constant proportional term [s-1]': 2.7778e-07,
        'Positive electrode OCP entropic change [V.K-1]': 0.0,
        "Positive electrode Paris' law constant b": 1.12,
        "Positive electrode Paris' law constant m": 2.2,
        "Positive electrode Poisson's ratio": 0.2,
        "Positive electrode Young's modulus [Pa]": 375000000000.0,
        'Positive electrode active material volume fraction': 0.665,
        'Positive electrode cracking rate': pp["Positive electrode cracking rate"],  # Function
        'Positive electrode critical stress [Pa]': 375000000.0,
        'Positive electrode initial crack length [m]': 2e-08,
        'Positive electrode initial crack width [m]': 1.5e-08,
        'Positive electrode number of cracks per unit area [m-2]': 3180000000000000.0,
        'Positive electrode partial molar volume [m3.mol-1]': 1.25e-05,
        'Positive electrode reference concentration for free of deformation [mol.m-3]': 0.0,
        'Positive electrode volume change': pp["Positive electrode volume change"],
        'Typical plated lithium concentration [mol.m-3]': 1000.0,
        'Exchange-current density for stripping [A.m-2]': pp["Exchange-current density for stripping [A.m-2]"],
        'Exchange-current density for plating [A.m-2]': pp["Exchange-current density for plating [A.m-2]"],
    }
    return degradation_parameters

def updated_parameters(k_cr_factor=1.0):
    # Get Chen2020 parameters and Okane2022 degradation parameters
    chen2020_params = get_chen2020_parameters()
    okane2022_degradation_params = get_okane2022_degradation_parameters()

    # Update the cracking rate for both positive and negative electrodes dynamically
    okane2022_degradation_params["Positive electrode cracking rate"] = graphite_cracking_rate_Ai2020(pybamm.Scalar(298.15), k_cr=k_cr_factor)
    okane2022_degradation_params["Negative electrode cracking rate"] = graphite_cracking_rate_Ai2020(pybamm.Scalar(298.15), k_cr=k_cr_factor)

    # Merge Okane2022 degradation parameters with Chen2020 parameters
    chen2020_params.update(okane2022_degradation_params, check_already_exists=False)

    return chen2020_params

# Main part of the code
k_cr_factors = [10]  # Factors for scaling k_cr
impedance_data_all = {}

# Loop over different k_cr factors
for k_cr_factor in k_cr_factors:
    print(f"Running simulation with k_cr factor = {k_cr_factor}")

    # Update parameters with the new k_cr scaling
    updated_parameters_dict = updated_parameters(k_cr_factor)
    
    # Define the experiment and model
    
    # Define model with configuration
    model = pybamm.lithium_ion.DFN(
    {   
        "thermal": "lumped", #lumped thermal submodel 
        "surface form": "differential", #Adding capacitance to the model 
        "SEI":"solvent-diffusion limited", #"reaction limited",
        "SEI porosity change": "true",
        "particle mechanics":"swelling and cracking",
        "SEI on cracks": "true",
    }
)

var_pts = {
    "x_n": 5,
    "x_s": 5,
    "x_p": 5,
    "r_n": 30,
    "r_p": 30,
}
cycle_numbers = [1]  # Example cycle numbers, modify as needed

for cycle_number in cycle_numbers:
    exp = pybamm.Experiment(
            [
                (
                    "Discharge at 1C until 2.5 V",  # ageing cycles
                    "Charge at 0.6C until 4.2 V (5 minute period)",
                    "Hold at 4.2 V until C/100 (5 minute period)",
                )
            ]
            * cycle_number
        )

    sim = pybamm.Simulation(
    model,
    parameter_values=updated_parameters_dict,
    experiment=exp,
    solver=pybamm.CasadiSolver(mode='safe without grid'),
    var_pts=var_pts,
    )

sol = sim.solve()

# Set initial conditions of aged_model from the solution and EIS Running
aged_model = model.set_initial_conditions_from(sol)

# Impedance calculation for each frequency
impedances_time = []
frequencies = np.logspace(-4, 2, 30)
I = 1e-3  # Current amplitude for EIS
number_of_periods = 20
samples_per_period = 20

# Current function for EIS
def current_function(t):
    return I * pybamm.sin(2 * np.pi * pybamm.InputParameter("Frequency [Hz]") * t)

updated_parameters_dict["Current function [A]"] = current_function

sim_aged = pybamm.Simulation(
        aged_model,
        parameter_values=updated_parameters_dict,
        solver=pybamm.CasadiSolver(mode='safe without grid'),
        var_pts=var_pts,
    )

# Impedance calculation for each frequency
impedances_time = []
for frequency in frequencies:
        period = 1 / frequency
        dt = period / samples_per_period
        t_eval = np.array(range(0, 1 + samples_per_period * number_of_periods)) * dt
        sol = sim_aged.solve(t_eval, inputs={"Frequency [Hz]": frequency}, initial_soc=0.9)

        # Extract final two periods of the solution
        current = sol["Current [A]"].entries[-3 * samples_per_period - 1 :]
        voltage = sol["Terminal voltage [V]"].entries[-3 * samples_per_period - 1 :]

        # FFT
        current_fft = fft(current)
        voltage_fft = fft(voltage)

        # Get index of first harmonic
        idx = np.argmax(np.abs(current_fft))
        impedance = -voltage_fft[idx] / current_fft[idx]
        impedances_time.append(impedance)

    # Store the impedance data for the current cycle
impedance_data_all[cycle_number] = impedances_time

# Plot Nyquist plots for all cycles
plt.figure(figsize=(8, 8))
for cycle_number, impedances_time in impedance_data_all.items():
    real_part = np.real(impedances_time)
    imaginary_part = np.imag(impedances_time)
    plt.plot(real_part, -imaginary_part, marker='o', linestyle='-', label=f'{cycle_number} cycles')

plt.title('Electrochemical Impedance Spectroscopy (EIS) - Nyquist Plot')
plt.xlabel('Real(Z) [Ohm]')
plt.ylabel('-Imaginary(Z) [Ohm]')
plt.legend()
plt.grid(True)
plt.show()

[rtimms]

You can make it a Parameter object and then treat it like any other parameter

def graphite_cracking_rate_Ai2020(T_dim):
    Eac_cr = 0  # Activation energy for cracking (can be updated if available)
    arrhenius = pybamm.exp(Eac_cr / pybamm.constants.R * (1 / T_dim - 1 / 298.15))
    return pybamm.Parameter("k_cr") * arrhenius  

[Akila1993]

Hi @rtimms,

Thank you for your quick response, I have made the changes to the code, and please check and let me know if it is the correct approach.

Thank you.

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft
import pybamm

# Set random seed
np.random.seed(0)

# Original function for degradation parameters
pp = pybamm.ParameterValues("OKane2022")

# Function to get Chen2020 parameters
def get_chen2020_parameters():
    parameter_values = pybamm.ParameterValues("Chen2020")
    return parameter_values

# Function to get OKane2022 degradation parameters with k_cr included
def get_okane2022_degradation_parameters():
    degradation_parameters = {
        'Dead lithium decay constant [s-1]': 1e-06,
        'Dead lithium decay rate [s-1]': pp["Dead lithium decay rate [s-1]"],
        'Initial plated lithium concentration [mol.m-3]': 0.0,
        'Lithium metal partial molar volume [m3.mol-1]': 1.3e-05,
        'Lithium plating kinetic rate constant [m.s-1]': 1e-09,
        'Lithium plating transfer coefficient': 0.65,
        'Negative electrode LAM constant exponential term': 2.0,
        'Negative electrode LAM constant proportional term [s-1]': 2.7778e-07,
        "Negative electrode Paris' law constant b": 1.12,
        "Negative electrode Paris' law constant m": 2.2,
        "Negative electrode Poisson's ratio": 0.3,
        "Negative electrode Young's modulus [Pa]": 15000000000.0,
        'Negative electrode cracking rate': pp["Negative electrode cracking rate"],  # Function
        'Negative electrode critical stress [Pa]': 60000000.0,
        'Negative electrode initial crack length [m]': 2e-08,
        'Negative electrode initial crack width [m]': 1.5e-08,
        'Negative electrode number of cracks per unit area [m-2]': 3180000000000000.0,
        'Negative electrode partial molar volume [m3.mol-1]': 3.1e-06,
        'Negative electrode reference concentration for free of deformation [mol.m-3]': 0.0,
        'Negative electrode volume change': pp["Negative electrode volume change"],
        'Positive electrode LAM constant exponential term': 2.0,
        'Positive electrode LAM constant proportional term [s-1]': 2.7778e-07,
        'Positive electrode OCP entropic change [V.K-1]': 0.0,
        "Positive electrode Paris' law constant b": 1.12,
        "Positive electrode Paris' law constant m": 2.2,
        "Positive electrode Poisson's ratio": 0.2,
        "Positive electrode Young's modulus [Pa]": 375000000000.0,
        'Positive electrode active material volume fraction': 0.665,
        'Positive electrode cracking rate': pp["Positive electrode cracking rate"],  # Function
        'Positive electrode critical stress [Pa]': 375000000.0,
        'Positive electrode initial crack length [m]': 2e-08,
        'Positive electrode initial crack width [m]': 1.5e-08,
        'Positive electrode number of cracks per unit area [m-2]': 3180000000000000.0,
        'Positive electrode partial molar volume [m3.mol-1]': 1.25e-05,
        'Positive electrode reference concentration for free of deformation [mol.m-3]': 0.0,
        'Positive electrode volume change': pp["Positive electrode volume change"],
        'Typical plated lithium concentration [mol.m-3]': 1000.0,
        'Exchange-current density for stripping [A.m-2]': pp["Exchange-current density for stripping [A.m-2]"],
        'Exchange-current density for plating [A.m-2]': pp["Exchange-current density for plating [A.m-2]"],
        "k_cr": pybamm.Parameter("k_cr"),  # Define as a parameter object
    }
    return degradation_parameters

# Graphite cracking rate function
def graphite_cracking_rate_Ai2020(T_dim):
    """
    Calculates the graphite cracking rate using an Arrhenius-like relationship.
    T_dim: Dimensional temperature (in Kelvin)
    """
    Eac_cr = 0  # Activation energy for cracking [J/mol], update if needed
    arrhenius = pybamm.exp(Eac_cr / pybamm.constants.R * (1 / T_dim - 1 / 298.15))
    return pybamm.Parameter("k_cr") * arrhenius

# Function to update parameters
def update_parameters():
    # Load Chen2020 parameters
    chen2020_params = get_chen2020_parameters()

    # Load OKane2022 degradation parameters
    okane2022_degradation_params = get_okane2022_degradation_parameters()

    # Update Chen2020 parameters with OKane2022 degradation parameters
    chen2020_params.update(okane2022_degradation_params, check_already_exists=False)

    # Add graphite cracking rate function
    chen2020_params["k_cr"] = 10*3.9e-20

    return chen2020_params

# Update the parameters
updated_parameters = update_parameters()

# Define the model with configuration
model = pybamm.lithium_ion.DFN(
    {   
        "thermal": "lumped",  # Lumped thermal submodel
        "surface form": "differential",  # Adding capacitance to the model
        "SEI": "solvent-diffusion limited",
        "SEI porosity change": "true",
        "particle mechanics": "swelling and cracking",
        "SEI on cracks": "true",
    }
)

# Variable points
var_pts = {
    "x_n": 5,
    "x_s": 5,
    "x_p": 5,
    "r_n": 30,
    "r_p": 30,
}


# Simulation parameters
cycle_numbers = [1]
frequencies = np.logspace(-4, 2, 30)


# Store results for Nyquist plot
impedance_data = {}

# Loop through each cycle count
for cycle_number in cycle_numbers:
    exp = pybamm.Experiment(
        [
            (
                "Discharge at 1C until 2.5 V",  # Ageing cycles
                "Charge at 0.6C until 4.2 V (5 minute period)",
                "Hold at 4.2 V until C/100 (5 minute period)",
            )
        ]
        * cycle_number
    )

    sim = pybamm.Simulation(
        model,
        parameter_values=updated_parameters,
        experiment=exp,
        solver=pybamm.CasadiSolver(mode='safe without grid'),
        var_pts=var_pts,
    )

    sol = sim.solve()

   

[rtimms]

BTW for EIS simulations I'd recommend taking a look at pybamm-eis