Three-dimensional peridynamic model to predict fracture evolution during lithiation

: Due to its large electric capacity, silicon has become one of the most promising electrode materials for lithium ion batteries. However, silicon undergoes large volumetric expansion and material stiffness reduction during the charging process. This can lead to fracture and failure of lithium-ion batteries. Damage formation and evolution inside the electrode are inﬂuenced by the lithium ion concentration and electrode material. High stress gradients induced by heterogeneous deformation can lead to massive migration of lithium ions towards high geometrical singularity regions, such as crack edge regions, which increases the lithium ion concentration. Fully coupled mechanical diffusion equations are important in describing the mechanics of this problem. In this study, the three-dimensional peridynamic theory is presented to solve the coupled ﬁeld problem. In addition, the newly developed peridynamic differential operator concept is utilized to convert partial differential equations into peridynamic form for the diffusion equation. Spherical and cylindrical shaped energy storage structures with different pre-existing penny-shaped cracks are considered to demonstrate the capability of the developed framework. It is shown that peridynamic theory is a suitable tool for predicting crack evolution during the lithiation process.


Introduction
In the marine industry, traditional propulsion systems with internal combustion engines cannot match the increasing strict standards for waste gas emissions (MARPOL Annex VI) [1]. Hence, the hybrid propulsion system, as an alternative approach to marine propulsion systems, has become popular in the shipping industry. Marine batteries, as primary energy storage devices, have been subject to excessive investigation [2]. The lithium ion battery is one of the most promising marine batteries due to its high energy to weight ratio, high energy density, low rate of self-discharge and long cycle life [3]. The performance of a lithium ion battery mainly depends on the material properties of its anode, cathode and electrolyte. The carbon group elements, carbon (C), silicon (Si), germanium (Ge), tin (Sn) and their alloys, can be selected as the anode material in lithium ion batteries [4,5]. In many commercial grade lithium ion batteries, graphitic carbon is the main component of the anode material due to its low expansion induced by lithiation during battery charging [4]. However, the major limitation of graphite carbon is its low electric capacity (372 mAhg −1 ) whereby it stores the lithium ion inside graphite sheets as (LiC 6 ) [4][5][6] which cannot satisfy the demand for high electric capacity in marine batteries. On the other hand, silicon (Si) and germanium (Ge) have relatively high theoretical capacities as compared to graphite carbon (3579 mAhg −1 and 1625 mAhg −1 , respectively) [6] which can be considered to be suitable candidates for the anode material in lithium ion batteries. However, both Si and Ge anodes experience large volumetric changes during battery cycling. The Si anode can experience a volume expansion of up to 400% during the charging process [7,8]. Similar to the Si anode, the Ge anode can expand 370% in volume during the lithiation process [9]. Frequent cycling of the lithium ion battery may lead to stress misdistribution, degradation and delamination of the battery components. This may then lead to failure or even pulverization of the battery anode.
Many efforts have been conducted to increase the electric capacity and to avoid the occurrence of defects inside the electrodes of lithium ion batteries. Winter et al. [10] suggested that some metals, such as aluminum (Al), tin (Sn) and antimony (Sb) can store more lithium ion than graphite carbon by forming alloys with lithium ions. However, the newly formed alloyed anode will experience a 500% volumetric expansion [4], which is the major barrier to becoming the anode material in rechargeable lithium ion batteries. Gao et al. [11] suggested carbon nanotubes as the anode material. Carbon nanotubes are tubular forms of graphite sheet with high conductivity, high tensile strength, high rigidity and low density. They increase the electric capacity of carbon by up to around 600 mAhg −1 without damage or pulverization. Metals or their alloys can also be added into the nanotube in composite form. Hence, the advantage of the high electric capacity of metal alloys and the low volumetric change in graphite carbon can be utilized in anode materials [12,13]. However, carbon nanotubes experience a loss of lithium ion capacity during cycling and a linear voltage drop during discharging [4].
In comparison with graphite carbon or carbon nanotube anodes, silicon has superior performance in terms of electrical capacity. However, as mentioned earlier, it suffers from pulverization due to large volumetric change during battery cycling. The fracture behaviour of silicon anodes has been subject to frequent investigation in recent years. Liu et al. [14] constructed a thin silicon film model to investigate lithiation-induced tensile stress and surface cracking by using analytical and finite element methods (FEM). A compression-traction transition zone was observed along the interface of electrode and electrolyte. The formation of the transition zone depends on the large volumetric change, plastic deformation and slow charging rate. Crack formation and propagation was led by the location of the transition zone [14,15]. The magnitude and profiles of tensile stress at the surface of the lithiated silicon zone depends on volumetric misfit strain, yield stress and modulus of unlithiated silicon. Ryu et al. [8] found that the large volume change during the normal cycling process is always accompanied by a pressure gradient. High hydrostatic stress will affect the lithium ion diffusion and fracture propagation inside the anode of a battery. Grantab and Shenoy [16] conducted a detailed investigation of the effect of the pressure gradient factor on crack propagation in silicon nanowires. The cohesive FEM was applied to model fracture evolution in silicon nanowires. Since localized stress around the crack tip is higher than the surrounding nanowire surface region, a large quantity of lithium ions moved into the crack tip region which caused a relatively large volume expansion. Therefore, the hydrostatic stress around the crack tip reduced. Zuo and Zhao [17] used the phase field method to study the stress evolution and crack propagation. Several damaged anode models with different crack numbers and different crack orientations were considered to illustrate fracture evolution in the battery anode. The pressure gradient was shown to be dependent on the elastic modulus, partial molar volume, concentration, Poisson's ratio and the localized lithium ion concentration. On the other hand, Gao and Zhou [18] applied the FEM and the J-integral method to study the softening effect during the lithiation process. They also captured the high lithium ion concentration at crack tip regions during the charging process, leading to relaxation of hydrostatic stress in later diffusion process. The fracture evolution in silicon nanowires also depends on the geometry of the structure. Ryu et al. [8] also concluded that silicon nanowires with diameters smaller than 300 nm will not fail during battery cycling, even pre-damaged nanowires. Due to their large pressure gradient and large volume expansion, lithium ions can rapidly diffuse into silicon nanowires. As an alternative approach, peridynamics [19][20][21][22][23][24][25] was used to investigate fracture evolution during the lithiation process by considering plane-stress conditions [26,27]. It was shown that peridynamics can capture accurate fracture propagation patterns for multiple cracks with random orientations inside the anode material.
In this study, the two-dimensional fully-coupled peridynamic model developed in Wang et al. [27] is extended to a three-dimensional peridynamic model. The effect of pressure gradient and material phase change during lithiation process (from Si to Li x Si) are taken into consideration. Coupled field equations are used to describe the lithiation process with the help of the newly developed peridynamic differential operator [26]. To demonstrate the capability of the model, the fracture evolution inside the cylindrical and spherical shaped lithium ion battery anodes is presented.

