Contour Integral Evaluation

Abaqus/Standard offers two different ways to evaluate the contour integral. The first approach is based on the conventional finite element method, which typically requires you to conform the mesh to the cracked geometry, to explicitly define the crack front, and to specify the virtual crack extension direction. Detailed focused meshes are generally required, and obtaining accurate contour integral results for a crack in a three-dimensional curved surface can be quite cumbersome. The extended finite element method (XFEM) alleviates these shortcomings. XFEM does not require the mesh to match the cracked geometry. The presence of a crack is ensured by the special enriched functions in conjunction with additional degrees of freedom. You must, however, generate a sufficient number of elements around the crack front to obtain path-independent contours, particularly in the region with high crack front curvature. This approach also removes the requirement for explicitly defining the crack front or specifying the virtual crack extension direction when evaluating the contour integral. The data required for the contour integral are determined automatically based on the level set signed distance functions at the nodes in an element (see Modeling Discontinuities as an Enriched Feature Using the Extended Finite Element Method).

Several contour integral evaluations are possible at each location along a crack. In a finite element model each evaluation can be thought of as the virtual motion of a block of material surrounding the crack tip (in two dimensions) or surrounding each node along the crack line (in three dimensions). Each block is defined by contours, where each contour is a ring of elements completely surrounding the crack tip or the nodes along the crack line from one crack face to the opposite crack face. These rings of elements are defined recursively to surround all previous contours.

Abaqus/Standard automatically finds the elements that form each ring from the regions defined as the crack tip or crack line. Each contour provides an evaluation of the contour integral. The possible number of evaluations is the number of such rings of elements for two-dimensional quadrilateral and three-dimensional brick elements. For tetrahedral elements, you must specify a small radius within which rings of elements are identified for fracture mechanics studies. A refined mesh is required to define the rings of elements around the crack front, especially in a region near the external free surfaces. In a case where the crack front intersects the external free surface in a model with tetrahedral elements at an angle not equal to 90°, you should specify surface normals at all the crack tip nodes that lie on the external free surfaces (see Normal Definitions at Nodes). This action ensures that the tangential directions of the crack front at those locations are estimated accurately for contour integral evaluation. The default value of the ring radius is twice the typical element characteristic length along the crack front, which works well for most problems. You must specify the number of contours to use in calculating contour integrals. In addition, you must specify the type of contour integral to calculate, as described below. By default, Abaqus/Standard calculates the J-integral.

You can assign a name to a crack that is used to identify the contour integral values in the data file and in the output database file. If you are using the conventional finite element method and do not specify a crack name, by default Abaqus/Standard generates crack numbers that follow the order in which the cracks are defined. If you are using XFEM, you must set the crack name equal to the name assigned to the enriched feature. Both the domain integral method and the line integral method are supported when you evaluate the contour integral using XFEM.

Domain Integral Method

Using the divergence theorem, the contour integral can be expanded into an area integral in two dimensions or a volume integral in three dimensions, over a finite domain surrounding the crack. This domain integral method is used to evaluate contour integrals in Abaqus/Standard. The method is quite robust in the sense that accurate contour integral estimates are usually obtained even with quite coarse meshes. The method is robust because the integral is taken over a domain of elements surrounding the crack and because errors in local solution parameters have less effect on the evaluated quantities such as J, $C_{t}$ , the stress intensity factors, and the T-stress.

Requesting Multiple Contour Integrals

Contour integrals at several different crack tips in two dimensions or along several different crack lines in three dimensions can be evaluated at any time by repeating the contour integral request as often as needed in the step definition. When you are using the conventional finite element method, you must specify the crack front and the direction of virtual crack extension (or the normal to the crack plane if this normal is constant) for each crack tip or crack line, as described below. When you are using XFEM, you do not need to specify the crack front or the virtual crack extension direction because they will be determined by Abaqus/Standard. However, you must set each crack name equal to the corresponding enriched feature, with each enriched feature consisting of only one crack. In addition, regardless of whether you are using either the conventional finite element method or XFEM, you must specify the number of contours to be calculated for each integral.

J-Integral

The J-integral is usually used in rate-independent quasi-static fracture analysis to characterize the energy release associated with crack growth. It can be related to the stress intensity factor if the material response is linear.

The J-integral is defined in terms of the energy release rate associated with crack advance. For a virtual crack advance $λ (s)$ in the plane of a three-dimensional fracture, the energy release rate is given by

\begin{matrix} \bar{J} & ​ = \int_{A} λ (s) n \cdot H \cdot q d A \end{matrix},

where $d A$ is a surface element along a vanishing small tubular surface enclosing the crack tip or crack line, $n$ is the outward normal to $d A$ , and $q$ is the local direction of virtual crack extension. $H$ is given by

H = (W I - σ \cdot \frac{\partial u}{\partial x}) .

For elastic material behavior W is the elastic strain energy; for elastic-plastic or elasto-viscoplastic material behavior W is defined as the elastic strain energy density plus the plastic dissipation, thus representing the strain energy in an “equivalent elastic material.” Therefore, the J-integral calculated is suitable only for monotonic loading of elastic-plastic materials.

