Determining the Length Scale of Transport Impedances in Li-Ion Electrodes: Li(Ni_0.33Mn_0.33Co_0.33)O₂

Figure 2a shows the configuration of the coin cell used in this study. Based on such configuration, two models were developed following the development by Knehr²² et al. and Brady et al.²³ In order to eliminate electrode-scale transport limitations, thin ${{\rm{LiNi}}}_{0.33}{{\rm{Mn}}}_{0.33}{{\rm{Co}}}_{0.33}{{\rm{O}}}_{2}\left({{\rm{NMC}}}_{111}\right)$ cathodes $\left(25\,\mu {\rm{m}},\,{\rm{mass}}\,{\rm{loading}}\,=\,2.84\,{\rm{mg}}\,{{\rm{cm}}}^{-2}\right)$ were used in the experimental study and in the present analysis. Both scaling arguments and simulations confirmed that electrode-scale effects were negligible.

Zoom In Zoom Out Reset image size

Nevertheless, a so-called pseudo 2D modeling paradigm²⁴ was used to most easily connect the smaller-length scale simulations to the parameters that are controlled during fabrication. Equations describing transport of lithium inside the NMC agglomerate or crystal are coupled to the electrode scale equations. As shown in Fig. 1c, in the crystal model, we assume that the transport through agglomerates is fast and do not contribute to voltage loss, thus the agglomerate scale transport loss is neglected. In contrast, crystal scale transport loss is neglected in the agglomerate model. All the coupled equations are solved simultaneously using the BAND(J) algorithms,²⁵ which gives values of dependent variables as a function of position and time. Detailed governing equations and boundary conditions are shown in Table I. The coefficients for the reversible potential in Table I, Eq. 11 were obtained by fitting the expression to experimental GITT (C/10) data taken after one hour of relaxation at 11 different state of discharge. While both models describe similar physics, there are differences. For example, the description of the interfacial reaction kinetics appears as a boundary condition in a crystal-scale model but in the governing equation of an agglomerate-scale model.

Table I.  Governing equations and boundary conditions for mathematical model of NMC cell.

	Electrode Scale Equation
	Governing Equation	Li Anode	Separator Edge	Current Collector
(1) Solid State Current (${i}_{1}$)	$\left(1-{\epsilon }\right)\sigma {{\rm{\nabla }}}^{2}{{\rm{\Phi }}}_{1}-a{i}_{in}\,=\,0$	${i}_{1}\,=\,0$	${\rm{\nabla }}\cdot {i}_{1}\,=\,0$	${i}_{1}\,=\,{i}_{applied}$
(2) Electrolyte Current (${i}_{2}$)	${\rm{\nabla }}\cdot \left(\kappa {\rm{\nabla }}{{\rm{\Phi }}}_{2}\right)\,+\,{\epsilon }F\left({z}_{+}{D}_{+}+{z}_{-}{D}_{-}\right){{\rm{\nabla }}}^{2}{c}_{0}\,+\,a{i}_{in}\,=\,0$	${{\rm{\Phi }}}_{2}\,=\,0$	${\rm{\nabla }}\cdot {i}_{2}\,=\,0$	${i}_{2}\,=\,0$
	$\left(\kappa \,=\,{\epsilon }{F}^{2}\left({z}_{+}^{2}{u}_{+}+{z}_{-}^{2}{u}_{-}\right){c}_{0}\right)$
(3) Electrolyte Concentration (${c}_{0}$)	${\epsilon }\tfrac{\partial {c}_{0}}{\partial t}\,=\,-{\rm{\nabla }}\cdot {N}_{+}+\tfrac{a{i}_{in}}{F}$	${N}_{+}\,=\,\tfrac{{i}_{applied}}{F}$	${\rm{\nabla }}\cdot {N}_{+}\,=\,0$	${N}_{+}\,=\,0$
	$\left({N}_{+}\,=\,-{D}_{0,eff}{\rm{\nabla }}{c}_{0}-\tfrac{F}{RT}{z}_{+}{D}_{+}{c}_{0}{\rm{\nabla }}{{\rm{\Phi }}}_{2}\right)$
(4) Solid State Concentration (${c}_{s}$)	$\left(1-{\epsilon }\right)\tfrac{\partial {c}_{s}}{\partial t}\,=\,-\tfrac{a\cdot {i}_{in}}{F}$	—	—	—
Crystal Scale equation
	Governing equation	Crystal Center	Crystal Edge
