A Gaussian Process-Based Crack Pattern Modeling Approach for Battery Anode Materials Design

Secondary (rechargeable) lithium ion batteries (LIBs) outperform other battery chemistries because of their high energy/power density and good cycle life performances [1,2], and they have become the key energy storage devices for portable electronic devices ranging from smart phones and tablet PCs to electric vehicles (EVs). Currently, commercial LIBs commonly use lithium metal oxides or phosphates (LiCoO₂ and LiMn₂O₄) as the cathode and graphite as the anode. However, these chemistries are approaching their performance limits, unable to meet increasing demands for high gravimetric/volumetric power/energy and cycle life, as well as lower cost and improved safety [3], especially for emerging EV systems.

To address these challenges, numerous efforts have been made to develop new electrodes with enhanced lithium storage capacity. Silicon is a promising candidate for the active anode material in LIBs. In contrast to the intercalation mechanism of graphite for Li-ion storage, i.e., Li-ion diffusion into the interstitial sites of the host lattice, Si will react with lithium, leading to bond-breaking between Si atoms and the formation of a Li_xSi alloy accompanied by dramatic structural changes [4]. With no constraints from the atomic framework of the host material, Si is able to store much more lithium compared to intercalation electrodes [5]. For instance, Si has the highest known theoretical specific capacity of 4200 mA h g⁻¹ for the Li₂₂Si₅ alloy phase; meanwhile, the theoretical capacity of graphite is ∼370 mA h g⁻¹ (with a fully lithiated graphite phase of LiC₆) [6].

However, due to this “alloying” lithium storage mechanism, the volume of silicon changes significantly (up to 300%) during lithiation/delithiation cycles, which creates large internal stresses and cracks, and leads to the pulverization of anode and reduction of capacity [5]. Various works have been conducted to study the crack initiation behaviors of the Si active material due to the significant volume changes during cycling. Li et al. [7] developed a modified spring-block model to capture the essential features of cracking patterns in amorphous Si thin film electrodes during the delithiation process. They found that the crack patterns were related to the thickness of the Si layer, where thicker films could lead to the straight cracks, whereas thinner films were more favorable for wiggles cracks. Chew et al. [8] investigated the crack initiation mechanisms of the amorphous Si film during the lithiation/delithiation cycles through a two-dimensional finite element (FE) model. A critical film thickness, below which the crack initiation could be mitigated, was predicted and showed a good agreement with experiments. Besides, the crack initiation and growth in the single-crystal Si electrodes through the electrode thickness direction during delithiation were analyzed using experiments and FE models [9–11], and the influences of anisotropic electrochemical/mechanical properties were explored. However, most of the previous works considered the Si anode with relatively straightforward structures, e.g., flat films, nanowires, etc. With widespread improvements to electrode fabrication methods, a number of novel structured Si anodes with designed 3D supporting scaffolds as current collectors are proposed [12–14]. Thus, fundamental understandings regarding various structural features of the sophisticated electrodes/scaffolds (i.e., curvature of the scaffold substrate, cohesive strength between the Si layer and the substrate, etc.) and their synergetic effects on the cracking initiation behaviors need to be investigated.

Aside from physics-based simulation methods, data-driven techniques [15] are also widely used for electrode design optimization. Compared to physics-based models, which are built to explore certain failure mechanisms of a specific battery, data-driven techniques can effectively and efficiently capture the general changing trends of the electrode properties affected by various design parameters. The drawbacks of data-driven models, however, are that they usually require massive experimental results as training data, which can be expensive to acquire; additionally, these models could not reveal the underlying physical mechanisms. This study presents an integrated platform that is mixing the physics-based simulation and the Gaussian process (GP) surrogate modeling technique to investigate the crack patterns within the anode design space. The GP surrogate modeling technique has been extensively utilized in engineering applications because of its accuracy and efficiency [16–18]. One of the major contributions of the GP surrogate model is that, compared to experimental methods or physics-based numerical simulation approaches, it can efficiently explore the multi-dimensional design parameter space of the anode, estimate its performances with respect to continuous or discrete design variables, and provide optimal design strategies with identification of the fidelity of the surrogate model.

