FEM-Summary

FEM simulation of 3D Deformable Solids. 这里主要是对FEM在三维物体变形的一个简单的介绍。

三维空间里的弹性物体

在这个章节主要关注三维空间里的弹性物体如何变形,并讨论如何对变形物体和力进行数学量化表示。

变形映射和变形梯度

首先,描述一个物体要放在一个特定的坐标系中。用$\Omega$表示物体,未变形状态称为reference configuration(或者undeformed configuration)。
用$\overrightarrow{X}\in \Omega$来表示物体$\Omega$的未变形状态时的一个顶点。当物体$\Omega$受力变形之后,每一个$\overrightarrow{X}$都会发生位移到新的位置$\overrightarrow{x}$。用变形映射来表示物体变形前后之间的关系$\overrightarrow{\phi}:R^{3}\rightarrow R^{3}$,它将每个点$\overrightarrow{X}$映射到对应的变形位置$\overrightarrow{x}=\overrightarrow{\phi}(\overrightarrow{x})$。

deformation map

从变形映射直接得出的一个重要的物体量就是变形梯度张量$F\in R^{3\times 3}$。如果$\overrightarrow{X}=(X_{1},X_{2},X_{3})^{T}$,$\overrightarrow{\phi}(\overrightarrow{X})=(\phi_{1}(\overrightarrow{X}),\phi_{2}(\overrightarrow{X}),\phi_{3}(\overrightarrow{X}))^{T}$,那么变形梯度就为
$$F = \frac{\partial{(\phi_{1},\phi_{2},\phi_{3})}}{\partial{(X_{1},X_{2},X_{3})}}=
\begin{pmatrix}
\partial\phi_{1}/\partial{X_{1}} & \partial\phi_{1}/\partial{X_{2}} &\partial\phi_{1}/\partial{X_{3}}\\
\partial\phi_{2}/\partial{X_{1}} & \partial\phi_{2}/\partial{X_{2}} & \partial\phi_{2}/\partial{X_{3}}\\
\partial\phi_{3}/\partial{X_{1}} & \partial\phi_{3}/\partial{X_{2}} & \partial\phi_{3}/\partial{X_{3}}
\end{pmatrix}
$$
或者用下标标记法 $F_{ij}=\phi_{i,j}$。从上面就知道 $F$其实就是变形映射Jacobian矩阵。需要注意的是,通常情况下$F$会随着时间变化。

对变形梯度的简单解释
Figure 2.1
如上图(a),表示的是物体在空间只进行平移,形状并没有任何改变。
$$\overrightarrow{x}=\phi_{(\overrightarrow{X})}=\overrightarrow{X}+\overrightarrow{t},\quad F=\partial\phi_{\overrightarrow{X}}/\partial{\overrightarrow{X}}=I$$
如上图(b),表示的是物体扩大了$\gamma$倍。
$$\phi_{\overrightarrow{X}}=\gamma\overrightarrow{X},\quad F = \gamma I$$
如上图(c),物体在水平方向缩小0.7倍,垂直方向放大了2倍。
$$\begin{pmatrix}x\\ y \end{pmatrix}=\phi\begin{pmatrix}X\\ Y \end{pmatrix}=\begin{pmatrix}0.7X\\ 2Y \end{pmatrix}, \quad F=\begin{pmatrix}0.7 & 0\\ 0 & 2 \end{pmatrix}$$
如上图(d),物体关于原点逆时针旋转了$45^{0}$
$$
\begin{pmatrix}
x\\
y
\end{pmatrix}=
\begin{pmatrix}
cos45^{0} & -sin45^{0}\\
sin45^{0} & cos45^{0}
\end{pmatrix}
\begin{pmatrix}
X\\
Y
\end{pmatrix}, \quad \quad
F= \begin{pmatrix}
cos45^{0} & -sin45^{0}\\
sin45^{0} & cos45^{0}
\end{pmatrix}$$

应变能和超弹性

一系列的变形之后,就会导致变形物体势能(potential energy)的积累,在变形体上下文中也即是应变能(strain energy)。用$E[\phi]$来表示应变能,也就说应变能完全由给定状态的变形映射决定。这里有个假设:物体变形产生的应变能只与变形初始和结束状态有关,与变形过程中的路径无关。也就说,不管你中间怎么着,我不管;我只看起点与终点物体的状态。这种与之前的变形历史无关的特性是超弹性(hyperelastic)材料的特性。这种特性与超弹性材料弹性力的保守性紧密相关。

保守性:变形路径上由弹性内力做的功只依赖于初始状态和最终状态,而与路径无关。这与高中物理重力做功形式,重力对物体做的功只与物体的初始位置与最终位置有关,而与物体中间行走的路径无关。所以重力场是一种保守场。

受力变形的物体不同部分的变形程度是不同的。所以,变形与应变能之间的关系应该定义在局部区域上。这也就引入了能量密度函数(energy density function) $\Psi[\phi;\overrightarrow{X}]$,它表示的是顶点$\overrightarrow{X}$附近无穷小区域$dV$上每个单元未变形体积(per unit undeformed volume)的应变能。那么$\Omega$域上的应变能可以通过能量密度函数积分得到:
$$E[\phi]=\int _{\Omega} \Psi[\phi;\overrightarrow{X}]d\overrightarrow{X}$$

在特定的位置$\overrightarrow{X}_{*}$,因为$\Psi[\phi;\overrightarrow{X}_{*}]$反应的是$\overrightarrow{X}_{*}$附近无穷小区域的变形。所以我们可以用一阶泰勒扩展式来估计这一小块区域的变形映射。
$$\phi(\overrightarrow{X})\approx \frac{\partial{\phi}}{\overrightarrow{X}_{*}}\mid_{\overrightarrow{X}_{*}}(\overrightarrow{X}-\overrightarrow{X}_{*}) = \overrightarrow{x}_{*} + F(\overrightarrow{X}_{*})(\overrightarrow{X}-\overrightarrow{X}_{*}) =F_{*}\overrightarrow{X}+\overrightarrow{t}$$
其中$F_{*}=F(\overrightarrow {X}_{*}) * X $, $\overrightarrow{t}=\overrightarrow{x}_{*}-F(\overrightarrow{X}_{*}) \overrightarrow{X}_{*}$
这个方程表明$\Psi[\phi;\overrightarrow{X}_{*}]$应该是$F_{*}$和$\overrightarrow{t}$的函数。更进一步,我们可以认为向量$\overrightarrow{t}$的值在这个式子中是无关的:这个参数的不同值表示的是不同的常量位移,因此它产生的效果不改变物体形状和应变能。

综上:应变能密度函数可以表示为:
$$\Psi[\phi;\overrightarrow{X}_{*}]=\Psi(F(\overrightarrow{X}_{*}))$$
它只与局部变形梯度有关。

但是,现在我们并没有给出应变能密度函数的具体表达形式。它与具体的材料属性有关,后面会讲不同材料的应变能密度函数怎么求,并推出相应材料的具体应变能密度函数的表达式。

力和牵引力

接下来要讲的一个物理量是物体发生形变时产生的弹性力
First:
举个栗子:一个质点(只有质量,没有体积)在保守场内运动,假设保守场就是重力场。质点在$\overrightarrow{x}=(x,y,z)$有势能$E(\overrightarrow(x))=m\overrightarrow{G}\overrightarrow{x}=mgz,\quad \quad \overrightarrow{G}=(0,0,g)$,是重力加速度($z-axis$为垂直方向)。此时可以得到重力就是质点在位置$\overrightarrow{x}$处势能梯度的负方向
$$\overrightarrow{f} =-\frac{\partial E(\overrightarrow(x))}{\partial \overrightarrow{x}}=(0,0,-mg)$$

Then:
我们也想可变性物体中力和能量之间这种关系,但是在对变形物体离散化之前,我们认为它是连续分布的,而不是一系列独立的质点的集合。所以使用内力密度$\overrightarrow{f}(\overrightarrow{X})$表示点$\overrightarrow{X}$附近无穷小区域内每单元体积受到的(force per unit undeformed volume)。此时,有限区域$A\in \Omega$上的力的总和可以由积分获得。
$$\overrightarrow{f}_{aggregate}(A)=\int _{A} \overrightarrow{f}(\overrightarrow{X})d\overrightarrow{X}$$

