岩土力学:Drucker-Prager 条形基础

本示例把两个岩土力学构件结合成一个非线性边值问题:TensorMesh 的公共 Drucker-Prager 材料模型和弹性条形基础设置。一个矩形土块由中心条形基础加载,基础压力按载荷步逐级增加,每个收敛步后提交逐求积点的 Drucker-Prager 历史。

本构模型是内置的 DruckerPragerPlasticity 装配器,配以 FrictionalMaterial。TensorMesh 内部固体力学约定保持应力拉为正。为岩土力学报告,沉降以向下为正显示,

\[s = -u_y.\]

这是一个紧凑的教学示例,而非基础设计方法。

问题

该模型是二维平面应变式土块,面外厚度为单位 1。本构模型是带线性各向同性硬化的小应变、关联 Drucker-Prager 塑性,内部以拉为正的屈服函数写成

\[f(\sigma, \alpha) = q + \eta I_1 - (k + H\alpha) \le 0, \qquad I_1 = \mathrm{tr}(\sigma), \qquad q = \sqrt{\tfrac{3}{2}\,s:s}.\]

边界条件

本示例复用弹性基础的滚动支座设置:

• 底部边界:竖向位移固定,``u_y = 0``;

• 左右边界:水平位移固定,``u_x = 0``;

• 顶部边界:除受载基础块外自由。

基础压力被集总到基础块内的顶面节点上。

求解器

由于问题现在是路径相关的,它通过带载荷步进的非线性能量最小化来求解。在每个载荷步,用 L-BFGS 在自由位移自由度上最小化总势能 internal - external。该示例遵循与 drucker_prager_triaxial 示例相同的逐求积点历史生命周期:

• 逐求积点历史变量存储在 self.history[etype] 中;

• 上一步的 eps_p 和 alpha 通过 element_data 传入;

• update_state(u) 在 torch.no_grad() 下、每个收敛载荷步后调用。

合理性检验

脚本报告最终基础沉降、最大沉降、提交的塑性历史的最大值与均值,以及塑性质心。

相应的测试检查:基础向下沉降、提交的塑性历史与沉降随载荷单调增长、塑性区在基础下方局部化、低载荷/高黏聚力情形基本保持弹性并与弹性基础沉降吻合,以及黏聚力更高的土体发展出更少的塑性。

图 58 drucker_prager_footing.py 的输出。左面板显示按沉降 -u_y 着色的变形土网格。中面板显示提交的 Drucker-Prager 塑性历史变量,它在基础下方发展。右面板显示载荷-沉降曲线,包括塑性累积时的非线性增长。为便于观察,变形被放大。

运行它

cd examples/solid/geomechanics/drucker_prager_footing
python drucker_prager_footing.py

若只想快速做纯数值运行而不写出图:

python drucker_prager_footing.py --no-plot --steps 6 --chara-length 0.60

核心实现

载荷步进的非线性求解器是本示例的核心。它使用内置 DruckerPragerPlasticity 装配器的 energy、element_data_from_history 和 update_state 方法。

def solve_drucker_prager_footing(
    problem: FootingProblem,
    material: FrictionalMaterial,
    n_steps: int = 10,
    dtype: torch.dtype = torch.float64,
    device: str | torch.device = "cpu",
    lbfgs_max_iter: int = 50,
) -> Dict[str, Any]:
    """Solve the footing problem with Drucker-Prager plasticity and load stepping.

    The total potential energy ``internal - external`` is minimized over the free
    displacement DOFs with L-BFGS at each load step.  After each converged step
    the per-quadrature Drucker-Prager history is committed with ``update_state``.
    """
    device = torch.device(device)
    mesh = build_mesh(problem, dtype=dtype, device=device)
    dp = DruckerPragerPlasticity.from_mesh(mesh, material=material)

    dim = mesh.dim
    n_points = mesh.n_points
    n_dof = n_points * dim

    fixed = make_boundary_mask(mesh, problem)
    rhs_full, loaded_nodes, total_vertical_load = make_load_vector(mesh, problem)
    rhs_full_flat = rhs_full.flatten()

    free_flat = (~fixed).flatten()
    free_idx = torch.where(free_flat)[0]
    base_zeros = torch.zeros(n_dof, dtype=dtype, device=device)

    free_u = torch.zeros(free_idx.numel(), dtype=dtype, device=device, requires_grad=True)

    def recover_full(free_values: torch.Tensor) -> torch.Tensor:
        """Scatter the free DOFs into a full [n_points, dim] displacement field."""
        u_flat = base_zeros.index_add(0, free_idx, free_values)
        return u_flat.reshape(n_points, dim)

    load_factors: List[float] = []
    pressures_kpa: List[float] = []
    footing_settlements: List[float] = []
    max_settlements: List[float] = []
    max_alphas: List[float] = []
    mean_alphas: List[float] = []

    for step in range(1, n_steps + 1):
        load_factor = step / n_steps
        f_ext_flat = load_factor * rhs_full_flat

        optimizer = torch.optim.LBFGS(
            [free_u],
            max_iter=lbfgs_max_iter,
            tolerance_grad=1.0e-10,
            tolerance_change=1.0e-14,
            history_size=50,
            line_search_fn="strong_wolfe",
        )

        def closure() -> torch.Tensor:
            optimizer.zero_grad()
            u_full = recover_full(free_u)
            internal = dp.energy(
                point_data={"displacement": u_full},
                element_data=dp.element_data_from_history(),
            )
            external = torch.dot(f_ext_flat, u_full.flatten())
            loss = internal - external
            loss.backward()
            return loss

        optimizer.step(closure)

        with torch.no_grad():
            u_full = recover_full(free_u)
            dp.update_state(u_full)

            settlement = -u_full[:, 1]
            load_factors.append(float(load_factor))
            pressures_kpa.append(float(load_factor * problem.footing_pressure / 1.0e3))
            footing_settlements.append(float(settlement[loaded_nodes].mean()))
            max_settlements.append(float(settlement.max()))
            max_alphas.append(float(dp.max_alpha()))
            mean_alphas.append(float(dp.mean_alpha()))

    with torch.no_grad():
        u_final = recover_full(free_u).detach()

    return {
        "mesh": mesh,
        "dp": dp,
        "u": u_final,
        "fixed": fixed,
        "loaded_nodes": loaded_nodes,
        "total_vertical_load_N_per_m": total_vertical_load,
        "load_factors": load_factors,
        "pressures_kpa": pressures_kpa,
        "footing_settlements_m": footing_settlements,
        "max_settlements_m": max_settlements,
        "max_alphas": max_alphas,
        "mean_alphas": mean_alphas,
        "n_nodes": int(n_points),
        "n_steps": int(n_steps),
    }

下一步

一个自然的后续是为公共 Drucker-Prager 模型加入非关联流动和孔压耦合,或随着更多内置岩土力学材料的加入,用此基础几何复用它们。