ABAQUS子程序UMAT里弹塑本构的实现

摘 要

前 言

有限元法是工程中广泛使用的一种数值计算方法。它是力学、计算方法和计算机技术相结合的产物。在工程应用中,有限元法比其它数值分析方法更流行的一个重要原因在于:相对与其它数值分析方法,有限元法对边界的模拟更灵活,近似程度更高。所以,伴随着有限元理论以及计算机技术的发展,大有限元软件的应用证变得越来越普及。

ABAQUS软件一直以非线性有限元分析软件而闻名,这也是它和ANSYS,Nastran等软件的区别所在。非线性有限元分析的用处越来越大,因为在所用材料非常复杂很多情况下,用线性分析来近似已不再有效。比方说,一个复合材料就不能用传统的线性分析软件包进行分析。任何与时间有关联,有较大位移量的情况都不能用线性分析法来处理。多年前,虽然非线性分析能更适合、更准确的处理问题,但是由于当时计算设备的能力不够强大、非线性分析软件包线性分析功能不够健全,所以通常采用线性处理的方法。

这种情况已经得到了极大的改善,计算设备的能力变得更加强大、类似ABAQUS这样的产品功能日臻完善,应用日益广泛。

非线性有限元分析在各个制造行业得到了广泛应用,有不少大型用户。航空航天业一直是非线性有限元分析的大客户,一个重要原因是大量使用复合材料。新一代波音 787客机将全部采用复合材料。只有像 ABAQUS这样的软件,才能分析包括多个子系统的产品耐久性能。在汽车业,用线性有限元分析来做四轮耐久性分析不可能得到足够准确的结果。分析汽车的整体和各个子系统的性能要求(如悬挂系统等)需要进行非线性分析。在土木工程业, ABAQUS能处理包括混凝土静动力开裂分析以及沥青混凝土方面的静动力分析,还能处理高度复杂非线性材料的损伤和断裂问题,这对于大型桥梁结构,高层建筑的结构分析非常有效。

瞬态、大变形、高级材料的碰撞问题必须用非线性有限元分析来计算。线性分析在这种情况下是不适用的。以往有一些专门的软件来分析碰撞问题,但现在ABAQUS在通用有限元软件包就能解决这些问题。所以,ABAQUS可以在一个软件完成线性和非线性分析。

ABAQUS给用户提供了强大二次开发接口,尤其是在材料本构方面,给用户开发符合实际工程的材料本构模型提供了强大帮助,本文将针对其用户材料子程序展开研究,总结常用材料模型的开发方法。

摘 要

ABAQUS软件功能强大,特别是能够模拟复杂的非线性问题,它包括了多种材料本构关系及失效准则模型,并具有良好的开放性,提供了若干个用户子程序接口,允许用户以代码的形式来扩展主程序的功能。

本文主要研究了ABAQUS用户子程序UMAT的开发方法,采用FORTRAN语言编制了各向同性硬化材料模型的接口程序,研究该类材料的弹塑性本构关系极其实现方法。

本文紧紧围绕UMAT的二次开发技术,首先对其接口原理做了详细介绍,然后针对非线性有限元增量理论中的常刚度法和切线刚度法的算法理论做了深入的剖析,推导出了常刚度法和切线刚度法的算法理论的具体表达式,然后分别编制了两种算法的UMAT程序,最后建立了一个具体的验算模型,通过与ABAQUS自带弹塑性本构关系的计算结果相比较,验证两者的正确性。

本文还对常刚度法和切线刚度法得算法效率做了对比,得出了在非线性程度较高时切线刚度法效率高于常刚度法的结论。

关键字: ABAQUS、UMAT、有限元、材料非线性、FORTRAN、切线刚度

ABSTRACT

ABAQUS software powerful, especially to simulate complex non-linear problem, which includes a wide range of material constitutive model and failure criteria, and has a good open, providing a number of user subroutine interface that allows users to code form to expand the functions of the main program.

This paper studies the user subroutine UMAT of ABAQUS development methods, the use of FORTRAN language isotropic hardening material model of the interface program, studied the effects of such material is extremely elastic-plastic constitutive relation method.

This article UMAT tightly around the secondary development of technology, the first principle of its interface detail, and then for the theory of nonlinear finite element incremental stiffness of the regular tangent stiffness method and the theory of algorithms to do an in-depth analysis of deduced a regular tangent stiffness and rigidity of the law of the specific expression of algorithm theory, and then the preparation of the two algorithms, respectively, of the UMAT program, and finally the establishment of a specific model checking, bringing with ABAQUS elasto-plastic constitutive relation of the calculated results compared to verify the correctness of the two.

This article also often stiffness and tangent stiffness method was to do a comparison of algorithm efficiency is obtained when a higher degree in the non-linear tangent stiffness method more efficient than the conclusions of law often stiffness.

KEY WORDS:ABAQUS、UMAT、Finite element、Material nonlinearity、FORTRAN、Tangent stiffness


1.​  绪论

1.1.​ 课题的研究背景

有限单元法基本思想的提出,可以追溯到克劳夫(R.W.Clough)在1943年的工作[1],他第一次尝试应用定义在三角形区域上的分片连续函数和最小位能原理相结合,来求解St. Venant扭转问题。1960年克劳夫进一步处理了平面弹性问题,并第一次提出了“有限单元法”的名称,使人们开始认识了有限单元法的功效。

四十多年来,随着电子计算机的广泛应用和发展,有限单元法的理论和应用都得到迅速的,持续不断的发展,其应用己由弹性力学平面问题扩展到空间问题、板壳问题,由静力学问题扩展到稳定问题、动力问题和波动问题。分析的对象从弹性材料扩展到塑性、粘弹性、粘塑性和复合材料等,从固体力学扩展到流体力学、传热学等连续介质力学领域。在工程分析中的作用已从分析和校核扩展到优化设计并和计算机辅助设计。

利用有限元软件解决工程和科学问题,是有限元理论应用于工程设计和科学研究实践的主要形式。由于工程设计的巨大市场需要,有限元软件的发展是很迅速的,目前常用的大型有限元软件常见的有Sap2000,ADINA,MSC/NASTRAN,MSC Marc,ANSYS,ABAQUS等,这些软件的共同特点是具有丰富的单元库和求解器,强大而可靠的分析功能,人们利用这些软件解决了很多工程建设和工业产品设计中遇到的问题,取得了巨大的经济技术效益。

由于工程问题的千差万别,不同的用户有不同的专业背景和发展方向,通用软件不免在具体的专业方面有所欠缺,针对这些不足,大部分的通用软件都提供了二次开发功能,以帮助用户减少重复性的编程工作、提高开发起点、缩短研发周期、降低开发成本,并能简化后期维护工作,给用户带来很多方便。基于通用软件平台进行开发,是目前研究的一个重要发展方向。

ABAQUS也提供了若干用户子程序(User Subroutines)接口,它是一个功能非

常强大且适用的分析工具,与命令行的程序格式相比,用户子程序的限制少得多,

从而使用更加灵活方便。针对ABAQUS所提供的本构关系模型种类有限,无法满足

工程应用需要的问题,用户子程序中的用户材料子程序(User-defined Materia

Mechanical Behavior,简称UMAT)接口可以帮助用户定义自己的材料本构模型和算

法,这是ABAQUS的独到之处。由于其操作方便,能被灵活地应用于各个领域中,

尤其受到用户的青睐。

1.2.​ 本文的研究内容和方法

ABAQUS中用户材料子程序UMAT的开发主要解决两方面的问题:本构模型的建立和积分算法的选择。

本文主要研究非线性材料的UMAT实现方法,并重点研究其迭代算法部分,目前,用户材料子程序UMAT的迭代算法主要是常刚度法,常刚度法的优点在于算法原理较简单,程序编写较方便,缺点是当遇到复杂非线性材料时,其迭代次数较多,收敛速度也较慢,在这个情况下,本文采取的是一种迭代次数较少且收敛速度较快的切线刚度法,具体就是采用FORTRAN语言编制了基于Von-Mises模型的接口程序,并采用切线刚度算法,通过与ABAQUS自带本构关系计算的结果相比较,验证其正确性。

本文的研究工作紧紧围绕UMAT的二次开发技术,首先根据有限元方法推导材料非线性问题算法的公式,然后参考UMAT接口规范设计程序的算法流程,继而编写出该程序,最后建立一个具体的本构和具体的模型做测试,验证程序的正确性,在这一过程中,调试是一个非常重要的过程,占用了大量的时间,在调试程序时采用了将中间变量输出到文本的方式,这样能明确跟进迭代过程,发现算法或程序的缺陷。

本文采用的本构关系是经过归纳和抽象的,也就是说本文的程序并不仅仅是只针对某个具体模型和问题,而是针对所有符合抽象出的各向同性硬化材料,这样做的好处是能保证程序的通用性和复用性,避免以后的重复劳动,当然,这也是符合ABAQUS软件设计UMAT接口的宗旨的。

2.​ 基于ABAQUS软件的二次开发

2.1.​ ABAQUS介绍

ABAQUS是一套功能强大的基于有限元法的工程模拟软件[2],其解决问题的 范围从相对简单的线性分析到最富有挑战性的非线性模拟问题。ABAQUS具备十分丰富的、可模拟任意实际形状的单元库。并与之对应拥有各种类型的材料模型库,可以模拟大多数典型工程材料的性能,其中包括金属、橡胶、高分子材料、复合材料、钢筋混凝土、可压缩弹性的泡沫材料以及岩石和土这样的地质材料。 作为通用的模拟分析工具,ABAQUS 不仅能解决结构分析中的问题,还能模拟和研究各种领域中的问题,如热传导、质量扩散、电子元器件的热控制(热一电耦合分析)、声学分析、土壤力学分析(渗流——应力耦合分析)和压电介质力学分析。

ABAQUS为用户提供了广泛的功能,且使用起来又十分简明。最复杂的问题也可以很容易地建立模型[3]。例如复杂的多部件问题可以通过对每个部件定义材料模型和几何形状,然后再把它们组装起来而构成。在大部分模拟分析问题中,甚至在高度非线性问题中,用户也只需要提供结构的几何形状、材料性能、边界条件和荷载工况这样的工程数据就可以进行分析。在非线性分析中,ABAQUS能自动选择合适的荷载增量和收敛精度。不仅能选择这些参数值,而且能在分析过程中不断地调整参数来保证有效地得到高精度的解,很少需用户去定义这些参数。

2.2.​ ABAQUS各模块简介

ABAQUS 有两个主要的分析模块:ABAQUS/Standard 和ABAQUS/Explicit 。ABAQUS/Standard还有两个特殊用途的附加分析模块:ABAQUS/Aqua和ABAQUS/Design。另外,还有ABAQUS 分别与ADAMS/Flex,C-MOLD和Mold flow的接口模块:ABAQUS/ADAMS,ABAQUS/C-MOLD和ABAQUS/ MOLDFLOW。ABAQUS/CAE是完全的ABAQUS工作环境模块,它包括了ABAQUS模型的构造,交互式提交作业、监控作业过程以及评价结果的能力。ABAQUS/Viewer是ABAQUS/CAE的子集,它具有后处理功能,这些模块之间的关系见图2- 1

ABAQUS子程序UMAT里弹塑本构的实现的图1图2-1

ABAQUS/Standard

ABAQUS/Standard是一个通用分析模块,在数值方法上采用有限元方法常用的

隐式积分。它能够求解广泛的线性和非线性问题,包括结构的静态、动态问题、热

力学场和电磁场问题等。对于通常同时发生作用的几何、材料和接触非线性可以采

用自动控制技术处理,也可以由用户自己控制。

ABAQUS/Explicit

ABAQUS/Explicit是一个在数值方法上采用有限元显式积分的特殊模块,它利用对时间的显式积分求解动态有限元方程。它适合于分析诸如冲击和爆炸这样短暂、瞬时的动态问题,同时对高度非线性问题如模拟加工成型过程中接触条件的改变等也非常有效。

ABAQUS/CAE

ABAQUS/CAE是ABAQUS进行有限元分析的前后处理模块,也是建模、分析和后处理的人机交互平台。该模块根据结构的几何图形生成网格,将材料和截面的特性分配到网格上,并施加载荷和边界条件。该模块可以进一步将生成的模型投入到分析模块中进行高效率的后台运行,并对运行情况进行监测,对计算结果进行后处理。ABAQUS/CAE的后处理支持ABAQUS分析模块的所有功能,并且对计算结果的描述和解释提供了范围很广的选择,除了通常的云图,等值线和动画显示之外,还可以用列表,曲线(包括部分常用运算)等其他常用工具来完成对结果数据的处理。该模块的许多独特功能与特点,例如CAD特征化建模、参数化建模、适应设计者要求的数据管理系统等极大的方便了ABAQUS的使用者。

ABAQUS/Aqua

ABAQUS/Aqua的一系列功能可以附加在ABAQUS/Standard中应用。它偏向于模拟海上结构,如海洋石油平台。它的功能包括模拟波浪,风载荷及浮力的 影响。在本指南中不讨论ABAQUS/Aqua。

ABAQUS/ADAMS

ABAQUS/ADAMS允许ABAQUS有限元模型作为柔性部件进入到MDIADAMS产品族中去进行分析。

ABAQUS/C-MOLD

ABAQUS/C-MOLD把注模分析软件C-MOLD中有限元网格、材料性质和初始应力数据转换成为ABAQUS 输入文件。

ABAQUS/Design

ABAQUS/Design 的一系列功能可附加在ABAQUS/Standard 中进行设计敏度计算。

ABAQUS/MOLDFLOW

ABAQUS/MOLDFLOW 模块把MOLDFLOW 分析软件中的有限元模型信息 转换成ABAQUVS 输入文件的一部分。

2.3.​ ABAQUS的二次开发平台

