WIPP-Drucker Model

Viscoplasticity is also modeled in FLAC3D or 3DEC by combining the viscoelastic WIPP model with the elasto-plastic Drucker-Prager model. Of the plasticity models currently embodied in FLAC3D/3DEC, the Drucker-Prager model is the most compatible with the WIPP model, because both models are formulated in terms of the second invariant of the deviatoric stress tensor. Viewed in the \(\pi\)-plane, both models exhibit responses that depend only on radial distance from the isotropic-stress locus. The response of the Mohr-Coulomb model, on the other hand, is not isotropic, because the intermediate principal stress does not enter into its formulation.

The following development is slightly different from that presented in the Drucker-Prager model, and is provided to demonstrate the compatibility of the Drucker-Prager formulation with the creep formulation given in the WIPP-reference creep model.

The shear yield function for the Drucker-Prager model is

(1)\[f^s=\tau+q_{\phi}\sigma_\circ-k_\phi\]

where \(f^s\) = 0 at yield, \(\sigma_\circ=\sigma_{kk}/3\), and \(\tau=\sqrt{J_2}\), where \(J_2\) is the second invariant of the deviatoric stress tensor. Parameters \(q_\phi\) and \(k_\phi\) are material properties.

Because

(2)\[J_2=\sigma_{ij}^d\sigma_{ij}^d/2\]

(e.g., see Malvern 1969, p. 93) \(\tau\) may be related to the stress magnitude, \(\bar\sigma\), defined by \(\bar\sigma=\sqrt{3\sigma^d_{ij}\sigma^d_{ij}/2}\):

(3)\[\bar\sigma=\sqrt{3}\,\tau\]

The plastic potential function in shear, \(g^s\), is similar to the yield function, with the substitution of \(q_\psi\) for \(q_\phi\) as a material property that controls dilation (similar to the Drucker-Prager model).

(4)\[g^s=\tau+q_\psi\sigma_\circ\]

If the yield condition (\(f^s\) = 0) is met, the following flow rules apply:

(5)\[\dot\epsilon_{ij}^{dp}=\lambda{\partial g^s\over\partial\sigma_{ij}^d}\]

(6)\[\dot\epsilon_\circ^p=\lambda{\partial g^s\over\partial\sigma_\circ}\]

where \(\lambda\) is a multiplier (not a material property) to be determined from the requirement that the final stress tensor must satisfy the yield condition. Superscript \(p\) denotes “plastic,” and \(d\) denotes “deviatoric.” By differentiating Equations (2) and (6), we obtain

(7)\[\dot\epsilon_{ij}^{dp}=\lambda{\sigma_{ij}^d\over2\tau}\]

(8)\[\dot\epsilon_\circ^p=\lambda q_\psi\]

In this elastic/plastic formulation, these equations are solved simultaneously with the condition \(f^s\) = 0 and the condition that the sum of elastic and plastic strain rates must equal the applied strain rate.

This implementation of the Drucker-Prager model also contains a tensile yield surface, with a composite decision function used near the intersection of the shear and tensile yield functions. The tensile yield surface is

(9)\[f^t=\sigma_\circ-\sigma^t\]

where \(\sigma^t\) is the tensile yield strength. The associated plastic potential function is

(10)\[g^t=\sigma_\circ\]

Using a similar approach to that for shear yield, the strain rates for tensile yield are

(11)\[\dot\epsilon_{ij}^{dp}=0\]

(12)\[\dot\epsilon_\circ^p=\lambda\]

where \(\lambda\) is determined from the condition that \(f^t\) = 0. Note that the tensile strength cannot be greater than the value of mean stress at which \(f^s\) becomes zero (i.e., \(\sigma^t<k_\phi/q_\phi\)).

When both creep and plastic flow occur, we assume that the associated strain rates act “in series”:

(13)\[\dot\epsilon_{ij}^d =\dot\epsilon_{ij}^{de}+\dot\epsilon_{ij}^{dv} +\dot\epsilon_{ij}^{dp}\]

where the terms represent elastic, viscous, and plastic strain rates, respectively. We first treat the case of shear yield, \(f^s\) > 0. Extending Equation (13), we obtain

(14)\[\dot\epsilon_{ij}^d={\dot\sigma_{ij}^d\over2G}+{\sigma_{ij}^d\over2\bar\sigma}\Bigl\{ 3\dot\epsilon + \sqrt{3}\lambda\Bigr\}\]

In contrast to the creep-only model, the volumetric response of the viscoplastic model is not uncoupled from the deviatoric behavior unless \(q_\psi\) = 0, so

(15)\[\dot\epsilon_{kk} = {3\dot\epsilon_\circ} = {\dot\sigma_{kk}\over3K}+\lambda q_\psi\]

