从布料到绳索的 Mass-Spring 物理模拟

发布于 2025-10-20 · 更新于2026-05-11 11:06阅读
Compute ShaderPhysics SimulationXPBD

引言:连续介质的离散化幻觉

在计算机图形学中,我们看到的飘逸裙摆、随风舞动的发丝,本质上是计算机对连续介质力学(Continuum Mechanics) 的一次”离散化幻觉”。真实世界的布料是无数纤维交织的连续体,其行为遵循复杂的偏微分方程(PDE);但在以 60 FPS 运行的实时管线里,直接求解 PDE 是不切实际的算力黑洞。

为了在”物理真实”与”计算效率”之间取得平衡,弹簧质点系统(Mass-Spring System) 应运而生。它将连续的拓扑面/线降维为离散的”点”(承载质量与状态)与”线”(提供约束与作用力)。

本文将系统性地拆解这一体系,从最朴素的胡克定律出发,一路打通经典力学、数值积分、PBD/XPBD、碰撞与摩擦,最终落到 GPU Compute Shader 的工业级并行实现。

一、拓扑学基础与经典力学建模 (Force-Based Formulation)

构建任何物理模拟器的第一步,是建立数学模型。我们需要对几何拓扑做图论级别的拆解,并为每条边赋予明确的物理意义。

1.1 质点网络的拓扑结构 (Topology)

设网格(布料)或曲线(绳索)包含 个顶点,顶点集 ,边集为 。

在 Mass-Spring 模型中,根据节点间拓扑距离的不同,弹簧被严格划分为三个正交的物理维度:

A. 结构弹簧 (Structural Springs)

  • 拓扑定义:连接图论距离为 1 的相邻节点。对 2D 布料,即 连接 与 ;对 1D 绳索,即 连接 。
  • 物理意义:承担拉伸与压缩刚度 (Tensile/Compressive Stiffness),是维持物体宏观尺寸的核心。结构弹簧失效即对应布料被撕裂的物理现象。

B. 剪切弹簧 (Shear Springs)

  • 拓扑定义:仅存在于 2D/3D 网格中,连接对角节点,即 连接 。
  • 物理意义:抵抗面内剪切力 (In-plane Shear)。若缺少剪切弹簧,正方形网格在结构弹簧长度不变的情况下,可以坍缩成极度压扁的菱形。

C. 弯曲弹簧 (Bending Springs)

  • 拓扑定义:跨越一个节点的连接,距离为 2。即 连接 、;绳索中即 连接 。
  • 物理意义:抵抗面外弯曲 (Out-of-plane Bending)。丝绸的弯曲弹簧极弱,皮革或粗麻绳的弯曲弹簧则极强。

下面这个交互式拓扑图可以让你直观看到三类弹簧在 5×5 网格上的连接形态:

交互演示 · 三类弹簧的拓扑结构

点击按钮切换不同类型弹簧的可见性。可以观察:仅有结构弹簧时,网格在剪切方向几乎无约束;启用剪切弹簧后,菱形坍缩被禁止;弯曲弹簧则负责抵抗弯折。

进阶补充:现代实现(如 Bridson 等人的工作)往往用 二面角约束 (Dihedral Bending Constraint) 替代距离弯曲弹簧,以解耦”压缩”与”弯曲”两种行为。本文为了行文连贯,沿用经典三类弹簧的简化模型。

1.2 内力场:广义胡克定律的向量推导

考虑通过弹簧连接的两个质点 ,位置 ,相对距离向量 。设原长 、刚度系数 。

胡克定律的标量形式为 ,但在 3D 空间中,我们需要其向量形式。从能量视角出发,弹簧存储的弹性势能为:

由”力是势能的负梯度”这一变分原理:

对 L2 范数求导,可得施加在质点 上的弹力:

工程陷阱:当 (两质点几乎重合)时,公式产生除零奇异性 (Singularity)。代码实现中必须引入极小值 :float dist = max(length(deltaX), 1e-6f);。

1.3 能量耗散:阻尼力场 (Damping Force)