Coupled Diffusion-Mechanical Deformation Formulation
Fick's Second Law describes the basic principles of ion diffusion [28]. Due to silicon being the anode material, larger volume expansion occurs as lithium ions diffuse inside the anode. Hence, the diffusion-induced stress during battery cycling should be taken into account in Fick's Second Law. The mechanical strain component, ε M ij , can be expressed by considering the influence of lithium ion diffusion as [16,17]: where ε ij refers to the total strain, α represents the coefficient of expansion, C avg represents the current average concentration value of the interacting material points, and δ ij is the Kronecker Delta. For three-dimensional cases, the stress-strain constitutive equation can be written as: By substituting Equation (2) into Equation (1), the coupled field normal stresses can be expressed as: where C is current value of the normalized lithium ion concentration, and C max is the maximum value of the lithium ion concentration. The total normal strain in all three dimensions can be expressed in terms of displacements as: Local stress will rise along with lithium ion diffusion at regions with high geometrical singularity. A high pressure gradient will lead to a large amount of lithium ion diffusion into these regions, which will increase the lithium ion concentration. As a result, stress will release as volume expands. Hence, the general Fick's Second Law should be modified by considering the pressure-gradient as: where M is the molecular mobility, k B is Boltzmann constant, T is the absolute temperature, N A is Avogadro's constant and σ is the hydrostatic stress. Zhang et al. [29] emphasized that silicon has different lithiated stages during the charging process. In the early stage, since lithium ion concentration is low, silicon will transform into partial lithiated silicon (Li x Si). As the lithium ion concentration increases, the partially lithiated silicon will further transform into fully lithiated silicon (Li 15 Si 4 ). Moreover, since some material properties of lithiated silicon, such as the elastic modulus and fracture toughness, are lower than those of pure silicon, the lithiation process can also be regarded as a material softening process [30]. However, the size of the partially lithiated silicon region in a battery anode is very small compared to the anode's geometry [29]. Therefore, in this study, only fully lithiated silicon was considered.

