From a43bce43a05b439099fd13f58377149e9d0d9361 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Fri, 11 Feb 2022 13:39:40 -0700 Subject: [PATCH 01/14] Working wrapper file --- SSRcode.py | 107 +++++++++++++++++++++++++++++++++++++++++++++++++++++ wrapper.py | 47 +++++++++++++++++++++++ 2 files changed, 154 insertions(+) create mode 100644 SSRcode.py create mode 100644 wrapper.py diff --git a/SSRcode.py b/SSRcode.py new file mode 100644 index 0000000..198a74d --- /dev/null +++ b/SSRcode.py @@ -0,0 +1,107 @@ + +#%% +import pandas as pd +import os + +#https://doi.org/10.1016/j.electacta.2014.09.074 for refdata +refpath = 'D:\projects\Data\ToProcess' #change this to file directory +path2 = os.path.abspath(os.getcwd())+'/outputs' +modelfile = 'output.csv' +comparison = 'Refdata.xlsx' +#%% + +def interpolate(voltage1, voltage2, capacity1, capacity2, capacity3): + voltage3 = (voltage2 - voltage1)*(capacity3-capacity1)/(capacity2-capacity1) +voltage1 + return voltage3 + +def locater1(refcapacity, refvoltage): + lowerneighbour_ind = RefData[RefData['Capacity'] < refcapacity]['Capacity'].idxmax() + higherneighbour_ind = lowerneighbour_ind + 1 + capacity1 = RefData['Capacity'][lowerneighbour_ind] + capacity2 = RefData['Capacity'][higherneighbour_ind] + capacity3 = refcapacity + voltage1 = RefData['Voltage'][lowerneighbour_ind] + voltage2 = RefData['Voltage'][higherneighbour_ind] + v3 = (interpolate(voltage1, voltage2, capacity1, capacity2, capacity3)) + SSR = (refvoltage - v3)**2 + return SSR + +def SSRmain(folder_name): + SumSR = 0 + CodeData = pd.read_csv(path2+"/" + folder_name + "/" + modelfile) + cyclestart = CodeData[CodeData['cycle'] > 0]['cycle'].idxmin() + 1 + chargestart = CodeData[CodeData['cycle'] > 1]['cycle'].idxmin() + Voltagesort = CodeData.filter(like='phi_ed') + Phi_ed_loc = Voltagesort.columns[-1] + Desiredcolumn = Voltagesort[Phi_ed_loc][cyclestart:chargestart] + Desiredcolumn = Desiredcolumn.reset_index() + for x, y in enumerate(CodeData['capacity'][cyclestart:chargestart]): + refvoltage = Desiredcolumn[Phi_ed_loc][x] + SumSR += locater1(y, refvoltage) + return SumSR + + +def locater2(dataset, refcapacity1, refvoltage1): + new = CodeData['capacity'][cyclestart:chargestart] + lowerneighbour_ind = new[new < refcapacity1].idxmax() + higherneighbour_ind = lowerneighbour_ind + 1 + capacity1 = new[lowerneighbour_ind] + capacity2 = new[higherneighbour_ind] + capacity3 = refcapacity + print(capacity1, capacity2, capacity3) + voltage1 = Desiredcolumn[Phi_ed_loc][lowerneighbour_ind] + voltage2 = Desiredcolumn[Phi_ed_loc][higherneighbour_ind] + v3 = (interpolate(voltage1, voltage2, capacity1, capacity2, capacity3)) + SSR = (refvoltage - v3)**2 + return SSR + +def SSRmain3(folder_name, dataset): + SumSR = 0 + CodeData = pd.read_csv(path2+"/" + folder_name + "/" + "output.csv") + cyclestart = CodeData[CodeData['cycle'] > 0]['cycle'].idxmin() + 1 + chargestart = CodeData[CodeData['cycle'] > 1]['cycle'].idxmin() + Voltagesort = CodeData.filter(like='phi_ed') + Phi_ed_loc = Voltagesort.columns[-1] + Desiredcolumn = Voltagesort[Phi_ed_loc][cyclestart:chargestart] + Desiredcolumn = Desiredcolumn.reset_index() + for x, y in enumerate(RefData[dataset]['Capacity']): + refvoltage1 = RefData[dataset]['Voltage'][x] + SumSR += locater2(dataset, y, refvoltage) + return SumSR + + + + +#%% +ID_key = {"CPCN04": "50CP50CNT0.4", "CP04":"CP0.4", "CPCN05":"50CP50CNT0.5", "CP05": "CP0.5", "CPCN06":"50CP50CNT0.6", "CP06":"CP0.6"} + +#%% +SSRarray = [] +altSSRarray = [] +for folder_name in os.listdir(path2): + if "1045" in folder_name: + string = folder_name + array = string.split("_") + dataset = ID_key[array[0]] + RefData = pd.read_excel(refpath + "/"+ comparison, sheet_name= dataset) + SSR_file = SSRmain(folder_name) + SSRarray.append(SSR_file) + +print(SSRarray) +# print(altSSRarray) + + +#result = [locater("50CP50CNT0.4", x, refcapcity) for x, refcapcity in enumerate(CodeData['Capacity'])] +#print(result) + +# print (Data['50CP50CNT0.4']['Capacity'][row]) +# lowerneighbour_ind = Data2[Data2['Capacity'] < Temp]['Capacity'].idxmax() +# higherneighbour_ind = Data2[Data2['Capacity'] > Temp]['Capacity'].idxmin() +# v3 = interpolate(Data2['Capacity'][higherneighbour_ind], voltage2,Data2['Capacity'][lowerneighbour_ind], Data2['Capacity'][higherneighbour_ind], capacity3) +# SSR += (Voltage_test[i] - v3)**2 + +# %% +# Official Playlist +# Small World - Bump of Chicken +# Taxi Driver - Rkomi +# 30 - UVERworld \ No newline at end of file diff --git a/wrapper.py b/wrapper.py new file mode 100644 index 0000000..ebb7d3b --- /dev/null +++ b/wrapper.py @@ -0,0 +1,47 @@ +# wrapper.py +#%% +path = os.path.abspath(os.getcwd())+'/outputs' +import sys +import os +import pandas as pd +from bat_can import bat_can +from SSRcode import SSRmain + +#%% + + +foldarray =os.listdir(path) +inputfile = "LiMetal_PorousSep_Air" +bat_can(inputfile) +foldarray2 = os.listdir(path) + +#%% +ID_key = {"CPCN04": "50CP50CNT0.4", "CP04":"CP0.4", "CPCN05":"50CP50CNT0.5", "CP05": "CP0.5", "CPCN06":"50CP50CNT0.6", "CP06":"CP0.6"} +comparison = 'Refdata.xlsx' +SSRarray = [] +refpath = 'D:\projects\Data\ToProcess' #change this to file directory +folder_name = list(set(foldarray).symmetric_difference(set(foldarray2)))[0] +array = folder.split("_") +dataset = ID_key[array[0]] +RefData = pd.read_excel(refpath + "/"+ comparison, sheet_name= dataset) +SSR_file = SSRmain(folder_name) +SSRarray.append(SSR_file) +print(SSRarray) + + + +#%% +# SSRarray = [] +# path = os.path.abspath(os.getcwd())+'/outputs' +# refpath = 'D:\projects\Data\ToProcess' #change this to file directory +# for folder_name in os.listdir(path): +# if "1045" in folder_name: +# string = folder_name +# array = string.split("_") +# dataset = ID_key[array[0]] +# RefData = pd.read_excel(refpath + "/"+ comparison, sheet_name= dataset) +# SSR_file = SSRmain(folder_name) +# SSRarray.append(SSR_file) + + +#%% From efc94c19c6b891c67fbab2647521b495ffd385d4 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Fri, 11 Feb 2022 13:45:47 -0700 Subject: [PATCH 02/14] Fixed branch --- SSRcode.py | 107 ----------------------------------------------------- wrapper.py | 47 ----------------------- 2 files changed, 154 deletions(-) delete mode 100644 SSRcode.py delete mode 100644 wrapper.py diff --git a/SSRcode.py b/SSRcode.py deleted file mode 100644 index 198a74d..0000000 --- a/SSRcode.py +++ /dev/null @@ -1,107 +0,0 @@ - -#%% -import pandas as pd -import os - -#https://doi.org/10.1016/j.electacta.2014.09.074 for refdata -refpath = 'D:\projects\Data\ToProcess' #change this to file directory -path2 = os.path.abspath(os.getcwd())+'/outputs' -modelfile = 'output.csv' -comparison = 'Refdata.xlsx' -#%% - -def interpolate(voltage1, voltage2, capacity1, capacity2, capacity3): - voltage3 = (voltage2 - voltage1)*(capacity3-capacity1)/(capacity2-capacity1) +voltage1 - return voltage3 - -def locater1(refcapacity, refvoltage): - lowerneighbour_ind = RefData[RefData['Capacity'] < refcapacity]['Capacity'].idxmax() - higherneighbour_ind = lowerneighbour_ind + 1 - capacity1 = RefData['Capacity'][lowerneighbour_ind] - capacity2 = RefData['Capacity'][higherneighbour_ind] - capacity3 = refcapacity - voltage1 = RefData['Voltage'][lowerneighbour_ind] - voltage2 = RefData['Voltage'][higherneighbour_ind] - v3 = (interpolate(voltage1, voltage2, capacity1, capacity2, capacity3)) - SSR = (refvoltage - v3)**2 - return SSR - -def SSRmain(folder_name): - SumSR = 0 - CodeData = pd.read_csv(path2+"/" + folder_name + "/" + modelfile) - cyclestart = CodeData[CodeData['cycle'] > 0]['cycle'].idxmin() + 1 - chargestart = CodeData[CodeData['cycle'] > 1]['cycle'].idxmin() - Voltagesort = CodeData.filter(like='phi_ed') - Phi_ed_loc = Voltagesort.columns[-1] - Desiredcolumn = Voltagesort[Phi_ed_loc][cyclestart:chargestart] - Desiredcolumn = Desiredcolumn.reset_index() - for x, y in enumerate(CodeData['capacity'][cyclestart:chargestart]): - refvoltage = Desiredcolumn[Phi_ed_loc][x] - SumSR += locater1(y, refvoltage) - return SumSR - - -def locater2(dataset, refcapacity1, refvoltage1): - new = CodeData['capacity'][cyclestart:chargestart] - lowerneighbour_ind = new[new < refcapacity1].idxmax() - higherneighbour_ind = lowerneighbour_ind + 1 - capacity1 = new[lowerneighbour_ind] - capacity2 = new[higherneighbour_ind] - capacity3 = refcapacity - print(capacity1, capacity2, capacity3) - voltage1 = Desiredcolumn[Phi_ed_loc][lowerneighbour_ind] - voltage2 = Desiredcolumn[Phi_ed_loc][higherneighbour_ind] - v3 = (interpolate(voltage1, voltage2, capacity1, capacity2, capacity3)) - SSR = (refvoltage - v3)**2 - return SSR - -def SSRmain3(folder_name, dataset): - SumSR = 0 - CodeData = pd.read_csv(path2+"/" + folder_name + "/" + "output.csv") - cyclestart = CodeData[CodeData['cycle'] > 0]['cycle'].idxmin() + 1 - chargestart = CodeData[CodeData['cycle'] > 1]['cycle'].idxmin() - Voltagesort = CodeData.filter(like='phi_ed') - Phi_ed_loc = Voltagesort.columns[-1] - Desiredcolumn = Voltagesort[Phi_ed_loc][cyclestart:chargestart] - Desiredcolumn = Desiredcolumn.reset_index() - for x, y in enumerate(RefData[dataset]['Capacity']): - refvoltage1 = RefData[dataset]['Voltage'][x] - SumSR += locater2(dataset, y, refvoltage) - return SumSR - - - - -#%% -ID_key = {"CPCN04": "50CP50CNT0.4", "CP04":"CP0.4", "CPCN05":"50CP50CNT0.5", "CP05": "CP0.5", "CPCN06":"50CP50CNT0.6", "CP06":"CP0.6"} - -#%% -SSRarray = [] -altSSRarray = [] -for folder_name in os.listdir(path2): - if "1045" in folder_name: - string = folder_name - array = string.split("_") - dataset = ID_key[array[0]] - RefData = pd.read_excel(refpath + "/"+ comparison, sheet_name= dataset) - SSR_file = SSRmain(folder_name) - SSRarray.append(SSR_file) - -print(SSRarray) -# print(altSSRarray) - - -#result = [locater("50CP50CNT0.4", x, refcapcity) for x, refcapcity in enumerate(CodeData['Capacity'])] -#print(result) - -# print (Data['50CP50CNT0.4']['Capacity'][row]) -# lowerneighbour_ind = Data2[Data2['Capacity'] < Temp]['Capacity'].idxmax() -# higherneighbour_ind = Data2[Data2['Capacity'] > Temp]['Capacity'].idxmin() -# v3 = interpolate(Data2['Capacity'][higherneighbour_ind], voltage2,Data2['Capacity'][lowerneighbour_ind], Data2['Capacity'][higherneighbour_ind], capacity3) -# SSR += (Voltage_test[i] - v3)**2 - -# %% -# Official Playlist -# Small World - Bump of Chicken -# Taxi Driver - Rkomi -# 30 - UVERworld \ No newline at end of file diff --git a/wrapper.py b/wrapper.py deleted file mode 100644 index ebb7d3b..0000000 --- a/wrapper.py +++ /dev/null @@ -1,47 +0,0 @@ -# wrapper.py -#%% -path = os.path.abspath(os.getcwd())+'/outputs' -import sys -import os -import pandas as pd -from bat_can import bat_can -from SSRcode import SSRmain - -#%% - - -foldarray =os.listdir(path) -inputfile = "LiMetal_PorousSep_Air" -bat_can(inputfile) -foldarray2 = os.listdir(path) - -#%% -ID_key = {"CPCN04": "50CP50CNT0.4", "CP04":"CP0.4", "CPCN05":"50CP50CNT0.5", "CP05": "CP0.5", "CPCN06":"50CP50CNT0.6", "CP06":"CP0.6"} -comparison = 'Refdata.xlsx' -SSRarray = [] -refpath = 'D:\projects\Data\ToProcess' #change this to file directory -folder_name = list(set(foldarray).symmetric_difference(set(foldarray2)))[0] -array = folder.split("_") -dataset = ID_key[array[0]] -RefData = pd.read_excel(refpath + "/"+ comparison, sheet_name= dataset) -SSR_file = SSRmain(folder_name) -SSRarray.append(SSR_file) -print(SSRarray) - - - -#%% -# SSRarray = [] -# path = os.path.abspath(os.getcwd())+'/outputs' -# refpath = 'D:\projects\Data\ToProcess' #change this to file directory -# for folder_name in os.listdir(path): -# if "1045" in folder_name: -# string = folder_name -# array = string.split("_") -# dataset = ID_key[array[0]] -# RefData = pd.read_excel(refpath + "/"+ comparison, sheet_name= dataset) -# SSR_file = SSRmain(folder_name) -# SSRarray.append(SSR_file) - - -#%% From 693603992db5ba8a35f0ddc6842c2ba3327e16b8 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Thu, 2 Feb 2023 13:04:33 -0700 Subject: [PATCH 03/14] Updated file for current batcan file --- inputs/LiMetal_PorousSep_Air.yaml | 89 ++++++++++++++++--------------- 1 file changed, 47 insertions(+), 42 deletions(-) diff --git a/inputs/LiMetal_PorousSep_Air.yaml b/inputs/LiMetal_PorousSep_Air.yaml index b4d30d3..8cd9ba6 100644 --- a/inputs/LiMetal_PorousSep_Air.yaml +++ b/inputs/LiMetal_PorousSep_Air.yaml @@ -30,7 +30,11 @@ cell-description: sigma_io: 0.272 # S/m for 1M LiTFSI in TEGDME at 25 deg C // Chen, 2019 phi_0: 2.91 # 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 @@ -42,9 +46,9 @@ cell-description: - species: 'ClO4-[elyt]' D_k: 8.79e-11 # m2 S-1 // Saito, 2017 - species: 'O2(e)' - D_k: 4.26e-11 # m2 s-1 // Gittleson, 2017 - - species: 'Li2O2[elyt]' - D_k: 2e-13 # can't find (assume low?) + D_k: 2.5283e-8 # m2 s-1 // Gittleson, 2017 + # - species: 'Li2O2[elyt]' + # D_k: 2e-13 # can't find (assume low?) cathode: class: 'air_electrode' n-points: 5 @@ -59,7 +63,7 @@ cell-description: thickness: 3.50e-4 # cathode thickness, m r_host: 6e-6 # Carbon Particle radius, m d_oxide: 1.0e-6 #Oxide diameter, m - th_oxide: 500e-9 #oxide thickness, m + th_product: 500e-9 #oxide thickness, m phi_0: 2.95 # Initial electric potential, V // Yin, 2017 sigma_el: 1.0E1 #electrical ondutivity of the host eps_host: 0.621138 # Solid phase volume fraction of carbon - @@ -75,8 +79,8 @@ cell-description: # Simulation parameters: parameters: - T: 23 # C - P: 101325 # Pa + T: 23 C + P: 101325 Pa # Describe simulation type, parameters, etc. simulations: - type: 'CC_cycle' # Constant current cycling @@ -94,39 +98,39 @@ parameters: outputs: show-plots: True # Show the plots and save them (True), or just save # them (False)? - savename: 'CPCN04' # Folder label for output files. - - type: 'CC_cycle' # Constant current cycling - # Specify only one of i_ext or C-rate. The other should be null: - i_ext: 0.0005 A/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: 1 # 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: 1e-12 # 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)? - savename: 'CPCN05' # Folder label for output files. - - type: 'CC_cycle' # Constant current cycling - # Specify only one of i_ext or C-rate. The other should be null: - i_ext: 0.0006 A/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: 1 # 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: 1e-12 # 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)? - savename: 'CPCN06' # Folder label for output files. + save-name: 'CPCN04' # Folder label for output files. + # - type: 'CC_cycle' # Constant current cycling + # # Specify only one of i_ext or C-rate. The other should be null: + # i_ext: 0.0005 A/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: 1 # 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: 1e-12 # 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: 'CPCN05' # Folder label for output files. + # - type: 'CC_cycle' # Constant current cycling + # # Specify only one of i_ext or C-rate. The other should be null: + # i_ext: 0.0006 A/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: 1 # 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: 1e-12 # 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: 'CPCN06' # Folder label for output files. # Cantera inputs: description: |- @@ -193,9 +197,9 @@ phases: - name: electrolyte thermo: ideal-condensed elements: [Li, Cl, C, H, O, E] - species: ['C10H22O5[elyt]', 'Li+[elyt]', 'ClO4-[elyt]', 'O2(e)', 'Li2O2[elyt]'] + species: ['C10H22O5[elyt]', 'Li+[elyt]', 'ClO4-[elyt]', 'O2(e)'] state: - X: {'C10H22O5[elyt]': 0.9007, 'Li+[elyt]': 0.049610, 'ClO4-[elyt]': 0.049610, 'O2(e)': 1.39E-4, 'Li2O2[elyt]': 1E-9} + X: {'C10H22O5[elyt]': 0.9007, 'Li+[elyt]': 0.049610, 'ClO4-[elyt]': 0.049610, 'O2(e)': 6.39E-5} #, 'Li2O2[elyt]': 1E-9} standard-concentration-basis: unity - name: air-elyte-interface @@ -416,4 +420,5 @@ air-electrolyte-reactions: #Official Soundtrack: #Day 6 - as of 2021 all albums/songs - # - highlight Finale, Congradulations, 예뻤어 \ No newline at end of file + # - highlight Finale, Congradulations, 예뻤어 + #Deathcab for Cutie - Pepper \ No newline at end of file From c7e42bf808acf7f714982a4a5db94bf7d300f612 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Thu, 2 Feb 2023 13:05:54 -0700 Subject: [PATCH 04/14] Units uncommented --- inputs/LiMetal_PorousSep_spmLCO_input.yaml | 12 ++++++------ inputs/LiMetal_Resistor_spmLCO_input.yaml | 8 ++++---- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/inputs/LiMetal_PorousSep_spmLCO_input.yaml b/inputs/LiMetal_PorousSep_spmLCO_input.yaml index ba11fa6..c384acf 100644 --- a/inputs/LiMetal_PorousSep_spmLCO_input.yaml +++ b/inputs/LiMetal_PorousSep_spmLCO_input.yaml @@ -16,7 +16,7 @@ cell-description: charge: 1 mobile-ion: 'Li+[elyt]' # Species name for Li+ in elyte. thickness: 8e-6 # Initial anode thickness, m - minimum-thickness: 1e-9 # minimum anode thickness, m + minimum-thickness: 1e-9 # minimum anode thickness, m phi_0: 0. # Initial electric potential, V A_surf_ratio: 1.0 # Interface area per unit geometric area C_dl: 9e-6 #F/m2 @@ -30,11 +30,11 @@ cell-description: transport: mobile-ion: 'Li+[elyt]' # Species name for Li+ in elyte. model: 'dilute-solution' - diffusion-coefficients: # Species names must match those in the phase + diffusion-coefficients: # Species names must match those in the phase #definition, below: - species: 'C3H4O3[elyt]' D_k: 1e-15 - - species: 'C4H6O3[elyt]' + - species: 'C4H6O3[elyt]' D_k: 1e-15 - species: 'Li+[elyt]' D_k: 2e-11 @@ -61,14 +61,14 @@ cell-description: # Simulation parameters: parameters: - T: 60 # C - P: 101325 # Pa + T: 60 C + P: 101325 Pa # Describe what to do with the results: outputs: show-plots: True # Do you want the plots shown and saved, or just saved? # Describe simulation type, parameters, etc. simulation: - + simulations: - type: 'CC_cycle' # Constant current cycling # Specify only one of i_ext or C-rate. The other should be null: diff --git a/inputs/LiMetal_Resistor_spmLCO_input.yaml b/inputs/LiMetal_Resistor_spmLCO_input.yaml index 834ec64..6da7023 100644 --- a/inputs/LiMetal_Resistor_spmLCO_input.yaml +++ b/inputs/LiMetal_Resistor_spmLCO_input.yaml @@ -16,7 +16,7 @@ cell-description: charge: 1 mobile-ion: 'Li+[elyt]' # Species name for Li+ in elyte. thickness: 15e-6 # anode thickness, m - minimum-thickness: 1e-9 # minimum anode thickness. + 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-3 #F/m2 @@ -48,8 +48,8 @@ cell-description: # Simulation parameters: parameters: - T: 60 # C - P: 101325 # Pa + T: 60 C + P: 101325 Pa # Describe what to do with the results: outputs: show-plots: True # Do you want the plots shown and saved, or just saved? @@ -61,7 +61,7 @@ parameters: C-rate: 0.01 # input C-rate n_cycles: 0.01 # Number of cycles. Currently must be 0.5 or an int. first-step: 'discharge' # Start with charge or discharge? - equilibrate: + equilibrate: enable: True # Begin with a hold at i_ext = 0? time: 100 # If true, how long is the hold, s phi-cutoff-lower: 2.25 # Simulation cuts off if E_Cell <= this value From 0a3f8fbcc49b12c2b0daf897ef0b942cf27d8c01 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Thu, 2 Feb 2023 13:06:10 -0700 Subject: [PATCH 05/14] corrected indexing --- separator_models/porous_separator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/separator_models/porous_separator.py b/separator_models/porous_separator.py index 93f1c6d..3a6b196 100644 --- a/separator_models/porous_separator.py +++ b/separator_models/porous_separator.py @@ -269,7 +269,7 @@ def output(self, axs, solution, an, ca, SV_offset, ax_offset): # Axis 5: Li+ concentration: Ck_elyte_an = solution[an.SVptr['C_k_elyte'][0]+SV_offset,:] - axs[ax_offset+1].plot(solution[0,:]/3600, + axs[ax_offset+1].plot(solution[0,:]/3600, Ck_elyte_an[an.index_Li_elyte,:], label="an interface") Ck_elyte_sep_ptr = \ @@ -283,7 +283,7 @@ def output(self, axs, solution, an, ca, SV_offset, ax_offset): Ck_elyte_ca = \ solution[ca.SV_offset+ca.SVptr['C_k_elyte'][j]+SV_offset,:] axs[ax_offset+1].plot(solution[0,:]/3600, - Ck_elyte_ca[ca.index_Li,:]) + Ck_elyte_ca[ca.index_Li_elyte,:]) axs[ax_offset+1].set_ylabel('Li+ concentration \n(kmol/m$^3$') From 33a35d9b6c3bf6ed00625411e2cd679b1a15300d Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Fri, 10 Mar 2023 11:04:43 -0700 Subject: [PATCH 06/14] implement graded cathode --- electrode_models/air_electrode.py | 234 +++++++++++++++--------------- 1 file changed, 120 insertions(+), 114 deletions(-) diff --git a/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index ad09bc7..03a4d1a 100644 --- a/electrode_models/air_electrode.py +++ b/electrode_models/air_electrode.py @@ -8,36 +8,36 @@ from math import tanh import numpy as np -class electrode(): +class electrode(): """ Create an electrode object representing the air 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.gas_obj = ct.Solution(input_file, inputs['gas-phase']) self.elyte_obj = ct.Solution(input_file, inputs['electrolyte-phase']) - self.gas_elyte_obj = ct.Interface(input_file, inputs['elyte-iphase'], + self.gas_elyte_obj = ct.Interface(input_file, inputs['elyte-iphase'], [self.gas_obj, self.elyte_obj]) self.host_obj = ct.Solution(input_file, inputs['host-phase']) self.product_obj = ct.Solution(input_file, inputs['product-phase']) - self.surf_obj = ct.Interface(input_file, inputs['surf-iphase'], + self.surf_obj = ct.Interface(input_file, inputs['surf-iphase'], [self.product_obj, self.elyte_obj, self.host_obj]) # Electrode thickness and inverse thickness: - self.n_points = inputs['n-points'] + self.n_points = inputs['n-points'] self.dy = inputs['thickness']/self.n_points self.dyInv = 1/self.dy - # Anode or cathode? Positive external current delivers positive charge - # to the anode, and removes positive charge from the cathode. For the - # cathode, the first node is at the separator (j=0), whereas the last + # Anode or cathode? Positive external current delivers positive charge + # to the anode, and removes positive charge from the cathode. For the + # cathode, the first node is at the separator (j=0), whereas the last # node has j=0. self.name = electrode_name if self.name=='anode': @@ -49,17 +49,24 @@ 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']) # Phase volume fractions - self.eps_host = inputs['eps_host'] - self.eps_product_init =inputs['eps_product'] - self.eps_elyte_init = 1 - self.eps_host - self.eps_product_init + eps_host = inputs['eps_host'] + if len(eps_host) == 1: + self.eps_host = np.repeat(eps_host[0], self.n_points) + elif len(eps_host) == self.n_points: + self.eps_host = np.asarray(eps_host) + else: + raise ValueError("Porosity must be either a single scalar, or must match the discretization of the air electrode.") + + self.eps_product_init = inputs['eps_product'] + self.eps_elyte_init = 1. - self.eps_host - self.eps_product_init self.sigma_el =inputs['sigma_el'] - # The following calculations assume spherical particles of a single + # The following calculations assume spherical particles of a single # radius, with no overlap. self.r_host = inputs['r_host'] self.th_product = inputs['th_product'] @@ -69,45 +76,45 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, if 'host_form' not in inputs or inputs['host_form'] == 'sphere': self.A_init = self.eps_host * 3./self.r_host elif inputs['host_form'] == 'cylinder': - self.A_init = self.eps_host * 2./self.r_host + self.A_init = self.eps_host * 2./self.r_host else: raise ValueError("Support host_form values include: 'sphere'", " or 'cylinder'.") if 'host_surf_frac' in inputs: self.A_init *= inputs['host_surf_frac'] - # For some models, the elyte thickness is different from that of the + # For some models, the elyte thickness is different from that of the # electrode, so we specify it separately: self.dy_elyte = self.dy # Inverse double layer capacitance, per unit interfacial area. self.C_dl_Inv = 1/inputs['C_dl'] - # 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_init**1.5 # where would we use this? - - # SV_offset specifies the index of the first SV variable for the + self.elyte_microstructure = self.eps_elyte_init[0]**1.5 #This is used by 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 # Determine the electrode capacity (Ah/m2) - # Max voume concentration of the product species (assuming all + # Max voume concentration of the product species (assuming all # electrolyte has been replaced by product species) stored_species = inputs['stored-species'] v_molar_prod = \ self.product_obj[stored_species['name']].partial_molar_volumes[0] self.capacity = (stored_species['charge']*ct.faraday - * self.eps_elyte_init * inputs['thickness'] + * self.eps_elyte_init[0] * inputs['thickness'] / (3600 * v_molar_prod)) - - # Minimum volume fraction for the product phase, below which product + + # Minimum volume fraction for the product phase, below which product # phase consumption reaction shut off: self.product_phase_min = inputs['product-phase-min'] - # Number of state variables: electrode potential, double layer - # potential, electrolyte composition, product phase volume fraction + # Number of state variables: electrode potential, double layer + # potential, electrolyte composition, product phase volume fraction self.n_vars = 3 + self.elyte_obj.n_species self.n_vars_tot = self.n_points*self.n_vars @@ -134,18 +141,18 @@ def initialize(self, inputs, sep_inputs): self.SVptr['phi_ed'] = np.arange(0, self.n_vars_tot, self.n_vars) self.SVptr['phi_dl'] = np.arange(1, self.n_vars_tot, self.n_vars) self.SVptr['eps_product'] = np.arange(2, self.n_vars_tot, self.n_vars) - self.SVptr['C_k_elyte'] = np.ndarray(shape=(self.n_points, - self.elyte_obj.n_species), dtype='int') + self.SVptr['C_k_elyte'] = np.ndarray(shape=(self.n_points, + self.elyte_obj.n_species), dtype='int') for i in range(self.n_points): - self.SVptr['C_k_elyte'][i,:] = range(3 + i*self.n_vars, + self.SVptr['C_k_elyte'][i,:] = range(3 + i*self.n_vars, 3 + i*self.n_vars + self.elyte_obj.n_species) - self.SVnames = (['phi_ed', 'phi_dl', 'eps_product'] + self.SVnames = (['phi_ed', 'phi_dl', 'eps_product'] + self.elyte_obj.species_names[:])*self.n_points - - # 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_tot) # Save the indices of any algebraic variables: @@ -159,7 +166,7 @@ def initialize(self, inputs, sep_inputs): SV[self.SVptr['phi_dl']] = sep_inputs['phi_0'] - inputs['phi_0'] #V SV[self.SVptr['eps_product']] = self.eps_product_init #Volume Fraction SV[self.SVptr['C_k_elyte']] = self.elyte_obj.concentrations # kmol/m3 - + return SV def residual(self, t, SV, SVdot, sep, counter, params): @@ -170,15 +177,15 @@ def residual(self, t, SV, SVdot, sep, counter, params): 1. The electric potential in the electrode phase 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 invariant (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 invariant (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) @@ -196,12 +203,12 @@ def residual(self, t, SV, SVdot, sep, counter, params): # Initialize the residual: resid = np.zeros((self.n_vars_tot,)) - # Save local copies of the solution vectors, pointers for this + # Save local copies of the solution vectors, pointers for this # electrode: SVptr = self.SVptr SV_loc = SV[SVptr['electrode']] SVdot_loc = SVdot[SVptr['electrode']] - + # Start at the separator boundary: j = self.nodes[0] @@ -210,18 +217,18 @@ def residual(self, t, SV, SVdot, sep, counter, params): phi_elyte = phi_ed + SV_loc[SVptr['phi_dl'][j]] 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 - - # 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 + eps_elyte = 1. - eps_product - self.eps_host[j] + + # 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 # correct sign. - N_k_in, i_io_in = (self.i_ext_flag*X for X in + N_k_in, i_io_in = (self.i_ext_flag*X for X in sep.electrode_boundary_flux(SV, self, params['T'])) - # No electronic current at the separator boundary: + # No electronic current at the separator boundary: i_el_in = 0 - + # Loop through the finite volumes: for j_next in self.nodes[1:]: # Set Cantera object properties: @@ -232,52 +239,52 @@ def residual(self, t, SV, SVdot, sep, counter, params): except: X_elyte = params['species-default'] self.elyte_obj.X = X_elyte - + # Set microstructure multiplier for effective diffusivities #TODO #48 self.elyte_microstructure = eps_elyte**1.5 - + # Read out state variables for 'next' node toward air boundary: phi_ed_next = SV_loc[SVptr['phi_ed'][j_next]] phi_elyte_next = phi_ed_next + SV_loc[SVptr['phi_dl'][j_next]] c_k_elyte_next = SV_loc[SVptr['C_k_elyte'][j_next]] eps_product_next = SV_loc[SVptr['eps_product'][j_next]] - eps_elyte_next = 1. - eps_product_next- self.eps_host - - # Load the node properties into dict structures + eps_elyte_next = 1. - eps_product_next- self.eps_host[j_next] + + # Load the node properties into dict structures # (1 = local, 2 = next node) - state_1 = {'C_k': c_k_elyte, 'phi':phi_elyte, + state_1 = {'C_k': c_k_elyte, 'phi':phi_elyte, 'T': params['T'], 'dy':self.dy, 'microstructure':eps_elyte**1.5} - state_2 = {'C_k': c_k_elyte_next, 'phi':phi_elyte_next, - 'T': params['T'], 'dy':self.dy, + state_2 = {'C_k': c_k_elyte_next, 'phi':phi_elyte_next, + 'T': params['T'], 'dy':self.dy, 'microstructure':eps_elyte_next**1.5} - # Calculate fluxes and currents out of this node, into the next + # Calculate fluxes and currents out of this node, into the next # node toward the current collector: N_k_out, i_io_out = sep.elyte_transport(state_1, state_2, sep) i_el_out = (self.sigma_el*(phi_ed - phi_ed_next)*self.dyInv) # Total current (ionic plus electronic) must be spatially invariant. 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 - eps_product/self.th_product + A_avail = self.A_init[j] - eps_product/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 - # drops below the user-specified minimum, any reactions that + # Multiplier to scale phase destruction rates. As eps_product + # drops below the user-specified minimum, any reactions that # consume the phase have their rates quickly go to zero: mult = tanh(eps_product / self.product_phase_min) - # Chemical production rate of the product phase: + # Chemical production rate of the product phase: # (kmol/m2-interface/s) sdot_product = (self.surf_obj.get_creation_rates(self.product_obj) - mult * self.surf_obj.get_destruction_rates(self.product_obj)) - + # Rate of change of the product phase volume fraction: resid[SVptr['eps_product'][j]] = \ - (SVdot_loc[SVptr['eps_product'][j]] - A_avail + (SVdot_loc[SVptr['eps_product'][j]] - A_avail * np.dot(sdot_product, self.product_obj.partial_molar_volumes)) # Production rate of the electron (moles / m2 interface / s) @@ -285,13 +292,13 @@ def residual(self, t, SV, SVdot, sep, counter, params): (mult * self.surf_obj.get_creation_rates(self.host_obj) - self.surf_obj.get_destruction_rates(self.host_obj)) - # Positive Faradaic current corresponds to positive charge created + # Positive Faradaic current corresponds to positive charge created # in the electrode: i_Far = -(ct.faraday * sdot_electron) # Double layer current has the same sign as i_Far i_dl = self.i_ext_flag*(i_io_in - i_io_out)/A_surf_ratio - i_Far - resid[SVptr['phi_dl'][j]] = (SVdot_loc[SVptr['phi_dl'][j]] - + resid[SVptr['phi_dl'][j]] = (SVdot_loc[SVptr['phi_dl'][j]] - i_dl*self.C_dl_Inv) # Species production rate for electrolyte species: @@ -301,14 +308,14 @@ def residual(self, t, SV, SVdot, sep, counter, params): # The double layer current acts as an additional chemical source/ # sink term: - sdot_elyte_host[self.index_Li_elyte] -= i_dl / ct.faraday + sdot_elyte_host[self.index_Li_elyte] -= i_dl / ct.faraday # Change in electrolyte species concentration, per unit time: - resid[SVptr['C_k_elyte'][j]] = (SVdot_loc[SVptr['C_k_elyte'][j]] - - (N_k_in - N_k_out+ sdot_elyte_host * A_surf_ratio) + resid[SVptr['C_k_elyte'][j]] = (SVdot_loc[SVptr['C_k_elyte'][j]] + - (N_k_in - N_k_out+ sdot_elyte_host * A_surf_ratio) * self.dyInv)/eps_elyte - - # Re-set "next" node state variables to be the new "current" node + + # Re-set "next" node state variables to be the new "current" node # state variables: phi_ed = phi_ed_next phi_elyte = phi_elyte_next @@ -316,14 +323,14 @@ def residual(self, t, SV, SVdot, sep, counter, params): eps_product = eps_product_next eps_elyte = eps_elyte_next - # Similarly, all fluxes out become fluxes in, when we move to the + # Similarly, all fluxes out become fluxes in, when we move to the # next node: N_k_in = N_k_out i_io_in = i_io_out i_el_in = i_el_out j = j_next - - """ The final node is at the air boundary: """ + + """ The final node is at the air boundary: """ # No ionic current or electrolyte flux to the air: i_io_out = 0 N_k_out = np.zeros_like(N_k_in) @@ -348,65 +355,65 @@ def residual(self, t, SV, SVdot, sep, counter, params): if params['boundary'] == 'current': # Total current (i_io + i_el) into the node equals i_ext: i_el_out = params['i_ext'] - resid[SVptr['phi_ed'][j]] = i_io_in + i_el_in - i_el_out - elif params['boundary'] == 'potential': + resid[SVptr['phi_ed'][j]] = i_io_in + i_el_in - i_el_out + elif params['boundary'] == 'potential': # Electrode potential = cell potential: - resid[SVptr['phi_ed'][j]] = (SV_loc[SVptr['phi_ed']] - - params['potential']) - + resid[SVptr['phi_ed'][j]] = (SV_loc[SVptr['phi_ed']] + - params['potential']) + # Calculate available surface area (m2 interface per m3 electrode): - A_avail = self.A_init - eps_product/self.th_product + A_avail = self.A_init[j] - eps_product/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 drops - # below the user-specified minimum, any reactions that consume the + # Multiplier to scale phase destruction rates. As eps_product drops + # below the user-specified minimum, any reactions that consume the # phase have their rates quickly go to zero: mult = tanh(eps_product / self.product_phase_min) # Chemical production rate of the product phase: (mol/m2 interface/s) sdot_product = (self.surf_obj.get_creation_rates(self.product_obj) - mult * self.surf_obj.get_destruction_rates(self.product_obj)) - + # Rate of change of the product phase volume fraction: - resid[SVptr['eps_product'][j]] = (SVdot_loc[SVptr['eps_product'][j]] + resid[SVptr['eps_product'][j]] = (SVdot_loc[SVptr['eps_product'][j]] - A_avail * np.dot(sdot_product, self.product_obj.partial_molar_volumes)) - + # Production rate of the electron (moles / m2 interface / s) sdot_electron = (mult * self.surf_obj.get_creation_rates(self.host_obj) - self.surf_obj.get_destruction_rates(self.host_obj)) - # Positive Faradaic current corresponds to positive charge created in + # Positive Faradaic current corresponds to positive charge created in # the electrode: i_Far = -(ct.faraday * sdot_electron) # Double layer current has the same sign as i_Far i_dl = self.i_ext_flag*(i_io_in - i_io_out)/A_surf_ratio - i_Far - resid[SVptr['phi_dl'][j]] = (SVdot_loc[SVptr['phi_dl'][j]] + resid[SVptr['phi_dl'][j]] = (SVdot_loc[SVptr['phi_dl'][j]] - i_dl*self.C_dl_Inv) - - - # Molar production rate of electrolyte species at the electrolyte-electrode interface (kmol / m2 of interface / s) + + + # Molar production rate of electrolyte species at the electrolyte-electrode interface (kmol / m2 of interface / s) sdot_elyte_host = (mult*self.surf_obj.get_creation_rates(self.elyte_obj) - self.surf_obj.get_destruction_rates(self.elyte_obj)) - # Double layer current represents an additional chemical source/sink + # Double layer current represents an additional chemical source/sink # term, for electrolyte chemical species: - sdot_elyte_host[self.index_Li_elyte] -= i_dl / ct.faraday + sdot_elyte_host[self.index_Li_elyte] -= i_dl / ct.faraday - # Molar production rate of electrolyte species at the electrolyte-air - # interface (kmol / m2 of interface / s) + # Molar production rate of electrolyte species at the electrolyte-air + # interface (kmol / m2 of interface / s) sdot_elyte_air = \ - self.gas_elyte_obj.get_net_production_rates(self.elyte_obj) + self.gas_elyte_obj.get_net_production_rates(self.elyte_obj) # Rate of change of electrolyte chemical species molar concentration: - resid[SVptr['C_k_elyte'][j]] = (SVdot_loc[SVptr['C_k_elyte'][j]] - - (N_k_in - N_k_out + sdot_elyte_air - + sdot_elyte_host * A_surf_ratio) + resid[SVptr['C_k_elyte'][j]] = (SVdot_loc[SVptr['C_k_elyte'][j]] + - (N_k_in - N_k_out + sdot_elyte_air + + sdot_elyte_host * A_surf_ratio) * self.dyInv)/eps_elyte return resid - + def voltage_lim(self, SV, val): """ Check to see if the voltage limits have been exceeded. @@ -414,15 +421,15 @@ 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 - # looks for instances where this value changes sign (i.e. where it - # crosses zero) + + # Calculate the current voltage, relative to the limit. The simulation + # looks for instances where this value changes sign (i.e. where it + # crosses zero) voltage_eval = SV_loc[SVptr['phi_ed'][-1]] - val # if voltage_eval <= 0.: # print("Voltage") - + return voltage_eval def species_lim(self, SV, val): @@ -439,7 +446,7 @@ 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: for j in range(self.n_points): Ck_loc = SV_loc[SVptr['C_k_elyte'][j,:]] - + local_eval = min(Ck_loc) - val species_eval = min(species_eval, local_eval) @@ -448,12 +455,12 @@ def species_lim(self, SV, val): break - # The simulation looks for instances where this value changes sign - # (i.e. where it equals zero) + # The simulation looks for instances where this value changes sign + # (i.e. where it equals zero) return species_eval def adjust_separator(self, sep): - """ + """ Sometimes, an electrode object requires adjustments to the separator object. This is not the case, for the SPM. """ @@ -463,9 +470,9 @@ def adjust_separator(self, sep): def output(self, axs, solution, SV_offset, ax_offset): """Plot the intercalation fraction vs. time""" for j in np.arange(self.n_points): - eps_product_ptr = (SV_offset + self.SV_offset + eps_product_ptr = (SV_offset + self.SV_offset + self.SVptr['eps_product'][j]) - axs[ax_offset].plot(solution[0,:]/3600, + axs[ax_offset].plot(solution[0,:]/3600, solution[eps_product_ptr, :]) axs[ax_offset].set_ylabel(self.name+' product \n volume fraction') @@ -473,9 +480,9 @@ def output(self, axs, solution, SV_offset, ax_offset): for name in self.plot_species: species_ptr = self.elyte_obj.species_index(name) for j in np.arange(self.n_points): - C_k_elyte_ptr = (SV_offset + self.SV_offset + C_k_elyte_ptr = (SV_offset + self.SV_offset + self.SVptr['C_k_elyte'][j, species_ptr]) - axs[ax_offset+1].plot(solution[0,:]/3600, + axs[ax_offset+1].plot(solution[0,:]/3600, 1000*solution[C_k_elyte_ptr,:]) axs[ax_offset+1].legend(self.plot_species) @@ -491,4 +498,3 @@ def output(self, axs, solution, SV_offset, ax_offset): #Jimmy Eat World - Chase the Light + Invented #Lay - Lit #George Ezra - Staying at Tamara's - \ No newline at end of file From a66fa95516eeef8526d3f30d7f87aaef6f93a13f Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Mon, 13 Mar 2023 11:33:40 -0600 Subject: [PATCH 07/14] Improvements to graded cathode design --- electrode_models/air_electrode.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index 03a4d1a..011be05 100644 --- a/electrode_models/air_electrode.py +++ b/electrode_models/air_electrode.py @@ -60,7 +60,8 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, elif len(eps_host) == self.n_points: self.eps_host = np.asarray(eps_host) else: - raise ValueError("Porosity must be either a single scalar, or must match the discretization of the air electrode.") + raise ValueError("Porosity must be either a single scalar, or" + " must match the discretization of the air electrode.") self.eps_product_init = inputs['eps_product'] self.eps_elyte_init = 1. - self.eps_host - self.eps_product_init @@ -107,7 +108,7 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, self.product_obj[stored_species['name']].partial_molar_volumes[0] self.capacity = (stored_species['charge']*ct.faraday - * self.eps_elyte_init[0] * inputs['thickness'] + * (np.average(self.eps_elyte_init))* inputs['thickness'] / (3600 * v_molar_prod)) # Minimum volume fraction for the product phase, below which product @@ -498,3 +499,4 @@ def output(self, axs, solution, SV_offset, ax_offset): #Jimmy Eat World - Chase the Light + Invented #Lay - Lit #George Ezra - Staying at Tamara's +# \ No newline at end of file From 51b9c1525cf3042e1fbfc9b939903d5ae680cfb0 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Fri, 17 Mar 2023 13:53:47 -0600 Subject: [PATCH 08/14] Packing up code to work remotely --- bat_can_fit.py | 122 +++---- data/Liu_O2/0.1mAhcm2.xlsx | Bin 0 -> 11404 bytes data/Liu_O2/0.2mAhcm2.xlsx | Bin 0 -> 10967 bytes data/Liu_O2/0.5mAhcm2.xlsx | Bin 0 -> 10667 bytes data/Liu_O2/1mAhcm2.xlsx | Bin 0 -> 10865 bytes electrode_models/air_electrode.py | 2 +- electrode_models/dense_electrode.py | 118 +++--- inputs/LiO2_Gittleson.yaml | 499 ++++++++++++++++++++++++++ inputs/LiO2_Liu_starter.yaml | 516 +++++++++++++++++++++++++++ separator_models/porous_separator.py | 3 +- submodels/bandwidth.py | 4 +- 11 files changed, 1141 insertions(+), 123 deletions(-) create mode 100644 data/Liu_O2/0.1mAhcm2.xlsx create mode 100644 data/Liu_O2/0.2mAhcm2.xlsx create mode 100644 data/Liu_O2/0.5mAhcm2.xlsx create mode 100644 data/Liu_O2/1mAhcm2.xlsx create mode 100644 inputs/LiO2_Gittleson.yaml create mode 100644 inputs/LiO2_Liu_starter.yaml 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 0000000000000000000000000000000000000000..a09fd47fa377a103e63e3fa5cf56d0b45b81a3dc GIT binary patch literal 11404 zcmeHtg&1e*`^DWMXmEEYxO;F74#C|uxVr@>xCDQhnVoNDGyDAo zyLF#iRsD4TPIp%w>vIlB8Uhj>01bcz006`Q!;=g%T`&M31`+^x3xEaJ61K5+G_iKn zQ+Bg8anNCK{b)s$3kgn@4FG?={(tSicm>K5M&x>cD5AF#kHX*SrB-SLVA+2751`Q~ z@OAda_m>!HWtyA6WQIPV2<2c{uvepvuXr$@kDAq3f2{Za)Ze0r816UNsiKR=&e%6_ zOx}To6KAifahi>eC%}ZGrx*Dl9gt$z(4)vIvMDH4hHHj_$>+PBu{5BHy58S}SfW|X zIk?=lA*Z^A$}#+98fAO(?Tm#vwTJHIHoYWdIElv!rn7DVh*({T+EB5{{ytf{-GQo+ zWmew^d<_^J(}uh{whrCW&Q2i9id*v~?5aewKgN=_!&BSSy@M(${0em6{2s$hk6tTN zz>nB}Y5@SU46ns3XSN>q>gC7ClNmYdCR%+U-+YVTyN+kvoEaTsEQd|Sq;cSzBnFX%U@2hEqALMTI>#AIF}KQ;yrxIK<1cC-Rc z&XW#&n4h5mfR`5t0O;SuvQ7m^diBa{GOtpH@G6#i_9j*ij10f^{}admVj29~tCz*e z%J%>fgHI(MLk4eWS7T6wq+J9hT8NdseI-}k)J5iz;jMO1;-e_z20)1Ww0S=cFR$`M z9u5-UY_XPwqoH$?*1MDkC*RvS!O&3HCyCjWZT6x%&tA{oq>4+qQ8~9o)08$A=SUB& zlZs8Bi&UeGGpOOhy)D2G#^g=))9jbkSTp=t0X8e7^rJkus)0TGFn&DUXFjR$0Ezz- zr_AwGD#oCLq3L3U=b#nw^&PH?k{PFYm0_km7m1smu~p}#NP5THCr^60&x492tY|-( z$HfQ8v#)(LYgx}nGTeIE;d)93&xZWN35(ZWE!01SgaB0?KIc_HuwO;!&FjjzS~0rV z*jpOf*jWB%xAK+OZPI}#-dVLTT^HP7@qw(x`NPqor86IA=js?^X$CmDXcxY{mn|uI zc^E^PDs(&0gnYmIZSVA?y=D6A$JjiNibxQ1SN>c5UUH=`F-MRu?P(sA*f`QAQV!&9 z{)CNPrg`I8a_my;^oh;F03um`7R!-^Tcim~MrT~;=nz52l? zAS_8ow$K(d@tjQ;s5DXEvrt{d)*KII)WO29eMrPyv(~>2sv`Tof`+#+42>v>1nVXT zy02LO{7#cbHkq@rjohnt)>iy11`u4u3nwEkj^S@^GBfN9ShA5IypspdiU~lmiZ(iX zpPyNI*JjJP@^xyGTSL^NY#%2Al2DA4s9ckl3Gy_qD5RL*X;)69M ztkGw!3i3M7l6gwASKwJ=`6G zGoOdA3F{I)98_416QAXK(56NW{n$>4SCpY`?B_vGub=h*#8MUE4aS4bjq z)if2%YB2_MVzzLYcsKm20Y_@T``!mC9hdvpHZ%~WQGeKZX%b$(fzI-GBb>a%fs^Dx zasy26&k)o8>=M0!y=Z;>AnP_}*_C&XxPxCRO;Ezd>ly4w6f!7I=tCXnu#=nGqH?1pKh%M%AeW&s2!6HNSpQ`FfQb4vEN4P$ZQ1{Njn_ zhCi*Wjzrp4S2KRH5l;bJY(YCsLyVt4lM6!1bAw`U7HRo}L-O77CP!@A#d$9Hd(qFfzKbY8~6d_O6sSqZzID&H4|786%y zK?t9U0O+GPJlw!aFV|&$vKwAPT2_I=lqeM}3sE2ver58ZY+EVFsdAeD%lOBc>{*M2 zN!mFnIepOAHO!oFLY}eytqrnum^6{nh>1~xR$l9 z!)^Nmp8Q0gyFoORB+jGkRNu(9nLAkrpSd$`+dWidjViLLcW*kGzoXZR4i(wiBM``K80 z_;%3dXho*r*syPnSJ~J1QfW2>od}Lp@GkLuihs!c?>m{0x@v(hjuNETx6hYXk@d_{ zL)7nvl7nr(13AJmw7h3|{q6YV+xvVb>}8tzaaA?1qLQbpn}PXC{`S^KS2yqb=e>2p zb|3eckT*6yZJr+6@`Tr$J^WjYm4rSowxjgG>(^p3g7*(x;hIFDl7sWPK{-m(BMyx`ogTBGwfgLdD#mH%m6$T(z8SbJcJ?F@_&Q%g zvzdEAabPDV&|*>HTn!Q)w7?CMu64h&KMgu-8XzXNuMm@637c9v$$*5xzC> zL`5g6rA0(k96Q+qY+y+;jR90mYy2H!pY`T1!*^l=l6xyrlw?D=_w9Z>66KQY%lPrH zCleWju^M0~v?W+kr+S1+b}K*7xT6wOK3tt9w;>ylgJ*GLkGM@XH}zAhQ{ zpgc=`xXTZ}M8XC|2hhndfw)eLW79@wl|mE4_~Tpx4C}Sa^0f~Z0K_n8g~e3woTp-l z$YTJ}_QES{Z@AIT#f5k&z^98zO~S%O8puYrYH`*TNe|9QvF3|a;RzyI+!;qeFucS{ z#$?UXoep%6c@`;#d5-%-_VIy7!CIAuA$d>;L_7xO%Doc3 z=IWv$X41^vUTo^(fdWo_%a`JIE#~+NXNwQKKCDSvsp-ubc#11f53G)&a@2Ge_E1p7ny|TfnTY)V09-D72eK zy3|U-bNE8M8AOU&YN5k^KKAgY`q&7A@9GF1=x7t6MPPB#r)v}j)*>sNF1+*wBWg!K zvbhO@XjTS+z5a5WvHKV=wm2_8^v=e3kL|!YBSt>C(pX`y^q>xt6l3dm9{PT=HLq#tdm7{>}v(RBY_gm^ybJ8K4z)Qbb8#hFY$$UUt1WguE1qOo& zO7h3guYCwXAaycz1^>+RDg$PfLa&J2@1ISS6O>|RkomPFtg~l=I2WfUzs^(XOzjXz znH)7|L8beXO_p;aRpw3sjp>={&F?!8*eXI{yEa+jrg`Y3emkuLIm05If*OSy?Oi8E z>_;BlgeSC>3?a3IP9ACv?%7p-20IaeG;8_nBv)Z~idGDoVhU3GsQ9q>Xr%9@q21_p{# zRPyNRlAOk+(~*h}(_80sz5G$obZ&~e5d=n?z8T>L2n#L=*+nhwkR30 z8)Pl4AA-1rDCX~C+d&6VjKfvW>vvS(xKm@!r15=(QqMenGO34R6Af00dIoRWuw~?E&E= z<>&%w`3_mSm6WD4Qn>hYTT#(%rnmd0&%s$xQ*v8He9DeH*5d#w4 zq*<5)`0oPNQZQ)?6w9$m<9_P;8-hH%n2AY`krJhRJifHbZ!*OSAw^KFtCduf=TvX^ z)CheF+E+YzN4#`lkjcbNu;npzONEB!ZFWcgw1?Heh&icBA+iIlA#9!f!rUaI@2NiN zaBp9b$Ly_sxd^WnCR?ocj4_#2hl_f|-(!Md6 zA*jy6rD+MnxO+e9?HXc$U|SsZxIzq4>9L}cj5-PXj%LIOGatYgT#-z2Yp4%Vegd|p zF`SBXqj#Bn8Wzna7>z=|QJES^O66xIpWkyKX;?jwIX_WHKU5u0@y6A#cSxrGkv4ml zorp1$*W>x@;4-9x*XL^gu1Bw?sg6e48irg+>0~xq@8#(UYfZ1?`RwqV{a}Xjyt}vK z@opr!*|x@v0$gR$M~;^H74Yu)R5Pd@tH(lY+Llr!Sc>a&VXqf?4-13%=4b1?pl zUZ7BKoB-;+bMeW@Ah@m;(P6h!)sOCFeD!{)R-O*a#b(0uCbFp1O0jCF6{84eY!|ML zlC+VtGuup=)%eS%(3a4NSj!UzsWQOr9!s{`mf~V$mB}tr7+Po{N7s~6N40RoyV4fn zPdoFx%v&;5p;sp;l>wio$%XCp$rTxDg$+NAbej+bE?tbSeXQ&$ww&elG$u88eGZ`a zr!PDGBW4T=$BLXHOl%6 z^uCXiw%?y>3dd%N<#k?&1o1d5$;1uzdQ%`1i^jSHNMdbsvl zkTgPlw@6E%{B&NTexYtOl!cq?@=M}U3$~JDL_fAi1D}KX zgaOOQUelx+iGJd{W2rLQ!7-KV&-!2{aV3^>)b`1+v?jOb2_LIWZS4GAF{wgP*Yodi zMOLMBICefOyJtip%WklN?+|_upQfIN+6Tt9bVwdW!fKc%yJMbs>e$kH+k73+*IVQL z(GH2@X@Ns@({MA2qG0LIV-O4rkBGJE{|Vz+TO^}Y)9NH4eIfn*8whkU)wZ|s-y{3R zLu)7>U4{ughL3aeA?zDc@SX!$O>JG6L1tBIM!)**X~k11w%vXr^qk z02?FunZ}kvw~|0gzimgsPszL#J%OkxY_i(T3ESCjm|idnbtpD z43coLYC=R3xU}QthM}2aQ#=xS>GkM+5|LOZ@pC`Ys`;YuX>O+p&=RKGp*V66x-xf> zng_39D`d(M`fVyz*glu^t-mWHTg2~fBaDo0gW%d@N6lAPXQYLR3iS;Sus{=-!(1lqdA;$j6VqB{*F|7Av{wT!=c_@UZz zDiS%qVH(7lo~y6-6+D|?Y`@7vJB|Cf?0|2CL{wv6#&6Wx&KAM!9yMLnIy7$w-J!87 zifhqcw%6N{nfqQ2kn(Ml*vwju{G?egtd`~vQ zBBQncQ140_XNi99A~IDZFE+x{Q-HaXjwsAb?DQH~Vbpkb!Hik&TD`9xDt;7yAfT`v zn2H>t8jf$B0Dean5J=G%JCnaZ{Sp?jR*1)ZP`TF!M~fpAi11m@ zG8izk0FyMB={aH~OQYv+kREC*yGRvtER?qXfykB*4$d$Cz#EpNUY(iL9&I^i2avjtnn9 zY&<+q9KHz=UU#CKFA&U}j}TiVS?dq4UOT0ec_47z@H<6(YDbDnIf+SYE`dSj0G_qx znBLVR_!Nuyw8H2N%uK78CSD3CsB*qRH%Y7jMj|J45Kw#VXzcKS5>&<@L^a+ymUTk> z?YLkz7I}b1ITd}p^fOFY&9GZGRcQu2pC!~y;-8e(GE|yd{GV2M7%@Q~Z7gdkjtfdW z;uX*snT}WB6%zw}?An7n_Bi`NahVaARN6ZlDKTw2(?YrA%BMLX{<(9OFx}1+5NFxk zLSkri%|uCnwCd_fCmT$K>1}DGta#zHtKa(r?^DfXm`fM~vTELSJBHx=VcnprwYwUj zV8IKuBKaC9S+>zyV~7nAVYneoovT=kYvYH(x|k!ks_dmC_$4RPw(sY-kRi2U`*i;Z^i4yE@ozz=?)m zo?{a5?}*V~106){;yg(;pZ97q__FPpa6Tq$)=LJU`vm7gwRVu7uBAnT2+H_Xl|-Gu zlm@HC-q-T4NFdA~8Q(Ta?%APEBsb6Cwn-1i99u~3aX5P=XR9Rj*WTo+0rkT12T%(c zGtB$0QuxPEeadhp`m=5qT%s33I3x$Qd2bACe1L=~SOqtJrcyrvEn5@>G2odqMUMw@ zu~h+osNEiWm41&(##&ZDdFb}{*u4;>SyOX3_72&O|Zz|*ip$Qu6ZuZcd z;~l2EAO<~?My=0l;FB6Dn_GF^o-&3il8DL9*yePQwstfyz$0gt9gP82Jq-blNhm^f z5dr9m9w|A8u0r^}Z3)4Ai8_RhC<`of*$&JBs9e>Oo@_Jd-8(9ot|HOgksy>)e;~Qu zBvdA1w`T#xYnRgxb$xGEIwc;3c*dnw$GlJi8MjRQWnNq>k?+D7#28agh7A*6$f>qwq!#v zp-9mg{~rHX4{FSm$||ZwY)tfW#(r?9%IZ%^4(i#Z$@fsQ7IaKEna0Pf!fyHD1Pr<{ z`Zyh~bV=gOruOO+JJQwa1fk#moDF^2F?=J&&}ZmI_Ld|B%CQB#bTm);D#5@luKJYl z8!4(7gKnJm&MQW=g6q#o9K%R&_<*4sx>I-J3>-oDj=~I_nj+ccREfUMGyP>lI`juM=to6>UkCFQAkUVUteedE7`?k%wV{ z)#Sq9;wX8t;4plfpN&|Y!W@Yo=u?TAJfiZjK0@0lX%X4+#F#5X?;Etpb<_2_RZnMb z(OYl!=6l51^6$Cfk)rxp5g6m&Lm|(1H|dn?Xn(@5X9zSO@W&=l^lj-)cLsiXBrQ^6 z^fMO^peiwbB7?Z()UsJZ4Aqn>w8r_wNaYcVMdHv48$VisIAzq@CAx5Nny4=cy&Kya zV8bL5nq4`gSBXQondmVxsiu&C*kyzFEhrYVdr@1`f`PfsM+m#!WE#=XpqadqFQ3T% zTeBSOknEyunk4*obk*RX8EWteZ4|))>au3%p-;8R`arhb!1te=PlaMFguZ7CTPBT? zXG#wv)nj}wQXL)8#wA(zVeNY|ksRU}l3W$2^fTlae?&Q%8y^l)S)z7CrjvMRz zkj^${av{9)p>3I`vZo%I-05rUFzx2)Utb!9&rUhO<;RbEZieh_Mj(-gv7AgT_d2RQ zj51n^SBlo5!5iGxy;^M6+>#a%@TU)(o@^KJZr)X5xh8q7yA-|A-LPx-u;HPx!A(+m z*GxI)(STiG=U+piCL8>-C5Eh#XA&-ND6s2UIB@8Xkn^wA1!_W%G)1r8$N1|X{a;>) znSs5Dv67>`xwYx<&WK2qApC32K=jt%1@>#h7X}eBNi2U#(5*0cz265#Zl&sm29m9d zw()?|61Bt8c(Kq_UaN34M7xai!Sax!omNg?R8i}yaT$USXvN*`J0df5^x+|LP8^!F zPU3I{gAnI6E0-xCEPO569*G4B@%%3hEj264fw5gPZid%Y5G;0{GLH9lwrByn4ehP69*;Ls&r8L{mU#DIss z1q=xnNdb0@^L82m5mYcVYs!XQEXL0fyiW6%J%P&4oOSTrcrIN1yr$B2C@{BkrCke7 zI&+C$5h)#=mD5Fg-m}W7UTWYd%eT|EgoV|^1>kp8Ml9{-apy~YS|+YRKHN0a_Yuk* z8(_w0dIcd7lsg6Z{9}}#Vm~LC%PHD(Q*zQ*9$>hYo<+O82V67tzl7c%fH&^IC=SPAWGk<#KSw1^)GLo zgFrjj#QJ@H7>Sx2JV!FP2m?_(!ZRmD1Ck<=e$#8afcF( zZG=po>d~Xq*Nd=7NTK!o5D28J#nF=Q2SPRbCJc}e@qu5 zS3}8tTp0$<&uz#8E9u>nOx?U0><0E}%j8-Da#HY2I5L;1qmUwdwPvRr2{Vh~HIMk(k2uln25pAMPpyCZo?Ei(N)n(+9u`Tq4Q z^AENK1E+sg_W!x=^UrPibNmlgpdjggC-CnTYJUQM8~>ex+f!Al3NWy^0C)f*005u@m>y+W8A1U7@vr~@4ge8aU&6uO z6=d&f^wP@_iu+6YlCu8YnZ<&$hOH%8t6n6wkx6;jO`*Sn=UH9kZ&nw`&NF9B5TX1qTgvX&RF7 zviA=j&~*|JB|7Wr9_Qkch;kAc8O2&;0@9rtd)1#wei4(eAhyE87Y^LYS{l^CS|4ae zEz>LEA6o9-P}W+*;u|TQ!Q7g{nYFQI^f5f&VwHmhQ~RvoyBij%Qt43V zvn5sdd6|O!-wzEp@pY*jnvwe27f&I4g_XK_!o9-UZ8`4}#~kjWqV= z)(<&4P2COP`T-9BJUzhxRR3m{^_m( zuSir<>E%ETKbC!n7`mQYjmH#M@DP)2rF!WfD7S)9A6r06vf9Z&iusZ_1V$#H-Tz@^ zc~vlWZ;0w@^H~KL8&`m)!J{%f?as*!fr;KZRobcIOCOf|+~wR=hK#%y(7ip5sl2Hq zPhohSMtbH{vIcX4O`8}Ar-(EhUnnC;Z$L?R&Ge=UYEE3^urj>5kvDfQX(BUVKDBrk zO(c?E@nAXwZ^*^eaMb8S84dYqk3Z76`B12GUnhzVgpX2$b1 zyN83bt*L{9?N8n+e7UZ&%z^2jQ}@(;=0hw{7bIO6*^op6K;)*?Eq z_`|UPgHql_m{{_H)A;Sf>Z*r_Z!H^&D|<|eUQzlLqO`bEmR!B!>2dk_1|rRV`iM?m zitK4f^^N((?aZlsW^40j4@5w+lGBkvN@(78GceYQ>D(oskkP*zGaBf%bLS%t7RBF4QTsKuPGbj2z!JFd9M zVB%m(ceBiJrB*nlV!4z*@a1shsRBWL2GXm=^iFAwr@@0vyZ3rt1;q<^mc;OpOT%>H zggWW*nc%vZ=BL1r*+1Nsd?P-Ps}lRN^ZrI(fhMfsOgJ(3|?Bj>pjz_F9!=0o zd#I1y9m6Z=z%sM8J+NygM8982+=1erlb=le1(GsG9K4vw3qF+a87oXtQ$!+)*ol0X zGuZKTzAta**wzbIB;Ge++^%x_neC5@%c+kZWDL+O8U*DSp=BKo3mjxF*W)tZv!HiM z?vs>qwaYt+(m3v!cr@F@@N&8Em)OeQRTn*8NezL}V|I=CE{cq)Cs_IgVv1P{T$-Ro z4mgmFRBe`ksA;?Hr=REf9Ss1|UK~bopOC|F0VThyRCyY{wuL@!x$^s4FS-abUKi{D|Q2 z$aKfUT6ATnKG599Kpd`PT&Abu_PhRI*!dU?qUm8Jd9Vx|4)(#`e)%B^`db&AX5qvkmTX+UZ#O_XqKJZGYRQY{ z>F{SZ;~Z{WHZGH7BM|>LTK6^ZE>z>7GPtg>k@790#r9LP#PStI(6!?ZHCHYuvWG;0|@0v|NF_j7P{f+?6t7NQNSSG0|l^KciLi_=~If z{Rm^{@)^tM#j(ZL;la6e$V&gy{kdx{D9S-x9}o}X0wH(&+x@v%fk3Vd`<3CJKA0-E&R9Yh`g=TDC}F4TE7vz`&1) ztynE{xN!wrF5*TS`7Xpv1@4KR?}K%2GEo>SUk!9&h*pE&Z%rMJIohrEqyd5`yKiDr zsa`Ug87Aj8J1irkj^lo|CO?@U%lJTg(zaNLr;uMV@hmO;2(h#^)C4S!oUT(dYutmr zDS8@sd#|qLuQIcV{qCE6O@&_19x@hpsrkZlRrxwuqi+K*gf(=o>5rVno};hBc8(hs z6NzX{tMl#KVJ*p8w7#1drKpD5+7x~9;fY#t1%$3u+tEeC)XY}Y+YE-IDnU$FdrZQd^g}zZ@smnE#nNl zQ4iy7^@rzYhZ(u59$x}lcVq<;?5N_}GYA#__{CwrRx&w#)t*EKHO#1g(r^AD$0V1~ zDqcA93-OMoUZ3=Ihu?&i7PD6o;|osc1^gP9RKbU#)6>Ja`T6S)6UP%OgU10sd~bf- zc11)qc4l|Be`!^5_PohFyxpETe)9G6^YMDT`f)1?zWK6yGi@3m#o>P$o)aKAb$U8D zb9xGTJQ%&SZ1cN{G;8>H^6(Upx*}Bcq?Ia5uPx7<`tCYzCS&Wu^&M;JX0nNVaen6` zb?f|r6Q-*AuoAOG$am4;G1=-L11MYVvrD2&OjcwSb6bNhDB0RY;iNu3K@QV0g_qBBUC5)wpSN26nOf$qrN@bnG^(%3+%$ zF&ZS0w74!^x`vXV1~;7@d#;RFAQiCvZCKMfuW1&4N7bMY6;lN?uTG^{3%6;UEi~@2 zQB(3~G~Hu`AtIz^`N+5G{5`1BK<~=~`#y36q&<8CW>B|>cKw~PR59c1BvZaTR=HKU zvB|7?4uW7qb3-Ynvq_b7ZqE?1AcCNwrg>ta9u#RG>$&=FHnO0*CR83Wnryf_staeB z_|fc^lwNogP11g4Xl$%HSuq@ilJLtzGydX}p)8-aT+F#Cj5(t5<%z)9H^Hdp&}Jsr zXu}#-lCBEu>eo%^4n~3WGC` zso4H@6C2031|)n^ZA>5nMaM-Rb*ZsXzR{$byzpMe(kDHT)k&?x6&5qE&R6=PHa^>n zdKR2;XlJ#CZwj`hRz_5kAeYup`otPdaAcyTPl={@pbeAj0hjG9{j8YL4voE06lfd0%A}W`WPn+M0NDC0&NQJMc;nhDNW&oRo*~a%t4hQ-mY^UcJy58+h zS8)Z#g$OrI9}qT6i-}{oUBD1)$VliH1EY4b;aHkJ+2bn`*EX4;Pqs*aoLGEU35c8l4z9|hfh4{)YTHE zO%3~WpW74Y+#uHtrMhg+TA4eAdMU z1#}DJOv^LFUhPt*eZ4Gy@>?LM6q>h0@$O=(HxL*@SiAxiC z|A6w7Cw{q@T1Ht#l}jZS`H4@=T0?S5%r%-^lib?sz!AlXz#4@!l&<_@d1Sd;I=b6t z4nt`hQ+XS#NXv@C_|zOD=&haQ#htY81Rw&-!s670pIeQ8(oZ)5W`{BY05CxLGjjBc zeq60Ub|Cg&?Jwpz&>aJl@Zxq7+=-yNcs%fI#nG*8OuSB7r8Un=CarHcP}Ab(N^B#+ zfu`fXR0kForHMNTB~OdOptWDbA<@?DQcN+7EzneM(`LSw*KFB<{uYJhL18(^-xD=JHy#k2Oc}*Cvf#eQYKxiwm32eHl<4PJ!C8A%*fmj8WL2 zHzLd8z2*LVF0PS&@4UjHp=#rWt*WG)b5l55SiOx$^AetU&(0g1HPjHX_C&@BwRp7h z19c5W9ctcfy_h2|VL%|XI<4I1a6hWbBw}q77|6fTw@f#Ki0u%LO>f-vk`YZ_^X7wa z(2ru8=mrk1!jxCWQ92~*SDvPQ!-@?L%(-*ClzNen z5$H5Dj^^Two*o|w){HuToa~+Q?#?ot_Vjf=+>WMoKApa}tlnI~WAE@gJKIenSogd9 zt`bMKv`k8tenJ&hb5ixi?3f;G5`=z#3MG=&$06RAD2lb?UUD=#1Z}8Kx!2=Xeb}>1 zsxtuBCfI4a*g|pIO#3F|gLDntidl?1p~s7jvbUorCyv>Qt4Zh0QLRyv3ARUU@)dyV z@7%fCo9c_P)u3-^(b!SNeBIL;oi!3MWaX_Ak=ye^T$_s3Q5Q#;A3`E$=p>wt>D1Zk zBupbmdq9+-OJ`$ib{~36Z0Cf0&1p;^zk|a+vUGc3)QU~*K%HNblZUJWZ2pzYIiOpH zzKWez(mR^*{LYLeRN>o;h)ZIzPGr7U*wzsLlFzDovu zM7pM6vGumfUg>g7tZ@14rRicgck^^#W-((Z`~fs>n{lo;aCvnpW|Pt88#F?He72LV z(pIz0O_h8!3fJE%?%@a9HWh=wQL?Po!q993HBNo=mM5KMC7FhG%ZXlufomyl5>-zE zwA%G9TJiFwCpUerSz&}5GhWVQ<-PK`XI$)zvSo8c%iB^ExWrpeMuS+1H2LBkJ)z^Z z%GrJYotGZEY4E_s8A2093lcRqLNz0Faa=>oP{YLxohetHDb>BSJUwZrZDbBRJUq|s zs|~64?Fg|Yv2Y90=UrRt+zpJ zWTJUydyBI}*gNkD-%K;BPmr(C+n@VwQOaIpf4I#jz7x!*U2hq7fE~C zDPrT=VP1Ua#VXX%VP{5o6BP&!ao!u_49cv=OW3~k*fbuons1Aq=1+b@^8&2uo}X@> z-NO@A^=V3lnSxMd@ojhBKAT{A8O&9?$|^6#+=vUj^*77nex|)y*QDJDfv;%gXA-lu zOhxvO%B;>+wrtJl*mGX;ZC_^nz%bn$tT~^!lyuvad!s*oR9mj&W^K3_3+pJ?@SYWWC(I{a)9k?Wls^$~68K>h`~B?u|<0ws+kkGYSO%x1>h zjDF@FfNPu9hNu4+{POf^bk>LDaJ0k#0JMM6uZyeKYmm#&SH3=jW&5q?n11Y&0hq8C zQ&X!A(97V3s3o|DUFghwzZB4-NdcH@O21X=sa;}EoU3Ch%9MKHl>*iTswFzca6)r| z%JQ?RjzOuYv)2{+yv0@hy2__RVjR*@J$jePj%M&i+)(2ZDAB(A1B@{ zyO5fpQPw(Fh?uo?@WgO=znQ6S8=iN9@6`P!MQqx{q9xtTw&|#rC~a-JP`t=HjS2D& zbIg?HqO{g&d;Pf`s4`C_?UGE{R-n_8D&3f)t!o&14r{QgXHRIYUv9Lp)L<1e+3Mr2 zj*I7&!G4OiL44~RmG8O!YW*aZ!&El#;W_Onx1#>eeS;@WqAl*Jhvf8o73opI-Xi?% zOjHRjDmTbWI=k)<5B&J`?$x`R;gScLdvbcm!Rgpx;7C&2B*#~DahYRt6Jl+ZPd%_> zD3av*;f(U`b9R|4xw}<43Rti?E6*ZnQUs6>P_OT{;styjr#f=yxP{ltPd>o3zfJ1k z#k!zn?cfX!W?dqn(vR?&53BUEaJAXTgaOa_Rv>UH@9o-!cbdGJPIzA&(!c;8-sSPQvy|@sE zCU`uZw;3x7Vy1}dLK$)C<}{O9ete#7O08LLT3C zHL;FZGL?;-(MANtP&WH*qzO`e|)6s+h7bbfXE!VE!?q7v0 zK|Q#ldVf#I8E?fX-D4@f40LOCuI-=Upx!qKs7Nt z1|Mj!t(zX-{%S@%pYG*xo~p>_Q|@TPu5?xYiq8uvQE|1CWg!Y$t4Cct2vwHX<*`aK z#WS8k&v*Tg^_CIN5lm=ngw~zd!V5ZSAYsJIG&a~gF)WRw8jMbULY)DEV4dcJM zNWi-^zaOfP-=~=ut&KAPztO21;yu~iaSIMHwtDO8;mT=*4wL=n+v4;4Q^Do0qLLF0 zyTQB~nGchigE=_~atn(yW&pE|N}(zBGD2zSm|6 z=Q?u|*`?_<$c5krgy+MxbdQvL?ex=BwS51l>x{$V#Up*qA?Fp#$71mI2*iy*3kN8UBk zbP7Nm5d?TBS8ob!by4?dq2M?MX29N4fUuu^%N@oz*yg+qV>1FZ>Azcpp3+VKvRTmM zt7xhzhnnV2Xw3ra=){D;vG2~ktvmS4NLQ3^3XVcYQWU|;mw_BiM_j{@Q(J;*8U=BK?DY$IZ9^WE*NRu}JcX?-3J}vVG z_2h+qJ=a1RY~nIi>AZ}7itw4aDtmkbQ_$C;Ttf;w&AXu_v3AiHdj#=sy!~8-btup} zjP+7fnnxJ<;hns`yj*5mMrP_qkgVcj{tVK;C+7F!_-uF1Q2BmU3hSjf)dOC#)J~Ha zQkyu#P$=ZFPp9&-nFVBj&OpPyNJGJmk@=A-9`v?JwJyaZ`B`w%^c> zLI&&!>s7R!kO8Bbrh%MhscO1K^pqQ2fdCpo>|S(e4V?tE3}%56G=lJ#cV^sf(XLeY ztQl0CK5q&L?BE?V^eLSLx9@0heI_ptPN)y!mXvf6*?EA&Zp6#4E&LW33bn%tc9 z4i0_3r}_OHtA;2>8M_QP3etm}})F?gq zVtb-UcA!rb0kumXV$xU@>aiE3SHI$O| zQAlU@njak(NUq2}5O}8gt$V!3Fx+tJ@Np1iav)CCBx_+9_h}?7atf}c*H#L@X^p3C z5}ZOR7!K|JArU(Kudu@Jd1`I(2d4Ls-Nyvvi~H|gh?R*m$XvtK+1lRn*JecWjTkZ{ z{w#GJ?1A{Pv6M}cR*oQ;LG@ZfpdrXYT|lFzv5|W7tbHQnxJ-L*EJ->lL+CXa8`UW* zbEq<6f4hx85KGFwdP0%h0=uNgds}jrg%up3?8c|Z>?VU$Gz4>6yKt%Xb4a4o^tLS>?;S51reFHAS5bqj@Kg z5c)r(8#u9jjRKI3moy|Fg9XXQm^+xMIy*SJu$wtJgZ^C=^S?9ar@WP z4g{XHHBP!akyvwPfK|95&FC7UK}v-xbkklVh1TIpPJ~bo|Mg&VgiYIw(izlv`R7^S z1fqObq!!c1s;N=aw`+NQ2==rNdN4doAu?v}ews)MQv1e7dMIqb6FFKUm9WYwdRPU{ zD^r-a_be|IakVRF1p5l~%|`nLj9}OhiA>Of6v+960U83BQ%H>cfpz-x7XU_%XSv<5 zapx?9Mr#1&dl3Zh?VeK2h2^Fh#v?lY=mNeb4@0zps{Puk1GByde*Y{Xq3iPE5z-dy zz=gb=@OuuhW)<#Ggiz+WU6SH+wl$_+=7(i9xhc!S9w-dTq#Gt73qj^l=d-4CYOVNu zMY#n0MsVCHDTSYbo|B z@tAF0+kA84j8-KYUz_|u3$#qh*YBdM)E9Y`UDDHKTcC(=iBCfvMGN@U_Xwdv%2Y2R z@dX(w44*%jIMmxL$&v8QBTkSTrX`G!;Ga{r^{M(ELWwunvt|q{2R7@A1w~w_R#Ki} z^qn)Fc5$4&;k$i3C2IH9A+)z>M%QTNxFR5L&>g2L=#-{?D^G78vqmw9Vp<4w)()Ai zW$kMjd}sfe#JS&8=jF^z-sOX7YR>j95d$}!NN1HlB-r~ewS|Iag?Rh_TtfKAGX2s2 zhk`;?g})2S_X3bVp=}VS{#!Z7@8G{@-2a3E0Fg+)!2eJF z{dYOPr(yq;)Cvj5|9?r@-=+K>P5)EM3+!J~evhsHF5vfgz@Gw|iGB(AGdl1)^mh{f z6DmUeFX-=-{=0;~v*n+7NMf1>0Qd)Y{to}U>-#G_n(i;~e>*}|1vtp+0szR6PauSz J#ehH0{tvL(PRjrQ literal 0 HcmV?d00001 diff --git a/data/Liu_O2/0.5mAhcm2.xlsx b/data/Liu_O2/0.5mAhcm2.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..677e9ce0b7f942cc97a5420bd53c00c3415553dc GIT binary patch literal 10667 zcmeHt1$P`t(rt^G(PCznEtYICGc&f(VrH_%%*9pVjuwk3;-;+wy>>@6VS#< zU&Y-H=%~x&X8oQd4-%Xv2LKMb{(swlaSxOyjL7w{po-o|JPPkHO0CqqgMD)xFn~^{ zDA3s(-(PC1on>M1k`?-ZDwK<5$ytp)zT(MtHfmmDV_hE*(%+(l816sVsj7$1$=o+^ zMA?D!D$YSm^CSlo{~hZqef>z&bU=!ILyyvHkxfCVay)Y^YyrRR%%uS}wDtZb#8Ryi z?!o1*4LP+nG_K)K)2Q2%7&Dd@w4Qnw+l-Qs;bfjG*e-g73Zxp!v_?uz4)@8@?T$2! z?6U^O;A<>FF>T1JW9!f@?VLoiukmU=goaO+ zz4Ipxm|6fRScTVOSFqWPd-uM@%9j~A?Iu}$pxneD=v~J*Y0ip{F_FV1vu4;R9LN%p z=jNmf@F^Y~YUJvaKQ<=tvMHK`*?|{}+z`#WDDozg`|E zE8oL{7<3}>7(940yBdQkB<(6F(L$=?<0rX-QWu#|fxp^8O@OL`7YHHl+vf8)yu8X6 zc`!(Nz4f|0937LFyxz4UDEZ#r8HSF^AxX@>e6ttLW%g?JI#pcCoyMgtny##|Bv*Q9 zom_1COr#ohoJk!I4x^AD2%A6EU#nkMbIs_s5^PpT`M4ses(~}-AbvdEcRs0TAL(rf zx6IL0D%PN*k=bIU*Wi27t2;bZWpi$eDx)k19x``*llPq$BIz9%PhO028G}k>uhDPW z#>EFHbFO@~YG0p?WV-iq!u6C5o(=_s6PK)kGSokV1fWz7mkSCITu=z1fO^L5J+rH= zgO!o3t<}%GRiLsgx6Oj;t#9$tb?%wQ_1a&opaf?^Z3$rsxLJtqK-h2=x1hP4D;nav zZa;-7rLPf!&O636U}ie@b$_vlVS;~Nj3(7+3V~%S?vu=*YM8ISt{XM(D`}t@0;PKx zabuTR&v>@nywo~lLbEV{L^jOPs;u`0X~bPpR$f&n0QkvCl9AzRm^1O)Yj*DM^$d?J zutpsPLR;9RleS$7BZUTDr5d^p7TTzz2$ufsL!v;Q+PF5T%7XjKO1`3DcoOV9)azWZ zz7mt6_mqi6;ovK4F|gLJS6KNd{78@Fz^A(u3k|eh)z$L*U29LlJjZuZn=nJB_u#f7 zeLAu7Vb(!*Eq8j9Q#(VZ_kkww3Z8_OD4Uc5@A9P9s@=xk+T2&dd~8l+%2U0{Okz@o z5`%hk&4-a>Mk^3hsCChX(p=eJN7ISdFNkaG+ZcaQDE3N!s<}#?WCxjujywv+aK|Np zMeO6QmeL5dl-9)N+~;V^%C4*(c+4oD^NzX+0(gW=sP=9NO(C3baZ3cV%@Sv?Q`NPe zHovCbr7XcEUPNf|ywDvi6K_6#6xcUpeGC?&b+mcda2=`e_9cYKDT$Z+n5;(bmMz;M zMCCK-If!8Lb~))75nMjXhH@IN$(t{9ym;Ce-96r5<>p>|XMs0VVm*Azj!l3#Qfzv% z4}&Xz`va}kvRAqTOy-ur&Q0fKqzw-OT%4B4_-mrT*T|JH*&W0z-DMc;`-Y4cDw3ff7VcF2z7WKr>{z_1LjL?&3+?} zJw)U$P=!U2;B}^AHL*Su;ZkW|(0S^uYGhMbr$Tu}$JwR|svlw}1N;|OcQ~LUtE~=t zd_lhYPjB`YNJzc_{n8BbXgmNc7|5Id#rl8s=|3_67-&KUii!X3tz1b~x|apD4dFSM z#Wmdp3vJPfne0gQ5CwLqmUfwnl-=iY9iOgQS7T0^8NxQm{bY2&<>n3Y76iss7kx<> zI;022`x9O$+tG_r2w0=z3ZejUC=B$YgTu081mrYFj3$A&L53_$u5Y(Al&^yGi^xPS zo2G(ZTaK|fvspS$5DtSFaHMv}@BV}GQAI#)Ljy?~t?BMdlkoC2basF{@#MuTI7vPv zcfjP{3@QC>m*_R@dF!L8toxW%SN=WXEU)ngw-MrKm)ssyp zp@(7Gj^%TPk;@a)ouPr*brW=k zh#Ru$A;A>ANxVl1@eVjfv@GXRg;5=?JgB9<`g0*04f*{b;cp zXuz)+gwIHx!+d>AmEea~pJ1vt5q(oUhD0e~r|bYk!gNB&T9d2%+;4S=-3|>43D=KP zi<~?3>LL$%)O()8z^#xWsz*=AV`HAug2!h4_lL?#DU zB8D$y6K_vdOGIp`-P_kvGCWBn$8$6O*7*B(f_6<}dEeIdc5hD~PamG=$BUET>}(V- zH?F6v$L$a60z;dd#P#2kQX~5g`P(5eW{&s5ZtXx+KcbaGBZ>5AFn_TC~jY z+_hG<7{iQFt~rNblFISp@vzYRD0WC5E^P_=`XP7`>Do-3g<4$FdO7VA1LrK+#t=NcyfaPJ}XcdDy7t5P~kUh_Nxemu_1!5qb0wNQMom z8k<}<)-14}Olz}4Nd+9JiRI%gW+zYzi>#Fo1+c@BQ%M|W5F=DfOs3n2VTXi6K6q5V zk?$|i@APfdSl0~X|1 zdT##&H>caB0E!B02^?L!PK=Cz1;!Ov@L*CYN9nlWl&;`xb1`|lbrTO$_52+v8j_(Z zl73S`*E81r0xrF%3{A>0#X-3P$T{!$D{ZOhX$g_i9XX9$FMiZ)k*Wx-O2_JGUJ>m4>)qmUFsa<|@LN+vC{O$Ca)cLbgp*PZ zeU_@m7f!KVj+}1X9~&G#o|kKJg^==mL4gE=Br#DJ>VI zaPjB1;^LbuACC*Cv#GWe5AhU)?zklVc7=2PjEIF4is>X@D=wPXUh)zN5oka-n>z%C zN9pBk9|tt+A=rThh&z$@xk<#5+{<*>eXv2Xjh&`k{_9i`0}|cj+1LXFgn?@**z|=; z6}aSaw|W6a3ZC98r_e&FDY4peW92MWH zAY#Pk4N-*8g0uoUJ;9l##b$@|Ihgv|J@eAXdI}AfRth4L4vj%fAL}e#o0hOly7!_m z)(`^)+u~@)6=RUfj+B&TG{`u2wIY772>|@Sl_(^)hWZfYCtzzD!)drTdY37uVbN`a z(5Va>RcMi+bD%ycH~tMM02 zp)H{ku~t8rq{;y|-`R82x0DtmtAO8-!q7vDxVom4JF0~v2+LZ8Lw4u+*|ub=LN9-y zehv(orWAHCpj2Y26*dYP=?0Q~SUMkFv;N#uVl~U}WkPNUdJYz!kWah)Bj!wsM@rlx ztZxY0!%cSB9DKXPsVbQ%L_ET1FYb*QK1hG#3BJM;?0`qyMt&gD?d6D*cG#P03dd!S z<#$<#RN!-5l8GDa^`Sy06^(Tbl*I9hUgQ`I0aON1DE*f7wSes*?D;`J|!--;3_*s^y7Lq2smm?7_yIiZ<d}tgFmy?E~DfX+qJ~3-0hlR;6^gb~98wGNX`X zH{O8n67Pgh)6PTfvBb4>NFGMQYMLc`VE^#awWIg3{W4&nzs7&u4*AN<@)gN-!}Tbt zqE!H&VGt}lBF<_+2-dTXNM@PV`yUDE3+Wsv5SU_WZ5Z)8k$vN#HPnx;!$huNj?hri z+Gh_5Jn1yR?;`$R;BJT+tUS!KRyytC${pw2RUP$aZ(kltO{et)J_E){D-Md&iAt9TdQ~3NfUmIK=ri~m%r{b!a3{A5*jcX`lim&L5tEi7GgsErtL2oftW70`uQsE#u zS4+%(7oNo$3(IZydPB5rH&}2uqWYmrR zR8-0ny}UMeb(XLVbmXtC?lAUnd*0hUW!Rp_nXkluaolh#JKW{XX`7#%9Dn3|B$XWgY*H)p)VA-Y%?v%j**YPQPtgnjv=o z{J{hcgYcDj%}ChwR&Hv>7a!X9Fso7vj0t##9jGdSO1WtDWglH|)kdDGbPk;9P&N!= z6?KrO5AVTr3@1cgEItp+GyT}!WgPv879K96{}rch^M0#cXVJ_41dSx~p>8-w)7r)I zqiG#K4P={K{>0H7&quuq?Gq%A)v8GXzSA~xB;=)88-+to_&SFMZCFc6pijPL2~46C zZ%G6GQ_d3!y|KYAt+z=#;Ks70{^dU^I&Pi~4%(o?gqk1#fb_5Q>*(bE9_aY*6$sof&_(k6T$#MV9;YDBi6 zD~&deXoie36x)<9zx;Z#eLyty{C&ALXHjLJW&*Q9J>Fr*Td20&%0#ZO*oJXtXBz8$ zh0hk}+bH?yCvmsSjs(U?BsC7@Z;e~q-$bx^L`_$<4$a#`cW8bS#WQMTP!nrn+OkuO z6SFW{C|cy4LIrw!v`ZIbBeBqEeLvSmBR@|n=9oazny=BEB-W6vuBrFo0#avH%Ldm% zyG(y!sop$dqQ%oi2@}gbmH7;51MkiuG|z2acl|Vy#Yn>Mk%MA{T}FHFq27%=&I1T+xwD@*Tv`Y%ccGiFZ#w8q5?O@OOj}_jgPL_wL5aF|4az|nfTvuFtV z#2X02gDX?ycm1}niuQvcT~6NtX_!wCVXMzl3a`OFnFk`b4gV9wr*@>Mlpir^&80BNTr8)pxn_6u2)-pEzO692 z12faAW{DT?6xFypFljYS@yQ%S`f zFUx=ls~L9Bp()E`6tIH2PW;2tTBd4qOTg0#A2YUswXIbx)lp%oXS^ahGwaa`yi(!^ zU;Fl;j_=(43UOHx*fcu3nkg}Dy3;~=lqx5=3UBAm*u(TXQxv$%=N1w}qiZHg1Etkg ze{{ZqsWiJOieami+b{lY_;TLDvL~nP+7_V3y zMeO6e$hDro*JKLhIIzC5PS&cI48-&e%7bd{pgdVii&h{if2*b}>I|klSS`j;`*uYF zVgAtMrcv^{J=#Qa^9){_^l;3PrPOyW7w_a8)ujH~>pXQ9{cwT-v?AtAi@wX0w_|9& z<#-bP**6QW(F?&`k^|fP*M_#fEX1fdh1dROQn!GXEh?fI@T{5Q$Njk2sz56~KZ=r5 zfD3*k%sY83qR!#QGXU%`AHY?vYKwonldM|}0mIHO74m@?i2nLp&JfDcF6-S#CVgO| zcE%d`q-M(IR(`jajFGA&VzLXa1p}m=Jsk|op$q%2=D=%x&39arP{bM{@1QGtq~sjC zis1XUB?Jp3>JT=fEO9VpJFo|!^3+Oua?GK>-O8Xgvvtf_9~`R@$o=n-41-VeqLw!Pyy zz=?_S@OBc=AVz98&`OkV8m8rjws+rDS!P;>r|W|esbr&{18P6ual5n3Szl0B6t9Y+ zsR)riVkL;~H44JD3Q-S!z?&@k(4aJZVCs#%jfve`-CYy9;;2E=mTUwj6e&6r&=U~r zNsFEGxr$~H7aOyJxgQ+rb9G3Pqef0y@;#KSB?IeqmdVkouzNu`5tCkw!K)58h9q$| zGY5@{UFm8KqR^c`T05zgRl?Pl)~WRsQI$s zFanzy#_Y~vPNWZvsidr)QTaI5(6-9jB=&qU7AnwthAnd44888v(^*@LHk-W#o^f_> zzq{j;qxrulG9kE!LZ0t#(yh?d2_dLw`d~2-fJ>z0*V3Er!V>aGUaZXQZy_E?Q)=== z0dc{tZM%dRswGuq^D2ay#xoR$%&`|XezX#C%DA;lbm9CY(LfUVTWo8fEvraq&gU8Z z&##C#6Fo;J)fE#EyKM2levHNLUeuAaWMXUc6~b)?P9qu_HdB5UC?Ik8+AIeK9FNSuyf1(R3z3y?03qv1#FZ&Rel($9us(x z>ga$rDb2nQYyX}i6ybcGN~9j7h^W5~Xu4DCP9ELRw^q)nem?1O+E{-J>0)ax7tB8& z+Lm=9d*YeJo4&RV({7Re<)u;h^n?pse*CECddR_c1QK}|$Jxy4dq=gWab`=&O7S{0 zc!T@8cZ=Y)|Nx1x>4_(i~ABO%V$NgS>RTp}sD+bLz#z9}qf6qe94IO|c%1#az zHfFyjBO+0P@ZCR);Q_9&Um89!iBL%51W+s72=mtan=0`tS2r|}ZJoD`2cDFwAB@I} zg{Jbq4@XC|&rBbz2tM3x<@Q4pwW%7HAu>fT>Gs$anPFfI50-Q0(xP`3hbtU}IICH? zNKwEc(5CN^Sdb7e_|(u+v$Fglwrj@S=&A~W-QG*a>E7NBJ#gBt8|<*;%NM+ir^e46 zyL#1y9rYu*r~k34%%|zc#{-&piGid@G*F4b#MW5B!Pd@^+1SR z4L!_QZc%N+rY*n?_eA7B-_2*(_{_k&hSQH#4-J%oPgK|8%>*@_!{9UKD&Y z38O)+t((grsR)_YAatPnQEx8XEsX4YDt5WW8diSuh?wjTFOWpoGlHAio z-}(Pk+CUl}>R%ns(9Z6E{s!uoe{UIa^75ckb>E4ESMZ?gL4q!7DM3C~)fOO62_eHQ zQ_^^IA=WHJOMqXFObR~Zn!`6M%h#vAd4+@NkR=YO9rZ$!iS3dqVa^u*tnt~?ft1bY zl(KQ78kLvA{!?`H3|kO3xpo1ArW>8e=n7K-g|ackYiX-&4%HAOJIigcpwb%X&w=jN z6JPj7FC0U)#q^Ebr+y&y?KSS*`pS+NcnElCoh(6?&t}&<&jE_H4?L}maun=czP#B> z$Y4^<)6R8|OmgK>cUUS-1PfLl(oo-kO?)IM(f%M%OSvGS<=nbAVHfghu>o`W)8I3h zVYG4xavsz;Tv>ADs6tB;*THvG{Sq+~B{sPJO$E8S>80nn3&EH$VMAM%fPy~Kt)&`a zkCG(RJD9Dz1UCa*ckv7_m$`S=X|^AF3a2&oS5C@(a|c`qEB(*N+j?@fW>KqU{7yrr`hF9rF?vJpOF7hfx*BTL1h1*m&N@t zOn)5z!y-8a>AxHJ`vRUnfIpAPpji0Ja-QFTzb|U|6WR*m^xu{^{0{zm4f#(f01yWE z3;h47DF1Hf_rl_zmf8^h?;-wCX8hgC?_&C&R#4D>S@~UB|J}gva=@PkqCn+NP+$Bh z4*U-NJ%|4Z#U%X~^!JSZyM@0e%Rlh|z%j_eKT_xK@V`^vU*UQbe}Vs-2o*;z=Ro8uw)Hz344gwMb01bcz005)_1AxDr|G(|O_z6^|3@diC041)aAH=p8WtZxOVL87C^`p_L z2!82F?khLZ%dxb4&WXGSisoZlan+%XEqSw@jaby%*){}4^tGxYMh6aj(J;X0V(#rf zq6A^%COPTqoaABP3$x-H8pgiO0%SNgcB^uTuZzf5;#pu~2?lItf9uyoUF~Z|EY~gL z8CdLGQ`B5RBrKpSH53^)|TJWR!u7Ci7mxax*AaCe>D>HCAnQx=WYqaHeTu zpLt^fzQPig(2l%3x(eOe!9}FNfmi<_>PNY5UxGC~$VcDD3q%ta{X_Y#WKHUD(5aD-AbzN7a+5Q@;SDU0)H7qO#oEK3x<&LZ})o`T3qIj zJs2Rp+Tf^+M#JDEZ*Z>)OTTk;g`uNzN|SV~T<<}3o4K61%9N7zq;YGHr>kfx%ah>wSDiQ{ZI+p~uYuDfF3bTBBIxNPOcL;X`ph8bMTa$Xb!>qUi7UY?AH4YRww zleMwEz4cGORiv(>w8#SV%l-5WycOs{w#-ANrDL3->kNAZ&77KW*c2;eY~COi4Se?9 zEuu~7V2MEE8|BFbwcZ~qD$bd$3pm(@B+-(=;CSk@rNTAKPn;Y^`p58_eE4dSPoXqH z&b8-qdO18OyMN!TdQAk7PA*vbMh6Gp?PQ?poJYshS+Cqugq7xnn<4A`F7v&!9+pUi zBSDI~^fVd8odAwi9WW*yvqFx^7NNjXgmIj!xlBG>v;O)^s`~0&>X^C?SdTtNm*A_s z0eT&_*RsX?ZGgtPqTMRaC02gOZlD`=pyf%KDoa(2?mR>QvhJk!6FLE~9zE#N4{;?N z1L~a+YbRPF)!i4`T5{+vZRUhk76s$g=owqWn%#F8u?z zClu!uKYE%uEpJeEQ8ZVdD^o~JQG)8Q)PVo%! zXk})$3SWU!nu?q`c0tEI>Zrkf@-cR|KTQhTt{cFwap~K2BP$?^cVFGZ>@&>C)ESiN z2DWcNoAA{~g#|^@)2&lQw(NvfY66MwiHVzWA&STG_zlvtdk>!@l(G0z?x@p{_X;_k zxP;X56W?~SQuy{TBobByrYK*#V+AGanlf|bobH@2KfK+u85=d4Te5>R%*1wVDsJq+ zWK*laxTSzBC2~)%=<}dw%L;C};|^l{-e|AS;QU73>#>6#ijem$?KdM1J6wge^X1xZ z8V`vNmW9GqnF0&axhWjxls~NJjI&6{m8g#q1xBQkdycsz5ujDwlNTfKP39(Ho?nRi zPbvorrleoI{L=D5XFLEb*b9~aW%B<@=s(Oq*h@6_q7eVxTcxUkTn`Jd9pNdQ#XZXn z6LrCbne0g85CwMd6YU}uDZAf~ReZXZ*V?mk%ng_9GV~5U|GItB8W6pwQ8d4h}26BOt$XMsF5O8equ5;NH2Rp~MX@EFlyB(L5Q( zVKvI)%4X#}PB`>p{)Rto1>A+G9aRN=YHTEVNBegBxmj%S3OYB)lX&6+7fyyB$rCWK zH%&@^(5=By0h>Oahsq*j;Nedz$eUwBz!=M*q=7ZdAlI9RpVscMf84% z7PNTIF#O}>?bcxb%<9WZ|C2I~wLN^>Pyhf2Iskz4a>l={pR?uN~`8Gx$Df z(5{OGL*hF14khh+jRKJlHUEG_-A$vk-O(T-G<^D1a`}^Yu2`U~jdEd+qs{J;`**&M zWh;rnU*AMvC*jMUR<`=F5}7iK?_MS&UE| z^E|%X`7+}B=UFHx`jt1OiN+n`sB@nB_wcD|6SVFp z6tp6+D=l{wn1$uWRd%95>;-T)aL8IlQD<=?zoM&tp_P>t5U34Aa0PmPxZl4{C;Xlh z=imOoyD??HBJe17cYb%f3tC-`tFZ^IfR2WGT;1nE4f8dyeM(Y3&)g}G7yU-h54?u= zZoaXjb2*M*ZyA=JpL2Url=eaArjk3oJ#t2oUGXR>i4fZ^Z?g3{qsnw45zL-$C@Q4B zCjv>|s^~YTC@-Vslask%C8|_OE1;Q7eDX#WKqHeD&z2JovDwOxleXA?;yx`Bu`Ej5 zO7mzT&Si2JT3BC1(;*ixLo=5YtPY8#snro-JeK@irXgpUs!h)=^-dV9R19p|F5(mw z)ekKZ{a8=73bzb~W1u~PXi%-nRQZS!DT+&+&(B$zl$>(-kMIn*~v6M7}8^o zawG}Uwr(Pvhm6!5&ZZd8T%zPGW}rt_H2x3`o~O{aY#3nH1;ktJ zd$7juJ5{o15tBe|54B#LM+?n-&bBbmZQW*zJG@WA}}p{UJu__q&oc00zfY16pcRAR~CDz$H2Q~8h(=;#9u2E^;i26`le zCxP^+%k1kIaNh7TRgyY3*^!~Csd^V4p=eQeO+uEiX%RV`rsdI|x^IZ|m5N4(cz`*= z=9;fm3X@~G)Dgjwd0~bzGQj1#9QeHr*ex_2p{X8l09YEO!yc`fvm7CE;@9F7IN}rl z02Ph%3g%{v932=}SSP;vt1Od<^l1@kH4#?TSTU94!$6kZbwqi2g`#o*M(x=wY$>TR zSxa$Hi6#Ws`xD3XP^ffxRdjYz>o{8z6N^(I6RH40gpEI>Zg);abw=U0Jae!mJ^KtB zte|l1$UP>oSHLm=6RW^n-upt5e59&A)w7~l%CP>~enp{okQ!APtI~DS?FEMJ;K#f* ze9|K7BJkswbIIFop@an#?hhj3;$|EjdF(FVLwVqJl)s{{K}eoDFI?S`}jD5ePm%#7QAlu6kF1e+LhbDS} zz(4^iJH}asS>ZU5dps!Xr1Ldd_(zIeqSBT~)7)LQaLYfsfux?K6eX?~V9?wKPE1Ex zdd@}^=kawVh{6d5*heK8Y(3KA*9>czFq`2FEH$TU^*n-^Dh=HZLkkoBtUG{L?jbVgoxAUo1`^$K>6^?X9qVc&fK-hOeh|ca0oSX5KD{Esm$;%X`*(n5{8;(>o+1ZlX z@X^63c`j9Hii*=k?FCXMg&~mIf5gL4eA*|TpdOhguiB=_vXRwwLyDGqYA-Fl&hhiQ za5A2 zFB>WO5mKtGzxRhWrFGUsQKT4}Rju+m%KW;`?t0OP&^^^-Leg*NMmel}L>u0d*EDEo zeipZskGt57%vckeRN~vvI%0Nt&uq=|Z+x^ToX__=ttLregLS=&4I#?(`!5e4bwLf2 zrC=5uL8&DT^SpUW1{zb!RCoXf*7)a#y;`);7z_`_ORT*j#>j4AZ9m3#9yOw+EtIP?|7V3GDD+A);` zq>3X|HF<3^u5I0zV>UrR0JtiJ%*J3ZqS82QeN!|I&sxtSzhB* z$=ShBs;M2%#2Y?8K4PyJf}Tzf&baocsn5E4Ko7UW>7eH`-pksJB~0cH-}CeRWb9Sn z%Uz{-!f%TNgc+x#k#(ol>n0~u(MEyD_h(>2**z?xJxRi-dv0aN!vo+3dL##3uC?F0 z774WbpxXFB)(b7fXU!CGnKhDiP)jB;ZaD6|Yvu2TPfs0k5#nj>2y#}loOnPe*g z*SqX_S{tehv9)G9NKt5!CET5pYM?r?7{ZEHv54(C0k#eK+Q=Wrz?$HQDM~S?H@HYfj1 zDXMB_3URL}+KW3Ah7h?O-tbF25fD7$Ch|Se>z-Fha!z}b&CxjQi2`o(vC90;-{g}9 zdisgl?;TpOg-;y;!ZCSlbp^msS8zro3i^*&!#&!v+ zjW+wn-!T7qB$Azv_u@n9w^kfAmzX{r??ypq?QtXa;oas5EwVSMgh#TK^aGq65ul275Wi6vEIcCi za!>^3lfHO%g|5wUO4fYVD-;L}NzHcjjd(#0=J67CB2_9Fi)XXFWBX4Q6hh@5|p#83@*xjap}3==NP+U5Z#`w)q4OQJtLc zU|32aXWVq85~61N7C(-t zD`vL*<$aX=rMd?@Rh~&vxGOzo?s(O`;)O?i?3AK;b7jlB4+v06w;r_mv0~}6C40I8 zCo5IchkkqN?mFqvGz(KWM)Gg*Ra|jY43R}K49r6e7BaObT(l>Y57P5>B^|fnS!^*e zJ+`mbB-*#bMZSrIy(LK5HNh6Og>Arnm5n3_BsRKNZP$lUrT2AG@DYhaWlRX`#~%UsZ+Le~Mh-ZeJp)GgJOM2-Jrd^Wv!}+$5D$d2G z>*C7wa`HJS6x?OhC|uxo4cWcOh{s0|gNZaA95kG1*-u5>7CZxU7Xk*5yt!1RsA zCEP4)`sbLAH@6u_LQ$im-s|1mG+;oaz@_flTb6afH`{^fq1T|8~foPWNK_2@6!ZN38fGLQQM zA%9FvEH{8JM$bongKF3Z&noauHCr$$j3%AXYn6Cz7dsGT>zIf%CYygPhdPF6j*K#x z*j%Wz$T88;FA;feQ>n*QQr)YQ!mQkYcL)-KYR|7u<^GCglw^LUz1my+WO=@cQiyhv zbhGG8V1h(a?^G#d($>Kl!{!w?RogZ==Likb*^$6AZeq}sY-ZYUP)U-sG@dV6;F<)Q zd4)P;NwSexYPZ?Uw$muhkxDwJkhB$Qx1>on=4$B}gj_)CFYDUjSn5?6&VOsLh#7D7 zc2mW`^vq;FLt4YTb&4$TSbe>E8p~oV9q{mqVwhcCZ|}argFML^7Vnb*>yLWbP721bHeS{1M{y-DU!x_v1uI-VD3oYQ<>{MEkqs z4ldLm6pS6LK|ze)uu1j8z2`!!eBZiQ9ReYuXM8GQz$}8z-oZ@Gk|J4~Jbi`q)?(j5 zs-9)1b_fzF6+mwqI4eY^WZi>&9d;hmxn}3!|JlLK{z&`8inL=EmIR`Dk%r8G_}Y^V zc3_Ot!)cqgsxW2(_frTh2Bi#NAs=#P+xn;}q%|=CwhjPW_wm-UG3*0huo)g)g^Gat z&b|ifF+-N3p%c;&zX-xcuazV=;$>~_C2`OK#OGm8s>5E?agq!qLBJ2AgRz!tm!139 z!Qa3hToAo>6LUvfQ9c|nlw1b5wmQ}KPO*?38ehC!yMG!#KnWLHb!C_<7Ri~5kz62K z>5HygIbo2$C-PVeJVAWyK#I#aPI%W+4uj0ia@v+}e%pZHUncI~2J^aqdP>7Q^+H%h zlLrOEEVY^?7C8k(MC-e)v(2xZqA>~~p+o3W`32%H$Hg;=$o+KcnHXag?_r|qhdlFW zDzX^`t)Z?`|B$qjtoaXZJ7fzUJ^XrOO z1*wuLkHA;^ekZz%Fc&aJ6mE4THcVIt>RrAqZs3Y;UKOd;09#oz|9 zUjInMyfnQZ_?&P^J||osuOA(!{b_*fbYssoDEN)VI~R8sRzqZnoVcBZSNCW9i(A6t zV-5R3TxwYl;~GU<^J}^TXJOVPgZS0GiyqDnYj6_L7iU;RLfevzmn_cWj!8b`x=*|H z*@Agathl!6x(zbH82(`eP;DT}la+Vz%0!hynrafRU}^((lCM4qElDHH9hzP@$?Q6! zj;FUwe)zbSBKWbe;H48$(J3N)9!0km#V5hZ}q#x#_^il=SO_iAw8t4MI2J#GuZ8h!&l118uc@MA?;=Fua1ht8g zI^O7}Dm4$$@<2O!uB$IHEyC0F!iZP1(a)Oc-Q)3ivdr3EP*;^MOQ5QYl0IOjNbEI< zz_p1|4}{=Nl!P>@P941U#oENc>Z$9hk6d!rCTUMM1`~~ym=5X=O7y11%BZQOS-`=< zsABE|hpMTINORWCt4P0tQm|rRy~;5?S{CyxiY8(*NO*$_@?c1lVl#Kr9^aO$(NFft^tQ@0oR%u}&ujd;upxcfSyTnbpI01hL#{@#S$HOhtb zo-vb@)jO^b+ZNhhO^?KpKfzKRde5j;v5TR{vu-M9gVAoir^q|WL1@<#pBy#NhRBrQ z4hng$tNC@+YyAj<2Br|p{vaG8)qvKXEH{>j2l7%i=0HoSV48B%M+%4w9zFYSh>^Oo zC3d(G%rxGS*ksNHZ9M+i|W|T z5RGxY$|TYXQ$aLb1vKBPcBPN38QQAl)jgebyR5AWLAu#nD25BnMYiXhD4clb@MW#6 z!gN?>eSB^bJ3ZkBR~kF&z8ZA0ABIF8!ge*c-UZcpn`F0^EtRf9gExAv`nKAyc&05N z5KJ93KRV3gUlG<|d!+fUx|gCDtT}dg+w;@eE>FC8qW!%1-tU&detEPG&zttd|P;|4Jk;D;|^h+P;efC-fx!NQh&p z@rQ|;n6CWJ#*F;SP|pWmBdbXRrhmDLV*P!ozn5wK{q=fVxRaN%>kqKmK}jtdIaQ3n zo!A5~z5MhU{+zr4&l0E#d?ZZXATc#BZ#5?|4-+Fxe27^&2Ne=asiqI(quydHZay$k z%@o`&K(1Jgui=Hu#@#o3i{PVEE+Lp6R2k+skTD~|;sUCtV3t=6$ABP{GLgV5WOJD? zP#aODlRvV{P^y}YKZjAN;4!D(=|}D+6h-8|2HWI@J>!rET))@+0`3)ks<|y20ec^E z-|>K8HS8!N!pUFP`x1`0y;_nTk$qTHfnPDBSOUKwDLAF$F>j+eG`?zjNB*wV6SdUy zeY}>ZY!oZ4q&?*ms}fA!fyqNIH|e#zRu4jeTQlgH|M8^_>Yw?vub1!716~;S^ujpc zUm0iQ;P5}fy)f8-P70+RKe=>^(bWMcNBW{(XR@mZ^-*kt6;hJgq~a&=N@ z>398+y1nB@rqL;s2uj5k8Lh=Brg%1iSeDvG6j*8)z0%A*{g@p4 z_vkAXTZ8j62+{?fYb-Ipevwkj@=Xu;+;YHcj^xl%{A8klh+#Eq=MJ5nbWtqpV$D%A zac;0HJF55j)!Hd9?jasQeiAW+YdA|t1311+MNYZ^9x^NWJ#()cRXs?Udo86P+@jOD z3(r)IeYToB?vAGi+X|SOA|GVQo#y<#tpVi13a{$OK!8Y<@QGV;N!v{tdl#*)h5UV}eVY^Se)4tuCG(GY+b$%Q4s&;Y^+RYfP1<9zhVZ}^xjG5141KP( zM(uX%4yLelzjqxUnMSd1$M^*{Y)ZEJ`)BN5p8vmC3kJ^kV(kBObKf7U^vC``wE8K_ z{awJ{+mQYM{@kZS004h!O8Oo6duzs@(6$$w{#&EQ@8G`|!~cW=0O4@I!2h4p`0sLl zuaN#JsRiNx-NZlYrN2x0J(>Qelo`}tQhraZ|1RM7bikhi)Ny|a_%k{1JM?!S{wMSs z>A#@Ad-UHW{M}jpi3b2ak^un!aGk%y|8DvI3a6y}3;f@9P+1P@Wpx1n_?J(>i>+Lw H{(1C&Fg6C8 literal 0 HcmV?d00001 diff --git a/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index 5a239dd..1901329 100644 --- a/electrode_models/air_electrode.py +++ b/electrode_models/air_electrode.py @@ -93,7 +93,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**1.5 #This is used by 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) diff --git a/electrode_models/dense_electrode.py b/electrode_models/dense_electrode.py index 9d89932..3380b48 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'] + self.eps_elyte = np.array([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..37d7211 --- /dev/null +++ b/inputs/LiO2_Liu_starter.yaml @@ -0,0 +1,516 @@ +#=============================================================================== +# 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: 1.42e-9 # m2 S-1 // Saito, 2017 + - species: 'TFSI-[elyt]' + D_k: 8.79e-11 # m2 S-1 // Saito, 2017 + - species: 'O2(e)' + D_k: 7.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: 5 + 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: 245.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. + # 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 + # - 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 + # - 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 + # - 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 + +fit-parameters: + - type: 'elyte-transport' + species: 'O2(e)' + guess-value: 7.30E-10 + min: 6.30e-10 + max: 10.30e-10 + - type: 'elyte-transport' + species: 'Li+[elyt]' + guess-value: 5.42e-10 + min: 1.42e-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: 245.e-9 + # min: 300.e-9 + # max: 145.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..609aaa7 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 @@ -210,7 +211,7 @@ def electrode_boundary_flux(self, SV, ed, T): state_1 = {'C_k': C_k_1, 'phi':phi_1, 'T':T, 'dy':self.dy, 'microstructure':self.elyte_microstructure} state_2 = {'C_k': C_k_2, 'phi':phi_2, 'T':T, 'dy':ed.dy_elyte, - 'microstructure':ed.elyte_microstructure} + 'microstructure':ed.elyte_microstructure[j_ed]} # Multiply by ed.i_ext_flag: fluxes are out of the anode, into the cathode. N_k_elyte, i_io = tuple(x*ed.i_ext_flag diff --git a/submodels/bandwidth.py b/submodels/bandwidth.py index 13ced6c..6c7634e 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 From 4516ff5a8c8adf93b825808ce6a05f07c0ed89b9 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Mon, 27 Mar 2023 12:01:26 -0600 Subject: [PATCH 09/14] troubleshooting --- electrode_models/air_electrode.py | 5 +++-- electrode_models/dense_electrode.py | 2 +- separator_models/porous_separator.py | 4 +++- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index 1901329..733721e 100644 --- a/electrode_models/air_electrode.py +++ b/electrode_models/air_electrode.py @@ -93,7 +93,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**1.5 #This is used by the + self.elyte_microstructure = self.eps_elyte_init**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) @@ -212,13 +212,14 @@ def residual(self, t, SV, SVdot, sep, counter, params): # Start at the separator boundary: j = self.nodes[0] - + print(self.eps_host) # Read out properties: phi_ed = SV_loc[SVptr['phi_ed'][j]] phi_elyte = phi_ed + SV_loc[SVptr['phi_dl'][j]] 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 diff --git a/electrode_models/dense_electrode.py b/electrode_models/dense_electrode.py index 3380b48..2eaa3be 100755 --- a/electrode_models/dense_electrode.py +++ b/electrode_models/dense_electrode.py @@ -56,7 +56,7 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, self.dy_elyte = inputs['dy_elyte'] # Electrolyte volume fraction in the separator: - self.eps_elyte = np.array([sep_inputs['eps_electrolyte']]) + self.eps_elyte = sep_inputs['eps_electrolyte'] # Microstructure-based transport scaling factor, based on Bruggeman # coefficient of -0.5: diff --git a/separator_models/porous_separator.py b/separator_models/porous_separator.py index 609aaa7..c026acf 100644 --- a/separator_models/porous_separator.py +++ b/separator_models/porous_separator.py @@ -210,8 +210,10 @@ def electrode_boundary_flux(self, SV, ed, T): # Create dictionaries to pass to the transport function: state_1 = {'C_k': C_k_1, 'phi':phi_1, 'T':T, 'dy':self.dy, 'microstructure':self.elyte_microstructure} + print(j_ed) + print(ed.elyte_microstructure) state_2 = {'C_k': C_k_2, 'phi':phi_2, 'T':T, 'dy':ed.dy_elyte, - 'microstructure':ed.elyte_microstructure[j_ed]} + 'microstructure':ed.elyte_microstructure} # Multiply by ed.i_ext_flag: fluxes are out of the anode, into the cathode. N_k_elyte, i_io = tuple(x*ed.i_ext_flag From f4cc52be30e04f7f8abc9ff4ae862e72ad6e4fdb Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Wed, 29 Mar 2023 08:16:30 -0600 Subject: [PATCH 10/14] prelim changes to fix cathode --- electrode_models/air_electrode.py | 11 ++- inputs/LiO2_Liu_starter.yaml | 130 +++++++++++++-------------- separator_models/porous_separator.py | 2 - 3 files changed, 72 insertions(+), 71 deletions(-) diff --git a/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index 733721e..635f067 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: @@ -93,8 +95,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**1.5 #This is used by the electrode_boundary_flux - + 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 @@ -212,7 +213,7 @@ def residual(self, t, SV, SVdot, sep, counter, params): # Start at the separator boundary: j = self.nodes[0] - print(self.eps_host) + # Read out properties: phi_ed = SV_loc[SVptr['phi_ed'][j]] phi_elyte = phi_ed + SV_loc[SVptr['phi_dl'][j]] @@ -220,7 +221,6 @@ def residual(self, t, SV, SVdot, sep, counter, params): 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 @@ -414,6 +414,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/inputs/LiO2_Liu_starter.yaml b/inputs/LiO2_Liu_starter.yaml index 37d7211..0e71793 100644 --- a/inputs/LiO2_Liu_starter.yaml +++ b/inputs/LiO2_Liu_starter.yaml @@ -25,7 +25,7 @@ cell-description: separator: class: 'porous_separator' thickness: 6.5e-4 # separator thickness, m - n_points: 5 # Number of finite volumes to discretize + n_points: 20 # 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 @@ -50,7 +50,7 @@ cell-description: # D_k: 2e-13 # can't find (assume low?) cathode: class: 'air_electrode' - n-points: 5 + n-points: 35 gas-phase: 'air' elyte-iphase: 'air-elyte-interface' electrolyte-phase: 'electrolyte' # Cantera phase name @@ -67,7 +67,7 @@ cell-description: th_product: 245.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_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. @@ -106,82 +106,82 @@ parameters: # them (False)? save-name: '0p1_liu' # Folder label for output files. validation: 'Liu_O2/0.1mAhcm2.xlsx' #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: 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 - # - 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 - # - 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 + - 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 + - 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 + - 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 fit-parameters: - type: 'elyte-transport' species: 'O2(e)' - guess-value: 7.30E-10 - min: 6.30e-10 + guess-value: 6.30E-10 + min: 5.00e-10 max: 10.30e-10 - type: 'elyte-transport' species: 'Li+[elyt]' guess-value: 5.42e-10 - min: 1.42e-10 + 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: 245.e-9 - # min: 300.e-9 - # max: 145.e-9 + - type: 'cathode-microstructure' + parameter: 'th_product' + guess-value: 245.e-9 + min: 250.e-9 + max: 105.e-9 diff --git a/separator_models/porous_separator.py b/separator_models/porous_separator.py index c026acf..0206e6b 100644 --- a/separator_models/porous_separator.py +++ b/separator_models/porous_separator.py @@ -210,8 +210,6 @@ def electrode_boundary_flux(self, SV, ed, T): # Create dictionaries to pass to the transport function: state_1 = {'C_k': C_k_1, 'phi':phi_1, 'T':T, 'dy':self.dy, 'microstructure':self.elyte_microstructure} - print(j_ed) - print(ed.elyte_microstructure) state_2 = {'C_k': C_k_2, 'phi':phi_2, 'T':T, 'dy':ed.dy_elyte, 'microstructure':ed.elyte_microstructure} From a113da5230caeefb432276796a6d47da0d93b753 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Sun, 9 Apr 2023 10:58:47 -0600 Subject: [PATCH 11/14] allow running separate radius --- bat_can.py | 1 - electrode_models/air_electrode.py | 12 ++++++++-- inputs/LiO2_Liu_starter.yaml | 38 +++++++++++++++++-------------- simulations/CC_cycle.py | 4 +++- submodels/bandwidth.py | 6 ++--- 5 files changed, 37 insertions(+), 24 deletions(-) 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/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index 635f067..381fd3c 100644 --- a/electrode_models/air_electrode.py +++ b/electrode_models/air_electrode.py @@ -57,6 +57,12 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, eps_host = np.array([inputs['eps_host']]) eps_host = eps_host[0] + # current = params['i_ext'] + # # Split the current from the units: + # i_ext, units = current.split() + self.surf_to_volume = 1.23224E-11 + + if len(eps_host) == 1: self.eps_host = np.repeat(eps_host[0], self.n_points) elif len(eps_host) == self.n_points: @@ -73,6 +79,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] @@ -271,7 +279,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 @@ -364,7 +372,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 diff --git a/inputs/LiO2_Liu_starter.yaml b/inputs/LiO2_Liu_starter.yaml index 0e71793..4be439f 100644 --- a/inputs/LiO2_Liu_starter.yaml +++ b/inputs/LiO2_Liu_starter.yaml @@ -41,16 +41,16 @@ cell-description: - species: 'C2H6O[elyt]' D_k: 1.43e-10 # m2 S-1 // Saito, 2017 - species: 'Li+[elyt]' - D_k: 1.42e-9 # m2 S-1 // Saito, 2017 + D_k: 6.42e-11 # m2 S-1 // Saito, 2017 - species: 'TFSI-[elyt]' - D_k: 8.79e-11 # m2 S-1 // Saito, 2017 + D_k: 6.79e-10 # m2 S-1 // Saito, 2017 - species: 'O2(e)' - D_k: 7.30E-10 #11 # m2 s-1 // Gittleson, 2017 for DMSO + 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: 35 + n-points: 50 gas-phase: 'air' elyte-iphase: 'air-elyte-interface' electrolyte-phase: 'electrolyte' # Cantera phase name @@ -60,13 +60,13 @@ cell-description: 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 + r_host: 0.75e-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. + host_surf_frac: 0.98 # Fraction of host surface area that is accessible. d_oxide: 1.0e-8 #Oxide diameter, m - th_product: 245.e-9 #185e-9 #oxide thickness, 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 ondutivity of the host + 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 @@ -106,8 +106,9 @@ parameters: # 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: + # 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. @@ -124,6 +125,7 @@ parameters: # 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) @@ -142,6 +144,7 @@ parameters: # 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) @@ -160,16 +163,17 @@ parameters: # 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: 6.30E-10 + guess-value: 24.30E-10 min: 5.00e-10 - max: 10.30e-10 + max: 20.30e-10 - type: 'elyte-transport' species: 'Li+[elyt]' - guess-value: 5.42e-10 + guess-value: 6.42e-11 min: 3.32e-10 max: 1.4e-9 # - type: 'cathode-kinetics' @@ -177,11 +181,11 @@ fit-parameters: # 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: 245.e-9 - min: 250.e-9 - max: 105.e-9 + # - type: 'cathode-microstructure' + # parameter: 'th_product' + # guess-value: 100.e-9 + # min: 60.e-9 + # max: 105.e-9 diff --git a/simulations/CC_cycle.py b/simulations/CC_cycle.py index 60ebb92..d748a4c 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 @@ -368,7 +370,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 6c7634e..d6bf24f 100644 --- a/submodels/bandwidth.py +++ b/submodels/bandwidth.py @@ -31,9 +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) + # 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) + # print('i > j', i, j, uband, lband) + # print(an.SVptr, sep.SVptr, ca.SVptr) return lband, uband From 59fd7c1309f8aa4ffd826982d7970fee938e0bc2 Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Tue, 18 Apr 2023 09:32:27 -0600 Subject: [PATCH 12/14] energy and power density --- inputs/LiO2_Liu_starter.yaml | 80 ++++++++++++++++++------------------ simulations/CC_cycle.py | 21 +++++++--- 2 files changed, 55 insertions(+), 46 deletions(-) diff --git a/inputs/LiO2_Liu_starter.yaml b/inputs/LiO2_Liu_starter.yaml index 4be439f..c42f30f 100644 --- a/inputs/LiO2_Liu_starter.yaml +++ b/inputs/LiO2_Liu_starter.yaml @@ -25,7 +25,7 @@ cell-description: separator: class: 'porous_separator' thickness: 6.5e-4 # separator thickness, m - n_points: 20 # Number of finite volumes to discretize + 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 @@ -50,7 +50,7 @@ cell-description: # D_k: 2e-13 # can't find (assume low?) cathode: class: 'air_electrode' - n-points: 50 + n-points: 20 gas-phase: 'air' elyte-iphase: 'air-elyte-interface' electrolyte-phase: 'electrolyte' # Cantera phase name @@ -126,44 +126,44 @@ parameters: 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 + # - 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' diff --git a/simulations/CC_cycle.py b/simulations/CC_cycle.py index d748a4c..a5f4c8b 100644 --- a/simulations/CC_cycle.py +++ b/simulations/CC_cycle.py @@ -93,6 +93,15 @@ def terminate_check(t, SV, SVdot, return_val, inputs): 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 + voltage = solution.values.y.T[-7] + power_density = voltage*i_data + energy_density = [0] + for item, value in enumerate(cycle_capacity): + if item == 0: + continue + else: + delta_capacity = value - cycle_capacity[item-1] + energy_density.append(voltage[item]*delta_capacity +energy_density[item-1]) # Append the current data array to any preexisting data, for output. # If this is the first step, create the output data array. @@ -100,7 +109,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: @@ -109,7 +118,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: @@ -249,14 +258,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], @@ -310,7 +319,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]) From e6201c68e1b7089274f23d224ec00c2a41276cfc Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Wed, 19 Apr 2023 11:19:00 -0600 Subject: [PATCH 13/14] more thought out power & energy density implement --- simulations/CC_cycle.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/simulations/CC_cycle.py b/simulations/CC_cycle.py index a5f4c8b..f30969f 100644 --- a/simulations/CC_cycle.py +++ b/simulations/CC_cycle.py @@ -92,16 +92,18 @@ 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 - voltage = solution.values.y.T[-7] - power_density = voltage*i_data - energy_density = [0] - for item, value in enumerate(cycle_capacity): + cycle_capacity = 0.1*solution.values.t*abs(i_data)/3600 #mAh/cm2 + + phi_ptr = SV_offset + 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: - continue + energy_density = [0] + power_density = currents[i]*voltage_value else: - delta_capacity = value - cycle_capacity[item-1] - energy_density.append(voltage[item]*delta_capacity +energy_density[item-1]) + energy_density.append(delta_t[item-1]*currents[i]*voltage_value+energy_density[item-1])/3600 #Wh/cm2 + power_density.append = 3600*energy_density[item]/sum(delta_t[:item]) #W/cm2 # Append the current data array to any preexisting data, for output. # If this is the first step, create the output data array. From f1fa384dff06d629c9514dd85f979901a953ae1b Mon Sep 17 00:00:00 2001 From: mglasser16 Date: Wed, 19 Apr 2023 12:18:53 -0600 Subject: [PATCH 14/14] unit fix --- electrode_models/air_electrode.py | 6 ------ simulations/CC_cycle.py | 9 ++++----- 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/electrode_models/air_electrode.py b/electrode_models/air_electrode.py index 381fd3c..c51eb0d 100644 --- a/electrode_models/air_electrode.py +++ b/electrode_models/air_electrode.py @@ -57,12 +57,6 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs, eps_host = np.array([inputs['eps_host']]) eps_host = eps_host[0] - # current = params['i_ext'] - # # Split the current from the units: - # i_ext, units = current.split() - self.surf_to_volume = 1.23224E-11 - - if len(eps_host) == 1: self.eps_host = np.repeat(eps_host[0], self.n_points) elif len(eps_host) == self.n_points: diff --git a/simulations/CC_cycle.py b/simulations/CC_cycle.py index f30969f..5286a25 100644 --- a/simulations/CC_cycle.py +++ b/simulations/CC_cycle.py @@ -93,17 +93,16 @@ def terminate_check(t, SV, SVdot, return_val, inputs): 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 #mAh/cm2 - - phi_ptr = SV_offset + ca.SV_offset+int(ca.SVptr['phi_ed'][-1]) + 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 + power_density = [currents[i]*voltage_value] else: - energy_density.append(delta_t[item-1]*currents[i]*voltage_value+energy_density[item-1])/3600 #Wh/cm2 - power_density.append = 3600*energy_density[item]/sum(delta_t[:item]) #W/cm2 + 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.