当前位置:首页 >> 其它课程 >>

有限元总复习


总复习
曹国华

1

1.2有限单元法简介
1. 有限元法
?

分析指导思想
化整为零,集零为整,裁弯取直,以简驭繁,变难为易

即把一个结构看成由若干通过结点相连的单元组
成的整体,先进行单元分析,然后再把这些单元组合 起来代表原来的结构进行整体分析。

1.3有限单元法分析步骤
有限元法分析问题的基本步骤:
1、结构的离散化

离散化就是将要分析的结构分割成有限个单元体,并在
单元的指定位置设置节点,使相邻单元的有关参数具有一定 的连续性,构成单元的集合体代替原来的结构。 结构离散化时,划分的单元大小和数目应根据计算精度 的要求和计算机的容量来决定。

? 选取坐标(右手法则)
? 选择合适的单元,离散结构物为有限个单元,并对单元、节 点进行编号
3

1.3有限单元法分析步骤
2、选择位移模式 为了能用节点位移表示单元体的位移、应变和应力, 在分析连续体问题时,必须对单元中位移的分布做出一定 的假设,一般假定位移是坐标的某种简单函数。

选择适当的位移函数是有限单元法中的关键
{f} — 单元内任意点的位移列矩阵

? f ? ? ? N ??? ?

e

[N ]

— 单元形函数矩阵
— 单元节点位移的列矩阵

?? ?

e

4

1.3有限单元法分析步骤

3、分析单元的力学特性
利用几何方程、本构方程和变分原理建立单元的刚度矩 阵和载荷矩阵

{F}e= [K]e {δ}e
{F}e — 单元节点力 [K]e — 单元刚度矩阵

5

1.3有限单元法分析步骤
4、集合所有单元平衡方程,得到整体结构的平衡方程 先将各个单元刚度矩阵集合成整体刚度矩阵[K]: {K } =Σ{K}e 然后将各单元的等效节点力列阵集合成总的载荷列阵 :

{F } =Σ{F}e
5、由平衡方程求解未知节点位移 {δ} 按照问题的边界条件修改总的平衡方程,并进行求解。

? K ??? ? ? ?F?
[K]整体刚度矩阵, {F }整体载荷 ,{δ}整体节点位移向量
6

1.3有限单元法分析步骤

6、单元应变和应力的计算 根据已知结点的位移利用弹性力学方程和位移插值函数算 出单元的应变和应力。

7

弹性微元

9

一、平衡方程 ?? x ?? xy
?x ?? yx ?x ? ? ?y ?? y ?y ?? zy

二、几何方程
? ? ?? xz ?z ?? yz ?z ?X ?0 ?Y ? 0

?? zx ?? z ? ? ?Z ?0 ?x ?y ?z

?u ?u ?v ? x ? , ? xy ? ? ?x ?y ?x ?v ?v ?w ? y ? , ? yz ? ? ?y ?z ?y ?w ?w ?u ? z ? , ? zx ? ? ?z ?x ?z

三、本构关系(物理方程)

1 2(1 ? ? ) ? x ? [? x ? ? (? y ? ? z )],? xy ? ? xy E E 1 2(1 ? ? ) ? y ? [? y ? ? (? z ? ? x )],? yz ? ? yz E E 1 2(1 ? ? ) ? z ? [? z ? ? (? x ? ? y )],? zx ? ? zx E E
10

对于平面问题

一、平衡方程
? ?? x ?? yx ? ? fx ? 0 ? ?y ? ?x ? ? ?? y ? ?? xy ? f ? 0 y ? ? y ? x ?

三、本构关系(物理方程)
对于平面应力问题

? ?u ? ? ? ?? x ? ? ?x ? ? ? ? ? ? ?v {? } ? ?? y ? ? ? ? ? ? ? ?y ? ? ? ? xy ? ? ?u ?v ? ? ? ? ? ?y ?x ?

二、几何方程

{? } ? [ D]{? }

? ?1 E ? D? ? 2 1? ? ? ?0 ?

?
1 0

? 0 ? 0 ? 1? ? ? ? 2 ?

对于平面应变问题将上式中的E,u置换为:

E ? , 2 1? ? 1? ?

杆单元

14

引 言
?结构离散 杆系的单元离散,各杆有局部坐标系,方向为 杆的轴线方向。结构的刚度方程是在统一的坐标系 (总体坐标系)中建立并求解,因此需将单元局部 坐标各量转换到总体坐标系中。 取杆件与杆件交点、集中力作用点、杆件与支 承的交点为节点。相邻两节点间的杆件段是单元。 节点编号时力求单元两端点号差最小。
15

单元分析
(1)位移模式和形函数 i ·
① 位移模式 px
l

ui

·u

j

x

j

LINK 因为只有 2 个结点,每个结点位移只有 1个自由度,因此 单元的位移模式可设为:

u ? a1 ? a 2 x

(1)

式中a1、a2为待定常数,可由结点位移条件 x=xi 时, u=ui x=xj 时, u=uj 确定。再将由此确定的a1、a2 其代入式(1),得

u ? (u i ?

u j ? ui l

xi ) ?

u j ? ui l

x
16

单元分析 u ? a1 ? a 2 x
u ? (u i ? u j ? ui l
a1

xi ) ?

u j ? ui l
a2

x

? ui ? 1 u ? [( x j ? x ) ? ( xi ? x )] ? ? l ?u j ?
1 N j ] ? [( x j ? x ) ? ( xi ? x )] l
17

形函数 [ N ] ? [ N i

单元分析
② 形函数

将位移函数改写为下列形式

u ? [ N ]{? }
式中形函数[N]为
[ N ] ? [ Ni

e

1 N j ] ? [( x j ? x ) ? ( xi ? x )] l

2 若令系数 ? ? ( x ? xc ), l

xc ?

x j ? xi 2

? 1 ? ? ? ?1

N1 ? 1 (1 ? ? ), 2

N 2 ? 1 (1 ? ? ) 2
18

(2)应变矩阵 一维铰接杆单元仅有轴向应变 e 代入 u ? [ N ]{? } du ?? dx 1

1 ? ? [ ?1 1]{? }e l
上式也可写为

[ N ] ? [( x j ? x ) ? ( xi ? x )] l

? ? [ B]{? }e

式中[B]为应变矩阵
1 [ B] ? [ Bi B j ] ? [?1 1] l
19

单元分析
(3)应力矩阵 由应力应变关系

? ? E?
? ? [ B]{? }
e

? ? E[B]{? }e ? [S ]{? }e
式中[S]为应力矩阵
E [ S ] ? [ ?1 1] l
20

单元分析
(4) 单元刚度矩阵

单元刚度矩阵由虚功原理推出
[k ]e ? ? ?B ? ?D ??B ?dv
T v

对于等截面铰接杆单元(截面积为A ) ,v=Adx,故有:

[k ] ? A? ?B ? ?D ??B ?dx
e T v

EA ? 1 ? 1? [k ] ? ? ? 1 1 l ? ? ?
e

1 [ B ] ? [?1 1] l E [ S ] ? [ ?1 1] l
21

平面桁架杆单元(2D LINK1)
y

?看成局部坐 标下的拉压杆

?2
?1 i j l

?4
?3 x

(1)单元坐标单元位移向量

?? ?

e

?? 1 ? ?? ? ? 2? ?? ? ?? 3 ? ? ?? 4 ? ?

y

?2
?1 i j

?4

?3

x

22

形函数

1 [ N ] ? [( x j ? x ) l

0 ? ( xi ? x ) 0]

1 应变矩阵 [ B ] ? [ Bi B j ] ? [?1 l

0

1

0]

应力矩阵

E [ S ] ? [?1 l

0

1

0]

单元刚度矩阵
?1 ?0 EA ? e [k ] ? l ?? 1 ? ?0 0 ?1 0 0 0 1 0 0 0? 0? ? 0? ? 0?
EA ? 1 ? 1? [k ] ? ? ? 1 1 l ? ? ?
e
23

3、空间杆单元(3D LINK8)
y ?2 ?1 i ?3 l ?6 ?5 j ? 4 x

z

(1)单元坐标单元位移向量

?? ?

e

? ?1 ? 2 ? 3 ? 4

?

?5 ?6 ?

T

24

(2) 形函数
1 [ N ] ? [( x j ? x ) l 0 0 ? ( xi ? x ) 0 0]

(3) 应变矩阵
1 [ B ] ? [ ?1 l 0 0 1 0 0]

(4) 应力矩阵

E [ S ] ? [?1 l (5) 等价节点力
e

0

0

1

0

0]

pl T ?1 0 0 1 0 0? {F } ? 2

25

(6) 单元坐标单元刚度矩阵
对于等截面铰接杆单元,
?1 ?0 ? EA ? 0 e [k ] ? ? l ?? 1 ?0 ? ?0

EA ? 1 ? 1? [k ] ? ? ? 1 1 l ? ? ?
e

0 0 ? 1 0 0? 0 0 0 0 0? ? 0 0 0 0 0? ? 0 0 1 0 0? 0 0 0 0 0? ? 0 0 0 0 0?
26

坐标变换
所讨论的[K]e,{R}e,{ε}, {σ} 都是在单元的局部坐 标下进行的。由于杆件系统在空间中各有自己的方 向,都建立了各自的坐标系,所以在分析中都必须 统一在一个整体坐标中进行,一定要进行坐标变换。
z
z
ui

v j uj
j

x

uj

vi
i O

?
ui
x
27

坐标变换
平面铰接杆单元(桁架元)
z
z

vj
vi ui
i

uj j

x

uj

?
ui
x

O
整体坐标与局部坐标的关系:

ui ? ui cos? ? vi sin ? vi ? ui (? sin ? ) ? vi cos ?
28

坐标变换
z
整体坐标与局部坐标的关系(矩 阵形式): z ui

v jv j

uj
j

x

uj

?ui ? ? cos ? ? ??? ? vi ? ? ? sin ?

sin ? ? ?ui ? ? ? ? cos ? ? ? vi ?

vi v i
i O

?
ui
x

? ui ? ? cos ? ?v ? ? ? i ? ? ? sin ? ?u ? ? ? j? ? 0 ?v j ? ? ? ? ? 0

sin ? cos ? 0 0

0 0 cos ? ? sin ?

0 ? ? ui ? ?v ? ? 0 ? i? ?? ? sin ? ? ?u j ? ?? ? cos ? ? ?v j ?
29

坐标变换
? ui ? ? cos ? ?v ? ? ? i ? ? ? sin ? ?u ? ? ? j? ? 0 ?v j ? ? ? ? ? 0 sin ? cos ? 0 0 0 0 cos ? ? sin ? 0 ? ? ui ? ? ? 0 ? ?vi ? ?? ? sin ? ? ?u j ? ? ? cos ? ? ? ?v j ?
z

?? ? ? T ?? ?

?? ? ? T ?? ?
?1

v j uj

z
vi

0 0 ? ?cos ? ? sin ? ? sin ? cos ? ? 0 0 ? T -1 ? ? 0 cos ? ? sin ? ? ? 0 ? ? 0 0 sin ? cos ? ? ?

ui

j

x uj

?
ui
x
30

i O

在整体坐标系和局部坐标系中节点力之间的 转换关系经推导后,有如下关系式
? Fix ? ? Fix ? ? ? ?F ? F ? iy ? ?1 ? iy ? F ? ? ? ? T ? ? ? T ?1F ? F jx ? ? F jx ? ? ? ? F jy ? ? ? F jy ? ?

z
z ui

vj

uj j

x

uj

vi
i O

?
ui
x
31

在局部坐标系下单元的平衡方程表示为 K? ? F 在整体坐标系下单元的平衡方程为

z

vj

uj j

x
uj

z
vi

ui i

K? ? F
将式转换矩阵代入式,得 O

?
ui

x

KT -1? ? T -1F
得 K = TKT ?1



TKT ?1? ? F

则桁架单元在整体坐标系下的单元刚度矩阵为

K = T KT = T KT
32

?1

T

z

vj

uj j

z
vi

x uj

ui

?
ui

i
O

x 将局部坐标系下的单元刚度矩阵和转换矩阵代入上式, 得平面桁架单元在整体坐标系下的单元刚度矩阵为
? cos 2 ? cos ? sin ? ? cos 2 ? ? cos ? sin ? ? ? ? 2 2 cos ? sin ? sin ? ? cos ? sin ? ? sin ? EA ? ? K ?e ? l ? ? cos 2 ? ? cos ? sin ? cos 2 ? cos ? sin ? ? ? ? 2 2 cos ? sin ? sin ? ? ? ? cos ? sin ? ? sin ?
33

桁架结构的有限元法分析简例
结构离散----单元分析----单元组装 ----形成方程----求解
2 y ② ⑤ 3 P



⑥ ④ 1 l 4 x

l
34



梁单元

1.无轴向变形的平面梁单元 1)单元位移模式 如图所示的梁单元为一个无轴向变形的等截面直杆, 共有两个节点,节点位移包括挠度和转角,节点力包 括剪力和弯矩。
y

