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 \]