Intercalation driven porosity effects in coupled continuum models for the electrical, chemical, thermal and mechanical response of battery electrode materials

A battery cell usually consists of porous, positive and negative electrodes, a separator and a current collector (see Figure 1). The simplified porous electrode model[25] is widely used, based on the theory set forth by Newman and Tiedemann[26]. The equations that follow correspond to those in Newman and Tiedemann[26]. They are modified to account first for the deformed configuration of the electrode. Since the evolution of porosity depends on the kinematics of deformation, these effects are treated jointly in Section 2.2. We note that the porosity was assumed to be constant by Doyle et al.[25]. All equations are expressed in the deformed configuration. We only present the final equations here. Detailed derivations appear in the Appendix.

Lithium intercalation/de-intercalation induces particle swelling/contraction which manifests itself as electrode deformation at the macro-scale. Additionally, the particle and separator undergo thermal expansion.^aaaThe electrolyte is assumed not to undergo thermal expansion. Therefore the decomposition of the deformation into elastic, swelling and thermal contributions does not apply to it. The kinematics of finite strain leads to the following decomposition:

where $\mbox{\boldmath$F$}=\mbox{\boldmath$1$}+\partial\mbox{\boldmath$u$}/\partial\mbox{\boldmath$X$}$, is the total deformation gradient tensor averaged over the constituents of the material. It is multiplicatively decomposed into $\mbox{\boldmath$F$}^{\text{e}}$, $\mbox{\boldmath$F$}^{\text{c}}$ and $\mbox{\boldmath$F$}^{\theta}$, which are, respectively, its elastic, chemical (induced by lithium intercalation) and thermal components. In the absence of a body force the strong form of the mechanics problem in the current configuration is

where $T$ is the Cauchy stress tensor and $W$ is the strain energy density function.

The electrodes are composed of solid particles, electrolyte-filled pores and binders as illustrated in Figure 3.

The sum of volume fractions gives

where $\epsilon_{\text{s}}$, $\epsilon_{\text{l}}$ have been introduced before as the volume fractions of solid particles and pores (assumed filled with electrolyte), and $\epsilon_{\text{b}}$ is the binder. The subscript $(\bullet)_{0}$ denotes the corresponding volume fraction in the reference configuration. For a volume element, we denote its total volume by $\delta V$ and the volume of solid particles in it by $\delta V_{\text{s}}$. Then, we have:

where $\mbox{\boldmath$F$}_{\text{s}}$ is the deformation gradient tensor over the particle, det denotes the determinant, and we have used $\text{det}\mbox{\boldmath$F$}=\delta V/\delta V_{0}$, $\text{det}\mbox{\boldmath$F$}_{\text{s}}=\delta V_{\text{s}}/\delta V_{\text{s}_{0}}$, and $\epsilon_{\text{s}_{0}}=\delta V_{\text{s}_{0}}/\delta V_{0}$. Similarly, we have

where $\mbox{\boldmath$F$}_{\text{l}}$ and $\mbox{\boldmath$F$}_{\text{b}}$ are, respectively, the deformation gradient tensors of the pore (assumed filled with electrolyte) and the binder material. From the theory of mixtures[31] the total stress is,

We assume that the electrolyte is an ideal, incompressible fluid, implying that $\mbox{\boldmath$T$}_{\text{l}}=P_{\text{l}}\mbox{\boldmath$1$}$, uniform in space and constant in time. Furthermore, the stress in the binder is modelled as isotropic, uniform in space and constant in time so that $\mbox{\boldmath$T$}_{\text{b}}=P_{\text{b}}\mbox{\boldmath$1$}$. This implies that the elastic deformation of the binder $\mbox{\boldmath$F$}^{\text{e}}_{\text{b}}$ is also uniform in space and constant in time. This represents a strong assumption for the binder, which we will seek to relax in future communications by taking account of the microstructure. For the electrolyte, isotropy of stress is reasonable, while its uniformity and constancy merit closer examination in microstructurally based treatments.

Assuming a strain energy density function $W(\mbox{\boldmath$F$}^{\text{e}})$ with volumetric term

