added Tesla NCA cell model from Marc's paper

configs/params_LiC6_Tesla.cfg

@@ -0,0 +1,25 @@
+# Anode parameters for isothermal benchmark from M. Torchio et al., J. Electrochem. Soc. 163, A1192 (2016).
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = diffn
+discretization = 2.e-7
+shape = sphere
+thickness = 20e-9
+
+[Material]
+muRfunc = Tesla_graphite
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+rho_s = 1.7438e28
+D = 1.0175453067862272e-13
+Dfunc = constant
+E_D = 5000
+
+[Reactions]
+rxnType = BV_mod01
+k0 = 0.6043
+E_A = 5000
+alpha = 0.5
+Rfilm = 0e-0

configs/params_NCA_Tesla.cfg

@@ -0,0 +1,27 @@
+# Cathode parameters for isothermal benchmark from M. Torchio et al., J. Electrochem. Soc. 163, A1192 (2016).
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = diffn
+discretization = 2.e-7
+shape = sphere
+thickness = 20e-9
+
+[Material]
+muRfunc = Tesla_NCA_Si
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+rho_s = 3.276e28
+D = 8.71599819329158e-14
+Dfunc = constant
+E_D = 5000
+
+[Reactions]
+# DZ change reaction
+rxnType = BV_mod01
+k0 = 73.708
+# Marc's parameters: 4.4381848254635963e-10[m^2.5/(mol^0.5*s)] * \rho_s
+E_A = 5000
+alpha = 0.5
+Rfilm = 0e-0

configs/params_system_Tesla.cfg

@@ -0,0 +1,112 @@
+# Cell parameters for isothermal benchmark from M. Torchio et al., J. Electrochem. Soc. 163, A1192 (2016).
+# See params_system.cfg for parameter explanations.
+
+[Sim Params]
+# Constant voltage or current or segments of one of them
+# Options: CV, CC, CCsegments, CVsegments
+profileType = CC
+# Battery (dis)charge c-rate (only used for CC), number of capacities / hr
+# (positive for discharge, negative for charge)
+Crate = -1
+# Optional nominal 1C current density for the cell, A/m^2
+# 1C_current_density = 30
+# Voltage cutoffs, V
+Vmax = 4.2
+Vmin = 2.7
+# Final time (only used for CV), [s]
+tend = 1.2e3
+# Number disc. in time
+tsteps = 200
+# Numerical Tolerances
+relTol = 1e-6
+absTol = 1e-6
+# Temperature, K
+T = 298
+# Random seed. Set to true to give a random seed in the simulation
+randomSeed = false
+# Value of the random seed, must be an integer
+seed = 0
+# Series resistance, [Ohm m^2]
+Rser = 0.
+# Cathode, anode, and separator numer disc. in x direction (volumes in electrodes)
+Nvol_c = 10
+Nvol_s = 10
+Nvol_a = 10
+# Number of particles per volume for cathode and anode
+Npart_c = 1
+Npart_a = 1
+
+[Electrodes]
+cathode = params_NCA_Tesla.cfg
+#cathode = params_LiCoO2_LIONSIMBA.cfg
+#anode = params_LiC6_LIONSIMBA.cfg
+anode = params_LiC6_Tesla.cfg
+# Rate constant of the Li foil electrode, A/m^2
+# Used only if Nvol_a = 0
+k0_foil = 1e0
+# Film resistance on the Li foil, Ohm m^2
+Rfilm_foil = 0e-0
+
+[Particles]
+# electrode particle size distribution info, m
+mean_c = 11e-6
+stddev_c = 0
+mean_a = 17e-6
+stddev_a = 0
+# Initial electrode filling fractions
+cs0_c = 0.8595
+cs0_a = 0.0142
+
+[Conductivity]
+# Simulate bulk cathode conductivity (Ohm's Law)?
+simBulkCond_c = true
+simBulkCond_a = true
+# Dimensional conductivity (used if simBulkCond = true), S/m
+sigma_s_c = 100
+sigma_s_a = 100
+# Simulate particle connectivity losses (Ohm's Law)?
+simPartCond_c = false
+simPartCond_a = false
+# Conductance between particles, S = 1/Ohm
+G_mean_c = 1e-14
+G_stddev_c = 0
+G_mean_a = 1e-14
+G_stddev_a = 0
+
+[Geometry]
+# Thicknesses, m
+L_c = 64e-6
+L_a = 83e-6
+L_s = 10e-6
+# Volume loading percents of active material (volume fraction of solid
+# that is active material)
+P_L_c = 0.7452
+P_L_a = 0.8277
+# Porosities (liquid volume fraction in each region)
+poros_c = 0.2298
+poros_a = 0.1473
+poros_s = 0.3589
+# Bruggeman exponent (tortuosity = porosity^bruggExp)
+BruggExp_c = -0.5
+BruggExp_a = -0.5
+BruggExp_s = -0.5
+
+[Electrolyte]
+# Initial electrolyte conc., mol/m^3
+c0 = 1200
+# Cation/anion charge number (e.g. 2, -1 for CaCl_2)
+zp = 1
+zm = -1
+# Cation/anion dissociation number (e.g. 1, 2 for CaCl_2)
+nup = 1
+num = 1
+# Electrolyte model,
+# Options: dilute, SM
+elyteModelType = SM
+# Stefan-Maxwell property set, see props_elyte.py file
+SMset = Tesla_electrolyte
+# Reference electrode (defining the electrolyte potential) information:
+# number of electrons transfered in the reaction, 1 for Li/Li+
+n = 1
+# Stoichiometric coefficient of cation, -1 for Li/Li+
+sp = -1

