基于Tersoff势的晶格中波动传播

周子清; 王鹏飞; 徐松林

doi:10.11883/bzycj-2024-0007

基于Tersoff势的晶格中波动传播

doi: 10.11883/bzycj-2024-0007

中国科学技术大学中国科学院材料力学行为和设计重点实验室，安徽合肥 230027

基金项目: 国家自然科学基金（12372372，11672286，11872361）；高压物理与地震科技联合实验室开放基金（2019HPPES01）；中石油与中科院重大战略合作项目（2015A-4812）；中央高校基本科研业务费专项资金（WK2480000008）

详细信息

作者简介:
周子清（1997－　），男，硕士研究生，zhzqing@mail.ustc.edu.cn

通讯作者:
王鹏飞（1985－　），男，博士，副研究员，pfwang5@ustc.edu.cn

中图分类号: O347.4
计量
- 文章访问数: 257
- HTML全文浏览量: 96
- PDF下载量: 91
- 被引次数: 0
出版历程
- 收稿日期: 2024-01-02
- 修回日期: 2024-03-13
- 网络出版日期: 2024-04-23
- 刊出日期: 2024-09-20

Wave propagation in lattices based on Tersoff potential

CAS Key Laboratory of Mechanical Behavior and Design of Materials, University of Scienceand Technology of China, Hefei 230027, Anhui, China

摘要

摘要: 在晶格间的Tersoff势作用下分别研究了单晶体系和多晶体系中的波动传播特性。首先，在微振动的情况下，分别基于晶格间线性作用、Tersoff势作用以及含缺陷的Tersoff势作用3种势能函数研究了单晶体系中格波的传播，得到了晶格中的色散关系以及格波波速的表达式。其次，分别以碳晶格和硅晶格为例，运用有限差分方法，研究了3种势能作用下单晶体系中的波动传播过程，对比了压缩和拉伸冲击下晶格的运动差异，并讨论了入射速度对位移峰值和受力峰值的影响，揭示了单晶体系中波动传播与连续介质中波动传播的差异。最后，分别以金刚石和碳化硅为例，采用分子动力学模拟方法，研究了多晶体系中的波动传播特性，讨论了不同空间位置原子的运动差异。结果表明：多晶体系中晶格结构更复杂，其中的波动传播特性与单晶体系存在差异；缺陷的存在对波动传播规律影响显著，这种影响在多晶体系中表现得更加突出。
- 晶格动力学 /
- 微振动 /
- Tersoff势 /
- 缺陷 /
- 分子动力学
Abstract: The propagation characteristics of waves are the basis for studying the dynamic behavior of materials, and the theoretical study of waves in continuous media at the macro scale has been well developed. With the widespread application of materials and structures at the micro- and nano- scales, the study of wave propagation characteristics at the lattice scale is receiving increasing attention. In this article, the Tersoff potential interaction between lattices is applied to study the wave propagation characteristics in single-crystal and polycrystalline systems. Firstly, in the case of micro-vibration, the propagation of lattice waves in a single-crystal system is studied based on three potential energy functions between lattices: linear interaction, Tersoff potential, and Tersoff potential with defects. The dispersion relationship in the lattice and the expression of lattice wave velocity are obtained. Secondly, taking carbon lattice and silicon lattice as examples, the finite difference method is applied to study the wave propagation process in the single-crystal system under three potential energies. The differences in lattice motion under compressive and tensile impacts are compared, and the influence of incident velocity on the peak displacement and peak force is discussed, which reveals the difference in wave propagation between single-crystal systems and continuous media. Finally, taking diamond and silicon carbide as examples, molecular dynamics simulations are used to study the wave propagation characteristics in polycrystalline systems, and the differences in atomic motion at different spatial positions are discussed. The results indicate that the lattice structure in polycrystalline systems is more complex, and the wave propagation characteristics in polycrystalline systems are different from those in single-crystal systems. The existence of defects has a significant impact on the propagation law of waves, which is more prominent in polycrystalline systems. This study has good reference significance for the study of material dynamics performance at the micro- and nano- scales.
- lattice dynamics /
- micro-vibration /
- Tersoff potential /
- defect /
- molecular dynamics

HTML全文

波动传播特性是研究材料动力学行为的基础。连续介质中的波动传播研究主要基于3类方程：几何关系或连续性方程、运动学方程和本构方程。较高的冲击压力会造成局部较高的温升，本构方程会替换成状态方程^[1]。此理论体系对于宏观介质的研究已经非常完善，取得了丰硕的研究成果^[2-5]。随着微纳尺度材料和结构的广泛应用^[6-9]，针对微纳尺度材料的动力学特性，尤其是波动传播规律的研究显得非常重要。微纳尺度下，连续介质假设失效，如何通过方程描述波动传播至关重要。

宏观尺度的颗粒介质考虑的是颗粒接触力，例如赫兹接触等^[10-11]，其运动学方程主要基于接触作用建立，由此进行载荷和变形的传递^[12-13]。Tang等^[14]研究了颗粒体系中颗粒形态对弹性波传播的影响。Alberdi等^[15]采用松弛微形态模型研究了由不同晶胞组成的异质元结构中的波传播。Waymel等^[16]研究了弹塑性波在由接触金属颗粒组成的二维有序颗粒介质中的传播，发现塑性变形导致振幅衰减。上述研究较好地描述了材料的宏观动力学行为。

与连续介质理论不同，微纳尺度的晶格间作用力可以为原子或分子提供稳固的结构，同时也可以进行原子或分子间的信息传递^[17-19]。晶格间以各种作用势定义的作用力发生相互作用，其运动学方程需要基于势作用力建立。Hu等^[20]研究了碳纳米管中的横波和扭转波，讨论了碳纳米管微观结构对波色散的影响。Agarwal等^[21]通过分子动力学方法研究了纳米晶Al微结构在原子尺度冲击压缩过程中的变形响应。Sam等^[22]采用动态结构因子和应变波动等方法评估了全硅沸石中的声学特性，讨论了声波在纳米多孔材料中的传播。但上述研究大多从分子动力学出发，尚未从理论上揭示清楚微纳尺度下晶格运动的基本规律以及波动传播的特性。

Tersoff势能够描述多个原子间共价键的作用，可通过考虑多体效应来实现共价键的形成和离解，其参数取决于所讨论的键的局部环境。这种势与薛定谔基于量子态波函数描述的粒子波动不同^[23-26]，与传统的分子间力场也不相同。本文中，拟从晶格尺度的运动方程出发对波动传播特性进行研究。首先，对线性作用、Tersoff势作用以及含缺陷的Tersoff势作用下单晶体系微振动中格波的传播进行研究，阐明格波与连续波之间的联系；其次，采用有限差分方法，分别以碳原子和硅原子为例，研究单晶体系中的波动传播特性；最后，采用分子动力学方法，以金刚石（C₄）和碳化硅（SiC）为例，研究多晶体系中的波动传播特性，以期探索晶格体系中的波传播规律。