Peridynamic Theory
An extensive number of studies about fracture mechanics based on classic continuum mechanics (CCM) can be found in the literature. The CCM applies the spatial partial differential equations to describe the motion of a material point. Within the framework of CCM, various tools, such as the FEM, boundary element method (BEM) and cohesive zone method (CZM) [31][32][33], have been applied in numerical fracture analysis. However, since partial differential equations require continuous material geometry, CCM, as a local theory, may have difficulty performing the numerical analysis of problems with discontinuous geometry, such as cracks and kinks. By regarding damage as newly formed boundaries, a remeshing process may be inevitable in numerical analysis [34].
Peridynamic theory (PD) was firstly introduced by Silling and Askari [35] as an alternative numerical approach to CCM. Unlike FEM, the PD uses spatial integral equations instead of partial differential equations to describe the motion of material points. On the other hand, PD is a non-local theory. Material points can build up interactions, called peridynamic bonds, with surrounding neighbors within a certain distance (δ), as shown in Figure 1. The collection of neighboring points is called the horizon (H x ). For neighboring points of material point x beyond the size of the horizon, it is assumed that the interactions are too weak and have little influence on point x, which can be ignored in a numerical simulation. By integrating the forces between material point x and its neighboring points, x', the equation of motion of material point x can be calculated and expressed as where b refers to the body force density and f is the pairwise force of each interaction (bond). Note that, in this study, the bond-based peridynamic theory is used for simplicity. However, extension of the formulation to the ordinary-state based formulation is straightforward [34,36].
In bond-based peridynamic theory, the bond force depends on the relative position and relative displacement of the material points associated with the bond which can be calculated as [37]: where the relative position of material points x and x ′ can be defined in the undeformed configuration as: Note that the bond force given in Equation (11) is calculated based on the deformed configuration. Therefore, the formulation has the capacity to capture large deformations. Moreover, the relative displacement can be expressed as: In Equation (11), α is the coefficient of expansion, s is the stretch of the bond and C avg represents the average lithium ion concentration of material points x and x ′ . µ is the failure parameter and c is the bond constant. For a three-dimensional isotropic material, the bond constant can be expressed as: Note that for the three-dimensional, bond-based peridynamic theory, the Poisson's ratio for a linear elastic isotropic material is equivalent to ν = 1/4.
In CCM, damage is specially treated by using remeshing or a pre-defined crack propagation path. On the other hand, the damage is represented as "bond breakage" in peridynamic theory. For brittle elastic material, once the bond stretch exceeds a critical stretch value, s c , after deformation, the bond will break and cannot recover by itself. The critical stretch depends on the critical strain energy release rate (G c ). For a three-dimensional model, the rate definition of the critical strain energy release in peridynamics has been derived by Silling and Askari [35] as: Hence, the critical stretch, s c , can be expressed as: The failure condition of each bond can be represented by a failure parameter as: Since the broken bond no longer carries any load, the total load associated with material point x will be redistributed to the remaining bonds. This leads to a change in mechanical conditions. The damage value of material point x can be defined as: which represents the percentage of broken bonds associated with material point x. The material point is totally damaged when the damage parameter is equivalent to 1, whereas it is undamaged when the damage parameter is equivalent to 0. Note that breakage of a single bond will not cause the emergence of a crack. Instead, all bonds passing through a crack surface should be broken for a crack to occur. In other words, as the number of broken bonds increases, cracks can emerge and eventually propagate inside the structure. Moreover, there is no need to define the crack path, as in some other existing approaches, such as the cracking particle method [38,39]. During the charging (lithiation) process, the anode will expand in volume, and high hydrostatic stresses will emerge inside the anode, especially at high singularity regions. For three-dimensional problems, the hydrostatic stress of a material point, x, can be calculated as [40]: The Cauchy stress tensor, σ ij , can be expressed in terms of the first Piola-Kirchoff stress, σ 0 , as [41]: where the deformation gradient, F, and its determinant, J, are defined as: and: In Equation (21), x ′ − x refers to the initial relative distance, and y ′ − y refers to the relative distance after deformation. The stress is not directly expressed by the bond-based peridynamic theory. However, the first Piola-Kirchoff stress can be calculated using peridynamic parameters as: Hence, the amount of hydrostatic stress depends on the displacement gradient and the lithium ion concentration, whereas the concentration change rate depends on the hydrostatic stress gradient and the lithium ion concentration gradient. Therefore, the coupled field equation is applied, as shown in Equation (9). In this study, the peridynamic form of partial differential equation given in Equation (9) was obtained by using peridynamic differential operator method developed by Madenci et al. [42].

Peridynamic Differential Operator Approach
The three-dimensional peridynamic differential operator can be derived by considering the Taylor series. Three-dimensional Taylor Series can be written by disregarding the higher order terms [41]: This is done by moving f (x) to the left-hand-side and multiplying each term in Equation (24) with a PD function, g p 1 p 2 p 3 3 (ξ), in which (p 1 , p 2 , p 3 = 0, 1, 2 and they cannot be equal to zero at the same time), and integrating throughout the horizon yields: The orthogonality property of the PD functions is invoked: and Equation (26) is substituted into Equation (25); then, the partial differential equations can be expressed in peridynamic form as: In Equation (27), the peridynamic functions, g p 1 p 2 p 3 3 (ξ), can be constructed by: ω 011 (|ξ|)ξ 2 ξ 3 (28) The weight function, ω q 1 q 2 q 3 (|ξ|) can be expressed as: where q 1 , q 2 , q 3 = 0, 1, 2 and they cannot be equal to zero at the same time.
According to Equation (28), the peridynamic functions, g p 1 p 2 p 3 3 (ξ), are composed of elements with an unknown coefficient matrix, a, weight function, ω and bond length in the x, y and z directions.
The unknown coefficient matrix, a, is derived from the peridynamic shape matrix, A, and the known coefficient matrix, b, via following relationship: The elements of the peridynamic shape matrix, A, can be expressed as: where the unit volume is dV x ′ = ξ 2 sin θdξdθdφ. Hence, the peridynamic shape matrix can be calculated as: The elements of the known matrix, b, can be expressed as: b p 1 p 2 p 3 n 1 n 2 n 3 = n 1 !n 2 !n 3 !δ n 1 p 1 δ n 2 p 2 δ n 3 p 3 Hence, the known coefficient matrix b can be written as: According to Equation (30), the unknown matrix a can be computed as: Substituting Equations (29) and (35) back to Equation (28), the peridynamic functions g p 1 p 2 p 3 3 (ξ) can be calculated as: where θ represents the angle between the projection of the bond over the x-y plane and x-axis. φ shows the angle between the bond and the z-axis, as shown in Figure 2 (red line shows the bond).  Substituting Equations (36)-(44) in Equation (27), the partial derivative expressions in peridynamic form can be obtained. Therefore, the extension of Fick's Second Law given in Equation (9) can be rewritten in peridynamic form as:

Numerical Studies
To demonstrate the capability of the current approach, spherical and cylindrical shaped anodes with different crack configurations were considered and the evolution of the damage during the lithiation process was investigated. In all cases, for the numerical solution, uniform spatial discretization was utilized with a grid size of 20 nm for the spherical model and 10 nm for the cylindrical model. Corresponding horizon sizes were specified as 60 nm and 30 nm, respectively. However, the recently introduced dual-horizon concept [20,22] can also be applied to improve the numerical efficiency of the solution. For the time integration, two different critical time step sizes for Equations (10) and (45) were determined as given in [34]. Since the critical time step size for each equation was different, we chose the smaller value, i.e., 2 × 10 −14 s, to ensure the stability of the solution.

Spherical Model
In the first case, a spherical energy storage particle with a pre-existing penny-shaped crack was considered. The materials of this particle, including pure silicon and lithiated silicon, were regarded as brittle materials. In addition, the concentration values were normalized by the maximum concentration as shown in Table 1. Before the charging started, the anode material remained as pure silicon. During the charging process, lithium ions at their maximum concentration were applied to the entire outer surface. Since the particle structure was free from any displacement constraints, it experienced free expansion during the charging process. For the three-dimensional bond-based peridynamic theory, the Poisson's ratio is limited to 1/4. However, the Poisson's ratio values of pure silicon and lithiated silicon, given in Table 1, are very close to the constrained value. Therefore, the limitation of bond-based peridynamic theory regarding Poisson's ratio would not be effective in this case. Other Poisson's ratio values can be specified using the ordinary-state based peridynamic formulation. Since the geometry of the spherical energy storage particle was symmetric and the penny-shaped crack was horizontally oriented, the sample planes were selected along the longitude and latitude of the sphere, shown in red and black, respectively, in Figure 3 to show the lithium ion concentration and mechanical deformation clearly inside the particle.  The penny-shaped crack was located in the central region of the anode structure. The diameter of the penny-shaped crack was half of the diameter of the anode and the crack was horizontally oriented, as shown in Figure 3. As lithiation progressed, the material at the anode surface region became lithiated silicon (Li 15 Si 4 ) first. Hence, the surface region experienced relatively large deformation while the central region remained undeformed. Due to this mechanical deformation state, compressive stresses emerged at the particle surface region, whereas tensile stresses formed at the particle central region, especially around the crack tip location. Hence, according to Equation (9), high hydrostatic stress affects the lithium ion distribution inside the spherical structure.
The results of lithium ion diffusion and diffusion-induced deformation during the charging process are shown in Figures 4 and 5. In order to have a clear understanding of the lithiation-induced damage in the three-dimensional structure, results are presented in both the x-z plane and x-y plane (crack surface plane). Due to the heterogeneous distribution of deformation in the spherical structure, high hydrostatic stresses arose at the edge of the penny-shaped crack where high geometrical singularity lies, as shown in Figures 4c and 5c. As the lithium ion diffused further into the spherical structure, the lithium ion concentration at the crack edge regions started to increase compared with the surrounding regions. However, since the penny-shaped crack had just started to propagate, as shown in Figures 4b and 5b, the change in lithium-ion concentration may not be very obvious from the concentration plot. The hydrostatic stress value around the crack edge region was around 1.5 GPa which led to a bond stretch far beyond the critical stretch of amorphous silicon. Hence, the penny-shaped crack would have continued to propagate until the stretch of the bonds at the crack tip region reduced below the critical stretch value given in Equation (16).

Cylindrical Model
In this section, several different pre-damaged models are considered to investigate the fracture evolution in the silicon nanowire. The silicon nanowire is represented by a cylindrical-shaped structure. The fracture analysis of cases with different crack orientations is separately discussed in the following sections.

Penny-Shaped Crack along the Horizontal Plane
A cylindrical-shaped structure with diameter L, shown in Figure 6, was considered to represent a silicon nanowire. A penny-shaped crack (marked as shadow lines in Figure 6), L/4 in diameter, was horizontally oriented at the center of the cylindrical structure. The material properties of the cylindrical anode were the same as for the spherical model in Table 1. Before charging the battery, the cylindrical material was pure silicon. During the charging process, lithium ions at maximum concentration was applied to the outer surface of the cylindrical structure. As a result, the material properties changed as the lithium ions diffused inside the electrode and the cylindrical structure deformed. Similar to the spherical model, the results are shown in two different plane views in order to present detailed fracture information. Due to the material phase change during the charging process, material points at the surface region expanded as pure silicon changed into lithiated silicon, whereas material points in the central region remained undeformed. Therefore, the hydrostatic stress in the outer surface region is shown as compression stress which is marked as blue in Figures 7c and 8c. The hydrostatic stress in the inner region of the cylindrical anode is shown as tensile stress, which is marked in red, especially at the crack tip region. In this case, the crack propagation is obvious. The penny-shaped crack edge lies close to the material phase boundary, as shown Figures 7b and 8b. Hence, according to Equation (9), the lithium ion concentration at the crack edge region is relatively higher than in the surrounding regions as shown in Figures 7d and 8d. Since the hydrostatic stress at the crack tip region is around 2.3 GPa, it leads to the bond stretching far beyond the critical stretch of the lithiated silicon. Therefore, the crack propagation will not stop until the stretch values of the bonds at the crack edge regions reduce below the critical stretch value.

