当前位置:首页 >> 建筑/土木 >>

有限元分析基础教程


有限元分析基础教程
Fundamentals of Finite Element Analysis





清华大学

2008-12

有限元分析基础教程

曾攀

有限元分析基础教程
Fundamentals

of Finite Element Analysis





(清华大学)

内容简介
全教程包括两大部分,共分 9 章;第一部分为有限元分析基本原理,包括第 1 章至第 5 章,内容有:绪论、有限元分析过程的概要、杆梁结构分析的有限元方法、连续体结构分析 的有限元方法、有限元分析中的若干问题讨论;第二部分为有限元分析的典型应用领域,包 括第 6 章至第 9 章,内容有:静力结构的有限元分析、结构振动的有限元分析、传热过程的 有限元分析、弹塑性材料的有限元分析。本书以基本变量、基本方程、求解原理、单元构建、 典型例题、MATLAB 程序及算例、ANSYS 算例等一系列规范性方式来描述有限元分析的力 学原理、程序编制以及实例应用;给出的典型实例都详细提供有完整的数学推演过程以及 ANSYS 实现过程。本教程的基本理论阐述简明扼要,重点突出,实例丰富,教程中的二部 分内容相互衔接,也可独立使用,适合于具有大学高年级学生程度的人员作为培训教材,也 适合于不同程度的读者进行自学;对于希望在 MATLAB 程序以及 ANSYS 平台进行建模分 析的读者,本教程更值得参考。 本基础教程的读者对象:机械、力学、土木、水利、航空航天等专业的工程技术人员、 科研工作者。

-1-

有限元分析基础教程

曾攀

作者简介
曾攀,男,1963 年生,海南省海口市人;1988 年在清华大学获博士学位,1988-1992 年先后在大连理工大学和西南交通大学从事两站博士后研究(领域为计算力学),为国家杰出 青年科学基金获得者(1998)、长江学者(2000)、德国“洪堡”学者(1994-1995)、“新世纪百千万 人才工程”国家级人选。现为清华大学机械工程系主任、教授、博士生导师,为《机械工程 学报》《工程力学》《塑性工程学报》等五个学术期刊的编委,为上海交通大学振动冲击噪 、 、 声国家重点实验室、 华中科技大学塑性成形模拟与模具国家重点实验室学术委员会委员, 先 后主持包括国家级重点基金项目、863 项目、霍英东基金项目等科研项目 30 多个,获教委 科技进步二等奖、机械部一等奖、北京高等教育教学成果二等奖各一项,获国家发明专利授 权两项; 已出版及翻译学术著作有: 《材料的概率疲劳损伤力学及现代结构分析原理》 (曾攀, 1993)、 《有限元分析及应用》(曾攀,2004)、 《工程中的有限元方法》(T.R.Chandrupatla,曾 攀译)(2006)、 《有限元方法第 1 卷:基本原理》(O.C.Zienkiewicz(第 5 版),曾攀译)(2008), 发表论文 100 多篇。主要从事计算力学、结构设计与分析、材料加工中的数值模拟等方面的 研究。

-2-

有限元分析基础教程

曾攀





有限元分析已经在教学、科研以及工程应用中成为重要而又普及的数值分析方法和工 具; 该基础教程力求提供具备现代特色的实用教程。 在教材的内容体系上综合考虑有限元方 法的力学分析原理、建模技巧、应用领域、软件平台、实例分析这几个方面,按照教科书的 方式深入浅出地叙述有限元方法,并体现出有限元原理“在使用中学习,在学习中使用”的 交互式特点,在介绍每一种单元的同时,提供完整的典型推导实例、MATLAB 实际编程以 及 ANSYS 应用数值算例,并且给出的各种类型的算例都具有较好的前后对应性,使学员在 学习分析原理的同时,也进行实际编程和有限元分析软件的操作,经历实例建模、求解、分 析和结果评判的全过程,在实践的基础上深刻理解和掌握有限元分析方法。 一本基础教材应该在培养学员掌握坚实的基础理论、 系统的专业知识方面发挥作用, 因 此,教材不但要提供系统的、具有一定深度的基础理论,还要介绍相关的应用领域,以给学 员进一步学习提供扩展空间, 本教程正是按照这一思路进行设计的; 全书的内容包括两个部 分,共分 9 章;第一部分为有限元分析基本原理,包括第 1 章至第 5 章,内容有:绪论、有 限元分析过程的概要、杆梁结构分析的有限元方法、连续体结构分析的有限元方法、有限元 分析中的若干问题讨论;第二部分为有限元分析的典型应用领域,包括第 6 章至第 9 章,内 容有:静力结构的有限元分析、结构振动的有限元分析、传热过程的有限元分析、弹塑性材 料的有限元分析。在基本原理方面,以基本变量、基本方程、求解原理、单元构建等一系列 规范的方式进行介绍;在阐述有限元分析与应用方面,采用典型例题、MATLAB 程序及算 例、ANSYS 算例的方式,以体现出分析建模的不同阶段和层次,引导学员领会有限元方法 的实质,还提供有大量的练习题。 本教程的重点是强调有限元方法的实质理解和融会贯通, 力求精而透, 强调学员综合能 力(掌握和应用有限元方法)的培养,为学员亲自参与建模、以及使用先进的有限元软件平台 提供较好的素材;同时,给学员进一步学习提供新的空间。 本教程力求体现以下特点。 (1)考虑教学适应性:强调对学员在数学原理、分析建模、软件应用几个方面的培养目 标要求,注重学员在工程数值方面的基础训练,培养学员“使用先进软件+分析实际问题” 的初步能力。 (2)考虑认知规律性:力求按照有限元分析方法的教学规律和认知规律,在教材中设计 了“基本变量、基本方程、求解原理、单元构建”这样的模块;并体现出有限元原理“在使 用中学习,在学习中使用”的交互式特点, 在介绍每一种单元的同时, 提供实用的 MATLAB 实际编程和数值实例;在每一章还进行要点总结,给出典型例题,以引导学员领会有限元方 法的实质,体现教材的启发性,有利于激发学员学习兴趣和便于自学。 (3)考虑结构完整性: 本教程提供完整的教材结构: 绪论、 正文、 典型例题、 基于 MATLAB 的编程算例与数值算例、具有一定深度的 ANSYS 算例、各章要点、习题、专业术语的英文 标注、关键词中文和英文索引、参考文献,便于学员查阅。 (4)内容上的拓展性:除基本内容外,还介绍了较广泛的应用领域,包括:静力结构分 析、结构振动分析、传热过程分析、弹塑性材料分析;提供了有关的典型问题的建模详细分 析过程,基本上反映了有限元分析在一些主要领域的应用状况及建模方法。 (5)编排上的逻辑性:本教程力求做到具有分明的层次和清楚的条理,在每一章中重点 突出有限元方法的思想、 数理逻辑及建模过程, 强调相应的工程概念, 提供典型例题及详解, 许多例题可作为读者进行编程校验的标准考题(Benchmark), 还提供了对应的 MATLAB 编程

-3-

有限元分析基础教程

曾攀

算例与 ANSYS 算例,特别是介绍了基于 APDL 参数化的 ANSYS 建模方法,并给出具体的 实例,力求反映有限元分析的内在联系及特有思维方式。 由于作者的水平有限,所写的教程一定会出现不少错误,敬请读者提出批评。

曾 攀 于清华大学 2008-12-11

-4-

有限元分析基础教程

曾攀

有限元分析基础教程
Fundamentals of Finite Element Analysis


第一部分 有限元分析基本原理
第 1 章 绪论
1.1 概况 1 1.2 有限元方法的历史 1 1.3 有限元分析的作用 5



第 2 章 有限元分析过程的概要 7
2.1 有限元分析的目的和概念 7 2.2 一维阶梯杆结构问题的求解 9 2.3 有限元分析的基本流程 17 2.4 有限元分析的特点 20 2.5 本章要点 22

第 3 章 杆梁结构分析的有限元方法

23

3.1 杆梁结构分析的工程概念 23 3.2 杆件有限元分析的标准化表征与算例 24 3.2.1 杆件分析的基本力学原理 24 3.2.2 局部坐标系中的杆单元描述 28 3.2.3 杆单元的坐标变换 32 3.2.4 杆单元分析的 MATLAB 程序 35 3.2.5 杆结构分析的算例 38 3.3 梁件有限元分析的标准化表征与算例 47 3.3.1 梁件分析的基本力学原理 48 3.3.2 局部坐标系中的平面梁单元 54 3.3.3 平面梁单元的坐标变换 62 3.3.4 空间梁单元及坐标变换 63 3.3.5 梁单元的常用等效节点载荷 66 3.3.6 梁单元分析的 MATLAB 程序 68 3.3.7 梁结构分析的算例 70 3.4 应用:桥梁结构的 ANSYS 参数化分析 77
I

有限元分析基础教程

曾攀

3.4.1 桥梁结构描述 77 3.4.2 基于 ANSYS 的桁架桥梁结构分析 78 3.5 本章要点 83 3.6 习题 83

第 4 章 连续体结构分析的有限元方法 89
4.1 连续体结构分析的工程概念 89 4.2 连续体结构分析的基本力学原理 89 4.3 平面问题有限元分析的标准化表征 95 4.3.1 平面问题的 3 节点三角形单元描述 95 4.3.2 平面问题的 4 节点矩形单元描述 101 4.3.3 平面问题 3 节点三角形单元的 MATLAB 程序 114 4.3.4 平面问题 4 节点矩形单元的 MATLAB 程序 116 4.4 轴对称问题有限元分析的标准化表征 118 4.4.1 轴对称问题的基本变量及方程 118 4.4.2 3 节点三角形轴对称单元(环形单元) 120 4.4.3 4 节点矩形轴对称单元(环形单元) 122 4.5 空间问题有限元分析的标准化表征 123 4.5.1 空间问题的 4 节点四面体单元描述 123 4.5.2 空间问题的 8 节点正六面体单元描述 126 4.5.3 空间问题 4 节点四面体单元的 MATLAB 程序 128 4.5.4 空间问题 8 节点正六面体单元的 MATLAB 程序 130 4.6 形状映射参数单元的一般原理和数值积分 133 4.6.1 两个坐标系之间的三个方面的变换 133 4.6.2 参数单元的三种类型 137 4.6.3 参数单元刚度矩阵计算的数值积分 137 4.7 平面问题分析的算例 143 4.7.1 平面 3 节点三角形单元分析的算例 143 4.7.2 平面 4 节点四边形单元分析的算例 151 4.8 空间问题分析的算例 155 4.8.1 空间 4 节点四面体单元分析的算例 155 4.8.2 空间 8 节点六面体单元分析的算例 161 4.9 本章要点 165 4.10 习题 166

第 5 章 有限元分析中的若干问题讨论 169
5.1 单元的节点编号与总刚度阵的存储带宽 169 5.2 单元形状函数矩阵与刚度矩阵的性质 170 5.2.1 形状函数矩阵的性质 170 5.2.2 刚度矩阵的性质 171 5.3 边界条件的处理与支反力的计算 177

II

有限元分析基础教程

曾攀

5.4 单元位移函数构造与收敛性要求 188 5.4.1 选择单元位移函数的一般原则 188 5.4.2 关于收敛性问题 189 5.4.3 位移函数构造的收敛性准则 190 5.5 C0 型单元与 C1 型单元 192 5.6 有限元分析结果的性质与节点应力的平均处理 193 5.6.1 有限元分析结果的下限性质 193 5.6.2 共用节点上应力的平均处理 195 5.7 高阶单元的构建 196 5.7.1 一维高阶单元 196 5.7.2 二维高阶单元 199 5.7.3 三维高阶单元 202 5.8 提高计算精度的 h 方法和 p 方法 204 5.9 本章要点 205 5.10 习题 205

第二部分 有限元分析的典型应用领域
第 6 章 静力结构的有限元分析 208
6.1 连续体平面问题的 MATLAB 有限元分析程序 208 6.1.1 程序原理 208 6.1.2 完整的 MATLAB 程序源代码 212 6.2 受均匀载荷方形板的有限元分析 216 6.3 自主程序开发与 ANSYS 前后处理器的衔接 222 6.4 工程应用:预应力万吨液压机机架的参数化建模与分析 228 6.4.1 模锻液压机的描述 228 6.4.2 8 万吨模锻液压机主牌坊的简化模型的有限元分析 230 6.5 习题 235

第 7 章 结构振动的有限元分析 237
7.1 结构振动分析的基本原理 237 7.1.1 结构振动分析的基本方程 237 7.1.2 结构振动的有限元分析列式 239 7.1.3 常用单元的质量矩阵 241 7.2 汽车悬挂系统的振动模态分析 243 7.3 带有张拉的绳索的振动模态分析 247 7.4 机翼模型的振动模态分析 251 7.5 习题 255

第 8 章 传热过程的有限元分析 258
III

有限元分析基础教程

曾攀

8.1 传热过程分析的基本原理 258 8.1.1 传热过程的基本方程 258 8.1.2 稳态传热过程的有限元分析列式 259 8.1.3 热应力问题的有限元分析列式 262 8.2 平面矩形板的稳态温度场分析 264 8.3 金属材料凝固过程的瞬态传热分析 267 8.4 温度变化下的结构热应力分析 271 8.5 习题 275

第 9 章 弹塑性材料的有限元分析 279
9.1 弹塑性材料分析的基本原理 279 9.1.1 弹塑性材料的物理方程 279 9.1.2 基于全量理论的有限元分析列式 282 9.1.3 基于增量理论的有限元分析列式 282 9.1.4 非线性方程求解的 Newton-Raphson(N-R)迭代法 283 9.2 三杆结构塑性卸载后的残余应力分析 284 9.3 悬臂梁在循环加载作用下的弹塑性分析 289 9.4 习题 294

参考文献 296 附录 A:MATLAB 程序基本操作 297 附录 B:ANSYS 程序基本操作 309 附录 C:常用材料的力学性能 316 附录 D:常用材料的热力学参数 317 附录 E:计量单位换算 318 中文索引 319 英文索引 323 单元及编程索引 327

IV

有限元分析基础教程

曾攀

典型例题、求解原理、MATLAB 算例、ANSYS 算例





第2章
【典型例题】2.1(1) 【典型例题】2.2(1) 【典型例题】2.2(2) 【典型例题】2.2(3) 【典型例题】2.3(1) 一个一维函数的两种展开方式的比较 1D 阶梯杆结构问题的材料力学求解 1D 阶梯杆结构的节点位移求解及平衡关系 1D 阶梯杆结构基于位移求解的通用形式 1D 三连杆结构的有限元分析过程

第3章
【基本变量】3.2.1(1) 1D 问题的基本变量 【基本方程】3.2.1(2) 1D 问题的基本方程 【求解原理】3.2.1(3) 1D 问题的直接求解 【求解原理】3.2.1(4) 1D 问题的虚功原理求解 【求解原理】3.2.1(5) 1D 问题的最小势能原理求解 【典型例题】3.2.1(6) 变截面杆单元的推导 【单元构造】3.2.2(1) 杆单元的描述 【MATLAB 程序】3.2.4(1) 1D 杆单元的有限元分析程序(Bar1D2Node) 【MATLAB 程序】3.2.4(2) 2D 杆单元的有限元分析程序(Bar2D2Node) 【典型例题】3.2.5(1) 四杆桁架结构的有限元分析 【MATLAB 算例】3.2.5(2) 四杆桁架结构的有限元分析(Bar2D2Node) 【ANSYS 算例】3.2.5(3) 四杆桁架结构的有限元分析 【基本变量】3.3.1(1) 平面梁的基本变量 【基本方程】3.3.1(2) 平面梁的基本方程 【求解原理】3.3.1(3) 简支梁的微分方程解 【求解原理】3.3.1(4) 简支梁的虚功原理求解 【求解原理】3.3.1(5) 简支梁的最小势能原理求解 【单元构造】3.3.2(1) 平面纯弯梁单元的描述 【单元构造】3.3.2(2) 一般平面梁单元的描述 【典型例题】3.3.2(3) 受均布载荷平面梁单元的等效节点载荷 【典型例题】3.3.2(4) 悬臂-简支平面连续梁的有限元分析 【MATLAB 程序】3.3.6(1) 1D 梁单元的有限元分析程序(Beam1D2Node) 【MATLAB 程序】3.3.6(2) 2D 梁单元的有限元分析程序(Beam2D2Node) 【典型例题】3.3.7(1) 三梁平面框架结构的有限元分析 【MATLAB 算例】3.3.7(2) 三梁平面框架结构的有限元分析(Beam2D2Node) 【ANSYS 算例】3.3.7(3) 三梁平面框架结构的有限元分析 【ANSYS 算例】3.4.2(1) 基于图形界面(GUI)的桁架桥梁结构分析 【ANSYS 算例】3.4.2(2) 基于命令流方式的桁架桥梁结构分析

i

有限元分析基础教程

曾攀

【ANSYS 算例】3.4.2(3)

基于参数化方式的桁架桥梁结构分析

第4章
【基本变量】4.2.1(1) 连续体问题的三大类变量 【基本方程】4.2.1(2) 连续体问题的三大类方程及边界条件 【求解原理】4.2.1(3) 直接法以及试函数法的求解思想 【求解原理】4.2.1(4) 连续体问题求解的虚功原理 【求解原理】4.2.1(5) 连续体问题求解的最小势能原理 【强度准则】4.2.1(6) 结构分析中的受力状态诊断(强度准则) 【单元构造】4.3.1(1) 平面问题的 3 节点三角形单元 【单元特征】4.3.1(2) 平面 3 节点三角形单元的位移坐标变换问题 【单元特征】4.3.1(3) 平面 3 节点三角形单元的常系数应变和应力 【单元构造】4.3.2(1) 平面问题的 4 节点矩形单元 【单元特征】4.3.2(2) 4 节点矩形单元的线性应变和应力 【典型例题】4.3.2(3) 三角形单元与矩形单元计算精度的比较 【ANSYS 算例】4.3.2(4) 三角形单元与矩形单元的精细网格的计算比较 【MATLAB 程序】4.3.3(1) 3 节点三角形单元的有限元分析程序(Triangle2D3Node) 【MATLAB 程序】4.3.4(1) 平面 4 节点矩形单元的有限元分析程序(Quad2D4Node) 【基本变量】4.4.1(1) 轴对称问题的三大类变量 【基本方程】4.4.1(2) 轴对称问题的三大类方程及边界条件 【单元构造】4.4.2(1) 3 节点三角形轴对称单元(环形单元) 【单元构造】4.4.3(1) 4 节点矩形轴对称单元(环形单元) 【单元构造】4.5.1(1) 空间问题的 4 节点四面体单元 【单元特征】4.5.1(2) 4 节点四面体单元的位移坐标变换问题 【单元特征】4.5.1(3) 4 节点四面体单元的常系数应变和应力 【单元构造】4.5.2(1) 空间问题的 8 节点正六面体单元 【单元特征】4.5.2(2) 8 节点正六面体单元的一次线性应变和应力 【MATLAB 程序】4.5.3(1) 4 节点四面体单元的有限元分析程序(Tetrahedron3D4Node) 【MATLAB 程序】4.5.4(1) 8 节点正六面体单元的有限元分析程序(Hexahedral3D8Node) 【基本原理】4.6.1(1) 两个坐标系之间的函数映射 【基本原理】4.6.1(2) 两个坐标系之间的偏导数映射 【基本原理】4.6.1(3) 两个坐标系之间的面(体)积元映射 【基本原理】4.6.2(1) 等参元、超参元以及亚参元 【基本原理】4.6.3(1) 数值积分的 Gauss 方法 【典型例题】4.6.3(2) 平面 4 节点四边形等参元的刚度矩阵的计算 【典型例题】4.7.1(1) 基于 3 节点三角形单元的矩形薄板分析 【MATLAB 算例】4.7.1(2) 基于 3 节点三角形单元的矩形薄板分析(Triangle2D3Node) 【ANSYS 算例】4.7.1(3) 基于 3 节点三角形单元的矩形薄板分析 【MATLAB 算例】4.7.2(1) 基于 4 节点四边形单元的矩形薄板分析(Quad2D4Node) 【ANSYS 算例】4.7.2(2) 基于 4 节点四边形单元的矩形薄板分析 【MATLAB 算例】4.8.1(1) 基于 4 节点四面体单元的空间块体分析(Tetrahedron3D4Node) 【ANSYS 算例】4.8.1(2) 基于 4 节点四面体单元的空间块体分析 【MATLAB 算例】4.8.2(1) 基于 8 节点六面体单元的空间块体分析(Hexahedral3D8Node)
ii

有限元分析基础教程

曾攀

【ANSYS 算例】4.8.2(2)

基于 8 节点六面体单元的空间块体分析

第5章
【基本原理】5.2.1(1) 单元形状函数性质 1:0/1 性质 【基本原理】5.2.1(2) 单元形状函数性质 2:和 1 性质 【基本原理】5.2.2(1) 单元刚度矩阵性质 1:对角线元素的 1/0 性质 【基本原理】5.2.2(2) 单元刚度矩阵性质 2:非对角线元素的 1/0 性质 【基本原理】5.2.2(3) 单元刚度矩阵性质 3:对称性质 【基本原理】5.2.2(4) 单元刚度矩阵性质 4:半正定性质 【基本原理】5.2.2(5) 单元刚度矩阵性质 5:奇异性质 【基本原理】5.2.2(6) 单元刚度矩阵性质 6:行(或列)的代数和为零的性质 【典型例题】5.2.2(7) 平面梁单元形状函数的性质 【基本原理】5.3(1) 处理边界条件的直接法 【基本原理】5.3(2) 处理边界条件的置“1”法 【基本原理】5.3(3) 处理边界条件的乘大数法 【基本原理】5.3(4) 支反力的计算 【基本原理】5.3(5) 处理耦合边界条件的拉格朗日(Lagrange)乘子法 【基本原理】5.3(6) 处理耦合边界条件的罚函数法 【典型例题】5.3(7) 平面问题斜支座的处理 【ANSYS 算例】5.3(8) 平面问题斜支座的处理 【基本原理】5.4.3(1) 收敛性准则 1:完备性要求(针对单元内部) 【基本原理】5.4.3(2) 收敛性准则 2:协调性要求(针对单元之间) 【典型例题】5.4.3(3) 平面单元位移函数选取的要求 【典型例题】5.4.3(4) 平面弯曲梁单元位移函数选取的要求 【典型例题】5.4.3(5) 平面 3 节点三角形单元的二次位移函数的选择与分析 【基本原理】5.5(1) C0 型单元的位移函数连续性 【基本原理】5.5(2) C1 型单元的位移函数连续性 【基本原理】5.6.1(1) 有限元位移结果的下限性质 【基本原理】5.6.1(2) 有限元模型的刚化性 【典型例题】5.6.1(3) 基于网格加密的求解精度估计 【基本原理】5.6.2(1) 共用节点上应力的直接平均 【基本原理】5.6.2(2) 共用节点应力的加权平均 【单元构造】5.7.1(1) 1D 高阶单元:二次杆单元 【单元构造】5.7.1(2) 1D 高阶单元:高次梁单元 【基本原理】5.7.2(1) (面积)自然坐标 【单元构造】5.7.2(2) 2D 高阶单元:6 节点三角形二次单元 【单元构造】5.7.2(3) 2D 高阶单元:矩形高阶 Lagrange 型单元 【单元构造】5.7.3(1) 3D 高阶单元:10 节点四面体二次单元 【单元构造】5.7.3(2) 3D 高阶单元:20 节点正六面体高阶单元 【基本原理】5.8(1) 提高计算精度的 h 方法(h-version 或 h-method) 【基本原理】5.8(2) 提高计算精度的 p 方法(p-version 或 p-method)

iii

有限元分析基础教程

曾攀

第6章
【MATLAB 程序】6.1.2(1) 平面问题有限元分析的通用程序 FEM2D.m 【MATLAB 算例】6.2(1) 受均匀载荷方形板的有限元分析(FEM2D.m) 【ANSYS 算例】6.2(2) 受均匀载荷方形板的有限元分析 【ANSYS 程序】6.3(1) ANSYS 前后处理器与自主程序的衔接 【ANSYS 算例】6.4.2(1) 8 万吨模锻液压机主牌坊的分析(GUI) 【ANSYS 算例】6.4.2(2) 8 万吨模锻液压机主牌坊的参数化建模与分析(命令流)

第7章
【基本变量】7.1.1(1) 结构振动的三大类变量 【基本方程】7.1.1(2) 结构振动的三大类方程及边界/初始条件 【求解原理】7.1.1(3) 结构振动求解的虚功原理 【单元构造】7.1.2(1) 结构振动分析的单元构造的基本表达式 【单元构造】7.1.3(1) 杆单元的质量矩阵 【单元构造】7.1.3(2) 梁单元的质量矩阵 【单元构造】7.1.3(3) 平面三节点三角形单元的质量矩阵 【ANSYS 算例】7.2(1) 汽车悬挂系统的振动模态分析(GUI) 【ANSYS 算例】7.2(2) 汽车悬挂系统的振动模态分析(命令流) 【ANSYS 算例】7.3(1) 带有张拉的绳索的振动模态分析(GUI) 【ANSYS 算例】7.3(2) 带有张拉的绳索的振动模态分析(命令流) 【ANSYS 算例】7.4(1) 机翼模型的振动模态分析(GUI) 【ANSYS 算例】7.4(2) 机翼模型的振动模态分析(命令流)

第8章
【基本方程】8.1.1(1) 传热过程的基本变量及方程 【求解原理】8.1.1(2) 传热过程分析的求解原理(求极值问题) 【单元构造】8.1.2(1) 稳态传热过程的单元构造基本表达式 【单元构造】8.1.2(2) 平面 3 节点三角形传热单元 【基本方程】8.1.3(1) 热应力问题中的物理方程 【求解原理】8.1.3(2) 热应力问题求解的虚功原理 【单元构造】8.1.3(3) 热应力问题分析的单元构造的基本表达式 【ANSYS 算例】8.2(1) 2D 矩形板的稳态热对流的自适应分析(GUI) 【ANSYS 算例】8.2(2) 2D 矩形板的稳态热对流的自适应分析(命令流) 【ANSYS 算例】8.3(1) 金属材料凝固过程的瞬态传热分析(GUI) 【ANSYS 算例】8.3(2) 金属材料凝固过程的瞬态传热分析(命令流) 【ANSYS 算例】8.4(1) 升温条件下杆件支撑结构的热应力分析(GUI) 【ANSYS 算例】8.4(2) 升温条件下杆件支撑结构的热应力分析(命令流)

iv

有限元分析基础教程

曾攀

第9章
【基本原理】9.1.1(1) 材料的弹塑性行为实验 【基本原理】9.1.1(2) 材料塑性行为的三方面准则 【单元构造】9.1.2(1) 基于全量理论的单元构造的基本表达式 【单元构造】9.1.3(1) 基于增量理论的单元构造的基本表达式 【求解原理】9.1.4(1) Newton-Raphson(N-R)迭代法的原理 【ANSYS 算例】9.2(1) 三杆结构塑性卸载后的残余应力计算(GUI) 【ANSYS 算例】9.2(2) 三杆结构塑性卸载后的残余应力计算(命令流) 【ANSYS 算例】9.3(1) 悬臂梁在循环加载作用下的弹塑性计算(GUI) 【ANSYS 算例】9.3(2) 悬臂梁在循环加载作用下的弹塑性计算(命令流)

v

有限元分析基础教程

曾攀

第一部分 有限元分析基本原理 第 1 章 绪论
1.1 概况

有限元方法(finite element method)或有限元分析(finite element analysis)[1][2],是求取复杂 微分方程近似解的一种非常有效的工具, 是现代数字化科技的一种重要基础性原理。 将它用 于在科学研究中,可成为探究物质客观规律的先进手段。将它应用于工程技术中,可成为工 程设计和分析的可靠工具。严格来说,有限元分析必须包含三个方面:(1)有限元方法的基 本数学力学原理,(2)基于原理所形成的实用软件,(3)使用时的计算机硬件。随着现代计算 机技术的发展,一般的个人计算机就能满足第(3)方面的要求;因此,本书的重点将在以上 的第(1)和第(2)方面, 将通过一些典型的实例来深入浅出地系统阐述有限元分析的基本原理, 并强调原理的工程背景和物理概念;基于 MATLAB 平台来系统演示基于有限元原理的编程 方法和过程;通过 ANSYS 分析平台来展示具体应用有限元方法的建模过程。

1.2

有限元方法的历史

