umat二次开发超弹性本构

合集下载

各向异性屈服准则的UMAT子程序二次开发研究

各向异性屈服准则的UMAT子程序二次开发研究乔顺成;吴建军;展学鹏【摘要】通过ABAQUS提供的UMAT子程序二次开发接口,结合完全隐式向后Euler图形返回算法,用Fortran语言编程,将Yld2004-18p各向异性屈服准则本构模型嵌入ABAQUS软件,并提出了一个通用的、柔性的本构模型二次开发结构模式.将Yld2004-18p屈服准则退化到Mises各向同性屈服准则,用于实例有限元分析,计算结果与ABAQUS自带的Mises屈服准则计算结果做比较,验证了二次开发过程的正确性,为后续嵌入更高级、更强健、更特殊的屈服准则提供研究思路和方法,使之用于精密塑性有限元分析.【期刊名称】《锻压装备与制造技术》【年(卷),期】2018(053)004【总页数】9页(P89-97)【关键词】UMAT子程序;各向异性;屈服准则;本构模型;二次开发【作者】乔顺成;吴建军;展学鹏【作者单位】中航工业西安飞机工业集团有限责任公司,陕西西安710072;西北工业大学机电学院,陕西西安710072;西北工业大学机电学院,陕西西安710072【正文语种】中文【中图分类】TP391.77ABAQUS以其强大的非线性迭代计算和前后处理功能而广泛应用于材料弹塑性有限元分析，它是现阶段在各个工程领域广泛使用的大型通用有限元软件之一。

ABAQUS[1]有大量的单元库和求解模型供用户使用，而且用户也能够通过这些模型求解绝大多数问题。

但是实际问题是相当复杂的，ABAQUS不可能直接处理所有可能的问题[2-4]，所以，有必要给用户提供二次开发的接口来解决实际中出现的新问题。

在用户自定义材料（UMAT即 User-defined Material Mechanical Behavior）研究中，很多学者通过UMAT子程序二次开发，将新的本构模型用于实际问题的有限元分析。

例如：李平[5]等通过数值分析模拟了普通碳钢连续热轧的过程，将率相关的各向同性硬化热变形本构模型通过UMAT子程序用于有限元仿真，分析了轧制过程中的温度和应力、应变之间的耦合关系。

umat二次开发超弹性本构()

APPENDIXNeo-Hookean Hyperelatic Material User SubroutineThis program is based on the derivation of hyperelastic material constitutive model inSection 4.4. A stress and strain relationship was derived from the neo-Hookean hyperelastic material constitutive model that is normally represented as the strain energywith strain invariants.subroutine vumat(C Read only (unmodifiable)variables -1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only (modifiable) variables -7 stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'Cdimension props(nprops), density(nblock), coordMp(nblock,*),1 charLength(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),8 defgradNew(nblock,ndir+nshr+nshr),9 fieldNew(nblock,nfieldv),1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),2 enerInternNew(nblock), enerInelasNew(nblock)Ccharacter*80 cmnameCif (cmname(1:6) .eq. 'VUMAT0') thencall VUMAT0(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,117else if (cmname(1:6) .eq. 'VUMAT1') thencall VUMAT1(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)end ifendCsubroutine vumat0 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -Cinclude 'vaba_param.inc'Cdimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0)real C10,D1,ak,twomu,amu,alamda,hydro,vonMises, maxShear,1 midStrain, maxPrincipalStrain118Cintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '//* 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1191 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two) scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1) stateNew(k,2) = stateOld(k,2) + strainInc(k,2) stateNew(k,3) = stateOld(k,3) + strainInc(k,3) stateNew(k,4) = stateOld(k,4) + strainInc(k,4) CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+ 1201 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C use 3/2 will get 2 (int) !CC write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4) CC Failure CriteriaCmidStrain = stateNew(k,1) + stateNew(k,2)maxShear = sqrt((stateNew(k,1) - midStrain)**2. +1 stateNew(k,4)**2.)if (midStrain .GE. 0.) thenmaxPrincipalStrain = midStrain + maxShearelsemaxPrincipalStrain = maxShear - midStrainend ifif (vonMises .GE. 10.8565e6) thenstateNew(k,5) = 0end ifend doend ifreturnendCsubroutine vumat1 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,121* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'dimension coordMp(nblock,*), charLength(nblock), props(nprops), 1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0)real C10,D1,ak,twomu,amu,alamda,hydro,vonMisesCintv(1) = ndirintv(2) = nshrif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '//* 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10122twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3) stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two) scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)123EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1) stateNew(k,2) = stateOld(k,2) + strainInc(k,2) stateNew(k,3) = stateOld(k,3) + strainInc(k,3) stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4)124 end do end if return end。

ABAQUS-二次开发资料-UMAT

各个楼层及内容索引2-------------------------------------什么是UMAT3-------------------------------------UMAT功能简介4-------------------------------------UMAT开始的变量声明5-------------------------------------UMAT中各个变量的详细解释6-------------------------------------关于沙漏和横向剪切刚度7-------------------------------------UMAT流程和参数表格实例展示8-------------------------------------FORTRAN语言中的接口程序Interface9-------------------------------------关于UMAT是否可以用Fortran90编写的问题10-17--------------------------------Fortran77的一些有用的知识简介20-25\30-32-----------------------弹塑性力学相关知识简介34-37--------------------------------用户材料子程序实例JOhn-cook模型压缩包下载38-------------------------------------JOhn-cook模型本构简介图40-------------------------------------用户材料子程序实例JOhn-cook模型完整程序+david详细注解[欢迎大家来看看,并提供意见,完全是自己的diy的,不保证完全正确,希望共同探讨,以便更正,带""部分,还望各位大师\同仁指教]1什么是UMAT1.1UMAT功能简介!!![-摘自庄茁老师的书UMAT子程序具有强大的功能，使用UMAT子程序：(1)可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进行计算，扩充程序功能。

ABAQUS材料用户子程序UMAT学习报告

具有友好的用户界面和易用的操作流程，方便用户进行学习和使用
UMAT子程序简介
UMAT子程序是 ABAQUS材料用户自定义模块，允许用户根据实际需求编写材料本构模型。
UMAT子程序采用C 语言编写，用户需要具备一定的编程基础。
UMAT子程序可以实现多种材料本构模型，如弹性、塑性、蠕变等。
UMAT子程序实现细节
编程语言和接口
A B A Q U S 材料用户子程序 U M AT 使用Fortran语言编写
U M AT 子程序中可以定义材料属性、本构关系等
添加标题
添加标题
添加标题
添加标题
U M AT 子程序通过 A B A Q U S 提供的接口与主程序进行交互
不足：使用门槛较高，需要用户具备一定的编程基础
未来展望：期待更多的用户参与开发，不断完善子程序功能
总结：UMAT子程序为用户提供了强大的材料模型描述能力，但使用过程中需要注意其局限性
在ABAQUS中的未来发展方向
开发更高效的材料模型集成人工智能和机器学习技术增强与CAD软件的集成扩展对多物理场模拟的支持
适用于金属材料
适用于复合材料
适用于橡胶材料
适用于陶瓷材料
参数的合理选择
参数选择需符合实际物理模型参数选择需考虑材料特性参数选择需经过实验验证参数选择需注意收敛性和稳定性
收敛性和稳定性问题
收敛性： U M AT 子程序在迭代过程中应满足收敛条件，否则可能导致计算失败或结果不准确。
边界条件和初始条件
边界条件：描述模型在边界上的行为，如位移、速度等初始条件：描述模型在初始时刻的状态，如温度、压力等

形状记忆合金本构的VUMAT二次开发

a a b
(1-3a)
卸载时：
σ b E A ε-ε b Ε Α εd Ε A ε εd Ε ε Α
其中: εc b a d
c b d c d
止应变； d 为再次按奥氏体弹性模量
卸载的开始应变。
加载段分为两段：
段 d-o 段，则再加载的第一加载段和第二加载段的分界点将是原来的 a 点。双旗型-SMA 超弹性本构模型为压缩和拉伸的本构关系对称的形式，如图 1-1 所示。

a
ΕΑ ΕΑ
PR OP S 物理性系系 1 2 3 4 5 6 7 8

-T
数

-T
a d
-T -T

-C
数

-C
a
-C
d
-C
A
ΕΑ
b
质
注：T -拉伸，C-压缩
c
d
A
a
SMA-VUMA T 子程序的主要编制思路，如下：
εb
o εd
c
1) 以应变是否大于零，来区分拉伸和压缩两种情况。 2) 以应变增量是否大于零，来区分加载和卸载两种情况。 3) 以关键点的应力值作为加载段 \ 卸载段的划分依据，对不同的加载段 \卸载段进行划分，根据当前传递过来的应力值，判断其所在的加载 \ 卸载段，定义新的应力值计算公式。最终编制成功的 SMA-VUMA T 子程序，经过单个单元和构件的双重测试，测试情况多样，全面，概况如下： 1) 单纯拉伸\压缩：加载进入第一加载段或第二加载段的各种加载再卸载情况。 2) 单纯拉伸\压缩：卸载进入第一加载段、第二加载段或第三卸载段的各种加载、卸载再加载的情况。拉伸与压缩各种混合加载\ 卸载情

UMAT以及depvar整理

UMAT以及depvar整理UMAT以及depvar整理by warmwormdk最近一直在论坛查资料，对自己感兴趣的一些问题专门进行了整理，希望对大家有所帮助，也希望能获得小小的奖励啊，哈哈1UMAT的状态变量问题Q:用DEPVAR定义的状态变量个数假设为10个。

是不是说一个积分点的状态变量是10个，单元的积分点是4的话，那么单元的状态变量就是40个。

也就是自己要存储单元变量的话，就得按40个状态变量来。

是不是呢？A:有人说跟我说，DEPVAR定义的状态变量个数指每个积分点的状态变量个数。

abaqus 会自动为每个单元的每个积分点开辟这样大小的状态变量数组，abaqus调用umat时能够自动根据单元好和积分点向umat提供状态变量，在此基础上umat修改状态变量。

2umat里的STATEV变量怎么输出到odb文件中Q:比如我想知道statev(10)的odb文件,怎么输出?又怎么打开.statev(10)是我自己定义的damage变量,请各位赐教A在element关键词中添加SDV，就像下面这样。

*Element Output,directions=YESALPHA,LE,PE,PEEQ,PEMAG,PS,S,STH,SE,SF,VE,VEEQ,VS,SDV3Q statev(?)请问这个状态变量（）内的数字代表什么含义？对应的变量是不是固定的？各自对应着哪些变量？A:括号中的数代表这个变量矩阵的维数，这个值等于depvar的值。

4umat中DEPVAR有几种定义方式?Q;UMAT中状态变量的定义方式,一般有两种形式,一是在inp文件中采用*initial conditions定义,二是特殊情况下可以采用SDVINI来定义;目前的疑问是,是否还有其他定义状态变量的方式?请专家针对上传的附件给于指点多谢!附件中的例子来自ABAQUS HELP中的例子!请问例子中INP文件中是如何定义状态变量的THANKS A:初值可以采用SDVINI来定义程序运行中用以下代码更新DO310K1=1,NTENSSTATEV(K1)=EELAS(K1)STATEV(K1+NTENS)=EPLAS(K1)310CONTINUESTATEV(1+2*NTENS)=EQPLASCRETURNEND5umat子程序定义问题？Q:在umat中定义的参数在材料中需不需再定义了？比如我umat 中已经给了密度了。

umat二次开发超弹性本构

APPENDIXNeo-Hookean Hyperelatic Material User SubroutineThis program is based on the derivation of hyperelastic material constitutive model in Section 4.4. A stress and strain relationship was derived from the neo-Hookean hyperelastic material constitutive model that is normally represented as the strain energy with strain invariants.subroutine vumat(C Read only (unmodifiable)variables -1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only (modifiable) variables -7 stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'Cdimension props(nprops), density(nblock), coordMp(nblock,*),1 charLength(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),8 defgradNew(nblock,ndir+nshr+nshr),9 fieldNew(nblock,nfieldv),1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),2 enerInternNew(nblock), enerInelasNew(nblock)Ccharacter*80 cmnameCif (cmname(1:6) .eq. 'VUMAT0') thencall VUMAT0(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)117else if (cmname(1:6) .eq. 'VUMAT1') thencall VUMAT1(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)end ifendCsubroutine vumat0 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'Cdimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0)real C10,D1,ak,twomu,amu,alamda,hydro,vonMises, maxShear,1 midStrain, maxPrincipalStrain118Cintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '//* 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1191 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two) scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1)stateNew(k,2) = stateOld(k,2) + strainInc(k,2)stateNew(k,3) = stateOld(k,3) + strainInc(k,3)stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1201 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C use 3/2 will get 2 (int) !CC write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4) CC Failure CriteriaCmidStrain = stateNew(k,1) + stateNew(k,2)maxShear = sqrt((stateNew(k,1) - midStrain)**2. +1 stateNew(k,4)**2.)if (midStrain .GE. 0.) thenmaxPrincipalStrain = midStrain + maxShearelsemaxPrincipalStrain = maxShear - midStrainend ifif (vonMises .GE. 10.8565e6) thenstateNew(k,5) = 0end ifend doend ifreturnendCsubroutine vumat1 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,121* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'dimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0)real C10,D1,ak,twomu,amu,alamda,hydro,vonMisesCintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '//* 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10122twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in localaxis)Cdet=stretchNew(k, 3)*1 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two)scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)123EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1)stateNew(k,2) = stateOld(k,2) + strainInc(k,2)stateNew(k,3) = stateOld(k,3) + strainInc(k,3)stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4) 124end doend ifreturnend。

