有限元2-弹性力学平面问题有限单元法(2.1三角形单元,2.2几个问题的讨论)综述

第2章 弹性力学平面问题有限单元法2.1 三角形单元(triangular Element)三角形单元是有限元分析中的常见单元形式之一，它的优点是:①对边界形状的适应性较好，②单刚形式及其推导比较简单,故首先介绍之。

一、结点位移和结点力列阵设右图为从某一结构中取出的一典型三角形单元。

在平面应力问题中，单元的每个结点上有沿x 、y 两个方向的力和位移，单元的结点位移列阵规定为： 相应结点力列阵为: (式2-1-1)二、单元位移函数和形状函数前已述及，有限单元法是一种近似方法，在单元分析中，首先要求假定（构造）一组在单元内有定义的位移函数作为近似计算的基础。

即以结点位移为已知量，假定一个能表示单元内部（包括边界）任意点位移变化规律的函数。

构造位移函数的方法是：以结点(i,j,m)为定点。

以位移(u i ,v i ,…u m v m )为定点上的函数值，利用普通的函数插值法构造出一个单元位移函数。

在平面应力问题中，有u,v 两个方向的位移,若假定单元位移函数是线性的,则可表示成:(,)123u u x y x y ααα==++546(,)v v x y x y ααα==++ (2-1-2)a式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标){}⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=m j i m ed d d d m j j i v u v u v u i {}ii j j m X Y X (2-1-1)Y X Y iej m m F F F F ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪==⎨⎬⎨⎬⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎪⎪⎩⎭确定。

将3个结点坐标(x i,y i ),(x j,y j ),(x m,y m )代入上式得如下两组线性方程:123i i i u x y ααα=++123j j j u x y ααα=++ (a)123m m m u x y ααα=++和546i i i v x y ααα=++546j j j v x y ααα=++ (b)546m m m v x y ααα=++利用线性代数中解方程组的克来姆法则,由(a)可解出待定常数1α 、2α 、3α :11A Aα=22A Aα=33A Aα=式中行列式:1i i i j j j m m m u x y A u x y u x y =2111i i j j m mu y A u y u y =3111i i j jm mx u A x u x u =2111i i j j m mAx y A x y x y ==A 为△ijm 的面积,只要A 不为0,则可由上式解出:11()2m m i ij j a u a u a u A α=++ 21()2m m i ij j bu b u b u A α=++ （C ）31()2m mi i j j c u c u c u A α=++式中:m m i j j a x y x y =- m m j i i a x y x y =- m i j j i a x y x y =-m i j b y y =- m j i b y y =- m i j b y y =- （d ）m i j c x x =- m j i c x x =- m j i c x x =-为了书写方便，可将上式记为: m m i j i a x y x y =- m ij by y =- (,,)i j m m i jc x x =-(,,)i j m 表示按顺序调换下标,即代表采用i,j,m 作轮换的方式便可得到(d)式。