1. 单晶中的微振动与格波

1.1 线性作用

晶体的微振动可按谐振近似处理。晶格的微振动模型如图1^[27]所示。考虑图1中由N个相同原子组成的简单立方晶体的一维纵向晶格的微振动。对于一维纵向运动，简单立方晶体可简化为一维原子链的振动。由于晶格的微振动主要受相邻晶格的影响，非相邻微粒的影响可以忽略不计，因此，每个原子只考虑受到相邻原子的作用。设晶格中的原子质量为m，平衡时的原子间距为a，以X_n表示第n个原子由于微振动而偏离平衡位置的位移，则第n个与第n+1个原子间的相对位移为：

图 1 一维原子链振动模型^[27]

Figure 1. Vibration model of one-dimensional atomic chain^[27]

下载: 全尺寸图片幻灯片

$\Delta X = {X_n} - {X_{n + 1}}$

(1)

通过假定原子的受力模型，可以分析晶格微振动过程中的运动规律。

对图1所示的微振动模型，原子间作用力考虑为线性作用，则第n个原子受到第n+1个原子的作用力^[28]可以表示为：

${F_{n,n + 1}} = - \alpha \left( {{X_n} - {X_{n + 1}}} \right)$

(2)

式中： $\alpha$ 为线性系数。定义作用力与原子间距的比值 $\varphi$ 为：

$\varphi = \frac{{{F_{n,n + 1}}}}{{{X_{n + 1}} - {X_n}}}$

(3)

当原子间假定为线性作用时，可以得到 $\varphi = \alpha$ ，即原子间作用力与原子间距的比值为常数。对于第n个原子，根据牛顿运动定律，其运动方程为：

$m{\ddot{X}}_{n}=\alpha \left({X}_{n+1}-2{X}_{n}+{X}_{n-1}\right)$

(4)

式(4)对于所有原子均成立，其解具有简谐波形式：

${X_n} = A{{\text{e}}^{{\text{i}}\left( {\omega t - 2{\text{π}}nak} \right)}}$

(5)

式中：A为振幅；ω为角频率；a为晶格常数，即原子间平衡间距；na为第n个原子相对于原点的平衡位置；k=1/λ为波数，λ为波长。

引入以下关系：

