UMAT全过程技术篇

UMAT全过程-"技术篇"[写在前面:

这篇文章是UMAT全过程-感想篇的姊妹篇,容许要给大家写的一篇帖子,同时也是为了记录自己的学习过程,与大家分享!首先指出,俺的"技术篇"--是加了引号的,因为确实称不上有多么大的技术含量,还望大家莫笑偶!只不过一是跟那个感想篇形成一个对照,同时主要内容为自己编子程序过程中涉及的技术边边上的小问题的一些解决方法,供仿友们参考!偶不是谦虚,也不是一个低调的人,大家谢谢和支持的话,我先行谢过啦!更希望大家能提出质疑或者别的更好的方法,大家相互交流,共同进步!]----------------------------------------------------------------------------*转*入*正*题*第一局部:相关知识[特别声明,这局部来自于华中科技大学杨曼娟同学的硕士学位论文,在此对作者表示感谢!--大家可以去知网下载]----------------------------------------------------------------------------1.ABAQUS中材料非线性问题的处理

ABAQUS中材料非线性问题用Newton-Raphson法来求解。首先将载荷分为假设干个微小增量，结构受到一个微小增量△P。ABAQUS用与初始结构位移相对应的初始刚度矩阵K0和荷载增量△P计算出结构的在这一步增量后的位移修正Ca、修正后的位移值Ua和相应的新的刚度矩阵Ka。ABAQUS用新的刚度矩阵计算结构的内力Ia，荷载P和Ia的差值为迭代的剩余力Ra，即Ra=P-Ia。如果Ra在模型内的每个自由度上的值都为零，如图2-2中的a点，那么结构处于平衡状态。但在非线性问题中，通常Ra是不可能为零，ABAQUS为此设置了一个剩余力容差。如果Ra小于这个数字，ABAQUS就认为结构的内外力是平衡的。一般这个缺省值取为平均内力的0.5％(如图2-2)。

在Ra满足小于剩余力容差的条件后，ABAQUS还要检查位移修正系数Ca是否远远小于结构的位移增量△Ua=Ua-Uo.如果大于位移增量的1％，ABAQUS将自动继续进行迭代。第二次迭代时，ABAQUS用刚度矩阵Ka和第一次迭代的剩余力Ra计算，得到一个新的位移修正值。同理，第二次迭代中ABAQUS用新的剩余力、新的位移系数，继续验证收敛性，直到收敛为止〔如图2-3〕。因此非线性问题的计算量要比线性问题的计算量大得多，计算的时候需要分配更多的内存和磁盘空间。

ABAQUS会自动调整荷载增量的大小，用户只需合理的设置一个第一次迭代的增量值即可，如果用户不给出第一次迭代的增量值，ABAQUS会在第一次迭代的过程中把设置的所有荷载都加到结构上去，然后不断的试算自动调整第一次迭代的增量大小。寻找一个增量值的收敛解迭代的次数取决于系统非线性的程度。ABAQUS中缺省的设置是，如果在一个增量值的作用下，迭代16次以后结果仍不收敛，ABAQUS将放弃这个增量值，把增量取为前一次增量值的25％再进行计算。如果ABAQUS在迭代次数小于5次的情况下就找到了收敛解，ABAQUS就自动把增量值增加50％再进行计算，得到更为合理的增量值。2.用户子程序接口

ABAQUS不仅提供标准的有限元分析程序，而且具有良好的开放性，可利用它

提供的用户子程序接口生成非标准的分析程序来满足用户的需要，在实际工程中得

到广泛应用。ABAQUS允许用户通过子程序以代码的形式来扩展主程序的功能，并

给用户提供了强大而又灵活的用户子程序接口和应用程序接口〔UtilityRoutine〕，

ABAQUS6.4一共有42个用户子程序接口，15个应用程序接口，用户可以定义包括

边界条件、荷载条件、接触条件、材料特性以及利用用户子程序和其它应用软件进

行数值交换等等。这些用户子程序接口使得用户解决一些问题时有很大的灵活性，

同时大大的扩充了ABAQUS的功能。例如：

通过用户定义单元接口，用户自定义的任何类型的线性或非线性单元都可以被引入模型中，对于线性单元刚度矩阵和质量矩阵可以直接确定。同时，用户子程序也可被用来定义这些单元的线性和非线性特性。通过用户材料子程序接口，用户可定义任何补充的材料模型，不但任意数量的

材料常数都可以作为资料被读取，而且ABAQUS对于任何数量的与解相关的状态变

量在每一材料计数点都提供了存储功能，以便在这些子程序中应用。3.用户子程序和主程序的结合

