1[1].13非协调元

表 13.1 悬臂梁计算结果
66
图 13.3 二次非协调位移模式
( ) 和 内部自 由 度 相 关的 附 加 位 移项 在 单 元 中引 起 的 位 移表 示 在 图 13.3 中。 α1 1− ξ 2 和
( ) ( ) ( ) α3 1− ξ 2 项在单元η = ±1 的边界上呈二次抛物线变化。 α2 1−η 2 和 α4 1−η 2 项则在单元
( ) ( ) 4
∑ v = Nivi + α 3 1 − ξ 2 + α 4 1 −η 2 i =1
其中
Ni
=
1 4
(1
+
ξ
i
ξ
)(1
+
η
iη
)
(i = 1,2,3,4)
(13.4)
α1 ,α2 ,α3 ,α4 为内部自由度。
将(13.4)式表示成矩阵形式
{u} = [N ]{δ}e + ⎡⎣N ⎤⎦{α}e
图 13.2 悬臂梁网格划分及载荷情况 现以图 13.2 所示悬臂梁为例，说明 4 结点单元增加内部自由度后对精度的改进。对于两种不同
载荷和网格划分，采用两种单元的计算结果列于表 13.1。在现在的载荷和网格化分情况下，原线性 单元，即无内部自由度单元的计算精度是相当差的。将改进后的，即具有内部自由度的单元和精确 解比较，可以看到计算精度有显著的提高。
为了检验采用非协调元的任意网格划分时能否达到上述连续性的要求，可以进行分片试验。若 能通过分片试验，则解的收敛性就能得到保证。
二、分片试验
图 13.4 单元片 考虑一任意的单元片如图 13.4 所示，其中至少有一个结点是完全被单元所包围的如图中的结点 i，结点 i 的平衡方程为
( ) ∑ { } m ⎣⎡k界上呈二次抛物线变化。α1 ,α2 ,α3 ,α4 是在单元内部定义的内部自由度，所以这些附加
位移项在单元与单元的交界面上是不保证协调的，也就是说由于单元内增加了附加位移项而致使单 元之间不能保证在交界面上位移的连续性。这些附加位移项称之为非协调项。引入非协调位移项的 单元称为非协调元。可以证明，对于 C。型问题，如果在单元尺寸不断缩小的极限情况(应变趋于常 应变)下，位移的连续性能得到恢复，则非协调元的解仍然趋于正确解。因此问题转换为检验非协调 元是否能描述常应变，以及在常应变条件下能否自动地保证位移的连续性。
{Ru}e = ∫ [N ]T { p}dV + ∫ [N ]T {q}dS
Ve
Sσe
{Rα }e = ∫ ⎡⎣N ⎤⎦T { p} dV + ∫ ⎡⎣N ⎤⎦T {q} dS
Ve
Sσe
从(13.7)式的第 2 式可以解出
(13.8)
( ) ( ) {α}e = [kαα ]e −1 [Rα ]e − [kαu ]e {δ }e
u = α1xy, v = 0
(13.3)
所以(13.1)式的第二式位移 v 的表达式实际上表示了利用一个双线性单元模拟纯弯曲应力状态时所
出现的误差。由(13.3)式表示的位移模式计算得到的剪应力τ xy 以及正应σ y 表示在图 13.1(d)和(e)。
导致误差的原因是缺少完全的二次多项式，即缺少 x2 项和 y2 项。
上式消去(13.7)式第 1 式中的内部自由度{α} ，则得到凝聚后的单元求解方程
(13.9)
[k ]e {α}e = {R}e
(13.10)
( ) 其中
[k ]e = [kuu ]e − [kuα ]e
[kαα ]e
[ ] k −1
e
αu
( ) {R}e = {Ru}e − [kuα ]e [kαα ]e { } −1 Rα e
(13.5)
其中：
64
{u}
=
⎧u ⎨⎩v
⎫ ⎬ ⎭
⎧u1 ⎫
{δ }e
=
⎪⎪ ⎨ ⎪
v1 #
⎪⎪ ⎬ ⎪
⎩⎪v4 ⎭⎪
⎧α1 ⎫
{α}e
=
⎪⎪α ⎨⎪α
2 3
⎪⎪ ⎬ ⎪
⎪⎩α4 ⎭⎪
[N ] = ⎡⎣[N1] [N2 ] [N3 ] [N4 ]⎤⎦
代入几何关系可以得到
⎡⎣
N
⎤⎦
=
⎡1 ⎢ ⎣
−ξ 0
e=1
(13.11)
其中 m 是单元片包含的单元数，结点 j 代表单元片内除 i 结点外的所有结点。 分片试验是指：当赋予单元片各个结点以与常应变状态相应的位移值和载荷值时，校验(13.11)
式是否满足。如能满足，则认为通过分片试验，也就是单元能满足常应变要求，因此当单元尺寸不 断缩小时，有限元解能够收敛于真正解。
ξ ,η ,ξη ,ξ 2 ,η 2 ,ξ 2η ,ξη 2 。这些插值函数中所包含的完全多项式分别只是一次的和二次的，它们所
要求的自由度分别是 3 和 6，即只需要单元的结点数是 3 和 6。就构成决定单元精度的完全多项式而 言，其有效的结点数分别只是 3 个和 6 个。从这个意义上说，二维等参元中有四分之一的结点自由 度对计算精度是不起作用的。另一方面，插值函数中非完全的高次项一般说来非但对改善精度不起 作用，而且还可能起相反的作用。所以从这个意义上来讲，等参元的精度在给定自由度的条件下是 不够理想的。上述缺点对三维等参元来说将更突出。因为在三维坐标中，一次完全多项式是 4 项：1，
第十三节 非协调元和分片试验
等参元有良好的适应性和表达格式的简明性，因而得到广泛的应用，但是从严格的意义上说， 它的精度和效率仍是不够高的。以二维单元为例：双线性单元有 4 个结点，对应的插值函数中包含
下 列 4 项 ： 1 ， ξ ,η ,ξη ； 二 次 单 元 有 8 个 结 点 ， 对 应 的 插 值 函 数 中 包 含 下 列 8 项 ： l ，
2
1−η2 0
0 1−ξ 2
0⎤
1
−η
2
⎥ ⎦
{ε} = [B]{δ}e + ⎡⎣B⎤⎦{α}e
将{u} ,{ε} 等代入位能泛函 Π 并按照通常的步骤，泛函变分为零可以得到
⎢⎡[kuu ]e ⎢⎣[kαu ]e
[kuα [kαα
]e ]e
⎤ ⎥ ⎥⎦
⎧⎪{δ }e ⎩⎨⎪{α}e
⎫⎪ ⎬ ⎭⎪
=
⎩⎧⎪⎨⎪{{RRαu
e=1
(13.14)
如果上式不成立，说明当单元片各结点具有与常应变状态相应的位移时，结点 i 不能保持平衡。 必须在结点 i 施加相应的外载荷，即给 i 结点以某种约束，平衡才能维持。这就说明这种非协调元 不能满足常应变的要求。为满足常应变的要求就必须施加必要的约束力，该约束力所作的功等于单 元交界面上位移不协调而引起的附加应变能。
分片试验的另一提法是：当分片的边界结点赋予和常应变相应的位移值时，求解方程(13.14)式，
得到分片的内部结点 i 的位移值{δi} ，如{δi} 和常应变状态相一致，则认为通过分片试验。如{δi}
与常应变状态不一致，则认为通不过分片试验。通常认为前一种提法应用更方便些，因为后一种提 法要涉及矩阵求逆的计算。
例如在平面问题中，与常应变相应的位移是线性位移。
u = α1 + α2 x + α3 y v = α4 + α5x + α6 y
赋予各结点以与常应变状态相应的位移，即令结点位移为
(13.12)
67
u j = α1 + α2 x j + α3 y j v j = α4 + α5x j + α6 y j
现在来研究 Wilson 的 4 结点平面非协调元通过分片试验的条件。这种非协调元的位移插值表示式
是(13.4)式。我们知道，当单元的位移插值不包含非协调项时(即α1 = α2 = α3 = α4 = 0 )，是满足
收敛条件的，当然也必定通得过分片试验。因此，当单元片各结点赋予与常应变相应的(13.13)式的
(13.13)
由平面问题的方程可知，与常应变(也相当于常应力)状态相应的载荷条件应有体积力为零，同时
在 i 结点上也不能作用集中力(包括面积力的等效力)，因此{Ri}e 也必须为零。所以此时，通过分片
试验的要求是：当赋予各结点以(13.13)式所示位移时，下式仍成立。
∑ { } m ⎡⎣kij ⎤⎦e δ j = 0
一、 Wilson 非协调元
为了改善二维线性单元的性质，提高其精度，Wilson 提出在单元的位移插值函数中附加内部无结
( ) ( ) 点的位移项。当单元是等参元，采用自然坐标时，此附加项为α1 1− ξ 2 和α2 1−η 2 。从形式上
看，这两项和(13.1)式第二式所包含的项次相同。而它们正是利用二维双线性单元模拟纯弯应力状
( ) {δ}e = {δ}e l
因此，通过分片试验的要求是
r
(13.15) (13.16)
( ) ( ) ∫ {α}e =
−
[kαα ]e
−1
⎛ ⎜⎜⎝
Ve
⎡⎣
B
⎤⎦T
[
D][
B]
dV
⎞ ⎟⎟⎠
{δ }e
l
=0
(13.17)
( ) [ ] [ ] 又因为 kαα e 是非奇异的(因为附加位移不存在刚体移动)，其逆
ξ ,η ,ζ ；二次完全多项式是 10 项，而三维线性单元和二次单元却分别具有 8 个和 20 个结点，也即
三维等参元中有二分之一的结点自由度对计算精度是无贡献的。因此，E．wilson 提出的二维和三 维的非协调等参元，对改进等参元的计算精度和提高计算效率是很有意义的，特别是对于三维问题 的有限元分析。
它对结点位移没有影响，而只对单元内部的位移起了调整作用。确定附加项前的待定系数α1 和α2 与
单元边界上的结点无关，这种仅在单元内部定义的待定参数称为内部自由度。 包含附加的无结点位移项的单元位移插值表示如下：
( ) ( ) 4
∑ u = Niui + α1 1 − ξ 2 + α 2 1 −η 2 i =1
α
1
a2
−
x2
1 2
α
1
μ
b2
−
y2
(13.1)
(a)弯曲应力σ x (b)实际位移 (c)4 结点矩形元的近似位移 (d)(e)有限元解的剪应力τ 及σ y
图 13.1 受纯弯作用的矩形单元
63
利用平面问题的几何关系和物理关系可以得到
σ x = α1Ey
σ y = τ xy = 0
(13.2)
单元处于纯弯应力状态，从而证明了(13.1)式的位移模式是精确满足纯弯应力状态的. 如果我们用一个线性矩形单元去模拟上述受力状态，得到的位移将如图 13. 1(c)所示,即
65
此即包含附加内部位移项的单元刚度矩阵和载荷列阵。它是在原单元刚度矩阵和载荷列阵内增加了 修正项而得到。消去内部自由度以及修正单元刚度矩阵和载荷列阵都是在单元分析过程中进行的， 此过程称为内部自由度的凝聚。经凝聚后，单元的自由度仍是原四边形单元的自由度，以后的分析
和计算步骤也和标准的解题步骤相同。顺带指出：在通常不存在体积力(即， { p} ≡ 0)的情况下， {Rα }e 中的第 2 项也一并略去，可减少计算工作量，而计算结果表明，对精度没有什么影响。
现在我们来讨论二维双线性单元在表示纯弯应力状态时出现的问题。由于二维双线性单元它的插
值函数中包含有非完全的二次项 ξη ，因此用它表示纯弯曲应力时，出现明显的误差。图 13.1(a)表
示受纯弯作用的矩形单元，其精确位移解答如图 13.1(b)所示，并可表达如下
u = α1xy
( ) ( ) v
=
1 2
位移值时，应有α1 = α2 = α3 = α4 = 0 。另一方面(13.9)式及(13.8)式可知
( ) ( ) ∫ {α}e = − [kαα ]e −1 [kαu ]e {δ }e
=−
[kαα ]e
−1
⎛ ⎜⎜⎝
Ve
⎡⎣
B
⎤⎦T
[
D][
B]
dV
⎞
⎟⎟⎠{δ
}e
( ) 将与常应变状态相应的结点位移记为 {δ}e ，则有 l
kαα
e
−1
存在。同时还因为和
常应变相应的应力是常应力，即
( ) [D][B] {δ}e = {σ}
l
l
是常应力状态。所以，通过分片试验的要求可以由(13.17)式简化为
(13.18)
∫ ∫ ∫ ⎡⎣B⎤⎦T dV ≡
态造成误差的原因所在。增加附加项后就有可能通过适当调整系数α1 和α2 使误差降到最小。从数
学上看，是通过引入 ξ 2 和η 2 项，使插值函数中的二次式趋于完全，从而达到提高计算精度的目的。
( ) ( ) 可以看到附加项α1 1− ξ 2 和α2 1−η 2 的位移，在二维线性单元的四个结点上都取零值，即
}e }e
⎫⎪ ⎬ ⎪⎭
其中
∫ [kuu ]e = [B]T [D][B]dV (是原 4 结点线性单元的刚度矩阵) Ve
( ) [kuα ]e = [kαu ]e T = ∫ [B]T [D]⎡⎣B⎤⎦ dV Ve
[kαα ]e = ∫ ⎡⎣B⎤⎦T [D] ⎡⎣B⎤⎦ dV Ve
(13.6) (13.7)