第四章平面问题的有限元法ppt课件.ppt

第四章平面问题的有限元法,第4章 平面问题的有限元法,平面问题的有限单元法不仅可以用于计算分析具有平面特征的实际机械结构，还可以通过它掌握有限单元法的基本思想和基本步骤。本章详细介绍平面三角形单元的基本原理、用平面三角形单元进行结构有限元分析的主要步骤，并给出了示例。,本章概述,4.1 平面三角形单元矩阵推导,平面弹性问题可以分为两类：,平面应力问题：在平面应力问题中，连续体的二维坐标尺寸远大于第三维尺寸（如板），平面的法向应力可以忽略。平面应变问题：在平面应变问题中，连续体的二维坐标尺寸远小于第三维尺寸，加载平面的法向应变可以假设为零,因此分析该连续体的应力和位移时，可以通过分析其法向应变为零的一个横断面来完成。,4.1 平面三角形单元的单元刚度矩阵推导,平面问题可以用最简单的平面三角形常应变单元加以分析。在这节中讨论平面三角形单元的构造方法，给出用平面三角形单元求解平面应力问题的详细过程。平面三角形单元刚度矩阵的推导包括如下6个步骤：,1.选择合适的单元，建立坐标系统，进行结构离散 2.选择合适的位移函数3.用节点位移表示单元内部各点位移4.用节点位移表达单元内任一点的应变5.用应变和节点位移表达单元内任一点的应力6.单元刚度矩阵的形成,注意矩阵的维数！,4.1 平面三角形单元的单元刚度矩阵推导,1.选择合适的单元，建立坐标系统，进行结构离散,用平面三角形分析时，可以只建立一个整体坐标系OXY。将所有作用在单元上的载荷，都按虚功等效的原则移置到节点上，成为等效节点载荷。由此得到平面问题的有限元计算模型，如图4-1所示。,图4-1 弹性体和离散化后的有限元计算模型,如图4-2所示，节点编号1、2、3 按逆时针顺序编排，三个节点的位置坐标分别是,和。,图4-2 直角坐标系下平面三角形单元的节点位移和节点力,对于平面问题，每个节点有x和y两个方向的自由度，即对应的位移是u和v。可以认为三角形单元共有6个自由度，即 T，相应的单元节点力分量分别为 T。,4.1 平面三角形单元的单元刚度矩阵推导,1.选择合适的单元，建立坐标系统，进行结构离散,三角形单元的6个节点位移分量用列阵表示为 三角形单元的节点载荷列阵表示为 单元节点载荷列阵和节点位移列阵之间的关系可用下式表示 其中，为单元刚度矩阵。对于平面三角形单元，节点位移列阵和节点载荷列阵都是6阶的，单元刚度矩阵 是一个66阶的矩阵。在有限元中，将利用函数插值、弹性力学几何方程和物理方程、最小势能原理，建立上式具体关系。,(4.1),(4.2),(4.3),4.1 平面三角形单元的单元刚度矩阵推导,2.选择合适的位移函数,考虑建立以单元节点位移表示的单元内各点位移的表达式，选择一个简单的单元位移模式，单元内各点的位移可按此位移模式由单元节点位移通过插值得到。设平面三角形单元的位移模式为 由于在x和y方向的位移都是线性的，从而保证了沿接触面方向相邻单元间任意节点位移的连续性。,(4.4),4.1 平面三角形单元的单元刚度矩阵推导,3.用节点位移表示单元内部各点位移,三角形单元的三个节点必定满足位移模式的要求。将单元三个节点的坐标和三个节点位移都代入位移模式方程可以求解。已知单元三个节点的坐标分别为，。对于节点1有 类似的，节点2、3也按上述方法处理，三个节点的位移模型表达式可以组成方程组 利用上式就可求出未知的多项式系数，即 可以求得，,(4.5),(4.6),(4.7),6个方程6个未知数,4.1 平面三角形单元的单元刚度矩阵推导,得到单元内任意一点（x，y）的位移为 式中,N为形函数矩阵(Shape function matrix)，。平面三角形单元的形函数矩阵具体表达式如下 其中,I为2阶单位矩阵,；为三角形单元的面积;,(4.8),(4.11),(4.12),单元内任意点的位移同节点位移之间的关系,矩阵维数？,(4.10),4.1 平面三角形单元的单元刚度矩阵推导,式中的其它各个系数为（1、2、3轮换）(4.13)式（4.8）经过整理可以写成如下展开形式,(4.14),4.1 平面三角形单元的单元刚度矩阵推导,三角形单元用于解决弹性力学平面问题，单元内任一点的应变列阵满足几何方程 式中，和 是线应变，是剪应变。上式中的分别用位移模式方程式(4.8)代入可求解应变分量。由于N是x,y的函数，对其进行偏微分处理后，可得,4.用节点位移表达单元内任一点的应变,(4.15),(4.16),4.1 平面三角形单元的单元刚度矩阵推导,其中,B为单元应变矩阵，其表达式为 由于和b1、b2、b3、c1、c2、c3 等都是常量，所以平面三角形单元的应变矩阵B中的诸元素都是常量，因而平面三角形单元中各点的应变分量也都是常量。,(4.17),上式可简记为,(4.18),说明求解精度不高,4.1 平面三角形单元的单元刚度矩阵推导,对于平面应力问题，一点的应力状态 可以用、这三个应力分量来表示，应力应变关系为 式中，D为弹性矩阵，其表达式为 其中，E为杨氏模量，为泊松比。把步骤(4)中导出的应变表达式(4.16)代入式(4.19)，可得用节点位移表示的应力。令S为应力矩阵，它是,5.用应变和节点位移表达单元内任一点的应力,(4.20),(4.21),(4.22),(4.19),4.1 平面三角形单元的单元刚度矩阵推导,用虚位移原理对图4-2中的单元建立节点力和节点位移之间的关系。该单元在等效节点力的作用下处于平衡。设单元节点载荷列阵为，在单元中有虚位移，相应的三个节点虚位移为。作用在单元体上的外力所做的虚功为 单元内任一点的虚位移也具有与真实位移相同的位移模式，即 因此，由式(4.17)，单元内的虚应变*为 于是，单元的应变能为这里假定单元厚度t为常量。引入单元应变的具体表达式(4.17)，注意到虚位移的任意性，可将 提到积分号的前面，有,6.单元刚度矩阵的形成,(4.23),(4.24),(4.25),(4.26),(4.27),4.1 平面三角形单元的单元刚度矩阵推导,根据虚位移原理，即去掉等号两边的，得 此为表征单元的节点力和节点位移之间关系的刚度方程，记作式中 就是单元刚度矩阵。对于材料是均质的单元，D的元素就是常量，并且对于平面三角形单元，B矩阵中的元素也是常量。当单元的厚度也是常量时，单元刚度矩阵可以简化为,(4.28),(4.29),(4.30),(4.31),(4.32),4.1 平面三角形单元的单元刚度矩阵推导,单元刚度矩阵的物理意义是，其任一列的元素分别等于该单元的某个节点沿坐标方向发生单位位移时，在各节点上所引起的节点力。单元的刚度取决于单元的大小、方向和弹性常数，而与单元的位置无关，即不随单元或坐标轴的平行移动而改变。单元刚度矩阵一般具有如下三个特性：对称性、奇异性和具有分块形式。对于平面三角形单元，按照每个节点两个自由度的构成方式，可以将单元刚度矩阵列写成33个子块、每个子块为22阶的分块矩阵的形式。,什么含义？,(4.32),4.1 平面三角形单元的单元刚度矩阵推导,附注 形函数矩阵N和应变矩阵B的程序推导Clearsyms x y x1 y1 x2 y2 x3 y3;F=1 x y;A=1 x1 y1;1 x2 y2;1 x3 y3;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=b1 0 b2 0 b3 0 0 c1 0 c2 0 c3 c1 b1 c2 b2 c3 b3/(2*AA),4.2 利用平面三角形单元进行结构整体分析,整体分析主要包括：单元的组集（整体有限元方程的建立）边界条件的引入及求解,4.2.1 单元的组集,对于平面问题，每个节点有x和y两个方向的自由度。首先，引入整个弹性体的节点位移列阵2n1，它由所有节点位移按节点整体编号顺序从小到大排列而成，即其中节点i的位移分量为 再者，确定结构整体载荷列阵。设某单元三个节点（对应的整体编号分别为i,j,m，每个单元三个节点的等效节点力分别记为，其中。将弹性体的所有单元的节点力列阵加以扩充，使之成为2n1阶的列阵，即,1.直接组集法形成有限元计算模型,(4.34),(4.35),(4.36),4.2.1 单元的组集,由于结构整体载荷列阵是由移置到节点上的等效节点载荷按节点号码对应叠加而成，相邻单元公共边内力引起的等效节点力在叠加过程中必然会全部相互抵消，所以结构整体载荷列阵只会剩下外载荷所引起的等效节点力，因此在结构整体载荷列阵中大量元素一般都为0值。,各单元的节点力列阵经过扩充之后就可以进行相加。把全部单元的节点力列阵叠加在一起，便可得到整个弹性体的载荷列阵R。结构整体载荷向量记为,其中节点i上的等效节点载荷是,(i=1,2,n),(4.37),(4.38),4.2.1 单元的组集,再者，直接集成结构的整体刚度矩阵。把平面三角形单元的6阶单元刚度矩阵 进行扩充，使之成为一个2n2n阶的方阵。单元三个节点（1，2，3节点）分别对应的整体编号i,j,m,即单元刚度矩阵 中的22阶子矩阵kij将处于扩展矩阵中的第i双行、第j双列中。扩充后的单元刚度矩阵 为,(4.39),整体刚度矩阵的维数！,节点自由度，节点数,单元刚度矩阵经过扩充以后，除了对应的i,j,m 双行和双列上的九个子矩阵之外，其余元素均为零。把上式对N个单元进行求和叠加，得到结构整体刚度矩阵，记为 结构整体的有限元方程也可以根据虚功原理建立起来。用整体刚度矩阵、节点位移向量和节点载荷向量表达的结构有限元方程为 K=R 这是一个关于节点位移的2n阶线性方程组。整体刚度矩阵K中每一列元素的物理意义为：欲使弹性体的某一节点在坐标轴方向发生单位位移、而其它节点都保持为零的变形状态，在各节点上所需要施加的节点力。具有如下性质:K中主对角元素总是正的，它是一个对称矩阵，它是一个带状稀疏矩阵。整体刚度矩阵K是一个奇异矩阵，在排除刚体位移之后，它是一个正定矩阵。,4.2.1 单元的组集,(4.40),(4.41),4.2.1 单元的组集,对于离散化的弹性体有限元计算模型，首先求得或列出的是各个单元的刚度矩阵、单元位移列阵和单元载荷列阵。在进行整体分析时，需把结构的各项矩阵表达成各个单元对应矩阵之和，同时要求单元各项矩阵的阶数和结构各项矩阵的阶数相同。为此，引入单元节点自由度对应扩充为结构节点自由度的转换矩阵G。设结构的节点总数为n，某平面三角形单元对应的整体节点序号为i,j,m，该单元节点自由度的转换矩阵为,2.形成整体刚度矩阵的转换矩阵法,(4.42),转换矩阵G的维数？,单元自由度，总刚的维数,4.2.1 单元的组集,也就是，在G矩阵中，单元三个节点对应的整体编号位置（i,j,m）所在的子块设为2阶单位矩阵，其它均为0。利用转换矩阵可以直接求和得到结构的整体刚度矩阵为 结构节点载荷向量为,(4.43),(4.44),4.2.1 单元的组集,规律？,4.2.2 边界条件的引入,1.整体刚度矩阵的性质,弹性体有限元的整体刚度矩阵具有如下性质:第一、整体刚度矩阵K中每一列元素的物理意义为：欲使弹性体的某一节点在坐标轴方向发生单位位移、而其它节点都保持为零的变形状态，在各节点上所需要施加的节点力。令节点1在坐标x方向的位移u1=1，而其余的节点位移v1=u2=v2=u3=v3=un=vn=0，可得到节点载荷列阵等于K的第一列元素组成的列阵，即 第二、整体刚度矩阵中主对角元素总是正的。例如，整体刚度矩阵中的元素k33 是表示节点2在x方向产生单位位移，而其它位移均为零时，在节点2的x方向上必须施加的力，很显然，力的方向应该与位移方向一致，故应为正号。第三、整体刚度矩阵是一个对称矩阵，即。,4.2.2 边界条件的引入,第四、整体刚度矩阵是一个稀疏矩阵。如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线附近呈带状。如前所述，总刚中第r双行的子矩阵Krs，有很多位置上的元素都等于零，只有当第二个下标s等于r或者s与r同属于一个单元的节点号码时才不为零，这就说明，在第r双行中非零子矩阵的块数，应该等于节点r周围直接相邻的节点数目加一。可见，K的元素一般都不是填满的，而是呈稀疏状（带状）。第五、整体刚度矩阵是一个奇异矩阵，在排除刚体位移之后，它是一个正定矩阵。,4.2.2 边界条件的引入,2.边界条件的引入,方法一：保持方程组为2n2n阶不变，仅对K和R进行修正。例如，若指定节点i在方向y的位移为vi，则令K中的元素k2i,2i 为1，而第2i行和第2i列的其余元素都为零。R中的第2i个元素则用位移vi 的已知值代入，R中的其它各行元素均减去已知节点位移的指定值和原来K中该行的相应列元素的乘积。例如一个只有4个方程的简单例子。假定该系统中节点位移u1 和u2分别被指定为u1=1，u2=2,(4.45),(4.46),只有在消除了整体刚度矩阵奇异性之后，才能联立方程组并求解出节点位移。整体刚度矩阵的奇异性需要通过引入边界约束条件、消除结构的刚体位移来实现。,4.2.2 边界条件的引入,当引入这些节点的已知位移之后，方程(4.45)就变成 利用这组维数不变的方程来求解所有的节点位移，显然，其解仍为原方程(4.45)的解。如果在整体刚度矩阵、整体位移列阵和整体节点力列阵中对应去掉边界条件中位移为0的行和列，将会获得新的减少了阶数的矩阵，达到消除整体刚度矩阵奇异性的目的，这样处理与本方法在原理和最终结果等方面都是一致的。,(4,47),矩阵的降维,4.2.2 边界条件的引入,方法二：将整体刚度矩阵K中与指定的节点位移有关的主对角元素乘上一个大数，如1015，将R中的对应元素换成指定的节点位移值与该大数的乘积。实际上，这种方法就是使K中相应行的修正项远大于非修正项。把此方法用于上面的例子，则方程(a)就变成 该方程组的第一个方程为因 故有 以此类推。,(4.48),(4.49),(4.50),(4.51),例4-1 如图4-3所示由两个三角形单元所组成的结构系统，试分别采用直接组集法及转换矩阵法进行刚度矩阵的组集，并进行边界条件的引入构件新的有限元方程。图4-3 两个三角形单元组成的结构系统,图4-3 两个三角形单元组成的结构系统,4.2.2 边界条件的引入,解：采用直接组集法进行整体刚度矩阵的组集 整个系统中共有4个节点，每个节点有两个自由度，则整体刚度矩阵的维数为8。三角形单元刚度矩阵的维数为8*8，因而需要将每个单元扩展成的矩阵再进行叠加。用分块矩阵表示，扩展完成的单元（1）和单元（2）的刚度矩阵可表示为,(a),4.2.2 边界条件的引入,矩阵中各元素的下标表示整体节点编号，和 则表示单元（1）和单元（2）对整个系统刚度矩阵的贡献。两者相加则得到总刚矩阵，表示为,(b),4.2.2 边界条件的引入,采用转换矩阵法进行整体刚度矩阵的组集 采用转换矩阵法进行整体刚度矩阵组集的关键是获得每个单元的转换矩阵。单元转换矩阵的行数为单元自由度数，转换矩阵的列数为整体刚度矩阵的维数。对于上述结构，每个转换矩阵的维数为，两个单元的转换矩阵分别为,(c),4.2.2 边界条件的引入,获得了每个单元的转换矩阵，则可按 进行整体刚度矩阵的求解。边界条件的引入图4-3所示结构，引入约束前有限元方程可表示为,(d),4.2.2 边界条件的引入,边界条件分析：由于1,2节点为固定约束，位移边界条件为u1=0,v1=0,u2=0,v2=0，另外，在节点4的x轴负方向有作用力，对应在整体刚度矩阵、整体位移列阵和整体节点力列阵中对应去掉边界条件中位移为0的行和列，并引入节点作用力，则最终引入边界条件的有限元方程为。,(e),上述方程已消除整体刚度矩阵的奇异性，可直接进行求解,4.2 利用平面三角形单元进行整体分析,4.3 有限元法的实施步骤,（1）结构的离散化 将分析的对象划分为有限个单元体，并在单元上选定一定数量的点作为节点，各单元体之间仅在指定的节点处相连。单元的划分，通常需要考虑分析对象的结构形状和受载情况。为了提高有限元分析计算的效率、达到一定的精度，还应该注意以下几个方面的问题：首先，在划分单元之前，有必要先研究一下计算对象的对称或反对称的情况，以便确定是取整个结构，还是部分结构作为计算模型。,图4-4 结构的对称性利用,共包括以下7个步骤：,4.3 有限元法的实施步骤,此外，节点的布置是与单元的划分互相联系的。通常集中载荷的作用点、分布载荷强度的突变点、分布载荷与自由边界的分界点、支承点等都应该取为节点。并且，当结构是由不同的材料组成时，厚度不同或材料不同的部分，也应该划分为不同的单元。另外，节点的多少及其分布的疏密程度（即单元的大小），一般要根据所要求的计算精度等方面来综合考虑。在进行节点编号时，应该注意要尽量使同一单元的相邻节点的号码差尽可能地小，以便最大限度地缩小刚度矩阵的带宽，节省存储、提高计算效率。,（2）单元位移模式的选择 根据分块近似的思想，选择一个简单的函数来近似地构造每一单元内的近似解。例如，若以节点位移为基本未知量，为了能用节点位移表示单元体的位移、应变和应力，在分析求解时，必须对单元中位移的分布做出一定的假设，即选择一个简单的函数来近似地表示单元位移分量随坐标变化的分布规律，这种函数称为位移模式。位移模式的选择是有限单元法分析中的关键。由于多项式的数学运算比较简单、易于处理，所以通常是选用多项式作为位移模式。多项式的项数和阶数的选择一般要考虑单元的自由度和解答的收敛性要求等。,4.3 有限元法的实施步骤,（3）单元刚度分析 通过分析单元的力学特性，建立单元刚度矩阵。首先利用几何方程建立单元应变与节点位移的关系式，然后利用物理方程导出单元应力与节点位移的关系式，最后由虚功原理或最小势能原理推出作用于单元上的节点力与节点位移之间的关系式，即单元刚度矩阵。（4）等效节点力计算 分析对象经过离散化以后，单元之间仅通过节点进行力的传递，但实际上力是从单元的公共边界上传递的。为此，必须把作用在单元边界上的表面力，以及作用在单元上的体积力、集中力等，根据静力等效的原则全都移置到节点上，移置后的力称为等效节点力。,4.3 有限元法的实施步骤,(5)整体结构平衡方程建立 建立整体结构的平衡方程也叫做结构的整体分析，也就是把所有的单元刚度矩阵集合并形成一个整体刚度矩阵，同时还将作用于各单元的等效节点力向量组集成整体结构的节点载荷向量。从单元到整体的组集过程主要是依据两点：一是所有相邻的单元在公共节点处的位移相等，二是所有各节点必须满足平衡条件。在本书中，组集整体刚度矩阵的方法包括直接组集法与转换矩阵法。(6)引入边界约束条件 在上述组集整体刚度矩阵时，没有考虑整体结构的平衡条件，所以组集得到的整体刚度矩阵是一个奇异矩阵，尚不能对平衡方程直接进行求解。只有在引入边界约束条件、对所建立的平衡方程加以适当的修改之后。具体引入边界条件的方法可参照4.3.2。,4.3 有限元法的实施步骤,（7）求解未知的节点位移及单元应力 静态有限元分析的计算结果主要包括位移和应力两方面。位移已经获得，而对于应力计算结果则需要进行如下整理。以平面三角形单元为例加以说明。如前所述，平面三角形单元是常应变单元，也就是常应力单元。计算得到的单元应力通常视为单元形心处的应力。为了能根据计算结果推算出结构任一点处的应力值，一般采用绕节点平均法或两单元平均法进行处理。所谓的绕节点平均法，就是将环绕某一节点的各单元常应力加以平均，用以表示该节点的应力。为了使求得的应力能较好地表示节点处的实际应力，环绕该节点的各个单元的面积不应相差太大。一般而言，绕节点平均法计算出来的节点应力，在内节点处较好，而在边界节点处则可能很差。因此，边界节点处的应力不宜直接由单元应力平均来获得，而应该由内节点的应力进行推算。另一种推算节点应力值的方法是两单元平均法，即把两个相邻单元中的常应力加以平均，用来表示公共边界中点处的应力。这种情况下，两相邻单元的面积也不应相差太大。,4.3 有限元法的实施步骤,4.3 有限元法的实施步骤,求解弹性力学平面问题的实施步骤,以三角形常应变单元为例，应用有限元法求解弹性力学平面问题的步骤如下：（1）将计算对象进行离散化，即把结构划分为许多三角形单元，并对节点进行编号。确定全部节点的坐标值。（2）对单元进行编号，并列出各单元三个节点的节点号。（3）计算外载荷的等效节点力，列写结构节点载荷列阵。（4）计算各单元的常数b1、c1、b2、c2、b3、c3 及行列式2，计算单元刚度矩阵。（5）组集结构整体刚度矩阵。（6）引入边界条件，处理约束，消除刚体位移。（7）求解线性方程组，得到节点位移。（8）整理计算结果，计算应力矩阵，求得单元应力，并根据需要计算主应力和主方向。,如图4-5所示的平面应力问题，cm，单元厚度 mm，弹性模量 Pa，泊松比。求解各节点位移、单元应力、单元应变。图 4-5 平面应力状态结构,图4-5 平面应力状态结构,4.4 平面三角形单元举例,作为示例，该平板分成6个平面三角形单元、8个节点。有限元求解的具体过程如下：(1)建立平面直角坐标系oxy，原点设在节点1处，水平为x轴，垂直为y轴。列写节点坐标值及单元编码如下：（这里采用国际单位制进行了单位的统一）,表4-1 节点坐标值,4.4 平面三角形单元举例,表4-2 单元编码(2)计算各单元的单元刚度矩阵并进行扩展。首先计算出所需的系数。根据式（4.13），对于单元1，三个节点对应的整体编码为，求得,4.4 平面三角形单元举例,类似地求得应变矩阵为,4.4 平面三角形单元举例,弹性力学平面问题的弹性矩阵为得到单元1的单元刚度矩阵为,4.4 平面三角形单元举例,单元1的单元刚度矩阵进行扩展，得到一个1616阶的方阵。只在上述(1,2,3)三个节点对应的元素上有值，其它元素上均为0。全部6个单元均按上述同样的过程进行计算。（3）对上述6个单元的扩展刚度矩阵进行叠加，得到该结构的整体刚度矩阵为,4.4 平面三角形单元举例,(4)在考虑位移约束条件的情况下列写结构节点位移列阵。在本示例中，节点1，2处均为全约束，即这两个节点的x,y方向对应的位移分量为0，即(5)考虑结构的外载荷，构造结构载荷列阵。在本示例中，只在节点8处作用水平和垂直载荷，因此可以得到,4.4 平面三角形单元举例,(6)引入边界条件，即根据约束情况修正结构有限元方程，特别是消除整体刚度矩阵的奇异性，得到考虑约束条件的、可解的有限元方程。(7)利用线性方程组的数值解法，对上述结构的有限元方程进行求解，得到所有各节点的位移向量。最后根据需要求解单元应力。,4.4 平面三角形单元举例,程序：,4.4 平面三角形单元举例,%基本数据clearNJ=8%节点总数Ne=6%单元总数XY=.%节点坐标0 00.08 00 0.040.08 0.040 0.080.08 0.080 0.120.08 0.12,Code=.%单元编码1 2 33 2 43 4 55 4 65 6 77 6 8;E=2.06e11%材料参数Nu=0.3t=0.001%计算单元刚度矩阵D=E/(1-Nu*Nu)*1 Nu 0;Nu 1 0;0 0(1-Nu)/2;Kz=zeros(2*NJ,2*NJ);,4.4 平面三角形单元举例,for e=1:Ne I=Code(e,1);J=Code(e,2);M=Code(e,3);x1=XY(I,1);x2=XY(J,1);x3=XY(M,1);y1=XY(I,2);y2=XY(J,2);y3=XY(M,2);A=0.5*det(1 x1 y1;1 x2 y2;1 x3 y3);b1=y2-y3;b2=y4-y1;b3=y1-y2;c1=-(x2-x3);c2=x1-x3;c3=x2-x1;B=.b1 0 b2 0 b3 0,0 c1 0 c2 0 c3c1 b1 c2 b2 c3 b3/(2*A);Ke=t*A*B*D*B;%单元刚度矩阵的扩展与叠加Kz(2*I-1:2*I,2*I-1:2*I)=Kz(2*I-1:2*I,2*I-1:2*I)+Ke(1:2,1:2);Kz(2*I-1:2*I,2*J-1:2*J)=Kz(2*I-1:2*I,2*J-1:2*J)+Ke(1:2,3:4);Kz(2*I-1:2*I,2*M-1:2*M)=Kz(2*I-1:2*I,2*M-1:2*M)+Ke(1:2,5:6);%=Kz(2*J-1:2*J,2*I-1:2*I)=Kz(2*J-1:2*J,2*I-1:2*I)+Ke(3:4,1:2);Kz(2*J-1:2*J,2*J-1:2*J)=Kz(2*J-1:2*J,2*J-1:2*J)+Ke(3:4,3:4);Kz(2*J-1:2*J,2*M-1:2*M)=Kz(2*J-1:2*J,2*M-1:2*M)+Ke(3:4,5:6);%=Kz(2*M-1:2*M,2*I-1:2*I)=Kz(2*M-1:2*M,2*I-1:2*I)+Ke(5:6,1:2);Kz(2*M-1:2*M,2*J-1:2*J)=Kz(2*M-1:2*M,2*J-1:2*J)+Ke(5:6,3:4);Kz(2*M-1:2*M,2*M-1:2*M)=Kz(2*M-1:2*M,2*M-1:2*M)+Ke(5:6,5:6);end,4.4 平面三角形单元举例,Kz%Kz 整体刚度矩阵:16X16,NJXNJ 子矩阵组成%F=Kz*U 节点力列阵F都是外力，支反力不计，因为支座位移为0。%F=0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 100 50;F=zeros(2*NJ,1);F(15)=100;F(16)=50;%引入约束条件：u1=v1=0;u2=v2=0相当于Kz(1,:)=0;Kz(:,1)=0;Kz(1,1)=1;Kz(2,:)=0;Kz(:,2)=0;Kz(2,2)=1;Kz(3,:)=0;Kz(:,3)=0;Kz(3,3)=1;Kz(4,:)=0;Kz(:,4)=0;Kz(4,4)=1;Kz%新的总体刚度矩阵%新的载荷列阵F(1)=0;F(2)=0;F(3)=0;F(4)=0;F%求解节点位移U=inv(Kz)*F,4.4 平面三角形单元举例,%后处理，计算单元应变应力 Strain=;Stress=;for e=1:Ne I=Code(e,1);J=Code(e,2);M=Code(e,3);x1=XY(I,1);x2=XY(J,1);x3=XY(M,1);y1=XY(I,2);y2=XY(J,2);y3=XY(M,2);A=0.5*det(1 x1 y1;1 x2 y2;1 x3 y3);b1=y2-y3;b2=y4-y1;b3=y1-y2;c1=-(x2-x3);c2=x1-x3;c3=x2-x1;,B=.b1 0 b2 0 b3 00 c1 0 c2 0 c3c1 b1 c2 b2 c3 b3/(2*A);%把当前单元的节点位移从总体位移列阵中提取出来dlta=U(2*I-1),U(2*I),U(2*J-1),U(2*J),U(2*M-1),U(2*M);Strain_e=B*dlta;Stress_e=D*Strain_e;Strain=Strain Strain_e;Stress=Stress Stress_e;endStress%Sx Sy TxyStrain,4.4 平面三角形单元举例,(1)节点位移（U=u1,v1,u2,v2,.,u8,v8）U=1.0e-005*0 0 0 0 0.0825 0.0538 0.0906-0.0323 0.1966 0.0857 0.2150-0.0453 0.3075 0.0977 0.3720-0.0471,(2)单元应力单元 1 2 3 4 5 6 1.0e+006*0.8871-0.3098 0.7472 0.2883 0.7024 1.7442 2.9569-1.7069 1.8163-0.5663 0.8140 0.4360 1.5861 0.9139 1.3674 1.1326 0.8721 1.6279,(3)单元应变单元 1 2 3 4 5 6 1.0e-004*0 0.0101 0.0101 0.0229 0.0229 0.0807 0.1345-0.0807 0.0796-0.0326 0.0302-0.0044 0.2062 0.1188 0.1778 0.1472 0.1134 0.2116,4.4 平面三角形单元举例,作为对照，用ANSYS求解该示例。利用ANSYS求解的命令流文件如下：,/prep7ET,1,plane42,3!选单元!定义材料特性 MP,EX,1,2E11 MP,PRXY,1,0.3 type,1r,1,0.001!绘制模型 n,1,0,0n,2,0.08n,3,0,0.04n,4,0.08,0.04n,5,0,0.08n,6,0.08,0.08n,7,0,0.12n,8,0.08,0.12,e,1,2,3e,3,2,4e,3,4,5e,5,4,6e,5,6,7e,7,6,8d,1,all!加约束d,2,allf,8,fx,100!加力f,8,fy,50/solusolve/post1!进入后处理PRNSOL,dof!列表显示位移、应力、应变PRESOL,s,compPRESOL,epel,comp,4.4 平面三角形单元举例,ANSYS计算结果如下：(1)节点位移,4.4 平面三角形单元举例,ANSYS计算结果如下：(2)单元应力,4.4 平面三角形单元举例,ANSYS计算结果如下：(3)单元应变,