有限元第5章-等参数单元

《有限元第5章-等参数单元》由会员分享，可在线阅读，更多相关《有限元第5章-等参数单元（74页珍藏版）》请在装配图网上搜索。

1、2021/8/1415 等参数单元等参数单元5-1 等参数单元的引入等参数单元的引入三角形单元内的应力为常量，不同单元的应力互三角形单元内的应力为常量，不同单元的应力互不相同，提高精度的方法：不相同，提高精度的方法：（1）减小单元尺寸；）减小单元尺寸；（2）提高单元插值函数的阶次。）提高单元插值函数的阶次。为了适应不规则边界，要求用曲边单元。为了适应不规则边界，要求用曲边单元。基于以上原因，引入等参数单元。基于以上原因，引入等参数单元。2021/8/1425-2 四节点四边形等参数单元四节点四边形等参数单元四节点四边形单元的位移插值函数可以写成（以四节点四边形单元的位移插值函数可以写成（以x方

2、向的位移插值函数为例）方向的位移插值函数为例）xyyxu4321 2021/8/143对于边界来讲，将带入上式，经简化可得上式中有三个待定系数，由所在单元的节点场变量值确定，但是不能由这个单元的这条边界的两个节点的场变量值唯一确定，因此相邻两单元在同一边界上的位移表达式并不一致，使相容性条件不能得到满足。这种情况该怎样处理？BAXYFEXDXu22021/8/144我们知道，矩形单元满足相容性条件。以图示为例。它有四个节点，各条边与总体坐标轴平行。单元内任意一点的位移插值函数可以包含四个待定系数 xy xyyxu43212021/8/145在矩形单元的任意一条边上，把该边的方程或带入上式，总可

3、以得到或上式只含有两个未知参数，由边界上的两个节点的位移值唯一确定。AY BX DCXuFEYu2021/8/146可见矩形单元的特点：（1） 矩形单元满足相容性条件。（2）含有一次项和常数项，故也满足收敛性条件。（3）单元插值函数含有交叉项xy，比三节点三角形单元的阶次要高。如果通过坐标变换，将任意四边形单元变换成矩形单元，只要在坐标变换中，任意四边形单元与矩形单元之间的点是一一对应的（称为坐标变换的几何相容性），而变换后的位移插值函数又满足解的收敛性条件，这两条合在一起，就能保证任意四边形在原坐标系中满足解的收敛性条件。2021/8/147即：使四节点四边形单元满足解的收敛性的途径是（1）

4、将四边形通过坐标变换，转化为矩形单元；（几何相容）（2）以四边形节点位移值作为矩形单元的节点位移值。（收敛性要求）以上两条结合，即可保证四节点四边形单元的几何相容性和有限元解的连续性。2021/8/148在建立四边形和矩形单元的坐标变换关系时应注意：四边形单元定义在总体坐标系中，而矩形单元定义在局部坐标系中。坐标系的变换是一个四边形单元到一个矩形单元的变换。矩形单元的局部坐标系，仅仅适用于每个要变换的单元。 12341, 1 1, 1 1 , 11 , 12021/8/149为此，首先讨论局部坐标系下的位移插值函数、形状函数和收敛性条件，然后再讨论具体的坐标变换。根据前述，矩形单元四个节点的位

5、移值，就是原四节点四边形单元的节点处的位移值。因此，局部坐标系下的矩形单元内任意一点的位移可以表示为或者87654321vu ,f2021/8/1410将和节点坐标 带入位移插值函数表达式，可得 443322114321vuvuvuvueeeee 1 , 1,1 , 1,1, 1,1, 1 Ae2021/8/1411其中解上面的方程从而其中 T87654321 e1 A eeNAf,14321432100000000,NNNNNNNNN Tevuvuvuvu443322112021/8/1412或其中iiiiiivNvuNu4141,1141,1141,1141,1141,4321NNNN20

