CN109725123B - 一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法 - Google Patents

Abstract

本发明涉及一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法，步骤为：(1)基于J‑C模型应力应变数据进行拟合，得到材料的位错演化模型参数；(2)利用ABAQUS有限元软件的VUSDFLD子程序进行编程，基于位错演化模型建立微细观参量与宏观参量的联系，完成程序编写；(3)基于ABAQUS软件进行喷丸强化数值仿真，获取强化之后表层晶粒尺寸与残余应力场；(4)根据相同材料不同应力下未喷丸的裂纹扩展试验数据，拟合得到N‑R模型参数，建立裂纹扩展模型；(5)将残余应力与晶粒细化考虑到N‑R裂纹扩展模型中，建立考虑喷丸诱导晶粒细化的裂纹扩展模型。

Description

一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法

技术领域

本发明是一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法，它是一种能够考虑喷丸强化之后微观结构变化对裂纹扩展影响的喷丸构件裂纹扩展寿命确定方法，属于航空航天发动机技术领域。

背景技术

先进航空发动机零部件结构复杂，工作环境苛刻，特别是涡轮部件长时间工作在高温、高压、高转速的极端服役环境中，在加工制造以及使用的过程中不可避免引入萌生裂纹。对于已经存在实际裂纹的结构，抑制裂纹扩展过程从而增加疲劳扩展寿命，是降低航空发动机在服役条件下的疲劳失效风险，保障航空发动机零部件长期安全使用的重要研究方向。喷丸强化技术是工业上常用的表面强化技术，喷丸之后的表面不仅仅会引入残余应力场，而且会使表层晶粒尺寸发生变化，最典型的就是晶粒会产生细化。晶粒细化会提高材料晶界的含量，晶界对裂纹扩展具有阻碍作用，降低裂纹扩展速率。现有的喷丸强化裂纹扩展评估方法往往仅仅考虑残余应力的影响，没有考虑晶粒细化带来的组织强化作用，预测结果相对于实际存在偏差。因此需要建立有效的喷丸强化晶粒细化模拟方法，综合考虑强化效果对裂纹扩展的影响。

现有文献“王建明,赵莉莉,吕鹤婷.喷丸残余应力对裂纹扩展疲劳寿命的影响[J].哈尔滨工程大学学报,2016,37(04):608-613.”基于ABAQUS软件建立含喷丸残余应力场的四点弯曲试样有限元模型和考虑裂纹闭合效应的裂纹扩展疲劳寿命预测模型，实现裂纹扩展寿命预测，但是没有建立喷丸强化过程晶粒细化的模拟方法，因此没有考虑其对裂纹扩展的影响。

发明内容

本发明技术解决方案：克服现有技术的不足，提供一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法，综合反映喷丸强化引入的残余应力与表层晶粒尺寸变化对裂纹扩展的影响，实现裂纹扩展寿命的准确预测，能够服务与支撑航空发动机构件喷丸之后裂纹扩展寿命分析。

本发明技术解决方案：一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法，基于位错密度演化模型实现喷丸强化过程晶粒尺寸变化的模拟，基于N-R小裂纹扩展模型将晶粒尺寸的影响考虑进来，实现考虑残余应力晶粒细化等喷丸强化综合效果的裂纹扩展寿命预测。实现步骤如下：

第一步，基于J-C模型应力应变数据进行拟合，得到材料的位错演化模型参数，用以计算材料表层的晶粒细化；设置位错演化模型需要拟合的参数的初始值，代入位错演化模型中计算得到应力。然后将计算得到的高应变率下的应力结果与J-C模型计算得到的结果进行对比，利用遗传算法将两种模型得到的结果误差减小到最小值，从而得到材料的位错演化模型参数。位错密度演化模型如下式：

式中，ρ_c、ρ_w分别为位错胞内、胞壁中的位错密度；α^*、β^*、k₀和n^*分别为控制位错密度演化率的参数；b为伯氏矢量；

分别为位错胞内、胞壁中的分切剪应变率，由于胞壁、胞内的边界需保证应变协调，通常认为这两个参量相同，即

为工艺的参考分切剪应变率；d为晶粒尺寸；f为位错胞壁所占的体积分数，通常可通过TEM试验检测获得。

