西安交通大学有限元原理与工程应用内部教材

第一章引言§1-1概述1、有限元方法（TheFiniteElementMethod,FEM）是计算机问世以后迅速发展起来的一种分析方法。众所周知，每•种自然现象的背后都有相应的物理规律，对物理规律的描述可以借助相关的定理或定律表现为各种形式的方程（代数、微分、或积分）。这些方程通常称为控制方程（Governingequation）o针对实际的工程问题推导这些方程并不十分困难，然而，要获得问题的解析的数学解却很困难。人们多采用数值方法给出近似的满足工程精度要求的解答。有限元方法就是•种应用十分广泛的数值分析方法。图1-1工程问题的求解思路有限元方法是处理连续介质问题的•种普遍方法，离散化是有限元方法的基础。然而，这种思想自古有之。齐诺（Zeno公元前5世纪前后古希腊埃利亚学派哲学家）曾说过：空间是有限的和无限可分的。故，事物要存在必有大小。亚里土多德（Aristotle古希腊大哲学家，科学家）也讲过：连续体山可分的元素组成。古代人们在计算圆的周长或面积时就采用了离散化的逼近方法：即采用内接多边形和外切多边形从两个不同的方向近似描述圆的周长或面积，当多边形的边数逐步增加时近似值将从这两个方向逼近真解。图1-2可以用来表示这一过程。

