CFD: 线弹性方程
大部分的固体应力分析都使用的是有限元。然而使用有限体积法计算应力应变以及位移也是可行的。同时有限体积法具有一些优良的优点：如守恒特性。在将力离散到网格单元面的时候，所有的力均被抵消仅仅在边界上存在受力。这与与流体的质量守恒类似。OpenFOAM中的solidDisplacementFoam则采用有限体积法来对位移矢量\(\bfD\)方程进行离散，进而获得体对称应变张量\(\boldsymbol{\sigma}\)。本文对其简述。注意本文的方程只适用于小变形，因为只有在这种情况下，不存在非线性的类似流体的对流非线性项。
控制方程与算法
在不考虑体积力的情况下，线弹性固体的控制方程为：
(1)\[
\frac{\p ^2 \bfD}{\p t ^2}-\nabla\cdot\boldsymbol{\sigma}=0
\]
其中其中\(\bfD\)表示位移矢量，体心的应变张量\(\boldsymbol{\sigma}\)可以写为：
(2)\[
\boldsymbol{\sigma}=\mu\left(\nabla\bfD+\nabla\bfD^T\right) +\lambda (\nabla\cdot\bfD)\bfI
\]
(3)\[
\mu=\frac{E}{2(1 + \nu)}
\]
(4)\[
\lambda=\frac{\nu E}{(1 + \nu)(1 - 2\nu)}
\]
\(E\)表示杨氏模量，\(\lambda\)表示拉梅系数的第1个参数，\(\mu\)表示拉梅系数的第2个参数（在流体力学中理解为粘度，在固体力学中理解为剪切模量）。杨氏模量越大，越不容易发生变形。比如钢的杨氏模量大约是e11量级，普通的麻绳杨氏模量大约是e9量级，比较硬的绳子可能是e10量级。\(\nu\)是泊松比，表示当一个材料被拉伸时，横向应变与轴向应变的比值。钢的泊松比在0.2-0.3量级。木头的泊松比基本为0。同时，钢的密度大约在8000左右。
将方程(2)代入到(1)，有控制方程：
(5)\[
\frac{\p ^2 \bfD}{\p t ^2}-\nabla\cdot\left( \mu(\nabla \bfD+\nabla \bfD^T) \right) + \nabla\cdot\left(\lambda (\nabla\cdot\bfD)\bfI \right)=0
\]
方程(5)求解的为位移矢量。因此，在给定边界条件的时候，需要给定位移矢量的边界条件。例如：边界处的固定变形是多少之类。
方程(1)也可以写为：
(6)\[
\frac{\p ^2 \bfD}{\p t ^2}-\frac{1}{\dV}\oint\boldsymbol{\sigma}_f\cdot\bfn_f\rd S=0
\]
\(\bfn=\frac{\bfS_f}{|\bfS_f|}\)。其中的\(\boldsymbol{\sigma}_f\)从方程(2)插值而来。\(\boldsymbol{\sigma}_f\)具有6个分量，但\(\boldsymbol{\sigma}_f\cdot\bfn_f\)具有3个分量（矢量），也即\(\boldsymbol{\sigma}_f\)的6个分量有3个不参与离散计算。方程(6)离散之后，需要对面上的\(\boldsymbol{\sigma}_f\cdot\bfn_f\)进行加和计算。其实\(\boldsymbol{\sigma}_f\cdot\bfn_f\)对应的是物体表面的受力。这里对应下文讨论的牵引力压力边界条件。
离散1
仔细观察方程(5)，其实际上表示3个位移矢量的分量的方程。每个位移矢量分量的方程里面，涉及到其他的位移矢量分量。传统的FEM算法，是采用耦合求解的方法，离散成一个大的方程组同时求解。在CFD里面也有类似的求解算法，即为耦合算法。
本文讨论分离式求解，即将三个位移矢量分量的方程离散为三个分量方程。这与速度矢量的方程离散是一样的。采用这种离散方式可以得到3个稀疏矩阵，稀疏矩阵的好处之一是可以降低内存存储，另外是可以调用迭代矩阵求解器。然而分离式求解一些变量需要做显性离散处理，因此分离式求解不仅需要用稀疏线性系统求解器来迭代求解矩阵，还需要迭代求解位移矢量分量在三个方程之间的耦合关系。
显而易见的，在求解方程(5)的时候，其中显性与隐性离散的处理关系可以写为：
(7)\[\underset{\mathrm{implicit}}{\underbrace{\frac{\p ^2 \bfD}{\p t ^2}-\nabla\cdot ( \mu\nabla \bfD)}}
\underset{\mathrm{explicit}}{\underbrace{-\nabla\cdot ( \mu\nabla \bfD^T)-\nabla\cdot\left(\lambda (\nabla\cdot\bfD)\bfI \right)}}
=0
\]
然而在实际过程中，方程(7)的离散过程收敛性非常差。这是因为显性项的贡献比隐性项的贡献还要大。因此，OpenFOAM将其改写为：
(8)\[\begin{split}
\begin{multline}
\underset{\mathrm{implicit}}{\underbrace{\frac{\p ^2 \bfD}{\p t ^2}-\nabla\cdot ( (2\mu+\lambda)\nabla \bfD)}}
\\
\underset{\mathrm{explicit}}{\underbrace{-\nabla\cdot ( \mu(\nabla\bfD+\nabla\bfD^T)-\nabla\cdot\left(\lambda (\nabla\cdot\bfD)\bfI \right)
+\nabla\cdot((2\mu+\lambda)\nabla\bfD)}}
=0
\end{multline}
\end{split}\]
依据方程(2)，OpenFOAM将进一步化简：
(9)\[
\underset{\mathrm{implicit}}{\underbrace{\frac{\p ^2 \bfD}{\p t ^2}-\nabla\cdot ( (2\mu+\lambda)\nabla \bfD)}}
\underset{\mathrm{explicit}}{\underbrace{-\nabla\cdot \boldsymbol{\sigma}
+\nabla\cdot((2\mu+\lambda)\nabla\bfD)}}
=0
\]
其中\(2\mu\)中的\(2\)可以任意改变，\(2\mu+\lambda\)其实可以任意的改变，其与收敛性有关联。
Note
方程(9)在一个时间步内要多次求解，来减少显性离散带来的滞后，这样更趋向于空间项全隐格式。
压力显性分离
对于连续介质，也可以将\(\boldsymbol{\sigma}\)处理为偏应力与正应力，前者处理剪切变形，后者处理法向挤压拉伸，有：
\[
\boldsymbol{\sigma}= \mathrm{dev}(\boldsymbol{\sigma}) - p\bfI
\]
在这种情况下，方程(5)可以写为：
(10)\[
\underset{\mathrm{implicit}}{\underbrace{\frac{\p ^2 \bfD}{\p t ^2}-\nabla\cdot ( (2\mu+\lambda)\nabla \bfD)}}
\underset{\mathrm{explicit}}{\underbrace{-\nabla\cdot \mathrm{dev}\left(\boldsymbol{\sigma}\right)
+\nabla\cdot((2\mu+\lambda)\nabla\bfD)}}
=-\nabla p
\]
其中的压力项显性计算。方程(9)与方程(10)无本质差异。类似的，方程(10)可以写为\(\bfU\)的形式：
(11)\[
\underset{\mathrm{implicit}}{\underbrace{\frac{\p \rho \bfU}{\p t }-\nabla\cdot ( \dt(2\mu+\lambda)\nabla \bfU)}}
\underset{\mathrm{explicit}}{\underbrace{- \nabla\cdot \mathrm{dev}\left(\boldsymbol{\sigma}\right)
+\nabla\cdot(\dt(2\mu+\lambda)\nabla\bfU)}}
=-\nabla p
\]
如果将压力显性离散，同时给定\(\bfU\)的边界条件，方程(11)可以用来求解。方程(11)与方程(9)无本质差异。
压力隐性离散
同时对于不可压缩固体，有（压力方程的第1种形式）：
(12)\[
\frac{p}{\lambda+\frac{2}{3}\mu}+\nabla\cdot\bfD=0
\]
压力方程也可以写为近似流体力学的连续性方程（压力方程的第2种形式）：
(13)\[
\frac{\p\rho}{\p t}+\nabla\cdot(\rho\bfU)=0
\]
如果将方程(10)或(11)中的压力隐性离散，其可以和方程(12)或方程(13)可以构成耦合系统，在指定压力边界条件后，可通过迭代方法进行求解。该方法与压力显性方法存在数值差异，结果会更精准。
边界条件
针对\(\bfD\)求解的情况下，需要给定\(\bfD\)的边界条件。\(\bfD\)如果给定固定值边界条件，则假定其固定位移，在这种情况下表面的应力可能会非均一分布。对于任意的固体表面，可给定牵引力\(\mathbf{T}_f\)的边界条件。牵引力是定义在面上的一个力，存在三个方向。给定具体的面，三个分量分别对应\(x,y,z\)方向的牵引力。牵引力会导致形变。除了牵引力外，还存在法向压力\(p^{udf}_f\)。其作用于面的法向，因此是一个标量。牵引力和压力都会导致物体发生形变。
压力\(p_f^{udf}\bfn\)（正值向内）与牵引力\(\mathbf{T}_f\)（正值向外）反向。\(p^{udf}_f\)与\(\mathbf{T}_f\)的存在会使物体发生形变。方程(9)可以用来构建边界条件，但其中并没有\(\bfD\)以及\(\nabla\bfD\)的存在，因此从方程(9)中构建第一二类边界条件并不是显而易见的。
另一方面，可以从物理的角度来构建边界条件，对于给定表面，\(p^{udf}\)与\(\mathbf{T}\)的定义为：
(14)\[
\mathbf{T}_f-p_f^{udf}\bfn_f=\boldsymbol{\sigma}_f\cdot\bfn_f
\]
其可以用来构建\(\bfD\)的边界条件：
(15)\[
\mathbf{T}_f-p_f^{udf}\bfn_f= \mu\left(\nabla\bfD+\nabla\bfD^T\right)_f\cdot\bfn_f +\lambda (\nabla\cdot\bfD)_f\bfI \cdot\bfn_f
\]
有：
(16)\[
(\nabla\bfD)_f\cdot\bfn_f=
\frac
{\mathbf{T}_f-p_f^{udf}\bfn_f -\lambda (\nabla\cdot\bfD)_f\bfI \cdot\bfn_f - \mu(\nabla\bfD^T)_f\cdot\bfn_f}
{\mu}
\]
方程(9)的另一种理解方法即为\(\boldsymbol{\sigma}\)分解为显性部分以及隐性部分：
(17)\[
\boldsymbol{\sigma}=\underset{\mathrm{implicit}}{\underbrace{ (2\mu+\lambda)\nabla \bfD}}
+
\underset{\mathrm{explicit}}{\underbrace{\boldsymbol{\sigma}
- (2\mu+\lambda)\nabla\bfD }}
\]
将方程(17)代入到方程(14)有
(18)\[
\mathbf{T}_f-p^{udf}_f\bfn_f= \left(\underset{\mathrm{implicit}}{\underbrace{ (2\mu+\lambda)\nabla \bfD}}
+
\underset{\mathrm{explicit}}{\underbrace{\boldsymbol{\sigma}
- (2\mu+\lambda)\nabla\bfD}}\right)_f\cdot\bfn_f
\]
(19)\[
\mathbf{T}_f-p_f^{udf}\bfn_f= \underset{\mathrm{implicit}}{\underbrace{ (2\mu+\lambda)(\nabla \bfD)_f\cdot\bfn_f}}
+
\underset{\mathrm{explicit}}{\underbrace{\boldsymbol{\sigma}_f\cdot\bfn_f
- (2\mu+\lambda)(\nabla\bfD)_f\cdot\bfn_f }}
\]
其中的隐性部分可以用来构建边界条件：
(20)\[
(\nabla\bfD)_f \cdot\bfn_f
=
\frac{
\mathbf{T}_f-p_f^{udf}\bfn_f- \boldsymbol{\sigma}_f \cdot\bfn_f
+(2\mu+\lambda)(\nabla\bfD)_f \cdot\bfn_f
}{2\mu+\lambda}
\]
从方程(20)可以看出，给定牵引力\(\mathbf{T}_f\)的时候，\(\bfD\)会形成一个法向梯度固定值边界条件，其值通过方程(20)计算而来。对于自由表面，牵引力\(\mathbf{T}_f\)为\((0, 0, 0)\)。对于固定牵引力边界，需要给定牵引力\(\mathbf{T}_f\)的值，然后方程(20)会自动计算出位移的法向梯度固定值边界条件。
Note
方程(16)与方程(20)均为边界条件，前者从物理推导而来，后者从数值推导而来。方程(16)只有在每个时间步收敛（空间项显性离散与隐性离散趋向于相等）的时候是相符的，方程(20)在每个时间步内的每次迭代都相等。二者的区别在稳态下基本无差异，在瞬态下可能会存在可见差异。更详细的请参考《无痛苦NS方程笔记》 中的OpenFOAM边界条件一节。
现在考虑压力的边界条件，对于第一类边界条件，需要给定压力的固定值。其可以从方程(14)推导而来（类似于prghPressure边界条件）：
\[
(\mathbf{T}_f-p\bfn_f)\cdot\bfn_f
=(\boldsymbol{\sigma}_f\cdot\bfn_f)\cdot\bfn_f
\]
\[
p_f=\mathbf{T}_f\cdot\bfn_f - (\boldsymbol{\sigma}_f\cdot\bfn_f)\cdot\bfn_f
\]
第二类边界条件，其更适合于自由运动的无牵引力的自由变形表面，在这种情况下，\(\mathbf{T}_f=0\)，\(\boldsymbol{\sigma}_f=0\)。有符合物理的（类似于buoyantPressure边界条件）：
\[
p\bfI
=
\mathrm{dev}\left(\boldsymbol{\sigma}\right)
\]
(21)\[
(\nabla p)_f\cdot\bfn_f
=
\left(\nabla\cdot \mathrm{dev}\left(\boldsymbol{\sigma}\right)\right)_f
\cdot\bfn_f
\]
另有更加符合迭代算法的压力梯度边界条件（类似于fixedFluxPressure边界条件）：
\[
\bfU_f=
\bfHbyA_f^{*}
+ \frac{\rho_f}{A_f}\mathrm{dev}\left(\nabla\cdot \boldsymbol{\sigma} \right)_f
-\frac{\rho_f}{A_f}(\nabla p)_f
\]
(22)\[
(\nabla p)_f\cdot\bfn_f
=
\frac
{
\bfHbyA_f^{*}
+ \frac{\rho_f}{A_f}\mathrm{dev}\left(\nabla\cdot \boldsymbol{\sigma} \right)_f
-\bfU_f
}
{\frac{\rho_f}{A_f}}
\cdot\bfn_f
\]
方程(22)在收敛情况下退化为方程(21)。
离散2
现在讨论具体的离散过程。将方程(11)移项有：
(23)\[
\underset{\mathrm{implicit}}{\underbrace{\frac{\p \rho \bfU}{\p t }-\nabla\cdot ( \dt(2\mu+\lambda)\nabla \bfU)}}
=-\nabla p
+\underset{\mathrm{explicit}}{\underbrace{ \nabla\cdot \mathrm{dev}\left(\boldsymbol{\sigma}\right)
-\nabla\cdot(\dt(2\mu+\lambda)\nabla\bfU)}}
\]
略去除掉压力之外的显性项，求解方程(23)可以获得
(24)\[
\mathbf{U}_\mathrm{P}^{*} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}
\sum_f
\left(
p_f^t\bfS_f
+ \left( \mathrm{dev}\left(\boldsymbol{\sigma}\right) + \dt(2\mu+\lambda)\nabla\bfU\right)_f^t\cdot\bfS_f
\right)
\]
其中\(\bfHbyA\)的定义可参考CFD: 不可压 + 稳态中的定义。定义：
\[
\phi_{g}^t=\left( \mathrm{dev}\left(\boldsymbol{\sigma}\right) + \dt(2\mu+\lambda)\nabla\bfU\right)_f^t\cdot\bfS_f
\]
有：
(25)\[
\mathbf{U}_\mathrm{P}^{*} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}
\sum_f
p_f^t\bfS_f
- \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}
\sum_f\phi_f^t
\]
定义新的变量：
(26)\[
\boldsymbol{\mathcal{U}}_\rP^{*}=\mathbf{HbyA}^{*}_\mathrm{P}-\frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}
\sum_f\phi_f^t
\]
将压力和密度的关系简化成正压模型：
(27)\[
\rho^{t+\dt} =\rho^0+\psi(p^{t+\dt}-p^0)
\]
其中\(\rho^0,p^0\)为参考密度参考压力。方程(13)可以重写为：
(28)\[
\frac{\p \psi p}{\p t}
+
\sum \rho_f
\left(
\boldsymbol{\mathcal{U}}_\rP^{*}
- \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}
\sum_fp_f^{t+\dt}\bfS_f
\right)\cdot\bfS_f
=0
\]
(29)\[
\frac{\p \psi p}{\p t}
+
\nabla\cdot\left( (\rho^0-\psi p^0) \boldsymbol{\mathcal{U}}_\rP^{*}\right)
+
\nabla\cdot (\psi \boldsymbol{\mathcal{U}}_\rP^{*} p)
-
\nabla\cdot\left(\frac{\rho^t}{A_\mathrm{P}^t}\nabla p\right)
=0
\]
求解后有更新的压力\(p^{*}\)。随后可以更新速度：
(30)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}
\sum_f
p_f^*\bfS_f
- \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}
\sum_f\phi_f^t
\]