Be 与CH 混合材料物理性质分子模拟研究

2022-05-18褚衍运祁建敏马纪敏李正宏

重庆理工大学学报(自然科学) 2022年4期

曹宇，褚衍运，祁建敏，王真，马纪敏，李正宏

(1.中国工程物理研究院核物理与化学研究所，四川绵阳 621900;2.中国核动力研究设计院反应堆工程研究所，成都 610213)

0 引言

铍(Be)元素由于具有简单的原子结构，被广泛应用于核能、物理与航空航天工业。在惯性约束聚变(inertial confinement fusion，ICF)中，铍元素由于具有高密度、低不透明度、高热导率及冲击压缩过程中稳定性较高的优点，被用作间接驱动靶丸的烧蚀层材料。由于靶丸的设计主要基于辐射-流体动力学模拟，其中靶丸材料在烧蚀过程中将经历不同的压力和温度条件，因此，铍金属在极端条件下的物理性质对于优化靶丸设计至关重要。采用实验方法难以精确测量Be 的物性参数，而采用辐射流体力学模拟则需要完整的状态方程(EOS)数据来涵盖相关的等离子体参数，因此，采用可靠的方法计算极端条件下Be 的状态方程具有重要意义。

基于密度泛函理论的分子动力学模拟(density functional theory molecular dynamics，DFT-MD)已被证明是研究温稠密等离子体多体混合系统的可靠方法［1-6］，DFT-MD 模拟可以很好地用于状态方程计算，如氘、CH、硅、碳、铍等元素。Ding 等［7］采用orbital-based Kohn-Sham 分子动力学(Kohn-Sham molecular dynamics，KSMD)和orbital-free 分子动力学(orbital-free molecular dynamics，OFMD)为铍元素建立了密度范围在0.01～500 g/cm3的第一性原理状态方程参数表，其研究结果表明，在靶丸压缩的最大值附近，温度和密度差异最大达到了18%，因此，理论计算出的状态方程对于靶丸烧蚀过程产生了重要影响。由于铍与CH 材料通常作为ICF 靶的烧蚀材料，其状态方程对于辐射流体力学模拟有着一定的影响，并且国内外学者的研究较少深入研究其混合物的物理性质，因此，本文采用第一性原理与反应分子动力学的方法计算Be 与CH 混合物的状态参数。

1 计算方法与模型

1.1 第一性原理分子动力学

采用VASP［8-10］程序计算Be 与CH 混合物材料的物理性质。在广义梯度近似(GGA)下，采用Perdew-Burke-Ernzerhof(PBE)交换关联泛函形式的投影缀加波(PAW)进行处理。其中C、H 与Be的赝势分别采用4 个(2s22p2)、1 个和4 个(1s22s2)价电子，芯半径分别设置为1.1、0.8、1.3 Bohr。K 网格选用单Gamma 点，平面波展开的波函数截断能设置为2 000 eV。模拟采用NVT 系综及Nose-Hoover［11］控温方式，时间步长设置为0.05～0.1 fs。

在温度范围为104～106K 时，采用聚苯乙烯材料(C160H162)与不同比例的Be 金属进行混合计算得到混合物的状态方程。图1 所示为聚苯乙烯与Be 金属混合物在环境温度下达到平衡时的模拟体系。模拟系统中混合物的总密度、Be 金属和聚苯乙烯所占的比例分别定义为 ρmix=。本文采用第一性原理分子动力学方法得到聚苯乙烯与Be 的混合物物理性质与光学性质［12］。

1.2 反应分子动力学

通过采用非平衡反应分子动力学的方法计算Be 与CH 混合物的导热系数，其中，模拟软件采用Materials Studio 8.0［13］进行建模计算，并应用Van等［14-15］提出的ReaxFF 反应力场，可以在高温和高压极端条件下研究金属及碳氢化合物的热物理性质［16-17］。

