A Cahn–Hilliard-type phase-field theory for species diffusion coupled with large elastic deformations: Application to phase-separating Li-ion electrode materials

Abstract

We formulate a unified framework of balance laws and thermodynamically-consistent constitutive equations which couple Cahn–Hilliard-type species diffusion with large elastic deformations of a body. The traditional Cahn–Hilliard theory, which is based on the species concentration c and its spatial gradient $\nabla c$, leads to a partial differential equation for the concentration which involves fourth-order spatial derivatives in c; this necessitates use of basis functions in finite-element solution procedures that are piecewise smooth and globally ${\mathcal{C}}^1\mathrm{-}\mathrm{continuous}$. In order to use standard ${\mathcal{C}}^0\mathrm{-}\mathrm{continuous}$ finite-elements to implement our phase-field model, we use a split-method to reduce the fourth-order equation into two second-order partial differential equations (pdes). These two pdes, when taken together with the pde representing the balance of forces, represent the three governing pdes for chemo-mechanically coupled problems. These are amenable to finite-element solution methods which employ standard ${\mathcal{C}}^0\mathrm{-}\mathrm{continuous}$ finite-element basis functions.

We have numerically implemented our theory by writing a user-element subroutine for the widely used finite-element program Abaqus/Standard. We use this numerically implemented theory to first study the diffusion-only problem of spinodal decomposition in the absence of any mechanical deformation. Next, we use our fully coupled theory and numerical-implementation to study the combined effects of diffusion and stress on the lithiation of a representative spheroidal-shaped particle of a phase-separating electrode material.

Introduction

Several Li intercalation compounds of interest in battery applications exhibit phase-separation (cf., e.g., Bruce et al., 2008, Tang et al., 2010) and central to the study of coupled diffusion–deformation problems in such electrode materials is the Cahn–Hilliard phase-field theory (Cahn and Hilliard, 1958, Cahn and Hilliard, 1959, Cahn, 1961). As background, first consider the classical Cahn–Hilliard theory for species diffusion and phase segregation which is uncoupled from the mechanical problem. Let $\overline{c}\in [0,1]$ denote a normalized species concentration, ${\widehat{\psi }}_{\mathrm{R}}(\overline{c},\nabla \overline{c})$ the free energy per unit reference volume, and let the functional

$$\Psi \left(\overline{c}\right)={\int }_{\mathrm{B}}{\widehat{\psi }}_{\mathrm{R}}(\overline{c},\nabla \overline{c}){\mathit{dv}}_{\mathrm{R}}$$

denote the total free energy of the region of space occupied by the body $\mathrm{B}$. Then classically, for theories in which the free energy depends on $\nabla \overline{c}$, the chemical potential μ is defined as the variational derivative of the functional $\Psi \left(\overline{c}\right)$:

$$\mu \stackrel{\mathrm{def}}{=}{\displaystyle \frac{\delta \Psi }{\delta \overline{c}}}=\frac{\partial {\widehat{\psi }}_{\mathrm{R}}(\overline{c},\nabla \overline{c})}{\partial \overline{c}}-\mathrm{Div}\left(\frac{\partial {\widehat{\psi }}_{\mathrm{R}}(\overline{c},\nabla \overline{c})}{\partial \nabla \overline{c}}\right).$$

A widely used specific form for the free energy is the following separable energy first proposed by Cahn and Hilliard:

$${\widehat{\psi }}_{\mathrm{R}}(\overline{c},\nabla \overline{c})={\psi }_{\mathrm{R}}^c\left(\overline{c}\right)+{\displaystyle \frac{1}{2}}{\lambda }_{\mathrm{CH}}|\nabla \overline{c}{|}^2,{\lambda }_{\mathrm{CH}}>0.$$

Here, ${\psi }_{\mathrm{R}}^c\left(\overline{c}\right)$ represents the coarse-grain chemical free energy, a double-well potential whose wells (the “binodal points”) define the phases; the second term which depends on the concentration gradient represents an interfacial energy, with ${\lambda }_{\mathrm{CH}}$ a gradient energy coefficient with units of energy per unit volume times a length squared. For the free energy (1.3), the chemical potential, using (1.2), is

$$\mu ={\displaystyle \frac{d{\psi }_{\mathrm{R}}^c\left(\overline{c}\right)}{d\overline{c}}}-{\lambda }_{\mathrm{CH}}\mathrm{\Delta }\overline{c},$$

where $\Delta$ is the Laplace operator. Next, species balance requires that

$$\dot{\overline{c}}=-\mathrm{Div}{\mathbf{j}}_{\mathrm{R}},$$

with the species flux ${\mathbf{j}}_{\mathrm{R}}$ related to the chemical potential μ through the constitutive equation:

$${\mathbf{j}}_{\mathrm{R}}=-m\left(\overline{c}\right)\nabla \mu ,$$

where $m\left(\overline{c}\right)$, the species mobility, is a nonlinear and positive function. Using (1.4), (1.6) in (1.5) yields the classical Cahn–Hilliard diffusion equation:

$$\dot{\overline{c}}=\mathrm{Div}\left(m\left(\overline{c}\right)\nabla \left({\displaystyle \frac{d{\psi }_{\mathrm{R}}^c\left(\overline{c}\right)}{d\overline{c}}}-{\lambda }_{\mathrm{CH}}\Delta \overline{c}\right)\right).$$

Two mechanisms dominate the evolution of a solution to the Cahn–Hilliard equation: a minimization of the chemical energy ${\psi }_{\mathrm{R}}^c\left(\overline{c}\right)$ drives the solution to binodal points and separates the phases, while minimization of the interface energy $(1/2){\lambda }_{\mathrm{CH}}|\nabla \overline{c}{|}^2$ effectively coarsens the phases.

The Cahn–Hilliard equation (1.7) involves fourth-order spatial derivatives of the concentration $\overline{c}$. For recent reviews and discussions of numerical solution techniques for the Cahn–Hilliard equation, cf., e.g., Wells et al. (2006), Kuhl and Schmidt (2007), and Gomez et al. (2008). As noted by Gomez et al. (2008), traditional numerical methods for dealing with higher-order operators on very simple geometries include finite differences and spectral approximations. However, for practical engineering applications, simple geometries are not very relevant, and more geometrically flexible techniques such as the finite-element method — which might allow for arbitrary two- and three-dimensional geometries — need to be utilized. In the context of the finite element method, the fourth order equation (1.7) typically requires basis functions that are piecewise smooth and globally ${\mathcal{C}}^1\mathrm{-}\mathrm{continuous}$, but there are only a very limited number of two-dimensional finite elements possessing ${\mathcal{C}}^1\mathrm{-}\mathrm{continuity}$ applicable to complex geometries, and none in three-dimensions. Instead of using ${\mathcal{C}}^1\mathrm{-}\mathrm{continuous}$ functions, the most common manner to solve this equation using finite elements, has been with split-methods (also known as mixed-methods) which in addition to the primal variable $\overline{c}$, introduce an extra degree of freedom (cf., e.g., Wodo and Ganapathysubramanian, 2011). We discuss next one such method proposed by Ubachs et al. (2004).¹

In contrast to the standard derivation of the Cahn–Hilliard equation discussed above, Ubachs et al. (2004) introduced another variable $\overline{\mathbb{c}}$, which they call a nonlocal species concentration, and took as the governing partial differential equation for their diffusion problem, the equation

$$\dot{\overline{c}}=\mathrm{Div}\left(m\right(\overline{c})\nabla \mu ),\mathrm{with}\mu ={\displaystyle \frac{d{\psi }_{\mathrm{R}}^c\left(\overline{c}\right)}{d\overline{c}}}+\beta (\overline{c}-\overline{\mathbb{c}}),\beta >0,$$

where the nonlocal field $\overline{\mathbb{c}}$ is evaluated by solving another Helmholtz-type partial differential equation:

$$\overline{\mathbb{c}}-{\ell }^2\mathrm{\Delta }\overline{\mathbb{c}}=\overline{c},\ell >0,$$

in which ℓ is a parameter with units of length. Note that using (1.9) we may eliminate the factor $(\overline{c}-\overline{\mathbb{c}})$ in (1.8)₂ to obtain

$$\mu ={\displaystyle \frac{d{\psi }_{\mathrm{R}}^c\left(\overline{c}\right)}{d\overline{c}}}-\lambda \mathrm{\Delta }\overline{\mathbb{c}},\mathrm{with}\lambda =\beta {\ell }^2>0.$$

This expression for the chemical potential is similar to the Cahn–Hilliard relation (1.4), except that in (1.10) one has the Laplacian of the nonlocal field $\mathrm{\Delta }\overline{\mathbb{c}}$ and two parameters $\beta$ and ℓ, instead of the Laplacian of the local field $\mathrm{\Delta }\overline{c}$ and a single parameter ${\lambda }_{\mathrm{CH}}$.

Eqs. (1.8), (1.9), being partial differential equations, require boundary conditions. The boundary conditions for (1.8) are standard. For (1.9) (Ubachs et al., 2004) introduced boundary conditions on either the value of $\overline{\mathbb{c}}$ or its normal derivative. In their calculations they used a homogeneous Neumann-type boundary condition of the form

$$(\nabla \overline{\mathbb{c}})·{\mathbf{n}}_{\mathrm{R}}=0\mathrm{on}\mathrm{the}\mathrm{external}\mathrm{boundary}\partial \mathrm{B}\mathrm{of}\mathrm{the}\mathrm{body},$$

