有限元分析基础教案(武汉理工)

有限元分析基础

第一章有限元法概述

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

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

§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.3 有限元分析的一般过程

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

一，结构离散化

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

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

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

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

二维单元：三角形单元、四边形单元（平面类问题）

三维单元：四面体单元、六面体单元等（空间问题）

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

二，单元分析

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

1，直接法。

2，虚功原理法（变分法）。

3，加权余数法。

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

变分法需要泛函的数学知识，其推导过程具有严谨的数学概念。

加权余数法适用于泛函不存在的应用范围。

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

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

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

由于杆的两个端点节点1、2是单元上的点，所以它们应该满足上述方程。 节点1，x 1=0，∴u 1=a 0+a 1×0=a 0

节点2，x 2=l ，∴u 2=a 0+a 1×l a 1=(u 2-u 1)/l

将求出的结果带入方程并整理，就得：

[][]{}e N u u N N u l x u l x x l u l

u u u δ=??????=+??? ??-=??? ??-+=2121211211

式中：N 1、N 2是形函数

[N]形函数矩阵

{δ}e 节点位移向量

由位移与应变的关系知道：

dx

du dx u du u =-+=ε 将上面推出的位移表达式代入，可得： []{}{}{}[]{}

e e e e B l l l l x l

x dx d N dx d dx du dx u du u δδδδε=??????-=??????-===-+=11 上式中的[B]称为应变矩阵或几何矩阵。

运用材料力学中的虎克定律，可以将应变和应力联系起来。单向应力状态的虎克定律为 []{}[]{}e e S u u l E l E B E E δδεσ=?

???????????-===21 ×× 其中[S]称为应力矩阵。

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

{}[]{}e e e K F δ=

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

[K]e 为单元刚度矩阵。后面将推导出它的计算公式为[][][][]dV B D B K v

T e ?=

[D]矩阵是弹性矩阵。对于一维单元来说，就是E 。

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

[]??????--=??

????-????

??????-=?111111110l EA Adx l l E l l K l e 求得单刚矩阵，也就完成了单元分析。

总结单刚矩阵推导的步骤，应该分为四步：

1） 假定单元内位移变化的近似规律，即选择位移模式。

2） 运用几何关系，推求位移与应变的关系。

3） 应用物理规律，把应变与应力联系起来。

4） 运用虚功方程的力与位移关系，求出单刚矩阵。

单元分析是整个有限元分析的核心。不同的单元因为其力学特性不同，而具有不同的单元刚度矩阵，我们这本教材就是要学习几种常用单元矩阵的推导和计算。了解各种单元的力学特性，为以后选择单元类型打好基础。

三，整体分析

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

[]??????--=1111111l EA K []??

????--=1111222l EA K 整个结构有三个节点，首先将单元刚度矩阵扩充为3X3的矩阵，移动各元素使之与单刚矩阵中的元素位置相对应，如下：

[]????

??????--=000011011111l EA K []????

??????--=110110000221l EA K 然后直接相加。 []?????????????????

?--+--=222222221111110110l EA l EA l EA l EA l EA l EA l EA l EA K 2， 把各单元的节点力向量组集成总的节点载荷向量。

??????????=00P R

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

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

????????????????????????????--+-

-=??????????3212222222211111

10110

00u u u l EA l EA l EA l EA l EA l EA l EA l EA P 修改以后（采用置“0，1”法）

??????????????????

??????????+-

-=??????????3212211111

1100011000u u u l EA l EA l EA l EA l EA P 4， 求解方程组，得出总的节点位移向量。

解得的解是：

????

??????????????+=??????????022*******EA Pl EA Pl EA Pl u u u

有了节点位移，再回代到前面单元推导过程中的公式×和××，就可以求得每个单元的应变和应力了。

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

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

第二章 平面问题

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

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

§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个，它们是：

{}[]T yz xz xy z y x τττσσσσ=

三，位移

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

{}[]T

w v u =δ 符号规定：沿坐标轴正向为正，反之为负。

四，应变

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

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

