物质点法运动学：从速度梯度到粒子变形

一个小方块转过了一个角度，它是否就产生了应变？如果计算结果显示旋转后的方块变大了，又该怀疑材料模型，还是数值算法？

回答这些问题之前，需要分清三个层次：物体的位置变化、材料内部的相对运动，以及数值更新对运动的近似。平移会改变位置，却不改变形状；刚体旋转会改变方向，却不拉伸材料；剪切则可能同时包含形状改变与旋转成分。

上一篇形函数文章解释了粒子—网格映射中的权重和梯度。本文接着讨论：网格上已有速度之后，怎样判断粒子附近的材料会如何运动、如何变形？推导沿用物质点法基础中的运动描述，重点放在几何意义，不展开本构、应力率或完整 cpGIMP。

读者只需具备偏导数、二维向量与矩阵乘法的基础。可先看第一节的四幅图，再按需要阅读公式；速度梯度与形变实验用于交互观察，不是阅读本文的前提。

一、先看四种运动，不急着写矩阵

下面始终取$x$向右、$y$向上。把小方块看成一块材料区域，内部细线连接的是同一批材料点，不是每一步重新生成的背景网格。橙色角点用于追踪方向，前后两图采用相同尺度。

图中尺寸用于几何示意，不代表具体材料参数；除了平移，后三种图示都以方块中心为运动中心。

1. 均匀平移：大家走得一样快

所有材料点的速度都相同。方块虽然移动了，但任意两点之间的距离和方向都没有变化。

存在速度，不意味着存在速度梯度。决定变形的不是“跑得多快”，而是“相邻位置跑得是否一样”。

2. 单轴伸长：两侧逐渐分开

左侧向左、右侧向右，离中心越远，水平速度的绝对值越大。水平材料线被拉长，竖向高度不变，面积因此增大。

这里指的是给定速度场中的几何伸长，不预设泊松收缩，也不根据材料参数求侧向响应。它不等同于对某种材料施加单轴应力的完整力学试验。

3. 简单剪切：上边相对下边错动

水平速度随高度变化，上边相对下边向右错动。原来的直角改变，竖向材料线也不再保持竖直，方块变成平行四边形。

这个例子说明：面积不变，也可以有变形。只检查体积或面积，无法识别剪切。

4. 刚体旋转：方向变了，形状没变

方块绕中心逆时针旋转。橙色角点的位置变了，所有材料线的方向都在转动，但边长、夹角和面积保持不变。

正方形有旋转对称性，只看轮廓容易误判；追踪同一个角点或一条材料线，才能看清它确实在转动。

运动	位置或方向	材料内部变化	面积
均匀平移	整体位置改变	距离、夹角不变	不变
单轴伸长	两侧分开	水平长度增大	增大
简单剪切	各高度相对错动	长度或夹角改变	不变
刚体旋转	整体方向改变	距离、夹角不变	不变

二、速度梯度描述的是相邻点的速度差

1. 从一段很短的材料线出发

设当前速度场为$\boldsymbol v(\boldsymbol x,t)$，$\boldsymbol x$是当前位置，$t$是时间。相邻两点之间有一个很小的间隔$\delta\boldsymbol x$，它们的速度差近似为

$$\delta\boldsymbol v \approx \boldsymbol L\,\delta\boldsymbol x, \qquad L_{ab}=\frac{\partial v_a}{\partial x_b}. \tag{1}\label{eq:kin-local-difference}$$

$\boldsymbol L$称为速度梯度，单位为$\mathrm{s}^{-1}$。本文约定：矩阵的行对应速度分量，列对应求导方向。在二维中，

$$\boldsymbol L= \begin{pmatrix} \partial v_x/\partial x & \partial v_x/\partial y\\ \partial v_y/\partial x & \partial v_y/\partial y \end{pmatrix}.$$

例如，$L_{xy}$表示水平速度沿竖直方向的变化率。上下两点即使位于同一条竖线上，也可能有不同的水平速度，这正是简单剪切的来源。

均匀平移时$\boldsymbol v=\boldsymbol b$，其中$\boldsymbol b$是不随位置变化的速度向量，因此$\boldsymbol L=\boldsymbol 0$。这与第一幅图的观察一致。

2. MPM 怎样从节点速度得到它

用$i$表示网格节点，$p$表示物质点。给定节点速度$\boldsymbol v_i$，用形函数$N_i$插值得到局部速度场：

$$\boldsymbol v^h(\boldsymbol x)=\sum_iN_i(\boldsymbol x)\boldsymbol v_i.$$

上标$h$表示离散近似。对当前物理坐标求导，并在粒子位置$\boldsymbol x_p$处取值，就得到

