基于分子动力学方法的复杂合金Grüneisen状态方程参数计算及验证：以GH4169合金为例

张圣来，丁 军，童 权，黄 霞，宋 鹍，路世青

（重庆理工大学 机械工程学院，重庆 400054）

1 研究进展

当前，航天航空等领域对性能优异的复杂合金的需求越来越大，对其使用环境要求也更加严格。例如，航空发动机在工作中一直处于高温高压状态，因此选择材料时必须重点考虑极端服役环境下的合金性能。

在超高压下，应力偏量几乎可以忽略，此时材料的本构模型退化为体积应力关系，即高压物理领域中的状态方程。Mie-Grüneisen状态方程作为应用最广泛的一种状态方程，通常用于描述物质在极端条件下的状态，如冲击波传播。在这些极端条件下，物质会处于高温高压状态，性质发生变化。近几十年来，等温和绝热状态方程研究取得了很大进展。等温状态方程有两种：一种是在0 K，也叫作冷曲线；另一种是在300 K。0 K等温状态方程一般根据300 K下的冲击波数据通过温度修正得到。300 K等温状态方程是静态状态方程，通过静态压缩保持温度不变。Bridgman［1］开发了一种静态压缩装置，用于探索材料在高压下的压缩特性。由于压缩过程非常缓慢，材料内部产生的热量能与外界完全交换，因此可将该压缩过程认为是等温压缩。在此基础上，许多学者对材料进行高压等温压缩。张剑等［2］采用原位高压同步辐射X射线衍射技术在室温、21.5 GPa的压力范围测定铝的等温状态方程，利用Murnaghan方程对实验结果进行数值拟合，得到铝的零压体弹模量及其1级压力导数。该结果同相关文献资料研究结果在误差范围内符合得很好。Syassen等［3］采用12 GPa的高压X射线衍射测定Ag和Al的等温压缩曲线，并选择Birch方程作为等温状态方程来拟合实验数据。Voˇcadlo等［4］和Shanker等［5］通过比较各种等温状态方程，研究了等温体积模量及其1阶压力导数。Wang等［6］使用经典平均场方法计算了压力高达1 TPa时5种金属Al、Cu、Ta、Mo和W的293 K等温线。Hrubiak等［7］采用双参数Birch-Murnaghan方程描述了Hf在室温下的等温线，得到了67 GPa压力下α相Hf的P（E，V）型状态方程。Wang等［8］用平均场势法对Al的静态和动态压缩及体积模量进行全面研究，检验了Hugoniot曲线和300 K等温线。Akahama等［9］通过X射线衍射实验，利用Birch-Murnaghan和Vinet状态方程描绘了Bi、Pt和Au的300 K等温线。Verma等［10］计算了Al、Ta、Mo和W的固态和液态Hugoniot曲线，讨论了0 K等温线和P－T曲线。

绝热状态方程通常基于Hugoniot数据得到，而Hugoniot数据则来源于通过轻气炮或者爆炸载荷驱动的飞片撞击实验。Al’tshuler等［11］研究了高压区Al、Cu和Pb的状态方程。Zhang等［12］研究了93 W合金的Hugoniot曲线、冲击温度和多项式状态方程。Dewale等［13］研究了Ta、Au、Pt、Al、Cu和W在94 GPa以上的静态和动态状态方程。Bringa等［14］和Agarwal等［15］对Mg、Cu单晶的绝热冲击线进行非平衡分子动力学（NEMD）研究，得到了不同晶向下冲击波速度与粒子速度的线性关系。Mackenchery等［16］通过分子动力学（MD）讨论了4种不同原子间势对单晶Ti在0.5～2.0 km／s碰撞速度下的线性关系和Hugoniot曲线的影响。

