Drucker-Prager model - DruckerPrager

The Drucker-Prager plasticity model¹ is an isotropic elasto-plastic model based on a yield function

$$f \left(\mbox{\boldmath$\sigma$}, \tau_{\mathrm {Y}}\right) = F\left(\mbox{\boldmath$\sigma$}\right) - \tau_{\mathrm {Y}}$$ (19)

with the pressure-dependent equivalent stress

$$F\left(\mbox{\boldmath$\sigma$}\right) = \alpha I_1 + \sqrt{J_2}$$ (20)

As usual, $\sigma$ is the stress tensor, $\tau_{\mathrm {Y}}$ is the yield stress under pure shear, and $I_1$ and $J_2$ are the first invariant and second deviatoric invariant of the stress tensor. The friction coefficient $\alpha$ is a positive parameter that controls the influence of the pressure on the yield limit, important for cohesive-frictional materials such as concrete, soils or other geomaterials. Regarding Mohr-Coulomb plasticity, relation to cohesion, $c$, and the angle of friction, $\theta$, exists for the Drucker-Prager model

$$\alpha = \frac{2\sin\theta}{(3-\sin\theta)\sqrt{3}}$$ (21)

$$\tau_{\mathrm {Y}}= \frac{6c\cos\theta}{(3-\sin\theta)\sqrt{3}}$$ (22)

The flow rule is derived from the plastic potential

$$g\left(\mbox{\boldmath$\sigma$}\right) = \alpha_{\psi}I_1 + \sqrt{J_2}$$ (23)

where $\alpha_{\psi}$ is the dilatancy coefficient. An associated model with $\alpha_{\psi}=\alpha$ would overestimate the dilatancy of concrete, so the dilatancy coefficient is usually chosen smaller than the friction coefficient. The model is described by the equations

$$\mbox{\boldmath$\sigma$} = \mbox{\boldmath$D$} : \left(\mbox{\boldmath$\varepsilon$}- \mbox{\boldmath$\varepsilon$}_{\mathrm{p}}\right)$$ (24)

$$\tau_{\mathrm {Y}} = h(\kappa)$$ (25)

$$\dot{\mbox{\boldmath$\varepsilon$}}_{\mathrm{p}}= \dot{\lambda} \frac{\partial g}{\partial \mbox{\boldmath$\sigma$}} = \dot{\lambda} \left( \alpha_{\psi}\mbox{\boldmath$\delta$} + \frac{\mbox{\boldmath$s$}}{2\sqrt{J_2}} \right)$$ (26)

$$\dot{\kappa} = \sqrt{\frac{2}{3}} \; \Vert \dot{\mbox{\boldmath$\varepsilon$}}_{\mathrm{p}}\Vert$$ (27)

$$\dot{\lambda} \ge 0, \;\;\; f \left(\mbox{\boldmath$\sigma$}, \ta... ...\dot{\lambda}\, f \left(\mbox{\boldmath$\sigma$}, \tau_{\mathrm {Y}}\right) = 0$$ (28)

which represent the linear elastic law, hardening law, evolution laws for plastic strain and hardening variable, and the loading-unloading conditions. In the above, $D$ is the elastic stiffness tensor, $\varepsilon$ is the strain tensor, $\varepsilon$$_{\mathrm{p}}$ is the plastic strain tensor, $\lambda$ is the plastic multiplier, $\delta$ is the unit second-order tensor, $s$ is the deviatoric stress tensor, $\kappa$ is the hardening variable, and a superior dot marks the derivative with respect to time. The flow rule has the form given in Eq. (26) at all points of the conical yield surface with the exception of its vertex, located on the hydrostatic axis.

For the present model, the evolution of the hardening variable can be explicitly linked to the plastic multiplier. Substituting the flow rule (26) into Eq. (27) and computing the norm leads to

$$\dot\kappa = k \dot \lambda$$ (29)

with a constant parameter $k = \sqrt{1/3 + 2 \alpha_{\psi}^2}$, so the hardening variable is proportional to the plastic multiplier. For $\alpha=\alpha_{\psi}=0$, the associated $J_2$-plasticity model is recovered as a special case.

In the simplest case of linear hardening, the hardening function is a linear function of $\kappa$, given by

$$h(\kappa)=\tau_0+H E\kappa$$ (30)

where $\tau_0$ is the initial yield stress, and $H$ is the hardening modulus normalized with the elastic modulus. Alternatively, an exponential hardening function

$$h(\kappa) = \tau_{\mathrm{limit}} + \left( \tau_0 - \tau_{\mathrm{limit}} \right) \rm e^{-\kappa/\kappa_{\mathrm{c}}}$$ (31)

can be used for a more realistic description of hardening.

The stress-return algorithm is based on the Newton-iteration. In plasticity, this is commonly called Closest-Point-Projection (CPP), and it generally leads to quadratic convergence. The implemented algorithm is convergent in any stress case, but in the vicinity of the vertex region, quadratic convergence might be lost because of insufficient regularity of the yield function.

The algorithmic tangent stiffness matrix is implemented for both the regular case and the vertex region. Generally, the error decreases quadratically (of course only asymptotically). Again, in the vicinity of the vertex region, quadratic convergence might be lost due to insufficient regularity. Furthermore, the tangent stiffness matrix does not always exist for the vertex case. In these cases, the elastic stiffness is used instead. It is generally safer (but slower) to use the elastic stiffness if you encounter any convergence problems, especially if your problem is tension-dominated.

Table 8: DP material - summary.

Description	DP material
Record Format	DruckerPrager num(in) # d(rn) # tAlpha(rn) # E(rn) # n(rn) # alpha(rn) # alphaPsi(rn) # ht(in) # iys(rn) # lys(rn) # hm(rn) # kc(rn) # [ yieldtol(rn) #]
Parameters	- num material model number
	- d material density
	- tAlpha thermal dilatation coefficient
	- E Young modulus
	- n Poisson ratio
	- alpha friction coefficient
	- alphaPsi dilatancy coefficient
	- ht hardening type, 1: linear hardening, 2: exponential hardening
	- iys initial yield stress in shear, $\tau_0$
	- lys limit yield stress for exponential hardening, $\tau_{\mathrm{limit}}$
	- hm hardening modulus normalized with E-modulus (!)
	- kc $\kappa_c$ for the exponential softening law
	- yieldtol tolerance of the error in the yield criterion, default value 1.e-14
	- newtonIter maximum number of iterations in $\lambda$ search, default value 30
Supported modes	3dMat, PlaneStrain, 3dRotContinuum