首

由于需要实现诸多deformable body（可形变物体）的交互处理，包括Softbody（柔体）、布料（Cloth）、绳索（Rope）以及基于他们之上的交互工作，在经过学习之后（也许是飘了），打算自己讲学习过程中遇到的思想、问题点，解决方案逐步汇总起来写一个学习系列的文章，这个方向确实有难度，也许写的不是很正确，本人同样一百分的欢迎小伙伴批评，这对我的工作也有很大的帮助。

另外附一个Unity相关的成熟deformable body库技术探讨研究的汇总篇，正在制作ing，完成后url会附上：

简介

PBD（Position Based Dynamics）算法是一种用于模拟物体的物理行为和运动的算法。其核心思想是离散化物体的运动方程，将其分解为一系列迭代步骤，每一步都用来更新物体的位置和速度。这些迭代步骤可以模拟多种物理现象，包括弹性、碰撞、液体流动等。

PBD算法是一种算法思想，为物理模拟提出了诸多的步骤去计算，算法也可以和诸多的计算分配做结合，比如利用Unity的DOTS框架使其做好多线程计算分配，或者利用CUDA计算单元+FEM（有限元）将问题矩阵化处理以实现大规模模拟。

相比直接计算力学，其也拥有两个优点

1. 直接基于位置而非受力情况来计算位置。计算的结果虽然并不是完全真实的物理模拟，但是可信度很高。
2. 性能更好，而且性能消耗的力度可控，可以根据实际应用中对真实性的需求强度来决定计算量。
3. 实时计算

算法过程

符号说明

仿真物体包括 \(N\)个节点和 \(M\)个约束。

每一个节点\(i \in [1, …, N]\)至少包含的参数有

• 质量\(m_i\)
• 位置\(p_i\) 也有很多人的文章喜欢用\(x_i\)表示，我喜欢用p，有position之意
• 速度\(v_i\)

每一个约束\(j \in [1, …, M]\)至少包含的参数有

• 约束基数\(n_j\) ：约束函数中例子的个数
• 约束函数\(C_j\)：三个维度（x,y,z）分别都有约束\(R^{3nj} -> R\)
• 约束参数
  • 刚度（stiffness）参数：\(k_j \in [0,1]\)
  • 阻尼（damping）参数：\(d_j \in [0,1]\)
  • and more……
• 约束类型
  • 等式表示 \(C_j(x_{i1}, x_{in2}, …, x_{inj}) = 0\) 
  • 不等式表示\(C_j(x_{i1}, x_{in2}, …, x_{inj}) >= 0\) 

算法流程

初始化阶段Start：

初始化节点参数

• 质量设置：\(m_i\)
• 位置设置：\(p_i\)
• 速度设置（这个一般初始化的时候都是设置为0）

初始化约束参数

• 刚度设置｜检查：\(k_j \in [0,1]\)
• 阻尼设置｜检查：\(d_j \in [0,1]\)
• and more

另外在真正的计算机模拟实现的时候，必须、一定要注意时间微分尺度的设置，不能无脑利用引擎的deltaTime( \(\text{d}t\) ) 或者\(\Delta t\) )，必须要形成自己的时间周期，这也是为什么很多成熟的柔体库会有一个Solve节点的原因，通过这个节点去约束计算时间并且汇总计算处理，甚至还要防止死亡循环（参考：https://johnaustin.io/articles/2019/fix-your-unity-timestep）\Delta

执行阶段Update：

按照执行顺序分为如下步骤

• 处理外力：\(v_{i next} = v_i + \Delta t \frac{F}{m}\)  F表示外力（显示欧拉）
• 速度阻尼：\(dampV elocities(v_1, …, v_N)\) 
• 下一个位置预测：\(p_{i next} =p_i + \Delta t v_i\) 
• 碰撞检测，生成碰撞约束：generateCollisionConstraints
• 诺干次迭代（计算结果）：每次迭代处理一遍所有约束projectConstraints\(C_j(x_{i1}, x_{in2}, …, x_{inj}) \)   注意这里最容易造成卡顿问题
• 速度、位置更新（数据应用）：update v and p
• 碰撞节点的速度更新（数据表达）：updateVelocity\(v_1, …, v_N\)

注意，迭代处理阶段可以采用以下一些方法的变种，复杂度约后越高，效果也越好，对算力要求也越高

1. 欧拉（包含显示欧拉Explicit（Forward）Euler’s Method｜隐式欧拉Implicit（Backward）Euler’s Method｜半隐式欧拉Semi-implicit Euler’s Method）
2. 二阶泰勒展开
3. Verlet积分（这个好处就在可以知道前后时间）
4. RK4（四阶伦哥库塔）Rounge-Kutta Method | 基于四阶泰勒展开的

详细过程

参数扩展

在开展详细的计算之前，由于物体有内外动静摩擦，粘性，阻尼等等因素，甚至现实生活中还要考虑空气阻力，因此开展更加详细的计算之前需要对参数扩展。

质心、角速度等参数