Qi i Mi
l

j

Qj

x
Mj

36

单元的节点位移向量表示为
? ? ?vi ?i v j ? j ?
e T

y

节点力向量表示为
F e ? ?Qi M i Q j M j ?
T

Qi i Mi
l

j

Qj

x
Mj

设单元的挠度 v的位移模式取为

v( x) ? ?1 ? ?2 x ? ?3 x2 ? ?4 x3
将单元上两个节点的坐标和相应的位移代入上式,可得单元 的位移模式为 N1 ? (l 3 ? 3lx 2 ? 2 x3 ) / l 3 ?

v( x) ? N e? e

形函数表示为 N e ? ? N1 N2 N3 N4 ?

? N 2 ? (l x ? 2lx ? x ) / l ? ? N 3 ? (3lx 2 ? 2 x3 ) / l 3 ? ? N 4 ? ( x3 ? lx 2 ) / l 2 ?
2 2 3 2

37

2)单元应变 由材料力学可知,梁在发生弯曲变形而引起梁的轴向变 形产生的应变称为梁的弯曲应变,可由下式计算

将位移代入,可得

d 2v ? ? ?y 2 dx

? e ? Be? e
式中: B e——单元应变矩阵,可分块表示为
Be ? ? B1 B2 B3 B4 ?

y ? 12 x ? 6 l ? ?? l3 ? y B2 ? ? 2 ? 6 x ? 4 ? ? ? l ? y B3 ? 3 ?12 x ? 6l ? ? ? l ? y B4 ? ? 2 ? 6 x ? 2l ? ? l ? B1 ? ?
38

3)单元应力 梁的弯曲应力的计算公式为
d 2v ? ? E? ? ? Ey 2 dx

将应变代入上式,可得

? e ? EBe? e ? S e? e
式中: S e——应力矩阵,可分块表示为

S e ? ? S1 S2 S3 S4 ?

Ey ? 12 x ? 6 l ?? 3 ? l ? Ey S2 ? ? 2 ? 6 x ? 4 ? ? ? l ? Ey S3 ? 3 ?12 x ? 6l ? ? ? l ? Ey S 4 ? ? 2 ? 6 x ? 2l ? ? l ? S1 ? ?
39

4)单元刚度矩阵 对于等截面梁单元,单元刚度矩阵为
K e ? ??? B eT DB e dV ? ? B eT DB e Adx
V

? 12 ? EI ? 6l ? 3 l ? ?12 ? ? 6l

?12 6l ? 4l 2 ?6l 2l 2 ? ? ?6l 12 ?6l ? ? 2l 2 ?6l 4l 2 ? 6l

I ——截面对主轴的惯性矩。 式中:

40

2.有轴向变形的平面刚架单元 1)单元位移模式 取节点为i和j之间的杆件为梁单元,在节点i和j上 所受到的节点力为轴力、剪力和弯矩,即节点力向量 为 T e
F ? ? Ni Qi M i N j Q j M j ?

与之相对应的节点位移向量为
? ? ?ui vi ?i u j v j ? j ?
e T

y

Qi i Ni Mi
l

j

Qj Nj x Mj

3)单元的刚度矩阵。由虚功原理得单元刚度矩阵为
? EA ? l ? 12 EI ? 0 l3 ? 6 EI 4 EI ? 0 2 ? l l e eT e k ? ? ? K ? ??? B DB dV ? ? EA ?? 0 0 V ? l 12 EI 6 EI ? 0 ? ? ? l3 l2 ? 6 EI 2 EI ? 0 l2 l ? 对 称 ? ? ? ? ? ? ? ? ? ? ? ? 4 EI ? ? l ?

EA l 0 0

12 EI l3 6 EI ? 2 l

平面三角形单元

3单元分析

-位移模式
F2 F y F1 F3

? 假定三角形单元的位移模式
o

x

y j

vj

?? ?
uj

e

??

?

T i

?

T j

?

T T m

? ? ?u
T

i

vi

uj

vj

um

vm

?

T

v(x,y)

vm m o um

(x,y)

.

vi
u(x,y)

?? i ? ? ?ui
ui

vi ? (i,j,m 轮换) (a)
? ?u ? x, y ? ? ?1 ? ? 2 x ? ? 3 y ? ? ?v ? x, y ? ? ? 4 ? ? 5 x ? ? 6 y

i

?u ? ?f??? ? ?v ?
x

44

3单元分析

-位移模式

根据克莱姆法则,可求出?1, ?2 , ?3

?1 ?
其中

A1 A

?2 ?

A2 A

?3 ?

A3 A

1 xi A ? 1 xj 1 xm

yi y j ? 2? ym

ui A1 ? u j um

xi xj xm

yi yj ym

1 ui A2 ? 1 u j 1 um

yi yj ym

1 xi A3 ? 1 x j 1 xm

ui uj um
45

展开后

3单元分析
A1 ? ui xj xm yj ym

-位移模式
?uj xi xm yi ym ? um xi xj yi yj

? ui ai ? u j a j ? umam 1 1 yj 1 yi A2 ? ?ui ?uj ? um 1 1 ym 1 ym

yi yj

? uibi ? u j bj ? umbm

1 xi A3 ? ui ?uj ? um 1 xj 1 xm 1 xm 1 xj 1 xi

? ui ci ? u j c j ? umcm
其中 ai ?
xj xm yj ym ? x j ym ? xm y j , bi ? (?1) ? y j ? ym 1 ym m
3

1

yj

ci ?

1 xj 1 xm

? xm ? x j

i

j
46

