UMESHMOTION子程序模拟接触磨损

1UMESHMOTION子程序模拟接触磨损1.背景UMESHMOTION是应用于Abaqus/Standard模块的一个子程序，与ALE自适应网格技术联合使用可以控制节点的运动。由于调整节点的过程并不会改变应力应变分布，所以该子程序常用于磨损或烧蚀模拟。以下是abaqus帮助文档中给出的UMESHMOTION子程序的接口：SUBROUTINEUMESHMOTION(UREF,ULOCAL,NODE,NNDOF,*LNODETYPE,ALOCAL,NDIM,TIME,DTIME,PNEWDT,*KSTEP,KINC,KMESHSWEEP,JMATYP,JGVBLOCK,LSMOOTH)CINCLUDE'ABA_PARAM.INC'CDIMENSIONULOCAL(NDIM),JELEMLIST(*)DIMENSIONALOCAL(NDIM,*),TIME(2)DIMENSIONJMATYP(*),JGVBLOCK(*)CusercodingtodefineULOCALand,optionallyPNEWDTRETURNEND所有的参数中只有变量PNEWDT、LSMOOTH以及数组ULOCAL是需要计算的。其中ULOCAL是节点移动的速度矢量（或位移矢量，这取决于ALE中的定义），二维时ULOCAL(2)为节点面法向的分量，三维时ULOCAL(3)为节点面法向的分量。由于一般磨损或烧蚀都应该是沿着面法线方向进行的，所以我们只需要定义ULOCAL的法向分量就行了。变量PNEWDT是用于调整增量步时间的参数，一般不需要计算，变量LSMOOTH是用于标识是否应用表面平滑算法的变量，LSMOOTH=1表示应用，LSMOOTH=0表示不应用。2.模型介绍本实例中模拟一个圆盘对矩形板的磨削作用，如下图所示，圆磨盘使用解析刚体建模，矩形板由4节点线性平面应力单元组成。对刚体施加一个200N的下压力，然后刚体与矩形板开始接触。我们假设在磨削过程中磨损的速度仅与等效塑形应变有关，这里为了方便分析我们假设法向移动速度与等效塑形应变的大小成正比：v=k*PEEQ。23.建模要点1）建立磨盘和矩形板部件由于磨盘变形较小所以作为刚体处理，这里使用解析刚体建模，解析刚体支持的最大圆弧不能超过180°所以这里我们仅建立一段圆弧，半径为10mm，然后还需要在圆心处建立一个参考点作为集中力的加载点，如图1所示。矩形板处于平面应力状态，所以我们在这里建成平面壳，如图2长为40mm宽为10mm。图1刚体磨盘建模图2矩形板建模2）分配材料解析刚体不需要分配材料，所以这里仅需要建立矩形板的材料，这里使用材料的弹性模量为100GPa，泊松比为0.3，屈服强度为50MPa。在建立材料截面时应该选择实体而不是壳，这是因为我们模拟的不是壳变形而是平面应力分析。3图3材料定义3）装配如图4对部件进行装配并调整好位置即可。图4装配4）建立分析步和ALE自适应网格由于并不需要考虑动力学效应，所以此处我们建立静力学分析步，为了观察磨削的过程我们调整初始增量步和最大增量步大小如图5所示。另外ALE自适应网格也需要在step模块中设置。首先在step模块中选择other菜单然后再选择“ALEadaptiveMeshControls”建立一个ALE自适应网格的控制属性，参数设置如图6所示，然后我们再从other菜单中选择“ALEadaptiveMeshDomain”来选择应用ALE自适应网格的分析步和区域，这里的区域应该选择整个矩形块，其他相关的设置可参看图7.UMESHMOTION子程序也应该在other菜单中设置，但是这里先做一个不使用UMESHMOTION子程序的模型作为参考，后面再设置对子程序引用。4图5分析步设置图6ALE自适应网格控制图7ALE自适应网格区域选择5）建立接触在相互作用模块中首先建立硬接触属性如图8所示，需要注意的是在使用UMESHMOTION子程序时接触分析会出现非常强的非线性，所以为了使计算能够收敛我们需要将默认的约束方法修改为“罚（penalty）”并将接触刚度缩放系数设置为0.01，如图8所示。然后我们对刚体面和矩形板上面建立接触对，如图9所示。5图8接触属性定义图9定义接触对6）定义边界条件和载荷本实例中共有三处边界条件，分别位于参考点、矩形板下部和矩形板中面处，三者均为对称约束。加载点位于圆盘的参考点上，载荷类型为集中力，大小为200N，方向垂直向下，如图10所示。6图10定义载荷和边界条件7）网格划分解析刚体不需要划分网格，这里仅对矩形板划分网格，首先设置单元尺寸为0.5mm，然后选择单元类型为四节点线性平面应力单元CPS4，如图11。如图12为最终网格。图11设置单元类型图12最终的单元网格78）设置对UMESHMOTION子程序的引用这里我们首先建立一个不含UMESHMOTION子程序的计算任务作为参考并提交计算。然后我们切换到step模块，在step模块中选择other菜单然后再选择“ALEadaptiveMeshConstraint”建立一个速度型ALE自适应网格的约束属性，将区域选择为矩形板上与磨盘接触的面，然后将“motion”类型修改为“User-defined”即可完成对UMESHMOTION子程序的引用设置。然后我们再建立一个新的任务并选择要使用的子程序，如图14，最后提交任务。图13完成对UMESHMOTION子程序的引用图14设置任务使用的子程序文件84.子程序解析SUBROUTINEUMESHMOTION(UREF,ULOCAL,NODE,NNDOF,*LNODETYPE,ALOCAL,NDIM,TIME,DTIME,PNEWDT,*KSTEP,KINC,KMESHSWEEP,JMATYP,JGVBLOCK,LSMOOTH)CINCLUDE'ABA_PARAM.INC'CDIMENSIONULOCAL(NDIM)DIMENSIONALOCAL(NDIM,*),TIME(2)DIMENSIONJMATYP(*),JGVBLOCK(*)DIMENSIONARRAY(15)DIMENSIONJELEMLIST(10),JELEMTYPE(10)doubleprecisionk,PEEQC定义与节点相连的最大单元数和移动速度系数kNELEMS=10k=1.0C获取与本节点相连的单元列表JELEMLISTCALLGETNODETOELEMCONN(NODE,NELEMS,JELEMLIST,JELEMTYPE,*JRCD,JGVBLOCK)C获取节点上的各塑形分量PE（来自积分点上场输出的差值）CALLGETVRMAVGATNODE(NODE,JTYP,'PE',ARRAY,JRCD,JELEMLIST,NELEMS,*JMATYP,JGVBLOCK)C等效塑形应变为PE数组中的第7项PEEQ=ARRAY(7)C若等效塑形应变大于0则计算节点的法向运动速度ULOCAL(NDIM)if(PEEQ.GT.0.0)thenULOCAL(NDIM)=-k*PEEQendifRETURNEND5.结果分析如图15为未使用UMESHMOTION子程序时的Y方向位移结果，图16位使用UMESHMOTION子程序后的Y方向位移结果，可以看出前者的最大值0.127，而后者的最大值为0.183。这就是磨削的结果，当然这里使用的磨损公式仅是为了演示方便而设，较为简单，实际的磨损公式应该参考相关文献。9图15未使用UMESHMOTION子程序时的计算结果图16使用UMESHMOTION子程序后的计算结果6.附录6.1未使用子程序的inp文件ale-moxiao-190306.inp6.2使用子程序的inp文件umeshmotion-ale-190306.inp6.3子程序文件umeshmotion-moxiao-190306.for