GENESIS框架与材料模型

Genesis 的核心是全新设计的物理引擎，整合了多个物理求解器到统一框架中。在此基础上，我们添加了生成代理系统，用于自动化数据生成，特别适合机器人等领域的应用。
《什么是 Genesis》

GENESIS中集成了多种物理引擎，包括物质点法(MPM)、有限单元法(FEM)、基于位置动力学法(PBD)、欧拉气体稳定流求解器(SF)和光滑粒子流体动力学(SPH)法。其中，MPM适合模型岩土材料等易断裂、大变形的物理过程；FEM适用于小变、高精度的物理仿真、PBD适用于模拟弹簧质点系统、SF和SPH适用于模拟流体、气体等的力学行为。

对于岩土工程来说，使用MPM法模拟土的物理状态具有显著的优势。但GENESIS框架中仅预制了几种常见的材料模型，包括弹性模型、一般弹塑性模型、砂模型、雪模型、液体模型、肌肉模型等，如下图所示，对岩土力学模型的支持效果较差。在岩土力学中，土的本构模型已经衍生出了几千种之多，不同的本构模型分别对应了材料的应力应变关系。在基坑工程分析中常用的本构模型有：摩尔库伦模型（MC）、剑桥模型（CC）、修正剑桥模型（MCC）、小应变模型（SS）、小应变硬化模型（HSS）等。其中，MC模型是较为经典的岩土工程模型，目前被广泛写入各类规范和教材。剑桥模型和修正剑桥模型适用于小变形软黏土分析，小应变模型和小应变硬化模型适用于岩土体与结构小变形相互作用的模拟。

为了提升 GENESIS 框架在岩土工程模拟中的适应性，我们计划扩展其 MPM 材料模型接口，以支持典型的岩土本构模型（如 Mohr-Coulomb、Cam-Clay、Modified Cam-Clay 等）。这一部分的核心目标，是在不破坏现有框架结构的前提下，实现可插拔、模块化的本构模型管理。

GENESIS MPM 材料模型机制解析

GENESIS 中的 MPM 材料模型被抽象为函数接口，由 MPMSolver 在 p2g() 阶段调用。其核心机制包括：

1. 材料注册机制：通过 add_material() 方法添加材料对象，并自动记录其本构行为：

   self._mats_update_F_S_Jp.append(material.update_F_S_Jp)
   self._mats_update_stress.append(material.update_stress)

   每种材料模型在加入后会被分配唯一的 mat_idx，与粒子的 particles_info[i].mat_idx 对应。

2. 材料调度逻辑：在 p2g() 中，每个活跃粒子都会根据 mat_idx 调用对应材料的 update_F_S_Jp 和 update_stress，实现差异化行为模拟。

3. 变量接口说明：

   • F_tmp: 当前时间步预测的形变梯度，由 (I + dt * C) @ F 计算得到。
   • F: 历史形变梯度；F_new 为更新后的结果。
   • U, S, V: 对 F_tmp 做 SVD 得到的分解结果，用于应变提取与旋转构造。
   • Jp: 体积压缩比，记录材料体积演化历史，常用于应变硬化。
   • J: 当前体积比（det(S)），用于判断材料变形程度与体积变化。
   • actu, m_dir: 可选的肌肉激活与方向参数，部分材料使用。

4. 应力张量构造约定：应力返回值必须是 Cauchy 应力（σ），且张量必须对称，以保证力矩守恒与动量守恒。

5. 存储机制与性能优化：材料模型本身无需存储状态，仅需依赖粒子场变量；这使得 GENESIS 可以在 GPU 上批量调用本构逻辑，实现高效并行计算。

6. 反向传播支持：GENESIS 支持 Taichi 自动微分。为保证材料模型能正确传播梯度，其 update_F_S_Jp 和 update_stress 函数中使用的操作应支持 Taichi autodiff（如 log, svd, matmul 等）。某些复杂逻辑可能需自定义反向函数或启用 ti.ad.grad_replaced()。

拓展材料的注意事项

1. 数值稳定性控制：本构模型中频繁涉及 SVD、log、逆矩阵等操作，应使用 ti.max, ti.min 等方式对其进行 clamp 或保护处理，避免数值发散。

2. 数值计算精度控制：在一些本构模型应力更新过程中涉及到塑性乘子$\Delta \lambda$的迭代计算，这些塑性乘子往往量级很小（$<=1e-6$），而更新塑性乘子过程中用到的雅格比行列式往往尺度极大（$>=1e12$），造成精度丢失等问题。编写此类模型时，执行初始化命令 gs.init() 应带上 precision=64 参数，同时考虑采用其他近似迭代方法绕过尺度差距过大的计算过程。