无阻尼弹簧系统是保守力场,动能与势能无限转换,绳索将永不停歇地抖动。我们需要引入粘性阻尼 (Viscous Damping),模拟材料内部纤维摩擦导致的能量热耗散。

设两点速度 ,相对速度 ,阻尼系数 。

关键设计:只有沿弹簧轴向的相对速度才产生内阻尼。正交方向的相对速度代表系统整体的旋转,不应被阻尼削减,否则会产生”非物理的太空阻滞感”。

记弹簧单位方向向量 ,则:

1.4 外力场建模 (External Forces)

A. 重力场 (Gravity Field)

最简单的恒力场:

B. 空气动力学阻力 (Aerodynamic Drag)

对高质量布料模拟,风的影响至关重要。空气阻力不仅与速度有关,还与布料的迎风投影面积 (Projected Area) 有关。

设三角形面片法线 、面积 、相对风速 ,则:

注意这里使用 而非简单的平方,是为了保留法向分量的符号——风从布料正面或背面吹过时,推力方向必须相反。在质点系统中,通常将面片的力均摊到其三个顶点上。

二、数值分析与积分器 (Numerical Integration)

至此我们已构建了二阶常微分方程组 (ODE):

由于函数高度非线性且相互耦合,解析解几乎不可能存在。我们必须依赖离散时间步 进行数值求积。

2.1 显式积分的死局 (The Curse of Explicit Integrators)

最朴素的前向欧拉法 (Forward Euler):

为什么实时引擎几乎从不直接使用它?

考察一个简谐振子 ,其中 是固有角频率。前向欧拉对应的更新矩阵的谱半径为:

也就是说,对任何时间步长 ,前向欧拉对无阻尼振子都是无条件发散的(系统能量在每一步严格递增)。当刚度 提高时, 增大,发散速度指数级飙升,系统在数帧内就会冲到 NaN。

为了把发散控制在视觉可接受范围内,必须极大缩小 (增加子步数),这会瞬间吃光 CPU 预算。

2.2 隐式欧拉 (Implicit Euler) 的降维打击

为了换取绝对稳定,离线渲染(Maya、Houdini)常用隐式欧拉:

注意公式右边用的是未来时刻的未知受力。它需要通过对力做一阶泰勒展开,构造一个巨大的稀疏矩阵,并使用共轭梯度法 (Conjugate Gradient, CG) 求解线性系统 。

  • 优点:无条件稳定, 可以开得很大,且自带数值阻尼。
  • 缺点:每帧求解大型稀疏方程组,计算量与内存压力都极重,移动端与多角色实时场景难以承受。

2.3 实时物理的王者:韦尔莱积分 (Verlet Integration)

韦尔莱算法源于分子动力学,是一种辛积分器 (Symplectic Integrator)——它在长时间积分下能近似保持相空间体积守恒(几何能量守恒)。

对位置进行前向、后向泰勒展开:

两式相加,奇数阶项相互抵消,得到基础 Verlet:

引入隐式速度(由前后位置差表示)和全局阻尼(模拟空气阻滞、提供数值稳定):

其中 为阻尼因子。隐式速度可由相邻两步位置反推:

优势:计算量与显式欧拉相当,但局部截断误差由 提升至 ,全局误差为 ,极大地推迟了系统的崩溃阈值。这也是 Jakobsen 在《Advanced Character Physics》(2001) 中将 Verlet 引入实时游戏物理的根本动力。

下面这个对比演示同时在两个 Canvas 上跑同一个简谐振子(初始位置 ,零初速度),只有积分器不同。拖动 滑块,可以亲眼看到前向欧拉的振幅随时间指数膨胀,而 Verlet 的轨迹始终被锁定在合理范围内:

交互演示 · Forward Euler vs Verlet 稳定性对比

Forward Euler能量发散
Verlet (Symplectic)能量守恒

同一组初始条件 (x₀=1, v₀=0),两种积分器并行运行 200 步。Δt 越大、k 越大,Forward Euler 的发散越剧烈;Verlet 始终保持有界振荡。