3单元分析 -位移模式
u ? x, y ? ? ?1 ? ?2 x ? ?3 y
1 ? ui ai ? u j a j ? um am ? ? x ? ui bi ? u j b j ? umbm ? ? y ? ui ci ? u j c j ? umcm ? ? 2? 1 ? ? ai ? bi x ? ci y ? ui ? ? a j ? b j x ? c j y ? u j ? ? am ? bm x ? cm y ? um 2?

? ?

?

?

1 令 N i ? x, y ? ? ? ai ? bi x ? ci y ? 2?
可得
u ? Niui ? N j u j ? Nmum ?
i , j ,m

(i, j, m)
i i

?Nu

Ni — 单元的形函数

v ? Ni vi ? N j v j ? Nmvm ?

i , j ,m

?Nv

i i
47

3单元分析
写成矩阵形式

-位移模式
? ui ? ?v ? ? i? 0 ?? ?u j ? ? ? ? Nm ? ? ?vj ? ?um ? ? ? ? ? vm ? ?
e
e

?u ? ? Ni 0 N j 0 N m ?f??? ?? ? 0 N 0 N 0 i j ?v ? ?

所以,单元的位移模式:

? f ? ? ? N ??? ?

[N] — 形函数矩阵
48

3单元分析
写成矩阵形式
?bi 1 ? 0 ?? ? ? ? 2? ?ci ? 0 ci bi

-单元应变
? ui ? ?v ? 0 ?? i ? e ?? ?u j ? ? cm ? ? ? ? ? B ??? ? vj ? ? ? bm ? ?um ? ? ? ? ? vm ? ?
ai ? xj xm yj ym ? x j ym ? xm y j

bj 0 cj

0 cj bj

bm 0 cm

矩阵[B]称为几何矩阵

? B? ? ? ? Bi

Bj

Bm ? ?

?bi 0 ? 其中 1 ? ? (i= i, j, m) B ? 0 c ? i? i? 2? ? ? ?ci bi ? ?

bi ? ( ?1)3 ci ?

1 yj ? y j ? ym 1 ym ? xm ? x j

1 xj 1 xm

因此单元内任一点的应变是节点位移的函数, [B] 是常数, 所以三角形单元是常应变单元。
49

3单元分析
对于平面应力情况

-单元应力
? ? ci ? ? ci ? 1? ? ? bi ? 2 ?

? ? bi ? E ? Si ? ? ? ? bi 2 2 ?1 ? ? ? ? ? 1? ? ci ? ? 2

(i=i, j, m)

对于平面应变情况

E?E

1? ?

2

,? ? ?

1? ?

? bi ? ? ? ? E ?1 ? ? ? bi ? Si ? ? ? 2 ?1 ? ? ??1 ? 2 ? ? ? ? 1 ? ? ? 1 ? 2? ci ? ? 2 ?1 ? ? ? ?

? ci ? 1? ? ? ? ci ? ? 1 ? 2? ? bi ? 2 ?1 ? ? ? ? ?

?

(i=i, j, m)

可知,三角形单元中的应力各处相等
50

4形函数的性质
? 形函数是坐标(x,y)的线性函数

? 形函数在节点上的值
形函数在节点i上的值=1 1 Ni ? xi , yi ? ? ? ai ? bi xi ? ci yi ? ? 1 2? 在其余二节点上的值等于零 1 Ni ( x j , y j ) ? (ai ? bi x j ? ci y j ) ? 0 2? 1 N i ( xm , ym ) ? (ai ? bi xm ? ci ym ) ? 0 2? 所以
?1 Ni ( x j , y j ) ? ? ij ? ? ?0
当 i? j
当 i? j

(i, j, m)

?在单元上任一点的三个形函数之和等于 1
51

4形函数的性质
? 在三角形单元任一边如 i j 边上的形函数,仅与该边的 两端节点坐标有关,与另外一个节点的坐标无关 y m x ? xi N i ( x, y ) ? 1 ? x j ? xi i
x ? xi N j ( x, y ) ? x j ? xi
j n

x o 1 N j ( x, y ) ? (a j ? b j x ? c j y ) yi ? y j 2? y?? ( x ? xi ) xi ? x j 1 b j cm ? bm c j ? [ ( x ? xi )] 2? cm bm

利用形函数的这一性质可以证明,相邻单元的位移分别进 行线性插值之后,在其公共边上将是连续的。
52

x ? xi ? x j ? xi

??

cm

( x ? xi ) ? yi

4形函数的性质
位移函数应该满足的条件: 1、包括常数项 2、包括一次项 3、尽量保证位移的连续性 4、几何各向同性 1、2,有限元解收敛的必要条件,称为完备条件,满足

条件的单元为完备单元;
3,有限元解收敛的充分条件,称为协调条件,满足条件 的单元为协调单元;

协调单元的有限元解一定是收敛的,但非协调单元的解
不一定收敛。
53

1单元平衡方程
[k ] ? ?? A[B] [D][B]tdxdy
T

? B? ? ? ? Bi

Bj

Bm ? ?

?bi 0 ? 1 ? ? B ? 0 c ? i? i? 2? ? ? ?ci bi ? ?

? ? 1 对 ? ? E ? ?D? ? u 1 称 ? 2 1? u ? 1? u ? ?0 0 ? 2 ? ?

其中 ai ?

xj xm

yj ym

? x j ym ? xm y j ,

bi ? (?1)
m i j

3

1

yj

1 ym

? y j ? ym

ci ?

1 xj 1 xm

? xm ? x j

54

1单元平衡方程
[k ] ? ?? A[B] [D][B]tdxdy
T

[k ] ? [B] [D][B]t ?? Adxdy
T

? [B]T[D][B] A ? kii ? ? ? kji ? ?k mi
其中

kij kjj k mj

kim ? ? kjm ? k mm ? ?
55

kii ? 为二阶子矩阵

1单元平衡方程

将几何矩阵和弹性矩阵代入,可得平面应力问题中的每个子 块的矩阵表达式:

K ? Br
e rs

eT

? K rx , sx DB t ? ? ? ? K ry , sx
e s

K rx ,sy ? K ry , sy ? ? 1? ?? ? ? ? br cs ? cr bs ? 2 ? 1? ?? cr cs ? br b s ? ? 2 ?
r, s ? i, j, m

1? ?? ? br bs ? cr cs ? E ?t 2 ? ? 2 1? ?? 4 ?1 ? ? ? ? ? ? ? ? cr bs ? br cs ? ? 2

56

2单元刚度矩阵及其特点
例: 求图3.4所示等腰直角三角形单元的刚度矩阵,设 ? ? 0 ,厚度为t。 解 根据1(0,a)、2(0,0)、3(a,0)的坐标可得
y

1

a x
2

a

3
57

2单元刚度矩阵及其特点
b1 ? y2 ? y3 ? 0 b3 ? y1 ? y2 ? a
2

c1 ? x3 ? x2 ? a c3 ? x2 ? x1 ? 0
a

y

b2 ? y3 ? y1 ? ?a c2 ? x1 ? x3 ? ?a
单元面积 ? ? a 2 ,将以上各数据代入几何矩阵为

1

x
2

a

3

?0 0 ?a 0 a 0 ? ?0 0 ?1 0 1 0 ? 1 ? ? 1 ?0 1 0 ?1 0 0 ? Be ? 2 ? 0 a 0 ? a 0 0 ? a? ? a ? ? ? ?0 0 ?a ?a 0 a ? ? ?0 0 ?1 ?1 0 1 ? ?

由于? ? 0 ,则平面应力与平面应变问题的弹性矩阵 D 相同,为 ?
?1 ? E? ? D? ?? 1 2 ? (1 ? ? ? ) ? 0 0 ? ? 0 ? ?1 0 0 ? ?0 1 0 ? 0 ? ? E ? ? ? ? 1? ?? ? ?0 0 0.5? ? 2 ? ?

58

2单元刚度矩阵及其特点
应力转换矩阵,如下 从以上求得的结果可以 看出,单元刚度矩阵是对 称的。 1 、 3 、 5 列相加等 于零, 2 、 4 、 6 列相加也 等于零,单元刚度矩阵是 奇异的。另外在几何矩阵 中所涉及到的元素均为 “坐标差量”,而与坐标 的具体值无关,从而也说 明了单元刚度矩阵不随单 元或坐标轴的平行移动而 改变。

?0 0 ?2 0 2 0 ? E ? ? e e S ? DB ? 0 2 0 2 0 0 ? 2a ? ? ?1 0 ?1 ?1 0 1 ? ?
则单元刚度矩阵为
? 1 0 ?1 ?1 0 1 ? ? 0 2 0 ?2 0 0 ? ? ? 3 1 ?2 ?1? Et ? ?1 0 e eT e K ? B DB t ? ? ? ? 3 0 ?1? 4 ? ?1 ?2 1 ? 0 0 ?2 0 2 0? ? ? 1 0 ? 1 ? 1 0 1 ? ?

59

2单元刚度矩阵及其特点
y

1

单元2

单元3

a
单元1
2

x a
3

2单元刚度矩阵及其特点
单元刚度矩阵具有如下特点:

1.K e 中的每一个元素都是一个刚度系数
假设单元的节点位移为? 到节点力为
e

? ?1 0 0 0 0 0?,由 F e
T

? K e? e ,得

? Fix ? ? K ix ,ix ? ?F ? ?K ? ? iy ? ? iy ,ix ? ? ? Fjx ? ? ? ? K jx ,ix ? ? e F ?? ??? ? F K ? jy ? ? jy ,ix ? ? Fmx ? ? K mx ,ix ? ? ? ? ? ? ? Fmy ? ? ? ? K my ,ix ? ?

