经典网文分享: FEM的数值精度与高阶单元

之前听过清华FEM的Mooc前九章，知道了FEM的最基本的概念，介绍了FEM的流程，简单的单元研究，但实际上开始用ABAQUS操作以后，发现在ABAQUS中，常常采用的单元，比C3D20R这种三维二阶减缩积分单元，分析结果是摸不清楚的，看了庄茁老师讲的实体单元，有了一些头绪，大致懂了什么叫剪切自锁，什么叫沙漏，但也并没有给出一些更具体的数学表达，于是回头看看曾老师讲的第十和十一章，有关数值问题的精度和高阶单元的介绍，补完这一系列的课程，再加一些相关的笔记，顺便推荐一本书《ABAQUS常见问题解答》。

在正式写笔记之前，先补一下参数单元的知识：

参数单元就是长得不规则的单元，因为真实物理性状复杂，不可能每个都长得四四方方。对于不同的参数单元，希望从规则单元通过变换来得到参数单元，以2D问题为例，首先引入两个坐标系，物理系（x,y）和基准系 $（\xi，\eta）$ ，基准系是建立在单元上的（可变），而物理系是建立在空间上的（不变），将建立在不同单元上的基准系得到的变量映射到建立在空间上的系中，可以增加不同形状单元得到变量的通用性，下面分析需要哪些变换：

由此可见，需要变换的有：

<1>坐标；<2>偏导数；<3>积分变量（面积微元）；

下面首先来说<1>坐标的映射：

学过斜角直线坐标的张量分析，就应该对坐标变换非常熟悉了，现在空间中有一个位置P，可以用一个向量P来表示，这个向量在不同的坐标系中有不同的坐标，即：

$P=x_{1}^{（1）}e_{1}^{（1）}+x_{2}^{（1）}e_{2}^{（1）}=x_{1}^{（2）}e_{1}^{（2）}+x_{2}^{（2）}e_{2}^{（2）}$

并且分量之间一定存在唯一确定的转换关系：

$x_{1}^{（1）}=x_{1}^{（1）}（x_{1}^{（2）}，x_{2}^{（2）}）$

$x_{2}^{（1）}=x_{2}^{（1）}（x_{1}^{（2）}，x_{2}^{（2）}）$

首先，这些都是理论上行的通的，然而，在这里我们不用精确的表达式算法，而用插值的方式来获得每一个基准系上的点在物理系中的坐标（以2D4Node为例）：

回忆一下形函数的构造，根据自由度数确定待定系数，插值函数从Pascal三角形中获得，这里是完全相同的道理，所以选取的坐标插值函数为：

接下来，通过四个节点处的八个坐标值来确定坐标插值函数中的待定系数：

其实仔细考察，便会发现，这和2D4Node的形函数是完全一致的，只不过一个是通过节点位移插值出位移场，一个是在坐标变换中通过节点坐标插值出整个场的坐标。

最后关于坐标变换，给出整个场的变换表达式：

关于<2>的偏导数变换，相信学过高等数学的人都会，引入雅各比矩阵，直接上图不费话了：

最后是<3>积分变量的变换，或者说面积微元的变换，这里需要用到一些行列式和外积的几何性质，以及微分几何的一点内容：

第一个表达式，可以从几何意义入手，叉乘的大小代表平行四边形的面积，这个不是绝对值，也不是行列式，是给一个向量取模。那么具体应该如何实施呢，通过将两个向量表达成以i和j为基地的空间中，再进行外积，我写一下我的理解，从几何角度分析的（如果不对请指正）：

$d\bar{\xi}=（d\bar{\xi}）_{x}\bar{i}+（d\bar{\xi}）_{y}\bar{j} =（d\bar{\xi}·\bar{i} ）\bar{i}+（d\bar{\xi}·\bar{j} ）\bar{j}$

$=d\xi·cos<d\bar{\xi},\bar{i}>\bar{i}+d\xi·cos<d\bar{\xi},\bar{j}>\bar{j}$

$=d\xi·\frac{\triangle x}{\triangle \xi} ·\bar{i}+d\xi·\frac{\triangle y}{\triangle \xi} ·\bar{j}$

$=d\xi·\frac{\partial x}{\partial \xi} ·\bar{i}+d\xi·\frac{\partial y}{\partial\xi} ·\bar{j}$

$dA=\left| d\bar{\xi}\times d\bar{\eta} \right|= \left| det \left( \begin{array}{ccc} \bar{i}& \bar{j} & \bar{k} \\ \frac{\partial x}{\partial \xi}· d \xi & \frac{\partial y}{\partial \xi}· d \xi & 0 \\ \frac{\partial x}{\partial \eta}· d \eta & \frac{\partial y}{\partial \eta}· d \eta & 0 \\ \end{array} \right) \right| = \left| J \right| d \xi d\eta$

到这里，三个代换就做完了，这里在第一个坐标变换中引入了坐标的插值，根据坐标的插值和位移的插值的相对大小，参数单元分为：

下次的补充内容再说高斯积分，下面是正式的课程记录笔记。