近代，这•方法首先在航空结构分析中取得了明显的效果：一种称为框架分析法(frameworkmethod)被用来分析平面弹性体(将平面弹性体描述为杆和梁的组合体)(1941,HrenikofT)；在采用三角形单元及最小势能原理研究St.Venant扭转问题时，分片连续函数被用来在子域中近似描述未知函数(1943,Courant)o此后，本方法在固体力学、温度场和温升应力、流体力学、流固耦合(水弹性)问题，以及航空、航天、建筑、水工、机械、核工程和生物医学等方面获得了广泛的应用。从而，促成了一个内容十分丰富的新兴分支——计算力学的出现，长期以来在力学中存在的求解手段落后于基本理论的现象得到了根本的扭转。由于拥有了强有力的分析手段，相比之下对物质世界本身(例如本构关系)的了解反而出现了一些新的薄弱环节。有限元方法的第二个关键时期出现于二十世纪六十年代中期，归功于Argyris,和Kelsey(1960)以及Turner,Clough,Martin和Topp(1956)。然而，"有限单元”是由Clough首次提出的(I960).在众多数学家的共同努力下，有限元方法的基本原理被揭示以后，这种方法摆脱了各种各样的工程背景而成为一种具有普遍意义的数学方法。这样就不仅极大地扩展了该方法的应用范围，而且拓宽了人们的思路，在构造方法时人们不再受工程直觉的束缚。2、众所周知，一个连续体有无限多个自由度(属于无限维空间)，有限元方法则是将它转化成一个有限自由度(属于有限维空间)，建立有限元方程，求其近似解。可以将有限元法理解为在子域内应用的瑞利―里兹法(Rayleigh—RitzMethod).在传统的瑞利―里兹法中，必须假定近似的位移函数和其各阶导数在整个求解区域内有良好连续性。然而，实际的工程结构往往比较复杂。例如，变压器的箱体可以看成是由板和梁的组合结构；管道系统中的阀门、接头和三通表现为集中质量。在数学的描述上，这些实际情况表现为间断点，在这些部位函数的导数(及应变)是不连续的。因此，瑞利―里兹法的工程应用受到了限制。另外，对于二维及三维的工程结构，如果其几何边界不规则，要寻找满足边界条件的连续的近似位移函数是极其困难的。在有限元方法中，由于利用了分片插值技术，连续体(区域)的形状可以不受任何限制。而这一难题正是以前其他分析方法所难以克服的。图1-3给出了有限元法与传统的有限差分法在描述同一对象时的比较。有限差分法有限单元法有限差分法建立有限元方程大体有三类方法：（1）直接方法这种方法是直接从结构力学引伸过来的，作为一种建立有限元方程的方法而言，只在简单情况下才能凑效。这种方法的优点在于简单、易于理解，一些基本概念和作法的物理意义清晰，对理解有限元方法的相关概念和具体作法十分有益。（2）变分方法这种方法是讨论有限元方法时最常用的一种形式。有限元方法最早的严格理论论证就是以这种形式给事的。变分方法主要用于线性问题，该方法要求被分析的问题存在一个“能量泛函”，由泛函取驻值建立有限元方程。对丁•线性弹性问题就表现为最小势能原理、最小余能原理或其他形式的广义变分原理。对•于某些非线性问题（弹塑性问题）的虚功方程也可归于这一类。本书主要利用这种形式。（3）加权残值法对于线性自共瓶形式方程，加权残值法可能和变分法得到相同的结果，这将使我们得到一个对称的刚度矩阵。对于那些“能量泛函”不存在的问题（主要是一些非线性问题和依赖于时间的问题）加权残值法是一种很有效的方法。伽辽金Galerkin）法（即，选形函数为权函数的加权残值法）属于这一类。3、计算机和应用程序要用有限元方法的理论来解决实际问题离不开计算机（硬件）和程序（软件），人们大体要完成以下四方面的工作：（1）数据储存应用有限元方法求解实际问题时，在计算过程中要存贮大量数据（原始数据、中间数据和最终结果）。对于一个中等规模以上的算题，数据量相当可观。例如，一个不到500个结点的板壳结构（中等规模）的算题，要占用18MB的外存空间，在处理规模稍大的算题时一定要有足够的外存空间。（2）数据管理为了充分利用存贮空间，编制程序时要注意到存贮空间的利用率。（3）数值计算计算成本主要取决于数值运算的时间，尽量选用先进的计算方法，提高求解效率。（4）前处理及后处理为了减少人工准备原始数据的工作量，程序要有尽可能完善的“自动生成“功能。由程序产生一部分原始数据。尽管如此，对于中等规模以上的算题来说，准备原始数据仍然是一件繁重的工作。数据是否没有差错，往往决定着一次上机的成败。分析计算结果也是一件繁重的工作，一个好的程序应有较完善的后处理功能。例如，将计算结果绘制成图形或曲线。有限元分析软件可分为通用软件和专用软件，下面的表格简略的介绍了一些近年来国内外著名的分析软件，供读者参考。编号软件名称二维杆单元=维梁单元平面应力单元平面成变单元薄膜单兀剪切板单元弯曲板单元薄壳单元厚壳单元旋转克中元轴对称体单元三维体单元边界单元间隙单元哑单元标量单元刚性单元1ANSYSLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGM2ADINALGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGM3MARCLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGM4NASTRANLGLGMLGMLGMLLGLGLGLGMLLLGMLLGM☆☆☆5ABAQUSLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGM6FENRJSLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGM7PAFECLGLGLGMLGMLGMLGMLGMLGMLGMLGMLGMLLGM8ASKALMLMLMLM1.MLMLMLMLMLMLMLM9EALLGLGLGLG1GLGLGLGLGLG10SAMCEFLGMLLGMLGMLLLLLGMLGMLGMLGMLGIILARSTRAN80LGMLGMLGMLGMLGMLGMLGMLGMLGMLGMLGM12HAJIF系列LGMLGMLGMLGM1GMLMLGMLGMLGMLLLGML注：L一线性、G—几何非线性、M—材料非线性表1-2有限元软件的应用领域编号软件名称线性静力分析同仃振动分析屈曲分析非线性静力分析后屈曲分析非线性接fell分析非线性振动分析非线性瞬态响应分析波的传播分析流体与构耦合分析热叮机械耦A口分析粘塑性分析动力粘塑性分析静力敏度分析固有特性敏度分析屈曲敏分析动态响应敏度分析iANSYS☆☆☆☆☆☆☆☆☆☆2ADINA☆☆☆☆☆☆☆☆☆☆☆☆☆3MARC☆☆☆☆☆☆☆☆☆☆☆☆4NASTRAN☆☆☆☆☆☆☆☆☆☆☆☆☆5ABAQUS☆☆☆☆☆☆☆☆☆☆☆☆☆6FENRJS☆☆☆☆☆☆☆☆☆☆7PAFEC☆☆☆☆☆☆☆☆☆☆☆☆8ASKA☆☆☆☆☆9EAL☆☆☆☆☆☆☆☆☆☆☆☆☆☆10SAMCEF☆☆☆☆☆☆☆☆☆11LARSTRAN80☆☆☆☆☆☆☆☆☆☆12HAJIF系列☆☆☆☆☆☆表1-3有限元软件的载荷计算功能编号软件名称集'I1找荷纹载荷对称我荷面我荷体积我荷力我荷初始变形热我荷离心我荷跟随力1活我循环我荷随机载荷陀螺载荷非比例我荷瞬态冲击我荷接触我荷IANSYS☆☆☆☆☆☆☆☆☆☆☆☆☆☆2ADINA☆☆☆☆☆☆☆☆☆☆☆☆☆☆3MARC☆☆☆☆☆☆☆☆☆☆☆☆☆4NASTRAN☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆5ABAQUS☆☆☆☆☆☆☆☆☆☆☆☆☆☆6FENRJS☆☆☆☆☆☆☆☆☆☆☆☆☆7PAFEC☆☆☆☆☆☆☆☆☆☆☆☆☆☆8ASKA☆☆☆☆☆☆☆☆☆☆☆☆☆9EAL☆☆☆☆☆☆☆☆☆☆☆☆☆☆10SAMCEF☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆11LARSTRAN80☆☆☆☆☆☆☆☆☆☆☆☆☆12HAJIF系列☆☆☆☆☆☆☆☆☆编号软件名称各向同性正交各向异性多层材料非均质材料热弹性热塑件线弹性非线性弹性弹邛性弹性力硬化粘弹性粘塑性高温端变1ANSYS☆☆☆☆☆☆☆☆☆☆☆☆☆2ADINA☆☆☆☆☆☆☆☆☆☆☆☆3MARC☆☆☆☆☆☆☆☆☆☆☆☆☆4NASTRAN☆☆☆☆☆☆☆☆☆☆☆☆☆5ABAQUS☆☆☆☆☆☆☆☆☆☆☆☆6FENR1S☆☆☆☆☆☆7PAFEC☆☆☆☆☆☆☆☆☆☆☆☆8ASKA☆☆☆☆☆☆☆☆☆9EAL☆☆☆☆☆☆10SAMCEF☆☆☆☆☆☆☆☆☆11LARSTRAN80☆☆☆☆☆☆☆☆☆☆☆12HAJIF系列☆☆☆☆☆☆☆☆☆表1-5有限元软件的模块化功能编号软件名称固定格式数据riIII格式数据面向问即的语言用户安排的数据顺序系统安排的数据顺序节点坐标自由生成单元连接信息自由生成约束对称性边界条件口动生成子结构连接信息自动生成相同部件的h动重复生成载荷fl动生成维单元信息的自动生成角形单元的自动生成四边形单元的II动生成旋转体壳单元的自动生成维实心单元的自动生成1ANSYS☆☆☆☆☆☆☆☆☆☆☆☆2ADINA☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆3MARC☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆4NASTRAN☆☆☆☆☆☆☆☆☆☆☆☆☆☆5ABAQUS☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆6FENRIS☆☆☆☆☆☆☆☆☆☆☆☆☆7PAFEC☆☆☆☆☆☆☆☆☆☆☆☆☆8ASKA☆☆☆☆☆☆☆☆☆☆☆☆☆9EAL☆☆☆☆☆☆☆☆☆☆☆☆☆10SAMCEF☆☆☆☆☆☆☆☆☆☆☆☆☆11LARSTRAN80☆☆☆☆☆☆☆☆12HAJIF系列☆☆☆☆☆☆☆☆☆☆☆☆表1-6有限元软件的前后处理功能编号软件名称维单元ri动生成相贯线自动生成数据检作的交国形技术未变形结构的网格图形'|1<示窗口显示隐藏线消除轴侧视图透视图剖面图节点与单元重新编号数据的列表检杳法固定格式的列表输出由用户确定的数组和顺年输出最大lj最小吊输出按节点分区的最•大值与最小值输出应力温度流量极限输出用于用户后处理的文件1ANSYS☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆2ADINA☆☆☆☆☆☆☆☆-V☆☆☆☆☆3MARC☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆4NASTRAN☆☆☆☆☆☆☆☆☆☆☆☆☆5ABAQUS☆☆☆☆☆☆☆☆☆☆☆☆☆6FENRIS☆☆☆☆☆☆☆☆☆☆☆☆☆☆7PAFEC☆☆☆☆☆☆☆☆☆☆☆☆☆8ASKA☆☆☆☆☆☆☆☆☆☆9EAL☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆10SAMCEF☆☆☆☆☆☆☆☆☆☆☆☆J1LARSTRAN80☆☆☆☆☆☆☆12HAJIF系列☆☆☆☆☆☆☆☆☆☆☆☆编号软件名称温度、流量与应力等直线绘图面函数绘图用户选择项目绘图1按元素或区域时间历程绘图参数规定奇异性检查错误修正对用户的服务热线电话咨询培训录像带信息报刊开发单位编程语言（程序规模）1ANSYS☆☆☆☆☆☆☆☆☆☆☆☆SwansonAnalysisSystem（美国）F0RTRAN77（150000行）2ADINA☆☆☆☆☆☆☆ADINA工程公司（美国）FORTRAN（150000行）3MARC☆☆☆☆☆☆☆☆☆MARC公司（美国）F0RTRAN4F0R66F0R77（100000行）4NASTRAN☆☆☆☆☆☆☆☆☆☆NASA（美）主持MSC公司（美）开发F0RTRAN4Assembler（600000行）5ABAQUS☆☆☆☆☆☆Hilbitt,Karlson,andSorensen公司（美）F0RTRAN77（140000行）6FENRIS☆☆☆☆☆☆☆NTH,SINTEF（挪威）F0RTRAN77（160000行）7PAFEC☆☆☆☆☆☆☆☆☆PAFEC公司（美）FORTRAN（400000行）8ASKA☆☆☆☆☆☆☆斯图加特大学静动力学研究所（西德）F0RTRAN4（600000行）9EAL☆☆☆☆☆☆☆☆E1SE公司（美）ANSIF0R66Assembler（150000行）10SAMCEF☆☆☆☆☆☆L.T.A.S（比利时）FORTRAN4（300000行）11LARSTRAN80☆☆☆☆斯图加特大学静动力学研究所（西德）F0RTRAN4FORTRAN77（200000行）12HAJIF系列☆☆航空部（中国）FORTRAN4（280000行）至今，有限单元方法已有50多年的历史，尽管研究工作目前仍在继续进行，应该说该方法已经是•种成熟的分析方法，而且已经十分普及。在许多高等院校，该方法已成为研究生和高年级本科生的选修课或必修课，众多的工程技术人员已经将该方法作为常用的分析工具。§1-2典型问题在有限元分析具体问题时，其基本未知量为某种场变量。场变量可以是标量（温度t）,也可以是矢量（位移3,E,日）。（最后处理的实际上都是标量）每个具体的工程问题总有一个确切的定义域空间：绝大多数具体的工程问题是二维或三维的（二维往往是对三维的筒化）。另外还可以将问题分为有限维问题和无限维问题，例如:有限维：曲轴、飞机构件、汽轮机叶片；无限维：电场、磁场、声场、水坝。适用有限元方法分析的工程问题中最具代表性的问题是力学问题，由于有限元方法的问世源于力学问题，因此许多概念带有明显的力学特征，其次，力学问题代表性强，力学量的矢量特性。共粗量（位移场、应力场）直观性强。（1）力学：位移场、应力场、速度场（2）传热学：温度场（3）电学：电场、磁场（4）流体力学：流场（5）声学：声场§1-3有限单元方法处理问题的基本步骤(文本或图形)

输出结果(文本或图形)

