MPM 形函数对比：双线性插值、二次 B 样条与网格穿越

一个物质点从网格节点左侧移动到右侧，位置只变化了一点，为什么组装到节点上的内力可能突然改变？问题不一定出在应力更新，也可能来自形函数：质量映射使用形函数的值，内力组装使用形函数的梯度；前者连续，不代表后者连续。

标准 MPM 的控制方程离散与计算流程已经介绍了弱式、粒子积分和双向映射。本文不重复完整时间步，而是对比规则笛卡尔网格上的线性形函数、二维双线性形函数和二次 B 样条，回答三个问题：权重与梯度分别做什么，粒子穿越网格时什么量发生变化，以及边界附近为什么不能随意删掉支持节点。

理解本文需要基本的导数、矩阵乘法和 MPM 粒子—网格映射概念。形函数选择及积分误差的讨论主要参考 Steffen、Kirby 与 Berzins（2008）；文中的数值表和手算例子是为说明这些关系构造的，不是该论文的算例复现。读者可用二维 MPM 组装实验核对局部映射结果，但不依赖交互页面也能完成推导。

一、先分清权重、梯度和粒子积分

用$p$表示物质点，$i$表示网格节点。物质点携带位置$\boldsymbol{x}_p$、质量$m_p$、速度$\boldsymbol{v}_p$、当前体积$V_p$和柯西应力$\boldsymbol{\sigma}_p$。在当前构形上定义

$$N_{ip}=N_i(\boldsymbol{x}_p), \qquad \boldsymbol{G}_{ip}=\nabla_{\boldsymbol{x}}N_i(\boldsymbol{x}_p). \tag{1}\label{eq:sf-value-gradient}$$

$N_{ip}$无量纲，表示粒子与节点之间的权重；$\boldsymbol{G}_{ip}$的单位为$\mathrm{m}^{-1}$，表示权重对物理空间位置的变化率。本文取$x$向右、$y$向上，拉应力为正、压应力为负。

采用集中质量和点式粒子积分时，节点质量、动量和内力为

$$\begin{aligned} m_i&=\sum_p m_pN_{ip},\\ \boldsymbol{p}_i&=\sum_p m_p\boldsymbol{v}_pN_{ip},\\ \boldsymbol{f}^{\mathrm{int}}_i &=-\sum_p V_p\boldsymbol{\sigma}_p\cdot\boldsymbol{G}_{ip}. \end{aligned} \tag{2}\label{eq:sf-assembly}$$

前两式回答“粒子的信息分给哪些节点、各分多少”；第三式来自弱式中的应力项，不能用$N_{ip}$替换梯度。这里的负号与拉正约定和弱式推导配套，不应为了使箭头符合直觉而单独翻转。

从网格速度场计算粒子处速度梯度，同样使用$\boldsymbol{G}_{ip}$：

$$\boldsymbol{L}_p =\sum_i\boldsymbol{v}_i\otimes\boldsymbol{G}_{ip}, \qquad \boldsymbol{v}_i=\frac{\boldsymbol{p}_i}{m_i}\quad(m_i>0). \tag{3}\label{eq:sf-velocity-gradient}$$

其中$\otimes$表示张量积，$\boldsymbol{L}_p$的分量为$\partial v_\alpha/\partial x_\beta$。因此，梯度的变化不仅影响内力，也会通过运动学更新影响后续应变和应力。

还要区分两种近似：形函数决定网格场如何表示，粒子积分决定空间积分如何求和。把网格基函数换得更光滑，并不等于把粒子积分改成了精确积分。

二、一维线性形函数：值连续，斜率分段切换

设一维网格节点为$x_i=ih$，$h>0$为均匀网格间距。令$r=(x-x_i)/h$，节点$i$的线性帽函数为

$$N_i(x)=\max(1-|r|,0). \tag{4}\label{eq:sf-linear}$$

它的支持区间为$[x_i-h,x_i+h]$。在不包含折点的开区间内，导数为

$$\frac{\mathrm{d}N_i}{\mathrm{d}x} =\begin{cases} 1/h,&-1<r<0,\\ -1/h,&0<r<1,\\ 0,&|r|>1. \end{cases} \tag{5}\label{eq:sf-linear-gradient}$$

$N_i$连续，但在$r=-1,0,1$处没有唯一的经典导数。因此称它为$C^0$连续，而不是$C^1$连续。

