【文献阅读】Peer2017 SPH弹性力

1 Overview

题目：An Implicit SPH Formulation for Incompressible Linearly
Elastic Solids

发表时间：2017年

期刊/会议：CGF

作者：Andreas Peer , Christoph Gissler, Stefan Band and Matthias Teschner
University of Freiburg, Germany

摘要：我们提出了一种新的用于可变形固体的光滑粒子流体动力学（SPH）公式。我们的方法的关键方面是隐式弹性力和变形梯度的自适应SPH公式，与之前的工作相比，该公式允许直接从SPH变形梯度中提取旋转。提出的隐式概念完全基于线性公式。当使用线性应变张量时，需要对变形梯度进行旋转感知计算。与现有工作相比，各自的旋转估计完全是在SPH概念内实现的，使用了一种新的公式，并结合了一阶一致性的核梯度校正。与显式形式相比，建议的隐式公式和自适应旋转估计允许显著更大的时间步长和更高的刚度。性能增益因子高达100。可变形固体的不可压缩性用ISPH压力求解器来解释。这进一步允许基于压力的边界处理和与SPH流体和刚体相互作用的可变形体的统一处理。自碰撞由压力解算器隐式处理

图片：

图：三个方块中，中间那个是本文的模拟结果。其他两个使用了St Venant-Kirchhoff模型。本文优势在于：剪切变形没有与旋转混淆，并且时间步要比St大五倍。

2 Methods

作者认为，在以往模拟中，系统总把旋转误认为剪切变形。本文的解决方案是：使用共旋(corotated)模型，即把旋转矩阵从形变梯度中抽离出来，然后融入参考构形里面，即让参考构形进行"共旋"，再重新计算形变梯度。

（注：构形configuration的意思就是坐标系，这是连续介质力学里面的概念。由于我不知道正确的翻译，这里“构型”和“构形”混用。初始构型/参考构型就是没进行变形前的坐标系，当前构型就是变形后的坐标系）

其算法概览如下

翻译：

1. 提取旋转：
   1.1 计算形变梯度F
   1.2 从F中提取旋转分量R
2. 利用提取的旋转计算弹性力
   2.1 通过R计算共旋的形变梯度F*
   2.2 通过F*计算应变张量$\epsilon$
   2.3 通过$\epsilon$计算应力张量P
   2.4 从P计算弹性力f

2.1形变梯度F

按照上面提到的算法流程，第一步是算出F

F是一个二阶张量，代表从初始构型到当前构型的位移（即变形）。所以它是一个变形矩阵。Remember that 矩阵的几何意义是线性变换。故**变形矩阵的物理意义是这样一种线性变换：它包含所有的剪切、旋转、拉伸。**只要物体中任意点的坐标左乘这一矩阵，就能得到这一点变形后的坐标。

上式含义：当前构型对初始构型的导数。

当前构型和初始构型是两种坐标系。

当前构型指的是：变形之后的坐标系。
初始构型指的是：变形之前的坐标系。

注意：初始构型是去掉了刚体平移的。但是！没有将刚体旋转从初始构型里面去掉，因为刚体旋转和剪切很容易混在一起。所以本文的创新点就在于采用共旋模型将刚体旋转给分离出来！

假如没有提出来，SPH是怎么样计算F的呢？

就是粒子离邻居粒子距离的梯度。由于SPH公式能将对任何物理量的算子挪到核函数上，故距离的梯度变为了距离对W梯度的并乘。

注1：对矢量的梯度就是nabla算子并乘矢量（并乘标志在这里写为$\otimes$，实际上单纯什么符号都不写也是并乘），这点请看张量分析

注2：这里距离指的是当前粒子对邻居粒子的距离。

注3：由于SPH公式需要邻域搜索，变形后邻域可能发生改变，但一个通常的假设是认为变形后邻域仍然不变。

注4：这里的上标指的是初始构型。

注5：V指的是粒子体积，SPH中通常用质量除以密度代替。这一项来源于SPH公式。

尝试从物理上解释一下这个公式：
简单来说，该公式就是邻居粒子与当前粒子的距离的梯度。为什么是梯度呢？想象一个流场，里面有空气在流动，空间的每个点都标了一个小箭头，这些小箭头大小就是流场的梯度大小，方向就是梯度方向。某一粒子的运动的动力是什么？就是周围的粒子推动和挤压。离得近了，就推他。离得远了，就拉他。当然，超过一定范围就没有作用了。