where ${\mathbf{n}}_{\mathrm{R}}$ denotes the outward unit normal to the external boundary $\partial \mathrm{B}$ of the body $\mathrm{B}$. The resulting problem then consists of two coupled second-order partial differential equations (1.8), (1.9), which are solved using finite element methods which employ standard ${\mathcal{C}}^0\text{<mi mathvariant="italic" is="true">-</mi>}\mathrm{continuous}$ finite element basis functions.

Such a theory, which uses a nonstandard variable $\overline{\mathbb{c}}$ in order to introduce nonlocal effects, has met with substantial operational success when modeling microstructure evolution in two-phase alloys. (cf., e.g., Ubachs et al., 2004). However, from a physical point of view the “derivation” of the important partial differential equation (1.9) for the variable $\overline{\mathbb{c}}$ in Ubachs et al. (2004) is not entirely satisfactory.

• In particular, it is not clear whether (1.9) is a balance law, a constitutive equation, or a combination of the two. Also, it is not clear how to arrive at such an equation based on classical thermodynamic arguments.

On the other hand — following the formulation of micromorphic theories of continua² — if from the outset one introduces $\overline{\mathbb{c}}$ as an additional kinematical variable in the theory, then it is possible to develop a thermodynamically consistent theory in a systematic fashion using the principle of virtual power, in which the important relation (1.9) represents a microforce balance supplemented by thermodynamically consistent constitutive equations. We discuss such a formal development of the partial differential equation (1.9) in what follows.

As stated previously, one of the objectives of this paper is to formulate a theory which couples species diffusion with large elastic deformations of the body, with an accounting for the swelling and phase segregation caused by the diffusing species. In the development of the theory in the sections to follow — for clarity of the physical units of the various quantities introduced — we find it useful to formulate our theory in terms of the actual species concentration $c$, given in terms of moles per unit reference volume, rather than in terms of a normalized species concentration $\overline{c}\in [0,1]$ used thus far in this introductory Section. Also, we employ a variable $\mathbb{c}$ with units of moles per unit volume, and its gradient $\nabla \mathbb{c}$ — with the latter being chosen to represent a measure of the inhomogeneity of the microscale species concentration. For want of a better terminology, we call $\mathbb{c}$ the “micromorphic concentration.”

In addition to the standard balance (1.5) (in terms of standard un-normalized quantities), viz.

$$\dot{c}=-\mathrm{Div}{\mathbf{j}}_{\mathrm{R}},$$

for the species concentration $c$, we develop the other balances in our theory by following the virtual-power approach of Germain (1973) and Gurtin (2002), which results in a “macroforce” balance and two “microforce balances,” for the forces associated with the rate-like kinematical descriptors in our theory. Briefly, with $\mathit{\chi }$ denoting the motion, and $\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^c$ the decomposition of the deformation gradient into elastic and chemical parts, we define a generalized virtual velocity to be a list

$$\mathcal{V}=(\tilde{\mathit{\chi }},{\tilde{\mathbf{F}}}^e,\tilde{c},\tilde{\mathbb{c}}),$$

and for any part $\mathrm{P}$ of the body $\mathrm{B}$ we take the virtual internal and external powers in the theory to be given by

$${\mathcal{W}}_{\mathrm{int}}(\mathrm{P},\mathcal{V})={\int }_{\mathrm{P}}({\mathbf{S}}^e:{\tilde{\mathbf{F}}}^e+\pi \tilde{c}+\mathbb{p}\tilde{\mathbb{c}}+\mathit{\xi }·\nabla \tilde{\mathbb{c}}){\mathit{dv}}_{\mathrm{R}},{\mathcal{W}}_{\mathrm{ext}}(\mathrm{P},\mathcal{V})={\int }_{\partial \mathrm{P}}{\mathbf{t}}_{\mathrm{R}}\left({\mathbf{n}}_{\mathrm{R}}\right)·\tilde{\mathit{\chi }}{\mathit{da}}_{\mathrm{R}}+{\int }_{\mathrm{P}}{\mathbf{b}}_{\mathrm{R}}·\tilde{\mathit{\chi }}{\mathit{dv}}_{\mathrm{R}}+{\int }_{\partial \mathrm{P}}\zeta \left({\mathbf{n}}_{\mathrm{R}}\right)\tilde{\mathbb{c}}{\mathit{da}}_{\mathrm{R}}.$$

