使用离散弹性棒（DER）模型模拟毛发

离散弹性棒（Discrete Elastic Rods，以下简称DER）模型[Bergou et al. 2008, Kaldor et al. 2010]是近年来被顶级特效和动画工作室（如，Weta Digital、Disney Animation Studio等）使用的一种毛发模型。这个模型可以建模不同种类的毛发，包括直发、大卷、小卷、发束、发辫等，并可以在全部种类的毛发上达到非常真实的效果。长久以来，DER一直缺乏中文资料。本文对DER进行一些简要介绍，目的在于使得未深入研究DER的同学也可以一定程度上理解DER模型。

Weta Digital使用离散弹性棒模型在电影《阿丽塔》中模拟的短发效果

DER模型是基尔霍夫弹性棒模型的离散化版本。Bergou等人于08年重新定义了离散版本的坐标架和曲率计算方式，开发了DER模型。起初版本的DER只使用空间平行输运，难以计算力的导数，因此只被用于计算弹性条扭曲以后的准静态；Kaldor等人以及Bergou等人于10年的两篇文章中，分别独立提出了时序平行输运的方法，由此，DER开始被特效工业用于大规模计算，如模拟毛发的运动等。

同时，由于计算准确度很高，在传统工程和科学领域，DER也开始崭露头角，被用于如海底线缆的运动[Jawed et al. 2014]、黄瓜卷须的盘旋[Gerbode et al. 2012]等场景的仿真，相对于传统有限元方法，DER在保证准确度的同时，大幅简化了计算量。

由于篇幅有限，以下主要介绍计算过程和方法。更详细的阐述，可以参考Jawed等人所著的书A primer on the kinematics of discrete elastic rods，其中通过总共8章介绍了DER的理论和推导过程。以下，我们将边上的变量索引以上标形式表示，同时将顶点上的变量索引以下标形式表示。

DER建模了毛发的三种不同的行为：伸缩、弯曲、和扭曲。为了达到这样的目标，尤其是建模弯曲和扭曲力，我们需要在离散成多边形线段（Polyline）的毛发上，定义多个坐标架（从而定义弯曲和扭曲角），并且需要设计这个坐标架在多边形线段不同位置上的计算方式。

首先，我们介绍DER中最重要的一个操作：平行输运。

平行输运考虑如何在一根毛发的线段和线段之间扭曲尽量小的情况下，将一根线段上的坐标架迁移到另一根线段的问题。

比如我们有两根线段，我们首先计算他们的切向量，然后便可以通过两根线段的切向量，将一根上面定义的坐标架转移到另一根上。

具体来说，如上图所示，比如我们要将向量 传输到下一个线段，那么我们首先计算两个线段的切线 和 的叉乘，单位化以后得到副法线 ，然后我们将 和 分别与 计算叉乘，得到法线 和 ，然后我们通过点乘计算 在 和 上面的截距，我们认为这个向量传输以后，在 和 的截距和在 和 的截距分别相等，因此我们用这两个截距插值 和 ，得到传输后的向量 。

在实际运算中，平行输运主要有两个地方需要用。

1. 用来初始化一根毛发所有线段上的参考坐标架（无扭曲状态时的坐标架）：我们在毛发根部，定义一个随机的（但与第0根线段垂直）坐标架，然后通过平行输运将其迭代传输到毛发尖部的线段上，这每个线段上的坐标架，就是毛发上该线段的初始参考坐标架（初始状态的参考坐标架也被称为Bishop坐标架），这个过程，也称为空间平行输运。
2. 我们将毛发的顶点位置从 时刻更新到 时刻（其中 是时间步长）之后，也要更新参考坐标架（否则坐标架将和新的线段不垂直），这里我们将毛发每根线段在 时刻的坐标架，根据 时刻切线和 时刻切线，平行输运到 时刻，这个过程，称为时序平行输运。

在每根线段上定义了参考坐标架之后，我们便可以定义扭曲角，如下图所示：

我们把经过了扭曲的坐标架（即参考坐标架以切线为轴旋转扭曲角后的坐标架）称为材料坐标架。

在DER中，基本的变量是顶点（上图中的黑球）的空间位置，和定义在每一条边（线段）上的扭曲角。因此对于 个顶点的毛发，总自由度是 个。

DER中力与力矩的计算流程

DER的计算过程可以概括为：

1. 通过顶点和扭曲角的变化，更新一些中间变量
2. 通过这些中间变量计算伸缩、弯曲、扭曲力，以及这些力对顶点和扭曲角的导数。
3. 对受力进行积分，更新顶点速度以及扭曲角的角速度
4. 对速度进行积分，更新顶点位置以及扭曲角
5. 回到1，进入下一时间步

DER的中间变量较多，如下图所示：

