MeanFieldHomogenization.jl

Appearance

Two routes to a crack: the COD tensor and the flattened ellipsoid ​

The whole page as one script

Every block of this page, in order, is scripts/17_crack_routes.jl. From a clone of the repository, julia --project=docs -e 'using Pkg; Pkg.instantiate()' once, then julia scripts/17_crack_routes.jl, runs it end to end.

Before this page

The theory page Crack opening displacement and compliance: §2 defines the COD tensor , normalized by the half-width of the crack, and §3 proves that the compliance contribution of an elliptical crack is       , with   the flatness of the ellipsoid it is the limit of. Every tensor on this page is in that normalization; the Conventions give the factor   that separates it from papers normalizing by .

The identity of §3 joins two independent computations of a crack, and this page runs both, first in an isotropic matrix where everything is a closed form, then in matrices drawn at random, to show that they meet.

  • From the crack plane. cod_tensor integrates the reduced kernel of the Green operator along the crack normal and around the crack front (Kanaun and Levin, 2009; Barthélémy, 2026): it returns , and follows from the identity (compliance_from_cod).

  • From a flattened ellipsoid. The crack is the limit of the ellipsoid of semi-axes as   (Barthélémy, 2009). The first-order term of its Hill tensor in gives by a computation on blocks of a matrix, and follows from the identity read the other way (cod_from_compliance).

julia
using MeanFieldHomogenization
using TensND
using LinearAlgebra
using ForwardDiff
using Random
using Printf
import SymPyPythonCall
import DECUHR, Integrals

const canonical = TensND.CanonicalBasis{3, Float64}()
const tol = (; abstol = 1.0e-12, reltol = 1.0e-10)
relerr(A, B) = norm(components_canon(A) - components_canon(B)) / norm(components_canon(B))
relerr (generic function with 1 method)

The flattened ellipsoid to first order ​

The Hill tensor of the ellipsoid of semi-axes , whose third axis is the crack normal , expands as

so that the second Hill tensor      expands as    with     . The leading term is singular: every strain   lies in its kernel, three dimensions which, in the Kelvin–Mandel components of the crack frame , are the components , and . In these blocks,

with the block of on the kernel of (Barthélémy, 2009). Only    is needed. The function below takes it as a Kelvin–Mandel matrix, in whatever arithmetic produced it, and returns in the crack frame .

julia
function H_from_first_order(P₁, C, ℬ)
    Q₁ = -mandel(C) * P₁ * mandel(C)
    Y = mandel(from_mandel(Q₁, canonical), ℬ)
    H = zeros(eltype(Y), 6, 6)
    H[3:5, 3:5] = inv(Y[3:5, 3:5])
    return from_mandel(H, ℬ)
end
H_from_first_order (generic function with 1 method)

Evaluating   at a small but finite flatness is the obvious alternative, and the theory page does it: its error decreases only like (the check of ), and two flatnesses are combined to cancel that first-order error (every route in a triclinic matrix). The block computation removes it at the source, with no flatness to choose.

In an isotropic matrix, exactly ​

With of Young's modulus and Poisson's ratio , both routes are closed forms, and SymPy carries them exactly. The COD tensor of a penny-shaped crack of radius is Sneddon's:

and cod_tensor returns it exactly (the difference simplifies to zero):

julia
E, ν, ρ = SymPyPythonCall.symbols("E ν ρ", positive = true)
C_sym = TensISO{3}(E / (1 - 2ν), E / (1 + ν))
B_sym = get_array(cod_tensor(PennyCrack(1), C_sym))
B_sneddon = 16 * (1 - ν^2) / (3 * SymPyPythonCall.PI * E) * Diagonal([2 / (2 - ν), 2 / (2 - ν), 1])
SymPyPythonCall.simplify.(B_sym - B_sneddon)

The flattened route takes the Hill tensor of the oblate spheroid , a closed form in , differentiates it with respect to and lets  . The block computation then gives