61

2单元刚度矩阵及其特点
K ix,ix 表示i节点在水平方向产生单位位移时,在 由此可得:

节点i的水平方向上需要施加的节点力; K iy,ix 表示i节点在水平 方向产生单位位移时,在节点i的垂直方向上需要施加的节点

力。选择不同的单元节点位移,可以得到单元刚度矩阵中每
个元素的物理含义。 因此, K e 中的每一个元素都是一个刚度系数,表示单位节点位 节点位移矢量对单元某个节点力矢量的贡献率。

e K 移分量所引起的节点力分量, 中的每一个子块表示单元某个

62

2单元刚度矩阵及其特点
2. K e是对称矩阵
K e 是对称矩阵。根据对称性,可以减少计算量和存储量。
e 3. K e 是奇异矩阵,即 K ? 0

根据这一性质,当已知单元节点位移时,可以求出单元 节点力。反之,由于单元刚阵奇异,不存在逆矩阵,因此, 当已知单元节点力时不能求出单元上的节点位移。
从物体变形的实际情况来说,单元刚阵奇异是必须的, 因为单元除了产生变形外,还会产生刚体位移,即仅仅依靠 节点力是无法唯一地确定刚体位移的。事实上,当单元的节 点力为零时,单元仍可做刚体运动。
63

2单元刚度矩阵及其特点
e 例如,假定单元产生了x方向的单位刚体位移即 ? ? ?1 0 1 0 1 0?

T

并假设此时对应的单元节点力为零,则由

F e ? K e? e 得

?0 ? ?1 ? ?0 ? ?0 ? ? ? ? ? ?0 ? ? ? e 1 ? ?? K ? ? ?0 ? ?0 ? ?0 ? ?1 ? ? ? ? ? ?0 ? ?0 ?
可以得到,在单元刚度矩阵中1,3,5列中对应行的系数相 K e ? 0。 加为零,由行列式的性质可知,

64

2单元刚度矩阵及其特点

4. 单元的刚度不随单元或坐标轴的平行移动而改变

单元的刚度取决于单元的形状、大小、方向和弹性系数,
而与单元的位置无关,即不随单元或坐标轴的平行移动而改 变。因此,只要单元的形状、大小、方向和弹性系数相同, 无论单元出现在任何位置均有相同的单元刚度矩阵,根据对 称性,可以减少计算量。

65

1 整体刚度矩阵
将各单元刚度方程右边相加,从而得到整体刚度矩阵,如下
K?K
(1)

?K

(2)

?K

(3)

?K

(4)

y

? K11 K12 ?K ? 21 K 22 ? ? K 31 K 32 ? ? K 41 K 42 ? ? K 51 K 52 (1) (2) ? K11 ? K11 ? (1) ? K 21 (1) (2) ? ? K 31 ? K 31 ? (2) ? K 41 ? 0 ?

K13 K 23 K 33 K 43 K 53

K14 K 24 K 34 K 44 K 54

K15 ? K 25 ? ? K 35 ? ? K 45 ? K 55 ? ?
(1) (2) K13 ? K13 (1) (3) K 23 ? K 23 (1) (2) (3) (4) K 33 ? K 33 ? K 33 ? K 33 (2) (4) K 43 ? K 43 (3) (4) K 53 ? K 53

4 ④ ② ① 1 2 3③

5

x

o
(2) K14 0 (2) (4) K 34 ? K 34 (2) (4) K 44 ? K 44 (4) K 54

(1) K12 (1) (3) K 22 ? K 22 (1) (3) K 32 ? K 32 0 (3) K 52

? 0 ? (3) K 25 ? (3) (4) ? K 35 ? K 35 ? (4) K 45 ? (3) (4) ? K 55 ? K 55 ?

66

1 整体刚度矩阵
整体刚度矩阵组装的基本步骤: 1)将单元刚度矩阵中的每个子块放到在整体刚度矩阵中的对应 位置上,得到单元的扩大刚度矩阵。注意对于单元刚度矩阵是 按照局部编码排列的,即对应单元刚度矩阵中的i、j、m;对于 整体刚度矩阵是按照整体编码排列的,即按节点号码以从小到 大的顺序排列。在组装过程中,必须知道单元节点的局部编码 与该节点在整体结构中的整体编码之间的关系,才能得到单元 刚度矩阵中的每个子块在整体刚度矩阵中的位置。将单元刚度 矩阵中的每个子块按总体编码顺序重新排列后,可以得到单元 的扩大矩阵。例如在图中,单元②的局部编码为i、j、m,对应 整体编码为1、3、4,然后将单元②刚度矩阵中的每个子块按总 体编码顺序重新排列后,可以得到单元的扩大矩阵。 2)将全部单元的扩大矩阵相加得到整体刚度矩阵。
67

1 整体刚度矩阵
通过以上组装过程可以得到组装整体刚度矩阵的一般规则: 1)结构中的等效节点力是相关单元结点力的叠加,整体 刚度矩阵的子矩阵是相关单元的单元刚度矩阵子矩阵的集成; 2)当整体刚度矩阵中的子矩阵 K rs中r=s时,该节点(节点r 或s)被哪几个单元所共有,则K rs 就是这几个单元的刚度矩阵 e 中的子矩阵 K rs 的相加。如 K 33 应该是单元①-④中对应子矩阵 (1) (2) (3) (4) 的集成,即 K33 ? K33 ? K33 ? K33 ? K33
y
4 ④ ② ① 1 2 3③ 5

x
68

o

1 整体刚度矩阵
3)当 K rs中 r ? s 时,若rs边是组合体的内边,则K rs就是共用该 e 边的两相邻单元刚度矩阵中的子矩阵K rs 的相加。如13边为单元 (1) (2) ①和②的共用边,则 K13 = K13 ? K13 4)当K 中 =0。如节点r=1和 rs r和s不同属于任何单元时,则 K rs s=5不同属于任何单元,此时K15 =0。
y
4 ④ ② ① 1 2 3③ 5

x
69

o

1 整体刚度矩阵
例: 如图所示有限元模型,弹性模量为 简化计算取 ? ? 0,求整体刚度矩阵。
y
1
a

,厚度为 t ,为 E

2
a

① ③ ②

3 ④ 5 6 a

x

4 a

70

2 整体刚度矩阵的特点
用有限元方法分析复杂工程问题时,节点的数目比较多,

整体刚度矩阵的阶数通常也是很高的。那么在进行计算时,
如果存储整体刚度矩阵的全部元素,将会浪费较大的资源、 降低计算效率。如果根据整体刚度矩阵的特点进行编写程序, 可以大大节省资源、并提高计算效率。因此有必要了解和掌 握整体刚度矩阵的特点,整体刚度矩阵具有以下几个显著的

特点:
1. K 是对称矩阵

2. K 中主对角元素总是正的 3. K 是稀疏矩阵,非零元素呈带状分布
4. K 是奇异矩阵,在排除刚体位移后,它是正定阵
71

2 整体刚度矩阵的特点
1. K 是对称矩阵
由单元刚度矩阵的对称性和整体刚度矩阵的组装过程, 可知整体刚度矩阵必为对称矩阵。利用对称性,在计算机 编写程序时,只存储整体刚度矩阵上三角或下三角部分即 可。 2. K 中主对角元素总是正的 例如,刚度矩阵 K 中的元素K 33表示节点2在x方向产生单位 位移,而其它位移均为零时,在节点2的x方向上必须施加的 力; K 44 表示节点2在y方向产生单位位移,而其它位移均为零 时,在节点2的y方向上必须施加的力。很显然在此情况下力 的方向应该与位移方向一致,故应为正号。

72

2 整体刚度矩阵的特点
3. K 是稀疏矩阵,非零元素呈带状分布 如果遵守一定的节点编号规则,就可使矩阵的非零元素都 集中在主对角线附近呈带状。整体刚度矩阵中的子矩阵 K rs 只 有当下标s等于r或者s与r同属于一个单元时才不为零,这就说 明,在第r双行中非零子矩阵的块数,应该等于节点r周围直接 相邻的节点数目加1。可见, K 中元素一般都不是填满的,是稀 疏矩阵,且非零元素呈带状分布。 以下图所示的单元网格为例,其整体刚度矩阵中的非零子 块(每个子块为2行2列)的分布情况如下图所示,图中阴影 部分表示该子块不为零,其它子块部位均为零。

73

2 整体刚度矩阵的特点
显然,带状刚度矩阵的带宽取决于单元网格 中相邻节点号码的最大差值 D 。把半个斜带形区 域中各行所具有的非零元素的最大个数叫做刚度 矩阵的半带宽(包括主对角元),用B表示,如下 B=2(D+1) 通常的有限元程序,一般都利用刚度矩阵的对 称和稀疏带状的特点,在计算求解中,只存储上 半带的元素,即所谓的半带存储。因此,在划分 完有限元网格进行节点编号时,要采用合理的编 码方式,使同一单元中相邻两节点的号码差尽可 能小,以便节省存储空间、提高计算效率。
74