Domain Dependence

The J-integral should be independent of the domain used provided that the crack faces are parallel to each other, but J-integral estimates from different rings may vary because of the approximate nature of the finite element solution. Strong variation in these estimates, commonly called domain dependence or contour dependence, typically indicates an error in the contour integral definition. Gradual variation in these estimates may indicate that a finer mesh is needed or, if plasticity is included, that the contour integral domain does not completely include the plastic zone. If the “equivalent elastic material” is not a good representation of the elastic-plastic material, the contour integrals will be domain independent only if they completely include the plastic zone. Since it is not always possible to include the plastic zone in three dimensions, a finer mesh may be the only solution.

If the first contour integral is defined by specifying the nodes at the crack tip, the first few contours may be inaccurate. To check the accuracy of these contours, you can request more contours and determine the value of the contour integral that appears approximately constant from one contour to the next. The contour integral values that are not approximately equal to this constant should be discarded. In linear elastic problems the first and second contours typically should be ignored as inaccurate.

For some three-dimensional models with an open crack front, the J-integral estimates may be inaccurate from the node sets (or elements in the case with XFEM) at the crack front ends. The resolution difficulty is compounded by the skewness of the outmost layer of elements. This accuracy loss is confined only to the contour integrals at the front ends and has no effect on the accuracy of the contour integral values at the neighboring node sets (or elements in the case with XFEM) along the crack front.

Including the Effect of a Residual Stress Field on J-Integral Evaluation

A residual stress field often occurs in a structure; for example, as a result of service loads that produce plasticity, a metal forming process in the absence of an anneal treatment, thermal effects, or swelling effects. When the residual stresses are significant, the standard definition of the J-integral as described above may lead to a path-dependent value. To ensure its path independence, the J-integral evaluation must include an additional term that accounts for the residual stress field. In Abaqus/Standard the problem with a residual stress field is treated as an initial strain problem. If the total strain is written as the sum of mechanical strain, $ε^{m}$ , and initial strain, $ε^{o}$ ; that is,

ε = ε^{m} + ε^{o},

a path-independent energy release rate in the presence of a residual stress field is given by

\begin{matrix} \bar{J} & = \int_{A} λ (s) n \cdot (W I - σ \cdot \frac{\partial u}{\partial x}) \cdot q d A + \int_{V} σ : \frac{\partial ε^{o}}{\partial x} \cdot q d V \end{matrix},

where V is the domain volume enclosing the crack tip or crack line, W is defined as the mechanical strain energy density only,

W = \int_{0}^{ε^{m}} σ : d ε^{m},

and $ε^{o}$ remains constant during the entire deformation.

The residual stress field can be specified by reading the stress data from a previous analysis step or by defining an initial condition (see Defining Initial Stresses). You specify the step number from which the stress data in the last available increment of the specified step will be considered as residual stresses. If the step number is set equal to zero (default), the residual stress field is defined by the initial condition definition. When XFEM is used, the residual stress field can be defined only with an initial condition definition.

The approach used to account for the residual stress effects was adopted from Lei, O'Dowd, and Webster (2000).

Allowing Non-Proportional Stressing

Usually, the proportional stressing condition is assumed to be satisfied when the contour integrals are calculated. In other words, this equation prevails:

σ = \frac{\partial W}{\partial ε^{m}} .

Otherwise, you should modify the path-independent energy release rate as follows:

\begin{matrix} \bar{J} & = \int_{A} λ (s) n \cdot (W I - σ \cdot \frac{\partial u}{\partial x}) \cdot q d A + \int_{V} (σ : \frac{\partial ε^{m}}{\partial x} - \frac{\partial W}{\partial x}) \cdot q d V . \end{matrix}

This approach used to account for the non-proportional stressing effects was adopted from Lei (2005).

C_t-Integral

The C_t-integral is supported with the conventional finite element method; however, it is not supported with XFEM.

The $C_{t}$ -integral can be used for time-dependent creep behavior, where it characterizes creep crack deformation under certain creep conditions, including transient crack growth. $C_{t}$ is, for example, proportional to the rate of growth of the crack-tip/crack-line creep zone for a stationary crack under small-scale creep conditions. Under steady-state creep conditions, when creep dominates throughout the specimen, $C_{t}$ becomes path independent and is known as $C^{*}$ . $C_{t}$ -integrals should be requested only in a quasi-static step.

The $C_{t}$ -integral is obtained by replacing the displacements with velocities and the strain energy density with the strain energy rate density in the J-integral expansion. The strain energy rate density is defined as

\dot{W} \overset{def}{=} \int_{0}^{\dot{ε}} σ : d \dot{ε} .

