有限元分析教案

第一章有限元法概述在机械设计中,人们常常运用材料力学、结构力学等理论知识分析机械零构件的强度、刚度和稳定性问题.但对一些复杂的零构件,这种分析常常就必须对其受力状态和边界条件进行简化.否那么力学分析将无法进行.但这种简化的处理常常导致计算结果与实际相差甚远,有时甚至失去了分析的意义.所以过去设计经验和类比占有较大比重.由于这个原因,人们也常常在设计中选择较大的平安系数.如此也就造成所设计的机械结构整体尺寸和重量偏大,而局部薄弱环节强度和刚度又缺乏的设计缺陷.近年来,数值计算机在工程分析上的成功运用,产生了一门全新、高效的工程计算分析学科一一有限元分析方法.该方法彻底改变了传统工程分析中的做法.使计算精度和计算领域大大改善.§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,虚功原理法〔变分法〕.3,加权余数法.直接法概念浅显,易于理解物理含义.变分法需要泛函的数学知识,其推导过程具有严谨的数学概念.加权余数法适用于泛函不存在的应用范围.本教材将运用虚功原理方程结合弹性力学和材料力学中的知识来推求几种常见单元的单刚计算公式.现在先看一个简单的阶梯轴的轴向拉伸问题例：如下图的变截面直杆,受拉力P,运用有限元方法分析其变形.取任意单元,长度为1,面积为A,单元内任意一点的轴向位移和其位置坐标成正比,即

u=a0+aix其中a.、a1为待定系数.由于杆的两个端点节点1、2是单元上的点,所以它们应该满足上述方程.节点1,x1=0,u1=ao+a1X0=a.节点2,x2=i,u2=ao+a1xla1=〔u2-u1〕/l将求出的结果带入方程并整理,就得：U=Uix=U=Uix=1--u1+-u2