ABAQUS-二次开发资料-UMAT

各个楼层与容索引2-------------------------------------什么是UMAT3-------------------------------------UMAT功能简介4-------------------------------------UMAT开场的变量声明5-------------------------------------UMAT中各个变量的详细解释6-------------------------------------关于沙漏和横向剪切刚度7-------------------------------------UMAT流程和参数表格实例展示8-------------------------------------FORTRAN语言中的接口程序Interface9-------------------------------------关于UMAT是否可以用Fortran90编写的问题10-17--------------------------------Fortran77的一些有用的知识简介20-25\30-32-----------------------弹塑性力学相关知识简介34-37--------------------------------用户材料子程序实例JOhn-cook模型压缩包下载38-------------------------------------JOhn-cook模型本构简介图40-------------------------------------用户材料子程序实例JOhn-cook模型完整程序+david详细注解[欢送大家来看看,并提供意见,完全是自己的diy的,不保证完全正确,希望共同探讨,以便更正,带"?"局部,还望各位大师\指教]1 什么是UMAT???1.1 UMAT功能简介!!![-摘自庄茁教师的书UMAT子程序具有强大的功能，使用UMAT子程序：(1)可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进展计算，扩大程序功能。

umat二次开发超弹性本构审批稿

u m a t二次开发超弹性本构YKK standardization office【 YKK5AB- YKK08- YKK2C- YKK18】APPENDIXNeo-Hookean Hyperelatic Material User SubroutineThis program is based on the derivation of hyperelastic material constitutive model in Section . A stress and strain relationship was derived from the neo-Hookean hyperelastic material constitutive model that is normally represented as the strain energy with strain invariants.subroutine vumat(C Read only (unmodifiable)variables -1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only (modifiable) variables -7 stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude ''Cdimension props(nprops), density(nblock), coordMp(nblock,*),1 charLength(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),8 defgradNew(nblock,ndir+nshr+nshr),9 fieldNew(nblock,nfieldv),1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),2 enerInternNew(nblock), enerInelasNew(nblock)Ccharacter*80 cmnameCif (cmname(1:6) .eq. 'VUMAT0') thencall VUMAT0(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)117else if (cmname(1:6) .eq. 'VUMAT1') thencall VUMAT1(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)end ifendCsubroutine vumat0 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude ''Cdimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = , one = , two = , three = ,* four = , half =real C10,D1,ak,twomu,amu,alamda,hydro,vonMises, maxShear, 1 midStrain, maxPrincipalStrain118Cintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented 'Q. zero) then do k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3) stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1191 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two) scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1) stateNew(k,2) = stateOld(k,2) + strainInc(k,2) stateNew(k,3) = stateOld(k,3) + strainInc(k,3) stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1201 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C use 3/2 will get 2 (int) !CC write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4) CC Failure CriteriaCmidStrain = stateNew(k,1) + stateNew(k,2)maxShear = sqrt((stateNew(k,1) - midStrain)**2. +1 stateNew(k,4)**2.)if (midStrain .GE. 0.) thenmaxPrincipalStrain = midStrain + maxShearelsemaxPrincipalStrain = maxShear - midStrainend ifif (vonMises .GE. thenstateNew(k,5) = 0end ifend doend ifreturnendCsubroutine vumat1 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength, * props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,121* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude ''dimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = , one = , two = , three = ,* four = , half =real C10,D1,ak,twomu,amu,alamda,hydro,vonMisesCintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented 'Q. zero) then do k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3) stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two) scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)123EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PRstressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1) stateNew(k,2) = stateOld(k,2) + strainInc(k,2) stateNew(k,3) = stateOld(k,3) + strainInc(k,3) stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4) 124end doend ifreturnend。

5.本构模型-应力更新专题-UMAT和VUMAT

3.前推后拉及Lie导数
Lie导数
后拉和前推的概念为定义张量的时间导数提供了数学上的一致性－Lie导数。如框5.17，Kirchhoff应力的Lie导数是其应力的后拉的时间导数的前推。
D D Lv τ * ( * τ) F (F1 τ) Dt Dt
Lagrangian矢量dX和Eulerian矢量dx定义的二阶张量
可以由后拉和前推运算给出E-L张量之间映射的统一描述。例如，L矢量dX由F前推到当前构形给出E矢量dx
dx F dX * dX
Eulerian-Lagrangian 前推运算
1 E矢量dx由 F 后拉到参考构形给出L矢量dX
Jaumann率发生有限剪切时慎用Jaumann率
Green-Naghdi率
2.几种客观率的关系
如何得到正确的结果？切线模量之间的关系
次弹性本构关系共同应用的形式为
σ J CsJ : D
σ T CsT : D
σ G CsG : D
对于各向同性材料Jaumann率的切线模量为
退化
Ω σ σ ΩT σG σ
=W 简化
Assume
W σ σ WT σJ σ
2.几种客观率的关系
总结一： Comparison of different objective stress rate
Truesdell •Difficult to implement •Not used in commercial software •Kinematically consistent with the rate of Cauchy stress •Must accurately determine the rotation tensor R •Relatively easy to implement •Produces symmetric tangent moduli

abaqus umat于计算力学之应用

abaqus umat是一种在计算力学中广泛应用的有限元分析软件。

它可以通过用户自定义的子程序（也称为umat）进行材料本构关系的定义，使得在模拟复杂材料行为时能够更加精确地描述材料的非线性和非均匀性等特性。

abaqus umat能够有效地模拟材料的机械性能，并在工程领域具有广泛的应用。

1. 什么是abaqus umat？abaqus umat是abaqus软件中用于用户自定义材料本构关系的子程序。

它可以实现对材料行为的精确描述，包括材料的非线性、非均匀性等特性。

通过abaqus umat，用户可以自定义材料的本构关系和材料参数，以满足对于各种材料行为的精确模拟需求。

2. abaqus umat的实现原理abaqus umat的实现依赖于有限元分析方法。

用户可以通过编写程序，在abaqus中调用该程序来定义材料的本构关系。

在有限元分析中，材料的本构关系是描述材料应力和应变之间关系的重要参数，通过用户自定义的umat程序，可以实现对材料行为的更为精确的描述。

3. abaqus umat的应用领域abaqus umat在工程领域有着广泛的应用。

例如在航空航天领域，abaqus umat可以用于模拟飞机结构的材料行为，预测飞机在不同载荷下的应力应变分布，进行疲劳分析等。

在汽车工业中，abaqus umat可以用于模拟汽车结构在碰撞时的材料行为，以及进行车身强度分析等。

abaqus umat还被广泛应用于建筑、船舶、能源等领域，在模拟复杂材料行为时发挥着重要作用。

4. abaqus umat的优势相较于其他有限元分析软件，abaqus umat的优势在于其灵活性和精确性。

用户可以通过编写自定义的umat程序，实现对材料行为的精确描述，满足各种复杂条件下的模拟需求。

