高应变率下温度对单晶铁中孔洞成核与生长影响的分子动力学研究*

王云天 曾祥国 杨鑫

(四川大学建筑与环境学院, 深地科学与工程教育部重点实验室, 成都 610065)

采用嵌入原子势的分子动力学模拟方法, 研究了5 × 109 s—1应变率下, 温度效应对单晶铁中孔洞成核与生长的影响, 并对 NAG (nucleation and growth)模型在单晶铁中的适用性进行了探讨.结果表明：随着温度的升高, 单晶铁的抗拉强度峰值降低, 1100 K温度下单晶铁抗拉强度峰值比100 K温度下降低了35.9%.在100—700 K温度下, 拉应力时程曲线表现出双峰值特点, 分析表明, 第一峰值是由于拉应力升高引起内部结构发生相变而产生, 第二峰值则是因发生孔洞成核与生长而产生; 900—1100 K温度下, 拉应力时程曲线表现为单峰值, 孔洞成核与生长是拉应力下降的主要原因.分析发现, 孔洞在高温下更容易成核, 高应变率下单晶铁中孔洞成核与生长和NAG模型有较好的符合度, 单晶铁中孔洞成核阈值与生长阈值都远高于低碳钢,并且孔洞成核阈值与生长阈值随着温度的升高而逐渐降低.研究结果可为建立高应变率下金属材料动态损伤演化模型提供借鉴.

1 引 言

金属材料在高速冲击下的动态损伤问题是力学中最困难、最复杂的问题之一, 且在国防与民用领域中都有非常重要的应用背景, 因此该问题广受关注[1].在高速冲击下, 材料内部加载稀疏波和自由面反射稀疏波相互作用在内部产生拉应力区域,引起材料内部微孔洞成核、生长和贯穿, 最终由于损伤的累积导致材料断裂[2-6].

金属材料动态损伤受多种因素影响, 例如材料的微结构、初始缺陷、温度、加载应变率等, 其中加载应变率的影响非常重要, 金属材料的断裂强度随加载应变率的增加而增大, 在极高的应变率下材料的断裂强度与理论强度十分接近[7].在实验研究中,常用的研究金属材料动态损伤的方法主要有霍普金森杆、轻气炮、激光加载等, 加载应变率在102—1010s—1之间, 利用上述加载手段, 已经对不同应变率下金属材料的动态损伤进行了大量的研究[8-15].Zaretsky和Kanel[14]采用轻气炮加载方法, 研究了应变率达到 105s—1量级时多晶铁样品在不同初始温度下的动态响应, 分析了温度对层裂强度的影响以及冲击引起的相变行为.Remington等[2]利用激光加载方法, 研究了不同晶粒尺寸与加载应变率 (106—107s—1)下钽的动态响应, 结果表明晶粒尺寸对层裂强度有显著影响, 高应变率下孔洞在晶界处更容易成核.

在实验研究的基础上, 建立了描述金属材料动态损伤的理论模型, 代表性的理论模型包括：美国洛斯阿拉莫斯国家实验室提出的Gurson模型[16]和Johnson模型[17], 以及Curran等[3]提出的NAG(nucleation and growth)模型, 这些理论模型描述了宏观尺度下金属材料的动态损伤行为.金属材料在高速冲击下的动态损伤是一个非常复杂的问题,需要考虑的因素众多, 包括孔洞的成核与生长、相变、热耗散等, 并且时间尺度从皮秒到毫秒, 空间尺度从纳米到毫米[18], 因此, 难以在一个理论模型中将所有因素考虑在内.

目前, 分子动力学模拟已成为实验研究的重要补充手段, 为高应变率下金属材料动态损伤研究提供了新的方法.分子动力学模拟可以在微观尺度揭示金属材料的损伤过程, 并且可以直接测量相关物理量的变化, 例如温度、拉应力、相变和孔洞演化等数据[18-24], 有助于从根本上了解金属材料在高速冲击下的动态损伤机制.利用分子动力学模拟,Mayer[25]对单晶铁在应变率为 5 × 109s—1、单轴拉伸加载条件下的动态剪切与拉伸强度进行了研究,结果表明位错的均匀成核不会导致剪切应力的松弛, 同时发现拉应力时程曲线中有两个峰值出现,Mayer认为拉应力的第一个峰值是由于位错成核导致的塑性变形而产生的, 第二个峰值是由孔洞成核与生长而产生的.Rawat和Raole[19]分析了应变率为 1.5 × 108—1.5 × 1010s—1、三轴拉伸加载条件下单晶铁中的孔洞演化情况, 结果表明在高应变率下孔洞更容易成核, 但孔洞在较低的应变率下有更快的生长速度.Shao等[26]研究了应变率为2 ×109s—1、压缩加载条件下温度和预置空洞对单晶铁中结构变化的影响, 结果表明在相变过程中密排六方结构 (hexagonal close packed crystal structure,HCP)相的剪切应力总是低于体心立方结构(bodycentered cubic crystal structure, BCC)相的剪切应力.马文等[27]对冲击压缩条件下纳米多晶铁的相变过程进行了研究, 结果表明晶粒尺寸对晶界变形及相变原子比例有显著影响.分子动力学模拟的结果也可以为理论模型提供参数, Sugandhi等[28]通过粒子群优化算法从分子动力学模拟的结果中获得NAG模型的最佳拟合参数.

Rawat等[29]采用分子动力学方法研究了不同温度下单晶铜在 5 × 109s—1应变率下的动态损伤情况, 结果表明温度对单晶铜中孔洞成核与生长有显著影响, 同时对NAG模型在微观尺度的适用性进行了分析.目前已经进行的关于铁的动态损伤研究中, 已有关于不同应变率、相变、晶粒尺寸等因素的研究, 但在高应变率下温度的影响尚少有涉及, 并且NAG模型在微观上能否描述单晶铁的动态损伤过程的研究也比较少.因此, 本文使用分子动力学模拟方法, 建立了单晶铁的三轴拉伸模型,该模型广泛应用于研究金属材料在高应变率下的动态损伤中[30-32].本文分析了温度对单晶铁在高应变率下孔洞成核与生长的影响, 并探讨了NAG模型在单晶铁中的适用性.研究结果可为建立高应变率下金属材料动态损伤演化模型提供参考.

2 分子动力学模型

本文采用Mendelev等[33]提出的嵌入原子法(embedded atom method, EAM)描述单晶铁原子间相互作用, 该势函数广泛应用于高应变率下单晶铁的分子动力学模拟研究中[34,35].分子动力学模拟采用Lammps软件[36]进行, 首先建立单晶铁的分子动力学模型, 模型尺寸为 28.55 nm × 28.55 nm ×28.55 nm, 模型x,y和z坐标分别沿 [100], [010],[001]方向, 模型共包含 2 × 106个原子, 如图1 所示.在3个方向施加周期性边界条件, 在加载之前采用NPT系综对模拟系统进行充分弛豫(30 ps),以确保模拟系统在加载前达到指定温度的平衡状态, 其控制温度分别为 100, 300, 500, 700, 900,1100 K.弛豫结束后, 使用速度标定法控制温度,采用NVE系综并在x,y和z坐标轴方向施加相等的应变率使系统处于三轴拉伸状态.

图1 单晶铁三轴拉伸模型Fig.1.Model of single crystal iron under triaxial tension(atoms are colored by CNA).

使用Ovito软件[37]对分子动力学模拟结果进行可视化分析, 采用Python语言编写的程序进行孔洞体积分析, 其基本思想是将整个模拟区域划分为小的立体方格并计算空方格的个数, 然后统计在任意时刻的空方格数量来计算孔洞体积, 这种孔洞体积统计方法已被用于计算BCC晶体与面心立方结构 (face-centered cubic crystal structure, FCC)晶体的孔洞体积分析中[29,32,38]; 模型内部结构变化使用位错分析 (dislocation analysis, DXA)[39]、径向分布函数 (radial distribution function, RDF)分析与共邻分析 (common neighbor analysis, CNA)[40]等方法来进行; 使用基于C语言编写的改进遗传算法程序确定NAG模型最佳拟合参数.

3 计算结果与讨论

本节讨论了在应变率为 5 × 109s—1、三轴拉伸加载下单晶铁在100—1100 K温度下的动态响应.分析了温度对拉应力、孔洞体积分数、系统势能的影响, 对拉应力在不同温度下的变化特点进行了探讨, 同时还讨论了NAG模型在高应变率下单晶铁孔洞成核与生长的适用性问题, 确定了NAG模型的最佳拟合参数, 并分析了温度对NAG模型参数的影响.

3.1 温度对拉应力的影响

通过分子动力学模拟计算了单晶铁在100—1100 K温度下拉应力随时间变化情况, 拉应力时程曲线见图2.分子动力学模拟得到单晶铁在300 K 温度下最大拉应力约为 18.3 GPa, Mayer[25]研究单晶铁在应变率为 5 × 109s—1、单轴拉伸下得到的最大拉应力约为 17.6 GPa, Ashitkov 等[13]在铁(样品厚度为50 mm)的层裂实验中测得的最大拉应力为 17.8—20.3 GPa.从图2 可以看出, 在100—700 K温度下, 拉应力时程曲线有两个峰值出现, 而在900—1100 K温度下拉应力时程曲线只有单峰值出现.

100—1100 K温度下拉应力峰值如图3所示.从图3可以看出, 随着温度的升高, 拉应力峰值逐渐降低, 最大拉应力为 100 K 时的 21.4 GPa, 最小拉应力为 1100 K 时的 13.7 GPa, 降低幅度为35.9%, 表明温度对单晶铁的拉应力峰值有明显的影响, 相对而言, 应变率的变化对单晶铁的拉应力峰值的影响并不明显[19].

图2 100—1100 K 温度下拉应力随时间的变化Fig.2.Tensile stress as a function of time at 100—1100 K.

图3 拉应力峰值随温度的变化Fig.3.Relationship of peak pressure and temperature.

双拉应力峰值时程曲线以100 K温度为例,从图2可以看出, 拉应力随时间持续增长, 在15.7 ps时出现第一个峰值 (约为 21.4 GPa)后开始下降, 在 16.5 ps 时, 拉应力停止降低, 随后有所增长 (约增长 0.03 GPa), 19.1 ps时达到第二个峰值 (约为 16.9 GPa), 19.1 ps之后, 由于孔洞成核与生长, 拉应力快速下降.在 5 × 109s—1、单轴拉伸下对单晶铁的研究中, Mayer[25]也观察到拉应力时程曲线的双峰值现象, 认为第一个拉应力峰值内结构发生塑性变形, 第二个拉应力峰值内孔洞成核.双拉应力峰值时程曲线不仅在BCC结构的单晶铁中出现, 在FCC结构的单晶铜中也观察到类似现象, Rawat等[29]在单晶铜的三轴拉伸分子动力学模拟中观察到单晶铜在1250 K温度下(接近铜的熔点)拉应力时程曲线有两个峰值出现,Rawat认为拉应力的第一个峰值是由于发生了相变, 第二个峰值是由于孔洞成核与生长, 这说明Mayer与Rawat等对拉应力时程曲线的双峰研究结果一致.不同的是, FCC结构的单晶铜在1250 K温度下出现拉应力双峰值现象, 而BCC结构的单晶铁是在100—700 K温度下出现, 这表明在三轴拉伸下温度对BCC和FCC晶体结构的微观损伤演化的影响是不同的.

为探究拉应力时程曲线出现双峰值的原因, 对模拟结果进行了如下分析：1)分析拉应力时程曲线与孔洞体积分数关系, 即确认孔洞的成核阈值与成核时间; 2)通过DXA, 辨别塑性变形发生时间;3)借助晶体结构分析, 利用RDF分析和CNA方法来分析结构变化.为便于对比分析, 在图2中以100 K温度为代表标示出了拉应力双峰值之间的四个特征点, 分别为：第一个峰值点(A点)、第一个峰值后拉应力快速下降点(B点), 第一个峰值点后拉应力再次增长点(C点)、第二个峰值点(D点), 100—700 K温度下四个特征点时间如表1所列.

表1 100—700 K 温度下四个特征点时间Table 1.Four characteristic points time at 100 —700 K.

从表1 可以看出, 在 100—500 K 温度下, 第一个拉应力峰值出现的时间逐渐提前, 500 K温度下第一个拉应力峰值出现的时间比100 K温度下提前了 2.9 ps, 700 K 温度下第一个拉应力峰值出现时间较500 K温度下有所延后, 但相差不大, 仅为 0.5 ps; 100—700 K 温度下拉应力再次增长的时间以及第二个拉应力峰值出现的时间随温度的升高逐渐提前, 700 K温度下拉应力第二个峰值较100 K温度下提前了1.5 ps.

3.2 温度对孔洞体积分数的影响

图4给出了100—1100 K温度下孔洞体积分数随时间的变化.从图4可以看出, 尽管100—700 K与900—1100 K温度下拉应力时程曲线有不同的特点, 但是孔洞体积分数随时间的变化趋势基本相同, 随着温度的升高, 孔洞出现的时间提前,孔洞体积分数增长速率有所加快, 100—1100 K温度下孔洞体积分数在 30 ps时基本相同, 约为29%.100—1100 K温度下孔洞成核时间如图5所示, 从图5 可以看出, 1100 K 温度下孔洞成核时间是 12.8 ps, 100 K 温度下孔洞成核时间为 17.4 ps,1100 K温度下孔洞成核时间比100 K温度下提前了26.4%, 这表明温度对孔洞成核时间有显著影响, 孔洞成核时间随温度升高而逐渐提前.图6给出了体积分数为0.1时, 100—1100 K温度下系统内部孔洞分布情况.从图6可以看出, 系统内部的孔洞总体积是多个孔洞的成核、生长和贯穿累积而来, 1100 K温度下孔洞数量明显高于100 K温度,这表明随着温度的升高系统内部孔洞成核数量更多, 导致孔洞体积分数增长速率随温度升高而加快.

图4 100—1100 K 温度下孔洞体积分数随时间的变化Fig.4.Void volume fraction as a function of time at 100—1100 K.

图5 100—1100 K 温度下孔洞成核时间Fig.5.Void nucleation time at 100—1100 K.

图6 100—1100 K 温度下内部孔洞分布 (孔洞体积分数为 0.1 时)Fig.6.Distribution of voids at 100—1100 K (void volume fraction is 0.1).

图7给出了孔洞体积分数随时间的变化情况,并与拉应力时程曲线进行了对比, 图中虚线对应孔洞成核的时间(以tn表示), 孔洞成核后体积分数开始增长.从图7可以看出, 拉应力时程曲线在特征点C之前孔洞体积分数一直保持为0; 在特征点C附近, 孔洞成核, 孔洞体积分数开始增长, 但增长速度较慢, 直到特征点D之后, 增长速率才明显加快.图8给出了100—700 K温度下拉应力时程曲线特征点D现, 即孔洞没有成核; 在第一个峰值之后, 拉应力下降到特征点C附近时的孔洞分布情况, 可以看出系统内部有明显的孔洞产生.因此, 100—700 K温度下拉应力达到第一个峰值时没有孔洞出时, 孔洞开始成核, 之后孔洞体积分数持续增长.

100—1100 K温度下孔洞体积分数的发展可以分为3个阶段, 以900和1100 K为例进行说明,如图7所示.从图7可以看出, 孔洞体积分数增长可以分为3个阶段.

1)初始增长阶段 (S1), 与 100—700 K 温度下相同, 在拉应力达到峰值之前已有孔洞成核, 此时孔洞体积增长主要是孔洞成核带来, 初始增长阶段孔洞体积分数增长缓慢, 表明孔洞还未开始大规模生长, 直到拉应力时程曲线峰值拐点之前.从图5可以看出, 随着温度的升高S1阶段逐渐提前,1100 K比100 K温度下孔洞出现时间提前了约4.6 ps, 这表明孔洞在高温下更容易成核.

2)快速增长阶段(S2), 当孔洞体积分数达到0.01后, 成核孔洞开始快速生长, 是孔洞体积分数迅速升高的主要原因.S2阶段孔洞体积分数增长速率最快, 孔洞体积分数变化趋势与指数函数特点相似.通过与拉应力时程曲线的对比, 可以发现孔洞体积分数快速增长与拉应力时程曲线的快速下降对应, 表明孔洞体积分数的快速增长导致了拉应力的快速下降.1100 K温度下S2阶段持续时间约为 3 ps, 100 K 温度下 S2 阶段约为 7.5 ps, 二者相差2.5倍, 这表明高温下孔洞生长的速度更快, 这也导致了高温下拉应力下降速率更快.

3)线性增长阶段(S3), 此时孔洞生长减慢, 孔洞贯穿带来体积增长, 增长速率适中, 在 30 ps时不同温度下孔洞体积分数数值相差不大, 约为29%.这表明单晶铁在 5 × 109s—1应变率三轴拉伸下, 温度对最终孔洞体积分数的影响并不明显, 这与应变率对孔洞体积分数的影响不同[19].

3.3 位错分析

图7 100-1100 K 温度下孔洞体积分数与拉应力的关系 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.7.Void volume fraction as a function of time at 100-700 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 K.

孔洞体积分数分析表明拉应力时程曲线的特征点C附近孔洞开始成核与生长, 拉应力时程曲线第一峰值内没有孔洞出现.使用Ovito软件对分子动力学模拟结果进行了DXA, 图9给出了100—1100 K温度下位错密度随时间的变化情况, 图中时间轴上的标记点(红色圆点)是拉应力时程曲线峰值出现的时间, 在100—700 K温度下是第二个峰值出现的时间.

DXA 表明, 单晶铁在应变率为 5 × 109s—1、三轴拉伸加载下有三种类型位错出现：1/2位错、位错、位错, 这些位错类型在单晶铁及其合金的实验中都有观测到[41].除了这三种类型的位错之外, 还有一部分无法识别的位错类型,定义为其他(Others).从图9可以看出, 在100—1100 K温度下, 单晶铁中 1/2位错是最主要的位错类型, 占比最大,位错次之, 占比略高于其余两种位错类型.100 K温度下位错密度在21 ps开始增长, 300—900 K 温度下位错密度在20 ps开始增长, 1100 K 温度下位错密度在18 ps开始增长.1100 K 温度下位错出现时间比100 K 温度提前了 3 ps左右, 这表明温度的升高使位错出现的时间有所提前.通过与拉应力时程曲线的对比分析, 可以发现位错密度增长时间比拉应力时程曲线峰值(100—700 K拉应力第二个峰值)出现时间滞后 2—3 ps左右.对于 100—700 K温度下拉应力时程曲线, 位错密度开始增长的时间位于拉应力时程曲线第二个峰值拐点之后, 此时拉应力开始快速下降, 这表明拉应力时程曲线中第一个峰值内没有出现因位错而引起的塑性变形, 这与Mayer[25]的结论不同.Rawat等[29]在单晶铜的三轴拉伸分子动力学模拟中认为拉应力时程曲线第一个峰值内发生了相变, 为了探讨拉应力时程曲线第一个峰值内是否发生了相变, 对分子动力学模拟结果进行了RDF分析和CNA.

图8 100-700 K 温度下拉应力时程曲线第二峰值点孔洞分布Fig.8.Void distribution on the second-peak of tensile stress at 100-700 K.

3.4 径向分布函数分析

径向分布函数通过给定某个粒子的坐标, 计算其他粒子在空间的分布概率随距离变化情况, 常用于研究物质内部的有序性.晶体具有规则的周期性结构, 原子在其晶格位置附近波动, 由于晶体有序的结构, 径向分布函数有长程的峰, 在晶体结构转变或被破坏后, 长程的峰消失, 因此可以用径向分布函数来反应模拟体系内部晶体结构变化.使用Ovito软件计算了100—1100 K温度下整个系统RDF 变化情况, 100—700 K 温度下选取了 0 ps、拉应力时程曲线的四个特征点对应时间, 900—1100 K温度下选取了0 ps、孔洞成核时间、拉应力峰值时间以及拉应力峰值后1 ps, 如图10所示.可以看出, 在 100—700 K 温度下 0 ps时系统内部长程有序, 远距离处依然有较明显的曲线波动, 在拉应力时程曲线第一个峰值时远距离处的曲线波动开始减弱; 特征点B过后, 远距离的曲线波动基本消失, 只有近距离内出现明显的曲线波动, 系统内部长程有序消失.在 900—1100 K 温度下, 0 ps时系统内部长程有序, 孔洞成核使系统内部长程有序减弱, 在拉应力达到峰值后, 内部长程有序完全消失.

RDF分析表明, 100—700 K温度下拉应力达到第一个峰值后系统内部长程有序开始消失, 表明在拉应力达到第一个峰值后系统内部发生了结构改变; 在 900—1100 K 温度下, 孔洞成核引起系统内部长程有序消失.

3.5 共邻分析

RDF分析表明, 拉应力时程曲线的第一个峰值内系统内部可能发生了相变, 使用Ovito软件中的自适应共邻分析方法, 对系统的晶体结构变化进行了分析.图11给出了100—1100 K温度下系统内部晶体结构变化情况, 图中FCC代表面心立方结构、BCC代表体心立方结构、HCP代表密排六方结构、Other代表未知结构.100—700 K 温度下选取拉应力时程曲线的四个特征点对应的时间,900—1100 K温度下选取孔洞成核时间与拉应力峰值时间, 将模型中间进行切片观察, 切片厚度为5 Å, 图中使用黑色圆圈将内部孔洞标示出来.

图9 100-1100 K 温度下位错密度变化情况 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.9.Evolution of dislocation density as a function of time at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K;(f) 1100 K.

从图11可以看出：在 100—700 K温度下拉应力达到第一个峰值时, 系统内部出现FCC,HCP和未知结构, 随着拉应力的降低, FCC,HCP和未知结构数量逐渐增多; 在拉应力达到特征点C时, 大部分BCC结构已经转变为FCC,HCP和未知结构, 只有少量BCC结构留存; 在拉应力达到第二个峰值时, 可以观察到多个孔洞出现, 孔洞被未知结构包围.RDF 分析表明, 拉应力时程曲线的第一个峰值内系统内部发生了相变,BCC结构转变为FCC, HCP和未知结构, 没有孔洞出现; 在第二个峰值时, 只有少量BCC结构留存, 孔洞在未知结构区域内成核.在 900—1100 K温度下系统内部未知结构数量增多, FCC结构数量有所减少, 特别是在 1100 K 温度下, 系统内部出现大量的未知结构, 在拉应力达到峰值时可以看到明显的孔洞出现, 孔洞被未知结构包围.在铁纳米管和纳米晶体铁的单轴拉伸分子动力学模拟中也出现了BCC结构到HCP结构的相变出现[42,43].本文中出现的BCC结构到HCP结构的转变可能是α-ε相变, 在铁的冲击加载实验和分子动力学模拟中都有观察到[44-47].

图10 100-1100 K 温度下径向分布函数变化情况 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.10.Radial distribution function of the system at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K;(f) 1100 K.

图11定性地给出了系统内部结构变化情况,为定量地描述系统的结构变化, 计算了100—1100 K温度下FCC, BCC, HCP和未知结构随时间的变化情况, 如图12 所示.从图12 可以看出, BCC 结构首先转变为未知结构, 在拉应力时程曲线达到第一个峰值后, FCC, HCP 结构开始出现; FCC,HCP结构在拉应力时程曲线达到第二个峰值后占比达到最大.在 100—300 K 温度下, 二者占比约为30%, 随着温度的升高, FCC结构占比下降,500 K 时占比约为 25%, 700 K 时占比约为 20%;BCC结构占比最低点出现在FCC, HCP占比峰值之后, 100 K 时约为 10%, 并且随着温度的升高BCC结构占比逐渐降低, 900 K温度下约为2%.从图7可以看出孔洞体积分数在第二个拉应力峰值后开始快速增长, 对比图12结构占比变化, 可以发现孔洞体积的快速增长导致FCC, HCP结构占比降低, BCC结构占比再次升高, 引起这种变化的机理将在以后的文章中进行进一步的分析.

对100—700 K温度下拉应力时程曲线的分析表明系统内部首先发生了相变, 随后孔洞成核与生长, 但是随着温度的升高, 在 900—1100 K 温度下拉应力时程曲线只有一个峰值出现.从图12可以看出, 在孔洞成核之前系统内部已经发生相变,FCC, HCP结构占比已经开始增长, 但与 100—700 K 温度下不同的是, BCC 到 FCC, HCP 的相变并未导致拉应力的降低, 原因可能是FCC,HCP相变发生的时间与孔洞成核的时间非常接近,孔洞成核时FCC, HCP相变的数量非常少, 因此由FCC, HCP相变引起的拉应力降低幅度很小,拉应力只出现一个峰值.结合图7孔洞体积分数变化情况, 在拉应力达到峰值后孔洞体积分数开始快速增长, 这表明孔洞成核与生长是导致900—1100 K温度下拉应力降低的主要原因[30,48].因此, 在100—700 K温度下, 拉应力下降的原因首先是系统内部发生相变, BCC结构向FCC, HCP和未知结构转变, 然后是孔洞成核与生长; 在 900—1100 K 温度下, 孔洞成核与生长是导致拉应力下降的主要原因.

图11 100-1100 K 温度下内部结构变化 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.11.Snapshots for the structural changes at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 K.

为了对比初始BCC结构与重生成的BCC结构之间的区别, 选取模型中部进行切片观察, 为了能够更加明显地辨识孪晶结构, 选择切片厚度为2 Å.图13 给出了 300 K 温度下系统内部结构变化情况, 红色箭头标示出孪晶出现的区域.从图13可以看出, FCC和HCP结构在拉应力第二个峰值之后转变为BCC结构, 孔洞被未知结构包围, 当孔洞周围的未知结构原子随着孔洞生长而移动时, 导致原本属于FCC和HCP结构的原子重新排列, 再次产生 BCC 结构.对比 24, 30 和 0 ps时的BCC结构, 可以发现重新生成的BCC结构与初始BCC结构并不相同, 重新出现的BCC结构可能是变形孪晶的结果[49], 与BCC结构Mo和BCC结构Ta的孔洞生长分子动力学模拟中出现的变形孪晶情形相似[50].

图12 100-1100 K 温度下内部晶体结构占比 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.12.Crystal structure fraction as a function of time at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K;(f) 1100 K.

通过孔洞体积分数, DXA, RDF, CNA 等分析方法对分子动力学模拟结果进行了分析, 结果表明拉应力时程曲线的第一个峰值内系统内部发生了相变, 单晶铁原有的BCC结构转变为FCC,HCP和未知结构, 没有出现位错引起的塑性变形;在拉应力时程曲线的第二个峰值内孔洞开始成核与生长, 孔洞体积的快速增长使FCC和HCP结构占比降低, 孔洞的生长导致BCC结构重新产生,但与初始的BCC结构不同, 新出现的BCC结构可能是变形孪晶的结果.

3.6 温度对原子键能及系统势能的影响

图14是铁原子键能在不同体系温度下的变化情况, 图中ε是摩尔键能, 单位 kJ/mol, 摩尔键能的计算参照 Eberhart和Horner[51]的文章提出的方法, 计算时间点为0 ps, 即模拟体系充分弛豫后的时刻, 计算所需数据来自美国国家标准与技术研究院(NIST)[52].从图14可以看出, 铁原子键能在100 K 时为 104.68 kJ/mol, 随着模拟体系温度的升高, 铁原子键能逐渐降低, 在 1100 K 时铁原子键能为 93.99 kJ/mol, 降低幅度为 10.21%.这表明温度的升高显著影响了铁原子键能, 在高温下铁原子键能明显降低.

图13 300 K 温度下系统内部结构变化Fig.13.Structural changes at different time at 300 K.

图14 300-1100 K 温度下铁原子键能Fig.14.Bond energy of iron at 300-1100 K.

图15 100-1100 K 温度下系统势能Fig.15.Potential energy of the system at 100-1100 K.

图15是不同温度下系统势能随时间变化情况.从图15可以看出, 系统总势能随着温度的升高而增加.在 100—700 K 温度下, 拉应力时程曲线第一个峰值(图中红色圆圈标示)导致系统势能减小, 但幅度不大, 并且随着温度升高这种影响越来越低, 此时系统内部并未有孔洞成核, 是相变引起的系统势能变化; 系统势能峰值出现在拉应力时程曲线第二个峰值拐点之后, 此时孔洞体积分数快速增长, 表明孔洞成核与生长引起势能降低.在900—1100 K温度下, 系统势能随时间持续增长,势能曲线的峰值出现在拉应力时程曲线拐点之后,随后快速降低, 表明在高温下孔洞成核与生长是导致势能降低的主要原因.当有孔洞成核时, 原子间的相互作用会变得更强, 同时键长缩短, 导致势能的快速下降.

3.7 NAG模型

Curran等[3]和Seaman等[53]通过对铝和铜等金属材料的研究, 揭示出其断裂是大量准球形空洞成核和生长的结果.通过研究发现单位体积的孔洞数量以及孔洞体积与拉应力呈指数关系, 从而提出一种考虑微孔洞成核和生长效应的损伤模型—NAG模型.NAG模型是一种微细观物理模型, 主要用于描述金属材料中孔洞的成核和生长而发生的损伤过程.通过3.2部分对孔洞体积的分析, 发现单晶铁在 5 × 109s—1应变率 100—1100 K 温度下孔洞体积非线性增长阶段与NAG模型相符合.因此计算了100—1100 K温度下NAG模型参数并且分析了温度对NAG模型参数的影响.

3.7.1 NAG 模型介绍

孔洞的总体积由两部分组成：孔洞成核体积和孔洞生长体积.

1)孔洞成核体积

当材料中的拉应力Ps超过材料成核阈值Pn0时, 新的孔洞生成, 成核率按下式计算：

在Δt时间间隔内孔洞成核体积由下式计算：

式中Rn为形核尺寸参数.在NAG模型的计算中,形核尺寸参数Rn的值为3.1 Å, 即1.01倍铁的晶格常数, 也是孔洞体积分析中立方体方格的尺寸.

2)孔洞生长体积

当拉应力Ps超过孔洞生长阈值Pg0时, 成核的孔洞就会生长, 这时新的孔洞体积由下式计算：

式中,Vv0为在每个时间间隔初始时的孔洞体积,η为材料黏度, 在Δt时间间隔内孔洞成核与生长的总体积为

每个数据点的相对误差φ定义如下：

式中,VMD为分子动力学 (molecular dynaimics,MD)模拟中得到的孔洞体积分数,VNAG为NAG模型得到的孔洞体积分数.

相对平方误差Σ定义如下：

式中M为数据点总数.

3.7.2 NAG 参数分析

对于任意给定的一组NAG参数, 在已知拉应力时程曲线的前提下, 可以得到孔洞总体积随时间变化的函数关系.参数拟合采用基于改进遗传算法编写的C语言程序进行, 得到了100—1100 K温度下 NAG 模型拟合参数使用拟合得到的NAG模型参数计算了孔洞的体积分数变化情况, 并与MD模拟的结果进行了对比,如图16所示.从图16可以看出, 不同温度下单晶铁的孔洞体积分数变化与NAG模型符合较好.

图16 100-1100 K 温度下 NAG 与 MD 的孔洞体积分数结果的对比 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.16.Comparison of void volume fraction between the NAG model and MD at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K;(d) 700 K; (e) 900 K; (f) 1100 K.

表2 100-1100 K 温度下 NAG 模型最佳拟合参数Table 2.Best-fit NAG parameters at 100-1100 K.

除了与MD结果进行对比之外, 还将100—1100 K 温度下 NAG 模型拟合参数 (Pn0,P1,,Pg0,η)与实验中得到的低碳钢 (mild steel) NAG模型参数[54]进行了对比, 详见表2.

从表2可以看出：

1)单晶铁的孔洞成核阈值(Pn0)远高于低碳钢中的成核阈值, 差距在10倍以上.这是由于低碳钢是多晶结构, 存在大量的缺陷, 孔洞在晶界处更容易成核[2], 而在单晶铁中不存在这些缺陷, 所以孔洞成核阈值更高.

2)单晶铁的孔洞成核拉应力灵敏度系数(P1)较低, 约为低碳钢的20%, 这表明单晶铁中的成核速率对拉应力的敏感性比对多晶材料的敏感性高得多.

3)与低碳钢相比, 单晶铁的材料黏度η非常低, 相差约1000倍;

4)随着温度的升高, 孔洞成核阈值Pn0和孔洞生长阈值Pg0都逐渐减小.高温下, 原子运动更加剧烈, 加速了原子间相互分离的趋势, 导致孔洞成核与生长阈值降低;

5)成核阈值与生长阈值的比值在单晶铁与低碳钢中相差不大, 单晶铁的成核阈值Pn0约为生长阈值Pg0的6倍, 低碳钢中成核阈值是生长阈值的5倍[54].

4 结 论

本文使用分子动力学模拟方法, 研究了100—1100 K 温度下单晶铁在 5 × 109s—1应变率三轴拉伸加载的动态响应, 通过拉应力变化、孔洞体积分数以及微结构分析研究了单晶铁的动态响应以及温度效应的影响, 并且对NAG模型在高应变率下单晶铁中的适用性进行了探讨, 得到以下主要结论.

1)在 100—700 K温度下, 拉应力时程曲线出现两个峰值, 采用 DXA, RDF, CNA 和孔洞体积分数等方法对结果进行了分析.分析表明, 100—700 K温度下拉应力时程曲线的第一个峰值内系统内部发生了相变, 原有的BCC结构转变为FCC, HCP和未知结构, 在第一个峰值之后孔洞开始成核与生长.在 900—1100 K 温度下, 拉应力时程曲线只有一个峰值出现, 孔洞成核与生长是导致拉应力下降的主要原因.

2)随着温度的升高, 系统内部拉应力峰值降低, 最大拉应力为 100 K 时的 21.4 GPa, 最小拉应力为 1100 K 时的 13.7 GPa, 降低幅度为 35.9%;系统总势能随温度升高而增加, BCC结构转变为FCC, HCP和未知结构会导致系统势能降低, 但降幅不大, 孔洞快速生长是导致势能快速降低的主要原因.

3)孔洞成核在发生HCP, FCC相变之后与拉应力达到峰值之前(100—700 K是拉应力时程曲线第二个峰值之前)发生, 孔洞在未知结构区域内成核, 孔洞的生长推动周边原子的移动, 使FCC和HCP结构排列转变为新的BCC结构, 新产生的BCC结构不同于初始的BCC结构, 可能是变形孪晶的结果.

4)孔洞体积分数的增长可以分为三个阶段,第一阶段由孔洞成核主导, 此时孔洞体积分数增长较慢, 第二阶段由孔洞生长主导, 此时孔洞体积分数增长最快, 呈指数型增长, 第三阶段由孔洞贯穿主导, 此时孔洞体积分数呈线性增长, 速度适中.孔洞在高温下更容易成核与生长, 100 K温度下孔洞体积分数快速增长阶段S2的时间是1100 K温度下的2.5倍.孔洞的快速生长阶段与拉应力的峰值相对应, 因此高温下拉应力下降速率随孔洞体积分数快速增长阶段时间缩短而加快.

5) NAG 模型可以用来描述单晶铁在 5 × 109s—1应变率100—1100 K温度下孔洞体积非线性增长,在100—1100 K温度下拟合结果与MD结果都符合较好, 误差在 20% 以内.与低碳钢相比, 单晶铁的拉应力敏感系数与黏度更低.随着温度的升高,单晶铁的孔洞成核和生长阈值逐渐降低.

本文的分析结果有助于更全面地了解单晶铁在高速冲击下的动态响应, 研究结果可以为建立考虑温度影响的单晶铁在高应变率下的损伤模型提供参考.本文研究只考虑了温度的影响, 然而金属材料在高应变率下的动态响应受多个因素影响, 温度只是重要因素之一, 因此后续研究将考虑多种因素共同作用, 例如晶粒尺寸、预置孔洞、应变率等因素与温度共同作用对单晶铁的动态响应的影响.