• 质心：\(p_c = \frac{\sum_i x_i m_i}{\sum_i m_i}\)
  • 取巧做法：直接提前设置一个特殊的质心节点，讲所需要参数提前输入处理
• 质心速度：\(v_c = \frac{\sum_i v_i m_i}{\sum_i m_i}\)
  • 取巧做法：同上，质心的坐标和位置丢去当作一般节点计算即可
• 系统角动量：\(L = \sum_i r_i (m_i v_i)\) 其中 \(r_i = p_i – p_c\) 每一个节点位置减去质心参数
• 系统转动惯量：\(I = \sum_i r_i^2\) 用这个就行了，不需要三重积分那种方法
• 系统角速度：\(\omega = I^{-1} L\)

阻尼参数

然后设定阻尼参数\(k_j \in [0,1]\)，对于每一个节点的速度\(v_i\)，需要做如下的变换：

\[\Delta v_i = v_c + \omega r_i – v_i\]

然后求的速度经过变换后的具体值

$$v_{i next} = v_i + k \Delta v_i$$

阻尼k只是一个0～1的系数，如果k=1，那么系统可以视作没形变（刚体），足够k=0，那就纯自由运动，不存在阻尼。

注意刚度系数和阻尼系数区别

刚度系数（Stiffness Coefficient）：

1. 刚度系数衡量了物体对外部力量的抵抗程度，即物体在受到外力时变形的难易程度。
2. 较高的刚度系数意味着物体更难变形，更坚硬。例如，弹簧的刚度系数越高，弹簧受力时变形越小。
3. 刚度系数通常与弹性材料的特性有关，它决定了物体回弹的速度和幅度。

阻尼系数（Damping Coefficient）：

1. 阻尼系数描述了物体在振动或运动过程中损失动能的速度，即减小振幅的速度。
2. 阻尼系数的存在可以减少振动的持续时间，使振动逐渐减小直至停止。
3. 较高的阻尼系数通常会导致振动衰减得更快，而较低的阻尼系数可能导致振动持续时间较长。

碰撞运动

碰撞运动包含碰撞检测、碰撞处理两个部分

碰撞检测

碰撞有三种状态，未碰撞、表面碰撞、穿透碰撞（卡到另一个物体的内部去了），需要处理的有两种，分别有如下约束方法（相对最简单的方法）

1. 表面碰撞(static collision)：\(C(p) = (p – q_s) n_s \)，刚度参数k=1
2. 穿透碰撞(continue collision)：\(C(p) = (p – q_c) n_c \)，刚度参数k=1

碰撞处理

碰撞发生后的处理方案无非也就是挤压推离开，最简单的办法就是根据每一个节点（一个小圆球）的半径做挤压，这也属于分子动力学的一种思想。

其他的处理还有比如SAT，EPA（GJK检测，不适用于可形变物体），几何反馈，深度图，八叉树空间构造等等……

迭代（逼近）

PBD的迭代本身属于牛顿迭代法思想

牛顿迭代法（Newton’s Method）是一种用于求解方程根的迭代数值方法，其核心思想是通过不断逼近函数的根来寻找方程的解。

 

1. 选择初始近似值： 首先，选择一个初始的近似解（通常记为x₀），这个值可以是方程解的任何估计值，但通常选择越接近真实解的值，迭代过程会更快收敛。
2. 计算函数值和导数值： 在当前近似解x₀处，计算方程左边的函数值f(x₀)和其导数（或梯度）值f'(x₀)。这是为了构建切线，并通过切线来估计下一个近似解。
3. 更新近似解： 使用以下公式来计算下一个近似解x₁：x₁ = x₀ – f(x₀) / f'(x₀)这个公式是通过在点(x₀, f(x₀))处使用切线来估计方程的根。新的近似解x₁通常更接近实际解。
4. 重复迭代： 重复步骤2和步骤3，不断更新近似解，直到满足一定的收敛条件，例如，达到预定的迭代次数或者近似解的变化很小。
5. 输出结果： 当满足收敛条件时，停止迭代，并将最终的近似解作为方程的根的估计值。

每一次迭代采用高斯算法的思想，独立处理每一个约束，比如求解\(Ax = b\)的问题时，一个个逐步求\(x_i\)，并且当前求得的\(x_i\)会用到后面的问题求解中，这样算法就逐步逐步迭代过去了

备注：

“当前求得的\(x_i\)会用到后面的问题求解中”，这句话的含义不是当前批次产生影响，这样会产生螺旋的，而是每一个批次Batch的结果输出会作为下一个批次Batch处理的输入，并且逐步趋向问题答案

迭代中约束的处理，需要分类两种处理

• 等式表示 \(C_j(x_{i1}, x_{in2}, …, x_{inj}) = 0\) 
• 不等式表示\(C_j(x_{i1}, x_{in2}, …, x_{inj}) >= 0\) 

等式约束