3. 保持力学一致性：本构模型必须返回对称应力张量，否则会破坏动量守恒；应变更新过程应保持物理约束（如不可逆变形的不可恢复性）。

4. 高阶模型性能开销：如 MCC、 HSS等模型内部含迭代过程，建议考虑近似解或 GPU-friendly 实现，避免过高计算开销。

5. 多材料混合场景：当不同粒子具有不同材料模型时，要确保 update_stress 和 update_F_S_Jp 能正确地被选取并稳定运行，必要时在外层逻辑中加入边界检测与模型切换判断。

6. 差异化自动微分支持：目前部分模型在 requires_grad=True 情况下无法稳定运行（如 CPIC），需谨慎验证其在反向传播中的梯度有效性。

7. 单位一致性与粒子体积缩放：注意 _p_vol_scale 的缩放效应，可能影响密度、质量、体积应变率等物理量的精度。

MPM Solver中的材料接口函数

MPM Solver中注册了两个材料更新方法：update_F_S_Jp 与 update_stress 函数，因此所有新定义的材料模型均需要在内部实现这两个方法。求解器中对这两个方法进一步规定了输入和输出参数，如下代码所示。

class MySoilModel:
    def update_F_S_Jp(self, J, F_tmp, U, S, V, Jp):
        # 根据本构关系计算更新后的 F, S, Jp
        # 比如：对 MCC 模型计算屈服面投影后的塑性应变增量
        return F_new, S_new, Jp_new

    def update_stress(self, U, S, V, F_tmp, F_new, J, Jp, actu, m_dir):
        # 构造基于当前应力状态的 Cauchy 应力张量
        # 比如使用弹塑性张量应力率公式：sigma = f(F, epsilon_p, Jp)
        return stress_tensor

这两个接口通过 MPMSolver.add_material() 进行注册，同时在 p2g() 中进行调度执行，确保材料行为与物理过程解耦。

修正剑桥模型的编码实现

修正剑桥模型介绍

剑桥模型是英国剑桥大学的Roscoe和Burland根据正常固结粘土和弱超固结粘土的三轴试验，采用状态边界面的概念，由塑性理论的流动法则和塑性势理论，采用简单曲线配合法，建立塑性与硬化定律的函数。它考虑了静水压力屈服特性、压硬性、剪缩性，但破坏面有尖角，该点的塑性应变方向不易确定。其假定的弹性墙内加载仍会产生塑性变形。

试验证明，对于正常固结粘土和弱固结的饱和重塑粘土，孔隙比e与外力p，q之间存在有唯一的关系，且不随应力路径而发生变化（这是通过p为常数的试验、固结排水试验和固结不排水试验，得出的试验破坏结果在p-q上的投影点连线斜率相等推导得出的），这个面叫作边界状态面。它可以看作由无数条不同应力比（p/q）的正常压缩曲线组成的。剑桥模型在e-p-q空间坐标的图形如图1。

而后他们又针对剑桥模型的弹头屈服面进行了修正得到修正剑桥模型(MCC)，使得屈服面的轨迹为椭圆，指出在完全状态边界面内，当剪应力增加时，虽不产生塑性体积变形，但产生塑性剪切变形。MCC模型主要适用于三轴加载的情况，对于软粘土的实验测量为该本构模型的建立提供了依据。该模型中孔隙比e（体积应变εv）是有效平均应力σ自然对数值的函数，如下公式所示：

$$ \lambda^*=\frac{\lambda}{1+e} $$ $$ \kappa^{*}=\frac{\kappa}{1+e} $$

其中：$\kappa=1,3\frac{1-\nu_{c}}{1+\nu}C_{z}$为回弹曲线斜率，$\lambda=\frac{C_c}{2,3}$为正常固结曲线（NCL）的斜率，$e$为当前孔隙比，$C_c$为侧限压缩实验得到的单向压缩指数，$C_s$为单向回弹指数。

MC与MCC模型在修正后的主要区别如下所示：

特征	原始剑桥模型	修正剑桥模型
屈服方程	$q = M p \sqrt{\frac{p_c}{p} - 1} )$	$\frac{q^2}{M^2 p^2} + p (p - p_c) = 0$
屈服面形状	非对称抛物线形	椭圆形，通过原点
几何特性	仅描述压缩行为，未覆盖拉伸区域	覆盖全应力空间（压缩与拉伸） （拉伸截断）
硬化参数	仅依赖塑性体积应变	依赖塑性体积应变与应力路径
临界状态线 (CSL)	未显式定义体积恒定条件	显式定义临界孔隙比 $e_c = e_\Gamma - \lambda \ln p$
临界状态应力比	固定值 $M$	引入状态参数 $\psi = e - e_c$，动态调整 $M$ （暂未在本示例代码中实现）
压缩指数 ($\lambda$)	仅用于正常固结线	结合临界状态线定义
回弹指数 ($\kappa$)	未显式分离弹性变形	显式分离弹性体积变化
适用土类	黏土（正常固结为主）	黏土、超固结土、部分砂土
数值稳定性	低应力水平下可能出现奇异	椭圆屈服面提高数值收敛性

上图包括一条正常固结曲线（NCL）和一组回弹曲线。在第一次加载时，原始岩土体沿着NCL曲线向下变化。然后，当岩土体固结压缩至某一应力水平（该应力水平又被称作先期固结压力pc）时对其进行卸载，并沿着当前回弹曲线上升。接下来，对岩土体进行再加载，使其沿着回弹曲线（再压缩曲线）下降，直到到达一个给定的应力状态pc，即卸载前岩土体的应力状态。这时，岩土体又开始沿着正常固结曲线向下变化（原始加载-压缩曲线）。

在离散的时间步上，屈服面的的硬化和软化规律由当前的先期固结压力$p_c$控制，如下图所示：

$$ p_{c}^{i+1}=p_{c}^{i}\exp\left[\frac{-\Delta\varepsilon_{v}^{pl}}{\lambda^{*}-\kappa^{*}}\right] $$

其中，$p_c^{i+1}$为当前的先期固结压力，$\Delta\varepsilon_{\nu}^{pl}$为塑型体积应变增量，$M_{cs}$为临界状态线的斜率，分别对于三轴压缩和三轴拉伸情况有不同的表达式：

$$ M_{cz}^{+30^{\circ}}\left(\varphi_{cv}\right)=\frac{2\sqrt{3}\sin\varphi_{cv}}{3-\sin\varphi_{cv}} $$ $$ M_{cz}^{-30^{\circ}}\left(\varphi_{cv}\right)=\frac{2\sqrt{3}.\sin\varphi_{cv}}{3+\sin\varphi_{cv}} $$

其中，$\phi_{cv}$是临界状态下的内摩擦角。

修正剑桥模型结构设计

在 GENESIS 框架中实现修正剑桥模型（Modified Cam-Clay, MCC）时，需围绕两个关键函数 update_F_S_Jp() 与 update_stress() 完成核心材料行为逻辑。这两个函数分别负责：

• update_F_S_Jp()：在每个时间步内，根据当前状态计算新的形变梯度 F、奇异值 S（用于应变）、体积压缩比 Jp，并执行应力点的投影与状态变量更新。
• update_stress()：根据更新后的形变与材料当前状态计算柯西应力张量，并在必要时执行临界状态修正、塑性体积修正等操作。

在 GENESIS 中，MCC 材料类需继承自 Base 材料基类，并在初始化时设定模型核心参数（如临界状态线斜率 M、固结曲线斜率 λ、回弹曲线斜率 κ、初始先期固结压力 pc0、初始孔隙比 e0 等）。此外，还需预定义若干控制参数以保证数值稳定性，例如：

• min_pc：最小先期固结压力，防止屈服面退化
• min_p_prime：最小平均有效应力，防止 SVD 退化
• tolerance：屈服面投影的残差容忍阈值
• max_iter：内点迭代的最大迭代次数

结构上，模型类应实现如下几个重要子函数，其中涉及到的变量解释如下表所示：

变量解释表

符号	物理意义
$p$	有效平均主应力（pressure）
$q$	偏应力（deviatoric stress），代表剪切强度
$M$	临界状态线斜率，定义临界状态条件 $ q = M p $
$p_c$	先期固结压力，控制屈服面的硬化或软化
$J$	当前总体积比，等于 $ \det(F) $
$J_p$	累积塑性体积比，记录体积压缩历史
$\epsilon_v^{pl}$	塑性体积应变增量
$\lambda$	正常固结曲线斜率（压缩曲线斜率）
$\kappa$	回弹曲线斜率
$e$	孔隙比（void ratio）
$e_c$	临界状态孔隙比
$E_e$	弹性 Green-Lagrange 应变张量
$F$, $F_p$, $F_e$	总变形、塑性变形与弹性变形梯度张量
$S$	SVD 分解的奇异值张量（主伸长率）
$K$	当前材料体积模量
$G$	当前材料剪切模量
$\Delta \lambda$	塑性乘子，用于控制塑性流动的强度

1. compute_moduli(p, e)

根据当前有效应力与孔隙比，计算弹性模量 K 与剪切模量 G：

$$ K = \frac{(1 + e) \cdot p}{\kappa}, \quad G = \frac{3K(1 - 2\nu)}{2(1 + \nu)} $$

2. yield_function(p, q, pc)

定义 MCC 模型的屈服面函数：

$$ f(p, q, p_c) = \frac{q^2}{M^2 p^2} + p(p - p_c) $$

3. df_dp, df_dq

屈服函数对 p 与 q 的偏导数：

$$ \frac{\partial f}{\partial p} = (2p - p_c) - \frac{2q^2}{M^2 p^3}, \quad \frac{\partial f}{\partial q} = \frac{2q}{M^2 p^2} $$

4. _binary_iteration()

常规的牛顿莱布尼茨迭代法一般是通过计算目标函数的梯度，后逐步迭代直到模型收敛，如下公式所示。

$$ \Delta \lambda^{(k+1)} = \Delta \lambda^{(k)} - \frac{f(\Delta \lambda^{(k)})}{J(\Delta \lambda^{(k)})} $$

其中$f$是残差函数（屈服函数），$\frac{\partial f}{\partial \Delta \lambda}$是雅可比矩阵。然而，在修正剑桥模型的塑性修正中，直接应用牛顿法可能遇到以下问题：

1. 雅可比矩阵的病态性：修正剑桥模型的雅可比矩阵$J$由三部分组成：

   $$ J=\underbrace{-K\left(\frac{\partial f}{\partial p^{\prime}}\right)^2}_{\text{体积项}}\underbrace{-3G\left(\frac{\partial f}{\partial q}\right)^2}_{\text{偏量项}}\underbrace{-\frac{(1+e_0)}{\lambda-\kappa}p^{\prime}p_c\frac{\partial f}{\partial p^{\prime}}}_{\text{硬化项}} $$

   当塑性应变增量较大时，硬化项可能主导雅可比矩阵，导致$J$的值剧烈变化，矩阵条件数恶化。若参数标定不当（如$\lambda \approx \kappa$，硬化项分母接近零，引发数值不稳定。

2. 导数计算的精度差：时间离散情况下导数所有导数计算过程均存在误差，而牛顿法需精确计算雅可比矩阵，而修正剑桥模型的导数项包含非线性耦合项（如$\partial f/\partial p^{\prime}$和$\partial f/\partial q$），公式复杂易错。

3. 局部收敛性：牛顿法仅在初始猜测接近真解时保证收敛。若试探应力$p_{\mathrm{trial}}^{\prime}$远离屈服面，初始$\Delta\lambda$猜测不当，可能陷入局部极小或发散。
4. 函数不连续与不可导：有限精度计算中，屈服函数$f$的导数因舍入误差出现不连续。且代码中的 if s_norm > 1e-10，在偏应力接近零时，导数$\partial f/\partial q$会突变，破坏牛顿法的连续性假设。

而直接使用二分迭代法能够直接解决以上1、3和4中涉及到的全部或部分问题。但MCC模型在材料硬化修正过程中材料的参数如压缩模量$K$和剪切模量$G$会直接随着有效主应力和孔隙比的连续变化，近似发生连续变化，而在二分迭代时同时考虑大量参数的随动变化将会使得搜索结果不断波动，造成难以收敛。在时间积分过程中，减小$\Delta t$的值能够较好地减小此类误差。

这里使用二分法迭代搜索满足屈服条件的 Δλ，使应力状态落在屈服面上，并根据 Δλ 更新：

$$ \Delta \varepsilon_v^{pl} = \Delta \lambda \cdot \frac{\partial f}{\partial p} $$ $$ p' = p - K \cdot \Delta \varepsilon_v^{pl}, \quad q' = q - 3G \cdot \Delta \lambda \cdot \frac{\partial f}{\partial q} $$ $$ p_c' = p_c \cdot \exp\left[\frac{(1 + e)\cdot \Delta \varepsilon_v^{pl}}{\lambda - \kappa}\right] $$

5. update_F_S_Jp()

该函数主要流程为：

• 基于当前 C 计算 F_tmp（预测形变梯度）
• 对 F_tmp 执行 SVD 得到 U, S, V
• 构造试应力状态（p_trial, q_trial）
• 判断是否进入塑性：若屈服函数 > 0，则执行 _binary_iteration 搜索 Δλ
• 更新变量：

$$ J_p^{\text{new}} = J_p \cdot \exp(\Delta \varepsilon_v^{pl}) $$ $$ S^{\text{new}} = S \cdot J_p^{-{1}/{3}}, \quad F^{\text{new}} = U \cdot S^{\text{new}} \cdot V^T $$

6. update_stress()

