The J-Integral and Fracture Mechanics

Overview

A part fails by fracture when a crack grows through it — and a crack can grow at a load far below the load that would yield the uncracked material. That is why fracture mechanics matters to engineers: the yield stress σy is not the whole story. Every real material contains flaws (from machining, welding, fatigue, or corrosion), and the question is whether the load you apply will drive one of them to grow.

Two numbers answer that question:

• The driving force — how hard the load pushes the crack. It is measured by the stress-intensity factor K (Pa·√m) or, for general situations, by the J-integral (energy released per unit crack advance, J/m²).
• The resistance — how much the material can take. It is measured by the fracture toughness K_Ic, a material property tabulated in handbooks alongside E and σy (Rasmah’s material library carries it for every metal).

For a linear-elastic solid the J-integral (Rice 1968) is a contour integral around the crack tip,

\[J = \oint_\Gamma \Bigl( W\, n_1 - \sigma_{ij} n_j \frac{\partial u_i}{\partial x_1} \Bigr) \, ds, \qquad W = \tfrac12 \sigma_{ij} \varepsilon_{ij},\]

whose defining property is path independence: any contour that encircles the tip gives the same value. For a mode-I (opening) field J equals the energy release rate J = K_I²/E′, with E′ = E in plane stress and E′ = E/(1−ν²) in plane strain.

In this tutorial we (1) verify those two claims numerically, (2) run the assessment an engineer actually makes — “given this working stress, how big a crack is fatal?” — for two real alloys, and (3) compute J from a full finite-element solve of a cracked plate and check it against the handbook stress-intensity factor.

using Rasmah

The mode-I K-field

Near a crack tip the displacement and stress have a universal 1/√r shape, parametrised by the stress-intensity factor K. Close enough to the tip the material does not care about the rest of the part — only about K and the material constants. Evaluate the analytic field at the instant of fracture in 7075-T6 aluminium, where K equals the material’s fracture toughness:

mat = aluminum_7075
K = fracture_toughness(mat)              # Pa·√m — the driving force at fracture
E, ν = youngs_modulus(mat), poisson_ratio(mat)
μ = shear_modulus(mat)

ux, uy = mode1_displacement(K, μ, ν, 0.01, 0.02; plane = :strain)
σ = mode1_stress(K, 0.01, 0.02)
println("u at (0.01, 0.02) = (", round(ux; digits = 8), ", ", round(uy; digits = 8), ") m")
println("σxx, σyy, σxy     = ", round.(σ; digits = 0), " Pa")

u at (0.01, 0.02) = (3.365e-5, 2.08e-5) m
σxx, σyy, σxy     = [3.1353241e7, 1.00274333e8, -3.107306e6] Pa

At 1 cm from the tip the opening stress is already ≈100 MPa — comparable to a design stress, from a load-free crack that was “just there”. Move ten times closer (r = 1 mm) and the stress grows by √10 ≈ 3× (1/√r), which is why the region where the elastic prediction exceeds σy — the plastic zone — is confined to a small disc at the tip. We quantify that disc below.

Verification — path independence and J = K²/E′

Evaluate the contour J-integral on circles of two different radii and check they agree with each other and with K²/E′:

E′ = E / (1 - ν^2)          # plane strain
J1 = j_integral_contour(K, E, ν; plane = :strain, R = 0.5)
J2 = j_integral_contour(K, E, ν; plane = :strain, R = 1.5)
println("J (R=0.5) = ", round(J1; digits = 1), "  J (R=1.5) = ", round(J2; digits = 1))
println("K²/E′     = ", round(K^2 / E′; digits = 1))

J (R=0.5) = 10452.1  J (R=1.5) = 10452.1
K²/E′     = 10452.1

The three numbers agree: the integral is path-independent and equals the energy release rate K²/E′. For aluminium that release rate at fracture is a deceptively small ≈10.5 kJ/m² — one Joule per square metre of crack advance is all it takes.

Fracture assessment of a cracked plate

Here is the calculation an engineer makes. For an edge crack of depth a in a wide plate under a remote tensile stress σ perpendicular to the crack, the stress-intensity factor is the classic handbook expression

\[K_I = \sigma \sqrt{\pi a}\, F(a/W),\]

with the geometry factor F tabulated for many configurations. For a shallow edge crack F ≈ 1.12; for deeper cracks F grows (the finite-width polynomial below is the one we verify against in the FEM section). Setting K_I = K_Ic and solving for a gives the critical crack size — the crack depth at which the part fails at that stress. Compute it for two alloys at a working stress of 150 MPa (a quarter to a third of their yield strengths):

σwork = 150.0e6                          # Pa — working stress
for m in (aluminum_7075, steel_4340)
    KIc = fracture_toughness(m)
    a_cr = (KIc / (1.12 * σwork))^2 / π  # K_Ic = 1.12·σ·√(π·a_cr) inverted
    println(m.name, ": KIc = ", round(KIc / 1e6; digits = 0),
            " MPa√m, σy = ", round(yield_strength(m) / 1e6; digits = 0),
            " MPa → critical crack a_cr = ", round(a_cr * 1e3; digits = 1), " mm")
