CFD: 线弹性方程

引言

大部分的固体应力分析都使用的是有限元。然而使用有限体积法计算应力应变以及位移也是可行的。同时有限体积法具有一些优良的优点：如守恒特性。在将力离散到网格单元面的时候，所有的力均被抵消仅仅在边界上存在受力。这与与流体的质量守恒类似。OpenFOAM中的solidDisplacementFoam则采用有限体积法来对位移矢量\(\bfD\)方程进行离散，进而获得体对称应变张量\(\boldsymbol{\sigma}\)。本文对其简述。注意本文的方程只适用于小变形。

控制方程与算法

在不考虑体积力的情况下，线弹性固体的控制方程，也即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个参数（在流体力学中理解为粘度，在固体力学中理解为剪切模量）。将方程(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}\]