形函数与等参单元

第五章形函数与等参单元1第1页内容回顾经过前面学习认识到了形函数在有限单元法中主要性，形函数作为单元内插函数，在单元位移模式中有主要作用，对单元刚度推导，外力节点等效起着关键作用；推导原理和过程明确，不过推导繁琐，只能适应简单少节点单元（常应变三角形单元等）；思变：形函数仅是单元内部位移插值函数（利用节点位移到达内部位移），可考虑利用数学中“插值函数”方法，直接给出形函数。从而避开繁琐推导假定简单代数多项式位移函数利用节点坐标和节点位移分量求得多项式的呆待定系数待定系数反代得到位移函数，推演得到形函数2第2页一．面积坐标

对于三角形单元，用面积坐标代替普通直角坐标，不但能够简化应力矩阵、刚度矩阵和载荷矩阵等运算，而且它不随三角形单元形状和方位改变，对于计算机应用也十分有利。如图所表示三角形单元ijm，任意一点P(x,y)位置，能够用以下三个比值来确定式中为三角形ijm面积，分别是三角形Pjm、Pmi、Pij面积。这三个比值称为P点面积坐标（节点对向三角形面积）。

（5-1）35.1面积坐标与自然坐标1．面积坐标定义第3页（i）三个面积坐标并不完全是独立，只有两个是独立所以（iv）

依据面积坐标定义，不难看出，在平行jm边直线上全部各点，都有相同Li坐标，而且这个坐标就等于“该直线至jm边距离”与“结点i至jm边距离”比值。图中示出了Li一些等值线。。（5-2）

（5-3）；42．面积坐标性质（ii）三角形三个顶点处（iii）三角形三条边上第4页3面积坐标和直角坐标之间关系于是（5-5a）（5-4）类似5-5b）5三角形Pjm面积是第5页

它们矩阵形式可写为

与(5-4)式对比，可见前述三角形常应变单元中形函数Ni、Nj、Nm就等于面积坐标Li、Lj、Lm。

（5-6）6

则三节点三角形单元位移函数能够写成面积坐标表示式第6页

将式（5-6）求逆，可得到同理有以及7以上三式就是面积坐标与直角坐标之间变换公式。设Li、Lj为独立变量，则Lm=1－Li－Lj，变换式能够把平面上任意三角形ijm变换为LiLj平面上三角形i1,j1,m1。第7页84面积坐标求导和积分当面积坐标函数对直角坐标求导时，能够应用以下公式可写为第8页

求面积坐标幂函数在三角形单元上积分时，能够应用积分公式

式中为整函数。求面积坐标幂函数在三角形某一边积分值时，能够应用积分公式式中l为该边长度。可得到

（5-7）（5-8）9第9页10二．四边形自然坐标（归一化处理）(-1,1)(-1,-1)O1(1,1)(1,-1)234(1,-1)(-1,-1)4(-1,1)321(1,1)OO母体单元实际单元以两对边中点建立局部坐标系，而且要求坐标值有效范围0→±1,正方形内部任意一点位置可由其局部坐标来定义,即四节点任意四边形,按照上述方法(两对边中点)定义局部坐标系，显然在这种局部坐标系下,二者是等同,二者之间有映射关系,这种基于四边形自然形状而定义坐标系成为“自然坐标系”二者之间可经过一定坐标转换公式进行转换第10页11优点该种坐标系下,不论四边形大小和形状怎样,其坐标特征是相同,所以可用统一表示式描述,能够推导统一有限元公式;以局部坐标系导出公式,有利于数值积分运算,可克服高精度单元单刚矩阵、等效节点力矩阵等因无法导出显式而必须进行积分所碰到困难；对高精度单元，单元边界能够使直线/曲线，能更加好迫近实际物体边界。第11页六节点三角形单元含有12个节点自由度，位移函数采取以下完全二次多项式：

（5-9）i(1,0,0)j(0,1,0)m(0,0,1)1

(0,1/2,1/2)2

(1/2,0,1/2)2

(1/2,1/2,0)xy1.位移函数12三六节点三角形单元第12页

待定常数a1～a6可由六个节点水平位移分量确定，a7～a12可由六个节点垂直位移分量确定。a1、a2、a3和a7、a8、a9反应了刚体位移和常应变。因为位移函数采取了二次曲线，所以在任意两个单元相邻边界上，位移分量是按抛物线改变。抛物线形状由三个系数确定，而边界上三个节点水平位移u与垂直位移v完全能够分别确定三个系数值，从而确保了边界上位移连续性。所以上述位移函数能够满足解答收敛性条件。

为了计算简便起见，把边界上节点取在三角形三条边中点，分别以1，2，3表示，如图所表示。位移函数可表示为

（5-10）13第13页式中Ni，…和N1是用面积坐标表示形函数角点：边中点：

（5-12）

（5-11）14第14页利用上述形函数可将位移函数用矩阵形式表示为式中（5-13）15第15页

