SEI Growth by Electron Diffusion & Migration

[MichaPhilipp]

Description

The feature contains SEI growth models based on the diffusion and migration of electrons through the SEI (L.v. Kolzenberg, A. Latz, B. Horstmann, ChemSusChem 2020, 13, 3901-3910). The model is dervied to describe SEI growth during storage and battery operation.

Motivation

This feature supports the simulation of SEI growth under different operating conditions. Thereby it extends the set of available SEI growth models on PyBaMM and enables the comparison of SEI growth for different operating protocols.

Possible Implementation

The simplest to implement is to add these lines into sei_growth.py

OCP_n = variables["Negative electrode open-circuit potential [V]"]
j0 = variables["Negative electrode exchange current density [A.m-2]"]

elif self.options["SEI"] == "electron-diffusion limited":
            # L. v. Kolzenberg, ChemSusChem 2020 (eq. 25)
            chem_pot_Li = 0.1803
            eta_int = 2/F_RT*pybamm.arcsinh(j/(2*j0))
            j_ED = -(phase_param.D_li * phase_param.c_li_0 * param.F / L_sei_outer) * pybamm.exp(-F_RT*(eta_int + OCP_n + chem_pot_Li))
            j_sei = j_ED

elif self.options["SEI"] == "electron-migration limited":
            # L. v. Kolzenberg, ChemSusChem 2020 (eq. 27a & 27b)
            chem_pot_Li = 0.1803
            eta_int = 2/F_RT*pybamm.arcsinh(j/(2*j0))
            j_migr = pybamm.minimum(param.F**2 * phase_param.D_li * phase_param.c_li_0 / (2 * param.R * T * phase_param.kappa_inner) *j*pybamm.exp(-F_RT*(eta_int + OCP_n + chem_pot_Li)),0)
            j_sei = j_migr

Additional context

In the following all cited equations refer to the equations from this work (L.v. Kolzenberg, A. Latz, B. Horstmann, ChemSusChem 2020, 13, 3901-3910).

The equations described in "Possible Implementation" origin from the general equation for the SEI current density under different operating conditions (eq. 19). Looking at different limits this equation reduces to the reaction limited (eq. 23), electron diffusion limited (eq. 25) and migration limited cases (eq. 27 a+b). Note that this "electron-migration limited" case differs from the mechanism already implemented in PyBaMM under the same name. Even though the basic equation for electromigration (see eq. 5) can yield both "migration" mechanisms, the major difference lies in the interpretation of conductivity and the origin of the gradient in the electrical potential in the SEI. Here the intercalation current of Li-ions is assumed to cause the gradient in the electrical potential. This leads then to SEI growth proportional to the intercalation current density and it goes linear with time (eq. 28).

Both mechanisms (electron diffusion & migration) depend on the overpotential of the SEI (eta_SEI), which is directly linked to the intercalation overpotential (eta_int). The latter is defined through the difference in potentials (eq. 10) and linked to the resulting intercalation current (eq. 9). In the above suggested implementation of the models we assume that j_SEI << j_int and therefore j_int ~ j. This let's us use eq. 9 directly to compute eta_int and insert it in the needed SEI overpotential (see eq. 18). Depending on the application this assumption may not be justified. Some attempts were made to keep the intercalation current a variable, and let the solver get a correct solution for j_int and j_SEI. However, I was not able to get numerical stability.

To fully capture the correct dependencies and smooth transitions from reaction to diffusion to migration limited SEI growth it would be necessary to implement eq. 19. In theory this captures the linear to square-root to linear growth regimes of the SEI, depending on transport parameters and current SEI thickness. Unfortunately i was not able to achieve any solutions due to solver problems.

Looking forward to thoughts on this.

[brosaplanella]

Hi @MichaPhilipp! Thanks for opening this. I have a few questions/comments:

• If we managed to make equation (19) work, would that mean that we only need that one and we do not need to implement (25) and (27) as two separate cases?
• Maybe the numerical issue are to do with the definition of the variables. I think the way eta_int is defined might be the cause of the problem, and instead we can access PyBaMM's already defined overpotential.
• The chem_pot_Li should not be hardcoded but instead passed as a parameter, but that's a matter of user interface so I suspect it would not be a problem numerically.

@kawaMANMI has been looking into this too and his input will also be very helpful.

Also, if you have the code already implemented, you can open a pull request so we can see the code and it will make discussions easier. Note that this will make the code available, so make sure you are not including any confidential files. Happy to help if you have never done this before.

[MichaPhilipp]

Thank you @brosaplanella for your comments!

• If it is possible to make eq. (19) work with numerical stability (also during switching between charge/discharge) this would model the transition from (23), (25) to (27). If this transition occurs and when depends then on the parameterization and also the initial conditions. However, I think it would be still benefitial to implement (25) & (27) to be able to look on the isolated limiting cases.
• This might solve the problem. It would be great to have the correct overpotential already in PyBaMM. I need to check that.
• I totally agree.

I'm excited for further input!

Since I'm busy at the moment, I'll try to open a pull request by next week. Some help would be great.

[MichaPhilipp]

In order to advance, I investigated the already implemented SEI models and the used/available overpotentials. Thereby I'd like to add some discussion points.

• the proposed electron diffusion model is an extension of an earlier developed model "interstitial-diffusion" (F. Single, A. Latz, B. Horstmann; ChemSusChem, 2018) by the influence of the intercalation overpotential (charging/discharging). In case of storage these models are basically identical. This simpler "interstitial-diffusion" model from F. Single then coincides in its form with the already implemented model "interstitial-diffusion limited" (S. Marquis, PyBaMM), depending on the definition of delta_phi. If I'm right, the used delta_phi includes a potential drop due to the resistivity of the SEI layer (~ j*L_sei*R_sei), which causes accelarated SEI growth once the SEI layer is thick enough. However, in the derivation by Single & Kolzenberg the overpotential is due to the formation of neutral Li atoms at the electrode-SEI interface, which then diffuse to the SEI-Elyte Interface. Therefore, the eta_SEI does not depend on the ohmic potential drop over the SEI. Unfortunately I'm not able to see the derivation of S. Marquis to investigate whether these descriptions are identical or not. I think it would be great to first check if these models coincide or should be differentiated (by using different overpotentials).
• to circumvent my above calculation of eta_int I found two possibilites that should yield the same overpotential needed for "electron diffusion". Firstly, use eta_SEI = delta_phi - U_sei - j*R_sei*L_sei defined in sei_growth.py with U_sei = 0. This already includes the OCP. Secondly, use the available variable "Negative electrode reaction overpotential [V]" as the intercalation overpotential eta_int. The reaction overpotential is given by eta_r = (2 * (param.R * T) / param.F / ne) * pybamm.arcsinh(j / (2 * j0 * u)) and should be identical to my previous definition. In comparison delta_phi - j*R_sei*L_sei ~ eta_r + OCP are pretty similar, however some minor differences occur. I'm not sure if these are due to additional terms in delta_phi accounting for porosities. I think it's correct and faster to use eta_r, however, I would appreciate some clarification on the difference to delta_phi.

[brosaplanella]

Apologies, this completely went off my radar.

• delta_phi is the difference between the electrode and the electrolyte potential. When substracting the relevant OCP then it becomes an overpotential. Scott Marquis' thesis (https://ora.ox.ac.uk/objects/uuid:8afdcc34-cc42-48ba-b316-96a6d0f33a45) should provide all the details needed.
• Note U_sei is a parameter, so you can set it to zero without modifying the code, just changing the parameter set you use. eta_r is the intercalation overpotential, which depending on the model of choice may or may not include the SEI resistance losses.

[kawaMANMI]