2 整体刚度矩阵的特点
4. K 是奇异矩阵,在排除刚体位移后,它是正定阵 无约束的弹性体(或结构物)的整体刚度矩阵是奇异的, 不存在逆矩阵,即关于位移的解不唯一。这是因为弹性体在 外力的作用下处于平衡,外力的分量应该满足三个静力平衡 方程。这反映在整体刚度矩阵中就意味着存在三个线性相关 的行或列,所以是个奇异阵,不存在逆矩阵。例如:设弹性 体在外力的作用下处于平衡,这时相应的解为 ? 1 ,然后在给 ?1 ? ? 2 仍是问 予弹性体以刚体位移而相应的节点位移 ? 2 ,这时, 题的解,因为刚体位移不会破坏平衡。 注:当排除刚体位移后,整体刚度矩阵是正定矩阵。

75

3 边界条件
1.代入法 该方法保持方程组仍为2n×2n系统,仅对整体刚度矩阵 和整体载荷列阵进行修正。下面以一个只有四个方程的简 单例子加以说明,方程如下
? K11 ?K ? 21 ? K 31 ? ? K 41 K12 K 22 K 32 K 42 K13 K 23 K 33 K 43 K14 ? ? u1 ? ? F1 ? ? v ? ?F ? K 24 ? ? 2? 1? ?? ? ? ? ? ? K 34 ? ?u2 ? ? F3 ? ? K 44 ? ? ? v2 ? ? ? ? F4 ? ?

假定系统中节点位移 u1 ? c1 、u2 ? c2 ,则当引入这些节 点的已知位移之后,方程就变成

76

3 边界条件
? K11 ?K ? 21 ? K 31 ? ? K 41 K12 K 22 K 32 K 42 K13 K 23 K 33 K 43 K14 ? ? u1 ? ? F1 ? ? v ? ?F ? K 24 ? ? 2? 1? ?? ? ? ? ? ? K 34 ? ?u2 ? ? F3 ? ? K 44 ? ? ? v2 ? ? ? ? F4 ? ? ?1 0 ?0 K 22 ? ?0 0 ? ?0 K 42 c1 0 ? ? u1 ? ? ? ? v ? ?F ? K c ? K c ? 0 K 24 ? ? 2 21 1 23 2 ? 1? ?? ? ??? ? c2 1 0 ? ?u2 ? ? ? ?? ? ? 0 K 44 ? ? v2 ? ? F4 ? K 41c1 ? K 43c2 ? ? 0

若 c1 ? c2 ? 0

,则

?1 0 ?0 K 22 ? ?0 0 ? ? 0 K 42

0 ? ? u1 ? ? 0 ? ? v ? ?F ? 0 K 24 ? ? 2? 1? ?? ? ? ? ? ? 1 0 ? ?u2 ? ? 0 ? ?? ? ? ? 0 K 44 ? ? v2 ? ? F4 ? 0

然后,用这组维数不变的方程来求解所有的节点位移。显 然,其解答仍为原方程的解答。 在手算时,可直接将零位移约束所对应的整体刚度矩阵中的 行和列直接划去,使得整体刚阵的维数变小,更便于手算。
77

3 边界条件
2.乘大数法 将 K 中与指定的节点位移相对应的主对角元素乘上一个 大数,同时将 F 中的对应元素换成指定的节点位移值、该大 数与节点位移相对应的主对角元素三者的乘积。若把此方法 用于上面的例子,则方程就变成
? K11 ?1015 ? ? K 21 ? K 31 ? ? ? K 41 K12 K 22 K 32 K 42 K13 K 23 K 33 ? 1015 K 43 K14 ? ? u1 ? ? c1 K11 ? 1015 ? ?? ? ? ? K 24 ? ? v1 ? ? F2 ? ? ??? ? ? K 34 ?u2 ? ?c2 K 33 ? 1015 ? ?? ? ? ? K 44 ? F4 ? ? ? v2 ? ?

该方程组的第一个方程为 K11 ?1015 u1 ? K12v1 ? K13u2 ? K14v2 ? c1K11 ?1015 解得 u1 ? c1 ,这种方法就是使 中相应行的修正项远大于非修正项。
78

5 收敛准则
对于有限元这种数值计算方法,一般总是希望随着网 格的逐步细分所得到的解能够收敛于问题的精确解。根据前 面的分析,可知在有限元分析中,一旦确定了单元的形状之 后,位移模式的选择将是非常关键的。由于载荷的移置、应 力矩阵和刚度矩阵的建立等等,都依赖于单元的位移模式, 所以,如果所选择的位移模式与真实的位移分布有很大的差 别,那么就很难获得良好的数值解。为了保证解答的收敛性, 要求位移模式必须满足以下三个条件,即 1)位移模式必须包含单元的刚体位移。

2)位移模式必须包含单元的常应变。
3)位移模式在单元内要连续、且在相邻单元之间的位移必 须协调。
79

5 收敛准则
? 经验证明,根据巴斯卡( Pascal )三角形来选
择二维多项式的各项。在二维多项式中,如果包

含有对称轴一边的某一项,那么就必须同时包含
有另一边的对称项。选择多项式位移模式时,还

应该要考虑的一个因素是,多项式中的项数必须
等于或稍大于单元边界上的外节点的自由度数。

通常是取项数与单元的外节点的自由度数相等,
取过多的项数是不恰当的。
80

5 收敛准则

C x x x x x
5 4 3 2

常数项

y xy y xy
2 2 2 2 2

线性项 二次项

x y
3

2

y xy
3 3

3

三次项

x y
4

x y
3 2

y xy
4

4

四次项

x y

x y

x y

y

5

五次项

??
对称轴
81

平面问题有限元分析 矩形单元-等参单元
曹国华

82

6.1四节点矩形单元位移函数
这里引入一个局部坐标系?、?,这样可以推出比较简洁 的结果。如图所示,取矩形单元的形心o为局部坐标系的原 点,?和?轴分别与整体坐标轴x和y平行,两坐标系存在有以 下的坐标变换关系

x ? xo? ? a?

y ? yo? ? b?
式中:xo? 、yo? ——矩形形心处坐标。 矩形形心处坐标以及矩形长、宽可由下式计算
x0? ? ( x1 ? x2 ) / 2 ? ( x3 ? x4 ) / 2 ? ? y0? ? ( y2 ? y3 ) / 2 ? ( y1 ? y4 ) / 2 ? ? a ? ( x2 ? x1 ) / 2 ? ( x3 ? x4 ) / 2 ? ? b ? ( y3 ? y2 ) / 2 ? ( y4 ? y1 ) / 2 ?
y 4 3

?
2b

o?

?

1

2a

2

o

x

83

四节点矩形单元位移函数
y 4 3

?
4(-1,1) 3(1,1)
?
2b

o?

?

o
2 x

?

1

2a

o

1(-1,-1)

2(1,-1)

在局部坐标系中,节点i的坐标是 (?i ,?i ) ,其值分别 为±1。如节点1在局部坐标系下的坐标为(-1,-1)。

84

四节点矩形单元位移函数
由于矩形有4个节点,共8个自由度,可以选择有8个待 定参数的位移模式,如下 常数项

C

x

y xy y2 xy 2 x2 y3 y3 xy 3 y4 xy 4 y5

线性项 二次项 三次项 四次项 五次项

u ? ?1 ? ? 2? ? ?3? ? ? 4?? v ? ?5 ? ? 6? ? ? 7? ? ?8??
x5 x4

x2 x3 x3 y x4 y x2 y x3 y 2

x2 y 2 ??
对称轴

该函数称为双线性函数。将节点的局部坐标值代入上 式,可列出四个节点处的位移分量,即两组四元联立方程 ,由此可求得位移模式中的8个未知参数?1,?2,…,?8
85

四节点矩形单元位移函数
?u1 ? ?v ? ? 1? ?u2 ? ? ? 0 ? ?v2 ? ? ? N4 ? ? ?u3 ? ?v3 ? ? ? ?u4 ? ?v ? ? 4?

?
4(-1,1) 3(1,1)

?u ? ? N1 e u ?? ??? ?v ? ? 0
4

0 N1

N2 0

0 N2
4

N3 0

0 N3

N4 0

o
1(-1,-1)

?

u ? ? Ni ?? ,? ? ui , v ? ? Ni ?? ,? ? vi
i ?1 i ?1

2(1,-1)

f ? N e? e
式中: Ni ——矩形单元的形函数,i=1,2,3,4;
e N —— 形函数矩阵;

? i ? ?ui ? ——单元节点位移列阵,
e

vi ? ,i=1,2,3,4。
T

86

四节点矩形单元位移函数
? N1 0 N 2 N ?? ? 0 N1 0
e

0 N2

N3 0

0 N3

N4 0

0? e ? ? N ? 1 N4 ? ?

e N2

e N3

e ? N4 ?

?1 0 ? N ? Ni I = Ni ? ? 0 1 ? ?
e i

(i=1,2,3,4)

形函数的表达式为
1 ? N1 ? ?1 ? ? ??1 ? ? ? ? 4 ? 1 N 2 ? ?1 ? ? ??1 ? ? ? ? ? 4 ? 1 N3 ? ?1 ? ? ??1 ? ? ? ? ? 4 ? 1 N 4 ? ?1 ? ? ??1 ? ? ? ? 4 ?
4(-1,1)

?
3(1,1)

o
1(-1,-1)

?

2(1,-1)

87

