# Convergence Check Phase-Field Scheme for Modelling of Brittle and Ductile Fractures
## 3 Numerical Implementation into Finite Element Method  
Details on the numerical implementation of the phase-field method into the widely used finite element method can be found in [48]. Therefore, only basic relations will be presented herein. Numerical implementation is based on a virtual work principle \(\partial W^{\text{ext}}-\partial W^{\text{int}}=0\). According to Equation (15), variation in the internal energy potential is equal to:  
\[\delta W^{\text{int}}=\delta\Psi=\int_{\Omega}\mathbf{\sigma}\delta\varepsilon \text{d}\Omega+\int_{\Omega}\bigg{\{}\frac{\text{d}g(\phi)}{\text{d}\phi}\big{(} \psi_{\text{e}}^{+}+\psi_{\text{p}}\big{)}+\hat{F}(\overline{\psi})G_{\text{c} }\frac{\partial\gamma(\phi,\
abla\phi)}{\partial\phi}\bigg{\}}\delta\phi \text{d}\Omega. \tag{20}\]  
By discretisation of the virtual work, as explained in [61], we obtain:  
\[(\mathbf{F}_{\text{ext}}^{\mathbf{v}}-\mathbf{F}_{\text{int}}^{\mathbf{v}}) \delta\mathbf{v}+\Big{(}\mathbf{F}_{\text{ext}}^{\phi}-\mathbf{F}_{\text{int} }^{\phi}\Big{)}\delta\phi=0, \tag{21}\]  
where \(\mathbf{F}_{\text{ext}}^{\mathbf{v}}\) and \(\mathbf{F}_{\text{ext}}^{\phi}=0\) denote external force vectors with regard to the displacement field and phase field, respectively. Vector \(\mathbf{v}\) denotes the displacement vector. In the phase-field method, there are no external nodal forces, while in Equations (7) and (8), the usual external nodal force vector in mechanics is obtained:  
\[\mathbf{F}_{\text{ext}}^{\mathbf{v}}=\int_{\Omega}\mathbf{N}^{\mathbf{v}} \text{b}\text{d}\Omega-\int\limits_{\partial\Omega_{\text{b}}}\mathbf{N}^{ \mathbf{v}}\text{t}\text{d}\partial\Omega. \tag{22}\]  
In Equation (22), \(\mathbf{N}^{\mathbf{v}}\) is a matrix consisting of shape functions related to the displacement field. \(\mathbf{F}_{\text{int}}^{\mathbf{v}}\) and \(\mathbf{F}_{\text{int}}^{\phi}\) denote internal nodal force vectors regarding to the displacement and phase field, respectively [26]:  
\[\begin{split}\mathbf{F}_{\text{int}}^{\mathbf{v}}&= \int\limits_{\Omega}(\mathbf{B}^{\mathbf{v}})^{\mathsf{T}}\sigma\text{d}\Omega,\\ \mathbf{F}_{\text{int}}^{\phi}&=-\int\limits_{\Omega} \Big{\{}l^{2}\big{(}\mathbf{B}^{\phi}\big{)}^{\mathsf{T}}\mathbf{B}^{\phi}\phi +\big{[}(1+\mathcal{H})\mathbf{N}^{\phi}\phi-\mathcal{H}\big{]}\mathbf{N}^{ \phi}\Big{\}}\text{d}\Omega.\end{split} \tag{23}\]  
The matrix \(\mathbf{B}^{\mathbf{v}}\), appearing in the first relation of (23), is the usual strain matrix, consisting of the interpolation polynomial spatial derivatives. \(\mathbf{N}^{\phi}\) is the shape function matrix of the phase-field variable. \(\mathbf{B}^{\phi}\) is a matrix of spatial gradients of the matrix \(\mathbf{N}^{\phi}\), and \(\phi\) is a phase-field vector of the degrees of freedom (DOF). The corresponding stiffness matrices are:  
\[\begin{split}\mathbf{K}^{\mathbf{v}}&=-\frac{ \partial\mathbf{F}_{\text{int}}^{\mathbf{v}}}{\partial\mathbf{v}}=\int \limits_{\Omega}(\mathbf{B}^{\mathbf{v}})^{\mathsf{T}}\mathbf{C}\mathbf{B}^{ \mathbf{v}}\text{d}\Omega,\\ \mathbf{K}^{\phi\phi}&=-\frac{\partial\mathbf{F}_{ \text{int}}^{\text{ext}}}{\partial\phi}=\int\limits_{\Omega}\Big{\{}l^{2}\big{(} \mathbf{B}^{\phi}\big{)}^{\mathsf{T}}\mathbf{B}^{\phi}+[1+\mathcal{H}|\mathbf{ N}^{\phi}\mathbf{N}^{\phi}\Big{\}}\text{d}\Omega.\end{split} \tag{24}\]  
\(\mathbf{K}^{\mathbf{v}\mathbf{v}}\) and \(\mathbf{K}^{\phi\phi}\) represent stiffness matrices related to displacement and phase-field DOFs, while \(\mathbf{C}\) denotes the degraded material tangent matrix. The quantities derived constitute the basic matrices and vectors required for solving the nonlinear finite element equation.  
As mentioned, there are many approaches to solving the problem of the discretised couple of equations derived. In the next section, a novel staggered approach will be explained.