which is     for the above, exactly, as §3 says:

julia
P_sym = mandel(hill_tensor(Ellipsoid(1, 1, ρ), C_sym))
P₁_sym = [SymPyPythonCall.limit(SymPyPythonCall.diff(p, ρ), ρ => 0) for p in P_sym]
H_sym = H_from_first_order(P₁_sym, C_sym, canonical)
SymPyPythonCall.simplify.(mandel(H_sym) - mandel(compliance_from_cod(Tens(B_sym), PennyCrack(1))))

In random anisotropic matrices ​

Each matrix below is a random symmetric positive definite Kelvin–Mandel matrix, fully anisotropic (21 independent constants), and each crack has a random orientation. is now the derivative, by ForwardDiff, of the Hill tensor computed by the nested cubature, taken at two small flatnesses and and combined so that its own first-order term cancels:    .

julia
function P₁_by_forwarddiff(a, b, R, C; ρ₀ = 1.0e-5)
    P(ρ) = mandel(hill_tensor(Ellipsoid(a, b, ρ * b, R), C; method = :nestedquadgk, tol...), canonical)
    P′(ρ) = ForwardDiff.derivative(P, ρ)
    return 2 * P′(ρ₀ / 2) - P′(ρ₀)
end

function random_stiffness(rng)
    G = randn(rng, 6, 6)
    return from_mandel(100 * (G * transpose(G) / 6 + I), canonical)
end

function random_rotation(rng)
    F = qr(randn(rng, 3, 3))
    R = Matrix(F.Q) * Diagonal(sign.(diag(F.R)))
    return det(R) > 0 ? R : -R
end
random_rotation (generic function with 1 method)

For each matrix and each aspect ratio  , the table compares from the flattening with from by the identity, and the reverse; then from the crack plane by its three integrators, the nested cubature being the reference: the residues and DECUHR.

julia
rng = Xoshiro(20261018)
@printf "%-6s %-5s  %-13s %-13s %-13s %-13s\n" "matrix" "η" "ℍ: flat/𝐁" "𝐁: flat/𝐁" "𝐁 residues" "𝐁 DECUHR"
for k in 1:4
    C = random_stiffness(rng)
    R = random_rotation(rng)
    for η in (1.0, 0.5, 0.2)
        crack = EllipticCrack(1.0, η, R)
        B = cod_tensor(crack, C; method = :nestedquadgk, tol...)
        H = H_from_first_order(P₁_by_forwarddiff(1.0, η, R, C), C, crack_basis(crack))
        @printf "%-6d %-5.2f  %-13.1e %-13.1e %-13.1e %-13.1e\n" k η relerr(
            H, compliance_from_cod(B, crack)
        ) relerr(cod_from_compliance(H, crack), B) relerr(
            cod_tensor(crack, C; method = :residues, tol...), B
        ) relerr(cod_tensor(crack, C; method = :decuhr, tol...), B)
    end
end
matrix η      ℍ: flat/𝐁     𝐁: flat/𝐁     𝐁 residues    𝐁 DECUHR
1      1.00   2.8e-09       2.9e-09       2.3e-15       5.0e-14
1      0.50   2.2e-10       2.6e-10       2.7e-15       1.4e-15
1      0.20   2.8e-08       2.2e-08       2.8e-15       2.4e-13
2      1.00   2.5e-09       2.6e-09       5.4e-15       6.0e-15
2      0.50   3.2e-10       3.2e-10       5.5e-15       2.2e-13
2      0.20   1.1e-08       9.8e-09       5.2e-15       1.6e-13
3      1.00   4.3e-10       4.7e-10       2.8e-15       4.0e-13
3      0.50   3.0e-10       3.5e-10       2.7e-15       1.5e-13
3      0.20   4.9e-09       4.0e-09       2.7e-15       2.4e-12
4      1.00   1.3e-09       1.2e-09       4.7e-15       7.7e-14
4      0.50   2.5e-10       2.2e-10       4.4e-15       2.3e-14
4      0.20   1.5e-09       1.6e-09       4.6e-15       1.9e-14

