## 3 Numerical Implementation of Phase Field Fracture Model Based on Finite Elements
### Numerical Implementation Process  
Since the total potential energy of the system exhibits non-convexity with respect to the coupled phase field and displacement fields but remains convex for individual field variables, a staggered solution algorithm is implemented to ensure numerical robustness [35; 36]. This decoupled approach is selected for its capability to address the non-convex energy landscape when solving displacement and phase field variables simultaneously, while leveraging convexity within each sub-problem. By alternately freezing one field and solving the other, the algorithm circumvents numerical instabilities inherent in fully coupled nonlinear systems. This strategy enhances computational efficiency by reducing equation system dimensions per iteration while maintaining solution accuracy. In practice, the governing equations are first converted to weak forms via the principle of virtual work. For instance, the displacement field \(u\) and phase field \(\phi\) are transformed into integral forms over the solution domain, which are subsequently discretized via finite element methods to derive linear equation systems. The discretized governing equation for the displacement field can be expressed as  
\[K^{u}u=f^{u} \tag{39}\]  
In this formula, \(K^{u}\) is the stiffness matrix of the displacement field and \(f^{u}\) is the load vector of the displacement field. The discretized governing equation for the phase field can be expressed as  
\[K^{\phi}\phi=f^{\phi} \tag{40}\]  
In this formula, \(K^{\phi}\) is the phase field stiffness matrix and \(f^{\phi}\) is the phase field load vector. Thus, the specific steps of the staggered solution algorithm are shown in Figure 2.  
At each time step, the elastic matrix is updated according to the current elastic strain and phase field distribution. Specifically, material damage severity is quantified via the phase field value, enabling the dynamic adjustment of elastic modulus and stiffness matrix parameters. Subsequently, the updated elastic matrix is utilized to recompute the stiffness matrix and load vector, propagating the solution iteratively until steady-state conditions are attained. To ensure numerical stability, convergence criteria enforce the following: residual force norm reduction below \(10^{-5}\times\) initial residual. Maximum absolute phase-field variation across nodes is less than \(10^{-6}\). The Newton-Raphson method addresses strong nonlinearities from plasticity-damage coupling, exploiting quadratic convergence near solutions to efficiently resolve challenges like crack propagation softening behavior. By linearizing residuals at each iteration, the method systematically minimizes errors between iterative and true solutions, effectively capturing abrupt stiffness changes during damage localization. For displacement field \(u\) and phase field \(\phi\), the authors of [24] establish a linearized system:  
\[ \begin{bmatrix}K^{\mu\mu}&0\\ 0&K^{\phi\phi}\end{bmatrix}\begin{Bmatrix}\delta u\\ \delta\phi\end{Bmatrix} =- \begin{Bmatrix}R^{u}\\ R^{\phi}\end{Bmatrix}  \tag{41}\]  
Figure 2: Flowchart of staggered solution algorithm. In this formula, \(K^{\mu\mu}\) and \(K^{\phi\phi}\) are the tangent stiffness matrices of the displacement field and the phase field, respectively, and they are the functions of the displacement field \(u\) and the phase field \(\phi\); \(R^{u}\) and \(R^{\phi}\) are the residual vectors, which reflect the difference between the current solution and the real solution. The basic principle of the Newton-Raphson iterative method is to linearly approximate the nonlinear function near the current solution and gradually approach the real solution by solving the system of linear equations. The specific iterative steps are shown in Figure 3.  
Anderson acceleration is integrated to expedite convergence during iterations. This method constructs optimized search directions through linear combinations of historical residuals, significantly accelerating convergence rates.