diff --git a/bin/mpetplot.py b/bin/mpetplot.py index a010fd37..7f10e3b0 100755 --- a/bin/mpetplot.py +++ b/bin/mpetplot.py @@ -36,6 +36,8 @@ ('cbar_a','average anode solid concentrations (movie)'), ('bulkp_c','macroscopic cathode solid phase potential(movie)'), ('bulkp_a','macroscopic anode solid phase potential (movie)'), + ('temp','temperature of full cell (movie)'), + ('max_temp','maximum temperature of the full cell'), ('text','convert the output to plain text (csv)') ]) diff --git a/bin/run_tests.py b/bin/run_tests.py index 763b3e39..91d55627 100644 --- a/bin/run_tests.py +++ b/bin/run_tests.py @@ -23,8 +23,7 @@ def run(test_outputs, testDir, tests=None): # Generate a list of tests from the directories in ref_outputs ref_outputs = osp.join(dirDict["suite"],"ref_outputs") _, directories, _ = next(walk(osp.join(ref_outputs))) - tests = directories - tests.sort() + tests = sorted(directories) runInfoAnalyt = { "testAnalytCylDifn": (defs.testAnalytCylDifn, defs.analytCylDifn), diff --git a/configs/params_system.cfg b/configs/params_system.cfg index ef5dca6e..d0b49ea6 100644 --- a/configs/params_system.cfg +++ b/configs/params_system.cfg @@ -49,8 +49,11 @@ tsteps = 200 relTol = 1e-6 # Absolute Tolerance absTol = 1e-6 -# Temperature, K +# Initial Temperature throughout electrode, K T = 298 +# Nonisothermal: true for heat generation throughout electrode, false +# for no heat generation +nonisothermal = true # Random seed. Set to true to give a random seed in the simulation # (affects noise, particle size distribution). Set to true exactly # reproducible results -- useful for testing. @@ -127,6 +130,26 @@ G_stddev_c = 0 G_mean_a = 1e-14 G_stddev_a = 0 +[Thermal Parameters] +# Heat capacity of anode, cathode, electrolyte material in J/(kg*K). +cp_c = 700 +cp_s = 700 +cp_a = 700 +# Mass density of anode, cathode, and electrolyte materialin kg/m^3. +rhom_c = 2500 +rhom_s = 1100 +rhom_a = 2500 +# Heat transfer coefficient with the cell separator (W/(m^2*K)) +h_h = 2e4 +# Thermal conductivity in battery cell (only used for dilute electrolyte model), (W/(m*K)). +# For Stefan-Maxwell concentrated electrolyte, input from props_elyte.py +k_h_c = 2.1 +k_h_a = 1.7 +k_h_s = 0.16 +# Includes temperature dependence of entropic heat generation if true, does not include +# if false +entropy_heat_gen = False + [Geometry] # Thicknesses, m L_c = 50e-6 diff --git a/configs/params_system_Fuller94.cfg b/configs/params_system_Fuller94.cfg index 8f111b6f..a9295151 100644 --- a/configs/params_system_Fuller94.cfg +++ b/configs/params_system_Fuller94.cfg @@ -6,7 +6,7 @@ profileType = CC # Battery (dis)charge c-rate (only used for CC), number of capacities / hr Crate = 1 -#Optional nominal 1C current density for the cell, A/m^2 +# Optional nominal 1C current density for the cell, A/m^2 1C_current_density = 40. # Voltage cutoffs, V Vmax = 5 diff --git a/configs/params_system_LIONSIMBA.cfg b/configs/params_system_LIONSIMBA.cfg index baa598b5..a0c4811c 100644 --- a/configs/params_system_LIONSIMBA.cfg +++ b/configs/params_system_LIONSIMBA.cfg @@ -102,7 +102,7 @@ num = 1 # Options: dilute, SM elyteModelType = SM # Stefan-Maxwell property set, see props_elyte.py file -SMset = LIONSIMBA_nonisothermal +SMset = LIONSIMBA_isothermal # Reference electrode (defining the electrolyte potential) information: # number of electrons transfered in the reaction, 1 for Li/Li+ n = 1 diff --git a/configs/params_system_LIONSIMBA_nonisothermal.cfg b/configs/params_system_LIONSIMBA_nonisothermal.cfg new file mode 100644 index 00000000..f5b4eed1 --- /dev/null +++ b/configs/params_system_LIONSIMBA_nonisothermal.cfg @@ -0,0 +1,127 @@ +#Cell parameters for nonisothermal benchmark from M. Torchio et al., J. Electrochem. Soc. 163, A1192 (2016). +# See params_system.cfg for parameter explanations. + +[Sim Params] +# Constant voltage or current or segments of one of them +# Options: CV, CC, CCsegments, CVsegments +profileType = CC +# Battery (dis)charge c-rate (only used for CC), number of capacities / hr +# (positive for discharge, negative for charge) +Crate = 1 +#Optional nominal 1C current density for the cell, A/m^2 +1C_current_density = 30 +# Voltage cutoffs, V +Vmax = 5 +Vmin = 2.5 +# Final time (only used for CV), [s] +tend = 1.2e3 +# Number disc. in time +tsteps = 200 +# Numerical Tolerances +relTol = 1e-6 +absTol = 1e-6 +# Temperature, K +T = 298 +nonisothermal = true +# Random seed. Set to true to give a random seed in the simulation +randomSeed = false +# Value of the random seed, must be an integer +seed = 0 +# Series resistance, [Ohm m^2] +Rser = 0. +# Cathode, anode, and separator numer disc. in x direction (volumes in electrodes) +Nvol_c = 10 +Nvol_s = 10 +Nvol_a = 10 +# Number of particles per volume for cathode and anode +Npart_c = 1 +Npart_a = 1 + +[Electrodes] +cathode = params_LiCoO2_LIONSIMBA.cfg +anode = params_LiC6_LIONSIMBA.cfg +# Rate constant of the Li foil electrode, A/m^2 +# Used only if Nvol_a = 0 +k0_foil = 1e0 +# Film resistance on the Li foil, Ohm m^2 +Rfilm_foil = 0e-0 + +[Particles] +# electrode particle size distribution info, m +mean_c = 2e-6 +stddev_c = 0 +mean_a = 2e-6 +stddev_a = 0 +# Initial electrode filling fractions +cs0_c = 0.4995 +cs0_a = 0.8551 + +[Conductivity] +# Simulate bulk cathode conductivity (Ohm's Law)? +simBulkCond_c = true +simBulkCond_a = true +# Dimensional conductivity (used if simBulkCond = true), S/m +sigma_s_c = 412.43 +sigma_s_a = 685.77 +# Simulate particle connectivity losses (Ohm's Law)? +simPartCond_c = false +simPartCond_a = false +# Conductance between particles, S = 1/Ohm +G_mean_c = 1e-14 +G_stddev_c = 0 +G_mean_a = 1e-14 +G_stddev_a = 0 + +[Geometry] +# Thicknesses, m +L_c = 8e-5 +L_a = 8.8e-5 +L_s = 2.5e-5 +# Volume loading percents of active material (volume fraction of solid +# that is active material) +P_L_c = 0.9593 +P_L_a = 0.9367 +# Porosities (liquid volume fraction in each region) +poros_c = 0.385 +poros_a = 0.485 +poros_s = 0.724 +# Bruggeman exponent (tortuosity = porosity^bruggExp) +BruggExp_c = -3 +BruggExp_a = -3 +BruggExp_s = -3 + +[Thermal Parameters] +# Heat capacity of anode, cathode, electrolyte material in J/(kg*K). +cp_c = 700 +cp_s = 700 +cp_a = 700 +# Mass density of anode, cathode, and electrolyte materialin kg/m^3. +rhom_c = 2500 +rhom_s = 1100 +rhom_a = 2500 +# Heat transfer coefficient with the cell separator (W/(m^2*K)) +h_h = 1 +k_h_c = 2.1 +k_h_a = 1.7 +k_h_s = 0.16 +entropy_heat_gen = False + +[Electrolyte] +# Initial electrolyte conc., mol/m^3 +c0 = 1000 +# Cation/anion charge number (e.g. 2, -1 for CaCl_2) +zp = 1 +zm = -1 +# Cation/anion dissociation number (e.g. 1, 2 for CaCl_2) +nup = 1 +num = 1 +# Electrolyte model, +# Options: dilute, SM +elyteModelType = SM +# Stefan-Maxwell property set, see props_elyte.py file +SMset = LIONSIMBA_nonisothermal +# Reference electrode (defining the electrolyte potential) information: +# number of electrons transfered in the reaction, 1 for Li/Li+ +n = 1 +# Stoichiometric coefficient of cation, -1 for Li/Li+ +sp = -1 diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 2f959267..52d92289 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -111,9 +111,14 @@ def _init_from_dicts(self): self.D_s = ParameterSet(None, 'system', self.path) # set which electrodes there are based on which dict files exist trodes = ['c'] + # set up types of materials (cathode, anode electrolyte) for thermal parameters + materials = ['c', 'l'] if os.path.isfile(os.path.join(self.path, 'input_dict_anode.p')): trodes.append('a') + materials.append('a') self['trodes'] = trodes + self['materials'] = materials + # create empty electrode parametersets self.D_c = ParameterSet(None, 'electrode', self.path) if 'a' in self['trodes']: @@ -137,9 +142,13 @@ def _init_from_cfg(self, paramfile): self.D_s = ParameterSet(paramfile, 'system', self.path) # the anode and separator are optional: only if there are volumes to simulate trodes = ['c'] + # set up types of materials (cathode, anode electrolyte) for thermal parameters + materials = ['c', 'l'] if self.D_s['Nvol_a'] > 0: trodes.append('a') + materials.append('a') self['trodes'] = trodes + self['materials'] = materials # to check for separator, directly access underlying dict of system config; # self['Nvol']['s'] would not work because that requires have_separator to # be defined already @@ -471,6 +480,7 @@ def _scale_system_parameters(self, theoretical_1C_current): # non-dimensional scalings self['T'] = self['T'] / constants.T_ref self['Rser'] = self['Rser'] / self['Rser_ref'] + self['h_h'] = self['h_h'] * self['L_ref'] / self['k_h_ref'] if self['Dp'] is not None: self['Dp'] = self['Dp'] / self['D_ref'] if self['Dm'] is not None: @@ -496,6 +506,17 @@ def _scale_electrode_parameters(self): self['L'][trode] = self['L'][trode] / self['L_ref'] self['beta'][trode] = self[trode, 'csmax'] / constants.c_ref self['sigma_s'][trode] = self['sigma_s'][trode] / self['sigma_s_ref'] + self['k_h'][trode] = self['k_h'][trode] / self['k_h_ref'] + self['cp'][trode] = self['cp'][trode] / \ + (self['k_h_ref'] * self['t_ref'] / (self['rho_ref']*self['L_ref']**2)) + self['rhom'][trode] = self['rhom'][trode] / self['rho_ref'] + if self['nonisothermal']: + if self['rhom'][trode] == 0 or self['cp'][trode] == 0 or self['k_h'][trode] == 0: + raise Exception( + "Please provide all nonisothermal parameters for the " + + trode + + " electrode and ensure they are nonzero") + if self[trode, 'lambda'] is not None: self[trode, 'lambda'] = self[trode, 'lambda'] / kT if self[trode, 'B'] is not None: @@ -508,6 +529,15 @@ def _scale_electrode_parameters(self): # scalings on separator if self['have_separator']: self['L']['s'] /= self['L_ref'] + self['k_h']['s'] = self['k_h']['s'] / self['k_h_ref'] + self['cp']['s'] = self['cp']['s'] / \ + (self['k_h_ref'] * self['t_ref'] / (self['rho_ref']*self['L_ref']**2)) + self['rhom']['s'] = self['rhom']['s'] / self['rho_ref'] + if self['nonisothermal']: + if self['rhom']['s'] == 0 or self['cp']['s'] == 0 or self['k_h']['s'] == 0: + raise Exception( + "Please provide all nonisothermal parameters for the " + + " separator and ensure they are nonzero") def _scale_macroscopic_parameters(self, Vref): """ diff --git a/mpet/config/constants.py b/mpet/config/constants.py index ed309c15..b89f122b 100644 --- a/mpet/config/constants.py +++ b/mpet/config/constants.py @@ -22,9 +22,9 @@ #: parameter that are defined per electrode with a ``_{electrode}`` suffix PARAMS_PER_TRODE = ['Nvol', 'Npart', 'mean', 'stddev', 'cs0', 'simBulkCond', 'sigma_s', 'simPartCond', 'G_mean', 'G_stddev', 'L', 'P_L', 'poros', 'BruggExp', - 'specified_psd'] + 'specified_psd', 'rhom', 'cp', 'k_h'] #: subset of ``PARAMS_PER_TRODE``` that is defined for the separator as well -PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp'] +PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp', 'k_h', 'cp', 'rhom'] #: parameters that are defined for each particle, and their type PARAMS_PARTICLE = {'N': int, 'kappa': float, 'beta_s': float, 'D': float, 'k0': float, 'Rfilm': float, 'delta_L': float, 'Omega_a': float, 'E_D': float, diff --git a/mpet/config/derived_values.py b/mpet/config/derived_values.py index fb18dd11..b111922b 100644 --- a/mpet/config/derived_values.py +++ b/mpet/config/derived_values.py @@ -148,6 +148,12 @@ def curr_ref(self): """ return 3600. / self.config['t_ref'] + def rho_ref(self): + """Reference mass density used for energy balances + """ + m_ref = 1000 # reference is 1000 kg + return m_ref + def sigma_s_ref(self): """Reference conductivity """ @@ -197,6 +203,12 @@ def D_ref(self): mpet_module=f"mpet.electrolyte.{SMset}") return elyte_function()[-1] + def k_h_ref(self): + """Reference heat transfer coefficient + """ + return constants.c_ref * constants.k * constants.N_A * \ + self.config['L_ref']**2 / (self.config['t_ref']) + def z(self): """Electrode capacity ratio """ @@ -232,9 +244,9 @@ def muR_ref(self, trode): solidType = self.config[trode, 'type'] if solidType in constants.two_var_types: - muR_ref = -muRfunc((cs0, cs0), (cs0bar, cs0bar), 0.)[0][0] + muR_ref = -muRfunc((cs0, cs0), (cs0bar, cs0bar), self.config['T'], 0.)[0][0] elif solidType in constants.one_var_types: - muR_ref = -muRfunc(cs0, cs0bar, 0.)[0] + muR_ref = -muRfunc(cs0, cs0bar, self.config['T'], 0.)[0] else: raise ValueError(f'Unknown solid type: {solidType}') return muR_ref diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py index 17e2ea76..16ed34eb 100644 --- a/mpet/config/schemas.py +++ b/mpet/config/schemas.py @@ -69,6 +69,7 @@ def tobool(value): 'relTol': And(Use(float), lambda x: x > 0), 'absTol': And(Use(float), lambda x: x > 0), 'T': Use(float), + Optional('nonisothermal', default=False): Use(tobool), 'randomSeed': Use(tobool), Optional('seed'): And(Use(int), lambda x: x >= 0), Optional('dataReporter', default='mat'): str, @@ -113,6 +114,17 @@ def tobool(value): 'BruggExp_c': Use(float), 'BruggExp_a': Use(float), 'BruggExp_s': Use(float)}, + 'Thermal Parameters': {Optional('cp_c', default=0): Use(float), + Optional('cp_a', default=0): Use(float), + Optional('cp_s', default=0): Use(float), + Optional('rhom_c', default=0): Use(float), + Optional('rhom_a', default=0): Use(float), + Optional('rhom_s', default=0): Use(float), + Optional('k_h_c', default=0): Use(float), + Optional('k_h_a', default=0): Use(float), + Optional('k_h_s', default=0): Use(float), + Optional('h_h', default=0): Use(float), + Optional('entropy_heat_gen', default=False): Use(tobool)}, 'Electrolyte': {'c0': Use(float), 'zp': Use(int), 'zm': And(Use(int), lambda x: x < 0), diff --git a/mpet/daeVariableTypes.py b/mpet/daeVariableTypes.py index be994d10..e7a253b8 100644 --- a/mpet/daeVariableTypes.py +++ b/mpet/daeVariableTypes.py @@ -10,3 +10,6 @@ elec_pot_t = dae.daeVariableType( name="elec_pot_t", units=dae.unit(), lowerBound=-1e20, upperBound=1e20, initialGuess=0, absTolerance=1.e-6) +temp_t = dae.daeVariableType( + name="temp_t", units=dae.unit(), lowerBound=0.01, + upperBound=1e20, initialGuess=1, absTolerance=1.e-6) diff --git a/mpet/electrode/materials/LiC6.py b/mpet/electrode/materials/LiC6.py index eca5eaf1..f905f01c 100644 --- a/mpet/electrode/materials/LiC6.py +++ b/mpet/electrode/materials/LiC6.py @@ -1,17 +1,17 @@ import numpy as np -def LiC6(self, y, ybar, muR_ref): +def LiC6(self, y, ybar, T, muR_ref): """ Ferguson and Bazant 2014 """ muRtheta = -self.eokT*0.12 muR1homog, muR2homog = self.graphite_2param_homog( - y, self.get_trode_param("Omega_a"), self.get_trode_param("Omega_b"), + y, T, self.get_trode_param("Omega_a"), self.get_trode_param("Omega_b"), self.get_trode_param("Omega_c"), self.get_trode_param("EvdW")) muR1nonHomog, muR2nonHomog = self.general_non_homog(y, ybar) muR1 = muR1homog + muR1nonHomog muR2 = muR2homog + muR2nonHomog - actR1 = np.exp(muR1/self.T) - actR2 = np.exp(muR2/self.T) + actR1 = np.exp(muR1/T) + actR2 = np.exp(muR2/T) muR1 += muRtheta + muR_ref muR2 += muRtheta + muR_ref return (muR1, muR2), (actR1, actR2) diff --git a/mpet/electrode/materials/LiC6_1param.py b/mpet/electrode/materials/LiC6_1param.py index 4aa74e5d..a8bc95c1 100644 --- a/mpet/electrode/materials/LiC6_1param.py +++ b/mpet/electrode/materials/LiC6_1param.py @@ -1,12 +1,12 @@ import numpy as np -def LiC6_1param(self, y, ybar, muR_ref): +def LiC6_1param(self, y, ybar, T, muR_ref, ISfuncs=None): muRtheta = -self.eokT*0.12 muRhomog = self.graphite_1param_homog_3( - y, self.get_trode_param("Omega_a"), self.get_trode_param("Omega_b")) + y, T, self.get_trode_param("Omega_a"), self.get_trode_param("Omega_b")) muRnonHomog = self.general_non_homog(y, ybar) muR = muRhomog + muRnonHomog - actR = np.exp(muR/self.T) + actR = np.exp(muR/T) muR += muRtheta + muR_ref return muR, actR diff --git a/mpet/electrode/materials/LiC6_2step_ss.py b/mpet/electrode/materials/LiC6_2step_ss.py index f8aa83fd..7b3b13e6 100644 --- a/mpet/electrode/materials/LiC6_2step_ss.py +++ b/mpet/electrode/materials/LiC6_2step_ss.py @@ -2,7 +2,7 @@ from mpet.props_am import step_up, step_down -def LiC6_2step_ss(self, y, ybar, muR_ref): +def LiC6_2step_ss(self, y, ybar, T, muR_ref): """ Fit function to the OCV predicted by the phase separating 2-variable graphite model (LiC6 function in this class). @@ -14,8 +14,8 @@ def LiC6_2step_ss(self, y, ybar, muR_ref): rEdge = 1 - edgeLen width = 1e-4 vshift = 1e-2 - lSide = -((np.log(y/(1-y)) - np.log(lEdge/(1-lEdge)) - vshift)*step_down(y, lEdge, width)) - rSide = -((np.log(y/(1-y)) - np.log(rEdge/(1-rEdge)) + vshift)*step_up(y, rEdge, width)) + lSide = -((np.log(y/(1-y))*T - np.log(lEdge/(1-lEdge)) - vshift)*step_down(y, lEdge, width)) + rSide = -((np.log(y/(1-y))*T - np.log(rEdge/(1-rEdge)) + vshift)*step_up(y, rEdge, width)) OCV = ( Vstd + Vstep*(step_down(y, 0.5, 0.013) - 1) diff --git a/mpet/electrode/materials/LiC6_LIONSIMBA.py b/mpet/electrode/materials/LiC6_LIONSIMBA.py index 564dd58a..af1d01a5 100644 --- a/mpet/electrode/materials/LiC6_LIONSIMBA.py +++ b/mpet/electrode/materials/LiC6_LIONSIMBA.py @@ -1,9 +1,8 @@ import numpy as np -def LiC6_LIONSIMBA(self, y, ybar, muR_ref): +def LiC6_LIONSIMBA(self, y, ybar, T, muR_ref): """ Torchio et al, 2016. """ - T = self.T Tref = 298 r1 = 0.7222 r2 = 0.1387 diff --git a/mpet/electrode/materials/LiC6_coke_ss.py b/mpet/electrode/materials/LiC6_coke_ss.py index edc7efa3..326685fd 100644 --- a/mpet/electrode/materials/LiC6_coke_ss.py +++ b/mpet/electrode/materials/LiC6_coke_ss.py @@ -1,7 +1,7 @@ import numpy as np -def LiC6_coke_ss(self, y, ybar, muR_ref): +def LiC6_coke_ss(self, y, ybar, T, muR_ref): """ Doyle, Newman, 1996 """ OCV = (-0.16 + 1.32*np.exp(-3.0*y) + 10.*np.exp(-2000.*y)) muR = self.get_muR_from_OCV(OCV, muR_ref) diff --git a/mpet/electrode/materials/LiC6_coke_ss2.py b/mpet/electrode/materials/LiC6_coke_ss2.py index e0599cf4..5535bd21 100644 --- a/mpet/electrode/materials/LiC6_coke_ss2.py +++ b/mpet/electrode/materials/LiC6_coke_ss2.py @@ -1,7 +1,7 @@ import numpy as np -def LiC6_coke_ss2(self, y, ybar, muR_ref): +def LiC6_coke_ss2(self, y, ybar, T, muR_ref): """ Fuller, Doyle, Newman, 1994 """ c1 = -0.132056 c2 = 1.40854 diff --git a/mpet/electrode/materials/LiC6_ss.py b/mpet/electrode/materials/LiC6_ss.py index 3b017f5d..bf96712a 100644 --- a/mpet/electrode/materials/LiC6_ss.py +++ b/mpet/electrode/materials/LiC6_ss.py @@ -1,7 +1,7 @@ import numpy as np -def LiC6_ss(self, y, ybar, muR_ref): +def LiC6_ss(self, y, ybar, T, muR_ref): """ Safari, Delacourt 2011 """ OCV = (0.6379 + 0.5416*np.exp(-305.5309*y) + 0.044*np.tanh(-(y - 0.1958)/0.1088) diff --git a/mpet/electrode/materials/LiC6_ss2.py b/mpet/electrode/materials/LiC6_ss2.py index 01a1f2de..ec6d62a8 100644 --- a/mpet/electrode/materials/LiC6_ss2.py +++ b/mpet/electrode/materials/LiC6_ss2.py @@ -1,7 +1,7 @@ from mpet.props_am import step_down -def LiC6_ss2(self, y, ybar, muR_ref): +def LiC6_ss2(self, y, ybar, T, muR_ref): """ Bernardi and Go 2011 """ p1, p2, p3, p4 = (0.085, 0.120, 0.210, 3.5) sfac = 0.3 diff --git a/mpet/electrode/materials/LiCoO2_LIONSIMBA.py b/mpet/electrode/materials/LiCoO2_LIONSIMBA.py index bd248f15..03343eaa 100644 --- a/mpet/electrode/materials/LiCoO2_LIONSIMBA.py +++ b/mpet/electrode/materials/LiCoO2_LIONSIMBA.py @@ -1,6 +1,5 @@ -def LiCoO2_LIONSIMBA(self, y, ybar, muR_ref): +def LiCoO2_LIONSIMBA(self, y, ybar, T, muR_ref): """ Torchio et al, 2016. """ - T = self.T Tref = 298 r1 = 4.656 r2 = 88.669 diff --git a/mpet/electrode/materials/LiFePO4.py b/mpet/electrode/materials/LiFePO4.py index 071d3876..7e184a5d 100644 --- a/mpet/electrode/materials/LiFePO4.py +++ b/mpet/electrode/materials/LiFePO4.py @@ -1,12 +1,12 @@ import numpy as np -def LiFePO4(self, y, ybar, muR_ref): +def LiFePO4(self, y, ybar, T, muR_ref): """ Bai, Cogswell, Bazant 2011 """ muRtheta = -self.eokT*3.422 - muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a")) + muRhomog = self.reg_sln(y, T, self.get_trode_param("Omega_a")) muRnonHomog = self.general_non_homog(y, ybar) muR = muRhomog + muRnonHomog - actR = np.exp(muR/self.T) + actR = np.exp(muR/T) muR += muRtheta + muR_ref return muR, actR diff --git a/mpet/electrode/materials/LiMn2O4_ss.py b/mpet/electrode/materials/LiMn2O4_ss.py index 81cabff5..c759c375 100644 --- a/mpet/electrode/materials/LiMn2O4_ss.py +++ b/mpet/electrode/materials/LiMn2O4_ss.py @@ -1,7 +1,7 @@ import numpy as np -def LiMn2O4_ss(self, y, ybar, muR_ref): +def LiMn2O4_ss(self, y, ybar, T, muR_ref): """ Doyle, Newman, 1996 """ # OCV in V vs Li/Li+ OCV = (4.19829 + 0.0565661*np.tanh(-14.5546*y + 8.60942) diff --git a/mpet/electrode/materials/LiMn2O4_ss2.py b/mpet/electrode/materials/LiMn2O4_ss2.py index ca7615cb..d36a7412 100644 --- a/mpet/electrode/materials/LiMn2O4_ss2.py +++ b/mpet/electrode/materials/LiMn2O4_ss2.py @@ -1,7 +1,7 @@ import numpy as np -def LiMn2O4_ss2(self, y, ybar, muR_ref): +def LiMn2O4_ss2(self, y, ybar, T, muR_ref): """ Fuller, Doyle, Newman, 1994 """ # OCV in V vs Li/Li+ OCV = (4.06279 + 0.0677504*np.tanh(-21.8502*y + 12.8268) diff --git a/mpet/electrode/materials/Li_ss.py b/mpet/electrode/materials/Li_ss.py index 25a3464b..3d90b07e 100644 --- a/mpet/electrode/materials/Li_ss.py +++ b/mpet/electrode/materials/Li_ss.py @@ -1,4 +1,4 @@ -def Li_ss(self, y, ybar, muR_ref, ISfuncs=None): +def Li_ss(self, y, ybar, T, muR_ref, ISfuncs=None): muR = 0.*y + muR_ref actR = 0.*y + 1. return muR, actR diff --git a/mpet/electrode/materials/NCA_ss1.py b/mpet/electrode/materials/NCA_ss1.py index 545ed5cf..3665dea5 100644 --- a/mpet/electrode/materials/NCA_ss1.py +++ b/mpet/electrode/materials/NCA_ss1.py @@ -1,13 +1,13 @@ import numpy as np -def NCA_ss1(self, y, ybar, muR_ref): +def NCA_ss1(self, y, ybar, T, muR_ref): """ This function was obtained from Dan Cogswell's fit of Samsung data. """ OCV = (3.86 + 1.67*y - 9.52*y**2 + 15.04*y**3 - 7.95*y**4 - - 0.06*np.log(y/(1-y))) + - 0.06*T*np.log(y/(1-y))) muR = self.get_muR_from_OCV(OCV, muR_ref) actR = None return muR, actR diff --git a/mpet/electrode/materials/NCA_ss2.py b/mpet/electrode/materials/NCA_ss2.py index 6344e710..40dd328f 100644 --- a/mpet/electrode/materials/NCA_ss2.py +++ b/mpet/electrode/materials/NCA_ss2.py @@ -1,7 +1,7 @@ import numpy as np -def NCA_ss2(self, y, ybar, muR_ref): +def NCA_ss2(self, y, ybar, T, muR_ref): """ Li_q Ni(0.8)Co(0.15)Al(0.05)O2 as a function of y. Here, y actually represents a practical @@ -10,7 +10,7 @@ def NCA_ss2(self, y, ybar, muR_ref): This function was obtained from a fit by Raymond B. Smith of Samsung data of a LiC6-NCA cell discharged at C/100. """ - OCV = (-self.kToe*np.log(y/(1-y)) + OCV = (-self.kToe*np.log(y/(1-y))*T + 4.12178 - 0.2338*y - 1.24566*y**2 + 1.16769*y**3 - 0.20745*y**4) muR = self.get_muR_from_OCV(OCV, muR_ref) diff --git a/mpet/electrode/materials/testIS_ss.py b/mpet/electrode/materials/testIS_ss.py index 545c6a91..0fe037f6 100644 --- a/mpet/electrode/materials/testIS_ss.py +++ b/mpet/electrode/materials/testIS_ss.py @@ -1,9 +1,9 @@ import numpy as np -def testIS_ss(self, y, ybar, muR_ref): +def testIS_ss(self, y, ybar, T, muR_ref): """Ideal solution material for testing.""" - OCV = -self.kToe*np.log(y/(1-y)) + OCV = -self.kToe*np.log(y/(1-y))*T muR = self.get_muR_from_OCV(OCV, muR_ref) actR = None return muR, actR diff --git a/mpet/electrode/materials/testRS.py b/mpet/electrode/materials/testRS.py index 3b61acab..5b418956 100644 --- a/mpet/electrode/materials/testRS.py +++ b/mpet/electrode/materials/testRS.py @@ -1,9 +1,9 @@ import numpy as np -def testRS(self, y, ybar, muR_ref): +def testRS(self, y, ybar, T, muR_ref): muRtheta = 0. - muR = self.reg_sln(y, self.get_trode_param("Omega_a")) - actR = np.exp(muR/self.T) + muR = self.reg_sln(y, T, self.get_trode_param("Omega_a")) + actR = np.exp(muR/T) muR += muRtheta + muR_ref return muR, actR diff --git a/mpet/electrode/materials/testRS_ps.py b/mpet/electrode/materials/testRS_ps.py index d638fb99..a3d3c7d4 100644 --- a/mpet/electrode/materials/testRS_ps.py +++ b/mpet/electrode/materials/testRS_ps.py @@ -1,11 +1,11 @@ import numpy as np -def testRS_ps(self, y, ybar, muR_ref, ISfuncs=None): +def testRS_ps(self, y, ybar, T, muR_ref, ISfuncs=None): muRtheta = -self.eokT*2. - muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a"), ISfuncs) + muRhomog = self.reg_sln(y, T, self.get_trode_param("Omega_a"), ISfuncs) muRnonHomog = self.general_non_homog(y, ybar) muR = muRhomog + muRnonHomog - actR = np.exp(muR/self.T) + actR = np.exp(muR/T) muR += muRtheta + muR_ref return muR, actR diff --git a/mpet/electrode/materials/testRS_ss.py b/mpet/electrode/materials/testRS_ss.py index 1b250198..a8d4b40e 100644 --- a/mpet/electrode/materials/testRS_ss.py +++ b/mpet/electrode/materials/testRS_ss.py @@ -1,7 +1,7 @@ from mpet.props_am import step_down, step_up -def testRS_ss(self, y, ybar, muR_ref, ISfuncs=None): +def testRS_ss(self, y, ybar, T, muR_ref, ISfuncs=None): """ Regular solution material which phase separates at binodal points, for modeling as a solid solution. For testing. @@ -9,7 +9,7 @@ def testRS_ss(self, y, ybar, muR_ref, ISfuncs=None): # Based Omg = 3*k*T_ref yL = 0.07072018 yR = 0.92927982 - OCV_rs = -self.kToe*self.reg_sln(y, self.get_trode_param("Omega_a"), ISfuncs) + OCV_rs = -self.kToe*self.reg_sln(y, T, self.get_trode_param("Omega_a"), ISfuncs) width = 0.005 OCV = OCV_rs*step_down(y, yL, width) + OCV_rs*step_up(y, yR, width) + 2 muR = self.get_muR_from_OCV(OCV, muR_ref) diff --git a/mpet/geometry.py b/mpet/geometry.py index c35dbba1..e8dd2c79 100644 --- a/mpet/geometry.py +++ b/mpet/geometry.py @@ -71,14 +71,16 @@ def calc_curv(c, dr, r_vec, Rs, beta_s, particleShape): return curv -def get_elyte_disc(Nvol, L, poros, BruggExp): +def get_elyte_disc(Nvol, L, poros, BruggExp, k_h): out = {} # Width of each cell out["dxvec"] = utils.get_dxvec(L, Nvol) # Distance between cell centers dxtmp = np.hstack((out["dxvec"][0], out["dxvec"], out["dxvec"][-1])) + out["dx"] = dxtmp out["dxd1"] = utils.mean_linear(dxtmp) + out["dxd2"] = dxtmp[1:-1] # for thermal finite differences # The porosity vector out["porosvec"] = utils.get_asc_vec(poros, Nvol) @@ -87,7 +89,12 @@ def get_elyte_disc(Nvol, L, poros, BruggExp): # Vector of Bruggeman exponents Brugg_pad = utils.pad_vec(utils.get_asc_vec(BruggExp, Nvol)) + # The porosity vector + khvec = utils.get_asc_vec(k_h, Nvol) + out["khvec"] = utils.pad_vec(khvec) + # Vector of posority/tortuosity (assuming Bruggeman) out["eps_o_tau"] = porosvec_pad/porosvec_pad**(Brugg_pad) + out["min_eps_o_tau"] = (1-porosvec_pad)**(1-Brugg_pad) return out diff --git a/mpet/mod_cell.py b/mpet/mod_cell.py index 46d244d2..7e7a68ec 100644 --- a/mpet/mod_cell.py +++ b/mpet/mod_cell.py @@ -18,7 +18,7 @@ import mpet.ports as ports import mpet.utils as utils from mpet.config import constants -from mpet.daeVariableTypes import mole_frac_t, elec_pot_t, conc_t +from mpet.daeVariableTypes import mole_frac_t, elec_pot_t, conc_t, temp_t # Dictionary of end conditions endConditions = { @@ -58,7 +58,9 @@ def __init__(self, config, Name, Parent=None, Description=""): self.phi_bulk = {} self.phi_part = {} self.R_Vp = {} + self.Q_Vp = {} self.ffrac = {} + self.T_lyte = {} for trode in trodes: # Concentration/potential in electrode regions of elyte self.c_lyte[trode] = dae.daeVariable( @@ -81,9 +83,17 @@ def __init__(self, config, Name, Parent=None, Description=""): "R_Vp_{trode}".format(trode=trode), dae.no_t, self, "Rate of reaction of positives per electrode volume", [self.DmnCell[trode]]) + self.Q_Vp[trode] = dae.daeVariable( + "Q_Vp_{trode}".format(trode=trode), dae.no_t, self, + "Rate of heat generation of positives per electrode volume", + [self.DmnCell[trode]]) self.ffrac[trode] = dae.daeVariable( "ffrac_{trode}".format(trode=trode), mole_frac_t, self, "Overall filling fraction of solids in electrodes") + self.T_lyte[trode] = dae.daeVariable( + "T_lyte_{trode}".format(trode=trode), temp_t, self, + "Temperature in the elyte in electrode {trode}".format(trode=trode), + [self.DmnCell[trode]]) if config['have_separator']: # If we have a separator self.c_lyte["s"] = dae.daeVariable( "c_lyte_s", conc_t, self, @@ -93,6 +103,10 @@ def __init__(self, config, Name, Parent=None, Description=""): "phi_lyte_s", elec_pot_t, self, "Electrostatic potential in electrolyte in separator", [self.DmnCell["s"]]) + self.T_lyte["s"] = dae.daeVariable( + "T_lyte_s", temp_t, self, + "Temperature in electrolyte in separator", + [self.DmnCell["s"]]) # Note if we're doing a single electrode volume simulation # It will be in a perfect bath of electrolyte at the applied # potential. @@ -104,6 +118,10 @@ def __init__(self, config, Name, Parent=None, Description=""): self.c_lyteGP_L = dae.daeVariable("c_lyteGP_L", conc_t, self, "c_lyte left BC GP") self.phi_lyteGP_L = dae.daeVariable( "phi_lyteGP_L", elec_pot_t, self, "phi_lyte left BC GP") + self.T_lyteGP_L = dae.daeVariable( + "T_lyteGP_L", temp_t, self, "T_lyte left BC GP") + self.T_lyteGP_R = dae.daeVariable( + "T_lyteGP_R", temp_t, self, "T_lyte left BC GP") self.phi_applied = dae.daeVariable( "phi_applied", elec_pot_t, self, "Overall battery voltage (at anode current collector)") @@ -195,6 +213,23 @@ def DeclareEquations(self): * self.particles[trode][vInd,pInd].dcbardt()) eq.Residual = self.R_Vp[trode](vInd) - RHS + # Define dimensionless R_Vp for each electrode volume + for trode in trodes: + for vInd in range(Nvol[trode]): + eq = self.CreateEquation( + "Q_Vp_trode{trode}vol{vInd}".format(vInd=vInd, trode=trode)) + # Start with no reaction, then add reactions for each + # particle in the volume. + RHS = 0 + # sum over particle volumes in given electrode volume + for pInd in range(Npart[trode]): + # The volume of this particular particle + Vj = config["psd_vol_FracVol"][trode][vInd,pInd] + RHS += -(config["beta"][trode] * (1-config["poros"][trode]) + * config["P_L"][trode] * Vj + * self.particles[trode][vInd,pInd].q_rxn_bar()) + eq.Residual = self.Q_Vp[trode](vInd) - RHS + # Define output port variables for trode in trodes: for vInd in range(Nvol[trode]): @@ -202,6 +237,10 @@ def DeclareEquations(self): "portout_c_trode{trode}vol{vInd}".format(vInd=vInd, trode=trode)) eq.Residual = (self.c_lyte[trode](vInd) - self.portsOutLyte[trode][vInd].c_lyte()) + eq = self.CreateEquation( + "portout_T_trode{trode}vol{vInd}".format(vInd=vInd, trode=trode)) + eq.Residual = (self.T_lyte[trode](vInd) + - self.portsOutLyte[trode][vInd].T_lyte()) eq = self.CreateEquation( "portout_p_trode{trode}vol{vInd}".format(vInd=vInd, trode=trode)) phi_lyte = self.phi_lyte[trode](vInd) @@ -286,18 +325,28 @@ def DeclareEquations(self): eq.Residual = self.c_lyte["c"].dt(0) - 0 eq = self.CreateEquation("phi_lyte") eq.Residual = self.phi_lyte["c"](0) - self.phi_cell() + eq = self.CreateEquation("T_lyte") + eq.Residual = self.T_lyte["c"].dt(0) - 0 else: - disc = geom.get_elyte_disc(Nvol, config["L"], config["poros"], config["BruggExp"]) + disc = geom.get_elyte_disc(Nvol, config["L"], config["poros"], config["BruggExp"], + config["k_h"]) cvec = utils.get_asc_vec(self.c_lyte, Nvol) dcdtvec = utils.get_asc_vec(self.c_lyte, Nvol, dt=True) phivec = utils.get_asc_vec(self.phi_lyte, Nvol) + Tvec = utils.get_asc_vec(self.T_lyte, Nvol) + dTdtvec = utils.get_asc_vec(self.T_lyte, Nvol, dt=True) Rvvec = utils.get_asc_vec(self.R_Vp, Nvol) + Qvvec = utils.get_asc_vec(self.Q_Vp, Nvol) + rhocp_vec = utils.get_thermal_vec(Nvol, config) # Apply concentration and potential boundary conditions # Ghost points on the left and no-gradients on the right ctmp = np.hstack((self.c_lyteGP_L(), cvec, cvec[-1])) + # temperature uses a constant boundary condition + Ttmp = np.hstack((self.T_lyteGP_L(), Tvec, self.T_lyteGP_R())) phitmp = np.hstack((self.phi_lyteGP_L(), phivec, phivec[-1])) - Nm_edges, i_edges = get_lyte_internal_fluxes(ctmp, phitmp, disc, config) + Nm_edges, i_edges, q_edges = get_lyte_internal_fluxes(ctmp, phitmp, Ttmp, disc, + config, Nvol) # If we don't have a porous anode: # 1) the total current flowing into the electrolyte is set @@ -312,6 +361,7 @@ def DeclareEquations(self): # We assume BV kinetics with alpha = 0.5, # exchange current density, ecd = k0_foil * c_lyte**(0.5) cWall = .5*(ctmp[0] + ctmp[1]) + TWall = .5*(Ttmp[0] + Ttmp[1]) ecd = config["k0_foil"]*cWall**0.5 # note negative current because positive current is # oxidation here @@ -323,7 +373,7 @@ def DeclareEquations(self): # phiWall = -eta + phi_cell [- T*ln(c)] phiWall = -eta + self.phi_cell() if config["elyteModelType"] == "dilute": - phiWall -= config["T"]*np.log(cWall) + phiWall -= TWall*np.log(cWall) eqP.Residual = phiWall - .5*(phitmp[0] + phitmp[1]) # We have a porous anode -- no flux of charge or anions through current collector @@ -331,8 +381,16 @@ def DeclareEquations(self): eqC.Residual = ctmp[0] - ctmp[1] eqP.Residual = phitmp[0] - phitmp[1] + # boundary equation for temperature variables. per volume + eqTL = self.CreateEquation("GhostPointT_L") + eqTR = self.CreateEquation("GhostPointT_R") + eqTL.Residual = Ttmp[0] - config["T"] + eqTR.Residual = Ttmp[-1] - config["T"] + dvgNm = np.diff(Nm_edges)/disc["dxvec"] dvgi = np.diff(i_edges)/disc["dxvec"] + dvgq = np.diff(q_edges)/disc["dxvec"] + q_ohm = get_ohmic_heat(ctmp, Ttmp, self.phi_lyte, self.phi_bulk, disc, config, Nvol) for vInd in range(Nlyte): # Mass Conservation (done with the anion, although "c" is neutral salt conc) eq = self.CreateEquation("lyte_mass_cons_vol{vInd}".format(vInd=vInd)) @@ -340,6 +398,17 @@ def DeclareEquations(self): # Charge Conservation eq = self.CreateEquation("lyte_charge_cons_vol{vInd}".format(vInd=vInd)) eq.Residual = -dvgi[vInd] + config["zp"]*Rvvec[vInd] + # Energy Conservation + if config['nonisothermal']: + # if heat generation is turned on. per volume. + eq = self.CreateEquation("lyte_energy_cons_vol{vInd}".format(vInd=vInd)) + eq.Residual = rhocp_vec[vInd] * dTdtvec[vInd] - \ + q_ohm[vInd] - Qvvec[vInd] + dvgq[vInd] + else: + # if heat generation is turned off + eq = self.CreateEquation("lyte_energy_cons_vol{vInd}".format(vInd=vInd)) + eq.Residual = dTdtvec[vInd] - 0 + # add heat generation from reaction later # Define the total current. This must be done at the capacity # limiting electrode because currents are specified in @@ -473,26 +542,28 @@ def DeclareEquations(self): setVariableValues=[(self.endCondition, 2)]) -def get_lyte_internal_fluxes(c_lyte, phi_lyte, disc, config): +def get_lyte_internal_fluxes(c_lyte, phi_lyte, T_lyte, disc, config, Nvol): zp, zm, nup, num = config["zp"], config["zm"], config["nup"], config["num"] nu = nup + num - T = config["T"] dxd1 = disc["dxd1"] eps_o_tau = disc["eps_o_tau"] # Get concentration at cell edges using weighted mean wt = utils.pad_vec(disc["dxvec"]) c_edges_int = utils.weighted_linear_mean(c_lyte, wt) + T_edges_int = utils.weighted_linear_mean(T_lyte, wt) + k_h = utils.weighted_linear_mean(disc["khvec"], wt) if config["elyteModelType"] == "dilute": # Get porosity at cell edges using weighted harmonic mean eps_o_tau_edges = utils.weighted_linear_mean(eps_o_tau, wt) Dp = eps_o_tau_edges * config["Dp"] Dm = eps_o_tau_edges * config["Dm"] Nm_edges_int = num*(-Dm*np.diff(c_lyte)/dxd1 - - Dm/T*zm*c_edges_int*np.diff(phi_lyte)/dxd1) + - Dm/T_edges_int*zm*c_edges_int*np.diff(phi_lyte)/dxd1) i_edges_int = (-((nup*zp*Dp + num*zm*Dm)*np.diff(c_lyte)/dxd1) - - (nup*zp**2*Dp + num*zm**2*Dm)/T*c_edges_int*np.diff(phi_lyte)/dxd1) + - (nup*zp**2*Dp + num*zm**2*Dm)/T_edges_int + * c_edges_int*np.diff(phi_lyte)/dxd1) elif config["elyteModelType"] == "SM": SMset = config["SMset"] elyte_function = utils.import_function(config["SMset_filename"], SMset, @@ -500,18 +571,64 @@ def get_lyte_internal_fluxes(c_lyte, phi_lyte, disc, config): D_fs, sigma_fs, thermFac, tp0 = elyte_function()[:-1] # Get diffusivity and conductivity at cell edges using weighted harmonic mean - D_edges = utils.weighted_harmonic_mean(eps_o_tau*D_fs(c_lyte, T), wt) - sigma_edges = utils.weighted_harmonic_mean(eps_o_tau*sigma_fs(c_lyte, T), wt) + D_edges = utils.weighted_harmonic_mean(eps_o_tau*D_fs(c_lyte, T_lyte), wt) + sigma_edges = utils.weighted_harmonic_mean(eps_o_tau*sigma_fs(c_lyte, T_lyte), wt) sp, n = config["sp"], config["n"] # there is an error in the MPET paper, temperature dependence should be # in sigma and not outside of sigma i_edges_int = -sigma_edges * ( np.diff(phi_lyte)/dxd1 - + nu*T*(sp/(n*nup)+tp0(c_edges_int, T)/(zp*nup)) - * thermFac(c_edges_int, T) + + nu*T_edges_int*(sp/(n*nup)+tp0(c_edges_int, T_edges_int)/(zp*nup)) + * thermFac(c_edges_int, T_edges_int) * np.diff(np.log(c_lyte))/dxd1 ) Nm_edges_int = num*(-D_edges*np.diff(c_lyte)/dxd1 - + (1./(num*zm)*(1-tp0(c_edges_int, T))*i_edges_int)) - return Nm_edges_int, i_edges_int + + (1./(num*zm)*(1-tp0(c_edges_int, T_edges_int))*i_edges_int)) + q_edges_int = -k_h*np.diff(T_lyte)/dxd1 + # replace boundary conditions since they are ghost points with convective BCs + q_edges_int[0] = config["h_h"]*(1 - T_lyte[1]) + q_edges_int[-1] = config["h_h"]*(T_lyte[-2] - 1) + return Nm_edges_int, i_edges_int, q_edges_int + + +def get_ohmic_heat(c_lyte, T_lyte, phi_lyte, phi_bulk, disc, config, Nvol): + eps_o_tau = disc["eps_o_tau"] + min_eps_o_tau = disc["min_eps_o_tau"] + dx = disc["dxd2"] + + wt = utils.pad_vec(disc["dxvec"]) + sigma_s = utils.get_asc_vec(config["sigma_s"], Nvol) + c_edges_int = utils.weighted_linear_mean(c_lyte, wt) + dphilytedx = utils.central_diff_lyte(phi_lyte, Nvol, dx) + dphibulkdx = utils.central_diff_bulk(phi_bulk, Nvol, dx) + c_mid = c_lyte[1:-1] + T_mid = T_lyte[1:-1] + + # initialize sigma + q_ohmic = 0 + + if config["elyteModelType"] == "dilute": + zp, zm, nup, num = config["zp"], config["zm"], config["nup"], config["num"] + + # Get porosity at cell edges using weighted harmonic mean + Dp = eps_o_tau[1:-1] * config["Dp"] + Dm = eps_o_tau[1:-1] * config["Dm"] + sigma_l = ((nup*zp*Dp + num*zm*Dm)*np.diff(c_edges_int)/dx) - \ + (nup*zp ** 2*Dp + num*zm**2*Dm)/T_mid*c_mid + elif config["elyteModelType"] == "SM": + SMset = config["SMset"] + elyte_function = utils.import_function(config["SMset_filename"], SMset, + mpet_module=f"mpet.electrolyte.{SMset}") + sigma_fs, thermFac, tp0 = elyte_function()[1:-1] + # Get diffusivity and conductivity at cell edges using weighted harmonic mean + sigma_l = eps_o_tau[1:-1]*sigma_fs(c_mid, T_mid) + q_ohmic = q_ohmic + 2*sigma_l*(1-tp0(c_mid, T_mid))*T_mid \ + * np.diff(np.log(c_edges_int))/dx*dphilytedx + # this is going to be dra + sigma_s = min_eps_o_tau[1:-1] * sigma_s + + q_ohmic = q_ohmic + sigma_s*dphibulkdx**2 + \ + sigma_l*dphilytedx**2 + # do we have to extrapolate these + return q_ohmic diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index bbb241db..c9d2be4d 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -53,6 +53,10 @@ def __init__(self, config, trode, vInd, pInd, "c2bar", mole_frac_t, self, "Average concentration in 'layer' 2 of active particle") self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") + self.dcbar1dt = dae.daeVariable("dcbar1dt", dae.no_t, self, "Rate of particle 1 filling") + self.dcbar2dt = dae.daeVariable("dcbar2dt", dae.no_t, self, "Rate of particle 2 filling") + self.q_rxn_bar = dae.daeVariable( + "q_rxn_bar", dae.no_t, self, "Rate of heat generation in particle") if self.get_trode_param("type") not in ["ACR2"]: self.Rxn1 = dae.daeVariable("Rxn1", dae.no_t, self, "Rate of reaction 1") self.Rxn2 = dae.daeVariable("Rxn2", dae.no_t, self, "Rate of reaction 2") @@ -73,6 +77,7 @@ def __init__(self, config, trode, vInd, pInd, "portInBulk", dae.eInletPort, self, "Inlet port from e- conducting phase") self.phi_lyte = self.portInLyte.phi_lyte + self.T_lyte = self.portInLyte.T_lyte self.c_lyte = self.portInLyte.c_lyte self.phi_m = self.portInBulk.phi_m @@ -89,7 +94,6 @@ def get_trode_param(self, item): def DeclareEquations(self): dae.daeModel.DeclareEquations(self) N = self.get_trode_param("N") # number of grid points in particle - T = self.config["T"] # nondimensional temperature r_vec, volfrac_vec = geo.get_unit_solid_discr(self.get_trode_param('shape'), N) # Prepare noise @@ -108,7 +112,7 @@ def DeclareEquations(self): # Figure out mu_O, mu of the oxidized state mu_O, act_lyte = calc_mu_O( - self.c_lyte(), self.phi_lyte(), self.phi_m(), T, + self.c_lyte(), self.phi_lyte(), self.phi_m(), self.T_lyte(), self.config["elyteModelType"]) # Define average filling fractions in particle @@ -128,26 +132,50 @@ def DeclareEquations(self): for k in range(N): eq.Residual -= .5*(self.c1.dt(k) + self.c2.dt(k)) * volfrac_vec[k] + # Define average rate of filling of particle for cbar1 + eq = self.CreateEquation("dcbar1dt") + eq.Residual = self.dcbar1dt() + for k in range(N): + eq.Residual -= self.c1.dt(k) * volfrac_vec[k] + + # Define average rate of filling of particle for cbar1 + eq = self.CreateEquation("dcbar2dt") + eq.Residual = self.dcbar2dt() + for k in range(N): + eq.Residual -= self.c2.dt(k) * volfrac_vec[k] + c1 = np.empty(N, dtype=object) c2 = np.empty(N, dtype=object) c1[:] = [self.c1(k) for k in range(N)] c2[:] = [self.c2(k) for k in range(N)] if self.get_trode_param("type") in ["diffn2", "CHR2"]: # Equations for 1D particles of 1 field varible - self.sld_dynamics_1D2var(c1, c2, mu_O, act_lyte, noises) + eta1, eta2, c_surf1, c_surf2 = self.sld_dynamics_1D2var(c1, c2, mu_O, act_lyte, + noises) elif self.get_trode_param("type") in ["homog2", "homog2_sdn"]: # Equations for 0D particles of 1 field variables - self.sld_dynamics_0D2var(c1, c2, mu_O, act_lyte, noises) + eta1, eta2, c_surf1, c_surf2 = self.sld_dynamics_0D2var(c1, c2, mu_O, act_lyte, + noises) + + # Define average rate of heat generation + eq = self.CreateEquation("q_rxn_bar") + if self.config["entropy_heat_gen"]: + eq.Residual = self.q_rxn_bar() - 0.5 * self.dcbar1dt() * \ + (eta1 - self.T_lyte()*(np.log(c_surf1/(1-c_surf1))-1/self.c_lyte())) \ + - 0.5 * self.dcbar2dt() * (eta2 - self.T_lyte() + * (np.log(c_surf2/(1-c_surf2))-1/self.c_lyte())) + else: + eq.Residual = self.q_rxn_bar() - 0.5 * self.dcbar1dt() * eta1 \ + - 0.5 * self.dcbar2dt() * eta2 for eq in self.Equations: eq.CheckUnitsConsistency = False def sld_dynamics_0D2var(self, c1, c2, muO, act_lyte, noises): - T = self.config["T"] c1_surf = c1 c2_surf = c2 (mu1R_surf, mu2R_surf), (act1R_surf, act2R_surf) = calc_muR( - (c1_surf, c2_surf), (self.c1bar(), self.c2bar()), self.config, + (c1_surf, c2_surf), (self.c1bar(), self.c2bar()), self.T_lyte(), self.config, self.trode, self.ind) eta1 = calc_eta(mu1R_surf, muO) eta2 = calc_eta(mu2R_surf, muO) @@ -159,11 +187,11 @@ def sld_dynamics_0D2var(self, c1, c2, muO, act_lyte, noises): eta2_eff += noise2(dae.Time().Value) Rxn1 = self.calc_rxn_rate( eta1_eff, c1_surf, self.c_lyte(), self.get_trode_param("k0"), - self.get_trode_param("E_A"), T, act1R_surf, act_lyte, + self.get_trode_param("E_A"), self.T_lyte(), act1R_surf, act_lyte, self.get_trode_param("lambda"), self.get_trode_param("alpha")) Rxn2 = self.calc_rxn_rate( eta2_eff, c2_surf, self.c_lyte(), self.get_trode_param("k0"), - self.get_trode_param("E_A"), T, act2R_surf, act_lyte, + self.get_trode_param("E_A"), self.T_lyte(), act2R_surf, act_lyte, self.get_trode_param("lambda"), self.get_trode_param("alpha")) eq1 = self.CreateEquation("Rxn1") eq2 = self.CreateEquation("Rxn2") @@ -174,10 +202,10 @@ def sld_dynamics_0D2var(self, c1, c2, muO, act_lyte, noises): eq2 = self.CreateEquation("dc2sdt") eq1.Residual = self.c1.dt(0) - self.get_trode_param("delta_L")*Rxn1[0] eq2.Residual = self.c2.dt(0) - self.get_trode_param("delta_L")*Rxn2[0] + return eta1[-1], eta2[-1], c1_surf[-1], c2_surf[-1] def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises): N = self.get_trode_param("N") - T = self.config["T"] # Equations for concentration evolution # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors Mmat = get_Mmat(self.get_trode_param('shape'), N) @@ -185,8 +213,9 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises): # Get solid particle chemical potential, overpotential, reaction rate if self.get_trode_param("type") in ["diffn2", "CHR2"]: - (mu1R, mu2R), (act1R, act2R) = calc_muR( - (c1, c2), (self.c1bar(), self.c2bar()), self.config, self.trode, self.ind) + (mu1R, mu2R), (act1R, act2R) = calc_muR((c1, c2), (self.c1bar(), self.c2bar()), + self.T_lyte(), self.config, self.trode, + self.ind) c1_surf = c1[-1] c2_surf = c2[-1] mu1R_surf, act1R_surf = mu1R[-1], act1R[-1] @@ -203,11 +232,11 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises): eta2_eff = eta2 + self.Rxn2()*self.get_trode_param("Rfilm") Rxn1 = self.calc_rxn_rate( eta1_eff, c1_surf, self.c_lyte(), self.get_trode_param("k0"), - self.get_trode_param("E_A"), T, act1R_surf, act_lyte, + self.get_trode_param("E_A"), self.T_lyte(), act1R_surf, act_lyte, self.get_trode_param("lambda"), self.get_trode_param("alpha")) Rxn2 = self.calc_rxn_rate( eta2_eff, c2_surf, self.c_lyte(), self.get_trode_param("k0"), - self.get_trode_param("E_A"), T, act2R_surf, act_lyte, + self.get_trode_param("E_A"), self.T_lyte(), act2R_surf, act_lyte, self.get_trode_param("lambda"), self.get_trode_param("alpha")) if self.get_trode_param("type") in ["ACR2"]: for i in range(N): @@ -235,7 +264,8 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises): noise1, noise2 = noises Flux1_vec, Flux2_vec = calc_flux_CHR2( c1, c2, mu1R, mu2R, self.get_trode_param("D"), Dfunc, - self.get_trode_param("E_D"), Flux1_bc, Flux2_bc, dr, T, noise1, noise2) + self.get_trode_param("E_D"), Flux1_bc, Flux2_bc, dr, self.T_lyte(), + noise1, noise2) if self.get_trode_param("shape") == "sphere": area_vec = 4*np.pi*edges**2 elif self.get_trode_param("shape") == "cylinder": @@ -263,6 +293,11 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises): eq1.Residual = LHS1_vec[k] - RHS1[k] eq2.Residual = LHS2_vec[k] - RHS2[k] + if self.get_trode_param("type") in ["ACR"]: + return eta1[-1], eta2[-1], c1_surf[-1], c2_surf[-1] + else: + return eta1, eta2, c1_surf, c2_surf + class Mod1var(dae.daeModel): def __init__(self, config, trode, vInd, pInd, @@ -285,6 +320,8 @@ def __init__(self, config, trode, vInd, pInd, "cbar", mole_frac_t, self, "Average concentration in active particle") self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") + self.q_rxn_bar = dae.daeVariable( + "q_rxn_bar", dae.no_t, self, "Rate of heat generation in particle") if config[trode, "type"] not in ["ACR"]: self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction") else: @@ -304,6 +341,7 @@ def __init__(self, config, trode, vInd, pInd, "portInBulk", dae.eInletPort, self, "Inlet port from e- conducting phase") self.phi_lyte = self.portInLyte.phi_lyte + self.T_lyte = self.portInLyte.T_lyte self.c_lyte = self.portInLyte.c_lyte self.phi_m = self.portInBulk.phi_m @@ -320,7 +358,6 @@ def get_trode_param(self, item): def DeclareEquations(self): dae.daeModel.DeclareEquations(self) N = self.get_trode_param("N") # number of grid points in particle - T = self.config["T"] # nondimensional temperature r_vec, volfrac_vec = geo.get_unit_solid_discr(self.get_trode_param('shape'), N) # Prepare noise @@ -334,7 +371,7 @@ def DeclareEquations(self): bounds_error=False, fill_value=0.) # Figure out mu_O, mu of the oxidized state - mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T, + mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), self.T_lyte(), self.config["elyteModelType"]) # Define average filling fraction in particle @@ -353,18 +390,25 @@ def DeclareEquations(self): c[:] = [self.c(k) for k in range(N)] if self.get_trode_param("type") in ["ACR", "diffn", "CHR"]: # Equations for 1D particles of 1 field varible - self.sld_dynamics_1D1var(c, mu_O, act_lyte, self.noise) + eta, c_surf = self.sld_dynamics_1D1var(c, mu_O, act_lyte, self.noise) elif self.get_trode_param("type") in ["homog", "homog_sdn"]: # Equations for 0D particles of 1 field variables - self.sld_dynamics_0D1var(c, mu_O, act_lyte, self.noise) + eta, c_surf = self.sld_dynamics_0D1var(c, mu_O, act_lyte, self.noise) + + # Define average rate of heat generation + eq = self.CreateEquation("q_rxn_bar") + if self.config["entropy_heat_gen"]: + eq.Residual = self.q_rxn_bar() - self.dcbardt() * \ + (eta - self.T_lyte()*(np.log(c_surf/(1-c_surf))-1/self.c_lyte())) + else: + eq.Residual = self.q_rxn_bar() - self.dcbardt() * eta for eq in self.Equations: eq.CheckUnitsConsistency = False def sld_dynamics_0D1var(self, c, muO, act_lyte, noise): - T = self.config["T"] c_surf = c - muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), self.config, + muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), self.T_lyte(),self.config, self.trode, self.ind) eta = calc_eta(muR_surf, muO) eta_eff = eta + self.Rxn()*self.get_trode_param("Rfilm") @@ -372,17 +416,17 @@ def sld_dynamics_0D1var(self, c, muO, act_lyte, noise): eta_eff += noise[0]() Rxn = self.calc_rxn_rate( eta_eff, c_surf, self.c_lyte(), self.get_trode_param("k0"), - self.get_trode_param("E_A"), T, actR_surf, act_lyte, + self.get_trode_param("E_A"), self.T_lyte(), actR_surf, act_lyte, self.get_trode_param("lambda"), self.get_trode_param("alpha")) eq = self.CreateEquation("Rxn") eq.Residual = self.Rxn() - Rxn[0] eq = self.CreateEquation("dcsdt") eq.Residual = self.c.dt(0) - self.get_trode_param("delta_L")*self.Rxn() + return eta[-1], c_surf[-1] def sld_dynamics_1D1var(self, c, muO, act_lyte, noise): N = self.get_trode_param("N") - T = self.config["T"] # Equations for concentration evolution # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors Mmat = get_Mmat(self.get_trode_param('shape'), N) @@ -392,9 +436,10 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise): if self.get_trode_param("type") in ["ACR"]: c_surf = c muR_surf, actR_surf = calc_muR( - c_surf, self.cbar(), self.config, self.trode, self.ind) + c_surf, self.cbar(), self.T_lyte(), self.config, self.trode, self.ind) elif self.get_trode_param("type") in ["diffn", "CHR"]: - muR, actR = calc_muR(c, self.cbar(), self.config, self.trode, self.ind) + muR, actR = calc_muR(c, self.cbar(), self.T_lyte(), + self.config, self.trode, self.ind) c_surf = c[-1] muR_surf = muR[-1] if actR is None: @@ -409,7 +454,7 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise): eta_eff = eta + self.Rxn()*self.get_trode_param("Rfilm") Rxn = self.calc_rxn_rate( eta_eff, c_surf, self.c_lyte(), self.get_trode_param("k0"), - self.get_trode_param("E_A"), T, actR_surf, act_lyte, + self.get_trode_param("E_A"), self.T_lyte(), actR_surf, act_lyte, self.get_trode_param("lambda"), self.get_trode_param("alpha")) if self.get_trode_param("type") in ["ACR"]: for i in range(N): @@ -432,10 +477,12 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise): f"mpet.electrode.diffusion.{Dfunc_name}") if self.get_trode_param("type") == "diffn": Flux_vec = calc_flux_diffn(c, self.get_trode_param("D"), Dfunc, - self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + self.get_trode_param("E_D"), Flux_bc, dr, + self.T_lyte(), noise) elif self.get_trode_param("type") == "CHR": Flux_vec = calc_flux_CHR(c, muR, self.get_trode_param("D"), Dfunc, - self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + self.get_trode_param("E_D"), Flux_bc, dr, + self.T_lyte(), noise) if self.get_trode_param("shape") == "sphere": area_vec = 4*np.pi*edges**2 elif self.get_trode_param("shape") == "cylinder": @@ -449,6 +496,11 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise): eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) eq.Residual = LHS_vec[k] - RHS[k] + if self.get_trode_param("type") in ["ACR"]: + return eta[-1], c_surf[-1] + else: + return eta, c_surf + def calc_eta(muR, muO): return muR - muO @@ -482,9 +534,9 @@ def calc_flux_diffn(c, D, Dfunc, E_D, Flux_bc, dr, T, noise): Flux_vec[-1] = Flux_bc c_edges = utils.mean_linear(c) if noise is None: - Flux_vec[1:N] = -D * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * np.diff(c)/dr + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * np.diff(c)/dr else: - Flux_vec[1:N] = -D * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ np.diff(c + noise(dae.Time().Value))/dr return Flux_vec @@ -535,10 +587,10 @@ def calc_mu_O(c_lyte, phi_lyte, phi_sld, T, elyteModelType): return mu_O, act_lyte -def calc_muR(c, cbar, config, trode, ind): +def calc_muR(c, cbar, T, config, trode, ind): muRfunc = props_am.muRfuncs(config, trode, ind).muRfunc muR_ref = config[trode, "muR_ref"] - muR, actR = muRfunc(c, cbar, muR_ref) + muR, actR = muRfunc(c, cbar, T, muR_ref) return muR, actR diff --git a/mpet/plot/outmat2txt.py b/mpet/plot/outmat2txt.py index d802466b..579d52e8 100644 --- a/mpet/plot/outmat2txt.py +++ b/mpet/plot/outmat2txt.py @@ -47,6 +47,7 @@ elyteiHdr = ("Electrolyte Current Density [A/m^2]\n" + RowsStr + FCStr) elytediviHdr = ("Electrolyte Divergence of Current Density [A/m^3]\n" + RowsStr + CCStr) +tempHdr = ("Electrolyte Temperature [K]\n" + RowsStr + CCStr) seeDiscStr = "See discData.txt for particle indexing information." partStr = "partTrode{l}vol{j}part{i}_" @@ -76,7 +77,7 @@ def main(indir, genData=True, discData=True, elyteData=True, - csldData=True, cbarData=True, bulkpData=True): + csldData=True, cbarData=True, bulkpData=True, tempData=True): config = plot_data.show_data( indir, plot_type="params", print_flag=False, save_flag=False, data_only=True) @@ -261,4 +262,11 @@ def get_trode_str(tr): np.savetxt(os.path.join(indir, fname), bulkp_cData, delimiter=dlm, header=bulkpHdr) + if tempData: + tempMat = plot_data.show_data( + indir, plot_type="temp", print_flag=False, + save_flag=False, data_only=True)[1] + np.savetxt(os.path.join(indir, "elyteTempData.txt"), + tempMat, delimiter=dlm, header=tempHdr) + return diff --git a/mpet/plot/plot_data.py b/mpet/plot/plot_data.py index c9550bc4..65cc3b77 100644 --- a/mpet/plot/plot_data.py +++ b/mpet/plot/plot_data.py @@ -300,6 +300,27 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu fig.savefig("mpet_current.png", bbox_inches="tight") return fig, ax + # Plot maximum temperature profile + if plot_type == "max_temp": + T_sep, T_anode, T_cath = pfx + 'T_lyte_s', pfx + 'T_lyte_a', pfx + 'T_lyte_c' + datay_T = utils.get_dict_key(data, T_cath, squeeze=False) + if config["have_separator"]: + datay_s_T = utils.get_dict_key(data, T_sep, squeeze=False) + datay_T = np.hstack((datay_s_T, datay_T)) + if "a" in trodes: + datay_a_T = utils.get_dict_key(data, T_anode, squeeze=False) + datay_T = np.hstack((datay_a_T, datay_T)) + Tmax = np.max(datay_T * Tref, axis=1) + if data_only: + return times*td, Tmax + fig, ax = plt.subplots(figsize=figsize) + ax.plot(times*td, Tmax) + ax.set_xlabel("Time [s]") + ax.set_ylabel("Maximum Temperature [K]") + if save_flag: + fig.savefig("mpet_max_temp.png", bbox_inches="tight") + return fig, ax + if plot_type == "power": current = utils.get_dict_key(data, pfx + 'current') * (3600/td) * (cap/3600) # in A/m^2 voltage = (Vstd - (k*Tref/e)*utils.get_dict_key(data, pfx + 'phi_applied')) # in V @@ -315,23 +336,26 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu return fig, ax # Plot electrolyte concentration or potential - elif plot_type in ["elytec", "elytep", "elytecf", "elytepf", - "elytei", "elyteif", "elytedivi", "elytedivif"]: + elif plot_type in ["elytec", "elytep", "elytecf", "elytepf", "elytei", + "elyteif", "elytedivi", "elytedivif", "temp"]: fplot = (True if plot_type[-1] == "f" else False) t0ind = (0 if not fplot else -1) datax = cellsvec - c_sep, p_sep = pfx + 'c_lyte_s', pfx + 'phi_lyte_s' - c_anode, p_anode = pfx + 'c_lyte_a', pfx + 'phi_lyte_a' - c_cath, p_cath = pfx + 'c_lyte_c', pfx + 'phi_lyte_c' + c_sep, p_sep, T_sep = pfx + 'c_lyte_s', pfx + 'phi_lyte_s', pfx + 'T_lyte_s' + c_anode, p_anode, T_anode = pfx + 'c_lyte_a', pfx + 'phi_lyte_a', pfx + 'T_lyte_a' + c_cath, p_cath, T_cath = pfx + 'c_lyte_c', pfx + 'phi_lyte_c', pfx + 'T_lyte_c' datay_c = utils.get_dict_key(data, c_cath, squeeze=False) datay_p = utils.get_dict_key(data, p_cath, squeeze=False) + datay_T = utils.get_dict_key(data, T_cath, squeeze=False) L_c = config['L']["c"] * config['L_ref'] * Lfac Ltot = L_c if config["have_separator"]: datay_s_c = utils.get_dict_key(data, c_sep, squeeze=False) datay_s_p = utils.get_dict_key(data, p_sep, squeeze=False) + datay_s_T = utils.get_dict_key(data, T_sep, squeeze=False) datay_c = np.hstack((datay_s_c, datay_c)) datay_p = np.hstack((datay_s_p, datay_p)) + datay_T = np.hstack((datay_s_T, datay_T)) L_s = config['L']["s"] * config['L_ref'] * Lfac Ltot += L_s else: @@ -339,8 +363,10 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu if "a" in trodes: datay_a_c = utils.get_dict_key(data, c_anode, squeeze=False) datay_a_p = utils.get_dict_key(data, p_anode, squeeze=False) + datay_a_T = utils.get_dict_key(data, T_anode, squeeze=False) datay_c = np.hstack((datay_a_c, datay_c)) datay_p = np.hstack((datay_a_p, datay_p)) + datay_T = np.hstack((datay_a_T, datay_T)) L_a = config['L']["a"] * config['L_ref'] * Lfac Ltot += L_a else: @@ -353,17 +379,23 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu elif plot_type in ["elytep", "elytepf"]: ylbl = 'Potential of electrolyte [V]' datay = datay_p*(k*Tref/e) - Vstd + elif plot_type in ["temp"]: + ylbl = 'Temperature [K]' + datay = datay_T * constants.T_ref elif plot_type in ["elytei", "elyteif", "elytedivi", "elytedivif"]: cGP_L = utils.get_dict_key(data, "c_lyteGP_L") pGP_L = utils.get_dict_key(data, "phi_lyteGP_L") + TGP_L = utils.get_dict_key(data, "T_lyteGP_L") cmat = np.hstack((cGP_L.reshape((-1,1)), datay_c, datay_c[:,-1].reshape((-1,1)))) pmat = np.hstack((pGP_L.reshape((-1,1)), datay_p, datay_p[:,-1].reshape((-1,1)))) + Tmat = np.hstack((TGP_L.reshape((-1,1)), datay_T, datay_T[:,-1].reshape((-1,1)))) disc = geom.get_elyte_disc( - Nvol, config["L"], config["poros"], config["BruggExp"]) + Nvol, config["L"], config["poros"], config["BruggExp"], config["k_h"]) i_edges = np.zeros((numtimes, len(facesvec))) for tInd in range(numtimes): + # no heat flux at boundary i_edges[tInd, :] = mod_cell.get_lyte_internal_fluxes( - cmat[tInd, :], pmat[tInd, :], disc, config)[1] + cmat[tInd, :], pmat[tInd, :], Tmat[tInd, :], disc, config, Nvol)[1] if plot_type in ["elytei", "elyteif"]: ylbl = r'Current density of electrolyte [A/m$^2$]' datax = facesvec diff --git a/mpet/ports.py b/mpet/ports.py index 6689fd1b..0c1686c6 100644 --- a/mpet/ports.py +++ b/mpet/ports.py @@ -10,6 +10,9 @@ def __init__(self, Name, PortType, Model, Description=""): self.c_lyte = dae.daeVariable( "c_lyte", mole_frac_t, self, "Concentration in the electrolyte") + self.T_lyte = dae.daeVariable( + "T_lyte", mole_frac_t, self, + "Temperature in the electrolyte") self.phi_lyte = dae.daeVariable( "phi_lyte", elec_pot_t, self, "Electric potential in the electrolyte") diff --git a/mpet/props_am.py b/mpet/props_am.py index 41dd6398..f7b1cb38 100644 --- a/mpet/props_am.py +++ b/mpet/props_am.py @@ -23,6 +23,7 @@ class muRfuncs(): muR -- chemical potential actR -- activity (if applicable, else None) """ + def __init__(self, config, trode, ind=None): """config is the full dictionary of parameters for the electrode particles, as made for the @@ -32,7 +33,6 @@ def __init__(self, config, trode, ind=None): self.config = config self.trode = trode self.ind = ind - self.T = config['T'] # nondimensional # eokT and kToe are the reference values for scalings self.eokT = constants.e / (constants.k * constants.T_ref) self.kToe = 1. / self.eokT @@ -69,32 +69,31 @@ def get_muR_from_OCV(self, OCV, muR_ref): # Helper functions ###### - def ideal_sln(self, y): + def ideal_sln(self, y, T): """ Helper function: Should not be called directly from simulation. Call a specific material instead. """ - T = self.T muR = T*np.log(y/(1-y)) return muR - def reg_sln(self, y, Omga): + def reg_sln(self, y, T, Omga): """ Helper function """ - muR_IS = self.ideal_sln(y) + muR_IS = self.ideal_sln(y, T) enthalpyTerm = Omga*(1-2*y) muR = muR_IS + enthalpyTerm return muR - def graphite_2param_homog(self, y, Omga, Omgb, Omgc, EvdW): + def graphite_2param_homog(self, y, T, Omga, Omgb, Omgc, EvdW): """ Helper function """ y1, y2 = y - muR1 = self.reg_sln(y1, Omga) - muR2 = self.reg_sln(y2, Omga) + muR1 = self.reg_sln(y1, T, Omga) + muR2 = self.reg_sln(y2, T, Omga) muR1 += Omgb*y2 + Omgc*y2*(1-y2)*(1-2*y1) muR2 += Omgb*y1 + Omgc*y1*(1-y1)*(1-2*y2) muR1 += EvdW * (30 * y1**2 * (1-y1)**2) muR2 += EvdW * (30 * y2**2 * (1-y2)**2) return (muR1, muR2) - def graphite_1param_homog(self, y, Omga, Omgb): + def graphite_1param_homog(self, y, T, Omga, Omgb): """ Helper function """ width = 5e-2 tailScl = 5e-2 @@ -106,7 +105,7 @@ def graphite_1param_homog(self, y, Omga, Omgb): muR = muLtail + muRtail + muLlin + muRlin return muR - def graphite_1param_homog_2(self, y, Omga, Omgb): + def graphite_1param_homog_2(self, y, T, Omga, Omgb): """ Helper function """ width = 5e-2 tailScl = 5e-2 @@ -124,7 +123,7 @@ def graphite_1param_homog_2(self, y, Omga, Omgb): muR = muLMod + muLtail + muRtail + muLlin + muRlin return muR - def graphite_1param_homog_3(self, y, Omga, Omgb): + def graphite_1param_homog_3(self, y, T, Omga, Omgb): """ Helper function with low hysteresis and soft tail """ width = 5e-2 tailScl = 5e-2 diff --git a/mpet/sim.py b/mpet/sim.py index c9130bc4..9fc31eed 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -70,6 +70,8 @@ def SetUpVariables(self): for i in range(Nvol[tr]): # Guess initial volumetric reaction rates self.m.R_Vp[tr].SetInitialGuess(i, 0.0) + # set initial temperature condition + self.m.T_lyte[tr].SetInitialCondition(i, config["T"]) # Guess initial value for the potential of the # electrodes if tr == "a": # anode @@ -116,22 +118,27 @@ def SetUpVariables(self): if not self.m.SVsim: self.m.c_lyteGP_L.SetInitialGuess(config["c0"]) self.m.phi_lyteGP_L.SetInitialGuess(0) + self.m.T_lyteGP_L.SetInitialGuess(config["T"]) + self.m.T_lyteGP_R.SetInitialGuess(config["T"]) # Separator electrolyte initialization if config["have_separator"]: for i in range(Nvol["s"]): self.m.c_lyte["s"].SetInitialCondition(i, config['c0']) + self.m.T_lyte["s"].SetInitialCondition(i, config['T']) self.m.phi_lyte["s"].SetInitialGuess(i, 0) # Anode and cathode electrolyte initialization for tr in config["trodes"]: for i in range(Nvol[tr]): self.m.c_lyte[tr].SetInitialCondition(i, config['c0']) + self.m.T_lyte[tr].SetInitialCondition(i, config['T']) self.m.phi_lyte[tr].SetInitialGuess(i, 0) # Set electrolyte concentration in each particle for j in range(Npart[tr]): self.m.particles[tr][i,j].c_lyte.SetInitialGuess(config["c0"]) + self.m.particles[tr][i,j].T_lyte.SetInitialGuess(config["T"]) else: dPrev = self.dataPrev @@ -153,6 +160,7 @@ def SetUpVariables(self): # Set the inlet port variables for each particle part.c_lyte.SetInitialGuess(data["c_lyte_" + tr][-1,i]) + part.T_lyte.SetInitialGuess(data["T_lyte_" + tr][-1,i]) part.phi_lyte.SetInitialGuess(data["phi_lyte_" + tr][-1,i]) part.phi_m.SetInitialGuess(data["phi_bulk_" + tr][-1,i]) @@ -178,12 +186,16 @@ def SetUpVariables(self): for i in range(Nvol["s"]): self.m.c_lyte["s"].SetInitialCondition( i, data["c_lyte_s"][-1,i]) + self.m.T_lyte["s"].SetInitialCondition( + i, data["T_lyte_s"][-1,i]) self.m.phi_lyte["s"].SetInitialGuess( i, data["phi_lyte_s"][-1,i]) for tr in config["trodes"]: for i in range(Nvol[tr]): self.m.c_lyte[tr].SetInitialCondition( i, data["c_lyte_" + tr][-1,i]) + self.m.T_lyte[tr].SetInitialCondition( + i, data["T_lyte_" + tr][-1,i]) self.m.phi_lyte[tr].SetInitialGuess( i, data["phi_lyte_" + tr][-1,i]) diff --git a/mpet/utils.py b/mpet/utils.py index d74a51b4..6843ddbc 100644 --- a/mpet/utils.py +++ b/mpet/utils.py @@ -28,6 +28,11 @@ def weighted_harmonic_mean(a, wt): return ((wt[1:]+wt[:-1])/(wt[1:]/a[1:]+wt[:-1]/a[:-1])) +def get_cell_Ntot(Nvol): + """Nvol is a dictionary containing the number of volumes in each simulated battery section.""" + return np.sum(list(Nvol.values())) + + def add_gp_to_vec(vec): """Add ghost points to the beginning and end of a vector for applying boundary conditions.""" out = np.empty(len(vec) + 2, dtype=object) @@ -82,6 +87,66 @@ def get_asc_vec(var, Nvol, dt=False): return out +def central_diff_bulk(array, Nvol, dx): + """Gets central diff for derivatives for use in thermal derivatives (which are split between + the individual electrodes) for bulk""" + varout = {} + for sectn in ["a", "c", "s"]: + # If we have information within this battery section + if sectn in ["a", "c"]: + if sectn in array.keys(): + # if it is one of the electrode sections + out = get_var_vec(array[sectn], Nvol[sectn]) + out = np.hstack((2*out[0]-out[1], out, 2*out[-1]-out[-2])) + varout[sectn] = (out[2:]-out[:-2]) + else: + varout[sectn] = np.zeros(0) + else: + # if anode does not exist + if sectn in Nvol: + varout[sectn] = np.zeros(Nvol[sectn]) + else: + varout[sectn] = np.zeros(0) + # sum solid + elyte poroisty + output = np.hstack((varout["a"], varout["s"], varout["c"]))/(2*dx) + return output + + +def central_diff_lyte(array, Nvol, dx): + """Gets central diff for derivatives for use in thermal derivatives (which are split between + the individual electrodes) for electrolyte""" + out = np.zeros(0) + for sectn in ["a", "s", "c"]: + # If we have information within this battery section + if sectn in array.keys(): + # if it is one of the electrode sections + out = np.append(out, get_var_vec(array[sectn], Nvol[sectn])) + else: + out = np.append(out, np.zeros(0)) + # now we have stacked everything + out = np.hstack((2*out[0]-out[1], out, 2*out[-1]-out[-2])) + output = (out[2:]-out[:-2])/(2*dx) + return output + + +def get_thermal_vec(Nvol, config): + """Get a numpy array for a variable spanning the anode, separator, and cathode.""" + varout = {} + for sectn in ["a", "s", "c"]: + # If we have information within this battery section + if sectn in Nvol: + # If it's an array of dae variable objects + out = config['rhom'][sectn] * config['cp'][sectn] + varout[sectn] = get_const_vec(out, Nvol[sectn]) + else: + # if anode does not exist + varout[sectn] = np.zeros(0) + + # sum solid + elyte poroisty + out = np.hstack((varout["a"], varout["s"], varout["c"])) + return out + + def get_dxvec(L, Nvol): """Get a vector of cell widths spanning the full cell.""" if "a" in Nvol: diff --git a/tests/compare_tests.py b/tests/compare_tests.py index 6bf675cf..e80afff2 100644 --- a/tests/compare_tests.py +++ b/tests/compare_tests.py @@ -53,10 +53,10 @@ def test_compare(Dirs, tol): varDataNew = newData[varKey][...] varDataRef = refData[varKey][...] diffMat = np.abs(varDataNew - varDataRef) - except ValueError: - assert False, "Fail from ValueError" - except KeyError: - assert False, "Fail from KeyError" + except ValueError as exception: + assert False, "Fail from ValueError: %s" % exception + except KeyError as exception: + assert False, "Fail from KeyError: %s" % exception # #Check absolute and relative error against tol assert np.mean(diffMat) < tol or \ diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_a.cfg b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_a.cfg new file mode 100644 index 00000000..b2af2ef7 --- /dev/null +++ b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_a.cfg @@ -0,0 +1,28 @@ +[Particles] +type = CHR +discretization = 2.e-7 +shape = sphere +thickness = 20e-9 + +[Material] +muRfunc = LiC6_LIONSIMBA +noise = false +noise_prefac = 1e-6 +numnoise = 200 +kappa = 4.0e-7 +B = 0.0 +rho_s = 1.839e28 +D = 3.9e-14 +Dfunc = constant +E_D = 5000 +dgammadc = 0e-30 +cwet = 0.98 + +[Reactions] +rxnType = BV_mod01 +k0 = 4.690 +E_A = 5000 +alpha = 0.5 +lambda = 6.26e-20 +Rfilm = 0e-0 + diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_c.cfg b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_c.cfg new file mode 100644 index 00000000..cbadb4c7 --- /dev/null +++ b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_c.cfg @@ -0,0 +1,28 @@ +[Particles] +type = CHR +discretization = 2.e-7 +shape = sphere +thickness = 20e-9 + +[Material] +muRfunc = LiCoO2_LIONSIMBA +noise = false +noise_prefac = 1e-6 +numnoise = 200 +kappa = 5.0148e-10 +B = 0.1916e9 +rho_s = 3.1036e28 +D = 1e-14 +Dfunc = constant +E_D = 5000 +dgammadc = 0e-30 +cwet = 0.98 + +[Reactions] +rxnType = BV_mod01 +k0 = 3.671 +E_A = 5000 +alpha = 0.5 +lambda = 6.26e-20 +Rfilm = 0e-0 + diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_system.cfg b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_system.cfg new file mode 100644 index 00000000..5c291691 --- /dev/null +++ b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/params_system.cfg @@ -0,0 +1,91 @@ +[Sim Params] +profileType = CC +Crate = 1 +1C_current_density = 30 +Vmax = 5 +Vmin = 2.5 +Vset = 0.12 +segments = [ + (0.3, 0.4), + (-0.5, 0.1), + ] +prevDir = false +tend = 1.2e3 +tsteps = 25 +relTol = 1e-6 +absTol = 1e-6 +T = 323 +randomSeed = false +nonisothermal = true +seed = 0 +Rser = 0. +Nvol_c = 10 +Nvol_s = 10 +Nvol_a = 10 +Npart_c = 1 +Npart_a = 1 + +[Electrodes] +cathode = params_c.cfg +anode = params_a.cfg +k0_foil = 1e0 +Rfilm_foil = 0e-0 + +[Particles] +mean_c = 2e-6 +stddev_c = 0 +mean_a = 2e-6 +stddev_a = 0 +cs0_c = 0.4995 +cs0_a = 0.8551 + +[Conductivity] +simBulkCond_c = true +simBulkCond_a = true +sigma_s_c = 412.43 +sigma_s_a = 685.77 +simPartCond_c = false +simPartCond_a = false +G_mean_c = 1e-14 +G_stddev_c = 0 +G_mean_a = 1e-14 +G_stddev_a = 0 + +[Geometry] +L_c = 8e-5 +L_a = 8.8e-5 +L_s = 2.5e-5 +P_L_c = 0.9593 +P_L_a = 0.9367 +poros_c = 0.385 +poros_a = 0.485 +poros_s = 0.724 +BruggExp_c = -3 +BruggExp_a = -3 +BruggExp_s = -3 + +[Thermal Parameters] +cp_c = 700 +cp_s = 700 +cp_a = 700 +rhom_c = 2500 +rhom_s = 1100 +rhom_a = 2500 +h_h = 10 +k_h_c = 2.1 +k_h_a = 1.7 +k_h_s = 0.16 +entropy_heat_gen = False + +[Electrolyte] +c0 = 1000 +zp = 1 +zm = -1 +nup = 1 +num = 1 +elyteModelType = SM +SMset = LIONSIMBA_nonisothermal +n = 1 +sp = -1 +Dp = 7.5e-10 +Dm = 7.5e-10 diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/commit.diff b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/commit.diff new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/commit.diff @@ -0,0 +1 @@ + diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/daetools_config_options.txt b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/daetools_config_options.txt new file mode 100644 index 00000000..7320f92e --- /dev/null +++ b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/daetools_config_options.txt @@ -0,0 +1,138 @@ +{ + "daetools": { + "core": { + "checkForInfiniteNumbers": "false", + "eventTolerance": "1E-7", + "logIndent": " ", + "pythonIndent": " ", + "checkUnitsConsistency": "true", + "resetLAMatrixAfterDiscontinuity": "true", + "printInfo": "false", + "nodes": { + "useNodeMemoryPools": "false", + "deleteNodesThreshold": "1000000", + "significantDecimalsForConstantsHash": "10" + }, + "equations": { + "info": [ + "If simplifyExpressions is true equation expressions will be simplified.", + "evaluationMode specifies the mode of evaluation: evaluationTree_OpenMP or computeStack_OpenMP.", + "computeStack_External is set by specifying the adComputeStackEvaluator_t object to the simulation.", + "If numThreads is 0 the default number of threads will be used (typically the number of cores in the system)." + ], + "simplifyExpressions": "false", + "evaluationMode": "computeStack_OpenMP", + "evaluationTree_OpenMP": { + "numThreads": "0" + }, + "computeStack_OpenMP": { + "numThreads": "0" + } + } + }, + "activity": { + "printHeader": "false", + "printStats": "false", + "timeHorizon": "100.0", + "reportingInterval": "1.0", + "reportTimeDerivatives": "false", + "reportSensitivities": "false", + "stopAtModelDiscontinuity": "true", + "reportDataAroundDiscontinuities": "true", + "objFunctionAbsoluteTolerance": "1E-8", + "constraintsAbsoluteTolerance": "1E-8", + "measuredVariableAbsoluteTolerance": "1E-8" + }, + "datareporting": { + "tcpipDataReceiverAddress": "127.0.0.1", + "tcpipDataReceiverPort": "50000", + "tcpipNumberOfRetries": "10", + "tcpipRetryAfterMilliSecs": "1000" + }, + "logging": { + "tcpipLogAddress": "127.0.0.1", + "tcpipLogPort": "51000" + }, + "minlpsolver": { + "printInfo": "false" + }, + "IDAS": { + "relativeTolerance": "1E-5", + "integrationMode": "Normal", + "reportDataInOneStepMode": "false", + "nextTimeAfterReinitialization": "1E-7", + "printInfo": "false", + "numberOfSTNRebuildsDuringInitialization": "1000", + "SensitivitySolutionMethod": "Staggered", + "SensErrCon": "false", + "sensRelativeTolerance": "1E-5", + "sensAbsoluteTolerance": "1E-5", + "MaxOrd": "5", + "MaxNumSteps": "1000", + "InitStep": "0.0", + "MaxStep": "0.0", + "MaxErrTestFails": "10", + "MaxNonlinIters": "4", + "MaxConvFails": "10", + "NonlinConvCoef": "0.33", + "SuppressAlg": "false", + "NoInactiveRootWarn": "false", + "NonlinConvCoefIC": "0.0033", + "MaxNumStepsIC": "5", + "MaxNumJacsIC": "4", + "MaxNumItersIC": "10", + "LineSearchOffIC": "false", + "gmres": { + "kspace": "30", + "EpsLin": "0.05", + "JacTimesVecFn": "DifferenceQuotient", + "DQIncrementFactor": "1.0", + "MaxRestarts": "5", + "GSType": "MODIFIED_GS" + } + }, + "superlu": { + "factorizationMethod": "SamePattern_SameRowPerm", + "useUserSuppliedWorkSpace": "false", + "workspaceSizeMultiplier": "3.0", + "workspaceMemoryIncrement": "1.5" + }, + "superlu_mt": { + "numThreads": "0" + }, + "intel_pardiso": { + "numThreads": "0" + }, + "BONMIN": { + "IPOPT": { + "print_level": "0", + "tol": "1E-5", + "linear_solver": "mumps", + "hessianApproximation": "limited-memory", + "mu_strategy": "adaptive" + } + }, + "NLOPT": { + "printInfo": "false", + "xtol_rel": "1E-6", + "xtol_abs": "1E-6", + "ftol_rel": "1E-6", + "ftol_abs": "1E-6", + "constr_tol": "1E-6" + }, + "deal_II": { + "printInfo": "false", + "assembly": { + "info": [ + "parallelAssembly can be: Sequential or OpenMP.", + "If numThreads is 0 the default number of threads will be used (typically the number of cores in the system).", + "queueSize specifies the size of the internal queue; when this size is reached the local data are copied to the global matrices." + ], + "parallelAssembly": "OpenMP", + "numThreads": "0", + "queueSize": "32" + } + } + } +} + diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_anode.p b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_anode.p new file mode 100644 index 00000000..d827e26d Binary files /dev/null and b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_anode.p differ diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_cathode.p b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_cathode.p new file mode 100644 index 00000000..9539343b Binary files /dev/null and b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_cathode.p differ diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_derived_values.p b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_derived_values.p new file mode 100644 index 00000000..0ba05e8e Binary files /dev/null and b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_derived_values.p differ diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_system.p b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_system.p new file mode 100644 index 00000000..c254b864 Binary files /dev/null and b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/input_dict_system.p differ diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/output_data b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/output_data new file mode 100644 index 00000000..e69de29b diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/output_data.mat b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/output_data.mat new file mode 100644 index 00000000..532c4c34 Binary files /dev/null and b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/output_data.mat differ diff --git a/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/run_info.txt b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/run_info.txt new file mode 100644 index 00000000..d7fac5f2 --- /dev/null +++ b/tests/ref_outputs/benchmark_LIONSIMBA_nonisothermal_with_heat_gen/sim_output/run_info.txt @@ -0,0 +1,17 @@ +mpet version: +0.1.8 + +branch name: +feature/temperature_effects + +commit hash: +f76687a + +to run, from the root repo directory, copy relevant files there, +edit input_params_system.cfg to point to correct material +params files, and: +$ git checkout [commit hash] +$ patch -p1 < commit.diff: +$ python[3] mpetrun.py input_params_system.cfg + +Total run time: 0.7441375255584717 s