CN110765688B - 一种基于abaqus的异质结构模拟方法 - Google Patents

Abstract

本发明提供一种基于ABAQUS的异质结构模拟方法，包括：建立均匀结构和异质结构的2D代表性体积单元；将滑移阻力模型和背应力模型修正并写入用户子程序UMAT，以定义晶体的弹塑形本构方程；利用三种不同结构的拉伸曲线对UMAT中的材料参数进行确定；施加周期性边界条件和外部载荷，利用ABAQUS计算每个晶粒的应力应变响应；对所有晶粒的应力应变平均化；改变模型微观结构和晶粒取向等参数，预测和分析对应的宏观应力应变响应、强化效果和微观变形云图，对比得出最优的微观结构。本发明的异质结构模拟方法考虑了异质结构的晶粒尺寸效应和额外背应力强化，可对异质结构进行模拟，具有直观、适用性强、精确度高的优点。

Description

一种基于ABAQUS的异质结构模拟方法

技术领域

本发明属于材料力学性能有限元数值计算领域，涉及一种基于ABAQUS的异质结构模拟方法。

背景技术

随着我国航空航天领域的快速发展，我们对材料的综合性能提出了更高的要求。大多数金属纳米和超细晶材料虽能满足强度要求，但低延展性限制了其在许多先进技术中的应用。近年来，已报道了在环境温度下具有良好强韧综合性能和在高温下的超塑性的金属异质结构材料，高强度和良好延性的独特组合使得异质结构金属成为结构和功能的“梦想工程材料”，并应用于交通、医药、能源等领域。

为了进一步发展纳米金属材料领域，揭示出可以产生和获得高延性金属异质结构材料的潜在机制是亟需解决的关键技术瓶颈之一。近些年，有限元软件的发展可以很好地满足人们对于复杂应力应变行为的理解以及提供了该状态下不同微观结构的综合力学性能。有限元软件ABAQUS不仅可以分析复杂的固定力学和结构力学系统，还可以通过强大的二次开发接口补充ABAQUS前后处理模块中不完善的功能。基于FORTRAN语言的用户自定义的子程序扩展了ABAQUS在本构方程的应用，并实现了损伤评估、力学性能预测等方面的功能。

然而为了考虑异质结构变形过程中存在的内部物理机制，传统未考虑晶粒尺寸的滑移阻力模型和背应力模型已无法满足要求，其中，现有的滑移模型没办法考虑跨尺度晶粒，背应力模型没办法考虑异质结构的特殊强化。

发明内容

本发明的目的在于提供一种基于ABAQUS的异质结构模拟方法，以使得模拟方法能够模拟跨尺度晶粒，并且直观、适用性强、精确度高。

为了实现上述目的，本发明提供了一种基于ABAQUS的异质结构模拟方法，包括：

S1：建立均匀结构和异质结构的ABAQUS代表性单元模型；

S2：将滑移阻力模型和背应力模型分别修正为多尺度的滑移阻力模型和考虑强度差的背应力模型并写入ABAQUS的用户子程序UAMT，通过用户子程序UMAT定义晶体的弹塑性本构方程，该弹塑性本构方程中包括多个材料参数；

S3：根据两组不同晶粒尺寸的均匀结构和一组异质结构的单轴拉伸试验得到的拉伸曲线的实验结果与所述步骤S2中的弹塑性本构方程进行拟合，对所述步骤S2中的材料参数进行确定；

S4：施加周期性边界条件和外部载荷，利用ABAQUS计算每个晶粒的应力应变响应；

S5：对所有晶粒的应力应变平均化，得到异质结构的宏观力学性能；

S6：改变所述步骤S1中的ABAQUS代表性单元模型的微观结构和晶粒取向，预测和分析对应的宏观应力应变响应、强化效应和微观变形云图，通过对比获得最优的微观结构；

其中，所述步骤S2包括：

S21：采用投影运算方法，求解步骤S1中的异质结构的ABAQUS代表性单元模型的坐标系与晶粒坐标系之间的转化矩阵，基于所述转化矩阵，通过张量运算，得到异质结构的ABAQUS代表性单元模型的坐标系下的施密特因子，随后求解作用在晶粒中的应力张量σ^g，并结合材料弹性张量C^-1，进而计算出每个晶粒的弹性应变；

S22：由热驱动流动准则算出塑性应变，包括：

S221：在滑移阻力模型中引入Hall-Petch准则，将滑移阻力模型中的初始滑移阻力参数s₀修正为多尺度的初始滑移阻力参数S₀；

S222；在背应力模型中材料参数h_b中引入强度差模型，得到考虑强度差的背应力模型；

S223：根据基于热驱动流动准则的弹塑性变形理论求解晶粒分配给每个潜在滑移系的切变量

进而计算每个晶粒的塑性应变：

S23：结合所述步骤S21中的弹性应变和所述步骤S22中的塑性应变，得到每个晶粒的应力应变响应；

在所述步骤S221中，每个滑移系的滑移阻力为：

h^αβ＝h_s[w+(1-w)δ_αβ]，

其中，h_s为静态硬化模量，w代表晶格硬化率，δ_αβ代表克罗内克函数，s₀、S_sat为传统晶粒尺寸对应的初始滑移阻力和饱和滑移阻力，S₀为多尺度的初始滑移阻力参数(多尺度体现在粗细晶区域使用不同的霍尔佩奇参数s₀,k_s，可见实施例中参数拟合)，k_s为霍尔佩奇系数，d为晶粒尺寸；

在所述步骤S222中，每个滑移系的背应力为：

其中，h_b代表硬化率，

代表参考硬化率，k_b、r_D为材料参数，d_n为当前计算晶粒的周围晶粒平均尺寸。

优选地，在所述步骤S223中，每个潜在滑移系的切变量

为：

其中，T为绝对温度，k为玻尔兹曼常数，B^α为第α滑移系的背应力，S^α为第α滑移系的滑移阻力，τ^α为第α滑移系的分解切应力，表示为τ^α＝σ^g:P^α，F₀为热激活位错越过障碍的激活能，τ₀、p、q为材料常数，

为参考剪应变率；

每个晶粒的塑性应变为：

其中，P^α为第α滑移系的方向因子张量，m^α为第α滑移系滑移方向的单位矢量，n^α为其滑移面法向的单位矢量。

优选地，在所述步骤S21中，所述异质结构的ABAQUS代表性单元模型的坐标系下的施密特因子α为：

α＝cosλcosβ，

其中，λ为滑移方向与外力F的夹角，β为滑移面与外力F中心轴的夹角；

所述弹性应变ε^e为：

ε^e＝C^-1:σ^g，

其中，ε^e为弹塑性本构方程中的弹性应变，σ^g为作用在晶粒中的应力张量，C^-1为材料弹性张量。

优选地，所述步骤S1包括：

S11：基于matlab语言根据Voronoi tessellation方法建立均匀结构的2D代表性体积单元，该2D代表性体积单元包括多个晶粒；

S12：根据步骤S11中的均匀结构的2D代表性体积单元的晶粒个数随机生成各个晶粒的取向参数，得到随机取向文件；

S13：通过ABAQUS的脚本程序读取所述步骤S11中的均匀结构的2D代表性体积单元和所述步骤S12中的随机取向文件，生成均匀结构的ABAQUS代表性单元模型，并作为异质结构的均匀细晶模型；

S14：根据需要研究的异质结构的微观结构的形貌图，在ABAQUS中对所述步骤S13中的均匀细晶模型进行手动修正，首先根据所述形貌图中的粗晶形貌将所述均匀细晶模型中的若干个晶粒合并成粗晶并删除内部多余线条，再删除每个合并的晶粒的取向参数并重新赋予粗晶取向，多次对比校正后，形成异质结构的ABAQUS代表性单元模型；

S15：使用ABAQUS的脚本程序计算所述步骤S14中的异质结构的ABAQUS代表性单元模型中每个晶粒的尺寸，并将计算结果赋予相应晶粒。

优选地，在所述步骤S12中，所述晶粒的取向参数采用三个欧拉角

φ、

表示。

优选地，在步骤S3中，所述材料参数是通过结合所述步骤S2中的弹塑性本构方程，采用试参法计算两组不同晶粒尺寸的均匀结构和一组异质结构的ABAQUS代表性单元模型的拉伸曲线的模拟结果，直至模拟结果与单轴拉伸试验所得的拉伸曲线的实验结果拟合，来确定的。

