第5章 杆单元和梁单元.ppt

第五章杆单元和梁单元 第5章杆单元和梁单元 本章主要介绍利用杆单元及梁单元进行结构静力学的有限元分析原理 首先介绍了杆单元的分析方法 详细给出了采用杆单元进行有限元分析的整个过程 紧接着介绍了平面梁单元 以一个平面悬臂梁力学模型为分析实例 分别采用材料力学 弹性力学解析计算以及有限元法进行了分析与求解 以加深读者对有限元法的理解 杆单元 桁架结构 梁单元 轴系 转子动力学 5 1杆件系统的有限元分析方法 一般情况下 认为杆件只承受轴向力 只有一个方向的受力和相应的变形 本节将采用有限元法来分析杆件系统 以下给出规范的有限元法中关于杆单元的推导过程 以及整个杆系的求解过程 如图5 1所示的杆件结构 左端铰支 右端作用一个集中力 相关参数如图 具体求解过程如下 图5 1杆件结构 待求解的问题 1 确定坐标系 单元离散 确定位移变量 外载荷及边界条件 5 1 1 一维杆单元 材料力学可轻易求解 5 1 1一维杆单元 要建立两种坐标系 单元坐标系 局部坐标系 整体坐标系 根据自然离散 坐标系建立成一维 单元划分为两个 给出相应的节点1 2 3以及相应的坐标值 见图5 1 在局部坐标系中 取杆单元的左端点为坐标原点 图5 2为任取的一个杆单元 图5 2杆单元 对于两个节点的杆单元 存在如下节点力和节点位移的关系式 5 1 其中 称为单元刚度矩阵 5 1 1一维杆单元 2 确定位移模式 假设单元位移场 取其线性部分 系数 可由节点位移 确定 称为位移插值模式 interpolationmodel 5 2 3 形函数矩阵的推导 由单元的节点条件 两个节点坐标为x1 x2 两个节点位移为 代入上式插值模式公式得 求解得到 5 1 1一维杆单元 这样 可以写成如下矩阵形式 导出 5 3 得到形函数矩阵 shapefunctionmatrix 5 4 记节点位移矢量 nodaldisplacementvector 是 5 5 5 1 1一维杆单元 因此 用形函数矩阵表达的单元内任一点的位移函数是 5 6 4 应变 由弹性力学的几何方程知1维杆单元满足 5 7 5 应力 由弹性力学的物理方程知 5 8 6 利用最小势能原理导出单元刚度矩阵 单元的势能表达式 B为应变矩阵 常应变 S为应力矩阵 常应力 5 1 1一维杆单元 上式记作如下矩阵形式 5 9 根据最小势能原理 可以得到 5 10 5 1 1一维杆单元 7 把所有单元按结构形状进行组集 assemblyofdiscreteelements 对于图5 1所示结构 第一个单元 其中 单元刚度矩阵 elementstiffnessmatrix 或称单元特性矩阵 elementcharacteristicmatrix 5 11 5 1 1一维杆单元 整体结构的总势能是所有单元的势能的和 即 第二个单元 在这里 把表达成整体位移矢量的函数 如下 5 1杆件系统的有限元分析方法 可记作 5 12 上式的即为整体刚度矩阵 即根据最小势能原理 由各单元刚度矩阵求出的整体刚度矩阵 下式是由整体刚度矩阵表达的系统方程 5 13 5 1杆件系统的有限元分析方法 8 引入边界条件 Treatmentofboundaryconditions 为获取许可位移场 需引入边界条件 5 14 由于 可划去它所对应的行和列 这样基于许可位移场的系统总势能为 5 1 1一维杆单元 由最小势能原理 势能函数对未知位移求变分 满足的条件是 得如下方程式 9 求解节点位移 由上式方程可以直接求解得到 注意到R2是内力 不做功 在求解过程中 可以视为0 也就是 5 15 5 1 1一维杆单元 10 求单元应变 5 16 5 17 11 各单元应力 利用物理方程 求单元的应力 5 18 5 1 1一维杆单元 12 各支点反力 各支反力公式是由单元最小势能原理得到的 即 5 19 为了清楚起见 将上述两杆结构代入具体数值 进行相应的单元应力计算 得到的结果如下 这里从坐标变换的角度出发来说明平面杆单元的建立方法 从势能的角度 杆单元的势能不会因坐标系变换而产生能量的变化 如图5 3所示 局部坐标系中杆单元1维位移 和 可以投影到整体坐标系中变成 两个节点的坐标变为 图5 2局部坐标与整体坐标的变换 5 1杆件系统的有限元分析方法 5 1 2平面杆单元 整体坐标系下的位移和局部坐标系下的位移的变换关系为 5 20 式中 坐标变换矩阵为 5 21 5 1 2平面杆单元 因此 平面杆单元节点位移矢量的变换关系记为 5 22 单元势能是一个标量 不会因坐标系的不同而改变 导出整体坐标系下的单元势能函数 5 23 从上式我们可以导出在整体坐标下平面杆单元的刚度矩阵 5 24 5 1 2平面杆单元 具本而言 在转换矩阵中 有 5 25 式中为单元的长度 这里令 则转换矩阵T可表示为因此 整体坐标系下的单元刚度矩阵为 5 26 5 27 5 1 2平面杆单元 5 1 2平面杆单元 例5 1 如图5 4所示的平面桁架结构 在点1处施加有10000N向下的力 材料的弹性模量为E 30 106Pa 所有杆的横截面面积A 2m2 点2 3和点4完全约束 试求解各点的位移 受力以及各个单元的应力情况 clearall clearmemory E modulusofelasticity A areaofcrosssection L lengthofbarE 30e6 A 2 EA E A generationofcoordinatesandconnectivitiesnumberElements 3 numberNodes 4 elementNodes 12 13 14 nodeCoordinates 00 0120 120120 1200 xx nodeCoordinates 1 yy nodeCoordinates 2 GDof 2 numberNodes GDof totalnumberofdegreesoffreedomdisplacements zeros GDof 1 force zeros GDof 1 appliedloadatnode2force 2 10000 0 120 120 5 1 2平面杆单元 Ge zeros 4 GDof numberElements K zeros GDof GDof Te fori 1 numberElementspos 1 elementNodes i 1 pos 2 elementNodes i 2 节点编号Le i sqrt xx pos 1 xx pos 2 2 yy pos 1 yy pos 2 2 Ke E A Le i 1 1 11 一维杆单元Te i xx pos 2 xx pos 1 Le i yy pos 2 yy pos 1 Le i 00 0 0 xx pos 2 xx pos 1 Le i yy pos 2 yy pos 1 Le i Ke b i Te i Ke Te i Ge 1 2 pos 1 1 i 1 组集矩阵Ge 2 2 pos 1 i 1 组集矩阵Ge 3 2 pos 2 1 i 1 组集矩阵Ge 4 2 pos 2 i 1 组集矩阵K K Ge i Ke b i Ge i end 5 1 2平面杆单元 K s K 1 2 1 2 添加边界后的刚度矩阵force force 1 2 外载荷化简x inv K s force 得到位移结果X x 0 0 0 0 0 0 扩展为完整的节点位移fore 1 numberElementspos 1 elementNodes e 1 pos 2 elementNodes e 2 sigma e E 11 Te e X 2 pos 1 1 X 2 pos 1 X 2 pos 2 1 X 2 pos 2 Le e end 利用该程序求得该桁架在各节点处的位移为 U 0 0041 0 0159000000 T各单元应力 1 0e03 3 96451 4645 1 0355 T 节点支反力 在5 1节及5 2节中 考虑了一维及平面杆单元的情况 接下来考虑空间杆单元的问题 如图所示 局部坐标系中杆单元1维位移 可以投影到三维的整体坐标系中 变成 两个节点的坐标变为 5 1 3空间杆单元 5 1杆件系统的有限元分析方法 两者之间存在的关系是 5 28 式中 分别为杆单元在整体坐标系中与各轴的夹角 图5 5局部坐标与整体坐标的变换 5 1 3空间杆单元 在整体坐标系下 空间杆单元刚度矩阵为 根据上一节内容 令 则 5 29 5 30 5 1 3空间杆单元 5 2平面梁单元 平面悬臂梁问题的解析分析 平面梁单元的分析与求解 新的一类单元 简化单元 有着广泛的应用 5 2平面梁单元 梁单元的应用 5 2平面梁单元 梁单元的应用 5 2 1平面悬臂梁问题的解析分析 作为对照 先用经典材料力学法和弹性力学法对平面悬臂梁进行分析求解 1 平面悬臂梁的材料力学求解 一端受载荷作用的悬臂梁如图5 6 a 所示 选取坐标系如图5 6 b 任意横截面上的弯矩为 a 结构示意图 b 力学模型图5 6平面悬臂梁力学模型 5 31 5 2 1平面悬臂梁问题的解析分析 受载荷作用后梁产生变形 在xy平面内梁的轴线将变成一条曲线 即挠曲线 根据材料力学有关假设 梁弯曲的挠曲线的近似微分方程为 5 32 由这两个公式可得挠曲线的微分方程为 积分得 5 33 5 2 1平面悬臂梁问题的解析分析 引入边界条件 左侧固定端A处的转角和挠度均等于零 即当x 0时 5 34 把边界条件式代入式 4 22 得 再将所得积分常数C和D代回前式 得转角方程和挠曲线方程分别为 5 35 5 2 1平面悬臂梁问题的解析分析 将悬臂梁的右端受载荷W处的横坐标x l代入以上两式 得右端受载荷截面的转角和挠度分别为 2 平面悬臂梁的弹性力学求解 5 36 末端受集中载荷作用的平面悬臂梁的位移场可以用以下多项式表示 x方向 y方向 5 2 1平面悬臂梁问题的解析分析 梁的中性面 y 0的面 上的挠曲为 受载荷作用的悬臂梁上任何位置处的转角为 5 37 左侧悬臂处 x 0 的挠曲为 右端处 x L 受到集中载荷作用 挠曲为 该结果与材料力学中的挠曲线公式相同 5 38 梁中性面 y 0 上的转角为 左端点 x 0 为悬臂点 转角为 5 2 1平面悬臂梁问题的解析分析 受载荷作用的悬臂梁的应力场可在应变场的基础上 由弹性力学物理方程直接求出 右端点 x L 为受集中载荷点 转角为 受载荷作用的悬臂梁的应变场可由弹性力学几何方程求出 5 39 5 40 注意平面梁仅存x方向应力 5 2 2平面梁单元的推导 1 建立坐标系 进行单元离散 坐标系包括结构的整体坐标系和单元局部坐标系 利用局部坐标系进行单元分析 2 建立平面梁单元的位移模式 设一个平面梁单元有两个节点 如图5 7所示 在局部坐标系内 平面梁单元定义有6个自由度 图5 7平面梁单元模型 5 41 x1 0 x2 L 节点坐标 5 2 2平面梁单元的推导 略去轴向位移 可以设平面梁单元有如下4个自由度 5 42 对于平面梁单元 其弯曲变形的位移场可以设为下式 5 43 因此 梁的斜率是 Hermite型 5 44 位移模式写成矩阵形式 5 45 5 2 2平面梁单元的推导 3 推导形函数矩阵 直接利用节点坐标使推导简单 代入节点位移和节点坐标 有 其中 L 梁单元的长度 得到 前两个方程直接解出和 代入后两个方程 解出和 具体如下 5 2 2平面梁单元的推导 上面的推导可以写成如下矩阵形式 5 2 2平面梁单元的推导 求得 5 46 将式 5 46 代入式 5 45 用节点的位移形式重新整理 得 得到的用形函数矩阵表达的单元内任一点的位移是 5 2 2平面梁单元的推导 5 47 其中 N x 平面梁单元的形函数 节点位移向量 对于 对应于挠度的形函数N x 的具体表达式是 5 48 5 2 2平面梁单元的推导 4 推导应变 应力 根据最小势能原理导出单元刚度矩阵 在这里直接根据瑞利法 也可以导出以节点位移的形式来表达梁单元的应变能 弯曲梁的应变能是 5 49 二阶导数可由方程 5 47 决定 利用形函数位移插值公式 表示为 5 50 5 2 2平面梁单元的推导 其中 5 51 代入梁单元应变能公式 同时假设对于该单元而言是常量 得单元应变能 5 52 节点位移向量不是x的函数 上式可以写成 应变能的一般形式可以表达成 5 2 2平面梁单元的推导 5 53 其中 平面梁单元的单元刚度矩阵 即 5 54 考虑到B是x的函数 上式所有项积分后得局部坐标系下 代入了坐标值 的平面梁单元的单元刚度矩阵 5 55 5 2 2平面梁单元的推导 5 整体刚度矩阵的组集与坐标变换 a 局部坐标系向整体坐标系的转换 局部坐标系 整体坐标系 两种坐标系下的节点载荷 节点位移和单元刚度矩阵的变换关系为 其中坐标变换矩阵为 5 57 式中 是轴相对于x轴的夹角 可以证明 转换矩阵T的逆矩阵等于它的转置矩阵 所以 在整体坐标系下的单元刚度矩阵为 5 58 5 56 5 2 2平面梁单元的推导 b 进行整体刚度矩阵的组集 可以采用直接刚度法 6 引入约束条件 7 求解系统方程 得到所有的节点位移 8 进而再求出单元的应力应变等 5 2 2平面梁单元的推导 例5 2平面梁单元应用举例 设一方形截面的悬臂梁 截面每边长为5cm 长度为10m 在左端约束固定 在右端施以一个沿y轴负方向的集中力w 100N 求其挠度与转角 图5 8平面梁单元实例图 5 2 2平面梁单元的推导 将整个梁分成两个平面梁单元 求出每个单元的刚度矩阵 然后将两个单元组集成总体刚度矩阵 引入边界条件后 再求解出各节点的挠度和转角 具体的计算过程可参见以下MATLAB程序 注意到其中的有关坐标变换部分利用了直接给出的转换矩阵G 其具体含义还可参考第三章中的有关内容 clearx1 0 x2 sym L x sym x j 0 3 v x jm 1x1x1 2x1 3012 x13 x1 21x2x2 2x2 3012 x23 x2 2 mm inv m N v mm N 1xx 2x 3 inv m B diff N 2 k transpose B B Ke int k 0 L Element1 E 4 0e11 I bh 3 12 5 2e 7EI 4 0e11 5 2e 7Ke1 EI subs Ke L 5 Ke2 Ke1T eye 4 4 Ke1 T Ke1 T Ke2 T Ke2 T 5 2 2平面梁单元的推导 K1 G1 Ke1 G1K2 G2 Ke2 G2K K1 K2F 0 0 0 0 100 0 u F inv K u v1 xta1 v2 xta2 v3 xta3 v1 0 xta1 0 K 1 0 K 1 0 K 2 0 K 2 0 KX K 3 6 3 6 F 1 1 0 F 2 1 0 FX F 3 6 1 u inv KX FX 利用该程序求得悬臂梁节点1 左端点 节点2 中间点 和节点3 右端点 处的挠曲和转角为 00 0 0501 0 0180 0 1603 0 0240 T 5 2 2平面梁单元的推导 ANSYS的命令流如下 CONFIG NRES 1E5 title analysisofthebeamelement prep7 选单元ET 1 beam3 每个节点有3个自由度的粱单元 定义材料特性MP EX 1 4 0e11MP DENS 1 7850MP PRXY 1 0 3type 1b 5e 2 截面的尺寸参数h bs b hI b h h h 12 截面的惯性矩r 1 s I h 绘制模型用的是4个节点的板单元n 1 0 0n 2 5 0n 3 10 0e 1 2 自动分配单元号e 2 3d 1 all 加约束f 3 fy 100 外力 soluSolve post1 进入后处理PRNSOL dof 列表显示节点位移 应力 应变PRNSOL s compPRNSOL epel comp 5 2 2平面梁单元的推导 ANSYS左端点沿y方向位移 挠曲 0左端点绕z轴的转角 0中间点沿y方向位移 挠曲 0 05中间点绕z轴的转角 0 018右端节点沿y方向位移 挠曲 0 16右端节点绕z轴的转角 0 024 利用matlab和ansys两种方法求得的结果基本一致 利用前面提到的材料力学公式求得右端点处的挠度值为 可见 用有限单元法 材料力学方法和用ANSYS计算得到的结果基本一致 5 3空间梁单元 对于具有两个节点的空间梁单元 设其节点坐标和相应的节点力如下 未完 节点 1 5 59 节点 2 5 60 5 3 1空间梁单元的节点坐标 5 3 2空间梁单元的坐标变换 整体坐标系记为OXYZ 梁单元的局部坐标系记为oxyz 其中ox轴正方向由i端截面形心指向j端面形心 y轴和z轴是梁截面的两个相互垂直的形心主轴 见图5 9 坐标变换公式具有如下形式 5 61 图5 9空间梁单元的坐标变换 5 3 2空间梁单元的坐标变换 由局部坐标向整体坐标的位移变换公式是 5 62 节点力的变换公式是 5 63 单元刚度矩阵变换公式是 5 64 在三维空间中 设x y z是局部坐标系 X Y Z是整体坐标系 单元局部坐标系的三个坐标轴的方向余弦分别如下式 5 65 5 3 2空间梁单元的坐标变换 坐标变换矩阵的具体求算方法包括如下步骤 1 局部坐标系x轴在整体坐标系中的方向余