abaqus_standard用户材料子程序实例_W

1、ABAQUS 软件 2003 年度用户年会论ABAQUS/Standard 用户材料子程序实例 Johnson-Cook 金属本构模型卢剑锋 庄茁* 张帆 清华大学工程力学系 北京 100084摘要：用户材料子程序是 ABAQUS 提供给用户定义自己的材料属性的 Fortran 程序接口，使用户能使用 ABAQUS 材料库中没有定义的材料模型。 ABAQUS 中自有的 Johnson-Cook 模型只能应用于显式 ABAQUS/Explicit 程序中，而我们希望能在隐式 ABAQUS/Standard 程序中更精确的实现本构积分，而且应用 Johnson-Cook 模型的修正形式。这就需要通

2、过 ABAQUS/Standard 的用户材料子程序 UMAT 编程实现。在 UMAT 编程中使用了率相关塑性理论以及完全隐式的应力更新算法。 1 Johnson-Cook 强化模型简介 Johnson-Cook（JC）模型用来模拟高应变率下的金属材料。JC 强化模型表示为三项的乘积， 分别反映了应变硬化，应变率硬化和温度软化。这里使用 JC 模型的修正形式： o = ( A + Be n ) + + e& (1 - T *m )1 C ln 1 e& 0 并使参考应变率e&0 = 1 ，这样公式中的 A 即为材料的静态屈服应力。公式中包含 A, B, n, C, m 五个参数，需要通过实验来

3、确定。 2 ABAQUS 用户材料子程序 用户材料子程序（ User-defined Material Mechanical Behavior，简称 UMAT）通过与ABAQUS 主求解程序的接口实现与 ABAQUS 的数据交流。在输入文件中，使用关键字“*USER MATERIAL”表示定义用户材料属性。 子程序概况与接口 UMAT 子程序具有强大的功能，使用 UMAT 子程序： (1) 可以定义材料的本构关系，使用 ABAQUS 材料库中没有包含的材料进行计算，扩充程序功能。 - 18 -(2) 几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予 ABAQUS 中的任何单元

4、； (3) 必须在 UMAT 中提供材料本构模型的雅可比（Jacobian）矩阵，即应力增量对应变增量的变化率。 (4) 可以和用户子程序“USDFLD”联合使用，通过“USDFLD”重新定义单元每一物质点上传递到 UMAT 中场变量的数值。 由于主程序与 UMAT 之间存在数据传递，甚至共用一些变量，因此必须遵守有关 UMAT 的书写格式，UMAT 中常用的变量在文件开头予以定义，通常格式为： SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,2 STRAN,DSTRAN,TIME,DTI

