Efficient formulation and simulation of solid-phase diffusion in silicon-like volume-expanding materials

Silicon-like electrode materials expand several-fold on lithiation, so solid-phase diffusion must be solved in a sphere whose radius, volume and host-site concentration change with the state of lithiation. The reduced-order solutions developed for fixed particles (polynomial profiles, eigenfunction Galerkin collocation and, most recently, Galerkin weak forms with symmetric polynomials for constant and concentration-dependent diffusivity) cannot be applied to such a particle directly, because the governing equation contains a lattice velocity, a dilution term and a moving boundary. We show that, for a homogeneous isotropic expansion of the particle, the lithium balance of a swelling particle reduces exactly to the standard problem of spherical diffusion in a fixed sphere when it is written for the scaled moles X = R_a³ C_a / (R_a0³ C^H_a0) (the moles per unit initial volume, that is, the referential concentration of finite-deformation diffusion theory, which is also the local occupancy of the host sites) in a coordinate that moves with the lattice. In this variable the time derivative is multiplied by τZ², where Z = R_a/R_a0 is the stretch of the particle, the surface flux is Z²τQ/3 and the average of X is the state of lithiation exactly. For a prescribed stretch history this holds for any diffusivity law; with the closure Z³ = f(X_avg) the radius and its rate of change, which appear explicitly in the concentration equation, are eliminated through the average of X, and for a concentration-dependent diffusivity the argument of the diffusivity carries the time-dependent stretch. Because the swelling problem takes the standard form, the weak-form equations and tables of the symmetric-polynomial formulation apply directly. We give them for constant diffusivity, for the moment form and for concentration-dependent diffusivity, in non-dimensional and in dimensional form, together with the relations needed to couple them to a pseudo-two-dimensional electrode model. For a volume law that reaches 3.75 at full capacity, charged to 2.47 times the initial volume, the weak form with three to six unknowns reaches errors of 10⁻⁴ to 10⁻⁸ in the surface occupancy after the first 100 s for constant diffusivity and 10⁻⁴ to 10⁻⁷ for a diffusivity D_ref·e^(6X), where a cell-centered finite-difference scheme needs 32 to more than 128 cells. Over the whole history, including the first seconds and the instants after a change of current, the low-order weak forms are less accurate than a fine mesh. The moles are conserved to round-off by all schemes.