diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 649c8ffa..73c1e8e2 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,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'] diff --git a/mpet/config/constants.py b/mpet/config/constants.py index 63a26c5b..b3a4f689 100644 --- a/mpet/config/constants.py +++ b/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, diff --git a/mpet/config/derived_values.py b/mpet/config/derived_values.py index 4bf88db9..5032c71c 100644 --- a/mpet/config/derived_values.py +++ b/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 diff --git a/mpet/config/parameterset.py b/mpet/config/parameterset.py index d14fa346..95e7e401 100644 --- a/mpet/config/parameterset.py +++ b/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}' diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py index 5c581bce..c305281b 100644 --- a/mpet/config/schemas.py +++ b/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, diff --git a/mpet/geometry.py b/mpet/geometry.py index 2b3fc98d..75020e45 100644 --- a/mpet/geometry.py +++ b/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 diff --git a/mpet/mod_cell.py b/mpet/mod_cell.py index 2a437ada..446bb6b8 100644 --- a/mpet/mod_cell.py +++ b/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 diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index f17718b9..67a448f7 100644 --- a/mpet/mod_electrodes.py +++ b/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 diff --git a/mpet/utils.py b/mpet/utils.py index 311344a6..192d9d48 100644 --- a/mpet/utils.py +++ b/mpet/utils.py @@ -85,6 +85,26 @@ def get_asc_vec(var, Nvol, dt=False): return out +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 ["a", "c"]: + # If it's an array of dae variable objects + out = config['rhom'][sectn]*(1-config["poros"][sectn])**(1-config["BruggExp"][sectn]) \ + * config['cp'][sectn] + config['rhom']['l'] * \ + config["poros"][sectn]**config["BruggExp"][sectn]*config['cp']['l'] + varout[sectn] = get_const_vec(out, Nvol[sectn]) + else: + out = config['rhom']['l'] * \ + config["poros"][sectn]**config["BruggExp"][sectn]*config['cp']['l'] + varout[sectn] = get_const_vec(out, Nvol[sectn]) + # 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: