Extended finite element method (XFEM)

For the purpose of fracture analysis, the enrichment functions typically consist of the near-tip asymptotic functions that capture the singularity around the crack tip and a discontinuous function that represents the jump in displacement across the crack surfaces. The approximation for a displacement vector function $\mathbf{u}$ with the partition of unity enrichment is

where $N_I\left(x\right)$ are the usual nodal shape functions. The first term on the right-hand side of the above equation, ${\mathbf{u}}_I$ , is the usual nodal displacement vector associated with the continuous part of the finite element solution. The second term is the product of the nodal enriched degree of freedom vector, ${\mathbf{a}}_I$ , and the associated discontinuous jump function $H\left(x\right)$ across the crack surfaces. The third term is the product of the nodal enriched degree of freedom vector, ${\mathbf{b}}_I^{\alpha }$ , and the associated elastic asymptotic crack-tip functions, $F_{\alpha }\left(x\right)$ . The first term on the right-hand side is applicable to all the nodes in the model; the second term is valid for nodes whose shape function support is cut by the crack interior; and the third term is used only for nodes whose shape function support is cut by the crack tip, as illustrated in Figure 1, where the circled nodes are enriched with the discontinuous jump function and the squared nodes with the asymptotic crack-tip functions.

Figure 1. Nodes in cracked elements enriched with discontinuous jump function or asymptotic crack-tip functions.

Figure 2. Normal and tangential coordinates for a smooth crack.

Figure 2 illustrates the discontinuous jump function across the crack surfaces, $H\left(x\right)$ , which is given by

where $\mathbf{x}$ is a sample (Gauss) point, ${\mathbf{x}}^{*}$ is the point on the crack closest to $\mathbf{x}$ , and $\mathbf{n}$ is the unit outward normal to the crack at ${\mathbf{x}}^{*}$ .

Figure 2 illustrates tangential and normal directions (with respect to the crack) at various points along the interior of the crack, as well as the crack tip. It also illustrates a local polar coordinate system at the crack tip, with respect to which the asymptotic crack tip functions in an isotropic elastic material $F_{\alpha }\left(x\right)$ are given by

where $\left(r,\theta \right)$ is the radial distance from the crack tip and $\theta =0$ is tangent to the crack at the tip.

As written, these functions span the asymptotic crack-tip function of elastostatics, and $\sqrt{r}\sin \frac{\theta }{2}$ takes into account the discontinuity across the crack face. Abaqus/Standard supports singular enrichment functions for isotropic linear elastic materials only. However, the use of asymptotic crack-tip functions is not restricted to crack modeling in an isotropic elastic material. The same approach can be used to represent a crack along a bimaterial interface, impinged on the bimaterial interface, or in an elastic-plastic power law hardening material. In each of these three cases different forms of asymptotic crack-tip functions are employed depending on the crack location and the extent of the inelastic material deformation, as discussed in Sukumar (2004), Sukumar and Prevost (2003), and Elguedj (2006).

Accurately modeling the crack-tip singularity requires constantly keeping track of where the crack propagates, which is cumbersome because the degree of crack singularity depends on the location of the crack in a nonisotropic material. Therefore, we consider the asymptotic singularity functions only when modeling stationary cracks in Abaqus/Standard. When modeling moving cracks, we always assume that the crack must propagate across an entire element at a time to avoid the need to model the near-tip asymptotic singularity. Therefore, only the displacement jump across a cracked element is considered when modeling moving cracks.

Phantom nodes, which are superposed on the original real nodes, are introduced to represent the discontinuity of the cracked elements, as illustrated in Figure 3. When the element is intact (not cracked yet), each phantom node is completely constrained to its corresponding real node. When the element is cut through by a crack, the cracked element splits into two parts. Each part is formed by a combination of some real and phantom nodes depending on the orientation of the crack. Each phantom node and its corresponding real node are no longer tied together and can move apart.

Figure 3. Principle of phantom node method.

The magnitude of the separation is governed by the cohesive law or the strain energy release rate until the traction on the cracked element is zero, after which the phantom and the real nodes move independently. To have a set of full interpolation bases, the part of the cracked element that belongs in the reference domain, ${\mathrm{\Omega }}_0$ , is extended to the phantom domain, ${\mathrm{\Omega }}_p$ . Then the displacement in the real domain, ${\mathrm{\Omega }}_0$ , can be interpolated by using the degrees of freedom for the nodes in the phantom domain, ${\mathrm{\Omega }}_p$ . The jump in the displacement field is realized by integrating the element quantities such as stiffness and stress only over the area from the side of the real nodes up to the crack; that is, ${\mathrm{\Omega }}_0^{+}$ and ${\mathrm{\Omega }}_0^{-}$ . This method provides an effective and attractive engineering approach and has been used for the simulation of the initiation and growth of multiple cracks in solids by Song (2006) and Remmers (2008). The equivalent polynomial methodology originally developed by Ventura and Benvenuti (2015) for a cracked element enriched with a Heaviside enrichment function has been extended in Abaqus/Standard to evaluate the stiffness matrix for the cracked element composed of real and phantom nodes. It has been proven to exhibit almost no mesh dependence if the mesh is sufficiently refined.

A key development that facilitates the treatment of cracks in an extended finite element analysis is the description of crack geometry because the mesh is not required to conform to the crack geometry. The level set method, which is a powerful numerical technique for analyzing and computing interface motion, fits naturally with the extended finite element method and makes it possible to model arbitrary crack growth without remeshing. The crack geometry is defined by two almost-orthogonal signed distance functions, as illustrated in Figure 4. The first, $\varphi$ , describes the crack surface, while the second, $\psi$ , is used to construct an orthogonal surface so that the intersection of the two surfaces gives the crack front. ${\mathbf{n}}^{+}$ indicates the positive normal to the crack surface, and ${\mathbf{m}}^{+}$ indicates the positive normal to the crack front. No explicit representation of the boundaries or interfaces is required because they are entirely described by the nodal data. Two signed distance functions per node are generally required to describe a crack geometry.

Figure 4. Representation of nonplanar crack in three dimensions by two signed distance functions $\varphi$ and $\psi$ .