ABAQUSFortran二次开发.docx

《ABAQUSFortran二次开发.docx》由会员分享，可在线阅读，更多相关《ABAQUSFortran二次开发.docx（81页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、ABAQUSFortran二次开发目 录摘 要I1.绪论11.1.课题的研究背景11.2.本文的研究内容和方法22.基于软件的二次开发32.1.介绍32.2.各模块简介32.3.的二次开发平台52.4.的二次开发语言63.用户材料子程序83.1.开发环境设置83.2.注意事项93.3.接口的原理103.4.的使用方法124.材料非线性问题144.1.材料的弹塑性本构关系144.2.非线性有限元算法理论174.3.增量理论常刚度法公式推导204.4.增量理论切线刚度法公式推导215.程序设计和编码255.1.本构关系描述255.2.常刚度法程序设计275.3.常刚度法程序编码295.4.切线刚度

2、法程序设计325.5.切线刚度法程序编码355.6.程序的调试386.程序验证406.1.问题描述406.2.本构关系416.3.自带材料模型计算416.4.常刚度法的验证436.5.切线刚度法的验证456.6.两种算法的比较分析477.结论及展望517.1.结论517.2.展望51致 谢53参考文献54附1：自带弹塑性材料验证的文件55附2：用于算法验证的文件6177 / 81摘 要 软件功能强大，特别是能够模拟复杂的非线性问题，它包括了多种材料本构关系及失效准则模型，并具有良好的开放性，提供了若干个用户子程序接口，允许用户以代码的形式来扩展主程序的功能。 本文主要研究了用户子程序的开发方法

3、，采用语言编制了各向同性硬化材料模型的接口程序，研究该类材料的弹塑性本构关系极其实现方法。 本文紧紧围绕的二次开发技术，首先对其接口原理做了详细介绍，然后针对非线性有限元增量理论中的常刚度法和切线刚度法的算法理论做了深入的剖析，推导出了常刚度法和切线刚度法的算法理论的具体表达式，然后分别编制了两种算法的程序，最后建立了一个具体的验算模型，通过及自带弹塑性本构关系的计算结果相比较，验证两者的正确性。 本文还对常刚度法和切线刚度法得算法效率做了对比，得出了在非线性程度较高时切线刚度法效率高于常刚度法的结论。关键字： 、有限元、材料非线性、切线刚度 , , a , a , a . , , . , ,

4、 a , , , , a , . a a . ：、 、 、 1. 绪论1.1. 课题的研究背景 有限单元法基本思想的提出，可以追溯到克劳夫()在1943年的工作1，他第一次尝试应用定义在三角形区域上的分片连续函数和最小位能原理相结合，来求解. 扭转问题。1960年克劳夫进一步处理了平面弹性问题，并第一次提出了“有限单元法”的名称，使人们开始认识了有限单元法的功效。 四十多年来，随着电子计算机的广泛应用和发展，有限单元法的理论和应用都得到迅速的，持续不断的发展，其应用己由弹性力学平面问题扩展到空间问题、板壳问题，由静力学问题扩展到稳定问题、动力问题和波动问题。分析的对象从弹性材料扩展到塑性、粘弹

5、性、粘塑性和复合材料等，从固体力学扩展到流体力学、传热学等连续介质力学领域。在工程分析中的作用已从分析和校核扩展到优化设计并和计算机辅助设计。 利用有限元软件解决工程和科学问题，是有限元理论应用于工程设计和科学研究实践的主要形式。由于工程设计的巨大市场需要，有限元软件的发展是很迅速的，目前常用的大型有限元软件常见的有2000， ，等，这些软件的共同特点是具有丰富的单元库和求解器，强大而可靠的分析功能，人们利用这些软件解决了很多工程建设和工业产品设计中遇到的问题，取得了巨大的经济技术效益。 由于工程问题的千差万别，不同的用户有不同的专业背景和发展方向，通用软件不免在具体的专业方面有所欠缺，针对这