利用式（5-11）与式（5-12），并将六个节点面积坐标依次代入式（5-10），将得出u分别等于ui，uj，um，u1，u2，u3；v分别等于vi，vj，vm，v1，v2，v3。同时由式（5-6）可见，面积坐标与直角坐标是线性相关，因而式（5-10）中形状函数即是面积坐标二次式，也是直角坐标二次式，都能在六个节点处给出六个节点位移，所以式（5-10）与式（5-9）是等同，一样都能确保解答收敛性。16第16页依据应变公式2．单元应变与应力17因为形函数是面积坐标函数，而面积坐标又是坐标x,y函数。所以上述求导需按照复合函数求导法则进行第17页

可导出用节点位移表示单元应变表示式式中18第18页从而得

（5-14）同理可得

（5-15）19第19页

（5-16）写成矩阵形式为20第20页式中（5-17）21第21页而22第22页

由上述各式可知，应变分量是面积坐标一次式，因而也是直角坐标一次式。故应变是按线性改变。单元应力为式中23（5-18）第23页而

因为应力矩阵[S]元素都是面积坐标一次式，也是直线坐标一次式，所以单元中应力是按线性改变，精度比常应变单元精度高。24第24页对于六节点三角形单元，节点力与节点位移之间关系依然可写成式中[k]为单元刚度矩阵，它是12×12方阵。其通式仍为

（5-19）3．单元刚度矩阵25第25页

将式（5-57）转置与式（5-58）相乘，依据式（5-47），应用以下面积坐标幂函数在单元上积分公式：(i,j,m)（5-20）(i,j,m)（5-21）(i,j,m)（5-22）26第26页

对其中各元素进行积分，再利用关系式bi

+bj

+bm=0和ci

+cj

+cm=0加以化简，最终得

（5-23）27第27页式中(i,j,m)(i,j,m)

对于平面应变问题，须在应力矩阵及刚度矩阵各个公式中，将E改换为

，将改换为。28第28页

因为位移函数是非线性，非节点荷载移置只能应用普遍公式求载荷列阵。

1.分布体积力设单元重力为W，其体积力分量为移置到各个结点体积力荷载列阵为（5-24）4．荷载移置29第29页

利用积分公式（5-8）、（5-11）、（5-12），可求得30第30页

将以上两式代入式（5-64），得载荷列阵为即Xi=Xj=Xm=X1=X2=X3=0，Yi=Yj=Ym=0，Y1=Y2=Y3=也就是说，只须向节点1，2，3（而不是i,j,m）分别移置三分之一重力。

（5-25）31第31页yxijm123lpα

2.分布面力

设单元在边界上作用有沿x方向线性改变面力，分布面力在i点密度为p，在j点为零，如图所表示。面积坐标Li在i点为1，Xi

=Xj

=Xm=X1=X2=X3=0并沿边界按线性改变。面力分量可表示为32第32页

移置到各个节点面力荷载列阵为（5-26）因为分布面力仅作用在边界上，所以；；N1=N2=Nm=0而33第33页代入式（5-26），并应用面积坐标幂函数在单元边界积分公式（5-8），对其中各元素进行积分得将积分结果代入式（5-26）得

（5-27）34第34页这就是说，只须把总面力移置到节点i，移置到节点3。依据这一标准，能够应用叠加原理求得边界上受任意线性分布面力荷载列阵。

整个弹性体平衡方程依然可取以下形式：式中[K]——总体刚度矩阵；{δ}——结点位移列阵；{R}——节点荷载列阵。采取六节点三角形单元进行计算，精度比简单三角形单元高，整理计算结果也较简单，但因为节点较多，总体刚度矩阵带宽较大，所以计算时间稍长。35第35页365.2等参数单元一等参单元基本概念在进行有限元分析时，单元离散化会带来计算误差，主要采取两种方法来降低单元离散化产生误差：1）提升单元划分密度，被称为h方法（h-method）；2）提升单元位移函数多项式阶次，被称为p方法（p-method）。。图5-1整体坐标系下四结点矩形单元同第四章方法类似，将单元位移模式选为（5-28）1整体坐标系下：

