大作业报告参考2有限元学习心得

有限元学习心得

吴清鸽车辆工程 50110802411

短短八周的有限元课已经结束。关于有限元，我一直停留在一个很模糊的概念。我知道这是一个各个领域都必须涉及的点，只要有关于CAE分析的，几乎都要涉及有限元。总体来说，这是一门非常重要又有点难度的课程。

有限元方法(finite element method) 或有限元分析(finite element analysis)，是

求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要

基础性原理。将它用于在科学研究中，可成为探究物质客观规律的先进手段。将

它应用于工程技术中，可成为工程设计和分析的可靠工具。本课程教学基本内容

有固体力学和结构力学简介；有限元法基础；桁架、梁、刚架、二维固体、板和

壳、三维固体的有限元法；建模技术；热传导问题的有限元分析；PATRAN软件

的使用.

通过有限元分析课程学习使我了解和掌握了一些有限元知识：

1．简要了解二维和三维固体以及桁架、梁和板结构的三组基本力学方程，即表示位移-应变关系的几何方程，表示应力-应变关系的本构方程和表示内力-外力关系的平衡方程。

2．了解利用能量法形成有限元离散系统方程的基本原理，即哈密尔顿原理。掌握有限元分

析的基本方法及步骤，包括域的离散、位移插值、构造形函数、单元有限元方程

的建立、坐标变换、整体有限元方程的组装、整体有限元方程的求解技术。

3．具体深入的了解并掌握桁架结构、梁结构、刚架结构、二维固体、板和壳结构、三维固体的有限元法分析技术，包括他们具体的形函数构造，应变矩阵，局部坐标系和整体坐标系中的单元矩阵。各种结构的实例研究。

4．了解并掌握建立高质量建模所涉及的各种关键技术。包括单元类型的选择，单元畸形的限制，不同阶数单元混用时网格的协调性问题，对称性的应用（平面对称、轴对称、旋转对称、重复对称），由多点约束方程形成刚域及应用（模拟偏移、不同自由度单元的连接、网格协调性的施加）等，以及多点约束方程的求解。以PATRAN有限元通用软件为例了解一般商业有限元软件的组成及结构。掌握PATRAN软件的基本使用。利用PATRAN软件上机实践完成两个上机练习：刚架结构有限元分析和三维固体有限元分析。

课程的具体学习内容：

内容：

1.三节点三角形单元：单元分析、总刚度矩阵组装、引入约束条件修正总刚度

矩阵、载荷移置、方程求解；

2.四边形单元分析、四节点四面体单元分析、八节点六面体单元分析；

3. 其他常用单元形函数、自由度。

1、三节点三角形单元 1.1. 单元分析

1.1.1 分析步骤

单元分析的任务是建立单元平衡方程，形成单元刚度矩阵。不失一般性，从图1-1三角形离散结构中任取一个单元，设单元编号为e ，单元节点按右手法则顺序编号为 i, j, m,在定义的坐标系xOy 中，节点坐标分别为(xi+yi)，(xj+yj)，(xm+ym)，节点位移和节点力表示如图1-1所示。

取结点位移作基本未知量。由结点位移求结点力：

其中，转换矩阵称为单元刚度矩阵。单元分析的主要目的就是要求出单元刚度矩阵。

1.1.2 位移模式和形函数

对于平面问题，单元任意一点的位移可用位移分量u, v 描述，他们是坐标x, y 的函数。假定三节点单元的位移函数为x, y 的线性函数，六个节点位移只能确定六个多项式的系数，所以平面问题的3结点三角形单元的位移函数如下：

所选用的这个位移函数，将单元内部任一点的位移定为座标的线性函数，位移模式很简单。

位移函数写成矩阵形式为：

{}?????

?????????????????=m m j j i i e

v u v u v u δ{}?????

??????

???????????=m m j j i i e

V U V U V U F {}[]{}

e

e e K F δ=?

?

?

++=++=y a x a a y a x a a u 654321v ??????21a a

将水平位移分量和结点坐写成矩阵： 代入位移函数第一式：

令 则有 A 为三角形单元

[T]的伴随矩阵为 令 则有

同样，将垂直位移分量与结点坐标代入位移插值公式：

最终确定六个待定系数 ：

m

