有限元分析基础

武汉理工大学教师备课专用纸

有限元分析基础 第一章 有限元法概述

在机械设计中，人们常常运用材料力学、结构力学等理论知识分析机械零构件的强度、刚度和稳定性问题。但对一些复杂的零构件，这种分析常常就必须对其受力状态和边界条件进行简化。否则力学分析将无法进行。但这种简化的处理常常导致计算结果与实际相差甚远，有时甚至失去了分析的意义。所以过去设计经验和类比占有较大比重。因为这个原因，人们也常常在设计中选择较大的安全系数。如此也就造成所设计的机械结构整体尺寸和重量偏大，而局部薄弱环节强度和刚度又不足的设计缺陷。

近年来，数值计算机在工程分析上的成功运用，产生了一门全新、高效的工程计算分析学科——有限元分析方法。该方法彻底改变了传统工程分析中的做法。使计算精度和计算领域大大改善。

§1.1 有限元方法的发展历史、现状和将来

一，历史

有限元法的起源应追溯到上世纪40年代（20世纪40年代）。1943年R.Courant从数学的角度提出了有限元法的基本观点。50年代中期在对飞机结构的分析中，诞生了结构分析的矩阵方法。1960年R.W.Clough在分析弹性力学平面问题时引入了“Finite Element Method”这一术语，从而标志着有限元法的思想在力学分析中的广泛推广。

60、70年代计算机技术的发展，极大地促进了有限元法的发展。具体表现在： 1）由弹性力学的平面问题扩展到空间、板壳问题。 2）由静力平衡问题——稳定性和动力学分析问题。 3）由弹性问题——弹塑性、粘弹性等问题。 二，现状

现在有限元分析法的应用领域已经由开始时的固体力学，扩展到流体力学、传热学和电磁力学等多个传统的领域。已经形成了一种非常成熟的数值分析计算方法。大型的商业化有限元分析软件也是层出不穷，如：

SAP系列的代表SAP2000（Structure Analysis Program） 美国安世软件公司的ANSYS大型综合有限元分析软件 美国航天航空局的NASTRAN系列软件

除此以外，还有MASTER、ALGO、ABIQUES、ADINA、COSMOS等。 三，将来

有限元的发展方向最终将和CAD的发展相结合。运用“四个化”可以概括其今后的发展趋势。那就是：可视化、集成化、自动化和网络化。

§1.2 有限元法的特点

机械零构件的受力分析方法总体说来分为解析法和数值法两大类。如大家学过的材料力学、结构力学等就是经典的解析力学分析方法。在这些解析力学方法中，弹性力学的分析方法在数学理论上是最为严谨的一种分析方法。

其解题思路是：从静力、几何和物理三个方面综合考虑，建立描述弹性体的平衡、应力、应变和位移三者之间的微分方程，然后考虑边界条件，从而求出微分方程的解析解。其最大的有点就是，严密精确。缺点就是微分方程的求解困难，很多情况下，无法求解。

数值方法是一种近似的计算方法。具体又分为“有限差分法”和“有限元法”。 “有限差分法”是将得到的微分方程离散成近似的差分方程。通过对一系列离散的差分

武汉理工大学教务处制 1

武汉理工大学教师备课专用纸

方程求解，得到最终的力学问题近似解。其优点就是：计算简单收敛性好。缺点是：计算程序无法标准化，在不能获得整个问题的微分方程时，该方法不能运用。由于其是将微分方程转为差分方程，所以它是一种数学近似。

“有限元法”的基本思想就是“先分后合”或者“化整为零，又积零为整”。与有限差分不同，它是在力学模型上进行近似处理，也就是（分块近似）。

具体做法：把连续体模型转为由有限个单元组成的离散体模型，离散体模型之间通过一些节点联系。对于每一个离散体个体选择简单的函数近似表示其中的物理变化规律（如位移等），运用力学方程推导单元的平衡方程组，然后集合所有的方程组形成表征整体结构的方程组，引入边界条件，求取最后问题的解。

优点：概念清晰、易于学习理解，适用性强，便于电算化。 缺点：计算精度受单元划分的影响较大。

§1.3 有限元分析的一般过程

为了能够了解有限元分析的全貌，我们就一个简单的例子，来分析一下有限元分析的三个过程：结构离散化、单元分析、整体分析。 一，结构离散化

在该阶段中，要完成把连续结构的力学模型转变为离散的力学模型。处理的好坏，直接影响到最后分析结果的正确与否、计算的精度和计算的效率。

根据模型的传力特性和分析的目标，正确选择单元类型。通常单元分为：一维单元、二维单元和三维单元。

所谓一维单元就是指所求物理量仅随一个坐标变量而变化的单元。如桁架、平面刚架和空间刚架单元。

一维单元：杆单元、梁单元。

二维单元：三角形单元、四边形单元（平面类问题） 三维单元：四面体单元、六面体单元等（空间问题）

计算精度和计算效率：取决于单元划分的形状、大小和分布状况。通常单元愈多、愈密集，计算精度愈高，但计算效率愈低。有限元分析工作就是要在精度和效率两者之间做到有机的统一。 二，单元分析

进行单元分析的目的是为了到处表征单元力学特性的“单元刚度矩阵”。一般说来该过程有三种方法：

1， 直接法。

2， 虚功原理法（变分法）。 3， 加权余数法。

直接法概念浅显，易于理解物理含义。

变分法需要泛函的数学知识，其推导过程具有严谨的数学概念。 加权余数法适用于泛函不存在的应用范围。

本教材将运用虚功原理方程结合弹性力学和材料力学中的知识来推求几种常见单元的单刚计算公式。

现在先看一个简单的阶梯轴的轴向拉伸问题

例：如图所示的变截面直杆，受拉力P，运用有限元方法分析其变形。

武汉理工大学教务处制 2

武汉理工大学教师备课专用纸 取任意单元，长度为l，面积为A， 单元内任意一点的轴向位移和其位置坐标成正比，即 u=a0+a1x 其中a0、a1为待定系数。 由于杆的两个端点节点1、2是单元上的点，所以它们应该满足上述方程。 节点1，x1=0，∴u1=a0+a1×0=a0

节点2，x2=l，∴u2=a0+a1×l a1=(u2-u1)/l 将求出的结果带入方程并整理，就得：

u??u?u?u1??2?1?x??1?l???l??N????ex?xu?u2??N1?1l?l?u?N2??1??u2?

式中：N1、N2是形函数 [N]形函数矩阵

{δ}e节点位移向量 由位移与应变的关系知道：

??u?du?udu ?dxdx将上面推出的位移表达式代入，可得：

??u?du?ududd?xe???N?????1??dxdxdxdx?lex?e?1l?e???????????l? 上式中?ll???B????的[B]称为应变矩阵或几何矩阵。

运用材料力学中的虎克定律，可以将应变和应力联系起来。单向应力状态的虎克定律为

??E??E?B????e????l?EE??u1?e?????S? ×× ???ul??2?其中[S]称为应力矩阵。

利用虚功方程可以建立力与位移之间的关系，也就是单刚方程。在后面我们将会推导出它的一般形式如下：

?F?e??K?e???e

式中：{F}e为单元节点力向量，对我们这个例子应为[U1 U2]T。

武汉理工大学教务处制 3

武汉理工大学教师备课专用纸

[K]e为单元刚度矩阵。后面将推导出它的计算公式为?K??eT??B??D??B?dV v[D]矩阵是弹性矩阵。对于一维单元来说，就是E。

所以我们这儿讨论的例题：

?K?e?1???l?EA?1?1??11? ???l?E??Adx????01l??11????ll??l?求得单刚矩阵，也就完成了单元分析。 总结单刚矩阵推导的步骤，应该分为四步：

1） 假定单元内位移变化的近似规律，即选择位移模式。 2） 运用几何关系，推求位移与应变的关系。 3） 应用物理规律，把应变与应力联系起来。

4） 运用虚功方程的力与位移关系，求出单刚矩阵。 单元分析是整个有限元分析的核心。不同的单元因为其力学特性不同，而具有不同的单元刚度矩阵，我们这本教材就是要学习几种常用单元矩阵的推导和计算。了解各种单元的力学特性，为以后选择单元类型打好基础。 三，整体分析

1， 由各单元刚度矩阵组集成整个结构的总刚度矩阵。

?K?1?EA1??l1?1?EA22 ??K??l2??11?1?1?1???11? ??整个结构有三个节点，首先将单元刚度矩阵扩充为3X3的矩阵，移动各元素使之与单

刚矩阵中的元素位置相对应，如下：

?K?1?1?10?EA1?? ?K?1?EA2??110?l1?l2??000??然后直接相加。

0??00?01?1? ????0?11???EA1??l1EA?K????1?l1??0???P???R??0?

?0???EA?1l1EA1EA2?l1l2EA2?l2?0??EA2? ??l2EA2??l2??2， 把各单元的节点力向量组集成总的节点载荷向量。

3， 根据边界条件，修改总刚度矩阵，获得总刚方程组。

武汉理工大学教务处制 4

武汉理工大学教师备课专用纸

边界条件修改之前的总刚方程：

?EA1??P??l1???EA1?0?????0??l1???0??EA1l1EA1EA2?l1l2EA2?l2?????u1?EA2?????u2? ?l2???u3?EA2??l2??0修改以后（采用置“0，1”法）

?EA1?l?P??1???EA1?0?????0??l1??0???EA1l1EA1EA2?l1l20??0???u1???0??u2? ???1??u3????4， 求解方程组，得出总的节点位移向量。

解得的解是：

Pl2??Pl1??EAEA?2u?1??1?Pl2????u??2??? ?u??EA2??3??0?????有了节点位移，再回代到前面单元推导过程中的公式×和××，就可以求得每个单元的

应变和应力了。

从这个简单的例子，我们了解了有限元法求解力学问题三大步骤中的内容，想必很多同学会说，这样复杂，如果运用材料力学的知识，我还来得快些。但是大家不要忘记，有限元的计算很多都是编程完成，而且现在很多的商业软件都已经完成了很多的工作。我们学习有限元主要是了解它的原理，并对常见单元的力学特性有所了解，这样对于以后运用有限元起到帮助作用。

所以下面章节的内容，就是围绕这个主题展开。要达到这个目的，我们还必须学习必要的弹性力学知识。对弹性力学知识的学习，也对我们以后把握问题的本质有帮助。

武汉理工大学教务处制 5

武汉理工大学教师备课专用纸

第二章 平面问题

平面问题在力学研究的课题中属弹性力学的范畴。该类问题不仅本身具有典型性，而且在机械零构件的分析中，也是应用得非常广泛。所以这类问题也称之为经典的力学问题。

