SPH- FEM耦合方法在切削加工数值模拟中的研究

王 睿

(安徽农业大学 工学院，合肥 230036)

SPH- FEM耦合方法在切削加工数值模拟中的研究

王 睿

(安徽农业大学 工学院，合肥 230036)

针对SPH与FEM的各自特点，采用SPH与FEM耦合的方法进行切削加工数值模拟。在变形比较小处采用FEM，以提高计算精度和效率；在变形较大处采用SPH，以避免网格畸变造成计算困难。在ABAQUS软件环境中，建立了弯颈刨刀对45号钢切槽的SPH- FEM耦合模型，模拟了刨削过程材料的应力、切削应力，验证了耦合算法的有效性。

切削；数值模拟；SPH- FEM耦合算法

切削加工在机械制造工艺中占据极其重要的地位。国内外学者对切削加工做了大量的研究，主要研究方法有解析理论法、实验法和有限元法(finite element method, FEM)等。由于切削加工过程是一个复杂的弹塑性变形和断裂的过程，解析理论法进行理论分析和推导比较困难。实验法则需要精密的实验仪器和检测设备，并且数据结果受到实验人员的操作水平的影响。有限元法是建立在弹塑性力学模型的基础上，通过变分方法，使得误差函数达到最小值并产生稳定解的数值模拟方法。它可以模拟切削中材料的变形过程，精确的显示出切削过程的瞬时应力、应变、温度等信息。通过有限元方法可以揭示切削机理，优化切削参数，对于降低生产成本，保证产品质量，提高生产效率具有重要意义。

但是依赖于网格的传统有限元法对于切削模拟并不是十分理想，原因是切削引起的大变形会导致网格严重畸变，使计算无法进行下去。目前解决的办法是采取网格重划分和删除畸变网格单元的技术。当网格畸变达到设定准则时，进行网格重划分，将原有的网格单元信息映射到新的网格单元，从而形成新的网格，对于畸变依然严重的网格单元则及时进行删除。这种方法虽然能够保证切削模拟的顺利进行，但是频繁地划分网格不仅占用了大量的计算时间，而且旧网格单元向新网格单元的映射以及删除网格单元不可避免地会引入误差，影响模拟精度，限制了有限元切削模拟的推广。

光滑粒子动力学(smoothed particle hydrodynamics, SPH)法是一种拉格朗日形式的无网格法。SPH法采用粒子划分，不依赖于网格，具有很好的自适应性，可以避免网格畸变，适合切削引起的大变形问题。然而SPH法对每个粒子进行计算时，需要搜索影响区域内近邻的粒子信息。粒子物理量计算和搜索信息都比较费时，

因此计算效率比普通的有限元法低，对于三维模型占用计算机资源较大。

针对SPH与FEM的各自特点，为提高计算效率并消除网格畸变，采用SPH与FEM耦合的方法解决切削数值模拟问题。在变形大的区域采用SPH，避免FEM的网格畸变过大造成计算困难。在变形小的区域采用FEM，以提高计算效率。SPH与FEM耦合算法分为固定耦合算法和自适应耦合算法。固定耦合算法在计算之前就已确定SPH区域和FEM区域。自适应耦合算法则在计算之前都是FEM网格，在计算过程中自动地将大变形的有限元网格单元转换为光滑粒子，并按SPH法计算物理量。本文在大型通用有限元软件ABAQUS 6.14的环境下，建立弯颈刨刀对45号钢切槽的SPH和FEM自适应耦合模型，对刨削过程的材料的应力、切削应力和切削机理进行了分析。

1 SPH法的基本原理

SPH法的核心是插值理论，它将连续介质离散成一系列具有质量的拉格朗日粒子，用积分近似表示场函数，再对域内的相邻粒子的值积分。SPH方程的构造，先进行核估计，再进行粒子近似。在SPH法中，对任意一个函数f(x)，其核估计的表达式为：

=∫Ωf(x')W(x-x',h)dx'

(1)

式中W(x-x′,h)为核函数，Ω为点x的支持域；x-x′为粒子间距离；h为SPH粒子的光滑长度。从(1)式看出函数可以通过其核函数光滑求得，核函数相当于权函数。函数f(x)导数的核估计表达式为：

<▽f(x)>=∫Ωf(x')▽W(x-x',h)dx'

(2)

图1 粒子近似示意图