有限元方法的思想最早可以追溯到古人的“化整为零”“化圆为直”的作法,如“曹 、 冲称象”的典故,我国古代数学家刘徽采用割圆法来对圆周长进行计算;这些实际上都体现 了离散逼近的思想,即采用大量的简单小物体来“冲填”出复杂的大物体。 早在 1870 年,英国科学家 Rayleigh 就采用假想的“试函数”来求解复杂的微分方程, 1909 年 Ritz 将其发展成为完善的数值近似方法,为现代有限元方法打下坚实基础。 20 世纪 40 年代,由于航空事业的飞速发展,设计师需要对飞机结构进行精确的设计和 计算,便逐渐在工程中产生了的矩阵力学分析方法;1943 年,Courant 发表了第一篇使用三 1956 年波音公司的 Turner, Clough, Martin 角形区域的多项式函数来求解扭转问题的论文[3]; [4] 和 Topp 在分析飞机结构时系统研究了离散杆、梁、三角形的单元刚度表达式 ;1960 年 Clough 在处理平面弹性问题,第一次提出并使用“有限元方法”(finite element method)的名 称[5][6];1955 年德国的 Argyris 出版了第一本关于结构分析中的能量原理和矩阵方法的书[7], 为后续的有限元研究奠定了重要的基础,1967 年 Zienkiewicz 和 Cheung 出版了第一本有关 有限元分析的专著;1970 年以后,有限元方法开始应用于处理非线性和大变形问题;我国 的一些学者也在有限元领域做出了重要的贡献,如胡海昌于 1954 提出了广义变分原理[8], 钱伟长最先研究了拉格朗日乘子法与广义变分原理之间关系[9],钱令希在 20 世纪五十年代 就研究了力学分析的余能原理,冯康在 20 世纪六十年代就独立地、并先于西方奠定了有限 元分析收敛性的理论基础。图 1-1 展示了有限元方法的发展过程。 随着计算机技术的飞速发展, 基于有限元方法原理的软件大量出现, 并在实际工程中发 挥了愈来愈重要的作用;目前,专业的著名有限元分析软件公司有几十家,国际上著名的通 用有限元分析软件有 ANSYS, ABAQUS, MSC/NASTRAN, MSC/MARC, ADINA, ALGOR, PRO/MECHANICA,IDEAS,还有一些专门的有限元分析软件,如 LS-DYNA,DEFORM, PAM-STAMP, AUTOFORM,SUPER-FORGE 等;国际上著名的主要有限元分析软件状况见 表 1-1。有关有限元分析的学术论文,每年也不计其数,学术活动非常活跃,表 1-2 列出的 是刊登有限元分析论文的常见学术期刊。

1

有限元分析基础教程

曾攀

图 1-1 有限元方法的发展过程

表 1-1 国际上著名的有限元分析软件状况 年份 1965 软件名称 ASKA (PERMAS) STRUDL 1966 1967 NASTRAN BERSAFE SAMCEF 1969 ASAS MARC PAFEC SESAM 1970 ANSYS SAP 开发者 IKOSS GmbH, (INTES),Germany MCAUTO, USA MacNeal-Schwendler Corp., USA CEGB, UK (restructured in 1990) Univer. of Liege, Belgium Atkins Res.&Devel., UK MARC Anal. Corp., USA PAFEC Ltd, UK now SER Systems DNV, Norway Swanson Anal. Syst., USA NISEE, Univ. of California, Berkeley, USA 1971 STARDYNE TITUS (SYSTUS) 1972 DIANA WECAN 1973 1975 GIFTS ADINA CASTEM Mech. Res. Inc., USA CITRA, France; ESI Group TNO, The Netherlands Westinghouse R&D, USA CASA/GIFTS Inc., USA ADINA R&D, Inc., USA CEA, France www.adina.com www.castem.org:8001/ HomePage.html FEAP NISEE, Univ. of California, Berkeley, USA 1976 1978 1979 1980 NISA DYNA2D, DYNA3D ABAQUS LUSAS Eng. Mech. Res. Corp., USA Livermore Softw. Tech. Corp., USA Hibbit, Karlsson & Sorensen, Inc., USA FEA Ltd., UK www.eerc.berkeley.edu/sof tware_and_data www.emrc.com www.lstc.com www.abaqus.com www.lusas.com www.dnv.no www.ansys.com www.eerc.berkeley.edu/sof tware_and_data www.reiusa.com www.systus.com www.diana.nl www.samcef.com www.wsasoft.com www.marc.com 网址 www.intes.de www.gtstrudl.gatech.edu www.macsch.com

2

有限元分析基础教程 1982 1984 COSMOS/M ALGOR Structural Res. & Anal. Corp., USA Algor Inc., USA

曾攀 www.cosmosm.com www.algor.com

表 1-2 刊登有限元分析论文的学术期刊 出版商 Academic Press Elsevier 服务商 IDEAL ScienceDirect 网址 www.idealibrary.com www.elsevier.nl 期刊名称 J. of Sound and Vibration Acta Materialia Advances in Eng. Software Applied Math. Modelling Composite Structures Composites A, B Composites Science and Technology Computational Materials Science Computer Meth. in Appl. Mech and Eng. Computers & Structures Eng. Analysis with Boundary Elements Eng. Failure Analysis Eng. Fracture Mechanics Eng. Structures European J. of Mechanics A, B Finite Elements in Analysis and Design Int. J. of Mechanical Sciences Int. J. of Solids and Structures Int. J. of Impact Eng. Int. J. of Fatigue Int. J. of Plasticity Int. J. of Non-Linear Mechanics J. of Constructional Steel Research J. of Material Processing Technology J. of the Mech. and Physics of Solids Materials Science and Eng. A Mechanics of Materials Mechanics Research Communications Nuclear Eng. And Design Probabilistic Eng. Mechanics Structural Safety Theoretical and Appl. Fracture Mech. Thin-Walled Structures IoP Publishing Kluwer Academic Publishing Electronic J. Kluwer Online www.iop.org/ EJ www.wkap.nl/ journals Smart Materials and Structures Applied Composite Materials Int. J. of Fracture Meccanica Mechanics of Composite Materials

3

有限元分析基础教程 MCB University Press Springer Verlag LINK Link.springer.de/ ol/eol/index.htm Archive of Applied Mechanics Computational Mechanics Engineering With Computers J. Wiley & Sons Interscience www.interscience. wiley.com Emerald Library www.mcb.co.uk/ portfolio.htm Engineering Computations

曾攀

Int. J. of Num. Meth. For Heat&Fluid Flow

Communicat. in Numer. Meth. in Eng. Int. J. for Numerical Methods in Eng. Mech. of Cohesive-Frictional Mater. Progress in Structural Eng. and Mater.

下面介绍几位有限元方法的创始人的情况。
Richard Courant(美国数学家,1888.01.08-1972.01.27,出生地:德国Lublinitz), 1910 年

在哥廷根大学获得博士学位,1928 年创办了应用数学研究所,并在 1928~1933 年担任 所长,开始应用数学的研究,1936 年去纽约大学创立了数学研究所,1964 年该研究所 被命名为Courant数学科学研究所,出版了数学物理方法教科书,成为最有影响的书籍 之一,还出版了科普读物“什么是数学?”,至今还再版发行,Courant的名字还与有 限元方法紧密联系在一起,是他给这种数值求解偏微分方程以坚实的理论基础, “Courant–Friedrichs–Lewy条件”以及“Courant最小值原理”都是以他的名字命名的重要定理。

John Argyris(德国人, 1913.08.19-2004.04.02,出生地:希腊Volos)是公认的计算科学、 航空航天、流体力学领域的杰出专家,被誉为有限元分析的发明人和创始者之一,所 创立的力学分析的矩阵位移方法为有限元方法的前身,被称为世界上最先进的数学工 具,先后因在有限元方法以及混沌理论方面的贡献而获得菲利普王子奖章、美国最高 荣誉的爱因斯坦奖,获得包括 5 个G8 国家在内的由 16 发达国家授予的杰出科技与工 程奖;Argyris于 1950 年代在斯图加特大学创立计算机应用研究所,并担任所长近 40 年,立至于有限元分析的研究,开发了几个著名的在欧洲广泛使用的有限元分析商业 化软件,在有限元方法的应用方面也作出了杰出的贡献。 Olgierd Cecil Zienkiewiz (英国人,1921.05.18-,出生地:英国Caterham)。 英国Swansea 大学的荣誉退休教授,是该校工程数值方法研究所的原主任,现在仍然是西班牙巴塞罗 那Calalunya技术大学工程数值方法的UNESCO主席。从 1961 至 1989 年,担任Swansea 大学土木工程系的主任,使该系成为有限元研究的重要中心之一。在 1968 年,创办了 International Journal for Numerical Methods in Engineering杂志并任主编,该杂志至今仍然 是该领域的主要刊物。他被授予 24 个荣誉学位和多种奖励。Zienkiewiez教授还是 5 所科 学院的院士,这是对他在有限元方法领域的奠基性发展和贡献的赞誉。1978 年,成为皇 家科学院和皇家工程院的院士;并先后被选为美国工程院的外籍院(1981),波兰科学院院士(1985),中国科 学院院士(1998)和意大利国家科学院院士(1999)。 1967 年, 他出版了有限元领域的第一本专著 《有限元方法》 , 该书目前也出版到第 5 版,可参见《有限元方法第 1 卷:基本原理》(O.C.Zienkiewicz, R.L.Taylor(第 5 版), 曾攀译)(清华大学出版社,2008)。

4

有限元分析基础教程

曾攀

1.3

有限元分析的作用

基于功能完善的有限元分析软件和高性能的计算机硬件对设计的结构进行详细的力学 分析, 以获得尽可能真实的结构受力信息, 就可以在设计阶段对可能出现的各种问题进行安 全评判和设计参数修改,据有关资料,一个新产品的问题有 60%以上可以在设计阶段消除, 甚至有的结构的施工过程也需要进行精细的设计, 要做到这一点, 就需要类似有限元分析这 样的分析手段。 下面举出几个涉及土木工程、车辆工程、航空工程以及生物工程的实例。 北京奥运场馆的鸟巢由纵横交错的钢铁枝蔓组成, 它是鸟巢设计中最华彩的部分, 见图 1-2,也是鸟巢建设中最艰难的。看似轻灵的枝蔓总重达 42000 吨,其中,顶盖以及周边悬 空部位重量为 14000 吨,在施工时,采用了 78 根支柱进行支撑,也就是产生了 78 个受力区 域, 在钢结构焊接完成后, 需要将其缓慢而又平稳地卸去, 让鸟巢变成完全靠自身结构支撑; 因而,支撑塔架的卸载,实际上就是对整个钢结构的加载,如何卸载?需要进行非常详细的 数值化分析,以确定出最佳的卸载方案。2006 年 9 月 17 日成功地完成了整体钢结构施工的 最后卸载。

图 1-2 北京奥运场馆鸟巢的钢铁枝蔓结构

图 1-3 列车车厢整体结构的有限元模型

图 1-3 给出的是现代列车车厢整体结构的有限元分析模型[10]; 1-4 表明的是空客 A350 图 [11] 后机身第 19 框的设计与有限元分析过程 ;图 1-5 为人体肩部区域的骨胳有限元分析模型 以及计算结果[12],并与实际测试结果进行了比较。所有这些表明了有限元方法的广泛适应 性以及在科学研究与重大工程中的重要性。

5

有限元分析基础教程

曾攀

图 1-2 空客 A350 后机身第 19 框的设计与有限元分析过程

图 1-5 人体肩部区域的骨胳有限元分析模型及计算结果

6

有限元分析基础教程

曾攀

第 2 章 有限元分析过程的概要
本章先通过一个简单的实例, 采用直接的推导方法, 逐步展示有限元分析的基本流程, 从中可以了解有限元方法的思路形成过程, 以及如何由具体的求解步骤归纳出一种通用的标 准求解方法。

2.1 有限元分析的目的和概念
任何具有一定使用功能的构件(称为变形体(deformed body))都是由满足要求的材料所 制造的,在设计阶段,就需要对该构件在可能的外力作用下的内部状态进行分析,以便核对 所使用材料是否安全可靠, 以避免造成重大安全事故。 描述可承力构件的力学信息一般有三 类: (1) 构件中因承载在任意位置上所引起的移动(称为位移(displacement)); (2) 构件中因承载在任意位置上所引起的变形状态(称为应变(strain)); (3) 构件中因承载在任意位置上所引起的受力状态(称为应力(stress)); 若该构件为简单形状,且外力分布也比较单一,如:杆、梁、柱、板就可以采用材料 力学的方法,一般都可以给出解析公式,应用比较方便;但对于几何形状较为复杂的构件却 很难得到准确的结果,甚至根本得不到结果。 有限元分析的目的:针对具有任意复杂几何形状变形体,完整获取在复杂外力作用下 它内部的准确力学信息,即求取该变形体的三类力学信息(位移、应变、应力)。 在准确进行力学分析的基础上,设计师就可以对所设计对象进行强度(strength)、刚度 (stiffness)等方面的评判,以便对不合理的设计参数进行修改,以得到较优化的设计方案;然 后,再次进行方案修改后的有限元分析,以进行最后的力学评判和校核,确定出最后的设计 方案。 图 2-1 给出一个针对大型液压机机架的设计过程以及采用有限元分析的状况。

Y

Z

X

(a)机架设计与受力状况

(b)有限元分析模型与得到的变形状况

图 2-1 大型液压机机架的设计过程与数字化分析

为什么采用有限元方法就可以针对具有任意复杂几何形状的结构进行分析,并能够得 到准确的结果呢?这时因为有限元方法是基于“离散逼近(discretized approximation)”的基 本策略,可以采用较多数量的简单函数的组合来“近似”代替非常复杂的原函数。 一个复杂的函数,可以通过一系列的基底函数(base function)的组合来“近似” ,也就是 函数逼近,其中有两种典型的方法:(1)基于全域的展开(如采用傅立叶级数展开),以及(2)
7

有限元分析基础教程

曾攀

基于子域(sub-domain)的分段函数(pieces function)组合(如采用分段线性函数的连接);下面, 仅以一个一维函数的展开为例说明全域逼近与分段逼近的特点。 【典型例题】2.1(1) 一个一维函数的两种展开方式的比较 设有一个一维函数 f ( x), x ∈ [ x0 , xL ] ,分析它的展开与逼近形式。 解答:首先考虑基于全域的展开形式,如采用傅立叶级数(Fourier series)展开,则有

f ( x) ≈ c0 ? ?0 ( x ∈ [ x0 , xL ]) + c1 ? ?1 ( x ∈ [ x0 , xL ]) + = ∑ ci?i ( x ∈ [ x0 , xL ])
i =0 n

(2-1)

其中 ?i ( x ∈ [ x0 , xL ]) 为所采用的基底函数,它定义在全域 [ x0 , xL ] 上, c0 , c1 , c2 系数。 第二种是基于子域 [ xi , xi +1 ] 上的分段展开形式,若采用线性函数,有

为展开的

f ( x) ≈ {a0 + b0 x( x ∈ [ x0 , x1 ])} + {a1 + b1 x( x ∈ [ x1 , x2 ])} + = ∑ {ai + bi x( x ∈ [ xi , xi +1 ])}
i =0 n

(2-2)

其中 ai + bi x( x ∈ [ xi , xi +1 ]) 为所采用的基底函数,它定义在子域 [ xi , xi +1 ] 上, a0 , b0 , a1 , b1 为展开的系数。 这两种函数展开的方式如图 2-2 所示。

(a)基于全域 [ x0 , xL ] 的函数展开与逼近

(b) 基于子域 [ xi , xi +1 ] 的函数展开与逼近 8

有限元分析基础教程 图 2-2 一个一维函数的两种展开方式

曾攀

比较以上两种方式的特点, 可以看出, 第一种方式所采用的基本函数 ?i ( x ∈ [ x0 , xL ]) 非 常复杂,而且是在全域上 [ x0 , xL ] 定义的,但它是高次连续函数,一般情况下,仅采用几个 基底函数就可以得到较高的逼近精度;而第二种方式所采用的基本函数

ai + bi x( x ∈ [ xi , xi +1 ]) 非常简单,而且是在子域上 [ xi , xi +1 ] 定义的,它通过各个子域组合出
全域 [ x0 , xL ] ,但它是线性函数,函数的连续性阶次较低,因此需要使用较多的分段才能得 到较好的逼近效果,则计算工作量较大。 对于第一种的函数逼近方式,就是力学分析中的经典瑞利-里兹方法(Rayleigh-Ritz principle)的思想, 而针对第二种的函数逼近方式, 就是现代力学分析中的有限元方法的思想, 其中的分段就是“单元”的概念。 基于分段的函数描述具有非常明显的优势:(1)可以将原函数的复杂性“化繁为简” , 使得描述和求解成为可能,(2)所采用的简单函数可以人工选取,因此,可取最简单的线性 函数,或取从低阶到高阶的多项式函数,(3)可以将原始的微分求解变为线性代数方程。但 分段的做法可能会带来的问题有:(1)因采用了“化繁为简” ,所采用简单函数的描述的能力 和效率都较低,(2)由于简单函数的描述能力较低,必然使用数量众多的分段来进行弥补, 因此带来较多的工作量。 综合分段函数描述的优势和问题,只要采用功能完善的软件以及能够进行高速处理的 计算机,就可以完全发挥“化繁为简”策略的优势,有限元分析的概念就在于此。

2.2 一维阶梯杆结构问题的求解
一维问题,即 1D(one dimension)问题,是最简单的分析对象,下面就以一个 1D 阶梯 杆结构为例,详细给出各种方法求解的过程,直观地引入有限元分析的基本思路,并以此逐 步介绍有限元分析的过程。

图 2-3

1D 阶梯杆结构

【典型例题】2.2(1) 1D 阶梯杆结构问题的材料力学求解 如图 2-3 所示为一个阶梯杆结构,已知相应的弹性模量和结构尺寸为

E1 = E2 = 2 × 107 Pa , A1 = 2 A2 = 2cm 2 , l1 = l2 = 10cm , F = 10N
用材料力学方法求解该问题。

9

有限元分析基础教程

曾攀

解答:首先对右端的杆件②进行力学分析,见图 2-4。

图 2-4 杆件分离体的受力分析

将两个杆件进行分解,并标出每一个关联节点处的受力状况,由于在 C 点处受有外力 F,则

Pc = F = 10 N
由杆件②的平衡关系可知,有
2 I B = Pc = F = 10 N 1 2 1 2

(2-3)

(2-4)

由于 I B 和 I B 是一对内力,也为作用力与反作用力,因此,有关系 I B ? I B = 0 ,则可计算出 所有作用力与内力的值为
1 PA = I B = I B = PC = F = 10 N
2

(2-5)

下面计算每根杆件的应力,这是一个等截面杆受拉伸的情况,则杆件①的应力 σ 1 为

σ1 =
杆件②的应力 σ 2 为

PA 10 N N = = 5 × 10 4 2 = 5 × 10 4 Pa 2 A1 2cm m

(2-6)

σ2 =

PC 10 N N = = 1×105 2 = 1× 105 Pa 2 A2 1cm m

(2-7)

由于材料是弹性的,由虎克定律(Hooke law)有

σ 1 = E1ε1 ? ? σ 1 = E2ε 2 ?
其中 ε1 和 ε 2 为杆件①和②的应变,则有

(2-8)

? 5 ×104 PA = 2.5 ×10?3 ? 7 E1 2 × 10 Pa ? ? 5 σ 2 1×10 PA ?3 ? = = 5 × 10 ε2 = ? E2 2 ×107 Pa ?

ε1 =

σ1

=

(2-9)

10

有限元分析基础教程

曾攀

由应变的定义可知,它为杆件的相对伸长量,即 ε = Δl / l ,因此, Δl = ε ? l ,具体对杆件 ①和②,有

Δl1 = ε1 ? l1 = 2.5 × 10?3 × 10 = 2.5 × 10?2 cm ? ? ? 2 ?3 Δl2 = ε 2 ? l2 = 5 × 10 × 10 = 5 × 10 cm ? ?

(2-10)

由于左端 A 为固定,则该点沿 x 方向的位移为零,记为 u A = 0 ,而 B 点的位移则为杆 件①的伸长量 Δl1 ,即 u B = Δl1 = 2.5 × 10 cm ,C 点的位移为杆件①和②的总伸长量,即
?2

uc = Δl1 + Δl2 = 7.5 × 10?2 cm ,则归纳为以上结果,有完整的解答:

σ 1 = 5 ×104 Pa, σ 2 = 1× 105 Pa

? ? ε1 = 2.5 ×10?3 , ε 2 = 5 × 10?3 ? ? u A = 0, uB = 2.5 × 10?2 cm, uC = 7.5 ×10?2 cm ?

(2-11)

以上的求解结果如图 2-5 所示。

图 2-5 阶梯杆结构的材料力学解答

讨论: 以上完全按照材料力学的方法,将对象进行分解来获得问题的解答,它所求解的基 本力学变量是力(或应力) ,由于以上问题非常简单,而且是静定问题,所以可以直接求出, 但对于静不定问题,则需要变形协调方程(compatibility equation),才能求解出应力变量,在 构建问题的变形协调方程时,则需要一定的技巧; 若采用位移作为首先求解的基本变量, 则可以使问题的求解变得更规范一些,下面就基于 A、B、C 三个点的位移 u A、u B、uC 来 进行以上问题的求解。 【典型例题】2.2(2) 1D 阶梯杆结构的节点位移求解及平衡关系 所处理的对象与【典型例题】2.2(1)相同,要求分别针对每个连接节点,基于节点的位 移来构建相应的平衡关系,然后再进行求解。 解答:考虑图 2-3 所示杆件的受力状况,分别画出每个节点的分离受力图,如图 2-6 所示。

11

有限元分析基础教程

曾攀

图 2-6 1D 阶梯杆结构的各种平衡关系

首先分析图 2-6(c)中杆①内部的受力及变形状况,它的绝对伸长量为 (u B ? u A ) ,则相 应伸长量 ε1 为

ε1 =
由虎克定律,它的应力 σ 1 为

uB ? u A l1

(2-12)

σ 1 = E1ε1 =
杆①的内力 I B 为
1 I B = σ 1 A1 =

E1 ( uB ? u A ) l1

(2-13)

1

E1 A1 ( uB ? u A ) l1
2

(2-14)

对于杆②进行同样的分析和计算,有它的内力 I B 为
2 I B = σ 2 A2 =

E2 A2 ( uC ? uB ) l2

(2-15)

由图 2-6(b)节点 A、B、C 的受力状况,分别建立它们各自的平衡关系如下。 对于节点 A,有平衡关系
1 ? PA + I B = 0

(2-16)

将式(2-14)代入,有

? PA +
对于节点 B,有平衡关系

E1 A1 ( uB ? u A ) = 0 l1

(2-17)

1 2 ?IB + IB = 0

(2-18)

12

有限元分析基础教程

曾攀

将式(2-14)和式(2-15)代入,有

?
对于节点 C,有平衡关系

E1 A1 E A ( uB ? u A ) + 2 2 ( uC ? uB ) = 0 l1 l2

(2-19)

2 PC ? I B = 0

(2-20)

将式(2-15)代入上式,有

PC ?

E2 A2 ( uC ? uB ) = 0 l2

(2-21)

将节点 A、B、C 的平衡关系写成一个方程组,有

? ?EA ? ?EA ? ? PA ? ? 1 1 ? u A + ? 1 1 ? uB + 0 = 0 ? ? l1 ? ? l1 ? ? ? ? E1 A1 ? ? E1 A1 E2 A2 ? ? E2 A2 ? ? + 0+? ? uA ? ? ? uB + ? ? uC = 0 ? l2 ? ? l1 ? ? l1 ? l2 ? ? ? ?E A ? ?E A ? ? PC ? 0 + ? 2 2 ? uB ? ? 2 2 ? uC = 0 ? l2 ? l2 ? ? ? ?
写成矩阵形式,有

(2-22)

? E1 A1 ? l ? ? PA ? ? 1 ? ? ? E1 A1 ? 0 ? ? ?? l 1 ? ? ? ? PC ? ? ? 0 ? ?

?

E1 A1 l1

E1 A1 E2 A2 + l1 l2 ? E2 A2 l2

? ? ? ?u A ? ?0 ? E2 A2 ? ? ? ? ? ? ? uB = 0 l2 ? ? ? ? ? ?uC ? ?0 ? ? ? ? ? E2 A2 ? ? l2 ? ? 0

(2-23)

将材料弹性模量和结构尺寸代入(2-23)方程中,有以下方程(采用国际单位)

? 4 × 104 ? 4 ? ?4 × 10 ? 0 ?

?4 × 104 6 ×104 ?2 × 10
4

0 ? ?u A ? ? ? PA ? ? ? ? ?2 × 104 ? ?uB ? = ? 0 ? ? ? 2 × 104 ? ?uC ? ? 10 ? ?? ? ? ?

(2-24)

由于左端点为固定,即 u A = 0 ,该方程的未知量为 u B、uC、PA ,求解该方程,有

u B = 2.5 ×10?4 m ? ? uC = 7.5 ×10?4 m ? ? PA = 10N ?

(2-25)

可以看出这里的 PA 就是支座反力(reaction force of support), 下面就很容易求解出杆①和②中 的其它力学量,即

13

有限元分析基础教程

曾攀

ε1 =

uB ? u A ? = 2.5 × 10?3 ? l1 ? uC ? uB ?3 ? ε2 = = 5 ×10 ? l2 ? σ 1 = E1ε1 = 5 ×104 Pa ? ? σ 2 = E2ε 2 = 1×105 Pa ? ?

(2-26)

这样得到的结果与【典型例题】2.2(1)所得到的结果完全一致。下面再对【典型例题】2.2(2) 的方法作进一步的讨论。 讨论:还可以将式(2-23)写成
( 3×1) ( 3×1) ( 3×1) ( 3×1)

P ? I =0

(2-27)

其中 P 称为外力列阵(load matrix), I 称为内力列阵(inner force matrix)或变形力列阵 (deformed force matrix),这里矩阵符号的下标表示行和列的维数,这两个列阵分别为

? ? PA ? ? ? P = 0 ? ( 3×1) ? ? PC ? ? ?

(2-28)

? E1 A1 ? ? l1 ? EA I = ?? 1 1 ( 3×1) ? l1 ? ? 0 ? ?

?

E1 A1 l1

E1 A1 E2 A2 + l1 l2 ? E2 A2 l2

? ? ? ?u A ? E2 A2 ? ? ? ? ? uB l2 ? ? ? ?uC ? ? ? E2 A2 ? ? l2 ? ? 0

(2-29)

式(2-27)的物理含义就是内力与外力的平衡关系, 由式(2-29)可知, 内力表现为各个节点上的 内力,并且可以通过节点位移(nodal displacement) u A、u B、uC 来获取。 由方程(2-23)可知,这是一个基于节点 A、B、C 描述的全结构的平衡方程,该方程的 特点为: (a)基本的力学参量为节点位移 u A、u B、uC 和节点力 PA、P 。 C (b)直接给出全结构的平衡方程,而不是象【典型例题】2.2(1)那样,需要针对每一个 杆件去进行递推。 (c)在获得节点位移变量 u A、u B、uC 后,其它力学参量(如应变和应力) ,都可以分别 求出(见式(2-26)) 。 为了将方程(2-23)写成更规范、更通用的形式,用来求解【典型例题】2.2(1)所示结构的 更一般的受力状况,下面在式(2-23)的基础上,直接推导出通用平衡方程。 【典型例题】2.2(3) 1D 阶梯杆结构基于位移求解的通用形式 由方程(2-23),可以发现它具有求解这类结构的通用形式,讨论它的物理含义,以总结 出相应的规律。

14

有限元分析基础教程

曾攀

解答:将式(2-23)写成

? E1 A1 ? ? l1 ? E1 A1 ?? ? l1 ? ? 0 ? ?

?

E1 A1 l1

E1 A1 E2 A2 + l1 l2 ? E2 A2 l2

? ? ? ?u A ? ? PA ? E2 A2 ? ? ? ? ? ? ? uB = PB l2 ? ? ? ? ? ?uC ? ? PC ? ? ? ? ? E2 A2 ? ? l2 ? ? 0

(2-30)

再将其分解为两个杆件之和,即写成

? E1 A1 ? l ? 1 ? E1 A1 ?? l 1 ? ? 0 ? ?

?

E1 A1 l1

E1 A1 l1 0

? ? 0? ? 0 0 ? ?u A ? ? ? E2 A2 ? 0 ? ?u B ? + ?0 ? ? l2 ? ?uC ? ? ? ? ? E2 A2 0? ?0 ? ? l2 ? ? ?

? ? 0 ? ?u A ? ? PA ? E2 A2 ? ? ? ? ? ? ? uB = PB l2 ? ? ? ? ? ?uC ? ? PC ? ? ? ? ? E2 A2 ? ? l2 ? ?

(2-31)

以上式(2-31)左端的第 1 项实质为

? E1 A1 ? l ? 1 ? E1 A1 ?? l ? 1
1 1

?

E1 A1 ? 1 l1 ? ?u A ? E1 A1 ?u A ? uB ? ? ? I B ? ?? ? = =? 1 ? ? ? E1 A1 ? ?uB ? l1 ?uB ? u A ? ? I B ? l1 ? ?

(2-32)

上式中的 ? I B 及 I B 见图 2-6 (c) 含义为杆件①中的左节点的内力和右节点的内力。 , 同样地, 式(2-31)左端的第 2 项实质为

? E2 A2 ? l ? 2 ? E2 A2 ?? l 2 ?
2 2

?

E2 A2 ? l2 ? ?uB ? E2 A2 ? = E2 A2 ? ?uC ? l2 ? ? ? l2 ?

2 ?uB ? uC ? ? ? I B ? =? 2 ? ?u ? u ? B? ? C ? IB ?

(2-33)

上式中的 ? I B 及 I B 见图 2-6(c),含义为杆件②中的左节点的内力和右节点的内力。 可以看出:方程(2-31)的左端就是杆件①的内力表达和杆件②的内力表达之和,这样就 将原来的基于节点的平衡关系, 变为通过每一个杆件的平衡关系来进行叠加。 这里就自然引 入单元的概念, 即将原整体结构进行“分段”, 以划分出较小的“构件”(component), 每一个“构 件”上具有节点,还可以基于节点位移写出该“构件”的内力表达关系,这样的“构件”就叫做 单元(element),它意味着在几何形状上、节点描述上都有一定普遍性(generalization)和标准 性(standardization),只要根据实际情况将单元表达式中的参数(如材料常数、几何参数)作 相应的代换,它就可以广泛应用于这一类构件(单元)的描述。 从式(2-32)和式(2-33)可以看出, 虽然它们分别用来描述杆件①和杆件②的, 但它们的表 达形式完全相同,因此本质上是一样,实际上,它们都是杆单元(bar element)。 可以将杆单元表达为如图 2-7 所示的标准形式。

15

有限元分析基础教程

曾攀

图 2-7 杆单元的表达

将单元节点位移写成

?u ? T q e = ? 1 ? = [u1 u2 ] ( 2×1) ?u2 ?
将单元节点外力写成

(2-34)

?P ? Pe = ? 1 ? = [ P 1 2×1) ( ? P2 ?
由式(2-32),该单元节点内力为

P2 ]

T

(2-35)

? EA ? ? I1 ? ? l ? I ? = ? EA ? 2 ? ?? ? l ?
e

?

EA ? l ? ? u1 ? ? EA ? ?u2 ? ? ? ? l ?

(2-36)

它将与单元的节点外力 P 相平衡,则有 P = ? I1 , P2 = I 2 因此,该方程可以写成 1
( 2×1)

? EA ? l ? ? ? EA ? l ?
进一步表达成

?

EA ? l ? ? u1 ? ? P ? 1 ?? ? = ? ? EA ? ?u2 ? ? P2 ? l ? ?

(2-37)

( 2×2 ) ( 2×1)

K e qe = Pe

( 2×1)

(2-38)

其中

EA ? ? EA ? ? l l ? ? K11 Ke = ? ?=? ( 2×2 ) ? ? EA EA ? ? K 21 ? l l ? ? ?

K12 ? K 22 ? ?

(2-39)

可以看出,方程(2-38)是单元内力与外力的平衡方程(equilibrium equation),它与单元的刚度 方 程 (stiffness equation) 是 相 同 的 ( 见 式 (3-41)) 。 K 叫 做 单 元 的 刚 度 矩 阵 (stiffness matrix), K11、K12、K 21、K 22 叫做刚度矩阵中的刚度系数(stiffness coefficient)。
e

16

有限元分析基础教程

曾攀

2.3 有限元分析的基本流程
下面以一个 1D 三连杆结构为例,展现有限元分析的全部过程。 【典型例题】2.3(1) 1D 三连杆结构的有限元分析过程 采用杆单元的方法,求解如图 2-8 所示结构的所有力学参量。相关的材料参量和尺寸为

E1 = E2 = E3 = 2 × 105 Pa , 3A1 = 2 A2 = A3 = 0.06m 2 l1 = l2 = l3 = 0.1m

图 2-8

1D 三连杆结构的受力状况

解答: 所谓基于单元的分析方法, 就是将原整体结构按几何形状的变化性质划分节点并进行 编号,然后将其分解为一个个小的构件(即:单元) ,基于节点位移,建立每一个单元的节 点平衡关系(叫做单元刚度方程) ,对于杆单元来说就是式(2-38);下一步就是将各个单元进 行组合和集成,类似于式(2-31),以得到该结构的整体平衡方程(也叫做整体刚度方程) ,按 实际情况对方程中一些节点位移和节点力给定相应的值(叫做处理边界条件) ,就可以求解 出所有的节点位移和支反力, 最后在得到所有的节点位移后, 就可以计算每一个单元的其它 力学参量(如应变、应力) ;下面给出该问题的有限元分析过程。 (1)节点编号和单元划分 由于该结构由三根不同几何尺寸的杆件组成,并且在一些杆件连接处还作用有集中载 荷,因此,需要在杆件连接处划分出节点,这样对于该结构就自动给出三个单元,其节点及 单元编号(numbering of element)见图 2-8,将每一个单元分离出来,并标出每一个节点的位 移和外力,如图 2-9 所示,注意,这里位移和力的方向都以 x 正方向来标注。

图 2-9

各个单元的节点位移和外力

(2)计算各单元的单元刚度方程 可以看出,图 2-9 所示的每一个单元都和图 2-7 的单元类似,则所对应的刚度方程也应 与式(2-37)类似,只需要将其中的各个参数进行代换,下面直接给出对应于单元①、②、 ③的单元刚度方程。 单元①的刚度方程为

17

有限元分析基础教程

曾攀

? E1 A1 ? l ? 1 ? E1 A1 ?? l 1 ?
单元②的刚度方程为

?

E1 A1 ? l1 ? ? u1 ? ? P (1) ? ? =? 1 ? E1 A1 ? ?u2 ? ? P2 (1) ? ? ? ? l1 ? E2 A2 ? l2 ? ?u2 ? ? P3(2) ? ? =? ? E2 A2 ? ? u3 ? ? P3(2) ? ? ? l2 ? ? E3 A3 ? l2 ? ?u3 ? ? P3(3) ? ? =? ? E3 A3 ? ?u3 ? ? P4 (3) ? ? ? l3 ? ?

(2-40)

? E2 A2 ? l ? 2 ? E2 A2 ?? l 2 ?
单元③的刚度方程为

?

(2-41)

? E3 A3 ? l ? 3 ? E3 A3 ?? l 3 ?

?

(2-42)

(3)组装各单元刚度方程 由于整体结构是由各个单元按一定连接关系组合而成的, 因此, 需要按照节点的对应位 置将以上方程(2-40)(2-41)(2-42)进行组装,以形成一个整体刚度方程,即 、 、

u1 ↓ ? E1 A1 ? l ? 1 ? E1 A1 ?? l 1 ? ? ? 0 ? ? ? 0 ? ? ?

u2 ↓ E1 A1 l1 ?

u3 ↓ 0 E2 A2 l2

u4 ↓ ? ? ? ? ? u1 ? ? P (1) ? 0 ? ? ? ? (1) 1 (2) ? ? ?u2 ? = ? P2 + P2 ? E A ? ? u ? ? P (2) + P (3) ? ? 3 3 ? ? 3 ? ? 3 (3) 3 ? l3 ? ?u4 ? ? P4 ? ? ? ? ? ? E3 A3 ? l3 ? ? 0

E1 A1 E2 A2 + l1 l2 ? E2 A2 l2 0

E2 A2 E3 A3 + l2 l3 ? E3 A3 l3

(2-43)

为表达更清楚,在式(2-43)的刚度矩阵上方标明了所对应的变量。上面的组装过程,实际 上就是将各个单元方程按照节点编号的位置进行集成。 式(2-43)中的 P 、P2 1
(1) (2)

+ P2(2)、P3(2) + P3(3)、P4(3) 就是节点 1、2、3、4 上的合成节点

力,即

P = P (1),P2 = P2(1) + P2(2),P3 = P3(2) + P3(3),P4 = P4(3) 1 1

(2-44)

对照图 2-8 可知, P = ?100 N,P2 = 0,P = 50 N ,而 P4 为支座的支反力。将该结构的材 1 3 料参数和几何尺寸参数代入式(2-43)中,则有

18

有限元分析基础教程

曾攀

? 4 × 104 ? 4 ? ?4 × 10 ? 0 ? ? 0 ?

?4 × 104 1× 105 ?6 × 104 0

0 ?6 × 104 1.8 × 105 ?1.2 × 105

? ? u1 ? ? P ? 1 ?? ? ? ? 0 ? ?u2 ? = ? P2 ? 5? ?1.2 × 10 ? u3 ? ? P3 ? ?? ? ? ? 1.2 ×105 ? ?u4 ? ? P4 ? ?? ? ? ? 0

(2-45)

上式中的 u1、u2、u3、u4 为节点位移, P、P2、P 、P4 为节点力,可以看出,图 2-9 中分别 1 3 就单元①、②、③写出了各自的节点力,如对于节点 2,即写出了单元①中节点 2 的节点力

P2(1) ,又给出了单元②中节点 2 的节点力 P2(2) ,可以看出,在单元组装后,实际上只需要合
成后的节点力(见 P、P2、P 、P4 ) ;因此,今后只需要对各个单元的刚度系数按对应节点 1 3 位移的位置进行组装,而节点力只需直接写出即可,见方程式(2-45) 。 (4)处理边界条件并求解 图 2-8 所示结构的位移边界条件为: u4 = 0 ,将已知的节点位移和节点力代入后,则 方程(2-45)变为

? 4 × 104 ? 4 ? ?4 × 10 ? 0 ? ? 0 ?

?4 × 104 1× 105 ?6 × 104 0

0 ?6 ×104 1.8 × 105 ?1.2 × 105

? ? u1 ? ? ?100 N ? ?? ? ? ? 0 ? ? u2 ? = ? 0 ? ?1.2 × 105 ? ? u3 ? ? 50 N ? ?? ? ? ? 1.2 ×105 ? ?u4 = 0 ? ? P4 ? ? ? ? ?? 0

(2-46)

由于 u4 = 0 ,则划掉上述刚度矩阵的第 4 列和第 4 行,则有

? 4 × 104 ? 4 ? ?4 ×10 ? 0 ?
求解上述方程,有

?4 ×104 1×105 ?6 ×10
4

0 ? ? u1 ? ? ?100 ? ? ?6 × 104 ? ?u2 ? = ? 0 ? ? ? ? ? 1.8 × 105 ? ? u3 ? ? 50 ? ? ?? ? ?
? ? ? ? ?

(2-47)

u1 = ?4.58333 × 10?3 m u2 = ?2.08333 ×10?3 m u3 = ?4.16667 ×10?4 m
(5)求支反力

(2-48)

在求得所有节点位移后,见式(2-48) ,并且 u4 = 0 ,由方程组(2-46)的最后一行方 程,可求出支反力 P4 为

P4 = ?1.2 × 105 × u3 = 50N
(6)求各个单元的其它力学量(应变、应力)

(2-49)

19

有限元分析基础教程

曾攀

由式(2-12) ,可以求出各个单元的应变,即

ε1(1) =

u2 ? u1 ( ?2.08333 + 4.5833) × 10 = l1 0.1

?3

= 2.49997 × 10 ?2
?3

(2-50)

ε

(2)

u ?u ( ?0.416667 + 2.08333) ×10 = 3 2 = l2 0.1

= 1.6667 ×10?2

(2-51)

ε (3) =

u4 ? u3 0 + 4.16667 ×10?4 = = 4.16667 ×10?3 l3 0.1

(2-52)

由式(2-13) ,可以求出各个单元的应力,即

σ (1) = E1ε (1) = 4.999 × 103 Pa σ (2) = E2ε (2) = 3.3333 × 103 Pa σ (3) = E3ε (3) = 8.3333 × 10 2 Pa

(2-53) (2-54) (2-55)

这样可以得到一种直观的有限元分析思路, 就是将复杂的几何和受力对象划分为一个一 个形状比较简单的标准“构件”,称为单元,然后给出单元节点的位移和受力描述,构建起单 元的刚度方程, 再通过单元与单元之间的节点连接关系进行单元的组装, 可以得到结构的整 体刚度方程,进而根据位移约束和受力状态,处理边界条件,并进行求解,基本流程的示意 见图 2-10。

图 2-10 有限元分析的基本流程图示

2.4 有限元分析的特点
有限元分析的最大特点就是标准化和规范化,这种特点使得大规模分析和计算成为可 能, 当采用了现代化的计算机以及所编制的软件作为实现平台时, 则复杂工程问题的大规模 分析就变为了现实。 实现有限元分析标准化和规范化的载体就是单元, 这就需要我们构建起各种各样的具有 代表性的单元,一旦有了这些单元,就好像建筑施工中有了一些标准的预制构件(如梁、楼 板等),可以按设计要求搭建出各种各样的复杂结构,如图 2-11 所示。

20

有限元分析基础教程

曾攀

图 2-11 在建筑中采用一些基本构件可以搭建出各种各样的复杂结构

图 2-11 中所示的基本构件,实际上就是有限元分析中的“单元” ,而搭建出的复杂结 构就是我们需要分析的“对象” ,下面用图示的方法列出有限元分析中常用的一些典型单元 (见图 2-12)。

图 2-12 常用的一些典型单元(ANSYS 平台中)

有限元分析的最主要内容,就是研究单元,即首先给出单元的节点位移和节点力,然后 基于单元节点位移与节点力的相互关系(见【典型例题】2.2(2))可以直接获得相应的刚度 系数,进而得到单元的刚度方程,实际上就是要得到针对单元节点的平衡方程,这就是单元

21

有限元分析基础教程

曾攀

的刚度方程,就可以针对实际的复杂结构,根据实际的连接关系,将单元组装为整体刚度方 程,这实际上也是得到整体结构的基于节点位移的整体平衡方程。 因此,有限元方法的主要任务就是对常用的各种单元(包括 1D、2D、3D 问题的单元) 构造出相应的单元刚度矩阵[13][14];当然,如果还采用如【典型例题】2.2(2)所示的直接法来 进行构造,会非常烦琐,而采用能量原理(如:虚功原理或最小势能原理)来建立相应的平 衡关系则比较简单[15],这种方法可以针对任何类型的单元进行构建,以得到相应的刚度矩 阵,这种推导单元刚度矩阵的方法的力学基础见第 3 章。

2.5 本章要点
函数的全域逼近与子域分段逼近的概念 将对象离散为一些标准的单元 将单元的节点位移作为基本力学变量 建立单元的节点的平衡关系 单元刚度矩阵以及直接求取的方法 有限元分析的标准流程(离散化、单元描述、整体组装、问题求解)

22

有限元分析基础教程

曾攀

第 3 章 杆梁结构的有限元分析方法
下面先从简单的杆梁结构入手全面介绍有限元方法,接着在后面的章节就连续体问题进行 讨论。

3.1 杆梁结构分析的工程概念
在建筑结构中,杆、梁、板是主要的承力构件,关于它们的计算分析对于建筑结构设计来 说具有非常重要的作用,对杆、梁、板的建模将充分考虑到实际结构的几何特征及连接方式, 并需要对其进行不同层次的简化,可以就某一特定分析目的得到相应的 1D、2D、3D 模型,由 于在设计时并不知道结构的真实力学性能(或许还没有实验结果,或许还得不到精确的解析解), 仅有计算分析的一些结果,因此,一种进行计算结果校核或验证的可能方法,就是对所分析对 象分别建立 1D、2D、3D 模型,来进行它们之间的相互验证和核对;图 3-1 给出一个建筑结构 中的杆梁框架以及建模简化过程。

图 3-1 建筑结构中的杆梁框架以及建模简化过程

23

有限元分析基础教程

曾攀

3.2 杆件有限元分析的标准化表征与算例
3.2.1 杆件分析的基本力学原理
杆件是最常用的承力构件,它的特点是连接它的两端一般都是铰接接头,因此,它主要是 承受沿轴线的轴向力,因两个连接的构件在铰接接头处可以转动,则它不传递和承受弯矩。 有一个左端固定的拉杆,其右端承受一外力 P。该拉杆的长度为 l,横截面积为 A,弹性模 量为 E,如图 3-2 所示,这是一个一维问题,下面讨论该问题的力学描述与求解。

E,A

P

x

l
图 3-2 一端固定的拉杆

【基本变量】3.2.1(1) 1D 问题的基本变量 由于该问题是沿 x 方向的一维问题,因此只有沿 x 方向的基本变量,即 定义沿 x 方向移动为位移:

u (x ) 定义沿 x 方向的相对伸长(或缩短)量为应变: ε x (x) 定义沿 x 方向的单位横截面上的受力为应力: σ x (x )

【基本方程】3.2.1(2) 1D 问题的基本方程 该问题的三大类基本方程和边界条件如下。 ① 取出杆件的任意一个截面,可得到平衡方程(无体力)为

σ x = c1 或
其中 c1 为待定的常数。 ②

dσ x =0 dx

(3-1)

取出杆件 x 位置处的一段长度 dx,设它的伸长为 du,则它的相对伸长量为

εx =


du dx

(3-2)

该方程称为几何方程。 由该材料的拉伸试验,可得到该材料的虎克定律为

εx =
该方程也称为物理方程。 边界条件(BC, boundary condition) 位移边界条件 BC(u)

σx
E

(3-3)



u ( x) x =0 = 0

(3-4)

24

有限元分析基础教程

曾攀

力边界条件 BC(p)

σ x ( x ) x =l =

P = px A

(3-5)

对于以上的力边界条件,只能作为一种近似,因为在 x = l 的端面, σ x (x) 不应是均匀分布的。由 圣维南原理(Saint-Venant’s principle),在远离 x = l 的截面,力的边界条件才较好地满足。 从求解思路来说,可以有两类方法来对该问题进行求解,即 (1)直接求解方法:该问题比较简单,因此,可以由 3 个方程来直接求解 3 个变量。 (2)基于试函数的间接方法:可以先选取一个变量(如位移)作为最基本的待求变量,将其它变量 都用它来表达,并采用间接的近似求解方法;具体的做法为:先对待求的位移变量假设一种事 先满足位移边界条件的可能解(其中有一些待定的系数),称为试函数(trail function),让该受力系 也可以让 统的势能取最小值来最后确定出可能解(试函数)中的那些待定系数(unknown constant); 该受力系统的内部变形虚功等于外部施加力的虚功,来求出试函数中的那些待定系数。 下面针对图 3-2 所示的一端固定的拉杆问题,分别讨论基于直接求解方法以及基于试函数 的间接方法的求解过程。 【求解原理】3.2.1(3) 1D 问题的直接求解 就所针对的问题,对方程(3-1)-(3-3)进行直接求解,可得到以下结果。

σ x ( x) = c1 ε x ( x) =
c1 E c1 u ( x) = x + c2 E

? ? ? ? ? ? ? ? ?

(3-6)

其中 c1 及 c2 为待定常数,由边界条件(3-4)和(3-5),可求出(3-6)中的常数为 c2 = 0 , c1 = 此,有最后的结果

P 。因 A

σ x ( x) =

P A P ε x ( x) = EA P u( x ) = x EA

? ? ? ? ? ? ? ? ?

(3-7)

下面讨论基于试函数的间接方法的求解过程,先介绍相应的求解原理,然后就图 3-2 所示 的一端固定的拉杆问题给出解答。 【求解原理】3.2.1(4) 1D 问题的虚功原理求解 先以一个简单的结构静力平衡问题来描述虚功原理的基本思想, 然后再具体求解一端固 定的拉杆问题。

25

有限元分析基础教程

曾攀

图 3-3 平衡力系由于微小扰动产生相应的虚位移

如图 3-3 所示的一个平衡力系,由于该系统处于平衡状态,则有

PA l B = PB l A

(3-8)

假想在该平衡力系上作用有微小的扰动(不影响原平衡条件),且外力所作用的位置产生了微小 的位移变化,即 Δ A , Δ B 。该假想的位移如果不影响原平衡条件,应满足以下几何关系

Δ B lB = Δ A lA

(3-9)

这就是任意扰动的位移应满足的条件,称为许可位移条件,我们把满足许可位移条件的、任意 微小的假想位移称为虚位移(假想的)(virtual displacement)。 将式(3-9)代入式(3-8)中,有

PA Δ A ? PB Δ B = 0

(3-10)

PA Δ A 为 A 点的虚功(virtual work)(即外力 PA 在 A 点虚位移 Δ A 上所做的功),? PB Δ B 为 B 点的
虚功(即外力 PB 在 B 点虚位移 Δ B 上做的功,注意负号表示力与位移的方向相反)。 实质上,关系式(3-10)表达了基于虚位移的虚功原理(principle of virtual work),即:对于一 个处于平衡状态的系统,作用于系统上的所有外力在满足许可位移条件的虚位移上所做的虚功 总和恒为零。 现在进一步讨论弹性力学中有关变形体的虚功原理,这时的虚功应包括外力虚功 δW 和内 力虚功 ? δU , δ U 叫做虚应变能(virtual strain energy)。由于弹性体在变形过程中,内力是抵抗 变形所产生的,其方向总是与变形的方向相反,所以内力虚功取负。由于虚功总和为零,则有

δW ? δU = 0
所以

(3-11)

(3-12) 弹性力学中的虚功原理可表述为:在外力作用下处于平衡状态的变形体,当给物体以微小 虚位移时,外力所做的总虚功等于物体的总虚应变能(即应力在由虚位移所产生虚应变上所作的 功)。注意这里的虚位移是指仅满足位移边界条件 BC(u)的许可位移。 下面应用虚功应力来具体求解如图 3-2 所示的一端固定的拉杆问题,设有满足位移边界条

δU = δW

26

有限元分析基础教程

曾攀

件的位移场

u ( x ) = cx

(3-13)

可以验证:它满足位移边界条件。这是一个待定函数,也称为试函数,所谓该函数是待定的, 就是因为它中间有一个待定系数 c ,这就需要通过一个原理来确认它,下面由虚功原理来进行 确认。基于式(3-13)的试函数,则它的应变、虚位移以及虚应变为

ε ( x) = c δ u ( x) = δ c ? x δε ( x) = δ c

? ? ? ? ?

(3-14)

其中 δ c 为待定系数的增量。计算如图 3-2 所示算例的虚应变能以及外力虚功为

δ U = ∫ σ xδε x dΩ = ∫
Ω

l

0



A

E ? ε xδε x dA ? dx = E ? c ? δ c ? A ? l

(3-15) (3-16)

δW = P ?δ u( x = l) = P ?δ c ? l
由虚功原理(3-12),有

E ? c ?δ c ? A ? l = P ?δ c ? l
消去 δ c 后,有解

(3-17)

c=
代回式(3-13)中,就可以得到该问题的解。

P EA

(3-18)

【求解原理】3.2.1(5)

1D 问题的最小势能原理求解

先介绍最小势能原理的基本表达式。设有满足位移边界条件 BC(u)的许可位移场 u ( x ) , 计算该系统的势能(potential energy)为

Π (u ) = U ? W
其中 U 为应变能, W 为外力功,对于如图 3-2 所示的算例,有

(3-19)

1 σ x (u ( x)) ? ε x (u ( x))dΩ 2 ∫Ω W = P ? u( x = l ) U=

? ? ? ? ?

(3-20)

对于包含有待定系数的试函数 u ( x ) 而言, 真实的位移函数 u ( x ) 应使得该系统的势能取极小值, 即
u ( x )∈BC ( u )

min

[Π (u ) = U ? W ]

(3-21)

27

有限元分析基础教程

曾攀

下面应用最小势能原理来具体求解如图 3-2 所示的一端固定的拉杆问题,如同样取满足位 移边界条件的位移场(3-13),则计算应力、应变为

du =c dx σ x ( x) = E ? ε x ( x) = E ? c

ε x ( x) =

? ? ? ? ?

(3-22)

由式(3-20)则计算该系统的势能为

Π (u ) = U ? W =
由式(3-21),应对式(3-23)求极值,即

1 E ? c2 ? A ? l ? P ? c ? l 2

(3-23)

?Π (u ) =0 ?c
则可以求出 c = P / EA ,与式(3-18)的结果相同。

(3-24)

由上面的计算可以看出,基于试函数的方法,包括虚功原理以及最小势能原理(principle of minimum potential energy),仅计算系统的能量,实际上就是计算积分(见式(3-15)以及(3-20)),然 后转化为求解线性方程,不需求解微分方程,这样就大大地降低了求解难度。同时,也可以看 出,试函数的方法的关键就是如何构造出适合于所求问题的位移试函数,并且该构造方法还应 具有规范性以及标准化,基于“单元”的构造方法就可以完全满足这些要求。

3.2.2 局部坐标系中的杆单元描述
最简单的标准单元就是杆单元,下面就要研究它的试函数描述,以及计算它的应变能和外 力功。 【单元构造】3.2.2(1) 杆单元的描述

单元的描述包括单元的几何及节点描述、位移场、应变场、应力场、势能,也就是要充 分利用描述问题的三大类变量以及三大类方程来计算单元的势能,然后,由最小势能原理(或虚 功原理)来得到单元的方程。实际上,单元内位移场的描述就是它的试函数的选取。 (1) 单元的几何及节点描述 图 3-4 所示为一个在局部坐标系(local coordinate system)中的杆单元(bar element),由于有两 个端节点(Node 1 和 Node 2) ,则基本变量为节点位移(向量)列阵 q
e

q e = [u1 u2 ]

T

(3-25)

将每一个描述物体位置状态的独立变量叫做一个自由度(DOF, degree of freedom),显然,以上的 节点位移为两个自由度。节点力(向量)列阵 P 为
e

28

有限元分析基础教程

曾攀

P e = [P 1

P2 ]

T

(3-26)

图 3-4 局部坐标系中的单元

若该单元承受有沿轴向的分布外载,可以将其等效到节点上,即表示为如式(3-26)所示的节 点力。利用函数插值、几何方程、物理方程以及势能计算公式,可以将单元的所有力学参量(即 场变量: u ( x ) , ε ( x ) , σ ( x ) 和 Π )用节点位移列阵 q 及相关的插值函数来表示。
e
e

(2) 单元位移场的表达 设该单元的位移场为 u ( x ) ,由 Taylor 级数,它可以表示为

u ( x) = a0 + a1 x + a2 x 2 +

(3-27)

该函数将由两个端节点的位移 u1 和 u2 来进行插值确定,可取式(3-27)的前两项来作为该单元的 位移插值模式(interpolation model):

u ( x) = a0 + a1 x
其中 a0 和 a1 为待定系数(unknowns)。 单元节点条件为

(3-28)

u ( x) |x =0 = u1 ? ? u ( x ) | x = l e = u2 ?
将节点条件(3-29)代入式(3-28),可以求得 a0 和 a1 为

(3-29)

a0 = u1 a1 = u2 ? u1 le

? ? ? ? ?

(3-30)

将其代入式(3-28)有

u ( x) = u1 + (

u2 ? u1 x x ) x = (1 ? e )u1 + ( e )u2 = N( x) ? q e e l l l

(3-31)

其中 N 叫做形状函数矩阵(shape function matrix),为

29

有限元分析基础教程

曾攀

x ? N ( x) = ?(1 ? e ) l ?

x? le ? ?

(3-32)

q e 叫做节点位移列阵(nodal displacement vector),即 q e = [u1 u2 ]T
(3) 单元应变场的表达 由弹性力学中的几何方程,有 1D 问题的应变

(3-33)

ε ( x) =
其中

du ( x) 1 = [? e dx l

1 ? u1 ? ] = B( x ) ? q e l e ?u 2 ? ? ?
1 ] le

(3-34)

B( x) =
叫做几何矩阵(strain-displacement matrix)。 (4) 单元应力场的表达

d 1 N ( x) = [? e dx l

(3-35)

由弹性力学中的物理方程,有 1D 问题的应力

σ ( x ) = E eε ( x ) = E e ? B ( x ) ? q e = S ( x ) ? q e
其中

(3-36)

S( x ) = E e ? B( x ) = [ ?
叫做应力矩阵(stress-displacement matrix)。 (5) 单元势能的表达 基于式(3-34)和式(3-36),有单元势能的表达式

Ee le

Ee ] le

(3-37)

Πe = U e ?W e = =

1 e e ∫Ωe σ ( x) ? ε ( x) ? dΩ ? ( P1 ? u1 + P2 ? u2 ) 2

1 le eT T q ? S ( x) ? B( x) ? q e ? Ae ? dx ? ( P e ? u1 + P2e ? u2 ) 1 2 ∫0
1 [u1 u2 ] 2 ? E e Ae ? e ? le e ? E A ?? l e ? E e Ae ? ? e ? l ? ? u1 ? ? e 1 ? ? ? ?P E e Ae ? ?u2 ? le ? ? ?u ? P2e ? ? 1 ? ? u ? 2?

=

T = 1 q eT ? K e ? q e ? Pe ? q e 2

(3-38)

30

有限元分析基础教程

曾攀

其中 K 叫做单元刚度矩阵(stiffness matrix of element),即
e

Ke =

Ee Ae ? 1 ?1? l e ??1 1 ? ? ?

(3-39)

P e 叫做节点力列阵(nodal force vector),即
? Pe ? P e = ? 1e ? ? P2 ?
(6) 单元的刚度方程 对式(3-38)取极小值,可以得到单元的刚度方程(stiffness equation of element)为
(2×2) (2×1)

(3-40)

K e ? qe = Pe

(2×1)

(3-41)

【典型例题】3.2.2(2)

变截面杆单元的推导

如图 3-5 所示,有一受轴载荷的线性变截面杆件,两端的截面积为 A1 和 A2,长度为 l,材 料的弹性模量为 E,试建立描述该杆件的一个杆单元。

(a) 线性变截面杆件 图 3-5 变截面杆件及单元

(b) 变截面杆单元

解答:建立的变截杆单元如图 3-5(b)所示,其中横截面积为

? x? ? x? A ( x ) = ?1 ? ? A1 + ? ? A2 ? l? ?l?
其节点位移列阵为 q e = [u1

(3-42)
T

u2 ] ,节点力列阵为 P e = [ P 1
T

P2 ]

由于为两个节点,则单元的位移模式取为

u ( x ) = a0 + a1 x = N ( x ) ? q e
由节点条件,可得到该单元的形状函数矩阵

(3-43)

31

有限元分析基础教程

曾攀

?? x ? ? x ? ? N ( x ) = ?? 1 ? ? ? ? ? ?? l ? ? l ? ?
由几何方程,可得到该单元的几何矩阵

(3-44)

? 1 1? B ( x ) = ?? ? l l? ?
则该单元的刚度方程为
( 2×2 ) ( 2×1)

(3-45)

K e ? qe = P

e

( 2×1)

(3-46)

其中刚度矩阵为
( 2×2 )

K e = ∫ BT ? E ? B ? d Ω
Ω

? 1? ?? l ? ? 1 1? = ∫ ? ? ? E ? ?? ? A ( x ) ? dx 0 ? l l? ? ?1? ? l ? ? ? l E ? 1 ?1? ?? x ? ? x? ? =∫ 2 ? ? ? ??1 ? l ? A1 + ? l ? A2 ? dx 0 l ? ? ? ? ? ?1 1 ? ??
l

=

E ( A1 + A2 ) ? 1 ?1? ? ?1 1 ? 2l ? ?

(3-47)
可以看出,与等截面杆单元相比,该线性变截面杆单元的截面取了平均面积。

3.2.3 杆单元的坐标变换
1. 平面杆单元的坐标变换 在工程实际中,杆单元可能处于整体坐标系(global coordinate system)中的任意一个位置, 如图 3-6 所示,这需要将原来在局部坐标系(local coordinate system)中所得到的单元表达等价地 变换到整体坐标系中, 这样, 不同位置的单元才有公共的坐标基准, 以便对各个单元进行集成(即 组装)。图 3-6 中的整体坐标系为( x o y ),杆单元的局部坐标系为( ox )。

32

有限元分析基础教程

曾攀

图 3-6 平面杆单元的坐标变换

局部坐标系中的节点位移为

q e = [u1 u2 ]
整体坐标系中的节点位移为

T

(3-48)

?_ _ _ q = ?u1 v1 u 2 ?
e
?

? v2 ? ?
_
?

T

(3-49)

如图 3-6 所示,在节点 1,整体坐标系下的节点位移 u 1 和 v1 ,其合成的结果应完全等效于局部 坐标系中的 u1;在节点 2,节点位移 u 2 和 v 2 合成的结果应完全等效于局部坐标系中的 u2,即存 在以下的等价变换关系
? ?

u1=u1 cos α + v1 sin α u2=u 2 cos α + v 2 sin α
写成矩阵形式

? ? ? ? ?

(3-50)

? u ? ?cosα qe = ? 1 ? = ? ? u2 ? ? 0

sin α 0

0 cosα

? u1 ? ? ? 0 ? ? v1 ? ? =Te ? q e sin α ? ?u 2 ? ?? ? ?v2 ? ? ?

(3-51)

其中 T 为坐标变换矩阵(transformation matrix),即
e

?cosα Te= ? ? 0

sin α 0

0 cosα

0 ? sin α ? ?

(3-52)

下面推导整体坐标系下的刚度方程;由于单元的势能是一个标量(能量),不会因坐标系的不同 而改变,因此,可将节点位移的坐标变换关系(3-51)代入原来基于局部坐标系的势能表达式中,

33

有限元分析基础教程

曾攀



1 Π e = q e T ? K e ? q e ? p eT ? q e 2 1 = q e T (Te T K e Te )q e ? (Te T ? P e )T ? q e 2 eT e e eT e = 1q K q ?P q 2
其中 K 为整体坐标系下的单元刚度矩阵, P 为整体坐标系下的节点力列阵,即
e e

(3-53)

K = Te T K e Te
e

(3-54) (3-55)
e

P = Te T P e
e

由最小势能原理(针对该单元),将式(3-53)对待定的节点位移列阵 q 取极小值,可得到整体坐标 系中的刚度方程
(4× 4) (4× 2) e e e K ?q = P (4×2)

(3-56)

对于如图 3-6 所示的杆单元,由式(3-54)具体给出

? cos 2 α ? E e Ae ? cos α sin α e K = e (4×4) l ? ? cos2 α ? ? ? cos α sin α ?
2. 空间杆单元的坐标变换

cos α sin α sin 2 α ? cos α sin α ? sin 2 α

? cos 2 α ? cos α sin α cos 2 α cos α sin α

? cos α sin α ? ? ? sin 2 α ? cos α sin α ? ? sin 2 α ? ?

(3-57)

空间问题中的杆单元如图 3-7 所示。

图 3-7 空间杆单元的坐标变换

该杆单元在局部坐标系下 ( Ox ) 的节点位移还是

34

有限元分析基础教程

曾攀

q e = [u1
而整体坐标系中 Ox y z 的节点位移列阵为

u2 ]

T

(3-58)

(

)

q e = ?u1 v1 ?
杆单元轴线在整体坐标系中的方向余弦为

w1 u 2

v2

w2 ? ?

T

(3-59)

cos(x, x) =

z 2 ? z1 y ? y1 x 2 ? x1 , cos( x, y ) = 2 , cos( x, z ) = l l l

(3-60)

其中 ( x1 , y1 , z1 ) 和 ( x 2 , y 2 , z 2 ) 分别是节点 1 和节点 2 在整体坐标系中的位置,是杆单元的长度, l 和平面情形类似, q e 与 q e 之间存在以下转换关系

? u ? ? cos( x, x) cos( x, y ) cos( x, z ) qe = ? 1 ? = ? 0 0 0 ? u2 ? ? ?

? u1 ? ? ? ? v1 ? ? ? ? ? w1 ? 0 0 0 ? cos( x, x) cos( x, y ) cos( x, z ) ? ? u2 ? ? ? ? ? v2 ? ? ? ? w2 ? ? ?
(3-61)

= Te ? q
e

e

其中 T 为坐标变换矩阵,即

? cos( x, x) cos( x, y ) cos( x, z ) Te = ? 0 0 0 ? ?
刚度矩阵和节点力的变换与平面情形相同,即为
(6×6)

? ? cos( x, x) cos( x, y ) cos( x, z ) ? ? 0 0 0
e

(3-62)

K = T
e e
(6×1)

eT

(6× 2)

(2× 2) (2×6)

K

T

e

(3-63) (3-64)

P = T

eT

(6×2)

(2×1)

P

e

3.2.4 杆单元分析的 MATLAB 程序
学习有限元方法的一个最佳途径,就是在充分掌握基本概念的基础上亲自编写有限元分析 程序,这就需要一个良好的编程环境或平台,MATLAB 软件就是这样一个平台[16],它以功能强 大、 编程逻辑直观、 使用方便见长, 从本节开始, 将提供有限元分析中主要单元完整的 MATLAB 程序,并给出详细的说明,有关 MATLAB 程序的基本操作可参见附录 A。

35

有限元分析基础教程

曾攀

【MATLAB 程序】3.2.4(1) 1D 杆单元的有限元分析程序(Bar1D2Node) 最简单的线性杆单元的程序应该包括单元刚度矩阵、单元组装、单元应力等几个基本计算 程序。下面给出编写的线性杆单元的四个 MATLAB 函数。
Bar1D2Node _Stiffness(E,A,L) 该函数计算单元的刚度矩阵,输入弹性模量 E,横截面积 A 和长度 L,输出单元刚度矩阵 k(2X2)。 Bar1D2Node _Assembly(KK,k,i,j) 该函数进行单元刚度矩阵的组装,输入单元刚度矩阵 k,单元的节点编号 i、j,输出整体刚度矩阵 KK。 Bar1D2Node _Stress(k,u,A) 该函数计算单元的应力, 输入单元刚度矩阵 k、 单元的位移列阵 u(2×1)以及横截面积 A 计算单元应力矢量, 输出单元应力 stress。 Bar1D2Node_Force(k,u) 该函数计算单元节点力矢量,输入单元刚度矩阵 k 和单元的位移列阵 u(2×1),输出 2×1 的单元节点力矢量 forces。

基于第 3.2.2 节的基本公式,写出具体实现以上每个函数的 MATLAB 程序如下。
%%%%%%%%%%% Bar1D2Node %% begin %%%%%%%%% function k=Bar1D2Node_Stiffness(E,A,L) %该函数计算单元的刚度矩阵 %输入弹性模量 E,横截面积 A 和长度 L %输出单元刚度矩阵 k(2X2) %--------------------------------------k=[E*A/L -E*A/L; -E*A/L E*A/L]; %%%%%%%%%%%%%%%%%%%%%%%%%% function z=Bar1D2Node_Assembly(KK,k,i,j) %该函数进行单元刚度矩阵的组装 %输入单元刚度矩阵 k,单元的节点编号 i、j %输出整体刚度矩阵 KK %----------------------------------DOF(1)=i; DOF(2)=j; for n1=1:2 for n2=1:2 KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2); end end z=KK; %-----------------------------------------------------------function stress=Bar1D2Node_Stress(k,u,A) %该函数计算单元的应力 %输入单元刚度矩阵 k, 单元的位移列阵 u(2×1) %输入横截面积 A 计算单元应力矢量 %输出单元应力 stress %----------------------------------stress=k*u/A;

36

有限元分析基础教程 %----------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%% function forces=Bar1D2Node_Force(k,u) %该函数计算单元节点力矢量 %输入单元刚度矩阵 k 和单元的位移列阵 u(2×1) %输出 2×1 的单元节点力分量 forces %----------------------------------------forces=k*u; %%%%%%%%%%% Bar1D2Node %% end %%%%%%%%%

曾攀

【MATLAB 程序】3.2.4(2) 2D 杆单元的有限元分析程序(Bar2D2Node) 编写如图 3-5 所示平面桁架单元的单元刚度矩阵、单元组装、单元应力的计算程序。编写 的平面桁架单元的四个 MATLAB 函数如下。
Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha) 该函数计算单元的刚度矩阵,输入弹性模量 E,横截面积 A,第一个节点坐标(x1,y1) ,第二个节点坐标 (x2,y2)和角度 alpha(单位是度) ,输出单元刚度矩阵 k(4X4)。 Bar2D2Node_Assembly(KK,k,i,j) 该函数进行单元刚度矩阵的组装,输入单元刚度矩阵 k,单元的节点编号 i、j,输出整体刚度矩阵 KK。 Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u) 该函数计算单元的应力,输入弹性模量 E,第一个节点坐标(x1,y1) ,第二个节点坐标(x2,y2) ,角度 alpha (单位是度)和单位节点位移矢量 u,返回单元应力标量。 Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u) 该函数计算单元的应力, 输入弹性模量 E, 横截面积 A, 第一个节点坐标 (x1,y1) 第二个节点坐标 , (x2,y2) , 角度 alpha(单位是度)和单元节点位移矢量 u,返回单元节点力。

基于第 3.2.3 节中的基本公式,可以编写出具体实现以上每个函数的 MATLAB 程序如下。
%%%%%%%%%%% Bar2D2Node %% begin %%%%%%%%%%%%%% function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha) %该函数计算单元的刚度矩阵 %输入弹性模量 E,横截面积 A %输入第一个节点坐标(x1,y1) ,第二个节点坐标(x2,y2) ,角度 alpha(单位是度) %输出单元刚度矩阵 k(4X4) %------------------------------------------------L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); x=alpha*pi/180; C=cos(x); S=sin(x); k=E*A/L*[C*C C*S -C*C -C*S; C*S S*S -C*S -S*S; -C*C -C*S C*C C*S; -C*S -S*S C*S S*S]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function z = Bar2D2Node_Assembly(KK,k,i,j) %该函数进行单元刚度矩阵的组装 %输入单元刚度矩阵 k,单元的节点编号 i、j %输出整体刚度矩阵 KK %-------------------------------------------------------DOF(1)=2*i-1;

37

有限元分析基础教程 DOF(2)=2*i; DOF(3)=2*j-1; DOF(4)=2*j; for n1=1:4 for n2=1:4 KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2); end end z=KK; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-function stress= Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u) %该函数计算单元的应力 %输入弹性模量 E,第一个节点坐标(x1,y1) ,第二个节点坐标(x2,y2) %输入角度 alpha(单位是度)和单位节点位移矢量 u %返回单元应力标量 stress %-----------------------------------------------L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); x=alpha*pi/180; C=cos(x); S=sin(x); stress=E/L*[-C -S C S]*u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function forces= Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u) %该函数计算单元的应力 %输入弹性模量 E,横截面积 A %输入第一个节点坐标(x1,y1) ,第二个节点坐标(x2,y2) ,角度 alpha(单位是度) %输入单元节点位移矢量 u %返回单元节点力 forces %------------------------------------------------------------L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); x=alpha*pi/180; C=cos(x); S=sin(x); forces= E*A/L*[-C -S C S]*u; %%%%%%%%%%% Bar2D2Node %% end %%%%%%%%%%%%%%

曾攀

3.2.5 杆结构分析的算例
【典型例题】3.2.5(1) 四杆桁架结构的有限元分析 如 图 3-8 所 示 的 结 构 , 各 杆 的 弹 性 模 量 和 横 截 面 积 都 为 E = 29.5 × 10 4 N/mm 2 ,

A = 100mm 2 ,试求解该结构的节点位移、单元应力以及支反力。

38

有限元分析基础教程

曾攀

图 3-8 四杆桁架结构

解答:对该问题进行有限元分析的过程如下。 (1) 结构的离散化与编号 对该结构进行自然离散,节点编号和单元编号如图 3-8 所示,有关节点和单元的信息见表 3-1 至表 3-3。
表 3-1 四杆桁架结构节点及坐标 节点 1 2 3 4 x 0 400 400 0 y 0 0 300 300

表 3-2 四杆桁架结构的单元编号及对应节点 单元 ① ② ③ ④ 节点 1 3 1 4 2 2 2 3 3

表 3-3 各单元的长度及轴线方向余弦 单元 ① ② ③ ④ l 400 300 500 400

nx
1 0 0.8 1

ny
0 -1 0.6 0

39

有限元分析基础教程

曾攀

(2)各个单元的矩阵描述 由于所分析的结构包括有斜杆,所以必须在总体坐标下对节点位移进行表达,所推导的单 元刚度矩阵也要进行变换,各单元经坐标变换后的刚度矩阵如下。

K (1)

u1 ↓ ?1 ? 29.5 ×104 ×100 ? 0 = ? ?1 400 ? ?0

v1 u2 ↓ ↓ 0 ?1 0 0 0 0 1 0

v2 ↓ 0? 0? ? 0? ? 0?

← u1 ← v1 ← u2 ← v2

(3-65)

K (2)

u3 v3 ↓ ↓ ?0 0 ? 29.5 ×104 ×100 ?0 1 = ?0 0 300 ? ?0 ?1

u2 v2 ↓ ↓ 0 0? 0 ?1? ? 0 0? ? 0 1?

← u3 ← v3 ← u2 ← v2

(3-66)

K (3)

u1 v1 u3 v3 ↓ ↓ ↓ ↓ ? 0.64 0.48 ?0.64 ?0.48? ? ? 29.5 ×104 × 100 ? 0.48 0.36 ?0.48 ?0.36 ? = ? ?0.64 ?0.48 0.64 0.48 ? 500 ? ? ? ?0.48 ?0.36 0.48 0.36 ? u4 ↓ ?1 ? 29.5 ×104 × 100 ? 0 = ? ?1 400 ? ?0 v4 u3 v3 ↓ ↓ ↓ 0 ?1 0 ? ? 0 0 0? 0 1 0? ? 0 0 0?

← u1 ← v1 ← u3 ← v3

(3-67)

K (4)

← u4 ← v4 ← u3 ← v3

(3-68)

(3) 建立整体刚度方程 将所得到的各个单元刚度矩阵按节点编号进行组装,可以形成整体刚度矩阵,同时将所有 节点载荷也进行组装。 刚度矩阵: K = K 节点位移: q = [u1
(1)

+ K (2) + K (3) + K ( 4)
v1 u 2 v2 u3 v3 u4 v4 ]
T

节点力: P = R + F = R x1

[

R y1

2 × 10 4

Ry2

0 ? 2.5 × 10 4

Rx 4

Ry4

]

T

40

有限元分析基础教程

曾攀

其中 ( Rx1 , R y1 ) 为节点 1 处沿 x 和 y 方向的支反力, y 2 为节点 2 处 y 方向的支反力, Rx 4 , R y 4 ) R ( 为节点 4 处沿 x 和 y 方向的支反力。 整体刚度方程为

u1

v1

u2

v2

u3

v3

u4 ↓ 0 0 0

v4

↓ ↓ ↓ ↓ 22.68 5.76 ?15.0 0 ? ? 5.76 4.32 0 0 ? ? ?15.0 0 15.0 0 ? 4 0 0 20.0 29.5 ×10 ×100 ? 0 ? ?7.68 ?5.76 0 0 6000 ? ?20.0 0 ??5.76 ?4.32 ? 0 0 0 0 ? 0 0 0 ? 0 ?

↓ ? 0? ? u1 ? ? Rx1 ? ? Rx1 ? ? ? v ? ? R ? ?R 0? ? 1 ? ? y1 ? ? y1 ? 4 ? 0? ?u2 ? ? Fx 2 ? ? 2 ×10 ? ?? ? ? ? ? ?20.0 0 0 0? ? v2 ? ? Ry 2 ? ? Ryz ? = = ? ? ?u3 ? ? Fx3 ? ?0 22.68 5.76 ?15.0 0 ?? ? ? ? ? 4? 5.76 24.32 0 0? ? v3 ? ? Fy 3 ? ? ?2.5 ×10 ? ? ?15.0 0 15.0 0? ?u4 ? ? Rx 4 ? ? Rx 4 ? ?? ? ? ? ? 0 0 0 0? ? v4 ? ? Ry 4 ? ? Ry 4 ? ?? ? ? ? ? ?
(3-69)

↓ ↓ ?7.68 ?5.76 ?5.76 ?4.32 0 0

(4) 边界条件的处理及刚度方程求解 边界条件 BC(u)为: u1 = v1 = v 2 = u 4 = v 4 = 0 ,代入方程(3-69)中,经化简后有

0 0 ? ?u2 ? ? 2 × 104 ? ?15 29.5 × 10 × 100 ? ? ? 0 22.68 5.76 ? ? u3 ? = ? 0 ? ? ?? ? 6000 ? 0 5.76 24.32 ? ? v3 ? ? ?2.5 × 104 ? ?? ? ? ? ?
4

(3-70)

对该方程进行求解,有

?u2 ? ? 0.2712 ? ?u ? = ? 0.0565 ? mm ? 3? ? ? ? v3 ? ? ?0.2225? ? ? ? ?
则所有的节点位移为

(3-71)

q = [ 0 0 0.2712 0 0.0565 ?0.2225 0 0] mm
T

(3-72)

(5) 各单元应力的计算

σ (1) = E ? B ? T ? q =

E [ ? 1 1] T ? q l ? 0 ? ? 0 ? 4 29.5 × 10 ? = 200 N/mm 2 = ?1 0 1 0 ] ? [ ? 0.2712 ? 400 ? ? ? 0 ?

41

有限元分析基础教程

曾攀

其中 T 为坐标转换矩阵;同理,可求出其它单元的应力。

σ (2) = ?218.8N/mm 2
(6) 支反力的计算

σ (3) = 52.08N/mm 2

σ (4) = 41.67N/mm 2

(3-73)

将节点位移的结果式(3-72)代入整体刚度方程(3-69)中,可求出
0 ? ? ? ? 0 ? 0? ? ? ? 0.2712 ? 0 ? ? ? 0 ? 0? ? ? 0.0565 ? ? 0? ? ? ? 0.2225 ? ?? 0? ? ? 0 ? ? 0 ? ? ? ?

? R x1 ? ? 22.68 ?R ? ? ? y1 ? 29.5 × 10 4 × 100 ? 5.76 ? Ry2 ? = ? 0 6000 ? ? ? ? Rx 4 ? ? 0 ? Ry4 ? ? 0 ? ? ? ? ? 15833.0 ? ? 3126.0 ? ? ? = ? 21879.0 ? N ? ? ? ? 4167.0 ? ? ? 0 ? ?

5.76 4.32 0 0 0

? 15.0 0 0 0 0

0 0 20.0 0 0

? 7.68 ? 5.76 0 ? 15.0 0

? 5.76 ? 4.32

0 0

0 ? 20.0 0 15.0 0 0

(3-74)
【MATLAB 算例】3.2.5(2) 四杆桁架结构的有限元分析(Bar2D2Node)

参见【典型例题】3.2.5(1)所给算例,如图 3-8 所示的结构,基于 MATLAB 平台求解该结 构的节点位移、单元应力以及支反力。 解答:对该问题进行有限元分析的过程如下。 (1) 结构的离散化与编号 对该结构进行自然离散,节点编号和单元编号如图 3-8 所示,有关节点和单元的信息见表 3-1 至表 3-3。 (2)计算各单元的刚度矩阵(基于国际标准单位) 建立一个工作目录, 将所编制的用于平面桁架单元分析的四个 MATLAB 函数放置于该工作 Bar2D2Node_Assembly, 目录中, 分别以各自函数的名称给出文件名, Bar2D2Node_Stiffness, 即: Bar2D2Node_Stress,Bar2D2Node_Forces。然后启动 MATLAB,将工作目录设置到已建立的目 录中,在 MATLAB 环境中,输入弹性模量 E、横截面积 A,各点坐标 x1,y1,x2,y2,x3,y3,x4,y4, 角度 alpha 1, alpha 2 和 alpha 3,然后分别针对单元 1, , 和 4, 2 3 调用四次 Bar2D2Node_Stiffness, 就可以得到单元的刚度矩阵。相关的计算流程如下。
>>E=2.95e11; >>A=0.0001; >>x1=0; >>y1=0; >>x2=0.4; >>y2=0; >>x3=0.4; >>y3=0.3;

42

有限元分析基础教程 >>x4=0; >>y4=0.3; >> alpha1=0; >> alpha2=90; >> alpha3=atan(0.75)*180/pi; >> k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2,alpha1) k1 = 73750000 0 -73750000 0 0 0 0 0 -73750000 0 73750000 0 0 0 0 0 >> k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3,alpha2) k2 = 1.0e+007 * 0.0000 0.0000 -0.0000 -0.0000 0.0000 9.8333 -0.0000 -9.8333 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -9.8333 0.0000 9.8333 >> k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3,alpha3) k3 = 1.0e+007 * 3.7760 2.8320 -3.7760 -2.8320 2.8320 2.1240 -2.8320 -2.1240 -3.7760 -2.8320 3.7760 2.8320 -2.8320 -2.1240 2.8320 2.1240 >> k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3,alpha1) k4 = 73750000 0 -73750000 0 0 0 0 0 -73750000 0 73750000 0 0 0 0 0

曾攀

(3) 建立整体刚度方程 由于该结构共有 4 个节点,因此,设置结构总的刚度矩阵为 KK(8×8),先对 KK 清零,然 后四次调用函数 Bar2D2Node _Assembly 进行刚度矩阵的组装。相关的计算流程如下。
>>KK=zeros(8,8); >>KK=Bar2D2Node_Assembly (KK,k1,1,2); >>KK=Bar2D2Node_Assembly (KK,k2,2,3); >>KK=Bar2D2Node_Assembly (KK,k3,1,3); >>KK=Bar2D2Node_Assembly (KK,k4,4,3) KK= 1.0e+008 * 1.1151 0.2832 -0.7375 0 0.2832 0.2124 0 0 -0.7375 0 0.7375 0.0000 0 0 0.0000 0.9833 -0.3776 -0.2832 -0.0000 -0.0000 -0.2832 -0.2124 -0.0000 -0.9833 0 0 0 0 0 0 0 0

-0.3776 -0.2832 -0.0000 -0.0000 1.1151 0.2832 -0.7375 0

-0.2832 -0.2124 -0.0000 -0.9833 0.2832 1.1957 0 0

0 0 0 0 -0.7375 0 0.7375 0

0 0 0 0 0 0 0 0

(4) 边界条件的处理及刚度方程求解 由图 3-8 可以看出,节点 1 的位移将为零,即 u1 = 0 , v1 = 0 ,节点 2 的位移 v2 = 0 ,节 点 4 的 u4 = 0 , v4 = 0 。节点载荷 F3 =10N。采用高斯消去法进行求解,注意:MATLAB 中的反

43

有限元分析基础教程

曾攀

斜线符号“\”就是采用高斯消去法。 该结构的节点位移为: q = [u1 而节点力为: P = R + F = R x1

v1 u 2

v2

u3

v3

u4

v4 ]

T

[

R y1

2 × 10 4

Ry2

0 ? 2.5 × 10 4

Rx 4

Ry4

]

T

其中 ( Rx1 , R y1 ) 为节点 1 处沿 x 和 y 方向的支反力, y 2 为节点 2 处 y 方向的支反力, Rx 4 , R y 4 ) R ( 为节点 4 处沿 x 和 y 方向的支反力。相关的计算流程如下。
>>k=KK([3,5,6],[3,5,6]) k =1.0e+008 * 0.7375 -0.0000 -0.0000 -0.0000 1.1151 0.2832 -0.0000 0.2832 1.1957 >> p=[20000;0;-25000]; >>u=k\p u = 1.0e-003 * 0.2712 0.0565

-0.2225

[这里将列排成了一行,以节省篇幅]

由此可以看出,所求得的结果 u2 = 0.2712mm, u3 = 0.0565mm, v3 = ?0.2225mm ,则所有节 点位移为

q = [ 0 0 0.2712 0 0.0565 ?0.2225 0 0] mm
T

(3-75)

与前面通过数学推导得到的结果相同,见式(3-72)。 (5)支反力的计算 在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力。将整体的 位移列阵 q(采用国际标准单位)代回原整体刚度方程,计算出所有的节点力 P,按上面的对应 关系就可以找到对应的支反力。相关的计算流程如下。
>> q=[0 0 0.0002712 0 0.0000565 -0.0002225 0 0]' q = 1.0e-003 * 0 0 0.2712 0 0.0565 -0.2225 0 0 [这里将列排成了一行,以节省篇幅] >>P=KK*q P = 1.0e+004 * -1.5833 0.3126 2.0001 2.1879 -0.0001 -2.5005 -0.4167 0 [这里将列排成了行]

按对应关系,可以得到对应的支反力为

(3-76)

44

有限元分析基础教程

曾攀

(6)各单元的应力计算 先从整体位移列阵 q 中提取出单元的位移列阵,然后,调用计算单元应力的函数 Bar2D2Node_Stress , 就 可 以 得 到 各 个 单 元 的 应 力 分 量 。 当 然 也 可 以 调 用 上 面 的 Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)函数来计算单元的集中力,然后除以面积求得单元 应力。相关的计算流程如下。
>>u1=[q(1);q(2);q(3);q(4)] u1 = 1.0e-003 * 0 0 0.2712 0 >> stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha1,u1) stress1 =2.0001e+008 >>u2=[q(3);q(4);q(5);q(6)] u2 = 1.0e-003 * 0.2712 0 0.0565

[这里将列排成了一行,以节省篇幅]

-0.2225

[这里将列排成了一行,以节省篇幅]

>> stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,alpha2,u2) stress2 = -2.1879e+008 >>u3=[q(1);q(2);q(5);q(6)] u3 = 1.0e-003 * 0 0 0.0565 -0.2225 >> stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,alpha3,u3) stress3 = -52097000 >>u4=[q(7);q(8);q(5);q(6)] u4 = 1.0e-003 * 0 0 0.0565 -0.2225 >> stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,alpha1,u4) stress4 = 41668750

[这里将列排成了一行,以节省篇幅]

[这里将列排成了一行,以节省篇幅]

可 以 看 出 : 计 算 得 到 的 单 元 1 的 应 力 为 σ = 2.0001× 108 Pa ; 单 元 2 的 应 力 为

σ = ?2.1879 ×108 Pa , 单 元 3 的 应 力 为 σ = ?5.2097 ×107 Pa , 单 元 4 的 应 力 为 σ = 4.167 ×107 Pa 。与前面通过数学推导得到的结果相同,见式(3-73)。
ANSYS 是大型的通用有限元分析系统,从本节开始也针对本书的典型例题提供完整的 ANSYS 操作流程,包括基于图形界面的操作以及基于命令流的操作。这样将使得以基于详细推 导的典型例题与基于 MATLAB 的编程实现、以及与基于 ANSYS 的分析都完整地结合起来,可 以更好的理解和使用有限元方法这一工具[20]。有关 ANSYS 软件的基本操作可参见附录 B。
【ANSYS 算例】3.2.5(3) 四杆桁架结构的有限元分析 下面针对【典型例题】3.2.5(1)的问题,在 ANSYS 平台上,完成相应的力学分析。 解答:对该问题进行有限元分析的过程如下。 以下为基于 ANSYS 图形界面(GUI, graphic user interface)的菜单操作流程; 注意: 符号 →” “

45

有限元分析基础教程

曾攀

表示针对菜单中选项的鼠标点击操作。 1 基于图形界面(GUI)的交互式操作(step by step)
(1) 进入 ANSYS(设定工作目录和工作文件) 程序 →ANSYS → ANSYS Interactive →Working directory (设置工作目录) →Initial jobname(设置工作 文件名): planetruss→Run → OK (2) 设置计算类型 ANSYS Main Menu: Preferences… → Structural → OK (3) 选择单元类型 ANSYS Main Menu: Preprocessor →Element Type→Add/Edit/Delete… →Add… →Link:2D spar 1 OK (返回到 Element Types 窗口) →Close (4) 定义材料参数 ANSYS Main Menu: Preprocessor →Material Props →Material Models→Structural →Linear →Elastic → Isotropic: EX:2.95e11 (弹性模量),PRXY: 0 (泊松比) → OK → 鼠标点击该窗口右上角的“ 窗口 (5) 定义实常数以确定单元的截面积 ANSYS Main Menu: Preprocessor →Real Constants… →Add/Edit/Delete →Add →Type 1→ OK→Real Constant Set No: 1 (第 1 号实常数), AREA: 1e-4 (单元的截面积) →OK→Close (6) 生成单元 ANSYS Main Menu: Preprocessor →Modeling →Creat→Nodes→ In Active CS→Node number 1 → X:0,Y:0,Z:0 →Apply→Node number 2 → X:0.4,Y:0,Z:0 →Apply→Node number 3 → X:0.4,Y:0.3,Z:0→ Apply→Node number 4 → X:0,Y:0.3,Z:0→OK ANSYS Main Menu: Preprocessor → Modeling → Create → Elements→Elem Attributes (接受默认值) →User numbered→Thru nodes→ OK→选择 node 1 和 node2→ Apply→选择 node 2 和 node3→ Apply→选择 node 1 和 node3→ Apply→选择 node 4 和 node3→ Apply→OK (7) 模型施加约束和外载 添加位移的约束,分别将 1 节点 X 和 Y 方向、2 节点 Y 方向、4 节点的 X 和 Y 方向位移约束。 ANSYS Main Menu: Solution → Define Loads → Apply →Structural → Displacement → On Nodes → 用鼠标选择节点 1→ Apply → Lab2 DOFs: UX,UY ,VALUE:0 →Apply→用鼠标选择节点 2→ Apply → Lab2 DOFs: UY , VALUE: →Apply→用鼠标选择节点 4→ Apply → Lab2 DOFs: UX,UY , 0 VALUE: → 0 OK 加载集中力 ANSYS Main Menu: Solution → Define Loads → Apply →Structural → Force/moment→ On Nodes →用鼠标选择结构节点 2→ Apply →FX, VALUE: 20000 FY,VALUE: -25000 →OK (9) 分析计算 ANSYS Main Menu: Solution → Solve → Current LS →OK → Should The Solve Command be Executed? Y→ Close (Solution is done! ) → 关闭文字窗口 (10) 结果显示 ANSYS Main Menu: General Postproc → Plot Results →Deformed Shape … → Def + Undeformed → OK (返回到 Plot Results) → Contour Plot → Nodal Solu … → DOF solution→Displacement vector sum (可 以看到位移云图) → Apply→用鼠标选择结构节点 3→ Apply → ”来关闭该 →

46

有限元分析基础教程

曾攀

ANSYS Main Menu: General Postproc →List Results →Nodal solution →DOF solution→Displacement vector sum (弹出的文本文件显示各个节点的位移) ANSYS Main Menu: General Postproc →List Results →Reaction Solu→ALL items→ OK (弹出的文本文 件显示各个节点反力) (11) 退出系统 ANSYS Utility Menu: File→ Exit …→ Save Everything→OK

2 完整的命令流 以下为命令流语句;注意:以“!”打头的文字为注释内容,其后的文字和符号不起运行作 用。关于命令流的调用方式见附录 B。 !%%%%%%%% [典型例题]3.2.5(1) %%%% begin %%%%%%
/ PREP7 /PLOPTS,DATE,0 ET,1,LINK1 UIMP,1,EX, , ,2.95e11, R,1,1e-4, N,1,0,0,0, N,2,0.4,0,0, N,3,0.4,0.3,0, N,4,0,0.3,0, E,1,2 E,2,3 E,1,3 E,4,3 FINISH /SOLU D,1,ALL D,2,UY, D,4,ALL F,2,FX,20000, F,3,FY,-25000, SOLVE FINISH /POST1 PLDISP,1 FINISH !进入前处理 !设置不显示日期和时间 !选择单元类型 !给出材料的弹性模量 !给出实常数(横截面积) !生成 1 号节点,坐标(0,0,0) !生成 2 号节点,坐标(0.4,0,0) !生成 3 号节点,坐标(0.4,0.3,0) !生成 4 号节点,坐标(0,0.3,0) !生成 1 号单元(连接 1 号节点和 2 号节点) !生成 2 号单元(连接 2 号节点和 3 号节点) !生成 3 号单元(连接 1 号节点和 3 号节点) !生成 4 号单元(连接 4 号节点和 3 号节点) !前处理结束 !进入求解状态(在该状态可以施加约束及外力) !将 1 号节点的位移全部固定 !将 2 号节点的 y 方向位移固定 !将 4 号节点的位移全部固定 !在 2 号节点处施加 x 方向的力(20000) !在 3 号节点处施加 y 方向的力(-25000) !进行求解 !结束求解状态 !进入后处理 !显示变形状况 !结束后处理

!=====设置单元、材料,生成节点及单元

!=====在求解模块中,施加位移约束、外力,进行求解

!=====进入一般的后处理模块

!%%%%%%%% [典型例题]3.2.5(1) %%%% end %%%%%%

3.3 梁件有限元分析的标准化表征与算例
47

有限元分析基础教程

曾攀