In this study, a multi-physics-based FE model is built to investigate the crack pattern of Si anode and the critical design parameters, including thickness of the Si layer T_Si, the cohesive strength between the Si layer and the substrate, and the curvature of the substrate. Subsequently, a GP-based surrogate model is developed with the simulated results as training data to predict the cracking behavior of the anodes within the anode design space. The GP surrogate model results are then used to analyze the general trends of the Si anode performance.

A multi-physics-based 3D FE model is established via comsol multiphysics to simulate the crack patterns of Si anodes with complex structures during the delithiation process. The model consists of two sub-models. In the first part, the electrochemical dynamics of the Si anode during delithiation, a lithium diffusion process, is simulated; subsequently, the lithium concentration profile is acquired and used to evaluate the crack pattern based on a standard cohesive zone modeling (CZM) technique. Finally, a GP-based surrogate model with adaptive fidelity enhancement is implemented to analyze the influences of the design parameters of the Si anode structure, such as the Si layer thickness, the curvature of the substrate, etc.

The CZM technique is utilized to investigate the crack initiation of Si during volume shrinkage. Recently, CZM methods have been widely used to analyze the crack initiation and growth in metals [22,23] and polymers [24,25]. In the following, the basic features of the model are outlined.

The governing equations are then implemented in the 3D FE model software (comsol multiphysics). The multi-physics simulation model consists of two sub-modules, the electrochemical (lithium battery) model and the mechanical (crack initiation) model, which are solved sequentially. The electrochemical model is used to calculate the lithium concentration, c, with respect to the delithiation process. Then, c is extracted from the simulation results and passed to the crack initiation model. The model simulates the Si half-cell battery, with pure lithium metal and Si as the two electrodes. One molar lithium hexafluorophosphate (LiPF₆) in 1:1 ethylene carbonate:diethyl carbonate, one of the most commonly used liquid electrolytes, is used in the model. The lithiation-induced expansion is simulated analogously to the thermal expansion process with lithium concentration, c, as a substitute for temperature in the model to acquire stresses and deformations of the anode. Table 1 lists the electrochemical and mechanical parameters used in the simulation. The lithium diffusion coefficients in electrolyte and silicon are kept constant. Silicon active material has a maximum capacity c^max of 3.367 × 10⁵ mol/m³ based on the theoretical specific capacity 3579 mA h/g of Li_3.75Si [28], which is considered as the fully lithiated state in this study. The open circuit voltage of Si versus lithium is measured as a function of the state of charge (SoC) from slow rate discharging of half cells.

To study the initiated crack patterns, a few assumptions are made. A Si thin layer is coated on a metal support substrate. The substrate is also used as a current collector of the electrode to transport electrons (e⁻) to the external circuit. The cracks are considered to be generated only during the delithiation process due to the volume shrinkage of the Si layer and propagate vertically (perpendicular to the substrate) [7,29]. In addition, it is assumed that the cracks do not change with the increase of cycle numbers since the volume change-induced stress is not large enough to generate new cracks [7]; therefore, the cyclic effect is not considered in this work. Figure 2 illustrates the schematic of the crack initiation model with a flat substrate. Note that the substrate can be changed to a curved one with designated curvature. The coated Si layer is divided into small rectangular volume cells with a length and width of 10 × 10 nm. The height of the cells represents the thickness of the Si coating layer and is treated as a controllable input parameter. Note that for a curved substrate, the Si cells exhibit trapezoidal-prism-like shapes with the side walls (height of the cells) perpendicular to the substrate, and top and bottom faces having identical curvatures as the underlying substrate. Each of the individual cells is bonded to adjacent ones obeying CZM to represent the Si material. Additionally, the Si layer is also adhered to the substrate. Periodic boundary conditions are applied to all four vertical surfaces to mimic bulk Si. During delithiation, cracks can form inside this bulk Si along the cohesive surfaces Γ_Si; meanwhile, delamination may also occur at the Si–substrate interface Γ_Interface, resulting in a slipping of Si cells. Therefore, the overall crack initiation and propagation is the result of competition between two crucial forces (fracture criteria): the cracking force of Si, F_c, and the slipping force of the Si–substrate interface, F_s, [7]. If F_c is exceeded by the delithiation-induced stress at a certain cohesive surface of Γ_Si, a crack is initiated; correspondingly, if F_s is exceeded at Γ_Interface, the Si layer will delaminate from substrate and move/slip along with adjacent cells. In the CZM model, F_c and F_s are evaluated via the cohesive strength, σ_max, and fracture toughness, G_1C, where the $σFc$ and $G1C−Fc$ of the cracking force are those of the yield strength and toughness of pure silicon; and $σFs$ and $G1C−Fs$ of the slipping force reflect the properties of the Si–substrate interface. The ratio $k=σFc/σFs$ is also considered as a controllable variable in this study to investigate the crack initiation properties. After full delithiation, cracks may be generated along either Γ_Si or Γ_Interface to form several Si islands. The areas of the Si islands separated by cracks are chosen as the main output and used to analyze the crack initiation and propagation performances.

