CN110222442A - 面心立方材料疲劳过程晶体塑性本构模型建立方法 - Google Patents

Abstract

本发明公开了一种面心立方材料疲劳过程晶体塑性本构模型建立方法。本发明解决了在计算材料损伤和疲劳失效过程中无法考虑介观尺度塑性演化和材料微观结构特征问题，提供了更准确的革新的介观晶体建模方法和晶体塑性本构模型及相关稳定高效的计算方法。本发明提出的晶体塑性本构模型与计算方法考虑了晶体滑移面法线方向和滑移方向，利用运动和等向硬化准则考虑背应力效应，最后建立晶体滑移和位错密度与塑性剪切应变速率函数关系的面心立方金属材料晶体塑性本构模型。可以从介观尺度研究材料微观特征对材料损伤与疲劳寿命的影响作用和贡献程度，为研究材料损伤演化和疲劳失效过程提供重要的理论指导和技术支撑，具有重要的科学意义和工程应用价值。

Description

面心立方材料疲劳过程晶体塑性本构模型建立方法

技术领域

本发明涉及材料疲劳失效技术领域，具体涉及一种面心立方材料疲劳过程晶体塑性本构模型建立方法。

背景技术

可靠性与长寿命是材料服役过程的核心问题。研究表明，疲劳是机械零部件和工业设备失效的主要原因之一。统计分析显示，超过80％的机械失效都是由疲劳引起的。疲劳失效具有潜伏性和突然性，并且对相邻零部件也会造成一定的影响。因此，疲劳失效是一种非常复杂的材料-力学相互影响的物理现象，通常材料的失效都是从晶粒或晶界处开始，经过一段塑性变形后在滑移带上形成可见的微观裂纹，大量微观裂纹相互融和形成长裂纹继续在材料内部进行扩展，因此整个疲劳失效过程包括裂纹萌生、裂纹扩展和最终断裂三个阶段,如图1所示。

疲劳裂纹萌生与材料微观结构有关。疲劳过程中，在材料局部区域会逐渐形成驻留滑移带，并沿着高角度晶粒边界从表面侵入挤出，造成该位置位错累积并伴随应力集中和能量升高，产生裂纹源。疲劳裂纹扩展可以用裂纹扩展率(da/dN)和应力强度范围(ΔK)的双对数函数来描述，一般整个曲线分为三个不同区域：1)门槛值区域；2)裂纹稳定扩展区域(Paris)；3)最终断裂区域，其中，裂纹扩展门槛值ΔK_th是描述材料疲劳失效的重要参数，可以定量表征长裂纹是否扩展，也是损伤容限设计中主要指标。研究表明，裂纹扩展门槛值也受材料微观结构的直接影响。

然而，传统的材料可靠性与疲劳寿命研究手段仅仅是试验测试，并且受限于试验环境和成本的限制往往很难对疲劳失效过程进行有效监测与研究。因此，急需一种高效可靠的研究方法或手段进行材料可靠性与疲劳寿命相关问题的研究。随着计算机技术的迅速发展和数值算法的提升，建立材料本构模型并开发对应的迭代算法进行材料疲劳失效行为研究成为一种可靠且可行的试验替代方法。但是，目前固体力学中所建立的材料本构模型只是从宏观层次来考虑各种材料外部因素对疲劳失效行为的影响，例如，载荷大小、加载方式与方向等，基本没有涉及从材料微观角度揭示疲劳失效机理，例如，疲劳裂纹萌生与扩展等，严重影响到对材料损伤和失效机理的有效研究及工程化应用。

发明内容

针对现有技术中的上述不足，本发明提供的一种面心立方材料疲劳过程晶体塑性本构模型建立方法解决了在计算材料损伤和疲劳失效过程中无法考虑介观塑性演化和材料微观结构特征问题的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种面心立方材料疲劳过程晶体塑性本构模型建立方法，包括以下步骤：

S1、通过面心立方材料的EBSD信息建立介观模型；

S2、在介观模型基础上建立材料的晶体塑性本构模型。

进一步地：所述步骤S1的具体步骤为：

S11、根据面心立方材料的EBSD信息利用图像术处理技术在空间生成晶粒晶核坐标信息；

S12、利用Voronoi算法将晶粒晶核坐标信息生成二维或三维晶体模型图；

S13、将晶体模型图信息导出并按一定顺序保存在文本文件中，通过Abaqus软件将文本文件生成晶体二维或三维几何模型；

S14、根据二维或三维几何模型的几何特征在Abaqus软件中进行结构/非结构网格划分得到网格，并对晶粒边界进行网格局部加密细化；

S15、将EBSD中的晶粒取向信息输入到Abaqus软件的网格中，完成考虑晶粒取向的介观模型建立；