我们知道，实际的机械零构件都是具有三维空间尺寸的物体，理应作为三维对象处理，但是当物体的几何形状和受力状态处于某些特定的情况下，近似地简化为平面问题，不仅可以大大简化计算的工作量，而且其精度也完全能够满足所要求。如：直齿圆柱齿轮可在垂直与孔轴线的截平面内作平面应力分析就足以了解整个齿轮的受力状态；大坝的横断面可作平面应变分析来了解整个大坝受力情况等。

本章是全书的重点，在这里不仅介绍弹性力学的基本知识。还将系统地讲解有限元的基本概念、原理和方法。是学习以后各章节的基础。

§2.1 外力、应力、应变和位移

外力、应力、应变和位移的概念在材料力学中已经学习过，由于这些概念在弹性力学、有限元法中具有和在材料力学中不同的规定，弹性力学中的规定和有限元法是完全相同的，所以在这里我们将按照弹性力学的习惯表达方法把他们集中的加以阐述。 一，外力

外界作用在物体上的作用力，可以分为两大类： 1） 体积力

分布在物体体积内的力。如：重力、惯性力和磁性力等。

单位体积的体力在坐标轴上的分量X、Y、Z，称为体力分量。符号规定为沿坐标轴正向的为正，沿负向的为负。 2） 面力

作用在物体表面的力。如轮压、水压等。 它又可细分为集中力与分布力。面力在坐标轴上的投影，表示为X、Y、Z。符号沿正轴为正，负轴为负。 二，应力

弹性体受到外力作用后，内部产生的抵抗变形的内力。

以弹性体中P点为定点的微单元体来考察。所谓微单元体，就是图中PA、PB、PC的边长分别为dx、dy和dz。以下简称这样的微单元体为微元体。

微元体每个面上的应力都可以分解为三

个应力分量。以图中红面为例，分别是σx、τxy、τxz。

应力命名的规则：以应力所在面垂直的坐标轴为第一个下标，应力指向为第二下标。如果下标相同就用一个下标表示。

符号规定：正面上的应力与坐标轴同向为正，反之为负。负面上的应力与坐标轴反向为正。反之为负。

所谓正面就是面的外法向与坐标轴同向为正。反之为负面。

作用在两个相互垂直的面上，并且垂直与该两面交线的剪应力互等。即： τxy=τyx；τyz=τzy；τxz=τzx

如此以来，代表P点应力状态的应力分量应有6个，它们是：

武汉理工大学教务处制 6

武汉理工大学教师备课专用纸

??????x?y?z?xy?xz?yz?T

三，位移

任一点的位置移动用u、v、w表示它在坐标轴上的三个投影分量。

?????uvw?

T符号规定：沿坐标轴正向为正，反之为负。 四，应变

弹性体内各点的位移在受力后一般是不相同的。各点之间距离的改变，从而使物体形状发生变化，即所谓的变形。而物体的形状总可以用它各部分的长度和角度表示。

长度的改变称为正应变ε，角度的改变称为剪应变γ。以微元体三个棱边的线伸长和角度的变化，就分别有和6个应力分量相对应的6个应变分量，即：

??????x?y?z?xy?xz?yz?T

为与前面符号规定一致，这里对剪应变的符号规定如下：

正应变伸长为正，缩短为负；剪应变使直角变小为正，变大为负。

§2.2 两类平面问题

前面我们讲过，实际受力物体都是三维的空间物体，作用在其上的外力，通常也是一个空间力系，其应力分量、应变分量和位移也都是x、y、z三个变量的函数。但是当所考察的物体具有某种特殊的形状和特殊的受力状态时，就可以简化为平面问题处理。

弹性力学中的平面问题有两类。 一，平面应力问题

当物体的长度与宽度尺寸，远大于其厚度（高度）尺寸，并且仅受有沿厚度方向均匀分布的、在长度和宽度平面内的力作用时，该物体就可以简化为弹性力学中的平面应力问题。

我们分析以下其应力特征。

当z=±t/2时，有σz=0、τzx=0、τzy=0。 由于板较薄（相对于长度和宽度尺寸），外力沿板厚又是均匀分布的，根据应力应连续的假定（弹性力

学中的基本假定），所以可以认为，整个板的各点均有σz=0、τzx=0、τzy=0。如此以来，描述空间问题的6个应力分量也就变为了3个，即

??????x?y?xy?T

而且这些应力分量仅是x、y两个变量的函数。 二，平面应变问题

当物体是一个很长、很长的柱形体，其横截面沿长度方向保持不变，物体承受平行于横截面且沿长度方向均匀分布的力时，该问题就可以简化为平面应变问题处理。

分析其应力特征。假定其长度方向为无限长，那么任一横截面都可以看作是物体的对称面，如此

武汉理工大学教务处制 7

武汉理工大学教师备课专用纸

则有该面上的点都有w=0，也就是横截面上的所有点都不会发生Z方向的位移。由这一点可以推出也就有εz=0、τzx=0、τzy=0。

和平面应力相比较，平面应变是εz=0，那么是否也就有σz=0呢？ 可能有同学想σ=Eε，当然也就有σz=0，这是错误的。

平面应变状态下σz≠0的。虽然不等于零，但它也不是一个独立的变量了，它由σx、σy的大小而决定。如此以来，独立的应力分量同平面应力问题一样也是3个：

??????x?y?xy?T

三，两类问题的比较

1， 几何特征

平面应力 厚度＜＜长度、宽度

平面应变 厚度＞＞长度、宽度 为便于说明可讲上述长度看作为厚度 2， 受力特征

外力都必须在其面内且不沿厚度方向变化 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

?????????yx?yxdy???x?xdx?dy?1??x?dy?1????dx?1??yx?dx?1?X?dydx?1?0?x?y????

展开化简得

???x????yx???X?0????????x???y????????同理可求得∑Fy=0满足得条件，? ?y???xy??Y?0??y???x?????由∑M=0，列出方程如下：

??xy???yx??dxdx?dydy??????dy?1???dy?1????dx?1???dx?1??0化简后得： xyxyyxyx?????x?22??y?22?武汉理工大学教务处制 8

武汉理工大学教师备课专用纸

?xy?1??xy1??yxdx??yx?dy

2?x2?y略去微量项，可得：?xy??yx。这就是前面所将的剪应力互等。

对于平面应变问题，微元体的前后面还有正应力σz，不过它们是互等的。对于推导出来的结果，没有任何影响。所以平面问题的平衡微分方程就是：

???x????yx???X?0????????x???y????y???y?????xy? ?????Y?0???x????写成矩阵形式

???x??00??y??y??x?????x??X????y?????0 ?????Y??xy?§2.4 几何方程

考察平衡微分方程，其中具有三个未知变量σx、σy、τxy，，而只有两个方程，方程具有无数个解。表明仅从静力学关系无法求解该方程。我们必须从其它方面寻求帮助。

弹性体在受到外力后，会发生位移和形变，从几何上描述弹性体各点位移于应变之间的关系，就是弹性力学中的又一个重要方程——几何方程。

仍然取截面的微元体ABCD，AB、CD边长为dx、dy，厚度为“1”。 位移u、v都是x、y的函数，即u(x,y)、v(x,y)，偏导数

?u?v、表示位移分量u、v沿?x?x坐标轴x的变化率，偏导数

?u?v、表示位移分量u、v沿坐标轴y的变化率，设A点的位?y?y移为u、v，那么B’点的位移就是：

uB'?u??u?vdx vB'?v?dx ?x?x武汉理工大学教务处制 9

武汉理工大学教师备课专用纸

同理的D’点的位移分量

uD'?u??u?vdy vD'?v?dy ?y?y由于α角在位移和形变很微小的情况下非常小，所以

A’B’≈A’B”

线段AB位移后的总伸长量为 A’B’-AB=A’B”-AB=uB’-uA=u?∴?x??u?udx-u=dx ?x?x

?v?v?u?udy/dy?，同理可得?y? dx/dx??y?y?x?xv?v??dx??v?x剪应变由α、β两个角度组成

B'B\??tg??'\?AB?u??dx??u?dx?u??x????v?x1??u?x

由于

?u?u?v ?1，所以??，同理可得???y?x?x?v?u? ?x?y∴?xy?????综合以上几何方程，并将它们写成矩阵形式：

??x????x?????y???0??????xy???y0??u??? ?y???v?????x?由以上方程可以看出，当弹性体的位移分量确定以后，由几何方程可以完全确定应变，反过来，已知应变却不能完全确定弹性体的位移。这是因为物体产生位移的原因有两点：

1） 变形产生的位移。 2） 因运动产生的位移。

因此弹性体有位移不一定有应变，有应变就一定有位移。

§2.5 物理方程

描述弹性体内应力与应变关系的方程，我们称之为物理方程，也叫材料的本构方程。 弹性力学通常研究的是各向同性材料，在三维应力状态下的应力应变关系。当弹性体处于小变形条件下，正应力只会引起微元体各棱边的伸长或缩短，而不会影响棱边之间角度的变化，剪应力只会引起角度的变化而不会引起各棱边的伸长或缩短。因此运用力的叠加原理、单向虎克定律和材料的横向效应（泊松效应），我们就可以很容易的推导出材料在三向应力状态下的虎克定律，也就是通常所说的广义虎克定律。

武汉理工大学教务处制 10

武汉理工大学教师备课专用纸

式中，E——材料线弹性模量 G——材料剪切弹性模量

μ——材料横向收缩系数，即泊松系数。 三者不是独立的。具有以下关系：

G?E

2?1???这些参数都是材料的固有属性系数，可以通过查材料手册获得。如：钢材的弹性模量E＝196～206GPa之间，通常取2.1×105MPa，μ＝0.24～0.28之间，也取为0.3进行计算，G＝79GPa。

将以上空间问题的物理方程运用到平面问题，其形式如下： 1） 平面应力问题的物理方程

前面分析已知，平面应力问题有σz=0、τzx=0、τzy=0。所以：

1??x???y? E1?y???y???x?

E1?z????x???y?

E12?1????zy??xy??xy GE?x?从以上物理方程，也论证了我们前面说的εz≠0的结论，但由于它是由x和y方向应

力产生的附加无约束变形，所以通常不予以考虑。

在有限元分析中更多的是运用应变表示的应力关系，所以我们将上式变形一下：

?x?E??x???y? 21??E??y???x?

1??2?y??xy?E?xy

2?1???武汉理工大学教务处制 11

武汉理工大学教师备课专用纸

以上方程的矩阵表达形式为：

??x?E????y??21??????xy???1????0???0???x???10???y? 简记为：?????D????

1????????xy?02?式中：{σ}、{ε}为该问题的应力、应变向量。

[D]为弹性矩阵。它是一个对称矩阵，且只与材料的弹性常数有关。 2） 平面应变问题的物理方程

因为εz=0所以由空间物理方程的第三式得：

?z????x??y?，代入（1）、（2）式得 1??2?x?E1??2?y?E????????x1??y?? ??????????y1??x?? ???zy?21?1?12?1??????xy??xy??xy EGE1??2??同理变型为应变表示应力的形式：

