利用热超声波晶体塑性分析方法研究12Cr2Ni4A齿轮钢在振动辅助磨削过程中的残余应力及表面完整性

《Journal of Materials Research and Technology》：Thermo-ultrasonic crystal plasticity analysis of residual stress and surface integrity in vibration-assisted grinding of 12Cr2Ni4A gear steel

　　邹润祥|唐金元|周卫华|周远生|叶世杰
中南大学极端服役性能精密制造国家重点实验室，长沙410083，中国

**摘要**
超声振动辅助磨削（UVG）可以提高硬化齿轮钢的表面完整性，但从振动引起的软化到残余应力保持的介观机制仍不够明确。本研究为12Cr2Ni4A齿轮钢开发了一种热超声晶体塑性有限元方法（CPFEM）框架。将工艺规模的磨削力和热流转换为多晶子模型的等效移动机械和热边界载荷，同时将超声和热软化引入基于位错密度的滑移律中，作为有效的、实验约束的项。模拟结果表明，超声振动降低了临界解析剪应力（CRSS），促进了多个{110}<111>滑移系的协调激活，并增加了近表面层的位错储存。这些变化提供了一种从III型到II型再到I型的应力传递路径：局部亚晶粒应力不均匀性在晶粒尺度上重新分布，最终以更深的宏观压缩残余应力形式保留下来。与传统磨削（CG）相比，UVG下预测的峰值压缩残余应力达到-369.13 MPa，而100 μm处的残余应力仍然比CG下的更压缩（-49.76 MPa vs. -23.48 MPa）。该框架使用热压缩数据进行了校准，并通过超声磨削实验、电子背散射衍射/KAM/GND表征和X射线衍射残余应力测量进行了验证。通过连接工艺载荷、晶体变形、位错演化和保留应力，本研究为高强度钢磨削中的表面完整性控制提供了基于材料的基础。

**1. 引言**
超声振动辅助磨削（UVG）因能降低切削阻力并提高表面质量而在高强度齿轮钢的精密加工中越来越受到应用[1]。随着超声设备的发展和对高精度齿轮需求的增加，UVG过程中的残余应力与微观结构控制已成为一个实际的材料加工问题[2]、[3]、[4]、[5]。
磨削产生的残余应力是由机械载荷和快速加热/冷却耦合形成的非平衡应力场，其大小和深度分布影响齿轮的接触疲劳强度、疲劳寿命和耐磨性[6]、[7]、[8]。同时，严重的近表面塑性变形会改变位错密度、晶粒结构、相分数和织构[9]、[10]。除了残余应力外，加工表面形态还会影响后续的润滑接触响应，因为点接触中的热混合弹性流体动力润滑对表面纹理和局部接触几何形状非常敏感[11]。因此，一个有用的UVG模型应能将工艺载荷与缺陷结构联系起来，进而与保留的应力状态关联起来。
从材料角度来看，UVG是一个受晶格旋转、滑移活性、位错储存和应力重新分布控制的表面完整性过程。关键问题不仅在于振动是否改变了磨削力，还在于它如何改变受影响层内的晶体变形路径。然而，当前在这一领域的研究仍存在以下局限性：

**（1）超声-热机械耦合建模的局限性**
大多数现有的UVG模型将工件视为均匀连续体，主要预测宏观力或温度。这些模型很有价值，但它们无法直接显示高频振动和局部温度升高如何在晶粒尺度上与位错滑移相互作用[12]、[13]。这限制了它们解释UVG和CG之间残余应力差异的能力。
电子背散射衍射（EBSD）、扫描电子显微镜（SEM）和X射线衍射（XRD）等实验方法提供了表面状态和残余应力的直接证据。例如，Maedeh[14]使用XRD研究了磨削参数对硬质合金涂层残余应力的影响。这些测量对于验证至关重要，尽管它们通常捕捉的是最终状态而非瞬态滑移过程。
有限元模型已被用于重建磨削温度、力、相变和残余应力场。Zhao等人[15]开发了一种用于超声振动辅助成形齿轮磨削的热机械耦合模型，并报告了主轴速度和超声振幅对磨削热量、磨削力和表面压缩残余应力的强烈影响。其他研究考虑了磨削引起的相变、位错密度、动态再结晶、晶粒尺寸或微裂纹形成[16]、[17]、[18]。这些模型提供了有用的工艺级趋势，但大多数模型未能解决晶粒取向依赖的变形问题。
在更小的尺度上，分子动力学、元胞自动机和相场模型可以描述原子流动、晶粒细化和微观结构转变[19]、[20]、[21]、[22]、[23]。相场方法特别适用于由热梯度和储存能量驱动的微观结构演变。然而，它们与UVG的直接耦合仍然困难，因为所需的热动力学参数和计算成本较高。
还提出了基于热弹性塑性理论、相变或位错密度演化的残余应力分析模型[24]、[25]、[26]、[27]。它们的优点在于物理清晰性，但封闭形式解所需的假设通常简化了晶界约束、取向效应和局部塑性应变不均匀性。
最近也有专门针对超声辅助加工的残余应力预测模型。Zhang等人[28]研究了轴承钢超声振动辅助磨削中的残余应力生成，Zhao等人[29]表明晶体取向影响微磨削残余应力，Feng等人[30]预测了超声振动辅助铣削中的残余应力分布。这些研究表明，超声激励可以纳入残余应力模型。然而，它们仍然无法充分揭示滑移系竞争、位错密度梯度以及晶粒间的应力传递。
最近的振动辅助加工研究进一步显示了当前问题的相关性。据报道，超声振动辅助的镍基单晶叶片榫头的轮廓磨削可以改善磨料轨迹、表面形态和亚表面完整性[31]。一篇关于航空发动机中高强度韧性合金的综述强调，齿轮、轴、叶片和盘需要能够控制高力、热量和表面完整性的切削技术[32]。对于连续纤维增强金属基复合材料，从能量角度分析了超声磨削以明确力和异质去除机制[33]。在γ-TiAl合金磨削中，超声振动辅助的高效深度磨削降低了磨削温度和力，同时改变了切屑形成和表面形态[34]。这些研究证实了振动辅助的工艺级优势，但它们仍未阐明振动修改的载荷如何通过滑移系活性和位错储存传递到高强度齿轮钢中。
因此，剩余的空白是介观尺度上的。实验确定了最终的微观结构，连续介质模拟提供了工艺场，分析模型阐明了趋势；但这些方法单独使用时都无法完全描述振动修改的热机械载荷如何激活滑移系并在多晶齿轮钢中重新分配残余应力。晶体塑性框架非常适合填补这一空白[35]、[36]。

**（2）介观机制的空白**
晶体塑性有限元方法（CPFEM）通过为有限元分配晶粒尺度的取向和滑移系，将晶体滑移、位错演化和宏观响应联系起来[37]、[38]、[39]。它已被用于分析不同合金的纳米压痕、断裂、孪生和空洞形成[40]、[41]、[42]、[43]，显示出其将局部变形机制与可测量的机械行为联系起来的能力。
最近的基于位错密度的CPFEM公式包括温度和应变率效应，并已扩展到相变辅助塑性和应变梯度效应[44]、[45]、[46]、[47]。这些发展支持将CPFEM用于同时发生高热梯度和严重塑性变形的表面加工问题。
对于金属材料，CPFEM可以比各向同性连续体模型更明确地表示晶体各向异性和局部应力集中[48]、[49]。这种能力对于磨削非常重要，因为即使在相同的外部载荷下，相邻晶粒也可能经历不同的解析剪应力。对焊料合金、钢、钛合金和双相钢的研究表明，CPFEM可以捕捉晶粒尺度的塑性和损伤相关特征[50]、[51]、[52]、[53]、[54]、[55]。
CPFEM在磨削中的应用仍然有限，尤其是在UVG方面。难点在于接触区的快速移动、同时的机械和热载荷以及额外的超声激励。现有研究描述了UVG的力、温度[56]、[57]和残余应力[58]、[59]，但从振动到滑移激活、位错储存和残余应力保持的晶粒尺度路径仍不清楚。
本研究的目的是为12Cr2Ni4A齿轮钢建立这一缺失的路径。基于位错密度的CPFEM与工艺得到的磨削力和热流相结合，从而将宏观的轮-工件相互作用转化为晶粒尺度的滑移和缺陷演变。将超声和热软化项引入滑移律中，以估计振动如何改变CRSS、滑移系激活、位错储存和近表面层的残余应力形成。与传统热机械残余应力模型相比，当前框架不仅输出应力大小，还确定了哪些晶粒取向和滑移系承担塑性变形，以及III型局部应力如何转化为II型晶粒应力，以及这些应力在卸载后如何保留为I型残余应力。

**2. 齿轮UVG的热超声晶体塑性有限元建模**
本节描述了UVG-CPFEM框架。在本文中，除非另有说明，UVG特指齿轮的超声振动辅助磨削。首先，使用齿轮磨削运动学计算等效力和热载荷。然后，将晶体塑性本构律扩展以包含超声和热软化项。材料参数来自EBSD信息和热压缩测试。最后，使用两个Neper生成的多晶子模型：一个用于局部应力和微观结构演变的表面模型，以及一个用于残余应力分布的深度模型。

**2.1. 超声辅助齿轮磨削力和热载荷计算**
图1显示了轴向超声振动辅助齿轮磨削的配置。在磨削过程中，砂轮以速度ns旋转，并以振幅A和频率f轴向振动，同时沿磨削方向以速度vw进给并径向切割至深度ar。由于渐开线齿面曲率不断变化，因此使用动态加工坐标系来描述砂轮和局部齿面接触区的相对位置。

下载：下载高分辨率图像（355KB）
下载：下载全尺寸图像
图1. 12Cr2Ni4A齿轮的轴向超声磨削示意图。
定义了一个笛卡尔加工坐标系xp-yp-zp，以有效渐开线轮廓的起点Op为原点。xp轴平行于齿轮轴，yp轴垂直于砂轮旋转轴，zp轴遵循右手规则。坐标变换后，该局部系统中的渐开线轮廓表示为：
(1)
其中rb是基圆半径，β是对应于基圆螺距的中心角的一半，φ=θk+αk。θk是渐开线展开角，αk是压力角。
齿面上的有效磨削深度an与径向磨削深度ar通过局部轮廓斜率相关：
(2)
有效砂轮半径随渐开线轮廓变化，因为局部齿面点与砂轮轴之间的距离不是恒定的。因此，有效砂轮半径表示为：
(3)
其中Rgmax是最大砂轮半径。
有效砂轮表面速度也随局部有效半径变化。它与砂轮旋转速度的关系为：
(4)
其中vs(φ)表示砂轮在渐开线磨削位置的线速度，ns是砂轮的旋转速度。
然后从局部砂轮半径和法向深度获得等效接触弧长：
(5)
在UVG中，间歇性分离和高频相对运动减少了磨料颗粒与工件之间的有效接触面积。因此，摩擦系数低于CG中的摩擦系数。根据参考文献[60]、[61]，超声和传统摩擦系数之比为：
(6)
其中ξ=vs/2πAf是局部砂轮表面速度与振动速度之比，θ是振动方向与宏观速度之间的角度。在轴向超声振动磨削（UVG）中，θ = π/2。μc 和 μu 分别是CG（常规磨削）和UVG的摩擦系数；α 和 β 是摩擦二项式表达式中的材料相关常数；ds 是活性磨料颗粒的等效直径；vs 是局部轮面速度；P0 是CG中磨料-工件界面的平均法向接触压力；vw 是工件进给速度。这些定义与方程（7）、（8）保持一致，其中Ft 表示切向磨削力，Qw 表示进入工件的热流。然后可以通过考虑磨削参数、局部齿面几何形状和修正后的摩擦状态来计算切向磨削力Ft。详细推导过程参考文献[62]：(7)其中a1 和 a2 是从实验中获得的材料系数。磨削热主要来自磨料-工件界面的摩擦和近表面材料的塑性变形。在本文使用的干磨条件下，产生的热量分为传导到齿轮中的热量、通过切屑带走的热量、散失到周围空气中的热量以及传导到磨料颗粒中的热量。空气对流遵循牛顿冷却定律[63]，工件/磨料的热量分配遵循经典的磨削热量分配理论[64]、[65]、[66]。单颗粒磨削中进入工件的热流表示为：(8)其中Hu 是UVG系数；Ft 是切向磨削力；Bw 是磨削宽度；f(t) 是时间移位函数；ρw 和 cw 分别是工件的密度和比热容；Tm 是熔点；λw 和 λg 分别是工件和磨料颗粒的热导率；ρg 和 cg 分别是磨料颗粒的密度和比热容。这个热源表达式为CPFEM子模型提供了瞬态热流，并使用了测量的齿轮材料属性以及CBN磨料的特性。在这个框架中，参考文献[62]、[66]提供了基线的磨削力和热量分配关系。超声波的贡献通过ξ=vs/2πAf中的A和f定量引入，这改变了方程（6）中的μu/μc。然后使用修正后的摩擦状态在方程（7）中计算Ft，并将Ft传递到方程（8）中以获得DLOAD/DFLUX边界条件下的Qw。因此，振动作为数值输入被包含在机械和热负荷中。

2.2. UVG晶体的塑性本构模型构建
晶体变形通过将总变形梯度分解为热变形、弹性变形和塑性变形来描述。在图2中，Fθ 表示由热膨胀引起的晶格变形，Fe 表示弹性晶格畸变，Fp 表示主要由滑移引起的塑性变形。因此，总变形梯度 F 可以表示为这三个张量的乘积[67]。

下载：下载高分辨率图像（348KB）
下载：下载全尺寸图像
图2. 晶体运动和变形的示意图。
晶体从参考配置A0开始，移动到热膨胀配置A1，经过塑性和弹性变形后达到当前配置A。滑移系统α的滑移方向和滑移面法线最初由m0(α)和s0(α)定义。经过晶格旋转和变形后，它们在中间配置变为mα和sα[67]：
(9)
其中m0(α)和s0(α)是滑移系统α的初始单位向量。更新后的向量被归一化以保持单位长度。
如果假设热变形是各向同性的，则当前配置中的速度梯度L可以分解为弹性部分、塑性部分和热部分[44]：
(10)
叠加的点表示材料的时间导数。Le、Lp 和 Lθ 分别是弹性速度梯度、塑性速度梯度和热速度梯度。αT 是热膨胀系数，ΔT 是相对于参考状态的温度变化，γα 是滑移系统α的塑性剪切率，Mα=mα?sα 是Schmid张量。
对于没有热旋转变的静水热膨胀（Wθ=0），L 分解为对称变形率张量D和反对称旋转变张量W：
(11)
对称变形率张量D包括弹性部分和塑性部分，表示为：
(12)
反对称旋转变张量W包括弹性部分和塑性部分，表示为：
(13)
这种运动学分解将可恢复的弹性变形与晶体滑移系统上的剪切率联系起来。
此外，在参考配置中定义的拉格朗日应变表示为：
(14)
通过对拉格朗日应变进行时间微分，并将其与参考配置中的第二Piola-Kirchhoff应力S共轭，可以得到格林应变率：
(15)
根据能量守恒定律，参考配置中的单位体积变形能量与当前配置中相应体积的变形能量相同，即：
(16)
热机械耦合的CPFEM通过UMAT子程序在ABAQUS/Standard中实现。选择ABAQUS是因为其UMAT接口允许隐式集成基于位错密度的晶体塑性定律和内部变量，同时其接触、有限旋转和大变形能力非常适合微观力学模拟。
首先使用经典的幂律关系描述滑移系统上的剪切应变率，该关系基于解析剪切应力τ(α)和当前滑移阻力g(α) [68]：
(17)
其中γ0 是参考剪切应变率（s-1），τ(α) 是解析剪切应力，g(α) 是当前滑移阻力，mt 是速率敏感指数。符号函数sgn(x)在x为负、零或正时分别给出-1、0或1。
对于UVG，超声波振动和温度变化会影响滑移和硬化响应。因此，通过引入超声波软化项Usoft和热软化项Tsoft来修改经典的剪切率定律。修改后的剪切应变率为：
(18)
其中Usoft 是超声波软化项，Tsoft 是热软化项。在本研究中，Usoft被用作振动辅助位错移动性的有效表示，而不是作为独立的位错动力学定律。

Usoft的物理基础。引入超声波项是为了表示高频振动导致的流动阻力降低。在位错尺度上，循环超声波应力可以促进从弱障碍物上的解离，减少滑移的时间平均障碍，并促进有利取向的滑移系统的早期激活。建立单独的振动-位错相互作用本构定律需要原位的超声波力学数据或位错动力学校准，而这些数据对于目前热处理的12Cr2Ni4A钢是不可用的。因此，Usoft被有意用作有效的、有文献支持的软化因子，并通过EBSD和XRD在实际UVG条件下的验证来评估其适用性。
在屈服后，塑性应变增加了位错密度并提高了进一步滑移的阻力。滑移阻力gα与位错密度ρα有关，通过总和所有滑移系统之间的相互作用来考虑这些相互作用[70]：
(19)
这里，τα0表示摩擦应力、晶界障碍和固溶强化；μ是剪切模量；Bur是Burgers矢量；ρα是滑移系统α上的位错密度；Gαβ是连接系统β上的位错密度与系统α上的滑移阻力的相互作用矩阵项。
方程（19）主要通过位错密度来表达滑移阻力。没有添加单独的Hall-Petch项，因为在多晶模型中平均晶粒尺寸是固定的，τ0是使用相同的材料状态拟合的。如果将来研究晶粒尺寸效应，则应明确添加Hall-Petch或GND基强化项。
位错密度在滑移过程中通过乘法和湮灭发生变化。使用Kocks-Mecking关系[71]来描述这种演变：
(20)
其中k1是位错乘积系数，k2是位错湮灭系数。
然后可以将位错密度的演变与滑移阻力率联系起来，得到硬化关系[45]：
(21)
方程（9）、（10）、（11）、（12）、（13）、（14）、（15）、（16）、（17）、（18）、（19）、（20）、（21）构成了齿轮UVG的CPFEM分析的本构基础。
耦合的热机械问题还需要线性动量和能量的平衡。在准静态假设下，动量平衡为：
(22)
其中σ是Cauchy应力张量，Vb是单位变形体积的体力，?·是空间散度算子。
热问题由能量平衡控制。在没有外部体积热源的情况下，傅里叶热传导表示为：
(23)
晶体塑性滑移改变了晶格取向，但不改变弹性本构形式。弹性关系表示为：
(24)
其中De是弹性变形率张量，I是四阶单位张量，L是弹性模量张量，σ?e是Cauchy应力的Jaumann率。Jaumann率是与晶格旋转相关的共旋应力率[72]：
(25)
因此
(26)
最终得到晶体塑性本构关系为：
(28)
这个关系将应力率、变形率和UMAT实现中使用的滑移剪切率联系起来。

2.3. 12Cr2Ni4A材料的晶体塑性模型参数确定
使用热压缩试验来确定12Cr2Ni4A钢的温度依赖性塑性流动行为。选择了Gleeble 3800热机械模拟器，因为它可以在压缩过程中控制温度、应变和应变率。合金成分列在表1中。圆柱形试样的直径为6毫米，高度为9毫米，首先在815°C下保持100分钟进行热处理，然后进行油淬火，在-80°C下进行100分钟的深冷处理，最后在150°C下进行180分钟的回火（图3）。
表1. 12Cr2Ni4A的化学成分
元素 CMn Si SP Cr Ni
百分比(%) 0.1～0.15 0.3～0.6 0.17～0.37 ≤0.015 ≤0.025 1.25～1.75 3.25～3.75
下载：下载高分辨率图像（459KB）
下载：下载全尺寸图像
图3. 热压缩试验：(a) 试验样品；(b) 样品热处理过程；(c) 热压缩过程。
热压缩过程如图3(c)所示。由于本研究的磨削温度接近300°C，因此将变形温度设置为225°C、325°C和425°C。总压缩变形为69%，应变率为0.01 s?1。样品以10°C/s的速度加热，在目标温度下保持2分钟以减少热梯度，压缩后立即水淬火。压缩后的均匀桶形表明摩擦条件得到了适当的控制。
图4(c)显示了12Cr2Ni4A钢在225°C至425°C之间、应变率为0.01 s?1时的平均应力-应变曲线。每种条件至少重复三次。重复之间的小散布表明数据适合参数识别。
下载：下载高分辨率图像（878KB）
下载：下载全尺寸图像
图4. 12Cr2Ni4A的真应力-应变曲线：(a) 12Cr2Ni4A的初始微观结构相图；(b) 热压缩实验的晶体塑性模型和模拟结果；(c) 不同温度下的12Cr2Ni4A的真应力-应变曲线；(d) 在5 μm振幅下进行超声波软化校正后的真应力-应变曲线；(e) 模拟和实验真应力-应变曲线的比较。
实验平台无法在压缩过程中同时施加超声波振动和高温。因此，使用参考文献[73]中提出的超声波软化校正将传统的热压缩本构响应转换为与UVG相关的有效响应。这种转换仅限于相同的钢材等级、相同的20 kHz频率以及参考文献覆盖的应变率范围。20 kHz频率也是本工作中使用的换能器-horn系统的共振频率，在该频率下振幅稳定性和能量传递最高。
这种处理并不声称基于文献的校正可以逐点再现完整的热超声波流变特性。它为当前的UVG条件提供了一个有效的本构桥梁。剩余的不确定性通过明确报告残余应力误差并通过将模型呈现为更适用于机制解释和趋势预测而不是通用绝对应力预测来处理。
在0.01 s?1的热压缩试验中，流动应力在20 kHz振动和5 μm振幅下以规律的方式减小（图4(d)）。这一趋势与报告的超声波软化效应一致，表明超声波能量有助于位错运动，特别是在高温下。
12Cr2Ni4A钢的CPFEM参数分为四组：
(1) 单晶弹性常数。BCC弹性矩阵由C11、C12和C44定义。初始材料的EBSD（图4(a)）显示钢主要是体心立方结构，含有少量的FCC相和渗碳体。CPFEM模型因此使用BCC晶格作为基体相，而在当前的介观尺度模型中，次要相没有分配单独的本构律。这种简化使模型专注于基体主导的滑移行为。少量的FCC分数和渗碳体颗粒可能会影响相界附近的局部绝对值，但由于EBSD图显示基体以BCC为主，并且验证重点关注平均残余应力和KAM/GND趋势，因此它们被视为次要因素。BCC晶体的弹性基体C表示为：(29)公式中的C11、C12和C44进一步定义为：(30)其中Ew是12Cr2Ni4A的弹性模量，ν是泊松比，这些数据来自单轴热压缩试验。(2)超声波和热软化因子。热软化项Tsoft遵循Johnson的温度依赖公式，而超声波软化项Usoft描述了在超声波能量作用下的流动阻力减小。这些项在方程(31)、(32)中给出：(31)(32)其中θtransition是转变温度，θmelt是熔化温度，pt是热软化指数，λsoft是超声波软化系数[73]。对于12Cr2Ni4A钢，λsoft=5.84(ε)-0.00847(E/σi)1/2。这个系数被视为有效的振动-位错相互作用项：周期性超声波加载有助于位错解离，并在给定的解析剪切应力下增加滑移的概率，但完整的振动-位错相互作用律需要额外的原位超声波力学试验，超出了现有的校准数据范围。(3)塑性流动参数。这些参数包括速率敏感指数mt和参考剪切应变率。在本研究中，将其设置为热压缩试验中使用的应变率。(4)塑性硬化参数。这些参数包括初始临界解析剪切应力τ0、伯格斯矢量Bur、初始位错密度ρ0、剪切模量μ以及位错增殖/湮灭系数k1和k2。像mt、τ0、k1和k2这样的参数无法直接以足够的可靠性进行测量。因此，通过最小化模拟和测量压缩曲线之间的均方误差（MSE）来确定这些参数。初始搜索范围是根据基于位错密度的CPFEM研究[69]、[72]、[74]选定的，最终的参数集对应于最低的MSE。(33)其中Nsim表示应力-应变曲线上的数据点数量；σi exp表示实验应力-应变数据点；σi sim表示模拟应力-应变数据点。校准后的本构参数列在表2中。C11、C12和C44是根据测量的杨氏模量和泊松比使用方程(30)计算得出的。Bur是根据BCC晶格参数计算得出的。其余参数，包括ρ0、τ0、k1和k2，是通过反演拟合热压缩曲线获得的。表2. 12Cr2Ni4A钢CPFEM模型中使用的晶体塑性参数符号描述参数值单位C11弹性常数280.12GPaC12弹性常数120.05GPaC44弹性常数75.41GPamt速率敏感指数10-参考剪切应变率0.01s-1ρ0初始位错密度7.72×107mm-2τ0临界解析剪切应力（CRSS）425MPaBur伯格斯矢量2.58×10-7mmμ剪切模量74980MPak1位错增殖系数6.54×105mm-1k2位错湮灭系数2.27-图4(e)比较了模拟和测量的压缩曲线。在测试的温度范围内，平均相对误差为1.8%，表明校准后的参数集能够捕捉用于模型识别的热压缩响应。因此，表2中的温度独立弹性和硬化常数应解释为225-425°C校准窗口内的有效参数，而不是适用于所有磨削温度的通用常数。2.4. 建立高强度合金钢12Cr2Ni4A的多晶模型为了使CPFEM计算可行，只对齿轮齿面的局部区域进行建模。由于建模区域处于微米级别，其大小远小于渐开线齿面的曲率半径。因此，将局部表面近似为平坦的。这种近似降低了模型复杂性，同时保留了感兴趣的晶粒尺度变形机制。多晶集合是使用Neper生成的。材料被建模为BCC 12Cr2Ni4A钢。BCC晶体具有{110}<111>、{112}<111>和{123}<111>滑移系；主要的{110}<111>系包含12个滑移系统，通常在BCC钢中主导塑性流动。为了在保留主要塑性机制的同时限制计算成本，本研究重点关注表3中列出的{110}<111>滑移系统。忽略{112}<111>和{123}<111>系可能会略微低估某些晶粒在高温或高应变率下的局部塑性应变，但不预期会改变此处讨论的滑移激活、位错存储和残余应力保留的整体趋势。表3. 12Cr2Ni4A的12组滑移系统滑移系统编号滑移面滑移方向1{1 1 0}<-1 1 1>2{1 1 0}<1 -1 1>3{1 -1 0}<1 1 1>4{1 -1 0}<1 1 -1>5{1 0 1}<1 1 -1>6{1 0 1}<-1 1 1>7{1 0 -1}<1 1 1>8{1 0 -1}<1 -1 1>9{0 1 1}<1 1 -1>10{0 1 1}<1 -1 1>11{0 1 -1}<-1 1 1>12{0 1 -1}<1 1 1>晶粒取向由晶体坐标系与宏观样品坐标系之间的关系描述。如图5(d)所示，RD、TD和ND定义了全局轴，而[100]、[010]和[001]定义了晶体轴。使用米勒指数为每个晶粒分配晶体学平面和方向：(34)下载：下载高分辨率图像（1MB）下载：下载全尺寸图像图5. 12Cr2Ni4A多晶模型的建立：(a)不同晶粒数量的多晶模型；(b)不同网格划分类型的多晶模型比较；(c)晶体结构类型；(d)晶格坐标系与参考坐标系之间的关系。每个晶粒都被分配了一个由Neper生成的唯一随机取向，并通过Python脚本导入ABAQUS。由于热处理后的12Cr2Ni4A钢在选定的建模区域内没有显示出强烈的初始优先取向，因此使用了随机纹理集合。初步的随机种子检查显示平均I型残余应力剖面的变化小于4%，表明集合大小适合当前的预测。Python脚本读取Neper生成的INP文件，将取向参数和晶体塑性常数编译成160×1的材料参数矩阵（表4），并将用户材料属性分配给ABAQUS中的每个晶粒。表4. 晶体塑性参数矩阵参数对应的矩阵位置材料参数名称Props(1)-Props(21)晶体弹性模量Props(25)-Props(56)滑移系统参数Props(57)-Props(72)晶体初始取向Props(73)-Props(96)滑移速率参数Props(97)-Props(144)基于位错的硬化参数Props(145)-Props(152)时间积分方案和有限变形参数Props(153)-Props(160)迭代方法参数使用四面体C3D4元素，因为它们比结构化的六面体元素更适合不规则的晶界。这种选择改善了晶界附近局部应力和应变梯度的表示，这对本研究至关重要。使用粗、中、细三种C3D4网格检查了网格敏感性，分别大约有1.5×10^4、4.8×10^4和1.2×10^5个元素。中等网格和细网格之间的最大冯·米塞斯应力差异小于3.2%，而细网格所需的计算量大约是前者的四倍。因此，在后续模拟中使用了中等网格。图6显示，在弹性阶段晶粒数量的影响很小，但会影响屈服和塑性流动。当晶粒数量超过约200个时，应力-应变曲线趋于收敛。因此，表面模型包含216个晶粒，这平衡了计算成本和统计稳定性。下载：下载高分辨率图像（195KB）下载：下载全尺寸图像图6. 在0.01 s-1的压缩应变率下不同晶粒数量的应力-应变曲线。图7总结了两个多晶子模型。表面模型尺寸为15 μm × 15 μm × 15 μm，包含216个平均尺寸约为2.5 μm的晶粒。它用于研究局部应力、滑移和位错演变。深度模型尺寸为100 μm × 50 μm × 3 μm，包含800个平均晶粒尺寸相同的晶粒。它用于计算沿深度方向的残余应力。由于3 μm的厚度接近平均晶粒尺寸，两个厚度方向的面受到对称位移条件的约束，以近似代表性的平面应变切片并抑制人为的自由表面松弛；没有使用周期性边界条件。两个模型都使用C3D4元素和随机晶粒取向，材料参数来自表2。边界条件如下：(1)约束：底部表面完全固定以表示工件支撑，在侧面表面施加对称位移约束以减少边界引起的应力奇异性。(2)初始应力：施加一个随深度衰减的压缩预应力场，表示为σyy(z)=-350×exp(-z/δ) MPa，其中δ = 100 μm。该场代表了通过XRD在预磨试样上测量的残余应力。采用指数衰减是因为磨削和预应力表面通常在表面附近显示高压缩应力，然后随深度衰减[75]、[76]。(3)加载：将过程衍生的机械载荷（DLOAD）和热通量（DFLUX）施加到顶部表面，并根据方程(8)计算。由于实际接触宽度远大于子模型尺寸，移动载荷宽度设置为模型长度的1.5倍。数值检查显示在稳态阶段应力误差小于5%，这对于局部子模型是可以接受的。下载：下载高分辨率图像（571KB）下载：下载全尺寸图像图7. 晶体模型尺寸和载荷应用的示意图：(a) 15 μm*15 μm*15 μm 表面多晶模型；(b) 100 μm*50 μm*3 μm 深度多晶模型。模拟在ABAQUS/Standard中运行。其他有限元软件也可以解决耦合的热机械问题，但选择了ABAQUS，因为其UMAT接口非常适合隐式晶体塑性积分，具有位错密度等内部变量。其有限旋转和大变形处理也支持严重表面载荷的介观尺度模拟。在CPFEM域中，轮子和单个磨料颗粒没有明确网格化。相反，它们的作用被转移到子模型中，作为等效的移动机械和热边界载荷。这种宏观到介观的策略保留了晶粒尺度分辨率，同时保持了三维计算的可行性。因此，该模型应被视为局部准稳态磨削子区域，而不是完整的轮子-工件接触模拟。UMAT在每个积分点接收应力、应变增量、状态变量、材料参数和晶粒取向信息。它根据上述本构方程更新滑移系统剪切、位错密度、滑移阻力和应力。实现顺序如图8所示，并遵循标准的隐式CPFEM积分程序[77]。下载：下载高分辨率图像（686KB）下载：下载全尺寸图像图8. 轮齿超声磨削的塑性模型计算和实现过程。3. 实验细节本节描述了用于验证UVG-CPFEM模型的实验。使用自开发的超声齿轮磨削平台在单因素参数变化下磨削12Cr2Ni4A齿轮。通过热电偶、XRD和EBSD分别测量磨削温度、残余应力和微观结构。3.1. UVG实验设计UVG设置如图9(a)所示。超声发生器、换能器和喇叭安装在垂直加工中心的主轴上，以对磨轮施加轴向振动。发生器提供20 kHz的信号，该信号由换能器转换为机械振动，并由喇叭放大。在磨削之前，使用激光振动计校准了轮子的振幅，并将其设置为5 μm以进行验证测试（图9(b)）。工件材料为12Cr2Ni4A钢，齿轮参数列在表5中。下载：下载高分辨率图像（903KB）下载：下载全尺寸图像图9. 超声波齿轮磨削实验平台：(a) 超声波齿轮磨削设备；(b) 磨轮超声振幅测量设备。表5. 齿轮特性参数。齿轮特性参数值模块（mm）3.8788齿数31压力角（°）28齿轮宽度（mm）45超声频率固定在20 kHz，因为换能器-喇叭系统在这个频率下工作在共振状态。在这种状态下，振幅稳定，可以通过发生器进行调整。因此，频率被视为系统参数，而不是自由变化的磨削参数。第4.1.4节中通过数值方法检查了频率的影响。使用单因素实验设计来测试磨削深度、轮速和进给速度的影响，同时保持其他变量不变。详细的参数矩阵显示在表6中。表6. 超声波振动磨削测试计划：'/'表示常规磨削条件，即没有施加超声振动振幅或频率。测试组编号磨削深度（μm）磨轮速度（r/min）进给速度（mm/min）超声振幅（μm）超声频率（kHz）15/10/15150010005202101000/1500/200010005203101500500/1000/150052041015001000//首先在数值敏感性分析（第4.1.4节）中检查了振幅和频率的范围。该分析表明，在测试范围内，振幅对残余应力的影响比频率更强，支持所选的实验设置。3.2. 磨削温度的测量磨削温度使用接触式热电偶（图10(a)）测量。在封闭的轮齿界面难以应用非接触方法。盲孔是通过电火花加工（EDM）从齿轮背面加工出来的，孔底距离加工表面在100 ± 50 μm范围内。热电偶使用高温陶瓷粘合剂固定在孔中，温度信号通过MC-USB-1616HS数据采集设备进行收集[78]。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像

图10. 超声波齿轮磨削和磨削温度及残余应力测量装置：(a) 磨削温度测量 (b) 电解抛光仪器 (c) 线切割后的齿轮样品及测量点 (d) 残余应力测量。

3.3 残余应力测量
磨削后，通过线切割电火花加工（EDM）获得包含加工齿面的块状样品（图10(c)）。表面残余应力通过XRD进行测量。该方法通过应力引起的衍射峰位移来确定残余应力，适用于磨削钢的近表面层。XRD测量使用了配备Cr靶材的Proto iXRD衍射仪。主要设置参数为25 kV、5 mA、波长2.291 ?、1 mm准直器直径、最大β角20°、3秒曝光时间和七次曝光。L1和R1探测器在156.41°的布拉格角下进行校准。样品在测量前用无水乙醇进行超声清洗。测量点位于齿尖、齿根和每个齿面的中间区域（图10(c)）。每个点至少测量三次并取平均值。由于渐开线轮廓导致齿面法线方向变化，因此沿齿宽方向测量残余应力，该方向上的表面更接近平面。

通过电抛光去除一层材料后，测量深度依赖的残余应力。材料以10 μm的步长从0去除到60 μm，然后以20 μm的步长从60 μm去除到100 μm。这一过程产生了深度达100 μm的残余应力分布。

4 结果与讨论
首先研究了UVG-CPFEM模型预测的晶粒尺度力学响应，包括应力不均匀性、滑移系激活、位错密度演变和残余应力传递。然后分析了磨削参数和超声参数的影响。最后，将预测结果与磨削实验中的EBSD和XRD测量结果进行了比较。

4.1 基于UVG-CPFEM的UVG残余应力形成机制
4.1.1 UVG的力学行为分析
UVG模拟显示12Cr2Ni4A钢中的应力场高度不均匀。在移动的DLOAD和DFLUX边界条件下，最大应力位置随磨削载荷变化，而局部应力水平强烈依赖于晶粒取向。在顶部0-15 μm层中，<110>取向晶粒的平均冯·米塞斯应力约为450 MPa，<100>取向晶粒的平均冯·米塞斯应力约为300 MPa，取向相关的波动约为35%（图11）。

图11. 表面多晶模型：UVG过程中的瞬态应力场等高线图。

图12显示了深度模型中冯·米塞斯应力的时间演变。在t=7.5×10-3 s时，表面应力迅速上升，相邻晶粒之间的应力差可达约300 MPa。随着载荷的移动，应力集中在晶界和三重晶界附近。卸载后，硬度较高的晶粒释放应力的速度较慢，留下不均匀的残余应力场。这种行为解释了宏观磨削载荷为何会产生晶粒尺度的应力分布而不是均匀应力状态。

图12. 深度多晶模型：UVG过程中的瞬态应力场等高线图。

局部响应还取决于晶界约束。如图13所示的例子中，在晶界附近约10%-20%的晶粒尺寸范围内出现应变梯度，局部峰值应力达到平均值的3.1倍。晶界应力比晶内平均值高40%-100%，在较大晶粒中这种效应更明显。这支持了细晶粒可以在磨削过程中减少局部应力集中的观点。

图13. 晶粒内部应力的不均匀分布：(a) 等应力等高线图；(b) 等应力箱形图。

微观应力-应变不均匀性有两个主要来源：
(1) 相邻晶粒之间的取向不匹配导致不同的弹塑性响应。较大的取向错配角增加了晶界间兼容性的需求，从而提高了晶界应力集中。
(2) 不同晶粒中的滑移系激活程度不同。不同的滑移活动改变了局部应变适应，并通过晶界约束将应力传递给相邻晶粒。

4.1.2 微观结构演变
滑移系的激活是晶体塑性的主要载体。Schmid因子（SF）描述了滑移面、滑移方向和加载方向之间的几何关系[79]。在图14中，外力P在A1面上产生解析剪应力。其中λ表示滑移方向与力的夹角，φ表示滑移面法线与力的夹角，解析剪应力为：
(35)
其中σ是施加的应力，μsf是Schmid因子。μsf的符号表示优先的滑移方向，接近0.5的值对应于位错运动的较低障碍。

图14. Schmid因子示意图。

在基于位错的晶体塑性框架中，塑性变形由活性滑移系上的剪切描述。12Cr2Ni4A钢模型中有12个{110}<111>滑移系。当绝对解析剪应力|τα|超过其当前强度gα时，滑移系变得活跃。累积的剪切应变γα记录了每个活性系统上的塑性剪切。

图15显示UVG过程中的滑移激活在空间和时间上是异质的。磨削开始时，只有局部高应力下的晶粒激活有限的滑移系。随着载荷的移动和更多晶粒进入塑性流动，高Schmid因子系统首先激活，活性系统的数量增加。一旦材料达到准稳态，活性系统的数量几乎稳定。图16进一步显示，不同晶粒在磨削后仍保持不同的滑移激活模式。

图15. 不同时间下{110}<111>滑移系的有效剪切应变图。

图16. 不同滑移系内的剪切应变图。

图17(a)和(b)比较了UVG和CG中的CRSS和塑性应变场。具有有利Schmid因子的滑移系4和9在早期加载期间显示出更快的CRSS演变和更高的剪切应变积累。低Schmid因子系统较晚激活。这一序列表明滑移竞争受晶体取向和移动动态载荷的共同控制。一些系统在卸载过程中也变得不活跃，这有助于最终残余应力场的形成。

图17(c)显示滑移系4和9占总激活事件的约30%。它们的主导性与滑移面法线与加载方向之间的小角度有关。这些系统的不均匀激活产生了晶粒间的塑性不均匀性，随后表现为位错密度梯度和不均匀残余应力。

图17(d)比较了UVG和CG下12个滑移系的累积剪切应变。UVG使{110}<111>滑移系的激活增加了9.36%。这种增加归因于超声软化作用，降低了有效CRSS，使更多滑移系，尤其是高Schmid因子晶粒中的滑移系参与变形。

图17(e)显示在相同的加载步长下，UVG中{110}<111>滑移系的激活比CG早约0.18 ms。CG中的平均CRSS从477 MPa降至UVG中的453 MPa，磨削过程中的最大减少量为5.03%（图17(f)）。这种减少有助于解释UVG下更高的多滑移活性和更大的位错增殖率。

这些结果表明，超声振动通过降低CRSS阈值和激活更多滑移系来改善塑性应变适应。这种机制减少了局部应力集中，并有助于在加工表面以下更深处保留压缩残余应力。

图18(a)显示了UVG表面层的总位错密度分布。排除异常节点后，密度范围从7.881×10^7到1.894×10^9 mm^-2。晶粒41的密度约为1.565×10^9 mm^-2，而相邻晶粒53的密度约为8.483×10^8 mm^-2。相邻晶粒之间大约两倍的差异显示了取向和滑移活动对缺陷存储的影响。

图18. 同时滑移系激活与位错密度演变之间的相关性：(a) 表面层的总位错密度分布；(b) 晶粒41中的滑移系激活；(c) 晶粒53中的滑移系激活。

选择晶粒41和53进行滑移激活比较。在晶粒41中，四个滑移系被强烈激活，增加了交叉滑移和位错缠结的机会，形成了高密度的位错结构。在晶粒53中，只有两个滑移系被弱激活，因此位错相互作用较少，存储的密度较低。这一结果符合Taylor的多滑移概念：多个滑移系的激活增加了位错增殖和加工硬化。这也解释了为什么加工表面处的残余应力是不均匀的。与均匀连续介质模型不同，当前的热超声CPFEM模型展示了降低CRSS和增加多滑移激活如何将超声效应传递到位错密度和残余应力[15]。

4.1.3 UVG中残余应力的形成机制
UVG中的残余应力是在耦合的机械和热载荷作用下通过非均匀塑性变形形成的。在CPFEM结果中，残余应力分为宏观尺度上的I型应力、晶粒平均的II型应力以及亚晶粒或积分点尺度上的III型应力（图19）。I型应力受整体力和热平衡控制，II型应力取决于晶粒取向和滑移活动，III型应力集中在晶界和位错缠结区附近。例如，在图11中，<110>取向晶粒的残余应力约为450 MPa，而<100>取向晶粒的残余应力约为300 MPa；晶界处的III型应力可达到约±600 MPa。

图19. 三种类型的残余应力示意图。

图20比较了在相同宏观条件下UVG和CG预测的x方向残余应力σxx。两个面板使用相同的等高线窗口（-300至0 MPa，25 MPa间隔），以便直接比较压缩应力场。UVG产生的压缩应力幅度更大，有效压缩层更深，而两个场的应力都随深度逐渐减小。

图20. x方向的残余应力（两个面板使用统一的颜色条，等高线窗口为-300至0 MPa（Δ = 25 MPa）。此窗口外的值以饱和端颜色显示，图例中的min/max表示绝对极值）：(a) UVG；(b) CG。

图21显示了从模拟中提取的I型、II型和III型残余应力。在表面下8 μm处，UVG的最大压缩值为-369.13 MPa，CG为-263.66 MPa，UVG的压缩幅度大约大40.0%。在100 μm处，UVG的残余应力为-49.76 MPa，CG为-23.48 MPa，表明UVG在更深的深度保持更大的压缩应力状态。I型应力梯度遵循晶粒级别的II型应力波动，而III型应力在约±400 MPa范围内分散，表明应力从亚晶粒不均匀性传递到晶粒和宏观尺度。

图21. 沿深度方向的三种残余应力分布：(a) UVG中的三种残余应力分布；(b) CG中的三种残余应力分布；(c) II型残余应力比较；(d) I型残余应力比较

Taylor因子将单晶滑移阻力与多晶体的宏观屈服响应联系起来[80]，[81]。它反映了不同晶粒取向需要不同的滑移系组合来适应相同的宏观变形。泰勒因子Ms定义为：(36)其中σyield是宏观屈服应力，τCRSS是单晶的临界分解剪切应力，Ms是泰勒因子，对于体心立方（BCC）金属，其范围通常在2.0-2.8之间。这里使用泰勒因子作为应力传递路径的分析后处理描述符，而不是在UMAT中求解的内部变量。因此，我们并不认为仅通过直接计算Ms的减小就能确定压缩残余应力的增加。机械证据来自于模拟的CRSS减小、早期滑移激活、活跃滑移系统的数量增加以及在超声振动研磨（UVG）下更深的I型压缩残余应力分布。为了清晰起见，关于泰勒因子的讨论应被视为对模拟滑移响应的解释，而不是一个额外的拟合变量。因此，结论是基于直接计算的CRSS减小、早期滑移激活和更多的活跃滑移系统。在本研究中，泰勒因子用于描述III型局部应力如何传递到II型晶粒尺度应力。如果晶界处的局部III型应力为σIII，则相邻晶粒中的相关II型应力可以示意性地表示为：(37)当局部应力梯度足够大时，位错增殖和滑移激活会重新分配相邻晶粒中的II型应力。II型应力随后通过晶界兼容性贡献于I型残余应力场。定量上，平均CRSS从CG的477 MPa降低到UVG的453 MPa（图17(f)），相当于单晶滑移抗力的5.03%减小。在可比的宏观载荷下，这种减小表示适应相同表面变形所需的泰勒因子需求降低。因此，整个过程可以总结为：位错滑移产生III型应力梯度；泰勒因子介导的传递改变了II型晶粒应力；晶界适应修改了I型残余应力。因此，该机制可以表示为超声软化 -> CRSS/泰勒因子需求降低 -> 多滑移激活 -> 改善的残余应力保持。这一结论得到了模拟CRSS减小、滑移激活增加9.36%以及在UVG下更深的压缩应力分布的支持。这种解释仅基于经过验证的观察结果，并不假设在没有单独进行疲劳测试的情况下直接获得疲劳寿命的提高。

4.1.4. 工艺参数对12Cr2Ni4A材料机械响应和位错密度演变的影响
研磨参数：ap、ns、vw
图22(a)和(d)显示了研磨深度ap的影响。随着ap从5微米增加到15微米，平均压缩残余应力从-276.51 MPa降低到-240.41 MPa，而平均位错密度从1.40×10^8个/平方毫米增加到1.54×10^8个/平方毫米。图22(b)和(e)显示了砂轮速度ns的影响。将ns从500转/分钟增加到1500转/分钟，平均压缩应力从-291.23 MPa降低到-250.41 MPa，且在1500转/分钟时表面位错密度达到1.56×10^8个/平方毫米。图22(c)和(f)显示，将进给速度vw从500毫米/分钟增加到1500毫米/分钟，残余压缩应力从-228.22 MPa增加到-274.45 MPa，而位错密度从3.34×10^8个/平方毫米减少到2.84×10^8个/平方毫米。这些趋势反映了机械压缩、热软化、相互作用时间和位错增殖之间的竞争。

(2) 超声参数：f、A
下载：下载高分辨率图像（721KB）
下载：下载全尺寸图像
图22. 研磨参数的影响：(a) 研磨深度ap对平均残余应力的影响；(b) 研磨轮速度ns对平均残余应力的影响；(c) 进给速度vw对平均残余应力的影响；(d) 研磨深度ap对平均位错密度的影响；(e) 研磨轮速度ns对平均位错密度的影响；(f) 进给速度vw对平均位错密度的影响。
超声振幅和频率区分了UVG和CG。图23(a)显示，将振幅从5微米增加到20微米，平均表面压缩残余应力从-283.82 MPa变化到-295.51 MPa。超过15微米后，振幅的增加效果变小，斜率从1.22 MPa/微米降低到0.44 MPa/微米。这个值应被视为当前参数窗口内的模拟得出的拐点，而不是一个普遍的实验临界振幅。
下载：下载高分辨率图像（191KB）
下载：下载全尺寸图像
图23. 超声参数的影响：(a) 超声振幅A的影响；(b) 超声频率f的影响。
在工艺尺度上，超声振动降低了有效摩擦系数并改变了热输入。在晶粒尺度上，它降低了CRSS并增加了多个滑移系统同时变形的可能性。图23(b)显示，将频率从15 kHz改变到30 kHz对残余应力的影响远小于改变振幅的影响。因此，对于当前的机床系统和材料，振幅是更有用的超声控制参数。模型建议的实际振幅范围是10-15微米，但在生产使用前应结合表面粗糙度、研磨温度、工具磨损和振幅分辨的UVG实验进行验证。

(3) 晶体塑性参数：τ0、ρ0、k1、k2
晶体塑性参数也影响预测的机械响应。图24(a)显示，较高的初始CRSS τ0主要提高了屈服应力，对后期硬化斜率的影响较小。这是预期的，因为较大的τ0需要在晶格滑移开始前施加更高的应力。相比之下，初始位错密度ρ0同时影响屈服和硬化。
下载：下载高分辨率图像（193KB）
下载：下载全尺寸图像
图24. 晶体塑性参数对屈服阶段的影响：(a) 临界剪切应力τ0；(b) 初始位错密度ρ0。
如图24(b)所示，将ρ0从7.7×10^6个/平方毫米增加到2.7×10^8个/平方毫米可以提高初始屈服强度。在硬化过程中，低ρ0允许快速的位错增殖，而高ρ0从一开始就增加了位错相互作用和运动阻力。因此，钢材的初始缺陷状态影响了表层对UVG的响应。
图25(a)显示，将增殖系数k1从6.0×10^5增加到7.5×10^5加速了位错存储并提高了硬化速率。图25(b)显示了湮灭系数k2的相反趋势：较大的k2促进了位错湮灭，减少了有效位错密度，并减弱了稳定硬化过程中的应力增加。
下载：下载高分辨率图像（167KB）
下载：下载全尺寸图像
图25. 位错系数对硬化阶段的影响：(a) 位错增殖系数k1；(b) 位错湮灭系数k2
这些参数研究表明，UVG应与材料的初始状态一起优化。较大的ap和ns可能提高效率，但由于更强的热效应可能会降低压缩残余应力。较高的vw可以通过缩短热暴露时间来增加压缩应力。振幅的拐点意味着简单地增加振动幅度并不总是有效的。初始位错密度和τ0也表明，在选择研磨参数时应考虑锻造和热处理历史。

4.2. 计算模型与实验结果的比较和验证
4.2.1. EBSD测试结果的分析
EBSD用于比较UVG和CG样品的近表面微观结构。图26显示了核平均错位（KAM）图。UVG的平均KAM值为0.58°，CG为0.53°，增加了9.4%。由于KAM反映了局部晶格畸变，并与几何必要位错（GND）密度相关，结果表明UVG后位错存储更强。
下载：下载高分辨率图像（2MB）
下载：下载全尺寸图像
图26. KAM图和统计：(a) UVG；(b) CG。
图27中的晶粒尺寸统计显示UVG后晶粒细化更为明显。平均等效晶粒直径从CG的1.24微米减小到UVG的0.70微米，面积加权平均值从3.21微米减小到1.50微米。这种细化与模型预测一致，即超声振动降低了CRSS并增加了多滑移活动，从而促进了表层中的位错积累和边界重排。
下载：下载高分辨率图像（239KB）
下载：下载全尺寸图像
图27. 晶粒等效直径的统计分析：(a) CG；(b) UVG。
图28给出了GND密度。预研磨状态的平均GND密度为0.37×10^14个/平方毫米。研磨后，UVG的GND密度为0.48×10^14个/平方毫米，CG为0.39×10^14个/平方毫米。因此，UVG使GND密度相对于预研磨状态增加了29.73%，相对于CG增加了23.08%。UVG下更均匀的表面GND分布与模拟的从局部滑移到协调多滑移变形的转变一致。
下载：下载高分辨率图像（4MB）
下载：下载全尺寸图像
图28. 几何必要位错（GND）密度图和统计：(a) 研磨前；(b) UVG；(c) CG。
图29比较了极图和逆极图。极图中的最大均匀密度倍数（MUD）从CG的5.39减小到UVG的3.16，减少了41.47%。在逆极图中，最大MUD从2.17减小到1.96。UVG后较弱的纹理表明取向响应更加分散，支持了振动促进多个滑移系统激活而非单一优选取向主导变形的模型结果。
下载：下载高分辨率图像（2MB）
下载：下载全尺寸图像
图29. UVG和CG的极图和逆极图的比较：(a) UVG极图；(b) CG极图；(c) UVG逆极图；(d) CG逆极图。
观察到的较高压缩残余应力、较高位错密度、晶粒细化和较弱的纹理都是表面完整性的有利指标。它们表明UVG不仅可以作为精加工过程，还可以作为齿轮的表面状态控制过程[8]。关于接触疲劳寿命、耐磨性或弯曲疲劳的直接声明仍需通过专门的服务性能测试来验证，因此目前的结果最好解释为潜在性能改进的机械证据。

4.2.2. 研磨残余应力实验结果的分析
通过两种方式将模型与实验结果进行比较和验证：不同工艺参数下的表面残余应力以及沿深度方向的残余应力。图30显示，预测和测量的表面应力趋势一致，平均误差为21.13%。尽管这个误差高于某些为特定情况校准的连续体残余应力模型，但CPFEM模型提供了关于滑移竞争、位错密度梯度和III型到II型应力传递的额外信息。因此，该模型主要用于解释机制和趋势，而不是在没有特定情况校准的情况下预测绝对残余应力值。
下载：下载高分辨率图像（310KB）
下载：下载全尺寸图像
图30. 不同参数下表面残余应力的验证：(a) ap；(b) ns；(c) vw。
模拟与测量之间的剩余差异对于当前的建模策略是合理的，但不应忽视。XRD平均的应力覆盖的区域远大于CPFEM的晶粒尺度域，且等效移动载荷无法再现每个磨料颗粒的随机接触。这些差异主要影响绝对应力水平，而预测的工艺趋势与测量结果一致。
当ap从5微米增加到15微米时（图30(a)），测量的表面压缩残余应力从-331.54 MPa降低到-298.46 MPa。当ns从1000转/分钟增加到2000转/分钟时（图30(b)），测量的压缩应力从-339.73 MPa降低到更低的值，趋势与模型相同。相比之下，将vw从500毫米/分钟增加到1500毫米/分钟时，压缩应力从-242.18 MPa增加到-316.97 MPa（图30(c)）。这些趋势反映了研磨过程中机械压缩、热输入和位错恢复的竞争作用。
图31比较了在ap = 10微米、ns = 1500转/分钟、vw = 1000毫米/分钟、A = 5微米和f = 20 kHz条件下UVG和CG的深度剖面。模拟和实验都显示，压缩残余应力在表面附近最高，并随深度减小。测量的UVG表面值为-328.37 MPa，而模型值为-248.14 MPa。深度上的衰减归因于塑性变形和热流从接触表面的衰减。
下载：下载高分辨率图像（213KB）
下载：下载全尺寸图像
图31. 沿层深度方向的残余应力验证：(a) UVG；(b) CG。
在前10微米内，测量和模拟残余应力之间的平均差异为22.3%。在100微米的深度剖面上，UVG预测的最大绝对误差为102.79 MPa，平均绝对误差为78.69 MPa，而CG的相应误差分别为122.14 MPa和85.48 MPa。在大约30微米以下，误差变得更大。主要原因是实际材料中的随机纹理和预应力变化、XRD对多个晶粒的体积平均处理、电抛光过程中可能发生的应力松弛、简化的等效超声移动载荷，以及使用从单轴压缩而不是多轴磨削特定载荷校准的参数。比较结果证实，热超声CPFEM能够捕捉到主要的残余应力趋势和观察到的微观结构变化。UVG测得的表面压缩残余应力为-328.37 MPa，比CG在验证条件下的值高出约22.8%，而模拟结果显示在100 μm深度处的压缩应力保留得更好。这一机制与通过多滑移激活将III型局部应力转化为II型晶粒尺度应力的过程一致（见图19）。该模型仅在测试的磨削深度、砂轮速度、进给速度、振幅和20 kHz激励范围内进行了验证。对于其他钢材或UVG系统，需要重新确定晶体塑性参数Usoft和Tsoft以及热摩擦边界条件。在测试范围之外的预测应视为定性结果，直到有新的验证数据可用。因此，这种本构校准应被视为当前UVG条件的有效参数化方法。未来的改进应包括高应变率或逆向磨削特定识别，以便更直接地表示多轴载荷和快速热循环。

总体而言，UVG-CPFEM框架通过识别磨削参数和超声输入如何影响滑移活动、位错密度和保留的残余应力，减少了依赖于试错实验的需求。其主要价值在于机制性：它展示了哪些微观结构路径能够提高12Cr2Ni4A齿轮钢的表面完整性。

5. 结论
本研究开发了一个用于12Cr2Ni4A齿轮钢热超声CPFEM的框架，并通过热压缩、EBSD和XRD数据进行了验证。主要结论如下：
(1) 基于位错密度的CPFEM通过将热软化和超声软化耦合到滑移律中得到了扩展。校准后的模型再现了热压缩应力-应变曲线，平均相对误差为1.8%，为后续的磨削模拟提供了材料参数基础。
(2) 与CG相比，UVG将有效CRSS降低了5.03%，并将{110}<111>滑移系统的激活提高了9.36%。这促进了协调的多滑移变形，并支持了从III型局部应力到II型晶粒应力，再到I型残余应力的应力传递。在CPFEM结果中，峰值压缩残余应力达到了-369.13 MPa，而在100 μm深度处的残余应力在UVG下更为压缩（-49.76 MPa vs. -23.48 MPa）。
(3) 工艺参数通过耦合的机械和热效应影响残余应力和位错密度。在模拟范围内，将磨削深度从5 μm增加到15 μm，平均压缩残余应力从-276.51 MPa变为-240.41 MPa，平均位错密度从1.40×10^8 mm^-2增加到1.54×10^8 mm^-2。提高砂轮速度会抑制位错密度的增长，因为热软化与加工硬化相互竞争。超声振幅显示出拐点效应，在当前模拟窗口内，当振幅超过约15 μm时，边际效益会降低。
(4) EBSD和XRD测量结果支持了预测的趋势。与CG相比，UVG使GND密度增加了23.08%，并产生了更细小的晶粒和更弱的纹理。XRD残余应力测量显示，在表面附近的平均相对误差为22.3%，沿100 μm UVG深度剖面的平均绝对误差为78.69 MPa。该模型比绝对应力幅值更可靠地捕捉到了趋势，因为绝对应力幅值受到测量平均、实际纹理和简化等效载荷的影响。
结果表明，超声振动通过改变晶体滑移和位错存储来改变表面完整性，而不仅仅是改变宏观磨削力。因此，该框架将加工、缺陷演变和残余应力联系起来，这对于高强度钢表面制造具有重要意义。除了这里研究的齿轮磨削案例外，只有在重新校准材料常数和载荷边界条件后，同样的工艺-结构-应力逻辑也可以用于其他机械和热驱动的表面处理过程。
该模型的局限性在于其校准范围和等效边界载荷处理方法。将其扩展到其他材料、振幅、频率或冷却条件将需要新的参数识别和验证。另一个局限性是C11、C12、C44和一些硬化参数被视为温度独立的。包括温度依赖的弹性常数和硬化参数应该能够更好地预测在陡峭磨削热梯度下的弹性回弹和最终残余应力。这一局限性对绝对应力值的影响大于本研究中识别的III型->II型->I型机制。

作者贡献
所有作者都参与了研究的构思和设计。材料制备、数据收集和分析由邹润祥、唐金元和周卫华完成。手稿的第一稿由邹润祥撰写，所有作者都对之前的手稿版本进行了评论。所有作者都阅读并批准了最终手稿。

利益冲突声明
作者声明他们没有已知的竞争性财务利益或个人关系，这些关系可能会影响本文报告的工作。