Modeling Discontinuities as an Enriched Feature Using the Extended Finite Element Method

Modeling discontinuities, such as cracks, as an enriched feature:

is commonly referred to as the extended finite element method (XFEM);
is an extension of the conventional finite element method based on the concept of partition of unity;
allows the presence of discontinuities in an element by enriching degrees of freedom with special displacement functions;
enables the modeling of discontinuities in the fluid pressure field as well as fluid flow within the cracked element surfaces as in hydraulically driven fracture;
can include the heat transport due to thermal conductance and radiation across the cracked element surfaces as well as thermal convection within the cracked element surfaces;
does not require the mesh to match the geometry of the discontinuities;
is a very attractive and effective way to simulate initiation and propagation of a discrete crack along an arbitrary, solution-dependent path without the requirement of remeshing in the bulk materials;
can be simultaneously used with the surface-based cohesive behavior approach (see Contact Cohesive Behavior) or the Virtual Crack Closure Technique (see Crack Propagation Analysis), which are best suited for modeling interfacial delamination;
can be performed using the static procedure (see Static Stress Analysis), the implicit dynamic procedure (see Implicit Dynamic Analysis Using Direct Integration), the general fatigue crack growth approach (see Linear Elastic Fatigue Crack Growth Analysis), the low-cycle fatigue analysis using the direct cyclic approach (see Low-Cycle Fatigue Analysis Using the Direct Cyclic Approach), the geostatic stress field procedure (see Geostatic Stress State), quasi-static analysis (see Quasi-Static Analysis), fully coupled thermal-stress analysis (see Fully Coupled Thermal-Stress Analysis), or coupled pore fluid diffusion/stress analysis (see Coupled Pore Fluid Diffusion and Stress Analysis);
can also be used to perform contour integral evaluations for an arbitrary stationary surface crack without the need to define the conforming mesh around the crack tip;
allows contact interaction of cracked element surfaces, including thermal interaction and pore fluid interaction, based on a small-sliding formulation or on a finite-sliding formulation within the general contact framework;
allows the application of distributed pressure loads or distributed heat fluxes to the cracked element surfaces;
allows the output of some surface variables on the cracked element surfaces;
allows both material and geometrical nonlinearity; and
is available only for first-order stress/displacement solid continuum elements, first-order displacement/pore pressure solid continuum elements, first-order displacement/temperature solid continuum elements, first-order displacement/pore pressure/temperature solid continuum elements, and second-order stress/displacement tetrahedral elements.

This page discusses:

Modeling Approach
Defining an Enriched Feature and Its Properties
Applying Cohesive Material Concepts to XFEM-Based Cohesive Behavior
Applying the VCCT Technique to the XFEM-Based LEFM Approach
Applying Distributed Pressure Loads and Heat Fluxes to Cracked Element Surfaces
Specifying the Initial Location of an Enriched Feature
Activating and Deactivating the Enriched Feature
Contour Integral
Procedures
Initial Conditions
Boundary Conditions
Loads
Predefined Fields
Material Options
Elements
Output
Limitations
Input File Template
References

Modeling Approach

Modeling stationary discontinuities, such as a crack, with the conventional finite element method requires that the mesh conforms to the geometric discontinuities. Creating a conforming mesh can be quite difficult. Modeling a growing crack is even more cumbersome because the mesh must be updated continuously to match the geometry of the discontinuity as the crack progresses.

The extended finite element method (XFEM) alleviates the need to create a conforming mesh. The extended finite element method was first introduced by Belytschko and Black (1999). It is an extension of the conventional finite element method based on the concept of partition of unity by Melenk and Babuska (1996), which allows local enrichment functions to be easily incorporated into a finite element approximation. The presence of discontinuities is ensured by the special enriched functions in conjunction with additional degrees of freedom. However, the finite element framework and its properties such as sparsity and symmetry are retained. XFEM does not alleviate the need for sufficient mesh refinement in the vicinity of the crack tip.

Introducing Nodal Enrichment Functions

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 $u$ with the partition of unity enrichment is

\begin{matrix} u \end{matrix} = \overset{N}{\sum_{I = 1}} N_{I} (x) [u_{I} + H (x) a_{I} + \overset{4}{\sum_{α = 1}} F_{α} (x) b_{I}^{α}],

where $N_{I} (x)$ are the usual nodal shape functions; the first term on the right-hand side of the above equation, $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, $a_{I}$ , and the associated discontinuous jump function $H (x)$ across the crack surfaces; and the third term is the product of the nodal enriched degree of freedom vector, $b_{I}^{α}$ , and the associated elastic asymptotic crack-tip functions, $F_{α} (x)$ . 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.

Figure 1 illustrates the discontinuous jump function across the crack surfaces, $H (x)$ , which is given by