abaqus umat还具有较强的兼容性和扩展性，可以与abaqus的其他模块结合使用，实现更为全面的分析和模拟。

5. 用户如何编写abaqus umat程序编写abaqus umat程序需要一定的编程和材料力学知识。

ABAQUS-UMAT弹塑本构二次开发地实现

word前言有限元法是工程中广泛使用的一种数值计算方法。

它是力学、计算方法和计算机技术相结合的产物。

在工程应用中，有限元法比其它数值分析方法更流行的一个重要原因在于：相对与其它数值分析方法，有限元法对边界的模拟更灵活，近似程度更高。

所以，伴随着有限元理论以与计算机技术的开展，大有限元软件的应用证变得越来越普与。

ABAQUS软件一直以非线性有限元分析软件而闻名，这也是它和ANSYS,Nastran等软件的区别所在。

非线性有限元分析的用处越来越大，因为在所用材料非常复杂很多情况下，用线性分析来近似已不再有效。

比方说，一个复合材料就不能用传统的线性分析软件包进展分析。

任何与时间有关联，有较大位移量的情况都不能用线性分析法来处理。

多年前，虽然非线性分析能更适合、更准确的处理问题，但是由于当时计算设备的能力不够强大、非线性分析软件包线性分析功能不够健全，所以通常采用线性处理的方法。

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

非线性有限元分析在各个制造行业得到了广泛应用，有不少大型用户。

航空航天业一直是非线性有限元分析的大客户，一个重要原因是大量使用复合材料。

新一代波音 787客机将全部采用复合材料。

只有像 ABAQUS这样的软件，才能分析包括多个子系统的产品耐久性能。

在汽车业，用线性有限元分析来做四轮耐久性分析不可能得到足够准确的结果。

分析汽车的整体和各个子系统的性能要求〔如悬挂系统等〕需要进展非线性分析。

在土木工程业， ABAQUS能处理包括混凝土静动力开裂分析以与沥青混凝土方面的静动力分析，还能处理高度复杂非线性材料的损伤和断裂问题，这对于大型桥梁结构，高层建筑的结构分析非常有效。

瞬态、大变形、高级材料的碰撞问题必须用非线性有限元分析来计算。

线性分析在这种情况下是不适用的。

以往有一些专门的软件来分析碰撞问题，但现在ABAQUS在通用有限元软件包就能解决这些问题。

ABAQUS材料用户子程序UMAT学习报告

NDI：某一点上直接应力组件数。 NSHR：某一点上剪切应力组件数。 NTENS：总应力分量个数，（=NDI+NSHR）。
NSTATV：存储与解有关的状态变量数组个数。 PROPS (NPROPS)：材料常数数组。 COORDS：当前积分点坐标。DROT(3,3) ：旋转增量矩阵。 CELENT：特征元素长度。
载入输入文件umat.inp，得到如下图形。
按下图所示操作载入umat.for子程序文件按下图所示创建作业（Job）
点一击段S时ub间m后it之，后结，果结如果下如图下所图示所。示。
最后查看可视化后处理，得到如下云图。
2、验证利用UMAT进行二次开发的实例
实例：最简单的杆件单轴拉伸，材料本构模型ARDMLO子A程D序；；特此殊外分，布读的取牵结
果引文力件的可U采T用RAUCRLODAFIDL ；子温程度序场。边界的 UTEMP 。
4、用户定义的单元
对于特殊类型的单元，可采用 UEL 子程序进行定义。
5、用户定义的材料特性和本构关系
几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予 ABAQUS 中的任何单元。
必须在 UMAT 中提供材料本构的雅可比（ Jacobian）矩阵，即应力增量对应变增量的变化率。
2、UMAT书写格式
定义了一些相关参数与变量的精确度，从 ABAQUS 安装目录下可找到
UMAT 中的应力矩阵、应变矩阵以及矩阵 DDSDDE、 DDSDDT、 DRPLDE 等，都是直接分量存储在前，剪切分量存储在后。直接分量有 NDI 个，剪切分量有 NSHR 个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写 UMAT 时要考虑到所使用单元的类别。

基于VUMAT的固体推进剂材料本构模型二次开发

基于VUMAT的固体推进剂材料本构模型二次开发易龙;彭云【期刊名称】《机械强度》【年(卷),期】2013(35)3【摘要】目前通用的大型商业软件都仅提供线性黏弹性本构模型,无法胜任固体推进剂材料的大变形分析。

基于Abaqus提供的用户材料子程序接口VUMAT(vectorized user defined material subroutine),选用适应性强的Swanson非线性黏弹性本构模型进行二次开发。

本构模型采用Fortran语言编写,可在Abaqus增量求解过程中作为子程序调用。

通过标准犬骨单轴拉伸算例,验证子程序的有效性。

所开发子程序考虑几何和材料双重非线性,能应用于大型固体火箭发动机药柱结构完整性分析,弥补Abaqus仅含线性黏弹性本构模型的不足。

详述材料子程序开发流程,可为用户扩充Abaqus的材料模型提供参考。

【总页数】4页(P391-394)【关键词】黏弹性;本构模型;用户子程序;固体推进剂【作者】易龙;彭云【作者单位】洪都航空工业集团第660研究所;南昌航空大学飞行器工程学院【正文语种】中文【中图分类】TB125;TB113【相关文献】1.基于能量守恒的复合固体推进剂粘弹性本构关系 [J], 龚建良;刘佩进;李强2.基于单胞体模型的复合固体推进剂松弛模量衰减过程数值模拟 [J], 何涛3.基于广义线性模型的固体推进剂贮存寿命评估 [J], 洪东跑;王英华;管飞;马小兵4.金属材料循环加载本构关系有限元模型的二次开发 [J], 杨珍;唐斌;陈涛5.基于材料基因工程的复合固体推进剂力学性能预估方法 [J], 张炜因版权原因，仅展示原文概要，查看原文内容请购买。

abaqus_二次开发详解

#开头的为注释行.第一步,建立建模环境,这一步中p y将从a b a q u s中导入建模所需的所有程序模块.f r o m p a r t i m p o r t*接下来定义草图环境m d b.m o d e l s['M o d e l A'].S k e t c h(n a m e='__p r o f i l e__',s h e e t S i z e=200.0)m d b.m o d e l s['M o d e lA'].s k e t c h e s['__p r o f i l e__'].s k e t c h O p t i o n s.s e t V a l u e s(c o n s t r u c t i o n G e o m e t r y=O N,d e c i m a l P l a c e s=2,d i m e n s i o n T e x t H e i g h t=5.0,g r i d=O N,g r i d F r e q u e n c y=2,g r i d S p a c i n g=5.0,s h e e t S i z e=200.0,v i e w S t y l e=A X I S Y M)上面的设定为大小200*200,格栅间距为5,文字标注高度为5.m d b.m o d e l s['M o d e l A'].s k e t c h e s['__p r o f i l e__'].O b l i q u e C o n s t r u c t i o n L i n e(p o i n t1=(0.0,-100.0),p o i n t2=(0.0,100.0))本句语句设定轴对称模型的对称轴线位置m d b.m o d e l s['M o d e l A'].s k e t c h e s['__p r o f i l e__'].r e c t a n g l e(p o i n t1=(0.0,0.0),p o i n t2=(40.0, -40.0))该语句绘制矩形,从点0,0至点40,-40m d b.m o d e l s['M o d e l A'].P a r t(d i m e n s i o n a l i t y=A X I S Y M M E T R I C,n a m e='B o d e n',t y p e=D E F O R M A B L E_B O D Y)定义模型为轴对称,名字为b o d e n,为可变形体m d b.m o d e l s['M o d e l A'].p a r t s['B o d e n'].B a s e S h e l l(s k e t c h=m d b.m o d e l s['M o d e lA'].s k e t c h e s['__p r o f i l e__'])d e l m d b.m o d e l s['M o d e l A'].s k e t c h e s['__p r o f i l e__']绘图完成丌要忘记收回建模环境所占的内存第二节:材料定义--------------------2楼第三节:装配--------------------3楼第四节:分析步定义--------------------4楼第五节:接触定义--------------------5楼第六节:荷载边界定义-----------------6楼第七节:网格划分控制------------------7楼第八节,任务提交及杂项功能--------8楼关于如何在p y t h o n中提交多个任务的问题9楼第二节,材料定义f r o m m a t e r i a l i m p o r t*f r o m s e c t i o n i m p o r t*从A B A Q U S提供的接口中导入材料库和组件库m d b.m o d e l s['M o d e l-A'].M a t e r i a l(n a m e='B o d e n')定义材料名m d b.m o d e l s['M o d e l A'].m a t e r i a l s['B o d e n'].D e n s i t y(t a b l e=((2000.0,),))定义材料密度m d b.m o d e l s['M o d e l A'].m a t e r i a l s['B o d e n'].E l a s t i c(t a b l e=((210546.3,0.3333),))定义材料线弹性模量和泊松比,其它的材料,如弹塑性,粘弹性材料均对应丌同的对象函数. m d b.m o d e l s['M o d e l A'].H o m o g e n e o u s S o l i d S e c t i o n(m a t e r i a l='B o d e n',n a m e='b o d e n',t h i c k n e s s=1.0)m d b.m o d e l s['M o d e lA'].p a r t s['B o d e n'].a s s i g n S e c t i o n(r e g i o n=R e g i o n(f a c e s=m d b.m o d e l s['M o d e lA'].p a r t s['B o d e n'].f a c e s[0:1]),s e c t i o n N a m e='b o d e n')设定组件为坐标无关性材料,厚度为单位厚度,并将属性附给所用的组件第三节,装配f r o m a s s e m b l y i m p o r t*首先,导入装配所用到的对象m d b.m o d e l s['M o d e lA'].r o o t A s s e m b l y.D a t u m C s y s B y T h r e e P o i n t s(c o o r d S y s T y p e=C Y L I N D R I C A L,o r i g i n=(0.0, 0.0,0.0),p o i n t1=(1.0,0.0,0.0),p o i n t2=(0.0,0.0,-1.0))定义坐标类型为柱坐标,原点0,0,0,另外两个为单位向量,确定该坐标轴的方向.m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.I n s t a n c e(n a m e='B o d e n-1',p a r t=m d b.m o d e l s['M o d e l A'].p a r t s['B o d e n'])生成草图对像b o d e n的实体,名字叨B o d e n-1.无偏移插入第四节,定义分析步f r o m s t e p i m p o r t*象其它步一样,先导入分析步要用到的模块m d b.m o d e l s['M o d e l A'].I m p l i c i t D y n a m i c s S t e p(i n i t i a l I n c=0.005,m a x N u m I n c=1024,n a m e='S t e p-1',n o S t o p=O F F,n o h a f=O F F,p r e v i o u s='I n i t i a l',t i m e I n c r e m e n t a t i o n M e t h o d=F I X E D,t i m e P e r i o d=5.12)定义对劢力隐式分析,时长为0.005*1024=5.12个时间单位,前一步为I n i t i a lm d b.m o d e l s['M o d e l A'].f i e l d O u t p u t R e q u e s t s['F-O u t p u t-1'].s e t V a l u e s(v a r i a b l e s=('U',))定义输出到O D B文件的数据,这里叧定义了位移输出m d b.m o d e l s['M o d e l A'].f i e l d O u t p u t R e q u e s t s['F-O u t p u t-1'].s e t V a l u e s(f r e q u e n c y=1)定义位移输出的频率为每步都输出m d b.m o d e l s['M o d e l A'].s t e p s['S t e p-1'].R e s t a r t(f r e q u e n c y=1,o v e r l a y=O N)定义重启劢析,每一步记录,叧记录最后一次的正确状态第五节,定义接触f r o m i n t e r a c t i o n i m p o r t*依然是先导入所用的模块m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.P a r t i t i o n E d g e B y P a r a m(e d g e s=(m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['B o d e n-1'].e d g e s[3],),p a r a m e t e r=0.975)在上部即第3面的97.5%的地方设定一个点,用于定义接触m d b.m o d e l s['M o d e l A'].C o n t a c t P r o p e r t y('I n t P r o p-1')定义接触属性名m d b.m o d e l s['M o d e l A'].i n t e r a c t i o n P r o p e r t i e s['I n t P r o p-1'].T a n g e n t i a l B e h a v i o r(f o r m u l a t i o n=F R I C T I O N L E S S)m d b.m o d e l s['M o d e l A'].i n t e r a c t i o n P r o p e r t i e s['I n t P r o p-1'].N o r m a l B e h a v i o r(a l l o w S e p a r a t i o n=O F F,a u g m e n t e d L a g r a n g e=O F F,p r e s s u r e O v e r c l o s u r e=H A R D)定义接触特性,为无摩擦硬接触丌允许分开m d b.m o d e l s['M o d e l A'].S u r f a c e T o S u r f a c e C o n t a c t S t d(a d j u s t M e t h o d=N O N E,c r e a t e S t e p N a m e='I n i t i a l',i n t e r a c t i o n P r o p e r t y='I n t P r o p-1',m a s t e r=R e g i o n(s i d e1E d g e s=m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['f u n-1'].e d g e s[0:1]) ,n a m e='I n t-1',s l a v e=R e g i o n(s i d e1E d g e s=m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['B o d e n-1'].e d g e s[4:5]),s l i d i n g=F I N I T E)这一句是建立接触对,分别为两个面上的一条边,这里边的定义由A B A Q U S内定,具体可以查阅参考手册第六节,荷载边界定义f r o m l o a d i m p o r t*m d b.m o d e l s['M o d e l A'].P e r i o d i c A m p l i t u d e(a_0=1.0,d a t a=((3.0,1.1),(3.2,1.7)),f r e q u e n c y=2.454,n a m e='F o u r i e r',s t a r t=0.0,t i m e S p a n=S T E P)定义f o u r i e r级数表示的荷载m d b.m o d e l s['M o d e l A'].P r e s s u r e(a m p l i t u d e='F o u r i e r',c r e a t e S t e p N a m e='S t e p-1',d i s t r i b u t i o n=U N I F O R M,m a g n i t u d e=50.0,n a m e='L o a d-1',re g i o n=R e g i o n(s i d e1E d g e s=m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['f u n-1'].e d g e s[2:3]))定义压强,设定加载的分析步,区域及放大系数m d b.m o d e l s['M o d e l A'].D i s p l a c e m e n t B C(a m p l i t u d e=U N S E T,c r e a t e S t e p N a m e='I n i t i a l',d i s t r i b u t i o n=U N I F O R M,l o c a l C s y s=N o n e,n a m e='B C-1',r e g i o n=R e g i o n(e d g e s=m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['B o d e n-1'].e d g e s[0:1]+\m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['B o d e n-1'].e d g e s[2:3]+\m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['f u n-1'].e d g e s[3:4]), u1=S E T,u2=U N S E T,u r3=U N S E T)m d b.m o d e l s['M o d e l A'].D i s p l a c e m e n t B C(a m p l i t u d e=U N S E T,c r e a t e S t e p N a m e='I n i t i a l',d i s t r i b u t i o n=U N I F O R M,l o c a l C s y s=N o n e,n a m e='B C-2',re g i o n=R e g i o n(e d g e s=m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['B o d e n-1'].e d g e s[1:2]),u1=U N S E T,u2=S E T,u r3=U N S E T)设定边界位移为0的边界条件,注意语法中对象相加的方法.第七节,网格划分控制f r o m m e s h i m p o r t*i m p o r t m e s h导入网格划分模块e l e m T y p e1=m e s h.E l e m T y p e(e l e m C o d e=C A X8,e l e m L i b r a r y=S T A N D A R D,s e c o n d O r d e r A c c u r a c y=O F F,h o u r g l a s s C o n t r o l=S T I F F N E S S,d i s t o r t i o n C o n t r o l=O F F)e l e m T y p e2=m e s h.E l e m T y p e(e l e m C o d e=C A X6M,e l e m L i b r a r y=S T A N D A R D)a1=m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l yf1=a1.i n s t a n c e s['B o d e n-1'].f a c e sf a c e s1=f1[0:1]r e g i o n s=(f a c e s1,)a1.s e t E l e m e n t T y p e(r e g i o n s=r e g i o n s,e l e m T y p e s=(e l e m T y p e1,e l e m T y p e2))定义其中一个物体的网格为二次8结点单元,如果其中有无法划分成四边形单元的情况,则用三角形二次6结点单元.e l e m T y p e1=m e s h.E l e m T y p e(e l e m C o d e=C A X4,e l e m L i b r a r y=S T A N D A R D)e l e m T y p e2=m e s h.E l e m T y p e(e l e m C o d e=C A X3,e l e m L i b r a r y=S T A N D A R D)a1=m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l yf1=a1.i n s t a n c e s['f u n-1'].f a c e sf a c e s1=f1[0:1]r e g i o n s=(f a c e s1,)a1.s e t E l e m e n t T y p e(r e g i o n s=r e g i o n s,e l e m T y p e s=(e l e m T y p e1,e l e m T y p e2))定义其中一个物体的网格为一次4结点单元,如果其中有无法划分成四边形单元的情况,则用三角形一次3结点单元.m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.s e e d P a r t I n s t a n c e(r e g i o n s=(m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['f u n-1'],),s i z e=0.5)m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.s e e d P a r t I n s t a n c e(r e g i o n s=(m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.i n s t a n c e s['B o d e n-1'],),s i z e=1)定义网格划分全局单元大小.m d b.m o d e l s['M o d e l A'].r o o t A s s e m b l y.g e n e r a t e M e s h(r e g i o n s=(m d b.m o d e l s['M o d e lA'].r o o t A s s e m b l y.i n s t a n c e s['B o d e n-1'],m d b.m o d e l s['M o d e lA'].r o o t A s s e m b l y.i n s t a n c e s['f u n-1']))按照定义划分单元第八节,任务提交及杂项功能m d b.m o d e l s.c h a n g e K e y(f r o m N a m e='M o d e l A',t o N a m e='F a l l-M u s t e r')修改模型名称m d b.M o d e l(n a m e='F a l l-015',o b j e c t T o C o p y=m d b.m o d e l s['F a l l-M u s t e r'])拷贝模型m d b.m o d e l s['F a l l-015'].m a t e r i a l s['B o d e n'].e l a s t i c.s e t V a l u e s(t a b l e=((210546.3,0.15),))修改模型中的材料属性m d b.J o b(c o n t a c t P r i n t=O F F,d e s c r i p t i o n='',e c h o P r i n t=O F F,e x p l i c i t P r e c i s i o n=S I N G L E,h i s t o r y P r i n t=O F F,m o d e l='F a l l-015',m o d e l P r i n t=O F F,m u l t i p r o c e s s i n g M o d e=T H R E A D S,n a m e='J o b-015',n o d a l O u t p u t P r e c i s i o n=S I N G L E, n u m C p u s=1,n u m D o m a i n s=1,p a r a l l e l i z a t i o n M e t h o d E x p l i c i t=L O O P,p a r a l l e l i z a t i o n M e t h o d S t a n d a r d=T R E E,p r e M e m o r y=1024.0,s c r a t c h='',s t a n d a r d M e m o r y=2048.0,s t a n d a r d M e m o r y P o l i c y=M O D E R A T E,t y p e=A N A L Y S I S, u s e r S u b r o u t i n e='')生成任务m d b.s a v e A s(p a t h N a m e='D:/t e m p/F u n d a m e n t S c h w i n g u n g S t u d i e')保存模型m d b.j o b s['J o b-015'].s u b m i t()提交任务关于如何在p y t h o n中提交多个任务的问题:如果您使用下面这样的命令做的提交m d b.j o b s['J o b-01'].s u b m i t()m d b.j o b s['J o b-02'].s u b m i t()m d b.j o b s['J o b-03'].s u b m i t()你就会看到,所有的任务是一次性提交的,多个任务在一起运行,这肯定丌是你想看到的结果,如何完成一个接着一个的提交呢,其实很简单,在每个任务后面加上一句m d b.j o b s['J o b-01'].w a i t F o r C o m p l e t i o n()就可以了.那么上面的语句就改为m d b.j o b s['J o b-01'].s u b m i t()m d b.j o b s['J o b-01'].w a i t F o r C o m p l e t i o n()m d b.j o b s['J o b-02'].s u b m i t()m d b.j o b s['J o b-02'].w a i t F o r C o m p l e t i o n()m d b.j o b s['J o b-03'].s u b m i t()m d b.j o b s['J o b-03'].w a i t F o r C o m p l e t i o n()......就一切O K了。

ABAQUS_Fortran二次开发

目录摘要............................................................... ABSTRACT .. (I)1.绪论 01.1.课题的研究背景 01。

2.............................................. 本文的研究内容和方法1 2。

基于ABAQUS软件的二次开发 (2)2。

1。

...................................................... ABAQUS介绍22.2。

................................................ ABAQUS各模块简介22。

3............................................. ABAQUS的二次开发平台4 2。

4。

........................................... ABAQUS的二次开发语言5 3。

用户材料子程序UMAT (7)3.1。

................................................ UMAT开发环境设置73.2.UMAT注意事项 (8)3。

3.................................................... UMAT接口的原理93.4。

................................................... UMAT的使用方法12 4。

材料非线性问题. (13)4。

1。

............................................ 材料的弹塑性本构关系144.2。

............................................. 非线性有限元算法理论174.3.增量理论常刚度法公式推导 (20)4.4.增量理论切线刚度法公式推导 (21)5.UMAT程序设计和编码 (25)5.1.本构关系描述 (25)5。

写UMAT或VUMAT概览目的书写UMAT或VUMAT所需要采取

第六讲写UMAT或VUMAT概览目的书写UMAT或VUMAT所需要采取的步骤UMAT接口例子VUMAT接口例子概览ABAQUS / Standard和ABAQUS /显有接口，使用户执行本构方程。

-在ABAQUS /标准用户定义的材料通过用户子程序UMAT模型实施。

-在ABAQUS /明确的用户定义的材料通过用户子程序VUMAT模型实施。

当在ABAQUS素材库中没有任何一种已经存在的材料，能够准确反映当前用来建模所用材料的特性时，就需要使用UMAT和VUMAT进行建模。

这些接口能够确定各种复杂本构关系的材料模型。

用户定义的材料模型可用于任何ABAQUS结构元素类型。

多用户材料可通过一个单一的UMAT或VUMAT或一起使用。

在这个讲座，在UMAT或VUMAT中的材料模型的实施将会被讨论并举例说明目的为了模拟实验结果而进行地高级的本构模行测试往往需要复杂的有限元模型。

-先进的结构element-复杂加载条件-热负荷-接触和摩擦条件-静态和动态分析如果本构模型模拟不稳定性和具有本地化现象的材料，特别的分析问题就会产生。

-准静态分析需要特别的解决方案。

-鲁棒元素配方应当提供。

-显式动态的解决方案以及强大矢量联系算法需要改进。

此外，强大的功能要求随时的可视化结果。

（就是可以动态的可视化结果）-轮廓和路径图的状态变量。

-函数地块。

-列表的结果。

材料模型的开发者应当只关注的材料模型的发展，而不是有限元软件的开发和维持。

-发展和材料没关系的建模方法-新系统的移植问题-用户开发的代码长期的计划维持书写UMAT或VUMAT所需要采取的步骤需要采取的步骤书面UMAT或VUMAT：这就需要定义一个适当的本构方程，如下：-明确定义应力（柯西应力大变形应用）-定义的应力变化率（仅在corotational框架）此外，它很可能需要：-时间，温度，或外地变量这些所依赖东西的定义-内部状态变量的定义，显式的或者用带有偏微分的函数。

umat二次开发超弹性本构

umat二次开发超弹性本构umat二次开发超弹性本构APPENDIXNeo-Hookean Hyperelatic Material User SubroutineThis program is based on the derivation of hyperelastic material constitutive model inSection 4.4. A stress and strain relationship was derived from the neo-Hookeanhyperelastic material constitutive model that is normally represented as the strain energywith strain invariants.subroutine vumat(C Read only (unmodifiable)variables -1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only (modifiable) variables -7 stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'Cdimension props(nprops), density(nblock), coordMp(nblock,*),1 charLength(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),8 defgradNew(nblock,ndir+nshr+nshr),9 fieldNew(nblock,nfieldv),1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),2 enerInternNew(nblock), enerInelasNew(nblock)Ccharacter*80 cmnameCif (cmname(1:6) .eq. 'VUMAT0') thencall VUMAT0(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)117else if (cmname(1:6) .eq. 'VUMAT1') then call VUMAT1(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)end ifendCsubroutine vumat0 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'Cdimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock),strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00,* four = 4.D00, half = 0.5D0)realC10,D1,ak,twomu,amu,alamda,hydro,vo nMises, maxShear,1 midStrain, maxPrincipalStrain118Cintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '// * 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k,1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNe w(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in localaxis)Cdet=stretchNew(k, 3)*1191 (stretchNew(k, 1)*stretchNew(k,2)-stretchNew(k, 4)**two)scale=det**(-ONE/THREE) stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B issymmetric)CBBar(k,1)=stretchNewBar(k,1)**two+stretchNewBar(k, 4)**twoBBar(k,2)=stretchNewBar(k,2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**two BBar(k,4)=stretchNewBar(k,1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3) EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/ Three) + PRstressNew(k,2)=EG*(BBar(k,2)-TRBBar/ Three) + PRstressNew(k,3)=EG*(BBar(k,3)-TRBBar/ Three) + PRstressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinates C stateNew(k,1) = stateOld(k,1) + strainInc(k,1) stateNew(k,2) = stateOld(k,2) + strainInc(k,2) stateNew(k,3) = stateOld(k,3) + strainInc(k,3) stateNew(k,4) = stateOld(k,4) +strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1201 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydro end dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2 end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C use 3/2 will get 2 (int) !CC write(6,*) totalTime,defgradNew(k,4),stretchNew(k,4)C1 ,defgradNew(k,3),defgradNew(k,4),def gradNew(k,5)C ,det,TRBBarC1 ,stressNew(k,1),stressNew(k,2),stress New(k,3),stressNew(k,4)CC Failure CriteriaCmidStrain = stateNew(k,1) + stateNew(k,2)maxShear = sqrt((stateNew(k,1) - midStrain)**2. +1 stateNew(k,4)**2.)if (midStrain .GE. 0.) then maxPrincipalStrain = midStrain + maxShearelsemaxPrincipalStrain = maxShear - midStrainend ifif (vonMises .GE. 10.8565e6) then stateNew(k,5) = 0end ifend doend ifreturnendCsubroutine vumat1 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,121* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'dimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00,* four = 4.D00, half = 0.5D0)realC10,D1,ak,twomu,amu,alamda,hydro,vo nMisesCintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '// * 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