Mie［17］和Grüneisen［18］提出了描述压力和体积之间关系的固体理论。目前，Mie-Grüneisen状态方程作为高压下的固体模型得到了广泛应用，许多仿真软件可以通过直接输入Mie-Grüneisen状态方程的参数来定义材料。沈天乐等［19］通过实验得到了3种不同密度的环氧树脂－SiO2空心微球复合材料的Hugoniot关系，讨论了基于Hugoniot关系得到的Mie-Grüneisen状态方程，并与P-α模型进行了比较，显示结果一致。Walsh等［20－21］通过爆炸系统测量了金属的状态方程，并估计了Grüneisen系数γ的体积依赖性。林华令等［22］给出了一种用冲击压缩数据计算Grüneisen物态方程参数的方法，计算了23种金属的Grüneisen系数，使Grüneisen系数沿实验冲击绝热线的计算值与晶格振动理论的计算值达到最优拟合。Mc-Queen等［23］通过平面波爆炸系统在100～200 GPa压力下获得了19种金属元素的状态方程。通过Mie-Grüneisen理论，Hugoniot关系（P，V，E）被扩展为更完整的状态方程（P，V，E，T）。吴绍曾等［24］利用刚性离子模型计算NaCl的Grüneisen参数及其随压力的变化，计算中对最近邻排斥势分别使用两种修正Born-Mayer势，并利用电子气理论的计算结果处理次近邻排斥势和关联势的影响，最后得到的结果与实验相符。Nellis等［25］通过80～440 GPa压力范围内的Hugoniot曲线测量Al、Cu和Ta的状态方程，并讨论了冲击温度和γ。Jin等［26］提出了基于绝热冲击线数据表征金属的0 K等温线的热力学公式，并使用Born-Meyer势预测了7种金属的0 K等温线和Grüneisen状态方程。

从上述研究中不难发现，Mie-Grüneisen状态方程可由等温线或Hugoniot曲线推导得到，因而两者均可作为参考线。虽然已推导出了Mie-Grüneisen状态方程的公式，但很难通过实验获得状态方程的所有参数。确定这些参数的手段除实验外还必须借助理论分析和分子动力学模拟方法。Yang等［27］使用分子动力学方法全面研究了Hugoniot曲线、等温线、等熵线及其内能，通过分子动力学模拟获得了单晶Al和Pb的Mie-Grüneisen状态方程，计算并讨论了沿Hugoniot曲线和等熵线的温度、沿Hugoniot曲线的熵增量和Grüneisen系数γ，研究结果与已有文献吻合较好。Chen等［28］通过分子动力学模拟研究了TiAl合金的Hugoniot曲线和Mie-Grüneisen状态方程，发现材料成分严重影响冲击波速与质点速度的线性关系、Hugoniot曲线和内能。Zhang等［29］使用分子动力学研究单晶镍的冲击行为。模拟得到的P－V曲线略高于实验数据，实验样品为多晶Ni且内含缺陷，因此导致了数据差异。Zhang等［30］提出了一种新的与温度相关的能热力学模型，用以研究材料的温度相关Grüneisen。通过分子动力学对Ni、Cu、Au进行分子动力学模拟，验证了新建立的热力学模型。

目前对于等温线、Hugoniot曲线和Mie-Grüneisen状态方程的研究较多，材料类型主要集中在单一元素和二元合金，对多元素的复杂合金的研究相对较少。众所周知，复杂合金的物理性质与其组成元素的物理性质有很大不同，且服役环境通常较复杂。因此，有必要研究复杂合金在极端条件下的动态力学性能。GH4169作为一种镍基高温合金，是航空发动机叶片的主要材料。该材料在高温下仍具有很好的力学性能，因此一直被研究人员所关注。

然而，由于缺少相应的势函数，分子动力学目前无法通过模拟直接得到复杂合金Mie-Grüneisen状态方程的各个参数，因此在Ji等［31］和Millett等［32］的工作基础上，提出一种用于计算复杂合金的Mie-Grüneisen状态方程参数，并以GH4169镍基合金为例计算了Mie-Grüneisen状态方程。