$\dot{W}$ is not uniquely defined if multiple deformation mechanisms contribute to the strain rate. However, the creep mechanism will dominate within a zone surrounding a crack tip or crack line, so elastic and plastic contributions to $\dot{W}$ are negligible. The size of that zone depends on the extent of creep relaxation: the zone is initially small but eventually encompasses the entire specimen when steady-state creep is reached. Abaqus/Standard considers only creep in the calculation of $\dot{W}$ . Neglecting elastic and plastic strain rates, the strain energy density for the power law creep model with time hardening form in Abaqus/Standard is

\dot{W} = \frac{n}{n + 1} q \dot{\bar{ε}},

where n is the power law exponent, q is the equivalent von Mises stress, and $\dot{\bar{ε}}$ is the equivalent uniaxial strain rate.

For the hyperbolic-sine law an analytical expression of $\dot{W}$ is not available. For this law $\dot{W}$ is obtained by numerical integration; a five-point Gauss quadrature scheme gives reasonable accuracy in the range of realistic creep strain rates.

The domain integral method is used for $C_{t}$ -integrals as described above for J-integrals.

For user-defined creep laws the strain energy rate density must be defined in user subroutine CREEP.

Domain Dependence

Prior to steady state $C_{t}$ -integral estimates will exhibit domain dependence, even if the finite element mesh is sufficiently refined, because of the assumption of creep dominance within the domain specified. These $C_{t}$ estimates should be extrapolated to zero radius to obtain an improved $C_{t}$ estimate corresponding to a contour shrunk onto the crack tip or crack line (see Ct-integral evaluation).

Stress Intensity Factors

The stress intensity factors $K_{I}$ , $K_{I I}$ , and $K_{I I I}$ are usually used in linear elastic fracture mechanics to characterize the local crack-tip/crack-line stress and displacement fields. They are related to the energy release rate (the J-integral) through

J = \frac{1}{8 π} K^{T} \cdot B^{- 1} \cdot K,

where $K = {⌊ K_{I}, K_{I I}, K_{I I I} ⌋}^{T}$ are the stress intensity factors and $B$ is called the pre-logarithmic energy factor matrix. For homogeneous, isotropic materials $B$ is diagonal, and the above equation simplifies to

J = \frac{1}{\bar{E}} (K_{I}^{2} + K_{I I}^{2}) + \frac{1}{2 G} K_{I I I}^{2},

where $\bar{E} = E$ for plane stress and $\bar{E} = E / (1 - ν^{2})$ for plane strain, axisymmetry, and three dimensions. For an interfacial crack between two dissimilar isotropic materials,

J = \frac{1 - β^{2}}{E^{*}} (K_{I}^{2} + K_{I I}^{2}) + \frac{1}{2 G^{*}} K_{I I I}^{2},

where

\frac{1}{E^{*}} = \frac{1}{2} (\frac{1}{{\bar{E}}_{1}} + \frac{1}{{\bar{E}}_{2}}), \frac{1}{G^{*}} = \frac{1}{2} (\frac{1}{G_{1}} + \frac{1}{G_{2}});

β = \frac{G_{1} (κ_{2} - 1) - G_{2} (κ_{1} - 1)}{G_{1} (κ_{2} + 1) + G_{2} (κ_{1} + 1)};

$κ = 3 - 4 ν$ for plane strain, axisymmetry, and three dimensions; and $κ = (3 - ν) / (1 + ν)$ for plane stress. Unlike their analogues in a homogeneous material, $K_{I}$ and $K_{I I}$ are no longer the pure Mode I and Mode II stress intensity factors for an interfacial crack. They are simply the real and imaginary parts of a complex stress intensity factor.

Although the energy release rate is calculated directly in Abaqus/Standard, it is usually not straightforward to compute stress intensity factors from a known J-integral for mixed-mode problems. Abaqus/Standard provides an interaction integral method to compute the stress intensity factors directly for a crack under mixed-mode loading. This capability is available for linear isotropic and anisotropic materials. The theory is described in detail in Stress intensity factor extraction.

In this case the J-integrals calculated from the stress intensity factors will also be output. These J-integral values may be slightly different from those estimated by requesting the J-integral directly, due to the different algorithms used for the calculations.

Domain Dependence

The stress intensity factors have the same domain dependence features as the J-integral.

Including the Effect of a Residual Stress Field on Stress Intensity Factor Evaluation

An additional term is included to account for the residual stress field when calculating the stress intensity factors, as described in Including the Effect of a Residual Stress Field on J-Integral Evaluation.

Crack Propagation Direction

For homogeneous, isotropic elastic materials the direction of cracking initiation can be calculated using one of the following three criteria: the maximum tangential stress criterion, the maximum energy release rate criterion, or the $K_{I I} = 0$ criterion. $K_{I I I}$ is not taken into account in any of these criteria.

Maximum Tangential Stress Criterion

Using either the condition $\partial σ_{θ θ} / \partial θ = 0$ or $τ_{r θ} = 0$ (where r and $θ$ are polar coordinates centered at the crack tip in a plane orthogonal to the crack line), we can obtain

