Study on numerical simulation of residual stresses induced by laser shock processing

激光冲击强化技术又称激光喷丸，与传统机械喷丸技术相比，它能在工件表面形成更深的残余压应力层^[1]，并且可控性强、适应性好，能够处理倒角、沟槽、小孔等难以处理的部位。基于这些优点，该技术已应用在航空发动机叶片、齿轮、核电站压力焊缝等众多领域^[2]。

数值模拟有助于获得特定应用场合最优的工艺参量，已逐渐成为研究激光冲击强化的重要手段。国内外学者针对激光冲击强化的建模及其优化做了大量的研究工作。BRAISTED和BROCKMAN采用ABAQUS软件建立了2维轴对称模型，对激光冲击处理形成的残余应力场进行了预测^[3]。NAM采用独立编程对激光冲击Ti-6Al-4V厚板及薄板诱导产生的残余应力场进行了模拟分析，模拟结果与实验结果具有较好的一致性^[4]。HU等人采用ANSYS/LS-DYNA软件对激光冲击过程中靶材的残余应力场进行了模拟研究^[5]。BHAMARE等人仅采用LS-DYNA软件进行动态应力计算，通过重复单次冲击模拟计算过程来获得多点搭接的残余应力场^[6]。ACHINTHA和NOWELL采用基于本征应变的模型来研究工件尺寸及冲击次数对残余应力场的影响，提高了模拟计算效率^[7]。JIANG等人采用ABAQUS软件模拟激光冲击铝合金小孔构件，研究了冲击区域范围对试样表面及孔壁残余应力场的影响^[8]。

以上模拟均是建立在单次激光冲击模型的基础上，为了获得更为准确的残余应力场模拟结果，建立更为有效的单次激光冲击强化模型显得尤为重要^[9]。但是从已发表的文献来看，不同研究者采用的模型网格单元尺寸及激光冲击波压力空间分布模型不同，得到的模拟结果也有一些差异，而针对模型网格单元尺寸及激光冲击波压力空间分布模型对残余应力场模拟结果影响的相关研究鲜有报道。

本文中建立了单次激光冲击35CD4合金钢厚板有限元模型，对激光冲击过程中靶材内部的能量变化、表面动态应力进行了分析，讨论了模型网格单元尺寸、冲击压力空间分布模型选取对残余应力场模拟结果的影响, 分析了单次激光冲击下材料的残余应力场的分布特点，并与已知的一些实验测量结果进行了对比。

采用ABAQUS进行激光冲击强化模拟主要分为两个步骤：(1)采用ABAQUS/explicit模块做显式动态分析计算，模拟靶材的动态响应过程；(2)将显式动态计算的应力场导入ABAQUS/standard模块，进行隐式分析计算，以确定最终稳定的残余应力场。如需进行多次激光冲击模拟，还需将以上的结果回传到ABAQUS/explicit中，作为后一次冲击的预定义场。

本文中模拟单面激光冲击35CD4合金钢厚板，载荷压力分布对称，载荷作用半径远小于靶材尺寸，因此为了提高计算效率，可以将3维几何模型简化为一个2维轴对称模型，如图 1所示。其中r表示到激光光斑中心的距离，激光冲击载荷作用半径为4mm，细网格区域6mm×6mm，采用四节点轴对称有限单元CAX4R模拟冲击区域的残余应力场。周边网格采用四节点轴对称无限单元CINAX4，该种单元为弹性单元，作为有限单元区域的无反射边界，以消除弹塑性波的反射。

Figure 1.  One-sided LSP model with axisymmetric boundary conditions

PEYRE等人^[10]对约束模型下的激光冲击波压力求解进行了半理论研究，其给出的单次激光冲击峰值压力p的估算公式如下：

式中, α为相互作用系数，一般取为0.1~0.2；Z₁和Z₂分别为目标靶材和约束层的声阻抗，单位为g·cm^-2·s^-1；I₀为激光功率密度，单位为GW·cm^-2，可以由公式I₀=E_l/(πr_p²τ)计算得到，其中E_l为激光器输出能量，r_p为光斑半径，τ为激光脉冲宽度。计算得到的冲击峰值压力p的单位为GPa。

根据参考文献[11]，Z₁=3.6×10⁶g·cm^-2·s^-1，Z₂=1.65×10⁵g·cm^-2·s^-1, α取为0.1, I₀取为10GW/cm²，计算得到的峰值压力为3.14GPa，而实验测得的峰值压力为2.8GPa^[12]，说明该模型估算的峰值压力比实验测得的峰值压力要大。

一般实验中测得的激光冲击载荷压力-时间分布为准高斯分布^[12]，在约束模式下, 冲击波压力的持续时间超过激光脉宽的2~3倍。BRAISTED等人^[3]采用三角形分布来近似准高斯分布。NAM^[4]基于实验给出了时间及空间上的载荷压力分布模型。

然而冲击载荷在空间上的分布很难通过实验来研究，模拟中一般采用均匀分布模型^[10]。ZHANG等人^[13]研究了激光冲击波压力在光斑径向上的分布，认为激光冲击波压力在空间服从高斯分布。本文中采用三角形瞬态分布模型(如图 2所示)，半峰全宽(full width at half maximum，FWHM)为50ns。在空间分布上采用图 3所示的3种分布模型^{[4, 13]}，其中p_max为最大峰值压力。

Figure 2.  Relationship between pressure and time of laser shock pressure

Figure 3.  Three kinds of spatial distribution models^{[4, 13]}

激光冲击材料产生的应变率超过10⁶s^-1，常用的描述材料在高应变率条件下的本构方程有理想弹塑性(elastic perfectly plastic，EPP)模型、Johnson-Cook(JC)模型、Khan-Huang-Liang(KHL)模型^[14]和Zerilli-Armstrong(ZA)模型^[15]等。

目前JC模型是在激光冲击强化数值模拟中使用最广泛的一种基于实验得到的材料本构模型，它考虑了应变速率效应和温度效应，忽略了冲击波压力和变形历史的影响，一般更适用于应变率小于的情形。KHL模型考虑了材料在高应变率下的硬化率下降的因素，是对JC模型的修正。ZA模型考虑了应变率和温度的相互作用，是基于材料晶体结构和位错机制的一种模型^[15]。但对于不同晶体结构，其流变应力关系式不同，而且该模型中有大量的参量需要通过实验来确定。

本文中采用的材料本构方程模型为EPP模型，它假设材料为各向同性理想弹塑性，材料遵循1维应变条件下的本构关系，当材料中的冲击波峰值压力大于材料的Hugoniot弹性极限(Hugoniot elastic limit, HEL)σ_HEL时，材料产生屈服，发生永久变形。材料的动态屈服强度为^[3]：

式中, ν为泊松比。该模型比较简单，参量易于获取，但是没有考虑材料的应变硬化等因素。根据参考文献[11]，热处理后硬度为30HRC的35CD4钢的机械性能参量为：密度ρ＝7800kg·m^-3，弹性模量E＝210GPa，泊松比ν＝0.29，Hugonoit弹性极限σ_HEL＝1.47GPa。

对于瞬态应力分析，网格密度对模拟结果影响很大。一般网格单元尺寸越细，分析结果越准确，但计算时间也会大大增长。为了精确捕捉到应力波而又尽可能减少计算时间，需要找到一个合适的网格密度。为了评估网格密度对结果的影响，采用最小有限元单元尺寸为0.10mm, 0.05mm, 0.03mm, 0.02mm分别进行模拟计算，如表 1所示。

此外，显式分析增量时间步长Δt对模拟结果的收敛性和准确性有很大的影响。如果增量时间步长大于稳定极限Δt_stable，模拟过程的不稳定性或许会导致无界解。一般而言，稳定极限难以精确测定，可以利用如下所示的计算公式来进行估算:

式中, L_e为最小网格单元尺寸；C_d为弹性波在材料中传播的波速，可以通过公式${C_{\rm{d}}} = \sqrt {\frac{{\left( {1 - \nu } \right)E}}{{\mathit{\rho }\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}} $得到，计算得到C_d=5.94×10³m/s。如果最小网格单元尺寸为L_e=0.03mm，则计算得到的增量时间步长Δt=Δt_stable=5.0ns。

在显式分析计算4500ns后，得到了不同网格单元尺寸下沿着径向的表面动态应力σ_d, 如图 4所示。可以看出, 网格单元尺寸为0.02mm和0.03mm的表面动态应力分布的计算结果基本接近，而与网格单元尺寸为0.10mm的结果相差较大。综合考虑收敛性和计算效率，后续模拟计算均选用有限元单元尺寸为0.03mm。

Figure 4.  Surface dynamic stresses with different FEM element length

采用均匀压力空间分布模型来进行模拟计算，参照参考文献[5]，初步选择的显式求解时间为4500ns，从能量和表面动态应力变化两个方面对激光冲击瞬间靶材的动态响应过程进行分析。

在单次激光冲击材料表面时，靶材能量历史变化曲线如图 5所示。激光冲击波对整个靶材表面产生的外力功W_t在100ns内迅速由0mJ增加到330mJ，并在随后转变为材料的内能W_i、动能W_k和黏性耗散能W_v。1000ns以后，动能和内能急剧减少，分别逐渐趋于0mJ和146mJ，而黏性耗散能则迅速增加到150mJ，最后稳定在185mJ左右。由前面计算可知，弹性波在该材料中传播的波速为5.94×10³m/s，由此可推算出弹性波在1000ns刚好传播至无反射边界处，导致各能量的急剧变化。

Figure 5.  Relationship between energy and time

弹性存储能W_e在90ns内急剧上升到127mJ，然后在1000ns以后急剧减少，最后稳定在17mJ左右。塑性耗散能W_p在300ns内急剧上升到125mJ，在2000ns以后稳定在128.8mJ左右，说明后续没有塑性变形发生，所以选取显式求解时间为4500ns是合理的。此外，伪应变能W_a在整个求解时间范围内不超过0.06mJ，说明采用的网格模型是合理的，在计算过程中不会因为采用了减积分单元CAX4R而引起“沙漏”问题。

由于存在弹性变形和塑性变形，靶材的应力状态只有在某个时间节点以后才能达到稳定状态。图 6所示是在不同时刻沿着轴向的表面动态应力σ_d的变化情况。在500ns时，应力波动幅度较大; 1500ns时, 在光斑中心处附近出现了超过800MPa的拉应力，这是由于光斑边缘产生的稀疏波传至光斑中心所致；而在4000ns以后，应力分布趋于稳定，再次验证了选取求解时间为4500ns的合理性。

Figure 6.  Relationship between surface dynamic stresses and distance from centerline

采用3种不同空间分布冲击载荷模型得到的残余应力模拟结果如图 7所示。其中均匀分布、高斯分布两种模型在距光斑中心0.27mm和4mm处均出现了由于很大应力变化而引起的应力尖峰，且前者的变化程度更剧烈，而NAM的模型分布没有这两个尖峰。相比其它两种分布，均匀分布在距中心1.5mm~4mm范围内残余压应力更趋于稳定。而在光斑中心沿不同深度处的残余应力σ_r分布上，三者的残余压应力影响层深度分别为0.86mm, 0.81mm和0.76mm，并分别在0.33mm, 0.39mm和0.13mm深度处达到最大残余应力-396MPa，-352MPa和-416MPa。这是由于3种空间分布冲击载荷模型在光斑边缘附近的压力大小及分布存在较大差别，这样使得在光斑边缘模拟产生的稀疏波和剪切波大小及特性不一样，最后导致靶材中形成的残余应力场存在一定差距。

Figure 7.  Distribution of residual stresses of different spatial distribution models of the loading

一般来说，对于直径较大的光斑(大于5mm)可以忽略激光脉冲能量空间分布带来的影响^[3]，故本文中采用均匀空间冲击压力模型。模拟得到的沿着表面和不同深度(距光斑中心距离r=3.5mm处)的残余应力σ_r分布如图 8a和图 8b所示，并与BALLARD^[12]的实验测量结果进行对比。BALLARD的实验条件为：直径30mm、厚15mm的35CD4合金钢圆棒试样，激光脉冲能量为150J左右，功率密度为10GW/cm²，激光半峰全宽约为25ns~30ns，光斑直径为8mm，采用黑漆涂层，水作为约束层。

Figure 8.  Residual stresses comparison between simulation and experiment

由图 8中可以看出，模拟结果和实验结果有着较好的一致性。从表面残余应力σ_r的分布来看，模拟得到的最大表面残余压应力为-402MPa，与实验值-360MPa相差11.7%。在距光斑中心4mm左右出现了一个残余应力剧烈变化区域，这是由于光斑边缘区域附近发生了剧烈的剪切变形而产生一定的拉应力，使该区域附近的残余压应力值减少。此外，在距离光斑中心半径0.5mm范围内残余应力剧烈变化，并于0.27mm处达到峰值残余应力179MPa，而在光斑中心处的残余应力值为-132MPa，表现出明显的残余应力缺失。这主要是由于材料表面受到冲击与材料弹性力作用产生振荡过程，冲击光斑边缘产生稀疏波的反向加载^[16]，引起反向塑性变形，从而形成上述的“残余应力洞”现象。

采用ABAQUS软件建立了单次激光冲击35CD4合金钢厚板轴对称有限元模型，验证了显式分析求解时间4500ns选取的合理性。对于2维轴对称模型，选用网格单元尺寸应为0.03mm左右，才能得到收敛的解。不同冲击压力空间分布模型对模拟结果影响较大，主要是由于不同冲击波空间分布在光斑边缘处产生的稀疏波和剪切波强度存在差异。模拟得到的残余应力场计算结果与参考文献中的实验测量结果吻合得较好，表明该模型能够有效预测单次激光冲击诱导的残余应力的大小和分布。