Penny-Shaped Crack along the Vertical Plane
In contrast to the spherical energy storage particle model presented earlier, the cylindrical silicon nanowire model may show different fracture behaviors in the radial direction and longitudinal direction. Hence, in this case study, the penny-shaped crack along the longitudinal (z-direction) direction is under investigation. A penny-shaped crack with a diameter of L/4 was located at the central position in the cylinder model, as shown in Figure 9. However, in contrast to the case study in Section 5.2.1, the crack was vertically oriented (or along the x-z mid-plane). Before charging the lithium ion battery, the cylinder was composed of pure silicon and the lithium ion concentration was zero throughout the whole structure. During the charging process, the lithium ion concentration at maximum value was applied to the entire outer surface of the cylinder. As the lithium ion concentration increased, pure silicon transformed into fully lithiated silicon which led to volume expansion. The stress and damage induced by this volume expansion are shown in Figures 9 and 10.  The results are presented in the x-z mid-plane view ( Figure 10) and x-y mid-plane view ( Figure 11). During the charging process, the material points in the outer surface region were first lithiated and then expanded whereas the material points in the central region remained undeformed. As a result, a deformation gradient from outer surface to central cylinder was produced, and hydrostatic stress increased, especially at the crack edge regions, in accordance with Equations (19)- (23). Figures 10c  and 11c, show the location of compression stress at the outer surface region, which is represented in blue in the figure. Tension stress exists in the central cylinder region, especially at crack edge regions where geometrical singularity exists, which is represented as red in the figure. As the crack propagated, the crack opened which caused relatively high tension stress near the crack surface. Since the crack propagated close to the material phase boundary, the lithium ion concentration was relatively higher than surrounding regions in accordance with to Equation (9), as shown in Figures 10d and 11d.

Twin Penny-Shaped Cracks along Horizontal Planes
The case studies in Sections 5.2.1 and 5.2.2 focused on the fracture analysis of a single crack on the symmetry plane of the cylinder structure. Hence, due to symmetric geometry and loading, the penny-shaped crack propagated along the crack surface plane. However, in reality, multiple cracks can occur at any position inside the battery electrode. Hence, in this case study, a silicon cylinder with two cracks not on the symmetry plane was investigated. As shown in Figure 12, twin cracks were located at the central region with diameters equal to L/4 and the distance between these cracks was also L/4. Both of the cracks were horizontally oriented. Before charging the battery, the cylindrical nanowire was composed of pure silicon only. During the charging process, lithium ions at maximum concentration were applied on the entire outer surface of the cylinder. As a result, the material points on the surface became lithiated and expanded while the material points in the central region remained undeformed. The heterogeneous deformation led to an increase in hydrostatic stress in regions with geometrical singularity and sped up the diffusion process, as shown in Figure 13.  Compared with the case in Section 5.2.1, the crack propagation behavior was different. At the initial stage of the charging process, both cracks propagated along the crack surface plane. However, as the charging process continued, the upper crack propagated upward from its edge and the lower crack propagated downwards from its edge. These cracks repelled each other, as shown in Figure 13b. Due to the heterogeneous volume expansion inside the cylinder, a deformation gradient formed from surface of the cylinder to the central cylinder. As a result, the hydrostatic stress increased, especially at the crack edge regions, in accordance with Equations (19)- (23). Since material points at the surface region experienced expansion during the charging process, the compression stresses were evident, as represented in Figure 13c. The expansion of these material points caused tensile stress for the material points in the central region, especially around the crack edge region, as shown in Figure 13c. Since the crack edge is close to the material phase boundary, lithium ions would have diffused into regions around the crack edge with higher priority in accordance with Equation (9). Hence, the lithium ion concentration was relatively higher than surrounding regions, as shown in Figure 13d.

