CN110597180A - 一种基于切削过程仿真的单晶铜微铣削力预测方法 - Google Patents

Abstract

本发明属于微小零件精密高效加工领域，特别涉及一种基于切削过程仿真的单晶铜微铣削力预测方法。由于单晶铜材料的特殊性，本发明首先综合考虑了晶体变形运动学、晶体塑性本构关系，率相关晶体的硬化规律等因素对于微铣削单晶铜的影响；然后对本构方程的增量形式进行了推导，采用Fortran语言编写了单晶铜VUMAT用户材料子程序，将晶体塑性本构引入到有限元仿真中；最后采用ABAQUS软件对单晶铜微铣削过程进行仿真建模，实现了单晶铜微铣削力预测。本发明以有限元切削仿真为基础，实现单晶铜微铣削力的准确预测，为研究单晶铜微铣削过程提供技术支撑，提高单晶铜微铣削加工精度与效率，具有实际应用价值。

Description

一种基于切削过程仿真的单晶铜微铣削力预测方法

技术领域

本发明属于微小零件精密高效加工领域，特别涉及一种基于切削过程仿真的单晶铜微铣削力预测方法。

技术背景

随着科学技术的进步，航空航天等领域微小结构/零件除了要满足加工工艺上的要求，还要满足一些特殊要求，如需要有优质的导电性能。单晶铜由于具有一致的结晶取向，且只有一个晶粒，内部没有晶界，其导电性能超过一般的金属且具有良好的热、疲劳和蠕变性能，使其在航空航天等领域应用广泛。微铣削技术是加工单晶铜微小零件的有效手段。研究微小零件微铣削技术，对提高其加工精度及效率具有重要意义。微铣削过程中，微铣削力作为一个重要过程物理参数，研究微铣削力建模方法对于优化微小零件微铣削加工工艺、提高加工质量具有重要指导作用。

目前，微铣削力建模方法主要有力学解析法、有限元仿真法及智能算法建模法等。力学解析法考虑了刀具参数、工件材料屈服强度和切削条件等因素的影响，但建模过程中做出的假设及简化会使得预测精度降低且建模过程复杂。有限元仿真法节省成本，可以综合考虑刀具及工件材料的影响，但计算成本高，预测精度差。智能算法建模法需要大量试验样本，且无法考虑实际加工过程情况，导致微铣削力预测精度低。尽管微铣削力建模方法较多，但由于单晶铜材料的各向异性，不同晶体取向的单晶铜微铣削加工机理也不同，目前还缺乏切削过程仿真的单晶铜微铣削力预测方法。

卢晓红等人在《一种基于刀具磨损效应的微铣削力建模方法》中提出了一种基于刀具磨损效应的微铣削力建模方法。应用DEFORM有限元软件建立微铣削三维仿真模型。采用Johnson-Cook材料本构模型表征材料属性，并将刀具磨损效应引入微铣削力建模过程，实现一种模型对于微铣削力的精准预测，提高模型鲁棒性。但文献并没有开展有关单晶铜这类各向异性材料的研究。

尽管微铣削力建模方法较多，但由于单晶铜材料的各向异性，不同晶体取向的单晶铜微铣削加工机理也不同，目前还缺乏切削过程仿真的单晶铜微铣削力预测方法。

发明内容

为解决上述问题，本发明提出一种基于切削过程仿真的单晶铜微铣削力预测方法。以有限元切削仿真为基础，实现对于单晶铜微铣削力的准确预测，为单晶铜微铣削过程提供技术支撑，提高单晶铜微铣削加工精度与效率，具有实际应用价值。

本发明采用的技术方案是：

一种基于切削过程仿真的单晶铜微铣削力预测方法，由于单晶铜材料的特殊性，首先综合考虑了晶体变形运动学、晶体塑性本构关系、率相关晶体的硬化规律等因素对于微铣削单晶铜的影响；然后对本构方程的增量形式进行推导，采用Fortran语言编写单晶铜VUMAT用户材料子程序，将晶体塑性本构引入到有限元仿真中；最后采用ABAQUS软件对单晶铜微铣削过程进行仿真建模，实现单晶铜微铣削力预测，方法的具体步骤如下：