6、21/8/1413写成统一的形式形函数的性质：（1） 保证位移在节点连续。又因为是双线性单元，故也保证在边界连续。（2） 保证单元包含刚体位移。这两条性质，保证了解的收敛性。 1 , 1,1 , 1,1, 1,1, 1,4 , 3 , 2 , 1,1141,44332211为：其中iiiiiiNjiNNjjiiii, 0, 1,1,41iiN2021/8/1414下面讨论坐标变换。可以证明：视整体坐标系下的四节点四边形单元的节点坐标值为“位移值”，采用与矩形单元内任意一点的插值函数完全相同的插值方式，就可以满足坐标变换的相容性（几何相容性），即iiiiiiyNyxNx4141,2021/8/1

7、415现证明如下：从四边形到矩形的坐标变换是点点对应，并能保证相邻单元的几何相容（前面的位移插值可以看成是位移相容）。所谓几何相容，即是指总体坐标系下的两四边形单元在转换到局部坐标系下的矩形单元后：（1）相邻单元的公共节点位置重合；（2）相邻单元的公共边界不开裂，不重叠，反之亦然。关于（1）因为 ，所以相邻单元的公共节点位置重合；关于（2）：局部坐标系下的矩形单元边界上的 或 保持常数，转换到总体坐标系下后，jiNNjjiiii, 0, 1,2021/8/1416边界为线性函数，该线性函数可以由边界上的两个节点坐标完全确定。因此，保证了相邻单元的公共边界既不开裂，也不重叠。对于矩形单元中的点也

8、可同样证明（提示：用通过矩形单元中任意点的水平或垂直的直线在总体坐标系和局部坐标系中的对应关系来证明）。我们看到：矩形单元的插值函数对于场变量和坐标变换完全一样，故称之为的等参数单元。如果两者不一样，就称为超参元或亚参元。在此不予介绍。2021/8/14175-3 等参数单元平面问题的有限元格式前述有限元求解的七个步骤中第13步：形成插值函数；第46步：求出单元刚度矩阵，并集成求解；第7步：用已知节点位移计算应力。对于等参元，已经得到四边形四节点的等参数单元的形状函数。下面主要讨论单元刚度矩阵的形成，即上述中的46步。2021/8/1418一 等参数单元刚度矩阵第4步：单元应变单元位移节点位移

9、的关系由平面问题几何方程和位移插值函数，有41414141,iiiiiiiiiiiixyyxvNxuNyvNyuNxxvyuyvxuyx2021/8/14194433221144332211443322114433221141414141,vxNvxNvxNvxNuyNuyNuyNuyNvyNvyNvyNvyNuxNuxNuxNuxNvNxuNyvNyuNxiiiiiiiiiiii2021/8/14204433221144332211432143214433221144332211443322114433221100000000vuvuvuvuxNyNxNyNxNyNxNyNyNyNyNyNx