Here, the stress ${\mathbf{S}}^e$ is power conjugate to ${\tilde{\mathbf{F}}}^e$; the scalar microstress π is conjugate to $\tilde{c}$; the scalar microstress $\mathbb{p}$ is conjugate to $\tilde{\mathbb{c}}$; and the vector microstress $\mathit{\xi }$ is conjugate to $\nabla \tilde{\mathbb{c}}$. Also, for each unit outward normal ${\mathbf{n}}_{\mathrm{R}}$ to the boundary $\partial \mathrm{P}$ of $\mathrm{P}$, ${\mathbf{t}}_{\mathrm{R}}\left({\mathbf{n}}_{\mathrm{R}}\right)$ is the macroscopic traction which expends power on $\tilde{\mathit{\chi }}$; ${\mathbf{b}}_{\mathrm{R}}$ is body force that expends power on $\tilde{\mathit{\chi }}$; and $\zeta \left({\mathbf{n}}_{\mathrm{R}}\right)$ is a scalar microscopic traction which expends power on $\tilde{\mathbb{c}}$. The principle of virtual power then consists of two basic requirements:

• Given any part P,

  $${\mathcal{W}}_{\mathrm{ext}}(\mathrm{P},\mathcal{V})={\mathcal{W}}_{\mathrm{int}}(\mathrm{P},\mathcal{V})\mathrm{for}\mathrm{all}\mathrm{generalized}\mathrm{virtual}\mathrm{velocities}\mathcal{V}.$$

• Given any part $\mathrm{P}$ and a rigid virtual velocity $\mathcal{V}$,

  $${\mathcal{W}}_{\mathrm{int}}(\mathrm{P},\mathcal{V})=0\mathrm{whenever}\mathcal{V}\mathrm{is}\mathrm{a}\mathrm{rigid}\mathrm{macroscopic}\mathrm{virtual}\mathrm{velocity.}$$

Consequences of the virtual power-principle are

• The stress ${\mathbf{T}}_{\mathrm{R}}\stackrel{\mathrm{def}}{=}{\mathbf{S}}^e{{\mathbf{F}}^c}^{-\top }$ satisfies ${\mathbf{T}}_{\mathrm{R}}{\mathbf{F}}^{\top }={\mathbf{FT}}_{\mathrm{R}}^{\top }$ and represents the classical Piola stress which is consistent with the macroscopic force balance

  $$\mathrm{Div}{\mathbf{T}}_{\mathrm{R}}+{\mathbf{b}}_{\mathrm{R}}=\mathbf{0},$$

  and the macroscopic traction condition:

  $${\mathbf{t}}_{\mathrm{R}}\left({\mathbf{n}}_{\mathrm{R}}\right)={\mathbf{T}}_{\mathrm{R}}{\mathbf{n}}_{\mathrm{R}}.$$

• With ${\mathbf{M}}^e\stackrel{\mathrm{def}}{=}{J}^{c-1}{{\mathbf{F}}^e}^{\top }{\mathbf{S}}^e$ defining a Mandel stress, the microstress $\pi$ satisfies the microforce balance:

  $$\pi =J^c{\mathbf{M}}^e:\mathbb{A},$$

  where $J^c=\det {\mathbf{F}}^c$, and $\mathbb{A}$ is a tensor defining the “direction” of chemical expansion/contraction; cf. Eq. (A.10).

• The microscopic stresses $\mathbb{p}$ and $\mathit{\xi }$ are consistent with the microscopic force balance:

  $$\mathrm{Div}\mathit{\xi }-\mathbb{p}=0,$$

  and the microtraction condition:

  $$\zeta \left({\mathbf{n}}_{\mathrm{R}}\right)=\mathit{\xi }·{\mathbf{n}}_{\mathrm{R}}.$$

These macro- and microforce balances, when supplemented with a set of thermodynamically consistent constitutive equations, provide the governing equations for the theory.

We have chosen to relegate a detailed development of the theory to Appendix A. In this Appendix, limiting our considerations to isothermal conditions, we develop a reasonably general theory, and then specialize the constitutive theory to describe isotropic materials. In the main body of the paper we begin by summarizing our isotropic constitutive theory in Section 2. In Section 3 we further specialize the constitutive equations, and discuss a specific set of equations which should be useful for applications. The governing partial differential equations and initial/boundary conditions are summarized in Section 4.

The resulting coupled diffusion–deformation theory is amenable to numerical implementation using standard ${\mathcal{C}}^0\mathrm{-}\mathrm{continuous}$ finite element basis functions. We have numerically implemented our theory by writing a user-element subroutine for the widely used finite element program Abaqus/Standard (2013). Using this numerical capability, in Section 6 we study the diffusion-only problem of spinodal decomposition with the specific task of characterizing the effect of the constant $\beta$ in the split-method used in this work. Then, in Section 7, we use our fully coupled theory and numerical implementation to study the combined effects of diffusion and stress on the lithiation of a representative spheroidal-shaped particle of a phase-separating electrode material.