Fracture Integrals


A finite element calculation in which a crack is represented in the mesh geometry can be used to calculate the displacement and stress fields in the vicinity of the crack. One straightforward way to evaluate the stress intensity from a finite element solution is through the -integral (Rice, 1968), which provides the mechanical energy release rate. If the crack is subjected to pure mode- loading, can be calculated from using the following relationship: (1) where is the Young's modulus and is the Poisson's ratio. The -integral was originally formulated as an integral over a closed curve around the crack tip in 2D. It can alternatively be expressed as an integral over a surface in 2D or a volume in 3D using the method of Li et al. (1985), which is more convenient for implementation within a finite element code.

Following the terminology of Healy et al. (2012), the integrated value of the -integral over a segment of a crack front, represented as , can be expressed as a summation of four terms: (2) where the individual terms are defined as: (3)



(6) In these equations, is the undeformed volume, is the combined area of the two crack faces, is the 1st Piola-Kirchoff stress tensor, is the displacement vector, is the undeformed coordinate vector, is the stress-work density, is the kinetic energy density, is time, is the material density, and is the vector of tractions applied to the crack face.

The vector field is a vector field that is oriented in the crack normal direction. This field has a magnitude of 1 between the crack tip and the inner radius of the ring over which the integral is to be performed, and drops from 1 to 0 between the inner and outer radius of that ring, and has a value of 0 elsewhere. In 3D, is evaluated by calculating the integral over a segment of the crack front. A separate field is formed for each segment along the crack front, and for each ring over which the integral is to be evaluated. Each of those fields is multiplied by a weighting function that varies from a value of 0 at the ends to a finite value in the middle of the segment. The value of at each point on the curve is evaluated by dividing by the integrated value of the weighting function over the segment containing that point.

The first term in , , represents the effects of strain energy in homogeneous material in the absence of thermal strains or inertial effects. The second term, , accounts for the effects of material inhomogenieties and gradients of thermal strain. For small strains, this term can be represented as: (7) where , and are the Cauchy stress, thermal conductivity, and difference between the current temperature and a reference temperature, respectively.

The third term in , , accounts for inertial effects in a dynamic analysis, and the fourth term, accounts for the effects of surface tractions on the crack face.

The MOOSE implementation of the -integral calculator can be used for arbitrary curved 3D crack fronts, and includes the and terms to account for the effects of quasistatic mechanically and thermally induced strains. The last two terms have not been implemented. For quasistatic loading conditions, there are no inertial effects, so . would be nonzero for pressurized cracks, and should be implemented in the future to consider such cases.

Pointwise values are calculated from the -integral over a crack front segment using the expression (8) To use the -integral capability in MOOSE, a user specifies a set of nodes along the crack front, information used to calculate the crack normal directions along the crack front, and the inner and outer radii of the set of rings over which the integral is to be performed. The code automatically defines the full set of functions for each point along the crack front, and outputs either the value of or at each of those points. In addition, the user can ask for any other variable in the model to be output at points corresponding to those where is evaluated along the crack front.

Interaction Integral

The interaction integral method is based on the -integral and makes it possible to evaluate mixed-mode stress intensity factors , and , as well as the T-stress, in the vicinity of three-dimensional cracks. The formulation relies on superimposing Williams' solution for stress and displacement around a crack (in this context called 'auxiliary fields') and the computed finite element stress and displacement fields (called 'actual fields'). The total superimposed can be separated into three parts: the of the actual fields, the of the auxiliary fields, and an interaction part containing the terms with both actual and auxiliary field quantities. The last part is called the interaction integral and for a fairly straight crack without body forces, thermal loading or crack face tractions, the interaction integral over a crack front segment can be written (Walters et al., 2005): (9) where is the stress, is the displacement, and is identical to the -functions used for -integrals. This is the formulation of the interaction integral used currently in MOOSE. It is not recommended to use the interaction integral method in analyses where body forces or crack face tractions are important, as the code does not include the necessary correction terms in the integral.

In the same way as for the -integral, the pointwise value at location is obtained from using: (10) Next, we relate the interaction integral to stress intensity factors. By writing , with actual and auxiliary fields superimposed, in terms of the mixed-mode stress intensity factors (11) the interaction integral part evaluates to (12) To obtain individual stress intensity factors, the interaction integral is evaluated with different auxiliary fields. For instance, by choosing and and computing the volume integral over the actual and auxiliary fields, can be solved for.

If all three interaction integral stress intensity factors are computed, there is an option to output an equivalent stress intensity factor computed by: (13) The -stress is the first second-order parameter in Williams' expansion of stress at a crack tip and is a constant stress parallel to the crack (Larsson and Carlsson, 1973) (Rice, 1974). Contrary to , -stress depends on geometry and size and can give a more accurate description of the stresses and strains around a crack tip than alone. The -stress characterizes the crack-tip constraint and a negative -stress is associated with loss of constraint and a higher fracture toughness than would be predicted from a one-parameter description of the load on the crack. -stresses can be calculated with the interaction integral methodology using appropriate auxiliary fields (Nakamura and Parks, 1992). The current implementation of the -stress computation is valid for two-dimensional and three-dimensional cracks under Mode I loading.


The MOOSE implementation of the capability to compute these fracture integrals is provided using a variety of MOOSE objects, which are quite complex to define manually, especially for 3D simulations. For this reason, the to compute -integrals or interaction integrals, the DomainIntegral Action should be used to set up all of the required objects.

  integrals = JIntegral
  boundary = 800
  crack_direction_method = CrackDirectionVector
  crack_direction_vector = '1 0 0'
  radius_inner = '4.0 5.5'
  radius_outer = '5.5 7.0'
  output_variable = 'disp_x'
  output_q = false
  incremental = true
  1. Toshio Nakamura and David M. Parks. Determination of elastic t-stress along three-dimensional crack fronts using an interaction integral. International Journal of Solids and Structures, 29(13):1597–1611, January 1992. doi:10.1016/0020-7683(92)90011-H.[BibTeX]
  2. Brian Healy, Arne Gullerud, Kyle Koppenhoefer, Arun Roy, Sushovan RoyChowdhury, Matt Walters, Barron Bichon, Kristine Cochran, Adam Carlyle, James Sobotka, Mark Messner, and Robert Dodds. Warp3d-release 17.3.1. Technical Report UILU-ENG-95-2012, University of Illinois at Urbana-Champaign, 2012.[BibTeX]
  3. J. R. Rice. A path independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics, 35(2):379, 1968. doi:10.1115/1.3601206.[BibTeX]
  4. J.R. Rice. Limitations to the small scale yielding approximation for crack tip plasticity. Journal of the Mechanics and Physics of Solids, 22(1):17–26, January 1974.[BibTeX]
  5. Matthew C. Walters, Glaucio H. Paulino, and Robert H. Dodds. Interaction integral procedures for 3-D curved cracks including surface tractions. Engineering Fracture Mechanics, 72(11):1635–1663, July 2005.[BibTeX]
  6. S.G. Larsson and A.J. Carlsson. Influence of non-singular stress terms and specimen geometry on small-scale yielding at crack tips in elastic-plastic materials. Journal of the Mechanics and Physics of Solids, 21(4):263–277, July 1973.[BibTeX]
  7. F. Z. Li, C. F. Shih, and A. Needleman. A comparison of methods for calculating energy release rates. Engineering Fracture Mechanics, 21(2):405–421, 1985. doi:10.1016/0013-7944(85)90029-3.[BibTeX]