应用有限元分析工程实例-1

应用有限元分析工程实例高级研修班讲义梁醒培2012年9月 致谢本讲义参考了不少文献资料 在此向有关作者一并表示感谢 有限元法发展迅速 成果丰富 相比之下本讲义所含内容有限 向没有能够反映其成果的作者表示歉意 河南工业大学梁醒培2012年9月 主要参考文献 1 王勖成 邵敏编著 有限元法基本原理和数值方法 M 北京 清华大学出版社 1999 2 美 K J Bathe著 付子智译 工程分析中的有限元法 M 北京 机械工业出版社 1991 3 英 O C 监凯维奇著 有限元法 M 北京 科学出版社 1985 4 赵经文 王宏钰编著 结构有限元分析 M 北京 科学出版社 2001 5 梁醒培 王辉编著 应用有限元分析 M 北京 清华大学出版社 2010 6 梁醒培 王辉著 基于有限元法的结构优化设计 原理与工程应用 M 北京 清华大学出版社 2010 7 博弈创作室编著 ANSYS9 0经典产品高级分析技术与实例详解 M 北京 中国水利水出版社 2005 8 ANSYS程序电子文档本讲义系为学习班专门编写 未经许可 不得翻印 梁醒培2012 8 25 第一讲1绪论 1有限元法简介1 1有限元法发展概况 1 起源 飞机的结构分析 用三角形单元求得了平面应力问题的解答中国科学院冯康教授1960年克拉夫 Clough 第一次提出 有限单元法 这一名称数学家和工程师共同奠定了有限元法的理论基础和框架 2 总体发展概况 国外20世纪70年代就有了商业有限元程序系统国内 1980年以前高校 科研院所 大型厂矿企业开始研究 特点是功能单一 重复 不能解决较为复杂的工程问题 1980 1990大型有限元程序的开发和引进同时开展国内自主开发 大连理工大学 北京大学 原机械部郑州机械研究所 航天部的研究所等 1979年引进了线性有限元程序SAP 1980年引进了非线性有限元程序ADINA ADINAT 目前我国有限元软件市场由国外商业有限元软件占据 3 技术发展概况到1990年 有限元理论体系 常用计算功能基本定型 1980年以后 商业软件的发展主要反映在前后处理功能 基于有限元系统本身 特殊功能的研发 1990年以后 随着三维CAD软件的发展 CAE与CAD结合 使有限元前处理功能极大地提高 缩短了有限元计算的周期 1 2有限元法的基本思路及其求解步骤 1 基本思路弹性力学经典解法可以求解的问题是简化的 有限的 实际工程问题一般是复杂的 例1 焊接接头变形和应力 例2钢板筒仓变形与应力计算 2 位移有限元法的求解步骤1 离散化2 选择单元位移函数3 单元特性分析4 外载荷处理5 建立节点上的力平衡方程6 处理边界条件 解算节点位移7 后处理 1 3国内外常用有限元软件 ANSYS ABAQUS ADINA 结构简介 2有限元分析基本原理 平面三角形单元 以三角形单元为例讲解有限元法的基本原理2 1平面线弹性问题基本方程两种平面问题 平面应力问题和平面应变问题2 1 1平面应力的基本方程1 定义构件厚度远小于其他两个方向的尺寸 载荷只作用在平行于xy的平面内 沿厚度均匀分布 在z t 2的表面上 板薄 含内部 2 待求的未知函数 位移 应力和应变 3 位移 应力和应变只的函数 与坐标无关 平面应力问题的基本方程位移分量2个 应力分量3个 应变分量3个 几何方程物理方程 平衡方程以应变表示应力的物理方程令 2 1 2平面应变问题的基本方程1 定义 1 纵向 z向 尺寸于远大于其他两个方向的尺寸 截面形状沿z轴不变 两端限制了其纵向位移 2 载荷沿z轴不变 且垂直于z轴 推断 除两端附近外 没有z向位移 因此 只需研究任一横截面 因为任一横截面的情况相同 为对称面 横截面内只有x y向位移 且只是的函数 因为所以 2 待求的未知函数 位移 应变 应力 3 位移 应力和应变只是的函数 与坐标无关 平面应变问题的基本方程位移分量2个 应力分量4个 应变分量3个 几何方程物理方程 平衡方程令 2 8 可写为 由胡克定律2 1 3平面线弹性问题小结1 对于两种平面问题 除弹性常数不同外 基本方程是相同的 所以 在有限元法中可以用同一个程序计算 只要注意区分两种平面问题的弹性常数 即弹性矩阵即可 2 这两种平面问题 都有8个独立的未知量 函数 2个位移分量 函数 3个应力分量 函数 3个应变分量 函数 3 用来求解的方程也是8个 即3个几何方程 3个物理方程 2个平衡方程 由上述8个方程 可以求解8个未知量 这实质上是微分方程的积分问题 4 要完全确定解的结果 还需要满足给定的位移边界条件和载荷边界条件 2 2平面问题结构离散平面有限元计算 将求解域划分为若干个单元 三角形单元 四边形等 2 2 1单元位移函数任意三角形单元节点逆时针排序每个节点2个自由度每个单元共6个自由度 单元节点位移列阵 通过插值方法 用单元的节点位移来表示单元内任意点的位移单元位移函数可选用坐标的一次多项式 即是待定系数 将节点的坐标代入 2 12 有 将做未知量 求解式 2 13 A是三角形单元的面积 同理 利用3个节点在y方向的位移条件 即将节点的节点坐标代入式 2 12 的第二个方程 可得将看做未知量 求解可得 在式 2 15 2 16 中 将带入式 2 12 经整理可得单元位移函数 式中 称为单元形状函数 简称形函数 它们是坐标的一次函数 是取决于三个节点坐标的常数 现在 三角形单元的面积又可以表示为将式 2 18 表示为矩阵形式 则有形函数矩阵 2 在单元中任意一点上各个形函数之和等于1 即3节点三角形单元 单元位移 形函数 是线性的 相邻单元在公共边界上位移是连续的 因为公共边界上的位移可由两个节点的位移唯一地确定 能保证全域位移的连续性 连续性 形函数性质 1 在其自身节点上 1 在其他节点上 0 2 2 2单元应变和单元应力由几何方程 2 2 求得单元内任意点处的应变其中 若令称为单元应变矩阵 单元应变可写为 参数由单元的节点坐标确定 因此 仅取决于单元形状 当单元的节点坐标确定后 这些参数都是常量 所以 3节点三角形单元的应变矩阵是常数矩阵 由物理方程 可得单元的应力其中 称为单元应力矩阵平面应力 平面应变 对于3节点平面三角形单元 应力矩阵也是常数矩阵 各点的应变和应力都是相同的 单元是常应变单元 也是常应力单元 在应力变化剧烈或应力梯度较大的部位 单元划分应该适当加密 2 2 3单元刚度矩阵位移有限元法是在节点上建立力的平衡方程 即建立 单元节点力与单元节点位移之间的关系 结构在载荷作用下产生变形和应力 于是在各单元之间就产生相互作用 实际上各单元之间的相互作用是通过相邻边界上 即单元的边 实际是面 的分布力而产生的 但是 按照有限元法 人为地把结构离散化为一个个单元后 单元之间的相互作用就由单元的节点力来实现 即用单元节点力等效代替相邻边界上的相互作用力 那么 节点力就与单元应力相关 而单元应力又与节点位移相关 因此 单元节点力与单元节点位移相关 虚位移原理 在外力作用下处于平衡状态的变形体 当发生约束允许的任意微小的虚位移时 外力在虚位移上所作之虚功等于整个体积内的应力在虚应变上所作之虚功 实位移 外力F产生的位移实应力 外力F产生的应力虚位移 约束允许的任意位移虚应变 外力F所做虚功 实应力所做虚功 由 单元节点位移 单元节点力 单元应变 应力 把一个单元作为分析对象 把节点力看作外力 单元节点力和单元节点位移之间的关系可由虚位移原理导出 节点虚位移 单元虚位移 单元虚应变 节点力 外力 在虚位移上所作之虚功单元应力在虚应变上所作的虚功 令 单元的虚功方程则有表示单元节点力与节点位移关系的单元刚度方程称为单元刚度矩阵 节点力不是结构上的外载荷 而是按虚位移原理把单元边界上的分布力近似等效到单元节点上的一种节点力 节点力在实际结构中是不存在的 式 2 48 和 2 49 推导过程中所基于的原理和方法具有普遍性 原则上说 式 2 48 是位移法有限元分析中普遍适用的单元刚度阵表达式 对于不同的单元 只是其中的具体计算细节不同 一般情况下 应变矩阵是坐标的函数矩阵 对于这里使用的三角形常应变单元 因为是常数矩阵 所以 如果材料是线性 均质的 矩阵也是常数矩阵 且单元厚度是常量 这样 三角形常应变单元的刚度矩阵就可以写成A是单元的面积 t是单元的厚度 单元刚度矩阵的特性 1 对称性 单元刚度矩阵是对称矩阵 2 奇异性 单元刚度矩阵是奇异矩阵 它不存在逆矩阵 3 主元恒正 单元刚度矩阵对角元素的数值恒大于零 这可由式看出 4 单元刚度矩阵的元素具有明确的物理意义 比如 对于图2 4所示的单元 其刚度矩阵第一列的6个元素的物理意义是 当单元的第1个节点位移 节点的 为1 而其他节点位移全为零时 需要在6个节点位移方向上施加的节点力的大小 三角形常应变单元 单元刚度矩阵是6 6的对称矩阵 其显式形式如下 例2 1平面应力直角三角形单元ijm 直角边长分别为a b 厚度为t 弹性模量为E 泊松比为 求单元的刚度阵 三角形面积为刚度阵中的 单元刚度矩阵 当单元为等腰直角三角形单元时 其单元刚度矩阵与单元的边长无关 对于平面问题 如果材料的物理特性相同 则几何形状相似的单元都具有相同的单元刚度矩阵 在程序中利用这一特点 可以减少计算量 当a b时 为等腰直角三角形单元 其刚度矩阵为 a b c d 讨论 a b c d 三角形ABC和三角形abc的单元刚度阵是否相同 四边形ABCD和四边形abcd的单元刚度阵是否相同 在图 a 中 若边AB ab固定 边BC和边bc仅限制其竖直方向位移 并在点C和点c施加同样的力F 假定两个三角形材料相同 并按线弹性考虑 则 1 点C和点c的水平位移是否相等 若不相等 哪个点的位移大 为什么 2 三角形ABC和三角形abc的位移场和应力场是否相同 2 2 4单元等效节点载荷作用在结构上的外载荷 集中力 表面力 体积力等 集中力作用点划分分节点 表面力和体积力采用虚位移原理移植到单元的节点上 移植后的节点载荷和移植前的载荷在约束允许的任意虚位移上所做的虚功相等 设单元的虚位移 1 体积力 自重 离心力等 特点 体积力分布在整个单元上 令单位体积的力令与体积力等效的节点载荷 等效节点载荷的虚功 移植前载荷的虚功 体积力可以是常数 也可以是坐标的函数 如果体积力是重力 且重力方向为负y方向 则 单元体积力等效到单元节点上的计算公式为 2 表面力表面力作用在单元边界上 令表面力矢量为 令与表面力等效的单元节点载荷为 等效节点载荷的虚功 移植前表面力的虚功单元表面力的等效节点载荷计算公式积分沿作用力边界 表面力可以是沿单元边界的均布力 非均布力均布力 以右图所示为例是边的长度 相当于将表面力均分到两个节点上 非均布力 以右图所示为例表面力可用插值函数表示 对于沿单元边界呈非线性变化的表面力 同样可以采用这种线性化的处理方法 因为当单元尺寸变小时 线性化的处理产生的误差也会减小 在有限元程序中 计算也由程序自动完成的 2 2 5总体平衡方程的建立二种方法 直接组装法 略 详见文献 5 按最小势能原理来建立有限元平衡方程在外力作用下 产生的实际位移使系统的势能的一阶变分为零 即使势能取驻值 即平面线弹性问题的最小势能原理 是系统的总势能 或称泛函总势能 应变能 体积力势能 外力势能 其中 是厚度 是弹性矩阵 势能原理适用于解析解 势能原理适用于解析解以拉伸杆为例总势能应变能外力势能总势能当求解域划分为若干单元时 求解域为全部单元之和 边界为有外力作用的边界之和 一个单元的应变能 是单元面积 全部单元的应变能则为 体积力和表面力的势能是总体节点载荷与总体节点位移之乘积 即 离散系统的总势能 节点位移的函数 求变分 转化为多元函数的极值条件 2 3 2稀疏性和带状分布总体刚度矩阵的几个主要特点和性质对称性 奇异性 稀疏性 主对角元素恒正 非零元素呈带壮分布 半带宽最大半带宽 一个节点的自由度 最大节点号差图a 图b 计算量与最大半带宽的平方成正比 2 3总体刚度矩阵的物理意义及特点2 3 1总体刚度矩阵的物理意义及特点各元素的物理意义与单元刚度矩阵类似 结构的第个节点位移 不是第个节点的位移 为1而其他节点位移全为零时 需要在第个节点位移 不是第个节点的位移 方向上施加的节点力 计算量与最大半带宽的平方成正比 对于上页图a和图b 图b的节点编号法 其计算量将是图a的1 56 102 82 倍 由于总体刚度矩阵是对称矩阵 在求解过程中 只有在带宽以内的元素发生变化 所以 只需存储包括对角线元素的对角线以上 或以左 的半带宽以内的元素 这时存储的总体刚度矩阵称为上 下 三角矩阵 对于图a所示的网格 需要保存106个元素 对于图b所示的网格 需要保存122个元素 同一个问题 虽然网格相同 但是编号方法不同 需要存储的元素个数也不相同 最终的计算时间也不相同 因此 图a所示的网格要优于图b 从图中可以看到总体刚度矩阵中元素的对称性 稀疏性 以及呈带状分布的特点 因为这里的计算模型只有10个节点 如果节点数增多 其稀疏性和呈带状分布的特点会更加突出 2 4位移边界条件的处理及总体平衡方程求解2 4 1位移边界条件的处理总体平衡方程 总体方程是节点上的力的平衡方程 形成总体平衡方程的目的是求解节点位移 在形成的总体平衡方程中 只处理了载荷边界条件 而并未涉及位移边界条件 组装后的总体刚度矩阵是奇异的 从力学概念上说 这时的总体平衡方程还没有排除刚体位移 为排除刚体位移 至少要给出消除刚体位移所需的位移边界条件 否则 对于静力问题 此时的总体平衡方程无法求解 在总体有限元平衡方程中 各个方程对应的是节点的位移 而单元已不再显现 所以 位移边界条件的处理也就是针对节点位移的 位移边界条件一般分为两种 即固定约束边界和已知位移边界 对于处在这两种边界上的节点 其位移是零位移或已知的非零位移 在总体有限元平衡方程中 为已知量 为待求量 在程序中 此时是一个空数组 1 零位移边界条件的处理 1 划行划列法把零位移对应的方程的行和列直接去掉 并指定其位移为0 比如 例2 2的方程 p 划去方程第5 8并令 2 对角线约束改1法 2 非零位移边界条件的处理 乘大数法如果与第个方程对应的节点位移是 即 可采用乘大数M法 比如取M 1020 2 4 2总体刚度方程的求解 当正确处理了位移边界条件以后 总体刚度矩阵的奇异性就被消除 求解代数方程组的方法很多 比如 Gaussian消元法 LU分解法 波前法等 对于静力问题 不管用以上哪种方法求解 只有消除了刚体位移以后 才能进行求解以求得节点位移 有了节点位移就可以按和求得单元的应变和应力 例2 2 设正方形的边长为a 弹性模量为E 波松比为u 厚度为t 将板划分为2个单元 见图 节点3 4处在固定边上 令其位移为零 令节点2的y方向位移为零 节点1的作用力 求节点1 2的位移 例2 3计算具有圆孔的平板在单向拉伸下的应力集中系数 边长为140mm 板厚为10mm的方板 圆孔直径20mm 泊松比为0 3 在方板对边施加单向均布拉力 求圆孔的应力集中系数 146个单元 92个节点文件 w2 3 txt理论 最大应力在节点1 理论应力集中系数为3 0 有限元算得单元1的应力为 应力集中系数为2 94 只要网格密度合适 三角形常应变单元也能较好的模拟应力集中处的应力 2 5有限元解的收敛条件简单讲 收敛条件是指 在什么条件下 当单元尺寸趋于零时 有限元的解趋于真解 协调性协调性是指单元位移函数要保证位移在单元内及在相邻单元之间的面上连续 即在单元内及在相邻单元之间不重叠 不开裂 选择多项式建立单元位移模式时 协调性在单元内总是能满足的 而在相邻单元之间并非是任意一个多项式就能够满足的 通常 当单元交界面上的位移取决于该交界面上的节点位移时 协调性是可以满足的 完备性完备性要求单元位移函数能表示刚体位移和常应变状态 刚体位移就是单元象刚体一样移动而不产生应力 单元位移函数必须具有常数项 从物理意义上讲 当单元尺寸无限缩小时 每个单元的应变都应该趋于常值 如果不含常应变 单靠缩小单元尺寸也不可能收敛于真解 所以 单元位移函数应包含常应变项 平面3节点三角形单元 其位移函数是完全的一次多项式 所以 满足完备性要求 是完备的协调性