RasmahRasmah

Cam-Clay Critical-State Soil Plasticity

Overview

Soils and “wet” clays fail gradually: they harden as they are compressed, but under sustained shear they approach a critical state at which they flow at constant volume. The modified Cam-clay model (Roscoe & Burland 1968; Wood 1990) captures this with a single yield surface in the mean–deviatoric stress invariants

\[f = q^2 + M^2 p\,(p - p_c) = 0\]

where p is the mean effective stress, q the deviatoric stress, M the critical-state stress ratio, and p_c the preconsolidation pressure (the hardening variable). An associated flow rule and a volumetric hardening law dp_c = (p_c/λ̃) dε_v^p (with λ̃ = λ − κ, the compression and swelling slopes) close the model. The key prediction: as dε_s grows, the stress ratio q/p → M and the volumetric strain stops changing — the critical state.

Rasmah integrates the Cam-clay return map as a scalar Newton solve on the plastic multiplier.

using Rasmah

A triaxial compression test

A single CamClayModel stores the current (p, q, p_c). Drive it through a sequence of volumetric + deviatoric strain increments — a drained triaxial test — and watch the stress path and the hardening variable evolve.

soil = CamClayModel(; M=1.2, λ=0.2, κ=0.05, K=1e7, G=1e7, p0=1e5, pc0=2e5)

dεs = 0.001
n = 300
ps = Float64[]; qs = Float64[]; ev = Float64[]
for _ in 1:n
    camclay_step!(soil, 0.0, dεs)          # undrained-per-step shear
    push!(ps, soil.p); push!(qs, soil.q)
end
println("final q/p = ", round(qs[end] / ps[end]; digits = 4), "  (critical-state M = 1.2)")
final q/p = 1.2  (critical-state M = 1.2)

The stress ratio approaches the critical-state value M:

println("M recovered = ", round(qs[end] / ps[end]; digits = 4), "  (M = 1.2)")
println("q = ", round(qs[end] / 1e3; digits = 1), " kPa")
M recovered = 1.2  (M = 1.2)
q = 120.0 kPa

Verification — the critical state line

The model must reproduce the two textbook signatures: (1) the stress ratio converges to M, and (2) the mean effective stress settles onto the critical state line q = M p (passing through the point where p stops evolving). Check that the final state satisfies q = M p:

M = 1.2
q_crit = M * ps[end]
println("q final      = ", round(qs[end]; digits = 1))
println("M·p final    = ", round(q_crit; digits = 1))
println("q/(M·p) → 1 : ", round(qs[end] / q_crit; digits = 4))
q final      = 120000.0
M·p final    = 120000.0
q/(M·p) → 1 : 1.0

Next steps