岩土力学：Drucker-Prager 三轴压缩¶
本示例在固体力学示例系列中介绍了一个小型岩土力学应用。它针对一种压力相关的土体或软岩材料，实现了一个仅作示例用的 Drucker-Prager 塑性装配器，并通过一条简单的三轴压缩应变路径进行驱动。
目标有意保持简单：在不增加公共 API 的前提下，演示 TensorMesh 在岩土力学中的使用约定。局部装配器位于 examples/solid/geomechanics/drucker_prager_triaxial/drucker_prager_triaxial.py。
问题¶
模型为小应变、关联 Drucker-Prager 塑性，并带有线性各向同性硬化。TensorMesh 内部保持固体力学约定：应力以受拉为正。为便于岩土力学报告，脚本将轴向应力和平均压力以受压为正的形式输出。
屈服函数在内部写为
\[f(\sigma, \alpha) = q + \eta I_1 - (k + H\alpha) \le 0,\]
其中
\[I_1 = \mathrm{tr}(\sigma), \qquad
q = \sqrt{\frac{3}{2}\,s:s}, \qquad
p = -\frac{I_1}{3}.\]
由于在受拉为正约定下受压会导致 I1 为负，更高的围压会降低 f 并延迟屈服。
历史变量¶
本示例遵循现有 J2 塑性示例的相同模式：
逐求积点的历史变量存储在
self.history[etype]中；上一步的
eps_p和alpha通过element_data传递；update_state(u)在每个收敛的载荷步结束后，于torch.no_grad()上下文中被调用。
本示例仅保留在脚本内部，以便在确立约定的示例经过审阅后，再讨论公共 API。
合理性检验¶
运行两种围压水平：
p0 = 0 kPa；p0 = 100 kPa。
脚本验证更高围压的情况更晚到达弹性试探屈服面，且提交的塑性历史变量是单调的。
运行¶
cd examples/solid/geomechanics/drucker_prager_triaxial
python drucker_prager_triaxial.py
若仅进行快速数值计算而不输出图像：
python drucker_prager_triaxial.py --no-plot --steps 16
默认运行会写入 drucker_prager_triaxial.png，其中包含轴向应力与塑性历史曲线。
图 45 Output of drucker_prager_triaxial.py. The left panel shows the
compression-positive axial stress response; the right panel shows the
committed plastic history variable. The higher-confinement case yields
later, matching the pressure-dependent Drucker-Prager sanity check.¶
核心实现¶
"""Example-only associated Drucker-Prager plasticity assembler.
This class is deliberately kept inside the example. It is not a public
TensorMesh API. The implementation follows the same high-level lifecycle
as J2Plasticity:
1. store per-quadrature history in ``self.history[etype]``;
2. pass previous-step state through ``element_data`` during energy calls;
3. call ``update_state(u)`` after each converged load step.
Notes
-----
TensorMesh uses tension-positive stress. With compression-positive mean
pressure ``p = -tr(sigma) / 3``, the yield function is written internally as
f = q + eta * I1 - (k + H * alpha) <= 0,
where ``I1 = tr(sigma)`` and ``q = sqrt(3/2 s:s)``. Because compression
gives negative ``I1``, confinement increases the yield stress.
"""
def __post_init__(self, params: DruckerPragerParameters | None = None):
if params is None:
params = DruckerPragerParameters()
self.params = params
self.E = float(params.E)
self.nu = float(params.nu)
self.cohesion = float(params.cohesion)
self.friction_angle_deg = float(params.friction_angle_deg)
self.H = float(params.H)
self.mu = self.E / (2.0 * (1.0 + self.nu))
self.bulk = self.E / (3.0 * (1.0 - 2.0 * self.nu))
phi = math.radians(self.friction_angle_deg)
sin_phi = math.sin(phi)
cos_phi = math.cos(phi)
# Triaxial-compression meridian fit to Mohr-Coulomb:
# q = M p + k, p compression-positive.
# With TensorMesh tension-positive stress, p = -I1 / 3, so
# f = q + (M/3) I1 - k.
self.M = 6.0 * sin_phi / (3.0 - sin_phi)
self.eta = self.M / 3.0
self.k = 6.0 * self.cohesion * cos_phi / (3.0 - sin_phi)
# Associated linear Drucker-Prager return denominator for q-based f.
self.return_denominator = 3.0 * self.mu + 9.0 * self.bulk * self.eta**2 + self.H
self.history: Dict[str, Dict[str, torch.Tensor]] = {}
for etype, trans in self.transformation.items():
n_elem = trans.n_elements
n_quad = trans.n_quadrature
eps_p = torch.zeros((n_elem, n_quad, 3, 3), device=self.device, dtype=self.dtype)
alpha = torch.zeros((n_elem, n_quad), device=self.device, dtype=self.dtype)
self.history[etype] = {"eps_p": eps_p, "alpha": alpha}
@staticmethod
def _small_strain_3d(graddisplacement: torch.Tensor) -> torch.Tensor:
"""Return the 3D small-strain tensor for 2D or 3D input gradients."""
dim = graddisplacement.shape[-1]
if dim == 2:
eps_2d = 0.5 * (graddisplacement + graddisplacement.transpose(-1, -2))
eps = torch.zeros((3, 3), device=graddisplacement.device, dtype=graddisplacement.dtype)
eps[:2, :2] = eps_2d
return eps
return 0.5 * (graddisplacement + graddisplacement.transpose(-1, -2))
def _elastic_stress(self, eps_e: torch.Tensor) -> torch.Tensor:
"""Tension-positive isotropic elastic stress from elastic strain."""
eye = torch.eye(3, device=eps_e.device, dtype=eps_e.dtype)
tr_eps = eps_e.diagonal(dim1=-2, dim2=-1).sum(-1)
dev_eps = eps_e - (tr_eps[..., None, None] / 3.0) * eye
return 2.0 * self.mu * dev_eps + self.bulk * tr_eps[..., None, None] * eye
@staticmethod
def _invariants(sigma: torch.Tensor) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
"""Return I1, deviatoric stress, and q for a tension-positive stress."""
eye = torch.eye(3, device=sigma.device, dtype=sigma.dtype)
I1 = sigma.diagonal(dim1=-2, dim2=-1).sum(-1)
s = sigma - (I1[..., None, None] / 3.0) * eye
s_contract_s = (s * s).sum(dim=(-2, -1))
q = torch.sqrt(torch.clamp(1.5 * s_contract_s, min=1.0e-30))
return I1, s, q
def element_energy(
self,
graddisplacement: torch.Tensor,
eps_p_n: torch.Tensor,
alpha_n: torch.Tensor,
) -> torch.Tensor:
"""Algorithmic incremental potential density at one quadrature point."""
eps = self._small_strain_3d(graddisplacement)
eps_e_trial = eps - eps_p_n
sigma_trial = self._elastic_stress(eps_e_trial)
I1_trial, _, q_trial = self._invariants(sigma_trial)
f_trial = q_trial + self.eta * I1_trial - (self.k + self.H * alpha_n)
dgamma = torch.clamp(f_trial, min=0.0) / self.return_denominator
tr_eps_e = eps_e_trial.diagonal(dim1=-2, dim2=-1).sum(-1)
eye = torch.eye(3, device=eps.device, dtype=eps.dtype)
dev_eps_e = eps_e_trial - (tr_eps_e / 3.0) * eye
elastic_energy = 0.5 * self.bulk * tr_eps_e**2 + self.mu * (dev_eps_e * dev_eps_e).sum()
return elastic_energy - 0.5 * self.return_denominator * dgamma**2
def update_state(self, u_vec: torch.Tensor) -> None:
"""Commit per-quadrature state after a converged load step."""
with torch.no_grad():
for etype, trans in self.transformation.items():
cells = trans.elements
u_elem = u_vec[cells]
grad_u = torch.einsum("bqkx,bku->bqux", trans.shape_grad, u_elem)
dim = grad_u.shape[-1]
if dim == 2:
eps = torch.zeros(grad_u.shape[:2] + (3, 3), device=u_vec.device, dtype=u_vec.dtype)
eps[..., :2, :2] = 0.5 * (grad_u + grad_u.transpose(-1, -2))
else:
eps = 0.5 * (grad_u + grad_u.transpose(-1, -2))
hist = self.history[etype]
eps_p_n = hist["eps_p"]
alpha_n = hist["alpha"]
eps_e_trial = eps - eps_p_n
sigma_trial = self._elastic_stress(eps_e_trial)
I1_trial, s_trial, q_trial = self._invariants(sigma_trial)
f_trial = q_trial + self.eta * I1_trial - (self.k + self.H * alpha_n)
dgamma = torch.clamp(f_trial, min=0.0) / self.return_denominator
q_safe = torch.clamp(q_trial, min=1.0e-30)
n_dev = 1.5 * s_trial / q_safe[..., None, None]
eye = torch.eye(3, device=u_vec.device, dtype=u_vec.dtype)
flow_dir = n_dev + self.eta * eye
active = (f_trial > 0.0).to(dtype=u_vec.dtype)
dgamma = dgamma * active
hist["eps_p"] += dgamma[..., None, None] * flow_dir
hist["alpha"] += dgamma
def element_data_from_history(self) -> Dict[str, Dict[str, torch.Tensor]]:
"""Return history in the element_data structure expected by energy()."""
return {
"eps_p_n": {etype: h["eps_p"] for etype, h in self.history.items()},
"alpha_n": {etype: h["alpha"] for etype, h in self.history.items()},
}
def mean_alpha(self) -> torch.Tensor:
"""Return the mean committed plastic multiplier across all quadrature points."""
values = [hist["alpha"].reshape(-1) for hist in self.history.values()]
return torch.cat(values).mean()
def max_alpha(self) -> torch.Tensor:
"""Return the maximum committed plastic multiplier."""
values = [hist["alpha"].reshape(-1) for hist in self.history.values()]
return torch.cat(values).max()
def mean_stress(self, u_vec: torch.Tensor) -> torch.Tensor:
"""Average committed Cauchy stress over elements and quadrature points."""
stresses = []
with torch.no_grad():
for etype, trans in self.transformation.items():
u_elem = u_vec[trans.elements]
grad_u = torch.einsum("bqkx,bku->bqux", trans.shape_grad, u_elem)
dim = grad_u.shape[-1]
if dim == 2:
eps = torch.zeros(grad_u.shape[:2] + (3, 3), device=u_vec.device, dtype=u_vec.dtype)
eps[..., :2, :2] = 0.5 * (grad_u + grad_u.transpose(-1, -2))
else:
eps = 0.5 * (grad_u + grad_u.transpose(-1, -2))
sigma = self._elastic_stress(eps - self.history[etype]["eps_p"])
stresses.append(sigma.reshape(-1, 3, 3))
return torch.cat(stresses, dim=0).mean(dim=0)