基于FEM-DEM的齿轮钢随机喷丸模型及残余应力仿真研究

严宏志，伊伟彬，朱鹏飞，吴顺兴

（中南大学 机电工程学院，长沙 410000）

0 引言

喷丸强化工艺广泛应用于航空航天、汽车制造以及工程机械等行业中各种尺寸和样式的工件，尤其对承受交变载荷齿轮、螺旋桨叶、曲轴以及轴承等金属制件能大幅提高疲劳强度。

喷丸强化数值仿真是研究喷丸工艺与机理的重要途径，相关研究已从最初的二维模型发展到三维模型，从单个弹丸喷射模型发展到考虑了覆盖率的多个规则位置弹丸模型。Johnson[1]等最早建立单弹丸模型，分析了喷丸强化工艺的残余压应力。Schiffner[2]、Guagliano[3]等建立了规则的多弹丸喷丸模型，并与单弹丸模型做了对比。近年来随机弹丸模型的研究得到学者们关注：盛湘飞[4]等对铝合金进行了随机喷丸模型研究，通过再次导入残余应力的方式将长时间的仿真分段来提高效率。YangF等[5]结合规则分布弹丸与随机弹丸，建立了周期性模型。S.M.H.Gangaraj[6]对比了规则弹丸模型和随机弹丸模型，发现在引入残余应力方面结果较接近，但规则弹丸在粗糙度和覆盖率的预测上相对弱势。在受喷对象方面，李智等[7]建立了面对曲面的喷丸模型，对比了不同喷射角度对残余应力的影响。王延忠[8]研究了面向TC4钛合金板件的喷丸仿真模型，通过提取弹丸坐标的方法预测喷丸覆盖率。不难看出，随机弹丸模型和不规则表面受喷是喷丸强化工艺仿真的发展方向，但由于要达到一定覆盖率时，既需要大量弹丸，同时受喷表面的网格也需要划分的十分精细，因此前人的喷丸随机模型不可避免的会带来巨大的计算量和分析时间。

随机弹丸模型可以更准确反映喷丸机理，但达到一定覆盖率时，既需要大量弹丸，同时受喷表面网格需要更精细，给喷丸随机模型的计算带来巨大了计算量和分析时间。

应用离散元模型可以准确的表征大量随机颗粒，模拟大量弹丸对受喷工件打击的集群效应和弹丸粒子之间的相互影响，还可大大减少了计算时间。而离散元模型难以表征连续体的应力应变，有限元模型则能准确反映受喷工件的弹塑性变形以及引入的残余应力，因此建立FEM-DEM耦合模型可较好的解决随机喷丸强化仿真问题。但以往的FEM-DEM模型需要应用EDEM软件与有限元软件的联合仿真，前后处理繁琐。

本文基于ABAQUS软件加入P3D3离散单元以及离散元模型，应用python语言通过对ABAQUS的前处理inp文件重新编写，建立了FEM-DEM耦合模型来表征喷丸强化的过程，较好的满足了随机喷丸数值仿真的高效率需求。

1 仿真模型建立

喷丸强化实质上是大量多数且高速撞击在受喷表面，在短时间内快速引起弹性和塑性变形。本文所建立的FEM-DEM模型可大致分为两部分，有限元分析部分是基于ABAQUS/Explicit内核建立，主要是受喷件齿轮模型几何模型和材料本构模型；另一部分是根据PARTICLE GENERATOR颗粒流模块建立的喷射弹丸的离散元模型。

1.1 齿轮几何模型

以往的受喷模型多简化为平面板件，本文则采用了株洲齿轮厂应用于某型号车桥中的齿轮模型，该齿轮使用20CrMnTi齿轮钢，根据表中参数使用SolidWorks建立三维模型。为了后续仿真简化模型，如图1所示，截取该齿轮模型中的一个齿根及左右两齿，并导入abaqus。

图1 截取的齿轮模型

1.2 材料模型

喷丸强化高度非线性且发生在瞬时，具有大变形和高应变速率。Johnson-Cook本构模型能够准确反映材料在塑性变形过程中应变硬化、应变速率硬化效应和温度软化效应[9]。

ε为等效应变；

T为温度，℃；

Tm为材料熔化温度，Tm=1350℃；

T0为温度参考值，通常取室温T0=25℃。

王佳斌对20CrMnTi钢进行动态压缩实验，获得了不同应变率下的动态应力-应变曲线，获得了修正后的J-C本构模型[10]，故本文参考该模型来描述材料形变特性。其主要参数及相应数值如表1所示。其中A为屈服应力，单位MPa；B为幂指系数，单位MPa；n为硬化指数；C为应变率敏感性系数；m为温度敏感性系数。

表1 20CrMnTi合金钢JC本构模型系数

