现代设计方法4-3 三角形三节点平面单元课件

1、4.3 平面问题的有限单元法,三角形三节点平面单元,1,学习交流PPT,结构离散化,单元分析,整体分析,有限元分析的基本步骤：,2,学习交流PPT,1 结构离散化,例：图示一悬臂梁，梁的厚度为t，设泊松比m =1/3，弹性模量为E，试用三节点三角形单元进行离散。,3,学习交流PPT,2 单元分析,单元分析的主要内容：由节点位移求内部任一点的位移，由节点位移求单元应变、应力和节点力。 单元分析的步骤：,节点 位移,单元内部 各点位移,单元 应变,单元 应力,节点 力,(1),(2),(3),(4),单元分析,4,学习交流PPT,(1) 单元位移模式 有限元法中，通常采用位移法进行计算，即取节点的

2、位移分量为基本未知量，单元中的位移、应变、应力等物理量，都和基本未知量相关联。 节点i的位移分量可写成 di = ui viT 单元节点位移向量d e可写成 d e = di dj dmT = ui vi uj vi um vmT,三角形三节点单元,2.1 由节点位移求单元内部任一点位移,5,学习交流PPT,六个位移分量需六个待定参数 a1，a2，a6 设单元内任一点位移分量是坐标(x,y)的线性函数： (4.22) 写成矩阵的形式为： (4.23),6,学习交流PPT,设节点i，j，m坐标分别是xi，yi；xj，yj；xm，ym。把三个节点的坐标及其水平位移代入式(4.22)中得： (a)

3、解得： (b) 其中，,(2) 由单元节点位移d e求位移参数a ,7,学习交流PPT,对v同理可列出a4、a5、a6的方程 ： vi = a4+a5 xi+a6 yi i, j, m 解得： a4、a5、a6 为了书写的简便与规格化，引入记号 ai, bi, ci分别为行列式A中各行各个元素的代数余子式,8,学习交流PPT,解出a1 a6结果：,a =Ld e,为三角形单元面积。,9,学习交流PPT,将a写成矩阵形式，有a =Ld e,由单元节点位移d e求单元内部任一点位移d (x,y) d (x,y)=m(x,y)a =m(x,y)Ld e,10,学习交流PPT,11,学习交流PPT,N

4、i, Nj, Nm是坐标的连续函数，反映单元内位移的分布状态，称为位移的形状函数，简称形函数(shape function)。矩阵N称为形函数矩阵(shape function matrix) 。 形函数物理意义：Ni(x,y) ，节点i单位位移，其它节点位移分量为0，单元内部产生位移分布形状,12,学习交流PPT,在单元任一点上，三个形函数之和等于1。 形函数Ni在i点的函数值为1，在j点及m点的函数值为零。 三角形单元i, j, m在ij边上的形函数与第三个顶点的坐标无关。,形函数的性质：,13,学习交流PPT,例：求图示单元和单元的形函数矩阵,14,学习交流PPT,（a）,(b),分别如

5、上图所示：,15,学习交流PPT,单元如图所示。设a=1m, b=2m.,(或直接由图形可知其面积),16,学习交流PPT,求系数ai, aj, am, bi, bj, bm, ci, cj, cm,17,学习交流PPT,求形函数矩阵,代入相关常数：,18,学习交流PPT,将a=1,b=2代入得：,19,学习交流PPT,求常数,单元如图所示=ab/2,20,学习交流PPT,求形函数矩阵,21,学习交流PPT,将a=1, b=2m代入上式得：,22,学习交流PPT,作业：,已知三角形三节点单元坐标如图示，设单元中一点A的坐标(0.5,0.2)，已知三角形三节点单元i节点位移(2.0,1.0)，

6、j节点位移(2.1,1.1)， m节点位移(2.15,1.05)， 1)写出单元的位移函数；2)求A点的位移分量。,A,23,学习交流PPT,2.2 由节点位移求单元的应变,几何方程,24,学习交流PPT,25,学习交流PPT,简记为e =Bd e,B可写成分块的形式： B=Bi Bj Bm,B称为应变矩阵，它的元素都只与单元的几何性质有关的常量。这种单元称为平面问题的常应变三角形单元。,单元应变与单元节点位移关系,(i, j, m),26,学习交流PPT,物理方程 s =De 而 e =Bde s = DBd e = S d e (求应力的表达式) 记 S=DB S应力矩阵: S=Si Sj