?x????????y?? 2?x1????????1???1????????????x?? 2?y1????????1???1?????E1??2E1??2?y??xy?E1??2???2?1??1???????xy

矩阵形式：

武汉理工大学教务处制 12

武汉理工大学教师备课专用纸

E??x?1??2????y??2???????xy?1???1?????????1????1?????0???1??10??0????x????0???y?

????xy????1?1????2?????D???? 也可简记为?平面应变问题的弹性矩阵不同于平面应力问题的弹性矩阵，比较可以发现只需将平面应

力问题弹性矩阵[D]中的材料常数E换为E/(1-μ2)，μ换为μ/(1-μ)就得到了平面应变问题的弹性矩阵。其实弹性矩阵的这种转换方法，是弹性力学中将平面应力结果，转换到平面应变问题结论的一般方法。因为在两种平面问题的描述方程中（平衡微分方程、几何方程和物理方程），只有物理方程是不同的。

§2.6 边界条件

求解弹性力学问题实际就是在确定边界条件下，求解8个基本方程（平面问题而言），以确定8个未知变量。所以从数学的角度看，就是求解偏微分方程的边值问题。边界条件的给出通常是各式各样的。大体可以分为三类： 1， 第一类边值问题

给定物体的体力和面力条件，确定弹性体的应力场和位移场。此类问题边界以力的形式给出，所以也称为应力边界条件。我们可以来考察一下应力边界的一般形

式：

上给出的力的

?jivj?Ti Ti是在Sσ面

分量。

平面问题如左图所示，设阴体弧长为ds，厚度为单元厚度

X轴的夹角为θ，由阴影部分微元体的平衡条件可以推出：

影部分的微元“1”，其法线与

?Xds?1??x?dscos??1??xy?dssin??1?0??? Yds?1???dssin??1???dscos??1?0yxy??化简后得：

??xcos???xysin??X??? 此即为平面问题应力边界方程。 ?sin???cos??Yxy?y?2， 第二类边值问题

给出弹性体的体力和物体表面各点得位移条件，确定弹性体得应力场和位移场。由于以位移给出已知得边界条件，所以也称为位移边界问题。

武汉理工大学教务处制 13

武汉理工大学教师备课专用纸

一般得位移边界条件为：ui?ui 在Su面上。

3， 第三类边值问题

给定弹性体得体力和一定边界上得面力，其余边界上的位移，确定其应力场和位移场。由于边界以力和位移两种形式给出，所以也称为混合边界问题。 针对不同的边界条件，弹性力学求解的方法也有所不同。

§2.7 弹性力学的解题方法（解析法）

1， 应力法

由于第一类问题的边界条件以应力形式给出，所以以应力作为基本的未知量的求解过程，就是人们通常所说的应力法。 由于平衡方程中有三个未知量，而只有两个平衡微分方程，必须找出另外一个包含应力分量的方程，才能求得方程解。

考虑到弹性体变形前是一个连续体，变形后也应是连续体的基本假设，所以要求微元体的变形一定要协调，才能使变形前、后，不会发生裂缝、重叠等现象。要使变形协调，就要研究几何方程。

前面介绍的平面问题几何方程如下：

?x??v?u?v?u? ?y? ?xy??y?x?y?x分别对εx、εy、求y、x的二阶偏导，然后相加：

2?2?x??y?2??222?y?x?y2??u?????2??x??x2??v??2??u?v???xy???y????x?y???y??x????x?y ????上式表明三个应变分量之间应满足的连续性条件，我们称之为形变协调方程（相容方程）。

通过物理方程，使上述的形变协调方程换成应力表示的形式，使之与平衡微分方程就构成了应力法中需要求解的方程组。具体我们来看看：

1） 利用物理方程消去相容方程中的形变分量（以平面应力为例）

?2?1?????xy1??2?2????????????? yyx???2x2E??yE?x?y?x?2?2?xy??2??2 ?2??x???y??2??y???x???2?1????x?y?y?x??2） 利用平衡微分方程，消去上述公式中的剪应力

??xy?y??xy?x????x1 ?X ○

?x??y?y2 ?Y ○

??1式对X求偏导，○2式对y求偏导，然后两者相加 ○

武汉理工大学教务处制 14

武汉理工大学教师备课专用纸

2?2?x?X??y?Y 2?????22?x?y?x?y?x?y?2?xy代入相容方程，化简

??2?x?X?2?y?Y??2?2??x???y??2??y???x???1?????2??2??2?x?y??y?x?y???x?22?2?x??y??y?2?x??X?Y??? ???????1???2222???x?y?x?y??x?y???2??X?Y??2??????????????1??y??x??y?? ??x2?y2?x????对于平面应变而言，运用前面讲过的物理方程的转换方法，只需将上式中的μ代以μ/(1-μ)

就可以了。

??2?2?1??X?Y????????????y??x??y?? ??x2?y2?x1??????3） 最终求解的方程组

平面应力问题 ???x????yx???X?0????????x???y????y???y?????xy? ??????x???Y?0?????2??X?Y??2??????????????1??y??x??y?? ??x2?y2?x????平面应变问题

???x???x???y???y? ????yx??????y???X?0???????xy? ??????x???Y?0?????2?2?1??X?Y?????????????? y????x2?y2?x1????x?y???三个微分方程，三个未知变量，再考虑边界条件，即可求得。

如果是单连体（只具有唯一的封闭边界）的对象，满足了以上方程组后就是实际的解。但对于多连体（具有多个封闭边界）的对象，中还包含有待定系数，这些待定系数会导致位移的解出现多值性。所以对于多连体的问题，还应考虑位移的单值条件，才能最终确定。该部分的内容可以参见徐芝纶编的《简明弹性力学教程》中圆环受均布压应力的情况（P87） 2， 位移法

位移法主要针对第二类边界条件问题求解。

武汉理工大学教务处制 15

武汉理工大学教师备课专用纸

解题步骤：

1）改写物理方程使之成为应变表示应力的形式 2）应用几何方程表示以上得到公式中的应变 3）将它们代入平衡微分方程

经整理最后得到的位移法求解平面应力问题方程为：

E1??2E1??2??2u1???2u1???2v????x2?2?y2?2?x?y???X?0 ????2v1???2v1???2u????y2?2?x2?2?x?y???Y?0 ??两个未知量，两个方程，再加以边界条件即可求得问题的解。

以上介绍的解析法中，应力法和位移法是求解弹性力学问题的基本方法。但都需要解联立的偏微分方程组。求解过程中的数学难度，常常导致这种求解是无法进行的。

由于应力法在体力为常量的情况下可以进一步简化为求解一个单独的微分方程的问题，所以应力法在解析法中相对应用较多。但即使这样，在应力法中，也常常采用逆解法或半逆解法。

§2.8 常体力情况下应力法的简化、应力函数及实例分析

我们前面讲述了弹性力学的三大方程，及应用这三大方程的应力法和位移法解题步骤。但是也说了要求解这些联立的偏微分方程在数学上是存在很大难度的。很多情况下，根本无法进行。那么弹性力学如何在实际中进行应用，它们和我们前面学过的材料力学区别究竟在哪里？我们将通过这一节的学习，一方面了解如何应用这些弹性力学的方程求解问题，另一方面加深对力学概念的理解，建立力学分析问题的直观感觉，为建立有限元模型打好基础。

我们知道在大多数情况下，我们分析的对象，体力是常数，它不随x、y坐标变化。如此以来，前面讲解的第三个方程（应力表示的相容方程），就可以简化为了：

??2?2???x2??y2??2??? 简记为：??x??y??0 ????0?xy??以上方程称为拉普拉斯微分方程，数学上也称之为调和方程，满足调和方程的函数称之

2为调和函数，及这里的?x??y。?是拉普拉斯算子。

??这样以来常体力情况下的应力法方程就是： ???x????yx???X?0????????x???y????y???y?????xy? ??????x???Y?0????2??x??y??0

以上方程都不含有材料常数E、μ，所以平面应力和平面应变两类问题具有相同的方程，

这表明：在单连体问题中，只要边界相同、受同样的分布外力，应力分布与材料无关；也与是平面应力还是平面应变的状态无关。

以上结论的意义：

1弹性力学平面解答的应用范围加宽。 ○

武汉理工大学教务处制 16

武汉理工大学教师备课专用纸

2为实验应力分析提供了理论依据（光弹实验） ○

下面我们考察平衡方程：

???其解由???x???其次微分方程的通解，加上任意一组特解组成。 ?yx??0??????x???y????y????xy? ???y??????x???0????特解我们可以很容易找到。如：?x??Xx ?y??Yy ?xy?0 所以现在关键是找其次方程的通解。

???由第一个方程，可得：???x?????yx??????x???y???????xy?，由数学微分理论，该式是一个函数全微?y分的充要条件。所谓全微分就是有一个函数

dA??xdy????xy?dx 且 ?x??A?A ??xy? ?y?x同理由第二式可得：?y??B?B ??xy?

?y?x由剪应力公式又知存在一个函数φ，可以使d??Bdy?Adx ∴A????? B? 故： ?y?x?2??2??2? ?x?2 ?y?2 ?xy???y?x?x?y由于应力与函数φ存在这样的关系，因此函数φ即是应力函数。

我们用应力函数来表示相容方程：

??2??2?2??2???2?2????x2??y2????x??y?????x2??y2?????x2??y2????0 ???????2?2???4??0

上式表明φ是重调和函数。

前面讲过在弹性力学中，常用逆解法和半逆解法。

所谓逆解法就是设定各种满足相容方程的应力函数，运用σx、σy与φ的关系，求得应力分量，再考察其满足何种边界条件，从而知晓这样的应力函数可以解决什么问题。

所谓半逆解法就是根据弹性体的边界形状和受力关系，设定部分应力分量为何种形式的函数，从而确定应力函数φ，在运用应力函数求出所有的应力分量，根据边界条件确定应力分量应具有的最终形式。

下面我们来看一个半逆解法的例子。

运用逆解法求简支梁受均布载荷的应力分布。

武汉理工大学教务处制 17

武汉理工大学教师备课专用纸 由材料力学知，弯曲应力主要由弯矩M引起，剪应力由剪力引起，而挤压应力由分布载荷q引起。现在q为不随x变化的常量。因此我们设σy不随x坐标变化，即?y?f?y?，

?2??f?y? 因此 ?y?2?x我们对x积分：

???xf?y??f1?y? ?xx2??f?y??xf1?y??f2?y?

2上式中f1(y)、f2(y)是待定函数。由于应力函数必须满足相容方程，所以：

?4??4??4?1 ???4?222?4?0 ○

