第一节三角形常应变单元(DOC)

第三章平面问题的有限元法本章通过三角形常应变单元，介绍有限元法应用于弹性体应力分析的基本原理和方法：包括弹性体的离散化，单元特性的分析，刚度矩阵的建立，等效节点力的计算，解答的收敛性以及实施步骤和注意事项，同时对形函数的性质也作简要的叙述。

第一节三角形常应变单元

一、结构离散化

用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。

(a) (b)

图3.1 弹性体和有限元模型

将整个结构（平板）划分成有限个三角形。这样的三角形称为单元（三角形单元）。

三角形单元的顶点取为节点。3节点三角形单元用边界节点间的直线段来近似板的曲线边界。

这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。

注：1. 全部节点和全部单元一般由1开始按自然顺序编号。2. 节点编码：

总码-----------用于整体分析，如1,2，…，按自然顺序编号局部码--------用于单元分析，i，j，m 要求按逆时针方向的顺序进行编码

每个单元的节点局部码i，j，m和节点总码有一一对应的关系

3. 单元间不能有重叠

4. 一个单元的任一顶点不许为另一单元任一边的内点

5. 所有作用在单元上的载荷，包括集中载荷、表面载荷和体积力，都按虚功等效的原则移置到节点上，成为等效节点载荷。

二、 位移模式

1. 单元节点位移列阵

i

u

图 3.2 平面三角形单元

设单元e 的节点号码为i ，j ，m 。由弹性力学平面可知，单元内任意一点有两个位移分量u ，v ，记为

{}T

f u v ????=

故每个节点也有两个位移分量，因此称节点自由度为2。3个节点得位移分量分别是 ,,,,,m m i i j j u v u v u v ，用列阵表示为

{}

i e

i i e j j j m m m u v u v u v δδδδ??

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

??????

== （3-1）

称单元自由度是6。其中任一子矩阵为

{}

T

i i

i u v δ??????

= （i ，j ，m 轮换）

2. 位移模式

结构承受载荷以后发生位移，但位移分布事先并不知道。位移法有限元采用节点位移为基本未知量。因此，在应用位移法有限元时，需要对单元内部的位移分布进行假定，使其满足节点的位移连续条件和单元边界的位移连续条件。

单元位移分布的假定一般选用代数多项式，多项式的系数待定。这是一种近似方法。代数多项式的次数可以选择很高，不过次数越高，分析越麻烦，但精确度越好。

这种假定的单元位移分布形式称为位移模式，它是坐标x 和y 的函数，所以也称为位移函数。

对于3节点三角形单元，选用的位移模式是把单元中任一点的位移u ，v 表示为坐标x 和y 的线性函数，即

123546u x y v x y αααααα???

??

=++=++ （3-2）

式中512346,,,,,αααααα为待定常数

设各节点坐标为（x i ,y i ）,(x j ,y j ),(x m ,y m )，同时设各节点位移为（u i ,v i ）,(u j ,v j ),(u m ,v m )代入式（3-2）得

i i i y x u 321ααα++=

i i i y x v 654ααα++= j j j y x u 321ααα++=

j j j y x v 654ααα++=

m m m y x u 321ααα++= m m m y x u 654ααα++=

由上式左边的三个方程可以求得

112i i i j j j m m m u x y u x y u x y α=?，211121i i

j j m m

u y u y u y α=?，311121i i j j m m x y x y x x α=?

其中

1211i i

j j m m

x y x y x y ?=

式中?为三角形面积，为了保证求得的面积为正值，三个节点i ，j ，m 必须按逆时针编排，如图3-2所示。

将321,,ααα代入式（3-2），经整理得

1[()()()]2m m m m i i i i j j j j u a b x c y u a b x c y u a b x c y u =++++++++?

其中

m m i j j m

i j m i j a x y x y b y y c x x ???

??

??

=-=-=- （i ，j ，m 轮换） （3-3）

同理得

1[()()()]2m m m m i i i i j j j j v a b x c y v a b x c y v a b x c y v =++++++++?

若令

1

()2i i i

i N a b x c y =

++? （i,j,m 轮换） （3-4） 则得位移模式为

m m i i j j m m i i j j u N u N u N u v N v N v N v ??

?

??

=++=++ （3-5）

也可写成矩阵形式

{}

{}00000

i i e

m i

j j i

j

j j m m u v N N N u u f N v N N N v u v δ????????????????????

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

???

?????????

=== （3-6）

式中,,m i j N N N 是坐标的线性函数，它们反映了单元的位移状态，所以称为形函数，且称

0000

00m i j m i j N

N N N N N N ??

????????

?

?

=

为形函数矩阵 。其中 {}

T

e

m m i i j j u v u v u v δ

?????

?

= 例题1 图示单元，已知各 节点的坐标（单位: m ）， 计算 ：1. 形函数的表达 式；13边中点A 的形函数。

2. 已知各节点的位移：

1（0，-0.001），2（0.002,0），

3（0,0），计算13边中点A 的位移。 图3.3 例题1

y

解：1. Δ=1

1202022011m m i j j m i j m i j a x y x y b y y c x x =-=?-==-=-=-=-=-=-

00202020m m j i i m j i m j i a x y x y b y y c x x =-=-?==-=-==-= 001

m i j j i m i j m j i a x y x y b y y c x x =-==-==-=

因 1()2i i i i N a b x c y =++?

（i,j,m 轮换）， 得

1(,)(22)