1.3 网格划分

网格划分的越精细结果就越准确，但同时也会导致计算量成倍增加。高应变速率的模型对网格尺寸精密程度更加敏感[11]，M.Frija建立了不同的最小单元尺寸模型，用仿真的残余应力结果对比理论计算结果，结论为最小单元的尺寸相比弹丸直径应至少小于1/10[12]。为了优化分析时间，采用单个积分点的C3D8R减缩积分单元。同时，在受喷区域表面做网格细化，使受喷区域的网格尺寸小于0.08mm。并在垂直表面方向做纵向梯度细化，以方便后处理时提取数据。划分网格后的模型如图2所示，共171360个单元。

图2 受喷齿的网格划分

1.4 边界条件与载荷

为防止在仿真过程中齿轮模型被弹丸撞击产生位移，对截取出来的齿轮模型底面做完全固定。

1.5 离散颗粒流模型

本文所建立的模型面向高强度齿轮，选用高强度CCW-φ0.8mm/G3高硬度钢丝切丸。

通过对inp文件的编写，定义弹丸与受喷表面之间的接触为hard接触，动摩擦因数μ为0.2。对于弹丸这种球形模型，根据Hertz接触理论来描述的法向接触力与法向位移的关系[13]，来定义弹丸之间的接触，最终离散元模型如图所示。其基本参数如表2所示。

表2 钢丝切丸参数

图3 随机弹丸离散模型

2 仿真模拟优化与结果分析

2.1 沙漏控制

由于减缩积分单元只有一个积分点，在某些时候单元有形变，而插值计算的应变为0进而导致内能为零，这种情况称为沙漏模式。避免使用减缩积分单元出现沙漏模式的根本方法是细化网格或是属于用多积分点单元，但这会增加计算时间，降低计算效率。人工引入“沙漏刚度”，可以一定程度上限制沙漏变形。仿真软件对沙漏刚度取值范围为0.2～3，当沙漏能与总内能的比值低于10%时，可视为该模型的沙漏模式得到比较好的控制，即该模型可用[14]。本文分别取沙漏刚度为0.2，0.5，1，2的沙漏控制进行分析。

图4所示为同一个模型下不同的沙漏刚度值产生的应变能与内能随时间的变化。当沙漏刚度为0.2，0.5，1时，沙漏能ALLAE随着仿真时间的增加增大，最大值（15.48mJ，18.82mJ，22.38mJ）与总内能ALLIE（136.59mJ）的比值大于了10%，这表明引入的沙漏刚度不足，产生了明显的沙漏现象。当沙漏刚度取到2时，沙漏能增幅度很大，但很快趋于稳定不变，这表明引入大的沙漏刚度有效控制了沙漏现象。

图4 不同沙漏刚度参数引起的沙漏能

2.2 表面覆盖率计算方法

表面覆盖率是喷丸工艺的一项基本工艺参数，它是指弹丸冲击工件表面后留下的凹坑面积占受喷面积的比率。在以往的仿真模型中，如图5所示，多采用规则的多弹丸模型，按照指定的位置冲击受喷工件，因而可以预设覆盖率。本文基于随机模型研究了一种新型表征覆盖率与弹丸数量的方法，以单弹丸模型出发，再进一步扩展到随机弹丸模型。

图5 传统喷丸仿真覆盖率表征方法

2.2.1 单弹丸模型的受喷区域界定研究

建立简化的单弹丸有限元模型，仿真后的塑性变形如图6所示。提取受喷工件沿深度方向的位移，以离弹坑的距离为横轴绘制图7（纵轴参考右轴），则距离弹坑中心第一个z向位移最大时的距离即是弹坑的半径，可得弹坑直径为3.38×10-4。

图6 单弹丸喷射有限元模型

图7 弹坑纵向变形和等效塑性变形曲线

图7中的红色曲线反映了弹丸冲击下沿深度方向的变形情况，曲线形状约为弹坑截面，弹坑面的Z向塑性应变值范围为-0.024mm～0.010mm，因此难以用该值界定临界值来区分喷与未喷。喷丸工艺中产生的弹丸位置是随机的，当要保证一定的覆盖率时，那必然某些微观表面会被多个弹丸撞击，产生的塑性变形也会累计受到弹丸冲击而产生累加的塑性变形，因此引入等效塑性变形（PEEQ）这一无量纲变量，如图所示绘制了以离弹坑的距离为变量，两段截面上的等效塑性变形曲线。不难得出，随着距离弹坑中心原来越近，PEEQ的值逐渐增大，到弹坑中心达到最大值，PEEQ的值与离心距成正比。上文得出弹坑直径3.38×10-4，由图7可知，在弹坑边缘，Z向位移最大时，PEEQ的值为0.115，因此弹坑内的PEEQ≥0.115，弹坑外PEEQ值＜0.115。由此可得，以等效塑性应变的值作为界定受喷区域时，PEEQ≥0.115的单元可视为受喷区域，反之则视为未喷区域。