end

Aluminum (7075-T6): KIc = 29.0 MPa√m, σy = 503.0 MPa → critical crack a_cr = 9.5 mm
Steel (4340): KIc = 120.0 MPa√m, σy = 860.0 MPa → critical crack a_cr = 162.4 mm

The contrast is the whole point of fracture mechanics: at the same working stress, 7075-T6 aluminium fails at a 9.5 mm edge crack, while 4340 steel tolerates a 162 mm one. Neither part would have yielded — 150 MPa is well below both yield strengths — yet the aluminium part fails first. This is why airframes are inspected for millimetre-scale cracks, and why the inspection interval (how often you must look) follows from a_cr and the crack growth rate.

rp = (K / yield_strength(aluminum_7075))^2 / (2π)
println("plastic zone rp ≈ ", round(rp * 1e3; digits = 2), " mm   vs   a_cr = 9.5 mm")

plastic zone rp ≈ 0.53 mm   vs   a_cr = 9.5 mm

a = 5.0e-3
F = 1.12 - 0.231 * a + 10.55 * a^2 - 21.72 * a^3 + 30.39 * a^4
K5 = σwork * sqrt(π * a) * F
println("K(5 mm crack) = ", round(K5 / 1e6; digits = 1), " MPa√m   safety factor KIc/K = ",
        round(K / K5; digits = 2))

K(5 mm crack) = 21.0 MPa√m   safety factor KIc/K = 1.38

The domain J-integral of a FEM solution

The contour form needs a smooth field to integrate along a path. A finite-element solution is piecewise linear, so instead of contours, fracture codes use the equivalent domain integral (Shih–Moran–Nakamura 1986): the same value, rewritten with a smooth weight q that is 1 at the tip and 0 on a surrounding ring, integrated over the elements between the two. Rasmah provides a cracked-plate mesh with duplicated nodes along the crack faces (so the faces can open) and the domain integral as a one-shot solve-and-evaluate:

m, bottom = cracked_grid_mesh(120, 60; crack_end = 0.1)
println("cracked mesh: ", size(m.vertices, 2), " nodes, ", size(m.faces, 2), " triangles")
vtk(m)

cracked mesh: 7393 nodes, 14400 triangles

The crack runs along y = ½ from the left edge to x = 0.1; the duplicated nodes let the two faces slide past each other, so the crack can actually open. Now load the plate: fix the bottom edge and pull the top edge in +y — tension perpendicular to the crack, which is what opens it in mode I:

n = size(m.vertices, 2)
σfar = 40.0e6                             # 40 MPa — keeps K below KIc of 7075
loads = zeros(2n)
fixed = Int[]
for i in 1:n
    y = m.vertices[2, i]
    y ≈ 0.0 && append!(fixed, 2i - 1, 2i)  # fix the bottom edge
    y ≈ 1.0 && (loads[2i] = σfar / 120)    # distribute σfar on the top edge
end

With the crack tip at (0.1, 0.5), define the q plateau (constant inside a 0.03 radius, smoothly falling to 0 at 0.2), solve the plane-stress elasticity problem, and evaluate the domain integral:

tip = (0.1, 0.5)
rin, rout = 0.03, 0.2
qgrad(x, y) = begin
    r = sqrt((x - tip[1])^2 + (y - tip[2])^2)
    dq = 0.0
    if rin < r < rout
        t = (r - rin) / (rout - rin)
        dq = -6t * (1 - t) / (rout - rin)
    end
    (dq * (x - tip[1]) / r, dq * (y - tip[2]) / r)
end

J_fem = fracture_j_integral(m, E, ν, fixed, loads, ones(size(m.faces, 2)), qgrad)
K_fem = sqrt(J_fem * E)                   # plane stress: J = K²/E

a = 0.1
F = 1.12 - 0.231 * a + 10.55 * a^2 - 21.72 * a^3 + 30.39 * a^4
K_ana = σfar * sqrt(π * a) * F
println("K_fem = ", round(K_fem / 1e6; digits = 2), " MPa√m,  K_ana = ",
        round(K_ana / 1e6; digits = 2), " MPa√m  (", round(100 * K_fem / K_ana; digits = 1), "%)")

K_fem = 25.67 MPa√m,  K_ana = 26.54 MPa√m  (96.7%)

The FEM J, converted to K through J = K²/E, reproduces the handbook stress-intensity factor σ√(πa)·F(a/W) to about 3% on this mesh — and refines toward it. At 40 MPa the plate is safe: K = 25.7 MPa√m is below the 29 MPa√m toughness of 7075-T6, though with only 1.1× margin — crack inspection would be mandatory.

Next steps

• Phase-field fracture — a regularised crack without a sharp tip.
• Linear elasticity — the stress analysis that feeds the J-integral.
• Materials — the built-in library that provides K_Ic, E, and σy.