Drucker-Prager/Cap model for geological materials

Strain rate decomposition

A linear strain rate decomposition is assumed, so

d ε = d ε^{e l} + d ε^{p l} + d ε^{c r},

where $d ε$ is the total strain rate, $d ε^{e l}$ is the elastic strain rate, $d ε^{p l}$ is the inelastic (plastic) time-independent strain rate, and $d ε^{c r}$ is the inelastic (creep) time-dependent strain rate.

Elastic behavior

The elastic behavior can be modeled as linear elastic or by using the porous elasticity model including tensile strength, described in Porous elasticity. If creep has been defined, the elastic behavior must be modeled as linear.

Plastic behavior

The yield/failure surfaces used with this model are written in terms of the three stress invariants: the equivalent pressure stress,

p = - \frac{1}{3} trace (σ);

the Mises equivalent stress,

q = \sqrt{\frac{3}{2} (S : S)};

and the third invariant of deviatoric stress,

r = {(\frac{9}{2} S \cdot S : S)}^{\frac{1}{3}},

where $S$ is the stress deviator, defined as

S = σ + p I .

We also define the deviatoric stress measure

t = \frac{q}{2} [1 + \frac{1}{K} - (1 - \frac{1}{K}) {(\frac{r}{q})}^{3}],

where $K (θ, f^{α})$ is a material parameter that may depend on temperature, $θ$ , and other predefined fields $f^{α}, α = 1, 2, \dots$ . This measure of deviatoric stress is used because it allows matching of different stress values in tension and compression in the deviatoric plane, thereby providing flexibility in fitting experimental results and a smooth approximation to the Mohr-Coulomb surface. Since $r / q = 1$ in uniaxial tension, $t = q / K$ in this case; since $r / q = - 1$ in uniaxial compression, $t = q$ in that case. When $K = 1$ , the dependence on the third deviatoric stress invariant is removed; and the Mises circle is recovered in the deviatoric plane: $t = q$ . Figure 2 shows the dependence of t on K. To ensure convexity of the yield surface, $0.778 \leq K \leq 1.0$ .

With this expression for the deviatoric stress measure, the Drucker-Prager failure surface is written as

F_{s} = t - p \tan β - d = 0,

where $β (θ, f^{α})$ is the material's angle of friction and $d (θ, f^{α})$ is its cohesion (see Figure 1).

The cap yield surface has an elliptical shape with constant eccentricity in the meridional (p–t) plane (Figure 1) and also includes dependence on the third stress invariant in the deviatoric plane (Figure 2). The cap surface hardens or softens as a function of the volumetric plastic strain: volumetric plastic compaction (when yielding on the cap) causes hardening, while volumetric plastic dilation (when yielding on the shear failure surface) causes softening. The cap yield surface is written as

F_{c} = \sqrt{{(p - p_{a})}^{2} + {[\frac{R t}{(1 + α - α / \cos β)}]}^{2}} - R (d + p_{a} \tan β) = 0,

where $R (θ, f^{α})$ is a material parameter that controls the shape of the cap, $α (θ, f^{α})$ is a small number that is defined below, and $p_{a}$ is an evolution parameter that represents the volumetric plastic strain driven hardening/softening. The hardening/softening law is a user-defined piecewise linear function relating the hydrostatic compression yield stress, $p_{b}$ , and the corresponding volumetric inelastic (plastic and/or creep) strain, $p_{b} = p_{b} (ε_{vol}^{i n}, θ, f^{α})$ (Figure 3), where $ε_{vol}^{i n} = ε_{vol}^{p l} + ε_{vol}^{c r}$ .

The evolution parameter, $p_{a}$ , is defined as

p_{a} = \frac{p_{b} - R d}{(1 + R \tan β)} .

The parameter $α$ is a small number (typically 0.01 to 0.05) used to define a smooth transition surface between the shear failure surface and the cap:

F_{t} = \sqrt{{(p - p_{a})}^{2} + {[t - (1 - \frac{α}{\cos β}) (d + p_{a} \tan β)]}^{2}} - α (d + p_{a} \tan β) = 0.

Flow rule

Plastic flow is defined by a flow potential that is associated on the cap and nonassociated on the failure yield surface and transition yield surfaces. The nonassociated nature of these surfaces stems from the shape of the flow potential in the meridional plane. The flow potential surface in the meridional plane is shown in Figure 4.

It is made up of an elliptical portion in the cap region that is identical to the cap yield surface:

G_{c} = \sqrt{{(p - p_{a})}^{2} + {[\frac{R t}{1 + α - α / \cos β}]}^{2}}

and another elliptical portion in the failure and transition regions that provides the nonassociated flow component in the model:

G_{s} = \sqrt{{[(p_{a} - p) \tan β]}^{2} + {[\frac{t}{1 + α - α / \cos β}]}^{2}} .

The two elliptical portions, $G_{c}$ and $G_{s}$ , form a continuous and smooth potential surface.

Nonassociated flow implies that the material stiffness matrix is not symmetric, so the unsymmetric solver should be invoked by the user. However, if the region of the model in which nonassociated inelastic deformation is occurring is confined, it is possible that a symmetric approximation to the material stiffness matrix will give an acceptable convergence rate: in such cases the unsymmetric solver may not be needed.

Creep model

Classical “creep” behavior of materials that also exhibit plastic behavior according to the modified Drucker-Prager/Cap model can be defined.

The creep behavior in such materials is intimately tied to the plasticity behavior (through the definition of creep flow potentials and test data), so it is necessary to define the modified Drucker-Prager/Cap plasticity and hardening behavior as well. The elastic part of the behavior must be linear.

The rate-independent part of the plastic behavior is limited by the following restrictions:

$α = 0$ —that is, no transition zone is allowed;
K=1—that is, no third stress invariant effects are taken into account.

In such a case, the deviatoric stress measure t is equal to the Mises equivalent stress, q, and the yield surface has a von Mises (circular) section in the deviatoric stress plane.

Creep behavior

The built-in Abaqus creep laws or uniaxial laws defined through user subroutine CREEP can be used. The integration of the creep strain rate is first attempted explicitly, as described in Rate-dependent metal plasticity (creep). The integration is done by the backward Euler method (as described in Rate-dependent metal plasticity (creep)) if the stability limit is exceeded, a geometrically nonlinear analysis is being performed, or plasticity becomes active.

In this model we assume the existence of two separate and independent creep mechanisms. One is a cohesion mechanism, which operates similarly to the Drucker-Prager creep model described in Models for granular or polymer behavior. The other is a consolidation mechanism, which operates similarly to the cap zone plasticity. We then have

d ε^{c r} = d ε_{s}^{c r} + d ε_{c}^{c r},

where $d ε_{s}^{c r}$ is the creep strain rate due to the cohesion mechanism and $d ε_{c}^{c r}$ is the creep strain rate due to the consolidation mechanism.

As described above, the cap surface hardens or softens as a function of the volumetric plastic strain and volumetric creep strain: volumetric inelastic compaction (when yielding on the cap or creeping through the consolidation mechanism) causes hardening, while volumetric plastic dilation (when yielding on the shear failure surface or creeping through the cohesion mechanism) causes softening. The separation between the two yield surfaces and the dominant regions for the two creep mechanisms are defined by the evolution parameter, $p_{a}$ , which relates to the user-defined hydrostatic compression yield stress, $p_{b} = p_{b} (ε_{vol}^{p l} + ε_{vol}^{c r}, θ, f^{α})$ (Figure 3).

The cohesion mechanism is active for all stress states that have a positive equivalent creep stress as explained below. The consolidation mechanism is active for all stress states in which the pressure is larger than $p_{a}$ . Figure 5 illustrates the active regions in this formulation.

We adopt the notion of the existence of creep isosurfaces (or equivalent creep surfaces) of stress points that share the same creep “intensity,” as measured by an equivalent creep stress. Consider the cohesion creep mechanism first. When the material plastifies, the equivalent creep surface should coincide with the yield surface; therefore, we define the equivalent creep surfaces by homogeneously scaling down the yield surface. In the p–q plane that translates into parallels to the yield surface, as depicted in Figure 6. Abaqus requires that cohesion creep properties be measured in a uniaxial compression test.

The equivalent creep stress, ${\bar{σ}}^{c r}$ , is determined as the intersection of the equivalent creep surface with the uniaxial compression curve. As a result,

\begin{matrix} (1) & {\bar{σ}}^{c r} = \frac{(q - p \tan β)}{(1 - \frac{1}{3} \tan β)}, \end{matrix}

where $β (θ, f^{α})$ is the material angle of friction. Figure 6 shows how the equivalent creep stress was determined. In uniaxial compression $p = \frac{1}{3} {\bar{σ}}^{c r}$ ; therefore, the uniaxial compression test line has a slope of 1/3. This approach has several consequences. One is that the cohesion creep strain rate is a function of both q and p. This allows the determination of realistic material properties in cases in which, due to high hydrostatic pressures, q is very high. If one looks at the yield strength of the material in this region to be a composite of cohesion strength and friction strength, this model corresponds to cohesion-determined creep. As a result, there is a cone in p–q space inside which there is no cohesion creep.

Next consider the consolidation creep mechanism. In this case we wish to make creep dependent on the hydrostatic pressure above a threshold value of $p_{a}$ , with a smooth transition to the areas in which the mechanism is not active ( $p < p_{a}$ ). Therefore, we define equivalent creep surfaces as constant pressure surfaces. In the p–q plane that translates into vertical lines. Abaqus requires that consolidation creep properties be measured in a hydrostatic compression test. The effective creep pressure, ${\bar{p}}^{c r}$ , is then the point on the p-axis with a relative pressure

\begin{matrix} (2) & {\bar{p}}^{c r} = p - p_{a} . \end{matrix}

This value is used in the uniaxial creep law. The equivalent volumetric creep strain rate produced by this type of law is defined as positive for a positive equivalent pressure. The internal tensor calculations in Abaqus will account for the fact that a positive pressure will produce negative (that is, compressive) volumetric creep components.

Creep flow rule

The creep flow rules are derived from creep potentials, $G^{c r}$ , in such a way that

\begin{matrix} (3) & d ε^{c r} = \frac{d {\bar{ε}}^{c r}}{f^{c r}} \frac{\partial G^{c r}}{\partial σ}, \end{matrix}

where $d {\bar{ε}}^{c r}$ is the equivalent creep strain rate, which must be work conjugate to the equivalent creep stress:

\begin{matrix} d {\bar{ε}}^{c r} = ​ & | d ε_{11}^{c r} | in
  the uniaxial compression
  case, \\ = & | d ε_{v o l}^{c r} | in
  the volumetric compression
  case . \end{matrix}

Since $d ε^{c r}$ is obviously work conjugate to $σ$ , $f^{c r}$ is a proportionality factor defined by

\begin{matrix} (4) & f^{c r} = \frac{1}{{\bar{σ}}^{c r}} σ : \frac{\partial G^{c r}}{\partial σ}, \end{matrix}

with

\begin{matrix} {\bar{σ}}^{c r} = ​ & | σ_{11} | in
  the uniaxial compression
  case, \\ = & p in
  the volumetric compression
  case . \end{matrix}

Cohesion creep

For the cohesion mechanism the creep potential is assumed to follow the same potential as the creep strain rate in the Drucker-Prager creep model (Models for granular or polymer behavior); that is, a hyperbolic function. This creep flow potential, which is continuous and smooth, ensures that the flow direction is always uniquely defined. The function approaches a parallel to the shear-failure yield surface asymptotically at high confining pressure stress and intersects the hydrostatic pressure axis at a right angle. A family of hyperbolic potentials in the meridional stress plane is shown in Figure 7:

\begin{matrix} (5) & G_{s}^{c r} = \sqrt{{(0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \tan β)}^{2} + q^{2}} - p \tan β, \end{matrix}

where d is the material cohesion.

The equivalent cohesion creep strain rate is then determined from the uniaxial law:

Δ {\bar{ε}}_{s}^{c r} = h_{s} ({\bar{σ}}^{c r}, {\bar{ε}}_{s}^{c r}, θ, f^{α}) .

Equation 1, Equation 3, and Equation 5 produce the flow rule for this mechanism

Δ ε_{s}^{c r} = \frac{Δ {\bar{ε}}_{s}^{c r}}{f^{c r}} (\frac{q}{\sqrt{{(0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \tan β)}^{2} + q^{2}}} n + \frac{1}{3} \tan β I),

where

