Modified Drucker-Prager/Cap Model

Yield Surface

The addition of the cap yield surface to the Drucker-Prager model serves two main purposes: it bounds the yield surface in hydrostatic compression, thus providing an inelastic hardening mechanism to represent plastic compaction; and it helps to control volume dilatancy when the material yields in shear by providing softening as a function of the inelastic volume increase created as the material yields on the Drucker-Prager shear failure surface.

The yield surface has two principal segments: a pressure-dependent Drucker-Prager shear failure segment and a compression cap segment, as shown in Figure 1. The Drucker-Prager failure segment is a perfectly plastic yield surface (no hardening). Plastic flow on this segment produces inelastic volume increase (dilation) that causes the cap to soften. On the cap surface plastic flow causes the material to compact. The model is described in detail in Drucker-Prager/Cap model for geological materials.

Failure Surface

The Drucker-Prager failure surface is written as

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

where $β (θ, f_{i})$ and $d (θ, f_{i})$ represent the angle of friction of the material and its cohesion, respectively, and can depend on temperature, $θ$ , and other predefined fields $f_{i}, i = 1, 2, 3 \dots$ . The deviatoric stress measure t is defined as

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

and

$p = - \frac{1}{3} trace (σ)$: is the equivalent pressure stress,
$q = \sqrt{\frac{3}{2} S : S}$: is the Mises equivalent stress,
$r = {(\frac{9}{2} S : S \cdot S)}^{\frac{1}{3}}$: is the third stress invariant, and
$S = σ + p I$: is the deviatoric stress.

$K (θ, f_{i})$ is a material parameter that controls the dependence of the yield surface on the value of the intermediate principal stress, as shown in Figure 2.

The yield surface is defined so that K is the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression. $K = 1$ implies that the yield surface is the von Mises circle in the deviatoric principal stress plane (the $Π$ -plane), so that the yield stresses in triaxial tension and compression are the same; this is the default behavior in Abaqus/Standard and the only behavior available in Abaqus/Explicit. To ensure that the yield surface remains convex requires $0.778 \leq K \leq 1.0$ .

Cap Yield Surface

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 inelastic strain: volumetric plastic and/or creep compaction (when yielding on the cap and/or creeping according to the consolidation mechanism, as described later in this section) causes hardening, while volumetric plastic and/or creep dilation (when yielding on the shear failure surface and/or creeping according to the cohesion mechanism, as described later in this section) causes softening. The cap yield surface is

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

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

p_{b} = p_{b} ({ε_{vol}^{i n} |}_{0} + ε_{vol}^{p l} + ε_{vol}^{c r}) .

The volumetric inelastic strain axis in Figure 3 has an arbitrary origin: ${ε_{vol}^{i n} |}_{0} ( = {ε_{vol}^{p l} |}_{0} + {ε_{vol}^{c r} |}_{0})$ is the position on this axis corresponding to the initial state of the material when the analysis begins, thus defining the position of the cap ( $p_{b}$ ) in Figure 1 at the start of the analysis. The evolution parameter $p_{a}$ is given 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 transition yield surface,

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

so that the model provides a smooth intersection between the cap and failure surfaces.

Defining Yield Surface Variables

You provide the variables d, $β$ , R, ${ε_{v o l}^{i n} |}_{0}$ , $α$ , and K to define the shape of the yield surface. In Abaqus/Standard $0.778 \leq K \leq 1.0$ , while in Abaqus/Explicit K= 1 ( $t = q$ ). If desired, combinations of these variables can also be defined as a tabular function of temperature and other predefined field variables.

Input File Usage

CAP PLASTICITY

Defining Hardening Parameters

The hardening curve specified for this model interprets yielding in the hydrostatic pressure sense: the hydrostatic pressure yield stress is defined as a tabular function of the volumetric inelastic strain, and, if desired, a function of temperature and other predefined field variables. The range of values for which $p_{b}$ is defined should be sufficient to include all values of effective pressure stress that the material is subjected to during the analysis. If the value of the volumetric inelastic strain becomes greater than the last user-specified value, Abaqus/Explicit extrapolates the stress-strain relationship based on the slope at the last point of the curve.

Input File Usage

CAP HARDENING

Plastic Flow

Plastic flow is defined by a flow potential that is associated in the deviatoric plane, associated in the cap region in the meridional plane, and nonassociated in the failure surface and transition regions in the meridional plane. The flow potential surface that we use 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 form a continuous and smooth potential surface.

Nonassociated Flow

Nonassociated flow implies that the material stiffness matrix is not symmetric and the unsymmetric matrix storage and solution scheme should be used in Abaqus/Standard (see Defining an Analysis). 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 gives an acceptable rate of convergence; in such cases the unsymmetric matrix scheme may not be needed.

Calibration

At least three experiments are required to calibrate the simplest version of the Cap model: a hydrostatic compression test (an odometer test is also acceptable) and either two triaxial compression tests or one triaxial compression test and one uniaxial compression test (more than two tests are recommended for a more accurate calibration).

The hydrostatic compression test is performed by pressurizing the sample equally in all directions. The applied pressure and the volume change are recorded.

The uniaxial compression test involves compressing the sample between two rigid platens. The load and displacement in the direction of loading are recorded. The lateral displacements should also be recorded so that the correct volume changes can be calibrated.

Triaxial compression experiments are performed using a standard triaxial machine where a fixed confining pressure is maintained while the differential stress is applied. Several tests covering the range of confining pressures of interest are usually performed. Again, the stress and strain in the direction of loading are recorded, together with the lateral strain so that the correct volume changes can be calibrated.

Unloading measurements in these tests are useful to calibrate the elasticity, particularly in cases where the initial elastic region is not well defined.

The hydrostatic compression test stress-strain curve gives the evolution of the hydrostatic compression yield stress, $p_{b} (ε_{v o l}^{p l})$ , required for the cap hardening curve definition.

The friction angle, $β$ , and cohesion, d, which define the shear failure dependence on hydrostatic pressure, are calculated by plotting the failure stresses of the two triaxial compression tests (or the triaxial compression test and the uniaxial compression test) in the pressure stress (p) versus shear stress (q) space: the slope of the straight line passing through the two points gives the angle $β$ and the intersection with the q-axis gives d. For more details on the calibration of $β$ and d, see the discussion on calibration in Extended Drucker-Prager Models.

R represents the curvature of the cap part of the yield surface and can be calibrated from a number of triaxial tests at high confining pressures (in the cap region). R must be between 0.0001 and 1000.0.

Abaqus/Standard Creep Model

Classical “creep” behavior of materials that exhibit plasticity according to the capped Drucker-Prager plasticity model can be defined in Abaqus/Standard. The creep behavior in such materials is intimately tied to the plasticity behavior (through the definitions of creep flow potentials and definitions of test data), so cap plasticity and cap hardening must be included in the material definition. If no rate-independent plastic behavior is desired in the model, large values for the cohesion, d, as well as large values for the compression yield stress, $p_{b}$ , should be provided in the plasticity definition: as a result the material follows the capped Drucker-Prager model while it creeps, without ever yielding. This capability is limited to cases in which there is no third stress invariant dependence of the yield surface ( $K = 1$ ) and cases in which the yield surface has no transition region ( $α = 0$ ). The elastic behavior must be defined using linear isotropic elasticity (see Defining Isotropic Elasticity).

Creep behavior defined for the modified Drucker-Prager/Cap model is active only during soils consolidation, coupled temperature-displacement, and transient quasi-static procedures.

Creep Formulation

This model has two possible creep mechanisms that are active in different loading regions: one is a cohesion mechanism, which follows the type of plasticity active in the shear-failure plasticity region, and the other is a consolidation mechanism, which follows the type of plasticity active in the cap plasticity region. Figure 5 shows the regions of applicability of the creep mechanisms in p–q space.

Equivalent Creep Surface and Equivalent Creep Stress for the Cohesion Creep Mechanism

Consider the cohesion creep mechanism first. We adopt the notion of the existence of creep isosurfaces of stress points that share the same creep “intensity,” as measured by an equivalent creep stress. Since it is desirable to have the equivalent creep surface coincide with the yield surface, we define the equivalent creep surfaces by homogeneously scaling down the yield surface. In the p–q plane the equivalent creep surfaces translate into surfaces that are parallel to the yield surface, as depicted in Figure 6.

Abaqus/Standard requires that cohesion creep properties be measured in a uniaxial compression test. The equivalent creep stress, ${\bar{σ}}^{c r}$ , is determined as follows:

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

Abaqus/Standard also requires that ${\bar{σ}}^{c r}$ be positive. Figure 6 shows such an equivalent creep stress. A consequence of these concepts is that there is a cone in p–q space inside which creep is not active. Any point inside this cone would have a negative equivalent creep stress.

Equivalent Creep Surface and Equivalent Creep Stress for the Consolidation Creep Mechanism

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 hydrostatic pressure surfaces (vertical lines in the p–q plane). Abaqus/Standard 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 of ${\bar{p}}^{c r} = p - p_{a}$ . 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/Standard account for the fact that a positive pressure will produce negative (that is, compressive) volumetric creep components.

Creep Flow

The creep strain rate produced by the cohesion mechanism is assumed to follow a potential that is similar to that of the creep strain rate in the Drucker-Prager creep model (Extended Drucker-Prager Models); that is, a hyperbolic function:

G_{s}^{c r} = \sqrt{{(0.1 \frac{d}{(1 - \frac{1}{3} \tan β)} \tan β)}^{2} + q^{2}} - p \tan β .

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. The cohesion creep potential is the von Mises circle in the deviatoric stress plane (the $Π$ -plane).

Abaqus/Standard protects for numerical problems that may arise for very low stress values. See Drucker-Prager/Cap model for geological materials for details.

The creep strain rate produced by the consolidation mechanism is assumed to follow a potential that is similar to that of the plastic strain rate in the cap yield surface (Figure 8):

G_{c}^{c r} = \sqrt{{[p - p_{a}]}^{2} + {[R q]}^{2}} .

The consolidation creep potential is the von Mises circle in the deviatoric stress plane (the $Π$ -plane). The volumetric components of creep strain from both mechanisms contribute to the hardening/softening of the cap, as described previously. For details on the behavior of these models refer to Verification of creep integration.

Nonassociated Flow

The use of a creep potential for the cohesion mechanism different from the equivalent creep surface implies that the material stiffness matrix is not symmetric, and the unsymmetric matrix storage and solution scheme should be used (see Defining an Analysis). If the region of the model in which cohesive inelastic deformation is occurring is confined, it is possible that a symmetric approximation to the material stiffness matrix will give an acceptable rate of convergence; in such cases the unsymmetric matrix scheme may not be needed.

Specifying Creep Laws

The definition of the creep behavior is completed by specifying the equivalent “uniaxial behavior”—the creep “laws.” In many practical cases the creep laws are defined through user subroutine CREEP because creep laws are usually of complex form to fit experimental data. Data input methods are provided for some simple cases. To avoid drawbacks of the time hardening and strain hardening forms, it is recommended that you use the time power law model rather than the time hardening form and the power law model rather than the strain hardening form, as discussed below.

User Subroutine CREEP

User subroutine CREEP provides a general capability for implementing viscoplastic models in which the strain rate potential can be written as a function of the equivalent stress and any number of “solution-dependent state variables.” When used in conjunction with these materials, the equivalent cohesion creep stress, ${\bar{σ}}^{c r}$ , and the effective creep pressure, ${\bar{p}}^{c r}$ , are made available in the routine. Solution-dependent state variables are any variables that are used in conjunction with the constitutive definition and whose values evolve with the solution. Examples are hardening variables associated with the model. When a more general form is required for the stress potential, user subroutine UMAT can be used.

Input File Usage

Use either or both of the following options:

CAP CREEP, MECHANISM=COHESION, LAW=USER
CAP CREEP, MECHANISM=CONSOLIDATION, LAW=USER

Time Hardening Form

With respect to the cohesion mechanism, the "time hardening" form is available

{\dot{\bar{ε}}}^{c r} = A {({\bar{σ}}^{c r})}^{n} t^{m},

where

${\dot{\bar{ε}}}^{c r}$: is the equivalent creep strain rate;
${\bar{σ}}^{c r}$: is the equivalent cohesion creep stress;
t: is the total or the creep time; and
A, n, and m: are user-defined creep material parameters specified as functions of temperature and field variables.

In using this form with the consolidation mechanism, ${\bar{σ}}^{c r}$ can be replaced by ${\bar{p}}^{c r}$ , the effective creep pressure, in the above relation.

Input File Usage

Use either or both of the following options:

CAP CREEP, MECHANISM=COHESION, LAW=TIME
CAP CREEP, MECHANISM=CONSOLIDATION, LAW=TIME

Strain Hardening Form

For the cohesion mechanism the "strain hardening" form is

{\dot{\bar{ε}}}^{c r} = {(A {({\bar{σ}}^{c r})}^{n} {[(m + 1) {\bar{ε}}^{c r}]}^{m})}^{\frac{1}{m + 1}} .

In using this form with the consolidation mechanism, ${\bar{σ}}^{c r}$ can be replaced by ${\bar{p}}^{c r}$ , the effective creep pressure, in the above relation.

For physically reasonable behavior A and n must be positive and $- 1 < m \leq 0$ .

Input File Usage

Use either or both of the following options:

CAP CREEP, MECHANISM=COHESION, LAW=STRAIN
CAP CREEP, MECHANISM=CONSOLIDATION, LAW=STRAIN

Time Power Law Model

The time power law model has the following form:

{\dot{\bar{ε}}}^{c r} = {\dot{ε}}_{0} {{(\frac{{\bar{σ}}^{c r}}{q_{0}})}^{n} ({\dot{ε}}_{0} t)}^{m},

where ${\dot{\bar{ε}}}^{c r}$ and ${\bar{σ}}^{c r}$ are defined above; and ${\dot{ε}}_{0}$ , $n$ , $m$ , and $q_{0}$ are material parameters.

The model is equivalent to the time hardening form. It is recommended that you use the time power law model when the value of the parameter $A$ is very small ( $A \leq 10^{- 27}$ ). In this case the equivalent time power law model is obtained by setting $q_{0} = {\dot{ε}}_{0}^{\frac{m + 1}{n}} A^{- \frac{1}{n}}$ , keeping the parameters $n$ and $m$ unchanged, and setting ${\dot{ε}}_{0}$ to an arbitrary value greater than zero (typically, ${\dot{ε}}_{0}$ is set to one).

Input File Usage

Use either or both of the following options:

CAP CREEP, MECHANISM=COHESION, LAW=TIME POWER LAW
CAP CREEP, MECHANISM=CONSOLIDATION, LAW=TIME POWER LAW

Power Law Model

The power law model has the following form:

{\dot{\bar{ε}}}^{c r} = {{\dot{ε}}_{0} [{{(\frac{{\bar{σ}}^{c r}}{q_{0}})}^{n} [(m + 1) {\bar{ε}}^{c r}]}^{m})]}^{\frac{1}{m + 1}},

where ${\dot{\bar{ε}}}^{c r}$ , ${\bar{ε}}^{c r}$ and ${\bar{σ}}^{c r}$ are defined above; and ${\dot{ε}}_{0}$ , $n$ , $m$ and $q_{0}$ are material parameters.

This model is equivalent to the strain hardening form. It is recommended that you use the power law model when the value of the parameter $A$ is very small ( $A \leq 10^{- 27}$ ). In this case the equivalent power law model is obtained by setting $q_{0} = {\dot{ε}}_{0}^{\frac{m + 1}{n}} A^{- \frac{1}{n}}$ , keeping the parameters $n$ and $m$ unchanged, and setting ${\dot{ε}}_{0}$ to an arbitrary value greater than zero (typically, ${\dot{ε}}_{0}$ is set to one).

Input File Usage

Define one or both of the following:

CAP CREEP, MECHANISM=COHESION, LAW=POWER LAW
CAP CREEP, MECHANISM=CONSOLIDATION, LAW=POWER LAW

Singh-Mitchell Law

A second cohesion creep law available as data input is a variation of the Singh-Mitchell law:

{\dot{\bar{ε}}}^{c r} = A e^{(α {\bar{σ}}^{c r})} {(t_{1} / t)}^{m},

where ${\dot{\bar{ε}}}^{c r}$ , t, and ${\bar{σ}}^{c r}$ are defined above and A, $α$ , $t_{1}$ , and m are user-defined creep material parameters specified as functions of temperature and field variables. For physically reasonable behavior A and $α$ must be positive, $0.0 < m \leq 1.0$ , and $t_{1}$ should be small compared to the total time.