优选地，所述步骤S4包括：结合所述步骤S1中的异质结构的ABAQUS代表性单元模型、步骤S3中确定的材料参数，并施加的周期性边界条件和外部载荷，在ABAQUS中划分模型网格，再结合所述步骤S2中的弹塑性本构方程，通过ABAQUS软件计算得到每个积分点的应力应变张量。

优选地，在步骤S6中，所述最优的微观结构通过最大化的激发微观结构能够诱发的背应力强化，进而以牺牲最小甚至不牺牲塑性的情况下提高强度，来获得。

优选地，所述异质结构包括异质双峰结构、层片结构和梯度结构。

(1)本发明的基于ABAQUS的异质结构模拟方法在现有的热驱动流动准则采用的描述晶体塑性应变的理论机制的基础上，与有限元方法相结合，采用基于晶粒组织的网格模型，并通过霍尔佩奇方程来考虑不同尺度下的晶粒尺寸效应，从而得到了多尺度的滑移阻力模型，能够准确预测不同异质微观结构的弹塑性转变过程和屈服应力，解决了现有模型无法预测不同晶粒尺寸引起的应变强化效应的问题。

(2)本发明的基于ABAQUS的异质结构模拟方法利用ABAQUS的有限元网格的单元和节点化来实现晶粒间的相互作用，并同时通过考虑异质结构中粗晶和细晶的强度差来考虑由于大的应变梯度引起的背应力强化效应，解决了现有模型无法考虑金属材料中由于组成相的强度差异大，产生应变梯度而引起的背应力强化的问题，从而可以更深层次的考虑异质结构变形过程中存在的内部物理机制和强化机制；

(3)本发明的基于ABAQUS的异质结构模拟方法采用有限元软件ABAQUS，其具有很强的直观性，能够直观的提取异质结构的背应力强化数值、变形滑移带和应力应变集中位置，并可充分发挥有限元模拟计算的优势，实现复杂的边界加载条件和多样的可视化后处理操作；采用本发明提供的方法可以计算得到晶粒组织有限元网格模型中每个晶粒的随着加载过程不断变化的应力应变响应，对再现材料组织演变过程、分析材料局部受力状态，以及预测材料失效形式等提供了更为准确的研究方法，具有直观、适用性强、精确度高的优点。

综上所述，本发明的基于ABAQUS的异质结构模拟方法在热驱动流动准则采用的描述晶体塑性应变的理论机制基础上与有限元方法相结合，不仅考虑了异质结构中的晶粒尺寸效应和由于微观结构不均匀引起的背应力强化，同时用有限元网络的单元和节点化来实现晶粒间的相互作用，可对异质结构-存在着跨尺度设计理念的微观结构进行模拟，弥补了传统模型只能模拟均匀晶粒大小微观组织的不足，并可充分发挥有限元模拟计算的优势，实现复杂的边界加载条件和多样的可视化后处理操作，具有直观、适用性强、精确度高的优点。

附图说明

图1为本发明实施例的基于ABAQUS的异质结构模拟方法的流程图；

图2为异质结构的拉伸曲线的实验结果、采用根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法的拉伸曲线的模拟结果以及采用现有的无约束背应力方法的模拟结果的对比图；

图3为室温下不同微观形态的异质双峰结构在采用多种模拟方法时分别得到的屈服应力预测对比图；

图4为根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法所得到的异质双峰结构及其两个尺度的晶粒的应力应变响应图，其示出了结构载荷传递和应变分配等强化机制；

图5为异质双峰结构和均匀结构的背应力强化效应图；

图6为根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法的不同细晶尺寸下异质结构的应力应变响应和应变硬化率对比图；

图7为根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法的不同粗晶尺寸下异质结构的应力应变响应和应变硬化率对比图；

图8为根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法的不同取向下异质结构应力应变响应和应变硬化率对比图。

具体实施方式

下面将结合说明书附图和具体的实施例对本发明所述的一种基于ABAQUS的异质结构模拟方法做进一步的详细说明。

图1为根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法的流程图。其中，异质结构为跨尺度晶粒尺寸设计的结构，包括异质双峰结构，层片结构和梯度结构等，具体是通过冷轧后退火形成局部晶粒再结晶获得的。

如图1所示，在本实施方式中，基于ABAQUS和晶体塑性本构对异质结构材料设计的模拟方法包括步骤：

S1：建立均匀结构和异质结构的ABAQUS代表性单元模型；

其中，在步骤S1具体包括：