第二步，利用ABAQUS有限元软件，将第一步中材料的位错演化模型参数带入位错演化模型，在有限元软件中建立细观参量与宏观参量的联系。利用SEM、TEM等技术对材料的晶粒尺寸与位错密度进行测量，将测量结果作为数值模拟的输入因素，根据位错胞内部位错密度的演化规律，如胞内、胞壁边界Frank-Read源产生的位错、胞内向胞壁转移的位错和由横向滑移导致的位错湮灭等，计算某时刻位错胞内、胞内位错密度，并按体积分数叠加计算出整体分解剪应力。利用可被ABAQUS有限元软件兼容的Fortran语言编程，实现有限元模拟计算。

第三步，基于ABAQUS软件建立有限元模型进行喷丸强化数值仿真，利用步骤二建立的位错演化模型，获取强化之后晶粒尺寸与残余应力场，作为后续裂纹扩展计算的输入参数。根据实际构件几何形状完成喷丸强化靶材有限元模型构建，根据实际情况设置边界条件，采用适于高应变率下的本构模型J-C模型。根据所需要模拟的喷丸覆盖率，利用Avrami公式预估所需要弹丸的个数，利用ABAQUS的Python脚本语言进行编程，给定弹丸生成位置并给定弹丸速度大小与方向，定义弹丸与靶材的接触。对靶材和弹丸进行网格划分，靶材边界设置成无限单元以控制应力场的稳定。计算时同时提交Fortran程序，完成常规计算的同时计算出晶粒尺寸的变化。提取喷丸强化之后的残余应力场与晶粒尺寸。根据有限元模拟得到的裂纹扩展路径上的残余应力与晶粒尺寸，拟合得到裂纹长度与残余应力、晶粒尺寸的关系曲线。

第四步，根据相同材料不同应力下未喷丸的裂纹扩展试验数据拟合得到N-R模型参数，建立裂纹扩展模型，用于后续裂纹扩展模拟计算。通过开展试验或者阅读文献的方式，获取未强化试验件裂纹长度与循环数的数据。根据模型分析结果计算出每个试验点对应的裂纹尖端塑性区位移与裂纹扩展速率，拟合得到模型参数，完成裂纹扩展模型构建。

第五步，根据步骤(3)得到的晶粒尺寸与残余应力分布结果，将残余应力与晶粒尺寸变化考虑到N-R模型中，根据步骤(4)得到的N-R模型结果建立考虑喷丸诱导的晶粒细化的裂纹扩展模型，以晶粒为单位，不断进行裂纹扩展计算，更新晶粒尺寸与应力场并计算循环数，得到相应的裂纹扩展寿命。重复上述过程直到达到给定裂纹长度为止。

第五步中，将残余应力与晶粒尺寸变化考虑到N-R模型中的方法是通过拟合残余应力、晶粒尺寸随裂纹长度的变化关系式，如下式所示，将关系式引入修正的N-R模型中，在完成上一个晶粒计算的时按照拟合的关系式不断更新下一个晶粒尺寸与应力值，

d＝f(a)

σ＝g(a)

其中d为晶粒尺寸，a为裂纹长度，σ为残余应力、g(a)为多项式的形式，取三次多项式，即：

f(a)/g(a)＝C₁a³+C₂a²+C₃a+C₄

其中C₁、C₂、C₃、C₄为多项式系数。

本发明与现有技术相比的优点在于：

(1)本发明基于位错演化模型构建喷丸强化多尺度模拟方法，可有效模拟弹丸冲击过程表层晶粒尺寸变化情况，现有喷丸强化数值模拟方法往往难以模拟晶粒尺寸变化的过程，因此无法将晶粒尺寸变化考虑进来。

(2)本发明利用N-R模型将喷丸强化的效果包括残余应力场与晶粒尺寸变化考虑到裂纹扩展模拟中，添加相应的修正项，实现裂纹扩展的模拟，更加符合喷丸强化实际情况，裂纹扩展寿命计算更加准确。现有裂纹扩展寿命分析方法往往仅仅考虑残余应力的影响，没有考虑晶粒尺寸的变化对裂纹扩展过程的影响，裂纹扩展寿命准确性不能满足需求。

附图说明

图1是本发明的一种考虑喷丸强化晶粒尺寸变化的裂纹扩展寿命分析方法流程图；

图2是进行喷丸强化模拟构建的两弹丸冲击有限元模型；

图3是进行喷丸强化多尺度模拟得到的晶粒尺寸云图，其中(a)是一次冲击结果，(b)是两次冲击结果；

图4是进行喷丸强化多尺度模拟结果，(a)是晶粒尺寸曲线(b)是残余应力曲线；

图5是开展裂纹扩展试验所用CT试验件图纸；

图6是计算得到的裂纹扩展速率曲线；

