ABAQUS_用户子程序UMAT例子.pdf-资料库

kou2080-7854119-16359647543974085497.pdf-第1页.png

第1页 / 共18页

kou2080-7854119-16359647543974085497.pdf-第2页.png

第2页 / 共18页

kou2080-7854119-16359647543974085497.pdf-第3页.png

第3页 / 共18页

kou2080-7854119-16359647543974085497.pdf-第4页.png

第4页 / 共18页

kou2080-7854119-16359647543974085497.pdf-第5页.png

第5页 / 共18页

kou2080-7854119-16359647543974085497.pdf-第6页.png

第6页 / 共18页

kou2080-7854119-16359647543974085497.pdf-第7页.png

第7页 / 共18页

kou2080-7854119-16359647543974085497.pdf-第8页.png

第8页 / 共18页

ABAQUS 软件 2003 年度用户年会论文集 ABAQUS/Standard 用户材料子程序实例－Johnson-Cook 金属本构模型卢剑锋庄茁* 张帆清华大学工程力学系北京 100084 摘要：用户材料子程序是 ABAQUS 提供给用户定义自己的材料属性的 Fortran 程序接口，使用户能使用 ABAQUS 材料库中没有定义的材料模型。 ABAQUS 中自有的 Johnson-Cook 模型只能应用于显式 ABAQUS/Explicit 程序中，而我们希望能在隐式 ABAQUS/Standard 程序中更精确的实现本构积分，而且应用 Johnson-Cook 模型的修正形式。这就需要通过 ABAQUS/Standard 的用户材料子程序 UMAT 编程实现。在 UMAT 编程中使用了率相关塑性理论以及完全隐式的应力更新算法。 1 Johnson-Cook 强化模型简介 Johnson-Cook（JC）模型用来模拟高应变率下的金属材料。JC 强化模型表示为三项的乘积，分别反映了应变硬化，应变率硬化和温度软化。这里使用 JC 模型的修正形式： σ = ( A B + n ε )  1    + C  ln 1   + ε & ε & 0        ( 1 − T m )* 并使参考应变率 0 , , 1ε =& ，这样公式中的 A 即为材料的静态屈服应力。公式中包含 , A B n C m 五个参 , 数，需要通过实验来确定。 2 ABAQUS 用户材料子程序用户材料子程序（User-defined Material Mechanical Behavior，简称 UMAT）通过与 ABAQUS 主求解程序的接口实现与 ABAQUS 的数据交流。在输入文件中，使用关键字“*USER MATERIAL”表示定义用户材料属性。子程序概况与接口 UMAT 子程序具有强大的功能，使用 UMAT 子程序： (1) 可以定义材料的本构关系，使用 ABAQUS 材料库中没有包含的材料进行计算，扩充程序功能。 - 1 -

