有限元第8讲_等参单元2

1、2011/4/192011/4/19奔阻单洁及钦件应用山东建筑大学201110材料非线性问題11几何非线性问题12热传导问题13有限元Fortran程序设计14 ANSYS有限元软件期末考试2011/4/192011/4/191有限元方法概述2数理力学基础3简单杆系结构有限元法4弹性力学平面有限元方法5等参元和高斯积分6空间问题有限元法7梁结构单元8板壳问题有限元法9结构动力问题有限元法等参数单元 5.1概述 5.2等参单元定义的给出 5.3四节点四边形等參数单元 5.4平面八节点等参数单元 5.5六面体尊参单元5.6数值积分2011/4/192011/4/195.5六面体等参单元32011/

2、4/19#2011/4/19多数弹性力学问晦要按照三维空间问题来求解.三维弹 性力学问题的有限元法的基本步骤与平面问题的步骤一样在 分析三堆问麵时，所选择的单元主鼻为四体单元和六面体单 元“毎个单元节点上定义有三个位移分*5 V.三峨问师有限元法有以下两个主鼻难点：(1) 亓J无法采用人工方法完成复杂三维实体的单元划分.褊鼻有功能0大的单元划分程序，从CAD橫型揍 生成离敢的单元网格现左的有限元软件可以读入IGES. STL 等格式的SQ影交换文件六面体元的计算精度比较高. 对于備实体无法实现六面体草元的自动划分.采用四面体单元能够实现单元自动划分.但是四面体单元的计算精度比 较低。 ANSY

3、S提供的Solid45单元就是六面体八节点等 参单元，每个节点有代表心八z三个方向位 移的三个自由度(DOF, Degree of Freedom), 可以退化为五面体棱柱和四面体单元，单元局 部坐标为r、s. tf六面体八结点等参单元的 基本单元如图所示.#2011/4/19#2011/4/19L-s|v-S (2)计算规模大-SL-E N1从三维问题的单元数目大，节点自由度多，导致 计算规模大，对计算机硬件的要求很高。为缩 短计算时间，有许多问题需要采用巨型计算机, 如CRAY.或并行计算机。常用的三单元有六体八节点粉单元等参单元的位移 模式和坐标变化式采用相同的形函数，如上#2011/4

4、/19EUm*ntSol id95的基本单元二十结点基本单元Sol id45的基本单元八结点基本单元六面体八节点等釜单元的基本单元如图所示，其形函数 为.与六面体八结点等参单元相比，六面体二十结点等参单 元能更好地适应不规则的形状，计算误差比较小，基本 单元如图所示，其形函数为#2011/4/19M = |(1 + 7+“)(< s I .8)"中，为纟古点的局部坐标#2011/4/19#2011/4/19ANSYS提供的Sol id95单元是六面体二十节点等参单元， 每个节点有代表从y、z三个方向位移的三个自由度，可 以退化为五面体棱柱.五面体金字塔形和四面体单元。Sol i

5、d95单元的基本单元如图所示。ANSYSffl 供的 Solid95 单元N, =(l-/72Xl+«Xl+<) (r =104248.20)N. =1-<'X1+X1+T7)0 = 1X1445.16)其中Bp J为单元结点在局部坐标系中的坐标单元刚度矩阵为，f=jn BTDBdxdydZ#2011/4/19#2011/4/1952011/4/19¥«") (l-i) 1-7)1-H)20 V H“"ieJ d丿.-.；”k” E小砂(/. /rV川豆叨；)&(1疔(17)(1片(”"，* J】一亍)U

6、.Q>U >#2011/4/J92011/4/J9按照上节介绍的等参单元分析的基本步骤可以得到 三维单元的单元刚度矩阵。雅可比矩阵为，比3?*勢*"6576而 6-英 at-舛冬奶al-%形函数对整体坐标的偏微分可以用雅可比矩阵表 示为形函数对局部坐标的偏微分，=U雅可比矩阵计算公式只有在难可比矩阵可逆 的情况下.可以求解岀 呢 Nylyl几 吗仏心亦一站-%刎一nzznsEnzE%£!:%些蛊莎券 处一黑一韶一兀 8=2011/4/J9利用雅可比矩阵的行列式，将整体坐标系下的积 分转换为在局部坐标系下的积分。在整体坐标系中的 体积微元为，嚴后，用高斯积分计算出

