24四边形单元

再将应变矩阵(用形函数与节点位移表示)代入应力矩阵中，有：
西
工 学
σ = D ⋅ B ⋅qe = S ⋅qe
院
(3×1) (3×3) (3×6) (6×1) (3×6) (6×1)
汽
其中，应力矩阵为：
车
工
程
S = D⋅B
系
(3×6) (3×3) (3×6)
单元势能的表达
以上将三大基本变量(u,ε ,σ ) 均用基于单元节点位移列阵来进
y)
⎥⎦
=
(3×2) (2×1)
∂⎥
∂x
⎥ ⎦
工 程
其中 [∂] 为几何方程的算子矩阵，即
系
⎡∂
⎢ ⎢
∂x
⎤ 0⎥
⎥
[ ]∂ =
⎢ ⎢
0
⎢ (3×2)
∂⎥
∂y
⎥ ⎥
⎢∂ ∂⎥
⎢ ⎣
∂y
∂x
⎥ ⎦
将节点位移代入上式，则有：
[ ] ε (x, y) = ∂ N (x, y) ⋅qe = B (x, y) ⋅qe
车 工 程
1 ci = 1
xj xm
= −x j + xm
系
a,b,c编号以顺序方式记即可。
将 β1, β2 , β3, β4 , β5, β6代入u(x,y), v(x,y)表达式中，有：
u(x, y) = Ni (x, y)ui + N j (x, y)u j + Nm (x, y)um v(x, y) = Ni (x, y)vi + N j (x, y)v j + Nm (x, y)vm
工
将各方向位移的位移模式代入节点位移表达式中，有：
程 系
⎧ui (xi , yi ) = β1 + β2 xi + β3 yi
⎪⎪vi (xi , yi ) = β4 + β5 xi + β6 yi
⎪⎪u j (x j , y j ) = β1 + β2 x j + β3 y j ⎨⎪v j (x j , y j ) = β4 + β5 x j + β6 y j
(6×6) (6×3) (3×3) (3×6)
⎢⎣kmi kmj kmm ⎥⎦
其中，
kii = BTi D Bi (2×3) (3×3) (3×2)
kij = BTi D B j (2×3) (3×3) (3×2)
kim = BTi D Bm ...... (2×3) (3×3) (3×2)
单元刚度方程
= 1 qeT K eqe − PeT qe 2
系
其中，Ke是单元刚度矩阵，即：
K ∫ B D B dΩ ∫ e =
T
(6×6)
Ωe (6×3) (3×3) (3×6)
= BT DB ⋅ dA⋅t Ae
t为平面问题单元的厚度。
广 西
由式
B (x,
(3×6)
y)
=
1 2A
⎡⎢⎢b0i ⎢⎣ci
0 ci bi
工 学 院
vm 节 点 m(xm,ym) um
汽 车
vj
vi
uj
工
节点j(xj,yj)
程
ui
系
y 节点i(xi,yi)
x
O
对于图示的平面3节点三角形单元，由于有3个节点，每一个节点有两个 位移，因此，共有6 个节点位移，即6个自由度。如果我们将x方向和y方 向的位移进行分别表达，则每个方向的位移由3个节点位移来确定，即 每个方向的位移可以设定3个待定系数。
⎢∂ ∂⎥
程 系
⎢ ⎣
∂y
∂x
⎥ ⎦
将
Ni
=
1 2A
(ai
+
bi x
+
ci
y)
(i, j, m)
代入上式，有：
B (x,
(3×6)
y)
=
1 2A
⎢⎢⎡b0i ⎢⎣ci
0 ci bi
bj 0 cj
0 cj bj
bm 0 cm
B 0
cm bm
⎤ ⎥ ⎥ ⎥⎦
=
⎡
⎢ ⎢⎣
(
i
3×2)
Bj
(3×2)
⎤
Bm
(3×2)
⎥ ⎥⎦
广 来自百度文库 工 学
其中，
Bi
(3×2)
=
1 2A
⎡bi
⎢ ⎢
0
⎢⎣ci
0⎤
ci
⎥ ⎥
(i, j, m)
bi ⎥⎦
院
汽 车
当单元的节点坐标确定后，这些参数都是常量(与坐标变量x,y
工
无关)，因此，B是常量矩阵。当单元的节点位移确定后，由
程 系
B转换求得的单元应变都是常量，也就是说在载荷作用下，
单元中各点具有同样的应变值，因此，3节点三角形单元也称
为常应变单元。
单元应力场的表达
由弹性力学中平面问题的物理方程，将其写成矩阵形式，为：
广 西 工 学 院
⎡
⎤
σ (x,
(3×1)
y,
z)
=
⎡σ ⎢⎢σ
x y
⎢⎣τ xy
⎤ ⎥ ⎥ ⎥⎦
=
E
1− μ2
⎢1
⎢⎢μ
⎢ ⎢0
μ
1
0
ε D 0
0
⎥ ⎥ ⎥
⎡ ⎢ ⎢
εx εy
⎤ ⎥ ⎥
=
⋅
1
−
μ
⎥ ⎥
⎢⎣γ
xy
bj 0 cj
0 cj bj
bm 0 cm
B 0
cm bm
⎤ ⎥ ⎥ ⎥⎦
=
⎡ ⎢⎢⎣(3×2i)
Bj
(3×2)
⎤
Bm
(3×2)
⎥ ⎥⎦
工
学
院 可知，B矩阵为常系数矩阵，因此，上式可以写成：
汽
车
工 程 系
K B D B t A e = T
=
⎡ ⎢ ⎢
kii k ji
kij k jj
kim k jm
⎤ ⎥ ⎥
x O
u2 节点2(x2,y2)
汽 车 工
u=0, v=0
三节点三角形单元如所示，三个节点的编
号为1,2,3，各自位置坐标为(xi , yi ), i = 1, 2,3
程
各个节点的位移(分别沿x方向和y方向)为
系
(ui , vi ), i = 1, 2,3
三角形单元共有6个节点位移自由度。将所有节点上的位移组成一个位 移列阵，记为qe，同样，将所有节点上的各个力也组成一个力列阵， 记为Pe，则有：
广
行表达了，将其代入单元的势能表达式中，有：
西
工 学 院
∫ ∫ ∫ ∏e = 1 σ T ⋅ε dΩ −[ b T ⋅udΩ + pT ⋅udA]
2 Ωe
Ωe
S
e p
汽
∫ ∫ ∫ = 1 qeT ( BT DBdΩ)qe − ( N TbdΩ + N T pT dA)qe
2
Ωe
Ωe
S
e p
车 工 程
⎥⎦
(3×3) (3×3)
⎣
2⎦
汽
其中平面应力问题的弹性系数矩阵为：
车
工
⎡
⎤
程
⎢1 μ 0 ⎥
系
D(3×3)
=
E
1− μ2
⎢⎢μ
⎢
⎢0
1 0
⎥ 0⎥
1− μ ⎥
⎥
⎣
2⎦
若为平面应变问题，则将上式中的系数换成平面应变问题的系数即可， 即将 (E, μ) 换成 ( E , μ )
1− μ2 1− μ
广
考虑到简单性、完备性、连续性及待定系数的唯一确定性分别选取单 元中各个方向的位移模式为：
广 西 工
⎧u ( x, ⎨⎩v( x,
y) y)
= =
β1 β4
+ +
β2 β5
x x
+ +
β3 β6
y y
学 院
同样，根据节点条件，即：
汽 车
⎧⎨⎩uv((xxkk
, ,
yk yk
) )
= =
uk vk
(k = i, j, m)
e p
工 程
= 1 qeT K eqe − PeT qe 2
系
⇒ K e qe = Pe (6×6) (6×1) (6×1)
广义坐标有限元位移模式的一般原则
广
A. 广义坐标是由结点场变量确定的，因此它的个数应与结点自由度
西 工
数相等。如平面三角形单元有6个结点自由度(结点位移)，广义坐
学
标个数应取6个，因此两个方向的位移u和v各取三项多项式。
qe = ⎡⎣ui vi u j v j um vm ⎤⎦T
6×1
Pe = ⎡⎣Pxi Pyi Pxj Pyj Pxm Pym ⎤⎦T 6×1
若该单元承受外载，则可以将外载等效到单元节点上，即等效节点
力。利用函数插值、几何方程、物理方程以及势能计算公式，可以
广 西
将单元的所有力学参量用节点位移列阵qe及相关的插值函数来表示。
项式的次数，以保证相邻单元的交界面上位移协调性。
不同形式广义坐标有限元的位移模式
广 西
单元形式
工 学
3结点三角形平面单元
院
6结点三角形平面单元
汽 车
4结点四边形平面单元
工 程
8结点四边形平面单元
系
4结点四面体三维单元
8结点六面体三维单元
位移模式 1xy
1 x y x2 xy y2 1 x y xy
1 x y x2 xy y2 x2y xy2 1xyz
1 x y z xy yz zx xyz
平面四边形单元
平面4节点矩形单元如图所示，单元的节点位移共有8个自由度，平
广
面的编号为1,2,3,4，各自的位置坐标为(xi,yi) i=1,2,3,4，各个节点的位
西 工
移(分别沿x方向与y方向)为(ui, yi) i=1,2,3,4。
yi yj
=
1 2A
(biui
+
bju
j
+
bmum
)
ym
A为三角形单元的面积
汽 车 工 程
1
β3
=
1 2A
1
xi xj
ui uj
=
1 2A
(ciui
+
c
ju
j
+
cmum
)
1 xm um
系
β4
=
1 2A
(aivi
+
ajvj
+
amvm )
β5
=
1 2A
(bivi
+ bjvj
+ bmvm )
β6
=
1 2A
汽
⎢⎣vm ⎥⎦
车
其中，N(x,y)为形状函数矩阵，即：
工
程 系
N (x,
( 2×6 )
y)
=
⎡Ni
⎢ ⎣
0
0 Ni
Nj 0
0 Nj
Nm 0
0⎤
N
m
⎥ ⎦
Ni
=
1 2A
(ai
+
bi x
+
ci
y)
(i, j, m)
单元应变场的表达
由弹性力学平面问题的几何方程(矩阵形式)
广 西
⎡ ∂u ⎤ ⎡ ∂
⎢
⎥⎢
⎪⎪um (xm , ym ) = β1 + β2 xm + β3 ym
⎪⎩vm (xm , ym ) = β4 + β5 xm + β6 ym
广
β1
=
1 2A
ui uj
xi xj
yi yj
=
1 2A
(aiui
+
a ju j
+
amum )
um xm ym
西 工 学 院
1
β2
=
1 2A
1
ui uj
1 um
学 院
v2 Node 2 u2
v1 u1
汽 车 工
(x2, y2) b
y,η x,ξ
O
Node 1 (x1, y1)
程 系
b
Nodev33 (x3, y3)
u3
v4 u4
Node 4
a
a
(x4, y4)
若采用无量纲坐标
v2 Node 2 u2 (x2, y2)
b
y,η x,ξ
O
v1 u1
Node 1 (x1, y1)
广
将单元的势能对节点位移取一阶极值，可得到单元的刚度方程：
西
工
学 院
∫ ∫ ∫ ∏e = 1 σ T ⋅ε dΩ −[ b T ⋅udΩ + pT ⋅udA]
2 Ωe
Ωe
S
e p
汽 车
∫ ∫ ∫ = 1 qeT ( BT DBdΩ)qe − ( N TbdΩ + N T pT dA)qe
2
Ωe
Ωe
S
广 西 工 学 院 汽 车 工 程 系
四边形单元与三角形单元的异同
广 西 工 学 院
汽
Node 3
车
工
程
系
Node 6
Node 9
Node 10
Node 4
Node 8
Node 12
Node 13
T1
①
T2
v3
节点3(x3,y3)
u3
v2
②③
广
④
西
⑤
工
⑦
学
⑧
⑥
院
⑩⑨
v1
u1 y 节点1(x1,y1)
(civi
+
cjvj
+
cmvm )
其中：
1 A= 1 1
2
xi xj
yi yj
=
1 2
(ai
+
aj
+
am )
=
1 2
(bic j
−
b j ci
)
广
1 xm ym
西
工 学
ai
=
xj xm
yj ym
= x j ym − xm y j
院 汽
1 bi = − 1
yj ym
= y j − ym
ai-cm
等9个参数的值均由单 元的节点坐标确定
⎤ 0⎥
工 学 院 汽 车
ε(
(3×1)
x,
y)
=
⎡ ⎢ ⎢
ε ε
x y
⎢⎣γ xy
⎤ ⎥ ⎥ ⎥⎦
=
⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
∂x ∂v ∂y ∂u + ∂y
∂v ∂x
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
=
⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
∂x 0
∂ ∂y
⎥
[∂] u ∂ ⎥ ⎡u(x, y)⎤
∂y
⎥ ⎥
⎢⎣v(
x,
广
(3×1)
(3×2) (2×6)
(6×1) (3×6)
(6×1)
西 工
其中几何函数矩阵B(x,y)为：
学 院
⎡∂
⎢ ⎢
∂x
⎤ 0⎥
⎥
汽 车 工
[ ] B (x, y) =
∂
N
=
⎢ ⎢
0
(3×6)
⎢ (3×2) (2×6)
∂ ⎥ ⎡Ni 0
∂y
⎥ ⎥
⎢ ⎣
0
Ni
Nj 0
0 Nj
Nm 0 ⎤
0
N
m
⎥ ⎦
同样，将所有节点上的各个节点力组成一个列阵，记作Pe，则有：
写成矩阵形式，有：
广 西
⎡ ui ⎤
⎢ ⎢
vi
⎥ ⎥
工 学 院
u N q (x,
( 2×1)
y)
=
⎡u(x, ⎢⎣ v( x,
y)⎤ y ) ⎥⎦
=
⎡ ⎢ ⎣
Ni 0
0 Ni
Nj 0
0 Nj
Nm 0
0 Nm
⎤ ⎥ ⎦
⎢ ⎢ ⎢
u v
j j
⎥ ⎥ ⎥
=
(x,
( 2×6 )
y)
⋅
e
(6×1)
⎢⎢um
⎥ ⎥
广 西
ξ = x,η = y
ab
b
Nodev33 (x3, y3)
u3
工 学
则单元4个节点的几何位置为：
a
院
⎧ξ1 = 1,η1 = 1
汽 车 工
⎪⎪⎨⎪ξξ32
= =
−1,η2 −1,η3
= =
1 −1
程 系
⎪⎩ξ4 = 1,η4 = −1
v4 u4
Node 4
a
(x4, y4)
将所有节点上的位移组成一个列阵，记作qe。
院
B. 选取多项式时，常数项和坐标的一次项必须完备。位移模式中的
汽 车
常数项和一次项反映了单元的刚体位移和常应变的特性。
工 程
C. 多项式的选取应由低阶到高阶，尽量选取完全多项式以提高单元
系
的精度。若由于项数限制不能选取完全多项式时，选择的多项式
应具有坐标的对称性。并且一个坐标方向的次数不应超过完全多