The two routes agree to between and  , whatever the anisotropy: the residual of the left in , larger for the narrower cracks. On the reduced kernel, the residues meet the nested cubature to rounding, and DECUHR, at the same tolerances, to   at worst: the nested cubature is the default, DECUHR a second opinion.

An analytic anchor in an anisotropic matrix ​

Neither route above is a closed form in an anisotropic matrix. One exists when the matrix is the image of an isotropic one by a linear map   ,       (elliptic orthotropy). The crack of radii maps to a crack of radii and normal in the isotropic matrix, whose is a closed form, and (Barthélémy, 2026) gives, in the normalization of the package,

julia
function transformed_crack(P, a, b, ℓ, m, n)
    g, h = P \ ℓ, P \ m
    area = norm(cross(g, h))
    s₊ = sqrt(a^2 * norm(g)^2 + b^2 * norm(h)^2 + 2a * b * area)
    s₋ = sqrt(a^2 * norm(g)^2 + b^2 * norm(h)^2 - 2a * b * area)
    ã, b̃ = (s₊ + s₋) / 2, (s₊ - s₋) / 2
    ñ = normalize(transpose(P) * n)
    w_ℓ = (a^2 * norm(g)^2 - b̃^2) * g + b^2 * dot(g, h) * h
    w_m = a^2 * dot(g, h) * g + (b^2 * norm(h)^2 - b̃^2) * h
    ℓ̃ = normalize(norm(w_ℓ) ≥ norm(w_m) ? w_ℓ : w_m)
    return ã, b̃, hcat(ℓ̃, cross(ñ, ℓ̃), ñ)
end

Pmap = [1.3 0.2 -0.1; 0.15 0.8 0.25; -0.05 0.3 1.1]
PP = mandel(Tens(Pmap) ⊠ˢ Tens(Pmap))
C̃ = TensISO{3}(3 * 20.0, 2 * 9.0)
C_EO = from_mandel(PP * mandel(C̃) * transpose(PP), canonical)
a, b = 1.0, 0.45
R = random_rotation(Xoshiro(2))
ℓ, m, n = eachcol(R)
ã, b̃, R̃ = transformed_crack(Pmap, a, b, ℓ, m, n)
B̃ = components_canon(cod_tensor(EllipticCrack(ã, b̃, R̃), C̃))
B_rule = (b̃ / b) / norm(transpose(Pmap) * n) * (transpose(Pmap) \ B̃ / Pmap)
crack = EllipticCrack(a, b, R)
H_EO = H_from_first_order(P₁_by_forwarddiff(a, b, R, C_EO), C_EO, crack_basis(crack))
for (label, B) in (
        ("crack plane, nested cubature", cod_tensor(crack, C_EO; method = :nestedquadgk, tol...)),
        ("crack plane, residues", cod_tensor(crack, C_EO; method = :residues, tol...)),
        ("crack plane, DECUHR", cod_tensor(crack, C_EO; method = :decuhr, tol...)),
        ("flattened ellipsoid", cod_from_compliance(H_EO, crack)),
    )
    @printf "%-30s %.1e\n" label norm(components_canon(B) - B_rule) / norm(B_rule)
end
crack plane, nested cubature   3.8e-15
crack plane, residues          8.7e-16
crack plane, DECUHR            2.8e-13
flattened ellipsoid            1.4e-10

Every route meets the closed form: the reduced kernel to rounding, the flattened ellipsoid to the accuracy of its first-order term. The same comparison, for transformed transversely isotropic matrices and ribbon-like cracks, is part of the test suite.

Where to go next ​

A single crack is all this page needs. Crack distributions: isotropic or parallel gives populations of them a density and an orientation rule, and follows the self-consistent estimates to percolation; Validating a finite-element crack computes the opening of one crack by finite elements and compares it with the closed-form .

This page was generated using Literate.jl.