m m j j j i i i y a x a a u y a x a a u y a x a a u 321321321++=++=++=??

???

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

??=??????????321111a a a y x

y x y x u u u m m j j i i

m j i []T 111=????

?

?????m m

j j i i y x y x

y x []1-123i j m a u a T u a u ????

????

=????????????

A

2T =[]T

*T ????

?

????

?---------=i j j

i i j j i m i i m m

i i m j m m j j

m m j x x y y y x y x x x y y y x y x x x y y y x y x ????

?????

?=??

????????=m j

i m j i m j

i m m

m

j j j

i i

i

c c c b b b a a a c b a c b a c b a T

*]T [??

??

?

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

?

?=???

???????m j i m j

i

m j i

m j i

u u u c c c b b b a a a A a a a 21321??

??

?

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

?

??=?????

?????m j i m j

i

m j i m j i v v v c c c b b b a a a A a a a 21654??

??

?

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

?=?????

?????m j i m j

i

m j i

m j i u u u c c c b b b a a a A a a a 21321??

??

?

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

?

??=?????

?????m j i m j

i

m j i m j i v v v c c c b b b a a a A a a a 21654])()()[(21

m m m m j j j j i i i i u y c x b a u y c x b a u y c x b a A

u ++++++++=

])()()[(21

m m m m j j j j i i i i v y c x b a v y c x b a v y c x b a A

v ++++++++=

??i u

令 （下标i ，j ，m 轮换）

[N]称为形态矩阵， N i 称为位移的形态函数

1.1.3 位移函数的收敛性

选择单元位移函数时，应当保证有限元法解答的收敛性，即当网格逐渐加密

时，有限元法的解答应当收敛于问题的正确解答。因此，选用的位移模式应当满足下列两方面的条件：

(1) 必须能反映单元的刚体位移和常量应变。 6个参数 到 反映了三个刚体位移和三个常量应变。 (2) 必须保证相邻单元在公共边界处的位移连续性。 (线性函数的特性)

1.1.4 应变矩阵和应力矩阵

利用几何方程、物理方程，实现用结点位移表示单元的应变和单元的应力。 用结点位移表示单元的应变的表达式为：

[B]矩阵称为几何矩阵

由物理方程，可以得到单元的应力表达式： 为应力矩阵

1.1.5 单元刚度矩阵

)(21

y c x b a A

N i i i i ++={}???

???????

????????????=??????????=m m j j i i m j i e

v u v u v u δδδδ1a

6a

{}????

????

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

??????=???

?????????????+??????=εm m j j i i m m

j j i i

m j i m j i v u v u v u b c b c b c c 0c 0c 0

0b 0b 0b A 21x v y u y v x u e

B }]{[}{δε=[][]

m

j i

B B B B =[]????

?

?????=i i i i i b c c b A B 0

021{}[]{}[][]{}e

B D D δεσ==[][][]

B D S =[][]

m

j i

S S S S =[][][]??

??

?

?

?

???????---==i i i i i i

i i b c c b c b A E B D S 212

1)1(22μμ

μμμ

讨论单元内部的应力与单元的结点力的关系，导出用结点位移表示结点

力的表达式。由应力推算结点力，需要利用平衡方程。用虚功方程表示出平衡方程。

考虑上图三角形单元的实际受力，任意虚设位移，节点位移结点力和内部应力为： 与内部应变为：

令实际受力状态在虚设位移上作虚功，外力虚功为

微小矩形的内力虚功为

根据虚功原理，得

这就是弹性平面问题的虚功方程，实质是外力与应力之间的平衡方程。

虚应变可以由结点虚位移求出： 代入虚功方程

接上式，将应力用结点位移表示出

有

令 {}

{}{}

{} dxdydz ζεF δT *T *???=

*{}

??

???

?????g εε=ε*xy *y *x *m

*m m *m j *j j *j i *i i *i V v U u V v U u V v U u +++++=T []

????

????

?

?????????=m m j j i i *m *

m *j *j *i

*

i V U V U V U v u v u

v u {}{}

e

eT *F δ=dy)(γtdx)(ηdy)(εtdx)(ζdx)(εtdy)(ζdU *xy xy *y y *x x ?+?+?=)tdxdy

ηγζεζ(εxy *xy y *y x *x ++=[]

tdxdy

ηζζ γεεxy y x *

xy *y *x ?????

?????={}{}{}

{}??σε=

δtdxdy

F T

*e

T e *{}

[]{}{}T

T

e

*T e *T

*[B]δ)δB (ε=={}{}{}

{}??=tdxdy

B F T T e

e

T

e

σδδ][**{}{}??=tdxdy

ζ[B]F T e {}[][]{}

e δB D ζ={}{}

??=e

T e δy [D][B]tdxd [B]F []??=y

[D][B]tdxd [B]K T e

则

建立了单元的结点力与结点位移之间的关系， 称为单元刚度矩阵。它

是6*6矩阵，其元素表示该单元的各结点沿坐标方向发生单位位移时引起的结点

力，它决定于该单元的形状、大小、方位和弹性常数，而与单元的位置无关，即不随单元或坐标轴的平行移动而改变。

1.2 总刚度矩阵组装

整体刚度矩阵 是单元刚度矩阵 的集成。 1、刚度集成法的物理概念：

刚度矩阵中的元素是刚度系数，即由单位结点位移引起的结点力。 2、刚度矩阵的集成规则： 先对每个单元求出单元刚度矩阵 ，然后将其中的每个子块 送到结构刚度矩阵中的对应位置上去，进行迭加之后即得出结构刚度矩阵[K]的子块，从而得出结构刚度矩阵[K]。

关键是如何找出 中的子块在[K]中的对应位置。这需要了解单元中的

结点编码与结构中的结点编码之间的对应关系。

结构中的结点编码称为结点的总码，各个单元的三个结点又按逆时针方向编为i,j,m,称为结点的局部码。

单元刚度矩阵中的子块是按结点的局部码排列的，而结构刚度矩阵中的子块是按结点的总码排列的。因此，在单元刚度矩阵中，把结点的局部码换成总码，并把其中的子块按照总码次序重新排列。

以单元②为例，局部码i,j,m 对应于总码5,2,4，因此 中的子块按照总码重新

排列后，得出扩大矩阵

为上图所示： 用同样的方法可得出其他单元的扩大矩阵 将各单元的扩大矩阵迭加，即得出结构刚度矩阵[K]：

集成规则包含搬家和迭加两个环节：

1、将单元刚度矩阵 中的子块搬家，得出单元的扩大刚度矩阵 。

2、将各单元的扩大刚度矩阵 迭加，得出结构刚度矩阵[K]。

{}[]{}

e

e e δK F =[]e

K []

K []

e

k []

e

k []

ij k []

e k

[](2)k [](2)

K [][][](4)(3)(1)K 、K 、

K [][][][][][]∑=+++=(e)(4)(3)(2)(1)K K K K K K []e

k []e

K

1.3 引入约束条件修正总刚度矩阵

整体刚度矩阵[K]求出后，结构的结点力{F}可表示为：

在无支杆的结点处，结点力就等于已知的结点载荷。在有支杆的结点处，则求结点力时，还应把未知的支杆反力考虑在内。如果用{P}表示结点载荷和支杆反力组成的向量，则结点的平衡方程为

根据支承条件对平衡方程加以处理。先考虑结点n 有水平支杆的情况。与结点n 水平方向对应的平衡方程是第2n-1个方程，

根据支承情况，上式应换成 ，即在[K]中，第2n-1行的对角线元素

应改为1，该行全部非对角线元素应改为0。在{P}中，第2n-1个元素应改为0。 此外，为了保持矩阵[K]的对称性，则第2n-1列全部非对角线元素也改为0。

同理，如果结点n 有竖向支杆，则平衡方程的第2n 个方程应改为 ，为此，在矩阵[K]中，第2n 行的对角线元素改为1，该行全部非对角线元素改为0，同时，第2n 列全部非对角线元素也改为0。在{P}中，第2n 个元素改为0。

[]e

K

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

P δ=[K]xn

n 1,2n 2n n 11,2n 2n 11,22n 11,12n P ...v K u K ...v K u K =+++++-----0u n =1

1,2n 2n K --xn

P 0v n =??

?

??????

?

?

??

?

??

???

???????????????????=????????????????????????????????????????????????????????????????????????????00P P P P v u v u v u 00000000000010000000000001000000000000000000y2x2y1x1

n n 2211

1.4 载荷移置 将载荷移置到结点上，必须遵循静力等效的原则。静力等效是指原载荷与结点载荷在任意虚位移上做的虚功相等。在一定的位移模式下，移置结果是唯一的，且总能符合静力等效原则。

单元的虚位移可以用结点的虚位移 表示为 令结点载荷为

集中力的移置

如图所示，在单元内任意一点作用集中力 由虚功相等可得

由于虚位移是任意的，则 体力的移置

令单元所受的均匀分布体力为

由虚功相等可得： 2、四边形、四节点四面体、八节点六面体单元分析 2.1 四边形单元分析

四边单元可取矩形作为研究对象，矩形单元也是一种常用的单元，它采用了比常应变三角形单元次数更高的位移模式，因而可以更好地反映弹性体中的位移状态和应力状态。

矩形单元1234如图所示，其边长分别为2a 和2b ，两边分别平行于x 、y 轴。若取该矩形的四个角点为节点，因每个节点位移有两个分量，所以矩形单元共有8个自由度。采用三角形单元中的方法，同样可以完成对这种单元的力学特性分析。然而，如果我们引入一个局部坐标系ξ、η，那么就可以推出比较简洁

的结果。

在局部坐标系中，节点i 的坐标是(

i , i )，其值分别为±1。取位移模式：

用与三角形单元相似的方法建立形函数，则位移模式可写成：

{}

e

*δ{}{}

e N

f **][δ={}?????

????

?

?

???????????=m m j j i i

e Y X Y X Y X R ???++=tdxdy

p N tds P N P N R T s

T T c e }{][}{][}{][}{{}?

?

?

???=y x P P P ()()}

{][}{}{}{**P N R T

T

e e

T

e δδ=}

{][}{P N R T e =?

??

???=y x p }{()()??=tdxdy

p N R T

T

e

T

e }{][}{}{}{*

*δδ??=tdxdy

p N R T e }{][}{ ) )

) )

