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.

import math

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

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:

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:

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

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

2a/W 0.1 | K/σ√(πa) 1.0053 | Feddersen 1.0062
2a/W 0.2 | K/σ√(πa) 1.0237 | Feddersen 1.0254
2a/W 0.3 | K/σ√(πa) 1.0565 | Feddersen 1.0594
2a/W 0.4 | K/σ√(πa) 1.1097 | Feddersen 1.1118
2a/W 0.5 | K/σ√(πa) 1.1868 | Feddersen 1.1892
2a/W 0.6 | K/σ√(πa) 1.3033 | Feddersen 1.3043
2a/W 0.7 | K/σ√(πa) 1.4891 | Feddersen 1.4841

Figure 1: The stress intensity factor of a central crack