6、些不足，大部分的通用软件都提供了二次开发功能，以帮助用户减少重复性的编程工作、提高开发起点、缩短研发周期、降低开发成本，并能简化后期维护工作，给用户带来很多方便。基于通用软件平台进行开发，是目前研究的一个重要发展方向。 也提供了若干用户子程序（ ）接口，它是一个功能非常强大且适用的分析工具，及命令行的程序格式相比，用户子程序的限制少得多，从而使用更加灵活方便。针对所提供的本构关系模型种类有限，无法满足工程应用需要的问题，用户子程序中的用户材料子程序（ ，简称）接口可以帮助用户定义自己的材料本构模型和算法，这是的独到之处。由于其操作方便，能被灵活地应用于各个领域中，尤其受到用户的青睐。1.2.

7、本文的研究内容和方法 中用户材料子程序的开发主要解决两方面的问题：本构模型的建立和积分算法的选择。 本文主要研究非线性材料的实现方法，并重点研究其迭代算法部分，目前，用户材料子程序的迭代算法主要是常刚度法，常刚度法的优点在于算法原理较简单，程序编写较方便，缺点是当遇到复杂非线性材料时，其迭代次数较多，收敛速度也较慢，在这个情况下，本文采取的是一种迭代次数较少且收敛速度较快的切线刚度法，具体就是采用语言编制了基于模型的接口程序,并采用切线刚度算法，通过及自带本构关系计算的结果相比较，验证其正确性。 本文的研究工作紧紧围绕的二次开发技术，首先根据有限元方法推导材料非线性问题算法的公式，然后参考接口

8、规范设计程序的算法流程，继而编写出该程序，最后建立一个具体的本构和具体的模型做测试，验证程序的正确性，在这一过程中，调试是一个非常重要的过程，占用了大量的时间，在调试程序时采用了将中间变量输出到文本的方式，这样能明确跟进迭代过程，发现算法或程序的缺陷。 本文采用的本构关系是经过归纳和抽象的，也就是说本文的程序并不仅仅是只针对某个具体模型和问题，而是针对所有符合抽象出的各向同性硬化材料，这样做的好处是能保证程序的通用性和复用性，避免以后的重复劳动，当然，这也是符合软件设计接口的宗旨的。 2. 基于软件的二次开发2.1. 介绍 是一套功能强大的基于有限元法的工程模拟软件2，其解决问题的 范围从相对

9、简单的线性分析到最富有挑战性的非线性模拟问题。具备十分丰富的、可模拟任意实际形状的单元库。并及之对应拥有各种类型的材料模型库，可以模拟大多数典型工程材料的性能，其中包括金属、橡胶、高分子材料、复合材料、钢筋混凝土、可压缩弹性的泡沫材料以及岩石和土这样的地质材料。 作为通用的模拟分析工具， 不仅能解决结构分析中的问题，还能模拟和研究各种领域中的问题，如热传导、质量扩散、电子元器件的热控制（热一电耦合分析）、声学分析、土壤力学分析（渗流应力耦合分析）和压电介质力学分析。 为用户提供了广泛的功能，且使用起来又十分简明。最复杂的问题也可以很容易地建立模型3。例如复杂的多部件问题可以通过对每个部件定义材

10、料模型和几何形状，然后再把它们组装起来而构成。在大部分模拟分析问题中，甚至在高度非线性问题中，用户也只需要提供结构的几何形状、材料性能、边界条件和荷载工况这样的工程数据就可以进行分析。在非线性分析中，能自动选择合适的荷载增量和收敛精度。不仅能选择这些参数值，而且能在分析过程中不断地调整参数来保证有效地得到高精度的解，很少需用户去定义这些参数。 2.2. 各模块简介 有两个主要的分析模块： 和 。还有两个特殊用途的附加分析模块：和。另外，还有 分别及，和 的接口模块：，和 。是完全的工作环境模块，它包括了模型的构造，交互式提交作业、监控作业过程以及评价结果的能力。是的子集，它具有后处理功能，这些

11、模块之间的关系见图2- 1图2-1 是一个通用分析模块，在数值方法上采用有限元方法常用的隐式积分。它能够求解广泛的线性和非线性问题，包括结构的静态、动态问题、热力学场和电磁场问题等。对于通常同时发生作用的几何、材料和接触非线性可以采用自动控制技术处理，也可以由用户自己控制。 是一个在数值方法上采用有限元显式积分的特殊模块，它利用对时间的显式积分求解动态有限元方程。它适合于分析诸如冲击和爆炸这样短暂、瞬时的动态问题，同时对高度非线性问题如模拟加工成型过程中接触条件的改变等也非常有效。 是进行有限元分析的前后处理模块，也是建模、分析和后处理的人机交互平台。该模块根据结构的几何图形生成网格，将材料和

12、截面的特性分配到网格上，并施加载荷和边界条件。该模块可以进一步将生成的模型投入到分析模块中进行高效率的后台运行，并对运行情况进行监测，对计算结果进行后处理。的后处理支持分析模块的所有功能，并且对计算结果的描述和解释提供了范围很广的选择，除了通常的云图，等值线和动画显示之外，还可以用列表，曲线(包括部分常用运算)等其他常用工具来完成对结果数据的处理。该模块的许多独特功能及特点，例如特征化建模、参数化建模、适应设计者要求的数据管理系统等极大的方便了的使用者。 的一系列功能可以附加在中应用。它偏向于模拟海上结构，如海洋石油平台。它的功能包括模拟波浪，风载荷及浮力的 影响。在本指南中不讨论。 允许有限

13、元模型作为柔性部件进入到产品族中去进行分析。 把注模分析软件中有限元网格、材料性质和初始应力数据转换成为 输入文件。 的一系列功能可附加在 中进行设计敏度计算。 模块把 分析软件中的有限元模型信息 转换成 输入文件的一部分。2.3. 的二次开发平台 的脚本语言接口非常友好，其自嵌的脚本语言是4，系国际上广泛使用、功能强大、具有良好开放性的一种面向对象程序设计语言。所以，应用在中进行二次开发也比较方便，且可移植性强。以基于的语法规则向二次开发者提供了许多库函数，这些库函数主要是用来增强的交互式()操作功能。用户可以通过的交互式()界面实现分析对象的特征造型、指定材料属性、完成网格剖分和控制、提交

14、并监控分析作业，也可以使用脚本语言越过的交互式()界面直接高效地向内核提交任务。使用可以进行参数化建模，修改交互式建立的模型，还可以一次提交多个作业。 出了脚本语言接口，还为用户提供了功能强大的用户子程序接口（ ），以帮助用户开发基于内核的程序，常用的用户子程序包括( ，用户单元子程序)，( a s ，用户材料子程序 )5，其中的使用最为广泛，它主要用于用户开发自己的材料模型，以弥补自带材料模型的不足，帮助用户完成各种材料分析，功能极为强大。 在国外，众多的有限元分析和研究者热衷于使用，一个很重要的原因就在于给用户提供了功能强大，使用方便的二次开发工具和接口，使得用户可以方便的进行富含个性化的

15、有限元建模、分析和后处理，满足特定工程问题的需要。通过用户材料子程序接口，用户可定义任何补充的材料模型，不但任意数量的材料常数都可以作为资料被读取，而且对于任何数量的及解相关的状态变量在每一材料计数点都提供了存储功能，以便在这些子程序中应用。2.4. 的二次开发语言 的二次开发语言主要有3种：，语言主要用于开发，语言主要用于用户子程序开发，而语言主要专注于其他高级开发部分。 本文主要是针对用户子程序的开发，所以采用语言，下面简要介绍一下该语言极其特点： 语言是世界上第一个被正式推广使用的高级语言6。它是1954年被提出来的，1956年开始正式使用，至今已有三十多年的历史，但仍历久不衰，它始终是

16、数值计算领域所使用的主要语言。 语言是 的缩写，意为“公式翻译”。它是为科学、工程问题或企事业管理中的那些能够用数学公式表达的问题而设计的，其数值计算的功能较强。 语言问世以来，根据需要几经发展，先后推出了不同的版本，主要版本有 77， 90， 95，采用 77，通常用固定格式编写代码。77语言同C语言一样，是一种结构化编程语言结构化程序设计方法规定，在结构化的程序中，只能有三种基本结构：(1)顺序结构这是一种最简单的基本结构形式，它的特点是，在这个结构内的各个功能模块或语句序列，是按其出现的先后顺序执行的，如赋值语句、输入/输出语句等。它有一个入口和一个出口，并在入口和出口之间包含着若干个功

17、能块，其中每一个功能块可以是一个非转移语句。因此，顺序基本结构块是由一系列的顺序执行语句组成的。(2)分支选择结构在给定的条件下，分支选择结构判断选择哪一条路径执行，不同路径完成的功能是不同的。实现分支选择结构主要由块语句、语句、 语句以及 语句组成的结构。(3)循环结构循环结构也称重复处理结构，即重复执行某一功能块，直到满足（或不满足）某一条件为止。实现循环结构的90语句主要是语句、块语句和逻辑语句的结合。以上三种基本结构，是组成结构化程序的基本结构形式。这里有两层意思：一是结构化的程序中，各个模块均由这三种基本结构组成；二是结构化程序本身，从宏观上也是这三种基本结构形式之一。3. 用户材料

18、子程序3.1. 开发环境设置 由于是采用语言编写，那么要运行就需要安装的开发环境， 同时还需要的支持，本文采用的版本为6.81，支持 9.1-10.1， 安装时又需要安装 的相应版本，经过比较，本文选用6.81 10.1 2005，相对于来说，开发环境的设置较为繁琐，这给子程序的使用带来诸多不便，为了解决这一问题，我用语言编制了子程序编译环境设置工具，只需要将安装文件解压到的安装目录，运行安装程序就可以了，整个过程不需要人工干预，也不需要安装庞大的 2005，如图3-1所示 图3-13.2. 注意事项 的用户子程序是根据提供的相应接口，按照语法，用户自己编写的代码。它是一个独立的程序单元，可以

19、独立的被存储和编译，也能被其它程序单元引用，因此，利用它可带回大量数据供引用程序使用，也可以用它来完成各种特殊的功能。它的一般结构形式是： S(x12,)（用于用户子程序中） )（用于用户子程序中）x1，x2，是提供的用户子程序的接口参数，有些参数是传到用户子程序中的，例如 中的、，有些是需要用户自己定义的，例如F，文件和随着软件的安装而包含在操作系统中，它们含有重要的参数，帮助主求解程序对用户子程序进行编译和链接。当控制遇到语句时便返回到引用程序单元中去，语句是用户子程序结束的标志。 在一个算例中，用户可以用到多个用户子程序，但必须把它们放在一个以为扩展名的文件中。运行带有用户子程序的算例同

20、时有两种方法：一是在中运行，在 菜单中的子菜单的 对话框中选择用户子程序所在的文件即可；另外是在中运行，语法如下： 编制用户子程序时应注意(1)用户子程序相互之间不能调用，但可以调用用户自己编写的子程序和应用程序。应用程序必须由用户子程序调用。当用户编写子程序时，建议子程序名以K开头，以免和内部程序冲突。(2)当用户在用户子程序中利用打开外部文件时，要注意以下两点：一是设备号的选择是有限制的，只能取1518和大于100的设备号，其余的都已被占用；二是用户需提供外部文件的绝对路径而不是相对路径。(3)对于不同的用户子程序调用的时间是不同的，有的是在每个的开始，有的是的结尾，有的是在每个的开始等等

21、。当调用用户子程序时，都会把当前的和利用用户子程序的两个实参和传给用户子程序，用户可把它们输出到外部文件中，这样就可清楚的知道何时调用该用户子程序。为保证用户子程序的正确执行，子程序的书写必须遵循的相关规定，下面以用户材料子程序为例详细说明。3.3. 接口的原理 用户材料子程序（ ，简称）是提供给用户定义自己的材料属性的程序接口78，使用户能使用材料库中没有定义的材料模型。用户材料子程序通过及主求解程序的接口实现及的资料交流。在输入文件中，使用关键词“* ”表示定义用户材料属性。子程序具有强大的功能，使用子程序：（1）可以定义材料的本构关系，使用材料库中没有包含的材料进行计算，扩充程序功能。（

22、2）几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予中的任何单元。（3）必须在中提供材料本构的雅可比（）矩阵，即应力增量对应变增量的变化率。由于主程序及之间存在数据传递，甚至共享一些变量，因此必须遵守有关的书写格式，中常用的变量在文件开头予以定义，通常格式 (,1 ,2 ,3 ,4 01)*80 ()(),1 ()()(),2 ()()(2)(1)(1),3 ()(3)(3,3)0(3,3)1(3,3) 中的应力矩阵、应变矩阵以及矩阵、等，都是直接分量存储在前，剪切分量存储在后。直接分量有个，剪切分量有个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写时要考虑到所使

23、用单元的类别。下面对中用到的一些变量进行说明：( )：一个的矩阵，称作矩阵是应力的增量，是应变的增量，(i，j)表示增量步结束时第j个应变分量的改变引起的第i个应力分量的变化。通常矩阵是一个对称矩阵，除非在“* ”语句中加入了“”参数。()：应力张量数组，对应个直接分量和个剪切分量。在增量步的开始，应力张量矩阵中的数值通过和主程序之间的接口传递到中，在增量步的结束将对应力张量矩阵更新。对于包含刚体转动的有限应变问题，一个增量步调用之前就已经对应力张量进行了刚体转动，因此中只需处理应力张量的共旋部分。中应力张量的度量为柯西应力。()：用于存储及解有关的状态变量的数组，在增量步开始时将数值传递到中

24、，也可在子程序或中先更新数据，然后增量步开始时将更新后的资料传递到中。在增量步的结束必须更新状态变量矩阵中的数据。和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的资料更新以外，及解有关的状态变量矩阵中的任何向量或者张量都必须通过旋转来考虑材料的刚体运动。状态变量矩阵的维数通过输入文件中的关键词“*”定义，关键词下面数据行的数值即为状态变量矩阵的维数。()：材料常数数组。材料常数的个数，等于关键词“* ”中“”常数设定的值。矩阵中元素的数值对应于关键词“ ”下面的数据行。，：分别定义每一增量步的弹性应变能，塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅作为能量输出。()：应变数

25、组。()：应变增量数组。：增量步的时间增量。：直接应力分量的个数。：剪切应力分量的个数。：总应力分量的个数，。由于子程序在单元的积分点上调用，增量步开始时，主程序路径将通过的接口进入，单元当前积分点必要变量的初始值将随之传递给的相应变量。在结束时，变量的更新值将通过接口返回主程序。3.4. 的使用方法 我们知道，有限元计算（增量方法）的基本问题7是：已知第n 步的结果（应力，应变等），然后给出一个应变增量，计算，要完成这一计算，并要计算()=。是应力增量矩阵，是应变增量矩阵，() 定义了第J 个应变分量的微小变化对第 I 个应力分量带来的变化。 该矩阵只影响收敛速度，不影响计算结果的准确性 （

26、当然，不收敛自然得不到结果）。 有限元计算的中心问题就是求得节点的位移 （进而应变、应力），以使内力和 外力达到平衡： （3-1）d 是节点位移矩阵，黑体字表示矩阵或矢量。除了小变形、线弹性问题，方程2-1是非性的，要用迭代的方法解出： （3-2） (3-3) i 表示一个增量步内的第i 次迭代，n表示第n个增量步。是切线刚度，由材料的 矩阵结合单元计算组装而得。刚度矩阵其实就是力对位移的梯度。要想快速收敛，位移增量应沿该梯度方向变化，也就是说，如果 矩阵不是那么准确，自然 也不怎么准确，那么满足3-1式的位移被找到的速度也就变慢，甚至发散，根本找不到。但收敛速度无论慢快，3-1式才是判断结果

27、准确及否的唯一标准。所以 矩阵不影响结果的准确性，只影响收敛速度的快慢。 以大变形、非性材料为例，整个计算步骤是这样的： 整个外力不是一次加上，而是一点点加上的，不然会发散得不到结果的。所以，每一个增量步开始时就是在原来的外力上加上一点点，得到。根据3-2得到位移增量，此时要知道力对位移的梯度，以尽快找到满足平衡条件的位移，由材料的 矩阵和单元结合起来组装得到(此处使用 提供的 矩阵)。然后可计算应变增量，调用，得到新的应力，进而得到新的内力，所以，程序不在乎新的应力是由增量方法得到，还是全量方法得到，而只在乎新应力是否准确。然后回到3-2，如此循环，直至3-2右端为0，也即满足3-1。这样第

28、1 步就完成了，然后开始第2 步，即 外力加上一点点，按同样的方法求解新的位移。直至整个外力全部施加并得到满足3-1的位移。4. 材料非线性问题 弹性力学作为精确理论，从本质上都是非线性的，早期, 和在这些方面都作出了重要贡献。后来又提出超弹性(即具有弹性势的)有限变形理论，由于理论方程的冗长而复杂，且工程应用也没有提出这方面要求而被搁置。20世纪40年代以后，由于橡胶材料、高分子合成材料的迅速发展和工业领域的大量应用，非线性弹性及超弹性的研究再次引起科学和工程界的重视，除了一般理论研究有了新的发展以外，工程应用计算方法也得到长足的发展。4.1. 材料的弹塑性本构关系弹塑性材料进入塑性的特征是

29、当荷载卸去后存在不可恢复的永久变形。所以，在卸载情况下，应力应变之间不再是唯一的对应关系。这是区别于非线性弹性材料的基本属性。只以加载时应力应变关系成非线性，还不足以判断材料是非线性弹性还是弹塑性。但是一经卸载就可以看出两者的区别。非线性弹性材料沿原路径返回，而弹塑性材料将依据不同的加载历史卸载后产生不同的永久变形。对大多数材料来说，在单调加载的情况下，存在一个明显的极限应力，当应力低于时，材料保持线弹性。而当应力达到以后，则材料开始进入弹塑性状态。如继续加载，然后在卸载，材料始终保持永久的塑性变形。如果应力达到后，应力不再增加，而材料变形可以继续增加，及变形处于不定的流动状态，则称材料为理想

30、弹塑性的。反之如果应力达到后，再增加变形，应力也必须增加，则材料是应变硬化的，这时应力是塑性应变的函数，可解析为： (4-1)本构关系反应着应力应变之间的关系。对于弹性材料变形是可以恢复的；而塑性材料变形是不可以恢复的。典型的弹塑性应变在卸载后要保持一个永久的变形。如图3-2图4-1 塑性应变有下列特性：（1）总应变分为弹性和塑性两部分，即 (4-2) 或者： (4-3) （2）塑性变形取决于加载路径，而应力应变之间没有一一对应的关系。所以必须确定二则之间的本构关系，这种本构关系可以用偏微分方程或者增量形式来描述。总之，弹塑性理论主要包括以下几个方面：（1）应变张量的分解；（2）应力空间的屈服

31、条件；（3）流动法则；（4）强化法则；（5）协调性条件。1：本构模型塑性力学的应力-应变曲线通常有5种简化模型8：（1）理想弹塑性模型，用于低碳钢或强化性质不明显的材料。（2）线性强化弹塑性模型，用于有显著强化性质的材料。（3）理想刚塑性模型，用于弹性应变比塑性应变小得多且强化性质不明显的材料。（4）线性强化刚塑性模型，用于弹性应变比塑性应变小得多且强化性质明显的材料。（5）幂强化模型，为简化计算中的解析式，可将应力-应变关系的解析式写为y（y）n，式中y为屈服应力，y为及y相对应的应变，n为材料常数。 图4-22：屈服条件 在复杂应力状态下，判断物体屈服状态的准则称为屈服条件9。屈服条件是各

32、应力分量组合应满足的条件。对于金属材料，最常用的屈服条件为最大剪应力屈服条件（又称屈服条件）和弹性形变比能屈服条件（又称 条件）。对于岩土材料则常用屈服条件、屈服条件和屈服条件。对于强化或软化材料，屈服条件将随塑性变形的增长而变化，改变后的屈服条件称为后继屈服条件。当已知主应力的大小次序时，使用屈服条件较为方便；若不知道主应力的大小次序，则使用 屈服条件较为方便。对于韧性较好的材料， 屈服条件及试验数据符合较好。 屈服准则具体形式是，对于各项同性材料，应力偏量第二不变量等于某一定值时，材料开始进入了塑性状态。 (4-4) 3：强化法则对理想的弹塑性材料而言，因无强化作用，所以，整个塑性变形过程