对于一个等式约束\(C(p) = C(p_1, …, p_n) = 0\)，对于当前的位置\(p\)我们需要计算一个\(\Delta p\)，使得质点在当前位置移动\(\Delta p\)后满足约束（\(p\)为0向量时，当前位置就是满足条件的）

$$C(p + \Delta p) \approx C(p) + \nabla _p C(p) \Delta p = 0$$

备注

复习一下泰勒级数，函数在a处的展开式如下：

$$f(a + h) = f(a) + f^`(a)h + o(h)$$

由于\(\Delta p\)在梯度的方向上，直接按照梯度方向设置步长方向设置系数\(\lambda\)去迭代目标值即可

$$ \Delta p = \lambda \nabla _p C(p) $$

带入上面的公式就可以得到

$$C(p + \Delta p) \approx C(p) + \nabla _p C(p) \lambda \nabla _p C(p) = 0$$

$$\lambda =  – \frac{C(p)}{\left\Vert \nabla _p C(p) \right\Vert ^2}$$

$$\Delta p = – \frac{C(p)}{\left\Vert \nabla _p C(p) \right\Vert ^2} \nabla _p C(p)$$

不等式约束

对于不等式约束，只需要在\(C(p) < 0\)时进行迭代变换即可，如果再考虑刚度、阻尼等参数，在处理约束函数的时候算上就行，另外也不要忽略了线形误差\(\Delta p(1 – k)\)

约束函数\(C(p)\)函数含义

距离约束

两个质点弹簧需要趋向恒定，需要双向处理，因此如果是距离约束的话\(C(p)\)是：

$$C(\overrightarrow{p1}, \overrightarrow{p2}) = \left\Vert \overrightarrow{p1} – \overrightarrow{p2} \right\Vert – d = 0$$

PS：好像任何弹簧-质点的可形变系统都是需要有距离约束，因此基础学科很重要呀

PPS：胡克定律在这个方面用的很重，\(F = kx\)，F表示力，K表示比例常数，x表示应变距离

弯曲约束

两个质点弹簧需要趋向恒定，需要双向处理，因此如果是距离约束的话\(C(p)\)

两个共轭顶点的Mesh发生弯曲，因此两个Mesh包含四个顶点组成一个约束，定义为p1,p2,p3,p4，其中p1,p2为两个Mesh的共轭顶点，约束函数为：

$$C_{bend}(p1, p2, p3, p4) = \arccos(\frac{(p2 – p1) \times (p3 – p1)}{\left\vert(p2 – p1) \times (p3 – p1)\right\vert} \ast \frac{(p2 – p1) \times (p3 – p1)}{\left\vert(p2 – p1) \times (p3 – p1)\right\vert}) – \psi _0 =0$$

太难了看不懂，直接截图吧

and more 约束…

还有很多约束，可以参考各种各样的论文

更新参数

直接更新参数即可

如果一个节点发生了碰撞，假设碰撞位置的法向量为 n \boldsymbol{n}n ，那么该节点速度垂直于 n \boldsymbol{n}n 的分量将因为摩擦力而减少，并且由于反弹会产生与 n \boldsymbol{n}n 同方向的速度分量。

扩展

动力学

除了PBD以外（PBD基于质点），需要补充一下粒子动力学（Particle Dynamics）和刚体动力学基本参数（Rigid Body Dynamics）的基本参数概念

粒子动力学

\(\overrightarrow{x}\)

\(\overrightarrow{v} = \frac{d \overrightarrow{x}}{dt} \)

\(\overrightarrow{a} = \frac{d \overrightarrow{v}}{dt} = \frac{d^2 \overrightarrow{x}}{dt^2} \)

\(M\)

\(p = M \overrightarrow{v} \)

\(\overrightarrow{F} = \frac{d \overrightarrow{p}}{dt} = M \overrightarrow{a}\)

刚体动力学

\(R\)

\(\overrightarrow{\omega} = \frac{\overrightarrow{v} \times \overrightarrow{r}}{\left\Vert \overrightarrow{v} \right\Vert}\)

\(\overrightarrow{\alpha} = \frac{d \overrightarrow{\omega}}{dt} = \frac{\overrightarrow{a} \times \overrightarrow{r}}{\left\Vert \overrightarrow{r}^2 \right\Vert}\)

\(I = RI_0R^T\)

\(L = I\overrightarrow{\omega}\)

\(\overrightarrow{\tau} = \frac{d \overrightarrow{L}}{dt}\)

复习一下牛顿定律

牛顿第一定律

一切物体保持其原来的状态，即保持静止的物体保持静止，匀速运动的物体保持匀速运动，除非有外力作用于它们。

牛顿第二定律

物体的加速度（a）与作用在物体上的总外力（F）之间成正比，且与物体的质量（m）成反比，表达式为：\(F = ma\)

本文参考

https://zhuanlan.zhihu.com/p/338474692

https://blog.csdn.net/dragonylee/article/details/121360710

https://johnaustin.io/articles/2019/fix-your-unity-timestep

https://www.xianlongok.site/post/d6327b48/#rigid-body-dynamics