用于计算应力张量：

• 弹性应变：$Ee = \frac{1}{2}(C_e - I), \quad C_e = F e^T F_e$

• 平均有效应力与偏应力：$p’ = K (J_e - 1), \quad s = 2G \left(E_e - \frac{1}{3} \text{tr}(E_e) I \right)$

• 柯西应力张量：$\sigma = p’I + s$

• 临界状态修正：

若：$\left| \frac{q}{p’} - M \right| < \varepsilon$，则执行：

$$ s \leftarrow s \cdot \frac{Mq}{q}, \quad \sigma = p'I + s $$

• 应力体积修正：$\sigma_{\text{true}} = \left(\frac{J_e}{J_{\text{total}}}\right) \cdot \sigma$

整体计算流程和结构如下图所示。通过上述结构，MCC 模型可嵌入 GENESIS 的 MPM 流程中实现状态变量更新与应力响应，适用于软土剪切与体积变形响应的精准建模。

flowchart TD
    A0["Start: 当前时间步"] --> A1["预测形变梯度 F_tmp = (I + dt*C) * F"]
    A1 --> A2["奇异值分解: F_tmp = U * S * Vᵗ"]
    A2 --> A3["试应变张量 Eₑ = 0.5*(Cₑ - I)"]
    A3 --> A4["计算 p', q_trial, f_trial"]
    A4 --> B1{"屈服函数 f_trial > tol?"}

    B1 -- 否 --> C1["弹性更新: F_new = F_tmp"]
    C1 --> C2["S_new = S; Jp_new = Jp"]
    C2 --> D1["更新 e, pc, εₚ"]

    B1 -- 是 --> E1["塑性投影: 启动二分迭代"]
    E1 --> E2["搜索 Δλ 使 f ≈ 0"]
    E2 --> E3["更新 p', q, pc"]
    E3 --> E4["计算 Δεᵥᵖ = Δλ*∂f/∂p"]
    E4 --> E5["Jp_new = Jp*exp(Δεᵥᵖ)"]
    E5 --> E6["S_new = S*Jp_new^(-1/3)"]
    E6 --> E7["F_new = U*S_new*Vᵗ"]
    E7 --> D1

    D1 --> Z2["调用 update_stress"]
    Z2 --> Z3["构造 σ = p'*I + s"]
    Z3 --> Z4{"是否 q/p' ≈ M?"}
    Z4 -- 是 --> Z5["修正 q = M*p', 缩放 s"]
    Z4 -- 否 --> Z6["保持当前 σ"]
    Z5 --> Z7["缩放 σ: σ = σ*(Jₑ / J_total)"]
    Z6 --> Z7
    Z7 --> End["返回 σ"]

修正剑桥模型的代码实现

基于以上的理论结构，完整的实现代码如下所示：

import taichi as ti
import genesis as gs
from math import pi as M_PI
from failureState import FailureState
from base import Base
import numpy as np

