Add lithium plating model

[Scottmar93]

Add the lithium plating model from:

Yang, X. G., Leng, Y., Zhang, G., Ge, S., & Wang, C. Y. (2017). Modeling of lithium plating induced aging of lithium-ion batteries: Transition from linear to nonlinear aging. Journal of Power Sources, 360, 28-40.

[felipe-salinas]

Hi @Scottmar93 and Pybamm team,

Thanks for the program, it is a great contribution!

I've been trying to include the equations in the referenced paper and faced two problems. First, the model requires a total current j_tot = j_n (already existing)+j_sei (which has a similar implementation than the one @Scottmar93 formulated in #834, but requires the definition of a new variable c_EC, the concentration of Ethylene Carbonate at the surface of graphite, and an algebraic equation that relates c_EC to j_sei and L_sei) + j_lpl (also similar to #834). This total current j_tot should be equal to the current in the solid phase i_s in the negative electrode. However, I can't find among the sub-models the interaction between the electrode (Ohms law) and interface (Butler-Volmer and associated kinetics), in which such modification should take place. I think that there was an attempt of this kind in the submodel "full_ohm.py" (pybamm/models/submodels/electrode/ohm/full_ohm.py) in the inclusion of the term sum_j (line 67), obtained from the dictionary "reactions" (apparently unused #933). What do you think that It would be the best way to include this modification?

The second issue is that both Butler-Volmer for j_n and Tafel equations for j_sei and j_lpl should include a term in the overpotential of the type L_sei * R_sei * j_tot, which implies an iterative solution. Has this already been implemented in any of the submodels (ex. actual Butler-Bolmer formulation)?

[rtimms]

Hi @felipe-salinas thanks for getting involved! There is currently a PR open (#948) which updates how the total current density is set (to allow for extra reactions etc.), so it might be worth checking there to see if that helps (@tinosulzer not sure how close this is to being finished?). I think after that you should be able to call up j_tot (or any of the other currents) more easily as an existing variable where you need it in the lithium plating submodel.

At the moment we haven't implemented having that resistance in the overpotential, but we were discussing it the other day. It requires updating the submodel to make j part of the state vector with a corresponding algebraic equation to solve for.

[valentinsulzer]

Hi @felipe-salinas , the PR I am working on is almost finished (there is just one example failing in the tests) so it should be merged in the next few days, hopefully tomorrow. What this PR will do is automatically add the current from any new reactions which implement _get_standard_whole_cell_interfacial_current_variables to some "sum" variables, then the other submodels automatically use those. For example, the "full_diffusion" model uses

sum_s_j = variables["Sum of electrolyte reaction source terms"]

while the "full_ohm" model uses

sum_j = variables[
            "Sum of " + self.domain.lower() + " electrode interfacial current densities"
        ]

(at the moment the exception is the particle models which still only use j_n but this should be changed, I will do it before merging the PR)

So if you wait a few days, or branch off the issue-933-remove-reactions-dict branch, then you will have that functionality

[valentinsulzer]

For the resistance term in Butler-Volmer, we were planning on doing this pretty soon, but happy to leave it to you if you want (it seems like you picked up the submodel structure pretty quickly)?

Two more comments:

1. Any new equations implemented will need to be nondimensionalised, see our paper "An asymptotic derivation of a single particle model with electrolyte" for details on this
2. You might find it easier to start off by making a copy of the BasicDFN model, where all the variables are defined in one place, and then adding lithium plating to that

[felipe-salinas]

Hi @tinosulzer and @rtimms, thanks for the feedback! As suggested, I will focus then on programming the current j as a state variable to include the effect of the film resistance. I hope to give some news in a couple of days.

[valentinsulzer]

New reactions formulation is now merged into develop

[felipe-salinas]

Thanks for the modifications! I cloned the branch from #834 and took as a reference the reaction_limited.py script. So far, I was able to change j_sei to a variable, created a new submodel that includes j_sei and the variable for EC surface concentration c_ec_s, programmed algebraic equations for both, where j_sei depends from the exponential term j_sum * L_sei * R_sei. I also modified the file "full_reaction_driven_porosity.py", and included it as an option for the porosity model in the dfn.py, so that there is a reduction in the negative electrode porosity proportional to j_sei, according to the paper of Yang et. al. All modifications were made following the existing scales to produce non-dimensional equations. The Lithium plating part should be another submodel and I'll explore that from now on. I wanted to share it using a pull request, but got stuck on a local problem when installing PyBaMM as a developer to run the necessary tests:

Complete output from command python setup.py egg_info:
Traceback (most recent call last):
File "", line 1, in
File "/Users/felipesalinasbarros/dep/PyBaMM/setup.py", line 7, in
import wheel.bdist_wheel as orig
ModuleNotFoundError: No module named 'wheel'

Any idea of what could be the problem? I followed the developer install.

Regarding the intercalation current in the negative electrode j_n, I tried to made modifications to include the influence of the film resistance but I think that I don't manage well enough the structure of the program, and changes I've been doing are unlikely to succeed in the near term. I think that it would be better if you make this modification, if you already planned on doing it :)

[valentinsulzer]

Sounds great, I think moving on to the lithium plating is a good plan.

For the install, you can try pip install wheel, or start again using virtualenv instead of python3 -m venv, which will install the wheel for you

Yeah, I thought the film resistance thing would be quite challenging, will take a look :)