?x?x?y?y4?4??4?d2f?y??0 22? 42?x?x?ydyd4f1?y?d4f2?y??4?x2d4f?y?1中 代入到式○??x?44442dy?ydydyd4f1?y?d4f2?y?1d4f?y?2d2f?y?x?x??2?0 44422dydydydy考察上式可以看出它是一个x的二次方程，所以一般情况下只有两个根。也就是说只有两个位置能够满足上式。但我们对相容方程的要求是绝对满足。也就是要求在整个梁的范围内都满足。所以只有该方程的系数项和自由项全部为零。即：

d4f1?y?d4f2?y?d2f2?y?d4f?y??0 ?0 ?2?0 4442dydydydy2 ∴f?y??Ay?By?Cy?D ○

32武汉理工大学教务处制 18

武汉理工大学教师备课专用纸

3 f1?y??Ey3?Fy2?Gy ○

d4f2?y?d2f?y???2??12Ay?4B

dy4dy2f2?y???A5B44 y?y?Hy3?Ky2 ○

1063、○4中分别省掉了常数项和一次项、常公式中的A.、B、C…K都是待定系数。公式○

数项。这是由于f1(y)和f2(y)分别是应力函数中x的一次项和常数项的原因，这样处理不会对应力分量产生影响。

最后求出的应力函数为：

x2AB??Ay3?By2?Cy?D?xEy3?Fy2?Gy?y3?y4?Hy3?Ky2

2106????由应力与应力函数的关系，可以求出各个应力分量：

?2?x2?x?2??6Ay?2B??x?6Ey?2F??2Ay3?2By2?6Hy?2K2?y?y?Ay3?By2?Cy?D

?xy??x?2Ay2?2By?C???3Ey2?2Fy?G?

由于以上求得的应力分量满足了平衡方程和相容方程，所以只需根据边界确定A…K的系

数，就求得了该问题的解。

根据对称性，知道为偶函数，为奇函数，所以有E=F=G=0

通常梁的跨度远大于梁的深度，因此上下边界是主要边界，它们必须满足。

???yy?h2?0 ??y?y??h??q ??xy?y??h?0

22将它们代入的表达式，并且考虑E=F=G=0

h3h2hA?B?C?D?0 842h3h2h?A?B?C?D??q 84232hA?hB?C?0 432hA?hB?C?0 4以上四个方程解四个未知数，求得：

A??有：

2q3qqC?D?? B=0 将他们代回到应力分量的表达式中，也就32h2h武汉理工大学教务处制 19

武汉理工大学教师备课专用纸

6q26q3xy?y?6Hy?2K 33hh2q3qq?y??3y3?y?

2h2h6q3q?xy?3x2y?x

2hh?x??左右两个边界，由于前面已经考虑了对称性，所以这个仅考虑优边界。

没有水平力。要x=l时，σx=0，考察σx的表达式，除非q=0。而这和已知条件相违背。所以在这个边界上我们只能要求部分满足。运用圣维南原理运用等效力系代替它。（这样产生的误差只在力作用点附近较大）。运用的等效力系就是合成力系为平衡力系：

?h2h2????h2xx?ldy?0 合力等于0

?h2????xx?lydy?0 合力矩等于0

ql2q由第一个条件得K=0，（奇函数在对称区间上的积分为零）；由第二个条件得H?3?

10hh可以证明剪应力的合力为-ql。即

h2?h2????xyx?ldy??ql

最终求得的结果，加以整理：

6q2y?y23?2 ?x?3?l?x?y?q?42????h?h5?hq?y??2y??y???1???1??2?h??h?

2?xy6q?h22????3x??y ??h?4?h3h2y2?由于厚度为“1”，此时其惯性矩I?，静矩S?（计算见图）

1282任意一点的弯矩M?ql?l?x??剪力Q??ql?q?l?x???qx

q?l?x?2?ql2?x2 22??武汉理工大学教务处制 20

武汉理工大学教师备课专用纸

所以上式中的应力分量可以改写为：

My?y23? ?x?y?q?42????Ih?h5?q?y??2y??y???1???1??2?h??h?

QS I2?xy?各项应力的分布

σx第一项为主应力项，与材料力学中的结果完全一致。第二项为应力修正项。当L/h＞4时，仅占主项的1/60；当L/h=2时，仅占主项的1/15。所以对于深梁的工程构件不容忽视第二项的影响。

§2.8 虚功方程

从前一节深梁的例子，我们可以看到，弹性力学解析求解的过程是非常复杂的。这样的求解对实际工程情况说来，很多情况根本是不可能的。所以长期以来，技术人员就一直探求数值求解的方法。有限元法是其中最成功的方法。为分析单元特性和简化分析过程，我们还需了解单元的能量关系。因为在力学上，很多时候从能量的角度分析，可以大大简化分析的过程。

一，应变能的概念

有材料力学知，弹性体在受到外力作用发生变形的过程中，弹性体内部会存储变形势能

武汉理工大学教务处制 21

武汉理工大学教师备课专用纸

——应变能。在单向应力场中，单位体积的应变能的计算可以表示为：

1dU???

2对于平面问题，有三个应力分量和与之对应的应变分量。由于在小变形情况下，正交力系互不影响，由力的叠加原理，所以该种情况的应变能为：

dU?1??x?x??y?y??xy?xy??1???T????1???T??? 222其中：??????x?y?xy?T ??????x?y?xy?T

11TT?????????dv ??dv???2?2vv整个弹性体的应变能：U?上式也表示应力在应变上所做的总功。

二，虚功原理和虚功方程

在理论力学中，我们曾经学到一个虚功原理，也称虚位移原理。其基本思想就是：假定加在物体上一个可能的、任意的、微小的位移，在平衡条件下，物体的约束反力所做虚功应等于外力所做虚功。因为能量必须守恒。

在这里所说的可能的虚位移是指位移必须满足的约束条件。任意的是指位移的方向和类型是随意的。

把这一原理运用到现在的弹性体中，衡量弹性体应满足的平衡能量关系就是：假定加在弹性体上一个可能的、任意的、微小的位移，在平衡条件下，弹性体内的应变能应等于外力所做虚功。同样是因为能量必须守恒。

运用这一原理，我们可以推到有限元中广泛用到的虚功方程。

假定弹性体发生u、v的虚位移，则由几何方程得：

**?x*?v*?u*?v*?u***?? ?y? ?xy?

?y?y?x?x现考察弹性体微元体和边界处微元体上的力所做虚功：

武汉理工大学教务处制 22

武汉理工大学教师备课专用纸

1内部微元体上的力所做虚功 ○

左面的应力虚功???xdy?1??u

*?*?u*???x??右面的应力虚功??x?dx?dy?1??u?dx?

?x?x????左、右两面上的虚功之和（略去高阶微量，并考虑?x*?u*?） ?x??x*??*?u?dv1 ??x?x?x??同理得剪应力的虚功之和

??u*??xy*????xy?y??yu??dv1 ??体力X的虚功Xudv1

同样的考虑Y方向的?y以及?yx以及体力Y的虚功，然后叠加成内部微元体上的虚功如下：

????x??xy?*???y??yx?*?**2边界上的○???dw1?????X?u????Y?v?dv1??x?x??y?*y??xy?xydv1??y?x????y?????x?*??微元体

设斜边中点处的虚位移u、v，形心处的应力为?x、?y、?xy那么在直角边上的应力和位移均有一个负的增量，如图所示，虚功计算为：

**?*?u*dx???xdx?????x??dy?1???u??x2?? ?x2????????*?dx?????xu*??xu*??x?x??dy?1（略去了高阶微量） ??x?2??同理的dy直角边上的剪应力虚功为：

????xy*?u*?dy?*???xyu????yu??xy?y??2?dx?1 ??????代入已有的几何关系：dx?dssin?，dy?dscos? 斜面上面力的虚功

体积力的虚功

武汉理工大学教务处制 23

武汉理工大学教师备课专用纸

同样的方法求得另一面上正应力与剪应力的虚功，全部相加即得斜边微元体上的虚功之和为：

????x??xy?*???y??yx?*?**?u???v?dv2??x?xdw2?????X??Y??y?*y??xy?xydv2??x????y?x????y?????????X??xcos???xysin??u*??Y??ysin???xycos??v*

支反力处的虚位移为零，所以支反力不做功，将dw1+dw2，并对整个体积积分，可以得到整

个弹性体内的总虚功：

????x??xy?*???y??yx?*??u???v?dv Wz??????X??Y??x????y?x???y??v????????sv??X???xcos???xysin??u*??Y??ysin???xycos??v*ds ?** ???x?x??y?*y??xy?xydv?根据平衡微分方程和静力边界条件，上式的第一、第二项都是零，所以弹性体的总虚功为：

***Wz???x?x??y?*y??xy?xydv???vv???????dv

T根据能量守恒，它应与外力的虚功相等

****** Wz???x?x??y?*y??xy?xydv??Xu?Yvdv??Xu?Yvdsvvs??????由于该等式引入了平衡方程和边界方程，所以上式虚功方程等价于静力平衡条件（内部和边

界微元体）。不同之处在于它是一种能量的表示形式。为了便于有限元中方便运用，引入广义力和物理方程，虚功方程变形为：

????F???????D????dv

**v综合以上推到过程，虚功方程表达的物理概念就是：“若弹性体处于平衡状态，那么外力在虚位移上作的虚功应等于应力在应变上作的虚应变功，或者说等于虚应变能”。

§2.9 平面问题单元划分

有限元法在平面问题进行分析时，才采用三角形单元和四边形单元、或者矩形单元，三角形单元的优点是简单且对结构的不规则边界逼近好，而矩形单元却更能反映实际弹性体内部的应力应变变化。这两点我们会逐渐向大家说明。

所以一般说来，有限元分析，单元划分的密度和单元种类选取，对计算结果起重要作用。一般单元划分越密集，结果越精确。单元多也导致求解的线性方程组阶数增高，要求计算机的内存也更大，计算的时间也越长，分析的效率就越低。解决这一矛盾的方法就是在应力集中区域单元划分密集一些，应力变化梯度小的位置，划稀疏些，这样就能兼顾精度与效率的关

武汉理工大学教务处制 24

武汉理工大学教师备课专用纸

系。

一般的原则是： 1） 根据结构的受力和支承特点，按对称和反对称的性质，简化分析模型，以减少

计算分析的规模。

2） 合理布局单元的密集程度，以使计算结果精度高而计算量小。 3） 在同一单元内，单元的特性数据和材质数据应保持一致。 4） 集中载荷的作用点和载荷密度突变处应有节点。 5） 在欲知道应力状态、内力情况和位移值的位置应有节点。 6） 单元的选取欲分析的目标密切相关。

模型的单元划分好后，把所有的单元和节点按一定的规律和顺序进行编号，选择适当的坐标系（直角、柱面和球面），以方便确定各节点的坐标值。

§2.10 节点位移、节点力和节点载荷