3.3.1 梁件分析的基本力学原理
设有一个受分布载荷作用的简支梁如图 3-9 所示,由于简支梁的宽度较小,外载沿宽度方 向无变化,该问题可以认为是一个 xoy 平面内的问题,可以有以下两种方法来建立基本方程。 方法之一是采用一般的建模及分析方法,即从对象取出 dxdy 微元体进行分析,建立最一般 的方程,见第 4.2 节中关于 2D 问题的基本变量及方程,这样,所用的变量较多,方程复杂,未 考虑到这一具体问题的特征。 方法之二是针对细长梁用 “特征建模” characterized modeling) 的简化方法来推导三大方程, ( 其基本思想是采用工程宏观特征量来进行问题的描述;图 3-9 所示问题的特征为:1梁为细长 , 2主要变形为垂直于 x 的挠度, 可只用挠度(deflection) 梁 long beam)因此可只用 x 坐标来刻画, ( 来描述位移场;针对这两个特征,可以对梁沿高度方向的变形做出以下设定:(1)变形后的直线 假定;(2)小变形假定。这两个假定对于细长梁的实际情况也是符合的。 【基本变量】3.3.1(1) 平面梁的基本变量 位移: 应力:

? v ( x , y = 0) (中性层的挠度)

σ (采用 σ x ,其它应力分量很小,不考虑),该变量对应于梁截面上的弯矩 M 应变: ε (采用 ε x ,沿高度方向满足直线假定)
下面取具有全高度梁的 dx “微段” 来推导三大方程(见图 3-10)。

图 3-9 受分布载荷作用的简支梁

图 3-10 梁问题的 dx“微段”及受力平衡

【基本方程】3.3.1(2) 平面梁的基本方程 该问题的三大类基本方程和边界条件如下。 (1) 平衡方程 针对图 3-10 中的“微段”,应该有三个平衡方程,首先由 x 方向的合力等效 ∑ X = 0 ,有

? M = ∫ σ x ? y ? dA
A

(3-77)

? 其中 y 是以梁的中性层(neutral layer)为起点的 y 坐标,M 为截面上的弯矩。
然后由 y 方向的合力平衡 ∑ Y = 0 ,有 dQ + p ( x) ? dx = 0 ,即

dQ +p = 0 dx

(3-78)

48

有限元分析基础教程

曾攀

其中 Q 为截面上的剪力,再由弯矩平衡 ∑ M = 0 ,有 dM ? Qdx = 0 ,即 0

Q=

dM dx

(3-79)

(2) 几何方程 梁问题 dx“微段”的弯曲变形分析见图 3-11。

图 3-11 梁的纯弯变形图

考虑梁的纯弯变形(pure bending deformation),如图 3-11 所示。由变形后的几何关系,可得

? 到位于 y 处纤维层的应变(即相对伸长量)为
? ε x ( y) = ? ( R ? y ) ? dθ ? R ? dθ R ? dθ ? y =? R
dθ dθ 1 = = ds R ? dθ R

(3-80)

其中 R 为曲率半径,而曲率 κ (curvature)与曲率半径 R 的关系为

κ=

(3-81)

? 对于梁的挠度函数 v ( x , y = 0) ,它的曲率 k 的计算公式为

κ =±
这里就图 3-11 所示的情形,应取为

v '' ( x) ≈ ± v ''( x ) (1 + v ' ( x) 2 )3/ 2

(3-82)

κ=
将式(3-81)和式(3-83)代入式(3-80)中,有

d 2v dx 2

(3-83)

? ε x ( x, y ) ? =?y
(3) 物理方程 利用虎克定律

d 2v dx 2

(3-84)

49

有限元分析基础教程

曾攀

σ x = Eε x
并对以上方程进行整理, 有描述平面梁弯曲问题的基本方程

(3-85)

? EI

d 4v + p( x) = 0 dx 4
A

(y 方向的平衡)

(3-86)

? M ( x ) = ∫ σ x yd A ? = ∫ ? y 2 Ev ''dA
A

(x 方向的平衡)

(3-87)

= ? EI

d 2v dx 2

? σ x ( x) = ? E y ? ε x ( x) = ? y
方程(3-87)中 I =

d 2v dx 2

(物理方程)

(3-88)

d 2v dx 2

(几何方程)

(3-89)



A

? y 2 d A 为梁截面的惯性矩(moment of inertia),可以看出:将原始基本变量定为

? 中性层的挠度 v ( x , y = 0) ,而其它力学参量都可以基于它来表达。
(4) 边界条件 图 3-11 所示简支梁的边界为梁的两端,由于在建立平衡方程时已考虑了分布外载 p ( x ) (见 方程(3-86)),因此不能再作为力的边界条件。 两端的位移边界:

BC(u):
两端的力(弯矩)边界:

v x = 0 = 0 , v x =l = 0

(3-90)

BC(p):

M

x =0

= 0,M

x =l

=0

(3-91)

由式(3-87),可将弯矩以挠度的二阶导数来表示,即

BC(p):

v ′′ |x =0 = 0

v ′′ |x =l = 0

(3-92)

【求解原理】3.3.1(3) 简支梁的微分方程解 若用基于 dxdy 微体所建立的原始方程 (即原平面应力问题中的三大类方程) 进行直接求解, 不仅过于繁琐,而且不易求解,若用基于以上“特征建模” 简化方法所得到的基本方程进行直 接求解则比较简单,对如图 3-9 所示的均匀分布外载的情况,其方程为:

50

有限元分析基础教程

曾攀

? EI
BC(u): BC(p):

d 4v + p0 = 0 dx 4

(3-93) (3-94) (3-95)

v x =0 = v x =l = 0 v ′′ x =0 = v ′′ x =l = 0

这是一个常微分方程,其解的形式为

v( x) =

1 p0 4 ( x + c3 x3 + c2 x 2 + c1 x + c0 ) EI 24

(3-96)

其中 c 0 … c3 为待定系数,可由四个边界条件求出,最后有结果

v( x) =
则位于中点处的挠度为

p0 ( x 4 ? 2lx3 + l 3 x) 24 EI

(3-97)

v( x = 1 l ) = 0.013020833 2

p0l 4 EI

(3-98)

计算平面梁弯曲问题有关能量方面的物理量如下。 应变能

U=

1 σ x ε x dΩ 2 ∫Ω ∧ d 2v 1 d 2v ? = ∫ (? E y 2 )(? y 2 )dAdx 2 Ω dx dx ? d 2v ? 1 = ∫ EI z ? 2 ? dx 2 l ? dx ?
2

(3-99)

外力功

W = ∫ p ( x) ? v( x)dx
l

(3-100)

势能

Π = U ?W =

1 d 2v EI z ( 2 )2 dx ? ∫ p ( x) ? v( x)dx l 2 ∫l dx

(3-101)

【求解原理】3.3.1(4) 简支梁的虚功原理求解

? 同样以如图 3-9 所示的简支梁为例, 假设有一个只满足位移边界条件 BC(u)的位移场 v ( x )


? v( x) = c1 ? sin

πx
l

(3-102)

51

有限元分析基础教程

曾攀

其中 c1 为待定系数。则它的微小变化,即虚位移场为

? δ v( x) = δ c1 ? sin

πx
l

(3-103)

δc1 为微小变化量,可以验证,式(3-102)满足位移边界条件 BC(u),将满足位移边界条件 BC(u)
的试函数叫做许可位移(admissible displacement)。 该简支梁的虚应变能为

δ U = ∫ σ xδε x dΩ
Ω l

=∫

0



A

E ? ε xδε x dA ? dx

(3-104)

其中 A 为梁的横截面,对于梁的弯曲问题,由式(3-84),有几何方程

? εx = ?y ?
将其代入式(3-104)中,则有

? d 2v 2 dx

(3-105)

δU = ∫ E
0

l

(
2

? d 2 v ? ? d 2δ v ? ? y 2dA ? ? 2 ? ? ? 2 ? ? dx ∫A ? dx ? ? dx ? ? ? ? ?

)

(3-106)

将式(3-102)和式(3-103)代入上式,有

πx ?π ? πx ?π ? δ U = ∫ EI ? ? ? c1 sin ? ? ? ? sin ? δ c1 ? dx 0 l ?l ? l ?l ?
l

2

EIl ? π ? = ? ? ? c1 ? δ c1 2 ?l ?
其中 I =

4

(3-107)



A

? y 2 dA 为截面惯性矩。

该简支梁的外力虚功为

? δ W = ∫ p0δ vdx
0

l

= p0 ? δ c1 ? ∫ sin
0

l

πx
l

dx
(3-108)

=
由虚功原理(3-12),即 δU = δW ,则

2lp0

π

? δ c1

2lp0
消去 δc 1 后,有

EIl ? π ? ? δ c1 = ? ? ? c1 ? δ c1 π 2 ?l ?

4

(3-109)

52

有限元分析基础教程

曾攀

c1 =

4l 4 p0 EI π 5 4l 4 πx p0 sin 5 EI π l

(3-110)

那么,由式(3-102)所表示的位移模式中,真实的一组为满足虚功原理时的位移,即

v ( x) =

(3-111)

【求解原理】3.3.1(5) 简支梁的最小势能原理求解 仍以如图 3-9 所示的平面简支梁的弯曲问题为例,为提高计算精度,可以选取多项函数的

? 组合,这里取满足位移边界条件 BC(u)的许可位移场 v ( x ) 为

? v( x) = c1 ? sin
其中 c1 和 c2 为待定系数;计算应变能 U 为

πx
l

+ c2 ? sin

3π x l

(3-112)

U=

1 σ x ? ε x ? dΩ 2 ∫Ω
2

? ? d 2v ? 1 l = ∫ EI ? ? 2 ? dx 0 2 ? dx ?
4 4 2 2 1 l ? 2 ?π ? πx 3π x ? ? π ? ? 3π ? 2?πx? 2 ? 3π ? 2 ? 3π x ? = ∫ EI ? c1 ? ? sin ? sin ?dx ? + c2 ? ? sin ? ? + 2c1c2 ? ? ? ? sin 2 0 ? ?l ? l l ? ? l ? ? l ? ? l ? ?l ? ? l ? ? ?

EI = 2

4 ? 2 ? π ?4 l ? 2 ? 3π ? l ?c1 ? ? + c2 ? ? ? ? l ? 2? ? ?l? 2 ? ?

(3-113)

相应的外力功 W 为
l πx 3π x ? ? + c2 sin W = ∫ p0 ? ? c1 sin ?dx 0 l l ? ? 2l ? ? 2l = p0 ? c1 + c2 3π ? ? π ?

(3-114)

则总势能为 Π = U ? W ,为使 Π 取极小值,则有
4 2l ?Π EI ? ? π ? l ? 2c1 ? ? ? ? p0 = 0 = ? π ?c1 2 ? ? l ? 2? ? ? 4 2l ?Π EI ? ? 3π ? l ? =0 = 2c2 ? ? ? ? p0 ? ?c2 2 ? ? l ? 2? 3π ? ?

? ? ? ? ? ? ?

(3-115)

53

有限元分析基础教程

曾攀

? 解出 c1 和 c2 后,有 v ( x ) 的具体表达
? v( x) = 4 p0l 4 4 p0l 4 ?πx? ? 3π x ? sin ? sin ? + ? ? 5 5 π EI ? l ? 243π EI ? l ?
(3-116)

可以看出,该方法得到的第一项与前面虚功原理求解出来的结果相同,与精确解(3-97)相比,该 结果比前面由虚功原理得到的结果更为精确,这时因为选取两项函数作为试函数,这也是提高 计算精度的重要途径。以上求解过程所用的试函数(3-112)为许可基底函数的线性组合,因此, 上述求解方法也是瑞利-里兹方法。 以上的【求解原理】3.3.1(4)以及【求解原理】3.3.1(5)都是基于试函数的能量方法(也称为 泛函极值方法),基本要点是不需求解原微分方程,但需要假设一个满足位移边界条件 BC(u)的 许可位移场。因此,如何寻找或构建满足所需要求的许可位移场是一个关键,并且,还期望这 种构建许可位移场的方法还应具有标准化和规范性。下面的重点将讨论通过基于“单元”的位 移函数的构建就可以满足这些要求。

3.3.2 局部坐标系中的平面梁单元
【单元构造】3.3.2(1) 平面纯弯梁单元的描述 图 3-12 所示为一局部坐标系中的纯弯梁单元(beam element),其长度为 l,弹性模量为 E, 横截面的惯性矩为 Iz。

图 3-12 局部坐标系中的梁单元

(1) 单元的几何及节点描述 设有两个端节点(Node 1 和 Node 2) ,节点位移列阵 q 为
e

q e = [ v1 θ1

v2 θ 2 ]
e

T

(3-117)

这表明该单元的节点位移有 4 个自由度(DOF)。节点力列阵 P 为

P e = [ Pv1

M1

Pv 2

M2]

T

(3-118)

其中 v1 , θ1 , v2 , θ 2 分别为各节点的挠度(deflection)和转角(slope),若该单元承受有分布外载,可 以将其等效到节点上,即也可以表示为如式(3-118)所示的节点力;和前面推导杆单元时的情形

54

有限元分析基础教程

曾攀

类似,利用函数插值、几何方程、物理方程以及势能计算公式,可以将单元的所有力学参量用 节点位移列阵 q 及相关的插值函数来表示。
e

(2) 单元位移场的表达 由于有 4 个位移节点条件,可假设纯弯梁单元的位移场挠度为具有四个待定系数的函数模 式,即

v( x) = a0 + a1 x + a2 x 2 + a3 x3
其中 a0,a1,a2,a3 为待定系数。由该单元的节点位移条件

(3-119)

v( x = 0) = v1 , v ' ( x = 0) = θ1 v( x = l ) = v2 , v ( x = l ) = θ 2
'

? ? ? ? ?

(3-120)

可求出式(3-119)中的 4 个待定系数,即

a0 = v1 a1 = θ1 a2 = 1 (?3v1 ? 2θ1l + 3v2 ? θ 2l ) l2 a3 = 1 (2v1 + θ1l ? 2v2 + θ 2l ) l3
将式(3-121)代入式(3-119)中,重写位移函数,有

? ? ? ? ? ? ? ? ?

(3-121)

v( x) = (1 ? 3ξ 2 + 2ξ 3 )v1 + l (ξ ? 2ξ 2 + ξ 3 )θ1 + (3ξ 2 ? 2ξ 3 )v2 + l (ξ 3 ? ξ 2 )θ 2 = N(ξ ) ? q e
其中

(3-122)

ξ = , N (ξ ) 叫做单元的形状函数矩阵,即

x l

N(ξ ) = ?(1 ? 3ξ 2 + 2ξ 3 ) l (ξ ? 2ξ 2 + ξ 3 ) (3ξ 2 ? 2ξ 3 ) l (ξ 3 ? ξ 2 )? ? ?
(3) 单元应变场的表达 由纯弯梁的几何方程,有梁的应变表达式
^ ^ ε ( x, y ) = ? y

(3-123)

d 2 v( x) dx 2 1 1 1 ? ^ ?1 (6ξ ? 4) ? 2 (12ξ ? 6) (6ξ ? 2) ? ? q e = ? y ? 2 (12ξ ? 6) l l l ?l ? = B(ξ ) ? q e

(3-124)

^ 其中 y 是以中性层为起点的 y 方向的坐标, B (ξ ) 叫做单元的几何矩阵,即

55

有限元分析基础教程
^ B(ξ ) = ? y [ B1

曾攀

B2

B3

B4 ]

(3-125)

其中

B1 = 1 (12ξ ? 6) , l2
(4) 单元应力场的表达 由梁的物理方程

1 1 B2 = (6ξ ? 4) , B3 = ? 2 (12ξ ? 6) , l l

1 B4 = (6ξ ? 2) l

^ ^ ^ ^ σ ( x, y ) = E ? ε ( x, y ) = E ? B ( x, y ) ? q e = S ( x, y ) ? q e

(3-126)

其中 E 为弹性模量, S ( x ) 叫做单元的应力矩阵。

(5) 单元势能 Π 的表达
e

该单元的势能为

Πe = U e ?W e
其中应变能

(3-127)

Ue =

1 l ^ ^ y) σ ( x, ? ε ( x, y) ? dA ? dx 2 ∫0 ∫A l 1 = q e T [ ∫ ∫ BT ? E ? B ? dA ? dx]q e 0 A 2 1 = qeT ? K e ? qe 2

(3-128)

其中 K 为刚度矩阵,具体地有
e

? B1 ? ?B ? l ^ K e = ∫ ∫ ( ? y ) ? 2 ? ? E ? [ B1 0 A ? B3 ? ? ? ? B4 ? ? ?

B2

B3

^ B4 ] ( ? y ) ? dA ? dx

? B12 ? l BB ^ 2 = ∫ ( ? y ) dA ? E ? ∫ ? 1 2 A 0 ?B B 1 3 ? ? B1 B4 ?

B1 B2 2 B2 B2 B3 B2 B4

B1 B3 B2 B3 B32 B3 B4

B1 B4 ? ? B2 B4 ? ? dx B3 B4 ? ? B42 ? ?

56

有限元分析基础教程

曾攀

6l ? 12 ? 6l 4l 2 E?I = 3z? l ? ?12 ?6l ? 2 ? 6l 2l

?12 ?6l 12 ?6l

6l ? ? 2l 2 ? ?6l ? ? 4l 2 ?

(3-129)

Iz 为惯性矩,式(3-127)中的外力功为 W e = Pv1 ? v1 + M 1θ1 + Pv 2 ? v2 + M 2θ 2 = P e T ? q e
其中

(3-130)

P e = [ Pv1
(6) 单元的刚度方程

M1

Pv 2

M2]

T

(3-131)

同样,由最小势能原理,将式(3-127)中的 Π 对 q 取极小值,有单元刚度方程
e
e

(4×4) 4×1) (4×1) (

K e ? qe = Pe

(3-132)

其中刚度矩阵 K 和力矩阵 P 分别见式(3-129)和式(3-130)。
e e

【单元构造】3.3.2(2) 一般平面梁单元的描述 为推导局部坐标系中的一般平面梁单元,在图 3-13 所示的纯弯梁的基础上叠加进轴向位 移(由于为线弹性问题,满足叠加原理),这时的节点位移自由度(DOF)共有 6 个,见图 3-13。

图 3-13 平面梁单元

图 3-13 所示平面梁单元的节点位移列阵 q 和节点力列阵 P 为
e
e

(6×1)

q e = [u1 v1 θ1 u1 v2 θ 2 ]
Pv1 M1 Pu 2 Pv 2

T

(3-133)
T

(6×1)

P e = [ Pu1

M2 ]

(3-134)

相应的刚度方程为
(6×6) 6×1) (6×1) (

K e ? qe = Pe

(3-135)

对应于图 3-13 的节点位移和式(3-133)式中节点位移列阵的排列次序,将杆单元刚度矩阵与纯弯 梁单元刚度矩阵进行组合,可得到式(3-135)中的单元刚度矩阵,即

57

有限元分析基础教程

曾攀

? EA ? l ? ? 0 ? ? ? 0 ? e K = (6×6) ? EA ? ? l ? ? ? 0 ? ? ? 0 ? ?

0 12 EI l3 6 EI l2 0 ? 12 EI l3 6 EI l2 ?

0 6 EI l2 4 EI l 0 6 EI l2 2 EI l

?

EA l 0 0 ?

0 12 EI l3 6 EI ? 2 l 0 12 EI l3 6 EI ? 2 l ?

0 6 EI l2 2 EI l 0 6 EI l2 4 EI l

EA l 0 0

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

(3-136)

【典型例题】3.3.2(3) 受均布载荷平面梁单元的等效节点载荷 工程中的梁构件,存在许多象图 3-14(a)所示的受均布载荷的情况,若就受均布载荷部分 的梁构件建立单元,则需要就所建立的梁单元给出相应的节点等效载荷,就该问题进行讨论。

(a) 几种受均布载荷作用的梁构件 图 3-14

(b) 节点等效载荷

受均布载荷作用的梁构件以及节点等效载荷

图 3-15 受均布载荷作用的直接静力等效的节点载荷(每个节点分一半)

解答:对于图 3-14(a)所示的几种受均布载荷的情况,若需要采用几个或一个梁单元,均可以建 立 如 图 3-14(b) 所 示 的 单 元 , 其 节 点 位 移 列 阵 为 q e = [ v1 θ1

v2 θ 2 ] , 节 点 力 列 阵 为
T

P e = [ Pv1

M1

Pv 2

M 2 ] ,单元的挠度位移场 v ( x ) 为 v ( x ) = N ( x ) ? q e ,其中
T

58

有限元分析基础教程

曾攀

N ( x ) = ?(1 ? 3ξ 2 + 2ξ 3 ) l (ξ ? 2ξ 2 + ξ 3 ) ?
而ξ =

( 3ξ

2

? 2ξ 3 ) l (ξ 3 ? ξ 2 ) ? ?

(3-137)

x ;计算该单元的外力功为 W 为 l
W = ∫ p ( x ) ? v ( x ) dx
l 0 l = ? ∫ ( ? p0 ) ? N ( x ) ? dx ? ? q e ? 0 ? ? ? eT e = P ?q

(3-138)

将式(3-137)的形状函数矩阵代入上式,可计算出节点力 P 为
e

P eT = ∫ ( ? p0 ) ? N ( x ) ? dx
l 0

= p0 ∫ ?(1 ? 3ξ 2 + 2ξ 3 ) l (ξ ? 2ξ 2 + ξ 3 ) 0?
l

( 3ξ

2

? 2ξ 3 ) l (ξ 3 ? ξ 2 ) ? dx ?
(3-139)

= ? ? p0 l 2 ? p0 l 2 12 ? p0 l 2 ?

p0 l 2 12 ? ?

讨论 1:若凭一种直觉,直接按照静力等效的方式来进行计算,即,每个节点各分一半进行静 力等效,见图 3-15,则计算出的节点等效力为

Pv1 = ?

p l2 p l2 p l2 p0l , M 1 = ? 0 , Pv 2 = ? 0 , M 2 = ? 0 2 8 2 8

(3-140)

将该结果与式(3-139)进行比较,显然这样计算出的 M1 和 M2 都是错误的。 讨论 2:该等效节点载荷是按照外力功(3-138)式进行计算的,是通用的均布载荷的节点等效 载荷,与节点的实际约束状态没有关系。也就是说,图 3-14(a)中的几种情况的节点等效载荷 。 都用式(3-139)

【典型例题】3.3.2(4)

悬臂-简支平面连续梁的有限元分析

有 一 个 连 续 平 面 悬 臂 简 支 梁 , 截 面 面 积 A = 6650mm 2 , 高 度 h = 317mm , 惯 性 矩

I = 118.6 ×106 mm4 , 如 图 3-16 所 示 , 梁 的 弹 性 模 量 E = 200 GPa , 承 受 的 均 布 荷 载
p0 = 25 000 N m 。试确定节点 3 的竖向位移、节点 2 和节点 3 的转角。同时计算节点 1 和节
点 2 的反力。

59

有限元分析基础教程 (a) 连续平面悬臂简支梁 (b)有限元分析模型

曾攀

图 3-16 连续平面悬臂简支梁以及有限元分析模型

解答:由于该梁在其中的一个位置有一个支撑,因此采用两个梁单元。则该结构的整体节点位 移列阵

q = [ v1 θ1 v2 θ 2 v3 θ3 ]T
载荷列阵为

(3-141)

P = F + R = [ Fv1 +Rv1

Fθ 1 +Rθ 1

Fv 2 +Rv 2

Fθ 2

Fv 3

Fθ 3 ]T

(3-142)

其中 F 为施加的外载,见式(3-148), R 为支反力。 由前面给出的梁的单元刚度矩阵的计算公式(3-129),代入单元(1)的属性值,有

v1

θ1

v2

θ2
(3-143)

k( )
1

6 × (5) 6 × ( 5) ? ?12 ? 12 ? 2 2? 200 × 109 × 1.186 × 10?4 ? 6 × ( 5 ) 4 × ( 5 ) ?6 × ( 5 ) 2 × ( 5 ) ? = ×? 12 ?12 ?6 × ( 5 ) ?6 × ( 5 ) ? 53 ? ? ? 6 × ( 5 ) 2 × ( 5 ) 2 ?6 × ( 5 ) 4 × ( 5 ) 2 ? ? ?

计算单元(2)的刚度矩阵,有

v2

θ2

v3

θ3
(3-144)

k ( 2)

?12 6 × ( 2.5 ) 6 × ( 2.5 ) ? ? 12 ? 2 2? 200 × 109 ×1.186 × 10?4 ?6 × ( 2.5 ) 4 × ( 2.5 ) ?6 × ( 2.5 ) 2 × ( 2.5 ) ? = ×? 3 12 ?12 ?6 × ( 2.5 ) ?6 × ( 2.5 ) ? ( 2.5) ? ? ?6 × ( 2.5 ) 2 × ( 2.5 )2 ?6 × ( 2.5 ) 4 × ( 2.5 )2 ? ? ?
( 2)

计算总刚度矩阵 K = k

+ k ( 2) ,有
v2

v1

θ1

θ2

v3

θ3

5 692 800 0 0 ? 2 277 120 5 692 800 ?2 277 120 ? ? 5 692 800 18 976 000 ?5 692 800 ? 9 488 000 0 0 ? ? ? ?2 277 120 ?5 692 800 20 494 080 17 078 400 ?18 216 960 22 771 200 ? K=? ? 56 928 000 22 771 200 18 216 960 ? ? 5 692 800 9 488 000 17 078 400 ? 0 0 ?18 216 960 ?22 771 200 18 216 960 ?22 771 200 ? ? ? 0 0 22 771 200 18 976 000 ?22 771 200 37 952 000 ? ? ? ?
(3-145) 参考【典型例题】3.3.2(3),对单元的分布载荷进行处理,分别计算单元(1)和单元(2)的载 荷列阵,即

60

有限元分析基础教程
3 ? wL ? ? 25 ×10 × 5 ? ? ? ? 2 ? ? 2 ? ? ? ? 2 3 2 ? ? wL ? ? 25 ×10 × 5 ? ? ?62 500 ? ? ? ?52 083 ? ? 12 ? ? 12 ? =? ?=? ?=? wL ? ? 25 ×103 × 5 ? ? ?62 500 ? ?? ? ? 52 083 ? ? 2 ? ? 2 ? ? ? ? ? ? wL2 ? ? 25 ×103 × 52 ? ? ? 12 ? ? ? ? ? ? ? 12

曾攀

F (1)

(3-146)

F ( 2)

3 ? wL ? ? 25 ×10 × 2.5 ? ? ? ? 2 ? ? 2 ? ? ? ? 2 3 2 ? ? wL ? ? 25 × 10 × 2.5 ? ? ?31 250 ? ? ? ?13 021? ? 12 ? ? 12 ? =? ?=? ?=? wL ? ? 25 × 103 × 2.5 ? ? ?31 250 ? ?? ? ? 13 021 ? ? 2 ? ? 2 ? ? ? ? ? ? wL2 ? ? 25 × 103 × 2.52 ? ? ? 12 ? ? ? ? ? ? 12 ?

(3-147)

合并上述两个荷载列阵,形成总载荷列阵,有

F = F (1) + F (2)

?62 500 ? ? ? ?62 500 ? ? ? ? ?52 083 ? ?52 083 ? ? ? ? ? ?62 200 ? 31 250 ? ? ?93 750 ? =? ?=? ? ? ?52 083 ? 13 021? ? 39 062 ? ? ? ? ?31 250 ? ? 31 250 ? ? ? ? 13 021 ? ? ? 13 021 ? ? ? ? ?

(3-148)

这样可以得到该结构的整体刚度方程为

Kq = P = F + R

(3-149)

考虑位移边界条件:在节点 1 处有 v1 = θ1 = 0 ,在节点 2 处有 v2 = 0 ,将上述刚度方程化为

? 56 958 000 ?22 771 200 18 976 000 ? ?θ 2 ? ? 39 062 ? ? ?22 771 200 18 216 960 ?22 771 200 ? ? v ? = ? ?31 250 ? ? ? ?? 3? ? ? 18 976 000 ?22 771 200 37 952 000 ? ?θ 3 ? ? 13 021 ? ? ? ?? ? ?
求解上述方程组,则可求得各节点的位移值

(3-150)

q = [ v1 θ1 v2 θ 2 v3 θ3 ]T = ? 0 0 0 ?0.001 3723 ( rad ) ?0.008 577 2 ( m ) ?0.004 117 ( rad ) ? ? ?
然后,根据下述关系求解得各节点反力和弯矩
T

(3-151)

R = Kq ? F
61

(3-152)

有限元分析基础教程

曾攀

其中 R 为反力矩阵。将相关值代入式(3-152),有

R = [ Rv1 Rθ 1 Rv 2 0 0 0]T = [54687.8(N) 39062.6(Nm) 132814.8(N) 0 0 0]T

(3-153)

3.3.3 平面梁单元的坐标变换
图 3-17 所示为一整体坐标系中的平面梁单元,它有两个端节点,梁的长度为 l,弹性模量 为 E,横截面的面积为 A,惯性矩为 Iz。

图 3-17 平面问题中梁单元的坐标变换

设局部坐标系下 ( Oxy ) 的节点位移列阵为

(6×1)

q e = [u1 v1 θ1 u 2

v2 θ 2 ]

T

(3-154)

整体坐标系中 ( Oxy ) 的节点位移列阵为

(6×1)

q e = u 1 v1 θ 1 u 2

[

v2 θ2

]

T

(3-155)

注意:转角 θ1 和 θ 2 在两个坐标系中是相同的。按照两个坐标系中的位移向量相等效的原则,可 推导出以下变换关系。

u1 = u1 cos α + v 1 sin α v1 = ?u1 sin α + v 1 cos α u2 = u 2 cos α + v 2 sin α v2 = ?u 2 sin α + v 2 cos α
写成矩阵形式有

? ? ? ? ? ? ?

(3-156)

(6×1) (6×6) (6×1)

q e = Te ? q e

(3-157)

62

有限元分析基础教程

曾攀

其中 T 为单元的坐标变换矩阵,即
e

? cos α ? ? sin α ? ? 0 e T =? (6×6) ? 0 ? 0 ? ? 0 ?

sin α cos α 0 0 0 0
e

0 0 1 0 0 0

0 0 0 cos α ? sin α 0
e e

0 0 0 sin α cos α 0

0 0 0 0 0 1

? ? ? ? ? ? ? ? ? ?

(3-158)

与平面杆单元的坐标变换类似,梁单元在整体坐标系中的刚度方程为
(6×6) 6×1) (6×1) (

K ?q = P

(3-159)

其中
(6×6) (6×6) 6×6) 6×6) ( (

K = TeT ? K e ? Te
e

(3-160) (3-161)

(6×1) (6×6) 6×1) (

P = TeT P e

e

3.3.4 空间梁单元及坐标变换
1. 空间梁单元 空间梁单元除承受轴力和弯矩外,还可能承受扭矩的作用,而且弯矩可能同时在两个坐标 面内存在. 图 3-18 所示为一局部坐标系中的空间梁单元,其长度为 l,弹性模量为 E,横截面的 惯性矩为 I z (绕并行于 z 轴的中性轴)和 I y (绕并行于 y 轴的中性轴),横截面的扭转惯性矩为 J。 对应于图 3-18 中梁单元,设有两个端节点,每一个节点的位移自由度有 6 个,单元共有 12 个自由度;其局部坐标系中的节点位移列阵 q 和节点力列阵 P 如下。
e
e

图 3-18 局部坐标系中的空间梁单元

q e = ? u1 ? (12×1)

v1

w1

θ x1 θ y1 θ z1

u2

v2

w2

θ x2 θ y2 θ z2 ? ?

T

(3-162)

63

有限元分析基础教程

曾攀

P e = ? Pu1 (12×1) ?

Pv1

Pw1

M x1

M y1

M z1

Pu 2

Pv 2

Pw 2

M x2

M y2

M z2 ? ?

T

(3-163)

下面,基于前面杆单元和平面梁单元的刚度矩阵分别写出对应于图 3-18 中各对应节点位移 的刚度矩阵,然后进行组合以形成完整的刚度矩阵。 (1) 对应于图 3-18 中的节点位移(u1,u2) 这是轴向位移,由式(3-160),有对应于杆单元的刚度矩阵为
e K u1u2 =

(2× 2)

EA ? 1 ?1? l ? ?1 1 ? ? ?

(3-164)

(2) 对应于图 3-18 中的节点位移( θ x1 , θ x 2 ) 这是杆受扭转时的情形,如果将扭转角类似于拉伸杆的轴向位移,则它的分析结果与拉伸 杆类似(见材料力学中的扭转问题),所推导出的刚度矩阵和(3-164)式相似,即
e Kθ x1θ x 2 =

(2× 2)

GJ l

? 1 ?1? ? ?1 1 ? ? ?

(3-165)

其中 J 为横截面的扭转惯性矩,G 为剪切模量。 (3) 对应于图 3-18 中 Oxy 平面内的节点位移(v1, θ z1 ,v2, θ z 2 ) 这是梁在 Oxy 平面内的纯弯曲情形,有对应的刚度矩阵

K eOxy ) (
(4× 4)

6l ? 12 ? 6l 4l 2 EI = 3Z ? l ? ?12 ?6l ? 2 ? 6l 2l

?12 6l ? ?6l 2l 2 ? ? 12 ?6l ? ? ?6l 4l 2 ?

(3-166)

其中 I z 为绕并行于 z 轴的中性轴的惯性矩。 (4) 对应于图 3-18 中 Oxz 平面内的节点位移(w1, θ y1 ,w2, θ y 2 ) 这是梁在 Oxz 平面内的纯弯曲情形,可得到与(3-166)式类似的刚度矩阵,但所对应的节点 位移是不同的。 (6) 将各分刚度矩阵进行组合以形成完整的单元刚度矩阵 对应于式(3-162)中的节点位移的次序,分别将上面各个部分的刚度阵的元素进行组合,可 形成局部坐标系中空间梁单元的完整刚度矩阵,即

64

有限元分析基础教程
u1 ? EA ? ? l ? ? 0 ? ? ? 0 ? ? ? 0 ? ? 0 ? ? ? 0 ? e K =? 12 ( ×12) ? ? EA ? l ? ? 0 ? ? ? 0 ? ? ? 0 ? ? ? 0 ? ? ? 0 ? ? ? v1 0 12EI z l3 0 0 0 6EI z l2 0 ? 12EI z l3 0 0 0 6EI z l2 ? ? ? w1 0 0 12EI y l3 0 6EI y l2 0 0 0 12EI y l3 0 6EI y l 0
2

曾攀
v2 w2 0 0 ? 12EI y l3 0 6EI y l2 0 0 0 12EI y l3 0 6EI y l
2

θx1
0 0 0 GJ l 0 0 0 0 0 ? GJ l 0 0 ?

θy1
0 0 6EI y l2 0 4EI y l 0 0 0 6EI y l2 0 2EI y l 0 ?

θz1
0 6EI z l2 0 0 0 4EI z l 0 6EI z l2 0 0 0 2EI z l ?

u2 EA l 0 0 0 0 0 EA l 0 0 0 0 0 ? ? ?

θx2
0 0 0 ? GJ l 0 0 0 0 0 GJ l 0 0 ?

θy2
0 0 6EI y l2 0 2EI y l 0 0 0 6EI y l2 0 4EI y l 0 ?

θz 2
0 6EI z l2 0 0 0 2EI z l 0 6EI z l2 0 0 0 4EI z l ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

0 12EIz l3 0 0 0 6EI z l2 0 12EI z l3 0 0 0 6EI z l2

0

(3-167)
2. 空间梁单元的坐标变换 空间梁单元坐标变换的原理和方法与平面梁单元的坐标变换相同,只要分别写出两个坐标 系中的位移向量的等效关系则可得到坐标变换矩阵,即 局部坐标系中空间梁单元的节点位移列阵为
(12×1)

q e = u1 v1

[

w1 θ x1 θ y1 θ z1 u 2

v2

w2 θ x 2 θ y 2 θ z 2

]

T

(3-168)

整体坐标系中的节点位移列阵为
(12×1)

q e = u1

[

v1

w1 θ x1 θ y1 θ z1

u2

v2

w2 θ x2 θ y 2 θ z 2

]

T

(3-169)

对应于式(3-168)中的各组位移分量,可分别推导相应的转换关系,具体地,对于端节点 1,有

? u1 ? ? u1 ? ? u1 cos( x, x ) + v1 cos( x, y ) + w1 cos( x, z ) ? ? v ? = ?u cos( y, x) + v cos( y, y ) + w cos( y, z ) ? = λ ? ? v ? 1 1 ? (3×3) ? 1 ? ? 1? ? 1 ? w1 ? ? u1 cos( z , x) + v1 cos( z , y ) + w1 cos( z , z ) ? ? w1 ? ? ? ? ? ? ?

(3-170)

?θ x1 ? ?θ x1 ? ? θ x1 cos( x, x) + θ y1 cos( x, y ) + θ z1 cos( x, z ) ? ? ? ?θ ? = θ cos( y, x) + θ cos( y, y ) + θ cos( y, z ) = λ ? ?θ ? y1 z1 ? (3×3) ? y1 ? ? y1 ? ? x1 ? ? ? ? ?θ z1 ? ? ? ? θ x1 cos( z , x) + θ y1 cos( z , y ) + θ z1 cos( z , z ) ? ?θ z1 ? ? ? ? ?

(3-171)

65

有限元分析基础教程

曾攀

同样,对于端节点 2,有以下转换关系

?u 2 ? ? u2 ? ?v ? = λ ? ?v2 ? ? 2 ? (3×3) ? ? ?w 2 ? ? w2 ? ? ? ? ?
?θ x 2 ? ?θ x 2 ? ?θ ? = λ ? ?θ y 2 ? ? y 2 ? (3×3) ? ? ?θ z 2 ? ?θ z 2 ? ? ? ? ?
以上的 λ 为节点坐标变换矩阵,即

(3-172)

(3-173)

? cos( x, x) cos( x, y ) cos( x, z ) ? ? ? λ =? cos( y, x) cos( y , y ) cos( y, z ) ? (3× 3) ? cos( z , x) cos( z , y ) cos( z , z ) ? ? ?

(3-174)

其中 cos( x, x ) , ……, cos( z , z ) 分别表示局部坐标轴(x,y,z)对整体坐标轴 x, y, z 的方向 余弦。 将式(3-170)~式(3-173)写在一起,有
( (12×1) (12×12) 12×1)

(

)

q e = Te ? q

e

(3-175)

其中 T 为坐标变换矩阵
e

? (3λ3) ? × ? 0 (3×3) e T =? (12×12) ? 0 ? (3×3) ? 0 ? (3×3) ?

(3×3) (3×3) (3×3) (3×3)

0 λ

(3×3) (3×3) (3×3) (3×3)

0

0

0 0

λ

0

0 ? ? 0 ? (3×3) ? 0 ? (3×3) ? λ ? (3×3) ? ?
(3×3)

(3-176)

有了坐标变换矩阵,就很容易写出整体坐标系下的刚度矩阵和刚度方程。

3.3.5 梁单元的常用等效节点载荷
表 3-4 列出了常用的梁单元在承受非节点载荷下的节点载荷等效值,该等效值是根据外力功 的计算公式得到的,因此,它与梁单元的边界条件没有关系(表 3-4 中的图示虽为固支,这些节 点载荷等效值也可以用在其它边界情况)。
表 3-4 支承与外载情况 梁单元的常见等效节点荷载 所等效节点荷载

66

有限元分析基础教程

曾攀

W e = ∫z p ( x ) ? v ( x ) dx
等效原理

= [ RA

? ? = ? ∫l p ( x ) ? N ( x ) dx ? MA RB

M B ] qe

RA = ? P / 2 RB = ? P / 2 M A = ? PL / 8 M B = PL / 8
RA = ? ( Pb 2 / L3 ) ( 3a + b ) RB = ? ( Pa 2 / L3 ) ( a + 3b )

M A = ? Pab 2 / L2 M B = Pa 2 b / L2 RA = ? p0 L / 2 RB = ? p0 L / 2 M A = ? p0 L2 /12 M B = p0 L2 /12 RA = ?3 p0 L / 20 RB = ?7 p0 L / 20 M A = ? p0 L2 / 30 M B = p0 L2 / 20
RB = ? ( p0 a 3 / 2 L3 ) ( 2 L ? a )
RA = ? ( p0 a / 2 L3 ) ( a 3 ? 2a 2 L + 2 L3 )

M A = ? ( p0 a 2 /12 L2 ) ( 3a 2 ? 8aL + 6 L2 )
M B = ( p0 a 3 /12 L2 ) ( 4 L ? 3a )

RA = ? p0 L / 4 RB = ? p0 L / 4 M A = ?5 p0 L2 / 96 M B = 5 p0 L2 / 96

67

有限元分析基础教程

曾攀

RA = ?6 M 0 ab / L3
M A = ? ( M 0b / L2 ) ( 3a ? L )

RB = 6M 0 ab / L3

M B = ? M 0 a / L2 (3b ? L )

(

)

3.3.6 梁单元分析的 MATLAB 程序
【MATLAB 程序】3.3.6(1) 1D 梁单元的有限元分析程序(Beam1D2Node) 编写如图 3-12 所示局部坐标系中的梁单元的刚度矩阵、单元组装、单元应力的计算程序。 解答:编写的 2 节点梁单元的五个 MATLAB 函数如下。
Beam1D2Node_Stiffness(E,I,L) 该函数计算单元的刚度矩阵,输入弹性模量 E,横截面的惯性矩 I,梁单元的长度 L,输出单元刚度矩阵 k(4×4)。 Beam1D2Node _Assembly(KK,k,i,j) 该函数进行单元刚度矩阵的组装,输入单元刚度矩阵 k,单元的节点编号 i、j、m,输出整体刚度矩阵 KK。 Beam1D2Node_ Strain(x,L,y) 该函数计算单元的几何矩阵,输入所测点距梁单元左节点的水平距离 x,输入所测点以中性层为起点的 y 方向的坐标,梁单元的长度 L,输出单元几何形状函数矩阵 B(1×4) 。 Beam1D2Node _Stress(E,B,u) 该函数计算单元内某点的应力,输入弹性模量 E,几何矩阵 B,节点位移列阵 u,输出单元的应力 stress Beam1D2Node_Deflection(x,L,u) 该函数计算单元内某点的挠度,输入所测点距梁单元左节点的水平距离 x,梁单元的长度 L,节点位移列 阵 u,输出该点的挠度 v。

基于第 3.3.2 节中的基本公式,可以编写出具体实现以上每个函数的 MATLAB 程序如下。
%%%%%%%%%% Beam1D2Node %%% begin %%%%%%%%%%%% function k =Beam1D2Node_Stiffness(E,I,L) %该函数计算单元的刚度矩阵 %输入弹性模量 E,横截面的惯性矩 I,梁单元的长度 L %输出单元刚度矩阵 k(4×4) %----------------------------------------k = E*I/(L*L*L)*[12 6*L -12 6*L ; 6*L 4*L*L -6*L 2*L*L ; -12 -6*L 12 -6*L ; 6*L 2*L*L -6*L 4*L*L]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function z = Beam1D2Node_Assembly(KK,k,i,j) %该函数进行单元刚度矩阵的组装 %输入单元刚度矩阵 k,单元的节点编号 i、j、m %输出整体刚度矩阵 KK %-----------------------------------------

68

有限元分析基础教程 DOF(1)=2*i-1; DOF(2)=2*i; DOF(3)=2*j-1; DOF(4)=2*j; for n1=1:4 for n2=1:4 KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2); end end z=KK; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function B= Beam1D2Node_Strain(x,L,y) %该函数计算单元的几何矩阵 %输入所测点距梁单元左节点的水平距离 x %输入所测点以中性层为起点的 y 方向的坐标,梁单元的长度 L %输出单元几何形状函数矩阵 B(1×4) %----------------------------------------e=x/L; B1=(12*e-6)/(L*L); B2=(6*e-4)/L; B3=-(12*e-6)/(L*L); B4=(6*e-2)/L; B=-y*[B1,B2,B3,B4]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function stress= Beam1D2Node_Stress(E,B,u) %该函数计算单元内某点的应力 %输入弹性模量 E,几何矩阵 B,节点位移列阵 u %输出单元的应力 stress %----------------------------------------stress = E*B*u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function v=Beam1D2Node _Deflection(x,L,u) %该函数计算单元内某点的挠度 %输入所测点距梁单元左节点的水平距离 x %输入梁单元的长度 L,节点位移列阵 u %输出该点的挠度 v %----------------------------------------e=x/L; N1=1-3*e*e+2*e*e*e; N2=L(e-2*e*e+e*e*e); N3=3*e*e-2*e*e*e; N4=L(e*e*e-e*e); N=[N1,N2,N3,N4]; v=N*u; %%%%%%%%%% Beam1D2Node %%% end %%%%%%%%%%%%

曾攀

【MATLAB 程序】3.3.6(2) 2D 梁单元的有限元分析程序(Beam2D2Node) 编写如图 3-13 所示平面梁单元的单元刚度矩阵、单元组装、单元应力的计算程序。 解答: 编写的平面 2 节点梁单元的三个 MATLAB 函数如下。
Beam2D2Node_Stiffness(E,I,A,L)

69

有限元分析基础教程

曾攀

以上函数计算单元的刚度矩阵,输入弹性模量 E,横截面积 A,惯性矩 I,长度 L,输出单元刚度矩阵 k(6 ×6)。 Beam2D2Node_Assemble(KK,k,i,j) 以上函数进行单元刚度矩阵的组装,输入单元刚度矩阵 k,单元的节点编号 i、j,输出整体刚度矩阵 KK。 Beam2D2Node_Forces(k,u) 以上函数计算单元的应力,输入单元刚度矩阵 k,节点位移 u,输出单元节点力 forces。

基于第 3.3.3 节的基本公式,可以编写出具体实现以上每个函数的 MATLAB 程序如下。
%%%%%%%%%% Beam2D2Node %% begin %%%%%%%%%%%% function k =Beam2D2Node_Stiffness(E,I,A,L) %以上函数计算单元的刚度矩阵 %输入弹性模量 E,横截面积 A,惯性矩 I,长度 L %输出单元刚度矩阵 k(6×6) %----------------------------------------k=[E*A/L,0,0,-E*A/L,0,0;0,12*E*I/(L^3),6*E*I/(L^2),0,-12*E*I/(L^3),6*E*I/(L^2);0,6*E*I/(L^2),4*E*I/L,0,-6*E*I/ (L^2),2*E*I/L;-E*A/L,0,0,E*A/L,0,0;0,-12*E*I/(L^3),-6*E*I/(L^2),0,12*E*I/(L^3),-6*E*I/(L^2);0,6*E*I/(L^2),2* E*I/L,0,-6*E*I/(L^2),4*E*I/L] %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y =Beam2D2Node_Assemble(KK,k,i,j) %以上函数进行单元刚度矩阵的组装 %输入单元刚度矩阵 k,单元的节点编号 i、j %输出整体刚度矩阵 KK %----------------------------------------DOF(1)=3*i-2; DOF(2)=3*i-1; DOF(3)=3*i; DOF(4)=3*j-2; DOF(5)=3*j-1; DOF(6)=3*j; for n1=1:6 for n2=1:6 KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2); end end y = KK; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function forces = Beam2D2Node_Forces(k,u) %以上函数计算单元的应力,输入单元刚度矩阵 k,节点位移 u, %输出单元节点力 forces %----------------------------------------forces = k * u; %%%%%%%%%% Beam2D2Node %% end %%%%%%%%%%%%

3.3.7 梁结构分析的算例

70

有限元分析基础教程

曾攀

【典型例题】3.3.7(1) 三梁平面框架结构的有限元分析 如图 3-19 所示的框架结构,其顶端受均布力作用,用有限元方法分析该结构的位移。结构 中各个截面的参数都为: E = 3.0 × 1011 Pa , I = 6.5 × 10 ?7 m 4 , A = 6.8 × 10 ?4 m 2 。

图 3-19 框架结构受一均布力作用

(a)节点位移及单元编号

(b)等效在节点上的外力

图 3-20 单元划分、节点位移及节点上的外载

解答:对该问题进行有限元分析的过程如下。 (1) 结构的离散化与编号 将该结构离散为三个单元,节点位移及单元编号如图 3-20(a)所示,有关节点和单元的信息 见表 3-5。
表 3-5 单元编号与节点编号

单元编号 ① ② ③ 1 3 4

节点编号 2 1 2

71

有限元分析基础教程

曾攀

节点位移列阵为

q = [u1

v1 θ 1

u2

v2 θ 2

u3

v3 θ 3

u4

v4 θ 4 ]

T

(3-177)

节点外载列阵为(均布外载的等效计算见表 4.1)

F = Fx1
支反力列阵为

[

F y1

Mθ1

0 Fy 2

Mθ 2

0 0 0 0 0 0

]

T

(3-178)

R = 0 0 0 0 0 0 Rx3

[

R y3

Rθ 3

Rx 4

Ry4

Rθ 4

]

T

(3-179)

其中: Rx 3 、 Ry 3 、 Rθ 3 为节点 3 的沿 x 方向支反力、沿 y 方向支反力和支反力矩, Rx 4 、 Ry 4 、

Rθ 4 为节点 4 的沿 x 方向支反力、沿 y 方向支反力和支反力矩,均为待求值。
总的节点载荷列阵为

P =F+R T = ?3000 ?3000 ?720 0 ?3000 720 Rx3 Ry3 Rθ 3 Rx 4 Ry 4 Rθ 4 ? ? ?
(2) 各个单元的描述 单元①的局部坐标与整体坐标是一致的,则可以由式(3-136)直接得到

(3-180)

K

(1)

u1 v1 θ1 u2 v2 θ2 ↓ ↓ ↓ ↓ ↓ ↓ ?141.7 0 0 0 0 ? ← u1 ? 141.7 ? 0 ?0.784 0.564 ? ← v1 0.784 0.564 0 ? ? ? 0 ?0.564 0.271 ? ← θ1 0.564 0.542 0 = 106 × ? ? 0 0 141.7 0 0 ? ← u2 ? ?141.7 ? 0 ?0.784 ?0.564 0 0.784 ?0.564 ? ← v2 ? ? 0.564 0.271 0 ?0.564 0.542 ? ← θ 2 ? 0 ? ?

(3-181)

单元②和单元③的情况相同,只是节点编号不同而已,其局部坐标系下的单元刚度矩阵为

? K (2)

0 0 ?212.5 0 0 ? ? 212.5 ? 0 2.645 1.270 0 ?2.645 1.270 ? ? ? ? 0 1.270 0.8125 0 ?1.270 0.4062 ? = 106 × ? ? 0 0 212.5 0 0 ? ? ?212.5 ? 0 ?2.645 ?1.270 0 2.645 ?1.270 ? ? ? 1.270 0.4062 0 ?1.270 0.8125 ? ? 0 ? ?

(3-182)

这两个单元轴线的方向余弦为 cos( x, x ) = 0 , cos( x, y ) = 1 ,有坐标转换矩阵

72

有限元分析基础教程

曾攀

?0 ? ?1 ? ?0 T=? ?0 ?0 ? ?0 ?

1 0 0 0 0 1 0 0 0 0 0 0

0 0? 0 0 0? ? 0 0 0? ? 0 1 0? ?1 0 0 ? ? 0 0 1? ? 0

(3-183)

则可以计算出整体坐标下的单元刚度矩阵(单元②和单元③)

K (2)

0 ?1.27 ?2.645 0 ?1.27 ? ? 2.645 ? 0 212.5 0 0 0 ? ?212.5 ? ? ? ?1.27 0 0.8125 1.27 0 0.4062 ? ? = TT ? K (2) ? T = 106 × ? ? 0 1.27 2.645 0 1.27 ? ? ?2.645 ? 0 0 0 212.5 0 ? ?212.5 ? ? 0 0.4062 1.27 0 0.8125? ? ?1.27 ? ?
对于单元②: [u 3 对于单元③: [u 4

(3-184)

注意这两个单元所对应的节点位移列阵分别为

v3 θ 3
v4 θ 4

u1
u2

v1 θ 1 ]

T

(3-185) (3-186)

v2 θ 2 ]

T

(3) 建立整体刚度方程 组装整体刚度矩阵并形成整体刚度方程

K ?q = P
其中刚度矩阵的组装关系为

(3-187)

K = K (1) + K (2) + K (3)
(4) 边界条件的处理及刚度方程求解 该问题的位移边界条件为

(3-188)

u 3 = v3 = θ 3 = u 4 = v 4 = θ 4 = 0
处理该边界条件后的刚度方程为

(3-189)

0 1.270 ?141.7 0 0 ? ? u1 ? ? 3000 ? ? 144.3 ? 0 213.3 0.564 0 ?0.784 0.564 ? ? v1 ? ? ?3000 ? ? ? ?? ? ? ? 1.270 0.564 1.3545 0 ?0.564 0.271 ? ?θ1 ? ? ?720 ? 106 × ? ? ?? ? = ? 0 0 144.3 0 1.270 ? ?u2 ? ? 0 ? ? ?141.7 ? 0 0 213.3 ?0.564 ? ? v2 ? ? ?3000 ? ?0.784 ?0.564 ? ? ?? ? ? 0.564 0.271 1.270 ?0.564 1.3545 ? ?θ 2 ? ? 720 ? ? 0 ? ? ?? ? ?

(3-190)

73

有限元分析基础教程

曾攀

求解后的结果为

u1 = 0.92mm v1 = ?0.0104mm θ1 = ?0.00139rad u2 = 0.901mm v2 = ?0.018mm θ 2 = 3.88 × 10?5 rad

? ? ? ? ? ? ?

(3-191)

【MATLAB 算例】3.3.7(2)

三梁平面框架结构的有限元分析(Beam2D2Node)

参见【典型例题】3.3.7(1)所给算例,对于如图 3-19 所示的结构,基于 MATLAB 平台求解 该结构的节点位移以及支反力。 解答:对该问题进行有限元分析的过程如下。 (1) 结构的离散化与编号 将该结构离散为三个单元,节点位移及单元编号如图 3-20 所示,有关节点和单元的信息见 表 3-5。 (2) 各个单元的描述 首先在 MATLAB 环境下,输入弹性模量 E、横截面积 A、惯性矩 I、长度 L,然后针对单 元 1 和单元 2 单元 3,分别二次调用函数 Beam2D2Node_Stiffness,就可以得到单元的刚度矩阵 k1(6×6)和 k2(6×6),且单元 2 和单元 3 的刚度矩阵相同。
>> E=3E11; >> I=6.5E-7; >> A=6.8E-4; >> L1=1.44; >> L2=0.96; >> k1=Beam2D2Node_Stiffness(E,I,A,L1); >> k2=Beam2D2Node_Stiffness(E,I,A,L2);

(3) 建立整体刚度方程 将单元 2 和单元 3 的刚度矩阵转换成整体坐标下的形式。由于该结构共有 4 个节点,则总 共的自由度数为 12,因此,结构总的刚度矩阵为 KK(12×12),对 KK 清零,然后两次调用函数 Beam2D2Node_Assemble 进行刚度矩阵的组装。
>> T=[0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1]; >> k3=T'*k2*T; >> KK=zeros(12,12); >> KK=Beam2D2Node_Assemble(KK,k1,1,2); >> KK=Beam2D2Node_Assemble(KK,k3,3,1); >> KK=Beam2D2Node_Assemble(KK,k3,4,2) KK = 1.0e+008 * 1.4431 0 0.0127 -1.4167 0 0 -0.0264 0 0.0127 0 0 2.1328 0.0056 0 -0.0078 0.0056 0 -2.1250 0 0 0.0127 0.0056 0.0135 0 -0.0056 0.0027 -0.0127 0 0.0041 0 -1.4167 0 0 1.4431 0 0.0127 0 0 0 -0.0264

0 0 0 0

0 0 0 0.0127

74

有限元分析基础教程 0 -0.0078 -0.0056 0 2.1328 -0.0056 0 0 0 0.0056 0.0027 0.0127 -0.0056 0.0135 0 0 -0.0264 0 -0.0127 0 0 0 0.0264 0 0 -2.1250 0 0 0 0 0 2.1250 0.0127 0 0.0041 0 0 0 -0.0127 0 0 0 0 -0.0264 0 -0.0127 0 0 0 0 0 0 -2.1250 0 0 0 0 0 0 0.0127 0 0.0041 0 0 0 0 -2.1250 0 0 -0.0127 0 0.0041 -0.0127 0 0 0 0 0 0 0 0.0081 0 0 0 0 0.0264 0 -0.0127 0 0 2.1250 0 0 -0.0127 0 0.0081

曾攀

(4) 边界条件的处理及刚度方程求解 该问题的位移边界条件为 u 3 = v3 = θ 3 = u 4 = v 4 = θ 4 = 0 。因此,将针对节点 1 和节点 2 的位移进行求解,节点 1 和节点 2 的位移将对应 KK 矩阵中的前 6 行和前 6 列,则需从 KK(12× 12)中提出,置给 k,然后生成对应的载荷列阵 p,再采用高斯消去法进行求解,注意:MATLAB 中的反斜线符号“\”就是采用高斯消去法。
>> k=KK(1:6,1:6); >> p=[3000;-3000;-720;0;-3000;720]; >> u=k\p u = 0.0009 -0.0000 -0.0014

0.0009

-0.0000

-0.0000

[将列排成了行]

(5)支反力的计算 在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力;先将上面 得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系),可以得到整体的位移列 阵 U(12×1),再代回原整体刚度方程,计算出所有的节点力 P(12×1),按式(3-180)的对应关系 就可以找到对应的支反力。
>> U=[u;0;0;0;0;0;0] U = 0.0009 -0.0000 0 0 >> P=KK*U P = 1.0e+003 * 3.0000 -3.0000 -0.6658 2.2012 -0.0014 0 0.0009 0 -0.0000 0 -0.0000 0 [将列排成了行] [将列排成了行]

-0.7200 0.6014

0.0000 -2.3342

-3.0000 0.7200 3.7988 1.1283

[将列排成了行] [将列排成了行]

【ANSYS 算例】3.3.7(3) 三梁平面框架结构的有限元分析 针对【典型例题】3.3.7(1)的问题,在 ANSYS 平台上,完成相应的力学分析。 解答:对该问题进行有限元分析的过程如下。 1 基于图形界面(GUI)的交互式操作(step by step)
(1) 进入 ANSYS(设定工作目录和工作文件) 程序 →ANSYS → ANSYS Interactive →Working directory (设置工作目录) →Initial jobname(设置工作 文件名): beam3→Run → OK (2) 设置计算类型 ANSYS Main Menu: Preferences… → Structural → OK (3) 选择单元类型 ANSYS Main Menu: Preprocessor →Element Type→Add/Edit/Delete… →Add… →beam:2D elastic 3

75

有限元分析基础教程 →OK (返回到 Element Types 窗口) →Close (4) 定义材料参数

曾攀

ANSYS Main Menu: Preprocessor →Material Props →Material Models→Structural →Linear →Elastic → Isotropic: EX:3e11 (弹性模量) → OK → 鼠标点击该窗口右上角的“ ”来关闭该窗口

(5) 定义实常数以确定平面问题的厚度 ANSYS Main Menu: Preprocessor →Real Constants… →Add/Edit/Delete →Add →Type 1 Beam3→ OK →Real Constant Set No: 1 (第 1 号实常数), Cross-sectional area:6.8e-4 (梁的横截面积) →OK →Close (6) 生成几何模型 生成节点 ANSYS Main Menu: Preprocessor →Modeling →Creat→Nodes→ In Active CS→Node number 1 → X:0,Y:0.96,Z:0 →Apply→Node number 2 → X:1.44,Y:0.96,Z:0 →Apply→Node number 3 → X:0,Y:0,Z:0→ Apply→Node number 4 → X:1.44,Y:0,Z:0→OK 生成单元 ANSYS Main Menu: Preprocessor → Modeling 成单元 3)→OK (7) 模型施加约束和外载 左边加 X 方向的受力 ANSYS Main Menu: Solution → Define Loads → Apply →Structural → Force/Moment → On Nodes → 选择节点 1→ apply →Direction of force: FX →VALUE:3000 → OK→ 上方施加 Y 方向的均布载荷 ANSYS Main Menu: Solution → Define Loads → Apply →Structural → Pressure →On Beams →选取 单元 1(节点 1 和节点 2 之间)→ apply →VALI:4167→VALJ:4167→OK 左、右下角节点加约束 ANSYS Main Menu: Solution → Define Loads → Apply →Structural → Displacement → On Nodes → 选取节点 3 和节点 4 → Apply → Lab:ALL DOF → OK (8) 分析计算 ANSYS Main Menu: Solution → Solve → Current LS → OK → Should The Solve Command be Executed? Y→ Close (Solution is done! ) → 关闭文字窗口 (9) 结果显示 ANSYS Main Menu: General Postproc → Plot Results →Deformed Shape … → Def + Undeformed → OK (返回到 Plot Results) (10) 退出系统 ANSYS Utility Menu: File→ Exit …→ Save Everything→OK (11) 计算结果的验证 与 MATLAB 支反力计算结果一致。

→ Create →Element → Auto Numbered → Thru