[felipe-salinas]

Thanks @tinosulzer, I was able to install properly the wheel. Regarding the issue, I've been able to code the SEI and plating reactions, and to include the film resistance in the Butler-Volmer equation, by adding a submodel that is initialized in the dfn.py if a sei option is called, only for the negative electrode, and implements the equation referenced in #984, and changes the base_interface.py function get_coupled_variables to an algebraic formulation to obtain the interfacial current, which is now a state variable.

I've run some simulations (some parameters and functions are missing in 2017 Yang, and I think I will contact the authors to obtain them) and although the variables are doing what I think is expected (j_sei is more or less a constant reaction were the rate is reduced with the increase of L_sei, due to the equation regarding the concentration of EC, j_lpl is a smaller current than j_sei that grows when the voltage increases during charge), I don't register any capacity loss. I was wondering if it would be necessary to implement an update of the state of charge, same as the approach taken by 2004 Ramadass, where the capacity lost to SEI (and in this particular case, to lithium plating) is subtracted from the state of charge of the cathode, and updated for the next cycle.

[valentinsulzer]

Great job @felipe-salinas !

I am also working on adding the film resistance so we should be careful to avoid duplication.
Regarding the capacity loss, I think it should be possible to see a cycle-to-cycle capacity loss without having to update externally, by making sure that the j_sei and j_lpl currents divert current away from the main reaction (so that the anode isn't charged as much as the cathode is discharged). This isn't currently implemented in the existing SEI models, so I am fixing that too as part of #984.

When you say you added SEI, is that the separate from the SEI submodels that were added last week?

I am nearly done with #984, so I think the easiest thing is if I get you to review that, and then add your lithium plating model after?

[felipe-salinas]

Thanks for the update @tinosulzer,

I just added the film resistance to understand how to do it, and to explore if that could be the problem that I have with (no)capacity loss. Maybe I made a mistake when adding the side reaction currents to the interface current. I'm following your program in #984 to develop the new submodels SEI limited by EC and Lithium plating. Regarding the former, there are a couple of differences with the existing SEI models, mainly, a new variable representing the concentration of EC in the surface of graphite with its corresponding algebraic equation, and the definition of certain constants.

[felipe-salinas]

Is it possible that there is an error in the definition of the scale used for j_sei? I think that the actual script of sei models (ex. reaction_limited.py) has two errors: one that the scale chosen does not go in line with the time derivative d/dt L_inner = -j_inner (mainly, t_discharge is not implied in this equation so it shouldn't be part of the scale for j_sei); and a second, that given that the scale is different than the one for the intercalation current, j_sei cannot be directly summed to j_n to form j_tot using the function get_standard_whole_cell_interfacial_current_variables. I guess that was the intention of the parameter Gamma_SEI_n defined in sei_parameters.py, but I don't see the use of this factor in any place.

[valentinsulzer]

I'll check those scalings over, they could be wrong. In the meantime, can you take a look at the branch issue-984-butler-volmer-resistance and run calendar_ageing.py and let me know if those results look sensible to you? It shows a loss of capacity in the negative electrode when zero current is applied.

Happy for you to add your EC-limited SEI model to the sei folder alongside the existing ones. By the way, what does "EC" stand for?

[felipe-salinas]

Hi @tinosulzer, I ran the script, and calendar aging showed a voltage reduction (which should be the case if there are side reactions) and the capacity is reduced during cycling. It seems to work as it should, well done!

I have another remark. In the calculation in the submodel base_kinetics.py where the current j included in the overpotential eta_r is defined as the negative electrode interfacial current density, I think it should be modified to the sum of negative electrode interfacial current density, to include the overpotential produced by j_sei.

The EC stands for Ethylene Carbonate :), and the algebraic equation proposed in Yang 2017 has the effect of reducing the diffusion constant of EC when L_sei increases, I think, producing a ''passivating effect" in the long term of SEI growth.

[valentinsulzer]

The way it's implemented in base_kinetics.py, phi_s and phi_e are algebraic variables (which get solved by Ohm's law and the MacInnes equation). For example:

i_s = -sigma * grad(phi_s)
div(i_s) = -(j + j_sei)

and

j = j(phi_s - phi_e - U - R_sei * L_sei * j)
j_sei = j_sei(phi_s - phi_e - U_sei - R_sei * L_sei * j)

All base_kinetics.py does is create the first j as a function of phi_s and phi_e , and later j gets plugged in to the div(i_s) equation

Regarding scales, you're right there is a missing factor of Gamma_sei. We're currently trying to decide whether it would be more appropriate to scale j_sei the same way as the intercalation current. The discharge timescale is fine though, because you want to scale all your equations with the same timescale, and the discharge timescale is the one that is used everywhere else

[felipe-salinas]

Thanks! I was missing the fact that all time derivatives have to include a time scale.

[valentinsulzer]

The branch with the SEI film resistance has been merged into develop. The nondimensionalisation is updated to fix the missing factor you pointed out (j_sei now scales like j). If you merge it into your lithium plating branch, we can go from there

[DrSOKane]

Hi Tino/Felipe,

My paper on reversible Li plating modelling has been published this week (https://dx.doi.org/10.1149/1945-7111/ab90ac) and I'm currently working on adding that model into PyBaMM. Ideally we want to have a minimum of three options: no plating, reversible plating and irreversible plating. It appears the two codes being developed could complement one another nicely.

Many thanks, Simon

[valentinsulzer]

Thanks @DrSOKane and congrats on the publication! I'll let you and Felipe coordinate adding the lithium plating model but let me know if you have any questions

[felipe-salinas]

Hi @DrSOKane,

thanks for joining and for sharing the parameters in the supplementary material of your article. You are very welcomed to add the plating model, because for the moment I'm focusing in merging the SEI model, and on validating a simulation to check that is implemented correctly. There are some variables and files that will be used in both models (like L_sei which will increase with the plating current as well, and the reduction in porosity from it) and maybe is useful to check first check what I did in #1009. I think that the most general approach is to implement only a submodel for reversible plating, as irreversible plating is a particular case where the backward reaction is set to 0.

[DrSOKane]

Thanks @felipe-salinas, this is good advice. I will push a separate LithiumPlating branch once the code is in better shape.

[maurosgroi]

Dear all,
I'm new to Pybamm and very interested to the Li-plating feature. Is there any news about it? Will it be added to the official release of the code?
Thanks a lot and best regards,
Mauro Sgroi.

[rtimms]

Hi @DrSOKane can you provide an update? thanks!

[DrSOKane]

The plating model is working but I've been having issues with my Linux subsystem that prevented me from adding it to PyBaMM. I will work on fixing this today.

[maurosgroi]

Dear DrSOKane,
thanks a lot!
Best regards,
Mauro Sgroi.

[valentinsulzer]

Closed by #1287