Infracture mechanics, theenergy release rate,${\displaystyle G}$, is the rate at whichenergy istransformed as amaterial undergoesfracture. Mathematically, the energy release rate is expressed as the decrease in totalpotential energy per increase in fracture surface area,^[1][2] and is thus expressed in terms of energy per unit area. Variousenergy balances can be constructed relating the energy released during fracture to the energy of the resulting new surface, as well as otherdissipative processes such asplasticity 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 andfatigue.

The energy release rate${\displaystyle G}$ is defined^[3] as the instantaneous loss of total potential energy${\displaystyle \mathrm{\Pi }}$ per unit crack growth area${\displaystyle s}$,

${\displaystyle G\equiv -\frac{\mathrm{\partial }\mathrm{\Pi }}{\mathrm{\partial }s},}$

where the total potential energy is written in terms of the total strain energy${\displaystyle \mathrm{\Omega }}$, surface traction${\displaystyle \mathbf{t}}$, displacement${\displaystyle \mathbf{u}}$, and body force${\displaystyle \mathbf{b}}$ by

${\displaystyle \mathrm{\Pi }=\mathrm{\Omega }-\left\{{\int }_{{\mathcal{S}}_t}\mathbf{t}\cdot \mathbf{u}dS+{\int }_{\mathcal{V}}\mathbf{b}\cdot \mathbf{u}dV\right\}.}$

The first integral is over the surface${\displaystyle S_t}$ of the material, and the second is over its volume${\displaystyle V}$.

The figure on the right shows the plot of an external force${\displaystyle P}$ vs. the load-point displacement${\displaystyle q}$, in which the area under the curve is the strain energy. The white area between the curve and the${\displaystyle P}$-axis is referred to as the complementary energy. In the case of alinearly-elastic material,${\displaystyle P(q)}$ is a straight line and the strain energy is equal to the complementary energy.

In the case of prescribed displacement, the strain energy can be expressed in terms of the specified displacement and the crack surface${\displaystyle \mathrm{\Omega }(q,s)}$, and the change in this strain energy is only affected by the change in fracture surface area:${\displaystyle \delta \mathrm{\Omega }=(\mathrm{\partial }\mathrm{\Omega }/\mathrm{\partial }s)\delta s}$. Correspondingly, the energy release rate in this case is expressed as^[3]

${\displaystyle G=-{\frac{\mathrm{\partial }\mathrm{\Omega }}{\mathrm{\partial }s}|}_q.}$

Here is where one can accurately refer to${\displaystyle G}$ as the strain energy release rate.

When the load is prescribed instead of the displacement, the strain energy needs to be modified as${\displaystyle \mathrm{\Omega }(q(P,s),s)}$. The energy release rate is then computed as^[3]

${\displaystyle G=-{\frac{\mathrm{\partial }}{\mathrm{\partial }s}|}_P\left(\mathrm{\Omega }-Pq\right).}$

If the material is linearly-elastic, then${\displaystyle \mathrm{\Omega }=Pq/2}$ and one may instead write

${\displaystyle G={\frac{\mathrm{\partial }\mathrm{\Omega }}{\mathrm{\partial }s}|}_P.}$

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,${\displaystyle \mathrm{\partial }s=B\mathrm{\partial }a}$. Therefore, the equation for computing${\displaystyle G}$ can be modified for the 2D case:

• Prescribed Displacement:${\displaystyle G=-{\frac{1}{B}\frac{\mathrm{\partial }\mathrm{\Omega }}{\mathrm{\partial }a}|}_q.}$
• Prescribed Load:${\displaystyle G=-{\frac{1}{B}\frac{\mathrm{\partial }}{\mathrm{\partial }a}|}_P\left(\mathrm{\Omega }-Pq\right).}$
• Prescribed Load, Linear Elastic:${\displaystyle G={\frac{1}{B}\frac{\mathrm{\partial }\mathrm{\Omega }}{\mathrm{\partial }a}|}_P.}$

One can refer to the example calculations embedded in the next section for further information. Sometimes, the strain energy is written using${\displaystyle U=\mathrm{\Omega }/B}$, an energy-per-unit thickness. This gives

