第4章-平面问题的有限元法

第四章平面问题的有限元法第4章平面问题的有限元法平面问题的有限单元法不仅可以用于计算分析具有平面特征的实际机械结构，还可以通过它掌握有限单元法的基本思想和基本步骤。本章详细介绍平面三角形单元的基本原理、用平面三角形单元进行结构有限元分析的主要步骤，并给出了示例。本章概述单元分析单元组集边界条件的引入4.1平面三角形单元矩阵推导平面弹性问题可以分为两类：平面应力问题：在平面应力问题中，连续体的二维坐标尺寸远大于第三维尺寸（如板），平面的法向应力可以忽略。平面应变问题：在平面应变问题中，连续体的二维坐标尺寸远小于第三维尺寸，加载平面的法向应变可以假设为零,因此分析该连续体的应力和位移时，可以通过分析其法向应变为零的一个横断面来完成。有限元方程4.1平面三角形单元的单元刚度矩阵推导平面问题可以用最简单的平面三角形常应变单元加以分析。在这节中讨论平面三角形单元的构造方法，给出用平面三角形单元求解平面应力问题的详细过程。平面三角形单元刚度矩阵的推导包括如下6个步骤：1.选择合适的单元，建立坐标系统，进行结构离散2.选择合适的位移函数3.用节点位移表示单元内部各点位移4.用节点位移表达单元内任一点的应变5.用应变和节点位移表达单元内任一点的应力6.单元刚度矩阵的形成注意矩阵的维数！！！！4.1平面三角形单元的单元刚度矩阵推导1.选择合适的单元，建立坐标系统，进行结构离散用平面三角形分析时，可以只建立一个整体坐标系OXY。将所有作用在单元上的载荷，都按虚功等效的原则移置到节点上，成为等效节点载荷。由此得到平面问题的有限元计算模型，如图4-1所示。xOy图4-1弹性体和离散化后的有限元计算模型如图4-2所示，节点编号1、2、3按逆时针顺序编排，三个节点的位置坐标分别是,和。),(11yx),(22yx),(33yx图4-2直角坐标系下平面三角形单元的节点位移和节点力对于平面问题，每个节点有x和y两个方向的自由度，即对应的位移是u和v。可以认为三角形单元共有6个自由度，即{}T，相应的单元节点力分量分别为{}T。332211,,,,,vuvuvu332211,,,,,yxyxyxFFFFFF单元数量及类型充分利用工程软件4.1平面三角形单元的单元刚度矩阵推导1.选择合适的单元，建立坐标系统，进行结构离散三角形单元的6个节点位移分量用列阵表示为三角形单元的节点载荷列阵表示为单元节点载荷列阵和节点位移列阵之间的关系可用下式表示其中，为单元刚度矩阵。对于平面三角形单元，节点位移列阵和节点载荷列阵都是6阶的，单元刚度矩阵是一个6×6阶的矩阵。在有限元中，将利用函数插值、弹性力学几何方程和物理方程、最小势能原理，建立上式具体关系。121122333{,,,,,}eTuvuvuvδδδδ123112233{,,}{,,,,,}eTTTTTxyxyxyRRRRFFFFFFeeeRkδekek(4.1)(4.2)(4.3)节点自由度和单元刚度矩阵维数的关系4.1平面三角形单元的单元刚度矩阵推导2.选择合适的位移函数考虑建立以单元节点位移表示的单元内各点位移的表达式，选择一个简单的单元位移模式，单元内各点的位移可按此位移模式由单元节点位移通过插值得到。设平面三角形单元的位移模式为由于在x和y方向的位移都是线性的，从而保证了沿接触面方向相邻单元间任意节点位移的连续性。123456uxyvxy(4.4)设单元内位移场位移模式：单元内部任一点的位移变化规律4.1平面三角形单元的单元刚度矩阵推导3.用节点位移表示单元内部各点位移三角形单元的三个节点必定满足位移模式的要求。将单元三个节点的坐标和三个节点位移都代入位移模式方程可以求解。已知单元三个节点的坐标分别为，，。对于节点1有类似的，节点2、3也按上述方法处理，三个节点的位移模型表达式可以组成方程组利用上式就可求出未知的多项式系数，即可以求得，),(11yx),(22yx),(33yx11111110000001xyxyδαAα123123{,,}[,,]eTTTTTTTTδδδδAαAAAααeδAα133221163322115333222111433221133322112333222111111121,11121,2111121,11121,21vxvxvxyvyvyvyxvyxvyxvuxuxuxyuyuyuyxuyxuyxu(4.5)(4.6)(4.7)6个方程6个未知数Matlab4.1平面三角形单元的单元刚度矩阵推导得到单元内任意一点（x，y）的位移为式中,N为形函数矩阵(Shapefunctionmatrix)，。平面三角形单元的形函数矩阵具体表达式如下其中,I为2阶单位矩阵,；为三角形单元的面积;1(,)(,)(,)eeuxyxyxyvδfαfAδNδ10000001xyxyf(x,y)1(,)xyNfA][321IIINNNN)3,2,1()(21iycxbaNiiii1001I3322111112yxyxyx(4.8)(4.11)(4.12)单元内任意点的位移同节点位移之间的关系矩阵维数？？(4.10)Matlab4.1平面三角形单元的单元刚度矩阵推导式中的其它各个系数为（1、2、3轮换）(4.13)式（4.8）经过整理可以写成如下展开形式32321323212332332211111xxxxcyyyybyxyxyxyxa33332222111121uycxbauycxbauycxbau33332222111121vycxbavycxbavycxbav(4.14)123uxy单独用水平方向位移能确定形函数吗？4.1平面三角形单元的单元刚度矩阵推导三角形单元用于解决弹性力学平面问题，单元内任一点的应变列阵满足几何方程式中，和是线应变，是剪应变。上式中的分别用位移模式方程式(4.8)代入可求解应变分量。由于N是x,y的函数，对其进行偏微分处理后，可得4.用节点位移表达单元内任一点的应变xvyuyvxuxyyxεeebcbcbccccbbbxNyNxNyNyNxNyNxNyNyNxNxNδδε33221132132133232321112100000021000000(4.15)xyxy(4.16)1122331122331212uNuNuNuvNvNvNv4.1平面三角形单元的单元刚度矩阵推导其中,B为单元应变矩阵，其表达式为由于和b1、b2、b3、c1、c2、c3等都是常量，所以平面三角形单元的应变矩阵B中的诸元素都是常量，因而平面三角形单元中各点的应变分量也都是常量。33221132132100000021bcbcbccccbbbBeBδε(4.17)上式可简记为(4.18)说明求解精度不高4.1平面三角形单元的单元刚度矩阵推导对于平面应力问题，一点的应力状态可以用、、这三个应力分量来表示，应力应变关系为式中，D为弹性矩阵，其表达式为其中，E为杨氏模量，为泊松比。把步骤(4)中导出的应变表达式(4.16)代入式(4.19)，可得用节点位移表示的应力。令S为应力矩阵，它是),(yxσxyxy),(),(yxyxDεσ2100010112EDeDBδσDBS5.用应变和节点位移表达单元内任一点的应力(4.20)(4.21)(4.22)(4.19)常应力4.1平面三角形单元的单元刚度矩阵推导用虚位移原理对图4-2中的单元建立节点力和节点位移之间的关系。该单元在等效节点力的作用下处于平衡。设单元节点载荷列阵为，在单元中有虚位移，相应的三个节点虚位移为。作用在单元体上的外力所做的虚功为单元内任一点的虚位移也具有与真实位移相同的位移模式，即因此，由式(4.17)，单元内的虚应变*为于是，单元的应变能为这里假定单元厚度t为常量。引入单元应变的具体表达式(4.17)，注意到虚位移的任意性，可将提到积分号的前面，有6.单元刚度矩阵的形成eReδeeTURδe*Nδδe*BδεddTWtxyεσddeTTeWtxyδBDBδ(4.23)(4.24)(4.25)(4.26)eTδ(4.27)4.1平面三角形单元的单元刚度矩阵推导根据虚位移原理，，即去掉等号两边的，得此为表征单元的节点力和节点位移之间关系的刚度方程，记作式中就是单元刚度矩阵。对于材料是均质的单元，D的元素就是常量，并且对于平面三角形单元，B矩阵中的元素也是常量。当单元的厚度也是常量时，，单元刚度矩阵可以简化为WUddeTeeTTetxyδRδBDBδeT*δddeTetxyRBDBδeeeδkRddeTtxykBDBekdxdytTeDBBk(4.28)(4.29)(4.30)(4.31)(4.32)4.1平面三角形单元的单元刚度矩阵推导单元刚度矩阵的物理意义是，其任一列的元素分别等于该单元的某个节点沿坐标方向发生单位位移时，在各节点上所引起的节点力。单元的刚度取决于单元的大小、方向和弹性常数，而与单元的位置无关，即不随单元或坐标轴的平行移动而改变。单元刚度矩阵一般具有如下三个特性：对称性、奇异性和具有分块形式。对于平面三角形单元，按照每个节点两个自由度的构成方式，可以将单元刚度矩阵列写成3×3个子块、每个子块为2×2阶的分块矩阵的形式。什么含义？？111213212223313233eeeeeeeeeekkkkkkkkkk(4.32)4.1平面三角形单元的单元刚度矩阵推导[附注]形函数矩阵N和应变矩阵B的程序推导Clearsymsxyx1y1x2y2x3y3;F=[1xy];A=[1x1y1;1x2y2;1x3y3];N=f*inv(A);simplify(factor(N));N1=N(1,1);N2=N(1,2);N3=N(1,3);b1=diff(N1,x);b2=diff(N2,x);b3=diff(N3,x);c1=diff(N1,y);c2=diff(N2,y);c3=diff(N3,y);AA=det(A)B=[b10b20b300c10c20c3c1b1c2b2c3b3]/(2*AA)4.2利用平面三角形单元进行结构整体分析整体分析主要包括：单元的组集（整体有限元方程的建立）边界条件的引入及求解4.2.1单元的组集对于平面问题，每个节点有x和y两个方向的自由度。首先，引入整个弹性体的节点位移列阵2n×1，它由所有节点位移按节点整体编号顺序从小到大排列而成，即其中节点i的位移分量为再者，确定结构整体载荷列阵。设某单元三个节点（对应的整体编号分别为i,j,m，每个单元三个节点的等效节点力分别记为，其中。将弹性体的所有单元的节点力列阵加以扩充，使之成为2n×1阶的列阵，即1.直接组集法形成有限元计算模型2112{}TTTTnnδδδδT{}(1,2,,n)iiiuviδemejeiRRR,,TeyiexieiFF}{RTTTT121emejeiennmjiRRRR(4.34)(4.35)(4.36)e61R待求的量4.2.1单元的组集由于结构整体载荷列阵是由移置到节点上的等效节点载荷按节点号码对应叠加而成，相邻单元公共边内力引起的等效节点力在叠加过程中必然会全部相互抵消，所以结构整体载荷列阵只会剩下外载荷所引起的等效节点力，因此在结构整体载荷列阵中大量元素一般都为0值。21211{}NeTnne