5、ME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)CINCLUDE ABA_PARAM.INCCCHARACTER*80 CMNAMEDIMENSION STRESS(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2 STRAN(NTENS),DSTRAN(NTE

6、NS),TIME(2),PREDEF(1),DPRED(1),3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDTRETURNENDUMAT 中的应力矩阵、应变矩阵以及矩阵 DDSDDE ，DDSDDT ，DRPLDE 等，都是直接分 量存储在前，剪切分量存储在后。直接分量有 NDI 个，剪切分量有 NS

7、HR 个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写 UMAT 时要考虑到所使用单元的类别。下面对UMAT 中用到的一些变量进行说明： DDSDDE ( NTENS, NTENS )是一个 NTENS 维的方阵，称作雅可比矩阵，D / D ，D 是应力的增量，D 是应变的增量， DDSDDE ( I , J ) 表示增量步结束时第 J 个应变分量的改变引起的第 I 个应力分量的变化。通常雅可比是一个对称矩阵，除非在“*USER MATERIAL”语句中加入了“UNSYMM”参数。 STRESS ( NTENS )应力张量矩阵，对应 NDI 个直接分量和 NSHR 个剪切分量。在增

8、量步的开始，应力张量矩阵中的数值通过 UMAT 和主程序之间的接口传递到 UMAT 中，在增量步的结束 UMAT 将对应力张量矩阵更新。对于包含刚体转动的有限应变问题，一个增量步调用 UMAT 之前就已经对应力张量的进行了刚体转动，因此在 UMAT 中只需处理应力张量的共旋部分。UMAT 中应力张量的度量为柯西（真实）应力。 STATEV ( NSTATEV )用于存储状态变量的矩阵，在增量步开始时将数值传递到 UMAT 中。也可在子程序 USDFLD 或 UEXPAN 中先更新数据，然后增量步开始时将更新后的数据传递到 UMAT 中。在增量步的结束必须更新状态变量矩阵中的数据。 和应力张量矩

9、阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。 状态变量矩阵的维数，等于关键字“*DEPVAR”定义的数值。状态变量矩阵的维数通过ABAQUS 输入文件中的关键字“*DEPVAR”定义，关键字下面数据行的数值即为状态变量矩阵的维数。 材料常数的个数，等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。 PROPS ( NPROPS )材料常数矩阵，矩阵中元素的数值对应于关键字“*USER MATERIAL”下面的数据行。 SSE ， SPD ， SCD分别定义每一增量步的弹性应变

10、能，塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅 作为能量输出。其他变量： STRAN ( NTENS ) ： 应 变 矩 阵 ； DSTRAN ( NTENS )：应变增量矩阵； DTIME ：增量步的时间增量； NDI ：直接应力分量的个数； NSHR ：剪切应力分量的个数； NTENS ：总应力分量的个数， NTENS = NDI + NSHR 。 使用 UMAT 时需要注意单元的沙漏控制刚度和横向剪切刚度。通常减缩积分单元的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。这些刚度基于材料初始剪切模量的值，通常在材料定义中通过“*ELASTIC”选项定义。

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

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

13、表 1 UMAT 材料常数 3 SHPB 实验的有限元模拟模型的简化与有限元网格 为了不使模型过于庞大，对模型进行了一些简化。首先，改变入力杆和出力杆的尺寸，长度由原来的 3040mm 减小为 1000mm，直径增加到 25mm，试件的长度和直径也分别变化为 22mm 和 18mm。这样不仅优化了网格的质量，还成倍地减小了模型的规模，其带来的影响就是试件能达到的应变将降低。另外，由于撞击杆仅仅起到产生应力脉冲的作用，在数值模型中没必要考虑撞击杆，取代的方法是直接在入力杆的输入端施加均布的应力脉冲。 考虑到实验装置的对称性，也做了一些简化。整个实验装置以及载荷等都是关于杆的中心线轴对称的，所以可

14、以使用轴对称单元进行二维分析。 二维轴对称模型如图 2 所示。在模型中，对试件以及入力杆，出力杆和试件接触的部分进行了局部网格加密，这样的网格划分可以取得比较经济的结果。 模 型 尺寸mm( F L ) 单元类型 单元个数 总节点 数 总单元数 二维模型 入力杆 251000CAX4530 1475 1220试件 1822CAX4160出力杆 251000CAX4530图 2二维轴对称有限元模型表 2 模型信息 材料定义 入力杆和出力杆使用线弹性材料，弹性模量和泊松比分别为 200GPa 和 0.3，密度为 7.85103 kg/m3。试件采用用户在 UMAT 中自定义材料，材料参数如表 3

15、所示，其中 Johnson-Cook 模型中参数的数值来源于前面的数值拟合程序。 表 3 试件的材料定义 杨氏模量 MPa 泊松比 Johnson-Cook 模型参数 性质 密度 Kg/m3A MPa B MPa n C M数 值 2.7103 68.0103 0.3366.562108.8530.238 0.029 0.5边界条件 为了保证 SHPB 实验的要求，在二维模型中施加了必要的边界条件。在对称轴上施加了对称性边界条件，同时保证压杆和试件可以沿轴线方向自由无约束的运动。压杆和试件之间的接触为硬接触，光滑无摩擦。 为了确定输入应力脉冲的时间，进行了简单的计算。弹性材料中纵波波速的计算公

16、式为： ErCd =其中 E 为材料弹性模量， r 为材料密度。由此可以计算输入应力波在压杆中的传播速度为 Cd = 5048 m/s。 要求在入力杆应力波的输入端不能出现入射波和反射波的重叠，也就是说在输入应力脉冲的时间内，应力波的传播距离不应超过两倍的杆长，即： sT 2L =2 4.0 10-4 (s)Cd5048sr根据这一估计，选择输入应力脉冲的持续时间T = 2.0 10-4 s，上升时间t = 3.0 10-5 s。 基 准 压 力 (1.7 10-4,200)(3.0 10-5,170)应 变 率(s-1)系数2500.902000.771000.52700.46经过若干次试算

17、，对输入应力脉冲的波形进行了适当的调整，使试件中产生较均匀的应变率。最后输入应力脉冲的波形如图 3 所示： 300250应 力 (MPa)2001501005000.00000.0001图 30.00020.00030.00040.0005时 间 (s) 输入应力脉冲 为了确定增量步的最大时间步长，需要先简单计算一下单元的稳定极限。基于一个单元的估 d算，稳定极限可以用单元特征长度 Le 和材料波速C 定义如下： CDLetstable =d压杆单元的特征单元长度 Le 10 mm，由此可以计算出应力波在压杆传递的稳定极限为 stableDt 2.010-6 (s)将它作为 ABAQUS 自动

18、增量控制里面的最大时间步长。二维动态分析 我们所对照的 SHPB 实验正是属于这一情况，所以可以将 ABAQUS/Standard 结合 UMAT 进行有限元模拟的结果和实验数据进行对比。下面是应变率 250 s-1 下的动态模拟过程。 在时间t = 1.9810-4 (s) 左右，应力波前沿到达试件，这一时间和前面使用弹性波波速计算的传播时间是相同的，此前试件上的 Mises 应力几乎为零，如图 4 所示。 图 4应力波前沿到达试件时的 Mises 应力 (t=1.9810-4 s)在时间t = 3.0 10-4 (s) ，试件经过应力波的上升时间后达到稳定变形的状态，一部分入射波 反射回入

19、力杆，一部分应力波经过试件进入出力杆，试件各点的变形都很均匀，如图 5 (a)所示。在图 5 (b)试件的放大图上可以看出，各点 Mises 应力相差不超过 1MPa，这个精度是相当可靠的。 (a) 全局视图 (b) 试件的放大视图 图 5 试件经历均匀变形时的 Mises 应力 (t=3.010-4 s)经过稳定变形阶段后，反射波和传递波分别向入力杆和出力杆扩散，试件上 Mises 应力逐渐减小到较低的水平，试件开始经历卸载，如图 6 所示。图中 Mises 应力云图的单位为 KPa。 图 6应力波消退后试件时的 Mises 应力 (t=4.210-4 s)实际上有限元模拟的应力应变曲线和恒

20、定应变率下实验的结果也能够很好的吻合。取出试件表面中间的一点，将应变率 250 s-1 和 200 s-1 下 ABAQUS 有限元模拟的结果与实验的结果对比 见图 7 和图 8。 160140120100500实验值 ABAQUS模拟应 变 率450400350应 变 率 (s-1)300应 力 (MPa)8025060200150401002050000.0000.0050.0100.0150.0200.0250.0300.0350.0400.045应 变 图 7 应力应变曲线的对比及模拟过程中真实应变率变化 (250 s-1)160140120力 (MPa)100500实验数据 ABAQ

21、US模拟 应 变 率450400350应 变 率 (s-1)30080250应60200150401002050000.0000.0050.0100.0150.0200.0250.0300.035应 变 图 8 应力应变曲线的对比及模拟过程中真实应变率变化 (200 s-1) 4 J-C UMAT 程序 SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTEN

22、S,3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)CINCLUDE ABA_PARAM.INCCCHARACTER*80 MATERLDIMENSION STRESS(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3 PROPS(NPROPS),

23、COORDS(3),DROT(3,3),4 DFGRD0(3,3),DFGRD1(3,3)CDIMENSION EELAS(6),EPLAS(6),FLOW(6)PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0, HALF =0.5d0) DATA NEWTON,TOLER/40,1.D-6/CC CUMAT FOR JOHNSON-COOK MODELC CPROPS(1) - YANGS MODULUSCPROPS(2) - POISSON RATIOCPROPS(3) - INELASTIC HEAT FRACTIONCPARAME