ABAQUS的脚本语言接口非常友好,其自嵌的脚本语言是Python[4],系国际上广泛使用、功能强大、具有良好开放性的一种面向对象程序设计语言。所以,应用Python在ABAQUS中进行二次开发也比较方便,且可移植性强。ABAQUS以基于Python的语法规则向二次开发者提供了许多库函数,这些库函数主要是用来增强ABAQUS的交互式(GUI)操作功能。用户可以通过ABAQUS的交互式(GUI)界面实现分析对象的特征造型、指定材料属性、完成网格剖分和控制、提交并监控分析作业,也可以使用ABAQUS脚本语言越过ABAQUS的交互式(GUI)界面直接高效地向ABAQUS内核提交任务。使用Python可以进行参数化建模,修改交互式建立的模型,还可以一次提交多个作业。

出了脚本语言接口,ABAQUS还为用户提供了功能强大的用户子程序接口(Abaqus User Subroutines ),以帮助用户开发基于ABAQUS内核的程序,常用的用户子程序包括UEL(User subroutine to define an element ,用户单元子程序),UMAT(User subroutine to define a material's mechanical behavior,用户材料子程序 )[5],其中UMAT的使用最为广泛,它主要用于用户开发自己的材料模型,以弥补ABAQUS自带材料模型的不足,帮助用户完成各种材料分析,功能极为强大。

在国外,众多的有限元分析和研究者热衷于使用ABAQUS,一个很重要的原因就在于ABAQUS给用户提供了功能强大,使用方便的二次开发工具和接口,使得用户可以方便的进行富含个性化的有限元建模、分析和后处理,满足特定工程问题的需要。通过用户材料子程序接口,用户可定义任何补充的材料模型,不但任意数量的材料常数都可以作为资料被读取,而且ABAQUS对于任何数量的与解相关的状态变量在每一材料计数点都提供了存储功能,以便在这些子程序中应用。

2.4.​ ABAQUS的二次开发语言

ABAQUS的二次开发语言主要有3种:Python,FORTRAN,C++

Python语言主要用于GUI开发,FORTRAN语言主要用于用户子程序开发,而c++语言主要专注于其他高级开发部分。

本文主要是针对用户子程序的开发,所以采用FORTRAN语言,下面简要介绍一下该语言极其特点:

FORTRAN语言是世界上第一个被正式推广使用的高级语言[6]。它是1954年被提出来的,1956年开始正式使用,至今已有三十多年的历史,但仍历久不衰,它始终是数值计算领域所使用的主要语言。

FORTRAN语言是Formula Translation的缩写,意为“公式翻译”。它是为科学、工程问题或企事业管理中的那些能够用数学公式表达的问题而设计的,其数值计算的功能较强。

FORTRAN语言问世以来,根据需要几经发展,先后推出了不同的版本,主要版本有FORTRAN 77,FORTRAN 90,FORTRAN 95,ABAQUS采用FORTRAN 77,通常用固定格式编写代码。

FORTRAN77语言同C语言一样,是一种结构化编程语言

结构化程序设计方法规定,在结构化的程序中,只能有三种基本结构:

(1)顺序结构

这是一种最简单的基本结构形式,它的特点是,在这个结构内的各个功能模块或语句序列,是按其出现的先后顺序执行的,如赋值语句、输入/输出语句等。它有一个入口和一个出口,并在入口和出口之间包含着若干个功能块,其中每一个功能块可以是一个非转移语句。因此,顺序基本结构块是由一系列的顺序执行语句组成的。

(2)分支选择结构

在给定的条件下,分支选择结构判断选择哪一条路径执行,不同路径完成的功能是不同的。实现分支选择结构主要由块IF语句、ELSE语句、END IF语句以及ELSE IF语句组成的IF-THEN-ELSE结构。

(3)循环结构

循环结构也称重复处理结构,即重复执行某一功能块,直到满足(或不满足)某一条件为止。实现循环结构的FORTRAN90语句主要是DO语句、块IF语句和逻辑IF语句的结合。

以上三种基本结构,是组成结构化程序的基本结构形式。这里有两层意思:一是结构化的程序中,各个模块均由这三种基本结构组成;二是结构化程序本身,从宏观上也是这三种基本结构形式之一。

3.​  用户材料子程序UMAT

3.1.​ UMAT开发环境设置

由于UMAT是采用FORTRAN语言编写,那么要运行UMAT就需要安装FORTRAN的开发环境, 同时还需要ABAQUS的支持,本文采用的ABAQUS版本为6.81,支持INTEL Fortran9.1-10.1,Intel Fortran安装时又需要安装Microsoft Visual Studio的相应版本,经过比较,本文选用ABAQUS6.81+Intel Fortran10.1+Microsoft VisualC++ 2005,相对于ABAQUS来说,UMAT开发环境的设置较为繁琐,这给子程序的使用带来诸多不便,为了解决这一问题,我用C#语言编制了ABAQUS子程序编译环境设置工具,只需要将安装文件解压到ABAQUS的安装目录,运行安装程序就可以了,整个过程不需要人工干预,也不需要安装庞大的VisualC++ 2005,如图3-1所示

ABAQUS子程序UMAT里弹塑本构的实现的图2 图3-1

3.2.​ UMAT注意事项

ABAQUS的用户子程序是根据ABAQUS提供的相应接口,按照Fortran语法,

用户自己编写的代码。它是一个独立的程序单元,可以独立的被存储和编译,也能

被其它程序单元引用,因此,利用它可带回大量数据供引用程序使用,也可以用它

来完成各种特殊的功能。它的一般结构形式是:

SUBROUTINE S(x1,x2,……,xn)

INCLUDE‘ABA_PARAM.INC’(用于ABAQUS/Standard用户子程序中)

OR INCLUDE‘VABA_PARAM.INC’)(用于ABAQUS/Explicit用户子程序中)

……

RETURN

END

x1,x2,……,xn是ABAQUS提供的用户子程序的接口参数,有些参数是ABAQUS传到用户子程序中的,例如SUBROUTINE DLOAD中的KSTEP、KINC、COORDS,有些是需要用户自己定义的,例如F,文件aba_param.inc和vaba_param.inc随着ABAQUS软件的安装而包含在操作系统中,它们含有重要的参数,帮助ABAQUS主求解程序对用户子程序进行编译和链接。当控制遇到RETURN语句时便返回到引用程序单元中去,END语句是用户子程序结束的标志。

在一个算例中,用户可以用到多个用户子程序,但必须把它们放在一个以.for为扩展名的文件中。运行带有用户子程序的算例同时有两种方法:一是在CAE中运行,在EDIT JOB菜单中的GENERAL子菜单的USER SUBROUTINE FILE对话框中选择用户子程序所在的文件即可;另外是在ABAQUS.COMMAND中运行,语法如下:

abaqus job=job-name user={source-file|object-file}

编制用户子程序时应注意

(1)用户子程序相互之间不能调用,但可以调用用户自己编写的Fortran子程序和ABAQUS应用程序。ABAQUS应用程序必须由用户子程序调用。当用户编写Fortran子程序时,建议子程序名以K开头,以免和ABAQUS内部程序冲突。

(2)当用户在用户子程序中利用OPEN打开外部文件时,要注意以下两点:一是设备号的选择是有限制的,只能取15~18和大于100的设备号,其余的都已被ABAQUS占用;二是用户需提供外部文件的绝对路径而不是相对路径。

(3)对于不同的用户子程序ABAQUS调用的时间是不同的,有的是在每个STEP

的开始,有的是STEP的结尾,有的是在每个INCREMENT的开始等等。当ABAQUS

调用用户子程序时,都会把当前的STEP和INCREMENT利用用户子程序的两个实

参KSTEP和KINC传给用户子程序,用户可把它们输出到外部文件中,这样就可清

楚的知道ABAQUS何时调用该用户子程序。

为保证用户子程序的正确执行,子程序的书写必须遵循ABAQUS的相关规定,

下面以用户材料子程序为例详细说明。

3.3.​ UMAT接口的原理

用户材料子程序(User-defined Material Mechanical Behavior,简称UMAT)是

ABAQUS提供给用户定义自己的材料属性的Fortran程序接口[7][8],使用户能使用

ABAQUS材料库中没有定义的材料模型。用户材料子程序UMAT通过与ABAQUS主求解程序的接口实现与ABAQUS的资料交流。在输入文件中,使用关键词“*USER MATERIAL”表示定义用户材料属性。

UMAT子程序具有强大的功能,使用UMAT子程序:

(1)可以定义材料的本构关系,使用ABAQUS材料库中没有包含的材料进行

计算,扩充程序功能。

(2)几乎可以用于力学行为分析的任何分析过程,几乎可以把用户材料属性赋

予ABAQUS中的任何单元。

(3)必须在UMAT中提供材料本构的雅可比(Jacobian)矩阵,即应力增量对

应变增量的变化率。

由于主程序与UMAT之间存在数据传递,甚至共享一些变量,因此必须遵守有

关UMAT的书写格式,UMAT中常用的变量在文件开头予以定义,通常格式

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

1 RPL,DDSDDT,DRPLDE,DRPLDT,

2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,

3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,

4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

INCLUDE‘ABA_PARAM.INC’

CHARACTER*80 CMNAME

DIMENSION STRESS(NTENS),STATEV(NSTATV),

1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),

3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)

user coding to define DDSDDE,STRESS,STATEV,SSE,SPD,SCD

and,if necessary,RPL,DDSDDT,DRPLDE,DRPLDT,PNEWDT

RETURN

END

UMAT中的应力矩阵、应变矩阵以及矩阵DDSDDE、DDSDDT、DRPLDE等,都是直接分量存储在前,剪切分量存储在后。直接分量有NDI个,剪切分量有NSHR个。各分量之间的顺序根据单元自由度的不同有一些差异,所以编写UMAT时要考虑到所使用单元的类别。下面对UMAT中用到的一些变量进行说明:

ABAQUS子程序UMAT里弹塑本构的实现的图3DDSDDE(NTENS NTENS):一个NTENS×NTENS的矩阵,称作Jacobian矩阵

ABAQUS子程序UMAT里弹塑本构的实现的图4是应力的增量,ABAQUS子程序UMAT里弹塑本构的实现的图5是应变的增量,DDSDDE(i,j)表示增量步结束时第j个应变分量的改变引起的第i个应力分量的变化。通常Jacobian矩阵是一个对称矩阵,除非在“*USER MATERIAL”语句中加入了“UNSYMM”参数。

STRESS(NTENS):应力张量数组,对应NDI个直接分量和NSHR个剪切分量。在增量步的开始,应力张量矩阵中的数值通过UMAT和主程序之间的接口传递到UMAT中,在增量步的结束UMAT将对应力张量矩阵更新。对于包含刚体转动的有限应变问题,一个增量步调用UMAT之前就已经对应力张量进行了刚体转动,因此UMAT中只需处理应力张量的共旋部分。UMAT中应力张量的度量为柯西应力

STATEV(NSTATEV):用于存储与解有关的状态变量的数组,在增量步开始时将数值传递到UMAT中,也可在子程序USDFLD或UEXPAN中先更新数据,然后增量步开始时将更新后的资料传递到UMAT中。在增量步的结束必须更新状态变量矩阵中的数据。和应力张量矩阵不同的是:对于有限应变问题,除了材料本构行为引起的资料更新以外,与解有关的状态变量矩阵中的任何向量或者张量都必须通过旋转来考虑材料的刚体运动。状态变量矩阵的维数通过ABAQUS输入文件中的关键词“*DEPVAR”定义,关键词下面数据行的数值即为状态变量矩阵的维数。

PROPS(NPROPS):材料常数数组。材料常数的个数,等于关键词“*USER MATERIAL”中“CONSTANTS”常数设定的值。矩阵中元素的数值对应于关键词“USER MATERIAL”下面的数据行。

SSE,SPD,SCD:分别定义每一增量步的弹性应变能,塑性耗散和蠕变耗散。

它们对计算结果没有影响,仅仅作为能量输出。

STRAN(NTENS):应变数组。

DSTRAN(NTENS):应变增量数组。

DTIME:增量步的时间增量。

NDI:直接应力分量的个数。

NSHR:剪切应力分量的个数。

NTENS:总应力分量的个数,NTENS=NDI+NSHR。

由于UMAT子程序在单元的积分点上调用,增量步开始时,主程序路径将通过

UMAT的接口进入UMAT,单元当前积分点必要变量的初始值将随之传递给UMAT

的相应变量。在UMAT结束时,变量的更新值将通过接口返回主程序。

3.4.​ UMAT的使用方法

我们知道,有限元计算(增量方法)的基本问题[7]是:已知第n 步的结果(应力,应变等)ABAQUS子程序UMAT里弹塑本构的实现的图6ABAQUS子程序UMAT里弹塑本构的实现的图7,然后给出一个应变增量ABAQUS子程序UMAT里弹塑本构的实现的图8,计算ABAQUS子程序UMAT里弹塑本构的实现的图9,UMAT要完成这一计算,并要计算DDSDDE(I,J)=ABAQUS子程序UMAT里弹塑本构的实现的图10ABAQUS子程序UMAT里弹塑本构的实现的图11是应力增量矩阵,ABAQUS子程序UMAT里弹塑本构的实现的图12是应变增量矩阵,DDSDDE(I,J) 定义了第J 个应变分量的微小变化对第 I 个应力分量带来的变化。 该矩阵只影响收敛速度,不影响计算结果的准确性 (当然,不收敛自然得不到结果)。

有限元计算的中心问题就是求得节点的位移 (进而应变、应力),以使内力和

外力达到平衡:

ABAQUS子程序UMAT里弹塑本构的实现的图13 (3-1)

d 是节点位移矩阵,黑体字表示矩阵或矢量。除了小变形、线弹性问题,方程2-1是非性的,要用迭代的方法解出:

ABAQUS子程序UMAT里弹塑本构的实现的图14 (3-2)

ABAQUS子程序UMAT里弹塑本构的实现的图15 (3-3)

i 表示一个增量步内的第i 次迭代,n表示第n个增量步。KT是切线刚度,由材料的Jacobian 矩阵结合单元计算组装而得。刚度矩阵其实就是力对位移的梯度。要想快速收敛,位移增量应沿该梯度方向变化,也就是说,如果Jacobian 矩阵不是那么准确,自然KT 也不怎么准确,那么满足3-1式的位移被找到的速度也就变慢,甚至发散,根本找不到。但收敛速度无论慢快,3-1式才是判断结果准确与否的唯一标准。所以Jacobian 矩阵不影响结果的准确性,只影响收敛速度的快慢。