步骤1：构建晶体变形运动学模型

晶体塑性理论认为晶体的变形是由晶格畸变和位错滑移组成。其中晶格畸变可以通过连续介质力学中的弹性部分来描述，位错滑移可以通过连续介质力学来表示。如图1所示，晶体的总变形梯度F被分解为：

F＝F_e·F_p (1)

其中，F_e为弹性变形梯度，F_p为塑性的变形梯度。

滑移方向：

m*^(α)＝F_em^(α) (2)

其中，m^(α)和m*^(α)分别是晶体的晶格发生畸变前的第α滑移系滑移方向的单位向量和晶体的晶格发生畸变后的第α滑移系滑移方向的单位向量。

滑移面法向：

n*^(α)＝(F_e ^-1)^Tn^(α) (3)

其中，n^(α)和n*^(α)分别是晶体的晶格发生畸变前的第α滑移面滑移方向的单位法向量和晶体的晶格发生畸变后的第α滑移面滑移方向的单位法向量。

当前状态的速度梯度张量L为：

其中，和分别为F、F_e和F_p的导数；L_e和L_p分别是弹性速度梯度张量和塑性速度梯度张量。

建立各滑移系的滑移剪应变与晶体塑性变形的关系：

其中，为在滑移系α上的滑移剪切应变率，对所有激活的滑移系进行求和。

另外，速度梯度张量L还可以分解为对称部分D和非对称部分W，D为变形率张量，W为旋率张量；变形率张量D和旋率张量W也可以按弹性部分和塑性部分，即：

L＝D+W，D＝D_e+D_p，W＝W_e+W_p (6)

其中，D_e为弹性变形率张量，D_p为塑性变形率张量，W_e为弹性旋转张量，W_p为塑性旋转张量。

并满足：

因此，变形率张量D与旋转张量W的表达式如下：

这些晶体变形学的公式，将晶体变形中的滑移剪切率与宏观上的变形率相联系起来。

步骤2：晶体塑性本构关系的建立

由于弹性势的存在，确定了晶格的弹性变形率张量D_e和柯西应力焦曼率之间存在着如下的关系：

其中，σ为当前的柯西应力；I表示单位二阶张量；L表示一组弹性模量张量，包括：L_ijkl，L_jikl，L_ijlk，L_klij；它们具有完全对称的特性；为晶格旋转轴上的共旋转应力率，与材料旋转轴上的共旋转应力率有如下的关系：

其中：

由此，建立了应力率、变形率和滑移剪切应变率的关系，而应变率的得出先求出各滑移系的剪切应变率，而它是根据硬化规律来计算的。

步骤3：率相关晶体的硬化规律分析

单晶体滑移的程度，取决于外加载荷所产生的切应力的大小、晶体结构的类型以及开动滑移面相对于切应力的取向。不同取向的单晶体产生滑移时所需的拉伸载荷是不同的。可以认为滑移开始时，这一数值与单晶体一般应力应变曲线上的屈服应力相当。

当滑移面上存在临界分切应力(CRSS)，当滑移方向的切应力达到临界值，就会出现滑移现象。临界分切应力是晶体的固有属性，由晶体构型、材料组成和温度决定，与外力无关。晶体在滑移系α上的滑移剪切应变率由分解剪应力τ_α决定：

其中，为了简化表达形式将简化为 表示滑移系α的参考应变率，g_α表示滑移系的强度，f_α是一个无量纲通用函数，描述应变率与应力的依赖关系。

采用幂分布来描述多晶蠕变f_α(x)：

f_α(x)＝x|x|^n-1 (14)

其中，n为率敏感指数；x为函数自变量。

由式(13)所示，还要得到滑移系的强度g_α的值，其演化方程可由下式确定：

其中，为g_α的导数；h_αβ为滑移硬化模量，表示滑移系中的滑移剪应变对滑移系所产生的硬化。当α＝β时，h_αα为自硬化模量，即相同滑移系自身变形对自身的硬化；当α≠β时，h_αβ为潜硬化模量，即不同滑移系间产生的硬化；为晶体在滑移系β上的滑移剪切应变率。

用幂函数来描述自硬化模量h_αα：

其中，h₀为初始硬化模量，τ₀为初始屈服应力，τ_s为剪切屈服强度，γ为滑移系累计剪切应变，其中γ表达如下：

其中，t是时间变量。

潜硬化模量为：

h_αβ＝qh(γ)(α≠β) (18)

其中，q是常数。

晶体材料中三个强化阶段的滑移硬化模量：

h_αβ＝qh_αα(β≠α) (20)

其中，h_s为易滑移时的硬化模量；γ_α为在滑移系α上的累计剪切应变；γ_β为在滑移系β上的累计剪切应变；函数G与交互式硬化有关：

其中，γ₀为滑移系间相互作用达到峰值强度后的滑移量，每个分量f_αβ表示特定滑移作用的强度值。

步骤4：本构模型的增量形式推导

有限元计算采用增量法，因此首先应将本构关系写成增量形式。

根据率相关的切线系数方法，定义第α个滑移系上的剪切应变增量Δγ^(α)为：

其中，和分别为t+Δt时刻的剪切应变量和t时刻的剪切应变量。

写成线性插值方式为：

其中，θ为积分参数，取值范围为[0,1]，当θ为零时，对应为简单欧拉时间积分格式，一般为了保证具有较高的精度和稳定性，推荐在[0.5,1]的范围内取；和分别为t+Δt时刻的剪切率和t时刻的剪切率。

滑移剪切应变率是关于分解剪应力和当前强度的函数，将公式(23)进行泰勒展开，如下：

由公式(22-24)可推知：

引入各滑移系的施密特因子和张量

由晶体滑移的硬化方程，可知当前滑移系的强度增量Δg_α为：

由此可知，剪应力的增量Δτ_α为：

其中，和为各滑移系的施密特张量。

可得共回转应力增量Δσ_ij为：

其中，和为各滑移系参考基向量上的位移分向量；σ_jk、σ_ij和σ_ik为各滑移系参考基向量上的柯西分应力；Δε_kk和Δε_ij各滑移系参考基向量上的应变增量。

由给定的应变增量得出剪切应变增量：

其中，δ_αβ为克罗内克尔记号。

求解方程组(30)可解的该增量步中的滑移剪切应力的增量，带入式(27)、(28)、(29)，其他未知量均可求得。

步骤5：ABAQUS/VUMAT用户材料子程序编写，输入步骤1-4得到的参数，构建晶体弹塑性本构模型

采用Fortran语言编写单晶铜VUMAT用户材料子程序，VUMAT用户材料子程序流程如图2所示。第一步，从ABAQUS中读取材料参数PROPS(n)，定义积分点序号I并且令I＝1。第二步，定义滑移系组数、积分系数和旋转增量；计算滑移方向滑移面法向、滑移变形张量、状态变量、剪切应变和累计剪切应变；计算滑移旋转张量、剪切应变率和其导数；求出自潜硬化模量；更新应力和应变；进而判断I是否是最后一个积分点序号，若不满足，则执行I＝I+1，循环第二步，若满足条件，则继续向下输出；最后，将更新好的应力和应变输出到ABAQUS中，结束程序。

步骤6：基于有限元方法获得单晶铜微铣削力预测

测绘微铣刀几何结构，建立微铣刀模型，导入ABAQUS软件中，将其设置为刚体。建立三维加工工件模型，设置为弹塑性体。对模型进行网格划分，刀具模型网格类型采用R3D3型网格，工件模型网格类型采用C3D8R，并将刀具及工件接触处切削区域的网格局部加密。