在整体坐标系中设四结点矩形单元，如图所表示，该矩形单元在x及y方向边长分别为2a和2b，厚度为t，单元编号如图（一）四节点平面矩形单元第36页37可得到（5-29）形态函数为（5-30）上述单元位移模式满足位移模式选择基本要求：反应了单元刚体位移和常应变，单元在公共边界上位移连续。在矩形单元边界上，坐标x和y其中一个取常量，所以在边界上位移是线性分布，由两个结点上位移确定。第37页38现选择以下局部坐标系，单元局部坐标系原点取在中心，并以（x,y)坐标作为位移插值函数，并将单元坐标函数归一化，采取无量纲坐标。则有局部坐标与整体坐标关系以下所以，不论原单元变长多大，其四个节点分别为2单元局部坐标系（a）注意（5-30）和（a）式关系第38页39与三结点三角形单元相比，四结点矩形单元位移模式是坐标二次函数，能够提升计算精度，但也有显著缺点，两种单元比较以下。单元类型优点缺点三结点三角形单元适应复杂形状，单元大小过渡方便计算精度低四结点矩形单元单元内应力、应变是线性改变，计算精度较高不能适应曲线边界和非正交直线边界假如任意形状四边形四结点单元仍采取矩形单元位移模式，则在公共边界上不满足位移连续性条件。为了既能得到较高计算精度，又能适应复杂边界形状，能够采取坐标变换，即经过坐标变换将局部坐标中与四边形方向无关边长为2正方形（下列图b），变换为整体坐标中四边形（下列图a）。3平面四节点等参单元第39页40图5-2a整体坐标下任意四结点四边形单元图5-2b局部坐标下四结点正方形单元如上图所表示任意四边形单元上，用等分四条边两族直线分割四边形，以两族直线中心为原点，建立局部坐标系沿和增大方向作为和轴。因为可将局部坐标系下正方形单元变换为整体坐标系四边形单元，对于局部坐标下正方形单元，可取以下位移函数a平面四节点等参单元形函数与坐标变换第40页41参考矩形单元，四结点正方形单元位移模式为其中式中是矩形单元四个节点局部坐标，其坐标值分别为其中（5-31）（5-32）（5-34）上式可简记为形态函数（5-33）母单元第41页42由式（5-31）位移表示式以及式（5-34）可知，在单元四条边上，分别有，位移线性改变，固边界上位移呈线性改变，确保了单元公共边界上位移连续性。所以给出任意四边形单元结点位移就能得到整个单元上位移，（5-31）位移模式就是所要找正确位移模式。矩形单元形函数含有与三角形单元形函数相类似特点，既含有以下关系在单元形心上第42页43对于母单元计算，他本身并没有太大使用价值，但能够利用他得到实际计算单元（整体坐标下单元），利用形函数（5-33）作以下变换（5-35）使得图5.2b中平面上几个角点分别映射为图5.2a中平面上四个角点。假如对实际计算单元位移函数依然采取母单元形函数，即能够证实他满足完备性和协调性要求，因为描述单元变形函数和描述单元几何形状函数相同，故称计算单元为等参单元。采取等参单元，使我们能够在局部坐标系中规则单元上进行单元分析，然后在映射笛卡尔坐标系实际单元上。等参单元同时含有计算精度高和适用性好特点，是有限元程序中主要采取单元形式（实际上就是等参变换）。第43页44将位移函数（5-35）代入几何方程，能够得到位移场b平面四节点等参单元几何矩阵式中第44页45依据复合函数求导法则记为J称为Jacobi矩阵可得到第45页46能够得到应力场单元刚度矩阵其子矩阵计算公式为第46页47二四边形八节点等参单元为了更加好地反应物体内应力改变，适应曲线边界，在弹性力学平面问题分析中经常使用四边形八节点等参单元。如图所表示，因为每条边上增加了一个结点，单元边是一条二次曲线，能够更加好地适应曲线边界。图四边形八结点单元图

八结点基本单元第47页48对于等参单元，先在图所表示八结点基本单元上进行分析。八结点单元一共有16个已知结点位移分量，基本单元中取以下位移模式：（5-36）该位移模式实际上是一个双二次函数，待定系数由结点位移分量确定。在单元每条边上，局部坐标位移是局部坐标二次函数，完全由边上三个结点位移值确定，所以这个位移模式满足位移连续性条件。实际单元内位移用形函数表示为，（5-37）第48页49其中形函数为将形函数归纳为（5-38）第49页50形函数在单元i结点上值为1，在其它结点上值均为0坐标变换式采取以下相同公式（5-39）代入公式（5-39），能够得到单元345边在整体坐标下参数方程：（5-40）可见在整体坐标系中，单元边是一条抛物线或退化为一条直线。第50页51如图所表示，ANSYS提供PLANE82单元是一个四边形八结点等参单元，局部坐标定义为s和t，如图所表示。PLANE82单元能够退化为三角形六结点单元。ANSYS提供Plane82单元Plane82基本单元ANSYS理论手册中给出PLANE82单元位移模式如图5-8所表示，位移模式与公式（5-18）展开后是一样。第51页52三等参单元单元分析以平面问题四边形八结点等参单元为例，介绍结构等参单元单元刚度矩阵基本过程。弹性力学平面问题单元刚度矩阵为单元应变为单元结点位移将形函数代入后，能够得到应变矩阵表示式（5-41）第52页53可得应变矩阵分块矩阵（5-42）因为等参单元形函数是局部坐标函数，所以应变矩阵[B]也是局部坐标函数。形成等参单元单元刚度矩阵需要在整体坐标系中对局部坐标函数进行积分，包含以下三个基本步骤：1）计算用局部坐标表示形函数对整体坐标x、y偏导数；2）将整体坐标系中面积积分转换为在局部坐标系中面积积分；3）用数值积分计算出单元刚度矩阵中元素。因为局部坐标与整体坐标之间存在坐标转换关系，所以形函数Ni

是局部坐标函数，同时也能够看作是整体坐标函数。由复合函数求导法则可得（一）计算形函数对整体坐标x，y偏导数第