24、TERS OF JOHNSON-COOK MODEL: CPROPS(4) - ACPROPS(5) - BCPROPS(6) - nCPROPS(7) - CCPROPS(8) - mC CIF (NDI.NE.3) THEN WRITE(6,1)1FORMAT(/,30X,*ERROR - THIS UMAT MAY ONLY BE USED FOR ,1ELEMENTS WITH THREE DIRECT STRESS COMPONENTS) ENDIFCCELASTIC PROPERTIES CEMOD=PROPS(1) ENU=PROPS(2)IF(ENU.GT.0.4999.AND

25、.ENU.LT.0.5001) ENU=0.499 EBULK3=EMOD/(ONE-TWO*ENU) EG2=EMOD/(ONE+ENU)EG=EG2/TWO EG3=THREE*EG ELAM=(EBULK3-EG2)/THREECCELASTIC STIFFNESS CDO 20 K1=1,NTENSDO 10 K2=1,NTENS DDSDDE(K2,K1)=0.010CONTINUE20CONTINUE CDO 40 K1=1,NDIDO 30 K2=1,NDI DDSDDE(K2,K1)=ELAM30CONTINUE DDSDDE(K1,K1)=EG2+ELAM40CONTINUE

26、DO 50 K1=NDI+1,NTENSDDSDDE(K1,K1)=EG50CONTINUE CCCALCULATE STRESS FROM ELASTIC STRAINS CDO 70 K1=1,NTENSDO 60 K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)60CONTINUE70CONTINUE CCRECOVER ELASTIC AND PLASTIC STRAINS CDO 80 K1=1,NTENS EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTEN

27、S)80CONTINUEEQPLAS=STATEV(1+2*NTENS)CCCALCULATE MISES STRESS CIF(NPROPS.GT.5.AND.PROPS(4).GT.0.0) THEN SMISES=(STRESS(1)-STRESS(2)*(STRESS(1)-STRESS(2) +1(STRESS(2)-STRESS(3)*(STRESS(2)-STRESS(3) +1(STRESS(3)-STRESS(1)*(STRESS(3)-STRESS(1) DO 90 K1=NDI+1,NTENSSMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)9

28、0CONTINUE SMISES=SQRT(SMISES/TWO)CCCALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESS CCCALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4) CDETERMINE IF ACTIVELY YIELDINGCIF (SMISES.GT.(1.0+TOLER)*SYIEL0) THENCCMATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTION CSHYDRO=(STRESS(1)+STRESS(2)+STR

29、ESS(3)/THREE ONESY=ONE/SMISESDO 110 K1=1,NDI FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)110CONTINUEDO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY120CONTINUE CCREAD PARAMETERS OF JOHNSON-COOK MODEL CA=PROPS(4) B=PROPS(5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8)CCNEWTON ITERATION CSYIELD=SYIEL0DEQPL=(SMISES-SYIELD)

30、/EG3 DSTRES=TOLER*SYIEL0/EG3 DEQMIN=HALF*DTIME*EXP(1.0D-4/C) DO 130 KEWTON=1,NEWTONDEQPL=MAX(DEQPL,DEQMIN)CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4) TVP=C*LOG(DEQPL/DTIME)TVP1=TVP+ONE HARD1=HARD*TVP1+SYIELD*C/DEQPL SYIELD=SYIELD*TVP1RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD1)IF(ABS(

31、RHS/EG3) .LE. DSTRES ) GOTO 140130CONTINUE WRITE(6,2) NEWTON2FORMAT(/,30X,*WARNING - PLASTICITY ALGORITHM DID NOT ,1CONVERGE AFTER ,I3, ITERATIONS)140CONTINUE EFFHRD=EG3*HARD1/(EG3+HARD1)CCCALCULATE STRESS AND UPDATE STRAINS CDO 150 K1=1,NDI STRESS(K1)=FLOW(K1)*SYIELD+SHYDROEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO150CONTINUEDO 160 K1=NDI+1,NTENS