从(2)式可以看出函数f(x)的导数可以转换成其核函数的导数，即函数f(x)的导数可以由其核函数导数光滑求得。

SPH法的粒子近似是将函数的解域离散成任意分布的粒子，在粒子上进行物理量计算，如图1所示。

函数及其导数的核估计的积分表达式转化为粒子的求和的离散形式的表达式为：

(3)

(4)

式中:mj为粒子j的质量，ρj为粒子j的密度；N为在粒子i影响域内所有粒子总数。从(3)和(4)式，可以看出函数及其导数可以由离散点上的函数值经过加权平均而获得，其权因子就是核函数或其导数。

核函数决定了SPH法的精度和效率，广泛使用的核函数是B样条函数，表达式为：

(5)

SPH法的计算过程是首先将连续介质离散为SPH粒子，进行近邻粒子搜索，确定每个粒子影响域中的粒子，根据动量守恒计算出每个粒子的速度及坐标。接着计算出粒子的应变率，根据质量守恒计算出粒子的密度，根据能量守恒计算出粒子的比内能。最后根据状态方程得到静水压力，根据弹塑性本构关系得到应力偏张量，依据应力张量是由静水压力和应力偏张量组成，从而得到某时刻的应力。

2 SPH-FEM的耦合处理

图2 SPH-FEM耦合算法流程图

本文采用的SPH-FEM耦合算法是自适应耦合算法，即计算之前都是FEM网格，计算中根据网格单元的变形情况决定是否转变为SPH粒子。自适应耦合算法包括SPH粒子的转换、接触滑移的计算、交界面SPH粒子的计算等内容。SPH粒子的转换可以把应变作为转换标准，当FEM网格单元应变超过设定的应变阀值时，则将其转为SPH粒子。转换过程是将超过应变阀值的网格单元删除，调整有限元的界面并计算相关的有限元节点质量，增加新的SPH粒子并计算其物理量。接触滑移的计算包括不同材料FEM单元间、FEM单元与SPH粒子间、不同材料SPH粒子间的接触滑移计算。接触滑移计算要求在接触面上的一方不能侵入另一方的内部。如果发生侵入现象，则按动量守恒、动量矩守恒及速度协调原则调节单元或粒子的速度和位置，以保障不发生侵入现象。交界面SPH粒子的计算是把交界面一侧处于粒子的影响域的FEM单元作为虚拟粒子参与SPH计算。但虚拟粒子不作为SPH粒子，而依然保留FEM特征，参与FEM计算。这样可以避免FEM一侧由于没有粒子而导致粒子分布不对称引起的计算误差。SPH-FEM耦合算法的流程图如图2所示。

3 计算模型

3.1 弯颈刨刀对45号钢切槽的SPH-FEM模型

图3 弯颈刨刀切槽的FEM模型

在初始时刻，弯颈刨刀对45号钢切槽全部用FEM建模，如图3所示。45号钢为尺寸为300mm×400mm×200mm的矩形块，在ABAQUS环境中绘制。弯颈刨刀在CREO软件绘制并导入到ABAQUS环境中。45号钢采用C3D8R网格单元划分，对切削部位采用比基体更密的网格。FEM网格转换SPH粒子的应变阀值设为0.7，即切削部位的网格应变超过0.7时，自动转换为SPH粒子。弯颈刨刀采用C3D4网格单元划分。

采用ABAQUS的EXPLICIT模块的显示算法进行切削模拟，分析时长设置为1s。定义弯颈刨刀与45号钢的接触类型为“通用接触”。弯颈刨刀沿Z轴直线运动，速度为400mm/s，约束其X、Y方向自由度及其全部转动自由度。45号钢的底部及左右两侧的全部自由度均被约束。

3.2 材料模型

由于在切削时，45号钢的塑性变形、温度和应变率均发生较大的变化，可采取Johnson- Cook本构模型作为材料模型。Johnson-Cook本构模型的表达式为：

(6)

表1 45号钢的本构参数

45号钢的损伤失效模型采取Johnson- Cook损伤模型，表达式为：

(7)

式中D为损伤参数，D∈[0、1]，初始时D=0。当D=1时材料破坏，在ABAQUS中标记相应单元的状态变量STATUS为1,在随后的切削变形中该单元即被删除。 为一个时间步的塑性应变增量， 为当前时间步的破坏应变，表达式为：

(8)

弯颈刨刀的材料为W18Cr4V，材料参数见表2。

表2 W18Cr4V参数

