[游戏物理探秘]Neo-Hookean模型模拟弹性物体及unity实现

KaaPexei · 发表于 2022-1-12 12:42

这篇文章主要介绍stvk模型的缺点，以及可压缩Neo-Hookean模型的显隐式方法，以及第二piola-kirchhoff应力。并且收录很多开源库的代码。
demo在
代码分别包括python实现一二三显隐式可压缩Neo-Hookean模型，以及unity实现三维隐式方法，仓库在
或者gitee上
接上篇
https://zhuanlan.zhihu.com/p/446370424
一，重新认识STVK

我们再来复习一下stvk模型。它的能量函数和第一piola-kirchhoff 应力长这个样子

$\Psi(\bold F) = \mu \bold E : \bold E + \frac{\lambda}{2}tr^2(\bold E) \qquad \bold P(\bold F) = \bold F(2\mu \bold E + \lambda tr(\bold E)\bold I) \tag{1.1}$
其中E 是green应变，用于衡量物体的变形程度。物体的变形程度越大，则变形梯度越大，那么green应变越大，进而能量也越大。

$\bold E = \frac{1}{2}(\bold F^T \bold F - \bold I) \tag{1.2}$
上面的公式写成代码大概长下面这个样子。
doubleProduct = 0 # 能量公式中的(E:E)
for i in range(3):
for j in range(3):
doubleProduct += E[i,j]*E[i,j]
trE = E[0,0] + E[1,1] + E[2,2]
psi = mu * doubleProduct + lambda / 2 * trE * trE这在之前已经介绍过了，详细的请看https://zhuanlan.zhihu.com/p/446363864 ，但现在我们将从不同的情况再观察这个stvk模型。
情况一：未变形

当物体处于未变形的状态的时候，变形梯度是单位矩阵，那么格林应变，能量，第一kirchhoff 分别计算如下

$\bold F = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \qquad \bold E = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad \Psi = 0 \qquad \bold P = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$
很好，没有变形的时候，应变，能量密度，以及第一piola kirchhoff 力都是零，没有奇怪的东西。
情况二：每种个方向都拉长

第二种情况是每个方向都拉伸成原来的两倍的时候，那么

$\bold F = \begin{bmatrix}2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{bmatrix} \qquad \bold E = \begin{bmatrix}1.5 & 0 & 0 \\ 0 & 1.5 & 0 \\ 0 & 0 & 1.5 \end{bmatrix} \qquad \Psi = \frac{27}{4}\mu + \frac{81}{8}\lambda \\ \bold P = \begin{bmatrix}6\mu + 9\lambda & 0 & 0 \\ 0 & 6\mu + 9\lambda & 0 \\ 0 & 0 & 6\mu + 9\lambda \end{bmatrix} \\$
看起来也是人畜无害相当正常。既然物体拉长了，那么就肯定产生了能量和节点力，会把物体拉回去，对吧？现在这是三个方向都拉长的情况，如果只拉长一个方向呢？
情况三：仅一个方向拉长

假设仅在x轴拉长两倍，那么变形梯度和piola应力如下

$\bold F = \begin{bmatrix}2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \qquad\bold E = \begin{bmatrix}1.5 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad \bold P = \begin{bmatrix}6\mu + 3\lambda & 0 & 0 \\ 0 & 1.5\lambda & 0 \\ 0 & 0 & 1.5\lambda \end{bmatrix} \\$
变形梯度看着没毛病，应变也只有一个方向不是零，而第一pk力则是三个方向都有数值？这意味着尽管y轴和z轴并没被拉长，但我们有限元模型为了让能量归尽可能小，仍然会让它缩短。换句话说，如果使用stvk模型，那么其中一个方向伸长时，其它的方向就会缩短。在参考[2]中，图书如下

这种特性怎么来的呢？再看一眼piola力的计算方式，我们把第一piola-kirchhoff力拆开，第一项包括一个常数2，一个拉梅系数mu，以及green应变。

$2\mu \bold E \tag{1.3}$
其中拉梅常数mu 也是剪切模量，如下