三、工业革命 —— PBD 与 XPBD 架构

2006 年,NVIDIA 的 Matthias Müller 提出了 Position Based Dynamics (PBD),彻底改变了实时柔体物理的格局。Unity 的 Cloth 组件、UE 的 Chaos 引擎、Havok 物理的布料模块,底层皆由此衍生。

3.1 PBD 的哲学:力是表象,约束是本质

PBD 的核心洞见是:计算弹簧力极易爆炸,那我们干脆不要力了。

我们将弹簧视为一个纯粹的几何约束 (Geometric Constraint) 函数:

算法流程被重构为三步:

  1. 外力预测 (Predict):仅用重力和 Verlet 推算出预测位置 。
  2. 约束投影 (Project Constraints):发现 之间距离不再是 ,通过 Gauss-Seidel 迭代 把两个点强行推/拉回满足 的位置。
  3. 更新速度 (Update Velocity):用修正后的最终位置减去上一帧位置,自动得到受约束影响后的新速度。

这套流程的精髓在于:位置永远不会越界,因此系统永远不会爆炸。

下面是 PBD 完整的每帧执行流:

PBD 算法流水线

下面这个交互演示是一根由 12 个粒子组成的绳索,顶端固定。鼠标可以拖动任意粒子,而滑块控制每帧 PBD 求解的迭代次数——你能直观感受到迭代次数对”刚度”的支配关系:迭代越多,绳索越接近刚体;迭代越少,绳索越像橡皮筋:

交互演示 · PBD 距离约束迭代

红色 = 固定点,绿色 = 自由粒子,橙色 = 正在拖动的粒子。注意:迭代 1 次时绳索可被严重拉伸(PBD 此时近似一个软弹簧);迭代 ≥ 10 次时,绳索近似刚体棒——这恰是传统 PBD 刚度依赖迭代次数的根本原因(也是 XPBD 要解决的问题)。

3.2 约束梯度的数学推导

给定标量约束 ,其中 是包含所有相关质点位置的超级向量。我们寻找一个极小位移 ,使得 。

一阶泰勒展开:

为了让位移在物理上合理(遵守动量守恒、按质量倒数加权),取:

其中 为质量倒数对角阵。代回展开式解出拉格朗日乘子 :

得到 PBD 通用位置修正公式:

代入距离约束 ,即可推导出常见代码中的对称推/拉形式。

将距离约束 完整代入 PBD 通用位置修正公式，推导过程和最终的对称推/拉形式如下：

  1. 计算约束函数的梯度
    首先，我们需要分别求约束函数对质点 和 的梯度。
    设两质点间的当前距离为 ，两点间的归一化方向向量为 。
  • 对 的梯度：

  • 对 的梯度（由于 在负号后面）：

  1. 计算分母部分
    将梯度代入公式的分母 中。由于 和 都是单位向量，它们的模长平方均为 1：

  1. 计算标量乘子
    将 和分母代入 的公式：

  1. 得到最终的位移修正量
    根据 ，分别代入 和 的权重与梯度，即可得到常用于代码实现中的对称形式：
    质点 1 的位移修正：

质点 2 的位移修正：

代码实现中的物理意义总结：

  • 方向： 决定了力的作用在两点的连线上。
  • 大小： 是当前的形变量（拉伸或压缩量）。
  • 权重分配： 确保了位移按照质量的倒数（）进行加权。质量越大的点（ 越小），位移修正量越小，完全符合动量守恒的物理直觉。如果某个点是固定点（质量无限大，），它的位移修正量就会直接变成 0。

3.3 传统 PBD 的致命缺陷与 XPBD 的诞生

PBD 的痛点在于它的”刚度”与真实物理参数解耦:

  • 迭代 5 次,布料像丝绸;迭代 20 次,布料像纸板。
  • 同样的参数下,网格越细,整体布料越软。

这让美术调参极其痛苦。