(5) Lithium balance (${{c}_{s}}^{x}$)	$\tfrac{\partial {{c}_{s}}^{x}}{\partial t}\,=\,{D}_{x}{{\rm{\nabla }}}^{2}{{c}_{s}}^{x}$	${\rm{\nabla }}{{c}_{s}}^{x}\,=\,0$	$-{D}_{x}{\rm{\nabla }}{{c}_{s}}^{x}\,=\,\tfrac{{i}_{in}}{F}$
Agglomerate Scale equation
	Governing equation	Agglomerate Center	Agglomerate Edge
(6) Solid State Current (${i}_{1}^{agg}$)	${\rm{\nabla }}\cdot {i}_{1}^{agg}\,=\,{a}^{agg}{i}_{in}^{agg}$	${i}_{1}^{agg}\,=\,0$	${i}_{1}^{agg}\,=\,{i}_{surface,agg}$
	$\left({i}_{surface,agg}\,=\,\tfrac{\partial {C}_{s}}{\partial t}F{V}_{agg}\left(1-{{\epsilon }}_{agg}\right)/{A}_{agg}\right)$
(7) Electrolyte Current (${i}_{2}^{agg}$)	${\rm{\nabla }}\cdot {i}_{2}^{agg}\,=\,{a}^{agg}{i}_{in}^{agg}$	${i}_{2}^{agg}\,=\,0$	${{\rm{\Phi }}}_{2}^{agg}\,=\,{{\rm{\Phi }}}_{2}$
(8) Electrolyte Concentration (${c}_{0}^{agg}$)	${{\epsilon }}^{agg}\tfrac{\partial {c}_{0}^{agg}}{\partial t}\,=\,-{\rm{\nabla }}\cdot {N}_{+}^{agg}+\tfrac{{a}^{agg}{i}_{in}^{agg}}{F}$	${N}_{+}^{agg}\,=\,0$	${c}_{+}^{agg}\,=\,{c}_{bulk}$
	$\left({N}_{+}^{agg}\,=\,-{{\epsilon }}^{agg}{D}_{agg}{\rm{\nabla }}{c}_{0}^{agg}-\tfrac{F}{RT}{z}_{+}{D}_{agg}{c}_{0}^{agg}{\rm{\nabla }}{{\rm{\Phi }}}_{2}^{agg}\right)$
(9) Solid State Concentration (${c}_{s}^{agg}$)	$\left(1-{{\epsilon }}^{agg}\right)\tfrac{\partial {c}_{s}^{agg}}{\partial t}\,=\,-\tfrac{{a}^{agg}{i}_{in}^{agg}}{F}$	—	—
Electrochemical Reaction Rate
(10) Current Density (${i}_{in}$)	${i}_{in}\,=\,{i}_{0}\left[\exp \left(\tfrac{{\alpha }_{a}F\eta }{RT}\right)-\exp \left(-\tfrac{{\alpha }_{c}F\eta }{RT}\right)\right]$	$\eta \,=\,{{\rm{\Phi }}}_{1}-{{\rm{\Phi }}}_{2}-U$
Exchange Current Density (${i}_{0}$)	${i}_{0}\,=\,F{k}_{rxn}{c}_{0}^{{\alpha }_{a}}{c}_{\alpha }^{{\alpha }_{c}}{\left({c}_{\alpha ,\max }-{c}_{\alpha }\right)}^{{\alpha }_{a}}$
Reversible Potential
(11) $\begin{array}{lll}U & = & {U}_{ref}+\tfrac{RT}{F}\,\mathrm{ln}\left[\left(\tfrac{c}{{c}_{0}}\right)\left(\tfrac{1-{\bar{c}}_{\max }}{{\bar{c}}_{\max }}\right)\right]\\ & & +\displaystyle {\sum }_{k=0}^{10}{A}_{k}\left[{\left(2{\bar{c}}_{\max }-1\right)}^{k+1}-\tfrac{2{\bar{c}}_{\max }k\left(1-{\bar{c}}_{\max }\right)}{{\left(2{\bar{c}}_{\max }-1\right)}^{1-k}}\right]\end{array}$					${\bar{c}}_{\max }\,=\,\tfrac{{c}_{s}-{c}_{s,0}}{{c}_{s,\max }-{c}_{s,0}}$
Parameter	Value	Parameter	Value	Parameter	Value
${c}_{s,\max }\left({\rm{mol}}\,{{\rm{cm}}}^{-3}\right)$	0.0262	A3	0.023162	A8	−0.259945
${U}_{ref}\left(V\right)$	3.868568	A4	−0.037789	A9	−0.589845
A0	−0.201805	A5	−0.330780	A10	0.052014
A1	0.112340	A6	0.239297
A2	−0.048369	A7	0.778712