ABAQUS的用户子程序是根据ABAQUS提供的相应接口，按照Fortran语法，用户自己编写的代码。它是一个独立的程序单元，可以独立的被存储和编译，也能被其它程序单元引用，因此，利用它可带回大量数据供引用程序使用，也可以用它来完成各种特殊的功能。它的一般结构形式是：SUBROUTINES(x1,x2,……,xn)

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

ORINCLUDE‘VABA_PARAM.INC’)〔用于ABAQUS/Explicit用户子程序中〕

……

RETURN

END

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

在一个算例中，用户可以用到多个用户子程序，但必须把它们放在一个以.for为

扩展名的文件中。运行带有用户子程序的算例同时有两种方法：一是在CAE中运行，

在EDITJOB菜单中的GENERAL子菜单的USERSUBROUTINEFILE对话框中选择用户子程序所在的文件即可；另外是在ABAQUS.COMMAND中运行，语法如下：

abaqusjob=job-nameuser={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的相关规定.4.用户材料子程序UMAT接口的原理

用户材料子程序〔User-definedMaterialMechanicalBehavior，简称UMAT〕是

ABAQUS提供应用户定义自己的材料属性的Fortran程序接口，使用户能使用ABAQUS材料库中没有定义的材料模型。用户材料子程序UMAT通过与ABAQUS主求解程序的接口实现与ABAQUS的资料交流.在输入文件中，使用关键词“*USERMATERIAL〞表示定义用户材料属性。[更详细相关知识,请看帖子ABAQUS-UMAT-自学知识整理贴[已经初步完成,不断完善更新,请跟帖讨论]]

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

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

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

ABAQUS调用UMAT过程及次数如图2-4所示。5.UMAT子程序流程

UMAT子程序采用Fortran语言编制，可以包括以下几个局部：子程序定义语句、ABAQUS定义的参数说明、用户定义的局部变量说明、用户编制的程序主体、子程序返回和结束语句。

主要求解过程：每一个增量加载步开始时，ABAQUS主程序在单元的积分点上调用UMAT子程序，并传入应变增量、时间步长及荷载增量，同时也传入当前状态的应力、应变及其它与求解过程相关的变量；UMAT子程序根据本构方程求解应力增量并更新应力及其它相关的变量，提供Jacobian矩阵给ABAQUS主程序以形成整体刚度矩阵；主程序结合当前荷载增量求解位移增量，继而进行平衡校核；如果不满足指定的误差，ABAQUS将进行迭代直到认为收敛，然后进行下一增量步的求解。

ABAQUS中应力σ采用Cauchy应力张量描述，剪应变分量εij按照工程剪应变

的定义存储。在求解的过程中，需记录保存与求解过程相关的变量，这可以储存在UMAT子程序的STATEV数组中，并且可以定义储存的数量。对本文所述的本构方程，定义了13个与求解过程相关的变量，即6个弹性应变，6个塑性应变，1个塑性参数增量。非非的ABAQUS用户子程序[zz]转自当用到某个用户子程序时，用户所关心的主要有两方面：一是ABAQUS提供的用户子程序的接口参数。有些参数是ABAQUS传到用户子程序中的，例如SUBROUTINEDLOAD中的KSTEP，KINC，COORDS；有些是需要用户自己定义的，例如F。二是ABAQUS何时调用该用户子程序，对于不同的用户子程序ABAQUS调用的时间是不同的。有些是在每个STEP的开始，有的是STEP结尾，有的是在每个INCREMENT的开始等等。当ABAQUS调用用户子程序是，都会把当前的STEP和INCREMENT利用用户子程序的两个实参KSTEP和KINC传给用户子程序，用户可编个小程序把它们输出到外部文件中，这样对ABAQUS何时调用该用户子程序就会有更深的了解。〔子程序中很重要的就是要知道由abaqus提供的那些参量的意义，如下〕首先介绍几个子程序：一．SUBROUTINEDLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,COORDS,JLTYP,SNAME)参数：

1．F为用户定义的是每个积分点所作用的荷载的大小；

2．KSTEP,KINC为ABAQUS传到用户子程序当前的STEP和INCREMENT值；

3．TIME(1),TIME(2)为当前STEPTIME和INCREMENTTIME的值；

4．NOEL,NPT为积分点所在单元的编号和积分点的编号；

5．COORDS为当前积分点的坐标；

6．除F外，所有参数的值都是ABAQUS传到用户子程序中的。功能：

1．荷载可以被定义为积分点坐标、时间、单元编号和单元节点编号的函数。

2．用户可以从其他程序的结果文件中进行相关操作来定义积分点F的大小。例1：这个例子在每个积分点施加的荷载不仅是坐标的函数，而且是随STEP变化而变化的。