$\qquad \mu = G = \frac{E}{2(1+\nu)} \tag{1.4}$
其中G是剪切模量，E是杨氏模量，nu是泊松比。剪切模量越大，那么当物体剪切越厉害，受到的力越多。例如当物体产生纯剪切形变时，情况可能如下

$\bold F = \begin{bmatrix}2 & 0 \\ 1 & 2 \end{bmatrix} \qquad 2\mu \bold E = 2\mu\begin{bmatrix}5 & 2 \\ 2 & 4 \end{bmatrix} \tag{1.5}$
这样第一piola-kirchhoff的第一项 2muE 矩阵中每个元素都不是零了，而我们之后可以算出正确的力，逐渐让2 mu E 成为全零矩阵，进而去掉物体的剪切形变。
而第一piola-kirchhoff的第二项，由单位矩阵，拉梅常数lambda，green应变的迹(trace) 以及单位矩阵组成

$\qquad \lambda tr(\bold E)\bold I \tag{1.6}$
其中拉梅系数lambda如下

$\qquad \lambda = \frac{E \nu}{(1+\nu)(1-2\nu)} \tag{1.7}$
简单分析可知道，green应变的迹只关心green 应变主对角线上的元素之和，也就是总拉伸程度。
我们观察之前算出的第一piola - kirchhoff力

$\bold P = \begin{bmatrix}6\mu + 3\lambda & 0 & 0 \\ 0 & 1.5\lambda & 0 \\ 0 & 0 & 1.5\lambda \end{bmatrix} \tag{1.8}$
当我们继续增加mu，那么也仅仅被伸长的x 轴受到的力变大，其它方向不变。
但如果继续增加lambda，则能同时增大xyz三个方向的力。也就是为了让在x轴被伸长的物体的体积恢复到原来的水平，不仅要压缩物体x轴的长度，也要压缩yz轴，也就是变短的同时变细了。
情况四：压缩50%

当三个轴压缩50%的时候，各部分情况如下

$\bold F = \begin{bmatrix}0.5 & 0 & 0 \\ 0 & 0.5 & 0 \\ 0 & 0 & 0.5 \end{bmatrix} \qquad \bold E = \begin{bmatrix}-0.375 & 0 & 0 \\ 0 &-0.375 & 0 \\ 0 & 0 & -0.375 \end{bmatrix} \qquad \Psi = 0.42\mu + 0.63\lambda \\ \bold P = \begin{bmatrix}-0.375\mu - 0.562\lambda & 0 & 0 \\ 0 & -0.375\mu - 0.56\lambda & 0 \\ 0 & 0 & -0.375\mu - 0.56\end{bmatrix} \\$
似乎很正常。
情况五：压缩40%

$\bold F = \begin{bmatrix}0.6 & 0 & 0 \\ 0 & 0.6 & 0 \\ 0 & 0 & 0.6 \end{bmatrix} \qquad \bold E = \begin{bmatrix}-0.32 & 0 & 0 \\ 0 &-0.32 & 0 \\ 0 & 0 & -0.32 \end{bmatrix} \qquad \Psi = 0.31\mu + 0.46\lambda \\ \bold P = \begin{bmatrix}-0.384\mu - 0.576\lambda & 0 & 0 \\ 0 & -0.384\mu - 0.576\lambda & 0 \\ 0 & 0 & -0.384\mu - 0.576\lambda\end{bmatrix} \\$
相比于压缩50%的时候，压缩40%的时候的变形应变和能量的绝对值都要小一点，但是第一piola kirchhoff 应力的绝对值反而大一些。实际上当压缩42%的时候，也就是变形梯度矩阵的主对角线为0.577左右的时候的第一piola kirchhoff 应力最大。这又是piola kirchhoff 应力的计算方式造成的。
假设物体是一维物体，变形梯度是x ，那么第一pk 应力大致的计算公式如下

$y = x((x^2-1) + \frac{1}{2}(x^2-1)) = \frac{3}{2}(x^3 - x) \tag{1.9}$
那么极值点可从以下情况推断出来

$\qquad y' = \frac{9}{2}x^2 - \frac{3}{2} \qquad x = \sqrt{\frac{1}{3}} \tag{1.10}$
对于这种情况，参考[1]是这么说的