将(c)式代入2-1-2中,整理后可得:(,)(,)(,)m m i i j ju N x y u N x y u N x y u =++同理: (,)(,)(,)m m i i j jv N x y v N x y v N x y v =++ （2-1-2）b 式中: 1()2ii i i N a b x c y A=++ (,,)i j m (2-1-3)将三角形单元的位移函数用矩阵表示：或：4)a -1-(2v u i v u v u m 0 j 0 i 0 0 m 0 j 0 ),(m j i i ⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎥⎦⎤⎢⎢⎣⎡=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=m N N N N N N v u f i y x 4)b-1-(2 }]{[}{e d N f v u =⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=三、单元的应变和应力１、应变──几何矩阵[B]由弹性力学知,弹性力学平面问题的几何方程： xu xε∂=∂； y v y ε∂=∂ ； xy u v y x ε∂∂=+∂∂用矩阵表示00x y xy x u v y y x εεε⎡⎤⎢⎥⎢⎥⎧⎫⎢⎥⎪⎪⎧⎫⎢⎥⎪⎪⎪⎪⎨⎬⎨⎬⎢⎥⎪⎪⎪⎪⎢⎥⎩⎭⎪⎪⎢⎥⎩⎭⎢⎥⎢⎥⎢⎥⎣⎦∂∂∂=∂∂∂∂∂或， {}{}H f ε⎡⎤⎣⎦= (2-1-5)[H]称为微分符矩阵,又称为微分算子,“[H]{f}”实际上不是一般的矩阵乘,可以称为微分符矩阵[H]作用在{f}上,其作用规律符合矩阵乘积规定，实际上是按[H]对{f}求导。

将2-1-4式的{f}=[N]{d}代入: {ε}=[H][N]{d}=[B]{d} 2-1-6式中: (2-1-7)000000000000012m i j m i j m i j m i j m m ii j j xN N N B H N y N N N y x b b b c c c c b c b c b A ⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥⎢⎥⎡⎤⎡⎤⎡⎤⎢⎥⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦∂∂∂==∂∂∂∂∂=称为几何矩阵，对于上述三角形单元，[B]是常量矩阵,因此常把这种三角形单元称为常应变单元２、应力矩阵[S]由弹性力学知,由应力求应变的物理方程是:()1x x y Eσμσε=-()1y y x Eσμσε=-2(1)xy xy xyEGτμγτ+==由上式解出应力,得到由应变表示应力的物理方程: ()21x y x E εμεσμ=+-()21y x y E εμεσμ=+-2121xyxy E μγτμ-=∙- 用矩阵表示:{}{}x y xy D σσσετ⎧⎫⎪⎪⎪⎪⎡⎤⎨⎬⎢⎥⎣⎦⎪⎪⎪⎪⎩⎭==21010(1)0021D E μμμμ⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥⎣⎦⎢⎥⎢⎥⎢⎥⎣⎦=-- 2-1-8称为弹性矩阵。

将2-1-6式{ε}=[B]{d}代入上式得: {}{}{}{}eeD D B d S d σε⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎣⎦=== 2-1-9 式中:221111*********m m i i j j m m i i j j m m i i j j S D B bc b c b c E b c b c b c A c b c b c b μμμμμμμμμμμμμ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎛⎫⎢⎥⎪⎝⎭⎢⎥⎢⎥⎣⎦==-------2-1-10上式称为应力矩阵,是结点位移与应力之间的关系矩阵,在上述三角形单元中它也是一个常量矩阵(常应力单元)。

四、单元刚度矩阵有了几何矩阵[B]和弹性矩阵[D]后,我们便可将其代入在§1-3中推导出的单元刚度矩阵的一般表达式: T vK B D B dv ⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎣⎦=⎰对于平面问题，积分dz t =⎰，是单元的厚度，并假定t 在单元内不变化(常数)，所以三角形单元的单刚:T K t B D B dxdy ⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎣⎦=⎰⎰积分式内的两个矩阵都是常量,矩阵乘后,积分得三角形单元的 单元刚度矩阵：ii ij im ji jj jm mm mi mj K K K K K K K K K K ⎡⎤⎢⎥⎢⎥⎡⎤⎣⎦⎢⎥⎢⎥⎣⎦= 2-1-11子块:21122114(1)22r s r s r s r s rs r s r s r s r s b b c c b c c b Et K A c b b c c c b b μμμμμμμ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦--++=---++ (r,s=i,j,m)上式常称为各向同性常应变单刚(stiffness matrix for isotropic constant strain triangle in plane stress)上面的[K]是在图1及2-1-1式的位移排列顺序和(2-1-2)a 下导得的,对此作些改变,便会获得不同的单元刚度矩阵形式,如将位移列阵的排列改为: {d}=[ui uj um vi vj vm ]T 。

（可自行推导,此种单元刚度矩阵的显式可参见《Finite Element Analysis Fundamentals 》）作业3：1．求形函数N i (x,y)在三角形形心（x c ,y c ）上的函数值。

2．设图中i 点有水平位移u i =1，试由单元刚度方程写出 各链杆的反力；并证明各水平、竖向反力之和为0。

3．求图示单元1、2的单元刚度矩阵和应力矩阵。

（结论：将单元逆时针转动180度，则单刚无影响而应力矩阵反号）2.2 三角形单元中几个问题的讨论一、形函数的物理意义 由(2-1-4)a{}{}000(,)000i i em i j j m j i j m m u v N N N u u f x y N dv v N N N u v ⎧⎫⎪⎪⎪⎪⎪⎪⎡⎤⎧⎫⎪⎪⎢⎥⎪⎪⎪⎪⎡⎤⎨⎬⎨⎬⎢⎥⎢⎥⎣⎦⎪⎪⎪⎪⎢⎥⎩⎭⎣⎦⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭==可以看出,当u ι =1,其它5个位移分量为0时,u(x,y)=Ni(x,y),或当v i =1其余为0时 v (x,y)=N i (x,y)(如图所示)。

故形函数Ni(x,y)表示当结点i 发生单位位移时,在单元内部产生的位移分布状态,函数Nj ,Nm 亦具有类似的性质,因此,Ni ,Nj ,Nm 称为位移的形状(态)函数,简称形函数,[N]即称形函数矩阵（shape function ）。

二、形函数的几何意义在单元分析中，很重要的一步是构造位移函数，得出以形函数和结点位移乘积表示的单元位移场。

{}{}000000i i e m ij j m j i j m m u v N N N u f N d v N N N u v ⎧⎫⎪⎪⎪⎪⎪⎪⎡⎤⎪⎪⎢⎥⎪⎪⎡⎤=⎨⎬⎢⎥⎢⎥⎣⎦⎪⎪⎢⎥⎣⎦⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭= 其中形函数: 12i i i i N a b x c y A ⎛⎫ ⎪⎝⎭=++ (),,i j m 而m i j b y y =-x m i j c x =-将其代回:jm m j i y x y x a -=1()()21121m m m m i j i j j j j m mN x y x y y y x x x y A x y x y Ax y ⎡⎤⎢⎥⎣⎦=-+-+-= 若从单元内任意点p(x,y)向各顶点引连线,将其分成三个小三角形,则上式中的行列式恰好是小三角形pjm 面积的2倍即: 1211221i x y AA i N x y i j j A A Ax y m m=== 同理可得:j j A N A = m m A N A=由此可得如下结论(三角形单元的几何性质)1. 任意一点形函数之和等于1,（ Ni +Nj +Nm =1）2. 形函数为≥0,且 ≤1的值, 0≤(Ni ,Nj ,Nm )≤13. 顶点坐标上的形函数值:当(x,y)坐标取在i 点时, Ni =1, Nj =Nm =0 当(x,y)坐标取在j 点时, Ni =1, Nj =1, Nm =0 当(x,y)坐标取在m 点时, Ni =Nj =0, Nm =1因此:Ni ,Nj ,Nm 又称为面积坐标面积坐标的概念在讲叙六结点三角形单元时将得到应用。