1、1,有限元程序设计,谷 音福州大学土木工程学院2010,梁单元，静力问题,2,1. 介绍.,框架结构，例如桁架、桥梁,受弯构件 flexural elements 梁,轴力构件 axial elements 杆,平面梁单元 plane beam element,3,2. 经典梁单元 (Bernoulli-Euler) Beam,平面梁假设 Plane-beam-assumption,中面法线在变形后仍保持和中面垂直的直法线假设,小变形理论,One-variable beam theory,几何关系,物理关系（应力应变关系）,梁在纯弯曲时的平面假设：梁的各个横截面在变形后仍保持为平 面，并仍垂直
2、于变形后的轴线，只是横截 面绕某一轴旋转了一个角度。,4,平衡方程,边界条件,or,or,where,k 曲率,M, Q 弯矩，剪力,I 惯性矩,5,最小势能原理,典型 C 1 连续问题,通常梁分析中常用2节点Hermite单元,6,其中,引入变形到最小 P , 得到,7,Pj 集中荷载;,Mj 弯矩力偶。,e.g. 对于均匀分布荷载,8,3. 铁木辛柯梁理论,对剪切变形的影响,3.1 理论,只考虑剪切变形,变形后轴线切向与变形前轴线之间的转角 ( x).,9,其中 (x) 为只考虑梁弯曲理论中的线性单元转角.,假设 : 截面上均匀分布剪应变,弯曲产生的位移：,( x) 相应给出沿着中线剪切角
3、 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-不会为零。这种现象
4、称为剪切闭锁。shear-locking,18,几种方法避免产生剪切闭锁,减缩积分 数值积分采用比精确积分要求少的积分点数 假设剪切应变 替代插值函数,举例说明,19,20,Timoshenko 梁 （采用精确积分）,21,采用缩减积分,22,23,结构离散,取杆件与杆件交点、集中力作用点、杆件与支承的交点为节点。相邻两节点间的杆件段是单元。节点编号时力求单元两端点号差最小。,24,坐标系,有限元中的坐标系有整体坐标系和局部坐标系。对于一个结构，整体坐标系一般只有一个；而局部坐标系有很多个，一个单元就有一个局部坐标。并且局部坐标系每一个单元的规定都是相同的，这样，同类型单元刚度矩阵相同。,25
5、,杆系结构单元主要有铰接杆单元和梁单元两种类型。它们都只有2个节点i、j。,约定：单元坐标系的原点置于节点i；节点i到j的杆轴（形心轴）方向为单元坐标系中x轴的正向。 y轴、z轴都与x轴垂直，并符合右手螺旋法则。对于梁单元， y轴和z轴分别为横截面上的两个惯性主轴。,26,平面桁架杆单元（2D LINK1） 空间杆单元（3D LINK8） 平面刚架，BEAM3 空间梁单元(BEAM4),27,2-D Elastic Beam,three degrees of freedom at each node,Ansys,28,BEAM3 is a uniaxial element with tensi
6、on, compression, and bending capabilities,BEAM23 2-D Plastic Beam,a uniaxial element with tension-compression and bending capabilities,This element allows a different unsymmetrical geometry at each end and permits the end nodes to be offset from the centroidal axis of the beam,BEAM54 2-D Elastic Tap
7、ered Unsymmetric Beam,29,30,3-D Elastic Beam,BEAM4 is a uniaxial element with tension, compression, torsion, and bending capabilities.,six degrees of freedom at each node,BEAM24 3-D Thin-walled Beam,The element has plastic, creep, and swelling capabilities in the axial direction as well as a user-de
8、fined cross-section.,This element allows a different unsymmetrical geometry at each end and permits the end nodes to be offset from the centroidal axis of the beam,BEAM44 3-D Elastic Tapered Unsymmetric Beam,31,32,33,BEAM188 3-D Linear Finite Strain Beam,BEAM188 is suitable for analyzing slender to
9、moderately stubby/thick beam structures. This element is based on Timoshenko beam theory. Shear deformation effects are included.,This element is well-suited for linear, large rotation, and/or large strain nonlinear applications.,34,35,BEAM189 3-D Quadratic Finite Strain Beam,BEAM189 is a quadratic
10、(3-node) beam element in 3-D.,For a description of the low-order beam, see BEAM188.,36,37,有限元程序设计方法简介,程序基本框图 1、输入基本数据（结构描述）： （1）控制数据：如结点总数、单元总数、约束条件总数等； （2）结点数据：如结点编号、结点坐标、约束条件等； （3）单元数据：如单元编号、单元结点序号、单元的材料特性、几何特性等； （4）载荷数据：包括集中载荷、分布载荷等。,38,2、单元分析,（1）各单元的bi,ci(i,j,m) ， 面积A； （2）应变矩阵B，应力矩阵S； （3）单元刚度矩阵k
11、； （4）单元等价载荷列向量F。,3、系统分析 （1）整体刚度矩阵K的组装； （2）整体载荷列阵P的形成；,K的存储；约束引入；求解,39,总刚存贮,全矩阵存贮法：不利于节省计算机的存贮空间，很少采用。Ki,j 对称三角存贮法：存贮上三角或下三角元素。 半带宽存贮法 ：存贮上三角形（或下三角形）半带宽以内的元素 。 一维压缩存贮法 ：半带宽存贮中仍包含了许多零元素。存贮每一行的第一个非零元素到主对角线元素。,40,等带宽形式,方阵形式,（1）半带宽存贮法,41,方阵存贮和半带宽存贮地址关系,半带宽计算：设结构单元网格中相邻结点编号的最大差值是d，则最大半带宽为UBW：,结点编号：欲使最大半带宽
12、UBW最小，必须注意结点编号方法，使直接联系的相邻节点的最大点号差最小。,42,举例,B = 2(4-1+1) = 8,B = 2(6-1+1) = 12,Advantages of 2D Storage 1)Space-saving; 2)Easy to be computerized Disadvantages of 2D Storage Enormous storage is required when local bandwidth is large.,43,例：计算下图半带宽。,结点数N=91，总刚K中的元素总数为：82（912）（91 2 ）=33124 最大半带宽UBW=(7+1
13、) 2=16，半带宽存储矩阵元素总数为182 16=2912，约方阵元素的8.8%。,44,（2） 变带宽存贮（一维压缩存贮）,等带宽存贮虽然已经节省了不少内存，但认真研究半带宽内的元素，还有相当数量的零元素。在平衡方程求解过程中，有些零元素只增加运算工作量而对计算结果不产生影响。如果这些零元素不存、不算，更能节省内存和运算时间，采用变带宽存贮可以实现（也称一维数组存贮） 。变带宽存贮编程技巧要求较高，程序较长。,45,对 称,方阵形式的刚度矩阵K,顶线以上零元素无须存贮，仅顶线以下元素。,46,一维数组A存贮刚度矩阵K,47,变带宽存贮：按列存贮方式。从左到右，逐列存放；对每一列，先存主对角
14、线元素，然后由下而上顺序存放，直到顶线下第一个元素为止。为避免混淆，我们把存贮K的一维数组称为A。实现变带宽存贮的关键问题是：总刚中元素Kij在一维数组A中的地址是什么？为此，需要知道主元Kii在A中的位置和相应列高hi。 主元位置：采用一个一维数组MAXA存主元在A中位置。 MAXA =1,2,4,6,10,12,16,18，22。,列高hj：第j行的左带宽。,48,从第j列的主对角线元素起到该列上方第一个非零元素为止，所含元素的个数称为第j列的列高，记为hj ；如果把第j列上方第1个非零元素的行号记为mj，则第j列的列高为hj = j - mj + 1 其实，hj就是第j行的左带宽，因而必
15、有UBW= max(hj)j=1,2, ,N,利用节点位移信息数组ID （去约束后节点位移自由度编码），可容易地确定刚度矩阵K任何一列的列高。,49,例：求图示框架结构h7=?。利用ID数组得各单元的连接数组LM（假定小号为i）,（1）ID数组 节点号：1 2 3 4,按列，遇1变0，遇0加1。,50,连接数组：1号单元：LM=0, 0, 1, 0, 0, 2 2号单元：LM=0, 0, 2, 3, 4, 53号单元：LM=3, 4, 5, 6, 7, 8 4单元：LM=0, 0, 1, 6, 7, 8,1 2 3 4,51,a) 如果ID（i, j ) = 0 则表明j号节点第i个自由度受有
16、约束。,b) 如果ID( i, j ) 0 则j号节点第i个自由度不受约束。并且， j号节点第i个位移分量在非约束节点位移列向量 f中的序号就是:ID( i, j ),52,主元在一维数组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 =7 1 + 1 = 7,53,MA
17、XA(j) = h1 + hj-2 - hj-1 + 1=（h1 + hj-2 + 1）+ hj-1= MAXA（j - 1）+ hj-1 因为永远有MAXA(1)= 1， MAXA(2)= 2 故计算主元地址的公式可写为MAXA(j+1)= MAXA(j)+ hj （2-1） 式中， j = 2，3，N；hj刚度矩阵K第j列的列高。,一维数组A的总长度(S) ，即刚度矩阵K按变带宽存贮的总存贮量S = MAXA(N+1)- MAXA(1),54,Ki,j在一维数组A中的地址,记Ki,j在一维数组A中的地址为AIJ。则由下图可知，,A I J = MAXA J + J I （2-2） 其中，I = mj，mj+1，J。,图5-12,j列,MAXA(J),AIJ,mj,
