Newton 柔性实体单元的代码位于 kernels_particle.py 的 eval_tetrahedra() 中。具体代码这里不完整展开，读者可以直接查看 Newton 的源码。这篇文章主要想回答一个问题：
Newton 到底用了什么方法实现柔性实体仿真？它和传统有限元有什么关系？
在传统有限元软件中，一个几何非线性实体单元通常按照如下流程计算：
\[\boldsymbol u
\rightarrow
\boldsymbol F
\rightarrow
\text{constitutive model}
\rightarrow
\boldsymbol \sigma
\rightarrow
\boldsymbol f_{\mathrm{int}}
\rightarrow
\boldsymbol K_{\mathrm T}.
\]
Newton 中的柔性四面体代码看起来却简单得多。在 eval_tetrahedra() 中几乎看不到传统有限元代码常见的 Gauss 积分、\(B\) 矩阵、材料 Jacobian 或几何刚度。它的核心计算过程可以概括为
\[\boxed{
\boldsymbol x
\rightarrow
\boldsymbol F
\rightarrow
(I_C,J)
\rightarrow
\boldsymbol P
\rightarrow
\boldsymbol f
}
\]
下面就从有限元的角度拆解这段代码，看看这些熟悉的有限元概念在 Newton 中是如何被组织起来的。
1. Linear Tet 与变形梯度
Newton 使用四节点线性四面体（Linear Tet）作为柔性实体单元。
![image-20261005101215710]()
记参考构形节点坐标为
\[\boldsymbol X_0,\quad
\boldsymbol X_1,\quad
\boldsymbol X_2,\quad
\boldsymbol X_3,
\]
当前构形坐标为
\[\boldsymbol x_0,\quad
\boldsymbol x_1,\quad
\boldsymbol x_2,\quad
\boldsymbol x_3.
\]
以节点 0 为基点，定义参考构形边矩阵
\[\boldsymbol D_m^0=
\begin{bmatrix}
\boldsymbol X_1-\boldsymbol X_0 &
\boldsymbol X_2-\boldsymbol X_0 &
\boldsymbol X_3-\boldsymbol X_0
\end{bmatrix},
\]
以及当前构形边矩阵
\[\boldsymbol D_s=
\begin{bmatrix}
\boldsymbol x_1-\boldsymbol x_0 &
\boldsymbol x_2-\boldsymbol x_0 &
\boldsymbol x_3-\boldsymbol x_0
\end{bmatrix}.
\]
对于 Linear Tet，单元内部的变形梯度 \(\boldsymbol F\) 为常数，并满足（推导见附录）
\[\boldsymbol F = \boldsymbol D_s(\boldsymbol D_m^0)^{-1}
\]
Newton 中对应的代码形式为
Ds = wp.matrix_from_cols(x10, x20, x30)
Dm = pose[tid]
F = Ds * Dm
其中 pose[tid] 保存的是参考构形边矩阵的逆 \(\boldsymbol D_m=(\boldsymbol D_m^0)^{-1}\)。
由于参考构形在仿真过程中保持不变，因此这个矩阵可以在初始化阶段预先计算。
从传统有限元角度看，这其实就是熟悉的 \(\boldsymbol F = \sum_a\boldsymbol x_a\otimes\nabla_X N_a\)。四节点线性 Tet 的形函数 \(N_a\) 关于参考坐标是线性的，因此 \(\nabla_XN_a=\text{constant}\)。
于是
\[\boxed{
\boldsymbol F=\text{constant}
}
\]
在整个单元内部成立。这也是为什么 Newton 的 Tet kernel 中看不到复杂的 Gauss 积分循环。
对于超弹性材料，内部虚功可以写成
\[\delta W_{\mathrm{int}} = \int_{\Omega_0} \boldsymbol P:\delta\boldsymbol F\,dV_0
\]
其中 \(\boldsymbol P\) 为第一 Piola-Kirchhoff 应力。
节点 \(a\) 对应的内力为
\[\boldsymbol f_a^{\mathrm{int}}
=
\int_{\Omega_0}
\boldsymbol P\nabla_XN_a\,dV_0.
\]
由于 Linear Tet 中 \(\boldsymbol F\)、\(\boldsymbol P\) 和 \(\nabla_XN_a\) 在单元内部都是常数，因此积分可以直接解析得到
\[\boldsymbol f = V_0\boldsymbol P\nabla_XN_a
\]
其中 \(V_0\) 为参考构形下的四面体体积。
也就是说，对于 Linear Tet，传统有限元中的体积分已经可以被压缩成非常紧凑的矩阵和向量运算。
2. Stable Neo-Hookean 模型
经典 compressible Neo-Hookean 常见形式之一是
\[\Psi = \frac{\mu}{2}(I_C-3) -\mu\ln J + \frac{\lambda}{2}(\ln J)^2
\]
Newton 代码中使用的则是代码注释中所说的
Neo-Hookean (with rest stability [Smith et al. 2018])
其能量形式可以写成
\[\Psi = \frac{\mu}{2}(I_C-3)
-\frac{\mu}{2}\ln(I_C+1)
+
\frac{\lambda}{2}(J-\alpha)^2
\]
其中 \(I_C=\boldsymbol F:\boldsymbol F,\quad J=\det\boldsymbol F\)。
与包含 \(\ln J\) 的经典 compressible Neo-Hookean 相比，这种形式避免了在 \(J\rightarrow0\) 时直接计算 \(\ln J\)，因此更适合可能出现严重压缩甚至单元翻转的图形学和机器人仿真场景。
对于超弹性材料，第一 Piola-Kirchhoff 应力由应变能对变形梯度求导得到：
\[\boxed{
\boldsymbol P
=
\frac{\partial\Psi}{\partial\boldsymbol F}
}
\]
首先考虑与 \(I_C\) 有关的部分：
\[\Psi_\mu = \frac{\mu}{2}(I_C-3) -\frac{\mu}{2}\ln(I_C+1)
\]
由于 \(I_C=\boldsymbol F:\boldsymbol F\)，对 \(\boldsymbol F\) 求导可以得到
\[\boxed{ \boldsymbol P_\mu = \mu\boldsymbol F \left( 1-\frac{1}{I_C+1} \right) }
\]
Newton 还在这一项中加入了与 \(\dot{\boldsymbol F}\) 有关的阻尼，代码形式类似于
P = F * k_mu * (1.0 - 1.0 / (Ic + 1.0)) \
+ dFdt * k_damp
它的作用是耗散能量、抑制柔性体振荡。
接下来考虑体积部分：
\[\Psi_{\mathrm{vol}}
=
\frac{\lambda}{2}(J-\alpha)^2.
\]
它对 \(\boldsymbol F\) 的贡献为
\[\boldsymbol P_{\mathrm{vol}}
=
\frac{\partial\Psi_{\mathrm{vol}}}{\partial\boldsymbol F}
=
\lambda(J-\alpha)
\frac{\partial J}{\partial\boldsymbol F}.
\]
利用
\[\frac{\partial J}{\partial\boldsymbol F} = \operatorname{cof}(\boldsymbol F) = J\boldsymbol F^{-T}
\]
得到
\[\boxed{ \boldsymbol P_{\mathrm{vol}} = \lambda(J-\alpha) \operatorname{cof}(\boldsymbol F) }
\]
代码中首先计算
J = wp.determinant(F)
f_volume = (J - alpha + act) * k_lambda
这里的 f_volume 并不是最终的节点力，也不是完整的 \(\boldsymbol P_{\mathrm{vol}}\)，而是体积能量关于 \(J\) 的导数 \(\texttt{f\_volume} = \frac{\partial\Psi_{\mathrm{vol}}}{\partial J} = \lambda(J-\alpha+\mathrm{act})\)。
代码中还有一个比较特殊的参数：
alpha = 1.0 + k_mu / k_lambda \
- k_mu / (4.0 * k_lambda)
即
\[\boxed{ \alpha = 1+\frac{3\mu}{4\lambda} }
\]
为什么这里不是简单地取 \(\alpha=1\)？
我们考虑未变形状态
\[\boldsymbol F=\boldsymbol I,
\qquad
I_C=3,
\qquad
J=1.
\]
此时 \(I_C\) 相关部分产生
\[\boldsymbol P_\mu = \mu\boldsymbol I \left( 1-\frac14 \right) = \frac34\mu\boldsymbol I
\]
也就是说，如果只考虑这一部分，参考构形并不是零应力状态。
体积部分在 \(J=1\) 时产生
\[\boldsymbol P_{\mathrm{vol}} = \lambda(1-\alpha)\boldsymbol I
\]
为了保证参考构形满足 \(\boldsymbol P(\boldsymbol I)=0\)， 需要 \(\lambda(1-\alpha) = -\frac34\mu\)，因此
\[\boxed{ \alpha = 1+\frac{3\mu}{4\lambda} }
\]
所以 \(\alpha\) 并不是一个任意的经验参数，而是一个 rest-state correction：它使完整材料模型在参考构形 \(\boldsymbol F=\boldsymbol I\) 时保持零应力状态。
3. 从能量得到 Tet 节点力
有了 \(\boldsymbol F\) 和材料模型之后，接下来就是把材料能量转化为节点力。
传统有限元通常写成
\[\boldsymbol f_a^{\mathrm{int}} = \int_{\Omega_0} \boldsymbol P\nabla_XN_a\,dV_0,
\]
或者在矩阵形式下写成
\[\boldsymbol f_{\mathrm{int}} = \int_{\Omega_0} \boldsymbol B^T\boldsymbol P\,dV_0.
\]
Newton 利用了 Linear Tet 的特点，把这一过程直接展开成对节点坐标的运算。
以体积能量为例：
\[\Psi_{\mathrm{vol}} = \frac{\lambda}{2}(J-\alpha)^2.
\]
根据链式法则，
\[\frac{\partial\Psi_{\mathrm{vol}}} {\partial\boldsymbol x_a} = \frac{\partial\Psi_{\mathrm{vol}}}{\partial J} \frac{\partial J}{\partial\boldsymbol x_a}.
\]
因此体积项对应的节点力与
\[\boxed{ \lambda(J-\alpha) \frac{\partial J} {\partial\boldsymbol x_a} }
\]
直接相关。
这里可以看到，前面代码中的
f_volume = (J - alpha + act) * k_lambda
只计算了 \(\frac{\partial\Psi_{\mathrm{vol}}}{\partial J}\)，而 \(\boldsymbol F\) 所包含的几何信息并没有消失，它继续通过 \(\frac{\partial J}{\partial\boldsymbol x_a}\) 进入最终的节点力。
因此，Newton 并没有改变有限元的基本原理，而是利用 Linear Tet 的结构，将传统的
\[\text{形函数}
\rightarrow
\text{Gauss 积分}
\rightarrow
\text{应力}
\rightarrow
\text{节点内力}
\]
压缩成了针对 Linear Tet 的直接矩阵和向量运算。
从整体上看，一个 Tet 的计算过程可以概括为
\[\boxed{ \begin{aligned} \{\boldsymbol x_a\} &\rightarrow \boldsymbol F \\ &\rightarrow I_C,\ J \\ &\rightarrow \boldsymbol P_\mu,\ \frac{\partial\Psi_{\mathrm{vol}}}{\partial J} \\ &\rightarrow \boldsymbol f_a^{\mathrm{elastic}} + \boldsymbol f_a^{\mathrm{damp}} \\ &\rightarrow \boldsymbol f_a. \end{aligned} }
\]
Newton 只显式计算 Tet 中三个节点的内力，第四个节点的内力利用单元内部力的平衡直接得到
f0 = -(f1 + f2 + f3)
为什么可以这样做呢？这来自形函数的 partition of unity：
\[\sum_{a=0}^3N_a=1,
\]
前面已经有 Linear Tet 的节点内力：
\[\boldsymbol f_a^{\mathrm{int}} = V_0\boldsymbol P\nabla_XN_a
\]
由 partition of unity 可得
\[\sum_{a=0}^3\nabla_XN_a=0.
\]
于是单元内部力满足
\[\boxed{
\sum_{a=0}^3\boldsymbol f_a=0
}.
\]
4. Newton 与传统有限元的区别
从连续介质力学角度看，Newton 并没有使用另一套有限元理论。真正不同的是软件实现路径。
传统通用隐式有限元通常采用
\[\boxed{
\boldsymbol u
\rightarrow
\boldsymbol F
\rightarrow
\text{stress}
\rightarrow
\text{material tangent}
\rightarrow
\boldsymbol f_{\mathrm{int}}
\rightarrow
\boldsymbol K_T
}
\]
这样的架构。这是因为 Abaqus、ANSYS 一类通用有限元软件需要支持大量材料模型、单元类型以及隐式 Newton-Raphson 求解，因此通常显式组织材料积分点、应力更新、材料 Jacobian 和单元切线刚度。
Newton 的 eval_tetrahedra() 则只是一个力计算 kernel，它实现的是\(\boldsymbol f(\boldsymbol x,\boldsymbol v)\)。其核心路径更加直接：
\[\boxed{
{\boldsymbol x_a}
\rightarrow
\boldsymbol F
\rightarrow
\Psi(\boldsymbol F)
\rightarrow
\boldsymbol P
\rightarrow
{\boldsymbol f_a}
}
\]
这种直接从节点位置计算单元力的形式计算路径短、数据局部性强，也非常适合 GPU 上以 Tet 为粒度进行并行计算。它与 Abaqus、ANSYS 等通用有限元软件完整的“材料积分点—材料 Jacobian—单元切线—全局 Newton”架构相比简单很多，但背后的连续介质力学基础仍然是有限变形超弹性和有限元离散。
附录
什么是边矩阵
在图形学和物理仿真的四面体实现中，经常会看到所谓的边矩阵（edge matrix）。它其实非常简单：选择四面体的一个节点作为基点，把从这个节点出发的三条边向量按列排列起来。
例如，对于四节点四面体，参考构形中的节点为
\[\boldsymbol X_0,\quad
\boldsymbol X_1,\quad
\boldsymbol X_2,\quad
\boldsymbol X_3.
\]
选择节点 \(0\) 作为基点，可以得到三条边：
\[\boldsymbol A_1=\boldsymbol X_1-\boldsymbol X_0,\\
\boldsymbol A_2=\boldsymbol X_2-\boldsymbol X_0,\\
\boldsymbol A_3=\boldsymbol X_3-\boldsymbol X_0.
\]
把这三个向量作为矩阵的三列：
\[\begin{bmatrix}
\boldsymbol X_1-\boldsymbol X_0 &
\boldsymbol X_2-\boldsymbol X_0 &
\boldsymbol X_3-\boldsymbol X_0
\end{bmatrix}
\]
这就是参考构形的边矩阵 \(\boldsymbol D_m^0\)。
变形之后，四个节点移动到
\[\boldsymbol x_0,\quad
\boldsymbol x_1,\quad
\boldsymbol x_2,\quad
\boldsymbol x_3.
\]
同样以节点 \(0\) 为基点，可以构造当前构形的边矩阵
\[\begin{bmatrix}
\boldsymbol x_1-\boldsymbol x_0 &
\boldsymbol x_2-\boldsymbol x_0 &
\boldsymbol x_3-\boldsymbol x_0
\end{bmatrix}
\]
根据变形梯度 \(\boldsymbol F\) 的定义 \(d\boldsymbol x=\boldsymbol F\,d\boldsymbol X\)，我们可以得到
\[\boldsymbol D_s = \boldsymbol F\boldsymbol D_m^0.
\]
因此
\[\boldsymbol F = \boldsymbol D_s(\boldsymbol D_m^0)^{-1}
\]
这就是 Newton 代码中使用边矩阵计算变形梯度的来源。
浙公网安备 33010602011771号