弹性体在承受外力作用后，其内力的传递实际是通过单元之间的边界来实现的。但我们把结构离散化后，如果单元划分得足够小时，可以看成为其内力的传递通过单元与单元之间的节点进行传递。对于平面问题而言，每个节点都有位移和力两个未知量，这两个量又都是x、y的函数，注意平面问题的节点是不能传递力矩的，为什么？ 一，节点位移

对三节点三角形单元而言，因有三个节点，每个节点的位移都有x，y两个分量，所以一共有6个自由度。单元节点位移向量可表示为：

???e??uiviujvjumvm

?T二，节点力

所谓节点力，就是单元对节点或节点对单元作用的力，它是弹性体内部的作用力，也就是我们常说的内力。

和上面相同的道理，节点力向量为：

?F?e??UiViUjVjUmVm

?T三，节点载荷

节点载荷是指作用在节点上的外力，包括： 1直接作用在节点上的外力 ○

2经等效处理后，移植到节点上的等效力。 ○

可用Xi、Yi.表示。由力平衡条件知，节点要保持平衡，那么作用在节点上的载荷应等于节点

武汉理工大学教务处制 25

武汉理工大学教师备课专用纸

内力的合力。即：

Xi??Ui，Yi??Vi，

ee所有的节点载荷向量表示为?Ri????F?

ie§2.11 三节点三角形单元的位移模式和形函数

弹性结构受外载荷后，内部各点的位移变化规律，一般都是很复杂的。很难找到一个函数来描述整个结构内部各点的位移变化规律。但当把整个结构离散以后，在一个相当小的单元内部，却可以用简单的函数来近似描述单元内部位移的这种变化规律。而不致造成很大的误差。这就好比一条复杂的曲线，用一个函数很难描述，但在把这条曲线分段以后，对于每一条分段曲线，却可以用直线或抛物线来描述。 一，三角形单元的位移模式

设单元内任意一点的位移分量u、v是其位置坐标x、y的线性函数，则：

?u?a1?a2x?a3y a1…a6是待定系数。 ?v?a?ax?ay456?改写方程组成为矩阵的形式：

?a1??u??1xy000?????????????? ??v??0001xy??a??6?单元的三个顶点i、j、m的坐标已知，分别为?xiyi?、?xjyj?和?xmym?，因为它们

也是单元上的点，所以应该满足以上假定的位移变化规律。代入上式：

ui?a1?a2xi?a3yi vi?a4?a5x?a6yi uj?a1?a2xj?a3yj vj?a4?a5xj?a6yj

um?a1?a2xi?a3yi vm?a4?a5xm?a6ym

解以上6个方程，求得6个待定系数。

uiuja1?um1xixjxmxiyiyjymyiyjym?1?xjym?xmyj?ui??xmyi?xiym?uj??xiyj?xjyi?um2A??

1xj1xm?1aiui?ajuj?amum 2A??武汉理工大学教务处制 26

武汉理工大学教师备课专用纸

1ui1uj1um1xi1xj1xmyiyjymyiyjym1?yj?ym?ui??ym?yi?uj??yi?y?um 2Aa2???? ?1biui?bjuj?bmum 2A??1xiuiujumyiyjym?1?xm?xj?ui??xi?xm?uj??xi?xi?um 2A1xja2?1xm1xi1xj1xm?? ?同理得

1ciui?cjuj?cmum 2A??1aiui?ajuj?amum 2A1a5?biui?bjuj?bmum

2A1a6?ciui?cjuj?cmum

2Aa4???????ai?xjym?xmyj??bi?yj?ym? ci??xj?xm??i,j,m顺序轮换，A是三角形的面积

??在单元节点的顺序号i,j,m必须是按逆时针排列，否则系数行列式是负值，而三角形的面积为负值，是不合理的。

求得的6个系数可以用以下矩阵表示：

?ai?b?i?a1???1?ci?????2A?0?a??6??0???0二，形函数

000aibiciajbjcj000000ajbjcjambmcm0000??ui??v?0???i?0???uj????? am??vj?bm??um????cm??vm????武汉理工大学教务处制 27

武汉理工大学教师备课专用纸

将所求得的6个待定系数代入位移模式表达式中：

?a1??u?1?1xy000????????????????v?2A?0001xy??a??6??ai0aj0am0??ui??b0b0b??0?jm?i??vi??uj??1?1xy000??ci0cj0cm0?????????2A?0001xy??0ai0aj0am??vj??0bi0bj0bm??um?????0c0c0c??ijm?vm??????ui??v??i?0ai?bix?ciy0ai?bix?ciy0??1?a?bx?ciy?uj????ii?? ?0ai?bix?ciy0ai?bix?ciy0ai?bix?ciy??vj?2A??um?????vm??1?ai?bix?ciy? 令Ni?2A?Ni????就有??00NiNj00Nji,j,m顺序轮换

??Nm0?u??ui?0??i????????N???? ?Nm????v??vm??m?上式就是假定位移模式下导出的单元内任一点位移表达式。

该式的数学意义就是单元内任一点的位移可以由单元节点的某种形式插值得到，其中的插值基函数就是Ni、Nj、Nm。对于我们目前假定的位移模式是线性函数，所以得出的插值基函数也是类似的线性函数。由此可以看出，插值基函数具有反映单元位移变化形态的特征，所以也称之为位移形态函数，简称形函数。[N]就是形函数矩阵。 三，形函数的性质

1单元内任一点的三个形函数之和恒等于1。N?N?N○ijm?1 2在单元的三个顶点处，有 ○

i节点处 Ni?1 Nj?0 Nm?0 j节点处 Nj?1 Ni?0 Nm?0 m节点处 Nm?1 Ni?0 Nj?0

以上这些，可以通过简单的数学运算进行证明。 四，位移模式收敛性的分析

由于位移模式的选取是人为假定的，这种假定只能近似模拟单元内位移的变化规律，由于单元刚度矩阵的推导是以假定的位移模式展开的，那么这种假定的位移模式能否使有限元数值解收敛与精确解，在很大程度上就取决于所选的位移模式，通过数学证明，可以找出位移模

武汉理工大学教务处制 28

武汉理工大学教师备课专用纸

式满足收敛性的几个条件： A）完备性

1） 位移模式必须包含单元的常应变状态。 将u、v的表达式代入几何方程，得：

?x??y??u???a1?a2x?a3y??a2 ?x?x?v???a4?a5x?a6y??a6 ?y?x?v?u?????a4?a5x?a6y???a1?a2x?a3y??a3?a5 ?x?y?x?y?xy?因为系数a2…a6都是常数，所以上面的应变分量也是常量。这也表明所选的位移模式中包含有弹性体的常应变状态。在上面的表达式中不含x、y的变量，说明单元的应变是常量，这也表明这种单元是一种常应力单元。 2） 位移模式必须包含单元的刚体位移。 弹性体位移一般包含两个部分，即变形位移和没有弹性变形的刚体位移。那么作为模拟单元

1中已经证明了弹性应变的位移状态的位移模式，也就应该能同时反映这两部分的位移。在○

变形，下名说明他也包含有刚体位移的特征。 改写位移表达式如下：

a5?a3a5?a3?u?a?ax?ay?a?ax?y?y12312?22 ?a5?a3a5?a3?v?a4?a5x?a6y?a4?a6y?x?x22?1知，a?a?a?a?0，所以上式： 当?x??y??xy?0时，由○2635a5?a3?u?a?y1?2 * ?a5?a3?v?a4?x2?我们再来看看一个点作刚体位移的运动方程。点M先

转到M1，再由M1移到M2，如左图所示：

?u?u0?r?cos??u0?wzy ** ??v?v0?r?sin??v0?wzx比较上面的*和**式，可以看出a5?u0，a4?v0，

a5?a3?wz 2由此得出在三角形线性位移模式中，也包含了单元的刚体位移。 B）协调性

武汉理工大学教务处制 29

武汉理工大学教师备课专用纸

3） 位移模式必须能够反映位移的连续性。 在弹性力学求解问题时，曾经讲到过变形协调方程，也说过它是弹性体变形后仍保持连续和不发生撕裂、侵入缺陷的条件。那么位移模式的选取，也应该保证单元之间不出现撕裂和侵入的缺陷。由于上面假定的位移模式是线性函数，而两点就能决定一条直线。由于相邻单元的公共边界上的位移由与之相连的两个节点插值获得，而相邻的单元具有两个公共的节点，所以通过这两个节点所得的插值值，不可能出现不同。也就是不可能出现下图a和图b的情况。

以上三个条件是选取位移模式必须考虑的。完备形条件是收敛的必要条件；协调条件是充分条件。在有限元中，满足完备条件的的单元是完备单元，满足协调条件的是协调单元。

§2.12 （三节点三角形）单元刚度矩阵的推导

上一节我们已经建立了三角形单元的位移插值模式，并求得了形函数的方程，这样就完成了单元内任一点的位移由单元节点位移表示（插值）的工作，接下来运用我们已学过的一系列知识，我们就可以完成单元刚度矩阵的推导了。 一，推导过程

1） 由位移插值函数导出单元应变的单元节点位移表达式

???e??x??????y?????xy?e?????x??0???????y?0?????u??? ??y?v?????x???????x??0???????y?0?????Ni??y??0????x??0NiNj00NjNm0?u?0??i???? Nm???v??m?武汉理工大学教务处制 30

武汉理工大学教师备课专用纸

??Ni???x???0???Ni???y将Ni?0?Ni?y?Ni?x?Nj?x0?Nj?y0?Nj?y?Nj?x?Nm?x0?Nm?y?0???ui??Nm???????y???

?v?Nm??m??x??1?ai?bix?ciy?代入上式，可得： 2A?ui??v?0??i????uj??cm???

vj???bm??um?????vm?????e?bi1??0?2A?ci?0cibibj0cj0cjbjbm0cm如按分块矩阵记忆，那么：

???e??BiBj?bi1?0Bm???e 其中?Bi??2A???ci?0?ci?? i= i, j, m bi??矩阵[B]称为应变矩阵，或称为几何矩阵。

由以上计算公式知，它与单元的节点坐标有关，但不随点的坐标变化，就是说在这一单元内所有点的应变是相同的。

2） 求得单元应力的单元节点位移的表达式 将几何矩阵代入单元的物理方程，就有：

???e??D????e??D??B????e

弹性矩阵[D]是由材料常数组成的矩阵。令[S]=[D][B]，代入平面应力的物理方程，就有

???e??S????e??1E???2?1???0????bi0?1?10?0??1??2A?c?0?i2?0cibibj0cj0cjbjbm0cm0??ui????cm?????v?bm???m?∴

??biE??b?S??i21??2A?1???ci?2??cici1??bi2bj?cjcjcj1??bj2bm???bj1??2?bm1??2cm??cm?cm?也可以

?1???bm?2?武汉理工大学教务处制 31

武汉理工大学教师备课专用纸