Annealed nickel is chosen as the current collector substrate. A Ni substrate endures elastoplastic deformation. The mechanical properties of nickel (annealed) is obtained from online resources (MatWeb) with Young’s modulus E_Ni = 200 GPa, Poisson’s ratio v_Ni = 0.30, and yield strength σ_Ni = 600 MPa, (assuming perfect plasticity). For the silicon active material, the modulus decreases in an approximately linear manner with increasing lithium concentration, leading to significant elastic softening. The non-constant Young’s modulus and Poisson’s ratio of silicon with respect to the lithiation state are shown in Fig. 3 [30]. Silicon bears elastoplastic deformation, with a yield strength $σy=σFc=170MPa$ and toughness $G=G1C−Fc=15J/m2$⁠. Silicon is assumed to have perfect plasticity with no hardening effect. The Si layer and Ni substrate are in contact with no inter-penetration.

Note that, by calibrating the mechanical parameters of the active material and substrate with experiments or theoretical calculation, the model can be widely applied to thin film LIB electrodes with other active materials.

The multi-physics-based FE model introduced above is computationally expensive, especially when the design space is high dimensional. To further improve efficiency without losing accuracy, a GP-based surrogate model is introduced. GP surrogate model is a statistical procedure that generates an estimated surface from a scattered set of points via a covariance governed Gaussian process interpolation method. It requires training samples for model construction and predicting responses of new sample points [31].

The fidelity of initial GP model results G(x) depends on the sampling scenario in the electrode design space and can be relatively low. In order to improve the fidelity of the model, the maximum expected information value (MEIV)-based adaptive sampling method is applied. Figure 4 illustrates the flowchart of developing the GP model with enhanced adaptive fidelity. In general, three consecutive steps are involved in the approach: (1) an initial set of sample points (X_S, G_S) are generated by the FE model to develop the initial GP model, where X_S are the model input parameters and G_S is the corresponding output performance. In this study, a random sampling strategy is used to select the input X_S. (2) In the next step, the simulation results (X_S, G_S) are stored in a sample data pool and used to construct the predictive models M employing the GP-based surrogate modeling technique [32–34]. (3) The most important sample points x* are then selected using the MEIV-based adaptive sampling technique. The performances of these points G* are evaluated with the FE model, and the new sample points (x*, G*) are imported into the data pool to update the surrogate model. (4) Cumulative confidence level (CCL) [32] is used to qualify the accuracy of reliability of the surrogate model. Eventually, the updated GP model acquires enhanced fidelity, which can be used for electrode performance prediction and optimization design.

Table 2 lists the main design variables and their ranges considered in this study. The ratio $k=σFc/σFs$ controls the two failure mechanisms of the anode structure, i.e., whether the Si layer will break and generate cracks (F_c) or the coated Si material will delaminate from the substrate (F_s). With $σFc$ as an intrinsic property of the Si material and kept as constant (170 MPa), k varies along with $σFs$ from 0.1 to 2.0, where smaller k indicates that Si is more vulnerable to fracture, and larger k implies that the anode has lower resistance to delamination. The influence of the thickness of the Si coating layer on crack initiation performance is also explored from 20 to 650 nm. Last but not least, the impact of flat and curved substrates is investigated, with the mean curvature of the substrate spanning from −0.6 to 0.6.

