关于有限元素----基础理论普及

miaorenfeng

1.对称边界条件

    有些结构由於具有某些对称性，我们可以在对称面上施加适当的对称边界条件，这样只要建立部分的模型，既省时又省力。
    这里讲的对称，不只是几何形状的对称，还包括边界条件、外力施加、材料性质的对称。如果仅有形状对称，其他条件只要有一样不对称，就不能用对称模型求解。
    举一个简单的例子。图一是一个含圆孔的平板，边界条件如图所示。我们知道图一的模型上下、左右对称，因此可以简化成图二。
    为什麽图二的边界条件可以设为如此？为回答这个问题，我们先将图一分为上下两部分。假设A点与B点是上下两部分相对称的点(图三)，则A点与B点的「y」方向位移方向相反，大小相等。想像A与B点越来越靠近对称面，直到几乎重合在一起(图四)，此时A、B两点仍可视为对称於对称面，两点的y方向位移仍为大小相等，方向相反。然而此时两点几乎重合，因此A、B两点的位移应该相同。为了满足上面的条件，唯一的可能就是在对称面上y方向位移等於零。图二中的另一边x方向位移等於零，可用相同方法解释。
    事实上，图二的两个对称边界，还有一个条件，那就是剪应力xy为零。这个道理的解释方法与位移类似。参考图三，A、B两点由於对称，剪应力也要一模一样；当A、B两点越来越靠近对称面(图四)，两点的剪应力就必须满足牛顿第三运动定律(作用力等於反作用力，方向相反)，此时唯一的可能就是剪应力为零。因此在对称面上，剪应力为零
    总之，在对称面上，垂直於对称面的位移以及作用於对称面上的剪应力皆为零
2.在破坏力学的应用

    破坏力学是固体力学的一个分支，这门学科主要在探讨四个部分：
1. 含裂纹结构的受到外力时的应力分布；
2. 含裂纹结构受到多大的外力，裂纹会成长；
3. 结构中，裂纹一旦成长，会往那个方向成长；
4. 各种工程结构抵抗裂纹成长的能力，这部分通常由实验决定。

    首先先介绍破坏力学的基本概念。图一为含裂纹结构的破坏模式，共分为第一型(mode I or opening mode)、第二型(mode II or sliding mode)以及第三型(mode III or tearing mode)。
    裂纹尖端为奇异点(应力正比於根号 r分之一(stress ~ 1/r^0.5)，r为结构中某一点与裂纹尖端的距离)，去探讨裂纹尖端(或附近)的应力有多大是没有意义的。有鉴於此，我们引入了一个参数：应力强度因子(stress intensity factor)。应力强度因子共有三个，一般写做KI、KII以及KIII，分别对应到三种不同的破坏模式。
    在线弹性破坏力学(linear elastic fracture mechanics，简称LEFM)中，外力的大小与应力强度因子成正比关系。此外，应力强度因子也和几何参数(例如裂纹长度、外力与裂纹的距离等)有关。对於同一结构，应力强度因子可视为含裂纹结构所受外力大小的一个指标；换言之，应力强度因子越大，结构越危险。当应力强度因子超过「破坏韧性(fracture toughness)，记为Kc」，裂纹开始成长。破坏韧性通常由制作标准试片由实验求得，不同的材料有不同的破坏韧性，所以破坏韧性可视为材料性质。
    至於裂缝会往那个方向成长？有好几个理论可以预测，例如能量释放率理论(energy release rate)、最大周向应力理论(maximum circumferential stress)、J积分理论(J-integral)、应变能密度理论(strain energy density theory)。以上几个理论，除了最大周向应力之外，基本上都和应变能有关，这个部分较为深入，不详述了。
    接着我们要探讨有限元素法在破坏力学的应用。前面讲过，探讨裂缝尖端附近的应力有多大是没有意义的。因此我们的重点在於，如何运用有限元素法求得裂纹尖端的应力强度因子。最常见的就是利用「四分之一节点(quarter point)」元素来模拟裂纹尖端，如图二。严格来说，四分之一节点元素不是「三角形元素」，它是由四边形元素退化而成，退化方式如图三。一般我们利用节点上的位移，求得裂纹尖端的开口位移(crack tip opening displacement)来反推得到应力强度因子。MARC及ANSYS皆用四分之一节点元素来算应力强度因子。
    以上所说的四分之一节点元素，其最原始的统御方程式和弹性力学并没有两样，只是将其中几个节点变了位置，因此形状函数也跟着改变，最後得到一个结果：这种元素内部的应力呈现根号r分之一的奇异性，因此可用来模拟裂纹。
    另外有些软体，例如NASTRAN，则利用另一种特殊元素，共有18个节点，如图四。这种元素最原始的统御方程式，和一般的弹性力学稍有不同，它一开始就假设应力正比於根号 r分之一。
    不管如何，裂纹尖端的元素与一般的元素不同，裂纹尖端网格的大小，由算出来的应力强度因子准不准来决定。如何知道准不准？可建立一个简单模型，有理论解的题目，即可比较，例如图四。以平面问题为例，若采用四分之一节点元素，元素边长大约在裂纹长度的百分之六左右；若采用NASTRAN的特殊元素，其长度更只要裂纹长度的百分之十至二十即可。当然以上的准则并非万用，不同的题目最好还是多试几个网格粗细，再做决定。
