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
\]
首先:
求解方程(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.