10、NxNxNxNvxNvxNvxNvxNuyNuyNuyNuyNvyNvyNvyNvyNuxNuxNuxNuxN2021/8/1421 eeeeeBBBBBvuvuvuvuxNyNxNyNxNyNxNyNyNyNyNyNxNxNxNxN43214321443322114433221143214321000000002021/8/1422第5步：单元应力应变节点位移的关系由平面问题的物理方程，有第6步：节点力节点位移间的关系由虚功原理，可得节点力于节点位移间的关系式对于平面问题有 eBDyxDyx, eVTeedVBDBF)()( eeeSTeKtdxdyBDBFe)()()(2021/8/142

11、3其中积分区域为四节点四边形单元四条边所围成的区域。 需要积分。遇到的问题：（1） 不是常量，不能提到积分号外面来；（2） 是基于局部坐标而建立的，而 涉及 。（3）上述积分区域是总体坐标系下的，不是局部坐标系下的。 )()(eSTetdxdyBDBK)(eK B N ByNxNii,2021/8/1424为了能够计算上式，可以（1）从局部坐标系向总体坐标系转换；也可以（2）从总体坐标系向局部坐标系转换。统一之后，再积分。比较方便的做法：把 的积分从总体坐标系向局部坐标系转换。下面讨论坐标变换。)(eK2021/8/1425二 等参数坐标变换现在的关键是（1）求解 在局部坐标系中的表达式；（2

12、）将 转换成 。已知：根据复合函数求导法则 BdxdyddiiiiiiyNyxNx4141,yyxxyyxx2021/8/1426写成矩阵形式其中称为雅可比矩阵。 yxJyxyxyx yxyxJ2021/8/1427由此可得其中 雅可比逆矩阵； 雅可比行列式。 xxyyJJyx11 1JJyxyxyxyxJ2021/8/1428利用 的表达式，可以将形状函数 对整体坐标x,y的偏导数，转换成对局部坐标 的偏导数。例如其中yx ,iN, iiiiiiNNxxyyJNNJyNxN1141414141,iiiiiiiiiiiixNxxNxyNyyNy2021/8/1429此外，整体坐标系与局部坐标系

13、的面积微分之间的关系为到此，完成 计算中，从总体坐标向局部坐标的变换。积分区域：从总体坐标系下的四节点四边形转换为局部坐标系下的正方形（边长为2）。即上式积分域简单，但是函数形式复杂，常采用数值积分的方法。ddJdxdy )(eK 1111)(ddJtBDBKTe2021/8/1430三 能进行等参数变换的条件在 的计算中，涉及：（1） 计算中 ，否则雅可比矩阵的逆矩阵不存在；（2） 中，也要求 。因此， 是 中将总体坐标转换成局部坐标的必要条件。问题：有使 的四边形吗？如果有，什么样的四边形会导致 ？怎样避免？)(eK B0JddJdxdy 0J0J)(eK0J0J2021/8/1431为此

14、，展开雅可比矩阵 414141414141414114141414,iiiiiiiiiiiiiiiiiiiiiiiiiiiiyxyxyNxNyNxNyxyxJ2021/8/1432BaAaBaAayyxxyyxxyxyxiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii4321414141414141414144444444141414142021/8/1433其中由此可得雅可比行列式为可见， 是 的线性函数。414413412411414141,4141,4141,41iiiiiiiiiiiiiiiiiiiiyaxayaxayBxA 342132413241Ba

15、AaAaBaaaaaAaBaBaAaJJ,2021/8/1434要使 在单元上不等于0，则只需 在四个节点上具有相同的正负即可，否则， 在单元上某些地方一定等于零。以节点1为例，将局部坐标 和整体坐标 带入 中，可得其中， 为总体坐标系中节点1、2边和节点1、4边的长度， 是这两条边的夹角。JJJ1, 1,114 , 3 , 2 , 1,iyxiiJ11414121243211, 143211, 1sin1412yyxxyyxxBaAaBaAaBaAaBaAaJ14,1212021/8/1435同理，在节点2、3、4， 的值分别为J41 , 131 , 121, 1sin4341sin3432

16、sin2321JJJ xy123412342021/8/1436由于四边形的内角之和为所以要使 全部为负是不可能的。全部为正是可能的。部分为正、部分为负是可能的，但是，是不能接受的（将使 ） 大于0的条件是： 上式意味着四边形为凸四边形，而不能有一个内角大于或等于 。可见，任意四边形的任意不包括：任意两条边通过延伸在单元内出现交点。243214 , 3 , 2 , 1sinii0J4 , 3 , 2 , 1sinii4 , 3 , 2 , 10ii2021/8/1437为了保证计算精度，在划分单元时，应尽量使四边形单元的形状与正方形相差不远。 1423交点交点2021/8/1438四 单元负荷

17、向量的形成单元上的外力不一定是集中力。如果是体力和面力，就涉及到把外载荷等效到节点上去的问题。即使是集中力，也有可能不落在节点上。等效方法：虚功原理。即体力和面力对单元的外力功，等于等效的节点外力对单元的功。上式是等效的节点外力对单元的虚功等于单元内储存的虚应变能。 int*int)(* ,)(WWdVBDBWFWexteVTeTeeTexte2021/8/14391 等效体力负荷节点向量 的形成设单位体积所受到的体积力，即体积力密度（二维）则该单元上的体积力所做的虚功为其中 分别为单元节点虚位移向量和单元形状函数。)(ePF yxppP (e)(e)(e)A*A*A*tdxdyPNtdxdy

18、PNtdxdyPWeTeTTeeText eeN,*2021/8/1440 是单元内任意一点的位移而等效的节点外力的虚功为 与体积力等效的节点外力。由此 * eeN* )(*)(*ePeTeeTextFFW)(ePF 1111A)(A*)(*(e)(e)dtdPNtdxdyPNFtdxdyPNFeTeTePeTeTePeT2021/8/14412 与表面力等效的节点向量 的形成设单元表面所受到的表面力密度（二维为单位长度上）为则该单元上的表面力所做的外力虚功为)(eQF yxqqQ eeeLeTeTLTeeLTextdlQNdlQNdlQW*2021/8/1442同样，等效的节点外力的虚功为

19、与表面力等效的节点外力。从而其中 )(*)(*eQeTeeTextFFW)(eQF 11)(*)(*dlQNdlQNFdlQNFeTLeTeQLeTeTeQeTee2222dydydxdxdydxdl2021/8/1443比如，在单元 和 边上， 则比如，在单元 和 边上， 则12340d22dydxdl14230d22dydxdl 12341, 1 1, 1 1 , 11 , 12021/8/14443 没有落在单元节点上的外力集中力向单元节点的等效，即 的形成设外力没有作用在节点上，而是作用在单元边界上某点C处，力的大小为 ，该力在虚位移上所做的功为式中 分别是外力集中力作用点处的虚位移和

20、形函数在该点的值。从而)(eRF)(eCF )(*)(*)(*eCeTCeTeCTeeCeCTCextFNFNFW eCCN,* )()()(*)(*eCeTCeReCeTCeTeReTFNFFNF2021/8/1445这样，单元节点的外力由等效的体积力节点向量、等效的表面力节点向量和等效的集中力节点向量合并成为得到了单元刚度矩阵和单元节点负荷向量，就可以按节点号叠加总体刚度矩阵和总体节点负荷向量，最后得到方程组置入边界条件后，即可求解。 )()()()(eReQePeFFFF RQPFFFFK2021/8/14465-4 八节点曲边四边形等参数单元四节点四边形等参数单元的不足：（1）阶次仍嫌

21、不高；（2）直边对曲边拟合精度不高。故提出八节点曲边四边形单元。一 平面八节点等参数单元位移插值函数还是完成两件事情：（1）以总体坐标系的节点位移，作为局部坐标系下的节点位移，建立插值函数；（2）以（1）的插值函数，建立总体坐标系下的节点坐标和局部坐标系下的节点坐标的关系。2021/8/1447先讨论第一个。局部坐标系的单元和节点配置：边长为2的正方形，节点位于正方形的角点和各边中点上。节点坐标为 0 , 1,1 , 1,1 , 0,1 , 1,0 , 1,1, 1,1, 0,1, 1:8 , 2 , 1,iii347856122021/8/1448位移插值函数取为下列多项式注意到在单元边界上

22、，最高次数为二次。这样的插值函数称为双二次插值函数，相应的插值称为双二次插值。显然，它比双线性插值的阶次要高。不出现边界上的三次项是为了相容性要求。边界上的二次函数，可以由边界上的三个节点的位移值唯一确定。因为相邻单元边界上的节点是公共节点，其位移唯一，所以相邻单元边界上任意一点的位移都相等。2827265243212021/8/1449下面来讨论第（2）个问题。如果用上述插值函数联系总体坐标系下的单元坐标和局部坐标系下的单元坐标，那么，总体坐标系中相邻单元的边界和节点值也是相等的，不会开裂，也不会重叠，这样，几何协调性（相容性）也得到了保证。二 八节点等参数单元的形状函数记节点处的位移值为

23、8 , 2 , 1,ivuiii2021/8/1450则单元内任意一点的位移可以写为 为形状函数，它可以利用形状函数的性质而得到。 iiiiiiiiieeNvNuNvu818181,8 , 2 , 1,iNiijiNNNjjiiiiii, 8 , 2 , 1, 0, 1, )2(1, ) 1 (812021/8/1451其中 是节点i的局部坐标。 是双二次函数。因为形状函数的值是唯一的，所以只要求出来，就是结果。以 为例。（1） 在节点28其值为0；（2） 在节点1的值为1。8 , 2 , 1,iii8 , 2 , 1,iNi,1N 34785612,1N,1N2021/8/1452节点28在

24、直线35、57和28上，这三条直线的方程为由此，可以构造 的表达式为它是双二次函数。再利用 ，即可以求得01, 01, 01,1N111,1CN1,111N414111,1, 1111CCN11141,1N2021/8/1453同理可得再求中点。以2为例，它在节点1和节点38为0。注意到三条直线17、35和57通过这些点于是可以构造11141,11141,11141,753NNN01, 01, 0111111,22DDN2021/8/1454它也是双二次函数。再利用 ，即可以求得同理可得1,222N21211,1, 02222DDN1121,22N2826241121,1121,1121,NN

25、N2021/8/1455将上述式子写成统一形式三 等参数变换及其实现条件局部坐标到总体坐标的坐标变换式,84 11212,6 11217 , 5 , 3 , 1 11141,22iiiNiiiiiiiiiiiiiyNyxNx8181,2021/8/1456上式可以满足坐标变换的几何相容性。雅可比矩阵相应地可以写出雅可比矩阵的逆矩阵 和雅可比行列式 。 iiiiiiiiiiiiyNxNyNxNyxyxyxJ81818181, 1JJ2021/8/1457从局部坐标系到总体坐标系的变换，局部坐标系下的直边，转换成了总体系下的二次曲线（特殊情况下成为直线）。由此，可以提高单元对于曲边边界的拟合精度。

26、为了计算 ，同样要保证 。类似于4节点直边四边形单元，对于8节点曲边四边形的节点，也要作一些限制。这些限制是（1）对边延长不能相交；（2）中间点在边界的中间部分。如果位于1/3分点，计算出现较大偏差；如果位于1/4分点，则计算结果会完全不正确，甚至方程出现奇异性，计算无法进行下去。)(eK0J2021/8/1458单元划分时应做到（1）尽量接近正方形；（2）中间节点尽量位于每边的1/2分点处。2021/8/1459四 等参数单元的特点 计算精度高。有资料表明，分别采用三节点三角形线性单元和8节点等参数单元求解，节点总数相同（200），三角形单元数300个左右，等参元只用了50个。三角形在应力集

27、中处计算结果不理想，等参元与实际情况基本相同。因此，现在有限元分析大多采用等参元计算。缺点：（1）等参数变换，使程序复杂；（2）数值积分，使计算时间增加。如果认为精度还不够高，还可以在每边再加一个节点，成为12节点等参元。实际情况是8节点等参数单元求解平面问题，精度通常足够了。2021/8/1460除了上述四节点或8节点等参数单元外，还可以在单元的四条边上配置不同数目的节点，形成其它类型的等参数单元。这样的单元可以实现单元分划时的疏密过渡。1554642112332345612345672021/8/14615-5 数值积分及其应用计算单元刚度矩阵被积函数 复杂，需要使用数值积分。采用常规数值

28、积分，要提高精度，就要增加计算时间，使用中必须考虑精度和计算开销的统一。为简单计，我们从一维函数积分 开始讨论数值积分方法。之后，再将其推广到二维和三维情况。 11111111)(,ddJtBDBddfKTe,f11,f2021/8/1462一般来讲，有两种数值积分方法。一 牛顿柯特斯求积法（Newton-Cotes）基本想法：在-1,1区间取n各分点：计算各点处的函数值 。用这些函数值去构造一个m次多项式，并用这个多项式的积分代替对 的积分。实用中，将m次多项式积分后，带入相应的节点坐标和节点处的函数即可。比如，m=0时，为矩形积分公式11121nn if f 11111nkkkkfdf20

29、21/8/1463m=1时，为梯形积分公式m=2时，为辛普森积分公式按照这样的办法，可以将类似的公式一直写下去。统一的形式为 11111121nkkkkkffdf 11112111431nkkkkkkfffdf nkkkfHdf1112021/8/1464 求积系数或积分系数。因为n个函数值最高可以确定一个(n-1)次多项式，所以牛顿柯特斯积分公式具有(n-1)次的代数精度，即公式的误差是 。二 高斯求积法如果不指定积分点的位置，只给定积分点的数目n，这样，需要求解2n个未知量，即n个积分点的位置和这n个积分点位置处的函数值2n个未知量求出后，可以构造一个最高次数达(2n-1)次的多项式，其误

30、差是 。kH nO nkfkk,2 , 1,和nO22021/8/1465这样做可以提高所构造的求积公式的精度。例如，仍采用公式按照上面这样的想法得到的数值积分公式，称为高斯求积公式。高斯积分中， 的求解比较复杂。例如，积分点应是勒让德多项式 的n个不同的实根，而积分系数 也与勒让德多项式的导数有关。 nkkkfHdf111kkH和 nLkH2021/8/1466工程中是利用已被求出的 所列成的表格来查取使用的。kkH和积分点坐标k n积分系数Hk23450.577 350 269 20.774 596 669 200.861 136 311 60.339 981 043 60.906 179

31、 845 90.538 469 310 1010.555 555 5560.888 888 8890.347 854 845 10.652 145 154 90.236 926 885 10.478 628 670 50.568 888 8892021/8/1467比较牛顿柯特斯法可见，高斯法求积精度高（积分点数目相同），或者积分点数少（精度相同）。对于等参元的有限元分析，被积函数的计算十分复杂，因此，计算次数少的高斯积分法更为适用。所以，等参元积分中的数值积分都采用高斯求积法。 JtBDBfT,2021/8/1468三 高维情况的高斯求积法对于二维情况，将二重积分化为二次积分即可。其中， 为

32、高斯积分点， 为相应的积分系数。此时，平面上的积分点数有 个。 ninjjiijnjjijniiiniiniiiefHHfHHdfHdfHddfddfK111111111111111111)(,njiji,2 , 1,和njiHHji,2 , 1,和nn2021/8/1469对于三维情况也是类似的此时，空间区域有积分点 个。 ninjnkkjiijkefHHHdddfK111111111)(,nnn2021/8/1470四 合适的高斯积分阶数数值积分很费时间，高精度数值积分尤其如此。例如，对于平面问题，当n=2时，每个单元有4个积分点；当n=3时，每单元有9个积分点；当n=4时，每单元有16个

33、积分点。二维是平方，三维是立方。所以，为了提高计算速度，在保证收敛和计算精度的允许范围内，应尽量降低数值积分的阶次。经验表明：二维情况下，四节点单元 n=2，八节点单元 n=3 较好。有时，八节点单元 n=2也会有较好的计算精度。2021/8/1471小结：（1）等参元是指单元内的插值函数的公式和坐标变换的公式一致。等参元的种类很多，可以根据实际需要选用。（2）以局部坐标系中的规则的正方形构造插值函数，继而完成从局部坐标系到总体坐标系的坐标变换和单元刚度矩阵从总体坐标系到局部坐标系的转换。整个讨论和计算都是在局部坐标系中的规则的正方形单元内进行，最后在总体坐标下叠加各单元刚度矩阵，然后求解。2

34、021/8/1472（3）等参元的优点：1）单元选择的自由度大；2）能很好地模拟曲线边界；3）计算精度高；4）输入数据少。这一点对于复杂区域的求解特别突出。（4）使用等参元：1）要坐标变换；2）数值积分；3）程序编制复杂；4）形成刚度矩阵时间长。但不对计算带来本质的影响。（5）为了保证计算精度和使计算切实可行，单元分割时要满足单元不过分歪斜的要求，中间节点尽量等分边长。2021/8/1473（6）等参数单元法只能满足在单元边界上函数值的相容性条件，其一阶导数在边界上不连续。因此，原则上只能适用于二阶微分方程所描述的问题，如应力分析，稳态温度场，电磁场分析等。对于板弯曲问题所必须满足的四阶微分方程，上述等参元方法一般不适用。部分资料从网络收集整理而来，供大家参考，感谢您的关注！