在单元$[x_j,x_{j+1}]$内，引入局部坐标$\xi=(x-x_j)/h$，便得到更适合编程的形式：

$$N_j=1-\xi,\quad N_{j+1}=\xi, \qquad N'_j=-1/h,\quad N'_{j+1}=1/h. \tag{6}\label{eq:sf-linear-cell}$$

一般位置只有两个非零权重。粒子恰好位于节点时，一个权重可能为零，但按照选定单元计算的相应梯度仍非零。不能仅以“权重为零”为条件，跳过内力所需的梯度贡献。

程序必须明确折点归属。例如实验室对内部网格节点采用右侧单元，在最右端点采用最后一个单元的左侧导数。这是离散实现约定，不会使数学上的左右导数变成相等。

三、二次 B 样条：让一阶导数也连续

本文采用中心位于网格节点上的均匀二次 B 样条。对于内部节点，或保留完整外延支持的均匀网格，写成

$$N_i(x)=B_2(r), \qquad B_2(r)=\begin{cases} \dfrac{3}{4}-r^2,&|r|<\dfrac{1}{2},\\[4pt] \dfrac{1}{2}\left(\dfrac{3}{2}-|r|\right)^2, &\dfrac{1}{2}\le |r|<\dfrac{3}{2},\\[4pt] 0,&|r|\ge\dfrac{3}{2}. \end{cases} \tag{7}\label{eq:sf-quadratic}$$

求导时必须保留从无量纲坐标到物理坐标的链式法则：

$$\frac{\mathrm{d}N_i}{\mathrm{d}x} =\frac{1}{h}B'_2(r), \qquad B'_2(r)=\begin{cases} -2r,&|r|<\dfrac{1}{2},\\[4pt] -\operatorname{sgn}(r)\left(\dfrac{3}{2}-|r|\right), &\dfrac{1}{2}\le |r|<\dfrac{3}{2},\\[4pt] 0,&|r|\ge\dfrac{3}{2}. \end{cases} \tag{8}\label{eq:sf-quadratic-gradient}$$

这里$\operatorname{sgn}$为符号函数。在$r=\pm1/2$处两段的值和一阶导数相等；到$r=\pm3/2$时，值和一阶导数都降到零。二次 B 样条因此具有$C^1$连续性，但二阶导数仍分段跳变。

与线性帽函数相比，它有三个值得注意的差异：

• 支持更宽：单个基函数覆盖长度为$3h$的区间；粒子一般与三个节点发生非零权重耦合。
• 不是节点插值型基函数：在$x=x_i$处，$N_i=3/4$，相邻两个节点各占$1/8$，而不是中间节点独占全部权重。
• 分段位置不同：这里的多项式分段点位于节点两侧的半整数网格位置；不能直接复用线性函数的两节点单元模板。

第二点意味着，B 样条展开中的网格自由度不能不加区分地解释为重构场在该节点处的点值。例如要对重构速度场施加严格的位移或速度边界条件，需要与所用基函数和边界构造相容，不能仅照搬节点插值型有限元的直觉。

四、从一维到二维：张量积与梯度方向

设规则二维网格在两个方向上的间距分别为$h_x,h_y$，节点坐标为$(x_i,y_j)$。节点$(i,j)$的形函数由两个一维函数相乘：

$$N_{ij}(x,y)=N_i^x(x)N_j^y(y), \qquad \nabla N_{ij} =\begin{bmatrix} (N_i^x)'N_j^y\\ N_i^x(N_j^y)' \end{bmatrix}. \tag{9}\label{eq:sf-tensor-product}$$

注意，$x$方向求导时只对$N_i^x$求导，$y$方向求导时只对$N_j^y$求导。梯度不是把两个一维导数直接相乘。

对于一个矩形单元，定义$\xi=(x-x_0)/h_x$、$\eta=(y-y_0)/h_y$，四个角点的双线性函数为

$$\begin{aligned} N_{00}&=(1-\xi)(1-\eta),&N_{10}&=\xi(1-\eta),\\ N_{01}&=(1-\xi)\eta,&N_{11}&=\xi\eta. \end{aligned} \tag{10}\label{eq:sf-bilinear}$$

“双线性”指分别固定另一个坐标时，对当前坐标为线性；它不是任意方向上的全局线性函数。每个粒子一般耦合四个角点，梯度跨单元边界不连续。二次 B 样条的二维张量积一般耦合九个节点，梯度连续。