下面我们分别介绍这些变量的计算，以及伸缩、弯曲、扭曲力的计算。

首先是伸缩力的计算。DER中，伸缩力实际上是每2个相邻顶点间的弹簧力，需要以下这些中间变量：

• 边向量：两个相邻顶点位置相减的向量，通常由接近根部的顶点指向远离根部的顶点。
• 边长：边向量的长度
• 切线：边向量除以边长的单位向量。

有了这些中间变量，我们便可以计算伸缩力。对于顶点 来说，加载在接近根部的顶点上的伸缩力来自第 个边和第 个边，其公式是：

其中是毛发半径， 是毛发的杨氏模量（抵抗形变能力的物理量，通常在1e9 Pa到6e9 Pa之间）， 是边长（上面加横杠的是初始边长）， 是该边的切线。

接下来，我们通过上面计算的中间变量和其他一些中间变量来计算每3个相邻顶点构成的弯曲力。这其中涉及的中间变量依计算顺序分为：

• 点长：与一个顶点相连的两个边的边长平均值，称为点长。
• 参考坐标架：如上一节所述，若是第0时间步，则通过根部初始化的随机坐标架通过空间平行输运计算而得；若不是，则从上一时间步的当前线段的参考坐标架通过时序平行输运而得。
• 材料坐标架：以切线为轴，将参考坐标架旋转扭曲角后而得。具体计算如下，其中 是扭曲角， 和 分别为第 个边的材料坐标架的法线和副法线， 和 分别是第i个边的参考坐标架的法线与副法线。

• 曲率副法线：由相邻两个边的切线经过点乘与叉乘而得，其方向即为顶点i的点副法线方向 ，大小等同于两个边夹角的半角正切值的两倍（这个标量也称为离散积分曲率，记为 ，注意和下面的材料曲率向量是不同的），方向是副法线方向。具体计算如下：

• 材料曲率向量：在DER中，材料曲率是由4个分量组成的一个向量，定义在除了根部和末尾的所有顶点上。对于顶点而言，它的材料曲率向量定义为曲率副法线和它前面的（更接近根部的） 号边的法线、副法线的点乘，以及与它之后的号边的法线、副法线的点乘，即：

以上的这些变量如图所示：

有了材料曲率向量，我们便可以计算弯曲力。在弯曲力中，顶点i的受力会来自附近的 和 顶点的状态，其公式为：

这其中， 是顶点的初始点长， 是初始时刻的 ， 这一项则为顶点的材料曲率向量对顶点位置的导数，是一个4行3列的矩阵。这一项推导较为复杂，由于篇幅所限，这里直接给出推导结果（具体推导可以参考本文作者于2019年发表的文章的补充材料第10节[Fei et al. 2019]）。

其中，以 时为例（另外情况可以将 换成 或 ），有：

另外，需要注意的是，弯曲力不仅仅对顶点的位置有影响，对扭曲角也有很大影响。对扭曲角施加的弯曲力力矩可以类似对位置的弯曲力计算，不同点是，需要把材料曲率向量的导数换成对扭曲角的导数，即：

最后，我们计算扭曲力。扭曲力需要的中间变量除了以上的部分中间变量以外，还需要：

• 参考旋度：由于我们对参考坐标架使用了时序平行输运，相邻线段间的参考坐标架实际上是独立的。为了建立他们之间的关系，我们需要计算参考坐标架之间的旋度，称为参考旋度。从毛发截面看过去，如下图所示：

在毛发扭曲比较严重的时候，参考旋度可能也比较大。为了不产生奇点，通常采用增量法计算，分为4步：a）将上一条边的参考坐标架法线通过空间平行输运到当前边；b）以切线为轴，将这个向量用上一时刻的参考旋度旋转得到旋转后的向量；c）计算这个向量和当前边经时序平行输运得到的参考坐标系法线之间的夹角；d）将这个夹角加到参考旋度上，得到当前时刻的参考旋度。

• 材料旋度：将参考旋度加上当前边扭曲角和上一个边扭曲角的差值，即得到材料旋度，记为 。

扭曲力则由当前材料旋度和初始材料旋度之间的差值决定。注意与弯曲力类似，扭曲力也会同时作用于顶点和扭曲角上，因此其计算也分成作用于顶点的力与作用于扭曲角的力矩两部分：

其中 是毛发的剪切模量（抵抗剪应力的物理量。对于毛发来说，通常是杨氏模量的约0.365倍）。其中，材料旋度的梯度可以由下式计算：

毛发的隐式积分

使用显式积分来更新毛发的位置和速度是很简单的。在对顶点和扭曲角计算了三种力和力矩之后，我们将力除以顶点的质量得到加速度，将力矩除以边的惯性矩得到角加速度。然后分别乘以时间步长后加到毛发的速度和角速度上，接下来将速度和角速度分别乘以时间步长后加到毛发的位置和扭曲角上。

