Search Fracture Mechanics Website

Determining K from Displacements

 home > LEFM > k from displacements
Prev Up Next

Introduction

This page addresses the challenge of determining the stress intensity factor, \(K\), from crack opening displacement data. The data can come from experimental measurements or analytical predictions. The idea is simple: take the equation for crack opening displacement as a function of \(K\) and invert it to obtain \(K\) as a function of displacement, and then insert displacement data into it. We will begin with some introductory exercises, and then focus on determining \(K\) from finite element analysis (FEA) results.

But why bother trying to estimate \(K\) from displacement data in the first place? Why not simply use \(K = \sigma_\infty \sqrt{\pi a}\) to compute it? The answer is that this equation only applies to an infinite plate with a center crack, and real world cases are never this simple. They usually involve shear, bending, and complex geometries that complicate the issue and lead to K-values different from \(K = \sigma_\infty \sqrt{\pi a}\).


Review of Pertinent Equations

The web page, displacements.html, developed the complete set of displacement equations for numerous scenarios: full displacement field, crack tip approximation, crack opening displacement, plane stress conditions, plane strain conditions, etc. But of all of these, only two equations are relevant to the topic at hand. They are (i) the exact solution for crack opening displacements derived from Westergaard's exact solution for stress

\[ d = {\sigma_\infty \over 2G} ( \kappa + 1 ) \, a \; \sqrt{1 - \left( x \over a \right)^2} \]
and (ii) Irwin's crack tip approximation for crack opening displacements in terms of \(K\), the stress intensity factor.

\[ d = {K \over G} \sqrt{ r \over 2 \pi} \; ( \kappa + 1 ) \]
The variables in the two equations are: \(x\), \(r\), \(a\), and \(d\), defined in the plot below, \(G\) is the material's shear modulus, \(\sigma_\infty\) is the far field stress, \(K\) is the stress intensity factor defined as \(K = \sigma_\infty \sqrt{\pi a}\), and \(\kappa\) is the Kolosov constant defined as

\[ \kappa \; = \; \left\{ \begin{eqnarray} 3 - 4 \nu & \qquad \text{Plane Strain} \\ \\ {3 - \nu \over 1 + \nu} & \qquad \text{Plane Stress} \end{eqnarray} \right. \]
The exact equation already contains \(\sigma_\infty\). We could easily multiply it by \((\sqrt{\pi a} / \sqrt{\pi a})\) to obtain \(\sigma_\infty \sqrt{\pi a}\) and thereby relate \(K\) to the exact crack opening displacement, \(d\). That's fine enough. But it only works when we have an exact theoretical solution to a problem. And we only have that for the center crack in an infinite plate. So this approach will not work in real-world problems with complicated geometries and loading conditions.

The key issue and the principal motivation for this page is the fact that as one moves away from the crack tip, the approximate solution for crack opening displacement diverges from the true, exact value. This is reflected in the plot below, which compares the two for the case of 70 MPa tension in a plate. This difference greatly complicates the task of accurately determining \(K\) from displacement data.

Crack Opening Example



Determining K from Crack Opening Displacements

The method of determining \(K\) from crack displacements is simple in concept. It is simply a matter of inverting the above equation for displacements in terms of \(K\) and solving for \(K\) in terms of \(d\).

\[ K = \sqrt{2 \pi} \; {G \over \kappa + 1} \left( {d \over \sqrt{r}} \right) \]
A notable feature of this equation is that only \(d\) and \(\sqrt{r}\) vary with \(r\). All other terms are constant. Actually, near the crack tip, the ratio of \(d\) to \(\sqrt{r}\) is constant as well. This reflects the fact that \(K\) can only have one single value for a single loading condition.

The challenge is that the equation is only accurate at the crack tip because it is the result of applying Irwin's near crack tip approximation. As one moves away from the tip, the predicted displacement by the approximate equation diverges from the exact solution. This leads to inaccuracies in estimating \(K\) from displacements as the example below will demonstrate.

Determining \(K\) at Two Positions

Let's return to the loading case from displacements.html, which consisted of a 60mm center crack (\(a\) = 30mm) in an infinite aluminum plate subjected to \(\sigma_\infty = \text{70 MPa}\) tension. Recall the material properties for aluminum are \(G = \text{26,000 MPa}\) and \(\nu = 1/3\).

Let's start by first computing \(K\) directly - and exactly - as follows.

\[ \begin{eqnarray} K & = & \sigma_\infty \sqrt{\pi a} \\ \\ & = & (70\,\text{MPa}) \sqrt{\pi\,(0.03\,\text{m})} \\ \\ & = & 21.49\,\text{MPa}\sqrt{\text{m}} \end{eqnarray} \]
So now we know \(21.49\,\text{MPa}\sqrt{\text{m}}\) is the exact, correct value. Now let's estimate \(K\) using displacement data, and do so for plane stress conditions. Therefore \(\kappa = 2\) and \(\kappa + 1 = 3\). We will perform the estimates at two different positions. The first, 10mm from the tip, so \(x / a = 66.7\text{%}\). The second, 3mm from the tip, so \(x / a = 90\text{%}\).

Start by inserting \(x / a = 0.667\) into the exact solution for crack opening displacement.

\[ \begin{eqnarray} d & = & {\sigma_\infty \over 2G} ( \kappa + 1 ) \, a \; \sqrt{1 - \left( x \over a \right)^2} \\ \\ & = & {70 \, \text{MPa} \over 2 (26,000 \, \text{MPa} ) } ( 3 ) \, ( 30 \, \text{mm} ) \sqrt{1 - 0.667^2} \\ \\ & = & 0.0903 \, \text{mm} \end{eqnarray} \]
Now let's insert this displacement into the equation for \(K\) in terms of \(d\) to obtain an estimate for it.

\[ \begin{eqnarray} K & = & \sqrt{2 \pi} \; {G \over \kappa + 1} \left( {d \over \sqrt{r}} \right) \\ \\ & = & \sqrt{2 \pi} \; {26,000 \, \text{MPa} \over 3} \left( {0.0000903 \, \text{m} \over \sqrt{0.01 \, \text{m}}} \right) \\ \\ & = & 19.62 \, \text{MPa} \sqrt{\text{m}} \end{eqnarray} \]
We have arrived at an estimate that is 8.7% lower than the exact value of \(21.49\,\text{MPa}\sqrt{\text{m}}\). Again, the source of this error is the fact that the equation relating \(K\) and \(d\) loses accuracy away from the crack tip. So even though we used the exact displacement value, we only had an approximate equation to work with.

Let's move closer to the crack tip and try again. Choose \(r = 3 \, \text{mm}\), giving \(x / a = 90\text{%}\). The exact displacement here is

\[ \begin{eqnarray} d & = & {\sigma_\infty \over 2G} ( \kappa + 1 ) \, a \; \sqrt{1 - \left( x \over a \right)^2} \\ \\ & = & {70 \, \text{MPa} \over 2 (26,000 \, \text{MPa} ) } ( 3 ) \, ( 30 \, \text{mm} ) \sqrt{1 - 0.90^2} \\ \\ & = & 0.0528 \, \text{mm} \end{eqnarray} \]
Now plug this displacement value into the equation for \(K\).

\[ \begin{eqnarray} K & = & \sqrt{2 \pi} \; {G \over \kappa + 1} \left( {d \over \sqrt{r}} \right) \\ \\ & = & \sqrt{2 \pi} \; {26,000 \, \text{MPa} \over 3} \left( {0.0000528 \, \text{m} \over \sqrt{0.003 \, \text{m}}} \right) \\ \\ & = & 20.94 \, \text{MPa} \sqrt{\text{m}} \end{eqnarray} \]
We are now within 2.6% of the exact value. So moving closer to the tip has improved the estimate as expected.

Expression for \(K\) Estimate Error

The above example can be generalized to obtain an equation giving the error in \(K\) from displacement values as a function of distance from the crack tip.

Recall Westergaard's exact solution for the crack opening displacement was

\[ d = {\sigma_\infty \over 2G} (\kappa + 1) \, a\, \sqrt{1 - (x/a)^2} \]
Substitute this into the expression for \(K\).

\[ K = \sqrt{2 \pi} \; {G \over \kappa + 1} \left( { {\sigma_\infty \over 2G} (\kappa + 1) \, a\, \sqrt{1 - (x/a)^2} \over \sqrt{r}} \right) \]
We obviously need to simplify this. But we also need to relate \(r\) to \(x\). That relationship is \(r = a - x\) for \(0 \le x \le a\). Substituting and simplifying leads to

\[ K = \sigma_\infty \sqrt{\pi \over 2} \; \sqrt{a + x} \]
And this can be further manipulated to obtain

\[ K = \sigma_\infty \sqrt{\pi a} \; \sqrt{1 + (x/a) \over 2} \]
Remarkably, \(\sigma_\infty \sqrt{\pi a}\) appears in this equation. It is the exact value for \(K\). So the equation can be written as

\[ {K \over K_\text{exact}} = \sqrt{1 + (x/a) \over 2} \]
which is a direct measure of the error in estimating \(K\) from displacement data measured at position, \(x\). Note that \(K \rightarrow K_\text{exact}\) as \(x \rightarrow a\), as it should.

The plot below shows the computed ratio of \(K\) to \(K_\text{exact}\) as a function of dimensionless position, \(x / a\), where the crack opening displacement is recorded, either by measurement or prediction.

Crack Opening Example

The graph here inspires an idea. We could take our estimated \(K\) values in the above example and scale them up by the amount the chart tells us they are underpredicted. For example, look at the case in the above example where we evaluated \(K\) 10mm away from the crack tip, at \(x/a = 0.667\). The chart tells us that \(K / K_\text{exact} = 0.913\). And recall we had computed a value of \(K = 19.62 \, \text{MPa} \sqrt{\text{m}}\) here. Dividing 19.62 by 0.913 gives us \(K = 21.49 \, \text{MPa} \sqrt{\text{m}}\). Exactly the correct value!


Estimating \(K\) by Extrapolation

Crack Opening Example
But while the above correction by scaling process may look foolproof, it is not. Consider the following case of a crack growing along a crooked path. In such a case, there could easily be three different values of crack length, \(a\), as shown in the sketch. In fact, there could be even more because the three values could be averaged or processed in some other way to obtain a still different effective value. So it is not at all clear what crack length to use to determine the correct scaling factor.

The resolution to this dilemma is to abandon the scaling factor approach all together. It just cannot be generalized to handle complex situations. Instead, the recommended procedure to use when estimating \(K\) from displacement data is to do something that is normally frowned upon. It is to extrapolate! Basically, take two measurements at different distances from the crack tip and extrapolate the computed \(K\) values back to \(r = 0\). Mathematically, this looks like

\[ {0 - r_1 \over r_2 - r_1} = {K_{(r=0)} - K_{r_1} \over K_{r_2} - K_{r_1}} \]
where \(K_{(r=0)}\) is the desired/extrapolated value at \(r=0\). Let's take the values from our example above and give it a go. Recall we obtained \(K = 19.62 \, \text{MPa} \sqrt{\text{m}}\) at \(r = 10\,\text{mm}\) and \(K = 20.94 \, \text{MPa} \sqrt{\text{m}}\) at \(r = 3 \,\text{mm}\). Inserting these values into the equation

\[ {0 - 3\,\text{mm} \over 10\,\text{mm} - 3\,\text{mm}} = {K_{(r=0)} - 20.94 \over 19.62 - 20.94} \]
gives \(K = 21.51 \, \text{MPa} \sqrt{\text{m}}\). Plenty close enough. And if the two evaluation locations were chosen to be closer to the crack tip to begin with, then the extrapolation would be even more accurate.

Note that this method works regardless of the complexity of the problem at hand. It relies only on measuring distance along the crack from its tip. In contrast, the earlier scaling method relied on having an analytical solution for the problem, something that is seldom available.

Finally, note that either location can be called \(r_1\) or \(r_2\) in the above equation as long as the corresponding \(K\) values are assigned to \(K_1\) and \(K_2\) appropriately. The exact same result will be obtained regardless.


Estimating \(K\) Using Finite Element Analysis

Finite element analyses (FEA) provide an ideal framework for performing the extrapolation steps presented above because of the many nodes present in a mesh. They provide ample opportunities to compute \(K\) at multiple locations and then extrapolate those values to \(r = 0\).

But the mesh must satisfy certain geometric guidelines in order to provide accurate predictions of the displacement, stress, and strain fields near the crack. Otherwise, the predictions can be wildly incorrect.

Proper LEFM FE Mesh
An ideal mesh for Linear Elastic Fracture Mechanics (LEFM) analysis is shown in the figure here. Its key property is that the triangular elements are quadratic with midside nodes placed at 25% of the element's edge. These midside nodes are called quarter-point nodes to emphasize that they are at one-quarter of the element's edge. The justification for this 25% node location was explained on the feamesh.html page. (Note - Quadratic wedge elements with quarter-point nodes are used in 3-D analyses.)

The extrapolation equation for this case is

\[ {0 - L/4 \over L - L/4} = {K_{(r=0)} - K_{L/4} \over K_L - K_{L/4}} \]
where \(L\) is the element edge length, \(K_L\) is evaluated at element vertex node, and \(K_{L/4}\) is evaluated at the quarter-point node. Rearranging the equation gives

\[ K_{(r=0)} = {4 \over 3} K_{L/4} - { 1 \over 3} K_L \]
This is the key result of this section. It is part of the generalized interpolation equation in terms of \(f\) with \(f = 4/3\).

\[ K_{(r=0)} = f * K_{L/4} + (1 - f) * K_L \]
We have \(0 \le f \le 1\) for conventional interpolations and \(f \ge 1\) for extrapolations.

Determining \(K\) from FEA

Let's return one last time to the loading case from displacements.html, which consisted of a 60mm crack (\(a\) = 30mm) in an aluminum plate subjected to \(\sigma_\infty = \text{70 MPa}\) tension. Recall the exact value of \(K\) for this case is \(21.49\,\text{MPa}\sqrt{\text{m}}\).

Assume the element size at the crack tip is 12mm. So \(L = 12 \, \text{mm}\) and \(L/4 = 3 \, \text{mm}\).

Let's assume for the sake of argument that the finite element analysis gives the exact displacement predictions given by Westergaard's solution. We will need \(x/a\) values for this. They are \(x/a = 0.6\) when \(L = 12\,\text{mm}\) and \(x/a = 0.9\) when \(L/4 = 3\,\text{mm}\).

Inserting these values into

\[ d = {\sigma_\infty \over 2G} ( \kappa + 1 ) \, a \; \sqrt{1 - \left( x \over a \right)^2} \]
gives \(d = 0.0969\,\text{mm}\) for \(L = 12\,\text{mm}\) and \(d = 0.0528\,\text{mm}\) for \(L/4 = 3\,\text{mm}\). These displacements are then inserted into the equation for \(K\) in terms of \(d\).

\[ K = \sqrt{2 \pi} \; {G \over \kappa + 1} \left( {d \over \sqrt{r}} \right) \]
This gives \(K_L = 19.22 \, \text{MPa} \sqrt{\text{m}}\) for \(L = 12\,\text{mm}\) and \(K_{L/4} = 20.94 \, \text{MPa} \sqrt{\text{m}}\) for \(L/4 = 3\,\text{mm}\).

Finally, insert these two values into the extrapolation equation to obtain

\[ \begin{eqnarray} K_{(r=0)} & = & {4 \over 3} K_{L/4} - {1 \over 3} K_L \\ \\ & = & {4 \over 3} (20.94) - {1 \over 3} (19.22) \\ \\ & = & 21.51 \, \text{MPa} \sqrt{\text{m}} \end{eqnarray} \]
Unsurprisingly, this is the same value that was obtained earlier by extrapolating from 10mm and 3mm away from the crack tip.

Summary

We have seen that the best method of determining \(K\) from crack opening displacement data is to extrapolate Irwin's approximate solution in terms of the square-root function to the crack tip at \(r = 0\). This requires that \(K\) first be computed at two different \(r\) values using

\[ K = \sqrt{2 \pi} \; {G \over \kappa + 1} \left( {d \over \sqrt{r}} \right) \]
Then the two \(K\) values are used to extrapolate back to \(r=0\).

A special case of this occurs in finite element analysis when node displacements are used in the calculations. In such cases, quadratic elements should be used and the midside node should be at 25% of the distance from the crack tip to the vertex node. In such cases, the proper extrapolation equation is

\[ K_{(r=0)} = {4 \over 3} K_{L/4} - {1 \over 3} K_L \]