写成分块矩阵的形式

??biE??b?Si??i221??A?1???c?2i???ci? i= i, j, m

?1???bi?2??ci??3） 用虚功方程导出单元刚度矩阵（单刚矩阵）

虚功方程

????F???????D????dv

*T*Tv假定单元的厚度为t，上式改写为单元的虚功方程形式，

??????F???T*ee?????D????tdxdy

*T虚应变也可以用几何方程表示

?????B???? ?????*e*eT*eeT*eAT*e??B??T????T*e???????B?代入上式

T*eT??????F??????????B??D??B????tdA

e由于虚位移元素是常量，所以可以提到积分号以外，并与左边的消去（为什么？）。于是上式变为：

?F?e???A?B?T?D??B????etdxdy

令?K??eT??B??A?D??B?tdxdy，虚功方程就成为了单刚方程

?F?e??K?e???e

由于[B]、[D]都是不含有x，y的常数矩阵，所以双重积分实际就是对面积积分了。

?K?e??B?T?D??B?tA A——是三角形面积

将前面求得的应变矩阵和弹性矩阵代入，然后作矩阵乘法。就得到我们要求得的矩阵计算公

式。我们在这里采用分块矩阵的方法记忆。

?K?e??Kii????Kji??Kmi???K??K?????K??K???

?K??K???ijimjjjmjmmm1???bb?crcs?rsEt2?Krs???24?1???A??bc?1??bcsrrs2??brcs?1???bscr?2? r,s=i,j,m 1??crcs?brbs?2?武汉理工大学教务处制 32

武汉理工大学教师备课专用纸

注意以上是平面应力状态的单刚矩阵，如果是平面应变问题呢？ 二，单元刚度矩阵的性质

单元刚度矩阵[K]e表示单元抵抗变形的能力。它与通常的弹簧刚度系数k的物理意义本质相同，只不过[K]e是一个6×6阶的矩阵。共有36个元素。这是因为三角形的节点力向量和节点位移向量都为6的缘故。 1） 物理意义

分块矩阵[Kij]表示的物理含义是：节点j处产生单位位移，而节点i、m被约束，此时在节点i处产生的节点力。我们写出它的4个分量元素：

?K?ij11?kij??21??kij12?kij 22?kij??在上面的矩阵中，元素的第一个下标表示产生节点力的节点，第二个下标表示产生节点位移

的节点。上标“1”表示水平分量，“2”表示竖直分量。 而且上标和下标的关系是对应的。也就是说第一个下标对应第一个上标；第二个下标对应第二个上标。如此就有：

12就表示节点j处产生竖直单位位移，在节点i处产生的水平方向的节点力。 kij2） 单元刚度矩阵是对称矩阵→?K???TTeTe??K?

这一点可以通过简单的数学证明如下：

??K?????B??D??B?tA?eTTTT??B??D??B???TT??B??D??B?tA??K?e

单元矩阵的对称性，从物理学角度反映出的道理就是，“功的互等”。也就是在节点j处产生某一位移引起节点i处的节点力，应等于在节点i处产生相同位移引起节点j处的节点力。 3） 单元刚度矩阵中的元素只与单元的材料性质、几何形状、尺寸大小有关，而与单元的位

置无关。

单元刚度矩阵中不含有ai、aj、am，上节对位移模式收敛性分析中，曾经说明了a1、a4分别表示单元的平移分量，而

1?aiui?ajuj?amum? 2A1?aivi?ajvj?amvm? a4?2Aa1?由上式知单元的平移运动与ai、aj、am有关，而[K]e又与ai、aj、am无关，所以说它不随坐标轴的平动而变。

4） 单元刚度矩阵是奇异矩阵，即?K??0

e从力学的角度理解单元的刚度方程?F???K????，当给定位移时，可以求得力；当给定

eee力时，却不能求得位移。因为[K]e不存在逆矩阵，在单元没有给出任何约束的情况下，除有应变的可能性外，还同时有刚体位移的可能性。所以方程无解。

武汉理工大学教务处制 33

武汉理工大学教师备课专用纸

§2.13 载荷的节点移置

前面对于有限元模型的分析时，曾经说过，单元之间的力传递是通过节点进行的。所以不在节点上的力，必须按静力等效原则，把它们移置到节点上。静力等效原则：原载荷在任何虚位移上所做的虚功，与移置到节点上的节点载荷所做的虚功相等。

这种处理方法，和我们前面讲到的圣维南原理相同。它们只会影响模型局部的应力分布，而不会影响整个结构的应力。下面根据力的类型，分类说明处理的方法。

1） 集中力的移置

如图所示的受力示意图，M处的力为

?P???Px?R?e??xi等就是：

*TPy?T，移置后的节点载荷为

yixj*Tyjxmym，虚功相

?T??????R??????P?

∵

?????N????，所以有

**eT??????R????N??????P????????N??P?

*TT*eT*e∴?R???N??P?

T由于我们现在一直局限于单元的载荷移置，所以上式中的{R}应记为{R}e。如果将它写成分量的形式：

?R?e??Xi?NiPxYiNiPyXjYjXmYm NjPyNmPxNmPy

?T?NjPx?T从上面公式可以清楚的看到，载荷移置结果与单元位移模式密切相关。 2） 体力移置 体积力密度为?p???XY?，将单元中微元体的体力?p?tdxdy看作集中力，那么

T?R?e???A?N??p?tdxdy

3） 面力的移置

设在单元的一个边界上作用有分布力，面力的密度为?p??X力也看作是集中力，则有

?Y?，将微面积tds上的面

T?R?e??s?N?T?p?tds

以上面力，体力的积分运算较为复杂，在线性模式下，可按照理论力学中静力学平行力分解

武汉理工大学教务处制 34

武汉理工大学教师备课专用纸

的原理，直接求得等效的节点载荷：

W的节点载荷。 32qlt2 i、j边受到集度为q的均布面力载荷，则i、j边节点受到○的等效节点载荷。

213 i、j边受到线性分布的面力，i处集度为零，j处集度为q其合力为R?○qlt，那么i

212处的节点载荷为R，j处的载荷为R。

331均质等厚的三角形单元，受W的体力，则三个节点分别受到○

这些简单的规则，对于今后实际中的应用，可以提高效率。

§2.14 整体分析

该过程是将离散分离的各个单元组集成离散的结构物，从而建立模型的总刚方程。 一，总刚方程和总刚矩阵的组集 1， 总刚度矩阵的组集原则

A， 整个离散结构变形后，各个单元在节点处仍然协调地相互连接。即环绕某个节点的n

个单元，在节点i处具有相同的位移。数学公式的描述：

???????1ii2?????i????i?

nB， 各个节点应满足静力平衡条件。即每个节点上的节点力合力应等于该节点的节点载荷。

数学公式描述：

??F???R?

iie此处的

?e代表环绕节点i的所有单元节点力求和。

2， 在该原则指导下，实例的组集过程

如图a所示的平面问题，采用图b的方法划分网格，单元分析完成后，现将它们组集成原问题的有限元离散网格，求该问题的总刚矩阵。

武汉理工大学教务处制 35

武汉理工大学教师备课专用纸

单元○1为i，j，m=1，2，3（i节点从最小号开始，然后逆时针排列节点号），单刚方程：??F1?1???K1??11??K12?1?K1??F2?1?????K121??13???K???1?1?1??F?1???K?1?3??31?22?1?K1?23?????2??? 将它们展开 K132??K33?1??????3?1???F1?1??K111???1?1??K12?1??2?1??K13?1??13?

?F12?1??K21?1??11???K122?1??2???K123???3? ?F11113???K131???1???K32???2?1??K33???3?1

完全相同的道理，得其他三个单元展开的单刚方程

单元○

2 节点号为i，j，m=1，3，4 ?F21?2??K11???1?2??K13?2??3?2??K14?2??24? ?F23?2??K31?2??21???K33?2??23???K234???4? ?F2224???K241???1???K243???3?2??K44???4?2

单元○

3 节点号为i，j，m=3，5，4 ?F33333???K333???3???K35???5?3??K34???4?3 ?F35?3??K53???3?3??K55?3??5?3??K54?3??34? ?F4?3??K343?3??33???K45?3??35???K344???4?

武汉理工大学教务处制 36

节点号

武汉理工大学教师备课专用纸

4 节点号为i，j，m=2，3，5 单元○

?F2?4??K22?4??2?4??K23?4??3?4??K25?4??5?5 ?F3?4??K32?4??2?4??K33?4??3?4??K35?4??5?4 ?F5?4??K52?4??2?4??K53?4??3?4??K55?4??5?4

实际上，刚度方程可以写成如下的一般形式：

?Fi????Kik???k? k表示单元的节点组成

k?i,j,m2，即各个节点力的合力等于节点载荷 首先我们运用规则○

??F???R?

eiie所以也就有：?F1???F1???R1?

12?F2?2??F2?4??R2?

?F3?1??F3?2??F3?4??F3?3??R3? ?F4?2??F4?3??R3? ?F5?3??F5?4??R5?

我们将上面得到的方程按这个规律相加：

?F1?1??F1?2??R1???K11?1??1?1??K12?1??2?1??K13?1??3?1

222222??K11???1???K13???3???K14???4?1，即对这个等式运用规则○

???????1ii2?????i????i?并整理有：

n?R1????K11?1??K11?2???1???K12?1??2????K13?1??K13?2???3???K14?2??4??R2???K21?1??1????K22?1??K22?4???2????K23?1??K23?4???3???K25?4??5??R3????K31?1??K31?2???1????K32?1??K32?4???2????K33?1??K33?2??K33?3??K33?4???3????K34?2??K34?3???4????K35?3??K35?4???5??R4???K41?2??1????K43?2??K43?3???3????K44?2??K44?3???4???K45?5??5??R5???K52?4??2????K53?3??K53?4???3???K54?4??4????K54?3??K55?4???5?4将上述的

这些方程写为矩阵的形式，就得到了该问题的总刚方程为：

武汉理工大学教务处制 37

武汉理工大学教师备课专用纸

?R1???K11?1?2?R??1??K221??????1?2?R3???K31??R???K?2?4??41???R5???0?K12?1?K13?1?2?K14?2?K22?1?4?K23?1?40?K32?1?4?K33?1?2?3?4?K34?2?3?K43?2?3?K44?2?30?K52?4?K53?3?4?K54?3???1?????K25?4???2???K35?3?4????3?

3???K45????4??3?4??K55?????5??0我们对这个总刚方程进行分析，可以得出以下的规律

3， 总刚矩阵的规律

总刚方程就是运用所有节点的位移向量和总刚矩阵相乘，得到结构的节点载荷，节点载荷向量（外界对弹性体作用的载荷）很容易求得，节点位移向量是未知量，关键就是总刚矩阵如何得出。如都按照上述的步骤，成千上万个节点，该如何作呢？所以有必要总结总刚度矩阵的规律。