<l)l-NiN2PUiU2=NU式中：m、也是形函数[N]形函数矩阵{8}e节点位移向量由位移与应变的关系知道：udu-ududxdxdx将上面推出的位移表达式代入,可得:udu-udu-u_du_£N；?一dxdxdx=B?：、/e—x——1--dxILl上式中的[B]称为应变矩阵或几何矩阵.运用材料力学中的虎克定律,可以将应变和应力联系起来.单向应力状态的虎克定律为EVEVu1i」J2二SL.：,e[ke=B[ke=BTDIBdvv[D]矩阵是弹性矩阵.对于一维单元来说,就是所以我们这儿讨论的例题：-11E.1mEA11Adx「111其中[S]称为应力矩阵.利用虚功方程可以建立力与位移之间的关系,也就是单刚方程.在后面我们将会推导出它的一般形式如下：式中：{F}e为单元节点力向量,对我们这个例子应为[UiU2]T.[K]e为单元刚度矩阵.后面将推导出它的计算公式为求得单刚矩阵,也就完成了单元分析.总结单刚矩阵推导的步骤,应该分为四步：1〕假定单元内位移变化的近似规律,即选择位移模式.2〕运用几何关系,推求位移与应变的关系.

3〕应用物理规律,把应变与应力联系起来.4〕运用虚功方程的力与位移关系,求出单刚矩阵.单元分析是整个有限元分析的核心.不同的单元由于其力学特性不同,而具有不同的单元刚度矩阵,我们这本教材就是要学习几种常用单元矩阵的推导和计算.了解各种单元的力学特性,为以后选择单元类型打好根底.三,整体分析1,由各单元刚度矩阵组集成整个结构的总刚度矩阵.ki=aki=a：1tli1-11一1-1il2[-11一整个结构有三个节点,首先将单元刚度矩阵扩充为3X3的矩阵,移动各元素使之与单刚矩阵中的元素位置相对应,如下:11-10K1=EA1-110I1000_然后直接相加.000k1=受01—1I2-0-11_EAEAI1—生EAEAI1—生I1I1EA.EA2I1l2EA2120EA2I2ea2I22,把各单元的节点力向量组集成总的节点载荷向量.riR=<0,03,根据边界条件,修改总刚度矩阵,获得总刚方程组.边界条件修改之前的总刚方程：EA1EA2I1EA1EA2I1I2EA2I20EA2I2EA2U1"2?

lu3,I2修改以后〔采用置“0,1〞法〕EAiiEAiiiEAi110EA尸EAi於7T0|Ui0<U211"J4,求解方程组,得出总的节点位移向量.解得的解是：Ui—Pll*Pl2eA1eAT.fUi—ea2J3,0有了节点位移,再回代到前面单元推导过程中的公式X和XX,就可以求得每个单元的应变和应力了.从这个简单的例子,我们了解了有限元法求解力学问题三大步骤中的内容,想必很多同学会说,这样复杂,如果运用材料力学的知识,我还来得快些.但是大家不要忘记,有限元的计算很多都是编程完成,而且现在很多的商业软件都已经完成了很多的工作.我们学习有限元主要是了解它的原理,并对常见单元的力学特性有所了解,这样对于以后运用有限元起到帮助作用.所以下面章节的内容,就是围绕这个主题展开.要到达这个目的,我们还必须学习必要的弹性力学知识.对弹性力学知识的学习,也对我们以后把握问题的本质有帮助.第二章平面问题平面问题在力学研究的课题中属弹性力学的范畴.该类问题不仅本身具有典型性,而且在机械零构件的分析中,也是应用得非常广泛.所以这类问题也称之为经典的力学问题.我们知道,实际的机械零构件都是具有三维空间尺寸的物体,理应作为三维对象处理,但是当物体的几何形状和受力状态处于某些特定的情况下,近似地简化为平面问题,不仅可以大大简化计算的工作量,而且其精度也完全能够满足所要求.如：直齿圆柱齿轮可在垂直与孔轴线的截平面内作平面应力分析就足以了解整个齿轮的受力状态；大坝的横断面可作平面应变分析来了解整个大坝受力情况等.本章是全书的重点,在这里不仅介绍弹性力学的根本知识.还将系统地讲解有限元的基本概念、原理和方法.是学习以后各章节的根底.§2.1外力、应力、应变和位移外力、应力、应变和位移的概念在材料力学中已经学习过,由于这些概念在弹性力学、有限元法中具有和在材料力学中不同的规定,弹性力学中的规定和有限元法是完全相同的,所以在这里我们将根据弹性力学的习惯表达方法把他们集中的加以阐述.一,外力外界作用在物体上的作用力,可以分为两大类：1〕体积力分布在物体体积内的力.如：重力、惯性力和磁性力等.单位体积的体力在坐标轴上的分量X、Y、Z,称为体力分量.符号规定为沿坐标轴正向的为正,沿负向的为负.2〕面力作用在物体外表的力.如轮压、水压等.面力在坐标轴上的投影,表示为X、Y面力在坐标轴上的投影,表示为X、Y、Z.符号沿正弹性体受到外力作用后,内部产生的反抗变形的内轴为正,负轴为负.力.二,应力力.以弹性体中P点为定点的微单元体来考察.所谓微单元体,就是图中PA、PB、PC的边长分别为dx、dy和dzo以下简称这样的微单元体为微元体.微元体每个面上的应力都可以分解为三个应力分量.以图中红面为例,分别是八、.xy、Txz0应力命名的规那么：以应力所在面垂直的坐标轴为第一个下标,应力指向为第二下标.如果下标相同就用一个下标表示.符号规定：正面上的应力与坐标轴同向为正,反之为负.负面上的应力与坐标轴反向为正.反之为负.所谓正面就是面的外法向与坐标轴同向为正.反之为负面.作用在两个相互垂直的面上,并且垂直与该两面交线的剪应力互等.即：6个,它们是:Txy=Pyx；Pyz=Tzy；P6个,它们是:如此以来,代表P点应力状态的应力分量应有J'-'二y二zxyxz-yzT三,位移任一点的位置移动用u、V、w表示它在坐标轴上的三个投影分量.符号规定：沿坐标轴正向为正,反之为负.四,应变弹性体内各点的位移在受力后一般是不相同的.各点之间距离的改变,从而使物体形状发生变化,即所谓的变形.而物体的形状总可以用它各局部的长度和角度表示.长度的改变称为正应变£,角度的改变称为剪应变丫.以微元体三个棱边的线伸长和角度的变化,就分别有和6个应力分量相对应的6个应变分量,即：L」；x；y；zxyxzyzT为与前面符号规定一致,这里对剪应变的符号规定如下：正应变伸长为正,缩短为负；剪应变使直角变小为正,变大为负.§2.2两类平面问题前面我们讲过,实际受力物体都是三维的空间物体,作用在其上的外力,通常也是一个空间力系,其应力分量、应变分量和位移也都是x、v、z三个变量的函数.但是当所考察的物体具有某种特殊的形状和特殊的受力状态时,就可以简化为平面问题处理.弹性力学中的平面问题有两类.,平面应力问题当物体的长度与宽度尺寸,远大于其厚度〔高度〕尺寸,并且仅受有沿厚度方向均匀分布的、在长度和宽度平面内的力作用时,该物体就可以简化为弹性力学中的平面应力问题.我们分析以下其应力特征.当2=±1/2时,有〔Tz=0、Pzx=0、Tzy=0o由于板较薄〔相对于长度和宽度尺寸〕,外力沿板厚又是均匀分布的,根据应力应连续的假定〔弹性力学中的根本假定〕,所以可以认为,整个板的各点均有bz=0、°zx=0、Tzy=0o如此以来,描述空间问题的6个应力分量也就变为了3个,即七:'-Lx二yxyT而且这些应力分量仅是x、y两个变量的函数.,平面应变问题和平面应力相比拟,平面应变是当物体是一个很长、很长的柱形体,其横截面沿长度方向保持不变,物体承受平行于横截面且沿长度方向均匀分布的力时,该问题就可以简和平面应力相比拟,平面应变是分析其应力特征.假定其长度方向为无限长,

那么任一横截面都可以看作是物体的对称面,如此那么有该面上的点都有w=0,也就是横截面上的所有点都不会发生Z方向的位移.由这一点可以推出也就有sz=0、Tzx=0、Tzy=0o£z=0,那么是否也就有bz=0呢可能有同学想d=E£,当然也就有bz=0,这是错误的.平面应变状态下bzW0的.虽然不等于零,但它也不是一个独立的变量了,它由

by的大小而决定.如此以来,独立的应力分量同平面应力问题一样也是3个:三,两类问题的比拟1,几何特征平面应力厚度〈V长度、宽度平面应变厚度＞＞长度、宽度为便于说明可讲上述长度看作为厚度2,受力特征外力都必须在其面内且不沿厚度方向变化3,应力特征平面应力6r0、Tzx=0、Tzy=0.&产0自由变形〔无约束〕平面应变,3zW0但不是自变量、Tzx=0、7rzy=0o£z=0总以上比拟可以看出,平面应力是真正的2维〔平面〕应力状态,而平面应变却不是,而是3维应力状态,只不过bz不是独立变量而是随横截面平面应力分量而定.独立变化的应力分量只有3个,类似于平面应力状态.衡微分方程弹性力学求解问题是从静力学、几何学和物理学三方面综合考虑的.所以我们首先微元体应该满足的平衡条件一一平衡微分方程.我们以平面问题为例推导,看看它应该具有什么形式.首先对平面问题的微元体进行受力分析图,如左所示.物体静力平衡的条件是：汇Fx〔y〕=0;汇M=0.先看EFx=0「x,±F1dxdy1-「xdy1-%,二Tyxdydx1-yxdx1,Xdydx1=0展开化简得同理可求得汇Fy=0满足得条件,色〕+色〕我=0N＞x由汇M=0,列出方程如下：卜熊町1喈十山咚-1管尸1岩-…岩=.化简后得:xy2fxdxyxy2fxdxyx2：:ydy略去微量项,可得：“y=3这就是前面所将的剪应力互等.对于平面应变问题,微元体的前后面还有正应力bz,不过它们是互等的.对于推导出来的结果,没有任何影响.所以平面问题的平衡微分方程就是：三.工.丫二0X=0写成矩阵形式何方程考察平衡微分方程,其中具有三个未知变量bx、by、Pxy,,而只有两个方程,方程具有无数个解.说明仅从静力学关系无法求解该方程.我们必须从其它方面寻求帮助.弹性体在受到外力后,会发生位移和形变,从几何上描述弹性体各点位移于应变之间的关系,就是弹性力学中的又一个重要方程一一几何方程.1坐标轴x的变化率,1坐标轴x的变化率,偏导数—：:y移为u、..,v,那么B点的位移就是:二uuR'=u——dxBex..'同理的D点的位移分量vB'二v=v—dx-:xuD'=uFu—dy

:y：:vvd''v—dy二y由于a角在位移和形变很微小的情况下非常小,所以AB'^AB线段AB位移后的总伸长量为仍然取截面的微元体ABCD,AB、CD边长为dx、dy,厚度为“1〞.位移u、v都是x、y的函数,即u(x,y)、v(x,y),偏导数出、0表示位移分量u、v沿x:x£v表示位移分量u、v沿坐标轴y的变化率,设A点的位-y

一,一一一〞一',cu,cu,AB-AB=AB-AB=ub-ua=u+—dx-u=—dx:x=—dy/dy=—：:y；:y；x=-udx/dx=-u=—dy/dy=—：:y；:y剪应变由a、3两个角度组成：：t：：tgABV-xdx-vdxu—dx-u

一：x由于里Y1,所以口=生,同理可得日=包一x二x二y:v：:u…xv=■=.xy--二x二y0Li.Cy综合以上几何方程,并将它们写成0Li.Cyyyxy由以上方程可以看出,当弹性体的位移分量确定以后,由几何方程可以完全确定应变,反过来,应变却不能完全确定弹性体的位移.这是由于物体产生位移的原因有两点：1〕变形产生的位移.2〕因运动产生的位移.因此弹性体有位移不一定有应变,有应变就一定有位移.理方程描述弹性体内应力与应变关系的方程,我们称之为物理方程,也叫材料的本构方程.弹性力学通常研究的是各向同性材料,在三维应力状态下的应力应变关系.当弹性体处于小变形条件下,正应力只会引起微元体各棱边的伸长或缩短,而不会影响棱边之间角度的变化,剪应力只会引起角度的变化而不会引起各棱边的伸长或缩短.因此运用力的叠加原理、单向虎克定律和材料的横向效应〔泊松效应〕,我们就可以很容易的推导出材料在三向应力状态下的虎克定律,也就是通常所说的广义虎克定律.式中,式中,EG材料线弹性模量-材料剪切弹性模量-材料横向收缩系数,即泊松系数.三者不是独立的.具有以下关系:

这些参数都是材料的固有属性系数,可以通过查材料手册获得.如：钢材的弹性模量=196〜206GPa之间,通常取2.1x105MPa,科=0.24〜0.28之间,也取为0.3进行计算,=79GPa.将以上空间问题的物理方程运用到平面问题,其形式如下：Zz=0、Tzx=0、Tzy=0.所以:1〕平面应力Zz=0、Tzx=0、Tzy=0.所以:1二-二x一二01=之x1二-七；x.」二yE121匕；zyTOC\o"1-5"\h\z=—T=TzycxyxyGE从以上物理方程,也论证了我们前面说的£zw0的结论,但由于它是由x和y方向应力产生的附加无约束变形,所以通常不予以考虑.在有限元分析中更多的是运用应变表示的应力关系,所以我们将上式变形一下：xyxy二E21」xy以上方程的矩阵表达形式为:ax,ayax,ayyxyE1-」201-口2yxy简记为：七〕-D,'.；〕式中:｛b｝、｛£｝为该问题的应力、应变向式中:[D]为弹性矩阵.它是一个对称矩阵,且只与材料的弹性常数有关.2〕平面应变问题的物理方程由于&z=0所以由空间物理方程的第三式得:〕,代入〔1〕、〔2〕式得一2CT1-」y11-」2P二一x1-」zy1TxyG21日一■xyE-xy同理变型为应变表示应力的形式：十—十——1-J-yE1-」22xyxy+z1」x1-UE1-j22E1-」2xy21+——1-U矩阵形式:□xVxy；x无\y转转换到平面应几何方程和物的.也可简记为七D「二比拟可以发现只需将平面应平面应变问题的弹性矩阵不同于平面应力问题的弹性矩阵,力问题弹性矩阵［D］中的材料常数E换为E/〔1-^2〕,科换为科/〔1-科〕比拟可以发现只需将平面应的弹性矩阵.其实弹性矩阵的这种转换方法,是弹性力学中将平面应力结果,变问题结论的一般方法.由于在两种平面问题的描述方程中〔平衡微分方程、理方程〕,只有物理方程是不同界条件求解弹性力学问题实际就c是在确定边界面问题而言〕,学的角度看,就融T条件下,求解8个根本c是在确定边界面问题而言〕,学的角度看,就融T是求解偏微分方程的边值问题.边界条件的给出通常是各式各样的.大体可以分为三类:1,第一类边值问题给定物体的体力和面力条件,确定弹性体的应力场和位移场.此类问题边界以力的形式给出,所以也称为应力边界条件.我们可以来考察一下应力边界的一般形式：ajivj=Ti五是在S,面上给出的力的分量.平面问题如左图所示,设阴影局部的微元体弧长为ds,厚度为单元厚度“1〞,其法线与X轴的夹角为由阴影局部微元体的平衡条件可以推出：Xds父1-crx•dscos6x1—e,dssine父1=0"人xyYdsx1-try,dssin6x1-zxy,dscos日x1=0'化简后得：◎xcos9+ixysin日=XJcy、_>此即为平面问题应力边界方程.仃ysine+%ycos@=Y2,第二类边值问题给出弹性体的体力和物体外表各点得位移条件,确定弹性体得应力场和位移场.由于以位移给出得边界条件,所以也称为位移边界问题.一般得位移边界条件为：5=也在Su面上.3,第三类边值问题给定弹性体得体力和一定边界上得面力,其余边界上的位移,确定其应力场和位移场.由于边界以力和位移两种形式给出,所以也称为混合边界问题.针对不同的边界条件,弹性力学求解的方法也有所不同.弹性力学的解题方法〔解析法〕1,应力法由于第一类问题的边界条件以应力形式给出,所以以应力作为根本的未知量的求解过程,就是人们通常所说的应力法.由于平衡方程中有三个未知量,而只有两个平衡微分方程,必须找出另外一个包含应力分量的方程,才能求得方程解.考虑到弹性体变形前是一个连续体,变形后也应是连续体的根本假设,所以要求微元体的变形一定要协调,才能使变形前、后,不会发生裂缝、重叠等现象.要使变形协调,就要研究几何方程.前面介绍的平面问题几何方程如下：分别对ex、ey、求y、x的二阶偏导,然后相加:fV二2fV二2xy十—I=.x:x:y上式说明三个应变分量之间应满足的连续性条件,我们称之为形变协调方程〔相容方程〕通过物理方程,使上述的形变协调方程换成应力表示的形式,使之与平衡微分方程就构成了应力法中需要求解的方程组.具体我们来看看：1〕利用物理方程消去相容方程中的形变分量〔以平面应力为例〕1f2.1.E守〞…y妥、…x2(1+N)3、xyEj1f2.1.E守〞…y妥、…x2(1+N)3、xyEjxjy_2_2「*-『『「yx=21・口：：2xy.:x；:y2〕利用平衡微分方程,消去上述公式中的剪应力xy-:y二二x-X_:xxxy:x/y—y-Y-:y二22"y；xyX求偏导,C2式对y求偏导,然后两者相加9.2_产二x:X二二y三Y———21f~~—x:x二y二y代入相容方程,化简二2——(er-kkr.2xy-y-2-

exd2OxdX32OyY.2一.2.x二x二y■:y-2_二:-x-2x.2C+—y

2~■：2c：__y-2x(仃y-N仃X)=I1+N)-1二十——ySo*ySo*y・二y11•日Y+——对于平对于平面应变而言,就可以了.运用前面讲过的物理方程的转换方法,只需将上式中的代以(1/(1-)1+CT1+CT)=y1-」.XFY-f-:x：y3〕最终求解的方程组平面应力问题三-三X二0一一Y二0,：y二x二x.二y=-1.」平平面应变问题三个微分方程,三个未知变量,再考虑边界条件,即可求得.如果是单连体〔只具有唯一的封闭边界〕的对象,满足了以上方程组后就是实际的解.但对于多连体〔具有多个封闭边界〕的对象,中还包含有待定系数,这些待定系数会导致位移的解出现多值性.所以对于多连体的问题,还应考虑位移的单值条件,才能最终确定.该局部的内容可以参见徐芝纶编的?简明弹性力学教程?中圆环受均布压应力的情况〔P87〕2,位移法位移法主要针对第二类边界条件问题求解.解题步骤：1〕改写物理方程使之成为应变表示应力的形式2〕应用几何方程表示以上得到公式中的应变3〕将它们代入平衡微分方程经整理最后得到的位移法求解平面应力问题方程为：E1-口E1-口222.uu十c2改1-J-2u1-2v+2cy22dxdy)E(d2v上1—□d2vJ+□d2u"1-N2(W22加2fixcy；两个未知量,两个方程,再加以边界条件即可求得问题的解.以上介绍的解析法中,应力法和位移法是求解弹性力学问题的根本方法.但都需要解联立的偏微分方程组.求解过程中的数学难度,常常导致这种求解是无法进行的.由于应力法在体力为常量的情况下可以进一步简化为求解一个单独的微分方程的问题,所以应力法在解析法中相对应用较多.但即使这样,在应力法中,也常常采用逆解法或半逆解法.体力情况下应力法的简化、应力函数及实例分析我们前面讲述了弹性力学的三大方程,及应用这三大方程的应力法和位移法解题步骤.但是也说了要求解这些联立的偏微分方程在数学上是存在很大难度的.很多情况下,根本无法进行.那么弹性力学如何在实际中进行应用,它们和我们前面学过的材料力学区别究竟在哪里我们将通过这一节的学习,一方面了解如何应用这些弹性力学的方程求解问题,另一方面加深对力学概念的理解,建立力学分析问题的直观感觉,为建立有限元模型打好根底.我们知道在大多数情况下,我们分析的对象,体力是常数,它不随x、y坐标变化.如此以来,前面讲解的第三个方程〔应力表示的相容方程〕,就可以简化为了：二2二2J+£y[&X+Oy〕=0简记为：▽2〔仃、+仃〕=0-2-2yxy：xn以上方程称为拉普拉斯微分方程,数学上也称之为调和方程,满足调和方程的函数称之为调和函数,及这里的〔仃*+byV2是拉普拉斯算子.

这样以来常体力情况下的应力法方程就是:二.吐X=0*：y三.三.丫一.:yY、'一二x•二y=0以上方程都不含有材料常数E、W,所以平面应力和平面应变两类问题具有相同的方程,这说明：在单连体问题中,只要边界相同、受同样的分布外力,应力分布与材料无关；也与是平面应力还是平面应变的状态无关.以上结论的意义：⑷弹性力学平面解答的应用范围加宽.②为实验应力分析提供了理论依据〔光弹实验〕下面我们考察平衡方程：其解由隹〕+匡=0其次微分方程的通解,加上任意一组特解组成.三■士工0特解特解我们可以很容易找到.如：:-x=-Xx二y=-Yyxy=0所以所以现在关键是找其次方程的通解.由第一个方程,可得：尚、2G已/1T"万=ly…分的充要条件.所谓全微分就是有一个函数由数学微分理论,该式是一个函数全微dA=：xdy+JdA=：xdy+J〞xydx且:二：:A-7一■xy：:Aex同理由第二式可得：;〕::Bexxy:By由剪由剪应力公式又知存在一个函数可以使d=BdyAdxB=—故::x由于应力与函数.存在这样的关系,因此函数.即是应力函数.我们用应力函数来表示相容方程：'J';=\4=0上式说明.是重调和函数.前面讲过在弹性力学中,常用逆解法和半逆解法.所谓逆解法就是设定各种满足相容方程的应力函数,运用bx、by与.的关系,求得应力分量,再考察其满足何种边界条件,从而知晓这样的应力函数可以解决什么问题.所谓半逆解法就是根据弹性体的边界形状和受力关系,设定局部应力分量为何种形式的函数,从而确定应力函数.,在运用应力函数求出所有的应力分量,根据边界条件确定应力分量应具有的最终形式.下面我们来看一个半逆解法的例子.运用逆解法求简支梁受均布载荷的应力分布.II由材料力学知,弯曲应力主要由弯矩M引起,剪应力由剪力引起,而挤压应力由分布载荷q引起.现在q为不随x变化的常量.因此我们设6y不随x坐标变化,即仃y=f(y),c2(P,、因此%=fyex我们对x积分:J=xf(y)+f1(y)；:x2=Yfyxfiyf2y上式中fi(y)、f2(y)是待定函数.由于应力函数必须满足相容方程,所以:.2..2.2二二4.x：14：d2fy一2一2一,2x二ydy代入到式O中独二/d4f(y)+d4fi(y)+d4f2(y)：:y42dy4dy4d代入到式O中4_.442_dfyx2xdfiydf2y2dfy=0dy4dy4dy4dy2考察上式可以看出它是一个x的二次方程,所以一般情况下只有两个根.也就是说只有两个位置能够满足上式.但我们对相容方程的要求是绝对满足.也就是要求在整个梁的范围内都满足.所以只有该方程的系数项和自由项全部为零.即：d4fy_nd4fiyd4f2y,d2f2y_nTOC\o"1-5"\h\z404042.20dydydydy•1-f(y)=Ay3+By2+Cy+D②fiy=Ey3Fy2Gy@d4f2d4f2y

dy4d2fydy2--12Ay-4BA5B432f2y=)yyHyKy⑷106公式中的A.、B、C…K都是待定系数.公式.3、O中分别省掉了常数项和一次项、常数项.这是由于fi(y)和f2(y)分别是应力函数中x的一次项和常数项的原因,这样处理不会对应力分量产生影响.最后求出的应力函数为：2=—Ay3By2CyDxEy3Fy2Gy--y3--y4Hy3Ky2TOC\o"1-5"\h\z2106由应力与应力函数的关系,可以求出各个应力分量：?2•:x232二x=r二一6Ay2Bx6Ey2F-2Ay3-2By26Hy2Ky232二y-Ay3By2CyDxy=-x2Ay22ByC-3Ey22FyG由于以上求得的应力分量满足了平衡方程和相容方程,所以只需根据边界确定A…K的系数,就求得了该问题的解.根据对称性,知道为偶函数,为奇函数,所以有E=F=G=0通常梁的跨度远大于梁的深度,因此上下边界是主要边界,它们必须满足.二yy4=0二yy#=-q,xyy：；h=0将它们代入的表达式,并且考虑E=F=G=0—A—BhCD-0842

,3,2hAhA—32-hAhBC=043h2A-hBC=04以上四个方程解四个未知数,求得:9将他们代回到9将他们代回到应力分量的表达式中,也就2A=-―3B=0C=-h32h有:■xy6qh32q26q34_■xy6qh32q26q34_xyq3-y6Hy2Kh36q2

二hyh3qq,—y--2h2四x2h左右两个边界,由于前面已经考虑了对称性,没有水平力.要x=l时,bx=0,考察bx所以这个仅考虑优边界.的表达式,除非q=0.而这和条件相违背.所以在这个边界上我们只能要求局部满足.运用圣维南原理运用等效力系代替它O的误差只在力作用点附近较大〕.运用的等效力系就是合成力系为平衡力系:h2〔这样产生.二xx'dy=0h~2h2.二xx^ydy=0h合力等于0合力矩等于0由第一个条件得K=0,〔奇函数在对称区间上的积分为零〕；由第二个条件得ql2

h3q

10h可以证实剪应力的合力为-ql.即h2,xyx'dy=「qlh最终求得的结果,加以整理：-x2yq-I

--xyTOC\o"1-5"\h\zhh\/由于厚度为“1〞,此时其惯性矩I=一,由于厚度为“1〞,此时其惯性矩1282任意一点的弯矩M=qll-x--l-x2=-l2-x222剪力Q--qlql-x--qx所以上式中的应力分量可以改写为:2y2yQSxy各项应力的分布bx第一项为主应力项,与材料力学中的结果完全一致.第二项为应力修正项.当L/h>4时,仅占主项的1/60;当L/h=2时,仅占主项的1/15.所以对于深梁的工程构件不容无视第二项的影响.§2.8虚功方程从前一节深梁的例子,我们可以看到,弹性力学解析求解的过程是非常复杂的.这样的求解对实际工程情况说来,很多情况根本是不可能的.所以长期以来,技术人员就一直探求数值求解的方法.有限元法是其中最成功的方法.为分析单元特性和简化分析过程,我们还需了解单元的能量关系.由于在力学上,很多时候从能量的角度分析,可以大大简化分析的过程.一,应变能的概念有材料力学知,弹性体在受到外力作用发生变形的过程中,弹性体内部会存储变形势能——应变能.在单向应力场中,单位体积的应变能的计算可以表示为：…1dU=一；二；2对于平面问题,有三个应力分量和与之对应的应变分量.由于在小变形情况下,正交力系互不影响,由力的叠加原理,所以该种情况的应变能为：TOC\o"1-5"\h\z1：,T：■1：,T：■dU=-x-x1y-y'-xyxy=-•22其中："‘」x二yxyT'：=x；yxyT整个弹性体的应变能：U=11二:=二卅二1J)Ti；,dv2v2v上式也表示应力在应变上所做的总功.二,虚功原理和虚功方程在理论力学中,我们曾经学到一个虚功原理,也称虚位移原理.其根本思想就是：假定加在物体上一个可能的、任意的、微小的位移,在平衡条件下,物体的约束反力所做虚功应等于外力所做虚功.由于能量必须守恒.在这里所说的可能的虚位移是指位移必须满足的约束条件.任意的是指位移的方向和类型是随意的.把这一原理运用到现在的弹性体中,衡量弹性体应满足的平衡能量关系就是：假定加在弹性体上一个可能的、任意的、微小的位移,在平衡条件下,弹性体内的应变能应等于外力所做虚功.同样是由于能量必须守恒.运用这一原理,我们可以推到有限元中广泛用到的虚功方程.

