控制方程与算法
在不考虑体积力的情况下，线弹性固体的控制方程，也即solidDisplacementFoam中求解方程为：
(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)求解的为位移矢量。因此，在给定边界条件的时候，需要给定位移矢量的边界条件。例如：边界处的固定变形是多少之类。
离散
仔细观察方程(5)，其实际上表示3个位移矢量的分量的方程。每个位移矢量分量的方程里面，涉及到其他的位移矢量分量。传统的FEM算法，是采用耦合求解的方法，离散成一个大的方程组同时求解。在CFD里面也有类似的求解算法，即为耦合算法。
本文讨论分离式求解，即将三个位移矢量分量的方程离散为三个分量方程。这与速度矢量的方程离散是一样的。采用这种离散方式可以得到3个稀疏矩阵，稀疏矩阵的好处之一是可以降低内存存储，另外是可以调用迭代矩阵求解器。然而分离式求解一些变量需要做显性离散处理，因此分离式求解不仅需要用稀疏线性系统求解器来迭代求解矩阵，还需要迭代求解位移矢量分量在三个方程之间的耦合关系。
显而易见的，在求解方程(5)的时候，其中显性与隐性离散的处理关系可以写为：
(6)\[\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
\]
然而在实际过程中，方程(6)的离散过程仅仅会部分收敛。这是因为显性项的贡献比隐性项的贡献还要大。因此，OpenFOAM将其改写为：
(7)\[\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将进一步化简：
(8)\[
\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\)其实可以任意的改变，其与收敛性有关联。
针对\(\bfD\)求解的情况下，需要给定\(\bfD\)的边界条件。\(\bfD\)如果给定固定值边界条件，则假定其固定位移，在这种情况下表面的应力可能会非均一分布。对于任意的固体表面，可给定牵引力\(\mathbf{T}\)的边界条件，如果其存在值，则表明存在一个牵引力来导致形变。如果不存在值，则表明不存在牵引力导致发生形变。对于给定表面，牵引力\(\mathbf{T}\)的定义为：
\[
\mathbf{T}=\boldsymbol{\sigma}\cdot\bfn,\; \bfn=\frac{\bfS_f}{|\bfS_f|}
\]
依据方程(2)，有：
\[
\mathbf{T}=
\left( (2\mu+\lambda)\nabla\bfD +\boldsymbol{\sigma} -(2\mu+\lambda)\nabla\bfD\right)\cdot \bfn
\]
整理后有：
\[
\mathbf{T}
-\left( \boldsymbol{\sigma} -(2\mu+\lambda)\nabla\bfD\right)\cdot\bfn
=
\left( (2\mu+\lambda)\nabla\bfD \right)\cdot\bfn
\]
\[
\frac{
\mathbf{T}
-\left( \boldsymbol{\sigma} -(2\mu+\lambda)\nabla\bfD\right)\cdot\bfn
}{2\mu+\lambda}
=
\nabla\bfD \cdot\bfn
\]
(9)\[
\nabla\bfD \cdot\bfn
=
\frac{
\mathbf{T}
+(2\mu+\lambda)\nabla\bfD \cdot\bfn
- \boldsymbol{\sigma} \cdot\bfn
}{2\mu+\lambda}
\]
从方程(9)可以看出，给定牵引力\(\mathbf{T}\)的时候，\(\bfD\)会形成一个法向梯度固定值边界条件，其值通过方程(9)计算而来。对于自由表面，牵引力\(\mathbf{T}\)为\((0, 0, 0)\)。对于固定牵引力边界，需要给定牵引力\(\mathbf{T}\)的值，然后方程(9)会自动计算出位移的法向梯度固定值边界条件。
在存在法向压力\(p\)的情况下，压力\(p\cdot\bfn\)（正值向内）与牵引力\(\mathbf{T}\)反向（正值向外），因此\(\bfD\)的边界条件变更为：
(10)\[
\nabla\bfD \cdot\bfn
=
\frac{
(\mathbf{T} - p\cdot\bfn)
+(2\mu+\lambda)\nabla\bfD \cdot\bfn
- \boldsymbol{\sigma} \cdot\bfn
}{2\mu+\lambda}
\]