模拟细节如下:首先，在环境温度下构造和优化随机的聚苯乙烯与Be 金属混合模型如图1 所示。然后，采用NVT 系综对模拟系统进行100 ps的平衡计算，在体系总能量达到最小值并维持稳定时获得优化后的结构模型。最后，通过在NVE系综下进行模拟获得特定的温度梯度，进而计算材料的热导率。在模拟过程中，采用了周期性边界条件，选择Berendsen 控温方式，热浴阻尼系数为0.1 ps，时间步长为0.1 fs，运动方程采用Verlet算法。为了计算体系的平均温度梯度(dT/dz)，将模拟系统z 向尺寸设置为X/Y 向尺寸的3～5 倍，然后沿z 向分割成9 个盒子，进而统计每个小盒子的平均温度，在1 000 K 到10 000 K(间隔为1 000 K)温度范围内采用非平衡反应分子动力学计算混合物的导热系数。

图1 在环境温度压力下平衡时Be 金属与聚苯乙烯混合分子模型示意图

2 结果与讨论

2.1 Be 与CH 混合材料状态方程

采用第一性原理分子动力学方法计算得到了聚苯乙烯与不同比例的Be 金属混合物状态参数。图2(a)所示为混合体系在密度为1 g/cm3时含有0～0.5 Be 金属的温度-压力曲线。从图2(a)中可知，在掺杂少量的Be 金属时，聚苯乙烯体系的压力发生了明显的降低，并且随着Be 含量的增加，压力不断减小，在106K 时，压力最大降低了13%。图2(b)所示为含有0.41 Be 金属的混合体系在不同密度下的温度-压力曲线。从图2(b)中可知，在混合物体系密度为1、2 g/cm3时，温度-压力关系近似呈线性分布，在密度为4、5 g/cm3且温度低于105K 时，存在一定的曲率。曲率的形成是由于聚合物发生热分解导致模拟体系分子种类的变化，从而使状态方程发生改变，因此，聚合物热解对于高密度的混合体系有较大影响，而对低密度的CH 材料影响较小。

图2 聚苯乙烯与Be 金属混合体系的温度-压力曲线

2.2 Be 与CH 混合材料导热系数

图3(a)所示为初始温度分别设置为1 000、3 000、5 000 K 时模拟体系的温度分布情况，由于将初始体系沿z 向分割成9 个模拟盒，因此，分别统计了每个模拟盒内的平均温度。从图3(a)中可以得出，模拟盒两侧的温度高于中间点的温度，因此，存在热量从两端向中间传递，可以通过计算得到模拟体系的平均温度梯度，图3(b)所示为混合模拟体系在不同温度达到热平衡时，平均温度梯度随模拟时间的变化情况。从图3(b)中可以得出，在模拟的初始阶段(t＜10 ps)，模拟体系的平均温度梯度存在较大的波动，但随着模拟时间的延长，平均温度梯度逐渐趋于稳定值，因此，在50 ps 的模拟计算中温度梯度达到了收敛。此外，平均温度梯度随着初始温度的升高而增加，说明初始温度设置越高，模拟盒内温度波动越大。

图3 聚苯乙烯与Be 金属混合体系在不同的初始条件下的温度及温度梯度分布曲线

图4 所示为混合体系在密度为3 g/cm3时体系的平均热流密度随时间的变化情况，从图4 中可以得出，随着温度升高，模拟盒内温度波动也不断增加，从而导致热流密度不断增加，在模拟的初始阶段(t＜10 ps)，由于模拟体系未达到热平衡，导致热流存在一定的波动，但是在模拟后期，热流密度的数值逐渐趋于稳定，因此，在50 ps 非平衡反应分子动力学模拟中计算出的热流数据是收敛且可靠的。

图4 聚苯乙烯与Be 金属混合体系在不同的初始温度下平均热流随时间的变化曲线