$$\boldsymbol L_p =\sum_i\boldsymbol v_i\otimes\boldsymbol G_{ip}, \qquad \boldsymbol G_{ip}=\nabla_{\boldsymbol x}N_i(\boldsymbol x_p). \tag{2}\label{eq:kin-grid-gradient}$$

这里的外积按$(\boldsymbol v_i\otimes\boldsymbol G_{ip})_{ab}=v_{i,a}G_{ip,b}$定义。节点速度单位为$\mathrm{m/s}$，形函数梯度单位为$\mathrm{m}^{-1}$，两者相乘正好得到$\mathrm{s}^{-1}$。

这一步容易出现两类错误：

• 把外积次序写反，导致非对角分量转置；均匀平移未必能暴露这个错误，简单剪切更适合检查。
• 直接使用单元局部坐标下的导数，没有换算成物理坐标梯度，使结果随网格尺寸出现错误缩放。

在本文对应的实验中，节点速度直接取自已知速度场，用于隔离速度梯度计算。完整 MPM 中，节点速度还受到动量映射、节点力、边界条件及时间更新顺序的影响，不能把这里的局部检查当作完整求解器验证。

三、把速度梯度拆成变形率与自旋

速度梯度同时包含拉伸、剪切和旋转的信息。将它拆成对称与反对称部分：

$$\boldsymbol D=\frac12(\boldsymbol L+\boldsymbol L^{\mathsf T}), \qquad \boldsymbol W=\frac12(\boldsymbol L-\boldsymbol L^{\mathsf T}), \qquad \boldsymbol L=\boldsymbol D+\boldsymbol W. \tag{3}\label{eq:kin-split}$$

$\boldsymbol D$是变形率张量，常称应变率张量；$\boldsymbol W$是自旋张量。两者单位均为$\mathrm{s}^{-1}$。它们描述瞬时变化率，不是已经累积的应变和转角。

1. 为什么材料线的长度变化由 D 决定

对无穷小材料线$\mathrm d\boldsymbol x$，随材料运动有$\frac{\mathrm D}{\mathrm Dt}\mathrm d\boldsymbol x=\boldsymbol L\,\mathrm d\boldsymbol x$。这里$\mathrm D/\mathrm Dt$表示物质导数，不是张量$\boldsymbol D$。于是

$$\frac{\mathrm D}{\mathrm Dt}\left|\mathrm d\boldsymbol x\right|^2 =2\,\mathrm d\boldsymbol x^{\mathsf T}\boldsymbol D\,\mathrm d\boldsymbol x. \tag{4}\label{eq:kin-line-length}$$

反对称矩阵满足$\mathrm d\boldsymbol x^{\mathsf T}\boldsymbol W\,\mathrm d\boldsymbol x=0$，因此$\boldsymbol W$不直接改变这段材料线的长度。$\boldsymbol D$则描述不同方向上的瞬时伸缩；把不同方向的材料线放在一起，也能识别剪切造成的夹角变化。

2. 刚体旋转：L 不为零，但 D 为零

设绕原点逆时针转动的角速度为常数$\Omega$，单位为$\mathrm{rad/s}$，速度场为$v_x=-\Omega y$、$v_y=\Omega x$。此时

$$\boldsymbol L= \begin{pmatrix}0&-\Omega\\\Omega&0\end{pmatrix}, \qquad \boldsymbol D=\boldsymbol 0, \qquad \boldsymbol W=\boldsymbol L.$$

速度随位置变化，所以$\boldsymbol L$不为零；但相邻材料点之间没有相对伸缩，$\boldsymbol D$为零。不能把整个$\boldsymbol L$都当成应变率，也不能看到材料线方向变化就认定它被剪切了。

3. 简单剪切：D 与 W 同时存在

设$v_x=\dot\gamma y$、$v_y=0$，其中常数$\dot\gamma$为工程剪切速率，单位为$\mathrm{s}^{-1}$。则

$$\boldsymbol L= \begin{pmatrix}0&\dot\gamma\\0&0\end{pmatrix}, \qquad \boldsymbol D= \begin{pmatrix}0&\dot\gamma/2\\\dot\gamma/2&0\end{pmatrix}, \qquad \boldsymbol W= \begin{pmatrix}0&\dot\gamma/2\\-\dot\gamma/2&0\end{pmatrix}.$$

因此，$D_{xy}$是工程剪切速率的一半。这不是少算了一项，而是张量分量与工程剪切量的定义不同。