I have a few concerns about implementing the model from the paper "Solid–Electrolyte Interphase During Battery Cycling: Theory of Growth Regimes - Kolzenberg - 2020 - ChemSusChem - Wiley Online Library"

I thought it might be helpful to raise them here:

1. I doubt it's appropriate to directly implement equation (19) as during rest, it will be undefined due to L_mig becoming infinity.
2. After equation (21), they mention that for realistic parameters, L_diff << L_mig, which makes sense. However, with the parameters provided in the supporting information, this condition is not satisfied. (The lithium-ion conductivity \kappa_{Li^+, SEI} can be changed to meet this condition, but what is a reasonable value?)
3. C_Li_0, the "Reference concentration of Li atoms at the anode-SEI interface", is given as 1000 mol/m^3 (in the supporting information), while in PyBaMM and other references in the literature, the 'Lithium interstitial reference concentration [mol.m-3]' is 15 mol/m^3. On the other hand, the original reference for this value is 15 mmol/m^3 "https://chemistry-europe.onlinelibrary.wiley.com/doi/abs/10.1002/cssc.201800077". I am not sure whether this is a typo or if other people are using it incorrectly.

[MichaPhilipp]

Thank you for your points, maybe I can help there:

1. I think this should not be a problem as L_mig only occurs in the denominator and therefore the equation simplifies during rest. However, I'm not sure about the numerical stability. Maybe this can cause troubles.
2. What values do you use for the currents j and j_0? If I stick to the values in the supporting information, this condition is satisfied.
3. In my opinion the problem originates from the fact, that one cannot identify these parameters one by one. As they state in the paper only the product of the concentration, the diffusivity and the effective area is effectively used to fit the parameters. Depending on what values you use for the electrode surface area and the diffusivity will have major impact on the value of c_LI_0.

[kawaMANMI]

Thank you for your responses, which are helpful.

1. Yes, theoretically there should be no problem, but I think in the implementation we will get a zero division error. However, I just found that I added a small number (1e-28) to the denominator of L_mig, and I think it works fine.
2. I do not set any values of J and J_0 explicitly as I use the full cell model and cycling protocol in PyBaMM. I am interested in the formation cycle, so the charge/discharge is expected to be slow, and I set it as C/20. With the "Chen2020" set of parameters, J will be around 0.075, but J_0 is used as "m_ref * arrhenius * c_e0.5 * c_s_surf0.5 * (c_s_max - c_s_surf) ** 0.5" (https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/lithium_ion/Chen2020.py), which is completely different from the paper. Yes, this was one of my troubles in mixing up PyBaMM parameters and the paper. I just double-checked, and it seems I forgot the factor F/RT in the exponential to find L_diff. So, yes, the condition is satisfied if parameters in the supporting information used.

Now, J_sei in the paper is added as one more option in the sei_growth, but my focus is on how I can set parameters in J_sei in order to meet the condition (L_diff << L_mig) . I also want the paraemeters (for now) to meet third mode of the model (L_mig << L_sei ) because my main focus will be on comparing several models in the first cycle with subsequent cycles, and other modes already exist in pyabmm. I assume that L_sei <= 1e-12.

[MichaPhilipp]

Hi @kawaMANMI!

1. Great that you found a work around.
2. Good to know that the condition is satisfied. However, I experienced myself some cases in which the condition is not satisfied, e.g. for using parameter sets in PyBaMM (Chen2020). Therefore, I think the interpretation of the condition is that once it is satisfied, you should see a transition from diffusion limited to migration limited SEI growth. If it is not satisfied (L_mig << L_diff) you should be in the migration dominated regime right from beginning.

Alternatively, you can adapt the value of j_SEI_0_0 to keep the condition satisfied. The third mode is describing the long-term SEI growth, so I'm not sure whether you really want to reach it (lower L_mig for lower \kappa_{Li+,SEI} or higher currents). Depending on your purpose it might be more valuable to simulate the different modes seperately (i.e. reaction limited, diffusion limited, migration limited).