LFMP implemented and tested

configs/params_LFMP.cfg

@@ -7,11 +7,9 @@ type = ACR2
 discretization = 1e-9
 shape = C3
 thickness = 25e-9
-surface_diffusion = false
 
 [Material]
 muRfunc = LiFeMnPO4
-logPad = false
 noise = false
 noise_prefac = 1e-6
 numnoise = 200
@@ -27,7 +25,7 @@ B2 = 0.2891e9
 rho_s = 1.341e28 
 D = 1e-18
 Dfunc = lattice
-dgammadc = 3.5e-29
+dgammadc = 0e-30
 cwet = 0.98
 
 [Reactions]

mpet/config/configuration.py

@@ -521,6 +521,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, 'B1'] is not None:
+                self[trode, 'B1'] = self[trode, 'B1'] / (kT * constants.N_A * self[trode, 'csmax'])
+            if self[trode, 'B2'] is not None:
+                self[trode, 'B2'] = self[trode, 'B2'] / (kT * constants.N_A * self[trode, 'csmax'])
             for param in ['Omega_a', 'Omega_b', 'Omega_c', 'EvdW']:
                 value = self[trode, param]
                 if value is not None:
@@ -759,6 +763,16 @@ def _indvPart(self):
                     # non-dimensional quantities
                     if self[trode, 'kappa'] is not None:
                         self[trode, 'indvPart']['kappa'][i, j] = self[trode, 'kappa'] / kappa_ref
+                    if self[trode, 'kappa1'] is not None:
+                        kappa_ref1 = (constants.k * constants.T_ref * constants.N_A
+                                      * self[trode, 'csmax'] * plen**2)
+                        self[trode, 'indvPart']['kappa1'][i, j] = (self[trode, 'kappa1']
+                                                                   / kappa_ref1)
+                    if self[trode, 'kappa2'] is not None:
+                        kappa_ref2 = (constants.k * constants.T_ref * constants.N_A
+                                      * self[trode, 'csmax'] * plen**2)
+                        self[trode, 'indvPart']['kappa2'][i, j] = (self[trode, 'kappa2']
+                                                                   / kappa_ref2)
                     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 \

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', 'k_h', 'cp', 'rhom']
 #: parameters that are defined for each particle, and their type
-PARAMS_PARTICLE = {'N': int, 'kappa': float, 'beta_s': float, 'D': float, 'k0': float,
+PARAMS_PARTICLE = {'N': int, 'kappa': float, 'kappa1': float, 'kappa2': float, 'beta_s': float,
+                   'D': float, 'k0': float,
                    'Rfilm': float, 'delta_L': float, 'Omega_a': float, 'E_D': float,
                    'E_A': float}

mpet/config/schemas.py

@@ -156,7 +156,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('kappa1', default=None): Use(float),
+                          Optional('kappa2', default=None): Use(float),
                           Optional('B', default=None): Use(float),
+                          Optional('B1', default=None): Use(float),
+                          Optional('B2', default=None): Use(float),
                           Optional('EvdW', default=None): Use(float),
                           'rho_s': Use(float),
                           'D': Use(float),

mpet/electrode/materials/LiFeMnPO4.py

@@ -1,31 +1,27 @@
 import numpy as np
 
 
-def LiFeMnPO4(self, y, ybar, muR_ref, ISfuncs=None):
-    """ LFMP chemical potential """
+def LiFeMnPO4(self, y, ybar, T, muR_ref):
+    """ Ombrini 2023 """
     muRtheta1 = -self.eokT*4.09
     muRtheta2 = -self.eokT*3.422
     stoich_1 = self.get_trode_param("stoich_1")
     stoich_2 = 1-stoich_1
-    if ISfuncs is None:
-        ISfuncs1, ISfuncs2 = None, None
-    else:
-        ISfuncs1, ISfuncs2 = ISfuncs
     y1, y2 = y
     Omga = self.get_trode_param('Omega_a')
     Omgb = self.get_trode_param('Omega_b')
     Omgc = (Omga+Omgb)*0.3
     # 2*ln(c/1-c) is used to account for electrons config entropy
-    muR1homog = 2*self.reg_sln(y1, (Omga/2)*stoich_1, ISfuncs1)
-    muR2homog = 2*self.reg_sln(y2, (Omgb/2)*stoich_2, ISfuncs2)
+    muR1homog = 2*self.reg_sln(y1, T, (Omga/2)*stoich_1)
+    muR2homog = 2*self.reg_sln(y2, T, (Omgb/2)*stoich_2)
     muR1nonHomog, muR2nonHomog = self.general_non_homog(y, ybar)
     muR1 = muR1homog + muR1nonHomog
     muR2 = muR2homog + muR2nonHomog
     # interaction between the two phases
     muR1 += stoich_2*Omgc*(1-2*y2)
     muR2 += stoich_1*Omgc*(1-2*y1)
-    actR1 = np.exp(muR1/self.T)
-    actR2 = np.exp(muR2/self.T) 
+    actR1 = np.exp(muR1/T)
+    actR2 = np.exp(muR2/T)
     muR1 += muRtheta1 + muR_ref
     muR2 += muRtheta2 + muR_ref
     return (muR1, muR2), (actR1, actR2)

mpet/electrode/reactions/BV_square.py

@@ -2,7 +2,7 @@
 
 
 def BV_square(eta, c_sld, c_lyte, k0, E_A, T, act_R=None,
-       act_lyte=None, lmbda=None, alpha=None):
+              act_lyte=None, lmbda=None, alpha=None):
     if act_R is None:
         act_R = c_sld/(1-c_sld)
     gamma_ts = (1./(1-c_sld))**2

