Coupling Diffusion and Finite Deformation in Phase Transformation Materials

Consider a body occupying a region $\Omega$ in three-dimensional Euclidean space $\mathbb{R}^{3}$ in the reference configuration. Let $\mathbf{x}$ be the position of a point on $\Omega$ in this reference configuration. We describe the deformation of this body as a function $\boldsymbol{\chi}:\Omega\to\mathbb{R}^{3}$, in which $\boldsymbol{\chi}(\mathbf{x},t)$ denotes the position of point $\mathbf{x}$ in the deformed configuration. Given any deformation $\boldsymbol{\chi}$, we describe the deformation gradient $\nabla\boldsymbol{\chi}$ as a matrix of partial derivatives with components,

We use $\mathbf{F}$ to denote this deformation gradient, i.e., a rank-2 tensor $\mathbf{F}=\nabla\boldsymbol{\chi}$ and the gradient operator $\nabla$ is calculated with respect to the reference configuration. Please note, in this work, we denote vectors and tensors using bold lower-case and upper-case Roman letters, respectively. We denote the components of these vectors and tensors, with respect to a Cartesian basis, using upper-case (lower-case) indices in the reference (deformed) configuration. A vector with a hat is a unit normal vector (e.g., $\hat{\mathbf{n}}$) in the reference configuration.

In the continuum theory of crystalline solids, the Cauchy-Born rule, relates the movement of atoms in a crystal to the overall deformation of the body. That is, consider the body $\Omega$ as described above, and at each point $\mathbf{x}\in\Omega$, there is a Bravais lattice that defines the crystalline arrangement of atoms. This Bravais lattice is an infinite set of points in three-dimensional space that can be generated by the translation of a single point $\mathbf{o}$ through three linearly independent lattice vectors $\{\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}\}$. In the reference configuration, the unit cell is defined by lattice vectors $\mathbf{e}^{\circ}_{i}$ and deforms according to the deformation gradient:

We define the Green-Lagrange strain tensor as

in which $\mathbf{I}$ is the identity matrix. The strain tensor in coordinate notation is given by $E_{IJ}=\frac{1}{2}(F_{kI}F_{kJ}-\delta_{{IJ}})$.

The Cauchy-Born rule gives an exact correspondence between these structural transformations of individual lattices and continuum microstructures at a material point. For example, in first-order phase transformation materials (i.e., materials undergoing displacive-type of transformation), we describe the structural transformation of lattices using a positive-definite symmetric matrix $\mathbf{U}$

Fig. 2(a) shows a cubic-to-tetragonal deformation of lattices with the three tetragonal variants described by $\mathbf{U}_{1}$, $\mathbf{U}_{2}$ and $\mathbf{U}_{3}$, respectively. In intercalation materials, such as Li₂Mn₂O₄, a similar displacive-type of lattice transformation is observed using in-situ TEM.

During phase transformation, the tetragonal variants rotate and fit compatibly with each other forming energy-minimizing deformations called twin interfaces. That is, two lattice variants described by transformation matrices $\mathbf{U}_{I}$ and $\mathbf{U}_{J}$ satisfy the Hadamard jump condition (or the Kinematic compatibility condition):

for a given rotation matrix $\mathbf{Q}$ and vectors $\mathbf{a}\neq 0$ and $\mathbf{\hat{n}}$, see Fig. 3. The solution to Eq. (5) then describes a twin plane that connects the two lattice variants $\mathbf{U}_{I}$ and $\mathbf{U}_{J}$ coherently. These twin interfaces have been observed in intercalation materials, such as Li₂Mn₂O₄, in-situ, during battery operation. These twin interfaces are rarely found in isolation. As illustrated in Fig. 3, finely twinned microstructures form when any two lattice variants, with deformation tensors $\mathbf{U}_{I}$ and $\mathbf{U}_{J}$, satisfy the kinematic compatibility condition for some scalar $0\leq f\leq 1$(Ball and James, 1987; Bhattacharya, 2003):

Here the reference cubic lattice is represented by $\mathbf{I}$, rotation matrices by $\mathbf{Q}$ and $\mathbf{Q}^{\prime}$, and vectors $\mathbf{b}\neq 0$ and $\mathbf{\hat{m}}$. These martensite-like microstructures describe energy-minimizing deformations that reduce coherency stresses at the phase boundary and, arise from the multi-well structure of the energy landscape. We discuss this in detail in section 6.1.