原文提到，提出共旋矩阵这一想法最早是Becker09做的，只不过它们用的是MLS，本文的贡献在于仅用SPH，无需MLS。因此就需要对SPH公式和核函数进行改造。改造如下。

2.2函数改造

为了利用SPH捕捉旋转，有一个条件：即核函数的梯度一阶相容(consistent)(TODO：相容具体是什么意思？)

因此他们借用Bonet99的做法，对核函数梯度进行了修正。修正方法就是左乘一个矩阵

这个矩阵L是

2.3旋转的提取

F中是包含旋转、剪切和拉伸的。要去掉旋转，有两种做法：
（1）将已经旋转了的F再转回去
（2）在计算F前就旋转，让初始构型本身包含旋转

其实（1）和（2）的区别就在于（1）是旋转物体，而（2）是旋转坐标系。

Becker09等人采取做法（1）。本文采取做法（2）。

Becker等人公式：

Becker等人公式翻译成SPH公式

从中看到，Becker等人就是用了相对位移（相对）的梯度（原来是距离的梯度）。

本文则将两端同时乘以R，得到

其中

就是说，改造以后的F*就是去掉了旋转的了。

至于R怎么算，要去看Muller16。（我猜是polar decomposition)。

（PS：共旋模型其实是连续介质力学很常见的一个模型，这个文章应用于了SPH上，即其创新点）

接下来，就是常规步骤了，即从变形梯度计算弹性力。

2.4从变形梯度计算弹性力

先计算应变

再计算应力（PK1)

最后对应力求散度，就是弹性力。

这里把i和j分开了，是为了保证SPH公式的对称性。

2.5 隐式时间积分

以上2.1-2.4节完整推导出了弹性力的公式。接下来，我们进行隐式时间积分的处理。首先回顾一下刚刚推导的弹性力公式：

其中PK1为

应变为

形变梯度为

如果采用的是显式时间积分的话，粒子位置使用当前时刻与初始位置。即：
$$\begin{gathered} \mathbf{d}={\mathbf{x}}^{\mathbf{t}} \\ {\mathbf{d}}^0={\mathbf{x}}^0 \end{gathered}$$

于是显式的下一时刻速度为

如果是隐式时间积分的话，粒子位置使用下一时刻和初始位置，即
$$\begin{gathered} \mathbf{d}={\mathbf{x}}^{\mathbf{t}+\mathbf{\Delta }\mathbf{t}} \\ {\mathbf{d}}^0={\mathbf{x}}^0 \end{gathered}$$

于是隐式的下一时刻速度为

将下一时刻速度换为

代入隐式公式，得到

将下一时刻弹性力换位

代入隐式公式，得到

这就是最终隐式时间积分的速度更新公式。

其中$0$是0向量
这个式子就是线性方程组
$$\mathbf{A}{\mathbf{v}}^{\mathbf{t}+\mathbf{\Delta }\mathbf{t}}=\mathbf{b}$$

其中未知量是下一时刻速度，而b就是当前时刻的量

3 Implements

实际上还采用了隐式时间积分和共轭梯度法加速求解，此处就不叙述了。

此外，压力求解是采用的IISPH。

3.1 初始构形预计算

仿照【BIT09】的做法，我们不存储粒子的绝对位置，而是存储粒子与其邻居之间的距离。每次插入新粒子的时候，将执行邻域搜索，并预计算修正矩阵，即下式：

3.2 求解线性系统

上一章节中我们推导了隐式积分，最后得到了一个线性系统。本节我们来求解该线性系统。我们采用共轭梯度法。为此，分为三个步骤：
（1）计算核函数梯度，即

（2）计算下式的等号右端项（即b）

（3）迭代求解线性系统

3.2.1 弹性力函数的计算

不论是等号右端，还是等号左端，都需要计算弹性力函数。这一函数即

在计算这个式子的时候，需要如下步骤，总结为如下算法流程：

解释如下：
总共需要两次对全体粒子的迭代。
在第一次迭代中，目的是求解PK1。为此，我们先求解修正后的形变梯度F*，然后计算应变$\epsilon$，最后计算应力第一柯西皮奥拉应力PK1。在存储P的时候，只需要6个量，因为该矩阵是对称的。
在第二次迭代中，目的是求解弹性力f。

3.2.2 修正的形变梯度F*

此外，对于F的求解，我们先重写F公式

因为

所以代入后得