性质	一维线性／二维双线性	二次 B 样条／其二维张量积
基函数连续性	$C^0$	$C^1$
一阶导数	分段定义，跨相关边界跳变	连续，分段变化
一维一般位置的非零权重数	2	3
二维一般位置的非零权重数	4	9
粒子正好位于内部网格节点	权重集中于该节点	权重分布于周围节点
主要取舍	支持紧凑、计算简单	更平滑，但耦合范围更大

这些节点数描述一般位置的非零权重；特殊位置可能出现零权重。实际梯度支持还必须按具体函数和折点约定处理，不能机械地按表格数量截断。

五、哪些性质应该用于自检

在保留完整支持的均匀网格上，两种基函数都满足单位分解和空间一阶再现：

$$\sum_iN_i(\boldsymbol{x})=1, \qquad \sum_i\boldsymbol{x}_iN_i(\boldsymbol{x})=\boldsymbol{x}. \tag{11}\label{eq:sf-completeness}$$

第一式意味着常数场能够被再现；第二式意味着坐标场以及一般仿射场能够被再现。对其求导，得到

$$\sum_i\nabla N_i=\boldsymbol{0}, \qquad \sum_i\boldsymbol{x}_i\otimes\nabla N_i=\boldsymbol{I}, \tag{12}\label{eq:sf-gradient-identities}$$

其中$\boldsymbol{I}$为单位张量。线性函数在折点处应理解为所选单元内的单侧梯度关系，而不是宣称折点存在经典导数。

将式 \eqref{eq:sf-completeness} 和式 \eqref{eq:sf-gradient-identities} 代入组装公式，可以检查

$$\sum_i m_i=\sum_p m_p, \qquad \sum_i\boldsymbol{p}_i=\sum_p m_p\boldsymbol{v}_p, \qquad \sum_i\boldsymbol{f}^{\mathrm{int}}_i=\boldsymbol{0}. \tag{13}\label{eq:sf-conservation}$$

这里的求和包括完整支持节点，检查的是施加边界约束、质量阈值筛除和时间更新之前的组装结果。

还应区分两个容易混淆的结论：

• 若所有粒子速度相同，质量加权映射后，每个正质量节点都应得到相同速度。
• 若直接给网格自由度赋值$\boldsymbol{v}_i=\boldsymbol{a}+\boldsymbol{A}\boldsymbol{x}_i$，则由一阶再现性可重构该仿射场及其梯度$\boldsymbol{A}$；但普通质量加权 P2G 不保证把任意粒子采样的仿射速度场恰好映射成这些节点值。

总内力为零不代表每个节点都处于平衡，更不代表完整算法没有能量误差。错误的局部分配也可能在总和中抵消；质量与动量守恒检查之外，仍需要节点级对照、解析解和收敛性验证。

六、一个可手算的二维内力例子

取$h_x=h_y=0.25\,\mathrm{m}$，考察单元$[0,0.25]\times[0,0.25]\,\mathrm{m}^2$。一个粒子位于单元中心，其状态为

$$\begin{aligned} \boldsymbol{x}_p&=(0.125,0.125)^{\mathsf T}\,\mathrm{m},\\ m_p&=1\,\mathrm{kg},\qquad \boldsymbol{v}_p=(0.8,0.3)^{\mathsf T}\,\mathrm{m\,s^{-1}},\\ V_p&=1000\,\mathrm{cm}^3=10^{-3}\,\mathrm{m}^3,\\ \boldsymbol{\sigma}_p&= \begin{bmatrix}10&5\\5&20\end{bmatrix}\,\mathrm{kPa}. \end{aligned} \tag{14}\label{eq:sf-example-state}$$

这里直接给定真实体积，不把二维面积默认为体积。若程序采用单位厚度，应明确厚度约定；只有“面积乘厚度”才具有$\mathrm{m}^3$的单位。该粒子也不被假定为填满整个背景单元。

采用双线性函数，四个节点权重均为$1/4$。根据式 \eqref{eq:sf-tensor-product}，可以得到：

节点坐标 / m	权重	梯度 / m⁻¹	该粒子贡献的内力 / N
$(0,0)$	$1/4$	$(-2,-2)$	$(30,50)$
$(0.25,0)$	$1/4$	$(2,-2)$	$(-10,30)$
$(0,0.25)$	$1/4$	$(-2,2)$	$(10,-30)$
$(0.25,0.25)$	$1/4$	$(2,2)$	$(-30,-50)$

