有限元 2-5弹性力学平面问题(非协调单元 ) [兼容模式]

本章只讨论二阶问题，主要包括：非协调元的构造和分析方法， 非协调元的理论基础（显然不能再利用最小势能原理），收敛判别 方法。这些结论对四阶问题同样适用。
有限单元法
土木工程学院
§1 Wilson 非协调元
η ４（-1,1）
３（1, 1）
1. 母体单元 形函数
ê
ξ
母体单元ê：边长为２的正方形, 自然坐
（3）Βιβλιοθήκη Baidu调性分析
有限单元法
土木工程学院
沿单元的一边，例如节点１、
２所在的边，η ＝－１。u,v是 y,v
ξ的二次函数，完全被u1, v1,α1,
4
和u2, v2,α3 所决定。但由于不
同单元的α1~α4 彼此独立，故不 能保证单元之间位移的协调性。
0
3 能否保证收敛到真实解 ？
v2
u2
４
η ３
e
2
Wilson 非协调元的基函数中与四节点等参元的基函数相同者记作φi ( i=1~n)。与单元内自由度有关的形函数记作ψl (l =1~2m)。其中φi 满足协调条
件。Ψl 不满足协调条件，穿过单元边界时ψl 有有限跳跃量。即
∂ψ l 、 ∂ψ l
∂x ∂y
为δ函数。
试探函数
n
m
∑ ∑ u =
uiϕi +
且
ψ l, j−1 = (1 − ξ 2 ), ψ l, j = (1 −η 2 ) (2-2)
对一般的非协调元来说，它的基函数或者一部分不满足协调条件。或者全部 不满足协调条件。这些基函数张成的有限元空间Sh 均不是容许空间（二阶问题 的容许空间为H1）的子空间。
2. 内积和模
泊松方程
∂2u + ∂2u = − f ∂x 2 ∂y 2
3. 单元内假设位移场
i=1
e
(x2, y2)
4
4
∑ u = Ni (ξ ,η)ui + α1 (1 − ξ 2 ) + α 2 (1 −η 2 ) i =1
(x4, y4)
1 (x1, y1) x,u
0
4
∑ v = Ni (ξ ,η)vi + α 3 (1 − ξ 2 ) + α 4 (1 −η 2 )
− We
=
1 2
uE
T
⎡ ⎢k
⎤ ⎥
uE
⎣~⎦
− uE
T
⎧⎫ ⎨r ⎬ ⎩~⎭
−
1 2
{rI
}T
[k
II
] {−1 rI
}
(1-13)
式（1-13）右端第三项与{uE}无关，不影响πPh 取驻值。第一项为{uE}的二
次型，⎢⎡⎣k~
⎤ ⎥⎦
为凝聚掉内自由度后的单元刚度矩阵。
⎡⎢⎣r~
⎤ ⎥⎦
为凝聚内自由度后的载荷向
有限单元法
补充： 非协调单元
土木工程学院
有一些单元，它们不满足协调条件，但仍可以收敛到真实解， 这类单元称为非协调单元，可以看成是对等参数单元的一种改进， 目的在于：在计算量增加不多的情况下，使单元的实际精度有 所改善。
对于四阶问题（例如板、壳），协调条件要求单元之间位移 和位移的一阶导数（转角）连续。实现上述协调条件不是件容易 的事，而且为此要增加相当大的计算量，因而人们在自编程序中 常常对非协调单元感兴趣。
e
−1 −1
单元变形能
[ ] Ve
=
1 ⎧uE ⎫T
2
⎨ ⎩
u
I
⎬ ⎭
k
⎧u E
⎨ ⎩
u
I
⎫ ⎬ ⎭
=
1 2
⎧u E
⎨ ⎩
u
I
⎫T ⎬ ⎭
⎡k EE
⎢ ⎣
k
IE
k EI k II
⎤ ⎥ ⎦
⎧u E ⎨⎩u I
⎫ ⎬ ⎭
由于 [k] 为对称阵，必有
有限单元法
[ ] [ ] kIE T = kEI
(1-8)
Ｍ v1
ξ
ê
1
u1
１
Mˆ
２
x,u
（ξ,-1）
图3
平面应力问题为例。设节点总数为n，单元总数为m。则总的自由度可区分为：
节点自由度 ui , vi （i = 1 ~ n）
单元内自由度
α(j １
)、α２( j
)、α３( j
)、α４( j
)
（i = 1 ~ m）
系统的总势能定义为
m
m
Vej , j 号单元的变形能
⎨ ⎩
u
I
⎫T ⎬ ⎭
⎩⎨⎧rrEI
⎫ ⎬ ⎭
(1-10)
∫ ∫ ⎨⎧⎩rrEI
⎫ ⎬ ⎭
=
1
1
⎡ ⎢
N1
⎣ −1−1 0
0 N1
N2 0
0 N2
N3 0
0 N3
N4 0
0 N4
(1 − ξ 2 ) 0
(1 −η 2 ) 0
0 (1 − ξ 2 )
0 ⎤T (1 −η 2 )⎥⎦
×
⎧ ⎨ ⎩
f f
x y
(α
( 1
j
)ψ
l , j−1
+
α
( 2
j
)ψ
l, j
)
i =1
j =1
(2-1)
有限单元法
土木工程学院
在穿过单元边界时也可能有有限跳跃量。故φi、ψl 所张成的有限元空间Ｓh 仅是 L2 的一个子空间（函数自身平方可积）。不是H1（协调单元的有限元空间）
的子空间（一阶导数平方可积）。基函数ψl,j-1、ψl,j 仅在第 j 号单元内的非零，
⎪∂x
[B]
=
⎪ ⎨
0
⎪
∂ ∂y
⎪⎪⎡ ⎬⎢ ⎪⎣
N1 0
0 N1
N2 0
0 N2
N3 0
0 N3
N4 0
0 N4
(1− ξ 2 ) 0
(1−η 2 ) 0
0 (1− ξ 2 )
0⎤ (1−η 2 )⎥⎦
⎪∂ ∂⎪
⎪⎩∂y ∂x ⎪⎭
(1-6)
（2） 单元刚度矩阵和体积力载荷向量
11
[k] = ∫∫ [B]T [E][B]⋅ tdxdy = ∫ ∫ [B]T [E][B]det J tdξdη (1-7)
量。
有限单元法
土木工程学院
⎡ ⎢k
⎤ ⎥
=
[k I
]−
[k EI
][k II
] [k −1 IE
]
(1-14)
⎣~⎦
（４）
⎡ ⎢r
⎤ ⎥
=
{rE
}−
[k EI
][k II
] {−1 rI
}
⎣~⎦
(1-15)
{ } 边界外载荷 px p y T 的等效节点力
∫ {rs
}
=
sP
⎡N1
⎢ ⎣
0
0 N1
N2 0
土木工程学院
{ } { } { } { } { } { } { } { } [ ] [ ] [ ] [ ] Ve
=
1 2 ( uE
T
kI
uE
+ u I T k IE
uE
+ uE T kEI
uI
+ u I T k II
uI )
{ } { } { } { } { } { } [ ] [ ] [ ] =
1 2
uE
T
kI
uE
+ u I T k IE
uE
+1 2
uI
T
k II
uI
(1-9)
{ } 体积力 f x , f y T 做功
∫∫ ∫ ∫ We
=
e
⎧u ⎨ ⎩v
⎫T ⎬ ⎭
⎧ ⎨ ⎩
f f
x y
⎬⎫tdσ
⎭
=
1 1 ⎧u⎫T ⎨⎬
−1 −1 ⎩v ⎭
⎧ ⎨ ⎩
fx fy
⎫ ⎬t ⎭
det
J
dξdη
=
⎧u E
∂π
h P
=
∂Vej
− ∂Wej
=0
(l = 1 ~ 4）
∂α
( l
j
)
∂α
( l
j
)
∂α
( l
j
)
(1-4)
在单元分析时可以先消去αl (j) （这一步骤称为静凝聚），只剩下ui, vi 进入总 体平衡方程。
４. 单元分析 静凝聚 单元的外自由度： 单元的内自由度：
有限单元法
土木工程学院
{ } { } uE = u1 v1 u2 v2 u3 v3 u4 v4 T
∑ ∑ ∑ πＰh = Vej − Wej − WS (1-2)
j =1
j =1
Wej , j 号单元体积力做功
4
4
∑ ∑ ∑ WS 为各边界外力在位移
N
i
u
、
i
N i vi
上做的功之和
i =1
i =1
不计算边界力在内自由度上的功！
有限元解：
有限单元法
土木工程学院
由 方程组：
∂π
h P
=0
,
∂π
有限单元法
§2 非协调元的基本理论
土木工程学院
（i）分单元假设的位移场（即试探函数）不完全满足协调条件；
（ii）形式上套用了协调元的具体作法。至于能否收敛到真实解，到目前 为止并不清楚。实际情况是：有时能保证收敛性，有时则不能。
1. 有限元空间
（在讨论非协调元的数学理论时，为了简单，以泊松方程为例,基 本未知量只有一个，变形能的表达式也比较简单 ）
(2-3)
在整个求解域内其变形能
∫∫Ω
1 2
⎡ ⎢⎣(
∂u ∂x
)
2
+
(
∂u ∂y
)
2
⎤ ⎥⎦dx
dy
(2-4)
有限单元法
土木工程学院
由于当Sh不是H1(Ω)的子空间，试探函数 u 在穿过单元边界时为δ─ 函数。
积分（2-4）在整个求解域内积分不存在。但是，在单元 e 内积分存在。
∫∫e
1 2
⎡ ⎢( ⎣
0 N2
N3 0
0 N3
N4 0
５. 组装及求解总体方程
0 N4
⎤ ⎥ ⎦
T
⎧ ⎨ ⎩
p p
x y
⎫ ⎬tds ⎭
(1-16)
[K ]{U} = {F} 具体作法与协调单元作法相同
静凝聚与非协调元是两个不同的概念。静凝聚的目的是消去内自由度，以减少总体 平衡方程的规模。不论是协调元还是非协调元，只要存在内自由度，都可以在单元分析 过程中将内自由度凝聚掉；不进行静凝聚，除了增加解总体方程的工作量外不会有其它 任何（好的或不好的）影响。静凝聚方法也广泛地用于组合单元和子结构，是一项很有 价值的技术。
标：ξ、η取四个角点为节点，在单元内的序 １（-1,-1）
２（1,-1）
号为１~4。
形函数 Ni (ξ ,η)
=
1 4
(1 +
ξiξ )(1 + ηiη) （i
=1
~
图１
4）
2. 实际单元 e
4
4
∑ ∑ F： eˆ → e x = Ni xi , y = Ni yi y,v
3 (x3, y3) 2
i=1
∂u ∂x
)
2
+
(
∂u ∂y
)
2
⎤ ⎥ ⎦
dx
dy
存在
内积：当 u、v ∈ Sh 时，u、v 的内积定义为
∑ ∫∫ Dh (u,v) =
m j =1
ej
(( ∂u )( ∂v ) + (∂u )( ∂v ))dxdy ∂x ∂y ∂y ∂y
“变形能”：
（2-5）
∑ ∫∫ 1
2
Dh (u, u)
=
1 2
图２
i =1
(1-1)
单元内的位移场精度有所改善，二次函数
有限单元法
土木工程学院
同四节点等参元相比，单元内假定的位移场多了四项：
α1 (1 − ξ 2 )、α 2 (1 −η 2 )、α 3 (1 − ξ 2 )、α 4 (1 −η 2 )
这四项有如下特性： （1）不影响节点处的位移值，故称 αl 为非节点自由度或单元
Dh (u,u) − ( f ,u)
(2-9)
π Sh 中使
h P
取驻值的元素即为非协调元的有限元解 uh 。或者换个提法：找一
个元素uh∈Sh，使得对任何δu ∈Sh 都有
h P
=0
(i = 1 ~ n）;
∂ui
∂vi
(1-3)
∂π
h P
=0
( j = 1 ~ m） (l = 1 ~ 4）
∂α
( l
j
)
求得的 ui, vi, αl (j) 以及由此求得的应力做为非协调单元的的有限元解。
在（1-3）中共有2n+4m个未知量。比四节点等参元多了４m个未知量。但 是α1 (j)、α2 (j)、α3 (j)、α4 (j) 仅属于第 j 号单元，故有
⎫ ⎬t ⎭
det
J
dξdη
(1-11)
有限单元法
土木工程学院
（3）静凝聚
{ } { } 内自由度： uI = α1 α 2 α 3 α 4 T
略去（1-4）中的单元编号 j ，以（1-9），（1-10）代入，（变形能对内部 自由度取偏导）
∂π
h P
=
∂Vej
− ∂Wej
=0
∂α
( l
j
)
∂α
( l
j
)
∂α
( l
j
)
(l = 1 ~ 4）
{ } { } { } [kIE ] uE + [kII ] uI − rI = 0
{ } { } { } [ ] [ ] [ ] uI
= − k II
k −1 IE
uE
+ k II −1 rI
(1-12)
将（1-12）代入（1-9）和（1-10）有：
{ } { } { } Ve
0 ⎤⎧uE ⎫
(1
−
η
2
)⎥⎦
⎨ ⎩
u
I
⎬ ⎭
（1）几何矩阵 应变
⎡∂
⎢
⎧ε ⎪⎨ε
x y
⎪⎩τ xy
⎫ ⎪ ⎬ ⎪ ⎭
=
⎢ ⎢ ⎢ ⎢ ⎢
∂x 0 ∂
⎢⎣ ∂y
⎤
0⎥
∂
∂y ∂
⎥⎥⎥⎥⎩⎨⎧uv ⎥
⎫ ⎬ ⎭
=
[B]⎩⎨⎧uuEI
⎫ ⎬ ⎭
∂x ⎥⎦
几何矩阵
有限单元法
土木工程学院
⎧∂ ⎫
⎪ 0⎪
{ }uI = {α1 α 2 α 3 α 4 }T
I for internal , E for external
（1-1）所定义的单元位 移场：
⎧u ⎨⎩v
⎫ ⎬ ⎭
=
⎡ ⎢ ⎣
N1 0
0 N1
N2 0
0 N2
N3 0
0 N3
N4 0
0 N4
(1− ξ 2 ) 0
(1−η 2 ) 0
0 (1− ξ 2 )
的“内自由度”。在计算单元变形能和单元体积力做功时计入这些 位移；但在计算边界外力做功（为了将边界力化为等效节点力）时 不计这些位移。即在计算边界外力做功时只计 Niui、Nivi 各项。
（2）补充这些项后，单元内的位移场是 ξ，η 的完全二次多项式。当实际 单元 e 为矩形时，单元内位移场将是 x、y 的完全二次多项式。
m j =1
ej
(( ∂u )2 ∂x
+ (∂u )2 )dxdy ∂y
（2-6）
能量模：
u= h
Dh (u,u)
(2-7)
3. 有限元解
有限单元法
土木工程学院
总势能函数
π
h P
=
1 2
Dh (u,u) − ( f
,u) − ( p,u)
(2-8)
当所有边界条件均为位移边界条件时：
π
h P
=
1 2