mpet/mod_electrodes.py

@@ -123,9 +123,14 @@ def DeclareEquations(self):
         for k in range(N):
             eq1.Residual -= self.c1(k) * volfrac_vec[k]
             eq2.Residual -= self.c2(k) * volfrac_vec[k]
-        stoich_1 = self.get_trode_param("stoich_1")
-        stoich_2 = 1-stoich_1
-        eq.Residual = self.cbar() - (stoich_1*self.c1bar() + stoich_2*self.c2bar())
+        eq_cbar = self.CreateEquation("cbar")
+        if self.get_trode_param("stoich_1") is not None:
+            stoich_1 = self.get_trode_param("stoich_1")
+            stoich_2 = 1-stoich_1
+        else:
+            stoich_1 = 0.5
+            stoich_2 = 0.5
+        eq_cbar.Residual = self.cbar() - (stoich_1*self.c1bar() + stoich_2*self.c2bar())
 
         # Define average rate of filling of particle
         eq = self.CreateEquation("dcbardt")
@@ -149,25 +154,25 @@ def DeclareEquations(self):
         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
+        if self.get_trode_param("type") in ["diffn2", "CHR2", "ACR2"]:
+            # Equations for 1D particles of 2 field varible
             eta1, eta2, c_surf1, c_surf2 = self.sld_dynamics_1D2var(c1, c2, mu_O, act_lyte,
                                                                     noises)
         elif self.get_trode_param("type") in ["homog2", "homog2_sdn"]:
-            # Equations for 0D particles of 1 field variables
+            # Equations for 0D particles of 2 field variables
             eta1, eta2, c_surf1, c_surf2 = self.sld_dynamics_0D2var(c1, c2, mu_O, act_lyte,
                                                                     noises)
 
         # Define average rate of heat generation
         eq = self.CreateEquation("q_rxn_bar")
         if self.config["entropy_heat_gen"]:
-            eq.Residual = self.q_rxn_bar() - 0.5 * self.dcbar1dt() * \
+            eq.Residual = self.q_rxn_bar() - stoich_1 * self.dcbar1dt() * \
                 (eta1 - self.T_lyte()*(np.log(c_surf1/(1-c_surf1))-1/self.c_lyte())) \
-                - 0.5 * self.dcbar2dt() * (eta2 - self.T_lyte()
-                                           * (np.log(c_surf2/(1-c_surf2))-1/self.c_lyte()))
+                - stoich_2 * self.dcbar2dt() * (eta2 - self.T_lyte()
+                                                * (np.log(c_surf2/(1-c_surf2))-1/self.c_lyte()))
         else:
-            eq.Residual = self.q_rxn_bar() - 0.5 * self.dcbar1dt() * eta1 \
-                - 0.5 * self.dcbar2dt() * eta2
+            eq.Residual = self.q_rxn_bar() - stoich_1 * self.dcbar1dt() * eta1 \
+                - stoich_2 * self.dcbar2dt() * eta2
 
         for eq in self.Equations:
             eq.CheckUnitsConsistency = False