{}[]T yz xz xy z y x γγγεεεε=

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

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

§2.2 两类平面问题

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

一，平面应力问题

当物体的长度与宽度尺寸，远大于其厚度（高度）

尺寸，并且仅受有沿厚度方向均匀分布的、在长度和

宽度平面内的力作用时，该物体就可以简化为弹性力

学中的平面应力问题。

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

当z=±t/2时，有σz =0、τzx =0、τzy =0。

由于板较薄（相对于长度和宽度尺寸），外力沿板

厚又是均匀分布的，根据应力应连续的假定（弹性力

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

{}[]T xy y x τσσσ=

而且这些应力分量仅是x 、y 两个变量的函数。

二，平面应变问题

当物体是一个很长、很长的柱形体，其横截面

沿长度方向保持不变，物体承受平行于横截面且沿

长度方向均匀分布的力时，该问题就可以简化为平

面应变问题处理。

分析其应力特征。假定其长度方向为无限长，

那么任一横截面都可以看作是物体的对称面，如此

则有该面上的点都有w=0，也就是横截面上的所

有点都不会发生Z 方向的位移。由这一点可以推

出也就有εz =0、τzx =0、τzy =0。

和平面应力相比较，平面应变是εz =0，那么是否也就有σz =0呢？

可能有同学想σ=E ε，当然也就有σz =0，这是错误的。

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

{}[]T xy y x τσσσ=

三，两类问题的比较

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

011111=??+??-????

? ????++??-???? ????+dydx X dx dx dy y dy dy dx x yx yx yx x x x τττσσσ

展开化简得

0=+???

? ????+??? ????X y x yx x τσ 同理可求得∑F y =0满足得条件，0=+???

? ????+???? ????Y x y xy y τσ 由∑M=0，列出方程如下：

021212121=??-?????? ????+-??+?????? ?

???+dy dx dy dx y dx dy dx dy x yx yx yx xy xy xy ττττττ化简后得： dy y

dx x yx yx xy xy ??+=??+ττττ2121 略去微量项，可得：yx xy ττ=。这就是前面所将的剪应力互等。

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

=

+

??

?

?

?

?

?

?

+

?

?

?

?

?

?

?

X

y

x

yx

x

τ

σ

=

+

??

?

?

?

?

?

?

+

??

?

?

?

?

?

?

Y

x

y

xy

y

τ

σ

写成矩阵形式

=

?

?

?

?

?

?

+

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

Y

X

xy

y

x

x

y

y

x

τ

σ

σ

§2.4 几何方程

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

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

仍然取截面的微元体ABCD，AB、CD边长为dx、dy，厚度为“1”。

位移u、v都是x、y的函数，即u(x,y)、v(x,y)，偏导数

x

u

?

?

、

x

v

?

?

表示位移分量u、v沿

坐标轴x的变化率，偏导数

y

u

?

?

、

y

v

?

?

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

dx

x

u

u

u

B?

?

+

=

'dx

x

v

v

v

B?

?

+

=

'

同理的D’点的位移分量

dy

y

u

u

u

D?

?

+

=

'dy

y

v

v

v

D?

?

+

=

'

由于α角在位移和形变很微小的情况下非常小，所以

A’B’≈A’B”

线段AB 位移后的总伸长量为

A ’

B ’-AB=A ’B ”-AB=u B ’-u A =dx x u u ??+

-u=dx x u ?? ∴x

u dx dx x u x ??=??=/ε，同理可得y v dy dy y v y ??=??=/ε 剪应变由α、β两个角度组成 ()x u x v x v u dx x u u dx v dx v B A B B tg ??????+=??????-??++-+==≈1"'"

'αα

由于1πx u ??，所以x

v ??=α，同理可得y u ??=β ∴y u x v xy ??+??=

+=βαγ 综合以上几何方程，并将它们写成矩阵形式：

??????????

??????=??????????????????v u x y y x xy y x 0

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

1） 变形产生的位移。

2） 因运动产生的位移。

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

§2.5 物理方程

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

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

G ——材料剪切弹性模量