2

1(,)22

1(,)2i j m N x y x y N x y x x N x y y

=--=?==

在13边中点A 有x=0, y=1 ，将其代入上式，的

1

(,)2

(,)01

(,)2

i j m N x y N x y N x y === 2. 单元节点位移

{}

000.00200

0.001T

T

e

m m i i j j

u v u v u v δ??

??

???

?

?

?

==-

由方程（3-5），得13边中点A 的位移

110000()

22

1100(0.001)0.0005()

22

m m i i j j m m i i j j u N u N u N u mm v N v N v N v mm =++=?++?==++=?++?-=-

三、 应变

有了单元位移模式（3-5），利用平面问题的几何方程

{}

x y xy u x v y u v y x εεεγ????

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

????????

???==???+?? 可以求得应变分量。

j m i x m

i j j m i y m

i j j j m m i i xy m m

i j i j N N u N u u u x x x x

N N v N v u v y y y y

N N N N u v N N u u u v v v y x y y y x x x

εεγ????==++????????==++????????????=+=+++++???????? 而

,22i i i i

N b N c x y ??==

??

?? （i ，j ，m 轮换） 所以

1()

21()21()

2x m m i i j j y m m i i j j xy m m m m i i j j i i j j bu b u b u c v c v c v c u c u c u b v b v b v εεγ=++?

=++?

=+++++?

写成矩阵形式

{}00010002i i m

i j e

j m i j j m

m i

i j j m m u v b b b u c c c v c b c b c b u v ε????????

???????????

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

????

=

? 简写成

{}

{}e

e

B εδ????

= (3-7)

将其写成分块矩阵形式

,,m i j B B B B ??????????

= (3-8)

而子矩阵

0102i i i i

i b

B c c b ??

??

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

=? （i ，j ，m 轮换） (3-9) 注：1. 式(3-9)是用节点位移表示单元应变的矩阵方程，其中矩阵[]B 称之为单元应变矩阵。

2. 由于,,,,,,m m i j i j b b b c c c ?等都是常数，所以矩阵[]B 中的元素都是常数，因而单元中各点的应变分量,,x y xy εεγ也都是常数。故这种单元称为常应变单元。

例题2 对于例1单元，试计算单元应变。 解：

{}

{}{}x e y m m i i j j xy B B B εεεδδδγ????????????????????????

?

?????????

??

??==++ {}{}

0i δ= 所以单元应变为

{}{}{}x e y m

m

j j

xy B B εεεδδγ??????????????????????

??==+

20000.00201100012200.00102100.0020.0005000110022j m j m m j m j m m j j b b u u c c v v c b c b ????

??

??????

??????

????

??

???????

??

??????????

?

??

??

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

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

???

?????

?????????

=+-=-=+??

四、 应力

求得应变后，利用物理方程

{}{}

D σε????= 便可导出以节点位移表示的应力关系式中。把式（3-7）代入上式，得

{}

{}e

e

D B σδ????????

= （3-10）

令

S D B ?

?????????

??= 则 {}

{}e

e

S σδ????

= （3-11）

上式表示的是应力与节点位移之间的关系。式中矩阵[]S 称之为单元应力矩阵，写成分块矩阵的形式

m m i

j i j S D B

B B S S S ??????????????

????

?

?

== （3-12） 对于平面应力问题，其弹性矩阵[]D 为

210101100

2E D μ

μμμ?

?

??

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

?

=--

代入式（3-12），得对应于平面应力问题的应力矩阵

22(1)1122i i i i i i i i b c E S D B b c c b μμμμμ??

??

??

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

?

?????

?

==-?--（i,j, m 轮换） (3-13)

对于平面应变问题，只要把平面应力问题的弹性矩阵中的E 换成

2

1E μ-，μ换成

1μ

μ

-，即得其弹性矩阵为

[]

10

1(1)

10(1)(12)11200

2(1)E D μ

μμμμμμ

μμ?

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

?

则对应于平面应变问题的应力矩阵为

[][][]1(1)2(1)(12)112122(1)2(1)i i i i i i i i b c E S D B b c c b μμμμμμμμμμμ????-??

-??

==??+-?-??

--????

--??

（i,j, m 轮换） (3-14)

注：1. 由(3-13) 和(3-14)式知，[]S 中的元素都是常数，所以每个单元中的应力分量也是常数。故这种单元也称为常应力单元。 2. 由于相邻单元将具有不同的应力和应变。这样越过公共边界，即从一个单元过渡到另一个与它相邻的单元，应力和应变的值都将有突变；但从一个单元过渡到另一个与它相邻的单元，位移是连续的。

3. 随着单元的缩小，可减小这种应力和应变的突变

例题3 对于例1单元，试计算单元应力。 解：{}

{}

{}{}

{}x e

e y m m i i

j j xy S S S S σσσδδδδτ????

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

?????

?

????

??===++ {}{}

0i δ=

{}{}

{}e

m m j j

S S σδδ??????????=+∴ 22

220202(1)2(1)0111220002(1)102j j j j j j j m b

c E E S b c c

b E S μμμμμμμμμμμ?

?

????

?

???????

??

?????????

???

??

??????

?

?

??

???????????

??????

?

==-?----=--

所以单元应力为

{}

2280200.0020200000.0012(1)101020.0040.0010.0040.0012(1)03.3670.210()

0e E E MPa μσμμμμμμμ????

???

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

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

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

?

?

??????

????

??????

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

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