@ti.data_oriented
class CamClay(Base):
    """Cam-Clay model for MPM
    """
    def __init__(self, E=100e6, nu=0.3, rho=1000, lam=None, mu=None, 
                 sampler="pbs",
                 lambda_c=0.25,
                 kappa=0.05,
                 M=1.2,
                 pc0=400e3,
                 alpha0=0,
                 beta=0.5,
                 n=1.5,
                 e_c=1.5,
                 e0=0.8,
                 tension_cutoff=0.0,
                 tolerance=1e-6, 
                 min_pc0=5,
                 min_K=1e3,
                 min_G=1e3,
                 min_pc=1,
                 min_p_prime=2,
                 max_iter=10):
        super().__init__(E, nu, rho, lam, mu, sampler)
        self._M = M
        self._lambda_c = lambda_c
        self._kappa = kappa
        self._e0 = e0
        self._alpha0 = alpha0
        self._beta = beta
        self._n = n
        self._e_c = e_c
        self._pc0 = pc0
        self._K = ti.field(gs.ti_float, ())
        self._G = ti.field(gs.ti_float, ())
        self._tolerance = tolerance
        self._max_iter = max_iter
        self._tension_cutoff = tension_cutoff
        self._stress_update = False
        self._min_pc0 = min_pc0
        self._min_K = min_K
        self._min_G = min_G
        self._min_pc = min_pc
        self._min_p_prime = min_p_prime

        self._pc = ti.field(gs.ti_float, ())
        self._epsilon_p = ti.field(gs.ti_float, ())
        self._e = ti.field(gs.ti_float, ())
        self._state_psi = ti.field(gs.ti_float, ())

        @ti.kernel
        def init_():
            self.init_state_variables()
            self.compute_moduli(p_prime=self._pc0, e_current=self._e0)        
        init_()

    @ti.func
    def init_state_variables(self):

        self._pc[None] = self._pc0
        self._epsilon_p[None] = 0
        self._e[None] = self._e0
        self._state_psi[None] = self._e[None] - self._e_c


    @ti.func
    def compute_yield_function(self, p, q, pc):
            term =  1.0 - self._beta * ti.pow(p / pc, self._n)
            f = q - self._M * p * term
            return f
    
    @ti.func
    def yield_function(self, p, q, pc):
        return (q * q) / (self._M * self._M * p * p) + p * (p - pc)
    
    @ti.func
    def df_dp(self, p, q, pc):
        return (2.0 * p - pc) - (2 * q * q) / (self._M * self._M * p * p * p)
    
    @ti.func
    def df_dq(self, p, q):
        return (2 * q) / (self._M * self._M * p * p)

    @ti.func
    def compute_moduli(self, p_prime, e_current):
        v = 1.0 + e_current
        # avoid all parametes are zero
        
        pc = max(p_prime, self._min_pc0)
        self._K[None] = max((v * pc) / self._kappa, self._min_K)
        self._G[None] = max((3 * self._K[None] * (1 - 2 * self.nu)) / (2 * (1 + self.nu)), self._min_G)

    @ti.func
    def _saturate_pc(self, pc):
        return ti.max(pc, self._tolerance)

    @ti.func
    def _binary_iteration(self, p_prime, q, pc, upper=0.0005, lower=1e-6):
        delta_lambda_min = lower
        delta_lambda_max = upper
        e = self._e[None]

        resi_upper = self._update_status_deltalambda(delta_lambda_max, p_prime, q, pc)
        resi_lower = self._update_status_deltalambda(delta_lambda_min, p_prime, q, pc)
        if resi_upper * resi_lower > 0:
            print("⚠️ 二分区间端点未包围根，建议检查初始试应力是否处于屈服面外")

        for i in range(self._max_iter):
            dl_mid = 0.5 * (delta_lambda_min + delta_lambda_max)
            res_mid = self._update_status_deltalambda(dl_mid, p_prime, q, pc)

            if ti.abs(res_mid) < 1e-3:
                break

            if self._update_status_deltalambda(delta_lambda_min, p_prime, q, pc) * res_mid < 0:
                delta_lambda_max = dl_mid
            else:
                delta_lambda_min = dl_mid

        delta_lambda = 0.5 * (delta_lambda_min + delta_lambda_max)
        delta_epsilon_p_vol = delta_lambda * self.df_dp(p_prime, q, pc)
        p_new = p_prime - self._K[None] * delta_epsilon_p_vol
        q_new = q - 3 * self._G[None] * delta_lambda * self.df_dq(p_prime, q)
        pc_new = pc * ti.exp((1 + self._e[None]) * delta_epsilon_p_vol / (self._lambda_c - self._kappa))

        return p_new, q_new, pc_new, delta_lambda, delta_epsilon_p_vol
    
    @ti.func
    def _update_status_deltalambda(self, delta_lambda, p_prime, q, pc):
        delta_epsilon_p_vol = delta_lambda * self.df_dp(p_prime, q, pc)
        p_new = p_prime - self._K[None] * delta_epsilon_p_vol
        q_new = q - 3 * self._G[None] * delta_lambda * self.df_dq(p_prime, q)
        pc_new = pc * ti.exp((1 + self._e[None]) * delta_epsilon_p_vol / (self._lambda_c - self._kappa))
        return self.yield_function(p_new, q_new, pc_new)

    @ti.func
    def update_F_S_Jp(self, J, F_tmp, U, S, V, Jp):
        if Jp == 0:
            Jp = 1.0
        J_e = F_tmp.determinant()
        C_e_trial = F_tmp.transpose() @ F_tmp
        E_e_trial = 0.5 * (C_e_trial - ti.Matrix.identity(gs.ti_float, 3))
        p_prime_trial = max(-E_e_trial.trace() / 3.0 * self._K[None], self._min_p_prime)

        self.compute_moduli(p_prime=p_prime_trial, e_current=self._e[None])

        dev_E_e_trial = E_e_trial - E_e_trial.trace() / 3 * ti.Matrix.identity(gs.ti_float, 3)
        s_trial = 2 * self._G[None] * dev_E_e_trial
        q_trial = ti.sqrt(1.5) * s_trial.norm()

        pc = self._pc[None]
        f_trial = (q_trial * q_trial) / (self._M* self._M * p_prime_trial * p_prime_trial
                                         ) + p_prime_trial * (p_prime_trial - pc)

        delta_epsilon_p_vol = 0.0
        F_new = F_tmp
        __, S_new, __ = ti.svd(F_tmp)
        Jp_new = Jp
        sigma_new = ti.Matrix.zero(gs.ti_float, 3, 3)
        if f_trial > self._tolerance:
            # binary iteration
            delta_lambda = 0.0
            q = q_trial
            p_prime = p_prime_trial
            p_prime, q, pc, delta_lambda, delta_epsilon_p_vol = self._binary_iteration(p_prime, q, pc)
            

            self._pc[None] = self._saturate_pc(pc)
            Jp_new = Jp * ti.exp(delta_epsilon_p_vol)
            Jp_new = ti.max(Jp_new, 1e-6)

            S_new = S * ti.pow(Jp_new, -1.0 / 3.0)
            F_new = U @ S_new @ V.transpose()
            self._epsilon_p[None] += delta_epsilon_p_vol

        delta_e = (self._lambda_c - self._kappa) * ti.log(pc / pc)
        self._e[None] = self._e0 - (1.0 + self._e0) * delta_epsilon_p_vol + delta_e
        self.update_psi(self._e[None])

        return F_new, S_new, Jp_new

    @ti.func
    def update_stress(self, U, S, V, F_tmp, F_new, J, Jp, actu, m_dir):

        F_e = F_new 
        F_p = self.get_plastic_deformation()
        
        C_e = F_e.transpose() @ F_e
        J_e = F_e.determinant() 
        E_e = 0.5 * (C_e - ti.Matrix.identity(ti.f32, 3))
        
        p_prime = self._K[None] * (J_e - 1.0)
        dev_E_e = E_e - (E_e.trace() / 3.0) * ti.Matrix.identity(ti.f32, 3)
        s = 2.0 * self._G[None] * dev_E_e 
        
        sigma = p_prime * ti.Matrix.identity(ti.f32, 3) + s
        
        q = ti.sqrt(1.5) * s.norm()

        if ti.abs(q / p_prime - self._M) < 1e-6 and p_prime > 0:
            q_target = self._M * p_prime
            if s.norm() > 1e-10:
                s_scaled = (q_target / q) * s  # 缩放偏应力
                sigma = p_prime * ti.Matrix.identity(ti.f32, 3) + s_scaled

        F_total = F_e @ F_p
        J_total = F_total.determinant()
        
        if ti.abs(J_total) > 1e-10:
            sigma *= (J_e / J_total)
        
        for i in ti.static(range(3)):
            if sigma[i, i] < -self._tension_cutoff:
                sigma[i, i] = -self._tension_cutoff
        
        return sigma
    
    @ti.func
    def get_plastic_deformation(self):
        F_p = ti.Matrix.identity(ti.f32, 3)
        F_p[0, 0] = ti.exp(-self._epsilon_p[None])
        F_p[1, 1] = ti.exp(-self._epsilon_p[None])
        F_p[2, 2] = ti.exp(2 * self._epsilon_p[None])  
        return F_p

