Feature/temperature effects

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)')
 ])
 

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

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

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

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

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):
         """

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,

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

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),

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)

mpet/electrode/materials/LiC6.py

@@ -1,17 +1,17 @@
 import numpy as np
 
 
-def LiC6(self, y, ybar, muR_ref, ISfuncs=(None, None)):
+def LiC6(self, y, ybar, T, muR_ref, ISfuncs=(None, None)):
     """ 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"), ISfuncs)
     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)

mpet/electrode/materials/LiC6_1param.py

@@ -1,12 +1,12 @@
 import numpy as np
 
 
-def LiC6_1param(self, y, ybar, muR_ref, ISfuncs=None):
+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"), ISfuncs)
+        y, T, self.get_trode_param("Omega_a"), self.get_trode_param("Omega_b"), 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