可直接拖动画布中的粒子；连线透明度表示形函数权重。
总质量：粒子 / 网格
总动量：粒子 / 网格
最大映射误差
粒子动能：更新前 / 后
更新后网格动能
已推进时间 / 步数
P2G：\(m_i=\sum_p N_i m_p\)，\(p_i=\sum_p N_i m_pv_p\)。G2P：\(v_p^{PIC}=\sum_iN_i v_i^{n+1}\)，\(v_p^{FLIP}=v_p^n+\sum_iN_i\Delta v_i\)。
模型边界：该实验展示一维线性形函数的 P2G—网格更新—G2P 循环，但没有应力、内力和本构更新；外部加速度仅用于观察传递算法与固定边界的影响。
拖动 A–D 或用位置控件调整。节点为小方块；彩色圆点标记粒子中心，圆大小无物理含义。uGIMP 模式下彩色正方形按实际域宽绘制，选中域边框加粗；小域的中心点缩小以免遮挡。空心虚线小方块为外延支持节点。粗灰框是物理区域，不施加边界条件；粒子域可以越界，不裁剪。
总质量：粒子 / 网格 · kg
x 动量：粒子 / 网格 · kg·m/s
y 动量：粒子 / 网格 · kg·m/s
质量误差 · kg
动量误差 (x, y) · kg·m/s
内力合力 (x, y) · N
最大 |ΣN − 1|
最大 ‖Σ∇N‖ · m⁻¹
选中粒子的 x 向一维权重与梯度
选中粒子的二维支持权重与梯度
保留零权重但非零梯度的节点；它们可以没有质量却仍有内力贡献。
| 节点 (ix, iy) | 类型 | N | ∂N/∂x / m⁻¹ | ∂N/∂y / m⁻¹ |
|---|
活跃节点：质量、动量、速度与内力
活跃表示至少一个粒子有非零权重或梯度贡献。合计包含外延节点；无质量节点的速度记为 0，力仍照实显示。
| 节点 (ix, iy) | 类型 | m / kg | px / kg·m/s | py / kg·m/s | vx / m/s | vy / m/s | fx / N | fy / N |
|---|
跨单元边界的数值对比：双线性、二次 B 样条与 uGIMP
在距选中粒子最近的内部网格节点左右各取 ε = 0.001h，逐节点比较一维 N 与 dN/dx；不把线性导数跳跃画成连续曲线。二次 B 样条与正域宽 uGIMP 的一阶导数连续。对比始终包含三种基函数；uGIMP 使用上方保留的半宽比。
| 基函数 | 采样位置 | x 节点 | Nx | dNx/dx / m⁻¹ |
|---|
连续跨界扫描：固定节点的权重、梯度与内力贡献
独立的一维剖面：固定观察中间网格节点，扫描粒子沿 x 方向穿过其左右两个单元，不修改上方四个粒子的状态。三条曲线共用上方网格间距与 uGIMP 半宽比；选择 uGIMP 后可调整半宽比，扫描始终同时比较三种方法。
固定体积 V = 1000 cm³ = 0.001 m³、拉应力 σ = +10 kPa；一维单粒子内力贡献为 fᵢₚ = −Vσ dNᵢ/dx = −10 dNᵢ/dx，单位 N。该量不是完整物体的节点合力。权重曲线使用一维函数，不是上方二维张量积权重。
蓝色实线：线性；青色虚线：二次 B 样条；橙色点划线：uGIMP。竖线表示当前粒子位置。线性梯度和内力在折点处分段绘制，不连成斜坡；区间端点表示单侧极限，当前值遵循原实验的单侧梯度约定。
| 方法 | 当前位置 Nᵢ | dNᵢ/dx / m⁻¹ | fᵢₚ / N |
|---|
这是指定位置扫描，不是动力学模拟：播放只改变扫描位置，不计算加速度、应力更新或时间积分。约 8 秒扫过两个单元只是演示速度，不是物理时间。到终点、折叠面板、切换实验或离开浏览器标签时自动暂停。
张量积形函数与节点组装（应力采用拉正约定）：\[ \begin{aligned} N_i(x_p,y_p)&=N_{i_x}(x_p)N_{i_y}(y_p),\quad \nabla N_i=(N'_{i_x}N_{i_y},\,N_{i_x}N'_{i_y}),\\ m_i&=\sum_p N_i m_p,\qquad \boldsymbol p_i=\sum_p N_i m_p\boldsymbol v_p,\\ \boldsymbol f_i^{\mathrm{int}}&=-\sum_p V_p\boldsymbol\sigma_p\nabla N_i,\qquad \boldsymbol\sigma_p=\begin{pmatrix}\sigma_{xx}&\sigma_{xy}\\\sigma_{xy}&\sigma_{yy}\end{pmatrix}. \end{aligned} \]
令 \(h=1/(n-1)\)、\(r=(x-x_i)/h\)。线性基函数为 \(N_i=\max(0,1-|r|)\)；二次 B 样条为 \[ N_i(x)=\begin{cases} \frac34-r^2,&|r|<\frac12,\\ \frac12(\frac32-|r|)^2,&\frac12\le |r|\le\frac32,\\ 0,&|r|>\frac32. \end{cases} \] 梯度对物理坐标求导，包含 \(1/h\)。完整支持上 \(\sum_iN_i=1\)、\(\sum_i\nabla N_i=\boldsymbol0\)，所以质量与动量守恒，内部组装力合力为零（浮点误差除外）；局部节点内力不必为零。
uGIMP 将线性帽函数在均匀粒子域上精确取平均。这里的 \(\ell_p\) 是半宽，\(0<\ell_p/h\le 0.5\)；两轴共用该半宽，构成边长 \(2\ell_p\) 的轴对齐正方形：\[ \begin{aligned} \bar N_i(x_p)&=\frac{1}{2\ell_p}\int_{x_p-\ell_p}^{x_p+\ell_p}N_i^{\mathrm{lin}}(x)\,\mathrm dx,\\ \bar N'_i(x_p)&=\frac{N_i^{\mathrm{lin}}(x_p+\ell_p)-N_i^{\mathrm{lin}}(x_p-\ell_p)}{2\ell_p},\\ N_{ip}^{\mathrm{uGIMP}}&=\bar N_{i_x}(x_p)\bar N_{i_y}(y_p),\\ \nabla N_{ip}^{\mathrm{uGIMP}}&=(\bar N'_{i_x}\bar N_{i_y},\,\bar N_{i_x}\bar N'_{i_y}). \end{aligned} \] 当 \(\ell_p/h=0.5\) 时恰为上述二次 B 样条；当正半宽趋于零时，权重趋于线性帽函数，梯度在远离折点处趋于线性梯度。参见
Steffen 等（2008），§3.3、§4.1：该文使用全宽，本页参数使用半宽。本实验是固定、均匀、未变形粒子域的 uGIMP，不是随形变更新域的 cpGIMP。
模型边界：双线性在网格线上采用所选半开单元的单侧梯度（x 或 y = 1 时取最后单元）；二次 B 样条与 uGIMP 在边界使用一层外延节点以保留完整支持，不截断或重新归一化。uGIMP 域也不裁剪到物理区域，其平滑长度独立于内力组装所用的粒子体积，不要求域铺满或互不重叠；这是诊断性选择，不是材料一致的时间模拟。外延节点只是数学支持，不是边界条件。本实验没有外力、接触、本构更新、网格更新、G2P 或时间积分，不是完整二维 MPM 时间推进求解器；预设只设置当前粒子的给定应力与速度，不代表已求解连续体边值问题。
Van Genuchten SeTarantino Seff孔隙率变化
m = 1 - 1/n
进气尺度 1/α
e / b / Sr(100 kPa)
Van Genuchten 与 Tarantino 模型：\[ \begin{aligned} S_e &= \left[1+(\alpha s)^n\right]^{-m}, \\ S_{eff} &= \left[1+\left(\frac{a s e^b}{\rho^l g}\right)^{\frac{1}{1-\lambda}}\right]^{-\lambda}, \\ e &= \frac{n}{1-n}, \qquad b=\frac{1-\lambda}{\lambda}. \end{aligned} \]
| 吸力 s / kPa | VG Se | VG θ | SWRC Seff | SWRC Sr | 相对渗透率 kr |
|---|
理论来源与边界：孔隙率相关曲线采用 Tarantino（2009）形式，并参照 Zhan 等（2023）式（12）—（13）实现；这里只模拟湿润主曲线，未考虑干湿滞回，参数不能替代具体土样标定。
地形操作
先选择栅格，再改变该单元高程。
低高程高高程选中路径起点
坡降：\(s_{ij}=(z_i-z_j)/d_{ij}\)。D8 将流向分配给八邻域中正坡降最大的单元。
模型边界：这是用于解释 D8 规则的简化网格，没有执行填洼、平坦区处理、汇流累积和河网提取，因此内部洼地可能终止路径。
初始状态 · t = 0
当前状态 · t = 0.00 s
一次 Euler 近似
各图使用相同且播放中不变的坐标范围；橙色角点用于追踪方向。1 秒的运动放慢为约 6 秒演示，到终点停止；这不是受力求解。
计算细节：参数、速度梯度与公式（展开时暂停）
梯度误差 ‖L − A‖F / s⁻¹
精确 J = det F（面积比）
一次 Euler J = det F（面积比）
固定采样与域中心（位置单位 m，速度单位 m/s）| 量 | x | y |
|---|
仿射参考、形变梯度与角点坐标
A 来自给定场，独立于形函数重构；矩阵按行列 xx、xy、yx、yy 排列。F 无量纲，初始与演化角点一一对应。
| 角点 | 初始 (x, y) / m | 精确 (x, y) / m | Euler (x, y) / m |
|---|
采样支持节点：给定速度、权重与梯度
包含零权重但非零梯度的节点，以及区域外的完整支持；节点速度直接取给定场，不是由粒子质量或动量 P2G 得到。
| 位置 (x, y) / m | 类型 | vx / m/s | vy / m/s | N | ∂N/∂x / m⁻¹ | ∂N/∂y / m⁻¹ |
|---|
给定稳态仿射速度场，以 \(\boldsymbol c=(0.5,0.5)\,\mathrm m\) 为中心，直接赋值给所有网格节点（包括外延节点），从而把运动学与 P2G 误差隔离：\[ \begin{aligned} \boldsymbol v(\boldsymbol x)&=\boldsymbol A(\boldsymbol x-\boldsymbol c)+\boldsymbol b,\qquad \boldsymbol v_p=\sum_iN_i(\boldsymbol x_p)\boldsymbol v_i,\\ L_{ab}&=\sum_i v_{i,a}\,\partial_bN_i,\qquad \boldsymbol D=\tfrac12(\boldsymbol L+\boldsymbol L^\mathsf T),\qquad \boldsymbol W=\tfrac12(\boldsymbol L-\boldsymbol L^\mathsf T). \end{aligned} \] 完整支持与仿射再现应使 \(\boldsymbol L=\boldsymbol A\)（浮点误差除外）。双线性在网格线上沿用半开单元的单侧梯度约定；二次 B 样条与 uGIMP 保留外延支持，不截断、不归一化。
初始 \(\boldsymbol F_0=\boldsymbol I\)。播放逐帧取解析时刻 \(t=\Delta t\)；可选的 Euler 图始终是从起点到当前时刻的一次大步近似，不是许多小步的累积积分：\[ \begin{aligned} \dot{\boldsymbol F}&=\boldsymbol L\boldsymbol F,\qquad \boldsymbol F_{\mathrm{exact}}=\exp(\boldsymbol A\Delta t),\qquad \boldsymbol F_{\mathrm{Euler}}=\boldsymbol I+\Delta t\,\boldsymbol L,\qquad J=\det\boldsymbol F,\\ \boldsymbol x_{\mathrm{exact}}&=\boldsymbol c+\boldsymbol F_{\mathrm{exact}}(\boldsymbol x_p-\boldsymbol c)+\boldsymbol b\Delta t,\qquad \boldsymbol x_{\mathrm{Euler}}=\boldsymbol x_p+\Delta t\,\boldsymbol v_p,\\ \boldsymbol q_{\mathrm{new}}&=\boldsymbol x_{\mathrm{new}}+\boldsymbol F_{\mathrm{new}}(\boldsymbol q_0-\boldsymbol x_p). \end{aligned} \] 这里的中心位置公式适用于本页四种预设：平移时 \(\boldsymbol A=0\)，其余预设 \(\boldsymbol b=0\)。小方域的四个初始角点 \(\boldsymbol q_0\) 与采样点相距每轴 ±0.06 m，使用同一个 F 变换，不代表已更新 uGIMP 的形函数支持。
四种预设（令 \(r\) 为速率、\(s=r\Delta t\)）：- 均匀平移：\(\boldsymbol A=0\)，\(\boldsymbol b=(0.2r,0.1r)\,\mathrm{m/s}\)，此时 r 是无量纲倍率。L、D、W 均为零；域只平移，精确与 Euler 相同，J = 1。
- 单轴伸长：\(\boldsymbol A=\begin{pmatrix}r&0\\0&0\end{pmatrix}\)，r 单位 s⁻¹。D = L、W = 0；精确 x 向伸长比为 \(e^s\)，Euler 为 \(1+s\)，二者也是各自面积比。
- 简单剪切：\(\boldsymbol A=\begin{pmatrix}0&r\\0&0\end{pmatrix}\)，r 单位 s⁻¹。D 与 W 均非零；\(\boldsymbol A^2=0\)，所以本预设精确 F 与单步 Euler F 相同，J = 1。
- 刚体旋转：\(\boldsymbol A=\begin{pmatrix}0&-r\\r&0\end{pmatrix}\)，r 单位 s⁻¹（角速度 rad/s）。D = 0、W = L；精确 F 是转角 s 的旋转矩阵，精确 J = 1，而 Euler J = \(1+s^2\)。Euler 面积增大是时间离散误差，不是材料可压缩性。
模型边界：这是给定仿射场的局部运动学诊断，不求解动量方程或连续体边值问题，不计算本构、应力、内外力或接触。uGIMP 仅沿用固定、轴对齐域的形函数与梯度；图中变形方域用于展示 F，不反馈到积分域，没有实现 cpGIMP，也不是完整 MPM 动力学求解器。
首次激活时运行默认算例。
位移误差 / 峰值 RMS
速度误差 / 峰值 RMS
应力误差 / 峰值 RMS
下方空间曲线：蓝实线与圆点为粒子数值解，橙虚线为解析解；横轴始终为固定参考坐标 X / m，不是变形后位置。粒子之间仅连线辅助阅读，应力并非连续重构。
能量图：蓝实线 K/E₀，青虚线 U/E₀，紫粗实线总能量/E₀，橙点线解析总能量/E₀ = 1；灰竖线为选中快照。横轴 t/T。
当前快照粒子数据（数值 / 解析）
| X / m | u / µm | u 解析 / µm | v / m/s | v 解析 / m/s | σ / kPa | σ 解析 / kPa |
|---|
网格收敛对照
尚未运行网格对照。
| 单元 / 粒子 | 步数 / Δt (s) | 终态 u 误差 % | 终态 v 误差 % | 终态 σ 误差 % | 终态总能量/E₀ |
|---|
这是相同终止相位下的网格对照；固定 CFL 时空间与时间同时细化。并不保证每个相位、每个误差指标都单调下降；PIC 的反复投影耗散也随步数变化，不能只据终态能量判定精度。
均匀杆参考问题（应力拉正）：\[ \rho\,u_{tt}=E\,u_{XX},\quad \sigma=E\,u_X,\quad u(0,t)=u(L,t)=v(0,t)=v(L,t)=0,\quad u(X,0)=0,\quad v(X,0)=v_0\sin(kX). \] \[ k=\pi/L,\quad c=\sqrt{E/\rho},\quad \omega=kc,\quad T=2L/c,\quad u=\frac{v_0}{\omega}\sin(kX)\sin(\omega t),\quad v=v_0\sin(kX)\cos(\omega t),\quad \sigma=\frac{Ev_0}{c}\cos(kX)\sin(\omega t). \] 这些解析式直接满足波动方程、初值和两端约束，仅用于初值与结果对照，不作为时间推进中的节点速度、力或应力。
固定 SI 数据：L = 1 m，E = 10⁶ Pa，ρ = 1000 kg/m³，A = 0.01 m²；c = √1000 ≈ 31.6228 m/s，v₀ = 0.001c ≈ 0.0316228 m/s，T ≈ 0.0632456 s。位移幅值 v₀/ω ≈ 318.310 µm，应力幅值 Ev₀/c = 1 kPa；v₀/c = 0.001 是应变幅值，符合小应变假设。解析总能量 E₀ = ρALv₀²/4 = 0.0025 J。
固定参考网格上的显式循环（h = L / 单元数）：粒子位于各单元均匀子区间中点，\(V_p=A h/n_{pc}\)、\(m_p=\rho V_p\)；\(X_p,N_i(X_p),N'_i(X_p),V_p,m_p\) 全程不变。\[ m_i=\sum_p m_pN_i,\quad \bar v_i=\frac{\sum_p m_pv_p^nN_i}{m_i},\quad f_i=-\sum_pV_p\sigma_p^nN'_i,\quad a_i=f_i/m_i,\quad v_i^{n+1}=\bar v_i+\Delta t\,a_i. \] 每步先将两端节点的映射速度 \(\bar v_i\) 和加速度 \(a_i\) 都置零，再更新网格速度。回传与材料更新为\[ v_p^{n+1}=\begin{cases}v_p^n+\Delta t\sum_iN_i a_i & \text{FLIP},\\ \sum_iN_i v_i^{n+1} & \text{PIC},\end{cases}\qquad u_p^{n+1}=u_p^n+\Delta t\sum_iN_i v_i^{n+1}, \] \[ \varepsilon_p^{n+1}=\varepsilon_p^n+\Delta t\sum_iN'_i v_i^{n+1},\qquad \sigma_p^{n+1}=E\varepsilon_p^{n+1}. \] 初始 u、ε、σ 均为零。FLIP 只回传已施加约束的加速度增量，不把未约束的端点映射速度差混入增量。PIC 每步覆盖粒子速度，通常有更强投影耗散；FLIP 并非严格能量守恒。没有反弹、位置裁剪或额外人工阻尼。
误差使用固定、非零的解析峰值 RMS 尺度，而不是瞬时解析范数：\[ e_q=\frac{\sqrt{\frac1{n_p}\sum_p(q_p-q_{\mathrm{exact}}(X_p,t))^2}}{q_{\mathrm{amp}}/\sqrt2},\qquad q_{\mathrm{amp}}\in\{v_0/\omega,\ v_0,\ Ev_0/c\}. \] 显示为 \(100e_q\%\)，因此位移、速度或应力过零时仍有意义。粒子能量为 \(K=\sum_p m_pv_p^2/2\)、\(U=\sum_pV_p\sigma_p^2/(2E)\)，以离散初始能量 \(E_0=\sum_p m_pv_p(0)^2/2\) 归一化；本例均匀中点积分与上述解析 E₀ 一致。
方法边界：这是固定参考网格的线性小应变物质点离散，不是随更新位置重算形函数的 updated-Lagrangian MPM，不研究网格穿越、有限应变、接触或不同基函数优劣，也没有实现完整 cpGIMP 求解器。只用线性形函数，是为了隔离已知解析解与端点边界处理。步数为 ceil(求解时长 / (CFL·h/c))，实际 Δt 缩短以精确到达终点；CFL ∈ [0.05, 0.5] 是本固定均匀算例的保守选区，不是通用 MPM 稳定性保证。
橙色虚线箭头：原始解析场；蓝色实线箭头：当前数值结果，重合时叠在参考箭头上。三图使用相同坐标范围和速度比例：当前 1 m/s 对应图中 0.72 m 的长度；比例由初始预设确定，不随往返轮数、粒子布局或显示阶段变化。不足 1 像素的箭头省略，数值仍见表格。x 向右、y 向上，位置单位 m，速度单位 m/s。
PIC · 覆盖
FLIP · 增量
APIC · 仿射
原始仿射场的再现误差
| 方法 | 网格速度误差 / % | 回传粒子速度误差 / % | 回传 C 误差 / s⁻¹ |
|---|
速度误差是全部分量的欧氏范数误差除以原始解析场范数，不是相对上一轮的变化。矩阵误差是四粒子 C−A 的合并 Frobenius 范数；PIC／FLIP 没有 C，显示“—”。当前三个预设的速度归一化分母均非零。
角动量：轨道项与 APIC 仿射项必须一起看
关于单元中心 (0.5,0.5) m 的出平面角动量，逆时针为正，单位 kg·m²/s。“初始”指第一次映射之前，“本轮输入”指选定轮的输入。APIC 总量包含仿射贡献，不能拿其绝对值与点粒子方法作相同初始总量的排名。
| 方法 | 初始总量 | 本轮输入总量 | 网格总量 | 回传总量 | 回传轨道项 | 回传仿射项 |
|---|
\[ J_{\mathrm{APIC}}=\sum_p m_p\left[(x_p-0.5)v_{p,y}-(y_p-0.5)v_{p,x}+B_{p,21}-B_{p,12}\right],\qquad \boldsymbol B_p=\boldsymbol C_p\boldsymbol D_p. \] PIC／FLIP 只统计轨道项。默认旋转、s = 0.25 m、解析 C 时，APIC 初始轨道项 0.5、仿射项 1.5，总量 2；PIC／FLIP 初始总量为 0.5。本例 FLIP 无网格速度增量，不等于任意动力学更新下严格守恒。
线动量与 APIC 回传矩阵
线动量单位 kg·m/s。C 按行排列为 xx、xy、yx、yy，单位 s⁻¹；它是传递状态，不是变形梯度 F。
| 方法 | 输入 Px | 输入 Py | 网格 Px | 网格 Py | 回传 Px | 回传 Py |
|---|
本轮网格与粒子速度数据（图形的文字替代）
位置单位 m，节点质量单位 kg，速度单位 m/s。网格没有受力或边界修正，P2G 后速度同时就是 G2P 使用的速度。
| 方法 | 节点 (x,y) | 质量 | vx | vy | 解析 vx | 解析 vy |
|---|
| 方法 | 粒子 (x,y) | 输入 vx | 输入 vy | 回传 vx | 回传 vy |
|---|
传递公式、解析预设与模型边界
单元边长 1 m，双线性权重 w 满足完整支撑、常量与线性再现。每轮位置不变，节点质量从粒子汇集。本页 s ∈ [0.1,0.4] 避开线性插值的退化位置，不代表通用奇异矩阵处理。
\[ \boldsymbol r_{ip}=\boldsymbol x_i-\boldsymbol x_p,\quad m_i=\sum_p m_pw_{ip},\quad \boldsymbol D_p=\sum_iw_{ip}\boldsymbol r_{ip}\boldsymbol r_{ip}^{\mathsf T}. \] APIC 的双向传递为 \[ m_i\boldsymbol v_i=\sum_p m_pw_{ip}(\boldsymbol v_p+\boldsymbol C_p\boldsymbol r_{ip}),\qquad \boldsymbol v_p^{\mathrm{new}}=\sum_iw_{ip}\boldsymbol v_i, \] \[ \boldsymbol B_p^{\mathrm{new}}=\sum_iw_{ip}\boldsymbol v_i\boldsymbol r_{ip}^{\mathsf T},\qquad \boldsymbol C_p^{\mathrm{new}}=\boldsymbol B_p^{\mathrm{new}}\boldsymbol D_p^{-1}. \] 本固定布局 D = s(1−s)I，以 m² 为单位；内核按权重和相对位置实际组装。PIC／FLIP 的 P2G 不含 C 项。PIC 覆盖回传，FLIP 因没有网格更新而保留本轮输入速度；APIC 的 C 每轮都重新从网格构造，不重复赋值解析矩阵。
平移：v = (1,−0.5) m/s；剪切：vx = γ(y−0.5 m)、vy = 0，γ = 1 s⁻¹；旋转：vx = −ω(y−0.5 m)、vy = ω(x−0.5 m)，ω = 1 s⁻¹。解析矩阵选项初始化 C = A = ∇v；零矩阵选项用于观察第一次映射的信息损失。
模型边界：每一轮只是固定几何上的 P2G／G2P，不推进物理时间或位置，不计算内力、压力、本构、应力、接触或边界约束。重复映射用于隔离投影损失，不是旋转动画或受力模拟。仿射再现与角动量守恒不等于任意场无耗散，也不能据此给三种方法作完整动力学精度排名。现有弹性杆实验仍仅提供 PIC／FLIP。
