[游戏物理探秘]Neo-Hookean模型模拟弹性物体及unity实现
这篇文章主要介绍stvk模型的缺点，以及可压缩Neo-Hookean模型的显隐式方法，以及第二piola-kirchhoff应力。并且收录很多开源库的代码。demo在
代码分别包括python实现一二三显隐式可压缩Neo-Hookean模型，以及unity实现三维隐式方法，仓库在
或者gitee上
接上篇
https://zhuanlan.zhihu.com/p/446370424
一，重新认识STVK
我们再来复习一下stvk模型。它的能量函数和第一piola-kirchhoff 应力长这个样子
https://www.zhihu.com/equation?tex=%5CPsi%28%5Cbold+F%29+%3D+%5Cmu+%5Cbold+E+%3A+%5Cbold+E+%2B+%5Cfrac%7B%5Clambda%7D%7B2%7Dtr%5E2%28%5Cbold+E%29+%5Cqquad+%5Cbold+P%28%5Cbold+F%29+%3D+%5Cbold+F%282%5Cmu+%5Cbold+E+%2B+%5Clambda+tr%28%5Cbold+E%29%5Cbold+I%29+%5Ctag%7B1.1%7D
其中E 是green应变，用于衡量物体的变形程度。物体的变形程度越大，则变形梯度越大，那么green应变越大，进而能量也越大。
https://www.zhihu.com/equation?tex=%5Cbold+E+%3D+%5Cfrac%7B1%7D%7B2%7D%28%5Cbold+F%5ET+%5Cbold+F+-+%5Cbold+I%29+%5Ctag%7B1.2%7D
上面的公式写成代码大概长下面这个样子。
doubleProduct = 0 # 能量公式中的(E:E)
for i in range(3):
for j in range(3):
doubleProduct += E*E
trE = E + E + E
psi = mu * doubleProduct + lambda / 2 * trE * trE这在之前已经介绍过了，详细的请看https://zhuanlan.zhihu.com/p/446363864 ，但现在我们将从不同的情况再观察这个stvk模型。
情况一：未变形
当物体处于未变形的状态的时候，变形梯度是单位矩阵，那么格林应变，能量，第一kirchhoff 分别计算如下
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D1+%26+0+%26+0+%5C%5C+0+%26+1+%26+0+%5C%5C+0+%26+0+%26+1+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+E+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+%5Cqquad+%5CPsi+%3D+0+%5Cqquad+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+
很好，没有变形的时候，应变，能量密度，以及第一piola kirchhoff 力都是零，没有奇怪的东西。
情况二：每种个方向都拉长
第二种情况是每个方向都拉伸成原来的两倍的时候，那么
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D2+%26+0+%26+0+%5C%5C+0+%26+2+%26+0+%5C%5C+0+%26+0+%26+2+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+E+%3D+%5Cbegin%7Bbmatrix%7D1.5+%26+0+%26+0+%5C%5C+0+%26+1.5+%26+0+%5C%5C+0+%26+0+%26+1.5+%5Cend%7Bbmatrix%7D+%5Cqquad+%5CPsi+%3D+%5Cfrac%7B27%7D%7B4%7D%5Cmu+%2B+%5Cfrac%7B81%7D%7B8%7D%5Clambda+%5C%5C+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D6%5Cmu+%2B+9%5Clambda+%26+0+%26+0+%5C%5C+0+%26+6%5Cmu+%2B+9%5Clambda+%26+0+%5C%5C+0+%26+0+%26+6%5Cmu+%2B+9%5Clambda+%5Cend%7Bbmatrix%7D++%5C%5C
看起来也是人畜无害相当正常。既然物体拉长了，那么就肯定产生了能量和节点力，会把物体拉回去，对吧？现在这是三个方向都拉长的情况，如果只拉长一个方向呢？
情况三：仅一个方向拉长
假设仅在x轴拉长两倍，那么变形梯度和piola应力如下
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D2+%26+0+%26+0+%5C%5C+0+%26+1+%26+0+%5C%5C+0+%26+0+%26+1+%5Cend%7Bbmatrix%7D+%5Cqquad%5Cbold+E+%3D+%5Cbegin%7Bbmatrix%7D1.5+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D6%5Cmu+%2B+3%5Clambda+%26+0+%26+0+%5C%5C+0+%26+1.5%5Clambda+%26+0+%5C%5C+0+%26+0+%26+1.5%5Clambda+%5Cend%7Bbmatrix%7D+%5C%5C
变形梯度看着没毛病，应变也只有一个方向不是零，而第一pk力则是三个方向都有数值？这意味着尽管y轴和z轴并没被拉长，但我们有限元模型为了让能量归尽可能小，仍然会让它缩短。换句话说，如果使用stvk模型，那么其中一个方向伸长时，其它的方向就会缩短。在参考中，图书如下
这种特性怎么来的呢？再看一眼piola力的计算方式，我们把第一piola-kirchhoff力拆开，第一项包括一个常数2，一个拉梅系数mu，以及green应变。
https://www.zhihu.com/equation?tex=2%5Cmu+%5Cbold+E++%5Ctag%7B1.3%7D
其中拉梅常数mu 也是剪切模量，如下
https://www.zhihu.com/equation?tex=%5Cqquad+%5Cmu+%3D+G+%3D+%5Cfrac%7BE%7D%7B2%281%2B%5Cnu%29%7D+%5Ctag%7B1.4%7D
其中G是剪切模量，E是杨氏模量，nu是泊松比。剪切模量越大，那么当物体剪切越厉害，受到的力越多。例如当物体产生纯剪切形变时，情况可能如下
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D2+%26+0+%5C%5C+1+%26+2+%5Cend%7Bbmatrix%7D+%5Cqquad+2%5Cmu+%5Cbold+E+%3D+2%5Cmu%5Cbegin%7Bbmatrix%7D5+%26+2+%5C%5C+2+%26+4+%5Cend%7Bbmatrix%7D+%5Ctag%7B1.5%7D
这样第一piola-kirchhoff的第一项 2muE 矩阵中每个元素都不是零了，而我们之后可以算出正确的力，逐渐让2 mu E 成为全零矩阵，进而去掉物体的剪切形变。
而第一piola-kirchhoff的第二项，由单位矩阵，拉梅常数lambda，green应变的迹(trace) 以及单位矩阵组成
https://www.zhihu.com/equation?tex=%5Cqquad+%5Clambda+tr%28%5Cbold+E%29%5Cbold+I++%5Ctag%7B1.6%7D
其中拉梅系数lambda如下
https://www.zhihu.com/equation?tex=%5Cqquad+%5Clambda+%3D+%5Cfrac%7BE+%5Cnu%7D%7B%281%2B%5Cnu%29%281-2%5Cnu%29%7D+%5Ctag%7B1.7%7D
简单分析可知道，green应变的迹只关心green 应变主对角线上的元素之和，也就是总拉伸程度。
我们观察之前算出的第一piola - kirchhoff力
https://www.zhihu.com/equation?tex=%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D6%5Cmu+%2B+3%5Clambda+%26+0+%26+0+%5C%5C+0+%26+1.5%5Clambda+%26+0+%5C%5C+0+%26+0+%26+1.5%5Clambda+%5Cend%7Bbmatrix%7D+%5Ctag%7B1.8%7D
当我们继续增加mu，那么也仅仅被伸长的x 轴受到的力变大，其它方向不变。
但如果继续增加lambda，则能同时增大xyz三个方向的力。也就是为了让在x轴被伸长的物体的体积恢复到原来的水平，不仅要压缩物体x轴的长度，也要压缩yz轴，也就是变短的同时变细了。
情况四：压缩50%
当三个轴压缩50%的时候，各部分情况如下
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D0.5+%26+0+%26+0+%5C%5C+0+%26+0.5+%26+0+%5C%5C+0+%26+0+%26+0.5+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+E+%3D+%5Cbegin%7Bbmatrix%7D-0.375+%26+0+%26+0+%5C%5C+0+%26-0.375+%26+0+%5C%5C+0+%26+0+%26+-0.375+%5Cend%7Bbmatrix%7D+%5Cqquad+%5CPsi+%3D+0.42%5Cmu+%2B+0.63%5Clambda+%5C%5C+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D-0.375%5Cmu+-+0.562%5Clambda+%26+0+%26+0+%5C%5C+0+%26+-0.375%5Cmu+-+0.56%5Clambda+%26+0+%5C%5C+0+%26+0+%26+-0.375%5Cmu+-+0.56%5Cend%7Bbmatrix%7D+%5C%5C
似乎很正常。
情况五：压缩40%
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D0.6+%26+0+%26+0+%5C%5C+0+%26+0.6+%26+0+%5C%5C+0+%26+0+%26+0.6+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+E+%3D+%5Cbegin%7Bbmatrix%7D-0.32+%26+0+%26+0+%5C%5C+0+%26-0.32+%26+0+%5C%5C+0+%26+0+%26+-0.32+%5Cend%7Bbmatrix%7D+%5Cqquad+%5CPsi+%3D+0.31%5Cmu+%2B+0.46%5Clambda+%5C%5C+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D-0.384%5Cmu+-+0.576%5Clambda+%26+0+%26+0+%5C%5C+0+%26+-0.384%5Cmu+-+0.576%5Clambda+%26+0+%5C%5C+0+%26+0+%26+-0.384%5Cmu+-+0.576%5Clambda%5Cend%7Bbmatrix%7D+%5C%5C
相比于压缩50%的时候，压缩40%的时候的变形应变和能量的绝对值都要小一点，但是第一piola kirchhoff 应力的绝对值反而大一些。实际上当压缩42%的时候，也就是变形梯度矩阵的主对角线为0.577左右的时候的第一piola kirchhoff 应力最大。这又是piola kirchhoff 应力的计算方式造成的。
假设物体是一维物体，变形梯度是x ，那么第一pk 应力大致的计算公式如下
https://www.zhihu.com/equation?tex=y+%3D+x%28%28x%5E2-1%29+%2B+%5Cfrac%7B1%7D%7B2%7D%28x%5E2-1%29%29+%3D+%5Cfrac%7B3%7D%7B2%7D%28x%5E3+-+x%29++%5Ctag%7B1.9%7D
那么极值点可从以下情况推断出来
https://www.zhihu.com/equation?tex=%5Cqquad+y%27+%3D+%5Cfrac%7B9%7D%7B2%7Dx%5E2+-+%5Cfrac%7B3%7D%7B2%7D+%5Cqquad+x+%3D+%5Csqrt%7B%5Cfrac%7B1%7D%7B3%7D%7D+%5Ctag%7B1.10%7D
对于这种情况，参考是这么说的
这个特点直接让stvk模型面对超过42%的压缩时无能为力，无法计算出更多的力来抵抗压缩。
情况六：压缩成一个点
第四种情况是，但是当物体坍缩成一个点的时候，stvk模型并不容易变回来，甚至可以说根本变不回来，因为
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+E+%3D+%5Cbegin%7Bbmatrix%7D-0.5+%26+0+%26+0+%5C%5C+0+%26+-0.5+%26+0+%5C%5C+0+%26+0+%26+-0.5+%5Cend%7Bbmatrix%7D+%5Cqquad+%5CPsi+%3D+%5Cfrac%7B3%7D%7B4%7D%5Cmu+%2B+%5Cfrac%7B9%7D%7B8%7D%5Clambda+%5C%5C+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+%5C%5C
既然物体长度缩短100%，但是能量并没有很大，而第一piola kirchhoff 应力甚至全是零，这导致当物体坍缩成点，就永远无法复原了。坍缩成面或线的时候也没法复原。
情况七：完全翻转过来
当物体完全反过来，也就是变形梯度为负的单位矩阵时，那么
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D-1+%26+0+%26+0+%5C%5C+0+%26+-1+%26+0+%5C%5C+0+%26+0+%26+-1+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+E+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+%5Cqquad+%5CPsi+%3D+0+%5Cqquad+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+%5C%5C
同样没有能量产生，没有节点力的作用。但此时物体每个轴都是完全反过来，并且毫无变回来的意思。不过这种情况不能完全算是缺点，也许某些时候正需要这些特性。
总结一下，stvk模型的特点如下
[*]不受旋转影响
[*]实现简单
[*]在小变形的时候效果不错
[*]在有压缩的时候，压缩至未变形的58%的时候产生的节点最大，再压缩产生的节点力反而会下降
[*]不能抗极限压缩，也就是当模型坍缩成面，线，点的时候，就没有应力了
[*]当物体完全反转过来的时候，应力无法让物体再返回最初始的状态了。
二，越压缩应力越小的原因
stvk模型中的能量函数用到了green 应变，这样可以去除旋转的影响。但是让green应力对变形梯度求导的结果仍然是变形梯度。那么当变形梯度越小，green应力对变形梯度求导的结果越小，由于stvk模型的公式定义，间接导致piola应力越小。所以，我们为了模拟可以压缩的材质，必须甩掉stvk模型，另寻新欢。
之前说过，能量公式都是人为设计的，都是尽量贴合现实中物体的具体变形。那么现在，就是我们当创世主，自己随心所欲创造物理规则的时候了。连物理规则都能自己创造，这自由度不甩teardown或者minecraft几个宇宙？
不过，有什么函数f(x)，在x的绝对值越小的时候，f(x)的绝对值反而越大呢？
简单起见，我们可以选择变形梯度的逆。这样当压缩到接近一个点的时候，应力就是无穷大了。
https://www.zhihu.com/equation?tex=%5Cbold+P%28%5Cbold+F%29+%3D+%5Cbold+F%5E%7B-1%7D+%5Ctag%7B2.1%7D
不过选择变形梯度的逆作为应力的一部分有个缺点，那就是无论变形梯度是压缩还是拉伸，应力的方向都相同，并且就算未变形的时候，也仍然有应力。要解决这个问题，我们简单将减去单位矩阵就行了
https://www.zhihu.com/equation?tex=%5Cbold+P%28%5Cbold+F%29+%3D+%28%5Cbold+F%5E%7B-1%7D+-%5Cbold+I%29+%5Ctag%7B2.2%7D
但是变形梯度的逆真正的缺点是，当拉伸程度越大，远远超过一的时候的时候，应力反而越小。这不合理，相当于拆东墙补西墙。
除了逆之外，我们还有一个选择，那就是log函数。
三，可压缩Neo-Hookean
简介
根据参考的6.4.3节，最常用的各项同性可压缩(isotropic compressible )Neo-Hookean模型长下面这个样子，它长下面这个样子。其中J是变形梯度的行列式(determinant)
https://www.zhihu.com/equation?tex=%5CPsi+%3D+%5Cfrac%7B%5Cmu%7D%7B2%7D%28%7C%7C%5Cbold+F%7C%7C%5E2+-+3%29+-+%5Cmu+%5Clog+J+%2B+%5Cfrac%7B%5Clambda%7D%7B2%7D%28%5Clog+J%29%5E2+%5Cqquad+J+%3D+%5Cdet%28%5Cbold+F%29+%5Ctag%7B3.1%7D
其中logJ 正是我们抵抗压缩所需要的。当变形梯度各项趋近于零，那么F 的行列式也趋近于零，而log J 趋向于负无穷，如此一来就产生很大的能量和应力了。请记住这个公式，因为Neo-Hookean模型有很多种，长得又很像。不过本篇仅仅讨论这个常用的可压缩Neo-Hookean模型。在参考中，计算第一piola-kirchhoff力可以如下计算
https://www.zhihu.com/equation?tex=%5Cbold+P%28%5Cbold+F%29+%3D++%5Cmu%28%5Cbold+F+-++%5Cbold+F%5E%7B-T%7D%29+%2B+%5Clambda+%5Clog%28J%29%5Cbold+F%5E%7B-T%7D+%5Ctag%7B3.2%7D
可压缩Neo-Hookean的能量公式和应力公式各有什么特性呢？再来逐一分析。
情况一：未变形
第一种情况仍然是未变形
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D1+%26+0+%26+0+%5C%5C+0+%26+1+%26+0+%5C%5C+0+%26+0+%26+1+%5Cend%7Bbmatrix%7D+%5Cquad+J+%3D+1+%5Cquad+%5CPsi+%3D+0+%5Cquad++%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0%26+0+%5C%5C+0+%26+0+%26+0%5Cend%7Bbmatrix%7D+%5C%5C
这很合理，没有变形的时候，就没有能量，也没用应力产生。
情况二：仅旋转
第一种情况是仅仅发生旋转，那么此时变形梯度就是个旋转矩阵。又因为旋转矩阵是正交矩阵，那么它的行列式det就是 1，并且它的逆等于它的转置。如此一来，能量和应力都是零了。
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbold+R+%5Cqquad+J+%3D+1+%5Cqquad+%7C%7C%5Cbold+F%7C%7C%5E2+%3D+3+%5Cqquad+%5Cbold+F+%3D+%5Cbold+F%5E%7B-T%7D+%5Cqquad+%5CPsi++%3D+0+%5Cqquad+%5Cbold+P+%3D+%5Cbold+0+%5C%5C
情况三：拉伸一个方向
假设仅有x轴伸长两倍，那么最后算出的应力，仍然会影响三个方向。
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D2+%26+0+%26+0+%5C%5C+0+%26+1+%26+0+%5C%5C+0+%26+0+%26+1+%5Cend%7Bbmatrix%7D+%5Cquad+J+%3D+4+%5Cquad+%5CPsi+%3D+0.1%5Cmu+%2B+%5Clambda+%5Cquad+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D1.5%5Cmu+%2B0.7%5Clambda+%26+0+%26+0+%5C%5C+0+%26+1.4%5Clambda+%26+0+%5C%5C+0+%26+0+%26+1.4%5Clambda+%5Cend%7Bbmatrix%7D+%5C%5C
情况四：压缩成一个点
这个时候，变形梯度中每个元素都零，那么它行列式也为零。而log 0 为无穷大。运用这种机制，这样无论物体怎么被压扁，总有力量复原。
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D+%5Cquad+J+%3D+0+%5Cquad+%5CPsi+%3D+%2Binf++%5C%5C
情况五：完全翻转过来
如果物体完全翻转过来，那么此时logJ处于未定义的状态，这说明使用了logJ 这个函数的Neo-Hookean模型无法处理物体完全翻转的情况。
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D-1+%26+0+%26+0+%5C%5C+0+%26+-1+%26+0+%5C%5C+0+%26+0+%26+-1+%5Cend%7Bbmatrix%7D+%5Cquad+J+%3D+-1+%5Cquad+%5Clog+J+%3D+NaN+%5Cquad++%5C%5C
总结一下，可压缩Neo-Hookean模型的特点如下
[*]不受旋转影响
[*]实现简单
[*]对于压缩变形有比较好的表现
[*]模型不能翻转过来，否则变形梯度的行列式小于零，将导致log J 未定义
将上面一堆东西，写成代码在https://github.com/clatterrr/PhysicsSimulationPlayground/blob/master/ImplictElasticFem/NeohookeanPython/NeoHookeanExplicit3D.py，核心部分如下如下
FinvT = np.linalg.inv(F).T
FtF = np.dot(F.T,F)
J = max(np.linalg.det(F),0.001)
logJ = np.log(J)
# 第一不变量
Ic = FtF + FtF + FtF
# 可压缩 neohookean 能量
energy = mu * 0.5 * (Ic - 3) - mu * logJ + la * 0.5 * logJ * logJ
# 第一 piola kirchhoff 应力
piola = mu * F - mu * FinvT + la * logJ * 0.5 * FinvT分析完毕，接下来就是去看看世界前几名的开源库是怎么实现可压缩Neo-Hookean的了。当然，我觉得单纯的公式实现并不值得看，值得看的是各个开源库在计算公式之前做了什么准备，计算完公式后，又把结果作了什么处理，送到哪里去了。这些公式会以各种各样形式的代码出现，不过想找出Neo-Hookean模型的代码，认准logJ就行了。
polyfem
在https://github.com/polyfem/polyfem/blob/master/src/assembler/NeoHookeanElasticity.cpp#L220对能量计算公式写的很清楚
const T log_det_j = log(polyfem::determinant(def_grad));
const T val = mu / 2 * ((def_grad.transpose() * def_grad).trace() - size() - 2 * log_det_j)
+ lambda / 2 * log_det_j * log_det_j;对应力计算公式写的很清楚
//stress = mu (F - F^{-T}) + lambda ln J F^{-T}
//stress = mu * (def_grad - def_grad^{-T}) + lambda ln (det def_grad) def_grad^{-T}
Eigen::MatrixXd stress_tensor = mu * (def_grad - FmT) + lambda * std::log(def_grad.determinant()) * FmT;FEbio
在https://github.com/febiosoftware/FEBio/blob/master/FEBioMech/FENeoHookean.cpp#L112 的能量公式计算写得很详细。
double sed = mu*((I1-3)/2.0 - lnJ)+lam*lnJ*lnJ/2.0;IPC
在https://github.com/ipc-sim/IPC/blob/master/src/Energy/Physics_Elasticity/NeoHookeanEnergy.cpp#L68、对能量计算公式写得很详细。
const double sigma2Sum = singularValues.squaredNorm();
const double sigmaProd = singularValues.prod();
const double log_sigmaProd = std::log(sigmaProd);
E = u / 2.0 * (sigma2Sum - dim) - (u - lambda / 2.0 * log_sigmaProd) * log_sigmaProd;那个奇异值和变形梯度的奇异值分解有关。所以我们暂时就把它当成变形梯度就行。之后计算应力时也用到了奇异值分解，所以到时候再写得详细一些。计算应力所用代码在
dE_div_dF = u * (F - FInvT) + lambda * std::log(J) * FInvT;IGSim
这是一个c++有限元库，在
https://github.com/milkpku/IGsim/blob/master/include/IGsim/neohookean_model.cpp#L52 对能量公式及应力公式写得很清楚。
/* energy(F) = \mu/2 * (I_1 - log I_3 - 3) + \lambda/8 * log^2 I_3 */
energy = 0.5 * Mu * (I1 - 2 * logJ - 3) + 0.5 * Lam * logJ * logJ;
/* P(F) = \mu * (F - F^{-T}) + \lambda/2 * log I_3 * F^{-T} */
DerivedF FinvT = F.inverse().transpose();
P = Mu * (F - FinvT) + Lam * logJ * FinvT;Q-Minh
这是一个positionBasedDynamics库，在https://github.com/Q-Minh/position-based-dynamics/blob/23fcf93bddd5daf425cdc3443da05760cc1343ec/src/xpbd/neohookean_elasticity_constraint.cpp#L136 中对能量公式和应力公式写得很清楚
// psi(I1, J) = (mu/2)*(I1 - 3) - mu*log(J) + (lambda/2)*log^2(J)
psi = static_cast<scalar_type>(0.5) * mu_ * (I1 - static_cast<scalar_type>(3.)) -
mu_ * logJ + static_cast<scalar_type>(0.5) * lambda_ * logJ * logJ;
// P(F) = mu*(F - mu*F^-T) + lambda*log(J)*F^-T
Piola = mu_ * (F - mu_ * FinvT) + lambda_ * logJ * FinvT;四，隐式方法
显式方法容易计算，但不稳定，因此可以使用隐式方法。这个我已经在https://zhuanlan.zhihu.com/p/446370424介绍了。我们希望算的最终结果是应力对定点位置求导，但是我们可将他们拆分成两部分，一部分是应力对变形梯度求导，另一部分是变形梯度对顶点位置求导。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+%5Cbold+P%7D%7B%5Cpartial+%5Cbold+x%7D+%3D+%5Cfrac%7B%5Cpartial+%5Cbold+P%7D%7B%5Cpartial+%5Cbold+F%7D%5Cfrac%7B%5Cpartial+%5Cbold+F%7D%7B%5Cpartial+%5Cbold+x%7D+%5Ctag%7B4.1%7D
变形梯度对顶点位置求导比较简单。而前者应力对变形梯度求导，就直接复制粘贴的公式了
将上面一堆写成代码在https://github.com/clatterrr/PhysicsSimulationPlayground/blob/master/ImplictElasticFem/NeohookeanPython/NeohookeanImplicit3D.py，核心部分如下
for i in range(12):
dF = np.dot(dD,Dminv)
dP = mu * dF
dP += (mu - la * logJ) * np.dot(np.dot(FinvT,dF.T),FinvT)
term = np.dot(Finv,dF)
dP += la * (term + term + term) * FinvT
dH = - volume * np.dot(dP,Dminv.T)其它开源库是怎么写的呢？
QuasiNewton
地址在https://github.com/ltt1598/Quasi-Newton-Metods-for-Real-time-Simulation-of-Hyperelastic-Materials/blob/master/GenPD/GenPD/source/constraint_tet.cpp#L796，核心代码如下
dPdF(i, j) = m_mu * deltaF;
const ScalarType& J0 = m_neohookean_clamp_value;
if (J > J0){
logJ = std::log(J);
dPdF(i, j) += (m_mu - m_lambda*logJ)*FinvT*deltaF.transpose()*FinvT
+ m_lambda*(Finv*deltaF).trace() * FinvT;
}MIT-6.807
这是个作业仓库。在https://github.com/BrandonPreble/MIT-6.807-Code/blob/f520f9b0699c6d36951b40b742b6c67cb2f19f1c/Common/FEM/neohookean_material.hpp#L50 写出了应力求导的公式
return mu * dF
+ (mu - lambda * log(J)) * F_inv_trans * dF.transpose()
* F_inv_trans
+ lambda * (F_inv * dF).trace() * F_inv_trans;ziran
地址在https://github.com/penn-graphics-research/ziran2020/blob/master/Lib/Ziran/Physics/ConstitutiveModel/NeoHookean.cpp#L62，主要代码在
T scale = lambda * s.logJ - mu;
dP = mu * dF + (lambda - scale) * s.FinvT.cwiseProduct(dF).sum() * s.FinvT;
dP += scale / s.J * dP;Otherlab
地址在https://github.com/otherlab/geode/blob/157cc904c113cc5e29a1ffe7c091a83b8ec2cf8e/geode/force/NeoHookean.cpp#L117，主要代码为
T mu_minus_lambda_logJ = mu()+lambda()*log(F_inverse.determinant());
SymmetricMatrix<T,3> F_inverse_outer = outer_product(F_inverse.to_vector());
DiagonalizedIsotropicStressDerivative<T,3> dP_dF;
dP_dF.x0000 = mu()+(lambda()+mu_minus_lambda_logJ)*F_inverse_outer.x00;
dP_dF.x1111 = mu()+(lambda()+mu_minus_lambda_logJ)*F_inverse_outer.x11;
dP_dF.x2222 = mu()+(lambda()+mu_minus_lambda_logJ)*F_inverse_outer.x22;
dP_dF.x1100 = lambda()*F_inverse_outer.x10;
dP_dF.x2200 = lambda()*F_inverse_outer.x20;
dP_dF.x2211 = lambda()*F_inverse_outer.x21;
dP_dF.x1010 = dP_dF.x2020 = dP_dF.x2121 = mu();
dP_dF.x1001 = mu_minus_lambda_logJ*F_inverse_outer.x10;
dP_dF.x2002 = mu_minus_lambda_logJ*F_inverse_outer.x20;
dP_dF.x2112 = mu_minus_lambda_logJ*F_inverse_outer.x21;五，对称的第二piola-kirchhoff应力
第一piola-kirohhoff 实际上就是能量对变形梯度求导得到
https://www.zhihu.com/equation?tex=%5Cbold+P+%3D+%5Cfrac%7B%5Cpartial+%5CPsi%7D%7B%5Cpartial+%5Cbold+F%7D+%3D+%5Cfrac%7B%5Cpartial+%5CPsi%7D%7B%5Cpartial+%5Cbold+C%7D%5Cfrac%7B%5Cpartial+%5Cbold+C%7D%7B%5Cpartial+%5Cbold+F%7D+%3D+2%5Cbold+F%5Cfrac%7B%5Cpartial+%5Cbold+%5CPsi%7D%7B%5Cpartial+%5Cbold+C%7D+%5Ctag%7B5.1%7D
其中C是右柯西格林变形张量(right Cauchy–Green deformation tensor)，表达式如下，也用于衡量物体的变形程度，并忽视旋转的影响
https://www.zhihu.com/equation?tex=%5Cbold+C+%3D+%5Cbold+F%5ET+%5Cbold+F+%5Ctag%7B5.2%7D
第一piola-kirchhoff应力最大的问题是它是不对称的。设发生带旋转的剪切，并且假设la = mu = 2，那么使用可压缩Neo-Hookean公式可得到下面的结果，
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D2+%26+0++%5C%5C+1+%26+2+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+P+%3D+%5Cbegin%7Bbmatrix%7D3.69+%26+0.15++%5C%5C+2+%26+3.69+%5Cend%7Bbmatrix%7D+%5Ctag%7B5.3%7D
也就是只要变形梯度不是对称的，那么第一piola kirchhoff 力就不是对称的。这样一来，各种对于对称矩阵的专用操作都没法用用了。除此之外，第一piola-kirchhoff应力还是作用于变形后的物体上的力。而第二piola kirchhoff 应力不仅是对称的，还是作用于变形前的物体上的力。
关于这个，参考是这么写的
The first Piola-Kirchhoff stress matrix can intuitively be viewed as the force in the deformed configuration per unit area of the undeformed configuration, while the second Piola-Kirchhoff stress can intuitively be viewed as the force in the undeformed configuration per unit area of the undeformed configuration.至于为什么是这样我没看懂，以后再慢慢理解。第二piola-kirohhoff 实际上就是能量对green应变求导得到
https://www.zhihu.com/equation?tex=%5Cbold+S+%3D+%5Cfrac%7B%5Cpartial+%5CPsi%7D%7B%5Cpartial+%5Cbold+E%7D+%3D+%5Cfrac%7B%5Cpartial+%5CPsi%7D%7B%5Cpartial+%5Cbold+C%7D%5Cfrac%7B%5Cpartial+%5Cbold+C%7D%7B%5Cpartial+%5Cbold+E%7D+%3D+2%5Cfrac%7B%5Cpartial+%5Cbold+%5CPsi%7D%7B%5Cpartial+%5Cbold+C%7D+%5Ctag%7B5.4%7D
推导过程可见https://www.bilibili.com/video/BV1zh411s7pR?t=1735.7。对于本篇介绍的可压缩Neo-Hookean模型来说，第二piola-kirohhoff应力如下
https://www.zhihu.com/equation?tex=%5Cbold+S+%3D+%5Cmu%28%5Cbold+I+-+%5Cbold+C%5E%7B-1%7D%29+%2B+%5Clambda+%28%5Clog+J%29%5Cbold+C%5E%7B-1%7D+%5Ctag%7B5.5%7D
那么这样，即使变形梯度不对称，算出来的第二piola kirchhoff 力也是对称的。
https://www.zhihu.com/equation?tex=%5Cbold+F+%3D+%5Cbegin%7Bbmatrix%7D2+%26+0++%5C%5C+1+%26+2+%5Cend%7Bbmatrix%7D+%5Cqquad+%5Cbold+S+%3D+%5Cbegin%7Bbmatrix%7D2.19+%26+-0.96+%5C%5C+-0.96+%26+2.24%5Cend%7Bbmatrix%7D+%5Ctag%7B5.6%7D
将上面这种应力写成代码如下
C = np.dot(F.T,F)
Cinv = np.linalg.inv(C)
J = np.linalg.det(F)
logJ = np.log(J)
# 第一不变量
Ic = C + C + C
# 可压缩 neohookean 能量
energy = mu * 0.5 * (Ic - 3) - mu * logJ + la * 0.5 * logJ * logJ
# 第二 piola kirchhoff 应力
S = mu * (np.identity(3) - Cinv) + la * logJ * Cinv实现了这种第二piola-kirchhoff的开源库也有很多，如下。
Akantu
地址，https://github.com/sunscale/Akantu/blob/master/src/model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.hh#L47
计算第二piola-kirchhoff 应力的代码如下
this->rightCauchy(F, C);
Real J = F.det() * sqrt(C33);
Cminus.inverse(C);
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
S(i, j) = (i == j) * mu + (lambda * log(J) - mu) * Cminus(i, j);AsFem
这也是有限元库，在B站上有账号。https://space.bilibili.com/100272198?spm_id_from=333.788.b_765f7570696e666f.2。其中https://github.com/M3Group/AsFem/blob/main/src/MateSystem/NeoHookeanMaterial.cpp#L58的代码为
_pk2=(_I-_Cinv)*mu+_Cinv*lambda*log(J);cubica
耶鲁大学图形组助理教授Theodore Kim的开源库Cubica。他的个人主页http://www.tkim.graphics/，Cubica主页http://www.tkim.graphics/cubica/。其中的主要代码为
MATRIX3 C = F.transpose() * F;
MATRIX3 Cinverse = C.inverse();
const Real J = det(F);
return _mu * (MATRIX3::I() - Cinverse) + _lambda * (log(J)) * Cinverse;fabsim
地址在https://github.com/DavidJourdan/fabsim/blob/master/src/NeoHookeanElement.cpp#L47，主要代码为
Matrix<double, 3, 2> F = deformationGradient(X);
Matrix2d C = F.transpose() * F;
double lnJ = log(C.determinant()) / 2;
return mu * Matrix2d::Identity() + (lambda * lnJ - mu) * C.inverse();FEBio
地址在https://github.com/febiosoftware/FEBio/blob/master/FEBioMech/FENeoHookean.cpp#L124，主要代码为
mat3ds C = I + ES*2;
mat3ds Ci = C.inverse();
double detF = sqrt(C.det());
double lndetF = log(detF);
// calculate stress
mat3ds S = (I - Ci)*mu + Ci*(lam*lndetF);Physika
地址在https://github.com/FeiZhu/Physika/blob/7887e1c88022a2065ff651b75d53f06f1e25b951/Physika_Src/Physika_Dynamics/Constitutive_Models/neo_hookean.cpp#L89，主要代码在
SquareMatrix<Scalar,Dim> identity = SquareMatrix<Scalar,Dim>::identityMatrix();
SquareMatrix<Scalar,Dim> inverse_c = (F.transpose()*F).inverse();
Scalar lnJ = log(F.determinant());
Scalar mu = this->mu_;
Scalar lambda = this->lambda_;
SquareMatrix<Scalar,Dim> S = mu*(identity-inverse_c)+lambda*lnJ*inverse_c;六，Cauchy应力
除了第一二piola-kirchhoff应力之外，另一个很重要的就是Cauchy应力。Cauchy应力的资料很多，这里就不再赘述了。根据参考的4.5节，可压缩Neo-Hookean的Cauchy应力可以如下计算
https://www.zhihu.com/equation?tex=%5Csigma+%3D+%5Cfrac%7B%5Cmu%7D%7BJ%7D%28%5Cbold+b+-+%5Cbold+I%29+%2B+%5Cfrac%7B%5Clambda%7D%7BJ%7D%28%5Clog+J%29%5Cbold+I+%3D+J%5E%7B-1%7D%5Cbold+F+%5Cbold+S+%5Cbold+F%5E%7BT%7D+%5Ctag%7B6.1%7D
其中左Cauhy-Green变形张量(left cauchy -green deformation tensor )如下
https://www.zhihu.com/equation?tex=%5Cbold+b+%3D+%5Cbold+F+%5Cbold+F%5ET+%5Ctag%7B6.2%7D
第一piola kirchhoff力与Cauchy应力的关系
https://www.zhihu.com/equation?tex=%5Cbold+P+%3D+J+%5Csigma+%5Cbold+F%5E%7B-T%7D+%5Cqquad+J+%3D+%5Cdet%28%5Cbold+F%29+%5Ctag%7B6.3%7D
第二piola kirchhoff力与Cauchy应力的关系
https://www.zhihu.com/equation?tex=%5Cbold+S+%3D+J+%5Cbold+F%5E%7B-1%7D+%5Csigma+%5Cbold+F%5E%7B-T%7D+%5Ctag%7B6.4%7D
第一piola kirchhoff力与第二piola kirchhoff力的关系
https://www.zhihu.com/equation?tex=%5Cbold+P+%3D+%5Cbold+F+%5Cbold+S+%5Ctag%7B6.5%7D
左右CauchyGreen张量的迹
https://www.zhihu.com/equation?tex=tr+%5Cbold+C+%3D+tr%28%5Cbold+F%5ET+%5Cbold+F%29+%3D+tr%28%5Cbold+F+%5Cbold+F%5ET%29+%3D+tr+%5Cbold+b+%5Ctag%7B6.6%7D
实现了这种Cauchy应力的开源库包括
Cubica
耶鲁大学图形组助理教授Theodore Kim的开源库Cubica。他的个人主页http://www.tkim.graphics/，Cubica主页http://www.tkim.graphics/cubica/。其中在Material.cpp的主要代码为
MATRIX3 MATERIAL::cauchyStress(MATRIX3& F, bool diagonal)
{
Real J = det(F);
MATRIX3 S = secondPiolaKirchhoff(F, diagonal);
return (1.0 / J) * F * S * F.transpose();
}FEMBio
在https://github.com/febiosoftware/FEBio/blob/master/FEBioMech/FENeoHookean.cpp#L40的主要的代码为
double detF = pt.m_J;
double detFi = 1.0/detF;
double lndetF = log(detF);
// calculate left Cauchy-Green tensor
mat3ds b = pt.LeftCauchyGreen();
// calculate stress
mat3ds s = (b - I)*(mu*detFi) + I*(lam*lndetF*detFi);fea-large
在https://github.com/fourier/fea-large/blob/master/solver-prototype/cylindric/materials/sigma_neohookean_compressible.m的主要代码为。
J = det(F);
b = F*F&#39;;
E = eye(3,3);
C = b-E;
% Sigma = lambda/J*(ln(J))E+mu/J*(b-E)
S = g_mu/J*C+g_lambda/J*log(J)*E;Neo-Hookean模型有很多种，除了可压缩Neo-Hookean外，另一种常用的是stable Neo-Hookean模型，将在另一篇介绍。
参考
FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction
Dynamic Deformables: Implementation and Production Practicalities
Heyliger, Paul R.. “Nonlinear continuum mechanics for finite element analysis (2nd edn). Javier Bonet and Richard D. Wood, Cambridge University Press, Cambridge, 2008.
Bower, Allan F.. “Applied Mechanics of Solids.” (2009).
https://www.bilibili.com/video/BV1zh411s7pR?t=1735.7
https://www.continuummechanics.org/deformationgradient.html
&#34;Invertible Finite Elements For Robust Simulation of Large Deformation &#34;
https://cn.comsol.com/blogs/why-all-these-stresses-and-strains/
https://engcourses-uofa.ca/books/introduction-to-solid-mechanics/stress/first-and-second-piola-kirchhoff-stress-tensors/
各种开源的github库 [赞][赞]太高产了 显式方法是把第一第二PK力以及Cauchy力求和作为显式伪代码里的P么 不是，应力求解可以用第一或第二pk应力，选择其中一个即可。显式方法只要把应力求解出来就行了，隐式方法还需要应力对顶点位置求导
页:
[1]