The lithium concentration profile during the delithiation process is simulated via the electrochemical sub-model. The delithiation process is assumed to be conducted under constant current with a slow delithiation rate (C/5) to minimize the influence of non-uniform lithium concentrations in the Si host material. Figure 5(a) illustrates the simulated SoC of the Si anode during the delithiation process. With constant delithiation current, the lithium concentration decreases in a linear manner with respect to time, driven by the lithium flux outwards from the electrode. The SoC of Si along the thickness direction of the thin anode layer after delithiation is illustrated in Fig. 5(b), where the thickness Ṫ_Si is a non-dimensionalized value for a clear label. As shown in the inset of Fig. 5(b), Ṫ_Si = 0 represents the Si–substrate contact interface Γ_Interface, and Ṫ_Si = 1 denotes the Si-electrolyte surface. It can be seen that the lithium concentration gradient along the thin anode layer is rather small. This is mainly because of the short lithium diffusion pathway inside the thin Si layer (at nanoscale), large diffusivity of lithium in Si, and slow discharging rate. Consequently, the small lithium concentration gradient could reduce the inner stress of Si, lead to a uniform stress distribution inside the Si layer, and thus improve the mechanical properties of the Si anode.

Figure 6 illustrates the representative simulation results of crack patterns (top and perspective view) in Si layers with different thicknesses coated on a flat substrate after delithiation/volume shrinkage. The red-dash-circled regions in Figs. 6(a) and 6(b) are the Si islands separated by cracks. Although, due to the modeling assumptions, the Si islands all have square shapes which are not realistically accurate comparing with experiments, the model can still capture the general trends of crack initiation behaviors and thus be used to analyze the impacts of critical design parameters. It is observed that, with the increase of the Si layer thickness, the area of the Si island also increases. The changes of the areas of the Si island are crucial to the reliability performances of Si anode: smaller Si island area indicates more cracks in the Si layer, which could introduce more sites for the formation of solid electrolyte interphase (SEI) and more intensive degradation of the anode.

Figure 7(a) shows the gap between adjacent Si cells upon delithiation, where the cells are adhered to the substrate and shrinking. It can be observed that the gap distance at the top surface is ∼2.3 nm and keeps decreasing downwards; at the bottom, the gap reduces to 0, meaning the adjacent cells are still connected at this timestep and the crack is not completely formed yet. This result suggests the process of the formation of cracks, as illustrated in Fig. 7(b): (1) during delithiation, the Si layer undergoes a tensile stress induced by the volume reduction; (2) in the beginning, as the deformation starts to grow, the tensile stress σ_Tensile increases correspondingly (Eq. 4); (3) when it exceeds the yield strength of Si $σy=σFc$⁠, the crack is initiated from the top surface; (4) as the crack approaches the bottom, the propagation of the gap is retarded by the adhesion interface $σFs$⁠; and (5) eventually, when deformation-induced stress σ_Tensile overcomes the synergistic effect of $σFc$ and $σFs$⁠, the Si layers will be thoroughly separated by the crack.

Figure 8 plots the area of the Si island with respect to the thickness of the coating layer, T_Si. A scaling relationship is observed, indicating that the area, A, is exponentially related to T_Si: $A∝TSin$⁠. n is approximated as 1.3 in this model. The exponential relationship between A and T_Si is also reported elsewhere [7]. The causes of this phenomenon can be explained as follows. With a slow deformation rate (discharge rate C/5), the crack propagation can be considered static. By applying the static force equilibrium condition (Fig. 7(b)), we have $σTensile×(TSi×WSi)=σFs×(LSi×WSi)$⁠, where T_Si, W_Si, and L_Si are the thickness, width, and length of the Si island, respectively. Therefore, along with the increase of T_Si, L_Si also grows, i.e., the area of the Si island increases. Besides, the power n is largely influenced by the mechanical data assigned in the model, including modulus of the active material, yield strength σ_y, and toughness of Si, as well as the adhesive properties at Si–substrate interfaces. Hence, with changes of these parameter values, the value of n also varies.