输出结果图is有限元方法处理问题的流程1、将给定的区域离散化为子区域（单元）的集合。离散化的目的是使问题的性质在每一单元内尽量简单。一般情况下，单元内部不能存在任何间断性。离散化的目的另一作用是使单元的几何形状尽可能的吻合实际问题的几何边界。（i）用预先选定的单元类型来划分求解域，创建有限元网格；（ii）给单元及节点编号；（iii）创建几何特性（例如：坐标系，横截面的面积等）。2、对有限元网格中现存的各种典型单元进行单元分析。利用各种方法形成单元的刚度矩阵、载荷矩阵及质量矩阵（动力分析需要质量矩阵），这里就需要选择近似的插值函数（位移模式），在直角坐标系中通常采用多项式函数，在圆柱坐标系中则常采用三角函数和多项式函数的混合形式。山于同类的单元可以采用相同的位移模式，因此只需对典型的单元进行单元分析。(i)对各典型单元创建与其微分方程等价的变分形式；(ii)假设典型的独立变量(Trialfunction)(例如〃)的形式为/=1(iii)将上式代入变分形式，并得到卜的=\f\(iv)推导或选择单元插值函数已计算相关的单元矩阵。3、将单元方程合并为总体方程组(i)给出局部自由度与总体自由度之间的关系(此关系反映了基本变量在单元之间的连续性或单元之间的连接性)：(ii)给出二阶变量之间的“平衡条件”(局部坐标系中的力分量与总体坐标系中的力分量之间的关系)；(iii)依据迭加性质及以上两步对单元方程进行合并。4、施加边界条件。5、求解总体方程。6、输出结果§1-4有限元方法的组成模块§1-5基本未知量的选择和有限元方法的分类1、位移单元以位移做为基本未知量，几何关系和弹性关系精确满足，平衡方程只能近似满足(近似解)。优点：未知量少；缺点：位移精度好，应力精度低。位移单元又分为两种：协调单元和非协调单元。位移单元是目前应用最广的一类单元。2、平衡单元以应力或应力函数为基本未知量，以应力协调关系(或最小余能原理)建立有限元方程。3、杂交(混合)单元同时以位移和应力为基本未知量。这种单元在处理板、壳问题时有显著的优点。§1-6结束语现代仃限元方法起源「矩阵结构分析方法，经过科学家和数学家的工作已经使其不断成熟完善，从而变成工程分析中一种处理偏微分方程边值问题的最有效的数值方法之一，对于工程中的许多场变量的定解问题，通过此方法可以得到满足工程要求的近似解。应用此方法在对于连续介质问题分析时，首先要将求解域离散化，然后的中心工作是单元分析(要用插值函数近似表达场变量)，即建立结点位移与结点力之间的关系(形成单元刚度矩阵)。此后的工作可以认为是程式化的工作，即组装总体方程和求解此方程，在解方程时要用到数值积分。变分原理在有限元方法中有重要的作用，在许多场合依据变分原理可以给出总体方程合理性的满意解释。创建新型单元的工作仍在不断进行，人们始终没有放弃构造新型的功能更强、性质更优的单元的努力。参考文献P-G•西阿莱，“有限元素法的数值分析”，蒋尔雄等译，上海科学技术出版(1978)。G•斯特朗，G•J费克斯，“有限元法分析”，崔俊芝等译，科学出版社(1983)。[3]姜礼尚，庞之恒，''有限元方法及其理论基础”。人民教育出版社(1979)。[4]陈传淼，''有限元方法及其提高精度的分析”，湖南科学技术出版社(1982)*[5]欧阳#，马文华等，“弹性塑性有限元”，湖南科学技术出版社(1983)R・DCook,“有限元分析的概念和应用”，程耿东译，科学技术出版社S•S•RAD,"THEFINITEELEMENTMETHODINENGINEERINGS1980)R,Glowinski,E,YRodin,O,CZienkiewicz,“EnergyMethodsinFiniteElementAnalysis”,。978)J,E,AKIN,“ApplicationandImplementationofFiniteElementMethods”(1982)I•霍拉德，K•贝尔，“有限单元法在应力分析中的应用"，凌复华译，国防工业出版社(1978).[11]董平，J•N•罗赛托斯，“有限单元法一一基本方法与实施”，张圣坤等译，国防工业出版社,(1979)G,P,Bazeley,Y•K•Cheung,B•M,Irons,C,C,Zienkiewicz,''Triangularelementsinplatebending conformingandnonconformingsolutions”,Wright----pattersonI(1965)H•G,Schaeffer,“AreviewoftheinternationalSymposiumonstructuralmechanicsSoftware”,ComputersandStructuresVol.8No.51978*[14]J.N.Reddy:ANINTRODUCTIONTOTHEFINITEELEMENTMETHOD.McGraw-HillBookCompany.1984.*[15]王勖成，邵敏：“《有限单元法基本原理和数值分析方法》清华大学出版社，1997[16]谭建国：《使用Ansys6.0进行有限元分析》北京大学出版社，2002.5.*[17]吴鸿庆，任侠：《结构有限元分析》中国铁道出版社，2000*[18]朱伯芳：《有限单元法原理与应用》中国水利水电出版社，1998.10注：* 主要参考书第二章结构矩阵分析由于有限元方法起源于力学中的结构分析，本章的作用是通过三个典型问题说明有限元方法应用于结构分析时的•般步骤，并借此了解有限元方法的一些基本概念。§2-1平面桁架（直接法，结构矩阵分析中常用的力法，处理静定问题，位移法，可处理静定&静不定）本节讨论的时象是图2-1所示的平面桁架。组成桁架的各杆为等截面直杆，外载荷p直接作用于杆的钱接点（结点）上。为简单起见不妨设各杆的截面积均为A,材料的弹性模量均为E。我们可按下述步骤求得桁架的变形和内力。图2图2—1图2—21、结构的离散化对结点及单元编号取组成桁架的每根杆为一个单元（该问题本身为一离散结构的力学问题），以①，②，③加以编号；取杆的较接点为结点，以1、2、3加以编号（总体结点序号）。如图2—2所示，即：我们所讨论的桁架包括三个单元、三个结点。各单元（杆）仅在结点处连接。2、建立总体坐标系并确定结点坐标和自由度为了描述结构的平衡需要建立一个坐标系，称为总体坐标系，以区别于以后出现的“局部坐标系”。总体坐标系的选择原则上不受限制，但希望使用方便。本节所选的总体坐标系示于图2-2,坐标原点与结点1重合。以u,v分别表示沿方向的位移分量，p.q分别表示力沿乂y轴的力分量（投影）。在总体坐标系中各结点的坐标为：（xi,yi）=（0,0（X2,y2）=（a,a）、（x3,y3）=（a,0） xupTOC\o"1-5"\h\z它们将作为程序的输入数据（几何参数）。 /每个结点有两个自由度，对结点1、2、3分别为T T T yVq /1{«/.V/} 、{U2>V2} 、 {W5.Vj}yvq/若暂时不考虑支承约束条件，整个结构的结点自由度 \ 7为 \/\u {tilV1U2v2U3V）}' y1« xup3、单元分析（建立结点力与结点位移之间的关系） *取一个一般性的单元，设它的两个结点在结构中的编号为（单元内部的结点序号）。由材料力学可知，杆 图2.3的轴向刚度为EA/L。其中L为杆的长度：

L= +(yj-yi)（1）单元局部坐标系现选取一典型单元对其进行单元分析，对所分析的单元按如下方式建立一个坐标系：原点：与结点i重合，x轴：沿ij方向,y轴：与x轴垂直。如图2・3所示。这个坐标系只属于一个单元，故称为单元局部坐标系，不同单元的单元局部坐标系一般是不相同的。在单元局部坐标系中可以规定：结点自由度 {〃jN：}'， {UJVy}T；单兀结点自由度 =V,iMj%}丁。（2）局部坐标系中的单元刚度矩阵在外载荷作用下，结构发生变形，单元必受到来自结点的作用力。桁架中的杆只承受轴向力S,大小与杆的轴向伸长△L成正比在局部坐标系中这种特性可以得到清楚的表述（这一点也是引入局部坐标系的理由之一）。若以pt,qi,pj,qj分别表示结点i,/作用于单元的力在x,y轴上的投影，由①号单元的静力平衡有（图2-3）有，_EA.,,、Pj= （"二％） 、*=q：=Q -p[=—s= ~u，*=q：=Q -用矩阵的形式可以写成若引入单元广义力矢量：{/}={p；q\*则上式可缩写为(2-1-2)0(2-1-2)00(2-1-3)0■10-1阳2°°°A-101000

