The strength of a sandstone: crushing grains and sliding contacts

Before this page

Friction of a granular medium from its contacts, the same polycrystal with rigid grains, and Layered sphere §7 for the strain carried by the opening of an interface.

A sandstone is a sand whose grains have been cemented together. Its strength involves two mechanisms: the grains can yield, and the cement between them can slide. This page reproduces the micromechanical model of Dormieux et al. (2010), in which the solid of the grains obeys a von Mises criterion,   with  , and the cemented contacts a Mohr–Coulomb one,    . The question is the macroscopic strength, in the plane of the mean stress    and of the deviatoric stress  .

The normalization of this page

This page keeps the normalization of Dormieux et al. (2010), with under the root, also for the deviatoric strain below, so that its formulas and figures are those of the article. The norm used in the rest of this documentation carries no factor (see Conventions): it is times the , and of this page. In that norm the grains obey   , and every ratio of this page, and the slope of the Drucker–Prager limit, is multiplied by .

The model follows the modified secant method: each nonlinear behavior is replaced by a linear one whose stiffness depends on the average strain it undergoes, and the strength is reached when these stiffnesses vanish. The linear problem behind it is therefore the heart of the model: an incompressible grain of shear modulus , bonded by an interface rigid in the normal direction and of tangential stiffness , in a self-consistent polycrystal with spherical pores of volume fraction .

1. The grain and its interface

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

brick(μs, kt) = LayeredSphere(
    (1.0,), (iso_stiffness(Inf, μs),);        # R = 1, incompressible grain
    interfaces = (SpringInterface(Inf, kt),), # rigid normally, kₜ tangentially
)

The average strain of the grain is read on the outer lip of its interface (their eq. 43): it includes the sliding of the contacts, which is what the concentration tensor of a LayeredSphere counts by default. With  , the behavior depends on the dimensionless contact stiffness   and on .

2. The self-consistent moduli in closed form

The self-consistent estimate      solves         , with the concentration tensors of the grain and of the pore computed in . Kept general, with a compressible grain and a compliant normal contact, the spherical equation is linear in :

using SymPyPythonCall
@syms k::positive μ::positive ks::positive μs::positive kn::positive kt::positive
@syms R::positive φ::positive x::positive κ::positive ϵ::positive

C = iso_stiffness(k, μ)
S_sym = LayeredSphere((R,), (iso_stiffness(ks, μs),); interfaces = (SpringInterface(kn, kt),))
αS, βS = get_data(gradient_gradient_concentration(S_sym, C, C))
aS, bS = get_data(flux_gradient_concentration(S_sym, C, C))
αP, βP = get_data(gradient_gradient_concentration(Ellipsoid(1.0), iso_stiffness(0, 0), C))

k_sc = simplify(solve(3k * ((1 - φ) * αS + φ * αP) - (1 - φ) * aS, k)[1])

which is eq. (49) of the article. Substituted in the deviatoric equation, it leaves a product of positive factors and of the cubic of eq. (50), the only factor of degree three in :

shear = 2μ * ((1 - φ) * βS + φ * βP) - (1 - φ) * bS
num, _ = sympy.fraction(sympy.together(subs(shear, k => k_sc)))
[sympy.degree(f, μ) for (f, _) in sympy.factor_list(num)[2]]

The article then lets the grain and the normal contact become incompressible,  . With   ,    and  , the lowest order in of the deviatoric equation is a positive multiple of the quadratic of eq. (52),            :

function lowest_order(e)
    n, _ = sympy.fraction(sympy.together(e))
    c = sympy.Poly(n, ϵ).all_coeffs()
    return c[findlast(c -> c != 0, c)]
end
eq52 = lowest_order(subs(shear, k => k_sc, ks => 1 / ϵ, kn => 1 / ϵ))
eq52 = subs(eq52, μ => x * μs, kt => κ * μs / R)
q52 = 16 * (5 + κ) * (3 - φ) * x^2 + ((9 + 77φ) * κ + 114 * (3φ - 1)) * x + 57κ * (2φ - 1)
factor(eq52 / q52)