H (x) = {\begin{matrix} 1 & if (x - x^{*}) . n \geq 0, \\ - 1 & otherwise, \end{matrix}

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

Figure 1 illustrates the asymptotic crack tip functions in an isotropic elastic material, $F_{α} (x)$ , which are given by

\begin{matrix} F_{α} (x) = [\sqrt{r} \sin \frac{θ}{2}, \sqrt{r} \cos \frac{θ}{2}, \sqrt{r} \sin θ \sin \frac{θ}{2}, \sqrt{r} \sin θ \cos \frac{θ}{2}], \end{matrix}

where $(r, θ)$ is a polar coordinate system with its origin at the crack tip and $θ = 0$ is tangent to the crack at the tip.

These functions span the asymptotic crack-tip function of elasto-statics, and $\sqrt{r} \sin \frac{θ}{2}$ takes into account the discontinuity across the crack face. 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. However, in each of these three cases different forms of asymptotic crack-tip functions are required depending on the crack location and the extent of the inelastic material deformation. The different forms for the asymptotic crack-tip functions are discussed by Sukumar (2004), Sukumar and Prevost (2003), and Elguedj (2006), respectively.

Accurately modeling the crack-tip singularity requires constantly keeping track of where the crack propagates and is cumbersome because the degree of crack singularity depends on the location of the crack in a non-isotropic material. Therefore, we consider the asymptotic singularity functions only when modeling stationary cracks in Abaqus/Standard. Moving cracks are modeled using one of the two alternative approaches described below.

Modeling Moving Cracks with the Cohesive Segments Method and Phantom Nodes

One alternative approach within the framework of XFEM is based on traction-separation cohesive behavior. This approach is used in Abaqus/Standard to simulate crack initiation and propagation. This is a very general interaction modeling capability, which can be used for modeling brittle or ductile fracture. The other crack initiation and propagation capabilities available in Abaqus/Standard are based on cohesive elements (Defining the Constitutive Response of Cohesive Elements Using a Traction-Separation Description) or on surface-based cohesive behavior (Contact Cohesive Behavior). Unlike these methods, which require that the cohesive surfaces align with element boundaries and the cracks propagate along a set of predefined paths, the XFEM-based cohesive segments method can be used to simulate crack initiation and propagation along an arbitrary, solution-dependent path in the bulk materials, since the crack propagation is not tied to the element boundaries in a mesh. In this case the near-tip asymptotic singularity is not needed, and only the displacement jump across a cracked element is considered. Therefore, the crack has to propagate across an entire element at a time to avoid the need to model the stress singularity.

Phantom nodes, which are superposed on the original real nodes, are introduced to represent the discontinuity of the cracked elements, as illustrated in Figure 2. When the element is intact, 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.

The magnitude of the separation is governed by the cohesive law until the cohesive strength of 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 real domain, $Ω_{0}$ , is extended to the phantom domain, $Ω_{p}$ . Then the displacement in the real domain, $Ω_{0}$ , can be interpolated by using the degrees of freedom for the nodes in the phantom domain, $Ω_{p}$ . The jump in the displacement field is realized by simply integrating only over the area from the side of the real nodes up to the crack; that is, $Ω_{0}^{+}$ and $Ω_{0}^{-}$ . This method provides an effective and attractive engineering approach and has been used for 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 to evaluate the stiffness matrix for the cracked element composed of real and phantom nodes in Abaqus/Standard. It has been proven to exhibit almost no mesh dependence if the mesh is sufficiently refined.

Modeling Hydraulically Driven Fracture

The cohesive segments method in conjunction with phantom nodes discussed above can also be extended to model hydraulically driven fracture. In this case additional phantom nodes with pore pressure degrees of freedom are introduced on the edges of each enriched element to model the fluid flow within the cracked element surfaces in conjunction with the phantom nodes that are superposed on the original real nodes to represent the discontinuities of displacement and fluid pressure in a cracked element. The phantom node at each element edge is not activated until the edge is intersected by a crack. The flow patterns of the pore fluid in the cracked elements are shown in Figure 3. The fluid is assumed to be incompressible. The fluid flow continuity, which accounts for both tangential and normal flow within and across the cracked element surfaces as well as the rate of opening of the cracked element surfaces, is maintained. The fluid pressure on the cracked element surfaces contributes to the traction-separation behavior of the cohesive segments in the enriched elements, which enables the modeling of hydraulically driven fracture.

Optionally, phantom nodes with both pore pressure and temperature degrees of freedom can be introduced on the edges of each enriched element, and phantom nodes with displacement, pore pressure, and temperature degrees of freedom can be superposed on the original real nodes. The hydraulically driven fracture functionality can be extended to include the heat transport due to thermal conductance and radiation across the cracked element surfaces as well as thermal convection within the cracked element surfaces.

Modeling Moving Cracks Based on the Principles of Linear Elastic Fracture Mechanics (LEFM) and Phantom Nodes

Another alternative approach to modeling moving cracks within the framework of XFEM is based on the principles of linear elastic fracture mechanics (LEFM). Therefore, it is more appropriate for problems in which brittle crack propagation occurs either upon reaching a critical value of a fracture criterion or associated with fatigue crack growth; for more discussion of fracture criteria and fatigue crack growth using the Paris law, see Fatigue Crack Growth Criterion, Linear Elastic Fatigue Crack Growth Analysis, and Low-Cycle Fatigue Analysis Using the Direct Cyclic Approach. The XFEM-based LEFM approach can be used to simulate crack propagation along an arbitrary, solution-dependent path in the bulk material.

Similar to the XFEM-based cohesive segments method described above, the near-tip asymptotic singularity is not considered, and only the displacement jump across a cracked element is considered. Therefore, the crack has to propagate across an entire element at a time to avoid the need to model the stress singularity. The strain energy release rate at the crack tip is calculated based on the modified Virtual Crack Closure Technique (VCCT), which has been used to model delamination along a known and partially bonded surface (see Crack Propagation Analysis).

Another similarity is the introduction of phantom nodes to represent the discontinuity of the cracked element when the fracture criterion is satisfied. The real node and the corresponding phantom node will separate when the equivalent strain energy release rate exceeds the critical strain energy release rate at the crack tip in an enriched element. The traction is initially carried as equal and opposite forces on the two surfaces of the cracked element. The traction is ramped down linearly over the separation between the two surfaces with the dissipated strain energy equal to either the critical strain energy required to initiate the separation or the critical strain energy required to propagate the crack depending on whether the VCCT or the enhanced VCCT criterion is specified.

Using the Level Set Method to Describe Discontinuous Geometry

A key development that facilitates 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, $ϕ$ , describes the crack surface, while the second, $ψ$ , is used to construct an orthogonal surface so that the intersection of the two surfaces gives the crack front. $n^{+}$ indicates the positive normal to the crack surface; $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.

Defining an Enriched Feature and Its Properties

An enriched feature is a region of a model, specified with an element set, in which the finite element interpolation functions are enhanced to include discontinuities. In other words, an enriched feature is a region in the model where you can use the extended finite element method (XFEM) to model a crack. You must specify an enriched feature and its properties to use XFEM. Because there is a computational cost associated with modeling discontinuities, limiting the extent of an enriched feature can enhance computational performance.

An enriched feature can model a stationary crack or a propagating crack, but not both. You must decide during model setup whether a particular enriched feature models a stationary crack or a propagating crack. For a stationary crack, the enrichment functions include the asymptotic elastic crack-tip fields (see Introducing Nodal Enrichment Functions). If you model a propagating crack, you must ensure that the enriched feature is large enough to include all areas of the part where the crack could potentially propagate. You should not include areas where the crack is unlikely to propagate in the enriched feature (if this information is known ahead of time). Including these areas increases the computational cost without having any impact on the accuracy of the model. Additional details regarding the criteria for initiation and propagation of a crack are discussed in Crack Initiation and Direction of Crack Extension.

The default XFEM implementation in Abaqus does not support interactions among multiple cracks or branching of a single crack. It is best suited for problems where only a single crack is modeled within an enriched feature. Real applications often require modeling multiple cracks within a single part. A part can contain no preexisting cracks or one or more preexisting cracks; it can also contain one or more cracks that initiate during the analysis. For such applications, you can use multiple enriched features within a part (one for each crack) or use a single enriched feature to model more than one crack, although the latter approach has limitations.

The choice of one approach versus the other depends on the ease of specifying multiple enrichment features at the time of model setup. This, in turn, depends on the proximity of the cracks to one another at initiation and during growth. You can choose a single enriched feature in situations where you do not have any prior knowledge regarding the potential locations of multiple cracks or when it is not straightforward to separate out two or more distinct enriched features in a complex part. The following paragraphs provide some guidelines to help you judge what types of situations you can model and to select the correct modeling approach for your problem, starting with the important concept of a level set conflict.

As discussed in Using the Level Set Method to Describe Discontinuous Geometry, an XFEM crack is represented in terms of nodal values of a level set function. A requirement of the level set representation is that the level set function at any given node must have a unique value. This requirement is violated when more than one crack cuts through a single element or adjacent elements of a given node, resulting in an attempt to assign multiple values of the level set function at the node. Such a condition is referred to as a level set conflict. If a level set conflict occurs during an analysis, the analysis ends with an error. Without any kind of prior experimental data, paths of propagating cracks are often difficult to predict. This can cause uncertainty prior to a simulation as to whether a level set conflict will occur if you model multiple cracks in a single enriched feature.

If there are multiple preexisting cracks in a part, you can define them within a single enriched feature as long as these cracks are relatively far away from one another to begin with and do not come close to one another during the crack growth process. This restriction ensures that a level set conflict does not occur at any node. Mesh refinement can help avoid a level set conflict.

By default, crack initiation exhibits the following behavior to avoid level set conflicts:

Crack initiation in an enriched feature cannot occur until all existing cracks propagate through the enriched feature boundary. As illustrated in Figure 5, Crack 2 does not initiate until Crack 1 propagates through the enriched feature boundary.
Crack initiation cannot occur near a preexisting crack (because this immediately causes a level set conflict).
More than one crack can initiate at nonpredetermined locations within an enriched feature only if two or more nonadjacent elements satisfy the damage initiation criterion in the same time increment and one of the following conditions is satisfied:
- There are no preexisting cracks in the enriched feature.
- All preexisting cracks in the enriched feature have propagated through the feature boundary.
A level set conflict occurs if the elements in which cracks initiate share nodes.
A newly initiated second crack cannot approach or enter an already cracked element.

You can use one of the following methods to change the default behavior described above:

Use multiple enriched features.
Use a single enriched feature that is capable of initiating more than one crack at different locations and at different time increments before existing cracks propagate through the feature boundary.

The sections that follow outline these two approaches.

Multiple Enriched Features

To use multiple enriched features, you must define the element sets belonging to each feature at the time of model setup. In other words, you must be able to determine the boundaries among multiple enriched features at the time of model setup. Figure 6 shows a part with two holes that result in stress concentrations. For the loading shown, two cracks would initiate and grow near these holes. For this part, it is easy to define multiple enriched features during model setup. If both cracks in Figure 7 are preexisting cracks and do not come close to each other as they grow, you can also use a single enriched feature. However, with a single enriched feature, performance suffers if most of the enriched feature remains uncracked. From a performance point of view, the better choice may still be to use two separate enriched features in the areas where the cracks grow.

Single Enriched Feature Capable of Initiating Multiple Cracks before Existing Cracks Propagate through the Feature Boundary

You may not be able to clearly define two or more enriched features at the time of model setup in all situations. Figure 7, which shows a plate with two closely spaced holes, illustrates an example of this situation. The potential crack initiation and growth paths are difficult to anticipate ahead of time. It may take some experience and modeling iterations to properly choose the boundaries that separate multiple enriched features. Using a single enriched feature may be more appropriate for this situation. XFEM cannot model cases in which cracks should coalesce during a simulation.

Another example where a clear separation of enriched features is difficult (maybe impossible) involves a composite panel with a hole and with +45/−45 plies. In this panel many parallel matrix cracks could nucleate and propagate in each ply, as illustrated in Figure 8. A single enriched feature per ply is more appropriate for this problem.

For the examples shown in Figure 7 and Figure 8, it is recommended that you use a single enriched feature (per ply), but activate the capability to initiate more than one crack at different locations and at different time increments even before preexisting cracks propagate through the feature boundary. This nondefault approach has fewer limitations compared to the default approach discussed earlier. During the model setup phase in this nondefault approach, you explicitly specify that multiple cracks can initiate within a single enriched feature. You also specify a small relative radius, $r_{s}$ , to define a crack initiation suppression zone near a crack tip within which the initiation of another crack is not allowed (to prevent a level set conflict). Multiple cracks are allowed to initiate in an enriched feature, but only outside the crack initiation suppression zone for each individual crack tip (see Activating and Deactivating the Enriched Feature). However, a level set conflict could still occur if two cracks propagate close to each other during the crack growth. If this happens, further mesh refinement in the region where the cracks approach each other may help run the analysis to completion.

Specifying the Enrichment

Abaqus activates the enriched degrees of freedom only when a crack cuts through an element. You can associate only stress/displacement, displacement/pore pressure, displacement/temperature, or displacement/pore pressure/temperature solid continuum elements with an enriched feature.

Defining the Type of Enrichment

You can choose to model an arbitrary stationary crack or a discrete crack propagation along an arbitrary, solution-dependent path. The former requires that the elements around the crack tips are enriched with asymptotic functions to catch the singularity and that the elements intersected by the crack interior are enriched with the jump function across the crack surfaces. The latter infers that crack propagation is modeled with either the cohesive segments method or the linear elastic fracture mechanics approach in conjunction with phantom nodes. However, the options are mutually exclusive and cannot be specified simultaneously in a model.

Assigning a Name to the Enriched Feature

You must assign a name to an enriched feature, such as a crack. This name can be used in defining the initial location of the crack surfaces, in identifying a crack for contour integral output, in activating or deactivating the crack propagation analysis, and in generating cracked element surfaces.

Identifying an Enriched Region

You must associate the enrichment definition with a region of your model. Only degrees of freedom in elements within these regions are potentially enriched with special functions. The region should consist of elements that are presently intersected by cracks and those that are likely to be intersected by cracks as the cracks propagate.

Defining a Crack Surface

As a crack propagates through the model, a crack surface representing both facets of cracked elements is generated on those enriched elements that are intersected by a crack during the analysis. You must associate the name of an enriched feature with the surface (see Assigning a Name to the Enriched Feature above).

Defining Contact of Cracked Element Surfaces Using a Small-Sliding Formulation

When an element is cut by a crack, the compressive behavior of the crack surfaces has to be considered. The formula that govern behavior are very similar to those used for surface-based small-sliding penalty contact (About Mechanical Contact Properties).

For an element intersected by a stationary crack or a moving crack with the linear elastic fracture mechanics approach, it is assumed that the elastic cohesive strength of the cracked element is zero. Therefore, compressive behavior of the crack surfaces is fully defined with the above options when the crack surfaces come into contact. For a moving crack with the cohesive segments method, the situation is more complex; traction-separation cohesive behavior as well as compressive behavior of the crack surfaces are involved in a cracked element. In the contact normal direction, the pressure-overclosure relationship governing the compressive behavior between the surfaces does not interact with the cohesive behavior, since they each describe the interaction between the surfaces in a different contact regime. The pressure-overclosure relationship governs the behavior only when the crack is “closed”; the cohesive behavior contributes to the contact normal stress only when the crack is “open” (that is, not in contact).

If the elastic cohesive stiffness of an element is undamaged in the shear direction, it is assumed that the cohesive behavior is active. Any tangential slip is assumed to be purely elastic in nature and is resisted by the elastic cohesive strength of the element, resulting in shear forces. If damage has been defined, the cohesive contribution to the shear stresses starts degrading with damage evolution. Once maximum degradation has been reached, the cohesive contribution to the shear stresses is zero. The friction model activates and begins contributing to the shear stresses.

Defining Contact of Cracked Element Surfaces Using General Contact

A general-purpose finite-sliding contact capability is available for XFEM-based crack surfaces within the general contact framework in Abaqus/Standard (About General Contact in Abaqus/Standard). Contact pairs do not support XFEM-based crack surfaces. Crack surfaces can be included in the general contact domain to participate in contact within the compressive regime. Contact between crack surfaces and contact between crack surfaces and other types of surfaces in the model can be considered. The general contact algorithm includes crack surfaces when the default all-inclusive exterior surface defines the contact domain. Alternatively, you can specify pairwise inclusions and exclusions to control more precisely the regions of the model that can potentially come into contact. In both cases the enrichment regions where cracks occur and participate in general contact have to be associated with at least one XFEM-based crack surface. In the former case the association of an enrichment region with an XFEM-based crack surface triggers general contact on cracks in that enrichment region. In the latter case the named XFEM-based crack surface can be explicitly used with contact inclusions and exclusions akin to other types of surfaces to precisely define the general contact domain (Defining the General Contact Domain).

General contact always takes precedence over the small-sliding contact formulation when both are included to model contact between opposing surfaces of a crack within the same analysis. The small-sliding formulation cannot be used to model contact between the crack surfaces and other surfaces; in that case, only general contact supports contact interactions. General contact is applicable only for a crack propagation analysis.

Thermal Contact Interaction at Cracked Element Surfaces

Optionally, you can specify thermal contact interaction properties (gap conductance, gap radiation, and gap heat generation) in a contact property definition for cracked element surfaces using the small-sliding formulation or using the general contact finite-sliding formulation. You can include all three types of thermal properties in the same contact property definition. General contact always takes precedence over the small-sliding contact formulation when both are included to model thermal contact between opposing surfaces of a crack within the same analysis.

Applying Cohesive Material Concepts to XFEM-Based Cohesive Behavior

The formulae and laws that govern the behavior of XFEM-based cohesive segments for a crack propagation analysis are very similar to those used for cohesive elements with traction-separation constitutive behavior (Defining the Constitutive Response of Cohesive Elements Using a Traction-Separation Description) and those used for surface-based cohesive behavior (Contact Cohesive Behavior). The similarities extend to the linear elastic traction-separation model, damage initiation criteria, and damage evolution laws.

Linear Elastic Traction-Separation Behavior

The available traction-separation model in Abaqus assumes initially linear elastic behavior followed by the initiation and evolution of damage. The elastic behavior is written in terms of an elastic constitutive matrix that relates the normal and shear stresses to the normal and shear separations of a cracked element.

The nominal traction stress vector, $t$ , consists of the following components: $t_{n}$ , $t_{s}$ , and (in three-dimensional problems) $t_{t}$ , which represent the normal and the two shear tractions, respectively. The corresponding separations are denoted by $δ_{n}$ , $δ_{s}$ , and $δ_{t}$ . The elastic behavior can then be written as

t = {\begin{matrix} t_{n} \\ t_{s} \\ t_{t} \end{matrix}} = [\begin{matrix} K_{n n} & 0 & 0 \\ 0 & K_{s s} & 0 \\ 0 & 0 & K_{t t} \end{matrix}] {\begin{matrix} δ_{n} \\ δ_{s} \\ δ_{t} \end{matrix}} = K δ .

The normal and tangential stiffness components will not be coupled: pure normal separation by itself does not give rise to cohesive forces in the shear directions, and pure shear slip with zero normal separation does not give rise to any cohesive forces in the normal direction.

The terms $K_{n n}$ , $K_{s s}$ , and $K_{t t}$ are calculated based on the elastic properties for an enriched element. Specifying the elastic properties of the material in an enriched region is sufficient to define both the elastic stiffness and the traction-separation behavior. For simplicity, we assume that $K_{n n} = K_{s s} = K_{t t}$ .

Damage Modeling

Damage modeling allows you to simulate the degradation and eventual failure of an enriched element. The failure mechanism consists of two ingredients: a damage initiation criterion and a damage evolution law. The initial response is assumed to be linear as discussed in the previous section. However, once a damage initiation criterion is met, damage can occur according to a user-defined damage evolution law. Figure 9 shows a typical linear and a typical nonlinear traction-separation response with a failure mechanism. The enriched elements do not undergo damage under pure compression.

Damage of the traction-separation response for cohesive behavior in an enriched element is defined within the same general framework used for conventional materials (see About Progressive Damage and Failure). However, unlike cohesive elements with traction-separation behavior, you do not have to specify the undamaged traction-separation behavior in an enriched element.

Crack Initiation and Direction of Crack Extension

Crack initiation refers to the beginning of degradation of the cohesive response at an enriched element. The process of degradation begins when the stresses or the strains satisfy specified crack initiation criteria. Crack initiation criteria are available based on the following Abaqus/Standard built-in models:

the maximum principal stress criterion,
the maximum principal strain criterion,
the maximum nominal stress criterion,
the maximum nominal strain criterion,
the quadratic traction-interaction criterion,
the quadratic separation-interaction criterion, and
the three-dimensional LaRC05 criterion.

In addition, a user-defined damage initiation criterion can be specified in user subroutine UDMGINI.

An additional crack is introduced or the crack length of an existing crack is extended after an equilibrium increment when the fracture criterion, f, reaches the value 1.0 within a given tolerance:

1.0 \leq f \leq 1.0 + f_{t o l} .

You can specify the tolerance $f_{t o l}$ . If $f > 1 + f_{t o l}$ , the time increment is cut back such that the crack initiation criterion is satisfied. The default value of $f_{t o l}$ is 0.05. To improve performance, a separate tolerance $f_{t o l}^{g}$ can be specified to control the crack growth of an existing crack while $f_{t o l}$ is used to control the nucleation of an additional crack. If it is not specified, the growth tolerance, $f_{t o l}^{g},$ is set equal to $f_{t o l} .$

Fracture of Multiple Elements in an Unstable Crack Growth Analysis

For an unstable crack growth problem, sometimes it is more efficient to allow multiple elements at and ahead of a crack tip to fracture without excessively cutting back the increment size when the fracture criterion is satisfied. This capability is activated automatically if you specify an unstable growth tolerance, $f_{t o l}^{u}$ . In this case if the fracture criterion, f, is within the given unstable growth tolerance:

1.0 + f_{t o l} \leq f \leq 1.0 + f_{t o l}^{u},

where

f_{t o l}

is the tolerance described earlier in this section, the time increment size by default is immediately reduced automatically to a very small value,

β ▵ t_{\min} .

Reducing the time increment size allows more elements to fracture until

f < 1

for all the elements ahead of the crack tip. You can, however, optionally specify the maximum number of cutbacks allowed,

n

, to be controlled by the regular tolerance,

f_{t o l}

, prior to the activation of the unstable growth tolerance in an increment. After this limit the time increment size is recovered automatically to a larger value,

α ▵ t_{p r e}

, where

▵ t_{\min}

is the minimum time increment allowed,

▵ t_{p r e}

is the time increment size prior to the unstable crack growth, and

α

β

and

n

are scaling parameters. The default values of

α

β

and

n

are 0.5, 2.0, and 0, respectively. If you do not specify a value for the unstable growth tolerance, the default value is infinity. In this case the fracture criterion, f, for unstable crack growth is not limited by any upper bound value in the above equation.

Specifying the Crack Direction

When the maximum principal stress or the maximum principal strain criterion is specified, the newly introduced crack is always orthogonal to the maximum principal stress/strain direction when the fracture criterion is satisfied. However, when one of the other Abaqus/Standard built-in crack initiation criteria is used, you have to specify if the newly introduced crack will be orthogonal to the element local 1-direction or orthogonal to the element local 2-direction (see Conventions) when the fracture criterion is satisfied. By default, the crack is orthogonal to the element local 1-direction. If a user-defined damage initiation criterion is specified, the normal direction to the crack plane or the crack line can be defined in user subroutine UDMGINI.

Maximum Principal Stress Criterion

The maximum principal stress criterion can be represented as

f = {\frac{⟨ σ_{m a x} ⟩}{σ_{m a x}^{o}}} .

Here, $σ_{m a x}^{o}$ represents the maximum allowable principal stress. The symbol $⟨ ⟩$ represents the Macaulay bracket with the usual interpretation (that is, $⟨ σ_{m a x} ⟩ = 0$ if $σ_{m a x} < 0$ and $⟨ σ_{m a x} ⟩ = σ_{m a x}$ if $σ_{m a x} \geq 0$ ). The Macaulay brackets are used to signify that a purely compressive stress state does not initiate damage. Damage is assumed to initiate when the maximum principal stress ratio (as defined in the expression above) reaches a value of one.

Maximum Principal Strain Criterion

The maximum principal strain criterion can be represented as

f = {\frac{⟨ ε_{m a x} ⟩}{ε_{m a x}^{o}}} .

Here, $ε_{m a x}^{o}$ represents the maximum allowable principal strain, and the Macaulay brackets signify that a purely compressive strain does not initiate damage. Damage is assumed to initiate when the maximum principal strain ratio (as defined in the expression above) reaches a value of one.

Maximum Nominal Stress Criterion

The maximum nominal stress criterion can be represented as

f = \max {\frac{⟨ t_{n} ⟩}{t_{n}^{o}}, \frac{t_{s}}{t_{s}^{o}}, \frac{t_{t}}{t_{t}^{o}}} .

The nominal traction stress vector, $t$ , consists of three components (two in two-dimensional problems). $t_{n}$ is the component normal to the likely cracked surface, and $t_{s}$ and $t_{t}$ are the two shear components on the likely cracked surface. Depending on what you specify (see Specifying the Crack Direction above), the likely cracked surface will be orthogonal either to the element local 1-direction or to the element local 2-direction. Here, $t_{n}^{o}$ , $t_{s}^{o}$ , and $t_{t}^{o}$ represent the peak values of the nominal stress. The symbol $⟨ ⟩$ represents the Macaulay bracket with the usual interpretation. The Macaulay brackets are used to signify that a purely compressive stress state does not initiate damage. Damage is assumed to initiate when the maximum nominal stress ratio (as defined in the expression above) reaches a value of one.

Maximum Nominal Strain Criterion

The maximum nominal strain criterion can be represented as

f = \max {\frac{⟨ ε_{n} ⟩}{ε_{n}^{o}}, \frac{ε_{s}}{ε_{s}^{o}}, \frac{ε_{t}}{ε_{t}^{o}}} .

Damage is assumed to initiate when the maximum nominal strain ratio (as defined in the expression above) reaches a value of one.

Quadratic Nominal Stress Criterion

The quadratic nominal stress criterion can be represented as

f = {\frac{⟨ t_{n} ⟩}{t_{n}^{o}}}^{2} + {\frac{t_{s}}{t_{s}^{o}}}^{2} + {\frac{t_{t}}{t_{t}^{o}}}^{2} .

Damage is assumed to initiate when the quadratic interaction function involving the stress ratios (as defined in the expression above) reaches a value of one.

Quadratic Nominal Strain Criterion

The quadratic nominal strain criterion can be represented as

f = {\frac{⟨ ε_{n} ⟩}{ε_{n}^{o}}}^{2} + {\frac{ε_{s}}{ε_{s}^{o}}}^{2} + {\frac{ε_{t}}{ε_{t}^{o}}}^{2} .

Damage is assumed to initiate when the quadratic interaction function involving the nominal strain ratios (as defined in the expression above) reaches a value of one.

Larc05 Three-Dimensional Criterion

The LaRC05 three-dimensional criterion can be applied generally to polymer-matrix fiber-reinforced composites. This criterion considers four different damage initiation mechanisms: matrix cracking, fiber kinking, fiber splitting, and fiber tension. For detailed information on the damage initiation criterion, see Larc05 Criterion.

The initiation criterion that first reaches a value of 1.0 determines the damage initiation.

User-Defined Damage Initiation Criterion

User subroutine UDMGINI provides a general capability for implementing a user-defined damage initiation criterion.

You can define several damage initiation mechanisms in user subroutine UDMGINI. You represent each damage initiation mechanism by a fracture criterion, $f_{i n d e x i}$ , and its associated normal direction to the crack plane or the crack line. Although you can define several damage initiation mechanisms, the actual damage initiation for an enriched element is governed by the most severe damage initiation mechanism:

f = \max {f_{i n d e x 1}, f_{i n d e x 2}, \dots f_{i n d e x n}} .

Damage is assumed to initiate when f, as defined in the expression above, reaches a value of one.

You must specify any material constants that are needed in user subroutine UDMGINI as part of a user-defined damage initiation criterion definition.

Limiting the Crack Propagation Direction

When the maximum principal stress, maximum principal strain, or user-defined damage initiation criterion is specified, you can limit the new crack propagation direction to within a certain angle (in degrees) of the previous crack propagation direction. The default is 85°.

Local Calculations of the Stress and Strain Fields Ahead of the Crack Tip

An accurate and efficient evaluation of the stress/strain fields ahead of the crack tip is important for both evaluating the crack initiation criterion and computing the crack propagation direction when needed. Abaqus/Standard offers several options for computing these fields.

Centroidal Values of Stress and Strain

By default, the stress/strain computed at the element centroid ahead of the crack tip is used to determine if the damage initiation criterion is satisfied and to determine the crack propagation direction. See Figure 10.

Computing the Stress and Strain Fields at the Crack Tip

With a sufficiently refined mesh, the centroidal approximation is accurate and economical. However, if the finite element mesh in the vicinity of the crack tip is coarse relative to the gradients in the stress/strain fields, the default centroidal approximation may not be sufficient. In such cases you can use the stress/strain extrapolated to the crack tip to determine if the damage initiation criterion is satisfied and to determine the crack propagation direction. See Figure 10.

Combining Crack Tip and Centroidal Calculations

You can also choose to combine the two previous alternatives: you can use the stress/strain values extrapolated to the crack tip to determine if the damage initiation criterion is satisfied, and you can use the stress/strain values at the element centroid to determine the crack propagation direction.

Nonlocal Averaging of the Stress/Strain Fields and Smoothing of the Crack Surface Normals to Improve the Accuracy of Crack Propagation Directions

The three options for evaluating the stress and strain fields discussed above are local calculations in the sense that the evaluated fields are local to the single element ahead of the crack tip. In the case of coarse and/or unstructured meshes a nonlocal averaging of the stress and strain fields ahead of the crack tip can lead to a more accurate evaluation of those fields, which can improve the accuracy of the computed propagation directions. In addition, a moving least-squares approximation by polynomials is used by default to obtain more accurate crack propagation directions. The least-squares approximation further smooths out the normals of the individual crack facets in elements along the crack front that satisfy the damage initiation criterion.

Specifying the Region of the Model Used for Nonlocal Averaging and Smoothing

To control the range of elements used for nonlocal averaging and smoothing in the crack direction calculations, you can specify a radius, $r_{c}$ , within which the elements ahead of the crack tip are included (see Figure 11). The default radius is three times the typical element characteristic length in the enriched region.

Smoothing the Stress/Strain Fields before Averaging

To further improve the nonlocal averaging, you can request an initial smoothing of the stress/strain fields ahead of the crack. In this case Abaqus/Standard averages the field values to element nodes and then interpolates the smoothed fields to the integration points. Once smoothing is complete, the nonlocal averaging is applied. No smoothing is applied by default.

Weighting Schemes for Nonlocal Averaging

Abaqus/Standard offers a number of weighting schemes for field smoothing that provide additional control over nonlocal averaging. For example, you may want to give a higher weighting to elements close to the crack tip. You can specify a weight function, $ω$ , to compute the average stress/strain based on the distance from the element integration points to the crack tip, $r$ . By default, a uniform weighting is applied to all elements used for averaging; alternatively, you can use a Gaussian function or a cubic spline function. You can also define a weight function with a user subroutine.

The Gaussian function is represented by:

ω (r) = \frac{1}{{(2 π)}^{3 / 2} r_{c}^{3}} e x p (\frac{- r^{2}}{2 r_{c}^{2}})

The cubic spline function is represented by:

ω (r) = {\begin{matrix} 4 (\frac{r}{r_{c}} - 1) {(\frac{r}{r_{c}})}^{2} + \frac{2}{3}, & 0 < r < \frac{r_{c}}{2} \\ \frac{4}{3} {(1 - - \frac{r}{r_{c}})}^{3}, & \frac{r_{c}}{2} \leq r \leq r_{c} \\ 0, & otherwise \end{matrix}

Smoothing the Normals of Individual Crack Facets Using Least-Squares Approximation

After the predicted crack propagation direction is obtained based on the nonlocal stress/strain averaging, a moving least-squares approximation by polynomials is used by default to further smooth out the crack normals. The least-squares approximation is applied to the normals of the individual facets in elements along the crack front that satisfy the damage initiation criterion, as highlighted in Figure 12. This approximation provides a smoother crack surface (as shown in Figure 13), leading to a more accurate crack propagation direction.

You can use linear, quadratic, or cubic polynomial approximation for the moving least-squares approximation to smooth out the crack normals. You specify the number of terms in the polynomial. You can also suppress the least-squares approximation. In this case, the predicted crack propagation direction is determined based only on the nonlocal stress/strain averaging.

Limiting the Elements Involved in Crack Normal Smoothing

At the beginning of the analysis, you can choose to include or exclude the preexisting crack facets in elements from the moving least-squares approximation to obtain the crack propagation direction. During the analysis, you can also limit the elements involved in the least-squares approximation. You can set the maximum allowed difference (in degrees) below which the normals of the crack facets are included in the moving least-squares approximation. The default is 70°.

Damage Evolution

The damage evolution law describes the rate at which the cohesive stiffness is degraded once the corresponding initiation criterion is reached. The general framework for describing the evolution of damage is conceptually similar to that used for damage evolution in surface-based cohesive behavior (Contact Cohesive Behavior).

A scalar damage variable, D, represents the averaged overall damage at the intersection between the crack surfaces and the edges of cracked elements. It initially has a value of 0. If damage evolution is modeled, D monotonically evolves from 0 to 1 upon further loading after the initiation of damage. The normal and shear stress components are affected by the damage according to

t_{n} = {\begin{matrix} (1 - D) T_{n}, & T_{n} \geq 0 \\ T_{n}, & otherwise (no damage to compressive stiffness); \end{matrix}

t_{s} = (1 - D) T_{s},

t_{t} = (1 - D) T_{t},

where $T_{n}$ , $T_{s}$ , and $T_{t}$ are the normal and shear stress components predicted by the elastic traction-separation behavior for the current separations without damage.

To describe the evolution of damage under a combination of normal and shear separations across the interface, an effective separation is defined as

δ_{m} = \sqrt{{⟨ δ_{n} ⟩}^{2} + δ_{s}^{2} + δ_{t}^{2}} .

Use with User-Defined Damage Initiation Criterion

A separate damage evolution law should be specified for each damage initiation criterion defined in user subroutine UDMGINI. Each combination of a damage initiation criterion and a corresponding damage evolution law is referred to as a failure mechanism. Damage will accumulate for only one failure mechanism per element, corresponding to the mechanism whose damage initiation criterion was achieved first.

Use with LaRC05 Criterion

You can specify four separate damage evolution laws, one for each of the four initiation mechanisms. Alternatively, you can specify fewer than four damage evolution laws. In this case, the initiation mechanisms that do not have a corresponding evolution law use the specified damage evolution law with the smallest failure index. Damage accumulates for only one failure mechanism per element, corresponding to the mechanism whose damage initiation criterion was achieved first.

Viscous Regularization in Abaqus/Standard

Models exhibiting various forms of softening behavior and stiffness degradation often lead to severe convergence difficulties in Abaqus/Standard. Viscous regularization of the constitutive equations defining cohesive behavior in an enriched element can be used to overcome some of these convergence difficulties. Viscous regularization damping causes the tangent stiffness matrix to be positive definite for sufficiently small time increments.

The approximate amount of energy associated with viscous regularization over the whole model is available using output variable ALLVD.

Defining the Constitutive Response of Fluid Flow within the Cracked Element Surfaces

The formulae and laws that govern the behavior of fluid flow within the XFEM-based cracked element surfaces are very similar to those used for fluid flow within the cohesive element gap (Defining the Constitutive Response of Fluid within the Cohesive Element Gap). The similarities extend to the traction-separation model, damage initiation criteria, damage evolution law, and the fluid flow behavior. The fluid constitutive response includes the tangential flow within the cracked element surfaces and the normal flow across the cracked element surfaces due to caking or fouling effects in the enriched elements.

Tangential Flow

The tangential flow within the cracked element surfaces can be modeled with either a Newtonian or power-law model. By default, there is no tangential flow of pore fluid within the cracked element surfaces. To allow tangential flow, define a gap flow property in conjunction with the pore fluid material definition.

In the case of a Newtonian fluid the volume flow rate density vector is given by the expression

q d = - k_{t} \nabla p,

where $k_{t}$ is the tangential permeability (the resistance to the fluid flow), $\nabla p$ is the pressure gradient along the cracked element surfaces, and $d$ is the opening of the cracked element surfaces.

Abaqus defines the tangential permeability, or the resistance to flow, according to Reynold's equation:

k_{t} = \frac{d^{3}}{12 μ},

where $μ$ is the fluid viscosity and $d$ is the opening of the cracked element surfaces. You can also specify an upper limit on the value of $k_{t}$ .

In the case of a power law fluid the constitutive relation is defined as

τ = K {\dot{γ}}^{α},

where $τ$ is the shear stress, $\dot{γ}$ is the shear strain rate, $K$ is the fluid consistency, and $α$ is the power law coefficient. Abaqus defines the tangential volume flow rate density as

q d = - (\frac{2 α}{1 + 2 α}) {(\frac{1}{K})}^{\frac{1}{α}} {(\frac{d}{2})}^{\frac{1 + 2 α}{α}} {∥ \nabla p ∥}^{\frac{1 - α}{α}} \nabla p,

where $d$ is the opening of the cracked element surfaces.

By default, the gap between the cracked element surfaces has an initial opening of 0.002 in both a Newtonian fluid and a power law fluid. However, you can specify this opening directly.

Normal Flow across the Cracked Element Surfaces

You can permit normal flow by defining a fluid leak-off coefficient for the pore fluid material. This coefficient defines a pressure-flow relationship between the phantom nodes located at the cracked element edges and cracked element surfaces. The fluid leak-off coefficients can be interpreted as the permeability of a finite layer of material on the cracked element surfaces, as shown in Figure 14.

The normal flow is defined as

q_{t} = c_{t} (p_{i} - p_{t})

and

q_{b} = c_{b} (p_{i} - p_{b}),

where $q_{t}$ and $q_{b}$ are the flow rates into the top and bottom surfaces of a cracked element, respectively; $p_{i}$ is the pressure at the phantom node located at the cracked element edge; and $p_{t}$ and $p_{b}$ are the pore pressures on the top and bottom surfaces of a cracked element, respectively. You can optionally define leak-off coefficients as functions of temperature and field variables.

Alternatively, you can use user subroutine UFLUIDLEAKOFF to define more complex leak-off behavior, including the ability to define a time accumulated resistance, or fouling, through the use of solution-dependent state variables.

Thermal Effect of Fluid Flow within the Cracked Surfaces

If the thermal effect of the fluid flow within the cracked element surfaces is considered, the heat transport response comprises the tangential heat flux convection within the cracked surfaces, the normal heat convection between the gap fluid and the fracture surfaces, and the normal heat flux convected by fluid flow across the cracked element surfaces.

Optionally, you can include the normal heat convection between the gap fluid and the fracture surfaces.

The normal heat convection is defined as

q_{t} = h_{t} (θ_{i} - θ_{t})

and

q_{b} = h_{b} (θ_{i} - θ_{b}),

where $q_{t}$ and $q_{b}$ are the heat fluxes into the top and bottom surfaces of a cracked element, respectively; $θ_{i}$ is the temperature at the phantom node located on the cracked element edge; and $θ_{t}$ and $θ_{b}$ are the temperatures on the top and bottom surfaces of a cracked element, respectively.

Applying the VCCT Technique to the XFEM-Based LEFM Approach

The formulae and laws that govern the behavior of the XFEM-based linear elastic fracture mechanics approach for crack propagation analysis are very similar to those used for modeling delamination along a known and partially bonded surface (see Crack Propagation Analysis), where the strain energy release rate at the crack tip is calculated based on the modified Virtual Crack Closure Technique (VCCT). However, unlike this method, the XFEM-based LEFM approach can be used to simulate crack propagation along an arbitrary, solution-dependent path in the bulk material with or without an initial crack. You complete the definition of the crack propagation capability by defining a fracture-based surface behavior and specifying the fracture criterion in enriched elements.

Crack Nucleation and Direction of Crack Extension

By definition, the XFEM-based LEFM approach inherently requires the presence of a crack in the model since it is based on the principles of linear elastic fracture mechanics. The crack can be preexisting, or it can nucleate during the analysis. If there is no preexisting crack for a given enriched region, the XFEM-based LEFM approach is not activated until a crack nucleates. The crack nucleation is governed by one of the six built-in stress- or strain-based crack initiation criteria or a user-defined crack initiation criterion discussed in Crack Initiation and Direction of Crack Extension above. After a crack is nucleated in an enriched region, subsequent propagation of the crack is governed by the XFEM-based LEFM criterion.

Specifying When a Preexisting Crack Will Extend

If there is a preexisting crack in an enriched region, the crack extends after an equilibrium increment when the fracture criterion, f, reaches the value 1.0 within a given tolerance:

1.0 \leq f \leq 1.0 + f_{t o l} .

You can specify the tolerance $f_{t o l}$ . If $f > 1 + f_{t o l}$ , the time increment is cut back such that the crack extension criterion is satisfied. The default value of $f_{t o l}$ is 0.2.

Fracture of Multiple Elements in an Unstable Crack Growth Analysis

1.0 + f_{t o l} \leq f \leq 1.0 + f_{t o l}^{u},

where

f_{t o l}

is the tolerance described earlier in this section, the time increment size by default is immediately reduced automatically to a very small value,

β ▵ t_{\min} .

Reducing the time increment size allows more elements to fracture until

f < 1

for all the elements ahead of the crack tip.

However, you can optionally specify the maximum number of cutbacks allowed, $n$ , to be controlled by the regular tolerance, $f_{t o l}$ , prior to the activation of the unstable growth tolerance in an increment. After this limit is reached, the time increment size is recovered automatically to a larger value, $α ▵ t_{p r e}$ , where $▵ t_{\min}$ is the minimum time increment allowed; $▵ t_{p r e}$ is the time increment size prior to the unstable crack growth; and $α$ , $β$ , and $n$ are scaling parameters. The default values of $α$ , $β$ , and $n$ are 0.5, 2.0, and 0, respectively. If you do not specify a value for the unstable growth tolerance, the default value is infinity. In this case the fracture criterion, f, for unstable crack growth is not limited by any upper bound value in the above equation.

Specifying the Crack Propagation Direction

You must specify the crack propagation direction when the fracture criterion is satisfied. The crack can extend at a direction normal to the direction of the maximum tangential stress, orthogonal to the element local 1-direction (see Conventions), or orthogonal to the element local 2-direction. By default, the crack propagates normal to the direction of the maximum tangential stress.

Limiting the Crack Propagation Direction

If the crack direction normal to the maximum tangential stress is specified, you can limit the new crack propagation direction to within a certain angle (in degrees) of the previous crack propagation direction. The default is 85°.

Nonlocal Smoothing of the Crack Normals

After the normals of the individual crack facets are obtained based on the fracture criterion defined above, you can use a moving least-squares approximation by polynomials to further smooth out the crack normals. The least-squares approximation is applied to the normals of the individual facets in elements along the crack front that satisfy the fracture criterion to obtain a more accurate crack propagation direction.

Specifying the Approximation Used in the Least-Squares Approximation

You can use linear, quadratic, or cubic polynomial approximation for the moving least-squares approximation to smooth out the crack normals. You specify the number of terms in the polynomial.

Specifying the Region of the Model Used for Nonlocal Smoothing of the Crack Normals

To control the range of elements used for nonlocal smoothing of the crack normals in the crack direction calculations, you can specify a radius, $r_{c}$ , within which the elements around the crack tip along the crack front are included. The default radius is three times the typical element characteristic length along the crack front in the enriched region.

Mixed-Mode Behavior

Abaqus provides three common mode-mix formulae for computing the equivalent fracture energy release rate $G_{e q u i v C}$ : the BK law, the power law, and the Reeder law models. The choice of model is not always clear in any given analysis; an appropriate model is best selected empirically.

BK Law

The BK law model is described in Benzeggagh and Kenane (1996) by the following formula:

G_{e q u i v C} = G_{I C} + (G_{I I C} - G_{I C}) {(\frac{G_{I I} + G_{I I I}}{G_{I} + G_{I I} + G_{I I I}})}^{η},

G_{e q u i v} = G_{I} + G_{I I} + G_{I I I} .

To define this model, you must provide $G_{I C}, G_{I I C},$ and $η$ . This model provides a power law relationship combining energy release rates in Mode I, Mode II, and Mode III into a single scalar fracture criterion.

Power Law

The power law model is described in Wu and Reuter (1965) by the following formula:

\frac{G_{e q u i v}}{G_{e q u i v C}} = {(\frac{G_{I}}{G_{I C}})}^{a_{m}} + {(\frac{G_{I I}}{G_{I I C}})}^{a_{n}} + {(\frac{G_{I I I}}{G_{I I I C}})}^{a_{o}} .

To define this model, you must provide $G_{I C}, G_{I I C}, G_{I I I C}, a_{m}, a_{n},$ and $a_{o}$ .

Reeder Law

The Reeder law model is described in Reeder et al. (2002) by the following formula:

G_{e q u i v C} = G_{I C} + (G_{I I C} - G_{I C}) {(\frac{G_{I I} + G_{I I I}}{G_{I} + G_{I I} + G_{I I I}})}^{η} +

(G_{I I I C} - G_{I I C}) (\frac{G_{I I I}}{G_{I I} + G_{I I I}}) {(\frac{G_{I I} + G_{I I I}}{G_{I} + G_{I I} + G_{I I I}})}^{η},

G_{e q u i v} = G_{I} + G_{I I} + G_{I I I} .

To define this model, you must provide $G_{I C}, G_{I I C}, G_{I I I C},$ and $η$ . The Reeder law is best applied when $G_{I I C} \neq G_{I I I C}$ ; when $G_{I I C} = G_{I I I C}$ , the Reeder law reduces to the BK law. The Reeder law applies only to three-dimensional problems.

Defining Variable Critical Energy Release Rates

You can define a VCCT criterion with varying energy release rates by specifying the critical energy release rates at the nodes.

If you indicate that the nodal critical energy rates will be specified, any constant critical energy release rates you specify are ignored and the critical energy release rates are interpolated from the nodes. The critical energy release rates must be defined at all nodes in the enriched region.

Enhanced VCCT Criterion

The formulae and laws governing the behavior of the enhanced VCCT criterion are very similar to those used for the VCCT criterion. However, unlike the VCCT criterion, the onset and growth of a crack can be controlled by two different critical fracture energy release rates: $G_{C}$ and $G_{C}^{P}$ . In a general case involving Mode I, II, and III fracture, when the fracture criterion is satisfied; i.e,

f = \frac{G_{e q u i v}}{G_{e q u i v C}} \geq 1.0,

the traction on the two surfaces of the cracked element is ramped down over the separation with the dissipated strain energy equal to the critical equivalent strain energy required to propagate the crack, $G_{e q u i v C}^{P}$ , rather than the critical equivalent strain energy required to initiate the separation, $G_{e q u i v C}$ . The formulae for calculating $G_{e q u i v C}^{P}$ are identical to those used for $G_{e q u i v C}$ for different mixed-mode fracture criteria.

Fatigue Crack Growth Criterion Based on the Principles of LEFM

If you specify the fatigue crack growth criterion, progressive crack growth at the enriched elements subjected to sub-critical cyclic loading can be simulated. This criterion can be used only in a fatigue crack growth analysis using the direct cyclic approach (Low-Cycle Fatigue Analysis Using the Direct Cyclic Approach) or the general fatigue crack growth approach (Linear Elastic Fatigue Crack Growth Analysis). A fatigue crack growth analysis step can be the only step, can follow a general static step, or can be followed by a general static step. You can include multiple fatigue crack growth analysis steps in a single analysis. If you perform a fatigue analysis in a model without a preexisting crack, you must precede the fatigue step with a static step that nucleates a crack, as discussed in Crack Nucleation and Direction of Crack Extension.

The onset and fatigue crack growth are characterized by using the Paris law, which relates the fracture energy release rate or the stress intensity factor to crack growth rates. The fracture energy release rates or the stress intensity factors at the crack tips in the enriched elements are calculated based on the above mentioned VCCT technique.

The Paris regime is bounded by the energy release rate threshold, $G_{t h r e s h}$ , below which there is no consideration of fatigue crack initiation or growth, and the energy release rate upper limit, $G_{p l}$ , above which the fatigue crack will grow at an accelerated rate. $G_{C}$ is the critical equivalent strain energy release rate calculated based on the user-specified mode-mix criterion and the fracture strength of the bulk material. The formulae for calculating $G_{C}$ have been provided above for different mixed-mode fracture criteria. You can specify the ratio of $G_{t h r e s h}$ over $G_{C}$ and the ratio of $G_{p l}$ over $G_{C}$ . The default values are $\frac{G_{t h r e s h}}{G_{C}} = 0.01$ and $\frac{G_{p l}}{G_{C}} = 0.85$ .

Onset of Fatigue Crack Growth

The onset of fatigue crack growth refers to the beginning of fatigue crack growth at the crack tip in the enriched elements. In a fatigue crack growth analysis the onset of the fatigue crack growth criterion is characterized by $Δ G$ , which is the relative fracture energy release rate when the structure is loaded between its maximum and minimum values. The fatigue crack growth initiation criterion is defined as

f = \frac{N}{c_{1} Δ G^{c_{2}}} \geq 1.0,

where $c_{1}$ and $c_{2}$ are material constants and $N$ is the cycle number. The enriched elements ahead of the crack tips will not be fractured unless the above equation is satisfied and the maximum fracture energy release rate, $G_{m a x}$ , which corresponds to the cyclic energy release rate when the structure is loaded up to its maximum value, is greater than $G_{t h r e s h}$ . If you do not specify the onset criterion, Abaqus/Standard assumes that the onset of fatigue crack growth is satisfied automatically.

Fatigue Crack Growth Using the Paris Law

Once the onset of the fatigue crack growth criterion is satisfied at the enriched element, the crack growth rate, $d a / d N$ , can be calculated based on the relative fracture energy release rate, $Δ G$ . The rate of the crack growth per cycle is given by the Paris law if $G_{t h r e s h} < G_{m a x} < G_{p l},$

\frac{d a}{d N} = c_{3} Δ G^{c_{4}},

where $c_{3}$ and $c_{4}$ are material constants.

Ratcliffe and Johnston (2014) and Deobald et al. (2017) proposed the following alternate form of the fatigue law which better accounts for mixed-mode fatigue crack growth:

\frac{d a}{d N} = c_{3} G_{T M a x}^{c_{4}} .

In the above expression, $G_{T M a x}$ is the total maximum strain energy release rate (as opposed to the strain energy release rate change over a cycle used in the original form of the Paris law), while $c_{3}$ and $c_{4},$ are material parameters that depend on mode-mix and stress ratios. Abaqus does not support the above form of the crack growth rate equation directly, but instead allows specification of $d a / d N$ as a tabular function of $G_{T M a x}$ , the mode-mix ratio, and the stress ratio.

In addition, user subroutine UMIXMODEFATIGUE provides a general capability for implementing a user-defined fatigue crack growth law.

At the end of cycle $N$ , Abaqus/Standard extends the crack length, $a_{N}$ , from the current cycle forward over an incremental number of cycles, $Δ N$ to $a_{N + Δ N}$ by fracturing at least one enriched element ahead of the crack tips. Given the material constants $c_{3}$ and $c_{4}$ , combined with the known element length and the likely crack propagation direction $Δ a_{N}_{j} = a_{N + Δ N} - a_{N}$ at the enriched elements ahead of the crack tips, the number of cycles necessary to fail each enriched element ahead of the crack tip can be calculated as $Δ N_{j}$ , where j represents the enriched element ahead of the $j$ th crack tip. The analysis is set up to advance the crack by at least one enriched element after the loading cycle is completed. The element with the fewest cycles is identified to be fractured, and its $Δ N_{m i n} = m i n (Δ N_{j})$ is represented as the number of cycles to grow the crack equal to its element length, $Δ a_{N}_{m i n} = m i n (Δ a_{N}_{j})$ . The most critical element is completely fractured with a zero constraint and a zero stiffness at the end of the completed cycle. As the enriched element is fractured, the load is redistributed and a new relative fracture energy release rate must be calculated for the enriched elements ahead of the crack tips for the next cycle. This capability allows at least one enriched element ahead of the crack tips to be fractured completely after each completed cycle and precisely accounts for the number of cycles needed to cause fatigue crack growth over that length.

If $G_{m a x} > G_{p l}$ , the enriched elements ahead of the crack tips will be fractured by increasing the cycle number count, $d N$ , by one only.

For information on how to accelerate the fatigue crack growth analysis and to provide a smooth solution for the crack front, see Controlling Element Fracture.

For linear elastic materials, the fracture energy release rate is related directly to the stress intensity factors by the following relationship:

G = \frac{K_{I}^{2}}{E^{'}} + \frac{K_{I I}^{2}}{E^{'}} + \frac{K_{I I I}^{2}}{2 μ},

where

E^{'}

is the Young's modulus, and

μ

is the shear modulus. A form based on the stress intensity factor,

K

, (equivalent to the fracture energy release rate–based form above) is available using

\frac{d a}{d N} = c_{3} Δ K_{e f f}^{c_{4}},

where

c_{3}

and

c_{4}

are material constants, and

Δ K_{e f f}

is the effective stress intensity factor range of a load cycle.

A mode-mix formula for computing the effective stress intensity factor is described in Irwin (1968) as follows:

Δ K_{e f f} = \sqrt{A * Δ K_{I}^{2} + B * Δ K_{I I}^{2} + C * \frac{1}{1 - ν} Δ K_{I I I}^{2}},

where A, B, and C are user-defined material constants; and

ν

is the Poisson's ratio.

For the fatigue crack growth criterion, the following forms based on the stress intensity factor are also available:

A tabular form to support multiple piecewise linear log-log ( $\frac{d a}{d N}$ versus $Δ K_{e f f}$ ) segments.
A user-defined crack growth criterion using user subroutine UMIXMODEFATIGUE.

Viscous Regularization for the XFEM-Based LEFM Approach

The simulation of structures with unstable propagating cracks is challenging and difficult. Nonconvergent behavior may occur from time to time. Localized damping is included for the XFEM-based LEFM approach by using the viscous regularization technique. Viscous regularization damping causes the tangent stiffness matrix of the softening material to be positive for sufficiently small time increments.

Applying Distributed Pressure Loads and Heat Fluxes to Cracked Element Surfaces

When an element is cut by a crack during the analysis, a XFEM-based crack surface is generated during the analysis (see Defining a Crack Surface above). You can apply a distributed pressure load or heat flux to the cracked element surfaces.

Specifying the Initial Location of an Enriched Feature

Because the mesh is not required to conform to the geometric discontinuities, the initial location of a preexisting crack must be specified in the model. The level set method is provided for this purpose. Two signed distance functions per node are generally required to describe a crack geometry. The first describes the crack surface, while the second is used to construct an orthogonal surface so that the intersection of the two surfaces gives the crack front (see Initial Conditions).

The first signed distance function must be either greater or less than zero and cannot be equal to zero. If an initial crack has to be defined at the boundaries of an element, a very small positive or negative value for the first signed distance function must be specified.

Activating and Deactivating the Enriched Feature

You can activate or deactivate the crack propagation capability within the step definition.

You can specify a small relative radius, $r_{s}$ , around the crack tip (as shown in Figure 15) to define a crack initiation suppression zone within which the elements in the enriched feature are excluded from crack nucleation. However, multiple crack nucleations are allowed to occur when the damage initiation criterion is satisfied in the elements lying outside the crack initiation suppression zone in the enriched feature. The default radius value is five times the typical element characteristic length in the enriched region.

Contour Integral

When you evaluate the contour integrals using the conventional finite element method (Contour Integral Evaluation), you must define the crack front explicitly and specify the virtual crack extension direction in addition to matching the mesh to the cracked geometry. Detailed focused meshes are generally required and obtaining accurate contour integral results for a crack in a three-dimensional curved surface can be cumbersome. The extended finite element in conjunction with the level set method alleviates these shortcomings. The adequate singular asymptotic fields and the discontinuity are ensured by the special enrichment functions in conjunction with additional degrees of freedom. In addition, the crack front and the virtual crack extension direction are determined automatically by the level set signed distance functions.

Specifying the Enrichment Radius

Although XFEM has alleviated the shortcomings associated with defining the conforming mesh in the neighborhood of the crack front, you must still generate a sufficient number of elements around the crack front to obtain path-independent contours, particularly in a region with high crack-front curvature. The group of elements within a small radius from the crack front are enriched and become involved in the contour integral calculations. The default enrichment radius is six times the typical element characteristic length of those elements along the crack front in the enriched area. You must include the elements inside the enrichment radius in the element set used to define the enriched region. Elements with a large aspect ratio should be avoided along the crack front.

Procedures

Modeling discontinuities as an enriched feature can be performed using any of the following:

static analysis (see Static Stress Analysis);
implicit dynamic analysis (see Implicit Dynamic Analysis Using Direct Integration);
low-cycle fatigue analysis using the direct cyclic approach (Low-Cycle Fatigue Analysis Using the Direct Cyclic Approach);
fatigue crack growth analysis (see Linear Elastic Fatigue Crack Growth Analysis);
fully coupled thermal-stress analysis (see Fully Coupled Thermal-Stress Analysis);
quasi-static analysis (see Quasi-Static Analysis);
geostatic stress field analysis (see Geostatic Stress State); or
coupled pore fluid diffusion/stress analysis (see Coupled Pore Fluid Diffusion and Stress Analysis).

Initial Conditions

Initial conditions to identify initial boundaries or interfaces of an enriched feature can be specified (see Initial Conditions).

Boundary Conditions

Boundary conditions can be applied to any of the displacement, temperature, or pore pressure degrees of freedom (see Boundary Conditions).

Loads

The following types of loading can be prescribed in a model with an enriched feature:

Concentrated nodal forces can be applied to the displacement degrees of freedom (1–3) or the pore pressure degree of freedom (8); see Concentrated Loads.
Distributed pressure forces or body forces can be applied; see Distributed Loads. The distributed load types available with particular elements are described in Abaqus Elements Guide.
Concentrated heat fluxes.
Body fluxes and distributed surface fluxes.
Node-based film and radiation conditions.
Average-temperature radiation conditions.
Element- and surface-based film and radiation conditions.

For more information on heat fluxes, film conditions, and radiation conditions, see Thermal Loads.

Predefined Fields

The following predefined fields can be specified in a model with an enriched feature, as described in Predefined Fields:

Nodal temperatures (although temperature is not a degree of freedom in stress/displacement elements). The specified temperature affects temperature-dependent critical stress and strain failure criteria.
The values of user-defined field variables. The specified value affects field-variable-dependent material properties.

Material Options

Any of the mechanical constitutive models in Abaqus/Standard, including user-defined materials (defined using user subroutine UMAT) can be used to model the mechanical behavior of the enriched element in a crack propagation analysis. See Abaqus Materials Guide. The inelastic definition at a material point must be used in conjunction with the linear elastic material model (Linear Elastic Behavior), the porous elastic material model (Elastic Behavior of Porous Materials), or the hypoelastic material model (Hypoelastic Behavior). Only isotropic elastic materials are supported when evaluating the contour integral for a stationary crack.

Elements

Only the following elements can be associated with an enriched feature:

First-order solid continuum stress/displacement elements.
First-order displacement/pore pressure solid continuum elements.
First-order displacement/temperature solid continuum elements.
First-order displacement/pore pressure/temperature solid continuum elements.
Second-order stress/displacement tetrahedral elements.

For propagating cracks these include

bilinear plane strain and plane stress elements,
bilinear axisymmetric elements,
linear brick elements,
linear tetrahedral elements, and
second-order tetrahedral elements.

For stationary cracks, these include

linear brick elements,
linear tetrahedral elements, and
second-order tetrahedral elements.

For an incompatible mode element, Abaqus/Standard discards the contribution due to the incompatible deformation mode immediately after the element is fractured under a tensile loading. Therefore, the stress level at the cracked element may not return completely to its originally unloaded state even when this cracked element is unloaded completely and the contact of the cracked element surfaces is reestablished.

Output

In addition to the standard output identifiers available in Abaqus (Abaqus/Standard Output Variable Identifiers), the following nodal, whole element, and surface variables have special meaning for a model with an enriched feature.

Nodal variables:

PHILSM: Signed distance function to describe the crack surface.
PSILSM: Signed distance function to describe the initial crack front.

Whole element variables:

STATUSXFEM: Status of the enriched element. (The status of an enriched element is 1.0 if the element is completely cracked and 0.0 if the element contains no crack. If the element is partially cracked, the value of STATUSXFEM lies between 1.0 and 0.0.)
ENRRTXFEM: All components of strain energy release rate when linear elastic fracture mechanics with the extended finite element method is used.
LOADSXFEM: Distributed pressure loads applied to the crack surface.
CYCLEINIXFEM: Minimum number of cycles needed to satisfy the condition for the onset of fatigue crack growth at an enriched element.

The following whole element output variables are available only when fluid flow is enabled within the cracked enriched element surfaces:

GFVRXFEM: Gap fluid volume rate of the enriched element.
CRDCUTXFEM: Crack midpoint coordinates at the element edges of the enriched element.
PFOPENXFEM: Fracture opening of the enriched element.
PFOPENXFEMCOMP: Fracture opening at the element edges of the enriched element.
PORPRES: Fluid pressure of the enriched element.
PORPRESCOMP: Fluid pressure at the element edges of the enriched element.
LEAKVRTXFEM: Leak-off flow rate at the top of the enriched element.
LEAKVRBXFEM: Leak-off flow rate at the bottom of the enriched element.
ALEAKVRTXFEM: Accumulated leak-off flow volume at the top of the enriched element.
ALEAKVRBXFEM: Accumulated leak-off flow volume at the bottom of the enriched element.

Surface variables (available only for propagating cracks modeled with first-order solid continuum elements):

CRKDISP: Crack opening and relative tangential motions on cracked surfaces in enriched elements.
CSDMG: Damage variable on cracked surfaces in enriched elements.
CRKSTRESS: Remaining residual pressure and tangential shear stresses on cracked surfaces in enriched elements.

The following surface output variables are available only when fluid flow is enabled within the cracked enriched element surfaces:

GFVR: Fluid volume rate within the cracked surfaces in the enriched element.
PORPRES: Pore pressure within the cracked surfaces in the enriched element.
PORPRESURF: Pore pressure on the cracked surfaces in the enriched element.
LEAKVR: Leak-off flow rate on the cracked surfaces in the enriched element.
ALEAKVR: Accumulated leak-off flow volume on the cracked surfaces in the enriched element.

The following surface output variable is available only when both fluid flow and thermal convection are enabled within the cracked enriched element surfaces:

PORTEMP: Temperature within the cracked surfaces in the enriched element.

Use of Unsymmetric Matrix Storage and Solution

When the pore pressure degrees of freedom are activated in the enriched elements, matrices are unsymmetric; therefore, unsymmetric matrix storage and solution may be needed to improve convergence (see Matrix Storage and Solution Scheme in Abaqus/Standard).

Limitations

The following limitations exist with an enriched feature:

An enriched element cannot be intersected by more than one crack.
A crack is not allowed to turn more than 90° in one increment during an analysis.
Only asymptotic crack-tip fields in an isotropic elastic material are considered for a stationary crack.
Adaptive remeshing is not supported.
Composite solid elements are not supported.
Import analysis is not supported.

Input File Template

The following is an example of modeling crack propagation with the XFEM-based cohesive segments method:

HEADING
...
NODE, NSET=ALL
...
ELEMENT, TYPE=C3D8, ELSET=REGULAR
ELEMENT, TYPE=C3D8, ELSET=ENRICHED
...
SOLID SECTION, MATERIAL=STEEL1, ELSET=REGULAR
SOLID SECTION, MATERIAL=STEEL12, ELSET=ENRICHED

ENRICHMENT, TYPE=PROPAGATION CRACK, ELSET=ENRICHED, 
NAME=ENRICHMENT, INTERACTION=INTERACTION
SURFACE, TYPE=XFEM, NAME=SURF_NAME
Data lines to specify the names of enriched features
MATERIAL, NAME=STEEL1
...
MATERIAL, NAME=STEEL2
DAMAGE INITIATION, CRITERION=MAXPS, TOLERANCE=0.05
DAMAGE EVOLUTION, TYPE=ENERGY
Data lines to specify the failure mechanism
...
SURFACE INTERACTION, NAME=INTERACTION
SURFACE BEHAVIOR
Data lines to specify the contact of cracked element surfaces
...
STEP
STATIC
...
END STEP
STEP
STATIC
...

ENRICHMENT ACTIVATION, TYPE=PROPAGATION CRACK, 
NAME=ENRICHMENT, ACTIVATE=OFF
...
END STEP

The following is an example of modeling crack propagation with the XFEM-based LEFM approach:

HEADING
...
NODE, NSET=ALL
...
ELEMENT, TYPE=C3D8, ELSET=REGULAR
ELEMENT, TYPE=C3D8, ELSET=ENRICHED
...
SOLID SECTION, MATERIAL=STEEL1, ELSET=REGULAR
SOLID SECTION, MATERIAL=STEEL12, ELSET=ENRICHED

ENRICHMENT, TYPE=PROPAGATION CRACK, ELSET=ENRICHED, 
NAME=ENRICHMENT, INTERACTION=INTERACTION
MATERIAL, NAME=STEEL1
...
MATERIAL, NAME=STEEL2
DAMAGE INITIATION, CRITERION=MAXPS, TOLERANCE=0.05
Data lines to specify the crack nucleation mechanism
...
SURFACE INTERACTION, NAME=INTERACTION
SURFACE BEHAVIOR
FRACTURE CRITERION, TYPE=VCCT, TOLERANCE=0.05,VISCOSITY=0.00001
Data lines to specify the crack propagation criterion
...
END STEP

The following is an example of calculating contour integrals in stationary cracks with the extended finite element method:

HEADING
...
NODE, NSET=ALL
...
ELEMENT, TYPE=C3D8, ELSET=REGULAR
ELEMENT, TYPE=C3D8, ELSET=ENRICHED
...
SOLID SECTION, MATERIAL=STEEL1, ELSET=REGULAR
SOLID SECTION, MATERIAL=STEEL12, ELSET=ENRICHED

ENRICHMENT, TYPE=STATIONARY CRACK, ELSET=ENRICHED, 
NAME=ENRICHMENT, ENRICHMENT RADIUS
MATERIAL, NAME=STEEL1
...
MATERIAL, NAME=STEEL2
...
STEP
STATIC
...
CONTOUR INTEGRAL, CRACK NAME=ENRICHMENT, XFEM
END STEP

References

Belytschko, T., and T. Black, “Elastic Crack Growth in Finite Elements with Minimal Remeshing,” International Journal for Numerical Methods in Engineering, vol. 45, pp. 601–620, 1999.
Benzeggagh, M., and M. Kenane, “Measurement of Mixed-Mode Delamination Fracture Toughness of Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending Apparatus,” Composite Science and Technology, vol. 56 439, 1996.
Deobald, L., G. Mabson, S. Engelstad, M. Rao, M. Gurvich, W. Seneviratne, S. Perera, T. O'Brien, G. Murri, J. Ratcliffe, C. Davila, N. Carvalho , and R. Krueger, “Guidelines for VCCT-Based Interlaminar Fatigue and Progressive Failure Finite Element Analysis,” NASA/TM-2017-219663, 2017.
Elguedj, T., A. Gravouil, and A. Combescure, “Appropriate Extended Functions for X-FEM Simulation of Plastic Fracture Mechanics,” Computer Methods in Applied Mechanics and Engineering, vol. 195, pp. 501–515, 2006.
Irwin, G. R., “Linear Fracture Mechanics Fracture Transition, and Fracture Control,” Engineering Fracture Mechanics, vol. 1, pp. 241–257, 1968.
Melenk, J., and I. Babuska, “The Partition of Unity Finite Element Method: Basic Theory and Applications,” Computer Methods in Applied Mechanics and Engineering, vol. 39, pp. 289–314, 1996.
Ratcliffe, J., and W. Johnston, “Influence of Mixed Mode I-Mode II Loading on Fatigue Delamination Growth Characteristics of a Graphite Epoxy Tape Laminate,” Proceedings of American Society for Composites 29th Technical Conference, 2014.
Reeder, J., S. Kyongchan, P. B. Chunchu, and D. R.. Ambur, “Postbuckling and Growth of Delaminations in Composite Plates Subjected to Axial Compression”43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Denver, Colorado, vol. 1746, p. 10, 2002.
Remmers, J. J. C., R. de Borst, and A. Needleman, “The Simulation of Dynamic Crack Propagation using the Cohesive Segments Method,” Journal of the Mechanics and Physics of Solids, vol. 56, pp. 70–92, 2008.
Song, J. H., P. M. A. Areias, and T. Belytschko, “A Method for Dynamic Crack and Shear Band Propagation with Phantom Nodes,” International Journal for Numerical Methods in Engineering, vol. 67, pp. 868–893, 2006.
Sukumar, N., Z. Y. Huang, J.-H. Prevost, and Z. Suo, “Partition of Unity Enrichment for Bimaterial Interface Cracks,” International Journal for Numerical Methods in Engineering, vol. 59, pp. 1075–1102, 2004.
Sukumar, N., and J.-H. Prevost, “Modeling Quasi-Static Crack Growth with the Extended Finite Element Method Part I: Computer Implementation,” International Journal for Solids and Structures, vol. 40, pp. 7513–7537, 2003.
Ventura, G., and E. Benvenuti, “Equivalent Polynomials for Quadrature in Heaviside Function Enriched Elements,” International Journal for Numerical Methods in Engineering, vol. 102, pp. 688–710, 2015.
Wu, E. M., and R. C. Reuter Jr., “Crack Extension in Fiberglass Reinforced Plastics,” T and M Report, University of Illinois, vol. 275, 1965.