APPENDIXNeo-Hookean Hyperelatic Material User SubroutineThis program is based on the derivation of hyperelastic material constitutive model in Section 4.4. A stress and strain relationship was derived from the neo-Hookean hyperelastic material constitutive model that is normally represented as the strain energy with strain invariants.subroutine vumat(C Read only (unmodifiable)variables -1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only (modifiable) variables -7 stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'Cdimension props(nprops), density(nblock), coordMp(nblock,*),1 charLength(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),8 defgradNew(nblock,ndir+nshr+nshr),9 fieldNew(nblock,nfieldv),1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),2 enerInternNew(nblock), enerInelasNew(nblock)Ccharacter*80 cmnameCif (cmname(1:6) .eq. 'VUMAT0') thencall VUMAT0(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)117else if (cmname(1:6) .eq. 'VUMAT1') thencall VUMAT1(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,7 stressNew, stateNew, enerInternNew, enerInelasNew)end ifendCsubroutine vumat0 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'Cdimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0)real C10,D1,ak,twomu,amu,alamda,hydro,vonMises, maxShear,1 midStrain, maxPrincipalStrain118Cintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '//* 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in local axis)Cdet=stretchNew(k, 3)*1191 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two) scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1)stateNew(k,2) = stateOld(k,2) + strainInc(k,2)stateNew(k,3) = stateOld(k,3) + strainInc(k,3)stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1201 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C use 3/2 will get 2 (int) !CC write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4) CC Failure CriteriaCmidStrain = stateNew(k,1) + stateNew(k,2)maxShear = sqrt((stateNew(k,1) - midStrain)**2. +1 stateNew(k,4)**2.)if (midStrain .GE. 0.) thenmaxPrincipalStrain = midStrain + maxShearelsemaxPrincipalStrain = maxShear - midStrainend ifif (vonMises .GE. 10.8565e6) thenstateNew(k,5) = 0end ifend doend ifreturnendCsubroutine vumat1 (C Read only -* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,* stepTime, totalTime, dt, cmname, coordMp, charLength,* props, density, strainInc, relSpinInc,* tempOld, stretchOld, defgradOld, fieldOld,121* stressOld, stateOld, enerInternOld, enerInelasOld,* tempNew, stretchNew, defgradNew, fieldNew,C Write only -* stressNew, stateNew, enerInternNew, enerInelasNew )Cinclude 'vaba_param.inc'dimension coordMp(nblock,*), charLength(nblock), props(nprops),1 density(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),9 defgradNew(nblock,ndir+nshr+nshr),1 fieldNew(nblock,nfieldv),2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),3 enerInternNew(nblock), enerInelasNew(nblock)Cdimension devia(nblock,ndir+nshr),1 BBar(nblock,4), stretchNewBar(nblock,4), intv(2)Ccharacter*80 cmnameparameter (zero = 0.D00, one = 1.D00, two = 2.D00, three = 3.D00, * four = 4.D00, half = 0.5D0)real C10,D1,ak,twomu,amu,alamda,hydro,vonMisesCintv(1) = ndirintv(2) = nshrCif (ndir .ne. 3 .or. nshr .ne. 1) thencall xplb_abqerr(1,'Subroutine VUMAT is implemented '//* 'only for plane strain and axisymmetric cases '//* '(ndir=3 and nshr=1)',0,zero,' ')call xplb_abqerr(-2,'Subroutine VUMAT has been called '//* 'with ndir=%I and nshr=%I',intv,zero,' ')call xplb_exitend ifCC10 = props(1)D1 = props(2)C C10=1.11619E6 D1=4.48E-8Cak=two/D1amu=two*C10122twomu=four*C10alamda=(three*ak-twomu)/threeCC if stepTime equals zero, assume pure elastic material and use initial elastic modulusCif(stepTime .EQ. zero) thendo k=1,nblocktrace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)stressNew(k,1) = stressOld(k,1)* + twomu * strainInc(k,1) + alamda * trace1stressNew(k,2) = stressOld(k,2)* + twomu * strainInc(k,2) + alamda * trace1stressNew(k,3) = stressOld(k,3)* + twomu * strainInc(k,3) + alamda * trace1stressNew(k,4) = stressOld(k,4)* + twomu * strainInc(k,4)C write(6,*) totalTime,k,defgradNew(k, 1),stretchNew(k,1),C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)end doelsedo k=1,nblockCC JACOBIAN OF STRETCH TENSOR (U is symmetric and in localaxis)Cdet=stretchNew(k, 3)*1 (stretchNew(k, 1)*stretchNew(k, 2)-stretchNew(k, 4)**two)scale=det**(-ONE/THREE)stretchNewBar(k, 1)=stretchNew(k, 1)*scalestretchNewBar(k, 2)=stretchNew(k, 2)*scalestretchNewBar(k, 3)=stretchNew(k, 3)*scalestretchNewBar(k, 4)=stretchNew(k, 4)*scaleCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)CBBar(k,1)=stretchNewBar(k, 1)**two+stretchNewBar(k, 4)**two BBar(k,2)=stretchNewBar(k, 2)**two+stretchNewBar(k, 4)**two BBar(k,3)=stretchNewBar(k, 3)**twoBBar(k,4)=stretchNewBar(k, 1)*stretchNewBar(k, 4) +1 stretchNewBar(k, 2)*stretchNewBar(k, 4)CC CALCULATE STRESS tensorCTRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)123EG=two*C10/detPR=two/D1*(det-one)stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR stressNew(k,4)=EG* BBar(k,4)CC Update the specific internal energyCstressPower = half * (1 ( stressOld(k,1)+stressNew(k,1) ) * strainInc(k,1) +2 ( stressOld(k,2)+stressNew(k,2) ) * strainInc(k,2) +3 ( stressOld(k,3)+stressNew(k,3) ) * strainInc(k,3) ) +4 ( stressOld(k,4)+stressNew(k,4) ) * strainInc(k,4) enerInternNew(k) = enerInternOld(k)1 + stressPower / density(k)CC Strains under corotational coordinatesCstateNew(k,1) = stateOld(k,1) + strainInc(k,1)stateNew(k,2) = stateOld(k,2) + strainInc(k,2)stateNew(k,3) = stateOld(k,3) + strainInc(k,3)stateNew(k,4) = stateOld(k,4) + strainInc(k,4)CC Calculate vonMisesChydro = (stressNew(k,1)+stressNew(k,2)+1 stressNew(k,3))/3.do k1=1,ndirdevia(k,k1) = stressNew(k,k1) - hydroend dodo k1=ndir+1,ndir+nshrdevia(k,k1) = stressNew(k,k1)end dovonMises = 0.do k1=1,ndirvonMises = vonMises + devia(k,k1)**2end dodo k1=ndir+1,ndir+nshrvonMises = vonMises + 2*devia(k,k1)**2end dovonMises = sqrt(3./2*vonMises)C write(6,*) totalTime,defgradNew(k, 4),stretchNew(k,4)C 1 ,defgradNew(k,3),defgradNew(k,4),defgradNew(k,5)C ,det,TRBBarC 1 ,stressNew(k,1),stressNew(k,2),stressNew(k,3),stressNew(k,4) 124end doend ifreturnend。