• Prescribed Displacement:${\displaystyle G=-{\frac{\mathrm{\partial }U}{\mathrm{\partial }a}|}_q.}$
• Prescribed Load:${\displaystyle G=-{\frac{\mathrm{\partial }}{\mathrm{\partial }a}|}_P\left(U-\frac{Pq}{B}\right).}$
• Prescribed Load, Linear Elastic:${\displaystyle G={\frac{\mathrm{\partial }U}{\mathrm{\partial }a}|}_P.}$

The energy release rate is directly related to thestress intensity factor associated with a given two-dimensional loading mode (Mode-I, Mode-II, or Mode-III) when the crack grows straight ahead.^[3] This is applicable to cracks underplane stress,plane strain, andantiplane shear.

For Mode-I, the energy release rate${\displaystyle G}$ rate isrelated to the Mode-I stress intensity factor${\displaystyle K_I}$ for a linearly-elastic material by

${\displaystyle G=\frac{K_I^2}{E^{\prime }},}$

where${\displaystyle E^{\prime }}$ is related toYoung's modulus${\displaystyle E}$ andPoisson's ratio${\displaystyle \nu }$ depending on whether the material is under plane stress or plane strain:

For Mode-II, the energy release rate is similarly written as

${\displaystyle G=\frac{K_{II}^2}{E^{\prime }}.}$

For Mode-III (antiplane shear), the energy release rate now is a function of theshear modulus${\displaystyle \mu }$,

${\displaystyle G=\frac{K_{III}^2}{2\mu }.}$

For an arbitrary combination of all loading modes, these linear elastic solutions may be superposed as

${\displaystyle G=\frac{K_I^2}{E^{\prime }}+\frac{K_{II}^2}{E^{\prime }}+\frac{K_{III}^2}{2\mu }.}$

Alternative representation of${\displaystyle G}$ in terms of${\displaystyle K}$[4]
Which can be seen to be equivalent to the previous representation through the relationship between Young's modulus and the shear modulus: ${\displaystyle E=2G(1+\nu )}$

Crack growth is initiated when the energy release rate overcomes a critical value${\displaystyle G_c}$, which is a material property,

${\displaystyle G\ge G_c,}$

Under Mode-I loading, the critical energy release rate${\displaystyle G_c}$ is then related to the Mode-Ifracture toughness${\displaystyle K_{IC}}$, another material property, by

${\displaystyle G_c=\frac{K_{IC}^2}{E^{\prime }}.}$

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.

The only method to calculate${\displaystyle G}$ 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:

• calculating the stress field resulting from the loading,
• calculating the strain energy in the material resulting from the stress field,
• calculating the work done by the external loads,

all in terms of the crack surface area.

Example calculation of${\displaystyle G}$ using the total potential energy

This problem is two-dimensional and has a fixed load, so with${\displaystyle s=aB}$, ${\displaystyle G=-\frac{1}{B}{\frac{\mathrm{\partial }}{\mathrm{\partial }a}|}_P(\mathrm{\Omega }-Pq).}$ Since the material is linearly-elastic,${\displaystyle \mathrm{\Omega }=Pq/2}$ and thus ${\displaystyle G=\frac{1}{B}\frac{\mathrm{\partial }\mathrm{\Omega }}{\mathrm{\partial }a}.}$ The stresses in the DCB are due to the bending stresses in each cantilever beam, ${\displaystyle {\sigma }_{11}(x,y)=-\frac{M(x)y}{I}=-\frac{12Pxy}{B{h}^3},}$ where${\displaystyle B}$ is the length into the page. Using${\displaystyle W={\sigma }_{11}^2/2E}$, one now has the strain energy ${\displaystyle \mathrm{\Omega }={\int }_{V}WdV=2B{\int }_{-h/2}^{h/2}{\int }_0^a\frac{1}{2E}{\left(\frac{12Pxy}{B{h}^3}\right)}^{2}dxdy,}$ where the factor of 2 out front is due to there being 2 cantilever beams. Solving, ${\displaystyle \mathrm{\Omega }=\frac{P^2}{B^3}\frac{4a^3}{E{h}^3},}$ then taking a derivative with respect to${\displaystyle a}$ and dividing by${\displaystyle B}$, ${\displaystyle G=\frac{P^2}{B^2}\frac{12a^2}{E{h}^3}.}$

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,${\displaystyle C}$^[3]

${\displaystyle C=\frac{q}{P}.}$

The corresponding strain energy${\displaystyle \mathrm{\Omega }}$ (area under the curve) is equal to^[3]

${\displaystyle \mathrm{\Omega }=\frac{1}{2}Pq=\frac{1}{2}\frac{q^2}{C}=\frac{1}{2}{P}^{2}C.}$

Using the compliance method, one can show that the energy release rate for both cases of prescribed load and displacement come out to be^[3]

${\displaystyle G=\frac{1}{2}{P}^2\frac{\mathrm{\partial }C}{\mathrm{\partial }s}.}$

Example calculation of${\displaystyle G}$ using the compliance method[3]

Consider a double cantilever beam (DCB) specimen as shown in the right figure. The displacement of a single beam is ${\displaystyle {\nu }_{max}=\frac{P{a}^3}{3EI},\text{with}I=\frac{B{h}^3}{12}.}$ The resulting load-point displacement${\displaystyle q}$ is therefore${\displaystyle 2{\nu }_{max}}$. Substitute into the equation for compliance and simplify: ${\displaystyle C=\frac{q}{P}=\frac{8a^3}{EB{h}^3}.}$ Now,${\displaystyle \mathrm{\partial }C/\mathrm{\partial }s}$ is computed as ${\displaystyle \frac{\mathrm{\partial }C}{\mathrm{\partial }s}=\frac{\mathrm{\partial }C}{\mathrm{\partial }a}\frac{\mathrm{\partial }a}{\mathrm{\partial }s}=\left(\frac{24a^2}{EB{h}^3}\right)\left(\frac{1}{B}\right).}$ Finally, the energy release rate of this DCB specimen can be expressed as ${\displaystyle G(P,a)=\frac{1}{2}{P}^2\frac{\mathrm{\partial }C}{\mathrm{\partial }s}=\frac{P^2}{B^2}\frac{12a^2}{E{h}^3}.}$ Note that alternatively, the energy release rate can be expressed in terms of${\displaystyle q}$ and${\displaystyle a}$: ${\displaystyle G(q,a)=\frac{3}{16}\frac{E{h}^3{q}^2}{a^4},}$ indicating that${\displaystyle G}$ decreases with the crack length${\displaystyle a}$ for the case of fixed displacement, and vice versa for the case of fixed load.

In the case of prescribed displacement, holding the crack length fixed, the energy release rate can be computed by^[3]

${\displaystyle G=-{\int }_0^q\frac{\mathrm{\partial }P}{\mathrm{\partial }s}dq,}$

while in the case of prescribed load,^[3]

${\displaystyle G={\int }_0^P\frac{\mathrm{\partial }q}{\mathrm{\partial }s}dP.}$

As one can see, in both cases, the energy release rate${\displaystyle G}$ times the change in surface${\displaystyle ds}$ returns the area between curves, which indicates the energy dissipated for the new surface area as illustrated in the right figure^[3]

${\displaystyle Gds=-ds{\int }_0^q\frac{\mathrm{\partial }P}{\mathrm{\partial }s}dq=ds{\int }_0^P\frac{\mathrm{\partial }q}{\mathrm{\partial }s}dP.}$

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^[3]

${\displaystyle G=\underset{\mathrm{\Delta }s\to 0}{lim}-\frac{1}{\mathrm{\Delta }s}{\int }_{\mathrm{\Delta }s}\frac{1}{2}{t}_i^0\left(\mathrm{\Delta }{u}_i^{+}-\mathrm{\Delta }{u}_i^{-}\right)dS,}$

where${\displaystyle \mathrm{\Delta }s}$ is the new fracture surface area,${\displaystyle t_i^0}$ are the components of the traction released on the top fracture surface as the crack grows,${\displaystyle \mathrm{\Delta }{u}_i^{+}-\mathrm{\Delta }{u}_i^{-}}$ 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${\displaystyle S}$.

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^[3]

${\displaystyle G=\underset{\mathrm{\Delta }a\to 0}{lim}\frac{1}{\mathrm{\Delta }a}{\int }_0^{\mathrm{\Delta }a}{\sigma }_{i2}(x_1,0)u_i(\mathrm{\Delta }a-x_1,\pi )d{x}_1,}$

where${\displaystyle \mathrm{\Delta }a}$ is the new crack length, and the displacement components are written as a function of the polar coordinates${\displaystyle r=\mathrm{\Delta }a-x_1}$ and${\displaystyle \theta =\pi }$.

Example calculation of${\displaystyle G}$ using the crack closure integral

Consider the crack in the DCB specimen shown in the figure. The nonzero stress and displacement components are given by[3] as ${\displaystyle {\sigma }_{22}(x_1,0)=\frac{K_I}{\sqrt{2\pi x_1}},u_{22}(\mathrm{\Delta }a-x_1,\pi )=\frac{8K_I}{E}\sqrt{\frac{\mathrm{\Delta }{a}_1-x_1}{2\pi }}.}$ The crack closure integral for this linearly elastic material, assuming the crack grows straight ahead, is ${\displaystyle G=\frac{2K_I^2}{\pi E}\underset{\mathrm{\Delta }a\to 0}{lim}\frac{1}{\mathrm{\Delta }a}{\int }_0^{\mathrm{\Delta }a}\sqrt{\frac{\mathrm{\Delta }{a}_1-x_1}{x_1}}d{x}_1.}$ Consider rescaling the integral using${\displaystyle \xi =x_1/\mathrm{\Delta }a}$ for ${\displaystyle G=\frac{2K_I^2}{\pi E}{\int }_0^1\sqrt{\frac{1-\xi }{\xi }}d\xi ,}$ where one computes the simpler integral to be ${\displaystyle {\int }_0^1\sqrt{\frac{1-\xi }{\xi }}d\xi =\frac{\pi }{2},}$ leaving the energy release rate as ${\displaystyle G=\frac{K_I^2}{E},}$ the expected relation. In this case, it is not straightforward to obtain${\displaystyle K_I}$ directly from the loading and geometry of the problem, but since the crack grows straight ahead and the material is linearly elastic, the energy release rate here should be the same as the energy release rate calculated using the other methods. This allows one to indirectly retrieve the stress intensity factor for this problem as ${\displaystyle K_I=\frac{P}{B}\frac{2a}{h}\sqrt{\frac{3}{h}}.}$

In certain situations, the energy release rate${\displaystyle G}$ can be calculated using theJ-integral, i.e.${\displaystyle G=J}$, using^[3]

${\displaystyle J={\int }_{\mathrm{\Gamma }}\left(W{n}_1-t_i\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_1}\right)d\mathrm{\Gamma },}$

where${\displaystyle W}$ is the elastic strain energy density,${\displaystyle n_1}$ is the${\displaystyle x_1}$ component of the unit vector normal to${\displaystyle \mathrm{\Gamma }}$, the curve used for the line integral,${\displaystyle t_i}$ are the components of the traction vector${\displaystyle \mathbf{t}=\mathit{\sigma }\cdot \mathbf{n}}$, where${\displaystyle \mathit{\sigma }}$ is the stress tensor, and${\displaystyle u_i}$ are the components of the displacement vector.

This integral is zero over a simple closed path and ispath independent, allowing any simple path starting and ending on the crack faces to be used to calculate${\displaystyle J}$.In order to equate the energy release rate to the J-integral,${\displaystyle G=J}$, the following conditions must be met:

• the crack must be growing straight ahead, and
• the deformation near the crack (enclosed by${\displaystyle \mathrm{\Gamma }}$) must be elastic (not plastic).

The J-integral may be calculated with these conditions violated, but then${\displaystyle G\ne J}$. 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^[3]

${\displaystyle G=J=\frac{K_I^2}{E^{\prime }}+\frac{K_{II}^2}{E^{\prime }}+\frac{K_{III}^2}{2\mu }.}$

Example calculation of${\displaystyle G}$ using the J-integral[3]

