research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Improved lithium-ion battery cti file and Matlab example

Browse source 
Extended and clarified comments, corrected density data, improved functionality
This commit is contained in:
wbessler 2019-05-06 10:44:00 +02:00 • committed by Ray Speth
parent ae792dde00
commit bc6dd4ddc0
2 changed files with 200 additions and 132 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

279
data/inputs/lithium_ion_battery.cti
View file

@ -1,21 +1,30 @@
#===================================================================================== #==============================================================================
# Cantera input file for an LCO/graphite lithium-ion battery # Cantera input file for an LCO/graphite lithium-ion battery
#
# This file includes a full set of thermodynamic and kinetic parameters of a
# lithium-ion battery, in particular:
# - Active materials: LiCoO2 (LCO) and LiC6 (graphite)
# - Organic electrolyte: EC/PC with 1M LiPF6
# - Interfaces: LCO/electrolyte and LiC6/electrolyte
# - Charge-transfer reactions at the two interfaces
#
# A MATLAB example using this file for simulating a discharge curve is
# samples/matlab/lithium_ion_battery.m
#
# Reference: # Reference:
# M. Mayur, S. DeCaluwe, B. L. Kee, W. G. Bessler, "Modeling # M. Mayur, S. DeCaluwe, B. L. Kee, W. G. Bessler, "Modeling
# thermodynamics and kinetics of intercalation phases for lithium-ion # thermodynamics and kinetics of intercalation phases for lithium-ion
# batteries in Cantera", Computer Physics Communications # batteries in Cantera", under review at Electrochimica Acta.
#===================================================================================== #==============================================================================
#==============================================================================
#=====================================================================================
# Bulk phases # Bulk phases
#===================================================================================== #==============================================================================
#------------------------------------------------------------------------------
#------------------------------------------------------------------
# Graphite (anode) # Graphite (anode)
# Thermodynamic data based on half-cell measurements by K. Kumaresan et al., J. Electrochem. Soc. 155, A164-A171 (2008) # Thermodynamic data based on half-cell measurements by K. Kumaresan et al.,
# Density: 5031.67 kg/m3 - used to calculate species molar volume as molecular weight (MW)/density # J. Electrochem. Soc. 155, A164-A171 (2008)
#------------------------------------------------------------------ #------------------------------------------------------------------------------
BinarySolutionTabulatedThermo( BinarySolutionTabulatedThermo(
name = "anode", name = "anode",
elements = "Li C", elements = "Li C",
@ -48,26 +57,11 @@ BinarySolutionTabulatedThermo(
1.92885E+01, 1.92876E+01, 1.92837E+01, 1.92769E+01, 1.92850E+01, 1.93100E+01, 1.93514E+01], 1.92885E+01, 1.92876E+01, 1.92837E+01, 1.92769E+01, 1.92850E+01, 1.93100E+01, 1.93514E+01],
"J/mol/K"))) "J/mol/K")))
# Lithium intercalated in graphite, MW: 79.0070 g/mol #------------------------------------------------------------------------------
species(
name = "Li[anode]",
atoms = "Li:1 C:6",
thermo = const_cp(h0 = (0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')), # these are dummy entries because the values are taken from the table
standardState = constantIncompressible(molarVolume = (79.0070/5.0317, 'cm3/gmol')))
# Vacancy in graphite, MW: 72.0660 g/mol. Note this species includes the carbon host matrix.
species(
name = "V[anode]",
atoms = "C:6",
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')), # values are set to 0 because this is the reference species for this phase
standardState = constantIncompressible(molarVolume = (72.0660/5.0317, 'cm3/gmol')))
#------------------------------------------------------------------
# Lithium cobalt oxide (cathode) # Lithium cobalt oxide (cathode)
# Thermodynamic data based on half-cell measurements by K. Kumaresan et al., J. Electrochem. Soc. 155, A164-A171 (2008) # Thermodynamic data based on half-cell measurements by K. Kumaresan et al.,
# Density: 2292 kg/m3 - used to calculate species molar volume as molecular weight (MW)/density # J. Electrochem. Soc. 155, A164-A171 (2008)
#------------------------------------------------------------------ #------------------------------------------------------------------------------
BinarySolutionTabulatedThermo( BinarySolutionTabulatedThermo(
name = "cathode", name = "cathode",
elements = "Li Co O", elements = "Li Co O",
@ -100,24 +94,21 @@ BinarySolutionTabulatedThermo(
-5.39817E+01, -5.45468E+01, -5.48125E+01, -5.51520E+01, -5.54526E+01, -5.52961E+01, -5.50219E+01, -5.46653E+01, -5.42305E+01], -5.39817E+01, -5.45468E+01, -5.48125E+01, -5.51520E+01, -5.54526E+01, -5.52961E+01, -5.50219E+01, -5.46653E+01, -5.42305E+01],
"J/mol/K"))) "J/mol/K")))
# Lithium cobalt oxide, MW: 97.8730 g/mol #------------------------------------------------------------------------------
species( # Carbonate based electrolyte
name = "Li[cathode]", # Solvent: Ethylene carbonate:Propylene carbonate (1:1 v/v)
atoms = "Li:1 Co:1 O:2", # Salt: 1M LiPF6
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')), # these are dummy entries because the values are taken from the table #------------------------------------------------------------------------------
standardState = constantIncompressible(molarVolume = (97.8730/2.292, 'cm3/gmol'))) IdealSolidSolution(
name = "electrolyte",
elements = "Li P F C H O E",
species = "C3H4O3[elyt] C4H6O3[elyt] Li+[elyt] PF6-[elyt]",
initial_state = state(mole_fractions = 'C3H4O3[elyt]:0.47901 C4H6O3[elyt]:0.37563 Li+[elyt]:0.07268 PF6-[elyt]:0.07268'),
standard_concentration = "unity")
# Vacancy in the cobalt oxide, MW: 90.9320 g/mol. Note this species includes the host matrix. #------------------------------------------------------------------------------
species(
name = "V[cathode]",
atoms = "Co:1 O:2",
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')), # values are set to 0 because this is the reference species for this phase
standardState = constantIncompressible(molarVolume = (90.9320/2.292, 'cm3/gmol')))
#------------------------------------------------------------------
# Electron conductor # Electron conductor
#------------------------------------------------------------------ #------------------------------------------------------------------------------
metal( metal(
name = "electron", name = "electron",
elements = "E", elements = "E",
@ -125,90 +116,172 @@ metal(
density = (1.0, 'kg/m3'), # dummy entry density = (1.0, 'kg/m3'), # dummy entry
initial_state = state( mole_fractions = "electron:1.0")) initial_state = state( mole_fractions = "electron:1.0"))
# Electron, MW: 0.000545 g/mol
#==============================================================================
# Species
#==============================================================================
#------------------------------------------------------------------------------
# Lithium intercalated in graphite, MW: 79.0070 g/mol.
# Note this species includes the carbon host matrix.
# Molar enthalpy and entropy are set to 0 because the values given in the
# BinarySolidSolutionTabulatedThermo class are used.
# Density of graphite: 2270 kg/m3 (W. M. Haynes et al, CRC Handbook of Chemistry
# and Physics, 94th edition, CRC press, Boca Raton, London, New York, 2013)
# (used to calculate species molar volume as molecular weight (MW)/density).
#------------------------------------------------------------------------------
species( species(
name = "electron", name = "Li[anode]",
atoms = "E:1", atoms = "Li:1 C:6",
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K'))) # dummy entries because chemical potential is set to zero for a "metal" phase thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (79.0070/2.270, 'cm3/mol')))
#------------------------------------------------------------------------------
# Vacancy in graphite, MW: 72.0660 g/mol.
# Note this species includes the carbon host matrix.
# Molar enthalpy and entropy are set to 0 because this is the reference species
# for this phase.
# Density of graphite: 2270 kg/m3 (W. M. Haynes et al, CRC Handbook of Chemistry
# and Physics, 94th edition, CRC press, Boca Raton, London, New York, 2013)
# (used to calculate species molar volume as molecular weight (MW)/density).
#------------------------------------------------------------------------------
species(
name = "V[anode]",
atoms = "C:6",
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (72.0660/2.270, 'cm3/mol')))
#------------------------------------------------------------------------------
# Lithium cobalt oxide, MW: 97.8730 g/mol.
# Note this species includes the cobalt oxide host matrix.
# Molar enthalpy and entropy are set to 0 because the values given in the
# BinarySolidSolutionTabulatedThermo class are used.
# Density of LCO: 4790 kg/m3 (E.J. Cheng et al., J. Asian Ceramic Soc. 5, 113,
# 2017) (used to calculate species molar volume as molecular weight/density).
#------------------------------------------------------------------------------
species(
name = "Li[cathode]",
atoms = "Li:1 Co:1 O:2",
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (97.8730/4.790, 'cm3/mol')))
#-------------------------------------------------------------------- #------------------------------------------------------------------------------
# Carbonate based electrolyte # Vacancy in the cobalt oxide, MW: 90.9320 g/mol.
# Solvent: Ethylene carbonate:Propylene carbonate (1:1 v/v) # Note this species includes the cobalt oxide host matrix.
# Salt: 1M LiPF6 # Molar enthalpy and entropy are set to 0 because this is the reference species
# Density: 1260.0 kg/m3 - used to calculate species molar volume as molecular weight (MW)/density # for this phase.
#-------------------------------------------------------------------- # Density of LCO: 4790 kg/m3 (E.J. Cheng et al., J. Asian Ceramic Soc. 5, 113,
IdealSolidSolution( # 2017) (used to calculate species molar volume as molecular weight/density).
name = "electrolyte", #------------------------------------------------------------------------------
elements = "Li P F C H O E", species(
species = "C3H4O3[elyt] C4H6O3[elyt] Li+[elyt] PF6-[elyt]", name = "V[cathode]",
initial_state = state(pressure = OneAtm, mole_fractions = 'C3H4O3[elyt]:0.47901 C4H6O3[elyt]:0.37563 Li+[elyt]:0.07268 PF6-[elyt]:0.07268'), atoms = "Co:1 O:2",
standard_concentration = "unity") thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (90.9320/4.790, 'cm3/mol')))
#------------------------------------------------------------------------------
# Ethylene carbonate, MW: 88.0630 g/mol # Ethylene carbonate, MW: 88.0630 g/mol
# Density of electrolyte: 1260 kg/m3 (used to calculate species molar volume
# as molecular weight (MW)/density)
# Molar enthalpy and entropy set to zero (dummy entries as this species does
# not participate in chemical reactions)
#------------------------------------------------------------------------------
species( species(
name = "C3H4O3[elyt]", name = "C3H4O3[elyt]",
atoms = "C:3 H:4 O:3", atoms = "C:3 H:4 O:3",
thermo = const_cp(h0 =(0.0, 'J/mol'), s0 = (0.0, 'J/mol/K')), # Dummy entries as this species does not participate in chemical reactions thermo = const_cp(h0 =(0.0, 'J/mol'), s0 = (0.0, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (88.0630/1.260, 'cm3/gmol'))) standardState = constantIncompressible(molarVolume = (88.0630/1.260, 'cm3/mol')))
#------------------------------------------------------------------------------
# Propylene carbonate, MW: 102.0898 g/mol # Propylene carbonate, MW: 102.0898 g/mol
# Density of electrolyte: 1260.0 kg/m3 (used to calculate species molar volume
# as molecular weight (MW)/density)
# Molar enthalpy and entropy set to zero (dummy entries as this species does
# not participate in chemical reactions)
#------------------------------------------------------------------------------
species( species(
name = "C4H6O3[elyt]", name = "C4H6O3[elyt]",
atoms = "C:4 H:6 O:3", atoms = "C:4 H:6 O:3",
thermo = const_cp(h0 =(0.0, 'J/mol'), s0 = (0.0, 'J/mol/K')), # Dummy entries as this species does not participate in chemical reactions thermo = const_cp(h0 =(0.0, 'J/mol'), s0 = (0.0, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (102.0898/1.260, 'cm3/gmol'))) standardState = constantIncompressible(molarVolume = (102.0898/1.260, 'cm3/mol')))
#------------------------------------------------------------------------------
# Lithium ion, MW: 6.940455 g/mol # Lithium ion, MW: 6.940455 g/mol
# Density of electrolyte: 1260.0 kg/m3 (used to calculate species molar volume
# as molecular weight (MW)/density)
# Molar enthalpy and entropy taken from Li+(aq) from P. Atkins "Physical
# Chemistry", Wiley-VCH (2006)
#------------------------------------------------------------------------------
species( species(
name = "Li+[elyt]", name = "Li+[elyt]",
atoms = "Li:1 E:-1", atoms = "Li:1 E:-1",
thermo = const_cp(h0 = (-278.49, 'kJ/mol'), s0 = (13.4, 'J/mol/K')), # Li+(aq) from P. Atkins "Physical Chemistry", Wiley-VCH (2006) thermo = const_cp(h0 = (-278.49, 'kJ/mol'), s0 = (13.4, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (6.940455/1.260, 'cm3/gmol'))) standardState = constantIncompressible(molarVolume = (6.940455/1.260, 'cm3/mol')))
#------------------------------------------------------------------------------
# Hexafluorophosphate ion, MW: 144.964745 g/mol # Hexafluorophosphate ion, MW: 144.964745 g/mol
# Density of electrolyte: 1260.0 kg/m3 (used to calculate species molar volume
# as molecular weight (MW)/density)
# Molar enthalpy and entropy set to zero (dummy entries as this species does
# not participate in chemical reactions)
#------------------------------------------------------------------------------
species( species(
name = "PF6-[elyt]", name = "PF6-[elyt]",
atoms = "P:1 F:6 E:1", atoms = "P:1 F:6 E:1",
thermo = const_cp(h0 = (0.0, 'J/mol'), s0 = (0.0, 'J/mol/K')), # Dummy entries as this species does not participate in chemical reactions thermo = const_cp(h0 = (0.0, 'J/mol'), s0 = (0.0, 'J/mol/K')),
standardState = constantIncompressible(molarVolume = (144.964745/1.260, 'cm3/gmol'))) standardState = constantIncompressible(molarVolume = (144.964745/1.260, 'cm3/mol')))
#------------------------------------------------------------------------------
# Electron, MW: 0.000545 g/mol
# Molar enthalpy and entropy set to zero (dummy entries because chemical
# potential is set to zero for a "metal" phase)
#------------------------------------------------------------------------------
species(
name = "electron",
atoms = "E:1",
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')))
#------------------------------------------------------------------------------
#===================================================================================== # Dummy species (needed for defining the interfaces)
# Interfaces for electrochemical reactions #------------------------------------------------------------------------------
#=====================================================================================
#--------------------------------------------------------------------
# Anode reaction
#--------------------------------------------------------------------
ideal_interface(
name = "edge_anode_electrolyte",
phases = "anode electron electrolyte",
reactions = "anode_*",
elements = "Li E C",
species = "(dummy)", # dummy entry for global kinetics
site_density = (1.0e-2, 'mol/cm2')) # dummy entry for global kinetics
edge_reaction("Li[anode] <=> Li+[elyt] + V[anode] + electron", [4, 0.0, (0, 'kJ/mol')], rate_coeff_type = "exchangecurrentdensity", beta = 0.5,id="anode_reaction")
#--------------------------------------------------------------------
# Cathode reaction
#--------------------------------------------------------------------
ideal_interface(
name = "edge_cathode_electrolyte",
phases = "cathode electron electrolyte",
reactions = "cathode_*",
elements = "Li E Co O",
species = "(dummy)", # dummy entry for global kinetics
site_density = (1.0e-2, 'mol/cm2')) # dummy entry for global kinetics
edge_reaction("Li+[elyt] + V[cathode] + electron <=> Li[cathode]", [100, 0.0, (0, 'kJ/mol')], rate_coeff_type = "exchangecurrentdensity", beta = 0.5,id="cathode_reaction")
# Dummy species
species( species(
name = "(dummy)", name = "(dummy)",
atoms = "", atoms = "",
thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K'))) thermo = const_cp(h0 = (0.0, 'kJ/mol'), s0 = (0.0, 'J/mol/K')))
#==============================================================================
# Interfaces for electrochemical reactions
#==============================================================================
#------------------------------------------------------------------------------
# Graphite/electrolyte interface
# Species and site density are dummy entries (as we do not consider surface-
# adsorbed species)
#------------------------------------------------------------------------------
ideal_interface(
name = "edge_anode_electrolyte",
phases = "anode electron electrolyte",
reactions = "anode_*",
species = "(dummy)",
site_density = (1.0e-2, 'mol/cm2'))
#------------------------------------------------------------------------------
# LCO/electrolyte interface
# Species and site density are dummy entries (as we do not consider surface-
# adsorbed species)
#------------------------------------------------------------------------------
ideal_interface(
name = "edge_cathode_electrolyte",
phases = "cathode electron electrolyte",
reactions = "cathode_*",
species = "(dummy)",
site_density = (1.0e-2, 'mol/cm2'))
#==============================================================================
# Electrochemical reactions
#==============================================================================
# Use Butler-Volmer kinetics (rate_coeff_type = "exchangecurrentdensity").
edge_reaction("Li[anode] <=> Li+[elyt] + V[anode] + electron", [4, 0.0, (0, 'kJ/mol')], rate_coeff_type = "exchangecurrentdensity", beta = 0.5,id="anode_reaction")
edge_reaction("Li+[elyt] + V[cathode] + electron <=> Li[cathode]", [100, 0.0, (0, 'kJ/mol')], rate_coeff_type = "exchangecurrentdensity", beta = 0.5,id="cathode_reaction")

53
samples/matlab/lithium_ion_battery.m
View file

@ -15,13 +15,7 @@
% M. Mayur, S. DeCaluwe, B. L. Kee, W. G. Bessler, "Modeling % M. Mayur, S. DeCaluwe, B. L. Kee, W. G. Bessler, "Modeling
% thermodynamics and kinetics of intercalation phases for lithium-ion % thermodynamics and kinetics of intercalation phases for lithium-ion
% batteries in Cantera", under review at Electrochimica Acta. % batteries in Cantera", under review at Electrochimica Acta.
%
% For the sake of simplicity, we're going to assume that the anode and
% cathode capacities are perfectly balanced (i.e. if the cathode lithium
% content is X percent of it's max possible (i.e. its capacity), then we
% will assume that the anode is at 1-X percent. Without loss of
% generality, we will define the anode composition:
% The routine below returns the cell voltage (in Volt) of a lithium-ion % The routine below returns the cell voltage (in Volt) of a lithium-ion
% cell for a given cell current and active material lithium stoichiometries. % cell for a given cell current and active material lithium stoichiometries.
% %
@ -32,18 +26,16 @@
% - externally-applied current I_app [A] % - externally-applied current I_app [A]
% - electrolyte resistance R_elyt [Ohm] % - electrolyte resistance R_elyt [Ohm]
X_Li_an = [0.005:0.025:0.995];
X_Li_ca = 1 - X_Li_an;
I_app = 0; % Input parameters
R_elyt = 0; SOC = 0:0.02:1; % [-] Input state of charge (0...1)
T = 300; X_Li_an = (0.75-0.01)*SOC+0.01; % anode balancing
P = oneatm; X_Li_ca = (0.99-0.49)*(1-SOC)+0.49; % cathode balancing
I_app = 0; % [A] Externally-applied current
global F R_elyt = 0; % [Ohm] Electrolyte resistance
% Parameters T = 300; % [K] Temperature
P = oneatm; % [Pa] Pressure
inputCTI = 'lithium_ion_battery.cti'; % cantera input file name inputCTI = 'lithium_ion_battery.cti'; % cantera input file name
F = 96485; % Faraday's constant [C/mol]
S_ca = 1.1167; % [m^2] Cathode total active material surface area S_ca = 1.1167; % [m^2] Cathode total active material surface area
S_an = 0.7824; % [m^2] Anode total active material surface area S_an = 0.7824; % [m^2] Anode total active material surface area
@ -56,14 +48,16 @@ anode_interface = Interface(inputCTI, 'edge_anode_electrolyte', anode, elde, ely
cathode_interface = Interface(inputCTI, 'edge_cathode_electrolyte', cathode, elde, elyt); cathode_interface = Interface(inputCTI, 'edge_cathode_electrolyte', cathode, elde, elyt);
% Set the temperatures and pressures of all phases % Set the temperatures and pressures of all phases
phases = [anode elde elyt cathode]; set(anode,'T',T,'P',P);
for ph = phases set(cathode,'T',T,'P',P);
set(ph,'T',T,'P',P); set(elde,'T',T,'P',P);
end set(elyt,'T',T,'P',P);
set(anode_interface,'T',T,'P',P);
set(cathode_interface,'T',T,'P',P);
% Calculate cell voltage, separately for each entry of the input vectors % Calculate cell voltage, separately for each entry of the input vectors
E_cell = zeros(length(X_Li_ca),1); E_cell = zeros(length(SOC),1);
for i = 1:length(X_Li_ca) for i = 1:length(SOC)
% Set anode electrode potential to 0 % Set anode electrode potential to 0
phi_s_an = 0; phi_s_an = 0;
@ -80,11 +74,12 @@ for i = 1:length(X_Li_ca)
E_cell(i) = phi_s_ca - phi_s_an; E_cell(i) = phi_s_ca - phi_s_an;
end end
% Let's plot the cell voltage, as a function of the cathode stoichiometry: % Let's plot the cell voltage, as a function of the state of charge:
plot(X_Li_ca,E_cell,'linewidth',2.5) figure(1);
plot(SOC*100,E_cell,'linewidth',2.5)
ylim([2.5,4.3]) ylim([2.5,4.3])
xlabel('Li Fraction in Cathode') xlabel('State of charge / %')
ylabel('Cell potential [V]') ylabel('Cell voltage / V')
set(gca,'fontsize',14) set(gca,'fontsize',14)
@ -107,7 +102,7 @@ function anCurr = anode_curr(phi_s,phi_l,X_Li_an,anode,elde,elyt,anode_interface
r = rop_net(anode_interface).*1e3; % [mol/m2/s] r = rop_net(anode_interface).*1e3; % [mol/m2/s]
% Calculate the current % Calculate the current
anCurr = r*F*S_an*1; anCurr = r*96485*S_an*1; % F = 96485 C/mol Faraday's constant
end end
% This function returns the Cantera calculated cathode current (in A) % This function returns the Cantera calculated cathode current (in A)
@ -126,5 +121,5 @@ function caCurr = cathode_curr(phi_s,phi_l,X_Li_ca,cathode,elde,elyt,cathode_int
r = rop_net(cathode_interface).*1e3; % [mol/m2/s] r = rop_net(cathode_interface).*1e3; % [mol/m2/s]
% Calculate the current % Calculate the current
caCurr = r*F*S_ca*(-1); caCurr = r*96485*S_ca*(-1); % F = 96485 C/mol Faraday's constant
end end