CoCrFeMnNi 高熵合金冲击波响应与层裂强度的分子动力学研究1)

杜 欣 袁福平 熊启林 张 波 阚前华 张 旭 ,2)

* (西南交通大学力学与航空航天学院,成都 610031)

† (中国科学院力学研究所非线性力学国家重点实验室,北京 100190)

** (华中科技大学航空航天学院,武汉 430074)

引言

高熵合金作为新型结构材料,其拥有高硬度、高强度、良好的塑性、抗辐照性能以及出色的热稳定性和抗腐蚀能力[1-5],这为高熵合金未来在航空航天、太空探索以及深海探测等领域的应用提供了巨大的应用前景.作为一种经典的高熵合金[6],CoCrFeMnNi 高熵合金在极端环境下的力学行为研究备受关注.然而,现有的表征手段无法观测到动态连续的微结构演变,使得高熵合金在纳米尺度的层裂行为缺少深入的理解.分子动力学(molecular dynamics,MD) 方法具有时间尺度短,空间尺寸小,应变率大的特点[7],可以在极短时间尺度下研究材料冲击动态力学响应,且容易捕捉受冲击过程中的微结构演化,例如孔洞形核过程[8].因此,利用纳米空间尺度和皮秒时间尺度的MD 方法分析材料的动态力学行为具有天然优势.

目前,MD 模拟方法已成为研究材料冲击动态力学行为的重要手段之一,并已有大量的研究报道.Li 等[9]利用MD 研究了缺陷可逆性对碳化硅动态拉伸强度的影响.结果表明,在低于临界压缩应变的情况下,SiC 在沿[001]晶向加载时产生可逆的变形孪晶,这种可逆性使得晶体拥有高强度.在超过临界应变时,该过程变得不稳定,并触发孪晶界运动,产生不可逆的缠结缺陷,导致强度显著降低.Jian等[10]通过MD 模拟研究了冲击引起CoCrNi 中熵合金非晶化的阈值和原因,并与采用平均原子势的模拟结果进行对比,发现中熵合金由于晶格畸变效应而产生了更多的无序结构,造成其层裂强度降低.Cekil 等[11]研究了碳化硼在动态载荷作用下的力学性能,在沿[0001]和两个方向分别以2.0 和3.0 km/s 的速度进行冲击时碳化硼发生了非晶化，C-B-C 链的屈曲以及晶格旋转是发生非晶化的主要原因，其中方向不发生晶格旋转并且其C-B-C链发生屈曲的临界冲击速度大于[0001]方向.Zhu等[12]通过MD 研究了Cu/Ni 纳米多层体系在不同冲击脉冲持续时间下的动态响应和层裂,模拟结果发现层裂只出现在铜层区域,而不会出现在镍层区域;并且随着冲击持续时间的增加,层裂机制从空隙在铜层内均匀成核转变为在Cu/Ni 界面处成核,导致整体层裂强度下降.Xie 等[13]使用蒙特卡洛以及MD 方法研究了单晶和多晶CoCrNi 在冲击加载下的力学行为、缺陷演化和层裂破坏特征.结果表明单晶中熵合金的Hugoniot 弹性极限表现出异常的各向异性,并且孔洞成核于层错相交处或者冲击压缩产生的非晶态区域.在多晶中,孔洞形核则位于晶界处的Ni 聚集区域,一般为三叉晶界处,并随着晶界分离逐渐扩大.Thürmera 等[14]利用实验和MD 模拟研究了CoCrFeMnNi 高熵合金冲击波作用下层裂现象,发现层裂面产生了纳米孪晶,MD 模拟中证实CoCrFeMnNi 高熵合金拥有较Cu,Fe 和Ta 更高的层裂强度.

以上研究虽已对传统合金,中熵合金和高熵合金的冲击动态损伤进行了报道,但现有文献中尚未考虑弹塑性双波分离对层裂强度的影响.此外,在CoCrFeMnNi 高熵合金中沿不同方向冲击时的微结构演化以及微结构对层裂强度的影响仍需进一步分析.利用MD 模拟软件Lammps[15],研究冲击方向以及速度对单晶CoCrFeMnNi 高熵合金冲击波响应以及层裂行为的影响,并利用后处理软件OVITO 分析二者对冲击时微观结构演化的影响[16-17].相关研究成果有助于从原子尺度上理解CoCrFeMnNi 高熵合金的动态力学行为,促进CoCrFeMnNi 高熵合金未来作为抗冲击材料在航空航天等领域内的应用.

1 原子模型与加载方式

为研究冲击方向和冲击速度对层裂强度以及冲击过程中微观结构演化的影响,利用Atomsk 软件[18]构建了取向分别为的三个单晶CoCrFeMnNi 高熵合金模型,三个模型的尺寸均为74.7 nm × 13.3 nm ×7.5 nm,晶格常数为a=3.595 Å.其中x方向为[100]的单晶CoCrFeMnNi 高熵合金模型如图1 所示.在冲击之前对模型进行弛豫,使系统达到平衡状态.在弛豫过程中,模型的x,y和z方向均采用周期性边界条件,并在NPT 系综下弛豫50 ps (1 ps=1×10-12s).在诱导冲击波产生和冲击波传播阶段采用NVE 系综,并且冲击方向(x方向)为自由边界条件.采用活塞法诱导冲击波的产生[19],即在模型的一端设置虚拟墙作为活塞并给予一定的速度进行冲击.在本文中,为探究冲击速度对层裂强度的影响,施加的冲击速度Up为0.6,0.9,1.2,1.5 和1.8 km/s.虚拟墙作为活塞进行冲击的过程持续5 ps.之后,在NVE 系综下冲击波在模型内传播并在自由面发生反射形成拉应力区,该过程持续时间为35 ps.在整个模拟过程中,时间步长为1 fs (1 fs=1×10-15s),采用的势函数为Choi等[20]提出的第二近邻修正嵌入原子(2 NN MEAM) 势.该势函数已被用来模拟CoCrFeMnNi 高熵合金的冲击动态力学行为[14,19].此外,该势函数也被用来模拟CoCrFeMnNi 高熵合金的拉伸、循环、纳米压痕以及辐照行为[21-24].因此,该势函数具有一定可靠性.

图1 x 方向为[100]取向的单晶CoCrFeMnNi 高熵合金原子模型.黄色和红色平面分别代表虚拟墙和自由面Fig.1 Atomistic model of single-crystal CoCrFeMnNi high-entropy alloy with [100] orientation in x direction.The yellow and red planes represent virtual wall and free surface,respectively

2 结果与分析

2.1 自由面速度

自由面速度可以反映层裂强度和动态损伤的发生.当冲击波到达样品自由面时,自由面速度增大到最大值.随后,冲击波在自由面发生反射,反射后的稀疏波在与衰减冲击波尾部的稀疏波相遇后,自由面发生卸载,造成自由面速度下降并且形成拉应力区域.当拉应力达到阈值后发生层裂,而反射后的稀疏波作用于层裂面时会发生再次反射,使得自由面速度增加,产生速度回跳[10,25-26].

沿[100],[110]和[111]方向进行冲击加载时的自由面速度演化如图2 所示.可以看出,在以0.6 km/s 的速度沿[100]方向进行冲击时自由面速度未发生回跳,表明在该速度下沿[100]方向冲击未出现动态损伤[27].但是,沿[110] 和[111] 方向以0.6 km/s 的速度冲击时自由面速度发生了多次回跳,表明沿[110]和[111]方向冲击时动态损伤发生的临界冲击速度小于[100]方向.

图2 自由面速度演化Fig.2 Evolutions of free surface velocity

在虚拟墙以0.6 km/s 的速度沿三个方向冲击5 ps 后,样品内部的塑性变形未随着冲击波的传播而增加,表明此时只产生了弹性波,未产生塑性波.但是,由于在CoCrFeMnNi 高熵合金中拉伸强度低于压缩强度,表现出了明显的拉-压非对称性[22],使得在[110]和[111]方向下冲击波反射后形成的拉应力区产生了塑性变形,如图3 所示.其中,绿色原子为面心立方(face-centered cubic,FCC)结构;红色原子为密排六方(hexagonal close-packed,HCP)结构,并且主要为层错;蓝色原子为体心立方(bodycentered cubic,BCC)结构;白色原子为无序结构[28].采用了共邻分析(common neighbor analysis,CNA)方法[29]分析.可以看出,在虚拟墙冲击5 ps 时,由于模型冲击表面原子受到虚拟墙的作用发生畸变,使得冲击表面附近产生了少量层错.在冲击波正向传播的过程中(10,15 ps),塑性变形并未随着冲击波的传播而增加.但是,当冲击波在自由面发生反射后(20 ps),由于拉-压非对称性,造成了在拉应力区发生塑性变形,产生了层错.此时,[110]方向在当前的二维视角中未产生贯穿模型的孔洞.但随着冲击波的进一步作用,在25 ps 时[110]方向出现了贯穿模型的孔洞.

图3 采用CNA 方法分析以0.6 km/s 的冲击速度沿[110]和[111]方向冲击时的微观组织演化Fig.3 Microstructure evolutions by CNA at the 0.6 km/s shock velocity along the [110] and [111] directions

此外,从图2(b)和图2(c)中可以看出,随着冲击速度的增加,在沿[110]和[111]方向进行冲击时,自由面速度达到峰值前经历了两次增加的过程.其中,第一次增加是由于弹性波作为先驱波首先到达自由面,第二次增加是由于塑性波到达自由面.由于沿[100]方向未发生弹塑性双波分离现象,所以其自由面速度在达到峰值前只经历了一次增加的过程.[100] 方向未出现双波分离现象与该方向拥有较[110]和[111]方向更多可开动滑移系有关,多个滑移系的开动造成弹性波波速较低,甚至低于塑性波波速[30].当冲击速度为0.9 km/s 时,在沿[110]和[111]方向冲击时样品内部产生了塑性波,但此时在[111]方向的自由面速度中未体现出弹塑性双波分离,并且[110]方向弹塑性双波分离现象不明显,只对自由面速度造成微小扰动.此外,在以1.2 km/s 的冲击速度沿[111]方向冲击时,自由面速度曲线反映出的弹塑性双波分离现象同样不明显.随着冲击速度的进一步增加,在以1.2 km/s 和1.5 km/s 分别沿[110]和[111]方向冲击时,自由面速度在达到峰值前明显经历了两次增加的过程,表明了弹塑性双波分离现象的出现.由于沿[111]方向加载时Schmid因子小于[110]方向(分别为0.272 和0.408),所以沿[111]方向进行冲击时塑性变形所需的外加应力更大[19],导致了[111]较[110]方向拥有更大的弹塑性双波分离现象临界冲击速度.之后,弹塑性双波分离现象随着冲击速度的增加而减弱.

为进一步分析弹塑性双波结构,得到沿[100],[110]和[111]方向以0.6～1.8 km/s 冲击5 ps 时的粒子速度剖面图,如图4 所示.可以看出,在[110]和[111]方向以0.6 km/s 的速度冲击时曲线并未出现台阶段,表明此时未出现弹塑性双波分离现象.在[110]和[111]方向以0.9 km/s 的速度冲击时,曲线中出现了台阶段,但是台阶段的粒子速度与冲击表面粒子速度相差不大,使得在图3 中[110]方向的自由面速度未表现出明显的弹塑性双波分离现象,而[111]方向未出现弹塑性双波分离现象.在1.2 km/s的冲击速度下,[110]和[111]方向出现了明显的台阶段,但是由于弹性波的单独扰动,使原子面在波前附近处弹性振动[31],导致[111]方向弹性波前附近粒子速度上升并与冲击表面附近粒子速度相差仍然不大,造成图3 中[111]方向的自由面速度未表现出明显的弹塑性双波分离现象.由于台阶段粒子速度与冲击表面速度相差不大,而此时已发生双波分离,使得自由面受到了分离的弹性波和塑性波的两次作用,造成了在图3 中自由面速度第一次回跳点滞后.随着冲击速度的进一步增加,塑性变形加剧,使塑性波传播速度增加,造成塑性波与弹性先驱波之间的速度差减小,从而使得台阶段的宽度逐渐减小,表明在发生明显的弹塑性双波分离现象后,弹塑性双波分离现象随着冲击速度的增加而减弱.在[111]方向发生双波分离现象后(Up≥1.2 km/s),[111]方向先驱波的波前附近粒子速度明显大于[110]方向.此外,在[100]方向下以0.6～1.8 km/s 速度冲击时,曲线中均未产生台阶段,表明[100]方向未出现双波分离现象.

图4 冲击5 ps 时粒子速度剖面图Fig.4 Stress profile at 5 ps

整体而言,随着冲击速度的增加,[110] 和[111]方向的弹塑性双波分离现象先增强后减弱,自由面速度多次回跳现象减弱,自由面峰值速度增加.此外,沿[100]方向进行冲击时产生动态损伤的临界冲击速度大于[110]和[111]方向,表明了[100]方向拥有更强的抗冲击能力.

2.2 层裂强度

式中,ρ0和cB分别为初始密度和体波声速,ΔU=umax-upull-back为自由面速度演化曲线中峰值速度与第一个回跳点速度之差.在这里,初始密度ρ0=7.856 g/cm3[14],体波声速cB=4.5 km/s[33].

然而,在冲击速度超过Hugoniot 弹性极限时,在冲击压缩过程中已发生塑性变形,此时需要考虑弹塑性响应和层裂片厚度的影响,修正后的层裂强度可表示为[26,34-35]

此外,层裂强度还可通过冲击方向最大拉应力计算,即[35-36]

为了进一步分析层裂强度与冲击方向和冲击速度的关系,采用经典和修正声学近似方法以及最大拉应力方法,得到[100],[110]和[111]方向下层裂强度与冲击速度的关系,如图5 所示.

图5 采用三种方法得到沿[100],[110]和[111]方向冲击时层裂强度与冲击速度的关系Fig.5 Three methods are used to obtain the relationship between spall strength and shock velocity when the shocking along the [100],[110] and[111] directions

从图5 中可以看出,在沿[100]方向进行冲击时以三种方法获得的层裂强度都表现出随冲击速度的增加而减小的变化趋势.在[110]和[111]方向下,当处于弹塑性双波分离的临界速度时([110] 和[111]方向下分别为0.9 和1.2 km/s),三种方法获得的层裂强度都发生了明显减小.其中,在双波分离现象出现时,由于弹性波前粒子速度与塑性波速度相差不大,造成自由面速度曲线中的峰值速度下降,从而使得通过经典和修正声学近似方法计算层裂强度时自由面速度演化曲线中峰值速度与第一个回跳点速度之差减小,造成层裂强度大幅减小,使得在该速度下通过经典和修正声学近似方法计算的层裂强度小于更快冲击速度下的层裂强度.此外,由于双波分离现象的出现,使得在塑性波前出现了热尖峰,如图6所示.塑性波前的热尖峰会在塑性波反射后再次对模型的温度产生影响,从而使得通过最大拉应力得到的层裂强度在双波分离现象出现时也表现出了明显的减小趋势.

图6 以0.6 和0.9 km/s 的速度沿[110]冲击时的温度剖面图Fig.6 Stress profile in the [110] direction at the shock velocities of 0.6 and 0.9 km/s

在采用经典和修正声学近似方法计算层裂强度时所需要的一些关键物理量,如体波声速、初始密度、纵波声速均不是动态变化的,因此只能近似地得到动态力学行为中的层裂强度.然而,MD 可以实时提取冲击方向的应力分布.因此,采用在冲击加载时模型冲击方向出现的最大局部拉伸应力作为层裂强度更准确.在这里,采用最大拉伸应力来描述层裂强度与冲击方向和冲击速度之间的关系,如图7 所示.

图7 采用最大拉伸应力方法得到沿[100],[110]和[111]方向冲击时层裂强度与冲击速度的关系Fig.7 The maximum tensile stress method is used to obtain the relationship between spall strength and shock velocity when shocked along the [100],[110] and [111] directions

可以看出,沿三个方向冲击时层裂强度都表现出了随着冲击速度增加而减小的变化趋势.弹塑性双波分离的出现使得[110]和[111]方向的层裂强度大幅减小.整体而言,[100]方向拥有更大的层裂强度.在塑性变形较小时(Up≤0.9 km/s),沿[111]方向冲击时层裂强度大于[110]方向.但是,随着冲击速度的增加产生较大的塑性变形后(Up≥1.2 km/s),沿[110]方向冲击时层裂强度大于[111]方向.

为探究冲击速度以及冲击方向带来的层裂强度的差异,采用CNA 方法[29]分析以0.9 km/s 和1.5 km/s 的冲击速度沿[100],[110]和[111]方向冲击时的微观组织演化,如图8 所示.从图8(a) 和8(b) 中可以看出,在沿[100] 方向的冲击阶段(5 ps)产生了BCC 中间相,并且在高速冲击下会产生更多的BCC 中间相.BCC 中间相的产生与冲击方向和压力相关,其过程为沿[100] 方向加载时,FCC 结构中的体心四方(body-centered tetragon,BCT)结构转变为BCC 结构,进而表现出FCC 结构向BCC 中间相的转变,该过程遵循Born 准则[19,37].BCC 中间相会抑制层错和无序结构的产生[19].此外,由于BCC 相的结构不稳定,以及FCC 和BCC 两相间能垒随着压缩应力的减小而减小[37],使得在冲击波到达自由面并发生反射后(15 ps),在稀疏波的卸载作用下BCC 中间相转变为FCC 结构.

由图8(c)～图8(f) 可以看出,在沿[110] 和[111]方向的冲击阶段主要产生了层错以及无序结构,并且产生的无序结构随着冲击速度的增加而增加.在层裂发生初期阶段(30 ps)，微孔洞形核于大量无序结构处[38],表明该区域发生了非晶化.此外,随着冲击速度的增加,产生了更多的微孔洞,并且微孔洞的分布更广.

图8 采用CNA 方法分析以0.9 km/s 和1.5 km/s 的冲击速度沿[100],[110]和[111]方向冲击时的微观组织演化Fig.8 Microstructure evolutions by CNA at the shock velocities of 0.9 km/s and 1.5 km/s along the [100],[110] and [111] directions

非晶化是由拉伸应力和温升共同导致的[13],为了探究非晶化与拉伸应力和温度的关系,得到沿[100]方向以0.9 km/s 的速度冲击时的应力剖面图与温度演化过程,分别如图9 和图10 所示.可以看出,在冲击过程中(5,10 ps)模型处于压应力状态,压应力和温升随着波前沿冲击方向的传播而出现在更深的区域.随后,冲击波到达自由表面并发生反射,反射后的稀疏波在与衰减冲击波尾部的稀疏波相遇后(20 ps),形成拉应力区域,并且在拉应力区中产生微孔洞以及发生非晶化.此外,在非晶化区域中拥有更高的温度.

图9 沿[100]方向以0.9 km/s 的速度冲击时的应力剖面图Fig.9 Stress profile in the [100] direction at the shock velocities of 0.9 km/s

图10 沿[100]方向以0.9 km/s 的速度冲击时的微结构和温度演化Fig.10 Microstructure and temperature evolutions in the [100]direction at the shock velocities of 0.9 km/s

层裂强度与非晶化程度密切相关[10].在这里,提取自由面速度达到峰值时的无序结构含量,如图11所示.可以看出,[100]方向产生的无序结构含量最少,这使得微孔洞较[110]和[111]方向更难以形核,所以[100] 方向拥有较大的层裂强度.[110] 和[111]方向产生了较多的无序结构,并且[111]方向产生的无序结构含量小于[110]方向.但是,只有在较低速度冲击(Up≤0.9 km/s)时[111]方向的层裂强度大于[110]方向.此外,随着冲击速度的增加,模型内温度增加(图6),甚至局部温度超过熔点,这使得在模型内发生了非晶化,导致产生的无序结构含量随着冲击速度的增加而增加.

图11 自由面速度达到峰值时无序结构含量Fig.11 The content of disordered structure at the peak free surface velocity stage

[110]和[111]方向层裂强度大小关系的转变与微孔洞形核区域的无序结构含量有关[10].提取了沿冲击方向宽度为10 nm 的微孔洞形核区域中无序结构含量,如图12 所示.可以看出,在Up≤0.9 km/s 时,[110] 方向在层裂发生初期的微孔洞形核区域产生了更多的无序结构.但是,在Up≥1.2 km/s 时,[111]方向产生了更多的无序结构.由于在[111]方向双波分离现象出现后,弹性先驱波的波前附近粒子速度明显大于[110] 方向(图4),导致[110] 和[111]方向下层裂初期微孔洞形核区域中无序结构含量大小关系的转变.此外,多主元使得高熵合金拥有晶格畸变效应并产生更多的无序结构,进而降低了层裂强度[10,13].

图12 沿[110]和[111]方向冲击时微孔洞形核区域的无序结构含量Fig.12 The content of disordered structure in the microvoid nucleation region along the [110] and [111] directions

3 结论

本文对单晶CoCrFeMnNi 高熵合金进行了冲击加载的MD 模拟,研究了冲击方向与冲击速度对冲击波响应、层裂强度以及微结构演化的影响,并从微观结构演化的角度分析了造成层裂强度改变的原因,获得以下主要结论.

(1) 沿[100]方向进行冲击时自由面速度未表现出弹塑性双波分离现象,而以一定速度沿[110]和[111]方向冲击时产生了双波结构,并且双波分离现象随着冲击速度的增加表现出先增强后减弱的变化趋势.

(2) 在层裂发生的初期,微孔洞的形核产生了大量无序结构.随着冲击速度的增加,模型内温度增加,在高冲击强度和高温下产生的无序结构更多,从而使得沿[100],[110]和[111]方向冲击时层裂强度随着冲击速度的增加而降低.此外,双波分离使得[110]和[111]方向的层裂强度下降明显.

(3) 由于BCC 中间相的存在阻碍了无序结构的产生,使得[100]方向拥有较高的层裂强度.在冲击速度较小时(Up≤0.9 km/s),沿[111]方向冲击的层裂强度大于[111] 方向.但是,在冲击速度较大时(Up≥1.2 km/s),沿[110]方向冲击的层裂强度略大于[111]方向.层裂初期微孔洞形核区域无序结构含量大小关系的转变,造成了[110]和[111]方向下层裂强度大小关系的转变.