有限元仿真中，为了保证仿真结果的准确性，需要选择合适的摩擦模型来反映实际加工中的摩擦。微铣削加工过程中切屑与前刀面的接触区域有着两种摩擦类型，分别为粘结区和滑移区。

摩擦特性表达如下：

τ_f＝τ_s(τ_f≥τ_s) (31)

τ_f＝μσ_n(τ_f<τ_s) (32)

其中，式(31)表示粘结区，式(32)表示滑动区，其中τ_f是摩擦剪切力，τ_s是工件的剪切屈服强度，μ是摩擦系数，σ_n是接触区域的正应力。根据摩擦剪切力与材料剪切屈服强度的比较，判断出当前处于的摩擦区。

采用单元删除的建模技术模拟材料去除。为实现单元删除，需要设置材料去除的准则。采用一种基于剪切应变的标准来解释材料去除过程中的各向异性：

max(γ-γ_cr,γ_sl,min-γ_sl,cr)≥0， γ_sl,min＝min(γ_α)，α＝1,2...N (33)

其中，γ_sl,cr和γ_cr分别是单滑移系统上剪切应变的临界值和所有滑移系统上的累积剪切应变，监测单个滑移系统上的剪切以及由于滑移系统引起的整体滑移；如果达到单个系统的临界值或累积滑移的临界值就删除该单元；γ_sl，min为滑移系α上累计剪切应变最小值。

设置主轴转速、每齿进给量、轴向切深等参数，即可对单晶铜微铣削力进行预测。

单晶铜微铣削仿真过程如图3所示，该参数下ABAQUS软件输出的三向切削力如图4所示，图中结果为采用ABAQUS后处理软件滤波之后的结果。

本发明的有益效果：本发明所述的基于切削过程仿真的单晶铜微铣削力预测方法无需大量实验而且对单晶铜这种各向异性材料适应性强，可实现单晶铜微铣削力预测，提高微铣削力的预测精度和效率，具有实际应用价值。

附图说明

图1是晶体变形运动示意图。

图2是VUMAT子程序流程图。

图3是单晶铜微铣削仿真过程。

图4是ABAQUS软件输出的切削力。

图5(a)是<100>晶向单晶铜微铣削Fx方向仿真值与实验值的对比图。

图5(b)是<100>晶向单晶铜微铣削Fy方向仿真值与实验值的对比图。

图5(c)是<100>晶向单晶铜微铣削Fz方向仿真值与实验值的对比图。

图6(a)是<110>晶向单晶铜微铣削Fx方向仿真值与实验值的对比图。

图6(b)是<110>晶向单晶铜微铣削Fy方向仿真值与实验值的对比图。

图6(c)是<110>晶向单晶铜微铣削Fz方向仿真值与实验值的对比图。

图7(a)是<111>晶向单晶铜微铣削Fx方向仿真值与实验值的对比图。

图7(b)是<111>晶向单晶铜微铣削Fy方向仿真值与实验值的对比图。

图7(c)是<111>晶向单晶铜微铣削Fz方向仿真值与实验值的对比图。

具体实施方式

下面结合附图和技术方案详细说明本发明的具体实施方式。

考虑到单晶铜微小零件微铣削过程中，微铣削力是一个重要的过程参量，因而研究单晶铜微铣削力建模方法对优化微小零件微铣削加工工艺、提高加工质量具有重要指导作用。此外，由于单晶铜材料的各向异性对于微铣削过程中的切削力变化有一定的影响。据此，针对单晶铜微小零件微铣削预测力难题，发明了一种基于切削过程仿真的单晶铜微铣削力预测方法。

以两刃平头铣刀微铣削直槽为例，采用的刀具是CrN涂层的超微粒子碳化钨两刃平头立铣刀，刀具的直径为0.4mm、刀具刃长为1.2mm、刃口圆弧半径1.9μm螺旋角20°、周刃前角9°、周刃后角12°，借助ABAQUS软件进行仿真，详细说明本发明实施过程。