In addition to T_Si, the impact of the adhesive properties of the Si–substrate interfaces is also investigated. Figure 9 shows the relationship between the Si island area and the fracture criteria ratio k, where the thickness of the Si layer is fixed at 200 nm, and $k=σFc/σFs$ is controlled by changing the cohesive strength $σFs$⁠. A second-order polynomial relationship is found between Si island area A and ratio k. With the decrease of the ratio k, the area of Si islands also decreases. This is because, smaller k means the cohesive strength $σFs$ is enhanced, so that the Si layer becomes more tightly adhered on the substrate; in other words, the Si layer becomes relatively vulnerable to cracking. As a result, to absorb the same amount of strain energy caused by the Si volume shrinkage, systems with small k need to form more cracks, leading to the reduction of the Si island area.

As mentioned above, with the improvement of electrode fabrication technology, much more sophisticated Si anode structures have been prepared in recent years. One of the common features among these novel electrodes is that the current collector is usually elaborately designed and fabricated into a three-dimensional bi-continuous supporting scaffold, so that a thin layer of the Si active material can be conformally coated on the scaffold [12,35]. The major advantages of this kind of design are that these architectures achieve high power density by simultaneously reducing the ion and electron transport lengths; meanwhile, they provide large porosity that allows the Si material to expand without causing large-scale pulverization. Among these electrodes, however, the Si layer must often be coated on a curved substrate rather than a simple flat one. Therefore, in this part, the influences of the substrate curvature are studied using the 3D finite element analysis (FEA) model.

Figure 10 illustrates the area of Si islands with respect to the changes in the curvature of the substrate from −0.6 to 0.6 µm⁻¹. The two insets are examples of the Si anodes with negative (−0.5 µm⁻¹) and positive curvatures (0.5 µm⁻¹), respectively. The thickness of the Si layer is fixed at 200 nm, and ratio k is kept as 0.5. Positive curvature leads to a reduction of the Si island area, suggesting an exacerbated fracture scenario. On the other hand, when the curvature is negative and continues decreasing, larger area of the Si island is found, indicating that substrates with negative curvature tend to resist the formation of cracks. This can be explained as follows. As mentioned previously in this model, the cracks are assumed to only propagate perpendicularly to the substrate. Therefore, as shown in Fig. 11, for the negatively curved substrate, the Si cells show a trend of reduced areas from bottom to top. In this case, the deformation-induced stress, especially in the beginning stages of delithiation, is very small and even becomes negative/compressive. Consequently, the formation of cracks, which initiates from the top surfaces of the Si layer as introduced in Fig. 7(a), is restrained, resulting in a larger area of the Si island. What is more, according to our previous works [18,36], it is found that with the same Si–substrate interfacial properties, i.e., the same cohesive strength and fracture toughness, the substrate with a negative curvature tends to be more vulnerable to delamination compared to the flat substrate, which indirectly diminishes the formation of cracks and benefits the integrality of the Si layer. On the contrary, with the positively curved substrate, the Si top surfaces would experience more intensified delithiation-induced deformation and larger tensile stresses, which leads to an aggravated crack initiation. The Si layer ends up with more cracks, dividing it into more islands with smaller areas.

The GP-based surrogate model is used to predict the Si island area within the design space. FEA-based simulation results with different Si anode design parameters, i.e., fracture ratio, k, Si layer thickness, T_Si, and curvature of the substrate, ρ, are used as training data to develop the surrogate model. Additionally, the CCL [32] is adopted to qualify the accuracy and reliability of the surrogate model. In order to have a clear comparison, the Si cracking performance is analyzed using two performance surfaces as illustrated in Figs. 12 and 13. Figure 12 shows the predicted cracking areas of the Si layer on a flat substrate (curvature ρ = 0) with respect to T_Si and ratio k. Within the whole parameter space, a monotonically increasing surface is generated from the point at (20 nm, 0.1) to the (600 nm, 2.0) corner. The performance surface is not linear with respect to the two variables: at the (400 nm, 1.5) point, a small concave region is observed. This is due to the non-linear relationships between the Si island area and the design variables, as illustrated in Figs. 8 and 9. Similarly, as illustrated in Fig. 13, a relatively smooth performance surface is observed. By increasing T_Si or curvature ρ, the Si island area also increases.

