常用单元的刚度矩阵
2(ru)2ru
2rr
由于各点在圆周方向上无位移,因而剪应变vr和vr均为
u
rru
零。将应变写成向量的形式,则rw
z
zrzuw
zr
根据上式,可推导出几何方程B(e)
zjk
N(r,z)ir10
其中几何矩阵B
2r
kj
00rkjzjk
zkj
Nj(r,z)r0rik
00rikzki
zij
Nk(r,z)r0rji
00rji zij
3.弹性方程和弹性矩阵[D]
依照广义虎克定律,同样可以写出在轴对称中应力和应变之间的弹性方程,其形式为
1
ru(z) E1
u(rz)
E1
zzu(r)
E
2(1)rrzrz
E
r
所以弹性方程为D 式中应力矩阵rzrzT
1
1E
弹性矩阵D1(1)(12)
000
00 12
2
4.单元刚度矩阵k(e)
与平面问题相同,仍用虚功原理来建立单元刚度矩阵,其积分式为
k(e)BTDBdV
V
在柱面坐标系中,dV将dV
2drdz代入k
(e)
2drdz
T
BDBdV
V
,则k(e)2BTDBrdrdz
即为轴对称问题求单元刚度矩阵的积分式。
与弹性力学平面问题的三角形单元不同,在轴对称问题中,几何矩阵[B]内有的元素(如
Ni(r,z)
等)是坐标r
r、z的函
数,不是常量。因此,乘积BTDB不能简单地从式
k(e)2BTDBrdrdz的积分号中提出。如果对该乘积逐项求
积分,将是一个繁重的工作。一般采用近似的方法:用三角形形心的坐标值代替几何矩阵[B]内的r和z的值。用B表示在形心(r,z)处计算出的矩阵[B]。其中
r
(rirjrk)
3
,z
(zizjzk)
3
只要单元尺寸不太大,经过这样处理引起的误差也不大。被积函数又成为常数,可以提出到积分号外面:
k(e)2BDBrdrdz2BDBr式中——三角形的面积。
T
T
由式k
(e)
2BDBrdrdz2BDBr可以看出,两轴对
T
T
称的三角形单元,当形状、大小及方位完全相同而位置不同时,其刚度矩阵也不相同。距离主轴线越远的单元,其刚度越大。这与平面问题不一样。
二、等参数的刚度矩阵
对一些由曲线轮廓的复杂结构,如果采用直角边单元进行离散,由于用直线代替了曲线,除非网格划分得很细,否则不能获得较高的精度;对另一些应力随坐标急剧变化的结构,采用简单的常应力单元离散时,也必须划分成大量的微小单元,以保证足够的精度。为此引入一种高精度的单元——等参数单元。它既能简化复杂单元划分的工作,又能在满足同样精度的要求时,大大减少使用的单元数。目前流行的大程序中较常用,它成功地解决了许多二维和三维的弹性力学问题。
为导出等参数单元的刚度矩阵,首先要建立根据每个单元的形状确定的自然坐标系,然后将位移模式和形状函数都写成自然坐标的函数。
一个单元在自然坐标系内的点余元整体坐标系内的点成一一对应的关系。通过映射,可以将整体坐标系中的图形转化为自然坐标系中的相应徒刑。例如可以将整体坐标系中的一个任意四边形(实际单元)映射到自然坐标系中成为一
个正方形(基本单元)。同样也可以将任意四面体、六面体(包括直边和曲边的)分别映射成正四面体和正六面体。
这里只介绍较简单的一种平面问题的情况,将整体坐标系中的一个任意四边形映射成自然坐标系中的一个正方体,并导出单元刚度矩阵。其它种单元的映射,可依次原理进行。不再叙述。
1. 位移模式和形状函数
图4-2中的任意四边形单元上,作连接对边中点的直线,取其交点为原点,这两条直线分别为和轴,并令四条边上的和值分别为1,建立一个新的坐标系,称之为该单元的自然坐标系。原坐标系XOY称为整体坐标系。在整体坐标系中,自然坐标系非正交,它由任意四边形的形状所确定。
图4-19
如果将自然坐标系改画成直角坐标系,那么图4-19(a)中的任意四边形单元就成为图4-19(b)所示的正方形。
上述两个四边形的点(包括顶点)一一对应,即它们之间相互映射。因此,需要写出整体坐标X、Y和自然坐标、之间的坐标转换式,即
X1234Y5678
*
四边形四个顶点的坐标值在XOY坐标系中分别为
X1,Y1,X2,Y2,X3,Y3,X4,Y4:在
o
坐标系中相应为
1,1,1,1,1,1,1,1。将有关数据代入*中的第一式,则有
X11234,X21234X31234,X41234
求解上述方程组得:
X1X2X3X4X1X2X3X4
,2
44
X1X2X3X4XX2X3X4
3,41
44
1
坐标变换方程*成为
X
1
1X11X21X31X44
同理
Y
1
1Y11Y21Y31Y44
当引入函数Ni,后,坐标变换方程成为
XNi,Xi
i144
YNi,Yi
i1
式中Ni,11i1i
4
变量、的正负号由相应节点的坐标值i、i决定。例如当i=4时,1,1,因此,N,11。
4
4
4
4
下面再来研究函数Ni,的特性。
对节点1X1,Y1,相应的自然坐标值为(-1,-1)。从式
Ni,
1
1i1i4
中很容易看出,除N1=1外,
N2=N3=N4=0。对其余各节点也一样。总而言之,对节点i(i=1,2,3,4),除Ni=1外,其余三个N值均为零。同时,不难看出
N1,N2,N3,N4,1,即四个节点的
Ni函数
之和等于1。
函数Ni,具备上章所介绍的形状函数应满足的条件,可作为本单元的形状函数。
采用Ni,做形状函数,其位移模式为
uNi,ui,vNi,vi
i1
i1
4
4
对比
XNi,Xi
i14
4
YNi,Yi
i1
和uNi,ui,vNi,vi可以看出:
i1
i1
44
在这种实际单元(任意四边形)中,坐标变换式和位移模式不仅采用了相同的形状函数Ni,,而且具有相同的数学模型。这种性质的实际单元称为等参数单元。
对用节点位移值ui(或vi等)求单元内某一点位移量u(或v等)的插值公式uNi,ui,只要将u(或v等)换成X(或Y
i14
等),便成为利用节点值Xi(或Yi等)求相应点坐标X(或Y等)的插值公式。相反也是这样。 2.几何矩阵[B]
由于几何矩阵[B]通过对位移求偏导数而得出,所以首先必须利用复合函数求导的规则得出下述公式
uuXuuX
XuYuuXY
或写成Ju XuYu
Yy
X
式中JX
Y
,此式称为雅可比矩阵。 Y
为了将几何矩阵[B]写成变量、的函数,必须将式
uu
X
uJu改写成
y
uvuv
XX11
uJu,同理,vJv yy
从表示单元内各点位移与其应变关系的几何方程可知:
X0
Y
0
1000
uu
0001
Yv0110
X
u
Y
vX
vY
T
uvuvXX11
将式uJu和vJv合并,则
yyuX
uY
vX
JvY0
T
1
0u
J1
u
v
v
T
对单元(e),任意一点的位移u,v对自然坐标,的偏导数可利用上式求出,写成矩阵形式为:
u
u
v
v
e Np
T
式中eu1
1p
p
1p
v1u1v1u1v1u1v1
2p
3p
T
4p
4p
N0N0N0N0
N0N0N0N0N
2p
3p
对于i=1,2,3,4 将
u
u
N
ip
Ni
1
Ni
T
u
v
v
T
uXv
uY
T
vX
JvY0
T
0u
J1
和入
v
eNp
代
X0
Y
0
1000uu
0001
Yv0110
X
u
YvX
vY
T
,则可得出表示
在整体坐标系中位移和应变关系的几何方程:eBee 式中的几何矩阵[B]是自然坐标,的函数:
1000
J1
B0001
00110
0
Np1J
T
也可利用Nip
Ni
Ni
求得的Nip以及
XNi,Xi
i14
4
和
YNi,Yi
i1
XJX
Y
X1
求出J,JN1pN2pN3pN4pYY1
X2
Y2X3Y3
X4Y4
T
。
3.单元刚度矩阵ke
设单元板厚为t,根据虚功方程有:keBTDBtdA,
A
此式中几何矩阵[B]和弹性矩阵[D]都已求出。因为几何矩阵[B]中的变量是自然坐标,,所以也要用自然坐标表示微分面积dA。
在实际单元中任取一点p,其整体坐标位X、Y,其相应
的自然坐标为,。过p点做,的等值线,同时做
d,d的等值线,围成一小块微分面积
dA,如图4-20
(a)所示。为便于分析,将四边形pqrs放大,如图4-20(b)所示。实际上,d,d取得很小,因此该四边形可视为平行四边形。若相邻的两边用向量a,b表示,则两者的乘积等于该平行四边形的面积dA。
图4-20
dAabsinab
若
aaxiayj,bbxibyj
axay
bxby
,则
dAaxiayjbxibyj
为了求出ax,ay,bx,by的值,要先写出a和b两端节点p、q、s的坐标值。
点p:XpX,,YpY, 点q:XqXd,,YqYd, 点s:XsX,d,YsY,d 利用泰勒技术展开并略去高阶项,可得
X
d
X
X,dX,d
Xd,X,
对Yd,,Y,d,也可写出相应的展开式。利用式
Xd
可得: X
X,dX,d
Xd,X,
XYd,ayYqYpdXY
bxXsXpd,byYsYpd
axXqXp
,
将此式代入式dAaxiayjbxibyj
ax
dA
ay
XbxXby
axay
bxby
得到:
Y
dd,简写为dAJdd Y
单元刚度矩阵为:
k
e
11
T
,用BDBtJdd,这个积分可以采用“数值方法”
11
高斯求积分公式很方便的求出,在此不作介绍。 例:求如图所示四边形的雅可比矩阵。
解:求雅可比矩阵可在整体坐标系中进行,也可以在实际单元的局部坐标系中进行。为便于计算,本例在局部坐标系中进行。
对单元(1):
将四个节点的自然坐标值(-1,-1)、(1,-1)、(1,1)、(-1,1)代入下式:
Ni,
1
1i1i,再将所得到的Ni,值及四个节点4
实际单元在局部坐标系的坐标值(-3,-2)、(3,-2)、(3,2)、(-3,2)代入下式计算:
XNi,Xi
i1
44,则X3,Y2
YNi,Yi
i1
X雅克比矩阵为JXY30。 Y02
对单元(2):
四个节点在局部坐标系中的坐标值分别为(-1,-3/4)、(1,-3/4)、(1,5/4)、(-1,1/4)。因此,
X1111111111111 4
Y1335111111111 44444因而雅可比矩阵为
J14
401 3
X1Y1X2Y2X3Y3X4Y4T也可利用JN1pN2pN3pN4p求雅可比矩阵,
其结果与上相同,同学们可自行验证。
图
11