XPBD (Extended Position Based Dynamics, 2016) 在约束求解方程中重新引入了真实的柔顺度参数 (Compliance) ,即刚度的倒数:

XPBD 的拉格朗日乘子增量更新:

如此一来,无论帧率如何变化、迭代多少次,材料的物理刚度表现始终保持恒定。这就是当今 3A 引擎柔体模拟的底层骨架。

四、空间交互 —— 碰撞检测与摩擦力

系统不仅要自洽,还必须与世界交互。碰撞检测分为两个阶段:宽相 (Broad-phase) 与 窄相 (Narrow-phase)。

4.1 宽相加速:空间哈希与 BVH

如果布料有 10000 个顶点,场景有 100 个胶囊体,暴力遍历是 次检测——直接卡死。我们必须使用空间划分机制。

空间哈希 (Spatial Hashing):将 3D 空间划分为均匀网格,用质点坐标计算其所在体素索引,再用大素数哈希函数映射到 1D 数组:

其中 为大素数(如 73856093, 19349663, 83492791)。只有哈希值相同的质点和碰撞体才进入精确求交阶段,从而把 降到接近 。

对于动态场景中带层级结构的碰撞体,BVH (Bounding Volume Hierarchy) 是更通用的选择,但构建成本更高。

4.2 窄相碰撞与 SDF (Signed Distance Field)

对于复杂角色(如带盔甲的人物),解析法球/胶囊求交不再适用。现代物理倾向使用 SDF (有向距离场)——一个 3D 体纹理 (Volume Texture),存储空间中每个点到模型表面的最短带符号距离 :

  • :质点在模型外部;
  • :质点已穿模进入内部;
  • 表面法线可直接由 SDF 梯度获得:。

在 PBD 框架下,SDF 碰撞自然地表示为不等式约束:

其中 为质点半径(可为 0)。若违反,则沿负梯度方向把质点推到表面之外:。

4.3 库仑摩擦力模型 (Coulomb Friction)

真实布料不会在角色身上无限滑动。在求解完碰撞约束后,我们已得到碰撞法线 与位置修正 。摩擦力作用于垂直于 的切平面 (Tangent Plane) 上。

记切向位移:

引入静摩擦系数 与动摩擦系数 :

  • 若 ,完全锁死(相当于把切向位移设为零);
  • 否则按动摩擦衰减:,或更严格地按 截断。

五、极致工程化 —— Unity GPU Compute Shader 并行架构

理论再丰满,也无法掩盖 CPU 串行执行的羸弱。我们要把成千上万根弹簧扔进 GPU。

5.1 数据竞争噩梦 (The Data Race Nightmare)

GPU 的计算模型是 SIMT (Single Instruction, Multiple Threads)。当数千个线程同时处理不同弹簧时:

  • 线程 A 处理弹簧 ;
  • 线程 B 处理弹簧 。

两个线程会在同一周期内读写同一块 的内存——这就是经典的 Data Race(数据竞争),后果是画面疯狂闪烁、布料随机撕裂、且不可复现。

5.2 解决方案:图着色 / 红黑排序

为了避免内存冲突,必须将弹簧分组(称为 Pass/Color),保证同一组内任意两条弹簧不共享任何顶点。

  • 1D 绳索(最简情形):
    • 偶 Pass(Red):求解
    • 奇 Pass(Black):求解
  • 2D 布料:结构 + 剪切 + 弯曲弹簧的组合通常需要 4 ~ 8 个颜色组 进行分批 Dispatch。颜色数越少,迭代收敛越快;颜色数过多则 Dispatch 调用开销上升,需要折中。

下面这个动画演示了 1D 红黑染色的并行调度过程——红色 Pass 中的 5 条弹簧同时被求解(无任何冲突),然后切换到黑色 Pass:

交互演示 · 红黑图着色并行调度

当前阶段:就绪

同一颜色组内的弹簧两两不共享顶点,可以无锁并行求解。GPU 线程组隐式同步后再切换到下一颜色组,即可避免 Data Race。

5.3 工业级 Compute Shader 核心源码

