## 2 Phase-Field Fracture Model for Brittle Fracture
### Constitutive Equations of Cohesive Interface Elements  
For multiaxial loading or mixed-mode fractures, material heterogeneities can induce mixed-mode crack propagation. Consequently, corresponding TSLs are used in finite element simulations [53]. As can be seen in Figure 3, a bilinear TSL is characterised by penalty stiffness \(K\), threshold stress \(\sigma_{\text{max}}\), and critical fracture energy \(G_{\text{C}}\), which is adopted in this study. Fracture energy follows the form  
\[G_{\text{C}}=\int_{0}^{\delta_{f}}\sigma\,\text{d}\delta\quad\text{for}\quad 0 \leq\delta\leq\delta_{f}. \tag{13}\] \begin{table}
\begin{tabular}{l c c} \hline \hline
**Parameters** & **Standard Phase-Field Fracture** & **Cohesive Phase-Field Fracture** \\ \hline \(\phi(x)\) & \(\begin{cases}\phi(x)=\left(\frac{|x|}{2\ell_{c}}-1\right)^{2}&\text{AT1 model}&\phi(x)=1-\sin\left(\frac{|x|}{L}\right)\\ \phi(x)=\exp\left(-\frac{|x|}{L_{c}}\right)&\text{AT2 model}&\\ \alpha(\phi)=\phi\quad\text{AT1 model}&\alpha(\phi)=2\phi-\phi^{2}\\ \alpha(\phi)=\phi^{2}\quad\text{AT2 model}&\\ \omega(\phi)\quad\quad\quad\quad\quad\quad\quad\quad\omega(\phi)=(1-\phi)^{2 }+k&\omega(\phi)=\frac{(1-\phi)^{2}}{(1-\phi)^{2}+a_{1}\phi+a_{1}a_{2}\phi^{2 }+a_{1}a_{2}a_{3}\phi^{3}}\end{cases}\) \\ \(c_{0}\) & \(\begin{cases}c_{0}=\frac{8}{3}&\text{AT1 model}&c_{0}=\pi\\ c_{0}=2&\text{AT2 model}&\\ H(x,t)&H(\mathbf{x},t)=\max\{\max\psi_{\text{eq}}(\epsilon(x,t))\}&\begin{cases}H( \mathbf{x},t)=\max\left\{\psi_{\text{eq}}^{0}\,,\max\psi_{\text{eq}}(\epsilon (x,t))\right\}\\ \psi_{\text{eq}}^{0}=\frac{1}{2L}\sigma_{u}^{2}\end{cases}\) \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Comparison between standard and cohesive phase-field fracture models. Before the initiation of damage, the normal and shear stress can be calculated as  
\[ \begin{cases}\sigma_{n}=k_{n0}\langle\delta_{n}\rangle,&\text{(normal direction)}\\ \sigma_{t}=k_{s0}\,\delta_{s},&\text{(shear direction)},\end{cases}  \tag{14}\]  
where \(k_{n0}\) and \(k_{s0}\) are penalty stiffnesses in the normal and shear directions. The Macaulay bracket \(\langle\bullet\rangle\) indicates that compressive deformation or stress states do not initiate damage in the model. A quadratic stress damage criterion is used for damage initiation. Damage is presumed to commence when the function involving nominal stress ratios attains a value of one:  
\[\left(\frac{\langle\sigma_{n}\rangle}{\sigma_{n}^{0}}\right)^{2}+\left(\frac{ \sigma_{s}}{\sigma_{s}^{0}}\right)^{2}=1, \tag{15}\]  
where \(\sigma_{n}^{0}\) and \(\sigma_{s}^{0}\) are the tensile and shear cohesive strength, respectively. In FE simulations, when Equation (15) is satisfied and applied displacement increases, both normal and shear stress diminish, signifying the onset of the softening region. A scalar damage variable \(D\) is employed to quantify this stiffness degradation, encapsulating all damage mechanisms. \(D\) monotonically evolves from 0 to 1, reaching 1 upon complete failure. The normal and shear stress in the softening region are calculated via  
\[ \begin{cases}\sigma_{n}=\left(1-D\right)k_{n0}\langle\delta_{n}\rangle&\text{ (normal direction)},\\ \sigma_{t}=\left(1-D\right)k_{s0}\,\delta_{s}&\text{(shear direction)} \end{cases}  \tag{16}\]  
The evolution of the damage variable \(D\) is computed using the expression proposed by [54], i.e.,  
\[D=\frac{\delta_{m}^{f}(\delta_{m}^{\text{max}}-\sigma_{m}^{o})}{\delta_{m}^{ \text{max}}(\delta_{m}^{f}-\delta_{m}^{o})}. \tag{17}\]  
where \(\delta_{m}^{0}\) and \(\delta_{m}^{f}\) are the effective displacements at damage onset and complete failure, respectively, and \(\delta_{m}^{max}\) is the peak effective displacement during loading. The effective displacements in Equation (17) are computed using an expression suggested by [54] as  
\[\delta_{m}= \sqrt{\left\langle\delta_{n}\right\rangle^{2}+\delta_{s}^{2}}. \tag{18}\]  
Figure 3: The quadratic nominal stress criterion and mixed-mode fracture criterion [52].