33、中，屈服函数值保持一个常量，强化定义了屈服面在应力空间的演化准则。 (4-5) 其中，是强化参数。通常采用的强化法则有以下几种：(1) 各向同性强化此法则规定材料进入塑性变形以后，加载曲面在各方向均匀的向外扩张，没有畸变。而其形状、中心及其在应力空间的方位均保持不变10。需要指出的是：各向同性强化法则主要适用于单调加载情况。如果用于卸载情况，它只适合反向屈服应力等于应力反转点的材料，而通常材料不具备这种性质，因此在塑性力学中还发展了其它强化准则。(2) 随动强化此法则规定材料进入塑性状态以后，加载曲面在应力空间作刚体移动而没有转动，因此初始屈服面的形状、大小和方向仍然保持不变。(3) 混合强化

34、把各向同性强化模型和随动强化模型加以组合，得到混合强化模型。它假定在塑性变形过程中，加载曲面不但作刚性平移，还同时在各个方向作均匀扩大。在以上几种强化模型中，各向同性强化模型应用最为广泛。本文也是采用该硬化法则，这一方面是由于它便于进行数学处理；另一方面，如果在加载过程中应力方向（或各个应力分量的比值）变化不大，采用各向同性强化模型的计算结果及实际情况也比要符合。随动强化模型可以考虑材料的包兴格（）效应，在循环加载或可能出现反向屈服的问题中，需要采用这种模型。 由于塑性变形及变形历史有关， 因此反映塑性应力-应变关系的本构关系用应变增量形式给出比较方便。用应变增量形式表示塑性本构关系的理论称为

35、塑性增量理论。增量理论的本构关系在理论上是合理的，但应用比较麻烦，因为要积分整个变形路径才能得到最后结果。因此，又发展出塑性全量理论，即采用全量应力和全量应变表示塑性本构关系的理论。在比例变形的条件下，可通过积分增量理论的本构关系获得全量理论的本构关系。当偏离比例变形条件不多时，全量理论的计算结果和实险结果比较接近。本文的程序都是基于增量理论。4.2. 非线性有限元算法理论 对于非线性问题，在有限元求解该问题时，对一个自由度总可以表达成，式中，为基本未矢量。如果是线性问题，及无关，而是一次项，显然这是一个线性方程。如果及相关，则方程的出现非一次项，变成非线性问题，在实际工程中，特别是塑性成型问

36、题，材料的几何方程，本构方程以及边界条件往往是非线性的也体现在中出现了，所以变为了非线性问题，要得到最基本的未知量，就必须求解非线性方程组1：直接迭代法又称常刚度法11，这是种最简.单的求解方法，在每次求解前，利用上次的解来求出这一次的值，然后利用和的倒数的乘积求出的当前值 (4-6) 表达为迭代形式 (4-7) 上式可以看出，这种方法首先需要有一个初始的值，以便开始迭代。另外，每一次求解都需要对求倒数，如果求解方程组，就是对刚度矩阵求逆，这种方法在求解中控制两次求解之差，当其值很小时，就认为接近真实值了，迭代结束 图4-32：方法方法的算法及常刚度法不同12，如果得近似表达式是不成立的，存在

37、着残余值，即，此式也可以作为近似值及真实值的差值量度，实际上在具体计算时，也可以控制其值，当极小时，就认为接近真实值了，当第次迭代的值是真实解，则可以按照级数展开得到 图4-43：切线刚度法在复杂非线性问题求解中，刚度及的大小是有一定关系的，在用增量法来求解这种问题时，就等于结构任一点处力及位移的曲线的局部梯度，称为切线刚度13，刚度矩阵的倒数很难用自变量显示表达，通过增量方式求解，在每一步荷载增量范围内把问题线性化，求解方法及方法相同总结以上可以得到： 以上几种算法中，通过比较，不难发现，直接迭代法采用了固定的刚度，适合解决非线性程度不高的本构关系，而切线刚度法采用了变化的刚度，在每一步上都