and the same form for the solid particle, this leads to

Computing the trace of Equation (27), using (29) and (30), we have

We write the swelling of the electrode’s volume element due to lithium intercalation and thermal expansion using empirically derived functions, $\beta$ and $\beta^{\theta}$:

Then, the volumetric deformation of the element and its solid component can be expressed as

Note that the models in Equations (32a–32d) represent a simplified approach as an alternative to rigorous homogenization. The lithium intercation-induced swelling of the solid component is due only to the active material, and not the binder. The effective swelling, applicable to the entire solid component has been represented by $\beta_{\text{s}}(C_{1})$. Thermal expansion also is assumed to occur only in the active material, but not in the binder; in this case, $\beta_{\text{s}}^{\theta}$ represents the effective thermal expansion over the solid component, since we do not consider the temperature field $\theta$ to vary over the spatial scale of active material and binder particles^bbbThis assumption is valid for the through plane temperature distribution due to the relatively thin electrode structure, but should be re-assessed for the case of large format cells which could have significant in-plane temperature distributions.. The functions, $\beta(C_{1})$ and $\beta_{\text{s}}(C_{1})$ have been fit to experimental data in this work (see Section 4.2). The more sophisticated treatment for these functions, using analytic or computational homogenization over the microstructure, will be presented in a future communication.

On substituting Equations (32a–32d) in Equation (24) we obtain

Assuming the binder to deform at constant volume during charging and discharging such that $\text{det}\mbox{\boldmath$F$}_{\text{b}}=1$, we have

For the thermal expansion functions, we have chosen

for constants $\Omega_{\text{s}}$ and $\Omega$.

Returning to Equation (18), we assume that although the particles undergo isotropic swelling in the electrode due to intercalation of lithium, there is unidirectional swelling along the normal to the separator; i.e., the $\mbox{\boldmath$e$}_{2}$-direction in Figure 4. Accordingly, we have

This assumption models the non-slip boundary condition that would be applied at the current collector-electrode interface, and which would provide a strong constraint against macroscopic expansion of the electrode in the $\mbox{\boldmath$e$}_{1}$ and $\mbox{\boldmath$e$}_{3}$ directions. Since we do not directly model the current collector, this assumption represents its mechanical effect. On a practical note for this work, $\beta(C_{1})$ and $\beta_{\text{s}}(C_{1})$ were fit to the observed expansion along $\mbox{\boldmath$e$}_{2}$ and used to obtain the numerical results of Section 4. Thermal expansion is modelled as isotropic[3]^cccThe aspect ratio of the computational domain has been chosen with the $\mbox{\boldmath$e$}_{2}$ dimension to be much smaller than the other two dimensions. Additionally, the boundary conditions do not allow displacement in the $\mbox{\boldmath$e$}_{1}$ and $\mbox{\boldmath$e$}_{3}$ directions on the surfaces of the computational domain that are perpendicular to $\mbox{\boldmath$e$}_{2}$ (see Section 4.1). The solution to the mechanics problem therefore renders the resulting total strain (driven by intercalation and thermal expansion) to also be primarily along the $\mbox{\boldmath$e$}_{2}$ direction..

Equations (18), (40), (41) and (19) allow the volume fractions to be calculated by also using Equations (35) and (37).

The separator can be treated as a medium in which the active material and binder are replaced by a porous polymer that does not undergo lithium intercalation-induced swelling, but does experience thermal expansion. Accordingly, the functions $\beta$ and $\beta_{\text{s}}$ vanish in the separator. The deformation of the solid component is $\text{det}(\mbox{\boldmath$F$}_{\text{s}})=\text{det}(\mbox{\boldmath$F$}^{\text{e}}_{\text{s}})\text{det}(\mbox{\boldmath$F$}^{\theta}_{\text{s}})$, which replaces Equation (34). The stress is $\mbox{\boldmath$T$}=\mbox{\boldmath$T$}_{\text{l}}+\mbox{\boldmath$T$}_{\text{s}}$, replacing Equation (27). The volume fraction expressions in the separator reduce to