@@ -207,8 +212,12 @@ def sld_dynamics_0D2var(self, c1, c2, muO, act_lyte, noises):
 
     def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
         N = self.get_trode_param("N")
-        stoich_1 = self.get_trode_param("stoich_1")
-        stoich_2 = 1-stoich_1
+        if self.get_trode_param("stoich_1") is not None:
+            stoich_1 = self.get_trode_param("stoich_1")
+            stoich_2 = 1-stoich_1
+        else:
+            stoich_1 = 0.5
+            stoich_2 = 0.5
         # 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)
@@ -223,6 +232,14 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
             c2_surf = c2[-1]
             mu1R_surf, act1R_surf = mu1R[-1], act1R[-1]
             mu2R_surf, act2R_surf = mu2R[-1], act2R[-1]
+        elif self.get_trode_param("type") in ["ACR2"]:
+            (mu1R, mu2R), (act1R, act2R) = calc_muR((c1, c2), (self.c1bar(), self.c2bar()),
+                                                    self.T_lyte(), self.config, self.trode,
+                                                    self.ind)
+            mu1R_surf, act1R_surf = mu1R, act1R
+            mu2R_surf, act2R_surf = mu2R, act2R
+            c1_surf = c1
+            c2_surf = c2
         eta1 = calc_eta(mu1R_surf, muO)
         eta2 = calc_eta(mu2R_surf, muO)
         if self.get_trode_param("type") in ["ACR2"]:
@@ -255,8 +272,10 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
 
         # Get solid particle fluxes (if any) and RHS
         if self.get_trode_param("type") in ["ACR2"]:
-            RHS1 = stoich_1*np.array([self.get_trode_param("delta_L")*self.Rxn1(i) for i in range(N)])
-            RHS2 = stoich_2*np.array([self.get_trode_param("delta_L")*self.Rxn2(i) for i in range(N)])
+            RHS1 = stoich_1*np.array([self.get_trode_param("delta_L")
+                                      * self.Rxn1(i) for i in range(N)])
+            RHS2 = stoich_2*np.array([self.get_trode_param("delta_L")
+                                      * self.Rxn2(i) for i in range(N)])
         elif self.get_trode_param("type") in ["diffn2", "CHR2"]:
             # Positive reaction (reduction, intercalation) is negative
             # flux of Li at the surface.
@@ -299,7 +318,7 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
             eq1.Residual = LHS1_vec[k] - RHS1[k]
             eq2.Residual = LHS2_vec[k] - RHS2[k]
 
-        if self.get_trode_param("type") in ["ACR"]:
+        if self.get_trode_param("type") in ["ACR","ACR2"]:
             return eta1[-1], eta2[-1], c1_surf[-1], c2_surf[-1]
         else:
             return eta1, eta2, c1_surf, c2_surf

mpet/props_am.py

@@ -196,7 +196,7 @@ def non_homog_round_wetting(self, y, ybar, B, kappa, beta_s, shape, r_vec):
             muR_nh = B*(y - ybar) - kappa*curv
         elif mod2var:
             stoich_1 = self.get_trode_param("stoich_1")
-            stoich_2 = 1 -stoich_1
+            stoich_2 = 1 - stoich_1
             ybar_avg = stoich_1*ybar[0]+stoich_2*ybar[1]
             y1 = y[0]
             y2 = y[1]
@@ -239,7 +239,16 @@ def general_non_homog(self, y, ybar):
                     muR_nh = self.non_homog_rect_fixed_csurf(
                         y, ybar, B, kappa, cwet)
                 elif mod2var:
-                    raise NotImplementedError("no 2param C3 model known")
+                    # only type using this is ACR2
+                    kappa1 = self.get_trode_param("kappa1")
+                    kappa2 = self.get_trode_param("kappa2")
+                    kappa = (kappa1,kappa2)
+                    B1 = self.get_trode_param("B1")
+                    B2 = self.get_trode_param("B2")
+                    B = (B1,B2)
+                    cwet = self.get_trode_param("cwet")
+                    muR_nh = self.non_homog_rect_fixed_csurf(
+                        y, ybar, B, kappa, cwet)
             elif shape in ["cylinder", "sphere"]:
                 beta_s = self.get_trode_param("beta_s")
                 r_vec = geo.get_unit_solid_discr(shape, N)[0]