喷丸残余应力场SPH/FEM 耦合建模与仿真

王利平，王建明，裴信超，赵莉莉

(山东大学 机械工程学院，山东 济南250061)

0 引言

喷丸强化是一种能提高工件表面疲劳寿命的工艺．它通过大量弹丸介质以20 ～100 m/s 的速度反复撞击工件表面，在工件表面形成均匀残余压应力层．该残余压应力层可抵消工件表面的拉应力，从而有效阻止裂纹的萌生和扩展，显著提高工件疲劳寿命［1-2］，被广泛应用于航空、军工、汽车等领域．因实验方法耗时耗力，且所用的阿尔门试片不能精确反应残余压应力的大小和分布． 随着计算机技术和有限元技术的发展，已出现许多基于有限元法研究喷丸残余应力的数值模型．Mori K 等［3］(1992)通过建立2D 单丸粒对称模型，对塑性变形工件进行了仿真． S． A． Meguid等［4-5］(1999)通过建立3D 单丸粒动力学模型，研究了弹丸和靶材尺寸对塑性应变及残余应力的影响．S． A． Meguid 等［5］和G． R． Majzoobi 等［6］利用3D 多丸粒对称模型对喷丸强化过程进行了系统的分析．凌祥等［7］(2006)建立了模拟喷丸残余应力场的三维有限元模型．路会龙等［8］(2009)分别建立单、多丸撞击模型，指出弹丸撞击速度存在最优值． 张洪伟等［9］(2010)利用LS-DYNA 有限元软件分别建立单、多丸对称胞元3D 模型，利用单丸模型研究弹丸速度、直径、入射角、摩擦系数和重复撞击次数对残余应力及塑性应变的影响;通过分别取4，6，8，9，13 个弹丸建立对称胞元模型，分析不同覆盖率对残余应力场的影响．上述有限元模型仅涉及单丸或数量有限的多丸，属微观模型，与采用大量弹丸反复撞击工件表面的实际喷丸过程存在明显差异．

笔者采用SPH 耦合FEM 方法针对高强度铝合金A12024-T3 通过显式动力分析软件LS-DYNA，进行了大量弹丸反复撞击工件表面的喷丸强化过程数值模拟;基于两类不同材料属性的SPH粒子模拟弹丸流，研究了弹丸速度、喷丸时间等对喷丸残余应力的影响以及工件表面各典型位置处残余应力的分布规律．

1 数值模型的建立

1．1 两类SPH 粒子均匀混合模型

笔者采用两种SPH 粒子分别表示空气和离散弹丸．其中，空气对应的SPH 粒子采用常规气体状态方程，弹丸对应的SPH 粒子采用离散固体的状态方程．两类粒子间的相互作用通过形函数的特殊定义实现，即在表示粒子形函数的相邻粒子支持域中可包含两类不同性质的SPH 粒子．

在SPH 法中，问题域的离散是由有限个具有质量并且占有一定空间的粒子组成． 粒子的质量可表示为

mi=ΔViρi． (1)

式中:ΔVi为粒子i 占有的体积;ρi为粒子i (i =1，2，…，N)的密度，N 为粒子支持域内的粒子数量．

根据弹丸尺寸及弹丸流中弹丸与空气的体积比，分别确定SPH 粒子大小和空气粒子与弹丸粒子的数量．为使两类SPH 粒子均匀混合，使用随机函数产生随机数，通过MATLAB 编程实现不同属性的两类SPH 粒子在弹丸流空间范围内随机均匀分布．将得到的两类SPH 粒子随机分布单元信息写入LS-DYNA 支持的K 文件中． 图1 所示为两类SPH 粒子均匀分布模型，其中左侧粒子代表弹丸，右侧粒子代表空气．

图1 两类SPH 粒子混合模型Fig．1 Hybrid model for two types of SPH particles

1．2 气体状态方程

在LS-DYNA 中，采用MAT_Null 材料建立空气SPH 粒子，使用状态方程EOS_LINEAR_ POLYNOMIAL 表示压力与密度变化间的关系． EOS_LINEAR_POLYNOMIAL 为线性多项式形式的状态方程，可表示为

Ρ=C0+C1μ +C2μ2+C3μ3+(C4μ4+C5μ5+C6μ6)Ε． (2)

式中:μ =ρ/ρ0-1，表示当前密度与参考密度间的比值，ρ0为初始参考密度，在材料模式MAT_Null中定义;C0～C6为各项系数，对于气体状态方程取