假定弹性体发生u、v的虚位移,那么由几何方程得:-**二-**二ux二Fx:y*cu二v,,,,,xy-Fy；:x现考察弹性体微元体和边界处微元体上的力所做虚功:⑷内部微元体上的力所做虚功左面的应力虚功-Jxdy11u右面的应力虚功左、右两面上的虚功之和〔略去高右面的应力虚功左、右两面上的虚功之和〔略去高阶微量,并考虑*x二xxudv1；x同理得剪应力的虚功之和二u二,xy*xy——udvi二y二y,rr一一,,1.r*体力X的虚功Xudv1同样的考虑Y方向的仃y以及Eyx以及体力Y的虚功,然后叠加成内部微元体上的虚功如下：dwi:-y；式yx***YvdviT.,.,x；x二y；y二y二x+Txy?xydvi⑸边界上的微元体**设斜边中点处的虚位移u、v,形心处的应力为仃xby、Txy那么在直角边上的应力和位移均有一个负的增量,如下图,虚功计算为::二xdx」,—xdy1：x2,*v**dx定一户xu-u+%%I-dyMi〔略去了局阶微量〕IL「：x2同理同理的dy直角边上的剪应力虚功为:r*一."xyui；*：xy*：:u*dyu•xydx1.n：：y2代入已有的几何关系:dx=dssin3,dy代入已有的几何关系:斜面上面力的虚功体积力的虚功同样的方法求得另一面上正应力与剪应力的虚功,全部相加即得斜边微元体上的虚功之和为：'x二.xydw2=；x7二:.'x二.xydw2=；x7,,YVdV27,1xJxxyxydV2Nfx+iX-：xcos、：-xysin-u*Y-0ysin：-.xycos-v*1支反力处的虚位移为零,所以支反力不做功,将dWi+dW2,并对整个体积积分,可以得到整个弹性体内的总虚功：**一।X-：xcos；1〞xysin「uY—'sin；一：；xycos「vdss+1V根据平衡微分方程和静力边界条件,上式的第一、第二项都是零,所以弹性体的总虚功为:***Wlz—x-xy-y-xyxyv根据能量守恒,它应与外力的虚功相等Wz：jHx；x二y；y,xy；ydv=Xu*Yv*dv-Ii：Xu*Yv*dsvvs由于该等式引入了平衡方程和边界方程,所以上式虚功方程等价于静力平衡条件〔内部和边界微元体〕.不同之处在于它是一种能量的表示形式.为了便于有限元中方便运用,引入广义力和物理方程,虚功方程变形为：；、,.*：*：F〕-〔；*书?：；加v综合以上推到过程,虚功方程表达的物理概念就是：“假设弹性体处于平衡状态,那么外力在虚位移上作的虚功应等于应力在应变上作的虚应变功,或者说等于虚应变能〞.面问题单元划分有限元法在有限元法在平面问题进行分析时,角形单元的优点是简单且对结构的不规那么边界逼近好,部的应力应变变化.这两点我们会逐渐向大家说明.或者矩形单元,三而矩形单元却更能反映实际弹性体内所以一般说来,有限元分析,单元划分的密度和单元种类选取,对计算结果起重要作用.一般单元划分越密集,结果越精确.单元多也导致求解的线性方程组阶数增高,要求计算机的内存也更大,计算的时间也越长,分析的效率就越低.解决才采用三角形单元和四边形单元、这一矛盾的方法就是在应力集中区域单元划分密集一些,应力变化梯度小的位置,划稀疏些,这样就能兼顾精度与效率的关系.一般的原那么是：根据结构的受力和支承特点,按对称和反对称的性质,简化分析模型,以减少计算分析的规模.合理布局单元的密集程度,以使计算结果精度高而计算量小.在同一单元内,单元的特性数据和材质数据应保持一致.集中载荷的作用点和载荷密度突变处应有节点.在欲知道应力状态、内力情况和位移值的位置应有节点.单元的选取欲分析的目标密切相关.模型的单元划分好后,把所有的单元和节点按一定的规律和顺序进行编号,选择适当的坐标系(直角、柱面和球面),以方便确定各节点的坐标值.点位移、节点力和节点载荷弹性体在承受外力作用后,其内力的传递实际是通过单元之间的边界来实现的.但我们把结构离散化后,如果单元划分得足够小时,可以看成为其内力的传递通过单元与单元之间的节点进行传递.对于平面问题而言,每个节点都有位移和力两个未知量,这两个量又都是x、y的函数,注意平面问题的节点是不能传递力矩的,为什么一,节点位移对三节点三角形单元而言,因有三个节点,每个节点的位移都有x,y两个分量,所以一共有6个自由度.单元节点位移向量可表示为：Je=LViUjVjUmVmT二,节点力所谓节点力,就是单元对节点或节点对单元作用的力,它是弹性体内部的作用力,也就是我们常说的内力.和上面相同的道理,节点力向量为：T?=UiViUjVjUmVmT■三,节点载荷■节点载荷是指作用在节点上的外力,包括：J⑷直接作用在节点上的外力②经等效处理后,移植到节点上的等效力.可用Xi、Yi.表示.由力平衡条件知,节点要保持平衡,那么作用在节点上的载荷应等于节点—n/.内力的合力.即:Xi=nUi,Yi=nVi,ee所有的节点载荷向量表示为：R；Fje§2.11三节点三角形单元的位移模式和形函数弹性结构受外载荷后,内部各点的位移变化规律,一般都是很复杂的.很难找到一个函数来描述整个结构内部各点的位移变化规律.但当把整个结构离散以后,在一个相当小的单元内部,却可以用简单的函数来近似描述单元内部位移的这种变化规律.而不致造成很大的误差.这就好比一条复杂的曲线,用一个函数很难描述,但在把这条曲线分段以后,对于每一条分段曲线,却可以用直线或抛物线来描述.一,三角形单元的位移模式设单元内任意一点的位移分量u、v是其位置坐标x、y的线性函数,那么：ai…a6是待定系数.u=ai…a6是待定系数.v=a4+a5x+a6y改写方程组成为矩阵的形式:单元的三个顶点i、j、m的坐标,分别为〔xiyj〔Xjyj怀口〔xmym〕,由于它们也是单元上的点,所以应该满足以上假定的位移变化规律.代入上式:UiUi=a〔a?xia3yiVi=a4a5xa6yiUjUj=a〔a?xja3yjVj=a'a§xja6yjUmUm:a-a?xia3yiVm=a4a5xma6ym解以上6个方程,求得6个待定系数.aiUi为yaiUi为yiujxjyjUmxmymxiyii।xjym-xmyj2AUi'"-.'xmyi_xiymUjxiyj_xjyiUm'1xjyj1xmym12AajUjamUmUiV\UjXiXjXmV」ViVj=2A-Vj-VmUiVm-ViUjVi-VUm-biUibjUjbmUm12Aa2同理得XiUiXjUjXmXiXjViV」2Vxm-XjUiXi-XmUjXi-XiUmXm_工一2AiuiCjUjCmUma4a5a62A

