added nonisothermal heat generation

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,10 +480,11 @@ def _scale_system_parameters(self, theoretical_1C_current):
         self['Dp'] = self['Dp'] / self['D_ref']
         self['Dm'] = self['Dm'] / self['D_ref']
         self['k_h'] = self['k_h'] / self['k_h_ref']
-        self['cp'] = self['cp'] / (self['k_h_ref'] * self['t_ref'] / self['L_ref']**3)
+        self['cp_l'] = self['cp_l'] / \
+            (self['k_h_ref'] * self['t_ref'] / (self['rho_ref']*self['L_ref']**2))
+        self['rhom_l'] = self['rhom_l'] / self['rho_ref']
         self['h_h'] = self['h_h'] * self['L_ref'] / self['k_h_ref']
-        self['sigma_lyte'] = self['sigma_lyte'] * self['k_h_ref'] * constants.T_ref / \
-            (constants.e**2 * self['D_ref'] * constants.N_A * constants.c_ref)
+        self['sigma_l'] = self['sigma_l'] / self['sigma_s_ref']
         self['c0'] = self['c0'] / constants.c_ref
         self['phi_cathode'] = 0.  # TODO: why is this defined if always 0?
         self['currset'] = self['currset'] / (theoretical_1C_current * self['curr_ref'])
@@ -504,6 +514,11 @@ def _scale_electrode_parameters(self):
                 if value is not None:
                     self[trode, param] = value / kT
 
+        for mat in self['materials']:
+            self['cp'][mat] = self['cp'][mat] / \
+                (self['k_h_ref'] * self['t_ref'] / (self['rho_ref']*self['L_ref']**2))
+            self['rhom'][mat] = self['rhom'][mat] / self['rho_ref']
+
         # scalings on separator
         if self['have_separator']:
             self['L']['s'] /= self['L_ref']

mpet/config/constants.py

@@ -20,9 +20,11 @@
 #: 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']
 #: subset of ``PARAMS_PER_TRODE``` that is defined for the separator as well
 PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp']
+# PARAMETERS THAT ARE NEEDED IN A THERMAL MODEL FOR ELECTROLYTE PROPERTIES
+PARAMS_ELYTE = ['cp', 'sigma', '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

@@ -147,6 +147,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 / self.config['L_ref']**3
+
     def sigma_s_ref(self):
         """Reference conductivity
         """
@@ -196,10 +202,8 @@ def D_ref(self):
     def k_h_ref(self):
         """Reference heat transfer coefficient
         """
-        if self.config['elyteModelType'] == 'dilute':
-            return self.config['k_h']
-        else:
-            return getattr(props_elyte, self.config['SMset'])()[-1]
+        return constants.c_ref * constants.k * \
+            self.config['L_ref']**2 / (self.config['t_ref'] * constants.N_A)
 
     def z(self):
         """Electrode capacity ratio

mpet/config/parameterset.py

@@ -2,7 +2,7 @@
 import configparser
 
 from mpet.config import schemas
-from mpet.config.constants import PARAMS_PER_TRODE, PARAMS_SEPARATOR
+from mpet.config.constants import PARAMS_PER_TRODE, PARAMS_SEPARATOR, PARAMS_ELYTE
 from mpet.exceptions import UnknownParameterError
 
 
@@ -84,6 +84,8 @@ def __getitem__(self, item):
             trodes = self['trodes'][:]  # make a copy here to avoid adding values to the original
             if item in PARAMS_SEPARATOR and self['have_separator']:
                 trodes.append('s')
+            if item in PARAMS_ELYTE:
+                trodes.append('l')
             for trode in trodes:
                 # get the value for this electrode/separator and store it
                 key = f'{item}_{trode}'

mpet/config/schemas.py

@@ -114,6 +114,15 @@ def tobool(value):
                        'BruggExp_c': Use(float),
                        'BruggExp_a': Use(float),
                        'BruggExp_s': Use(float)},
+          'Thermal Parameters': {Optional('cp_c', default=1e8): Use(float),
+                                 Optional('cp_a', default=1e8): Use(float),
+                                 Optional('cp_l', default=1e8): Use(float),
+                                 Optional('rhom_c', default=0.2): Use(float),
+                                 Optional('rhom_a', default=0.2): Use(float),
+                                 Optional('rhom_l', default=0.2): Use(float),
+                                 Optional('k_h', default=0.2): Use(float),
+                                 Optional('h_h', default=500): Use(float),
+                                 Optional('sigma_l', default=500): Use(float)},
           'Electrolyte': {'c0': Use(float),
                           'zp': Use(int),
                           'zm': And(Use(int), lambda x: x < 0),
@@ -124,11 +133,7 @@ def tobool(value):
                           'n': Use(int),
                           'sp': Use(int),
                           'Dp': Use(float),
-                          'Dm': Use(float),
-                          Optional('cp', default=1e8): Use(float),
-                          Optional('sigma_lyte', default=0.2): Use(float),
-                          Optional('k_h', default=0.2): Use(float),
-                          Optional('h_h', default=500): Use(float)}}
+                          'Dm': Use(float)}}
 
 #: Electrode parameters, per section
 electrode = {'Particles': {'type': lambda x: check_allowed_values(x,

mpet/geometry.py

@@ -85,6 +85,7 @@ def get_elyte_disc(Nvol, L, poros, BruggExp):
 
     # Distance between cell centers
     dxtmp = np.hstack((out["dxvec"][0], out["dxvec"], out["dxvec"][-1]))
+    out["dx"] = dxtmp
     out["dxd1"] = utils.mean_linear(dxtmp)
 
     # The porosity vector

mpet/mod_cell.py

@@ -59,6 +59,7 @@ 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:
@@ -83,6 +84,10 @@ 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")
@@ -210,6 +215,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]):
@@ -312,15 +334,19 @@ def DeclareEquations(self):
             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)
+            phibulkvec = utils.get_asc_vec(self.phi_bulk, 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]))
+            phibulktmp = np.hstack((self.phi_cell(), phibulkvec, config["phi_cathode"]))
 
             Nm_edges, i_edges, q_edges = get_lyte_internal_fluxes(ctmp, phitmp, Ttmp, disc, config)
 
@@ -374,7 +400,7 @@ def DeclareEquations(self):
             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(cvec, Tvec, i_edges, disc, config)
+            q_ohm = get_ohmic_heat(ctmp, Ttmp, phibulktmp, phitmp, 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))
@@ -386,8 +412,8 @@ def DeclareEquations(self):
                 if config['nonisothermal']:
                     # if heat generation is turned on. per volume.
                     eq = self.CreateEquation("lyte_energy_cons_vol{vInd}".format(vInd=vInd))
-                    eq.Residual = disc["porosvec"][vInd]*config["cp"] * \
-                        dTdtvec[vInd] - dvgq[vInd] - q_ohm[vInd]
+                    eq.Residual = rhocp_vec[vInd] * dTdtvec[vInd] - \
+                        dvgq[vInd] - q_ohm[vInd] - Qvvec[vInd]
                 else:
                     # if heat generation is turned off
                     eq = self.CreateEquation("lyte_energy_cons_vol{vInd}".format(vInd=vInd))
@@ -575,20 +601,34 @@ def get_lyte_internal_fluxes(c_lyte, phi_lyte, T_lyte, disc, config):
     return Nm_edges_int, i_edges_int, q_edges_int
 
 
-def get_ohmic_heat(c_lyte, T_lyte, i_edges_int, disc, config):
-    eps_o_tau = disc["eps_o_tau"][1:-1]
+def get_ohmic_heat(c_lyte, T_lyte, phi_lyte, phi_bulk, disc, config, Nvol):
+    eps_o_tau = disc["eps_o_tau"]
+    dx = disc["dx"][1:-1]
 
-    # get average current
-    i_cent = utils.mean_linear(i_edges_int)
+    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)
+    phi_lyte_int = utils.weighted_linear_mean(phi_lyte, wt)
+    phi_bulk_int = utils.weighted_linear_mean(phi_bulk, wt)
+    c_mid = c_lyte[1:-1]
+    T_mid = T_lyte[1:-1]
 
     # initialize sigma
-    sigma = 0
+    q_ohmic = 0
 
     if config["elyteModelType"] == "dilute":
-        sigma = eps_o_tau * config["sigma_lyte"]
+        sigma_l = eps_o_tau * config["sigma_l"]
     elif config["elyteModelType"] == "SM":
+        tp0 = getattr(props_elyte,config["SMset"])()[-3]
+
         sigma_fs = getattr(props_elyte,config["SMset"])()[1]
         # Get diffusivity and conductivity at cell edges using weighted harmonic mean
-        sigma = eps_o_tau*sigma_fs(c_lyte, T_lyte)
-    q_ohmic = i_cent**2/sigma
+        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)) * \
+            np.diff(np.log(c_edges_int))/dx*np.diff(phi_lyte_int)/dx
+    # this is going to be dra
+    sigma_s = (1-eps_o_tau[1:-1]) * sigma_s
+    q_ohmic = q_ohmic + sigma_s*(np.diff(phi_bulk_int)/dx)**2 + \
+        sigma_l*(np.diff(phi_lyte_int)/dx)**2
+    # do we have to extrapolate these
     return q_ohmic

mpet/mod_electrodes.py

@@ -55,6 +55,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")
@@ -138,16 +142,37 @@ 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, ISfuncs, noises)
+            eta1, eta2, c_surf1, c_surf2 = self.sld_dynamics_1D2var(c1, c2, mu_O, act_lyte,
+                                                                    ISfuncs, 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, ISfuncs, noises)
+            eta1, eta2, c_surf1, c_surf2 = self.sld_dynamics_0D2var(c1, c2, mu_O, act_lyte,
+                                                                    ISfuncs, noises)
+
+        # Define average rate of heat generation
+        eq = self.CreateEquation("q_rxn_bar")
+        eq.Residual = self.q_rxn_bar() - 0.5 * self.dcbar1dt() * (-eta1
+                                                                  - (np.log(c_surf1/(1-c_surf1))+1/self.c_lyte())) \
+            - 0.5 * self.dcbar2dt() * (-eta2
+                                       - (np.log(c_surf2/(1-c_surf2))+1/self.c_lyte()))
 
         for eq in self.Equations:
             eq.CheckUnitsConsistency = False
@@ -183,6 +208,7 @@ def sld_dynamics_0D2var(self, c1, c2, muO, act_lyte, ISfuncs, 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, eta2, c1_surf[-1], c2_surf[-1]
 
     def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, ISfuncs, noises):
         N = self.get_trode_param("N")
@@ -274,6 +300,11 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, ISfuncs, 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,
@@ -296,6 +327,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:
@@ -368,10 +401,15 @@ 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.ISfuncs, self.noise)
+            eta, c_surf = self.sld_dynamics_1D1var(c, mu_O, act_lyte, self.ISfuncs, 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.ISfuncs, self.noise)
+            eta, c_surf = self.sld_dynamics_0D1var(c, mu_O, act_lyte, self.ISfuncs, self.noise)
+
+        # Define average rate of heat generation
+        eq = self.CreateEquation("q_rxn_bar")
+        eq.Residual = self.q_rxn_bar() - self.dcbardt() * \
+            (-eta - (np.log(c_surf/(1-c_surf))+1/self.c_lyte()))
 
         for eq in self.Equations:
             eq.CheckUnitsConsistency = False
@@ -394,6 +432,8 @@ def sld_dynamics_0D1var(self, c, muO, act_lyte, ISfuncs, noise):
         eq = self.CreateEquation("dcsdt")
         eq.Residual = self.c.dt(0) - self.get_trode_param("delta_L")*self.Rxn()
 
+        return eta, c_surf[-1]
+
     def sld_dynamics_1D1var(self, c, muO, act_lyte, ISfuncs, noise):
         N = self.get_trode_param("N")
         # Equations for concentration evolution
@@ -462,6 +502,11 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, ISfuncs, 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