Consider the double cantilever beam specimen shown in the figure, where the crack centered in the beam of height${\displaystyle 2h}$ has a length of${\displaystyle a}$, and a load${\displaystyle P}$ is applied to open the crack. Assume that the material is linearly-elastic and that the crack grows straight forward. Consider a rectangular path shown in the second figure: start on the top crack face, (1) go up to the top at${\displaystyle h}$, (2) go to the right past the crack tip, (3) go down to the bottom at${\displaystyle -2h}$, (4) go along the bottom to the left, and (5) go back up to the bottom crack face. The J-integral is zero along many parts of this path. The material is effectively unloaded behind the crack, so both the strain energy density and traction are zero along (1) and (5), and hence the J-integral. Along (2) and (4) one has${\displaystyle n_1=0}$ as well as${\displaystyle \mathbf{t}=\mathbf{0}}$ (no traction on the free surface), so the J-integral is zero on (2) and (4) as well. This leaves only (3); assuming one is far enough from the crack on (3), the traction term is zero since${\displaystyle u_1=0}$ and${\displaystyle {\sigma }_{12}=0}$ far from the crack, leaving ${\displaystyle J={\int }_{{\mathrm{\Gamma }}_3}W{n}_{1}d{\mathrm{\Gamma }}_3.}$ ${\displaystyle n_1=1}$ along (3), and${\displaystyle W}$ is due to bending stress for a cantilever beam ${\displaystyle {\sigma }_{11}(y)=-\frac{My}{I}=-\frac{12Pay}{B{h}^3},}$ where${\displaystyle B}$ is the length into the page. Using${\displaystyle W={\sigma }_{11}^2/2E}$, one now has ${\displaystyle J=2{\int }_{-h/2}^{h/2}\frac{1}{2E}{\left(\frac{12Pay}{B{h}^3}\right)}^{2}dy,}$ where the factor of 2 out front is due to there being 2 cantilever beams. Solving, ${\displaystyle G=J=\frac{P^2}{B^2}\frac{12a^2}{E{h}^3}.}$

A handful of methods exist for calculating${\displaystyle G}$ with finite elements. Although a direct calculation of theJ-integral is possible (using the strains and stresses outputted byFEA), 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.

If the crack is growing straight, theenergy release rate can be decomposed as a sum of 3 terms${\displaystyle G_i}$ associated with the energy in each 3 modes. As a result, the Nodal Release method (NR) can be used to determine${\displaystyle G_i}$ 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${\displaystyle \mathrm{\Delta }a}$. First, we calculate the displacement variation at the node of interest${\displaystyle \mathrm{\Delta }\overrightarrow{u}={\overrightarrow{u}}^{(t+1)}-{\overrightarrow{u}}^{(t)}}$(before and after the crack tip node is released). Secondly, we keep track of the nodal force${\displaystyle \overrightarrow{F}}$ outputted by FEA. Finally, we can find each components of${\displaystyle G}$ using the following formulas:

$${\displaystyle G_1^{\text{NR}}=\frac{1}{\mathrm{\Delta }a}{F}_2\frac{\mathrm{\Delta }{u}_2}{2}}$$$${\displaystyle G_2^{\text{NR}}=\frac{1}{\mathrm{\Delta }a}{F}_1\frac{\mathrm{\Delta }{u}_1}{2}}$$$${\displaystyle G_3^{\text{NR}}=\frac{1}{\mathrm{\Delta }a}{F}_3\frac{\mathrm{\Delta }{u}_3}{2}}$$Where${\displaystyle \mathrm{\Delta }a}$ 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${\displaystyle G=\underset{\mathrm{\Delta }a\to 0}{lim}{G}^{\text{NR}}}$. Note that the equations above are derived using the crack closure integral.

If the energy release rate exceeds 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^[5]).

Similar to the Nodal Release Method, the Modified Crack Closure Integral (MCCI) is a method for calculating theenergy release rate utilizingFEA nodal displacements${\displaystyle (u_i^j)}$ and forces${\displaystyle (F_i^j)}$.^[6][7] Where${\displaystyle i}$ represents the direction corresponding to theCartesian basis vectors with origin at the crack tip, and${\displaystyle j}$ 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${\displaystyle (\mathrm{\Delta }a)}$ along the crack face in the${\displaystyle x_1}$-direction. Additionally, this method requires sufficient discretization such that over the length of one element stress fields areself-similar. This implies that${\displaystyle K(a+\mathrm{\Delta }a)\approx K(a)}$ as the crack propagates. Below are examples of the MCCI method with two types of common finite elements.