以左下角为例，先把应力换算为 Pa：

$$\boldsymbol{f}^{\mathrm{int}}_{00,p} =-10^{-3} \begin{bmatrix}10000&5000\\5000&20000\end{bmatrix} \begin{bmatrix}-2\\-2\end{bmatrix} =\begin{bmatrix}30\\50\end{bmatrix}\,\mathrm{N}. \tag{15}\label{eq:sf-example-force}$$

四个力向量之和为零。每个节点得到该粒子的质量贡献$0.25\,\mathrm{kg}$、动量贡献$(0.2,0.075)\,\mathrm{kg\,m\,s^{-1}}$。固定其他量而把体积加倍，内力加倍；质量不会自动加倍，因为本例把质量与体积作为独立给定输入。这是组装自检，不是描述材料真实变形过程。

在这个特殊位置，两个方向都恰好位于相邻节点中点，二次 B 样条的有效权重和梯度与上述四节点结果相同。因此，不能只用这个中心点例子来判断两种形函数的差异；还应把粒子移离中心，尤其观察接近节点和跨越相关分段点时的变化。

在实验室中复现时，选择$5\times5$网格和双线性函数，将粒子 A 设置为式 \eqref{eq:sf-example-state} 的状态，并保持其余粒子应力为零。此时节点内力可直接与表格比较；质量与动量表仍包含其他三个粒子的贡献，不应拿总表与单粒子贡献逐项混比。

该例没有求解边值问题：非零应力不是由加载或本构计算得到的，内力也未与外力平衡。它验证的是符号、单位和组装运算，不是材料响应。

七、网格穿越时，究竟什么量发生跳变

取内部节点$x_i$，令粒子位于$x_i-\varepsilon$或$x_i+\varepsilon$，其中$0<\varepsilon\ll h$。只观察同一个节点$i$的基函数：

基函数	$N_i(x_i-\varepsilon)$	$N_i(x_i+\varepsilon)$	左侧梯度	右侧梯度
线性	$1-\varepsilon/h$	$1-\varepsilon/h$	$1/h$	$-1/h$
二次 B 样条	$3/4-(\varepsilon/h)^2$	$3/4-(\varepsilon/h)^2$	$2\varepsilon/h^2$	$-2\varepsilon/h^2$

线性函数的权重在两侧几乎没有区别，但梯度跳跃幅度为$2/h$。对于固定体积和固定标量应力的一维粒子，单粒子节点内力贡献$-V_p\sigma_pN'_i(x_p)$会相应跳变。

二次 B 样条两侧的梯度都随着$\varepsilon\to0$趋于零，没有这样的有限跳跃。在半整数分段位置，式 \eqref{eq:sf-quadratic-gradient} 的左右极限也相等；不是把原有梯度跳跃简单移到了另一条线上。

但一次局部贡献跳变，不能直接等同于完整模拟的误差大小。节点内力本来是积分

$$\boldsymbol{f}^{\mathrm{int}}_i =-\int_{\Omega}\boldsymbol{\sigma}(\boldsymbol{x})\cdot \nabla N_i(\boldsymbol{x})\,\mathrm{d}V \approx-\sum_pV_p\boldsymbol{\sigma}_p\cdot\nabla N_i(\boldsymbol{x}_p), \tag{16}\label{eq:sf-quadrature}$$

其中$\Omega$为当前物体区域。实际误差还取决于粒子分布、粒子代表的积分体积、应力变化和边界。

Steffen 等人的分析强调：常规有限元可以按单元分开积分，从而让积分区域尊重基函数梯度的不连续位置；MPM 的积分点随材料移动，通常不再保持这种对应关系。更平滑的 B 样条能减轻相应的内力积分误差，但粒子积分近似仍然存在。

因此，实验室中的跨界梯度表只能展示这一误差机制的局部来源。要比较完整算法，应进一步设置均匀应力体或有解析解的动力学问题，独立改变网格间距、每单元粒子数和时间步，分别检查空间误差、时间误差与能量演化，不能只凭动画是否平滑下结论。

八、边界支持不能靠删节点解决

考虑一维物理区间$[0,1]$，取$h=0.25\,\mathrm{m}$，粒子正好位于左端点$x_p=0$。按均匀二次 B 样条计算，完整支持为

