<100>拉伸载荷下钽单晶弹性特性的温度和压力依赖性：分子动力学研究

背景

一般来说，钽（Ta）在常温下属于 BCC 结构。目前已有大量文献证明Ta单晶在高压下表现出优异的相稳定性[1,2,3]。更重要的是，Ta 在环境压力下具有非常高的熔化温度 3269 K，高于大多数其他金属 [4]。由于其优异的性能，Ta是许多技术应用的理想材料，例如微/纳米电子学中的扩散屏障、磨损保护涂层和高温超合金。

最近，在实验 [2, 5,6,7] 和理论 [8,9,10,11,12,13,14] 领域的大量努力已经投入到塔。德瓦勒等人。 [5] 研究了压力对金刚石砧室 (DAC) 中 Ta 屈服强度的影响，最高可达 93 GPa，DAC 实验还表明 BCC 结构在高达 135 GPa 时仍保持稳定 [2]。此外，Shigeaki [8] 使用 DFT 模拟了高达 100 GPa 和 3000 K 的 Ta 的状态方程 (EOS)。吴等人。 [9] 研究了 Ta 在高达 350 GPa 的高压下的弹性和热力学特性。与此同时，Škoro 等人。 [6, 7] 分别在高达 2250 和 2500 K 的极高温度下测量了 Ta 的屈服强度和杨氏模量。顾等人。 [10] 使用第一性原理方法对高达 500 GPa 的立方体 Ta 的高压结构和弹性性能进行了研究。发现作为压力函数的弹性常数以及Ta的体积、杨氏和剪切模量都随着压力的增加而增加。

除了DAC实验和DFT计算，在MD模拟领域也有很多在高温高压下的研究[15,16,17,18]。刘等人。 [15] 使用扩展的 Finnis-Sinclair (EFS) 势并研究了热 EOS 以及 Ta 在高达 400 GPa 的压力下的熔化特性。此外，Tramontina 等人。 [16, 17] 研究了晶体取向对高压塑性机制的影响。他们还讨论了冲击强度和冲击上升时间对其微观结构的影响。此外，Ruestes 等人。 [18] 使用三种不同的原子间势对 BCC Ta 进行了压痕模拟，并提出了导致塑性变形区产生和扩展的缺陷机制​​。

尽管进行了大量研究，但还没有使用 MD 模拟对拉伸载荷下 Ta 的动态响应进行系统的原子模拟研究。本工作的主要目的是研究 Ta 单晶在 <100> 拉伸载荷下的弹性特性，同时考虑尺寸、应变速率、温度和压力的影响。此外，了解是否可以通过<100>拉伸载荷诱导相变是这项工作的另一个目的。

方法/实验

物理建模

如图 1 所示，本文研究的 Ta 立方是通过沿 <001>、<010> 和 <100> 方向重复 BCC 晶胞生成的，晶格参数为 a =b =c =3.301 Å，分别。四种不同边长的三次模型，包括 12a (3.96 纳米), 18a (5.94 纳米), 24a (7.92 nm) 和 30a (9.90 nm)，构建。对应的原子数分别为 3456、11664、27648 和 54000。图 1 显示了边长为 3.96 nm 的 Ta 立方体的示意图，这是我们目前模拟中的原始结构。

边长为 3.96 nm 的 Ta 立方体的 BCC 结构和示意图。草图是我们目前模拟中的原始结构

MD 模拟详情

我们在本文中报告的 MD 模拟是使用大规模原子/分子大规模并行模拟器 (LAMMPS) [19] 进行的。在MD模拟过程中，根据这些原子之间的原子间势函数计算施加在Ta原子上的外力。在这项研究中，考虑了两种不同的原子间势：Zhou 等人的 EAM 势。 [20] 和 Ravelo 等人。 [21]，分别。两种势均应用于再现表 1 中列出的 Ta 的一些弹性特性。为简单起见，将它们称为 Ravelo-EAM 和 Zhou-EAM 势。

<图>

三次对称结构的弹性行为完全由它们的弹性常数 C 描述 11、C 12、C 44. <100>、<110> 和<111> 单晶的依赖于取向的弹性模量通过几个方程计算[18]。在我们的工作中，进行了 MD 模拟以研究 <100> 拉伸载荷下对弹性模量的影响。因此，我们关注 <100> 取向的弹性模量并取 E 的弹性模量 100 计入。因此，我们考虑弹性常数 C 11 和 C 12 以及以下等式[22]：

$$ {E}_{100}=\frac{\left({C}_{11}-{C}_{12}\right)\left({C}_{11}+2{C}_ {12}\right)}{\left({C}_{11}+{C}_{12}\right)} $$ (1)

通常，从 MD 计算弹性常数的方法有 3 种，包括应力波动法、应变波动法和直接法。在目前的工作中，类似于高等人的直接方法。 [23] 用于计算弹性常数C 11 和 C 12 对于两种类型的 EAM 电位，如图 2 所示。

C 弹性常数的压力相关性 11 和 C 12 基于所研究的电位。 一 周-EAM。 b 拉夫洛-EAM。在目前的工作中，我们利用直接法计算弹性常数C 11 和 C 12 对于两种类型的 EAM 电位。两种电位的 C11/C12 压力曲线如图 a 所示 和 b , 分别

从图 2a 中，可以发现获得的 C 曲线 11 和 C 12 基于 Zhou-EAM 势的压力增加到 140 GPa 不能保持平稳。在 ~ 40 GPa 的压力下有一个交叉点，即 C 11 =C 12，其中 E 的弹性模量 根据等式，100 将为零。 (1).此外，E 100 当压力高于 40 GPa 时会呈现负值，这很可疑，与理论和实验结果相冲突 [24, 25]。因此，Zhou-EAM 潜力在此处考虑的范围内表现不佳。然后，让我们从 C 的 MD 模拟结果来检验 Ravelo-EAM 势的可行性 11 和 C 12 基于图 2b 所示的 Ravelo-EAM 电位。数值结果表明，压力越高，C的值越大 11 和 C 图 12，这与通过 DFT 计算得出的弹性常数随压力变化的趋势非常吻合 [9, 25]。同时，使用 Ravelo-EAM 势计算的结果与 Ruestes 等人报告的值显着相同。 [18]。 Ravelo-EAM 势在高压下表现良好，同时它还可以描述动态变形下 Ta 的弹性和机械性能 [26]。因此，我们将在以下部分中基于 Ravelo-EAM 势进行模拟。

在几何构建之后，我们进行了一系列相关的 MD 模拟。在 MD 模拟期间，周期性边界条件 (PBC) 用于三次模型的所有三个方向。时间步长设置为 1 fs，系统温度设置为 1、300、600、900、1200 和 1500 K，以探索 Ta 的弹性特性的温度依赖性。首先，使用规范集成 (NVT) MD 模拟将模型松弛约 50 ps 松弛过程，以使系统处于局部电位最小值。然后，它使用等温-等压 (NPT) MD 模拟来确保指定的静水压力范围为 0 到 140 GPa，以研究压力对 Ta 弹性性能的影响 [27]。最后，应变率范围为 5 × 10 ⁸ 的拉伸载荷 s ^{− 1} 到 7.5 × 10 ⁹ s ^{− 1} [28, 29] 应用于 Ta 立方的 x 方向。同时，NPT 模拟在 y 和 z 方向上以第二步中施加的相同压力进行。因此，此处计算的弹性模量适用于 <100> 取向。对于所有 MD 模拟，模型将通过 <100> 拉伸载荷在 x 方向拉伸至 15%。

结果与讨论

拉伸过程

在拉伸过程中，使用科学软件包 Open Visualization Tool (OVITO) [30] 对模拟配置进行可视化。零压力下<100>单轴拉伸应变下Ta的应力-应变曲线及不同应变下对应的原子构型如图3所示。

Ta在零压力下的应力-应变曲线以及不同应变下的相应配置。注：蓝、绿、白球分别对应BCC、FCC等晶格结构

如图 3 所示，从应力-应变曲线可以得出结论，表面在配置 (IV) 附近断裂。在单轴拉伸应变开始时，应力通常随应变线性变化，图 3I 显示原子构型保持 BCC 结构。随着应变的增加，从 BCC 到 FCC 结构的相变开始于~ 7.4% 的应变，并在~ 9.8% 的应变时完成，分别如图 3II、II 所示。并且这些FCC结构在第一次表面断裂之前保持最大值。当应变为~ 13.1% 时，y 和 z 方向的边缘长度突然减小，导致表面断裂。同时，值得注意的是，集群发生的时间很短，如图 3IV 所示。在连续单轴变形下，原子构型保持条状直到应变~ 13.3%，如图3V所示。

在本文中，我们专注于拉伸特性对模型尺寸、应变率、温度和压力的依赖性，如本节所述。理论上，模型在弹性变形阶段被线性拉伸，弹性模量定义为应力应变曲线线性部分的斜率。可以发现，所有模型的拉伸过程相似，应力-应变曲线具有相似的变化趋势。因此，我们使用相同的方法获得了不同模型尺寸和应变率下Ta的弹性模量。

依赖于大小和应变率

表2列出了不同模型尺寸在1K温度和5 ×10 ⁸ 应变速率下的弹性模量和屈服应力 s ^{− 1} .可以很容易地得出结论，模型尺寸对 Ta 的弹性模量和屈服应力没有影响。这很容易解释，弹性模量是描述原子之间的相互作用，而弹性模量不随模型大小而变化。从表 2 可以看出，弹性模量为 ~ 139 GPa，与 140 GPa 的模拟结果显着相同 [18]。

<图>

根据现有的参考文献。 [28, 29]，大部分应变率范围在10 ⁸ s ^{− 1} 到 10 ¹⁰ s ^{− 1} .本文选取4种应变率进行拉伸模拟，包括5.0 × 10 ⁸ s ^{− 1} , 7.5 × 10 ⁸ s ^{− 1} , 5.0 × 10 ⁹ s ^{− 1} , 和 7.5 × 10 ⁹ s ^{− 1} .表 3 列出了在 300 K 温度和不同应变率下的弹性模量和屈服应力。很容易得出结论，应变速率对弹性模量和屈服应力没有明显影响。

<图>

同时，我们还模拟了模型尺寸和应变速率对不同温度和压力下的弹性模量和屈服应力的影响。这些模拟得出相同的结论。因此，我们将使用相同的模型尺寸 3.96 nm 和相同的应变率 5 × 10 ⁸ s ^{− 1} 用于以下模拟。

对温度的依赖

图 4 显示了高达 1500 K 的不同温度下的应力-应变曲线。可以看出，这些曲线的斜率表示 <100> 方向的弹性模量，即 E 100，在弹性拉伸期间，屈服应力随着温度的升高而逐渐减小。根据热力学理论[31]，体系所有原子的总动能一般满足以下方程：

$$ {E}_{\mathrm{k}}=\sum \limits_{i=1}^N\frac{1}{2}{mv_i}^2=\frac{3}{2}{Nk} _BT $$ (2)

其中 E k 是系统的总动能； N 是原子总数； K B 是玻尔兹曼常数； T 是热力学温度。因此，可以得出结论，温度越高，系统包含的总动能越大，原子运动得越快。从热力学的角度来看，原子变得更加活跃，原子的运动更加强烈，这意味着其平衡位置的振幅更大。在拉伸过程中，原子之间的吸引力相对减小，原子容易从平衡位置逸出，因此在相同应变下x方向的应力减小。因此，较高温度下的弹性模量将小于较低温度下的弹性模量。此外，这些曲线的趋势与Ta的早期理论和实验结果相吻合[6, 7, 32]。

3.96 nm Ta模型在应变速率为5 × 10 ⁸ 时的应力-应变曲线 s ^{− 1} 以及从 1 到 1500 K 的不同温度。从图 4 中可以看出，这些曲线的斜率表示 <100> 方向的弹性模量，即 E 100，在弹性拉伸期间，屈服应力随着温度的升高而逐渐降低

为便于观察，表4列出了Ta在不同温度下的弹性模量和屈服应力。随着温度从 1 升到 1500 K，弹性模量将从 136.49 到 76.67 GPa 降低 ~ 42.3%，屈服应力从 ~ 8 到 ~ 4 GPa 降低 ~ 51%。

<图>

根据表 4，我们可以进一步推导出弹性模量 (E 100) 结果显示为图 5a 中的实线。参数化如下：

$$ {E}_{100}=a+{b}^{\ast }T $$ (3)

其中 E 100 在 (GPa) 和 T 在 (K ); 一 =138.07 ± 0.92111，b =- 0.04094 ± 0.00101。这个等式表明 E 100 随着温度的升高而线性下降，建议在 0 到 1500 K 的温度范围内使用。 (3) 容易得到E 在 T 的温度下，100 将达到 0 GPa 临界 =−a /b =3372 K，非常接近 Ta 的熔化温度 [15]。

一 弹性模量和b Ta 的屈服应力与温度的关系。在图 5 中，我们还展示了弹性模量 (E 100) 结果

对于屈服应力，推荐的参数化为

$$ {Y}_{\mathrm{stress}}=a+{b}^{\ast }T+{c}^{\ast {T}^2 $$ (4)

其中 Y 压力在 (GPa) 和 T 在 (K ); 一 =7.99610 ± 0.0415，b =− 0.0039 ± 1.30126 × 10 ^{− 4} , 和 c =7.97381 × 10 ^{− 7} ± 8.32307 × 10 ^{− 8} .从方程。 (4) 可以发现，屈服应力可能随着温度的升高而降低，符合二次多项式模型，如图 5b 中的实线所示。

对压力的依赖

正如“引言”部分所述，人们对高压条件下 Ta 的热弹性特性进行了大量的理论和实验研究。与静态方法不同，我们采用动态方法通过 <100> 拉伸加载来检查 E 弹性模量的压力依赖性 100在不同的静水压力下。

图 6 显示了 E 的弹性模量曲线 在 1 到 1500 K 的不同温度下，100 对高达 140 GPa 的压力。虽然所有不同颜色的实线都是在不同温度下通过动态方法获得的，但带有方形标记的虚线是通过静态方法通过使用方程获得的。 (1) 基于C的值 11 和 C 12 at 0 K. 很明显，1 K 处的曲线，即带圆圈标记的红色实线，与通过静态方法（虚线）获得的曲线几乎重叠，这表明在目前的工作适用于高达 140 GPa 的高压。

Ta 在不同温度和压力下 <100> 取向的弹性模量。值得指出的是，这些曲线的弹性模量为E 100

如图 6 所示，E 的弹性模量 当压力从 20 GPa 增加到 60 GPa 时，100 在不高于 600 K 的温度下显示出向下的凹形截面。 Ruestes 等。 [18] 通过使用 MD 模拟报告了作为压力函数的弹性常数，高达 60 GPa，结果与计算的 C 非常吻合 11 和 C 12 在目前的工作。作为回报，E 的计算弹性模量 100 也显示出与我们的结果相同的趋势。但是E的计算弹性模量 DFT 计算 [33] 中的 100 随着压力的增加而逐渐增加，曲线中没有发现下凹部分。是什么导致了这种不一致？通常，MD 模拟中使用的电位是通过拟合 DFT 计算和实验结果获得的。从这个意义上说，DFT 计算比 MD 方法具有更高的精度。 Ravelo-EAM 势 [21] 是通过将高压特性引入 Ta 单晶的拟合 DFT 状态方程 (EOS) 曲线来构建的。在拟合过程中，EOS 的冷曲线被扩展为包括晶格常数的高阶（三次和四次）项，这使得 MD EOS 对晶格常数的高阶项非常敏感。换句话说，这种不一致可能是由于 Ravelo-EAM 势不能准确描述在 20 到 60 GPa 压力下 Ta 的 EOS。综合来看，弹性模量-压力曲线在不同温度下具有相似的趋势，且弹性模量在~ 40 GPa以上随着压力的增加而逐渐增加。

结论

在本文中，进行了 MD 模拟以研究通过 <100> 拉伸载荷的 Ta 单晶弹性特性的温度和压力依赖性。首先，我们对两种类型的 EAM 势，包括 Zhou-EAM 和 Ravelo-EAM，根据 Ta 在 0 K 和不同静水压力下的弹性特性进行了比较研究。结果表明，在不同压力下，Ravelo-EAM 势比 Zhou-EAM 势表现更好。然后，基于 Ravelo-EAM 势对 Ta 单晶的拉伸行为进行 MD 模拟。观察表明，Ta在<100>拉伸载荷下断裂前会经历BCC-FCC相变，模型尺寸和应变速率对Ta单晶的拉伸行为没有明显影响。此外，E 的弹性模量 随着温度从 1 到 1500 K 的升高，100 会从~ 136 线性下降到~ 79 GPa，随着温度的升高，屈服应力从~ 8 到~ 4 GPa，符合二次多项式公式。最后，弹性性能的压力依赖性在 0 到 140 GPa 之间进行，观察表明弹性模量随着整体压力的增加而增加。 MD 模拟结果还表明 Ravelo-EAM 势能在较高压力下表现良好以及 E 的计算公式 100 使用 C 11 和 C 12 在低于 140 GPa 的压力下。

更改历史

缩写

密件抄送：

体心立方

DAC：

金刚石砧座

DFT：

密度泛函理论

EAM：

嵌入原子法

EFS：

扩展的芬尼斯-辛克莱

EOS：

状态方程

FCC：

面心立方

LAMMPS：

大规模原子/分子大规模并行模拟器

MD：

分子动力学

NPT：

等温-等压系综

NVT：

规范合奏

OVITO：

打开可视化工具

PBC：

周期性边界条件

Ta：

钽