03杆梁结构.ppt

第三章杆梁结构有限元,3.1有限元分析的基本力学原理3.2局部坐标系中的杆单元描述3.3FEA的基本步骤3.4杆单元坐标变换3.5杆单元分析的MATLAB程序3.6梁单元及坐标变换3.7梁单元分析的MATLAB程序3.8单元形状函数矩阵与刚度矩阵的性质,在杆件系统中根据单元受力的特点分成两大类：杆和梁。两端铰接，只受轴向力的基本结构称为杆单元。受轴向力和弯矩、扭矩、剪力共同作用的基本结构称为梁单元。根据其截面的一个主轴与它所承受的载荷是否在同一平面内，杆件单元可分成平面单元或空间单元。,3.1有限元分析的基本力学原理,1.弹塑性力学基本假设,小变形假设：变形后物体内各点的位移远小于原尺寸，可忽略变形引起的几何变化。无初应力假设：外力作用前各点应力为零。载荷分类：体积力和表面力,材料构造模型材料力学性质模型结构计算模型,材料构造模型,连续性假设：固体材料是连续介质应力、应变和位移都可用连续函数来描述。各向同性及均匀性假设：物体由同一类型的均匀性材料组成物体内任意一点在各方向上的物理性质响应相同（各向同性）；物体每一部分具有相同的性质，物理常数不随位置的变化而变化（均匀性）。,力学分析模型的三个层次：,结构计算模型,材料力学性质模型,弹性材料塑性材料黏性材料,材料构造模型和结构计算模型是弹塑性力学问题的共同基础，材料力学性质模型的选取需根据材料本身的力学性质、工作环境及限定的研究范围确定。,2.基本变量,位移u(ui),应变（ij）,应力（ij）,物体变形后的形状,物体变形的程度,物体的受力状态,3.弹性力学基本方程,（1）平衡方程（受力情况描述）,i，j自由下标，i，jx，y，z,（2）几何方程（变形程度描述）,指标形式：,（3）本构关系（材料描述）,（4）协调方程（Saint-Venant方程）,D弹性矩阵,（5）边界条件,应力边界条件,位移边界条件,4.能量概念,功能定理,在等温或绝热过程中，物体所获得的应变能等于外力所做的功。,5.虚位移原理,弹性体处于平衡状态的必要与充分条件：对于任意的、满足相容条件的虚位移，外力所做的功等于弹性体所接受的总虚变形功（外力在虚位移上的虚功等于应力在虚应变上的虚功）。,总虚变形功：,对于平面问题：,虚位移原理,总外力虚功：,虚位移：满足许可位移条件的，任何微小的假想位移,6.最小势能原理,总势能：,形变势能：,外力势能：,形变势能变分:,外力势能变分:,即：形变势能的变分表达式与虚变形功的表达式完全相同。,即：外力势能的变分表达式与外力虚功负值的表达式完全相同。,在几何可能的一切容许位移和形变中，真正的位移和形变使总势能取最小值；反之，使总势能取最小值者也必是真正的位移和形变。,最小势能原理,有限元法的理论基础是变分原理。最常用的变分原理有最小势能原理、最小余能原理和混合变分原理。采用不同的变分原理，将得到不同的未知场变量：采用最小势能原理，需假设单元内位移场函数的形式位移法；采用最小余能原理，需假设单元内应力场函数的形式应力法；采用混合变分原理，需假设单元内某些位移和某些应力场函数的形式混合法。,6.杆件分析的基本力学原理应用,基本变量,沿x方向移动为位移ux(x),沿x方向的相对伸长(或缩短)量为应变x(x),沿x方向的单位横截面上的受力为应力x(x),基本方程,平衡方程（体力为0）,几何方程,物理方程,边界条件,求解方法:直接求解方法：简单问题可以由3个方程来直接求解3个变量。基于试函数的间接方法：可以先选取一个变量(如位移)作为最基本的待求变量，将其它变量都用它来表达，并采用间接的近似求解方法。具体的做法为：先对待求的位移变量假设一种事先满足位移边界条件的可能解(其中有一些待定的系数)，称为试函数(trailfunction)，基于虚功原理或最小势能原理求出试函数中的那些待定系数。,位移边界条件,力边界条件,（1）直接求解法,（2）虚功原理求解,设有满足位移边界条件的位移场,u(x)=cx,虚应变能,外力虚功,由虚功原理,（3）最小势能原理求解,设有满足位移边界条件的位移场,u(x)=cx,该系统的势能,基于试函数的方法，包括虚功原理以及最小势能原理，仅计算系统的能量，实际上就是计算积分，然后转化为求解线性方程，不需求解微分方程，这样就大大地降低了求解难度。试函数的方法的关键就是如何构造出适合于所求问题的位移试函数，并且该构造方法还应具有规范性以及标准化，基于“单元”的构造方法就可以完全满足这些要求。,3.2局部坐标系中的杆单元描述,单元的描述包括单元的几何及节点描述、位移场、应变场、应力场、势能的描述；充分利用描述问题的三大类变量以及三大类方程来计算单元的势能，然后，由最小势能原理(或虚功原理)来得到单元的方程。单元内位移场的描述就是它的试函数的选取。,(1)单元的几何及节点描述,基本变量为节点位移(向量)列阵,节点力(向量)列阵,设单元的位移场为u(x),(2)单元位移场的表达,由Taylor级数,取前两项来作为该单元的位移插值模式,单元节点条件,形状函数矩阵(shapefunctionmatrix),节点位移列阵(nodaldisplacementvector),(3)单元应变,几何矩阵,(4)单元应力,应力矩阵,单元位移场,局部坐标系下单元刚度矩阵,单元节点外载,KEY！,(5)单元势能,(6)单元的刚度方程,E(1)=E(2)=2107Pa，A(1)2A(2)2cm2，l(1)l(2)10cm,节点位移向量,刚度矩阵,节点外载,单元1,节点位移向量,刚度矩阵,节点外载,单元2,集成以得到结构的总体势能,边界条件处理,u10,求解：由最小势能定理,各单元应变,各单元应力,节点1的反力,讨论,u10,先处理边界条件再应用最小势能原理先应用最小势能原理再处理边界条件,3.3FEA的基本步骤,1.物体几何的离散,具有相同特征的个体单元,2.单元分析：所有力学信息表示为节点位移,(1)单元节点位移,(2)单元位移模式（唯一确定性原则，完备性原则）,(3)力学量的表达（节点位移的函数）,形函数,形函数矩阵,单元的平衡关系,3.装配集成,整体的平衡关系,4.边界条件处理并求解：得到满足边界条件的许可位移场,u未知节点位移k已知节点位移,Pu未知节点载荷Pk已知节点载荷,5.其他力学量：,1.平面问题杆单元坐标变换,局部坐标系下单元节点位移,总体坐标系下单元节点位移,相互关系,3.4杆单元坐标变换,坐标变换矩阵,总体坐标系下单元刚度矩阵,总体坐标系下单元节点载荷,总体坐标系下单元刚度方程,2.空间问题杆单元坐标变换,局部坐标系下单元节点位移,总体坐标系下单元节点位移,杆单元轴线在总体坐标系中的方向余弦,坐标变换矩阵,总体坐标系下单元刚度矩阵,总体坐标系下单元节点载荷,总体坐标系下单元刚度方程,3.5杆单元分析的MATLAB程序,学习有限元方法的一个最佳途径，就是在充分掌握基本概念的基础上亲自编写有限元分析程序;MATLAB软件功能强大、编程逻辑直观、使用方便。,2D杆单元的有限元分析程序(Bar2D2Node),MATLAB中的反斜线符号“”就是采用高斯消去法,2D杆单元的有限元分析程序(Bar2D2Node),Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)该函数计算单元的刚度矩阵，输入弹性模量E，横截面积A，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2）和角度alpha（单位是度），输出单元刚度矩阵k(44)。2.Bar2D2Node_Assembly(KK,k,i,j)该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j，输出整体刚度矩阵KK。3.Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)该函数计算单元的应力，输入弹性模量E，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）和单位节点位移矢量u，返回单元应力标量。4.Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)该函数计算单元的应力，输入弹性模量E，横截面积A，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）和单元节点位移矢量u，返回单元节点力。,functionk=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)L=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);x=alpha*pi/180;C=cos(x);S=sin(x);k=E*A/L*C*CC*S-C*C-C*S;C*SS*S-C*S-S*S;-C*C-C*SC*CC*S;-C*S-S*SC*SS*S;,Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha),输入参数,2.Bar2D2Node_Assembly(KK,k,i,j),functionz=Bar2D2Node_Assembly(KK,k,i,j)DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;forn1=1:4forn2=1:4KK(DOF(n1),DOF(n2)=KK(DOF(n1),DOF(n2)+k(n1,n2);endendz=KK;,单元刚度矩阵和全体刚度矩阵的关系,3.Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u),functionstress=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)L=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);x=alpha*pi/180;C=cos(x);S=sin(x);stress=E/L*-C-SCS*u;,4.Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u),functionforces=Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)L=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);x=alpha*pi/180;C=cos(x);S=sin(x);forces=E*A/L*-C-SCS*u;,四杆桁架结构的有限元分析,各杆的弹性模量和横截面积都为2.95104/mm2,A=100mm2，试求解该结构的节点位移、单元应力以及支反力。,（1）结构的离散化与编号,（2）计算各单元的刚度矩阵(基于国际标准单位),建立一个工作目录，将所编制的用于平面桁架单元分析的四个MATLAB函数放置于该工作目录中，分别以各自函数的名称给出文件名，即：Bar2D2Node_Stiffness，Bar2D2Node_Assembly，Bar2D2Node_Stress，Bar2D2Node_Forces。,然后启动MATLAB，将工作目录设置到已建立的目录中，在MATLAB环境中，输入弹性模量E、横截面积A，各点坐标x1,y1,x2,y2,x3,y3,x4,y4,角度alpha1,alpha2和alpha3,分别针对单元1，2，3和4，调用四次Bar2D2Node_Stiffness，就可以得到单元的刚度矩阵。,（3）建立整体刚度方程,由于该结构共有4个节点，因此，设置结构总的刚度矩阵为KK(88)，先对KK清零，然后四次调用函数Bar2D2Node_Assembly进行刚度矩阵的组装。,（4）边界条件的处理及刚度方程求解,（5）支反力的计算,在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。将整体的位移列阵q（采用国际标准单位）代回原整体刚度方程，计算出所有的节点力P，即可找到对应的支反力。,（6）各单元的应力计算,先从整体位移列阵q中提取出单元的位移列阵，然后，调用计算单元应力的函数Bar2D2Node_Stress，就可以得到各个单元的应力分量。当然也可以调用上面的Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)函数来计算单元的集中力，然后除以面积求得单元应力。,3.6梁单元及坐标变换,1.局部坐标系中的纯弯梁单元,节点载荷向量,节点位移向量,单元位移场,（1）单元位移场,形函数矩阵,（2）单元应变场,基于中性层的坐标,单元几何函数矩阵,（3）单元应力场,单元应力函数矩阵,（4）单元势能,局部坐标系下单元刚度矩阵,局部坐标系下单元刚度矩阵,外力功,（5）单元刚度方程,Ke,2.局部坐标系中的平面梁单元,在纯弯单元上直接叠加轴向位移,节点载荷向量,节点位移向量,单元刚度方程,qe,Ke,杆单元刚度矩阵,纯弯梁单元刚度矩阵,3.平面问题梁单元的坐标变换,局部坐标系下单元节点位移,总体坐标系下单元节点位移,相互关系,坐标变换矩阵,qe,qe,总体坐标系下单元刚度矩阵,总体坐标系下单元节点载荷,总体坐标系下单元刚度方程,4.空间梁单元及其坐标变换,节点载荷向量,节点位移向量,单元刚度方程,qe,z2,单元刚度矩阵,结合杆单元和平面梁单元的刚度矩阵,空间梁单元的刚度矩阵,（1）对应于u1、u2的轴向位移,（2）对应于x1、x2的扭转角位移,（3）对应于xoy平面内位移v1、z1、v2、z2,（4）对应于zoy平面内位移w1、y1、w2、y2,坐标变换,局部坐标系下单元节点位移,总体坐标系下单元节点位移,相互关系（对节点1）,qe,qe,坐标变换矩阵,总体坐标系下单元刚度矩阵,总体坐标系下单元节点载荷,总体坐标系下单元刚度方程,3.7梁单元的常用等效节点载荷,非节点载荷作用下的梁单元须转化为节点载荷作用下的梁单元，即需得到非节点载荷的节点载荷等效值。,等效原则,常用梁,外力在单元位移场所做的功=等效节点力在节点位移上所做的功,3.8梁单元分析的MATLAB程序,平面问题梁单元,Beam2D2Node_Stiffness(E,I,A,L)计算单元的刚度矩阵，输入弹性模量E，横截面积A，惯性矩I，长度L，输出单元刚度矩阵k(66)。2.Beam2D2Node_Assemble(KK,k,i,j)进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j，输出整体刚度矩阵KK。3.Beam2D2Node_Forces(k,u)计算节点的力，输入单元刚度矩阵k，节点位移u，输出单元节点力forces。,u2,functionk=Beam2D2Node_Stiffness(E,I,A,L)k=E*A/L,0,0,-E*A/L,0,0;0,12*E*I/(L3),6*E*I/(L2),0,-12*E*I/(L3),6*E*I/(L2);0,6*E*I/(L2),4*E*I/L,0,-6*E*I/(L2),2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/(L3),-6*E*I/(L2),0,12*E*I/(L3),-6*E*I/(L2);0,6*E*I/(L2),2*E*I/L,0,-6*E*I/(L2),4*E*I/L,Beam2D2Node_Stiffness(E,I,A,L),functiony=Beam2D2Node_Assemble(KK,k,i,j)DOF(1)=3*i-2;DOF(2)=3*i-1;DOF(3)=3*i;DOF(4)=3*j-2;DOF(5)=3*j-1;DOF(6)=3*j;forn1=1:6forn2=1:6KK(DOF(n1),DOF(n2)=KK(DOF(n1),DOF(n2)+k(n1,n2);endendy=KK;,2.Beam2D2Node_Assemble(KK,k,i,j),单元刚度矩阵和全体刚度矩阵的关系,functionforces=Beam2D2Node_Forces(k,u)forces=k*u;,3.Beam2D2Node_Forces(k,u),实践课题1,用matlab编制平面2D梁单元的计算程序。用编制的matlab程序计算图示2D框架结构固定点C、D的弯矩，点A、B的弯矩及转角，载荷作用点的位移；并用材料力学理论求解，比较两者的结果。E=20000kgf/mm2，A1=2000mm2，A2=1500mm2，I1=108mm4，I2=107mm4.,P=4000N,A,B,C,D,5m,6m,I1，A1,I2，A2,I2，A2,题1图,题2图,要求：1.按照报告模板提交纸质文档及电子文档；2.提交截止日期：2011.10.13,functionk=Element_Stiffness(E,I,A,L,alpha)k=E*A/L,0,0,-E*A/L,0,0;0,12*E*I/(L3),6*E*I/(L2),0,-12*E*I/(L3),6*E*I/(L2);0,6*E*I/(L2),4*E*I/L,0,-6*E*I/(L2),2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/(L3),-6*E*I/(L2),0,12*E*I/(L3),-6*E*I/(L2);0,6*E*I/(L2),2*E*I/L,0,-6*E*I/(L2),4*E*I/L;x=alpha*pi/180;C=cos(x);S=sin(x);T=C,S,0,0,0,0;-S,C,0,0,0,0;0,0,1,0,0,0;0,0,0,C,S,0;0,0,0,-S,C,0;0,0,0,0,0,1;k=T*k*T;,functionforces=Forces(k,u)forces=k*u;,functiony=Global_Assemble(KK,k,i,j)DOF(1)=3*i-2;DOF(2)=3*i-1;DOF(3)=3*i;DOF(4)=3*j-2;DOF(5)=3*j-1;DOF(6)=3*j;forn1=1:6forn2=1:6KK(DOF(n1),DOF(n2)=KK(DOF(n1),DOF(n2)+k(n1,n2);endendy=KK;,E=2e11;A1=2e-3;A2=1.5e-3;I1=1e-4;I2=1e-5;L2=5;L1=3;k1=Element_Stiffness(E,I1,A1,L1,0);k2=Element_Stiffness(E,I2,A2,L2,90);KK=zeros(15,15);KK=Global_Assemble(KK,k1,1,5);KK=Global_Assemble(KK,k1,5,2);KK=Global_Assemble(KK,k2,3,1);KK=Global_Assemble(KK,k2,4,2);k=KK(1,2,3,4,5,6,13,14,15,1,2,3,4,5,6,13,14,15);p=0;0;0;0;0;0;0;-4000;0;u=kp;U=u(1:6);0;0;0;0;0;0;u(7:9);P=KK*U;p3=k2*0;0;0;U(1:3);p4=k2*0;0;0;U(4:6);,实践课题1-（2）的输入步骤：,C点:M3=-289.2NMD点:M4=289.2NMA点:M1=-580.8NM1=-0.000363radB点:M2=580.8NM2=0.000363radP点位移:u5=0v5=-0.0008298M,实践课题1-（2）的结果：,3.单元形状函数矩阵与刚度矩阵的性质,形状函数矩阵与刚度矩阵在有限元方法中占有最重要的位置，同时它们也具有非常明确的物理意义，分析和了解它们的性质对于更深层次地掌握有限元方法具有重要的作用。以一维杆单元为例进行讨论，其结论完全可以推广到一般单元。,一维杆单元,u1、u2为节点位移，N1、N2为对应于节点1和节点2的形状函数,（1）考虑单元一端发生单位位移，而另一端固定时的情形：u1=0，u2=1,单元形状函数性质1,0/1性质,Ni表示在i点的节点位移为1，其它节点位移为0时的单元位移场函数,1.单元形状函数矩阵,（2）考虑单元发生刚体位移的情形,设单元有刚体位移0u，由于是刚体位移，则单元的位移场函数及节点位移,单元形状函数性质2,和1性质,单元的形状函数满足：在单元的任意点处,其中n为单元的节点数，该性质表明形状函数能够描述单元的刚体位移,考察：既具有线位移又具有角位移的梁单元作刚体运动的过程,case1：沿垂直方向的刚体平动,设梁单元在垂直方向的刚体平动量为0v，此时的梁单元的位移函数为,case2：刚体转动,若绕节点1有刚体转动，此时的梁单元的位移函数为,不满足杆单元问题中,case3：一般性刚体运动（平动和转动）,如果绕节点1有刚体转动,并且在垂直方向的刚体平动量为，此时的梁单元的位移函数为,平面梁单元的形状函数分成两个部分：,对应于连续部分的形状函数,对应于连续部分的形状函数,C0型单元与C1型单元,C0型单元：单元的节点参数中仅包含场函数的节点值，即在泛函(势能)中位移函数出现的最高阶导数是1阶，在单元交界面上具有0阶的连续导数，即节点上只要求位移连续。一般的杆单元、平面问题单元、空间问题单元都是C0型单元。,C型单元：单元的节点参数中除场函数的节点值外，还包含场函数的导数的节点值，即泛函(势能)中位移函数出现的最高阶导数是2阶，在单元交界面上具有1阶的连续导数，即节点上除要求位移连续外，还要求1阶导数连续。梁单元、板单元、壳单元都是C型单元。,2.单元刚度矩阵,一维2节点杆单元,（1）单元刚度矩阵性质1：对角线元素恒为正性质,k11为使节点2的位移为零，使节点1产生单位位移而需要在节点1上所施加的力。,单元刚度矩阵的对角线元素kii表示要使单元的第i个节点产生单位位移（ui=1），而其它节点位移为0时，需在节点i所施加的力,考虑单元左端发生单位位移，而右端固定时的情形：u1