称为局部坐标系中的单元刚度矩阵，它只与杆的几个参数E、A、L有关，与杆的方位无关。(3)坐标变换局部坐标系中的单元刚度矩阵公式简捷。但不同单元的局部坐标系一般不同，为了研究结构整体的平衡，必须将结点给单元的力以及相应的单元刚度矩阵转换到统的坐标系——总体坐标系。在总体坐标系中单元结点自由度 ｛0｝=｛斯V；UjVj｝T结点给单元的力 ｛工｝=｛科QiPj%｝T在图2-3中，x'轴与x轴的夹角为a结点的位移分量的坐标变换为sina=八%>=cosasina>=[4Uj一sinacosaVjY,Lu结点的位移分量的坐标变换为sina=八%>=cosasina>=[4Uj一sinacosaVjY,Luj\/>―cosasina■=Huj一sinacosakJ单元的位移分量的坐标变换为%cosasina0o—%/匕<r»=-sinacosa00匕■=[t｝匕%00cosasinaUjuj—00一sinacosaVJ.vj.(2-1-4)或缩写为怔｝=［施｝类似，｛[｝与位｝之间的转换关系为怔｝=[选｝ (*5)由于[]cosasina—sinacosa (2-1-6)是正交矩阵，因此[T]=(2-1-7)也是正交矩阵。所以有将(2-1-4).(2-1-5)代入(2-1-2)有一比—旭｝从上式可得到.=[T]Tk]T旭｝=炖｝ (2-1-8)其中M=[t]T[^It] (2-1-9)称为单元在总体坐标系中的单元刚度矩阵。以后将会看到，(2-1-9)是一个具有普遍意义的公式。它表明，当单元的自由度由一种形式换成另一种形式时，单元刚度矩阵只需进行一次相似变换。对于平面桁架单元，将(2-1-3)、(2-1-6)、(2-1-7)代入(2-1-9)可得到更便于应用的单元刚度矩阵公式cos2a

cos2a

cosasina

-cos2a

-cosasinacosasina

sin2a

一cosasina

一sin2a2-cosa一cosasinacos2a

cosasina-cosasina

一sin2a

