梁的有限元分析原理.ppt

1,有限元程序设计,谷音福州大学土木工程学院2012,梁单元，静力问题,2,1.介绍.,框架结构，例如桁架、桥梁,受弯构件flexuralelements梁,轴力构件axialelements杆,平面梁单元planebeamelement,3,2.经典梁单元(Bernoulli-Euler)Beam,平面梁假设Plane-beam-assumption,中面法线在变形后仍保持和中面垂直的直法线假设,小变形理论,One-variablebeamtheory,几何关系,物理关系（应力应变关系）,梁在纯弯曲时的平面假设：梁的各个横截面在变形后仍保持为平面，并仍垂直于变形后的轴线，只是横截面绕某一轴旋转了一个角度。,4,平衡方程,边界条件,or,or,where,k曲率,M,Q弯矩，剪力,I惯性矩,5,最小势能原理,典型C1连续问题,通常梁分析中常用2节点Hermite单元,6,其中,引入变形到最小P,得到,7,Pj集中荷载;,Mj弯矩力偶。,e.g.对于均匀分布荷载,8,3.铁木辛柯梁理论,对剪切变形的影响,3.1理论,只考虑剪切变形,变形后轴线切向与变形前轴线之间的转角(x).,9,其中(x)为只考虑梁弯曲理论中的线性单元转角.,假设:截面上均匀分布剪应变,弯曲产生的位移：,(x)相应给出沿着中线剪切角xz,10,内部力,其中假设,11,实际上xz采用以下形式：,其中变量与z相关。,为了确定截面的不均匀剪应力分布，引入因素k修正剪应力：,12,其中k为与截面及泊松比相关的函数,可从弹性理论推导得到,假设变形场的整体势能为：,13,14,铁木辛柯梁单元采用两个独立变量,3.2离散公式,挠度w,截面曲率，不考虑剪切,每个单元的节点数量,Lagrange插值函数,15,16,17,挠度与转动采用了同阶的插值表示式。dw/dx与不同阶，因此，泛函中的第二项中的dw/dx-的积分，对于柔性梁（l/n趋于无穷大时）会被严重放大。除非是常数（没有弯曲变形），否则，dw/dx-不会为零。这种现象称为剪切闭锁。shear-locking,18,几种方法避免产生剪切闭锁,减缩积分数值积分采用比精确积分要求少的积分点数假设剪切应变替代插值函数,举例说明,这种高斯积分阶数低于被积函数所有项次精确积分所需要阶数的积分方案称之为减缩积分。实际计算表明：采用缩减积分往往可以取得较完全积分更好的精度。这是由于：精确积分常常是由插值函数中非完全项的最高方次要求，而决定有限元精度的是完全多项式的方次。这些非完全的最高方次项往往不能提高精度，反而可能带来不好的影响。取较低阶的高斯积分，使积分精度正好保证完全多项式方次的要求，而不包括更高次的非完全多项式的要求，其实质是相当用一种新的插值函数替代原来的插值函数，从而一定情况下改善了单元的精度。,19,基于最小位能原来基础上建立的位移有限元，其解答具有下限性质。即有限元的计算模型具有较实际结构偏大的整体刚度。选取缩减积分方案将使有限元计算模型的刚度有所降低，因此可能有助于提高计算精度。另外，这种缩减积分方案对于泛函中包含罚函数的情况也常常是必须的，用以保证和罚函数相应的矩阵的奇异性（见相应教程），否则将可能导致完全歪曲了的结果。,20,21,22,Timoshenko梁（采用精确积分）,23,采用缩减积分,24,25,结构离散,取杆件与杆件交点、集中力作用点、杆件与支承的交点为节点。相邻两节点间的杆件段是单元。节点编号时力求单元两端点号差最小。,26,坐标系,有限元中的坐标系有整体坐标系和局部坐标系。对于一个结构，整体坐标系一般只有一个；而局部坐标系有很多个，一个单元就有一个局部坐标。并且局部坐标系每一个单元的规定都是相同的，这样，同类型单元刚度矩阵相同。,27,杆系结构单元主要有铰接杆单元和梁单元两种类型。它们都只有2个节点i、j。,约定：单元坐标系的原点置于节点i；节点i到j的杆轴（形心轴）方向为单元坐标系中x轴的正向。y轴、z轴都与x轴垂直，并符合右手螺旋法则。对于梁单元，y轴和z轴分别为横截面上的两个惯性主轴。,28,平面桁架杆单元（2DLINK1）空间杆单元（3DLINK8）平面刚架，BEAM3空间梁单元(BEAM4),29,2-DElasticBeam,threedegreesoffreedomateachnode,Ansys,30,BEAM3isauniaxialelementwithtension,compression,andbendingcapabilities,BEAM232-DPlasticBeam,auniaxialelementwithtension-compressionandbendingcapabilities,Thiselementallowsadifferentunsymmetricalgeometryateachendandpermitstheendnodestobeoffsetfromthecentroidalaxisofthebeam,BEAM542-DElasticTaperedUnsymmetricBeam,31,32,3-DElasticBeam,BEAM4isauniaxialelementwithtension,compression,torsion,andbendingcapabilities.,sixdegreesoffreedomateachnode,BEAM243-DThin-walledBeam,Theelementhasplastic,creep,andswellingcapabilitiesintheaxialdirectionaswellasauser-definedcross-section.,Thiselementallowsadifferentunsymmetricalgeometryateachendandpermitstheendnodestobeoffsetfromthecentroidalaxisofthebeam,BEAM443-DElasticTaperedUnsymmetricBeam,33,34,35,BEAM1883-DLinearFiniteStrainBeam,BEAM188issuitableforanalyzingslendertomoderatelystubby/thickbeamstructures.ThiselementisbasedonTimoshenkobeamtheory.Sheardeformationeffectsareincluded.,Thiselementiswell-suitedforlinear,largerotation,and/orlargestrainnonlinearapplications.,36,37,BEAM1893-DQuadraticFiniteStrainBeam,BEAM189isaquadratic(3-node)beamelementin3-D.,Foradescriptionofthelow-orderbeam,seeBEAM188.,38,39,有限元程序设计方法简介,程序基本框图1、输入基本数据（结构描述）：（1）控制数据：如结点总数、单元总数、约束条件总数等；（2）结点数据：如结点编号、结点坐标、约束条件等；（3）单元数据：如单元编号、单元结点序号、单元的材料特性、几何特性等；（4）载荷数据：包括集中载荷、分布载荷等。,40,2、单元分析,（1）各单元的bi,ci(i,j,m)，面积A；（2）应变矩阵B，应力矩阵S；（3）单元刚度矩阵k；（4）单元等价载荷列向量F。,3、系统分析（1）整体刚度矩阵K的组装；（2）整体载荷列阵P的形成；,K的存储；约束引入；求解,41,总刚存贮,全矩阵存贮法：不利于节省计算机的存贮空间，很少采用。Ki,j对称三角存贮法：存贮上三角或下三角元素。半带宽存贮法：存贮上三角形（或下三角形）半带宽以内的元素。一维压缩存贮法：半带宽存贮中仍包含了许多零元素。存贮每一行的第一个非零元素到主对角线元素。,42,等带宽形式,方阵形式,（1）半带宽存贮法,43,方阵存贮和半带宽存贮地址关系,半带宽计算：设结构单元网格中相邻结点编号的最大差值是d，则最大半带宽为UBW：,结点编号：欲使最大半带宽UBW最小，必须注意结点编号方法，使直接联系的相邻节点的最大点号差最小。,44,举例,B=2(4-1+1)=8,B=2(6-1+1)=12,Advantagesof2DStorage1)Space-saving;2)EasytobecomputerizedDisadvantagesof2DStorageEnormousstorageisrequiredwhenlocalbandwidthislarge.,45,例：计算下图半带宽。,结点数N=91，总刚K中的元素总数为：82（912）（912）=33124最大半带宽UBW=(7+1)2=16，半带宽存储矩阵元素总数为18216=2912，约方阵元素的8.8%。,46,（2）变带宽存贮（一维压缩存贮）,等带宽存贮虽然已经节省了不少内存，但认真研究半带宽内的元素，还有相当数量的零元素。在平衡方程求解过程中，有些零元素只增加运算工作量而对计算结果不产生影响。如果这些零元素不存、不算，更能节省内存和运算时间，采用变带宽存贮可以实现（也称一维数组存贮）。变带宽存贮编程技巧要求较高，程序较长。,47,对称,方阵形式的刚度矩阵K,顶线以上零元素无须存贮，仅顶线以下元素。,48,一维数组A存贮刚度矩阵K,49,变带宽存贮：按列存贮方式。从左到右，逐列存放；对每一列，先存主对角线元素，然后由下而上顺序存放，直到顶线下第一个元素为止。为避免混淆，我们把存贮K的一维数组称为A。实现变带宽存贮的关键问题是：总刚中元素Kij在一维数组A中的地址是什么？为此，需要知道主元Kii在A中的位置和相应列高hi。主元位置：采用一个一维数组MAXA存主元在A中位置。MAXA=1,2,4,6,10,12,16,18，22。,列高hj：第j行的左带宽。,50,从第j列的主对角线元素起到该列上方第一个非零元素为止，所含元素的个数称为第j列的列高，记为hj；如果把第j列上方第1个非零元素的行号记为mj，则第j列的列高为hj=j-mj+1其实，hj就是第j行的左带宽，因而必有UBW=max(hj)j=1,2,N,利用节点位移信息数组ID（去约束后节点位移自由度编码），可容易地确定刚度矩阵K任何一列的列高。,51,例：求图示框架结构h7=?。利用ID数组得各单元的连接数组LM（假定小号为i）,（1）ID数组节点号：1234,按列，遇1变0，遇0加1。,52,连接数组：1号单元：LM=0,0,1,0,0,22号单元：LM=0,0,2,3,4,53号单元：LM=3,4,5,6,7,84单元：LM=0,0,1,6,7,8,1234,53,a)如果ID（i,j)=0则表明j号节点第i个自由度受有约束。,b)如果ID(i,j)0则j号节点第i个自由度不受约束。并且，j号节点第i个位移分量在非约束节点位移列向量f中的序号就是:ID(i,j),54,主元在一维数组A中的地址,数组MAXA的长度是K的行或列数加1（N+1）。K的任何一个主对角元在一维数组A中的地址：第j列主对角线元素Kjj在一维数组A中的地址等于前(j-1)列的列高之和加1，即,确定第j列列高的办法是：从1号单元起，对所有单元逐个进行检查。其中，与7号位移分量同在一个连接数组中的最小非零号码就是m7。显然有m7=1第7列的列高为：h7=j-mj+1=71+1=7,55,MAXA(j)=h1+hj-2-hj-1+1=（h1+hj-2+1）+hj-1