下面是一套经过内存对齐、支持图着色 PBD 的高性能 GPU 核心(HLSL):

HLSL
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#pragma kernel Integrate
#pragma kernel SolveDistanceConstraints
#pragma kernel CollisionSDF

// 32 字节,自然对齐到 16,迎合 GPU 合并访存 (Coalesced Memory Access)
struct Particle {
    float3 position;     // 当前位置
    float  invMass;      // 质量倒数,0 表示固定点 (Pinned)
    float3 oldPosition;  // 上一帧位置 (Verlet 隐式速度依据)
    uint   isPinned;     // 冗余标记,便于分支预测
};

// 由 C# 端经过图着色预先分组后传入
struct Spring {
    uint  p1;
    uint  p2;
    float restLength;
};

RWStructuredBuffer<Particle> _Particles;
StructuredBuffer<Spring>     _Springs;

// ====== Uniforms ======
float  dt;
float3 gravity;
float  velocityRetain;  // 显式语义:保留多少速度,= 1 - damping
float  stiffness;       // PBD 松弛系数 [0,1],或 XPBD 的 alpha 编码

// -------------------------------------------------------------
// Kernel 1:外力预测 + Verlet 积分
// -------------------------------------------------------------
[numthreads(256, 1, 1)]
void Integrate(uint3 id : SV_DispatchThreadID)
{
    uint idx = id.x;
    Particle p = _Particles[idx];

    if (p.invMass == 0.0) return;  // 固定点不参与预测

    float3 velocity = p.position - p.oldPosition;
    p.oldPosition = p.position;

    // Verlet:x(t+1) = x(t) + v * (1 - damping) + a * dt^2
    p.position += velocity * velocityRetain + gravity * (dt * dt);

    _Particles[idx] = p;
}

// -------------------------------------------------------------
// Kernel 2:PBD 距离约束求解(C# 侧按颜色分批 Dispatch)
// -------------------------------------------------------------
[numthreads(256, 1, 1)]
void SolveDistanceConstraints(uint3 id : SV_DispatchThreadID)
{
    Spring s = _Springs[id.x];

    Particle p1 = _Particles[s.p1];
    Particle p2 = _Particles[s.p2];

    float3 delta = p2.position - p1.position;
    float  currentDist = length(delta);

    // 奇异性保护
    if (currentDist < 1e-6) return;

    float w1 = p1.invMass;
    float w2 = p2.invMass;
    float sumW = w1 + w2;
    if (sumW == 0.0) return;  // 两端都是固定点

    float3 dir   = delta / currentDist;
    float  error = currentDist - s.restLength;

    // PBD 修正,stiffness ∈ [0,1] 充当松弛因子
    float3 corr = dir * (error / sumW) * stiffness;

    // 由于上层经过图着色保证了无冲突,此处可直接写入,
    // 完全无需 InterlockedAdd / atomics
    p1.position += corr * w1;
    p2.position -= corr * w2;

    _Particles[s.p1] = p1;
    _Particles[s.p2] = p2;
}

// -------------------------------------------------------------
// Kernel 3:解析法球体碰撞 + 库仑摩擦
// -------------------------------------------------------------
float3 colliderCenter;
float  colliderRadius;
float  frictionCoefficient;  // [0,1],切向速度衰减比例

[numthreads(256, 1, 1)]
void CollisionSDF(uint3 id : SV_DispatchThreadID)
{
    uint idx = id.x;
    Particle p = _Particles[idx];

    if (p.invMass == 0.0) return;

    float3 delta = p.position - colliderCenter;
    float  distSq = dot(delta, delta);
    float  radiusSq = colliderRadius * colliderRadius;

    if (distSq < radiusSq)  // 发生穿透
    {
        float  dist   = sqrt(max(distSq, 1e-12));
        float3 normal = delta / dist;

        // 1) 位置挤出:沿法线推到表面外侧,+epsilon 避免反复抖动
        float3 correctedPos = colliderCenter + normal * (colliderRadius + 1e-3);

        // 2) 摩擦力:通过修正 oldPosition 来衰减切向"速度"
        float3 vPred       = correctedPos - p.oldPosition;
        float  vNormalMag  = dot(vPred, normal);
        float3 vNormal     = normal * vNormalMag;
        float3 vTangent    = vPred - vNormal;

        // 切向衰减(简化的库仑摩擦)
        vTangent *= max(0.0, 1.0 - frictionCoefficient);

        p.oldPosition = correctedPos - (vNormal + vTangent);
        p.position    = correctedPos;
    }

    _Particles[idx] = p;
}