3.2.3 迭代求解器

我们采用共轭梯度法求解线性系统。其中必然要计算矩阵向量乘积$\mathbf{A}\mathbf{p}$
但是我们不显式地构建矩阵A，而是一行行地构建$\mathbf{A}\mathbf{p}$

3.2.4 计算中可跳过形变梯度

其中关键在于计算弹性力f，为此就需要先计算形变梯度F*，然后计算应变，接着是应力，最后就得到了f。

形变梯度

应变

我们发现，在形变梯度中的单位阵代入应变中后，恰好被消掉了。所以在很多图形学论文中，不是使用形变梯度，而是用形变梯度中的第二项（即去掉单位阵那一项），这一项又恰好是速度的梯度（也是个二阶张量）。

所以应变可以由速度梯度来表示，而不是形变梯度，即

3.3 共轭梯度法

我们采用共轭梯度法求解矩阵。该方法是针对对称矩阵的。而我们的矩阵恰好就是对称矩阵。论证如下：我们为每个粒子赋予相同的静止密度${\rho }_0$，并且以间距（particle spacing）h(译者：我猜是相邻粒子中心的距离)均匀排布我们的粒子。质量、静止密度和间距的关系为$m_0={\rho }_0{h}^3$。给定静止密度和粒子间距，通过这个式子计算出粒子质量。于是粒子体积也是相同的$V_0=m_0/{\rho }_0$。由于有这些假设存在，所以我们的矩阵是对称的。

迭代的停止判据是粒子的平均绝对误差，设定为$1\times 10^{-3}$。为了快速收敛，我们要给一个合适的初值，因此我们用中间速度${\mathbf{v}}^{\mathbf{t}}+\Delta t_c{\mathbf{a}}^{\ast }$来进行热启动（warm start）。

3.4 热启动以加速收敛

计算上节所说的热启动的算法如下所示：

解释如下：

本文的弹性力求解是嵌入到IISPH中的。上述算法绝大部分就是IISPH，只做了一些修改：即我们增加了${\mathbf{V}}_{elas}$的计算。由于弹性力的引入，造成加速度变化，可能会导致原本的CFL条件不再满足，因此要重新计算时间步。为此，

1. 先从现有的时间步中估算出候选时间步$\Delta t_c$。
2. 然后计算其他力的加速度（重力、粘性力、摩擦力）。
3. 提取旋转
4. 计算附加弹性力修正后的速度${\mathbf{V}}_{elas}$
5. 更新加速度a**
6. 从${\mathbf{V}}_{elas}$再次计算时间步$\Delta t$
7. 从$v^t$、$a^{\ast \ast }$和$\Delta t$计算压力和压力梯度力，从而得到$a^{\ast \ast \ast }$
8. 更新速度
9. 更新位置

注意，与原本的线性系统

相类似，给出了弹性力所对应的速度

两式子相比：原本是想要直接求解下一时刻的速度，现在是先计算出一个中间量（${\mathbf{v}}_{elast}$），然后再计算下一时刻速度。

4 Results

首先，看动量守恒。与圣维南模型（也是去除旋转的一个模型）基本相等，可以保证动量的守恒。

上图左边是用了修正的核函数的，右图则没有。左图的方块可以很好地翻转过来。右图的方块翻不过去，大部分旋转被误认为了剪切变形。

上图中，最近的是采用了圣维南模型，中间和最远的都是本文模型。其中最远的时间步调大了五倍。可见大时间步下仍然能够保证不崩溃。并且效果和圣维南模型一样好。

最后是各个场景的测试。包括多相流。相变。布料自碰撞等。

4.1 scences

本文采用的两个固体力学材料参数是剪切模量G和体积模量K。
还可以有其他的组合方式，详见WIKI(下图）

4.1.1 scences: 相变

这里只讲一下相变 phase transition

图7显示了一个场景，兔子溶解成液体，剪切模量和体积模量设置为零。然后，当剪切模量和体积模量增加时，再重新组装。我们的隐式公式能够处理溶解产生的大粒子位移。图8中的场景进一步说明了这一点，四个犰狳在一个盒子里溶解成液体，然后依次重新组装。然而表1中的迭代计数显示，对于大多数其他场景，压力求解器不会受到挑战，它需要在当前场景中增加迭代次数，因为它负责分离犰狳。由于封闭的盒子空间有限，犰狳无法完全展开到休息状态，这进一步加剧了这种情况。