The GP surrogate modeling can be particularly useful to investigate the Si anodes with complex structures. Recently, with emerging electrode fabrication technologies, various complicated designs of the Si anode have been proposed, such as nanowires [37,38], porous inverse opal structures [12], etc. It is challenging and time consuming to study the Si cracking phenomenon in these sophisticated designs and find the optimal design through experimental approaches or physics-based numerical simulation methods. The GP-based surrogate model, therefore, is an effective alternative method to evaluate the crack generations in those Si anodes. First, researchers may acquire the statistical distributions of the critical geometric features of the Si anode through experimental measurements, such as the thickness of the Si layer, the curvature of the substrate, etc. Then, one can refer to the GP surrogate model results as lookup tables to quickly predict the Si island area distribution and the crack generation in the anodes. Through this approach, enormous experimental efforts and computational costs can be saved, which will largely benefit the design of the novel Si anode. In the future work, we will use the GP surrogate model results to investigate the Si cracking in the novel inverse opal structured Si anode.

Figure 14 illustrates the reliability (CCL value) of the GP surrogate model with increasing of data points used to train the model. It can be seen that the reliability of the model gradually improves with the increase in numbers of the sample points and reaches as high as 0.997 when 20 training data points are used. This demonstrates that the developed GP surrogate model can rather precisely predict the performance of the Si anode within the design space range and can be further used for design optimization.

In this work, we studied the crack initiation performances of Si anode during the delithiation process. Multi-physics-based simulations coupled with the GP-based surrogate model are utilized to investigate the formation of cracks in a Si layer and the critical factors that largely influence the cracking behaviors. By applying a constant current discharging profile with a slow delithiation rate, the lithium concentration in the Si active material decreases linearly, and the concentration gradient along the Si layer thickness direction is small, which is thus neglected by assuming a uniform lithium concentration in Si. The thickness of the Si coating layer, T_Si, the yield strength $σFc$ of the Si material, the cohesive strength between Si and the substrate, $σFs$⁠, and the curvature of the substrate, ρ, are studied as the factors that affect the formation of cracks. During the simulation of delithiation, cracks initiate on the surface of Si due to its volume contraction, propagate toward the substrate, and eventually divide the Si layer into separated Si islands. The area of the Si island is selected as the main output and analyzed, because it is directly related to the reliability performances of Si anode: smaller Si island area indicates more cracks in the Si layer, which could introduce more sites for the formation of SEI and more intensive degradation of the anode. The area of the Si island A is found to have a scaling relationship with T_Si: $A∝TSin$⁠. Besides, the cohesive strength $σFs$ of the Si–substrate interface shows large impact on the formation of cracks: an increase of cohesive strength $σFs$ could lead to larger Si island area. The curvature of the substrate plays an important role in the crack initiation behaviors of Si. A substrate with negative mean curvature could prevent the generation of crack; a positively curved substrate, on the other hand, exacerbates crack formation. The GP-based surrogate model is further implemented to predict the cracking properties in the design space. GP surrogate model is found to be a powerful tool to effectively and efficiently explore the design space and assist further optimization design. It is worth mentioning that, by adjusting and calibrating the parameters of the models, this finite element simulation-coupled GP surrogate model framework can be extensively used to analyze the properties of the electrodes of other active materials and with complicated geometry structures.

This research is partially supported by the National Science Foundation through the Faculty Early Career Development award (CMMI-1813111), the NSF award (CMMI-1802489) for the Gaussian Process-based surrogate modeling, and the Office of Naval Research (ONR) through the Navy and Marine Corps Department of Defense University Research-to-Adoption (DURA) Initiative (N00014-18-S-F004) for the Si-based battery anode study.