## Abstract

Silicon-based anodes are one of the promising candidates for the next generation high-power/energy density lithium ion batteries (LIBs). However, a major drawback limiting the practical application of the Si anode is that Si experiences a significant volume change during lithiation/delithiation, which induces high stresses causing degradation and pulverization of the anode. This study focuses on crack initiation within a Si anode during the delithiation process. A multi-physics-based finite element (FE) model is built to simulate the electrochemical process and crack generation during delithiation. On top of that, a Gaussian process (GP)-based surrogate model is developed to assist the exploration of the crack patterns within the anode design space. It is found that the thickness of the Si coating layer, TSi, the yield strength of the Si material, σFc, the cohesive strength between Si and the substrate, σFs, and the curvature of the substrate, ρ, have large impacts on the cracking behavior of Si. This coupled FE simulation-GP surrogate model framework is also applicable to other types of LIB electrodes and provides fundamental insights as building blocks to investigate more complex internal geometries.

## 1 Introduction

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 (LiCoO2 and LiMn2O4) 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 LixSi 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−1 for the Li22Si5 alloy phase; meanwhile, the theoretical capacity of graphite is ∼370 mA h g−1 (with a fully lithiated graphite phase of LiC6) [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 [911], 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 [1214]. 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 [1618]. 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 TSi, 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.

## 2 Modeling Approach

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.

### 2.1 The Governing Equations.

The delithiation process of LixSi is considered as a single-phase reaction in this model. The lithium concentration profile of the anode is determined by Fick’s law (1). When considering a thin Si layer as a host material, the effects of the hydrostatic stress gradient on the diffusion flux, J, have also been taken into consideration:
$J=−D(∇c−ηcRT∇p)$
(1)
where c is the concentration of lithium, D is the diffusion coefficient, η is the expansion coefficient, R is the gas constant, T is the absolute temperature, ∇ is gradient operator, and p = tr(σ) is the hydrostatic pressure. Further combining Fick’s law (1) with the mass conservation equation, $∂c/∂t+∇⋅J=0$, we have
$∂c∂t−D∇(∇c−ηcRT∇p)=0$
(2)
The boundary condition of Eq. (2) at the Si-electrolyte surface is controlled by the electric current density in
$J⋅n=in/F$
(3)
where n is the normal vector of the anode-electrolyte surface. Generally, periodic boundary conditions are applied to other surfaces for approximating a large Si anode structure.
Afterwards, the lithium concentration profile is used to simulate the delithiation-induced volume reduction and the crack initiation process. The volumetric shrinkage of Si caused by delithiation is analogous to the 3D thermal contraction process [1921], so that the stress and deformation of Si during delithiation are correspondingly treated as the thermal contraction of a material surrounded by a low temperature media. An isotropic volume contraction of amorphous Si is assumed. Therefore, the total deformation of Si is composed of the elastic deformation that obeys Hooke’s law, the diffusion-induced term, and perfectly plastic deformation term:
$εij=εel+εLi+εpl=1E[(1+ν)σij−νσkkδij]+cΩ3δij+εpl$
(4)
where ɛij and σij are the components of strain and stress tensors, respectively. δij is the Kronecker delta, ɛel, ɛLi, and ɛpl are the elastic, lithiation-induced, and plastic strain, respectively. Ω is the partial molar volume (m3/mol). E and v are Young’s modulus and Poisson’s ratio, respectively. The lithium diffusion (Eq. (2)) and the stress–strain (Eq. (4)) are coupled through the lithium concentration term c.

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.

Figure 1 shows the schematic of the cohesive model. The CZM presets a cohesive surface Γc that permits the material planes adjacent to the surface to separate and uses a rate-independent, history-dependent, and irreversible traction-separation law that corelates the cohesive traction vector, T, and the displacement jump vector, Δ, along the cohesive surface to evaluate the resistance to crack nucleation and growth. In the tensile separation mode (mode I), which is the focus of the work, this cohesive model has the form
$T=ζ1−ζΔΔCσmax$
(5)
where T and Δ are the cohesive traction and crack opening displacement vectors, respectively. σmax is the tensile cohesive failure strength. ΔC is the critical value of the opening displacement jump, beyond which complete delamination is regulated. ΔC is codetermined by σmax and the fracture toughness G1C of the cohesive surface by
$G1C=12σmaxΔC$
(6)
Fig. 1
Fig. 1
Close modal
ζ is a coefficient within [0, 1), is updated every timestep, and is defined as
$ζ=min(ζp,⟨1−Δ/ΔC⟩)$
(7)
where ζp denotes the ζ value from the previous timestep. 〈a〉 represents a function where 〈a〉 = a if a > 0, 〈a〉 = 0 otherwise. ζ monotonically decreases over time and is used to quantify the real-time damage state and prevent the healing of the de-bonded/failed interfaces over subsequent loading/unloading cycles. The initial value of ζ is set at around 1 (set as 0.98 in this study).

### 2.2 Finite Element Model Implementation.

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 (LiPF6) 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 cmax of 3.367 × 105 mol/m3 based on the theoretical specific capacity 3579 mA h/g of Li3.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.

Table 1

Electrochemical parameters used in the model

ParametersValue
Lithium diffusion coefficient in electrolyte DLi2.6 × 10−10 m2/s [26]
Lithium transference number t+0.363 [26]
Lithium diffusion coefficient in amorphous Si DLi–Si1 × 10−16 m2/s [27]
Si maximum capacity cmax3.367 × 105 mol/m3 [28]
ParametersValue
Lithium diffusion coefficient in electrolyte DLi2.6 × 10−10 m2/s [26]
Lithium transference number t+0.363 [26]
Lithium diffusion coefficient in amorphous Si DLi–Si1 × 10−16 m2/s [27]
Si maximum capacity cmax3.367 × 105 mol/m3 [28]

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, Fc, and the slipping force of the Si–substrate interface, Fs, [7]. If Fc is exceeded by the delithiation-induced stress at a certain cohesive surface of ΓSi, a crack is initiated; correspondingly, if Fs is exceeded at ΓInterface, the Si layer will delaminate from substrate and move/slip along with adjacent cells. In the CZM model, Fc and Fs are evaluated via the cohesive strength, σmax, and fracture toughness, G1C, 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.

Fig. 2
Fig. 2
Close modal

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 ENi = 200 GPa, Poisson’s ratio vNi = 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.

Fig. 3
Fig. 3
Close modal

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.

## 3 Gaussian Process Surrogate Modeling

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].

