# Three Dimensional CS-FEM Phase-Field Modeling Technique for Brittle Fracture in Elastic Solids
## 2 Methods
### Three-Dimensional Cell Based Smoothed-Finite Elements
\begin{array}{ccc}\tilde{b}_{I1}(\mathbf{x}_{C})&0&0\\ 0&\tilde{b}_{I2}(\mathbf{x}_{C})&0\\ 0&0&\tilde{b}_{I3}(\mathbf{x}_{C})\\ \tilde{b}_{I2}(\mathbf{x}_{C})&\tilde{b}_{I1}(\mathbf{x}_{C})&0\\ 0&\tilde{b}_{I2}(\mathbf{x}_{C})&\tilde{b}_{I3}(\mathbf{x}_{C})\\ \tilde{b}_{I1}(\mathbf{x}_{C})&0&\tilde{b}_{I3}(\mathbf{x}_{C})\end{array} \right] \tag{23}\]  
where  
\[\tilde{b}_{I1}(\mathbf{x}_{C})=\frac{1}{V_{C}}\int_{\Gamma_{C}}N_{I}(\mathbf{x})n_{k}( \mathbf{x})d\mathbf{\Gamma} \tag{24}\]  
\(N_{I}\) is the regular shape functions for a eight-node hexahedral finite element. Since we are using eight-noded hexahedral elements, one Gauss point at the center is sufficient to integrate along each surface boundary. The above equation then reduces to  
\[\tilde{b}_{I1}(\mathbf{x}_{C})=\frac{1}{V_{C}}\sum_{I=1}^{M}N_{I}(\mathbf{x}_{i}^{GP}) n_{ki}^{C}A_{i}^{C} \tag{25}\]  
where \(x_{i}\) is the midpoint or the intersection of the diagonals of the boundary segment and \(\mathbf{n_{i}}\) and \(A_{i}\) are the corresponding outward unit normal and surface area, respectively. The smoothed stiffness matrix of an element is thus given by the assembly of the stiffness of each smoothing cell  
\[\tilde{K}_{e}=\sum_{C}\tilde{B}_{C}^{T}C\tilde{B}_{C}V_{C} \tag{26}\]  
The overall smoothed stiffness matrix of the domain is calculated by the summation of sub-stiffness matrices for nodes \(I\) in relation to node \(J\) as in FEM, except that the summation here is performed over the smoothing domains, not elements:  
\[\tilde{K}_{IJ}^{CS-FEM}=\sum_{i=1}^{N_{c}}\sum_{m=1}^{N_{c}}\int_{\Omega_{i,m} ^{S}}\tilde{B}_{I}^{T}C\tilde{B}_{J}d\mathbf{\Omega}=\sum_{i=1}^{N_{c}}\sum_{m=1} ^{N_{c}}\tilde{B}_{I}^{T}C\tilde{B}_{I}V_{i,m}^{s}=\sum_{k=1}^{N_{c}}\tilde{B }_{I}^{T}C\tilde{B}_{I}V_{k}^{s} \tag{27}\]  
Thus, it is an algebraic summation over the stiffness of each element and the equations required to calculate the stiffness matrix of an element are the Equations (23) and (24) where we do not need the spatial derivative of the shape functions, thereby eliminating the Jacobian inversion problem. Only the values of shape function at certain points (Gauss point/center of each surface of a smoothing cell) are needed (Figure 2). The calculation of displacements or force vector, based on the nature of the problem is similar to classical FEM, by solving the following algebraic equation  
\[\mathbf{\tilde{K}}^{CS-FEM}\mathbf{\tilde{d}}=\mathbf{\tilde{f}} \tag{28}\]  
Once the displacements are calculated, we deviate from the standard stress/strain calculations used in [1; 2; 38]. This is because ABAQUS does not let users obtain information from surrounding elements; the user element (UEL) is written such that it only specifies the formulation of a single element. In a typical CS-FEM setting [1; 2], the stress/strains at a node is calculated by the weighted average of the values of the surrounding smoothing domains of the node. In this case, due to our inability to access data from surrounding elements and because of the way the UEL is formulated, we calculate the elemental strains as  
\[\epsilon=\bar{B}\mathbf{u} \tag{29}\]  
where \(\bar{B}\) is the smoothed strain displacement matrix obtained at the center of the smoothing cell. The corresponding elemental stress is given by  
\[\sigma_{ij}=\mathbf{C}_{ijk}e_{kl} \tag{30}\]  
ABAQUS then allows its internal algorithms to calculate the nodal stress/strain values. Similarly, the elemental displacement values are given as  
\[\mathbf{u_{i}}=N_{ij}\mathbf{U_{j}} \tag{31}\]  
where \(U_{j}\) is the \(24\times 1\) nodal displacement vector and \(u_{i}\) is the \(3\times 1\) elemental displacement vector.