S11：基于matlab语言根据Voronoi tessellation方法建立均匀结构的2D代表性体积单元，该2D代表性体积单元包括多个晶粒；

其中，此2D代表性体积单元的数量为一个且施加以周期性边界条件，此2D代表性体积单元的体积尺寸是最小，且2D代表性体积单元内包含了足够多微观组织结构组成物的几何信息、晶体学取向信息、分布信息与相场信息，并能在统计学意义上(统计平均性质)代表材料微观组织结构的基本特征。2D代表性体积单元的晶粒个数和晶粒尺寸与实验结构相吻合或根据用户需要设置。

S12：根据步骤S11中的均匀结构的2D代表性体积单元的晶粒个数并使用matlab程序随机生成2D代表性体积单元的各个晶粒的取向参数，得到随机取向文件；其中，晶粒的取向参数采用三个欧拉角

φ、

表示；

S13：通过ABAQUS的脚本程序读取所述步骤S11中的均匀结构的2D代表性体积单元和所述步骤S12中的随机取向文件，生成均匀结构的ABAQUS代表性单元模型。特别地，对于异质结构，此模型作为均匀细晶模型，在S14中以此为基础合成异质结构中的粗晶，进而形成异质结构的代表性单元模型。

S14：根据需要研究的异质结构的微观结构的形貌图，在ABAQUS中对所述步骤S13中的均匀细晶(小晶粒)模型进行手动修正，首先根据所述形貌图中的粗晶形貌将所述均匀细晶模型中的若干个晶粒合并成粗晶(大晶粒)并删除内部多余线条，再删除每个合并的晶粒已经赋予的取向参数并重新赋予粗晶取向，多次对比校正后，形成异质结构的ABAQUS代表性单元模型。对于均匀结构，则忽略此步骤S14。

S15：使用ABAQUS的脚本程序计算所述步骤S14中的异质结构的ABAQUS代表性单元模型中每个晶粒的尺寸，并将计算结果赋予相应晶粒。

所述晶粒的尺寸d为：

其中，d为晶粒的尺寸，S为晶粒的面积。

S2：将滑移阻力模型和背应力模型分别修正为多尺度滑移阻力模型和考虑强度差的背应力模型并写入ABAQUS的用户子程序UAMT，通过用户子程序UMAT定义晶体的弹塑性本构方程，该弹塑性本构方程中包括多个材料参数，其中，所述材料参数包括三个弹性常数C₁₁,C₁₂,C₄₄,和十二个流动准则参数

F₀,

r_D,h_s,w,τ₀,p,q,k_s,s₀,k_b，各参数含义见下文具体模型。

其中，在步骤S2中，弹塑性本构方程包括：背应力模型的方程、滑移阻力模型的方程和基于热驱动流动准则的方程，步骤S2具体包括以下步骤：

S21：采用投影运算方法，求解所述步骤S1中的异质结构的ABAQUS代表性单元模型的坐标系与晶粒坐标系之间的转化矩阵，基于所述转化矩阵，通过张量运算，得到异质结构的ABAQUS代表性单元模型的坐标系下的施密特因子α，随后求解作用在晶粒中的应力张量σ^g，并结合材料弹性张量C^-1，进而计算出每个晶粒的弹性应变ε^e；

其中，如上文所述，晶粒的取向参数采用三个欧拉角

φ、

定义，即将晶粒坐标系绕其z轴，在xy平面内转过

角；再绕其x轴，在yz平面内转过φ角；最后再一次绕z轴，在xy平面内转过

角，最终与有限元网格模型坐标系完全重合，由此得到每个晶粒的三个空间取向欧拉角。

因此，所述步骤S1中的异质结构的ABAQUS代表性单元模型的坐标系与晶粒坐标系之间的转化矩阵为：

因此，ABAQUS代表性单元模型的坐标系与晶粒坐标系之间的转换矩阵位置关系可表示为：C_c＝gC_s，

其中，C_c为晶粒坐标系，C_s为ABAQUS代表性单元模型的坐标系。

对于张量运算，作用在晶粒的Cauchy应力张量分量可表示为：σ_ij＝F_ijA_i ^-1(i,j＝1,2,3),

其中，F_ij为i方向所受的沿X_j方向的力，A_i为当前晶粒X_i的横截面积，X₁，X₂和X₃方向为ABAQUS代表性单元模型坐标系的三个坐标轴。

