物质点法速度回传：PIC 与 FLIP 的耗散差异

同一根无阻尼弹性杆，采用相同网格、粒子和时间步，只切换速度回传方式：两个振动周期后，FLIP 的总能量约为初始值的 0.994379，PIC 则约为 0.120633。材料没有变软，也没有加入黏性，振动为何衰减得如此不同？

关键在于：网格算完一步以后，是用网格速度替换粒子速度，还是把网格速度的增量加回粒子？两者看似只差一行，却决定了粒子原有速度信息能否越过一次往返映射。

本文先读实际能量曲线，再用一个可手算的单元拆开原因。阅读只需了解速度、动量、动能，以及形函数的权重含义；完整循环可参照物质点法基本流程。这里不比较所有 MPM 变体，也不把“耗散较小”直接等同于“算法更好”。

一、先看两个周期里发生了什么

固定参考网格弹性杆在两个周期内的归一化总能量：FLIP 保持在初始值附近，PIC 明显衰减

横轴为 $t/T$ ，即时间除以解析振动周期；纵轴为 $E/E_0$ ，即数值总能量除以初始总能量。曲线来自固定两端弹性杆实验实际推进所得的数据，不是预设的衰减函数。解析无阻尼解的归一化总能量恒为 $1$ 。

两组计算均使用 $16$ 个单元、每单元 $4$ 个粒子、CFL 数 $0.2$ ，共推进 $320$ 步至 $2T$ 。固定参数为：

量	取值
杆长 $L$ 、截面积 $A$	$1\,\mathrm{m}$ 、 $0.01\,\mathrm{m^2}$
杨氏模量 $E_Y$ 、密度 $\rho$	$10^6\,\mathrm{Pa}$ 、 $1000\,\mathrm{kg/m^3}$
波速 $c=\sqrt{E_Y/\rho}$	$31.6228\,\mathrm{m/s}$
初速度幅值 $v_0$	$0.001c\approx0.0316228\,\mathrm{m/s}$
周期 $T=2L/c$	$0.0632456\,\mathrm{s}$
实际时间步 $\Delta t$	$3.9528471\times10^{-4}\,\mathrm{s}$
初始能量 $E_0$	$0.0025\,\mathrm{J}$

初始位移、应变和应力为零，粒子初速度为 $v(X,0)=v_0\sin(\pi X/L)$ ，其中 $X$ 是沿杆向右的参考坐标。图中的总能量按粒子求和：

$E=K+U, \qquad K=\sum_p\frac12m_pv_p^2, \qquad U=\sum_p\frac{V_p\sigma_p^2}{2E_Y}.$

$p$ 表示粒子， $m_p$ 为质量， $V_p$ 为体积， $v_p$ 为速度， $\sigma_p$ 为拉正应力； $K$ 、 $U$ 分别是动能与应变能，单位均为焦耳。这里用 $E_Y$ 表示材料模量，避免与图中总能量 $E$ 混淆。

这是一维固定参考网格上的线性小应变物质点离散。粒子参考坐标、线性形函数权重、梯度、质量和体积全程不变；位移会更新，但不会据此移动映射位置。它不是随粒子当前位置重算形函数的 updated-Lagrangian MPM，因此这张图不能证明网格穿越、有限应变或不同基函数下的优劣。

在这个限定模型里，PIC 的反复覆盖带来了明显更强的耗散；FLIP 曲线则在 $1$ 附近波动，终点也不恰好等于 $1$ 。这已经提示：低耗散不是严格能量守恒。

二、回传之前，先分清节点速度的时间状态

设 $n$ 表示当前时间步， $i$ 表示节点， $w_{ip}=N_i(X_p)$ 是无量纲权重。对质量非零的节点，粒子到网格的动量映射给出

$m_i=\sum_p m_pw_{ip}, \qquad \widehat v_i^n=\frac{\sum_p m_pw_{ip}v_p^n}{m_i}.$

$\widehat v_i^n$ 是刚映射、尚未处理速度边界的节点速度，不是上一时间步留下来的节点速度。对它施加本步速度约束后，记为 $\bar v_i^n$ ；再用已处理约束的节点加速度 $a_i^n$ 更新：

$v_i^{n+1}=\bar v_i^n+\Delta t\,a_i^n. \tag{1}\label{eq:pic-flip-grid-update}$

速度单位为 $\mathrm{m/s}$ ，加速度为 $\mathrm{m/s^2}$ ， $\Delta t$ 为秒。本文中的 $v_i^{n+1}$ 特指本步网格更新后的状态；下一步通常会重新映射，并不保证它等于下一步的映射速度。