7、 Sm,2.3 由单元节点位移求单元的应力,27,学习交流PPT,1) 平面应力问题：,代入D及B 得: S =Si Sj Sm. 对于平面应力： (i, j, m ),28,学习交流PPT,将上式中以 E/(1-m2) 代 E，以 m/(1-m) 代 m 则子矩阵： (i, j, m),2)对于平面应变问题：,29,学习交流PPT,考虑单元平衡，节点力是作用在单元上的外力，与单元应力平衡。有限元法中以虚功方程代替平衡方程。 节点力列阵及单元内应力列阵： 单元节点力是指单元和节点相连接的内力； 考虑节点平衡，节点力为外力，与节点外载荷平衡；,2.4 由节点位移求单元节点力(求单元刚度矩阵),3

8、0,学习交流PPT,节点虚位移列阵及虚应变：,31,学习交流PPT,令实际受力状态在虚位移状态上做虚功，虚功方程： 由 e =Bd e 知 e *=Bd *e 则 e *T=(d *e)TBT 由于d *e中的元素为常量，提至前积分号前，故： (对于三角形三节点单元，B和s 为常量，单元厚度t也是常量；,为三角形单元面积，用表示）,32,学习交流PPT,单元刚度矩阵,物理方程s = De ,简记为 Fe = ked e,ke = BTDBt 或 = BTSt,几何方程e = Bd e,33,学习交流PPT,对于平面应力问题：,(r =i, j, m; s =i, j, m).,对于平面应变问题

9、：,34,学习交流PPT,求例4.2（p84）单元的单元刚度矩阵,(i, j, m),解：（1）求矩阵B,35,学习交流PPT,（2）求矩阵S,36,学习交流PPT,（3）求矩阵ke,ke=BTDBt =BTSt,37,学习交流PPT,代入a=1, b=2m 得： 可算出，当a=b时单元刚度矩阵与尺寸a,b无关,38,学习交流PPT,在单元节点力列阵Fe、单元应力列阵s e、单元应变列阵e e和单元节点位移列阵d e的四个列阵之间，存在五个转换关系，可得五个转换矩阵。,平衡条件,物理关系,几何关系,39,学习交流PPT,单元刚度矩阵的特性：,(1)单元刚度矩阵的物理意义： 单元刚度矩阵ke表示

10、了单元抵抗变形的能力，即表示了节点位移d e与节点力Fe之间的关系。 kij表示节点j发生单位位移时，其它节点位移分量均为零时，在节点i上产生的节点力。,40,学习交流PPT,(2)分块性质:单元刚度矩阵可以分块运算。,(4.51),按节点进行分块，则单元刚度矩阵的分块形式可写为：,(4.52),41,学习交流PPT,（3）对称性 单元刚度矩阵是一个对称矩阵 式(4.51)中： knl=kln (n=16; l=16) 分块形式中： krs=ksrT (r =i, j, m; s =i, j, m) （4）奇异性 单元刚度矩阵是一个奇异矩阵 |ke|=0，表明其逆矩阵不存在。即，如果给定了单元

11、节点位移可以得出唯一的节点力: Fe=ked e。反之，如果给出节点力却无法求出确定的节点位移。因为这时的单元未考虑所受约束时，可能存在不引起单元应力和节点力的刚性位移，这部分刚体位移由节点力是无法唯一确定的。 三节点三角形单元每行元素之和为零,42,学习交流PPT,例：证明图示单元刚度矩阵: kI=kIII,证明：由于单元刚度矩阵 ke=BTDBt 可知：当两个三角形单元几何尺寸相同时，t值和单元面积值均相同；当两个单元的材料性质相同时，弹性矩阵D也时相同的。故ke是否相同，取决于矩阵B是否相同。,43,学习交流PPT,不难验证,I、III单元的上述br, cr (i, j, m)值均相等。

12、 结论： 两个单元刚度矩阵ke相等的条件为：只要两单元的形状、大小，方向和单元弹性常数均相同，并且编号的方式也相同(如按逆时针方向编号为 i, j, m，直角顶点编号为m)，则两个单元的刚度矩阵时相等的。,44,学习交流PPT,2.5 单元载荷的移置. （离散时每个单元受载作用于节点上）,(a)原则：将单元载荷向节点处移置，按照虚功等效的原则进行。对于变形体(包括弹性体)，虚功等效是指原载荷与节点载荷载在任何虚位移上做的虚功相等。当位移模式确定后，载荷移置(或分解)其结果是唯一的。 虚功等效包含了刚体体系的静力等效，当虚位移为刚体位移时，虚功等效即为静力等效，静力等效是虚功等效的特例。,离散,

13、45,学习交流PPT,(b)载荷移置公式,(1)集中力 设单元i, j, m中任一点M(x, y)处受有集中力P=Px PyT，移置到该单元各节点处载荷列阵为 Re=Xi Yi Xj Yj Xm YmT,(1) 集中力P=Px PyT,46,学习交流PPT,(1)集中力,假设该单元发生一微小虚位移，M点相应的虚位移为 f *，该单元各节点处相应虚位移为d * ，由静力等效原理，载荷与节点等效载荷在虚位移上所作虚功相等： (d *e)TRe= f *TP 将 f *=N d *e 代入上式： 有 (d *e)TRe = f *TP=(d *e)TNTP 则 Re=NTP Re=Xi Yi Xj

14、Yj Xm YmT =NiPx NiPy NjPx NjPy NmPx NmPyT,47,学习交流PPT,设单元i, j, m的一边受有分布的面力P=Px Py将微元面积tds上的面力合力Ptds当作集中载荷，可得面力的移置公式:,(2)面力,(2) 面力,48,学习交流PPT,(3)体力,设单元i, j, m受有分布体力G=Px PyT将微分体积tdxdy上的体力合力Gtdxdy当作集中载荷dG同理可得：,(3) 分布体力G=X Y T,49,学习交流PPT,(4) 三节点三角形单元上同时有体力、面力和集中力等,(a) 集中力P=Px PyT,(b) 分布体力G=X Y T,(c) 面力,(

15、d) 单元虚位移,d * =u*i v*i u*j v*j u*m v*mT,50,学习交流PPT,单元各节点处载荷列阵为 Re=Xi Yi Xj Yj Xm YmT 应用虚功等效原则： 将 f *=Nd *e 代入上式： 虚位移是任意的，从而矩阵d *eT 也是任意的，故：,51,学习交流PPT,例：设三角形单元i, j, m的ij边作用有线形分布的法向载荷，i和j两点的压力集度分别为qi和qj，试用公式 求其等效节点载荷。单元厚度为t，节点坐标如图示。,52,学习交流PPT,解：计算常数 ai=xj ymxm yj=0； bi = yj ym= ym； ci = xmxj = xm； aj

16、=xm yi xi ym=xm yj； bj = ym yi=ym yi； cj=xi xm= xm； am=xi yj xj yi=0； bm=yi yj = yi； cm=xj xi= 0. 计算形函数,53,学习交流PPT,由 得： 在边界jm和mi上的面力为零，故上式积分中后两项为0，在ij边上的面力分量可表示为： 代入上式中得：,计算等效节点载荷,54,学习交流PPT,= 0,积分沿逆时针方向，有ds=-dy 所以,在ij边上 x = 0，代入 Ni Nj Nm中：,55,学习交流PPT,引入支承条件,解方程求位移,求单元应力,3 整体分析,建立整体刚度矩阵,任务：建立整个结构的总刚

17、度方程；引入边界条件解方程。,56,学习交流PPT,3.1 建 立 整 体 刚 度 矩 阵,3,P/2,P/2,4,1,2,y,i,i,j,j,m,m,I,II,x,57,学习交流PPT,单元节点受力分析,作用于节点的集中力,单元作用在节点上的力,作用在节点上等效载荷,58,学习交流PPT,由节点的平衡条件，有平衡方程： 环绕节点2的所有单元求和： F2=U2e V2eT 单元e作用在2节点上的节点力 R2=X2 Y2T 绕节点2的各单元作用在2节点上的等效节点载荷及直接作用于2的集中力之和。,59,学习交流PPT,每一个单元可用节点位移表示节点力，采用分块P89 4.46,考虑节点2，对应I

18、单元的局部编号为i即I单元i节点，故I单元i节点上的节点力： 代入 得,60,学习交流PPT,结构的整体刚度矩阵,结构的节点位移列阵,结构的节点载荷列阵,每个节点均可得到类似的方程，每个节点可写出两个平衡方程，按节点的序号排列，用矩阵表示：,61,学习交流PPT,结构的节点载荷列阵 R = X1 Y1 X2 Y2 T 故 其中 l=14 l, k=14,62,学习交流PPT,整体刚度矩阵的集成规则：,(1)先求出每个单元的刚度矩阵ke； (2)将ke的每个方块kije进行换号，换成对应的整体编号; (3)将换号后的子块填入整体刚度矩阵上对应的位置。 (4)若在同一位置上有几个单元的相应子块填入

19、同一位置，则进行叠加。,63,学习交流PPT,图4.24 刚度矩阵,64,学习交流PPT,i j m 2 4 1,i j m 4 2 3,i 2 j 4 m 1,i 4 j 2 m 3,65,学习交流PPT,解: (1)求各单元刚度矩阵(P90已求)P90 4-50 求得单元的刚度矩阵为： kI=kII (I II单元刚度矩阵相等),例：4.2。,66,学习交流PPT,(2)整体刚度矩阵k,换号的单元刚度矩阵：,67,学习交流PPT,(3)将、的单元刚度矩阵填入 (P101 式4.62),如果I II都有就相加便于计算机运算,68,学习交流PPT,总体刚度矩阵的特性：,(1) 对称性 总体刚度

20、矩阵是一个对称矩阵。因单元刚度矩阵升阶后对称性不变，由之合成的总体刚度矩阵自然是对称矩阵。 (2) 奇异性 总体刚度矩阵行列式的值 |K|=0,69,学习交流PPT,(3) 稀疏性 总体刚度矩阵是一个稀疏矩阵；即矩阵中的绝大多数元素为0，非0元素只占元素总数的很小的一部分。因为只有当节点ln相关时Kln才不是0，P101图4.26。 (4) 带状分布规律 分布在以主对角线为中心的带状区域内。,总体刚度矩阵的特性：,70,学习交流PPT,3.2 节点载荷列阵,集成法求整体结构的节点载荷列阵步骤： (1)求出每个单元的等效节点载荷； (2)单元等效节点子向量(k=i, j, m)换号，换成对应整体

21、编号； (3)换号后等效节点载荷子向量送到整体节点载荷列阵对应位置； (4)同一位置上若有多个单元的等效节点载荷子向量，叠加； (5)节点上有直接作用的节点载荷，按整体节点号进行叠加； (6)若节点k具有水平和垂直方向的支承，支承反力为未知量，可暂设为Rke=Qxk QykT。,71,学习交流PPT,图4.24中：,R=Qx1 Qy1 X2 Y2 X3 Y3 Qx4 Qy4 T =Qx1 Qy1 0 -P/2 0 -P/2 Qx4 Qy4 T,72,学习交流PPT,3.3 引入支承条件,(4.66),例：,73,学习交流PPT,基本未知量为节点2的位移，只需在方程中抽出第三、四行即可： 矩阵形

22、式 (4.67),74,学习交流PPT,为了便于编程，修改后的矩阵仍保持原矩阵K的阶数、排列次序及矩阵对称性。式(4.67)扩大成：,等效：u1=v1=u3=v3= 0 与支承条件及前述矩阵一致,75,学习交流PPT,多个节点：节点n的水平位移un=0，则Kd =R改为： K中的2n-1行与列中主对角元素为1，其他为0 载荷向量R中2n-1元素置0 若vn=0 则对应K中2n行和列作上述修改。,76,学习交流PPT,仍以课本P98 图4.24 u1= v1= u4= v4= 0,77,学习交流PPT,3.4 解方程 求节点位移 Kd =R,将整体刚度矩阵K代入支承条件 图4.24,78,学习交

23、流PPT,设m =1/3得：,得：,于是节点的位移向量为 d =u1 v1 u2 v2 u3 v3 u4 v4T =P/Et0 0 -1.50 -8.42 1.88 -8.99 0 0T,79,学习交流PPT,单元1：当a=1, b=2, m=1/3时，由例4.3,3.5 求单元应力s e=Sed e,80,学习交流PPT,同理可求单元,81,学习交流PPT,3.6求节点力及支承反力 Fe=ked e,单元的节点力 a=1,b=2, m=1/3,82,学习交流PPT,单元的节点力,同理可得FII 也进行验算; 由受力图可得,验算,83,学习交流PPT,总结：,有限元求解弹性力学平面问题步骤如下

24、： (1)整理原始数据，结构离散化，对单元和节点编号。 (2)求单元刚刚度矩阵ke (3)用刚度集成法，形成结构整体刚度矩阵K (4)求节点等效载荷，写出载荷列阵 R (5)引入支承条件 (6)解K d =R 求出节点位移d (7)求单元应力s e=S ed e (8)求单元节点力 Fe=Ked e (9)整理结果，作节点位移图及应力图。,84,学习交流PPT,一、matlab 基础,在 matlab 的提示符 符“”下输入命令. 3*4+5 ans= 17 cos(30*pi/180) ans= 0.8660 x=4 x=4 2/sqrt(3+x) ans= 0.7559,4 线性三角形弹性

25、力学平面问题的matlab程序：,(德)卡坦-韩来彬 清华大学出版社,85,学习交流PPT,若不让 matlab 输出运算数据，在命令行的结尾输入分号：, y=32; z=5; x=2*y-z; w=3*y+4*z w=116,matlab区分大小写,x=1 x=1 X=2 X=2 x x=1,使用help命令可以获得所有matlab命令的详细用法 help inv,86,学习交流PPT,下面的例子显示如何输入矩阵并实现简单的矩阵运算。 x=1 2 3;4 5 6;7 8 ,1 2 3 4 5 6 7 8 9,x=, y=2 ; 0 ; -3,y = 2 0 -3, w=x*y w= -7 -

26、10 -13,87,学习交流PPT,求解下面的联立方程组 采用高斯消去法解方程组: A= 2 -1 3 0 ;1 5 -2 4;2 0 3 -2; 1 2 3 4, b= 3; 1; -2 ; 2 b= 3 1 -2 2, x=Ab. x= 1.9259 -1.8148 -0.8889 1.5926,A= 2 -1 3 0 1 5 -2 4 2 0 3 -2 1 2 3 4,88,学习交流PPT,采用如下方法也可求 x=inv(A)*b x = 1.9259 -1.8148 -8.8889 1.5926 时间长，矩阵大尤其长,D=1 2 3 4 5 ; 2 4 6 8 9 ; 2 4 6 2

27、4 ; 1 7 2 3 -2 ; 9 0 2 3 1 ,D= 1 2 3 4 5 2 4 6 8 9 2 4 6 2 4 1 7 2 3 -2 9 0 2 3 1,89,学习交流PPT,可以以矩阵中提取 24 行，35列作子矩时：, E=D(2:4, 3:5) E= 6 8 9 6 2 4 2 3 -2,提取D的第3列： F=D (1:5, 3),F = 3 6 6 2 2,提取D的第2行作为子矩阵 G=D (2, 1:5) G = 2 4 6 8 9,90,学习交流PPT,提取D中的4行 3列 H=D (4, 3) H=2 绘制y=f(x) 定义 x. y. 再用plot(x,y) 绘图。

28、x= 1 2 3 4 5 6 7 8 9 10 y=x.2 plot(x,y),91,学习交流PPT,二、 线性三角元 1.基本方程. Linear triangular element 线性形函数.,系数： 弹性模量E 泊松比m 厚度t,单元刚度阵：,k= BTDBt,2= xi(yj-ym)+xj(ym-yi)+xm(yi-yj),92,学习交流PPT,ai bi ci (P79 4.26).,平面应力,平面应变,显然线性三角形元有6个自由度,每个节点2个自由度。,93,学习交流PPT,对一个有n个节点的结构，整体K为2n*2n Kd =F 边界条件手动赋值 采用高斯消去求解 s=Sd 求

29、得单元应力矩,2，用到的Matlab函数 线性三角形元用到的5个matlab函数分别为： (1)LinearTriangleElementArear(xi，yi，xj，yj，xm , ym) 该函数根据给出的i, j, m的节点坐标返回单元的面积。,94,学习交流PPT,(2)LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p) 该函数用于计算弹性模量为E，泊松比为m，厚度为t，以及i, j, m节点坐标已知的线性三角形单元的刚度矩阵。p=1表明函数用于平面应力情况。P=2表明用于平面应变情况。返回6*6的单元刚度矩阵k。,(3)

30、 LinearTriangleElementAssemble(K, k,i, j, m) 将连接i ,j ,m线性三角形元的单元刚度矩阵k集成到整体刚度矩阵K。返回2n*2n的整体刚度矩阵 K。,(4) LinearTriangleElementStresses(E,NU,t,xi,yi,xj,yjxm,ym,p,u) u为单位位移矢量, 返回单元应力,(5) LinearTriangleElementPStresses (sigma) 计算单元主应力，返回3*1. sigma1,sigma2,thetaT,95,学习交流PPT,例：求解如图所示的受均布载荷作用的薄平板结构， 将平板离散比化为两个线性三角形元，如右图示:,E=210MPa m=0.3 t=