## 3 Numerical Implementation
### Finite Element Discretization  
\[\mathbf{u}=\sum_{i}^{m}N_{i}^{u}\mathbf{u}_{i},\ T=\sum_{i}^{m}N_{i}^{T}T_{i},\ d= \sum_{i}^{m}N_{i}^{d}d_{i} \tag{26}\]  
where \(N_{i}^{u}=N_{i}^{T}=N_{i}^{d}=N_{i}\) is the shape function of node \(i\). \(\mathbf{u}_{i}\), \(T_{i}\), \(d_{i}\) are the values of the displacement field, temperature field, and the phase-field at node \(i\), respectively. The corresponding gradient is given by:  
\[\mathbf{\epsilon}=\sum_{i}^{m}\mathbf{B}_{i}^{u}\mathbf{u}_{i},\ \
abla T=\sum_{i}^{m}\mathbf{B}_{i}^{T}T_{i},\ \
abla d=\sum_{i}^{m}\mathbf{B}_{i}^{d}d_{i} \tag{27}\]where \(\mathbf{B}_{i}^{u}\), \(\mathbf{B}_{i}^{T}\), \(\mathbf{B}_{i}^{d}\) is the derivative of the corresponding shape function, respectively.  
Moreover, the test functions and their derivatives inherit identical discretization properties to the trial functions.  
\[\delta\mathbf{u}=\sum_{i}^{m}N_{i}^{u}\delta\mathbf{u}_{i},\ \delta T=\sum_{i}^{m}N_{i}^{T} \delta T_{i},\ \delta d=\sum_{i}^{m}N_{i}^{d}\delta d_{i} \tag{28}\]  
\[\delta\mathbf{\epsilon}=\sum_{i}^{m}\mathbf{B}_{i}^{u}\delta\mathbf{u}_{i},\ \
abla\delta T=\sum_{i}^{m}\mathbf{B}_{i}^{T} \delta T_{i},\ \
abla\delta d=\sum_{i}^{m}\mathbf{B}_{i}^{d}\delta d_{i} \tag{29}\]  
The weak form residuals of the coupled system are derived as:  
\[\left\{\begin{array}{l}\begin{aligned} \text{R}_{i}^{u}=& \underbrace{\int\Omega\mathbf{b}\cdot\delta\mathbf{u}\text{d}\,\Omega+\int \Omega_{i}\overline{\mathbf{t}}\cdot\delta\mathbf{u}\text{d}S}_{\mathbf{F}_{i}^{out}}- \underbrace{\int\Omega\mathbf{\sigma}\cdot\delta\mathbf{\epsilon}\text{d}\Omega}_{ \mathbf{F}_{i}^{out}}-\underbrace{\int\Omega\mathbf{\rho}\bar{u}\mathbf{\omega}\text{d} \,\Omega}_{\mathbf{F}_{i}^{in}\mathbf{\text{e}}^{in}}\\ \text{R}_{i}^{T}=&\underbrace{\int\Omega\mathcal{Q} \delta T\text{d}\,\Omega+\int\Omega_{i}\overline{J}\delta T\text{d}S}_{\mathbf{F} _{i}^{out}}-\underbrace{\int\Omega\
abla T\cdot\
abla\delta T\text{d}\,\Omega }_{\mathbf{F}_{i}^{T,int}}-\underbrace{\int\Omega\mathbf{\rho}C_{p}\bar{T}\delta T \text{d}\,\Omega}_{\mathbf{F}_{i}^{T,loc}}\\ \text{R}_{i}^{d}=&-\underbrace{\int\alpha\bigg{\{}\frac{G_{ c}}{l_{0}}\Big{(}d\delta d-l_{0}^{2}\
abla d\cdot\
abla\delta d\Big{)}-2(1-d)H( \mathbf{\epsilon}_{e})\delta d\bigg{\}}\,\text{d}\Omega}_{\mathbf{F}_{i}^{d,int}}\end{aligned} \end{aligned}\right. \tag{30}\]  
where \(R_{i}^{u}\), \(R_{i}^{T}\), \(R_{i}^{d}\) represents the residual terms corresponding to the displacement field, the temperature field, and the phase field, respectively. \(\mathbf{F}_{i}^{u,ext}\), \(\mathbf{F}_{i}^{u,int}\), \(\mathbf{F}_{i}^{u,line}\) are the external force terms. \(\mathbf{F}_{i}^{Text}\), \(\mathbf{F}_{i}^{T,int}\), \(\mathbf{F}_{i}^{u,loc}\) are the external force term, internal force term, and local internal energy increase term of the temperature field, respectively. \(\mathbf{F}_{i}^{d,int}\) is the internal force term of the phase field.  
\[\left\{\begin{array}{l}\begin{aligned} K_{ij}^{uu}= \frac{\partial F_{i}^{unit}}{\partial u_{i}}=&\int\Omega\big{(}\mathbf{B} _{i}^{u}\big{)}^{\mathsf{T}}\mathbb{C}_{\mathbf{c}}\mathbf{B}_{j}^{u}\text{d}\Omega\\ K_{ij}^{TT}=&\frac{\partial F_{i}^{T,int}}{\partial u_{i}}=& \int\Omega\big{(}\mathbf{B}_{i}^{T}\big{)}^{\mathsf{T}}k\mathbf{B}_{j}^{\mathsf{T}} \text{d}\Omega\text{d}\Omega\\ K_{ij}^{dd}=&\frac{\partial F_{i}^{d,int}}{\partial u_{j}}=& \int\Omega\Big{\{}N_{i}\Big{[}\frac{\delta_{c}}{l_{0}}+2H(\mathbf{\epsilon}_{e}) \Big{]}+\big{(}\mathbf{B}_{i}^{T}\big{)}^{\mathsf{T}}G_{c}l_{0}\mathbf{B}_{j}^{T} \Big{\}}\end{aligned}\right.\end{aligned}\right. \tag{31}\]  
where \(K_{ij}^{uu}\), \(K_{ij}^{TT}\), \(K_{ij}^{dd}\) represents the stiffness matrices of the displacement field, the temperature field, and the phase field, respectively.