38、做了实时的修正，对非线性程度较高本构关系任然有效，在效率和迭代精度方面，切线刚度法采用的修正更符合非线性材料的应力应变关系，具有较大的优势，这也是本文采用切线刚度法计算的原因，当然，非线性有限元算法还有很多，切线刚度法也不见得就是最好的能解决所有问题的算法，但是它是在程序开发难度不高和精度方面较高的条件下相对来说最好的 本文采用的本构关系是同性硬化弹塑性模型14，采用屈服准则，下面将根据J2理论15，分别推导常刚度法和切线刚度法计算该问题的的算法公式4.3. 增量理论常刚度法公式推导由应力应变关系得： (4-8) (4-9) 其中是弹性矩阵，它的表达式为是塑性应变增量，它的表达式为其中为等效塑

39、性应变增量，它的表达式为为切线模量对于3维空间问题，流动方向：等效应力 应力偏量 4.4. 增量理论切线刚度法公式推导本文采用的是一种切线刚度法，其应力应变关系为 (4-10) (4-11) 3-14两端左右同乘 ，得到：对于强化材料，准则代入可得其中，为等效塑性应变增量，它的表达式为由于：由流动法则可知 应用上式得所以得到 (4-12) 这就是切线刚度法的矩阵表达式其中为弹塑性矩阵，在里面称为雅可比矩阵，它的表达式为其中是弹性矩阵，它的表达式为为塑性矩阵，它的表达式为对于3维空间问题流动方向 所以最后推导可得，弹塑性矩阵的表达式其中 为切线模量，对本构关系求导得到 总结推导过程：上面的推导过

40、程看似复杂，其实核心的问题只有一个，即两者对于塑性阶段应力更新的算法不同。常刚度法采用的是：而切线刚度法采用的是：把握好了两者的本质上的区别，对于两者的算法设计和程序开发问题便迎刃而解5. 程序设计和编码 本章将严格按照前一章推导的公式展开程序设计和编码，为了便于编程，本文将本构关系做了抽象化处理，即将其描述成一个含参数的表达式，改变参数即可应用于不同的模型，这样做的好处是能保证程序的复用性，这也是本文反复强调的使用的原则。5.1. 本构关系描述 本文采用各向同性硬化弹塑性材料，材料参数如下： (5-1) 弹性部分：，弹性模量200000，泊松比0.3 图5-1 弹性部分本构关系塑性部分：，为

41、了研究方便，取700，0.5，400图5-2 塑性部分本构关系将2个曲线统一到同一个坐标系（为方便显示，x轴标注时扩大了1000倍）图5-3 本构关系由此可以求出两条线的交点即初始屈服点的应力0=400（注意两条曲线相差了一个屈服应变，因为两者其实不是一个坐标系）综上定义的材料常数见表5-1：表5-1弹性模量E200000泊松比0.3屈服应力0400A700B0.5C400注意：上面的取值只是为了便于理解和分析本文的材料模型，为了保证程序的通用性，本文的参数在程序中的一律用变量表示。5.2. 常刚度法程序设计算法设计1：定义程序需要用到的常数和变量2：读取定义的材料常数和状态变量(这里只定义了一个状态变量)，材料常数为，弹性模量E，泊松比，屈服应力0，参数,并且计算出剪切模量G,状态变量为等效塑性应变3：读取应力分量，计算平均应力，应力偏量以及等效应力平均应力： 应力偏量：等效应力：4：根据3计算的等效应力和2读取的屈服应力0比较，如果等效应力小于屈服应力，表明此时材料未屈服，那么转到5，否则转到65：雅可比矩阵,初始化为0，计算弹性矩阵，按照弹性理论更新应力6：雅可比矩阵,初始化为01） ：计算切线模量H,注意到当等效塑性应时对应于本构关系的屈服点，此时的H不能通过上式计算，可以取此时的H为弹性模