PIC（Particle-in-Cell）从更新后的网格重新构造粒子速度：

$v_{p,\mathrm{PIC}}^{n+1}=\sum_iw_{ip}v_i^{n+1}. \tag{2}\label{eq:pic-flip-pic}$

FLIP（Fluid-Implicit-Particle）保留旧粒子速度，只接收网格变化量：

$v_{p,\mathrm{FLIP}}^{n+1} =v_p^n+\sum_iw_{ip}\Delta v_i. \tag{3}\label{eq:pic-flip-flip}$

若增量的起点选为式 $\eqref{eq:pic-flip-grid-update}$ 的 $\bar v_i^n$ ，则 $\Delta v_i=v_i^{n+1}-\bar v_i^n=\Delta t\,a_i^n$ 。于是区别很直接：PIC 问“网格现在认为我该多快”，FLIP 问“这一步网格让我改变多少”。

写程序时，必须保留增量两端的状态，或单独保存加速度。不能先原地覆盖旧节点速度，再用同一数组相减求增量；同样，计算两种候选回传结果或混合结果时，不能先丢掉 $v_p^n$ 。不过，保存这些数学状态不一定要求复制整套数组：像本例直接存储加速度，就不必为差分额外保存旧网格速度。

三、边界处理也是增量定义的一部分

本例由当前应力组装内力：

$f_i^n=-\sum_pV_p\sigma_p^nN_i'(X_p), \qquad a_i^n=\frac{f_i^n}{m_i}.$

$N_i'=\mathrm{d}N_i/\mathrm{d}X$ 的单位为 $\mathrm{m^{-1}}$ ，内力 $f_i$ 的单位为牛顿。两个固定端点的映射速度和加速度均设为零，内节点则使用上式。因而实际 FLIP 分支只回传

$v_p^{n+1}=v_p^n+\Delta t\sum_iw_{ip}a_i^n.$

如果另一个程序把增量定义为“最终节点速度减去最初未约束的 $\widehat v_i^n$ ”，端点还会包含 $0-\widehat v_i^n$ 这项投影修正。那是另一种增量口径，不等于这里的受约束加速度回传。不能为了套用通用差分公式，悄悄把这项修正加进本实验，也不能据此承诺不同边界实现的代码等价。

无论选择 PIC 还是 FLIP，本实验都用同一份更新后网格速度推进位移与应变：

$\begin{aligned} u_p^{n+1}&=u_p^n+\Delta t\sum_iw_{ip}v_i^{n+1},\\ \varepsilon_p^{n+1}&=\varepsilon_p^n+\Delta t\sum_iN_i'(X_p)v_i^{n+1},\\ \sigma_p^{n+1}&=E_Y\varepsilon_p^{n+1}. \end{aligned}$

$u_p$ 的单位为米， $\varepsilon_p$ 为无量纲小应变。这里不直接用 FLIP 粒子速度推进位移；速度回传仍会通过下一步动量映射影响后续位移和应力。速度梯度为何进入应变更新，可接着阅读物质点法运动学。解析解只用于初始化与对照，没有被代入这些推进式。

四、一个单元就能看见投影耗散

暂时拿掉弹性力和边界约束。取一个孤立的一维单元，两端节点位于 $0$ 、 $h$ ， $h>0$ 是单元长度。两个粒子质量均为 $1\,\mathrm{kg}$ ，分别位于 $h/4$ 、 $3h/4$ ，速度为 $+1$ 、 $-1\,\mathrm{m/s}$ 。

线性形函数的权重矩阵以粒子为行、节点为列：

$W=\begin{pmatrix} 3/4&1/4\\ 1/4&3/4 \end{pmatrix}.$

两个节点质量均为 $1\,\mathrm{kg}$ ；映射后速度为

$v_{g,0}=\frac34(1)+\frac14(-1)=\frac12\,\mathrm{m/s}, \qquad v_{g,1}=-\frac12\,\mathrm{m/s}.$

没有力，也没有速度约束，因此网格增量为零。PIC 仍会把网格速度插回粒子，得到

$v_{1,\mathrm{PIC}}=\frac34\frac12+\frac14\left(-\frac12\right) =\frac14\,\mathrm{m/s}, \qquad v_{2,\mathrm{PIC}}=-\frac14\,\mathrm{m/s}.$

FLIP 加回的增量是零，所以两个粒子仍以 $+1$ 、 $-1\,\mathrm{m/s}$ 运动。对比动量与动能：

