Note

本文的算法收录于LCO算法编程课中。离散过程不来自于任何公开发表的文献。在一定程度上适合作为算法编程范例在课堂上介绍。

CFD: 线弹性方程

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

\(\bfD\)算法

在不考虑体积力的情况下,线弹性固体的控制方程为:

(1)\[ \frac{\p ^2 \rho \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)}, \lambda=\frac{\nu E}{(1 + \nu)(1 - 2\nu)} \]

\(E\)表示杨氏模量,\(\lambda\)表示拉梅系数的第1个参数,\(\mu\)表示拉梅系数的第2个参数(在流体力学中理解为粘度,在固体力学中理解为剪切模量)。杨氏模量越大,越不容易发生变形。比如钢的杨氏模量大约是e11量级,普通的麻绳杨氏模量大约是e9量级,比较硬的绳子可能是e10量级。\(\nu\)是泊松比,表示当一个材料被拉伸时,横向应变与轴向应变的比值。钢的泊松比在0.2-0.3量级。木头的泊松比基本为0。同时,钢的密度大约在8000左右。如果\(\nu=0.5\),其表示的为不可压缩固体。在这种情况下\(\lambda\)趋向于无穷大因此\(\boldsymbol{\sigma}\)不可计算。解决方法是可以将\(\boldsymbol{\sigma}\)写成另外一种形式(压力分离),将在后续讨论。

将方程(2)代入到(1),有控制方程:

(4)\[ \frac{\p ^2 \rho\bfD}{\p t ^2}-\nabla\cdot\left( \mu(\nabla \bfD+\nabla \bfD^T) \right) - \nabla\cdot\left(\lambda (\nabla\cdot\bfD)\bfI \right)=0 \]

方程(4)求解的为位移矢量。因此,在给定边界条件的时候,需要给定位移矢量的边界条件。例如:边界处的固定变形是多少之类。

仔细观察方程(4),其实际上表示3个位移矢量的分量的方程。每个位移矢量分量的方程里面,涉及到其他的位移矢量分量。传统的FEM算法,是采用耦合求解的方法,离散成一个大的方程组同时求解。在CFD里面也有类似的求解算法,即为耦合算法。本文讨论分离式求解,即将三个位移矢量分量的方程离散为三个分量方程。这与速度矢量的方程离散是一样的。采用这种离散方式可以得到3个稀疏矩阵,稀疏矩阵的好处之一是可以降低内存存储,另外是可以调用迭代矩阵求解器。然而分离式求解一些变量需要做显性离散处理,因此分离式求解不仅需要用稀疏线性系统求解器来迭代求解矩阵,还需要迭代求解位移矢量分量在三个方程之间的耦合关系。

显而易见的,在求解方程(4)的时候,理想情况下,空间项应该全部隐性离散:

(5)\[ \frac{\p ^2 \rho\bfD}{\p t ^2}- \underset{\mathrm{implicit}}{\underbrace{\nabla\cdot\left( \mu(\nabla \bfD+\nabla \bfD^T) \right) - \nabla\cdot\left(\lambda (\nabla\cdot\bfD)\bfI \right)}} =0 \]

但一般将其写为可植入的形式:

(6)\[ \frac{\p ^2 \rho\bfD}{\p t ^2}- \underset{\mathrm{implicit}}{\underbrace{\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 \]

Note

具体原因请参考无痛苦NS方程笔记中的 \(\nabla\cdot((\nabla\cdot\bfU)\bfI)\)以及\(\nabla\cdot(\mu\nabla\bfU^T)\)的隐性离散 一节。

进一步的,方程(6)因为是用于描述固体而不是流体,其显性项的贡献比隐性项的贡献还要大,因此其收敛性比较差[JW00]。因此,OpenFOAM将其改写为:

(7)\[ \frac{\p ^2 \rho\bfD}{\p t ^2}- \underset{\mathrm{implicit}}{\underbrace{\nabla\cdot ( (2\mu+\lambda)\nabla \bfD)}} + \underset{\mathrm{explicit}}{\underbrace{ \nabla\cdot((2\mu+\lambda)\nabla\bfD)}} =\nabla\cdot{\bff} \]
(8)\[ \bff= \mu(\nabla\bfD+\nabla\bfD^T) + \lambda (\nabla\cdot\bfD)\bfI \]

其中\(2\mu\)中的\(2\)可以任意改变,\(2\mu+\lambda\)其实可以任意的改变,其与收敛性有关联。之所以默认为\(2\mu+\lambda\)是因为如果全部隐性离散,那么其对角线系数即为\(2\mu+\lambda\)。等号右边的\(\bff\)为必然显性离散的体积力。求解方程(7)即为OpenFOAM中solidDisplacementFoam的求解流程。

Note

方程(7)在一个时间步内要多次求解,来减少显性离散带来的滞后,这样更趋向于空间项全隐格式。

\(\bfU\)算法

首先,位移矢量与速度存在下列关系(欧拉隐性):

(9)\[ \bfD^{t+\dt}=\bfD^{t} + \int_t^{t+\dt}\bfU\rd t=\bfD^{t} +\bfU^{t+\dt}\dt \]

因此对于方程(5)中意图隐性离散的\(\nabla\cdot\left( \mu(\nabla \bfD+\nabla \bfD^T) \right) - \nabla\cdot\left(\lambda (\nabla\cdot\bfD)\bfI \right)\),有:

(10)\[\begin{split} \begin{multline} \underset{\mathrm{implicit}}{\underbrace{ -\nabla\cdot\left( \mu(\nabla \bfD+\nabla \bfD^T) \right) - \nabla\cdot\left(\lambda (\nabla\cdot\bfD)\bfI \right) }} = \\ \underset{\mathrm{implicit}}{\underbrace{ -\nabla\cdot\left( \dt\mu(\nabla \bfU+\nabla \bfU^T) \right) - \nabla\cdot\left(\dt\lambda (\nabla\cdot\bfU)\bfI \right) }} - \nabla\cdot\bff \end{multline} \end{split}\]

其中\(\bff\)的定义参考方程(8)。在这种情况下有\(\bfU\)的控制方程:

(11)\[ \frac{\p \rho\bfU}{\p t }- \underset{\mathrm{implicit}}{\underbrace{\nabla\cdot\left( \dt\mu(\nabla \bfU+\nabla \bfU^T) \right) - \nabla\cdot\left(\dt\lambda (\nabla\cdot\bfU)\bfI \right)}} =\nabla\cdot\bff \]

方程(11)(5)一样同样面临隐性离散的问题。参考\(\bfD\)的数值处理方法,有最终\(\bfU\)的控制方程:

(12)\[ \underset{\mathrm{implicit}}{\underbrace{\frac{\p \rho\bfU}{\p t }-\nabla\cdot ( \dt(2\mu+\lambda)\nabla \bfU)}} +\underset{\mathrm{explicit}}{\underbrace{ \nabla\cdot(\dt(2\mu+\lambda)\nabla\bfU)}} =\nabla\cdot\bff+\nabla\cdot\bffU \]
\[ \bffU= \dt\mu(\nabla\bfU+\nabla\bfU^T) + \dt\lambda (\nabla\cdot\bfU)\bfI \]

其中\(\bffU\)显性离散。求解方程(12)有另外一种形式的固体控制方程。

Note

\(\bffU\)在同一个时间步内可以多次更新来提高精准性,\(\bff\)在一个时间步内只能采用旧值。

\(\bfU-p\)耦合算法

另外,也可参考流体力学的方法,\(\boldsymbol{\sigma}\)可以写为另外的形式[Pap08]

(13)\[ \boldsymbol{\sigma}=\mu\left(\nabla\bfD+\nabla\bfD^T\right) -\frac{2}{3}\mu (\nabla\cdot\bfD)\bfI -p\bfI \]

方程(13)中不存在\(\lambda\),因此其适用于不可压缩固体。在这种情况下,方程(1)可以写为:

(14)\[ \frac{\p ^2 \rho\bfD}{\p t ^2} -\nabla\cdot\left( \mu(\nabla \bfD+\nabla \bfD^T) - \frac{2}{3}\mu(\nabla\cdot\bfD)\bfI\right) = -\nabla p \]

其构成\(\bfD-p\)耦合方程。同理,参考上文中\(\bfU\)方程的推导,类似的有\(\bfU-p\)耦合方程:

(15)\[ \frac{\p \rho\bfU}{\p t } -\nabla\cdot\left( \dt\mu(\nabla \bfU+\nabla \bfU^T) - \frac{2}{3}\dt\mu(\nabla\cdot\bfU)\bfI \right) = -\nabla p+\nabla\cdot\bff' \]
(16)\[ \bff'= \mu(\nabla \bfD+\nabla \bfD^T) - \frac{2}{3}\mu (\nabla\cdot\bfD)\bfI \]

如果同时调用连续性方程,以及压力与密度的关系(状态方程,在弱可压缩固体中可以简化为正压模型),则形成固体的控制方程系统:

(17)\[ \frac{\p\rho}{\p t}+\nabla\cdot(\rho\bfU)=0 \]
(18)\[ \rho =\rho^0+\psi(p-p^0) \]

其中\(\rho^0,p^0\)表示参考密度以及参考压力。方程(15)(17)(18)可用于封闭求解。

现在讨论具体的离散过程。类似方程(11)(5),方程(15)在处理离散项的时候,需要类似的特殊处理:

(19)\[ \frac{\p \rho \bfU}{\p t }- \underset{\mathrm{implicit}}{\underbrace{\nabla\cdot ( \dt(2\mu+\lambda)\nabla \bfU)}} + \underset{\mathrm{explicit}}{\underbrace{\nabla\cdot ( \dt(2\mu+\lambda)\nabla \bfU)}} =-\nabla p + \nabla\cdot\bff' + \nabla\cdot\bffU' \]
(20)\[ \bffU'= \dt\mu(\nabla \bfU+\nabla \bfU^T) - \frac{2}{3}\dt\mu (\nabla\cdot\bfU)\bfI \]

首先:

  • 显性离散体积力显性项\(\bff'\)以及\(\bffU'\)

  • 将压力显性离散;

求解方程(19)可以获得(也即动量预测)\(\mathbf{HbyA}^{*}_\mathrm{P} \)\(\bfHbyA\)的定义可参考CFD: 不可压 + 稳态中的定义。在考虑显性体积力,以及压力修正后,然后有SIMPLE/PISO循环的更新速度:

(21)\[ \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_g^t \]
(22)\[ \frac{1}{V_\rP}\sum_f\phi_g^t = \frac{1}{V_\rP}\sum_f \left( \bff'+\bffU' \right)_f^t\cdot\bfS_f \]

定义新的变量:

(23)\[ \boldsymbol{\mathcal{U}}_\rP^{*}=\mathbf{HbyA}^{*}_\mathrm{P} +\frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f\phi_g^t \]

有:

(24)\[ \mathbf{U}_\mathrm{P}^{**} = \boldsymbol{\mathcal{U}}_\rP^{*} - \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP}\sum_fp_f^t\bfS_f \]

\(\mathbf{U}_\mathrm{P}^{**} \)需要满足连续性方程。将压力和密度的关系简化成正压模型:

(25)\[ \rho^{*} =\rho^0+\psi(p^{*}-p^0) \]

其中\(\rho^0,p^0\)为参考密度参考压力。方程(17)可以重写为(\(\rho_f^*\)在每个时间步最初进行更新,在这里是一个已知量):

(26)\[ \frac{\p \psi p}{\p t} + \sum \rho_f^* \boldsymbol{\mathcal{U}}_\rP^{*}\cdot\bfS_f + \sum \rho_f^* \left( - \frac{1}{{{A^t_\mathrm{P}}}} \frac{1}{V_\rP} \sum_fp_f^{*}\bfS_f \right)\cdot\bfS_f =0 \]
(27)\[ \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^*}{A_\mathrm{P}^t}\nabla p\right) =0 \]

求解后有更新的压力\(p^{*}\)。随后可以更新速度:

(28)\[ \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_g^t \]

然后进入下一次PISO迭代。

边界条件

针对\(\bfD\)求解的情况下,需要给定\(\bfD\)的边界条件。\(\bfD\)如果给定固定值边界条件,则假定其固定位移,在这种情况下表面的应力可能会非均一分布。对于任意的固体表面,可给定牵引力\(\mathbf{T}_f\)的边界条件。牵引力是定义在面上的一个力,存在三个方向。给定具体的面,三个分量分别对应\(x,y,z\)方向的牵引力。除了牵引力外,还存在用户施加的法向压力\(p^{udf}_f\)。其作用于面的法向,因此是一个标量。牵引力和压力都会导致物体发生形变。压力\(p_f^{udf}\bfn\)(正值向内)与牵引力\(\mathbf{T}_f\)反向。

方程(7)可以用来构建边界条件,但其中并没有\(\bfD\)以及\(\nabla\bfD\)的存在,因此从方程(7)中构建第一二类边界条件并不是显而易见的。

另一方面,可以从物理的角度来构建边界条件,对于给定表面,在任意的情况下均存在下述力平衡:

(29)\[ \mathbf{T}_f-p_f^{udf}\bfn_f=\boldsymbol{\sigma}_f\cdot\bfn_f \]

通过\(\boldsymbol{\sigma}_f\)的定义,其可以用来构建\(\bfD\)的边界条件:

(30)\[ \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 \]

有:

(31)\[ (\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} \]

方程(7)的另一种理解方法即为\(\boldsymbol{\sigma}\)分解为:

(32)\[ \boldsymbol{\sigma}=\underset{\mathrm{implicit}}{\underbrace{ (2\mu+\lambda)\nabla \bfD}} +\boldsymbol{\sigma}- \underset{\mathrm{explicit}}{\underbrace{(2\mu+\lambda)\nabla\bfD }} \]

将方程(32)代入到方程(29)

(33)\[ \mathbf{T}_f-p^{udf}_f\bfn_f= \left(\underset{\mathrm{implicit}}{\underbrace{ (2\mu+\lambda)\nabla \bfD}} +\boldsymbol{\sigma}- \underset{\mathrm{explicit}}{\underbrace{ (2\mu+\lambda)\nabla\bfD}}\right)_f\cdot\bfn_f \]
(34)\[ \mathbf{T}_f-p_f^{udf}\bfn_f= \underset{\mathrm{implicit}}{\underbrace{ (2\mu+\lambda)(\nabla \bfD)_f\cdot\bfn_f}} +\boldsymbol{\sigma}_f\cdot\bfn_f- \underset{\mathrm{explicit}}{\underbrace{(2\mu+\lambda)(\nabla\bfD)_f\cdot\bfn_f }} \]

移项有:

(35)\[ \underset{\mathrm{implicit}}{\underbrace{(\nabla\bfD)_f \cdot\bfn_f}} = \frac{ \mathbf{T}_f-p_f^{udf}\bfn_f- \boldsymbol{\sigma}_f \cdot\bfn_f +(2\mu+\lambda) \underset{\mathrm{explicit}}{\underbrace{(\nabla\bfD)_f \cdot\bfn_f}} }{2\mu+\lambda} \]

从方程(35)可以看出,给定牵引力\(\mathbf{T}_f\)的时候,\(\bfD\)会形成一个法向梯度固定值边界条件,其值通过方程(35)计算而来。对于自由表面,牵引力\(\mathbf{T}_f\)\((0, 0, 0)\)。对于固定牵引力边界,需要给定牵引力\(\mathbf{T}_f\)的值,然后方程(35)会自动计算出位移的法向梯度固定值边界条件。

Note

方程(31)与方程(35)均为边界条件,前者从物理推导而来,后者从数值推导而来。方程(31)只有在每个时间步收敛(空间项显性离散与隐性离散趋向于相等)的时候是相符的,方程(35)在每个时间步内的每次迭代都相等。二者的区别在稳态下基本无差异,在瞬态下可能会存在可见差异。更详细的请参考《无痛苦NS方程笔记》 中的OpenFOAM边界条件一节。

同理,对于\(\bfU\),方程(35)可以用于推导\(\bfU\)的边界条件:

(36)\[ \underset{\mathrm{implicit}}{\underbrace{(\nabla\bfU)_f \cdot\bfn_f}} = \frac { \mathbf{T}_f-p_f^{udf}\bfn_f- \boldsymbol{\sigma}_f \cdot\bfn_f } {\dt(2\mu+\lambda)} + (\nabla\bfD)_f \cdot\bfn_f \]

现在考虑压力的边界条件,对于第一类边界条件,需要给定压力的固定值。其可以从方程(13)推导而来

(37)\[ p\bfn_f = \left(\mu\left(\nabla\bfD+\nabla\bfD^T\right) -\frac{2}{3}\mu (\nabla\cdot\bfD)\bfI\right)_f\cdot\bfn_f -\mathbf{T}_f+p_f^{udf}\bfn_f \]

左右两边同时乘以\(\bfn_f\)有:

(38)\[ p = \bfn_f\cdot\left(\left(\mu\left(\nabla\bfD+\nabla\bfD^T\right) -\frac{2}{3}\mu (\nabla\cdot\bfD)\bfI\right)_f\cdot\bfn_f\right) -\bfn_f\cdot\mathbf{T}_f+p_f^{udf} \]

除了第一类边界条件,还存在第二类边界条件,依据符合迭代算法的压力梯度边界条件(类似于fixedFluxPressure边界条件),依据方程(21)

\[ \bfU_f= \bfHbyA_f^{*} + \frac{1}{A_f} \nabla\cdot\left(\bff'+\bffU' \right)_f -\frac{1}{A_f}(\nabla p)_f \]
(39)\[ (\nabla p)_f\cdot\bfn_f = \frac { \bfHbyA_f^{*} + \frac{1}{A_f} \nabla\cdot\left(\bff'+\bffU' \right)_f -\bfU_f } {\frac{1}{A_f}} \cdot\bfn_f \]
[JW00]

Hrvoje Jasak and HG Weller. Application of the finite volume method and unstructured meshes to linear elasticity. International journal for numerical methods in engineering, 48(2):267–287, 2000.

[Pap08]

George Papadakis. A novel pressure–velocity formulation and solution method for fluid–structure interaction problems. Journal of computational physics, 227(6):3383–3404, 2008.