## 2 Problem and Equations
### The XFEM Spatial and Time Discretization  
Applying the trial and the test functions, the discretized forms of integral Equations (16) and (17) can be obtained according to the Bubnov-Galerkin technique as follows.  
\[\begin{array}{c}K\widehat{\overline{U}}-W\widehat{\overline{T}}-f^{ext}=0 \\ R\dot{\overline{T}}+C\dot{\overline{T}}-g^{ext}=0\end{array} \tag{27}\]  
where \(\dot{\overline{U}}^{T}=\left\langle\overline{u}^{T},\overline{a}^{T},\, \overline{b}^{T}\right\rangle\) and \(\dot{\overline{T}}^{T}=\left\langle\overline{T}^{T},\overline{e}^{T},\, \overline{f}^{T}\right\rangle\) are the complete set of the standard and enriched DOF of the displacement and temperature fields respectively. The matrices \(\mathbf{K}\), \(\mathbf{W}\), \(\mathbf{R}\), \(\mathbf{C}\) and external force vectors \(\mathbf{f}\) and \(\mathbf{g}\) are defined as follows.  
\[\begin{array}{c}\mathbf{K}_{\alpha\beta}=\int\limits_{\Omega}\left(B_{\mu}^ {\alpha}\right)^{T}D\beta_{\mu}^{\beta}d\Omega\\ \mathbf{W}_{\alpha x}=\int\limits_{\Omega}\left(B_{\mu}^{\alpha}\right)^{T}D \frac{1}{3}\beta mN_{p}^{\gamma}d\Omega\\ \mathbf{f}_{\alpha}^{ext}=\int\limits_{\Omega}\left(N_{\mu}^{\alpha}\right)^{ T}\varphi bd\Omega+\int\limits_{\Gamma_{t}}\left(N_{\mu}^{\alpha}\right)^{T} \overline{t}\,d\Gamma\end{array} \tag{28}\]and  
\[\begin{array}{l}\mathbf{R}_{\lambda\kappa}=\int\limits_{\Omega}\left(N_{T}^{ \Lambda}\right)^{T}(\rho c_{E})N_{T}^{\kappa}d\Omega\\ \mathbf{C}_{\lambda\kappa}=\int\limits_{\Omega}\left(\
abla N_{T}^{\Lambda} \right)^{T}\!k\
abla N_{T}^{\kappa}d\Omega\\ \mathbf{g}_{\lambda}^{ext}=\int\limits_{\Omega}\left(N_{T}^{\Lambda}\right)^{ T}\!q_{s}d\Omega+\int\limits_{\Gamma_{q}}\left(N_{T}^{\Lambda}\right)^{T}\!\overline{q}d\Gamma\end{array} \tag{29}\]  
where \((\alpha,\beta)\in(u^{std},a^{Hev},b^{tip})\) denote the \"standard\", \"Heaviside\", and \"asymptotic tip\" functions of the displacement field. \((\lambda,\kappa)\in(T^{std},e^{Hev},f^{tip})\) denote the \"standard\", \"Heaviside\", and \"asymptotic tip\" functions of the temperature field. In these definitions, \(\mathbf{m}\) is the vector of the delta Dirac function defined as \(\mathbf{m}=\left\{110\right\}^{T}\).  
In order to complete the numerical solution of XFEM equations, it is necessary to integrate the differential Equation (27) in time. The Newmark method is used for the elastic part of the problem [15; 17]. Meanwhile, in the numerical integration of the XFEM formulations, the crack discontinuity must be properly taken into account when the integration is performed, or this may lead to poor numerical results. So, the crack discontinuity must be taken as an internal boundary of the domain of the element, and the integration must be carried out by dividing the domain of the element into two sub-domains to perform the integration on these two sub-domains. The XFEM formulations involve the discontinuity along the crack and the singularity at the crack tip. In order to implement the standard integration method in the XFEM formulation, the element cut by the crack surface must be first divided into several sub-polygons, and the numerical integration is then performed over each sub-polygon. The numerical integration method of discontinuity is shown in Figure 4, in which the crack does not coincide with the rectangular edges that cause some approximation in the numerical simulation [15]. The discontinuity is divided into sub-triangles. It must be noted that more triangular sub-divisions are necessary in the tip element.