DOI：10.19476/j.ysxb.1004.0609.2019.05.10
基于位错理论的TC4钛合金的动态本构模型与数值模拟
李云飞1，曾祥国2
(1. 中国工程物理研究院 总体工程研究所，绵阳 621999；
2. 四川大学 建筑与环境学院，成都 610065)
摘 要：为描述Ti-6Al-4V(TC4)两相钛合金在高应变率、高温载荷条件下的复杂力学行为，基于细观塑性变形机理和位错动力学理论，从细观尺度构建一种两相钛合金粘塑性本构模型，并阐释各本构参数与微结构特征量之间的关联及其表征的物理意义。为确定本构参数并提高参数的识别效率与精度，提出一种基于拉丁超立方抽样、Spearman秩相关分析的参数敏感度整体分析方法，并在参数敏感度分析结果和基本遗传算法的基础上，建立基于改进小生境算法、可疑峰值点判断策略和局域精确搜索技术的改进遗传算法，得到了TC4两相钛合金的本构参数。采用应力补偿更新算法，通过显式用户子程序接口VUMAT将两相钛合金本构模型嵌入ABAQUS有限元软件中，实现了钛合金在高应变率、高温条件下动态本构行为的数值模拟。对比模拟结果与实验数据发现，所构建的本构模型描述材料高应变率条件下力学行为的准确性优于Johnson-Cook模型的。
关键词：Ti-6Al-4V；本构模型；位错动力学；VUMAT用户子程序
文章编号：1004-0609(2019)-05-0972-11 中图分类号：TG113.25 文献标志码：A
随着科学技术的迅猛发展，各工业化领域对材料尤其是金属材料的要求越来越高，需要了解其在各种极端条件下的性能才可充分、有效利用材料的潜能[1]。钛合金因优良的力学性能和耐蚀性能，在航空航天、造船、石油化工等领域得到广泛应用。钛合金主要分为α相、β相钛合金和α+β两相钛合金，目前两相钛合金的应用最为广泛，其中Ti-6Al-4V最具代表性。
由于对应变率、温度加载条件高度敏感，钛合金在高应变率下的动态力学性能与准静态下相比有很大差异[2-4]。早期诸多研究从材料的宏观动态力学行为出发，总结其动态行为的变化规律提出了经验型本构，如Johnson-Cook模型、Khan-Huang-Liang模型[5]等。在高应变率(1×103 s-1以上)条件下，上述经验型本构方程的流变应力预测误差较大，不能够理想地描述材料的塑性变形行为。高应变率下材料的力学响应行为与其微结构特征及其演化形式密切相关，因此，从微观力学的层面出发，基于塑性变形的物理基础有望建立更加准确地宏观本构模型。
目前，对钛合金本构模型的研究主要集中在α相、β相这两种晶体结构相对简单的钛合金，而针对具备相对复杂晶体结构的α+β两相钛合金本构研究还较少。NEMAT-NASSER等[6]通过归纳推广其FCC模型，提出了一种Ti-6A1-4V的本构模型。ZERILLI等[7]基于HCP晶体兼具体心立方(BCC)和面心立方(FCC)两种晶体结构的特征的论点，提出了一种针对HCP金属的本构模型(Z-A模型)。SEO等[8]通过引入一个温度相关函数对传统J-C模型进行改进，从而考虑了动态再结晶对流变应力的软化效应。
本文基于细观塑性变形机理和热激活位错动力学理论，综合考虑相关宏观因素(如实验温度、应变速率、应变量等)和部分微观因素(如晶粒大小、点阵类型、原始位错分布和种类、动态回复和再结晶、变形激活能、金属类型等)的共同影响，在前人工作基础上针对Ti-6Al-4V提出一种新的黏塑性动态本构模型，并阐释了本构参数与微结构特征量之间的关联及其表征的物理意义。所提出的两相钛合金本构模型共包含12个待定参数，为了确保本构参数识别结果的可靠性，减少参数识别的计算量，提出了一种基于拉丁超立方抽样、Spearman秩相关分析的参数敏感度整体分析方法，实现了对上述12个待定参数的敏感度分析，确定了各参数搜索区间合理的离散疏密程度。然后在参数敏感度分析结果和基本遗传算法的基础上，建立了一种基于改进小生境算法、可疑峰值点判断策略和局域精确搜索技术的改进遗传算法，对Ti-6Al-4V两相钛合金本构模型12个参数进行了优化识别。通过显式用户材料子程序接口VUMAT，采用应力补偿更新算法将该两相本构模型以及优化识别的本构参数导入ABAQUS有限元软件中，验证子程序的准确性，同时比较计算结果与实验数据，验证该本构模型在高应变率下的对两相钛合金本构行为预测的准确性。
1 基于位错理论的钛合金本构模型
金属钛具有同素异构现象，当温度在882 ℃以下为密排六方晶体结构(HCP)，简称α-Ti；在882 ℃和熔点1678 ℃之间，金属钛为体心立方形式晶格结构(BCC)，简称β-Ti[9]。利用上述两种结构的不同特点添加适当的合金元素，改变其相变温度和各相的含量形成不同形式的钛合金，主要分为α型、β型、α+β型，工业上相应称为TA、TB、TC型钛合金。
ZERILLI和ARMSTRONG[7]提出HCP晶体兼具体心立方(BCC)和面心立方(FCC)两种晶体结构的特征，其变形机制介于BCC和FCC晶体之间。MACDOUGALL等[10]研究发现，纯钛(BCC)的Zerilli- Armstrong本构模型可用于描述Ti-6Al-4V(HCP)的本构行为，这意味着HCP和BCC之间可能存在某些共同的结构特性；将HCP与FCC对比，HCP与FCC两种晶体的原子堆垛方式非常相像，因此可推断HCP的本构方程与BCC、FCC的本构方程相关，即可写作BCC与FCC本构方程的线性组合，则钛合金的本构方程为
(1)
式中：
、
分别为BCC与FCC的本构方程；
为塑性应变，
为实际应变率，T为试验温度；k1、k2为线性系数，分别代表FCC和BCC结构在HCP结构中所占的比例。
REGAZZONI等[11]认为金属的塑性变形在常规应变率下主要为位错热激活机制所控制。热激活机制是指运动位错必须克服不同势垒的阻碍。流变应力的大小实质上就是反映材料对位错运动的阻碍能力。以流变应力在0 K时的值
(力学阀值应力Mechanical Threshold Stress，MTS)[12]为基准，可将流变应力写成非热应力项
与热应力项
之和。根据BCC和FCC结构位错热激活机制的不同，其流变应力可做进一步探究。
对BCC而言，随着温度增加，热能会降低短程势垒高度，即热应力项
减小，而非热应力项
不变；对FCC而言，其热激活机制依赖于应变，其势垒主要是森林位错相互之间形成的阻碍作用。当T=0 K时尚没有热激活能，短程势垒高度最大；而随着温度的增加，热能会增加原子振动的振幅，降低势垒高度[13]，从而帮助位错克服势垒，热应力项
也相应减小，但非热应力项
不发生变化。因此，BCC和FCC的本构方程可分别写为式(2a)与(2b)：
(2a)
(2b)
式中：和分别为BCC结构与FCC结构的非热应力项；
、
分别为描述应变率和温度效应的影响因子(＜1)，称为热激活函数；与分别为BCC结构与FCC结构的热应力项。
由式(2a)、(2b)可看出，BCC和FCC的本构方程均包含3个部分：非热应力项、热激活函数、热应力项，下面对其分别讨论。
1) 非热应力项、
对BCC而言，其热激活面积与塑性应变无关，其势垒主要peierls内应力[14]，应变硬化项和反映应变率与温度效应的项之间是不耦合的，将其计入非热应力项；对FCC而言，需考虑材料颗粒尺寸对屈服应力的影响作用，其形式根据Hall-Petch关系得出。BCC 和FCC 本构方程中的非热应力项可分别写为式(3a)、(3b)：
(3a)
(3b)
式中：σG代表初始缺陷和杂质等带来的应力；d为颗粒直径；K·(εp)n为应变硬化项；K为应变硬化系数；n为应变硬化指数；C为常数。
2) 热激活函数、
金属内部的位错是由外力作用下产生的应力引起，OROWAN[15]最早提出了塑性应变率
与平均位错运动速度v之间的关系：
(4)
式中：
为Schmidt方向因子；b为Burgers矢量的模量，表示晶体因位错滑移引起的偏移；ρm为移动位错密度。
GILMAN等[16]根据实验与理论分析，提出位错运动速度
可进一步写为
(5)
式中：v0为参考位错速度；k为Boltzman常数；T为温度；而G0表示参考热激活能；r和s为决定势垒形状的一对相关联的常参数。
结合公式(4)~(5)并考虑热激活函数＜1的物理性质可知，BCC与FCC的热激活函数均可写为如式(6)所示的相同形式：
(6)
式中：
定义为参考应变率；p、q为决定势垒形状的一对相关联的常参数。
3) 热应力项、(均为阀值热应力)：
FOLLANSBEE等[12]发现，当应变硬化随着动态回复过程逐渐减小直至最后达到饱和状态时，饱和阀值应力
可用下述描述：
(7)
式中：
为饱和阀值应力在T=0 K时的参考值；a0为常数；μ为材料的剪切模量。
对BCC而言，其热应力与应变无关。且阀值热应力
与其在应变硬化达到饱和状态时的饱和值
相等，则有以下关系：
(8)
式中：
为饱和阀值热分量(在T=0 K时的参考值)。
对FCC而言，定义应变硬化率
，定义饱和阀值热应力为
(即应变硬化率为0时的阀值热应力)，FOLLANSBEE等[12]根据
—
和—的曲线关系图得出
(9)
分别将式(8)和(9)代入式(7)中得到BCC和FCC本构方程中的热应力项：
(10a)
(10b)
至此，BCC和FCC的本构方程中所包含的3个部分(即非热应力项、热激活函数、热应力项)的表达式已全部求出。联立式(3a)、(6)和(10a)并合并常系数，可得BCC的本构方程为
(11)
联立式(3b)、(6)和(10b)并合并常系数，可得FCC的本构方程为
(12)
将式(11)和(12)代入式(1)，并合并常系数，整理后可得两相钛合金的本构方程为
(13)
式中：
，K3=k2，K2=
，
，
，K1=k1K。需要通过相关实验数据确定的本构参数：K1、K2、K3、
、
、
、p、q、m、n、
、
，共12个待定本构参数。
2 本构参数优化选取
2.1 本构参数敏感度分析
Ti-6Al-4V两相钛合金本构模型共12个待定参数，每个参数对目标函数产生的影响程度不同。实际分析中若对所有参数都予以考虑，计算量将非常大，需对所有参数作筛选及优化处理，重点分析影响程度大的参数，以减少计算量。
传统的参数敏感度分析方法为单因素分析法，该方法忽略了本构参数间的相互作用对计算结果的影响。而文中的两相本构模型呈非线性，且各参数相互关联，所有参数会共同影响该模型对材料本构关系的拟合精度，因此需要进行参数敏感度的整体性分析，即当所有参数在指定范围内同时变化时，考察各参数对目标函数的影响程度。为综合分析所有变化的参数，首先需在整个参数空间中抽样。拉丁超立方抽样(LHS)[17]是一种适用于复杂多维参数空间抽样的方法，其基本实现过程如下，过程示意如图1所示。
1) 设有m个本构参数构成的参数集{x1, x2, …, xm}，将每个参数的取值范围分为等概率的n个互不重叠的区间。
2) 在每个参数xi的各个区间内随机抽取一个样本代表
，即每个参数xi含n个样本值。
3) 将每个参数xi的n个样本值均随机进行排列，共生成m个随机排列。
4) 每次从m个排列中顺序提取各参数的一个样本值，即形成一个样本参数集Xk；依次提取n次，即得n个样本参数集{X1, X2, …, Xn}。
完成参数总体抽样后，采用非参数统计方法的Spearman秩相关分析法进行随机参数之间的相关性分析，其基本实现过程如下。
1) 设有样本集{A1, A2, …, An}，把它们按大小顺序排列，构成一个有序序列A；若Ai在该序列A中占第Ri个位置，则称Ri为Ai在{A1, A2, …, An}中的“秩”。
2) 将n个样本参数集{X1, X2, …, Xn}代入式(1)可计算得目标函数解集{W1, W2, …, Wn}，即：每个参数xi的样本集
所对应的目标函数解集为{W1, W2, …, Wn}。
3) 设参数
在中的“秩”为
，
在{W1, W2, …, Wn} 中的“秩”为
，则参数
的随机输入样本集与其对应的输出结果集{W1, W2, …, Wn}的Spearman秩相关系数的绝对值即为参数xi的敏感度ri，其表达式为
(14)
将各本构参数xi的敏感度ri按大小进行排列，敏感度越大说明该参数对输出结果的影响越大，在后续的参数识别分析时需重点关注。利用该方法计算得到的参数敏感度分析结果为确定本构模型参数离散区间的疏密程度提供了理论基础，可极大地减少后续参数识别的计算工作量，确保参数识别结果的可靠性。
2.2 多参数非线性优化方法
本构参数优化识别的目的是寻找一组解，使模拟结果与观测结果之间的误差为最小，即寻找一组材料本构参数，使如下目标函数成立：
(15)
式中：W(X)为目标函数；X为离散设计变量；
为塑性应变
处所对应的流变应力实验测量值；为通过式(13)计算得到的塑性应变处所对应的流变应力模拟计算值；
和
分别表示参数xi的上限、下限值。此时目标函数W(X)取得极小值，相应模拟计算用的材料参数
即为两相钛合金本构模型参数的最优解。
图1 LHS方法的基本实现过程示意图
Fig. 1 Realization process schematic drawing of LHS method
基于参数敏感度分析结果，需采用遗传算法快速实现对本文建立的两相钛合金粘塑性本构模型12个待定参数的识别，但基本遗传算法(SGA)在应用中存在诸多缺陷，例如无法同时满足精度、可靠性和计算时间三方面的要求，且容易产生早熟现象，局部寻优能力较差等。本文针对基本遗传算法的诸多缺陷做了有效的改进，采用改进小生境算法、可疑峰值点判断策略和局域精确搜索技术，提高遗传算法的全局搜索能力，其关键实现要点如下。
1) 对小生境算法的改进
建立最优个体保存策略。各小生境独立进化，每代进化完毕后保存当前最优个体。保存策略为：在每代进化开始前首先复制当前最优个体作为复本，然后让所有个体平等地参与交叉、变异运算。该代进化完毕后，若该最优个体仍保留在群体中，则不做任何处理，否则就用当前最优个体的复本替换掉当代最差个体。在各小生境间建立数据交换机制。每进化一代都会在各小生境间进行最优个体的交换，即用第一个小生境的最优结果替换第二个小生境的最差结果，用第二个的最优个体替换第三个的最差个体，以此类推。该方法能在保证群体多样性的前提下提高优等个体的比例，加快收敛速度。
2) 建立可疑峰值点判断策略
每个小生境独立进化完毕后，均收敛于一个峰值点处。当小生境数目较大时，求得的所有峰值点必含全局最优点和局部最优点。在最终确定全局最优点前，对求得的所有峰值点均视为“可疑峰值点”。所有可疑峰值点中的最优点A必为全局最优点，记其目标函数值为F*，对应的变量为
。设某一可疑峰值点B的目标函数为F，其对应的变量为{x1, x2, …, xm}，则B点为相异于A点的全局最优点的判断条件为
(16)
式中：为0~1之间的常数，指最优点附近的目标函数值域相对搜索范围；βi为判断两个体是否“相邻”的常数；di为两个体的海明距离，若di≤βi，则说明A、B两个个体的相隔距离较近，经精确搜索后必收敛于同一个峰值点。式(16)的第一个不等式用于确定该点为全局最优点，第二个不等式用于确定B点为相异于A点的全局最优点。
3) 对所有相异的全局最优点作局域精确搜索
设某个全局最优点对应的变量集为，其各个变量的取值范围为xil≤xi≤xiu，则将各个变量的上下限边界条件改为
(17)
式中：
为变量xi附近的搜索邻域。
然后对各变量在缩小后的范围内，重新进行更精密的遗传进化处理，获得精度更高的该多峰值问题的全局最优峰值点。改进遗传算法的程序流程图如图2所示。
2.3 本构参数优化结果
采用文献中不同应变率、温度加载条件下的实验数据[5]，对Ti-6Al-4V两相钛合金动态本构模型待定参数的敏感度进行求解。得到不同取值范围、抽样次数下12个本构参数的敏感度计算结果如表1所示。将敏感度根据归一化敏感值大小分为3类，A代表非常敏感(敏感值范围0.65~1)；B代表一般敏感(敏感值范围0.2~0.649)；C代表不敏感(敏感值范围0~0.199)。
根据表1可知，各参数敏感度与其自身的取值范围、抽样次数有关。当抽样次数较小时，计算得的各参数敏感度的排序出现较大的振荡，杂乱无规律；随着抽样次数增大，各参数敏感度的排序逐渐趋于稳定。在抽样次数足够大的情况下，无论各参数设置的取值范围是否包含最优解，各参数的敏感程度都是较稳定的。因此，对抽样次数大的若干算例结果进行综合比较可得，各本构参数的敏感度排序不完全相同，但存在明显的规律：参数K2、的敏感度值相对较低，即K2、的变化在整个计算过程中对输出值的影响较小；参数K3、n、q、、m的敏感度值相对较大，即参数K3、n、q、、m对输出值的影响程度较大。
根据各参数的物理意义，确定了参数的初始搜索范围上下限分别为：K1
[650, 700]、K2 [1000, 1500]、K3[2000, 2500]、n[0.01, 0.99]、α[10-7, 10-5]、β[10-7, 10-5]、p[0.01, 0.99]、q[1.01, 1.99]、
[55, 65]、m[0.05, 0.06]、
[0.9×1011, 1.1× 1011]、
[2×109, 4× 109]。并基于各参数敏感度结果，增加敏感度较大参数的离散搜索精度，减小敏感度较小参数的离散搜索精度。采用上述改进遗传算法并基于参数敏感度分析结果与实验数据，优化识别了Ti-6Al-4V两相钛合金本构模型的12个本构参数值如表2所示。
图2 改进遗传算法的程序流程图
Fig. 2 Program flow diagram of advanced GA
3 两相动态本构的数值模拟
3.1 VUMAT与应力补偿更新算法
本文通过ABAQUS有限元软件二次开发，将建立的Ti-6Al-4V钛合金两相本构模型应用于数值模拟中。ABAQUS有两个求解器：隐式分析模块Standard和显式分析模块Explicit，都提供了FORTRAN语言编写的专门用于定义材料模型的用户子程序接口：UMAT和VUMAT。二者的数据结构是相似的，不同之处在于UMAT中需要形成Jacobian矩阵，即应力增量对应变增量的变化率，而在VUMAT中则不要求Jacobian矩阵，此外在子程序与主程序的接口处二者也有细微差别。本文采用显式用户材料子程序VUMAT接口对钛合金两相本构模型进行定义，适用于冲击、高速率应变等瞬时载荷下的结构响应有限元模拟，可有效缩短计算时间，提高计算效率。
表1 Ti-6Al-4V双相钛合金本构参数敏感度分析结果
Table 1 Sensitivities of constitutive parameters for Ti-6Al-4V
表2 Ti-6Al-4V钛合金两相本构参数优化值
Table 2 Optimized values of dual-phase constitutive based on experimental data
应力张量更新算法的选择是VUMAT子程序编写过程中的核心问题。传统粘塑性应力更新算法中若材料点已进入塑性变形阶段，则应变增量由弹性部分和黏塑性部分组成
，弹性部分应力增量按广义Hooke定律计算
，黏塑性部分应力增量为
(18)
式中：Dep为黏塑性张量，且Dep=De-Dp；De为弹性张量；Dp为塑性张量。塑性变形阶段的应力增量为Δσe与Δσep之和。该算法中每一增量步必须首先对Dep进行迭代更新，计算量大效率低，且易累积误差，影响计算结果精度。因此采用基于应力补偿的更新算法[18]，算法原理如图3所示。首先判断当前状态下材料点是否屈服，若屈服则应变增量
。假定应变张量全为弹性的，有广义Hooke定律计算试探应力增量
以及试探应力
。然后在试探应力的基础上减去多算的塑性应变增量对应的应力增量部分，对试探应力进行负补偿：
(19)
由式(19)可知，该算法没有迭代积分，计算中只需要恒定的弹性张量De，可以显著减小计算量，提高计算效率。
3.2 本构模型的子程序开发
ABAQUS用户子程序VUAMT的接口提供了依赖于结果的状态变量，用于存储在计算时用户定义的需要不断更新的状态变量。实现Ti-6Al-4V钛合金两相本构模型的子程序计算流程图和计算步骤如图4所示。
图3 应力补偿更新算法示意图
Fig. 3 Schematic diagram of stress compensation update algorithm
步骤(1) 从ABAQUS子程序接口读入本增量步的应变增量Δε，增量步开始时的应力张量σt和等效塑性应变
；
步骤(2) 计算试探应力
；
步骤(3) 调用子程序，计算初始流动应力σf；
步骤(4) 将试探应力代入屈服准则，判断是否屈服；
步骤(5) 若没有屈服，则
转到步骤(8)；
步骤(6) 若屈服，计算本增量步的塑性应变增量Δεp，采用径向返回应力补偿算法更新本增量步结束时的应力；
步骤(7) 更新内能、塑性功、等效塑性应变及各状态变量的值；
步骤(8) 结束，返回主程序。
4 结果与讨论
基于应力补偿更新方法，按照ABAQUS用户子程序接口规范，通过FORTRAN语言编写Ti-6Al-4V两相钛合金本构模型的VUMAT用户子程序，本构参数采用表1中的优化选取结果。考虑到材料行为与结构形态无关，本文采用单个八节点六面体等参元(C3D8R)计算Ti-6Al-4V钛合金在单轴加载条件下的力学响应。
图4 VUMAT子程序计算流程图
Fig. 4 Program flow diagram of VUMAT subroutine
图5 不同应变率、温度下Ti-6Al-4V模拟结果与实验值对比
Fig. 5 Comparison between simulation and experimental data of Ti-6Al-4V alloy under different strain rates and different temperatures
通过在ABAQUS/Explicit中调用VUMAT用户子程序，得到Ti-6Al-4V两相钛合金本构模型在高应变率(1900和2700 s-1)和不同温度(77、296和598 K)下的应力-应变曲线。本文中的流变应力预测(除应变率为0.001 s-1时)均为绝热过程，图中所标示温度为塑性变形过程中的初始温度。数值模拟结果与实验数据[19]对比如图5所示，二者吻合良好验证了用户子程序的准确性。由图5可知，在不同应变率且温度跨度较大的情况下，该本构模型可以较好地描述两相钛合金的本构行为。
图6所示为室温时在6000 s-1高应变率条件下，本文模型模拟结果和Johnson-Cook模型预测结果基于实验数据[19]的比较，其中Johnson-Cook模型本构参数优化识别结果如表3所示。表中：X为初始屈服应力；Y为应变硬化系数；Z为应变率强系数；n为应变硬化指数；m为温度软化系数。结果表明，高应变率条件下本文模型模拟结果与实验数据吻合相对较好，平均误差小于5%。而Johnson-Cook模型预测结果在应变大于15%时与实验数据的偏差逐渐增大，已不能很好地描述材料的本构行为。
图6 室温、6000 s-1应变率下Ti-6Al-4V不同模型模拟结果与实验数据对比
Fig. 6 Comparisons between simulation data and experimental data of Ti-6Al-4V alloy at 6000 s-1 and room temperatures
表3 Johnson-Cook模型本构参数优化值
Table 3 Optimized results of constitutive parameters for Johnson-Cook model
同时还对该本构模型准静态加载条件下的响应进行了数值模拟，并与Johnson-Cook模型预测结果以及实验数据[19]进行对比。由图7可见两相本构模拟结果与实验数据在小塑性变形段吻合良好，但应变大于15%时模拟结果小于实验数据并逐渐偏离。Johnson- Cook模型模拟结果比实验数据稍大，但趋势与实验数据一致并逐渐靠近，比两相本构的结果较接近实验数据。结合图5~7的数据比对结果可知，本文模型更适用于描述高应变率(1×103 s-1以上)加载条件下的本构行为。
图7 室温、准静态条件下的Ti-6Al-4V不同模型模拟结果与实验值对比
Fig. 7 Comparison between simulation data, experimental data and prediction of Johnson-Cook model under quasi-static condition
为了进一步验证用户子程序的准确性与可行性，将本文模型的模拟结果与其他实验数据(KHAN等[5]和HARDING等[10]的各一组实验数据)进行对比。如图8所示，本文模型模拟结果与KHAN等的实验数据吻合良好，表明其具备较好预测能力，而与HARDING等的实验数据则存在一定偏差。该现象主要因为HARDING等的拉伸与扭转实验中采用了Ti-6Al-4V薄壁管试件，而不是Hopkinson实验中的圆杆试件。
图8 Ti-6Al-4V本构模型模拟结果与不同实验结果对比
Fig. 8 Comparisons between simulation data and others’ experimental data of Ti-6Al-4V alloy
5 结论
1) 针对Ti-6Al-4V两相钛合金，基于细观塑性变形机理和位错动力学理论，从细观尺度提出了一种新的粘塑性本构模型。
2) 提出了一种基于拉丁超立方抽样、Spearman秩相关分析的参数敏感度整体分析方法，并在参数敏感度分析结果和基本遗传算法的基础上，建立了基于改进小生境算法、可疑峰值点判断策略和局域精确搜索技术的改进遗传算法，优化识别了Ti-6Al-4V两相钛合金的本构参数。
3) 利用径向返回应力补偿更新算法，通过VUMAT显式用户子程序将本构模型嵌入ABAQUS中，对钛合金不同应变率、不同温度下的动态响应进行数值模拟，模拟结果与实验数据吻合良好，通过对比发现所构建的本构模型在描述钛合金高应变率、温度跨度大的加载条件下的力学行为的准确性优于Johnson-Cook模型的。
REFERENCES
[1] 张长清, 谢兰生, 陈明和, 商国强. 高应变率下TC4-DT钛合金的动态力学性能及塑性本构关系[J]. 中国有色金属学报, 2015, 25(2): 323-329.
ZHANG Chang-qing, XIE Lan-sheng, CHEN Ming-he, SHANG Guo-qiang. Dynamic mechanical property and plastic constitutive relation of TC4-DT Ti alloy under high strain rate[J]. The Chinese Journal of Nonferrous Metals, 2015, 25(2): 323-329.
[2] MAJORELL A, STRIVATSA S, PICU R C. Mechanical behavior of Ti-6Al-4V at high and moderate temperatures— Part I: Experimental results[J]. Materials Science and Engineering A, 2002, 326(2): 297-305.
[3] 王冰洁, 郭鹏程, 李世康, 叶 拓, 曹淑芬, 李落星. 应变速率对AM80镁合金压缩变形行为的影响[J]. 中国有色金属学报, 2015, 25(3): 560-567.
WANG Bing-jie, GUO Peng-cheng, LI Shi-kang, YE Tuo, CAO Shu-fen, LI Luo-xing. Influence of strain rate on compression deformation behavior of AM80 magnesium alloy[J]. The Chinese Journal of Nonferrous Metals, 2015, 25(3): 560-567.
[4] CHIOU S T, TSAI H L, LEE W S. Impact mechanical response and micro structural evolution of Ti alloy under various temperatures[J]. Journal of Materials Processing Technology, 2009, 209(5): 2282-2294.
[5] KHAN A S, SUH Y S, KAZMI R. Quasi-static and dynamic loading responses and constitutive modeling of titanium alloys[J]. International Journal ofPlasticity, 2004, 20(12): 2233-2248.
[6] NEMAT-NASSER S, GUO W G, NESTERENKO V F, INFRAKANTI S S, GU Y B. Dynamic response of conventional and hot isostatically pressed Ti-6Al-4V alloys: Experiments and modeling[J]. Mechanics of Materials, 2001, 33(8): 425-439.
[7] ZERILLI F J, ARMSTRONG R W. High strain rate effects on polymer, metal and ceramic matrix composites and other advanced materials[C]//Proceeding of ASME International Mechanical Engineering Congress and Exposition. San Francisco: American Society of Mechanical Engineers, 1995: 121-126.
[8] SEO S, MIN O, YANG H. Constitutive equation for Ti-6Al-4V at high-temperatures measured using the SHPB technique[J]. International Journal of Impact Engineering, 2005, 31(6): 735- 754.
[9] ZHAN H, WANG G, KENT D, DARGUSCH M. Constitutive modeling of the flow behavior of 
titanium alloy at high strain rates and elevated temperatures using the Johnson–Cook and modified Zerilli–Armstrong models[J]. Materials Science and EngineeringA, 2014, 612(1): 71-79.
[10] MACDOUGALL D A S, HARDING J. A constitutive relation and failure criterion for Ti6Al4V alloy at impact rates of strain[J]. Journal of the Mechanics and Physics of Solids, 1999, 47(5): 1157-1185.
[11] REGAZZONI G, Kocks U F, Follansbee P S. Dislocation kinetics at high strain rates[J]. ActaMetallurgica, 1987, 35(2): 2865-2875.
[12] Follansbee P S, Kocks U F. A constitutive description of the deformation of copper based on the use of mechanical threshold stress as an internal state variable[J]. Acta Metallurgica, 1988, 36(1): 81-93.
[13] GAO C Y, Zhang L C, Yan H X. A new constitutive model for HCP metals[J]. Materials Science and Engineering A, 2011, 528(13-14): 4445-4452.
[14] Peierls P E. The size of a dislocation[J]. Proceedings of the PhysicalSociety, 1940, 52(1): 34-37.
[15] Orowan E. Problems of plastic gliding[J]. Proceedings of the PhysicalSociety, 1940, 52(1): 8-22.
[16] Johnson W G, Gilman J J. Dislocation velocities, dislocation densities, and plastic flow in lithium fluoride crystals[J]. Journal ofAppliedPhysics, 1959, 30(2): 129- 144.
[17] Simpson T W, Lin D K. Sampling strategies for computer experiments: design and analysis[J]. International Journal of Reliability and Applications, 2001, 2(3): 209-240.
[18] GAO C Y, Liu P H, Zhang L C. An explicit integration algorithm for introducing user-defined Thermo-Viscoplastic constitutive models in FE Simulations[C]//Proceeding of Materials Science Forum. Switzerland: Trans Tech Publications, 2011: 204-207.
[19] Follansbee P S, Gray G T. An analysis of the low temperature, low and high-rate deformation of Ti-6Al-4V[J]. Metallurgical Transaction A, 1989, 20(5): 863-874.
Dislocation-mechanics-based constitutive model of TC4 alloy and its numerical simulation
LI Yun-fei1, ZENG Xiang-guo2
(1. Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang 621999, China;
2. College of Architecture and Environment, Sichuan University, Chengdu 610065, China)
Abstract: A new physically-based constitutive model of dual-phase Ti-6Al-4V alloy (TC4) was developed based on the microscopic mechanism of plastic deformation and the theory of thermally activated dislocation motion. The relationships among constitutive parameters and microstructure characteristics, as well as the physical meanings were illustrated. In order to determine the constitutive parameters and improve the efficiency and precision of parameters identification, an overall analysis method of parameter sensitivity based on Latin Hypercube Sampling (LHS) method and spearman rank correlation method was presented. Furthermore, the constitutive parameters of the material were determined by a global genetic algorithm composed of an improved niche genetic algorithm, a global peak determination strategy and local accurate searching techniques. By adopting a stress compensation updating algorithm, a subroutine VUAMT of the proposed constitutive model of dual-phase titanium alloy was programmed on plat of ABAQUS/Explicit. Then finite element simulation of dynamic responses of TC4 alloy under loading conditions of high strain rate and high temperature was investigated. The comparison results indicate that the proposed model is more suitable and accurately to describe the mechanical behavior of material subjected to high strain rate than Johnson-Cook model.
Key words: Ti-6Al-4V; constitutive model; dislocation dynamics; VUMAT user subroutine
Foundation item: Project(U1430119) supported by the National Natural Science Foundation of China Academy of Engineering Physics and Jointly Set up “NSAF” Joint Fund
Received date: 2017-08-15; Accepted date: 2017-12-20
Corresponding author: LI Yun-fei; Tel: +86-816-2482384; E-mail: yunfei_ise@163.com
(编辑 龙怀中)
基金项目：国家自然科学基金委员会与中国工程物理研究院联合基金资助项目(U1430119)
收稿日期：2017-08-15；修订日期：2017-12-20
通信作者：李云飞，助理研究员，博士；电话：0816-2482384；E-mail: yunfei_ise@163.com