四节点矩形单元应变与应力矩阵
x ? xo? ? a? y ? yo? ? b?
? ? Ni ? ? ?x ? e Bi ? ? 0 ? ? ? Ni ? ??y ? ? 1 ? Ni 0 ? ? ? ? a ?? ? Ni ? ? ??? 0 ?y ? ? ? Ni ? ? 1 ? Ni ? ?x ? ? ? b ?? ? 0 ? ? 0 ? b?i ?1 ? ?i? ? ? 1 ? Ni ? 1 ? ? ? 0 a ? 1 ? ? ? ? ? i i ? ? b ?? ? 4 ab ? ? a?i ?1 ? ?i? ? b?i ?1 ? ?i? ? ? ? ? 1 ? Ni ? a ?? ? ?

1 N i ? ?1 ? ?0 ??1 ? ?0 ? 4

?0 ? ?i?
?0 ? ?? i
88

四节点矩形单元刚度矩阵
矩形单元刚度矩阵的推导过程与三节点三角形单元类似, 即 K e ? ?? ? B eT DB e dV,由前文可知 K e的推导过程与形函数的具体 表达形式、节点个数均无关,该表达式具有普遍意义。 若单元厚度t 为常量,则K 可以进一步表示为
K e ? t ?? B eT DB e dxdy
e

将单元刚度矩阵写成子块的形式,如下
e ? K11 ? e K 21 e ? K ? e ? K 31 ? e ? ? K 41 e K12 e K 22 e K 32 e K 42 e K13 e K 23 e K 33 e K 43 e ? K14 e ? K 24 ? e ? K 34 e ? K 44 ? ?

89

四节点矩形单元刚度矩阵
? K1 E ?t K ? ? 4ab ?1 ? ? ?2 ? ? K 3
e rs

K2 ? K4 ? ?

式中:
? ?r? s ? 1 ? ? ? 2 ? ? r? s ? ? K1 ? b ? r? s ?1 ? ? a ? ? r s ?1 ? ? ? 3 2 3 ?? ? ? ? ? ? 1? ?? ? ? K 2 ? ab ? ? ?? r? s ? ? r? s ? ? 2 ? ? ? ? ? 1 ? ? ? ? ? (r、s=1,2,3,4) K 3 ? ab ? ? ?? s?r ? ? s? r ? ? 2 ? ? ? ? ? ? ? 1? ?? 2 ? ?? ? K 4 ? a 2?r? s ?1 ? r s ? ? b ? r? s ? 1 ? r s ? ? 3 ? 2 3 ?? ? ? ?
2

90

坐标变换与等参元的概念
? x ? ? N1 e x ?? ??? ? y? ? 0 0 N1 N2 0 0 N2 N3 0 0 N3 N4 0

? x1 ? ?y ? ? 1? y ? x2 ? ? ? 0 ? ? y2 ? ? ? N4 ? ? ? x3 ? ? y3 ? ? ? ? x4 ? ?y ? 0 ? 4?

整体坐标
( x2 , y2 )

?
? ?1

? ? ?1

4 1

3

( x1 , y1 )

o
? ? ?1

?

? ?1
( x4 , y4 )

( x3 , y3 )

2

x

x ? ? N i ?? ,? ? xi , y ? ? N i ?? ,? ? yi
i ?1 i ?1

4

4

局 部 坐 标 自 然 坐 标

?
4(-1,1) 3(1,1)

o
1(-1,-1)

?

/

2(1,-1)

等参单元-任意直边四边形单元(4节点)
x ? ? Ni ?? ,? ? xi , y ? ? Ni ?? ,? ? yi
i ?1 i ?1 4 4

式(Δ)是形状函数,与位移函数一样,在i节点,Ni=1,在 其它节点,Ni=0,该形状函数与位移函数一样求解。
为简便起见,将(Δ)式写成如下通用公式 1 N i ? (1 ? ? i? )(1 ? ?i? ) (i ? 1,2,3,4) 4
92

1 ? (1 ? ? )(1 ? ? ) ? 4 ? 1 N 2 ? (1 ? ? )(1 ? ? ) ? ? 4 ? 1 N 3 ? (1 ? ? )(1 ? ? ) ? ? 4 ? 1 N 4 ? (1 ? ? )(1 ? ? ) ? ? 4 N1 ?

(?)

等参单元-任意直边四边形单元(4节点)
位移函数:

u ? ?1 ? ? 2? ? ?3? ? ? 4??
?u1 ? ?v ? ? 1? ?u2 ? ? ? N ( ? , ? ) 0 N ( ? , ? ) 0 N ( ? , ? ) 0 N ( ? , ? ) 0 u ? ?v2 ? ? ? ? 1 2 3 4 e u ?? ??? ? ?u ? 0 N ( ? , ? ) 0 N ( ? , ? ) 0 N ( ? , ? ) 0 N ( ? , ? ) v ? ? ? 1 2 3 4 ?? 3? ?v3 ? ? ? ?u4 ? ?v ? 4 4 ? 4?

v ? ?5 ? ? 6? ? ? 7? ? ?8??

u ? ? Ni ?? ,? ? ui , v ? ? Ni ?? ,? ? vi
i ?1 i ?1

93

等参单元-任意直边四边形单元(4节点)
u ? ? Ni ?? ,? ? ui , v ? ? Ni ?? ,? ? vi x ? ? Ni ?? ,? ? xi , y ? ? Ni ?? ,? ? yi
i ?1 i ?1 i ?1 4 i ?1 4 4 4

对比

节点位移个数与节点坐标个数一样,在选取单元的位移函 数时,可取单元的自然坐标作为自变量,取图形(坐标)变换 式的形状函数为位移函数的形状函数。 坐标变换函数与位移函数采用相同的形状函数
94

1 ? (1 ? ? )(1 ? ? ) ? 4 ? 1 N 2 ? (1 ? ? )(1 ? ? ) ? ? 4 ? 1 N 3 ? (1 ? ? )(1 ? ? ) ? ? 4 ? 1 N 4 ? (1 ? ? )(1 ? ? ) ? ? 4 N1 ?

(?)

等参单元刚度(4节点)
根据形函数,几何矩阵为
? ? Ni ? ? ?x ? ? 0 ? ? ? Ni ??y ? ? 0 ? ? ? Ni ? ? y ? (i=1,2,…,4) ? ? Ni ? ?x ? ?

Bie ?

N1 ? N2 ? N3 ? N4 ?

1 ? (1 ? ? )(1 ? ? ) ? 4 ? 1 (1 ? ? )(1 ? ? ) ? ? 4 ? 1 (1 ? ? )(1 ? ? ) ? ? 4 ? 1 (1 ? ? )(1 ? ? ) ? ? 4
95

等参单元刚度(4节点)
N1 N2 N3 N4 1 ? ? (1 ? ? )(1 ? ? ) ? 4 ? Ni ? Ni ? 1 4 ? (1 ? ? )(1 ? ? ) ? ?y ?x x ? ? N i ?? ,? ? xi ? 4 ? i ?1 1 4 ? (1 ? ? )(1 ? ? ) ? ? 4 y ? ? N i ?? ,? ? yi ? i ?1 1 ? (1 ? ? )(1 ? ? ) ? ? 4 N i 是ξ、η的函数,因此必须用坐标变换式来转 在上式中,

换导数关系,根据复合函数的求导法则,有 ? ? N i ? ?? x ? y ? ?? N ? ?? Ni ? i ? ?? ? ? ? ?x ? ? ? ?? ? ? ? ?? ?? ? x ? ? ? ? ? ?? ? ? ?? N ? ? J ? ? (i=1,2,…,4) ? N ? x ? y i ? N ? ? ? i ? ?? i ? ? ? ? ?? ? ? ? ??y ? ? ? ?? ?? ? ?? ? ??y ? ?
96

等参单元刚度(4节点)
? ? Ni ? ? ? x ? ?? ? ? ?? ? ? ? ? ??? ? ? Ni ? ? ? x ? ? ?? ? ?? ? ? ?

? y ? ? ? Ni ? ? ? Ni ? ?? ? x ? ? ?x ? ? ?? ? ? ? ? ? ?? N ? ? J ?? N ? ? y ?? i ? i? ? ? ?? ? ??y ? ? ??y ? ? ??

? Ni ?x

? Ni ?y

式中的J称为 雅可比(Jacobi)矩阵

?? x ? ?? ? J=? ?? x ? ? ??

? y? ?? ? ? ? ? y? ?? ? ?

? ? Ni ? ? ? Ni ? ? ? ? ? ?? ? ?x ? ? ?1 ? ?? N ? ? J ? ? ? i? ? ? Ni ? ? ? ??y ? ? ? ?? ? ?

关键

求出

J

?1

97

等参单元刚度(4节点)

将单元应变代入平面问题的物理方程式,就得到平面4 节点等参单元的应力列阵,为

Sie ? DBie (i=1,2,…,4)
等参元刚度矩阵仍可用下式表示

K e ? ?? ? B eT DB e dV ? t ? ? B eT DB e dxdy
K ? t?
e 1 ?1 ?1

?

1

BeT DBe J d? d?

98

等参单元刚度(4节点)
K ? t?
e 1 ?1 ?1

?

1

BeT DBe J d? d?