However:
上述方法并不适合表示施加在可变性物体边界($\partial \Omega$)上的力。假设一个弹性体被挤压(e.g. $\phi(\overrightarrow{X}=\alpha \overrightarrow{X},\quad \alpha<1$)。那么弹性体表面$S\subset \partial \Omega$应该会有一个恢复力(restorative force)来反推造成挤压的物体。同样,定义牵引力(traction)$\overrightarrow\tau (\overrightarrow{X})$为表面力密度函数。用它来表示物体表面上$\overrightarrow{X}$附近无穷小区域单位面积力(force per unit undeformed area)。那么,有限表面区域$B\subset \partial \Omega$上力的总和可以由积分获得。
$$\overrightarrow{f}_{aggregate}(B)=\oint _{B} \overrightarrow\tau (\overrightarrow{X})dS$$

The First Piola-Kirchhoff 应力张量

内力密度与表面牵引力的不同表明:他们中任何一个都不能很好的描述弹性体所有方面的弹性响应性质。但是,这两种力都可以引出一种基本的力描述符:应力张量(Stress Tensor)。其实,有很多种”应力”描述符。但我们只讨论 1st Piola-Kirchhoff stress tensor $P$,它是一个$3 \times 3$的矩阵,它有如下性质:

  • 表面位置$\overrightarrow{X} \in \partial\Omega$的内部牵引力为
    $$\overrightarrow\tau (\overrightarrow{X}) = -P \cdot \overrightarrow{N}$$
    $\overrightarrow{N}$是未变形状态表面单位法向,方向向外。
    由上可得应力张量$P$的定义:对于变形体内任一点$\overrightarrow{X} \in \Omega \setminus \partial\Omega$,可以假设将物体一刀切开,切的轨迹经过$\overrightarrow{X}$并且与$\overrightarrow{N}$相互垂直。那么,$P$就是一个与$\overrightarrow\tau,\overrightarrow{N}$有关如上式的一个矩阵。
  • 由$P$可以得出内力密度函数
    $$\overrightarrow{f}(\overrightarrow{X})=div_{\overrightarrow{X}}P(\overrightarrow{X})$$
    或者下标形式$$f_{i}=\sum_{j=1}^{3}P_{ij,j}=\frac{\partial P_{i1} }{\partial X_{1}}+\frac{\partial P_{i2} }{\partial X_{2}}+\frac{\partial P_{i3} }{\partial X_{3}}$$
  • 对于超弹性材料,$P$仅仅是关于变形梯度的函数,并且与应变能之间的关系很简单:
    $$P(F)=\frac {\partial \Psi(F)}{\partial F}$$

如上所述,由1st Piola-Kirchhoff stress tensor可以表示内力(force)和张力(tension)。而且可以很简单从应变能密度函数求得。

有两种(等价)方法表示超弹性材料的材料属性:
(a) 显示方法:以$F$为参数的$\Psi$
(b) 隐式方法:以$F$为参数的$P$

下面谈论的所有材料,我们都会给出两种这两种方式的表示方法。


材料的本构模型

在一个部分,主要介绍在仿真领域内比较流行的一些材料,并展示如何表示各种材料的物理属性,本构模型(给定材料物理属性的数学描述)和特定刺激(e.g. 变形)下材料的反应方程(比如:力,应力,能量)。
如前面所讲,每个本构方程有两种表现形式:
(a) Piola 应力 $P$ 关于变形梯度$F$的函数的表示形式
(b) 能量密度函数$\Psi$ 关于变形梯度$F$ 的函数的表示形式
为了简单,我们讨论各向同性(isotropic)材料。

各种同性材料:弹性物体对于变形的反应与变形的方向无关。

应变度量(Strain measures)

原则上来说,直接用$F$和$\Psi$(或$F$和$P$)表示本构方程就足够了。但是呢,这种方式并不直观。所以引入了从$F$得出的一个中间物理属性,比如:应变度量(strain measures)或者不变量(invariants)。

应变度量是对变形程度的定量描述符。所以,虽然应变度量是从变形梯度得来,但是还能获得变形程度的信息。
Green Strain tensor $E \in R^{3\times 3}$:
$$E=\frac{1}{2}(F^{T}F-I)$$

当弹性体为未变形状态的时候:比如$\overrightarrow{\phi}(\overrightarrow{X})=\overrightarrow{X}\quad \rightarrow \quad F=I \rightarrow \quad E=0$。
当弹性物体仅仅旋转和平移的时候:$\overrightarrow{\phi}(\overrightarrow{X})=R\overrightarrow{X}+\overrightarrow{t}\quad$(R为旋转矩阵)。那么$F=R\quad \rightarrow E=0$因为$R^{T}R=I$
一般情况:变形梯度可以进行极分解,$F=RS$,$R$为旋转矩阵,$S$为半正定对称矩阵。三维旋转阵$R$包含3个自由度,$S$有6个独立自由度。此时
$$E=\frac{1}{2}[(RS)^{T}(RS)-I]=\frac{1}{2}(S^{2}-I)$$
上面的式子是非线性(二次)函数,增加了本构模型的复杂度,也会导致离散化节点力(nodal forces)和节点位置(nodal positions)时候的非线性问题。为了解决非线性问题,通过泰勒展开式在$F+I$处,对$E$进行了线性化估计。
$$E(F)\approx E(I)+\frac{\partial E}{\partial F} \mid_{F=I} :(F-I)$$
其中$E(I)=0$,用$\delta E$计算微分$\partial E / \partial F$
$$\frac{\partial E}{\partial F}:\delta F = \delta E = \frac{1}{2}(\delta F^{T}+F^{T}\delta F)$$
所以在$F=I$时:
$$\frac{\partial E}{\partial F} \mid_{F=I} :(F-I)=\frac{1}{2}[(F-I)^{T}I-I^{T}(F-I)]=\frac{1}{2}(F+F^{T})-I$$
由此,得到对$E(F)$的线性估计:小应变张量(small strain tensor)或者无穷下应变张量(infinitesimal strain tensor)表示为$\varepsilon$:
$$\varepsilon = \frac{1}{2}(F+F^{T})-I$$

优势: 离散化后的节点力和节点位置之间的映射是线性的。
局限:只能应用于小变形,不能用于表示大变形。

线性弹性(Linear Elasticity)

最简单的本构模型就是线弹性(Linear Elasticity)模型,对应的应变能密度函数为:
$$\Psi(F) = \mu \varepsilon: \varepsilon + \frac{\lambda}{2}tr^{2}(\varepsilon)$$
其中$\varepsilon$为小应变张量,$\mu$和$\lambda$为$Lam\acute{e} coefficients$,都是材料的属性Young's modulus$k$和Poisson's ratios$v$有关。

$$\mu = \frac{k}{2(1+v)},\quad \quad \lambda = \frac{kv}{(1+v)(1-2v)}$$

the Piola stress $P$ 和 $F$ 之间的关系:
$$P(F) = \mu(F+F^{T}-2I)+ \lambda tr(F-I)I$$

该模型的缺点:只能用于对小形变的建模,不使用于大型变。在有大形变的情况下,会产生体积膨胀的现象。

St. Venant-Kirchhoff model

不再使用对Green strain的线性估计进行计算,而是使用非线性Green-Lagrange应变$E=\frac {1}{2}(F^{T}F-I)$,这样做可以消除大旋转(large rotation)带来的缺陷。
$$\Psi(F) = \mu E: E + \frac{\lambda}{2}tr^{2}(E)\\
P(F)=F[2\mu E + \lambda tr(E) I]$$
StVK的弹性能量函数是四次(4阶)多项式。因此,内力和刚度矩阵分别是三次和二次多项式。

StVK可能是最简单的非线性模型了,因为它使用一个线性映射来对stress-strain关系进行建模。不像一般的线性材料可能产生的各向异性(anisotropic)问题,StVK是各向同性(isotropic)的,所以也称StVK为几何各向同性的非线性材料模型(isotropic geometrically nonlinear material model)。

优点

  • 积分可以在预处理时进行
  • 多项式系数的存储只需要少量存储空间
  • 内力是三次多项式,可以计算出内力之后求出其一阶导数(刚度矩阵)和二阶导数(Hessian of the internal forces)。这些高阶导数在优化和控制变形过程中很重要
  • 三次的内力在model reduction的时候很有用

缺点

  • 在大挤压情况下,材料会塌陷甚至翻转
  • 耗时

Corotated Linear Elasticity

极分解之后$F=RS$,由此得出的新的应变度量为$\varepsilon_{c} = S-I$,只是线性依赖于对称张量$S$,而与旋转成分$F$无关。
此时:
$$\Psi(F) = \mu \varepsilon_{c}: \varepsilon_{c} + \frac{\lambda}{2}tr^{2}(E)=\mu \left | S-I \right | _{F}^{2} +(\lambda/2)tr^{2}(S-I)$$

$$P(F)=F[2\mu \varepsilon_{c} + \lambda tr(\varepsilon_{c}) I] = 2\mu (F-R)+\lambda tr(R^{T}F-I)R$$

各向同性材料和不变量

Neohookean Elasticity

文章目录
  1. 1. 三维空间里的弹性物体
    1. 1.1. 变形映射和变形梯度
    2. 1.2. 应变能和超弹性
    3. 1.3. 力和牵引力
    4. 1.4. The First Piola-Kirchhoff 应力张量
  2. 2. 材料的本构模型
    1. 2.1. 应变度量(Strain measures)
    2. 2.2. 线性弹性(Linear Elasticity)
    3. 2.3. St. Venant-Kirchhoff model
    4. 2.4. Corotated Linear Elasticity
    5. 2.5. 各向同性材料和不变量
    6. 2.6. Neohookean Elasticity
,