## 3 Brittle Fracture Phase Field Model
### Governing Equations  
According to [PERSON] and [PERSON] [17], the fracture surface energy is obtained from the elastic energy, which determines the damage variable evolution. The elastic energy is decomposed into tensile and compressive components to ensure that only the tensile component of such density is governing the phase field variable evolution. This assumption cannot predict the crack initiation in the case of uniaxial specimen compression as a concrete compression test. It also discards the influence of the crack roughness [12] on the propagation phenomena. Here, the decomposition method described by [PERSON] et al. [27] was adopted and the decomposition of the strain tensor is as follows:  
\[ \begin{cases}\epsilon_{-}=\sum\limits_{a=1}^{d}\left\langle\epsilon_{a}\right\rangle _{-}\mathbf{n}_{a}\otimes\mathbf{n}_{a}\\ \epsilon_{+}=\sum\limits_{a=1}^{d}\left\langle\epsilon_{a}\right\rangle_{+} \mathbf{n}_{a}\otimes\mathbf{n}_{a}\end{cases}  \tag{5}\]  
where \(\epsilon_{-}\) and \(\epsilon_{+}\) are the compressive and tensile strain tensors, respectively. The operators \(\left\langle\cdot\right\rangle_{+}\) and \(\left\langle\cdot\right\rangle_{-}\) are the positive and the negative parts, respectively. \(\epsilon_{a}\) is the principal strain and \(\mathbf{n}_{a}\) is its direction vector. Using the spectral decomposition of the strain tensor, the elastic energy density is illustrated by the following equation:  
\[ \begin{cases}\psi_{0}^{+}=\frac{1}{2}\lambda\left\langle tr(\epsilon)\right\rangle _{+}^{2}+\mu\ tr(\epsilon_{+}^{2})\\ \psi_{0}^{-}=\frac{1}{2}\lambda\left\langle tr(\epsilon)\right\rangle_{-}^{2 }+\mu\ tr(\epsilon_{-}^{2})\end{cases}  \tag{6}\]  
To model the stiffness reduction due to damage phase propagation, the following equation was used [42]:  
\[\psi_{0}(\epsilon)=((1-\alpha)^{2}+k)\psi_{0}^{+}(\epsilon)+\psi_{0}^{-}( \epsilon) \tag{7}\]  
Figure 4: Body with: (**a**) sharp fracture; (**b**) diffuse fracture. where \(k\) is a numerical parameter that prevents numerical singularities if the phase field variable \(\alpha\) tends towards 1. Furthermore, this parameter has to be positive and very small compared to 1. In dynamic fracture propagation, the kinetic energy of the body is added to the energy balance. The kinetic energy of a solid is given by the following equation:  
\[E_{cin}(\dot{\mathbf{u}})=\frac{1}{2}\!\!\int\limits_{\Omega}\dot{\rho}\dot{\mathbf{u}} \dot{\mathbf{u}}\dot{\mathbf{u}}\dot{\mathbf{d}}\Omega \tag{8}\]  
where \(\dot{\mathbf{u}}=\frac{\partial\mathbf{u}}{\partial t}\) is the velocity vector and \(\rho\) is the mass density of the material. The total Lagrange energy functional is given by the following equation:  
\[ \begin{array}{l}L=-\frac{1}{2}\!\!\int\limits_{\Omega}\dot{\rho}\dot{\mathbf{u}} \dot{\mathbf{u}}\dot{\mathbf{d}}\Omega-\int\limits_{\Omega}\left[((1-\alpha)^{2}+k) \psi_{0}^{+}(\epsilon)+\psi_{0}^{-}(\epsilon)\right]\!\mathrm{d}\Omega-\int \limits_{\Omega}G_{c}\!\left[\frac{\alpha^{2}}{4\dot{\mathbf{u}}}+\eta|\
abla\alpha |^{2}\right]\!d\Omega\\ +\!\!\int\limits_{\Omega}\mathbf{b}.\mathbf{u}\ d\Omega+\int\limits_{\partial\Omega} \mathbf{f}\mathbf{u}\ d\bar{S}\end{array}  \tag{9}\]  
The functional \(L\)'s variation may be calculated, and its first variation must be 0, which leads us to the governing Equation (10) given by [PERSON] et al. [38]:  
\[\left\{ \begin{array}{l}\frac{\partial\sigma}{\dot{\alpha}\dot{\mathbf{x}}}+\mathbf{b} =\rho\ddot{\mathbf{u}}\\ -4\eta^{2}\Delta\alpha+\alpha=\frac{4\dot{\mathbf{u}}}{\dot{\mathbf{c}}_{c}}(1-\alpha) H\\ \forall\mathbf{x}:\ \psi_{0}^{+}<\psi_{0}^{-}\Rightarrow\alpha=0\end{array} \right. \tag{10}\]  
where \(\ddot{\mathbf{u}}=\frac{\partial^{2}\mathbf{u}}{\partial t^{2}}\) is the acceleration vector, \(\sigma\) the stress tensor and \(I\) is the unit tensor \(\mathrm{R}^{\mathrm{d}\times\mathrm{d}}\). The law of damaged material behavior given by the thermodynamic of continuums using Equations (6) and (7) is:  
\[