\hat{θ} = \cos^{- 1} (\frac{3 K_{I I}^{2} + \sqrt{K_{I}^{4} + 8 K_{I}^{2} K_{I I}^{2}}}{K_{I}^{2} + 9 K_{I I}^{2}}),

where the crack propagation angle $\hat{θ}$ is measured with respect to the crack plane and $\hat{θ} = 0$ represents the crack propagation in the “straight-ahead” direction. $\hat{θ} < 0$ if $K_{I I} > 0$ while $\hat{θ} > 0$ if $K_{I I} < 0$ The crack propagation angle is measured from $q$ to $n$ ; that is, it is measured about the direction $- t$ , or counterclockwise measured from $q$ in Figure 1.

The crack propagation angle $\hat{θ}$ will be output.

Maximum Energy Release Rate Criterion

This criterion postulates that a crack initially propagates in the direction that maximizes the energy release rate.

The crack propagation angle $\hat{θ}$ will be output.

K_II = 0 Criterion

This criterion assumes that a crack initially propagates in the direction that makes $K_{I I} = 0$ .

The crack propagation angle $\hat{θ}$ will be output.

T-Stress

The T-stress component represents a stress parallel to the crack faces at the crack tip. Its magnitude can alter not only the size and shape of the plastic zone but also the stress triaxiality ahead of the crack tip. It is, therefore, a useful indicator of whether measures of the strength of the crack-tip singularity (such as the J-integral or the stress intensity factors) are useful in characterizing a crack under a particular loading. In a linear elastic analysis the T-stress should be calculated using loads equal to the loads in the elastic-plastic analysis. See T-stress extraction for more information.

Domain Dependence

In general, the T-stress has larger domain dependence or contour dependence than the J-integral and the stress intensity factors. Numerical tests suggest that the estimates from the first two rings of elements abutting the crack tip or crack line generally do not provide accurate results. Sufficient contours extending from the crack tip or crack line should be chosen so that the T-stress can be determined to be independent of the number of contours, within engineering accuracy. Particularly for axisymmetric models, the closer the crack tip is to the symmetry axis, the more refined the mesh in the domain should be to achieve path independence of the contour integral.

Including the Effect of a Residual Stress Field on T-Stress Evaluation

An additional term is included to account for the residual stress field when calculating the T-stress, as described in Including the Effect of a Residual Stress Field on J-Integral Evaluation.

Defining the Data Required for a Contour Integral with the Conventional Finite Element Method

To request contour integral output with the conventional finite element method, you must define the crack front and specify the virtual crack extension direction.

Defining the Crack Front

You must specify the crack front; that is, the region that defines the first contour. Abaqus/Standard uses this region and one layer of elements surrounding it to compute the first contour integral. An additional layer of elements is used to compute each subsequent contour.

The crack front can be equivalent to the crack tip in two dimensions or the crack line in three dimensions; or it can be a larger region surrounding the crack tip or crack line, in which case it must include the crack tip or crack line.

If blunted crack tips are modeled, the crack front should include all the nodes going from one crack face to the other that would collapse onto the crack tip if the radius of the blunted tip were reduced to zero. Otherwise, the contour integral value will depend on the path until the contour region reaches the parallel crack faces.

Defining the Crack Tip or Crack Line

By default, Abaqus/Standard defines the crack tip as the first node specified for the crack front and the crack line as the sequence of first nodes specified for the crack front. The first node is the node with the smallest node number, unless the node set is generated as unsorted. Alternatively, you can specify the crack-tip node or crack-line nodes directly. This specification plays a critical role for a three-dimensional crack with a blunt crack tip.

Defining a Closed-Loop Crack Line

Sometimes a crack line may form a closed loop (for example, when modeling a full penny-shaped crack without invoking symmetry conditions). In such cases the finite element mesh in the crack-tip region can be created with or without seams; that is, linear constraint equations (Linear Constraint Equations) or multi-point constraints (General Multi-Point Constraints) may or may not be used to tie two layers of nodes together.

If a crack line forms a closed loop, the starting node set of the crack front can be chosen arbitrarily and the other node sets defining the crack front must go around the crack front sequentially. The last node set defining the crack front must be the same as the first node set. If a closed loop is formed by creating coincident nodes that are then tied together by linear constraint equations and multi-point constraints, the node sets must be specified in order starting from one of the node sets involved in the constraint equation or multi-point constraint and terminating with the other node set.

Specifying the Virtual Crack Extension Direction

You must specify the direction of virtual crack extension at each crack tip in two dimensions or at each node along the crack line in three dimensions by specifying either the normal to the crack plane, $n$ , or the virtual crack extension direction, $q$ .

If the virtual crack extension direction is specified to point into the material (parallel to the crack faces), the J-integral values calculated will be positive. Negative J-integral values are obtained when the virtual crack extension direction is specified in the opposite direction.

Specifying the Normal to the Crack Plane