The 4-node square linear elements seen in Figure 2 have a distance between nodes${\displaystyle j}$ and${\displaystyle j+1}$ equal to${\displaystyle \mathrm{\Delta }a.}$ Consider a crack with its tip located at node${\displaystyle j.}$ Similar to the nodal release method, if the crack were to propagate one element length along the line of symmetry (parallel to the${\displaystyle x_1}$-axis) the crack opening displacement would be the displacement at the previous crack tip, i.e.${\displaystyle {\mathit{u}}^{\mathit{j}}}$ and the force at the new crack tip${\displaystyle (j+1)}$ would be${\displaystyle {\mathit{F}}^{j+1}.}$ Since the crack growth is assumed to be self-similar the displacement at node${\displaystyle j}$ after the crack propagates is equal to the displacement at node${\displaystyle j-1}$ before the crack propagates. This same concept can be applied to the forces at node${\displaystyle j+1}$ and${\displaystyle j.}$ Utilizing the same method shown in the nodal release section we recover the following equations for energy release rate:

${\displaystyle G_1^{\text{MCCI}}=\frac{1}{2\mathrm{\Delta }a}{F}_2^j\mathrm{\Delta }{u}_2^{j-1}}$

${\displaystyle G_2^{\text{MCCI}}=\frac{1}{2\mathrm{\Delta }a}{F}_1^j\mathrm{\Delta }{u}_1^{j-1}}$

${\displaystyle G_3^{\text{MCCI}}=\frac{1}{2\mathrm{\Delta }a}{F}_3^j\mathrm{\Delta }{u}_3^{j-1}}$

Where${\displaystyle \mathrm{\Delta }{u}_i^{j-1}=u_i^{(+)j-1}-u_i^{(-)j-1}}$(displacement above and below the crack face respectively). Because we have a line of symmetry parallel to the crack, we can assume${\displaystyle u_i^{(+)j-1}=-u_i^{(-)j-1}.}$

Thus,${\displaystyle \mathrm{\Delta }{u}_i^{j-1}=2u_i^{(+)j-1}.}$

The 8-node rectangular elements seen in Figure 3 have quadraticbasis functions. The process for calculating G is the same as the 4-node elements with the exception that${\displaystyle \mathrm{\Delta }a}$ (the crack growth over one element) is now the distance from node${\displaystyle j}$ to${\displaystyle j+2.}$ Once again, making the assumption of self-similar straight crack growth the energy release rate can be calculated utilizing the following equations:

${\displaystyle G_1^{\text{MCCI}}=\frac{1}{2\mathrm{\Delta }a}\left(F_2^j\mathrm{\Delta }{u}_2^{j-2}+F_2^{j+1}\mathrm{\Delta }{u}_2^{j-1}\right)}$

${\displaystyle G_2^{\text{MCCI}}=\frac{1}{2\mathrm{\Delta }a}\left(F_1^j\mathrm{\Delta }{u}_1^{j-2}+F_1^{j+1}\mathrm{\Delta }{u}_1^{j-1}\right)}$

${\displaystyle G_3^{\text{MCCI}}=\frac{1}{2\mathrm{\Delta }a}\left(F_3^j\mathrm{\Delta }{u}_3^{j-2}+F_3^{j+1}\mathrm{\Delta }{u}_3^{j-1}\right)}$

Like with the nodal release method the accuracy of MCCI is highly dependent on the level of discretization along the crack tip, i.e.${\displaystyle G=\underset{\mathrm{\Delta }a\to 0}{lim}{G}^{\text{MCCI}}.}$ 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^[8] in the mesh.

The J-integral may be calculated directly using the finite element mesh and shape functions.^[9] We consider a domain contour as shown in figure 4 and choose an arbitrary smooth function${\displaystyle \tilde{q}(x_1,x_2)=\sum _i{N}_i(x_1,x_2){\tilde{q}}_i}$ such that${\displaystyle \tilde{q}=1}$ on${\displaystyle \mathrm{\Gamma }}$ and${\displaystyle \tilde{q}=0}$ on${\displaystyle {\mathcal{C}}_1}$.