2.2.2 多弹丸模型覆盖率与弹丸个数之间的关系

当弹丸直径、速度与材料等条件一定时，喷丸时间越大，喷射到工件表面的弹丸越多，覆盖率越高。如图8所示，本文取部分齿根和齿底直径6mm的圆作为总受喷区域，将该部分离散为1015个单元，再提取受喷的齿根区域1015个单元的PEEQ值，根据覆盖率的定义：凹坑面积比总受喷面积，则覆盖率Cr的值为：

图8 离散受喷区域

式中，Np为PEEQ≥0.115的单元个数；

Nt表示受喷区域内总的单元数，1015。

图9～图12分别是加载100弹丸、220弹丸、340弹丸和460弹丸的受喷面离散点的等效塑性应变值，图中红线是判别是否为受喷区域的分界线，则根据单弹丸模型的结论，统计表面1015个单元中PEEQ≥0.115的单元数目，计算得到不同弹丸冲击下的覆盖率分别为：100弹丸冲击下覆盖率为34.3%，220弹丸冲击下覆盖率为72.7%，340弹丸冲击下覆盖率为93.6%，460弹丸冲击下覆盖率为98.4%。

图9 100弹丸冲击下表面离散点PEEQ值

图10 220弹丸冲击下表面离散点PEEQ值

图11 340弹丸冲击下表面离散点PEEQ值

图12 460弹丸冲击下表面离散点PEEQ值

图13是加载弹丸数目从100增加到520所对应的覆盖率变化曲线。不难看出，随着施加的弹丸数从100增加至520，覆盖率逐渐从34.3%增加到99.4%，但整体增加幅度不断减小。当弹丸数达到460个时，覆盖率达到98.4%，由于边际递减效应，一般情况下认为覆盖率达到98%时即达到100%。因此在该模型中，弹丸数达到460个时覆盖率达到100%。

图13 不同弹丸数目对应的覆盖率

2.3 覆盖率对残余应力的影响

结合上文结论，取弹丸直径为0.8mm，弹丸速度为50m/s，入射角为90°，分别模拟了50%、100%、200%、300%、400%喷丸覆盖率下的齿根的受喷情况，得到的应力云图如图14～图17所示。而根据模拟结果，提取出沿层深方向的残余应力曲线，结果如图18所示。

图14 覆盖率100%表面残余应力云图

图15 覆盖率200%表面残余应力云图

图16 覆盖率300%表面残余应力云图

图17 覆盖率400%表面残余应力云图

图18 不同覆盖率残余应力曲线

由图14～图17可知，100%覆盖率时，弹坑还较明显，部分区域应力较大，随着覆盖率增加，受喷的齿底与齿根区域表面残余应力趋于均匀。由图18可以看到，当覆盖率从50%增加到100%时，表面残余应力由-448.22MPa增加到-845.89MPa，最大残余应力由-753.10MPa增加到-967.04MPa，同时最大残余应力层深和引入的残余应力层深也有所增加；当覆盖率增加到200%时，表面残余应力略微减少了45.34MPa，而最大残余应力略微增加了59.36MPa，层深变化不大；覆盖率继续增加，残余应力曲线介于100%和200%之间，这表明覆盖率超过一定值后，喷丸强化产生硬化现象，随着覆盖率增加，难以产生更多的塑性变形，残余应力曲线趋于稳定。而喷丸表面完整性的研究同样表明[15]，高覆盖率并不强化残余应力，而是从组织强化的角度来进一步提升疲劳强度。

4 结语

本文应用ABAQUS软件建立的FEM-DEM模型，通过对ABAQUS软件的二次开发，高效的实现了面向齿根齿底喷丸强化的模拟。解决了一般的FEM-DEM模型需要应用EDEM软件和有限元软件联合仿真、前处理和后处理繁琐的问题。

本文研究了在使用减缩积分单元时，不同的沙漏刚度对沙漏模式产生的影响，当取沙漏刚度为2时可以有效控制沙漏模式持续产生。

基于单弹丸模型研究了受喷区域的判定问题，当PEEQ≥0.115时应视为受喷区域；在此基础上基于随机模型研究了当喷嘴直径为4mm时，覆盖率与弹丸个数之间的关系，结果表明当弹丸数达到460时覆盖率可达到100%。

研究了覆盖率对残余应力场的影响，当覆盖率低于100%时，随着覆盖率增加残余应力会增大；覆盖率超过200%后持续的喷丸强化对残余应力场影响甚微。