4 结果与分析

当弯颈刨刀与45号钢发生接触进行切削时，由图4可以看出，刀尖部位的应力最大，并且离刀尖的距离越远应力越小。由于此时刨削处于弹性变形阶段，FEM单元没有太大的畸变，因此没有SPH粒子。当切削进行到中间阶段时，由图5可以看出，此时应力最大处是紧邻刀尖处的材料，刀尖处应力次之，而已切削的材料则沿刨刀面向上运动形成切屑。这是因为在刀尖的应力作用下，刀尖前部材料达到强度极限产生裂纹，所以在裂纹尖端处材料应力最大。由于切削会产生残余内应力，因此切屑和已切好的沟槽应力虽然减少但不为零。由于切削产生较大的变形，促使一些变形大的FEM单元转变为SPH粒子分布在切屑和沟槽中，从而避免了网格畸变使计算出错，以及删除网格单元所造成的误差，更加符合切削的真实过程。切削结束时，由图6可以看出，沟槽应力均衡并较小，这是由于沟槽此时均为残余内应力。从图6还可看出在沟槽的开始处有圆角，在结束处有毛刺。这是由于开始时材料产生塑性变形，造成圆角，结束时裂纹拉长造成毛刺。从图7可以看出主切削应力在刚开始切削时，由零迅速增大，然后在一定范围内进行上下震荡，在切削结束时，又减少至零。这是因为在切削初始时，刨刀与工件接触发生弹性变形，应力迅速增大。随着切削进行，材料达到屈服极限，发生塑性变形，然后再达到强度极限，材料产生裂纹并扩展，导致刀尖不与材料直接接触，切削应力有所下降。刀具继续向前，刀尖又与材料接触，又发生弹性变形，切削应力又有所上升。如此循环，切削应力在一定范围上下震荡，直至切削完成，切削应力降为零。

图4 0.07秒时切削应力图

图5 0.3秒时切削应力图

图6 0.8秒时切削应力图

图7 主切削应力变化曲线

通过弯颈刨刀对45号钢切槽的数值模拟，证明了采用SPH- FEM耦合方法进行切削加工的仿真是可行的。SPH- FEM耦合的方法，吸收了FEM法计算效率高和SPH法模拟大变形能力强的优点，为大变形的材料切削加工提供一种高效、准确的途径。

[1] 卞梁,王肖钧,章杰. SPH/FEM 耦合算法在陶瓷复合靶抗侵彻数值模拟中的应用[J].高压物理学报,2010,24(3): 161- 167.

[2] 于峻宏. 土壤切削大变形的三维数值仿真[J]. 农业机械学报,2007,38(4):118- 121.

[3] 宿崇,唐亮,侯俊铭，等. 基于FEM与SPH耦合方法的金属切削仿真研究[J].系统仿真学报,2009，21(16):5002- 5005.

[4] 王建明,余丰,刘飞宏，等. SPH和FEM耦合法模拟磨料水射流中单磨粒加速过程[J].山东大学学报(工学版), 2011,41(5):114- 119.

[5] 郭晓光,魏延军,张小翼，等.基于FEM- SPH耦合方法的AISI4340钢超精密切削过程仿真[J].大连理工大学学报, 2013,53(4):526- 531．

[6] 林晓东,卢义玉,汤积仁，等.基于SPH- FEM 耦合算法的磨料水射流破岩数值模拟[J].振动与冲击,2014,33(18):170- 176.

[责任编辑：张永军]

Simulation Research of Cutting Based on FEM and SPH

WANG Rui

(School of Engineering, Anhui Agricultural University, Hefei 230036, China)

According to the characteristics of SPH and FEM, the process of cutting is simulated by using the coupled method of SPH and FEM. FEM is employed in the lower distortion regions in order to improve the accuracy and efficiency of calculation. SPH is employed in the higher distortion regions in order to avoid the problem of calculation for mesh distortion. In ABAQUS software environment, the coupled model of SPH and FEM for 45 steel grooving by using bent neck plane knife is established. The simulation of cutting show the stress of material and the cutting stress，and accordingly prove the effectiveness of the coupled method of SPH and FEM.

cutting; simulation; SPH- FEM coupling

2016-09-28

2017-02-20

王 睿(1973— )，男，安徽合肥人，安徽农业大学工学院讲师；研究方向：模具CAD及数值仿真。

TP391.9

A

2096-2371(2017)02-0084-05