Fracture Resistance of Chemo-Mechanically Coupled Solid Solutions

Solid-state diffusion occurs ubiquitously in many material systems, such as rechargeable battery electrodes, fuel cells for energy conversion, and metal alloys. In these solid solutions, mass transport and diffusion often generate substantial stress, which can lead to the fracture and ultimate failure of the materials. For the rational design and development of these material systems, the combined effects of chemical and mechanical driving forces on fracture must be thoroughly understood.

A general thermodynamic framework accounting for stress–diffusion coupling in multicomponent solids was established by Larché and Cahn [1]. Using a stress-dependent chemical potential as one of the key ingredients, this framework was initially applied to metallic systems [2–5] and later extended to study ionic crystals [6,7], nonstoichiometric oxide ceramics [8–11], and oxide-metal interfaces [12,13]. More recently, the chemomechanics of lithium-ion battery (LIB) electrodes has been a topic of intensive research, driven by the escalating electrification in various sectors. The operation of a LIB involves repeated insertion and extraction of mobile lithium (Li) ions in electrode materials, which causes appreciable volume changes. In particular, many alloying- and conversion-type electrode materials, such as silicon (Si), germanium (Ge), tin, and oxide and sulfide compounds, undergo large-volume-change transformations (up to 300%) during lithiation [14–17]. If left unaddressed, this large deformation can lead to mechanical failure and consequently rapid capacity fading of the electrodes which compromises their high-capacity advantage [18–20].

Extensive experimental and modeling work has recently been conducted to gain the fundamental understanding of chemo-mechanical interplay in the context of LIB electrodes. The multibeam optical sensor (MOS) and nanoindentation techniques have been employed to characterize the stress evolution [21], stress-potential coupling [22], and stress–diffusion coupling of lithiated thin-film electrodes [23,24]. Christensen and Newman are among the first to model the coupling between mechanical stresses and Li diffusion using a mathematical model for diffusion-induced stresses [25,26]. A variety of computational continuum models were subsequently developed to account for more general material behaviors, such as finite deformation, elastoplasticity, and phase segregation [27–35]. Atomistic mechanisms underpinning diffusive transport and plastic deformation of Li_xSi nanostructures have also been revealed using ab initio calculations [36,37] and molecular dynamics (MD) simulations [38–40]. In addition, analytical approaches have been taken to derive closed-form relationships to characterize stress buildup in electrode materials [41–43].

Fracture caused by diffusion-induced stresses within solid electrodes has been analyzed with several classical fracture criteria, including Griffith’s criterion [44,45], stress intensity factor criterion [46–48], and nonlinear energy release rate (ERR) (J-integral) criterion [49–51]. Grantab et al. [52] and Bower et al. [53] used a cohesive zone approach to model crack nucleation and propagation in electrode materials by assuming composition-independent fracture properties. Employing large-scale atomistic simulations, Ding et al. investigated the fracture behavior of amorphous lithiated Si [54]. Their work revealed a brittle fracture at low Li concentrations due to void nucleation and coalescence and ductile fracture at high Li concentrations through extensive shear banding. The fracture toughness of lithiated Si versus composition and charging/discharging history was further studied by Khosrownejad et al. through K-controlled MD simulations [55] and was found to be consistent with the Gurson model for ductile fracture. More recently, Xu et al. investigated the corrosive fracture of electrodes in LIBs by implementing a theory of coupled diffusion and large deformation [56]. The theory was applied to model intergranular fracture in LiNi_xMn_yCo_zO₂ (NMC) cathode particles during delithiation and Li embrittlement was found to accelerate the initiation and growth of intergranular cracks in NMC particles.