The parameter estimation method follows the development of Brady et al.²⁶ A designated parameter space was sampled using Sobol sequences from Python module sobol_seq.²⁷ After simulations are generated using different sets of parameters from sobol sampling, the residual sum of squares (RSS) is calculated to describe how well the simulations emulate the experimental data. The RSS value of each parameter set is then fed into a Markov-chain monte-carlo (MCMC)²⁸ method to show statistical distributions of parameters. The parameter values are assumed to form a normal distribution, in which the mean value is the most likely parameter to emulate real electrochemical performance, and the standard deviation indicates the uncertainty in this parameter estimation.

In our study, a three-parameter space comprised of either the agglomerate or crystal scale diffusion coefficient (${D}_{agg}$ or ${D}_{x}$), reaction rate constant ${k}_{rxn}$ of lithium intercalation inside NMC materials and contact resistance (${R}_{ct}$) is sampled. The assumption that anode impedances can be neglected was tested to confirm the modeling approach. Results indicate that the anode overpotential is on the order of 1–2 mV for 0.85 mA cm⁻² (2C) discharge. It was concluded that the estimation of the agglomerate-scale diffusion coefficient is not impacted by inclusion of the anode overpotential.

Cathodes were cast on ∼18 μm Al foil with doctor blades using a mixture of $LiN{i}_{0.33}M{n}_{0.33}C{o}_{0.33}{O}_{2}$ (MSE Supplies LLC), carbon black (Timcal), and PVDF (Arkema Kynar 761) at a mass ratio of 9:0.5:0.5, dissolved in N-Methyl-2-Pyrrolidone (NMP) (Sigma-Aldrich). The as-casted electrode was heated on a hot plate to 110 °C to evaporate NMP solvent and cathode with mass loading ∼2.84 mg cm⁻² could be obtained. Such 25 μm NMC₁₁₁ cathode contains NMC agglomerates with average diameter of 10 μm and NMC crystals with an average diameter of 400 nm.

CR2032 Coin cells were made inside an argon filled glovebox. A commercial electrolyte, 1 M $LiP{F}_{6}$ in EC/DEC (5:5, w/w) (Gotion Inc.) was used as liquid electrolyte. 250 μm lithium metal foils were used as counter electrode. After assembly, two formation cycles at a C/10 rate were conducted prior to the galvanostatic and GITT experiments, which were performed employing a Landt battery tester.

Mass transport processes occur at a minimum of at least three different length scales inside the NMC electrode: electrode scale, agglomerate scale and crystal scale. In order to characterize the agglomerate scale and crystal scale mass transport, it is helpful to design experiments in which variations on the electrodes scale are negligible. Intuitively, this implies the use of thin electrodes, and Brady et al.²⁶ showed that mock experiments can quantitively show when a parameter is not important. Specifically, in this case, when the simulated uncertainty ${\sigma }_{{D}_{0}}$ of the electrode-scale effective diffusion coefficient is large and ${\sigma }_{{D}_{agg}}$ is small, the electrode scale is not important. Figure 3 shows the normalized parameter estimation results for cathodes with different thickness (mock experiments). Comparing the uncertainty of electrode scale diffusion coefficient (${\sigma }_{{D}_{0}}$) and agglomerate scale diffusion coefficient (${\sigma }_{{D}_{agg}}$), the result suggests that for thick cathodes ($\gt 100\,\mu {\rm{m}}$), the electrode scale mass transport is more important. While for thin electrodes, the electrode scale transport impedance is negligible. Thus, as we want to eliminate the electrode effect, a thin cathode (25 μm) was used to obtain experimental data for further study.

Zoom In Zoom Out Reset image size