式中:Cp为定压比热;Cv为定容比热．

上述状态方程相关参数使用AUTODYNAVersion6．1 材料库中的参数，如表1 所示．

表1 气体材料和状态方程相关参数Tab．1 Relative parameters of the air material and state equation

1．3 丸粒状态方程

丸粒SPH 粒子使用MAT_Null 模拟，考虑弹丸流具有离散介质的特性，根据文献［10］推导弹丸SPH 粒子的状态方程，其为特殊形式的Gruneisen 方程:

式中:Cs为声音在弹丸材料中的传播速度，对于钢丸其值为5．2 km/s;ρ0为弹丸材料初始密度．

故弹丸SPH 粒子也可以使用线性多项式形式的状态方程EOS_LINEAR_POLYNOMIAL，此时有

弹丸材料和状态方程参数如下:密度ρ0=7．8 g·cm-3;多项式系数C1=2．096．

1．4 工件材料本构方程

工件材料为A12024-T3，力学性能参数如下:弹性模量E=71．3 GPa;密度ρ0=2．78 g·cm-3;泊松比γ=0．3．

由于撞击过程速度较高且高应变率，需要考虑温度对应力流动的影响． 工件材料采用Johnson-Cook 黏塑性本构方程模拟，屈服应力σy与等效塑性应变的关系式为

式中:A、B、C、n、m 是通过实验获得与应变、应变率、温度相关的常数． Dabboussi 等人提出的相关参数使得仿真结果与实验结果具有较高的一致性［11］，见表为等效塑性应变=等于1 s-1．

式中:T 为任一时刻的瞬时温度;Tm为熔解温度，Tr为参考温度取292K．

表2 Johnson-Cook 相关参数Tab.2 Relative parameters of Johnson-Cook

1．5 SPH 耦合FEM 的喷丸模型

采用LS-DYNA 软件建立喷丸模型如图1 所示．其中弹丸流使用两类SPH 粒子建模，共包含320 个SPH 粒子;工件使用有限元建模，为减少计算时间，工件尺寸取6 mm×6 mm×4 mm，由于工件中部区域与弹丸发生接触，故对其网格进行细化，由53 760 个solid 单元组成．工件底面采用固定位移边界，4 个侧面采用非反射边界以避免应力波在其边界发生反射．

通过定义接触算法耦合FEM 与SPH 模型．在LS-DYNA 中通过关键字Contact_Automatic_Nodes_to_Surface 定义有限单元模型上表面与SPH 粒子间的接触关系，其中将FEM 单元节点定义为主节点，SPH 粒子定义为从节点． SPH 耦合FEM 的过程如图2 所示．

图2 SPH 耦合FEM 的过程Fig．2 The process of SPH coupled FEM

2 仿真结果分析值模型

图3 为大量弹丸以速度50 m/s 撞击工件时不同时刻的撞击效果图． 基于此模型可模拟大量弹丸撞击工件的喷丸过程．

图3 大量弹丸撞击工件效果图Fig．3 The process of shots impacting the target

2．1 弹丸速度对残余应力的影响

由于本文仿真模型中，弹丸流由大量随机分布的丸粒组成，可通过弹坑附近的残余压应力研究大量丸粒撞击下的喷丸结果． 图4 给出弹丸速度分别为25，50，75 和100 m/s 的大量弹丸经过120 μs 撞击后，弹坑底部残余应力分布图． 对应上述4 种弹丸速度，其最大残余应力分别为-632．2，-692．71，-729．0，-738．9 MPa，发生在距离弹坑表面0．18 ～0．24 mm 的次表层．最大残余压应力值和其深度均随速度的增加而增加．因此在实际生产中，应根据最大残余压应力及应力层深度要求合理选择喷丸速度．

2．2 喷丸时间对残余应力的影响

基于本文含大量丸粒的宏观喷丸模型，可研究残余应力场随喷丸时间的变化． 随着喷丸时间的增加，撞击工件表面的弹丸数量、工件表面的有效受喷面积及覆盖率均随之增加． 图5 所示为速度为50 m/s 的弹丸分别经过26．6，79．8，102．6，120 μs 时间后，弹坑底部残余应力分布图． 图示弹坑处受多个不同粒子撞击，其中26．6 μs 对应弹丸首次撞击后的残余应力结果，而79．8 μs 则对应弹丸多次撞击同一位置后的残余应力结果．两者对应的最大残余压应力值分别为-658．86 MPa 和-692．17 MPa，其最大残余压应力值有所提高．撞击时间继续增加至102．6 μs 和120 μs时，其对应的残余应力变化甚微，最大变化幅度小于10%，即喷丸残余应力存在饱和状态．

2．3 不同位置点的残余应力分布

图4 喷丸速度对残余应力的影响Fig．4 Effect of impact velocity on residual stress

图5 喷丸时间对残余应力的影响Fig．5 Effect of impact time on residual stress

图6 典型弹坑位置点Fig．6 Typical positions of craters

图6 所示为一典型四弹坑分布图． C、E、G、I分别为4 个弹坑的中心，B、D、H、F 分别为4 个弹坑两两重叠区域，A 点未受到撞击． ABCD 区域可视为一典型区域．图7 分别表示A、B、C、D 4 个位置点的残余应力分布情况． C 点为粒子撞击中心点，其最大残余压应力为-653．45 MPa，发生在深度为0．20 ～0．22 mm 的次表层;B 点和D 点位于两弹坑重叠处，该处的最大残余压应力值为-343．56 MPa，位于深度为0．12 ～0．18 mm 的次表层;A 点的最大残余压应力值为-450 MPa，发生在表面处．

图7 不同位置残余应力分布Fig．7 Distribution of the residual stress vs． depth at different positions

3 结论

针对有限元法难以模拟大量弹丸群集撞击的喷丸过程，采用SPH 耦合FEM 法进行了喷丸过程数值建模与分析．仿真结果如下:

(1)最大残余压应力值及所在次表层深度均随弹丸速度的增加而增加． 故应根据最大残余压应力及应力层的深度要求合理选择喷丸速度．

(2)开始时，随着喷丸时间的增加，最大残余应力值增大;达到一定时间后，继续增加喷丸时间，残余应力变化甚微，最大残余压应力的变化幅度小于10%，达到喷丸饱和状态．

(3)弹坑中心点处最大残余压应力值较高且位于次表层;两坑交汇处的最大残余压应力值较中心点低且其次表层深度较浅;未受撞击点处的最大残余压应力发生在表面上．

［1］ KAZUYUKI O． Fatigue life enhancement of aluminum alloy for aircraft by Fine Particle Shot Peening(FPSP)［J］． Journal of Materials Processing Technology，2011，211(8):1395 -1399．

［2］ 高玉魁．喷丸强化对23CO14Ni12Cr3MoE 钢疲劳性能的影响［J］． 材料热处理学报，2007，28(6):75 -77．

［3］ MORI K，OSAKADA K，MATSUOKA N． Finite element analysis of peening process with plastically deforming shot［J］． Mater Process Technol，1994，45(1 -4):607 -612．

［4］ MEGUID S A，SHAGAL G，STRANART J C． Finite element modelling of shot-peening residual stresses［J］． Journal of Materials Processing Technology，1999，92/93:401 -404．

［5］ MEGUID S A，SHAGAL G，STRANART J C，et al．Three-dimensional dynamic finite element analysis of shot-peening induced residual stresses［J］． Finite Elements in Analysis and Design，1999，31(3):179-191．

［6］ MAJZOOBI G H，AZIZI R，NIA A A． A three-dimensional simulation of shot peening process using multiple shot impacts ［J］． Journal of Materials Processing Technology，2005:164/165:1226 -1234．

［7］ 凌翔，彭薇薇，倪红芳． 喷丸三维残余应力场的有限元模拟［J］． 机械工程学报，2006，42 (8):182 -189．

［8］ 路会龙，姚平喜，刘海英．基于ANSYS/LS-DYNA 的受控喷丸工艺过程仿真［J］．机械设计与制造，2009(2):214 -216．

［9］ 张洪伟，张以都，吴琼．喷丸强化残余应力场三维数值分析［J］． 航空动力学报，2010，25(3):603-609．

［10］ GRUJICIC M，PANDURANGAN B，QIAO R，et al．Parameterization of the porous-material model for sand with different levels of water saturation［J］． Soil Dynamics and Earthquake Engineering，2008，28(1):20 -35．

［11］ DABBOUSSI W，NEMES J A． Modeling of ductile fracture using the dynamic punch test ［J］． International Journal of Mechanical Sciences，2005，47(8):1282 -1299．