In a solid solution, the coupling of crack-tip mechanical and chemical fields can have a profound influence on fracture. Stress concentration near a crack tip modifies the gradient of chemical potential which drives solute diffusion. The resulting change in solute concentration, in turn, can alter the crack-tip deformation and stress fields. Despite its importance, the intriguing interaction between mass diffusion and crack growth in solid solutions remains understudied and poorly understood. In this work, we develop a comprehensive computational framework for modeling the fracture of chemo-mechanically coupled solids. In Sec. 2, we present this framework which features an integration of a composition-dependent cohesive zone model with a two-way, chemo-mechanical coupling continuum model. Numerical implementation of this framework is explained in Sec. 3. In Sec. 4 we apply our computational framework to investigate the crack growth resistance of solid solutions and the fracture behavior of lithiated nanostructures. Finally, we conclude in Sec. 5.

The above finite element scheme is implemented through the user element (UEL) interface of ABAQUS. The bulk material is modeled using three-dimensional (3D) 20-node brick elements (Fig. 1(c)). The eight corner nodes have displacement, concentration, and chemical potential DOFs, and the 12 edge nodes only have displacement DOFs. Second-order isoparametric shape functions are chosen for displacement interpolation, and linear isoparametric shape functions are used to interpolate the concentration and chemical potential fields. The cohesive zone is modeled using 16-node cohesive zone elements, which have zero initial thickness and share nodes with the adjacent bulk elements. Like the bulk elements, the cohesive zone elements combine quadratic interpolation for the crack opening displacement and linear interpolation for the other two fields. This hybrid interpolation scheme is used here to provide an optimal balance between computational accuracy and efficiency.

The computational framework outlined in the preceding sections has been applied to investigate a model problem, where a plane strain specimen containing an initial edge crack is subjected to mode I loading (Fig. 1(a)). Through this model problem, we seek to gain fundamental insights into how the fracture behavior of a solid solution is controlled by the interplay between the mechanical and chemical fields around a crack tip. The specimen has a width of W, a height of H = 0.3W, an initial crack length of a₀ = 0.3W, and a uniform mesh size of L = W/200. The solid solution is initially stress free with a homogeneous distribution of guest species concentration c₀, which corresponds to an initial intrinsic fracture energy of $Gc0=Gc(c0)$⁠. For simplicity, an isotropic linear elastic material behavior is assumed and Young’s modulus and Poisson’s ratio are taken to be independent of the concentration. Due to the geometrical and loading symmetry about the crack plane, only the upper half of the specimen is simulated. Mode I fracture loading is applied by prescribing a pair of opening displacements with constant speed v_app on the upper and lower vertical surfaces to the left of the crack. As the fracture loading increases, a cohesive crack initiates and propagates with a nonlinear fracture process zone at the crack tip, whose constitutive behavior is governed by the traction-separation law described earlier. The path-independent chemo-mechanical J-integral calculation (Sec. 2.4) is performed on the cracked specimen to evaluate the energy release rates at various stages of loading, which are then used to construct the fracture resistance curve (R-curve) of J_∞ versus Δa of the material.

In quasi-static fracture simulations, snap-back instability in the load–displacement path is often observed because of the excessive rate of change in crack-tip ERR (i.e., dJ_tip/da > dG_c/da), leading to numerical convergence issues. With the displacement-controlled loading condition in the current model, the ERR decreases with crack growth, which helps maintain crack stability. However, as will be shown later, the resistance curve of the material is strongly influenced by chemo-mechanical coupling, which may lead to crack instability even under displacement control. To overcome this issue, a small artificial viscosity is introduced in the cohesive zone law [61,62]. This viscosity term provides an energy dissipation mechanism for capturing post-instability behaviors, while making a negligible contribution to the fracture resistance during stable crack propagation.

Systematic parametric studies have been carried out to study the effects of two controlling parameters, the crack growth velocity v_c and the concentration dependence of intrinsic fracture energy dG_c/dc. Our dimensional analysis yields a characteristic length $l0=Gc0/E$⁠, a characteristic time $t0=(Gc0)2/(E2D)$⁠, and a characteristic velocity $v0=ED/Gc0$⁠, which are used for dimensionless representation of the relevant quantities. Also, the concentration dependence of intrinsic fracture energy is normalized by its initial value $Gc0$ to be expressed as a dimensionless parameter $χ=(1/Gc0)(dGc/dc)$⁠.