2 理论框架

2.1 复杂合金的C0和λ

由于目前已有的势函数中，包含4种元素以上的势函数较少，因此使用单一的势函数无法直接计算出复杂合金us和up的关系。Ji等［31］在对累计滚焊Al／Ni复合材料的冲击诱导反应特性进行研究时，通过纯Al、纯Ni的参数C0和λ，以质量占比为权重，拟合Al／Ni复合材料的C0和λ，从而得到Al／Ni复合材料的Grüneisen状态方程。结果证明理论结果与实验结果符合良好，验证了所提方法的正确性。Millett等［32］通过平板撞击实验得到TiAl合金Ti-46.5Al-2Nb-2Cr的Hugoniot数据，并通过合金中各元素的参数C0和λ，以质量占比为权重拟合了Ti-46.5Al-2Nb-2Cr合金的C0和λ。结果表明，通过这种计算方式获得的Hugoniot参数与实验测得的参数误差较大。鉴于结论存在矛盾，本文中提出一种新的思路：以质量分数为权重，通过合金中各稳定相的Hugoniot参数，对合金的Hugoniot参数进行拟合：

其中：mPi代表第i种相的质量分数；下标“A”和“Pi”分别代表合金和第i种相；n是合金中存在的相的总数。

2.2 Hugoniot曲线及其内能

Hugoniot曲线由质量、动量和能量守恒方程决定，其表达式为：

其中：ρ为密度；V＝1／ρ为受冲击材料的比体积；up为粒子速度；us为冲击波波速；P为压力；E为比内能。下标“0”和“H”分别表示这些量处于未受冲击的初始态和Hugoniot状态。

常温（300 K）下，在某个很宽泛的压力范围内，us和up的关系是线性的：

其中：C0是压力为0时的体积声速；λ是材料常数。这些量是通过拟合us和up得到的。

由式（3）～（6）可以推导出Hugoniot曲线PH（V）和内能EH（V）：

其中：cV是定容比热；T是冲击温度。当EH远大于E0时，可以忽略E0。

2.3 冷压和冷能

0 K下，us0K与up0K之间线性关系和Hugoniot曲线PH0K（V）的表达式为：

通常，式（6）是在常温（300 K）下通过实验得到，因此式（10）中的C′0和λ′需要基于C0和λ进行温度修正。

在0～300 K范围内进行温度修正时，体积膨胀系数αv可以近似地看作一个常数。V0K、ρ0K、C′0和λ′可以写成：

其中：γ（V0）是300 K时的Gruneisen系数。从式（11）～（15）可以算出C′0和λ′，此时可以求得式（10）的Hugoniot曲线PH0K。

Born-Meyer势可以精准描述离子晶体和金属的可压缩性，因此常用于计算冷压Pc－BM和冷能Ec－BM，其表达式为：

其中：Q和q是材料常数；ρ0K是0K时的材料密度；δ＝V0K／V。

Morse势同样可以描述高压0 K条件下原子间的相互作用，因此Morse势的冷压Pc－M和冷能Ec－M可以表达为：

其中，A和B是材料常数。

一般，假设PH0K、等熵线PS（V）和Pc（V）在初始点（P＝0，V＝V0K）的1阶导和2阶导近似相等：

由式（8）～（21）可得材料常数Q、q、A和B的表达式为：

2.4 Grüneisen状态方程

Grüneisen状态方程是高压物理和爆炸力学领域描述固体行为最常用的状态方程之一。作为特征方程，它描述了晶格热振动的贡献。Grüneisen状态方程最常用的参考曲线是Hugoniot曲线，其经典形式为：

经典形式以冷压Pc和冷能Ec作为参考线，但常温下获得的Hugoniot曲线PH（V）及其内能EH（V）同样可以作为参考曲线。替换参考曲线后的状态方程为：

Grüneisen系数量纲为1，表达式为［33］：

其中：t代表表征γ的不同模型，当t＝0、1、2，式（28）分别表示Slate模型γS，Dugdale-MacDonald模型γDM和自由体积模型γf。

PS在初始点（P＝0，V＝V0）近似满足式（28），因此式（28）也可以写成：

当t＝0、1和2时，有：

假设PS与PH的1阶导和2阶导在初始点（P＝0，V＝V0）处恒等，而P′H（V0）和P″H（V0）又可直接由实验或者分子动力学模拟得到的PH求导得到，那么当获得式（2）中的参数，可以得到：

将式（33）～（35）代入式（29）可得：

Grüneisen系数γ仅是V的函数，且γ＞0。γ的表达式一般是通过实验数据拟合得到的经验公式，而由于本文中研究在无实验条件下仅通过数值模拟得到复杂合金的Grüneisen状态方程，因此可用的经验公式为：

基于式（1）～（2）、（6）～（8）、（17）和（40）～（41），最终可以得到复杂合金的Grüneisen状态方程。

3 计算方法

采用的势函数有Mishin等［34］开发的Ni的嵌入原子法（EAM）势，Zhang等［35］开发的Ni-Nb的嵌入原子法（EAM）势，Ko等［36］开发的Ni-Ti的修正嵌入原子法（MEAM）势，Purja等［37］开发的Ni-Al的嵌入原子法（EAM）势，Howells等［38］开发的Cr的角相关势（ADP），Chamati等［39］开发的Fe的嵌入原子法势和Zhou等［40］开发的嵌入法原子势。

所有模型尺寸均为20a×20a×20a，a为单个晶胞的晶格常数。在3个方向上均使用周期性边界条件。在加载冲击前，首先在NPT系统中对模型在300 K进行40 ps的弛豫，使其达到可以忽略系统中残余应力的状态。弛豫后通过MSST在X轴方向对模型施加冲击波，Ni的波速范围为11～15 km／s，Ni3Nb的波速范围为10～18 km／s，Ni3Ti的波速范围为6～10 km／s，Ni3Al的波速范围为6～10 km／s，Cr的波速范围为11～15 km／s，Fe的波速范围为6～10 km／s，Mo的波速范围为6～10 km／s。当模拟持续足够长的时间，粒子速度、压强、温度和其他物理量均达到稳定后，从稳定状态下的数据中提取需要的数据。

4 结果与讨论

4.1 GH4169镍基高温合金的C0和λ

GH4169镍基高温合金的化学成分复杂，且许多微量金属元素的质量分数小于1%，因此对GH4169的成分进行简化：首先，保留质量分数大于1%的金属元素，满足此条件的元素有Ni、Cr、Fe、Nb、Mo和Ti；其次，由于GH4169镍基合金中具有γ′相，对增强GH4169的高温性能具有重要意义，因此所有γ′相包含的金属元素也要保留，满足此条件的元素有Ni、Nb、Ti和Al。经简化后的GH4169镍基高温合金各主要元素的质量分数如表1所示。该成分来源由上海康晟航材的GH4169高温合金成分表简化而来，对所选元素中Ni以外的元素，基于成分的含量范围进行质量分数取整后，余量即为Ni元素的质量分数。

表1 GH4169镍基合金各元素的质量分数 和原子分数 %

GH4169高温镍基合金中，γ″相为Ni3Nb，是主要强化相；γ′相包括Ni3Al和Ni3Ti，为辅助强化相；Ni是γ相，为基体。所涉及的相中无δ相，根据上海康晟航材给出的热处理制度，无δ相GH4169的热处理方式为（1 010～1 065）℃±10℃，1 h，油冷、空冷或水冷＋720℃±5℃，8 h，以50℃／h炉冷至620℃±5℃，8 h，空冷。以原子占比为依据进行计算，可得各相的理想质量分数如表2所示。此外，通过分子动力学的多尺度模拟技术（MSST）获得的各相的C0和λ如表2中所示。

表2 GH4169镍基合金中各相的质量分数 和Hugoniot参数

表3为通过拟合得到的GH4169镍基高温合金的C0和λ。通过与文献［41］的研究成果进行对比发现，所提出的理论计算得到的C0和λ与通过一级轻气炮进行平板撞击实验得到的C0和λ非常接近，证明所提出的计算理论是正确的。

表3 GH4169镍基高温合金的Hugoniot参数 和相对误差

4.2 Hugoniot曲线及其内能

得到C0和λ后代入式（7）（8）可以得到GH4169合金的Hugoniot曲线及其内能曲线。

图1 GH4169合金的Hugoniot曲线

图2 GH4169合金的内能

红色方块为使用文献［41］方法得到的数据计算结果，可以发现该结果与本文方法计算结果一致，表明本文中所提出方法的计算结果可靠。

4.3 0 K等温线及其内能

通过分子动力学模拟升温过程，将获得的体积数据代入式（11）（12），得到GH4169合金的体积膨胀系数αv为251×10－6K－1。通过式（13）～（15），（22）～（25）可以计算得到Born-Meyer势和Morse势的0 K等温线及其内能曲线。材料常数Q、q、A和B分别为70.55 GPa、11.55、151.52 GPa和4.45。

图3、4为GH4169合金的0 K等温线Pc及其内能Ec。蓝色方块为使用文献［41］方法得到的数据计算结果，黑线和红线分别表示Born-Meyer势和Morse势，这两种势均可以用来描述0 K等温线。

图3 GH4169合金的冷压曲线

图4 GH4169合金的冷能曲线

可以看到，0 K等温线Pc及其内能Ec都随着V0K／V的增大而增大，且V0K／V越大，压力Pc和内能Ec的增长速率越大，其冷曲线特征与Cu和Pb一致［42］。当V0K／V小于1.5时，Born-Meyer势与Morse势的等温线、内能曲线几乎是重合在一起的，且与刘焦等［41］的结果吻合良好；当V0K／V大于1.5时，Born-Meyer势与Morse势的等温线、内能曲线逐渐分离，此时不论是等温线还是内能，Morse势都更接近于文献［41］的结果。这表明当V0K／V小于1.5时，Born-Meyer势和Morse势都是比较准确的，而当V0K／V大于1.5时，Morse势更加准确。

4.4 Grüneisen状态方程

通过式（1）～（2），（6）～（8），（27）和（40）～（41），最终可以得到Grüneisen状态方程的3D曲面图，如图5、6所示。

图5、6分别是以式（40）、（41）作为Grüneisen系数获得的状态方程。虽然从正面看，2个Grüneisen状态方程的形状均为一个凹面，且形状相似，但可以发现最大压力是不相同的。这是由Grüneisen系数方程不同造成的，因此可以推断，Grüneisen系数是影响Grüneisen状态方程压力值的重要因素。

图5 P-V-E形式的Grüneisen状态方程曲面，以式（40）作为Grüneisen系数

图6 P-V-E形式的Grüneisen状态方程曲面，以式（41）作为Grüneisen系数

5 结论

1）使用所提出的计算复杂合金Hugoniot参数的理论对GH4169镍基高温合金的C0和λ进行计算，结果与实验结果吻合度良好，证明该理论可行。

2）GH4169镍基高温合金的Born-Meyer势材料常数Q、q和Morse势材料常数A、B分别为70.55、11.55、151.52 GPa和4.45。当V0K／V小于1.5时，Born-Meyer势和Morse势都可以准确地描述0 K等温线，而当V0K／V大于1.5时，Morse势对于0 K冷曲线的描述更加准确。

3）Grüneisen状态方程在P-V-E空间中表现为一个凹面，Grüneisen系数的模型是影响状态方程最大压力值的重要因素。