x

ξη

αηαξααξηαηαξαα87654321+++=+++=v u {}[]{}∑=?

??

??

?=i i N v u f δ

式中：

由几何方程可以求得单元的应变

可推出：

式中

由虎克定律我们可以得出用节点位移表示的单元应力，即

式中： 对于平面应力问题 有： 若将单元刚度矩阵写成分块形式：

则其中的子矩阵可按下式进行计算：

四边形单元的节点位移与单元节点力之间的关系仍为 ： 其中载荷列阵{R }e 与上节中的(c)式相同，仍可按上式计算等效节点力。但是，需要注意的是，矩形单元有四个节点（1，2，3，4），所以{R }e 具有8个元素，即： 2.2 四节点四面体单元分析

2.2.1 单元划分及位移模式 如图1所示的四面体单元，单元结点的编码为i, j, m, n 。每个结点的位移具有三个分量u , v , w 。这样单元结点的位移列阵可表示成：

[][][]{})4321( , 1001

, ，，，=?

??

???=????

??==i v u I I N N i i i i

i

δ{}?????

?

??

????

??????+=????

??????????????+=??????????????????+=??????????=?ξ??η??η??ξ??ξ??η??η??ξ?????????g εεεv b u a v a u

b ab v a u b v b u a x v y u y v x u xy y x 11111{}[]{}

e

B B B B δε43

21

=[]()()()()???????

???++++=?

??

???

?

??

????????

?=0000111001410

01ηξξηξηηξ?ξ??η??η??ξ?i i i i i i i i i b a a b ab N b N a N a N b ab

B {

}[]{}[]{}

e

S S S S D δεσ4321==[][][]

i i B D S =[]()

()()()()()()??????

?

???????+-+-++++-=00000

02121121111114ηξμξημξηηξμξημηξμi i i i i i i b a a b a b ab E

S []?????

???????=4443

42

41

3433323124232221

14131211

k k k k k k k k k k k k k k k k k [][][][]??=tdxdy

B D B k j

T

i

ij

[]{}{}e e R

k =δ{}[]T e V U V U V U V U R 44332211

=

{}[]

δδδδδe

i j m n i

i i i i i i i i i i

i T

u v w u v w u v w u v w =????????????

??=

单元的位移模式采用线性多项式：

u x y z

v x y z w x y z =+++=+++=+++αααααααααααα123456789101112

式中，为待定系数，由单元结点的位移和坐标决定。将四个结点的坐标（x i , y i , z i ）、（x j , y j , z j ）、（x m , y m , z m ）、（x n , y n , z n ）和结点位移（u i , v i , w i ）、（u j , v j , w j ）、（u m , v m , w m ）、（u n , v n , w n ）代入（2）式可得12个联立方程，解方程组便可求出。将这十二个系数回代到式，则得到由结点位移和形函数表示的单元内任一点的位移表达式：

u N u N u N u N u v N v N v N v N v w N w N w N w N w i i j j m m n n i i j j m m n n i i j j m m n n

=+++=+++=+++

式中

()()

()()N V a b x c y d z N V a b x c y d z

N V a b x c y d z N V a b x c y d z i i i i i j j j j j m m

m m m n n n n n =

+++=-+++=+++=-+++1

61

61

61

6

N i ，N j ，N m ，N n 为四面体单元的形函数

位移模式可以用矩阵形式表示：

{}[

]

{}

[]{}

f u v w N N N N N N N N N N N N N I N I

N I

N I N i j m n i j m n i

j

m

n i j m n i j m n e

e

=????????

?

?=

????????????????????????==0000000000000000

000

δδδδδδ

式中，[I ]为三阶单位阵，[N ]为形函数矩阵。上式即为单元结点位移和单元任意点位移之间的关系。 2.2.2 单元应变和应力

知道单元内任意一点位移后，可利用几何方程确定单元内该点的应变。将位移矩阵式代入空间问题几何方程得：

{}[]{}[

]

{}

εδδ==--B B B B B e i

j

m

n e

其中

上式表明几何矩阵[B ]中的元素都是常量，因此单元中的应变也是常量。也就是说，采用线性位移模式的四面体单元是常应变单元。

将上式代入物理方程，就得到单元的应力列阵：

式中：[S ]为四面体单元的应力矩阵，其分块形式为：

[][][]S D B A V

b A

c A

d A b c A d A b A c d A c A b A d A c A d A b i

i

i

i i i i i i i i i i i i i

i

==

6000

3111111222222

其中

()

()()A A A E 12311221136112=

-=

--=

-+-μ

μ

μμμμμ()

2.2.3 单元刚度矩阵 对于四面体单元，利用虚功原理，采用类似平面问题的处理方法可以得到其单刚矩阵。

{}[][][]{}[]{}R B D B dxdydz K e T e e e

=

=???

δδ

其中：[K ]e

为单元刚度矩阵

[][][][][][][]K B D B dxdydz B D B V

e T T ==??

写成分块形式为：

[]K k k k k k

k k k k k k k k k k k e

ii

ij im in ji

jj jm jn mi mj mm mn ni

nj

nm

nn =--------??????

??

?

??

?

式中子矩阵[ K rs ]由下式计算：

[][]

[][]k B D B V

A V b b A c c d d A b c A c b A b b A d d A c b A b c c c A d d b b A c d A d c A d b A b d A d c A c d d d rs

r

T

s r s r s r s r s r s

r s r s

r s r s

r s r s r s r s r s

r s r s

r s r s

r s ==+++++++++++321212122121212()()A b b c c r s r s 2()+????????

?

?

{ } [ ] { } [ ] [ ] { } [ ] { } [

]

{ }

σ ε δ δ δ = = = = - - D D B S S S S S e e

i

j m n e

可以看出， 单元刚度矩阵是由单元结点的坐标和单元材料的弹性常数所决定的，是一个常数矩阵。 2.3 八节点六面体单元分析 2.3.1 形函数与坐标变换 形函数

坐标变换

2.3.2 位移插值函数与几何矩阵

简记为 ： 2.3.3 单元刚度矩阵与等效节点载荷向量

单元刚度矩阵可以表示为：

将上式中的 替换为 则有：

写成高斯积分形式为：

()()()t t s s r r N i i i i +?+?+=

1118

1

()()()???

?

??

????=?=?=∑∑∑===818

18

1

,,,,,,i i i i i

i i i

i z t s r N z y t s r N y x t s r N x ()()()????

??

??

?

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

?=??????????e e e e N N N N N N N N N x y x w x y x v x y x u 243218882221

11

00000

0000

00000,,,,,,δδδδ

{}[]{}

e

N u δ?=[][][]????

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

?????????????????=??=8882221

11

00000

000000

00000

000

0000

N N N N N N N N N

x z

y z

x y z y

x N B ????????????

????

??[][][][][][][]dxdydz

B D B dv B D B K T

v T

v e

e

e ??=??=??????z y x ,,t s r ,,[][][][]drdsdt

J B D B K T

e

??=???---111111{}[]{}k

j i

k

j

i

k

j

i

b

t

n

i n

j k

j

i

n

k e

b

h

h h t s r J t s r F t s r N P ?=∑∑∑===),,(),,(),,(111

3、其他常用单元形函数、自由度

3.1 轴对称单元

轴对称结构体可以看成由任意一个纵向剖面绕着纵轴旋转一周而形成。此旋转轴即为对称轴，纵向剖面称为子午面，如图表示一圆柱体的子午面abcd 被分割为若干个三角形单元，再经过绕对称轴旋转，圆柱体被离散成若干个三棱圆环单元，各单元之间用圆环形的铰链相连接。对于轴对称问题，采用圆柱坐标较为方便。以弹性体的对称轴为z 轴，其约束及外载荷也都对称于z 轴，因此弹性体内各点的各项应力分量、应变分量和位移分量都与环向坐标θ无关，只是径向坐标r 和轴向坐标z 的函数。也就是说，在任何一个过z 轴的子午面上的位移、应变和应力的分布规律都相同。因此轴对称问题可把三维问题简化为以（z ，r ）为自变量的二维问题。

由于轴对称性，弹性体内各点只可能存在径向位移u 和轴向位移w 。此时，位

移u 、

w 只是r 、z 的函数，而环向位移v =0。即：

轴对称问题的物理方程可写为：

由于轴对称性，我们只需分析任意一个子午面上的位移、应力和应变情况。其有限元分析计算步骤和平面问题相似。首先进行结构区域的有限元剖分。采用的单元是三角形、矩形或任意四边形环绕对称轴z 旋转一周而得到的整圆环，通常采用的单元是三角形截面的整圆环。在单元类型确定之后，单元剖分可以在子午面内进行，如图5-1表示的abcde 子午面被分割为若干个三角形，绕对称轴z 旋转后即形成若干个三棱圆环单元。

相邻的单元由圆环形的铰链相连接。单元的棱边都是圆，故称为结圆。每个结圆与rz 平面的交点称为结点。

),(),(===v z r w w z r u u {}[][]{}εg εεεt σσσσθθD D rz z r rz z r =???

???????????=??????????????=

这样，各单元在子午面rz 平面上形成三角形网格，就如同平面问题中在xy 平面上的网格一样。采用位移法有限元分析，其基本未知量为结点位移。单元的结点位移列阵如下：

对于每一个环形单元，需要假定其位移模式。仿照平面三角形单元，取线

性位移模式：

类似于平面三角形单元的推导，即将单元的结点坐标及结点位移 代入式中,可以解出六个待定系数。再将这些待定系数回代到式中，就可以得到由结点位移和形函数所表示的单元内任一点的位移表达式：

其中形函数：

形函数矩阵的表达式为：

有上面分析可知，轴对称单元自由度有六个。

以上就是关于课程总结的全部内容。通过这八周的学习，我已经对有限元的基础有了一个大致掌握，关于用有限元进行具体分析也掌握了一些最基本的方法。其中应用到很多矩阵变换之中的知识，我会加强这方面知识的巩固。在以后的研究方向重，我也会对有限元分析的方法勤加练习。这门课对我以后的课题方向和分析方法有着举足轻重的作用。感谢雷老师严谨认真的教学，把理论课学习与上机练习紧密结合起来，是我们更加容易掌握要点，更加容易记住方法。在此表示衷心的感谢。

{}[]

[

]

T

m

m j j i i

T

T m T j T i e w u w u w u ==δδδδz

r z r w w z

r z r u u 654321),(),(αααααα++==++==m m j j i i m m j j i i w N w N w N w u N u N u N u ++=++=()()

m j i z c r b a N i i i i ,,2?

++=

[]??

??

??=m j

i

m j i

N N N N N N N 0

000