ABAQUS-Standard用户材料子程序实例-Johnson-Cook金属本构模型卢剑锋庄茁

ABAQUS软件2003年度用户年会论文集-1-ABAQUS/Standard用户材料子程序实例�Johnson-Cook金属本构模型卢剑锋庄茁*张帆清华大学工程力学系北京100084摘要�用户材料子程序是ABAQUS提供给用户定义自己的材料属性的Fortran程序接口�使用户能使用ABAQUS材料库中没有定义的材料模型。ABAQUS中自有的Johnson-Cook模型只能应用于显式ABAQUS/Explicit程序中�而我们希望能在隐式ABAQUS/Standard程序中更精确的实现本构积分�而且应用Johnson-Cook模型的修正形式。这就需要通过ABAQUS/Standard的用户材料子程序UMAT编程实现。在UMAT编程中使用了率相关塑性理论以及完全隐式的应力更新算法。1Johnson-Cook强化模型简介Johnson-Cook�JC�模型用来模拟高应变率下的金属材料。JC强化模型表示为三项的乘积�分别反映了应变硬化�应变率硬化和温度软化。这里使用JC模型的修正形式�()()*01ln11nmABCTεσεε����=+++−����������&&并使参考应变率01ε=&�这样公式中的A即为材料的静态屈服应力。公式中包含,,,,ABnCm五个参数�需要通过实验来确定。2ABAQUS用户材料子程序用户材料子程序�User-definedMaterialMechanicalBehavior�简称UMAT�通过与ABAQUS主求解程序的接口实现与ABAQUS的数据交流。在输入文件中�使用关键字“*USERMATERIAL”表示定义用户材料属性。子程序概况与接口UMAT子程序具有强大的功能�使用UMAT子程序�(1)可以定义材料的本构关系�使用ABAQUS材料库中没有包含的材料进行计算�扩充程序功能。BackABAQUS软件2003年度用户年会论文集-2-(2)几乎可以用于力学行为分析的任何分析过程�几乎可以把用户材料属性赋予ABAQUS中的任何单元�(3)必须在UMAT中提供材料本构模型的雅可比�Jacobian�矩阵�即应力增量对应变增量的变化率。(4)可以和用户子程序“USDFLD”联合使用�通过“USDFLD”重新定义单元每一物质点上传递到UMAT中场变量的数值。由于主程序与UMAT之间存在数据传递�甚至共用一些变量�因此必须遵守有关UMAT的书写格式�UMAT中常用的变量在文件开头予以定义�通常格式为�SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1RPL,DDSDDT,DRPLDE,DRPLDT,2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)CINCLUDE'ABA_PARAM.INC'CCHARACTER*80CMNAMEDIMENSIONSTRESS(NTENS),STATEV(NSTATV),1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)usercodingtodefineDDSDDE,STRESS,STATEV,SSE,SPD,SCDand,ifnecessary,RPL,DDSDDT,DRPLDE,DRPLDT,PNEWDTRETURNENDUMAT中的应力矩阵、应变矩阵以及矩阵DDSDDE�DDSDDT�DRPLDE等�都是直接分BackABAQUS软件2003年度用户年会论文集-3-量存储在前�剪切分量存储在后。直接分量有NDI个�剪切分量有NSHR个。各分量之间的顺序根据单元自由度的不同有一些差异�所以编写UMAT时要考虑到所使用单元的类别。下面对UMAT中用到的一些变量进行说明�(),DDSDDENTENSNTENS是一个NTENS维的方阵�称作雅可比矩阵�/∂∆∂∆σε�∆σ是应力的增量�∆ε是应变的增量�(),DDSDDEIJ表示增量步结束时第J个应变分量的改变引起的第I个应力分量的变化。通常雅可比是一个对称矩阵�除非在“*USERMATERIAL”语句中加入了“UNSYMM”参数。()STRESSNTENS应力张量矩阵�对应NDI个直接分量和NSHR个剪切分量。在增量步的开始�应力张量矩阵中的数值通过UMAT和主程序之间的接口传递到UMAT中�在增量步的结束UMAT将对应力张量矩阵更新。对于包含刚体转动的有限应变问题�一个增量步调用UMAT之前就已经对应力张量的进行了刚体转动�因此在UMAT中只需处理应力张量的共旋部分。UMAT中应力张量的度量为柯西�真实�应力。()STATEVNSTATEV用于存储状态变量的矩阵�在增量步开始时将数值传递到UMAT中。也可在子程序USDFLD或UEXPAN中先更新数据�然后增量步开始时将更新后的数据传递到UMAT中。在增量步的结束必须更新状态变量矩阵中的数据。和应力张量矩阵不同的是�对于有限应变问题�除了材料本构行为引起的数据更新以外�状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。状态变量矩阵的维数�等于关键字“*DEPVAR”定义的数值。状态变量矩阵的维数通过ABAQUS输入文件中的关键字“*DEPVAR”定义�关键字下面数据行的数值即为状态变量矩阵的维数。材料常数的个数�等于关键字“*USERMATERIAL”中“CONSTANTS”常数设定的值。()PROPSNPROPS材料常数矩阵�矩阵中元素的数值对应于关键字“*USERMATERIAL”下面的数据行。SSE�SPD�SCD分别定义每一增量步的弹性应变能�塑性耗散和蠕变耗散。它们对计算结果没有影响�仅仅BackABAQUS软件2003年度用户年会论文集-4-作为能量输出。其他变量�()STRANNTENS�应变矩阵�()DSTRANNTENS�应变增量矩阵�DTIME�增量步的时间增量�NDI�直接应力分量的个数�NSHR�剪切应力分量的个数�NTENS�总应力分量的个数�NTENSNDINSHR=+。使用UMAT时需要注意单元的沙漏控制刚度和横向剪切刚度。通常减缩积分单元的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。这些刚度基于材料初始剪切模量的值�通常在材料定义中通过“*ELASTIC”选项定义。但是使用UMAT的时候�ABAQUS对程序输入文件进行预处理的时候得不到剪切模量的数值。所以这时候用户必须使用“*HOURGLASSSTIFFNESS”选项来定义具有沙漏模式的单元的沙漏控制刚度�使用“*TRANSVERSESHEARSTIFFNESS”选项来定义板、壳、梁单元的横向剪切刚度。编程基于上面所述的率相关材料公式和应力更新算法�参照ABAQUS用户材料子程序的接口规范�进行UMAT的编程。有限元模拟结果将在下一节给出�在最后一节中还给出了相应的程序源代码。由于UMAT在单元的积分点上调用�增量步开始时�主程序路径将通过UMAT的接口进入UMAT�单元当前积分点必要变量的初始值将随之传递给UMAT的相应变量。在UMAT结束时�变量的更新值将通过接口返回主程序。整个UMAT的流程如图1所示。一共有8个材料常数需要给定�并申请一个13维的状态变量矩阵�它们表示的物理含义如表1所示。下一步将使用建立的UMAT结合ABAQUS/Standard进行霍布金森冲击杆(SHPB)实验的有限元模拟�并对结果进行比较。BackABAQUS软件2003年度用户年会论文集-5-图1UMAT流程图表1UMAT材料常数PROPS12345678物理性质杨氏模量泊松比塑性耗散比ABnCMSTATEV1-67-1213变量意义弹性应变塑性应变等效塑性应变3SHPB实验的有限元模拟模型的简化与有限元网格为了不使模型过于庞大�对模型进行了一些简化。首先�改变入力杆和出力杆的尺寸�长度由原来的3040mm减小为1000mm�直径增加到25mm�试件的长度和直径也分别变化为22mm和18mm。这样不仅优化了网格的质量�还成倍地减小了模型的规模�其带来的负面影响就是试件能达到的应变将降低。另外�由于撞击杆仅仅起到产生应力脉冲的作用�在数值模型中没必要考虑撞击杆�取代的方法是直接在入力杆的输入端施加均布的应力脉冲。考虑到实验装置的对称性�也做了一些简化。整个实验装置以及载荷等都是关于杆的中心线轴对称的�所以可以使用轴对称单元进行二维分析。二维轴对称模型如图2所示。在模型中�对试件以及入力杆�出力杆和试件接触的部分进行了局部网格加密�这样的网格划分可以取得比较经济的结果。BackABAQUS软件2003年度用户年会论文集-6-图2二维轴对称有限元模型表2模型信息模型尺寸[mm](Φ×L)单元类型单元个数总节点数总单元数入力杆25×1000CAX4530试件18×22CAX4160二维模型出力杆25×1000CAX453014751220材料定义入力杆和出力杆使用线弹性材料�弹性模量和泊松比分别为200GPa和0.3�密度为7.85×103kg/m3。试件采用用户在UMAT中自定义材料�材料参数如表3所示�其中Johnson-Cook模型中参数的数值来源于前面的数值拟合程序。表3试件的材料定义Johnson-Cook模型参数性质密度[Kg/m3]杨氏模量[MPa]泊松比A[MPa]B[MPa]nCM数值2.7×10368.0×1030.3366.562108.8530.2380.0290.5BackABAQUS软件2003年度用户年会论文集-7-边界条件为了保证SHPB实验的要求�在二维模型中施加了必要的边界条件。在对称轴上施加了对称性边界条件�同时保证压杆和试件可以沿轴线方向自由无约束的运动。压杆和试件之间的接触为硬接触�光滑无摩擦。为了确定输入应力脉冲的时间�进行了简单的计算。弹性材料中纵波波速的计算公式为�dECρ=其中E为材料弹性模量�ρ为材料密度。由此可以计算输入应力波在压杆中的传播速度为5048dC=m/s。要求在入力杆应力波的输入端不能出现入射波和反射波的重叠�也就是说在输入应力脉冲的时间内�应力波的传播距离不应超过两倍的杆长�即�4224.010()5048sdLTsC−=≈×根据这一估计�选择输入应力脉冲的持续时间42.010sT−=×s�上升时间53.010rt−=×s。经过若干次试算�对输入应力脉冲的波形进行了适当的调整�使试件中产生较均匀的应变率。最后输入应力脉冲的波形如图3所示�0.00000.00010.00020.00030.00040.0005050100150200250300应力(MPa)时间(s)基准压力(3.0×10-5,170)(1.7×10-4,200)应变率(s-1)系数2500.902000.771000.52700.46图3输入应力脉冲为了确定增量步的最大时间步长�需要先简单计算一下单元的稳定极限。基于一个单元的估算�稳定极限可以用单元特征长度eL和材料波速dC定义如下�BackABAQUS软件2003年度用户年会论文集-8-establedLtC∆=压杆单元的特征单元长度10eL≈mm�由此可以计算出应力波在压杆传递的稳定极限为62.010()stablets−∆≈×将它作为ABAQUS自动增量控制里面的最大时间步长。二维动态分析我们所对照的SHPB实验正是属于这一情况�所以可以将ABAQUS/Standard结合UMAT进行有限元模拟的结果和实验数据进行对比。下面是应变率250s-1下的动态模拟过程。在时间