## 1 Introduction  
The capability of industrial equipment to resist fractures during its whole service life is directly related to the solid materials used in industrial manufacturing. Numerical simulation is an economical and efficient way to study the relationship between the mechanical properties of such materials and fracture problems in engineering applications. With the demand for industrial equipment for wide applications, there is an increasing need for structural components operating in extremely high-temperature environments. However, the potential risk of accidents brings significant challenges in assessing the ability of structures such as engine nozzles, metallic pressure vessels and ceramic insulation coatings to withstand fractures under thermal shock. Therefore, it is of great importance to develop an accurate computational method to simulate the ductile or brittle fracture behaviors under thermomechanical deformation and predict the crack propagation arising from thermal stress. So far, some methods based on classical continuum mechanics have been proposed to investigate the thermomechanical behavior involving fractures, such as the boundary element method (BEM) [1; 2] and the finite element method (FEM) [3]. Incorporating damage mechanics, these methods have been utilized to simulate the multiple crack patterns of brittle solids in quenching tests. However, the implementation of the above methods is inevitably faced with some difficulties in describing the crack initiation and propagation, due to over-simplified crack modeling and mesh dependence. To overcome such limitations and model arbitrary thermomechanical crack patterns, a number of novel methods were developed, including the extended finite element method (XFEM) [4], the extended element free Galerkin method (XEFG) [5], the radial point interpolation method (RPIM) [6], and smoothed particle hydrodynamics (SPH) [7; 8]. These methods allow for better tracking of crack paths. Moreover, the cohesive zone model [9] and the phase-field (PF) model [10] have also been extended to simulate fracture problems with thermal effects [11; 12], further providing an effective approach for researchers. Although the aforementioned numerical methods and models can be applied for crack analysis under thermal shocks, some extra computational efforts such as element delete techniques [13] and damage equations [14] still introduce numerical difficulty due to the discontinuous crack surfaces. Peridynamics (PD), as an alternative non-local model in solid mechanics, was formulated in terms of space integration instead of differentiation in the governing equations by [PERSON] [15]. The main advantage is that the integration equation allows crack propagation at any paths without a special crack growth assumption. The original version of the PD model is called the bond-based PD (BB-PD) model [16; 17], where the assumption of equal pairwise force states results in a limitation on Poisson's ratio. To overcome this constraint, the state-based PD (SB-PD) model [18] was further proposed, including the ordinary state-based PD (OSB-PD) model [19; 20; 21] and the non-ordinary state-based PD (NOSB-PD) model [22; 23; 24]. As a more general form of the SB-PD model, the NOSB-PD model has widely attracted attention because of its advantage of easy application of classical continuum mechanics (CCM) to study complex material models, including plastic [25; 26], visco-plasticity [27], and hyperelastic [28]. For thermomechanical problems, some excellent works based on the PD model have been conducted. Early, the 1-D and 2-D PD formulations based on the BB-PD model were developed by [PERSON] and [PERSON] [29; 30]. [PERSON] et al. [31] extended the heat conduction equation with the SB-PD model based on the Lagrangian formalism. Subsequently, many studies of thermoelastic, thermoplastic, and thermodynamic problems [32; 33; 34; 35; 36] were investigated on the basis of the theoretical foundation of modeling thermal problems.