SUBROUTINEDLOAD(P,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,COORDS,

1JLTYP,SNAME)

INCLUDE'ABA_PARAM.INC'C

DIMENSIONTIME(2),COORDS(3)

CHARACTER*80SNAME

PARAMETER(PLOAD=100.E4)

IF(KSTEP.EQ.1)THEN！当STEP=1时的荷载大小

P=PLOAD

ELSEIF(KSTEP.EQ.2)THEN！当STEP＝2时的荷载大小

P=COORDS(1)*PLOAD！施加在积分点的荷载P是坐标的函数

ELSEIF(KSTEP.EQ.3)THEN！当STEP=3时的荷载大小

P=COORDS(1)**2*PLOAD

ELSEIF(KSTEP.EQ.4)THEN！当STEP=4时的荷载大小

P=COORDS(1)**3*PLOAD

ELSEIF(KSTEP.EQ.5)THEN！当STEP=5时的荷载大小

P=COORDS(1)**4*PLOAD

ENDIF

RETURN

END

UMAT子程序具有强大的功能，使用UMAT子程序：

(1)

可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进行计算，扩充程序功能。

(2)几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予ABAQUS中的任何单元；(3)必须在UMAT中提供材料本构模型的雅可比〔Jacobian〕矩阵，即应力增量对应变增量的变化率。

(4)可以和用户子程序“USDFLD〞联合使用，通过“USDFLD〞重新定义单元每一物质点上传递到UMAT中场变量的数值。由于主程序与UMAT之间存在数据传递，甚至共用一些变量，因此必须遵守有关UMAT的书写格式，UMAT中常用的变量在文件开头予以定义，通常格式为：

SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

1RPL,DDSDDT,DRPLDE,DRPLDT,

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

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

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

INCLUDE'ABA_PARAM.INC'

CHARACTER*80CMNAME

DIMENSIONSTRESS(NTENS),STATEV(NSTATV),

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

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

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

usercodingtodefineDDSDDE,STRESS,STATEV,SSE,SPD,SCD

and,ifnecessary,RPL,DDSDDT,DRPLDE,DRPLDT,PNEWDT

RETURN

END

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

DDSDDE〔NTENS，NTENS〕是一个NTENS维的方阵，称作雅可比矩阵，，是应力的增量，是应变的增量，DDSDDE〔I，J〕表示增量步结束时第J个应变分量的改变引起的第I个应力分量的变化。通常雅可比是一个对称矩阵，除非在“*USERMATERIAL〞语句中参加了“UNSYMM〞参数。

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

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

PROPS〔NPROPS〕材料常数矩阵，矩阵中元素的数值对应于关键字“*USERMATERIAL〞下面的数据行。

SSE，SPD，SCD分别定义每一增量步的弹性应变能，塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅作为能量输出。其他变量：

STRAN(NTENS)：应变矩阵；DSTRAN(NTENS)：应变增量矩阵；

DTIME：增量步的时间增量；

NDI：直接应力分量的个数；

NSHR：剪切应力分量的个数；

NTENS：总应力分量的个数，NTENSNDINSHR=+。使用UMAT时需要注意单元的沙漏控制刚度和横向剪切刚度。通常减缩积分单元的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。这些刚度基于材料初始剪切模量的值，通常在材料定义中通过“*ELASTIC〞选项定义。但是使用UMAT的时候，ABAQUS对程序输入文件进行预处理的时候得不到剪切模量的数值。所以这时候用户必须使用“*HOURGLASSSTIFFNESS〞选项来定义具有沙漏模式的单元的沙漏控制刚度，使用“*TRANSVERSESHEARSTIFFNESS〞选项来定义板、壳、梁单元的横向剪切刚度。几个关于子程序的问题及相应解答

Q:本人在用umat作本构模型时，

*static,

1,500,0.000001,0.1此时要求的增量步很多，即每次增量要很小，

*static

1,500时，在弹性向塑性过度时，出现错误，增量过大，出现尖点.?

A:YOUCANTRYASFOLLOWS:

*STEP,EXTRAPOLATION=NO,INC=2000000

*STATIC

0.001,500.0,0.00001,0.1。

Q:在abaqus中，如果采用umat,利用自己的本构，如何让abaqus明白这种材料的弹塑性应变，也就是说，如何让程序返回弹性应变与塑性应变，好在output中输出，我曾想用最笨地方法，在uvarm中定义输出，利用getvrm获取材料点的值，但无法获取增量应力，材料常数等，研究了帮助中的例子，umatmst3.inp,umatmst3.for,他采用misesJ2流动理论，我在outputhistory显示他已进入塑性状态，但他的PE仍然为0！！?