3.素的「方向」

    元素有方向性。有些元素如果方向搞错，则跑不出来；有些虽然跑得出来，但结果却有问题。

    以平面元素(包括平面应变、平面应力以及轴对称)为例，四个节点的编号必须是「逆时钟方向」。如果是顺时钟方向，则在有限元素的定义中，这个元素的面积小於零，就跑不出来了。

    板(或壳)元素也有方向性。板元素的方向由右手定则决定，也就是元素的正向方向，由元素编号顺序的方向决定，如此也定义了板元素的上表面或下表面。当板受到bending时，上表面和下表面分别会受到张应力以及压应力(当然也可能颠倒，视bending的方向而定)。图一为两个相邻的板元素受到bending作用，元素1的方向和元素2的方向不同，因此它们的上、下表面也不一致(参考图一)。在受到图一所示的bending时，元素1的上表面为张应力，下表面为压应力；元素2的上表面则为压应力，下表面为张应力。有些有限元素软体(如MARC)在计算板元素的「上表面」节点应力时，是根据两个元素的「上表面」高斯点上应力经外差到节点上再平均求得，但注意图一中两个元素的上下表面定义颠倒，所以平均後的上表面节点应力就很不准了。

    再解释清楚一点，我们先将两个元素分开，同时画上厚度，比较清楚，如图二，图一中的A点在图二则可分为A1至A4四点。如果单独看元素1，A1点应该受到张应力，A2为压应力；单独看元素2，A3为张应力，A4为压应力。当有限元素法在计算A点「上表面」应力时，是将元素1以及元素2的「上表面」应力平均(即将A1及A4的应力平均)，本来应该上表面是张应力的，被这麽一搞，应力下降，误差就来了。

    一般有限元素法软体，在做automesh时，元素的方向与surface(或area)的方向一致，因此同一个surface(或area)做出来的网格，方向都相同。但如果两个surface(或area)相邻，就要特别注意方向有没有一致。如果是「手动」建立网格，更要特别注意建立出来的元素的方向，以免造成困扰。
4.形状函数
形状函数在有限元素法中是非常重要的一个概念，它定义了元素内部位移的分布。以一维线性元素为例(图一)
此种元素有两个节点，分别以节点1及2表示。节点1的座标为xi = -1，节点2则为xi = -2。元素中有几个节点，就有几个形状函数。因此我们有两个形状函数。形状函数有个特点：考虑第n个形状函数，若代入第n个节点的座标，其函数值为1；若代入其他节点的座标，则函数值为0。由於图一为两节点元素，所以形状函数为线性，其表示式为： 其他种类的元素的形状函数，一般有限元素法书籍(或商用有限元素法软体的使用手册)皆有提到，这里不多写了。
    元素的特性主要由形状函数所主宰。例如三节点之三角形元素，沿着元素内部任一方向，位移皆为线性分布。为什麽我们知道它是线性分布？只要查查该元素的形状函数立刻就明白了。由於应变为位移的对空间的一次微分，所以在三节点三角形元素内部，应变是保持不变的，因此该元素又称为等应变元素(constant strain triangle)。应变既然不变，应力也不变。同样道理，四节点金字塔3D元素，应力及应变也是保持不变的。而四边形四节点元素，沿着边长，位移为线性，但若沿着其他方向，位移则为二次曲线分布。
    当我们充分掌握各种元素的特性後，利用有限元素法解题，就比较能够掌握下列几个问题：