cosasina*2sina(2-1-10)(4)具体结果由(2-1-10)可求得各单元的刚度矩阵的具体形式如下：单元①：单元自由度｛孙V；u2v2｝T,a=45°单元刚度矩阵为[4=-a[4=-aVI4vi4v2-4后一4-723T一4近一4¥¥¥¥V2-4V2-4V24724单元②：单元自由度{〃2旷2 丫3}Ta=-90°单元刚度矩阵为000O'r.1EA010-1阳2=a0000(2-1-12)0-101单元③：单元自由度{UjVju3\屋3}T，a=0°单元刚度矩阵为■10-10r.iEA0000l勺3=a-1010(2-1-13)0000请注意，单元刚度矩阵与单元自由度中位移分量的排列次序有关。如果改动这种排列次序,例如对①号单兀，将单兀自由度次序山｛〃/V/〃2W｝T改为｛〃2也V/｝T»必然导致刚

度矩阵(2-1-11)元素位置的变动。(5)单元刚度矩阵的物理意义和特点设平面桁架单元在总体坐标系中刚度矩阵的一般形式为同=2131k.k.k.131444同=2131k.k.k.131444由(2-1-8),当单元结点位移为｛10)丁时，在单元各结点上施加的力刚好为单元刚度矩阵中的第一列：｛ku心k3.k“｝\对［k］的其他各列也可做出类似的解释。即单元刚度矩阵的每一列相当于一组特定位移下的结点力，如表27所示。由图24可以获得更为直观的理解。表2T平面桁架上台元刚度矩阵的物理意义单元结点位移作用于单元的结点力{1000}T{女〃k21kuk4i}T{0100}T{^12女22女32自2}T{0010}T[如3223%33自3}{0001}T{k[4左24234如}T图2-4对图2-4中的各种情况，据平面力系的平衡条件应有3+*=ok2s+%4s=0(X,fh-(匕-xh=。"(S=l〜4)

这三个关系说明，［k］的四个行向量中只有一个线性独立（四个元素有三个约束方程）。从以上分析可以看出，一般的单元刚度矩阵均具备以下两个特征。（对平面桁架单元而言，从（2-1-10）也可以得出这些结论）（i）单元刚度矩阵是对称矩阵，这是线性系统互易定理的具体体现。由于对称性，对行向量或列向量两者之一得到的结论，对另一个也适用。（ii）单元刚度矩阵是奇异矩阵。它的行向量（或列向量）线性相关，具有零特征值，det［k］=0.对平面桁架的单元刚度矩阵而言，它的四个行向量（或列向量）中只有一个线性独立，而［k］有三个零特征值。这三个零特征值对应的特征向量相当于三种独立的刚体位移模式：两个平移，一个旋转。这是我们在单元分析中不考虑位移约束条件的自然结果。4、总体刚度矩阵的组装 总体平衡方程将图2-1所示的桁架中的支承约束以约束反力代替，如图2-5所示。下面来建立平衡问题的有限元方程。（1）结点平衡条件作用于图2-5每个结点上的外载荷、支座反力以及来自单元的力应处于平衡。以p产、q!m,示结点i作用于单元m的力在x,y以p产、q!m,的力在x,丁轴上的投影应为一p”、一%的。对结点1：对结点1：对结点2：对结点3；可以合并成仇⑴:0图2—50/⑶&YP2⑴►+«Pl(2)0Pq”0夕2⑵P3⑵»TV0P3⑶00(2-1-14)04、.犬3y.式（2-1-14）的右边为外载荷和支反力。左边则为单元给结点的力，它们是未知的，但可以借助单元刚度矩阵以结点位移来表示。（2）单元刚度矩阵的扩充为了表示（2-1-14）左边的各个列向量，设想将每个单元的自由度扩充到与结构总体自由度相同（本例为6）,并在单元刚度矩阵中补充零元素，由（2-1-11）、（2-1-12）,（2-1-13）和（2-1-8）可以用结点位移表示（2-1-14）左边的各列向量。由单元①(2-1-15)ooooooooooooV24vi4vi-4^-4ooV24&4VI-4VI-4Oo收一4a一4vi4VI40oV24V2一4->/2丁-V240on⑴〕P\c⑴qIP2⑴EA>= 夕2⑴a00由单元②(2-1-16)000000'"1"10000004K000000>=[k]2,M,*2200010-1匕匕000000町"3000-10I_V3.EAa小⑵。3⑵S'"由单元③(回+皿+吼乂w,(回+皿+吼乂w,viU2%"3.V3.pF■I000-Io-s<7i<3>000000vlvl0EA000000u2r.iu20►—1，a000000V2>=团3,V2> (2-1-17)(3)P3一100010"3“3(3)[<?3000000V3.(3)组装总体刚度矩阵将(2-1-15)、(2-1-16),(2-P-17)代入(2-1-14)得到&x'&丫

P00.R3y.其中%&x匕段丫"2*=«P>0匕“30&Y.(2-1-18)1+2+3=（2-M9）上式称为没有考虑位移约束条件情况卜的总体刚度矩阵（求和对所有单元进行）。时本节所分析的平面桁架有4-14-1000-叵-亚1n4J54V乙00V,R44y\V200u2<p>=<> (2-1-20)44V20V2V2.——+10-1"30440010I'3J3yJ0-101_要形成总体平衡方程（2-1-18）,只需组装出它的右端项和总体刚度矩阵［©就足够了，这只是同一件事的两种不同的提法而已。（2-1-19）用来“书写”组装总刚度矩阵的过程是简单而明了的。但事实上不能这样做,将所有［k］m补充零元素扩充为［©m极大地浪费了宝贵的存贮空间，这些零元素仅起到使白益的元素在总刚度矩阵中就位的作用。实际上采用的是另一种方法。如果已选定各结点位移在结构总体自由度中的排列次序为{川V,U2V2U3V3}T,对每个单元在形成单元刚度矩阵的同时还形成了个定位数组LM,它将指出单元自由度的各分量在总体自由度中的序号，如下表：单元号单元自由度LM(1)LM(2)LM(3)LM(4)①{UiV1U2V2了1234②T{U2V2U3V3}3456③{uiviU3V3)11256单元刚度矩阵中第s行第t列的元素kst加到总刚度矩阵的第LM（s）行LM⑴列即可。这一组装总体刚度矩阵的方法被形象的称为“对号入座二从（2-1-20）可以看出面是奇异阵，它的六个行向量（或列向量）中只有三个线性独立。这是尚未考虑位移约束条件，结构的刚体位移未受到限制的必然结果。

(4)引入位移约束条件由图2-1,平面桁架的位移约束条件为(2-1-21)U］=V|=V3=0(2-1-21)代入方程(2-1-18),得到EAa子后JJL。7EAa子后JJL。7

斗4¥也4。。V2一4后一4V2-4或-40oooO-on一ooo1io0R\x0u2»=«Pv20“300.R3y.它显然可以分成两个方程组oO1V2ZV板一4V2一40oO1V2ZV板一4V2一40-血4企(2-1-23)(2-1-22)就是平常所说的平衡问题的有限元方程，一般把它写成［kRu｝=｛e｝ (2-1-24)［K］为考虑了约束条件(2-1-21)后的总体刚度矩阵。｛U｝为非约束自由度位移向量，｛F｝为作用于这些自由度上的载荷向量。矩阵［K］相当从矩阵因中划去了第1、2、6行和第1、2、6歹IJ。这说明(2-1-21)这种零位移约束条件，可以通过对总刚度矩阵划行划列的方法来实现二至于其他形式的位移约束条件的实现方法，后面将有专门的一章加以讨论。实际上，在组装总刚度矩阵的过程中，只要不组装应划去的行和列即可直接得到［K］和｛F｝,即方程(2-1-24),至于方程(2-1-23),尽管求得结点位移后可以由它求得支承反力，但是有限元分析中一般不形成这个方程。如确有必要求支座反力，将另找其他途径。

5、解有限元方程对于平衡问题，归结为解一个代数方程。本例的自由度较少，从(2-1-22)可以方便地求得非约束自由度的位移.、(2-\/2+1) 4 EA* aP《V、》《 >被约束自山度根据(2-1-21)可直接赋零6、求单元内力桁架单元的内力只有一个轴力S。由(2-1-1)S的大小和正负与pj相同。由(2-1-2)和(2-1-4)也｝=卜恒｝=[k[T版｝取出它的第三行即得到(2-1-24)cEAr . .］匕(2-1-24)S=—I—cosa-sinacosasina\L Uj3对单元①、②、③求得内力S于下表单元号单元结点位移内力(以拉为正)①(00(272+1)---}TEAEA41P②{(2V2+1)---00}TEAEA一产③{0000}T0由于本节所讨论的桁架是一个十分简单的静定桁架，用理论力学的知识即可得到各杆的内力，结果相同。但是，若桁架为静不定桁架且杆的数目有上千个，那么本节所讨论的方法原则上不会遇到任何困难，我们要解决的课题将转为如何管理有关的大量数据和如何解•个数下阶的代数方程组这样一些技术问题。

§2-2平面框架本节将讨论另一个典型问题——平面框架，示于图2—6。构成它的元件是梁。为讨论简单，不妨设外我荷只作用于梁的刚结点上。我们用同上节类似的步骤求得框架的变形和内力。1、结构的离散化取每根梁为一个单元，编号为①、②、③。取梁之间的刚结点为结点，以1、2、3、4对结点进行编号。2、总体坐标系结点坐标和自由度轴在框架平面内，原点与结点1重合，z轴与框架平面垂直，如图2—6所示。"/分别表示各结点沿x,y方向的位移，。为绕z轴的转角。p,q;分别表示力在x,y轴上的投影，M为绕z轴之力偶矩。不难定出框架四个结点的坐标为（阳，必）=（。，0）。3，%）=（2凡。）（》4，、4）=（2。,0）暂时不考虑支承约束，整个结构的结点总自由度为考虑1、4结点处的约束条件=匕=4=〃4=二°(2-2-1)考虑1、4结点处的约束条件=匕=4=〃4=二°(2-2-1)(2-2-2)V2(2-2-2)V202〃3匕结构的非约束结点自由度为3、单元分析取一个一般性的单元进行分析，设它的两个结点在结构中的编号为i、j»见图2—7。（1）建立单元局部坐标系原点：与结点i重合；x轴：沿i,，方向；V轴：在xy平面内，与x轴垂直;Z轴：与Z轴一致。单元局部坐标系中，各结点自由度为单元结点自由度可表示为[={〃：/4UjM打(2)局部坐标系中的单元刚度矩阵根据单元刚度矩阵的物理意义，它的元素相当单元结点位移在六种位移模式{10000o}T-{o00001}T情况下的结点力。如果我们约定：(i)单元内弯曲刚度EI为常数；(ii)单元只在端点受到外力；(iii)不计剪切变形。则”将是x的线性函数，v是x的三次函数，8=骼。利用材料力学知识不难求得梁的变形曲线和结点力，如表2—2所示。设单元长度为L,截面积为A。图中实线为现实构形，虚线为参照构形。

由此求得平面梁单元在局部坐标系中的单元刚度矩阵［R■EA~L00EA~L00\2EI6EI\2EI6EI001}I?I?片6E14EI6EI2EI001=EAL2LEAL2L(2-2-3)0000LL12EI6EI\2E16E100I?I?I?L-6E12EI6EI4£70ITL0[2L.它是一个对称阵；也是一个奇阵。在单元分析中没有考虑位移约束条件，允许单元作刚体运动。我们还可以看到，局部坐标系为单元分析提供了方便。(3)坐标变换为了讨论结构的平衡，需要把局部坐标系中形成的单元刚度矩阵转换到总体坐标系。若图2—7中x轴与x的夹角为a,则两个坐标系中位移分量(以结点i为例)的转换关系为cosasina0‘外一-sinacosa0匕•=H%(2-2-4)0 0 1同同每个单元有两个结点，单元结点自由度的转换关系为W｝=j［也｝=T应｝ (2-2-5)利用与2-1节中类似的步骤，可推出单元刚度矩阵的变换公式为同向卜什］ (2-2-6)形式上与(2・1・9)完全相同。但请注意，矩阵口］的定义不同，因而［T］的具体表达式也不同。

利用(2-2-3)和(2-2-6)可求得单元①、②、③在总体坐标系中的单元刚度矩阵［kh、［k］2>［k］3o4、总体刚度矩阵和载荷向量的组装本节侧重讨论在组装刚度矩阵的过程中实现划行划列，即由因直接组装［K］的方法。对每个单元，根据(2-2-1)和(2-2-2)可形成一个数组LM,如下表所示：单元号单元自由度LM(I)LM(2)LM(3)LM(4)LM(5)LM(6)①{"1V,d"2 %000123②{"2V202 "3匕/}'123456③{%匕“4V4456000每个LM数组有六个元素，分别与单元自由度中的各位移分量对应。当这个位移分量根据(2-2-1)被约束时(例如»/,v/(0,.u4,v4,%)对应的LM元素取零；若这个位移属于非约束自由度，对应的LM元素取该位移在(2-2-1)中的序号，利用LM数组即可实现由［k］直接组装［K］。对于［k］的第s行第t列的元素却若LM®和LM(t)中至少一个为零，则对该元素不予组装(相当于划行划列)；否则将右迭加入总刚度矩阵［K］的第LM(s)行LM⑴第列。非约束自由度中只有“2的自由度上作用外载荷P，故总体载荷向量为/}={p0000o}T最后得到有限元方程[kRu}={丹最后得到有限元方程[kRu}={丹(2-2-7)图2图2—85、有限元方程(2-2-7)得到非约束自由度结点位移，被约束自由度上的位移，据(2-2-1)直接赋零。6、由结点位移求单元内力。例如：轴力、剪力、弯矩，对空间梁单元还可求扭矩。§2-3平面应力问题常应变三角形图2-8为一边长为a、厚度为t的正方形薄板。其中AB边固定，BC、CD边自由，AD边作用均布压力外对这一问题，有限元分析的步骤是：1、将ABCD划分(离散)为8个三角形(单元)，编号①1⑧。各单元仅在顶点(结点)较接，结点编号1—9。建立坐标系后，不难定出各结点的坐标2、单元分析任取一个一般性的单元，如图2-9所示。三个结点的编号为i.j,鼠结点位移为(Ub片)、(Uj,Vj),(uk,vk)单元结点位移为

8}={%匕"jVjukvJT(1)假定单元内位移场",V是x,y的一次函数〃=«+%x+a3yv=cr4〃=«+%x+a3yv=cr4+cx5x+a6y(2-3-1)名〜&6为待定常数，在结点处应有可解出:,1y, , 1y>2A=2A=detDA=a,+aj+ak当i,_/,左的位置为逆时针排列时，2△恒正，且等于三角形单元面积的两倍。将这些结果代入(2-3-1)有〃=白(4+3+C/H+20+与x+c/)。+白(4+bkx+cky)uk=Ni(x,y)ui+N/x,y)“j+Nk(x,y)uk

类似可得到v=Ni(x,y)vi+Nj(x,y)vj+Nk(x,y)vk可以合并成o0*lvd0NjN：00Nk*类似可得到v=Ni(x,y)vi+Nj(x,y)vj+Nk(x,y)vk可以合并成o0*lvd0NjN：00Nk* 0也｝(2-3-2)(2)单元的应变、应力利用(2-3-1)不难求得oa¥Aaxa-dxo}.也.- _oMsoON，5ooMMo_一oa苏a-ax其中bj0bj0c1.0J瓦Cj040c,0ckJ nbj。卜bk{«}=历应}(2-3-3)J 0 c*号 c* bk(2-3-4)在假定单元内位移场是xj的一次函数的前提下，单元内的应变和应力将是常数,故这种单元又称为常应变三角元。

(2-3-4)图(2—10)(3)为了在单元内构成均匀应力场，必须在单元的各边施加均布载荷，它们的合力一定作用在各边的中点，如图2-10(a)所示。再将各边上的合力平分到这边的两点结点，由图2-10(b)不难得出Pi=；(凡+心)=:配•一匕4+&-Xj)rxy-(yk-匕-(x,-xk)rxy]=；阮-H沅+&-xj卜孙12 -类似可求得%、0、%、Pk、qk并可合并写成Pi瓦0ci%0Cib.Pj鸟0ci%20CjbJPk4°q4.0ckbk根据单元刚度矩阵闵的直观意义，常应变三角元的单元刚度矩阵即为同=时团团(2-3-5)(4)单元①和②的边上作用着均布的外载荷，可以把它们的合力平分到两结点，如图2-11所示。f\= -X7) Zz=54(X1一工4)

图（2—11）(5)为了组装总体刚度矩阵，每个单元还应形成一个数组LM。元素LM(1)〜LM(6)分别为Ui、Uj、V,'姝、%在总体自由度中的序号。由于单元的结点位移和结点力都是用总体坐标描述的，单元刚度矩阵不必再进行坐标变换。3、组装总体刚度矩阵和载荷向量。为了实现位移约束条件〃7='==丫9=0这些自由度对应的行和列可以不必组装。总体平衡方程为国］｛。｝=｛尸｝ (2-3-6)其中,非约束自山度位移为｛。｝=｛"1匕〃2丫2"3匕〃4V4U5V5UbV6｝T相应自由度的载荷向量为伊｝=｛0-/200000-(/,+/2)0000｝T解方程(2-3-6)可得到各非约束自由度的位移，再山(2-3-3)可求得各单元的应力。§2-4结束语1、本章所讨论的三具体问题尽管力学背景不同，但分析步骤大体相同。区别仅在于结点参数不同，单元分析公式不同。而组装总体矩阵的方法、处理约束条件的方法以及代数方程组的求解方法完全通用。2、从建立有限元方程的方法看，本章属于直接方法一类。这种方法只在简单情况下有效，但是物理意义清晰。尽管今天有限元方法已发展成为一种抽象的数学分析方法，本章提供的直观解释仍然常常被用来说明一些问题。3、对于第一、二两个问题(平面桁架和平面框架)本章求的将是精确解(真实解)。而对第三个问题(平面应力问题)，本章求得的只是近似解(位移场是假设的)，这样就自然引出以下两个问题：(1)所求的近似解与真实解的接近程度如何？或者说，怎样才能接近得更好一些？(2)在进行单元分析时，曾作了两条关键性的假设：(i)单元内位移u、v是坐标的一次函数：(ii)各单元仅在结点处较接。这样的假定有没有依据？是否必需？这些问题将在后续的章节中讨论。第三章最小势能原理和分片插值最小势能原理和分片插值是有限单元方法的核心内容之一,在这一章中我们将介绍最小势能原理的具体应用，同一力学问题的几种不同的表达方式及它们之间的联系；Ritz方法在单元内的应用；两种最常用的插值形式(Lagrange型和Hermite型)。协调的位移型单元的收敛条件。§3-1最小势能原理一个平衡问题，可以至少用以下叁种不同的方式加以描述:(i)平衡方程(ii)虚位移原理(iii)总势能取驻值(函数的极值问题)1、有限自由度系统质点系图3-1(a)为两个重分别为已，Pe的小球，由不计重量，弹性系数为％的弹簧相连，放置在光滑的曲面尸优切=0上。该系统的平衡问题可由以下三种方法来描述：(1)平衡方程(在平衡位置描述平衡)分别取月、6为研究对象，建立一直角坐标系，如图3-1(b)所示。由？1、3两球组成的系统的平衡条件为：A：XX=0B：ZX=。Zy=0Xr=0其中，弹簧力以压力为正，即：T=k(l0-AB)(2)虚位移原理(在平衡位置附近描述平衡)系统的平衡位置可描述为以,rB或XA>yA,kb、ys如图3-1(c)所示。两质点的虚位移分别为：8rB或8xa,6%,6xb,8yB由虚位移原理可得pAa+QbM-亍3(AB)=0上式也称为微分型变分原理。进一步改写为-。血-P血+他-45»(/3)=0一)PM+PM+5H43-/。丁=0(3)总势能取驻值(在全范围内描述平衡(寻找平衡位置))引入总势能函数万p=gk(AB-lJ2+pm+PbYb由总势能函数取驻值，即：6兀p=0寻找平衡位置。以上三种描述系统平衡的方法，对于有限自由度系统而言是等价的。然而,三种方法的特殊性是：平衡方程是在平衡位置上描述平衡问题；虚位移原理是在平衡位置附近描述平衡问题；而最小势能原理是在全局范围上描述(寻找)平衡问题。最小势能原理可以表述为：在所有满足给定位移边界条件和协调条件的位移中，满足平衡条件的位移使总势能取驻值，若驻值是最小值，则平衡是稳定的。作为•条基本原理，它的正确性应由事实加以检验，不必从理论上给以证明。但为了便于理解“位移边界条件“和“协调条件”的含义，下面通过三个例子对最小势能原理和平衡方程的等价关系加以验证。2、无限自由度系统弹性体(1)轴向受拉的直杆设杆长为入，截面积为A,轴向分布载荷f(x)^x=0端固定，户L端受端点集中力P»设位移"(X)满足：(i)〃仰=0(位移边界条件)(ii)“㈤在上连续(协调条件)

(iii)调条件)

(iii)使取最小值。[L L7tp(w)=—公一j/udx-Pu(L)o o若u(x)+巾〃㈤为不同于u(x)的另外一种位移分布函数，也满足上述的位移边界条件和协调条件，山[w(x)4-血(x)[=o=w(0)+质«0)血(0)=0血(0)=0(3-1-2)不难求得：L L .L兀 +3u)_兀p(u)=^EAu8nfdx-^f・企欣-尸加(£)+-j£4(式/「治o o 2oL=加〃+-,/(血')2公o因np(u)取最小值，即兀p(w+8u)—兀p(w)>0

的充分必要条件是：对任意满足(3-1-2)的5〃㈤有加「=\EAu'dii'dx-J/-8udx-Pdu{L)=0 (3-1-3)0 0式(3-1-3)即势能驻值条件。若补充假定“'㈤存在、连续，则对(3-1-3)分部积分一次，并利用(3-1-2),可得到+fdudx++fdudx+\EAu'(L)-P]8n(L)=0(3-1-4)(3-1-4)式对任意Su(x)都成立的充分必要条件是：/(E4/)+/=0 (平衡方程)EAu'(L)-P=0 (力边界条件)由上述过程不难看出：①由势能取驻值可以推出平衡方程。反之也对，说明两种描述方法在力学上等价。②两种描述方法对边界条件的要求不同。用微分方程描述时，〃必须满足：(£1//)+f—0 (平衡方程)u(0)=0 (位移边界条件)EAu\L)-P=O (力边界条件)用最小势能原理描述时，要求函数满足位移边界条件“(0)=0而力边界条件将作为势能取驻值的自然结果。当"④使"P精确取驻值时，平衡方程和力边界条件将精确地得到满足；当"㈤使"P近似取驻值时，平衡方程和力边界条件只能近似的得到满足。在这种意义下，力边界条件又称为自然边界条件(NaturalBoundaryCondition),而位移边界条件由称为强制边界条件(EssentialBoundaryCondition)»③两种描述方法对函数的光滑程度(即可微性)要求不同。用微分方程描述时要求“同有连续的二阶导数(记作"GC2(0,L))。而用最小势能原理描述时，为了保证变形能存在,要求"々平方可积(记作"GH'(0,L))(2)平面应力问题%("，v)=；JJ{"lE胆}以Mv」正方形区域边长为a,厚度为t,受到体积力区工),边界48固定。边界BC%("，v)=；JJ{"lE胆}以Mv」B o cn 图3—31H1为一阶导数平方可积的函数组成的函数空间，C?为二阶导数连续的函数组成的函数空间。

UIAB=VIAB=O利用唱｛叫础£｝卜|候尸团｛£｝+｛卬历倏｝)唱｛叫础£｝卜|候尸团｛£｝+｛卬历倏｝)=｛时展｝=㈤小f»<»k