这个特点直接让stvk模型面对超过42%的压缩时无能为力，无法计算出更多的力来抵抗压缩。
情况六：压缩成一个点

第四种情况是，但是当物体坍缩成一个点的时候，stvk模型并不容易变回来，甚至可以说根本变不回来，因为

$\bold F = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad \bold E = \begin{bmatrix}-0.5 & 0 & 0 \\ 0 & -0.5 & 0 \\ 0 & 0 & -0.5 \end{bmatrix} \qquad \Psi = \frac{3}{4}\mu + \frac{9}{8}\lambda \\ \bold P = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \\$
既然物体长度缩短100%，但是能量并没有很大，而第一piola kirchhoff 应力甚至全是零，这导致当物体坍缩成点，就永远无法复原了。坍缩成面或线的时候也没法复原。
情况七：完全翻转过来

当物体完全反过来，也就是变形梯度为负的单位矩阵时，那么

$\bold F = \begin{bmatrix}-1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{bmatrix} \qquad \bold E = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad \Psi = 0 \qquad \bold P = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \\$
同样没有能量产生，没有节点力的作用。但此时物体每个轴都是完全反过来，并且毫无变回来的意思。不过这种情况不能完全算是缺点，也许某些时候正需要这些特性。
总结一下，stvk模型的特点如下

不受旋转影响
实现简单
在小变形的时候效果不错
在有压缩的时候，压缩至未变形的58%的时候产生的节点最大，再压缩产生的节点力反而会下降
不能抗极限压缩，也就是当模型坍缩成面，线，点的时候，就没有应力了
当物体完全反转过来的时候，应力无法让物体再返回最初始的状态了。

二，越压缩应力越小的原因

stvk模型中的能量函数用到了green 应变，这样可以去除旋转的影响。但是让green应力对变形梯度求导的结果仍然是变形梯度。那么当变形梯度越小，green应力对变形梯度求导的结果越小，由于stvk模型的公式定义，间接导致piola应力越小。所以，我们为了模拟可以压缩的材质，必须甩掉stvk模型，另寻新欢。
之前说过，能量公式都是人为设计的，都是尽量贴合现实中物体的具体变形。那么现在，就是我们当创世主，自己随心所欲创造物理规则的时候了。连物理规则都能自己创造，这自由度不甩teardown或者minecraft几个宇宙？
不过，有什么函数f(x)，在x的绝对值越小的时候，f(x)的绝对值反而越大呢？
简单起见，我们可以选择变形梯度的逆。这样当压缩到接近一个点的时候，应力就是无穷大了。

$\bold P(\bold F) = \bold F^{-1} \tag{2.1}$
不过选择变形梯度的逆作为应力的一部分有个缺点，那就是无论变形梯度是压缩还是拉伸，应力的方向都相同，并且就算未变形的时候，也仍然有应力。要解决这个问题，我们简单将减去单位矩阵就行了

$\bold P(\bold F) = (\bold F^{-1} -\bold I) \tag{2.2}$
但是变形梯度的逆真正的缺点是，当拉伸程度越大，远远超过一的时候的时候，应力反而越小。这不合理，相当于拆东墙补西墙。
除了逆之外，我们还有一个选择，那就是log函数。
三，可压缩Neo-Hookean

简介

根据参考[3]的6.4.3节，最常用的各项同性可压缩(isotropic compressible )Neo-Hookean模型长下面这个样子，它长下面这个样子。其中J是变形梯度的行列式(determinant)

$\Psi = \frac{\mu}{2}(||\bold F||^2 - 3) - \mu \log J + \frac{\lambda}{2}(\log J)^2 \qquad J = \det(\bold F) \tag{3.1}$
其中logJ 正是我们抵抗压缩所需要的。当变形梯度各项趋近于零，那么F 的行列式也趋近于零，而log J 趋向于负无穷，如此一来就产生很大的能量和应力了。请记住这个公式，因为Neo-Hookean模型有很多种，长得又很像。不过本篇仅仅讨论这个常用的可压缩Neo-Hookean模型。在参考[1]中，计算第一piola-kirchhoff力可以如下计算

