有限基础分析 2

FINITEELEMENTMETHODFUNDAMENTALS第4章有限元法基础从结构离散化到整体方程求解与应力分析工学院·有限元方法课程01基本概念有限元思想、基本概念、分析步骤与虚位移原理02离散与位移模式单元选择、网格离散与形函数03单元分析与等效节点载荷应变、应力、单元刚度与载荷等效04整体分析与解法整体组装、边界条件与有限元方程05应力结果分析由节点位移得到应变与应力结果06计算实例平面三角形单元的完整求解流程CONTENTS本章内容4.1基本概念建立有限元思想、基本术语与分析流程4.1基本概念：学习导引建立有限元思想、基本术语与分析流程核心思想将复杂结构离散为有限个容易分析的单元，单元之间通过有限个节点相互连接。基本概念节点自由度描述独立运动参数；约束限制自由度；载荷需等效到邻近节点。分析步骤有限元分析主要包括结构离散化、单元特性分析和整体分析。有限元简介有限元法的基本思想是将结构离散化，用有限个容易分析的单元来表示复杂的对象，单元之间通过有限个结点相互连接，然后根据变形协调条件综合求解。由于单元的数目是有限的，结点的数目也是有限的，所以称为有限元法(FEM，FiniteElementMethod)。有限元法是最重要的工程分析技术之一。它广泛应用于弹塑性力学、断裂力学、流体力学、热传导等领域。有限元法是60年代以来发展起来的新的数值计算方法，是计算机时代的产物。虽然有限元的概念早在40年代就有人提出，但由于当时计算机尚未出现，它并未受到人们的重视。随着计算机技术的发展，有限元法在各个工程领域中不断得到深入应用，现已遍及宇航工业、核工业、机电、化工、建筑、海洋等工业，是机械产品动、静、热特性分析的重要手段。早在70年代初期就有人给出结论：有限元法在产品结构设计中的应用，使机电产品设计产生革命性的变化，理论设计代替了经验类比设计。有限元法的孕育过程及诞生和发展牛顿(Newton)莱布尼茨(LeibnizG.W.)大约在300年前，牛顿和莱布尼茨发明了积分法，证明了该运算具有整体对局部的可加性。虽然，积分运算与有限元技术对定义域的划分是不同的，前者进行无限划分而后者进行有限划分，但积分运算为实现有限元技术准备好了一个理论基础。在牛顿之后约一百年，著名数学家高斯提出了加权余值法及线性代数方程组的解法。这两项成果的前者被用来将微分方程改写为积分表达式，后者被用来求解有限元法所得出的代数方程组。高斯(Gauss)拉格朗日(LagrangeJ.)在18世纪，另一位数学家拉格朗日提出泛函分析。泛函分析是将偏微分方程改写为积分表达式的另一途径。在19世纪末及20世纪初，数学家瑞利和里兹（RayleighRitz）首先提出可对全定义域运用展开函数来表达其上的未知函数。

瑞利(Rayleigh)1915年，数学家伽辽金(Galerkin)提出了选择展开函数中形函数的伽辽金法，该方法被广泛地用于有限元。1943年，数学家库朗德第一次提出了可在定义域内分片地使用展开函数来表达其上的未知函数。这实际上就是有限元的做法。（对象、变量、方程、求解途径）各力学学科分支的关系非变形体（刚体）变形体(1)桥梁隧道问题研究对象：任意形状的变形体任意变形体力学分析的基本变量及方程圆形隧道三维模型(2)中华和钟(3)矿山机械(4)压力容器的成形变形体及受力情况的描述思路：以计算机为工具，分析任意变形体以获得所有力学信息，并使得该方法能够普及、简单、高效、方便，一般人员可以使用。实现办法：有限元方法的思路及发展过程技术路线：发展过程：如何处理对象的离散化过程......点(质量)线(弹簧，梁，杆，间隙)面(薄壳,二维实体,轴对称实体)二次体(三维实体)线性二次..线性..............................常用单元的形状点单元线单元一维波传导问题点单元线单元面单元受垂直载荷的托架线性单元/二次单元更高阶的单元模拟曲面的精度就越高。低阶单元更高阶单元体单元复杂问题的建模简化与特征等效软件的操作技巧（单元、网格、算法参数控制）计算结果的评判二次开发工程问题的研究误差控制有限元分析的作用应力应变位移弹性力学各个量之间的关系平衡方程物理方程几何方程外力有限元基本思想只要位移场确定，就可得到应变、应力。应力应变位移平衡方程放弃物理方程几何方程外力能量原理有限元的基本思想：在弹性体内选取足够多、有限个点，假定这些点的位移已知，再用这些假定的位移量描述其它位置点的位移，就得到了用特定点位移表示的弹性体的位移场。这些选定的有代表性的点——结点，（node）结点：代表性——尖点、拐角、截面改变处等集中载荷作用、位移约束位置等。位移场：某个点（非结点）位移不是由所有结点位移来表述的，而是划分成小区域/小块上的结点来表示的，这些小区域/小块——单元。有限元处理问题的方法——连续体剖分小块（单元），即离散体。单元网格划分中每一个小的块体节点确定单元形状、单元之间相互联结的点节点力单元上节点处的结构内力载荷作用在单元节点上的外力（集中力、分布力）约束限制某些节点的某些自由度弹性模量（杨式模量）E泊松比（横向变形系数）μ密度单元单元载荷节点节点力约束有限元单元模型中几个重要概念有限单元法的分析步骤如下：物体离散化单元特性分析单元组集，整体分析求解未知节点的位移由节点的位移求解各单元的位移和应力有限元法特点：概念浅显，容易掌握，可以在不同程度上理解与应用通用性强，应用广泛，几乎所有领域；计算格式统一，便于编程计算；大型通用程序成熟商业化，无需专门知识编程先进的前处理，网格自动划分，完善的后处理，可视或动态显示，直观形象。误差难估计结构离散化

将结构分成有限个小的单元体，单元与单元、单元与边界之间通过节点连接。结构的离散化是有限元法分析的第一步，关系到计算精度和效率，包括以下三个方面：单元类型的选择。选定单元类型，确定单元形状、单元节点数、节点自由度数等。单元划分。网格划分越细，节点越多，计算结果越精确，但计算量越大。网格加密到一定程度后计算精度提高就不明显，对应应力变化平缓区域不必要细分网格。节点编码。

注意：有限元分析的结构已不是原有的物体或结构物，而是由同样材料、众多单元以一定方式连接成的离散物体。用有限元分析计算所获得的结果是近似的（满足工程要求即可）。有限元单元法分析步骤单元特性分析选择未知量模式选择节点位移作为基本未知量时，称为位移法；选节点力作为基本未知量时，称为力法；取一部分节点位移和一部分节点力作为未知量，称为混合法。分析单元力学性质根据单元材料性质、形状、尺寸、节点数目、位置等，找出单元节点力和节点位移关系式，应用几何方程和物理方程建立力和位移的方程式，从而导出单元刚度矩阵。计算等效节点力作用在单元边界上的表面力、体积力或集中力都需要等效地移到节点上去，即用等效力来替代所有作用在单元上的力。有限元单元法分析步骤整体分析集成整体节点载荷矢量F。结构离散化后，单元之间通过节点传递力，作用在单元边界上的表面力、体积力或集中力都需要等效地移到节点上去，形成等效节点载荷。将所有节点载荷按照整体节点编码顺序组集成整体节点载荷矢量。组成整体刚度矩阵K

，得到总体平衡方程：引进边界约束条件，解总体平衡方程求出节点位移。通过上述分析可以看出有限单元法的基本思想是“一分一合”，分是为了进行单元分析，合是为了对整体的结构进行综合分析。有限元单元法分析步骤有限单元法(FEM)中，为了简洁清晰地表示各个基本量以及它们之间的关系，也为了便于编制程序利用计算机进行计算，广泛采用矩阵表示和矩阵运算。平面问题中，物体受体力，可用体力列阵表示：物体受面力，可用面力列阵表示：3个应力分量的应力列阵表示3个形变分量的应变列阵表示2个位移分量的位移列阵表示有限单元法中基本量的矩阵表示几何方程的矩阵表示为：

物理方程矩阵表示为：

利用应力列阵和应变列阵得：

其中矩阵

只与弹性常数E及μ有关，称为平面问题的弹性矩阵。用u*和v*表示虚位移，用表示与该虚位移相应的虚应变。根据虚功方程：处于平衡状态的变形体，外力在虚位移上所做的虚功等于应力在虚应变上所做的虚功。对于厚度为t的薄板，虚功方程可用矩阵表示为：其中，分别为体力列阵，面力列阵和应力列阵。为虚位移列阵有限单元法中，作用于弹性体的各种外力常以作用于某些点的等效集中力来代替。在厚度为t的薄板上，设作用于i点的集中力沿x及y方向的分量为Fix,Fiy,作用于j点的力为Fjx,Fjy等。这些集中力以及它们相应的虚位移用列阵表示为：为虚应变列阵虚位移原理代入虚功方程，得：

上式为集中力作用下的虚功方程。集中力列虚位移列阵外力在虚位移上所做的功为：4.2离散与位移模式从连续体到有限自由度模型4.2离散与位移模式：学习导引从连续体到有限自由度模型离散化将连续体的无限自由度问题转化为离散体的有限自由度问题。单元选择单元类型需匹配结构维度、形状、节点数量及节点自由度。位移模式描述单元内部任意点位移分量与节点位移之间的关系，也称插值函数。单元的形式是多样的实体单元模型离散单元类型维数

主要应用杆单元2-D承受轴向力作用，平面桁架结构3-D空间桁架、网架等梁单元2-D主要承受横向载荷作用，即承受弯矩；也可横向+轴向3-D固体2-D平面应力、平面应变问题：三角形、四边形3-D一般三维（空间）问题：四面体、八面体板壳2-D主要承受横向载荷作用三角形、四边形3-D单元类型与作用杆单元梁单元二维单元线性单元二次单元三维单元线性单元二次单元板壳单元离散化应注意的问题：首要的问题是根据结构的几何特点、受力特征选择合理的单元形式。对称性的利用，在划分单元之前，有必要先研究一下计算对象的对称或反对称的情况，以便确定是取整个物体，还是部分物体作为计算模型。取四分之一作为计算模型（以平面三角形单元为例）共边：覆盖求解区域，单元间既不允许相互重叠，也不允许相互脱开；共点：任意三角形的顶点必须是相邻单元的顶点，不能为相邻单元的内点。边长接近：单元的边长尽可能接近，采用锐角三角形数目与精度兼顾：单元划分细，计算精度越高，但结点数增加，计算时间加长。单元大小过渡，应力梯度大的区域单元尺寸小，应力变化小的区域，单元可以划分大些。或在初步计算的基础上对于高应力区，在进一步细化网格，进行二次分析。适当简化。不可以可以较差较好节点编号顺序

在进行节点编号时，应该注意要尽量使同一单元的相邻节点的号码差尽可能地小，以便最大限度地缩小刚度矩阵的带宽，节省存储、提高计算效率。平面问题的半带宽为

B=2(d+1)单元编码节点i节点j节点m11232243…nijm节点坐标x坐标y10221.60…mxy代码实现clearformatshorteNE=;%单元个数

NP=;%节点个数E=;%弹性模量u=;%泊松比T=;%厚度NN=3;%单元节点数LN=[];%单元节点号相应为单元节点号按逆时针顺序输入

CD=[];%节点坐标（1）取三角形单元的节点位移为基本未知量：

其中，

称为单元的节点位移列阵；（2）应用插值公式，由单元节点位移求出单元的位移函数：

其中，N称为形函数矩阵；（3）应用几何方程，由单元的节点位移求出单元的应变：

其中，B是表示与之间关系的矩阵；三角形单元离散化结构分析步骤

其中，Fe

是单元的节点力，k称为单元刚度列阵； 对三角形板单元，节点力为：

（5）应用虚功方程，导出单元节点力与节点位移之间的关系。对右图中的i节点：节点对单元的作用力为节点力，作用于单元上。（4）应用物理方程，由单元的节点位移求出单元的应力： 其中，S称为应力转换矩阵； Fe是作用于单元的外力，此外，单元内部还作用有应力。根据虚功方程，从而得到节点力的公式：（7）列出各节点的平衡方程，组成整个结构的平衡方程组。由于节点i受有环绕节点的单元移置而来的节点载荷和节点力因而i节点的平衡方程为： （i=1,2,…,n）

（6）应用虚功方程，将单元中的外力载荷向节点移置，化为节点载荷（即求出单元的节点载荷）：

将（f）代入（h），整理得：

其中，K称为整体刚度矩阵，FL是整体节点载荷列阵，δ是整体节点位移列阵。在上述求解步骤中，（2）至（6）是针对每个单元进行的，称为单元分析；（7）是针对整个结构进行的称为整体分析。要求：i、j、m按逆时针排序单元的结点位移向量用来描述单元内各点位移变化规律的函数，称为位移模式位移模式三角形单元的位移模式假定为位移模式：位移模式矩阵表达：位移模式通式——单元内任一点的位移；——单元的结点位移向量；——单元的形函数矩阵。面积坐标形函数与面积坐标的关系形函数的性质三角形的面积位移模式反映了单元内任意一点的位移与结点位移之间的关系。是有限元计算精度的关键。有限单元法中，应力转换矩阵和劲度矩阵的建立以及载荷的移置等，都依赖于位移模式。代码实现fori=1:NE%对所有单元进行循环%计算当前单元的面积A=det([1CD(LN(i,1),1)CD(LN(i,1),2);1CD(LN(i,2),1)CD(LN(i,2),2);1CD(LN(i,3),1)CD(LN(i,3),2)])/2;end4.3单元分析与等效节点载荷建立单元力学方程与载荷等效关系4.3单元分析与等效节点载荷：学习导引建立单元力学方程与载荷等效关系单元分析通过建立单元的力学方程，确定单元刚度特性和受力状态。应变与应力由几何方程和节点位移得到应变，再由物理方程得到应力。载荷等效等效节点载荷与原外力在任意虚位移上所做的虚功相等。单元上任意一点的应变单元上任意一点的应力单元的能量单元刚度矩阵的性质单元分析单元上任意一点的应变几何方程或写成通式[B]矩阵叫做单元几何矩阵，反映了单元内任意一点的应变分量与结点位移之间的关系几何矩阵[B]中的每个元素，均为常数，它们由结点坐标确定。单元内任意一点P的应变分量与坐标（x，y）无关，说明单元中应变是常量。单元上任意一点的应力物理方程[D]——弹性矩阵平面应力问题平面应变问题[S]叫做应力矩阵代码实现fori=1:NE…%生成应变矩阵Bforj=0:2b(j+1)=CD(LN(i,(rem((j+1),3))+1),2)-CD(LN(i,(rem((j+2),3))+1),2);c(j+1)=-CD(LN(i,(rem((j+1),3))+1),1)+CD(LN(i,(rem((j+2),3))+1),1);endB=[b(1)0b(2)0b(3)0;0c(1)0c(2)0c(3);c(1)b(1)c(2)b(2)c(3)b(3)]/(2*A);B1(:,:,i)=B;…end%生成弹性矩阵D平面应力D=[1u0;u10;00(1-u)/2]*E/(1-u^2);%求应力矩阵S=D*B;ES=B'*S*T*A;%求解单元刚度矩阵1、单元的应变能一维问题应变能密度为平面问题应变能密度为

单元的能量

称为单元刚度矩阵，简称单刚，它反映了单元应变能与单元结点向量之间的关系。

(1)、体力势能(2)、面力势能(3)、集中力势能2、外力势能3、单元的总势能某单元的刚度矩阵，仔细看看，会发现该矩阵有哪些特点？单元刚度矩阵的性质单元刚度矩阵的性质1、对称性单元刚度矩阵是对称方阵，其元素都对称于主对角线。2、奇异性单元刚度矩阵中任意一行或列元素之和为零。其物理意义是在没有给单元施加任何约束时，单元可有刚体运动，位移不能唯一的确定。3、主对角线元素恒为正值主对角线元素是正值说明结点位移方向与施加结点荷载的方向是一致的。

4、单元刚度矩阵与单元位置无关单元刚度矩阵与单元位置无关，也就是单元在平移时，[K]e不变；单元结点排列顺序不同时，[K]e中元素大小不变，而排列顺序相应改变。

弹性体所受外力包括体积力、表面力、集中力。分别作用在弹性体内部、物体表面上、物体的一个点上。载荷列阵{R}，是由弹性体的全部单元的等效节点力集合而成，是将全部载荷转移到单元的节点上，它们的作用位置发生了变化——载荷移置。它们的作用效果是等效的，故称等效节点力向量{R}e

。各种载荷分别移置到节点上，再逐点加以合成求得单元的等效结点载荷。等效载荷1、体力等效结点载荷自重情况下：y0xijmy0xijm2、面力等效结点载荷代码实现NF=;%结点荷载个数F=[10-1];%结点力数组（受力结点编号,x方向,y方向），此例中只有一个受力点%生成荷载向量ASL(1:2*NP)=0;%总体荷载向量置零fori=1:NFASL((F(i,1)*2-1):F(i,1)*2)=F(i,2:3);end4.4整体分析与解法由单元方程组装整体有限元方程4.4整体分析与解法：学习导引由单元方程组装整体有限元方程整体组装将各个离散单元组合成整体结构，建立整体结构的力学方程。总刚度矩阵各单元刚度矩阵按节点总码扩充并叠加，形成结构总刚度矩阵。方程求解施加位移边界条件后，求解有限元方程获得节点位移。结构的结点位移向量假设弹性体被划分为N个单元和n个节点，整个弹性体的节点位移向量{

}2n×1整个弹性体的载荷列阵{R}2n×1整体分析矢量——有方向性，外力、应力，不能直接相加标量——没有方向，只有大小，可以相加。弹性体的能量是标量，可以直接相加。结构的总势能

单刚的扩充为了实现上述运算扩展——结构的总刚度矩阵——结构总的体力列阵

——结构总的面力列阵

——结构总的集中力列阵

结构的总势能

123421q组装总刚[K]的一般规则：1.

当[Krs]中r=s时，该点被哪几个单元所共有，则总刚子矩阵[krs]就是这几个单元的刚度矩阵子矩阵[krs]e的相加。2.

当[krs]中rs时，若rs边是组合体的内边，则总体刚度矩阵[krs]就是共用该边的两相邻单元单刚子矩阵[krs]e的相加。3.当[krs]中r和s不同属于任何单元时，则总体刚度矩阵[krs]=[0]。下面，我们考查一个组装总刚的实例：1.整体刚度矩阵及载荷列阵的组集

根据叠加原理，整体结构的各个刚度矩阵的元素显然是由有关单元的单元刚度矩阵的元素组集而成的，为了便于理解，现结合图5说明组集过程。整体刚度矩阵形成方法图中有两种编码：一是节点总码：1、2、3、4；二是节点局部码，是每个单元的三个节点按逆时针方向的顺序各自编码为1，2，3。图中两个单元的局部码与总码的对应关系为：单元e的刚度矩阵分块形式为：

单元①：1，2，3 1，2，3

单元②

：1，2，33，4，1或：单元①

：1，2，3 1，2，3

单元②

：1，2，31，3，4整体刚度矩阵分块形式为：其中每个子块是按照节点总码排列的。

通常，采用刚度集成法或直接刚度法来组集整体结构刚度矩阵。刚度集成法分两步进行。

第一步，把单元刚度矩阵扩大成单元的贡献矩阵，使单元刚度矩阵的四个子块按总体编号排列，空白处作零子块填充。

第二步，以单元②

为例，局部码1，2，3

对应于总码3，4，1，按照这个对应关系扩充后，可得出单元②

的贡献矩阵

总码1234

2 3 431

2

局部码用同样的方法可得单元1的贡献矩阵。

第三步，把各单元的贡献矩阵对应行和列的子块相叠加，即可得出整体结构的刚度矩阵，如（42）式。

在这里应该指出，整体刚度矩阵中每个子块为2*2阶矩阵，所以若整体结构分为n个节点，则整体刚度矩阵的阶数是。

总码1234

123（42）

123局部码

至于整体结构的节点载荷列阵的组集，只需将各单元的等效节点力列阵扩大成2n行的列阵，然后按各单元的节点位移分量的编号，对应相叠加即可⒈刚度矩阵[K]中每一列元素的物理意义为：欲使弹性体的某一节点在坐标轴方向发生单位位移，而其它节点都保持为零的变形状态，在各节点上所需要施加的节点力。⒉正定性，刚度矩阵[K]中主对角元素总是正的。⒊刚度矩阵[K]是一个对称矩阵，即[Krs]=[Ksr]T。⒋刚度矩阵[K]是一个稀疏矩阵。如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线附近呈带状。5.奇异性。刚度矩阵[K]是一个奇异矩阵，在排除刚体位移后，它是正定阵。整体刚度矩阵的性质半带存储半带宽B=（相邻节点号的最大差值D+1）*2代码实现AS=zeros(2*NP,2*NP);%生成特定大小总体刚度矩阵并置0

fori=1:NE…a=LN(i,:);%临时向量,用来记录当前单元的节点编号forj=1:3fork=1:3AS((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=AS((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ES(j*2-1:j*2,k*2-1:k*2);%根据节点编号对应关系将单元刚度分块叠加到总体刚度矩阵中endendend结构的总势能最小势能原理，对于线弹性体，某一变形可能位移状态为真实位移状态的必要和充分条件是，此位移状态的变形体势能取最小值。结构总势能泛函对结点位移的变分为0.结构有限元方程它是一个2n阶的线性代数方程组。因为该方程中[K]是结构的总刚度矩阵，{F}是外荷载列阵，都通过计算求得，因此可以根据有限元方程可以确定结点位移。位移边界条件的处理由于总体刚度矩阵是奇异的，物理意义是结构中存在刚体位移，不能直接求解。必须引入限制结构刚体位移的位移边界条件，即位移约束条件，消除总体刚度矩阵的奇异性，才能求解结构有限方程。位移边界条件是指结构的某些区域位移已知，对于离散体来说，位移约束条件是某些结点的位移分量受到限制，包括位置限制和方向限制两个方面。具体哪些结点受到限制，受限制结点哪个方向位移分量受到限制，要根据结构受力后变形特征来确定。

处理的方法，主要有三种：

降阶法（紧缩法）置大数法改1法

解法1.降阶法降阶法也称紧缩法或直接代入法，该法是将结构有限元方程中已知结点位移的自由度全部消去，得到一组降阶的修正方程，用以求解其它未知的结点位移。如果给定的位移均为零位移，则只需将总刚[K]、荷载列阵{F}中与该位移所对应的行和列全部划去即可。如果给定的位移不为零位移，也只保留了待定的结点位移作为未知量，但需对右端荷载列阵进行相应的修正。2.置大数法将结构总刚度矩阵中与被约束的位移分量相对应的主对角线元素赋予一个大数A，如取A=10e30或更大。再将右端荷载列阵对应的荷载值换成已知的位移值与该大数的乘积。设结点位移分量r为已知，则有限元方程变为：经过修改后第r个方程的为方程两边同时除以A，除第r项外，其余各项均为微小量可略去。3.对角元素改1法当给定的位移值为零时，将总刚中与之相对应主对角线元素改为1，相对应的行和列中其余所有元素改为0，荷载列阵对应的元素也改为0即可。代码实现NV=;%受约束边界点数FD=[];%约束信息（约束点，x约束，y约束）有约束为1，无约束为0%将约束信息加入总体刚度矩阵（对角元素改一法）fori=1:NVifFD(i,2)==1AS(:,(FD(i,1)*2-1))=0;%一列为零AS((FD(i,1)*2-1),:)=0;%一行为零AS((FD(i,1)*2-1),(FD(i,1)*2-1))=1;%对角元素为1end%生成单元刚度矩阵并组成总体刚度矩阵ifFD(i,3)==1AS(:,FD(i,1)*2)=0;%一列为零AS(FD(i,1)*2,:)=0;%一行为零AS(FD(i,1)*2,FD(i,1)*2)=1;%对角元素为1endend4.5应力结果分析从节点位移恢复应变场与应力场4.5应力结果分析：学习导引从节点位移恢复应变场与应力场节点位移有限元方程求解后，得到结构全部节点的位移。应变场利用几何矩阵与单元节点位移，可计算单元应变场。应力场利用应力矩阵得到单元应力场，并结合工程目的进行显示和评价。计算模型中：位移场已经确定，就可得到应变、应力。应力应变位移物理方程几何方程外力有限元方程应力结果有限元方程求解之后，得到了所有结点的位移，单元应力计算对每个单元循环。对于任一单元①根据结点i、j、m的实际编号，从结构结点位移向量中选出单元结点位移向量②计算单元的应变分量，③计算单元的应力分量：单元应力计算步骤以上分析得到了所有单元的应力分量，为了强度分析，进一步计算主应力或等效应力。主应力取“＋”号为最大应力，取“－”号为最小应力最大应力与x轴的夹角MISES应力由应力分量表示的三维MISES应力由主应力表示的三维MISES应力由应力分量表示的二维MISES应力应力显示x应力mises应力代码实现%求解应力ASD=AS\ASL';%计算节点位移向量ELEDISP(1:6)=0;%当前单元节点位移向量fori=1:NEforj=1:3ELEDISP(j*2-1:j*2)=ASD(LN(i,j)*2-1:LN(i,j)*2);%取出当前单元的节点位移向量endSTRESS(:,:,i)=D*B1(:,:,i)*ELEDISP';%求应力end确定根据工程实际情况确定问题的力学模型，并按一定比例绘制结构图、注明尺寸、载荷和约束情况等。将计算对象进行离散化，即弹性体划分为许多三角形单元，并对节点进行编号。确定全部节点的坐标值，对单元进行编号，并列出各单元三个节点的节点号。④计算载荷的等效节点力。③单元分析，由各单元的相关参数，计算单元的几何矩阵、刚度矩阵。有限元分析的实施步骤处理约束，消除刚体位移，求解线性方程组，得到节点位移。计算应力矩阵，求得单元应力，并根据需要计算主应力和主方向。整理计算结果（后处理部分）。组集整体刚度矩阵，即形成总刚的非零子矩阵。组装各单元的等效结点载荷，形成总的外载荷向量。4.6计算实例以正方形薄板演示完整有限元计算流程4.6计算实例：学习导引以正方形薄板演示完整有限元计算流程问题模型均质正方形薄板上下承受均匀拉力，属于平面应力问题。对称简化利用几何与载荷对称性，可取四分之一结构进行研究。离散求解将四分之一结构离散为两个三角形单元，依次完成单刚、组装、约束和应力计算。图1所示为一厚度t=1cm的均质正方形薄板，上下受均匀拉力q=106N/m，材料弹性模量为E，泊松比，不记自重，试用有限元法求其应力分量。123421xy2myxq=106N/m计算实例1解：１.力学模型的确定２.结构离散由于此结构长、宽远大于厚度，而载荷作用于板平面内，且沿板厚均匀分布，故可按平面应力问题处理。考虑到结构和载荷的对称性，可取结构的1/4来研究。该1/4结构被离散为两个三角形单元，节点编号，单元划分及取坐标如图2所示，其各节点的坐标值见表1。节点坐标１２３４xy００１０１１０１３.求单元的刚度矩阵计算单元的节点坐标差及单元面积单元面积单元１（i、j、m

1，2，3）2）组装单元的几何矩阵3）计算单元的应力矩阵弹性矩阵应力矩阵单刚计算先计算用到的常数单元的刚度矩阵中各个子矩阵单元1的刚度矩阵为:123123（i、j、m=1，2，3）单元2：若按i、j、m=3、4、1顺序，对应单元1的123排码时，则这两个单元刚度矩阵内容完全一样，故有:3413414）组集整体刚度矩阵由于[Krs]=[Ksr]T，又单元1和单元2的节点号按123对应341，则可得：按刚度集成法可得整体刚度矩阵为：所以组集的整体刚度矩阵为：5.引入约束条件，修改刚度方程并求解根据约束条件：u1

=v1=0；v2=0；u4=0和等效节点力列阵： ，并代入刚度方程： ，划去[K]中与0位移相对应的1，2，4，7的行和列，则刚度方程变为：求解上面方程组可得出节点位移为：所以先求出各单元的应力矩阵[S]1、[S]2，然后再求得各单元的应力分量：6.计算各单元应力矩阵，求出各单元应力单元应力可看作是单元形心处的应力值。代码实现functionk_ele=TriangleElementStiffness(E,miu,t,node_ele)%TriangleElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforaTriangle(CST)%elementwithmodulusofelasticityE,%Poission'sratiomiu,constantthicknesst,%node_elethenodecoordinateofelement.%Thesizeoftheelementstiffness%matrixis6x6.%---------nodecoordinate------x1=node_ele(1,1);y1=node_ele(1,2);x2=node_ele(2,1);y2=node_ele(2,2);x3=node_ele(3,1);y3=node_ele(3,2);%-------------------------------A=(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2;%单元面积a1=x2*y3-y2*x3;a2=y1*x3-x1*y3;a3=x1*y2-y1*x2;b1=y2-y3;b2=y3-y1;b3=y1-y2;c1=x3-x2;c2=x1-x3;c3=x2-x1;B=1/2/A*[b10b20b30;0c10c20c3;c1b1c2b2c3b3];D=E/(1-miu^2)*[1miu0;miu10;

00(1-miu)/2];%strss/strainmatrixforplanestress%D=E/(1-miu^2)*[1miu0;%miu10;%00(1-miu)/2];%strss/strainmatrixforplanestraink_ele=t*A*B'*D*B;%的单元刚度矩阵单元刚度子函数代码实现functionstr=TriangleElementStress(E,miu,node_ele,u1,p)%TriangleElementStressThisfunctionreturnstheelement%stressmatrixforaTriangle%elementwithmodulusofelasticityE,%Poission'sratiomiu,%node_elethenodecoordinateofelement，planestressorplanestrainoption.

%---------nodecoordinate------x1=node_ele(1,1);y1=node_ele(1,2);x2=node_ele(2,1);y2=node_ele(2,2);x3=node_ele(3,1);y3=node_ele(3,2);%-------------------------------A=(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2;%单元面积a1=x2*y3-y2*x3;a2=y1*x3-x1*y3;a3=x1*y2-y1*x2;b1=y2-y3;b2=y3-y1;b3=y1-y2;c1=x3-x2;c2=x1-x3;c3=x2-x1;B=1/2/A*[b10b20b30;0c10c20c3;c1b1c2b2c3b3];ifp==1D=E/(1-miu^2)*[1miu0;miu10;00(1-miu)/2];%strss/strainmatrixforplanestresselseifp==2D=E/(1+miu)/(1-2*miu)*[1-miumiu0;

miu1-miu0;00(1-2*miu)/2];%strss/strainmatrixforplanestrainendstr=D*B*u1;%单元应力

单元应力子函数代码实现functionk_t=assemTriangle(k_t,k_ele,node1,node2,node3)%assemTriangleThisfunctionassemblestheelementstiffness%matrixkoftheplaneTriangleelementwithnodes%iandjintotheglobalstiffnessmatrixK.%Thisfunctionreturnstheglobalstiffness%matrixKaftertheelementstiffnessmatrix%kisassembled.d(1:2)=2*node1-1:2*node1;d(3:4)=2*node2-1:2*node2;d(5:6)=2*node3-1:2*node3;forii=1:6forjj=1:6k_t(d(ii),d(jj))=k_t(d(ii),d(jj))+k_ele(ii,jj);endend整刚组装子函数代码实现主程序clcclear;node=[1000;2100;3110;4010;];%节点信息，第一列为节点编号，2~4列分别为x,y,z方向坐标ele=[1123;2341;];%单元信息，第一列为单元编号，后面各列为单元上的节点号码num_ele=size(ele,1);%单元数%---物理参数------------------E=2.1e11;%弹性模量MPat=0.01;%单元厚度mmiu=1/3;%泊松比q=1e6;%均布载荷N/m%-------------------------------n_ele=length(ele(:,1));%单元数%组装总体刚度矩阵dof=length(node(:,1))*2;%自由度数f=ones(dof,1)*1e8;%结构整体外载荷矩阵，整体坐标系下f_loc=zeros(6,1);%单元外载荷矩阵，局部坐标系下u=ones(dof,1)*1e6;%位移矩阵K=zeros(dof);%总体刚度矩阵stress=zeros(n_ele,1);%单元应力矩阵fori=1:n_elek_ele=TriangleElementStiffness(E,miu,t,node(ele(i,2:4),2:4));K=assemTriangle(K,k_ele,ele(i,2),ele(i,3),ele(i,4));end代码实现主程序%力边界条件f(6)=q/2;%3节点垂向力f(8)=q/2;%4节点垂向力f(3)=0;f(5)=0;%位移边界条件u(1)=0;u(2)=0;u(4)=0;u(7)=0;%求解未知自由