图5(a)所示为采用反应分子动力学方法计算得到的聚苯乙烯与0.5 Be 金属混合体系热导率随温度的变化情况，从图中可以得出:聚苯乙烯与Be 金属混合物热导率随着温度升高呈现波动上升的趋势，在9 000 K时，热导率最高达到了98 W/m·K。在较低的温度条件下(＜5 000 K)，聚苯乙烯与Be 金属混合体系热导率随着密度升高而增加，但是高温条件下聚合物会发生一定程度的热分解，从而导致体系分子结构的改变，因此，在较高温度下呈现热导率先减少后增加的趋势。

图5 聚苯乙烯与Be 金属混合体系的热导率随温度变化的情况

图5(b)所示为聚苯乙烯在含有不同比例的Be 金属构成混合物体系时热导率随温度的变化情况，从图5(b)中可以得出:纯聚苯乙烯体系热导率随温度上升不断增加，并在9 000 K 时增加到4.2 W/m·K。在添加0.5 Be 金属后，由于体系原子结构发生了改变，导致热导率降低，但是混合体系的热导率随温度增加仍然波动上升，在7 000 K时，达到最大值2.6 W/m·K。纯Be 金属体系的热导率相比混合体系显著降低，在9 000 K 的温度下，热导率降低了近40%。总体而言，纯聚苯乙烯材料在掺杂Be 金属后，热导率存在不同程度的降低，但是总体随温度呈现上升的趋势。

2.3 Be 与CH 混合材料光学性质

为了研究材料的光学性质，选择折射率N 描述电磁波通过真空和其他材料传播时的差异:

式中N 在真空里为实数，数值等于1。N 的虚部k与吸收系数有关。

式中c 代表光速。

材料的反射率定义为:

式中n0代表环境的折射率。在计算光学性质时，通常需要计算介电函数，然后推导出其他性质。介电函数ε 的定义如下:

可以获得介电函数与折射率的实部和虚部之间的关系为:

因此，可以根据材料的介电函数计算折射率，得到吸收系数与反射率等光学性质参数，进而研究Be 金属对聚苯乙烯光学性质的影响，类似的计算方法已广泛应用于建立CH［18-19］及氘［20］材料的光学参数。基于密度泛函理论计算出介电函数虚部的数值ε2，然后通过Kramas-Kronig 关系式同时计算出实部ε1［21］。介电函数部分中，ε2出现的峰值部分与第一布里渊区能带间电子跃迁直接相关。图6(a)所示为温度在10 000 K 时的聚苯乙烯和含0.5 Be 金属的混合物体系在不同密度下的介电函数。从图6(a)中可以得出，混合物密度越高，介电函数虚部数值越大，并且实部的峰值位置也发生了改变，这说明在不同密度条件下，聚合物热分解产物是有区别的，因此导致体系的分子结构发生变化，体现在介电函数存在一定的差异。图6(b)所示为纯聚苯乙烯与含0.5 Be 金属的混合物体系在温度与密度分别为3 000 K 与1 g/cm3时的折射率，其中折射率的实部与虚部分别定义为n 与k。从图6(b)中可以得出:由于Be 与CH材料具有不同的原子结构，因此添加0.5 Be 金属后聚苯乙烯材料的折射率发生了改变，在光子能量处在较低范围时(hv＜10 eV)，添加Be 金属对聚苯乙烯光学性质影响较大，但是随着光子能量不断增加，混合体系与纯聚苯乙烯体系折射率趋于一致。

图6 聚苯乙烯与Be 金属混合体系在不同密度下的介电函数与折射率曲线

