From 66993583fe2db982605996f2d09bf6f12ed64411 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Wed, 30 Nov 2022 15:31:07 +0000 Subject: [PATCH 01/27] start of 2d lfp --- configs/params_LFP2D.cfg | 31 +++++ mpet/config/configuration.py | 2 +- mpet/config/constants.py | 2 +- mpet/mod_cell.py | 4 +- mpet/mod_electrodes.py | 224 +++++++++++++++++++++++++++++++++++ mpet/sim.py | 11 +- 6 files changed, 270 insertions(+), 4 deletions(-) create mode 100644 configs/params_LFP2D.cfg diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg new file mode 100644 index 00000000..47c8699c --- /dev/null +++ b/configs/params_LFP2D.cfg @@ -0,0 +1,31 @@ +# Default parameters for simulating LFP in 1D using the ACR model. +# See params_electrodes.cfg for parameter explanations. + +[Particles] +type = ACR2D +discretization = 1e-9 +shape = C3 +thickness = 20e-9 + +[Material] +muRfunc = LiFePO4 +noise = false +noise_prefac = 1e-6 +numnoise = 200 +Omega_a = 1.8560e-20 +kappa = 5.0148e-10 +B = 0.1916e9 +rho_s = 1.3793e28 +D = 5.3e-19 +Dfunc = lattice +dgammadc = 0e-30 +cwet = 0.98 + +[Reactions] +rxnType = BV +k0 = 1.6e-1 +E_A = 13000 +alpha = 0.5 +# Fraggedakis et al. 2020, lambda = 8.3kBT +lambda = 3.4113e-20 +Rfilm = 0e-0 diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index f83aa316..078d4942 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -634,7 +634,7 @@ def _distr_part(self): # For particles with internal profiles, convert psd to # integers -- number of steps solidDisc = self[trode, 'discretization'] - if solidType in ['ACR']: + if solidType in ['ACR', 'ACR2D']: psd_num = np.ceil(raw / solidDisc).astype(int) psd_len = solidDisc * psd_num elif solidType in ['CHR', 'diffn', 'CHR2', 'diffn2']: diff --git a/mpet/config/constants.py b/mpet/config/constants.py index ed309c15..9579ec5c 100644 --- a/mpet/config/constants.py +++ b/mpet/config/constants.py @@ -13,7 +13,7 @@ #: General particle classification (1 var) two_var_types = ["diffn2", "CHR2", "homog2", "homog2_sdn"] #: General particle classification (2 var) -one_var_types = ["ACR", "diffn", "CHR", "homog", "homog_sdn"] +one_var_types = ["ACR", "diffn", "CHR", "homog", "homog_sdn", "ACR2D"] #: Reference concentration, mol/m^3 = 1M c_ref = 1000. #: Reaction rate epsilon for values close to zero diff --git a/mpet/mod_cell.py b/mpet/mod_cell.py index 46d244d2..ab456690 100644 --- a/mpet/mod_cell.py +++ b/mpet/mod_cell.py @@ -137,7 +137,9 @@ def __init__(self, config, Name, Parent=None, Description=""): dae.eOutletPort, self, "Bulk electrode port to particles") solidType = config[trode, "type"] - if solidType in constants.two_var_types: + if solidType == "ACR2D": + pMod = mod_electrodes.Mod2D + elif solidType in constants.two_var_types: pMod = mod_electrodes.Mod2var elif solidType in constants.one_var_types: pMod = mod_electrodes.Mod1var diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index bbb241db..203b9b34 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -23,6 +23,216 @@ from mpet.daeVariableTypes import mole_frac_t +class Mod2D(dae.daeModel): + def __init__(self, config, trode, vInd, pInd, + Name, Parent=None, Description=""): + super().__init__(Name, Parent, Description) + + self.config = config + self.trode = trode + self.ind = (vInd, pInd) + + # Domain + self.Dmn = dae.daeDomain("discretizationDomainX", self, dae.unit(), + "discretization domain in x direction") + self.Dmny = dae.daeDomain("discretizationDomainY", self, dae.unit(), + "discretization domain in y direction") + + # Variables + self.cy = {} + self.cybar = {} + self.cx = dae.daeVariable( + "cx", mole_frac_t, self, + "Concentration in x direction of active particle", [self.Dmn]) + + Nx = 50 + for k in range(Nx): + self.cy[k] = dae.daeVariable( + "cy_{k}".format(k=k), mole_frac_t, self, + "Concentration in y direction of element {k}".format(k=k), + [self.Dmny]) + for k in range(Nx): + self.cybar[k] = dae.daeVariable( + "cybar_{k}".format(k=k), mole_frac_t, self, + "Average concentration in y of element {k}".format(k=k)) + + self.cbar = dae.daeVariable( + "cbar", mole_frac_t, self, + "Average concentration in active particle") + # self.cxbar = dae.daeVariable( + # "cxbar", mole_frac_t, self, + # "Average concentration in x of active particle") + + + self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") + self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn]) + + # Get reaction rate function + rxnType = config[trode, "rxnType"] + self.calc_rxn_rate = utils.import_function(config[trode, "rxnType_filename"], + rxnType, + f"mpet.electrode.reactions.{rxnType}") + + # Ports + self.portInLyte = ports.portFromElyte( + "portInLyte", dae.eInletPort, self, "Inlet port from electrolyte") + self.portInBulk = ports.portFromBulk( + "portInBulk", dae.eInletPort, self, + "Inlet port from e- conducting phase") + self.phi_lyte = self.portInLyte.phi_lyte + self.c_lyte = self.portInLyte.c_lyte + self.phi_m = self.portInBulk.phi_m + + def get_trode_param(self, item): + """ + Shorthand to retrieve electrode-specific value + """ + value = self.config[self.trode, item] + # check if it is a particle-specific parameter + if item in self.config.params_per_particle: + value = value[self.ind] + return value + + def DeclareEquations(self): + dae.daeModel.DeclareEquations(self) + # Nx = self.get_trode_param("N") # number of grid points in particle + Nx = 50 + Ny = 10 + T = self.config["T"] # nondimensional temperature + # r_vec, volfrac_vec_x = geo.get_unit_solid_discr(self.get_trode_param('shape'), Nx) + # r_vec, volfrac_vec_y = geo.get_unit_solid_discr(self.get_trode_param('shape'), Ny) + + self.noise = None + + mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T, + "SM") + + # Define average filling fractions in particle + # eqx = self.CreateEquation("cxbar") + # for i in range(Nx): + # eqy = self.CreateEquation("cybar_{i}".format(i=i)) + # eqy.Residual = self.cx(i) + # for j in range(Ny): + # eqy.Residual -= self.cy[i](j)/Ny + + + for i in range(Nx): + eq = self.CreateEquation("cx_{i}".format(i=i)) + eq.Residual = self.cx(i) + for j in range(Ny): + eq.Residual -= self.cy[i](j)/Ny + + eq = self.CreateEquation("cbar") + eq.Residual = self.cbar() + for k in range(Nx): + eq.Residual -= self.cx(k)/Nx + + # Define average rate of filling of particle + eq = self.CreateEquation("dcbardt") + eq.Residual = self.dcbardt() + for k in range(Nx): + eq.Residual -= self.cx.dt(k)/Nx + + # cx = np.empty(Nx, dtype=object) + # cx[:] = [self.cx(k) for k in range(Nx)] + c_mat = np.empty((Ny, Nx), dtype=object) + for k in range(Nx): + for j in range(Ny): + c_mat[j, k] = self.cy[k](j) + + self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) + + for eq in self.Equations: + eq.CheckUnitsConsistency = False + + def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): + Ny = np.size(c_mat, 0) + Nx = np.size(c_mat, 1) + T = self.config["T"] + Dfunc_name = self.get_trode_param("Dfunc") + Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), + Dfunc_name, + f"mpet.electrode.diffusion.{Dfunc_name}") + # Equations for concentration evolution + # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors + + Mmat = get_Mmat("C3", Nx) + dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), Ny) + + + # muR, actR = calc_muR(c_mat[:,0], self.cbar(), self.config, self.trode, self.ind) + # muR_mat = np.empty((Nx,Ny), dtype=object) + # actR_mat = np.empty((Nx,Ny), dtype=object) + # muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), + # self.config, self.trode, self.ind) + # for k in range(Nx): + # muR_mat[k,:], actR_mat[k,:] = calc_muR(c_mat[k,:], self.cx(k), + # self.config, self.trode, self.ind) + + # print(muR_mat[1,:]) + + c_surf = c_mat[-1,:] + muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), + self.config, self.trode, self.ind) + + eta = calc_eta(muR_surf, muO) + eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm") + for i in range(Nx)]) + 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("lambda"), self.get_trode_param("alpha")) + + for i in range(Nx): + eq = self.CreateEquation("Rxn_{i}".format(i=i)) + eq.Residual = self.Rxn(i) - Rxn[i] + + RHSx = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(Nx)]) + dcdt_vec_x = np.empty(Nx, dtype=object) + dcdt_vec_x[0:Nx] = [self.cx.dt(k) for k in range(Nx)] + Mmatx = get_Mmat("C3", Nx) + LHS_vec_x = MX(Mmatx, dcdt_vec_x) + + for k in range(Nx): + eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) + eq.Residual = LHS_vec_x[k] - RHSx[k] + + # Flux_mat = np.empty((Nx,Ny), dtype=object) + RHS_mat = np.empty((Ny,Nx), dtype=object) + area_vec = 1 + LHS_mat = np.empty((Ny,Nx), dtype=object) + Mmaty = get_Mmat("C3", Ny) + # print(c_mat[1,:]) + for k in range(Nx): + Flux_bc = -self.Rxn(k) + # muR_vec = muR_mat[k,:] + # c_v = c_mat[:,k] + # c_vec = utils.get_var_vec(c_mat[:,k],Ny) + muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cx(k), + self.config, self.trode, self.ind) + + # Flux_mat[k,:] = calc_flux_CHR(c_mat[k,:], muR_vec, self.get_trode_param("D"), Dfunc, + # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + Flux_vec = calc_flux_C3ver(c_mat[k,:], muR_vec, self.get_trode_param("D"), Dfunc, + self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + # c_vec = utils.get_var_vec(c_mat[:,k],Ny) + + + # RHS_mat = -np.diff(Flux_vec * area_vec) + RHS_vec = -np.diff(Flux_vec * area_vec) + + dcdt_vec_y = np.empty(Ny, dtype=object) + dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] + + LHS_vec_y = MX(Mmaty, dcdt_vec_y) + for j in range(Ny): + eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) + eq.Residual = LHS_vec_y[j] - RHS_vec[j] + LHS_mat[:,k] = LHS_vec_y + + + + class Mod2var(dae.daeModel): def __init__(self, config, trode, vInd, pInd, Name, Parent=None, Description=""): @@ -502,6 +712,19 @@ def calc_flux_CHR(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): np.diff(mu + noise(dae.Time().Value))/dr return Flux_vec +def calc_flux_C3ver(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): + N = len(c) + Flux_vec = np.empty(N+1, dtype=object) + Flux_vec[0] = 0 # Symmetry at r=0 + Flux_vec[-1] = Flux_bc + c_edges = utils.mean_linear(c) + if noise is None: + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu)/dr + else: + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ + np.diff(mu + noise(dae.Time().Value))/dr + return Flux_vec + def calc_flux_CHR2(c1, c2, mu1_R, mu2_R, D, Dfunc, E_D, Flux1_bc, Flux2_bc, dr, T, noise1, noise2): N = len(c1) @@ -546,6 +769,7 @@ def MX(mat, objvec): if not isinstance(mat, sprs.csr.csr_matrix): raise Exception("MX function designed for csr mult") n = objvec.shape[0] + print(n) if isinstance(objvec[0], dae.pyCore.adouble): out = np.empty(n, dtype=object) else: diff --git a/mpet/sim.py b/mpet/sim.py index c9130bc4..d5f16102 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -53,8 +53,12 @@ def SetUpParametersAndDomains(self): self.m.DmnPart[tr].CreateArray(config["Npart"][tr]) for i in range(config["Nvol"][tr]): for j in range(config["Npart"][tr]): + # self.m.particles[tr][i, j].Dmn.CreateArray( + # int(config["psd_num"][tr][i,j])) self.m.particles[tr][i, j].Dmn.CreateArray( - int(config["psd_num"][tr][i,j])) + 50) + self.m.particles[tr][i, j].Dmny.CreateArray( + 10) def SetUpVariables(self): config = self.config @@ -85,6 +89,11 @@ def SetUpVariables(self): solidType = self.config[tr, "type"] if solidType in constants.one_var_types: part.cbar.SetInitialGuess(cs0) + if solidType == "ACR2D": + for k in range(50): + part.cx.SetInitialCondition(k, cs0) + for j in range(10): + part.cy[k].SetInitialCondition(j, cs0) for k in range(Nij): part.c.SetInitialCondition(k, cs0) elif solidType in constants.two_var_types: From d41f275d5633f7bdff4f665c8b5abb6575c51f4a Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 13:28:30 +0000 Subject: [PATCH 02/27] sometimes it's working --- configs/params_LFP2D.cfg | 3 +- mpet/config/configuration.py | 20 +- mpet/config/schemas.py | 3 +- mpet/geometry.py | 4 +- mpet/mod_electrodes.py | 67 +-- mpet/mod_electrodes_old.py | 785 +++++++++++++++++++++++++++++++++++ mpet/sim.py | 16 +- 7 files changed, 853 insertions(+), 45 deletions(-) create mode 100644 mpet/mod_electrodes_old.py diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index 47c8699c..e9f418ce 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -5,6 +5,7 @@ type = ACR2D discretization = 1e-9 shape = C3 +discretization_ver = 3e-9 thickness = 20e-9 [Material] @@ -16,7 +17,7 @@ Omega_a = 1.8560e-20 kappa = 5.0148e-10 B = 0.1916e9 rho_s = 1.3793e28 -D = 5.3e-19 +D = 5.3e-18 Dfunc = lattice dgammadc = 0e-30 cwet = 0.98 diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 078d4942..cdd2b503 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -615,6 +615,14 @@ def _distr_part(self): if not np.all(self['specified_psd'][trode]): # If PSD is not specified, make a length-sampled particle size distribution # Log-normally distributed + if self[trode,'type'] == 'ACR2D': + np.random.seed(self.D_s['seed']) + mean_t = self[trode,'thickness'] + stddev_t = mean_t*0.1 + var_t = stddev_t**2 + mu_t = np.log((mean_t**2) / np.sqrt(var_t + mean_t**2)) + sigma_t = np.sqrt(np.log(var_t/(mean_t**2) + 1)) + raw_t = np.random.lognormal(mu_t, sigma_t, size=(Nvol, Npart)) mean = self['mean'][trode] stddev = self['stddev'][trode] if np.allclose(stddev, 0., atol=1e-12): @@ -634,9 +642,14 @@ def _distr_part(self): # For particles with internal profiles, convert psd to # integers -- number of steps solidDisc = self[trode, 'discretization'] - if solidType in ['ACR', 'ACR2D']: + if solidType in ['ACR']: psd_num = np.ceil(raw / solidDisc).astype(int) psd_len = solidDisc * psd_num + elif solidType in ['ACR2D']: + solidDisc_ver = self[trode, 'discretization_ver'] + psd_num = np.ceil(raw/ solidDisc).astype(int) + psd_len = solidDisc * psd_num + psd_num_ver = np.ceil(raw_t/solidDisc_ver).astype(int) elif solidType in ['CHR', 'diffn', 'CHR2', 'diffn2']: psd_num = np.ceil(raw / solidDisc).astype(int) + 1 psd_len = solidDisc * (psd_num - 1) @@ -674,6 +687,9 @@ def _distr_part(self): self['psd_area'][trode] = psd_area self['psd_vol'][trode] = psd_vol self['psd_vol_FracVol'][trode] = psd_frac_vol + # store values to config 2d + if self[trode,'type'] in 'ACR2D': + self['psd_num_ver'][trode] = psd_num_ver def _G(self): """ @@ -727,6 +743,8 @@ def _indvPart(self): kappa_ref = constants.k * constants.T_ref * cs_ref_part * plen**2 # J/m gamma_S_ref = kappa_ref / plen # J/m^2 # non-dimensional quantities + if self[trode,'type'] in 'ACR2D': + self[trode, 'indvPart']['N_ver'][i,j] = self['psd_num_ver'][trode][i,j] if self[trode, 'kappa'] is not None: self[trode, 'indvPart']['kappa'][i, j] = self[trode, 'kappa'] / kappa_ref if self[trode, 'dgammadc'] is not None: diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py index 17e2ea76..1f607d91 100644 --- a/mpet/config/schemas.py +++ b/mpet/config/schemas.py @@ -133,7 +133,8 @@ def tobool(value): 'discretization': Use(float), 'shape': lambda x: check_allowed_values(x, ["C3", "sphere", "cylinder", "homog_sdn"]), - Optional('thickness'): Use(float)}, + Optional('discretization_ver', default = None):Use(float), + Optional('thickness'): Use(float)}, 'Material': {Optional('muRfunc_filename', default=None): str, 'muRfunc': str, 'noise': Use(tobool), diff --git a/mpet/geometry.py b/mpet/geometry.py index c35dbba1..e5b952f2 100644 --- a/mpet/geometry.py +++ b/mpet/geometry.py @@ -9,7 +9,9 @@ def get_unit_solid_discr(Shape, N): r_vec = None volfrac_vec = np.ones(1) elif Shape == "C3": - r_vec = None + # r_vec = None + Rs = 1. # (non-dimensionalized by itself) + r_vec = np.linspace(0, Rs, N) # For 1D particle, the vol fracs are simply related to the # length discretization volfrac_vec = (1./N) * np.ones(N) # scaled to 1D particle volume diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 203b9b34..45f9e39c 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -39,22 +39,21 @@ def __init__(self, config, trode, vInd, pInd, "discretization domain in y direction") # Variables - self.cy = {} - self.cybar = {} self.cx = dae.daeVariable( "cx", mole_frac_t, self, "Concentration in x direction of active particle", [self.Dmn]) - Nx = 50 + Nx = self.get_trode_param("N") + self.cy = {} for k in range(Nx): self.cy[k] = dae.daeVariable( - "cy_{k}".format(k=k), mole_frac_t, self, + "cy{k}".format(k=k), mole_frac_t, self, "Concentration in y direction of element {k}".format(k=k), [self.Dmny]) - for k in range(Nx): - self.cybar[k] = dae.daeVariable( - "cybar_{k}".format(k=k), mole_frac_t, self, - "Average concentration in y of element {k}".format(k=k)) + + # self.cybar[k] = dae.daeVariable( + # "cybar{k}".format(k=k), mole_frac_t, self, + # "Average concentration in y of element {k}".format(k=k)) self.cbar = dae.daeVariable( "cbar", mole_frac_t, self, @@ -96,7 +95,7 @@ def get_trode_param(self, item): def DeclareEquations(self): dae.daeModel.DeclareEquations(self) # Nx = self.get_trode_param("N") # number of grid points in particle - Nx = 50 + Nx = self.get_trode_param("N") Ny = 10 T = self.config["T"] # nondimensional temperature # r_vec, volfrac_vec_x = geo.get_unit_solid_discr(self.get_trode_param('shape'), Nx) @@ -116,11 +115,11 @@ def DeclareEquations(self): # eqy.Residual -= self.cy[i](j)/Ny - for i in range(Nx): - eq = self.CreateEquation("cx_{i}".format(i=i)) - eq.Residual = self.cx(i) + for k in range(Nx): + eq = self.CreateEquation("avgs_cy_cx{k}".format(k=k)) + eq.Residual = self.cx(k) for j in range(Ny): - eq.Residual -= self.cy[i](j)/Ny + eq.Residual -= self.cy[k](j)/Ny eq = self.CreateEquation("cbar") eq.Residual = self.cbar() @@ -137,8 +136,7 @@ def DeclareEquations(self): # cx[:] = [self.cx(k) for k in range(Nx)] c_mat = np.empty((Ny, Nx), dtype=object) for k in range(Nx): - for j in range(Ny): - c_mat[j, k] = self.cy[k](j) + c_mat[:,k] = [self.cy[k](j) for j in range(Ny)] self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) @@ -156,10 +154,8 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): # Equations for concentration evolution # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors - Mmat = get_Mmat("C3", Nx) - dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), Ny) - - + dr, edges = geo.get_dr_edges("C3", Ny) + # dr = 1 # muR, actR = calc_muR(c_mat[:,0], self.cbar(), self.config, self.trode, self.ind) # muR_mat = np.empty((Nx,Ny), dtype=object) # actR_mat = np.empty((Nx,Ny), dtype=object) @@ -187,34 +183,40 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq = self.CreateEquation("Rxn_{i}".format(i=i)) eq.Residual = self.Rxn(i) - Rxn[i] - RHSx = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(Nx)]) - dcdt_vec_x = np.empty(Nx, dtype=object) - dcdt_vec_x[0:Nx] = [self.cx.dt(k) for k in range(Nx)] - Mmatx = get_Mmat("C3", Nx) - LHS_vec_x = MX(Mmatx, dcdt_vec_x) + # RHSx = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(Nx)]) + # dcdt_vec_x = np.empty(Nx, dtype=object) + # dcdt_vec_x[0:Nx] = [self.cx.dt(k) for k in range(Nx)] + # Mmatx = get_Mmat("C3", Nx) + # LHS_vec_x = MX(Mmatx, dcdt_vec_x) - for k in range(Nx): - eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) - eq.Residual = LHS_vec_x[k] - RHSx[k] + # for k in range(Nx): + # eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) + # eq.Residual = LHS_vec_x[k] - RHSx[k] # Flux_mat = np.empty((Nx,Ny), dtype=object) - RHS_mat = np.empty((Ny,Nx), dtype=object) + # RHS_mat = np.empty((Ny,Nx), dtype=object) area_vec = 1 - LHS_mat = np.empty((Ny,Nx), dtype=object) + # LHS_mat = np.empty((Ny,Nx), dtype=object) Mmaty = get_Mmat("C3", Ny) # print(c_mat[1,:]) for k in range(Nx): + # print(k, ' --- \n') Flux_bc = -self.Rxn(k) + # eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) + # eq.Residual = LHS_vec_x[k] - self.get_trode_param("delta_L")*self.Rxn(k) # muR_vec = muR_mat[k,:] # c_v = c_mat[:,k] # c_vec = utils.get_var_vec(c_mat[:,k],Ny) - muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cx(k), + muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cbar(), self.config, self.trode, self.ind) # Flux_mat[k,:] = calc_flux_CHR(c_mat[k,:], muR_vec, self.get_trode_param("D"), Dfunc, # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) - Flux_vec = calc_flux_C3ver(c_mat[k,:], muR_vec, self.get_trode_param("D"), Dfunc, + Flux_vec = calc_flux_C3ver(c_mat[:,k], muR_vec, self.get_trode_param("D"), Dfunc, self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + # Flux_vec = calc_flux_diffn(c_mat[:,k], self.get_trode_param("D"), Dfunc, + # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + # c_vec = utils.get_var_vec(c_mat[:,k],Ny) @@ -228,7 +230,7 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): for j in range(Ny): eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) eq.Residual = LHS_vec_y[j] - RHS_vec[j] - LHS_mat[:,k] = LHS_vec_y + # LHS_mat[:,k] = LHS_vec_y @@ -769,7 +771,6 @@ def MX(mat, objvec): if not isinstance(mat, sprs.csr.csr_matrix): raise Exception("MX function designed for csr mult") n = objvec.shape[0] - print(n) if isinstance(objvec[0], dae.pyCore.adouble): out = np.empty(n, dtype=object) else: diff --git a/mpet/mod_electrodes_old.py b/mpet/mod_electrodes_old.py new file mode 100644 index 00000000..0a00c2f9 --- /dev/null +++ b/mpet/mod_electrodes_old.py @@ -0,0 +1,785 @@ +"""These models define individual particles of active material. + +This includes the equations for both 1-parameter models and 2-parameters models defining + - mass conservation (concentration evolution) + - reaction rate at the surface of the particles +In each model class it has options for different types of particles: + - homogeneous + - Fick-like diffusion + - Cahn-Hilliard (with reaction boundary condition) + - Allen-Cahn (with reaction throughout the particle) +These models can be instantiated from the mod_cell module to simulate various types of active +materials within a battery electrode. +""" +import daetools.pyDAE as dae +import numpy as np +import scipy.sparse as sprs +import scipy.interpolate as sintrp + +import mpet.geometry as geo +import mpet.ports as ports +import mpet.props_am as props_am +import mpet.utils as utils +from mpet.daeVariableTypes import mole_frac_t + + +class Mod2D(dae.daeModel): + def __init__(self, config, trode, vInd, pInd, + Name, Parent=None, Description=""): + super().__init__(Name, Parent, Description) + + self.config = config + self.trode = trode + self.ind = (vInd, pInd) + + # Domain + self.Dmn = dae.daeDomain("discretizationDomainX", self, dae.unit(), + "discretization domain in x direction") + self.Dmny = dae.daeDomain("discretizationDomainY", self, dae.unit(), + "discretization domain in y direction") + + # Variables + self.cx = dae.daeVariable( + "cx", mole_frac_t, self, + "Concentration in x direction of active particle", [self.Dmn]) + + Nx = self.get_trode_param("N") + self.cy = {} + for k in range(Nx): + self.cy[k] = dae.daeVariable( + "cy{k}".format(k=k), mole_frac_t, self, + "Concentration in y direction of element {k}".format(k=k), + [self.Dmny]) + + # self.cybar[k] = dae.daeVariable( + # "cybar{k}".format(k=k), mole_frac_t, self, + # "Average concentration in y of element {k}".format(k=k)) + + self.cbar = dae.daeVariable( + "cbar", mole_frac_t, self, + "Average concentration in active particle") + # self.cxbar = dae.daeVariable( + # "cxbar", mole_frac_t, self, + # "Average concentration in x of active particle") + + + self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") + self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn]) + + # Get reaction rate function + rxnType = config[trode, "rxnType"] + self.calc_rxn_rate = utils.import_function(config[trode, "rxnType_filename"], + rxnType, + f"mpet.electrode.reactions.{rxnType}") + + # Ports + self.portInLyte = ports.portFromElyte( + "portInLyte", dae.eInletPort, self, "Inlet port from electrolyte") + self.portInBulk = ports.portFromBulk( + "portInBulk", dae.eInletPort, self, + "Inlet port from e- conducting phase") + self.phi_lyte = self.portInLyte.phi_lyte + self.c_lyte = self.portInLyte.c_lyte + self.phi_m = self.portInBulk.phi_m + + def get_trode_param(self, item): + """ + Shorthand to retrieve electrode-specific value + """ + value = self.config[self.trode, item] + # check if it is a particle-specific parameter + if item in self.config.params_per_particle: + value = value[self.ind] + return value + + def DeclareEquations(self): + dae.daeModel.DeclareEquations(self) + # Nx = self.get_trode_param("N") # number of grid points in particle + Nx = self.get_trode_param("N") + Ny = 10 + T = self.config["T"] # nondimensional temperature + # r_vec, volfrac_vec_x = geo.get_unit_solid_discr(self.get_trode_param('shape'), Nx) + # r_vec, volfrac_vec_y = geo.get_unit_solid_discr(self.get_trode_param('shape'), Ny) + + self.noise = None + + mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T, + "SM") + + # Define average filling fractions in particle + # eqx = self.CreateEquation("cxbar") + # for i in range(Nx): + # eqy = self.CreateEquation("cybar_{i}".format(i=i)) + # eqy.Residual = self.cx(i) + # for j in range(Ny): + # eqy.Residual -= self.cy[i](j)/Ny + + + for k in range(Nx): + eq = self.CreateEquation("avgs_cy_cx{k}".format(k=k)) + eq.Residual = self.cx(k) + for j in range(Ny): + eq.Residual -= self.cy[k](j)/Ny + + eq = self.CreateEquation("cbar") + eq.Residual = self.cbar() + for k in range(Nx): + eq.Residual -= self.cx(k)/Nx + + # Define average rate of filling of particle + eq = self.CreateEquation("dcbardt") + eq.Residual = self.dcbardt() + for k in range(Nx): + eq.Residual -= self.cx.dt(k)/Nx + + # cx = np.empty(Nx, dtype=object) + # cx[:] = [self.cx(k) for k in range(Nx)] + c_mat = np.empty((Ny, Nx), dtype=object) + for k in range(Nx): + c_mat[:,k] = [self.cy[k](j) for j in range(Ny)] + + self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) + + for eq in self.Equations: + eq.CheckUnitsConsistency = False + + def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): + Ny = np.size(c_mat, 0) + Nx = np.size(c_mat, 1) + T = self.config["T"] + Dfunc_name = self.get_trode_param("Dfunc") + Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), + Dfunc_name, + f"mpet.electrode.diffusion.{Dfunc_name}") + # Equations for concentration evolution + # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors + + dr, edges = geo.get_dr_edges("C3", Ny) + # dr = 1 + # muR, actR = calc_muR(c_mat[:,0], self.cbar(), self.config, self.trode, self.ind) + # muR_mat = np.empty((Nx,Ny), dtype=object) + # actR_mat = np.empty((Nx,Ny), dtype=object) + # muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), + # self.config, self.trode, self.ind) + # for k in range(Nx): + # muR_mat[k,:], actR_mat[k,:] = calc_muR(c_mat[k,:], self.cx(k), + # self.config, self.trode, self.ind) + + # print(muR_mat[1,:]) + + c_surf = c_mat[-1,:] + muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), + self.config, self.trode, self.ind) + + eta = calc_eta(muR_surf, muO) + eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm") + for i in range(Nx)]) + 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("lambda"), self.get_trode_param("alpha")) + + for i in range(Nx): + eq = self.CreateEquation("Rxn_{i}".format(i=i)) + eq.Residual = self.Rxn(i) - Rxn[i] + + RHSx = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(Nx)]) + dcdt_vec_x = np.empty(Nx, dtype=object) + dcdt_vec_x[0:Nx] = [self.cx.dt(k) for k in range(Nx)] + Mmatx = get_Mmat("C3", Nx) + LHS_vec_x = MX(Mmatx, dcdt_vec_x) + + for k in range(Nx): + eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) + eq.Residual = LHS_vec_x[k] - RHSx[k] + + # Flux_mat = np.empty((Nx,Ny), dtype=object) + # RHS_mat = np.empty((Ny,Nx), dtype=object) + area_vec = 1 + # LHS_mat = np.empty((Ny,Nx), dtype=object) + Mmaty = get_Mmat("C3", Ny) + # print(c_mat[1,:]) + for k in range(Nx): + # print(k, ' --- \n') + Flux_bc = -self.Rxn(k) + # muR_vec = muR_mat[k,:] + # c_v = c_mat[:,k] + # c_vec = utils.get_var_vec(c_mat[:,k],Ny) + muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cbar(), + self.config, self.trode, self.ind) + + # Flux_mat[k,:] = calc_flux_CHR(c_mat[k,:], muR_vec, self.get_trode_param("D"), Dfunc, + # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + Flux_vec = calc_flux_C3ver(c_mat[:,k], muR_vec, self.get_trode_param("D"), Dfunc, + self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + # Flux_vec = calc_flux_diffn(c_mat[:,k], self.get_trode_param("D"), Dfunc, + # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + + # c_vec = utils.get_var_vec(c_mat[:,k],Ny) + + + # RHS_mat = -np.diff(Flux_vec * area_vec) + RHS_vec = -np.diff(Flux_vec * area_vec) + + dcdt_vec_y = np.empty(Ny, dtype=object) + dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] + + LHS_vec_y = MX(Mmaty, dcdt_vec_y) + for j in range(Ny): + eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) + eq.Residual = LHS_vec_y[j] - RHS_vec[j] + # LHS_mat[:,k] = LHS_vec_y + + + + +class Mod2var(dae.daeModel): + def __init__(self, config, trode, vInd, pInd, + Name, Parent=None, Description=""): + super().__init__(Name, Parent, Description) + + self.config = config + self.trode = trode + self.ind = (vInd, pInd) + + # Domain + self.Dmn = dae.daeDomain("discretizationDomain", self, dae.unit(), + "discretization domain") + + # Variables + self.c1 = dae.daeVariable( + "c1", mole_frac_t, self, + "Concentration in 'layer' 1 of active particle", [self.Dmn]) + self.c2 = dae.daeVariable( + "c2", mole_frac_t, self, + "Concentration in 'layer' 2 of active particle", [self.Dmn]) + self.cbar = dae.daeVariable( + "cbar", mole_frac_t, self, + "Average concentration in active particle") + self.c1bar = dae.daeVariable( + "c1bar", mole_frac_t, self, + "Average concentration in 'layer' 1 of active particle") + self.c2bar = dae.daeVariable( + "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") + 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") + else: + self.Rxn1 = dae.daeVariable("Rxn1", dae.no_t, self, "Rate of reaction 1", [self.Dmn]) + self.Rxn2 = dae.daeVariable("Rxn2", dae.no_t, self, "Rate of reaction 2", [self.Dmn]) + + # Get reaction rate function + rxnType = config[trode, "rxnType"] + self.calc_rxn_rate = utils.import_function(config[trode, "rxnType_filename"], + rxnType, + f"mpet.electrode.reactions.{rxnType}") + + # Ports + self.portInLyte = ports.portFromElyte( + "portInLyte", dae.eInletPort, self, "Inlet port from electrolyte") + self.portInBulk = ports.portFromBulk( + "portInBulk", dae.eInletPort, self, + "Inlet port from e- conducting phase") + self.phi_lyte = self.portInLyte.phi_lyte + self.c_lyte = self.portInLyte.c_lyte + self.phi_m = self.portInBulk.phi_m + + def get_trode_param(self, item): + """ + Shorthand to retrieve electrode-specific value + """ + value = self.config[self.trode, item] + # check if it is a particle-specific parameter + if item in self.config.params_per_particle: + value = value[self.ind] + return value + + 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 + self.noise1 = self.noise2 = None + if self.get_trode_param("noise"): + numnoise = self.get_trode_param("numnoise") + noise_prefac = self.get_trode_param("noise_prefac") + tvec = np.linspace(0., 1.05*self.config["tend"], numnoise) + noise_data1 = noise_prefac*np.random.randn(numnoise, N) + noise_data2 = noise_prefac*np.random.randn(numnoise, N) + self.noise1 = sintrp.interp1d(tvec, noise_data1, axis=0, + bounds_error=False, fill_value=0.) + self.noise2 = sintrp.interp1d(tvec, noise_data2, axis=0, + bounds_error=False, fill_value=0.) + noises = (self.noise1, self.noise2) + + # 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.config["elyteModelType"]) + + # Define average filling fractions in particle + eq1 = self.CreateEquation("c1bar") + eq2 = self.CreateEquation("c2bar") + eq1.Residual = self.c1bar() + eq2.Residual = self.c2bar() + for k in range(N): + eq1.Residual -= self.c1(k) * volfrac_vec[k] + eq2.Residual -= self.c2(k) * volfrac_vec[k] + eq = self.CreateEquation("cbar") + eq.Residual = self.cbar() - .5*(self.c1bar() + self.c2bar()) + + # Define average rate of filling of particle + eq = self.CreateEquation("dcbardt") + eq.Residual = self.dcbardt() + for k in range(N): + eq.Residual -= .5*(self.c1.dt(k) + 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) + 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) + + 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, + self.trode, self.ind) + eta1 = calc_eta(mu1R_surf, muO) + eta2 = calc_eta(mu2R_surf, muO) + eta1_eff = eta1 + self.Rxn1()*self.get_trode_param("Rfilm") + eta2_eff = eta2 + self.Rxn2()*self.get_trode_param("Rfilm") + noise1, noise2 = noises + if self.get_trode_param("noise"): + eta1_eff += noise1(dae.Time().Value) + 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("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("lambda"), self.get_trode_param("alpha")) + eq1 = self.CreateEquation("Rxn1") + eq2 = self.CreateEquation("Rxn2") + eq1.Residual = self.Rxn1() - Rxn1[0] + eq2.Residual = self.Rxn2() - Rxn2[0] + + eq1 = self.CreateEquation("dc1sdt") + 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] + + 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) + dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), N) + + # 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) + c1_surf = c1[-1] + c2_surf = c2[-1] + mu1R_surf, act1R_surf = mu1R[-1], act1R[-1] + mu2R_surf, act2R_surf = mu2R[-1], act2R[-1] + eta1 = calc_eta(mu1R_surf, muO) + eta2 = calc_eta(mu2R_surf, muO) + if self.get_trode_param("type") in ["ACR2"]: + eta1_eff = np.array([eta1[i] + + self.Rxn1(i)*self.get_trode_param("Rfilm") for i in range(N)]) + eta2_eff = np.array([eta2[i] + + self.Rxn2(i)*self.get_trode_param("Rfilm") for i in range(N)]) + else: + eta1_eff = eta1 + self.Rxn1()*self.get_trode_param("Rfilm") + 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("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("lambda"), self.get_trode_param("alpha")) + if self.get_trode_param("type") in ["ACR2"]: + for i in range(N): + eq1 = self.CreateEquation("Rxn1_{i}".format(i=i)) + eq2 = self.CreateEquation("Rxn2_{i}".format(i=i)) + eq1.Residual = self.Rxn1(i) - Rxn1[i] + eq2.Residual = self.Rxn2(i) - Rxn2[i] + else: + eq1 = self.CreateEquation("Rxn1") + eq2 = self.CreateEquation("Rxn2") + eq1.Residual = self.Rxn1() - Rxn1 + eq2.Residual = self.Rxn2() - Rxn2 + + # Get solid particle fluxes (if any) and RHS + if self.get_trode_param("type") in ["diffn2", "CHR2"]: + # Positive reaction (reduction, intercalation) is negative + # flux of Li at the surface. + Flux1_bc = -0.5 * self.Rxn1() + Flux2_bc = -0.5 * self.Rxn2() + Dfunc_name = self.get_trode_param("Dfunc") + Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), + Dfunc_name, + f"mpet.electrode.diffusion.{Dfunc_name}") + if self.get_trode_param("type") == "CHR2": + 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) + if self.get_trode_param("shape") == "sphere": + area_vec = 4*np.pi*edges**2 + elif self.get_trode_param("shape") == "cylinder": + area_vec = 2*np.pi*edges # per unit height + RHS1 = -np.diff(Flux1_vec * area_vec) + RHS2 = -np.diff(Flux2_vec * area_vec) +# kinterlayer = 1e-3 +# interLayerRxn = (kinterlayer * (1 - c1) * (1 - c2) * (act1R - act2R)) +# RxnTerm1 = -interLayerRxn +# RxnTerm2 = interLayerRxn + RxnTerm1 = 0 + RxnTerm2 = 0 + RHS1 += RxnTerm1 + RHS2 += RxnTerm2 + + dc1dt_vec = np.empty(N, dtype=object) + dc2dt_vec = np.empty(N, dtype=object) + dc1dt_vec[0:N] = [self.c1.dt(k) for k in range(N)] + dc2dt_vec[0:N] = [self.c2.dt(k) for k in range(N)] + LHS1_vec = MX(Mmat, dc1dt_vec) + LHS2_vec = MX(Mmat, dc2dt_vec) + for k in range(N): + eq1 = self.CreateEquation("dc1sdt_discr{k}".format(k=k)) + eq2 = self.CreateEquation("dc2sdt_discr{k}".format(k=k)) + eq1.Residual = LHS1_vec[k] - RHS1[k] + eq2.Residual = LHS2_vec[k] - RHS2[k] + + +class Mod1var(dae.daeModel): + def __init__(self, config, trode, vInd, pInd, + Name, Parent=None, Description=""): + super().__init__(Name, Parent, Description) + + self.config = config + self.trode = trode + self.ind = (vInd, pInd) + + # Domain + self.Dmn = dae.daeDomain("discretizationDomain", self, dae.unit(), + "discretization domain") + + # Variables + self.c = dae.daeVariable("c", mole_frac_t, self, + "Concentration in active particle", + [self.Dmn]) + self.cbar = dae.daeVariable( + "cbar", mole_frac_t, self, + "Average concentration in active particle") + self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") + if config[trode, "type"] not in ["ACR"]: + self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction") + else: + self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn]) + + # Get reaction rate function + rxnType = config[trode, "rxnType"] + self.calc_rxn_rate = utils.import_function(config[trode, "rxnType_filename"], + rxnType, + f"mpet.electrode.reactions.{rxnType}") + + # Ports + self.portInLyte = ports.portFromElyte( + "portInLyte", dae.eInletPort, self, + "Inlet port from electrolyte") + self.portInBulk = ports.portFromBulk( + "portInBulk", dae.eInletPort, self, + "Inlet port from e- conducting phase") + self.phi_lyte = self.portInLyte.phi_lyte + self.c_lyte = self.portInLyte.c_lyte + self.phi_m = self.portInBulk.phi_m + + def get_trode_param(self, item): + """ + Shorthand to retrieve electrode-specific value + """ + value = self.config[self.trode, item] + # check if it is a particle-specific parameter + if item in self.config.params_per_particle: + value = value[self.ind] + return value + + 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 + self.noise = None + if self.get_trode_param("noise"): + numnoise = self.get_trode_param("numnoise") + noise_prefac = self.get_trode_param("noise_prefac") + tvec = np.linspace(0., 1.05*self.config["tend"], numnoise) + noise_data = noise_prefac*np.random.randn(numnoise, N) + self.noise = sintrp.interp1d(tvec, noise_data, axis=0, + 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, + self.config["elyteModelType"]) + + # Define average filling fraction in particle + eq = self.CreateEquation("cbar") + eq.Residual = self.cbar() + for k in range(N): + eq.Residual -= self.c(k) * volfrac_vec[k] + + # Define average rate of filling of particle + eq = self.CreateEquation("dcbardt") + eq.Residual = self.dcbardt() + for k in range(N): + eq.Residual -= self.c.dt(k) * volfrac_vec[k] + + c = np.empty(N, dtype=object) + 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) + 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) + + 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, + self.trode, self.ind) + eta = calc_eta(muR_surf, muO) + eta_eff = eta + self.Rxn()*self.get_trode_param("Rfilm") + if self.get_trode_param("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("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() + + 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) + dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), N) + + # Get solid particle chemical potential, overpotential, reaction rate + 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) + elif self.get_trode_param("type") in ["diffn", "CHR"]: + muR, actR = calc_muR(c, self.cbar(), self.config, self.trode, self.ind) + c_surf = c[-1] + muR_surf = muR[-1] + if actR is None: + actR_surf = None + else: + actR_surf = actR[-1] + eta = calc_eta(muR_surf, muO) + if self.get_trode_param("type") in ["ACR"]: + eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm") + for i in range(N)]) + else: + 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("lambda"), self.get_trode_param("alpha")) + if self.get_trode_param("type") in ["ACR"]: + for i in range(N): + eq = self.CreateEquation("Rxn_{i}".format(i=i)) + eq.Residual = self.Rxn(i) - Rxn[i] + else: + eq = self.CreateEquation("Rxn") + eq.Residual = self.Rxn() - Rxn + + # Get solid particle fluxes (if any) and RHS + if self.get_trode_param("type") in ["ACR"]: + RHS = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(N)]) + elif self.get_trode_param("type") in ["diffn", "CHR"]: + # Positive reaction (reduction, intercalation) is negative + # flux of Li at the surface. + Flux_bc = -self.Rxn() + Dfunc_name = self.get_trode_param("Dfunc") + Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), + Dfunc_name, + 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) + 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) + if self.get_trode_param("shape") == "sphere": + area_vec = 4*np.pi*edges**2 + elif self.get_trode_param("shape") == "cylinder": + area_vec = 2*np.pi*edges # per unit height + RHS = -np.diff(Flux_vec * area_vec) + + dcdt_vec = np.empty(N, dtype=object) + dcdt_vec[0:N] = [self.c.dt(k) for k in range(N)] + LHS_vec = MX(Mmat, dcdt_vec) + for k in range(N): + eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) + eq.Residual = LHS_vec[k] - RHS[k] + + +def calc_eta(muR, muO): + return muR - muO + + +def get_Mmat(shape, N): + r_vec, volfrac_vec = geo.get_unit_solid_discr(shape, N) + if shape == "C3": + Mmat = sprs.eye(N, N, format="csr") + elif shape in ["sphere", "cylinder"]: + Rs = 1. + # For discretization background, see Zeng & Bazant 2013 + # Mass matrix is common for each shape, diffn or CHR + if shape == "sphere": + Vp = 4./3. * np.pi * Rs**3 + elif shape == "cylinder": + Vp = np.pi * Rs**2 # per unit height + vol_vec = Vp * volfrac_vec + M1 = sprs.diags([1./8, 3./4, 1./8], [-1, 0, 1], + shape=(N, N), format="csr") + M1[1,0] = M1[-2,-1] = 1./4 + M2 = sprs.diags(vol_vec, 0, format="csr") + Mmat = M1*M2 + return Mmat + + +def calc_flux_diffn(c, D, Dfunc, E_D, Flux_bc, dr, T, noise): + N = len(c) + Flux_vec = np.empty(N+1, dtype=object) + Flux_vec[0] = 0 # Symmetry at r=0 + 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 + else: + Flux_vec[1:N] = -D * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ + np.diff(c + noise(dae.Time().Value))/dr + return Flux_vec + + +def calc_flux_CHR(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): + N = len(c) + Flux_vec = np.empty(N+1, dtype=object) + Flux_vec[0] = 0 # Symmetry at r=0 + Flux_vec[-1] = Flux_bc + c_edges = utils.mean_linear(c) + if noise is None: + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu)/dr + else: + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ + np.diff(mu + noise(dae.Time().Value))/dr + return Flux_vec + +def calc_flux_C3ver(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): + N = len(c) + Flux_vec = np.empty(N+1, dtype=object) + Flux_vec[0] = 0 # Symmetry at r=0 + Flux_vec[-1] = Flux_bc + c_edges = utils.mean_linear(c) + if noise is None: + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu)/dr + else: + Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ + np.diff(mu + noise(dae.Time().Value))/dr + return Flux_vec + + +def calc_flux_CHR2(c1, c2, mu1_R, mu2_R, D, Dfunc, E_D, Flux1_bc, Flux2_bc, dr, T, noise1, noise2): + N = len(c1) + Flux1_vec = np.empty(N+1, dtype=object) + Flux2_vec = np.empty(N+1, dtype=object) + Flux1_vec[0] = 0. # symmetry at r=0 + Flux2_vec[0] = 0. # symmetry at r=0 + Flux1_vec[-1] = Flux1_bc + Flux2_vec[-1] = Flux2_bc + c1_edges = utils.mean_linear(c1) + c2_edges = utils.mean_linear(c2) + if noise1 is None: + Flux1_vec[1:N] = -D/T * Dfunc(c1_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu1_R)/dr + Flux2_vec[1:N] = -D/T * Dfunc(c2_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu2_R)/dr + else: + Flux1_vec[1:N] = -D/T * Dfunc(c1_edges) * np.exp(-E_D/T + E_D/1) * \ + np.diff(mu1_R+noise1(dae.Time().Value))/dr + Flux2_vec[1:N] = -D/T * Dfunc(c2_edges) * np.exp(-E_D/T + E_D/1) * \ + np.diff(mu2_R+noise2(dae.Time().Value))/dr + return Flux1_vec, Flux2_vec + + +def calc_mu_O(c_lyte, phi_lyte, phi_sld, T, elyteModelType): + if elyteModelType == "SM": + mu_lyte = phi_lyte + act_lyte = c_lyte + elif elyteModelType == "dilute": + act_lyte = c_lyte + mu_lyte = T*np.log(act_lyte) + phi_lyte + mu_O = mu_lyte - phi_sld + return mu_O, act_lyte + + +def calc_muR(c, cbar, config, trode, ind): + muRfunc = props_am.muRfuncs(config, trode, ind).muRfunc + muR_ref = config[trode, "muR_ref"] + muR, actR = muRfunc(c, cbar, muR_ref) + return muR, actR + + +def MX(mat, objvec): + if not isinstance(mat, sprs.csr.csr_matrix): + raise Exception("MX function designed for csr mult") + n = objvec.shape[0] + if isinstance(objvec[0], dae.pyCore.adouble): + out = np.empty(n, dtype=object) + else: + out = np.zeros(n, dtype=float) + # Loop through the rows + for i in range(n): + low = mat.indptr[i] + up = mat.indptr[i+1] + if up > low: + out[i] = np.sum( + mat.data[low:up] * objvec[mat.indices[low:up]]) + else: + out[i] = 0.0 + return out diff --git a/mpet/sim.py b/mpet/sim.py index d5f16102..00d75ce9 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -53,12 +53,11 @@ def SetUpParametersAndDomains(self): self.m.DmnPart[tr].CreateArray(config["Npart"][tr]) for i in range(config["Nvol"][tr]): for j in range(config["Npart"][tr]): - # self.m.particles[tr][i, j].Dmn.CreateArray( - # int(config["psd_num"][tr][i,j])) self.m.particles[tr][i, j].Dmn.CreateArray( - 50) - self.m.particles[tr][i, j].Dmny.CreateArray( - 10) + int(config["psd_num"][tr][i,j])) + if self.config[tr, "type"] in 'ACR2D': + self.m.particles[tr][i, j].Dmny.CreateArray( + int(config["psd_num_ver"][tr][i,j])) def SetUpVariables(self): config = self.config @@ -90,12 +89,13 @@ def SetUpVariables(self): if solidType in constants.one_var_types: part.cbar.SetInitialGuess(cs0) if solidType == "ACR2D": - for k in range(50): + for k in range(Nij): part.cx.SetInitialCondition(k, cs0) for j in range(10): part.cy[k].SetInitialCondition(j, cs0) - for k in range(Nij): - part.c.SetInitialCondition(k, cs0) + else: + for k in range(Nij): + part.c.SetInitialCondition(k, cs0) elif solidType in constants.two_var_types: part.c1bar.SetInitialGuess(cs0) part.c2bar.SetInitialGuess(cs0) From 4fcce5fe2f1d2a388c9ecde69a08e9c697cdeb24 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 13:35:18 +0000 Subject: [PATCH 03/27] vertical discretization from config --- configs/params_LFP2D.cfg | 2 +- mpet/config/configuration.py | 2 ++ mpet/config/constants.py | 2 +- mpet/mod_electrodes.py | 2 +- mpet/sim.py | 3 ++- 5 files changed, 7 insertions(+), 4 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index e9f418ce..a4e4099e 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -5,7 +5,7 @@ type = ACR2D discretization = 1e-9 shape = C3 -discretization_ver = 3e-9 +discretization_ver = 1.5e-9 thickness = 20e-9 [Material] diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index cdd2b503..5ba4532d 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -605,6 +605,8 @@ def _distr_part(self): self['psd_area'] = {} self['psd_vol'] = {} self['psd_vol_FracVol'] = {} + # if self[trode,'type'] in 'ACR2D': + self['psd_num_ver'] = {} for trode in self['trodes']: solidType = self[trode, 'type'] diff --git a/mpet/config/constants.py b/mpet/config/constants.py index 9579ec5c..375b23b4 100644 --- a/mpet/config/constants.py +++ b/mpet/config/constants.py @@ -26,6 +26,6 @@ #: subset of ``PARAMS_PER_TRODE``` that is defined for the separator as well PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp'] #: parameters that are defined for each particle, and their type -PARAMS_PARTICLE = {'N': int, 'kappa': float, 'beta_s': float, 'D': float, 'k0': float, +PARAMS_PARTICLE = {'N': int, 'N_ver':int, 'kappa': float, 'beta_s': float, 'D': float, 'k0': float, 'Rfilm': float, 'delta_L': float, 'Omega_a': float, 'E_D': float, 'E_A': float} diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 45f9e39c..1ee5f1af 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -96,7 +96,7 @@ def DeclareEquations(self): dae.daeModel.DeclareEquations(self) # Nx = self.get_trode_param("N") # number of grid points in particle Nx = self.get_trode_param("N") - Ny = 10 + Ny = self.get_trode_param("N_ver") T = self.config["T"] # nondimensional temperature # r_vec, volfrac_vec_x = geo.get_unit_solid_discr(self.get_trode_param('shape'), Nx) # r_vec, volfrac_vec_y = geo.get_unit_solid_discr(self.get_trode_param('shape'), Ny) diff --git a/mpet/sim.py b/mpet/sim.py index 00d75ce9..44de3866 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -89,9 +89,10 @@ def SetUpVariables(self): if solidType in constants.one_var_types: part.cbar.SetInitialGuess(cs0) if solidType == "ACR2D": + N_ver_ij = config["psd_num_ver"][tr][i,j] for k in range(Nij): part.cx.SetInitialCondition(k, cs0) - for j in range(10): + for j in range(N_ver_ij): part.cy[k].SetInitialCondition(j, cs0) else: for k in range(Nij): From 003925b8d940e5a2a245aab53f26d8577f48a07c Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 14:10:07 +0000 Subject: [PATCH 04/27] fix ref D --- mpet/config/configuration.py | 9 +- mpet/mod_electrodes.py | 71 +--- mpet/mod_electrodes_old.py | 785 ----------------------------------- 3 files changed, 9 insertions(+), 856 deletions(-) delete mode 100644 mpet/mod_electrodes_old.py diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 5ba4532d..370c135e 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -607,6 +607,7 @@ def _distr_part(self): self['psd_vol_FracVol'] = {} # if self[trode,'type'] in 'ACR2D': self['psd_num_ver'] = {} + self['psd_len_ver'] = {} for trode in self['trodes']: solidType = self[trode, 'type'] @@ -652,6 +653,7 @@ def _distr_part(self): psd_num = np.ceil(raw/ solidDisc).astype(int) psd_len = solidDisc * psd_num psd_num_ver = np.ceil(raw_t/solidDisc_ver).astype(int) + psd_len_ver = solidDisc_ver * psd_num_ver elif solidType in ['CHR', 'diffn', 'CHR2', 'diffn2']: psd_num = np.ceil(raw / solidDisc).astype(int) + 1 psd_len = solidDisc * (psd_num - 1) @@ -692,6 +694,7 @@ def _distr_part(self): # store values to config 2d if self[trode,'type'] in 'ACR2D': self['psd_num_ver'][trode] = psd_num_ver + self['psd_len_ver'][trode] = psd_len_ver def _G(self): """ @@ -753,7 +756,11 @@ def _indvPart(self): nd_dgammadc = self[trode, 'dgammadc'] * cs_ref_part / gamma_S_ref self[trode, 'indvPart']['beta_s'][i, j] = nd_dgammadc \ / self[trode, 'indvPart']['kappa'][i, j] - self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] * self['t_ref'] / plen**2 + if self[trode,'type'] in 'ACR2D': + pthick = self['psd_len_ver'][trode][i,j] + self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] * self['t_ref'] / pthick**2 + else: + self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] * self['t_ref'] / plen**2 self[trode, 'indvPart']['E_D'][i, j] = self[trode, 'E_D'] \ / (constants.k * constants.N_A * constants.T_ref) self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \ diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 1ee5f1af..e6672854 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -50,18 +50,10 @@ def __init__(self, config, trode, vInd, pInd, "cy{k}".format(k=k), mole_frac_t, self, "Concentration in y direction of element {k}".format(k=k), [self.Dmny]) - - # self.cybar[k] = dae.daeVariable( - # "cybar{k}".format(k=k), mole_frac_t, self, - # "Average concentration in y of element {k}".format(k=k)) self.cbar = dae.daeVariable( "cbar", mole_frac_t, self, - "Average concentration in active particle") - # self.cxbar = dae.daeVariable( - # "cxbar", mole_frac_t, self, - # "Average concentration in x of active particle") - + "Average concentration in active particle") self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn]) @@ -94,26 +86,14 @@ def get_trode_param(self, item): def DeclareEquations(self): dae.daeModel.DeclareEquations(self) - # Nx = self.get_trode_param("N") # number of grid points in particle Nx = self.get_trode_param("N") Ny = self.get_trode_param("N_ver") T = self.config["T"] # nondimensional temperature - # r_vec, volfrac_vec_x = geo.get_unit_solid_discr(self.get_trode_param('shape'), Nx) - # r_vec, volfrac_vec_y = geo.get_unit_solid_discr(self.get_trode_param('shape'), Ny) self.noise = None mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T, "SM") - - # Define average filling fractions in particle - # eqx = self.CreateEquation("cxbar") - # for i in range(Nx): - # eqy = self.CreateEquation("cybar_{i}".format(i=i)) - # eqy.Residual = self.cx(i) - # for j in range(Ny): - # eqy.Residual -= self.cy[i](j)/Ny - for k in range(Nx): eq = self.CreateEquation("avgs_cy_cx{k}".format(k=k)) @@ -151,21 +131,8 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), Dfunc_name, f"mpet.electrode.diffusion.{Dfunc_name}") - # Equations for concentration evolution - # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors dr, edges = geo.get_dr_edges("C3", Ny) - # dr = 1 - # muR, actR = calc_muR(c_mat[:,0], self.cbar(), self.config, self.trode, self.ind) - # muR_mat = np.empty((Nx,Ny), dtype=object) - # actR_mat = np.empty((Nx,Ny), dtype=object) - # muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), - # self.config, self.trode, self.ind) - # for k in range(Nx): - # muR_mat[k,:], actR_mat[k,:] = calc_muR(c_mat[k,:], self.cx(k), - # self.config, self.trode, self.ind) - - # print(muR_mat[1,:]) c_surf = c_mat[-1,:] muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), @@ -183,57 +150,21 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq = self.CreateEquation("Rxn_{i}".format(i=i)) eq.Residual = self.Rxn(i) - Rxn[i] - # RHSx = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(Nx)]) - # dcdt_vec_x = np.empty(Nx, dtype=object) - # dcdt_vec_x[0:Nx] = [self.cx.dt(k) for k in range(Nx)] - # Mmatx = get_Mmat("C3", Nx) - # LHS_vec_x = MX(Mmatx, dcdt_vec_x) - - # for k in range(Nx): - # eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) - # eq.Residual = LHS_vec_x[k] - RHSx[k] - - # Flux_mat = np.empty((Nx,Ny), dtype=object) - # RHS_mat = np.empty((Ny,Nx), dtype=object) area_vec = 1 - # LHS_mat = np.empty((Ny,Nx), dtype=object) Mmaty = get_Mmat("C3", Ny) - # print(c_mat[1,:]) for k in range(Nx): - # print(k, ' --- \n') Flux_bc = -self.Rxn(k) - # eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) - # eq.Residual = LHS_vec_x[k] - self.get_trode_param("delta_L")*self.Rxn(k) - # muR_vec = muR_mat[k,:] - # c_v = c_mat[:,k] - # c_vec = utils.get_var_vec(c_mat[:,k],Ny) muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cbar(), self.config, self.trode, self.ind) - - # Flux_mat[k,:] = calc_flux_CHR(c_mat[k,:], muR_vec, self.get_trode_param("D"), Dfunc, - # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) Flux_vec = calc_flux_C3ver(c_mat[:,k], muR_vec, self.get_trode_param("D"), Dfunc, self.get_trode_param("E_D"), Flux_bc, dr, T, noise) - # Flux_vec = calc_flux_diffn(c_mat[:,k], self.get_trode_param("D"), Dfunc, - # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) - - # c_vec = utils.get_var_vec(c_mat[:,k],Ny) - - - # RHS_mat = -np.diff(Flux_vec * area_vec) RHS_vec = -np.diff(Flux_vec * area_vec) - dcdt_vec_y = np.empty(Ny, dtype=object) dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] - LHS_vec_y = MX(Mmaty, dcdt_vec_y) for j in range(Ny): eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) eq.Residual = LHS_vec_y[j] - RHS_vec[j] - # LHS_mat[:,k] = LHS_vec_y - - - class Mod2var(dae.daeModel): def __init__(self, config, trode, vInd, pInd, diff --git a/mpet/mod_electrodes_old.py b/mpet/mod_electrodes_old.py deleted file mode 100644 index 0a00c2f9..00000000 --- a/mpet/mod_electrodes_old.py +++ /dev/null @@ -1,785 +0,0 @@ -"""These models define individual particles of active material. - -This includes the equations for both 1-parameter models and 2-parameters models defining - - mass conservation (concentration evolution) - - reaction rate at the surface of the particles -In each model class it has options for different types of particles: - - homogeneous - - Fick-like diffusion - - Cahn-Hilliard (with reaction boundary condition) - - Allen-Cahn (with reaction throughout the particle) -These models can be instantiated from the mod_cell module to simulate various types of active -materials within a battery electrode. -""" -import daetools.pyDAE as dae -import numpy as np -import scipy.sparse as sprs -import scipy.interpolate as sintrp - -import mpet.geometry as geo -import mpet.ports as ports -import mpet.props_am as props_am -import mpet.utils as utils -from mpet.daeVariableTypes import mole_frac_t - - -class Mod2D(dae.daeModel): - def __init__(self, config, trode, vInd, pInd, - Name, Parent=None, Description=""): - super().__init__(Name, Parent, Description) - - self.config = config - self.trode = trode - self.ind = (vInd, pInd) - - # Domain - self.Dmn = dae.daeDomain("discretizationDomainX", self, dae.unit(), - "discretization domain in x direction") - self.Dmny = dae.daeDomain("discretizationDomainY", self, dae.unit(), - "discretization domain in y direction") - - # Variables - self.cx = dae.daeVariable( - "cx", mole_frac_t, self, - "Concentration in x direction of active particle", [self.Dmn]) - - Nx = self.get_trode_param("N") - self.cy = {} - for k in range(Nx): - self.cy[k] = dae.daeVariable( - "cy{k}".format(k=k), mole_frac_t, self, - "Concentration in y direction of element {k}".format(k=k), - [self.Dmny]) - - # self.cybar[k] = dae.daeVariable( - # "cybar{k}".format(k=k), mole_frac_t, self, - # "Average concentration in y of element {k}".format(k=k)) - - self.cbar = dae.daeVariable( - "cbar", mole_frac_t, self, - "Average concentration in active particle") - # self.cxbar = dae.daeVariable( - # "cxbar", mole_frac_t, self, - # "Average concentration in x of active particle") - - - self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") - self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn]) - - # Get reaction rate function - rxnType = config[trode, "rxnType"] - self.calc_rxn_rate = utils.import_function(config[trode, "rxnType_filename"], - rxnType, - f"mpet.electrode.reactions.{rxnType}") - - # Ports - self.portInLyte = ports.portFromElyte( - "portInLyte", dae.eInletPort, self, "Inlet port from electrolyte") - self.portInBulk = ports.portFromBulk( - "portInBulk", dae.eInletPort, self, - "Inlet port from e- conducting phase") - self.phi_lyte = self.portInLyte.phi_lyte - self.c_lyte = self.portInLyte.c_lyte - self.phi_m = self.portInBulk.phi_m - - def get_trode_param(self, item): - """ - Shorthand to retrieve electrode-specific value - """ - value = self.config[self.trode, item] - # check if it is a particle-specific parameter - if item in self.config.params_per_particle: - value = value[self.ind] - return value - - def DeclareEquations(self): - dae.daeModel.DeclareEquations(self) - # Nx = self.get_trode_param("N") # number of grid points in particle - Nx = self.get_trode_param("N") - Ny = 10 - T = self.config["T"] # nondimensional temperature - # r_vec, volfrac_vec_x = geo.get_unit_solid_discr(self.get_trode_param('shape'), Nx) - # r_vec, volfrac_vec_y = geo.get_unit_solid_discr(self.get_trode_param('shape'), Ny) - - self.noise = None - - mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T, - "SM") - - # Define average filling fractions in particle - # eqx = self.CreateEquation("cxbar") - # for i in range(Nx): - # eqy = self.CreateEquation("cybar_{i}".format(i=i)) - # eqy.Residual = self.cx(i) - # for j in range(Ny): - # eqy.Residual -= self.cy[i](j)/Ny - - - for k in range(Nx): - eq = self.CreateEquation("avgs_cy_cx{k}".format(k=k)) - eq.Residual = self.cx(k) - for j in range(Ny): - eq.Residual -= self.cy[k](j)/Ny - - eq = self.CreateEquation("cbar") - eq.Residual = self.cbar() - for k in range(Nx): - eq.Residual -= self.cx(k)/Nx - - # Define average rate of filling of particle - eq = self.CreateEquation("dcbardt") - eq.Residual = self.dcbardt() - for k in range(Nx): - eq.Residual -= self.cx.dt(k)/Nx - - # cx = np.empty(Nx, dtype=object) - # cx[:] = [self.cx(k) for k in range(Nx)] - c_mat = np.empty((Ny, Nx), dtype=object) - for k in range(Nx): - c_mat[:,k] = [self.cy[k](j) for j in range(Ny)] - - self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) - - for eq in self.Equations: - eq.CheckUnitsConsistency = False - - def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): - Ny = np.size(c_mat, 0) - Nx = np.size(c_mat, 1) - T = self.config["T"] - Dfunc_name = self.get_trode_param("Dfunc") - Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), - Dfunc_name, - f"mpet.electrode.diffusion.{Dfunc_name}") - # Equations for concentration evolution - # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors - - dr, edges = geo.get_dr_edges("C3", Ny) - # dr = 1 - # muR, actR = calc_muR(c_mat[:,0], self.cbar(), self.config, self.trode, self.ind) - # muR_mat = np.empty((Nx,Ny), dtype=object) - # actR_mat = np.empty((Nx,Ny), dtype=object) - # muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), - # self.config, self.trode, self.ind) - # for k in range(Nx): - # muR_mat[k,:], actR_mat[k,:] = calc_muR(c_mat[k,:], self.cx(k), - # self.config, self.trode, self.ind) - - # print(muR_mat[1,:]) - - c_surf = c_mat[-1,:] - muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), - self.config, self.trode, self.ind) - - eta = calc_eta(muR_surf, muO) - eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm") - for i in range(Nx)]) - 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("lambda"), self.get_trode_param("alpha")) - - for i in range(Nx): - eq = self.CreateEquation("Rxn_{i}".format(i=i)) - eq.Residual = self.Rxn(i) - Rxn[i] - - RHSx = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(Nx)]) - dcdt_vec_x = np.empty(Nx, dtype=object) - dcdt_vec_x[0:Nx] = [self.cx.dt(k) for k in range(Nx)] - Mmatx = get_Mmat("C3", Nx) - LHS_vec_x = MX(Mmatx, dcdt_vec_x) - - for k in range(Nx): - eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) - eq.Residual = LHS_vec_x[k] - RHSx[k] - - # Flux_mat = np.empty((Nx,Ny), dtype=object) - # RHS_mat = np.empty((Ny,Nx), dtype=object) - area_vec = 1 - # LHS_mat = np.empty((Ny,Nx), dtype=object) - Mmaty = get_Mmat("C3", Ny) - # print(c_mat[1,:]) - for k in range(Nx): - # print(k, ' --- \n') - Flux_bc = -self.Rxn(k) - # muR_vec = muR_mat[k,:] - # c_v = c_mat[:,k] - # c_vec = utils.get_var_vec(c_mat[:,k],Ny) - muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cbar(), - self.config, self.trode, self.ind) - - # Flux_mat[k,:] = calc_flux_CHR(c_mat[k,:], muR_vec, self.get_trode_param("D"), Dfunc, - # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) - Flux_vec = calc_flux_C3ver(c_mat[:,k], muR_vec, self.get_trode_param("D"), Dfunc, - self.get_trode_param("E_D"), Flux_bc, dr, T, noise) - # Flux_vec = calc_flux_diffn(c_mat[:,k], self.get_trode_param("D"), Dfunc, - # self.get_trode_param("E_D"), Flux_bc, dr, T, noise) - - # c_vec = utils.get_var_vec(c_mat[:,k],Ny) - - - # RHS_mat = -np.diff(Flux_vec * area_vec) - RHS_vec = -np.diff(Flux_vec * area_vec) - - dcdt_vec_y = np.empty(Ny, dtype=object) - dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] - - LHS_vec_y = MX(Mmaty, dcdt_vec_y) - for j in range(Ny): - eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) - eq.Residual = LHS_vec_y[j] - RHS_vec[j] - # LHS_mat[:,k] = LHS_vec_y - - - - -class Mod2var(dae.daeModel): - def __init__(self, config, trode, vInd, pInd, - Name, Parent=None, Description=""): - super().__init__(Name, Parent, Description) - - self.config = config - self.trode = trode - self.ind = (vInd, pInd) - - # Domain - self.Dmn = dae.daeDomain("discretizationDomain", self, dae.unit(), - "discretization domain") - - # Variables - self.c1 = dae.daeVariable( - "c1", mole_frac_t, self, - "Concentration in 'layer' 1 of active particle", [self.Dmn]) - self.c2 = dae.daeVariable( - "c2", mole_frac_t, self, - "Concentration in 'layer' 2 of active particle", [self.Dmn]) - self.cbar = dae.daeVariable( - "cbar", mole_frac_t, self, - "Average concentration in active particle") - self.c1bar = dae.daeVariable( - "c1bar", mole_frac_t, self, - "Average concentration in 'layer' 1 of active particle") - self.c2bar = dae.daeVariable( - "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") - 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") - else: - self.Rxn1 = dae.daeVariable("Rxn1", dae.no_t, self, "Rate of reaction 1", [self.Dmn]) - self.Rxn2 = dae.daeVariable("Rxn2", dae.no_t, self, "Rate of reaction 2", [self.Dmn]) - - # Get reaction rate function - rxnType = config[trode, "rxnType"] - self.calc_rxn_rate = utils.import_function(config[trode, "rxnType_filename"], - rxnType, - f"mpet.electrode.reactions.{rxnType}") - - # Ports - self.portInLyte = ports.portFromElyte( - "portInLyte", dae.eInletPort, self, "Inlet port from electrolyte") - self.portInBulk = ports.portFromBulk( - "portInBulk", dae.eInletPort, self, - "Inlet port from e- conducting phase") - self.phi_lyte = self.portInLyte.phi_lyte - self.c_lyte = self.portInLyte.c_lyte - self.phi_m = self.portInBulk.phi_m - - def get_trode_param(self, item): - """ - Shorthand to retrieve electrode-specific value - """ - value = self.config[self.trode, item] - # check if it is a particle-specific parameter - if item in self.config.params_per_particle: - value = value[self.ind] - return value - - 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 - self.noise1 = self.noise2 = None - if self.get_trode_param("noise"): - numnoise = self.get_trode_param("numnoise") - noise_prefac = self.get_trode_param("noise_prefac") - tvec = np.linspace(0., 1.05*self.config["tend"], numnoise) - noise_data1 = noise_prefac*np.random.randn(numnoise, N) - noise_data2 = noise_prefac*np.random.randn(numnoise, N) - self.noise1 = sintrp.interp1d(tvec, noise_data1, axis=0, - bounds_error=False, fill_value=0.) - self.noise2 = sintrp.interp1d(tvec, noise_data2, axis=0, - bounds_error=False, fill_value=0.) - noises = (self.noise1, self.noise2) - - # 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.config["elyteModelType"]) - - # Define average filling fractions in particle - eq1 = self.CreateEquation("c1bar") - eq2 = self.CreateEquation("c2bar") - eq1.Residual = self.c1bar() - eq2.Residual = self.c2bar() - for k in range(N): - eq1.Residual -= self.c1(k) * volfrac_vec[k] - eq2.Residual -= self.c2(k) * volfrac_vec[k] - eq = self.CreateEquation("cbar") - eq.Residual = self.cbar() - .5*(self.c1bar() + self.c2bar()) - - # Define average rate of filling of particle - eq = self.CreateEquation("dcbardt") - eq.Residual = self.dcbardt() - for k in range(N): - eq.Residual -= .5*(self.c1.dt(k) + 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) - 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) - - 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, - self.trode, self.ind) - eta1 = calc_eta(mu1R_surf, muO) - eta2 = calc_eta(mu2R_surf, muO) - eta1_eff = eta1 + self.Rxn1()*self.get_trode_param("Rfilm") - eta2_eff = eta2 + self.Rxn2()*self.get_trode_param("Rfilm") - noise1, noise2 = noises - if self.get_trode_param("noise"): - eta1_eff += noise1(dae.Time().Value) - 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("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("lambda"), self.get_trode_param("alpha")) - eq1 = self.CreateEquation("Rxn1") - eq2 = self.CreateEquation("Rxn2") - eq1.Residual = self.Rxn1() - Rxn1[0] - eq2.Residual = self.Rxn2() - Rxn2[0] - - eq1 = self.CreateEquation("dc1sdt") - 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] - - 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) - dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), N) - - # 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) - c1_surf = c1[-1] - c2_surf = c2[-1] - mu1R_surf, act1R_surf = mu1R[-1], act1R[-1] - mu2R_surf, act2R_surf = mu2R[-1], act2R[-1] - eta1 = calc_eta(mu1R_surf, muO) - eta2 = calc_eta(mu2R_surf, muO) - if self.get_trode_param("type") in ["ACR2"]: - eta1_eff = np.array([eta1[i] - + self.Rxn1(i)*self.get_trode_param("Rfilm") for i in range(N)]) - eta2_eff = np.array([eta2[i] - + self.Rxn2(i)*self.get_trode_param("Rfilm") for i in range(N)]) - else: - eta1_eff = eta1 + self.Rxn1()*self.get_trode_param("Rfilm") - 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("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("lambda"), self.get_trode_param("alpha")) - if self.get_trode_param("type") in ["ACR2"]: - for i in range(N): - eq1 = self.CreateEquation("Rxn1_{i}".format(i=i)) - eq2 = self.CreateEquation("Rxn2_{i}".format(i=i)) - eq1.Residual = self.Rxn1(i) - Rxn1[i] - eq2.Residual = self.Rxn2(i) - Rxn2[i] - else: - eq1 = self.CreateEquation("Rxn1") - eq2 = self.CreateEquation("Rxn2") - eq1.Residual = self.Rxn1() - Rxn1 - eq2.Residual = self.Rxn2() - Rxn2 - - # Get solid particle fluxes (if any) and RHS - if self.get_trode_param("type") in ["diffn2", "CHR2"]: - # Positive reaction (reduction, intercalation) is negative - # flux of Li at the surface. - Flux1_bc = -0.5 * self.Rxn1() - Flux2_bc = -0.5 * self.Rxn2() - Dfunc_name = self.get_trode_param("Dfunc") - Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), - Dfunc_name, - f"mpet.electrode.diffusion.{Dfunc_name}") - if self.get_trode_param("type") == "CHR2": - 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) - if self.get_trode_param("shape") == "sphere": - area_vec = 4*np.pi*edges**2 - elif self.get_trode_param("shape") == "cylinder": - area_vec = 2*np.pi*edges # per unit height - RHS1 = -np.diff(Flux1_vec * area_vec) - RHS2 = -np.diff(Flux2_vec * area_vec) -# kinterlayer = 1e-3 -# interLayerRxn = (kinterlayer * (1 - c1) * (1 - c2) * (act1R - act2R)) -# RxnTerm1 = -interLayerRxn -# RxnTerm2 = interLayerRxn - RxnTerm1 = 0 - RxnTerm2 = 0 - RHS1 += RxnTerm1 - RHS2 += RxnTerm2 - - dc1dt_vec = np.empty(N, dtype=object) - dc2dt_vec = np.empty(N, dtype=object) - dc1dt_vec[0:N] = [self.c1.dt(k) for k in range(N)] - dc2dt_vec[0:N] = [self.c2.dt(k) for k in range(N)] - LHS1_vec = MX(Mmat, dc1dt_vec) - LHS2_vec = MX(Mmat, dc2dt_vec) - for k in range(N): - eq1 = self.CreateEquation("dc1sdt_discr{k}".format(k=k)) - eq2 = self.CreateEquation("dc2sdt_discr{k}".format(k=k)) - eq1.Residual = LHS1_vec[k] - RHS1[k] - eq2.Residual = LHS2_vec[k] - RHS2[k] - - -class Mod1var(dae.daeModel): - def __init__(self, config, trode, vInd, pInd, - Name, Parent=None, Description=""): - super().__init__(Name, Parent, Description) - - self.config = config - self.trode = trode - self.ind = (vInd, pInd) - - # Domain - self.Dmn = dae.daeDomain("discretizationDomain", self, dae.unit(), - "discretization domain") - - # Variables - self.c = dae.daeVariable("c", mole_frac_t, self, - "Concentration in active particle", - [self.Dmn]) - self.cbar = dae.daeVariable( - "cbar", mole_frac_t, self, - "Average concentration in active particle") - self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") - if config[trode, "type"] not in ["ACR"]: - self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction") - else: - self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn]) - - # Get reaction rate function - rxnType = config[trode, "rxnType"] - self.calc_rxn_rate = utils.import_function(config[trode, "rxnType_filename"], - rxnType, - f"mpet.electrode.reactions.{rxnType}") - - # Ports - self.portInLyte = ports.portFromElyte( - "portInLyte", dae.eInletPort, self, - "Inlet port from electrolyte") - self.portInBulk = ports.portFromBulk( - "portInBulk", dae.eInletPort, self, - "Inlet port from e- conducting phase") - self.phi_lyte = self.portInLyte.phi_lyte - self.c_lyte = self.portInLyte.c_lyte - self.phi_m = self.portInBulk.phi_m - - def get_trode_param(self, item): - """ - Shorthand to retrieve electrode-specific value - """ - value = self.config[self.trode, item] - # check if it is a particle-specific parameter - if item in self.config.params_per_particle: - value = value[self.ind] - return value - - 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 - self.noise = None - if self.get_trode_param("noise"): - numnoise = self.get_trode_param("numnoise") - noise_prefac = self.get_trode_param("noise_prefac") - tvec = np.linspace(0., 1.05*self.config["tend"], numnoise) - noise_data = noise_prefac*np.random.randn(numnoise, N) - self.noise = sintrp.interp1d(tvec, noise_data, axis=0, - 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, - self.config["elyteModelType"]) - - # Define average filling fraction in particle - eq = self.CreateEquation("cbar") - eq.Residual = self.cbar() - for k in range(N): - eq.Residual -= self.c(k) * volfrac_vec[k] - - # Define average rate of filling of particle - eq = self.CreateEquation("dcbardt") - eq.Residual = self.dcbardt() - for k in range(N): - eq.Residual -= self.c.dt(k) * volfrac_vec[k] - - c = np.empty(N, dtype=object) - 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) - 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) - - 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, - self.trode, self.ind) - eta = calc_eta(muR_surf, muO) - eta_eff = eta + self.Rxn()*self.get_trode_param("Rfilm") - if self.get_trode_param("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("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() - - 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) - dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), N) - - # Get solid particle chemical potential, overpotential, reaction rate - 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) - elif self.get_trode_param("type") in ["diffn", "CHR"]: - muR, actR = calc_muR(c, self.cbar(), self.config, self.trode, self.ind) - c_surf = c[-1] - muR_surf = muR[-1] - if actR is None: - actR_surf = None - else: - actR_surf = actR[-1] - eta = calc_eta(muR_surf, muO) - if self.get_trode_param("type") in ["ACR"]: - eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm") - for i in range(N)]) - else: - 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("lambda"), self.get_trode_param("alpha")) - if self.get_trode_param("type") in ["ACR"]: - for i in range(N): - eq = self.CreateEquation("Rxn_{i}".format(i=i)) - eq.Residual = self.Rxn(i) - Rxn[i] - else: - eq = self.CreateEquation("Rxn") - eq.Residual = self.Rxn() - Rxn - - # Get solid particle fluxes (if any) and RHS - if self.get_trode_param("type") in ["ACR"]: - RHS = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(N)]) - elif self.get_trode_param("type") in ["diffn", "CHR"]: - # Positive reaction (reduction, intercalation) is negative - # flux of Li at the surface. - Flux_bc = -self.Rxn() - Dfunc_name = self.get_trode_param("Dfunc") - Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), - Dfunc_name, - 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) - 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) - if self.get_trode_param("shape") == "sphere": - area_vec = 4*np.pi*edges**2 - elif self.get_trode_param("shape") == "cylinder": - area_vec = 2*np.pi*edges # per unit height - RHS = -np.diff(Flux_vec * area_vec) - - dcdt_vec = np.empty(N, dtype=object) - dcdt_vec[0:N] = [self.c.dt(k) for k in range(N)] - LHS_vec = MX(Mmat, dcdt_vec) - for k in range(N): - eq = self.CreateEquation("dcsdt_discr{k}".format(k=k)) - eq.Residual = LHS_vec[k] - RHS[k] - - -def calc_eta(muR, muO): - return muR - muO - - -def get_Mmat(shape, N): - r_vec, volfrac_vec = geo.get_unit_solid_discr(shape, N) - if shape == "C3": - Mmat = sprs.eye(N, N, format="csr") - elif shape in ["sphere", "cylinder"]: - Rs = 1. - # For discretization background, see Zeng & Bazant 2013 - # Mass matrix is common for each shape, diffn or CHR - if shape == "sphere": - Vp = 4./3. * np.pi * Rs**3 - elif shape == "cylinder": - Vp = np.pi * Rs**2 # per unit height - vol_vec = Vp * volfrac_vec - M1 = sprs.diags([1./8, 3./4, 1./8], [-1, 0, 1], - shape=(N, N), format="csr") - M1[1,0] = M1[-2,-1] = 1./4 - M2 = sprs.diags(vol_vec, 0, format="csr") - Mmat = M1*M2 - return Mmat - - -def calc_flux_diffn(c, D, Dfunc, E_D, Flux_bc, dr, T, noise): - N = len(c) - Flux_vec = np.empty(N+1, dtype=object) - Flux_vec[0] = 0 # Symmetry at r=0 - 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 - else: - Flux_vec[1:N] = -D * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ - np.diff(c + noise(dae.Time().Value))/dr - return Flux_vec - - -def calc_flux_CHR(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): - N = len(c) - Flux_vec = np.empty(N+1, dtype=object) - Flux_vec[0] = 0 # Symmetry at r=0 - Flux_vec[-1] = Flux_bc - c_edges = utils.mean_linear(c) - if noise is None: - Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu)/dr - else: - Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ - np.diff(mu + noise(dae.Time().Value))/dr - return Flux_vec - -def calc_flux_C3ver(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): - N = len(c) - Flux_vec = np.empty(N+1, dtype=object) - Flux_vec[0] = 0 # Symmetry at r=0 - Flux_vec[-1] = Flux_bc - c_edges = utils.mean_linear(c) - if noise is None: - Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu)/dr - else: - Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.exp(-E_D/T + E_D/1) * \ - np.diff(mu + noise(dae.Time().Value))/dr - return Flux_vec - - -def calc_flux_CHR2(c1, c2, mu1_R, mu2_R, D, Dfunc, E_D, Flux1_bc, Flux2_bc, dr, T, noise1, noise2): - N = len(c1) - Flux1_vec = np.empty(N+1, dtype=object) - Flux2_vec = np.empty(N+1, dtype=object) - Flux1_vec[0] = 0. # symmetry at r=0 - Flux2_vec[0] = 0. # symmetry at r=0 - Flux1_vec[-1] = Flux1_bc - Flux2_vec[-1] = Flux2_bc - c1_edges = utils.mean_linear(c1) - c2_edges = utils.mean_linear(c2) - if noise1 is None: - Flux1_vec[1:N] = -D/T * Dfunc(c1_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu1_R)/dr - Flux2_vec[1:N] = -D/T * Dfunc(c2_edges) * np.exp(-E_D/T + E_D/1) * np.diff(mu2_R)/dr - else: - Flux1_vec[1:N] = -D/T * Dfunc(c1_edges) * np.exp(-E_D/T + E_D/1) * \ - np.diff(mu1_R+noise1(dae.Time().Value))/dr - Flux2_vec[1:N] = -D/T * Dfunc(c2_edges) * np.exp(-E_D/T + E_D/1) * \ - np.diff(mu2_R+noise2(dae.Time().Value))/dr - return Flux1_vec, Flux2_vec - - -def calc_mu_O(c_lyte, phi_lyte, phi_sld, T, elyteModelType): - if elyteModelType == "SM": - mu_lyte = phi_lyte - act_lyte = c_lyte - elif elyteModelType == "dilute": - act_lyte = c_lyte - mu_lyte = T*np.log(act_lyte) + phi_lyte - mu_O = mu_lyte - phi_sld - return mu_O, act_lyte - - -def calc_muR(c, cbar, config, trode, ind): - muRfunc = props_am.muRfuncs(config, trode, ind).muRfunc - muR_ref = config[trode, "muR_ref"] - muR, actR = muRfunc(c, cbar, muR_ref) - return muR, actR - - -def MX(mat, objvec): - if not isinstance(mat, sprs.csr.csr_matrix): - raise Exception("MX function designed for csr mult") - n = objvec.shape[0] - if isinstance(objvec[0], dae.pyCore.adouble): - out = np.empty(n, dtype=object) - else: - out = np.zeros(n, dtype=float) - # Loop through the rows - for i in range(n): - low = mat.indptr[i] - up = mat.indptr[i+1] - if up > low: - out[i] = np.sum( - mat.data[low:up] * objvec[mat.indices[low:up]]) - else: - out[i] = 0.0 - return out From 89f826a6281896c60d28737a37d1ac96186a3121 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 14:14:13 +0000 Subject: [PATCH 05/27] normal LFP works again --- mpet/config/configuration.py | 6 +++--- mpet/sim.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 370c135e..d3731206 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -692,7 +692,7 @@ def _distr_part(self): self['psd_vol'][trode] = psd_vol self['psd_vol_FracVol'][trode] = psd_frac_vol # store values to config 2d - if self[trode,'type'] in 'ACR2D': + if self[trode,'type'] == 'ACR2D': self['psd_num_ver'][trode] = psd_num_ver self['psd_len_ver'][trode] = psd_len_ver @@ -748,7 +748,7 @@ def _indvPart(self): kappa_ref = constants.k * constants.T_ref * cs_ref_part * plen**2 # J/m gamma_S_ref = kappa_ref / plen # J/m^2 # non-dimensional quantities - if self[trode,'type'] in 'ACR2D': + if self[trode,'type'] == 'ACR2D': self[trode, 'indvPart']['N_ver'][i,j] = self['psd_num_ver'][trode][i,j] if self[trode, 'kappa'] is not None: self[trode, 'indvPart']['kappa'][i, j] = self[trode, 'kappa'] / kappa_ref @@ -756,7 +756,7 @@ def _indvPart(self): nd_dgammadc = self[trode, 'dgammadc'] * cs_ref_part / gamma_S_ref self[trode, 'indvPart']['beta_s'][i, j] = nd_dgammadc \ / self[trode, 'indvPart']['kappa'][i, j] - if self[trode,'type'] in 'ACR2D': + if self[trode,'type'] == 'ACR2D': pthick = self['psd_len_ver'][trode][i,j] self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] * self['t_ref'] / pthick**2 else: diff --git a/mpet/sim.py b/mpet/sim.py index 44de3866..ee0d2544 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -55,7 +55,7 @@ def SetUpParametersAndDomains(self): for j in range(config["Npart"][tr]): self.m.particles[tr][i, j].Dmn.CreateArray( int(config["psd_num"][tr][i,j])) - if self.config[tr, "type"] in 'ACR2D': + if self.config[tr, "type"] == 'ACR2D': self.m.particles[tr][i, j].Dmny.CreateArray( int(config["psd_num_ver"][tr][i,j])) From 38a6c388b7e4dde496daf98c5263693ebdc7741b Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 14:17:47 +0000 Subject: [PATCH 06/27] style --- mpet/mod_electrodes.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index e6672854..2d6ea78d 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -32,7 +32,7 @@ def __init__(self, config, trode, vInd, pInd, self.trode = trode self.ind = (vInd, pInd) - # Domain + # Domain self.Dmn = dae.daeDomain("discretizationDomainX", self, dae.unit(), "discretization domain in x direction") self.Dmny = dae.daeDomain("discretizationDomainY", self, dae.unit(), @@ -111,9 +111,7 @@ def DeclareEquations(self): eq.Residual = self.dcbardt() for k in range(Nx): eq.Residual -= self.cx.dt(k)/Nx - - # cx = np.empty(Nx, dtype=object) - # cx[:] = [self.cx(k) for k in range(Nx)] + c_mat = np.empty((Ny, Nx), dtype=object) for k in range(Nx): c_mat[:,k] = [self.cy[k](j) for j in range(Ny)] @@ -161,8 +159,8 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): RHS_vec = -np.diff(Flux_vec * area_vec) dcdt_vec_y = np.empty(Ny, dtype=object) dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] - LHS_vec_y = MX(Mmaty, dcdt_vec_y) - for j in range(Ny): + LHS_vec_y = MX(Mmaty, dcdt_vec_y) + for j in range(Ny): eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) eq.Residual = LHS_vec_y[j] - RHS_vec[j] From 25e9f94e062fb3ed56f523487c3a0d4cb03b3820 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 15:29:36 +0100 Subject: [PATCH 07/27] style check mod_electrodes --- mpet/mod_electrodes.py | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 2d6ea78d..ad845dc0 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -36,7 +36,7 @@ def __init__(self, config, trode, vInd, pInd, self.Dmn = dae.daeDomain("discretizationDomainX", self, dae.unit(), "discretization domain in x direction") self.Dmny = dae.daeDomain("discretizationDomainY", self, dae.unit(), - "discretization domain in y direction") + "discretization domain in y direction") # Variables self.cx = dae.daeVariable( @@ -46,15 +46,14 @@ def __init__(self, config, trode, vInd, pInd, Nx = self.get_trode_param("N") self.cy = {} for k in range(Nx): - self.cy[k] = dae.daeVariable( - "cy{k}".format(k=k), mole_frac_t, self, - "Concentration in y direction of element {k}".format(k=k), - [self.Dmny]) - + self.cy[k] = dae.daeVariable("cy{k}".format(k=k), mole_frac_t, self, + "Concentration in y direction of element {k}".format(k=k), + [self.Dmny]) + self.cbar = dae.daeVariable( "cbar", mole_frac_t, self, - "Average concentration in active particle") - + "Average concentration in active particle") + self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling") self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn]) @@ -94,10 +93,10 @@ def DeclareEquations(self): mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T, "SM") - + for k in range(Nx): eq = self.CreateEquation("avgs_cy_cx{k}".format(k=k)) - eq.Residual = self.cx(k) + eq.Residual = self.cx(k) for j in range(Ny): eq.Residual -= self.cy[k](j)/Ny @@ -127,18 +126,18 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): T = self.config["T"] Dfunc_name = self.get_trode_param("Dfunc") Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), - Dfunc_name, - f"mpet.electrode.diffusion.{Dfunc_name}") + Dfunc_name, + f"mpet.electrode.diffusion.{Dfunc_name}") dr, edges = geo.get_dr_edges("C3", Ny) c_surf = c_mat[-1,:] - muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), - self.config, self.trode, self.ind) + muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), + self.config, self.trode, self.ind) eta = calc_eta(muR_surf, muO) eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm") - for i in range(Nx)]) + for i in range(Nx)]) 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, @@ -152,10 +151,10 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): Mmaty = get_Mmat("C3", Ny) for k in range(Nx): Flux_bc = -self.Rxn(k) - muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cbar(), - self.config, self.trode, self.ind) + muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cbar(), + self.config, self.trode, self.ind) Flux_vec = calc_flux_C3ver(c_mat[:,k], muR_vec, 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, T, noise) RHS_vec = -np.diff(Flux_vec * area_vec) dcdt_vec_y = np.empty(Ny, dtype=object) dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] @@ -164,6 +163,7 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) eq.Residual = LHS_vec_y[j] - RHS_vec[j] + class Mod2var(dae.daeModel): def __init__(self, config, trode, vInd, pInd, Name, Parent=None, Description=""): @@ -643,6 +643,7 @@ def calc_flux_CHR(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): np.diff(mu + noise(dae.Time().Value))/dr return Flux_vec + def calc_flux_C3ver(c, mu, D, Dfunc, E_D, Flux_bc, dr, T, noise): N = len(c) Flux_vec = np.empty(N+1, dtype=object) From ae88a1ef2504de1fea7757a0dcb856ac170548b7 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 15:32:28 +0100 Subject: [PATCH 08/27] configuration style check --- mpet/config/configuration.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index d3731206..087d9c70 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -650,7 +650,7 @@ def _distr_part(self): psd_len = solidDisc * psd_num elif solidType in ['ACR2D']: solidDisc_ver = self[trode, 'discretization_ver'] - psd_num = np.ceil(raw/ solidDisc).astype(int) + psd_num = np.ceil(raw / solidDisc).astype(int) psd_len = solidDisc * psd_num psd_num_ver = np.ceil(raw_t/solidDisc_ver).astype(int) psd_len_ver = solidDisc_ver * psd_num_ver @@ -757,10 +757,12 @@ def _indvPart(self): self[trode, 'indvPart']['beta_s'][i, j] = nd_dgammadc \ / self[trode, 'indvPart']['kappa'][i, j] if self[trode,'type'] == 'ACR2D': - pthick = self['psd_len_ver'][trode][i,j] - self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] * self['t_ref'] / pthick**2 + pthick = self['psd_len_ver'][trode][i,j] + self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] \ + * self['t_ref'] / pthick**2 else: - self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] * self['t_ref'] / plen**2 + self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] \ + * self['t_ref'] / plen**2 self[trode, 'indvPart']['E_D'][i, j] = self[trode, 'E_D'] \ / (constants.k * constants.N_A * constants.T_ref) self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \ From 6444910499eafd2c9da6752282353f8886d8cfe9 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 1 Dec 2022 16:27:11 +0100 Subject: [PATCH 09/27] nx now works only if it is 10 nm and 1 discr --- mpet/config/configuration.py | 2 +- mpet/mod_electrodes.py | 18 +++++++++--------- mpet/sim.py | 2 +- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 087d9c70..1da05a7b 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -621,7 +621,7 @@ def _distr_part(self): if self[trode,'type'] == 'ACR2D': np.random.seed(self.D_s['seed']) mean_t = self[trode,'thickness'] - stddev_t = mean_t*0.1 + stddev_t = mean_t*0 var_t = stddev_t**2 mu_t = np.log((mean_t**2) / np.sqrt(var_t + mean_t**2)) sigma_t = np.sqrt(np.log(var_t/(mean_t**2) + 1)) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index ad845dc0..7e841566 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -39,8 +39,8 @@ def __init__(self, config, trode, vInd, pInd, "discretization domain in y direction") # Variables - self.cx = dae.daeVariable( - "cx", mole_frac_t, self, + self.c = dae.daeVariable( + "c", mole_frac_t, self, "Concentration in x direction of active particle", [self.Dmn]) Nx = self.get_trode_param("N") @@ -96,20 +96,20 @@ def DeclareEquations(self): for k in range(Nx): eq = self.CreateEquation("avgs_cy_cx{k}".format(k=k)) - eq.Residual = self.cx(k) + eq.Residual = self.c(k) for j in range(Ny): eq.Residual -= self.cy[k](j)/Ny eq = self.CreateEquation("cbar") eq.Residual = self.cbar() for k in range(Nx): - eq.Residual -= self.cx(k)/Nx + eq.Residual -= self.c(k)/Nx # Define average rate of filling of particle eq = self.CreateEquation("dcbardt") eq.Residual = self.dcbardt() for k in range(Nx): - eq.Residual -= self.cx.dt(k)/Nx + eq.Residual -= self.c.dt(k)/Nx c_mat = np.empty((Ny, Nx), dtype=object) for k in range(Nx): @@ -143,11 +143,11 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): self.get_trode_param("E_A"), T, actR_surf, act_lyte, self.get_trode_param("lambda"), self.get_trode_param("alpha")) - for i in range(Nx): - eq = self.CreateEquation("Rxn_{i}".format(i=i)) - eq.Residual = self.Rxn(i) - Rxn[i] + for k in range(Nx): + eq = self.CreateEquation("Rxn_{k}".format(k=k)) + eq.Residual = self.Rxn(k) - Rxn[k] - area_vec = 1 + area_vec = 1. Mmaty = get_Mmat("C3", Ny) for k in range(Nx): Flux_bc = -self.Rxn(k) diff --git a/mpet/sim.py b/mpet/sim.py index ee0d2544..dbb41756 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -91,7 +91,7 @@ def SetUpVariables(self): if solidType == "ACR2D": N_ver_ij = config["psd_num_ver"][tr][i,j] for k in range(Nij): - part.cx.SetInitialCondition(k, cs0) + part.c.SetInitialCondition(k, cs0) for j in range(N_ver_ij): part.cy[k].SetInitialCondition(j, cs0) else: From 7fb18f937baf2d56e32289b0a45e737682a3b5a5 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Fri, 2 Dec 2022 13:05:12 +0000 Subject: [PATCH 10/27] solve error init and start fully 2D --- mpet/electrode/materials/LiFePO42D.py | 32 +++++++++ mpet/mod_electrodes.py | 95 +++++++++++++++++++++------ mpet/sim.py | 2 +- 3 files changed, 109 insertions(+), 20 deletions(-) create mode 100644 mpet/electrode/materials/LiFePO42D.py diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py new file mode 100644 index 00000000..e348af30 --- /dev/null +++ b/mpet/electrode/materials/LiFePO42D.py @@ -0,0 +1,32 @@ +import numpy as np + + +def LiFePO42D(self, c_mat, ybar, muR_ref): + """ Bai, Cogswell, Bazant 2011 """ + muRtheta = -self.eokT*3.422 + Ny = np.size(c_mat, 0) + dys = 1./Ny + Nx = np.size(c_mat, 1) + dxs = 1./Nx + muR_mat = np.empty((Ny,Nx), dtype=object) + actR_mat = np.empty((Ny,Nx), dtype=object) + shape = self.get_trode_param("shape") + kappa = self.get_trode_param("kappa") + B = self.get_trode_param("B") + cwet = self.get_trode_param("cwet") + # ytmp = np.empty(N+2, dtype=object) + # muRnonHomog_mat = self.general_non_homog(c_mat, ybar) + for i in range(Ny-2): + for j in range(Nx-2): + y = c_mat[i,j] + muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a")) + muR_mat[i,j] = muRhomog + + curvx = (c_mat[i,j] - c_mat[i,j+1]) - (c_mat[i,j+1] - c_mat[i,j+2]) + curvy = (c_mat[i,j] - c_mat[i+1,j]) - (c_mat[i+1,j] - c_mat[i+2,j]) + curv_en = -kappa*curvx/(dxs**2) -kappa*curvy/(dys**2) + muR_mat[i,j] += curv_en + muR_mat[i,j] += B*(y - ybar) + actR_mat[i,j] = np.exp(muR_mat[i,j]/self.T) + muR_mat[i,j] += muRtheta + muR_ref + return muR_mat, actR_mat diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 7e841566..d2ea7119 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -95,7 +95,7 @@ def DeclareEquations(self): "SM") for k in range(Nx): - eq = self.CreateEquation("avgs_cy_cx{k}".format(k=k)) + eq = self.CreateEquation("avgscy_isc{k}".format(k=k)) eq.Residual = self.c(k) for j in range(Ny): eq.Residual -= self.cy[k](j)/Ny @@ -109,13 +109,16 @@ def DeclareEquations(self): eq = self.CreateEquation("dcbardt") eq.Residual = self.dcbardt() for k in range(Nx): - eq.Residual -= self.c.dt(k)/Nx + for l in range(Ny): + eq.Residual -= (self.cy[k].dt(l)/Ny)/Nx + # eq.Residual -= self.c.dt(k)/Nx c_mat = np.empty((Ny, Nx), dtype=object) for k in range(Nx): c_mat[:,k] = [self.cy[k](j) for j in range(Ny)] - self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) + # self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) + self.sld_dynamics_2Dfull(c_mat, mu_O, act_lyte, self.noise) for eq in self.Equations: eq.CheckUnitsConsistency = False @@ -128,33 +131,82 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), Dfunc_name, f"mpet.electrode.diffusion.{Dfunc_name}") - dr, edges = geo.get_dr_edges("C3", Ny) + area_vec = 1. + Mmaty = get_Mmat("C3", Ny) + for k in range(Nx): + c_vec = c_mat[:,k] + muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), + self.config, self.trode, self.ind) + c_surf = c_mat[-1,k] + muR_surf = muR_vec[-1] + actR_surf = actR_vec[-1] + eta = calc_eta(muR_surf, muO) - c_surf = c_mat[-1,:] - muR_surf, actR_surf = calc_muR(c_surf, self.cbar(), - self.config, self.trode, self.ind) + eta_eff = eta + self.Rxn(k)*self.get_trode_param("Rfilm") - eta = calc_eta(muR_surf, muO) - eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm") - for i in range(Nx)]) - 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("lambda"), self.get_trode_param("alpha")) + 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("lambda"), self.get_trode_param("alpha")) - for k in range(Nx): eq = self.CreateEquation("Rxn_{k}".format(k=k)) - eq.Residual = self.Rxn(k) - Rxn[k] + eq.Residual = self.Rxn(k) - Rxn + + Flux_bc = -self.Rxn(k) + + Flux_vec = calc_flux_C3ver(c_vec, muR_vec, self.get_trode_param("D"), Dfunc, + self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + + RHS_vec = -np.diff(Flux_vec * area_vec) + dcdt_vec_y = np.empty(Ny, dtype=object) + dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] + LHS_vec_y = MX(Mmaty, dcdt_vec_y) + for j in range(Ny): + eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) + eq.Residual = LHS_vec_y[j] - RHS_vec[j] + + def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): + Ny = np.size(c_mat, 0) + Nx = np.size(c_mat, 1) + T = self.config["T"] + Dfunc_name = self.get_trode_param("Dfunc") + Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), + Dfunc_name, + f"mpet.electrode.diffusion.{Dfunc_name}") + dr, edges = geo.get_dr_edges("C3", Ny) area_vec = 1. Mmaty = get_Mmat("C3", Ny) + print('call calc_muR2D') + muR_mat , actR_mat = calc_muR2D(c_mat, self.cbar(), + self.config, self.trode, self.ind) for k in range(Nx): + c_vec = c_mat[:,k] + muR_vec = muR_mat[:,k] + actR_vec = actR_mat[:,k] + # muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), + # self.config, self.trode, self.ind) + c_surf = c_mat[-1,k] + muR_surf = muR_vec[-1] + actR_surf = actR_vec[-1] + eta = calc_eta(muR_surf, muO) + + eta_eff = eta + self.Rxn(k)*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("lambda"), self.get_trode_param("alpha")) + + eq = self.CreateEquation("Rxn_{k}".format(k=k)) + eq.Residual = self.Rxn(k) - Rxn + Flux_bc = -self.Rxn(k) - muR_vec, actR_vec = calc_muR(c_mat[:,k], self.cbar(), - self.config, self.trode, self.ind) - Flux_vec = calc_flux_C3ver(c_mat[:,k], muR_vec, self.get_trode_param("D"), Dfunc, + + Flux_vec = calc_flux_C3ver(c_vec, muR_vec, self.get_trode_param("D"), Dfunc, self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + RHS_vec = -np.diff(Flux_vec * area_vec) dcdt_vec_y = np.empty(Ny, dtype=object) dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] @@ -696,6 +748,11 @@ def calc_muR(c, cbar, config, trode, ind): muR, actR = muRfunc(c, cbar, muR_ref) return muR, actR +def calc_muR2D(c_mat, cbar, config, trode, ind): + muRfunc = props_am.muRfuncs(config, trode, ind).muRfunc + muR_ref = config[trode, "muR_ref"] + muR_mat, actR_mat = muRfunc(c_mat, cbar, muR_ref) + return muR_mat, actR_mat def MX(mat, objvec): if not isinstance(mat, sprs.csr.csr_matrix): diff --git a/mpet/sim.py b/mpet/sim.py index dbb41756..ac950878 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -91,7 +91,7 @@ def SetUpVariables(self): if solidType == "ACR2D": N_ver_ij = config["psd_num_ver"][tr][i,j] for k in range(Nij): - part.c.SetInitialCondition(k, cs0) + part.c.SetInitialGuess(k, cs0) for j in range(N_ver_ij): part.cy[k].SetInitialCondition(j, cs0) else: From 13aea3e8dd39a1e950318baf67ed016df7e19fe4 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Fri, 2 Dec 2022 13:07:52 +0000 Subject: [PATCH 11/27] update lfp2d --- configs/params_LFP2D.cfg | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index a4e4099e..51516aea 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -5,19 +5,21 @@ type = ACR2D discretization = 1e-9 shape = C3 -discretization_ver = 1.5e-9 -thickness = 20e-9 +discretization_ver = 5e-9 +thickness = 50e-9 [Material] -muRfunc = LiFePO4 +muRfunc = LiFePO42D noise = false noise_prefac = 1e-6 numnoise = 200 Omega_a = 1.8560e-20 +# Omega_a = 0.5560e-20 kappa = 5.0148e-10 B = 0.1916e9 +# B = 1e9 rho_s = 1.3793e28 -D = 5.3e-18 +D = 5.3e-17 Dfunc = lattice dgammadc = 0e-30 cwet = 0.98 From efc909d2cfc3eda45c40fbcd0f6010ab84041a7a Mon Sep 17 00:00:00 2001 From: Ombrini Date: Fri, 2 Dec 2022 17:10:29 +0100 Subject: [PATCH 12/27] small steps --- configs/params_LFP2D.cfg | 2 +- configs/params_system.cfg | 18 +++++------ mpet/config/schemas.py | 2 +- mpet/electrode/materials/LiFePO42D.py | 9 +++--- mpet/mod_electrodes.py | 45 +++++++++++++++++++++------ 5 files changed, 51 insertions(+), 25 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index 51516aea..e32435c4 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -9,7 +9,7 @@ discretization_ver = 5e-9 thickness = 50e-9 [Material] -muRfunc = LiFePO42D +muRfunc = LiFePO4 noise = false noise_prefac = 1e-6 numnoise = 200 diff --git a/configs/params_system.cfg b/configs/params_system.cfg index ef5dca6e..a53ce44d 100644 --- a/configs/params_system.cfg +++ b/configs/params_system.cfg @@ -13,7 +13,7 @@ Crate = 1 # 1C_current_density = 12.705 # Voltage cutoffs, V Vmax = 3.6 -Vmin = 2.0 +Vmin = 3.0 # Battery applied voltage (only used for CV), V Vset = 0.12 # Battery constant power (only used for CP), W/m^2 @@ -70,16 +70,16 @@ Rser = 0. # - If Nvol_c = 1 & Nvol_a = 0 & Nvol_s = 0, simulate a single volume with no # separator and infinitely fast counter electrode # - If Nvol_a = 0, simulate a Li foil electrode -Nvol_c = 10 -Nvol_s = 5 -Nvol_a = 10 +Nvol_c = 1 +Nvol_s = 0 +Nvol_a = 0 # Number of particles per volume for cathode and anode -Npart_c = 2 -Npart_a = 2 +Npart_c = 1 +Npart_a = 0 [Electrodes] # The name of the parameter file describing the cathode particles -cathode = params_LFP.cfg +cathode = params_LFP2D.cfg # The name of the parameter file describing the anode particles # Used only if Nvol_a > 0 anode = params_graphite_1param.cfg @@ -95,8 +95,8 @@ Rfilm_foil = 0e-0 # sphere or cylinder -- radius # If using stddev = 0, set Npart = 1 for that electrode to avoid # wasting computational effort. -mean_c = 100e-9 -stddev_c = 1e-9 +mean_c = 40e-9 +stddev_c = 0e-9 mean_a = 100e-9 stddev_a = 1e-9 # Use a specific set of particle sizes for the distribution diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py index 1f607d91..e75bb224 100644 --- a/mpet/config/schemas.py +++ b/mpet/config/schemas.py @@ -133,7 +133,7 @@ def tobool(value): 'discretization': Use(float), 'shape': lambda x: check_allowed_values(x, ["C3", "sphere", "cylinder", "homog_sdn"]), - Optional('discretization_ver', default = None):Use(float), + Optional('discretization_ver', default=None):Use(float), Optional('thickness'): Use(float)}, 'Material': {Optional('muRfunc_filename', default=None): str, 'muRfunc': str, diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py index e348af30..e3be2068 100644 --- a/mpet/electrode/materials/LiFePO42D.py +++ b/mpet/electrode/materials/LiFePO42D.py @@ -2,18 +2,18 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): - """ Bai, Cogswell, Bazant 2011 """ muRtheta = -self.eokT*3.422 + Ny = np.size(c_mat, 0) dys = 1./Ny Nx = np.size(c_mat, 1) dxs = 1./Nx muR_mat = np.empty((Ny,Nx), dtype=object) actR_mat = np.empty((Ny,Nx), dtype=object) - shape = self.get_trode_param("shape") + # shape = self.get_trode_param("shape") kappa = self.get_trode_param("kappa") B = self.get_trode_param("B") - cwet = self.get_trode_param("cwet") + # cwet = self.get_trode_param("cwet") # ytmp = np.empty(N+2, dtype=object) # muRnonHomog_mat = self.general_non_homog(c_mat, ybar) for i in range(Ny-2): @@ -21,10 +21,9 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): y = c_mat[i,j] muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a")) muR_mat[i,j] = muRhomog - curvx = (c_mat[i,j] - c_mat[i,j+1]) - (c_mat[i,j+1] - c_mat[i,j+2]) curvy = (c_mat[i,j] - c_mat[i+1,j]) - (c_mat[i+1,j] - c_mat[i+2,j]) - curv_en = -kappa*curvx/(dxs**2) -kappa*curvy/(dys**2) + curv_en = -kappa*curvx/(dxs**2) - kappa*curvy/(dys**2) muR_mat[i,j] += curv_en muR_mat[i,j] += B*(y - ybar) actR_mat[i,j] = np.exp(muR_mat[i,j]/self.T) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index d2ea7119..f69f3661 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -109,8 +109,8 @@ def DeclareEquations(self): eq = self.CreateEquation("dcbardt") eq.Residual = self.dcbardt() for k in range(Nx): - for l in range(Ny): - eq.Residual -= (self.cy[k].dt(l)/Ny)/Nx + for h in range(Ny): + eq.Residual -= (self.cy[k].dt(h)/Ny)/Nx # eq.Residual -= self.c.dt(k)/Nx c_mat = np.empty((Ny, Nx), dtype=object) @@ -166,7 +166,6 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) eq.Residual = LHS_vec_y[j] - RHS_vec[j] - def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): Ny = np.size(c_mat, 0) Nx = np.size(c_mat, 1) @@ -176,18 +175,18 @@ def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): Dfunc_name, f"mpet.electrode.diffusion.{Dfunc_name}") dr, edges = geo.get_dr_edges("C3", Ny) + # C3 shape has constant area along the depth area_vec = 1. Mmaty = get_Mmat("C3", Ny) - print('call calc_muR2D') - muR_mat , actR_mat = calc_muR2D(c_mat, self.cbar(), - self.config, self.trode, self.ind) + muR_mat, actR_mat = calc_muR2D(c_mat, self.cbar(), + self.config, self.trode, self.ind) for k in range(Nx): c_vec = c_mat[:,k] muR_vec = muR_mat[:,k] actR_vec = actR_mat[:,k] # muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), # self.config, self.trode, self.ind) - c_surf = c_mat[-1,k] + c_surf = c_vec[-1] muR_surf = muR_vec[-1] actR_surf = actR_vec[-1] eta = calc_eta(muR_surf, muO) @@ -748,12 +747,40 @@ def calc_muR(c, cbar, config, trode, ind): muR, actR = muRfunc(c, cbar, muR_ref) return muR, actR + def calc_muR2D(c_mat, cbar, config, trode, ind): - muRfunc = props_am.muRfuncs(config, trode, ind).muRfunc + eokT = 38.92714646 + T = 1. muR_ref = config[trode, "muR_ref"] - muR_mat, actR_mat = muRfunc(c_mat, cbar, muR_ref) + B = config[trode,'B'] + Omega = config[trode,'Omega_a'] + kappa = config[trode,'kappa'] + muRtheta = -eokT*3.422 + Ny = np.size(c_mat, 0) + dys = 1./Ny + Nx = np.size(c_mat, 1) + dxs = 1./Nx + muR_mat = np.empty((Ny,Nx), dtype=object) + actR_mat = np.empty((Ny,Nx), dtype=object) + for i in range(Ny-2): + for j in range(Nx-2): + y = c_mat[i,j] + muRhomog = T*np.log(y/(1-y)) + Omega*(1-2*y) + muR_mat[i,j] = muRhomog + curvx = (c_mat[i,j] - c_mat[i,j+1]) - (c_mat[i,j+1] - c_mat[i,j+2]) + curvy = (c_mat[i,j] - c_mat[i+1,j]) - (c_mat[i+1,j] - c_mat[i+2,j]) + curv_en = -kappa*curvx/(dxs**2) - kappa*curvy/(dys**2) + muR_mat[i,j] += curv_en + muR_mat[i,j] += B*(y - cbar) + actR_mat[i,j] = np.exp(muR_mat[i,j]/T) + muR_mat[i,j] += muRtheta + muR_ref + actR_mat[-2,-2] = actR_mat[-3,-3] + actR_mat[-1,-1] = actR_mat[-3,-3] + actR_mat[-2,-1] = actR_mat[-3,-3] + actR_mat[-1,-2] = actR_mat[-3,-3] return muR_mat, actR_mat + def MX(mat, objvec): if not isinstance(mat, sprs.csr.csr_matrix): raise Exception("MX function designed for csr mult") From ad37a0ed60da9372c2b7c9b2699a0ef904dfd0ae Mon Sep 17 00:00:00 2001 From: Ombrini Date: Wed, 7 Dec 2022 18:38:06 +0100 Subject: [PATCH 13/27] lfp2d now its in the correct position --- configs/params_LFP2D.cfg | 2 +- mpet/config/derived_values.py | 1 - mpet/electrode/materials/LiFePO42D.py | 47 ++++++++++++++----------- mpet/mod_electrodes.py | 49 ++++++--------------------- mpet/props_am.py | 1 + 5 files changed, 41 insertions(+), 59 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index e32435c4..51516aea 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -9,7 +9,7 @@ discretization_ver = 5e-9 thickness = 50e-9 [Material] -muRfunc = LiFePO4 +muRfunc = LiFePO42D noise = false noise_prefac = 1e-6 numnoise = 200 diff --git a/mpet/config/derived_values.py b/mpet/config/derived_values.py index fb18dd11..e27e7eda 100644 --- a/mpet/config/derived_values.py +++ b/mpet/config/derived_values.py @@ -229,7 +229,6 @@ def muR_ref(self, trode): muRfunc = props_am.muRfuncs(self.config, trode).muRfunc cs0bar = self.config['cs0'][trode] cs0 = np.array([cs0bar]) - solidType = self.config[trode, 'type'] if solidType in constants.two_var_types: muR_ref = -muRfunc((cs0, cs0), (cs0bar, cs0bar), 0.)[0][0] diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py index e3be2068..c142aefb 100644 --- a/mpet/electrode/materials/LiFePO42D.py +++ b/mpet/electrode/materials/LiFePO42D.py @@ -3,29 +3,38 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): muRtheta = -self.eokT*3.422 - + T = 1. + cbar = ybar + B = self.get_trode_param("B") + Omega = self.get_trode_param("Omega_a") + kappa = self.get_trode_param("kappa") + # for allowing muR_ref to be calculated + if np.size(c_mat) == 1: + c_mat_tmp = c_mat*np.ones((4,4)) + c_mat = c_mat_tmp Ny = np.size(c_mat, 0) - dys = 1./Ny Nx = np.size(c_mat, 1) + dys = 1./Ny dxs = 1./Nx + ywet = 0.98*np.ones(Ny+2, dtype=object) + c_mat_tmp = np.zeros((Ny+2,Nx+2), dtype=object) muR_mat = np.empty((Ny,Nx), dtype=object) actR_mat = np.empty((Ny,Nx), dtype=object) - # shape = self.get_trode_param("shape") - kappa = self.get_trode_param("kappa") - B = self.get_trode_param("B") - # cwet = self.get_trode_param("cwet") - # ytmp = np.empty(N+2, dtype=object) - # muRnonHomog_mat = self.general_non_homog(c_mat, ybar) - for i in range(Ny-2): - for j in range(Nx-2): - y = c_mat[i,j] - muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a")) - muR_mat[i,j] = muRhomog - curvx = (c_mat[i,j] - c_mat[i,j+1]) - (c_mat[i,j+1] - c_mat[i,j+2]) - curvy = (c_mat[i,j] - c_mat[i+1,j]) - (c_mat[i+1,j] - c_mat[i+2,j]) + c_mat_tmp[2:,1:-1] = c_mat + c_mat_tmp[:,-1] = ywet + c_mat_tmp[:,0] = ywet + c_mat_tmp[0,2:] = c_mat[0,:] + c_mat_tmp[1,2:] = c_mat[0,:] + for i in range(2,Ny+2): + for j in range(1,Nx+1): + y = c_mat_tmp[i,j] + muRhomog = T*np.log(y/(1-y)) + Omega*(1-2*y) + muR_mat[i-2,j-1] = muRhomog + curvx = (c_mat_tmp[i,j-1] - c_mat_tmp[i,j]) - (c_mat_tmp[i,j] - c_mat_tmp[i,j+1]) + curvy = (c_mat_tmp[i-2,j] - c_mat_tmp[i-1,j]) - (c_mat_tmp[i-1,j] - c_mat_tmp[i,j]) curv_en = -kappa*curvx/(dxs**2) - kappa*curvy/(dys**2) - muR_mat[i,j] += curv_en - muR_mat[i,j] += B*(y - ybar) - actR_mat[i,j] = np.exp(muR_mat[i,j]/self.T) - muR_mat[i,j] += muRtheta + muR_ref + muR_mat[i-2,j-1] += curv_en + muR_mat[i-2,j-1] += B*(y - cbar) + actR_mat[i-2,j-1] = np.exp(muR_mat[i-2,j-1]/T) + muR_mat[i-2,j-1] += muRtheta + muR_ref return muR_mat, actR_mat diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index f69f3661..7bd89205 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -178,19 +178,25 @@ def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): # C3 shape has constant area along the depth area_vec = 1. Mmaty = get_Mmat("C3", Ny) - muR_mat, actR_mat = calc_muR2D(c_mat, self.cbar(), - self.config, self.trode, self.ind) + print('cmat shape ', np.shape(c_mat)) + muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), + self.config, self.trode, self.ind) + print('muR_mat size', np.shape(muR_mat)) for k in range(Nx): c_vec = c_mat[:,k] muR_vec = muR_mat[:,k] + print('c_vec size ', np.shape(c_vec)) + print('muR_vec size ', np.shape(muR_vec)) + # print(muR_vec) actR_vec = actR_mat[:,k] - # muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), - # self.config, self.trode, self.ind) + c_surf = c_vec[-1] muR_surf = muR_vec[-1] + print('muR surf size ' , np.shape(muR_surf)) + print('c_surf surf size ' , np.shape(c_surf)) + # print(muR_surf) actR_surf = actR_vec[-1] eta = calc_eta(muR_surf, muO) - eta_eff = eta + self.Rxn(k)*self.get_trode_param("Rfilm") Rxn = self.calc_rxn_rate( @@ -748,39 +754,6 @@ def calc_muR(c, cbar, config, trode, ind): return muR, actR -def calc_muR2D(c_mat, cbar, config, trode, ind): - eokT = 38.92714646 - T = 1. - muR_ref = config[trode, "muR_ref"] - B = config[trode,'B'] - Omega = config[trode,'Omega_a'] - kappa = config[trode,'kappa'] - muRtheta = -eokT*3.422 - Ny = np.size(c_mat, 0) - dys = 1./Ny - Nx = np.size(c_mat, 1) - dxs = 1./Nx - muR_mat = np.empty((Ny,Nx), dtype=object) - actR_mat = np.empty((Ny,Nx), dtype=object) - for i in range(Ny-2): - for j in range(Nx-2): - y = c_mat[i,j] - muRhomog = T*np.log(y/(1-y)) + Omega*(1-2*y) - muR_mat[i,j] = muRhomog - curvx = (c_mat[i,j] - c_mat[i,j+1]) - (c_mat[i,j+1] - c_mat[i,j+2]) - curvy = (c_mat[i,j] - c_mat[i+1,j]) - (c_mat[i+1,j] - c_mat[i+2,j]) - curv_en = -kappa*curvx/(dxs**2) - kappa*curvy/(dys**2) - muR_mat[i,j] += curv_en - muR_mat[i,j] += B*(y - cbar) - actR_mat[i,j] = np.exp(muR_mat[i,j]/T) - muR_mat[i,j] += muRtheta + muR_ref - actR_mat[-2,-2] = actR_mat[-3,-3] - actR_mat[-1,-1] = actR_mat[-3,-3] - actR_mat[-2,-1] = actR_mat[-3,-3] - actR_mat[-1,-2] = actR_mat[-3,-3] - return muR_mat, actR_mat - - def MX(mat, objvec): if not isinstance(mat, sprs.csr.csr_matrix): raise Exception("MX function designed for csr mult") diff --git a/mpet/props_am.py b/mpet/props_am.py index 41dd6398..cd2f35eb 100644 --- a/mpet/props_am.py +++ b/mpet/props_am.py @@ -149,6 +149,7 @@ def non_homog_rect_fixed_csurf(self, y, ybar, B, kappa, ywet): ytmp[1:-1] = y ytmp[0] = ywet ytmp[-1] = ywet + print(ytmp) dxs = 1./N curv = np.diff(ytmp, 2)/(dxs**2) muR_nh = -kappa*curv + B*(y - ybar) From c2b19ee90cf8e56b02d3357922c5ea95a22d1d0c Mon Sep 17 00:00:00 2001 From: Ombrini Date: Wed, 7 Dec 2022 18:47:32 +0100 Subject: [PATCH 14/27] check direction --- mpet/mod_electrodes.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 7bd89205..663ed868 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -178,24 +178,22 @@ def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): # C3 shape has constant area along the depth area_vec = 1. Mmaty = get_Mmat("C3", Ny) - print('cmat shape ', np.shape(c_mat)) muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), self.config, self.trode, self.ind) - print('muR_mat size', np.shape(muR_mat)) for k in range(Nx): c_vec = c_mat[:,k] + c_vec = np.reshape(c_vec, (Ny,1)) muR_vec = muR_mat[:,k] - print('c_vec size ', np.shape(c_vec)) - print('muR_vec size ', np.shape(muR_vec)) + muR_vec = np.reshape(muR_vec, (Ny,1)) # print(muR_vec) actR_vec = actR_mat[:,k] + actR_vec = np.reshape(actR_vec, (Ny,1)) - c_surf = c_vec[-1] - muR_surf = muR_vec[-1] - print('muR surf size ' , np.shape(muR_surf)) - print('c_surf surf size ' , np.shape(c_surf)) + c_surf = c_vec[-1,:] + muR_surf = muR_vec[-1,:] + print(muR_surf) # print(muR_surf) - actR_surf = actR_vec[-1] + actR_surf = actR_vec[-1,:] eta = calc_eta(muR_surf, muO) eta_eff = eta + self.Rxn(k)*self.get_trode_param("Rfilm") From 72ad007b35c32bdba2433042865f06f88407ba27 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 8 Dec 2022 17:14:23 +0100 Subject: [PATCH 15/27] 2d lfp no mechanics --- mpet/electrode/materials/LiFePO42D.py | 80 +++++++++++++++++---------- mpet/mod_electrodes.py | 38 +++++++------ mpet/props_am.py | 1 - 3 files changed, 71 insertions(+), 48 deletions(-) diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py index c142aefb..91934885 100644 --- a/mpet/electrode/materials/LiFePO42D.py +++ b/mpet/electrode/materials/LiFePO42D.py @@ -7,34 +7,56 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): cbar = ybar B = self.get_trode_param("B") Omega = self.get_trode_param("Omega_a") - kappa = self.get_trode_param("kappa") + kappax = self.get_trode_param("kappa") + kappay = kappax*10 + # beta_s = self.get_trode_param("beta_s") # for allowing muR_ref to be calculated if np.size(c_mat) == 1: - c_mat_tmp = c_mat*np.ones((4,4)) - c_mat = c_mat_tmp - Ny = np.size(c_mat, 0) - Nx = np.size(c_mat, 1) - dys = 1./Ny - dxs = 1./Nx - ywet = 0.98*np.ones(Ny+2, dtype=object) - c_mat_tmp = np.zeros((Ny+2,Nx+2), dtype=object) - muR_mat = np.empty((Ny,Nx), dtype=object) - actR_mat = np.empty((Ny,Nx), dtype=object) - c_mat_tmp[2:,1:-1] = c_mat - c_mat_tmp[:,-1] = ywet - c_mat_tmp[:,0] = ywet - c_mat_tmp[0,2:] = c_mat[0,:] - c_mat_tmp[1,2:] = c_mat[0,:] - for i in range(2,Ny+2): - for j in range(1,Nx+1): - y = c_mat_tmp[i,j] - muRhomog = T*np.log(y/(1-y)) + Omega*(1-2*y) - muR_mat[i-2,j-1] = muRhomog - curvx = (c_mat_tmp[i,j-1] - c_mat_tmp[i,j]) - (c_mat_tmp[i,j] - c_mat_tmp[i,j+1]) - curvy = (c_mat_tmp[i-2,j] - c_mat_tmp[i-1,j]) - (c_mat_tmp[i-1,j] - c_mat_tmp[i,j]) - curv_en = -kappa*curvx/(dxs**2) - kappa*curvy/(dys**2) - muR_mat[i-2,j-1] += curv_en - muR_mat[i-2,j-1] += B*(y - cbar) - actR_mat[i-2,j-1] = np.exp(muR_mat[i-2,j-1]/T) - muR_mat[i-2,j-1] += muRtheta + muR_ref - return muR_mat, actR_mat + y = c_mat + muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a")) + muRnonHomog = self.general_non_homog(y, ybar) + muR = muRhomog + muRnonHomog + actR = np.exp(muR/self.T) + muR += muRtheta + muR_ref + return muR, actR + else: + Ny = np.size(c_mat, 1) + Nx = np.size(c_mat, 0) + dys = 1./Ny + dxs = 1./Nx + ywet = 0.98*np.ones(Ny, dtype=object) + c_mat_tmp = np.zeros((Nx+2,Ny), dtype=object) + muR_mat = np.zeros((Nx,Ny), dtype=object) + actR_mat = np.zeros((Nx,Ny), dtype=object) + # from second to second to last row we have the original conc + c_mat_tmp[1:-1,:] = c_mat + # first and last row is 0.98 + c_mat_tmp[-1,:] = ywet + c_mat_tmp[0,:] = ywet + curvy = np.empty(Ny, dtype=object) + # print(muR_mat) + # plan: produce the matrix muR without curvy + # then another for loop for the curvy + # if then works try to optmize the loops + for i in range(1,Nx+1): + curvy[0] = 2*c_mat_tmp[i,1] - 2*c_mat_tmp[i,0] + curvy[-1] = 2*c_mat_tmp[i,-2] - 2*c_mat_tmp[i,-1] \ + # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) + curvy[1:Ny-1] = np.diff(c_mat_tmp[i,:],2) + for j in range(Ny): + y = c_mat_tmp[i,j] + muRhomog = T*np.log(y/(1-y)) + Omega*(1-2*y) + muR_mat[i-1,j] = muRhomog + curvx = (c_mat_tmp[i-1,j] - c_mat_tmp[i,j]) - (c_mat_tmp[i,j] - c_mat_tmp[i+1,j]) + # curvy[i-1,j] = (c_mat_tmp[i-1,j-1] - c_mat_tmp[i-1,j]) - (c_mat_tmp[i-1,j] + # - c_mat_tmp[i-1,j+1]) + # curvy[i-1,-1] = 2*c_mat_tmp[i-1,-2] - 2*c_mat_tmp[i-1,-1] + # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) + # curv_en = -kappax*curvx/(dxs**2) - kappay*curvy/(dys**2) + curv_en = -kappax*curvx/(dxs**2) - kappay*curvy[j]/(dys**2) + muR_mat[i-1,j] += curv_en + muR_mat[i-1,j] += B*(y - cbar) + # print(muR_mat) + actR_mat = np.exp(muR_mat/T) + muR_mat += muRtheta + muR_ref + return muR_mat, actR_mat diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 663ed868..0a27af52 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -113,9 +113,9 @@ def DeclareEquations(self): eq.Residual -= (self.cy[k].dt(h)/Ny)/Nx # eq.Residual -= self.c.dt(k)/Nx - c_mat = np.empty((Ny, Nx), dtype=object) + c_mat = np.empty((Nx, Ny), dtype=object) for k in range(Nx): - c_mat[:,k] = [self.cy[k](j) for j in range(Ny)] + c_mat[k,:] = [self.cy[k](j) for j in range(Ny)] # self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) self.sld_dynamics_2Dfull(c_mat, mu_O, act_lyte, self.noise) @@ -124,8 +124,8 @@ def DeclareEquations(self): eq.CheckUnitsConsistency = False def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): - Ny = np.size(c_mat, 0) - Nx = np.size(c_mat, 1) + Ny = np.size(c_mat, 1) + Nx = np.size(c_mat, 0) T = self.config["T"] Dfunc_name = self.get_trode_param("Dfunc") Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), @@ -135,10 +135,10 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): area_vec = 1. Mmaty = get_Mmat("C3", Ny) for k in range(Nx): - c_vec = c_mat[:,k] + c_vec = c_mat[k,:] muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), self.config, self.trode, self.ind) - c_surf = c_mat[-1,k] + c_surf = c_mat[k,-1] muR_surf = muR_vec[-1] actR_surf = actR_vec[-1] eta = calc_eta(muR_surf, muO) @@ -167,8 +167,8 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq.Residual = LHS_vec_y[j] - RHS_vec[j] def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): - Ny = np.size(c_mat, 0) - Nx = np.size(c_mat, 1) + Ny = np.size(c_mat, 1) + Nx = np.size(c_mat, 0) T = self.config["T"] Dfunc_name = self.get_trode_param("Dfunc") Dfunc = utils.import_function(self.get_trode_param("Dfunc_filename"), @@ -178,22 +178,24 @@ def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): # C3 shape has constant area along the depth area_vec = 1. Mmaty = get_Mmat("C3", Ny) + # print(c_mat) muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), self.config, self.trode, self.ind) for k in range(Nx): - c_vec = c_mat[:,k] - c_vec = np.reshape(c_vec, (Ny,1)) - muR_vec = muR_mat[:,k] - muR_vec = np.reshape(muR_vec, (Ny,1)) + c_vec = c_mat[k,:] + # print(c_vec) + # c_vec = np.reshape(c_vec, (Ny,1)) + muR_vec = muR_mat[k,:] + # muR_vec = np.reshape(muR_vec, (Ny,1)) # print(muR_vec) - actR_vec = actR_mat[:,k] - actR_vec = np.reshape(actR_vec, (Ny,1)) + actR_vec = actR_mat[k,:] + # actR_vec = np.reshape(actR_vec, (Ny,1)) - c_surf = c_vec[-1,:] - muR_surf = muR_vec[-1,:] - print(muR_surf) + c_surf = c_vec[-1] + # print(c_surf) + muR_surf = muR_vec[-1] # print(muR_surf) - actR_surf = actR_vec[-1,:] + actR_surf = actR_vec[-1] eta = calc_eta(muR_surf, muO) eta_eff = eta + self.Rxn(k)*self.get_trode_param("Rfilm") diff --git a/mpet/props_am.py b/mpet/props_am.py index cd2f35eb..41dd6398 100644 --- a/mpet/props_am.py +++ b/mpet/props_am.py @@ -149,7 +149,6 @@ def non_homog_rect_fixed_csurf(self, y, ybar, B, kappa, ywet): ytmp[1:-1] = y ytmp[0] = ywet ytmp[-1] = ywet - print(ytmp) dxs = 1./N curv = np.diff(ytmp, 2)/(dxs**2) muR_nh = -kappa*curv + B*(y - ybar) From fc40d20d10696ccf4051bd43162c3d98f291aff4 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 8 Dec 2022 17:34:56 +0100 Subject: [PATCH 16/27] style schemas fix --- mpet/config/schemas.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py index e75bb224..142cc868 100644 --- a/mpet/config/schemas.py +++ b/mpet/config/schemas.py @@ -133,8 +133,8 @@ def tobool(value): 'discretization': Use(float), 'shape': lambda x: check_allowed_values(x, ["C3", "sphere", "cylinder", "homog_sdn"]), - Optional('discretization_ver', default=None):Use(float), - Optional('thickness'): Use(float)}, + Optional('discretization_ver', default=None):Use(float), + Optional('thickness'): Use(float)}, 'Material': {Optional('muRfunc_filename', default=None): str, 'muRfunc': str, 'noise': Use(tobool), From 5faf4c1ef79b4004a55ac08324eca78ee0044cf8 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Mon, 12 Dec 2022 09:28:48 +0000 Subject: [PATCH 17/27] base of mechanical variables --- mpet/mod_electrodes.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 0a27af52..7f962169 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -49,6 +49,18 @@ def __init__(self, config, trode, vInd, pInd, self.cy[k] = dae.daeVariable("cy{k}".format(k=k), mole_frac_t, self, "Concentration in y direction of element {k}".format(k=k), [self.Dmny]) + # check unit of measures + self.uy = {} + for k in range(Nx): + self.uy[k] = dae.daeVariable("uy{k}".format(k=k), mole_frac_t, self, + "Displacement in y direction of element {k}".format(k=k), + [self.Dmny]) + self.ux = {} + Ny = np.size(self.Dmny) + for k in range(Ny): + self.ux[k] = dae.daeVariable("ux{k}".format(k=k), mole_frac_t, self, + "Displacement in x direction of element {k}".format(k=k), + [self.Dmn]) self.cbar = dae.daeVariable( "cbar", mole_frac_t, self, @@ -113,12 +125,21 @@ def DeclareEquations(self): eq.Residual -= (self.cy[k].dt(h)/Ny)/Nx # eq.Residual -= self.c.dt(k)/Nx + c_mat = np.empty((Nx, Ny), dtype=object) for k in range(Nx): c_mat[k,:] = [self.cy[k](j) for j in range(Ny)] + # to check !! + u_y_mat = np.empty((Nx, Ny), dtype=object) + for k in range(Nx): + u_y_mat[k,:] = [self.uy[k](j) for j in range(Ny)] + u_x_mat = np.empty((Nx, Ny), dtype=object) + for k in range(Ny): + u_x_mat[:,k] = [self.ux[k](j) for j in range(Nx)] + # self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) - self.sld_dynamics_2Dfull(c_mat, mu_O, act_lyte, self.noise) + self.sld_dynamics_2Dfull(c_mat, u_x_mat, u_y_mat, mu_O, act_lyte, self.noise) for eq in self.Equations: eq.CheckUnitsConsistency = False @@ -166,7 +187,7 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq = self.CreateEquation("dcydt_{k}_{j}".format(k=k, j=j)) eq.Residual = LHS_vec_y[j] - RHS_vec[j] - def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): + def sld_dynamics_2Dfull(self, c_mat, u_x_mat, u_y_mat, muO, act_lyte, noise): Ny = np.size(c_mat, 1) Nx = np.size(c_mat, 0) T = self.config["T"] @@ -181,6 +202,8 @@ def sld_dynamics_2Dfull(self, c_mat, muO, act_lyte, noise): # print(c_mat) muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), self.config, self.trode, self.ind) + # muR_mat, actR_mat = calc_muR_mech(c_mat, u_x_mat, u_y_mat, self.cbar(), + # self.config, self.trode, self.ind) for k in range(Nx): c_vec = c_mat[k,:] # print(c_vec) From 982035979f3cc391df915ab4ae02d4699bd3cfbd Mon Sep 17 00:00:00 2001 From: Ombrini Date: Mon, 12 Dec 2022 13:00:33 +0100 Subject: [PATCH 18/27] elastic constant in the code chem pot --- mpet/electrode/materials/LiFePO42D.py | 50 ++++++++++++++++++++------- 1 file changed, 37 insertions(+), 13 deletions(-) diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py index 91934885..b04c6b96 100644 --- a/mpet/electrode/materials/LiFePO42D.py +++ b/mpet/electrode/materials/LiFePO42D.py @@ -3,14 +3,43 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): muRtheta = -self.eokT*3.422 - T = 1. - cbar = ybar - B = self.get_trode_param("B") Omega = self.get_trode_param("Omega_a") kappax = self.get_trode_param("kappa") kappay = kappax*10 # beta_s = self.get_trode_param("beta_s") # for allowing muR_ref to be calculated + # FePo4 elastic constants (GPa) + c11 = 175.9 + c22 = 153.6 + c33 = 135.0 + c44 = 38.8 + c55 = 47.5 + c66 = 55.6 + c13 = 54.0 + c12 = 29.6 + c23 = 19.6 + + Cij = np.zeros((6,6)) + Cij[0,0] = c11 + Cij[1,1] = c22 + Cij[2,2] = c33 + Cij[3,3] = c44 + Cij[4,4] = c55 + Cij[5,5] = c66 + Cij[1,0] = c12 + Cij[0,1] = c12 + Cij[2,0] = c13 + Cij[0,2] = c13 + Cij[1,2] = c23 + Cij[2,1] = c23 + # strain + e1 = 0.0517 + e2 = 0.0359 + e3 = -0.0186 + ej = np.array([e1, e2, e3, 0, 0, 0]) + sigma_i = np.dot(Cij,ej) + print(sigma_i) + if np.size(c_mat) == 1: y = c_mat muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a")) @@ -40,23 +69,18 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): # if then works try to optmize the loops for i in range(1,Nx+1): curvy[0] = 2*c_mat_tmp[i,1] - 2*c_mat_tmp[i,0] - curvy[-1] = 2*c_mat_tmp[i,-2] - 2*c_mat_tmp[i,-1] \ - # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) + curvy[-1] = 2*c_mat_tmp[i,-2] - 2*c_mat_tmp[i,-1] + # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) curvy[1:Ny-1] = np.diff(c_mat_tmp[i,:],2) for j in range(Ny): y = c_mat_tmp[i,j] - muRhomog = T*np.log(y/(1-y)) + Omega*(1-2*y) + muRhomog = self.T*np.log(y/(1-y)) + Omega*(1-2*y) muR_mat[i-1,j] = muRhomog curvx = (c_mat_tmp[i-1,j] - c_mat_tmp[i,j]) - (c_mat_tmp[i,j] - c_mat_tmp[i+1,j]) - # curvy[i-1,j] = (c_mat_tmp[i-1,j-1] - c_mat_tmp[i-1,j]) - (c_mat_tmp[i-1,j] - # - c_mat_tmp[i-1,j+1]) - # curvy[i-1,-1] = 2*c_mat_tmp[i-1,-2] - 2*c_mat_tmp[i-1,-1] - # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) - # curv_en = -kappax*curvx/(dxs**2) - kappay*curvy/(dys**2) curv_en = -kappax*curvx/(dxs**2) - kappay*curvy[j]/(dys**2) muR_mat[i-1,j] += curv_en - muR_mat[i-1,j] += B*(y - cbar) + muR_mat[i-1,j] += self.get_trode_param("B")*(y - ybar) # print(muR_mat) - actR_mat = np.exp(muR_mat/T) + actR_mat = np.exp(muR_mat/self.T) muR_mat += muRtheta + muR_ref return muR_mat, actR_mat From 64ddfda8b2ce8a867c2a56afd2737025efc16111 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Mon, 12 Dec 2022 13:01:07 +0100 Subject: [PATCH 19/27] config changes --- configs/params_LFP2D.cfg | 4 ++-- configs/params_system.cfg | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index 51516aea..6d7335c0 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -5,8 +5,8 @@ type = ACR2D discretization = 1e-9 shape = C3 -discretization_ver = 5e-9 -thickness = 50e-9 +discretization_ver = 1e-9 +thickness = 10e-9 [Material] muRfunc = LiFePO42D diff --git a/configs/params_system.cfg b/configs/params_system.cfg index a53ce44d..721eaaa1 100644 --- a/configs/params_system.cfg +++ b/configs/params_system.cfg @@ -13,7 +13,7 @@ Crate = 1 # 1C_current_density = 12.705 # Voltage cutoffs, V Vmax = 3.6 -Vmin = 3.0 +Vmin = 3.3 # Battery applied voltage (only used for CV), V Vset = 0.12 # Battery constant power (only used for CP), W/m^2 @@ -62,7 +62,7 @@ seed = 0 # printing internal variable concentrations) files. hdf5 files # are better for cycling, as they store less information and there is less # opening/rewriting of files. Default is mat -dataReporter = hdf5 +dataReporter = mat # Series resistance, [Ohm m^2] Rser = 0. # Cathode, anode, and separator numer disc. in x direction (volumes in electrodes) @@ -95,8 +95,8 @@ Rfilm_foil = 0e-0 # sphere or cylinder -- radius # If using stddev = 0, set Npart = 1 for that electrode to avoid # wasting computational effort. -mean_c = 40e-9 -stddev_c = 0e-9 +mean_c = 10e-9 +stddev_c = 10e-9 mean_a = 100e-9 stddev_a = 1e-9 # Use a specific set of particle sizes for the distribution From 9b9f8b81ed1c614bf2e42770da9827bffe91d699 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Mon, 12 Dec 2022 14:06:00 +0000 Subject: [PATCH 20/27] start of a divergence calculation --- mpet/electrode/materials/LiFePO42D.py | 32 ------- mpet/mod_electrodes.py | 115 +++++++++++++++++++++++--- 2 files changed, 103 insertions(+), 44 deletions(-) diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py index b04c6b96..9da50bc8 100644 --- a/mpet/electrode/materials/LiFePO42D.py +++ b/mpet/electrode/materials/LiFePO42D.py @@ -8,38 +8,6 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): kappay = kappax*10 # beta_s = self.get_trode_param("beta_s") # for allowing muR_ref to be calculated - # FePo4 elastic constants (GPa) - c11 = 175.9 - c22 = 153.6 - c33 = 135.0 - c44 = 38.8 - c55 = 47.5 - c66 = 55.6 - c13 = 54.0 - c12 = 29.6 - c23 = 19.6 - - Cij = np.zeros((6,6)) - Cij[0,0] = c11 - Cij[1,1] = c22 - Cij[2,2] = c33 - Cij[3,3] = c44 - Cij[4,4] = c55 - Cij[5,5] = c66 - Cij[1,0] = c12 - Cij[0,1] = c12 - Cij[2,0] = c13 - Cij[0,2] = c13 - Cij[1,2] = c23 - Cij[2,1] = c23 - # strain - e1 = 0.0517 - e2 = 0.0359 - e3 = -0.0186 - ej = np.array([e1, e2, e3, 0, 0, 0]) - sigma_i = np.dot(Cij,ej) - print(sigma_i) - if np.size(c_mat) == 1: y = c_mat muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a")) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 7f962169..010d9259 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -47,20 +47,19 @@ def __init__(self, config, trode, vInd, pInd, self.cy = {} for k in range(Nx): self.cy[k] = dae.daeVariable("cy{k}".format(k=k), mole_frac_t, self, - "Concentration in y direction of element {k}".format(k=k), + "Concentration in y direction of element in row {k}".format(k=k), [self.Dmny]) # check unit of measures self.uy = {} for k in range(Nx): self.uy[k] = dae.daeVariable("uy{k}".format(k=k), mole_frac_t, self, - "Displacement in y direction of element {k}".format(k=k), + "Displacement in y direction of element in row {k}".format(k=k), [self.Dmny]) self.ux = {} - Ny = np.size(self.Dmny) - for k in range(Ny): + for k in range(Nx): self.ux[k] = dae.daeVariable("ux{k}".format(k=k), mole_frac_t, self, - "Displacement in x direction of element {k}".format(k=k), - [self.Dmn]) + "Displacement in x direction of element in row {k}".format(k=k), + [self.Dmny]) self.cbar = dae.daeVariable( "cbar", mole_frac_t, self, @@ -125,18 +124,17 @@ def DeclareEquations(self): eq.Residual -= (self.cy[k].dt(h)/Ny)/Nx # eq.Residual -= self.c.dt(k)/Nx - c_mat = np.empty((Nx, Ny), dtype=object) for k in range(Nx): c_mat[k,:] = [self.cy[k](j) for j in range(Ny)] - # to check !! u_y_mat = np.empty((Nx, Ny), dtype=object) for k in range(Nx): u_y_mat[k,:] = [self.uy[k](j) for j in range(Ny)] + u_x_mat = np.empty((Nx, Ny), dtype=object) - for k in range(Ny): - u_x_mat[:,k] = [self.ux[k](j) for j in range(Nx)] + for k in range(Nx): + u_x_mat[k,:] = [self.ux[k](j) for j in range(Ny)] # self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) self.sld_dynamics_2Dfull(c_mat, u_x_mat, u_y_mat, mu_O, act_lyte, self.noise) @@ -202,8 +200,10 @@ def sld_dynamics_2Dfull(self, c_mat, u_x_mat, u_y_mat, muO, act_lyte, noise): # print(c_mat) muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), self.config, self.trode, self.ind) - # muR_mat, actR_mat = calc_muR_mech(c_mat, u_x_mat, u_y_mat, self.cbar(), - # self.config, self.trode, self.ind) + muR_el = calc_muR_el(c_mat, u_x_mat, u_y_mat, + self.config, self.trode, self.ind) + + muR_mat += muR_el for k in range(Nx): c_vec = c_mat[k,:] # print(c_vec) @@ -776,6 +776,97 @@ def calc_muR(c, cbar, config, trode, ind): muR, actR = muRfunc(c, cbar, muR_ref) return muR, actR +def mech_tensors(): + # FePo4 elastic constants (GPa) + c11 = 175.9 + c22 = 153.6 + c33 = 135.0 + c44 = 38.8 + c55 = 47.5 + c66 = 55.6 + c13 = 54.0 + c12 = 29.6 + c23 = 19.6 + + Cij = np.zeros((6,6)) + Cij[0,0] = c11 + Cij[1,1] = c22 + Cij[2,2] = c33 + Cij[3,3] = c44 + Cij[4,4] = c55 + Cij[5,5] = c66 + Cij[1,0] = c12 + Cij[0,1] = c12 + Cij[2,0] = c13 + Cij[0,2] = c13 + Cij[1,2] = c23 + Cij[2,1] = c23 + # strain + e01 = 0.0517 + e02 = 0.0359 + e03 = -0.0186 + e0 = np.array([e01, e02, e03, 0, 0, 0]) + + return Cij, e0 + +def calc_muR_el(c_mat, u_x, u_y, config, trode, ind): + Ny = np.size(c_mat, 1) + Nx = np.size(c_mat, 0) + dys = 1./Ny + dxs = 1./Nx + max_conc = config[trode, "rho_s"] + muR_el = np.zeros((Nx,Ny), dtype=object) + Cij, e0 = mech_tensors() + + e1_mat = np.zeros((Nx,Ny), dtype=object) + e2_mat = np.zeros((Nx,Ny), dtype=object) + duxdy = np.zeros((Nx,Ny), dtype=object) + duydx = np.zeros((Nx,Ny), dtype=object) + + # check boundaries !! + for j in range(Ny): + e1_mat[:,j] = np.diff(u_x[:,j])/dxs + duydx[:,j] = np.diff(u_y[:,j])/dxs + for i in range(Nx): + e2_mat[i,:] = np.diff(u_y[i,:])/dys + duxdy[i,:] = np.diff(u_x[i,:])/dys + e12_mat = 0.5*(duxdy+duydx) + + e_vec_of_ij = np.zeros(6, dtype=object) + sigma_vec_of_ij = np.zeros(6, dtype=object) + + sigma_11_mat = np.zeros((Nx,Ny), dtype=object) + sigma_22_mat = np.zeros((Nx,Ny), dtype=object) + sigma_33_mat = np.zeros((Nx,Ny), dtype=object) + sigma_12_mat = np.zeros((Nx,Ny), dtype=object) + sigma_13_mat = np.zeros((Nx,Ny), dtype=object) + sigma_23_mat = np.zeros((Nx,Ny), dtype=object) + div_sigm_1_mat = np.zeros((Nx,Ny), dtype=object) + div_sigm_2_mat = np.zeros((Nx,Ny), dtype=object) + div_sigm_3_mat = np.zeros((Nx,Ny), dtype=object) + for i in range(Nx): + for j in range(Ny): + e_vec_of_ij = np.array([e1_mat[i,j]- e0[0]*c_mat[i,j], + e2_mat[i,j]- e0[1]*c_mat[i,j], + 0, + e12_mat[i,j], + 0, + 0]) + sigma_vec_of_ij = np.dot((Cij,e_vec_of_ij)) + sigma_11_mat[i,j] = sigma_vec_of_ij[0] + sigma_22_mat[i,j] = sigma_vec_of_ij[1] + sigma_33_mat[i,j] = sigma_vec_of_ij[2] + sigma_12_mat[i,j] = sigma_vec_of_ij[3] + sigma_13_mat[i,j] = sigma_vec_of_ij[4] + sigma_23_mat[i,j] = sigma_vec_of_ij[5] + muR_el[i,j] = np.dot((sigma_vec_of_ij,e0))/max_conc + for j in range(Ny): + dsigma1dx = np.diff(sigma_11_mat[:,j])/dxs + dsgima12dx = np.diff(sigma_12_mat[:,j])/dxs + + + return muR_el, div_stress_mat + def MX(mat, objvec): if not isinstance(mat, sprs.csr.csr_matrix): From 5bc760922267e3e83a3c377ca3eb65aa2396e1ed Mon Sep 17 00:00:00 2001 From: Ombrini Date: Mon, 12 Dec 2022 18:36:28 +0000 Subject: [PATCH 21/27] missing boundary condi on sigma --- configs/params_LFP2D.cfg | 2 +- configs/params_system.cfg | 2 +- mpet/mod_electrodes.py | 108 ++++++++++++++++++++++++++------------ mpet/sim.py | 2 + 4 files changed, 77 insertions(+), 37 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index 6d7335c0..a23a255e 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -6,7 +6,7 @@ type = ACR2D discretization = 1e-9 shape = C3 discretization_ver = 1e-9 -thickness = 10e-9 +thickness = 2e-9 [Material] muRfunc = LiFePO42D diff --git a/configs/params_system.cfg b/configs/params_system.cfg index 721eaaa1..00997d71 100644 --- a/configs/params_system.cfg +++ b/configs/params_system.cfg @@ -95,7 +95,7 @@ Rfilm_foil = 0e-0 # sphere or cylinder -- radius # If using stddev = 0, set Npart = 1 for that electrode to avoid # wasting computational effort. -mean_c = 10e-9 +mean_c = 4e-9 stddev_c = 10e-9 mean_a = 100e-9 stddev_a = 1e-9 diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 010d9259..c94b9093 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -200,10 +200,19 @@ def sld_dynamics_2Dfull(self, c_mat, u_x_mat, u_y_mat, muO, act_lyte, noise): # print(c_mat) muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), self.config, self.trode, self.ind) - muR_el = calc_muR_el(c_mat, u_x_mat, u_y_mat, + muR_el, div_stress_mat = calc_muR_el(c_mat, u_x_mat, u_y_mat, self.config, self.trode, self.ind) muR_mat += muR_el + actR_mat = np.exp(muR_mat) + + for i in range(Nx): + for j in range(Ny): + eq1 = self.CreateEquation("divsigma1_{i}_{j}_equal0".format(i=i, j=j)) + eq1.Residual = div_stress_mat[i,j,0] + eq2 = self.CreateEquation("divsigma2_{i}_{j}_equal0".format(i=i, j=j)) + eq2.Residual = div_stress_mat[i,j,1] + for k in range(Nx): c_vec = c_mat[k,:] # print(c_vec) @@ -818,53 +827,82 @@ def calc_muR_el(c_mat, u_x, u_y, config, trode, ind): muR_el = np.zeros((Nx,Ny), dtype=object) Cij, e0 = mech_tensors() - e1_mat = np.zeros((Nx,Ny), dtype=object) - e2_mat = np.zeros((Nx,Ny), dtype=object) + # the middle of the particle does not move along x + N_hal_x = int(Nx/2) + u_x_tmp = np.zeros((Nx+1,Ny+1), dtype=object) + u_y_tmp = np.zeros((Nx+1,Ny+1), dtype=object) + u_x_tmp[:N_hal_x,1:] = u_x[:N_hal_x,:] + u_x_tmp[N_hal_x+2:,1:] = u_x[N_hal_x+1:,:] + u_y_tmp[:N_hal_x,1:] = u_y[:N_hal_x,:] + u_y_tmp[N_hal_x+2:,1:] = u_y[N_hal_x+1,:] + + e1 = np.zeros((Nx,Ny), dtype=object) + e2 = np.zeros((Nx,Ny), dtype=object) duxdy = np.zeros((Nx,Ny), dtype=object) duydx = np.zeros((Nx,Ny), dtype=object) - # check boundaries !! + e1 = np.diff(u_x_tmp, axis = 1)/dxs + e2 = np.diff(u_y_tmp, axis = 0)/dys + # diff shape + duydx = np.diff(u_y_tmp, axis = 1)/dxs + duxdy = np.diff(u_x_tmp, axis = 0)/dys + + # for j in range(Ny+1): + # e1[:,j] = np.diff(u_x_tmp[:,j])/dxs + # duydx[:,j] = np.diff(u_y_tmp[:,j])/dxs + # for i in range(Nx+1): + # e2[i,:] = np.diff(u_y_tmp[i,:])/dys + # duxdy[i,:] = np.diff(u_x_tmp[i,:])/dys + e12 = np.zeros((Nx,Ny), dtype=object) for j in range(Ny): - e1_mat[:,j] = np.diff(u_x[:,j])/dxs - duydx[:,j] = np.diff(u_y[:,j])/dxs + e12[:,j] = 0.5*utils.mean_linear(duydx[:,j]) for i in range(Nx): - e2_mat[i,:] = np.diff(u_y[i,:])/dys - duxdy[i,:] = np.diff(u_x[i,:])/dys - e12_mat = 0.5*(duxdy+duydx) + e12[i,:] += 0.5*utils.mean_linear(duxdy[i,:]) - e_vec_of_ij = np.zeros(6, dtype=object) - sigma_vec_of_ij = np.zeros(6, dtype=object) + e_mat = np.zeros((Nx,Ny,6), dtype=object) + sigma_mat = np.zeros((Nx,Ny,6), dtype=object) - sigma_11_mat = np.zeros((Nx,Ny), dtype=object) - sigma_22_mat = np.zeros((Nx,Ny), dtype=object) - sigma_33_mat = np.zeros((Nx,Ny), dtype=object) - sigma_12_mat = np.zeros((Nx,Ny), dtype=object) - sigma_13_mat = np.zeros((Nx,Ny), dtype=object) - sigma_23_mat = np.zeros((Nx,Ny), dtype=object) - div_sigm_1_mat = np.zeros((Nx,Ny), dtype=object) - div_sigm_2_mat = np.zeros((Nx,Ny), dtype=object) - div_sigm_3_mat = np.zeros((Nx,Ny), dtype=object) for i in range(Nx): for j in range(Ny): - e_vec_of_ij = np.array([e1_mat[i,j]- e0[0]*c_mat[i,j], - e2_mat[i,j]- e0[1]*c_mat[i,j], + e_mat[i,j,:] = np.array([e1[i,j]- e0[0]*c_mat[i,j], + e2[i,j]- e0[1]*c_mat[i,j], 0, - e12_mat[i,j], + e12[i,j], 0, 0]) - sigma_vec_of_ij = np.dot((Cij,e_vec_of_ij)) - sigma_11_mat[i,j] = sigma_vec_of_ij[0] - sigma_22_mat[i,j] = sigma_vec_of_ij[1] - sigma_33_mat[i,j] = sigma_vec_of_ij[2] - sigma_12_mat[i,j] = sigma_vec_of_ij[3] - sigma_13_mat[i,j] = sigma_vec_of_ij[4] - sigma_23_mat[i,j] = sigma_vec_of_ij[5] - muR_el[i,j] = np.dot((sigma_vec_of_ij,e0))/max_conc + sigma_mat[i,j,:] = np.dot(Cij,e_mat[i,j,:]) + muR_el[i,j] = np.dot(sigma_mat[i,j,:],e0)/max_conc + + # now that ew found the chem pot + # it is time to create the div(sigma) so that out of the function + # can be posed = 0 + dsigma1dx_mat = np.zeros((Nx-1,Ny), dtype=object) + dsigma2dy_mat = np.zeros((Nx,Ny-1), dtype=object) + dsigma12dx_mat = np.zeros((Nx-1,Ny), dtype=object) + dsigma12dy_mat = np.zeros((Nx,Ny-1), dtype=object) + div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) + # check boundaries ! + for i in range(Nx): + dsigma2dy_mat[i,:] = np.diff(sigma_mat[i,:,1])/dys + dsigma12dy_mat[i,:] = np.diff(sigma_mat[i,:,2])/dys for j in range(Ny): - dsigma1dx = np.diff(sigma_11_mat[:,j])/dxs - dsgima12dx = np.diff(sigma_12_mat[:,j])/dxs - - + dsigma1dx_mat[:,j] = np.diff(sigma_mat[:,j,0])/dxs + dsigma12dx_mat[:,j] = np.diff(sigma_mat[:,j,2])/dxs + for i in range(int(Nx/2)): + for j in range(Ny-1): + div_stress_mat[i,j+1,:] = np.array([ + dsigma1dx_mat[i,j] + dsigma12dy_mat[i,j], + dsigma2dy_mat[i,j] + dsigma1dx_mat[i,j] + ] + ) + for i in range((int(Nx/2)+1),Nx,1): + for j in range(Ny-1): + div_stress_mat[i,j+1,:] = np.array([ + dsigma1dx_mat[i-1,j] + dsigma12dy_mat[i-1,j], + dsigma2dy_mat[i-1,j] + dsigma1dx_mat[i-1,j] + ] + ) + print(div_stress_mat) return muR_el, div_stress_mat diff --git a/mpet/sim.py b/mpet/sim.py index ac950878..518974c7 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -94,6 +94,8 @@ def SetUpVariables(self): part.c.SetInitialGuess(k, cs0) for j in range(N_ver_ij): part.cy[k].SetInitialCondition(j, cs0) + part.ux[k].SetInitialCondition(j, 0) + part.uy[k].SetInitialCondition(j, 0) else: for k in range(Nij): part.c.SetInitialCondition(k, cs0) From 87fb8bd2c6b16d4b77a6f9a4f5ad78ddde0f181f Mon Sep 17 00:00:00 2001 From: Ombrini Date: Mon, 12 Dec 2022 18:47:01 +0000 Subject: [PATCH 22/27] problem in the immition of 0 in the boundaries --- mpet/mod_electrodes.py | 49 +++++++++++++++++++++++++----------------- 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index c94b9093..944fab36 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -206,8 +206,8 @@ def sld_dynamics_2Dfull(self, c_mat, u_x_mat, u_y_mat, muO, act_lyte, noise): muR_mat += muR_el actR_mat = np.exp(muR_mat) - for i in range(Nx): - for j in range(Ny): + for i in range(Nx-1): + for j in range(Ny-1): eq1 = self.CreateEquation("divsigma1_{i}_{j}_equal0".format(i=i, j=j)) eq1.Residual = div_stress_mat[i,j,0] eq2 = self.CreateEquation("divsigma2_{i}_{j}_equal0".format(i=i, j=j)) @@ -873,35 +873,44 @@ def calc_muR_el(c_mat, u_x, u_y, config, trode, ind): sigma_mat[i,j,:] = np.dot(Cij,e_mat[i,j,:]) muR_el[i,j] = np.dot(sigma_mat[i,j,:],e0)/max_conc - # now that ew found the chem pot + # now that ew found the chem pot # it is time to create the div(sigma) so that out of the function # can be posed = 0 dsigma1dx_mat = np.zeros((Nx-1,Ny), dtype=object) dsigma2dy_mat = np.zeros((Nx,Ny-1), dtype=object) dsigma12dx_mat = np.zeros((Nx-1,Ny), dtype=object) dsigma12dy_mat = np.zeros((Nx,Ny-1), dtype=object) - div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) + # div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) + div_stress_mat = np.zeros((Nx-1,Ny-1,2), dtype=object) # check boundaries ! - for i in range(Nx): - dsigma2dy_mat[i,:] = np.diff(sigma_mat[i,:,1])/dys - dsigma12dy_mat[i,:] = np.diff(sigma_mat[i,:,2])/dys - for j in range(Ny): - dsigma1dx_mat[:,j] = np.diff(sigma_mat[:,j,0])/dxs - dsigma12dx_mat[:,j] = np.diff(sigma_mat[:,j,2])/dxs - for i in range(int(Nx/2)): + dsigma2dy_mat = np.diff(sigma_mat[:,:,1], axis = 1)/dys + dsigma12dy_mat = np.diff(sigma_mat[:,:,2], axis = 1)/dys + dsigma1dx_mat = np.diff(sigma_mat[:,:,0], axis = 0)/dxs + dsigma12dx_mat = np.diff(sigma_mat[:,:,2], axis = 0)/dxs + for i in range(Nx-1): for j in range(Ny-1): - div_stress_mat[i,j+1,:] = np.array([ + div_stress_mat[i,j,:] = np.array([ dsigma1dx_mat[i,j] + dsigma12dy_mat[i,j], - dsigma2dy_mat[i,j] + dsigma1dx_mat[i,j] - ] - ) - for i in range((int(Nx/2)+1),Nx,1): - for j in range(Ny-1): - div_stress_mat[i,j+1,:] = np.array([ - dsigma1dx_mat[i-1,j] + dsigma12dy_mat[i-1,j], - dsigma2dy_mat[i-1,j] + dsigma1dx_mat[i-1,j] + dsigma2dy_mat[i,j] + dsigma12dx_mat[i,j] ] ) + + # it does not like div(sigma) = np.zeros + # for i in range(int(Nx/2)): + # for j in range(Ny-1): + # div_stress_mat[i,j+1,:] = np.array([ + # dsigma1dx_mat[i,j] + dsigma12dy_mat[i,j], + # dsigma2dy_mat[i,j] + dsigma12dx_mat[i,j] + # ] + # ) + + # for i in range((int(Nx/2)+1),Nx,1): + # for j in range(Ny-1): + # div_stress_mat[i,j+1,:] = np.array([ + # dsigma1dx_mat[i-1,j] + dsigma12dy_mat[i-1,j], + # dsigma2dy_mat[i-1,j] + dsigma12dx_mat[i-1,j] + # ] + # ) print(div_stress_mat) return muR_el, div_stress_mat From c70a194dde6bfcd8a9e3e148e3e360da5bb17804 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Tue, 13 Dec 2022 11:35:38 +0000 Subject: [PATCH 23/27] the simulation starts but it fails, check numbers --- .devcontainer/Dockerfile | 16 ++++++ .devcontainer/devcontainer.json | 24 +++++++++ .devcontainer/noop.txt | 3 ++ mpet/mod_electrodes.py | 90 ++++++++++++++++++++------------- mpet/sim.py | 4 +- 5 files changed, 99 insertions(+), 38 deletions(-) create mode 100644 .devcontainer/Dockerfile create mode 100644 .devcontainer/devcontainer.json create mode 100644 .devcontainer/noop.txt diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile new file mode 100644 index 00000000..4b8cdce1 --- /dev/null +++ b/.devcontainer/Dockerfile @@ -0,0 +1,16 @@ +FROM --platform=linux/amd64 mcr.microsoft.com/devcontainers/miniconda:0-3 + +# Copy environment.yml (if found) to a temp location so we update the environment. Also +# copy "noop.txt" so the COPY instruction does not fail if no environment.yml exists. +COPY environment.yml* .devcontainer/noop.txt /tmp/conda-tmp/ +RUN if [ -f "/tmp/conda-tmp/environment.yml" ]; then umask 0002 && /opt/conda/bin/conda env update -n base -f /tmp/conda-tmp/environment.yml; fi \ + && rm -rf /tmp/conda-tmp + +# [Optional] Uncomment to install a different version of Python than the default +# RUN conda install -y python=3.6 \ +# && pip install --no-cache-dir pipx \ +# && pipx reinstall-all + +# [Optional] Uncomment this section to install additional OS packages. +# RUN apt-get update && export DEBIAN_FRONTEND=noninteractive \ +# && apt-get -y install --no-install-recommends diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 00000000..43d26c59 --- /dev/null +++ b/.devcontainer/devcontainer.json @@ -0,0 +1,24 @@ +// For format details, see https://aka.ms/devcontainer.json. For config options, see the +// README at: https://github.com/devcontainers/templates/tree/main/src/miniconda +{ + "name": "Miniconda (Python 3)", + "build": { + "context": "..", + "dockerfile": "Dockerfile" + } + + // Features to add to the dev container. More info: https://containers.dev/features. + // "features": {}, + + // Use 'forwardPorts' to make a list of ports inside the container available locally. + // "forwardPorts": [], + + // Use 'postCreateCommand' to run commands after the container is created. + // "postCreateCommand": "python --version", + + // Configure tool-specific properties. + // "customizations": {}, + + // Uncomment to connect as root instead. More info: https://aka.ms/dev-containers-non-root. + // "remoteUser": "root" +} diff --git a/.devcontainer/noop.txt b/.devcontainer/noop.txt new file mode 100644 index 00000000..abee1954 --- /dev/null +++ b/.devcontainer/noop.txt @@ -0,0 +1,3 @@ +This file is copied into the container along with environment.yml* from the +parent folder. This is done to prevent the Dockerfile COPY instruction from +failing if no environment.yml is found. \ No newline at end of file diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index 944fab36..dcdfccb4 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -206,8 +206,8 @@ def sld_dynamics_2Dfull(self, c_mat, u_x_mat, u_y_mat, muO, act_lyte, noise): muR_mat += muR_el actR_mat = np.exp(muR_mat) - for i in range(Nx-1): - for j in range(Ny-1): + for i in range(Nx): + for j in range(Ny): eq1 = self.CreateEquation("divsigma1_{i}_{j}_equal0".format(i=i, j=j)) eq1.Residual = div_stress_mat[i,j,0] eq2 = self.CreateEquation("divsigma2_{i}_{j}_equal0".format(i=i, j=j)) @@ -818,28 +818,46 @@ def mech_tensors(): return Cij, e0 -def calc_muR_el(c_mat, u_x, u_y, config, trode, ind): +def calc_muR_el(c_mat, u_x, u_y, conf, trode, ind): Ny = np.size(c_mat, 1) Nx = np.size(c_mat, 0) dys = 1./Ny dxs = 1./Nx - max_conc = config[trode, "rho_s"] - muR_el = np.zeros((Nx,Ny), dtype=object) + max_conc = conf[trode, "rho_s"] + T_ref = 298 + k = 1.381e-23 + N_A = 6.022e23 + kT = k * T_ref + norm_stress = kT * N_A * max_conc + muR_el = np.zeros((Nx+1,Ny+1), dtype=object) Cij, e0 = mech_tensors() + ywet = 0.98*np.ones(Ny+2, dtype=object) + c_mat_tmp = np.zeros((Nx+2,Ny+2), dtype=object) + c_mat_tmp[1:-1,1:-1] = c_mat + # first and last row is 0.98 + c_mat_tmp[-1,:] = ywet + c_mat_tmp[0,:] = ywet + c_mat_tmp[1:-1,0] = c_mat[:,0] + c_mat_tmp[1:-1,-1] = c_mat[:,-1] + # the middle of the particle does not move along x N_hal_x = int(Nx/2) - u_x_tmp = np.zeros((Nx+1,Ny+1), dtype=object) - u_y_tmp = np.zeros((Nx+1,Ny+1), dtype=object) - u_x_tmp[:N_hal_x,1:] = u_x[:N_hal_x,:] - u_x_tmp[N_hal_x+2:,1:] = u_x[N_hal_x+1:,:] - u_y_tmp[:N_hal_x,1:] = u_y[:N_hal_x,:] - u_y_tmp[N_hal_x+2:,1:] = u_y[N_hal_x+1,:] - - e1 = np.zeros((Nx,Ny), dtype=object) - e2 = np.zeros((Nx,Ny), dtype=object) - duxdy = np.zeros((Nx,Ny), dtype=object) - duydx = np.zeros((Nx,Ny), dtype=object) + u_x_tmp = np.zeros((Nx+2,Ny+2), dtype=object) + u_y_tmp = np.zeros((Nx+2,Ny+2), dtype=object) + u_x_tmp[2:,2:] = u_x + u_x_tmp[2:,0] = u_x[:,0] + u_x_tmp[2:,1] = u_x[:,0] + # u_x_tmp[N_hal_x+2:,1:] = u_x[N_hal_x+1:,:] + u_y_tmp[2:,2:] = u_y + u_y_tmp[0,2:] = u_y[0,:] + u_y_tmp[1,2:] = u_y[0,:] + # u_y_tmp[N_hal_x+2:,1:] = u_y[N_hal_x+1,:] + + e1 = np.zeros((Nx+1,Ny+1), dtype=object) + e2 = np.zeros((Nx+1,Ny+1), dtype=object) + duxdy = np.zeros((Nx+1,Ny+1), dtype=object) + duydx = np.zeros((Nx+1,Ny+1), dtype=object) e1 = np.diff(u_x_tmp, axis = 1)/dxs e2 = np.diff(u_y_tmp, axis = 0)/dys @@ -853,42 +871,42 @@ def calc_muR_el(c_mat, u_x, u_y, config, trode, ind): # for i in range(Nx+1): # e2[i,:] = np.diff(u_y_tmp[i,:])/dys # duxdy[i,:] = np.diff(u_x_tmp[i,:])/dys - e12 = np.zeros((Nx,Ny), dtype=object) - for j in range(Ny): + e12 = np.zeros((Nx+1,Ny+1), dtype=object) + for j in range(Ny+1): e12[:,j] = 0.5*utils.mean_linear(duydx[:,j]) - for i in range(Nx): + for i in range(Nx+1): e12[i,:] += 0.5*utils.mean_linear(duxdy[i,:]) - e_mat = np.zeros((Nx,Ny,6), dtype=object) - sigma_mat = np.zeros((Nx,Ny,6), dtype=object) - - for i in range(Nx): - for j in range(Ny): - e_mat[i,j,:] = np.array([e1[i,j]- e0[0]*c_mat[i,j], - e2[i,j]- e0[1]*c_mat[i,j], + e_mat = np.zeros((Nx+1,Ny+1,6), dtype=object) + sigma_mat = np.zeros((Nx+1,Ny+1,6), dtype=object) + + for i in range(Nx+1): + for j in range(Ny+1): + e_mat[i,j,:] = np.array([e1[i,j]- e0[0]*c_mat_tmp[i,j], + e2[i,j]- e0[1]*c_mat_tmp[i,j], 0, e12[i,j], 0, 0]) sigma_mat[i,j,:] = np.dot(Cij,e_mat[i,j,:]) - muR_el[i,j] = np.dot(sigma_mat[i,j,:],e0)/max_conc - + muR_el[i,j] = np.dot(sigma_mat[i,j,:],e0)/norm_stress + muR_el = muR_el[1:,1:] # now that ew found the chem pot # it is time to create the div(sigma) so that out of the function # can be posed = 0 - dsigma1dx_mat = np.zeros((Nx-1,Ny), dtype=object) - dsigma2dy_mat = np.zeros((Nx,Ny-1), dtype=object) - dsigma12dx_mat = np.zeros((Nx-1,Ny), dtype=object) - dsigma12dy_mat = np.zeros((Nx,Ny-1), dtype=object) + dsigma1dx_mat = np.zeros((Nx,Ny), dtype=object) + dsigma2dy_mat = np.zeros((Nx,Ny), dtype=object) + dsigma12dx_mat = np.zeros((Nx,Ny), dtype=object) + dsigma12dy_mat = np.zeros((Nx,Ny), dtype=object) # div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) - div_stress_mat = np.zeros((Nx-1,Ny-1,2), dtype=object) + div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) # check boundaries ! dsigma2dy_mat = np.diff(sigma_mat[:,:,1], axis = 1)/dys dsigma12dy_mat = np.diff(sigma_mat[:,:,2], axis = 1)/dys dsigma1dx_mat = np.diff(sigma_mat[:,:,0], axis = 0)/dxs dsigma12dx_mat = np.diff(sigma_mat[:,:,2], axis = 0)/dxs - for i in range(Nx-1): - for j in range(Ny-1): + for i in range(Nx): + for j in range(Ny): div_stress_mat[i,j,:] = np.array([ dsigma1dx_mat[i,j] + dsigma12dy_mat[i,j], dsigma2dy_mat[i,j] + dsigma12dx_mat[i,j] @@ -911,7 +929,7 @@ def calc_muR_el(c_mat, u_x, u_y, config, trode, ind): # dsigma2dy_mat[i-1,j] + dsigma12dx_mat[i-1,j] # ] # ) - print(div_stress_mat) + # print(div_stress_mat) return muR_el, div_stress_mat diff --git a/mpet/sim.py b/mpet/sim.py index 518974c7..f10d855b 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -94,8 +94,8 @@ def SetUpVariables(self): part.c.SetInitialGuess(k, cs0) for j in range(N_ver_ij): part.cy[k].SetInitialCondition(j, cs0) - part.ux[k].SetInitialCondition(j, 0) - part.uy[k].SetInitialCondition(j, 0) + part.ux[k].SetInitialGuess(j, 0) + part.uy[k].SetInitialGuess(j, 0) else: for k in range(Nij): part.c.SetInitialCondition(k, cs0) From a665f11d0d6e13433012ee425e50c22020763c3e Mon Sep 17 00:00:00 2001 From: Ombrini Date: Thu, 16 Mar 2023 18:29:28 +0100 Subject: [PATCH 24/27] 2D LFP simplified mechanics --- configs/params_LFP.cfg | 2 +- configs/params_LFP2D.cfg | 16 +- configs/params_system.cfg | 16 +- mpet/config/configuration.py | 8 + mpet/config/constants.py | 3 +- mpet/config/schemas.py | 4 + mpet/electrode/materials/LiFePO42D.py | 41 ++-- mpet/mod_electrodes.py | 337 +++++++++++++------------- mpet/sim.py | 4 +- 9 files changed, 233 insertions(+), 198 deletions(-) diff --git a/configs/params_LFP.cfg b/configs/params_LFP.cfg index 1dc1a982..0d0a8c5d 100644 --- a/configs/params_LFP.cfg +++ b/configs/params_LFP.cfg @@ -5,7 +5,7 @@ type = ACR discretization = 1e-9 shape = C3 -thickness = 20e-9 +thickness = 40e-9 [Material] muRfunc = LiFePO4 diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index a23a255e..8923c8ad 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -5,8 +5,9 @@ type = ACR2D discretization = 1e-9 shape = C3 -discretization_ver = 1e-9 -thickness = 2e-9 +discretization_ver = 10e-9 +# remember the it is half the thickness, check the normalizations +thickness = 100e-9 [Material] muRfunc = LiFePO42D @@ -15,18 +16,21 @@ noise_prefac = 1e-6 numnoise = 200 Omega_a = 1.8560e-20 # Omega_a = 0.5560e-20 -kappa = 5.0148e-10 -B = 0.1916e9 +kappa_x = 5.0148e-10 +kappa_y = 15e-10 +Bx = 0.1916e9 +By = 0.5e9 # B = 1e9 rho_s = 1.3793e28 -D = 5.3e-17 +D = 5e-15 Dfunc = lattice dgammadc = 0e-30 cwet = 0.98 [Reactions] +# surface_diff = true rxnType = BV -k0 = 1.6e-1 +k0 = 1.6e0 E_A = 13000 alpha = 0.5 # Fraggedakis et al. 2020, lambda = 8.3kBT diff --git a/configs/params_system.cfg b/configs/params_system.cfg index 00997d71..89fe8fb9 100644 --- a/configs/params_system.cfg +++ b/configs/params_system.cfg @@ -8,12 +8,12 @@ profileType = CC # Battery (dis)charge c-rate (only used for CC), number of capacities / hr # (positive for discharge, negative for charge) -Crate = 1 +Crate = 5 # Optional nominal 1C current density for the cell, A/m^2 # 1C_current_density = 12.705 # Voltage cutoffs, V Vmax = 3.6 -Vmin = 3.3 +Vmin = 3.0 # Battery applied voltage (only used for CV), V Vset = 0.12 # Battery constant power (only used for CP), W/m^2 @@ -25,9 +25,9 @@ power = 1 # This can have as many rows (segments) as desired. # Note: It's okay to leave commented lines within the segments list segments = [ - (0.3, 0.4), - # (0, 0.2), - (-0.5, 0.1), + (10, 3), + (0, 30), + (5, 10), ] # Continuation directory. If false, begin a fresh simulation with the # specified input parameters here. Otherwise, this should be the @@ -55,7 +55,7 @@ T = 298 # (affects noise, particle size distribution). Set to true exactly # reproducible results -- useful for testing. # Options: true, false -randomSeed = false +randomSeed = true # Value of the random seed, must be an integer seed = 0 # Data reporter: choice of mat (MATLAB), hdf5 (hdf5), or hdf5Fast (hdf5, without @@ -95,8 +95,8 @@ Rfilm_foil = 0e-0 # sphere or cylinder -- radius # If using stddev = 0, set Npart = 1 for that electrode to avoid # wasting computational effort. -mean_c = 4e-9 -stddev_c = 10e-9 +mean_c = 50e-9 +stddev_c = 0e-9 mean_a = 100e-9 stddev_a = 1e-9 # Use a specific set of particle sizes for the distribution diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 1da05a7b..6254b371 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -500,6 +500,10 @@ def _scale_electrode_parameters(self): self[trode, 'lambda'] = self[trode, 'lambda'] / kT if self[trode, 'B'] is not None: self[trode, 'B'] = self[trode, 'B'] / (kT * constants.N_A * self[trode, 'cs_ref']) + if self[trode, 'Bx'] is not None: + self[trode,'Bx'] = self[trode, 'Bx'] / (kT * constants.N_A * self[trode, 'cs_ref']) + if self[trode, 'By'] is not None: + self[trode,'By'] = self[trode, 'By'] / (kT * constants.N_A * self[trode, 'cs_ref']) for param in ['Omega_a', 'Omega_b', 'Omega_c', 'EvdW']: value = self[trode, param] if value is not None: @@ -752,6 +756,10 @@ def _indvPart(self): self[trode, 'indvPart']['N_ver'][i,j] = self['psd_num_ver'][trode][i,j] if self[trode, 'kappa'] is not None: self[trode, 'indvPart']['kappa'][i, j] = self[trode, 'kappa'] / kappa_ref + if self[trode, 'kappa_x'] is not None: + self[trode,'indvPart']['kappa_x'][i, j] = self[trode,'kappa_x'] / kappa_ref + if self[trode, 'kappa_y'] is not None: + self[trode,'indvPart']['kappa_y'][i, j] = self[trode,'kappa_y'] / kappa_ref if self[trode, 'dgammadc'] is not None: nd_dgammadc = self[trode, 'dgammadc'] * cs_ref_part / gamma_S_ref self[trode, 'indvPart']['beta_s'][i, j] = nd_dgammadc \ diff --git a/mpet/config/constants.py b/mpet/config/constants.py index 375b23b4..ebdc5f51 100644 --- a/mpet/config/constants.py +++ b/mpet/config/constants.py @@ -26,6 +26,7 @@ #: subset of ``PARAMS_PER_TRODE``` that is defined for the separator as well PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp'] #: parameters that are defined for each particle, and their type -PARAMS_PARTICLE = {'N': int, 'N_ver':int, 'kappa': float, 'beta_s': float, 'D': float, 'k0': float, +PARAMS_PARTICLE = {'N': int, 'N_ver':int, 'kappa': float, 'kappa_x': float, 'kappa_y': float, + 'beta_s': float, 'D': float, 'k0': float, 'Rfilm': float, 'delta_L': float, 'Omega_a': float, 'E_D': float, 'E_A': float} diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py index 142cc868..1e66aa8d 100644 --- a/mpet/config/schemas.py +++ b/mpet/config/schemas.py @@ -144,7 +144,11 @@ def tobool(value): Optional('Omega_b', default=None): Use(float), Optional('Omega_c', default=None): Use(float), Optional('kappa', default=None): Use(float), + Optional('kappa_x', default=None): Use(float), + Optional('kappa_y', default=None): Use(float), Optional('B', default=None): Use(float), + Optional('Bx', default=None): Use(float), + Optional('By', default=None): Use(float), Optional('EvdW', default=None): Use(float), 'rho_s': Use(float), 'D': Use(float), diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py index 9da50bc8..bd15406d 100644 --- a/mpet/electrode/materials/LiFePO42D.py +++ b/mpet/electrode/materials/LiFePO42D.py @@ -4,8 +4,8 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): muRtheta = -self.eokT*3.422 Omega = self.get_trode_param("Omega_a") - kappax = self.get_trode_param("kappa") - kappay = kappax*10 + kappax = self.get_trode_param("kappa_x") + kappay = self.get_trode_param("kappa_y") # beta_s = self.get_trode_param("beta_s") # for allowing muR_ref to be calculated if np.size(c_mat) == 1: @@ -30,24 +30,33 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): # first and last row is 0.98 c_mat_tmp[-1,:] = ywet c_mat_tmp[0,:] = ywet - curvy = np.empty(Ny, dtype=object) - # print(muR_mat) - # plan: produce the matrix muR without curvy - # then another for loop for the curvy - # if then works try to optmize the loops - for i in range(1,Nx+1): - curvy[0] = 2*c_mat_tmp[i,1] - 2*c_mat_tmp[i,0] - curvy[-1] = 2*c_mat_tmp[i,-2] - 2*c_mat_tmp[i,-1] - # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) - curvy[1:Ny-1] = np.diff(c_mat_tmp[i,:],2) - for j in range(Ny): + + # curvy = np.empty(Ny, dtype=object) + + curvy = np.empty_like(c_mat_tmp) + + curvy[:,0] = 2 * c_mat_tmp[:,1] - 2 * c_mat_tmp[:,0] + curvy[:,-1] = 2 * c_mat_tmp[:,-2] - 2 * c_mat_tmp[:,-1] + curvy[:,1:-1] = np.diff(c_mat_tmp, n=2, axis=1) + + curvx = np.diff(c_mat_tmp, n=2, axis=0) + # muRhomog = self.T * np.log(c_mat_tmp / (1 - c_mat_tmp)) + Omega * (1 - 2 * c_mat) + y_vert_avg = np.average(c_mat_tmp, axis=1) + y_oriz_avg = np.average(c_mat_tmp, axis=0) + for i in np.arange(1,Nx+1): + # curvy[0] = 2*c_mat_tmp[i,1] - 2*c_mat_tmp[i,0] + # curvy[-1] = 2*c_mat_tmp[i,-2] - 2*c_mat_tmp[i,-1] + # # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) + # curvy[1:Ny-1] = np.diff(c_mat_tmp[i,:],2) + for j in np.arange(Ny): y = c_mat_tmp[i,j] muRhomog = self.T*np.log(y/(1-y)) + Omega*(1-2*y) muR_mat[i-1,j] = muRhomog - curvx = (c_mat_tmp[i-1,j] - c_mat_tmp[i,j]) - (c_mat_tmp[i,j] - c_mat_tmp[i+1,j]) - curv_en = -kappax*curvx/(dxs**2) - kappay*curvy[j]/(dys**2) + # curvx = (c_mat_tmp[i-1,j] - c_mat_tmp[i,j]) - (c_mat_tmp[i,j] - c_mat_tmp[i+1,j]) + curv_en = -kappax*curvx[i-1,j]/(dxs**2) - kappay*curvy[i,j]/(dys**2) muR_mat[i-1,j] += curv_en - muR_mat[i-1,j] += self.get_trode_param("B")*(y - ybar) + muR_mat[i-1,j] += self.get_trode_param("Bx")*(y - y_oriz_avg[j]) + muR_mat[i-1,j] += self.get_trode_param("By")*(y - y_vert_avg[i]) # print(muR_mat) actR_mat = np.exp(muR_mat/self.T) muR_mat += muRtheta + muR_ref diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index dcdfccb4..f432f697 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -46,20 +46,25 @@ def __init__(self, config, trode, vInd, pInd, Nx = self.get_trode_param("N") self.cy = {} for k in range(Nx): - self.cy[k] = dae.daeVariable("cy{k}".format(k=k), mole_frac_t, self, - "Concentration in y direction of element in row {k}".format(k=k), - [self.Dmny]) - # check unit of measures - self.uy = {} - for k in range(Nx): - self.uy[k] = dae.daeVariable("uy{k}".format(k=k), mole_frac_t, self, - "Displacement in y direction of element in row {k}".format(k=k), - [self.Dmny]) - self.ux = {} - for k in range(Nx): - self.ux[k] = dae.daeVariable("ux{k}".format(k=k), mole_frac_t, self, - "Displacement in x direction of element in row {k}".format(k=k), + self.cy[k] = dae.daeVariable("cy{k}".format(k=k), + mole_frac_t, self, + "Conc in ver direction of element in row {k}".format(k=k), [self.Dmny]) + # # check unit of measures + # self.uy = {} + # for k in range(Nx): + # self.uy[k] = dae.daeVariable("uy{k}".format(k=k), mole_frac_t, + # self, + # "Displacement in y direction of element in row {k}". + # format(k=k), + # [self.Dmny]) + # self.ux = {} + # for k in range(Nx): + # self.ux[k] = dae.daeVariable("ux{k}".format(k=k), mole_frac_t, + # self, + # "Displacement in x direction of element in row {k}". + # format(k=k), + # [self.Dmny]) self.cbar = dae.daeVariable( "cbar", mole_frac_t, self, @@ -128,16 +133,17 @@ def DeclareEquations(self): for k in range(Nx): c_mat[k,:] = [self.cy[k](j) for j in range(Ny)] - u_y_mat = np.empty((Nx, Ny), dtype=object) - for k in range(Nx): - u_y_mat[k,:] = [self.uy[k](j) for j in range(Ny)] + # u_y_mat = np.empty((Nx, Ny), dtype=object) + # for k in range(Nx): + # u_y_mat[k,:] = [self.uy[k](j) for j in range(Ny)] - u_x_mat = np.empty((Nx, Ny), dtype=object) - for k in range(Nx): - u_x_mat[k,:] = [self.ux[k](j) for j in range(Ny)] + # u_x_mat = np.empty((Nx, Ny), dtype=object) + # for k in range(Nx): + # u_x_mat[k,:] = [self.ux[k](j) for j in range(Ny)] # self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) - self.sld_dynamics_2Dfull(c_mat, u_x_mat, u_y_mat, mu_O, act_lyte, self.noise) + # self.sld_dynamics_2Dfull(c_mat, u_x_mat, u_y_mat, mu_O, act_lyte, self.noise) + self.sld_dynamics_2D1var(c_mat, mu_O, act_lyte, self.noise) for eq in self.Equations: eq.CheckUnitsConsistency = False @@ -152,11 +158,15 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): f"mpet.electrode.diffusion.{Dfunc_name}") dr, edges = geo.get_dr_edges("C3", Ny) area_vec = 1. - Mmaty = get_Mmat("C3", Ny) + Mmaty = get_Mmat(self.get_trode_param('shape'), Ny) + muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), + self.config, self.trode, self.ind) for k in range(Nx): c_vec = c_mat[k,:] - muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), - self.config, self.trode, self.ind) + muR_vec = muR_mat[k,:] + actR_vec = actR_mat[k,:] + # muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), + # self.config, self.trode, self.ind) c_surf = c_mat[k,-1] muR_surf = muR_vec[-1] actR_surf = actR_vec[-1] @@ -174,8 +184,8 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): Flux_bc = -self.Rxn(k) - Flux_vec = calc_flux_C3ver(c_vec, muR_vec, self.get_trode_param("D"), Dfunc, - self.get_trode_param("E_D"), Flux_bc, dr, T, noise) + Flux_vec = calc_flux_CHR(c_vec, muR_vec, self.get_trode_param("D"), Dfunc, + self.get_trode_param("E_D"), Flux_bc, dr, T, noise) RHS_vec = -np.diff(Flux_vec * area_vec) dcdt_vec_y = np.empty(Ny, dtype=object) @@ -200,18 +210,18 @@ def sld_dynamics_2Dfull(self, c_mat, u_x_mat, u_y_mat, muO, act_lyte, noise): # print(c_mat) muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), self.config, self.trode, self.ind) - muR_el, div_stress_mat = calc_muR_el(c_mat, u_x_mat, u_y_mat, - self.config, self.trode, self.ind) + # muR_el, div_stress_mat = calc_muR_el(c_mat, u_x_mat, u_y_mat, + # self.config, self.trode, self.ind) - muR_mat += muR_el + # muR_mat += muR_el actR_mat = np.exp(muR_mat) - for i in range(Nx): - for j in range(Ny): - eq1 = self.CreateEquation("divsigma1_{i}_{j}_equal0".format(i=i, j=j)) - eq1.Residual = div_stress_mat[i,j,0] - eq2 = self.CreateEquation("divsigma2_{i}_{j}_equal0".format(i=i, j=j)) - eq2.Residual = div_stress_mat[i,j,1] + # for i in range(Nx): + # for j in range(Ny): + # eq1 = self.CreateEquation("divsigma1_{i}_{j}_equal0".format(i=i, j=j)) + # eq1.Residual = div_stress_mat[i,j,0] + # eq2 = self.CreateEquation("divsigma2_{i}_{j}_equal0".format(i=i, j=j)) + # eq2.Residual = div_stress_mat[i,j,1] for k in range(Nx): c_vec = c_mat[k,:] @@ -243,7 +253,7 @@ def sld_dynamics_2Dfull(self, c_mat, u_x_mat, u_y_mat, muO, act_lyte, noise): Flux_vec = calc_flux_C3ver(c_vec, muR_vec, self.get_trode_param("D"), Dfunc, self.get_trode_param("E_D"), Flux_bc, dr, T, noise) - + area_vec = 1.*edges RHS_vec = -np.diff(Flux_vec * area_vec) dcdt_vec_y = np.empty(Ny, dtype=object) dcdt_vec_y[0:Ny] = [self.cy[k].dt(j) for j in range(Ny)] @@ -617,7 +627,6 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise): # Mass matrix, M, where M*dcdt = RHS, where c and RHS are vectors Mmat = get_Mmat(self.get_trode_param('shape'), N) dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), N) - # Get solid particle chemical potential, overpotential, reaction rate if self.get_trode_param("type") in ["ACR"]: c_surf = c @@ -785,133 +794,133 @@ def calc_muR(c, cbar, config, trode, ind): muR, actR = muRfunc(c, cbar, muR_ref) return muR, actR -def mech_tensors(): - # FePo4 elastic constants (GPa) - c11 = 175.9 - c22 = 153.6 - c33 = 135.0 - c44 = 38.8 - c55 = 47.5 - c66 = 55.6 - c13 = 54.0 - c12 = 29.6 - c23 = 19.6 - - Cij = np.zeros((6,6)) - Cij[0,0] = c11 - Cij[1,1] = c22 - Cij[2,2] = c33 - Cij[3,3] = c44 - Cij[4,4] = c55 - Cij[5,5] = c66 - Cij[1,0] = c12 - Cij[0,1] = c12 - Cij[2,0] = c13 - Cij[0,2] = c13 - Cij[1,2] = c23 - Cij[2,1] = c23 - # strain - e01 = 0.0517 - e02 = 0.0359 - e03 = -0.0186 - e0 = np.array([e01, e02, e03, 0, 0, 0]) - - return Cij, e0 - -def calc_muR_el(c_mat, u_x, u_y, conf, trode, ind): - Ny = np.size(c_mat, 1) - Nx = np.size(c_mat, 0) - dys = 1./Ny - dxs = 1./Nx - max_conc = conf[trode, "rho_s"] - T_ref = 298 - k = 1.381e-23 - N_A = 6.022e23 - kT = k * T_ref - norm_stress = kT * N_A * max_conc - muR_el = np.zeros((Nx+1,Ny+1), dtype=object) - Cij, e0 = mech_tensors() - - ywet = 0.98*np.ones(Ny+2, dtype=object) - c_mat_tmp = np.zeros((Nx+2,Ny+2), dtype=object) - c_mat_tmp[1:-1,1:-1] = c_mat - # first and last row is 0.98 - c_mat_tmp[-1,:] = ywet - c_mat_tmp[0,:] = ywet - c_mat_tmp[1:-1,0] = c_mat[:,0] - c_mat_tmp[1:-1,-1] = c_mat[:,-1] - - # the middle of the particle does not move along x - N_hal_x = int(Nx/2) - u_x_tmp = np.zeros((Nx+2,Ny+2), dtype=object) - u_y_tmp = np.zeros((Nx+2,Ny+2), dtype=object) - u_x_tmp[2:,2:] = u_x - u_x_tmp[2:,0] = u_x[:,0] - u_x_tmp[2:,1] = u_x[:,0] - # u_x_tmp[N_hal_x+2:,1:] = u_x[N_hal_x+1:,:] - u_y_tmp[2:,2:] = u_y - u_y_tmp[0,2:] = u_y[0,:] - u_y_tmp[1,2:] = u_y[0,:] - # u_y_tmp[N_hal_x+2:,1:] = u_y[N_hal_x+1,:] - - e1 = np.zeros((Nx+1,Ny+1), dtype=object) - e2 = np.zeros((Nx+1,Ny+1), dtype=object) - duxdy = np.zeros((Nx+1,Ny+1), dtype=object) - duydx = np.zeros((Nx+1,Ny+1), dtype=object) - - e1 = np.diff(u_x_tmp, axis = 1)/dxs - e2 = np.diff(u_y_tmp, axis = 0)/dys - # diff shape - duydx = np.diff(u_y_tmp, axis = 1)/dxs - duxdy = np.diff(u_x_tmp, axis = 0)/dys - - # for j in range(Ny+1): - # e1[:,j] = np.diff(u_x_tmp[:,j])/dxs - # duydx[:,j] = np.diff(u_y_tmp[:,j])/dxs - # for i in range(Nx+1): - # e2[i,:] = np.diff(u_y_tmp[i,:])/dys - # duxdy[i,:] = np.diff(u_x_tmp[i,:])/dys - e12 = np.zeros((Nx+1,Ny+1), dtype=object) - for j in range(Ny+1): - e12[:,j] = 0.5*utils.mean_linear(duydx[:,j]) - for i in range(Nx+1): - e12[i,:] += 0.5*utils.mean_linear(duxdy[i,:]) - - e_mat = np.zeros((Nx+1,Ny+1,6), dtype=object) - sigma_mat = np.zeros((Nx+1,Ny+1,6), dtype=object) - - for i in range(Nx+1): - for j in range(Ny+1): - e_mat[i,j,:] = np.array([e1[i,j]- e0[0]*c_mat_tmp[i,j], - e2[i,j]- e0[1]*c_mat_tmp[i,j], - 0, - e12[i,j], - 0, - 0]) - sigma_mat[i,j,:] = np.dot(Cij,e_mat[i,j,:]) - muR_el[i,j] = np.dot(sigma_mat[i,j,:],e0)/norm_stress - muR_el = muR_el[1:,1:] - # now that ew found the chem pot - # it is time to create the div(sigma) so that out of the function - # can be posed = 0 - dsigma1dx_mat = np.zeros((Nx,Ny), dtype=object) - dsigma2dy_mat = np.zeros((Nx,Ny), dtype=object) - dsigma12dx_mat = np.zeros((Nx,Ny), dtype=object) - dsigma12dy_mat = np.zeros((Nx,Ny), dtype=object) - # div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) - div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) - # check boundaries ! - dsigma2dy_mat = np.diff(sigma_mat[:,:,1], axis = 1)/dys - dsigma12dy_mat = np.diff(sigma_mat[:,:,2], axis = 1)/dys - dsigma1dx_mat = np.diff(sigma_mat[:,:,0], axis = 0)/dxs - dsigma12dx_mat = np.diff(sigma_mat[:,:,2], axis = 0)/dxs - for i in range(Nx): - for j in range(Ny): - div_stress_mat[i,j,:] = np.array([ - dsigma1dx_mat[i,j] + dsigma12dy_mat[i,j], - dsigma2dy_mat[i,j] + dsigma12dx_mat[i,j] - ] - ) +# def mech_tensors(): +# # FePo4 elastic constants (GPa) +# c11 = 175.9 +# c22 = 153.6 +# c33 = 135.0 +# c44 = 38.8 +# c55 = 47.5 +# c66 = 55.6 +# c13 = 54.0 +# c12 = 29.6 +# c23 = 19.6 + +# Cij = np.zeros((6,6)) +# Cij[0,0] = c11 +# Cij[1,1] = c22 +# Cij[2,2] = c33 +# Cij[3,3] = c44 +# Cij[4,4] = c55 +# Cij[5,5] = c66 +# Cij[1,0] = c12 +# Cij[0,1] = c12 +# Cij[2,0] = c13 +# Cij[0,2] = c13 +# Cij[1,2] = c23 +# Cij[2,1] = c23 +# # strain +# e01 = 0.0517 +# e02 = 0.0359 +# e03 = -0.0186 +# e0 = np.array([e01, e02, e03, 0, 0, 0]) + +# return Cij, e0 + +# def calc_muR_el(c_mat, u_x, u_y, conf, trode, ind): +# Ny = np.size(c_mat, 1) +# Nx = np.size(c_mat, 0) +# dys = 1./Ny +# dxs = 1./Nx +# max_conc = conf[trode, "rho_s"] +# T_ref = 298 +# k = 1.381e-23 +# N_A = 6.022e23 +# kT = k * T_ref +# norm_stress = kT * N_A * max_conc +# muR_el = np.zeros((Nx+1,Ny+1), dtype=object) +# Cij, e0 = mech_tensors() + +# ywet = 0.98*np.ones(Ny+2, dtype=object) +# c_mat_tmp = np.zeros((Nx+2,Ny+2), dtype=object) +# c_mat_tmp[1:-1,1:-1] = c_mat +# # first and last row is 0.98 +# c_mat_tmp[-1,:] = ywet +# c_mat_tmp[0,:] = ywet +# c_mat_tmp[1:-1,0] = c_mat[:,0] +# c_mat_tmp[1:-1,-1] = c_mat[:,-1] + +# # the middle of the particle does not move along x +# N_hal_x = int(Nx/2) +# u_x_tmp = np.zeros((Nx+2,Ny+2), dtype=object) +# u_y_tmp = np.zeros((Nx+2,Ny+2), dtype=object) +# u_x_tmp[2:,2:] = u_x +# u_x_tmp[2:,0] = u_x[:,0] +# u_x_tmp[2:,1] = u_x[:,0] +# # u_x_tmp[N_hal_x+2:,1:] = u_x[N_hal_x+1:,:] +# u_y_tmp[2:,2:] = u_y +# u_y_tmp[0,2:] = u_y[0,:] +# u_y_tmp[1,2:] = u_y[0,:] +# # u_y_tmp[N_hal_x+2:,1:] = u_y[N_hal_x+1,:] + +# e1 = np.zeros((Nx+1,Ny+1), dtype=object) +# e2 = np.zeros((Nx+1,Ny+1), dtype=object) +# duxdy = np.zeros((Nx+1,Ny+1), dtype=object) +# duydx = np.zeros((Nx+1,Ny+1), dtype=object) + +# e1 = np.diff(u_x_tmp, axis = 1)/dxs +# e2 = np.diff(u_y_tmp, axis = 0)/dys +# # diff shape +# duydx = np.diff(u_y_tmp, axis = 1)/dxs +# duxdy = np.diff(u_x_tmp, axis = 0)/dys + +# # for j in range(Ny+1): +# # e1[:,j] = np.diff(u_x_tmp[:,j])/dxs +# # duydx[:,j] = np.diff(u_y_tmp[:,j])/dxs +# # for i in range(Nx+1): +# # e2[i,:] = np.diff(u_y_tmp[i,:])/dys +# # duxdy[i,:] = np.diff(u_x_tmp[i,:])/dys +# e12 = np.zeros((Nx+1,Ny+1), dtype=object) +# for j in range(Ny+1): +# e12[:,j] = 0.5*utils.mean_linear(duydx[:,j]) +# for i in range(Nx+1): +# e12[i,:] += 0.5*utils.mean_linear(duxdy[i,:]) + +# e_mat = np.zeros((Nx+1,Ny+1,6), dtype=object) +# sigma_mat = np.zeros((Nx+1,Ny+1,6), dtype=object) + +# for i in range(Nx+1): +# for j in range(Ny+1): +# e_mat[i,j,:] = np.array([e1[i,j]- e0[0]*c_mat_tmp[i,j], +# e2[i,j]- e0[1]*c_mat_tmp[i,j], +# 0, +# e12[i,j], +# 0, +# 0]) +# sigma_mat[i,j,:] = np.dot(Cij,e_mat[i,j,:]) +# muR_el[i,j] = np.dot(sigma_mat[i,j,:],e0)/norm_stress +# muR_el = muR_el[1:,1:] +# # now that ew found the chem pot +# # it is time to create the div(sigma) so that out of the function +# # can be posed = 0 +# dsigma1dx_mat = np.zeros((Nx,Ny), dtype=object) +# dsigma2dy_mat = np.zeros((Nx,Ny), dtype=object) +# dsigma12dx_mat = np.zeros((Nx,Ny), dtype=object) +# dsigma12dy_mat = np.zeros((Nx,Ny), dtype=object) +# # div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) +# div_stress_mat = np.zeros((Nx,Ny,2), dtype=object) +# # check boundaries ! +# dsigma2dy_mat = np.diff(sigma_mat[:,:,1], axis = 1)/dys +# dsigma12dy_mat = np.diff(sigma_mat[:,:,2], axis = 1)/dys +# dsigma1dx_mat = np.diff(sigma_mat[:,:,0], axis = 0)/dxs +# dsigma12dx_mat = np.diff(sigma_mat[:,:,2], axis = 0)/dxs +# for i in range(Nx): +# for j in range(Ny): +# div_stress_mat[i,j,:] = np.array([ +# dsigma1dx_mat[i,j] + dsigma12dy_mat[i,j], +# dsigma2dy_mat[i,j] + dsigma12dx_mat[i,j] +# ] +# ) # it does not like div(sigma) = np.zeros # for i in range(int(Nx/2)): @@ -930,7 +939,7 @@ def calc_muR_el(c_mat, u_x, u_y, conf, trode, ind): # ] # ) # print(div_stress_mat) - return muR_el, div_stress_mat + # return muR_el, div_stress_mat def MX(mat, objvec): diff --git a/mpet/sim.py b/mpet/sim.py index f10d855b..42183d7c 100644 --- a/mpet/sim.py +++ b/mpet/sim.py @@ -94,8 +94,8 @@ def SetUpVariables(self): part.c.SetInitialGuess(k, cs0) for j in range(N_ver_ij): part.cy[k].SetInitialCondition(j, cs0) - part.ux[k].SetInitialGuess(j, 0) - part.uy[k].SetInitialGuess(j, 0) + # part.ux[k].SetInitialGuess(j, 0) + # part.uy[k].SetInitialGuess(j, 0) else: for k in range(Nij): part.c.SetInitialCondition(k, cs0) From e0b638abfe0b409e9dbe9971fac07bbedc712e45 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Fri, 17 Mar 2023 10:19:02 +0100 Subject: [PATCH 25/27] fix normalation rxn 2d lfp --- mpet/mod_electrodes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index f432f697..fa15e226 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -182,7 +182,7 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq = self.CreateEquation("Rxn_{k}".format(k=k)) eq.Residual = self.Rxn(k) - Rxn - Flux_bc = -self.Rxn(k) + Flux_bc = -self.Rxn(k)*self.get_trode_param("delta_L") Flux_vec = calc_flux_CHR(c_vec, muR_vec, self.get_trode_param("D"), Dfunc, self.get_trode_param("E_D"), Flux_bc, dr, T, noise) From 616a6e680619919fb5cff09ade41b42bc84480a6 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Fri, 17 Mar 2023 13:52:18 +0100 Subject: [PATCH 26/27] surface diffusion --- configs/params_LFP2D.cfg | 13 ++++++++----- mpet/config/configuration.py | 5 +++++ mpet/config/constants.py | 4 ++-- mpet/config/schemas.py | 3 +++ mpet/mod_electrodes.py | 32 +++++++++++++++++++++++++++++--- 5 files changed, 47 insertions(+), 10 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index 8923c8ad..d309fd20 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -7,7 +7,7 @@ discretization = 1e-9 shape = C3 discretization_ver = 10e-9 # remember the it is half the thickness, check the normalizations -thickness = 100e-9 +thickness = 50e-9 [Material] muRfunc = LiFePO42D @@ -22,15 +22,18 @@ Bx = 0.1916e9 By = 0.5e9 # B = 1e9 rho_s = 1.3793e28 -D = 5e-15 +D = 1e-14 +E_D = 0 Dfunc = lattice dgammadc = 0e-30 cwet = 0.98 +D_surf = 1e-16 +E_D_surf = 0 [Reactions] -# surface_diff = true -rxnType = BV -k0 = 1.6e0 +surface_diffusion = false +rxnType = CIET +k0 = 30 E_A = 13000 alpha = 0.5 # Fraggedakis et al. 2020, lambda = 8.3kBT diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index 6254b371..ec830d64 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -768,11 +768,16 @@ def _indvPart(self): pthick = self['psd_len_ver'][trode][i,j] self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] \ * self['t_ref'] / pthick**2 + self[trode,'indvPart']['D_surf'][i, j] = self[trode,'D_surf'] \ + * self['t_ref'] / plen**2 else: self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] \ * self['t_ref'] / plen**2 self[trode, 'indvPart']['E_D'][i, j] = self[trode, 'E_D'] \ / (constants.k * constants.N_A * constants.T_ref) + if self[trode, 'E_D_surf'] is not None: + self[trode, 'indvPart']['E_D_surf'][i, j] = self[trode, 'E_D_surf'] \ + / (constants.k * constants.N_A * constants.T_ref) self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \ / (constants.e * F_s_ref) self[trode, 'indvPart']['E_A'][i, j] = self[trode, 'E_A'] \ diff --git a/mpet/config/constants.py b/mpet/config/constants.py index ebdc5f51..126d2287 100644 --- a/mpet/config/constants.py +++ b/mpet/config/constants.py @@ -27,6 +27,6 @@ PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp'] #: parameters that are defined for each particle, and their type PARAMS_PARTICLE = {'N': int, 'N_ver':int, 'kappa': float, 'kappa_x': float, 'kappa_y': float, - 'beta_s': float, 'D': float, 'k0': float, + 'beta_s': float, 'D': float,'D_surf': float, 'k0': float, 'Rfilm': float, 'delta_L': float, 'Omega_a': float, 'E_D': float, - 'E_A': float} + 'E_D_surf': float, 'E_A': float} diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py index 1e66aa8d..97594025 100644 --- a/mpet/config/schemas.py +++ b/mpet/config/schemas.py @@ -152,15 +152,18 @@ def tobool(value): Optional('EvdW', default=None): Use(float), 'rho_s': Use(float), 'D': Use(float), + Optional('D_surf', default=0.): Use(float), Optional('Dfunc_filename', default=None): str, 'Dfunc': str, Optional('E_D', default=0.): Use(float), + Optional('E_D_surf', default=0.): Use(float), Optional('dgammadc', default=None): Use(float), Optional('cwet', default=None): Use(float)}, 'Reactions': {Optional('rxnType_filename', default=None): str, 'rxnType': str, 'k0': Use(float), Optional('E_A', default=0.): Use(float), + Optional('surface_diffusion', default=None): Use(tobool), 'alpha': Use(float), Optional('lambda', default=None): Use(float), 'Rfilm': Use(float)}} diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index fa15e226..e61574ad 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -161,13 +161,19 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): Mmaty = get_Mmat(self.get_trode_param('shape'), Ny) muR_mat, actR_mat = calc_muR(c_mat, self.cbar(), self.config, self.trode, self.ind) + + if self.get_trode_param("surface_diffusion"): + surf_diff_vec = calc_surf_diff(c_mat[:,-1], muR_mat[:,-1], + self.get_trode_param("cwet"), + self.get_trode_param("D_surf"), self.config["T"], + self.get_trode_param("E_D_surf")) for k in range(Nx): c_vec = c_mat[k,:] muR_vec = muR_mat[k,:] actR_vec = actR_mat[k,:] # muR_vec, actR_vec = calc_muR(c_vec, self.cbar(), # self.config, self.trode, self.ind) - c_surf = c_mat[k,-1] + c_surf = c_vec[-1] muR_surf = muR_vec[-1] actR_surf = actR_vec[-1] eta = calc_eta(muR_surf, muO) @@ -182,7 +188,10 @@ def sld_dynamics_2D1var(self, c_mat, muO, act_lyte, noise): eq = self.CreateEquation("Rxn_{k}".format(k=k)) eq.Residual = self.Rxn(k) - Rxn - Flux_bc = -self.Rxn(k)*self.get_trode_param("delta_L") + if self.get_trode_param("surface_diffusion"): + Flux_bc = -self.Rxn(k)*self.get_trode_param("delta_L") + surf_diff_vec[k] + else: + Flux_bc = -self.Rxn(k)*self.get_trode_param("delta_L") Flux_vec = calc_flux_CHR(c_vec, muR_vec, self.get_trode_param("D"), Dfunc, self.get_trode_param("E_D"), Flux_bc, dr, T, noise) @@ -714,6 +723,23 @@ def get_Mmat(shape, N): return Mmat +def calc_surf_diff(c_surf, muR_surf, cwet, D, T, E_D_surf): + N = np.size(c_surf) + dxs = 1./N + muR_surf_long = np.empty(N+2, dtype=object) + muR_surf_long[1:-1] = muR_surf + muR_surf_long[0] = muR_surf[0] + muR_surf_long[-1] = muR_surf[-1] + c_surf_long = np.empty(N+2, dtype=object) + c_surf_long[1:-1] = c_surf + c_surf_long[0] = cwet + c_surf_long[-1] = cwet + c_surf_short = utils.mean_linear(c_surf_long) + D_eff = D/T*np.exp(-E_D_surf/T + E_D_surf/1) + surf_diff = D_eff*(np.diff(c_surf_short*(1-c_surf_short)*np.diff(muR_surf_long)))/(dxs**2) + return surf_diff + + def calc_flux_diffn(c, D, Dfunc, E_D, Flux_bc, dr, T, noise): N = len(c) Flux_vec = np.empty(N+1, dtype=object) @@ -861,7 +887,7 @@ def calc_muR(c, cbar, config, trode, ind): # u_y_tmp[2:,2:] = u_y # u_y_tmp[0,2:] = u_y[0,:] # u_y_tmp[1,2:] = u_y[0,:] -# # u_y_tmp[N_hal_x+2:,1:] = u_y[N_hal_x+1,:] +# # u_y_tmp[N_hal_x+2:,1:] = u_y[N_hal_x+1,:] # e1 = np.zeros((Nx+1,Ny+1), dtype=object) # e2 = np.zeros((Nx+1,Ny+1), dtype=object) From d44fbfff03a15c29fea1b4a3d91c49c83d258f81 Mon Sep 17 00:00:00 2001 From: Ombrini Date: Wed, 22 Mar 2023 15:22:15 +0100 Subject: [PATCH 27/27] tentative of surface diffusion --- configs/params_LFP2D.cfg | 16 ++++----- configs/params_system.cfg | 19 +++++----- mpet/config/configuration.py | 4 ++- mpet/electrode/materials/LiFePO42D.py | 50 ++++++++------------------- mpet/mod_electrodes.py | 10 +++--- 5 files changed, 41 insertions(+), 58 deletions(-) diff --git a/configs/params_LFP2D.cfg b/configs/params_LFP2D.cfg index d309fd20..4e56cca3 100644 --- a/configs/params_LFP2D.cfg +++ b/configs/params_LFP2D.cfg @@ -5,9 +5,9 @@ type = ACR2D discretization = 1e-9 shape = C3 -discretization_ver = 10e-9 +discretization_ver = 0.75e-9 # remember the it is half the thickness, check the normalizations -thickness = 50e-9 +thickness = 20e-9 [Material] muRfunc = LiFePO42D @@ -15,12 +15,10 @@ noise = false noise_prefac = 1e-6 numnoise = 200 Omega_a = 1.8560e-20 -# Omega_a = 0.5560e-20 kappa_x = 5.0148e-10 -kappa_y = 15e-10 -Bx = 0.1916e9 +kappa_y = 10e-10 +Bx = 0.19e9 By = 0.5e9 -# B = 1e9 rho_s = 1.3793e28 D = 1e-14 E_D = 0 @@ -31,9 +29,9 @@ D_surf = 1e-16 E_D_surf = 0 [Reactions] -surface_diffusion = false -rxnType = CIET -k0 = 30 +surface_diffusion = true +rxnType = BV +k0 = 0.5 E_A = 13000 alpha = 0.5 # Fraggedakis et al. 2020, lambda = 8.3kBT diff --git a/configs/params_system.cfg b/configs/params_system.cfg index 89fe8fb9..3cbefc4c 100644 --- a/configs/params_system.cfg +++ b/configs/params_system.cfg @@ -8,12 +8,12 @@ profileType = CC # Battery (dis)charge c-rate (only used for CC), number of capacities / hr # (positive for discharge, negative for charge) -Crate = 5 +Crate = 0.01 # Optional nominal 1C current density for the cell, A/m^2 # 1C_current_density = 12.705 # Voltage cutoffs, V Vmax = 3.6 -Vmin = 3.0 +Vmin = 3.1 # Battery applied voltage (only used for CV), V Vset = 0.12 # Battery constant power (only used for CP), W/m^2 @@ -25,9 +25,10 @@ power = 1 # This can have as many rows (segments) as desired. # Note: It's okay to leave commented lines within the segments list segments = [ - (10, 3), + (0.2,150), + # (5, 6), (0, 30), - (5, 10), + (3, 15), ] # Continuation directory. If false, begin a fresh simulation with the # specified input parameters here. Otherwise, this should be the @@ -46,9 +47,9 @@ tend = 1.2e3 # spacing between initial and final times with tsteps total outputs. tsteps = 200 # Relative Tolerance -relTol = 1e-6 +relTol = 1e-10 # Absolute Tolerance -absTol = 1e-6 +absTol = 1e-10 # Temperature, K T = 298 # Random seed. Set to true to give a random seed in the simulation @@ -95,7 +96,7 @@ Rfilm_foil = 0e-0 # sphere or cylinder -- radius # If using stddev = 0, set Npart = 1 for that electrode to avoid # wasting computational effort. -mean_c = 50e-9 +mean_c = 60e-9 stddev_c = 0e-9 mean_a = 100e-9 stddev_a = 1e-9 @@ -106,7 +107,7 @@ specified_psd_c = False specified_psd_a = False # Initial electrode filling fractions # (for disch, anode starts full, cathode starts empty) -cs0_c = 0.01 +cs0_c = 0.02 cs0_a = 0.99 [Conductivity] @@ -129,7 +130,7 @@ G_stddev_a = 0 [Geometry] # Thicknesses, m -L_c = 50e-6 +L_c = 30e-6 L_a = 50e-6 L_s = 25e-6 # Volume loading percents of active material (volume fraction of solid diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py index ec830d64..7c5fd8a0 100644 --- a/mpet/config/configuration.py +++ b/mpet/config/configuration.py @@ -759,7 +759,9 @@ def _indvPart(self): if self[trode, 'kappa_x'] is not None: self[trode,'indvPart']['kappa_x'][i, j] = self[trode,'kappa_x'] / kappa_ref if self[trode, 'kappa_y'] is not None: - self[trode,'indvPart']['kappa_y'][i, j] = self[trode,'kappa_y'] / kappa_ref + pthick = self['psd_len_ver'][trode][i,j] + self[trode,'indvPart']['kappa_y'][i, j] = self[trode,'kappa_y'] \ + / (kappa_ref * (pthick**2/plen**2)) if self[trode, 'dgammadc'] is not None: nd_dgammadc = self[trode, 'dgammadc'] * cs_ref_part / gamma_S_ref self[trode, 'indvPart']['beta_s'][i, j] = nd_dgammadc \ diff --git a/mpet/electrode/materials/LiFePO42D.py b/mpet/electrode/materials/LiFePO42D.py index bd15406d..26bee711 100644 --- a/mpet/electrode/materials/LiFePO42D.py +++ b/mpet/electrode/materials/LiFePO42D.py @@ -6,6 +6,7 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): Omega = self.get_trode_param("Omega_a") kappax = self.get_trode_param("kappa_x") kappay = self.get_trode_param("kappa_y") + ywet = self.get_trode_param("cwet") # beta_s = self.get_trode_param("beta_s") # for allowing muR_ref to be calculated if np.size(c_mat) == 1: @@ -21,43 +22,22 @@ def LiFePO42D(self, c_mat, ybar, muR_ref): Nx = np.size(c_mat, 0) dys = 1./Ny dxs = 1./Nx - ywet = 0.98*np.ones(Ny, dtype=object) - c_mat_tmp = np.zeros((Nx+2,Ny), dtype=object) + ywet = ywet*np.ones((1,Ny), dtype=object) muR_mat = np.zeros((Nx,Ny), dtype=object) actR_mat = np.zeros((Nx,Ny), dtype=object) - # from second to second to last row we have the original conc - c_mat_tmp[1:-1,:] = c_mat - # first and last row is 0.98 - c_mat_tmp[-1,:] = ywet - c_mat_tmp[0,:] = ywet - - # curvy = np.empty(Ny, dtype=object) - - curvy = np.empty_like(c_mat_tmp) - - curvy[:,0] = 2 * c_mat_tmp[:,1] - 2 * c_mat_tmp[:,0] - curvy[:,-1] = 2 * c_mat_tmp[:,-2] - 2 * c_mat_tmp[:,-1] - curvy[:,1:-1] = np.diff(c_mat_tmp, n=2, axis=1) - - curvx = np.diff(c_mat_tmp, n=2, axis=0) - # muRhomog = self.T * np.log(c_mat_tmp / (1 - c_mat_tmp)) + Omega * (1 - 2 * c_mat) - y_vert_avg = np.average(c_mat_tmp, axis=1) - y_oriz_avg = np.average(c_mat_tmp, axis=0) - for i in np.arange(1,Nx+1): - # curvy[0] = 2*c_mat_tmp[i,1] - 2*c_mat_tmp[i,0] - # curvy[-1] = 2*c_mat_tmp[i,-2] - 2*c_mat_tmp[i,-1] - # # + 2*dys*beta_s*c_mat_tmp[i-1,-1]*(1-c_mat_tmp[i-1,-1]) - # curvy[1:Ny-1] = np.diff(c_mat_tmp[i,:],2) - for j in np.arange(Ny): - y = c_mat_tmp[i,j] - muRhomog = self.T*np.log(y/(1-y)) + Omega*(1-2*y) - muR_mat[i-1,j] = muRhomog - # curvx = (c_mat_tmp[i-1,j] - c_mat_tmp[i,j]) - (c_mat_tmp[i,j] - c_mat_tmp[i+1,j]) - curv_en = -kappax*curvx[i-1,j]/(dxs**2) - kappay*curvy[i,j]/(dys**2) - muR_mat[i-1,j] += curv_en - muR_mat[i-1,j] += self.get_trode_param("Bx")*(y - y_oriz_avg[j]) - muR_mat[i-1,j] += self.get_trode_param("By")*(y - y_vert_avg[i]) - # print(muR_mat) + curvy = np.empty_like(c_mat) + curvy[:,0] = (2 * c_mat[:,1] - 2 * c_mat[:,0])/(dys**2) + curvy[:,1:Ny-1] = (np.diff(c_mat, n=2, axis=1))/(dys**2) + curvy[:,Ny-1] = (2 * c_mat[:,-2] - 2 * c_mat[:,-1])/(dys**2) + curvx = np.diff(np.concatenate((ywet,c_mat,ywet), axis=0), n=2, axis=0)/(dxs**2) + y_vert_avg = np.average(c_mat, axis=1) + y_oriz_avg = np.average(c_mat, axis=0) + muR_mat = self.T*np.log(c_mat/(1-c_mat)) + Omega*(1-2*c_mat) + muR_mat += -kappax*curvx - kappay*curvy + muR_mat += self.get_trode_param("Bx")*np.subtract(c_mat,y_oriz_avg) + for i in np.arange(Nx): + y = c_mat[i,:] + muR_mat[i,:] += self.get_trode_param("By")*(y - y_vert_avg[i]) actR_mat = np.exp(muR_mat/self.T) muR_mat += muRtheta + muR_ref return muR_mat, actR_mat diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py index e61574ad..3bac5844 100644 --- a/mpet/mod_electrodes.py +++ b/mpet/mod_electrodes.py @@ -704,9 +704,9 @@ def calc_eta(muR, muO): def get_Mmat(shape, N): r_vec, volfrac_vec = geo.get_unit_solid_discr(shape, N) - if shape == "C3": + if shape == "ACR": Mmat = sprs.eye(N, N, format="csr") - elif shape in ["sphere", "cylinder"]: + elif shape in ["sphere", "cylinder","C3"]: Rs = 1. # For discretization background, see Zeng & Bazant 2013 # Mass matrix is common for each shape, diffn or CHR @@ -714,6 +714,8 @@ def get_Mmat(shape, N): Vp = 4./3. * np.pi * Rs**3 elif shape == "cylinder": Vp = np.pi * Rs**2 # per unit height + elif shape == "C3": + Vp = Rs vol_vec = Vp * volfrac_vec M1 = sprs.diags([1./8, 3./4, 1./8], [-1, 0, 1], shape=(N, N), format="csr") @@ -728,8 +730,8 @@ def calc_surf_diff(c_surf, muR_surf, cwet, D, T, E_D_surf): dxs = 1./N muR_surf_long = np.empty(N+2, dtype=object) muR_surf_long[1:-1] = muR_surf - muR_surf_long[0] = muR_surf[0] - muR_surf_long[-1] = muR_surf[-1] + muR_surf_long[0] = muR_surf[0] - (np.diff(muR_surf)[0]/dxs)*dxs + (np.diff(muR_surf,2)[0]/dxs**2)*dxs**2 + muR_surf_long[-1] = muR_surf[-1] + (np.diff(muR_surf)[-1]/dxs)*dxs + (np.diff(muR_surf,2)[-1]/dxs**2)*dxs**2 c_surf_long = np.empty(N+2, dtype=object) c_surf_long[1:-1] = c_surf c_surf_long[0] = cwet