以大变形、非性材料为例,整个计算步骤是这样的:

整个外力不是一次加上,而是一点点加上的,不然会发散得不到结果的。所以,每一个增量步开始时就是在原来的外力上加上一点点,得到ABAQUS子程序UMAT里弹塑本构的实现的图16。根据3-2得到位移增量ABAQUS子程序UMAT里弹塑本构的实现的图17,此时要知道力对位移的梯度KT,以尽快找到满足平衡条件的位移,由材料的Jacobian 矩阵和单元结合起来组装得到(此处使用UMAT 提供的Jacobian 矩阵)。然后可计算应变增量ABAQUS子程序UMAT里弹塑本构的实现的图18,调用UMAT,得到新的应力,进而得到新的内力,所以,程序不在乎新的应力是由增量方法得到,还是全量方法得到,而只在乎新应力是否准确。然后回到3-2,如此循环,直至3-2右端为0,也即满足3-1。这样第n+1 步就完成了,然后开始第n+2 步,即 外力加上一点点,按同样的方法求解新的位移。直至整个外力全部施加并得到满足3-1的位移。

4.​  材料非线性问题

弹性力学作为精确理论,从本质上都是非线性的,早期Cauchy, Green,Kirchhoff和Kelvin在这些方面都作出了重要贡献。后来又提出超弹性(即具有弹性势的)有限变形理论,由于理论方程的冗长而复杂,且工程应用也没有提出这方面要求而被搁置。20世纪40年代以后,由于橡胶材料、高分子合成材料的迅速发展和工业领域的大量应用,非线性弹性与超弹性的研究再次引起科学和工程界的重视,除了一般理论研究有了新的发展以外,工程应用计算方法也得到长足的发展。

4.1.​ 材料的弹塑性本构关系

弹塑性材料进入塑性的特征是当荷载卸去后存在不可恢复的永久变形。所以,在卸载情况下,应力应变之间不再是唯一的对应关系。这是区别于非线性弹性材料的基本属性。只以加载时应力应变关系成非线性,还不足以判断材料是非线性弹性还是弹塑性。但是一经卸载就可以看出两者的区别。非线性弹性材料沿原路径返回,而弹塑性材料将依据不同的加载历史卸载后产生不同的永久变形。

对大多数材料来说,在单调加载的情况下,存在一个明显的极限应力ABAQUS子程序UMAT里弹塑本构的实现的图19,当应力低于ABAQUS子程序UMAT里弹塑本构的实现的图20时,材料保持线弹性。而当应力达到ABAQUS子程序UMAT里弹塑本构的实现的图21以后,则材料开始进入弹塑性状态。如继续加载,然后在卸载,材料始终保持永久的塑性变形。如果应力达到ABAQUS子程序UMAT里弹塑本构的实现的图22后,应力不再增加,而材料变形可以继续增加,及变形处于不定的流动状态,则称材料为理想弹塑性的。反之如果应力达到ABAQUS子程序UMAT里弹塑本构的实现的图23后,再增加变形,应力也必须增加,则材料是应变硬化的,这时应力ABAQUS子程序UMAT里弹塑本构的实现的图24是塑性应变ABAQUS子程序UMAT里弹塑本构的实现的图25的函数,可解析为:

ABAQUS子程序UMAT里弹塑本构的实现的图26 (4-1)

本构关系反应着应力应变之间的关系。对于弹性材料变形是可以恢复的;而塑性材料变形是不可以恢复的。典型的弹塑性应变在卸载后要保持一个永久的变形。如图3-2

ABAQUS子程序UMAT里弹塑本构的实现的图27

图4-1

塑性应变有下列特性:

(1)总应变分为弹性和塑性两部分,即

ABAQUS子程序UMAT里弹塑本构的实现的图28 (4-2)

或者:

ABAQUS子程序UMAT里弹塑本构的实现的图29 (4-3)

(2)塑性变形取决于加载路径,而应力应变之间没有一一对应的关系。所以必须确定二则之间的本构关系,这种本构关系可以用偏微分方程或者增量形式来描述。

总之,弹塑性理论主要包括以下几个方面:

(1)应变张量的分解;

(2)应力空间的屈服条件;

(3)流动法则;

(4)强化法则;

(5)协调性条件。

1:本构模型

塑性力学的应力-应变曲线通常有5种简化模型[8]

(1)理想弹塑性模型,用于低碳钢或强化性质不明显的材料。

(2)线性强化弹塑性模型,用于有显著强化性质的材料。

(3)理想刚塑性模型,用于弹性应变比塑性应变小得多且强化性质不明显的材料。

(4)线性强化刚塑性模型,用于弹性应变比塑性应变小得多且强化性质明显的材料。

(5)幂强化模型,为简化计算中的解析式,可将应力-应变关系的解析式写为

σ=σy(ε/εy)n,式中σy为屈服应力,εy为与σy相对应的应变,n为材料常数。

ABAQUS子程序UMAT里弹塑本构的实现的图30

图4-2

2:屈服条件

在复杂应力状态下,判断物体屈服状态的准则称为屈服条件[9]。屈服条件是各应力分量组合应满足的条件。对于金属材料,最常用的屈服条件为最大剪应力屈服条件(又称Tresca屈服条件)和弹性形变比能屈服条件(又称Von Mises条件)。对于岩土材料则常用Tresca屈服条件、Drucker-Prager屈服条件和Mohr-Coulomb屈服条件。对于强化或软化材料,屈服条件将随塑性变形的增长而变化,改变后的屈服条件称为后继屈服条件。当已知主应力的大小次序时,使用Tresca屈服条件较为方便;若不知道主应力的大小次序,则使用Von Mises屈服条件较为方便。对于韧性较好的材料,Von Mises屈服条件与试验数据符合较好。

Von Mises屈服准则具体形式是,对于各项同性材料,应力偏量第二不变量ABAQUS子程序UMAT里弹塑本构的实现的图31等于某一定值时,材料开始进入了塑性状态。

ABAQUS子程序UMAT里弹塑本构的实现的图32 (4-4)

3:强化法则

对理想的弹塑性材料而言,因无强化作用,所以,整个塑性变形过程中,屈服函数值保持一个常量,强化定义了屈服面在应力空间的演化准则。

ABAQUS子程序UMAT里弹塑本构的实现的图33 (4-5)

其中,ABAQUS子程序UMAT里弹塑本构的实现的图34是强化参数。

通常采用的强化法则有以下几种:

(1) 各向同性强化

此法则规定材料进入塑性变形以后,加载曲面在各方向均匀的向外扩张,没有畸变。而其形状、中心及其在应力空间的方位均保持不变[10]。需要指出的是:各向同性强化法则主要适用于单调加载情况。如果用于卸载情况,它只适合反向屈服应力等于应力反转点的材料,而通常材料不具备这种性质,因此在塑性力学中还发展了其它强化准则。

(2) 随动强化

此法则规定材料进入塑性状态以后,加载曲面在应力空间作刚体移动而没有转动,因此初始屈服面的形状、大小和方向仍然保持不变。

(3) 混合强化

把各向同性强化模型和随动强化模型加以组合,得到混合强化模型。它假定在塑性变形过程中,加载曲面不但作刚性平移,还同时在各个方向作均匀扩大。

在以上几种强化模型中,各向同性强化模型应用最为广泛。本文也是采用该硬化法则,这一方面是由于它便于进行数学处理;另一方面,如果在加载过程中应力方向(或各个应力分量的比值)变化不大,采用各向同性强化模型的计算结果与实际情况也比要符合。随动强化模型可以考虑材料的包兴格(Bauschinger)效应,在循环加载或可能出现反向屈服的问题中,需要采用这种模型。

由于塑性变形与变形历史有关, 因此反映塑性应力-应变关系的本构关系用应变增量形式给出比较方便。用应变增量形式表示塑性本构关系的理论称为塑性增量理论。增量理论的本构关系在理论上是合理的,但应用比较麻烦,因为要积分整个变形路径才能得到最后结果。因此,又发展出塑性全量理论,即采用全量应力和全量应变表示塑性本构关系的理论。在比例变形的条件下,可通过积分增量理论的本构关系获得全量理论的本构关系。当偏离比例变形条件不多时,全量理论的计算结果和实险结果比较接近。本文的程序都是基于增量理论。

4.2.​ 非线性有限元算法理论

对于非线性问题,在有限元求解该问题时,对一个自由度总可以表达成ABAQUS子程序UMAT里弹塑本构的实现的图35,式中,ABAQUS子程序UMAT里弹塑本构的实现的图36为基本未矢量。如果是线性问题,ABAQUS子程序UMAT里弹塑本构的实现的图37ABAQUS子程序UMAT里弹塑本构的实现的图38无关,而ABAQUS子程序UMAT里弹塑本构的实现的图39是一次项,显然这是一个线性方程。如果ABAQUS子程序UMAT里弹塑本构的实现的图40ABAQUS子程序UMAT里弹塑本构的实现的图41相关,则方程的ABAQUS子程序UMAT里弹塑本构的实现的图42出现非一次项,变成非线性问题,在实际工程中,特别是塑性成型问题,材料的几何方程,本构方程以及边界条件往往是非线性的也体现在ABAQUS子程序UMAT里弹塑本构的实现的图43中出现了ABAQUS子程序UMAT里弹塑本构的实现的图44,所以变为了非线性问题,要得到最基本的未知量,就必须求解非线性方程组

1:直接迭代法

又称常刚度法[11],这是种最简.单的求解方法,在每次求解前,利用上次的ABAQUS子程序UMAT里弹塑本构的实现的图45解来求出这一次的ABAQUS子程序UMAT里弹塑本构的实现的图46值,然后利用ABAQUS子程序UMAT里弹塑本构的实现的图47ABAQUS子程序UMAT里弹塑本构的实现的图48的倒数的乘积求出ABAQUS子程序UMAT里弹塑本构的实现的图49的当前值

ABAQUS子程序UMAT里弹塑本构的实现的图50 (4-6)

表达为迭代形式

ABAQUS子程序UMAT里弹塑本构的实现的图51 (4-7)

上式可以看出,这种方法首先需要有一个初始的ABAQUS子程序UMAT里弹塑本构的实现的图52值,以便开始迭代。另外,每一次求解都需要对ABAQUS子程序UMAT里弹塑本构的实现的图53求倒数,如果求解方程组,就是对刚度矩阵求逆,这种方法在求解中控制两次求解ABAQUS子程序UMAT里弹塑本构的实现的图54之差,当其值很小时,就认为接近真实值了,迭代结束

ABAQUS子程序UMAT里弹塑本构的实现的图55 图4-3

2:Newton-Raphson方法

Newton-Raphson方法的算法与常刚度法不同[12],如果ABAQUS子程序UMAT里弹塑本构的实现的图56得近似表达式ABAQUS子程序UMAT里弹塑本构的实现的图57是不成立的,存在着残余值,即ABAQUS子程序UMAT里弹塑本构的实现的图58,此式也可以作为近似值与真实值的差值量度,实际上在具体计算时,也可以控制其值,当ABAQUS子程序UMAT里弹塑本构的实现的图59极小时,就认为ABAQUS子程序UMAT里弹塑本构的实现的图60接近真实值了,当第ABAQUS子程序UMAT里弹塑本构的实现的图61次迭代的ABAQUS子程序UMAT里弹塑本构的实现的图62

ABAQUS子程序UMAT里弹塑本构的实现的图63是真实解,则可以按照Taylor级数展开得到

ABAQUS子程序UMAT里弹塑本构的实现的图64

ABAQUS子程序UMAT里弹塑本构的实现的图65

图4-4

3:切线刚度法

在复杂非线性问题求解中,刚度ABAQUS子程序UMAT里弹塑本构的实现的图66ABAQUS子程序UMAT里弹塑本构的实现的图67的大小是有一定关系的,在用增量法来求解这种问题时,ABAQUS子程序UMAT里弹塑本构的实现的图68就等于结构任一点处力与位移的曲线的局部梯度,称为切线刚度[13],刚度矩阵的倒数很难用自变量显示表达,通过增量方式求解,在每一步荷载增量范围内把问题线性化,求解方法与Newton-Raphson方法相同

总结以上可以得到:

以上几种算法中,通过比较,不难发现,直接迭代法采用了固定的刚度,适合解决非线性程度不高的本构关系,而切线刚度法采用了变化的刚度,在每一步上都做了实时的修正,对非线性程度较高本构关系任然有效,在效率和迭代精度方面,切线刚度法采用的修正更符合非线性材料的应力应变关系,具有较大的优势,这也是本文采用切线刚度法计算的原因,当然,非线性有限元算法还有很多,切线刚度法也不见得就是最好的能解决所有问题的算法,但是它是在程序开发难度不高和精度方面较高的条件下相对来说最好的

本文采用的本构关系是同性硬化弹塑性模型[14],采用Mises屈服准则,下面将根据J2理论[15],分别推导常刚度法和切线刚度法计算该问题的的算法公式

4.3.​ 增量理论常刚度法公式推导

由应力应变关系得:

ABAQUS子程序UMAT里弹塑本构的实现的图69 (4-8)

ABAQUS子程序UMAT里弹塑本构的实现的图70 (4-9)

其中ABAQUS子程序UMAT里弹塑本构的实现的图71是弹性矩阵,它的表达式为

ABAQUS子程序UMAT里弹塑本构的实现的图72

ABAQUS子程序UMAT里弹塑本构的实现的图73是塑性应变增量,它的表达式为

ABAQUS子程序UMAT里弹塑本构的实现的图74

其中ABAQUS子程序UMAT里弹塑本构的实现的图75为等效塑性应变增量,它的表达式为

ABAQUS子程序UMAT里弹塑本构的实现的图76

ABAQUS子程序UMAT里弹塑本构的实现的图77

ABAQUS子程序UMAT里弹塑本构的实现的图78

ABAQUS子程序UMAT里弹塑本构的实现的图79为切线模量

对于3维空间问题,流动方向:

ABAQUS子程序UMAT里弹塑本构的实现的图80

等效应力

ABAQUS子程序UMAT里弹塑本构的实现的图81

应力偏量

ABAQUS子程序UMAT里弹塑本构的实现的图82

4.4.​ 增量理论切线刚度法公式推导

本文采用的是一种切线刚度法,其应力应变关系为

ABAQUS子程序UMAT里弹塑本构的实现的图83 (4-10)

ABAQUS子程序UMAT里弹塑本构的实现的图84 (4-11)

3-14两端左右同乘 ABAQUS子程序UMAT里弹塑本构的实现的图85,得到:

ABAQUS子程序UMAT里弹塑本构的实现的图86

对于强化材料,Mises准则

ABAQUS子程序UMAT里弹塑本构的实现的图87

代入可得

ABAQUS子程序UMAT里弹塑本构的实现的图88

其中,ABAQUS子程序UMAT里弹塑本构的实现的图89为等效塑性应变增量,它的表达式为

ABAQUS子程序UMAT里弹塑本构的实现的图90

由于:

ABAQUS子程序UMAT里弹塑本构的实现的图91

ABAQUS子程序UMAT里弹塑本构的实现的图92

由流动法则可知

ABAQUS子程序UMAT里弹塑本构的实现的图93

应用上式得

ABAQUS子程序UMAT里弹塑本构的实现的图94

所以得到

ABAQUS子程序UMAT里弹塑本构的实现的图95 (4-12)

这就是切线刚度法的矩阵表达式

其中ABAQUS子程序UMAT里弹塑本构的实现的图96为弹塑性矩阵,在ABAQUS里面称为雅可比矩阵,它的表达式为ABAQUS子程序UMAT里弹塑本构的实现的图97

其中ABAQUS子程序UMAT里弹塑本构的实现的图98是弹性矩阵,它的表达式为

ABAQUS子程序UMAT里弹塑本构的实现的图99

ABAQUS子程序UMAT里弹塑本构的实现的图100为塑性矩阵,它的表达式为

ABAQUS子程序UMAT里弹塑本构的实现的图101

对于3维空间问题

流动方向

ABAQUS子程序UMAT里弹塑本构的实现的图102

所以

ABAQUS子程序UMAT里弹塑本构的实现的图103

最后推导可得,弹塑性矩阵的表达式

ABAQUS子程序UMAT里弹塑本构的实现的图104

其中

ABAQUS子程序UMAT里弹塑本构的实现的图105 ABAQUS子程序UMAT里弹塑本构的实现的图106

ABAQUS子程序UMAT里弹塑本构的实现的图107为切线模量,对本构关系求导得到

ABAQUS子程序UMAT里弹塑本构的实现的图108

总结推导过程:上面的推导过程看似复杂,其实核心的问题只有一个,即两者对于塑性阶段应力更新的算法不同。

常刚度法采用的是:

ABAQUS子程序UMAT里弹塑本构的实现的图109

而切线刚度法采用的是:

ABAQUS子程序UMAT里弹塑本构的实现的图110

把握好了两者的本质上的区别,对于两者的算法设计和程序开发问题便迎刃而解

5.​  UMAT程序设计和编码

本章将严格按照前一章推导的公式展开程序设计和编码,为了便于编程,本文将本构关系做了抽象化处理,即将其描述成一个含参数的表达式,改变参数即可应用于不同的模型,这样做的好处是能保证程序的复用性,这也是本文反复强调的使用UMAT的原则。

5.1.​ 本构关系描述

本文采用各向同性硬化弹塑性材料,材料参数如下:

ABAQUS子程序UMAT里弹塑本构的实现的图111 (5-1)

弹性部分:ABAQUS子程序UMAT里弹塑本构的实现的图112,弹性模量E=200000Mpa,泊松比Mu=0.3

ABAQUS子程序UMAT里弹塑本构的实现的图113

图5-1 弹性部分本构关系

塑性部分:ABAQUS子程序UMAT里弹塑本构的实现的图114,为了研究方便,取A=700,B=0.5,C=400

ABAQUS子程序UMAT里弹塑本构的实现的图115

图5-2 塑性部分本构关系

将2个曲线统一到同一个坐标系(为方便显示,x轴标注时扩大了1000倍)

ABAQUS子程序UMAT里弹塑本构的实现的图116

图5-3 本构关系

由此可以求出两条线的交点即初始屈服点的应力Yield0=400Mpa(注意两条曲线相差了一个屈服应变,因为两者其实不是一个坐标系)

综上定义的材料常数见表5-1:

表5-1

弹性模量E

200000

泊松比Mu

0.3

屈服应力Yield0

400

A

700

B

0.5

C

400

注意:上面A,B,C的取值只是为了便于理解和分析本文的材料模型,为了保证程序的

通用性,本文的参数在程序中的A,B,C一律用变量表示。

5.2.​ 常刚度法程序设计

算法设计

1:定义程序需要用到的常数和变量

2:读取ABAQUS定义的材料常数和状态变量(这里只定义了一个状态变量),材料常数为,弹性模量E,泊松比Mu,屈服应力Yield0,参数A,B,C,并且计算出剪切模量G,状态变量为等效塑性应变EQPLAS

3:读取应力分量,计算平均应力,应力偏量以及Mises等效应力

平均应力:ABAQUS子程序UMAT里弹塑本构的实现的图117

应力偏量:ABAQUS子程序UMAT里弹塑本构的实现的图118

Mises等效应力:ABAQUS子程序UMAT里弹塑本构的实现的图119

4:根据3计算的Mises等效应力和2读取的屈服应力Yield0比较,如果Mises等效应力小于屈服应力,表明此时材料未屈服,那么转到5,否则转到6

5:雅可比矩阵,初始化为0,计算弹性矩阵,按照弹性理论更新应力

ABAQUS子程序UMAT里弹塑本构的实现的图120

ABAQUS子程序UMAT里弹塑本构的实现的图121

6:雅可比矩阵,初始化为0

1:计算切线模量H'

ABAQUS子程序UMAT里弹塑本构的实现的图122,注意到当等效塑性应ABAQUS子程序UMAT里弹塑本构的实现的图123时对应于本构关系的屈服点,此时的H不能通过上式计算,可以取此时的H为弹性模量

2:计算等效塑性应变增量并更新

ABAQUS子程序UMAT里弹塑本构的实现的图124

EQPLAS=DEQPLAS+EQPLAS更新状态变量

3)计算流动方向

ABAQUS子程序UMAT里弹塑本构的实现的图125

4计算塑性应变增量

ABAQUS子程序UMAT里弹塑本构的实现的图126

5)更新应力

ABAQUS子程序UMAT里弹塑本构的实现的图127

算法流程图

ABAQUS子程序UMAT里弹塑本构的实现的图128

图5-4 常刚度法算法流程图

5.3.​ 常刚度法程序编码

根据算法流程,用FORTRAN77 固定格式编制了常刚度法的计算程序如下:

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,

1 DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,

2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,

3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

INCLUDE 'ABA_PARAM.INC'

CHARACTER*8 CMNAME

DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),

1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),

2 PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),

3 DFGRD0(3,3),DFGRD1(3,3),DPLAS(6)

C---------------------------- 常刚度法

CCCCCCCCCCCCCCCCCCCCCCCCCC1:定义变量

c 定义常数

PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)

C 定义材料常数

DOUBLE PRECISION E,Mu,Yield0,A,B,C

C 定义中间变量

DOUBLE PRECISION H,W,CEGMA_EQ,CEGMA_X,CEGMA_Y,CEGMA_Z,TAO_XY,

1TAO_YZ,TAO_ZX,CEGMA_CP,CEGMA_SX,CEGMA_SY,CEGMA_SZ,TAO_SXY,TAO_SYZ

2,TAO_SZX

C 定义状态变量

DOUBLE PRECISION EQPLAS

C 定义更新变量

DOUBLE PRECISION DEQPLAS

CCCCCCCCCCCCCCCCCCCCCCCCCC2:读取变量

C 弹性模量

E=PROPS(1)

C 泊松比

Mu=PROPS(2)

C 剪切模量

G=E/TWO/(ONE+Mu)

C 屈服应力

Yield0=PROPS(3)

C 参数

A=PROPS(4)

B=PROPS(5)

C=PROPS(6)

C 读取状态参数

EQPLAS=STATEV(1)

Yield=A*(EQPLAS+Yield0/E)**B+C

CCCCCCCCCCCCCCCCCCCCCCCCCC3:读取应力分量,计算平均应力,应力偏量以及Mises等效应力

C 6个应力分量

CEGMA_X= STRESS(1)

CEGMA_Y= STRESS(2)

CEGMA_Z= STRESS(3)

TAO_XY= STRESS(4)

TAO_ZX= STRESS(5)

TAO_YZ= STRESS(6)

C 平均应力

CEGMA_CP=(CEGMA_X+CEGMA_Y+CEGMA_Z)/THREE

C 6个应力偏量

CEGMA_SX=CEGMA_X-CEGMA_CP

CEGMA_SY=CEGMA_Y-CEGMA_CP

CEGMA_SZ=CEGMA_Z-CEGMA_CP

TAO_SXY=TAO_XY

TAO_SZX=TAO_ZX

TAO_SYZ=TAO_YZ

c Mises等效应力

CEGMA_EQ=SQRT(THREE/TWO*(CEGMA_SX**2+CEGMA_SY**2+CEGMA_SZ**2+

1TWO*(TAO_SXY**2+TAO_SYZ**2+TAO_SZX**2)))

C 雅可比矩阵,初始化默认为

DO k1=1,NTENS

DO k2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO

c 计算弹性矩阵

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K1,k2)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)*(Mu/(ONE-Mu))

END DO

DDSDDE(K1,K1)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)*(ONE-TWO*Mu)/TWO/

1(ONE-Mu)

END DO

CCCCCCCCCCCCCCCCCCCCCCCCCC4:根据计算的Mises等效应力和读取的屈服应力Yield0比较,

CCCCCCCCCCCCCCCCCCCCCCCCCC如果Mises等效应力小于屈服应力,表明此时材料未屈服,那么转到,否则转到

IF(CEGMA_EQ.LT.Yield0) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCC5:没有屈服,按照弹性理论计算

C 按照弹性理论更新应力

DO K1=1,NTENS

DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

END DO

END DO

ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCC6:屈服发生,按照塑性理论计算

C 切线模量H,根据本构关系求导

IF(EQPLAS.EQ.0) THEN

H=E

ELSE

H=A*B*EQPLAS**(B-ONE)

END IF

C 求W

W=9.D0*G/TWO/CEGMA_EQ**2/(H+THREE*G)

C 等效塑性应变增量

DEQPLAS=(6.D0*G*CEGMA_EQ/(TWO*CEGMA_EQ**2*H+9.D0*G*(CEGMA_SX

1**2+CEGMA_SY**2+CEGMA_SZ**2+TWO*TAO_XY**2+TWO*TAO_ZX**2+TWO*TAO_YZ

2**2)))*(CEGMA_SX*DSTRAN(1)+CEGMA_SY*DSTRAN(2)+CEGMA_SZ*

3DSTRAN(3)+TAO_XY*DSTRAN(4)+TAO_ZX*DSTRAN(5)+TAO_YZ*DSTRAN(6))

C write(10,*) "DEQPLAS",DEQPLAS,EQPLAS,H

C 更新状态变量

STATEV(1)=EQPLAS+DEQPLAS

C 计算塑性应变增量

DPLAS(1)= DEQPLAS*THREE*CEGMA_SX/TWO/CEGMA_EQ

DPLAS(2)= DEQPLAS*THREE*CEGMA_SY/TWO/CEGMA_EQ

DPLAS(3)= DEQPLAS*THREE*CEGMA_SZ/TWO/CEGMA_EQ

DPLAS(4)= DEQPLAS*THREE*TAO_XY/CEGMA_EQ

DPLAS(5)= DEQPLAS*THREE*TAO_ZX/CEGMA_EQ

DPLAS(6)= DEQPLAS*THREE*TAO_ZY/CEGMA_EQ

C 按照塑性理论更新应力

DO K1=1,NTENS

DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*(DSTRAN(K1)-DPLAS(K1))

END DO

END DO

END IF

CCCCCCCCCCCCCCCCCCCCCCCCCC计算完成

RETURN

END

5.4.​ 切线刚度法程序设计

算法设计

1:定义程序需要用到的常数和变量

2:读取ABAQUS定义的材料常数和状态变量(这里只定义了一个状态变量),材料常数为,弹性模量E,泊松比Mu,屈服应力Yield0,参数A,B,C,并且计算出剪切模量G,状态变量为等效塑性应变EQPLAS

3:读取应力分量,计算平均应力,应力偏量以及Mises等效应力

平均应力:ABAQUS子程序UMAT里弹塑本构的实现的图129

应力偏量:ABAQUS子程序UMAT里弹塑本构的实现的图130

Mises等效应力:ABAQUS子程序UMAT里弹塑本构的实现的图131

4:根据3计算的Mises等效应力和2读取的屈服应力Yield0比较,如果Mises等效应力小于屈服应力,表明此时材料未屈服,那么转到5,否则转到6

5:雅可比矩阵,初始化为0,然后计算弹性矩阵,按照弹性理论更新应力

ABAQUS子程序UMAT里弹塑本构的实现的图132

ABAQUS子程序UMAT里弹塑本构的实现的图133

6:雅可比矩阵,初始化为0

1计算切线模量H

ABAQUS子程序UMAT里弹塑本构的实现的图134

注意到当等效塑性应ABAQUS子程序UMAT里弹塑本构的实现的图135时对应于本构关系的屈服点,此时的H不能通过上式计算,可以取此时的H为弹性模量

2计算w

根据前一章推导的公式

ABAQUS子程序UMAT里弹塑本构的实现的图136

3计算等效塑性应变增量DEQPLAS,并更新状态变量

根据前一章推导的公式:

ABAQUS子程序UMAT里弹塑本构的实现的图137

ABAQUS子程序UMAT里弹塑本构的实现的图138

ABAQUS子程序UMAT里弹塑本构的实现的图139

带入后可得等效塑性应变增量DEQPLAS

ABAQUS子程序UMAT里弹塑本构的实现的图140

然后EQPLAS=DEQPLAS+EQPLAS更新状态变量

4计算雅可比矩阵

首先初始化默认为0,然后用下式计算雅可比矩阵

ABAQUS子程序UMAT里弹塑本构的实现的图141

注意:abaqus的剪应力方向跟弹性力学规定的方向不一致,所以上式的最后两行的ABAQUS子程序UMAT里弹塑本构的实现的图142应交换

5更新应力

ABAQUS子程序UMAT里弹塑本构的实现的图143

算法流程图

ABAQUS子程序UMAT里弹塑本构的实现的图144

图5-5 切线刚度法算法流程图

5.5.​ 切线刚度法程序编码

根据算法流程,用FORTRAN77 固定格式编制了切线刚度法的计算程序如下:

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,

1 DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,

2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,

3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

INCLUDE 'ABA_PARAM.INC'

CHARACTER*8 CMNAME

DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),

1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),

2 PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),

3 DFGRD0(3,3),DFGRD1(3,3)

CCCCCCCCCCCCCCCCCCCCCCCCCC1:定义变量

c 定义常数

PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)

C 定义材料常数

DOUBLE PRECISION E,Mu,Yield0,A,B,C

C 定义中间变量

DOUBLE PRECISION H,W,CEGMA_EQ,CEGMA_X,CEGMA_Y,CEGMA_Z,TAO_XY,

1TAO_YZ,TAO_ZX,CEGMA_CP,CEGMA_SX,CEGMA_SY,CEGMA_SZ,TAO_SXY,TAO_SYZ

2,TAO_SZX

C 定义状态变量

DOUBLE PRECISION EQPLAS

C 定义更新变量

DOUBLE PRECISION DEQPLAS

CCCCCCCCCCCCCCCCCCCCCCCCCC2:读取变量

C 弹性模量

E=PROPS(1)

C 泊松比

Mu=PROPS(2)

C 剪切模量

G=E/TWO/(ONE+Mu)

C 屈服应力

Yield0=PROPS(3)

C 参数

A=PROPS(4)

B=PROPS(5)

C=PROPS(6)

C 读取状态参数

EQPLAS=STATEV(1)

Yield=A*(EQPLAS+Yield0/E)**B+C

CCCCCCCCCCCCCCCCCCCCCCCCCC3:读取应力分量,计算平均应力,应力偏量以及Mises等效应力

C 6个应力分量

CEGMA_X= STRESS(1)

CEGMA_Y= STRESS(2)

CEGMA_Z= STRESS(3)

TAO_XY= STRESS(4)

TAO_ZX= STRESS(5)

TAO_YZ= STRESS(6)

C 平均应力

CEGMA_CP=(CEGMA_X+CEGMA_Y+CEGMA_Z)/THREE

C 6个应力偏量

CEGMA_SX=CEGMA_X-CEGMA_CP

CEGMA_SY=CEGMA_Y-CEGMA_CP

CEGMA_SZ=CEGMA_Z-CEGMA_CP

TAO_SXY=TAO_XY

TAO_SZX=TAO_ZX

TAO_SYZ=TAO_YZ

c Mises等效应力

CEGMA_EQ=SQRT(THREE/TWO*(CEGMA_SX**2+CEGMA_SY**2+CEGMA_SZ**2+

1TWO*(TAO_SXY**2+TAO_SYZ**2+TAO_SZX**2)))

CCCCCCCCCCCCCCCCCCCCCCCCCC4:根据3计算的Mises等效应力和2读取的屈服应力Yield0比较,

CCCCCCCCCCCCCCCCCCCCCCCCCC如果Mises等效应力小于屈服应力,表明此时材料未屈服,那么转到5,否则转到6

IF(CEGMA_EQ.LT.Yield0) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCC5:没有屈服,按照弹性理论计算

C 雅可比矩阵,初始化默认为0

DO k1=1,NTENS

DO k2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO

c 计算弹性矩阵

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K1,k2)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)*(Mu/(ONE-Mu))

END DO

DDSDDE(K1,K1)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)*(ONE-TWO*Mu)/TWO/

1(ONE-Mu)

END DO

C 按照弹性理论更新应力

DO K1=1,NTENS

DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

END DO

END DO

ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCC6:屈服发生,按照塑性理论计算

C 切线模量H,根据本构关系求导

IF(EQPLAS.EQ.0) THEN

H=E

ELSE

H=A*B*EQPLAS**(B-ONE)

END IF

C 求W

W=9.D0*G/TWO/CEGMA_EQ**2/(H+THREE*G)

C 等效塑性应变增量

DEQPLAS=(6.D0*G*CEGMA_EQ/(TWO*CEGMA_EQ**2*H+9.D0*G*(CEGMA_SX

1**2+CEGMA_SY**2+CEGMA_SZ**2+TWO*TAO_XY**2+TWO*TAO_ZX**2+

2TWO*TAO_YZ**2)))*(CEGMA_SX*DSTRAN(1)+CEGMA_SY*DSTRAN(2)+CEGMA_SZ*

3DSTRAN(3)+TAO_XY*DSTRAN(4)+TAO_ZX*DSTRAN(5)+TAO_YZ*DSTRAN(6))

C 更新状态变量

STATEV(1)=EQPLAS+DEQPLAS

C 雅可比矩阵,初始化默认为0

DO k1=1,NTENS

DO k2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO

C 雅可比矩阵更新

DDSDDE(1,1)=E/(ONE+Mu)*((ONE-Mu)/(ONE-TWO*Mu)-W*CEGMA_SX**2)

DDSDDE(2,1)=E/(ONE+Mu)*((Mu)/(ONE-TWO*Mu)-W*CEGMA_SX*CEGMA_Y)

DDSDDE(3,1)=E/(ONE+Mu)*((Mu)/(ONE-TWO*Mu)-W*CEGMA_SX*CEGMA_Z)

DDSDDE(4,1)=E/(ONE+Mu)*(-W*CEGMA_SZ*TAO_SXY)

DDSDDE(5,1)=E/(ONE+Mu)*(-W*CEGMA_SZ*TAO_SZX)

DDSDDE(6,1)=E/(ONE+Mu)*(-W*CEGMA_SZ*TAO_SYZ)

DDSDDE(2,2)=E/(ONE+Mu)*((ONE-Mu)/(ONE-TWO*Mu)-W*CEGMA_Y**2)

DDSDDE(3,2)=E/(ONE+Mu)*((Mu)/(ONE-TWO*Mu)-W*CEGMA_SY*CEGMA_SZ)

DDSDDE(4,2)=E/(ONE+Mu)*(-W*CEGMA_SY*TAO_SXY)

DDSDDE(5,2)=E/(ONE+Mu)*(-W*CEGMA_SY*TAO_SZX)

DDSDDE(6,2)=E/(ONE+Mu)*(-W*CEGMA_SY*TAO_SYZ)

DDSDDE(3,3)=E/(ONE+Mu)*((ONE-Mu)/(ONE-TWO*Mu)-W*CEGMA_SZ**2)

DDSDDE(4,3)=E/(ONE+Mu)*(-W*CEGMA_SZ*TAO_SXY)

DDSDDE(5,3)=E/(ONE+Mu)*(-W*CEGMA_SZ*TAO_SZX)

DDSDDE(6,3)=E/(ONE+Mu)*(-W*CEGMA_SZ*TAO_SYZ)

DDSDDE(4,4)=E/(ONE+Mu)*(ONE/TWO-W*TAO_SXY**2)

DDSDDE(5,4)=E/(ONE+Mu)*(-W*TAO_SXY*TAO_SZX)

DDSDDE(6,4)=E/(ONE+Mu)*(-W*TAO_SXY*TAO_SYZ)

DDSDDE(5,5)=E/(ONE+Mu)*(ONE/TWO-W*TAO_SZX**2)

DDSDDE(6,5)=E/(ONE+Mu)*(-W*TAO_SYZ*TAO_SZX)

DDSDDE(6,6)=E/(ONE+Mu)*(ONE/TWO-W*TAO_SYZ**2)

DDSDDE(1,2)=DDSDDE(2,1)

DDSDDE(1,3)=DDSDDE(3,1)

DDSDDE(1,4)=DDSDDE(4,1)

DDSDDE(1,5)=DDSDDE(5,1)

DDSDDE(1,6)=DDSDDE(6,1)

DDSDDE(2,3)=DDSDDE(3,2)

DDSDDE(2,4)=DDSDDE(4,2)

DDSDDE(2,5)=DDSDDE(5,2)

DDSDDE(2,6)=DDSDDE(6,2)

DDSDDE(3,4)=DDSDDE(4,3)

DDSDDE(3,5)=DDSDDE(5,3)

DDSDDE(3,6)=DDSDDE(6,3)

DDSDDE(4,5)=DDSDDE(5,4)

DDSDDE(4,6)=DDSDDE(6,4)

DDSDDE(5,6)=DDSDDE(6,5)

C 按照塑性理论更新应力

DO K1=1,NTENS

DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

END DO

END DO

END IF

CCCCCCCCCCCCCCCCCCCCCCCCCC计算完成

RETURN

END

5.6.​ 程序的调试

程序设计是否合理,程序编码是否正确并不容易从程序的代码中简单的看出,而是通过在调试中发现的,可以说,程序的调试同程序编码同等重要,通常我们在做UMAT程序的时候,是写好程序直接到ABAQUS提交,这样的方法本身并没有什么错误,但是从程序调试的角度看,这是一种不太好的方法,因为这样做出错的可能性极大,不管是语法错误还是算法错误,一旦出错,我们很难发现到底错误在哪里,到底为什么会出错,所以在这种情况下,我们迫切需要采用一种好的调试手段,以能快速发现程序中的BUG,修正设计,下面介绍一下本文采用的调试方法:

调试方法主要分为2个方面

(1程序语法检查

(2程序算法检查

对于语法检查主要是通过编译器检查错误的方式进行的,但是事实上,如果直接将UMAT子程序拿来给FORTRAN编译器编译,不管是否有语法错误,都是不能编译通过的,这是因为每个UMAT程序都有一句INCLUDE 'ABA_PARAM.INC',这句的意思是需要包含一个名为ABA_PARAM.的头文件,但是事实上,翻遍ABAQUS的目录,我们都不可能找到这个文件,这是因为ABA_PARAM这个头文件的真正文件名并不是ABA_PARAM.INC,而是Site目录下的ABA_PARAM.INC_DP.INC和ABA_PARAM.INC_SP.INC的其中的一个,将2个文件任意一个改名为ABA_PARAM.INC,然后复制到UMAT子程序的同一个目录下,就可以编译通过了,通过这种方法,凡是语法错误,我们可以在编译阶段发现和修改,避免在ABAQUS里反复提交的繁琐手段。

对于程序没有语法错误却有算法错误,这种情况也是很常见的,对程序的算法做检查,本文采用的主要方式是通过输出变量到文本的方法,具体就是在程序内用OPEN语句打开一个文本文件,然后将程序中关心的中间变量用WRITE语句输到该文本,提交给ABAQUS以后,ABAQUS执行到该处就会将其输出到文本,这样,只需要查看最后的文本文件,就可以发现算法中是否纯在错误,修正算法。

通过以上的方法,程序错误就可以在调试中快速解决。

6.​  程序验证

任何有关材料模型的开发都需要经过实践的检验和验证,尤其对于工程问题,通过验证可以发现程序中的BUG,并做相应的修改,本文的验证主要采用综合比较法,整个验证过程如下:

(1)确定一个具体的模型和它的本构关系。

(2)由于ABAQUS自带的弹塑性材料本构关系是经过大量工程实践检验并被认为是一种正确的材料模型,所以对这个模型采用ABAQUS自带的弹塑性材料本构关系来计算,得出其Mises等效应力值的范围以及应力应变关系,该结果作为一个真实正确的解。

(3)采用常刚度算法的UMAT子程序计算,得出其Mises等效应力值的范围以及应力应变关系,然后同ABAQUS计算的相应结果做比较,看两者是否一致,验证常刚度法的正确性。

(4)采用切线刚度算法的UMAT子程序计算,得出其Mises等效应力值的范围以及应力应变关系;然后同ABAQUS计算的相应结果做比较,验证切线刚度法的正确性。

(5)提高非线性的程度,分别用常刚度法和切线刚度法计算,比较两者的效率,得出相应的结论。

6.1.​ 问题描述

用ABAQUS建立一个左端固结的悬臂梁尺寸为20x20x80,简单拉伸,右端作用一个均布力420MPa,材料采用前面所述的内容,弹性模量200000MPa,泊松比0.3,如图6-1所示

ABAQUS子程序UMAT里弹塑本构的实现的图145

图 6-1 验证模型图

6.2.​ 本构关系

采用Mises材料模型,屈服应力400MPa,塑性阶段的应力应变关系为前文所述:ABAQUS子程序UMAT里弹塑本构的实现的图146,其中

A=700 B=0.5 C=400

6.3.​ ABAQUS自带材料模型计算

塑性阶段应力应变关系用WPS表格软件描点得到,然后采用表格形式输入到ABAQUS中,具体取值如表6-1所示(由于篇幅有限,这里没有提供所有数据,详细请见附录的INP文件)

表6-1应力应变关系表

塑性应变

等效应力

塑性应变

等效应力

400

0

419.17028951

0.00075

404.94974747

5E-05

419.79898987

0.0008

407

0.0001

420.40833163

0.00085

408.5732141

0.00015

421

0.0009

409.89949494

0.0002

421.57544901

0.00095

411.06797181

0.00025

422.13594362

0.001

412.12435565

0.0003

422.68259244

0.00105

413.09580085

0.00035

423.21637353

0.0011

414

0.0004

423.73815494

0.00115

414.8492424

0.00045

424.24871131

0.0012

415.65247584

0.0005

424.74873734

0.00125

416.41645516

0.00055

425.23885893

0.0013

417.1464282

0.0006

425.7196423

0.00135

417.8465683

0.00065

426.19160171

0.0014

418.52025918

0.0007

426.65520587

0.00145

采用ABAQUS自带的弹塑性材料应力计算结果如图6-2所示

ABAQUS子程序UMAT里弹塑本构的实现的图147

图6-2 自带弹塑性材料Mises等效应力

图6-2表明Mises等效应力值介于308.3MPa—421.3MPa之间。

应力应变关系如图6-3所示:

ABAQUS子程序UMAT里弹塑本构的实现的图148

图6-3 自带弹塑性材料应力应变关系

图6-3表明,材料的屈服点为400Mpa,塑性阶段按幂函数增长,这个结果与参考材料的应力应变关系完全一致。

6.4.​ 常刚度法的UMAT验证

模型与前文的6.2一致,材料定义如下:

状态变量数:1

用户自定义材料见表6-2:

表6-2 用户自定义材料表

弹性模量E

200000

泊松比Mu

0.3

屈服应力Yield0

400

A

700

B

0.5

C

400

采用常刚度法应力计算结果如图6-4所示

ABAQUS子程序UMAT里弹塑本构的实现的图149

图6-4 常刚度法Mises等效应力

图6-4表明:常刚度法计算的Mises等效应力介于312.1MPa—421.1MPa之间

将常刚度法与ABAQUS计算的比较结果列于表6-3所示:

表6-3

ABAQUS

常刚度法

相差

Mises应力最小值

308.3MPa

312.1MPa

3.8MPa

Mises应力最大值

421.3MPa

422.1MPa

0.8MPa

计算结果比较:等效应力最大值较ABAQUS自带的材料略大0.8Mpa为422.1Mpa,最小值略大3.8Mpa为312.1Mpa,考虑到ABAQUS得计算时采用的表格形式,而本文采用的是直接由本构关系计算的应力应变,两者之间肯定存在误差,同时又考虑到在应力强度达到近420Mpa的情况下,4Mpa以内的波动属于正常范围,故本文的计算结果和ABAQUS自带的材料结果是完全一致的。

应力应变关系:

ABAQUS子程序UMAT里弹塑本构的实现的图150

图6-5 常刚度法应力应变关系

图6-5表明:常刚度法的应力应变关系同前面的ABAQUS计算的应力应变关系以及本文采用材料的应力应变关系均一致。

验证结论:常刚度法UMAT子程序的算法理论和程序设计完全正确。

6.5.​ 切线刚度法的UMAT验证

模型同前文6.2,材料定义如下:

状态变量数:1

用户自定义材料见表6-4:

表6-4 用户自定义材料表

弹性模量E

200000

泊松比Mu

0.3

屈服应力Yield0

400

A

700

B

0.5

C

400

采用切线刚度法应力计算结果如图6-6所示

ABAQUS子程序UMAT里弹塑本构的实现的图151

图6-6,切线刚度法Mises等效应力

图6-6表明:切线刚度法计算的Mises等效应力介于307.9MPa—422.1MPa之间

ABAQUS比较结果如表6-5所示:

表6-5

ABAQUS

切线刚度法

相差

Mises应力最小值

308.3MPa

307.9MPa

0.4MPa

Mises应力最大值

421.3MPa

422.1MPa

0.8MPa

计算结果比较:切线刚度法计算的等效应力最大值较ABAQUS自带的材料略大0.8Mpa为422.1Mpa,最小值略小0.4Mpa为307.9Mpa,两者相差极小,考虑到误差问题,我们可以看出,切线刚度法的计算结果几乎和ABAQUS自带的材料模型计算结果完全一致。

应力应变关系:

ABAQUS子程序UMAT里弹塑本构的实现的图152

图6-7,切线刚度法应力应变关系

图6-7表明:切线刚度法的盈利应变关系与前面的ABAQUS计算的应力应变关系以及本文采用材料的应力应变关系均一致。

验证结论:切线刚度法UMAT子程序的算法理论和程序设计完全正确。

6.6.​ 两种算法的比较分析

接下来,将改变本构关系中的常数,对比常刚度法和切线刚度法在非线性程度较高的情况下的算法效率,将模型的A值修改为2000,荷载修改为-470Mpa,重新运行分析

ABAQUS子程序UMAT里弹塑本构的实现的图153

图6-8 常刚度法Mises等效应力

ABAQUS子程序UMAT里弹塑本构的实现的图154

图6-9 常刚度法应力应变关系

ABAQUS子程序UMAT里弹塑本构的实现的图155

图6-9 切线刚度法Mises等效应力

ABAQUS子程序UMAT里弹塑本构的实现的图156

图6-11 切线刚度法应力应变关系

可以看出两种算法得到的Mises计算结果完全正确,它们的应力应变关系与材料的应许力应变关系也是完全一致得。

在算法的效率方面,比较结果见表6-7

表6-7 算法时间比较

常刚度法

21.3秒

切线刚度法

18.1秒

注意:我们考虑到整个计算时间包含了弹性部分,所以严格的说,切线刚度法在塑性阶段的优势更加明显,通过我的做的大量计算分析来看,常刚度法在塑性阶段的迭代时间大约是切线刚度法的2-3倍左右,也就是说,切线刚度法在计算非线性程度较高的材料非线性问题时具有明显的优势。

结论如下:

1:常刚度UMAT子程序的算法理论和公式推导均正确,程序编码无误

2:切线刚度法UMAT子程序的算法理论和公式推导均正确,程序编码无误

3:切线刚度法在处理非线性程度较高的材料模型时相对于常刚度法,能体现出明显的优势,这种优势包括,计算时间更少,收敛效率更高,通常在计算该类问题时应该采用切线刚度法。

4:本验算模型实际上计算规模上是比较小的,如果非线性程度更高,模型更复杂,那么采用切线刚度法将大大节省计算和分析的时间和成本。

5:本文采用的分析和研究手段可以推广到其他UMAT材料模型的开发中。

7.​  结论与展望

7.1.​ 结论

本文首先探讨了ABAQUS软件的二次开发接口,对其中的UMAT子程序的原理做了详细的阐述,然后从材料的弹塑性本构关系入手,研究了ABAQUS用户子程序UMAT的具体开发方法,并对常刚度法和切线刚度法的算法理论和程序设计做了深入的剖析,具体有:

(1)深入探讨了ABAQUS软件极其二次开发平台。

(2)对用户子程序UMAT的编译环境设置,接口原理,使用规范,以及开发方法做了详尽的阐述,并总结了常用的开发过程。

(3)对材料非线性问题展开了的探讨,并着重对其中的迭代算法部分进行了深入的分析和比较,对切线刚度法和其他算法做了比较,指出了其优势所在,最后针对各向同性硬化材料的常刚度法和切线刚度法的计算公式做了推导,详尽的阐述了其计算原理。

(4)针对弹塑性材料本构关系,对其进行了抽象,设计了计算该类材料模型的常刚度法和切线刚度法的程序。

(5)为了验证程序,建立了一个单轴拉伸模型,分别用ABAQUS自带弹塑性材料本构关系和本文开发的UMAT程序进行了对比和分析,验证了程序的正确性,在研究过程中主要运用了综合比较分析法。

7.2.​ 展望

通过本文的研究工作已经完成,但是我认为可以从以下2个方面继续对用户材料子程序UMAT展开研究:

(1)本文所编制的UMAT子程序中,采用的是Mises屈服准则,该准则主要可以应用于各向同性硬化材料,比如延性金属材料(钢筋),事实上Mises模型自身是存在一定局限性的,如果本模型拿来做混凝土材料的分析会不够准确,所以如果是做混凝土材料的本构模型开发的话,需要用到更符合实际情况的Drucker-Prager模型,在以后的开发中,可以通过UMAT用户材料子程序,实现该模型,来做混凝土材料的研究。

(2本文没有考虑时间因素,时间其实同样是UMAT子程序里的一个重要的变量,我们完全可以用考虑时间的方式将本构关系的表达式表示成率的形式,然后重新编写本文的代码,结果肯定是一样的,事实上,如果考虑时间效应来开发材料的本构模型,我们就可以对工程问题进行粘弹性分析或粘塑性分析,比如研究混凝土材料的收缩徐变,动态损伤与断裂,岩体的长期稳定性等复杂问题。

致 谢

经过两个月的努力,我完成了基于ABAQUS用户子程序UMAT的二次开发工作,在此期间得到老师和同学的关心和帮助,在这表示由衷的感谢。

感谢指导教师周东教授,在编程过程中,他给了我大力的支持和悉心指导。周老师严谨的治学作风,渊博的知识无时无刻不在影响着我。

感谢在大学4年中给过我无私关怀的王家林教授,陈世民副教授,以及力学系的其他老师们,你们实事求是的工作作风和崇高的敬业精神,感染着我,激励着我不断前进。

再次,感谢所有曾经关心和帮助过我的人,感谢家人的理解和支持,使我顺利完成了这次的毕业论文。

最后,谨向百忙之中审阅我程序设计和毕业论文的老师们以及参加答辩的每一位老师表示最诚挚的谢意!

参考文献

[1]王勖成、邵敏,有限单元法基本原理和数值方法(第2版)[M],清华大学出版社,1997

[2]ABAQUS/Standard Theory Manual,Hibbitt,Karlsson,Sorensen,Inc.2002

[3]Abaqus/CAE User's Manual[C]

[4]石亦平、周玉蓉,ABAQUS有限元分析实例详解[M],机械工业出版社,2006

[5]Abaqus User Subroutines Reference Manual[C]

[6]谭世语、FORTRAN程序设计[M], 重庆大学出版社,2002

[7]朱以文,蔡元奇等译,ABAQUS/Standard 有限元软件入门指南

[8]庄茁,张帆,岑松,ABAQUS非线性有限元分析与实例[M],科学出版社,2005

[9]D.R.J.欧文,塑性力学有限元理论与应用[M],兵器工业出版社,1989

[10]薛守义,弹塑性力学[M],中国建筑工业出版社,2005.7

[11]陈惠发、A.F.萨里普著,余天庆、王勋文、刘再华译,弹性与塑性力学[M],中国建筑工业出版社,2004.6

[12]殷有泉,非线性有限元基础[M],北京大学出版社,2007.12

[13]朱伯芳,有限单元法原理与应用(第2版)[M],中国水利水电出版社,1998

[14]吕西林,钢筋混凝土结构非线性有限元理论与应用[M],同济大学出版社,1997

[15]张汝清,詹先义,非线性有限元分析[M],重庆大学出版社,1990

附1:ABAQUS自带弹塑性材料验证的INP文件

*Heading

** Job name: beam6 Model name: Model-1

** Generated by: Abaqus/CAE Version 6.8-1

*Preprint, echo=NO, model=NO, history=NO, contact=NO

**

** PARTS

**

*Part, name=Part-1

*Node

1, 0., 20., 80.

2, 0., 13.333333, 80.

3, 0., 6.66666651, 80.

4, 0., 0., 80.

5, 0., 20., 72.

6, 0., 13.333333, 72.

7, 0., 6.66666651, 72.

8, 0., 0., 72.

9, 0., 20., 64.

10, 0., 13.333333, 64.

11, 0., 6.66666651, 64.

12, 0., 0., 64.

13, 0., 20., 56.

14, 0., 13.333333, 56.

15, 0., 6.66666651, 56.

16, 0., 0., 56.

17, 0., 20., 48.

18, 0., 13.333333, 48.

19, 0., 6.66666651, 48.

20, 0., 0., 48.

21, 0., 20., 40.

22, 0., 13.333333, 40.

23, 0., 6.66666651, 40.

24, 0., 0., 40.

25, 0., 20., 32.

26, 0., 13.333333, 32.

27, 0., 6.66666651, 32.

28, 0., 0., 32.

29, 0., 20., 24.

30, 0., 13.333333, 24.

31, 0., 6.66666651, 24.

32, 0., 0., 24.

33, 0., 20., 16.

34, 0., 13.333333, 16.

35, 0., 6.66666651, 16.

36, 0., 0., 16.

37, 0., 20., 8.

38, 0., 13.333333, 8.

39, 0., 6.66666651, 8.

40, 0., 0., 8.

41, 0., 20., 0.

42, 0., 13.333333, 0.

43, 0., 6.66666651, 0.

44, 0., 0., 0.

45, -6.66666651, 20., 80.

46, -6.66666651, 13.333333, 80.

47, -6.66666651, 6.66666651, 80.

48, -6.66666651, 0., 80.

49, -6.66666651, 20., 72.

50, -6.66666651, 13.333333, 72.

51, -6.66666651, 6.66666651, 72.

52, -6.66666651, 0., 72.

53, -6.66666651, 20., 64.

54, -6.66666651, 13.333333, 64.

55, -6.66666651, 6.66666651, 64.

56, -6.66666651, 0., 64.

57, -6.66666651, 20., 56.

58, -6.66666651, 13.333333, 56.

59, -6.66666651, 6.66666651, 56.

60, -6.66666651, 0., 56.

61, -6.66666651, 20., 48.

62, -6.66666651, 13.333333, 48.

63, -6.66666651, 6.66666651, 48.

64, -6.66666651, 0., 48.

65, -6.66666651, 20., 40.

66, -6.66666651, 13.333333, 40.

67, -6.66666651, 6.66666651, 40.

68, -6.66666651, 0., 40.

69, -6.66666651, 20., 32.

70, -6.66666651, 13.333333, 32.

71, -6.66666651, 6.66666651, 32.

72, -6.66666651, 0., 32.

73, -6.66666651, 20., 24.

74, -6.66666651, 13.333333, 24.

75, -6.66666651, 6.66666651, 24.

76, -6.66666651, 0., 24.

77, -6.66666651, 20., 16.

78, -6.66666651, 13.333333, 16.

79, -6.66666651, 6.66666651, 16.

80, -6.66666651, 0., 16.

81, -6.66666651, 20., 8.

82, -6.66666651, 13.333333, 8.

83, -6.66666651, 6.66666651, 8.

84, -6.66666651, 0., 8.

85, -6.66666651, 20., 0.

86, -6.66666651, 13.333333, 0.

87, -6.66666651, 6.66666651, 0.

88, -6.66666651, 0., 0.

89, -13.333333, 20., 80.

90, -13.333333, 13.333333, 80.

91, -13.333333, 6.66666651, 80.

92, -13.333333, 0., 80.

93, -13.333333, 20., 72.

94, -13.333333, 13.333333, 72.

95, -13.333333, 6.66666651, 72.

96, -13.333333, 0., 72.

97, -13.333333, 20., 64.

98, -13.333333, 13.333333, 64.

99, -13.333333, 6.66666651, 64.

100, -13.333333, 0., 64.

101, -13.333333, 20., 56.

102, -13.333333, 13.333333, 56.

103, -13.333333, 6.66666651, 56.

104, -13.333333, 0., 56.

105, -13.333333, 20., 48.

106, -13.333333, 13.333333, 48.

107, -13.333333, 6.66666651, 48.

108, -13.333333, 0., 48.

109, -13.333333, 20., 40.

110, -13.333333, 13.333333, 40.

111, -13.333333, 6.66666651, 40.

112, -13.333333, 0., 40.

113, -13.333333, 20., 32.

114, -13.333333, 13.333333, 32.

115, -13.333333, 6.66666651, 32.

116, -13.333333, 0., 32.

117, -13.333333, 20., 24.

118, -13.333333, 13.333333, 24.

119, -13.333333, 6.66666651, 24.

120, -13.333333, 0., 24.

121, -13.333333, 20., 16.

122, -13.333333, 13.333333, 16.

123, -13.333333, 6.66666651, 16.

124, -13.333333, 0., 16.

125, -13.333333, 20., 8.

126, -13.333333, 13.333333, 8.

127, -13.333333, 6.66666651, 8.

128, -13.333333, 0., 8.

129, -13.333333, 20., 0.

130, -13.333333, 13.333333, 0.

131, -13.333333, 6.66666651, 0.

132, -13.333333, 0., 0.

133, -20., 20., 80.

134, -20., 13.333333, 80.

135, -20., 6.66666651, 80.

136, -20., 0., 80.

137, -20., 20., 72.

138, -20., 13.333333, 72.

139, -20., 6.66666651, 72.

140, -20., 0., 72.

141, -20., 20., 64.

142, -20., 13.333333, 64.

143, -20., 6.66666651, 64.

144, -20., 0., 64.

145, -20., 20., 56.

146, -20., 13.333333, 56.

147, -20., 6.66666651, 56.

148, -20., 0., 56.

149, -20., 20., 48.

150, -20., 13.333333, 48.

151, -20., 6.66666651, 48.

152, -20., 0., 48.

153, -20., 20., 40.

154, -20., 13.333333, 40.

155, -20., 6.66666651, 40.

156, -20., 0., 40.

157, -20., 20., 32.

158, -20., 13.333333, 32.

159, -20., 6.66666651, 32.

160, -20., 0., 32.

161, -20., 20., 24.

162, -20., 13.333333, 24.

163, -20., 6.66666651, 24.

164, -20., 0., 24.

165, -20., 20., 16.

166, -20., 13.333333, 16.

167, -20., 6.66666651, 16.

168, -20., 0., 16.

169, -20., 20., 8.

170, -20., 13.333333, 8.

171, -20., 6.66666651, 8.

172, -20., 0., 8.

173, -20., 20., 0.

174, -20., 13.333333, 0.

175, -20., 6.66666651, 0.

176, -20., 0., 0.

*Element, type=C3D8R

1, 45, 46, 50, 49, 1, 2, 6, 5

2, 46, 47, 51, 50, 2, 3, 7, 6

3, 47, 48, 52, 51, 3, 4, 8, 7

4, 49, 50, 54, 53, 5, 6, 10, 9

5, 50, 51, 55, 54, 6, 7, 11, 10

6, 51, 52, 56, 55, 7, 8, 12, 11

7, 53, 54, 58, 57, 9, 10, 14, 13

8, 54, 55, 59, 58, 10, 11, 15, 14

9, 55, 56, 60, 59, 11, 12, 16, 15

10, 57, 58, 62, 61, 13, 14, 18, 17

11, 58, 59, 63, 62, 14, 15, 19, 18

12, 59, 60, 64, 63, 15, 16, 20, 19

13, 61, 62, 66, 65, 17, 18, 22, 21

14, 62, 63, 67, 66, 18, 19, 23, 22

15, 63, 64, 68, 67, 19, 20, 24, 23

16, 65, 66, 70, 69, 21, 22, 26, 25

17, 66, 67, 71, 70, 22, 23, 27, 26

18, 67, 68, 72, 71, 23, 24, 28, 27

19, 69, 70, 74, 73, 25, 26, 30, 29

20, 70, 71, 75, 74, 26, 27, 31, 30

21, 71, 72, 76, 75, 27, 28, 32, 31

22, 73, 74, 78, 77, 29, 30, 34, 33

23, 74, 75, 79, 78, 30, 31, 35, 34

24, 75, 76, 80, 79, 31, 32, 36, 35

25, 77, 78, 82, 81, 33, 34, 38, 37

26, 78, 79, 83, 82, 34, 35, 39, 38

27, 79, 80, 84, 83, 35, 36, 40, 39

28, 81, 82, 86, 85, 37, 38, 42, 41

29, 82, 83, 87, 86, 38, 39, 43, 42

30, 83, 84, 88, 87, 39, 40, 44, 43

31, 89, 90, 94, 93, 45, 46, 50, 49

32, 90, 91, 95, 94, 46, 47, 51, 50

33, 91, 92, 96, 95, 47, 48, 52, 51

34, 93, 94, 98, 97, 49, 50, 54, 53

35, 94, 95, 99, 98, 50, 51, 55, 54

36, 95, 96, 100, 99, 51, 52, 56, 55

37, 97, 98, 102, 101, 53, 54, 58, 57

38, 98, 99, 103, 102, 54, 55, 59, 58

39, 99, 100, 104, 103, 55, 56, 60, 59

40, 101, 102, 106, 105, 57, 58, 62, 61

41, 102, 103, 107, 106, 58, 59, 63, 62

42, 103, 104, 108, 107, 59, 60, 64, 63

43, 105, 106, 110, 109, 61, 62, 66, 65

44, 106, 107, 111, 110, 62, 63, 67, 66

45, 107, 108, 112, 111, 63, 64, 68, 67

46, 109, 110, 114, 113, 65, 66, 70, 69

47, 110, 111, 115, 114, 66, 67, 71, 70

48, 111, 112, 116, 115, 67, 68, 72, 71

49, 113, 114, 118, 117, 69, 70, 74, 73

50, 114, 115, 119, 118, 70, 71, 75, 74

51, 115, 116, 120, 119, 71, 72, 76, 75

52, 117, 118, 122, 121, 73, 74, 78, 77

53, 118, 119, 123, 122, 74, 75, 79, 78

54, 119, 120, 124, 123, 75, 76, 80, 79

55, 121, 122, 126, 125, 77, 78, 82, 81

56, 122, 123, 127, 126, 78, 79, 83, 82

57, 123, 124, 128, 127, 79, 80, 84, 83

58, 125, 126, 130, 129, 81, 82, 86, 85

59, 126, 127, 131, 130, 82, 83, 87, 86

60, 127, 128, 132, 131, 83, 84, 88, 87

61, 133, 134, 138, 137, 89, 90, 94, 93

62, 134, 135, 139, 138, 90, 91, 95, 94

63, 135, 136, 140, 139, 91, 92, 96, 95

64, 137, 138, 142, 141, 93, 94, 98, 97

65, 138, 139, 143, 142, 94, 95, 99, 98

66, 139, 140, 144, 143, 95, 96, 100, 99

67, 141, 142, 146, 145, 97, 98, 102, 101

68, 142, 143, 147, 146, 98, 99, 103, 102

69, 143, 144, 148, 147, 99, 100, 104, 103

70, 145, 146, 150, 149, 101, 102, 106, 105

71, 146, 147, 151, 150, 102, 103, 107, 106

72, 147, 148, 152, 151, 103, 104, 108, 107

73, 149, 150, 154, 153, 105, 106, 110, 109

74, 150, 151, 155, 154, 106, 107, 111, 110

75, 151, 152, 156, 155, 107, 108, 112, 111

76, 153, 154, 158, 157, 109, 110, 114, 113

77, 154, 155, 159, 158, 110, 111, 115, 114

78, 155, 156, 160, 159, 111, 112, 116, 115

79, 157, 158, 162, 161, 113, 114, 118, 117

80, 158, 159, 163, 162, 114, 115, 119, 118

81, 159, 160, 164, 163, 115, 116, 120, 119

82, 161, 162, 166, 165, 117, 118, 122, 121

83, 162, 163, 167, 166, 118, 119, 123, 122

84, 163, 164, 168, 167, 119, 120, 124, 123

85, 165, 166, 170, 169, 121, 122, 126, 125

86, 166, 167, 171, 170, 122, 123, 127, 126

87, 167, 168, 172, 171, 123, 124, 128, 127

88, 169, 170, 174, 173, 125, 126, 130, 129

89, 170, 171, 175, 174, 126, 127, 131, 130

90, 171, 172, 176, 175, 127, 128, 132, 131

*Nset, nset=_PickedSet2, internal, generate

1, 176, 1

*Elset, elset=_PickedSet2, internal, generate

1, 90, 1

** Section: Section-1

*Solid Section, elset=_PickedSet2, controls=EC-1, material=Material-1

*End Part

**

**

** ASSEMBLY

**

*Assembly, name=Assembly

**

*Instance, name=Part-1-1, part=Part-1

*End Instance

**

*Nset, nset=_PickedSet5, internal, instance=Part-1-1

1, 2, 3, 4, 45, 46, 47, 48, 89, 90, 91, 92, 133, 134, 135, 136

*Elset, elset=_PickedSet5, internal, instance=Part-1-1

1, 2, 3, 31, 32, 33, 61, 62, 63

*Elset, elset=Set-1, instance=Part-1-1

28,

*Elset, elset=__PickedSurf6_S5, internal, instance=Part-1-1

28, 29, 30, 58, 59, 60, 88, 89, 90

*Surface, type=ELEMENT, name=_PickedSurf6, internal

__PickedSurf6_S5, S5

*End Assembly

**

** ELEMENT CONTROLS

**

*Section Controls, name=EC-1, hourglass=ENHANCED

1., 1., 1.

**

** MATERIALS

**

*Material, name=Material-1

*Elastic

200000., 0.3

*Plastic

400., 0.

404.95, 5e-05

407., 0.0001

408.573, 0.00015

409.899, 0.0002

411.068, 0.00025

412.124, 0.0003

413.096, 0.00035

414., 0.0004

414.849, 0.00045

415.652, 0.0005

416.416, 0.00055

417.146, 0.0006

417.847, 0.00065

418.52, 0.0007

419.17, 0.00075

419.799, 0.0008

420.408, 0.00085

421., 0.0009

421.575, 0.00095

422.136, 0.001

422.683, 0.00105

423.216, 0.0011

423.738, 0.00115

424.249, 0.0012

424.749, 0.00125

425.239, 0.0013

425.72, 0.00135

426.192, 0.0014

426.655, 0.00145

427.111, 0.0015

427.559, 0.00155

428., 0.0016

428.434, 0.00165

428.862, 0.0017

429.283, 0.00175

429.698, 0.0018

430.108, 0.00185

430.512, 0.0019

430.911, 0.00195

431.305, 0.002

431.694, 0.00205

432.078, 0.0021

432.458, 0.00215

432.833, 0.0022

433.204, 0.00225

433.571, 0.0023

433.934, 0.00235

434.293, 0.0024

434.648, 0.00245

435., 0.0025

435.348, 0.00255

435.693, 0.0026

436.035, 0.00265

436.373, 0.0027

436.708, 0.00275

437.041, 0.0028

437.37, 0.00285

437.696, 0.0029

438.02, 0.00295

438.341, 0.003

438.659, 0.00305

438.974, 0.0031

439.287, 0.00315

439.598, 0.0032

439.906, 0.00325

440.212, 0.0033

440.515, 0.00335

440.817, 0.0034

441.116, 0.00345

441.413, 0.0035

441.707, 0.00355

442., 0.0036

442.291, 0.00365

442.579, 0.0037

442.866, 0.00375

443.151, 0.0038

443.434, 0.00385

** ----------------------------------------------------------------

**

** STEP: Step-1

**

*Step, name=Step-1

*Static, direct

0.01, 1.,

**

** BOUNDARY CONDITIONS

**

** Name: BC-1 Type: 位移/转角

*Boundary

_PickedSet5, 1, 1

_PickedSet5, 2, 2

_PickedSet5, 3, 3

_PickedSet5, 4, 4

_PickedSet5, 5, 5

_PickedSet5, 6, 6

**

** LOADS

**

** Name: Load-1 Type: Pressure

*Dsload

_PickedSurf6, P, -420.

**

** OUTPUT REQUESTS

**

*Restart, write, frequency=0

**

** FIELD OUTPUT: F-Output-1

**

*Output, field, variable=PRESELECT

**

** HISTORY OUTPUT: H-Output-1

**

*Output, history

*Element Output, elset=Set-1

E33, S33

*End Step

附2:用于算法验证的INP文件

*Heading

** Job name: beam3 Model name: Model-1

** Generated by: Abaqus/CAE Version 6.8-1

*Preprint, echo=NO, model=NO, history=NO, contact=NO

**

** PARTS

**

*Part, name=Part-1

*Node

1, 0., 20., 80.

2, 0., 13.333333, 80.

3, 0., 6.66666651, 80.

4, 0., 0., 80.

5, 0., 20., 72.

6, 0., 13.333333, 72.

7, 0., 6.66666651, 72.

8, 0., 0., 72.

9, 0., 20., 64.

10, 0., 13.333333, 64.

11, 0., 6.66666651, 64.

12, 0., 0., 64.

13, 0., 20., 56.

14, 0., 13.333333, 56.

15, 0., 6.66666651, 56.

16, 0., 0., 56.

17, 0., 20., 48.

18, 0., 13.333333, 48.

19, 0., 6.66666651, 48.

20, 0., 0., 48.

21, 0., 20., 40.

22, 0., 13.333333, 40.

23, 0., 6.66666651, 40.

24, 0., 0., 40.

25, 0., 20., 32.

26, 0., 13.333333, 32.

27, 0., 6.66666651, 32.

28, 0., 0., 32.

29, 0., 20., 24.

30, 0., 13.333333, 24.

31, 0., 6.66666651, 24.

32, 0., 0., 24.

33, 0., 20., 16.

34, 0., 13.333333, 16.

35, 0., 6.66666651, 16.

36, 0., 0., 16.

37, 0., 20., 8.

38, 0., 13.333333, 8.

39, 0., 6.66666651, 8.

40, 0., 0., 8.

41, 0., 20., 0.

42, 0., 13.333333, 0.

43, 0., 6.66666651, 0.

44, 0., 0., 0.

45, -6.66666651, 20., 80.

46, -6.66666651, 13.333333, 80.

47, -6.66666651, 6.66666651, 80.

48, -6.66666651, 0., 80.

49, -6.66666651, 20., 72.

50, -6.66666651, 13.333333, 72.

51, -6.66666651, 6.66666651, 72.

52, -6.66666651, 0., 72.

53, -6.66666651, 20., 64.

54, -6.66666651, 13.333333, 64.

55, -6.66666651, 6.66666651, 64.

56, -6.66666651, 0., 64.

57, -6.66666651, 20., 56.

58, -6.66666651, 13.333333, 56.

59, -6.66666651, 6.66666651, 56.

60, -6.66666651, 0., 56.

61, -6.66666651, 20., 48.

62, -6.66666651, 13.333333, 48.

63, -6.66666651, 6.66666651, 48.

64, -6.66666651, 0., 48.

65, -6.66666651, 20., 40.

66, -6.66666651, 13.333333, 40.

67, -6.66666651, 6.66666651, 40.

68, -6.66666651, 0., 40.

69, -6.66666651, 20., 32.

70, -6.66666651, 13.333333, 32.

71, -6.66666651, 6.66666651, 32.

72, -6.66666651, 0., 32.

73, -6.66666651, 20., 24.

74, -6.66666651, 13.333333, 24.

75, -6.66666651, 6.66666651, 24.

76, -6.66666651, 0., 24.

77, -6.66666651, 20., 16.

78, -6.66666651, 13.333333, 16.

79, -6.66666651, 6.66666651, 16.

80, -6.66666651, 0., 16.

81, -6.66666651, 20., 8.

82, -6.66666651, 13.333333, 8.

83, -6.66666651, 6.66666651, 8.

84, -6.66666651, 0., 8.

85, -6.66666651, 20., 0.

86, -6.66666651, 13.333333, 0.

87, -6.66666651, 6.66666651, 0.

88, -6.66666651, 0., 0.

89, -13.333333, 20., 80.

90, -13.333333, 13.333333, 80.

91, -13.333333, 6.66666651, 80.

92, -13.333333, 0., 80.

93, -13.333333, 20., 72.

94, -13.333333, 13.333333, 72.

95, -13.333333, 6.66666651, 72.

96, -13.333333, 0., 72.

97, -13.333333, 20., 64.

98, -13.333333, 13.333333, 64.

99, -13.333333, 6.66666651, 64.

100, -13.333333, 0., 64.

101, -13.333333, 20., 56.

102, -13.333333, 13.333333, 56.

103, -13.333333, 6.66666651, 56.

104, -13.333333, 0., 56.

105, -13.333333, 20., 48.

106, -13.333333, 13.333333, 48.

107, -13.333333, 6.66666651, 48.

108, -13.333333, 0., 48.

109, -13.333333, 20., 40.

110, -13.333333, 13.333333, 40.

111, -13.333333, 6.66666651, 40.

112, -13.333333, 0., 40.

113, -13.333333, 20., 32.

114, -13.333333, 13.333333, 32.

115, -13.333333, 6.66666651, 32.

116, -13.333333, 0., 32.

117, -13.333333, 20., 24.

118, -13.333333, 13.333333, 24.

119, -13.333333, 6.66666651, 24.

120, -13.333333, 0., 24.

121, -13.333333, 20., 16.

122, -13.333333, 13.333333, 16.

123, -13.333333, 6.66666651, 16.

124, -13.333333, 0., 16.

125, -13.333333, 20., 8.

126, -13.333333, 13.333333, 8.

127, -13.333333, 6.66666651, 8.

128, -13.333333, 0., 8.

129, -13.333333, 20., 0.

130, -13.333333, 13.333333, 0.

131, -13.333333, 6.66666651, 0.

132, -13.333333, 0., 0.

133, -20., 20., 80.

134, -20., 13.333333, 80.

135, -20., 6.66666651, 80.

136, -20., 0., 80.

137, -20., 20., 72.

138, -20., 13.333333, 72.

139, -20., 6.66666651, 72.

140, -20., 0., 72.

141, -20., 20., 64.

142, -20., 13.333333, 64.

143, -20., 6.66666651, 64.

144, -20., 0., 64.

145, -20., 20., 56.

146, -20., 13.333333, 56.

147, -20., 6.66666651, 56.

148, -20., 0., 56.

149, -20., 20., 48.

150, -20., 13.333333, 48.

151, -20., 6.66666651, 48.

152, -20., 0., 48.

153, -20., 20., 40.

154, -20., 13.333333, 40.

155, -20., 6.66666651, 40.

156, -20., 0., 40.

157, -20., 20., 32.

158, -20., 13.333333, 32.

159, -20., 6.66666651, 32.

160, -20., 0., 32.

161, -20., 20., 24.

162, -20., 13.333333, 24.

163, -20., 6.66666651, 24.

164, -20., 0., 24.

165, -20., 20., 16.

166, -20., 13.333333, 16.

167, -20., 6.66666651, 16.

168, -20., 0., 16.

169, -20., 20., 8.

170, -20., 13.333333, 8.

171, -20., 6.66666651, 8.

172, -20., 0., 8.

173, -20., 20., 0.

174, -20., 13.333333, 0.

175, -20., 6.66666651, 0.

176, -20., 0., 0.

*Element, type=C3D8R

1, 45, 46, 50, 49, 1, 2, 6, 5

2, 46, 47, 51, 50, 2, 3, 7, 6

3, 47, 48, 52, 51, 3, 4, 8, 7

4, 49, 50, 54, 53, 5, 6, 10, 9

5, 50, 51, 55, 54, 6, 7, 11, 10

6, 51, 52, 56, 55, 7, 8, 12, 11

7, 53, 54, 58, 57, 9, 10, 14, 13

8, 54, 55, 59, 58, 10, 11, 15, 14

9, 55, 56, 60, 59, 11, 12, 16, 15

10, 57, 58, 62, 61, 13, 14, 18, 17

11, 58, 59, 63, 62, 14, 15, 19, 18

12, 59, 60, 64, 63, 15, 16, 20, 19

13, 61, 62, 66, 65, 17, 18, 22, 21

14, 62, 63, 67, 66, 18, 19, 23, 22

15, 63, 64, 68, 67, 19, 20, 24, 23

16, 65, 66, 70, 69, 21, 22, 26, 25

17, 66, 67, 71, 70, 22, 23, 27, 26

18, 67, 68, 72, 71, 23, 24, 28, 27

19, 69, 70, 74, 73, 25, 26, 30, 29

20, 70, 71, 75, 74, 26, 27, 31, 30

21, 71, 72, 76, 75, 27, 28, 32, 31

22, 73, 74, 78, 77, 29, 30, 34, 33

23, 74, 75, 79, 78, 30, 31, 35, 34

24, 75, 76, 80, 79, 31, 32, 36, 35

25, 77, 78, 82, 81, 33, 34, 38, 37

26, 78, 79, 83, 82, 34, 35, 39, 38

27, 79, 80, 84, 83, 35, 36, 40, 39

28, 81, 82, 86, 85, 37, 38, 42, 41

29, 82, 83, 87, 86, 38, 39, 43, 42

30, 83, 84, 88, 87, 39, 40, 44, 43

31, 89, 90, 94, 93, 45, 46, 50, 49

32, 90, 91, 95, 94, 46, 47, 51, 50

33, 91, 92, 96, 95, 47, 48, 52, 51

34, 93, 94, 98, 97, 49, 50, 54, 53

35, 94, 95, 99, 98, 50, 51, 55, 54

36, 95, 96, 100, 99, 51, 52, 56, 55

37, 97, 98, 102, 101, 53, 54, 58, 57

38, 98, 99, 103, 102, 54, 55, 59, 58

39, 99, 100, 104, 103, 55, 56, 60, 59

40, 101, 102, 106, 105, 57, 58, 62, 61

41, 102, 103, 107, 106, 58, 59, 63, 62

42, 103, 104, 108, 107, 59, 60, 64, 63

43, 105, 106, 110, 109, 61, 62, 66, 65

44, 106, 107, 111, 110, 62, 63, 67, 66

45, 107, 108, 112, 111, 63, 64, 68, 67

46, 109, 110, 114, 113, 65, 66, 70, 69

47, 110, 111, 115, 114, 66, 67, 71, 70

48, 111, 112, 116, 115, 67, 68, 72, 71

49, 113, 114, 118, 117, 69, 70, 74, 73

50, 114, 115, 119, 118, 70, 71, 75, 74

51, 115, 116, 120, 119, 71, 72, 76, 75

52, 117, 118, 122, 121, 73, 74, 78, 77

53, 118, 119, 123, 122, 74, 75, 79, 78

54, 119, 120, 124, 123, 75, 76, 80, 79

55, 121, 122, 126, 125, 77, 78, 82, 81

56, 122, 123, 127, 126, 78, 79, 83, 82

57, 123, 124, 128, 127, 79, 80, 84, 83

58, 125, 126, 130, 129, 81, 82, 86, 85

59, 126, 127, 131, 130, 82, 83, 87, 86

60, 127, 128, 132, 131, 83, 84, 88, 87

61, 133, 134, 138, 137, 89, 90, 94, 93

62, 134, 135, 139, 138, 90, 91, 95, 94

63, 135, 136, 140, 139, 91, 92, 96, 95

64, 137, 138, 142, 141, 93, 94, 98, 97

65, 138, 139, 143, 142, 94, 95, 99, 98

66, 139, 140, 144, 143, 95, 96, 100, 99

67, 141, 142, 146, 145, 97, 98, 102, 101

68, 142, 143, 147, 146, 98, 99, 103, 102

69, 143, 144, 148, 147, 99, 100, 104, 103

70, 145, 146, 150, 149, 101, 102, 106, 105

71, 146, 147, 151, 150, 102, 103, 107, 106

72, 147, 148, 152, 151, 103, 104, 108, 107

73, 149, 150, 154, 153, 105, 106, 110, 109

74, 150, 151, 155, 154, 106, 107, 111, 110

75, 151, 152, 156, 155, 107, 108, 112, 111

76, 153, 154, 158, 157, 109, 110, 114, 113

77, 154, 155, 159, 158, 110, 111, 115, 114

78, 155, 156, 160, 159, 111, 112, 116, 115

79, 157, 158, 162, 161, 113, 114, 118, 117

80, 158, 159, 163, 162, 114, 115, 119, 118

81, 159, 160, 164, 163, 115, 116, 120, 119

82, 161, 162, 166, 165, 117, 118, 122, 121

83, 162, 163, 167, 166, 118, 119, 123, 122

84, 163, 164, 168, 167, 119, 120, 124, 123

85, 165, 166, 170, 169, 121, 122, 126, 125

86, 166, 167, 171, 170, 122, 123, 127, 126

87, 167, 168, 172, 171, 123, 124, 128, 127

88, 169, 170, 174, 173, 125, 126, 130, 129

89, 170, 171, 175, 174, 126, 127, 131, 130

90, 171, 172, 176, 175, 127, 128, 132, 131

*Nset, nset=_PickedSet2, internal, generate

1, 176, 1

*Elset, elset=_PickedSet2, internal, generate

1, 90, 1

** Section: Section-1

*Solid Section, elset=_PickedSet2, controls=EC-1, material=Material-1

*End Part

**

**

** ASSEMBLY

**

*Assembly, name=Assembly

**

*Instance, name=Part-1-1, part=Part-1

*End Instance

**

*Nset, nset=_PickedSet5, internal, instance=Part-1-1

1, 2, 3, 4, 45, 46, 47, 48, 89, 90, 91, 92, 133, 134, 135, 136

*Elset, elset=_PickedSet5, internal, instance=Part-1-1

1, 2, 3, 31, 32, 33, 61, 62, 63

*Elset, elset=Set-1, instance=Part-1-1

28,

*Elset, elset=__PickedSurf6_S5, internal, instance=Part-1-1

28, 29, 30, 58, 59, 60, 88, 89, 90

*Surface, type=ELEMENT, name=_PickedSurf6, internal

__PickedSurf6_S5, S5

*End Assembly

**

** ELEMENT CONTROLS

**

*Section Controls, name=EC-1, hourglass=ENHANCED

1., 1., 1.

**

** MATERIALS

**

*Material, name=Material-1

*Depvar

15,

*User Material, constants=6

200000., 0.3,400.,700., 0.5,400.

**

** BOUNDARY CONDITIONS

**

** Name: BC-1 Type: 位移/转角

*Boundary

_PickedSet5, 1, 1

_PickedSet5, 2, 2

_PickedSet5, 3, 3

_PickedSet5, 4, 4

_PickedSet5, 5, 5

_PickedSet5, 6, 6

** ----------------------------------------------------------------

**

** STEP: Step-1

**

*Step, name=Step-1

*Static, direct

0.01, 1.,

**

** LOADS

**

** Name: Load-1 Type: Pressure

*Dsload

_PickedSurf6, P, -420.

**

** OUTPUT REQUESTS

**

*Restart, write, frequency=0

**

** FIELD OUTPUT: F-Output-1

**

*Output, field, variable=PRESELECT

**

** HISTORY OUTPUT: H-Output-1

**

*Output, history

*Element Output, elset=Set-1

E33, PEEQ, PEEQT, S33

*End Step

(4条)
默认 最新
评论 点赞 1
厉害,
评论 点赞

查看更多评论 >

点赞 18 评论 5 收藏 158
关注