12A

12AiUiiUiiUiaiajUjbjUj,CjUj=XjVm-XmYjb=Vj-VmCi=—Xj-Xm'amUmbmUmCmUm—fijm顺序轮换,a是三角形的面积公在单元节点的顺序号i,j,m必须是按逆时针排列,否那么系数行列式是负值,而三角形的面积为负值,是不合理的.求得的6个系数可以用以下矩阵表示:■ai0aj0am01‘Uir1bi0bj0bm0Vi忸11Ci0Cj0Cm0Uj2A0ai0aj0amVja6UJ0bi0bj0bmUm0Ci0Cj0CmLVm二,形函数

将所求得的6个待定系数代入位移模式表达式中:■LJa0aj0am0"[|Ui'b0bj0bn0Vixy000];0cj0cm0uUjI001xy0ai0aj0am1Vj0b0>0bmUm0ci0cj0cmVmJ11x012Ay000也01xya6UiV1

UiV1

2Aabxcy

00bxcy0abxciyat^xciy0——fi,j,m顺序轮换.人一1,令Ni=——aibixciy2A0NiNm01：K:0Nm」lVm,ui卜人＞上式就是假定位移模式下导出的单元内任一点位移表达式.其中的插值该式的数学意义就是单元内任一点的位移可以由单元节点的某种形式插其中的插值基函数就是Ni、Nj、Nm.对于我们目前假定的位移模式是线性函数,所以得出的插值基函数也数也是类似的线性函数.由此可以看出,以也称之为位移形态函数,简称形函数.三,形函数的性质插值基函数具有反映单元位移变化形态的特征,[N]就是形函数矩阵.③单③单元内任一点的三个形函数之和恒等于1.NiNjNm=1CD在单元的三个顶点处,有i节点处i节点处Ni=1Nj=0Nm=0j节点处j节点处Nj=1Ni=0Nm=0m节点处m节点处Nm=1Ni=0Nj=0以上这些,可以通过简单的数学运算进行证实.四,位移模式收敛性的分析由于位移模式的选取是人为假定的,这种假定只能近似模拟单元内位移的变化规律,由于单元刚度矩阵的推导是以假定的位移模式展开的,那么这种假定的位移模式能否使有限元数值解收敛与精确解,在很大程度上就取决于所选的位移模式,通过数学证实,可以找出位移模式满足收敛性的几个条件：A〕完备性1〕位移模式必须包含单元的常应变状态.