For linear elastic cracks growing straight ahead,${\displaystyle G=J}$. The energy release rate can then be calculated over the area bounded by the contour using an updated formulation:${\displaystyle J={\int }_{\mathcal{A}}({\sigma }_{ij}{u}_{i,1}{\tilde{q}}_{,j}-W{\tilde{q}}_{,1})d\mathcal{A}}$

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 refinement)

Derivation of the J-integral for domain integral approach

The J-intregral may be expressed over the full contour as follows: ${\displaystyle G=J={\int }_{\mathrm{\Gamma }}(W{n}_1-t_i\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_1})d\mathrm{\Gamma }={\int }_{{\mathcal{C}}_1+{\mathcal{C}}_{+}+{\mathcal{C}}_{-}-\mathcal{C}}(W{n}_1-{\sigma }_{ij}{n}_j\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_1})\tilde{q}d\mathrm{\Gamma }}$ With${\displaystyle \mathcal{C}={\mathcal{C}}_1+{\mathcal{C}}_{+}+{\mathcal{C}}_{-}-\mathrm{\Gamma }}$.${\displaystyle \overrightarrow{n}=-\overrightarrow{m}}$,${\displaystyle \tilde{q}=0}$ on${\displaystyle {\mathcal{C}}_1}$ and the work and stresses cancel out on${\displaystyle {\mathcal{C}}_{+}}$ and${\displaystyle {\mathcal{C}}_{-}}$, hence by application of thedivergence theorem this leads to: ${\displaystyle J={\int }_{\mathcal{C}}(-W{m}_1+{\sigma }_{ij}{m}_i\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_1})d\mathcal{C}={\int }_{\mathcal{A}}(\frac{\mathrm{\partial }{\sigma }_{ij}{u}_{i,1}\tilde{q}}{\mathrm{\partial }{x}_j}-\frac{\mathrm{\partial }W\tilde{q}}{\mathrm{\partial }{x}_1})d\mathcal{A}}$ Finally, by noting that${\displaystyle \frac{\mathrm{\partial }W}{\mathrm{\partial }{x}_1}={\sigma }_{ij}\frac{{\mathrm{\partial }}^2{u}_i}{\mathrm{\partial }{x}_j\mathrm{\partial }{x}_1}}$ and using the equilibrium equation: ${\displaystyle J={\int }_{\mathcal{A}}({\sigma }_{ij}{u}_{i,1}{\tilde{q}}_{,j}-W{\tilde{q}}_{,1})d\mathcal{A}}$

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 performed by utilizing quarter-point elements around the crack tip.^[10] 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.^{[11][12][13][14]} Assuming perfect elasticity the stress fields will experience a${\displaystyle \frac{1}{\sqrt{r}}}$ crack tip singularity.

The 8-node quadratic element is described by Figure 5 in both parent space with local coordinates${\displaystyle \xi }$ and${\displaystyle \eta ,}$ and by the mapped element in physical/global space by${\displaystyle x}$ and${\displaystyle y.}$ The parent element is mapped from the local space to the physical space by the shape functions${\displaystyle N_i(\xi ,\eta )}$ and the degree of freedom coordinates${\displaystyle (x_i,y_i).}$ The crack tip is located at${\displaystyle \xi =-1,\eta =-1}$ or${\displaystyle x=0,y=0.}$

${\displaystyle x(\xi ,\eta )=\sum _{i=1}^8{N}_i(\xi ,\eta )x_i}$

${\displaystyle y(\xi ,\eta )=\sum _{i=1}^8{N}_i(\xi ,\eta )y_i}$

In a similar way, displacements (defined as${\displaystyle u\equiv u_1,v\equiv u_2}$) can also be mapped.

${\displaystyle u(\xi ,\eta )=\sum _{i=1}^8{N}_i(\xi ,\eta )u_i}$

${\displaystyle v(\xi ,\eta )=\sum _{i=1}^8{N}_i(\xi ,\eta )v_i}$

A property of shape functions in the finite element method iscompact support, specifically theKronecker delta property (i.e.${\displaystyle N_i=1}$ at node${\displaystyle i}$ and zero at all other nodes). This results in the following shape functions for the 8-node quadratic elements:^[8]