由此式可计算出单元的刚度矩阵,此处的积分一般得不到 显式形式,因此要采用数值积分法(如Gauss积分法),该积 分对任何形状的四边形单元具有相同的积分表达式和相同的积 分限。上式积分运算是在母单元(自然坐标)内进行的,但节 点力与节点位移的方向均沿整体坐标x,y方向。所以[k]e是整 体坐标内的刚度矩阵。 由以上分析可知,在划分单元时,只需确定单元节点的整 体坐标值,而不必画出等参元的具体形状,因为在计算中实际 使用的只有单元四个节点在整体坐标系下的位移值。 等参元变换的条件为 J ? 0 ,因此在有限元网格划分时,要 特别注意这一点。
99

高斯积分法
被积函数f一般很复杂,往往不能得出它的显式, 在有限元计算中常采用数值积分。在单元内选出某些 积分点(2~10),算出被积函数f在这些积分点的值, 然后用加权系数乘上这些积分值,再求出总和作为近 似的积分值

高斯积分法用同样数目的积分点具有较高精度
一维高斯积分公式

?

1

?1

f (? )d? ?

? W f (? )
i i i ?1

n

式中,n为积分的点数,?i为积分点i的局部坐标 ,Wi 为加权系数 ,f(?i)为被积函数在积分点处的函数值

?

1

?1

f (? )d? ?

? W f (? )
i i i ?1
f

n

f (? )

f1

?1

1 ?

I ? ? f (? )d? ? 2 f1
?1
101

1

f (? ) ? c0 ? c1? ? c2? ? c3?
2

3

f

f (? )

c2 I ? ? f (? )d ? ? 2(c0 ? ) ?1 3
1

f (?1)

f (?2)

?

1

?1

f (? )d? ?

?W f (? )
i ?1 i i

2

?1 ?1

?2 1 ?

c2 Wi f (?i ) ? 2(c0 ? ) ? 3 i ?1
102

2

f (? ) ? c0 ? c1? ? c2? ? c3?
2

f
3

f (? )

c2 Wi f (?i ) ? 2(c0 ? ) ? 3 i ?1
2 3 3 1 2 3

2

f (?1)

f (?2)

c2 W1 (c0 ? c1?1 ? c2?1 ? c ? ) ? W2 (c0 ? c1?2 ? c2?2 ? c3?2 ) ? 2(c0 ? ) 3
?W1 ? W2 ? 2 ?W ? ? W ? ? 0 1 1 2 2 ? ? 2 ? 2 2 ?W1?1 ? W2?2 ? 3 ? 3 3 W ? ? W ? ? ? 1 1 2 2 ?0

?1 ?1

?2 1 ?

?W1 ? W2 ? 1 ? 1 ? ??1 ? ?2 ? ? 3 ? ?0.57735 ?
103

高斯积分法
积分点数的选取与被积函数f(?)有关,当f(?)是m次 多项式时,取 m ?1 n? 2 当f(?)不是多项式时,需要通过一些试算来判断选取 适当的n值。积分点数n(2~10)不能取得过大,否则将 随着积分点的增多计算量急剧增加

??

1

1

?1 ?1

f (? ,? )d?d? ?

?? W W
i
n

n

m

j

f (? i ,? j )
j k

i ?1 j ?1
m

???

1

1

1

?1 ?1 ?1

f (? ,?, ? )d?d?d? ?

??? W W W
i i ?1 j ?1 k ?1

q

f (? i ,? j , ? k )

Gauss积分
高斯求积公式中的积分点坐标和权系数P126表6-1 n 1 2 3 4
±? i
Hi

0 +/- 0.5773502692 +/- 0.7745966692 0.0000000000

2 1 0.5555555556 0.8888888889

+/- 0.8611363116 +/- 0.3399810436

0.3478548451 0.6521451549

热传导问题的有限元法
泛函与变分

变分法预备定理P133 7-5

结论P135 7-9

设函数F(x)在[x1, x2]连续,对于δy(x),如果有 则 F ?x ? ? 0, ?x1 , x2 ? 。 δy(x)是y的变分。

?

x2

x1

F ?x ??ydx ? 0

设泛函

I ?y?x?? ? ? F ?x, y?x?, y' ?x??dx
x2 x1
x2
1

泛函取极值的必要条件 ?I ? ?x 即

? ?F d ? ?F ?? ? ? ? ? ??ydx ? 0 ? ? ? ?y dx ? ?y ' ??

?F d ? ?F ? ? ? ? ?0 ? ? ?y dx ? ?y' ?

第8章 场(流体流动)问题的有限元法
1. 流场问题及加权余量法 2. 二维流体流动的有限元格式

3. 流场有限元分析的几个问题

二 加权余量法

加权余量法的基本思想:通过使试探函数与真值
的加权误差在求解域内的总和为零,以求得满足

微分方程的近似解。
解释:如果已知求解问题的微分方程,假定一个 近似解,把近似解代入到微分方程中,这种运算 必定产生误差项(余数)。该余数定义为R,虽 然不能使这些误差(余数)为0,但可以让该余 数项的加权积分为0。

由于试函数的不同,余量可有如下三种情况,依此加 权余量法可分为: 1.内部法 试函数满足边界条件,此时消除余量的条件成为:

? W R dV ? 0
V Ii I

(i ? 1, 2,? , n)

(8.1.6)

2.边界法 试函数满足控制方程,此时消除余量的条件为:

?W
S

Bi

RB dS ? 0

(i ? 1, 2,? , n)

(8.1.7)

3.混合法 试函数不满足控制方程和边界条件,此时用式 (8.1.5)来消除余量。

8.1.2

基本方法概述

下面以内部法为例,介绍按权函数分类时加权余量的2种 基本方法。对内部法来说,消除余量的统一格式是:

? W R dV ? 0
V Ii I

(i ? 1, 2,? , n)

1.子域法(Subdomain Method) 2. 配点法(Collocation Method) 3.最小二乘法(Least Square Method) 4.伽辽金法(Galerkin Method) 5.矩法(Method of Moment)

3.最小二乘法(Least Square Method) 本法通过使在整个求解域上余量的平方和取极小来建

立消除余量的条件。
若记余量平方和为I(C),即

I (C ) ? ? RI2 dV ? ? RIT RI dV
V V

?RI T ?I (C ) 则极值条件为: ?C ? 2?V ( ?C ) RI dV ? 0

?RI 由此可见,本法权函数为: WIi ? ?C i

(i ? 1, 2,?, n)

4.伽辽金法(Galerkin Method)
本法是使余量与每一个基函数正交,也即以基

函数作为权函数

WIi ? Ni

(i ? 1,2,?, n)

当试函数 包含整个完备函数集时,用本法必 可求得精确解。

? 4. 例题解析
例题1
求解二阶常微分方程:

d 2u ?u ? x ? 0 2 dx
边界条件: 当x=0时,u=0 当x=1时,u=0

?0 ? x ? 1?



② ③

取它的近似解为

u ? x?1 ? x ??a1 ? a2 x ? ??



? 4. 例题解析
例题1

u ? x?1 ? x ??a1 ? a2 x ? ??
其中ai为待定参数;

N1 ? x ?1? x ? , N2 ? x ?1? x ? x....... 为试探函数
显然近似解满足边界条件,但是不满足微分方程,所
以会产生余量。余量的加权积分为零:
j

W Rdx ? 0 j ?
0

为方便起见,我们只讨论一项和两项的近似解: 一项近似解,n=1:

u1 ? a1x ?1 ? x ?
代入,得余量为:



R1 ?x? ? x ? a1 ? 2 ? x ? x2
两项近似解,n=2:

?

?



u2 ? x ?1? x ?? a1 ? a2 x ?
余量为:



R2 ?x? ? x ? a1 ? 2 ? x ? x2 ? a2 2 ? 6x ? x2 ? x3

?

? ?

?



——最小二乘法
将余量的二次方 R 2 在域中积分:

I ? ? R 2 d?
选择近似解的待定系数ai,使余量在全域的积分值达到极 小。为此必须有: 对ai求导,得到:
?

?I ?0 ?a i

i ? 1,2,?, n i ? 1,2,?, n

?R R d? ? 0 ? ?ai ?

由此得到n个方程,由此求解n个待定系数ai,将上式与 式④比较可得,最小二乘法的权函数选择为:

?R Wi ? ?ai

i ? 1,2,?, n

一项近似解:

R1 ?x? ? x ? a1 ? 2 ? x ? x2
?R1 ? ?2 ? x ? x 2 ?a1

?

?

代入求导式得到:
1 1

解得: a1 ? 0.2723 一项近似解为:

?R1 2 2 R dx ? x ? a ? 2 ? x ? x ? 2 ? x ? x dx ? 0 1 1 ? ? ?a1 0 0

?

?

???

?

u1 ? 0.2723x?1 ? x?

两项近似解:
W1 ?

R2 ?x? ? x ? a1 ? 2 ? x ? x2 ? a2 2 ? 6x ? x 2 ? x3
?R2 ? ?2 ? x ? x 2 W ? ?R2 ? 2 ? 6 x ? x 2 ? x3 2 ?a1 ?a
2

?

? ?

?

代入求导式得到两个方程:
1 1

?R2 2 2 3 2 R dx ? x ? a ? 2 ? x ? x ? a 2 ? 6 x ? x ? x ? 2 ? x ? x dx ? 0 2 1 2 ? ? ?a1 0 0

?

?

? ?

???

?

?R2 2 2 3 2 3 R dx ? x ? a ? 2 ? x ? x ? a 2 ? 6 x ? x ? x 2 6 x ? x ? x dx ? 0 2 1 2 ? ? ?a2 0 0

1

1

?

?

? ?

???

?

解得: a1 ? 0.1875