还要注意：自旋不等于每一条材料线的角速度，也不能直接积分成任意有限变形下的总转角。在这个剪切场里，水平材料线始终水平，竖向材料线却会倾斜；它们的方向变化并不相同。有限旋转需要结合完整运动或形变梯度分析。

四、形变梯度 F 把初始材料线带到当前构形

1. L 看“此刻怎样变”，F 看“已经变成什么样”

用$\boldsymbol X$表示材料点在参考构形中的坐标，运动写成$\boldsymbol x=\boldsymbol\chi(\boldsymbol X,t)$。形变梯度定义为

$$\boldsymbol F=\frac{\partial\boldsymbol x}{\partial\boldsymbol X}, \qquad \mathrm d\boldsymbol x=\boldsymbol F\,\mathrm d\boldsymbol X. \tag{5}\label{eq:kin-F}$$

$\boldsymbol F$无量纲，用来变换材料线的相对向量；它不包含整体平移。若要画出材料域，还需要知道它的中心位置。

沿材料运动对$\eqref{eq:kin-F}$求导，利用链式法则可得

$$\dot{\boldsymbol F}=\boldsymbol L\boldsymbol F. \tag{6}\label{eq:kin-F-rate}$$

点号表示物质时间导数。注意乘法顺序是$\boldsymbol L\boldsymbol F$，不是一般意义上可以互换的两个矩阵。

对于本文的空间均匀、时间恒定的$\boldsymbol L$，且初始$\boldsymbol F(0)=\boldsymbol I$，有精确解

$$\boldsymbol F(t)=\exp(\boldsymbol L t).$$

$\boldsymbol I$是单位矩阵，$\exp$是矩阵指数，不是把每个矩阵元素分别取指数。若$\boldsymbol L$随运动和时间变化，一般不能直接套用这个常系数公式。

2. 面积变化只描述变形的一部分

定义$J=\det\boldsymbol F$。在三维中，$J$是当前体积与参考体积之比；在本文二维示意中，它是面积比。对可逆的$\boldsymbol F$，

$$\dot J=J\,\operatorname{tr}(\boldsymbol L) =J\,\operatorname{tr}(\boldsymbol D).$$

$\operatorname{tr}$表示矩阵的迹，即对角元素之和。简单剪切和刚体旋转都可以有$J=1$，但前者改变形状，后者不改变形状。检查$J$不能代替对边长、夹角和方向的检查。

对仿射运动，若材料域的初始中心为$\boldsymbol x_{c0}$、当前中心为$\boldsymbol x_c(t)$，初始角点为$\boldsymbol q_0$，则

$$\boldsymbol q(t)=\boldsymbol x_c(t) +\boldsymbol F(t)(\boldsymbol q_0-\boldsymbol x_{c0}).$$

这也是文中示意图和实验室里“方块随运动改变”的基本做法。一般非均匀速度场中的有限大小材料域还会发生非仿射变化，不能总用中心处一个$\boldsymbol F$精确表示整个区域。

五、为什么一次 Euler 更新会把旋转的方块变大

用显式 Euler 方法离散$\eqref{eq:kin-F-rate}$，得到

$$\boldsymbol F^{n+1} =(\boldsymbol I+\Delta t\,\boldsymbol L^n)\boldsymbol F^n,$$

其中上标$n$表示时间层，$\Delta t$为步长。对于从$\boldsymbol F^0=\boldsymbol I$出发的一次刚体旋转，令$s=\Omega\Delta t$，精确结果与一次 Euler 近似分别是

$$\boldsymbol F_{\mathrm{exact}}= \begin{pmatrix}\cos s&-\sin s\\\sin s&\cos s\end{pmatrix}, \qquad \boldsymbol F_{\mathrm{Euler}}= \begin{pmatrix}1&-s\\s&1\end{pmatrix}.$$

精确旋转保持长度与面积，而 Euler 近似满足

$$\det\boldsymbol F_{\mathrm{Euler}}=1+s^2, \qquad \boldsymbol F_{\mathrm{Euler}}^{\mathsf T}\boldsymbol F_{\mathrm{Euler}} =(1+s^2)\boldsymbol I. \tag{7}\label{eq:kin-euler-dilation}$$

因此，它不只是“转得不准确”，还把所有方向的长度都乘上了$\sqrt{1+s^2}$。取$\Omega=1\ \mathrm{rad/s}$、$\Delta t=1\ \mathrm{s}$，精确转角为$1\ \mathrm{rad}$，面积比为$1$；Euler 图形的面积比却达到$2$，边长增至约$1.414$倍。

这是时间离散引入的虚假伸长，不是材料突然变得可压缩。即使节点速度与$\boldsymbol L$完全正确，也可能出现这一误差。