ABAQUS 软件 2003 年度用户年会论文集 (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.INC' C 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 END UMAT 中的应力矩阵、应变矩阵以及矩阵 DDSDDE ，DDSDDT ，DRPLDE 等，都是直接分 - 2 -

ABAQUS 软件 2003 年度用户年会论文集量存储在前，剪切分量存储在后。直接分量有 NDI 个，剪切分量有 NSHR 个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写 UMAT 时要考虑到所使用单元的类别。下面对 UMAT 中用到的一些变量进行说明： ) DDSDDE NTENS NTENS ( , 是一个 NTENS 维的方阵，称作雅可比矩阵， / ∂∆ ∂∆σ ε ，∆σ 是应力的增量，∆ε 是应变的增量， DDSDDE I J 表示增量步结束时第 J 个应变分量的改变引起的第 I 个应力分量的变化。通常雅可 , ( ) 比是一个对称矩阵，除非在“*USER MATERIAL”语句中加入了“UNSYMM”参数。 ) STRESS NTENS ( 应力张量矩阵，对应 NDI 个直接分量和 NSHR 个剪切分量。在增量步的开始，应力张量矩阵中的数值通过 UMAT 和主程序之间的接口传递到 UMAT 中，在增量步的结束 UMAT 将对应力张量矩阵更新。对于包含刚体转动的有限应变问题，一个增量步调用 UMAT 之前就已经对应力张量的进行了刚体转动，因此在 UMAT 中只需处理应力张量的共旋部分。UMAT 中应力张量的度量为柯西（真实）应力。 ) STATEV NSTATEV ( 用于存储状态变量的矩阵，在增量步开始时将数值传递到 UMAT 中。也可在子程序 USDFLD 或 UEXPAN 中先更新数据，然后增量步开始时将更新后的数据传递到 UMAT 中。在增量步的结束必须更新状态变量矩阵中的数据。和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。状态变量矩阵的维数，等于关键字“*DEPVAR”定义的数值。状态变量矩阵的维数通过 ABAQUS 输入文件中的关键字“*DEPVAR”定义，关键字下面数据行的数值即为状态变量矩阵的维数。材料常数的个数，等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。 ) PROPS NPROPS ( 材料常数矩阵，矩阵中元素的数值对应于关键字“*USER MATERIAL”下面的数据行。 SSE ， SPD ， SCD 分别定义每一增量步的弹性应变能，塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅 - 3 -

作为能量输出。其他变量： ABAQUS 软件 2003 年度用户年会论文集 STRAN NTENS ：应变矩阵； ( ) DSTRAN NTENS ：应变增量矩阵； ( ) DTIME ：增量步的时间增量； NDI ：直接应力分量的个数； NSHR ：剪切应力分量的个数； NTENS ：总应力分量的个数， NTENS NDI NSHR 使用 UMAT 时需要注意单元的沙漏控制刚度和横向剪切刚度。通常减缩积分单元的沙漏控。 = + 制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。这些刚度基于材料初始剪切模量的值，通常在材料定义中通过“*ELASTIC”选项定义。但是使用 UMAT 的时候， ABAQUS 对程序输入文件进行预处理的时候得不到剪切模量的数值。所以这时候用户必须使用 “*HOURGLASS STIFFNESS” 选项来定义具有沙漏模式的单元的沙漏控制刚度，使用 “*TRANSVERSE SHEAR STIFFNESS”选项来定义板、壳、梁单元的横向剪切刚度。编程基于上面所述的率相关材料公式和应力更新算法，参照 ABAQUS 用户材料子程序的接口规范，进行 UMAT 的编程。有限元模拟结果将在下一节给出，在最后一节中还给出了相应的程序源代码。由于 UMAT 在单元的积分点上调用，增量步开始时，主程序路径将通过 UMAT 的接口进入 UMAT，单元当前积分点必要变量的初始值将随之传递给 UMAT 的相应变量。在 UMAT 结束时，变量的更新值将通过接口返回主程序。整个 UMAT 的流程如图 1 所示。一共有 8 个材料常数需要给定，并申请一个 13 维的状态变量矩阵，它们表示的物理含义如表 1 所示。下一步将使用建立的 UMAT 结合 ABAQUS/Standard 进行霍布金森冲击杆(SHPB)实验的有限元模拟，并对结果进行比较。 - 4 -

ABAQUS 软件 2003 年度用户年会论文集 7 8 C M 5 B 6 n 13 等效塑性应变图 1 UMAT 流程图表 1 UMAT 材料常数 1 2 3 PROPS 4 物理性质杨氏模量泊松比塑性耗散比 A STATEV 变量意义弹性应变塑性应变 7-12 1-6 3 SHPB 实验的有限元模拟模型的简化与有限元网格为了不使模型过于庞大，对模型进行了一些简化。首先，改变入力杆和出力杆的尺寸，长度由原来的 3040mm 减小为 1000mm，直径增加到 25mm，试件的长度和直径也分别变化为 22mm 和 18mm。这样不仅优化了网格的质量，还成倍地减小了模型的规模，其带来的负面影响就是试件能达到的应变将降低。另外，由于撞击杆仅仅起到产生应力脉冲的作用，在数值模型中没必要考虑撞击杆，取代的方法是直接在入力杆的输入端施加均布的应力脉冲。考虑到实验装置的对称性，也做了一些简化。整个实验装置以及载荷等都是关于杆的中心线轴对称的，所以可以使用轴对称单元进行二维分析。二维轴对称模型如图 2 所示。在模型中，对试件以及入力杆，出力杆和试件接触的部分进行了局部网格加密，这样的网格划分可以取得比较经济的结果。 - 5 -

ABAQUS 软件 2003 年度用户年会论文集图 2 二维轴对称有限元模型表 2 模型信息单元类型单元个数总节点数总单元数 CAX4 CAX4 CAX4 530 160 530 1475 1220 尺寸[mm] ( Φ × L ) 25×1000 18×22 25×1000 模型入力杆二维试件模型出力杆材料定义入力杆和出力杆使用线弹性材料，弹性模量和泊松比分别为 200GPa 和 0.3，密度为 7.85×103 kg/m3。试件采用用户在 UMAT 中自定义材料，材料参数如表 3 所示，其中 Johnson-Cook 模型中参数的数值来源于前面的数值拟合程序。密度 [Kg/m3 杨氏模量 ] [MPa] 泊松比 2.7×103 68.0×103 0.33 性质数值表 3 试件的材料定义 Johnson-Cook 模型参数 A [MPa ] 66.56 2 B [MPa] 108.85 3 - 6 - n C M 0.23 8 0.029 0.5

ABAQUS 软件 2003 年度用户年会论文集边界条件为了保证 SHPB 实验的要求，在二维模型中施加了必要的边界条件。在对称轴上施加了对称性边界条件，同时保证压杆和试件可以沿轴线方向自由无约束的运动。压杆和试件之间的接触为硬接触，光滑无摩擦。为了确定输入应力脉冲的时间，进行了简单的计算。弹性材料中纵波波速的计算公式为： C d = E ρ 其中 E 为材料弹性模量， ρ为材料密度。由此可以计算输入应力波在压杆中的传播速度为 dC = 5048 m/s。要求在入力杆应力波的输入端不能出现入射波和反射波的重叠，也就是说在输入应力脉冲的时间内，应力波的传播距离不应超过两倍的杆长，即： T s < L 2 C d = 2 5048 ≈ s 4.0 10 ( ) × 4 − 根据这一估计，选择输入应力脉冲的持续时间 sT = 2.0 10 × 4 − s，上升时间 rt = 3.0 10 × 5 − s。经过若干次试算，对输入应力脉冲的波形进行了适当的调整，使试件中产生较均匀的应变率。最后输入应力脉冲的波形如图 3 所示： 300 250 200 150 100 50 ) a P M ( 力应基准压力 (1.7× 10-4,200) (3.0× 10-5,170) 应变率(s-1) 系数 250 0.90 200 0.77 100 0.52 70 0.46 0 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 时间 (s) 图 3 输入应力脉冲为了确定增量步的最大时间步长，需要先简单计算一下单元的稳定极限。基于一个单元的估算，稳定极限可以用单元特征长度 eL 和材料波速 dC 定义如下： - 7 -

ABAQUS 软件 2003 年度用户年会论文集 t ∆ stable = L e C d 压杆单元的特征单元长度 10 eL ≈ mm，由此可以计算出应力波在压杆传递的稳定极限为 t ∆ stable ≈ s− 2.0 10 ( ) 6 × 将它作为 ABAQUS 自动增量控制里面的最大时间步长。二维动态分析我们所对照的 SHPB 实验正是属于这一情况，所以可以将 ABAQUS/Standard 结合 UMAT 进行有限元模拟的结果和实验数据进行对比。下面是应变率 250 s-1 下的动态模拟过程。在时间 t = s 1.98 10 ( ) × 4 − 左右，应力波前沿到达试件，这一时间和前面使用弹性波波速计算的传播时间是相同的，此前试件上的 Mises 应力几乎为零，如图 4 所示。图 4 应力波前沿到达试件时的 Mises 应力 (t=1.98×10-4 s) 在时间 t = s 3.0 10 ( ) × 4 − ，试件经过应力波的上升时间后达到稳定变形的状态，一部分入射波反射回入力杆，一部分应力波经过试件进入出力杆，试件各点的变形都很均匀，如图 5 (a)所示。在图 5 (b)试件的放大图上可以看出，各点 Mises 应力相差不超过 1MPa，这个精度是相当可靠的。 - 8 -

资料库

ABAQUS_用户子程序UMAT例子.pdf

相关推荐

开发技术

热门标签

最新资料