梁的有限元分析原理

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 总刚存贮 全矩阵存贮法 不利于节省计算机的存贮空间 很少采用 K i 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 91 2 91 2 33124最大半带宽UBW 7 1 2 16 半带宽存储矩阵元素总数为182 16 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 2 2号单元 LM 0 0 2 3 4 5 3号单元 LM 3 4 5 6 7 8 4单元 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 7 1 1 7 55 MAXA j h1 hj 2 hj 1 1 h1 hj 2 1