1. 该用哪种元素？
2. 网格大小？
3. 预期会有何种结果？
4. 如何判读结果？
当然要回答以上的问题，不止要对元素特性有充分了解，还要对材料力学、弹性力学等基础理论有一定的熟悉。以2D元素模拟工程梁问题为例。选用哪种元素比较好？我们知道梁的问题，基本上是受到bending作用，因此轴向应力沿着横向的分布为线性。因此我们不会选用三节点三角形元素以及四节点四边形元素，采用八节点四边形元素是比较恰当的选择。当然我们可以选择线性元素，然後网格密一点。这对简单几何形状的模型而言，虽然行得通；但若针对较复杂的模型，在前处理建网格时，会比较花时间。一般而言，要看使用者的情形而定。
    针对同一个题目，使用相同大小的网格，但一个采用低阶元素，另一个则采用高阶元素。一般而言，低阶元素的模型位移较小，高阶元素的位移较大。因为高阶元素采用较高阶的形状函数，元素在变形方面，显得较为自由，容易变形。这个特性也会反映在模态分析上。一般而言，越硬的东西，自然频率较高，因此，高阶元素的模型，算出来的自然频率较高。以上两个模型，若元素数量够多，基本上结果不会有太大差异(即元素够多，计算结果收敛)。
    除此之外，形状函数还会影响到分布力、体力(body force)在有限元素法中的输入、高斯积分等
5.Shear Locking 误差

Shear Locking 这种误差, 许多有限元素的初学者可能没有听过,
但是, 若你的 FEM model 没做好, 或元素使用不当,
这种误差就会出现, 而且答案差了十万八千里.
最糟的是, 可能算出了错误答案还不知道

shear locking

shear locking 是 FEM 造成的数值误差, 发生於细长结构的分析(尤其bending),
现象: 算出 之 shear strain energy 过大, 大到不合理.
细长结构bending 之正常现象应是: shear strain energy << bending strain energy.

例: 使用low order linear element (如4-node plane element, 8-node solid element)於 bending 分析, 会发生以上现象.
这是因为 low order element 变形会 overstiffness.
相对的, 使用high order element (有 mid-side node) 较不会发生 shear locking.
但若 high order element 的 aspect ratio 过大, 仍有可能 shear locking.

解决对策:
   1. reduced integration . (用於 low order 或 high order element )
   2. incompatible mode (extra shape function, ANSYS 有使用)
       用於 low order element .

例如 ansys 的 plane 42 & solid 45 element,
内定为 incompatible mode , 以防止 shear locking.
所以, 於 细长梁厚度方向, 使用 少量 的 plane 42来分析bending cantilever beam , 也能得到好的结果.

reduced integration 用於 low order element , 会有 hourglass 困扰,
所以 incompatible mode 用於 low order element 是较好的方案,

据书上的说法, 以板元素来说,
薄板理论 (K-L plate) 因为没考虑 shear stress, 所以无 shear locking 问题.
但厚板理论 (Mindlin plate) 有考虑 shear stress进来, 所以可能会有 shear locking 发生, 解决方法(可能)同上.
注意:

本问题主要的现象为 bending.
shear locking 对於 bending 问题影响颇大 ~~!
以上各位可看到, 有两种元素误差过大, 主要起因於 shear locking.
分析应避免这两类元素,
也就是说:

  对於 3D solid elements,
   1. 使用 低阶的 hex (６面体) element, 如8节点的 , 要注意软体是否加入 incompatible mode, 也就是 extra shape function. 若没有, 就可能不准.
      对策: 使用 8-nodes hex + incompatible mode 或 20-nodes hex.

   2. 使用 tet (四面体或称三角锥) element, 注意, 4-nodes tet 的 shear locking 严重,
      因为 ANSYS 的 4-nodes tet 没有 extra shape function. (其它软体未知)
      且, 4-nodes tet 是 constant strain tet , 刚性比较大.
      避免用 4-nodes tet.
      对策: 使用 10-nodes tet.