Figure 2(a) shows a plot of crack advance versus time for three distinct loading velocities (v_app = 0.000265, 0.0265, and 2.65) and a representative value of χ = 3.33. The crack length is seen to increase linearly with time in the low- and high-velocity cases, but this linear relationship does not hold very well for the intermediate velocity. Linear fitting of each curve yields an average crack growth velocity v_c, which takes the place of v_app as a more intrinsic fracture condition. The R-curves for the three cases are shown in Fig. 2(b). For the high and low crack velocities (v_c/v₀ = 24.7 and 0.0019), J_∞ is found to increase abruptly and reach a plateau (maximum) value with little crack extension. In sharp contrast, the J_∞ curve for the intermediate crack velocity (v_c/v₀ = 0.19) fluctuates as the crack first starts to grow, and eventually converges to a plateau. In the following, the plateau value of each J_∞ curve will be interpreted as the effective fracture resistance $Gceff$ under the specific fracture condition.

To elucidate the fracture mechanism underlying the different resistance curves in Fig. 2(b), we resort to the steady-state mechanical and chemical fields around the crack tip, which are obtained when the plateau values of J_∞ are reached. Figure 3 shows the spatial distributions of concentration change Δc = c − c₀, hydrostatic stress σ_h, and crack opening stress component σ_yy. When the crack propagates at the low velocity (v_c/v₀ = 0.0019), the gradient in hydrostatic stress near the crack tip drives the guest species to diffuse toward the crack tip and cause solute accumulation around the tip, as evidenced by the contour plot a1 in Fig. 3(a) and the line plot in Fig. 3(b). The influx of the guest species can influence the effective fracture resistance in two separate ways. First, it causes the local area in the wake of the crack tip to expand. Such dilatational deformation lowers the local stress level of the crack-tip region, thus producing an effective shielding from the far-field mechanical loading (see a2 in Figs. 3(a) and 3(c)). We note that this toughening mechanism closely resembles the concept of stress-induced phase transformation in shape memory alloys [63] and zirconia-containing ceramics [64]. From an energetic point of view, this toughening effect can also be attributed to the energy dissipation by diffusion in the near-crack-tip region. Second, since the fracture energy of the solid solution is in general dependent on the chemical composition, the change in crack-tip concentration directly alters the intrinsic resistance to fracture.

For fast crack growth (FCG) (v_c/v₀ = 24.7), the time taken for crack growth is much smaller than the characteristic time for solute diffusion. Therefore, no appreciable amount of diffusion takes place in the crack-tip region, leading to no significant change in solute concentration (see c1 in Figs. 3(a) and 3(b)). In this case, the solute is essentially immobile and the fracture process is purely mechanically driven. This renders an effective fracture resistance equal to the intrinsic fracture energy of the material.

For the intermediate crack velocity (v_c/v₀ = 0.19), the crack growth and diffusion processes occur on comparable time scales, which give rise to a more complex fracture behavior than the two previous cases. Figure 4 shows the contour and line plots of Δc and σ_h at four successive time instants that are marked as A–D on the R-curve in Fig. 2(b). Upon initial loading, the stress singularity at the crack tip creates a localized region of high solute concentration (instant A). When the crack starts to propagate, the crack growth velocity is so large that solute diffusion cannot catch up, causing the high-concentration site to lag behind the crack tip (instant B). This hysteretic phenomenon is more evident as the crack further propagates, and a second high-concentration site forms at the moving crack tip and the original one starts to fade away (instant C). The negotiation of the propagating crack with solute diffusion results in the nonconstant crack velocity shown in Fig. 2(a) and the oscillating crack growth resistance curve as seen in Fig. 2(b). Beyond time instant D, a dynamic equilibrium between the chemical and mechanical driving forces for diffusion is reached and maintained, giving rise to a plateau fracture resistance value between those of the low- and high-crack-velocity cases.