and since the binder is absent,

Since we assume that lithium intercalation makes the spherical particles swell isotropically, the particle radius needs to be updated as follows to determine the specific area, $a_{\text{p}}=S_{\text{p}}/V_{\text{p}}=3/R_{\text{p}}$. Here, the particle volume is $V_{\text{p}}=V_{\text{p}_{0}}\det{\mbox{\boldmath$F$}_{\text{s}}}$. From equations (31), (33) and (34), this gives

from which the spherical radius can also be written as

As the particles swell/contract due to lithium intercalation/de-intercalation, the maximum concentration of lithium that they can accommodate varies. We define $C^{\text{max}_{0}}_{1}$ to be the maximum concentration of lithium in the undeformed particle. This quantity is a material constant. Furthermore, as the particle swells/contracts, the maximum number of lithium atoms that it can accommodate remains fixed. This implies that

from which, using Equation (44), $C^{\text{max}}_{1}$ is given by

We note that Equation (2) is derived by using the particle volume-averaged concentration, $C_{1}$, but ignoring diffusion of lithium in the solid particle. It is traditional to assume a spherically symmetric distribution of lithium within the particle due to uniform boundary conditions on the particle surface, with a parabolic profile along the particle radius[32]. Letting $c_{1}$ denote the intra-particle lithium concentration, its relation to the volume-averaged concentration, $C_{1}$, is

Where $R_{\text{p}}$ is the time varying particle radius due to swelling given by Equation (45). The detailed derivation appears in the Appendix.

To obtain the Galerkin weak form for each strong form (2), (4), (5-11) and (19), we multiply the strong form by the corresponding weighting function, integrate by parts and apply boundary conditions appropriately. The final weak form expressed in terms of residual equations appears below.

where $Q=Q_{\text{rxn}}+Q_{\text{rev}}+Q_{\text{ohm}}$ in the electrodes as in Equation (10) and $Q=Q_{\text{ohm}}$ in the separator as in Equation (11). The fields $w_{c_{1}}$, $w_{c_{2}}$, $w_{\phi_{\text{S}}}$, $w_{\phi_{\text{E}}}$, $w_{\theta}$ and $\mbox{\boldmath$w$}_{u}$ are weighting functions for the corresponding primary variables. Time integration is achieved by the Backward-Euler algorithm.

The analytical linearization of the residual equations (50-55) to obtain the Jacobian matrix is tedious and is fraught with the danger of algebraic mistakes. Symbolic differentiation is an option, but its speed can be a limitation for complicated nonlinear and coupled problems such as those in the present communication. An easy alternative is the use of numerical differentiation tools built into many standard solver packages. However, for a highly non-linear set of equations, numerical differentiation is inaccurate and ultimately unstable. An effective and efficient alternative is the use of algorithmic (or automatic) differentiation (AD), which works by application of the chain rule to algebraic operations and functions (polynomial, trigonometric, logarithmic, exponential or reciprocal) in the code. AD thus works to machine precision at a computational cost that is comparable to the cost of evaluation of the original equations. We use AD in this work to linearize Equations (50-55), and compute the Jacobian matrix. Specifically, we use the Sacado package, which is part of the open-source Trilinos project[36, 37].

