1、第5章 单元形函数的讨论,在有限单元的基本理论中，形函数不仅可以用作单元的内插函数，把单元内任一点的位移用节点表示，还可以做为加权函数，将分布力等效为节点上的集中力和力矩。根据形函数的思想，首先将单元位移场函数表示为多项式形式，然后利用节点条件将多项式的待定参数表示成场函数的节点值和单元几何参数的函数，从而将场函数表示成节点值差值形式的表达式。,本章概述,第5章 单元形函数的讨论,本章重点讨论几种典型单元的形函数插值函数的构造方式，然后以三角形单元为例，讨论形函数的性质，在此基础上分析,有限元的收敛准则。,5.1 形函数构造的一般原理,单元的类型和形状决定于结构总体求解域的几何特点、问题类型和
2、求解精度。根据单元形状,可分为一维、二维、三维单元。单元插值形函数主要取决于单元的形状、节点类型和单元的节点数目。节点的类型可以是只包含场函数的节点值，也可能还包含场函数导数的节点值。是否需要场函数导数的节点值作为节点变量一般取决于单元边界上的连续性要求，如果边界上只要求函数值保持连续，称为C0型单元，若要求函数值及其一阶导数值都保持连续，则是C1型单元。,5.1 形函数构造的一般原理,在有限元中，单元插值形函数均采用不同阶次的幂函数多项式形式。对于C0型单元，单元内的未知场函数的线性变化仅用角（端）节点的参数来表示。节点参数只包含场函数的节点,值。,而对于C1型单元，节点参数中包含场函数及其
3、一阶导数的节点值。与此对应，形函数可分为Lagrange型（不需要函数在节点上的斜率或曲率）和Hermite型（需要形函数在节点上的斜率或曲率）两大类，而形函数的幂次则是指所采用的多项式的,具有一次、二次、三次、或更高次等。,幂次，可能,5.1 形函数构造的一般原理,另外，有限元形函数N是坐标x、y、z的函数，而节点 位移不是x、y、z的函数，因此静力学中的位移对坐标微分,只对节点位移向量起作用。,时，只对形函数N作用，而在动力学中位移对时间t微分,时，,5.1 .1常用单元的形函数,1. 一维一次两节点单元（杆单元）,图5-1一维一次两节点单元模型,如图5-1所示，设位移函数u(x)沿x轴呈
4、线性变化, 即,设两个节点的坐标为 ；两节点的位移分别为 ，可以代入上式并解出 ，得,写成向量形式为：,(5.1),(5.2),(5.3),5.1.1 常用单元的形函数,(5.4),(5.5),借助于Matlab软件，可以很方便的推导出上述单元形函数，具体代码如下：,位移函数u(x)记作形函数与节点参数乘积的形式,得到形函数,1. 一维一次两节点单元（杆单元）,5.1 .1常用单元的形函数,clearx1=sym(x1);x2=sym(x2);x=sym(x);,j=0:1;,v=x.j; % v=1 x;,m=1,x1;1,x2,mm=inv(m),N=v*mm,simplify(facto
5、r(N),1. 一维一次两节点单元（杆单元）,5.1 .1常用单元的形函数,2. 二维一次三节点单元（平面三角形单元）,在总体坐标系统下，任一点的某一方向的位移是设三个节点的坐标是 ， 为三个节点在某方向上的位移，具有如下关系,(5.6),(5.7),5.1 .1常用单元的形函数,2. 二维一次三节点单元（平面三角形单元）,上述推导可用如下MATLAB程序：,Clear,v=sym(1, x,y),m=sym(1,x1,y1;1,x2,y2;1,x3,y3),mm=inv(m),N=v*mm,simplify(factor(N),得到形函数矩阵如下式,(5.8),5.1 .1常用单元的形函数,
6、在总体坐标系统下，任一点的某一方向的位移是按相似的方法可以得到,(5.9),(5.10),3.三维一次四节点单元（三位四面体单元）,5.1 .1常用单元的形函数,形函数矩阵如下式,3.三维一次四节点单元（三位四面体单元）,(5.11),5.1 .1常用单元的形函数,4.一维二次三节点单元（高次单元）,图5-2一维二次三节点单元模型,(5.12),(5.13),如图5-2设位移函数为,用节点位移,代入并求解,5.1 .1常用单元的形函数,(5.14),得到,上式等号右端第一项矩阵即为形函数。,4.一维二次三节点单元（高次单元）,5.1 .1常用单元的形函数,5. 一维三次四节点单元（Lagran
7、ge型 ),(5.15),图5-3一维三次四节点单元模型,如图5-3，位移函数为三次方程为,需要四个节点参数才能唯一地确定其中的常系数。这四个节点可以分别取两个端点和两个三分点。类似地，可以得到如下形函数方程,5.1 .1常用单元的形函数,(5.16),(5.17),其中形函数中的各元素为,5.1 .1常用单元的形函数,(5.18),6. 一维三次二节点单元（Hermite型 ),图5-4一维三次二节点单元模型,如图5-4所示的一维三次两节点单元，这类单元的位移插函数为,对应的转角方程为,(5.19),6. 一维三次二节点单元（Hermite型 ),5.1 .1常用单元的形函数,用节点参数,代
8、入求解,(5.20),即,5.1 .1常用单元的形函数,得到,(5.21),其中形函数矩阵中各元素为,(5.22),5.1 .1常用单元的形函数,上述结果可用MATLAB程序进行验证：,Clear,x=sym(x);,j=0:3;,v=x.j % v=1 x x2 x3;,m=sym(1,x1,x12,x13;1,x2,x22,x23;,0,1,2*x1,3*x12;0,1,2*x2,3*x22,mm=inv(m),N=v*mm;,simplify(factor(N),),5.1 .1常用单元的形函数,(5.23),7. 二维一次四节点单元（平面四边形单元或矩形单元）,用形函数表达的位移方程如
9、下,其中形函数矩阵的元素为,i=i, j, k, l,(5.24),5.1 .1常用单元的形函数,(5.25),8. 三维一次八节点单元,在三维一次单元形函数中，函数值沿三坐标轴（x、y、z轴）呈线性变化。假设位移函数沿各坐标轴的线性变化,可写成,假设在i节点的位移值为ui，并将数值代入上式，其他各节点(j,k,l,m,n,p,q)亦类推，共有8个式子，其中第1式如下,(5.26),5.1 .1常用单元的形函数,(5.27),可是以求得系数解,则,(5.28),5.1 .1常用单元的形函数,5.1 .1常用单元的形函数,(5.28),最后得到形函数的表达式为,5.1.2 形函数的构造规律帕斯卡
10、三角形,上述位移函数的构造规律，可以根据帕斯卡三角形加以确定，同时，这样制定的位移模式，还能够满足有限元的收敛性要求。以下是几种典型情况。例如,一维两节点单元的情况，,如图5-6：,图5-5一维二节点单元的形函数组成,图5-6一维三节点单元的形函数组成,5.1.2 形函数的构造规律帕斯卡三角形,图5-8二维六节点单元的形函数组成,二维高阶单元的情况，见图5-7、8、9、10：,图5-7二维四节点单元的形函数组成,5.1.2 形函数的构造规律帕斯卡三角形,图5-10二维九节点单元的形函数组成,图5-9二维八节点单元的形函数组成,5.1.2 形函数的构造规律帕斯卡三角形,斯卡三角形上圈定相应区域,
11、（3）对应写出位移函数的插值公式，即形函数。,可以看出，形函数时可以按照帕斯卡三角形构造，具体方法 （1）按照所研究问题的维数绘制坐标轴，一维对应一个坐标轴，二维对应两个坐标轴，三维对应三个坐标轴。（2）按照所选单元的节点数，用三角形、矩形或长方体在帕,三维单元的情况，见图5-11,图5-11三维单元的形函数组成,5.2 形函数的性质,下面以平面三角形单元为例讨论形函数的一些性质（见图5-12）。平面三角形单元的形函数为,（i =1, 2 , 3）,(5.30),图5-12三角形单元,其中，,为三角形单元的面积，,为与节点坐标,有关的系数，它们分别等,于,公式中的行列式的有关代数余子式,前面已
12、经介绍了，这里不再详述,5.2 形函数的性质,对于任意一个行列式, 其任一行（或列）的元素与其相应的代数余子式的乘积之和等于行列式的值，而任一行（或列）的元素与其他行（或列）对应元素的代数余子式乘积之和为零。因此有：,第一，形函数在各单元节点上的值，具有“本点是1、它点为零”的性质，即在单元节点1上，满足,(5.31),在节点2、3上，有,(5.32),(5.33),5.2 形函数的性质,类似地有,(5.34),第二，在单元的任一节点上，三个形函数之和等于1，即,(5.35),简记为,(5.36),这说明，三个形函数中只有二个是独立的,5.2 形函数的性质,第三，三角形单元任意一条边上的形函数
13、，仅与该边的两端节点坐标有关、而与其它节点坐标,无关。,例如，在2-3 边上有,(5.36),根据形函数的这一性质可以证明，相邻单元的位移分别进行线性插值之后，在其公共边上将是连续的。例如，单元1-2-3和1-2-4具有公共边1-2。由上式可知，在1-2边上两个单元的第三个形函数都等于0，即,(5.37),5.2 形函数的性质,不论按哪个单元来计算，公共边1-2上的位移均由下式表示,(5.38),可见，在公共边上的位移u、v 将完全由公共边上的两个节点1、2的位移所确定，因而相邻单元的位移是保持连续,的,为了能够更好地理解形函数的概念，这里引入面积坐标。在如图5-13所示的三角形单元ijm中，
14、任意一点P(x , y)的位置可以用以下三个比值来确定。,5.3 用面积坐标表示的形函数,图5-13平面三角形单元的面积坐标,(5.40),式中，三角形单元ijm的面积，i 、j 、m 为三角形Pjm、Pmi、Pij的面积。Li，Lj，Lm叫做P点的面积坐标。显然，这三个面积坐标不是完全独立的，这是由于,5.3 用面积坐标表示的形函数,i +j +m =,(5.41),所以有,Li +Lj +Lm =1,(5.42),对于三角形Pjm，其面积为,(5.43),故有,(5.44),类似地有,(5.45),(5.46),5.3 用面积坐标表示的形函数,可见，前面讲述的平面三角形单元的形函数Ni 、
15、Nj 、Nm 等于面积坐标Li 、Lj 、Lm 。,容易看出，单元三个节点的面积坐标分别为,节点 i： Li =1 Lj =0 Lm =0 节点 j： Li =0 Lj =1 Lm =0 节点m： Li =0 Lj =0 Lm =1,根据面积坐标的定义，平行于jm边的某一直线上的所有各点都有相同的坐标Li，并且等于该直线至jm边的距离与节点i至jm边的距离之比，图5-13中给出了Li的一些等值线。平行于其它边的直线也有类似的情况。,5.3 用面积坐标表示的形函数,(5.47),当面积坐标的函数对直角坐标求导时，有下列公式,(5.48),不难验证，面积坐标与直角坐标之间还存在以下变换关系：,5.
16、3 用面积坐标表示的形函数,求面积坐标的幂函数在三角形单元上的积分时，有,(5.49),式中, 、 为整常数。求面积坐标的幂函数在三角形某一边上的积分值时，有,(5.50),式中, l为该边的长度。,5.4 有限元的收敛准则,对于一个数值计算方法，一般总希望随着网格的逐步细分所得到的解答能够收敛于问题的精确解。根据前面的分析，在有限元中，一旦确定了单元的形状，位移模式的选择将是非常关键的。由于载荷的移置、应力矩阵和刚度矩阵的建立都依赖于单元的位移模式，所以，如果所选择的位移模式与真实的位移分布有很大的差别，会将很难获得良好的数值解。,可以证明，对于一个给定的位移模式，其刚度系数的数值比精确值要
17、大。所以，在给定的载荷之下，有限元计算模型的变形将比实际结构的变形小。因此细分单元网格，位移近似解将由下方收敛于精确解，即得到真实解的下界。,5.4 有限元的收敛准则, 位移模式必须包含单元的刚体位移。也就是，当节点位移由某个刚体位移引起时，弹性体内将不会产生应变。所以，位移模式不但要具有描述单元本身形变的能力，而且还要具有描述由于其它单元形变而通过节点位移引起单元刚体位移的能力。例如，平面三角形单元位移模式的常数项1、4 就是用于提供刚体位移的。, 位移模式必须能包含单元的常应变。每个单元的应变一般包含两个部分：一部分是与该单元中各点的坐标位置有关的应变，另一部分是与位置坐标无关的应变（即所
18、谓的常应变）。,为了保证解答的收敛性，位移模式要满足以下三个条件，即,5.4 有限元的收敛准则, 位移模式在单元内要连续、且在相邻单元之间的位移必须协调。当选择多项式来构成位移模式时，单元内的连续性要求总是得到满足的，单元间的位移协调性，就是要求单元之间既不会出现开裂也不会出现重叠的现象。通常，当单元交界面上的位移取决于该交界面上节点的位移时，就可以保证位移的协调性。,从物理意义上看，当单元尺寸无限缩小时，每个单元中的应变,单元中的常应变的。,应趋于常量。因此，在位移模式中必须包否则就不可能使数值,元的位移模式中，与2、3、5、6 有关的线性项就是,收敛于正确解。很显然，在平面三角形含有这些常
19、应变，单,解,提供,5.4 有限元的收敛准则,在有限单元法中，把能够满足条件1和2的单元，称为完备单元；满足条件3的单元，叫做协调单元或保续单元。前面讨论过的三角形单元和矩形单元，都属于完备的协调单元。在某些梁、板及壳体分析中，要使单元满足条件3会比较困难，实践中有时也出现一些只满足条件1和2的单元，其收敛性往往也能够令人满意。放松条件3的单元，即完备而不协调的单元。不协调单元的缺点主要是不能事先确定其刚度与真实刚度之间的大小关系。但不协调单元一般不像协调单元那样刚硬（即比较柔软），因此有可能会比协调单元收敛得快。,5.4 有限元的收敛准则,在选择多项式作为单元的位移模式时，其阶次的确定要考虑
20、解答的收敛性，即单元的完备性和协调性要求。实践证明，虽然这两项确实是所要考虑的重要因素，但并不是唯一的因素。选择多项式位移模式阶次时，需要考虑的另一个因素是，所选的模式应该与局部坐标系的方位无关，这一性质称为几何各向同性。对于线性多项式，各向同性的要求通常就等价于位移模式必须包含常应变状态。对于高次位移模式，就是不应该有一个偏移的坐标方向，也就是位移形式不应该随局部坐标的更换而改变。经验证明，实现几何各向同性的一种有效方法是，可以根据巴斯卡三角形来选择二维多项式的各项。在二维多项式中，如果包含有对称轴一边的某一项，就必须同时包含有另一边的对称项。,5.4 有限元的收敛准则,选择多项式位移模式时
21、，还应考虑多项式中的项数必须等于或稍大于单元边界上的外节点的自由度数。通常取项数与单元的外节点的自由度数相等，取过多的项数是不恰当的。,5.5 等效节点载荷列阵,在结构有限元整体分析时，结构的载荷列阵R是由结构的全部单元的等效节点力集合而成，而其中单元的等效节点力Re 则是由作用在单元上的集中力、表面力和体积力分别移置到节点上，再逐点加以合成求得。本节以平面三角形单元为例，讨论集中力、表面力和体积力的等效移置方法以及如何形成结构等效载荷列阵，并与静力等效进行了对比,5.5.1 单元载荷的移置,根据虚位移原理，等效节点力所做的功与作用在单元上的集中力、表面力和体积力在任何虚位移上所做的功相等，由
22、此确定等效节点力的大小。对于平面三角形单元，有,(5.51),等号左边表示单元的等效节点力Re 所做的虚功；等号右边第一项是集中力G所做的虚功，等号右边第二项是面力q所做的虚功，积分沿着单元的边界进行；等号右边第三项表示体积力p所做的虚功，积分遍及整个单元；t为单元的厚度，假定为常量。,5.5.1 单元载荷的移置,用形函数矩阵表示的单元位移模式方程为,(5.52),代入式(5.51)，注意到节点虚位移列阵*e可以提到积分号的外面，于是有,(5.53),上式右端括号中的第一项与节点虚位移相乘等于集中力所,虚功，它是单元上的集中力移置到节点上所得到的等效节,点力,做的,5.5.1 单元载荷的移置,
23、它是一个61阶的列阵，记为Fe。同理，式(5.53)右端括号中的第二项是单元上的表面力移置到节点上所得到的等效节点力，记为Qe；第三项是单元上的体积力移置到节点上所得到的等效节点力，记为Pe。,注意到( *e ) T 的任意性，上式化简为,R e = F e +Q e +P e,(5.54),F e = N T G,其中,(5.55),(5.57),(5.56),5.5.2 结构整体载荷列阵的形成,结构载荷列阵由所有单元的等效节点载荷列阵叠加得到。注意到叠加过程中相互联接的单元之间存在大小相等方向相反的作用力和反作用力，它们之间相互抵消，因此，结构载荷列阵中只有与外载荷有关的节点有值。下面逐项
24、进行讨论。,（1）集中力的等效载荷列阵,逐点合成各单元的等效节点力，并按节点号码的顺序进行排列，组成结构的集中力等效载荷列阵F，即,(5.58),5.5.2 结构整体载荷列阵的形成,式中,(5.60),(Ni )c 、(Nj )c 、(Nm )c 为形函数在集中力作用点处的值。,（2）表面力的等效载荷列阵Q,把作用在单元边界上的表面力移置到节点上，得到各单元的表面力的等效节点力。按照节点号码的顺序进行排列，逐个节点叠加合成后，组成结构表面力的等效载荷列阵Q，即,(5.61),上式中，单元e的集中力的等效节点力为(记单元节点局部编号为i, j, m),(5.59),5.5.2 结构整体载荷列阵的
25、形成,式中，,(5.62),由于作用在单元边界上的内力在合成过程中已相互抵消，上式中的节点力只由作用在结构边界上的表面力所引起。,（3）体积力的等效载荷列阵P,与表面力类似，体积力的等效载荷列阵也是由单元体积力的,5.5.2 结构整体载荷列阵的形成,等效节点力按节点号码顺序排列，在各节点处合成得到,(5.63),式中，单元e的体积力的等效节点力为,(5.64),（4）结构整体载荷列阵R,根据式（5.58）（5.61）（5.63）求和得到R=F+Q+P,(5.65),5.5.3 载荷移置与静力等效关系,上述基于形函数的载荷等效所得到的结果与按照静力学的平行力分解原理得到的结果完全一致。,例如，如
26、图5-14所示的单元e，在ij边上作用有表面力。假设ij边的长度为l，其上任一点P距节点i的距离为s。根据面积坐标的概念，有,(5.66),代入式(5.57)，求得单元表面力的等效节点力,5.5.3 载荷移置与静力等效关系,可见，求得的结果与按照静力等效原理将表面力q向节点i及j分解所得到的分力完全相同。,图5-14 表面力等效示意图,再如，从图5-15所示的单元e的A点处取体积微元tdxdy，作用在其上的体积力为ptdxdy，为便于分析，认为力的作用方向与单元平面垂直。根据平行力分解原理，对jm边取力矩，求得节点i处的分力为,5.5.3 载荷移置与静力等效关系,(5.67),整个单元e的体积
27、力在节点i处的分力为,(5.68),类似地，分别对im及ij边取力矩，可得到节点j和节点m处的分力,(5.69),(5.70),图5-15 体积力等效示意图,5.5.3 载荷移置与静力等效关系,因此，对于平面三角形单元，按照静力学中平行力的分解原理所得到的节点力与按照虚功原理求得的节点力完全一致，在实际计算等效节点力时，可以直接应用静力学中有关平行力分解的结果。例如，对均质等厚度的三角形单元所受的重力，只要把1/3的重量直接加到每个节点上，对于作用在长度为l的ij边上强度为q的均布表面力，可以直接把 (qtl) /2 移置到节点i和j上。,5.6 根据形函数基本原理进行三维实体分析,如图5-1
28、6所示的一个长方体，长宽高分别为40cm，10cm，30cm，分成3个单元，每个单元有8个节点且几何尺寸相同,都为,。长方体顶部及底部为固定约束，在节点,。试对该,6处沿z方向作用F=1000N的力,结构进行静力学分析。,图 5-16 三维实体,5.6 根据形函数基本原理进行三维实体分析,（1）划分单元，输入坐标值,选择3维一次8节点单元，将结构分为3个单元。,（2）求单元的形函数,该单元的位移模式为,参照5.1.1节求形函数的方法，通过编制matlab程序，求得本实例的形函数，以单元1为例，形函数为,5.6 根据形函数基本原理进行三维实体分析,（3）求解单元的刚度矩阵,求得单元的形函数后，单
29、元内任意点的位移，可以用节点位移来表示，即,三维实体中有6个应变分量，其中3个正应变，3个剪应变，按照几何方程，有,5.6 根据形函数基本原理进行三维实体分析,设,为单元的应变矩阵，而单元的弹性矩阵为,5.6 根据形函数基本原理进行三维实体分析,由此可以确定单元的刚度矩阵，求解式为,三个单元的几何形状相同又是同一种单元，所以有,即3个单元的单元刚度矩阵相同。,5.6 根据形函数基本原理进行三维实体分析,（4）单元刚度矩阵的组集,组集单元刚度矩阵可以直接组集，也可以利用转换矩阵，因为本实例中共有16个节点，而每个节点有3个自由度，因此，组集成的总刚矩阵K为,维的方阵。,（5）形成求解方程式,对于
30、静力学问题，求解方程式为,（6）引入约束条件进行求解,节点1，2，3，4，13，14，15，16为固定约束，因此，对应的每个节点,方向的位移为0，按此可以对总刚矩阵进行修正，只,取总刚矩阵的13-36行和13-36列的元素，形成新的刚度矩阵,另外，因为在节点6处沿z轴正向作用有1000N的外力，所以设,其余,元素为0，对应于前面的总刚矩阵，,这里也,取总载荷列阵的13-36行的元素，形成新的载荷列阵,5.6 根据形函数基本原理进行三维实体分析,按以上分析，可以求得5-12节点的位移列阵,即,结果为,1.0e-006 *-0.1160-0.0205-0.19270.17800.00840.87330.1182-0.01990.1287,-0.16270.01240.0470-0.0816-0.02060.09660.1007-0.01790.1426,0.08590.03710.1780-0.0985-0.04700.0112,5.6 根据形函数基本原理进行三维实体分析,（7）求各节点外力、应力和应变,求出了各节点位移，紧接着可以求出各节点力。利用几何方程求出应变，利用物理方程求出应立，从而完成静力学分析。,