将u、v的表达式代入几何方程,得：；:u；x=——=—a1a2xa3y=a2Fx；x.:v-y.:v-y=—a4a5xa6y=a6.x.:vFu一：：xy=———=—a4a5xa6y—a〔a?xa3y=a3axy:xfy：:xfy由于系数a2-a6都是常数,所以上面的应变分量也是常量.这也说明所选的位移模式中包含有弹性体的常应变状态.在上面的表达式中不含x、y的变量,说明单元的应变是常量,这也说明这种单元是一种常应力单元.2〕位移模式必须包含单元的刚体位移.弹性体位移一般包含两个局部,即变形位移和没有弹性变形的刚体位移.那么作为模拟单元位移状态的位移模式,也就应该能同时反映这两局部的位移.在③中已经证实了弹性应变的变形,下名说明他也包含有刚体位移的特征.改写位移表达式如下：u=au=a1a2xa3y=a1v=a4a5xa6y=a4a5-a3a5a3a2xy22a5-a32当%当%=%=/y=0时,由O知,a2=a6=a3+a§=0,所以上式:u二a1v—a4a5-u二a1v—a4a5-a32a5-a3x2我们再来看看一个点作刚体位移的运动方程.点M先转到Mi,再由Mi移到M2,如左图所示：u=u0-r切cosa=u0-wzy