The virtual crack extension direction can be defined by specifying the normal, $n$ , to the crack plane. In this case Abaqus/Standard will calculate a virtual crack extension direction, $q$ , that is orthogonal to the crack front tangent, $t$ , and the normal, $n$ . As shown in Figure 1, $q = t \times n$ for a three-dimensional crack; for a two-dimensional crack, we simply have $q_{x} = - n_{y}$ and $q_{y} = n_{x}$ . Specifying the normal implies that the crack plane is flat since only one value of $n$ can be given per contour integral.

Specifying the Virtual Crack Extension Direction

Alternatively, the virtual crack extension direction, $q$ , can be specified directly. In three dimensions the virtual crack extension direction, $q$ , will be corrected to be orthogonal to any normal defined at a node or in other cases to the tangent to the crack line itself. The tangent, $t$ , to the crack line at a particular point is obtained by parabolic interpolation through the crack front for which the virtual crack extension vector is defined and the nearest node sets on either side of this region. Abaqus/Standard will normalize the virtual crack extension direction, $q$ .

Defining Surface Normals

In a case where the crack front intersects the external surface of a three-dimensional solid, where there is a surface of material discontinuity in the model, or where the crack is in a curved shell, the virtual crack extension direction, $q$ , must lie in the plane of the surface for accurate contour integral evaluation. Surface normals should be specified at all nodes that lie on such surfaces within the contours requested for this purpose (these nodes are printed out under the “Contour Integral” information in the data file). For shell element models the normals can be specified with the nodal coordinates if the normals calculated by Abaqus/Standard are not adequate. For solid element models the normals can be specified either directly (see Normal Definitions at Nodes and A plate with a part-through crack: elastic line spring modeling) or using the nodal coordinates (the fourth–sixth coordinates).

If surface normals are not specified for the nodes on the crack surfaces and the external surfaces at the ends of a crack line, Abaqus/Standard can calculate the normals for these nodes to correct any inadequate virtual crack extension directions, $q$ except for a model with tetrahedral elements. For a model with tetrahedral elements, you must always specify the surface normals to improve the accuracy of contour integral evaluation at the free surface. For large models, requesting that Abaqus/Standard calculate the surface normals on free or crack surfaces can increase the preprocessing time.

Defining the Data Required for a Contour Integral with XFEM

If you are using XFEM to evaluate the contour integral, both the crack front and the virtual crack extension direction are determined by Abaqus/Standard.

Symmetry with the Conventional Finite Element Method

If the crack is defined on a symmetry plane, only half the structure needs to be modeled. The change in potential energy calculated from the virtual crack front advance is doubled to compute the correct contour integral values.

Constructing a Fracture Mechanics Mesh for Small-Strain Analysis with the Conventional Finite Element Method

Sharp cracks (where the crack faces lie on top of one another in the undeformed configuration) are usually modeled using small-strain assumptions. Focused meshes, as shown in Figure 1, should typically be used for small-strain fracture mechanics evaluations. However, for a sharp crack the strain field becomes singular at the crack tip. This result is obviously an approximation to the physics; however, the large-strain zone is very localized, and most fracture mechanics problems can be solved satisfactorily using only small-strain analysis.

The crack-tip strain singularity depends on the material model used. Linear elasticity, perfect plasticity, and power-law hardening are commonly used in fracture mechanics analysis. Power-law hardening has the form

\frac{\bar{ε}}{ε_{0}} = α {(\frac{\bar{σ}}{σ_{0}})}^{n},

where $\bar{ε}$ is the equivalent total strain, $ε_{0}$ is a reference strain, $\bar{σ}$ is the von Mises stress, $σ_{0}$ is the initial yield stress, n is the power-law hardening exponent (typically in the range of 3 to 8; $n > 10$ is very close to perfect plasticity for large $ε$ ), and $α$ is a material constant (typically in the range 0.5 to 1.0).

Results for pure power-law nonlinear elastic materials in a body under traction loading are proportional to the load to some power. Therefore, the fracture parameters for one geometry under a particular load can be scaled to any other load of the same distribution but different magnitude.

If the loading is proportional (the direction of the stress increase in stress space is approximately constant) and monotonically increasing, power-law hardening deformation plasticity and incremental plasticity are essentially equivalent. However, deformation plasticity is a nonlinear elastic material for which more analytical results are available. Abaqus uses the Ramberg-Osgood form of deformation plasticity (see Deformation Plasticity); this model is not a pure power law model, which must be considered.

Creating the Singularity

In most cases the singularity at the crack tip should be considered in small-strain analysis (when geometric nonlinearities are ignored). Including the singularity often improves the accuracy of the J-integral, the stress intensity factors, and the stress and strain calculations because the stresses and strains in the region close to the crack tip are more accurate. If r is the distance from the crack tip, the strain singularity in small-strain analysis is

\begin{matrix} ε \propto r^{- 1 / 2} & for
  linear
  elasticity, \\ ε \propto r^{- 1} & for
  perfect plasticity,
  and \\ ε \propto r^{- \frac{n}{n + 1}} & for
  power-law
  hardening. \end{matrix}

Modeling the Crack-Tip Singularity in Two Dimensions