An ordinary GP model generally has a form of
$G(x)≈GGP(x)=fT(x)⋅α+S(x)$
(8)
where G(x) and GGP(x) are the true and the GP prediction of the performance of the system at point x, respectively. fT(x) = [f1(x), …, fb(x)] is a basis function, α = [α1, …, αb] is a regression coefficient vector, and S(x) is a Gaussian stochastic process with mean equal to zero and variance equal to σ. The covariance function between two inputs xi and xj is expressed as
$Cov(i,j)=σ2⋅exp[−∑p=1kap|xip−xjp|bp]$
(9)
where ap and bp are parameters of the GP model, and k is the number of the input parameters x. The GP model is then trained with the imported simulation results (XS, GS) to generate the initial battery performance function G(x). With the GP model, the performance $G^(x′)$ of a new point x′ can be estimated via maximizing the likelihood function [32] as
$G^(x′)=μ+rTR−1(G−Aμ)$
(10)
$e^(x′)=σ2[1−rTR−1r+(1−ATR−1r)2/ATR−1A]$
(11)
where $e^(x′)$ is the variation and r is the correlation vector between x′ and the input parameters x.

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 (XS, GS) are generated by the FE model to develop the initial GP model, where XS are the model input parameters and GS is the corresponding output performance. In this study, a random sampling strategy is used to select the input XS. (2) In the next step, the simulation results (XS, GS) are stored in a sample data pool and used to construct the predictive models M employing the GP-based surrogate modeling technique [3234]. (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.

Fig. 4
Fig. 4
Close modal

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 (Fc) or the coated Si material will delaminate from the substrate (Fs). 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.

Table 2

Main design variables of the silicon anode and their ranges

Design variableInitial designDesign range
k = Fc/Fs0.50.1 to 2.0
TSi (nm)20020 to 650
Mean curvature (µm−1)0 (Flat)−0.6 to 0.6
Design variableInitial designDesign range
k = Fc/Fs0.50.1 to 2.0
TSi (nm)20020 to 650
Mean curvature (µm−1)0 (Flat)−0.6 to 0.6

## 4 Results and Discussion

### 4.1 Electrochemical Modeling Results.

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.

Fig. 5
Fig. 5
Close modal

### 4.2 Crack Initiation on the Flat Substrate.

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.

Fig. 6
Fig. 6
Close modal

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.

Fig. 7
Fig. 7
Close modal

Figure 8 plots the area of the Si island with respect to the thickness of the coating layer, TSi. A scaling relationship is observed, indicating that the area, A, is exponentially related to TSi: $A∝TSin$. n is approximated as 1.3 in this model. The exponential relationship between A and TSi 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 TSi, WSi, and LSi are the thickness, width, and length of the Si island, respectively. Therefore, along with the increase of TSi, LSi 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.

Fig. 8
Fig. 8
Close modal

In addition to TSi, 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.

Fig. 9
Fig. 9
Close modal

### 4.3 The Impact of the Curved Substrate.

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−1. The two insets are examples of the Si anodes with negative (−0.5 µm−1) and positive curvatures (0.5 µm−1), 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.

Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal

### 4.4 Gaussian Process Surrogate Model Predictions.

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, TSi, 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 TSi 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 TSi or curvature ρ, the Si island area also increases.

Fig. 12
Fig. 12
Close modal
Fig. 13
Fig. 13
Close modal

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.

Fig. 14
Fig. 14
Close modal

## 5 Conclusion

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, TSi, 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 TSi: $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.

## Acknowledgment

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.

## References

1.
Tarascon
,
J. M.
, and
Armand
,
M.
,
2001
, “
Issues and Challenges Facing Rechargeable Lithium Batteries
,”
Nature
,
414
(
6861
), pp.
359
367
. 10.1038/35104644
2.
McDowell
,
M. T.
,
Lee
,
S. W.
,
Nix
,
W. D.
, and
Cui
,
Y.
,
2013
, “
25th Anniversary Article: Understanding the Lithiation of Silicon and Other Alloying Anodes for Lithium-Ion Batteries
,”
,
25
(
36
), pp.
4966
4984
3.
Armand
,
M.
, and
Tarascon
,
J. M.
,
2008
, “
Building Better Batteries
,”
Nature
,
451
(
7179
), pp.
652
657
. 10.1038/451652a
4.
Limthongkul
,
P.
,
Jang
,
Y.-I.
,
Dudney
,
N. J.
, and
Chiang
,
Y.-M.
,
2003
, “
Electrochemically-Driven Solid-State Amorphization in Lithium-Silicon Alloys and Implications for Lithium Storage
,”
J. Acta Mater.
,
51
(
4
), pp.
1103
1113
. 10.1016/S1359-6454(02)00514-1
5.
Kasavajjula
,
U.
,
Wang
,
C.
, and
Appleby
,
A. J.
,
2007
, “
Nano- and Bulk-Silicon-Based Insertion Anodes for Lithium-Ion Secondary Cells
,”
J. Power Sources
,
163
(
2
), pp.
1003
1039
. 10.1016/j.jpowsour.2006.09.084
6.
Boukamp
,
B. A.
,
Lesh
,
G. C.
, and
Huggins
,
R. A.
,
1981
, “
All-Solid Lithium Electrodes With Mixed-Conductor Matrix
,”
J. Electrochem. Soc.
,
128
(
4
), pp.
725
729
. 10.1149/1.2127495
7.
Li
,
J.
,
Dozier
,
A. K.
,
Li
,
Y.
,
Yang
,
F.
, and
Cheng
,
Y.-T.
,
2011
, “
Crack Pattern Formation in Thin Film Lithium-Ion Battery Electrodes
,”
J. Electrochem. Soc.
,
158
(
6
), pp.
A689
A694
. 10.1149/1.3574027
8.
Chew
,
H. B.
,
Hou
,
B. Y.
,
Wang
,
X. J.
, and
Xia
,
S. M.
,
2014
, “
Cracking Mechanisms in Lithiated Silicon Thin Film Electrodes
,”
Int. J. Solids Struct.
,
51
(
23–24
), pp.
4176
4187
. 10.1016/j.ijsolstr.2014.08.008
9.
Shi
,
F. F.
,
Song
,
Z. C.
,
Ross
,
P. N.
,
Somorjai
,
G. A.
,
Ritchie
,
R. O.
, and
Komvopoulos
,
K.
,
2016
, “
Failure Mechanisms of Single-Crystal Silicon Electrodes in Lithium-Ion Batteries
,”
Nat. Commun.
,
7
, p.
11886
.
10.
Yang
,
H.
,
Fan
,
F.
,
Liang
,
W.
,
Guo
,
X.
,
Zhu
,
T.
, and
Zhang
,
S.
,
2014
, “
A Chemo-Mechanical Model of Lithiation in Silicon
,”
J. Mech. Phys. Solids
,
70
, p.
349
361
. 10.1016/j.jmps.2014.06.004
11.
Liu
,
X. H.
,
Zheng
,
H.
,
Zhong
,
L.
,
Huang
,
S.
,
Karki
,
K.
,
Zhang
,
L. Q.
,
Liu
,
Y.
,
Kushima
,
A.
,
Liang
,
W. T.
, and
Wang
,
J. W.
,
2011
, “
Anisotropic Swelling and Fracture of Silicon Nanowires During Lithiation
,”
Nano Lett.
,
11
(
8
), pp.
3312
3318
. 10.1021/nl201684d
12.
Zhang
,
H. G.
, and
Braun
,
P. V.
,
2012
, “
Three-Dimensional Metal Scaffold Supported Bicontinuous Silicon Battery Anodes
,”
Nano Lett.
,
12
(
6
), pp.
2778
2783
. 10.1021/nl204551m
13.
Li
,
J.
,
Wang
,
J.
,
Yang
,
J.
,
Ma
,
X.
, and
Lu
,
S.
,
2016
, “
Scalable Synthesis of a Novel Structured Graphite/Silicon/Pyrolyzed-Carbon Composite as Anode Material for High-Performance Lithium-Ion Batteries
,”
J. Alloys Compd.
,
688
, pp.
1072
1079
. 10.1016/j.jallcom.2016.07.148
14.
Yang
,
J.
,
Wang
,
Y.-X.
,
Chou
,
S.-L.
,
Zhang
,
R.
,
Xu
,
Y.
,
Fan
,
J.
,
Zhang
,
W.-X.
,
Liu
,
H. K.
,
Zhao
,
D.
, and
Dou
,
S. X. J. N. E.
,
2015
, “
Yolk-Shell Silicon-Mesoporous Carbon Anode With Compact Solid Electrolyte Interphase Film for Superior Lithium-Ion Batteries
,”
Nano Energy
,
18
, pp.
133
142
.
15.
Ng
,
S. S. Y.
,
Xing
,
Y.
, and
Tsui
,
K. L.
,
2014
, “
A Naive Bayes Model for Robust Remaining Useful Life Prediction of Lithium-Ion Battery
,”
Appl. Energy
,
118
, pp.
114
123
. 10.1016/j.apenergy.2013.12.020
16.
Meng
,
Z.
,
Zhang
,
D.
,
Liu
,
Z.
, and
Li
,
G.
,
2018
, “
An Adaptive Directional Boundary Sampling Method for Efficient Reliability-Based Design Optimization
,”
ASME J. Mech. Des.
,
140
(
12
), p.
121406
. 10.1115/1.4040883
17.
Fan
,
X.
,
Wang
,
P.
, and
Hao
,
F.
,
2019
, “
Reliability-Based Design Optimization of Crane Bridges Using Kriging-Based Surrogate Models
,”
J. Struct. Multidiscipl. Optim.
,
59
, pp.
993
1005
.
18.
Zheng
,
Z.
,
Chen
,
B.
,
Gurumukhi
,
Y.
,
Cook
,
J.
,
Ates
,
M. N.
,
Miljkovic
,
N.
,
Braun
,
P. V.
, and
Wang
,
P.
,
2019
, “
Surrogate Model Assisted Design of Silicon Anode Considering Lithiation Induced Stresses
,”
2019 IEEE International Reliability Physics Symposium (IRPS)
,
Monterey, CA
,
Mar. 31–Apr. 4
.
19.
Park
,
J.
,
Lu
,
W.
, and
Sastry
,
A. M.
,
2011
, “
Numerical Simulation of Stress Evolution in Lithium Manganese Dioxide Particles Due to Coupled Phase Transition and Intercalation
,”
J. Electrochem. Soc.
,
158
(
2
), pp.
A201
A206
. 10.1149/1.3526597
20.
Verbrugge
,
M. W.
, and
Cheng
,
Y. T.
,
2009
, “
Stress and Strain-Energy Distributions Within Diffusion-Controlled Insertion-Electrode Particles Subjected to Periodic Potential Excitations
,”
J. Electrochem. Soc.
,
156
(
11
), pp.
A927
A937
. 10.1149/1.3205485
21.
Cheng
,
Y. T.
, and
Verbrugge
,
M. W.
,
2009
, “
Evolution of Stress Within a Spherical Insertion Electrode Particle Under Potentiostatic and Galvanostatic Operation
,”
J. Power Sources
,
190
(
2
), pp.
453
460
. 10.1016/j.jpowsour.2009.01.021
22.
De-Andrés
,
A.
,
Pérez
,
J.
, and
Ortiz
,
M.
,
1999
, “
Elastoplastic Finite Element Analysis of Three-Dimensional Fatigue Crack Growth in Aluminum Shafts Subjected to Axial Loading
,”
Int. J. Solids Struct.
,
36
(
15
), pp.
2231
2258
. 10.1016/S0020-7683(98)00059-6
23.
Deshpande
,
V. S.
,
Needleman
,
A.
, and
Van der Giessen
,
E.
,
2001
, “
A Discrete Dislocation Analysis of Near-Threshold Fatigue Crack Growth
,”
Acta Mater.
,
49
(
16
), pp.
3189
3203
. 10.1016/S1359-6454(01)00220-8
24.
Maiti
,
S.
, and
Geubelle
,
P. H.
,
2006
, “
Cohesive Modeling of Fatigue Crack Retardation in Polymers: Crack Closure Effect
,”
Eng. Fract. Mech.
,
73
(
1
), pp.
22
41
. 10.1016/j.engfracmech.2005.07.005
25.
Bower
,
A. F.
, and
Guduru
,
P. R.
,
2012
, “
A Simple Finite Element Model of Diffusion, Finite Deformation, Plasticity and Fracture in Lithium Ion Insertion Electrode Materials
,”
Modell. Simul. Mater. Sci. Eng.
,
20
(
4
), p.
045004
. 10.1088/0965-0393/20/4/045004
26.
Xu
,
K.
,
2004
, “
Nonaqueous Liquid Electrolytes for Lithium-Based Rechargeable Batteries
,”
Chem. Rev.
,
104
(
10
), pp.
4303
4417
. 10.1021/cr030203g
27.
Ding
,
N.
,
Xu
,
J.
,
Yao
,
Y. X.
,
Wegner
,
G.
,
Fang
,
X.
,
Chen
,
C. H.
, and
Lieberwirth
,
I.
,
2009
, “
Determination of the Diffusion Coefficient of Lithium Ions in Nano-Si
,”
Solid State Ionics
,
180
(
2–3
), pp.
222
225
. 10.1016/j.ssi.2008.12.015
28.
Cho
,
Y. H.
,
Booh
,
S.
,
Cho
,
E.
,
Lee
,
H.
, and
Shin
,
J.
,
2017
, “
Theoretical Prediction of Fracture Conditions for Delithiation in Silicon Anode of Lithium Ion Battery
,”
APL Mater.
,
5
(
10
), p.
106101
. 10.1063/1.4997978
29.
Beaulieu
,
L. Y.
,
Eberman
,
K. W.
,
Turner
,
R. L.
,
Krause
,
L. J.
, and
Dahn
,
J. R.
,
2001
, “
Colossal Reversible Volume Changes in Lithium Alloys
,”
Electrochem. Solid State Lett.
,
4
(
9
), pp.
A137
A140
. 10.1149/1.1388178
30.
Feifei
,
F.
,
Shan
,
H.
,
Hui
,
Y.
,
Muralikrishna
,
R.
,
Dibakar
,
D.
,
Vivek
,
B. S.
,
,
C. T. V. D.
,
Sulin
,
Z.
, and
Ting
,
Z.
,
2013
, “
Mechanical Properties of Amorphous Li x Si Alloys: A Reactive Force Field Study
,”
Modell. Simul. Mater. Sci. Eng.
,
21
(
7
), p.
074002
. 10.1088/0965-0393/21/7/074002
31.
Rubinstein
,
R. Y.
, and
Kroese
,
D. P.
,
2016
,
Simulation and the Monte Carlo Method
, Vol.
10
,
John Wiley & Sons
,
Hoboken, NJ
.
32.
Wang
,
Z. Q.
, and
Wang
,
P. F.
,
2014
, “
A Maximum Confidence Enhancement Based Sequential Sampling Scheme for Simulation-Based Design
,”
ASME J. Mech. Des.
,
136
(
2
), p.
021006
. 10.1115/1.4026033
33.
Wang
,
P. F.
,
Wang
,
Z. Q.
, and
Almaktoom
,
A. T.
,
2014
, “
Dynamic Reliability-Based Robust Design Optimization With Time-Variant Probabilistic Constraints
,”
Eng. Optim.
,
46
(
6
), pp.
784
809
. 10.1080/0305215X.2013.795561
34.
Wang
,
Z. Q.
, and
Wang
,
P. F.
,
2015
, “
A Double-Loop Adaptive Sampling Approach for Sensitivity-Free Dynamic Reliability Analysis
,”
Reliab. Eng. Syst. Saf.
,
142
, pp.
346
356
. 10.1016/j.ress.2015.05.007
35.
Pikul
,
J. H.
,
Braun
,
P. V.
, and
King
,
W. P.
,
2017
, “
Performance Modeling and Design of Ultra-High Power Microbatteries
,”
J. Electrochem. Soc.
,
164
(
11
), pp.
E3122
E3131
. 10.1149/2.0151711jes
36.
Zheng
,
Z. Y.
,
Chen
,
B.
,
Fritz
,
N.
,
Gurumukhi
,
Y.
,
Cook
,
J.
,
Ates
,
M. N.
,
Miljkovic
,
N.
,
Braun
,
P. V.
, and
Wang
,
P. F.
,
2019
, “
Lithiation Induced Stress Concentration for 3D Metal Scaffold Structured Silicon Anodes
,”
J. Electrochem. Soc.
,
166
(
10
), pp.
A2083
A2090
. 10.1149/2.1031910jes
37.
Chan
,
C. K.
,
Peng
,
H.
,
Liu
,
G.
,
McIlwrath
,
K.
,
Zhang
,
X. F.
,
Huggins
,
R. A.
, and
Cui
,
Y.
,
2008
, “
High-Performance Lithium Battery Anodes Using Silicon Nanowires
,”
Nat. Nanotechnol.
,
3
(
1
), p.
31
35
. 10.1038/nnano.2007.411
38.
Kim
,
H.
, and
Cho
,
J.
,
2008
, “
Superior Lithium Electroactive Mesoporous Si@Carbon Core−Shell Nanowires for Lithium Battery Anode Material
,”
Nano Lett.
,
8
(
11
), pp.
3688
3691
. 10.1021/nl801853x