Crack opening displacement and compliance

A crack has no volume, so it cannot be described by a volume fraction and a stiffness. It is described instead by how much it opens under a given remote stress. That is the crack opening displacement tensor $\boldsymbol{B}$, and everything else on this page follows from it.

\[\boldsymbol{B} \;\longrightarrow\; \mathbb{H} \;\longrightarrow\; \Delta\mathbb{S} \qquad\text{and}\qquad \boldsymbol{B} \;\longrightarrow\; \underline{K},\ \underline{N}\]

Read left to right: the COD tensor $\boldsymbol{B}$ (order 2) gives the compliance contribution $\mathbb{H}$ (order 4), which gives the dilute correction $\Delta\mathbb{S}$ to the effective compliance; independently, $\boldsymbol{B}$ gives the stress and displacement intensity factors $\underline{K}$, $\underline{N}$ at the crack front.

MeanFieldHom computes $\boldsymbol{B}$ first — the intensity factors need it, and rebuilding $\mathbb{H}$ from it avoids inverting a rank-deficient order-4 tensor. Its normalization is not unique in the literature, hence the Conventions section below.

Geometry

A flat crack is the limit of a flat spheroidal inclusion. Keeping the shape tensor of Notation,

\[\boldsymbol{A} = \underline{\ell}\otimes\underline{\ell} + \eta\,\underline{m}\otimes\underline{m} + \omega\,\underline{n}\otimes\underline{n}, \qquad \eta = \frac{b}{a}, \qquad \omega = \frac{c}{a}\to 0,\]

with $(\underline{\ell},\underline{m})$ the in-plane unit vectors along the major and minor semi-axes $a\ge b$, and $\underline{n}$ the unit normal. Two families matter, and they are genuinely different objects, not two regimes of one:

The COD tensor $\boldsymbol{B}$

Under a remote stress $\boldsymbol{\Sigma}$, the two crack faces separate by the displacement jump $[\![\underline{u}]\!]$. By linearity and the superposition principle, that jump depends on the loading only through the resolved traction $\boldsymbol{\Sigma}\cdot\underline{n}$ on the crack plane. Its average over the crack surface $\mathcal{I}$ defines $\boldsymbol{B}$:

\[\boxed{\; \frac{\bigl\langle [\![\underline{u}]\!] \bigr\rangle_{\mathcal{I}}}{b} = \boldsymbol{B}\cdot\boldsymbol{\Sigma}\cdot\underline{n} \;}\]

The normalization is by the in-plane half-width $b$ — the minor semi-axis of the ellipse, the half-width of the ribbon. This makes $\boldsymbol{B}$ size-independent: it depends on the crack shape (through $\eta$) and on its orientation, never on how big it is. This is the convention of [30], following [31], [32], and it is the one MeanFieldHom implements (cod_tensor, alias B_tensor).

The shape coefficient $\chi$

The opening profile is an ellipsoidal cap, so the average jump is a fixed fraction of the maximum jump $\underline{\beta}$ at the crack centre:

\[\bigl\langle [\![\underline{u}]\!] \bigr\rangle_{\mathcal{I}} = \chi\,\underline{\beta}, \qquad \chi^{\mathcal{E}} = \frac{1}{\pi a b}\!\!\int_{\frac{x^2}{a^2}+\frac{y^2}{b^2}\le 1}\!\! \sqrt{1-\frac{x^2}{a^2}-\frac{y^2}{b^2}}\;\mathrm{d}x\,\mathrm{d}y = \frac{2}{3}, \qquad \chi^{\mathcal{R}} = \frac{1}{2b}\!\int_{-b}^{b}\!\sqrt{1-\frac{y^2}{b^2}}\;\mathrm{d}y = \frac{\pi}{4}.\]

$\chi$ is worth introducing explicitly, because every numerical factor on this page comes from it. It is the one quantity that distinguishes the elliptic geometry from the ribbon one.

The 3-D → 2-D limit: a trap