状态	总动量 $P$ / $\mathrm{kg\,m/s}$	动能 $K$ / $\mathrm{J}$
初始粒子	$0$	$1$
映射后的网格	$0$	$0.25$
PIC 回传后的粒子	$0$	$0.0625$
FLIP 回传后的粒子	$0$	$1$

两次平均把 PIC 速度幅值变成原来的四分之一，动能则变为十六分之一。这里既没有材料阻尼，也没有时间积分力的误差：能量损失来自映射本身。网格动能只是中间表示的诊断量，不能与粒子动能相加；表格也不表示 FLIP 在网格阶段凭空补回了能量——它始终保存着原来的粒子速度。

五、为什么动量不丢，动能却会丢

上述结论有明确前提：粒子质量为正，权重非负且满足 $\sum_iw_{ip}=1$ ，计算包含完整支撑域，映射与回传使用相容权重；这里还排除了外力和边界速度修改。

节点速度是质量加权平均。平方函数的凸性给出

$\frac12m_i\bar v_i^2 \leq\sum_p\frac12m_pw_{ip}(v_p^n)^2.$

对所有节点求和，可得网格动能不超过原粒子动能；对 PIC 的节点到粒子平均再用一次相同论证，得到

$K_{\mathrm{PIC}}\leq K_{\mathrm{grid}}\leq K_{\mathrm{particles}}^n.$

而动量是速度的一次量。无力、无约束的往返映射满足

$\sum_pm_pv_{p,\mathrm{PIC}} =\sum_im_i\bar v_i =\sum_pm_pv_p^n.$

所以“动量守恒”与“动能减少”并不矛盾。若被平均的速度恰好相同，平均也未必损失动能；不能把这个不等式说成每一步必定严格下降。

这些是无力投影步骤的结论，不是完整弹性求解器的总能量定理。力会使动能与应变能相互交换，时间离散会引入误差，约束还会交换动量。固定端点的杆不是无外界作用的孤立系统，不能直接拿上面的总动量等式要求它始终保持初始动量。理想静止支座不做功，也不意味着离散边界处理自动守能。

六、FLIP 留下来的信息不一定都应该删掉

FLIP 不反复用网格平均值覆盖粒子速度，因此也会保留网格分辨不了的速度分量。例如，同一单元内放三个等质量粒子，位置为 $h/4$ 、 $h/2$ 、 $3h/4$ ，速度为 $[1,-2,1]\,\mathrm{m/s}$ 。两个节点的映射动量都为零：

$\frac34(1)+\frac12(-2)+\frac14(1)=0, \qquad \frac14(1)+\frac12(-2)+\frac34(1)=0.$

网格只看见静止，但粒子并非静止。无力时，PIC 会清除这个分量，FLIP 则保留它。它属于该映射的零空间模式，说明网格表示不能唯一还原粒子速度。

网格未解析的信息，不应一概称为噪声：它可能是希望保留的细节，也可能是数值振荡或不受网格力有效控制的运动。判断需要看空间尺度、初值与误差演化。FLIP 的低耗散是一项特点，不是自动正确的保证；PIC 的平滑能抑制某些振荡，也可能抹掉真正的振动幅值。

常见折中是混合两种候选速度：

$v_p^{n+1}=(1-\alpha)v_{p,\mathrm{PIC}}^{n+1} +\alpha v_{p,\mathrm{FLIP}}^{n+1}, \qquad 0\leq\alpha\leq1.$

这里 $\alpha=0$ 为纯 PIC， $\alpha=1$ 为纯 FLIP；阅读其他代码时应先确认其系数约定。混合可以调节平滑程度，却没有一个通用于所有问题的最佳比例。本弹性杆界面只提供纯 PIC 和纯 FLIP，没有混合系数控件。

七、在实验里怎样比较，才不会只看能量下结论

打开弹性杆实验，按以下顺序复现开头的对照：

设单元数为 $16$ 、每单元粒子数为 $4$ 、CFL 为 $0.20$ 、求解时长为 $2.00T$ ，选择 FLIP，点击“运行验证”。核对已运行算例为 $320$ 步，再查看终态总能量比。
只把速度回传改为 PIC，重新点击“运行验证”。参数变化不会自动替换已保存结果；末端能量比分别应约为 $0.994379$ 与 $0.120633$ 。
用快照查看位移、速度和应力曲线相对解析解的偏差，并检查完整能量历史。单个终点接近 $1$ ，不能排除途中波动或相位偏移。
选定一种回传方式后再做网格收敛对照。固定 CFL 时，网格加密也会缩小时间步、增加回传次数；这不是只改变空间误差的实验。