The square root and $1 / r$ singularity can be built into a finite element mesh using standard elements. The crack tip is modeled with a ring of collapsed quadrilateral elements, as shown in Figure 2.

To obtain a mesh singularity, generally second-order elements are used and the elements are collapsed as follows:

Collapse one side of an 8-node isoparametric element (CPE8R, for example) so that all three nodes—a, b, and c—have the same geometric location (on the crack tip).
Move the midside nodes on the sides connected to the crack tip to the 1/4 point nearest the crack tip. You can create “quarter point” spacing with second-order isoparametric elements when you generate nodes for a region of a mesh; see Creating Quarter-Point Spacing.

This procedure will create the strain singularity

ε \to \frac{A}{r} + \frac{B}{r^{1 / 2}} as r \to 0.

The $r^{- \frac{n}{n + 1}}$ singularity cannot be created using Abaqus elements, but the combination of the $r^{- 1}$ and $r^{- 1 / 2}$ terms can provide a reasonable approximation for $r^{- \frac{n}{n + 1}}$ .

If 4-node isoparametric elements (for example, CPE4R) are used, one side of the element is collapsed, and the two coincident nodes are free to displace independently, a $1 / r$ singularity is created.

If the crack region is meshed with linear elements, the position specified for the midside nodes is ignored.

Creating a Square Root Singularity

If nodes a, b, and c are constrained to move together, $A = 0$ and the strains and stresses are square root singular (suitable for linear elasticity).

Creating a 1/r Singularity

If the midside nodes remain at the midside points rather than being moved to the 1/4 points and nodes a, b, and c are allowed to move independently, only the $1 / r$ singularity in strain is created (suitable for perfect plasticity).

Creating a Combined Square Root and 1/r Singularity

If the midside nodes are moved to the 1/4 points but nodes a, b, and c are allowed to move independently, the singularity created is a combination of the square root and $1 / r$ singularities. This combination is usually best for a power-law hardening material. However, since the $1 / r$ singularity dominates, moving the midside nodes to the 1/4 points gives only slightly better results than if the nodes are left at the midside points. Since creating a mesh with the midside nodes moved to the quarter points can be difficult, it is often best to simply use the $1 / r$ singularity.

Modeling the Crack-Tip Singularity in Three Dimensions

To create singular fields, 20-node bricks and 27-node bricks can be used with a collapsed face (see Figure 3).

The planes of the three-dimensional elements perpendicular to the crack line should be planar for the best accuracy. If they are not planar, the element Jacobian may become negative at some integration points when the midside nodes are moved to the 1/4 points. To correct this problem, move the midside nodes slightly away from the 1/4 points toward the midpoint position (the distance moved is not critical).

Creating a Square Root Singularity

To obtain a square root singularity, constrain the nodes on the collapsed face of the edge planes to move together and move the nodes to the 1/4 points.

If the nodes at the midplane of a collapsed 20-node brick are constrained to move together, $ε \propto r$ ; therefore, the singularity is not the same on the midplane as on an edge plane. This difference causes local oscillations in the solution about the crack tip along the crack line, although normally the oscillations are not significant.

If all midface nodes and the centroid node are included in a 27-node brick and the midside and midface nodes are moved to the 1/4 points closest to the crack line, the oscillation in the local stress and strain fields can be reduced.

Creating a 1/r Singularity

To obtain a $1 / r$ singularity, allow the three nodes on the collapsed face to displace independently and keep the midside nodes at the midpoints.

Creating a Combined Square Root and 1/r Singularity

To obtain a combined square root and $1 / r$ singularity, allow the nodes on the collapsed face to displace independently and move the midside nodes to the 1/4 points. As in the two-dimensional case, if it is difficult to create the mesh with the nodes moved to the 1/4 points, simply use the $1 / r$ singularity.

Mesh Refinement

The size of the crack-tip elements influences the accuracy of the solutions: the smaller the radial dimension of the elements from the crack tip, the better the stress, strain, etc. results will be and, therefore, the better the contour integral calculations will be.

The angular strain dependence is not modeled with the singular elements. Reasonable results are obtained if typical elements around the crack tip subtend angles in the range of 10° (accurate) to 22.5° (moderately accurate).

Since the crack tip causes a stress concentration, the stress and strain gradients are large as the crack tip is approached. Path dependence in the evaluation of the J-integral may be an indication that the mesh is not sufficiently refined, but path independence does not prove mesh convergence. The finite element mesh must be refined in the vicinity of the crack to get accurate stresses and strains; however, accurate J-integral results can frequently be obtained even with a relatively coarse mesh.

In many cases if sufficiently fine meshes are used, accurate contour integral values can be obtained without using singular elements.

Modeling the Crack-Tip Region in Shells

Focused meshes can be used, but not all of the three-dimensional shell elements in Abaqus/Standard can be collapsed. Elements S8R and S8RT cannot be degenerated into triangles; element types S4, S4R, S4R5, S8R5, and S9R5 can.

