Introduction
This page addresses 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.
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.
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
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.
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
\]