Its positive root is  , so that   and, from the limit of (49),    (eqs. 51, 53). Perfect contacts ( ) give back the self-consistent porous polycrystal,     (eq. 54): as grows, the quadratic is dominated by its coefficient of , which vanishes at that value.

simplify(subs(sympy.Poly(q52, κ).coeff_monomial(κ), x => 3 * (1 - 2φ) / (3 - φ)))

3. The numerical estimate

The self-consistent iteration keeps every modulus above times its starting value, a guard against the null fixed point; above   the moduli fall to the order of as  , so the iteration starts from that scale when it is the smaller.

function moduli(μs, kt, φ)
    rve = RVE()
    add_phase!(rve, :S, brick(μs, kt), Dict(:C => iso_stiffness(1.0, 1.0)); fraction = 1 - φ)
    add_phase!(rve, :P, Ellipsoid(1.0), Dict(:C => iso_stiffness(0.0, 0.0)); fraction = φ)
    start = ForwardDiff.value(min(μs, kt))    # the scale of the moduli, also under a Dual
    sc = SelfConsistent(init = iso_stiffness(start, start), abstol = 1.0e-15, reltol = 1.0e-13, maxiters = 2000)
    return collect(k_mu(homogenize(rve, sc, :C)))
end

M52(κ, φ) = (a = 16 * (5 + κ) * (3 - φ); b = (9 + 77φ) * κ + 114 * (3φ - 1); c = 57κ * (2φ - 1);
    (-b + sqrt(b^2 - 4a * c)) / (2a))
for (κv, φv) in ((0.1, 0.35), (1.0, 0.25), (10.0, 0.30))
    kv, μv = moduli(1.0, κv, φv)
    @printf("κ = %4.1f, φ = %.2f: μ/μs = %.8f (eq. 52: %.8f), k/μ = %.6f (eq. 51: %.6f)\n",
        κv, φv, μv, M52(κv, φv), kv / μv, 4 * (1 - φv) / (3φv))
end

κ =  0.1, φ = 0.35: μ/μs = 0.06999429 (eq. 52: 0.06999429), k/μ = 2.476190 (eq. 51: 2.476190)
κ =  1.0, φ = 0.25: μ/μs = 0.32903820 (eq. 52: 0.32903820), k/μ = 4.000000 (eq. 51: 4.000000)
κ = 10.0, φ = 0.30: μ/μs = 0.40057945 (eq. 52: 0.40057945), k/μ = 3.111111 (eq. 51: 3.111111)

4. From the averages to the strength

Two averages drive the criteria, and both come from derivatives of the effective moduli at fixed load (their eqs. 23 and 33):

• the quadratic mean of the deviatoric strain in the grains,        with    and   ;

• the quadratic mean of the tangential traction on the contacts,         with    and   .

The mean normal traction on the contacts is     (eq. 38), whatever the scheme. The four derivatives are taken by automatic differentiation through the self-consistent iteration, at  :

function coefficients(κv, φv)
    k, μ = moduli(1.0, κv, φv)
    J = ForwardDiff.jacobian(p -> moduli(p[1], p[2], φv), [1.0, κv])
    return J[1, 1] / k^2, J[2, 1] / μ^2, κv^2 * J[1, 2] / k^2, κv^2 * J[2, 2] / μ^2
end

Grains alone. When the grains yield, their secant shear modulus behaves as  , and the first average gives       with and evaluated at  : the moduli are proportional to , so the criterion does not depend on it. With perfect contacts this is the ellipse of the von Mises porous solid (eq. 75),       .

Grains and contacts. When the contacts slide too, the Mohr–Coulomb law holds on the averages,    , and is no longer fixed: it parameterizes the envelope (eq. 85). Eliminating gives the ratio  , and the grain criterion then gives :

for  , wherever  . Two ends close the curve. As   the contacts have no tangential stiffness left; above   the curve then tends to a finite point, and the contacts alone fail along the Drucker–Prager line through the origin whose slope is that limit of (their eq. 81), the first part of the envelope. Below ,   and the curve starts on the deviatoric axis. At the other end, either vanishes at a finite and the curve comes down to the mean-stress axis, or it stays positive and the curve meets the ellipse as  . The first end is located by bisection between two values of the grid.

κs = 10.0 .^ range(-8, 6; length = 141)
H²(c, φv, αv) = (3αv^2 / ((1 - 2φv) * (1 - φv)^2) - c[3]) / c[4]
function point(c, φv, αv)
    h² = max(H²(c, φv, αv), 0.0)
    Sm = sqrt((1 - φv) / (c[1] + c[2] * h²))
    return (Sm, Sm * sqrt(h²))
end
function envelope(φv, αv, coefs)
    pts = φv > 1 / 3 ? [(0.0, 0.0)] : Tuple{Float64, Float64}[]   # contacts alone
    for i in eachindex(κs)
        if H²(coefs[i], φv, αv) ≥ 0
            push!(pts, point(coefs[i], φv, αv))
        elseif i > 1                       # H² vanishes between κs[i-1] and κs[i]
            lo, hi = log10(κs[i - 1]), log10(κs[i])
            for _ in 1:40
                mid = (lo + hi) / 2
                H²(coefficients(10.0^mid, φv), φv, αv) ≥ 0 ? (lo = mid) : (hi = mid)
            end
            push!(pts, point(coefficients(10.0^lo, φv), φv, αv))
            break
        end
    end
    return pts
end
ellipse_radii(φv) = (c = coefficients(1.0e8, φv); (sqrt((1 - φv) / c[1]), sqrt((1 - φv) / c[2])))
plt = plot(; aspect_ratio = :equal, xlabel = "−Σm / K", ylabel = "Σd / K", legend = :topright,
    title = "Strength of a sandstone (K = 1)")
for (i, φv) in enumerate((0.35, 0.25))
    a, b = ellipse_radii(φv)
    t = range(0, π / 2; length = 100)
    plot!(plt, a .* cos.(t), b .* sin.(t); color = :red, lw = 2,
        label = i == 1 ? "grains alone (eq. 75)" : "")
    coefs = [coefficients(κv, φv) for κv in κs]
    for (j, αv) in enumerate((0.16, 0.20, 0.25, 0.30))
        pts = envelope(φv, αv, coefs)
        plot!(plt, first.(pts), last.(pts); color = :blue, lw = 1.5, marker = :cross, ms = 2,
            label = (i == 1 && j == 1) ? "grains and contacts (eq. 85)" : "")
    end
end
annotate!(plt, [(0.35, 0.05, text("φ = 0.35", 9)), (1.0, 0.05, text("φ = 0.25", 9))])
plt

For each porosity, the red ellipse is the strength of the porous solid when only the grains yield; the blue curves, for   from the bottom up, are the strength when the contacts slide as well. The weaker the friction, the lower the curve. Below   a polycrystal whose contacts have lost all tangential stiffness still resists shear,   (their eq. 56), so the curves leave the deviatoric axis at a finite stress; above it does not, and the envelope starts from the origin, with the straight segment of the contacts failing alone, as in a sand. Frictional contacts keep the curve up to the ellipse, which it meets as  , where the contacts stop sliding; weakly frictional ones bring it down to the mean-stress axis first, so that a hydrostatic compression alone breaks the rock before the grains reach their own limit.

The friction coefficient enters only through , so the whole family of envelopes at   is one surface. Each envelope is built as in the plane figure, its end included, then resampled at points evenly spaced along its length, so that the edges of the surface are the envelopes and their ends rather than the cells of the grid:

include(joinpath(pkgdir(MeanFieldHomogenization), "scripts", "common", "docviz.jl"))
φv = 0.25
coef = [coefficients(κv, φv) for κv in κs]
## n points evenly spaced along the polyline pts.
function resample(pts, n)
    s = vcat(0.0, cumsum([hypot((pts[k] .- pts[k - 1])...) for k in 2:length(pts)]))
    return map(range(0, s[end]; length = n)) do t
        k = clamp(searchsortedlast(s, t), 1, length(pts) - 1)
        w = s[k + 1] > s[k] ? (t - s[k]) / (s[k + 1] - s[k]) : 0.0
        (1 - w) .* pts[k] .+ w .* pts[k + 1]
    end