The quarter-point technique (moving the midside nodes to the quarter points to give a $1 / \sqrt{r}$ singularity for elastic fracture mechanics applications) can be used with S8R5 and S9R5 elements but not with S8R(T) elements. When the quarter-point technique is used with S9R5 elements, the midface node should be moved to the quarter-point position along with the two midside nodes.

If S8R(T) elements are used, a keyhole should be introduced at the crack tip.

Flaws lying in the plane through the thickness of a shell can be modeled using line spring elements; see Line Spring Elements for Modeling Part-through Cracks in Shells. In many cases line spring elements provide accurate J-integral and stress intensity values, but these elements are limited to modeling small strain and rotations. Limited modeling of plasticity is also allowed with line springs.

Constructing a Fracture Mechanics Mesh for Finite-Strain Analysis with the Conventional Finite Element Method

In large-strain analysis (when geometric nonlinearities are included) singular elements should not normally be used. The mesh must be sufficiently refined to model the very high strain gradients around the crack tip if details in this region are required. Even if only the J-integral is required, the deformation around the crack tip may dominate the solution and the crack-tip region will have to be modeled with sufficient detail to avoid numerical problems.

Physically, the crack tip is not perfectly sharp. Therefore, it is normally modeled as a blunted notch with a radius of $ \sim 10^{- 3} r_{p}$ , where $r_{p}$ is a characteristic dimension of the plastic zone ahead of the crack tip. The notch must be small enough that, at the loads of interest, the deformed shape of the notch no longer depends on the original geometry. Typically, the notch must blunt out to more than four times its original radius for the deformed shape to be independent of the original geometry. The size of the elements around the notch should be about 1/10 the notch-tip radius to obtain accurate results.

If a crack is modeled as sharp, the finite elements near the crack tip may not be able to approximate the high gradients, resulting in convergence problems. The stress and strain results around the crack tip will probably be inaccurate even if convergence is achieved. However, if the solution converges, the contour integral results should be reasonably accurate. The convergence difficulties will probably be greater in three dimensions than in two dimensions.

In situations involving finite rotations but small strains, such as bending of slender structures, a small “keyhole” around the crack tip should be modeled. If the hole is small, the results will not be affected significantly and problems in dealing with the singular strains at the crack tip will be avoided.

Using Constraints with the Conventional Finite Element Method

General multi-point constraints and linear constraint equations (About Kinematic Constraints) should not be used on nodes in the mesh regions where contour integrals are calculated unless the nodes involved in the constraint are located at the same point. The nodes at the crack tip of a focused mesh can be tied together using multi-point constraints without adversely affecting the contour integral calculations. Tying these nodes will change the singularity at the crack tip, but path independence of the contour integral will be maintained. In addition, path independence of the contour integrals will not be affected if two faces of a model are joined using MPC type TIE or a linear constraint equation, provided that all nodes of the two faces are coincident. Using multi-point constraints for mesh refinement or for applying symmetry/antisymmetry boundary conditions within the contour integral region will result in path dependence of the contour integrals. No warning or error messages are provided if this rule is violated.

Procedures

You can request contour integrals in fracture mechanics problems that were modeled using the following procedures:

static (Static Stress Analysis) with both XFEM and the conventional finite element methods;
quasi-static (Quasi-Static Analysis) with the conventional finite element method only;
steady-state transport (Steady-State Transport Analysis) with the conventional finite element method only;
coupled thermal-stress procedures (Fully Coupled Thermal-Stress Analysis) with the conventional finite element method only; and
crack propagation (Crack Propagation Analysis) with the conventional finite element method only.

Contour integrals can be requested only in general analysis steps: they are not calculated in linear perturbation analyses (General and Perturbation Procedures).

A crack analysis with pressure applied on the crack surfaces might give inaccurate contour integral values if geometric nonlinearity is included in a step. Similarly, the calculated results of the stress intensity factors and T-stress might not be accurate if geometric nonlinearity is included in a step.

Loads

Contour integral calculations include the following distributed load types:

thermal loads;
distributed loads, including crack face pressure and traction loads on continuum elements as well as those applied using user subroutine DLOAD and UTRACLOAD;
distributed loads, including surface traction loads and crack face edge loads on shell elements as well as those applied using user subroutine UTRACLOAD;
uniform and nonuniform body forces; and
centrifugal loads on continuum and shell elements.

Contributions to the contour integral due to concentrated loads in the domain are not included; instead, the mesh must be modified to include a small element and a distributed load must be applied to this element.

Contributions due to contact forces are not included.

Material Options

J-integral calculations are valid for linear elastic, nonlinear elastic, and elastic-plastic materials. Plastic behavior can be modeled as nonlinear elastic (Deformation Plasticity), but the results are generally best if the material is modeled by incremental plasticity and is subject to proportional, monotonic traction loading.

If unloading has taken place in the plastic zone around the crack tip, the J-integral will not be valid except in very limited cases.