In Fig. 5(a) we present a parametric study of the dependence of effective fracture resistance $Gceff$ on crack growth velocity and χ. The characteristic velocity v₀, which is set by the ratio of diffusivity to intrinsic fracture energy, serves as an important parameter for dictating the effective fracture resistance. Depending on the crack velocity relative to this characteristic velocity, there are three distinct fracture mechanism regimes: the slow crack growth (SCG) regime, the FCG regime, and the intermediate crack growth (ICG) regime in between. For a fixed value of χ, $Gceff$ is found to maintain a relatively constant value in the SCG regime (v_c/v₀ < 0.002). As v_c/v₀ increases, all the $Gceff$ curves first rise in the ICG regime, and then converge nonmonotonically to the same asymptotic value in the FCG regime (v_c/v₀ > 2) that is equal to the initial intrinsic fracture energy $Gc0$⁠. To understand the complex behavior of these $Gceff$ curves, we consider the decomposition of $Gceff$ into the crack-tip fracture energy G_c and the diffusion-induced energy dissipation ΔJ, which are presented in Figs. 5(b) and 5(c). In the SCG and ICG regimes, the solute atoms have sufficient time to diffuse and redistribute around the crack tip. The diffusion may lead to enhanced $(Gc>Gc0)$ or reduced $(Gc<Gc0)$ crack-tip fracture energy, depending on the sign of the χ value. However, the diffusion-induced ΔJ is always positive and therefore contributes to fracture toughening as discussed earlier. It is worth noting that the ΔJ curves in Fig. 5(c) are maximized in the ICG regime and decay toward zero for both low- and high-crack velocities, giving rise to the nonmonotonic dependence of $Gceff$ on crack velocity. The trend of ΔJ → 0 as v_c → 0 may seem counterintuitive, considering that a low crack velocity promotes diffusion in the near-crack-tip region. Nevertheless, this result can be rationalized by the fact that diffusion-induced dilatational eigenstrains around a stationary crack do not affect the energy release rate (i.e., ΔJ = 0). Such assertion has previously been proved in the context of dilatant transformation toughening in ceramics [65].

The fracture mechanism map in Fig. 5(a) has profound implications for materials design and characterization of chemo-mechanically coupled solid solutions. In a solid solution where the solute concentration induces toughening, crack growth velocities in the ICG regime are favorable for attaining high effective damage tolerance of the material. On the other hand, if a material is susceptible to solute embrittlement, crack velocities beyond the SCG regime can mitigate the embrittlement effect. For fracture testing of solid solutions, crack growth velocities maintained within the FCG regime can prevent an appreciable amount of solute diffusion around the crack tip and therefore allow for the accurate measurement of intrinsic fracture properties. For example, typical diffusivity values for most solid solutions range from 10⁻¹² to 10⁻¹⁶ m²/s. Assuming representative values of E = 100 GPa and $Gc0=10J/m2$⁠, we obtain characteristic crack velocities ranging from 0.001 to 10 mm/s. While crack growth velocities higher than these are relatively easy to attain for macroscopic fracture specimens, fracture testing at the nano- and micro-scale is usually performed at fairly slow loading rates due to instrumental constraints. Therefore, care must be taken in these fracture experiments regarding the FCG requirement to properly interpret measurement results.