n = \frac{\partial q}{\partial σ} = \frac{3}{2} \frac{S}{q},

and

f^{c r} = \frac{{\bar{σ}}^{c r}}{\sqrt{{(0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \tan β)}^{2} + {({\bar{σ}}^{c r})}^{2}}} - \frac{1}{3} \tan β .

The proportionality factor, $f^{c r}$ , is not a constant in this model. Its expression indicates that it will become negative if

{\bar{σ}}^{c r} < 0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \frac{\frac{1}{3} {(\tan β)}^{2}}{\sqrt{1 - {(\frac{1}{3} \tan β)}^{2}}} .

It turns out that below this stress level, which for typical materials will be very low, the stress vector and the normal to the creep potential are pointing in opposite directions:

σ : \frac{\partial G_{s}^{c r}}{\partial σ} < 0,

which is equivalent to

q - p \tan β < q - \frac{q^{2}}{\sqrt{{(0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \tan β)}^{2} + q^{2}}} .

Thus, there is a small zone just outside the “no creep” cone for which this is the case. Consequently, creep data obtained within this zone should show a creep strain rate in the opposite direction from the applied stress at very low stress levels, which will usually not be the case. To overcome this difficulty, Abaqus will modify the creep data entered such that $f^{c r} \geq 0.1$ . Therefore, you would not expect correspondence between calculated creep strains and measured creep properties in a region defined by

{\bar{σ}}^{c r} < 0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \frac{(0.1 + \frac{1}{3} \tan β) \tan β}{\sqrt{1 - {(0.1 + \frac{1}{3} \tan β)}^{2}}} .

This modification is usually not significant, since typical creep analyses have loads that are applied quickly, followed by long-term creep. Hence, the stress level for most of the analysis will usually be well beyond the modified zone.

An example of “slow” loading in which the approximation is visible is included in Verification of creep integration. As is clear in the example, the effect of the approximation is small in spite of the fact that the load is ramped up over the step.

The equivalent cohesion creep strain rate is a function of both q and p through ${\bar{σ}}^{c r}$ . The creep potential is the von Mises circle in the deviatoric stress plane (the $Π$ -plane). Although creep flow is associated in the deviatoric stress plane, the use of a creep potential different from the equivalent creep surface implies that creep flow is nonassociated.

Consolidation creep

For the consolidation mechanism the creep potential is derived from the plastic potential of the cap zone (Figure 8):

\begin{matrix} (6) & G_{c}^{c r} = \sqrt{{(p - p_{a})}^{2} + {(R q)}^{2}} . \end{matrix}

Recall that this mechanism is active only if $p \geq p_{a}$ .

The equivalent consolidation creep strain rate is then determined from the uniaxial law

Δ {\bar{ε}}_{c}^{c r} = h_{c} ({\bar{p}}^{c r}, {\bar{ε}}_{c}^{c r}, θ, f^{α}) .

Equation 3 and Equation 4 produce $f^{c r} = 1$ ; and Equation 2, Equation 3, and Equation 6 produce the flow rule for this mechanism:

Δ ε_{c}^{c r} = \frac{Δ {\bar{ε}}_{c}^{c r}}{G_{c}^{c r}} (R^{2} q n - \frac{1}{3} (p - p_{a}) I) .

Note that there is an equivalent pressure stress, $\bar{p}$ , work conjugate of the equivalent consolidation creep strain, which is different from the effective creep pressure, ${\bar{p}}^{c r}$ . Such equivalent pressure stress is given by

\bar{p} = \frac{R^{2} q^{2} + p (p - p_{a})}{G_{c}^{c r}}

and has the characteristic that it reduces to the pressure in a hydrostatic compression test.

The creep potential is the von Mises circle in the deviatoric stress plane (the $Π$ -plane). Creep flow is nonassociated in this mechanism.

This formulation is quite simplistic and ignores the effects of q on the creep function, $h_{c}$ . The two creep mechanisms operate independently from each other. This implies that $h_{s}$ does not depend on ${\bar{ε}}_{c}^{c r}$ and that $h_{c}$ does not depend on ${\bar{ε}}_{s}^{c r}$ . The only cross effects between both mechanisms are obtained through the dependency of $p_{a}$ on the volumetric creep from any of them.