$\bold P(\bold F) = \mu(\bold F - \bold F^{-T}) + \lambda \log(J)\bold F^{-T} \tag{3.2}$
可压缩Neo-Hookean的能量公式和应力公式各有什么特性呢？再来逐一分析。
情况一：未变形

第一种情况仍然是未变形

$\bold F = \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad J = 1 \quad \Psi = 0 \quad \bold P = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0& 0 \\ 0 & 0 & 0\end{bmatrix} \\$
这很合理，没有变形的时候，就没有能量，也没用应力产生。
情况二：仅旋转

第一种情况是仅仅发生旋转，那么此时变形梯度就是个旋转矩阵。又因为旋转矩阵是正交矩阵，那么它的行列式det就是 1，并且它的逆等于它的转置。如此一来，能量和应力都是零了。

$\bold F = \bold R \qquad J = 1 \qquad ||\bold F||^2 = 3 \qquad \bold F = \bold F^{-T} \qquad \Psi = 0 \qquad \bold P = \bold 0 \\$
情况三：拉伸一个方向

假设仅有x轴伸长两倍，那么最后算出的应力，仍然会影响三个方向。

$\bold F = \begin{bmatrix}2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad J = 4 \quad \Psi = 0.1\mu + \lambda \quad \bold P = \begin{bmatrix}1.5\mu +0.7\lambda & 0 & 0 \\ 0 & 1.4\lambda & 0 \\ 0 & 0 & 1.4\lambda \end{bmatrix} \\$
情况四：压缩成一个点

这个时候，变形梯度中每个元素都零，那么它行列式也为零。而log 0 为无穷大。运用这种机制，这样无论物体怎么被压扁，总有力量复原。

$\bold F = \begin{bmatrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad J = 0 \quad \Psi = +inf \\$
情况五：完全翻转过来

如果物体完全翻转过来，那么此时logJ处于未定义的状态，这说明使用了logJ 这个函数的Neo-Hookean模型无法处理物体完全翻转的情况。

$\bold F = \begin{bmatrix}-1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{bmatrix} \quad J = -1 \quad \log J = NaN \quad \\$
总结一下，可压缩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[0,0] + FtF[1,1] + FtF[2,2]
# 可压缩 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介绍了。我们希望算的最终结果是应力对定点位置求导，但是我们可将他们拆分成两部分，一部分是应力对变形梯度求导，另一部分是变形梯度对顶点位置求导。

$\frac{\partial \bold P}{\partial \bold x} = \frac{\partial \bold P}{\partial \bold F}\frac{\partial \bold F}{\partial \bold x} \tag{4.1}$
变形梯度对顶点位置求导比较简单。而前者应力对变形梯度求导，就直接复制粘贴[1]的公式了

将上面一堆写成代码在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[0,0] + term[1,1] + term[2,2]) * 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 实际上就是能量对变形梯度求导得到

$\bold P = \frac{\partial \Psi}{\partial \bold F} = \frac{\partial \Psi}{\partial \bold C}\frac{\partial \bold C}{\partial \bold F} = 2\bold F\frac{\partial \bold \Psi}{\partial \bold C} \tag{5.1}$
其中C是右柯西格林变形张量(right Cauchy–Green deformation tensor)，表达式如下，也用于衡量物体的变形程度，并忽视旋转的影响

$\bold C = \bold F^T \bold F \tag{5.2}$
第一piola-kirchhoff应力最大的问题是它是不对称的。设发生带旋转的剪切，并且假设la = mu = 2，那么使用可压缩Neo-Hookean公式可得到下面的结果，

$\bold F = \begin{bmatrix}2 & 0 \\ 1 & 2 \end{bmatrix} \qquad \bold P = \begin{bmatrix}3.69 & 0.15 \\ 2 & 3.69 \end{bmatrix} \tag{5.3}$
也就是只要变形梯度不是对称的，那么第一piola kirchhoff 力就不是对称的。这样一来，各种对于对称矩阵的专用操作都没法用用了。除此之外，第一piola-kirchhoff应力还是作用于变形后的物体上的力。而第二piola kirchhoff 应力不仅是对称的，还是作用于变形前的物体上的力。
关于这个，参考[9]是这么写的

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应变求导得到