节点坐标 / m	权重	梯度 / m⁻¹
$-0.25$	$1/8$	$-2$
$0$	$3/4$	$0$
$0.25$	$1/8$	$2$

完整支持满足权重和为 1、梯度和为 0、一阶矩为 0。如果因为节点$-0.25$在物理区间外就删除它，权重和变为$7/8$，梯度和变为$2\,\mathrm{m}^{-1}$。此时质量丢失和虚假的总内力不是材料现象，而是支持被截断的结果。

把剩余权重除以$7/8$，只能修复这一位置的单位分解；其重构坐标会变成

$$\sum_i x_i\widetilde N_i =0\cdot\frac{6}{7}+0.25\cdot\frac{1}{7} =\frac{1}{28}\,\mathrm{m}\ne0. \tag{17}\label{eq:sf-renormalized-moment}$$

一阶再现已经被破坏。如果采用位置相关的归一化$\widetilde N_i=N_i/S$，其中$S=\sum_{j\in\mathcal R}N_j$、$\mathcal R$为保留节点集合，那么在该集合固定且$S>0$的区域内，正确导数还必须包含商法则：

$$\nabla\widetilde N_i =\frac{S\nabla N_i-N_i\nabla S}{S^2}. \tag{18}\label{eq:sf-renormalization-gradient}$$

正确求导可以恢复归一化函数的梯度和关系，却不能自动恢复一阶再现；只归一化权重而保留旧梯度则更加不一致。

实验室采用简单、可核对的处理：保留一层外延节点，并把它们纳入守恒总和。外延节点只补足数学支持，不代表已经施加固定边界、自由表面或接触条件。

实际求解器可以采用经过设计的边界基函数、外延自由度处理或其他相容方案。Steffen 等人原文还讨论了专门的边界 B 样条构造；本文展示的均匀外延方式并不是对该论文边界实现的复现。不同方案必须分别验证单位分解、再现性和边界条件，不能把内部节点公式直接裁剪后当作边界公式。

九、实现时保留一份小而明确的检查清单

实现形函数时，建议把“一维权重与导数”作为独立计算，再通过张量积形成二维支持。检查应覆盖内部一般位置、网格节点、半整数分段点和物理边界，而不只检查单元中心。

1. 坐标与单位：使用当前物理坐标；无量纲导数分别乘$1/h_x$、$1/h_y$，应力和体积使用相容单位。
2. 支持完整性：线性折点的单元归属明确；保留内力需要的零权重、非零梯度项；边界不静默截断支持。
3. 局部恒等式：检查非负权重、单位分解、空间一阶再现、梯度和，以及坐标—梯度张量恒等式。
4. 组装对照：先用单粒子检查节点贡献，再检查多粒子叠加；让体积和应力分别缩放，确认内力按比例变化。
5. 完整求解中的额外处理：小质量节点阈值、速度边界和接触投影会影响更新，不能把组装阶段的恒等式直接当成这些步骤之后仍成立的保证。
6. 验证层次：代数恒等式检查实现一致性，手算例子检查单位和索引，解析解与收敛性算例才用于判断完整方法的精度。

二次 B 样条主要改变网格基函数的光滑性；GIMP 通过粒子特征函数及粒子域平均构造有效权重，APIC 则改变粒子与网格之间传递的速度信息。它们作用于不同环节，不能因为都能改善某些数值现象，就把名称或公式互换。

选择形函数时，值得首先问的不是“阶数越高是否越好”，而是：当前误差来自梯度不连续、粒子积分、边界截断，还是时间更新？只有把这些问题分开，形函数的光滑性、支持范围和计算成本才有明确的比较依据。

参考资料

1. Steffen, M., Kirby, R. M., & Berzins, M. Analysis and reduction of quadrature errors in the material point method (MPM). International Journal for Numerical Methods in Engineering, 2008, 76: 922–948. DOI: 10.1002/nme.2360；作者提供的全文。第 3.2 节讨论网格基函数及边界构造，第 4–5 节讨论粒子积分和内力误差。
2. 物质点法基础（二）：标准 MPM 的控制方程离散与计算流程。用于补充本文省略的弱式推导、粒子积分与时间更新背景。
3. 计算实验室：二维 MPM 组装。用于检查本文的局部权重、梯度和节点内力；不包含二维本构更新或时间推进。

本文含有 AI 生成内容