In our recent work (Zhang and Balakrishna, 2023), we used compatibility conditions in Eqs. (5) and (6) to analytically construct the twin interfaces in Li_2xMn₂O₄. As shown in Fig. 3, our analytical solutions of both the twin plane $\mathbf{K}$ and the volume fraction $f$ for the austenite-martensite microstructure are consistent with experimental measurements (Erichsen et al., 2020). These observations further support our efforts to derive a diffusion-deformation model to predict crystallographic microstructures in Li₂Mn₂O₄.

Consider any spatial region $\mathcal{P}$ inside the reference body $\Omega$, with an outward normal $\hat{\mathbf{n}}$ and a boundary $\partial\mathcal{P}$. The diffusion of species (e.g., Li-ions in batteries) across this boundary $\partial\mathcal{P}$ is accompanied by a change in the species concentration $c(\mathbf{x},t)$ and characterized by a flux $\mathbf{j}(\mathbf{x},t)$. That is, the rate of change of the chemical species across $\mathcal{P}$ is given by

Using the Divergence theorem over the integral we have,

Note that the choice of $\mathcal{P}$ was arbitrary, and therefore the local mass balance law is as follows

The principle of virtual power states that there exists a fundamental power balance between the external power $W_{\mathrm{ext}}(\mathcal{P})$ expended on $\mathcal{P}$, and the internal power $W_{\mathrm{int}}(\mathcal{P})$ expended within $\mathcal{P}$. To obtain the expressions for external and internal powers, we follow Gurtin (2002) and Anand (2012) and use “rate-like” kinematic descriptors — $\dot{\boldsymbol{\chi}}$, $\dot{\mathbf{F}}$, $\nabla\dot{\mathbf{F}}$, $\dot{c}$, and $\nabla\dot{c}$ — to derive the associated balance of forces by the principle of virtual power. These rates are not independent but are constrained by:

In Table LABEL:Tab:1 we identify individual force systems expending power externally on $\mathcal{P}$ and expending power internally within $\mathcal{P}$. We use these individual force systems to construct the total external $W_{\mathrm{ext}}(\mathcal{P})$ and internal power $W_{\mathrm{int}}(\mathcal{P})$ as follows:

Please note that in Eq. (2.5) the stresses $\mathbf{T}_{\mathrm{R}}$, $\mathbf{Y}$, and microscopic forces $\pi$ and $\boldsymbol{\xi}$ are defined over the body at all times. We define a line force $\mathbf{l}$ across the boundary edge $\zeta^{L}$ of the region $\mathcal{P}$. We assume that at any given time, the fields $\boldsymbol{\chi}$, $\mathbf{F}$, and $c$ are known, and we consider the fields $\dot{\boldsymbol{\chi}}$, $\dot{\mathbf{F}}$, and $\dot{c}$ as virtual velocities. We specify each of these virtual velocities independently and keep them consistent with Eq. (10). Next, we introduce virtual fields of the form $\widetilde{\boldsymbol{\chi}}$, $\widetilde{\mathbf{F}}$, and $\widetilde{c}$, to distinguish these fields from the previously described virtual velocities (as associated with the real-time evolution of the body) and we ensure that virtual fields satisfy:

We define a generalized virtual velocity list $\mathcal{V}=(\widetilde{\boldsymbol{\chi}},\widetilde{\mathbf{F}},\widetilde{c})$ that is consistent with Eq. (12). Following Anand (2012) and Gurtin (2002), we assume that under a change in the frame, the fields comprising the generalized virtual velocity convert as their nonvirtual counterparts. That is,

The external and internal power expenditures, respectively, are:

From Eqs. (14) and (15), the principle of virtual power comprises two essential requirements:

Power Balance: Given any spatial region $\mathcal{P}$,

Frame-Indifference: Given any spatial region $\mathcal{P}$ and any generalized virtual velocity $\mathcal{V}$,

To deduce the outcomes of the principle of virtual power, we assume that Eqs. (16) and (17) are satisfied. We note that in applying the virtual balance described by Eq. (16), we are free to choose any $\mathcal{V}$ that is consistent with the constraint Eq. (12).

We next derive the free energy imbalance under isothermal conditions. Let $\psi$ represent the Helmholtz free energy density of the system per unit volume and $\mu$ represent the chemical potential of the diffusing species in the reference configuration. By neglecting the kinetic energy of the system, we have:

Using the power balance property $W_{\mathrm{ext}}(\mathcal{P})=W_{\mathrm{int}}(\mathcal{P})$, Eq. (2.5) and the divergence theorem, we have:

Using the mass balance from Eq. (9) and noting that Eq. (36) holds for all spatial regions $\mathcal{P}$, we write the local form of the free energy imbalance as:

In Eq. (37), we introduce a net chemical potential,

At this point, we define a dissipation density $\mathcal{D}$ per unit volume per unit time:

Please note that all quantities in Eq. (37) and Eq. (39) are invariant under a change in frame.

We next consider the constitutive forms for the free energy density $\psi$, the first Piola-Kirchhoff stress tensor $\mathbf{T}_{\mathrm{R}}$, the third-order stress tensor $\mathbf{Y}$, the net chemical potential $\mu_{net}$, and the vector microscopic force $\boldsymbol{\xi}$, and the species flux $\mathbf{j}$. Following the free energy imbalance derived in Eq. (37), we note that each term can be expressed as a functional of the deformation gradient $\mathbf{F},\nabla\mathbf{F}$, species concentration $c,\nabla c$, and the chemical potential $\nabla\mu$:

Substituting the constitutive forms in Eq. (2.7) into the free-energy imbalance in Eq. (37), and using the chain rule, we have

This above inequality holds for all values of $\mathbf{F}$, $\nabla\mathbf{F}$, $c$, and $\nabla c$, which, in Eq. (41) appear in a linear form. The corresponding coefficients of these fields must vanish, so that $\dot{\mathbf{F}}$, $\nabla\dot{\mathbf{F}}$, $\dot{c}$, and $\nabla\dot{c}$ can not be chosen to violate the free energy imbalance in Eq. (41). This argument, consequently, leads to the following thermodynamic restriction in which the first Piola-Kirchhoff stress $\mathbf{T}_{\mathrm{R}}$, the third-order stress $\mathbf{Y}$, the chemical potential $\mu$, and the vector microscopic force $\boldsymbol{\xi}$ are defined as derivatives of the free energy function $\hat{\psi}$:

In Eq. (2.7), please note that the symmetry condition for third-order stress $\mathbf{Y}$ is $Y_{iJK}=Y_{iKJ}$. The dissipation inequality introduced in Eq. (39) reduces to

Recalling the definition of $\mu_{net}=\mu+\pi$ in Eq. (38) and by using the constitutive form for $\mu$ in Eq. (2.7) and the microforce balance for $\pi$ in Eq. (34), we derive the constitutive equation for the chemical potential:

We note that the form of the chemical potential in Eq. (44), derived from a microforce balance and thermodynamically consistent constitutive relations, is the same as the chemical potential obtained from a variational derivative of the energy functional (Rudraraju et al., 2016; Zhang and Kamlah, 2018).

Finally, we write the constitutive form of the flux of the diffusing species $\mathbf{j}$ in Eq. (2.7) as:

in which $\hat{\mathbf{M}}$ is a mobility tensor. Substituting Eq. (45) in Eq. (43), we express the dissipation inequality as:

The inequality in Eq. (46) requires the mobility tensor $\hat{\mathbf{M}}$ to be a positive-semidefinite quantity.

We now propose a free energy density for intercalation materials undergoing a symmetry-lowering lattice transformation. This energy density is applicable to other chemo-mechanically coupled phase transformation materials undergoing a displace-type of transformation.

Let us begin with the free energy density $\hat{\psi}(\mathbf{F},\nabla\mathbf{F},c,\nabla c)$ that depends on the lattice deformation gradient $\mathbf{F}$, $\nabla\mathbf{F}$ and the concentration of diffusing species $c$, $\nabla c$. We assume that this free energy density is defined and continuous for all $\mathbf{F}\in\mathrm{M}^{3\times 3}$ and $\mathrm{det}\mathbf{F}>0$ and for all values of the species concentration $c$ during composition evolution. Here $\mathrm{M}^{3\times 3}$ denotes a set of real $m\times n$ matrices. We assume that the free energy is Galilean invariant: for all $\mathbf{F}\in\mathrm{M}^{3\times 3}~{}\mathrm{and}~{}\mathrm{det}\mathbf{F}>0$, for all values of $c$ around the critical point, and each orthogonal rotations $\mathbf{R}$ with $\mathrm{det}\mathbf{R}=1$, we have:

Or equivalently, we describe the free energy in terms of the symmetric Green-Lagrange strain tensor from Eq. (3) $\mathbf{E}=\frac{1}{2}(\mathbf{F}^{\top}\mathbf{F}-\mathbf{I})$, such that:

Furthermore, we assume that the free energy density is an isotropic function of its arguments, and consequently depends only on the magnitude of the gradient terms $\nabla\mathbf{E}$ and $\nabla{c}$,

We note that this free energy density satisfies frame-indifference and material symmetry, i,e., for all rotations $\mathbf{R}\in\mathcal{P}(\mathbf{e}^{\circ}_{i})$ and $\mathcal{P}(\mathbf{e}^{\circ}_{i})$ is a finite point group of the undistorted crystalline lattice. The free energy density $\psi$ therefore satisfies:

The free energy density in Eq. (50) is a higher-order polynomial of its arguments and has a multi-well landscape. This together with the the higher-rank tensors in Eq. (50) introduces nonlinearities making it a challenge to numerically solve the free energy functional. We overcome some of these challenges, by writing the free energy density in terms of a symmetry-adapted strain measure vector $\mathbf{e}=(e_{1},e_{2},e_{3},e_{4},e_{5},e_{6})^{\top}$:

Following (Barsch and Krumhansl, 1984; Thomas and Van der Ven, 2017), we use $\mathbf{e}$ to describe the symmetry-breaking structural transformations of the lattices. The components of this strain measure are in turn a linear combination of the Green-Lagrange strain tensor components. That is, the strain measures $\mathbf{e}=\{e_{1},e_{2},\dots,e_{6}\}^{T}$ are defined using the Green-Lagrange strain tensor components in three dimensions as follows:

We note that within the limits of a small deformation framework, $e_{1}$ represents a stretch-like deformation, and $e_{4}$, $e_{5}$ and $e_{6}$ describe the shear-like deformation of lattices. The strain measures $e_{2}$ and $e_{3}$ in Eq. (3.1) not only assume different values in relation to a high-symmetry cubic lattice and a lower-symmetry tetragonal lattice, but also assume separate values to differentiate among the three tetragonal lattice variants. This construction of $e_{2}$ and $e_{3}$, therefore, is most suitable as structural order parameters describing the cubic-to-tetragonal lattice transformations in intercalation materials.

It follows from Eqs. (50) and (51) that the free energy density satisfying both frame-indifference and material-symmetry is

and holds for all $\mathbf{R}\in\mathcal{P}(\mathbf{e}^{\circ}_{i})$. In this work, we assume a zero vector $\mathbf{e}=0$ minimizes the free energy density in the reference state and a non-zero tensor $\mathbf{e}\neq 0$ minimizes the free energy in the intercalated state. The invariant condition in Eq. (53) implies that if a pair $(\mathbf{e},c)$ minimizes $\psi$, so is $(\mathbf{Re},c)$ for each $\mathbf{R}\in\mathcal{P}(\mathbf{e}^{\circ}_{i})$. These minimizers correspond to the different variants of the lower-symmetry phase and contribute to the multi-well energy landscape of the free energy function.

We next prescribe the specific form of the free energy function. We start by decomposing the free energy density into bulk and gradient energy terms. The bulk energy, in turn, includes contributions from the thermodynamic, elastic, and chemo-mechanically coupled terms:

We construct the thermodynamic energy contribution $\psi_{\mathrm{ther}}(c)$ as a sum of an ideal solution entropy and an excess free energy—representing a deviation from thermodynamic ideality—using a Redlich-Kirster polynomial series (Redlich and Kister, 1948):

Eq. (55) describes an energy landscape, which distinguishes between the reference and intercalated phases, as a function of the normalized concentration $\bar{c}$ (scaled with the maximum species concentration $c_{0}$ as $\bar{c}=c/c_{0}$). The constants, $R$ and $T_{0}$ correspond to the gas constant and the reference temperature, respectively. The Redlich-Kirster coefficients $\alpha_{i}$ are obtained from the least square method and represent the excess energy contribution.

The elastic energy term penalizes the bulk and shear deformations:

The coefficients $\mathit{K}$ and $\mathit{G}$ correspond to the bulk and shear modulus, respectively, and $\Delta V$ represents the volume change associated with the cubic to tetragonal structural transformation of the host lattices.

Following Shchyglo et al. (2012); Barsch and Krumhansl (1984); Ahluwalia et al. (2006), we next construct the coupled chemo-mechanical energy term $\psi_{\mathrm{coup}}(\mathbf{e},\bar{c})$ that describes a multi-well energy landscape in terms of order parameters $\bar{c},e_{2},e_{3}$. This energy landscape governs the cubic to tetragonal lattice transformation:

The coefficient $\beta_{1}(\bar{c})$ represents the concentration-dependent deviatoric modulus governing the cubic to tetragonal transformations. The coefficients $\beta_{2}$ and $\beta_{3}$ are the nonlinear elastic constants. The third-order terms accompanying concentration-dependent $\beta_{2}$ coefficient are required to describe a first-order type of phase-transformation (Cowley, 1976). Fig. 4(a-b) shows three-dimensional plots of the free energy as a function of strain measure components $e_{2}$, $e_{3}$ and normalized concentration $\bar{c}$. The energy well at $(e_{2},e_{3},\bar{c})=(0,0,0.5)$ corresponds to the higher-symmetry cubic phase (denoted by Identity tensor $\mathbf{I}$) and the energy wells at $(e_{2},e_{3},\bar{c})\in\{(0,-0.1,1.0),(-0.1,0.1,1.0),(0.1,0.1,1.0)\}$ corresponds to the tetragonal variants (denoted by stretch tensors $\mathbf{U}_{1},\mathbf{U}_{2},\mathbf{U}_{3}$) at the lower-symmetry phase.

Finally, the gradient energy term penalizes changes in concentration and/or strain and is given by:

In Eq. (58) we only include energy contributions from the quadratic gradient terms that also account for the mixed terms between $\nabla\bar{c}$ and $\nabla\mathbf{e}$. These gradient terms describe nonlocal elastic and composition behavior, and their corresponding coefficients, in their general forms, are tensors that can be functions of both concentration and strain. Specifically, the concentration gradient energy coefficient is denoted by $\boldsymbol{\lambda}$ a symmetric tensor, the strain gradient energy coefficients are denoted by tensors $\boldsymbol{\kappa}^{ij}$ for each combination of $i,j=1,\dots,6$, and mixed gradient energy coefficient(s) $\boldsymbol{\gamma}^{i}$ is a tensor for $i=1,\dots,6$. Eqs. (55)-(58) collectively construct the total free energy density.

We next derive the constitutive equations for the stresses, chemical potential, and microscopic forces for the form of the free energy function in Eq. (54). Using Eqs. (54) and (2.7a) the first Piola-Kirchhoff stress tensor is given by:

We write the third-order stress tensor using Eqs. (2.7b) and (54) as:

We derive the chemical potential expression using Eq. (44) and (54) as follows:

Finally, we construct the vector microscopic force using Eq. (54) and Eq. (2.7d):

We discretize these constitutive equations and implement them in a finite element framework in section 4.

During intercalation a guest-species, such as lithium, is inserted into the host-material lattices. This insertion is accompanied by an electrochemical reaction, at the electrode/electrolyte interface, which can be described for any guest species ‘A’ as follows:

We incorporate this electrochemical reaction in our constitutive model by using a phenomenological Butler-Volmer equation (Bai et al., 2011; Mykhaylov et al., 2019; Ganser et al., 2019) and describing the reaction rate. The Butler-Volmer reaction relates the current density $i$ to surface overpotential $\eta$ at the interface between electrode and electrolyte as:

In Eq. (64) $\beta$ is the electron-transfer symmetry factor, and $F$ is the Faraday constant. The exchange current density $i_{0}$ given by:

with $k_{0}$ denoting the reaction rate constant (units of $\mathrm{mol/m^{2}s}$) and $\mu_{+}$ denoting the chemical potential in the electrolyte. We assume that the guest-species moves much faster in the electrolyte than in the active host material and therefore set $\mu_{+}=0$ (Bai et al., 2011). We define the surface overpotential $\eta$ as

where $\Delta\phi$ is the voltage drop across the interface between electrode and electrolyte. This voltage drop serves as a driving force for species insertion (or removal) into the electrode.

By combining Eqs. (64)-(66) we write the species flux at the electrode/electrolyte interface as: