0 前言

喷丸强化技术是一种金属材料表面处理工艺，喷丸设备通过加速大量微小弹丸冲击材料表面，诱导材料表面产生塑性变形和残余应力，产生加工硬化效应^[1]。金属材料在实际服役过程中发生频率最高、危害性最大的断裂模式是疲劳断裂，喷丸强化工艺主要用于改善材料的疲劳性能，广泛应用于航空航天、车辆船舶等领域^[2]。喷丸强化机理主要是残余应力强化机制和组织结构强化机制^[2]，前者是由表层发生不均匀弹塑性变形产生残余应力，残余压应力可以抑制和阻碍表面裂纹的萌生和扩展，后者是由弹丸高速冲击使材料表层发生循环塑性变形，造成材料表层显微组织结构发生变化^[2-3]。 POUR-ALI，KUMAR 等^[4-5]在喷丸处理的材料表面观察到位错滑移和晶粒细化，利用剧烈喷丸技术制备了梯度分布的纳米晶表面层。KUMAR 等^[5]研究发现，表面晶粒细化为纳米级会改变材料的断裂力学行为，从而改善材料的疲劳断裂、微动疲劳、磨损和腐蚀等性能。因此许多学者通过试验和模拟进行了喷丸引起晶粒细化的研究^[6-8]。

已有研究表明，喷丸处理诱导材料表层晶粒细化的效果与喷丸材料、大小、速度、角度、覆盖率等工艺参数有关^[9-11]。与通过试验方法调整喷丸工艺从而获得所需的纳米级表面层^[8]相比，数值模拟方法可以大大节省成本和时间，便于研究人员优化工艺参数。HASSANI-GANGARAJ 等^[12]研究了喷丸过程中晶粒细化机制，建立了一个多尺度数值模型模拟 AISI4340 钢在喷丸过程中的微观组织演变规律，该数值模型预测的 AISI4340 钢位错密度演化过程和晶粒细化结果与试验相吻合。WANG 等^[11]进行了随机多喷丸冲击的参数化建模，将多尺度本构模型与有限元模拟相结合，详细研究了喷丸工艺参数对 AISI4340 钢晶粒细化的影响，当喷丸尺寸、速度和覆盖率增大时，晶粒细化程度增加。

7050 铝合金被广泛用作航天航空领域大型构件，相比于传统铝合金具有较好的强度和抗应力腐蚀性能^[6]，目前对于 7050 铝合金喷丸处理后组织形貌的研究主要是试验方法。宋颖刚等^[6]采用喷丸技术在 7050 铝合金表面制备出纳米结构表面层，利用 XRD、TEM、EBSD 等技术研究了喷丸后 7050 铝合金的表层纳米结构，研究表明表面纳米化机理是由于表面局部强烈塑性变形导致位错聚集。李永辉等^[13]通过增压喷丸技术在 7050 铝合金表面制备了纳米晶层，研究了喷丸强化对 7050 铝合金微观组织的影响，表层纳米化可以提高表层显微硬度和耐磨性能。但对于喷丸过程中 7050 铝合金晶粒细化与位错密度演化的关系鲜有研究，且关于喷丸工艺参数与 7050 铝合金纳米化表层之间的联系不够明确。本文采用有限元模拟和算法相结合的方法，建立一个将喷丸冲击的宏观三维有限元模型和微观位错密度演化模型相结合的多尺度数值模型，用以研究喷丸过程中 7050 铝合金晶粒细化与工艺参数之间的关系，预测将 7050 铝合金表层晶粒细化为纳米尺度对应的喷丸工艺参数，进一步揭示 7050 铝合金喷丸晶粒细化机理，指导 7050 铝合金喷丸表面纳米化工艺。

1 数值模拟方法

利用 Abaqus 软件分别建立单个和多个喷丸冲击的三维有限元模型，并将有限元模拟与位错密度演化联系起来，以预测喷丸冲击诱导的晶粒细化。建立单个喷丸冲击模型时靶材属性采用 JC 本构模型，输出单元的应变和应变率等，代入 Matlab 求解位错密度本构模型方程组，利用遗传算法求得两种本构模型的流动应力差值最小的位错密度本构模型的参数集。建立多个喷丸冲击模型时利用 VUMAT 用户材料子程序作为靶材属性，采用上述位错密度本构模型参数集。研究喷丸工艺参数对晶粒细化的影响，为优化工艺参数制备纳米级表层提供参考依据。参数集优化及求解流程图如图1 所示。

1.1 喷丸有限元模型

模拟所用 7050 铝合金化学成分（质量分数）见表1^[14]。7050 铝合金的固溶热处理温度为 471～482℃，在冷水中淬火，转移时间不大于 15 s，时效处理制度为加热温度 116～126℃，保温 3～5 h，加热温度 158～168℃，保温 24～26 h^[15]。喷丸采用常规常温喷丸强化。

1.1.1 单个喷丸冲击

采用 Abaqus 有限元软件，建立单个喷丸冲击的三维有限元模型如图2 所示。待喷靶材尺寸为 3 mm ×3 mm×1.5 mm 的立方体，靶材底面约束所有自由度，喷丸为半径 0.3 mm 的球体。采用 C3D8R 八节点线性六面体减缩积分和沙漏控制的单元。中心区域加密网格，受喷区域网格尺寸与凹坑直径之比约 1∶20^[16]，最小网格尺寸为 14 μm。使用罚算法定义接触属性，库仑摩擦因数为 0.2。在预定义场中定义喷丸速度为 65 m / s，采用动力显式分析，结果输出中包含时间增量、应力、塑性应变及应变率等。弹丸定义为弹塑性材料，密度 7 850 kg / m³ ，屈服应力 1.55 GPa，弹性模量 210 GPa，泊松比 0.3。为了简单起见，将喷丸撞击靶材的角度设置为 90°，因为垂直撞击是产生均匀残余应力的推荐入射角度，偏离垂直冲击可能诱导局部多向塑性变形^[12]。

喷丸冲击过程中，应变率决定了残余应力大小、分布以及塑性变形区状态^[17-18]。塑性变形时应力与应变的关系称为本构方程。JOHNSON 和 COOK^[19]提出一个流动应力的计算模型，其中包括应变硬化、应变速率硬化和热软化的影响，JC 本构方程如下：

式中，σ 是材料的流动应力，ε 是材料变形时的等效塑性应变，$\dot{\epsilon }$ 和 ${\dot{\epsilon }}_0$ 分别表示材料的应变率和参考应变率，T、T_room、T_melt 分别为材料发生变形时的温度、参考温度和熔化温度，A 为材料初始屈服应力， B 和 n 是描述硬化系数和指数的常数，C 和 m 分别是描述流动应力对应变速率和温度敏感性的常数。

采用 JC 本构模型定义靶材的属性。由于喷丸强化过程中温度几乎不发生变化，所以忽略了热效应的影响。FU 等^[14]研究了 7050 铝合金高应变率冲击压缩中流动应力与应变、应变率的关系，得到 7050 铝合金JC方程参数如表2所示。JC本构方程可使用Abaqus 软件自带的材料属性模块设置并求解出来。

Note: A is the initial yield stress of the material; B is the hardening modulus of the material; n is the hardening index of the material; C is the constant describing the sensitivity of the flow stress to the strain rate; m is the constant describing the sensitivity of the flow stress to temperature; ${\dot{\epsilon }}_0$ is the reference strain rate; T_melt is the melting point temperature of the material

1.1.2 多个喷丸冲击

建立规则排列的弹丸冲击和大量随机弹丸冲击的有限元模型。规则排列的弹丸冲击由 1.1.1 节中单个弹丸阵列装配而成，模型如图3 所示，用以研究靶材同一位置受力变形和晶粒细化随弹丸撞击次数的变化，以及邻近弹丸冲击的影响。

为更真实地模拟喷丸强化过程，建立一个基于概率控制多弹丸冲击靶面位置的三维喷丸有限元模型。为消除剧烈喷丸时应力波在靶材边界产生反射回弹影响结果精度，在靶材四周和底部装配一层无限元 CIN3D8 单元，使得模型具有无限边界^[20-21]，如图4 所示。

喷丸覆盖率是喷丸区域与总区域的比值，以百分比表示，用 Avrami 方程来评估喷丸覆盖率^[22]：

式中，C 表示喷丸覆盖率，r 表示喷丸凹坑的半径， R 为每个单位区域凹坑的生成率，t 表示喷丸时间。一般将 98%喷丸覆盖率视为 100%完全覆盖，200% 喷丸覆盖率定义为达到 100%喷丸覆盖率所需时间的两倍，并依次类推。

在实际喷丸过程中，弹丸在空间上并不是完全随机分布的，因此，使用概率控制弹丸随机分布，即弹丸之间彼此的时空位置不能重叠，弹丸分布函数如式（3）所示^[23]：

式中，P 为凹坑生成的概率，l 为随机生成的凹坑中心与之前凹坑中心的距离。如果 l＝0，即新生成的凹坑与之前的凹坑完全重叠，则概率为零（P＝0），即禁止新生成凹坑。如果 l＞2r，即新生成的凹坑完全超出之前的凹坑，则概率为 1（P＝1），允许生成新的凹坑。否则，如果为 0＜l≤2r，新的凹坑产生的概率由分布函数式（3）控制。

将上述随机概率控制过程使用 Python 脚本对 Abaqus 进行二次开发可完成模型创建，具体开发过程见文献^[23-24]。由于弹丸数量多，在脚本中增加了 Abaqus 中弹丸与靶材相互作用时录制的宏，将弹丸数量、半径、速度、冲击区域范围等写为用户输入数据，可完成自动化模拟，便于更改工艺参数时一键操作完成整个模拟。

1.2 位错密度演化模型

1.2.1 理论模型

金属的塑性变形依赖于应变率，它被普遍认为是与金属的晶格缺陷，主要是位错的产生、运动、相互作用和累积有关，所以位错与各种障碍及其自身的相互作用影响着金属的流动应力^[25]。金属的塑性变形程度依赖于应变率和应变速率，ESTRIN 等^[26] 研究了多种应力水平下的晶体的应变硬化过程，提出一个位错密度模型，由两个微分方程耦合组成，用于预测晶胞内和晶胞壁位错密度的演化。 HASSANI-GANGARAJ 等^[12]对上述位错密度模型进行改进，将位错演化与材料的流动应力结合起来，建立位错密度本构模型中初始参数优化的方法，并通过 TEM 验证该模型预测高应变率硬化下的位错密度和晶粒尺寸的可靠性。采用该模型对喷丸过程中的晶粒细化进行评价。

在该模型中，位错胞结构包括胞壁和胞内两部分。胞壁 ρ_w 和胞内 ρ_c 的位错密度的演化遵循不同的规则^[26]：

式（4）、（5）右边第一项表示由 Frank-Read 源的激活而引起的位错产生，第二项表示胞内的位错转移到胞壁，最后一项表示塑性变形时位错的湮灭导致位错密度动态平衡。参数 α^*、 β^*、k_c、k_w 为位错形成和湮灭相关的常数，n_c、n_w 为温度相关的参数。由于喷丸是一种冷处理过程，因此 n_c、n_w 假定为常数。参数 b 表示伯氏矢量的模长， r^γ 是分解的剪切应变，${\dot{\gamma }}_0$ 为分解的参考剪切应变率，d 为晶粒大小，f 为位错胞壁的体积分数。考虑了位错壁的体积分数，总位错密度可以表示为胞壁和胞内位错的加权和^[10，26]：

参数 f₀ 和 f_∞ 分别为 f 的初始值和饱和值。$\tilde{\gamma }$ 是一个衡量 f 下降率的常数。平均晶粒尺寸与总位错的关系表示为^[26]：

式中，K 是一个取决于累积塑性应变的参数，随总位错密度的增加而迅速减少^[26]：

式中，K₀ 和 K_∞ 分别为 K 的初始状态值和饱和值。

为了将上述位错相互作用与材料的力学行为联系起来，将材料的流动应力表示为两项之和^[12]：

式中，σ₁ 表示一个与应变无关的应力，该应力源于位错滑移阻力，与位错间相互作用无关。该应力可根据变形前材料的屈服应力来估算。事实上，这个 σ₁ 与式（1）中的 A 起着相同的作用，下一节将详细讨论式（10）与式（1）的一致性，介绍位错密度演化模型中参数优化的方法^[12]。σ₂ 来源于位错间相互作用，且与应变和应变率相关。通过泰勒系数，分解剪切应力和剪切应变率可与σ₂ 和等效应变率相关^[26]：

式中，分解剪应力τ^r 描述位错胞结构的力学行为，与胞壁和胞内的剪应力有关。剪应力是通过应用混合规则得到的^[26-27]：

将这些分解剪应力与其相应的位错密度联系起来，就形成一组方程^[26-27]：

式中，α为常数，G 为剪切模量，1 / m^* 为应变速率敏感性指数。

因此，该模型可以将塑性应变与位错演化联系起来，并考虑塑性变形时应变、应变率和位错密度对材料流动应力的贡献。

1.2.2 最佳输入参数集

式（4）～（15）构成包含隐式微分方程的复杂非线性方程组，使用 Matlab 软件可求解。将式（4）～（15）写为 Matlab 函数程序，求解大致过程如下：首先利用 1.1.1 节所建模型从 Abaqus 软件中提取到单元的应变和应变率等随时间变化，提取出 ${\dot{\gamma }}_0$ 约为 0.35 ×10⁶， r^γ 为 0.03，$\dot{\mathcal{E}}$为 0.2×10⁶，ε 为 0.07；然后利用Matlab软件四阶五级Runge-Kutta-Felhberg算法^[28-29] 求出 ρ_c、ρ_w、ρ_t 和 d 的数值解，该数值解为 n 行 4 列的矩阵（n 与定义的求解步长总时间有关）；利用 1.1.1 节所建模型从 Abaqus 软件中提取单元塑性变形所用时长约为 1.6μs，定位该总时长所在的该行矩阵为所要求解的ρ_c、ρ_w、ρ_t 和 d 的终值。

b、M、α、G、f₀、 f_∞、 ${\dot{\gamma }}_0$、K₀、 K_∞、$\tilde{\gamma }$、β、ρ_c^t=0、ρ_w^t=0 作为材料的初始参数值；α^*、β^*、k_c、 k_w、n_c、n_w、m^* 作为使有限元 JC 本构求得的流动应力值（式（1））与位错密度演化本构求得的流动应力值（式（10））相一致的参数集^[12]。将两种本构得出的流动应力值的差值的绝对值作为目标函数，使用Matlab软件遗传算法进行全局最小化优化来实现最佳的输入参数集，最终得出的最佳参数集见表3 和表4。综合考虑运算效率与精确性，采用 100 次遗传迭代计算最优解，该参数集最终使得流动应力的差值绝对值约 2×10⁻²⁵ MPa，如图5 所示为每代最佳适应度值变化图。

Note: α^* is the proportion of moving dislocation sources in the total dislocation sources, β^* is the proportion of intracellular dislocations transferred to the cell wall when they occur, k_c and k_w are the relevant parameters of energy loss rate associated with the dynamic process of dislocation cross slip in the cell and cell wall during dislocation annihilation, n_c and n_w are the relevant parameters of the sensitivity of dislocation stack fault energy values in the cell and cell wall to temperature, and m^* is the reciprocal of the strain rate sensitivity index.

Note: b is the size of the Bernstein vector, M is the Taylor factor, α is the empirical constant that relates the stress to the square root of the dislocation density, G is the shear modulus, f_o is the initial value of the cell wall volume fraction f before the start of the shear strain, f_∞ is the saturation value that the cell wall volume fraction f decreases to when the shear strain increases, $\mathrm{}{\dot{\gamma }}_0$ is the reference decomposition shear strain rate, use K to represent the proportional coefficient that relates the dislocation density to the cumulative plastic strain, K₀ is the initial value of the proportional coefficient K before the onset of plastic strain, K_∞ is the saturation value reached when the proportional coefficient K decreases with the development of plastic strain, $\stackrel{-}{\gamma }$ is the constant representing the decline rate of f, β is the best fitting parameter in the fitting function relationship between K and decomposed shear strain, andρ_c^t=0, ρ_w^t=0 are the initial dislocation densities in cells and cell walls before deformation.

1.2.3 位错密度演化在靶材中的分布

剧烈喷丸冲击过程中应变率时程提取计算较难，为了直观地在 Abaqus 中显示出不同深度处晶粒尺寸和位错密度随冲击过程变化的结果云图，开发用户材料子程序 VUMAT，将参数优化后的位错密度演化本构方程嵌入到 Abaqus 显式分析中。

2 结果与讨论

2.1 单个喷丸冲击的晶粒细化

半径为 0.3 mm 的弹丸以 65 m / s 速度垂直单次冲击 7050 铝合金引起的残余应力、等效塑性应变、总位错密度和晶粒尺寸的分布如图6 所示。根据所选的初始位错密度参数，计算出初始晶胞尺寸约为 12.2 μm。单次冲击后晶粒细化为约 2.6 μm，位错密度增加约两个数量级。最大等效塑性应变、最大位错密度、最小晶粒约在同一位置，即弹丸冲击产生的凹坑表面下方约 40 μm 处。由此可见，单次冲击产生了显著的晶粒细化，但单次冲击不会产生均匀的纳米结构层。

四阶五级 RKF 算法求解位错密度本构方程组时，求解出不定步长时间对应的晶粒尺寸，图7 显示了单个弹丸 65 m / s 冲击速度时最小晶粒尺寸随时间的演化，可以看出晶粒细化随碰撞时间不断发展。图8 显示了单个弹丸冲击后的最小晶粒尺寸与速度关系，可以看出晶粒细化程度随着速度的增加而增加。

2.2 多个喷丸冲击的晶粒细化

图9显示了撞击同一位置的弹丸个数n、尺寸 r、速度 v 与最小晶粒尺寸 d 的关系，可以看出，撞击次数越多，速度越大，弹丸尺寸越大，晶粒细化越显著，且速度是对晶粒细化影响最显著的因素。当弹丸如图3 呈 4 列每列 15 个撞击时，半径 0.3 mm，速度 80 m / s 时得到最小晶粒尺寸 493 nm，影响区平均晶粒尺寸约 1.4 μm，而相同尺寸和速度的弹丸呈单列共 15 个撞击时，获得 533 nm 的最小晶粒尺寸和 1.6 μm 的平均晶粒尺寸，说明邻近弹丸冲击的影响会提高晶粒细化作用，且邻近弹丸冲击的作用使得晶粒细化影响区域宽度变大，如图10 所示。如图11 所示，当 4 列每列 15 个弹丸的撞击速度增大到 100 m / s 时，靶材表面平均晶粒尺寸进入纳米级，约为 290 nm，晶粒细化明显的区域深度约为 1.3 mm，表层 300 μm 深度范围内的平均晶粒尺寸约为 0.6 μm，由此可见，喷丸速度增大会加深晶粒细化的范围。

弹丸的反复冲击会提高晶粒细化程度，在随机弹丸冲击模型中，喷丸覆盖率会对晶粒细化产生影响，图12 显示了喷丸半径 0.3 mm，速度 50 m / s 时，不同覆盖率下的最小晶粒尺寸，随着喷丸覆盖率增大，最小晶粒尺寸减小。

如图13 所示，弹丸半径 0.3 mm，速度 80 m / s，当喷丸覆盖率在 100%时，晶粒细化影响深度约为 50 μm，凹坑表面晶粒尺寸虽然细化为纳米级，但影响深度内平均晶粒尺寸仍为微米级，没有形成均匀纳米结构层。当喷丸覆盖率在 200%，速度 50 m / s 时，晶粒细化影响深度约为 100 μm，影响深度内平均晶粒尺寸仍为微米级。

当喷丸速度为 70 m / s，覆盖率在 300%以上时，喷丸表面和次表面的平均晶粒尺寸细化为纳米级，纳米级结构表层深度约 70 μm。如图14 所示，当喷丸速度为 50 m / s，覆盖率在 500%以上时，喷丸表面和次表面的平均晶粒尺寸细化为纳米级，形成均匀的纳米表层，纳米级表层深度约 150 μm。沿靶材深度方向提取路径，统计晶粒尺寸梯度分布，如图15 所示，随着表面到基体深度的增加，晶粒尺寸逐渐增大。这与宋颖刚等^[6]通过喷丸工艺制备梯度纳米结构表面层的结果相符合。

3 结论

（1）JC 本构方程与位错密度本构方程均能描述高应变率材料的塑性应力应变关系，基于两者对流动应力预测的一致性，可实现位错密度本构方程的优化调参。

（2）使用基于位错密度理论建立的喷丸强化 7050 铝合金的混合数值模型，可直观显示出靶材的动态力学行为（应力、应变率等）变化及微观结构 （位错密度、晶粒尺寸等）演变。位错演化机理与塑性应变过程中流动应力变化有关，位错演化与塑性应变相互适应。

（3）喷丸冲击造成靶材的高应变率塑性变形，诱导靶材表面和次表面产生了晶粒细化和位错密度增加。改变喷丸工艺参数会改变晶粒尺寸和分布，较高的喷丸强度和覆盖率可得到梯度纳米结构层。

参考文献