[交通运输]21ABAQUS用户材料子程序.doc

21 ABAQUS用户材料子程序(UMAT)虽然ABAQUS为用户提供了大量的单元库和求解模型，使用户能够利用这些模型处理绝大多数的问题；但是实际问题毕竟非常复杂，ABAQUS不可能直接求解所有可能出现的问题。所以ABAQUS提供了大量的用户自定义子程序（User Subroutine），允许用户在找不到合适模型的情况下自行定义符合自己问题的模型。这些用户子程序涵盖了建模、载荷到单元的几乎各个部分。用户子程序具有以下的功能和特点：（1）如果ABAQUS的一些固有选项模型功能有限，用户子程序可以提高ABAQUS中这些选项的功能；（2）通常用户子程序是用FORTRAN语言的代码写成；（3）它可以以几种不同的方式包含在模型中；（4）由于它们没有存储在restart文件中，如果需要的话，可以在重新开始运行时修改它；（5）在某些情况下它可以利用ABAQUS允许的已有程序。要在模型中包含用户子程序，可以利用ABAQUS执行程序，在执行程序中应用user选项指明包含这些子程序的FORTRAN源程序或者目标程序的名字。提示：ABAQUS的输入文件除了可以通过ABAQUS/CAE的作业模块中提交运行外，还可以在ABAQUS Command窗口中输入ABAQUS执行程序直接运行：ABAQUS job=输入文件名 user=用户子程序的Fortran文件名ABAQUS/Standard和ABAQUS/Explicit都支持用户子程序功能，但是他们所支持的用户子程序种类不尽相同，读者在需要使用时请注意查询手册。在接下来的两章里，我们将讨论两种常用的用户子程序用户材料子程序和用户单元子程序。本章将通过在ABAQUS/Standard中创建Johnson-Cook的材料模型，介绍编写ABAQUS/Standard的用户材料子程序UMAT。在ABAQUS/Explicit中编写用户材料子程序VUMAT与之相似，但是由于隐式和显式两种方法本身的差异，它们之间也有一些不同，请读者在具体使用前仔细查阅ABAQUS手册中的相关内容。21.1 引言用户材料子程序是ABAQUS提供给用户自定义材料属性的Fortran程序接口，它使用户能使用ABAQUS材料库中没有的材料模型。ABAQUS中自有的Johnson-Cook模型只能应用于显式ABAQUS/Explicit程序中，而我们希望能在隐式ABAQUS/Standard程序中更精确地实现本构积分，而且应用Johnson-Cook模型的修正形式。这就需要通过ABAQUS/Standard的用户材料子程序UMAT编程实现。在UMAT编程中使用了率相关塑性理论以及完全隐式的应力更新算法。21.2 模型的数学描述21.2.1 Johnson-Cook强化模型简介Johnson-Cook模型用来模拟在冲击载荷作用下的变形。Johnson-Cook强化模型（JC）表示为三项的乘积，分别反映了应变硬化、应变率硬化和温度软化。这里使用JC模型的修正形式：(21-1)模型中包含五个参数，需要通过实验来确定。使参考应变率，这样公式中的即为材料的静态屈服应力。公式中的为无量纲化的温度其中为室温，为材料的熔点。Johnson-Cook模型在温度从室温到材料熔点温度的范围内都是有效的。高应变率的变形经常伴有温升现象，这是因为材料变形过程中塑性功转化为热量。对于大多数金属，90-100的塑性变形将耗散为热量。所以JC模型中温度的变化可以用如下的公式计算：(21-2)其中，为温度的增量，为塑性耗散比，表示塑性功转化为热量的比例，为材料的比热，为材料的密度。公式（21-2）考虑的是一个绝热过程，即认为温度的升高完全起因于塑性耗散。21.2.2 率相关塑性的基本公式Johnson-Cook本构模型考虑率相关塑性，塑性变形是关联的，即塑性流动沿着屈服面的法线方向，并采用Mises屈服面。将应变的增量分解为弹性部分和塑性部分：(21-3)将上式两端同时对时间的增量微分得到率形式：(21-4)在率相关塑性中，材料的塑性反应取决于加载率，以率的形式给出材料的弹性反应：(21-5)为了发生塑性变形，率相关塑性必须满足或者超过屈服条件，塑性应变率为：，(21-6)上式即为流动法则，其中为塑性流动势能，为塑性率参数，为运动硬化时的背应力。对于各项同性硬化，不存在背应力，因此有，此时有：(21-7) (21-8)与率无关塑性不同的是，率相关塑性中等效塑性应变率不能通过一致性条件获得，而是直接通过经验定律给出，成为过应力模型： (21-9)式中是过应力，为粘性。在过应力模型中，等效塑性应变率取决于超过了多少屈服应力。上面一维的率相关塑性公式可以很方便地推广到三维情况。对于小应变的情形，应力度量之间无需区分，这里采用Cauchy应力，塑性率参数由应力和内变量的经验函数给出。对照一维情况，三维情况下分解应变率为弹性和塑性部分：(21-10)应力率和弹性应变率之间的关系为：(21-11)塑性流动法则和内变量的演化方程为：，(21-12)塑性率参数为：(21-13)对于J2流动理论，Perzyna（1971）中提出了典型的过应力模型为： (21-14)式中为Macualay括号，如果，则；如果，则。为Mises等效应力，为等效应变，为率敏感系数。对于Johnson-Cook模型，可以得到等效塑性应变率的表达式为：(21-15)其中为静态屈服应力。(21-16)21.2.3 完全隐式的应力更新算法对率形式的本构方程进行积分的算法称为应力更新算法。在完全隐式的算法中，在步骤结束时计算塑性应变和内变量的增量，同时强化屈服条件。积分算法写为： (21-17a) (21-17b)(21-17c)(21-17d)(21-17e)在时刻给出一组和应变增量，公式（21-17）是一组关于求解的非线性代数方程。将公式（21-17b）代入（21-17d）得到：(21-18)式中是弹性预测的试应力，而数值是塑性修正量，它沿着结束点塑性流动的方向。对于J2流动理论，塑性流动的方向为：(21-19)它是屈服面的法向，即。在偏应力空间，von Mises屈服面为环状，所以屈服面的法向通过圆心，如图21-1所示。图21-1 径向返回算法从图21-1可以看出，在弹性预测阶段，塑性应变和内变量保持固定；而在塑性修正阶段，总体应变保持不变。在迭代开始时，程序对应力和应变设初始值为：，(21-20)应力在第次迭代时为：(21-21)定义屈服面的单位法向矢量为：，(21-22)且在整个算法的塑性修正状态过程中始终保持不变，因此塑性应变的更新是的线性函数。在次迭代时将检查屈服条件：(21-23)若收敛，则迭代完毕，增量步结束。否则将计算塑性参数的增量： (21-24)并对塑性应变和内变量进一步更新：， (21-25a)(21-25b) (21-25c)(21-25d)(21-25e)然后将更新的变量返回屈服条件进行检查，整个过程将重复直至收敛为止。这就是增量步中应力更新的过程。21.3 ABAQUS用户材料子程序用户材料子程序（User-defined Material Mechanical Behavior，简称UMAT）通过与ABAQUS主求解程序的接口实现与ABAQUS的数据交流。在输入文件中，使用关键字“*USER MATERIAL”表示定义用户材料属性。21.3.1 子程序概况与接口UMAT子程序具有强大的功能，使用UMAT子程序：(1) 可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进行计算，扩充程序功能；(2) 几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予ABAQUS中的任何单元；(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,DTIME,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)C INCLUDE ABA_PARAM.INCC CHARACTER*80 CMNAME DIMENSION 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),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, PNEWDT RETURN ENDUMAT中的应力矩阵、应变矩阵以及矩阵，等，都是直接分量存储在前，剪切分量存储在后。直接分量有个，剪切分量有个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写UMAT时要考虑到所使用单元的类别。下面对UMAT中用到的一些变量进行说明：是一个维的方阵，称作雅可比矩阵，即，是应力的增量，是应变的增量，表示增量步结束时第个应变分量的改变引起的第个应力分量的变化。通常雅可比是一个对称矩阵，除非在“*USER MATERIAL”语句中加入了“UNSYMM”参数。应力张量矩阵，对应个直接分量和个剪切分量。在增量步的开始，应力张量矩阵中的数值通过UMAT和主程序之间的接口传递到UMAT中；在增量步的结束，UMAT将对应力张量矩阵更新。对于包含刚体转动的有限应变问题，一个增量步调用UMAT之前就已经对应力张量进行了刚体转动，因此在UMAT中只需处理应力张量的共旋部分。UMAT中应力张量的度量为柯西（真实）应力。用于存储状态变量的矩阵，在增量步开始时将数值传递到UMAT中。也可在子程序USDFLD或UEXPAN中先更新数据，然后在增量步开始时将更新后的数据传递到UMAT中。在增量步结束时必须更新状态变量矩阵中的数据。和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。状态变量矩阵的维数，等于关键字“*DEPVAR”中定义的数值。状态变量矩阵的维数通过ABAQUS输入文件中的关键字“*DEPVAR”定义，关键字下面数据行的数值即为状态变量矩阵的维数。材料常数的个数，等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。材料常数矩阵，矩阵中元素的数值对应于关键字“*USER MATERIAL”下面的数据行。，分别定义每一增量步的弹性应变能，塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅作为能量输出。其他变量：应变矩阵：应变增量矩阵：增量步的时间增量：直接应力分量的个数：剪切应力分量的个数：总应力分量的个数，使用UMAT时需要注意单元的沙漏控制刚度和横向剪切刚度。通常减缩积分单元的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。这些刚度基于材料初始剪切模量的值，通常在材料定义中通过“*ELASTIC”选项定义。但是使用UMAT时，ABAQUS对程序输入文件进行预处理的时候得不到剪切模量的数值。所以这时候用户必须使用“*HOURGLASS STIFFNESS”选项来定义具有沙漏模式的单元的沙漏控制刚度，使用“*TRANSVERSE SHEAR STIFFNESS”选项来定义板、壳、梁单元的横向剪切刚度。21.3.2 编程基于上面所述的率相关材料公式和应力更新算法，参照ABAQUS用户材料子程序的接口规范，进行UMAT的编程。有限元模拟结果将在下一节给出，在最后一节中还给出了相应的程序源代码。由于UMAT在单元的积分点上调用，增量步开始时，主程序路径将通过UMAT的接口进入UMAT，单元当前积分点必要变量的初始值将随之传递给UMAT的相应变量。在UMAT结束时，变量的更新值将通过接口返回主程序。整个UMAT的流程如图21-2所示。图21-2 UMAT流程图一共有8个材料常数需要给定，并申请了一个13维的状态变量矩阵，它们表示的物理含义如表21-1所示。表21-1 UMAT材料常数PROPS12345678物理性质杨氏模量泊松比塑性耗散比ABnCMSTATEV1-67-1213变量意义弹性应变塑性应变等效塑性应变下一步将使用建立的UMAT结合ABAQUS/Standard进行分离式Hopkinson压杆（SHPB）实验的有限元模拟，并对结果进行比较。21.4 SHPB实验的有限元模拟下面将建立SHPB实验的有限元模型，并把前面所建立的UMAT接入ABAQUS/Standard进行有限元模拟。进行有限元模拟的目的不仅是为了再现SHPB实验的过程，同时也是为了对选择的数值模型和建立的UMAT进行评价。21.4.1 分离式Hopkinson压杆（SHPB）实验分离式Hopkinson压杆（Split Hopkinson Pressure Bar，简称SHPB）实验是从经典Hopkinson实验基础之上发展而来的一种实验技术，用来测量材料的动态应力应变行为。该实验技术的理论基础是一维应力波理论，通过测量两根压杆上的应变来推导试件上的应力应变关系。分离式Hopkinson压杆实验示意图见图21-3。图21-3 分离式Hopkinson压杆装置21.4.2 有限元建模有限元模型主要是参照前面介绍的SHPB实验装置，通过载荷、边界条件等的定义在有限元中模拟SHPB实验的环境，尽量在较少的机时耗费下达到更高的精度。模型的简化与有限元网格为了不使模型过于庞大，对模型进行了一些简化。首先，改变入力杆和出力杆的尺寸，长度由原来的3040mm减小为1000mm，直径增加到25mm，试件的长度和直径也分别变化为22mm和18mm。这样不仅优化了网格的质量，还成倍地减小了模型的规模，其带来的负面影响就是试件能达到的应变将降低。另外，由于撞击杆仅仅起到产生应力脉冲的作用，在数值模型中没必要考虑撞击杆，取代的方法是直接在入力杆的输入端施加均布的应力脉冲。考虑到实验装置的对称性，也做了一些简化。整个实验装置以及载荷等都是关于杆的中心线轴对称的，所以可以使用轴对称单元进行二维分析。另外也建立了四分之一横截面的三维模型作为补充。二维轴对称模型和三维模型分别如图21-4和图21-5所示。在模型中，对试件以及入力杆，出力杆和试件接触的部分进行了局部网格加密，这样的网格划分可以取得比较经济的结果。图21-4 二维轴对称有限元模型图21-5 三维有限元模型单元类型上，选择一阶常规单元，由于没有使用减缩积分单元，所以使用UMAT时无需指定单元的沙漏控制刚度。最后的模型中，二维网格单元总数为1220，三维模型网格的单元总数为17160。关于单元的详细信息如表21-2所示。表21-2 模型信息模 型尺寸mm() 单元类型单元个数总节点数总单元数二维模型入力杆251000CAX453014751220试件1822CAX4160出力杆251000CAX4530三维模型入力杆251000C3D875002104917160试件1822C3D82160出力杆251000C3D87500在二维模型局部网格存在疏密连接的部位，一个单元边要同时和两个单元连接，这在通常的有限元网格中是不好实现的。本例在这里使用了ABAQUS的多点约束技术（Multi-Points Constrain，简称MPC）来解决，线性多点约束方式如图21-6所示。图21-6 线性多点约束技术图中点使用线性多点约束后，其节点自由度均由旁边的节点和线性插值得到。所以使用多点约束方式可以很好连接模型中网格疏密不同的部位，划分出比较精练的网格。材料定义入力杆和出力杆使用线弹性材料，弹性模量和泊松比分别为200GPa和0.3，密度为7.85103 kg/m3。试件采用用户在UMAT中的自定义材料，材料参数如表21-3所示，其中Johnson-Cook模型中参数的数值来源于对实验数据的拟合。表21-3 试件的材料定义性质密度Kg/m3杨氏模量MPa泊松比Johnson-Cook模型参数AMPaBMPanCM数值2.710368.01030.3366.562108.8530.2380.0290.5边界条件为了保证SHPB实验的要求，在二维模型和三维模型中均施加了必要的边界条件。在对称轴或对称面上施加了对称性边界条件，同时保证压杆和试件可以沿轴线方向自由无约束的运动。压杆和试件之间的接触为硬接触，光滑无摩擦。为了确定输入应力脉冲的时间，进行了简单的计算。弹性材料中纵波波速的计算公式为：(21-26)其中为材料弹性模量，为材料密度。由此可以计算输入应力波在压杆中的传播速度为m/s。要求在入力杆应力波的输入端不能出现入射波和反射波的重叠，也就是说在输入应力脉冲的时间内，应力波的传播距离不应超过两倍的杆长，即：(21-27)根据这一估计，选择输入应力脉冲的持续时间s，上升时间s。经过若干次试算，对输入应力脉冲的波形进行了适当的调整，使试件中产生较均匀的应变率。最后输入应力脉冲的波形如图21-7所示：为了确定增量步的最大时间步长，需要先简单计算一下单元的稳定极限。基于一个单元的估算，稳定极限可以用单元特征长度和材料波速定义如下：(21-28)图21-7 输入应力脉冲压杆单元的特征单元长度Le = 10 mm，由此可以计算出应力波在压杆传递的稳定极限为 (s)(21-29)将它作为ABAQUS自动增量控制里面的最大时间步长。21.4.3 二维动态分析下面将讨论二维动态分析的结果。为了便于比较，进行了三类的分析，首先是无温度影响时的强化模型，然后是考虑温度影响的强化模型，最后是考虑单元失效的模型。无温度影响强化模型分析所进行的SHPB实验正是属于这一情况，可以将ABAQUS/Standard结合UMAT进行有限元模拟结果与实验数据进行对比。下面是应变率250 s-1下的动态模拟过程。在时间左右，应力波前沿到达试件，这一时间和前面使用弹性波波速计算的传播时间是相同的，此前试件上的Mises应力几乎为零，如图21-8所示。图21-8 应力波前沿到达试件时的Mises应力 (t=1.9810-4 s)在时间，试件经过应力波的上升时间后达到稳定变形的状态，一部分入射波反射回入力杆，一部分应力波经过试件进入出力杆，试件各点的变形都很均匀，如图21-9 (a) 所示。在图21-9 (b) 试件的放大图上可以看出，各点Mises应力相差不超过1MPa，这个精度是相当可靠的。(a) 全局视图(b) 试件的放大视图图21-9 试件经历均匀变形时的Mises应力 (t=3.010-4 s)经过稳定变形阶段后，反射波和传递波分别向入力杆和出力杆扩散，试件上Mises应力逐渐减小到较低的水平，试件开始经历卸载，如图21-10所示。图中Mises应力云图的单位为kPa，如不作特别说明，下面各种应力云图中应力单位都为kPa。图21-10 应力波消退后试件时的Mises应力 (t=4.210-4 s)在离压杆两端0.2m处各取一个单元作为输出，其应力历史如图21-11所示，从中可直观地看出压杆上应力的传播过程。图21-11 压杆上的应力输出（实际输出）取出试件同一横截面的三个单元以及试件表面长度方向的三个单元，将不同点的应力应变历史比较如下：图21-12 试件横截面应变、应力及应变率历史比较图21-13 试件长度方向应变、应力及应变率历史比较从试件各点的应力应变分布上看，图21-12中应变、应力及应变率历史曲线基本重合，同一横截面内各点的变化历史基本一致。图21-13中试件长度方向各点的变形历史也基本一致，除了初始阶段和最后卸载阶段曲线出现较小波动外，其他阶段曲线也基本重合。出现小幅度波动的原因可能是应力波刚刚到达试件不能立即达到平衡状态，另外应力波经历试件的弹塑性变形后沿其他方向也会出现扩散。可以认为总体上试件的变形是均匀的，试件各点的应变率基本保持在250 s-1。其他三种应变率下（分别为70、100、200 s-1）的有限元模拟结果与250 s-1的情况大体相同，试件的变形均匀，各点的应变率基本保持在设计的水平。这给研究恒定应变率下的应力应变曲线提供了有力的前提条件。实际上有限元模拟的应力应变曲线和恒定应变率下实验的结果也能够很好的吻合。取出试件表面中间的一点，将应变率250 s-1和200 s-1下ABAQUS有限元模拟结果与实验数据对比见图21-14和图21-15。四种应变率（分别为70、100、200、250 s-1）下使用ABAQUS进行有限元模拟结果的对比见图21-16到图21-20。图21-14 应力应变曲线的对比及模拟过程中真实应变率变化 (250 s-1)图21-15 应力应变曲线的对比及模拟过程中真实应变率变化 (200 s-1)图21-16 四种应变率下的应变时间曲线图21-17 四种应变率下的应力时间曲线图21-18 四种应变率下的应力应变曲线图21-19 四种设计应变率下的真实应变率时间曲线图21-20 四种应变率下的应变率应变曲线在不同的应变率下，材料表现出不同的应力应变关系，体现出本构模型中应变率的作用。而且应力应变水平相对于应变率变化较大，在较高的应变率作用下，试件发生了较大的变形，承受的应力也较大。以上结果说明，选用Johnson-Cook强化模型能很好地反映不同应变率作用下材料的力学行为，对于金属材料的冲击问题是适用的。在此基础上建立的用户材料子程序是可靠的。考虑温度影响的强化模型分析下面将在UMAT中考虑塑性耗散向热量转化所引起的温度的改变，从而影响本构模型中的应力应变关系。这一过程是非耦合的，看作是绝热过程，即热量在单元中产生，引起单元温度的变化，但是相邻单元之间没有热量的传递。在应变率250 s1下试件最终的温度场如图21-21所示。(a) t=1.9810-4 s时的温度场(b) t=3.010-4 s时的温度场(c) t=3.010-4 s时的温度场图21-21 试件上温度场的变化取塑性耗散比为0.95，即95的塑性耗散将转化为热量。这一比例对于大多数金属是合适的。从图21-21中看出，早期应力波没有到达试件时，试件上的温度没有变化，经历塑性变形后，塑性耗散开始向热量转化，引起单元温度的变化。试件上的温升较为均匀，变形结束时试件上的最高温升约为1.9。为了使温度对屈服应力的影响明显，取Johnson-Cook模型中参数（通常对于大多数金属值在1.0左右），这相当于夸大了模型中的温度软化效应。将考虑温度效应时单元的应力应变与前面没有考虑温度效应的情况做一下比较，试件表面中间单元的温度曲线如图21-22所示。曲线中发生塑性变形后，单元的温度逐渐上升。图21-22 单元温度的变化在温度的影响下单元的应力、应变历史曲线如图21-23所示，并同前面不考虑温度作用的情况进行了对比。图中随着单元温度的升高，应力较前面不考虑温度的情况要略低，而应变则略大。这说明单元表现出温度软化的特征。应力应变曲线表示如下（图21-24），在单元进入塑性后，温度升高的同时，屈服强度和硬化率要比前面不考虑温度的情况要小。Johnson-Cook模型中最后一项关于温度软化的情况在有限元模拟中得到了验证，同时也说明UMAT能正确地考虑温度对材料本构关系的影响。(a) 应变时间曲线(b) 应力时间曲线图21-23 单元应力、应变历史图21-24 考虑温度效应与不考虑温度效应的比较21.4.4 三维动态分析本例题进行的三维动态分析也相当成功，试件在模拟过程中各处应力都比较均匀，最大应力和最小应力的差别不超过1MPa。图21-25为试件的Mises应力分布。图21-25 三维有限元模拟中试件的Mises应力分布变形过程中单元的应力、应变以及应变率历史都与二维的情况非常吻合，如图21-26所示。图21-26 三维有限元模拟的应力、应变率历史下面是应力应变曲线：图21-27 三维有限元模拟中的应力应变曲线这是对编写的UMAT用于三维实体单元的一个验证。UMAT虽然是从一维Johnson-Cook模型中建立起来，但是在UMAT中率相关塑性公式将它扩展到了三维情况，所以能用于三维实体单元也是意料之中。虽然取得的结果是一致的，但是计算时间上却有很大差别，在个人计算机上，使用轴对称二维网格完成一次计算只需要10分钟，然后完成一次三维网格的计算需要3个小时。相比而言使用轴对称二维模型要经济得多。21.5 UMAT的Fortran程序21.5.1 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,NTENS, 3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)C INCLUDE ABA_PARAM.INCC CHARACTER*80 MATERL DIMENSION 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),COORDS(3),DROT(3,3), 4 DFGRD0(3,3),DFGRD1(3,3)C DIMENSION 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 -C UMAT FOR JOHNSON-COOK MODELC -C PROPS(1) - YANGS MODULUSC PROPS(2) - POISSON RATIOC PROPS(3) - INELASTIC HEAT FRACTIONC PARAMETERS OF JOHNSON-COOK MODEL:C PROPS(4) - AC PROPS(5) - B C PROPS(6) - nC PROPS(7) - CC PROPS(8) - mC -C IF (NDI.NE.3) THEN WRITE(6,1) 1 FORMAT(/,30X,*ERROR - THIS UMAT MAY ONLY BE USED FOR , 1 ELEMENTS WITH THREE DIRECT STRESS COMPONENTS) ENDIFCC ELASTIC PROPERTIESC EMOD=PROPS(1) ENU=PROPS(2) IF(ENU.GT.0.4999.AND.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)/THREECC ELASTIC STIFFNESSC DO 20 K1=1,NTENS DO 10 K2=1,NTENS DDSDDE(K2,K1)=0.0 10 CONTINUE 20 CONTINUEC DO 40 K1=1,NDI DO 30 K2=1,NDI DDSDDE(K2,K1)=ELAM 30 CONTINUE DDSDDE(K1,K1)=EG2+ELAM 40 CONTINUE DO 50 K1=NDI+1,NTENS DDSDDE(K1,K1)=EG 50 CONTINUECC CALCULATE STRESS FROM ELASTIC STRAINSC DO 70 K1=1,NTENS DO 60 K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 60 CONTINUE 70 CONTINUECC RECOVER ELASTIC AND PLASTIC STRAINSC DO 80 K1=1,NTENS EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTENS) 80 CONTINUE EQPLAS=STATEV(1+2*NTENS)CC CALCULATE MISES STRESSC IF(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,NTENS SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) 90 CONTINUE SMISES=SQRT(SMISES/TWO)CC CALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESSCC CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4)C DETERMINE IF ACTIVELY YIELDINGC IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THENCC MATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTIONC SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3)/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) 110 CONTINUE DO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY 120 CONTINUECCREAD PARAMETERS OF JOHNSON-COOK MODELC A=PROPS(4) B=PROPS(5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8)CC NEWTON ITERATIONC SYIELD=SYIEL0 DEQPL=(SMISES-SYIELD)/EG3 DSTRES=TOLER*SYIEL0/EG3 DEQMIN=HALF*DTIME*EXP(1.0D-4/C) DO 130 KEWTON=1,NEWTON DEQPL=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*TVP1 RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD1) IF(ABS(RHS/EG3) .LE. DSTRES ) GOTO 140 130 CONTINUE WRITE(6,2) NEWTON 2 FORMAT(/,30X,*WARNING - PLASTICITY ALGORITHM DID NOT , 1 CONVERGE AFTER ,I3, ITERATIONS) 140 CONTINUE EFFHRD=EG3*HARD1/(EG3+HARD1)CC CALCULATE STRESS AND UPDATE STRAINSC DO 150 K1=1,NDI STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO 150 CONTINUE DO 160 K1=NDI+1,NTENS STRESS(K1)=FLOW(K1)*SYIELD EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL 160 CONTINUE EQPLAS=EQPLAS+DEQPL SPD=DEQPL*(SYIEL0+SYIELD)/TWO RPL = PROPS(3)*SPD/DTIMECC JACOB