In [ ]:
%load_ext autoreload
%autoreload 2
In [2]:
import os
from pathlib import Path
from pprint import pprint
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cellpy
from cellpy import prms
from cellpy import prmreader
from cellpy import cellreader
from cellpy.utils import ocv_rlx
import holoviews as hv
%matplotlib inline
hv.extension('bokeh')
In [3]:
######################################################################
## ##
## development ##
## ##
######################################################################
if os.name == 'nt':
# Use these when working on my work PC:
raw_data_path = r"C:\Scripting\MyFiles\development_cellpy\dev_data\gitt"
out_data_path = r"C:\Scripting\MyFiles\development_cellpy\out"
else:
# Use these when working on my MacBook:
raw_data_path = "/Users/jepe/scripting/cellpy/dev_data/gitt"
out_data_path = "/Users/jepe/scripting/cellpy/dev_data/out"
raw_data_path = Path(raw_data_path)
out_data_path = Path(out_data_path)
print(" SETTING SOME PRMS ".center(80, "="))
prms.Paths["db_filename"] = "cellpy_db.xlsx"
prms.Paths["cellpydatadir"] = out_data_path
prms.Paths["outdatadir"] = out_data_path
prms.Paths["rawdatadir"] = raw_data_path
prms.Paths["db_path"] = out_data_path
prms.Paths["filelogdir"] = out_data_path
pprint(prms.Paths)
pd.set_option('display.max_rows', 30)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)
In [4]:
fn = "20190403_cen59_04_rateGITT_01"
mass = 0.3 # mg
resn = fn + ".res"
cellpyn = fn + ".h5"
filename = prms.Paths["rawdatadir"] / resn
cellpyname = prms.Paths["cellpydatadir"] / cellpyn
In [5]:
cell = cellreader.get(filename, mass=mass, logging_mode="INFO")
cell.save(cellpyname)
cell = cellreader.get(cellpyname)
In [6]:
# print(cell)
In [7]:
steptable = cell.make_step_table(all_steps=True).cell.steps
In [8]:
steptable[steptable.cycle==5].head(4)
Out[8]:
In [9]:
# retrieving the voltage vs time curve
t, v = cell.get_timestamp(), cell.get_voltage()
all_cycs = hv.Curve((t,v), ('t', 'time (sec)'), ('v', 'voltage (v)'), label="voltage-time").opts()
In [10]:
# creating labels
cycle_label_df = steptable.drop_duplicates("cycle")
cycle_labels = hv.Labels((cycle_label_df.test_time_first, cycle_label_df.voltage_first, cycle_label_df.cycle))
In [11]:
# creating a dictionary of curves (for each cycle)
cycs_dict = dict()
for c in cell.get_cycle_numbers():
t = cell.get_timestamp(cycle=c)
t = t - np.amin(t.values) # setting first point to t=0
curve = hv.Curve((t, cell.get_voltage(cycle=c)), ("time", "time (seconds)"), ("voltage", "voltage (v vs. Li/Li+)"))
cycs_dict[c] = curve
In [12]:
# creating a holomap object
hmap = hv.HoloMap(cycs_dict,"cycle")
In [13]:
%%opts Curve [width=800, xformatter="%6.0f"]
(all_cycs * cycle_labels + hmap).cols(1)
Out[13]:
In [14]:
from bokeh.models import HoverTool
hover = HoverTool(tooltips=[("index", "$index")])
As we have seen, only some of the cycles contains the GITT curves. We need to select one of those. Lets start with cycle 5
In [15]:
dfdata = cell.cell.raw
In [16]:
cyc5 = dfdata.loc[dfdata.Cycle_Index == 5, :]
gt5 = steptable.loc[steptable.cycle==5, :]
In [17]:
cyc5curve = hv.Curve(cyc5, ("Test_Time", "time"), ("Voltage", "voltage")).opts(color="grey", alpha=0.5)
cyc5points = hv.Scatter(cyc5, ("Test_Time", "time"), ("Voltage", "voltage")).opts(size=5, fill_alpha=0.2)
In [18]:
slabels5 = hv.Labels((gt5.test_time_first, gt5.voltage_first, gt5.ustep))
In [19]:
%%opts Curve [width=1000, height=600, tools=['hover']]
cyc5curve * cyc5points * slabels5
Out[19]:
In [20]:
gt5.head()
Out[20]:
In [21]:
# approximating electrode-electrolyte contact area
diameter_si = 200.0 * 10e-7 # cm
rho_si = 2.32 # g/cm3
_mass = mass/1000
area = (2*3/diameter_si) * (_mass/rho_si) # CHECK THIS
print(f"Diameter: {diameter_si:8.6f} cm")
print(f"Mass: {mass:8.6f} mg")
print(f"Calculated contact area (ideal case): {area:8.6f} cm2")
# For reference (v, a, m is pr Si particle)
v = (4/3)*np.pi*((diameter_si/2)**3)
m = rho_si * v
n = mass / m
a = 4 * np.pi * (diameter_si/2)**2
In [22]:
# estimate some values
# electrolyte-electrode area (cm2)
radius = 0.75 # cm
roughness = 1.0 # a dimensionless parameter ranging from 1 to a lot.
print("Should probably use BET to find out this?")
area = (1*roughness) * np.pi * (radius**2)
print(f"S: {area:6.2f} cm2")
print("You should compare this number to what was calculated in the cell above (the theoretical surface area for the set of same sized nanoparticles)")
# number of moles
# 1 mol Si weighs 28.0855 g
_mass = mass / 1000 # convert mass from mg to g
number_of_moles = _mass / 28.0855
print(f"n_m: {number_of_moles:8.6f} mol")
# molar volume
# The 2006 CODATA recommended value for the molar volume of silicon is 12.0588349(11)×10−6 m3/mol, with a relative standard uncertainty of 9.1×10−8
molar_volume = 12.06 # cm3/mol
print(f"V_m: {molar_volume:6.2f} cm3/mol")
In [ ]:
In [ ]:
The following equation to estimate the Diffusion Coefficients appears in a paper:
\begin{equation} D = \frac{4}{\pi\cdot\tau} \cdot ( \frac {m \cdot V_m}{M \cdot S})^2 \cdot ( \frac {\Delta E_s}{ \Delta E_t})^2 \: for \: \tau << L^2 / D \end{equation}$ \Delta E_s $ = change of steady state voltage at the end of two sequential open-circuit relaxation periods
$ \Delta E_t $ = total change in the cell voltage $ E $ during the current pulse, neglecting the IR drop
$ V_s $ = mole volume of active material (cm^3/mol) (molar mass divided by the materials density)
$ M $ = ?
$ S $ = Surface area of the electrode (cm^2)
$ \tau $ = time of the current pulse
$ L $ = diffustion length
Andrzej P. Nowak et al., Procedia Engineering 98 (2014) 8
In [23]:
# Simple calculations by reading values from the graph
m = mass/1000
molar_mass = 28.0855 # g/mol
density = 2.32 # g/cm3
Vm = molar_mass / density # cm3/mol
S = 1.7 # cm2
t = 2880 # s
M = 1
Es1 = 0.378
Es2 = 0.344
DEs = Es2-Es1
Et11 = 0.378 # starting the current pulse
Et12 = 0.309 # after (possible) IR drop
Et2 = 0.28
DEt1 = Et2 - Et11
DEt2 = Et2 - Et12
p1 = (4 / (np.pi * t)) * (m * Vm / (M * S))**2
p21 = (DEs / DEt1)**2
p22 = (DEs / DEt2)**2
D1 = p1*p21
D2 = p1*p22
print(f"Assuming no IR drop: {D1}\nWith IR drop: {D2}")
print(f"A: {p1}")
print(f"DEt: {DEt1}")
print(f"DEs: {DEs}")
Comparing another reference (below) I assume that M is 1/number_of_moles ?
In [ ]:
This function is based on the Metrohm Autolab document. It will need some adjustments.
In [24]:
def calc_A(n_m=0.0011, V_m=12.06, S=1.77):
"""
D = 4 /((pi)(tau)) * (n_m * V_m / S)^2 * (DEs / DEt)^2 = A /tau * (DEs / DEt)^2
A = 4 /(pi) * (n_m * V_m / S)^2
tau: duration of the current pulse (s)
n_m: number of moles (mol)
V_m: molar volume of electrode (cm3/mol)
S: electrode-electrolyte contact area (cm2)
DEs: steady state voltage change due to the current pulse (V)
DEt: voltage change during the constant current pulse (eliminating the iR drop) (V)
Ref.: application note from Metrohm Autolab b.v. pdf (BAT03)
"""
A = (4 / np.pi) * (n_m * V_m / S)**2
return A
In [25]:
def auto_calc_D(steptable, cycle_number, A=1.0, tau=None, ustep_first=None, ustep_last=None):
"""Function for extracting diffusion constant(s) and inserting into steptable.
D = 4 /((pi)(tau)) * (n_m * V_m / S)^2 * (DEs / DEt)^2 = A /tau * (DEs / DEt)^2
A = 4 /(pi) * (n_m * V_m / S)^2
OBS! not corrected for IR drop yet.
"""
st = None
if cycle_number is None:
# This function is intended to only work on a pr. cycle basis
# to prevent users from "poluting" the steptable for "non_GITT" experiments.
print("no cycle number given")
return
st = steptable[steptable.cycle==cycle_number]
st = st[st.type.isin(["charge", "discharge", "ocvrlx_up", "ocvrlx_down"])]
if st.empty:
print("the given cycle is not found")
return
if ustep_first is not None:
st = st[st.ustep>=ustep_first]
if ustep_last is not None:
st = st[st.ustep<=ustep_last]
# used for finding DE
n3 = st["voltage_last"].shift(periods=-3)
n2 = st["voltage_last"].shift(periods=-2)
n1 = st["voltage_last"].shift(periods=-1)
n0 = st["voltage_last"]
# st["n3"] = n3
# st["n2"] = st["voltage_last"].shift(periods=-2)
# st["n1"] = n1
# st["n0"] = st["voltage_last"].shift(periods=0)
# used for finding tau
if tau is None:
tau = st["step_time_last"] - st["step_time_first"]
# used for validating if proper GITT step
t0 = st["type"]
t1 = st["type"].shift(periods=-1)
t2 = st["type"].shift(periods=-2)
t3 = st["type"].shift(periods=-3)
st["valid_D"] = (t0==t2) & (t1.str.contains("ocv"))
# calculating
st["DEt"] = st["voltage_last"] - st["voltage_first"] # should remove the IR drop here (maybe calculate it based on the measured IR time the current?)
st["DEs"] = n3 - n1
st["D"] = A / tau * (st["DEs"] / st["DEt"])**2
return st
In [26]:
# Lets just run it and see what happens
A = calc_A(n_m=0.0011, V_m=12.06, S=1.77)
print(f"A: {A}")
x = auto_calc_D(steptable, 5, A=A, ustep_first=27)
discharge = x.loc[(x.type=="discharge") & (x.valid_D)]
charge = x.loc[(x.type=="charge") & (x.valid_D)]
# discharge.head(25)
In [27]:
%%opts Curve [width=600 xformatter="%6.0f" tools=["hover"]]
X = discharge["discharge_avr"]*1_000_000/mass
Y = discharge["D"]
discharge_diffcurve = hv.Scatter((X,Y), ("capacity"), ("diffusion coef. (cm2 s-1)"), label="discharge").opts(size=12) * hv.Curve((X,Y)).opts(alpha=0.4)
X = charge["charge_avr"]*1_000_000/mass
Y = charge["D"]
charge_diffcurve = hv.Scatter((X,Y), ("capacity"), ("diffusion coef. (cm2 s-1)"), label="charge").opts(size=12) * hv.Curve((X,Y)).opts(alpha=0.4)
discharge_diffcurve * charge_diffcurve
Out[27]:
In [ ]:
In [28]:
cycle = 5
titcurve_ustep = 28
# filtering wrt cycle number to get the Data_Points for the step
datapoints = x.loc[(x.type=="discharge") & (x.ustep==titcurve_ustep), ["point_first", "point_last"]]
datapoints
Out[28]:
In [29]:
# Using the data points to select the data from dfdata
first = datapoints.iloc[0, 0]
last = datapoints.iloc[0, 1]
dftit = dfdata.loc[first:last, :]
dftit.head()
Out[29]:
In [30]:
dfcurve = hv.Curve(dftit, ("Step_Time", "time"), ("Voltage", "voltage")).opts(color="grey", alpha=0.5).opts(width=1000, xformatter="%6.0f")
min_time = hv.VLine(100.0)
max_time = hv.VLine(1000.0)
a = hv.Arrow(x=500, y = 0.31, text="fit", direction= "v")
dfcurve * min_time * max_time * a
Out[30]:
In [31]:
# plotting dE vs sqrt(t) and doing a linear regression to find the slope
dE = dftit["Voltage"]
t = np.sqrt(dftit["Step_Time"])
In [32]:
decurve = hv.Scatter((t, dE), "sqrt-time", "voltage").opts(width=1000, xformatter="%6.0f")
decurve
Out[32]:
In [ ]:
Eq. 1.1.
D = (4 / pi) (i*V_m / (Z_A*F*S)^2 * ((dE/dd)/(dE/d(sqrt(t)))^2
i: current (A)
Z_A: charge number
F: Faraday´s constant (96_458 C/mol)
V_m: molar volume of electrode (cm3/mol)
S: electrode-electrolyte contact area (cm2)
(dE/dd): the slope of the coulometric titration curve, found by plotting the
steady state voltages E (V) measured after each titration step δ
steady state voltage change due to the current pulse (V)
(dE/d(sqrt(t)): the slope of the linearized plot of the potential E (V) during the current pulse of duration t (s).Need to find the slope of E vs δ. This is the same as finding dV/dQ. One way to do this is to fit the "relaxed" voltage curve (the OCV_RLX points) to a polynomial and diffing it.
Then we also need to find the slope of E vs time during the pulse. For this we can e.g. use the polytfit function from scipy. The first 100 seconds should not be used as they are typically from charge-transfere resistance and the concentration overpotential.
Now it is time to get all the other prms (where the most difficult one is the electrode-electrolyte contact area). It should ideally be found using BET. But I guess the best is to calculate it using a sensible model, or from some scattering experiments (SANS?).
In [ ]:
Based on application note from Metrohm Autolab b.v. pdf (BAT03)
Eq. 1.1.
D = (4 / pi) (i*V_m / (Z_A*F*S)^2 * ((dE/dd)/(dE/d(sqrt(t)))^2
i: current (A)
Z_A: charge number
F: Faraday´s constant (96_458 C/mol)
V_m: molar volume of electrode (cm3/mol)
S: electrode-electrolyte contact area (cm2)
(dE/dd): the slope of the coulometric titration curve, found by plotting the
steady state voltages E (V) measured after each titration step δ
steady state voltage change due to the current pulse (V)
(dE/d(sqrt(t)): the slope of the linearized plot of the potential E (V) during the current pulse of duration t (s).Eq. 1.2.
D = 4 /((pi)(tau)) * (n_m * V_m / S)^2 * (DEs / DEt)^2
tau: duration of the current pulse (s)
n_m: number of moles (mol)
V_m: molar volume of electrode (cm3/mol)
S: electrode-electrolyte contact area (cm2)
DEs: steady state voltage change due to the current pulse (V)
DEt: voltage change during the constant current pulse (eliminating the iR drop) (V)
In [ ]:
D = (A / tau) * (DEs / DEt)**2
# during charge
DEt = ustep_52.voltage[-1] - ustep_52.voltage[0] # charge step
DEs = ustep_55.voltage[-1] - ustep_53.voltage[0] # ocv steps
# during discharge
DEt = ustep_27.voltage[-1] - ustep_27.voltage[0] # discharge step
DEs = ustep_30.voltage[-1] - ustep_28.voltage[0] # ocv steps
Using step-table
DEt = ustep[n].voltage_last - ustep[n].voltage_first # charge step
DEs = ustep[n+3].voltage_last - ustep[n+1].voltage_first # ocv steps