根据式(3)中材料折射率与反射率之间的关系，可以得到聚苯乙烯与Be 金属混合体系的反射率。图7(a)所示为含0.5 Be 金属的混合体系在不同密度下的反射率，从图7(a)中可以得出:在光子能量范围为0～5 eV 区间内，密度对混合物体系的反射率影响较小，随着光子能量增加，计算出的反射率在不同密度下差别逐渐变大，并且呈现密度越高、反射率越大的规律。总体而言，由于聚苯乙烯在不同密度下的热解程度不同，将导致体系生成不同的分子碎片，从而影响体系的反射率。为了研究在相同的密度与混合比的条件下，不同温度对混合物体系的影响，计算得到了在1 000、5 000、10 000 K 温度下，含0.5 Be 金属混合物材料的反射率，如图7(b)所示。从图7(b)中可以得出:温度对混合体系反射率存在一定的影响，其反射率随着温度升高而降低，在温度从1 000 K 变化到5 000 K 时，反射率最高减少了17%，但是随着温度继续升高，反射率的变化较小，在10 000 K时，反射率较5 000 K 仅减少了5%，这可能是由于聚苯乙烯材料在5 000 K 以上的温度下基本热解成原子体系，此时在不同温度下混合体系原子结构相似，因此其光学性质也非常接近。

图7 聚苯乙烯与Be 金属混合体系的反射率变化情况

为了探究含Be 金属对聚苯乙烯反射率的影响，图8 给出了纯聚苯乙烯与含0.5 Be 金属的混合体系在3 000、10 000 K 温度下的反射率。从图8 中可以得出:在较低的光子能量区间段(0～5 eV)，Be 金属对聚苯乙烯反射率影响较小，但是随着光子能量不断增加，其影响逐渐变大。温度对纯聚苯乙烯的反射率影响较大，随着温度增加，其反射率从0.75 减少到0.15 左右。含有0.5 Be金属时，混合体系反射率随着温度升高最多减少了60%。

图8 纯聚苯乙烯与含0.5 Be 金属混合体系的反射率在不同温度下的变化曲线

根据式(2)中材料折射率与吸收系数之间的关系，可以得到聚苯乙烯与Be 金属混合体系的吸收系数。图9(a)所示为含0.5 Be 金属的混合体系在密度为1～5 g/cm3时的反射率，从图9(a)中可以得出:在光子能量范围为0～1 eV 区间内，此时密度对混合物体系的吸收系数影响较小，随着光子能量增加，计算出的吸收系数在不同密度下差别逐渐变大，并且随着密度增加而增大。总体而言，由于聚苯乙烯在不同密度下的热解程度不同，将导致体系生成不同种类的分子与离子，体系结构的不同会改变混合物的吸收系数。为了探究体系在相同的密度与混合比的条件下，不同温度对混合物材料的影响，计算得到了在1 000、5 000、10 000 K 温度下，含0.5 Be 金属混合物材料的吸收系数，如图9(b)所示。从图9(b)中可以得出:温度对混合体系吸收系数存在一定的影响，其大小随着温度升高而降低，在温度从1 000 K 变化到5 000 K 时，吸收系数最高减少了4%，但是随着温度继续升高，其吸收系数在不同温度下变化较小，在10 000 K 时，混合物体系的吸收系数与纯聚苯乙烯体系非常接近，误差最大为2.3%，这可能是由于在高温下聚合物完全分解形成了相似的原子结构，进而导致吸收系数高度一致。

图9 聚苯乙烯与含0.5 Be 金属混合体系的吸收系数变化曲线

为了探究含Be 金属对聚苯乙烯吸收系数的影响，图10 给出了纯聚苯乙烯与含0.5 Be 金属的混合体系在3 000、10 000 K 温度下的吸收系数。从图10 中可以得出:在较低的光子能量区间段(0～1 eV)，Be 金属对聚苯乙烯吸收系数影响较小，但是随着光子能量不断增加，其影响逐渐变大，在光子能量为2.5 eV 时差别最高达13%。此外，由于Be 与CH 材料原子结构存在差异，因此吸收系数的峰值出现的位置也发生了变化，含0.5 Be金属混合物体系吸收系数的峰值位置较纯聚苯乙烯材料左移了0.05 eV。