有限元分析及应用第六章

e
−1 −1
(6-1-7)
单元变形能
第 3 页 共 10 页
有限元分析与应用
霍战鹏
[ ] Ve
=
1 uΙ T
2
u
Π
k
uΙ u Π
=
1 2
uΙ u Π
T
k ΙΙ k ΠΙ
k ΙΠ k ΠΠ
uΙ u Π
由于[k]为对称阵，必有
[ ] kΠΙ T = [kΙΠ ]
以及
来自百度文库
(6-1-8)
{ } { } { } { } { } { } { } { } [ ] [ ] [ ] [ ] Ve
=
1 2
(
uΙ
T
k ΙΙ
uΙ
+ uΠ T kΠΙ
uΙ
+ uΙ T kΙΠ
uΠ
+ uΠ T kΠΠ
uΠ )
{ } { } { } { } { } { } [ ] [ ] [ ] =
1 2
uΙ
T
k ΙΙ
uΙ
+ uΠ T kΠΙ
uΙ
+1 2
uΠ
T
k ΠΠ
uΠ
(6-1-9)
{ } 体积力 f x , f y T 做功
系统的总势能定义为
m
m
∑ ∑ ∑ π
h Ｐ
=
Vej −
Wej −
WS
j =1
j =1
(6-1-2)
∑ 其中 Vej , j 号单元的变形能；Wej , j 号单元体积力做功； WS
4
4
∑ ∑ N
i
u
、
i
N i vi 上做的功之和。
i =1
i =1
我们将把由
∂π
h P
=0
(i = 1 ~ n）
∂ui
为各边界外力在位移
][kΠΠ
] {−1 rΠ
}
~
第 4 页 共 10 页
(6-1-15)
有限元分析与应用
霍战鹏
为了避免对一个 4 × 4的方阵[kⅡ]求逆，静凝聚可以分四次进行，每次只凝聚一个内自由度。
{ } （４） 边界外载荷 p x p y T 的等效节点力
∫ {rs
}=
sP
N1 0
0 N1
N2 0
0 N2
N3 0
v2
η
y,v
2 u2
４
３
4
e
v1 Ｍ
1
u1
0
ξ ê
１
Mˆ
２
x,u
（ξ,-1）
图６－3
（３） 协调性分析
沿单元的一边，例如节点１、２所在的边，η＝－１。u,v 是ξ的二次函数，完全被 u1, v1,α1, 和 u2, v2,α3 所决定。但由于不同单元的α1~α4 彼此独立，故不能保证单元之间位移 的协调性。对于不满足协调条件的单元，显然不能再用最小势能原理。尽管如此，人们仍然
有限元分析与应用
霍战鹏
第六章 非协调单元
第二章至第五章的讨论以最小势能原理为基础，要求在单元内假设的位移场（试探函数） 满足协调条件（在不同的单元内可以假设不同的的位移场）。满足协调条件的单元，它们的 收敛性等问题已在第四章中做了研究。等参数单元就是目前处理二阶问题时应用最广的一种 协调单元。
此外，还有一些单元，它们不满足协调条件，但仍可以收敛到真实解，这类单元称为非 协调单元，可以看成是对等参数单元的一种改进，目的在于：在计算量增加不多的情况下， 使单元的实际精度有所改善。对于四阶问题（例如板、壳），协调条件要求单元之间位移和 位移的一阶导数（转角）连续。在第七章中将会看到，实现上述协调条件不是件容易的事， 而且为此要增加相当大的计算量，因而人们在自编程序中常常对非协调单元感兴趣。
至此，我们对非协调元的印象是：（i）分单元假设的位移场（即试探函数）不完全满足 协调条件；（ii）形式上套用了协调元的具体作法。至于能否收敛到真实解，到目前为止并不 清楚。实际情况是：有时能保证收敛性，有时则不能。为了明确回答这些问题，必须涉及有 关非协调元的数学理论。
§6-2 非协调元的基本理论
为了简单，以泊松方程为例，基本未知量只有一个，变形能的表达式也比较简单。 １、有限元空间 Wilson 非协调元的基函数中与四节点等参元的基函数相同者记作 φi (i=1~n)。与单元内 自由度有关的形函数记作 ψl (l=1~2m)。其中 φi 满足协调条件。ψl 不满足协调条件，穿过 单元边界时 ψl 有有限跳跃量。即
{ } { } 单元的外自由度为： uΙ = u1 v1 u2 v2 u3 v3 u4 v4 T { } { } 单元的内自由度为： uΙΙ = α1 α 2 α 3 α 4 T
（6-1-1）所定义的单元位移场可写成
u v
=
N01
0 N1
N2 0
0 N2
N3 0
0 N3
N4 0
0 N4
(1−ξ2) 0
照搬协调单元的一些具体做法（这样，对协调元编制的程序基本上可以照搬）。至于能否保
证收敛到真实解等问题，留待下面几节加以讨论。
现假如我们研究某个具体的平面应力问题。节点总数为 n，单元总数为 m。则总的自由
度可区分为：
节点自由度 ui , vi
（ i = 1 ~ n）
单元内自由度 α１( j)、α２( j)、α３( j)、α４( j) （ i = 1 ~ m）
α2(j)、α3(j)、α4(j)仅属于第 j 号单元，故有
∂π
h P
∂α
( l
j
)
=
∂Vej
∂α
( l
j
)
−
∂Wej
∂α
( l
j
)
=0
(l = 1 ~ 4）
(6-1-4)
这样，就可以在单元分析时先消去αl(j)（这一步骤称为静凝聚），只剩下 ui, vi 进入总体平衡 方程。
４、单元分析 静凝聚
在单元分析中，节点仍取它在单元内的序号，并约定：
[K ]{U }= {F }
可得到节点位移 ui、vi (i=1－n) 具体作法与协调单元作法相同。 最后强调一点：静凝聚与非协调元是两个不同的概念。静凝聚的目的是消去内自由度，
以减少总体平衡方程的规模。不论是协调元还是非协调元，只要存在内自由度，都可以在 单元分析过程中将内自由度凝聚掉；不进行静凝聚，除了增加解总体方程的工作量外不会 有其它任何（好的或不好的）影响。静凝聚方法也广泛地用于组合单元和子结构，是一项 很有价值的技术。
(1−η2) 0
0 (1−ξ2)
0 uΙ
(1−η2)uΙΙ
（１） 几何矩阵
应变
ε ε
x y
τ xy
=
∂ ∂x
0
∂
∂y
0
∂
∂y ∂
uv
=
[B
] uΙ
u ΙΙ
∂x
(6-1-5)
其中[B]为几何矩阵。
∂∂x
[B] = 0
0 ∂∂yN01
0 N1
N2 0
0 N2
N3 0
0 T
0
0 (1−ξ 2 ) (1−η2 )
×
f f
x y
t
det
J
dξdη
(6-1-11)
（３） 静凝聚
略去（6-1-4）中的单元编号 j，以（6-1-9），（6-1-10）代入，（变形能对内部自由度取偏
导）并注意到
则有 解得
{ } { } uΙΙ = α1 α 2 α 3 α 4 T { } { } { } [kΠΙ ] uΙ + [kΠΠ ] uΠ − rΠ = 0
积）。基函数 ψl，j-1、ψl,j 仅在第 j 号单元内的非零，且
ψ l, j−1 = (1 − ξ 2 ) ψ l, j = (1 − η 2 )
(6-2-2)
对一般的非协调元来说，它的基函数或者一部分不满足协调条件。或者全部不满足协调 条件。这些基函数张成的有限元空间 Sh 均不是容许空间（二阶问题的容许空间为 H1）的子
本章只讨论二阶问题，主要包括：非协调元的构造和分析方法，非协调元的理论基础（显 然不能再利用最小势能原理），收敛判别方法。这些结论对四阶问题同样适用。
从关于非协调元的讨论中，读者可以看到，有限元方法有了坚实的数学基础以后，在构 造方法时思路可以开阔很多。
§6-1Wilson 非协调元
Wilson 非协调元可以看成是由等参数单元演变来的单元，现以二维情况为例。
∂π
h P
=0
(i = 1 ~ n）
∂vi
∂π
h P
∂α
( l
j
)
=0
( j = 1 ~ m）
、(l = 1 ~ 4）
(6-1-3)
第 2 页 共 10 页
有限元分析与应用
霍战鹏
求得的 ui, vi, αl(j)以及由此求得的应力做为非协调单元的的有限元解。 在（6-1-3）中共有 2n+4m 个未知量。比四节点等参元多了４m 个未知量。但是α1(j)、
0 N3
N4 0
0 N4
T
p p
x y
tds
(6-1-16)
{ } 积分沿单元受到载荷的边界 SP 进行。其中 p x p y T 为作用在单位面积上的边界载荷。
５、组装及求解总体方程
引入由强[制k~]位组移装约总束体条刚件度。矩求阵解[~k]总。体由平{~r衡}、方{r程s}组装总体载荷向量{F}。在上述组装过程中可以
{ } { } { } [ ] [ ] [ ] uΠ = − kΠΠ −1 kΠΙ uΙ + kΠΠ −1 rΠ
(6-1-12)
将（6-1-12）代入（6-1-9）和（6-1-10）有
{ } { } { } Ve
− We
=
1 2
uΙ
T
k
uΙ
~
−
uΙ
T
r
~
−
1 2
{rΠ
}T
[k ΠΠ
] {−1 rΠ
图６－２
(6-1-1)
第 1 页 共 10 页
有限元分析与应用
霍战鹏
同四节点等参元相比，单元内假定的位移场多了四项：
α1 (1 − ξ 2 )、α 2 (1 − η 2 )、α 3 (1 − ξ 2 )、α 4 (1 − η 2 )
它们有如下特性： （１） 不影响节点处的位移值，故称αl 为非节点自由度或单元的“内自由度”。在计
}
(6-1-13)
（6-1-13）右端第三项与{uΙ}无关，不影响πPh 取驻值。第一项为{uΙ}的二次型，[k]为凝聚
掉内自由度后的单元刚度矩阵。{r}为凝聚内自由度后的载荷向量。它们的计算式为~
~
k
=
[kΙΙ
]
−
[kΙΠ
][kΠΠ
] [−1 kΠΙ
]
(6-1-14)
~
r
=
{rΙ
}−
[kΙΠ
0 N3
N4 0
0 N4
(1−ξ2) 0
(1−η2) 0
0 (1−ξ2)
0 (1−η2)
∂ ∂
∂y ∂x
(6-1-6)
（２） 单元刚度矩阵和体积力载荷向量
单元刚度矩阵，该矩阵为一12 ×12 阶的方阵
11
[k ] = ∫∫ [B]T [E][B]⋅ tdxdy = ∫ ∫ [B]T [E][B]det J tdξdη
∫∫ ∫ ∫ We
=
e
u v
T
f f
x y
tdσ
=
1 1 u T
−1
−1
v
f f
x y
t
det
J
dξdη
=
uΙ u Π
T
rrΠΙ
其中体积力载荷向量
(6-1-10)
∫ ∫ rrΠΙ
=
−11−11N01
0 N1
N2 0
0 N2
N3 0
0 N3
N4 0
0 N4
(1−ξ 2 ) (1−η2 ) 0
∑ u = N i (ξ ,η)ui + α1 (1 − ξ 2 ) + α 2 (1 − η 2 )
i =1
0
4
∑ v = N i (ξ ,η)vi + α 3 (1 − ξ 2 ) + α 4 (1 − η 2 )
i =1
4 (x4, y4)
3 (x3, y3)
2
e
(x2, y2)
1 (x1, y1) x,u
２ （1,-1）
F： eˆ → e
利用上面定义的形函数，坐标的变换可写成
图６－１
4
∑ x = N i xi i =1
4
∑ y = N i yi
y,v
i =1
其中（xi, yi）为实际单元中节点的坐标。 至此，还看不出 Wilson 非协调单元与
上一章介绍的等参数单元之间的差别。
３、单元内假设位移场
4
第 5 页 共 10 页
有限元分析与应用
霍战鹏
空间。 ２、内积和模 泊松方程
在整个求解域内其变形能为
∂2u ∂x 2
+
∂2u ∂y 2
=
−
f
(6-2-3)
∫∫Ω
1 2
(
∂u ∂x
)
2
+
(
∂u ∂y
)
2
dxdy
(6-2-4)
由于当 Sh 不是 H1(Ω)的子空间时，试探函数 u 在穿过单元边界时为 δ──函数。积分（6-2-4） 在整个求解域内积分不存在。但是，在单元 e 内积分
∂ψ l 、 ∂ψ l ∂x ∂y
为 δ 一函数。试探函数
n
m
∑ ∑ u =
uiϕi +
(α1(
ψ j ) l , j−1
+
α
( 2
j
)ψ
l
,
j
)
i =1
j =1
(6-2-1)
在穿过单元边界时也可能有有限跳跃量。故 φi、ψl 所张成的有限元空间Ｓh 仅是 L2 的一个子 空间（函数自身平方可积）。不是 H1（协调单元的有限元空间）的子空间（一阶导数平方可
算单元变形能和单元体积力做功时计入这些位移；但在计算边界外力做功（为了将边界力 化为等效节点力）时不计这些位移。即在计算边界外力做功时只计 Niui、Nivi 各项。
（２） 补充这些项后，单元内的位移场是ξ、η的完全二次多项式。当实际单元 e 为矩形时，单元内位移场将是 x、y 的完全二次多项式。
3
１、母体单元 形函数 母体单元 ê：边长为２的正方形 自然坐标：ξ、η 取四个角点为节点，在单元内的序号为１~4。
η
４ （-1,1）
３ （1, 1）
形函数
ξ
Ni
(ξ ,η)
=
1 4
(1 +
ξ iξ
)(1 + ηiη)
（i = 1 ~ 4）
ê
２、实际单元 e 可看成母体单元 ê 经变换 F 得到
１ （-1,-1）