To compare the simulation results derived from the two alternative hypotheses, the agglomerate model and crystal model were trained using galvanostatic experiments. The cell was operated between 3.0 V and 4.3 V under 0.2 C, 0.33 C, 0.5 C, 1 C and 2 C, in which 1 C = 150 mA g⁻¹. For model training, models were fitted to discharge voltage profiles under the 5 current rates, and a single set of parameters was obtained. For the agglomerate model, the diffusion coefficient ${D}_{agg}\,=\,\left(1.2\pm 0.2\right)\,\times \,{10}^{-9}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1},$ which is within a reasonable range compared with others' reports.^29,30 This low diffusion coefficient might result from the close packing of NMC crystals inside the agglomerate, which can be observed in Fig. 1b. For the crystal scale model, a value of ${D}_{x}\,=\,\left(5.2\pm 2.6\right)\,\times \,{10}^{-12}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}$ is obtained. The diffusion coefficient obtained from a crystal scale model is characterized by a significantly higher uncertainty. The higher uncertainty is because the parameter estimation algorithm suggests that the transport impedance is less dominant in the crystal-scale model.

Figure 4 shows the experimental voltage profiles and simulation results for both models. Simulation results do not extend to the experimental cutoff (2.8 V) because of numerical limitations. Both models agree with experimental data for C-rates below 1 C. For the 2 C rate, the agglomerate model still agrees well with experimental performance, while the crystal model is less successful at the large depth of discharge.

Zoom In Zoom Out Reset image size

Another method for model selection or hypothesis testing is to compare parameter estimates of the two models when trained with different types of experiments. Ideally, the derived parameters are consistent. Here, models were trained by (1) galvanostatic discharge experiments and (2) relaxation (voltage recovery) experiments. Figure 5 show the probability distribution function for the diffusion coefficient obtained by the two models when trained with discharge data and relaxation data. The experimental variation (${s}_{\exp }$) is estimated to be 20 mV. It can be seen that for the agglomerate model, discharge and relaxation experiments yield very similar estimates. This is seen in Fig. 5c, where ${D}_{discharge}/{D}_{relaxation}\,=\,0.85\pm 0.2.$ In contrast, the mean crystal-scale state diffusion coefficient (${D}_{x}$) estimated from discharge and relaxation is $5.2\,\times \,{10}^{-12}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}$ and $1.7\,\times \,{10}^{-12}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}.$ In Fig. 5c, we see ${D}_{discharge}/{D}_{relaxation}\,=\,3.1\pm 1.7$ for the crystal model.

Zoom In Zoom Out Reset image size

For an ideal system with experimental measurements with no variance, the ratio of diffusion coefficients should be one. In our case, for the agglomerate model, the distribution of ${D}_{discharge}/{D}_{relaxation}$ is narrow and has a mean value relatively close to 1. To be more quantitative, the agglomerate model has a 97% chance that ${D}_{discharge}/{D}_{relaxation}$ is between 0.6 and 1.4. For the crystal model, there is only a 17% chance that the ratio falls in this range. Thus, by comparing parameter values from different experimental sets for model selection, the agglomerate model has a much higher probability to be correct.

Another means of model selection is to evaluate its prediction ability. Similar to leave-one-out cross validation,³¹ the model is trained on a combination of data sets to obtain parameters. Then the trained model is used to predict the experimental performance of a test set that is not included in the training sets. In this study, the training sets and test sets were selected from 0.1C, 0.2C, 0.33C, 0.5C, 1C and 2C discharge experiments. To start, the models were trained by the lowest rate (0.1C) to predict the performance at 0.2C. Then the models were trained using a combination of 0.1C and 0.2C to predict 0.33C performance. For the last step, all other current rates was used to predict the performance of 2C discharge. This is named "forward prediction." For "backward prediction," The models are trained by 2C to predict 1C, etc. Figure 6 shows the average prediction error, which is the average deviation of prediction from experimental data (${e}_{sim}$) normalized by the experimental variation (${s}_{\exp }\,=\,20\,{\rm{mV}}$):

$$\begin{eqnarray}&&{\bar{e}}_{sim}\,=\,\displaystyle \frac{{e}_{sim}}{{s}_{\exp }}\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

For both forward and backward prediction, the prediction errors are similar for C-rates less than 0.5C. However, the error starts to differentiate at the higher discharge rates. For 1C prediction, the backward method for the crystal model has a high normalized error (${\bar{e}}_{sim}\,\sim \,3$), while the error for the agglomerate model is relatively low (${\bar{e}}_{sim}\,\sim \,1.2$). For 2C prediction, the forward prediction by the agglomerate model has much lower error than crystal model. This further suggests that the agglomerate scale model is in better agreement with experiment than the crystal-scale model. The deviation at 2C is also seen in Fig. 4, but perhaps in a less convincing manner.