A:用uvar()勉强成功。

Q:偶在umat中调用求主应力函数

CALLSPRINC(STRESS,PS,LSTR,NDI,NSHR)后，存储主应力得数组PS中各个主应力排列顺序是什么？

PS1>PS2>PS3?

PS1<PS2<PS3?

PS1>PS3>PS2?

A:第二个。个人觉得：

umat实现自己的本构没有固定的方法，对于不同的本构有可能必须采用不同的方法。这要靠自己不断地摸索。有可能一种方法对于简单加载问题还行，但有可能对于复杂问题并不收敛。最重要一点，就是umat中采用的算法必须consistent.再就是ddsdde必须正确，〔如果采用back_Euler方法等一些算法，ddsdde错误有时不影响结果〔对于简单加载问题没有影响，能收敛，〕，但对于复杂问题不收敛。

uptonow，你这个算法对于Mises,hill,J2，J2d等一类的屈服函数是正确的，但具体的本构还要灵活运用，这我也正学习，正在摸索。有时，umat需要很强的有限元根底，并且对采用的本构要很熟悉，不要在一颗树上吊死才好。首先要确认自己的umat没有错误，如果没有，但就是不收敛〔在不断减小加载步长的情况下，当然最好对步长不敏感，特别是对于粘弹性，粘塑性，内变量一类的材料，有的本构取决于背应力的计算〕。那就应该考虑换一种算法。一点体会，请大家探讨。

Q:abaqus-uamt的老问题，缺少'ABA_PARAM.INC'文件?

A:在cvf6.5调试时，显示缺少'ABA_PARAM.INC'文件！这个没有任何关系的，这个错误将在ABAQUS调用UMAT的时候自动会找到，仅仅有这个错误将没有任何影响的。也就是说，ABAQUS中调用的时候，实际并不存在这个错误。FT，忘了说一句了，你把或拷到你的程序工作空间后，应该将把或的后缀.sp或.dp去掉，即将ABA_PARAM.INC.dp或改名为ABA_PARAM.INC。呵呵，他的意思是在VisualFortran中调试其子程序，我觉得这是一个好方法，我当时也是这麽办的，毕竟在ABAQUS中调试是非常麻烦的，只有当你的UMAT没有语法或者明显的逻辑错误，你在ABAQUS中调试才能事半功倍。

Q：uvarm可以输出到哪里？.odb可以么？另外那个strav？？就是自己定义用于umat的那个数组里的数可以输出到.odb里么？

A：在umat中，statev是不能用在output中的，statev只是作为一个解的状态变量，说来惭愧，我是在umat中定义peeq,(peeq)的求解一般在弹塑性力学书上有〔等效塑性应变〕，用write()写入一个临时文件，〔〔切记：这个文件unit号不要与abaqus中的重合，因为他有一些系统默认的文件号，〕〕然后在uvarm中读取，以uvarm输出，因为uvarm可以以odb的形式输出，支持output,field,output,history

Q：在本版看了一个一维固结的例子，其中含有用户子程序，如下

SUBROUTINEUFIELD(FIELD,KFIELD,NSECPT,KSTEP,KINC,TIME,NODE,COORDS,TEMP,DTEMP)

INCLUDE'ABA_PARAM.INC'

DIMENSIONFIELD(NSECPT),TIME(2),COORDS(3),TEMP(NSECPT),DTEMP(NSECPT)

KFIELD=1

FIELD(1)=COORDS(2)

RETURN

END我在VF6.5中进行调试，提示找不到ABA_PARAM.INC请问大侠这如何解决。还有我因为是初次接触用户子程序，我查阅了本版所有的相关贴子，都讲的不太详细，我将问题总结一下，大侠们能不能详细的讲解一下，

1子程序格式(程序后缀是.f;.f90;.for;.obj??)

2CAE中如何调用，command下如何调用？

3假设有多个子程序同时存在，如何处理

4我对VF不是很熟，是否可以用VC，C＋＋编写子程序？

A:假设要在vf中调试，那么应该根据需要把SITE文件夹中的ABA_PARAM_DP.INC〔双精度〕或ABA_PARAM_SP.INC〔单精度〕拷到相应的位置，并改名为ABA_PARAM.INC即可。

1。我试过，.for格是应该是不可以的，至少6.2和6.3版本应该是不行，其他的没用过，没有发言权。在Abaqus中，运行abaqusj=jobnameuser=username时，默认的用户子程序后缀名是.for〔.f,.f90应该都不行的，手册上也有讲过〕，只有在username.for文件没有找到的情况下，才会去搜索username.