## 4 Stress Intensity Factors (SIFs)  
The SIFs are extracted from the XFEM solution with the help of the crack tip path-independent integrals, which have been extended to the thermal crack problems [10; 11; 12; 15]. The interaction integral method is applied to calculate the SIFs [18]. The \(J\) integral over a path \(\Gamma_{J}\) around the crack tip is defined as (see Figure 5)  
\[J=\int_{\Gamma_{J}}\left[W\delta_{1i}-\sigma_{ij}\frac{\partial u_{j}}{ \partial x_{1}}\right]n_{i}d\Gamma \tag{30}\]  
Figure 4: The triangular sub-divisions scheme of the numerical integration of crack discontinuity in the extended finite element method (XFEM). \(W\) is the mechanical strain energy density defined as  
\[W=\frac{1}{2}\sigma:(\varepsilon-\varepsilon_{T}) \tag{31}\]  
\(\sigma_{ij}\) denotes the stress tensor, \(u_{j}\) is the displacement field, \(n_{i}\) is the unit outward normal vector to the contour integral \(\Gamma_{I}\), \(\delta\) is the Kronecker delta, and \(T_{i}=\sigma_{ij}n_{j}\) is the traction on the contour integral \(\Gamma_{J}\). In order to establish a domain form of the integral that is suitable for the extraction of the SIFs from a numerical solution, a second path \(\Gamma_{e}\) surrounding \(\Gamma_{J}\) is considered and the region \(A\) between \(\Gamma_{e}\) and \(\Gamma_{J}\) is defined. The region \(A\) is bounded by \(\Gamma_{J}+\Gamma_{+}+\Gamma_{e}+\Gamma_{-}\) where \(\Gamma_{+}\) and \(\Gamma_{-}\) are paths located on each of the crack faces. Equation (30) can be recast into  
\[J=\int_{\Gamma_{J}+\Gamma_{+}+\Gamma_{e}+\Gamma_{-}}\left[\sigma_{ij}\frac{ \partial u_{j}}{\partial x_{1}}-W\delta_{1i}\right]qmd\Gamma-\int_{\Gamma^{+} -\Gamma^{-}}\sigma_{2j}\frac{\partial u_{j}}{\partial x_{1}}qd\Gamma \tag{32}\]  
where \(m\) is the outward normal to \(A\). \(q\) is a weighting function defined over the domain of integration. The crack faces are assumed to be traction free. Applying the divergence theorem to Equation (32) and using the equilibrium and strain-displacement equations, the domain form of the \(J\) integral is then obtained by  
\[J=\int_{A}\left\{\left[\sigma_{ij}\frac{\partial u_{j}}{\partial x_{1}}-(W) \delta_{1i}\right]\frac{\partial q}{\partial x_{i}}+\left[\alpha\sigma_{kk} \frac{\partial T}{\partial x_{1}}\right]q\right\}dA \tag{33}\]  
where \(\alpha\) is the expansion coefficient. Two states of a cracked body are considered. The interaction integral is formulated by superimposing the actual and auxiliary fields on the path independent \(J\)-integral. State one corresponds to the actual fields. State two is an auxiliary field which will be chosen as the asymptotic field for modes I, II or III. The \(J\)-integral for the sum of the two states can be written as  
\[J^{(1,2)}=J^{(1)}+J^{(2)}+M^{(1,2)} \tag{34}\]  
where \(M^{(1,2)}\)is called the interaction integral between states one and two. The auxiliary fields are independent of the thermal loading, they are exactly the same as in non-thermoelasticity. The interaction integral is obtained as  
\[M^{(1,2)}=\int_{A}\left\{\left[\sigma_{ij}^{(1)}\frac{\partial u_{j}^{(2)}}{ \partial x_{1}}+\sigma_{ij}^{(2)}-W^{(1,2)}\delta_{1i}\right]\frac{\partial q }{\partial x_{i}}+\left[\alpha\sigma_{kk}^{(2)}\frac{\partial T}{\partial x_{ 1}}\right]q\right\}dA \tag{35}\]  
From the second term of the integrand of Equation (35), the thermal loading has a similar effect in the interaction integrals between the actual fields and some auxiliary fields, which is suitable for the extraction of the individual SIFs. Figure 5: \(J\) integral contour around the crack tip. For general mixed-mode loading, the part of the crack inside \(\Gamma_{J}\) is straight, it is independent of the path, and the value of the \(J\)-integral related to SIFs is obtained as  
\[J=\frac{1}{E^{*}}\Big{(}K_{\mathrm{I}}^{2}+K_{\mathrm{II}}^{2}\Big{)}+\frac{1}{2 \mu}K_{\mathrm{III}}^{2},\ E^{*}=\left\{ \begin{array}{ll}E&\text{plane stress}\\ E/\big{(}1-p^{2}\big{)}&\text{plane strain}\end{array} \right.