Nodes → 选择节点 1、2(生成单元 1)→ apply → 选择节点 1、 3(生成单元 2)→ apply →选择节点 2、 4(生

2 完全的命令流

76

有限元分析基础教程 !%%%%%%%%%% / PREP7 ET,1,beam3 R,1,6.5e-7,6.8e-4 MP,EX,1,3e11 N,1,0,0.96,0 N,2,1.44,0.96,0 N,3,0,0,0 N,4,1.44,0,0 E,1,2 E,1,3 E,2,4 D,3,ALL D,4,ALL F,1,FX,3000 SFBEAM,1,1,PRESS,4167 FINISH /SOLU SOLVE FINISH /POST1 PLDISP,1 FINISH !%%%%%%%%%% !将 3 号节点的位移全部固定 !将 4 号节点的位移全部固定 !在 1 号节点处施加 x 方向的力(3000) !施加均布压力 !结束前处理状态 !进入求解模块 !求解 !结束求解状态 !进入后处理 !显示变形状况 !结束后处理 [典型例题]3_3_7(3) %%% end %%%%% !生成单元(连接 1 号节点和 2 号节点) ,以下类似 [典型例题]3_3_7(3) !进入前处理 !选择单元类型 !给出实常数(横截面积、惯性矩) !给出材料的弹性模量 !生成 4 个节点,坐标(0,0.96,0),以下类似 %%% begin %%%%%

曾攀

3.4 应用:桥梁结构的 ANSYS 参数化分析
3.4.1 桥梁结构描述 桥梁是现代建筑中最重要的结构形式之一,它对一个区域的政治、经济、文化有着重要的 影响,现代桥梁正向着长距离、大跨度方向发展;目前的大跨度桥梁主要有斜拉桥以及悬索桥 两种形式;例如,法国于 1995 年建成的 Normandie 大桥就是大跨度斜拉桥,有 184 根拉索,主 跨长度达到 864m,桥面高度为 51m,使用了 1 万 9 千吨钢材,7 万立方米的混凝土,见图 3-21。

(a) 法国 Normandie 大桥全景

(b) 斜拉桥的施工过程

图 3-21 现代斜拉桥的结构

现代桥梁设计的一个重要特点就是需要采用精确的有限元方法对桥梁的静力、 振动、 风载、 地震、四季温差等状况进行分析,并提出优化方案,并且还要就桥梁施工过程中因结构变化的

77

有限元分析基础教程

曾攀

各种受力状态进行细致的分析,从而确定出可行的施工方案,可以毫不夸张地说,有限元分析 在现代结构的设计和分析中正发挥着不可替代的作用。 下面以一个简单桁架桥梁为例,以展示有限元分析的全过程;背景素材选自位于密执安的 "Old North Park Bridge" (1904-1988),见图 3-22;该桁架桥由型钢组成,顶梁及侧梁,桥身弦杆, 底梁分别采用三种不同型号的型钢,结构参数见表 3-6。桥长 L=32m,桥高 H=5.5m。桥身由 8 段桁架组成,每段长 4m。该桥梁可以通行卡车,若这里仅考虑卡车位于桥梁中间位置,假设卡 车的质量为 4000kg,若取一半的模型,可以将卡车对桥梁的作用力简化为 P1 ,P2 ,P3 ,其 中 P1= P3=5000 N, P2=10000N,见图 3-23。

图 3-22 位于密执安的"Old North Park Bridge" (1904-1988)

图 3-23 桥梁的简化平面模型(取桥梁的一半) 表 3-6 桥梁结构中各种构件的几何性能参数 构件 惯性矩 横截面积

顶梁及侧梁(beam1)

3.83 ×10?6 m 4 1.87 ×10?6 m 4

2.19 ×10?3 m 2 1.185 ×10?3 m 2

桥身弦梁(beam2)

底梁(beam3)

8.47 ×10?6 m 4

3.031×10?3 m 2

下面分别给出该桥梁的三种分析方法:基于图形界面(GUI)的交互式操作、基于命令流方式 以及基于参数化方式,这三种方式都是在 ANSYS 平台上完成的。 3.4.2 基于 ANSYS 的桁架桥梁结构分析

78

有限元分析基础教程

曾攀

【ANSYS 算例】3.4.2(1) 基于图形界面(GUI)的桁架桥梁结构分析 解答:以下为基于 ANSYS 图形界面(GUI, Graphic User Interface)的菜单操作流程。
(1) 进入 ANSYS(设定工作目录和工作文件)   程序 → ANSYS → ANSYS Interactive → Working directory(设置工作目录)→ Initial jobname(设 置工作文件名):TrussBridge → Run → OK (2) 设置计算类型    ANSYS Main Menu:Preferences… → Structural → OK (3) 定义单元类型 ANSYS Main Menu: Preprocessor → Element Type → Add/Edit/Delete... → Add…→ Beam: 2d elastic 3 → OK(返回到 Element Types 窗口)→ Close (4) 定义实常数以确定梁单元的截面参数 ANSYS Main Menu: Preprocessor → Real Constants…→ Add/Edit/Delete → Add…→ select Type 1 Beam 3 → OK → input Real Constants Set No. : 1 , AREA: 2.19E-3,Izz: 3.83e-6(1 号实常数用于顶梁和侧梁) → Apply → input Real Constants Set No. : 2 , AREA: 1.185E-3,Izz: 1.87E-6 (2 号实常数用于弦杆) → Apply → input Real Constants Set No. : 3, AREA: 3.031E-3,Izz: 8.47E-6 (3 号实常数用于底梁) → OK (back to Real Constants window) → Close (the Real Constants window) (5) 定义材料参数 ANSYS Main Menu: Preprocessor → Material Props → Material Models → Structural → Linear → Elastic → Isotropic → input EX: 2.1e11, PRXY: 0.3(定义泊松比及弹性模量) → OK → Density (定义材料密 度) → input DENS: 7800, → OK → Close(关闭材料定义窗口) (6) 构造桁架桥模型 生成桥体几何模型 ANSYS Main Menu :Preprocessor → Modeling → Create → Keypoints → In Active CS → NPT Keypoint number:1,X,Y,Z Location in active CS:0,0 → Apply → 同样输入其余 15 个特征点坐标(最 左端为起始点, 坐标分别为 (4,0), (8,0), (12,0), (16,0), (20,0), (24,0), (28,0), (32,0), (4,5.5), (8,5.5), (12,5.5), (16.5.5), (20,5.5), (24,5.5), (28,5.5))→ Lines → Lines → Straight Line → 依次分别连接特征点 → OK 网格划分 ANSYS Main Menu: Preprocessor → Meshing → Mesh Attributes → Picked Lines → 选择桥顶梁及侧 梁 → OK → select REAL: 1, TYPE: 1 → Apply → 选择桥体弦杆 → OK → select REAL: 2, TYPE: 1 → Apply → 选择桥底梁 → OK → select REAL: 3, TYPE:1 → OK → ANSYS Main Menu:Preprocessor → Meshing → MeshTool → 位于 Size Controls 下的 Lines: → Element Size on Picked → Pick all → Apply Set → NDIV:1 → OK → Mesh → Lines → Pick all → OK (划分网格) (7) 模型加约束 ANSYS Main Menu: Solution → Define Loads → Apply → Structural→ Displacement → On Nodes → 选取桥身左端节点 → OK → select Lab2: All DOF(施加全部约束) → Apply → 选取桥身右端节点 → OK → select Lab2: UY(施加 Y 方向约束) → OK (8) 施加载荷 ANSYS Main Menu: Solution → Define Loads → Apply → Structural → Force/Moment → On Keypoints → 选取底梁上卡车两侧关键点(x 坐标为 12 及 20) → OK → select Lab: FY,Value: -5000 →

79

有限元分析基础教程 Apply → 选取底梁上卡车中部关键点(x 坐标为 16)→ OK → select Lab: FY,Value: -10000 → OK → ANSYS Utility Menu:→ Select → Everything (9) 计算分析 ANSYS Main Menu:Solution → Solve → Current LS → OK (10) 结果显示

曾攀

General Postproc → Plot Results → Deformed shape → Def shape only → OK (返 ANSYS Main Menu: 回到 Plot Results) → Contour Plot → Nodal Solu → DOF Solution, Y-Component of Displacement → OK (显示 Y 方向位移 UY)(见图 3-24(a)) 定义线性单元 I 节点的轴力 ANSYS Main Menu → General Postproc → Element Table → Define Table → Add → Lab: [bar_I], By sequence num: [SMISC,1] → OK → Close 定义线性单元 J 节点的轴力 ANSYS Main Menu → General Postproc → Element Table → Define Table → Add → Lab: [bar_J], By sequence num: [SMISC,1] → OK → Close 画出线性单元的受力图(见图 3-24(b)) ANSYS Main Menu → General Postproc → Plot Results → Contour Plot → Line Elem Res → LabI: [ bar_I], LabJ: [ bar_J], Fact: [1] → OK (11) 退出系统 ANSYS Utility Menu:File → Exit → Save Everything → OK

(a)桥梁中部最大挠度值为 0.003374m

(b)桥梁中部轴力最大值为 25380N

图 3-24 桁架桥挠度 UY 以及单元轴力计算结果

【ANSYS 算例】3.4.2(2) 基于命令流方式的桁架桥梁结构分析
!%%%%% [ANSYS 算例]3_4_2(2) %%%%% begin %%%%%% !------注:命令流中的符号$,可将多行命令流写成一行-----/prep7 /PLOPTS,DATE,0 !=====设置单元和材料 ET,1,BEAM3 R,1,2.19E-3,3.83e-6, , , , , R,2,1.185E-3,1.87e-6,0,0,0,0, !定义单元类型 !定义 1 号实常数用于顶梁侧梁 !定义 2 号实常数用于弦杆 !进入前处理 !设置不显示日期和时间

80

有限元分析基础教程 R,3,3.031E-3,8.47E-6,0,0,0,0, MP,EX,1,2.1E11 MP,PRXY,1,0.30 MP,DENS,1,,7800 !-----定义几何关键点 !定义 3 号实常数用于底梁 !定义材料弹性模量 !定义材料泊松比 !定义材料密度

曾攀

K,1,0,0,, $ K,2,4,0,, $ K,3,8,0,, $K,4,12,0,, $K,5,16,0,, $K,6,20,0,, $K,7,24,0,, $K,8,28,0,, $K,9,32,0,, $K,10,4,5.5,, $K,11,8,5.5,, $K,12,12,5.5,, $K,13,16,5.5,, $K,14,20,5.5,, $K,15,24,5.5,, $K,16,28,5.5,, !-----通过几何点生成桥底梁的线 L,1,2 $L,2,3 $L,3,4 $L,4,5 $L,5,6 $L,6,7 $L,7,8 $L,8,9 !------生成桥顶梁和侧梁的线 L,9,16 $L,15,16 $L,14,15 $L,13,14 $L,12,13 $L,11,12 $L,10,11 $L,1,10 !------生成桥身弦杆的线 L,2,10 $L,3,10 $L,3,11 $L,4,11 $L,4,12 $L,4,13 $L,5,13 $L,6,13 $L,6,14 $L,6,15 $L,7,15 $L,7,16 $L,8,16 !------选择桥顶梁和侧梁指定单元属性 LSEL,S,,,9,16,1, LATT,1,1,1,,,, !-----选择桥身弦杆指定单元属性 LSEL,S,,,17,29,1, LATT,1,2,1,,,, !-----选择桥底梁指定单元属性 LSEL,S,,,1,8,1, LATT,1,3,1,,,, !------划分网格 AllSEL,all LESIZE,all,,,1,,,,,1 LMESH,all /solu NSEL,S,LOC,X,0 D,all,,,,,,ALL,,,,, AllSEL,all NSEL,S,LOC,X,32 D,all,,,,,,,UY,,,, ALLSEL,all !------基于几何关键点施加载荷 FK,4,FY,-5000 /replot Allsel,all solve !=====进入一般的后处理模块 /post1 !后处理 $FK,6,FY,-5000 $FK,5,FY,-10000 !重画图形 !选择所有信息(包括所有节点,单元,载荷等) !求解 !根据几何位置选择节点 !对所选择的节点施加位移约束 !再恢复选择所有对象 !根据几何位置选择节点 !对所选择的节点施加位移约束 !再恢复选择所有对象 !再恢复选择所有对象 !对所有对象进行单元划分前的分段设置 !对所有几何线进行单元划分

!=====在求解模块中,施加位移约束、外力,进行求解

81

有限元分析基础教程 PLNSOL, U,Y, 0,1.0 PLNSOL, U,X, 0,1.0 !------显示线单元轴力-----ETABLE,bar_I,SMISC, 1 ETABLE,bar_J,SMISC, 1 PLLS,BAR_I,BAR_J,0.5,1 finish !画出轴力图 !结束 !显示 Y 方向位移 !显示 X 方向位移

曾攀

!%%%%% [ANSYS 算例]3_4_2(2) %%%%% end %%%%%%

【ANSYS 算例】3.4.2(3) 基于参数化方式的桁架桥梁结构分析
!%%%%%%% [ANSYS 算例]3_4_2(3) %%% begin /prep7 /PLOPTS,DATE,0 !=====进行参数化建模的参数设置 !(1)将桥梁总长度设为 L,每段水平桁架长度设为 DL,桥梁高度设为 H L=32 DL=L/8 H=5.5 !(2)将桥梁钢截面的面积设为 A,惯性矩设为 I A1=2.19E-3 I1=3.83E-6 A2=1.185E-3 I2=1.87E-6 A3=3.031E-3 I3=8.47E-6 !(3)将弹性模量和泊松比设为参数 e_modu=2.1e11 prxy_Poi=0.3 !(4)将载荷值设为参数 P1=-5000 ET,1,BEAM3 R,1,A1,I1, , , , , R,2,A2,I2,0,0,0,0, R,3,A3,I3,0,0,0,0, MP,EX,1,e_modu MP,PRXY,1,prxy_Poi MP,DENS,1,,7800 !-----定义关键点 !定义材料密度 $P2=-10000 $P3=-5000 !======下面开始有限元的建模和分析 !定义单元类型 !定义 1 号实常数用于顶梁侧梁 !定义 2 号实常数用于弦杆 !定义 3 号实常数用于底梁 %%%%%%% !------注:以下命令流中的符号$,表示可将多行命令流写成一行----------!进入前处理 !设置不显示日期和时间

!-------定义钢的弹性模量和泊松比

82

有限元分析基础教程 K,1,0,0,, K,7,6*DL,0,, $K,2,DL,0,, $K,3,2*DL,0,, $K,4,3*DL,0,, $K,10,DL,H,, $K,16,7*DL,H,, $K,5,4*DL,0,, $K,11,2*DL,H,, $K,6,5*DL,0,, $K,12,3*DL,H,,

曾攀

$K,8,7*DL,0,, $K,9,8*DL,0,,

K,13,4*DL,H,, $K,14,5*DL,H,, $K,15,6*DL,H,, !-----通过几何点生成桥底梁的线 L,1,2 L,9,16 L,2,10

$L,2,3 $L,3,4 $L,4,5 $L,5,6 $L,6,7 $L,7,8 $L,8,9 $L,15,16 $L,14,15 $L,13,14 $L,12,13 $L,11,12 $L,10,11 $L,1,10 $L,3,10 $L,3,11 $L,4,11 $L,4,12 $L,4,13 $L,5,13 $L,6,13 $L,6,14 $L,6,15 $L,7,15 $L,7,16 $L,8,16

!------生成桥顶梁和侧梁 的线 !------生成桥身弦杆的线 !------选择桥顶梁和侧梁指定单元属性 LSEL,S,,,9,16,1, LATT,1,1,1,,,, !-----选择桥身弦杆指定单元属性 LSEL,S,,,17,29,1, LATT,1,2,1,,,, !-----选择桥底梁指定单元属性 LSEL,S,,,1,8,1, LATT,1,3,1,,,, !------划分网格 AllSEL,all LESIZE,all,,,1,,,,,1 LMESH,all /solu NSEL,S,LOC,X,0 D,all,,,,,,ALL,,,,, AllSEL,all NSEL,S,LOC,X,L D,all,,,,,,,UY,,,, ALLSEL,all !-----施加载荷 FK,4,FY,P1 /replot Allsel,all solve !=====进入一般的后处理模块 /post1 PLNSOL, U,Y, 0,1.0 PLNSOL, U,X, 0,1.0 !------显示线单元轴力 ETABLE,bar_I,SMISC, 1 !后处理 !显示 Y 方向位移 !显示 X 方向位移 !选择所有信息(包括所有节点,单元,载荷等) !求解 $FK,6,FY,P3 $FK,5,FY,P2 !选择 x=L 处的节点 !施加 UY 方向固定的约束 !选择 x=0 处的节点 !施加完全固定的约束 !再恢复选择所有对象 !对所有对象进行单元划分前的分段设置 !对所有几何线进行单元划分

!=====在求解模块中,施加位移约束、外力,进行求解

83

有限元分析基础教程 ETABLE,bar_J,SMISC, 1 PLLS,BAR_I,BAR_J,0.5,1 finish !结束 end %%%%% !%%%%%%% [ANSYS 算例]3_4_2(3) %%

曾攀

3.5 本章要点
·有限元分析的基本方程(刚度矩阵、载荷列阵、单元的刚度方程、组装后的刚度方程) ·1D 问题(杆及梁)的虚功原理求解 ·1D 问题(杆及梁)的最小势能原理求解 ·杆单元的位移模式、刚度矩阵及坐标变换 ·梁单元的位移模式、刚度矩阵及坐标变换 ·平面梁单元的等效节点载荷的计算方法

3.6 习题
3-1 试用虚功原理推导杆 1D 单元的单元方程。

? x? ? x? 3-2 推导二次函数变截面的杆单元的刚度矩阵 A( x) = ?1 ? ? A1 + ? ? A2 。 ? l? ?l?
3-3 对于沿轴向承受均布载荷的杆单元,推导相应的节点等效载荷。 3-4 对于沿轴向承受三角形载荷分布的杆单元,推导相应的节点等效载荷。 3-5 有一个变截面的杆结构如图所示, 材料的弹性模量为 E = 2 × 10 MPa , 试采用两个单元(其
5

2

中一个为变截面单元)来分析该问题,并得到整个结构所有的力学信息。

(a) 一个变截面的杆结构 题图 3-5

(b)两个单元的计算模型

3-6 有一个杆结构,受有沿轴向的分布载荷,杆的横截面积为 A = 2 × 10 cm ,材料的弹性模
2

?2

量为 E = 2 × 10 MPa ,分析该问题的建模方案,并得到整个结构所有的力学信息。
5

84

有限元分析基础教程

曾攀

题图 3-6

3-7 有一个变截面的杆结构, 受有沿轴向的均匀分布载荷, 试采用如图 3-7(b)所示的多种建模方 式来进行分析,得到整个结构所有的力学信息,并比较几种建模方式的结果及特点。

(a) 一个变截面的杆结构 题图 3-7

(b)几种计算模型

3-8 有一个桁架结构,如题图 3-8 所示,每根杆件的材料以及面积都相同,试采用有限元方法 分析该结构的受力状态,并求出支反力。

题图 3-8

3-9 试用虚功原理推导 2 节点平面梁单元的单元方程。 3-10 有一个变截面的梁单元,假设梁的长度为 l,梁的横截面积为 A( x) = ?1 ?

? ?

x? ? x? ? A1 + ? ? A2 , l? ?l?

85

有限元分析基础教程

曾攀

惯性矩为 I ( x) = ?1 ?

? ?

x? ? x? ? I1 + ? ? I 2 ,试推导该纯弯梁单元的刚度矩阵。 l? ?l?

3-11 对于一个梁单元,承受有三角形分布的垂直载荷,推导该单元的相应节点等效载荷。 3-12 有一个变截面的梁结构,还受有承受有梯形分布的垂直载荷,如题图 3-12 所示,试采用如 题图 3-12(b)所示的多种建模方式来进行分析,得到整个结构所有的力学信息,并比较几种建模 方式的结果及特点。

(a) 一个变截面的梁结构 题图 3-12

(b)几种计算模型

3-13 有一个梁结构如题图 3-13 所示,试采用梁单元来分析该结构,并求取支反力。

题图 3-13

3-14 有一个悬臂简支梁如图所示, 在其中的一段受有均布载荷, 试分析该结构, 并求取支反力。

86

有限元分析基础教程

曾攀

题图 3-14

3-15 一 座 两 层 建 筑 可 以 简 化 为 如 题 图 3-14 所 示 的 模 型 , 梁 的 参 数 为 A = 4 × 10 m ,
2

?2

I = 6.75 ×10?4 m 4 ,材料的弹性模量为 E = 2 ×105 MPa ,试分析该结构,并求取支反力。

题图 3-15

L ,如题图所 3-16 有一个细长的弹性曲梁,当它未受力时,它由一个常曲率半径 R,并且 R 示,该曲梁的左端点固支,并与刚性水平面相切;现在它的右端点施加一个垂直集中力 F,若 试使用多个梁单元进行分析, 该梁仍保持线弹性, 多大的 F 才能使曲梁的右端点的曲率变为零。 并比较单元的个数与计算精度之间的关系。

题图 3-16

87

有限元分析基础教程

曾攀

3-17 在【ANSYS 算例】3.4.2(3)的基础上,当卡车从左端开到右端一系列过程情况下,即将载 荷 P1_P2_P3 从左端到右端进行一系列平移(每步移动五段中的一段),分别计算桥梁中最大位移 及应力的变化。

题图 3-17

88

有限元分析基础教程

曾攀

第 4 章 连续体结构分析的有限元方法
杆梁结构由于有自然的连接关系, 可以凭一种直觉将其进行自然的离散, 而连续体则不 同,它的内部没有自然的连接节点,必须完全通过人工的方法进行离散;有人说有限元方法 的真正魅力在于它成功地处理了连续体(场)问题,人们公认的有限元方法的开拓者之一 Courant,就是在 1943 年处理连续体问题时使用了三角形区域的分片连续函数和最小势能原 理;而“有限单元”的名称,是 Clough 于 1960 年在处理平面连续体问题时正式提出的。 本章将就连续体问题的有限元方法进行全面的讨论和研究。

4.1 连续体结构分析的工程概念
最简单的 1D 几何图形就是直线,最简单的 2D 几何图形就是三角形,最简单的 3D 几 何图形就是四面体,杆单元以及梁单元从局部坐标系来看可以说是 1D 单元,本章将重点讨 论用于连续体结构离散的 2D 单元及 3D 单元。 第 3 章虽然仅讨论了杆单元以及梁单元, 但从中已经给出了有限元分析的具有标准化分 析过程的核心要素:单元 本章除了使读者接受将连续的几何域划分为一系列的 2D 单元或 3D 单元这一基本思想 外,其它的过程与第 3 章完全相同,将再一次说明了有限元分析的标准化及规范性的特征。 图 4-1 展示了将连续的一个圆离散为有限个三角片的过程, 正是按照这一原理, 同样也 可以将一个任意复杂的几何形状离散为一系列(有限个)标准几何图形(单元)的组合, 如图 4-2 所示。

图 4-1 将连续的一个圆离散为有限个三角片

图 4-2 将连续的几何域离散为由有限个三角片组成的逼近域

4.2 连续体结构分析的基本力学原理
与第 3 章类似,连续体结构的力学分析也包括:基本变量、基本方程、求解原理,下 面分别进行描述,有关的详细内容见参考文献[17]及[18]。

89

有限元分析基础教程

曾攀

【基本变量】4.2.1(1) 连续体问题的三大类变量 就基本变量而言,2D 和 3D 问题将在 1D 问题上进行另外方向上的推广延伸,注意:除 了在主方向上进行延伸外,还存在有每两个坐标轴之间的交叉(夹角)项,对于应力来说就是 剪应力,对于应变来说就是剪应变;几种情况的基本变量列于表 4-1 中。
表 4-1 连续体问题的三大类变量的各个分量(直角坐标系下) 1D 问题
位移分量 应变分量

2D 问题

3D 问题

u (x)

u ( x, y ), v( x, y )

u ( x, y, z ), v( x, y, z ), w( x, y, z )

ε xx ( x)
σ xx ( x)

ε xx ( x, y ), ε yy ( x, y ), γ xy ( x, y )

ε xx ( x, y, z ), ε yy ( x, y, z ), ε zz ( x, y, z ), γ xy ( x, y, z ), γ yz ( x, y, z ), γ xz ( x, y, z )

应力分量

σ xx ( x, y ), σ yy ( x, y ), τ xy ( x, y )

σ xx ( x, y, z ), σ yy ( x, y, z ), σ zz ( x, y, z ), τ xy ( x, y, z ), τ yz ( x, y, z ), τ xz ( x, y, z )

注:σ xx ( x) 中的第一个下标表明应力的方向,第二个下标表明应力所作用面的法线方向;对于应变也如此。 一般情况下,将下标相同的正应力 σ xx ( x, y, z ) 或正应变 ε yy ( x, y, z ) 表达成 σ x ( x, y, z ) 或正应变 ε y ( x, y, z ) , 即只写一个下标,y 方向以及 z 方向的情况类似。

【基本方程】4.2.1(2) 连续体问题的三大类方程及边界条件 对于基本方程而言,2D 和 3D 问题同样将在 1D 问题上进行另外方向上的推广延伸; 几种情况的基本方程列于表 4-2 中。
表 4-2 连续体问题的三大类方程的各个表达式(直角坐标系下)

1D 问题
平衡方程

2D 问题
?σ xx ?τ xy + =0 ?x ?y ?τ xy ?σ yy + =0 ?x ?y

3D 问题
?σ xx ?τ xy ?τ xz + + =0 ?x ?y ?z ?τ xy ?σ yy ?τ yz + + =0 ?x ?y ?z ?τ xz ?τ yz ?σ zz + + =0 ?x ?y ?z
ε xx = γ xy
?u ?v ?w , ε yy = , ε zz = , ?x ?y ?z ?v ?u ?w ?v ?w ?u = + , γ yz = + , γ zx = + ?x ?y ?y ?z ?x ?z

dσ xx =0 dx

几何方程

ε xx =

du dx

?u ?v , ε yy = ?x ?y ?v ?u γ xy = + ?x ?y

ε xx =

物理方程

ε xx =

σ xx
E

ε xx = ε yy =

?σ xx ? μσ yy ? ? E? ?σ yy ? μσ xx ? ? E?
1 G 1

1

ε xx = ε yy = ε zz =

?σ xx ? μ (σ yy + σ zz ) ? ? E?
1 ?σ yy ? μ (σ xx + σ zz ) ? ? E?

1

γ xy = τ xy

?σ zz ? μ (σ xx + σ yy ) ? ? E?
1 1 1 G G G

1

γ xy = τ xy , γ yz = τ yz , γ zx = τ zx

注: E , μ , G 为材料的弹性模量,泊松比,剪切模量。

90

有限元分析基础教程

曾攀

对于一般的力学问题, 还有两类边界条件(Boundary Condition), 位移边界条件 BC(u) 即, 以及力边界条件 BC(p);几种情况的边界条件列于表 4-3 中。
表 4-3 连续体问题的两大类边界条件的各个表达式(直角坐标系下) 1D 问题 几何边 界条件 BC(u) 外力边 界条件 BC(p) 2D 问题 3D 问题

u ( x) x = x = u
0

u ( x, y ) x = x , y = y = u
0 0

u ( x, y , z ) x = x , y = y , z = z = u
0 0 0

v ( x, y ) x = x , y = y = v
0 0

v ( x, y , z ) x = x , y = y , z = z = v
0 0 0

w( x, y, z ) x = x , y = y , z = z = w
0 0 0

σ xx ( x) x = x = px
0

nxσ xx ( x0 , y0 ) + nyτ xy ( x0 , y0 ) = px

nxσ xx ( x0 , y0 , z0 ) + nyτ xy ( x0 , y0 , z0 ) + nzτ xz ( x0 , y0 , z0 ) = px nxτ xz ( x0 , y0 , z0 ) + nyτ yz ( x0 , y0 , z0 ) + nzσ zz ( x0 , y0 , z0 ) = pz

nxτ xy ( x0 , y0 ) + nyσ yy ( x0 , y0 ) = p y nxτ xy ( x0 , y0 , z0 ) + nyσ yy ( x0 , y0 , z0 ) + nzτ xz ( x0 , y0 , z0 ) = p y

注: x0 , y0 , z0 为边界上的几何坐标; nx , n y , nz 边界外法线上的方向余弦; u , v , w 为给定的对应方 向上的位移值;

px , p y , pz 为给定的对应方向上的边界分布力。

【求解原理】4.2.1(3) 直接法以及试函数法的求解思想 针对一个具体对象,目标是通过三大类方程(平衡、几何、物理)来求解出三大类变量(位 移、应变、应力);对于 1D 问题实际就是通过三个方程来求解三个变量,一般可以直接进 行联立求解;而对与 2D 问题,变量及方程的个数都较多,而且方程为偏微分方程,除一些 具有简单几何形状(如矩形、圆形)外,一般较难直接求解。必须寻找能够求解任意复杂几何 形状的通用方法;一个好的方法应该具有以下几个特征: 几何形状的适应性(通用方法) 数学力学原理上的标准化(不需要较多的技巧,都较容易掌握) 处理过程的规范性(具有统一的处理流程) 操作实施的可行性(具备相应的软件及硬件条件,可处理大规模线性方程组) 分析误差的可控制性(可确定计算的最佳效率,方法具有良好的收敛性) 针对弹性问题的三大类变量和方程,从第一层次(求解方法论)的思路上说,总体上有 两大类求解方法:直接法(直接求解原始微分方程)、以及试函数法(基于假设解的调试方法)。 从第二层次(求解策略)的思路上说,对于三大类变量,要同时进行联立求解一般比较 困难,需要将变量先进行代换,最好是代换成一种类型的变量先进行求解,然后再求解另外 类型的变量。 对于直接法,也就是解析方法,若将要求解的基本变量先确定为应力的话,这就是应 力方法,最常用的有 Airy 应力函数(stress function)方法,还有一些诸如逆解法、半逆解法的 方法,但这些方法虽然可以获得解析解,但仅针对一些简单几何形状,而大量的实际问题目 前还不能获得解析解;而且,对求解的数学技巧也要求较高。因此,这类方法大部分不具备 上面提到的几个特征。 而对于试函数法,可以追溯到 Rayleigh(1870 年)、以及 Ritz(1909 年),它首先假设一 个可能的解(即试函数,其中包含一些待定系数),将试函数再带入到原方程中,通过确定相 应误差函数的最小值来获得其中的待定系数,这样就求出以试函数(确认了待定系数)为结果 的解,这种方法大大降低了求解的难度和技巧,方法也具有标准化和规范性,但由于最早的 试函数是基于整个几何全域来选择的,因此,它的几何形状的适应性也是受到一定限制的, 加上过去还没有合适的计算工具,它实施的可行性也不具备。 进入 20 世纪 50 年代, 随着计算机技术的飞速发展, 复杂力学问题的处理有了两个实质 性的突破,其一:大规模计算成为可能,其二:将基于整个几何全域的试函数变为基于分片
91

有限元分析基础教程

曾攀

(子区域)的多项式函数表达(可以很好适应任意复杂的几何形状),然后再将分片的函数进行 集成组合得到全域的试函数。这种分片就是“单元” ,分片的过程就是将整体区域进行离散 的过程,将这种基于分片函数描述的试函数方法叫做有限元方法[19],相关的示意见图 4-3。

图 4-3 基于整体区域离散的有限元方法

【求解原理】4.2.1(4) 连续体问题求解的虚功原理 对于一般弹性问题(以 2D 问题为例),在几何域 Ω 中,受有体积力( b x , b y ),在外力边界

S p 上,受有施加的分布力( px , p y )。设有满足位移边界条件的位移场 (u, v) ,这就是试函数
(其中有一些待定的系数),则它的虚位移为 (δ u , δ v) ,虚应变为 (δε xx , δε yy , δγ xy ) 。
相应的虚应变能为

δ U = ∫ (σ xxδε xx+σ yyδε yy+τ xyδγ xy )dΩ
Ω

(4-1)

而外力虚功为

δ W = ∫ (b x ? δ u + b y ? δ v)dΩ + ∫ ( px ? δ u + p y ? δ v)dA
Ω Sp

(4-2)

那么,虚功原理 δU = δW 可以表达为



Ω

(σ xxδε xx+σ yyδε yy+τ xyδγ xy )dΩ = ∫ (b xδ u + b yδ v)dΩ + ∫ ( pxδ u + p yδ v)dA
Ω Sp

(4-3)
通过以上方程就可以确定出试函数中的待定系数。 【求解原理】4.2.1(5) 连续体问题求解的最小势能原理 同样对于 2D 问题,设有满足位移边界条件的位移场 (u , v) ,即试函数(其中有一些待定 的系数),确定其中待定系数的方法,就是使得该系统的势能取极小值,即
( u , v )∈BC ( u )

min

Π (u , v)

(4-4)

其中

92

有限元分析基础教程

曾攀

Π = U ?W 1 = ∫ (σ xx ? ε xx+σ yy ? ε yy+τ xy ? γ xy )dΩ 2 Ω ? ∫ (b x ? u + b y ? v)dΩ + ∫ ( px ? u + p y ? v)dA
Ω Sp

(4-5)

实际上,虚功原理与最小势能原理是等价的。 【强度准则】4.2.1(6) 结构分析中的受力状态诊断(强度准则) 实际的工程设计往往需要进行多次的力学分析、优化修改才能完成,强度准则在整个 设计过程中具有重要的作用,它是判断受力状态是否满足需要的主要依据,图 4-4 展示了一 个结构件的设计过程[21]。

图 4-4 一个结构的力学分析、状态诊断以及优化改进的主要过程

对于不同的材料,由于它的承载破坏的机理不同,需要采用针对性的强度判断准则。 下面给出常用的几种强度准则。一个问题具体应采用何种准则,应根据材料的受力状态、环 境要求、设计要求,甚至还需要通过一系列的实际状况的实验来确定。 (1)最大拉应力准则(max. tensile stress criterion):若材料发生脆性断裂失效,其原因是材料内 所承受的最大拉应力达到了所能承受的极限 (一般用于脆性材料) 。 已知危险点的应力状态 σ ij ,首先通过斜面分解方法求出最大的拉应力 σ 1 (实际上就是 第一主应力)和所在面的主方向,然后进行应力失效校核和判断

σ 1 ≤ [σ ]
上式中 [σ ] 为材料的许用应力,由材料的单向拉伸试验和安全系数确定,即 [σ ] =

(4-6)

σs
n

,其中

σ s 为单向拉伸试验得到的屈服应力, n 为安全系数。

93

有限元分析基础教程

曾攀

(2)最大剪应力准则(max. shearing stress criterion):若材料发生屈服(或剪断),其原因是材料 内所承受的最大剪应力达到了所能承受的极限(一般用于韧性材料)。 已知危险点的应力状态 σ ij ,首先通过斜面分解方法求出最大的剪应力 τ max

τ max =
然后进行剪应力失效校核和判断

σ1 ? σ 3
2

(4-7)

τ max ≤ [τ ]

(4-8)

其中 [τ ] 为材料的许用剪应力,由材料的单向拉伸试验和安全系数确定。对于材料的单向拉 伸试验,有 σ 3 = 0 ,因此

[τ ] =

[σ ] 2

(4-9)

将(4-7)和(4-9)式代入(4-8)中,有以主应力形式来表达的最大剪应力准则

σ 1 ? σ 3 ≤ [σ ]

(4-10)

(3)最大畸变能准则(max. distortion energy criterion):若材料发生屈服(或剪断),其原因是材 料内的畸变能密度达到了所能承受的极限(一般用于韧性材料,也称为 Mises 等效应力强度 准则)。 已知危险点的应力状态 σ ij ,按下式计算该点的畸变应变能密度

′ Ud =

1+ μ ?(σ 1 ? σ 2 ) 2 + (σ 2 ? σ 3 ) 2 + (σ 1 ? σ 3 ) 2 ? ? 6E ? ′ ′ U d ≤ [U d ]

(4-11)

然后进行畸变应变能密度校核和判断

(4-12)

′ 其中 [U d ] 为材料的临界畸变应变能密度,即由材料的单向拉伸试验和安全系数确定。对于
′ 材料的单向拉伸试验,有 σ 2 = 0, σ 3 = 0 ,因此,由式(4-11),有 [U d ] = 1+ μ [σ ]2 , 3E

将该关系以及式(4-11)代入(4-12)中,有以主应力形式来表达的最大畸变能准则

1 ?(σ 1 ? σ 2 ) 2 + (σ 2 ? σ 3 ) 2 + (σ 1 ? σ 3 ) 2 ? ≤ [σ ] ? 2?
或写成更一般的形式

(4-13)

1 2 2 2 ? (σ xx ? σ yy ) 2 + (σ yy ? σ zz ) 2 + (σ xx ? σ zz ) 2 + 6(τ xy + τ yz + τ xz ) ? ≤ [σ ] (4-14) ? ? 2
若定义

94

有限元分析基础教程

曾攀

σ eq =

1 ? (σ 1 ? σ 2 ) 2 + (σ 2 ? σ 3 ) 2 + (σ 1 ? σ 3 ) 2 ? ? 2?

(4-15)

1 2 2 2 ? (σ xx ? σ yy ) 2 + (σ yy ? σ zz ) 2 + (σ xx ? σ zz ) 2 + 6(τ xy + τ yz + τ xz ) ? = ? 2?

则称 σ eq 为 Mises 等效应力,也叫做应力强度,由上式可以看出,该等效应力反映了材料受 力变形的畸变能的平方根。

4.3 平面问题有限元分析的标准化表征
4.3.1 平面问题的 3 节点三角形单元描述
平面问题 3 节点单元具有几何特征简单、描述能力强的特点,是平面问题有限元分析 中最基础的单元,也是最重要的单元之一。 【单元构造】4.3.1(1) 平面问题的 3 节点三角形单元

(1) 单元的几何和节点描述 3 节点三角形单元(3-node triangular element) 如图 4-5 所示。3 个节点的编号为 1、2、3, 各自的位置坐标为(xi, yi), i=1,2,3,各个节点的位移(分别沿 x 方向和 y 方向)为(ui, vi), i=1,2,3。

图 4-5

平面 3 三节点三角形单元

如图 4-5 所示,该单元共有 6 个节点位移自由度(DOF)。将所有节点上的位移组成一个 列阵,记作 q ;同样,将所有节点上的各个力也组成一个列阵,记作 P ,那么
e e

(6×1)

q e = [u1 v1 u2
Py1 Px 2

v2
Py 2

u3

v3 ]
Px 3

T

(4-16)
T

P e = ? Px1 ? (6×1)

Py 3 ? ?

(4-17)

若该单元承受分布外载,可以将其等效到节点上,即也可以表示为如图 4-5 所示的节点力。 利用函数插值、几何方程、物理方程以及势能计算公式,可以将单元的所有力学参量(即场 变量: u ( x, y ) , ε ( x, y ) ,σ ( x, y ) 和 Π )用节点位移列阵 q 及相关的插值函数来表示;下
e
e

面进行具体的推导。 (2) 单元位移场的表达
95

有限元分析基础教程

曾攀

就如图 4-5 所示的平面 3 节点三角形单元,由于有 3 个节点,每一个节点有两个位移, 因此共有 6 个节点位移,考虑到简单性、完备性、连续性及待定系数的唯一确定性原则,分 别选取单元中各个方向的位移模式为

u ( x, y ) = a0 + a1 x + a2 y ? ? ? v ( x, y ) = b0 + b1 x + b2 y ? ?
由节点条件,在(x=xi,y=yi)处,有

(4-18)

u ( xi , yi ) = ui v ( xi , yi ) = vi

? ? ? i=1,2,3 ? ?

(4-19)

将式(4-19)代入节点条件式(4-20)中,可求解出式(4-18)中的待定系数,即

u1 1 a0 = u2 2A u3 1 u1 1 a1 = 1 u2 2A 1 u3

x1 x2 x3 y1

y1 y2 = y3 y2 = y3 u1 u2 = u3

1 ( a1u1 + a2u2 + a3u3 ) 2A

(4-20)

1 ( b1u1 + b2u2 + b3u3 ) 2A

(4-21)

1 x1 1 a2 = 1 x2 2A 1 x3

1 ( c1u1 + c2u2 + c3u3 ) 2A

(4-22)

b0 =

1 ( a1v1 + a2v2 + a3v3 ) 2A 1 b1 = ( b1v1 + b2v2 + b3v3 ) 2A 1 b2 = ( c1v1 + c2v2 + c3v3 ) 2A
1 x1 1 A = 1 x2 2 1 x3 y1 y2 = y3

(4-23) (4-24) (4-25)

在式(4-20)~式(4-25)中

1 1 ( a1 + a2 + a3 ) = ( b1c2 ? b2c1 ) 2 2

(4-26)

a1 =

? y2 = x2 y3 ? x3 y2 ? y3 ? ? 1 y2 ? b1 = ? = y2 ? y3 ? 1 y3 ? ? 1 x2 c1 = = ? x2 + x3 ? 1 x3 ? ? x2 x3

(1,2,3)

(4-27)

上式中的符号(1,2,3)表示下标轮换,如 1→2,2→3,3→1 同时更换。 将系数式(4-20)~式(4-25)代入式(4-18)中,重写位移函数,并以节点位移的形式进行表
96

有限元分析基础教程

曾攀

示,有

u ( x, y ) = N1 ( x, y ) ? u1 + N 2 ( x, y ) ? u2 + N 3 ( x, y ) ? u3 v ( x, y ) = N1 ( x, y ) ? v1 + N 2 ( x, y ) ? v2 + N 3 ( x, y ) ? v3
写成矩阵形式,有

(4-28) (4-29)

?u ( x, y ) ? ? N1 u ( x, y ) = ? ?= ( 2×1) v ( x, y ) ? ? 0 ? ?

0 N1

N2 0

0 N2

N3 0

? u1 ? ?v ? ? 1? 0 ? ? u2 ? e N ? ? v ? = (2×6) ( x, y ) ? q N3 ? ? 2 ? (6×1) ? u3 ? ? ? ? v3 ? ? ?

(4-30)

其中 N( x, y ) 为形状函数矩阵,即

?N N ( x, y ) = ? 1 (2×6) ?0


0 N1

N2 0

0 N2

N3 0

0? N3 ? ?

(4-31)

Ni =
其中的系数 ai , bi , ci 见式(4-27)。 (3) 单元应变场的表达

1 ( ai + bi x + ci y ) 2A



i=1,2,3

(4-32)

由弹性力学平面问题的几何方程(矩阵形式)

? ?u ? ? ? ? ? ? ?ε xx ? ? ?x ? ? ?x ? ? ? ?v ? ? ε ( x, y ) = ?ε yy ? = ? ?=? 0 ( 3×1) ?y ? ? ?γ xy ? ? ? ? ? ?u ?v ? ? ? ? + ? ? ? ?y ?x ? ? ?y ? ? ?
其中 [ ? ] 为几何方程的算子矩阵(operator matrix),即

? 0? ? ? ? ? u ( x, y ) ? ? = [? ] u ?? ?y ? ? v ( x, y ) ? ( 3×2) ( 2×1) ?? ? ?x ? ?

(4-33)

?? ? ? ?x ? [? ] = ? 0 ? ?? ? ? ?y ?
将式(4-30)代入式(4-33)中,有

? 0? ? ?? ? ?y ? ?? ? ?x ? ?

(4-34)

97

有限元分析基础教程

曾攀

( 3×1)

ε ( x, y ) = [ ? ] N ( x , y ) ? q e = B ( x, y ) ? q e
( 3×2 ) ( 2×6 ) ( 6×1)
(3×6)

(4-35)

(6×1)

其中几何矩阵 B( x, y ) 为

?? ? ? ?x ? B ( x, y ) = [ ? ] N = ? 0 (3×6) ( 3×2 ) ( 2×6 ) ? ?? ? ? ?y ?
将式(4-32)代入上式,有

? 0? ? ? ? ? N1 ?y ? ? 0 ?? ?? ? ?x ? ?

0 N1

N2 0

0 N2

N3 0

0? N3 ? ?

(4-36)

?b1 0 b2 1 ? 0 c1 0 B ( x, y ) = (3×6) 2A ? ?c1 b1 c2 ?
其中

0 c2 b2

b3 0 c3

0? ? c3 ? = ? B1 ? (3×2) b3 ? ? ?

( 3×2)

B2

? B3 ? ( 3×2 ) ?

(4-37)

?bi 1 ? Bi = 0 2A ? (3×2) ?ci ?
(4) 单元应力场的表达

0? ci ? , ? bi

相关文章:
有限元分析基础教程
有限元分析基础教程_数学_自然科学_专业资料。有限元分析基础教程 前言 有限元分析已经在教学、科研以及工程应用中成为重要而又普及的数值分析方法和工具;该基 础教程...
1 有限元分析基础教学大纲
七、教材及主要参考书 教材:石伟,有限元分析基础与应用教程,机械工业出版社,2010. 主要参考书: [1] 王新荣, 初旭宏. ANSYS 有限元基础教程. 电子工业出版社....
有限元分析基础教学大纲
有限元分析基础 Fundamentals of Finite Element Analysis 一、课程基本信息学分:...七、教材与参考书 1、 《有限单元法原理及应用简明教程》 ,高秀华、张小江、...
UG有限元分析入门与提高实例教程G
UG有限元分析入门与提高实例教程G_计算机软件及应用_IT/计算机_专业资料。UG有限元分析实例培训大纲 UG 有限元分析入门与提高实例教程 第一讲: UG 高级仿真模块...
有限元分析的基础教程chapter3
有限元分析基础教程chapter3_机械/仪表_工程科技_专业资料。有限元分析基础教程,从有限元的基础理论3 弹性力学平面问题的有限元法本章包括以下的内容: 3.1 弹性...
有限元分析的基础教程chapter5
有限元分析基础教程chapter5_机械/仪表_工程科技_专业资料。有限元分析基础教程,从有限元的基础理论5.等参单元本章包括以下内容: 5.1 等参单元的基本概念 5...
edu_ecologychuanke1477652465
视频教程,艺人教育全套教学,在线学习软件技术课程,UG NX 8.0动力学与有限元分析入门到精通视频下载
有限元分析开题报告
水利电力出版社,1993 5 有限元分析及应用 [2] 曾攀。有限元分析基础教程。清华大学,2008 [3] 石伟。有限元分析基础与应用教程。机械工业出版社,2010 [4] ...
西工大有限元法研究生课程008
有限元分析基础教程(曾攀) 343页 免费 有限元法讲义-研究生 175页 5财富值 弹性力学及有限元讲义-研究... 110页 免费 西工大有限元法研究生课程... 4页 10...
授课内容-基础知识
3页 免费 有限元分析基础教程(曾攀) 343页 免费如要投诉违规内容,请到百度文库投诉中心;如要提出功能问题或意见建议,请点击此处进行反馈。 ...
更多相关标签:
有限元基础教程 | 有限元分析 | 有限元分析教程 | 有限元分析公司 | 有限元分析培训班 | 有限元分析培训 | 有限元分析基础 | 有限元分析视频教程 |