首先确定选用的工件材料为单晶铜，其材料参数性能为密度：8,960kg/m³；弹性常数C₁₁:168,150Mpa；弹性常数C₁₂:121,400Mpa；弹性常数C₄₄:75,400Mpa；比热容8,960J/kg·℃；热膨胀系数0.000016/℃；热传导系数：385W/m·K；熔化温度1,356K。

将通过扫描电镜得到的图片导入到AutoCAD中进行精确描摹，得到其底面轮廓，将之导入Creo软件中，建立微铣刀的三维几何模型，采用有限元分析前处理软件Hypermesh对微铣刀模型进行网格划分，网格类型采用R3D3型网格，为了保证微铣刀的形状精度，在微铣刀的切削刃部分划分了较为密集的网格，刀具网格划分了23967个，最后导入至ABAQUS软件中设置为刚体。模型工件为阶梯状半环形体，设置为弹塑性材料。网格类型采用C3D8R，在与刀具接触的切削区域使用网格局部加密，网格数量为180540，网格最小尺寸为1.5μm。利用Fortran语言编写单晶铜VUMAT用户材料子程序，将晶体塑性本构引入到有限元仿真中。定义材料参数时带入上述材料参数以及率敏感指数n＝50，初始硬化模量h₀＝180MPa，易滑移时的硬化模量h_s＝24MPa，剪切屈服强度τ_s＝52MPa，初始屈服应力τ₀＝4MPa，常数q＝1。

ABAQUS软件中：在Property模块中，剪切摩擦系数设置为0.48，设置材料去除的准则，公式(33)中单滑移系统上剪切应变的临界值γ_sl,cr值为6.0，所有滑移系统上的累积剪切应变γ_cr的值为0.068。在Assembly模块中，调入工件和刀具模型。在模型中，通过调整铣刀和工件的相对位置确定切削深度和进给距离。定义边界条件，将工件侧面及底面自由度严格约束。在Step模块中，建立初始步后依次插入微铣削加工分析步、退刀分析步、约束转换分析步，过程类型选Dynamics，Explicit，场变量输出中勾选力，以便读取切削力值。选用了三种单晶铜(<100>,<110>,<111>)，上表面分别是(100)、(110)和(111)晶面，进给方向分别为<001>，<001>，在Job模块中创建数据，检查任务，设置刀具主轴转速为40,000rpm；设置每齿进给量为3μm/z；轴向切削深度为40μm。数据检查无误后，递交任务，进行有限元分析。在仿真完成后，输出X、Y、Z三个方向的切削力。

表1<100>晶向单晶铜微铣削实验与仿真的切削力对比

从图5(a)，5(b)，5(c)中可以看出，<100>晶向单晶铜仿真得到的微铣削力与实验得到的微铣削力基本吻合，但仍然存在偏差，如表1所示为<100>晶向单晶铜微铣削力的实验值和仿真值的对比，从表中可以看出最大相对误差为28.6％，平均相对误差为14.8％。因此验证了所建立的<100>晶向的单晶铜微铣削仿真模型的准确性。

表2<110>晶向单晶铜微铣削实验与仿真的切削力对比

如图6(a)，6(b)，6(c)所示，<110>晶向单晶铜在进给方向力Fx，径向力Fy和轴向力Fz三个方向力的仿真值与实验值对比。结果表明，仿真得到的微铣削力与实验得到的微铣削力基本吻合，但仍然存在偏差，如表2所示为<110>晶向单晶铜微铣削力的实验值和仿真值的对比，从表中可以看出最大相对误差为29.4％，平均相对误差为19.5％。因此验证了所建立的<110>晶向的单晶铜微铣削仿真模型的准确性。

表3<111>晶向单晶铜微铣削实验与仿真的切削力对比

如图7(a)，7(b)，7(c)所示，<111>晶向单晶铜在进给方向力Fx，径向力Fy和轴向力Fz三个方向力的仿真值与实验值对比。结果表明，仿真得到的微铣削力与实验得到的微铣削力基本吻合，但仍然存在偏差，如表3所示为<111>晶向单晶铜微铣削力的实验值和仿真值的对比，从表中可以看出最大相对误差为28.0％，平均相对误差为15.8％。因此验证了所建立的<111>晶向的单晶铜微铣削仿真模型的准确性。