In using this variation of the Singh-Mitchell law with the consolidation mechanism, ${\bar{σ}}^{c r}$ can be replaced by ${\bar{p}}^{c r}$ , the effective creep pressure, in the above relation.

Input File Usage

Use either or both of the following options:

CAP CREEP, MECHANISM=COHESION, LAW=SINGHM
CAP CREEP, MECHANISM=CONSOLIDATION, LAW=SINGHM

Time-Dependent Behavior

In the time hardening form, the time power law model, and the Singh-Mitchell law model, the total time or the creep time can be used. The total time is the accumulated time over all general analysis steps. The creep time is the sum of the times of the procedures with time-dependent material behavior. If the total time is used, it is recommended that small step times compared to the creep time be used for any steps for which creep is not active in an analysis; this is necessary to avoid changes in hardening behavior in subsequent steps.

Input File Usage

Use one of the following options:

CAP CREEP, TIME=TOTAL (default)
CAP CREEP, TIME=CREEP

Numerical Difficulties

Depending on the choice of units for the creep laws described above, the value of A may be very small for typical creep strain rates. If A is less than 10⁻²⁷, numerical difficulties can cause errors in the material calculations. Therefore, use another system of units or use the time power law or power law model to avoid such difficulties in the calculation of creep strain increments.

Creep Integration

Abaqus/Standard provides both explicit and implicit time integration of creep and swelling behavior. The choice of the time integration scheme depends on the procedure type, the parameters specified for the procedure, the presence of plasticity, and whether or not a geometric linear or nonlinear analysis is requested, as discussed in Rate-Dependent Plasticity: Creep and Swelling.

Initial Conditions

The initial stress at a point can be defined (see Defining Initial Stresses). If such a stress point lies outside the initially defined cap or transition yield surfaces and under the projection of the shear failure surface in the p–t plane (illustrated in Figure 1), Abaqus will try to adjust the initial position of the cap to make the stress point lie on the yield surface and a warning message will be issued. If the stress point lies outside the Drucker-Prager failure surface (or above its projection), an error message will be issued and execution will be terminated.

Elements

The modified Drucker-Prager/Cap material behavior can be used with plane strain, generalized plane strain, axisymmetric, and three-dimensional solid (continuum) elements. This model cannot be used with elements for which the assumed stress state is plane stress (plane stress, shell, and membrane elements).

Output

In addition to the standard output identifiers available in Abaqus (Abaqus/Standard Output Variable Identifiers and Abaqus/Explicit Output Variable Identifiers), the following variables have special meaning in the cap plasticity/creep model:

PEEQ: $p_{b}$ , the cap position.
PEQC: Equivalent plastic strains for all three possible yield/failure surfaces (Drucker-Prager failure surface - PEQC1, cap surface - PEQC2, and transition surface - PEQC3) and the total volumetric inelastic strain (PEQC4). For each yield/failure surface, the equivalent plastic strain is ${\bar{ε}}^{p l} = \int_{0}^{t} \sqrt{\frac{2}{3} {\dot{ε}}^{p l} : {\dot{ε}}^{p l}} d t,$ where ${\dot{ε}}^{p l}$ is the corresponding rate of plastic flow. The total volumetric inelastic strain is defined as $ε_{vol}^{i n} = \int_{0}^{t} {\dot{ε}}_{k k}^{p l} d t + \int_{0}^{t} {\dot{ε}}_{k k}^{c r} d t .$
CEEQ: Equivalent creep strain produced by the cohesion creep mechanism, defined as $\int \frac{σ : d ε^{c r}}{{\bar{σ}}^{c r}},$ where ${\bar{σ}}^{c r} = \frac{(q - p \tan β)}{(1 - \frac{1}{3} \tan β)}$ is the equivalent creep stress.
CESW: Equivalent creep strain produced by the consolidation creep mechanism, defined as $\int \frac{σ : d ε^{c r}}{\bar{p}}$ , where $\bar{p} = \frac{R^{2} q^{2} + p (p - p_{a})}{G_{c}^{c r}}$ is the equivalent creep pressure.