$\left\{ {\begin{array}{*{20}{l}} {\rho = {m}/{{{a^3}}}} \\ {E = {\alpha }/{a}} \\ {{c_0} = \sqrt {{E}/{\rho }} = a\sqrt {{\alpha }/{m}} } \end{array}} \right.$

(6)

式中：ρ为晶体密度，E为杨氏模量，c₀为一维应力弹性波速。

将式(5)代入式(4)，得到：

$\omega {\text{ = 2}}\sqrt {\frac{\alpha }{m}} \sin \left( {{\text{π}}ak} \right) = {\omega _0}\sin \left( {{\text{π}}ak} \right)$

(7)

${\omega _0} = \frac{\omega }{{\sin \left( {{\text{π}}ak} \right)}} = 2\sqrt {\frac{\alpha }{m}}$

(8)

式中：ω₀为角频率峰值。式(5)表明，各原子在平衡位置微振动时，以简谐波形式在晶格中传播，称为格波。式(7)描述了晶格间波传播依赖于频率的弥散现象，也称为色散关系。

当 $ak \ll 1$ 时，即在长波情况下， $\sin \left( {{\text{π}}ak} \right) \approx {\text{π}}ak$ ，波速的表达式为：

$c=\frac{1}{k}\frac{\omega }{2\text{π}}=a\sqrt{\frac{\alpha }{m}}\frac{\mathrm{sin}\left(\text{π}ak\right)}{\text{π}ak}\approx a\sqrt{\frac{\alpha }{m}}=\sqrt{\frac{E}{\rho }}={c}_{0}$

(9)

式(9)表明格波是一个波速依赖于晶格特征尺寸（晶格常数a）与波长λ（λ=1/k）之比的弥散波。钢的晶格常数a约为3 Å，平均振动时间约为10^–13 s，其格波波速约为3 km/s，接近其弹性波波速（3.5 km/s）^[2]。同时，式(7)表明：长波情况下格波的波速与波数无关，与材料特性相关。当波长远大于晶格常数时，晶体可以看作连续介质。

1.2 Tersoff势作用

Tersoff多体势函数适用于描述共价键结合的多个原子之间的相互作用，可以很好地表示三维空间分子表面的重构能，能够准确地描述C与Si晶格中的作用。

线性作用不能很好地描述原子间作用：在原子间距减小时，原子间斥力急剧增大，而原子间距增大时，原子间引力缓慢增大。本文中，以C和Si晶格为研究对象，并通过Tersoff势来描述原子间的相互作用。Tersoff势的表达式^[29]为：

$\left\{ \begin{array}{*{20}{l}} E = \dfrac{1}{2}\displaystyle\sum\limits_i {\displaystyle\sum\limits_{j \ne i} {{V_{ij}}} } \\ {V_{ij}} = {f_{\rm{C}}}( {{r_{ij}}})\left[ {{f_{\rm{R}}}( {{r_{ij}}} ) + {b_{ij}}{f_{\rm{A}}}( {{r_{ij}}} )} \right] \\ {f_{\rm{C}}}\left( r \right) = \left\{ \begin{array}{*{20}{l}} 1 & \;\;\;{r \leqslant R - D} \\ {\dfrac{1}{2} - \dfrac{1}{2}\sin \left( {\dfrac{{\text{π}} }{2}\dfrac{{r - R}}{D}} \right)}&\;\;\; {R - D {\text{＜}} r {\text{＜}} R + D}\\ 0 & \;\;\;{r \geqslant R + D} \end{array} \right. \\ {{f_{\rm{R}}}\left( r \right) = {A_1}\exp \left( { - {\lambda _1}r} \right)} \\ {{f_{\rm{A}}}\left( r \right) = - B\exp \left( { - {\lambda _2}r} \right)} \\ {{b_{ij}} = {{\left( {1 + {\beta ^{{n_1}}}\zeta _{ij}^{{n_1}}} \right)}^{{-} \tfrac{1}{{2{n_1}}}}}} \\ {{\zeta _{ij}} = \displaystyle\sum\limits_{k \ne i,j} {{f_{\rm{C}}}( {{r_{ij}}} )g( {{\theta _{ijk}}} )\exp \left[ {\lambda _3^{{m_1}}{{( {{r_{ij}} - {r_{ik}}} )}^{{m_1}}}} \right]} } \\ {g\left( \theta \right) = {\gamma _{ijk}}\left[ {1 + \dfrac{{{C^2}}}{{{d^2}}} - \dfrac{{{C^2}}}{{ {{d^2} + {{\left( {\cos \theta - \cos {\theta _0}} \right)}^2}} }}} \right]} \end{array} \right.$

(10)

式中：m₁、γ、λ₃、C、d、cosθ₀、n₁、β、λ₂、B、R、D、λ₁、A₁为材料参数。Tersoff势包含两部分：斥力项和引力项，对于不同的共价键，参数的取值不同，具体参数如表1^[30]所示。

表 1 Tersoff势函数中的材料参数^[30]

Table 1. Material parameters in Tersoff potential^[30]

共价键种类	m₁	γ	λ₃/Å^–1	C	d	cosθ₀	n₁
C―C	3.0	1.0	0	38 049	4.348 4	–0.570 58	0.727 51
Si―Si	3.0	1.0	0	100 390	16.217 0	–0.598 25	0.787 34

共价键种类	β	λ₂/Å^–1	B/eV	R/Å	D/Å	λ₁/Å^–1	A₁/eV
C―C	1.572 4×10^–7	2.211 90	346.70	1.95	0.15	3.487 9	1 393.6
Si―Si	1.100 0×10^–6	1.732 22	471.18	2.85	0.15	2.479 9	1 830.8

下载: 导出CSV

| 显示表格

表 2 2种晶格在平衡位置处的参数

Table 2. Parameters of two types of lattices at equilibrium positions

原子种类	摩尔质量/(g·mol⁻¹)	平衡距离/Å	晶格常数/Å
C	12	1.54	3.57
Si	28	2.35	5.43

下载: 导出CSV

| 显示表格

考虑图1的一维原子链模型，取θ = π，则b_ij转化为常数b₀。在平衡位置附近， ${f_{\rm{C}}}\left( r \right) = 1$ ，Tersoff势作用下每个原子的势能和受力的表达形式可简化为：

$u( r ) = {A_1}\exp \left( { - {\lambda _1}r} \right) - {b_{\text{0}}}B\exp \left( { - {\lambda _2}r} \right)$

(11)

$f\left( r \right) = - \frac{{{\rm{d}}u( r )}}{{{\rm{d}}r}} = {\lambda _1}{A_1}\exp \left( { - {\lambda _1}r} \right) - {b_{\text{0}}}{\lambda _2}B\exp \left( { - {\lambda _2}r} \right)$

(12)

对于图1的微振动模型，第n个原子受到第n+1个原子的作用力为：

${F_{n,n + 1}} = f\left( {{X_0}} \right) - f\left( {{X_{n + 1}} - {X_n} + {X_0}} \right)$

(13)

此时，第n个原子微振动的控制方程为：

$m{\ddot{X}}_{n}=-f\left({X}_{n+1}-{X}_{n}+{X}_{0}\right)+f\left({X}_{n}-{X}_{n-1}+{X}_{0}\right)$

(14)

采用简谐解 ${X_n} = A{{\text{e}}^{{\text{i}}\left( {\omega t - 2{\text{π}}nak} \right)}}$ ，由于微振动过程 $X{}_{n + 1} - {X_n} \ll {X_0} = a$ ，因此，可以得到：

$\omega \approx 2\sqrt {\frac{1}{m}\left( {\lambda _1^2{A_1}{{\text{e}}^{ - {\lambda _1}a}} - {b_0}\lambda _2^2B{{\text{e}}^{ - {\lambda _2}a}}} \right)\left[ {1 + 2{\lambda _1}{{\sin }^2}\left( {{\text{π}}ak} \right){X_n}} \right]} \sin \left( {{\text{π}}ak} \right)$

(15)

式(15)为Tersoff势作用下的ω-k关系。可以看出，在考虑Tersoff势的情况下，角频率峰值受X_n影响。当 ${X_n} = 0$ 时，原子处于平衡位置，此时角频率峰值为ω₀ ，其表达式为：

${\omega _0} = 2\sqrt {\frac{1}{m}\left( {\lambda _1^2{A_1}{{\text{e}}^{ - {\lambda _1}a}} - {b_0}\lambda _2^2B{{\text{e}}^{ - {\lambda _2}a}}} \right)}$

(16)

考虑微振动过程，由式(15)可以近似得到：

$\varphi = \left( {\lambda _1^2{A_1}{{\text{e}}^{ - {\lambda _1}a}} - {b_0}\lambda _2^2B{{\text{e}}^{ - {\lambda _2}a}}} \right)\left[ {1 + 2{\lambda _1}{{\sin }^2}\left( {{\text{π}}ak} \right){X_n}} \right]$

(17)

式(17)描述了原子间作用力与原子间距的比值 $\varphi$ 随X_n的变化。特别地，在X_n=0时，Tersoff势作用转化为线性作用，此时 $\varphi$ 为Tersoff势在平衡位置处的一阶导数。

图2给出了原子间作用为Tersoff势时的角频率ω与波数k之间的对应关系。可以看出：当X_n=0时，ω与k的关系与线性作用时相同；当X_n>0，即原子间距离减少时，ω随X_n增大而增大；而当X_n<0，即原子间距离增大时，ω随X_n减小而减小。

图 2 Tersoff势下角频率与波数的关系

Figure 2. Relationship between angular frequency and wavenumber under Tersoff potential

下载: 全尺寸图片幻灯片

当 $ak \ll 1$ 时，即在长波情况下， $\sin \left( {{\text{π}}ak} \right) \approx {\text{π}}ak$ ，对应的波速表达式为：

$c=\frac{1}{k}\frac{\omega }{2\text{π}}\approx a\sqrt{\frac{1}{m}\left({\lambda }_{1}^{2}{A}_{1}{\text{e}}^{-{\lambda }_{1}a}-{b}_{0}{\lambda }_{2}^{2}B{\text{e}}^{-{\lambda }_{2}a}\right)\left(1+2{\lambda }_{1}{\text{π}}^{\text{2}}{a}^{\text{2}}{k}^{\text{2}}{X}_{n}\right)}={c}_{0}\sqrt{1+2{\lambda }_{1}{\text{π}}^{\text{2}}{a}^{\text{2}}{k}^{\text{2}}{X}_{n}}$

(18)

可以看出，格波波速随着X_n变化而变化。

1.3 含缺陷的Tersoff势作用

晶格中存在位错等缺陷时，由于局部晶粒的缺失，粒子间的相互作用发生变化。无位错时，晶格考虑滑移过程，在平衡位置时势能最小，结合力为零；滑移半个原子间距时，需要克服的势能最大，此时对应的势能也最大；滑移一个原子间距后，恢复到初始状态，势能降到最小。此过程中，势能可取为周期性函数^[28]：

$u( r ) = - 2{A_2}\cos {\frac{{2{\text{π}}r}}{a}}$

(19)

式中：A₂为势能峰值。

若存在一个刃型位错，即在位错滑移面上方有N+1个原子，其下方只有N个原子，则滑移过程中晶格势能和作用力可表示为：

$u( r ) = - {A_2}\cos \left( {\dfrac{{2{\text{π}}r}}{a}\dfrac{{N + 1}}{{N + {1}/{2}}}} \right) - {A_2}\cos \left( {\dfrac{{2{\text{π}}r}}{a}\dfrac{N}{{N + {1}/{2}}}} \right)$

(20)

$f\left( r \right) = - \dfrac{{{\rm{d}}u( r )}}{{{\rm{d}}r}} = - \dfrac{{2{\text{π}}{A_2}}}{a}\dfrac{{N + 1}}{{N + {1}/{2}}}\sin \left( {\dfrac{{2{\text{π}}r}}{a}\dfrac{{N + 1}}{{N + {1}/{2}}}} \right) - \dfrac{{2{\text{π}}{A_2}}}{a}\dfrac{N}{{N + {1}/{2}}}\sin \left( {\dfrac{{2{\text{π}}r}}{a}\dfrac{N}{{N + {1}/{2}}}} \right)$

(21)

式(19)～(21)表明：由于位错的存在，晶体振动过程中原子间作用力减小。此时，无位错和存在位错的作用力f (r)的峰值分别为 $\dfrac{{4{\text{π}}{A_2}}}{a}$ 和 $\dfrac{{4{\text{π}}{A_2}}}{a}\cos \dfrac{{\text{π}}}{{4N + 2}}$ 。

微振动模型中存在的位错可近似等效为原子间作用力减小，若原子间作用为Tersoff势，第p个原子处受到位错影响，则第p个原子的运动方程为：

$m{\ddot{X}}_{p}=\eta \left[-f\left({X}_{p+1}-{X}_{p}+{X}_{0}\right)+f\left({X}_{p}-{X}_{p-1}+{X}_{0}\right)\right]\begin{array}{cc}\end{array}$

(22)

式中：η为损失因子，代表位错引起的作用力衰减， $0 {\text{＜}} \eta {\text{＜}} 1$ 。

依据式(19)～(21)，在位错处原子间作用受损失因子η 影响，则：

$\eta = \cos \frac{{\text{π}}}{{4N + 2}}$

(23)

采用简谐解 ${X_n} = A{{\text{e}}^{{\text{i}}\left( {\omega t - 2{\text{π}}nak} \right)}}$ ，可以得到以下关系：

$\omega \approx 2\sqrt {\frac{\eta }{m}\left( {\lambda _1^2{A_1}{{\text{e}}^{ - {\lambda _1}a}} - {b_0}\lambda _2^2B{{\text{e}}^{ - {\lambda _2}a}}} \right)\left[ {1 + 2{\lambda _1}{{\sin }^2}\left( {{\text{π}}ak} \right){X_n}} \right]} \sin \left( {{\text{π}}ak} \right)$

(24)

${\omega _0} = 2\sqrt {\frac{\eta }{m}\left( {\lambda _1^2{A_1}{{\text{e}}^{ - {\lambda _1}a}} - {b_0}\lambda _2^2B{{\text{e}}^{ - {\lambda _2}a}}} \right)}$

(25)

$\varphi = \eta \left( {\lambda _1^2{A_1}{{\text{e}}^{ - {\lambda _1}a}} - {b_0}\lambda _2^2B{{\text{e}}^{ - {\lambda _2}a}}} \right)\left[ {1 + 2{\lambda _1}{{\sin }^2}\left( {{\text{π}}ak} \right){X_n}} \right]$

(26)

当 $ak \ll 1$ 时，即在长波情况下， $\sin \left( {{\text{π}}ak} \right) \approx {\text{π}}ak$ ，此时波速的表达式为：

$c=\frac{1}{k} \frac{\omega }{2\text{π}}\approx {c}_{0}\sqrt{\eta \left(1+2{\lambda }_{1}{\text{π}}^{\text{2}}{a}^{\text{2}}{k}^{\text{2}}{X}_{n}\right)}$

(27)

式(21)～(25)表明：存在位错的振动模型，一方面与Tersoff势作用呈现相似的规律，ω、 $\varphi$ 和c随着X_n变化而变化，另一方面，位错造成了波的减弱，导致ω、 $\varphi$ 和c减小。图3给出了存在位错缺陷的情况下，作用为Tersoff势时的角频率ω与波数k之间的对应关系。ω的峰值随X_n变化而改变，对比无位错情况，角频率的峰值及其上下极限均呈现减小的趋势。

图 3 含缺陷的Tersoff势下角频率与波数的关系

Figure 3. Relationship between angular frequency and wavenumber under Tersoff potential with defects

下载: 全尺寸图片幻灯片

2. 单晶中波传播分析

2.1 波传播模型及差分方法

晶格处于势力场中，受到其他晶格的作用。其中，非相邻晶格间作用远小于相邻原子间作用。为了探究单晶体系中波的传播过程，采用有限差分方法对运动过程进行分析，即对运动方程进行时间和空间差分。以C和Si晶格为例，分析单晶下的波传播过程，其参数如表1~2所示。晶格波动计算中，初始晶格均处于平衡位置。

一般地，第n个晶格的运动方程为：

$m{\ddot{X}}_{n}={F}_{n,n+1}-{F}_{n-1,n}$

(28)

对式(28)的时间项进行中心差分，其具有二阶精度：

$\frac{m}{{\Delta {t^2}}}\left( {X_n^{t + 1} - 2X_n^t + X_n^{t - 1}} \right) = {F_{n,n + 1}} - {F_{n - 1,n}}$

(29)

当晶格间作用为线性作用时，可以得到：

$\frac{m}{{\Delta {t^2}}}\left( {X_n^{t + 1} - 2X_n^t + X_n^{t - 1}} \right) = \alpha \left( {X_{n + 1}^t - 2X_n^t + X_{n - 1}^t} \right)$

(30)

当晶格间作用为Tersoff势时，可以得到：

$\frac{m}{{\Delta {t^2}}}\left( {X_n^{t + 1} - 2X_n^t + X_n^{t - 1}} \right) = \left. {\frac{{{\rm{d}}u( r )}}{{{\rm{d}}r}}} \right|_{r = X_{n + 1}^t - X_n^t + a}^{} - \left. {\frac{{{\rm{d}}u( r )}}{{{\rm{d}}r}}} \right|_{r = X_n^t - X_{n - 1}^t + a}^{}$

(31)

当晶格间存在缺陷，且作用为Tersoff势时，在位错处可以表达为：

$\frac{m}{{\Delta {t^2}}}\left( {X_n^{t + 1} - 2X_n^t + X_n^{t - 1}} \right) = \eta \left. {\frac{{{\rm{d}}u( r )}}{{{\rm{d}}r}}} \right|_{r = X_{n + 1}^t - X_n^t + a}^{} - \left. {\eta \frac{{{\rm{d}}u( r )}}{{{\rm{d}}r}}} \right|_{r = X{}_n^t - X{}_{n-1}^t + a}^{}$

(32)

图4为3种势作用下碳晶格中的相互作用与晶格间距的关系。其中，在无缺陷晶格中，平衡位置处线性作用与Tersoff势作用相同；晶格间距减小时，Tersoff势下晶格间斥力急速增大；晶格间距增大时，Tersoff势下晶格间引力缓慢增大。缺陷的存在会导致晶格间作用减弱。

图 4 碳晶格中势力与晶格间距的关系

Figure 4. Relationship between force and lattice spacing in carbon lattices

下载: 全尺寸图片幻灯片

通过有限差分方法，可以得到不同情况下的运动差分方程，进而可以得到各个晶格的运动过程，研究波的传播过程。时间步长 $\Delta t$ 应满足 $\Delta t \ll {a}/{c_0}$ ，其中 ${a}/{c_0}$ 为波传过一个晶格间平衡间距的时间。

2.2 碳晶格中的波传播特性

在波的传播过程中，晶格间分别选取线性作用、Tersoff势作用以及含缺陷的Tersoff势作用。其中，将第3个与第4个晶格间的初始间距调整为1.1a来实现缺陷。开始时，赋予冲击区域（第1个晶格）初速度v₀，v₀分别取1 km/s（压缩过程）和 –1 km/s（拉伸过程），损失因子η取0.8，时间步长 $\Delta t$ 取10^–16 s。特别地，为了便于对比压缩和拉伸过程，在拉伸时位移取相反数。首先考虑3种势作用下碳晶格在压缩过程中的波传播过程，取冲击区域方向的前5个晶格进行研究，X₁、X₂、X₃、X₄、X₅依次为5个晶格偏离平衡位置的位移，选取不同时刻来探究波的起始阶段。

图5给出了碳晶格在压缩冲击下的波传播过程。可以看出，在无缺陷的模型中，各个晶格初始处于平衡状态；随着运动的开始，冲击区域速度缓慢减小，第2个晶格在一定的时间间隔后开始运动，速度快速升高，直至一个峰值，相对地，冲击区域速度快速减小，直至到达第2个平衡位置；第2个晶格开始运动相同的时间间隔后，第3个晶格开始运动，速度快速升高，并呈现与第2个晶格相同的规律，而第2个晶格速度快速减小。随着波的传播，各个晶格依次开始运动，呈现加速-减速2个阶段，直至到达第2个平衡位置，随后在平衡位置处小幅度震荡。在这一过程中，波速为2个晶格的位移差与2个原子开始运动的时间差的比值。

图 5 不同势作用下碳晶格中波的传播特性

Figure 5. Propagation characteristics of waves in carbon lattice under different potentials

下载: 全尺寸图片幻灯片

缺陷的种类繁多，会对波的传播过程产生不同的影响。本文中，通过增大晶格间初始间距来实现缺陷，可以看出，缺陷前后晶格的位移峰值差异巨大，同时晶格到达第2个平衡位置后震荡更剧烈。在运动开始时，缺陷两端的晶格存在相反的运动趋势。

另一方面，3种势作用下碳晶格的运动过程存在差异。图6给出了线性作用和Tersoff势作用下碳晶格在压缩和拉伸冲击下波形的对比。可以看出，线性作用下，每个晶格只受到相邻晶格的作用，开始运动时，各晶格依次运动。此时，晶格间作用与晶格间距的比值不变，压缩与拉伸过程相比，同一时刻各晶格的位移大小相等，方向相反。Tersoff势作用于全场，开始运动时，所有晶格均开始运动，但由于非相邻晶格的相互作用远小于相邻晶格的作用，因而后方的晶格在开始时的位移几乎忽略不计。此外，Tersoff势作用可以更真实地反映晶格间距对晶格间作用的影响，图4反映了Tersoff势和线性作用的区别：相比于线性作用，在Tersoff势作用下，压缩过程中晶格间距减小，晶格间斥力急速增大，而拉伸过程中晶格间距增大，晶格间引力缓慢增大。碳晶格在压缩过程中晶格间距减小，晶格的加速度增大，相同时间下各晶格的位移随之减小；拉伸过程中，晶格间距增大，晶格的加速度变小，相同时间下晶格的位移随之增大。图7给出了含缺陷的Tersoff势作用下碳晶格在压缩和拉伸冲击下波形的对比。当存在缺陷时，缺陷处存在内力，缺陷两端的晶格初始具有相反的运动趋势，此时第3和第4个晶格初始的运动方向相反，而两者到达相同位移的时间差显著增大，这也导致了缺陷处波速的减小。

图 6 压缩和拉伸过程中波形的对比

Figure 6. Comparison of waveforms during compression and stretching processes

下载: 全尺寸图片幻灯片

图 7 含缺陷时压缩和拉伸过程中波形的对比

Figure 7. Comparison of waveforms during compression and stretching processes with defects

下载: 全尺寸图片幻灯片

2.3 硅晶格中的波传播特性

C和Si属于不同周期的同族元素，晶格结构类似，具有相似的性质。图8(a)～(b)分别给出了Tersoff势作用及含缺陷的Tersoff势作用下C晶格与Si晶格中波传播过程的对比，选取不同时刻来探究波的起始阶段。

图 8 C晶格和Si晶格中波的传播特性

Figure 8. Propagation characteristics of waves in carbon lattice and silicon lattice

下载: 全尺寸图片幻灯片

Si晶格的运动过程呈现与C晶格相似的规律。图8(a)为无缺陷时Tersoff势作用下C和Si晶格中波的传播过程。在Si晶格中，共价键较长，Si晶格的位移峰值较高，后方晶格出现明显位移的时间会更晚。图8(b)为含缺陷的Tersoff势作用下的晶格中波的传播过程，相比于C晶格，Si晶格中缺陷两端的晶格出现明显位移的时间更晚，到达相同位移的时间差更大。

入射速度的选取会对晶格的位移峰值和受力峰值产生影响。图9(a)和(b)分别给出了单晶体系中不同入射速度的C和Si晶格中前3个晶格的位移峰值和受力峰值的变化。可以看出，在2种晶格中，位移和受力峰值都与入射速度线性相关。这与宏观连续介质的波传播过程中应力与速度线性相关的结论一致。此外，不同晶格的位移峰值和受力峰值存在细微差别，相同入射速度下，Si晶格中的位移峰值和受力峰值都高于C晶格；提升相同的入射速度，Si晶格中的位移和受力峰值增长更快。

图 9 单晶中入射速度的影响

Figure 9. Influence of incident velocity in single crystals

下载: 全尺寸图片幻灯片

3. 多晶中波传播的分子动力学模拟

3.1 模型构建与模拟方法

基于广泛使用的LAMMPS平台，采用分子动力学方法研究波在金刚石和碳化硅中的传播，同时使用OVITO进行可视化分析，对波传播过程中原子的位移进行分析。碳和硅属于同族元素，而金刚石与碳化硅的结构相似，其晶格常数分别为0.357和0.436 nm。采用Tersoff势函数分别描述C－C键作用和C－Si键作用。

建立一个尺寸为30 nm×1 nm×1 nm的盒子，在中间20 nm×0.2 nm×0.2 nm的圆柱形区域中填充碳原子，得到金刚石模型，模型中包含478个碳原子，如图10(a)所示。金刚石模型中C－C键长为1.54 Å。模型分为两部分，分别为冲击区域（红色）与传播区域（蓝色），其中，冲击区域包含10个碳原子。通过赋予冲击区域初速度实现冲击。在施加初速度前，对模型进行能量最小化处理，消除模型的初始震荡。在冲击过程中，为了消除y轴和z轴作用的影响，采用三维方向自由边界，同时固定所有原子只沿x向运动。通过对冲击区域赋予1 和 –1 km/s的初速度，探究金刚石压缩和拉伸过程中原子的运动和波的传播。运动过程中系统处于NVE系综，时间步长为1 fs，运行1 000步。

图 10 金刚石和碳化硅的分子动力学模型

Figure 10. Molecular dynamics models of diamond and silicon carbide

下载: 全尺寸图片幻灯片

碳化硅晶体中Si－C键长大于C－C键长，相同结构时碳化硅的尺寸大于金刚石。建立同样大的盒子，在中间20 nm×0.3 nm×0.3 nm的圆柱形区域中建立碳化硅模型，模型中包含221个碳原子和176个硅原子，如图10(b)所示。碳化硅模型中Si－C键长为1.89 Å。采用与金刚石模型中相同的设置，其中，冲击区域包含6个碳原子和4个硅原子。通过赋予冲击区域初速度，来研究碳化硅压缩和拉伸过程中原子的运动和波的传播。

3.2 金刚石晶体中波的传播过程

多晶体系不同于单晶体系，存在更复杂的结构。为探究多晶体系中波的传播特性，考虑金刚石及含缺陷的金刚石在冲击过程中原子位移的变化，其中，通过删除金刚石中的2个原子来实现晶体中的缺陷，如图10(a)所示。图11给出了0.04、0.08和0.16 ps等3个时刻压缩冲击下不同原子的位移随时间的变化，选取冲击区域方向中间一列的前5个原子进行研究。当无缺陷时，随着波的传播，各个原子依次开始运动，呈现加速-减速的2个阶段，直至到达第2个平衡位置，随后在平衡位置处小幅度震荡，这与单晶体系中的结论一致。特别地，在每个时刻，近似有3个原子向第2个平衡位置运动，即传播过程中的波长近似是冲击区域长度的2倍。另外，多晶体系中原子到达第2个平衡位置后的震荡较小。

图 11 金刚石中波的传播特性

Figure 11. Propagation characteristics of waves in diamond

下载: 全尺寸图片幻灯片

图12给出了金刚石及含缺陷的金刚石中不同原子运动过程的对比。分别选取模型冲击区域方向的顶端一列的前5个原子和中间一列的前5个原子进行研究。无缺陷情况下，顶端原子和中间原子呈现相似的运动趋势，不同的是，对于中间一列原子，第1个原子与第2个原子开始时的初速度和位移近似相同，运动一段时间后才产生差异。另外，缺陷破坏了结构的规则性，造成了能量的衰减，导致初始时刻缺陷两端原子的运动方向相反。顶端一列原子受缺陷影响较大，缺陷两端原子的偏离位移较大；而中间一列原子受缺陷影响较小，缺陷两端原子的偏离位移较小。

图 12 金刚石中间与顶端原子的运动过程

Figure 12. Motion process of the middle and top atoms in diamond

下载: 全尺寸图片幻灯片

图13给出了压缩和拉伸冲击下金刚石及含缺陷金刚石的中间一列原子的运动过程，拉伸时位移取相反数。可以看出，压缩过程中原子受力更大，原子的位移峰值小于拉伸过程。特别地，对于含缺陷的金刚石，在压缩和拉伸冲击下波的传播呈现较大差异：缺陷处原子间作用减小，而在压缩过程中原子间斥力快速增大，会弱化缺陷的影响，使缺陷处原子达到相同位移的时间差微弱增大；在拉伸过程中，原子间引力缓慢增大，缺陷会进一步减小原子间作用，原子的速度变化更慢，原子的位移峰值急剧增大，缺陷前后原子达到相同位移的时间差大幅增加。

图 13 压缩和拉伸过程中金刚石中的波传播特性

Figure 13. Wave propagation characteristics of diamond during compression and stretching processes

下载: 全尺寸图片幻灯片

3.3 碳化硅晶体中波的传播过程

图14(a)给出了无缺陷金刚石和碳化硅的中间一列原子的运动过程。在波的传播过程中，二者总体的运动规律相似，不同的是，碳化硅中键长更长，波传播过程中原子的位移峰值更大，相邻原子到达相同位移的时间更长，而波速更小。图14(b)给出了存在缺陷的金刚石和碳化硅中波的传播过程。碳化硅相对质量较大，但弹性极限较低，此时缺陷对碳化硅的影响更显著，缺陷前后原子与无缺陷情况相比，偏离的位移更大。

图 14 金刚石与碳化硅中波传播过程的对比

Figure 14. Comparison of wave propagation processes in diamond and silicon carbide

下载: 全尺寸图片幻灯片

此外，考虑入射端初速度对多晶体系中波传播的影响，图15(a)和(b)分别给出了金刚石及碳化硅中入射速度对不同原子位移峰值和受力峰值的的影响。取中间一列的前3个原子进行研究，可以看出：金刚石中位移和受力峰值都与入射速度线性相关；碳化硅中键长更长，相同入射速度下原子的位移峰值和受力峰值都高于金刚石，但碳化硅强度更低，在入射速度高于1 km/s后呈现非线性规律，位移峰值增长放缓，受力峰值急剧增加。

图 15 多晶体系中入射速度对波传播的影响

Figure 15. Influences of incident velocity on wave propagation in polycrystalline systems

下载: 全尺寸图片幻灯片

4. 结　论

从晶格的微振动出发，基于3种势作用，得到了格波波速、角频率与波数的关系。基于有限差分方法，以C和Si晶格为例，分析了单晶中波的传播过程，讨论了入射速度对波传播中位移峰值和受力峰值的影响。结合分子动力学方法，研究了金刚石和碳化硅两种多晶体系中波的传播过程。得到以下主要结论。

(1) 在格波的传播中，原子间在不同的势作用下，格波波速、角频率与波数的关系存在差异。在线性作用下，格波波速和角频率与波数的关系不受原子位移的影响；在Tersoff势作用下，格波波速和角频率与原子的位移呈正相关；当存在位错时，体系能量减小，格波波速和角频率也随之降低。

(2) 在单晶体系的波传播中，晶格的运动依次呈现加速-减速2个阶段。Tersoff势作用相比于线性作用，能够更准确地描述晶格间距对晶格间作用力的影响；位错的存在导致内力产生，也引发了位错两端晶格初始的相反运动，同时造成波速的减小。相比于压缩过程，拉伸过程中晶格受力更小，位移峰值更大。相比于碳晶格，硅晶格位移峰值较大，而波速较低。在2种晶格中，位移峰值和受力峰值都与入射速度线性相关，与宏观连续介质中波传播的规律类似。

(3) 在多晶体系的波传播中，原子的运动呈现与单晶体系中相似的规律：原子的运动呈现加速-减速2个阶段，拉伸过程的位移峰值大于压缩过程。无缺陷时顶端一列原子与中间一列原子的运动规律相似，而存在缺陷时顶端一列原子受影响较大。压缩和拉伸冲击下，含缺陷的金刚石中波的传播呈现较大差异：压缩会弱化缺陷的影响，而拉伸会强化缺陷的影响，拉伸冲击下原子的最大位移大幅增大。在入射速度低于1 km/s时，金刚石和碳化硅中位移和受力峰值都与入射速度线性相关；入射速度高于1 km/s时，碳化硅中呈现非线性规律，位移峰值增长缓慢，受力峰值急剧增大。

图 1 一维原子链振动模型^[27]

Figure 1. Vibration model of one-dimensional atomic chain^[27]

下载: 全尺寸图片幻灯片

图 2 Tersoff势下角频率与波数的关系

Figure 2. Relationship between angular frequency and wavenumber under Tersoff potential

下载: 全尺寸图片幻灯片

图 3 含缺陷的Tersoff势下角频率与波数的关系

Figure 3. Relationship between angular frequency and wavenumber under Tersoff potential with defects

下载: 全尺寸图片幻灯片

图 4 碳晶格中势力与晶格间距的关系

Figure 4. Relationship between force and lattice spacing in carbon lattices

下载: 全尺寸图片幻灯片

图 5 不同势作用下碳晶格中波的传播特性

Figure 5. Propagation characteristics of waves in carbon lattice under different potentials

下载: 全尺寸图片幻灯片

图 6 压缩和拉伸过程中波形的对比

Figure 6. Comparison of waveforms during compression and stretching processes

下载: 全尺寸图片幻灯片

图 7 含缺陷时压缩和拉伸过程中波形的对比

Figure 7. Comparison of waveforms during compression and stretching processes with defects

下载: 全尺寸图片幻灯片

图 8 C晶格和Si晶格中波的传播特性

Figure 8. Propagation characteristics of waves in carbon lattice and silicon lattice

下载: 全尺寸图片幻灯片

图 9 单晶中入射速度的影响

Figure 9. Influence of incident velocity in single crystals

下载: 全尺寸图片幻灯片

图 10 金刚石和碳化硅的分子动力学模型

Figure 10. Molecular dynamics models of diamond and silicon carbide

下载: 全尺寸图片幻灯片

图 11 金刚石中波的传播特性

Figure 11. Propagation characteristics of waves in diamond

下载: 全尺寸图片幻灯片

图 12 金刚石中间与顶端原子的运动过程

Figure 12. Motion process of the middle and top atoms in diamond

下载: 全尺寸图片幻灯片

图 13 压缩和拉伸过程中金刚石中的波传播特性

Figure 13. Wave propagation characteristics of diamond during compression and stretching processes

下载: 全尺寸图片幻灯片

图 14 金刚石与碳化硅中波传播过程的对比

Figure 14. Comparison of wave propagation processes in diamond and silicon carbide

下载: 全尺寸图片幻灯片

图 15 多晶体系中入射速度对波传播的影响

Figure 15. Influences of incident velocity on wave propagation in polycrystalline systems

下载: 全尺寸图片幻灯片

表 1 Tersoff势函数中的材料参数^[30]

Table 1. Material parameters in Tersoff potential^[30]

共价键种类	m₁	γ	λ₃/Å^–1	C	d	cosθ₀	n₁
C―C	3.0	1.0	0	38 049	4.348 4	–0.570 58	0.727 51
Si―Si	3.0	1.0	0	100 390	16.217 0	–0.598 25	0.787 34

共价键种类	β	λ₂/Å^–1	B/eV	R/Å	D/Å	λ₁/Å^–1	A₁/eV
C―C	1.572 4×10^–7	2.211 90	346.70	1.95	0.15	3.487 9	1 393.6
Si―Si	1.100 0×10^–6	1.732 22	471.18	2.85	0.15	2.479 9	1 830.8

下载: 导出CSV

表 2 2种晶格在平衡位置处的参数

Table 2. Parameters of two types of lattices at equilibrium positions

原子种类	摩尔质量/(g·mol⁻¹)	平衡距离/Å	晶格常数/Å
C	12	1.54	3.57
Si	28	2.35	5.43

下载: 导出CSV

参考文献(30)

[1]	徐松林, 刘永贵, 席道瑛. 岩石物理与动力学原理 [M]. 北京: 科学出版社, 2019: 1–21. XU S L, LIU Y G, XI D Y. Rock physics and dynamics principle [M]. Beijing: Science Press, 2019: 1–21.
[2]	王礼立. 应力波基础 [M]. 2版. 北京: 国防工业出版社, 2005: 1–28. WANG L L. Foundation of stress waves [M] 2nd ed. Beijing: National Defense Industry Press, 2005: 1–28.
[3]	袁良柱, 陆建华, 苗春贺, 等. 基于分数阶模型的牡蛎壳动力学特性研究 [J]. 爆炸与冲击, 2023, 43(1): 011101. DOI: 10.11883/bzycj-2022-0318. YUAN L Z, LU J H, MIAO C H, et al. Dynamic properties of oyster shells based on a fractional-order model [J]. Explosion and Shock Waves, 2023, 43(1): 011101. DOI: 10.11883/bzycj-2022-0318.
[4]	袁良柱, 苗春贺, 单俊芳, 等. 冲击下混凝土试样应变率效应和惯性效应探讨 [J]. 爆炸与冲击, 2022, 42(1): 013101. DOI: 10.11883/bzycj-2021-0114. YUAN L Z, MIAO C H, SHAN J F, et al. On strain-rate and inertia effects of concrete samples under impact [J]. Explosion and Shock Waves, 2022, 42(1): 013101. DOI: 10.11883/bzycj-2021-0114.
[5]	CHEN M D, XU S L, YUAN L Z, et al. Influence of stress state on dynamic behaviors of concrete under true triaxial confinements [J]. International Journal of Mechanical Sciences, 2023, 253: 108399. DOI: 10.2139/ssrn.4332010.
[6]	JANG S, RABBANI M, OGRINC A L, et al. Tribochemistry of diamond-like carbon: interplay between hydrogen content in the film and oxidative gas in the environment [J]. ACS Applied Materials & Interfaces, 2023, 15(31): 37997–38007. DOI: 10.1021/acsami.3c05316.
[7]	HU L F, ZHAI X Y, LI J G, et al. Improving the mechanical properties and tribological behavior of sulfobetaine polyurethane based on hydrophobic chains to be applied as artificial meniscus [J]. ACS Applied Materials & Interfaces, 2023, 15(25): 29801–29812. DOI: 10.1021/acsami.3c02940.
[8]	WANG D Y, WANG P F, WU Y F, et al. Temperature and rate-dependent plastic deformation mechanism of carbon nanotube fiber: experiments and modeling [J]. Journal of the Mechanics and Physics of Solids, 2023, 173: 105241. DOI: 10.1016/j.jmps.2023.105241.
[9]	薛晓. 碳纳米管纤维的动静态力学性能研究 [D]. 合肥: 中国科学技术大学, 2020: 53–61. DOI: 10.27517/d.cnki.gzkju.2020.000565. XUE X. Investigation of dynamic and quasi-static mechanical properties of carbon nanotube fibers [D]. Hefei: University of Science and Technology of China, 2020: 53–61. DOI: 10.27517/d.cnki.gzkju.2020.000565.
[10]	MACHADO M, MOREIRA P, FLORES P, et al. Compliant contact force models in multibody dynamics: evolution of the Hertz contact theory [J]. Mechanism and Machine Theory, 2012, 53: 99–121. DOI: 10.1016/j.mechmachtheory.2012.02.010.
[11]	BORODICH F M. The Hertz-type and adhesive contact problems for depth-sensing indentation [J]. Advances in Applied Mechanics, 2014, 47: 225–366. DOI: 10.1016/b978-0-12-800130-1.00003-5.
[12]	YANG F, XIE W H, MENG S H. Impact and blast performance enhancement in bio-inspired helicoidal structures: a numerical study [J]. Journal of the Mechanics and Physics of Solids, 2020, 142: 104025. DOI: 10.1016/j.jmps.2020.104025.
[13]	PENG Q, LIU X M, WEI Y G. Elastic impact of sphere on large plate [J]. Journal of the Mechanics and Physics of Solids, 2021, 156: 104604. DOI: 10.1016/j.jmps.2021.104604.
[14]	TANG X, YANG J. Wave propagation in granular material: what is the role of particle shape? [J]. Journal of the Mechanics and Physics of Solids, 2021, 157: 104605. DOI: 10.1016/j.jmps.2021.104605.
[15]	ALBERDI R, ROBBINS J, WALSH T, et al. Exploring wave propagation in heterogeneous metastructures using the relaxed micromorphic model [J]. Journal of the Mechanics and Physics of Solids, 2021, 155: 104540. DOI: 10.1016/j.jmps.2021.104540.
[16]	WAYMEL R F, WANG E, AWASTHI A, et al. Propagation and dissipation of elasto-plastic stress waves in two dimensional ordered granular media [J]. Journal of the Mechanics and Physics of Solids, 2018, 120: 117–131. DOI: 10.1016/j.jmps.2017.11.007.
[17]	LI S F, WANG G. Introduction to micromechanics and nanomechanics [M]. Singapore: World Scientific Publishing Co. Pre. Ltd., 2008.
[18]	ZUNDEL L, MALONE K, CERDÁN L, et al. Lattice resonances for thermoplasmonics [J]. ACS Photonics, 2023, 10(1): 274–282. DOI: 10.1021/acsphotonics.2c01610.
[19]	CERDÁN L, ZUNDEL L, MANJAVACAS A. Chiral lattice resonances in 2.5-dimensional periodic arrays with achiral unit cells [J]. ACS Photonics, 2023, 10(6): 1925–1935. DOI: 10.1021/acsphotonics.3c00369.
[20]	HU Y G, LIEW K M, WANG Q, et al. Nonlocal shell model for elastic wave propagation in single-and double-walled carbon nanotubes [J]. Journal of the Mechanics and Physics of Solids, 2008, 56(12): 3475–3485. DOI: 10.1016/j.jmps.2008.08.010.
[21]	AGARWAL G, VALISETTY R R, DONGARE A M. Shock wave compression behavior and dislocation density evolution in Al microstructures at the atomic scales and the mesoscales [J]. International Journal of Plasticity, 2020, 128: 102678. DOI: 10.1016/j.ijplas.2020.1026.
[22]	SAM A, ÁLVAREZ M B, VENEGAS R, et al. Multiscale acoustic properties of nanoporous materials: from microscopic dynamics to mechanics and wave propagation [J]. The Journal of Physical Chemistry C, 2023, 127(15): 7471–7483. DOI: 10.1021/acs.jpcc.3c00060.
[23]	薛定谔. 薛定谔讲演录 [M]. 2版. 范岱年, 胡新和, 译. 北京: 北京大学出版社, 2019: 7–8. SCHRÖDINGER E. Lectures of Schrödinger [M]. 2nd ed. Translated by FAN D N, HU X H. Beijing: Peking University Press, 2019: 7–8.
[24]	TOLOS L, CENTELLES M, RAMOS A. The equation of state for the nucleonic and hyperonic core of neutron stars [J]. Publications of the Astronomical Society of Australia, 2017, 34: e065. DOI: 10.1017/pasa.2017.60.
[25]	AARABI M, SARKA J, PANDEY A, et al. Quantum dynamical investigation of dihydrogen-hydride exchange in a transition-metal polyhydride complex [J]. The Journal of Physical Chemistry A, 2023, 127(31): 6385–6399. DOI: 10.1021/acs.jpca.3c01863.
[26]	HO W W, CHOI S. Exact emergent quantum state designs from quantum chaotic dynamics [J]. Physical Review Letters, 2022, 128(6): 060601. DOI: 10.1103/PhysRevLett.128.060601.
[27]	黄昆. 固体物理学 [M]. 北京: 人民教育出版社, 1966: 35–43.
[28]	王礼立, 胡时胜, 杨黎明, 等. 材料动力学 [M]. 合肥: 中国科学技术大学出版社, 2017: 33–158.
[29]	BRENNER D W. Tersoff-type potentials for carbon, hydrogen and oxygen [J]. MRS Online Proceedings Library (OPL), 1988, 141: 59. DOI: 10.1557/proc-141-59.
[30]	TERSOFF J. Modeling solid-state chemistry: interatomic potentials for multicomponent systems [J]. Physical Review B, 1989, 39(8): 5566. DOI: 10.1103/PhysRevB.39.5566.