We model a three dimensional cell (Figure 4) of size $120,000\mu\text{m}\times 128\mu\text{m}\times 85,000\mu\text{m}$. For ease of interpretation of the coupled physics in this first communication, the boundary conditions and distribution of coefficients have been chosen to render the problems to be effectively one-dimensional. The fields do not vary along the $\mbox{\boldmath$e$}_{1}$ and $\mbox{\boldmath$e$}_{3}$ directions. Therefore a single element sufficed along these directions. An element width of 1$\mu\text{m}$ was used along the $\mbox{\boldmath$e$}_{2}$ direction. Trilinear hexahedral elements were used. Because of the one-dimensional nature of the problems considered, the high element aspect ratios did not manifest in numerical ill-conditioning. In order to overcome the stiff dynamics using the Backward Euler algorithm, initial time steps of $0.1$s were used, gradually increasing to $10$s for computations at the 1 C rate, and $1$s for the 10 C rate.^dddA 1 C rate means that a battery rated at N Ah should supply N Amperes for 1 hour. Details of all parameters and initial conditions are summarized in the Appendix. The partial differential equations solved are either parabolic (Equations (51) and (54)) or elliptic (Equations (52), (53) and (55)), and so boundary conditions must be specified on each surface. Conservation of lithium ions in the electrolyte translates to zero flux boundary conditions for $C_{2}$. Since $\phi_{\text{S}}$ is only defined over the electrodes, it has boundary conditions specified on the domain boundaries, and interface conditions at the electrode-separator interfaces. Its boundary conditions are $\phi_{\text{S}}=0$ at $x_{2}=0$ (reference potential), and $-\sigma_{\text{eff}}\nabla\phi_{\text{S}}\cdot\mbox{\boldmath$n$}=I/F$ at $x_{2}=128\mu\text{m}$ (applied current). The boundary conditions on the field in the electrolyte are $\nabla\phi_{\text{E}}\cdot\mbox{\boldmath$n$}=0$ on all boundaries. Boundary conditions for the thermal field, $\theta$, correspond to conductive heat transfer to the ambient air which is assumed to be at $25^{\circ}\text{C}$. For mechanics, displacement boundary conditions are prescribed, $\mbox{\boldmath$u$}=\mbox{\boldmath$0$}$ at $x_{2}=0$, and $\mbox{\boldmath$u$}=u_{0}\mbox{\boldmath$e$}_{2}$ at $x_{2}=128\mu\text{m}$. These mechanical boundary conditions compress the cell to a fixed total strain as is common in automotive battery packs. The remaining surfaces are taken to be traction-free. The boundary and interface conditions are summarized in Table 1, where the surfaces referred to are depicted in Figure 4.

We computed a full cycle (discharging and charging) under 10 C and 1 C current rates for the cases outlined in (i) and (ii) above. For cases with decoupled porosity, the initial porosity, $\epsilon_{\text{s}_{0}}$, from which discharging happens, has been set higher than in the coupled porosity case to ensure that when fully discharged, both these models have the same porosity. This leads to smaller initial volume fractions of solid particles for cases with decoupled porosity. However, the total number of moles of lithium in the electrode is identical for both cases when they start at the same state of charge (SOC). The parameter $C^{\text{max}}_{1}$ therefore changes as expressed by Equation (47). All initial porosity distributions were taken to be uniform. For studies of the effect of porosity we considered two cases: (a) $\beta_{\text{s}}\neq 0$, $\beta=0$; i.e., particle swelling is accommodated within the electrode free volume and (b) $\beta_{\text{s}}\neq 0$, $\beta\neq 0$; particle swelling causes electrode swelling. Here, SOC was calculated by coulomb counting. Since constant current rates were applied, the SOC was defined as $\text{SOC}=1-\frac{\text{current time}}{\text{total time}}$ for discharging, or $\text{SOC}=\frac{\text{current time}}{\text{total time}}$ for charging, where the total time was for fully discharging or fully charging the cell. These definitions are equivalent to a volume averaged notion of SOC, given by $\text{SOC}=(\bar{\eta}-C_{\text{1soc}_{0}})/(C_{\text{1soc}_{100}}-C_{\text{1soc}_{0}})$, where $\bar{\eta}$ is the electrode volume average of $\eta$, recalling the definition $\eta=C_{1}/C_{1}^{\text{max}}$, and $C_{\text{1soc}_{y}}$ is the fraction of the theoretical electrode maximum lithium content at the given SOC. Thus defined, SOC is equivalent for both electrodes in the absence of side reactions or loss of lithium.

For our prototypical Li ion battery cell, the swelling function $\beta$ was fit to data of Oh et al.[38]–see Equation (56) and Figure 5–and $\beta_{\text{s}}$ was fit to data of Takami et al.[39]–see Equation (57) and Figure 6. The swelling functions were redefined as, $\beta(C_{1})=\hat{\beta}(\eta)$ and $\beta_{\text{s}}(C_{1})=\hat{\beta}_{\text{s}}(\eta)$, with

