吐血分享：研究生非线性有限元(双语)作业.doc

非线性有限元分析非线性有限元分析 NONLINEAR FINITE ELEMENT ANALYSIS 作业题目 平面应变问题四节点四边形单元求解 宋盈盈 S12080203065 任忠凯 S12080203061 黎 原 S12080203030 郝守义 S12080203020 杨海涛 S12080203099 苏永林 S12080203067 2013 年年 6 月月 非线性有限元分析非线性有限元分析 平面应变问题四节点四边形单元求解 硕士研究生：黎 原 宋盈盈 任忠凯 郝守义 杨海涛 苏永林 指导教师：孙静娜 学科专业：机械设计及理论 所 在 单 位：机械工程学院 Nonlinear finite element analysis FOUR-NODE QUADRILATERAL ELEMENT SOLUTION FOR PLANE STRAIN PROBLEM By Group Li Supervisor: Professor Sun Jingna Yanshan University June, 2013 摘 要 I 摘摘 要要 对于大多数工程问题的数值求解，有限元方法已经成为一种有效的求解工 具，其所涉及和应用的领域及其广泛，不仅包括汽车、飞机、建筑、桥梁等的 结构变形和应力分析，还包含热流、热体流、磁通量、渗流等流动问题的场分 析。随着计算机技术和 CAD 技术的发展，许多复杂问题的求解通过计算机的 建模分析已经很容易就能实现了。有限元的基本过程是：将被研究物体的复杂 几何区域离散为简单的限单元，而单元内的材料性质和控制方程通过单元节点 的未知量来进行表达，再通过单元集成、载荷和约束条件的处理，得到方程组， 求解该方程组就可以得到该物体行为的近似表达式。本文根据平面应变问题的 特点，采用四节点四边形单元，对所研究问题的截面进行离散和单元划分，通 过设定形函数，列出问题的载荷矩阵以及相应的刚度矩阵，借助于 MATLAB 工具及其相应程序，编辑输入文件，从而得出问题的近似解。此外，本文还针 对所研究的问题，运用 ANSYS 有限元分析软件，对其进行有限元分析，通过 建立几何模型，划分网格，施加载荷和约束设定求解模式，得出比较精确的求 解办法。基于编程计算结果和软件分析结果，本文对以上两种求解过程以及对 应的求解结果进行了对比和分析，从而进一步加强了对有限元求解问题过程的 理解和掌握。 关键词：关键词：有限元分析；四节点；四边形单元；平面应变； 非线性有限元分析 II Abstract For most of the numerical solution of engineering problems, finite element method has become an effective tool for solving, which involves a wide application area, including not only structural deformation and stress analysis for automobiles, airplanes, buildings, bridges and other, but also field analysis of heat, hot body float, flux, seepage flow. With the development of computer technology and CAD technology, many complex issues solution has been easily achieved through computer modeling analysis. The basic process of finite element is as follow: discrete the complex geometry area of the object being studied into the simple finite element zone, material properties and control equations in the element are expressed via the unknown quantity for the node. The equations are obtained through integration units, loads, and constraints treatment, solving the equations，we can get the approximate expression of the body act. Based on the characteristics of plane strain problem, we use the four-node quadrilateral element to discrete and mesh the sections of problem under study, by setting the shape function and listing loading matrix and the corresponding stiffness matrix of problem, by means of MATLAB tools and corresponding program, we edit the input file so as to arrive approximate solution. Additionally, for the research question, in this article, we also use ANSYS, finite element analysis software, analysis the problem, by establishing the geometric model, meshing, applying loads and constraints and setting the solving mode to obtain more precise way to the solution. Based on the programming results and analysis results by software, the two solution process and the corresponding results above were compared solving and analyzed in the paper, thus, we can understand and master the further finite element problem solving process. Keywords: finite element analysis; four nodes; quadrilateral element; plane strain; 目 录 III 目目 录录 摘 要.1 Abstract .II 第一章 MATLAB 的编程求解.1 1.1 MATLAB 简介.1 1.2 四边形单元求解原理.1 1.2.1 形状函数.1 1.2.2 刚度矩阵.2 1.2.3 载荷矩阵.3 1.3 求解题目及建模.4 1.4 利用 QUAD 程序求解.5 第二章 ANSYS 的分析求解 .9 2.1 ANSYS 软件简介 .9 2.2 ANSYS 软件的优势 .9 2.3 问题求解过程.9 2.3.1 建立几何模型.9 2.3.2 网格划分.10 2.3.3 施加载荷和约束.11 2.3.4 选定求解类型及求解模式.12 2.3.5 查看结果.12 总 结.14 致 谢.15 参考文献.16 MATLAB 的编程求解 1 第一章第一章 MATLAB 的编程求解的编程求解 1.1 MATLAB 简介 MATLAB 是由美国 mathworks 公司发布的主要面对科学计算、可视化以及交 互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化 以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗 环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提 供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言 （如 C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平，在 数学类科技应用软件中在数值计算方面首屈一指。 MATLAB 可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、 连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通 讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB 的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形 式十分相似，故用 MATLAB 来解算问题要比用 C，FORTRAN 等语言完成相同的事 情简捷得多，并且 MATLAB 也吸收了像 Maple 等软件的优点，使 MATLAB 成为一 个强大的数学软件。 1.2 四边形单元求解原理 1.2.1 形状函数 首先讨论规则的基准单元的形状函数，在 坐标系中定义形状为正方 形的基准单元。定义 Lagrange 形状函数为，即在节点 i 处，为 1，在其 i N i N 他节点处为零（i=1，2，3，4） 。如形状函数，在节点 1 处=1，而在节 i N 1 N 1 N 点 2，3，4 处要求=0。求得所有的四个形状函数为 1 N （1 11 4 1 11 4 1 11 4 1 11 4 1 41 3 2 1 N N N N 1） 非线性有限元分析 2 图 1.1 坐标系下的四边形单元（基准单元） 基于单元节点值，可以描述单元的位移场。为点处的位移分量，则, =NqNq。其中 q q 表示单元节点的位移列阵，N N 。 432 432 1 1 00 000 0 0 0NNN NNN N N J J 为 Jacobian 矩阵， J J yx yx J=J= 43214321 43214321 111)1 (111)1 ( 111)1 (111)1 ( 4 1 yyyyxxxx yyyyxxxx （1 2221 1211 JJ JJ 2） 1.2.2 刚度矩阵 求得的维数为 88 单元刚度矩阵为 （1 1 1 1 1 detdJdDBBtk Te 3） 其中 t 为单元厚度； MATLAB 的编程求解 3 B=AG， (14) 12221121 1121 1222 00 00 det 1 JJJJ JJ JJ J A G= （1 10101010 10101010 01010101 01010101 4 1 5） 1.2.3 载荷矩阵 本题中不考虑体力，也不受集中载荷，所以只考虑面力 T 的作用，面力为 作用在直线上的线性均布压力载荷。假如和是施加在直线 12 上的法向压 1 p 2 p 力，则， ，其中和 11 cpTx 22 cpTx 11 spTy 22 spTy 21 21 l xx s 。 21 21 l yy c 面力载荷矩阵为 e T （1 T yyxxyyxx e TTTTTTTT tl T 21212121 21 2222 6 6） 本题中边 412 受线性载荷，将其分为 48 和 812 两段直线来分析。其 中 4-8 边 29 5 29 2 84 48 84 84 l yy c l xx s 非线性有限元分析 4 MPaT MPaT MPaT MPaT y y x x 29 150 75 29 2 29 100 50 29 2 29 375 75 29 5 29 250 50 29 5 8 4 8 4 所以边 4-8 的载荷矩阵为 e T 84 边 8-12 T e tttt T 6 13000 6 100000 6 3500 6 87500 84 29 5 29 2 128 812 128 128 l yy c l xx s MPaT MPaT MPaT MPaT y y x x 29 200 100 29 2 29 150 75 29 2 29 500 100 29 5 29 375 75 29 5 12 8 12 8 边 8-12 的载荷矩阵为 e T 128 T e tttt T 6 55000 6 137500 6 50000 6 125000 128 将厚度 t=2000mm 代入，并将两矩阵进行叠加得最终的载荷矩阵 T 为 （1 783 . 1 758 . 4 00000071 . 275 . 700 0000717 . 1 792 . 2 000000 eeee eeT 7） 1.3 求解题目及建模 本问题要求采用四节点四边形单元（QUAD）求解如下平面应变问题。其中 E=70GPa，=0.33，厚度 t=2000mm。 MATLAB 的编程求解 5 图 1.2 求解问题的描述 针对该问题，将二维求解域离散为 6 个直边四节点四边形单元，如图 1.3 所示。 图 1.3 有限元离散模型 单元的节点连接信息如表 1.1 所示。 表 1.1 单元的节点连接信息 四个节点编号单元编号 1234 11265 非线性有限元分析 6 22376 33487 456109 5671110 6781211 1.4 利用 QUAD 程序求解 利用已有的程序 QUAD.m，需编制输入文件。利用 1.3 节中建立的有限元模 型和式（17）可以编制输入文件如下所示。 PROBLEM 7.4 12 6 1 2 4 2 8 6 0 1 0 0 2 250 0 3 500 0 4 1000 0 5 0 500 6 250 500 7 500 500 8 800 500 9 0 1000 10 250 1000 11 500 1000 12 600 1000 1 1 2 6 5 1 2000 0 2 2 3 7 6 1 2000 0 3 3 4 8 7 1 2000 0 4 5 6 10 9 1 2000 0 MATLAB 的编程求解 7 5 6 7 11 10 1 2000 0 6 7 8 12 11 1 2000 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 7 -2.92e7 8 -1.17e7 15 -7.5e7 16 -2.1e7 23 -4.58e7 24 -1.83e7 1 70000 .33 12E-6 B1 i B2 j B3 (Multi-point constr. B1*Qi+B2*Qj=B3) 其中 NN节点数； NE单元数； NMmat 数； MDIM节点自由度数； NEN单元节点数； NDN节点自由度； ND边界位移约束的数目； NL边界载荷的数目； NMPC是否具有多点约束，有取 1，没有取 0。 求解的结果如下所示 节点位移为 Node# X-Displ Y-Displ 非线性有限元分析 8 1 -6.4468E-006 -1.1483E-005 2 -2.6974E-006 -7.0127E-006 3 -6.7647E-006 1.2464E-006 4 -9.5912E-006 8.5794E-006 5 -1.6944E+000 -9.5758E-001 6 -1.5757E+000 -3.2290E-001 7 -1.6108E+000 -3.9874E-003 8 -1.8640E+000 5.1072E-001 9 -3.2439E+000 -1.1600E+000 10 -3.2656E+000 -4.5804E-001 11 -3.3815E+000 1.6052E-001 12 -3.4702E+000 4.0050E-001 位移约束节点的支反力 DOF# Reaction 1 3.7922E+007 2 6.7547E+007 3 1.5867E+007 4 4.1251E+007 5 3.9792E+007 6 -7.3320E+006 7 5.6419E+007 8 -5.0467E+007 单元的等效应力 ELEM# von Mises Stresses at 4 Integ_points 1 1.5592E+002 1.3173E+002 8.1479E+001 1.1541E+002 2 1.3459E+002 1.3405E+002 1.0054E+002 1.0146E+002 3 1.4225E+002 1.6360E+002 1.4472E+002 1.0813E+002 4 4.2109E+001 4.4549E+001 3.0822E+001 2.7810E+001 5 8.6471E+001 9.2103E+001 6.4856E+001 5.7359E+001 6 9.8635E+001 1.0097E+002 9.3342E+001 8.8395E+001 MATLAB 的编程求解 9 非线性有限元分析 10 第二章第二章 ANSYS 的分析求解的分析求解 本章运用 ANSYS 有限元分析软件对问题进行了模拟和分析，并得出了分析 的图像结果，意在与上一章的编程计算结果进行对比，从中进行总结，从而进 一步加深对有限元求解方法、求解过程以及软件求解的理解和掌握。 2.1 ANSYS 软件简介 ANSYS 有限元软件是一个多用途的有限元法计算机设计程序，可以用来求 解结构、流体、电力、电磁场及碰撞等问题，应用于航空航天、汽车工业、生 物医学、桥梁、建筑、电子产品、重型机械、微机电系统、运动器械等众多领 域。 ANSYS 软件主要包括三个部分：前处理模块，分析计算模块和后处理模块。 前处理模块提供了一个强大的实体建模及网格划分工具，用户可以方便地构造 有限元模型。分析计算模块包括结构分析、流体动力学分析、电磁场分析、声 场分析、压电分析以及多物理场的耦合分析，可模拟多种物理介质的相互作用， 具有灵敏度分析及优化分析能力。后处理模块可将计算结果以彩色等值线显示、 梯度显示、矢量显示、粒子流迹显示、立体切片显示、透明及半透明显示等图 形方式显示出来，也可将计算结果以图表、曲线形式显示或输出。 2.2 ANSYS 软件的优势 ANSYS 的技术涵盖多个学科领域。不论是需要结构分析、流体、热力、电 磁学、显式分析、系统仿真还是数据管理，ANSYS 的产品均能为各个行业的企 业提供巨大帮助。ANSYS 的成套产品极具灵活性。不论是为企业中新手还是能 手使用；是单套部署还是企业级部署；是首次通过还是复杂分析；是桌面计算、 并行计算还是多核计算，这一工程设计的高扩展性均能满足当前与未来的需求。 ANSYS 是唯一一家能提供客户所需能力水平的仿真软件供应商，而且能随此类 需求的发展无限扩展。工程设计与开发可使用多种 CAD 产品、内部开发代码、 物料库、第三方求解器、产品数据管理流程等其他工具。与那些刻板、僵化的 系统不同，ANSYS 的软件具有开放性和适应性特性，能实现高效的工作流程。 此外，其产品数据管理可使知识和经验在工作组间与企业内的实现共享。 2.3 问题求解过程 2.3.1 建立几何模型 如图所示，截面为梯形，上底为 600mm，下底为 1000mm，高为 1000mm ANSYS 的分析求解 11 2.3.2 网格划分 根据边长进行网格划分，如下图所示： 非线性有限元分析 12 2.3.3 施加载荷和约束 梯形的底边固定，梯形的腰部施加非均匀分布载荷：从 50Mpa 到 100Mpa， 加载情况如图所示： 载荷以约束的施加情况如下图： ANSYS 的分析求解 13 2.3.4 选定求解类型及求解模式 材料的弹性模量 E=70Gpa，泊松比为：=0.33，采用静力分析如下图所示： 2.3.5 查看结果 总体的变形图为： Mises 应力云图为： 非线性有限元分析 14 Mises 应变云图为： ANSYS 的分析求解 15 总总 结结 本文运用有限元四节点四边形单元的理论，对平面应变问题进行分析。首 先，把平面应变问题划分单元，建立模型，并运用 MATLAB 的相关程序，编辑输 入文件，通过运行程序最终得到结果。然后，再通过 ANSYS 软件进行实际建模， 设定网格划分，运用其自带的程序进行运算求解。经过上述两种求解过程得到 结果的差异还是很明显的。 用 MATLAB 编程所得到的数值结果是解析解，如果不考虑结果数值的精度， 该解应该是在现有模型网格划分基础上的精确解，如果扩大划分网格的数目， 则所得的结果将是更精确的近似解。与此相比，用 ANSYS 软件进行分析所得到 的结果是更加接近与实际的解一方面，软件进行网格化的数目比较大，网格越 多就越接近最终的实际解；另一方面，软件本身具有参数以及运算模式设定的 优势，ANSYS 会根据问题的性质以及用户的一些自定义设定，同时也考虑到实 际情况的影响因素进行优化处理，因此，得到的结果更加精确一些。 总的来说，MATLAB 编程得到的解精确，但与实际问题的真实解差别比较大。 增大网格划分数量可以改进真实度，但相应的程序计算效率将会降低。用 ANSYS 分析得到的结果与实际结果更加接近一些，但由于软件自身参数设定的 原因，有可能与实际结果产生叫的偏差。增加网格数量，细化参数的设定将开 始结果更理想，但这样对计算机的内存以及处理能力是个挑战。 非线性有限元分析 16 致致 谢谢 短短六周的非线性有限元分析课程结束了。在这些宝贵的六周的课程学习 期间，我们收获了很多知识，过得很充实也很快乐。 感谢孙静娜老师在有限元授课期间对我们的关心和帮助。我们都非常喜欢 孙老师的讲课风格和授课方式。孙老师总是那么和蔼可亲，平易近人，不要求 出勤率，只要求自己尽最大的努力给来听课的学生讲好每一堂课。夏天的天气 非常干燥闷热，但孙老师每堂课都讲的那么努力、耐心、细心，丝毫没有抱怨。 每次上课之前，孙老师爬楼梯进教室，身上都是汗。由于这堂课是临时决定由 孙老师来讲授，所以临危受命，老师需要加紧备课，做课件，熟悉内容，可以 说，老师尽职尽责，兢兢业业。老师的表现深深的感动了每一个听课的学生， 也让每一个对这门课程感兴趣的同学不仅收获了知识，还收获了老师的那份爱。 因此，借此课程结束的时机，您的学生黎原，代表我们小组，也代表所有喜欢 孙老师的同学，对孙老师这一学期在非线性有限元分析授课阶段所做出的 努力和付出致以最真挚的感谢和最崇高的敬意，老师您辛苦啦！真诚祝福老师： 在以后的工作中，顺利、开心、快乐；在以后的生活中，美满、如意、幸福。 致 谢 17 参考文献参考文献 1 Tirupathi R.Chandrupatla Ashok D. Belegundu, 曾攀. 工程中的有限元方法M. 第 3 版. 清华大学出版社 2 于润伟, 朱晓慧. MATLAB 基础及应用M. 第 2 版. 机械工业出版社. 3 郑阿奇, 曹弋编.