## 2 Methodology  
### Brief Description of RFPA3D  
Through continuum mechanics, RFPA3D is programmed to simulate discontinuous mechanical problems by reducing parameters after element failure. This code reproduces nonlinear deformation characteristics by introducing the heterogeneity of material properties into the model and assuming that the mechanical parameters conform to the distribution statistics (Weibull distribution function, distribution function, or normal distribution function). This paper assumes that the strength and elastic modulus follow the Weibull distribution function [14; 19; 29; 30], as:  
\[P(u)=\frac{m}{u_{0}}\bigg{(}\frac{u}{u_{0}}\bigg{)}^{m-1}\exp\bigg{(}\frac{u}{ u_{0}}\bigg{)}^{m}, \tag{1}\]  
where \(u_{0}\) is the mean value of the elements (such as strength or elastic modulus) put into the numerical simulation program, and \(u\) is the mechanical parameter of the element. In addition, the parameter m describes the shape of the distribution function and reflects the degree of homogeneity of the material; it is called the homogeneity index.  
It is initially linearly elastic until the specified damage threshold is reached. This study has calculated the damage thresholds using the maximum tensile strain criterion and the Mohr-Coulomb criterion. As expressed in Equation (2), the maximum tensile strain criterion primarily determines whether a crack has been initiated. The second damage threshold, defined in Equation (3), is the Mohr-Coulomb criterion, which determines whether an element is damaged in the shear mode in the absence of tensile damage.  
The maximum tensile strain criterion is:  
\[\varepsilon_{3}\leq\frac{\sigma_{t0}}{E}, \tag{2}\]  
The Mohr-Coulomb criterion is:  
\[F=\sigma_{1}-\sigma_{3}\frac{1+\sin\phi}{1-\sin\phi}\geq\sigma_{c0}, \tag{3}\]  
Until the element stress satisfies the above criteria, the elastic modulus remains constant and the same as before loading. According to the strain equivalence assumption, the elastic modulus of an element degrades gradually with damage progress, as follows:  
\[E=(1-D)E_{0}, \tag{4}\]  
It is characterized by an elastic-brittle damage constitutive relation, with a given residual strength (illustrated in Figure 1). As shown in the negative half of the axis of Figure 1, there are three kinds of elements determined by the damage variable, namely, elastic elements (\(D=0\)), damaged elements (\(0<D<1\)), and complete failure elements (\(D=1\)). Under multiaxial stress states, we can express the parameter \(D\) of an element as [31; 32; 33]:  
\[D=\left\{\begin{array}{rl}0&\overline{\varepsilon}>\varepsilon_{t0}\\ 1-\frac{\sigma_{tt}}{E_{0}}&\varepsilon_{tu}<\overline{\varepsilon}\leq \varepsilon_{t0}\\ 1&\overline{\varepsilon}\leq\varepsilon_{tu}\\ \end{array}\right., \tag{5}\]  
\[\varepsilon_{tu}=\eta\varepsilon_{t0}, \tag{6}\]  
where \(\overline{\varepsilon}\) is the equivalent principal strain. It is defined as follows:  
\[\overline{\varepsilon}=-\sqrt{\langle-\varepsilon_{1}\rangle^{2}+\langle- \varepsilon_{2}\rangle^{2}+\langle-\varepsilon_{3}\rangle^{2}}, \tag{7}\]where \(\varepsilon_{1}\), \(\varepsilon_{2}\), and \(\varepsilon_{3}\) are the three principal strains and \(\langle\rangle\) is a function defined as follows:  
\[\langle x\rangle=\left\{\begin{array}{cc}x&x\geq 0\\ 0&x<0,\end{array}\right. \tag{8}\]  
As shown in the positive half of the axis in Figure 1, for a shear mode element under multiaxial stress states, the damage variable \(D\) can be described as follows [31; 32; 33]:  
\[D=\left\{\begin{array}{cc}0&\varepsilon_{1}<\varepsilon_{0}\\ 1-\frac{\sigma_{tr}}{\varepsilon_{1}\varepsilon_{0}}&\varepsilon_{c0}\leq \varepsilon_{1},\end{array}\right. \tag{9}\]  
According to the constitutive law of elastic damage, the element gradually deteriorates. It is important to emphasize that the finite element analysis will fail if the modulus equals zero. Therefore, the limit of the elastic modulus is specified as a relatively small number, i.e., \(10^{-5}\).