7、单元刚度矩阵。同样， 用上节中类似的公式就可以在局部坐标下完成单元的 载荷移置。体力移置的公式为2011/4/J9(IV = dx x dz) = dxdydzw = £ “/刖川如“2011/4/J9微矢虽在局部坐标系中表示为,在gh的面上受到面力作用，面力移置的公式为:w -1 “口戸叩L.奶妬其中在点4 0. no. go集中力移苣的公式为:町讪羸川其中强为局部坐标系中E Tlg方向上的单位向量.5.6数值积分等参单元刚度矩阵的每个元素都是局部 坐标的函数，等参数变换后具有非常复 杂的形式，在有限元程序中不用解析的 办法来计算局部坐标系中的积分，而采 用数值积分方法。通常采用高

8、斯积分方法计算单元刚度矩 阵中的元素及等效节点载荷列阵的元素。K卜门：頤6別I丿1心52011/4/J92011/4/J9基hx间内/!个貉r“,a=h2,.，刀人格笔项式#upu为srnm"插（ft多项无2011/4/J92011/4/J92011/4/J9其中f 丿 *n-IfifTlMgriuigelA tfi 旳 C(一个需姿粉文衣力可“通过八金住遇的缶斂仅驰加权侦令余热序基本思路是：左单元上选择某些特征点（积分点），求 出峻朽屈数存这些积分点二的取伯.然壬用一些&肩菽 乘这些函数值，最后求和就可得到近似积分值。有限元 分析中.最常用的高斯数值积分法.则川扯X近似代f

9、t/UxW分式变为I=r w=£ *）武=r s【me 川 专仍*）/专”至少具有,7次 代数稱度2011/4/J92011/4/J9如果n个结点的插值型求唸 Newton-cotes求积公式具有n“次代数精度 几个常用求积公式等距分布.则前面 切防艮£0x02求积公式。梯形公式.n-1Cf(x)di * 气勺/Xa” fb)数值积分的基本思想对于一个定积分“（7（4%构造一个多项式砍臥佚11在区间内*1、点（“=门刖上败。切'0相同.即0心丿）则用败0来近似代酚（勺,枳分式变为问息是:如何构造姜项式败0使其对f（0自处好的逼近？2011/4/J9132011/4/

10、J9写成统一的形式EM为求识系数豉称为积分系数牛顿-柯斯特积分公式粕度Gauss-Legendre 求积公式 个播值结点非暮距分布结点初积分权系数可以査表nmggn毂牙点n农分收糸皴人10.00000 ooooo 000002 0<»00 OOOOO OOOOO2± 0. 57755 02691 996261 oo»oo ooooo oocoo3i 0. rfM 66692 4148J0.00000 ooooo ooooo0 SS555 55555 555560 S8®« KMCfli MWHX4±0.36113 63113 M

11、O53 ± 0. M998 10455 84S860.547X5 4S451 374540.65214 51S« 37454高斯求积法妁呉不事先叔定税分点的位而是允许这一些点位于能得!»度最好的 枳分值之处.在给Jfc枳分点败目的条件下.这样做可以提禹解构墩的求枳 公或的楮度.i-1如黒现定可以取0个分点.我幻必须求出2“个未知A<t）U = 1,2,，丿对它迓行箱确积分.并用积分结果代替原国数的积分.其误差是iiM*积公式具有代皴帕度2一1这样求出的数值积分公式嫁为离斯枳分公式.>i匚他）銘文HjgAi则右边四个n如見左边的左«1代入釈分,

12、也只含有四个M.则它们Z旬的关KAtlft完全确定Z = 2c0 + 2c2/3=Hg）七 HJ3 =叭+*； +%丿+£5+ 心 +（ /；）勺任童给定.以上两式祁要相等.因此对虫的系敷应満足以下关1L#2011/4/J9152011/4/J9高斯积分方法预先定义了,求岀被积分的函数在指定积分点上的数值，加权后求和，就得到了该函数的积分C髙斯积分方法具有最高的计算精度采用n个 积分点的高斯积分可以达到阶的精度，也 就是说，如果被积分的函数是次多项式. 用n个积分点的高斯积分可以得到精确的积分 结果。（2讪如奧际超过三次多项式，则高斯积分得岀的是梢碗似 超过三次多项式，高斯积分得到的

13、只是近似值【例】试用2点高斯积分计算I二【解】说为积極僦（-L 1）厂所以 无需积分限规格化：2点高斯积分表达式£(2-r2>/r- HiAqj+Hj/C/s)由表4.1査得：积分点:5 = 0.57735 02691 89626 I/=-0 57735 0269L 89626 j#2011/4/J9172011/4/J9权函数:枳分点函数值:垠后得 =精度比较理论解I二杆对误差头”5) = 2了=0.3368602药If1 (2r - r2)drI7=1.49563314= 2° - = 1.158772915j：aT 如 |忌y)=1497375951."