Twin Penny Shape Cracks along Vertical Planes
After the analysis of twin cracks along horizontal planes, a case with twin cracks along vertical planes was also considered. In this case study, twin vertical oriented cracks were arranged in the center region of a silicon cylinder, as shown in Figure 14. The diameters of both cracks were L/4 and the distance between these cracks was also L/4. Before charging the battery, the silicon cylinder was made of pure silicon. During the charging process, lithium ions at maximum concentration were applied on the entire outer surface of the cylinder. Hence, the material particles inside the outer surface region became lithiated and expanded. This influenced the lithium ion diffusion and hydrostatic stress inside the cylinder, as shown in Figure 15.  Figure 15a shows the initial damage in the view of the x-z mid-plane. Early in the charging process, penny shape cracks propagated along their crack surfaces. However, as the charging process continued, the left crack propagated towards the left and the right crack propagated towards the right. Generally, these cracks repelled each other, as shown in Figure 15b. Heterogeneous volume expansion led to the formation of a deformation gradient from the surface region to the central region. Therefore, hydrostatic stress rose inside the cylinder in accordance with Equations (19)- (23). The heterogeneous expansion during lithiation introduced compression stress on material points at the cylinder surface region, which is marked in blue in Figure 15c. As a consequence, the material points at the central cylinder region would have suffered from the tensile stress, marked in red. As the crack edge reached the material phase boundary, the lithium ion would have diffused into material points in the crack edge region first. Hence, a relatively high lithium ion concentration at the crack edge region as compared with its surrounding region was observed, as shown in Figure 15d.

Discussion
The results presented show that lithiation can influence the fracture behavior of a battery electrode. In the perspective of peridynamic theory, crack propagation is calculated based on the bond stretch which depends on the deformation of the structure. Hence, crack propagation generally depends on deformation which is induced by lithiation during battery charging. Since the pure silicon transforms into lithiated silicon, the material properties also change which leads to change in the critical bond stretch value. The comparison between fracture behavior with and without material phase change has been discussed in previous works [26,27]. In this article, since the material phase change during the charging process was also under consideration, the critical bond stretch was also influenced by the material properties.
For an electrode structure with symmetric geometry, the method of crack propagation is different depending on the position of the initial crack. If the initial cracks lie on the symmetry plane, cracks will propagate along the crack surface plane, since the geometry and loading are symmetric. However, if the initial cracks do not lie on the symmetry plane, the newly formed cracks may not align with the initial crack. For the twin cracks cases, the two-dimensional penny-shaped cracks turned into three-dimensional bowl-shaped cracks after the battery as charged. This is because of the material softening phenomenon during lithiation, as described in [32]. Due to the increase in lithium ion concentration, pure silicon at the crack edge regions will transform into lithiated silicon (Li 15 Si 4 ). As a consequence, the critical bond stretch is reduced. Hence, the bonds at crack edge region reach a critical stretch value earlier which leads the crack edge to propagate towards the high lithium ion concentration region. High hydrostatic stresses exist between the twin cracks, as shown in Figures 13c and 15c, and the bonds located in these regions may not reach the critical value. However, in situations of high hydrostatic stress, there is potential for these two cracks to merge into one larger crack as shown in Figure 15b.

Conclusions
In this study, peridynamic theory, in conjunction with the newly developed three-dimensional peridynamic differential operator approach, was utilized as a new method in three-dimensional fracture analysis during the lithiation process. Two differently shaped structures, a spherical energy storage particle and cylindrical nanowire, with different pre-existing cracks, were considered. According to numerical results, crack propagation is usually influenced by high hydrostatic stresses which yield bond breakage. In addition, the material phase change from pure silicon to lithiated silicon also influences the crack propagation. If the model geometry and loading are symmetric, the direction of crack propagation also depends on the initial crack position. High hydrostatic stresses lie between the cracks in twin crack cases, which means twin cracks have a potential to merge into one larger crack.
As a summary, it was shown that peridynamic theory is a suitable candidate for prediction of three-dimensional damage evolution in lithium ion battery electrodes. In peridynamic theory, the damage formation and evolution can be simplified without remeshing or applying sophisticated damage criteria. By using peridynamics, one can have better understanding of the failure mechanisms during the cyclic operation of lithium ion batteries.