综上，验证了所建立的<100>，<110>，<111>晶向的单晶铜微铣削仿真模型的准确性。

Claims (1)

1.一种基于切削过程仿真的单晶铜微铣削力预测方法，其特征在于，具体步骤如下：

步骤1：构建晶体变形运动学模型

晶体的总变形梯度F被分解为：

F＝F_e·F_p (1)

其中，F_e为弹性变形梯度，F_p为塑性的变形梯度；

滑移方向：

m*^(α)＝F_em^(α) (2)

其中，m^(α)和m*^(α)分别是晶体的晶格发生畸变前的第α滑移系滑移方向的单位向量和晶体的晶格发生畸变后的第α滑移系滑移方向的单位向量；

滑移面法向：

n*^(α)＝(F_e ^-1)^Tn^(α) (3)

其中，n^(α)和n*^(α)分别是晶体的晶格发生畸变前的第α滑移面滑移方向的单位法向量和晶体的晶格发生畸变后的第α滑移面滑移方向的单位法向量；

当前状态的速度梯度张量L为：

其中，和分别为F、F_e和F_p的导数；L_e和L_p分别是弹性速度梯度张量和塑性速度梯度张量；

建立各滑移系的滑移剪应变与晶体塑性变形的关系：

其中，为在滑移系α上的滑移剪切应变率，对所有激活的滑移系进行求和；

另外，速度梯度张量L可分解为对称部分D和非对称部分W，D为变形率张量，W为旋率张量；变形率张量D和旋率张量W可按弹性部分和塑性部分，即：

L＝D+W，D＝D_e+D_p，W＝W_e+W_p (6)

其中，D_e为弹性变形率张量，D_p为塑性变形率张量，W_e为弹性旋转张量，W_p为塑性旋转张量；

并满足：

因此，变形率张量D与旋转张量W的表达式如下：

步骤2：晶体塑性本构关系的建立

由于弹性势的存在，确定晶格的弹性变形率张量D_e和柯西应力焦曼率之间存在着如下的关系：

其中，σ为当前的柯西应力；I表示单位二阶张量；L表示一组弹性模量张量，包括：L_ijkl，L_jikl，L_ijlk，L_klij，具有完全对称的特性；为晶格旋转轴上的共旋转应力率，与材料旋转轴上的共旋转应力率有如下的关系：

其中：

步骤3：率相关晶体的硬化规律分析

晶体在滑移系α上的滑移剪切应变率由分解剪应力τ_α决定：

其中，将简化为 表示滑移系α的参考应变率，g_α表示滑移系的强度，f_α是一个无量纲通用函数，描述应变率与应力的依赖关系；

采用幂分布来描述多晶蠕变f_α(x)：

f_α(x)＝x|x|^n-1 (14)

其中，n为率敏感指数；x为函数自变量；

滑移系的强度g_α，其演化方程由下式确定：

其中，为g_α的导数；h_αβ为滑移硬化模量，表示滑移系中的滑移剪应变对滑移系所产生的硬化；当α＝β时，h_αα为自硬化模量，即相同滑移系自身变形对自身的硬化；当α≠β时，h_αβ为潜硬化模量，即不同滑移系间产生的硬化；为晶体在滑移系β上的滑移剪切应变率；

用幂函数来描述自硬化模量h_αα：

其中，h₀为初始硬化模量，τ₀为初始屈服应力，τ_s为剪切屈服强度，γ为滑移系累计剪切应变，其中γ表达如下：

其中，t是时间变量；

潜硬化模量为：

h_αβ＝qh(γ)(α≠β) (18)

其中，q是常数；

晶体材料中三个强化阶段的滑移硬化模量：

h_αβ＝qh_αα(β≠α) (20)