14、; . 1二维高斯积分公式! -J式中积分点和权函数仍按上表采用OPnAvid力零效结点力X hnM*J体积力等效结点力:侃严另£ HHN«!那元刚度矩阵幻=£另 施弘时0刃|丿“m-J#2011/4/J9#2011/4/J9三、高维情况的高斯求积法#2011/4/J9由上可见高斯积分能达到较高的精度。2 阶高斯积分的精度可以达到3阶插值函数 的积分精度。从几何上讲，高斯积分是 用低阶曲线通过选积分点达到较高精度。f j:/6加的訂:f几仙耐Mk l-J)l-J二另 /1 另 Hjfy®)|=英 H，m zj7对于三却况也是矣似的#2011/4/J919

15、2011/4/J9#2011/4/J9#2011/4/J9等参元中积分阶次的选择#2011/4/J9#2011/4/J9积分阶次的选择直接影响计算的精度和 计算工作量。积分阶次的选择必须保证积分的精度.很多情况下，实际选取的高斯积分点数 低于精确积分的要求，往往可以取得较 完全精确积分更好的精度。（*籀®今）用数值积分代替精确积分，无疑计算时会产生 误差。为尽量减小这一误差，希望采用尽可能 高的数值积分阶次。但是数值积分的计算时非常费时的，组集刚度 矩阵时相当多的时间花费在计算数值积分上. 高阶次的数值积分计算将使计算成本大幅度提 高。#2011/4/19对于二维4结点双线性单元，它

16、的插值函数中 包含X gn. g n项.在假设单元的IJ是常 数（单元形状为矩形或平行四边形）的情况下, 刚度矩阵的被积函数中包含1, 4,42>n2. 4 n 项。由于被积函数存E和n方向的晶髙方次为2. 所以要达到精确积分，应采用2X2阶高斯积分 （2n-1=3,被积函数为3次或3次以下的多项式. 积分将给岀精确解）。精确积分常常是由插值函数中非完全项的最高 方次要求，而决定有限元精度的是完全多项式 的方次。这些非完全的最高方次项往往不能提 高精度，反而可能带来不好的影响。取较低阶的髙斯积分，使积分精度正好保证完 全务项式方次的要求.而不包括更高次的非完 全多项式的要求.其实质是相当

17、用一种新的插 值函数替代原来的插值函数，从而一定情况下 改善了单元的精度。#2011/4/19#2011/4/19如果单元I衣常数，则需要选取更多的积分 点对于二维8结点单元也可作类似的分析。 结论是：为精确积分单元刚度矩阵，在I二 常数条件下.应采用3X3阶高斯积分。如果 1丰常数，则需要采用更高阶的高斯积分。这种高斯积分阶数等于被积函数所有项次精确 积分所需要阶数的积分方案称之为 或完全积分.正如前面已指岀的，在对单元刚度矩阵进行精 确积分的条件下，将保证当单元尺寸h不断减 小时，有限元解单调地收敛于精确解。基于最小位能原来基础上建立的位移有限元， 其解答具有下限性质。即有限元的计算模型具

18、 有较实际结构偏大的整体刚度。选取缩减积分 方案将使有限元计算模型的刚度有所降低，因 此可能有助于提高计算精度。另外，这种缩减积分方案对于泛函中包含罚函 数的情况也常常是必须的.用以保证和罚函数 相应的矩阵的奇异性（见相应教程），否则将 可能导致完全歪曲了的结果#2011/4/19212011/4/19三、完全积分当单元具有規则形状时所用的高斯积分点可以对单元刚度 矩阵中的多項式进行轿确地积分对六面体和四面体单元而宫所谓“规则积分”是播单元的边 相交成直角并且任何边上的节点是边中点线性单元要宪全积分. 则衣每一方向燮有两个积分点相应的线性平面单元變有2X2个 积分点线性空间单元要有2X2X2个

19、积分点.同理二次单元要完全积分则在每一方向要有三个积分点相 应的二次平面单元要有3X3=9个积分点二次空间单元要有 3X3X3=27个积分点.二维等參元推荐采用的积分阶数1 * *J4>'o > * ( °sd线性平面单元要有2X2个积分点 二次平面草元有3X3个积分点常用积分阶数嚴高积分阶数单元减缩积分(Reduced I ntegrat i on)只有平面四边形和空间六面体才可以用减缩积分 其他的楔形.棱锥和棱柱形实体单元必须采用完全积 分.减缩积分是扌旨在积分时，比完全积分降低一阶(每个方向少一个积分点)，如图所表示理论上讲，高斯积分的阶次选择要根据被积函