1当r=s时，即子块矩阵[Krs]是总刚矩阵主对角线上的子块矩阵。由环绕节点r（s）各单元○

刚度矩阵的相应对角线子块相加。

2当r≠s时，但r、s是相邻单元的公共边上的节点时，子块矩阵[Krs]等于两相邻单元刚矩○

阵的相应子块矩阵相加。如上例中的[K13]、 [K31]等。

3当r≠s时，且r、s只是一个单元的边上节点时，子块矩阵[Krs]就等于该单元刚矩阵的相○

应子块矩阵。如上例中的[K12]、 [K21]等。

2当r≠s时，且r、s不属于一个单元的边上的节点时，子块矩阵[Krs]等于零。如上例中的○

[K15]、 [K51]等。

利用上述规律不仅可以检验总刚组集的正确与否，而且可以直接组集总刚矩阵（手工组集）。但这种方法不利于计算机编程。下面介绍计算机组集总刚矩阵的方法 4， 总刚矩阵组集的步骤

1扩大各单元刚度矩阵，使之成为与总刚矩阵相同的阶数。 ○

2按照总体节点的编号，将各单元刚度矩阵的各个子块移到相应的位置上。其余位置充零。○ 3把各个改造过的矩阵直接相加，就得到总刚矩阵。 ○

4为例讲解步骤○1、○2 我们以单元○

单元4的i、j、m=2、5、3，所以单元刚度方程是

??F2?4???K22?4??4?1??F5?????K25???F?4???K?4?3??32单刚矩阵

?K25?4?K23?4????2?4??4????K55?4?K52?4????5? ????4??K35?4?K33?4???3?4?K?4??K22??1???K25???K?4?324?K25??K55?4?K35?4?K23???K52?4??4??K33??4???K22?4?4??K52???K?4?32?0?0??K25?4?K23?4?K55?4?K53?4?K35?4?K33?4000000??00?00??00?00???武汉理工大学教务处制 38

武汉理工大学教师备课专用纸

?0?0??0??0??00404?K22??K23??K32?4?K33?400?K52?4?K53?400000??K25???4 ?K35???0??K55?4??40二，总刚度矩阵的性质

总刚度矩阵由单元刚度矩阵组集而成，所以也具有单元刚度矩阵的一些性质，如相同的物理意义，位置无关性、对称性和奇异性等，还具有以下性质： 1） 稀疏性

由规律4知，当rs，且r、s不属于同一单元的两个节点是[Krs]=0，表明互不相关的节点数愈多，零子块矩阵就愈多。一般说来相关节点数不超过9，而整个分析对象常常成百上千、上万或几十万。如果整体有100个节点，那么百分之九是非零子块，而百分之九十多都是零子块。所以在总刚度矩阵中非零的子块是很稀少的。 2） 带状性

所谓带状性，就是指总刚矩阵中的非零子块，集中分布在主对角线的两侧，呈带状分布。 半带宽值就是计算这带状宽度的数值。它是以排列元素最长的一行，从第一个非零元素起至主对角线元素止，所有元素的个数。其数值可以由节点的总体编号算出。 B=（相邻节点号的最大差值+1）×（单个节点的自由度个数） 对于平面问题说来，就是：

B=（相邻节点号的最大差值+1）×2

下图中的两种编号方式，可以的分别计算其半带宽值： 图a，B=（2+1）×2=6 图b，B=（6+1）×2=14

武汉理工大学教务处制 39

武汉理工大学教师备课专用纸

半带宽值影响计算机存取总刚矩阵所需的内存大小。愈小愈好。由于半带宽值是直接受节点总体编号的影响。所以我们在建立有限元模型时，应慎重考虑，采用优化的方法。（现在通常在程序中都配有这样的优化程序）。 三，总刚矩阵的压缩存取技术

由于总刚矩阵具有对称性，所以我们只需存入主对角线上半带或下半带的元素，就可以完成解方程的运算。此即为总刚矩阵的半带宽存储方法。 假定取主对角线上半带元素存储，具体做法就是每行元素以主对角线上的元素开始，存储每行半带宽数值的元素个数。如此各元素的行号不变，改变的只是列号。新列号和原列号的关系式如下：

新列号=原列号-行号+1

前面图示的例子，如果采用图示表示总刚矩阵的存储，就是：

武汉理工大学教务处制 40

武汉理工大学教师备课专用纸

如下图所示的情况。可以看出原总刚矩阵存储需要24×24的2维数组，改为半带宽存储后，就成为了24×6的数组了。 另有一维的压缩存储方法，比这种二维的半带宽存储方法压缩更多。在这儿我们就不做介绍了。有兴趣的同学可以看相关的参考书籍。

讲到这儿，实际我们已经知道，每个元素在总刚矩阵中的位置，在节点编号完成后，就完全确定下来了。所以计算机对总刚矩阵的计算，是一部完成的。即“边对号移置、边改列号、边累加”。

四，总体边界条件的处理

前面介绍总刚矩阵的性质时，说明了总刚矩阵是奇异矩阵。即?K??0。就是说总刚矩阵不存在逆矩阵。要求得节点位移的位移解，还必须引入边界条件。 1，边界条件的类型

1节点固定，即u?v?0 ○ij2给定节点位移值。即u?u，v?v ○jjii2， 处理方法 1置“0、1”法

即首先在总刚矩阵[K]中与已知位移分量相对应的行和列元素改为0，但主对角线上的元素改为1，然后在节点载荷向量的列阵中，与已知对应元素的位移用代替，其余元素减去已知位移分别乘[K]中相应元素。

??10?1??R?10?1，已知u1?c1、v2?c4，按照上面的方例：某结构的总刚方程为?K?10?10?法修改如下：

?10?0k22??0k32??00??????0k1020k23k330?k1030??u1??c1??v??Y?kc?kc?0?k210?211244???1??10k310???u2????X2?k31c1?k34c4????????

1?0??v2??c4???????????????????0?k1010?vY?kc?kc10111044???5??5?0?如果展开上式，马上就有u1?c1、v2?c4

第二行 k22v1?k23u2?k25u3???k210v5?Y1?k21c1?k24c4 移项 k21c1?k22v1?k23u2?k24c4?k25u3???k21v05?Y1

从中可以看出，这样处理后不仅可以直接得到u1?c1、v2?c4，而且它们产生的效果也计入到了其他方程中。

如果是固定位移情况，那么载荷向量中的变化是什么呢？（好好想想）

武汉理工大学教务处制 41

武汉理工大学教师备课专用纸

2乘大数法

在总刚矩阵[K]中，把与已知位移位移相对应的行与列主对角线上的元素乘一个很大的数1010，然后把载荷向量中的对应元素代以给定位移乘以相应主对角线上的元素，再同样乘以一个很大的数。 同上例。

?k11?1010??k21?k31??k41?????k101k12k22k32k42?k102k13k23k33k43?k103k14k24k34k44?1010?0k110??u1??c1k11?1010???????k210??v1??Y1??k310??X2?u2???????????

?k410??v2??c4k44?1010?????????????????k1010?Y5?v5????????考察第一个方程：

k11?1010u1?k12v1???k110v5?c1k11?1010

两边同除以k11?10，由于k11?10＞＞k1j，所以

u1=c1

同理可得v2=c4，（思考该方法不能用于固定位移的情况，为什么？） 必须说明的是，当总刚矩阵采用半带宽存储的时候，以上边界条件的引入也应当按照半带宽存储的格式修改，否则就是错误的。

1；当第二种情况时，用方法○2。 实际应用中，当是第一种情况时，采用方法○

五，应力计算

总刚方程在引入足够的边界条件后，就可以进行求解了。所谓足够的边界条件就是要使系统是一个几何不变的系统。数值解方程的方法大多采用高斯消元法。当求解出位移向量以后，就可以通过单元的几何方程求解应变，通过物理方程求解应力了。这个过程称之为回代过程。 1求解应变 ○

1010???e??B????e

2求解应力 ○

???e??D????e??D??B????e

由于三角形单元是常应力与常应变单元，所以由上

述方法求得的应力或应变看作是形心处的应力或应变。而且很明显，这样求得的应力或应变是不连续的。为了推算弹性体内某一点的接近实际应力，我们通常采用以下两种方法，来平滑应力的突变。

1绕节点平均法——就是将绕某一节点的各单元形心应○

力加以平均，来表示该节点的应力。 如左图所示。

武汉理工大学教务处制 42

武汉理工大学教师备课专用纸

??x?1?1??Ax??Bx? ??x?221CDEF ??Ax??Bx??x??x??x??x2??该方法在各单元面积相差不大，单元的内部节点时较为准确。而在外部边界节点上，则不理想。所以对于边界节点多采用三点插值的方法求得。也就是用内部的节点来推算。 2两单元平均法——把相邻两单元的应力相加平均，表示两单元公共边中点的应力值。 ○

??x?2?1??Ax??Bx?

2边界边上的应力值也采用插值方法得到。

至此我们对于有限元求解平面问题的方法全部讲解完毕了。下面通过一个例子来看看，具体如何计算的。

§2.15 例题计算

一，有限元法求解平面问题的步骤：

1） 根据分析对象和给定的条件，按照一定的比例绘出计算简图，该图应有各部位的尺寸、

外力和支承情况。

2） 选择合适的坐标系，划分单元准备数据。数据包括节点数、编号和坐标，单元数、编号

和节点组成；材料特性常数；载荷大小及位置；边界条件。 3） 计算单元的面积、应变矩阵、弹性矩阵和单元的刚度矩阵。 4） 组集总刚矩阵，处理载荷移置，形成载荷向量。

5） 处理边界条件，修改总刚矩阵，得到最后的刚度方程。 6） 求解线性方程组，得到位移向量。 7） 求解应力值。 三，例题

设有如图所示的正方形薄板，在对角线上作用有沿厚度均匀分布的载荷，其合力为2kN，板厚为1个单元厚度，为使计算简便，假定材料的μ=0，求该板的变形。

分析：该薄板对称于对角线，受力也是如此，所以可取其1/4来计算。由于是取1/4，所以该模型的两个直角边，都是对称线。X的对称线不能有x方向的位移，y轴的对称线不

武汉理工大学教务处制 43

武汉理工大学教师备课专用纸

能有y方向的位移。直角的顶点既是x轴对称线上的点，又是y轴对称线上的点，所以该点在x、y方向的位移都应该约束。最后得出的计算模型如下。 解：1准备的数据 A单元节点数据 坐标 X Y 1 0 2 2 0 1 3 1 1 4 0 0 5 1 0 6 2 0 B单元信息数据 I J M 1 ○2 ○3 ○4 ○1 2 3 2 4 5 2 5 3 3 5 6 所有单元的板厚是“1” C材料常数

弹性模量E，泊送比μ=0 D载荷数据