异质结构的ABAQUS代表性单元模型的坐标系下的施密特因子α为：

α＝cosλcosβ (1)；

其中，λ为滑移方向与外力F的夹角，β为滑移面与外力F中心轴的夹角。

所述弹性应变ε^e为：

ε^e＝C^-1:σ^g (2)；

其中，ε^e为弹塑性本构方程中的弹性应变，σ^g为作用在晶粒中的应力张量，C^-1为材料弹性张量。

S22：由热驱动流动准则算出塑性应变，具体包括：

S221：在滑移阻力模型中引入Hall-Petch准则，将滑移阻力模型中的参数S₀修正为考虑晶粒尺寸的多尺度的初始滑移阻力参数，因此，S₀为对应晶粒尺寸的初始滑移阻力；

在热驱动流动滑移准则中，首先需要对现有技术中的滑移模型进行修正。由于异质结构例如双峰结构中存在着跨尺度的结构设计，所以必须考虑小尺寸晶粒的高密度晶界对位错运动的阻碍作用，因此，在滑移阻力模型中引入Hall-Petch准则以考虑不同晶粒尺度对应的不同强度，将滑移阻力模型中的初始滑移阻力参数s₀修正为考虑晶粒尺寸的多尺度的初始滑移阻力参数S₀，具体包括：

异质结构中，每个滑移系的滑移阻力为：

h^αβ＝h_s[w+(1w)δ_αβ] (4)；

其中，h_s为静态硬化模量，w代表晶格硬化率，δ_αβ代表克罗内克函数，s₀、S_sat为传统晶粒尺寸对应的初始滑移阻力和饱和滑移阻力，S₀为多尺度的初始滑移阻力参数(多尺度体现在粗细晶区域使用不同的霍尔佩奇参数s₀,k_s，可见实施例中参数拟合)，k_s为霍尔佩奇系数，d为晶粒尺寸；

S222：在背应力模型中材料参数h_b中引入强度差模型，得到考虑强度差的背应力模型；

由于异质结构例如双峰结构的粗细晶强度差，在界面处存在应变梯度，需要大量的几何必须位错来容纳，引发了额外的背应力强化，所以在背应力模型中考虑晶粒尺寸来描述这种额外的强化，具体为在背应力模型的材料参数h_b中，引入强度差模型。

在异质结构中，每个滑移系的背应力为：

其中，h_b代表硬化率，

代表参考硬化率，k_b、r_D为材料参数，d_n为当前计算晶粒的周围晶粒平均尺寸；

S223：根据基于热驱动流动准则的弹塑性变形理论求解晶粒分配给每个潜在滑移系的切变量

进而计算每个晶粒的塑性应变：

每个潜在滑移系的切变量

为：

其中，T为绝对温度，k为玻尔兹曼常数，B^α为第α滑移系的背应力，S^α为第α滑移系的滑移阻力，τ^α为第α滑移系的分解切应力，表示为τ^α＝σ^g:P^α，F₀为热激活位错越过障碍的激活能，τ₀、p、q为材料常数，

为参考剪应变率。符号运算<.>表示：当x≥0时，<x>≡x；当x≤0.时，<x>≡0，

每个晶粒的塑性应变为：

其中，P^α为第α滑移系的方向因子张量，m^α为第α滑移系滑移方向的单位矢量，n^α为其滑移面法向的单位矢量。

S23：结合所述步骤S21中的弹性应变和所述步骤S22中的塑性应变，得到每个晶粒的应力应变响应；

S3：根据两组不同晶粒尺寸的均匀结构和一组异质结构一共三种不同微观结构的单轴拉伸试验得到的拉伸曲线的实验结果与所述步骤S2中的弹塑性本构方程进行拟合，对所述步骤S2中的材料参数进行确定；

在步骤S3中，所述材料参数是通过结合所述步骤S2中的弹塑性本构方程，采用试参法计算两组不同晶粒尺寸的均匀结构和一组异质结构的ABAQUS代表性单元模型的拉伸曲线的模拟结果，直至模拟结果与单轴拉伸试验所得的拉伸曲线的实验结果拟合，来确定的。

S4：施加周期性边界条件和外部载荷，利用ABAQUS计算每个晶粒的应力应变响应；

所述步骤S4具体包括：结合所述步骤S1中的异质结构的ABAQUS代表性单元模型、步骤S3中确定的材料参数，并施加的周期性边界条件和外部载荷，在ABAQUS中划分模型网格，再结合所述步骤S2中的弹塑性本构方程，通过ABAQUS软件计算得到每个积分点的应力应变张量，

S5：利用fortran脚本对所有晶粒的应力应变平均化，得到异质结构的宏观力学性能；

S6：改变所述步骤S1中的ABAQUS代表性单元模型的微观结构和晶粒取向，预测和分析对应的宏观应力应变响应、强化效应和微观变形云图，通过对比获得最优的微观结构，进而指导材料设计。

在步骤S6中，所述最优的微观结构通过最大化的激发微观结构能够诱发的背应力强化，进而以牺牲最小甚至不牺牲塑性的情况下提高强度，来获得。

实验结果

为了更好地阐述本案的基于ABAQUS的异质结构模拟方法的预测效果，下面将采用低碳钢异质双峰结构作为待测的异质结构进行验证。

验证所采用的低碳钢异质双峰结构为典型的异质结构，具有晶粒尺寸差距较大的粗晶和细晶组成，由于晶粒尺寸的差异，能够引起较强的背应力强化，促使强度和延性综合提升。

图2为异质结构的拉伸曲线的实验结果、采用根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法的拉伸曲线的模拟结果以及采用现有的无约束背应力方法的模拟结果的对比图。

如图2所示，通过试参法调整单轴拉伸预测曲线和实验曲线能够较好地吻合，其中，对用户子程序UMAT中的材料参数确定如下，此实施例中，定义晶粒尺寸小15μm的为细晶，细晶参数：E＝200GPa,C₁₁＝269.2MPa,C₁₂＝115.4MPa,C₄₄＝76MPa v＝0.3,

F₀＝286kJ mol^-1,T＝293K,p＝1.0,q＝2.0,τ₀＝90MPa,s₀＝2MPa,h_s＝160MPa,k_s＝18.25MPa.μm^-1/2,w＝2,r_D＝6MPa,

k_b＝27MPa.μm^-1/2，晶粒尺寸大15μm的定义为粗晶，粗晶参数：s₀＝28.9MPa,k_s＝9.346MPa.μm^-1/2，其他参数与细晶相同。

通过得到的参数对异质双峰结构进行材料设计，能够获得如下结果：

1、拉伸曲线的演化。

图2示出了室温下两种微观形态的均匀结构和异质双峰结构的拉伸曲线的实验结果和多种模拟结果的对比，可以看出忽略约束导致的背应力强化会明显低估异质双峰结构的流动应力。

2、异质双峰结构的屈服应力预测

图3为室温下不同微观形态的异质双峰结构在采用多种模拟方法时分别得到的屈服应力预测对比图，可以看出使用基于数值方法的混合法则对异质双峰结构屈服强度的预测误差带在25％之内。而使用现阶段广泛使用的基于统一尺度霍尔佩奇公式的滑移阻力对屈服强度的预测的误差带在20％之内。而使用本方法，则能将异质双峰结构的屈服强度预测误差带控制在10％之内，证明此数值模拟方法具有较高的可靠性。

3、异质双峰结构中两个尺度的晶粒应力应变响应

本模型支持将异质结构中不同范围尺寸的晶粒应力应变状态分别提取以供分析，图4为根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法所模拟得到的异质双峰结构及其两个尺度的晶粒的应力应变响应图，从图中可以体现异质结构的许多强化机制，例如应变梯度强化，背应力强化，载荷传递和应变分配。软相提前进入屈服阶段，但是此时仍比处于弹性阶段的硬相更容易变形，由于协调变形要求，此时软相虽受硬相约束又需承担更多正向应变，几何必须位错快速形成塞积群，由此在向界面产生长程背应力。在共变形阶段，两相间仍具有一定塑性不相容性，因而该过程中依然存在应力和应变分配，仅是分配因子q＝（σ_hard-σ_soft)/ε_soft-ε_hard)逐渐下降，如图4所示。

4、异质双峰结构背应力强化数值

图5是异质双峰结构的背应力强化效应图。从图中可以体现，在整体应变为13％时，异质双峰结构的背应力相比于普通均匀结构多出30％左右，双峰结构比均匀结构多出的额外强化效应基本都是由于粗细晶相互作用,也就是约束诱导产生的背应力主导的。

5、不同微观形态的应变硬化率预测

图6和图7预测了不同微观形态的异质结构拉伸力学性能，其中，图6示出了不同细晶尺寸下异质结构的应力应变响应和应变硬化率，图7示出了不同粗晶尺寸下异质结构的应力应变响应和应变硬化率，可以看出当粗细晶的尺寸相差越大时，由于粗细晶应变梯度的增大，材料由于微观结构而得到的背应力强化和应变硬化率也越大，但是其中粗晶的尺寸越大，也会使得材料的屈服强度有所下降。

6、不同取向的异质双峰结构力学性能预测

图8为根据本发明的一个实施例的基于ABAQUS的异质结构模拟方法的不同取向下异质结构应力应变响应和应变硬化率对比图，其预测了不同取向粗晶的拉伸力学性能。可以从图中看出当粗晶的取向选择较硬取向(932)时，整体材料的屈服强度和应变硬化率都得到了较大的提高。

通过对比不同微观结构和取向的双峰结构预测模拟结果，可以更清楚的了解异质结构中的强化机制，并可以针对所需的力学性能调控相应的异质微观结构，指导材料设计。

需要说明的是，本发明的保护范围中现有技术部分并不局限于本申请文件所给出的实施例，所有不与本发明的方案相矛盾的现有技术，包括但不局限于在先专利文献、在先公开出版物，在先公开使用等等，都可纳入本发明的保护范围。

此外，本案中各技术特征的组合方式并不限本案权利要求中所记载的组合方式或是具体实施例所记载的组合方式，本案记载的所有技术特征可以以任何方式进行自由组合或结合，除非相互之间产生矛盾。

还需要注意的是，以上所列举的实施例仅为本发明的具体实施例。显然本发明不局限于以上实施例，随之做出的类似变化或变形是本领域技术人员能从本发明公开的内容直接得出或者很容易便联想到的，均应属于本发明的保护范围。

Claims (9)

1.一种基于ABAQUS的异质结构模拟方法，其特征在于，包括：

S1：建立均匀结构和异质结构的ABAQUS代表性单元模型；

S2：将滑移阻力模型和背应力模型分别修正为多尺度的滑移阻力模型和考虑强度差的背应力模型并写入ABAQUS的用户子程序UAMT，通过用户子程序UMAT定义晶体的弹塑性本构方程，该弹塑性本构方程中包括多个材料参数；

S3：根据两组不同晶粒尺寸的均匀结构和一组异质结构的单轴拉伸试验得到的拉伸曲线的实验结果与所述步骤S2中的弹塑性本构方程进行拟合，对所述步骤S2中的材料参数进行确定；

S4：施加周期性边界条件和外部载荷，利用ABAQUS计算每个晶粒的应力应变响应；

S5：对所有晶粒的应力应变平均化，得到异质结构的宏观力学性能；

S6：改变所述步骤S1中的ABAQUS代表性单元模型的微观结构和晶粒取向，预测和分析对应的宏观应力应变响应、强化效应和微观变形云图，通过对比获得最优的微观结构；

其中，所述步骤S2包括：

S21：采用投影运算方法，求解步骤S1中的异质结构的ABAQUS代表性单元模型的坐标系与晶粒坐标系之间的转化矩阵，基于所述转化矩阵，通过张量运算，得到异质结构的ABAQUS代表性单元模型的坐标系下的施密特因子，随后求解作用在晶粒中的应力张量σ^g，并结合材料弹性张量C^-1，进而计算出每个晶粒的弹性应变；

S22：由热驱动流动准则算出塑性应变，包括：

S221：在滑移阻力模型中引入Hall-Petch准则，将滑移阻力模型中的初始滑移阻力参数s₀修正为多尺度的初始滑移阻力参数S₀；

S222；在背应力模型中材料参数h_b中引入强度差模型，得到考虑强度差的背应力模型；

S223：根据基于热驱动流动准则的弹塑性变形理论求解晶粒分配给每个潜在滑移系的切变量

进而计算每个晶粒的塑性应变：

S23：结合所述步骤S21中的弹性应变和所述步骤S22中的塑性应变，得到每个晶粒的应力应变响应；

在所述步骤S221中，每个滑移系的滑移阻力为：

h^αβ＝h_s[w+(1-w)δ_αβ]，

其中，h_s为静态硬化模量，w代表晶格硬化率，δ_αβ代表克罗内克函数，s₀、S_sat为传统晶粒尺寸对应的初始滑移阻力和饱和滑移阻力，S₀为多尺度的初始滑移阻力参数，k_s为霍尔佩奇系数，d为晶粒尺寸，在粗细晶区域使用不同的霍尔佩奇参数s₀,k_s，以考虑异质结构中涉及的跨尺度设计准则，降低模拟误差，S^α为第α滑移系的滑移阻力；

在所述步骤S222中，每个滑移系的背应力为：

其中，h_b代表硬化率，

代表参考硬化率，k_b、r_D为材料参数，d_n为当前计算晶粒的周围晶粒平均尺寸，

为每个潜在滑移系的切变量，B^α为第α滑移系的背应力。

2.根据权利要求1所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，在所述步骤S223中，每个潜在滑移系的切变量

为：

其中，T为绝对温度，k为玻尔兹曼常数，B^α为第α滑移系的背应力，S^α为第α滑移系的滑移阻力，τ^α为第α滑移系的分解切应力，表示为τ^α＝σ^g:P^α，F₀为热激活位错越过障碍的激活能，τ₀、p、q为材料常数，

为参考剪应变率；

每个晶粒的塑性应变为：

其中，P^α为第α滑移系的方向因子张量，m^α为第α滑移系滑移方向的单位矢量，n^α为其滑移面法向的单位矢量。

3.根据权利要求1所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，在所述步骤S21中，所述异质结构的ABAQUS代表性单元模型的坐标系下的施密特因子α为：

α＝cosλcosβ，

其中，λ为滑移方向与外力F的夹角，β为滑移面与外力F中心轴的夹角；

所述弹性应变ε^e为：

ε^e＝C^-1:σ^g，

其中，ε^e为弹塑性本构方程中的弹性应变，σ^g为作用在晶粒中的应力张量，C^-1为材料弹性张量。

4.根据权利要求1所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，所述步骤S1包括：

S11：基于matlab语言根据Voronoi tessellation方法建立均匀结构的2D代表性体积单元，该2D代表性体积单元包括多个晶粒；

S12：根据步骤S11中的均匀结构的2D代表性体积单元的晶粒个数随机生成各个晶粒的取向参数，得到随机取向文件；

S13：通过ABAQUS的脚本程序读取所述步骤S11中的均匀结构的2D代表性体积单元和所述步骤S12中的随机取向文件，生成均匀结构的ABAQUS代表性单元模型，并作为异质结构的均匀细晶模型；

S14：根据需要研究的异质结构的微观结构的形貌图，在ABAQUS中对所述步骤S13中的均匀细晶模型进行手动修正，首先根据所述形貌图中的粗晶形貌将所述均匀细晶模型中的若干个晶粒合并成粗晶并删除内部多余线条，再删除每个合并的晶粒的取向参数并重新赋予粗晶取向，多次对比校正后，形成异质结构的ABAQUS代表性单元模型；

S15：使用ABAQUS的脚本程序计算所述步骤S14中的异质结构的ABAQUS代表性单元模型中每个晶粒的尺寸，并将计算结果赋予相应晶粒。

5.根据权利要求4所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，在所述步骤S12中，所述晶粒的取向参数采用三个欧拉角

φ、

表示。

6.根据权利要求1所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，在步骤S3中，所述材料参数是通过结合所述步骤S2中的弹塑性本构方程，采用试参法计算两组不同晶粒尺寸的均匀结构和一组异质结构的ABAQUS代表性单元模型的拉伸曲线的模拟结果，直至模拟结果与单轴拉伸试验所得的拉伸曲线的实验结果拟合，来确定的。

7.根据权利要求1所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，所述步骤S4包括：结合所述步骤S1中的异质结构的ABAQUS代表性单元模型、步骤S3中确定的材料参数，并施加的周期性边界条件和外部载荷，在ABAQUS中划分模型网格，再结合所述步骤S2中的弹塑性本构方程，通过ABAQUS软件计算得到每个积分点的应力应变张量。

8.根据权利要求1所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，在步骤S6中，所述最优的微观结构通过最大化的激发微观结构能够诱发的背应力强化，进而以牺牲最小甚至不牺牲塑性的情况下提高强度，来获得。

9.根据权利要求1所述的一种基于ABAQUS的异质结构模拟方法，其特征在于，所述异质结构包括异质双峰结构、层片结构和梯度结构。

Images