DafaliasManzari

The modified Dafalias--Manzari sand model.

The original model 10.1061/(ASCE)0733-9399(2004)130:6(622) is modified slightly. Readers can also refer to the corresponding section in Constitutive Modelling Cookbook for details on the theory.

The modifications can be summaries as follows.

1. The continuum mechanics convention (tension positive) is used. All volumetric strain and hydrostatic stress related quantities shall flip their signs.
2. The Lode angle dependency is removed, which is equivalent to set $c=1$ in Eq. (19).
3. Constants such as $2/3$ and $\sqrt{2/3}$ are removed. They can be combined with the model parameters.
4. The hardening related parameter $h$ as defined in Eq. (24) causes numerical issues under small cyclic loads. Hence, it is changed to a similar form.

The algorithm is correct but the theory needs validation. Your collaboration would be much appreciated.

Syntax

material DafaliasMazanri (1) (2-18) [19]
# (1) int, unique material tag
# (2) double, reference shear modulus G_0 ==> 125
# (3) double, poissons ratio \nu
# (4) double, \alpha^c ==> 1.5
# (5) double, slope of critical state line \lambda_c ==> 0.02
# (6) double, initial void ratio e_0
# (7) double, exponent \xi ==> 0.7
# (8) double, initial yield surface size m ==> 0.01
# (9) double, hardening constant h_0
# (10) double, hardening constant h_1 ==> 0.1
# (11) double, hardening constant c_h ==> 0.9
# (12) double, n^b ==> 1.1
# (13) double, dilatancy constant A_0 ==> -0.7
# (14) double, n^d ==> 3.5
# (15) double, z_{max} ==> 4
# (16) double, c_z ==> 600
# (17) double, atmospheric pressure p_{at} ==> -130
# (18) double, threshold G_r ==> 0.1
# [19] double, density, default: 0.0

Theory

Hyperelasticity

The hyperelastic response is defined as

$$G=G_0{\displaystyle \frac{{(2.97-e)}^2}{1+e}}\sqrt{p{p}_{at}},K={\displaystyle \frac{2}{3}}{\displaystyle \frac{1+\nu }{1-2\nu }}G.$$

To improve numerical stability, $G$ is bounded by $G_r{G}_0|p_{at}|$ where $G_r$ is a constant can be chosen as for example $0.1$. This is equivalent to define an elastic response for $|p|<0.01|p_{at}|$.

The void ratio can be associated to strain so that

$$e=e_0+(1+e_0)\mathrm{tr}{\epsilon }^{tr}.$$

The strain increment can be decomposed into elastic and plastic parts.

$${\epsilon }^{tr}={\epsilon }_n+\mathrm{\Delta }\epsilon ={\epsilon }_n+\mathrm{\Delta }{\epsilon }^e+\mathrm{\Delta }{\epsilon }^p.$$

As such, the stress increment can be expressed accordingly,

$$\sigma ={\sigma }_n+\mathrm{\Delta }\sigma ={\sigma }_n+2G(\mathrm{\Delta }\mathbf{e}-\mathrm{\Delta }{\mathbf{e}}^p)+K(\mathrm{\Delta }{\epsilon }_v-\mathrm{\Delta }{\epsilon }_v^p)\mathbf{I}.$$

In deviatoric and spherical components,

$$\begin{gathered} \sigma =\mathbf{s}+p\mathbf{I}, \\ p=p_n+K(\mathrm{\Delta }{\epsilon }_v-\mathrm{\Delta }{\epsilon }_v^p), \\ \mathbf{s}={\mathbf{s}}_n+2G(\mathrm{\Delta }\mathbf{e}-\mathrm{\Delta }{\mathbf{e}}^p), \end{gathered}$$

with

$$\mathrm{\Delta }\epsilon =\mathrm{\Delta }\mathbf{e}+{\displaystyle \frac{1}{3}}\mathrm{\Delta }{\epsilon }_v\mathbf{I},$$

where $\mathbf{s}=\mathrm{dev}\sigma$ is the deviatoric stress, $p={\displaystyle \frac{1}{3}}\mathrm{tr}\sigma$ is the hydrostatic stress.

Critical State

The critical state parameter is chosen as

$$\psi =e-e_0+{\lambda }_c{\left({\displaystyle \frac{p}{p_{at}}}\right)}^{\xi }.$$

The dilatancy surface is defined as

$${\alpha }^d={\alpha }^c\exp \left(n^d\psi \right).$$

The bounding surface is defined as

$${\alpha }^b={\alpha }^c\exp (-n^b\psi ).$$

Yield Function

A wedge-like function is chosen to be the yield surface.

$$F=|\mathbf{s}+p\alpha |+mp=|\eta |+mp,$$

where $\alpha$ is the so called back stress ratio and $m$ characterises the size of the wedge. For simplicity, $m$ is assumed to be a constant in this model.

By denoting $\eta =\mathbf{s}+p\alpha$, the directional unit tensor is defined as

$$\mathbf{n}={\displaystyle \frac{\eta }{|\eta |}}.$$

Flow Rule

A non-associated plastic flow is used, the corresponding flow rule is defined as follows.

$$\mathrm{\Delta }{\epsilon }^p=\gamma (\mathbf{n}+{\displaystyle \frac{1}{3}}D\mathbf{I}),$$

where $D$ is the dilatancy parameter.

$$D=A_d({\alpha }^d-m-\alpha :\mathbf{n})=A_0(1+⟨\mathbf{z}:\mathbf{n}⟩)({\alpha }_d-m-\alpha :\mathbf{n}).$$

Due to the change of sign convention, a negative $D$ gives contractive response, thus $A_0$ often needs to be negative.

Hardening Rule

The evolution rate of the back stress ratio $\alpha$ is defined in terms of a proper distance measure from the bounding surface,

$$\mathrm{\Delta }\alpha =\gamma h(({\alpha }^b-m)\mathbf{n}-\alpha ),$$

where $h$ controls the hardening rate,

$$h=b_0\exp (h_1({\alpha }_{in}-\alpha ):\mathbf{n}).$$

The parameter $b_0$ is defined as a function of current state,

$$b_0=G_0{h}_0(1-c_{h}e)\sqrt{{\displaystyle \frac{p_{at}}{p}}}.$$

${\alpha }_{in}$ is updated whenever load reversal occurs.

Fabric Effect

The fabric tensor changes when $\mathrm{\Delta }{\epsilon }_v^p$ is positive,

$$\mathrm{\Delta }\mathbf{z}=c_z⟨\mathrm{\Delta }{\epsilon }_v^p⟩(z_m\mathbf{n}-\mathbf{z}).$$