20、数的阶次决定.阶次选得越高积 分精度越高, 但所H的计算时间越长.实际计算表明.高 阶数低于被积函数所有项次精确积分 所II的积分阶数.往往可以取得较完全精确积 分更好的精度.#2011/4/J9爭仝单元数值积分计算步鼻与算例 前面分析可知，等鑫单元较规则的三角形单 元或矩形单元解题要复杂一些.除了需要进行 坐标变换相关的导数计算.矩阵计算及其求逆 计算外，还要进行数值积分坛算，荷载计算也 要应用数值积分.程序开发难度大大增加 好在计算机计算并不费事，可以充分利用对 高斯点的循环特点.但是等参单元大大减轻的 网格剖分难度，避免复杂边界带来的精度误差, 受到广泛应用.笹蟄单元数值积分计算步骤无论

21、是程序计算还是程序开发，了解筹参单元 的计算过程.对于结构分析.特别是结构非线 性分析都是很重要的1应用数值枳分计算单元刚度矩阵 步骤如下：（1）设定枳分点坐标的3、0、幺 的值：J232011/4/J9J#2011/4/J9尽管等参单元精度较高，在网格剖分时 要把握分析重点，注意具有足够疏密的单 元和适当的疏密过渡，尽量接近结构形状, 尽量采取规则单元，避免使用大钝角、边 长尺寸相差较大、严重扭曲等单元。（2）计算Jacobi矩阵的系数，行列 式|J|的值，Jacobi矩阵 的逆矩阵匹（3）计算B矩阵元素；（4）计算积分点上的03 ；J#2011/4/19#2011/4/193应用数值积分计

22、算单元应力 步骤如下：（1）取数值枳分点坐标：（2）计算Jacobi矩阵逆矩阵：（5）对髙斯点k循环nw骤1-4, 并累加I,求出01;（6）对单元循环，计算QI；2应用数值积分计算荷载列阵 步骤如下#2011/4/19#2011/4/19） itJ；伦列阵（1）设定形住;嚳协、的值包 折边界而（线冇:的高斯点坐标。（2）计算乃cobi行列式值，包括A等;（3）计算形函数并乘以Jacobi行列式值：（5）按高斯点k循环,累加汁算。（3）计算B矩阵兀素laiBiB:（4）计算单元应力矩阵圆：（5）取出与本单元结点位移：（6）做矩阵乘法amas得枳分点应 力53。#2011/4/19252011/

23、4/19【解】数值积分采用n=2来计算。 整体编号下的单元结点编号按右 手顺序为(1,2,4,3),局部编号 按(442),对应的高斯点如图 4.12。等参单元数值积分算例【例】所示平面任意四边形，只划 分成一个四结点尊参单元，其结点 坐标为 1 (0,0) , 2 (2,0) ,3(0.2) , 4 (3,3)已知，取P2ZZZ积分阶试用数值 积分计算平面应力I魂卅单元刚度 矩阵.#2011/4/19#2011/4/19积分权系数(1)由表4. 2查得积分点坐标I (-0. 5774. -0. 5774)II (0. 5774. -0. 5774)III (0. 5774, 0. 5774)

24、(-0. 5774. 0. 5774)#2011/4/19272011/4/19<2）按高斯点循环，求单元刚度矩阵。1）以第I个高斯点，枳分点坐标I（-0.5774,-0.5774）.计JfJacobi矩阵系数，将BL0I代入卜式得-0 5774赳av；期旺(r=L-4)0"-03M4 03W4 aiQ1057-0 3HI -01057 01057 0 沏1 1057 Q105 J01057 11057(bl)| J |= 1.21130.9128 -00872-0 0872 0 9127(b2)2）计算B矩阵元素讥 dx dN, dy(f-1-4)16av.09128 -0.