To further illustrate the utility of our computational framework, we have also applied it to study the fracture behavior of Si nanoparticles (SiNPs) and Ge nanoparticles (GeNPs) upon lithiation. In our earlier works [66,67], the fracture properties of electrochemically lithiated Si and Ge at various Li concentrations were quantitatively measured by nanoindentation. Electrodes for nanoindentation testing were fabricated in a thin-film form and electrochemically lithiated to various Li concentrations. To measure the fracture toughness, the Li_ξSi and Li_ξGe electrodes were fracture tested with a cube-corner indenter tip. The indentation load and crack sizes were used to evaluate the fracture toughness of Li_ξSi and Li_ξGe with the Morris nanoindentation model [68,69]. Figure 6(a) shows the measured fracture energies of Li_ξSi and Li_ξGe as functions of Li concentration, which clearly reveal that the fracture energies of Li_ξSi and Li_ξGe alloys depend sensitively on the Li concentration. As the Li concentration increases, the fracture energy of Li_ξSi first decreases slightly, indicating lithiation-induced embrittlement. However, upon further lithiation beyond ξ = 0.31, the measured fracture energy increases substantially with increasing Li concentration. Unlike Li_ξSi, the obtained fracture energy of Li_ξGe is shown to increase monotonically with the increase of Li concentration. These experimentally characterized fracture properties of Li_ξSi and Li_ξGe are incorporated into the computational model in this work to improve its predictive capability.

Figure 7 presents the simulated fracture behavior of lithiated Si and Ge nanoparticles using our computational framework. An initial crack is placed at the edge of each nanoparticle configuration. Three particle diameters, 100 nm, 250 nm, and 500 nm, are used for both Si and Ge nanoparticles in the comparative study, with the initial crack length being 4% of the particle diameter in each simulation. Due to the geometry and loading symmetry about the crack plane, only the upper half of each nanoparticle is modeled. A linear rule of mixtures is used to determine the composition-dependent mechanical properties, including the yield strength and elastic modulus, of lithiated Si and Ge. The lithiation of Si and Ge has been experimentally observed to proceed with the movement of a sharp phase boundary between pristine and lithiated regions [70–72], as shown in Fig. 6(b). This two-phase lithiation mechanism is implemented into our model by adding a double-well potential [73] to the chemical potential described in Eq. (18) and assuming ξ₁ = 0 and ξ₂ = 2.33 as two favorable states (Fig. 6(c)). As shown in Fig. 7, SiNPs exhibit size-dependent fracture upon lithiation, with no cracking observed for a small particle diameter of 100 nm (Fig. 7(a)) but obvious fracture for particle diameters of 250 nm and 500 nm (Figs. 7(b) and 7(c)). In distinct contrast to SiNPs, GeNPs undergo no visible cracking for all the three diameters investigated (Figs. 7(d)–7(f)). These simulation results conform quantitatively to the in situ transmission electron microscopy (TEM) observations reported by previous works [74,75], indicating that GeNPs are more robust than their Si counterparts and can be a promising electrode candidate for durable, high-rate rechargeable batteries.

In this work, we integrate a composition-dependent cohesive zone model with a diffusion–deformation coupling continuum model to investigate the fracture mechanics of chemo-mechanically coupled solid solutions. We apply our computational framework to analyze the crack growth resistance of a model solid-solution system. We show that the dynamic interaction between stress-induced diffusion and crack propagation dictates the fracture response of the system. A chemo-mechanical J-integral method allows the effective and crack-tip fracture energies to be calculated separately. For sufficiently high-crack velocities, both energy terms converge to the intrinsic fracture energy. However, when the crack grows at lower velocities, stress-induced diffusion can lead to enhancement or weakening in the crack-tip fracture resistance depending on the composition dependence of the intrinsic fracture energy. Diffusion-induced dilatational eigenstrains are found to contribute positively to the effective fracture resistance. Such a contribution is most significant when the crack velocity is comparable to a characteristic velocity set by the ratio of diffusivity to intrinsic fracture energy. We also use our computational method to study the fracture behavior of lithiated Si and Ge nanoparticles. The results reveal that SiNPs exhibit size-dependent fracture, while GeNPs experience no obvious cracking under similar lithiation conditions. The demonstrated mechanical robustness of GeNPs make them a promising electrode candidate for durable rechargeable batteries.

The authors gratefully acknowledge the support of the National Science Foundation (CMMI-1300458 and CMMI-1554393).

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.