然而，毛发通常来说都很硬，这导致使用显式积分需要较小的（通常在1e-6秒的量级）时间步长，由此，计算一帧耗时非常大。

因此，计算毛发的运动通常采用单步隐式积分（在图形学中也称为半隐式积分）或者全隐式积分。在全隐式积分中，每一步模拟可以看成一个非线性优化问题，然后使用牛顿法求解（具体可以参考Kaufman等人[2014]的文章）。在很多应用中，使用单步隐式积分的模拟结果几乎可以媲美全隐式积分的结果，对于对精度要求不高的应用足矣。

在单步隐式积分中，每一步计算受力的负梯度（力的梯度也称为力的雅克比矩阵） ，然后求解线性方程

其中 是要求解的、包含了毛发的速度和角速度的自由度为 向量（若有 个顶点）， 是对角项包含了毛发的质量和惯性矩（对于顶点是三个同样的质量值，对于边是一个惯性矩）的大小为 的对角阵。 是上一时刻的速度， 是包含了重力的所有受力和力矩， 是时间步长。

力的负梯度矩阵通过对上面三种力求导计算。以弯曲力为例，其负梯度矩阵中第i个顶点对第k个顶点的分块通过下式计算：

其中，第二项涉及了材料曲率向量的二阶导，后者是一个三阶张量，计算相当复杂（具体可以参考本文作者2019年文章的附录）。并且，第二项在多种情况下，可能会导致线性方程的左端矩阵不定，会带来严重的数值稳定性问题。一个非常有效的近似是忽略第二项，直接使用第一项来近似计算（即单步的高斯-牛顿算法）。

在实际使用中，在时间步长不是特别大的时候，使用近似矩阵得到的模拟结果和使用完整矩阵的结果差别很小。并且，使用近似矩阵的好处有两点：

1. 大量简化了计算复杂度，以及降低了高度并行环境中高效实现的编码难度。
2. 近似负梯度矩阵是一个对称半正定矩阵，加上对角质量阵 以后变成了一个对称正定阵，很容易求解。

因此，对于三种力，可以全部使用近似的负梯度矩阵来计算。

另外，对于对精度需要较高的使用全隐式积分或者其他需要求解非线性优化问题的积分器来说，使用近似的负梯度矩阵（其中伸缩力的负梯度矩阵因计算简单，可以使用完整版，但钳制于半正定）通常也是较优的选择：可以在迭代次数差不多的情况下，避免方程左端矩阵不定或者负定导致的牛顿法找不到下降方向的情况。

对于实时应用来说，模拟的毛发通常不会太多（几百根），剩下的毛发通常是由模拟毛发插值而来。这里我们可以参考NVIDIA之前做毛发插值的方法。我们在初始化时，对每一根插值毛发找到最近的3根模拟毛发，然后根据距离计算权重。同时我们也记录最近的1根毛发。在模拟时，我们顺着毛发找到第一个与头部的碰撞点，在碰撞点前用3根毛发插值位移得到位置，在碰撞点后用最近的一根插值位移得到位置，中间会做一些混合操作。这样的插值方案可以减少插值毛发穿透头部的现象。

以下是用上面介绍的模型对3万2千根毛发（其中512根模拟）的仿真结果。在2080Ti上可以做到50FPS。

DER的CPU代码可以在GitHub上找到开源版本：nepluno/libwetcloth

参考文献：

Bergou, Miklós, et al. "Discrete elastic rods." ACM SIGGRAPH 2008 papers. 2008. 1-12.

Bergou, Miklós, et al. "Discrete viscous threads."ACM Transactions on graphics (TOG)29.4 (2010): 1-10.

Kaldor, Jonathan M., Doug L. James, and Steve Marschner. "Efficient yarn-based cloth with adaptive contact linearization." ACM SIGGRAPH 2010 papers. 2010. 1-10.

Gerbode, Sharon J., et al. "How the cucumber tendril coils and overwinds."Science 337.6098 (2012): 1087-1091.

Jawed, Mohammad K., et al. "Coiling of elastic rods on rigid substrates."Proceedings of the National Academy of Sciences 111.41 (2014): 14663-14668.

Jawed, M. Khalid, Alyssa Novelia, and Oliver M. O'Reilly. A primer on the kinematics of discrete elastic rods. Springer International Publishing, 2018.

Fei, Yun, et al. "A multi-scale model for coupling strands with shear-dependent liquid." ACM Transactions on Graphics (TOG) 38.6 (2019): 1-20.

Kaufman, Danny M., et al. "Adaptive nonlinearity for collisions in complex rod assemblies." ACM Transactions on Graphics (TOG) 33.4 (2014): 1-12.