?R???0?10000000000?

TE位移边界

?????0v10v2u31x1v3y100u50u60?

T2计算单元刚度矩阵 单元1 面积A?11x221x3y2?y31101?0.5 2113102b1?y2?y3?1?1?0 b2?y3?y1?1?2??1 b3?y1?y2?2?1?1 c1??x2?x3?0?1?1 c2??x3?x1??1?0??1 c3??x1?x2?0?0?0 形成的几何矩阵

?b1?B?1?1?02???c1弹性矩阵

0c1b1b20c20c2b2b30c30??00?1010?1??

c3??000?100?2?1??2?b3???10?1?101????1?D??E???1???0?????0?100??10??E?010?

?1?1???00?0??2??2?武汉理工大学教务处制 44

武汉理工大学教师备课专用纸

单元1的刚度矩阵

0?0.50.500.5??0.5?0?10?100????0.501.50.5?1?0.5??K?1??B?T?D??B?tA?E??

1.50?0.5?2??0.5?10.5?00?1010????0.50?0.5?0.500.5????同理可求得?K?、?K?和?K?。需要说明的是可以利用位置无关性直接得出?K?和?K?，

23424由于它们是1单元平移得到，所以单刚和?K?一样。

1?K?300?1?00.50.5?0.50.5E?0??002?0??1?0.5?0.5???0?0.5?0.50?0?0.5?0.5??0?0.5?0.5??

10?1?01.50.5??1?0.51.5??0?13组集总刚

??K11?1?1??K21???K?1?K??E?31?0?0???00?0.25?00.5???0.250???0.25?0.5?00?0.50?E??00?0?0?00??00?0?0?0?0?K12?1?K13?1000???K22?1?2?3?K23?1?3?K24?2?K25?2?30??K32?1?3?K33?1?3?4?K35?3?4?K36?4? 0??K42?2?K44?2000??K52?2?3?K53?3?4?K55?2?3?400??44?K63??K66??000??0.25?0.250?0.5???1.50.25?1?0.25?0.25?0.2500.2500??0.251.5?0.25?0.50?0.500.2500??1?0.251.50.2500?0.5?0.2500.25?4引入边

??0.25?0.50.251.500?0.25?100??0.250000.750.25?0.5?100???0.25?0.5000.250.750?0.2500?00.25?0.5?0.25?0.501.50.250.5?0.25??0.250?0.25?1?0.25?0.250.251.50?0.25??000000?0.500.50?000.25000?0.25?0.2500.25??000.250000000000000界条件，采用划行划列（由第一中方法演变而来）就是将已知边界条件的节点位移所在的行和列全部划掉。得到的就是需要求解的位移向量，最后得到的总刚方程为：

武汉理工大学教务处制 45

武汉理工大学教师备课专用纸

?0.5??1??0.5?0???0.51.5?????0.25?0???0?????0.5?0??0?0??0?0.25????0?0????0?v1???3.252??v???1.252?2??????u3??1???0.088??????? ?v3?E??0.372??u5??0.176?????????u0.176???6?按照每个单元的节点组成，从总节点向量位移中挑选出每个单元的节点位移向量，运用前面

单元分析中得到的公式，回代就能计算应变、应力了。（是哪个公式？）

§2.16 四节点矩形单元

前面我们运用三节点三角形单元，获得单元的应变和应力是不变的。这和实际情况有着明显的差异（实际中应力应变在单元中是连续变化的）。产生这种差异的原因，就是我们的位移模式过于简单。下面我们将就如何选择位移模式，提高有限元计算精度，进行说明。同时也可以就如何选择单元，有个初步的认识。首先讲解四节点的矩形单元。 一，矩形单元位移模式

如图所示的矩形单元，四个节点分别为i、j、m、p。

0?0.251.50.25?0.500?0.50.251.5?0.25000.25?0.5?0.251.5?0.50??v1??v?0???2?0???u3????? 0??v3??0.5??u5????0.5????u6??解以上这个方程，可以得到以下的解。

单元节点位移向量????uie?viujvjumvmupvp

?T对应的会有8个节点力分量。由于矩形单元的节点位移向量是8个，那么根据前面我们学过

的三角形推导过程知道，可以建立8个方程，求解8个待定系数。所以在位移模式中的待定系数可以是8个。即：

武汉理工大学教务处制 46

武汉理工大学教师备课专用纸

?u?a1?a2x?a3y?a4xy ?v?a?ax?ay?axy5678?由于我们选择的是多项式插值。所以对于2元函数多项式项元的选择，可以用帕斯卡三角形确定。如左图。

从该三角形可以看出，2元多项式的二次项有三项，即x2、y2、xy。如果选取前面两项中的任何一项，都会造成位移模式偏惠与你所选择的那个方向。从而使位移出现不对称的情况。只有增加xy这一项，才能避免这种情况出现。也就满足了我们材料各向同性的要求。

选取局部坐标系，将四个节点的坐标代入该公式，并求解待定系数，最后可以得出和前面三角形单元类似的关系：

??????Ni?00NiNj00NjNm00NpNm0?ui??ui?0??????????N???? ?Np????v?vp???p?其中的Ni?1?x??y?1?x??y??1???1?? Nj??1???1?? 4?a??b?4?a??b?1?x??y?1?x??y??1???1?? Nm??1???1?? 4?a??b?4?a??b? Np?在上面公式中，代表节点坐标的是什么？

位移模式的讨论。 1，应变和刚体位移 由几何方程得

?x??y??xy??u???a1?a2x?a3y?a4xy??a2?a4y ?x?x?v???a5?a6x?a7y?a8xy??a2?a8x?y?y?v?u?????a5?a6x?a7y?a8xy???a1?a2x?a3y?a4xy??a3?a4x?a6?a8y ?x?y?x?y同三角形单元位移模式一样，常应变和刚体位移包括在位移模式之中。同时我们可以看到，应变中还有随x、y的变量应变。所以在矩形单元不再是常应变或常应力单元了。

2，在局部坐标系中，当x=±a（y=±b）时，位移模式是x（或y）的线性函数（将x=±a（y=±b）代入位移表达式中即可）。这表明在边界上位移按线性规律变化，在公共边界上的位移是连续的，满足收敛性的充分条件。 二，矩形单元的刚度矩阵

??x??????y?????xy?e???e?????x??0???????y?0?0b?y0b?y0??b?y?0??ui????b?y?e????u?1????0??a?x?0??a?x?0a?x0a?x?????????y??v?4ab???a?xb?ya?x??b?y?????a?x???b?y???a?x?b?y??vm?????x??应力向量

???e??D??B???e?。分别代入平面应力和平面应变的弹性矩阵，就可以得出两类问题的应

武汉理工大学教务处制 47

武汉理工大学教师备课专用纸

力解。

刚度矩阵?K??eT相乘后逐项积分即得单元刚度??B??A?D??B?tdxdy代入[B]和[D]矩阵，

矩阵，见教材中的具体计算公式。需要说明的是，由于[B]不再是常数矩阵（还有x、y），所

以积分运算较三角形单元复杂些。具体见如下矩阵。

武汉理工大学教务处制 48

武汉理工大学教师备课专用纸

矩形单元的载荷向节点移置的方法同三角形单元总刚矩阵组集的原则和方法、边界条件的处理等也同三角形单元。其中不同的位置仅是单元节点位移和节点载荷的数目，不再是6个而是8个。 三，讨论

矩形单元的优缺点：

1）其位移模式推导出的应力、应变不再是常量，分布更接近实际物体中的分布。从而使这种单元的计算精度更高。

2）对于斜边性边界和曲线边界的拟合性差，且不便于在不同部位采用不同大小的单元。实际应用中常需要三角形单元在大小不同的矩形单元之间进行过渡。见下图。

§2.17 六节点的三角形单元

从以上对矩形单元的推导，可以看出，增加单元节点的数目可以提高单元位移模式的插值次数，进而可以提高单元计算精度。沿用这样的思路，我们对三角形单元也可以采用相同的办法。不过增加的节点是在每条边的中点。这样就得到了一个六节点的三角形单元。如图所示。

一，位移模式

由于这样的单元有6个节点，所以节点位移向量的个数就是2×6=12个，那么在位移模式中可以有12个待定系数。参看帕斯卡三角形，我们选用以下的位移模式：

?u?a1?a2x?a3y?a4x2?a5xy?a6y2 ?22v?a?ax?ay?ax?axy?ay789101112?将每个单元公共边上的方程代入，可以求得一个二次的抛物线方程，而公共边上有三个公共节点，所以可以唯一确定这条抛物线。这表明该插值模式满足位移连续的收敛性充分条件。

仍然可以按照前面三节点三角形单元的方法，分别将6个节点位移代入，然后解联立方程组，求12个待定系数。这样的计算，由于待定系数过多，计算过程也过于繁杂。所以实际中是采用面积坐标的形式进行计算。 二，面积坐标

三角形三个顶点i、j、m，p为三角形中任意一点，

其在三角形中的位置，可以用Li定。其中：

?LjLm?来确

Li?AjAiA Lj? Lm?m

AAAA——三角形的面积。

Ai——是三角形pjm的面积。 Aj——是三角形pim的面积。

武汉理工大学教务处制 49

武汉理工大学教师备课专用纸

Am——是三角形pij的面积。

当P在图中虚线上任一点是，Li是相同的（三角形面积为底×高的一半），在该三角形的三个顶点，分别有： i，Li=1、Lm=Lj=0 j，Lj=1、Lm=Li=0 m，Lm=1、Li=Lj=0

而且有Li?Lj?Lm?1，将节点的面积坐标和前面我们推导出的形函数Ni进行比较，可以知道，形函数实际就是面积坐标。 三，六节点三角形单元的位移模式

?u?Niui?Njuj?Nmum?N1u1?N2u2?N3u3? v?Nv?Nv?Nv?Nv?Nv?Nviijjmm112233?推求如下：

Ni在i节点处应为1，在j、3节点处应为零，面积坐标Li虽在i节点处为1，但在2、3节点处却为1/2，所以还是用Li作为Ni就不行了。为此我们构造一个函数?Li?Li?该函数是否满足形函数的性质。 由于Li在2、3节点处，Li???1??，考察2?1所以很显然上式等于零。 2??1??=1，可以计算出来） 2?在i节点处要为1，那么β=2。（令?Li?Li?1??? 所以Ni?2Li?Li?2??同理可得Nj?2Lj?Lj???1?1??? Nm?2Lm?Lm??

2?2??对于N1应满足j、m节点处为零，所以构造函数?LjLm， 令?LjLm=1（在节点1处应为1），可以求得β=4 所以N1?4LjLm，同理得N2?4LiLm

N3?4LiLj

综合以上推导，六节点三角形单元的形函数如下：

武汉理工大学教务处制 50

本文来源：https://www.bwwdw.com/article/osxw.html