CN109558692B - 预测微粒冲击金属件残余应力和马氏体相变的有限元方法 - Google Patents

Abstract

本发明公开了一种预测微粒高速冲击下金属构件残余应力和马氏体相变的有限元方法，包括以下步骤：(1)建立微粒冲击金属构件的有限元模型；(2)建立综合考虑了塑性应变强化和马氏体相变强化的率相关弹塑性本构模型；(3)编写ABAQUS‑VUMAT用户动态材料子程序，运用应力补偿更新算法实现提出的率相关弹塑性本构模型；(4)将子程序嵌入到ABAQUS中，对微粒流冲击金属构件的过程进行计算，进一步获得金属构件内部的残余应力以及马氏体相的百分比含量。本发明可以用于任意数量的微粒对任意结构形式的金属构件进行冲击的残余应力和马氏体相变预测。

Description

预测微粒冲击金属件残余应力和马氏体相变的有限元方法

技术领域

本发明涉及金属构件残余应力和组织相变预测，旨在提供预测微粒高速冲击下金属构件残余应力和马氏体相变的有限元方法。

背景技术

微粒流冲击采用大量微粒反复冲击金属构件，引入有益的残余压应力，残余压应力可促使金属表面裂纹闭合，抑制裂纹扩展，提高疲劳裂纹扩展寿命，改善金属构件的抗疲劳性能。而且表层组织发生塑性变形，发生形变诱发马氏体相变，进而产生马氏体相变强化。因此，研究残余应力产生和马氏体相变对推进微粒冲击技术的应用与发展非常关键。由于实验研究费时费力，而且在实际的工程应用中还会受到一定的制约，因此，有限元法成为主要的研究手段。

构建微粒高速冲击金属构件的数值模型时，准确描述金属材料的本构关系至关重要。近年来研究中材料的本构关系多采用Johnson-Cook模型，该模型考虑了应变硬化和应变率硬化的影响，但是没有考虑到马氏体相变强化的影响，残余应力的预测结果并不准确，更无法对高速冲击下金属材料马氏体相变进行研究。因此，提出一种可以准确预测微粒高速冲击下金属构件内残余应力和马氏体相变的材料率相关弹塑性本构模型尤为必要。

发明内容

本发明的主要目的在于克服现有技术中的不足，提供一种能准确预测微粒高速冲击下金属构件残余应力和马氏体相变的有限元方法。

为解决上述技术问题，本发明的解决方案是：一种预测微粒高速冲击下金属构件残余应力和马氏体相变的有限元方法，包括以下步骤：

(1)建立微粒冲击金属构件的有限元模型；

(2)建立综合考虑了塑性应变强化和马氏体相变强化的率相关弹塑性本构模型；

(3)编写ABAQUS-VUMAT用户动态材料子程序，运用应力补偿更新算法实现提出的率相关弹塑性本构模型；

(4)将子程序嵌入到ABAQUS中，对微粒流冲击金属构件的过程进行计算，进一步获得金属构件内部的残余应力以及马氏体相的百分比含量。

进一步的，所述步骤(1)具体为：

基于ABAQUS建立单个或多个球形微粒、金属构件，分别设置材料属性和划分网格，再进行ASSEMBLY对其组装之后设置分析步和通用接触属性。

进一步的，所述步骤(2)具体包括下述步骤：

步骤(2.1)建立金属材料的本构关系：

Δσ＝D：Δε^e

其中，Δσ为应力增量张量，D为弹性模量张量，Δε^e为弹性应变增量张量；

步骤(2.2)建立塑性模型，具体建立方式为：

考虑塑性应变强化和马氏体相变强化的流动准则为：

其中，f为屈服方程；σ_e为Mises等效应力；所述K和分别为各向同性强化应力和等效塑性应变；所述σ₁₁，σ₂₂，σ₃₃，σ₁₂，σ₂₃和σ₃₁为六个方向的应力分量。

各向同性强化应力：

k_w＝h_wwⁿ

其中，k_ε为塑性应变强化应力，k_w为马氏体相变强化应力，C为材料参数，为等效塑性应变率，为参考应变率，σ₀为初始屈服应力，h_ε为应变强化系数，a为材料参数，h_w为相变强化模量，n为强化指数，w为马氏体的百分比含量，其增量方程由Santacreu模型给出：

其中，w_max为马氏体相的最大含量，D₀，D₁和m为材料参数，为应力三轴度。

进一步的，所述步骤(3)具体包括下述步骤：

步骤(3.1)：通过用户子程序VUMAT的SDV定义第n+1增量步开始时的初始状态变量w_n值，同时也是第n步结束时的状态变量值，在第n+1增量步开始时读入；其中，所述n是指第n增量步，w_n是指第n增量步结束时的金属构件马氏体的百分比含量；

步骤(3.2)：VUMAT由应变增量驱动，首先假设应变增量全为弹性，由步骤(2.1)中公式Δσ＝D：Δε^e计算试探应力：将试探应力代入到屈服准则中，该公式是步骤(2.2)中的公式在第n+1步的特定计算；

如果f_n+1≤0，则材料尚处于弹性变形阶段，将试探应力更新为n+1增量步应力，马氏体相含量更新为n+1增量步状态变量w_n+1；

如果f_n+1＞0，则材料屈服，根据应力补偿显式算法更新增量步结束时的应力，即在试探应力的基础上减去多算的塑性应变增量对应的应力增量部分，对试探应力进行负补偿；

等效塑性应变增量运用牛顿迭代法求解：

其中，所述i表示第i次迭代，表示第i次迭代所得的等效塑性应变增量；所述i-1 表示第i-1次迭代，表示第i-1次迭代所得的等效塑性应变增量；更新等效塑性应变增量，直至满足精度要求，得到马氏体相的百分比含量w_n+1，更新第n+1增量步应力，更新内能，塑性功。

进一步的，所述步骤(4)具体为：

将步骤(1)建立的模型主文件和步骤(3)建立的ABAQUS-VUMAT用户子程序联合，使用ABAQUS/EXPLICT方法对微粒高速冲击金属构件进行计算，获得马氏体相变含量；将ABAQUS/EXPLICIT中得到的动态应力状态导入到ABAQUS/STANDARD模块中进行隐式分析，确定静态平衡下的残余应力场，即完成微粒流高速冲击下金属构件内残余应力和马氏体相变的预测。

与现有技术相比，本发明的有益效果是：

(1)建立的材料率相关弹塑性本构模型同时考虑了塑性应变强化和马氏体相变强化的影响。

(2)利用ABAQUS-VUMAT用户子程序来数值实现所建立的材料本构模型，能准确预测微粒高速冲击下金属构件残余应力和马氏体相的百分比含量。

(3)基于通用有限元技术，本发明可以用于任意数量的微粒对任意结构形式的金属构件进行冲击的残余应力和马氏体相变预测。

附图说明

图1为本发明对某一具体实例所采用的有限元模型；

图2为本发明对所提出的率相关弹塑性本构模型的VUMAT数值实现流程图；

图3为对图2中实例进行预测得到的残余应力分布结果图；

图4为对图2中实例进行预测得到的马氏体相含量结果图。

具体实施方式

首先需要说明的是，本发明是计算机技术在残余应力和组织相变预测领域的一种应用。在本发明的实现过程中，会涉及到多个软件功能模块的应用。申请人认为，如在仔细阅读申请文件、准确理解本发明的实现原理和发明目的以后，在结合现有公知技术的情况下，本领域技术人员完全可以运用其掌握的软件编程技能实现本发明。凡本发明申请文件提及的均属此范畴，申请人不再一一列举。

下面结合附图与具体实施方式对本发明作进一步详细描述：

根据微粒和金属板材的对称性，在ABAQUS/CAE中建立包括四个微粒、一个金属板材构件的高速冲击四分之一对称模型，如图1所示，微粒将依次撞击金属板材的四个拐角区域。金属板材的大小为0.5×0.5×2.0mm，球形微粒的直径为1.0mm。固定金属板材的底面，在其四个侧面分别施加相应的对称约束。微粒和金属面之间设置库伦摩擦，摩擦系数为0.3。金属材料为AISI348，密度为7800kg/m³，弹性模量E＝210GPa，泊松比μ＝0.3。采用减缩积分三维八节点实体单元C3D8R对金属板材和微粒进行离散，金属板材的单元尺寸为0.05mm，微粒的单元尺寸为0.03mm。微粒材料为二氧化硅，其强度和硬度都高于金属板材，因此微粒的变形可以忽略不计，将微粒模拟成刚体，参考点位于其球心，施加初始速度v＝100m/s在参考点上。本实例中所用材料参数如表1所示。

表1实例中所用的材料参数

建立材料AISI348的本构关系如下：

Δσ＝D：Δε^e

其中，Δσ为应力增量张量，D为弹性模量张量，Δε^e为弹性应变增量张量；

建立塑性模型，具体建立方式为：

考虑塑性应变强化和马氏体相变强化的流动准则为：

其中，f为屈服方程；σ_e为Mises等效应力；所述K和分别为各向同性强化应力和等效塑性应变；所述σ₁₁，σ₂₂，σ₃₃，σ₁₂，σ₂₃和σ₃₁为六个方向的应力分量。

各向同性强化应力：

k_w＝470*w^0.97

其中，k_ε为塑性应变强化应力，k_w为马氏体相变强化应力，为等效塑性应变率，为参考应变率，w为马氏体的百分比含量，其增量方程由Santacreu模型给出：

其中，为应力三轴度。

利用ABAQUS/EXPLICIT计算模拟微粒冲击金属板材的过程，如图2利用用户子程序VUMAT判断材料点是否进入塑性，运用应力补偿更新算法实现提出的率相关弹塑性本构模型，计算得到马氏体相变含量。具体过程如下：通过用户子程序VUMAT的SDV定义第n+1 增量步的开始时的初始状态变量w_n值，同时也是第n步结束时的状态变量值，在第n+1增量步开始时读入；其中，所述n是指第n增量步，w_n是指第n增量步结束时的金属板材内的马氏体的百分比含量。

VUMAT由应变增量驱动，首先假设应变增量全为弹性，由公式Δσ＝D：Δε^e计算试探应力：将试探应力代入到屈服准则该公式是步骤(2.2)中的公式在第n+1步的特定计算；

如果f_n+1≤0，则材料尚处于弹性变形阶段，将试探应力更新为n+1增量步应力，马氏体相含量更新为n+1增量步状态变量w_n+1；如果f_n+1＞0，则材料屈服，根据应力补偿显式算法更新增量步结束时的应力，即在试探应力的基础上减去多算的塑性应变增量对应的应力增量部分，对试探应力进行负补偿；

等效塑性应变增量运用牛顿迭代法求解：

其中，所述i表示第i次迭代，表示第i次迭代所得的等效塑性应变增量；所述i-1 表示第i-1次迭代，表示第i-1次迭代所得的等效塑性应变增量；更新等效塑性应变增量，直至满足精度要求，得到马氏体相的百分比含量w_n+1，更新第n+1增量步应力，更新内能，塑性功。

最后，将ABAQUS/EXPLICIT中得到的动态应力状态导入到ABAQUS/STANDARD模块中进行隐式分析，确定静态平衡下金属板材的残余应力场。

图3为沿金属板材沿厚度方向的残余应力分布，随着深度的增加残余压应力先增大后减小；图4为金属板材表面中心点的马氏体相变含量随时间的变化，可见马氏体相变含量随着冲击次数的增加而增多。所以本发明提出的率相关弹塑性本构模型可以准确预测微粒流高速冲击下金属构件内部的残余应力和马氏体相变。

本发明提出的率相关弹塑性本构模型同时考虑了塑性变形和材料相变的影响，结合该本构模型，利用FORTRAN语言编写子程序并嵌入到ABAQUS中，能准确预测微粒高速冲击下金属构件残余应力和马氏体相变，为深入阐明残余应力产生和马氏体相变的机理提供了技术支撑，推进微粒流冲击在提高关键材料和零部件疲劳性能中的应用。

最后，需要注意的是，以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有很多变形。本领域的普通技术人员能从本发明公开的内容中直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (4)

1.一种预测微粒冲击金属件残余应力和马氏体相变的有限元方法，其特征在于，包括以下步骤：

(1)建立微粒冲击金属构件的有限元模型；

(2)建立综合考虑了塑性应变强化和马氏体相变强化的率相关弹塑性本构模型；

所述步骤(2)具体包括下述步骤：

步骤(2.1)建立金属材料的本构关系：

Δσ＝D：Δa^e

其中，Δσ为应力增量张量，D为弹性模量张量，Δa^e为弹性应变增量张量；

步骤(2.2)建立塑性模型，具体建立方式为：

考虑塑性应变强化和马氏体相变强化的流动准则为：

其中，f为屈服方程；σ^e为Mises等效应力；所述K和分别为各向同性强化应力和等效塑性应变；所述σ₁₁，σ₂₂，σ₃₃，σ₁₂，σ₁₃和σ₃₁为六个方向的应力分量；

各向同性强化应力：

k_w＝h_wwⁿ

其中，k_ε为塑性应变强化应力，k_w为马氏体相变强化应力，C为材料参数，为应变率，为参考应变率，σ₀为初始屈服应力，h_ε为应变强化系数，a为材料参数，h_w为相变强化模量，n为强化指数，w为马氏体的百分比含量，其增量方程由Santacreu模型给出：

其中，w_max为马氏体相的最大含量，D₀，D₁和m为材料参数，为应力三轴度；

(3)编写ABAQUS-VUMAT用户动态材料子程序，运用应力补偿更新算法实现提出的率相关弹塑性本构模型；

(4)将子程序嵌入到ABAQUS中，对微粒流冲击金属构件的过程进行计算，进一步获得金属构件内部的残余应力以及马氏体相的百分比含量。

2.根据权利要求1所述的一种预测微粒冲击金属件残余应力和马氏体相变的有限元方法，其特征在于，所述步骤(1)具体为：

基于ABAQUS建立单个或多个球形微粒、金属构件，分别设置材料属性和划分网格，再进行ASSEMBLY对其组装之后设置分析步和通用接触属性。

3.根据权利要求1所述的一种预测微粒冲击金属件残余应力和马氏体相变的有限元方法，其特征在于，所述步骤(3)具体包括下述步骤：

步骤(3.1)：通过用户子程序VUMAT的SDV定义第n+1增量步开始时的初始状态变量w_n值，同时也是第n步结束时的状态变量值，在第n+1增量步开始时读入；其中，所述n是指第n增量步，w_n是指第n增量步结束时的金属构件马氏体的百分比含量；

步骤(3.2)：VUMAT由应变增量驱动，首先假设应变增量全为弹性，由步骤(2.1)中公式Δσ＝D：Δε^e计算试探应力：将试探应力代入到屈服准则中，该公式是步骤(2.2)中的公式在第n+1步的特定计算；

如果f_n+1≤0，则材料尚处于弹性变形阶段，将试探应力更新为n+1增量步应力，马氏体相含量更新为n+1增量步状态变量w_n+1；

如果f_n+1＞0，则材料屈服，根据应力补偿显式算法更新增量步结束时的应力，即在试探应力的基础上减去多算的塑性应变增量对应的应力增量部分，对试探应力进行负补偿；

等效塑性应变增量运用牛顿迭代法求解：

其中，所述i表示第i次迭代，表示第i次迭代所得的等效塑性应变增量；所述i-1表示第i-1次迭代，表示第i-1次迭代所得的等效塑性应变增量；更新等效塑性应变增量，直至满足精度要求，得到马氏体相的百分比含量w_n+1，更新第n+1增量步应力，更新内能，塑性功。

4.根据权利要求1所述的一种预测微粒冲击金属件残余应力和马氏体相变的有限元方法，其特征在于，所述步骤(4)具体为：

将步骤(1)建立的模型主文件和步骤(3)建立的ABAQUS-VUMAT用户子程序联合，使用ABAQUS/EXPLICT方法对微粒高速冲击金属构件进行计算，获得马氏体相变含量；将ABAQUS/EXPLICIT中得到的动态应力状态导入到ABAQUS/STANDARD模块中进行隐式分析，确定静态平衡下的残余应力场，即完成微粒流高速冲击下金属构件内残余应力和马氏体相变的预测。