S16、当介观模型验证合格时，输出该介观模型，否则，根据面心立方材料的EBSD信息修改晶粒晶核坐标，并返回步骤S11。

进一步地：所述步骤S11中晶粒晶核坐标信息的计算公式为：

上式中，I^g(x,y)为晶粒晶核坐标信息，M和N分别为所选微观区域x和y坐标最大值，p为图片像素值，Ω_mic为微观区域。

进一步地：所述二维晶体模型图为：

上式中，D_i和B_ij分别为二维几何区域和几何边界位置，i和j分别为x和y两个方向的变量；

所述三维晶体模型图为：

上式中，V(p_i)和B'_ij分别为三维集合区域与几何边界位置，i、j和k分别为x、y和z三个方向的变量，I_n为y方向边界位置最大值。

进一步地：所述步骤S2的具体步骤为：

S21、通过面心立方金属材料晶体特征，在介观模型的基础上基于晶体特征和计算塑性理论分别建立滑移和位错密度与塑性剪切应变速率函数关系的本构模型；

S22、获得并定义本构模型的材料和滑移系相关参数；

S23、根据材料参数建立局部坐标下的弹性矩阵，通过取向矩阵对局部坐标下的弹性矩阵进行全局转化，并通过全局弹性矩阵计算滑移系强度；

S24、判断初始滑移系强度是否为0，当为0时通过滑移系相关参数计算滑移系方向和滑移系数量，并通过滑移系方向和滑移系数量计算施密特因子对称部分，否则，直接调用已有的施密特因子对称部分；

S25、通过滑移系相关参数对称部分计算施密特因子反对称部分，并通过施密特因子对称部分和施密特因子反对称部分计算滑移系算子；

S26、通过滑移系相关参数和滑移系算子计算滑移系剪切应变率；

S27、通过滑移系剪切应变率计算公式获得n元显示线性方程迭代计算显示线性方程的剪应变率，然后进行对剪切应变率进行非线性近似求解，获得第n步的剪切应变增量，n为迭代法中线性方程的个数；

S28、迭代计算第n+1步的剪切应变增量；

S29、利用第n+1步的剪切应变增量迭代计算一致切线刚度矩阵，并判断其是否收敛，若是，则结束迭代计算，更新本构模型中的相关参数，建立晶体塑性本构模型，否则，返回步骤S27。

进一步地：所述步骤S25滑移系算子的计算公式为：

C＝P^a+W^aσ+σW^a

上式中，C为滑移系算子，P^a为施密特因子对称部分，W^a为施密特因子反对称部分，σ为应力。

进一步地：所述步骤S26中滑移系剪切应变率的计算公式为：

上式中，为参考应变率，每个滑移系均相同；τ^α为第α滑移系的分解切应力；χ^α为第α滑移系的背应力或随动硬化阻力；g^α为第α滑移系各向同性硬化阻力；m是控制材料粘性的应变率敏感指数。

进一步地：所述步骤S27中剪切应变增量的计算公式为：

上式中，Δγ^α为剪切应变增量，Δt为时间增量，为t时刻剪切应变率，为t+1时刻剪切应变率，θ为积分算法控制参数。

进一步地：所述步骤S29中一致切线刚度矩阵的计算公式为：

上式中，C^alg为一致切线刚度矩阵，l_ijkl为弹性模量张量，为Schmid因子，为滑移方向与滑移面法向叉乘之和，σ_jk为应力张量，为滑移方向与滑移面法向叉乘之和，σ_ik为应力张量，γ^α为剪切应变增量，Δε为应变增量，a为滑移系，N为滑移系最大数量。

本发明的有益效果为：本发明解决了在计算材料损伤和疲劳失效过程中无法考虑介观塑性演化和材料微观结构特征问题，提供了更准确的革新的介观晶体建模方法和晶体塑性本构模型及相关稳定高效的计算方法。本发明提出的介观晶界模型是基于材料真实微观结构特征，利用Voronoi算法和Matlab图像处理技术来实现，具有良好的兼容性和可移植性。同时，本发明提出的晶体塑性本构模型与计算方法考虑了晶体滑移面法线方向和滑移方向，利用运动和等向硬化准则考虑背应力效应，最后建立晶体滑移和位错与塑性剪切应变速率函数关系的面心立方金属材料晶体塑性本构模型。相关本构模型与计算方法的建立，可以从介观角度研究材料微观特征对材料损伤与疲劳寿命的影响作用和贡献程度，并建立相互作用的定量表达式，为研究材料损伤演化和疲劳失效过程提供重要的理论指导和技术支撑，具有重要的科学意义和工程应用价值。

附图说明

图1为本发明背景技术中疲劳失效不同阶段定义示意图；

图2为本发明流程图；

图3为本发明中的介观尺度模型示意图；

图4为本发明中的介观尺度模型塑性示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图2所示，一种面心立方材料疲劳过程晶体塑性本构模型建立方法，包括以下步骤：

S1、通过面心立方材料的EBSD(电子背散射衍射)信息建立介观模型；具体步骤为：

S11、根据面心立方材料的EBSD信息利用图像术处理技术在空间生成晶粒晶核坐标信息；晶粒晶核坐标信息的计算公式为：

上式中，I^g(x,y)为晶粒晶核坐标信息，M和N分别为所选微观区域x和y坐标最大值，p为图片像素值，Ω_mic为微观区域。

S12、利用Voronoi算法将晶粒晶核坐标信息生成二维或三维晶体模型图；所述二维晶体模型图为：

上式中，D_i和B_ij分别为二维几何区域和几何边界位置，i和j分别为x和y两个方向的变量；

所述三维晶体模型图为：

上式中，V(p_i)和B'_ij分别为三维集合区域与几何边界位置，i、j和k分别为x、y和z三个方向的变量，I_n为y方向边界位置最大值。

S13、将晶体模型图信息导出并按一定顺序保存在文本文件中，通过Abaqus软件将文本文件生成晶体二维或三维几何模型；

S14、根据二维或三维几何模型的几何特征在Abaqus软件中进行结构/非结构网格划分得到网格，并对晶粒边界进行网格局部加密细化；

S15、将EBSD中的晶粒取向信息输入到Abaqus软件的网格中，完成考虑晶粒取向的介观模型建立；

S16、当介观模型验证合格时，输出该介观模型，否则，根据面心立方材料的EBSD信息修改晶粒晶核坐标，并返回步骤S11。

S2、在介观模型的基础上建立材料的晶体塑性本构模型。具体步骤为：

S21、通过面心立方金属材料晶体特征，在介观模型的基础上基于晶体特征和计算塑性理论分别建立滑移和位错密度与塑性剪切应变速率函数关系的本构模型；

S22、获得并定义本构模型的材料和滑移系相关参数；

S23、根据材料参数建立局部坐标下的弹性矩阵，通过取向矩阵对局部坐标下的弹性矩阵进行全局转化，并通过全局弹性矩阵计算滑移系强度；

S24、判断初始滑移系强度是否为0，当为0时通过滑移系相关参数计算滑移系方向和滑移系数量，并通过滑移系方向和滑移系数量计算施密特因子对称部分，否则，直接调用已有的施密特因子对称部分；

S25、通过滑移系相关参数对称部分计算施密特因子反对称部分，并通过施密特因子对称部分和施密特因子反对称部分计算滑移系算子；滑移系算子的计算公式为：

C＝P^a+W^aσ+σW^a

上式中，C为滑移系算子，P^a为施密特因子对称部分，W^a为施密特因子反对称部分，σ为应力。

S26、通过滑移系相关参数和滑移系算子计算滑移系剪切应变率；滑移系剪切应变率的计算公式为：

上式中，为参考应变率，每个滑移系均相同；τ^α为第α滑移系的分解切应力；χ^α为第α滑移系的背应力或随动硬化阻力；g^α为第α滑移系各向同性硬化阻力；m是控制材料粘性的应变率敏感指数。

S27、通过滑移系剪切应变率计算公式获得n元显示线性方程迭代计算显示线性方程的剪应变率，然后进行对剪切应变率进行非线性近似求解，获得第n步的剪切应变增量，n为迭代法中线性方程的个数；

剪切应变增量的计算公式为：

上式中，Δγ^α为剪切应变增量，Δt为时间增量，为t时刻剪切应变率，为t+1时刻剪切应变率，θ为积分算法控制参数。

S28、迭代计算第n+1步的剪切应变增量；

S29、利用第n+1步的剪切应变增量迭代计算一致切线刚度矩阵，并判断其是否收敛，若是，则结束迭代计算，更新本构模型中的相关参数，建立晶体塑性本构模型，否则，返回步骤S27。

一致切线刚度矩阵的计算公式为：

上式中，C^alg为一致切线刚度矩阵，l_ijkl为弹性模量张量，为Schmid因子，为滑移方向与滑移面法向叉乘之和，σ_jk为应力张量，为滑移方向与滑移面法向叉乘之和，σ_ik为应力张量，γ^α为剪切应变增量，Δε为应变增量，a为滑移系，N为滑移系最大数量。

根据本发明建立的二维几何和网格模型分别如图3(a)和(b)所示。图3(a)为介观尺度几何模型，图3(b)为介观尺度网格模型图。

面向面心立方材料的晶体塑性算法，在介观尺度模型上来进行塑性演化相关机理的计算，结果如图4所示。从图4中可以看出基于介观模型和晶体塑性算法，在材料受载过程中可以很好地预测出晶粒边界的塑性应力分布情况，同时也可以从介观尺度了解到塑性的演化过程，从而为材料失效机理的研究提供理论依据。

Claims (9)

1.一种面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，包括以下步骤：

S1、根据面心立方材料的EBSD信息建立介观模型；

S2、在介观模型的基础上建立材料的晶体塑性本构模型。

2.根据权利要求1所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述步骤S1的具体步骤为：

S11、根据面心立方材料的EBSD信息利用图像术处理技术在空间生成晶粒晶核坐标信息；

S12、利用Voronoi算法将晶粒晶核坐标信息生成二维或三维晶体模型图；

S13、将晶体模型图信息导出并按一定顺序保存在文本文件中，通过Abaqus软件将文本文件生成晶体二维或三维几何模型；

S14、根据二维或三维几何模型的几何特征在Abaqus软件中进行结构/非结构网格划分得到网格，并对晶粒边界进行网格局部加密细化；

S15、将EBSD中的晶粒取向信息输入到Abaqus软件的网格中，完成考虑晶粒取向的介观模型建立；

S16、当介观模型验证合格时，输出该介观模型，否则，根据面心立方材料的EBSD信息修改晶粒晶核坐标，并返回步骤S11。

3.根据权利要求2所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述步骤S11中晶粒晶核坐标信息的计算公式为：

上式中，I^g(x,y)为晶粒晶核坐标信息，M和N分别为所选微观区域x和y坐标最大值，p为图片像素值，Ω_mic为微观区域。

4.根据权利要求2所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述二维晶体模型图为：

上式中，D_i和B_ij分别为二维几何区域和几何边界位置，i和j分别为x和y两个方向的变量；

所述三维晶体模型图为：

上式中，V(p_i)和B'_ij分别为三维集合区域与几何边界位置，i、j和k分别为x、y和z三个方向的变量，I_n为y方向边界位置最大值。

5.根据权利要求1所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述步骤S2的具体步骤为：

S21、通过面心立方金属材料晶体特征，在介观模型的基础上基于晶体特征和计算塑性理论分别建立滑移和位错密度与塑性剪切应变速率函数关系的本构模型；

S22、获得并定义本构模型的材料和滑移系相关参数；

S23、根据材料参数建立局部坐标下的弹性矩阵，通过取向矩阵对局部坐标下的弹性矩阵进行全局转化，并通过全局弹性矩阵计算滑移系强度；

S24、判断初始滑移系强度是否为0，当为0时通过滑移系相关参数计算滑移系方向和滑移系数量，并通过滑移系方向和滑移系数量计算施密特因子对称部分，否则，直接调用已有的施密特因子对称部分；

S25、通过滑移系相关参数对称部分计算施密特因子反对称部分，并通过施密特因子对称部分和施密特因子反对称部分计算滑移系算子；

S26、通过滑移系相关参数和滑移系算子计算滑移系剪切应变率；

S27、通过滑移系剪切应变率计算公式获得n元显示线性方程迭代计算显示线性方程的剪应变率，然后进行对剪切应变率进行非线性近似求解，获得第n步的剪切应变增量，n为迭代法中线性方程的个数；

S28、迭代计算第n+1步的剪切应变增量；

S29、利用第n+1步的剪切应变增量迭代计算一致切线刚度矩阵，并判断其是否收敛，若是，则结束迭代计算，更新本构模型中的相关参数，建立晶体塑性本构模型，否则，返回步骤S27。

6.根据权利要求5所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述步骤S25滑移系算子的计算公式为：

C＝P^a+W^aσ+σW^a

上式中，C为滑移系算子，P^a为施密特因子对称部分，W^a为施密特因子反对称部分，σ为应力。

7.根据权利要求6所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述步骤S26中滑移系剪切应变率的计算公式为：

上式中，为参考应变率，每个滑移系均相同；τ^α为第α滑移系的分解切应力；χ^α为第α滑移系的背应力或随动硬化阻力；g^α为第α滑移系各向同性硬化阻力；m是控制材料粘性的应变率敏感指数。

8.根据权利要求5所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述步骤S27中剪切应变增量的计算公式为：

上式中，Δγ^α为剪切应变增量，Δt为时间增量，为t时刻剪切应变率，为t+1时刻剪切应变率，θ为积分算法控制参数。

9.根据权利要求5所述的面心立方材料疲劳过程晶体塑性本构模型建立方法，其特征在于，所述步骤S29中一致切线刚度矩阵的计算公式为：

上式中，C^alg为一致切线刚度矩阵，l_ijkl为弹性模量张量，为Schmid因子，为滑移方向与滑移面法向叉乘之和，σ_jk为应力张量，为滑移方向与滑移面法向叉乘之和，σ_ik为应力张量，γ^α为剪切应变增量，Δε为应变增量，a为滑移系，N为滑移系最大数量。