图7是计算得到的裂纹长度与寿命关系曲线。

具体实施方式

下面结合附图，通过举例的方式，对本发明一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命分析方法做进一步说明。

如图1所示，本发明实现过程如下：

第一步，基于J-C模型应力应变数据进行拟合，得到材料的位错演化模型参数，用以计算材料表层的晶粒细化；J-C模型是常用的描述高应变率下材料应力应变关系的本构模型，J-C模型的表达式为：

式中，σ_eq为流动应力，A、B、C和m为材料常数，n是加工硬化系数，

为标准化的有效塑性应变率，T_melt为熔点，T_r为室温，T_r＝298K，ε为等效塑性应变，

是参考应变率，取值为1。

对于预研究材料GH4169，通过查阅文献的方式，确定J-C模型的参数A＝1241MPa，B＝622，C＝0.0134，m＝1.05，n＝0.6522。给定应变、应变率与温度项，可以得到材料的流动应力。

设置位错演化模型需要拟合的参数的初始值，代入模型中计算得到应力。然后将计算得到的高应变率下的应力结果与J-C模型计算得到的结果进行对比，设置最大允许的拟合误差，利用遗传算法不断更新位错演化模型的参数，将两种模型得到的结果误差减小到允许值，从而得到材料的位错演化模型参数。

第二步，利用ABAQUS有限元软件，将第一步中材料的位错演化模型参数带入位错演化模型，在有限元软件中建立细观参量与宏观参量的联系。利用SEM、TEM等技术对材料的晶粒尺寸与位错密度进行测量，将测量结果作为数值模拟的输入因素，根据位错胞内部位错密度的演化规律，如胞内、胞壁边界Frank-Read源产生的位错、胞内向胞壁转移的位错和由横向滑移导致的位错湮灭等，计算某时刻位错胞内、胞内位错密度，如下式所示：

式中，ρ_c、ρ_w分别为位错胞内、胞壁中的位错密度；α^*、β^*、k₀和n^*分别为控制位错密度演化率的参数；b为伯氏矢量；

分别为位错胞内、胞壁中的分切剪应变率，由于胞壁、胞内的边界需保证应变协调，通常认为这两个参量相同，即

为工艺的参考分切剪应变率；d为晶粒尺寸；f为位错胞壁所占的体积分数，通常可通过TEM试验检测获得。

按体积分数叠加计算出整体分解剪应力。利用ABAQUS有限元软件中适合动态分析的VUSDFLD子程序将上述公式编程，实现有限元模拟计算。

第三步，基于ABAQUS软件建立有限元模型进行喷丸强化数值仿真，利用第二步建立的位错演化模型，获取强化之后晶粒尺寸与残余应力场，作为后续裂纹扩展计算的输入参数。以最简单的单丸粒冲击与双丸粒冲击模型为例说明上述过程。模型如图2所示，靶材表面为正方形。弹丸的尺寸均为1.0mm，入射速度v＝70m/s，在靶材上的投影位于靶材正中央。其中最外围区域的网格为无限单元CIN3D8，以消除冲击造成的应力波激荡。靶材其余区域与弹丸为六面体减缩积分单元C3D8R，靶材底面进行位移全约束。冲击区域设置沙漏(hourglass)控制，对弹丸冲击接触区域进行网格细分，网格尺寸为0.05mm×0.05mm×0.05mm。弹丸与靶材之间设置为面面接触，接触面之间的摩擦系数设定为0.3。

采用Johnson-Cook(J-C)模型，可有效描述材料在高应变率状态下的应力应变关系。利用ABAQUS的Python脚本语言对上述过程进行编程。计算时同时提交Fortran程序，完成常规计算的同时计算出晶粒尺寸的变化。提取喷丸强化之后的残余应力场与晶粒尺寸。一次冲击与二次冲击之后的弹丸的晶粒尺寸分布情况如图3所示，可以看出喷丸冲击之后表层的晶粒发生明显细化，两次冲击的晶粒尺寸细化程度比一次冲击的晶粒尺寸细化程度更加明显。

对获得的晶粒尺寸分布与残余应力分布进行后处理，提取沿靶材厚度方向的分布情况。由于一次冲击与二次冲击计算过程相近，因此以一次冲击为例说明计算过程。

冲击之后靶材中心沿深度方向的晶粒尺寸分成3部分，分别进行幂次拟合，拟合结果如下式，拟合曲线如图4(a)，晶粒尺寸随深度增大逐渐增大直到不发生变化：

其中d是晶粒尺寸，单位是um，x是裂纹长度，单位是mm。

对深度方向的残余应力分成3部分，分别进行拟合，结果如下式，拟合曲线如图4(b)，残余应力随深度增大先增大后减小：

其中σ_res是残余应力，单位是MPa，x是裂纹长度，单位是mm。

由此建立了裂纹扩展路径上晶粒尺寸与残余应力随裂纹长度变化的关系式。

第四步，根据相同材料不同应力下未喷丸的裂纹扩展试验数据拟合得到N-R模型参数，建立裂纹扩展模型，用于后续裂纹扩展模拟计算。基于国家标准GB/T6398-2000《金属材料疲劳裂纹扩展速率试验方法》的要求，开展裂纹扩展标准CT件试验，试验件材料为GH4169，试验件图纸如图5所示。试验在两种不同的载荷下进行，最大载荷分别是4kN与6kN，应力比均为0.1，获取裂纹长度与循环数数据点，利用割线法分别计算得到每个载荷下每个点对应的裂纹扩展速率。

根据N-R模型，按下式迭代计算每个点对应的n_c，在这里取迭代差值为0.0001。由于n_c的值未知，因此首先假设初始的n_c＝1。

其中n是裂纹尖端无量纲位置，数学表达式为n＝a/c，其中a是裂纹长度，c是裂纹长度加上裂纹尖端塑性区的长度。n_c是n在一个晶粒的最大值，σ为外加应力，σ_comp是大于外加应力σ的比较应力，K为应力强度因子，K_th是应力强度因子门槛值

计算出n_c之后，按下式计算n_s。

其中n是裂纹尖端无量纲位置，n_c是n在一个晶粒的最大值，n_s是n在一个晶粒的最小值；σ为外加应力，σ_FL是疲劳门槛应力，K为应力强度因子，K_th是应力强度因子门槛值。

然后将n_s作为试验点对应的n，代入下式计算出裂尖塑性区位移

其中n是裂纹尖端无量纲位置，σ为外加应力，K为应力强度因子，b是伯氏矢量，A₁的计算方法如下式。

其中G是剪切模量，v是泊松比。

计算出每个试验点对应的裂纹尖端塑性区位移

后，假定裂纹扩展速率da/dN与

满足下式。

da/dN＝f_aφ

利用最小二乘法进行拟合直线，分别对两个载荷下的数据进行上述操作，获得da/dN与

的两个系数f_a1、f_a2，按照下式分别计算p₁、p₂。

其中Δσ为外加应力，p₁、p₂为两个常数。

由此根据裂纹扩展试验数据，拟合得到N-R模型参数，建立了裂纹扩展模型。

第五步，根据第三步得到的晶粒尺寸与残余应力分布结果，将残余应力与晶粒尺寸变化考虑到N-R模型中，根据第四步得到的N-R模型结果建立考虑喷丸诱导的晶粒细化的裂纹扩展模型，得到相应的裂纹扩展寿命。N-R模型的基本思路是以晶粒为单位，从第一个晶粒开始，计算裂纹扩展到每个晶粒时所对应的n_c与n_s。在n_s与n_c之间均分出8个点，加上n_s与n_c，共有10个点。对于给定的晶粒，c已知，由可以计算出每个点对应的裂纹长度a，进而得到每个点对应的裂纹尖端塑性区位移

再乘以比例因子f，得到每个点对应的裂纹扩展速率。将各个裂纹扩展速率点连接起来，即可做出在一个晶粒内的裂纹扩展速率曲线。按照这个步骤，继续计算裂纹在下一个晶粒扩展时的c、n_s、n_c并对n进行均分。

喷丸强化之后引入残余应力，使得裂纹扩展过程的实际应力变化。假设实际应力满足叠加原理，即：

σ_t＝σ+σ_res

其中σ_t是实际应力，σ是外加应力，σ_res为残余应力。

根据第三步建立的裂纹扩展路径上晶粒尺寸与残余应力随裂纹长度变化的关系式，利用MATLAB编程计算，完成一个晶粒的计算之后，计算已经计算的晶粒的总长度，代入拟合公式求解下一个晶粒位置的残余应力与晶粒尺寸。将残余应力项代入修正此时的实际应力项与求解得到的晶粒尺寸一起代入，求解裂纹扩展速率，求解得到的一次冲击、二次冲击以及未强化的裂纹扩展速率随裂纹长度的变化如图6所示，其中横坐标为裂纹长度，纵坐标为裂纹扩展速率。可以看出未冲击的裂纹扩展速率最快，二次冲击相对一次冲击裂纹扩展速率更小。

完成裂纹扩展速率分析之后即可对裂纹扩展寿命进行分析，根据模型假设有：

其中da/dN是裂纹扩展速率，f是比例因子，v是泊松比，G是剪切模量，n是裂纹尖端无量纲位置参数，c是裂纹长度加裂纹尖端塑性区总长度，Δσ是应力幅值。

因为在塑性区被一个晶粒内的晶界阻拦时，c保持不变。因此有：

因此有：

在一个晶粒内扩展时，n值从n_s增长到n_c。将等式左边从n_s到n_c进行积分，就可以得到在裂纹在一个晶粒内扩展所需要的循环数ΔN，如下式所示：

进而计算得到给定裂纹长度下的循环数。计算得到裂纹长度与循环数的曲线如图7所示，其中横坐标为循环数，纵坐标为裂纹长度，可以看出未冲击的裂纹扩展寿命最短，二次冲击相对于一次冲击强化效果更好。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (4)

1.一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法，其特征在于：包括步骤如下：

步骤(1)：基于J-C模型应力应变数据进行拟合，得到材料的位错演化模型参数，用以计算材料表层的晶粒细化；所述应力应变数据是J-C模型刻画的应力应变曲线上取点获得的数据；所述位错演化模型参数是指基于位错理论建立的描述位错演化的模型中未知的参数；

步骤(2)：利用ABAQUS有限元软件，将步骤(1)中材料的位错演化模型参数带入位错演化模型，在有限元软件中建立细观参量与宏观参量的联系；所述细观参量是指晶粒尺寸；所述宏观参量是指应力场应变场参数；

步骤(3)：基于ABAQUS软件建立有限元模型进行喷丸强化数值仿真，利用步骤(1)的位错演化模型，获取强化之后晶粒尺寸与残余应力场，作为后续裂纹扩展计算的输入参数；所述晶粒尺寸与残余应力场是指喷丸强化模拟之后获得的表层晶粒尺寸与残余应力场分布情况；

步骤(4)：根据相同材料不同应力下未喷丸的裂纹扩展试验数据拟合得到N-R模型参数，建立裂纹扩展模型，用于后续裂纹扩展模拟计算；所述相同材料指和欲进行寿命预测的喷丸之后的试验件进行完全相同热处理的没有喷丸的试验件；所述裂纹扩展试验数据是指通过试验获得的裂纹长度与循环数数据点，通过材料手册进行获取；所述N-R模型参数是指N-R模型中表征裂纹扩展速率与裂尖塑性位移的两个与材料相关的常数；

步骤(5)：根据步骤(3)得到的晶粒尺寸与残余应力场分布，将残余应力场分布与晶粒尺寸变化考虑到N-R模型中，根据步骤(4)得到的N-R模型建立考虑喷丸诱导的晶粒细化的裂纹扩展模型，以晶粒为单位，不断进行裂纹扩展计算，得到相应的裂纹扩展寿命。

2.根据权利要求1所述的一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法，其特征在于：所述步骤(1)中，得到材料的位错演化模型参数实现如下：

设置位错演化模型需要拟合的参数的初始值，代入位错演化模型中计算得到应力，然后将计算得到的高应变率下的应力结果与J-C模型计算得到的结果进行对比，利用遗传算法将位错演化模型与J-C模型得到的结果误差减小到最小值，从而得到材料的位错演化模型参数。

3.根据权利要求1所述的一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命确定方法，其特征在于：所述步骤(2)中，建立位错密度演化模型所需要的输入参数方法为利用SEM或TEM技术对材料的晶粒尺寸与位错密度进行测量，将测量结果作为数值模拟的输入参数。

4.根据权利要求1所述的一种考虑喷丸强化表层晶粒细化的裂纹扩展寿命分析方法，其特征在于：步骤(5)中，将残余应力与晶粒尺寸变化考虑到N-R模型中的方法是通过拟合残余应力、晶粒尺寸随裂纹长度的变化关系式，如下式所示，将关系式引入修正的N-R模型中，在完成上一个晶粒计算时按照拟合的关系式不断更新下一个晶粒尺寸与应力值，

d＝f(a)

σ＝g(a)

其中d为晶粒尺寸，a为裂纹长度，σ为残余应力、g(a)为多项式的形式，取三次多项式，即：

f(a)/g(a)＝C₁a³+C₂a²+C₃a+C₄

其中C₁、C₂、C₃、C₄为多项式系数。

Images