5.4 C# 调度侧逻辑流 (Dispatcher)

  1. Awake():基于 Mesh 拓扑生成 Particle[] 与多组着色后的 Spring[],分别上传到独立的 ComputeBuffer。
  2. FixedUpdate():
    • 设置全局参数(dt、gravity、velocityRetain、stiffness);
    • Dispatch("Integrate") —— 外力预测;
    • Substep 子迭代循环(增加刚度的关键):
      • Dispatch("SolveDistanceConstraints") —— 颜色 0(GPU 自然同步);
      • Dispatch("SolveDistanceConstraints") —— 颜色 1;
      • … 直到所有颜色组都被求解;
    • Dispatch("CollisionSDF") —— 边界与摩擦处理。
  3. 渲染管线对接:将最终的粒子位置 Buffer 直接绑定到自定义 Shader,通过 SV_VertexID 在 Vertex Stage 中索引粒子位置进行顶点变换,实现真正的 Zero CPU Overhead 渲染管线(无需 readback)。

5.5 实时布料Demo

最后,我们用 Three.js 在浏览器里复刻一遍上述完整体系:Verlet 积分 + 多次 PBD 距离约束迭代 + 球体 SDF 碰撞 + 鼠标交互。鼠标可以拖拽布料的任意顶点,布料挂在两个角点,自由垂落并与红色球体发生碰撞:

交互演示 · 自由落体、碰撞与空气动力学

FPS: --

玩法：调节风力拉杆，布料会受到基于法线的真实空气动力学升力。点击【剪断/自由落体】可观察布料掉落并包裹球体的精密碰撞过程！

六、高阶视觉化 —— 法线重计算与自阴影

物理对了,但若渲染光影错了,作品依然是失败的。

布料顶点在 GPU 中被移动后,Mesh 原有的法线 (Normals) 与切线 (Tangents) 数据立刻失效。使用旧法线渲染,被风吹起的布料依然呈现”平铺地面”的光照分布,违和感极强。

我们需要再编写一个 Compute Shader Kernel RecalculateNormals:

  1. 遍历所有三角面片(Triangle ID 作为线程索引);
  2. 根据三个顶点的新位置 计算面法线:
  3. 利用原子操作(InterlockedAdd,通常需将法线分量编码为定点整数避免并发写冲突)将面法线累加到顶点法线累加器;
  4. 在第二个 Kernel 中遍历顶点,归一化得到平滑后的顶点法线。

对于高级布料/毛发渲染(如 Kajiya-Kay 或 Marschner 各向异性高光模型),光照对法线与切线的精度极度敏感。法线一旦正确,配合 PBR 光照模型,即可呈现 3A 级次世代的柔体表现。

总结:代码与物理的交响曲

从一根最朴素的弹力系数 ,到积分发散的数值黑洞,再到拯救一切的基于位置动力学 (PBD/XPBD);从 CPU 苦涩的串行循环,到 GPU 千军万马式的红黑着色调度——这数以千计的代码行,只是为了在虚拟的像素宇宙中,还原一阵微风拂过少女发丝时那一毫秒的真实。

计算机图形学的魅力正在于此:我们用最严谨冰冷的数学矩阵,绘制出最具生命感与温度的动态世界。

愿这篇深入到底层的剖析,能成为你在物理渲染道路上的一座指路灯塔。

📚 参考文献与延伸阅读

游戏物理引擎(五) Mass-Spring |知乎|小张