# Numerical Modeling of Temperature Effect on Tensile Strength of Granitic Rock
## 2 Numerical Methods
### Rock Fracture Model  
Considering a 2D body discretized with constant strain triangular (CST) finite elements having (possibly) intersecting cracks, i.e., displacement discontinuities, the discontinuities split some of the finite elements into two or more parts, as illustrated in Figure 1a. Two intersecting discontinuities are needed in the present application of modeling the effect of thermally induced cracks on the tensile strength of a rock specimen. More precisely, heat treatment of the numerical rock induces a substantial number of cracks with somewhat random orientations. As the fixed crack concept is adopted here (cracks do not rotate), an element with a crack unfavorably (close to parallel) oriented to the consequent uniaxial tension direction would generate spurious stresses without introducing the second crack. Under the small strain assumption, the displacement and strain fields for such an element can be written as:  
\[\mathbf{u}(\mathbf{x})=N_{i}(\mathbf{x})\mathbf{u}_{i}^{c}+\sum_{i=1}^{2}M_{\Gamma_{ \text{d}i}}(\mathbf{x})\mathbf{\alpha}_{\text{d}i}\ \ \text{ with }M_{\Gamma_{\text{d}i}}(\mathbf{x})=H_{\Gamma_{\text{d}i}}(\mathbf{x})- \varphi_{\Gamma_{\text{d}i}}(\mathbf{x}),\ \ (i=1,2) \tag{1}\]  
\[\epsilon(\mathbf{x})=(\
abla N_{i}\otimes\mathbf{u}_{i}^{c})^{sym}-\sum_{i=1}^{2} \left(\left(\
abla\varphi_{\Gamma_{\text{d}i}}(\mathbf{x})\otimes\mathbf{\alpha}_ {\text{d}i}\right)^{sym}+\delta_{\Gamma_{\text{d}i}}(\mathbf{n}_{i}\otimes \mathbf{\alpha}_{\text{d}i})^{sym}\right), \tag{2}\]  
where \(\mathbf{\alpha}_{\text{d}i}\) is the displacement jump (crack opening) for crack \(i\), and \(N_{i}\) and \(\mathbf{u}_{i}^{c}\) are the standard interpolation functions for the CST element and nodal displacements (\(i=1,\,\dots\,\), 3 with summation on repeated indices), respectively. Moreover, \(H_{\Gamma_{\text{d}i}}\) and \(\delta_{\Gamma_{\text{d}i}}\) denote, respectively, the Heaviside function and its gradient, the Dirac delta function. Consistent with the constant strain of the linear triangle element, the displacement jump is likewise an assumed elementwise constant, which means that \(\
abla\mathbf{\alpha}_{\text{d}i}\equiv\mathbf{0}\), and thus (2) follows in a straightforward manner by taking the gradient of (1). Moreover, the terms containing the Dirac's delta function, \(\delta_{\Gamma_{\text{d}i}}(\mathbf{n}_{i}\otimes\mathbf{\alpha}_{\text{d}i})^{sym}\), in (2), are non-zero only when \(\mathbf{x}\,\in\,\Gamma_{\text{d}i}\). Outside the discontinuity, this term is zero and can thus be neglected at the global level when solving the discretized equations of motion. Function \(M_{\Gamma_{\text{d}i}}\) in (1) restricts the effect of \(\mathbf{\alpha}_{\text{d}i}\) inside the corresponding finite element. This substantially facilitates the finite element implementation of the kinematics as there is no need for special treatment of the essential boundary conditions. The ramp function \(\varphi_{\Gamma_{\text{d}i}}\) appearing in \(M_{\Gamma_{\text{d}i}}\) is chosen, from among the combinations of the nodal interpolation functions so that its gradient is as parallel as possible to the crack normal \(\mathbf{n}_{i}\):  
\[\
abla\varphi_{\Gamma_{\text{d}i}}=\arg\left(\max_{k=1,2}\frac{\left|\sum_{n= 1}^{k}\
abla N_{n}\cdot\mathbf{n}_{i}\right|}{\left\|\sum_{i=1}^{k}\
abla N_{ n}\right\|}\right),\ (i=1,2) \tag{3}\]  
The displacement decomposition (1) and the related functions in the 1D case are illustrated in Figure 1b with a single two-node bar element under tension.