A ribbon is the limiting shape of an ellipse as $\eta\to 0$, so one expects $\boldsymbol{B}^{\mathcal{R}}$ to be $\lim_{\eta\to 0}\boldsymbol{B}^{\mathcal{E}}$. It is not [30]:

\[\boxed{\; \boldsymbol{B}^{\mathcal{R}}(\underline{m},\underline{n}) = \frac{\chi^{\mathcal{R}}}{\chi^{\mathcal{E}}}\, \lim_{\eta\to 0}\boldsymbol{B}^{\mathcal{E}}(\underline{m},\underline{n},\eta) = \frac{3\pi}{8}\, \lim_{\eta\to 0}\boldsymbol{B}^{\mathcal{E}}(\underline{m},\underline{n},\eta). \;}\]

The reason is exactly $\chi$. The geometry does converge, and so does the pointwise opening profile; what changes is the relation between the average opening and the maximum one, because averaging over a shrinking ellipse is not the same operation as averaging across the width of a strip ($2/3 \ne \pi/4$). Since $\boldsymbol{B}$ is defined from the average, it picks up the ratio $\chi^{\mathcal{R}}/\chi^{\mathcal{E}} = 3\pi/8$.

This matters in practice: the intensity factors of an elliptic crack are computed from the ribbon tensor of the tangent ribbon at each front point (see Intensity factors at the crack front below), so both objects appear in the same formula and must not be confused.

From $\boldsymbol{B}$ to the compliance $\mathbb{H}$

The extra strain a crack contributes, per unit volume of the embedding ellipsoid, is the average of the displacement jump over the crack surface $S$, spread over the volume $V$:

\[\boldsymbol{\varepsilon}^{\text{extra}} = \frac{1}{V}\int_{S} [\![\underline{u}]\!] \stackrel{s}{\otimes}\underline{n}\,\mathrm{d}S = \frac{S}{V}\, \bigl\langle [\![\underline{u}]\!] \bigr\rangle_{\mathcal{I}} \stackrel{s}{\otimes}\underline{n} = \frac{S\,b}{V}\; \underline{n}\stackrel{s}{\otimes}\boldsymbol{B}\stackrel{s}{\otimes}\underline{n} :\boldsymbol{\Sigma},\]

using the definition of $\boldsymbol{B}$. Identifying this with $\mathbb{Q}^{-1}:\boldsymbol{\Sigma}$, where $\mathbb{Q} = \mathbb{C}-\mathbb{C}:\mathbb{P}:\mathbb{C}$ is the second Hill tensor, gives $\mathbb{Q}^{-1} = (Sb/V)\,\underline{n}\stackrel{s}{\otimes}\boldsymbol{B} \stackrel{s}{\otimes}\underline{n}$, and hence

\[\boxed{\; \mathbb{H} \;=\; \lim_{c/b\,\to\,0}\ \frac{c}{b}\,\mathbb{Q}^{-1} \;=\; \frac{c\,S}{V}\; \underline{n}\stackrel{s}{\otimes}\boldsymbol{B}\stackrel{s}{\otimes}\underline{n}. \;}\]

The limit is finite: the components $(\mathbb{Q}^{-1})_{nijk}$ diverge like $b/c$, so $(c/b)\,\mathbb{Q}^{-1}$ stays bounded. Evaluating the geometric factor $cS/V$ for each family gives the two numbers used in the code:

Both factors are independent of $\eta$, which is the point of normalizing the limit by $b$: one definition covers the 2-D and the 3-D geometry, and $\mathbb{H}$ inherits $\boldsymbol{B}$'s size-independence.

The bridge in both directions is compliance_from_cod / cod_from_compliance (src/Cracks/cod_H_bridge.jl), dispatching on the crack type to apply the correct factor.

Conventions: what $\boldsymbol{B}$ and $\mathbb{H}$ are normalized by

The $3/4$ above looks like it disagrees with the literature. It does not — the three sources normalize the limit differently, while agreeing on $\boldsymbol{B}$. Naming them once avoids a great deal of confusion:

So the two elliptic compliances differ by exactly $\eta$:

\[\mathbb{H}^{\mathcal{E}}_{\texttt{MeanFieldHom}} = \frac{1}{\eta}\;\mathbb{H}^{\mathcal{E}}_{\text{Echoes}}, \qquad \mathbb{H}^{\mathcal{R}}_{\texttt{MeanFieldHom}} = \mathbb{H}^{\mathcal{R}}_{\text{Echoes}} .\]

They coincide for the penny crack $\eta=1$, which is why the discrepancy is invisible on the most common test case — and why a penny-only cross-check cannot detect a wrong $\eta$-dependence.

Measured, not inferred from the papers — Echoes' crack_compliance against compliance_contribution on a flat triaxial ellipsoid ($a = 1$, $b = \eta$, $c = 10^{-5}$):

On the MeanFieldHom side the $3/4$ is $\eta$-independent to machine precision ($\mathbb{H}_{3333}/B_{33} = 0.750000$ for every $\eta$).

Checking $\mathbb{H}$ without going through $\boldsymbol{B}$

The closed forms of the next section give $\boldsymbol{B}$, and $\mathbb{H}$ follows from it by the $3/4$. Comparing $\mathbb{H}$ to those closed forms would therefore verify nothing: it is the same formula read twice.

Read from right to left the boxed definition is a recipe: take a flat ellipsoid, compute $\mathbb{P}$, assemble $\mathbb{Q} = \mathbb{C}-\mathbb{C}:\mathbb{P}:\mathbb{C}$, invert, scale by the flatness $\omega = c/b$. That path never touches the crack machinery, so agreement with compliance_contribution is evidence rather than tautology:

using MeanFieldHom, TensND, Plots
gr()  # headless backend; GKSwstype is set to "100" in make.jl

E, ν = 210.0, 0.3
C₀ = TensISO{3}(E / (1 - 2ν), E / (1 + ν))   # (3k, 2μ)

## ω-family of the definition, evaluated at finite flatness
function H_from_hill(a, b, ω, C₀)
    P = hill_tensor(Ellipsoid(a, b, ω * b), C₀)
    Q = C₀ - C₀ ⊡ P ⊡ C₀
    return ω * inv(Q)
end

relerr(A, B) =
    maximum(abs(A[i,j,k,l] - B[i,j,k,l]) for i in 1:3, j in 1:3, k in 1:3, l in 1:3) /
    maximum(abs(B[i,j,k,l]) for i in 1:3, j in 1:3, k in 1:3, l in 1:3)

ωs = exp10.(range(-1, -3; length = 9))
plt = plot(; xscale = :log10, yscale = :log10, legend = :topleft,
           xlabel = "flatness ω = c/b", ylabel = "relative error on ℍ")
for η in (1.0, 0.5, 0.3)
    crack = η == 1.0 ? PennyCrack(1.0) : EllipticCrack(1.0, η)
    H_ref = compliance_contribution(crack, C₀)
    plot!(plt, ωs, [relerr(H_from_hill(1.0, η, ω, C₀), H_ref) for ω in ωs];
          marker = :circle, label = "η = $η")
end
plot!(plt, ωs, ωs; linestyle = :dash, color = :black, label = "slope 1")
plt

The three curves are straight lines of slope 1: the error decays like $\omega$, which is the order of the Taylor term that resolves the limit ([34]). That slope is the real content of the check — a single $\omega$ would not distinguish a true limit from a coincidence.

Closed forms of $\boldsymbol{B}$

Isotropic matrix

For $\mathbb{C} = 3k\,\mathbb{J}+2\mu\,\mathbb{K}$ with Young's modulus $E$ and Poisson ratio $\nu$, in the crack frame $(\underline{\ell},\underline{m},\underline{n})$:

\[\boldsymbol{B} = B_{nn}\,\underline{n}\otimes\underline{n} + B_{mm}\,\underline{m}\otimes\underline{m} + B_{\ell\ell}\,\underline{\ell}\otimes\underline{\ell},\]