end
αs = range(0.12, 0.34; length = 23)
envs = [envelope(φv, αv, coef) for αv in αs]
pts = resample.(envs, 61)
grid = [pts[i][j] for i in eachindex(pts), j in 1:61]
Xs, Zs = first.(grid), last.(grid)
Ys = [αv for αv in αs, _ in 1:61]
surface = surface_trace(Xs, Ys, Zs; color = "#2c7fb8", opacity = 0.9)
rim = line_trace(
    hcat(permutedims(hcat(Xs[1, :], Ys[1, :], Zs[1, :])), permutedims(hcat(Xs[:, end], Ys[:, end], Zs[:, end])),
        permutedims(hcat(reverse(Xs[end, :]), reverse(Ys[end, :]), reverse(Zs[end, :]))),
        permutedims(hcat(reverse(Xs[:, 1]), reverse(Ys[:, 1]), reverse(Zs[:, 1]))));
    color = FIG_VECTOR, width = 4)

Two families of planes cut it. A plane   gives back the envelope of that friction coefficient, as in the plane figure above. A plane   gives the strength in shear as a function of the friction at a fixed confinement, read on each envelope where it crosses that mean stress; it is empty where the rock breaks under that hydrostatic compression alone, whatever the shear.

xm, zm = maximum(Xs), maximum(Zs)
## Plane α = const: the envelope itself.
αc = 0.12:0.02:0.34
envc = [envelope(φv, αv, coef) for αv in αc]
friction = CutFamily(
    "Plane α = const", "α", string.(αc),
    [plane_patch([xm / 2, αv, zm / 2], [1, 0, 0], [0, 0, 1], 0.55xm, 0.55zm) for αv in αc],
    [vcat(first.(e)', fill(αv, 1, length(e)), last.(e)') for (αv, e) in zip(αc, envc)],
    [vcat(first.(e)', last.(e)') for e in envc],
    ("−Σm / K", "Σd / K"); active = 5,
)
## Plane −Σm = c: Σd where each envelope crosses that mean stress.
function crossing(e, c)
    for k in 1:(length(e) - 1)
        (x0, z0), (x1, z1) = e[k], e[k + 1]
        (x0 - c) * (x1 - c) ≤ 0 && x1 != x0 && return z0 + (c - x0) / (x1 - x0) * (z1 - z0)
    end
    return NaN
end
cs = 0.1:0.1:floor(xm; digits = 1)
shear = [[crossing(e, c) for e in envs] for c in cs]
confinement = CutFamily(
    "Plane Σm = const", "−Σm / K", string.(cs),
    [plane_patch([c, (first(αs) + last(αs)) / 2, zm / 2], [0, 1, 0], [0, 0, 1], 0.55 * (last(αs) - first(αs)), 0.55zm) for c in cs],
    [vcat(fill(c, 1, length(αs)), collect(αs)', d') for (c, d) in zip(cs, shear)],
    [vcat(collect(αs)', d') for d in shear],
    ("friction α", "Σd / K"); active = 5, equal = false,
)
section_figure([surface, rim], [friction, confinement];
    uid = "sandstone-envelopes-3d", height = 620, aspectmode = "cube",
    xlabel = "−Σm / K", ylabel = "friction α", zlabel = "Σd / K",
    title = "Strength envelopes of a sandstone, φ = 0.25")

The buttons above the surface choose a family of planes, or none, and the slider below moves the plane; on the right, the section is drawn in coordinates of its plane.

Where to go next

The rigid-grain limit of this model, with frictional contacts only, is Friction of a granular medium from its contacts. The secant reasoning on the energy derivatives is the same as in Quasi-brittle strength of cement paste and mortar, where it is applied to a brittle solid. The next application, Viscoelastic complex modulus of a bituminous mixture, turns to time-dependent behavior.