有限元分析基础

有限元分析根底有限元法概述在机械设计中，人们常常运用材料力学、结构力学等理论知识分析机械零构件的强度、刚度和稳定性问题。但对一些复杂的零构件，这种分析常常就必须对其受力状态和边界条件进行简化。否那么力学分析将无法进行。但这种简化的处理常常导致计算结果与实际相差甚远，有时甚至失去了分析的意义。所以过去设计经验和类比占有较大比重。因为这个原因，人们也常常在设计中选择较大的平安系数。如此也就造成所设计的机械结构整体尺寸和重量偏大，而局部薄弱环节强度和刚度又缺乏的设计缺陷。近年来，数值计算机在工程分析上的成功运用，产生了一门全新、高效的工程计算分析学科——有限元分析方法。该方法彻底改变了传统工程分析中的做法。使计算精度和计算领域大大改善。§1.1有限元方法的开展历史、现状和将来历史有限元法的起源应追溯到上世纪40年代〔20世纪40年代〕。1943年R.Courant从数学的角度提出了有限元法的根本观点。50年代中期在对飞机结构的分析中，诞生了结构分析的矩阵方法。1960年R.W.Clough在分析弹性力学平面问题时引入了“FiniteElementMethod”这一术语，从而标志着有限元法的思想在力学分析中的广泛推广。60、70年代计算机技术的开展，极大地促进了有限元法的开展。具体表现在：1〕由弹性力学的平面问题扩展到空间、板壳问题。2〕由静力平衡问题——稳定性和动力学分析问题。3〕由弹性问题——弹塑性、粘弹性等问题。现状现在有限元分析法的应用领域已经由开始时的固体力学，扩展到流体力学、传热学和电磁力学等多个传统的领域。已经形成了一种非常成熟的数值分析计算方法。大型的商业化有限元分析软件也是层出不穷，如：SAP系列的代表SAP2000〔StructureAnalysisProgram〕美国安世软件公司的ANSYS大型综合有限元分析软件美国航天航空局的NASTRAN系列软件除此以外，还有MASTER、ALGO、ABIQUES、ADINA、COSMOS等。将来有限元的开展方向最终将和CAD的开展相结合。运用“四个化”可以概括其今后的开展趋势。那就是：可视化、集成化、自动化和网络化。§1.2有限元法的特点机械零构件的受力分析方法总体说来分为解析法和数值法两大类。如大家学过的材料力学、结构力学等就是经典的解析力学分析方法。在这些解析力学方法中，弹性力学的分析方法在数学理论上是最为严谨的一种分析方法。其解题思路是：从静力、几何和物理三个方面综合考虑，建立描述弹性体的平衡、应力、应变和位移三者之间的微分方程，然后考虑边界条件，从而求出微分方程的解析解。其最大的有点就是，严密精确。缺点就是微分方程的求解困难，很多情况下，无法求解。数值方法是一种近似的计算方法。具体又分为“有限差分法”和“有限元法”。“有限差分法”是将得到的微分方程离散成近似的差分方程。通过对一系列离散的差分方程求解，得到最终的力学问题近似解。其优点就是：计算简单收敛性好。缺点是：计算程序无法标准化，在不能获得整个问题的微分方程时，该方法不能运用。由于其是将微分方程转为差分方程，所以它是一种数学近似。“有限元法”的根本思想就是“先分后合”或者“化整为零，又积零为整”。与有限差分不同，它是在力学模型上进行近似处理，也就是〔分块近似〕。具体做法：把连续体模型转为由有限个单元组成的离散体模型，离散体模型之间通过一些节点联系。对于每一个离散体个体选择简单的函数近似表示其中的物理变化规律〔如位移等〕，运用力学方程推导单元的平衡方程组，然后集合所有的方程组形成表征整体结构的方程组，引入边界条件，求取最后问题的解。优点：概念清晰、易于学习理解，适用性强，便于电算化。缺点：计算精度受单元划分的影响较大。§1.3有限元分析的一般过程为了能够了解有限元分析的全貌，我们就一个简单的例子，来分析一下有限元分析的三个过程：结构离散化、单元分析、整体分析。结构离散化在该阶段中，要完成把连续结构的力学模型转变为离散的力学模型。处理的好坏，直接影响到最后分析结果的正确与否、计算的精度和计算的效率。根据模型的传力特性和分析的目标，正确选择单元类型。通常单元分为：一维单元、二维单元和三维单元。所谓一维单元就是指所求物理量仅随一个坐标变量而变化的单元。如桁架、平面刚架和空间刚架单元。一维单元：杆单元、梁单元。二维单元：三角形单元、四边形单元〔平面类问题〕三维单元：四面体单元、六面体单元等〔空间问题〕计算精度和计算效率：取决于单元划分的形状、大小和分布状况。通常单元愈多、愈密集，计算精度愈高，但计算效率愈低。有限元分析工作就是要在精度和效率两者之间做到有机的统一。单元分析进行单元分析的目的是为了到处表征单元力学特性的“单元刚度矩阵”。一般说来该过程有三种方法：直接法。虚功原理法〔变分法〕。加权余数法。直接法概念浅显，易于理解物理含义。变分法需要泛函的数学知识，其推导过程具有严谨的数学概念。加权余数法适用于泛函不存在的应用范围。本教材将运用虚功原理方程结合弹性力学和材料力学中的知识来推求几种常见单元的单刚计算公式。现在先看一个简单的阶梯轴的轴向拉伸问题例：如下图的变截面直杆，受拉力P，运用有限元方法分析其变形。取任意单元，长度为l，面积为A，单元内任意一点的轴向位移和其位置坐标成正比，即u=a0+a1x其中a0、a1为待定系数。由于杆的两个端点节点1、2是单元上的点，所以它们应该满足上述方程。节点1，x1=0，∴u1=a0+a1×0=a0节点2，x2=l，∴u2=a0+a1×la1=(u2-u1)/l将求出的结果带入方程并整理，就得：式中：N1、N2是形函数[N]形函数矩阵{δ}e节点位移向量由位移与应变的关系知道：将上面推出的位移表达式代入，可得：上式中的[B]称为应变矩阵或几何矩阵。运用材料力学中的虎克定律，可以将应变和应力联系起来。单向应力状态的虎克定律为××其中[S]称为应力矩阵。利用虚功方程可以建立力与位移之间的关系，也就是单刚方程。在后面我们将会推导出它的一般形式如下：式中：{F}e为单元节点力向量，对我们这个例子应为[U1U2]T。[K]e为单元刚度矩阵。后面将推导出它的计算公式为[D]矩阵是弹性矩阵。对于一维单元来说，就是E。所以我们这儿讨论的例题：求得单刚矩阵，也就完成了单元分析。总结单刚矩阵推导的步骤，应该分为四步：假定单元内位移变化的近似规律，即选择位移模式。运用几何关系，推求位移与应变的关系。应用物理规律，把应变与应力联系起来。运用虚功方程的力与位移关系，求出单刚矩阵。单元分析是整个有限元分析的核心。不同的单元因为其力学特性不同，而具有不同的单元刚度矩阵，我们这本教材就是要学习几种常用单元矩阵的推导和计算。了解各种单元的力学特性，为以后选择单元类型打好根底。整体分析由各单元刚度矩阵组集成整个结构的总刚度矩阵。整个结构有三个节点，首先将单元刚度矩阵扩充为3X3的矩阵，移动各元素使之与单刚矩阵中的元素位置相对应，如下：然后直接相加。把各单元的节点力向量组集成总的节点载荷向量。根据边界条件，修改总刚度矩阵，获得总刚方程组。边界条件修改之前的总刚方程：修改以后〔采用置“0，1”法〕求解方程组，得出总的节点位移向量。解得的解是：有了节点位移，再回代到前面单元推导过程中的公式×和××，就可以求得每个单元的应变和应力了。从这个简单的例子，我们了解了有限元法求解力学问题三大步骤中的内容，想必很多同学会说，这样复杂，如果运用材料力学的知识，我还来得快些。但是大家不要忘记，有限元的计算很多都是编程完成，而且现在很多的商业软件都已经完成了很多的工作。我们学习有限元主要是了解它的原理，并对常见单元的力学特性有所了解，这样对于以后运用有限元起到帮助作用。所以下面章节的内容，就是围绕这个主题展开。要到达这个目的，我们还必须学习必要的弹性力学知识。对弹性力学知识的学习，也对我们以后把握问题的本质有帮助。平面问题平面问题在力学研究的课题中属弹性力学的范畴。该类问题不仅本身具有典型性，而且在机械零构件的分析中，也是应用得非常广泛。所以这类问题也称之为经典的力学问题。我们知道，实际的机械零构件都是具有三维空间尺寸的物体，理应作为三维对象处理，但是当物体的几何形状和受力状态处于某些特定的情况下，近似地简化为平面问题，不仅可以大大简化计算的工作量，而且其精度也完全能够满足所要求。如：直齿圆柱齿轮可在垂直与孔轴线的截平面内作平面应力分析就足以了解整个齿轮的受力状态；大坝的横断面可作平面应变分析来了解整个大坝受力情况等。本章是全书的重点，在这里不仅介绍弹性力学的根本知识。还将系统地讲解有限元的根本概念、原理和方法。是学习以后各章节的根底。§2.1外力、应力、应变和位移外力、应力、应变和位移的概念在材料力学中已经学习过，由于这些概念在弹性力学、有限元法中具有和在材料力学中不同的规定，弹性力学中的规定和有限元法是完全相同的，所以在这里我们将按照弹性力学的习惯表达方法把他们集中的加以阐述。外力外界作用在物体上的作用力，可以分为两大类：体积力分布在物体体积内的力。如：重力、惯性力和磁性力等。单位体积的体力在坐标轴上的分量X、Y、Z，称为体力分量。符号规定为沿坐标轴正向的为正，沿负向的为负。面力作用在物体外表的力。如轮压、水压等。它又可细分为集中力与分布力。面力在坐标轴上的投影，表示为X、Y、Z。符号沿正轴为正，负轴为负。应力弹性体受到外力作用后，内部产生的抵抗变形的内力。以弹性体中P点为定点的微单元体来考察。所谓微单元体，就是图中PA、PB、PC的边长分别为dx、dy和dz。以下简称这样的微单元体为微元体。微元体每个面上的应力都可以分解为三个应力分量。以图中红面为例，分别是σx、τxy、τxz。应力命名的规那么：以应力所在面垂直的坐标轴为第一个下标，应力指向为第二下标。如果下标相同就用一个下标表示。符号规定：正面上的应力与坐标轴同向为正，反之为负。负面上的应力与坐标轴反向为正。反之为负。所谓正面就是面的外法向与坐标轴同向为正。反之为负面。作用在两个相互垂直的面上，并且垂直与该两面交线的剪应力互等。即：τxy=τyx；τyz=τzy；τxz=τzx如此以来，代表P点应力状态的应力分量应有6个，它们是：位移任一点的位置移动用u、v、w表示它在坐标轴上的三个投影分量。符号规定：沿坐标轴正向为正，反之为负。应变弹性体内各点的位移在受力后一般是不相同的。各点之间距离的改变，从而使物体形状发生变化，即所谓的变形。而物体的形状总可以用它各局部的长度和角度表示。长度的改变称为正应变ε，角度的改变称为剪应变γ。以微元体三个棱边的线伸长和角度的变化，就分别有和6个应力分量相对应的6个应变分量，即：为与前面符号规定一致，这里对剪应变的符号规定如下：正应变伸长为正，缩短为负；剪应变使直角变小为正，变大为负。§2.2两类平面问题前面我们讲过，实际受力物体都是三维的空间物体，作用在其上的外力，通常也是一个空间力系，其应力分量、应变分量和位移也都是x、y、z三个变量的函数。但是当所考察的物体具有某种特殊的形状和特殊的受力状态时，就可以简化为平面问题处理。弹性力学中的平面问题有两类。平面应力问题当物体的长度与宽度尺寸，远大于其厚度〔高度〕尺寸，并且仅受有沿厚度方向均匀分布的、在长度和宽度平面内的力作用时，该物体就可以简化为弹性力学中的平面应力问题。我们分析以下其应力特征。当z=±t/2时，有σz=0、τzx=0、τzy=0。由于板较薄〔相对于长度和宽度尺寸〕，外力沿板厚又是均匀分布的，根据应力应连续的假定〔弹性力学中的根本假定〕，所以可以认为，整个板的各点均有σz=0、τzx=0、τzy=0。如此以来，描述空间问题的6个应力分量也就变为了3个，即而且这些应力分量仅是x、y两个变量的函数。平面应变问题当物体是一个很长、很长的柱形体，其横截面沿长度方向保持不变，物体承受平行于横截面且沿长度方向均匀分布的力时，该问题就可以简化为平面应变问题处理。分析其应力特征。假定其长度方向为无限长，那么任一横截面都可以看作是物体的对称面，如此那么有该面上的点都有w=0，也就是横截面上的所有点都不会发生Z方向的位移。由这一点可以推出也就有εz=0、τzx=0、τzy=0。和平面应力相比拟，平面应变是εz=0，那么是否也就有σz=0呢？可能有同学想σ=Eε，当然也就有σz=0，这是错误的。平面应变状态下σz≠0的。虽然不等于零，但它也不是一个独立的变量了，它由σx、σy的大小而决定。如此以来，独立的应力分量同平面应力问题一样也是3个：两类问题的比拟几何特征平面应力厚度＜＜长度、宽度平面应变厚度＞＞长度、宽度为便于说明可讲上述长度看作为厚度受力特征外力都必须在其面内且不沿厚度方向变化应力特征平面应力σz=0、τzx=0、τzy=0。εz≠0自由变形〔无约束〕平面应变σz≠0但不是自变量、τzx=0、τzy=0。εz=0总以上比拟可以看出，平面应力是真正的2维〔平面〕应力状态，而平面应变却不是，而是3维应力状态，只不过σz不是独立变量而是随横截面平面应力分量而定。独立变化的应力分量只有3个，类似于平面应力状态。§2.3平衡微分方程弹性力学求解问题是从静力学、几何学和物理学三方面综合考虑的。所以我们首先微元体应该满足的平衡条件——平衡微分方程。我们以平面问题为例推导，看看它应该具有什么形式。首先对平面问题的微元体进行受力分析图，如左所示。物体静力平衡的条件是：∑Fx(y)=0；∑M=0。先看∑Fx=0展开化简得同理可求得∑Fy=0满足得条件，由∑M=0，列出方程如下：化简后得：略去微量项，可得：。这就是前面所将的剪应力互等。对于平面应变问题，微元体的前后面还有正应力σz，不过它们是互等的。对于推导出来的结果，没有任何影响。所以平面问题的平衡微分方程就是：写成矩阵形式§2.4几何方程考察平衡微分方程，其中具有三个未知变量σx、σy、τxy，，而只有两个方程，方程具有无数个解。说明仅从静力学关系无法求解该方程。我们必须从其它方面寻求帮助。弹性体在受到外力后，会发生位移和形变，从几何上描述弹性体各点位移于应变之间的关系，就是弹性力学中的又一个重要方程——几何方程。仍然取截面的微元体ABCD，AB、CD边长为dx、dy，厚度为“1”位移u、v都是x、y的函数，即u(x,y)、v(x,y)，偏导数、表示位移分量u、v沿坐标轴x的变化率，偏导数、表示位移分量u、v沿坐标轴y的变化率，设A点的位移为u、v，那么B’点的位移就是：同理的D’点的位移分量由于α角在位移和形变很微小的情况下非常小，所以A’B’≈A’B”线段AB位移后的总伸长量为A’B’-AB=A’B”-AB=uB’-uA=-u=∴，同理可得剪应变由α、β两个角度组成由于，所以，同理可得∴综合以上几何方程，并将它们写成矩阵形式：由以上方程可以看出，当弹性体的位移分量确定以后，由几何方程可以完全确定应变，反过来，应变却不能完全确定弹性体的位移。这是因为物体产生位移的原因有两点：变形产生的位移。因运动产生的位移。因此弹性体有位移不一定有应变，有应变就一定有位移。§2.5物理方程描述弹性体内应力与应变关系的方程，我们称之为物理方程，也叫材料的本构方程。弹性力学通常研究的是各向同性材料，在三维应力状态下的应力应变关系。当弹性体处于小变形条件下，正应力只会引起微元体各棱边的伸长或缩短，而不会影响棱边之间角度的变化，剪应力只会引起角度的变化而不会引起各棱边的伸长或缩短。因此运用力的叠加原理、单向虎克定律和材料的横向效应〔泊松效应〕，我们就可以很容易的推导出材料在三向应力状态下的虎克定律，也就是通常所说的广义虎克定律。式中，E——材料线弹性模量G——材料剪切弹性模量μ——材料横向收缩系数，即泊松系数。三者不是独立的。具有以下关系：这些参数都是材料的固有属性系数，可以通过查材料手册获得。如：钢材的弹性模量E＝196～206GPa之间，通常取2.1×105MPa，μ＝0.24～0.28之间，也取为0.3进行计算，G＝79GPa。将以上空间问题的物理方程运用到平面问题，其形式如下：平面应力问题的物理方程前面分析，平面应力问题有σz=0、τzx=0、τzy=0。所以：从以上物理方程，也论证了我们前面说的εz≠0的结论，但由于它是由x和y方向应力产生的附加无约束变形，所以通常不予以考虑。在有限元分析中更多的是运用应变表示的应力关系，所以我们将上式变形一下：以上方程的矩阵表达形式为：简记为：式中：{σ}、{ε}为该问题的应力、应变向量。[D]为弹性矩阵。它是一个对称矩阵，且只与材料的弹性常数有关。平面应变问题的物理方程因为εz=0所以由空间物理方程的第三式得：，代入〔1〕、〔2〕式得同理变型为应变表示应力的形式：矩阵形式：也可简记为平面应变问题的弹性矩阵不同于平面应力问题的弹性矩阵，比拟可以发现只需将平面应力问题弹性矩阵[D]中的材料常数E换为E/(1-μ2)，μ换为μ/(1-μ)就得到了平面应变问题的弹性矩阵。其实弹性矩阵的这种转换方法，是弹性力学中将平面应力结果，转换到平面应变问题结论的一般方法。因为在两种平面问题的描述方程中〔平衡微分方程、几何方程和物理方程〕，只有物理方程是不同的。§2.6边界条件求解弹性力学问题实际就是在确定边界条件下，求解8个根本方程〔平面问题而言〕，以确定8个未知变量。所以从数学的角度看，就是求解偏微分方程的边值问题。边界条件的给出通常是各式各样的。大体可以分为三类：第一类边值问题给定物体的体力和面力条件，确定弹性体的应力场和位移场。此类问题边界以力的形式给出，所以也称为应力边界条件。我们可以来考察一下应力边界的一般形式：是在Sσ面上给出的力的分量。平面问题如左图所示，设阴影局部的微元体弧长为ds，厚度为单元厚度“1”，其法线与X轴的夹角为θ，由阴影局部微元体的平衡条件可以推出：化简后得：此即为平面问题应力边界方程。第二类边值问题给出弹性体的体力和物体外表各点得位移条件，确定弹性体得应力场和位移场。由于以位移给出得边界条件，所以也称为位移边界问题。一般得位移边界条件为：在Su面上。第三类边值问题给定弹性体得体力和一定边界上得面力，其余边界上的位移，确定其应力场和位移场。由于边界以力和位移两种形式给出，所以也称为混合边界问题。针对不同的边界条件，弹性力学求解的方法也有所不同。§2.7弹性力学的解题方法〔解析法〕应力法由于第一类问题的边界条件以应力形式给出，所以以应力作为根本的未知量的求解过程，就是人们通常所说的应力法。由于平衡方程中有三个未知量，而只有两个平衡微分方程，必须找出另外一个包含应力分量的方程，才能求得方程解。考虑到弹性体变形前是一个连续体，变形后也应是连续体的根本假设，所以要求微元体的变形一定要协调，才能使变形前、后，不会发生裂缝、重叠等现象。要使变形协调，就要研究几何方程。前面介绍的平面问题几何方程如下：分别对εx、εy、求y、x的二阶偏导，然后相加：上式说明三个应变分量之间应满足的连续性条件，我们称之为形变协调方程〔相容方程〕。通过物理方程，使上述的形变协调方程换成应力表示的形式，使之与平衡微分方程就构成了应力法中需要求解的方程组。具体我们来看看：利用物理方程消去相容方程中的形变分量〔以平面应力为例〕利用平衡微分方程，消去上述公式中的剪应力eq\o\ac(○,1)eq\o\ac(○,2)eq\o\ac(○,1)式对X求偏导，eq\o\ac(○,2)式对y求偏导，然后两者相加代入相容方程，化简对于平面应变而言，运用前面讲过的物理方程的转换方法，只需将上式中的μ代以μ/(1-μ)就可以了。最终求解的方程组平面应力问题平面应变问题三个微分方程，三个未知变量，再考虑边界条件，即可求得。如果是单连体〔只具有唯一的封闭边界〕的对象，满足了以上方程组后就是实际的解。但对于多连体〔具有多个封闭边界〕的对象，中还包含有待定系数，这些待定系数会导致位移的解出现多值性。所以对于多连体的问题，还应考虑位移的单值条件，才能最终确定。该局部的内容可以参见徐芝纶编的《简明弹性力学教程》中圆环受均布压应力的情况〔P87〕位移法位移法主要针对第二类边界条件问题求解。解题步骤：1〕改写物理方程使之成为应变表示应力的形式2〕应用几何方程表示以上得到公式中的应变3〕将它们代入平衡微分方程经整理最后得到的位移法求解平面应力问题方程为：两个未知量，两个方程，再加以边界条件即可求得问题的解。以上介绍的解析法中，应力法和位移法是求解弹性力学问题的根本方法。但都需要解联立的偏微分方程组。求解过程中的数学难度，常常导致这种求解是无法进行的。由于应力法在体力为常量的情况下可以进一步简化为求解一个单独的微分方程的问题，所以应力法在解析法中相对应用较多。但即使这样，在应力法中，也常常采用逆解法或半逆解法。§2.8常体力情况下应力法的简化、应力函数及实例分析我们前面讲述了弹性力学的三大方程，及应用这三大方程的应力法和位移法解题步骤。但是也说了要求解这些联立的偏微分方程在数学上是存在很大难度的。很多情况下，根本无法进行。那么弹性力学如何在实际中进行应用，它们和我们前面学过的材料力学区别究竟在哪里？我们将通过这一节的学习，一方面了解如何应用这些弹性力学的方程求解问题，另一方面加深对力学概念的理解，建立力学分析问题的直观感觉，为建立有限元模型打好根底。我们知道在大多数情况下，我们分析的对象，体力是常数，它不随x、y坐标变化。如此以来，前面讲解的第三个方程〔应力表示的相容方程〕，就可以简化为了：简记为：以上方程称为拉普拉斯微分方程，数学上也称之为调和方程，满足调和方程的函数称之为调和函数，及这里的。是拉普拉斯算子。这样以来常体力情况下的应力法方程就是：以上方程都不含有材料常数E、μ，所以平面应力和平面应变两类问题具有相同的方程，这说明：在单连体问题中，只要边界相同、受同样的分布外力，应力分布与材料无关；也与是平面应力还是平面应变的状态无关。以上结论的意义：eq\o\ac(○,1)弹性力学平面解答的应用范围加宽。eq\o\ac(○,2)为实验应力分析提供了理论依据〔光弹实验〕下面我们考察平衡方程：其解由其次微分方程的通解，加上任意一组特解组成。特解我们可以很容易找到。如：所以现在关键是找其次方程的通解。由第一个方程，可得：，由数学微分理论，该式是一个函数全微分的充要条件。所谓全微分就是有一个函数且同理由第二式可得：由剪应力公式又知存在一个函数φ，可以使∴故：由于应力与函数φ存在这样的关系，因此函数φ即是应力函数。我们用应力函数来表示相容方程：上式说明φ是重调和函数。前面讲过在弹性力学中，常用逆解法和半逆解法。所谓逆解法就是设定各种满足相容方程的应力函数，运用σx、σy与φ的关系，求得应力分量，再考察其满足何种边界条件，从而知晓这样的应力函数可以解决什么问题。所谓半逆解法就是根据弹性体的边界形状和受力关系，设定局部应力分量为何种形式的函数，从而确定应力函数φ，在运用应力函数求出所有的应力分量，根据边界条件确定应力分量应具有的最终形式。下面我们来看一个半逆解法的例子。运用逆解法求简支梁受均布载荷的应力分布。由材料力学知，弯曲应力主要由弯矩M引起，剪应力由剪力引起，而挤压应力由分布载荷q引起。现在q为不随x变化的常量。因此我们设σy不随x坐标变化，即，因此我们对x积分：上式中f1(y)、f2(y)是待定函数。由于应力函数必须满足相容方程，所以：eq\o\ac(○,1)代入到式eq\o\ac(○,1)中考察上式可以看出它是一个x的二次方程，所以一般情况下只有两个根。也就是说只有两个位置能够满足上式。但我们对相容方程的要求是绝对满足。也就是要求在整个梁的范围内都满足。所以只有该方程的系数项和自由项全部为零。即：∴eq\o\ac(○,2)eq\o\ac(○,3)eq\o\ac(○,4)公式中的A.、B、C…K都是待定系数。公式eq\o\ac(○,3)、eq\o\ac(○,4)中分别省掉了常数项和一次项、常数项。这是由于f1(y)和f2(y)分别是应力函数中x的一次项和常数项的原因，这样处理不会对应力分量产生影响。最后求出的应力函数为：由应力与应力函数的关系，可以求出各个应力分量：由于以上求得的应力分量满足了平衡方程和相容方程，所以只需根据边界确定A…K的系数，就求得了该问题的解。根据对称性，知道为偶函数，为奇函数，所以有E=F=G=0通常梁的跨度远大于梁的深度，因此上下边界是主要边界，它们必须满足。将它们代入的表达式，并且考虑E=F=G=0以上四个方程解四个未知数，求得：B=0将他们代回到应力分量的表达式中，也就有：左右两个边界，由于前面已经考虑了对称性，所以这个仅考虑优边界。没有水平力。要x=l时，σx=0，考察σx的表达式，除非q=0。而这和条件相违背。所以在这个边界上我们只能要求局部满足。运用圣维南原理运用等效力系代替它。〔这样产生的误差只在力作用点附近较大〕。运用的等效力系就是合成力系为平衡力系：合力等于0合力矩等于0由第一个条件得K=0，〔奇函数在对称区间上的积分为零〕；由第二个条件得可以证明剪应力的合力为-ql。即最终求得的结果，加以整理：由于厚度为“1”，此时其惯性矩，静矩〔计算见图〕任意一点的弯矩剪力所以上式中的应力分量可以改写为：各项应力的分布σx第一项为主应力项，与材料力学中的结果完全一致。第二项为应力修正项。当L/h＞4时，仅占主项的1/60；当L/h=2时，仅占主项的1/15。所以对于深梁的工程构件不容无视第二项的影响。§2.8虚功方程从前一节深梁的例子，我们可以看到，弹性力学解析求解的过程是非常复杂的。这样的求解对实际工程情况说来，很多情况根本是不可能的。所以长期以来，技术人员就一直探求数值求解的方法。有限元法是其中最成功的方法。为分析单元特性和简化分析过程，我们还需了解单元的能量关系。因为在力学上，很多时候从能量的角度分析，可以大大简化分析的过程。应变能的概念有材料力学知，弹性体在受到外力作用发生变形的过程中，弹性体内部会存储变形势能——应变能。在单向应力场中，单位体积的应变能的计算可以表示为：对于平面问题，有三个应力分量和与之对应的应变分量。由于在小变形情况下，正交力系互不影响，由力的叠加原理，所以该种情况的应变能为：其中：整个弹性体的应变能：上式也表示应力在应变上所做的总功。虚功原理和虚功方程在理论力学中，我们曾经学到一个虚功原理，也称虚位移原理。其根本思想就是：假定加在物体上一个可能的、任意的、微小的位移，在平衡条件下，物体的约束反力所做虚功应等于外力所做虚功。因为能量必须守恒。在这里所说的可能的虚位移是指位移必须满足的约束条件。任意的是指位移的方向和类型是随意的。把这一原理运用到现在的弹性体中，衡量弹性体应满足的平衡能量关系就是：假定加在弹性体上一个可能的、任意的、微小的位移，在平衡条件下，弹性体内的应变能应等于外力所做虚功。同样是因为能量必须守恒。运用这一原理，我们可以推到有限元中广泛用到的虚功方程。假定弹性体发生、的虚位移，那么由几何方程得：现考察弹性体微元体和边界处微元体上的力所做虚功：eq\o\ac(○,1)内部微元体上的力所做虚功左面的应力虚功右面的应力虚功左、右两面上的虚功之和〔略去高阶微量，并考虑〕同理得剪应力的虚功之和体力X的虚功同样的考虑Y方向的以及以及体力Y的虚功，然后叠加成内部微元体上的虚功如下：eq\o\ac(○,2)边界上的微元体设斜边中点处的虚位移、，形心处的应力为、、那么在直角边上的应力和位移均有一个负的增量，如下图，虚功计算为：〔略去了高阶微量〕同理的dy直角边上的剪应力虚功为：代入已有的几何关系：，斜面上面力的虚功体积力的虚功同样的方法求得另一面上正应力与剪应力的虚功，全部相加即得斜边微元体上的虚功之和为：支反力处的虚位移为零，所以支反力不做功，将dw1+dw2，并对整个体积积分，可以得到整个弹性体内的总虚功：根据平衡微分方程和静力边界条件，上式的第一、第二项都是零，所以弹性体的总虚功为：根据能量守恒，它应与外力的虚功相等由于该等式引入了平衡方程和边界方程，所以上式虚功方程等价于静力平衡条件〔内部和边界微元体〕。不同之处在于它是一种能量的表示形式。为了便于有限元中方便运用，引入广义力和物理方程，虚功方程变形为：综合以上推到过程，虚功方程表达的物理概念就是：“假设弹性体处于平衡状态，那么外力在虚位移上作的虚功应等于应力在应变上作的虚应变功，或者说等于虚应变能”。§2.9平面问题单元划分有限元法在平面问题进行分析时，才采用三角形单元和四边形单元、或者矩形单元，三角形单元的优点是简单且对结构的不规那么边界逼近好，而矩形单元却更能反映实际弹性体内部的应力应变变化。这两点我们会逐渐向大家说明。所以一般说来，有限元分析，单元划分的密度和单元种类选取，对计算结果起重要作用。一般单元划分越密集，结果越精确。单元多也导致求解的线性方程组阶数增高，要求计算机的内存也更大，计算的时间也越长，分析的效率就越低。解决这一矛盾的方法就是在应力集中区域单元划分密集一些，应力变化梯度小的位置，划稀疏些，这样就能兼顾精度与效率的关系。一般的原那么是：根据结构的受力和支承特点，按对称和反对称的性质，简化分析模型，以减少计算分析的规模。合理布局单元的密集程度，以使计算结果精度高而计算量小。在同一单元内，单元的特性数据和材质数据应保持一致。集中载荷的作用点和载荷密度突变处应有节点。在欲知道应力状态、内力情况和位移值的位置应有节点。单元的选取欲分析的目标密切相关。模型的单元划分好后，把所有的单元和节点按一定的规律和顺序进行编号，选择适当的坐标系〔直角、柱面和球面〕，以方便确定各节点的坐标值。§2.10节点位移、节点力和节点载荷弹性体在承受外力作用后，其内力的传递实际是通过单元之间的边界来实现的。但我们把结构离散化后，如果单元划分得足够小时，可以看成为其内力的传递通过单元与单元之间的节点进行传递。对于平面问题而言，每个节点都有位移和力两个未知量，这两个量又都是x、y的函数，注意平面问题的节点是不能传递力矩的，为什么？节点位移对三节点三角形单元而言，因有三个节点，每个节点的位移都有x，y两个分量，所以一共有6个自由度。单元节点位移向量可表示为：节点力所谓节点力，就是单元对节点或节点对单元作用的力，它是弹性体内部的作用力，也就是我们常说的内力。和上面相同的道理，节点力向量为：节点载荷节点载荷是指作用在节点上的外力，包括：eq\o\ac(○,1)直接作用在节点上的外力eq\o\ac(○,2)经等效处理后，移植到节点上的等效力。可用Xi、Yi.表示。由力平衡条件知，节点要保持平衡，那么作用在节点上的载荷应等于节点内力的合力。即：，，所有的节点载荷向量表示为§2.11三节点三角形单元的位移模式和形函数弹性结构受外载荷后，内部各点的位移变化规律，一般都是很复杂的。很难找到一个函数来描述整个结构内部各点的位移变化规律。但当把整个结构离散以后，在一个相当小的单元内部，却可以用简单的函数来近似描述单元内部位移的这种变化规律。而不致造成很大的误差。这就好比一条复杂的曲线，用一个函数很难描述，但在把这条曲线分段以后，对于每一条分段曲线，却可以用直线或抛物线来描述。三角形单元的位移模式设单元内任意一点的位移分量u、v是其位置坐标x、y的线性函数，那么：a1…a6是待定系数。改写方程组成为矩阵的形式：单元的三个顶点i、j、m的坐标，分别为、和，因为它们也是单元上的点，所以应该满足以上假定的位移变化规律。代入上式：解以上6个方程，求得6个待定系数。同理得顺序轮换，A是三角形的面积在单元节点的顺序号i,j,m必须是按逆时针排列，否那么系数行列式是负值，而三角形的面积为负值，是不合理的。求得的6个系数可以用以下矩阵表示：形函数将所求得的6个待定系数代入位移模式表达式中：令顺序轮换就有上式就是假定位移模式下导出的单元内任一点位移表达式。该式的数学意义就是单元内任一点的位移可以由单元节点的某种形式插值得到，其中的插值基函数就是Ni、Nj、Nm。对于我们目前假定的位移模式是线性函数，所以得出的插值基函数也是类似的线性函数。由此可以看出，插值基函数具有反映单元位移变化形态的特征，所以也称之为位移形态函数，简称形函数。[N]就是形函数矩阵。形函数的性质eq\o\ac(○,1)单元内任一点的三个形函数之和恒等于1。eq\o\ac(○,2)在单元的三个顶点处，有i节点处j节点处m节点处以上这些，可以通过简单的数学运算进行证明。位移模式收敛性的分析由于位移模式的选取是人为假定的，这种假定只能近似模拟单元内位移的变化规律，由于单元刚度矩阵的推导是以假定的位移模式展开的，那么这种假定的位移模式能否使有限元数值解收敛与精确解，在很大程度上就取决于所选的位移模式，通过数学证明，可以找出位移模式满足收敛性的几个条件：A〕完备性位移模式必须包含单元的常应变状态。将u、v的表达式代入几何方程，得：因为系数a2…a6都是常数，所以上面的应变分量也是常量。这也说明所选的位移模式中包含有弹性体的常应变状态。在上面的表达式中不含x、y的变量，说明单元的应变是常量，这也说明这种单元是一种常应力单元。位移模式必须包含单元的刚体位移。弹性体位移一般包含两个局部，即变形位移和没有弹性变形的刚体位移。那么作为模拟单元位移状态的位移模式，也就应该能同时反映这两局部的位移。在eq\o\ac(○,1)中已经证明了弹性应变的变形，下名说明他也包含有刚体位移的特征。改写位移表达式如下：当时，由eq\o\ac(○,1)知，，所以上式：*我们再来看看一个点作刚体位移的运动方程。点M先转到M1，再由M1移到M2，如左图所示：**比拟上面的*和**式，可以看出，，由此得出在三角形线性位移模式中，也包含了单元的刚体位移。B〕协调性位移模式必须能够反映位移的连续性。在弹性力学求解问题时，曾经讲到过变形协调方程，也说过它是弹性体变形后仍保持连续和不发生撕裂、侵入缺陷的条件。那么位移模式的选取，也应该保证单元之间不出现撕裂和侵入的缺陷。由于上面假定的位移模式是线性函数，而两点就能决定一条直线。由于相邻单元的公共边界上的位移由与之相连的两个节点插值获得，而相邻的单元具有两个公共的节点，所以通过这两个节点所得的插值值，不可能出现不同。也就是不可能出现以下图a和图b的情况。以上三个条件是选取位移模式必须考虑的。完备形条件是收敛的必要条件；协调条件是充分条件。在有限元中，满足完备条件的的单元是完备单元，满足协调条件的是协调单元。§2.12〔三节点三角形〕单元刚度矩阵的推导上一节我们已经建立了三角形单元的位移插值模式，并求得了形函数的方程，这样就完成了单元内任一点的位移由单元节点位移表示〔插值〕的工作，接下来运用我们已学过的一系列知识，我们就可以完成单元刚度矩阵的推导了。推导过程由位移插值函数导出单元应变的单元节点位移表达式将代入上式，可得：如按分块矩阵记忆，那么：其中i=i,j,m矩阵[B]称为应变矩阵，或称为几何矩阵。由以上计算公式知，它与单元的节点坐标有关，但不随点的坐标变化，就是说在这一单元内所有点的应变是相同的。求得单元应力的单元节点位移的表达式将几何矩阵代入单元的物理方程，就有：弹性矩阵[D]是由材料常数组成的矩阵。令[S]=[D][B]，代入平面应力的物理方程，就有∴也可以写成分块矩阵的形式i=i,j,m用虚功方程导出单元刚度矩阵〔单刚矩阵〕虚功方程假定单元的厚度为t，上式改写为单元的虚功方程形式，虚应变也可以用几何方程表示代入上式由于虚位移元素是常量，所以可以提到积分号以外，并与左边的消去〔为什么？〕。于是上式变为：令，虚功方程就成为了单刚方程由于[B]、[D]都是不含有x，y的常数矩阵，所以双重积分实际就是对面积积分了。A——是三角形面积将前面求得的应变矩阵和弹性矩阵代入，然后作矩阵乘法。就得到我们要求得的矩阵计算公式。我们在这里采用分块矩阵的方法记忆。r,s=i,j,m注意以上是平面应力状态的单刚矩阵，如果是平面应变问题呢？单元刚度矩阵的性质单元刚度矩阵[K]e表示单元抵抗变形的能力。它与通常的弹簧刚度系数k的物理意义本质相同，只不过[K]e是一个6×6阶的矩阵。共有36个元素。这是因为三角形的节点力向量和节点位移向量都为6的缘故。物理意义分块矩阵[Kij]表示的物理含义是：节点j处产生单位位移，而节点i、m被约束，此时在节点i处产生的节点力。我们写出它的4个分量元素：在上面的矩阵中，元素的第一个下标表示产生节点力的节点，第二个下标表示产生节点位移的节点。上标“1”表示水平分量，“2”表示竖直分量。而且上标和下标的关系是对应的。也就是说第一个下标对应第一个上标；第二个下标对应第二个上标。如此就有：就表示节点j处产生竖直单位位移，在节点i处产生的水平方向的节点力。单元刚度矩阵是对称矩阵→这一点可以通过简单的数学证明如下：单元矩阵的对称性，从物理学角度反映出的道理就是，“功的互等”。也就是在节点j处产生某一位移引起节点i处的节点力，应等于在节点i处产生相同位移引起节点j处的节点力。单元刚度矩阵中的元素只与单元的材料性质、几何形状、尺寸大小有关，而与单元的位置无关。单元刚度矩阵中不含有ai、aj、am，上节对位移模式收敛性分析中，曾经说明了a1、a4分别表示单元的平移分量，而由上式知单元的平移运动与ai、aj、am有关，而[K]e又与ai、aj、am无关，所以说它不随坐标轴的平动而变。单元刚度矩阵是奇异矩阵，即从力学的角度理解单元的刚度方程，当给定位移时，可以求得力；当给定力时，却不能求得位移。因为[K]e不存在逆矩阵，在单元没有给出任何约束的情况下，除有应变的可能性外，还同时有刚体位移的可能性。所以方程无解。§2.13载荷的节点移置前面对于有限元模型的分析时，曾经说过，单元之间的力传递是通过节点进行的。所以不在节点上的力，必须按静力等效原那么，把它们移置到节点上。静力等效原那么：原载荷在任何虚位移上所做的虚功，与移置到节点上的节点载荷所做的虚功相等。这种处理方法，和我们前面讲到的圣维南原理相同。它们只会影响模型局部的应力分布，而不会影响整个结构的应力。下面根据力的类型，分类说明处理的方法。集中力的移置如下图的受力示意图，M处的力为，移置后的节点载荷为，虚功相等就是：∵，所以有∴由于我们现在一直局限于单元的载荷移置，所以上式中的{R}应记为{R}e。如果将它写成分量的形式：从上面公式可以清楚的看到，载荷移置结果与单元位移模式密切相关。体力移置体积力密度为，将单元中微元体的体力看作集中力，那么面力的移置设在单元的一个边界上作用有分布力，面力的密度为，将微面积tds上的面力也看作是集中力，那么有以上面力，体力的积分运算较为复杂，在线性模式下，可按照理论力学中静力学平行力分解的原理，直接求得等效的节点载荷：eq\o\ac(○,1)均质等厚的三角形单元，受W的体力，那么三个节点分别受到的节点载荷。eq\o\ac(○,2)i、j边受到集度为q的均布面力载荷，那么i、j边节点受到的等效节点载荷。eq\o\ac(○,3)i、j边受到线性分布的面力，i处集度为零，j处集度为q其合力为，那么i处的节点载荷为，j处的载荷为。这些简单的规那么，对于今后实际中的应用，可以提高效率。§2.14整体分析该过程是将离散别离的各个单元组集成离散的结构物，从而建立模型的总刚方程。总刚方程和总刚矩阵的组集总刚度矩阵的组集原那么整个离散结构变形后，各个单元在节点处仍然协调地相互连接。即环绕某个节点的n个单元，在节点i处具有相同的位移。数学公式的描述：各个节点应满足静力平衡条件。即每个节点上的节点力合力应等于该节点的节点载荷。数学公式描述：此处的代表环绕节点i的所有单元节点力求和。在该原那么指导下，实例的组集过程如图a所示的平面问题，采用图b的方法划分网格，单元分析完成后，现将它们组集成原问题的有限元离散网格，求该问题的总刚矩阵。单元eq\o\ac(○,1)节点号为i，j，m=1，2，3〔i节点从最小号开始，然后逆时针排列节点号〕，单刚方程：将它们展开完全相同的道理，得其他三个单元展开的单刚方程单元eq\o\ac(○,2)节点号为i，j，m=1，3，4单元eq\o\ac(○,3)节点号为i，j，m=3，5，4单元eq\o\ac(○,4)节点号为i，j，m=2，3，5实际上，刚度方程可以写成如下的一般形式：k表示单元的节点组成首先我们运用规那么eq\o\ac(○,2)，即各个节点力的合力等于节点载荷所以也就有：我们将上面得到的方程按这个规律相加：对这个等式运用规那么eq\o\ac(○,1)，即并整理有：将上述的这些方程写为矩阵的形式，就得到了该问题的总刚方程为：我们对这个总刚方程进行分析，可以得出以下的规律总刚矩阵的规律总刚方程就是运用所有节点的位移向量和总刚矩阵相乘，得到结构的节点载荷，节点载荷向量〔外界对弹性体作用的载荷〕很容易求得，节点位移向量是未知量，关键就是总刚矩阵如何得出。如都按照上述的步骤，成千上万个节点，该如何作呢？所以有必要总结总刚度矩阵的规律。eq\o\ac(○,1)当r=s时，即子块矩阵[Krs]是总刚矩阵主对角线上的子块矩阵。由环绕节点r〔s〕各单元刚度矩阵的相应对角线子块相加。eq\o\ac(○,2)当r≠s时，但r、s是相邻单元的公共边上的节点时，子块矩阵[Krs]等于两相邻单元刚矩阵的相应子块矩阵相加。如上例中的[K13]、[K31]等。eq\o\ac(○,3)当r≠s时，且r、s只是一个单元的边上节点时，子块矩阵[Krs]就等于该单元刚矩阵的相应子块矩阵。如上例中的[K12]、[K21]等。eq\o\ac(○,2)当r≠s时，且r、s不属于一个单元的边上的节点时，子块矩阵[Krs]等于零。如上例中的[K15]、[K51]等。利用上述规律不仅可以检验总刚组集的正确与否，而且可以直接组集总刚矩阵〔手工组集〕。但这种方法不利于计算机编程。下面介绍计算机组集总刚矩阵的方法总刚矩阵组集的步骤eq\o\ac(○,1)扩大各单元刚度矩阵，使之成为与总刚矩阵相同的阶数。eq\o\ac(○,2)按照总体节点的编号，将各单元刚度矩阵的各个子块移到相应的位置上。其余位置充零。eq\o\ac(○,3)把各个改造过的矩阵直接相加，就得到总刚矩阵。我们以单元eq\o\ac(○,4)为例讲解步骤eq\o\ac(○,1)、eq\o\ac(○,2)单元4的i、j、m=2、5、3，所以单元刚度方程是单刚矩阵总刚度矩阵的性质总刚度矩阵由单元刚度矩阵组集而成，所以也具有单元刚度矩阵的一些性质，如相同的物理意义，位置无关性、对称性和奇异性等，还具有以下性质：稀疏性由规律4知，当rs，且r、s不属于同一单元的两个节点是[Krs]=0，说明互不相关的节点数愈多，零子块矩阵就愈多。一般说来相关节点数不超过9，而整个分析对象常常成百上千、上万或几十万。如果整体有100个节点，那么百分之九是非零子块，而百分之九十多都是零子块。所以在总刚度矩阵中非零的子块是很稀少的。带状性所谓带状性，就是指总刚矩阵中的非零子块，集中分布在主对角线的两侧，呈带状分布。半带宽值就是计算这带状宽度的数值。它是以排列元素最长的一行，从第一个非零元素起至主对角线元素止，所有元素的个数。其数值可以由节点的总体编号算出。B=〔相邻节点号的最大差值+1〕×〔单个节点的自由度个数〕对于平面问题说来，就是：B=〔相邻节点号的最大差值+1〕×2以下图中的两种编号方式，可以的分别计算其半带宽值：图a，B=〔2+1〕×2=6图b，B=〔6+1〕×2=14半带宽值影响计算机存取总刚矩阵所需的内存大小。愈小愈好。由于半带宽值是直接受节点总体编号的影响。所以我们在建立有限元模型时，应慎重考虑，采用优化的方法。〔现在通常在程序中都配有这样的优化程序〕。总刚矩阵的压缩存取技术由于总刚矩阵具有对称性，所以我们只需存入主对角线上半带或下半带的元素，就可以完成解方程的运算。此即为总刚矩阵的半带宽存储方法。假定取主对角线上半带元素存储，具体做法就是每行元素以主对角线上的元素开始，存储每行半带宽数值的元素个数。如此各元素的行号不变，改变的只是列号。新列号和原列号的关系式如下：新列号=原列号-行号+1前面图示的例子，如果采用图示表示总刚矩阵的存储，就是：如以下图所示的情况。可以看出原总刚矩阵存储需要24×24的2维数组，改为半带宽存储后，就成为了24×6的数组了。另有一维的压缩存储方法，比这种二维的半带宽存储方法压缩更多。在这儿我们就不做介绍了。有兴趣的同学可以看相关的参考书籍。讲到这儿，实际我们已经知道，每个元素在总刚矩阵中的位置，在节点编号完成后，就完全确定下来了。所以计算机对总刚矩阵的计算，是一部完成的。即“边对号移置、边改列号、边累加”。总体边界条件的处理前面介绍总刚矩阵的性质时，说明了总刚矩阵是奇异矩阵。即。就是说总刚矩阵不存在逆矩阵。要求得节点位移的位移解，还必须引入边界条件。1，边界条件的类型eq\o\ac(○,1)节点固定，即eq\o\ac(○,2)给定节点位移值。即，处理方法1置“0、1”即首先在总刚矩阵[K]中与位移分量相对应的行和列元素改为0，但主对角线上的元素改为1，然后在节点载荷向量的列阵中，与对应元素的位移用代替，其余元素减去位移分别乘[K]中相应元素。例：某结构的总刚方程为，、，按照上面的方法修改如下：如果展开上式，马上就有、第二行移项从中可以看出，这样处理后不仅可以直接得到、，而且它们产生的效果也计入到了其他方程中。如果是固定位移情况，那么载荷向量中的变化是什么呢？〔好好想想〕2乘大数法在总刚矩阵[K]中，把与位移位移相对应的行与列主对角线上的元素乘一个很大的数1010，然后把载荷向量中的对应元素代以给定位移乘以相应主对角线上的元素，再同样乘以一个很大的数。同上例。考察第一个方程：两边同除以，由于＞＞k1j，所以u1=c1同理可得v2=c4，〔思考该方法不能用于固定位移的情况，为什么？〕必须说明的是，当总刚矩阵采用半带宽存储的时候，以上边界条件的引入也应当按照半带宽存储的格式修改，否那么就是错误的。实际应用中，当是第一种情况时，采用方法eq\o\ac(○,1)；当第二种情况时，用方法eq\o\ac(○,2)。应力计算总刚方程在引入足够的边界条件后，就可以进行求解了。所谓足够的边界条件就是要使系统是一个几何不变的系统。数值解方程的方法大多采用高斯消元法。当求解出位移向量以后，就可以通过单元的几何方程求解应变，通过物理方程求解应力了。这个过程称之为回代过程。eq\o\ac(○,1)求解应变eq\o\ac(○,2)求解应力由于三角形单元是常应力与常应变单元，所以由上述方法求得的应力或应变看作是形心处的应力或应变。而且很明显，这样求得的应力或应变是不连续的。为了推算弹性体内某一点的接近实际应力，我们通常采用以下两种方法，来平滑应力的突变。eq\o\ac(○,1)绕节点平均法——就是将绕某一节点的各单元形心应力加以平均，来表示该节点的应力。如左图所示。该方法在各单元面积相差不大，单元的内部节点时较为准确。而在外部边界节点上，那么不理想。所以对于边界节点多采用三点插值的方法求得。也就是用内部的节点来推算。eq\o\ac(○,2)两单元平均法——把相邻两单元的应力相加平均，表示两单元公共边中点的应力值。边界边上的应力值也采用插值方法得到。至此我们对于有限元求解平面问题的方法全部讲解完毕了。下面通过一个例子来看看，具体如何计算的。§2.15例题计算一，有限元法求解平面问题的步骤：根据分析对象和给定的条件，按照一定的比例绘出计算简图，该图应有各部位的尺寸、外力和支承情况。选择适宜的坐标系，划分单元准备数据。数据包括节点数、编号和坐标，单元数、编号和节点组成；材料特性常数；载荷大小及位置；边界条件。计算单元的面积、应变矩阵、弹性矩阵和单元的刚度矩阵。组集总刚矩阵，处理载荷移置，形成载荷向量。处理边界条件，修改总刚矩阵，得到最后的刚度方程。求解线性方程组，得到位移向量。求解应力值。例题设有如下图的正方形薄板，在对角线上作用有沿厚度均匀分布的载荷，其合力为2kN，板厚为1个单元厚度，为使计算简便，假定材料的μ=0，求该板的变形。分析：该薄板对称于对角线，受力也是如此，所以可取其1/4来计算。由于是取1/4，所以该模型的两个直角边，都是对称线。X的对称线不能有x方向的位移，y轴的对称线不能有y方向的位移。直角的顶点既是x轴对称线上的点，又是y轴对称线上的点，所以该点在x、y方向的位移都应该约束。最后得出的计算模型如下。解：1准备的数据A单元节点数据坐标123456X001012Y211000B单元信息数据eq\o\ac(○,1)eq\o\ac(○,2)eq\o\ac(○,3)eq\o\ac(○,4)I1223J2455M3536所有单元的板厚是“1”C材料常数弹性模量E，泊送比μ=0D载荷数据E位移边界2计算单元刚度矩阵单元1面积形成的几何矩阵弹性矩阵单元1的刚度矩阵同理可求得、和。需要说明的是可以利用位置无关性直接得出和，由于它们是1单元平移得到，所以单刚和一样。3组集总刚4引入边界条件，采用划行划列〔由第一中方法演变而来〕就是将边界条件的节点位移所在的行和列全部划掉。得到的就是需要求解的位移向量，最后得到的总刚方程为：解以上这个方程，可以得到以下的解。按照每个单元的节点组成，从总节点向量位移中挑选出每个单元的节点位移向量，运用前面单元分析中得到的公式，回代就能计算应变、应力了。〔是哪个公式？〕§2.16四节点矩形单元前面我们运用三节点三角形单元，获得单元的应变和应力是不变的。这和实际情况有着明显的差异〔实际中应力应变在单元中是连续变化的〕。产生这种差异的原因，就是我们的位移模式过于简单。下面我们将就如何选择位移模式，提高有限元计算精度，进行说明。同时也可以就如何选择单元，有个初步的认识。首先讲解四节点的矩形单元。一，矩形单元位移模式如下图的矩形单元，四个节点分别为i、j、m、p。单元节点位移向量对应的会有8个节点力分量。由于矩形单元的节点位移向量是8个，那么根据前面我们学过的三角形推导过程知道，可以建立8个方程，求解8个待定系数。所以在位移模式中的待定系数可以是8个。即：由于我们选择的是多项式插值。所以对于2元函数多项式项元的选择，可以用帕斯卡三角形确定。如左图。从该三角形可以看出，2元多项式的二次项有三项，即x2、y2、xy。如果选取前面两项中的任何一项，都会造成位移模式偏惠与你所选择的那个方向。从而使位移出现不对称的情况。只有增加xy这一项，才能防止这种情况出现。也就满足了我们材料各向同性的要求。选取局部坐标系，将四个节点的坐标代入该公式，并求解待定系数，最后可以得出和前面三角形单元类似的关系：其中的在上面公式中，代表节点坐标的是什么？位移模式的讨论。1，应变和刚体位移由几何方程得同三角形单元位移模式一样，常应变和刚体位移包括在位移模式之中。同时我们可以看到，应变中还有随x、y的变量应变。所以在矩形单元不再是常应变或常应力单元了。2，在局部坐标系中，当x=±a〔y=±b〕时，位移模式是x〔或y〕的线性函数〔将x=±a〔y=±b〕代入位移表达式中即可〕。这说明在边界上位移按线性规律变化，在公共边界上的位移是连续的，满足收敛性的充分条件。二，矩形单元的刚度矩阵应力向量。分别代入平面应力和平面应变的弹性矩阵，就可以得出两类问题的应力解。刚度矩阵代入[B]和[D]矩阵，相乘后逐项积分即得单元刚度矩阵，见教材中的具体计算公式。需要说明的是，由于[B]不再是常数矩阵〔还有x、y〕，所以积分运算较三角形单元复杂些。具体见如下矩阵。矩形单元的载荷向节点移置的方法同三角形单元总刚矩阵组集的原那么和方法、边界条件的处理等也同三角形单元。其中不同的位置仅是单元节点位移和节点载荷的数目，不再是6个而是8个。三，讨论矩形单元的优缺点：1〕其位移模式推导出的应力、应变不再是常量，分布更接近实际物体中的分布。从而使这种单元的计算精度更高。2〕对于斜边性边界和曲线边界的拟合性差，且不便于在不同部位采用不同大小的单元。实际应用中常需要三角形单元在大小不同的矩形单元之间进行过渡。见以下图。§2.17六节点的三角形单元从以上对矩形单元的推导，可以看出，增加单元节点的数目可以提高单元位移模式的插值次数，进而可以提高单元计算精度。沿用这样的思路，我们对三角形单元也可以采用相同的方法。不过增加的节点是在每条边的中点。这样就得到了一个六节点的三角形单元。如下图。一，位移模式由于这样的单元有6个节点，所以节点位移向量的个数就是2×6=12个，那么在位移模式中可以有12个待定系数。参看帕斯卡三角形，我们选用以下的位移模式：将每个单元公共边上的方程代入，可以求得一个二次的抛物线方程，而公共边上有三个公共节点，所以可以唯一确定这条抛物线。这说明该插值模式满足位移连续的收敛性充分条件。仍然可以按照前面三节点三角形单元的方法，分别将6个节点位移代入，然后解联立方程组，求12个待定系数。这样的计算，由于待定系数过多，计算过程也过于繁杂。所以实际中是采用面积坐标的形式进行计算。二，面积坐标三角形三个顶点i、j、m，p为三角形中任意一点，其在三角形中的位置，可以用来确定。其中：A——三角形的面积。Ai——是三角形pjm的面积。Aj——是三角形pim的面积。Am——是三角形pij的面积。当P在图中虚线上任一点是，Li是相同的〔三角形面积为底×高的一半〕，在该三角形的三个顶点，分别有：i，Li=1、Lm=Lj=0j，Lj=1、Lm=Li=0m，Lm=1、Li=Lj=0而且有，将节点的面积坐标和前面我们推导出的形函数Ni进行比拟，可以知道，形函数实际就是面积坐标。三，六节点三角形单元的位移模式推求如下：Ni在i节点处应为1，在j、3节点处应为零，面积坐标Li虽在i节点处为1，但在2、3节点处却为1/2，所以还是用Li作为Ni就不行了。为此我们构造一个函数，考察该函数是否满足形函数的性质。由于Li在2、3节点处，所以很显然上式等于零。在i节点处要为1，那么β=2。〔令=1，可以计算出来〕所以同理可得对于N1应满足j、m节点处为零，所以构造函数，令=1〔在节点1处应为1〕，可以求得β=4所以，同理得综合以上推导，六节点三角形单元的形函数如下：1、2、3、i、j、m循环四，单元刚度矩阵推导过程类似于三角形单元的刚度矩阵，不同之处在于在对形函数求偏导时，要运用复合函数求导的方法〔因为插值基N是面积坐标的函数，面积坐标才是x、y的函数〕如：在求得[B][D]矩阵后，运用就可以求得单刚矩阵。但是由于[B]矩阵的复杂性，所以求解[K]e也不再是很容易的一项工作了。必须运用一些相关的积分公式。具体的计算公式可以见教材的相关内容。五，讨论在单元数目相同的情况下，六节点的三角形单元计算精度远比三节点三角形单元高，也比矩形单元高。但由于一个节点的相关节点数目，在6节点的三角形单元中大大增加，所以总刚度矩阵的的带宽也较三节点三角形单元宽的多。所需的内存也相应增加。从理论上说来，运用这种增加节点，改善单元计算精度的方法，可以不断的运用下去。如4节点矩形单元可以变为8节点矩形单元等等，但是实际运用中，我们可以看到，节点数目的增加，导致单元计算复杂性的大大增加，有时可能求解不出单刚计算公式〔在现有的计算分析理论上〕。所以在有限元程序中，再高次数的单元运用就比拟少了。6节点三角形单元为二次单元，而三节点三角形单元为一次单元。§2.18平面问题的计算实例例1一直齿齿轮，齿数为20，模数是3，齿厚20mm，压力角为20°，试分析其齿的受力状况。例2变厚度圆筒的压力容器，受有内压，试用有限元法计算圆筒内外壁的应力。第二章内容小结一，根底理论1，外力、应力、应变和位移的概念2，两类平面问题的区分3，弹性力学的3大方程4，位移相容方程的理解5，虚功方程二，有限元概念1，单元划分的原那么2，位移模式、形函数及其性质、保证有限元解收敛的条件3，单元刚度矩阵的推导过程4，非节点载荷的等效移置5，总刚矩阵的组集原那么、规律和方法6，总刚矩阵的半带宽压缩存取方法7，边界条件的类型及处理方法8，应力求解及平滑处理的措施三，三种单元的不同之处1，三节点三角形单元2，四节点矩形单元3，六节点三角形单元第三章轴对称问题、空间问题和等参数元问题工程中有许多零件是对称于一条轴线的回转体，如飞轮、螺杆和发动机的汽缸套等。当它们受到的载荷和约束也是对称于他的轴线时，弹性力学就将其归属于轴对称问题求解。与平面问题类似，也是将空间问题，简化为二维问题处理。本章的主要内容轴对称问题的根本方程轴对称问题的有限元方法空间问题的根本方程空间问题的有限元方法等参数元的概念§3.1轴对称问题的根本方程根据轴对称问题的特点，一般分析时采用圆柱坐标系，而不是直角坐标系，所以首先我们应该清楚两者之间的转化关系：直角坐标与柱面坐标的关系：自变量子午面：就是通过roz的任一平面，它是一对称面。轴对称的弹性体变形后，其子午面上任意一点P只在该平面上发生位移，而与θ无关。即P点的应力、应变和位移都只是r、z的函数，所以轴对称问题也是二维问题。弹性体变形后，不发生歪曲，所以任一点只有径向位移u和轴向位移w。如左图所示的微元体，其应力分量相应的应变分量其中的与u、w类似于平面问题中的关系，即现在我们要推求的就是轴对称问题中不同的，看左图可以推出：综合起来，就是轴对称问题的几何方程：其物理方程参照直角坐标系中的可以写出：转化为应变表示应力的矩阵形式：虚功方程的柱面坐标表达式：〔因为σ中不含θ的变量，所以〕§3.2轴对称问题的单元分析轴对称问题的单元是圆环形单元，其断面常常采用三角形或矩形，单元之间的连接不是点，而是园线，称作节园。虽然其单元与平面问题中的单元形式不同，但因其可以取任一子午面进行分析，所以在子午面上的单元界面却构成类似于平面问题的三角形网格，因此可以采用平面问题的分析方法，所不同之处是：1单元是圆环体2单元之间由节园铰接3节点力与节点载荷是施加在节园上的均布力4单元的边界是回转面下面推导单元的刚度矩阵一，单元位移模式节点位移向量节点力向量选线性插值函数作为位移模式如果用插值基函数表示，就是：其中：i、j、m循环二，单元刚度矩阵代位移函数到几何方程中其它局部类似于平面问题，只有其中的轴向应变是其特色。所以公式中的i、j、m循环要引起注意的是，fi是r、z的函数，说明其对应的应变分量εθ是随r、z坐标变化，不是如平面问题中的那样是常量。应力×其中应力矩阵的分块矩阵为i、j、m循环上式和一起代入虚功方程令由于[B]中的fi、fj、fm是r、z的函数，所以上式应该积分运算。但为了消去f在r=0处的奇异性，同时简化积分运算，我们取单元形心处的代替公式中的r、z，如此以来[B]中就没有r、z的变量了。这种近似处理，在网格划分较密集时，完全能够满足工程计算的精度要求。是单元形心处的周长，相当于平面问题中的厚度t。仍然采用分块矩阵的方法记忆单刚：r、s=i、j、m§3.3载荷的移置作用在环形单元上的体积力、面力和集中力都应分别移置到单元的节园上，形成节园载荷。移置后的单元载荷表示为：移置的原那么仍然是虚功相等的原那么。在这里我们应注意与平面问题的不同之处。集中力的节园载荷平面问题的是体力的节园载荷面力的节园载荷一，单元自重的移置如果自重沿z轴向下作用，那么可以推导出节点载荷计算公式：具体推导过程可见教材。二，单元离心力的移置离心力三，面力的移置如果作用有如下图的均布载荷，当是均布载荷时那么可以推导出：§3.4计算实例例1运用有限元法计算螺栓螺母之间受力的分布例2轮毂热套在轴上紧配合，工作时由于轮缘温升，导致轮毂与轴之间的压力分布变化，要求计算这种变化。§3.5空间问题一，弹性力学空间问题的三大方程1，平衡微分方程2，几何方程3，物理方程二，空间问题的有限元法一般采用四节点的四面体单元，如左图所示其节点位移向量其节点力向量单元共12个自由度，取12个待定系数位移模式为：代入节点坐标，求得的插值基函数：i、j、m、p循环其中：V——四面体体积，二，单刚矩阵推导过程完全同平面问题的过程，只是方程的维数增加了。得到的单刚矩阵，其分块矩阵为其中的§3.6等参数单元一，四边形等参数单元1，平面单元的位移模式及坐标转换如左图所示的任意四边形单元，仍采用同矩形单元相同的位移插值模式，即：现在考察其边界的位移情况，由于任意四边形的边界不再象矩形单元的边界，可以表示为x=±a〔y=±b〕的形式，而是的直线方程，代入其上位移模式方程：同理得以上方程说明在任意四边形的边界，其位移不再是线性方程了，而是一个抛物线方程，所以相邻单元之间的位移连续条件不能满足。〔两点不能唯一决定一条抛物线〕，所以这种位移模式不满足位移收敛的充分条件。从数学的观点来看待矩形和任意四边形的关系，可以发现任意四边形是矩形的畸变“映像”。就好比园与椭圆之间的关系。如果能够找到一种映射关系，确保四边形与矩形之间是一一对应的关系，那么就可以在局部坐标系中采用矩形单元分析〔这样可以保证位移连续的条件〕，然后用坐标转换的映射关系，转换到总体坐标系中成为任意的四边形单元，由于这种映射是一一对应的，所以也就保证了任意四边形单元边界的连续性，这样就回避了采用更高次的位移插值模式的问题。这就是等参数元的根本思想。矩形单元只是分析问题的过度单元，称作根本单元；任意四边形单元是实际单元，称作实际单元。为了描述问题的清楚方便，命局部坐标系为，那么根据前面的推导，得：其中：设矩形单元与任意四边形单元的转换关系如：将四个节点的局部坐标值和总体坐标值代入，就可以确定8个待定系数，如果我们也写成型函数的形式，就有：*其中的插值基函数同上。这种位移以节点位移值插值，位置以节点坐标值插值，而且插值参数相同，插值基函数也相同，具有这样性质的单元就称为“等参数单元”，简称“等参数元”。以上的转换关系，在数学上是可以证明其正确性的。我们在这儿只对其进行一下说明。当时，是矩形单元的jm边，将该值代入式*，**该式是以节点j、m坐标值为参数的直线方程，显然它也是实际单元的jm边。以代入式**，得、，是实际单元的m节点。以代入式**，得、，是实际单元的j节点。二等参数单元的单刚矩阵是ξ、η的函数，而ξ、η与整体坐标x、y之间又有坐标转换关系式，所以它是复合函数求导，即：代入形函数Ni令该矩阵称为坐标转换的Jacobian矩阵其中以上的推导过程说明了，原本形函数对整体坐标x、y求偏导的运算，转变为求局部坐标ξ、η的偏导。原本对实际单元的运算转为对根本单元的运算，附加求雅可比矩阵的运算。逆矩阵中的四个元素：由于[B]矩阵都是ξ、η的函数，所以单元刚度矩阵也应转换到局部坐标系下。首先看左图。沿局部坐标的ξ、η作微向量、，由于在ξ方向只是ξ变化，而η不变，所以同理得微向量在两坐标轴上的投影由它们组成的四边形面积：而，所以如果写为分块矩阵的形式：以上表达式是一个十分沉繁的积分运算，一般都只能用高斯数值积分运算。四边形单元的边界逼近性较矩形单元大大提高了，但仍然防止不了折线代替曲线的情况，如果要构造曲线边界的单元就必须增加单元边界上的节点数，使之称为由三个节点插值才能求得的情况。如此就可以构造以下的几种等参数元。其推导过程可以参看教材或其它参考书二高斯数值积分等参元中，形如以下的积分运算：由于f的复杂性，且常常不能列出显式，因此在有限元中常采用数值积分方法，即在积分区间上选取某些点，计算函数f在该处的函数值，然后用加权系数，乘以这些函数值，再求和，用该和值近似代替积分值的方法。高斯积分法的积分点值和该点的权重值，如下表所示：积分点数N积分点的坐标值xi权系数ωi234那么一维问题二维问题三维问题例：第四章梁杆结构问题所谓梁杆结构是指其长度比横截面尺寸大很多的梁和杆件、以及由它们组成的系统，这一类结构的应力、应变和位移都是一个坐标的函数，所以属于一维单元问题。§4.1梁干结构的单元划分一，分类1，平面桁架特点：杆件位于一个平面内，杆件间用铰节点连接，作用力也在该平面内。单元特性：只承受拉力或压力。单元划分：常采用自然单元划分。即以两个铰接点之间的杆件作为一个单元。2，平面刚架特点：与平面桁架的不同在于杆件之间以刚性铰接点作为连接。单元特性：除承受拉、压轴向载荷外，还要承受剪力和弯矩。单元划分：除采用自然单元划分外。即以两个刚接点之间的杆件作为一个单元。同时也要考虑在集中载荷作用位置、截面突变位置和分布载荷突变的位置增设节