${\displaystyle N_1=\frac{-(\xi -1)(\eta -1)(1+\eta +\xi )}{4}}$

${\displaystyle N_2=\frac{(\xi +1)(\eta -1)(1+\eta -\xi )}{4}}$

${\displaystyle N_3=\frac{(\xi +1)(\eta +1)(-1+\eta +\xi )}{4}}$

${\displaystyle N_4=\frac{-(\xi -1)(\eta +1)(-1+\eta -\xi )}{4}}$

${\displaystyle N_5=\frac{(1-{\xi }^2)(1-\eta )}{2}}$

${\displaystyle N_6=\frac{(1+\xi )(1-{\eta }^2)}{2}}$

${\displaystyle N_7=\frac{(1-{\xi }^2)(1+\eta )}{2}}$

${\displaystyle N_8=\frac{(1-\xi )(1-{\eta }^2)}{2}}$

When considering a line in front of the crack that is co-linear with the${\displaystyle x}$- axis (i.e.${\displaystyle N_i(\xi ,\eta =-1)}$) all basis functions are zero except for${\displaystyle N_{1,2,5}.}$

${\displaystyle N_1(\xi ,-1)=-\frac{\xi (1-\xi )}{2}}$

${\displaystyle N_2(\xi ,-1)=\frac{\xi (1+\xi )}{2}}$

${\displaystyle N_5(\xi ,-1)=(1-{\xi }^2)}$

Calculating the normal strain involves using thechain rule to take the derivative of displacement with respect to${\displaystyle x.}$

${\displaystyle {\gamma }_{xx}=\frac{\mathrm{\partial }u}{\mathrm{\partial }x}=\sum _{i=1,2,5}\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }\xi }\frac{\mathrm{\partial }\xi }{\mathrm{\partial }x}{u}_i}$

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${\displaystyle (\frac{L}{4})}$ of the element closer to the crack tip as seen in figure 5, the mapping from${\displaystyle \xi \to x}$ becomes:

${\displaystyle x(\xi )=\frac{\xi (1+\xi )}{2}L+(1-{\xi }^2)\frac{L}{4}}$

Solving for${\displaystyle \xi }$ and taking the derivative results in:

${\displaystyle \xi (x)=-1+2\sqrt{\frac{x}{L}}}$

${\displaystyle \frac{\mathrm{\partial }\xi }{\mathrm{\partial }x}=\frac{1}{\sqrt{xL}}}$

Plugging this result into the equation for strain the final result is obtained:

${\displaystyle {\gamma }_{xx}=\frac{4}{L}\left(\frac{u_2}{2}-u_5\right)+\frac{1}{\sqrt{xL}}\left(2u_5-\frac{u_2}{5}\right)}$

By moving the mid-nodes to a quarter position results in the correct${\displaystyle \frac{1}{\sqrt{r}}}$ crack tip singularity.

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${\displaystyle \frac{1}{\sqrt{r}}}$ 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.^[15] This is done by collapsing the nodes on the line${\displaystyle \xi =-1}$ to the mid-node location and shifting the mid-nodes on${\displaystyle \eta =\pm 1}$ 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:

${\displaystyle u=u_3+\sqrt{\frac{x}{L}}\left[4u_6-3u_3-u_1\right]+\frac{x}{L}\left[2u_1+2u_3-4u_6\right]}$

${\displaystyle {\gamma }_{xx}=\frac{\mathrm{\partial }u}{\mathrm{\partial }x}=\frac{1}{\sqrt{xL}}\left[-\frac{u_1}{2}-\frac{3u_3}{2}+2u_6\right]+\frac{1}{L}\left[2u_1+2u_3-4u_6\right]}$

This method reproduces the first two terms of the Williams solutions^[16] with a constant and singular term.

An advantage of the 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 errors if that crack tip propagates with a large degree of curvature.^[17]

• Fracture mechanics
• Stress intensity factor
• Fracture toughness
• J-integral

• Nonlinear Fracture Mechanics Notes by Prof. John Hutchinson (from Harvard University)
• Griffith's Strain Energy Release Rate onwww.fracturemechanics.org

• Fracture mechanics
• Solid mechanics
• Mechanics

• CS1 maint: location missing publisher
• Articles with short description
• Short description is different from Wikidata