The iteration procedure embodied in the creep solution scheme can be extended to include plastic strain increments. We have

(16)\[\sigma_{ij}^{d'}=\sigma_{ij}^{d^\circ}+2G\Delta t \left[\dot\epsilon_{ij}^d -{\sigma_{ij}^d\over2\bar\sigma} \bigl(3\dot\epsilon+\sqrt{3}\lambda\bigr) \right]\]

And Equation (15) becomes

(17)\[\sigma'_\circ=\sigma_\circ^\circ+(\dot\epsilon_{kk} -\lambda q_\psi)K\Delta t\]

The value of \(\lambda\) can be adjusted in each iteration so that the solution converges to \(f^s\) = 0. Using Newton’s method for roots,

(18)\[\lambda'=\lambda^\circ-{f^s\over({\partial f^s\over\partial\lambda})}\]

Note that \(f^s\) is evaluated with “new” stress components, \(\sigma_{ij}^{d'}\). The derivative in Equation (18) can be evaluated:

(19)\[{\partial f^s \over \partial\lambda} ={\partial f^s \over \partial\sigma_{ij}^{d'}} {\partial\sigma_{ij}^{d'} \over \partial\lambda} +{\partial f^s \over \partial\sigma'_\circ} {\partial\sigma'_\circ \over \partial\lambda}\]

Hence,

(20)\[{\partial f^s \over \partial\lambda} = -G\Delta t -Kq_\phi q_\psi\Delta t\]

assuming that the mean stress components (\(\sigma_{ij}^d\) and \(\bar\sigma\)) are constant.

For tensile yield, \(\sigma_\circ>\sigma^t\). Further, if the shear stress is nonzero, the following function is used to decide if shear or tensile yield is occurring:

\[h=\tau-\tau_p-\alpha_p(\sigma_\circ-\sigma^t)\]

where

\[\tau_p=k_\phi-q_\phi\sigma^t\]

\[\alpha_p=\sqrt{1+{q_\phi}^2}-q_\phi\]

Tensile yield is declared if \(h\) < 0; otherwise, shear yield occurs. In the former case, the last (plastic) term of Equation (13) is zero, and the value of \(\lambda\) is such that Equation (17) reduces to

\[\sigma'_\circ=\sigma^t\]

In order to include softening behavior, an accumulated plastic strain, \(\epsilon^{dp}\), is computed, based on the second invariant of the deviatoric strain-increment tensor, as follows:

(21)\[\epsilon^{dp}:=\epsilon^{dp} +\Delta t\sqrt{\dot\epsilon_{ij}^{dp}\dot\epsilon_{ij}^{dp}/2}\]

There is no built-in support for softening tables, but a FISH function that scans the grid every few steps and recomputes properties based on the current value of \(\epsilon^{dp}\) may be written (see WIPP-Drucker Model: Compression Test Showing Localization for an illustration of this technique).

wipp-drucker model properties

Use the following keywords with the zone property (FLAC3D) or zone property (3DEC) command to set these properties of the WIPP-Drucker model.

wipp-drucker

activation-energy f

activation energy, \(Q\)

bulk f

bulk modulus, \(K\)

cohesion-drucker f

material parameter, \(k_ϕ\)

constant-a f

WIPP model constant, \(A\)

constant-b f

WIPP model constant, \(B\)

constant-d f

WIPP model constant, \(D\)

constant-gas f

gas constant, \(R\)

creep-rate-critical f

critical steady-state creep rate, \(ε̇^*_{ss}\)

dilation-drucker f

material parameter, \(q_k\)

exponent f

WIPP model exponent, \(n\)

friction-drucker f

material parameter, \(q_ϕ\)

temperature f

zone temperature, \(T\)

poisson f

Poisson’s ratio, \(v\)

shear f

shear modulus, \(G\)

tension f

tension limit, \(σ^t\). The default is 0.0.

young f

Young’s modulus, \(E\)

creep-strain-primary f (r)

primary creep strain, \(ε̇_s\)

creep-rate-primary f (r)

primary creep rate, \(ε_s\)

strain-shear-plastic f (r)

accumulated plastic shear strain

strain-tension-plastic f (r)

accumulated plastic tensile strain

Key

(r) Read-only property.

This property cannot be set by the user. Instead, it can be listed, plotted, or accessed through FISH.

Notes:

• Only one of the two options is required to define the elasticity: bulk modulus \(K\) and shear modulus \(G\), or Young’s modulus \(E\) and Poisson’s ratio \(v\).
• The tension cut-off is \(σ^t\) = min (\(σ_t\), \(k_ϕ\)/\(q_ϕ\))
• The creep behavior is triggered by deviatoric stress, while the volumetric behavior does not consider creep.