+用户材料子程序实例-Johnson-Cook+金属本构模型

1、ABAQUS软件2003年度用户年会论文集中国CAE联盟 工 www.CAE 资料收集 M版权归届原作者ABAQUS/Standard用户材料子程序实例Johnson-Cook金属本构模型卢剑锋庄茁"张帆清华大学工程力学系北京100084摘要：用户材料子程序是ABAQUS提供给用户定义Fl己的材料属性的Fortran程序接口，使用 户能使用ABAQUS材料库中没有定义的材料模型。ABAQUS中H有的Johnson-Cook模熨只能应用丁显式ABAQUS/ExpUcit程序中，而我们 希與能在隐式ABAQUS/Standard程序屮更椿确的实现本构积分，而且应用Jolmson-Cook

2、模型 的修正形式。这就需更通过ABAQUS/Standard的用户材料子程序UMAT编程实现。在UMAT 编程中使用了率相关蜩性理论以及完全隐式的W力更新算法。1 Jolmson-Cook强化模型简介Jolmson-Cook (JC)模型用來模拟高应变率下的金属材料。JC强化模型农示为三项的乘积, 分别反映了应变硬化，应变率硬化和温度软化。这里使用JC模熨的修1E形式：a =（A+Bsn） 1 + Cln（l +彳（1 一广）并使参考应变率金=1,这样公式中的力即为材料的莎态屈服应力。公式中包含A,E,n,C,m五个参 数，需要通过实验來确定。2 ABAQUS用户材料子程序用户材料子程序(Us

3、eT-defined Material Mechaiucal Behavior,简称 UM AT)通过与ABAQUS ±求解程序的接口实现与ABAQUS的数据交流。在输入文件中，使用关键字'“USERMATERIAL"表示定义用户材料属性。子程序概况与接口UMAT子程序具有强大的功能，使用UMAT子程序：可以定义材料的木构关系，使用ABAQUS材料咋中没有包含的材料进行计如扩充程序功能。-1-BackABAQUS软件2003年度用户年会论文集(2) 儿乎可以用丁力学行为分析的任何分析过程，几乎可以把用户材料属性赋予ABAQUS中 的任何单元：(3) 必须在UMAT中

4、提供材料本构模型的雅可比(Jacobian)矩阵，即应力增址对应变增量 的变化率。(4) 可以和用户子程序'PSDFLD联A使用，通过rjSDFLD'亜新定义单元每一物质点上传 递到UMAT中场变最的数值。由于主程序与UMAT之间存在数据传递，其至共用一些变帚:，因此必须遵守有关UMAT的 书写格式，UMAT中常用的变応在文件开头予以泄义，通常格式为：SUBROUTINE UMAT(STRESS,ST ATEV,DDSDDESEZSPD,SCD,1 RPL,DDSDDT,DRPLDE,DRPLDTZ2 STRAN,DSTRAN,TIlvIE,DTIME,TElvIPzDTElv

5、lP,PIxEDEF,DPRED/CMN AME,3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDSQROT,PNEWDT,4 CELENT,DFGRDOZDFGRD1ZNOELZNPT,L AYER,KSPTZKSTEP,IONC)CINCLUDE 'ABA-PARAMINC*CCHARACTER*80 CMNAMEDIMENSION STRESS(NTENS)STATEV(NSTATV),1 DDSDDE(NTENS,NTENS ),DDSDDT(NTENS ),DRPLD 玖 NTENS),2 STRAN(NTENS)ZDSTCAN(NTENS)

6、ZTIME(2),PREDEF( 1 ),DPRED( 1),3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRDO(3,3),DFGRD1(3,3)user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCDand, if necessary, RPLZ DDSDDT, DRPLDE, DRPLDT, PNEWDTRETURNENDUMAT中的应力矩阵、应变矩阵以及矩阵DDSQDE, DDSDDT , DRPLDE等,都是H接分 最存储在前，剪切分帚储在后。W接分最有ND/个，剪切分最冇临饭个。各分最之间的

7、顺序 根据单元H由度的不同有一些差异，所以编写UMAT时耍考虑到所使用单元的类別。下面对 UMAT中用到的一些变虽进行说明：DDSDDE (NTENS, NTENS)是一个NTENS维的方阵，称作雅可比矩阵，沁IE, <!是应力的増鼠，山是应变的增暈, DDSDDE表示增帚:步结束时第丿个应变分帚:的改变引起的第I个应力分帚的变化。通常雅可 比是一个对称矩阵，除非在"USER MATERIAL"语句中加入:T 'TJNSUNT参数。STRESS (NTENS)应力张帚矩阵，对应NZV个氏接分帚和冋饭个剪切分量。在增帚步的开始，应力张最矩阵 屮的数值通过UMAT

8、和主程序Z间的接口传递到UMAT屮，在增量步的结束UMAT将对应力 张最矩阵更新。对丁包含刚体转动的有限应变问题，一个增駐步调用UMAT之前就已经对应力 张量的进行了刚体转动，因此在UMAT中只需处理应力张最的共旋部分。UMAT中应力张帚:的 度最为柯西(真实)应力。STATEV(NSTATEV)用丁存储状态变最的矩阵，在増暈步开始时将数值传递到UMAT中。也可在子程序USDFLD 或UEXPAN中先更新数据，然后增最步开始时将更新后的数据传递到UMAT中。在增最步的结 朿必须虫新状态变屋矩阵中的数据°和应力张最矩阵不同的是：对丁有限应变问题，除了材料木构行为引起的数据更新以外，状

9、态变量矩阵屮的任何矢戢或者张量都必须通过旋转來考虑材料的刚体运动。状态变砒矩阵的维数，等丁关键字'“DEPVAR'定义的数值。状态变量矩阵的维数通过 ABAQUS输入文件中的关键字'“DEPVAR，定义，关键字下面数据行的数值即为状态变帚矩阵的 维数。材料常数的个数，等丁关键'j"uSERMATERIAL"中"CONSTANTS'常数设定的值。PROPS(NPROPS)材料常数矩阵，矩阵中元素的数值对应丁咲键字'“USERMATERIAL下而的数据行。SSE , SPD, SCD分别定义每一增最步的弹性应变能，型性耗

10、散和疑变耗散。它们对计算结果没有影响，仅仅 3BackABAQUS软件2003年度用丿'*年会论文集作为能最输出。其他变最：STRAN NTENS):应变矩阵；DSTRAM(NTENS)： hY变增最矩阵；DTIME：増帚步的时间增帚:：NDI： X接应力分量的个数：NSHR：剪切应力分最的个数；NTENS：总应力分斌的个数，NTENS = NDI + NSHR。使用UMAT时碍要注意单元的沙漏控制刚皮和横向剪切刚皮。通常减缩积分单元的沙漏控 制刚度和板、売、梁单元的横向剪切刚度是通过材料屈性中的弹性性质定义的。这些刚皮基丁材 料初始剪切模量的值，通常在材料定义小通过"ELA

11、STIU选项定义。但是使用UMAT的时候， ABAQUS对程序输入文件进行预处理的时候得不到剪切模鼠的数值。所以这时候用户必须使用 "HOURGLASS STIFFNESS"选项來定义具有沙漏模式的单元的沙漏控制刚度，使用 TRANSVERSE SHEAR STIFFNESS"选项來定义板、壳、梁单元的横向剪切刚度。编程基丁上而所述的率相关材料公式和应力更新算法，参照ABAQUS用户材料子程序的接口规 范，进行UMAT的编程。有限元模拟结果将在下一节给出，在垠后一节中还给出了相应的程序 源代码。由于UMAT在单.元的积分点上调用，增最步开始时，主程序路径将通过UM

12、AT的接口进入 UMAT,单元当询积分点必要变量的初始值将随Z传递给UMAT的相应变量。在UMAT结束时, 变最的更新值将通过接口返回主程序。幣个UMAT的流程如图1所示。一共有8个材料常数需要给定，并申谙一个13维的状态变最矩阵，它们表不的物理含义如表 1所示。下一步将使用建立的UMAT结介ABAQUS/Standard进行袍布金森冲击杆(SHPB)实验的有 限元模拟，并对结果进行比较。图1UMAT流程图表1UMAT材料常数PROPS12345678物理性质杨氏模量泊松比塑性耗散比ABnCMSTATEV1-67-1213变最意义弹性应变塑性应变等效塑性应变3SHPB实验的有限元模拟模型的简化

13、与有限元网格为了不使模烈过于庞大，对模型进行了一些简化。首先，改变入力杆和出力杆的尺寸，长度 由原來的3040nmi减小为1 OOOmm, 径增加到25mm,试件的长度和H径也分别变化为22mm 和18nmio这样不仅优化了网格的质暈，还成倍地减小了模型的规模，其带來的负面影响就是试 件能达到的应变将降低。另外，由于撞击杆仅仅起到产生应力脉冲的作用，在数值模型中没必耍 考虑撞击杆，取代的方法是IT接在入力杆的输入端施加均布的应力脉冲。考虑到实验装置的对称性，也做了一些简化。整个实验装置以及载荷等都是关丁杆的中心线 轴对称的，所以可以使用轴对称单元进行二维分析。二维轴对称模型如图2所示。在模型中

14、，对试件以及入力杆，出力杆和试件接触的部分进行 了局部网格加密，这样的网格划分可以取得比较经济的结果。-#-BackABAQUS软件2003年度用户年会论文集入力杆出力杆试件图2二维轴对称有限元模型表2模型信息模熨尺寸inm（ xL）单元类些单元个数总节点数总单元数二维模型入力杆25x1000CAX453014751220试件18x22CAX4160出力杆25x1000CAX4530材料定义入力杆和出力杆使用线弹性材料，弹件模帚:和泊松比分别为200GPa和0.3,密度为7.85x103 kg/m3o试件采用用户在UMAT中门泄义材料，材料参数如表3所示，其屮Johnson-Cook模型 中参

15、数的数值來源丁前而的数值拟介程序。衣3试件的材料建义性质密度Kg/irf杨氏模 量 MPa泊松比Jolmson-Cook模型参数AMPaBMPanCM数2.7x10368.0x1030.3366.562108.8530.2380.0290.5边界条件为了保证SHPB实验的要求，在二维模型中施加了必要的边界条件。在对称轴上施加了对称 性边界条件，同时保证床杆和试件可以沿轴线方向IH由无约束的运动。床杆和试件之间的接触为 硬接触，光滑无摩擦。为了确定输入应力脉冲的时间，进行了简单的计算。弹性材料中纵波波速的计算公式为：q =其中E为材料弹性模量，Q为材料密度。由此可以计算输入应力波在压杆中的传播速

16、度为Cd = 5048 m/So耍求在入力杆应力波的输入端不能出现入射波和反射波的重廉，也就是说在输入应力脉冲的时间内，应力波的传播距离不应超过两倍的杆长,即:Tt < =?« 4.0x10 (s) Cd 5048根据这一佔计，选择输入应力脉冲的持续时间7； =2.0x10'4s,上升时间fr =3.0x10-5so经过若干次试算，对输入应力脉冲的波形进行适当的调格，使试件中产生较均匀的应变率。最后输入应力脉冲的波形如图3所示：图3输入应力脉冲为了饥定用战步的最大时间步长，品耍先简单计算一卜单元的稳定极限。基一个单尤的估 算，稳定极限可以用单元特征长度厶'和材料

17、波速Q定义如下：斥杆单元的特征单元长度r«10nmi.由此可以计算出应力波金斥杆传递的稳定极限为S宀0x10%)将它作为ABAQUS门动增磺控制里而的最大时间步匕 二维动态分析我们所对照的SHPB实验止是屈于这一情况，所以可以将ABAQUS/Standard结合UMAT 进行有限元模拟的结果和实验数据进行对比。下面是应变率250 s-i下的动态模拟过程。在时间f = 1.98x107($)左右，应力波前沿到达试件，这一时间和前面使用弹性波波速计算的传播时间是相同的，此前试件上的Mises应力儿乎为零，如图4所示。S, Mises|Ave Cert ： fil7e*05 66 60S

18、514e*05 3 6 心05 211etO5 060&-05 08be*04 571e*04 057e+04 543e*CI4 029rr*04 514e*04 801C-25图4应力波前沿到达试件时的Mises应力(t=1.98xl0-4 s)在时间心3.0x107(s)，试件经过应力波的上升时间后达到稳定变形的状态，一部分入射波 反射冋入力杆，一部分应力波经过试件进入出力杆，试件各点的变形都很均匀，如图5(a)所示。 在图5(b)试件的放犬图上可以看出，各点Mises应力相差不超过IMPa,这个精度是相当可靠的。(Ave. Cr it.: 75)H. 277-05*!. L70e

19、-*-051.06*et05 -*9.574e*04 _ M 一 f383H04 5 32*04M.255e*04 >3. L91e+04 2 12 8*04 *L04e+04 9. 595e 11(a)全局视图S, Mises(Ave Cir 丄t ：75 % |M.277e-tO5M.276etO5 M.275e-K)5 M.275e-tO5 r "274bOS r 1.273et-O5 k M.273e-K)5 k M.272e*O5 一 *1.272et-O5 i-1.271e-i-O5 *1.27Oe-t-O5 tl.27OetO5 M.2 69e-t-O5(b)试件

20、的放大视图图5试件经历均匀变形时的Mises应力(t=3.0xl0-4 s)经过稳定变形阶段后，反射波和传递波分别向入力杆和出力杆扩散，试件上Mises应力逐渐 减小到较低的水平，试件开始经历卸载，如图6所示。图中Mises应力云图的单位为KPa。1. 189*05 1.09D-fD5 9.9X3O047.938004 +6950d 3.962e*04 H.975e*01 3.98e*04 t3.000«t04 2.O12e*04 1.02 5e>D4 3.7n4etD2图6应力波消退后试件时的Mises应力(t=4.2xl0-4s)实际上有限元模拟的应力一应变曲线和恒定应变率

21、下实验的结果也能够很好的吻介。取出试件表面中间的一点，将应变率250 s-i和200sJ下ABAQUS有限元模拟的结果与实验的结果对比 9BackABAQUS软件2003年度用户年会论文集见图7和图8。160n140120-100-80-60-20-00.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040应变740-应变率实验值-e- ABAQUS 模拟450400350300250200«1501005000.045-,500图7应力一应变曲线的对比及模拟过程中真实应变率变化(250 sJ)图8应力一丿应变曲线的对比及模拟过程中冀

22、实应变率变化(200 S-1)WRITER 6,1)BackABAQUS软件2003年度用户年会论文集WRITER 6,1)BackABAQUS软件2003年度用户年会论文集4J-CUMAT 程序SUBROUTINE UMAT(STRESS,ST ATnEVQDSDDESESPDSCD1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRANZ2HME,DTIME；rEMPQTEMP,PREDEFQPRED,MATERL,NDI,NSHR,NTENS,3NSTATVfROPS,NPROPS,COORDSQROT,PNEWDT；CELENT,4DFGRDO,DFGRD1,NO

23、EL/NPT,KSLAY/KSPT,KSTEP/KINC)CINCLUDE 'ABA-PARAM-INC*CCHARACTER*80 MATERLDIMENSION STRESS(NTENS ),ST ATE V(NST ATV),1 DDSDD E(NTENS,NTENS)QDSDDT(NTENS ),DRPLD 玖 NTENS),2 STRAN(NTENS)ZDSTR AN(NTENS ),HME(2),PREDEF( 1 ),DPRED( 1),3 PROPS(NPROPS),COORDS(3),DROT(3,3),4 DFGRDO(3,3)QFGRD1(3,3)CDIMENSI

24、ON EEL AS( 6)ZEPL AS( 6),FLOW(6)PARAMETER (ONE=1.ODO,TWO=2.0D0zTHREE=3.0D0,SIX=6.0D0, HALF =0.5d0)DATA NEWTON,TOLER/40,1 .D-6/CCCUMAT FOR JOHNSON-COOK MODELCCPROPS( 1)-YANG'S MODULUSCPROPS(2) - POISSON RATIOCPROPS(3) - INELASTIC HEAT FRACTIONCP ARAMETERS OF JOHNSON-COOK MODEL:CPROPS(4) ACPROPS(5

25、)BCPROPS(6)nCPROPS(7) CCPROPS(8) mCCIF (NDI.NE.3) THEN1FORMAT(/,30X；FERROR - THIS UMAT MAY ONLY BE USED FOR 1'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')ENDIFCC ELASTIC PROPERTIESCEMOD=PROPS(1)ENU=PROPS(2)IF(ENU.GT.0.4999.AND.ENU.LT.0.5001)ENU=0.499EBULK3=EMOD/(ONE-TWCX ENU)EG2=EMOD/(ONE+

26、ENU)EG=EG2/TWOEG3=THREE*EGEL AM=(EBULK3-EG2)/THREECC ELASTIC STIFFNESSCDO 20 Kl=l,MIENSDO 10 K2=1ZNTENSDDSDDE(K2,K1 )=0.010 CONTINUE20 CONTINUECDO 40 K1=1,NDIDO 30 K2=1,NDIDDSDDE(K2,K1)=ELAM30 CONTINUEDDSDDE(K1,K1 )=EG2+EL AM40 CONTINUEDO 50 K1=NDI+1,NTENSWRITER 6,1)BackABAQUS软件2003年度用户年会论文集DDSDDE(K

27、1ZK1)=EG50 CONTINUECCC ALCULATE STRESS FROM EL ASTIC STRAINSCDO 70 K1=1,NTENSDO 60 K2=1,NTENSSTRESS(K2)=STRESS(K2)+DDSDDE(K2,KirDSTRAN(Kl)60 CONTINUE70 CONTINUECC RECOVER ELASTIC AND PLASTIC STRAINSCDO 80 K1=1,NTENSEELAS(K1)=STATEV(K1)+DSTRAN(K1)EPLAS(K1)=STATEV(K1+NTENS)80 CONTINUEEQPL AS=ST ATE V(

28、 1+2*NTENS)CC C ALCUL ATE MISES STRESSCIF(NPROPS.GT,5. AND.PROPS(4).GT.O.O) THENSMISES=(STRESS(1)-STRESS(2)*(STRESS(1>STRESS(2)4-1(STRESS(2STPESS(3)r(STRESS(2STRESS(3) +1(STRESS(3>STRESS( 1 )*(STRESS(3)-STRESS( 1)DO 90 K1=NDI+1,NTENSSMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)90 CONTINUESMISES=SQRT

29、(SMISES/TWO)C #BackABAQUS软件2003年度用户年会论文集CC ALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESSCCCALL USERH ARD(S YIELOZH ARD,EQPL AS,PROPS( 4)C DETERMINE IF ACTIVELY YIELDINGCIF (SMISES.GT.(1.0+TOLER)*SYIEL0) THENCCMATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTIONCSHWRO=(STRESS( 1 )+ST

30、RESS(2)+STBESS(3)/THREEONESY=ONE/SMISESDO 11OK1=1,NDIFLOW(K1 )=ONES Y*(STRESS(K1 >SHYDRO)110 CONTINUEDO 120 K1=NDI+1,NTENSFLOW(K1 )=SIPESS(K1 )*ONESY120 CONTINUECC READ PARAMETERS OF JOHNSON-COOK MODELCA=PROPS(4)B=PROPS(5)EN=PROPS(6)C=PROPS(7)EM=PROPS(8)CCNEWTON ITERATIONSYIELD=SYIEL0C 15BackABAQ

31、US软件2003年度用丿'*年会论文集DEQPL=(SNnSES-SYIELD)/EG3DSTRES=TOLER*S YIEL0/EG3DEQMIN=H ALFDTTMEDOX1. OD-4/C)DO 130 KEWTON=1,NEVVTONDEQPL=MAX(DEQPLzDEQNflN)C ALL USERHARD(SYIELD,HARD,EQPL AS+DEQPL,PROPS(4)tvpylocxdeqfl/dume)TVP1=TVP+ONEHARD1=HARD*TVP1+SYIELD*C/DEQPLSYIELD=S YIELD”TVP 1RHS=SMISES-EG3*DEQPL-

32、S1ELDDEQFL=DEQPL+RHS/(EG3+HARD1)IF(ABS(RHS/EG3) ,LE. DSTRES ) GOTO 140130 CONTINUEVVRITE(6,2)NEWTON2FORMAT(/,30X/*"WARNING - PLASTICITY ALGORITHM DID NOT1'CONVERGE AFTER ',13； ITERATIONS，)140CONTINUEEFFHRD=EG3*HARD1/(EG3+HARD1)CC CALCULATE STRESS AND UPDATE STRAINSCDO 150 K1=1,NDISTRESS(K1)=FLOW(K1)*STELD+SHYDROEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWOEELAS(Kl)=EELAS(Kl)-THREEFLOW(KirDEQPL/TVVO150 CONTINUEDO 160