（1）先说第十章的前三节，和文章题目关系不大，简单做一下笔记。分别是讲了带宽（见下图），形函数与刚度矩阵的性质，边界条件的处理，研究总体刚度矩阵的带宽是为了节省计算机的存储空间的，我目前并不怎么关系怎么存储；边界条件的处理，是讲除了我们在矩阵位移法里面对于零边界划掉行列以外，为了保持原有刚度矩阵各行各列的编号，采用了一些利于计算机编程的处理思路，比如置一法，乘大数法，罚函数法，并且他们可以处理非零边界，这个在有限元软件上是已经实现过的，我目前不自己编程，也不怎么需要了解；另外的就是刚度矩阵的性质，其实考试前早都背过了，再写一写，先说单元刚度矩阵，性质有一下几条：Kii代表在i处产生单位一的位移时，在i节点所施加的力的大小；Kij（非对角线元素）表示在j处产生单位一的位移，在i处所施加的力的大小；害怕第二条记混淆吧，没关系，第三条告诉我们们Kij=Kji，即刚度矩阵是对称阵；然而这个矩阵也是奇异的，从求解刚度矩阵前需要带入边界条件就可以得知，刚度矩阵这个平衡方程是允许产生刚体位移的；既然允许产生刚体位移，这个矩阵就不能是正定的，因为产生的如果是刚体位移，其应变能应该依然为零，可见，刚度矩阵具有半正定的特征；最后回到刚度矩阵的由来，推导它通过的是最小势能原理，最小是能原理是平衡方程在变分条件下转化而来的，刚度矩阵因此具有平衡的特性，在每一行对刚度矩阵的元素求和，结果均为零，也不难通过这个性质得到它的行列式为零，又反推出奇异性。另外还有形函数的性质：Ni的含义是i处产生单位一的位移，其他位置不产生位移时整个的位移场；对Ni进行求和，结果为常数一，感觉像是权函数对吧，实际上是通过带入刚体位移得到的这个性质。

（2）第四节重点来了，关于收敛性性讨论。收敛的意思是指"单元的尺寸趋向于0时，FEM所给出的解趋向于真实解，不仅要讨论单元内部，还要讨论单元与单元之间"，另一方面，我们应该尽量使得结果单调收敛。

思路很明确，我们从收敛这个结果，一步一步的推出它的必要条件，第一个是单元尺寸趋向于零时，势能的表达式依然存在，首先对势能的表达式进行展开：

$\Pi\left( u \right)=\Sigma（\frac{1}{2}\int_{\Omega}^{} \varepsilon^{T} D \varepsilon d \Omega - \int_{\Omega}^{} f u d \Omega - \int_{S}^{} \bar{f}u d S）$

$\Pi\left( u \right)=\Sigma（\frac{1}{2}\int_{\Omega}^{} (\varepsilon^{T}_{0} D \varepsilon_{0}+......) d \Omega - \int_{\Omega}^{} (f u_{0}+......) d \Omega - \int_{S}^{} （\bar{f}u_{0}+......） d S）$

位移函数必须连续，当体积趋向于零时，省略的高阶部分都趋向于零，为了使得势能存在，应该存在常应变项与常位移项，另外，位移函数在单元之间保持连续，才能保证应变的存在。

然后就有两个收敛准则要给出来：

<1>完备性要求：如果在势能泛函中所出现位移函数的最高阶导数为m阶，则单元内部位移场函数应该为m阶完备的多项式；

<2>协调性要求：如果在势能泛函中所出现位移函数的最高阶导数为m阶，则单元边界位移场应该具有m-1阶连续导数；

其实，之前在我们确定位移模式的时候已经默认了这些要求，只不过没有具体的提出来，现在给出位移模式中函数的选取规则：

<1>唯一确定性原则：节点的自由度个数等于待定系数个数，从低阶到高阶进行选取，选取完备的多项式提高精度（参考帕斯卡三角形）；

<2>单元内部的完备性：包含常数项（常位移）和完备的一次项（常应变）；

<3>单元之间的协调性：对于 $C_{1}$ 连续性单元（存在一阶连续导数）较难保证；

（3）第五节，第六节，关于 $C_{0}$ 单元与 $C_{1}$ 单元，说一下我的理解，这里说的 $C_{i}$ 型单元，是针对协调性要求来讲的，由模型的势能泛函决定。对于平面问题和空间问题，势能泛函中出现了应变，而应变为位移的一阶导数，所以只要求在边界上是零阶导数连续，也就是函数本身连续，所以平面和空间问题都是 $C_{0}$ 型单元，一阶单元和二阶单元都满足 $C_{0}$ 连续性。对于梁单元，板和壳单元，势能泛函中出现了位移的二阶导数，则在边界上要求一阶导数连续，梁，板，壳属于 $C_{1}$ 型单元。单元的完备性往往容易满足，但是对于 $C_{1}$ 型问题的协调性不太容易满足。满足完备性和协调性要求的单元称之为协调单元，协调单元在尺度趋向于零时，FEM分析的结果趋向于真实的结果。