Hardening Soil (MN, Bricks)

This model extends the Hardening Soil (MN) model with the BRICK small-strain stiffness formulation of Cudny & Truty (2020)¹, itself based on the nested strain-space yield surface (multi-surface / "field of work-hardening moduli") concept of Simpson (1992)[@simpson1992]. All plasticity mechanisms — the Matsuoka–Nakai shear cone, the elliptical volumetric cap, the flow rules, the vertex cut-off and the return-mapping sequence — are identical to the base model and are not repeated here; this page only documents the additions. Where a symbol or equation is not redefined below, it carries the meaning given in the base model documentation.

Upcoming Release

The Hardening-Soil Bricks model will be available with the upcoming release.

Overview

The base Hardening Soil model uses a stress-dependent but otherwise constant unloading/reloading modulus \(E_{ur}\). Physically, soils are considerably stiffer than \(E_{ur}\) over a small range of recent strain history and degrade towards \(E_{ur}\) only as the strain amplitude grows (Atkinson & Sallfors 1991; Benz 2007²). The BRICK extension reproduces this behaviour by replacing the constant reference shear modulus \(G_{ur}^{\text{ref}}\) of the elastic operator with a strain-history-dependent tangent modulus \(G_{\text{ref},t}\) that

equals a much stiffer small-strain reference value \(G_0^{\text{ref}}\) immediately after any strain reversal, and
degrades smoothly, following a Hardin–Drnevich-type S-shaped curve, towards \(G_{ur}^{\text{ref}}\) as the strain moves away from the point of the last reversal.

The degradation is tracked using a discretised (multi-surface) analogue of the S-curve: NBRICKS = 10 nested yield surfaces ("bricks") in strain space, dragged behind a reference point (the "man") that represents the current total strain. This is the same mechanical analogy as the overlay/multi-surface plasticity models of Iwan and Mróz, here applied purely to the elastic stiffness rather than to a stress-space yield surface.

In addition to stiffening the elastic operator, the model enhances the hardening rate of both plasticity mechanisms while the material is inside the small-strain range (see Hardening Enhancement), following the approach used in the HS-Small model (Benz 2007²).

Strain-Space Bricks

Man, Bricks and String Lengths

The current total strain state is represented by a point in strain space, the man, denoted \(\boldsymbol{\varepsilon}^{\text{man}}\). Ten bricks, each attached to the man by a string of a fixed length \(L_j\) (\(j=1,\dots,10\)), are dragged along whenever the man moves far enough to pull a string taut. A brick is slack (not dragged) while the distance between the man and the brick's anchor point \(\boldsymbol{\varepsilon}^{\text{b}}_j\) is smaller than \(L_j\); it is active (dragged) once that distance would otherwise exceed \(L_j\).

For a strain (sub-)increment \(\Delta\boldsymbol{\varepsilon}^{\text{man}}\), each brick is updated independently:

\[ \boldsymbol{\varepsilon}^{\text{man}}_{n+1} = \boldsymbol{\varepsilon}^{\text{man}}_{n} + \Delta\boldsymbol{\varepsilon}^{\text{man}}, \qquad d_j = \left\lVert \boldsymbol{\varepsilon}^{\text{man}}_{n+1} - \boldsymbol{\varepsilon}^{\text{b}}_{j} \right\rVert_{\gamma}, \]

where \(\lVert \cdot \rVert_\gamma\) is the shear strain invariant defined below. If \(d_j > L_j\), brick \(j\) is dragged along the connecting direction until its string is taut again,

\[ \boldsymbol{\varepsilon}^{\text{b}}_{j} \leftarrow \boldsymbol{\varepsilon}^{\text{b}}_{j} + \frac{d_j - L_j}{d_j}\left(\boldsymbol{\varepsilon}^{\text{man}}_{n+1} - \boldsymbol{\varepsilon}^{\text{b}}_{j}\right); \]

otherwise \(\boldsymbol{\varepsilon}^{\text{b}}_{j}\) is left unchanged. Each brick \(j\) carries a fixed stiffness proportion \(\Delta\omega_j\) (see below); the actual reference tangent shear modulus after the update is

\[ G_{\text{ref},t} = G_0^{\text{ref}} \left(1 - \sum_{j:\ \text{active}} \Delta\omega_j\right), \]

i.e. every brick that gets dragged permanently removes its share of stiffness from \(G_0^{\text{ref}}\) until the man reverses direction far enough for that brick to go slack again. Immediately after a full reversal all strings are slack, no bricks are dragged, and \(G_{\text{ref},t} = G_0^{\text{ref}}\); after monotonic straining well beyond the degradation range, all ten bricks are dragged and \(G_{\text{ref},t} = G_{ur}^{\text{ref}}\).

Both the man strain and the brick anchor strains are stored using tensorial shear components (i.e. \(\varepsilon_{12}=\tfrac12\gamma_{12}\), not the engineering shear strain used elsewhere in dstrain); the Voigt component order is \(11,22,33,12,13,23\).

Shear Strain Invariant

The distance metric in strain space is the shear strain invariant

\[ \gamma = \sqrt{\varepsilon_{11}^2+\varepsilon_{22}^2+\varepsilon_{33}^2-\varepsilon_{11}\varepsilon_{22}-\varepsilon_{22}\varepsilon_{33}-\varepsilon_{11}\varepsilon_{33} + 3\left(\varepsilon_{12}^2+\varepsilon_{13}^2+\varepsilon_{23}^2\right)}, \]

equivalent to \(\gamma=\sqrt{\tfrac32\,\mathbf{e}:\mathbf{e}}\) with \(\mathbf{e}\) the deviatoric strain tensor. Being purely deviatoric, \(\gamma\) is insensitive to the volumetric part of the strain increment: a purely isotropic strain path never drags any brick. The radicand is clamped at zero as a numerical safeguard.

Discretisation of the Stiffness Degradation Curve

The string lengths \(L_j\) and stiffness proportions \(\Delta\omega_j\) are derived once per material point (in material_collect_properties) from the two HS-Small-type reference parameters \(\gamma_{0.7}\) and \(G_0^{\text{ref}}\), following the modified Hardin–Drnevich (Santos & Correia) secant degradation curve

\[ \frac{G_{\text{sec}}}{G_0} = \frac{1}{1+a\,\gamma/\gamma_{0.7}},\qquad a = \frac{3}{7}, \]

which by definition passes through \(G_{\text{sec}}/G_0=0.722\) at \(\gamma=\gamma_{0.7}\). The stiffness range

\[ \text{Rng} = \frac{G_0^{\text{ref}}-G_{ur}^{\text{ref}}}{G_0^{\text{ref}}} \]

is split into NBRICKS equal steps \(\Delta\text{Rng} = \text{Rng}/10\) (so \(\Delta\omega_j = \Delta\text{Rng}\) for every brick), and the string length of brick \(j\) is the tangent-curve inversion evaluated at the mid-height of its step:

\[ L_j = \frac{7}{3}\gamma_{0.7}\left(\sqrt{\dfrac{1}{1-j\,\Delta\text{Rng} + \tfrac12\Delta\text{Rng}}}-1\right), \qquad j=1,\dots,10. \]

For \(G_0^{\text{ref}} = G_{ur}^{\text{ref}}\), \(\text{Rng}=0\), all \(L_j\) and \(\Delta\omega_j\) vanish, and the BRICK extension degenerates exactly to the base model (see Degeneration to the Base Model).

Elastic Operator

The Young's modulus entering the elastic Jacobian is no longer the stress-dependent \(E_{ur}\) of the base model, but a strain- and stress-dependent modulus

\[ E_t = 2\,(1+\nu_{ur})\,G_{\text{ref},t}\,\left(\frac{\sigma_3^\ast+c\cot\varphi}{p_{\text{ref}}+c\cot\varphi}\right)^m, \qquad G_{ur}^{\text{ref}} \le G_{\text{ref},t} \le G_0^{\text{ref}}, \]

i.e. the barotropy factor of the base model is retained unchanged and simply applied to the current brick-degraded shear modulus instead of the fixed \(G_{ur}^{\text{ref}}=E_{ur}^{\text{ref}}/2(1+\nu_{ur})\). Poisson's ratio \(\nu_{ur}\) is assumed constant (only the shear modulus degrades; the bulk-to-shear ratio implied by \(\nu_{ur}\) is preserved at every degradation state).

The yield-function moduli remain barotropic

\(E_i\) and \(E_{ur}\) entering the shear yield function \(f_s\) (see base model) are not replaced by the brick-degraded modulus; they retain their standard stress-dependent (but strain-history-independent) definition, following Sect. 4.1 of Cudny & Truty (2020)¹. Only the elastic operator \(D^e\) (and, through the hardening enhancement below, the hardening rate) is affected by the BRICK state.

Sub-Incremented Predictor

Because \(G_{\text{ref},t}\) is a piecewise-constant function of strain (it changes only when a brick becomes active or slack), the elastic predictor is sub-incremented so that each sub-increment sees an (approximately) constant tangent modulus. The strain increment dstrain (converted to tensorial shear components) is split into

\[ n_{\text{sub}} = \max\!\left(1,\ \left\lceil \frac{\gamma(\Delta\boldsymbol{\varepsilon})}{\gamma_{\text{sub}}^{\max}} \right\rceil\right), \qquad \gamma_{\text{sub}}^{\max} = 5\times10^{-5}, \]

equal parts. For each sub-increment: the bricks are moved first (update_bricks), the resulting \(G_{\text{ref},t}\) is used to form \(E_t\), and the stress is advanced with the corresponding elastic Jacobian, \(\sigma \leftarrow \sigma + D^e(E_t):\Delta\boldsymbol{\varepsilon}_{\text{sub}}\). The barotropy factor \(\left(\frac{\sigma_3^\ast+c\cot\varphi}{p_{\text{ref}}+c\cot\varphi}\right)^m\) itself is frozen at the stress state passed in by the caller for the whole predictor call (consistent with the semi-explicit treatment of barotropy in the base model), and is only re-evaluated between the outer sub-increments of the base model's own friction-based sub-stepping (see Integration Algorithm Notes).

The stiffness ratio

\[ G_m = \min_{\text{loading history}} \left(\frac{G_{\text{ref},t}}{G_{ur}^{\text{ref}}}\right) \ge 1 \]

is updated after every brick sub-increment as \(G_m \leftarrow \min(G_m,\ G_{\text{ref},t}/G_{ur}^{\text{ref}})\). Being a running minimum, \(G_m\) can only decrease: once the material has been strained far enough to reach \(G_m=1\) (full degradation, all bricks dragged), a subsequent unloading–reloading cycle regains elastic stiffness (\(G_{\text{ref},t}\) rises again as bricks go slack) but does not reset \(G_m\), so the hardening enhancement below is not reapplied on a cycle that stays within a strain amplitude already visited. This mirrors the intent of the HS-Small hardening modification, which is meant to suppress spurious plastic softening only within the (not yet visited) small-strain range.

Hardening Enhancement

Both hardening laws of the base model are scaled by a factor

\[ H_i = G_m^{\ 1+\Large\frac{E_{ur}^{\text{ref}}}{2E_{50}^{\text{ref}}}}, \qquad 1 \le H_i \le \left(\frac{G_0^{\text{ref}}}{G_{ur}^{\text{ref}}}\right)^{1+\Large\frac{E_{ur}^{\text{ref}}}{2E_{50}^{\text{ref}}}}, \]

evaluated once per (sub-)increment from the current \(G_m\) (Eqs. 22–24 in Cudny & Truty 2020¹). \(H_i \ge 1\) increases the hardening rate of both mechanisms while the material is inside the small-strain range and reverts to \(H_i = 1\) once \(G_m = 1\) (fully degraded, base-model behaviour). Concretely:

Shear (cone) mechanism: the fixed hardening factor \(3/2\) of the base model (relating the cone's plastic multiplier to the increment of \(\gamma^p\), and appearing in the consistency-condition denominator of the local Newton solve) is replaced by \(\tfrac32 H_i\) in both the cone and the mixed-surface return mappings.
Cap mechanism: the cap hardening modulus of the base model is scaled by the same factor,

\[ H_c = 2 H_{pp}\left(\frac{\sigma_3^\ast+c\cot\varphi}{p_{\text{ref}}+c\cot\varphi}\right)^m p^\prime\, H_i. \]

The yield-function shapes and the flow-rule directions (\(\mathbf{m}\), \(\mathbf{n}\) for the cone; \(\mathbf{w}\), \(\mathbf{u}\) for the cap) are unaffected by \(H_i\) — only the rate at which the respective yield surface expands is modified. A stiffer, more slowly hardening small-strain response and a correspondingly larger apparent hardening rate combine so that, for a given accumulated plastic strain, less plastic shear strain (and hence less stress relaxation towards the base-model backbone) accrues within the small-strain range; the two monotonic curves converge again as strain grows and \(G_m \to 1\).

Degeneration to the Base Model

Setting \(G_0^{\text{ref}} = G_{ur}^{\text{ref}} = E_{ur}^{\text{ref}}/2(1+\nu_{ur})\) makes \(\text{Rng}=0\), so all string lengths and stiffness proportions vanish, \(G_{\text{ref},t} \equiv G_{ur}^{\text{ref}}\) for any strain path, \(G_m \equiv 1\), and \(H_i \equiv 1\) identically. In this limit the model is, to machine precision, the base Hardening Soil (MN) model — this equivalence is checked as part of the model's element-test regression suite.

Integration Algorithm Notes

Two independent sub-incrementation mechanisms coexist in this model:

The brick sub-incrementation described above, triggered whenever the shear strain invariant of the full increment (or of a sub-increment of item 2) exceeds \(\gamma_{\text{sub}}^{\max}=5\times10^{-5}\). This refines only the elastic predictor and the resulting \(E_t\)/\(H_i\); it is active regardless of whether the increment is elastic or plastic.
The base model's friction-based sub-stepping, triggered when the trial mobilised friction angle exceeds \(\varphi\) (see base model), which re-integrates the entire increment (elastic predictor and return mapping together) in smaller strain sub-increments. Within this outer loop, the barotropy factor and the brick state are both re-evaluated at the start of every outer sub-increment, and the brick predictor (item 1) is called once per outer sub-increment — since the outer sub-increment size is typically well below \(\gamma_{\text{sub}}^{\max}\), this call almost always resolves to a single inner brick step.

The BRICK state (man strain, brick anchor strains, and \(G_m\)) is driven by the total strain increment at all times, including during plastic sub-increments — this is consistent with the HS-Small/HS-Brick model definition, where the small-strain stiffness history tracks total (not elastic) strain.

When the outer friction-based sub-stepping is triggered, the whole-increment brick trial performed while deciding whether to sub-step is discarded and the brick state is reset to its value at the start of the increment before the sub-stepping loop begins, so that the bricks see the strain increment exactly once, split into the outer sub-increments.

As in the base model, the return-mapping subroutines do not modify the elastic operator passed to them (intent(in)); the tangent returned to the caller is the elastic operator built from \(E_t\) at the last evaluated (sub-)increment, not a consistent algorithmic tangent.

Internal Optimisation

The optimisation routine optimize_hs_bricks_internal_constants(...) computes \((\alpha, H_{pp})\) using the identical virtual-oedometer procedure as the base model (see base model), but the return-mapping calls inside the objective function are evaluated with \(H_i = 1\) fixed, i.e. the optimisation reproduces the base HS model's oedometric response. Consequently \(\alpha\) and \(H_{pp}\) retain their base-model meaning, and \(E_{oed}^{ref}\)/\(K_0^{nc}\) are matched at the fully degraded (large-strain) stiffness state — independently of \(\gamma_{0.7}\) and \(G_0^{\text{ref}}\).

Material Parameters

The parameters must be provided in the props array in the following order. Parameters 1–14 are identical to the base model; two additional parameters are appended:

Index	Symbol	Unit	Description
1	\(E_{50}^{ref}\)	F/A	Secant stiffness in standard drained triaxial test.
2	\(E_{oed}^{ref}\)	F/A	Tangent stiffness for primary oedometer loading (target for optimisation).
3	\(E_{ur}^{ref}\)	F/A	Unloading/reloading stiffness.
4	\(m\)	-	Power law exponent for stress-dependent stiffness.
5	\(c\)	F/A	Cohesion.
6	\(\varphi\)	deg	Friction angle (internally converted to radians).
7	\(\psi\)	deg	Dilatancy angle (internally converted to radians).
8	\(\nu_{ur}\)	-	Poisson's ratio for unloading/reloading.
9	\(p_{\text{ref}}\)	F/A	Reference mean stress level.
10	\(K_0^{nc}\)	-	At-rest earth pressure coefficient (normally consolidated).
11	\(R_f\)	-	Failure ratio.
12	\(E_i^{ref}\)	F/A	Initial tangent stiffness.
13	\(\alpha\)	-	Cap aspect ratio parameter.
14	\(H_{pp}\)	-	Cap hardening parameter.
15	\(\gamma_{0.7}\)	-	Shear strain at which \(G_{\text{sec}}/G_0=0.722\) (Hardin–Drnevich threshold).
16	\(G_0^{\text{ref}}\)	F/A	Small-strain reference shear modulus at \(p_{\text{ref}}\).

Automatic determination of cap parameters

As in the base model, if \(\alpha = 0\) and/or \(H_{pp} = 0\), numgeo determines the missing parameter(s) from \(E_{oed}^{ref}\) and \(K_0^{nc}\) (see Internal Optimisation).

Parameter validity

At least 16 entries are required in props. \(G_0^{\text{ref}}\) must not be smaller than \(G_{ur}^{\text{ref}}=E_{ur}^{\text{ref}}/2(1+\nu_{ur})\); if \(G_0^{\text{ref}} > G_{ur}^{\text{ref}}\), \(\gamma_{0.7}\) must be strictly positive. Violations of either condition terminate the analysis with a diagnostic error stop.

State Variables

At least 73 state variables must be provided (nstatev >= 73); this is checked and enforced with error stop. State variables 1–5 are identical to the base model; the BRICK state is stored in variables 6–73:

Index	Symbol	Description
1	`void_ratio`	Void ratio (only meaningful if initialised at analysis start).
2	`gammapss`	Accumulated plastic shear strain measure \(\gamma^p\).
3	`pp`	Preconsolidation stress \(p_p\).
4	`p`	Current mean stress \(p\).
5	`q`	Current deviatoric stress \(q\).
6	`Gm`	Stiffness ratio \(G_m\) (monotonically non-increasing, \(\ge 1\)).
7	`n_bricks`	Number of currently active (dragged) bricks, diagnostic only (\(0\)–\(10\)).
8–13	`sn`	Man strain, tensorial shear components (order \(11,22,33,12,13,23\)).
14–73	`snb`	Brick anchor strains, six tensorial shear components per brick, bricks \(1,\dots,10\) in order (same component order as `sn`).

Initialisation. On the first call of an analysis (istep == 1 .and. iinc == 1), the BRICK state is initialised to the virgin state: all strings slack, man and brick strains at zero, and \(G_m = G_0^{\text{ref}}/G_{ur}^{\text{ref}}\) (no degradation recorded yet). A state read back with \(G_m < 1\) (which cannot occur by construction once initialised) is treated as uninitialised and reset the same way, as a safeguard for restart or geostatic-step conventions where statev may not have been written by this model yet.

Base model

Hardening Soil (MN)
Reference manual

How to use it

References

Marcin Cudny and Andrzej Truty. Refinement of the Hardening Soil model within the small strain range. Acta Geotechnica, 15(8):2031–2051, August 2020. doi:10.1007/s11440-020-00945-5. ↩↩↩
Thomas Benz. Small-Strain Stiffness of Soils and Its Numerical Consequences. PhD thesis, Inst. für Geotechnik, 2007. URL: https://www.igs.uni-stuttgart.de/dokumente/Mitteilungen/55_Benz.pdf. ↩↩