mpet/electrode/materials/Tesla_NCA_Si.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+
+def Tesla_NCA_Si(self, y, ybar, muR_ref):
+    """ Berliner et al., 2022.
+    OCV for graphite [V]
+    theta_p is the solid particle surface concentration normalized by its maximum concentration [-]
+    """
+    a = np.array([0.145584993881910,2.526321858618340,172.0810484337340,1.007518156438100,
+                 1.349501707184530,0.420519124096827,2.635800979146210,
+                 3.284611867463240]).reshape([1,-1])
+    b = np.array([0.7961299985689542, 0.2953029849791878, -1.3438627370872127, 0.6463272973815986,
+                 0.7378056244779166, 0.948857021183584, 0.5372357238527894,
+                 0.8922020984716097]).reshape([1,-1])
+    c = np.array([0.060350976183950786, 0.20193410562543265, 0.7371221766768185,
+                 0.10337785458522612, 0.09513470475980132, 0.0422930728072207, 0.1757549310633964,
+                 0.1413934223088055]).reshape([1,-1])
+    y = y.reshape([-1,1])
+
+    OCV = np.sum(a*np.exp(-((y - b)/c)**2), axis=1)
+    muR = self.get_muR_from_OCV(OCV, muR_ref)
+    actR = None
+    return muR, actR

mpet/electrode/materials/Tesla_graphite.py

@@ -0,0 +1,22 @@
+import numpy as np
+
+
+def Tesla_graphite(self, y, ybar, muR_ref):
+    """ Berliner et al., 2022.
+    OCV for graphite [V]
+    theta_p is the solid particle surface concentration normalized by its maximum concentration [-]
+    """
+    a0 = -48.9921992984694
+    a = np.array([29.9816001180044, 161.854109570929, -0.283281555638378,
+                 - 47.7685802868867, -65.0631963216785]).reshape([1,-1])
+    b = np.array([0.005700461982903098, -0.1056830819588037, 0.044467320399373095,
+                 - 18.947769999614668, 0.0022683366694012178]).reshape([1,-1])
+    c = np.array([-0.050928145838337484, 0.09687316296868148, 0.04235223640014242,
+                 7.040771011524739, 0.0011604439514018858]).reshape([1,-1])
+    y = y.reshape([-1,1])
+
+    OCV = a0 + np.squeeze(a[0,0]*np.exp((y - b[0,0])/c[0,0])) + \
+        np.sum(a[0,1:]*np.tanh((y - b[0,1:])/c[0,1:]), axis=1)
+    muR = self.get_muR_from_OCV(OCV, muR_ref)
+    actR = None
+    return muR, actR

mpet/electrolyte/Tesla_electrolyte.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+from mpet.config import constants
+
+
+def Tesla_electrolyte():
+    """ Set of parameters from Marc Berliner
+    """
+
+    def tp0(c, T):
+        return 0.455
+
+    def sigma(c, T):
+        """ K_e evaluates the conductivity coefficients for the electrolyte phase [S/m]
+        """
+        # we use the dimensional form in the equations
+        c = c*1000
+        T = T*constants.T_ref
+        return (-0.5182386306736273 +
+                - 0.0065182740160006376 * c +
+                + 0.0016958426698238335 * T +
+                + 1.4464586693911396e-6 * c**2 +
+                + 3.0336049598190174e-5 * c*T +
+                + 3.046769609846814e-10 * c**3 +
+                - 1.0493995729897995e-8 * c**2*T)
+
+    def D(c, T):
+        """
+        D_e evaluates the diffusion coefficients for the electrolyte phase [m^2/s]
+        """
+        c = c*1000
+        T = T*constants.T_ref
+        a = np.array([1.8636977199162228e-8, -1.3917476882039536e-10, 3.1325506789441764e-14,
+                      -7.300511963906146e-17, 5.119530992181613e-20, -1.1514201913496038e-23,
+                      2.632651793626908e-13, -1.1262923552112963e-16, 2.614674626287283e-19,
+                      -1.8321158900930459e-22, 4.110643579807774e-26])
+        return a[0] + a[1]*T + a[2]*c*T + a[3]*c**2*T + a[4]*c**3*T + a[5]*c**4*T + \
+            a[6]*T**2 + a[7]*c*T**2 + a[8]*c**2*T**2 + a[9]*c**3*T**2 + a[10]*c**4*T**2
+
+    def therm_fac(c, T):
+        return 1.
+
+    Dref = D(constants.c_ref/1000, 1)
+
+    def D_ndim(c, T):
+        return D(c, T) / Dref
+
+    def sigma_ndim(c, T):
+        return sigma(c, T) * (
+            constants.k*constants.T_ref/(constants.e**2*Dref*constants.N_A*(constants.c_ref)))
+    return D_ndim, sigma_ndim, therm_fac, tp0, Dref