$\bold S = \frac{\partial \Psi}{\partial \bold E} = \frac{\partial \Psi}{\partial \bold C}\frac{\partial \bold C}{\partial \bold E} = 2\frac{\partial \bold \Psi}{\partial \bold C} \tag{5.4}$
推导过程可见https://www.bilibili.com/video/BV1zh411s7pR?t=1735.7。对于本篇介绍的可压缩Neo-Hookean模型来说，第二piola-kirohhoff应力如下

$\bold S = \mu(\bold I - \bold C^{-1}) + \lambda (\log J)\bold C^{-1} \tag{5.5}$
那么这样，即使变形梯度不对称，算出来的第二piola kirchhoff 力也是对称的。

$\bold F = \begin{bmatrix}2 & 0 \\ 1 & 2 \end{bmatrix} \qquad \bold S = \begin{bmatrix}2.19 & -0.96 \\ -0.96 & 2.24\end{bmatrix} \tag{5.6}$
将上面这种应力写成代码如下
C = np.dot(F.T,F)
Cinv = np.linalg.inv(C)
J = np.linalg.det(F)
logJ = np.log(J)
# 第一不变量
Ic = C[0,0] + C[1,1] + C[1,1]
# 可压缩 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应力的资料很多，这里就不再赘述了。根据参考[3]的4.5节，可压缩Neo-Hookean的Cauchy应力可以如下计算

$\sigma = \frac{\mu}{J}(\bold b - \bold I) + \frac{\lambda}{J}(\log J)\bold I = J^{-1}\bold F \bold S \bold F^{T} \tag{6.1}$
其中左Cauhy-Green变形张量(left cauchy -green deformation tensor )如下

$\bold b = \bold F \bold F^T \tag{6.2}$
第一piola kirchhoff力与Cauchy应力的关系

$\bold P = J \sigma \bold F^{-T} \qquad J = \det(\bold F) \tag{6.3}$
第二piola kirchhoff力与Cauchy应力的关系

$\bold S = J \bold F^{-1} \sigma \bold F^{-T} \tag{6.4}$
第一piola kirchhoff力与第二piola kirchhoff力的关系

$\bold P = \bold F \bold S \tag{6.5}$
左右CauchyGreen张量的迹

$tr \bold C = tr(\bold F^T \bold F) = tr(\bold F \bold F^T) = tr \bold b \tag{6.6}$
实现了这种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模型，将在另一篇介绍。
参考

[1]FEM Simulation of 3D Deformable Solids: A practitioner’s guide to theory, discretization and model reduction
[2]Dynamic Deformables: Implementation and Production Practicalities
[3] Heyliger, Paul R.. “Nonlinear continuum mechanics for finite element analysis (2nd edn). Javier Bonet and Richard D. Wood, Cambridge University Press, Cambridge, 2008.
[4]Bower, Allan F.. “Applied Mechanics of Solids.” (2009).
[5]https://www.bilibili.com/video/BV1zh411s7pR?t=1735.7
[6]https://www.continuummechanics.org/deformationgradient.html
[7]&#34;Invertible Finite Elements For Robust Simulation of Large Deformation &#34;
[8]https://cn.comsol.com/blogs/why-all-these-stresses-and-strains/
[9]https://engcourses-uofa.ca/books/introduction-to-solid-mechanics/stress/first-and-second-piola-kirchhoff-stress-tensors/
[10]各种开源的github库

JoshWindsor · 发表于 2022-1-12 12:47

[赞][赞]太高产了

Ilingis · 发表于 2022-1-12 12:48

显式方法是把第一第二PK力以及Cauchy力求和作为显式伪代码里的P么

xiangtingsl · 发表于 2022-1-12 12:54

不是，应力求解可以用第一或第二pk应力，选择其中一个即可。显式方法只要把应力求解出来就行了，隐式方法还需要应力对顶点位置求导

		自动登录	找回密码
密码			立即注册

[笔记] [游戏物理探秘]Neo-Hookean模型模拟弹性物体及unity实现

本帖子中包含更多资源