Lithium intercalation and de-intercalation were modelled only in the negative electrode. In the positive electrode we assume $\beta(C_{1})=0$ and $\beta_{\text{s}}(C_{1})=0$, motivated by experimental observations[12] that the positive electrode has insignificant swelling. While $P_{\text{l}}$ and $P_{\text{b}}$ are not evaluated numerically in our computations, we assume that $P_{\text{l}}\ll\kappa_{\text{s}}$, and $P_{\text{b}}\ll\kappa_{\text{s}}$. Therefore, the actual values chosen for $P_{\text{l}}$ and $P_{\text{b}}$ do not influence the porosity strongly, and these parameters were set to zero in the numerical examples. All other parameter values are listed in the Appendix.

With $\beta_{\text{s}}=0,\beta=0$, there is no expansion/contraction of the particle and composite electrode due to lithium intercalation/de-intercalation. For $\beta_{\text{s}}\neq 0,\beta=0$, the particle expands with Li intercalation; however, the composite electrode does not itself expand but accommodates the expanding particle in the microstructure’s free volume. For $\beta_{\text{s}}\neq 0,\beta\neq 0$, both the particle and composite electrode expand with Li intercalation. Figures 7 and 8, respectively, show the porosity distribution through the thickness of the electrodes and separator (along the $x_{2}$ coordinate) at the indicated state of charge for the 10 C rate, under isothermal and non-isothermal (governed by the heat equation) conditions. The compressive boundary conditions cause porosity to decrease rapidly in the separator as the free volume is compressed, while the porosity decrease in the stiffer electrodes is relatively small. This is seen in the black (initial) and purple (very short times after start of the computation) curves of Figure 7. The porosity has fallen further to $\epsilon_{\text{l}}\sim 0.32$ by SOC = 91% . During discharge, lithium undergoes de-intercalation from the solid particles in the negative electrode, causing its porosity to increase by contraction. When the cell is fully discharged, the porosity has increased to about $\epsilon_{\text{l}}=0.35$ in the negative electrode. We note that the negative electrode also undergoes mechanical contraction due to the de-intercalation. When fully discharged, the volume of the negative electrode therefore contracts by about 9$\%$. This results in a smaller increase in porosity compared with the case of no expansion of the composite electrode ($\beta_{\text{s}}\neq 0,\beta=0$) for which $\epsilon_{\text{l}}\sim 0.36$. We assume there is no intercalation strain in the positive electrode as justified above (see Section 4.2). Therefore its deformation is only due to stress transmitted from the negative electrode’s expansion/contraction, and due to its own thermal strain. However, the small thermal expansion results in insignificant dimensional change, and since the separator is one order of magnitude less stiff than the electrodes, the stress state induced under the applied displacement boundary condition causes insignificant strain in the positive electrode. As a result, the porosity of the positive electrode remains fairly constant.

From the computations, we note that temperature significantly affects the coefficients of the electro-chemical equations; i.e., the transport, conductivity and reaction parameters. Under isothermal conditions, the lithium and lithium ion concentrations are much more non-uniform in the neighborhood of the separator, as seen in Figures 9 and 10, and the porosity is also highly non-uniform by the end of the cycle, as seen in Figure 7. With evolving porosity, the changes in solid particle volume fraction, $\epsilon_{\text{s}}$, during discharging/charging exaggerate the nonlinearity of the reaction rate. The reaction is further confined to the region adjacent to the separator at high currents due to the reduced porosity. The non-linear and highly coupled set of equations (2-9) which describe the transport and reaction of lithium ions across the cell remains a challenge for numerical solvers as discussed by Ramadesigan et al.[33] Our simulations failed to converge after a few time steps at high current rates (10 C), as the local value of $C_{1}\rightarrow C_{1}^{\text{max}}$ in the negative electrode during charging, and the open circuit potential $U$ approaches zero (see Equations (7-9) and Figure 2). In this limit, the surface over-potential $\phi_{\text{S}}-\phi_{\text{E}}-U$ takes on large positive values. Consequently, the positive exponential term in Equation (7) grows, leading to high rates in the Butler-Volmer equation. Combined with the nonlinearly coupled equations, this results in a stiffness due to which the discretized Jacobian matrix is not positive definite even with very small time steps, and the direct solver (UMFPACK) fails. Consequently, we report that at fixed temperature, the computation with constant porosity can be advanced up to SOC = 76% during charging, which is greater than the maximum SOC of the cases with evolving porosity, as shown in Figure 9. The above trends in the negative electrode are reversed during discharging: $C_{1}\rightarrow 0$, $U$ approaches its maximum, $\phi_{\text{S}}-\phi_{\text{E}}-U$ takes on large negative values, and the negative exponential term in term in Equation (7) grows, leading to high rates in the Butler-Volmer equation. In this regime also, the direct solver fails and the numerical solution did not converge after some time.