oAayAax

AaxoAay以及格林公式(其中Lx，5为区域Q之边界广的外法线〃的方向余弦)。势能的驻值条件为:Adxoa¥oAblaxAdxoa¥oAblax(3-1-5)注意到沿边界r,外法线〃的方向余弦为ABBCCDDALx-10I00-101以及沿线ZB：8u=8v=0oa-砂a-dxAaxoa-oa-砂a-dxAaxoa-力沿BC：。.尸t沿CD：o=tx=O沿AD：。，q,txFP关于边界条件的讨论以及从v可微性的讨论与上例相似，这里从略。(3)梁的平面弯曲O对于图34所示的弯曲问题，总势能和强制边界条件为OL-jq-vdx-Q-v(Z)+M•/(£)ov(O)=v，(O)=O势能驻值条件b兀p= •加必・凉(L)+A/•苏'(£)=0 (3・1»6)对上式分部积分两次，并注意到由于必须满足强制边界条件5『㈤=3/如=0则有6兀p=J(£7/)-q3vdx+[£7v*(£)4-M^v\L)。|_」一(E")' +Q(5v(L)=0 (3-1-7)_ x=L_使(3-1-7)对任意Sv(x)都成立的充分必要条件是：〃(EIv")-q=O(平衡方程)EIv\L)+A/=0 ]⑻行+Q=OJ(自然边界条件)x-L微分方程的阶数为4。关于一、J的边界条件为自然边界条件，关于V.v的边界条件为强制边界条件。当用微分方程描述时要求3市四阶的连续导数IvGC,QL)]。用最小势能原理描述时，为保证变形能存在，只要求「㈤平方可积[vGH2(0,L)|。本节讨论的三个例题，可作为维数不同，阶数不同的问题的代表。现把一些重要结论归纳如下表。其中，“微分方程阶数”是以位移为基本未知量来计算，“可微性要求”是对最小势能原理而言的。问题微分方程阶数强制边界条件协调条件可微性要求杆的拉伸2关于U的边界条件〃连续〃平方可积平面问题2关于U,V的边界条件U,V连续U,V平方可积梁的弯曲4关于V.V的边界条件V,V连续V,平方可积§3-2Ritz法(有限元方法的基础之一)由于有限单元方法可以理解为在单元(子域)内应用的Ritz法。Ritz法是一种求近似解的常用方法，它的基本步骤是：(i)选一组满足强制边界条件、协调条件和可微性要求的函数⑺、夕？、……中”这组函数被称为基函数。(ii)假定近似解(试探函数trialfimction)的形式为 +-小(iii)将试探函数作为近似解代入描述问题的能量泛函中，由泛函取驻值，即圾=0,鬼=0…其=0阳da2'四定出系数四〜时。从而得到近似解。下面以图3-5所示的简支梁为例，求解在集中力尸作用下的变形，分别采用不同的基函数来介绍Ritz法。 *解法1：取多项式形式的基函数： 次(px(x)=x(L-x)*2(X)=X2(L—X)显然，"八夕/、。2、92连续，<Pl、(P2平方可积，且(P\(0)=(P\⑷=%(0)=%⑷=o取v=a1^1(x)+a2^2(x)为试探函数,泛函可写成、 Bx7CL >图3—56，名为待定常数，则描述这一问题的能量兀p=1\EI(yydx+Pv(^)=2EIL(a^由能量泛函的驻值条件 、*0,da,可得到一代数方程组4%+2La?—2%+4Lc(2—解代数方程组可得试探函数中的待定常数7PL(X,一 ,128EIC点位移的近似值为v(^)=-0.009523 3+cXy£+a~2£~)+ a2'应=0da23PL16EI3PL64EI1Pa)一■64EIPI}~El解法2：取正弦函数为基函数（px（x）=sin——, %（x）=sin—^一这样选择的基函数同样满足协调性、微性以及强制边界豪件的要求。可计算出:42a\42a\+=-EIL4+^-axP-a1-P也=0, "=。阴oa2可求得axV2PI]1PI}%可求得axV2PI]1PI}%=——7 2 8万4EIc点位移兀4EIPT}=-0.011549—EI该问题的精确解为3PI? PI? 3PI? PI? =-0.011719, 256EIEI两种方法求得的c点位移绝对值小于精确值。正弦的基函数，使支座处弯矩为零的条件（不属于强制边界条件）也得到满足。尽管Ritz法本身并不要求这一点，但是第二个近似解的精度显然比第一个要好得多。不难看到，当选取多项式或三角函数为基函数时,由于这些函数可以求导任意次，协调条件和可微性要求很容易得到满足，我们要做的事仅为设法满足强制性边界条件。这在一维情况下很容易做到，但在二维、三维情况下，除一些形状规则的区域（例如矩形、园、扇形、长方体等）外却很难实现。因而古典的Ritz法的应用受到了限制。§3-3分片插值形式的基函数和试探函数（解的收敛性与插值函数的选取关系很大）一维Lagrange型插值图3-6为一轴向受拉的直杆，截面积A和轴向分布载荷/可以是x的函数。因而轴向位移“㈤可能是x的复杂函数。将区间回〃分成若干（这里是三个）子区间（单元），编号为①〜③。取端点和分点为结点，编号为1〜4。坐标为X/〜X/.(1)基函数定义基函数G(1)基函数定义基函数G〜卬4。<Pi(Xj)=8t满足1当户i0当尸了设基函数在单元内是X的一次函数。试探函数的形式取为基函数的线性组4合1=1根据0,的定义显然有：〃㈤是X的分段线性函数，且“(xj=»,系数"，恰好代表结点i的位移值，相互之间是独立的。这样分段(片)定义的试探函数的一个显著优点是：强制边界条件很容易得到满足。例如u(0)=0的条件只要简单地令结点1的位移"产0即可以实现。而且允许我们在任何方便的时候(例如组装总体刚度矩阵时)引入这些边界条件。由于强制边界条件问题已经有了妥善的解决办法，我们的注意力将转向协调条件和可微性问题。(3)协调性和可微性上面定义的如㈤和〃④都存在着“尖点“，光滑程度不高。但是：0,和〃㈤在单元内连续，在结点处也连续；打和“'㈤在单元内连续，在结点处可能不连续。但只有有限的跳跃量。在区间［0,L］上平方可积。内㈤和“㈤属于同一类型的函数。对于轴向受拉杆(二阶问题),“CW满足最小势能原理对协调性和可微性的要求。由可求得A的=o,加3U4的值,0,从而得到一个近似解。(4)Lagrange插值可求得A的=o,加3U4的值,0,从而得到一个近似解。(4)Lagrange插值如㈤、〃㈤都涉及这样一个问题：由两个结点上的函数值在单元内确定一个线性变化的函数。图3-7为一个一般性的单元，两个结点八)的坐标为X”Xj,假定单元内"㈤是X的线性函数u(x)=%+a2x由u{xj)=M(,U(X)=u,定出a\>a2,可得到x—x

w(x)= -ui+x,-X,x-xi=+Nj(x)uj其中Ni、2称为形函数，它们在单元内是x的线性函数，且满足n.=0M(x,)=lNG/)=ON/(xJ=O Nj(x)=\每个形函数由分子和分母两部分组成，分子保证了一个结点的形函数在其他结点处为0,而分母的选择则恰好使得这个形函数在自己的结点个取值为1。掌握了这些特点就可用“凑”的方法得到形函数的表达式。M(x)=[x-X-\x-X/)M(x)=N,(x)=(X—X/\X—Xj)(x-xjx-xj如果在每个单元内在增设一个结点M(x)=[x-X-\x-X/)M(x)=N,(x)=(X—X/\X—Xj)(x-xjx-xj试探函数m（x）=NjUj+NjUj+N/U,图3-8是它们的图形。用插值点的函数值构造的插值函数通常称为Lagrange插值。一维Hermite型插值图3-9为一根梁，横向载荷q和截面惯性矩I可以是X的函数。因而挠度V是X的复杂函数。梁的弯曲是四阶问题，试探函数V及V应在［0,〃上连续。将［0,〃分为若干（这里仍是三个）子区间（单元），编号为①〜③。取端点和分点为结点，编号为1〜4。(1)基函数定义基函数处60〜%。入M(x)〜«4(X)满足(i)(p/x)、M 在单元内是x的三次函数(ii)r1当j=i/(X/)=%={ d(X/)=010当jW/C1当j=i，(xj)=B={ %G)=o0当产了(2)试探函数 4 4/=1/=|根据3，、。，的定义可知，v㈤是x的分段三次函数，且满足v(xi)=V,./(x,)=V；系数匕、V,恰为结点处V、V之值。这些值相互之间是独立的。这样定义的试探函数保证了V、J在单元内连续，在结点处也连续，满足四阶问题所要求的协调条件。由于试探函数是分段定义的强制边界条件v(0)=v(0)=0可以简单地令V!(0)=Vl(0)=0得以实现。(3)可微性W，(x)、”,(x)和v㈤是分段定义的三次函数，函数本身及一阶导数在［0,L］上连续。二阶导数是分段线性函数，在结点处不连续，但只有有限跳跃量，在［0,L］上平方可积。满足最小势能原理对试探函数可微性的要求。如图3-9(续)所示。图3-9(续)(4)Hermite插值定义3i(X)、都涉及这样•个问题：由函数及一阶导数在结点处的值确定单元内的一个三次函数。图3-10为一个一般性单元。它的长度为L,结点为八八i为局部坐标原点，孑与x平行。设梁的挠度函数(试探函数)为v(J)=%+a古+a3g之+a/3

-=a2+2a34+3a/2以结点坐标£=。、/尸/代入则有匕=«!v\—(x2vi=a]+a2L+a3L2+aAI?v'j=a2+2a3L+3a4L2定出4个常数代入挠度函数并整理成/=1/=!的形式，即啥)=。^亿+2a.匕L+延727+。3八2加LL।~j)/+F—，匕=A；(x)v,+/7,(x)v.+%(x)%+H/(x)v'j兀、珏、hj、Hj的图形已示于图3-10。显然基函数的非零部分正是由它们拼成的。用插值点的函数值及导数值构造的插值函数通常称为hermite插值。§3-4常应变三角元的理论依据现在重新考虑图2-8所示的平面应力问题。仍把所研究的区域Q分为8个单元。单元和节编号方案不变。1、基函数定义基函数/i(x,y)〜<p9(K,y)满足 1当j=i(i)0，例,力卜^小0当产i(ii)0,〃,力在每个单元内是x、y的线性函数。2、试探函数“(x,y)=工""v(x,y)=Z/。，图2-8根据0/为,力，的定义，在结点上有

图2-8u(xi,yiu(xi,yi)^uiv(x,，N)=v：在每个单元内，小v是x、y的完全一次多项式，可以由三个结点上的函数值唯一确定,可由第二章的（2-3-2）表示。显然在单兀内心v连续，且其导数dudu 3v'dy、Hx、Hy为常数。沿单元边界（例如3/边），〃、v按线性变化，完全由这条边上两个结点上的函数值小、片、砂、竹所决定，故穿过单元边界时"'V连续（如仙图3-11所示）。但穿过单元边界时其导数、、dudu9v3vdx'dy'dx'dy却一般不连续，有有限的跳跃量，但在Q内它们平方可积。所选择的试探函数满足平面应力问题（二阶问题）对协调性和可微性的要求。可用Ritz法求近似解。k3、总势能将假定的位移场写成u\viu2.J% o 内 o % o （p9 0］匕，=［①Ru}［vj0 夕। 0 % 0 必 0 %：%丫9,万尸=；JJW团和}以叨-1

2Q AD"仿蚓外力-J“.Ni=\.. An"i=\eiADI:

0'

r.tdx>tdx万尸=；JJW团和}以叨-1

2Q AD"仿蚓外力-J“.Ni=\.. An"i=\eiADI:

0'

r.tdx>tdxoAayAayAaxoAax回oAayAayAaxoAax回一画一oa-aya苏Aaxoa一axDe8Z/=1Is,yrHT--[<1>fjw=抠｝］£叫以一邮型［「J/=匏｝「依应｝-的但｝ —I)4、势能取驻值W［kM-W｛£｝=0IK也｝=因强制边界条件/7=匕="8=%="9=丫9=°可以通过对区］、［£］划行划列或在组装［K］,［F］的过程中加以实现。5、单元刚度矩阵［①］在一个单元内有12列全零，故［kh有12行、12列全零。［©i就是第二章中“补零扩充"后的单元刚度矩阵。(3-4-1)的推导过程给出了由单元刚度矩阵IKh组装总体刚度矩阵的另一种解释：总变形能等于单元变形能之和。去掉［©i中全零的12行和12歹U，可得到一个6X6的方阵［k］。若单元三个结点的编号为3人k,当单元内位移按线性分布时,由(2-3-3)有匕回=<£［=［碎。卜历应｝£匕I1以其中［B］的表达式为(2-3-4)o在单元内为常数。单元变形能%=；0缶/历上》公办=fjl5K［EfB｝dxdy｛u\e e=*｝「(时间叫她｝=；也｝也｝

单元刚度矩阵[k]=[BY[E\B^与（2-3-5）完全相同6.等效结点载荷载荷q的等效结点力归结为计算下列两个积分单元①有限元方法可以看成采用分片插值形式的Ritz法。由于试探函数采用分片插值形式。即使在区域形状比较复杂的情况下，强制边界条件也很容易得到满足。但所选择的试探函数必须满足协调性和可微性要求。这是最早出现的关于有限元方法的理论论证。§3-5收敛条件一般说来，用分片插值形式定义的试探函数很难做到与问题本身的真实解（精确解）完全吻合，因而有限元解一般都是近似解。我们希望在网格逐步加密、单元尺度无限变小时有限元解能收敛到真实解。为了保证收敛性，各单元内假定的位移场（试探函数）应满足以下条件（1）假定的位移场在单元内连续。这是最小势能原理所要求的，也是很容易做到的。（2）能够描述任何一种常应变状态（常曲率）。归纳出这一条件的理由是：当单元无限细分时，真实解在每个单元内都将接近某种常应变状态。如果有限元解能够无限接近总势能的真实最小值，随之也就逼近了真实解。满足这一条件的另一个好处是，如果真实解在整个区域上为•常应变状态或者应力变化比较平缓，那么大尺度的单元即给出精确解或较好的近似解。（3）包括足够的刚体位移模式。在假定单元位移场时不考虑位移约束条件（强制边界条件），这些条件在组装总体矩阵时再加以考虑。如果单元位移场缺少某种刚体位移模式，实际意味着加上了额外的约束反力。实现条件（2）和（3）是不困难的。对于轴向受拉杆，只要包含完全一次多项式"=«+ot-^x —=a,dx即可。其中明为刚体平移，为常应变项。对于平面应力问题，只需包含x、y的完全一次多项式即可，即1(a,+a.

y\+\a2x+-yx+abya.-a,u-a}+a1(a,+a.

y\+\a2x+-yx+abyv=a4+a5x+aby-

其中第一个括号内为三个刚体型位移：二个平移一个旋转。第二个括号内为常应变项;加 due——=0C-y. £=—=（X, -+-=a?+a”xdx2vdy6