\[\begin{aligned} B_{nn} &= \frac{8\,\eta\,(1-\nu^{2})}{3E}\,\frac{1}{\mathcal{E}_\eta},\\[4pt] B_{mm} &= \frac{8\,\eta\,(1-\nu^{2})}{3E}\, \frac{1-\eta^{2}} {\bigl(1-(1-\nu)\eta^{2}\bigr)\mathcal{E}_\eta - \nu\,\eta^{2}\,\mathcal{K}_\eta},\\[4pt] B_{\ell\ell} &= \frac{8\,\eta\,(1-\nu^{2})}{3E}\, \frac{1-\eta^{2}} {\bigl(1-\nu-\eta^{2}\bigr)\mathcal{E}_\eta + \nu\,\eta^{2}\,\mathcal{K}_\eta}, \end{aligned}\]

with $\mathcal{K}_\eta = \mathcal{K}(\sqrt{1-\eta^{2}})$ and $\mathcal{E}_\eta = \mathcal{E}(\sqrt{1-\eta^{2}})$ the complete elliptic integrals of the first and second kind [11], provided by ell_K and ell_E.

The three components are the three fracture modes: $B_{nn}$ opens the crack (mode I), $B_{mm}$ and $B_{\ell\ell}$ shear it (modes II and III).

Penny crack $\eta=1$, where $\mathcal{K}_1 = \mathcal{E}_1 = \pi/2$:

\[B_{nn} = \frac{16\,(1-\nu^{2})}{3\pi E}, \qquad B_{mm} = B_{\ell\ell} = \frac{B_{nn}}{1-\nu/2}.\]

Implementation: src/Cracks/cod_analytical.jl, selected by method = :auto when $\mathbb{C}_0$ is a TensISO{4,3}.

Transversely isotropic matrix

When the matrix is transversely isotropic with its axis aligned with the crack normal $\underline{n}$, $\boldsymbol{B}$ is still analytical. The closed forms use the engineering parameters $(E,\nu_1,\nu_2,H,\Gamma)$ defined on the compliance $\mathbb{S} = \mathbb{C}^{-1}$ [35], [36], [37], and reduce to the isotropic case for $\nu_1=\nu_2=\nu$, $H=\Gamma=1$. The auxiliary coefficients are documented inline in src/Cracks/cod_analytical.jl.

The more general cases — a TI axis not aligned with the crack normal, or an elliptic-orthotropic matrix — are treated in [33] and [30] but are not yet exposed here.

Arbitrary anisotropy — numerical

No closed form exists in general. Following [34], the limit $\omega\to 0$ is resolved by extracting the first-order term of the Taylor expansion of $\mathbb{P}$ in $\omega$; that term has an integral representation on the unit circle of the crack plane, evaluated by either algorithm trait:

• DECUHR — adaptive cubature [15], ForwardDiff-safe (src/Cracks/green_decuhr.jl);
• Residue — Cauchy-residue reduction to a 1-D quadrature, as in [14] adapted to the crack kernel, Float64 only (src/Cracks/green_residue.jl).

method = :auto picks Residue for anisotropic Float64 inputs and falls back to DECUHR for symbolic or ForwardDiff.Dual scalars.

Dilute correction to the effective compliance

compliance_contribution(crack, C₀) returns $\mathbb{H}$ itself — the size-independent contribution, not the dilute correction. Cracks have no volume fraction, so the amount of cracking is measured by a Budiansky crack density [38], [32], and reintroduced by delta_compliance:

Implementation: src/Cracks/compliance.jl, dispatching on the crack shape.

Intensity factors at the crack front

At a point $\underline{x}^{\star}_{0}$ of the crack front, with in-plane outer normal $\underline{\nu}$ and tangent $\underline{\tau} = \underline{n}\wedge\underline{\nu}$, the asymptotic expansions of the jump and the traction read [39], [40], [41]:

\[[\![\underline{u}]\!](\underline{x}^{\star}_{0}+r\underline{\nu}) \underset{r\to 0^{-}}{\sim} 8\sqrt{\frac{-r}{2\pi}}\;\underline{N}, \qquad \underline{t}(\underline{x}^{\star}_{0}+r\underline{\nu}) \underset{r\to 0^{+}}{\sim} \frac{\underline{K}}{\sqrt{2\pi r}} .\]