In some cases, it may be desirable to predict the changes in voltage rather than the voltage, especially when an experimental uncertainty (such as a contact resistance or anode impedance in a cathode study) is difficult to eliminate. In such cases, models can be trained on for example a dV/dQ curve, which is the derivative of the voltage with respect to discharge capacity. The two models were trained using dV/dQ profiles for 0.33C, 0.5C, 1C and 2C and the resulted simulation results are shown in Fig. 7. Just as in Fig. 3, by fitting to dV/dQ curve, the two models are only distinct at 2C. The agglomerate model agrees better with experimental dV/dQ. Results gave ${D}_{agg}\,=\,\left(1.2\pm 0.2\right)\,\times \,{10}^{-9}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1},$ identical to values obtained by fitting V vs time. For the crystal scale model, ${D}_{x}\,=\,\left(5.9\pm 2.5\right)\,\times \,{10}^{-12}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1},$ which are also close to values obtained by fitting V vs time.

Zoom In Zoom Out Reset image size

Often, experimental investigations are not closely integrated with modeling development. Thus, a broadly applied theory is often used to estimate diffusion coefficient^32,33:

$$\begin{eqnarray}&&\tilde{D}\,=\,\displaystyle \frac{4}{\pi \tau }{\left(\displaystyle \frac{{m}_{B}{V}_{M}}{{M}_{B}S}\right)}^{2}{\left(\displaystyle \frac{{\rm{\Delta }}{E}_{s}}{{\rm{\Delta }}{E}_{t}}\right)}^{2}\,\left(\tau \,\ll \,{L}^{2}/\tilde{D}\right)\end{eqnarray} \tag{ 2 }$$

where $\tau \left(s\right)$ is the time duration for each discharge step, ${m}_{B}$(g) is the mass of NMC, ${M}_{B}$ (g mol⁻¹) is the atomic weight of NMC, ${V}_{M}$ (cm³ mol⁻¹) is the molar volume of the NMC materials, S (cm²) is the area of electrode. ${\rm{\Delta }}{E}_{s}\left(V\right)$ and ${\rm{\Delta }}{E}_{t}\left(V\right)$ are the change of voltage between relaxations and during discharge.

Note that Eq. 2 is derived assuming a crystal scale model, and the inequality constraint arises because the physics used to derive Eq. 2 assumes semi-infinite diffusion.³⁴ The use of Eq. 2 may lead to, for example, a conclusion that the diffusion coefficient is a function of the state of charge (SOC). While that is theoretically possible if not probable, whether it leads to improved agreement with experiment is a valid question.

To study the apparent SOC dependency of the diffusion coefficient, voltage recovery data from GITT experiments are used to determine the diffusion coefficient values from the two models, with results shown in Table II. Instead of absolute diffusion coefficients, values of D/L² are shown, where L is the diffusion length. Also, D/L² values calculated from experimental GITT data by Eq. 2 are listed. Figure 8 shows D/L² from agglomerate model fitting, crystal model fitting and from Eq. 2 calculation at 1C.

Table II.  ${\rm{D}}/{{\rm{L}}}^{2}$ values from fitting agglomerate model (Agg) and crystal model (Xtal) to relaxation experiments at different SOC, and ${\rm{D}}/{{\rm{L}}}^{2}$ values calculated by Eq. 2.

SOC	C/2			1C			2C
	Equation 2	Agg	Xtal	Equation 2	Agg	Xtal	Equation 2	Agg	Xtal
83.3%	0.97	6.97	3.94
80%				1.56	4.70	4.72
66.7%	0.88	5.82	2.27				1.68	6.17	3.68
60%				1.42	4.69	3.32
50%	0.77	7.03	1.64
40%				1.08	4.67	1.96
33.3%	0.4	6.16	0.87				0.55	5.33	2.06
20%				0.42	4.47	1.96
$\mu$	0.76	6.50	2.17	1.12	4.63	2.99	1.12	5.75	2.87
$\sigma$	0.22	0.52	1.13	0.44	0.09	1.14	0.56	0.42	0.81

Zoom In Zoom Out Reset image size

As shown in Table II and Fig. 8, for the agglomerate model, the D/L² value is relatively constant with state of charge. In short, a constant diffusion-coefficient model is consistent with experiment. The crystal scale model suggests that D_x increases with state of charge (this corresponds to decreasing D_x with increasing solid-state Li concentrations). Note however, that the value of D_x remains inconsistent with values estimated from the discharge curve (cf., Fig. 4). Furthermore, diffusion coefficients fit through agglomerate and crystal models are consistent with values obtained by Eq. 2 only to an order of magnitude.