《弹性力学与有限元》第6章平面问题的高次单元

（6-5）
要说明的是，矩形单元的插值函数仍然满足第五章所讲的性质，即
⎧1 N i ( xi , y j ) = ⎨ ⎩0
i=j i≠ j
∑N
i =1
4
i
=1
根据平面问题，单元的应变为
⎫ ⎧ ∂u ⎫ ⎧ ∂N i ui ⎪ ⎪ ⎪∑ ⎪ ∂x ⎪ ⎪ ∂x ⎪ ⎪ ⎪ ∂v ⎪ ⎪ ∂N i ⎪ e vi {ε } = ⎨ ⎬ = ⎨∑ ⎬ = [ B ]{a} ∂y ⎪ ∂y ⎪ ⎪ ⎪ ⎪ ∂u ∂v ⎪ ⎪ ∂N i ∂N i ⎪ ui + ∑ vi ⎪ ⎪ + ⎪ ⎪∑ ∂y ∂x ⎭ ⎩ ∂y ∂x ⎭ ⎩
（6-3）
（6-4）
热动 3 同学修改
王正伟 13601363209
其中
ξ = ,η =
x a
y b
写成矩阵形式：
{u} = ⎨
T ⎧u ⎫ ⎡ N i 0 N j 0 N m 0 N p 0 ⎤ ⎥ ( ui , vi , um , vm , u p , v p ) ⎬=⎢ ⎢0 Ni 0 N j 0 N m 0 N p ⎦ ⎥ ⎩v ⎭ ⎣
⎧ui = α1 + α 2 a + a3b + abα 4 , ⎪u = α − α a + a b − abα , 1 3 4 2 ⎪ j ⎨ ⎪um = α1 − α 2 a − a3b + abα 4 , ⎪u p = α1 + α 2 a − a3b − abα 4 , ⎩
vi = α 5 + aα 6 + α 7b + abα 8 v j = α 5 − aα 6 + α 7 b − abα 8 v m = α 5 − aα 6 − α 7 b + abα 8 v p = α 5 + aα 6 − α 7b − abα 8
(6-16)
利用数学知识可以证明： 当求面积坐标的幂函数在三角形单元上的面积分值时，可以用下列积分公式：
∫
V
β γ Lα i L j Lm dV =
α ! β !γ ! 2tA (α + β + γ + 2)!
(6-17)
当求面积坐标的幂函数在三角形单元某一边上的线积分值时，可应用下列积分公式：
∫
Sij
[K ]
e
= ∫∫ [ B ] [ D ][ B ] tdxdy
T
（6-9）
[K ]
e
⎡ kii ⎢ ⎢ k ji =⎢ k ⎢ mi ⎢ ⎣ k pi
T
kij kim kip ⎤ ⎥ k jj k jm k jp ⎥ kmj kmm kmp ⎥ ⎥ k pj k pm k pp ⎥ ⎦
(6-10)
6.1 矩形单元
6.1.1 单元位移模式 图 6-1 是一个边长为 2a 和 2b 的矩形单元，把坐标原点放在单元形心，取单元的位移函数为
⎧u = a1 + a2 x + a3 y + a4 xy ⎨ ⎩v = a5 + a6 x + a7 y + a8 xy
（6-1）
图 6-1 矩形单元 ： 按第五章的方法，把 i、j、m、p 坐标代入式（6-1）
因此 Li 、 L j 、 Lm 中只有两个是独立的变量。 ② 根据定义，面积坐标是一种大小不超过 1 的无因次坐标。当 p 在结点 l 时， Ll ＝1，
(l = i, j , m) ；当 p 在结点 l 所对的边线上时 Ll =0 ；
③ 当 p 在与 ij 边平行的直线上移动时， Am = const ， Lm = const ，因此称这些直线为等
矩形单元的应力为：
(6-7b)
{σ }
e
= [ D]{ε } = [ D][ B]{a} = [ S ]{a}
e e
e
(6-7c)
⎡ ⎤ ⎢1 µ 0 0 ⎥ ⎥ E0 ⎢ [ D] = ⎢ µ0 1 0 ⎥ 1 − µ0 ⎢ 1-µ 0 ⎥ ⎢0 0 ⎥ 2 ⎦ ⎣
计算平面应力时：
(6-7d)
E0 = E 、 µ 0 = µ
[B]的具体数值如
（6-7a）
b-y 0 ⎤ ⎡b + y 0 -(b+y) 0 -(b-y) 0 1 ⎢ a -x 0 a+x 0 0 -(a-x) 0 -(a + x ) ⎥ [ B] = ⎥ 4ab ⎢ ⎢ ⎣ a + x b + y a -x -(b+y) -(a-x) -(b-y) -(a + x) b-y ⎥ ⎦
（6-2）
，整理得 由 6-2 式解得 α1 ~ α 8 ，代入式（6-1）
⎧ ⎪u = N i ui + N j u j + N mum + N p u p ⎨ ⎪ ⎩v = N i vi + N j v j + N m vm + N p v p
这里
1 ⎧ ⎪ N i = 4 (1 + ξ )(1 + η ) ⎪ ⎪ N = 1 (1 − ξ )(1 + η ) ⎪ j 4 ⎨ ⎪ N = 1 (1 − ξ )(1 − η ) ⎪ m 4 ⎪ 1 ⎪ N p = (1 + ξ )(1 − η ) ⎩ 4
对 6-12 式整理得：
热动 3 同学修改
王正伟 13601363209
⎧1 ⎫ ⎡1 1 1 ⎤ ⎧ Li ⎫ ⎥⎪ ⎪ ⎪ ⎪ ⎢ ⎨ x ⎬ = ⎢xi x j xm ⎥ ⎨Lj ⎬ ⎪ ⎪ ⎢ ⎥⎪ ⎪ ⎩ y ⎭ ⎣ yi y j y m ⎦ ⎩ Lm ⎭
这就是直角坐标用面积坐标表示的形式。 6.2.3 面积坐标的微积分运算 建立了坐标系间的变换关系，则按求导法则可得：
下面我们讨论如何构造 [ N ] ，这里则采用了一种简便的求形函数的方法，它根据
Ni = δ iP 的条件用几何方法构造 Ni 函数，然后再用位移模式的完备性和协调性条件作为校
核。 假设单元内具有 d 个结点， P 为任意一个结点号 （P=1， 2， …， i， …， d） 。 欲求 N i ( x, y ) ， 可做一组（ m 条）不通过 i 结点而通过其它所有结点的不可约代数曲线 Lk ( x, y ) = 0
第六章
平面问题的高次单元
在第五章， 我们介绍了有限元法中最简单的平面单元----常应变的三角形单元。 实际上 一个物体的应力场是随坐标变化的， 甚至有些地方变化急剧。 如果采用常应变三角形单元分 析，则必须划分很多单元才能得到较高的精度，这样势必需要求解庞大的方程组。本章我们 将介绍几种平面问题的高次单元， 在单元内部位移函数采用二次以上多项式， 应力应变都是 沿坐标变化的，这样，即使采用较少单元也能得到较好的计算结果。
β Lα i L j dS =
α !β ! tlij (ij → jm → mi ) (α + β + 1)!
（6-18）
式中 α 、 β 、 γ 为面积坐标的幂， lij 为 ij 边的长度， t 为平面问题板厚。
6.3 六节点三角形单元
6.3.1 单元位移函数 1、 六节点三角形单元的位移函数 三节点三角形是常应变单元，它的插值函数是坐标的一次函数，为了提高精度，下面我 们就讨论六节点的三角形单元。 图 6-3 是一个六节点三角形单元。6-3a 是网格中的六节点三角形单元，6-3b 是六节点三 角形单元的编号方式和各节点的面积坐标。 单元的位移模式可以仿照三节点形式的构造， 即
坐标线，且它的大小是从 0 到 1。 6.2.2 面积坐标与直角坐标之间的关系 设 p 点坐标为（x，y） ，则由数学可知，三块小面积为：
1 x y 1 1 1 Ai = 1 x j y j = ⎡ x j ym − y j xm ) + ( y j − ym ) x + ( xm − x j ) y ⎤ = ( ai + bi x + ci y ) ⎣ ⎦ 2 2 2 1 x m ym
（6-15）
∂ ∂ ∂Li ∂ ∂L j ∂ ∂Lm 1 ⎛ = + + = ⎜ bi ∂x ∂Li ∂x ∂L j ∂x ∂Lm ∂x 2 A ⎜ ⎝ ∂ ∂ ∂Li ∂ ∂L j ∂ ∂Lm 1 ⎛ = + + = ⎜ ci ∂y ∂Li ∂y ∂L j ∂y ∂Lm ∂y 2 A ⎜ ⎝
∂ ∂ ∂ ⎞⎫ + bj + bm ⎟⎪ ∂Li ∂L j ∂Lm ⎟ ⎠⎪ ⎬ ∂ ∂ ∂ ⎞⎪ + cj + cm ⎟⎪ ∂Li ∂L j ∂Lm ⎟ ⎠⎭
⎧ Li ⎫ ⎡ a i bi ci ⎤ ⎧1 ⎫ ⎪ ⎪ 1 ⎢ ⎥⎪ ⎪ a j bj cj ⎥ ⎨x ⎬ ⎨Lj ⎬ = ⎢ ⎪ ⎪ 2A ⎢ ⎥⎪ ⎪ ⎣a m bm cm ⎦ ⎩ y ⎭ ⎩ Lm ⎭
由于
（6-13）
1 x i yi 1 A = 1 x j yj 2 1 x m ym
(6-14)
图 6-2 分别记这三个小三角形面积为
Ai = Apjm
Aj = Apmi
Am = Apij
则显然存在如下恒等关系： A = Ai + A j + Am
热动 3 同学修改
王正伟 13601363209
若令
Ll =
Al A
(l = i, j , m)
（6-11）
因此 Ll (l = i, j , m) 那么 p 点的坐标就可以用量纲为 1 的参数 Li 、L j 、Lm 中的两个来确定， 可以作为确定点位置的一种坐标。因为 Ll 是根据面积来定义的，故称 Li 、 L j 、 Lm 为面积 坐标。 由上述定义可见：面积坐标是一种固定于单元的局部坐标，它具有如下性质： ① 这三个面积坐标并不是完全独立的，由于 A = Ai + A j + Am ，所以有 Li + L j + Lm = 1 ；
对于任一子块
[ K rs ] = ∫∫ [ Br ] [ D ][ Bs ] tdxdy
(r , s = i，j，m，p)
6.2 面积坐标及其应用
6.2.1 面积坐标的定义 思考我们讨论三节点三角形单元时， 如果三角形的边上节点数增加， 如何计算单元位移 模式中的 α i ？事实上在计算等效节点载荷时，在 xoy 坐标系就显得比较复杂。那么有没有 可能选择一种坐标，使计算简化一些？答案时肯定的，下面我们引入面积坐标。 设 p 为三角形中的一点，与三个顶点相连将 ∆ijm 分割成 3 小块（图 6-2） 。
u = ∑ N i ui
i =1 6
6
v = ∑ N i vi
i =1
（6-19）
热动 3 同学修改
王正伟 13601363209
图 6-3
{u} = ⎨
⎧u ⎫ e ⎬ = [ N ]{a} ⎩v ⎭
⎡ N1 0 N 2 0 N3 0 N 4 0 N 5 0 N 6 0 ⎤ T =⎢ ( u1 , v1 , u2 , v2 , u3 , v3 , u4 , v4 , u5 , v5 , u6 , v6 ) ⎥ ⎣0 N1 0 N 2 0 N 3 0 N 4 0 N 5 0 N 6 ⎦
(6-6)
⎡ ∂N i ⎤ ∂N j ∂N p ∂N m 0 0 0 0 ⎥ ⎢ ∂x ∂x ∂x ⎢ ∂x ⎥ ⎢ ∂N j ∂N p ⎥ ∂Ni ∂N m 0 0 0 [ B] = ⎢ 0 ⎥ ∂y ∂y ∂y ∂y ⎥ ⎢ ⎢ ∂N ∂N ∂N j ∂N j ∂N ∂N ∂N p ∂N p ⎥ i m m ⎢ i ⎥ ∂x ∂y ∂x ⎦ ⎢ ∂y ∂x ∂y ∂x ∂y ⎥ ⎣
热动 3 同学修改
王正伟 13601363209
计算平面应变时：
源自文库
E0 =
E µ 、 µ0 = 2 1− µ 1− µ
(6-8)
[ S ] = [ D][ B]
[B]、 [S]都不是常数矩阵， 而是坐标 x、 y 的函数， 就是说，σ x、σ y、τ xy 及 ε x、ε y、γ xy 不再是常量，是 x、y 的一次函数。 6.1.2 单元刚度矩阵 利用最小势能原理，对平面问题 矩形单元刚度矩阵展开为：
( k = 1, 2, ⋅⋅⋅, m ) ，并按下式确定 Ni ( x, y ) ：
N i ( x, y ) =
∏ L ( x, y ) ∏L (x , y )
k =1 k i i k =1 m k
m
注意，当结点 P ≠ i 时，P 结点是位于某条（设其编号为 k）上述所作代数曲线上，故 有 Lk ( x, y ) = 0 ； 而 i 结点是不通过上述一条代数曲线的， 故 Lk ( x, y ) ≠ 0 由此可知，可以满足条件：
脚标轮换规则为 i → j → m → i 因此：
Li =
Ai 1 = ( ai + bi x + ci y ) A 2A
（6-12）
不难看出，6-11 式的右边就是前面讲过的 3 节点三角形的插值函数 N i , 可见，面积坐标 Li , Lj, Lm 就是 Ni, Nj, Nm。即 Li = N i , L j = N j , Lm = N m ，将 6-11 表示成矩阵形式是：