1、Chap 5 有限单元法应用中的若干实际考虑,建立有限元计算模型应遵循的一般原则。采用基于最小位能原理的位移元进行有限元分析所得应力结果的性质及其近似性的表现常用的几种改善应力结果的方法。,重点和应掌握的内容,Wilson非协调元的特点和分片试验的意义及实施方法子结构方法的特点、使用条件和实施步骤有限元建模中有效利用结构对称性和周期性的方法和实施步骤。,分析过程的有效性计算结果的可靠性,5.1 单元选取与网格划分,有限元方法的两大核心,一.单元类型和形状的选择,需要考虑：,问题的维数：一、二、三维,单元类型：,实体单元,结构单元（梁、板、壳单元）,根据工程中的问题，解决的办法：,单元形状,单元
2、阶次,三角形单元比较适合不规则形状 四边形比较适合规则性状,与求解域内应力变化特点有关,例题：如图是一悬臂梁， E=10MPa ，=0.3，L/B=100，端部有一集中力P=10KN作用。计算精度与单元网格划分及网格畸变间的关系。,此悬臂梁在端部的垂直位移需要考虑横向剪切的影响，可按弹性力学解出：,T3、Q4的自由度为16，T6的自由度为42，Q8的自由度为36。,网格划分：M1至M5， 三角形单元数是矩形单元数的2倍 网格M5，T3、Q4的自由度为16，T6的自由度为42，Q8的自由度为36。,NDLT：总自由度数,网格划分与计算结果,二、网格的划分的考虑,1. 网格疏密的布置,对于解有局部
3、集中现象，加密网格,在原网格中进行重分析。即，局部分析。,采用自适应分析方法,即：对前一次分析进行误差估计。,若误差超过规定，再由程序自动加密网格， 或提高单元阶次后进行重分析，直至满足精 度要求。,三项关键技术：,误差判断 单元升阶 网格再生成,不连续处的网格自然划分,(c) 板厚度的突变,(d) 材料性质的突变,三、疏密网格的过渡,(a) 采用形状不规则的单元过渡,不同密度划分网格过渡,缺点：可能因单元形状不好而影响局部精度。,(b) 采用三角形单元过渡,不同密度划分网格过渡,不足之处：可能因引入不同形式的单元而带来不便。,(c) 采用多点约束方程过渡,不同密度划分网格过渡,需引入约束方程
4、：,区域III,例如：图示，区域I已经用5个4边形4节点单元离散，区域III已经用2个4边形7节点Serendipity单元离散。要求采用适当的单元离散方法将区域II与区域I和区域III正确地连接起来。要求说明所采用单元的类型。,方案1：采用三角形与四边形等参单元过渡,方案2：采用形状不规则的四边形等参单元过渡,方案3：采用四边形等参元附加多点约束方程过渡。其中：表示4边形5节点Serendipity过渡单元,约束方程：ub=(ua+uc)/2; ud=(ue+uc)/2; ug=(uf+uh)/2,四、无限域问题,例如：设备的基础问题，是半无限域的地基的局部承载设备的自重和工作载荷。,此问题
5、可以采用无限元来分析。,5.2 应力近似解的性质与处理,分析：,真位移实,位移近似解,故，得到：,需考虑 的极小值问题。,同样有：,不变量,应变、应力近似解的性质：,它们是真实应变和真实应力在加权最小二乘 意义上的近似解。,5. 子结构 (Substructures),子结构可以看作一个超级单元(Sup-element)，它本身可以是一些有限元组成，也可以是由下一级的子结构构成。,有限元,基本子结构,一级子结构,整体结构,同级子结构之间交界面上的结点自由度：,外部自由度 内部自由度,例如，对一个实际物体按照3维问题进行离散处理,当实际物体太大，可以考虑先离散成子结构(Substructures
6、)，用超级单元处理,一.子结构法的理论根据：,1. 最小位能原理：,不失一般性，考虑仅有一级子结构， 上述位能泛函可以写成：,转化矩阵：,表示了单元内结点编号和子结构内结点编号的关系。,表示了子结构内结点编号和整体结构内结点编号的关系。,首先可以对子结构集成方程：,2.子结构集成方程的讨论,将子结构中自由度 分离为 和,相应的方程可写为：,i) 与内部自由度有关的项,已集成完毕。,内部结点与子结构以外的任何单元不相联,所以相应的项不受其它结构的影响.,ii) 尚未集成完毕,iii) 由前面讨论过的Gauss消去法的特点，可以边集成边消元。,a) 元素集成完毕的行可以作为消元行。,b) 被消元行
7、元素可以没有集成完毕，集成与消元可以同时进行，这里我们可以先用的有关行进行消元。,我们叫做内部自由度的凝聚,二.内部自由度的凝聚,由第一组方程，可解出：,代入第二组方程得：,或：,( * ),1. 对总体方程的贡献仅外部自由度有关项,因为内部自由度已经消元完毕，没有必要 再在总刚中出现，因此大大降低了总刚的 自由度数。,由总体方程可解出 得到每一个子结构的,2. 由(*)式可以计算 ，并计算各单元应力。,（注：需要对每个子结构保留一些信息 外存）,3. 程序实现不采用 求逆的方法，而是采用消元 回代方法。,消元：将内部自由度有关各行作为消元行。,b) 、 正是集成总体方程所需的项。,c) 为回
8、代所需项，要保留 外存。,d) 总体方程求解完后，,i) 对每个子结构取出,ii) 由外存读入,iii) 由下往上回代得,注：适用范围，整体分析、分叉结构,5. 非协调元与分片试验,问题提出：经典的矩形单元，具有双线形位移模式。位移模式是具有完全一次式，非完全二次式。然而，有限元的计算精度取决于插值函数的完全项的阶数因此，导致了经典的矩形单元产生其固有的缺陷，如：剪切锁定和奇异能量模式,一. Wilson平面非协调单元,Wilson平面非协调单元的思路：引入内部位移自由度，使位移模式具有完全二次式的形式 。然而，内部位移自由度并不和实际的节点相联系，因此可以理解为增加了虚拟的内部节点这种节点不
9、能保证两个单元的界面位移连续。因此，称为“非协调单元”,在矩形双线性单元基础上，增加一位移修正项：,显然在四个结点处修正项等于零，因此它只影响单元内部位移，可见e是单元内部位移参数。所求得的结点位移将仍直接是实际问题的近似值。,Wilson平面非协调单元的构造,根据这一位移模式，可得应变矩阵为,经虚位移原理做单元列式，可得,双线性,修正项,相关项,相关项,结点力,双线性等效载荷,修正等效载荷,各项的具体计算公式见下页。,双线性,修正项等效,相关项,修正项,双线性等效,由于e是单元内部位移参数，与单元集合体（整体）的其他单元无关，因此在集成之前可做如下消去(称:静态凝聚),由此可见做了上述消去后
10、，单元刚度方程形式上和双线性单元完全一样。,协调元是肯定收敛的，非协调元收敛性不能保证。下一章将介绍如何检查非协调元的收敛性。这里只给出结论：,1) 当单元网格为矩形或平行四边形时，分析将试收敛的。 2) 任意四边形时，以= =0点的J作为单元的各点的J，这样计算的修正单元刚度也能使分析收敛。,二. 三维非协调元,1.位移模式,以三维8结点，6面体为例,其中：,其位移模式包含：,一次项已完全了，二次项缺：,选择什么样的附加插值函数为好？,其好处是在各角结点， 均为0，这样保证了：,在结点处位移是协调的。,协调性：在边界面上，位移的分布是双向二次分布，但仅有4个边界结点，而内部附加自由度是与其它
11、单元无关的。在各单元内值不变。,所以，在边界上位移不协调，此类单元称为非协调元。,2.Wilson非协调元的有限元一般格式：,代入 可得：,近似取为零,可以得到：,类似子结构法中凝聚内部自由度,可得：,内部自由度与其它单元无关，不出现在总体 方程，不增加总刚的自由度数。,二.Wilson非协调元收敛性与分片试验：,非协调,收敛性不能保证？,分片试验：,赋予单元片上各点线性位移(常应变) 检查结点i的平衡方程。,满足,单元包围 一个结点i,单元尺寸变小时，有限元解,真解,二维Wilson非协调元,因为常应变(常应力),(无体力，分布力集中力),证明：,设：,而,当J为常矩阵,平行四边形,被积函数
12、中包含的项为,分片试验区通过。,当单元形状为矩形或平行四边形时,小结：,Wilson单元为平行四边形(六面体)，J为常矩阵时，能通过分片试验，收敛性可得证。,2) 其它情况：不能通过分片试验，收敛性不能证明，单元扭曲不大的情况下，结果还是不错。,3) 对扭曲不大的单元，取J为形心值代替，则收敛性可得证。, 5.4 应力计算结果的性质和处理,一、有限元位移解的下限性质,物理意义：单元本身是连续体的一部分，应具有无穷多个自由度。当设定单元的位移函数后，自由度限制以结点位移表示的有限自由度。即：位移函数u(x,y)对单元变形进行了限制和约束。,所以，位移元得到的位移解总体不大于真解。,注：,离散后的
13、刚度 实际的刚度,二、应力精度的改进,应力计算公式,其中 是自然坐标的系数，原则上可以计算单元 内任一点的 。,实际效果：角结点处最差，其次边中点（面中心点）内部点较好。应力在界面跳跃，（不满足力边界条件）假设仅场函数连续。,三、应力抹平的概念,实际效果，角结点处最差，边中点、面中点、内部点较好。,(i)应力计算公式：,(ii)应力在单元内分布：,对p次完全的位移模式，应力分布是p-1次， 最佳应力点：相应的Gauss积分点上应力精度最好。,位移元的有限单元法,得到各结点的位移值,利用,所以应变、应力的精度比位移低一阶。,四、应力近似解的性质,应力解的近似性：,a) 单元内部一般不满足平衡方程
14、；,b) 单元与单元之间一般不连续；,c) 在力的边界条件上一般不满足力的边界条件；,2、 应力在单元内分布,对于p完备的位移模式，应力为p-1次分布。,3、 最佳应力点: 相应的Gauss积分点应力精度最好。,定义：,为p次插值， p-1次，采用P11阶Gauss积分得到2p1精度。,4、等参元的最佳应力点,位移近似解u是p次多项式，若L是m阶微分算子， 应变近似解或应力近似解是n=p-m次多项式,精确积分至少应采用n+1阶Gauss积分，达到2n+1次多项式的精度。,插值函数中完全 多项式次数 p,最佳应力点应力 精度 p-m+1,插值函数中完全 多项式次数 p,最佳应力点应力 精度 p-
15、m+1,1,1,2,2,可视单元内应力平均值，或形心处应力。所以，平均应力 =,5、单元应力抹平技术 （最简单的一种方法 ，3结点三角形）,（b）加权抹平,（a）结点应力抹平： （最简单的一种方法 ，3结点三角形）,面积加权,(i)构造一个改进的应力解，,由于位移元得到的应力场不连续，用整体应力抹平改进计算结果，得到应力场全域连续。步骤：,ne是单元结点数， 是插值函数矩阵。,是改进后的结点应力， 是改进后的应力场函数。,（c）整体抹平,因此,可以表示为：,( i=1, 2, ., N ),应力插值的单元结点可以取结点位移，与位移取相同插值函数,求得的结点位移所得的应力解,得到各个结点改进的应力值，如果 是协调， 单元交界面上是连续。,显然整体抹平工作量很大，相对于进行第二 次有限元计算。,应力在单元内抹平，减少工作量。,原全域应力场,现单元应力场,单元足够小,区别不大,（ii）单元抹平,或 只算一个方向上的应力；,单元应力抹平求解方程组只有n个，,其中：,变分：,对于等参元：Gauss积分点的应力值求得后，,二维8结点等参元，利用四个Gauss积分点 改进四个角结点的应力值。,将22Gauss积分点应力代入,可以解出： 角结点处的应力值,
