Drucker-Prager Elastoplasticity
The Drucker-Prager plasticity model is widely used for simulating granular materials such as sand and soil. It generalizes the von Mises model by incorporating a friction angle, which governs how much shear stress the material can sustain relative to the normal stress. Physically, this corresponds to Coulomb-like friction between particles: the material yields when the shear stress exceeds a friction-dependent bound based on pressure.
In stress space, the Drucker-Prager yield surface takes the shape of a cone, with three distinct cases to handle:
- Case I (Elastic): The stress lies strictly inside the cone, and no plasticity occurs.
- Case II (Expansion): The stress corresponds to a configuration where the material expands volumetrically (positive trace), and no resistance is applied—this maps to the cone tip.
- Case III (Shearing): The stress lies outside the cone but with compressive pressure and must be projected back to the cone surface.
This model is best implemented in the log-strain (Hencky strain) space using the SVD of the deformation gradient, which we have already introduced in the previous section.
Example 31.1.1 (Drucker-Prager Yield Criterion, Log-Strain Formulation). We define the logarithmic strain as:
The deviatoric part is:
The plastic multiplier is computed as:
Here, is the Drucker-Prager friction coefficient derived from the friction angle .
Then we apply the return mapping:
- If (Case II), we project to the cone tip: set .
- If , we are inside the cone (Case I): no change.
- Otherwise (Case III), we project back to the cone surface:
Finally, we compute the updated singular values:
and reconstruct the elastic deformation:
Example 31.1.2 (Drucker-Prager Plasticity with Volume Correction).
In granular materials like sand, volumetric expansion can result in non-physical volume gain if not properly handled. The standard Drucker-Prager projection maps stress to the tip of the cone under expansion (positive trace), which corresponds to a stress-free state. However, this may unrealistically "lock in" expanded configurations as new rest shapes.
This effect can cause persistent volume inflation when a particle experiences elastic expansion followed by plastic projection. Any future compression then incurs artificial elastic penalties, resulting in incorrect material response.
To correct this, we follow the volume correction treatment described by [Tampubolon et al. 2017] by introducing a per-particle scalar accumulator that tracks log-volume changes induced by plastic projection:
This correction is naturally integrated into the log-strain formulation by adjusting the strain before return mapping:
where is the spatial dimension. This allows future compression to neutralize previous volume gain rather than being resisted elastically. In the code below,
diff_log_J
provides this volume correction term, computed as the accumulation of log-difference of determinants.
Implementation 31.1.1 (Drucker-Prager Elastoplasticity Return Mapping, simulator.py).
@ti.func
def Drucker_Prager_return_mapping(F, diff_log_J):
dim = ti.static(F.n)
sin_phi = ti.sin(friction_angle_in_degrees/ 180.0 * ti.math.pi)
friction_alpha = ti.sqrt(2.0 / 3.0) * 2.0 * sin_phi / (3.0 - sin_phi)
U, sig_diag, V = ti.svd(F)
sig = ti.Vector([ti.max(sig_diag[i,i], 0.05) for i in ti.static(range(dim))])
epsilon = ti.log(sig)
epsilon += diff_log_J / dim # volume correction treatment
trace_epsilon = epsilon.sum()
shifted_trace = trace_epsilon
if shifted_trace >= 0:
for d in ti.static(range(dim)):
epsilon[d] = 0.
else:
epsilon_hat = epsilon - (trace_epsilon / dim)
epsilon_hat_norm = ti.sqrt(epsilon_hat.dot(epsilon_hat)+1e-8)
delta_gamma = epsilon_hat_norm + (dim * lam + 2. * mu) / (2. * mu) * (shifted_trace) * friction_alpha
epsilon -= (ti.max(delta_gamma, 0) / epsilon_hat_norm) * epsilon_hat
sig_out = ti.exp(epsilon)
for d in ti.static(range(dim)):
sig_diag[d,d] = sig_out[d]
return U @ sig_diag @ V.transpose()
The return mapping is enforced during particle state update:
Implementation 31.1.2 (Particle State Update, simulator.py).
@ti.kernel
def update_particle_state():
for p in range(N_particles):
# trial elastic deformation gradient
F_tr = F[p]
# apply return mapping to correct the trial elastic state, projecting the stress induced by F_tr
# back onto the yield surface, following the direction specified by the plastic flow rule.
new_F = Drucker_Prager_return_mapping(F_tr, diff_log_J[p])
# track the volume change incurred by return mapping to correct volume, following https://dl.acm.org/doi/10.1145/3072959.3073651 sec 4.3.4
diff_log_J[p] += -ti.log(new_F.determinant()) + ti.log(F_tr.determinant())
F[p] = new_F
# advection
x[p] += dt * v[p]