Intercalation kinetics in multiphase layered materials
Raymond B. Smith, Edwin Khoo, Martin Z. Bazant
TL;DR
This paper develops a phase-field model to simulate intercalation kinetics in layered materials like graphite, capturing staging dynamics and voltage profiles in lithium-ion batteries, highlighting the tradeoff between detailed and simplified modeling approaches.
Contribution
It introduces a thermodynamically consistent phase-field model for layered intercalation materials, including a reduced model for computational efficiency and better experimental data fitting.
Findings
Model captures voltage staircase due to staging dynamics
Explicit layered models reveal internal domain structures
Reduced model efficiently fits experimental voltage profiles
Abstract
Many intercalation compounds possess layered structures or inter-penetrating lattices that enable phase separation into three or more stable phases, or "stages," driven by competing intra-layer and inter-layer forces. While these structures are often well characterized in equilibrium, their effects on intercalation kinetics and transport far from equilibrium are typically neglected or approximated by empirical solid solution models. Here, we formulate a general phase-field model with thermodynamically consistent reaction kinetics and cooperative transport to capture the dynamics of intercalation in layered materials. As an important case for Li-ion batteries, we model single particles of lithium intercalated graphite as having a periodic two-layer structure with three stable phases, corresponding to zero, one, or two layers full of lithium. The electrochemical intercalation reaction is…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Intercalation kinetics in multiphase layered materials
Raymond B. Smith
Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Edwin Khoo
Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Martin Z. Bazant
Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Department of Mathematics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Abstract
Many intercalation compounds possess layered structures or inter-penetrating lattices that enable phase separation into three or more stable phases, or “stages,” driven by competing intra-layer and inter-layer forces. While these structures are often well characterized in equilibrium, their effects on intercalation kinetics and transport far from equilibrium are typically neglected or approximated by empirical solid solution models. Here, we formulate a general phase-field model with thermodynamically consistent reaction kinetics and cooperative transport to capture the dynamics of intercalation in layered materials. As an important case for Li-ion batteries, we model single particles of lithium intercalated graphite as having a periodic two-layer structure with three stable phases, corresponding to zero, one, or two layers full of lithium. The electrochemical intercalation reaction is described by a generalized Butler-Volmer equation with thermodynamic factors to account for the flexible structure of the graphene planes. The model naturally captures the “voltage staircase” discharge curves as a result of staging dynamics with internal “checkerboard” domains, which cannot be described by solid-solution models based on Fickian diffusion. On the other hand, the two-layer model is computationally expensive and excludes low-density stable phases with longer-range periodicity, so we also present a reduced model for graphite, which captures the high-density stages while fitting the low-density voltage profile as an effective solid solution. The two models illustrate the general tradeoff between the explicit modeling of periodic layers or lattices and the needs for computational efficiency and accurate fitting of experimental data.
I Introduction
Staging, the phenomenon of pattern formation by intercalants or defects across a layered structure, occurs in a wide range of materials including high temperature superconductors Wells et al. (1996); Mohottala et al. (2006), clays and layered double hydroxides Ijdo and Pinnavaia (1998); Iyi et al. (2002); Pisson et al. (2003), MXenes Naguib et al. (2012); Wang et al. (2015), metal borocarbides Joshi et al. (2015), layered metal oxides Van der Ven et al. (1998); Chen et al. (2002), and the prototypical case of graphite Dresselhaus and Dresselhaus (1981). The properties of these materials can depend strongly on the microscopic staging structure. In high temperature superconductors, the layered segregation of oxygen can lead to regions of insulating behavior Mohottala et al. (2006). Staging phase-separation phenomena can also arise in three-dimensional materials with interpenetrating lattices of intercalation sites, as with hydrogen trapping in metals Di Leo and Anand (2013). Another example could be lithium intercalation in anatase (TiO2) de Klerk et al. (2016), which has been suggested to have multiple types of intercalation sites Lafont et al. (2010); Shin et al. (2011); Shen et al. (2014) and demonstrates staircase voltage curves Zachau-Christiansen et al. (1988), similar to graphite Ohzuku et al. (1993).
In graphite, the staging process is common for many intercalants, including lithium which forms LiC6, the most common anode material in Li-ion batteries. The staging of the lithium-graphite system is directly related to the observed voltage in batteries. The dynamics of staging during (de)lithatiation of graphite is important in consideration of Li-ion battery rate capability, because graphite is (de)intercalated with Li during (dis)charging. Li-ion batteries are ubiquitous as energy storage devices for portable electronics, with growing use in the (hybrid) electric vehicle market. In almost all commercial cells, graphite is used as the negative electrode (anode) Scrosati and Garche (2010). Li-ion batteries work by shuttling Li ions back and forth between the two electrodes Tarascon and Armand (2001), a process which involves transport within the solid and electrolyte phases as well as electrochemical reaction to transfer the lithium between these phases. The rates of these processes (and those of competing processes) fundamentally govern the power capabilities of the batteries.
Complicating our understanding of these processes is the fact that graphite tends to phase separate into Li-rich and Li-poor regions Dahn (1991), which affects both the reaction and the transport processes Bazant (2013). As Li is transferred into graphite, it intercalates between graphene planes, tending to form structures with filled layers separated by a number of empty layers. These structures are referred to by a stage number, which refers to the number of planes of graphene between filled layers. For example, stage 1 graphite is nearly full (Figure 1 (c)), whereas high stage number graphite has many empty layers between each full layer (e.g. Figure 1 (b) for stage 2).
In the dilute limit, there is no internal structure, which is referred to as stage 1’ Dresselhaus and Dresselhaus (1981) (Figure 1 (a)). Upon formation of the low-stage number structures (stages 1’, 2, and 1), graphite undergoes distinct color transformations, which allows direct visual probing of the local state of charge of graphite electrodes Maire et al. (2008); Harris et al. (2010). In addition, theoretical Daumas and Herold (1969); Safran (1980a); Kirczenow (1982); Krishnan et al. (2013) and experimental Clarke et al. (1979); Dimiev et al. (2013) studies suggest intercalants within graphite tend to phase separate within layered planes such that internal domains are formed (Figure 1 (d)).
Phase separation is common in many battery electrode materials such as lithium iron phosphate Padhi et al. (1997), spinel manganese oxides Ariyoshi et al. (2004), and lithium titanate Ohzuku et al. (1995), and it can dramatically affect the internal dynamics of both single particles Lim et al. (2016); Tang et al. (2010); Bai et al. (2011); Cogswell and Bazant (2012, 2013); Zeng and Bazant (2014); Welland et al. (2015) and porous electrodes Dreyer et al. (2010); Chueh et al. (2013); Li et al. (2014); Ferguson and Bazant (2012, 2014); Orvananos et al. (2014, 2015). As such, in order to develop predictive models of the behavior of such materials with complex phase behavior, it is critical that these models naturally capture the phase separating dynamics. Some models of phase separating electrode materials approximate them as solid solutions described by Fickian diffusion Fuller et al. (1994); Verbrugge and Koch (2003); Bernardi and Go (2011); Safari and Delacourt (2011); Srinivasan and Newman (2004a); Baker and Verbrugge (2012). Others approximate the phase separation behavior using artificially imposed equilibrium interfaces, which move via Stefan conservation conditions and use the equilibrium voltage curve to determine the functional form of the chemical potential in the solid Funabiki et al. (1999); Srinivasan and Newman (2004b); Baker and Verbrugge (2012); Heß and Novák (2013). Others have attempted to describe the volume-averaged phase behavior using nucleation models Gallagher et al. (2012). However, such approaches cannot capture the formation of interfaces by spinodal decomposition or other processes of non-equilibrium thermodynamics.
Phase-field models, pioneered by Van der Waals Rowlinson (1979) and Cahn and Hilliard Cahn and Hilliard (1958); Cahn (1961), have been widely applied to phase separating systems to describe the dynamics of their approach to equilibrium states. This approach is based on linear irreversible thermodynamics Groot and Mazur (1962), has been extended to electrical systems Garcıá et al. (2004), and has been self-consistently connected to electrochemical reactions to describe phase transformations in chemical and electrochemical systems Bazant (2013). Here, focusing on staging of lithiated graphite, we implement a Cahn-Hilliard style model in a multi-layer free energy framework to capture the phase behavior of graphite systems, similar to the work of Hawrylak and Subbaswamy Hawrylak and Subbaswamy (1984). We also directly relate the phase separation behavior to electrochemical kinetics via a thermodynamically consistent reaction boundary condition Bazant (2013), which has recently been validated by in operando scanning X-ray transmission microscopy of individual lithium iron phosphate nanoparticles Lim et al. (2016). A previous work Guo et al. (2016) demonstrates good agreement of the two-layer graphite model with an experiment visually observing the phase separation dynamics during a single particle discharge. In a companion paper Thomas-Alyea et al. (2016), we apply a one-variable reduced model, also presented below, in porous electrode simulations to capture electrochemical and optical measurements on a porous graphite battery electrode. The porous electrode experiments and modeling suggest that the modifications of the basic framework presented here are important to describe the non-equilibrium thermodynamics of Li-ion batteries, as we will discuss in Section V. In this paper, however, we focus on the development and predictions of the general model, which could be applied to other intercalation materials exhibiting staging associated with layered crystal structures or inter-penetrating lattices.
II General Theory
Layered intercalation materials tend to exhibit stable phases with periodicity in the transverse direction, which could be described by multiple concentration variables. We assume the layers of the material are stacked and perfectly overlapping, i.e. they look like those in Figure 1 with all layers present throughout the particle and to the edges. Then, we describe the concentration of intercalated species within each layer in only 2 dimensions, effectively depth averaging along the thickness within each individual layer. This allows us to represent the 3-dimensional layered structure as discrete, overlapping, 2-dimensional planes. Thus, we will present a 2-dimensional model with multiple concentration variables corresponding to species in different layers. As noted above, the same model could also apply to materials with multiple interpenetrating sub-lattices, represented by separate concentration variables.
II.1 Multi-Layer Free Energy Model
Mathematically, we begin by defining the total free energy of the system, , as an integral of the local free energy density, , over the system volume, ,
[TABLE]
The choice for determines the physical description of the system. We will focus on systems with two-dimensional symmetry in which multiple concentration variables correspond to species occupying different vertically stacked layers. However, the treatment is analogous to descriptions of three-dimensional materials in which intercalated species can occupy different types of sites, such as hydrogen trapping in metals Di Leo and Anand (2013) or lithium intercalation in anatase Zachau-Christiansen et al. (1988); Lafont et al. (2010); Shin et al. (2011); Shen et al. (2014); de Klerk et al. (2016).
For a general system with layers or site types, each occupying some part of the total system volume, ,
[TABLE]
For an -layer structure, in which each 2-dimensional layer has area and occupies a fraction, , of the height of the unit, , the total free energy is
[TABLE]
Next, we assume individual layer free energy densities are defined by some function, which depends only on the configuration of layer , as well as a series expansion of mixing enthalpies of particle clusters interacting between layers. Thus
[TABLE]
The reference concentration, , can be chosen to be some suitable scale, the local scaled concentrations are denoted by , and the parameters define the enthalpic interaction energies of multi-body particle clusters.
II.2 Intra-layer Free Energy
We first choose a model for the intra-layer free energy, , which will not depend on the layer, assuming that each region between adjacent planes is structurally similar. As the simplest model that allows phase separation within layers, we consider a regular solution, a lattice gas model Derosa and Balbuena (1999); Ledovskikh and Wagemaker (2016), which describes the entropic configurations of particles and vacancies on a lattice along with an enthalpic energy of mixing between particles and vacancies. We add a gradient penalty term to capture the effect of interfacial energy between the intra-layer domains Cahn and Hilliard (1958),
[TABLE]
where is Boltzmann’s constant, is the absolute temperature, is the regular solution parameter for interaction of particles and vacancies within a layer, is the reference chemical potential, and is the gradient energy penalty, assumed here to be isotropic within each plane. Of note, we use the gradient, , and divergence, , operators here to indicate their evaluation only within the plane of a single layer, and the non-dimensional concentrations are scaled such that they are filling fractions of individual layers. When , this model defines a free energy curve with two minima near the extreme filling fractions, leading to a tendency for internal phase separation to regions with high and low concentrations. The function is plotted in Figure 2 using thermodynamic parameters from Section III, at , choosing , and neglecting contributions from the gradient penalty term.
II.3 Periodic Bilayer Model Applied to Graphite
The simplest example of the multilayer framework describes a periodic solid with repeated pairs of identical layers for intercalation, which is a natural approximation for lithiated graphite that captures the major staging transitions. A bilayer model allows description of stage 1, 2, and 1’ structures, neglecting higher stage number structures in favor of model simplicity and reduced computational cost. Each modeled layer represents a plane of Li sites between graphene planes in (de)intercalated graphite, and the two modeled layers repeat in the direction normal to the modeled planes. Thus, we have a discretely “stacked” two-dimensional continuum description of the intercalant within the graphite. We do not assign any clear physical separation between layers in the vertical direction other than through energetic interactions, which act as local approximations of the misfit strain energy in graphite Schiffer et al. (2016), without accounting for nonlocal coherency strain Cogswell and Bazant (2012). From the choice of structurally similar layers, each of which occupies half the total volume of the system, the reference concentration, , where is the concentration of Li in LiC6.
For the inter-layer interactions, we propose a model incorporating 2- and 4-body interactions,
[TABLE]
where subscripts 1 and 2 represent the repeated layers 1 and 2 respectively. For graphite, the parameter is related to a repulsive interaction between cross-plane lithium Safran and Hamann (1979), and it allows for prediction of the “staircase” voltage plateau measured in near-equilibrium (dis)charge curves for graphite electrodes Ohzuku et al. (1993); Ferguson and Bazant (2014) by increasing the free energy of the full stage 1 state compared to that of the stage 1’ and stage 2 states (see Figure 3). The term can be interpreted as a particle-vacancy cross-plane mixing enthalpy, which Ferguson and Bazant found to be necessary to push particles toward stage 2 in the filling process rather than along the direction, although we will reexamine this in Sections IV.2.2 and IV.2.3. The primary change of this free energy model beyond that of Ferguson and Bazant Ferguson and Bazant (2014) is the gradient penalty in Eq. 5, which allows us to predict intra-layer phase separation. The thermodynamic treatment is similar to that in previous work describing staging in graphite Hawrylak and Subbaswamy (1984); Safran (1980b); Kirczenow (1984). For a review of other theoretical work on graphite staging, the reader is directed to the informative ref. Kirczenow (1990).
Here, it is helpful to examine the basic thermodynamic structure of the model, which is presented in Figures 2 and 3. Using the thermodynamic parameters from Section III (and choosing ), the free energy density of the system, , is plotted in Figure 3 (a) and (b) as a function of the filling fractions in each layer, and , assuming minimal energetic contribution from concentration gradients. Because the predicted free energy density is not a convex function, the system will tend to phase separate. For a simpler, single-variable free energy model such as the regular solution model applied to individual layers (Figure 2), the equilibrium phase separating behavior can be described by creating a common tangent line between the energy wells. When the concentration of a homogeneous system lies between the two intersection points, the system can lower its free energy to the common tangent by separating into two phases, each of which has equilibrium concentration set by the intersection points. The overall ratio of phase amounts is set by the lever rule. For the two-variable free energy model developed here, we instead construct common tangent planes, pictured in green in Figure 3 (a). Similar to the simpler case, the intersections of these planes with the homogeneous free energy surface, depicted as red dots, correspond to the stable equilibrium concentrations of the phases that will develop from overall system concentrations lying within the triangles formed by the three red dots of a given plane. This predicts that the equilibrium concentrations of stage 2 phases depend on whether the system is all in stage 2, in equilibrium with state 1’ (lower plane), or in equilibrium with stage 1 (upper plane).
In Figure 3 (a) and (b), three paths are plotted, which correspond to selected filling/emptying paths as the graphite is (de)intercalated. The orange curve represents a fully homogeneous (and energetically unfavorable) filling path in which . The red curve represents the minimum energy path for filling with individual layers remaining fully homogeneous (no intra-layer phase separation). The blue curve represents an equilibrium path allowing intra-layer phase separation. Following these paths motivates the definition of a quantity which we will refer to as an effective chemical potential,
[TABLE]
This quantity represents the change in free energy along the given filling path, , to insert a particle in either of the layers; it is dependent on the state of the system. is related to predicted equilibrium voltage curves Ferguson and Bazant (2014), as we will see in Section IV.1.1. Examples of its values over the three selected filling paths are presented in Figure 3 (c).
Returning to the formulation of the model, we now focus on the kinetic processes. The interstitial lithium flux, , occurs within each modeled plane of sites and is driven by gradients in the diffusional chemical potential Groot and Mazur (1962), , defined as the variational derivative of the free energy with respect to the concentration within the layer Cahn (1962),
[TABLE]
[TABLE]
where .
II.4 Concentrated Solid Diffusion
From the assumptions of linear irreversible thermodynamics Groot and Mazur (1962), the flux can be expressed as
[TABLE]
where is the tracer diffusivity, related to the mobility by the Einstein relation, . Whereas Hawrylak and Subbaswamy chose a constant flux prefactor investigating a similar model Hawrylak and Subbaswamy (1984), we propose a flux relationship describing diffusion on a lattice Bazant (2013). First, we begin with tracer diffusivity given by
[TABLE]
where is the tracer diffusivity in the dilute limit, is the activity coefficient of the diffusing species, and is that of the diffusing transition state. Following the assumptions outlined in Appendix A, we find
[TABLE]
Here, we have assumed that the transition state only differs thermodynamically from lithium in a stable lattice state in the entropic contribution. However, the enthalpic terms could differ, and extra, local stresses could also affect the energetics of the diffusing transition state via an activation volumeAziz et al. (1991); Bazant (2013).
Using these flux relations, we can solve for the time and space evolution of concentration profiles within the layers using species conservation,
[TABLE]
where is a volumetric source/sink term, which could describe, e.g., inter-layer reaction in which Li can move from layer 1 to layer 2 in the bulk via defects (Li hopping through pristine graphene is unlikely because of large energetic barriers Meunier et al. (2002); Yao et al. (2012)). We will set here except as noted in Section IV.2.4.
II.5 Thermodynamically Consistent Reaction Boundary Conditions
The system is closed by specifying boundary conditions on the spatially 4th order partial differential equation for each layer, according to the Cahn-Hilliard reaction model for heterogenous reactions Bazant (2013); Zeng and Bazant (2014). The concentration at the surface of the particle is governed by the natural boundary condition Cahn (1977); Bazant (2013), \widehat{\textbf{n}}\cdot\mathopen{}\mathclose{{}\left(\kappa\boldsymbol{\nabla}\widetilde{c}_{i}}\right)_{S}=\frac{\partial\gamma_{S}}{\partial\widetilde{c}_{i}}, where is an outward facing unit normal vector, the subscript indicates evaluation at the surface, and is the surface energy. Although surface “wetting” has been shown to be important in nanoparticle dynamics Cogswell and Bazant (2013), we set \widehat{\textbf{n}}\cdot\mathopen{}\mathclose{{}\left(\kappa\boldsymbol{\nabla}\widetilde{c}_{i}}\right)_{S}=0 here. The remaining boundary condition is a flux specification at the surface, which is determined by the dynamic process being examined. For example, in order to study spinodal decomposition, no flux boundary conditions can be imposed Cahn (1962). Here, because we are interested in examining battery operation, we specify the total current, , into the particle as the sum of integrated current densities, , into the perimeter of the particle ,
[TABLE]
where is the height of the particle and are the fractions of the surface area occupied by each layer. The current density into each layer is related to the flux from the electrochemical reaction associated with intercalation into each of the repeating layers,
[TABLE]
We describe the reaction using a thermodynamically consistent modification Bazant (2013) to the Butler-Volmer equation Bard and Faulkner (2001), which depends on the diffusional electrochemical potentials of the oxidized () and reduced (Li) states as well as the activity coefficient of the reaction transition state,
[TABLE]
where subscripts indicate variables corresponding to layer , is the reaction rate constant defined per unit reaction area of the particle, is the elementary charge, is the activity of the oxidized state (assumed to be unity here), is the activity of the reduced state in layer (), is the reaction symmetry factor, is the activity coefficient of the transition state for reaction into/out of layer , and is the activation overpotential, or reaction driving force, for layer (). is the diffusional electrochemical potential of the oxidized state,
[TABLE]
where is the reference chemical potential in the oxidized state, is the electric potential in the electrolyte, and is the electric potential in the graphite. Thus, we see that the total integrated current will be directly related to the applied interfacial potential, , which will allow simulations with either specified current or specified voltage. We note that use of the Butler-Volmer equation ignores recent evidence of departures from its predictions, especially at large driving forces Chidsey (1991); Bai and Bazant (2014); Laborda et al. (2013). Here, it is convenient to group the standard potentials, , to be the standard half-cell potential relative to the Li/Li+ metal standard half-cell potential.
The reaction transition state activity coefficient is used to describe physics of the reaction process. For example, for Li intercalation into the rigid ion channels of the positive electrode material FePO4 (iron phosphate), Bai et al. assumed the reaction transition state requires (excludes) a single empty lattice site, such that Bai et al. (2011). To describe reaction at the surface single crystal of graphite, ref. Guo et al. (2016) used
[TABLE]
which is based on the assumption that a reaction event in layer requires both a vacant site and a nearby filled site for the reaction to occur. This requires that an empty site be present near the surface for intercalation. In addition, when there is no nearby filled site, the graphene layers can collapse, making the reaction unlikely. Thus Eq. 18 predicts that the reaction resistance for the layer, which scales with , diverges at empty and full filling fractions, as presented in Figure 4.
This implicitly assumes that the reaction transition state behaves like an ideal solution except for entropic exclusion effects, unlike the diffusion transition state (Eq. 31), which we assume also depends on all the enthalpic terms of the diffusional chemical potential because it occurs in the bulk of the system rather than at the interface between one phase and another. We use this reaction model here except where noted to demonstrate some of its predictions. However, because of an inability to carefully experimentally control for other kinetic limitations in ref. Guo et al. (2016), more detailed consideration of this reaction model would be worthwhile. For example, Thomas-Alyea et al. Thomas-Alyea et al. (2016) find when describing secondary graphite particles with a variant of this model (developed in Section V) that a reaction model dominated by film resistance enabled much better agreement with the experimental results. Thus, we vary this model where noted below.
III Model Implementation
The energetic parameters in the free energy are taken as those fit to thermodynamic data in previous work Ferguson and Bazant (2014): , , and with , and the system is simulated at room temperature, . Here, we have adjusted their fit value of from 0.1366 V to 0.12 V in order to account for the fact that our particles can phase separate within the layers, which changes the observed voltage plateau. The difference corresponds to the dimensional length of the black dashed line in Figure 3 (c) as explained previously Levi and Aurbach (1997a); Dreyer et al. (2010) and demonstrated for a different material model in the same paper Ferguson and Bazant (2014). With the intra-layer phase separating particles, the value of corresponds to the voltage of the low-filling-fraction plateau, which we find to be 0.12 V from near-open-circuit experimental data Ohzuku et al. (1993). We choose to match the experimental and theoretically calculated density of lithiated graphite Wissler (2006); Kganyago and Ngoepe (2003).
The values of and are taken from ref. Guo et al. (2016) in which we used this model and found excellent agreement with optical measurements of the concentration profiles of a single graphite crystal in time and space. From this, and J/m, although the fit was relatively insensitive to the choice of , and the value of agrees well with ab initio calculations Persson et al. (2010). The reaction symmetry factor, , is set to 0.5, which implies a symmetric electrochemical reaction Bard and Faulkner (2001). Because the single particle electrode experiment did not allow unique determination of the reaction rate constant, , we do not take the value used in ref. Guo et al. (2016) and instead vary this parameter here as stated.
Inspired by our previous work Guo et al. (2016), we simulate a cylindrical particle in which intercalation occurs around the circumference of the particle into the circular layers. For simplicity, we assume axial symmetry about an axis normal to the simulated plane in the center of the particle, allowing us to reduce the model to solving in time, , and the radial coordinate, . We use a mass-conserving variant of finite volume discretization in space following ref. Zeng and Bazant (2014), ensuring that the grid spacing is less than the scale for the interfacial width, \lambda_{b}\sim\mathopen{}\mathclose{{}\left(\kappa/\mathopen{}\mathclose{{}\left(c_{\mathrm{ref}}\Omega_{a}}\right)}\right)^{1/2}=5.8\times 10^{-8}\ \mathrm{m}. Time integration is carried out using DAE Tools Nikolić (2016), which wraps the SUNDIALS integration suite Hindmarsh et al. (2005) with the ADOL-C automatic differentiation library Griewank et al. (1996) to facilitate the non-linear solver involved in the implicit time stepping. Single particle simulations took between a few minutes and a few hours on one CPU. In the following two sections, we validate the simulation by implementing the same model in COMSOL Multiphysics and comparing outputs, then using the COMSOL Multiphysics implementation in 2D to examine the validity of the assumption of axial symmetry in context of the experiment simulations from ref. Guo et al. (2016).
III.1 Simulation Validation
Beginning from the simulations in ref. Guo et al. (2016), we first validate the 1D finite volume simulation by implementing the same model in COMSOL Multiphysics, which uses the finite element method, to verify that both methods give similar results. In each case, we use the parameters given in ref. Guo et al. (2016): , , , with respect to Li/Li+. The low value of could be thought of as adjusting the value of , as the two compensate for each other (changing , , and by the same factor leads to the same simulation outputs). Because the results in ref. Guo et al. (2016) were not sensitive to the value of , we simply use their reported value of directly here with the value for in Section III for the remainder of this work. Using a grid spacing of , both simulations reproduce the previous results and agree with each other in Figure 5. We note that the model is symmetric in and . Therefore, in both the COMSOL and finite volume simulations, and profiles can be flipped, and the occurrence of this flipping is random and depends on numerical noise. Full movies are provided in the supplement.
III.2 Axial Symmetry Assumption
We examine the assumption of axial symmetry by relaxing it and implementing the resulting 2D model in COMSOL Multiphysics, again using parameters from ref. Guo et al. (2016). We choose the simulation domain to be the largest possible 2D slice of a circular disk such that the simulation finishes within a few hours on one CPU. This periodic simulation domain is replicated as many times as necessary to give the full circular disk. With this constraint, we use a sector with a central angle of and a grid spacing of . Taking a slice imposes an extra non-physical length scale on the simulation, related to the arc length of the simulated slice, and small enough slices may prevent formation of phase interfaces with normal components in the -direction and artificially encourage axial symmetry. Our simulated arc length of is several times the interface width for this simulation, so we expect this not to prevent breaking of axial symmetry at least near the particle surface. We see in Figure 6 that the axial symmetry assumption is broken at early times before complete phase separation occurs, but as the concentration front propagates inward from the particle surface, the particle regains axial symmetry and profiles become nearly identical to that predicted by the 1D model. In Appendix B, we consider a more detailed comparison of the 2D concentration profiles with a simulation of the same system in 1D. Movies are also provided in the supplement. Overall, at early times, the 1D and 2D cases show moderate differences, but at later times, the simulations become quite similar, so these simulations support the validity of the axial symmetry assumption.
IV Results
Because the general intercalation framework is consistently connected to electrochemical reactions, we are able to investigate typical experimental procedures such as constant current and constant voltage (dis)charge. We present here a few example predictions made by the model to highlight some of its key features, with focus on those which deviate strongly from models based on Fickian diffusion. Given initial conditions and a specified current (voltage) profile, the model predicts the concentration profiles and as well as the output voltage (current) profile. Here, we also define the local average filling fraction,
[TABLE]
which we will relate to the stage number.
In the following sections, we explore a number of model predictions. In Section IV.1, we focus on (near-)equilibrium properties of the model and relate them to the free energy model by examining the open circuit voltage of a single particle, the equilibrium structures and concentrations, and the phase diagram predicted by the thermodynamic model. In Section IV.2, we examine various predictions the model makes about transient processes. First, in IV.2.1, we examine a constant voltage intercalation process, which highlights the overall phase behavior the model predicts, including propagating phase boundaries and diffusive profiles in the stable phases between phase boundaries, as well as the formation of checkerboard patterns during the filling process. In Section IV.2.2, we study the process of spinodal decomposition of the particle from an initially homogeneous state to stable phases. This again highlights the checkerboard patterns that emerge because of the free energy structure of this model which could not be captured by some other models which are commonly applied to phase separating systems. In Section IV.2.3, we examine how the single particle behaves under a constant current discharge with particular attention given to the impact of the reaction rate model on the predicted system voltage. We look at discharge voltage curves both with the reaction model presented in Eq. 16 and an alternate version to highlight the importance of this component of the model in the predictions it makes. We show that the model in Eq. 16 (used in ref. Guo et al. (2016)) does not give predictions qualitatively matching experimental graphite discharge curves, whereas a simple alternative choice works better. We also give some discussion about the significance of this observation in context of applying the model to battery electrodes as in ref. Thomas-Alyea et al. (2016). In Section IV.2.4, we repeat the constant current discharge but allow for bulk inter-layer exchange of lithium and note the differences. This exchange is not relevant for pristine graphite, but demonstrates qualitative differences in predictions for systems in which exchange between lattices is possible, including graphite with defects Yao et al. (2012). We conclude the results with Section IV.3, in which we make some comparisons of this model with a model based on Fickian diffusion (solid solution) with a voltage curve fit to the predicted open circuit voltage (OCV) of the layered phase separating model. We choose to compare the model predictions of the solid solution and phase separating models when exposed to a current pulse within a two-phase plateau, where the models give particularly different predictions of both concentration profiles and predicted voltage.
IV.1 Equilibrium Behavior
IV.1.1 Open Circuit Voltage
First, we simulate a single particle with radius under a constant current discharge at a rate of C/10,000, (where C corresponds to the rate to (dis)charge the particle in hours). With a rate constant of , the characteristic reaction time is given by
[TABLE]
and the characteristic time for species transport is
[TABLE]
We note that another time scale could be formed for transport over the length of a phase boundary, , but that is a significantly shorter time scale. Thus, because the imposed discharge process time scale (10,000 h) is significantly longer than any inherent time scale in the system, we expect the system to be near equilibrium during the entire process. Because the system remains near equilibrium, the system voltage should be determined by a Nernstian relationship,
[TABLE]
where is defined in Eq. 7 and corresponds to the blue paths in Figure 3 in which intra-layer phase separation is allowed. This is what we see in Figure 7, which also presents experimental data of a very slow graphite electrode discharge Ohzuku et al. (1993). As the simulation proceeds, lithium initially fills either layer 1 or layer 2. Each layer phase separates internally into high/low concentration regions, until the system is at an overall filling fraction of 0.5 and in a stage 2 structure everywhere. This process corresponds to moving along the lower common tangent plane in Figure 3 (a). Then, lithium proceeds to fill regions into a full, stage 1 structure, which corresponds to moving along the upper common tangent plane in Figure 3 (a). Thus, in this simulation, the system follows along common tangent planes, similar to the blue paths in Figure 3, which predicts the “staircase” behavior of the voltage curve. Of note, the initial “overshoot” until a filling fraction near is caused by the system entering a metastable region before falling to the stable-equilibrium plateau. A movie of the simulated concentration profiles is included in the supplement.
The disagreement between model and data in Figure 7 at low filling fraction is related to the model assumption of only having two repeating, structurally similar layers, whereas real graphite can form high stage number structures with different energies Dahn (1991). Thus, because each layer is treated identically for physical reasons, the first half and second half of the equilibrium discharge voltage curve should look translationally similar, and we are unable to capture the details of the low-filling voltage curve with a 2-layer model. As we have only implemented the model with two layers here, we can only conjecture that adding more layers to the model may enable us to capture the higher stage structures and also the shape of the low-filling open circuit voltage curve.
IV.1.2 Equilibrium Structures
As described in Section II, we can calculate equilibrium phase concentrations by making the homogeneous free energy density function convex using common tangent planes. Calculation of the planes was carried out numerically and provides values for the expected results of equilibrium structure calculations. In order to compare to results using the full model, we perform simulations to examine the simple cases of stage 1’ - stage 2 equilibrium (Figure 8 (a)) and stage 1 - stage 2 equilibrium (Figure 8 (b)). These equilibrium structures are calculated by initializing the system with uniform high/low concentration in one layer and a step function for concentration in the other and simulating with no current until stage 2 - stage 1 (1’) equilibrium is achieved. The transients of the process involve (1) forming the correct interface shape (with width ), (2) reaching the equilibrium concentrations in the phases, and (3) shifting the interface left/right according to the lever rule to preserve the total average concentration. Simulated final concentrations in each phase agree with those calculated directly from the common tangent construction, and show qualitative agreement with experimental images of lithiated graphite phase equilibrium Song et al. (1996).
IV.1.3 Phase Diagram
One major advantage of models constructed from simpler thermodynamic models like the regular solution is that they make predictions about temperature dependence, unlike solid solution models in which the chemical potential as a function of filling fraction at a particular temperature is directly fit to an open circuit voltage at that temperature. Using the free energy presented here, we can construct phase diagrams following a similar procedure used to add the common tangent planes in Figure 3. We first find common tangent planes (if any) of the free energy surface at a fixed temperature. These indicate the presence of inter-layer phase separation, as demonstrated in Figure 3. Then, in regions outside of the planes, common tangent lines are constructed along slices of the free energy at constant average filling fraction (e.g. see the stage 2 region in Figure 3). The bounds of the stage 2 regions correspond to either (a) entering regions of intra-layer phase separation (entering a region with a common tangent plane), or (b) the largest/smallest filling fractions with a non-convex free energy slice possessing a common tangent.
Doing this process over a range of temperatures, we obtain the diagrams in Figure 9, in which we present the result both for the thermodynamic model parametrized as in Section III and also for the same model neglecting the particle-vacancy mixing enthalpy (captured by the term), which leads to similar predictions as made by Safran Safran (1980a). Ferguson and Bazant Ferguson and Bazant (2014), who originally introduced the term, did not carefully fit this parameter, as its value had little influence on their results. Here we can see that although the phase diagrams are similar at room temperature, this term strongly influences the predicted equilibrium behavior as temperature increases. Curiously, for the case in Figure 9 (a) near , the common tangent planes intersect the free energy only at points off the line, leaving regions of stage 2 both at intermediate filling fractions and also at low and high filling fractions between the intra-layer phase separation and the fully homogeneous region.
Comparing these diagrams to that of the near-room temperature LiC6 system Dahn (1991); Fischer et al. (1983); Kirczenow (1990), we see again that the model cannot capture the details of the low average filling (high stage number) phases. Also, neither model captures the details of the top line of the phase diagram including the high-temperature stable stage 2 region Fischer and Kim (1988); Woo et al. (1983); Kirczenow (1990), although the case with more closely approximates the temperature values at which order disappears (top line of the diagram). Although neither model accurately captures the phase diagram, the result highlights the value of free energy based models as starting points to make predictive temperature-dependent models to capture electrochemical behavior in non-isothermal systems Reynier et al. (2004). Natural modifications like inclusion of stresses or more careful representation of higher order cluster terms in the free energy would enable a better representation of the phase diagram, and including more layers with longer-range interactions leads to a rich set of phase diagram predictions Millman and Kirczenow (1982, 1983). For consistency with Ferguson and Bazant Ferguson and Bazant (2014), we will use the model as presented in Section III, as the focus of this work is at room temperature where the differences between the two variants in Figure 9 are not significant. Importantly, this diagram is a reduction from the complete space for visual simplicity. For example, within the “stage 2” region, the system is actually phase separated with , and the model predicts the values of these equilibrium concentrations, but they are not shown here.
IV.2 Predicted Dynamics
IV.2.1 Constant Voltage Intercalation
Here we apply a constant voltage to a particle with radius and set . In this simulation, we choose with respect to lithium metal, leading to a strong driving force for intercalation. This corresponds to the applied voltage simulated in ref. Guo et al. (2016). Note that we are ignoring the competing reaction for lithium plating which is also thermodynamically favored under these conditions. To connect the concentration profiles to the microstructure and the associated visual particle colors Harris et al. (2010), we find it informative to depict the simulated concentration profiles on a cross-sectional slice of the simulated cylindrical particle. For example, we plot a snapshot of this simulation in Figure 10, demonstrating the basic features of the model under these operating conditions.
The concentration profiles on the bottom correspond to the layers of green/yellow between the black graphene planes in the particle slice above. In addition, the particle surface (top) is colored as it would be seen experimentally, according to the internal stage number Harris et al. (2010). This snapshot also demonstrates some key characteristics of the model. First, we see the anticipated behavior of phase separation within layers caused by the intra-layer regular solution model of Eq. 5. Second, in the intermediate region (), the system forms a stage 2 structure with lithium organized in alternating full/empty layers, caused by the repulsive energetic parameters in Eq. 6. Interestingly, in this stage 2 region, the simulation naturally predicts the internal “checkerboard” domains of stage 2 structures Krishnan et al. (2013); Dimiev et al. (2013). The overall concentration profile going from full near the particle edge to empty near the center is a result of transport limitations as the lithium is inserted from the edge of the particle. The gradual sloping of the concentration profiles within each stage region (including the half-full stage 2) demonstrates the diffusive profiles within the solid solution regimes of the particle between the phase interfaces. A movie of the simulated process is included in the supplement.
IV.2.2 Spinodal Decomposition
Spinodal decomposition is the process of transitioning from a high energy, linearly unstable, homogeneous system to a phase separated system. It is characterized by a gradual process of domain coarsening in which initially small phase regions grow and coalesce in order to minimize energetically expensive interfaces between phases Balluffi et al. (1954). By setting the rate constant to zero (, which imposes zero current into each layer) and beginning with a randomly perturbed filling fraction of 0.5 in each layer, we can examine the spinodal decomposition process for this coupled, two-layer model. At early times, both layers undergo internal spinodal decomposition, much as predicted by typical Cahn-Hilliard models Balluffi et al. (1954). However, because the layers are coupled via the overall free energy relation in Eq. 6 such that it is energetically unfavorable for both layers to be lithium-rich at the same position and time, the layers coordinate to form a checkerboard pattern as shown in Figure 11. Then, domain coarsening occurs during which the internal domains expand in size as the most central domain shrinks and disappears to minimize the circumferential interface between the domains. Throughout the process, interfaces move in concert to maintain near-equilibrium concentrations within each domain while preserving the average filling fraction within each layer. The final structure in Figure 11 (f) with a single interface in each layer is stable because the layers are unable to exchange lithium, forcing the average concentration within each layer to remain constant. Of note, the depicted rings of stage 1 (green) and stage 1’ (yellow) in Figure 11 (b) - (f) are related to the large penalty for intermediate concentrations in both layers from the term in Eq. 6. Because of the strong penalty, the domains “shift” relative to each other on the adjacent planes such that they either overlap slightly or are separated in space to avoid both having an interface (and intermediate filling fractions) at the same location (Figure 11 (b) - (f)). When is set to zero in an otherwise identical simulation, the rings do not appear.
IV.2.3 Constant Current Intercalation
When (dis)charging a particle at a fixed current, we can observe multiple regimes. At low enough currents (defined by comparing system time scales), the observed voltage is set by the effective chemical potential of the solid, as discussed in Section IV.1.1. However, at larger currents, other losses in the system such as internal gradients and reaction losses contribute to the measured voltage. In particular, we find that when the imposed time scale from the specified current density approaches that of the reaction time scale, we enter a regime in which the behavior of the observed voltage is dominated by the reaction losses. As an example, we present the results of a discharge for a particle with radius and . The value of was semi-arbitrarily chosen such that the reaction time scale was of similar magnitude to the imposed time scale of 1 hr. The reaction time scale for this system is . In Figure 12, we consider two cases to compare the effect of the chosen model for the reaction transition state. In Figure 12 (a), we set the value, , whereas in Figure 12 (b), we retain the model originally proposed in Eq. 18. In the first, we obtain a voltage curve similar to that in Figure 7 but shifted downward from reaction resistance and showing features of transients from transport losses between the surface and bulk of the particle. The mild increase in voltage over the second plateau is a result of a growing annulus of stage 1 graphite, causing the surface concentration to gradually approach a stable equilibrium value from a metastable value (see movies of concentration profiles in supplement). However, in Figure 12 (b), we see the system has a voltage curve that is nearly uncorrelated with the effective chemical potential along the equilibrium path. As mentioned above, deviations from equilibrium voltage curves can come from both transport and reaction losses. Reaction losses arise from slow reaction kinetics leading to large overpotentials, , which indicate the departure of the interfacial voltage from the equilibrium value Bazant (2013). Thus, we examine the reaction rate prefactor, called the exchange current density, which corresponds to a reaction conductance. We find the voltage is directly related to an effective reaction resistance, defined by taking two reaction resistances in parallel for the two layers. Noting that each reaction resistance scales as (from Eq. 16), we define the quantity
[TABLE]
which differs functionally from the actual reaction resistance only because the reaction driving force, , can be different for each layer. We compare this quantity with the system voltage in Figure 12 (b) and observe that the two quantities are correlated: large reaction resistance, , corresponds to lower system voltage, indicating that this quantity (not changes in at the surface) is the primary contribution to departures of the observed voltage from the equilibrium value at the overall filling fraction. Surprisingly, the voltage has a sharp increase near half-filling rather than a decrease as in the equilibrium profile in Figure 7, something that could not be predicted by a simple diffusion-based model with imposed phase boundaries. Here, the voltage increases near half filling because the current transitions from being split by the two layers, both of which are enlarging internal domains, to being carried largely by layer 2. The associated steeper internal concentration gradient near the surface of layer 2 leads to values of closer to 0.5 and smaller overall reaction resistance, despite the fact that the current is carried by only half the surface area (see concentration profiles in Figure 13). Throughout, the concentration at the surface is the result of coupled transport and reaction processes, and in particular the steep voltage drop off near complete filling is related to transport losses causing the surface filling fraction to approach unity which also leads to diverging reaction resistance. The large spikes in the voltage are related to surface concentration transitions from low to high values, leading to temporary low resistance at the intermediate filling fractions. It is worth mentioning that some of these features would be smoothed and not visible macroscopically in an ensemble of coupled particles as in a porous electrode.
Thus, we see that the choice of the reaction model can have significant impact on the predicted macroscopic properties. These observations underscore the importance of careful consideration of the reaction model when attempting to relate model predictions with experimental data, as this model makes strong connections between the microstructure and the single-particle discharge voltage profile. The single-particle discharge curve predicted in Figure 12 (b) bears little resemblance to experimental discharge curves of graphite electrodes Safari and Delacourt (2011), which suggests that, at least for the secondary graphite particles used in electrodes Wissler (2006), Eq. 18 may not correctly describe the particle-electrolyte interface, and a reaction model more like that used in Figure 12 (a) or Eq. 28 may be more appropriate, as we find in ref. Thomas-Alyea et al. (2016). We also note that in a porous electrode with many particles, some of the details of single-particle voltage curves can be masked by particle-particle interactions and system-scale transport losses Dreyer et al. (2010); Ferguson and Bazant (2014).
Because the term in the free energy involving was originally chosen for a very different model and only more strongly penalizes states of intermediate filling in both layers (which is already unfavorable even with ), we chose this as a representative simulation to study the model without that term. In this case, the simulated concentration profiles look qualitatively similar (see video in supplement), but the simulation neglecting the term ran considerably faster, so it may be reasonable to neglect this term in full spatially resolved models of graphite based on this two-layer framework.
IV.2.4 Inter-layer (Homogeneous) Reaction
As discussed above, various physical phenomena could enable exchange between the two concentration fields. In the case of the graphite model, a spatially regular presence of defects could allow internal exchange between layers of intercalated lithium Yao et al. (2012) to be modeled as a homogeneous reaction term. In order to properly formulate the reaction term in the context of non-equilibrium thermodynamics, simple mass action kinetics cannot apply, and the diffusional chemical potential must be the driving force for the reaction. Following Bazant Bazant (2013) (Eq. 7), we note that the rate of inter-layer exchange from the first to second layer can be written as
[TABLE]
where is the activity coefficient for the inter-layer reaction transition state, and we have made use of the fact that is identical for each layer to absorb reference chemical potential factors into the inter-layer rate constant, . For asymmetric lattices, the above should be modified to include the reference chemical potentials for each layer, biasing the reaction relative to equal activity in both lattices. Again, we are left a choice in the activity coefficient of the transition state. Here, we consider a reaction in which a species hops from a site on lattice 1 to a corresponding site on lattice 2. A natural choice for the transition state, then, is to assume that it excludes a single site in both layers, such that it is given simply by
[TABLE]
Other choices could include terms weighting enthalpic contributions from initial and final states (similar to the model proposed in Eq. 31); however, the above seems to be a simple reasonable model of the reaction, such that
[TABLE]
Thus, we have a slightly modified version of the conservation equation in each layer, such that instead of having zero homogeneous reaction term in Eq. 13, we have expressions for the source term, . In Figure 14, we compare the concentration profiles predicted during constant current intercalation with for the same set of conditions as presented in Figure 13 in which was set to zero. This corresponds to a Damköhler number of
[TABLE]
indicating relatively slow inter-layer reactions compared to the bulk transport. Nevertheless, even with this small exchange rate, we see qualitative differences in the model predictions.
With the inter-layer reaction, because the layers can exchange species, the system is able to fill a single layer completely before the other layer fills, thus minimizing the number of interfaces to follow a lower energy path during the filling. Physically, this occurs because as the lithium enters at the surface of either layer, it raises the diffusional chemical potential (and activity) in that layer. However, as soon as layer 2 phase separates, it is able to support increasing amounts of lithium without a significant change in its diffusional chemical potential by moving the phase boundary. In layer 1, as the concentration increases, it can go beyond the binodal point into the metastable region and reach a higher diffusional chemical potential than that in layer 2 (see Figure 3 (c)), thus causing species from layer 1 to transfer to the lower activity layer 2, which simply advances the position of the interface. Once layer 2 is completely full, layer 1 fills independently, moving a phase boundary from the edge of the particle toward the center, much like the profiles predicted in similar systems with only one concentration variable Zeng and Bazant (2014).
IV.3 Comparison to Solid Solution Particles
Although graphite phase separates upon intercalation with lithium, it is commonly modeled in battery electrodes using various simplifications to avoid the challenges of properly capturing phase separation dynamics. For example, it is common to approximate the transport of lithium in the material as a solid solution Yu et al. (1999); Verbrugge and Koch (2003); Srinivasan and Newman (2004a); Safari and Delacourt (2011); Bernardi and Go (2011); Baker and Verbrugge (2012), which has the advantage of being computationally very straightforward. By employing a concentration-dependent diffusivity, these models can even predict steep concentration gradients similar to phase boundaries Baker and Verbrugge (2012), but these gradients relax to uniform concentrations under zero current, even within spinodal regions. It is also possible to describe the intercalation as multiple processes involving solid-solution diffusion and phase boundary motion which enables use of electroanalytical techniques with simple expressions to determine values of transport coefficients Bard and Faulkner (2001); Aurbach et al. (1998); Levi and Aurbach (1997a), particularly in the solid solution regimes Funabiki et al. (1998), but doing so can lead to results in which the chemical diffusivity remains positive even in spinodal regions Levi et al. (1997); Levi and Aurbach (1997b) where uphill diffusion should occur. This can be interpreted as an “effective” transport parameter.
Another approach involves explicitly solving for motion of a phase boundary via a Stefan condition assuming phase equilibrium at the interface and solid-solution diffusion elsewhere Funabiki et al. (1999); Baker and Verbrugge (2012); Heß and Novák (2013). However, when solving numerically, the creation, elimination, and tracking of arbitrary numbers of phase boundaries becomes cumbersome. Moving boundary models also neglect physical contributions of interfacial energy and stresses (not modeled here but naturally a part of this framework Bazant (2013)), both of which influence the location of phase interfaces Cogswell and Bazant (2012). They also would not capture the uphill diffusion process involved in spinodal decomposition (Section IV.2.2), metastable concentrations reached when minimal gradients are present (e.g. the overshoot in Figure 7), or the effects of surface (de)wetting properties Zeng and Bazant (2014). Using a free energy approach like that developed here, open circuit voltages are emergent properties of the free energy, and phase interfaces develop naturally and do not need to be created artificially and tracked numerically. Nevertheless, for near-equilibrium systems without significant surface (de)wetting, the moving boundary and free energy approaches can lead to similar results.
Levi et al. compared the solid solution and moving boundary approaches in ref. Levi et al. (2005) and found that for many cases, the two lead to similar fit transport parameters. This similarity in predicted electrochemical outputs could help explain the success of solid solution models of graphite when applied to porous battery electrodes and their ability to fit macroscopic electrochemical data Verbrugge and Koch (2003); Safari and Delacourt (2011); Bernardi and Go (2011); Baker and Verbrugge (2012); Thomas-Alyea et al. (2016). In each case, transport parameters fit using the above models could be interpreted as describing “effective” transport within the active material particles, perhaps a combination of transport within primary particles and along grain boundaries of highly polycrystalline graphite secondary particle agglomerates Wissler (2006).
Although a comprehensive comparison of solid solution and phase separating models is beyond the scope of this work, we develop and present a simple example which highlights the differences in the predictions of a simple solid solution model of graphite and the phase separating model developed here – a current pulse-relaxation process in which the average concentration remains within the miscibility gap. We will not present the details of solid solution models, as that has been extensively studied elsewhere Newman and Thomas-Alyea (2004); Doyle et al. (1993); Fuller et al. (1994). Because the structures of the models are different, some care must be taken to obtain comparisons with few enough differences to reasonably interpret the outputs. First, we consider the reaction model. In both models, the driving force is the thermodynamic driving force, , which depends on the diffusional chemical potential at the surface of the solid. In the solid solution models, this depends only on the concentration at the surface, whereas the free-energy based models permit dependence on additional factors such as gradients of concentration. To match exchange current densities in the simplest way, we replace Eq. 16 with a form using a constant exchange current density with arbitrarily set ,
[TABLE]
with for each layer in the phase separating model and unity in the solid solution model. This could be interpreted as a particular functional form assumed for . The solid solution diffusional chemical potential is defined entirely by the open circuit voltage, and we use a function which matches the blue dashed curve in Figure 16 (b) to match the stable-equilibrium open circuit voltage of the phase separating model.
For the transport processes, we want to compare situations with similar diffusive behavior, at least outside of the miscibility gap. For the solid solution model, we use a constant and uniform solid solution chemical diffusivity, such that the solid solution flux is given simply by . As presented in Appendix C, rather than matching the flux prefactors, and , we approximately match the concentration dependence (but not the magnitude) of the diffusive behavior in the solid solution regimes of the particles (far from the interfaces in the phase separating particle). To do this we retain the flux expression in Eq. 12 with and use . The value of is insignificant to the interpretation of the results, as described below. We simulate both particles as cylinders of radius with a uniform mesh of 400 volumes.
In Figure 15, we compare simulations of the solid solution and phase separating particles exposed to a 2 C current pulse for 300 s beginning after a low-current filling to reach an overall filling fraction of 0.6. We immediately see differences between the initial states, the lightest orange lines in Figure 15 (a) and (b). The solid solution model begins this pulse from a uniform 0.6 filling fraction, whereas the phase separating model begins with a stable phase interface and a high surface concentration. The concentration profiles for the solid solution and phase separating particles (Figure 15 (a) and (b) respectively) are qualitatively different. The solid solution model leads to gradually sloping concentration profiles, which relax toward a uniform concentration again after the current pulse ends. The phase separating model fills by propagating a phase interface in toward the center of the particle and retains a stable interface after relaxation.
The simulated voltage profiles in Figure 15 (c) show that the phase separating particle has transport losses during the current pulse because the surface diffusional chemical potential can rise. In contrast, the solid solution particle voltage is fixed only by the reaction losses because the equilibrium voltage is only a function of concentration and is within the voltage plateau over the range of filling fractions studied. As a result, the predicted voltage profile for the solid solution model with this current pulse is relatively insensitive to the solid solution diffusivity (as long as the surface concentration remains within the miscibility gap). That is, for this current pulse, any value of larger than the simulated value would lead to more uniform concentration profiles and give an identical voltage profile. Of note, when using the phase separating flux prefactor defined by Eq. 37, as we use in ref. Thomas-Alyea et al. (2016), in an otherwise similar simulation the voltage profiles are more similar because that form leads to larger flux prefactors in the stable phases (and smaller gradients), but the phase separating voltage profile still retains some decreasing slope. In multi-particle porous electrode simulations with particle-particle interactions Chung et al. (2014); Stephenson et al. (2007), both models could predict similar transport losses because of non-uniform particle filling and other losses Newman and Thomas-Alyea (2004); Garcıá et al. (2005), but in the case of single particles, the models make qualitatively distinct predictions, both of the concentration profiles in the particles and their associated voltage profiles.
Here, we should reemphasize that solid solution models of graphite have had good success capturing macroscopic porous electrode behavior. The present model is designed to accurately capture the single crystal dynamics, which may not strongly affect the macroscopic current and voltage measurements of full porous electrodes with complex secondary electrode particles. Nevertheless, consistently coupling a phase separating single-crystal model to a hierarchy of kinetic processes including grain boundary diffusion should more consistently predict overall electrode lithium distribution. This would lead to more accurate surface concentration predictions, which are important to describe accurately when developing models to predict lithium plating risk Harris et al. (2010) as well as a more consistent coupling to stresses Di Leo et al. (2014); Schiffer et al. (2016), which we have neglected here but are often studied using solid solution models Christensen and Newman (2006); Cheng and Verbrugge (2010); Grantab and Shenoy (2011); Woodford et al. (2010).
V Reduced Model for Graphite Battery Electrodes
By enforcing two-layer periodicity, our free-energy model is able to capture the thermodynamics of lithium intercalation at high filling fractions, but it cannot describe the plethora of stable or metastable phases in graphite at low filling fractions (not only for lithium), which exhibit longer-range periodicity across three or more layers. The complex phase behavior is reflected in the irregular rise of the “voltage staircase” at low filling fractions, which cannot be easily described by the two-variable model (Figure 7). Moreover, even with these physical limitations, the model is much more computationally expensive than phase-field intercalation models with a single concentration variable, such as that studied by Zeng and Bazant Zeng and Bazant (2014). More realistic models with three or more concentration variables to capture low-density phase behavior or including other effects such as elastic coherency strain Cogswell and Bazant (2012) would be even more costly. Similar concerns would apply to the modeling of other multilayer multiphase materials mentioned in the Introduction. There is a trade-off between explicitly describing the microscopic phase behavior with multiple concentration variables (for periodically repeating layers or interpenetrating lattices) and accurately fitting the chemical potential or voltage with a computationally efficient model.
To illustrate the construction of such a reduced model, we develop a single-variable free-energy model for lithium intercalation in graphite, which fits the open circuit voltage at low filling fractions as an effective solid solution, while still capturing the two primary voltage plateaus at high filling fractions, albeit without allowing for symmetry breaking (e.g. checkerboard patterns) or any prediction of temperature dependence of the staging behavior. Because of the large energetic barrier to fill along the fully homogeneous yellow paths in Figure 3, particles typically follow along the lower energy blue path. To simplify the model, we began with a function which follows the red free energy path to capture the energetics of the unstable phases as a function of only an average concentration. We further made the assumption that below a filling fraction of , the high stage numbers can be effectively treated as a solid solution, which allowed us to modify the homogeneous free energy function (or equivalently, the homogeneous diffusional chemical potential) to match the voltage in this low-filling-fraction region. Doing so ignores any higher order staging and stage transitions, but leads to a simple model which is better suited to scale up within multi-particle simulations, e.g. for porous, phase separating battery simulation as we explore in a companion work Thomas-Alyea et al. (2016). We also adjusted the high-filling behavior to have a more gradual drop off, as observed in various experiments Verbrugge and Koch (2003); Thomas and Newman (2003); Dahn et al. (1990). Finally, we adjusted the decreasing regions of diffusional chemical potential (increasing equilibrium voltage) to reduce the metastability region which causes open circuit voltage hysteresis Levi and Aurbach (1997a); Dreyer et al. (2010). The diffusional chemical potential function we used which satisfies the criteria above is presented in Appendix D and in Figure 16 noting that . The free energy permits construction of two common tangents between filling fractions near and and another between filling fractions near and which leads to the two clear voltage plateaus. We emphasize that, although this approach predicts phase separation based on an average filling fraction, it cannot capture checkerboard microstructures or make predictions about temperature dependence of the phase diagram (Section IV.1.3), which are possible with the two-variable model above.
Simulations using this free energy lead to predicted concentration profiles with stable average filling fractions near and in the high-filling states as in the two-parameter model, but the low-filling state has an average filling fraction near , which is clearly demonstrated in a constant current discharge like those presented in Figure 17 (compare to Figures 13 and 14). The predicted average concentration profiles look similar to those with inter-layer reaction but miss the formation of internal rings of stage 1 () predicted by the original 2-parameter model.
We apply this model within a porous electrode simulation in ref. Thomas-Alyea et al. (2016) and demonstrate good agreement with both macroscopic voltage transients and visual indicators of lithium distribution within a graphite electrode responding to current pulses. We also find that using Eq. 12 and the reaction model described in Eqs. 16 and 18 does not provide reasonable fits to the experimental data, but we are able to fit data using the flux expression defined in Eqs. 10, 11, and 37. The reaction model we use to fit the data is a Butler-Volmer expression dominated by film resistance. The general framework developed here is inspired by simple thermodynamic arguments developed for single crystals, whereas practical graphite electrode particles are secondary agglomerates composed of many smaller crystals joined by grain boundaries Wissler (2006). This suggests that these modifications to the general model are capturing effective properties of the secondary particles and provide a better starting point to describe graphite electrodes.
VI Summary and Conclusions
We have developed and studied a model for the kinetics of intercalation for layered, phase-separating compounds that tend to form staged structures. Using the particular example of lithium in graphite, the most common lithium-ion battery anode, we specialized the model to describe the most visually distinct low-stage number structures of lithiated graphite as well as the staircase character of its equilibrium voltage. Having applied the model to a single particle in ref. Guo et al. (2016), we explored the generalization of the model as well as some of the unique predictions it makes which differ from other models commonly employed to describe lithiated graphite. We found the specifics of the choice for reaction model to be critically important in affecting practical macroscopic simulation outputs such as predicted voltage under constant current intercalation (discharge). In particular, the reaction model used for the single-crystal simulations and experiment studied in ref. Guo et al. (2016) does not lead to qualitative agreement with typical graphite electrode discharge curves, which supports the idea that for secondary graphite electrode particles as we study in ref. Thomas-Alyea et al. (2016), the reaction model described in Section V is more appropriate. Simulations directly comparing the model to more commonly employed solid solution models highlight some of the key differences between the models, including concentration profiles of particles under current pulses and the associated voltage profiles. In addition, we developed a simplified model designed to accurately match observed graphite open circuit voltage profiles while retaining the overall phase separating characteristics of the 2-layer model at a substantially reduced computational cost. This facilitated the work in ref. Thomas-Alyea et al. (2016), in which we show that the model can reproduce electrochemical data as well as visual indications of concentration carried out in situ during current pulses. Natural extensions of the work include more carefully capturing the temperature dependence of the phase diagram by refining the free energy model. We have also neglected here the effects of elastic contributions to the energetics and dynamics of the materials, and incorporating them is a natural next step to make the model more physically descriptive of graphite Dahn et al. (1982); Qi et al. (2010); Qi and Harris (2010); Bower et al. (2011); Sethuraman et al. (2012); Qi et al. (2014) in a way that consistently couples the stresses to the phase separating concentration profiles.
While much of the focus here was on lithiated graphite, the general form of the model could be useful for other materials exhibiting staging or similar behavior, such as some layered double hydroxides Pisson et al. (2003), or metal borocarbides Joshi et al. (2015). MXenes Wang et al. (2015); Naguib et al. (2014) have also shown behavior analogous to staging, indicating the general framework developed here may help capture their interactions with intercalants. Because many intercalants in graphite exhibit similar behavior Dresselhaus and Dresselhaus (1981), the model could likely be adapted to describe other graphitic electrodes such as those intercalated with aluminum Lin et al. (2015) or sodium Kim et al. (2015). It may also be extensible to other sodium-based layered systems with non-trivial phase behavior Wu et al. (2015). Beyond the layered framework developed here, the model’s general structure of spatially overlapping lattices may be useful for similar phenomena such as that of different occupational sites available for Li in the electrode material anatase TiO2 Zachau-Christiansen et al. (1988); Lafont et al. (2010); Shin et al. (2011); Shen et al. (2014); de Klerk et al. (2016) or different interstitial sites available for hydrogen in metals Di Leo and Anand (2013).
The choice to develop this model with a thermodynamically consistent connection to electrochemical reaction kinetics makes it particularly suited to describe systems with internal phase separation as well as staging and the associated implications on battery and supercapacitor performance, but the generality of the model enables adaptation to describe intercalation kinetics of other related behaviors. Because the model is based on physical understanding of system microstructure, it is well suited to parameterization by ab initio calculations and can be used to infer microscopic physical mechanisms for transport and intercalation from experimental data. This allows it to provide meso-scale understanding of system behavior while retaining a clear connection to both microscopic modeling and macroscopic experiments with the generality to be specialized and applied to a wide range of systems.
VII Acknowledgements
R. B. Smith and M. Z. Bazant acknowledge support from the Samsung Advanced Institute of Technology. E. Khoo acknowledges support from the National Science Scholarship (PhD) funded by Agency for Science, Technology and Research, Singapore (A*STAR). The authors gratefully thank K. E. Thomas-Alyea for helpful discussions throughout, guidance toward the one-parameter model reduction, and very thoughtful criticisms of the manuscript. We also thank E. Rejovitzky for insightful criticisms of the manuscript.
Appendix A Transport coefficient for diffusion on a lattice
The activity coefficients are related to excess chemical potentials, , so using Eq. 9,
[TABLE]
Then, we postulate a diffusion transition state activity coefficient in which the diffusing species excludes two sites and maintains the bulk enthalpic contributions to its diffusional chemical potential Bazant (2013),
[TABLE]
which, combined with Eq. 11 leads to Eq. 12.
Appendix B Simulation details for examining axial symmetry assumption
In the interest of keeping the 2D simulation time to a few hours, we perform the simulation up to instead of as in ref. Guo et al. (2016). We impose no-flux boundary conditions on the straight edges of the sector and lithium intercalates along the arc of the sector. In Figure 18, we plot at where axial symmetry is broken, and at where axial symmetry is not broken. We provide plots of against at and to examine how axial symmetry is broken or not broken. In addition, we compare these plots with those obtained from a 1D simulation in COMSOL Multiphysics using the same grid spacing of to examine how well the 1D model approximates the 2D model. This grid spacing is larger than the approximate value of the interface width in these simulations (), which is a compromise on accuracy to speed up the simulations and enable a larger central angle for the 2D sector. We retain the same grid spacing in both cases in Figure 18 to focus on the comparison rather than accurate model predictions. Full movies are provided in the supplement.
In Figure 18, comparing plots (a), (c) and (e) for , we observe that axial symmetry in the profile is broken; for instance, at the arc of the sector, increases with increasing . On the other hand, comparing plots (b), (d) and (f) for , does not vary as a function of and axial symmetry is not broken. More generally, even though we initialize and to be uniformly, we observe that axial symmetry is broken immediately, turns “less broken” as increases, first becomes unbroken at (not shown in Figure 18, but can be seen in the full videos in the supplement), and then remains unbroken till . Although we did not perform 2D simulations up to as was done in ref. Guo et al. (2016), we expect that axial symmetry remains unbroken from to . Comparing plots (c) and (e) with (g), the broken axial symmetry at causes some deviations between the profiles for , and the true 1D simulation. However, comparing plots (d) and (f) with (h), these profiles agree reasonably well with each other when axial symmetry is not broken at ; this also holds true when axial symmetry first becomes unbroken at . Crucially, the four internal domains for in the 2D simulation are also reproduced by the 1D simulation. Therefore, the 1D simulation approximates the 2D simulation reasonably well except at early times where there are minor deviations that do not significantly affect the accuracy of the 1D simulation at later times.
Appendix C Matching diffusive behavior in solid solution and phase separating simulations
To approximately match the transport behavior of the two models, it is helpful to rearrange the flux expression in Eqs. 10 and 11 in terms of a chemical diffusivity, which is valid in regions where the gradient term is negligible (far from phase interfaces),
[TABLE]
where
[TABLE]
Here, we are constrained on our choice for the ratio by the second law of thermodynamics Groot and Mazur (1962), so we cannot choose it such that we have perfectly constant . We note that the choice originally proposed in Appendix A and simplified in Eq. 12 gives us Guo et al. (2016)
[TABLE]
Another natural choice may be to assume constant prefactor in front of the , or
[TABLE]
leading to
[TABLE]
which has been used in a similar model Hawrylak and Subbaswamy (1984) and has diverging chemical diffusivity near full and empty states. We also find this model specified by Eq. 37 to match porous electrode data in our companion paper in which the particle models are describing effective properties of secondary (polycrystalline) graphite particles Thomas-Alyea et al. (2016). However, to retain the most similar behavior of chemical diffusivity between the solid solution and phase separating models, we focus here on the model proposed in Eq. 12.
Appendix D Diffusional chemical potential for single-variable model reduction
In the reduced thermodynamic model for lithium intercalation in graphite, the function describing the diffusional chemical potential as a function of a single total filling fraction variable is given by
[TABLE]
where
[TABLE]
and step up and step down functions respectively defined by
[TABLE]
The homogeneous free energy can be computed from an integral of the homogeneous contribution to the diffusional chemical potential.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Wells et al. (1996) B. O. Wells, R. J. Birgeneau, F. C. Chou, Y. Endoh, D. C. Johnston, M. A. Kastner, Y. S. Lee, G. Shirane, J. M. Tranquada, and K. Yamada, Zeitschrift für Physik B Condensed Matter 100 , 535 (1996) . · doi ↗
- 2Mohottala et al. (2006) H. E. Mohottala, B. O. Wells, J. I. Budnick, W. A. Hines, C. Niedermayer, L. Udby, C. Bernhard, A. R. Moodenbaugh, and F.-C. Chou, Nature Materials 5 , 377 (2006) . · doi ↗
- 3Ijdo and Pinnavaia (1998) W. L. Ijdo and T. J. Pinnavaia, Journal of Solid State Chemistry 139 , 281 (1998).
- 4Iyi et al. (2002) N. Iyi, K. Kurashima, and T. Fujita, Chemistry of Materials 14 , 583 (2002) . · doi ↗
- 5Pisson et al. (2003) J. Pisson, C. Taviot-Gueho, Y. Israëli, F. Leroux, P. Munsch, J.-P. Itie, V. Briois, N. Morel-Desrosiers, and J.-P. Besse, The Journal of Physical Chemistry B 107 , 9243 (2003).
- 6Naguib et al. (2012) M. Naguib, J. Come, B. Dyatkin, V. Presser, P.-L. Taberna, P. Simon, M. W. Barsoum, and Y. Gogotsi, Electrochemistry Communications 16 , 61 (2012) . · doi ↗
- 7Wang et al. (2015) X. Wang, S. Kajiyama, H. Iinuma, E. Hosono, S. Oro, I. Moriguchi, M. Okubo, and A. Yamada, Nature Communications 6 , 6544 (2015) . · doi ↗
- 8Joshi et al. (2015) R. P. Joshi, B. Ozdemir, V. Barone, and J. E. Peralta, The Journal of Physical Chemistry Letters 6 , 2728 (2015) . · doi ↗