若把固定总时长$T$均分成$N$个 Euler 小步，在该恒定旋转场中，总面积比为

$$J_N=\left[1+\left(\frac{\Omega T}{N}\right)^2\right]^N \longrightarrow 1\qquad(N\to\infty).$$

缩小步长能改善这个例子的误差，但有限步长下仍不严格保面积。更适合旋转的更新方法需要另外讨论；不能通过强行把$J$改回$1$，就认为其他运动误差也得到解决。

另一个对照是简单剪切：此时$\boldsymbol L^2=\boldsymbol 0$，矩阵指数恰好截断为$\boldsymbol I+t\boldsymbol L$。所以这个特殊场的一次 Euler 图可以与精确图重合，不代表 Euler 对所有运动都准确。

六、在计算实验室完成三个观察任务

打开速度梯度与形变实验。默认只显示初始与当前状态；播放将 1 秒的运动放慢为约 6 秒，视野保持不变。以下任务按默认中心$(0.5,0.5)\ \mathrm m$和默认速率$1$进行；若此前修改过参数，可在“计算细节”中恢复。

任务一：先凭图判断，再看 D 与 W

分别选择均匀平移、简单剪切和刚体旋转，播放并跟踪橙色角点。先回答：位置变了吗？形状变了吗？面积变了吗？

然后展开“计算细节”。预期结果是：平移的$\boldsymbol L$为零；旋转的$\boldsymbol D$为零、$\boldsymbol W$不为零；剪切的$\boldsymbol D$与$\boldsymbol W$都不为零。展开细节会暂停播放，便于对照当前状态。

任务二：让旋转误差显形

选择刚体旋转，勾选“对比一次 Euler 近似”，把进度移到$t=1\ \mathrm s$。检查精确面积保持$100\%$，而一次 Euler 近似显示$200\%$。

再把进度移到$t=0.1\ \mathrm s$，Euler 面积比应变成$1.01$。这里改变的是从起点到当前时刻的一次大步长度，不是在固定总时长内增加积分步数；不要把进度滑块误认为求解器的时间步细化试验。

任务三：检查形函数，而不是寻找“更漂亮”的运动

展开细节，在双线性、二次 B 样条和 uGIMP 之间切换，也可把初始采样点移到边界。实验给定的是仿射速度场，使用完整支持时，三种形函数都应重构同一个$\boldsymbol L$，只剩浮点量级差异。

如果只因为切换形函数，均匀旋转就明显变成拉伸，应优先检查梯度方向、物理尺度换算和边界支持是否被截断。这个任务检查的是仿射再现能力，不能据此推断三种方法在一般 MPM 动力学问题中的误差完全相同。

七、如何把这些关系放回 MPM 程序

在实际实现里，可以沿着下面的链条检查，而不是从一张变形图直接判断程序是否正确：

1. 节点速度是否对应预期时间层与边界条件？ 本文不规定所有算法都采用同一种应力更新顺序。
2. 粒子处的$\boldsymbol L_p$是否由物理梯度计算？ 用平移检验零梯度，用简单剪切检验非对角分量。
3. 分解是否满足$\boldsymbol L=\boldsymbol D+\boldsymbol W$？ 刚体旋转应有$\boldsymbol D=0$。
4. 更新后的$\boldsymbol F$是否同时通过长度、方向与面积检查？ 不只检查行列式。
5. 画出的材料域与实际积分域是否一致？ 为展示$\boldsymbol F$画出一个旋转方块，不意味着 uGIMP 的固定轴对齐积分域已经随之更新，更不等于实现了 cpGIMP。

至此，“粒子如何变形”可以分成两个相接但不同的问题：形函数梯度把节点速度变成局部运动信息，时间更新再把这种瞬时信息累积成材料的几何变化。前一环节正确，并不自动保证后一环节正确。

应力如何响应这些变化，是本构关系与应力更新的问题；PIC 和 FLIP 如何回传速度，则会影响后续时间步的节点速度与能量表现。这些问题都建立在运动学关系之上，但不需要在同一篇文章中一次讲完。

相关基础与延伸阅读

• 物质点法基础（一）：基本原理与控制方程：参考构形、现时构形与连续体运动描述。
• 标准 MPM 的控制方程离散与计算流程：将局部运动学放回完整的粒子—网格更新流程。
• MPM 形函数对比：双线性插值、二次 B 样条与网格穿越：权重、物理梯度与边界支持。

本文四种运动、角点示意与 Euler 旋转对照是围绕上述运动学关系构造的教学例子，不是某篇论文的算例复现，也不构成对完整工程求解器的精度认证。

本文含有 AI 生成内容