$\underline{N}$ is the displacement intensity factor (DIF) and $\underline{K}$ the stress intensity factor (SIF), normalized so that the local energy release rate is simply $G = \underline{K}\cdot\underline{N}$ [42], [43].

The central result of the anisotropic theory [44], [45], [36] is that SIF and DIF are purely local and are exchanged by the COD tensor of the ribbon crack tangent to the real crack at the observation point:

\[\boxed{\; \underline{K} = \pi\,\bigl(\boldsymbol{B}^{\mathcal{R}}(\underline{\nu},\underline{n})\bigr)^{-1} \cdot\underline{N} \;}\]

This holds whatever the matrix anisotropy and whatever the remote loading — and it is why $\boldsymbol{B}^{\mathcal{R}}$, with its $3\pi/8$ factor, cannot be dispensed with even when studying elliptic cracks.

Elliptic crack

Parametrize the front by an angle $\theta_y$:

\[\underline{y}^{\star}_{0} = \cos\theta_y\,\underline{\ell} + \sin\theta_y\,\underline{m}, \qquad \underline{\nu} = \frac{\boldsymbol{S}^{\dagger}\!\cdot\underline{y}^{\star}_{0}} {\|\boldsymbol{S}^{\dagger}\!\cdot\underline{y}^{\star}_{0}\|},\]

with $\boldsymbol{S}$ the in-plane semi-axis tensor and $\boldsymbol{S}^{\dagger}$ its pseudo-inverse. Then

\[\underline{N}^{\mathcal{E}} = \tfrac{3}{8}\sqrt{\pi b}\;\sqrt{\varrho}\; \boldsymbol{B}^{\mathcal{E}}(\underline{m},\underline{n},\eta) \cdot\boldsymbol{\Sigma}\cdot\underline{n}, \qquad \underline{K}^{\mathcal{E}} = \tfrac{3}{8}\pi^{3/2}\sqrt{b}\;\sqrt{\varrho}\; \bigl(\boldsymbol{B}^{\mathcal{R}}(\underline{\nu},\underline{n})\bigr)^{-1} \cdot\boldsymbol{B}^{\mathcal{E}}(\underline{m},\underline{n},\eta) \cdot\boldsymbol{\Sigma}\cdot\underline{n},\]

where the dimensionless front factor is

\[\varrho = b\,\|\boldsymbol{S}^{\dagger}\!\cdot\underline{y}^{\star}_{0}\| = \sqrt{\eta^{2}\cos^{2}\theta_y + \sin^{2}\theta_y}.\]

Ribbon crack

Here $\underline{\nu} = \pm\underline{m}$ and $\varrho = 1$:

\[\underline{N}^{\mathcal{R}} = \sqrt{\frac{b}{\pi}}\; \boldsymbol{B}^{\mathcal{R}}(\underline{m},\underline{n}) \cdot\boldsymbol{\Sigma}\cdot\underline{n}, \qquad \underline{K}^{\mathcal{R}} = \sqrt{\pi b}\;\boldsymbol{\Sigma}\cdot\underline{n}.\]

The SIF of an infinite ribbon crack is independent of the matrix stiffness — the $\boldsymbol{B}^{\mathcal{R}}$ of the DIF cancels against its inverse in the exchange relation.

Modes

\[K_{I} = |\underline{K}\cdot\underline{n}|, \qquad K_{II} = |\underline{K}\cdot\underline{\nu}|, \qquad K_{III} = |\underline{K}\cdot\underline{\tau}| .\]

Evaluation: sif and dif (src/Cracks/sif.jl).

Dispatch

Entry points: cod_tensor / B_tensor for $\boldsymbol{B}$, compliance_contribution for $\mathbb{H}$, delta_compliance for $\Delta\mathbb{S}$, sif / dif for the front quantities. The transport counterpart — a scalar COD and a rank-1 resistivity contribution — is treated in Thermal cracks, with the same geometric factors $3/4$ and $2/\pi$.