25、08721 (-0 39441 f -00872 09L27j-0 3944" |-0.9128 -00872(0 3944 1(0 3692 1-00872 0.9127 J - 01057(0.13O8|#2011/4/19J292011/4/19则有 ” 96S7-11.1615-27 9036-11 1615-66 9687-27 903675 940812X55(58-11.2151-4 4861-26 916331 642017.94422 99077.47682.990717 94427.4768-26 9163-4 486131 642012 656875 940S-

26、11.21513）计算高斯点刚度矩阵。 平面应力问题的弹性矩阵第I个高斯点代表的单元区域产生 的刚度矩阵S7-I117 15MHB -2SS2SS-10?OS6-100244 -M1277 -1S6S5 -OSi-O56S5 3573笑-70211 6S591 Q1S2S -1T0512 7-1575 -1CT0M -1865 5 -70211 1&4154 2S69-I 0也、1-S605 -170512 -100244 -M1277 68)912S6W26S6011060Q4W301-11277 -1002-M 012504995HOM26S602MW6891-1M55-10708

27、6-170512 1-«05 0-IW5 2859M JW1M -70211 -085-253238 曲片、-110512 O152S 68)91 -70211 35755PJ59MS 57-1117-Oj6S5 -1S6S5 41277 -100244-KF0S6-25S2：J#2011/4/1939SSC5-16657-1500®S3Q2-6»4IT£JS120S1-23017-i«r-C1MC-now-4)eop-50®571»-4500®168rrw-C7.WM-<753-1S21-B7X00&

28、;班-C1MC-row夕沪661L«C-X4P67 DOOJ【兀.Y4-12十T7巧lime10 5>、airMTT8»-a-LM09ieoi2S4T6、畑2LWD-a-rw76H5E01JiOJSS-»2i-e30001.TM7-O47M4 8805-20096cson57W*7 JDOO-IS2b6TISS676H5-20QP6YS807J对于其它3个高斯点，重复1-3步，得到 其它离斯点代表的单元区域产生的刚度 矩阵J312011/4/19J#2011/4/19-0663)-Qi53-678 卩貝25SJ7QP 1S19023SF051 -6WT-01

29、53-0663)46BQ 23SP25»5 -6SSS6-1S5J-<1153 365T-42017M1235 -ilKBT 11J012 TJSFH10YFT -251 25»5125152515SJQQ-21KB7 -&SSS6£弓乂 -67ST YSSS6 -21XB* 10-BM65T43工AJry 11X)12 -21XB?2MK VR ttIMS 4腮勺-6S»X" 12515 2S»5446MJ -120J7M-HOW 275D612SI5253S52734）计算单元刚度矩阵单元的刚度矩 阵应该是各高斯点代

30、表区域的刚度矩阵 之和，即心汰;Z-1这n2rh式（F） （1）叠加起來得J#2011/4/J9对于其它单元，只是形函数M不同. 对于三维问题，只是由二维变成三 维即由两个变量増加为（x,y,z）, 局部坐标由（gQ）变为（gtaia）, 其它计算步骤一样.332011/4/J9#2011/4/J9【例】图示单元，己知结点位移为 ul=vl=u2=v2=0 u3=0 v3=- 0209X10-5m> u4=0> v4=-0.62X 10-5m 若采用2点高斯枳分，试讣算单元积分点 应力及结点应力。从计算可知，在每一积分点上都要计算 个单元刚度矩阵元素IQ3； i二I IV）, 这些

31、 仍删称的，而且每行（列）元 素和为0,它们迭加得到的单元刚度矩阵 也是对称的，且每行（列）元素和也为0. 单元刚度矩阵还要按整体坐标扩充，贡 献到总体坐标中去。#2011/4/195）取给宦的结点位移=0.0.0.0.0.-0 62.0.-0 209TX105m【解】有限元解答的应力是按单元 给出的，3结点三角形单元和矩形单 元给出的是结点应力，而等参单元 给出的是高斯点应力。（1）单元高斯积分点应力计算.等 参单元最佳应力点为高斯枳分点应力， 取n2,则枳分点有4个。352011/4/19#2011/4/19积分点坐标-0 5774 871 = -0 5 46）计JTCTWKT应力？将式（

32、k）及给左 的结点位移式（1）代入.矩阵相乘后得#2011/4/192)讹算JdCObi矩阵逆矩阵见式（b2） I3）计算B矩阵元素 R卩式（d） 4)汁算单元应力矩阵间目国|（e）转垃.得到第I个為斯点的应力矩阵咗=5Avx4叫._ 26.99国-2 29 H（ml）#2011/4/J91057740.57740.46750.4675205774-0.57741.74410.58953057740.57742.19942.19944-0.57740.57740.58951.7441高箭点 序号仿照上述过程，将其余积分点坐 标（EE）代入计算过程，计算出 来的高斯点应力为（单位陋）：表4. 2