其中，h_s为易滑移时的硬化模量；γ_α为在滑移系α上的累计剪切应变；γ_β为在滑移系β上的累计剪切应变；函数G与交互式硬化有关：

其中，γ₀为滑移系间相互作用达到峰值强度后的滑移量，每个分量f_αβ表示特定滑移作用的强度值；

步骤4：本构模型的增量形式推导

根据率相关的切线系数方法，定义第α个滑移系上的剪切应变增量Δγ^(α)为：

其中，和分别为t+Δt时刻的剪切应变量和t时刻的剪切应变量；

写成线性插值方式为：

其中，θ为积分参数，取值范围为[0,1]；和分别为t+Δt时刻的剪切率和t时刻的剪切率；

滑移剪切应变率是关于分解剪应力和当前强度的函数，将公式(23)进行泰勒展开，如下：

由公式(22-24)可推知：

引入各滑移系的施密特因子和张量

由晶体滑移的硬化方程，可知当前滑移系的强度增量Δg_α为：

则，剪应力的增量Δτ_α为：

其中，和为各滑移系的施密特张量；

得到共回转应力增量Δσ_ij为：

其中，和为各滑移系参考基向量上的位移分向量；σ_jk、σ_ij和σ_ik为各滑移系参考基向量上的柯西分应力；Δε_kk和Δε_ij各滑移系参考基向量上的应变增量；

由给定的应变增量得出剪切应变增量：

其中，δ_αβ为克罗内克尔记号；

步骤5：ABAQUS/VUMAT用户材料子程序编写，输入步骤1-4得到的参数，构建晶体弹塑性本构模型

第一步，从ABAQUS中读取材料参数PROPS(n)，定义积分点序号I并且令I＝1；第二步，定义滑移系组数、积分系数和旋转增量；计算滑移方向滑移面法向、滑移变形张量、状态变量、剪切应变和累计剪切应变；计算滑移旋转张量、剪切应变率和其导数；求出自潜硬化模量；更新应力和应变；进而判断I是否是最后一个积分点序号，若不满足，则执行I＝I+1，循环第二步，若满足条件，则继续向下输出；最后，将更新好的应力和应变输出到ABAQUS中，结束程序；

步骤6：基于有限元方法获得单晶铜微铣削力预测

测绘微铣刀几何结构，建立微铣刀模型，导入ABAQUS软件中，将其设置为刚体；建立三维加工工件模型，设置为弹塑性体；对模型进行网格划分，刀具模型网格类型采用R3D3型网格，工件模型网格类型采用C3D8R，并将刀具及工件接触处切削区域的网格局部加密；

微铣削加工过程中切屑与前刀面的接触区域有着两种摩擦类型，分别为粘结区和滑移区；

摩擦特性表达如下：

τ_f＝τ_s(τ_f≥τ_s) (31)

τ_f＝μσ_n(τ_f<τ_s) (32)

其中，式(31)表示粘结区，式(32)表示滑动区，其中τ_f是摩擦剪切力，τ_s是工件的剪切屈服强度，μ是摩擦系数，σ_n是接触区域的正应力；

根据摩擦剪切力与材料剪切屈服强度的比较，判断出当前处于的摩擦区；

采用单元删除的建模技术模拟材料去除；为实现单元删除，需要设置材料去除的准则，采用一种基于剪切应变的标准来解释材料去除过程中的各向异性：

max(γ-γ_cr,γ_sl,min-γ_sl,cr)≥0，γ_sl,min＝min(γ_α)，α＝1,2...N (33)

其中，γ_sl,cr和γ_cr分别是单滑移系统上剪切应变的临界值和所有滑移系统上的累积剪切应变，监测单个滑移系统上的剪切以及由于滑移系统引起的整体滑移；如果达到单个系统的临界值或累积滑移的临界值就删除该单元；γ_sl，min为滑移系α上累计剪切应变最小值；

设置主轴转速、每齿进给量、轴向切深参数，即可对单晶铜微铣削力进行预测。