## 2 SPH Formulation  
In the SPH method, particle approximation is employed to discretize governing equations (mass, momentum and energy conservation equations) and yield SPH equations. With different discretization procedures, various SPH equations can be derived. The SPH equations for mass, momentum and energy conservation used in this work are given as follows:  
\[\frac{\mathrm{d}\rho_{i}}{\mathrm{d}t}=\rho_{i}\sum_{j=1}^{N}(v_{i}^{\alpha}- v_{j}^{\alpha})\frac{\partial W_{ij}}{\partial x_{i}^{\alpha}}\frac{m_{j}}{ \rho_{j}} \tag{1}\]  
\[\frac{\mathrm{d}v_{i}^{\alpha}}{\mathrm{d}t}=\sum_{j=1}^{N}m_{j}\left(\frac{ \sigma_{i}^{\alpha\beta}-Q_{ij}\delta^{\alpha\beta}}{\rho_{i}^{2}}+\frac{ \sigma_{j}^{\alpha\beta}-Q_{ji}\delta^{\alpha\beta}}{\rho_{j}^{2}}\right)\frac {\partial W_{ij}}{\partial x_{i}^{\beta}}+g_{i}^{\alpha} \tag{2}\]  
\[\frac{\mathrm{d}u_{i}}{\mathrm{d}t}=\frac{1}{2}\sum_{j=1}^{N}m_{j}\left(\frac{ \sigma_{i}^{\alpha\beta}-Q_{ij}\delta^{\alpha\beta}}{\rho_{i}^{2}}+\frac{ \sigma_{j}^{\alpha\beta}-Q_{ji}\delta^{\alpha\beta}}{\rho_{j}^{2}}\right)(v_{j} ^{\alpha}-v_{i}^{\alpha})\frac{\partial W_{ij}}{\partial x_{i}^{\beta}} \tag{3}\]  
where \(t\) is the time, \(\rho\) is the density, \(m\) is the mass, \(x\) is the coordinate, \(v\) is the velocity, \(\sigma\) is the stress, \(u\) is the internal energy per unit mass, \(g\) is the gravitational acceleration, the superscripts \(\alpha\) and \(\beta\) are tensor indices (ranging from 1 to 3 in 3D space), for which [PERSON] summation convention holds. The subscripts \(i\) and \(j\) denote particle number, \(\delta\) is the Kronecker delta function, \(W\) is the kernel function, and \(Q\) is the artificial viscosity. One can refer to the reference [19] for a detailed derivation of Equations (1)-(3) in the absence of artificial viscosity. To solve Equations (1)-(3) for the density, the velocity and the internal energy per unit mass of particles, the Leap-Frog algorithm is adopted to carry out the time integration of these equations, and a fixed time step is used. The time step is subject to the Courant-Friedrichs-Levy condition for the stability of time integration. It is estimatedas \(\kappa\cdot\min(h_{i}/c_{i})\), where \(h\) is the smoothing length, \(c\) is the sound speed, and \(\kappa\) is a safety coefficient, generally no more than 1.0. The coefficient \(\kappa\) is problem dependent. It can be taken as values close to 1.0 for some problems, such as the problems in Sections 4.1 and 4.2, when only considering the stability of time integration. However, conservative values are used in our practical calculations for better accuracy. The coefficient \(\kappa\) is in the range of 0.1 to 0.3 for the problems considered in this work. For the kernel function, the commonly used B-spline smoothing function, which was originally presented by [PERSON] and [PERSON] (2010), is used, and the symmetrization of the kernel function is employed. Then, the value of the kernel function is calculated for each interaction particle pair. It is calculated as  
\[W_{ij}=W(d_{ij},\overline{h}_{ij})=\alpha_{d}\times\left\{ \begin{array}{ll}2/ 3-\eta^{2}+\eta^{3}/2,&0\leq\eta<1\\ (2-\eta)^{3}/6,&1\leq\eta<2\\ 0,&\eta\geq 2\end{array} \right., \tag{4}\]  
where \(\eta=d_{ij}/\overline{h}_{ij}\), \(d_{ij}\) is the distance between particles \(i\) and \(j\), \(\overline{h}_{ij}=\left(h_{i}+h_{j}\right)/2\) is the average smoothing length, \(h_{i}\) and \(h_{j}\) are the smoothing lengths of particles \(i\) and \(j\), respectively, and \(\alpha_{d}\) is a constant, which is \(3/2\pi h^{3}\), for 3D problems considered in the present work. The smoothing length is an important parameter for SPH simulation. For large deformation problems, particle spacing generally changes significantly as deformation develops. If the smoothing length is kept constant, some particles may have few neighboring particles or even no neighboring particle in certain directions due to the large particle spacing formed in the deformation process. This can result in a numerical fracture.