Energy release rate (fracture mechanics)

In fracture mechanics, the energy release rate, formula_1, is the rate at which energy is transformed as a material undergoes fracture. Mathematically, the energy release rate is expressed as the decrease in total potential energy per increase in fracture surface area, and is thus expressed in terms of energy per unit area. Various energy balances can be constructed relating the energy released during fracture to the energy of the resulting new surface, as well as other dissipative processes such as plasticity and heat generation. The energy release rate is central to the field of fracture mechanics when solving problems and estimating material properties related to fracture and fatigue.
Definition.
The energy release rate formula_1 is defined as the instantaneous loss of total potential energy formula_3 per unit crack growth area formula_4,
where the total potential energy is written in terms of the total strain energy formula_6, surface traction formula_7, displacement formula_8, and body force formula_9 by
The first integral is over the surface formula_11 of the material, and the second over its volume formula_12.
The figure on the right shows the plot of an external force formula_13 vs. the load-point displacement formula_14, in which the area under the curve is the strain energy. The white area between the curve and the formula_13-axis is referred to as the complementary energy. In the case of a linearly-elastic material, formula_16 is a straight line and the strain energy is equal to the complementary energy.
Prescribed displacement.
In the case of prescribed displacement, the strain energy can be expressed in terms of the specified displacement and the crack surface formula_17, and the change in this strain energy is only affected by the change in fracture surface area: formula_18. Correspondingly, the energy release rate in this case is expressed as 
Here is where one can accurately refer to formula_1 as the strain energy release rate.
Prescribed loads.
When the load is prescribed instead of the displacement, the strain energy needs to be modified as formula_21. The energy release rate is then computed as 
If the material is linearly-elastic, then formula_23 and one may instead write
G in two-dimensional cases.
In the cases of two-dimensional problems, the change in crack growth area is simply the change in crack length times the thickness of the specimen. Namely, formula_25. Therefore, the equation for computing formula_1 can be modified for the 2D case: 
One can refer to the example calculations embedded in the next section for further information. Sometimes, the strain energy is written using formula_30, an energy-per-unit thickness. This gives
Relation to stress intensity factors.
The energy release rate is directly related to the stress intensity factor associated with a given two-dimensional loading mode (Mode-I, Mode-II, or Mode-III) when the crack grows straight ahead. This is applicable to cracks under plane stress, plane strain, and antiplane shear.
For Mode-I, the energy release rate formula_1 rate is related to the Mode-I stress intensity factor formula_35 for a linearly-elastic material by
where formula_37 is related to Young's modulus formula_38 and Poisson's ratio formula_39 depending on whether the material is under plane stress or plane strain:
For Mode-II, the energy release rate is similarly written as
For Mode-III (antiplane shear), the energy release rate now is a function of the shear modulus formula_42,
For an arbitrary combination of all loading modes, these linear elastic solutions may be superposed as
Relation to fracture toughness.
Crack growth is initiated when the energy release rate over comes a critical value formula_45, which is a material property,
Under Mode-I loading, the critical energy release rate formula_45 is then related to the Mode-I fracture toughness formula_48, another material property, by
Calculating "G".
There are a variety of methods available for calculating the energy release rate given material properties, specimen geometry, and loading conditions. Some are dependent on certain criteria being satisfied, such as the material being entirely elastic or even linearly-elastic, and/or that the crack must grow straight ahead. The only method presented that works arbitrarily is that using the total potential energy. If two methods are both applicable, they should yield identical energy release rates.
Total potential energy.
The only method to calculate formula_1 for arbitrary conditions is to calculate the total potential energy and differentiate it with respect to the crack surface area. This is typically done by:
all in terms of the crack surface area.
Compliance method.
If the material is linearly-elastic, the computation of its energy release rate can be much simplified. In this case, the Load vs. Load-point Displacement curve is linear with a positive slope, and the displacement per unit force applied is defined as the compliance, formula_51 
The corresponding strain energy formula_6 (area under the curve) is equal to 
Using the compliance method, one can show that the energy release rate for both cases of prescribed load and displacement come out to be 
Multiple specimen methods for nonlinear materials.
In the case of prescribed displacement, holding the crack length fixed, the energy release rate can be computed by 
while in the case of prescribed load,
As one can see that, in both cases, the energy release rate formula_1 times the change in surface formula_59 returns the area between curves, which indicates the energy dissipated for the new surface area as illustrated in the right figure 
Crack closure integral.
Since the energy release rate is defined as the negative derivative of the total potential energy with respect to crack surface growth, the energy release rate may be written as the difference between the potential energy before and after the crack grows. After some careful derivation, this leads one to the crack closure integral 
where formula_62 is the new fracture surface area, formula_63 are the components of the traction released on the top fracture surface as the crack grows, formula_64 are the components of the crack opening displacement (the difference in displacement increments between the top and bottom crack surfaces), and the integral is over the surface of the material formula_65.
The crack closure integral is valid only for elastic materials, but is still valid for cracks that grow in any direction. Nevertheless, for a two-dimensional crack that does indeed grow straight ahead, the crack closure integral simplifies to 
where formula_67 is the new crack length, and the displacement components are written as a function of the polar coordinates formula_68 and formula_69.
J"-"integral.
In certain situations, the energy release rate formula_1 can be calculated using the 
J-integral, i.e. formula_71, using 
where formula_73 is the elastic strain energy density, formula_74is the formula_75component of the unit vector normal to formula_76, the curve used for the line integral, formula_77are the components of the traction vector formula_78, where formula_79 is the stress tensor, and formula_80are the components of the displacement vector.
This integral is zero over a simple closed path and is path independent, allowing any simple path starting and ending on the crack faces to be used to calculate formula_81.
In order to equate the energy release rate to the J-integral, formula_71, the following conditions must be met:
The J-integral may be calculated with these conditions violated, but then formula_84. When they are not violated, one can then relate the energy release rate and the J-integral to the elastic moduli and the stress intensity factors using 
Computational methods in fracture mechanics.
A handful of methods exist for calculating formula_1 with finite elements. Although a direct calculation of the J-integral is possible (using the strains and stresses outputted by FEA), approximate approaches for some type of crack growth exist and provide reasonable accuracy with straightforward calculations. This section will elaborate on some relatively simple methods for fracture analysis utilizing numerical simulations.
Nodal release method.
If the crack is growing straight, the energy release rate can be decomposed as a sum of 3 terms formula_87 associated with the energy in each 3 modes. As a result, the Nodal Release method (NR) can be used to determine formula_87 from FEA results. The energy release rate is calculated at the nodes of the finite element mesh for the crack at an initial length and extended by a small distance formula_67. First, we calculate the displacement variation at the node of interest formula_90(before and after the crack tip node is released). Secondly, we keep track of the nodal force formula_91 outputted by FEA. Finally, we can find each components of formula_1 using the following formulas:
formula_93formula_94formula_95Where formula_67 is the width of the element bounding the crack tip. The accuracy of the method highly depends on the mesh refinement, both because the displacement and forces depend on it, and because formula_97. Note that the equations above are derived using the crack closure integral.
If the energy release rate exceed a critical value, the crack will grow. In this case, a new FEA simulation is performed (for the next time step) where the node at the crack tip is released. For a bounded substrate, we may simply stop enforcing fixed Dirichlet boundary conditions at the crack tip node of the previous time step (i.e. displacements are no longer restrained) . For a symmetric crack, we would need to update the geometry of the domain with a longer crack opening (and therefore generate a new mesh).
Modified crack closure integral.
Similar to the Nodal Release Method, the Modified Crack Closure Integral (MCCI) is a method for calculating the energy release rate utilizing FEA nodal displacements formula_98 and forces formula_99. Where formula_100 represents the direction corresponding to the Cartesian basis vectors with origin at the crack tip, and formula_101 represents the nodal index. MCCI is more computationally efficient than the nodal release method because it only requires one analysis for each increment of crack growth.
A necessary condition for the MCCI method is uniform element length formula_102 along the crack face in the formula_103direction. Additionally, this method requires sufficient discretization such that over the length of one element stress fields are self-similar. This implies that formula_104as the crack propagates. Below are examples of the MCCI method with two types of common finite elements.
4-node elements.
The 4-node square linear elements seen in Figure 2 have a distance between nodes formula_101 and formula_106 equal to formula_107 Consider a crack with its tip located at node formula_108 Similar to the nodal release method, if the crack were to propagate one element length along the line of symmetry (parallel to the formula_75-axis) the crack opening displacement would be the displacement at the previous crack tip, i.e. formula_110and the force at the new crack tip formula_111 would be formula_112 Since the crack growth is assumed to be self-similar the displacement at node formula_101 after the crack propagates is equal to the displacement at node formula_114 before the crack propagates. This same concept can be applied to the forces at node formula_106 and formula_108 Utilizing the same method shown in the nodal release section we recover the following equations for energy release rate:formula_117
formula_118
formula_119Where formula_120(displacement above and below the crack face respectively). Because we have a line of symmetry parallel to the crack, we can assume formula_121
Thus, formula_122
8-node elements.
The 8-node rectangular elements seen in Figure 3 have quadratic basis functions. The process for calculating G is the same as the 4-node elements with the exception that formula_67 (the crack growth over one element) is now the distance from node formula_101 to formula_125 Once again, making the assumption of self-similar straight crack growth the energy release rate can be calculated utilizing with the following equations:formula_126
formula_127
formula_128Like with the nodal release method the accuracy of MCCI is highly dependent on the level of discretization along the crack tip, i.e. formula_129 Accuracy also depends on element choice. A mesh of 8-node quadratic elements can produce more accurate results than a mesh of 4-node linear elements with the same number of degrees of freedom in the mesh.
Domain integral approach for J.
The J-integral may be calculated directly using the finite element mesh and shape functions. We consider a domain countour as shown in figure 4 and choose an arbitrary smooth function formula_130 such that formula_131 on formula_76 and formula_133 on formula_134.
For linear elastic cracks growing straight ahead, formula_71. The energy release rate can then be calculated over the area bounded by the contour using an updated formulation:
formula_136
The formula above may be applied to any annular area surrounding the crack tip (in particular, a set of neighboring elements can be used). This method is very accurate, even with a coarse mesh around the crack tip (one may choose an integration domain located far away, with stresses and displacement less sensitive to mesh refinment)
2-D crack tip singular elements.
The above-mentioned methods for calculating energy release rate asymptotically approach the actual solution with increased discretization, but fail to fully capture the crack tip singularity. More accurate simulations can be preformed by utilizing quarter-point elements around the crack tip. These elements have a built in singularity which more accurately produces stress fields around the crack tip. The advantage of the quarter-point method is that it allows for coarser finite element meshes and greatly reduces computational cost. Furthermore, these elements are derived from small modifications to common finite elements without requiring special computational programs for analysis. For the purposes of this section elastic materials will be examined, although this method can be extended to elastic-plastic fracture mechanics. Assuming perfect elasticity the stress fields will experience a formula_137 crack tip singularity.
8-node isoparametric element.
The 8-node quadratic element are described by Figure 5 in both parent space with local coordinates formula_138 and formula_139 and by the mapped element in physical/global space by formula_140 and formula_141 The parent element is mapped from the local space to the physical space by the shape functions formula_142 and the degree of freedom coordinates formula_143 The crack tip is located at formula_144 or formula_145formula_146formula_147In a similar way, displacements (defined as formula_148) can also be mapped.formula_149formula_150A property of shape functions in the finite element method is compact support, specifically the Kronecker delta property (i.e. formula_151 at node formula_152 and zero at all other nodes). This results in the following shape functions for the 8-node quadratic elements:formula_153
formula_154
formula_155
formula_156
formula_157
formula_158
formula_159
formula_160When considering a line in front of the crack that is co-linear with the formula_140- axis (i.e. formula_162) all basis functions are zero except for formula_163formula_164
formula_165
formula_166Calculating the normal strain involves using the chain rule to take the derivative of displacement with respect to formula_167formula_168
If the nodes are spaced evenly on the rectangular element then the strain will not contain the singularity. By moving nodes 5 and 8 position to a quarter of the length formula_169 of the element closer to the crack tip as seen in figure 5, the mapping from formula_170 becomes:formula_171Solving for formula_172 and taking the derivative results in:formula_173formula_174Plugging this result into the equation for strain the final result is obtained:formula_175By moving the mid-nodes to a quarter position results in the correct formula_137 crack tip singularity.
Other element types.
The rectangular element method does not allow for singular elements to be easily meshed around the crack tip. This impedes the ability to capture the angular dependence of the stress fields which is critical in determining the crack path. Also, except along the element edges the formula_177 singularity exists in a very small region near the crack tip. Figure 6 shows another quarter-point method for modeling this singularity. The 8-node rectangular element can be mapped into a triangle. This is done by collapsing the nodes on the line formula_178 to the mid-node location and shifting the mid-nodes on formula_179 to the quarter-point location. The collapsed rectangle can more easily surround the crack tip but requires that the element edges be straight or the accuracy of calculating the stress intensity factor will be reduced.
A better candidate for the quarter-point method is the natural triangle as seen in Figure 7. The element's geometry allows for the crack tip to be easily surrounded and meshing is simplified. Following the same procedure described above, the displacement and strain field for the triangular elements are:formula_180formula_181This method reproduces the first two terms of the Williams solutions with a constant and singular term.
An advantage of quarter-point method is that it can be easily generalized to 3-dimensional models. This can greatly reduce computation when compared to other 3-dimensional methods but can lead to error if that crack tip propagates with a large degree of curvature.