v=v0+r6sina=v0+wzx**比拟上面的*和**式,可以看出a5=u0,a4=v0,a5-a32=Wz由此得出在三角形线性位移模式中,也包含了单元的刚体位移.B〕协调性3〕位移模式必须能够反映位移的连续性.在弹性力学求解问题时,曾经讲到过变形协调方程,也说过它是弹性体变形后仍保持连续和不发生撕裂、侵入缺陷的条件.那么位移模式的选取,也应该保证单元之间不出现撕裂和侵入的缺陷.由于上面假定的位移模式是线性函数,而两点就能决定一条直线.由于相邻单元的公共边界上的位移由与之相连的两个节点插值获得,而相邻的单元具有两个公共的节点,所以通过这两个节点所得的插值值,不可能出现不同.也就是不可能出现下列图~a和图b的情况.以上三个条件是选取位移模式必须考虑的.完备形条件是收敛的必要条件；协调条件是充分条件.在有限元中,满足完备条件的的单元是完备单元,满足协调条件的是协调单元.〔三节点三角形〕单元刚度矩阵的推导上一节我们已经建立了三角形单元的位移插值模式,并求得了形函数的方程,这样就完成了单元内任一点的位移由单元节点位移表示〔插值〕的工作,接下来运用我们已学过的一系列知识,我们就可以完成单元刚度矩阵的推导了.一,推导过程1〕由位移插值函数导出单元应变的单元节点位移表达式「£

「£

除0的0Ni曳|ex0Ni0NNxN曳|ex0Ni0NNxNjx0-Nj0Nj：:y.Njx:NNm：:y：Nm二y：NmUi、vmj,1,、••一将Nj=—(ai+bix+ciyM弋入上式,可得:2AUibi0Uibi0bj0bmtie10zf0——0c0cj02A!jCibCjbjCm0]cm'bmViujvjum1Vm01Cibii=i,j,m如01Cibii=i,j,m-bi{/=BiBjBm^}e其中fei0j2A矩阵[B]称为应变矩阵,或称为几何矩阵.由以上计算公式知,它与单元的节点坐标有关,但不随点的坐标变化,就是说在这一单元内所有点的应变是相同的.2〕求得单元应力的单元节点位移的表达式将几何矩阵代入单元的物理方程,就有：J)e二D；)e=D1B:e弹性弹性矩阵[D]是由材料常数组成的矩阵.令[S]=[D][B],代入平面应力的物理方程,就有七产"St1e=—E71-」21sL2152A■七产"St1e=—E71-」21sL2152A■bi1七

1-NrC写成分块矩阵的形式rbi%1一」Ci001-」2」CCi

1-<

-Vbi」CiC

1-\

Tbi12Abj」bj1—J-bi10:G—Cj2jbj0CjCj1—Jbji=i,j,m0Cjbjbm

0CmbmJbm1-JC2Ui1%Cm1-\bm2一也可以3〕用虚功方程导出单元刚度矩阵〔单刚矩阵〕虚功方程%*,.｝=1｛屋yd及｝dvV假定单元的厚度为t,上式改写为单元的虚功方程形式,

