## 4 Numerical Examples
### Three-Point Bending Test  
A classical three-point bending test was employed to validate the PFM's capability in simulating fracture behavior under mechanically driven loading conditions. Figure 5 represents the geometry and boundary conditions of the symmetric three-point bending beam. The parameters of the material are as follows: \(G_{c}=42.6\) N/m, \(l_{0}=5\) mm, \(E=37\) GPa, and \(v=0.2\). The geometric model is discretized by quadrilateral elements, and the mesh is locally refined in the upper region of the notch. There are a total of 39,680 quadrilateral elements. Figure 6 displays the load-CMOD curve obtained from the three-point bending beam test, showing good agreement with the corresponding results reported by [PERSON] et al. [46]. Moreover, Figure 7 illustrates the phase-field simulated evolution of fracture propagation in the three-point bending specimen under mechanical loading. Figure 4: Comparison between numerical and analytical solutions. water or liquid nitrogen) to create steep thermal gradients at rock-fluid interfaces. This section employs a TM-PFM to investigate the underlying mechanisms. We chose a classic example, the 2D brittle solid quenching test, to simulate thermal shock. First, 99% Al\({}_{2}\)O\({}_{3}\) powder was melted into a 50 mm \(\times\) 100 mm \(\times\) 1 mm sheet. Then, the brittle solid sheet was heated to an initial temperature, which varied from 300 \({}^{\circ}\)C to 600 \({}^{\circ}\)C. Finally, the high-temperature brittle solid sheet was abruptly placed in a 20 \({}^{\circ}\)C water bath. Leveraging symmetry, we model the upper-left quadrant of the brittle plate (Figure 8). In the displacement field, the upper left is the free boundary condition, and the lower right is the displacement-constrained boundary condition. Moreover, in the temperature field, the upper left boundary is the convection boundary condition, and the lower right boundary is the thermal insulation boundary condition. The total calculation time is 100 ms, and the time step is 0.1 ms. The parameters of the material are as follows: \(\rho=9398\) kg/m\({}^{3}\), \(E=370\) GPa, \(v=0.3\), \(G_{c}=42.47\) N/m, \(l_{0}=0.05\) mm, \(\alpha_{\mathrm{T}}=7.5\times 10^{-6}\) /K, \(C_{p}=880\) J/(kg\(\cdot\) K), \(k=31\) W/(m\(\cdot\) K), and heat transfer coefficient \(h=5\times 10^{5}\) W/(m\({}^{2}\cdot\) K). To ensure grid independence, the element size was chosen as \(h_{\mathrm{size}}=l_{0}/4\). In the final post-processing of the numerical results, we perform two mirroring operations on the two axes of symmetry to obtain information on the entire brittle solid sheet. When subjected to thermal shock, high-temperature samples tend to develop a multitude of thermal cracks. Figure 9 shows that the numerical simulation results of the thermal crack pattern are in high agreement with the experimental findings, thereby substantiating the capacity of the proposed TM-PFM to accurately capture complex thermal fracture patterns. Furthermore, it is noticeable that the cracks exhibit a distinct hierarchical morphology. Based on the combined experimental and numerical findings, all cracks are categorized into three types according to their length: long, medium-long, and short. As depicted in Figure 10a, there is a gradual increase in the total number of cracks within the brittle solid sheet following thermal shock, as the initial temperature rises. Concurrently, the crack density becomes more concentrated. When comparing the numerical and experimental results, the counts of long and medium-long cracks align closely. However, the numerical results show a significantly higher number of short cracks than the experimental findings (see Figure 10b-d). This discrepancy is predominantly attributed to the difficulty in visually identifying short cracks in experiments, as they are not visible in experimental images. In contrast, numerical methods can overcome this limitation.