The $C_{t}$ -integral is valid for problems involving creep (Rate-Dependent Plasticity: Creep and Swelling).

The stress intensity factor calculation is valid for cracks in homogeneous, linear elastic materials. It is also valid for an interfacial crack between two different isotropic linear elastic materials. It is not valid for any other types of materials, including user-defined materials.

The crack propagation direction is valid only for homogeneous, isotropic linear elastic materials.

The T-stress is valid only for homogeneous, isotropic linear elastic materials. Although the T-stress is calculated using the linear elastic material properties of the body with a crack, it is usually used with the J-integral calculated using the elastic-plastic material properties of the body (see T-stress extraction).

If there is material discontinuity, the normal to the material discontinuity line must be specified for all nodes on the material discontinuity that will lie in a contour integral domain. The normal can be specified by defining user-specified normals (see Normal Definitions at Nodes) for the elements on both sides of the discontinuity or by using nodal normal coordinates for the nodes on the discontinuity. Contour integral calculations cannot be performed for a crack with a material discontinuity line passing through its tip (except for an interfacial crack between two different materials). Therefore, you should be careful when specifying a normal that is not perpendicular to the virtual crack extension direction, $q$ , for the nodes at the crack tip.

Elements

When used with XFEM, the contour integral can be evaluated only in first-order or second-order tetrahedral and first-order brick elements. The following paragraphs apply only to the conventional finite element method.

The contour integral evaluation capability in Abaqus/Standard assumes that the elements that lie within the domain used for the calculations are quadrilateral elements in two-dimensional or shell models or are bricks or second-order tetrahedral elements in continuum three-dimensional models. You should not use triangles or wedges in the mesh that is included in the contour integral regions. The elements within the contour domain should be of the same type.

In shell structures the contour integrals calculated by Abaqus/Standard are contour independent only if the deformation mode around the crack tip is primarily membrane. If there are significant bending or transverse shear effects in the domain, the contour integrals may not be contour independent and contour integral values should be obtained directly from the displacements and/or the stresses.

Generalized plane strain elements, generalized axisymmetric elements with twist, asymmetric-axisymmetric elements, membrane elements, and cylindrical elements should not be used in the contour integral regions.

The contribution of rebar is included only in the calculations of the J-integral and the $C_{t}$ -integral for shell elements defined with a shell section integrated during the analysis (see Using a Shell Section Integrated during the Analysis to Define the Section Behavior).

Output

The domain associated with each contour is calculated automatically. The nodes belonging to each domain can be printed in the data file; see Controlling the Amount of analysis input file processor Information Written to the Data File. If you are using the conventional contour integral method, for each domain Abaqus/Standard creates a new node set in the output database to include these nodes. In addition, new node sets are created in the output database for nodes on crack surfaces and on free surfaces whose nodal normals are calculated by Abaqus/Standard.

Contour integrals cannot be recovered from the restart file as described in About Output.

You should not request element output extrapolated to the nodes (Element Output) for second-order elements with one collapsed side in two dimensions or one collapsed face in three dimensions.

Default Contour Integral Output

By default, the contour integral values are written to the data file and to the output database file. The following naming convention is used for contour integrals written to the output database:

integral-type: abbrev-integral-type at history-output-request-name_crack-name_internal-crack-tip-node-set-name__Contour_contour-number

where integral-type can be

Crack propagation direction (Cpd)
J-integral (J)
J-integral estimated from Ks (JKs)
Stress intensity factor K1 (K1)
Stress intensity factor K2 (K2)
T-stress (T)

For example,

J-integral: J at JINT_CRACK_CRACKTIP-1__Contour_1

Writing the Contour Integrals to the Results File

You can choose to write the contour integral values to the results file in addition to the data file.

Controlling the Output Frequency

You can control the output frequency, in increments, of contour integrals. By default, the crack-tip location and associated quantities will be printed every increment. Specify an output frequency of 0 to suppress contour integral output.

The output frequency for contour integral output to the output database is controlled by the larger of the frequency values specified for history output to the output database or for contour integral output. If you specify an output frequency of 0 for the history output to the output database, contour integral values will not be written to the output database.

Requesting Field Output of the Contour Integral

For the conventional contour integral method, if $n$ contours are specified, you can request the averaged value of contour integrals over $\frac{n}{2}$ contours (starting from contour number $\frac{n}{2}$ ) as nodal field output written to the output database.

Nodal field output is not available for the XFEM-based contour integral method.

The following nodal field output variables are available:

J: Averaged value of the contour integral.
K: Averaged values of the stress intensity factors.
TSTRESS: Averaged value of the T-stress.

References

Lei, Y., “J-Integral Evaluation for Cases Involving Non-proportional Stressing,” Engineering Fracture Mechanics, vol. 72, pp. 577–596, 2005.
Lei, Y., N. P. O'Dowd, and G. A. Webster, “Fracture Mechanics Analysis of a Crack in a Residual Stress Field,” International Journal of Fracture, vol. 106, pp. 195–216, 2000.

Contour Integral Evaluation