有限元基本概念和原理

有限元基本概念和原理 有限元基本概念和原理 有限元分析（FEA，Finite Element Analysis）的基本概念是用较简单的问题代替复杂问题后再求解。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的（较简单的）近似解，然后推导求解这个域总的满足条件(如结构的平衡条件），从而得到问题的解。这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。 有限元是那些集合在一起能够表示实际连续域的离散单元。有限元的概念早在几个世纪前就已产生并得到了应用，例如用多边形（有限个直线单元）逼近圆来求得圆的周长，但作为一种而被提出，则是最近的事。有限元法最初被称为矩阵近似方法，应用于航空器的结构强度计算，并由于其方便性、实用性和有效性而引起从事力学研究的科学家的浓厚兴趣。经过短短数十年的努力，随着计算机技术的快速发展和普及，有限元方法迅速从结构工程强度分析计算扩展到几乎所有的科学技术领域，成为一种丰富多彩、应用广泛并且实用高效的数值分析方法。 有限元方法与其他求解边值问题近似方法的根本区别在于它的近似性仅限于相对小的子域中。20世纪60年代初首次提出结构力学计算有限元概念的克拉夫（Clough）教授形象地将其描绘为：“有限元法=Rayleigh Ritz法＋分片函数”，即有限元法是Rayleigh Ritz法的一种局部化情况。不同于求解（往往是困难的）满足整个定义域边界条件的允许函数的Rayleigh Ritz法，有限元法将函数定义在简单几何形状（如二维问题中的三角形或任意四边形）的单元域上（分片函数），且不考虑整个定义域的复杂边界条件，这是有限元法优于其他近似方法的原因之一。 对于不同物理性质和数学模型的问题，有限元求解法的基本步骤是相同的，只是具体公式推导和运算求解不同。有限元求解问题的基本步骤通常为： 第一步：问题及求解域定义：根据实际问题近似确定求解域的物理性质和几何区域。 第二步：求解域离散化：将求解域近似为具有不同有限大小和形状且彼此相连的有限个单元组成的离散域，习惯上称为有限元网络划分。显然单元越小（网络越细）则离散域的近似程度越好，计算结果也越精确，但计算量及误差都将增大，因此求解域的离散化是有限元法的核心技术之一。 第三步：确定状态变量及控制方法：一个具体的物理问题通常可以用一组包含问题状态变量边界条件的微分方程式表示，为适合有限元求解，通常将微分方程化为等价的泛函形式。 第四步：单元推导：对单元构造一个适合的近似解，即推导有限单元的列式，其中包括选择合理的单元坐标系，建立单元试函数，以某种方法给出单元各状态变量的离散关系，从而形成单元矩阵（结构力学中称刚度阵或柔度阵）。 为保证问题求解的收敛性，单元推导有许多原则要遵循。对工程应用而言，重要的是应注意每一种单元的解题性能与约束。例如，单元形状应以规则为好，畸形时不仅精度低，而且有缺秩的危险，将导致无法求解。 第五步：总装求解：将单元总装形成离散域的总矩阵方程（联合方程组），反映对近似求解域的离散域的要求，即单元函数的连续性要满足一定的连续条件。总装是在相邻单元结点进行，状态变量及其导数（可能的话）连续性建立在结点处。 第六步：联立方程组求解和结果解释：有限元法最终导致联立方程组。联立方程组的求解可用直接法、选代法和随机法。求解结果是单元结点处状态变量的近似值。对于计算结果的质量，将通过与准则提供的允许值比较来并确定是否需要重复计算。 简言之，有限元分析可分成三个阶段，前处理、处理和后处理。前处理是建立有限元模型，完成单元网格划分；后处理则是采集处理分析结果，使用户能简便提取信息，了解计算结果。 对称边界条件 有些结构由於具有某些对称性，我们可以在对称面上施加适当的对称边界条件，这样只要建立部分的模型，既省时又省力。 这里讲的对称，不只是几何形状的对称，还包括边界条件、外力施加、材料性质的对称。如果仅有形状对称，其他条件只要有一样不对称，就不能用对称模型求解。 举一个简单的例子。图一是一个含圆孔的平板，边界条件如图所示。我们知道图一的模型上下、左右对称，因此可以简化成图二。 为什麽图二的边界条件可以设为如此？为回答这个问题，我们先将图一分为上下两部分。假设A点与B点是上下两部分相对称的点(图三)，则A点与B点的「y」方向位移方向相反，大小相等。想像A与B点越来越靠近对称面，直到几乎重合在一起(图四)，此时A、B两点仍可视为对称於对称面，两点的y方向位移仍为大小相等，方向相反。然而此时两点几乎重合，因此A、B两点的位移应该相同。为了满足上面的条件，唯一的可能就是在对称面上y方向位移等於零。图二中的另一边x方向位移等於零，可用相同方法解释。 事实上，图二的两个对称边界，还有一个条件，那就是剪应力xy为零。这个道理的解释方法与位移类似。参考图三，A、B两点由於对称，剪应力也要一模一样；当A、B两点越来越靠近对称面(图四)，两点的剪应力就必须满足牛顿第三运动定律(作用力等於反作用力，方向相反)，此时唯一的可能就是剪应力为零。因此在对称面上，剪应力为零 总之，在对称面上，垂直於对称面的位移以及作用於对称面上的剪应力皆为 在破坏力学的应用　 破坏力学是固体力学的一个分支，这门学科主要在探讨四个部分： 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的特殊元素，其长度更只要裂纹长度的百分之十至二十即可。当然以上的准则并非万用，不同的题目最好还是多试几个网格粗细，再做决定。 方向 元素有方向性。有些元素如果方向搞错，则跑不出来；有些虽然跑得出来，但结果却有问题。 以平面元素(包括平面应变、平面应力以及轴对称)为例，四个节点的编号必须是「逆时钟方向」。如果是顺时钟方向，则在有限元素的定义中，这个元素的面积小於零，就跑不出来了。 板(或壳)元素也有方向性。板元素的方向由右手定则决定，也就是元素的正向方向，由元素编号顺序的方向决定，如此也定义了板元素的上表面或下表面。当板受到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)相邻，就要特别注意方向有没有一致。如果是「手动」建立网格，更要特别注意建立出来的元素的方向，以免造成困扰。 形状函数 形状函数在有限元素法中是非常重要的一个概念，它定义了元素内部位移的分布。以一维线性元素为例(图一) 此种元素有两个节点，分别以节点1及2表示。节点1的座标为xi = -1，节点2则为xi = -2。元素中有几个节点，就有几个形状函数。因此我们有两个形状函数。形状函数有个特点：考虑第n个形状函数，若代入第n个节点的座标，其函数值为1；若代入其他节点的座标，则函数值为0。由於图一为两节点元素，所以形状函数为线性，其表示式为：其他种类的元素的形状函数，一般有限元素法书籍(或商用有限元素法软体的使用手册)皆有提到，这里不多写了。 元素的特性主要由形状函数所主宰。例如三节点之三角形元素，沿着元素内部任一方向，位移皆为线性分布。为什麽我们知道它是线性分布？只要查查该元素的形状函数立刻就明白了。由於应变为位移的对空间的一次微分，所以在三节点三角形元素内部，应变是保持不变的，因此该元素又称为等应变元素(constant strain triangle)。应变既然不变，应力也不变。同样道理，四节点金字塔3D元素，应力及应变也是保持不变的。而四边形四节点元素，沿着边长，位移为线性，但若沿着其他方向，位移则为二次曲线分布。 当我们充分掌握各种元素的特性後，利用有限元法解题，就比较能够掌握下列几个问题： 1. 该用哪种元素？ 2. 网格大小？ 3. 预期会有何种结果？ 4. 如何判读结果？ 当然要回答以上的问题，不止要对元素特性有充分了解，还要对材料力学、弹性力学等基础理论有一定的熟悉。以2D元素模拟工程梁问题为例。选用哪种元素比较好？我们知道梁的问题，基本上是受到bending作用，因此轴向应力沿着横向的分布为线性。因此我们不会选用三节点三角形元素以及四节点四边形元素，采用八节点四边形元素是比较恰当的选择。当然我们可以选择线性元素，然後网格密一点。这对简单几何形状的模型而言，虽然行得通；但若针对较复杂的模型，在前处理建网格时，会比较花时间。一般而言，要看使用者的情形而定。 针对同一个题目，使用相同大小的网格，但一个采用低阶元素，另一个则采用高阶元素。一般而言，低阶元素的模型位移较小，高阶元素的位移较大。因为高阶元素采用较高阶的形状函数，元素在变形方面，显得较为自由，容易变形。这个特性也会反映在模态分析上。一般而言，越硬的东西，自然频率较高，因此，高阶元素的模型，算出来的自然频率较高。以上两个模型，若元素数量够多，基本上结果不会有太大差异(即元素够多，计算结果收敛)。 除此之外，形状函数还会影响到分布力、体力(body force)在有限元素法中的输入、高斯积分等 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 发生, 解决方法(可能)同上 ANSYS小知识（经典） 1．绘制等值线 最近俺写文章，期刊上大都不用彩色，所以打出的云图一片模糊，无法识别，这时候可以选择出等值线图，但是等值线图也是彩色的，如何把它转成黑白的呢？开始是抓图后用Photoshop处理，太麻烦，ansys自己行不行呢？小弟琢磨了一阵，终于弄出来了，现贴出来给大家参考一下。 1 用命令jpgprf,500,100,1将背景变为白色； 2 plotctrls>device option中，把vector mode改为on,画出等值线图； 3 plotctrls>style>contour>contour labeling, 将key vector mode contour labels设为on every Nth ele,对N输入一个数值，值越大，图中的label越少； 4 plotctrls>style>colors>contour colors,将所有的系列都改为黑色; 5 如果不喜欢ansys给出的MX,MN标志，可以用plotctrls>window controls>window options把它们去掉，将MINM 后的Mix-Min Symbols改为off就可以了。 这时候，一幅清晰的黑白等值线图就出来了。 2．改彩色为灰度显示 可以自己合成颜色，我自己调配的一个颜色文件给你看看，在默认的9段表示时，可以达到灰度由白到黑，当然你还可以调整使它由黑到白（文件中的4到12颠倒即可）。plotctrls>style>color>load color map加载我上面提供的附件即可 *** COLOR MAP CREATED FOR THE WIN32 DRIVER *** 0 0 0 0 1 7 0 60 2 100 0 100 3 62 0 100 4 95 95 95 5 88 88 88 6 80 80 80 7 70 70 70 8 60 60 60 9 51 50 50 10 41 41 41 11 30 30 30 12 20 20 20 13 62 62 62 14 78 78 78 15 100 100 100 3．黑白等值线 对体和面来说，ANSYS默认的结果输出格式是云图格式，而这种彩色云图打印为黑白图像时对比很不明显，无法表达清楚，这对于发表文章来说是非常不便的。发文章所用的结果图最好是等值线图，并且最好是黑白的等值线图。笔者原来进行这项工作时一般借用photoshop等第三方软件，很麻烦，并且效果不好。现通过摸索，发现通过灵活运用ansys本身也能实现这项功能。现将步骤写给大家，感谢simwe对我的帮助。 （1）将要输出的结果调出，这时为彩色云图； （2）将云图转换为等值线图的形式 GUI：plotCtrls—>Device Options—>[/DEVI]中的vector mode 选为on 命令：/DEVICE,VECTOR,1 这时结果为彩色等值线，若直接输出，打印为黑白图像时仍然不清晰，为此需进行以下几步将图像转换为黑白形式； （3） 将背景变为白色 命令：jpgprf,500,100,1 /rep （4）对等值线中的等值线符号（图中为A,B,C等）的疏密进行调整 GUI：plotCtrls—>Style—>Contours—> Contours Labeling 在Key Vector mode contour label 中选中on every Nth elem，然后在N= 输入框中输入合适的数值，例如5，多试几次，直到疏密合适 命令：/clabel,1,5 （5）将彩色等值线变为黑色 GUI：plotCtrls—>Style—>Colors—>Contours Colors 将Items Numbered 1，Items Numbered 2等复选框中的颜色均选为黑色，图像即可变为黑白等值线图像 命令：/color,cntr,whit,1 等等 （6）最后一步：出图 GUI：plotCtrls—>Capture Image 希望对大家能有所帮助。 /DEVICE,VECTOR,1 JPGPRF,500,100,1 /CLABEL,1,5 /COLOR,CNTR,WHIT,1