弹性力学的有限元求法

1 1
xi
A 1 2
xj
1 xm
yi
1a 1
yj
1 2
0
ym 1 0
0 0 1 a2
2 a
例2.2、3结点三角形单元如图所示。结点位移如下， ui=0.1a, vi=0.05a, uj=0.15a, vj=0.1a, um=0.05a, vm=0。单元内 一点p的坐标为(0.25a,0.5a)，求p点的位移。
虚功相等，
{ *}e T{R}e
*
e
T
[N]T {q}tds
s
{R}e [N]T{q}tds s
例题2.5、在均质、等厚的三角形单元ijm的ij边上作用有 沿x方向按三角形分布的载荷，求移置后的结点载荷。
Rxi Ni
Ryi
0
Rxj Ryj
s
Nj 0
Rxm
m
Rym 0
0
Ni
0 Px
N
j
Py
0
Nm
Ni (0.25a,0.5a) 0.25 N j (0.25a,0.5a) 0.5 Nm (0.25a,0.5a) 0.25
2）体力的移置
令单元所受的均匀分布体力为，
{q}
qx
q
y
虚功相等，
{ *}e T{R}e { *}e T [N]T{q}tdxdy
u a1 a2x a3 y a4x2 a5xy a6 y2 ... v b1 b2x b3 y b4x2 b5xy b6 y2 ...
3结点三角形单元内位移由结点的6个位移分量来 确定。六个位移分量只能确定六个多项式的系数， 所以3结点三角形单元的位移函数如下，
u v
(3╳6)
Ke [ B]T[ D][ B]tA
(6╳6)
4.4.1 单元刚度矩阵
ci xm x j
i,j,m坐标轮换
aj bj
am bm
c j cm
(2-10)
aa12 a3
1 2A
abii ci
aj bj cj
am bm
uuij
cm um
(2-11)
同样，将垂直位移分量与结点坐标代入公式（2-6）中 的第二式，可得，
a4 a5 a6
1 2A
ai bi ci
vm
单元内的位移函数可以简写成，
f N e
(2-14)
把[N]称为形态矩阵，Ni称为形态函数。
例2.1、3结点三角形单元如图所示，求其形态矩 阵[N]。
结点坐标为，
i(a,0), j(0, a),m(0,0)
ai x j ym xm y j 0 0 0 a 0 bi y j ym a 0 a
f ( p) N( p) e
u( p) Ni ( p)ui N j ( p)u j Nm ( p)um
0.25a 0.1a 0.5a 0.15a (1 0.25a 0.5a )0.05a
a
a
aa
0.025 a 0.075 a 0.0125 a 0.1125 a
v( p) Ni ( p)vi N j ( p)v j Nm ( p)vm
f * [N] * e
令结点载荷为
Rxi
Ryi
Re
Rxj Ryj
Rxm
Rym
1）集中力的移置
在单元内任意一点p作用集中力
{ f p} [N]
{ *}e
( xp, yp)
P
Px Py
{ *}e T {R}e { f P}T {P}
{ *}e T {R}e { *}e T ([N ]
§6.1 有限元网格(Finite Element[Nk] etwork)
1. 常用单元 (1) 线性单元(linear elements)：三角形、矩形或其他四边形。形 函数是线性函数，即单元内任一点的坐标可用单元节点坐标的线性 函数来表示。
图5.3 线形单元及节点位置
(2) 抛物形单元(parabola elements)：除了角上有节点（主外节 点），边缘上也有节点（副外节点）的单元。有时在内部也有节点 （内节点），直边或曲边均可。单元内任一点的坐标，可用一个抛 物线内插法来求得。
uuij um
1 1 1
xi xj xm
yi yj ym
aaa132
(2-7)
1 xi 1 x j 1 xm
yi
y
j
T
ym
aa12
T1
ui uj
(2-8)
a3
um
[T]1 [T]* T
T 2A A为三角形单元的面积
[T]的伴随矩阵为，
T
T
*
x j ym xm yi
假定结点位移： ui 1, vi 0 u j 1, v j 0 um 1, vm 0
u(x, y) Ni (x, y)ui N j (x, y)u j Nm (x, y)um
按照上面给出的结点位移，单元实际上产生了刚体位移， 对于单元中的任意一个点u(x,y)=1。所以，
1 Ni (x, y) N j (x, y) Nm (x, y)
{R}e [N]T {q}tdxdy
(2-17)
例题2.4、设有均质等厚的三角形单元ijm，受 到沿y方向的均布载荷qy的作用。求均布体力 移置到各结点的载荷。
Rxi Ni
Ryi
0
Rxj Ryj
N 0
j
Rxm
N
m
Rym 0
0
Ni
0 0
N
j
q
y
t
dxLeabharlann Baidu
dy
0
Nm
0.25a 0.05a 0.5a 0.1a (1 0.25a 0.5a )0
a
a
aa
0.0125 a 0.05a 0.0625 a
形态函数Ni具有以下性质： 1）在单元结点上形态函数的值为1或为0。 2）在单元中的任意一点上，三个形态函数之和
等于1。
对于任意一点p(x,y)，
Ni (x, y) N j (x, y) Nm (x, y) 1
6.2单元位移函数
1
常数项1
xy
线性项3
x2 xy y2
二次项6
x3 x2y xy2 y3 立方项10
x4 x3y x2y2 xy3 y4 四次项 15
弹性体内的实际位移分布可以用单元内的位移分
布函数来分片近似地表示。在单元内的位移变化 可以假定一个函数来表示，这个函数称为单元位 移函数、或单元位移模式。单元位移函数可以用 多项式表示，
xm y j xi ym
y j ym ym yi
xm xi
xj xm
(2-9)
xi y j x j yi yi y j x j xi
ai x j ym xm y j bi y j ym
ai [T]* a j
bi bj
ci cj
T
abii
am bm cm ci
a1 a4
a2 a5
x x
a3 a6
y y
(2-6)
将3个结点上的坐标和位移分量代入公式（2-6） 就可以将六个待定系数用结点坐标和位移分量 表示出来。首先计算位移分量u的系数，
ui a1 a2 xi a3 yi u j a1 a2 x j a3 y j um a1 a2 xm a3 ym
aj bj cj
am vi
bm
v
j
cm vm
(2-12)
将（2-11）、（2-12）代回（2-6）整理后可得，
u
1 2A
[(ai
bi x
ci
y)ui
(a j
bj x
cj
y)u
j
(am
bm x
cm
y)um ]
v
1 2A
[(ai
bi x
ci
y)vi
(a j
bjx
cj
y)v j
(am
bm x
cm
y)vm ]
令
Ni
1 2A
(ai
bi x
ci
y)
（下标i，j，m轮换）
ui
vi
u
v
Ni
0
0 Ni
Nj 0
0 Nj
Nm 0
0 Nm
u v
j j
(2-13)
um
vm
单元内的位移记为 单元的结点位移记为
f
u v
ui
e
i j
m
vi
u
v
j j
um
)T {P}
( xp, yp)
由于虚位移是任意的，
{R}e ([N ]
)T {P}
( xp, yp)
(2-16)
例2.3、在均质，等厚的三角形单元ijm的一点p （0.25a, 0.5a）上作用有集中载荷Px, Py。
Rxi Ni
Ryi
0
Rxj Ryj
Nj 0
Rxm
N
位移函数构造和收敛性要求
▪ 单元中的位移模式采用待定系数的有限多项式做 为近似函数，有限多项式选取的原则：
－待定系数是由结点位移条件确定的，三角形6结 点位移，矩形单元待定系数8结点位移。
－在选择多项式时，必须要选择常数项和完备的一 次项。因为这两项可以反映单元刚体位移和常应 变的特性。
－选择多项式应由低阶到高阶，并具有坐标对称性。
N
m
Rym 0
0
Ni
0 N
j
q0x
tds
0
Nm
取局部坐标s，在i点s=0，在j点s=l，L为ij边的长度。
在ij边上，以局部坐标表示的插值函数为，
Ni
1
s L
s Nj L
载荷为
qx
q
s L
Nm 0
Rxi
L 0
(1
s )q L
s L
tds
qt( s2 2L
s3 3L2
)
L 1 qtL 06
Rxj
L 0
s L
q
s L
tds
qt
s3 3L2
L 1 qtL 03
设ij边的长度为L，先把分布面力 等效为作用在距 i结点2/3L处P点 的集中力，再移置到结点上。
4.4 单元刚度矩阵的求解过程
Fe
(6) BTtA
(6╳3)
s
(3)
e D (3)
(3╳3)
e
B
(3╳6)
S [ D][ B]
三角形积为 形态函数为
A a2 2
Ni
1 2A
(ai
bi x
ci y)
1 a2
(0 ax
0)
x a
N
j
1 2A
(a j
bj x
cj
y)
1 a2
(0 0
ay)
y a
Nm
1 2A
(am
bm x
cm y)
1 a2
(a2
ax
ay)
1
x a
y a
形态矩阵为
x
[N
]
a
0
0 x
y a 0
图5.4 抛物线形单元及节点位置
2. 单元划分要注意的几个问题 (1) 相邻两个单元的节点要与节点重合（外节点与外节点、内节点 与内节点），不能与无节点边重合。
( ,)
(2) 单元不必是相同尺寸，应力有突变的地方，单元划分应较小。 (3) 任何一个单元必须只能在一种材料区，即它不能跨越两种材料
区。 (4) 同一单元的各个边长，一般不要相差太大。
载荷移置要满足静力等效原则。静力等效是指 原载荷与结点载荷在任意虚位移上做的虚功 (virtual work)相等。
处于平衡状态的弹性体，真实位移分量发生了 位移边界条件所允许的微小改变，这个微小的 改变就是虚位移(virtual displacement)。
结点的虚位移
* e
单元的虚位移
三角形中的一点P可以用子三角形面 积定义的自然坐标来确定。面积坐 标定义为，
Li
Ai A
Lj
Aj A
Lm
Am A
点P表示为，P(Li , Lj , Lm )
面积坐标在三角形全面积上的积分为
A
Lai Lbj Lcmdxdy
(a
a!b!c! bc
2)!
A
3）分布面力的移置
设在单元的边上分布有面力， q [qx , qy ]T
ci xm x j 0 0
a j xm yi xi ym 0 0 a 0 0 bj ym yi 0 0 0
c j xi xm a 0 a
am xi y j x j yi a a 0 0 a2
bm yi y j 0 a a
cm x j xi 0 a a
0 y
1 x y aa 0
0
1
x
y
a
a
a a
用来计算三角形面积时，要注意单元结点的排列顺序。
当三个结点i，j，m取逆时针顺序时，用行列式计算出
的三角形面积为正值。
1 A 1 1
2
xi xj
1 xm
yi yj
1 11
2
a 0
0 a 1 a2
2
ym 1 0 0
如果把三个结点按顺时针方向排列，即i(a,0)、j(0,0)、 m(0,a)，行列式的计算结果为负值。
形态函数的几何意义
1 Ni 2A (ai bi x ci y)
1x y
Ni
11 2A
xj
yj
1 xm ym
Ni
SPJM SMIJ
Nj
SPMI SIJM
Nm
SPIJ SIJM
任意一点P的形态函数Ni是点P与结点I的对边所构 成的三角形面积与整个单元面积之比。
2.3单元载荷移置
有限元法的求解对象是单元的组合体，待求解 的未知变量都定义在单元的结点上，因此作用 在弹性体上的外力，需要移置到相应的单元结 点上成为结点载荷。
Rxi 0, Rxj 0, Rxm 0
Ryi Niqytdxdy qyt Nidxdy
N i dx dy
1 2A
(ai
bi
x
ci
y)dxdy
1 2A [ai
A
bi
Axc
ci
Ayc ]
A
1 2A
(ai
bi xc
ci
yc )
1 3
A
Ryi
1 3
qy At
Ryj
1 3
q
y
At
Rym
1 3 qy At