也*?T{f}e=J[Q*TIdLltdxdy虚应变也可以用几何方程表示qf=ib七eq*tT=〔fefc*ri=〔fe*y>fei代入上式隹*FT{FF=H(%*}eTlB『bIfek}etdAA由于虚位移元素是常量,所以可以提到积分号以外,并与左边的消去〔为什么〕.于是上式变为：『产=BTbIBLVtdxdy,A\令ke=以BTb】Btdxdy,虚功方程就成为了单刚方程由于［B］、［D］都是不含有x,y的常数矩阵,所以双重积分实际就是对面积积分了.ke=BTbIBtAa——是三角形面积将前将前面求得的应变矩阵和弹性矩阵代入,然后作矩阵乘法.就得到我们要求得的矩阵计算公式.我们在这里采用分块矩阵的方法记忆.2kj^^41-J22kj^^41-J2AHbsCr1-」—2"CrCs1-brCs2,,1-<JbrCs——bsCr21-」--2-brbsr,s=i,j,m注意以上是平面应力状态的单刚矩阵,如果是平面应变问题呢二,单元刚度矩阵的性质单元刚度矩阵［K］e表示单元反抗变形的水平.它与通常的弹簧刚度系数k的物理意义本质相同,只不过［K］e是一个6X6阶的矩阵.共有36个元素.这是由于三角形的节点力向量和节点位移向量都为6的缘故.1〕物理意义分块矩阵［kj］表示的物理含义是：节点j处产生单位位移,而节点i、m被约束,此时在节点i处产生的节点力.我们写出它的4个分量元素：12ij22在上面的矩阵中,元素的第一个下标表示产生节点力的节点,第二个下标表示产生节点位移的节点.上标“1〞表示水平分量,“2〞表示竖直分量.而且上标和下标的关系是对应的.也就是说第一个下标对应第一个上标；第二个下标对应第

二个上标.如此就有:12ki；2就表示节点j处产生竖直单位位移,在节点i处产生的水平方向的节点力.2〕单元刚度矩阵是对称矩阵一〔keT二k?这一点可以通过简单的数学证实如下：〔kH=〔BTb】B[aT=BTb1伯TT二bTDIB[a二k7单元矩阵的对称性,从物理学角度反映出的道理就是,“功的互等〞.也就是在节点j处产生某一位移引起节点i处的节点力,应等于在节点i处产生相同位移引起节点j处的节点力.3〕单元刚度矩阵中的元素只与单元的材料性质、几何形状、尺寸大小有关,而与单元的位置无关.单元刚度矩阵中不含有araj、am,上节对位移模式收敛性分析中,曾经说明了ana4分别表示单元的平移分量,而1a；=--aiUi-ajUjamum2A1a4二一aiViajVjamVm2A由上式知单元的平移运动与ah小am有关,而[K]e又与ar电、am无关,所以说它不随坐标轴的平动而变.4〕单元刚度矩阵是奇异矩阵,即KT=0从力学的角度理解单兀的刚度方程-K升疗,当给定位移时,可以求得力；当给定力时,却不能求得位移.由于[K]e不存在逆矩阵,在单元没有给出任何约束的情况下,除有应变的可能性外,还同时有刚体位移的可能性.所以方程无解.荷的节点移置前面对于有限元模型的分析时,曾经说过,单元之间的力传递是通过节点进行的.所以不在节点上的力,必须按静力等效原那么,把它们移置到节点上.静力等效原那么：原载荷在任何虚位移上所做的虚功,与移置到节点上的节点载荷所做的虚功相等.这种处理方法,和我们前面讲到的圣维南原理相同.它们只会影响模型局部的应力分布,而不会影响整个结构的应力.下面根据力的类型,分类说明处理的方法.1〕集中力的移置如下图的受力示意图,M处的力为MbhPyT,移置后的节点载荷为{r*=XiYiXjYjXmym】T,虚功相等就是：{・=TM,=「：*P〕•••忖〕=N'ide,所以有〔和*]£r}=〔NKsMe"p}=〔fe*}eTInI{p}・.：R」NHP〕由于我们现在一直局限于单元的载荷移置,所以上式中的{R}应记为{R}eo如果将它写成分量的形式：{R}e=XYXjYjXmYm『=NiRNiPyNjBNjPyNmPxNmPy1从上面公式可以清楚的看到,载荷移置结果与单元位移模式密切相关.2〕体力移置体积力密度为{p}=反Y1,将单元中微元体的体力{pXdxdy看作集中力,那么=N■ptdxdy一,A3〕面力的移置设在单元的一个边界上作用有分布力,面力的密度为{p}=!XY『,将微面积tds上的面力也看作是集中力,那么有R=sN1ptds以上面力,体力的积分运算较为复杂,在线性模式下,可根据理论力学中静力学平行力分解的原理,直接求得等效的节点载荷：③均质等厚的三角形单元,受W的体力,那么三个节点分别受到W的节点载荷.3⑶i、j边受到集度为q的均布面力载荷,那么i、j边节点受到邺的等效节点载荷.21③i、j边受到线性分布的面力,i处集度为零,j处集度为q其合力为R=-qlt,那么i2一一一1一一2处的节点载荷为1R,j处的载荷为2R.33这些简单的规那么,对于今后实际中的应用,可以提升效率.§2.14整体分析该过程是将离散别离的各个单元组集成离散的结构物,从而建立模型的总刚方程.一,总刚方程和总刚矩阵的组集1,总刚度矩阵的组集原那么A,整个离散结构变形后,各个单元在节点处仍然协调地相互连接.即环绕某个节点的n个单元,在节点i处具有相同的位移.数学公式的描述：W[=如¥="•・=也F=也}B,各个节点应满足静力平衡条件.即每个节点上的节点力合力应等于该节点的节点载荷.

数学公式描述:此处的£代表环绕节点i的所有单元节点力求和.e2,在该原那么指导下,实例的组集过程如图a所示的平面问题,采用图b的方法划分网格,单元分析完成后,现将它们组集成原问题的有限元离散网格,求该问题的总刚矩阵.3为i,j,m=1,23为i,j,m=1,2,111单元G〕节点号〔i节点从最小号开始,然后逆时针排列节点号〕,单刚方程：kkk将它们展开311汨：1二忆/⑥)1•K2k2』•k3k3,下2、卜211Y・卜22匚2』IN"」22112222335j31aJ•K32k2?-K331C.3?完全相同的道理,得其他三个单元展开的单刚方程单元Q节点号为i,j,m=1,3,4KJKdYEJKlJJkHJkF、〕2单元单元o节点号为i,j,m=3,5,4g3=K33k3：K35k513K34泞/『5：'3=K53HG3-K55Hf3・K54k4?下d二K43HdK45HJK44N/单元Q节点号为i,j,m=2,3,5X二K22"、"K23S"K25k.5)5M=K32k2：4K33I〞/K35跳J(F514=K52"、LK"、/K"/实际上,刚度方程可以写成如下的一般形式：{Fj=K工卜［瓦}k表示单元的节点组成k4,j,m首先我们运用规那么.2,即各个节点力的合力等于节点载荷“,括:e=Me所以也就有：；FJ,岸二口){f2}2+k2}4={r2}匕明岸『3~g3二代)下4)2•尸4节=R:'仁手•Is)4=g:'我们将上面得到的方程按这个规律相加：匕汜：2=4：'=匕1工：屋k2)1413忆3)1—2."忆3匕3〞忆4忆4)2对这个等式运用规那么.i,即%了=也F=…=电}n=%}并整理有：叱Kii1KI•〔KM%".2"：,代)=忆1〞「靖工仁：'"小14.)％"一代：'=^311-K312''ll卜321,K32I'J)k331・K3312k333k334",K342K3433-^35-K354X代)二k1忆JK4312K43M3"K4412K4414；忆585)

应｝=k52M&｝+收537+卜531近｝+-543J+fKJ+1^55/砧5｝将上述的这些方程写为矩阵的形式,就得到了该问题的总刚方程为:RiR2」RRiR2」R3,=R4R5「ki产Ik2ilIkJ-0ki21

卜22F

卜32「0

k524kJ2kJk23140k33l2-23匕3卜3-2324k54『卜35/IS'1〞与花3・占40J我们对这个总刚方程进行分析,可以得出以下的规律3,总刚矩阵的规律得到结构的节点载荷,节点载总刚方程就是运用所有节点的位移向量和总刚矩阵相乘,荷向量〔外界对弹性体作用的载荷〕很容易求得,节点位移向量得到结构的节点载荷,节点载阵如何得出.如都根据上述的步骤,成千上万个节点,该如何作呢所以有必要总结总刚度矩阵的规律..当r=s时,即子块矩阵［Krs］是总刚矩阵主对角线上的子块矩阵.由环绕节点r〔s〕各单元刚度矩阵的相应对角线子块相加.②当rws时,但r、s是相邻单元的公共边上的节点时,子块矩阵［Krs］等于两相邻单元刚矩阵的相应子块矩阵相加.如上例中的［K13］、［K31］等.③当rws时,且r、s只是一个单元的边上节点时,子块矩阵［Krs］就等于该单元刚矩阵的相应子块矩阵.如上例中的［K12］、［K21］等.②当rws时,且r、s不属于一个单元的边上的节点时,子块矩阵［Krs］等于零.如上例中的［K15］、［K51］等.利用上述规律不仅可以检验总刚组集的正确与否,而且可以直接组集总刚矩阵〔手工组集〕.但这种方法不利于计算机编程.下面介绍计算机组集总刚矩阵的方法4,总刚矩阵组集的步骤③扩大各单元刚度矩阵,使之成为与总刚矩阵相同的阶数.②根据总体节点的编号,将各单元刚度矩阵的各个子块移到相应的位置上.其余位置充零.⑶把各个改造过的矩阵直接相加,就得到总刚矩阵.我们以单元.4为例讲解步骤.1、Q单元4的i、j、m=2、5、3,所以单元刚度方程是K25K55K3K25K55K35h's"伟2川K4k52K4K33电｝4值｝43-4.1,.4・222532kkkTOC\o"1-5"\h\z000000k224k2340K2J0K32TK3310K35700000-0K524K5340叱557_二,总刚度矩阵的性质总刚度矩阵由单元刚度矩阵组集而成,所以也具有单元刚度矩阵的一些性质,如相同的物理意义,位置无关性、对称性和奇异性等,还具有以下性质：1〕稀疏性由规律4知,当rs,且r、s不属于同一单元白两个节点是[K「s]=0,说明互不相关的节点数愈多,零子块矩阵就愈多.一般说来相关节点数不超过9,而整个分析对象常常成百上千、上万或几十万.如果整体有100个节点,那么百分之九是非零子块,而百分之九十多都是零子块.所以在总刚度矩阵中非零的子块是很稀少的.2〕带状性所谓带状性,就是指总刚矩阵中的非零子块,集中分布在主对角线的两侧,呈带状分布.半带宽值就是计算这带状宽度的数值.它是以排列元素最长的一行,从第一个非零元素起至主对角线元素止,所有元素的个数.其数值可以由节点的总体编号算出.B=〔相邻节点号的最大差值+1〕X〔单个节点的自由度个数〕对于平面问题说来,就是：B=〔相邻节点号的最大差值+1〕X2下列图中的两种编号方式,可以的分别计算其半带宽值:图a,B=〔2+1〕x2=6图b,B=〔6+1〕X2=141士―彳国身彳国身凶1厂71工半带宽值影响计算机存取总刚矩阵所需的内存大小.愈小愈好.由于半带宽值是直接受节点总体编号的影响.所以我们在建立有限元模型时,应慎重考虑,采用优化的方法.〔现在通常在程序中都配有这样的优化程序〕.三,总刚矩阵的压缩存取技术由于总刚矩阵具有对称性,所以我们只需存入主对角线上半带或下半带的元素,就可以完成解方程的运算.此即为总刚矩阵的半带宽存储方法.假定取主对角线上半带元素存储,具体做法就是每行元素以主对角线上的元素开始,存储每行半带宽数值的元素个数.如此各元素的行号不变,改变的只是列号.新列号和原列号的关系式如下：新列号=原列号-行号+1前面图示的例子,如果采用图示表示总刚矩阵的存储,就是：如下列图所示的情况.可以看出原总刚矩阵存储需要24X24的2维数组,改为半带宽存储后,就成为了24X6的数组了.