a2 ? 0.1695

两项近似解为: u ? x?1? x??0.1875? 0.1695x? 2

——伽辽金法

取近似函数作为权函数
一项近似解: 取权函数:
1

u1 ? N1a1 ? a1x ?1 ? x ?

W1 ? N1 ? x?1 ? x ?
1 2 1

? N R ?x ?dx ? ? x?1 ? x ??x ? a ?? 2 ? x ? x ??dx ? 0
1 1 0 0

解得:

a1 ? 5 / 18

5 一项近似解: u1 ? x ?1 ? x ? 18

两项近似解:
W1 ?

R2 ?x? ? x ? a1 ? 2 ? x ? x2 ? a2 2 ? 6x ? x 2 ? x3
?R2 ? ?2 ? x ? x 2 W ? ?R2 ? 2 ? 6 x ? x 2 ? x3 2 ?a1 ?a
2

?

? ?

?

代入求导式得到两个方程:
1 1

?R2 2 2 3 2 R dx ? x ? a ? 2 ? x ? x ? a 2 ? 6 x ? x ? x ? 2 ? x ? x dx ? 0 2 1 2 ? ? ?a1 0 0

?

?

? ?

???

?

?R2 2 2 3 2 3 R dx ? x ? a ? 2 ? x ? x ? a 2 ? 6 x ? x ? x 2 6 x ? x ? x dx ? 0 2 1 2 ? ? ?a2 0 0

1

1

?

?

? ?

???

?

解得: a1 ? 0.1875

a2 ? 0.1695

两项近似解为: u ? x?1? x??0.1875? 0.1695x? 2

第九章 薄板弯曲问题的有限元法
1 薄板弯曲问题的力学描述 2 矩形单元薄板弯曲问题有限元方程

3 三角形单元薄板弯曲问题有限元方程

3 三角形单元薄板弯曲问题有限元方程
一.整体坐标系下的位移插值函数
单元节点位移用矩阵表示为

?? ? ? ? A??? ?
e
可惜在此情况下,对于二个边界

?? ? ? ? A? ??
?1

e

?

y 3(0,c)

分别平行于x轴和y轴的等腰三角
形单元,矩阵[A]是奇异的.

1(0,0)

2(c,0)

x

?1 0 ?0 0 ? ?0 ?1 ? ?1 c ? A? ? ?0 0 ? ?0 ?1 ?1 0 ? ?0 0 ?0 ?1 ?

0 1 0 0 1 c 1 0

0 0 0 c 0
2

0 0 0 0 c 0 0 0 ?c

0 0 0 0 0 0 c2 2c 0

0 0 0 c 0
3

0 0 0 0 c2 0 0 0 ?c
2

0 ?2c 0 0 0

?3c 2 0 0 0

0 ? ? 0 ? 0 ? ? 0 ? 0 ? ? 0 ? c3 ? ? 2 3c ? ? 0 ?

在实际意义中,避免这种情况发生,通常采用面积坐标的方法。

1、面积坐标
三角形单元ijm中,任一 点P(x,y)的位置,可以用 如下的三个比值来确定:

Ai Li ? A

Am Lj ? Lm ? A A

Aj

我们以原三角形边所对的角码来命名此三个子三角形的面积,即

?Pjm

面积为Ai;

? Pmi

面积为Aj;

P点的位置可由三个比值来唯一确定

?Pij 面积为Am P ? Li , L j , Lm ?

5.基于面积坐标的插值函数 三角形单元的插值函数

1. 六结点三角形单元计算公式
3

(0,0,1)

(0,0,1)

1? ?1 , 0, ? ? 2? ?2

1

(1,0,0)

(1,0,0)

? 1 1? ? 0, , ? ? 2 2?

2

(0,1,0)

?1 1 ? ? , ,0? ?2 2 ?

(0,1,0)

对于三个边中点,可选择非此边中点所在的其他二条边作为直线方程,则可得

N4=4L1L2

N5=4L2L3

N6=4L1L3

1. 十结点三角形单元计算公式
2? ?1 , 0, ? ? 3? ?3 1? ?2 , 0, ? ? 3? ?3
1(1,0,0)
?2 1 ? ? , ,0? ?3 3 ? ?1 2 ? ? , ,0? ?3 3 ?
? 2 1? 6 ? 0, , ? ? 3 3?

3(0,0,1)

? 1 2? 7 ? 0, , ? ? 3 3? ?1 1 1? 10 ? , , ? ?3 3 3?

2(0,1,0)

对于角节点

Ni ?

Li ? 2

Li ? 1 L 1 3? 3 ? i ? ? 3L ? 1?? 3L ? 2 ? L i i i 1 2 1 2 3 3

习题1:用试凑法求出如图所示10节点三角形单元每个
结点的形函数(用面积坐标L1、L2、L3表示),1-3边长

为L,作用三角载荷,最大值为q,将三角载荷表示成L1
的线性形式,并求9节点等效力。
3
8 9

a !b ! a b L1 L3ds ? L ? (a ? b ? 1)! 1?3
9 解: N 9 ? L1 L3 ? 3L1 ? 1? 2
1

7
10

6

4

5

2

qx ? L1q

F9 ?

1? 3

? N q ds ? .......
9 x

第二章 平面问题的基本理论
§2-1 弹性力学基本概念 §2-2 弹性力学的基本方程

§2-3 平面应变和平面应力问题
§2-4 能量原理与变分法

弹性理论问题需要解一系列偏微分方程组,并满足 边界条件,这在数学上往往遇到困难。因此需要寻求近似

的解法。变分法的近似解法是常用的一种方法。在数学上, 变分问题是求泛函的极限问题。在弹性力学里,泛函就是

弹性问题中的能量(功),变分法是求能量(功)的极值, 在求极值时得到弹性问题的解,变分问题的直接法使我们

比较方便地得到近似解。

显然,实际存在的位移,除了满足位移边
界条件以外,还应当满足位移表示的平衡方程 和应力边界条件; 现在又看到,实际存在的位移,除了满足 位移边界条件外,还满足位移变分方程。而且, 通过运算,还可以从位移变分方程导出用位移 表示的平衡微分方程和应力边界条件。于是可 见:位移变分方程可以代替平衡微分方程和应 力边界条件。

例题1 试用里兹法求如图所示梁的位移和固定端弯矩。 (用Ritz法、最小势能原理、虚位移原理求解)

q
A
EI , l

F
B
x

注: 对一维杆(梁)系结构:

? W? ? ? ? ( x)? ( x)dx
0

l

? ? W ? ? ? M ?? dx
e

M 为弯矩,

??

虚曲率

小挠度情形下

1

M(x)


?(x)

EI

2v d k ( x) ? =+ ?( x) - d x 2

1

此即弹性曲线的小挠度微分方程

M ( x) ? EIv??
? We ? ? M ?? dx ? ? EIv??? v??dx

非常感谢同学们一 个学期以来对我教 学工作的支持!

预祝同学们期末考 得好成绩,实现自 己的人生目标!


相关文章:
有限元期末复习
有限元复习 有限元考试主要围绕平面力系、杆件结构、等参元出题,分简答,小计算,计算三大题 老师没说具体出几道题,让大家好好看看印的那本书 一、简答、填空:...
有限元复习题及答案.pdf
有限元复习题及答案.pdf_从业资格考试_资格考试/认证_教育专区。有限元课程习题 1、试简要阐述有限元分析的 基 本步骤主要有哪些。有限元分析的主要步骤主要 有:...
有限元复习问题
有限元复习问题_理学_高等教育_教育专区。1.有限元程序设计的基本原理是什么? ...(即总的平衡方程),包括将刚度集成总刚,以及将单元的等效结点力列阵集成总的...
有限元复习重点总结
有限元复习重点总结_理学_高等教育_教育专区。1. 有限元法采用:加强余量法(或...70. 总刚的性质:1)对称性;2)奇异性;3)稀蔬;4)非零元素呈带状分布;5)主...
有限元总复习_图文
有限元总复习_理学_高等教育_教育专区。适合北交岩土专业硕士学生,内容详细 解决力学和物理问题的方法 :1 解析方法 精确解;2 数值分析方法 近似解; 直接求解基本...
有限元复习重点
有限元复习重点_信息与通信_工程科技_专业资料。●有限元起源于 20 世纪 50 ...单元贡献矩阵 总刚度矩阵为 刚度矩阵后,再形成载荷列阵,即可得整体刚度方程,经...
华中科技大学硕士 有限元 复习
华中科技大学硕士 有限元 复习_工学_高等教育_教育专区。有限元分析有限...总体刚度矩阵组装原则及总刚阵 特点 18、 固有频率与特征向量(振型)定义 及理解...
有限元法复习提纲
有限元复习提纲 有限元复习提纲第一章 绪论 有限元法实质 1. 有限元法...3)整体刚度矩阵的特点 ) a)K 是对称矩阵 b)K 中的主对角元素总是正的 c...
有限元 总复习
有限元复习宝典 11页 免费有​限​元​ ​ ​总​复​习 ...2. ANSYS求解模块是程序用来完成对已经生成的有限元模型进行分析和求解。在此...
有限元分析复习资料打印版
有限元分析复习资料打印版_理学_高等教育_教育专区。有限元复习资料 1.简述有限...当采用位移法时,物体或 结构离散化之后,就可把单元总的一些物理量如位移,应变...
更多相关标签: