diff --git a/bat_can.py b/bat_can.py index 4511fb1..063e1eb 100644 --- a/bat_can.py +++ b/bat_can.py @@ -162,7 +162,6 @@ def model_run(sim): for sim in parameters['simulations']: model = importlib.import_module('.'+sim['type'], package='simulations') - solution = model.plot(an, sep, ca, parameters, sim) #=========================================================================== diff --git a/bat_can_fit.py b/bat_can_fit.py index 3335536..8dbe842 100644 --- a/bat_can_fit.py +++ b/bat_can_fit.py @@ -8,7 +8,7 @@ """ # Import modules -from datetime import datetime +from datetime import datetime from functools import partial import importlib # allows us to import from user input string. import itertools @@ -23,7 +23,7 @@ from submodels.fitting import voltage_capacity as fit from bat_can_init import initialize -# This is the main function that runs the model. We define it this way so it +# This is the main function that runs the model. We define it this way so it # is called by "main," below: def bat_can(input, cores, print_flag): if input is None: @@ -55,32 +55,32 @@ def bat_can(input, cores, print_flag): #=========================================================================== # CREATE ELEMENT CLASSES AND INITIAL SOLUTION VECTOR SV_0 #=========================================================================== - # For each element (anode 'an', separator 'sep', cathode 'ca') the 'class' - # variable from the inputs tells what kind of anode, separator, or cathode - # it is, and points to a '.py' file in this directory. We import that - # module, and then run its 'initialize' routine to create an intial + # For each element (anode 'an', separator 'sep', cathode 'ca') the 'class' + # variable from the inputs tells what kind of anode, separator, or cathode + # it is, and points to a '.py' file in this directory. We import that + # module, and then run its 'initialize' routine to create an intial # solution vector and an object that stores needed parameters. # import single_particle_electrode as an_module_0 - - an_module = importlib.import_module('electrode_models.' + + an_module = importlib.import_module('electrode_models.' + an_inputs['class']) - an = an_module.electrode(input_file, an_inputs, sep_inputs, ca_inputs, + an = an_module.electrode(input_file, an_inputs, sep_inputs, ca_inputs, 'anode', parameters, offset=0) - - sep_module = importlib.import_module('separator_models.' + + sep_module = importlib.import_module('separator_models.' + sep_inputs['class']) - sep = sep_module.separator(input_file, sep_inputs, parameters, + sep = sep_module.separator(input_file, sep_inputs, parameters, offset=an.n_vars) - + # Check to see if the anode object needs to adjust the separator properties: sep = an.adjust_separator(sep) - ca_module = importlib.import_module('electrode_models.' + ca_module = importlib.import_module('electrode_models.' + ca_inputs['class']) - ca = ca_module.electrode(input_file, ca_inputs, sep_inputs, an_inputs, + ca = ca_module.electrode(input_file, ca_inputs, sep_inputs, an_inputs, 'cathode', parameters, offset= an.n_vars+sep.n_vars*sep.n_points) - # Check to see if the cathode object needs to adjust the separator - # properties: + # Check to see if the cathode object needs to adjust the separator + # properties: sep = ca.adjust_separator(sep) # Initialize the solution vector: @@ -96,23 +96,23 @@ def bat_can(input, cores, print_flag): #=========================================================================== # RUN THE SIMULATION #=========================================================================== - # The inputs tell us what type of experiment we will simulate. Load the + # The inputs tell us what type of experiment we will simulate. Load the # module, then call its 'run' function: # If the user requests a specific initailization routine, run that first: if 'initialize' in parameters and parameters['initialize']['enable']: if parameters['initialize']['type'] == 'open-circuit': - model = importlib.import_module('.'+'CC_cycle', + model = importlib.import_module('.'+'CC_cycle', package='simulations') - - sim = {'i_ext': None, 'C-rate': 0.0, 'n_cycles': 0, - 'first-step': 'discharge', 'equilibrate': - {'enable': True, 'time': 5}, 'phi-cutoff-lower': 2.0, + + sim = {'i_ext': None, 'C-rate': 0.0, 'n_cycles': 0, + 'first-step': 'discharge', 'equilibrate': + {'enable': True, 'time': 5}, 'phi-cutoff-lower': 2.0, 'phi-cutoff-upper': 4.8, 'init':True} - + solution = model.run(SV_0, an, sep, ca, algvars, parameters, sim) - - # Save final state as the initial state for all subsequent + + # Save final state as the initial state for all subsequent # simulation steps: SV_0 = model.final_state(solution) @@ -121,17 +121,17 @@ def bat_can(input, cores, print_flag): def run_model(x_guess, final_flag = False): print('x = ',x_guess) - + # Set the current fitting parameter values: for i, x in enumerate(x_guess): if fit_params[i]['type'] == 'elyte-transport': sep.D_k[sep.elyte_obj.species_index(fit_params[i]['species'])] \ = x elif fit_params[i]['type'] == 'cathode-kinetics': - ca.surf_obj.set_multiplier(x, + ca.surf_obj.set_multiplier(x, fit_params[i]['reaction-index']) elif fit_params[i]['type'] == 'air-kinetics': - ca.gas_elyte_obj.set_multiplier(x, + ca.gas_elyte_obj.set_multiplier(x, fit_params[i]['reaction-index']) elif fit_params[i]['type'] == 'cathode-microstructure': setattr(ca, fit_params[i]['parameter'], x) @@ -139,35 +139,35 @@ def run_model(x_guess, final_flag = False): if final_flag: SSR_net = 0 icolor = 0 - fit_fig, fit_axs = plt.subplots(1, 1, sharex=True, + fit_fig, fit_axs = plt.subplots(1, 1, sharex=True, gridspec_kw = {'wspace':0, 'hspace':0}) - + fit_fig.set_size_inches((5.0, 2.25)) # pool = mp.Pool(processes = int(cores)) # SSR_net = pool.map(run, (list(parameters['simulations']))) for sim in parameters['simulations']: - SSR_net += run(sim, SV_0, algvars, parameters, an, ca, sep, + SSR_net += run(sim, SV_0, algvars, parameters, an, ca, sep, final_flag, fit_fig, fit_axs, icolor) icolor += 1 else: with Pool(processes = int(cores)) as pool: - SSR_net = pool.map(partial(run, SV_0=SV_0, algvars=algvars, - parameters=parameters, an=an, sep=sep, ca=ca), + SSR_net = pool.map(partial(run, SV_0=SV_0, algvars=algvars, + parameters=parameters, an=an, sep=sep, ca=ca), list(parameters['simulations'])) print('SSR = ', np.sum(SSR_net)) - + if final_flag: fit_axs.annotate(f"SSR = 1.0", xy=(0,0))#, # xytext=(0.5, 0.5), textcoords='axes fraction')#, fontsize = 8) - + if len(parameters['simulations']) == 1: - savename = (parameters['output'] +'_' + savename = (parameters['output'] +'_' + sim['outputs']['save-name'] ) else: savename = (parameters['output']) - + fit_fig.savefig(savename + '/fit.pdf') plt.show() @@ -183,16 +183,16 @@ def run_model(x_guess, final_flag = False): for sim in parameters['simulations']: sim['ref_data'] = pd.read_excel('data/' + sim['validation']) - + if print_flag: - # Just print the fit of the starting guess - do not run the fitting + # Just print the fit of the starting guess - do not run the fitting # routine: result = run_model(x_start, final_flag=True) else: # Fit to the reference data: - result = spo.minimize(run_model, x_start, bounds = x_bounds, + result = spo.minimize(run_model, x_start, bounds = x_bounds, options={'disp': True}) - + print(result) print("Best fit = ", result.x) @@ -200,56 +200,56 @@ def run_model(x_guess, final_flag = False): run_model(result.x, final_flag=True) stop = timeit.default_timer() - print('Time: ', stop - start) + print('Time: ', stop - start) -def run(sim, SV_0, algvars, parameters, an, ca, sep, final_flag=False, +def run(sim, SV_0, algvars, parameters, an, ca, sep, final_flag=False, fit_fig=None, fit_axs=None, icolor=0): - + sim_local = sim try: - model = importlib.import_module('.'+sim['type'], + model = importlib.import_module('.'+sim['type'], package='simulations') sim['init'] = False - solution = model.run(SV_0, an, sep, ca, algvars, + solution = model.run(SV_0, an, sep, ca, algvars, parameters, sim) # Read out the results: if final_flag: sim_local['outputs']['show-plots'] = False - results = model.output(solution, an, sep, ca, parameters, + results = model.output(solution, an, sep, ca, parameters, sim_local, plot_flag=True, return_flag=True, save_flag=True) # for sim in parameters['simulations']: # model = importlib.import_module('.'+sim['type'], package='simulations') solution = model.plot(an, sep, ca, parameters, sim_local) else: - results = model.output(solution, an, sep, ca, parameters, + results = model.output(solution, an, sep, ca, parameters, sim, plot_flag=False, return_flag=True, save_flag=False) phi_sim = results['phi_ed'].to_numpy()[:,-1] - + sim_data = np.array((results['capacity'].to_numpy(), phi_sim)) # This function calculates the SSR for this simulation: # Scale to convert capacity units in validation data to mAh/cm2: - UnitsScale = 46.968 #todo #72 - ssr_calc = fit.SSR(sim['ref_data'].to_numpy(), sim_data.T, + UnitsScale = 1 #todo #72 46.968 + ssr_calc = fit.SSR(sim['ref_data'].to_numpy(), sim_data.T, units_scale = UnitsScale) print('SSR = ', ssr_calc) - + ndata = len(parameters['simulations']) + 1 cmap = plt.get_cmap('plasma') - + color_ind = np.linspace(0,1,ndata) colors = list() - + for i in np.arange(ndata): colors.append(cmap(color_ind[i])) - + if final_flag: - fit_axs, fit_fig = fit.plot(sim['ref_data'].to_numpy(), - sim_data.T, fit_axs, fit_fig, units_scale = UnitsScale, + fit_axs, fit_fig = fit.plot(sim['ref_data'].to_numpy(), + sim_data.T, fit_axs, fit_fig, units_scale = UnitsScale, color = colors[icolor]) - + except: # Assign a large penalty for failed parameter sets: ssr_calc = 1e23 @@ -265,12 +265,12 @@ def run(sim, SV_0, algvars, parameters, an, ca, sep, final_flag=False, if __name__ == '__main__': import argparse - # Currently, the only command line keyword enabled is --input, to specify + # Currently, the only command line keyword enabled is --input, to specify # the input file location: parser = argparse.ArgumentParser() parser.add_argument('--input') parser.add_argument('--cores') parser.add_argument('--print', action='store_true') args = parser.parse_args() - + bat_can(args.input, args.cores, args.print) diff --git a/data/Liu_O2/0.1mAhcm2.xlsx b/data/Liu_O2/0.1mAhcm2.xlsx new file mode 100644 index 0000000..a09fd47 Binary files /dev/null and b/data/Liu_O2/0.1mAhcm2.xlsx differ diff --git a/data/Liu_O2/0.2mAhcm2.xlsx b/data/Liu_O2/0.2mAhcm2.xlsx new file mode 100644 index 0000000..a50be85 Binary files /dev/null and b/data/Liu_O2/0.2mAhcm2.xlsx differ diff --git a/data/Liu_O2/0.5mAhcm2.xlsx b/data/Liu_O2/0.5mAhcm2.xlsx new file mode 100644 index 0000000..677e9ce Binary files /dev/null and b/data/Liu_O2/0.5mAhcm2.xlsx differ diff --git a/data/Liu_O2/1mAhcm2.xlsx b/data/Liu_O2/1mAhcm2.xlsx new file mode 100644 index 0000000..f71de3f Binary files /dev/null and b/data/Liu_O2/1mAhcm2.xlsx differ diff --git a/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index 5a239dd..c51eb0d 100644 --- a/electrode_models/air_electrode.py +++ b/electrode_models/air_electrode.py @@ -55,6 +55,8 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, # Phase volume fractions eps_host = np.array([inputs['eps_host']]) + eps_host = eps_host[0] + if len(eps_host) == 1: self.eps_host = np.repeat(eps_host[0], self.n_points) elif len(eps_host) == self.n_points: @@ -71,6 +73,8 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, # radius, with no overlap. self.r_host = inputs['r_host'] self.th_product = inputs['th_product'] + + self.V_host = 4./3. * np.pi * (self.r_host)**3 # Volume of a single carbon / host particle [m3] self.A_host = 4. * np.pi * (self.r_host)**2 # Surface area of a single carbon / host particle [m2] # m2 of host-electrolyte interface / m3 of total volume [m^-1] @@ -93,8 +97,7 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, # Microstructure-based transport scaling factor, based on Bruggeman # coefficient of -0.5: - self.elyte_microstructure = self.eps_elyte_init[0]**1.5 #This is used by the - + self.elyte_microstructure = self.eps_elyte_init[0]**1.5 #This is used by the electrode_boundary_flux # SV_offset specifies the index of the first SV variable for the # electode (zero for anode, n_vars_anode + n_vars_sep for the cathode) self.SV_offset = offset @@ -219,7 +222,7 @@ def residual(self, t, SV, SVdot, sep, counter, params): c_k_elyte = SV_loc[SVptr['C_k_elyte'][j]] eps_product = SV_loc[SVptr['eps_product'][j]] eps_elyte = 1. - eps_product - self.eps_host[j] - + self.elyte_microstructure = eps_elyte**1.5 # Read electrolyte fluxes at the separator boundary. No matter the # electrode, the function returns a value where flux to the electrode # is considered positive. We multiply by `i_ext_flag` to get the @@ -270,7 +273,7 @@ def residual(self, t, SV, SVdot, sep, counter, params): resid[SVptr['phi_ed'][j]] = i_io_in - i_io_out + i_el_in - i_el_out # Calculate available surface area (m2 interface per m3 electrode): - A_avail = self.A_init[j] - eps_product/self.th_product + A_avail = self.A_init[j] - eps_product/self.thickness_temp#self.th_product # Convert to m2 interface per m2 geometric area: A_surf_ratio = A_avail*self.dy # Multiplier to scale phase destruction rates. As eps_product @@ -363,7 +366,7 @@ def residual(self, t, SV, SVdot, sep, counter, params): - params['potential']) # Calculate available surface area (m2 interface per m3 electrode): - A_avail = self.A_init[j] - eps_product/self.th_product + A_avail = self.A_init[j] - eps_product/self.thickness_temp#self.th_product # Convert to m2 interface per m2 geometric area: A_surf_ratio = A_avail*self.dy @@ -413,6 +416,9 @@ def residual(self, t, SV, SVdot, sep, counter, params): + sdot_elyte_host * A_surf_ratio) * self.dyInv)/eps_elyte + eps_elyte = 1. - eps_product - self.eps_host[0] + self.elyte_microstructure = eps_elyte**1.5 + return resid def voltage_lim(self, SV, val): diff --git a/electrode_models/dense_electrode.py b/electrode_models/dense_electrode.py index 9d89932..2eaa3be 100755 --- a/electrode_models/dense_electrode.py +++ b/electrode_models/dense_electrode.py @@ -12,20 +12,20 @@ class electrode(): """ Create an electrode object representing the dense electrode """ - def __init__(self, input_file, inputs, sep_inputs, counter_inputs, + def __init__(self, input_file, inputs, sep_inputs, counter_inputs, electrode_name, params, offset): """ Initialize the model. - """ + """ # Import relevant Cantera objects. self.bulk_obj = ct.Solution(input_file, inputs['bulk-phase']) self.elyte_obj = ct.Solution(input_file, inputs['electrolyte-phase']) self.conductor_obj = ct.Solution(input_file, inputs['conductor-phase']) - self.surf_obj = ct.Interface(input_file, inputs['surf-phase'], + self.surf_obj = ct.Interface(input_file, inputs['surf-phase'], [self.bulk_obj, self.elyte_obj, self.conductor_obj]) - # Anode or cathode? Positive external current delivers positive charge + # Anode or cathode? Positive external current delivers positive charge # to the anode, and removes positive charge from the cathode. self.name = electrode_name if self.name=='anode': @@ -35,13 +35,13 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, else: raise ValueError("Electrode must be an anode or a cathode.") - # Store the species index of the Li ion in the Cantera object for the + # Store the species index of the Li ion in the Cantera object for the # electrolyte phase: self.index_Li_elyte = self.elyte_obj.species_index(inputs['mobile-ion']) # Electrode thickness self.dy = inputs['thickness'] - # The electrode consumption rate quickly goes to zero, below a + # The electrode consumption rate quickly goes to zero, below a # user-specified minimum thickness: self.min_thickness = inputs['minimum-thickness'] @@ -51,35 +51,35 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, # Inverse of the double layer capacitance, per unit interface area: self.C_dl_Inv = 1/inputs['C_dl'] - # Thickness of separator node considered as part of the anode domain. + # Thickness of separator node considered as part of the anode domain. # This is "subtracted" off from the total separator thickness. self.dy_elyte = inputs['dy_elyte'] - + # Electrolyte volume fraction in the separator: self.eps_elyte = sep_inputs['eps_electrolyte'] - # Microstructure-based transport scaling factor, based on Bruggeman + # Microstructure-based transport scaling factor, based on Bruggeman # coefficient of -0.5: self.elyte_microstructure = self.eps_elyte**1.5 - - # SV_offset specifies the index of the first SV variable for the + + # SV_offset specifies the index of the first SV variable for the # electode (zero for anode, n_vars_anode + n_vars_sep for the cathode) self.SV_offset = offset - - # Dense Li is not capacity-limiting, in and of itself. Rather the - # total amount of Li in the system is the limit. This is done in a - # separate routine, at a later time. Provide a large placeholder number + + # Dense Li is not capacity-limiting, in and of itself. Rather the + # total amount of Li in the system is the limit. This is done in a + # separate routine, at a later time. Provide a large placeholder number # here, so that it will not be the minimum, when evaluated later: self.capacity = 1e21 - - # Mumber of state variables: electrode potential, electrolyte + + # Mumber of state variables: electrode potential, electrolyte # potential, thickness, electrolyte composition (n_species) self.n_vars = 3 + self.elyte_obj.n_species # This model produces one plot, for the electrode thickness: self.n_plots = 1 - # Set the Cantera object state. + # Set the Cantera object state. self.bulk_obj.electric_potential = inputs['phi_0'] # If the user provided an initial composition, use that, here: if 'X_0' in inputs: @@ -99,15 +99,15 @@ def residual(self, t, SV, SVdot, sep, counter, params): 1. The electric potential is an algebraic variable. In the anode, phi = 0 is the reference potential for the system. - In the cathode, the electric potential must be such that the ionic current is spatially in_variant (i.e. it is constant and equal to the external applied current, for galvanostatic simulations). + In the cathode, the electric potential must be such that the ionic current is spatially in_variant (i.e. it is constant and equal to the external applied current, for galvanostatic simulations). The residual corresponding to these variables (suppose an index 'j') are of the form: resid[j] = (epression equaling zero) 2. All other variables are governed by differential equations. - - We have a means to calculate dSV[j]/dt for a state variable SV[j] (state variable with index j). - + + We have a means to calculate dSV[j]/dt for a state variable SV[j] (state variable with index j). + The residuals corresponding to these variables will have the form: resid[j] = SVdot[j] - (expression equalling dSV/dt) @@ -119,7 +119,7 @@ def residual(self, t, SV, SVdot, sep, counter, params): - counter: the object representing the electrode counter to the current electrode - params: dict of battery simulation parameters. """ - + # Initialize the residual array: resid = np.zeros((self.n_vars,)) @@ -136,8 +136,8 @@ def residual(self, t, SV, SVdot, sep, counter, params): self.bulk_obj.electric_potential = phi_ed self.conductor_obj.electric_potential = phi_ed self.elyte_obj.electric_potential = phi_elyte - - # Multiplier on the electrode removal reaction. Quickly goes to zero, + + # Multiplier on the electrode removal reaction. Quickly goes to zero, # for thicknesses below a user-specified minimum: mult = tanh(SV_loc[SVptr['thickness']]/self.min_thickness) @@ -146,28 +146,28 @@ def residual(self, t, SV, SVdot, sep, counter, params): (mult*self.surf_obj.get_creation_rates(self.conductor_obj) - self.surf_obj.get_destruction_rates(self.conductor_obj)) - # Molar production rate of electrode species (kmol/m2/s). Here, we scale + # Molar production rate of electrode species (kmol/m2/s). Here, we scale # the destruction rate by our multiplier. sdot_electrode = (self.surf_obj.get_creation_rates(self.bulk_obj) - mult*self.surf_obj.get_destruction_rates(self.bulk_obj)) - # Faradaic current density is positive when electrons are consumed + # Faradaic current density is positive when electrons are consumed # (Li transferred to the anode) i_Far = -ct.faraday*sdot_electron - # Flux of electrolyte species between the separator and the electrolyte + # Flux of electrolyte species between the separator and the electrolyte # in the current electrode domain: N_k_sep, i_io = sep.electrode_boundary_flux(SV, self, params['T']) # Electrode electric potential if self.name=='anode': - # For the anode, the electric potential is an algebraic variable, + # For the anode, the electric potential is an algebraic variable, # always equal to zero: resid[SVptr['phi_ed']] = SV_loc[SVptr['phi_ed']] elif self.name=='cathode': - # For the cathode, the potential of the cathode must be such that - # the electrolyte electric potential (calculated as phi_ca + + # For the cathode, the potential of the cathode must be such that + # the electrolyte electric potential (calculated as phi_ca + # dphi_dl) produces the correct ionic current between the separator # and cathode: if params['boundary'] == 'current': resid[SVptr['phi_ed']] = i_io - params['i_ext'] @@ -177,34 +177,34 @@ def residual(self, t, SV, SVdot, sep, counter, params): # Cell potential must equal phi: resid[SVptr['phi_ed']] = SV_loc[SVptr['phi_ed']] - phi - - # Double layer current has the same sign as i_Far, and is based on + + # Double layer current has the same sign as i_Far, and is based on # charge balance in the electrolyte phase: i_dl = self.i_ext_flag*i_io/self.A_surf_ratio - i_Far - + # Differential equation for the double layer potential difference: - resid[SVptr['phi_dl']] = (SVdot_loc[SVptr['phi_dl']] + resid[SVptr['phi_dl']] = (SVdot_loc[SVptr['phi_dl']] - i_dl*self.C_dl_Inv) - + # Change in thickness per time: dH_dt = np.dot(sdot_electrode, self.bulk_obj.partial_molar_volumes) resid[SVptr['thickness']] = SVdot_loc[SVptr['thickness']] - dH_dt - - # Set time derivatives for electrolyte species concentrations to zero + + # Set time derivatives for electrolyte species concentrations to zero # (temporary) - # Molar production rate of electrode species (kmol/m2/s). Here, we scale + # Molar production rate of electrode species (kmol/m2/s). Here, we scale # the destruction rate by our multiplier. sdot_electrolyte = \ (mult*self.surf_obj.get_creation_rates(self.elyte_obj) - self.surf_obj.get_destruction_rates(self.elyte_obj)) - # Double layer current removes Li from the electrolyte. Add this to + # Double layer current removes Li from the electrolyte. Add this to # sdot_electrolyte: sdot_electrolyte[self.index_Li_elyte] -= i_dl / ct.faraday dCk_elyte_dt = \ - (sdot_electrolyte * self.A_surf_ratio - + self.i_ext_flag * N_k_sep) / self.dy_elyte + (sdot_electrolyte * self.A_surf_ratio + + self.i_ext_flag * N_k_sep) / self.dy_elyte resid[SVptr['C_k_elyte']] = SVdot_loc[SVptr['C_k_elyte']] - dCk_elyte_dt return resid @@ -219,24 +219,24 @@ def initialize(self, inputs, sep_inputs): self.SVptr['phi_ed'] = np.array([0]) self.SVptr['phi_dl'] = np.array([1]) self.SVptr['thickness'] = np.array([2]) - self.SVptr['C_k_elyte'] = np.arange(3, + self.SVptr['C_k_elyte'] = np.arange(3, 3 + self.elyte_obj.n_species) - self.SVnames = (['phi_ed', 'phi_dl', 'thickness'] + self.SVnames = (['phi_ed', 'phi_dl', 'thickness'] + self.elyte_obj.species_names[:]) - + # There is only one node, but give the pointer a shape so that SVptr # ['C_k_elyte'][j] accesses the pointer array: self.SVptr['C_k_elyte'].shape = (1, self.elyte_obj.n_species) - # A pointer to where the SV variables for this electrode are, within + # A pointer to where the SV variables for this electrode are, within # the overall solution vector for the entire problem: - self.SVptr['electrode'] = np.arange(self.SV_offset, + self.SVptr['electrode'] = np.arange(self.SV_offset, self.SV_offset+self.n_vars) # Save the SV indices of any algebraic variables: self.algvars = self.SV_offset + self.SVptr['phi_ed'][:] - + # Load intial state variable values: SV[self.SVptr['phi_ed']] = inputs['phi_0'] SV[self.SVptr['phi_dl']] = sep_inputs['phi_0'] - inputs['phi_0'] @@ -252,8 +252,8 @@ def voltage_lim(self, SV, val): # Save local copies of the solution vector and pointers for this electrode: SVptr = self.SVptr SV_loc = SV[SVptr['electrode']] - - # Calculate the current voltage, relative to the limit. The simulation + + # Calculate the current voltage, relative to the limit. The simulation # looks for instances where this value changes sign (i.e. crosses zero) voltage_eval = SV_loc[SVptr['phi_ed']] - val @@ -269,9 +269,9 @@ def species_lim(self, SV, val): # For each electrode point, find the minimum species concentration, and # compare to the user-provided minimum. Save only the minimum value: species_eval = min(SV_loc[SVptr['C_k_elyte'][0]]) - val - - # The simulation looks for instances where this value changes sign - # (i.e. where it crosses zero) + + # The simulation looks for instances where this value changes sign + # (i.e. where it crosses zero) return species_eval def adjust_separator(self, sep): @@ -279,17 +279,17 @@ def adjust_separator(self, sep): The electrode domain considers the electrode object plus a thin layer of the separator, adjacent to the self. We subtract this thickness from the total separator thickness, so that we do not inadvertently increase the total transport resistance through the separator. """ sep.dy -= self.dy_elyte - - # Reduce the number of points in the separator by one, unless the - # separator already only contains one point (which is the case for the + + # Reduce the number of points in the separator by one, unless the + # separator already only contains one point (which is the case for the # `ionic_resistor` model. In this case, leave sep.n_points at 1.) sep.n_points = max(sep.n_points - 1, 1) - + return sep def output(self, axs, solution, SV_offset, ax_offset): - axs[ax_offset].plot(solution[0,:]/3600, + axs[ax_offset].plot(solution[0,:]/3600, 1e6*solution[SV_offset+int(self.SVptr['thickness'])]) axs[ax_offset].set_ylabel(self.name+' Thickness \n($\mu$m)') diff --git a/inputs/LiO2_Gittleson.yaml b/inputs/LiO2_Gittleson.yaml new file mode 100644 index 0000000..bed7a0b --- /dev/null +++ b/inputs/LiO2_Gittleson.yaml @@ -0,0 +1,499 @@ +#=============================================================================== +# input_template.yaml +# +# User inputs for bat-can model creation and running. +#=============================================================================== + +cell-description: +#This section can have any attributes just do whatever you want + anode: + class: 'dense_electrode' + bulk-phase: 'lithium_metal' # Cantera phase name + surf-phase: 'lithium_electrolyte' # Cantera phase name + electrolyte-phase: 'electrolyte' # Cantera phase name + conductor-phase: 'electron' # Cantera phase name + stored-ion: # Details on the stored species. Used for capacity calc. + name: Li[metal] + charge: 1 + mobile-ion: 'Li+[elyt]' # Species name for Li+ in elyte. + thickness: 2.5e-4 # anode thickness, m + minimum-thickness: 1e-9 # minimum anode thickness. + phi_0: 0. # Initial electric potential, V + A_surf_ratio: 1.0 # Interface area per unit geometric area + C_dl: 6e-4 #F/m2 + dy_elyte: 2e-6 #maybe should remove this? + separator: + class: 'porous_separator' + thickness: 2.6e-4 # separator thickness, m + n_points: 5 # Number of finite volumes to discretize + electrolyte-phase: 'electrolyte' # Cantera phase name + sigma_io: 0.74 # S/m for 1M LiTFSI in TEGDME at 25 deg C // Chen, 2019 + phi_0: 3.117 # Initial electric potential, V // Yin, 2017 + eps_electrolyte: 0.65 # Electrolyte volume fraction, - + transport: + flag_lithiated: 0 + D_scale_coeff: 1e-11 + diffusion-scaling: 'ideal' + mobile-ion: 'Li+[elyt]' + model: 'dilute-solution' + diffusion-coefficients: # Species names must match those in the phase + #definition, below: # fix diffusion coefficients later... + - species: 'C2H6OS[elyt]' + D_k: 2.50e-10 # m2 S-1 // Saito, 2017 + - species: 'Li+[elyt]' + D_k: 6.30e-10 # m2 S-1 // Saito, 2017 + - species: 'TFSI-[elyt]' + D_k: 1.40e-10 # m2 S-1 // Saito, 2017 + - species: 'O2(e)' + D_k: 5.67e-9 #11 # m2 s-1 // Gittleson, 2017 for DMSO + # - species: 'Li2O2[elyt]' + # D_k: 2e-13 # can't find (assume low?) + cathode: + class: 'air_electrode' + n-points: 5 #15 + gas-phase: 'air' + elyte-iphase: 'air-elyte-interface' + electrolyte-phase: 'electrolyte' # Cantera phase name + host-phase: 'conductor' # Cantera phase name + surf-iphase: 'cathode-surf' # Cantera phase name + product-phase: 'li-oxide' + product-phase-min: 1e-6 # Minimum volume fraction for product phase + mobile-ion: 'Li+[elyt]' # Species name for Li+ in elyte. + thickness: 2.350e-4 # cathode thickness, m + r_host: 1.e-6 # Carbon Particle radius, m + host_form: 'cylinder' # host morphology. Options are 'cylinder' and 'sphere' + host_surf_frac: 0.95 # Fraction of host surface area that is accessible. + d_oxide: 1.0e-8 #Oxide diameter, m + th_product: 220.e-9 #185e-9 #oxide thickness, m + phi_0: 3.15 # Initial electric potential, V // Yin, 2017 + sigma_el: 100 #electrical ondutivity of the host + eps_host: 0.2 # Solid phase volume fraction of carbon - + eps_product: 0.00 # solid phase volume fraction of Li2O2 - + C_dl: 6e-2 #F/m2 + stored-species: # Details on the stored species. Used for capacity calc. + name: Li2O2[cathode] + charge: 2 + MW: 45.881 + plot-species: # any extra species in the electrolyte that you want to plot. + - name: 'O2(e)' + + +# Simulation parameters: +parameters: + T: 25 C + P: 3.5 atm + # Describe simulation type, parameters, etc. + species-default: {'C2H6OS[elyt]': 0.8007, 'Li+[elyt]': 1.E-12, 'TFSI-[elyt]': 0.109610, 'O2(e)': 1.E-12} # Replacement mole fractions if elyte composition goes to NaN. + # initialize: # Calculate a common inititial condition for all steps? + # enable: True + # type: 'open-circuit' + # time: 5 # s + simulations: + - type: 'CC_cycle' # Constant current cycling + # Specify only one of i_ext or C-rate. The other should be null: + i_ext: 2.546 A/m # input current density, format: XX units. Units are 'current/area', with no carat for exponents (e.g. A/cm2) + C-rate: null # input C-rate + n_cycles: 0.5 # Number of cycles. + first-step: 'discharge' # Start with charge or discharge? + phi-cutoff-lower: 2.0 # Simulation cuts off if E_Cell <= this value + phi-cutoff-upper: 4.8 # Simulation cuts off if E_Cell >= this value + species-cutoff: 1.e-18 # Simulation cuts off if C_k <= this value + equilibrate: + enable: True # Start with a hold at open circuit? (True or False) + time: 5 # If true, how long is the hold, s + outputs: + show-plots: True # Show the plots and save them (True), or just save + # them (False)? + save-name: '5p4_Gittleson' # Folder label for output files. + validation: 'Gittleson_LiO2/Gittleson_5p42.xls' #Location of validation data + - type: 'CC_cycle' # Constant current cycling + # Specify only one of i_ext or C-rate. The other should be null: + i_ext: 10.19 A/m2 #0.01 A/cm2 # input current density, format: XX units. Units are 'current/area', with no carat for exponents (e.g. A/cm2) + C-rate: null # input C-rate + n_cycles: 0.5 # Number of cycles. + first-step: 'discharge' # Start with charge or discharge? + phi-cutoff-lower: 2.0 # Simulation cuts off if E_Cell <= this value + phi-cutoff-upper: 4.8 # Simulation cuts off if E_Cell >= this value + # species-cutoff: 1.e-14 # Simulation cuts off if C_k <= this value + species-default: {'C2H6OS[elyt]': 0.9007, 'Li+[elyt]': 1.E-12, 'TFSI-[elyt]': 0.049610, 'O2(e)': 1.E-12} # Replacement mole fractions if elyte composition goes to NaN. + equilibrate: + enable: False # Start with a hold at open circuit? (True or False) + time: 10 # If true, how long is the hold, s + outputs: + show-plots: True # Show the plots and save them (True), or just save + # them (False)? + save-name: '21p7_Gittleson' # Folder label for output files. + validation: 'Gittleson_LiO2/Gittleson_21p7.xls' #Location of validation data + # - type: 'CC_cycle' # Constant current cycling + # # Specify only one of i_ext or C-rate. The other should be null: + # i_ext: 20.384 A/m2 #0.01 A/cm2 # input current density, format: XX units. Units are 'current/area', with no carat for exponents (e.g. A/cm2) + # C-rate: null # input C-rate + # n_cycles: 0.5 # Number of cycles. + # first-step: 'discharge' # Start with charge or discharge? + # phi-cutoff-lower: 1.75 # Simulation cuts off if E_Cell <= this value + # phi-cutoff-upper: 4.8 # Simulation cuts off if E_Cell >= this value + # # species-cutoff: 0.000565 # Simulation cuts off if C_k <= this value + # species-default: {'C2H6OS[elyt]': 0.9007, 'Li+[elyt]': 0.049610, 'TFSI-[elyt]': 0.049610, 'O2(e)': 1.E-12} # Replacement mole fractions if elyte composition goes to NaN. + # equilibrate: + # enable: False # Start with a hold at open circuit? (True or False) + # time: 5 # If true, how long is the hold, s + # outputs: + # show-plots: True # Show the plots and save them (True), or just save + # # them (False)? + # save-name: '43p4_Gittleson' # Folder label for output files. + # validation: 'Gittleson_LiO2/Gittleson_43p4.xls' #Location of validation data + + +fit-parameters: + - type: 'elyte-transport' + species: 'O2(e)' + guess-value: 1.67e-9 + min: 1.67e-9 + max: 4.70e-9 + # - type: 'elyte-transport' + # species: 'Li+[elyt]' + # guess-value: 1.10e-10 + # min: 1.82e-11 + # max: 1.50e-10 + # - type: 'cathode-kinetics' + # reaction-index: 0 + # guess-value: 0.46210493 # This is a multiplier to scale the value in the mech. + # min: 0.1 + # max: 10 + # - type: 'cathode-microstructure' + # parameter: 'th_product' + # guess-value: 220.e-9 + # min: 100.e-9 + # max: 420.e-9 + + + +# Cantera inputs: +description: |- + Cantera input file for an Li-Oxygen battery + + This file includes a full set of thermodynamic and kinetic parameters of a + lithium-oxygen battery, in particular: + - Active materials: Li2O2 and Li (li metal) + - Organic electrolyte: TEGDME with 1M LiPF6 + - Interfaces: LCO/electrolyte and Li/electrolyte + - Cathode: The cathode is comprised of catalyst, carbon, (will binder be a later addition.) + - Charge-transfer reactions at the two interfaces + + Bulk phases + =========== + Lithium (anode) + + Ether based electrolyte (electrolyte) + Solvent: Tetraethylene glycol dimethyl ether + Salt: 0.25M LiTFSI + + Interface phases + ================ + + lithium/electrolyte interface (lithium_electrolyte) + Species and site density are dummy entries (as we do not consider surface- + adsorbed species) + + LCO/electrolyte interface (cathode_electrolyte) + Species and site density are dummy entries (as we do not consider surface- + adsorbed species) + + Builds on prior code development (LeBar, 2019) +generator: cti2yaml +cantera-version: 2.5.0 +date: Wed, 11 Dec 2019 16:59:07 -0500 +input-files: [lithium_ion_battery.cti] + +# ================================== +# PHASE DEFINITIONS +# ================================== +phases: + +- name: lithium_metal + thermo: ideal-condensed + species: ['Li[metal]'] + density: 0.534 g/cm^3 #verified + state: + T: 300.0 + P: 1 atm + +- name: electron + thermo: electron-cloud + elements: [E] + species: [electron] + state: + X: {electron: 1.0} + density: 1.0 kg/m^3 + +- name: lithium_electrolyte + thermo: ideal-surface + species: [(dummy)] + kinetics: surface + reactions: [lithium_electrolyte-reactions] + site-density: 0.01 mol/cm^2 #Question: should I change this?? + +- name: electrolyte + thermo: ideal-condensed + elements: [Li, Cl, C, H, O, E, S, F, N] + species: ['C2H6OS[elyt]', 'Li+[elyt]', 'TFSI-[elyt]', 'O2(e)'] #, 'Li2O2[elyt]'] + state: + X: {'C2H6OS[elyt]': 0.8207, 'Li+[elyt]': 0.089610, 'TFSI-[elyt]': 0.089610, 'O2(e)': 7.7E-5} #,'O2(e)': 1.39E-4} 'Li2O2[elyt]': 1E-9} + standard-concentration-basis: unity + +- name: air-elyte-interface + thermo: ideal-surface + species: [(dummy)] + kinetics: surface + reactions: [air-electrolyte-reactions] + site-density: 0.01 mol/cm^2 + +- name: air + thermo: ideal-gas + elements: [O, N, Ar] + species: [O2(g), N2(g), AR(g)] + kinetics: gas + reactions: none + transport: mixture-averaged + state: {T: 300.0, P: 1 atm, X: {O2(g): 0.21, N2(g): 0.78, AR(g): 0.01}} + +- name: conductor + thermo: electron-cloud + elements: [E] + species: [electron] + state: + X: {electron: 1.0} + density: 1.0 kg/m^3 + +- name: cathode-surf + thermo: ideal-surface + species: [(dummy)] + kinetics: surface + reactions: [cathode-surf-reactions] + state: + T: 300.0 + P: 1 atm + site-density: 2.50e-5 mol/cm^2 + +- name: li-oxide + composition: {Li: 2, O: 2} + thermo: fixed-stoichiometry + species: ['Li2O2[cathode]'] + density: 2310 kg/m^3 # https://pubchem.ncbi.nlm.nih.gov/compound/Lithium-peroxide-_Li2_O2#section=Physical-Description + state: + T: 300.0 + P: 1 atm + +# ================================== +# SPECIES DEFINITIONS +# ================================== +species: +- name: Li[metal] + composition: {Li: 1} + thermo: + model: constant-cp + h0: 19.50 kJ/mol # Question: unverified - source? + s0: 29.1 J/mol/K # verified NIST webbook + equation-of-state: + model: constant-volume + molar-volume: 12.998 cm^3/mol #verified + +- name: electron + composition: {E: 1} + thermo: + model: constant-cp + h0: 0.0 kJ/mol + s0: 0.0 J/mol/K + note: |- + 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) + +- name: (dummy) + composition: {} + thermo: + model: constant-cp + h0: 0.0 kJ/mol + s0: 0.0 J/mol/K + note: Dummy species (needed for defining the interfaces) + +- name: C2H6OS[elyt] + composition: {C: 2, H: 6, O: 1, S: 1} + thermo: + model: constant-cp + h0: -203.4 kJ/mol # https://webbook.nist.gov/cgi/cbook.cgi?ID=C67685&Mask=2 + s0: 0.0 J/mol/K # unknown + equation-of-state: + model: constant-volume + molar-volume: 71.001454016721192293711377680843 cm^3/mol + note: |- + Dimethyl sulfoxide, MW: 78.13 g/mol + Density of tetraethylene glycol dimethyl ether: 1100.4 kg/m3 (used to calculate species molar volume + as molecular weight (MW)/density) + Taken from: https://www.sigmaaldrich.com/US/en/sds/sigma/d8418 + Molar enthalpy and entropy set to zero (dummy entries as this species does + not participate in chemical reactions) + + +- name: Li+[elyt] + composition: {Li: 1, E: -1} + thermo: + model: constant-cp + h0: -278.49 kJ/mol + s0: 13.4 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 5.508297619047619 cm^3/mol + note: |- #Question: why is this using the density of electrolyte + 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) + +- name: TFSI-[elyt] + composition: {C: 2, F: 6, N: 1, O: 4, S: 1, E: 1} + thermo: + model: constant-cp + h0: 0.0 J/mol + s0: 0.0 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 78.92857142857143 cm^3/mol + note: |- #Question: why is this using the density of electrolyte + bis(trifluoromethylsulphonyl)imide ion, MW: 280.147 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) + +- name: Li2O2[elyt] + composition: {Li: 2, O: 2} + thermo: + model: constant-cp + h0: 0.0 J/mol # Question: should I put this here https + s0: 0.0 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 19.344262295081968261 cm^3/mol + note: |- + Li2O2 particle, MW: 38.94 g/mol + Density of electrolyte: 2013 kg/m3 (used to calculate species molar volume + as molecular weight (MW)/density) + https://www.chemspider.com/Chemical-Structure.145811.html + Molar enthalpy and entropy set to zero (dummy entries as this species does + not participate in chemical reactions) + +- name: O2(e) + composition: {O: 2} + thermo: + model: constant-cp + h0: -40.375 kJ/mol + s0: 0. J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 28.0445 cm^3/mol + note: |- + An O2 molecule dissolved into electrolyte. Thermo tuned to get correct solubility at room temperature. + Molar volume estimated from liquid O2 values. https://en.wikipedia.org/wiki/Liquid_oxygen#cite_note-3 + +- name: O2(g) + composition: {O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, + -1063.94356, 3.65767573] + - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, + -1088.45772, 5.45323129] + note: TPIS89 + transport: + model: gas + geometry: linear + well-depth: 107.4 + diameter: 3.458 + polarizability: 1.6 + rotational-relaxation: 3.8 +- name: N2(g) + composition: {N: 2} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.298677, 1.4082404e-03, -3.963222e-06, 5.641515e-09, -2.444854e-12, + -1020.8999, 3.950372] + - [2.92664, 1.4879768e-03, -5.68476e-07, 1.0097038e-10, -6.753351e-15, + -922.7977, 5.980528] + note: '121286' + transport: + model: gas + geometry: linear + well-depth: 97.53 + diameter: 3.621 + polarizability: 1.76 + rotational-relaxation: 4.0 +- name: AR(g) + composition: {Ar: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366] + - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366] + note: '120186' + transport: + model: gas + geometry: atom + well-depth: 136.5 + diameter: 3.33 + +- name: Li2O2[cathode] + composition: {Li: 2, O: 2} + thermo: + model: constant-cp + h0: -600. kJ/mol #-634.3 kJ/mol #Haynes, 2016 + s0: 0.0 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 19.861904761904753514 cm^3/mol + note: |- + h0 value tuned to give correct open circuit potential, relative to li metal anode. + Lithium peroxide, MW: 45.881 g/mol. + Density of L2O2: 2310 kg/m3 + (Haynes, 2016) + +lithium_electrolyte-reactions: +- equation: Li[metal] <=> Li+[elyt] + electron + id: lithium_faradaic_reaction + rate-constant: {A: 6.0e+20, b: 0.0, Ea: 0.0} #{A: 6.0e+09, b: 0.0, Ea: 0.0} + beta: 0.5 + +cathode-surf-reactions: +- equation: 2 Li+[elyt] + O2(e) + 2 electron <=> Li2O2[cathode] # Reaction 1 + id: cathode_reaction + rate-constant: {A: 9.242e-1, b: 0.0, Ea: 58.0 kJ/mol} #This is already scaled by 0.001 from fit params. #need a new reaction rate #5.629e+6, b: 0.0, Ea: 58.0 kJ/mol} #need a new reaction rate + exchange-current-density-formulation: true + beta: 0.5 + +air-electrolyte-reactions: +- equation: O2(g) <=> O2(e) # Reaction 2 + id: cathode_gas_reaction + rate-constant: {A: 2.e-2, b: 0.0, Ea: 0 kJ/mol} + +#"""Citations""" +#Chen, J.; et al., ACS Omega, vol. 4, no. 20708-20714, 2019, DOI: 10.1021/acsomega.9b02941. + +#Gittleson, F; et al, Energy Environ. Sci., 2017,10, 1167-1179, DOI: 10.1039/C6EE02915A. + +#Haynes, W.M.; CRC handbook of chemistry and physics, 2016, 97th Edition / Boca Raton, Florida: CRC Press. + +#LeBar, A., Colorado School of Mines, 2019, URL: https://mountainscholar.org/handle/11124/173373 + +#Yin, Y.; et al., J. Phys. Chem. C, vol. 121, no. 36, pp. 19577–19585, 2017, doi: 10.1021/acs.jpcc.7b05224. + +#Official Soundtrack: +#Day 6 - as of 2021 all albums/songs +# - highlight Finale, Congradulations, 예뻤어 +# - the Magnetic Fields, 69 Love Songs, vol III +#Julia Rocka - Blaza \ No newline at end of file diff --git a/inputs/LiO2_Liu_starter.yaml b/inputs/LiO2_Liu_starter.yaml new file mode 100644 index 0000000..c42f30f --- /dev/null +++ b/inputs/LiO2_Liu_starter.yaml @@ -0,0 +1,520 @@ +#=============================================================================== +# input_template.yaml +# +# User inputs for bat-can model creation and running. +#=============================================================================== + +cell-description: +#This section can have any attributes just do whatever you want + anode: + class: 'dense_electrode' + bulk-phase: 'lithium_metal' # Cantera phase name + surf-phase: 'lithium_electrolyte' # Cantera phase name + electrolyte-phase: 'electrolyte' # Cantera phase name + conductor-phase: 'electron' # Cantera phase name + stored-ion: # Details on the stored species. Used for capacity calc. + name: Li[metal] + charge: 1 + mobile-ion: 'Li+[elyt]' # Species name for Li+ in elyte. + thickness: 2.5e-4 # anode thickness, m + minimum-thickness: 1e-9 # minimum anode thickness. + phi_0: 0. # Initial electric potential, V + A_surf_ratio: 1.0 # Interface area per unit geometric area + C_dl: 6e-4 #F/m2 + dy_elyte: 2e-6 #maybe should remove this? + separator: + class: 'porous_separator' + thickness: 6.5e-4 # separator thickness, m + n_points: 5 # Number of finite volumes to discretize + electrolyte-phase: 'electrolyte' # Cantera phase name + sigma_io: 0.74 # S/m for 1M LiTFSI in TEGDME at 25 deg C // Chen, 2019 + phi_0: 3.117 # Initial electric potential, V // Yin, 2017 + eps_electrolyte: 0.65 # Electrolyte volume fraction, - + transport: + flag_lithiated: 0 + D_scale_coeff: 1e-11 + diffusion-scaling: 'ideal' + mobile-ion: 'Li+[elyt]' + model: 'dilute-solution' + diffusion-coefficients: # Species names must match those in the phase + #definition, below: # fix diffusion coefficients later... + - species: 'C2H6O[elyt]' + D_k: 1.43e-10 # m2 S-1 // Saito, 2017 + - species: 'Li+[elyt]' + D_k: 6.42e-11 # m2 S-1 // Saito, 2017 + - species: 'TFSI-[elyt]' + D_k: 6.79e-10 # m2 S-1 // Saito, 2017 + - species: 'O2(e)' + D_k: 24.30E-10 #11 # m2 s-1 // Gittleson, 2017 for DMSO + # - species: 'Li2O2[elyt]' + # D_k: 2e-13 # can't find (assume low?) + cathode: + class: 'air_electrode' + n-points: 20 + gas-phase: 'air' + elyte-iphase: 'air-elyte-interface' + electrolyte-phase: 'electrolyte' # Cantera phase name + host-phase: 'conductor' # Cantera phase name + surf-iphase: 'cathode-surf' # Cantera phase name + product-phase: 'li-oxide' + product-phase-min: 1e-6 # Minimum volume fraction for product phase + mobile-ion: 'Li+[elyt]' # Species name for Li+ in elyte. + thickness: 2.350e-4 # cathode thickness, m + r_host: 0.75e-6 # Carbon Particle radius, m + host_form: 'cylinder' # host morphology. Options are 'cylinder' and 'sphere' + host_surf_frac: 0.98 # Fraction of host surface area that is accessible. + d_oxide: 1.0e-8 #Oxide diameter, m + th_product: 87.e-9 #185e-9 #oxide thickness, m + phi_0: 3.15 # Initial electric potential, V // Yin, 2017 + sigma_el: 100 #electrical condutivity of the host + eps_host: [0.2] # Solid phase volume fraction of carbon -0.2, 0.4, 0.6, 0.7, 0.8 + eps_product: 0.00 # solid phase volume fraction of Li2O2 - + C_dl: 6e-2 #F/m2 + stored-species: # Details on the stored species. Used for capacity calc. + name: Li2O2[cathode] + charge: 2 + MW: 45.881 + plot-species: # any extra species in the electrolyte that you want to plot. + - name: 'O2(e)' + + +# Simulation parameters: +parameters: + T: 25 C + P: 3.5 atm + # Describe simulation type, parameters, etc. + # initialize: # Calculate a common inititial condition for all steps? + # enable: True + # type: 'open-circuit' + # time: 3600 # s + simulations: + - type: 'CC_cycle' # Constant current cycling + # Specify only one of i_ext or C-rate. The other should be null: + i_ext: 0.1 mA/cm2 #0.01 A/cm2 # input current density, format: XX units. Units are 'current/area', with no carat for exponents (e.g. A/cm2) + C-rate: null # input C-rate + n_cycles: 0.5 # Number of cycles. + first-step: 'discharge' # Start with charge or discharge? + phi-cutoff-lower: 2.0 # Simulation cuts off if E_Cell <= this value + phi-cutoff-upper: 4.8 # Simulation cuts off if E_Cell >= this value + species-cutoff: 1.e-18 # Simulation cuts off if C_k <= this value + species-default: {'C2H6O[elyt]': 0.9007, 'Li+[elyt]': 1.E-12, 'TFSI-[elyt]': 0.049610, 'O2(e)': 1.E-18} # Replacement mole fractions if elyte composition goes to NaN. + equilibrate: + enable: True # Start with a hold at open circuit? (True or False) + time: 5 # If true, how long is the hold, s + outputs: + show-plots: True # Show the plots and save them (True), or just save + # them (False)? + save-name: '0p1_liu' # Folder label for output files. + validation: 'Liu_O2/0.1mAhcm2.xlsx' #Location of validation data + thickness: 180.e-9 + - type: 'CC_cycle' # Constant current cycling + # Specify only one of i_ext or C-rate. The other should be null: + i_ext: 0.2 mA/cm2 #0.01 A/cm2 # input current density, format: XX units. Units are 'current/area', with no carat for exponents (e.g. A/cm2) + C-rate: null # input C-rate + n_cycles: 0.5 # Number of cycles. + first-step: 'discharge' # Start with charge or discharge? + phi-cutoff-lower: 2.0 # Simulation cuts off if E_Cell <= this value + phi-cutoff-upper: 4.8 # Simulation cuts off if E_Cell >= this value + # species-cutoff: 1.e-14 # Simulation cuts off if C_k <= this value + species-default: {'C2H6O[elyt]': 0.9007, 'Li+[elyt]': 1.E-12, 'TFSI-[elyt]': 0.049610, 'O2(e)': 1.E-18} # Replacement mole fractions if elyte composition goes to NaN. + equilibrate: + enable: False # Start with a hold at open circuit? (True or False) + time: 10 # If true, how long is the hold, s + outputs: + show-plots: True # Show the plots and save them (True), or just save + # them (False)? + save-name: '0p2_liu' # Folder label for output files. + validation: 'Liu_O2/0.2mAhcm2.xlsx' #Location of validation data + thickness: 108.E-9 + # - type: 'CC_cycle' # Constant current cycling + # # Specify only one of i_ext or C-rate. The other should be null: + # i_ext: 0.5 mA/cm2 #0.01 A/cm2 # input current density, format: XX units. Units are 'current/area', with no carat for exponents (e.g. A/cm2) + # C-rate: null # input C-rate + # n_cycles: 0.5 # Number of cycles. + # first-step: 'discharge' # Start with charge or discharge? + # phi-cutoff-lower: 1.75 # Simulation cuts off if E_Cell <= this value + # phi-cutoff-upper: 4.8 # Simulation cuts off if E_Cell >= this value + # # species-cutoff: 0.000565 # Simulation cuts off if C_k <= this value + # species-default: {'C2H6O[elyt]': 0.9007, 'Li+[elyt]': 1.E-12, 'TFSI-[elyt]': 0.049610, 'O2(e)': 1.E-18} # Replacement mole fractions if elyte composition goes to NaN. + # equilibrate: + # enable: False # Start with a hold at open circuit? (True or False) + # time: 5 # If true, how long is the hold, s + # outputs: + # show-plots: True # Show the plots and save them (True), or just save + # # them (False)? + # save-name: '0p5_liu' # Folder label for output files. + # validation: 'Liu_O2/0.5mAhcm2.xlsx' #Location of validation data + # thickness: 70.e-9 + # - type: 'CC_cycle' # Constant current cycling + # # Specify only one of i_ext or C-rate. The other should be null: + # i_ext: 1.0 mA/cm2 #0.01 A/cm2 # input current density, format: XX units. Units are 'current/area', with no carat for exponents (e.g. A/cm2) + # C-rate: null # input C-rate + # n_cycles: 0.5 # Number of cycles. + # first-step: 'discharge' # Start with charge or discharge? + # phi-cutoff-lower: 1.75 # Simulation cuts off if E_Cell <= this value + # phi-cutoff-upper: 4.8 # Simulation cuts off if E_Cell >= this value + # # species-cutoff: 0.000565 # Simulation cuts off if C_k <= this value + # species-default: {'C2H6O[elyt]': 0.9007, 'Li+[elyt]': 1.E-12, 'TFSI-[elyt]': 0.049610, 'O2(e)': 1.E-18} # Replacement mole fractions if elyte composition goes to NaN. + # equilibrate: + # enable: False # Start with a hold at open circuit? (True or False) + # time: 5 # If true, how long is the hold, s + # outputs: + # show-plots: True # Show the plots and save them (True), or just save + # # them (False)? + # save-name: '1_liu' # Folder label for output files. + # validation: 'Liu_O2/1mAhcm2.xlsx' #Location of validation data + # thickness: 109.e-9 + +fit-parameters: + - type: 'elyte-transport' + species: 'O2(e)' + guess-value: 24.30E-10 + min: 5.00e-10 + max: 20.30e-10 + - type: 'elyte-transport' + species: 'Li+[elyt]' + guess-value: 6.42e-11 + min: 3.32e-10 + max: 1.4e-9 + # - type: 'cathode-kinetics' + # reaction-index: 0 + # guess-value: 0.46210493 # This is a multiplier to scale the value in the mech. + # min: 0.1 + # max: 10 + # - type: 'cathode-microstructure' + # parameter: 'th_product' + # guess-value: 100.e-9 + # min: 60.e-9 + # max: 105.e-9 + + + +# Cantera inputs: +description: |- + Cantera input file for an Li-Oxygen battery + + This file includes a full set of thermodynamic and kinetic parameters of a + lithium-oxygen battery, in particular: + - Active materials: Li2O2 and Li (li metal) + - Organic electrolyte: TEGDME with 1M LiPF6 + - Interfaces: LCO/electrolyte and Li/electrolyte + - Cathode: The cathode is comprised of catalyst, carbon, (will binder be a later addition.) + - Charge-transfer reactions at the two interfaces + + Bulk phases + =========== + Lithium (anode) + + Ether based electrolyte (electrolyte) + Solvent: Tetraethylene glycol dimethyl ether + Salt: 0.25M LiTFSI + + Interface phases + ================ + + lithium/electrolyte interface (lithium_electrolyte) + Species and site density are dummy entries (as we do not consider surface- + adsorbed species) + + LCO/electrolyte interface (cathode_electrolyte) + Species and site density are dummy entries (as we do not consider surface- + adsorbed species) + + Builds on prior code development (LeBar, 2019) +generator: cti2yaml +cantera-version: 2.5.0 +date: Wed, 11 Dec 2019 16:59:07 -0500 +input-files: [lithium_ion_battery.cti] + +# ================================== +# PHASE DEFINITIONS +# ================================== +phases: + +- name: lithium_metal + thermo: ideal-condensed + species: ['Li[metal]'] + density: 0.534 g/cm^3 #verified + state: + T: 300.0 + P: 1 atm + +- name: electron + thermo: electron-cloud + elements: [E] + species: [electron] + state: + X: {electron: 1.0} + density: 1.0 kg/m^3 + +- name: lithium_electrolyte + thermo: ideal-surface + species: [(dummy)] + kinetics: surface + reactions: [lithium_electrolyte-reactions] + site-density: 0.01 mol/cm^2 #Question: should I change this?? + +- name: electrolyte + thermo: ideal-condensed + elements: [Li, F, C, H, O, E, N, S] + species: ['C2H6O[elyt]', 'Li+[elyt]', 'TFSI-[elyt]', 'O2(e)'] #, 'Li2O2[elyt]'] + state: + X: {'C2H6O[elyt]': 0.8007, 'Li+[elyt]': 0.09350, 'TFSI-[elyt]': 0.09350, 'O2(e)': 19.7E-5} + standard-concentration-basis: unity + +- name: air-elyte-interface + thermo: ideal-surface + species: [(dummy)] + kinetics: surface + reactions: [air-electrolyte-reactions] + site-density: 0.01 mol/cm^2 + +- name: air + thermo: ideal-gas + elements: [O, N, Ar] + species: [O2(g), N2(g), AR(g)] + kinetics: gas + reactions: none + transport: mixture-averaged + state: {T: 300.0, P: 1 atm, X: {O2(g): 0.21, N2(g): 0.78, AR(g): 0.01}} + +- name: conductor + thermo: electron-cloud + elements: [E] + species: [electron] + state: + X: {electron: 1.0} + density: 1.0 kg/m^3 + +- name: cathode-surf + thermo: ideal-surface + species: [(dummy)] + kinetics: surface + reactions: [cathode-surf-reactions] + state: + T: 300.0 + P: 1 atm + site-density: 2.50e-5 mol/cm^2 + +- name: li-oxide + composition: {Li: 2, O: 2} + thermo: fixed-stoichiometry + species: ['Li2O2[cathode]'] + density: 2310 kg/m^3 # https://pubchem.ncbi.nlm.nih.gov/compound/Lithium-peroxide-_Li2_O2#section=Physical-Description + state: + T: 300.0 + P: 1 atm + +# ================================== +# SPECIES DEFINITIONS +# ================================== +species: +- name: Li[metal] + composition: {Li: 1} + thermo: + model: constant-cp + h0: 19.50 kJ/mol # Question: unverified - source? + s0: 29.1 J/mol/K # verified NIST webbook + equation-of-state: + model: constant-volume + molar-volume: 12.998 cm^3/mol #verified + +- name: electron + composition: {E: 1} + thermo: + model: constant-cp + h0: 0.0 kJ/mol + s0: 0.0 J/mol/K + note: |- + 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) + +- name: (dummy) + composition: {} + thermo: + model: constant-cp + h0: 0.0 kJ/mol + s0: 0.0 J/mol/K + note: Dummy species (needed for defining the interfaces) + +- name: C2H6O[elyt] + composition: {C: 2, H: 6, O: 1} + thermo: + model: constant-cp + h0: -184.1 kJ/mol # https://webbook.nist.gov/cgi/cbook.cgi?ID=C115106&Mask=7 + s0: 0.0 J/mol/K # unknown + equation-of-state: + model: constant-volume + molar-volume: 106.29931972789115646258503401361 cm^3/mol + note: |- + Dimethyl ether, MW: 46.069 g/mol + Density of tetraethylene glycol dimethyl ether: 735 kg/m3 (used to calculate species molar volume + as molecular weight (MW)/density) + Taken from: https://www.sigmaaldrich.com/US/en/sds/sigma/d8418 + Molar enthalpy and entropy set to zero (dummy entries as this species does + not participate in chemical reactions) + + +- name: Li+[elyt] + composition: {Li: 1, E: -1} + thermo: + model: constant-cp + h0: -278.49 kJ/mol + s0: 13.4 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 5.508297619047619 cm^3/mol + note: |- #Question: why is this using the density of electrolyte + 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) + +- name: TFSI-[elyt] + composition: {C: 2, F: 6, N: 1, O: 4, S: 1, E: 1} + thermo: + model: constant-cp + h0: 0.0 J/mol + s0: 0.0 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 78.92857142857143 cm^3/mol + note: |- #Question: why is this using the density of electrolyte + bis(trifluoromethylsulphonyl)imide ion, MW: 280.147 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) + +- name: Li2O2[elyt] + composition: {Li: 2, O: 2} + thermo: + model: constant-cp + h0: 0.0 J/mol # Question: should I put this here https + s0: 0.0 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 19.344262295081968261 cm^3/mol + note: |- + Li2O2 particle, MW: 38.94 g/mol + Density of electrolyte: 2013 kg/m3 (used to calculate species molar volume + as molecular weight (MW)/density) + https://www.chemspider.com/Chemical-Structure.145811.html + Molar enthalpy and entropy set to zero (dummy entries as this species does + not participate in chemical reactions) + +- name: O2(e) + composition: {O: 2} + thermo: + model: constant-cp + h0: -40.87 kJ/mol + s0: 0. J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 28.0445 cm^3/mol + note: |- + An O2 molecule dissolved into electrolyte. Thermo tuned to get correct solubility at room temperature. + Molar volume estimated from liquid O2 values. https://en.wikipedia.org/wiki/Liquid_oxygen#cite_note-3 + +- name: O2(g) + composition: {O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, + -1063.94356, 3.65767573] + - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, + -1088.45772, 5.45323129] + note: TPIS89 + transport: + model: gas + geometry: linear + well-depth: 107.4 + diameter: 3.458 + polarizability: 1.6 + rotational-relaxation: 3.8 +- name: N2(g) + composition: {N: 2} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.298677, 1.4082404e-03, -3.963222e-06, 5.641515e-09, -2.444854e-12, + -1020.8999, 3.950372] + - [2.92664, 1.4879768e-03, -5.68476e-07, 1.0097038e-10, -6.753351e-15, + -922.7977, 5.980528] + note: '121286' + transport: + model: gas + geometry: linear + well-depth: 97.53 + diameter: 3.621 + polarizability: 1.76 + rotational-relaxation: 4.0 +- name: AR(g) + composition: {Ar: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366] + - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366] + note: '120186' + transport: + model: gas + geometry: atom + well-depth: 136.5 + diameter: 3.33 + +- name: Li2O2[cathode] + composition: {Li: 2, O: 2} + thermo: + model: constant-cp + h0: -640.3 kJ/mol #-634.3 kJ/mol #Haynes, 2016 + s0: 0.0 J/mol/K + equation-of-state: + model: constant-volume + molar-volume: 19.861904761904753514 cm^3/mol + note: |- + h0 value tuned to give correct open circuit potential, relative to li metal anode. + Lithium peroxide, MW: 45.881 g/mol. + Density of L2O2: 2310 kg/m3 + (Haynes, 2016) + +lithium_electrolyte-reactions: +- equation: Li[metal] <=> Li+[elyt] + electron + id: lithium_faradaic_reaction + rate-constant: {A: 6.0e+20, b: 0.0, Ea: 0.0} #{A: 6.0e+09, b: 0.0, Ea: 0.0} + beta: 0.5 + +cathode-surf-reactions: +- equation: 2 Li+[elyt] + O2(e) + 2 electron <=> Li2O2[cathode] # Reaction 1 + id: cathode_reaction + rate-constant: {A: 5.242e+3, b: 0.0, Ea: 58.0 kJ/mol} #This is already scaled by 0.001 from fit params. #need a new reaction rate + exchange-current-density-formulation: true + beta: 0.5 + +air-electrolyte-reactions: +- equation: O2(g) <=> O2(e) # Reaction 2 + id: cathode_gas_reaction + rate-constant: {A: 2.e-5, b: 0.0, Ea: 0 kJ/mol} + +#"""Citations""" +#Chen, J.; et al., ACS Omega, vol. 4, no. 20708-20714, 2019, DOI: 10.1021/acsomega.9b02941. + +#Gittleson, F; et al, Energy Environ. Sci., 2017,10, 1167-1179, DOI: 10.1039/C6EE02915A. + +#Haynes, W.M.; CRC handbook of chemistry and physics, 2016, 97th Edition / Boca Raton, Florida: CRC Press. + +#LeBar, A., Colorado School of Mines, 2019, URL: https://mountainscholar.org/handle/11124/173373 + +#Yin, Y.; et al., J. Phys. Chem. C, vol. 121, no. 36, pp. 19577–19585, 2017, doi: 10.1021/acs.jpcc.7b05224. + +#Official Soundtrack: +#Day 6 - as of 2021 all albums/songs +# - highlight Finale, Congradulations, 예뻤어 +# - the Magnetic Fields, 69 Love Songs, vol III +#Julia Rocka - Blaza \ No newline at end of file diff --git a/separator_models/porous_separator.py b/separator_models/porous_separator.py index 3a6b196..0206e6b 100644 --- a/separator_models/porous_separator.py +++ b/separator_models/porous_separator.py @@ -4,6 +4,7 @@ Class file for porous separator methods """ +#why won't you come to the top for commits? import cantera as ct import numpy as np diff --git a/simulations/CC_cycle.py b/simulations/CC_cycle.py index 60ebb92..5286a25 100644 --- a/simulations/CC_cycle.py +++ b/simulations/CC_cycle.py @@ -37,6 +37,8 @@ def run(SV_0, an, sep, ca, algvars, params, sim): if 'species-default' in sim: params['species-default'] = sim['species-default'] + ca.thickness_temp=sim['thickness'] + # Figure out which steps and at what currents to run the model. This # returns a tuple of 'charge' and 'discharge' steps, and a tuple with a # current for each step. 'equil' is a flag to indicate whether there is an @@ -90,7 +92,17 @@ def terminate_check(t, SV, SVdot, return_val, inputs): # Create an array of currents, one for each time step: i_data = currents[i]*np.ones_like(solution.values.t) cycle_number = int(i+1-equil)*np.ones_like(solution.values.t) - cycle_capacity = 0.1*solution.values.t*abs(i_data)/3600 + cycle_capacity = 0.1*solution.values.t*abs(i_data)/3600 #mAh/cm2 + phi_ptr = ca.SV_offset+int(ca.SVptr['phi_ed'][-1]) + voltage = solution.values.y.T[phi_ptr] + delta_t = solution.values.t[1:] - solution.values.t[:-1] #s + for item, voltage_value in enumerate(voltage): + if item == 0: + energy_density = [0] + power_density = [currents[i]*voltage_value] + else: + energy_density.append((delta_t[item-1]*currents[i]*voltage_value/3600+energy_density[item-1])) #Wh/m2 + power_density.append(3600*energy_density[item]/sum(delta_t[:item])) #W/m2 # Append the current data array to any preexisting data, for output. # If this is the first step, create the output data array. @@ -98,7 +110,7 @@ def terminate_check(t, SV, SVdot, return_val, inputs): # Stack the times, the current at each time step, and the solution # vector at each time step into a single data array. SV = np.vstack((solution.values.t+data_out[0,-1], cycle_number, - i_data, cycle_capacity, solution.values.y.T)) + i_data, cycle_capacity, power_density, energy_density, solution.values.y.T)) data_out = np.hstack((data_out, SV)) # Use SV at the end of the simualtion as the new initial condition: @@ -107,7 +119,7 @@ def terminate_check(t, SV, SVdot, return_val, inputs): # Stack the times, the current at each time step, and the solution # vector at each time step into a single data array. SV = np.vstack((solution.values.t, cycle_number, i_data, - cycle_capacity, solution.values.y.T)) + cycle_capacity, power_density, energy_density, solution.values.y.T)) data_out = SV # Use SV at the end of the simualtion as the new initial condition: @@ -247,14 +259,14 @@ def output(solution, an, sep, ca, params, sim, plot_flag=True, n_plots = 2 + an.n_plots + ca.n_plots + sep.n_plots # There are 4 variables stored before the state variables: (1) time (s), - # (2) cycle number, (3) current density(A/cm2) , and (4) Capacity (mAh/cm2) - SV_offset = 4 + # (2) cycle number, (3) current density(A/cm2) , (4) Capacity (mAh/cm2) and (5) power density and (6) energy density + SV_offset = 6 # Pointer for cell potential: phi_ptr = SV_offset + ca.SV_offset+int(ca.SVptr['phi_ed'][-1]) # Save the solution as a Pandas dataframe: - labels = (['cycle', 'current', 'capacity'] + an.SVnames + sep.SVnames + labels = (['cycle', 'current', 'capacity','power density','energy density'] + an.SVnames + sep.SVnames + ca.SVnames) solution_df = pd.DataFrame(data = solution.T[:,1:], index = solution.T[:,0], @@ -308,7 +320,7 @@ def plot(an, sep, ca, params, sim): # There are 4 variables stored before the state variables: (1) time (s), # (2) cycle number, (3) current density(A/cm2) , and (4) Capacity (mAh/cm2) - SV_offset = 4 + SV_offset = 6 # Pointer for cell potential: phi_ptr = SV_offset + ca.SV_offset+int(ca.SVptr['phi_ed'][-1]) @@ -368,7 +380,7 @@ def plot(an, sep, ca, params, sim): cycle_fig.set_size_inches((4.0,2.0)) - # iterate over cycles: + # iterate over cycles: for i in range(int(solution[1,-1])): cycle = solution_df[solution_df.iloc[:,0] == i+1] # All times are relative to the start of the step: diff --git a/submodels/bandwidth.py b/submodels/bandwidth.py index 13ced6c..d6bf24f 100644 --- a/submodels/bandwidth.py +++ b/submodels/bandwidth.py @@ -31,7 +31,9 @@ def calc_resid(SV): if abs(dF[i]) > 0: if j > i and abs(i - j) > uband: uband = abs(i - j) + # print('j > i', i, j, uband, lband) elif i > j and abs(i - j) > lband: lband = abs(i - j) - + # print('i > j', i, j, uband, lband) + # print(an.SVptr, sep.SVptr, ca.SVptr) return lband, uband