另有一维的压缩存储方法,比这种二维的半带宽存储方法压缩更多.在这儿我们就不做介绍了.有兴趣的同学可以看相关的参考书籍.讲到这儿,实际我们已经知道,每个元素在总刚矩阵中的位置,在节点编号完成后,就完全确定下来了.所以计算机对总刚矩阵的计算,是一部完成的.即“边对号移置、边改列号、边累加〞.四,总体边界条件的处理前面介绍总刚矩阵的性质时,说明了总刚矩阵是奇异矩阵.即|k】二0.就是说总刚矩阵不存在逆矩阵.要求得节点位移的位移解,还必须引入边界条件.1,边界条件的类型①节点固定,即Ui=Vj=00给定节点位移值.即Ui=Ui,Vj=Vj2,处理方法1置“0、1〞法即首先在总刚矩阵［K］中与位移分量相对应的行和列元素改为0,但主对角线上的元素改为1,然后在节点载荷向量的列阵中,与对应元素的位移用代替,其余元素减去位移分别乘［K］中相应元素.例：某结构的总刚方程为卜］0乳0缶｝10刈=010M,3=5、v2=c4,根据上面的方法修改如下：如果展开上式,马如果展开上式,马上就有U1=c1、V2=C4第二行k22V1k23U2k25U3k210V5二丫1一k21c1一k24c4一1000…0]lrU1L0k22k230k210V1丫1—k21c1—k24c40k32k330k310U2>=<X2-k31c1-k34c4>0001…0V2c4iaaasi-a——■■■■■■■■■■■■■■J0k102k1030k10101V5,丫5—k101c1—k104c4.移项k21.X22V1k23U2k24c4k25U3k21V5=丫1从中可以看出,这样处理后不仅可以直接得到U1=c1、V2=c4,而且它们产生的效果也计入到了其他方程中.如果是固定位移情况,那么载荷向量中的变化是什么呢〔好好想想〕2乘大数法在总刚矩阵［K］中,把与位移位移相对应的行与列主对角线上的元素乘一个很大的数101°,然后把载荷向量中的对应元素代以给定位移乘以相应主对角线上的元素,再同样乘以一个很大的数.同上例./-1/-10,,,,「父10kI2k〔3k〔4…ku.C[k[1M10k21k22k23k24…卜210V1Yk31k32k33k34卜310U2X2,,,,“10,<>=<k41k42k43k44M10…k410V2C4k44父10amm*mm■.ak101k102k1030k1010V5丫51010、、考察第一个方程:k1k12V1k110V5=c1k111010⑤求解⑤求解应力上乒=D.Lje-DI〔b,’e两边同除以"黑1010,由于k11x1010>>k1j,所以U1=c1同理可得V2=C4,〔思考该方法不能用于固定位移的情况,为什么〕必须说明的是,当总刚矩阵采用半带宽存储的时候,以上边界条件的引入也应当根据半带宽存储的格式修改,否那么就是错误的.实际应用中,当是第一种情况时,采用方法.1；当第二种情况时,用方法O2.五,应力计算总刚方程在引入足够的边界条件后,就可以进行求解了.所谓足够的边界条件就是要使系统是一个几何不变的系统.数值解方程的方法大多采用高斯消元法.当求解出位移向量以后,就可以通过单元的几何方程求解应变,通过物理方程求解应力了.这个过程称之为回代过程.@求解应变行=B㈠由于三角形单元是常应力与常应变单元,所以由上述方法求得的应力或应变看作是形心处的应力或应变.而且很明显,这样求得的应力或应变是不连续的.为了推算弹性体内某一点的接近实际应力,我们通常采用以下两种方法,来平滑应力的突变.①绕节点平均法一一就是将绕某一节点的各单元形心应力加以平均,来表示该节点的应力.如左图所示.二x1J二Ax.二Bx12二--1i0AxBx,■.C,■.D,■.E,1.Fx22xxxx该方法在各单元面积相差不大,单元的内部节点时较为准确.而在外部边界节点上,那么不理想.所以对于边界节点多采用三点插值的方法求得.也就是用内部的节点来推算.,两单元平均法一一把相邻两单元的应力相加平均,表示两单元公共边中点的应力值.

Ox2二；二晨,二Bx边界边上的应力值也采用插值方法得到.卜面通过一个例子来看看,具体至此我们对于有限元求解平面卜面通过一个例子来看看,具体如何计算的.题计算,有限元法求解平面问题的步骤:1〕根据分析对象和给定的条件,根据一定的比例绘出计算简图,该图应有各部位的尺寸、外力和支承情况.2〕选择适宜的坐标系,划分单元准备数据.数据包括节点数、编号和坐标,单元数、编号和节点组成；材料特性常数；载荷大小及位置；边界条件.3〕计算单元的面积、应变矩阵、弹性矩阵和单元的刚度矩阵.4〕组集总刚矩阵,处理载荷移置,形成载荷向量.5〕处理边界条件,修改总刚矩阵,得到最后的刚度方程.6〕求解线性方程组,得到位移向量.7〕求解应力值.三,例题设有如下图的正方形薄板,在对角线上作用有沿厚度均匀分布的载荷,其合力为2kN,板厚为1个单元厚度,为使计算简便,假定材料的科=0,求该板的变形.分析：该薄板对称于对角线,受力也是如此,所以可取其1/4来计算.由于是取1/4,所以该模型的两个直角边,都是对称线.X的对称线不能有x方向的位移,y轴的对称线不能有y方向的位移.直角的顶点既是x轴对称线上的点,又是y轴对称线上的点,所以该点在x、y方向的位移都应该约束.最后得出的计算模型如下.解：1准备的数据A单元节点数据坐标123456X001012Y211000

B单元信息数据①鸣⑷⑷I1223J2455M3536所有单元的板厚是“1〞C材料常数弹性模量E,泊送比科=0D载荷数据00000000TE位移边界1、-10Vl0v2u3V3E位移边界1、-10Vl0v2u3V300u50u601T2计算单元刚度矩阵1xiyi一.i.单兀1面积A=-1x2y221X3y3b1=y2-y3=1-1=0C1=-X2X3=01=1形成的几何矩阵b1B1=1020C1b1C1弹性矩阵10单元1的刚度矩阵111210201=0.513b2=y3-y1=1-2=一1b3=y1-y2=2-1=1c2--x3-x1--1-0--1C3=一x1x2=00=0b20b3000-1cc1Lcc0c20c3=0000..2屋一八,c2b2c3b310-1010-100-1010.0.50-0.5010k1=B「b】B1A=E-0.501.52|-0.5-10.500-10.500.5-1000.5-1-0.51.50-0.5010「0「0.50-0.5-0.500.5同理可求得k2、k产和k4.需要说明的是可以利用位置无关性直接得出kI2和k了,1000.5k3=E.0.5200由于它们是i1000.5k3=E.0.5200TOC\o"1-5"\