How to add experimental protocol to custom models

[suryankolluri]

Hello,
How to add experimental protocol, If I am writing custom battery model and simulating through pybamm casadi solver?

A sample Single Particle Model with parabolic approximation is included for reference. Thanks, Surya.

`
"""
@author: Suryanarayana Kolluri
Date Mon Jun 27 2022
Model: Isothermal Single Particle Model with Parabolic Profile Approximation for solid phase concentration
"""

import pybamm as pb
import tests
import numpy as np
import os
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.colors as mcolors
import matplotlib.patches as mpatch
import seaborn as sns
import math
from pprint import pprint
import time
import csv
plt.close("all")

1. PyBAMM settings

1.1 loading base model

model = pb.BaseModel()

1.2 Defining variables that are present in DAE's

cspavg = pb.Variable("cspavg") # 1. Positive electrode Average solid phase concentation
csnavg = pb.Variable("csnavg") # 2. Negative electrode Average solid phase concentation
Q = pb.Variable("Q") # 3. Charge stored (while charging) or charge withdrawn (while discharging)
csp = pb.Variable("csp") # 4. Positive electrode solid phase surface concentration
csn = pb.Variable("csn") # 5. Negative electrode solid phase surface concentration
phip = pb.Variable("phip") # 6. Positive electrode potential
phin = pb.Variable("phin") # 7. Negative electrode potential
v = pb.Variable("v") # 8. Cell Potential
iapp=pb.Variable("iapp") # 9. Input current density (A/m^2)

model.variables = {"cspavg":cspavg,
"csnavg":csnavg,
"Q":Q,
"csp":csp ,
"csn":csn ,
"phip":phip ,
"phin":phin ,
"v":v,
"iapp":iapp}

2. Parameters

F = 96487 # Faraday's constant
R = 8.314472 # Universal Gas Constant
T = 273.15+25 # Cell Temperature

Rpp = 2E-6 # Positive particle radius [m]
Rpn = 2E-6 # Negative particle radius [m]

ep = 0.385 # Positive electrode porosity
es = 0.724 # Separator porosity
en = 0.485 # Negative electrode porosity

efp = 0.025 # Positive electrode filler fraction
efn = 0.0326 # Negative electrode filler fraction

lp = 80E-6 # Positive electrode thickness [m]
ls = 25E-6 # Separator thickness [m]
ln = 88E-6 # Negative electrode thickness [m]

ctp = 51554 # Maximum concentration in positive electrode [mol.m-3]
ctn = 30555 # Maximum concentration in negative electrode [mol.m-3]
c0 = 1000 # Typical electrolyte concentration [mol.m-3]

t1 = 0.38 # Cation transference number

kp = 2.334e-11 # Positive electrode reaction rate constant [m2.5/mol0.5 s]
kn = 5.031e-11 # Negative electrode reaction rate constant [m2.5/mol0.5 s]

Dsp = 1.0E-14 # Positive electrode diffusivity [m2.s-1]
Dsn = 3.9E-14 # Negative electrode diffusivity [m2.s-1]

socp = 25751/51554.0 # Initial soc in positive electrode
socn = 26128/30555.0 # Initial soc in negative electrode

Particle surface area per unit volume

ap=3/Rpp*(1-ep-efp) # Positive electrode
an=3/Rpn*(1-en-efn) # Negative electrode

Calculating 1C-rate current density

socpmax = 0.98
iapp_pos = (1-efp-ep)ctp(socpmax-socp)Flp/3600
socnmin = 0.01
iapp_neg = (1-efn-en)ctn(socn-socnmin)Fln/3600
iapp_1c = min(iapp_pos,iapp_neg)

Open circuit potential of positive electrode

def OCV_pos(sto):
Up = (-4.656 + 88.669sto**2 - 401.119sto4 + 342.909*sto6 - 462.4710sto**8 + 433.434sto10)/(-1.0 + 18.933*sto2 - 79.532sto**4 + 37.311sto6 - 73.083*sto8 + 95.960*sto**10)
return Up

Open circuit potential of negative electrode

def OCV_neg(sto):
Un = 0.7222 + 0.13870sto + 0.029sto0.5 - 0.0172/sto + 0.0019/(sto1.5) + 0.2808pb.exp(0.9 - 15sto) - 0.7984pb.exp(0.44650sto - 0.4108)
return Un

Overpotential Terms in both electrodes

etap = phip-OCV_pos(csp)
etan = phin-OCV_neg(csn)

Input current

crate = -1

iapp = iapp_1c*crate

Time

tmax= 3600/abs(crate) # 3600/abs(crate) # 3600/abs(crate) # s, max time

# nt_init= 101 # time steps in exponential ramp of current to curr

dt = 5

nt = 200
t_eval = np.linspace(0, tmax, nt)

Butler Volmer Kinetics

jp = kpc0**0.5ctp*(1-csp)0.5csp0.5*(pb.exp(0.5Fetap/(RT))-pb.exp(-0.5Fetap/(RT)))
jn = kn*c00.5ctn(1-csn)**0.5csn0.5*(pb.exp(0.5Fetan/(RT))-pb.exp(-0.5Fetan/(RT)))

4. Defining ODE's and AE's

model.rhs = {cspavg:-3iapp/ap/lp/F/ctp/Rpp,
csnavg:3iapp/an/ln/F/ctn/Rpn,
Q:iapp/3600}

model.algebraic = {csp:(Dspctp/Rpp)(csp-cspavg) + iapp/5/ap/lp/F,
csn:(Dsnctn/Rpn)(csn-csnavg) - iapp/5/an/ln/F,
phip:jp-iapp/ap/lp/F,
phin:jn+iapp/an/ln/F,
v:v-phip+phin,
iapp:iapp-iapp_1c*crate}

model.initial_conditions = {cspavg:socp,
csnavg:socn,
Q:0,
csp:socp,
csn:socn,
phip:OCV_pos(socp),
phin:OCV_neg(socn).value,
v:OCV_pos(socp) - OCV_neg(socn).value,
iapp:iapp_1c*crate}

model.events +=[pb.Event("voltagep cutoff", v-2.8,
event_type=pb.EventType.TERMINATION)]

Since we are sovling set of DAE's related to isothermal SPM, and solid phase concentration in radial direction is approximated using parabolic approximation.

1.1 Create solver

dae_solver = pb.CasadiSolver(mode="safe",rtol=1e-4,atol=1e-6,dt_max=100) # use safe mode so that solver stops at events
disc = pb.Discretisation()
disc.process_model(model);
start_simtime = time.time()
sol = dae_solver.solve(model, t_eval);
end_simtime = time.time()
print("Simulation Time : %f sec \n" %(end_simtime - start_simtime))

Plotting simulation variables

ttime = sol.t
def extract_simdata(tt,var):
vart = var(tt)
return vart

Positive Electrode related variables

cs_pos = extract_simdata(ttime,sol["csp"])
avgcs_pos = extract_simdata(ttime,sol["cspavg"])
phi_pos = extract_simdata(ttime,sol["phip"])

Negative Electrode related variables

cs_neg = extract_simdata(ttime,sol["csn"])
avgcs_neg = extract_simdata(ttime,sol["csnavg"])
phi_neg = extract_simdata(ttime,sol["phin"])

Cell Voltage

cellpot = extract_simdata(ttime,sol["v"])

Charge Withdrawn / Charge Stored

Q_cell = extract_simdata(ttime,sol["Q"])

Iapp = extract_simdata(ttime,sol["iapp"])

plt.figure(figsize=(12, 9))
plt.subplot(2,2,1)
plt.plot(ttime,cs_pos,color='#66CC00',label="Positive Electrode")
plt.plot(ttime,cs_neg,color='#660033',label="Negative Electrode")
plt.xlabel("Time (s)", fontsize=12) #, fontsize=16
plt.ylabel("Surface Concentration (scaled)", fontsize=12) #, fontsize=16
plt.legend(loc="center right")
plt.grid(axis='both')

plt.subplot(2,2,2)
plt.plot(ttime,avgcs_pos,color='#66CC00',label="Positive Electrode")
plt.plot(ttime,avgcs_neg,color='#660033',label="Negative Electrode")
plt.xlabel("Time (s)", fontsize=12) #, fontsize=16
plt.ylabel("Avg Surface Concentration (scaled)", fontsize=12) #, fontsize=16
plt.legend(loc="center right")
plt.grid(axis='both')

plt.subplot(2,2,3)
plt.plot(ttime,cellpot,color='#660033')
plt.xlabel("Time (s)", fontsize=12) #, fontsize=16
plt.ylabel("Cell Potential (V)", fontsize=12) #, fontsize=16
plt.grid(axis='both')

plt.subplot(2,2,4)
plt.plot(ttime,Iapp,color='#660033')
plt.xlabel("Time (s)", fontsize=12) #, fontsize=16
plt.ylabel("Current (A/m2)", fontsize=12) #, fontsize=16
plt.grid(axis='both')
`

[valentinsulzer]

You can't do experimental protocols with custom models. But SPM({"particle": "quadratic profile"}) should be a similar model to your one

[YannickNoelStephanKuhn]

Technically, you can make your equations work with experimental protocols, but the way to do it would make your code incompatible with future releases of PyBaMM.
First, you would need to make your model extend pybamm.BaseSubModel, let's call that YourModel. That involves adding a lot of (superfluous to you) variable names.
Second, you extend pybamm.BaseBatteryModel, much like https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/models/full_battery_models/lithium_ion/spm.py. You only need self.set_external_circuit_submodel() to enable the compatibility with experimental protocols. Your equations then go into self.submodels[<your_model_name>] = YourModel(...).
If tinosulzer's model suggestion works well enough, I would not bother with my suggestion.

[valentinsulzer]

that's a good idea

[suryankolluri]

Hello Valentin Sulzer, Understood!! Thanks for the response. Surya e

…

Sent from my iPhone

On Jun 27, 2022, at 11:08 PM, Valentin Sulzer ***@***.***> wrote:  You can't do experimental protocols with custom models. But SPM({"particle": "quadratic profile"}) should be a similar model to your one — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.