---
jupyter:
  description: 'The stress intensity factor of a crack with planex2d and NGSolve: a central crack in a plate under
    tension from the energy that a longer crack releases, against Griffith''s infinite plate and Feddersen''s finite
    width correction.'
---

# The stress intensity factor of a crack

How strongly does a crack concentrate the stress at its tips? Close to a crack tip, the stress grows as $1/\sqrt{r}$
with the distance $r$ from it, and the stress intensity factor $K$ is the strength of that growth: ahead of the tip,
$\sigma_{yy} \approx K / \sqrt{2\pi r}$. A crack grows once $K$ reaches the fracture toughness of the material. For a
crack of length $2a$ across a wide plate under the tension $\sigma$,

$$K = \sigma \sqrt{\pi a}.$$

In a plate of finite width $W$, the free edges raise it, and Feddersen's correction gives
$K = \sigma \sqrt{\pi a} \sqrt{\sec(\pi a / W)}$.

Here, NGSolve finds $K$ from the energy that the crack releases as it grows, for cracks across a tenth to seven tenths of a
plate's width, and checks it against these values.

![The stress around a crack in a plate pulled at its ends, the crack's opening drawn larger](../_static/gallery/crack-stress-intensity.png){.result}

```python
import math

import ngsolve as ngs
import numpy as np
import planex as px
```

```python nbsphinx="hidden"
from docs_figures import Chart
```

## A quarter of the plate

The plate is 100 wide and 300 long, pulled at its ends with $\sigma = 1$, and the crack runs across its middle. By
symmetry, a quarter of the plate is enough. Its bottom edge lies on the crack's line: the crack itself is the part
of it that is free, and the ligament beyond the tip is held on the line, so the crack has no width at all. The edges
are named while drawing, and the two along the crack's line are meshed finely:

```python
W, H, E, nu = 100.0, 300.0, 1.0, 0.3
mu, lam = E / (2 * (1 + nu)), E * nu / (1 - nu**2)  # plane stress


def quarter(a):
    sk = px.QuickSketch((0, 0)).Line((a, 0), name="crack").Line((W / 2, 0), name="ligament")
    sk.Line((W / 2, H / 2), name="side").Line((0, H / 2), name="top").Close(name="axis")
    face = sk.Face()
    for edge in face.edges:
        if edge.name in ("crack", "ligament"):
            edge.maxh = a / 40
    return face
```

## The energy release rate

Under a fixed load, a crack that grows by $\delta A$ raises the stored strain energy $U$ by $G\,\delta A$, and the
energy release rate $G$ gives the stress intensity factor in plane stress as $K = \sqrt{E G}$. The full plate stores
four times the quarter's energy, and its crack grows by $2\,\delta a$ at two tips, so $G = 2\, dU_q/da$. `strain_energy`
solves the quarter, held along its axis and on the ligament, and returns half the work of the load:

```python
def strain_energy(a):
    mesh = ngs.Mesh(px.OCCGeometry(quarter(a)).GenerateMesh(maxh=W / 10, grading=0.2))
    V = ngs.VectorH1(mesh, order=5, dirichletx="axis", dirichlety="ligament")
    u, v = V.TnT()

    def stress(w):
        strain = ngs.Sym(ngs.Grad(w))
        return 2 * mu * strain + lam * ngs.Trace(strain) * ngs.Id(2)

    stiffness = ngs.BilinearForm(ngs.InnerProduct(stress(u), ngs.Sym(ngs.Grad(v))) * ngs.dx).Assemble()
    load = ngs.LinearForm(v[1] * ngs.ds("top")).Assemble()
    displacement = ngs.GridFunction(V)
    displacement.vec.data = stiffness.mat.Inverse(V.FreeDofs(), inverse="umfpack") * load.vec
    return ngs.InnerProduct(load.vec, displacement.vec) / 2


def intensity(a, step=0.005):
    G = 2 * (strain_energy(a * (1 + step)) - strain_energy(a * (1 - step))) / (2 * a * step)
    return math.sqrt(E * G)
```

## Against the formulas

For cracks across a tenth to seven tenths of the plate's width, $K / (\sigma\sqrt{\pi a})$ follows Feddersen's
correction to within half a percent. The wider the crack, the further it rises above Griffith's 1 for the wide plate
(Figure 1):

```python
ratios = [0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35]
factors = {}
for ratio in ratios:
    a = ratio * W
    factors[ratio] = intensity(a) / math.sqrt(math.pi * a)
    print("2a/W", 2 * ratio, "| K/σ√(πa)", round(factors[ratio], 4), "| Feddersen",
          round(math.sqrt(1 / math.cos(math.pi * ratio)), 4))
```

```python nbsphinx="hide-input"
fine = np.linspace(0.02, 0.37, 50)
Chart("crack across the width, 2a/W", "K / σ√(πa)",
      dots={"finite elements": ([2 * r for r in ratios], list(factors.values()))},
      lines={"Feddersen, √sec(πa/W)": (2 * fine, np.sqrt(1 / np.cos(np.pi * fine))), "Griffith, the wide plate": ([0, 0.74], [1, 1])},
      title="The stress intensity factor of a central crack")
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch, Polygon
from docs_figures import INK, MAP
from gallery import save_figure
a = 0.25 * W
mesh = ngs.Mesh(px.OCCGeometry(quarter(a)).GenerateMesh(maxh=W / 10, grading=0.2))
V = ngs.VectorH1(mesh, order=5, dirichletx="axis", dirichlety="ligament")
u, v = V.TnT()
strain = lambda w: ngs.Sym(ngs.Grad(w))
sig = lambda w: 2 * mu * strain(w) + lam * ngs.Trace(strain(w)) * ngs.Id(2)
stiffness = ngs.BilinearForm(ngs.InnerProduct(sig(u), strain(v)) * ngs.dx).Assemble()
load = ngs.LinearForm(v[1] * ngs.ds("top")).Assemble()
disp = ngs.GridFunction(V)
disp.vec.data = stiffness.mat.Inverse(V.FreeDofs(), inverse="umfpack") * load.vec
syy = sig(disp)[1, 1]


def turned(face, plot=True, **_):
    """σyy for PlotField on the middle of the plate turned a quarter, the quarter's values mirrored to all four"""
    gx, gy = np.meshgrid(np.linspace(-75, 75, 451), np.linspace(-50, 50, 301))
    across = np.clip(np.abs(gy.ravel()), 1e-6, W / 2 - 1e-6)
    along = np.clip(np.abs(gx.ravel()), 1e-6, H / 2)
    return {"points": list(zip(gx.ravel(), gy.ravel())), "values": np.asarray(syy(mesh(across, along))).ravel(),
            "label": "σyy / σ"}


fig, ax = plt.subplots(figsize=(6.4, 4))
px.PlotField(turned, px.QuickSketch((-75, -50)).Rectangle(150, 100).Face(), ax=ax, colormap=MAP.reversed(), min=0,
             max=3)
s = np.linspace(-a, a, 121)
uy = np.array([disp(mesh(min(abs(q), a - 1e-6), 1e-9))[1] for q in s])
k = 3.0 / uy.max()  # the opening drawn 3 wide
lens = [(k * w, q) for q, w in zip(s, uy)] + [(-k * w, q) for q, w in zip(s[::-1], uy[::-1])]
ax.add_patch(Polygon(lens, facecolor="white", edgecolor=INK, lw=1.2, zorder=5))
for side in (-1, 1):
    for q in (-32, 0, 32):
        ax.add_patch(FancyArrowPatch((side * 78, q), (side * 92, q), arrowstyle="-|>", mutation_scale=14, color=INK,
                                     lw=1.6))
ax.set_aspect("equal", adjustable="box")
ax.set_xlim(-95, 95)
ax.set_ylim(-53, 53)
save_figure(fig, "crack-stress-intensity", bare=True, bars={1: "σyy / σ"})
```