单元测试

笔者为材料模型中的每个函数都设计了相应的单元测试：

import math
import unittest
import taichi as ti
import numpy as np
import genesis as gs  # gs.EPS 等常量
from cam_clay import CamClay
from failureState import FailureState

gs.init(precision="64")

class TestModifiedCamClay(unittest.TestCase):


    def setUp(self):
        ti.init(arch=ti.cpu, debug=True, default_fp=ti.f64)  # 初始化Taichi

        self.material = CamClay(
            E=1e6, nu=0.3, rho=1000, 
            M=1.2, lambda_c=0.15, kappa=0.05,
            pc0=100, beta=0.8, n=1.5, e_c=1.2
        )
        self.tolerance = 1e-6
    
    def test_initialization(self):
        """测试参数初始化是否正确"""
        self.assertAlmostEqual(self.material._M, 1.2)
        self.assertAlmostEqual(self.material._lambda_c, 0.15)
        self.assertAlmostEqual(self.material._pc0, 100)
        self.assertAlmostEqual(self.material._tolerance, 1e-6)

    def test_state_variables_initialization(self):
        """测试状态变量初始化"""
        @ti.kernel
        def test_init_vars() -> gs.ti_float:
            self.material.init_state_variables()
            return self.material._pc[None]

        pc = test_init_vars()
        self.assertAlmostEqual(pc, 100.0)
    
    def test_yield_function_calculation(self):
        """测试屈服函数计算"""
        @ti.kernel
        def compute_f(p: gs.ti_float, q: gs.ti_float, pc: gs.ti_float) -> gs.ti_float:
            return self.material.compute_yield_function(p, q, pc)
        
        # 测试弹性状态 (f < 0)
        f1 = compute_f(100.0, 50.0, 200.0)
        self.assertLess(f1, 0.0)
        
        # 测试屈服状态 (f = 0)
        term = 1.0 - 0.8 * pow(150/200, 1.5)
        expected_f = 180.0 - 1.2*150.0*term
        f2 = compute_f(150.0, 180.0, 200.0)
        self.assertAlmostEqual(f2, expected_f, delta=0.0001)

    def test_moduli_calculation(self):
        """测试弹性模量计算"""
        @ti.kernel
        def test_moduli(p: gs.ti_float, e: gs.ti_float):
            self.material.compute_moduli(p, e)
        
        test_moduli(100.0, 1.0)
        v = 1.0 + 1.0
        expected_K = (v * 100) / 0.05
        expected_G = (3 * expected_K * (1 - 2*0.3)) / (2*(1+0.3))
        
        self.assertAlmostEqual(self.material._K[None], expected_K, delta=self.tolerance)
        # self.assertAlmostEqual(self.material._G[None], expected_G, delta=self.tolerance)
        self.assertAlmostEqual(self.material._G[None], expected_G, delta=1e-4)

    def test_elastic_response(self):
        """测试弹性响应情况"""
        @ti.kernel
        def test_update() -> ti.types.matrix(3,3,gs.ti_float):
            F_tmp = ti.Matrix.identity(gs.ti_float, 3)
            S = ti.Matrix.identity(gs.ti_float, 3)
            V = ti.Matrix.identity(gs.ti_float, 3)  # 代替 None
            U = ti.Matrix.identity(gs.ti_float, 3)  # 代替 None
            return self.material.update_F_S_Jp(1.0, F_tmp, U, S, V, 0.0)[1]
        
        S = test_update()
        # 弹性响应时应变无变化
        self.assertTrue(np.allclose(S.to_numpy(), np.eye(3), atol=self.tolerance))

    def test_plastic_response(self):
        """测试塑性响应情况"""
        @ti.kernel
        def test_update(F_tmp: ti.types.matrix(3,3,gs.ti_float)) -> ti.types.matrix(3,3,gs.ti_float):
            U, S, V = ti.svd(F_tmp)
            _, S_new, _ = self.material.update_F_S_Jp(1.0, F_tmp, U, S, V, 0.0)
            return S_new

        # 构造产生塑性变形的变形梯度
        F_plastic = ti.Matrix([[1, 0.05, 0.05], [0.0, 1, 0.0], [0.0, 0.0, 1]])
        S = test_update(F_plastic)

        self.assertFalse(np.allclose(S.to_numpy(), np.zeros((3,3)), atol=self.tolerance))
        self.assertTrue(np.allclose(S.to_numpy(), S.to_numpy().T, atol=self.tolerance))

    def test_stress_update(self):
        """测试应力更新逻辑"""
        @ti.kernel
        def test_stress() -> ti.types.matrix(3,3,gs.ti_float):
            U = ti.Matrix.identity(gs.ti_float, 3)
            S = ti.Matrix.zero(gs.ti_float, 3, 3)
            V = ti.Matrix.identity(gs.ti_float, 3)
            F_tmp = ti.Matrix([[1.05, 0, 0], [0, 0.95, 0], [0, 0, 0.98]])
            return self.material.update_stress(U, S, V, F_tmp, F_tmp, 1.0, 0.0, 0.0, 0)
        
        stress = test_stress().to_numpy()
        # 验证主应力分量符号
        self.assertGreater(stress[0,0], 0.0)  # 压缩为正
        self.assertEqual(stress[1,1], 0.0)      # 拉伸截断

if __name__ == '__main__':
    unittest.main()

小结

本文系统梳理了 GENESIS 框架中材料模型的调用机制、扩展路径与实现规范，并以修正剑桥模型（Modified Cam-Clay, MCC）为例，从本构理论、求解器机制、函数接口、数值实现与可视化流程图等多个角度，对该模型在物质点法（MPM）中的嵌入与运行过程进行了深入分析与逐步实现。作者理论和实践水平相当有限，还请各位老师同学们多多指正。

通过代码中的 update_F_S_Jp() 与 update_stress() 两大核心接口，本文构建了一个既满足塑性硬化规律、又适应 GPU 并行和自动微分的 MCC 材料模块，并成功融合于 GENESIS 的多物理场仿真架构中。然而，MCC 模型的其他流动准则、迭代收敛性、参数标定灵敏性与张量稳定性，计算框架的计算精度不足、实际计算任务的显存占用过大仍是实际应用中的挑战。