During discharge, the contracting electrode undergoes stress relaxation. Stress equilibrium under the applied displacement boundary condition (Table 1) then causes the separator to expand. Thereby, its porosity increases from $\epsilon_{\text{l}}=0.62$ to $\sim 0.64$ as seen in the blue curves in Figure 11. This porosity gain is surrendered during charging.

As explained above, the high concentration of lithium forces a numerical divergence, which we trace to the growth of exponential terms in the Butler-Volmer model as $C_{1}\rightarrow C_{1}^{\text{max}}$ during charging, and numerical stiffness of the system of equations causing the direct solver to fail. This divergence happens before the cell is fully charged to its initial state from which discharging begins. However, this holds only for computations in the isothermal case. During operation, the battery heats up to $\sim 45^{\circ}$ C at the 10 C rate (Figure 12), leading to higher diffusivities and a more uniform distribution of lithium concentration as seen by comparing Figure 9 with Figure 13, as well as more uniform lithium ion concentrations, as seen by comparing Figure 10 and Figure 14. The Butler-Volmer model does not cause a numerical divergence of the solution before being fully charged back to its initial state under these conditions. The porosity also is largely recovered, although at a slightly more non-uniform distribution than before discharge (Figure 8). From the computations we note that the evolving porosity most strongly affects the potential in the electrolyte at states that are close to fully discharged, as seen in Figures 15 and 16. Figure 17 shows the solid phase potential profile where porosity changes show little effect. However, when the porosity undergoes large changes, such as if the function $\beta_{\text{s}}(C_{1}/C^{\text{max}}_{1})\sim 2$ as has been assumed for the non-isothermal case computed in Figure 18, it has a strong influence on the potential profile. Although such large values of $\beta_{\text{s}}$ are not the focus of this communication, this scenario may be instructive for studying materials such as tin oxide which undergoes $250\%$ volume expansion[19]. The non-symmetric potential profile with respect to discharging/charging in Figure 18 reflects the fundamental irreversibilities in battery operation, here arising due to the transport-reaction and heat conduction/generation phenomena.

During operation, the electrodes deform due to lithium intercalation, thermal expansion and external traction. Since the stress induced by cell deformation is uniform through the cell thickness in these effectively one-dimensional phenomena, we have shown the surface traction force evolving with time in Figure 19. The displacement boundary condition imposes an initial compressive stress on the cell. This compressive stress relaxes since the electrode contracts during discharging, and increases again as the electrode swells due to intercalation during charging. Figure 19 also shows that the force induced by lithium intercalation is large, but that generated by thermal expansion is small and evolves steadily.

The same initial and boundary value problem was also run at the 1 C current rate under non-isothermal conditions (Figures 20-23). Under this low current rate, however, thermal effects are insignificant, the kinetics are slow and thermodynamic dissipation rates are also low. Consequently the battery’s performance is much more symmetric in a discharging$\to$charging cycle. Additionally, porosity evolution is uniform, as shown in Figure 20, and the lithium concentration is also uniform as shown in Figure 21. The porosity evolution in the separator in Figure 22 and the solid phase potential in Figure 23 reflect this symmetry.

The information, data, or work presented herein was funded in part by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.