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

.,1,4.3平面问题的有限单元法,三角形三节点平面单元,.,2,结构离散化,单元分析,整体分析,有限元分析的基本步骤：,.,3,1结构离散化,例：图示一悬臂梁，梁的厚度为t，设泊松比m=1/3，弹性模量为E，试用三节点三角形单元进行离散。,.,4,2单元分析,单元分析的主要内容：由节点位移求内部任一点的位移，由节点位移求单元应变、应力和节点力。单元分析的步骤：,节点位移,单元内部各点位移,单元应变,单元应力,节点力,(1),(2),(3),(4),单元分析,.,5,(1)单元位移模式有限元法中，通常采用位移法进行计算，即取节点的位移分量为基本未知量，单元中的位移、应变、应力等物理量，都和基本未知量相关联。节点i的位移分量可写成di=uiviT单元节点位移向量de可写成de=didjdmT=uiviujviumvmT,三角形三节点单元,2.1由节点位移求单元内部任一点位移,.,6,六个位移分量需六个待定参数a1，a2，a6设单元内任一点位移分量是坐标(x,y)的线性函数：(4.22)写成矩阵的形式为：(4.23),.,7,设节点i，j，m坐标分别是xi，yi；xj，yj；xm，ym。把三个节点的坐标及其水平位移代入式(4.22)中得：(a)解得：(b)其中，,(2)由单元节点位移de求位移参数a,.,8,对v同理可列出a4、a5、a6的方程：vi=a4+a5xi+a6yii,j,m解得：a4、a5、a6为了书写的简便与规格化，引入记号ai,bi,ci分别为行列式A中各行各个元素的代数余子式,.,9,解出a1a6结果：,a=Lde,为三角形单元面积。,.,10,将a写成矩阵形式，有a=Lde,由单元节点位移de求单元内部任一点位移d(x,y)d(x,y)=m(x,y)a=m(x,y)Lde,.,11,.,12,Ni,Nj,Nm是坐标的连续函数，反映单元内位移的分布状态，称为位移的形状函数，简称形函数(shapefunction)。矩阵N称为形函数矩阵(shapefunctionmatrix)。形函数物理意义：Ni(x,y)，节点i单位位移，其它节点位移分量为0，单元内部产生位移分布形状,.,13,在单元任一点上，三个形函数之和等于1。形函数Ni在i点的函数值为1，在j点及m点的函数值为零。三角形单元i,j,m在ij边上的形函数与第三个顶点的坐标无关。,形函数的性质：,.,14,例：求图示单元和单元的形函数矩阵,.,15,（a）,(b),分别如上图所示：,.,16,单元如图所示。设a=1m,b=2m.,(或直接由图形可知其面积),.,17,求系数ai,aj,am,bi,bj,bm,ci,cj,cm,.,18,求形函数矩阵,代入相关常数：,.,19,将a=1,b=2代入得：,.,20,求常数,单元如图所示=ab/2,.,21,求形函数矩阵,.,22,将a=1,b=2m代入上式得：,.,23,作业：,已知三角形三节点单元坐标如图示，设单元中一点A的坐标(0.5,0.2)，已知三角形三节点单元i节点位移(2.0,1.0)，j节点位移(2.1,1.1)，m节点位移(2.15,1.05)，1)写出单元的位移函数；2)求A点的位移分量。,A,.,24,2.2由节点位移求单元的应变,几何方程,.,25,.,26,简记为e=Bde,B可写成分块的形式：B=BiBjBm,B称为应变矩阵，它的元素都只与单元的几何性质有关的常量。这种单元称为平面问题的常应变三角形单元。,单元应变与单元节点位移关系,(i,j,m),.,27,物理方程s=De而e=Bdes=DBde=Sde(求应力的表达式)记S=DBS应力矩阵:S=SiSjSm,2.3由单元节点位移求单元的应力,.,28,1)平面应力问题：,代入D及B得:S=SiSjSm.对于平面应力：(i,j,m),.,29,将上式中以E/(1-m2)代E，以m/(1-m)代m则子矩阵：(i,j,m),2)对于平面应变问题：,.,30,考虑单元平衡，节点力是作用在单元上的外力，与单元应力平衡。有限元法中以虚功方程代替平衡方程。节点力列阵及单元内应力列阵：单元节点力是指单元和节点相连接的内力；考虑节点平衡，节点力为外力，与节点外载荷平衡；,2.4由节点位移求单元节点力(求单元刚度矩阵),.,31,节点虚位移列阵及虚应变：,.,32,令实际受力状态在虚位移状态上做虚功，虚功方程：由e=Bde知e*=Bd*e则e*T=(d*e)TBT由于d*e中的元素为常量，提至前积分号前，故：(对于三角形三节点单元，B和s为常量，单元厚度t也是常量；,为三角形单元面积，用表示）,.,33,单元刚度矩阵,物理方程s=De,简记为Fe=kede,ke=BTDBt或=BTSt,几何方程e=Bde,.,34,对于平面应力问题：,(r=i,j,m;s=i,j,m).,对于平面应变问题：,.,35,求例4.2（p84）单元的单元刚度矩阵,(i,j,m),解：（1）求矩阵B,.,36,（2）求矩阵S,.,37,（3）求矩阵ke,ke=BTDBt=BTSt,.,38,.,39,代入a=1,b=2m得：可算出，当a=b时单元刚度矩阵与尺寸a,b无关,.,40,在单元节点力列阵Fe、单元应力列阵se、单元应变列阵ee和单元节点位移列阵de的四个列阵之间，存在五个转换关系，可得五个转换矩阵。,平衡条件,物理关系,几何关系,.,41,单元刚度矩阵的特性：,(1)单元刚度矩阵的物理意义：单元刚度矩阵ke表示了单元抵抗变形的能力，即表示了节点位移de与节点力Fe之间的关系。kij表示节点j发生单位位移时，其它节点位移分量均为零时，在节点i上产生的节点力。,.,42,(2)分块性质:单元刚度矩阵可以分块运算。,(4.51),按节点进行分块，则单元刚度矩阵的分块形式可写为：,(4.52),.,43,（3）对称性单元刚度矩阵是一个对称矩阵式(4.51)中：knl=kln(n=16;l=16)分块形式中：krs=ksrT(r=i,j,m;s=i,j,m)（4）奇异性单元刚度矩阵是一个奇异矩阵|ke|=0，表明其逆矩阵不存在。即，如果给定了单元节点位移可以得出唯一的节点力:Fe=kede。反之，如果给出节点力却无法求出确定的节点位移。因为这时的单元未考虑所受约束时，可能存在不引起单元应力和节点力的刚性位移，这部分刚体位移由节点力是无法唯一确定的。三节点三角形单元每行元素之和为零,.,44,例：证明图示单元刚度矩阵:kI=kIII,证明：由于单元刚度矩阵ke=BTDBt可知：当两个三角形单元几何尺寸相同时，t值和单元面积值均相同；当两个单元的材料性质相同时，弹性矩阵D也时相同的。故ke是否相同，取决于矩阵B是否相同。,.,45,不难验证,I、III单元的上述br,cr(i,j,m)值均相等。结论：两个单元刚度矩阵ke相等的条件为：只要两单元的形状、大小，方向和单元弹性常数均相同，并且编号的方式也相同(如按逆时针方向编号为i,j,m，直角顶点编号为m)，则两个单元的刚度矩阵时相等的。,.,46,2.5单元载荷的移置.（离散时每个单元受载作用于节点上）,(a)原则：将单元载荷向节点处移置，按照虚功等效的原则进行。对于变形体(包括弹性体)，虚功等效是指原载荷与节点载荷载在任何虚位移上做的虚功相等。当位移模式确定后，载荷移置(或分解)其结果是唯一的。虚功等效包含了刚体体系的静力等效，当虚位移为刚体位移时，虚功等效即为静力等效，静力等效是虚功等效的特例。,离散,.,47,(b)载荷移置公式,(1)集中力设单元i,j,m中任一点M(x,y)处受有集中力P=PxPyT，移置到该单元各节点处载荷列阵为Re=XiYiXjYjXmYmT,(1)集中力P=PxPyT,.,48,(1)集中力,假设该单元发生一微小虚位移，M点相应的虚位移为f*，该单元各节点处相应虚位移为d*，由静力等效原理，载荷与节点等效载荷在虚位移上所作虚功相等：(d*e)TRe=f*TP将f*=Nd*e代入上式：有(d*e)TRe=f*TP=(d*e)TNTP则Re=NTPRe=XiYiXjYjXmYmT=NiPxNiPyNjPxNjPyNmPxNmPyT,.,49,设单元i,j,m的一边受有分布的面力P=PxPy将微元面积tds上的面力合力Ptds当作集中载荷，可得面力的移置公式:,(2)面力,(2)面力,.,50,(3)体力,设单元i,j,m受有分布体力G=PxPyT将微分体积tdxdy上的体力合力Gtdxdy当作集中载荷dG同理可得：,(3)分布体力G=XYT,.,51,(4)三节点三角形单元上同时有体力、面力和集中力等,(a)集中力P=PxPyT,(b)分布体力G=XYT,(c)面力,(d)单元虚位移,d*=u*iv*iu*jv*ju*mv*mT,.,52,单元各节点处载荷列阵为Re=XiYiXjYjXmYmT应用虚功等效原则：将f*=Nd*e代入上式：虚位移是任意的，从而矩阵d*eT也是任意的，故：,.,53,例：设三角形单元i,j,m的ij边作用有线形分布的法向载荷，i和j两点的压力集度分别为qi和qj，试用公式求其等效节点载荷。单元厚度为t，节点坐标如图示。,.,54,解：计算常数ai=xjymxmyj=0；bi=yjym=ym；ci=xmxj=xm；aj=xmyixiym=xmyj；bj=ymyi=ymyi；cj=xixm=xm；am=xiyjxjyi=0；bm=yiyj=yi；cm=xjxi=0.计算形函数,.,55,由得：在边界jm和mi上的面力为零，故上式积分中后两项为0，在ij边上的面力分量可表示为：代入上式中得：,计算等效节点载荷,.,56,=0,积分沿逆时针方向，有ds=-dy所以,在ij边上x=0，代入NiNjNm中：,.,57,引入支承条件,解方程求位移,求单元应力,3整体分析,建立整体刚度矩阵,任务：建立整个结构的总刚度方程；引入边界条件解方程。,.,58,3.1建立整体刚度矩阵,3,P/2,P/2,4,1,2,y,i,i,j,j,m,m,I,II,x,.,59,单元节点受力分析,作用于节点的集中力,单元作用在节点上的力,作用在节点上等效载荷,.,60,由节点的平衡条件，有平衡方程：环绕节点2的所有单元求和：F2=U2eV2eT单元e作用在2节点上的节点力R2=X2Y2T绕节点2的各单元作用在2节点上的等效节点载荷及直接作用于2的集中力之和。,.,61,每一个单元可用节点位移表示节点力，采用分块P894.46,考虑节点2，对应I单元的局部编号为i即I单元i节点，故I单元i节点上的节点力：代入得,.,62,结构的整体刚度矩阵,结构的节点位移列阵,结构的节点载荷列阵,每个节点均可得到类似的方程，每个节点可写出两个平衡方程，按节点的序号排列，用矩阵表示：,.,63,结构的节点载荷列阵R=X1Y1X2Y2T故其中l=14l,k=14,.,64,整体刚度矩阵的集成规则：,(1)先求出每个单元的刚度矩阵ke；(2)将ke的每个方块kije进行换号，换成对应的整体编号;(3)将换号后的子块填入整体刚度矩阵上对应的位置。(4)若在同一位置上有几个单元的相应子块填入同一位置，则进行叠加。,.,65,图4.24刚度矩阵,.,66,ijm241,ijm423,i2j4m1,i4j2m3,.,67,解:(1)求各单元刚度矩阵(P90已求)P904-50求得单元的刚度矩阵为：kI=kII(III单元刚度矩阵相等),例：4.2。,.,68,(2)整体刚度矩阵k,换号的单元刚度矩阵：,.,69,(3)将、的单元刚度矩阵填入(P101式4.62),如果III都有就相加便于计算机运算,.,70,总体刚度矩阵的特性：,(1)对称性总体刚度矩阵是一个对称矩阵。因单元刚度矩阵升阶后对称性不变，由之合成的总体刚度矩阵自然是对称矩阵。(2)奇异性总体刚度矩阵行列式的值|K|=0,.,71,(3)稀疏性总体刚度矩阵是一个稀疏矩阵；即矩阵中的绝大多数元素为0，非0元素只占元素总数的很小的一部分。因为只有当节点ln相关时Kln才不是0，P101图4.26。(4)带状分布规律分布在以主对角线为中心的带状区域内。,总体刚度矩阵的特性：,.,72,3.2节点载荷列阵,集成法求整体结构的节点载荷列阵步骤：(1)求出每个单元的等效节点载荷；(2)单元等效节点子向量(k=i,j,m)换号，换成对应整体编号；(3)换号后等效节点载荷子向量送到整体节点载荷列阵对应位置；(4)同一位置上若有多个单元的等效节点载荷子向量，叠加；(5)节点上有直接作用的节点载荷，按整体节点号进行叠加；(6)若节点k具有水平和垂直方向的支承，支承反力为未知量，可暂设为Rke=QxkQykT。,.,73,图4.24中：,R=Qx1Qy1X2Y2X3Y3Qx4Qy4T=Qx1Qy10-P/20-P/2Qx4Qy4T,.,74,3.3引入支承条件,(4.66),例：,.,75,基本未知量为节点2的位移，只需在方程中抽出第三、四行即可：矩阵形式(4.67),.,76,为了便于编程，修改后的矩阵仍保持原矩阵K的阶数、排列次序及矩阵对称性。式(4.67)扩大成：,等效：u1=v1=u3=v3=0与支承条件及前述矩阵一致,.,77,多个节点：节点n的水平位移un=0，则Kd=R改为：K中的2n-1行与列中主对角元素为1，其他为0载荷向量R中2n-1元素置0若vn=0则对应K中2n行和列作上述修改。,.,78,仍以课本P98图4.24u1=v1=u4=v4=0,.,79,3.4解方程求节点位移Kd=R,将整体刚度矩阵K代入支承条件图4.24,.,80,设m=1/3得：,得：,于是节点的位移向量为d=u1v1u2v2u3v3u4v4T=P/Et00-1.50-8.421.88-8.9900T,.,81,单元1：当a=1,b=2,m=1/3时，由例4.3,3.5求单元应力se=Sede,.,82,同理可求单元,.,83,3.6求节点力及支承反力Fe=kede,单元的节点力a=1,b=2,m=1/3,.,84,单元的节点力,同理可得FII也进行验算;由受力图可得,验算,.,85,总结：,有限元求解弹性力学平面问题步骤如下：(1)整理原始数据，结构离散化，对单元和节点编号。(2)求单元刚刚度矩阵ke(3)用刚度集成法，形成结构整体刚度矩阵K(4)求节点等效载荷，写出载荷列阵R(5)引入支承条件(6)解Kd=R求出节点位移d(7)求单元应力se=Sede(8)求单元节点力Fe=Kede(9)整理结果，作节点位移图及应力图。,.,86,一、matlab基础,在matlab的提示符符“”下输入命令.3*4+5ans=17cos(30*pi/180)ans=0.8660x=4x=42/sqrt(3+x)ans=0.7559,4线性三角形弹性力学平面问题的matlab程序：,(德)卡坦-韩来彬清华大学出版社,.,87,若不让matlab输出运算数据，在命令行的结尾输入分号：,y=32;z=5;x=2*y-z;w=3*y+4*zw=116,matlab区分大小写,x=1x=1X=2X=2xx=1,使用help命令可以获得所有matlab命令的详细用法helpinv,.,88,下面的例子显示如何输入矩阵并实现简单的矩阵运算。x=123;456;78,123456789,x=,y=2;0;-3,y=20-3,w=x*yw=-7-10-13,.,89,求解下面的联立方程组采用高斯消去法解方程组:A=2-130;15-24;203-2;1234,b=3;1;-2;2b=31-22,x=Ab.x=1.9259-1.8148-0.88891.5926,A=2-13015-24203-21234,.,90,采用如下方法也可求x=inv(A)*bx=1.9259-1.8148-8.88891.5926时间长，矩阵大尤其长,D=12345;24689;24624;1723-2;90231,D=1234524689246241723-290231,.,91,可以以矩阵中提取24行，35列作子矩时：,E=D(2:4,3:5)E=68962423-2,提取D的第3列：F=D(1:5,3),F=36622,提取D的第2行作为子矩阵G=D(2,1:5)G=24689,.,92,提取D中的4行3列H=D(4,3)H=2绘制y=f(x)定义x.y.再用plot(x,y)绘图。x=12345678910y=x.2plot(x,y),.,93,二、线性三角元1.基本方程.Lineartriangularelement线性形函数.,系数：弹性模量E泊松比m厚度t,单元刚度阵：,k=BTDBt,2=xi(yj-ym)+xj(ym-yi)+xm(yi-yj),.,94,aibici(P794.26).,平面应力,平面应变,显然线性三角形元有6个自由度,每个节点2个自由度。,.,95,对一个有n个节点的结构，整体K为2n*2nKd=F边界条件手动赋值采用高斯消去求解s=Sd求得单元应力矩,2，用到的Matlab函数线性三角形元用到的5个matlab函数分别为：(1)LinearTriangleElementArear(xi，yi，xj，yj，xm,ym)该函数根据给出的i,j,m的节点坐标返回单元的面积。,.,96,(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)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,.,97,例：求解如图所示的受均布载荷作用的薄平板结构，将平板离散比化为两个线性三角形元，如右图示:,E=210MPam=0.3