含能材料冲击起爆的计算机模拟方法

摘要

本发明请求保护一种含能材料在冲击波引导条件下起爆的高效且相对高精度的量子分子动力学模型。该模型是一种直接动力学模拟方法，即不预先构建势能面，体系的能量和作用于原子核的力采取随用随算模式。原子之间的作用力以及电子结构信息的计算采用电荷自洽的基于紧束缚近似的密度泛函方法(DFTB)，对于冲击波的描述选择多尺度近似方法。模拟结果将输出原子的运动规矩，及物理状态信息，比如温度、体积、压力等。该方法将为现有含能材料的改进，新型含能材料的筛选、设计、优化等提供先期的理论预测和科学依据。

说明书

技术领域

本发明属于材料学，计算机科学领域，具体涉及一种含能材料在冲击条件下起爆 的计算机模拟算法。

背景技术

随着安全问题在当代含能材料技术中发挥越来越重要的作用，对暴露于极端环境 中的固体猛炸药的相对安全性的关注日益增加，冲击波引爆炸药现象的研究是其中一个重 要方面。

冲击波通过受体炸药时，炸药承受压缩和绝热加热，引起化学反应，炸药分子分 解，释放化学能，同时使得波由降速变为加速，形成高速爆轰波。整个冲击起爆过程发生在 皮秒至纳秒时间尺度内，对实验研究提出了非常高的要求，测量系统必须具有极强的快速 响应能力和很高分辨率。随着实验技术的发展，科学家虽然可以直接研究爆轰过程中含能 材料的化学性质^[1]，但是，这样的实验非常昂贵，且很难真正达到时间尺度范围之内^[2]。理 论上常用经典分子动力学(ClassicalMolecularDynamics,CMD)模拟方法研究含能材料 的性质，包括平衡状态下的热力学、结构和力学性质等^[3]。CMD模拟通过求解经典运动方程 得到每一时刻材料中原子的位置和速度，从而使得我们能够在原子水平上直接观察含能材 料的许多超快动态性质。由于CMD模拟计算需要预先知道体系中粒子间相互作用的势函数， 为研究含能材料的化学性质，国际上已经发展出不同的分子力场来描述有机分子在高温高 压时的化学反应性能，影响较大的有反应性经验键级(ReactiveEmpiricalBondOrder, REBO)势函数^[4]和反应性分子力场(ReactiveForceFields,ReaxFF)^[5]。国内外多个研究 团队应用这些力场函数研究了RDX、PETN等材料在冲击波作用下的化学反应及爆轰波中能 量转换和传递机制^[6-9]。这些研究结果对了解含能材料的化学性质有着重要的作用。

2006年，美国陆军研究实验室(U.S.ArmyResearchLaboratory)的Rice等人^[10]提出用量子分子动力学(QuantumMolecularDynamics,QMD)方法模拟冲击波引起的含能 材料的化学和物理性质的变化。QMD计算应用量子力学方法处理电子的运动，原子核的运动 则由经典力学近似处理。虽然忽略了核运动的某些量子效应，但是QMD方法能够正确的处理 含能材料分子的电子结构，特别是化学键的断裂和形成。其选择的量子化学计算方法为密 度泛函理论(densityfunctionaltheory,DFT)。Rice^[11]等人应用DFT近似模拟了高能量 密度材料PETN和cg-N在冲击波加载下的动力学行为。他们的模拟结果从分子水平上描述了 这两种材料在能量释放时化学和物理性质的变化，明确地显示出两种不同的反应机理。最 近，南京理工大学和上海交通大学也开始运用从头算分子动力学模拟(AbInitio MolecularDynamicsSimulations)方法对含能材料体系的研究进行了有益的尝试^[12-13]。 可惜的是，由于从头算方法巨大的计算量，模拟的过程仅能到达几个皮秒。这与含能材料爆 炸过程化学反应的时间相差甚远，因而不能揭示引爆时含能材料性质变化的全貌。譬如，在 Rice等人工作中，体系包含2303个原子，模拟计算在多达512个CPU的超级计算机上进行，仅 完成了4000个时间步长(每一个时间步长设为1飞秒，共计4个皮秒)。针对这一问题，美国劳 伦斯利弗莫尔国家实验室的Reed等人提出应用基于紧束缚近似的密度泛函(Density FunctionalbasedTightBonding，DFTB)，该方法已经被证明其计算所得的分子结构和能 量数据能够比拟用abinitial方法取得的结果，并与分子体系和凝聚态的实验数据一致。

基于冲击Hugoniot的分子动力学方法可通过固定系统体积的同时约束温度对材 料进行微观模拟，但是该方法不能俘获过程中涉及的化学反应路径等^[14]。非平衡分子动力 学(NEMD)能够揭示冲击波结构中的很多微观细节^[15]，但所建立的模型长度需要满足真实 爆轰形成的尺度要求，并且模型的尺度还依赖于冲击波速度的大小，因此需要大量的计算 资源。多尺度冲击技术(MSST)^[16-18]采用可压缩流Navier-Stokes方程，较采用NEMD方法模 拟冲击波传播所需的模型尺寸显著降低，并且与NEMD模拟结果具有较好的一致性^[19]。

通过对当前含能材料计算机模拟相关技术研究成果的归纳与总结可以发现，目 前主要存在如下几个问题。

1.目前对含能材料爆轰的计算机模拟还大量采用的是经典的分子动力学技术。国 内外多个研究团队应用力场函数研究了RDX、PETN等材料在冲击波作用下的化学反应及爆 轰波中能量转换和传递机制^[20-23]。然而分子动力学模拟方法的一个重要缺陷是其对于体 系中粒子相互作用的势函数的依赖性。粒子相互作用的势函数是通过对实验数据或者高水 平的量子化学计算结果拟合而得到。模拟结果的可信性取决于所采用的经验势函数是否能 够准确的描述含能材料的化学性质。事实上，用参数化的经验势函数是很难准确描述高温 高压等极端条件下分子内的相互作用，完善现有的模型最多能够使其适用于特定的体系和 特定的条件。化学反应中存在化学键的断裂和重新形成过程均涉及电子的相互作用，因此 精确的含能材料的化学性质应该用量子化学的方法处理。

2.基于紧束缚近似的密度泛函计算所得的分子结构和能量数据虽然能够比拟用 abinitial方法取得的结果，但是不能处理极端条件下原子的运动轨迹。

3.多尺度冲击技术虽然能够处理冲击波条件下原子的运动轨迹，但是不能处理含 能材料分子的电子结构，特别是化学键的断裂和形成。

参考文献

[1]TeipelU(主编)，欧育湘(主译)，含能材料.北京:国防工业出版社，2009,471- 480.

[2]HolianBL,GermannTC,StrachanA,andMailletJB,"Non-Equilibrium MolecularDynamicsStudiesofShockandDetonationProcessesinEnergetic Materials,"inChemistryatExtremeConditions.chapter.9,ManaaMR,Ed.,New York:Elsevier,2005,269-298.

[3]RiceBMandSewellTD,"EquilibriumMolecularDynamics Simulations,"inStaticCompressionofEnergeticMaterials.chapter.7,Berlin: Springer,2008.

[4]BrennerDW,ShenderovaOA,HarrisonJA,StuartSJ,NiB,andSinnott SB,"Asecond-generationreactiveempiricalbondorder(REBO)potentialenergy expressionforhydrocarbons,"JournalofPhysics-CondensedMatter,2002,14,783- 802.

[5]vanDuinACT,DasguptaS,LorantF,andGoddardWA,"ReaxFF:Areactive forcefieldforhydrocarbons,"JournalofPhysicalChemistryA,2001,105,9396- 9409.

[6]StrachanA,vanDuinACT,ChakrabortyD,DasguptaS,andGoddardWA," Shockwavesinhigh-energymaterials:Theinitialchemicaleventsinnitramine RDX,"PhysicalReviewLetters,2003,91.

[7]NomuraKI,KaliaRK,NakanoA,andVashishtaP,"Ascalableparallel algorithmforlarge-scalereactiveforce-fieldmoleculardynamics simulations,"ComputerPhysicsCommunications,2008,178,73-87.

[8]LiuLC,LiuY,ZybinSV,SunH,andGoddardWA,"ReaxFF-/g:Correction oftheReaxFFReactiveForceFieldforLondonDispersion,withApplicationsto theEquationsofStateforEnergeticMaterials,"JournalofPhysicalChemistry A,2011,115,11016-11022.

[9]XuJ,ZhaoJ,andSunL,"ThermaldecompositionbehaviourofRDXby first-principlesmoleculardynamicssimulation,"MolecularSimulation,2008,34, 961-965.

[10]RomeroNA,MattsonWD,andRiceBM,"UsingQuantumMechanicsto PredictShockPropertiesofExplosives,"U.S.ArmyResearchLaboratory,Aberdeen ProvingGround,MD21005-5066,01NOV2006.

[11]MattsonWD,BaluR,andRiceBM,"DirectQuantumMechanical SimulationsofShockedEnergeticMaterials"2008DoDHPCMPUsersGroup Conference,2008.

[12]ChangJ,LianP,WeiDQ,ChenXR,ZhangQM,andGongZZ,"Thermal DecompositionoftheSolidPhaseofNitromethane:AbInitioMolecularDynamics Simulations,"PhysicalReviewLetters,2010,105,188302.

[13]朱卫华，肖鹤鸣，“从头算分子动力学模拟冲击波引爆炸药的机理”，第十一届 全国量子化学会议，合肥，2011.

[14]MailletJB,MareschalM,SoulardL,etal.UniaxialHugoniostat:A methodforatomisticsimulationsofshockedmaterials[J].PhysicalReviewE, 2000,63(1):016121.

[15]HeimAJ,JensenNG,KoberEM,GermannTC2008Phys.Rev.E78046710

[16]ReedEJ,FriedLE,JoannopoulosJD.Amethodfortractable dynamicalstudiesofsingleanddoubleshockcompression[J].Physicalreview letters,2003,90(23):235503.

[17]ReedEJ,FriedLE,ManaaMR,etal.Amulti-scaleapproachto moleculardynamicssimulationsofshockwaves[J].ChemistryatExtreme Conditions,2005:297-326.

[18]ReedEJ,FriedLE,HenshawWD,etal.Analysisofsimulation techniqueforsteadyshockwavesinmaterialswithanalyticalequationsof state[J].PhysicalReviewE,2006,74(5):056706.

[19]ReedEJ,MaitiA,FriedLE.Anomaloussoundpropagationandslow kineticsindynamicallycompressedamorphouscarbon[J].PhysicalReviewE, 2010,81(1):016607.

[20]ManaaMR,ReedEJ,FriedLE,andGoldmanN,"Nitrogen-Rich HeterocyclesasReactivityRetardantsinShockedInsensitiveExplosives," JournaloftheAmericanChemicalSociety,2009,131,5483-5487.

[21]MaxwellR.DissolvingMoleculestoImproveTheirPerformance [Online].Available:https://str.llnl.gov/June09/maxwell.html

[22]DouYS,TorralvaBR,andAllenRE,"SemiclassicalElectron- Radiation-IonDynamics(SERID)andCis-TransPhotoisomerizationofButadiene," JournalofModernOptics,2003,50,2615-2643.

[23]TandAllenRE,"Femtosecond-scaleresponseofGaAsto ultrafastlaserpulses,"PhysicalReviewB,2002,66,081202.

发明内容

针对以上现有技术的不足，提出了一种方法。本发明的技术方案如下：一种

含能材料冲击起爆的计算机模拟方法，其包括以下步骤；

步骤1、用MaterialsStudio软件建立含能材料冲击起爆模拟体系的三维模型，并 导出所有原子的三维坐标，

步骤2、确定模拟参数，包括初始温度、初始压力、迭代时间、时间步长、冲击速度， 模拟计算并输出数据，具体包括：

101、根据玻尔兹曼分布随机抽样获得扰动初始速度；

102、步骤101获得原子的初始速度后，采用基于紧束缚近似的密度泛函SCC-DFTB 方法计算迭代的第n步所有原子间的作用力来计算原子的速度；

103、完成步骤102计算所有原子间的作用力后，根据公式(5)计算第n+1/2步所有 原子的速度；

104、根据公式r_i＝r_i+v_iΔt计算第n+1步的三维坐标位置，其中r_i表示第i个粒子 的位置，v_i为上一步计算得到的速度，Δt为设定的时间步长；

105、根据公式(5)计算第n+1/2步所有原子的速度；

106、根据公式(2)、公式(6)以及公式(7)计算并输出体系的压力、体积、能量、温 度。

进一步的，步骤101中扰动初始速度为：

令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得 到。

进一步的，步骤102中原子间作用力的计算具体为：

原子之间的相互作用力由电荷自洽的基于紧束缚近似的密度泛函SCC-DFTB确定；

$E=\underset{i}{\overset{occ}{\Sigma }}<{\phi }_i\left|\widehat{H^0}\right|{\phi }_i>+\frac{1}{2}\int \int (\frac{1}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}+\frac{{\delta }^2{E}_{ex}}{{\mathrm{\delta \rho \delta \rho }}^{\prime }}){\mathrm{\delta \rho \delta \rho }}^{\prime }-\frac{1}{2}\int \int \frac{{\rho }_0^{\prime }\left(\overrightarrow{r^{\prime }}\right){\rho }_0\left(\overrightarrow{r}\right)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}-\int V_{xc}\left[{\rho }_0\right]{\rho }_0+E_{core}---\left(1\right)$

公式(1)中，φ_i表示Kohn-Sham的轨道函数，代表参考密度ρ₀有效的KS哈密顿 量，V_xc和E_ex两项分别表示交换相关势和交换相关能，E_core表示原子核之间的斥力能量。其中 电荷密度起伏δρ＝ρ-ρ₀，其在二级项中表示为进一步的，步骤103、104、105中 原子速度以及位置的计算为：

建立分子动力学模拟体系的哈密顿量H^MD：

$H^{MD}=\underset{l,\alpha }{\Sigma }\frac{p_{l,\alpha }^2}{2m_l{a}_{\alpha }^2}+\phi \left(As\right)+\frac{p_{\alpha x}^{2}M}{2Q{\left(a_y{a}_z\right)}^2}-\frac{1}{2}{\mathrm{MV}}_s^2{(1-\frac{a_x}{a_{x,0}})}^2-p_0{a}_y{a}_z(a_{x,o}-a_x)---\left(2\right)$

公式(2)中，粒子i在a＝{x,y,z}方向的动量是计算体系是正交体 系，体系的三个边长分别是a_x，a_y，和a_z；计算体系晶胞向量a_x的动量是公式(4)中，p是压力，v_s是冲击波传播速度，是体系的总质量，Q是与模拟体系的质 量相关的参量，所有下标为0的参量表示冲击波以前的性，s是经过矩阵A变化后的位置向 量，在公式(3)中，右边第一项表示原子的动能，第一项表示原子的势能，第三项表示模拟体 系体积变化的动能，最后两项是模拟体系受外界压力后变化的势能。

原子速度的计算方程如下：

$V=V+(\frac{F}{m}+\frac{\alpha }{V_{sum}m})\Delta t+\frac{1}{2}[{\left(\frac{\alpha }{V_{sum}m}\right)}^2+\frac{F\alpha }{V_{sum}{m}^2}]{\mathrm{\Delta t}}^2---\left(5\right)$

公式(5)中，V表示原子的速度，F表示原子受到的力，m表示原子的质量，α表示材料 的黏度系数，V_sum表示所有原子的速度之和，Δt表示设定的步长时间。

以横冲击速度在材料内传播的波阵面前后的能量关系式为其中E是能量，P是在该冲击方向上的负的应力张量的对角分量，V是体积，下标0表示冲击之 前的量，无下标表示冲击后的量。用于约束压力的方程为P-P₀＝U⁷ρ₀(1-ρ₀/ρ)(7)，其中，U是 冲击速度，ρ是密度。

本发明的优点及有益效果如下：

1.采用SCC-DFTB近似计算原子之间的作用力，计算速度比完全的DFT计算快几个 数量级，而且计算所得的分子结构和能量数据与分子体系和凝聚态的实验数据一致，从而 有效地提高了模拟的计算速度。

2.采用多尺度近似来描述模拟过程中冲击波的传播，既不增加模拟体系的体积， 又解决了冲击波在遇到模拟体系边界时的反射问题。

3.应用量子分子动力学的方法研究含能材料引爆时的化学反应机理。与分子动力 学模拟方法相比，量子分子动力学方法不需要提前构建势能函数，因此可以得出未经人为 预测的过程，从而能真正从原子水平上了解含能材料在极端条件下的化学反应。

附图说明

图1是本发明提供优选实施例分子动力学工作方框图；

图2：本发明的分子动力学算法流程图；

图3：RDX模拟体系的三维图；

图4：不同冲击波加载后的温度随时间变化图；

图5：不同冲击波加载后的压力随时间变化图；

图6：不同冲击波加载后的体积随时间变化图；

图7：RDX在8km/s冲击下的断键过程。

具体实施方式

以下结合附图，对本发明作进一步说明：

如图1、2所示，本发明提供了一种含能材料在冲击波引导条件下起爆的高效且相 对高精度的量子分子动力学模型，算法过程如下。

1.1如图3所示，采用MaterialsStudio软件建立RDX模拟体系，确定原子的初始位 置；

1.2设定初始温度，初始压力，迭代时间，时间步长，冲击速度；

1.3计算扰动初始速度；

1.4计算第n步所有原子间的作用力；

1.5计算第n+1/2步所有原子的速度；

1.6计算第n+1步的三维坐标位置；如图6所示为模拟体系在一定冲击速度下的三 维快照图，从图中可以看出有化学键的断裂。

1.7计算第n+1/2所有原子的速度；

1.8计算并输出体系的温度(图3)、压力(图4)、体积(图5)随时间变化图；

1.9重复第4至7步直到达到迭代的时间。

令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得 到。分子动力学模拟体系的哈密顿量H^MD描述为：

公 式中，粒子i在a＝{x,y,z}方向的动量是计算体系是正交体系，体系的三个边 长分别是a_x，a_y，和a_z；计算体系晶胞向量a_x的动量是在这些公式中，p是 压力，v_s是冲击波传播速度，是体系的总质量，Q是与模拟体系的质量相关的参量 (该参量的单位是(质量)²(长度)^-4。所有下标为0的参量表示冲击波以前的性，s是经过矩阵 A变化后的位置向量。在公式(1)中，右边第一项表示原子的动能，第一项表示原子的势能， 第三项表示模拟体系体积变化的动能，最后两项是模拟体系受外界压力后变化的势能。

原子速度的计算方程如下：

$V=V+(\frac{F}{m}+\frac{\alpha }{V_{sum}m})\Delta t+\frac{1}{2}[{\left(\frac{\alpha }{V_{sum}m}\right)}^2+\frac{F\alpha }{V_{sum}{m}^2}]{\mathrm{\Delta t}}^2---\left(2\right)$

公式(2)中，V表示原子的速度，F表示原子受到的力，m表示原子的质量，α表示材料 的黏度系数，V_sum表示所有原子的速度之和，Δt表示设定的步长时间。

原子之间的相互作用力可以由电荷自洽的基于紧束缚近似的密度泛函(self- consistentchargedensityfunctionalbasedtightbonding，SCC-DFTB)确定。

$E=\underset{i}{\overset{occ}{\Sigma }}<{\phi }_i\left|\widehat{H^0}\right|{\phi }_i>+\frac{1}{2}\int \int (\frac{1}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}+\frac{{\delta }^2{E}_{ex}}{{\mathrm{\delta \rho \delta \rho }}^{\prime }}){\mathrm{\delta \rho \delta \rho }}^{\prime }-\frac{1}{2}\int \int \frac{{\rho }_0^{\prime }\left(\overrightarrow{r^{\prime }}\right){\rho }_0\left(\overrightarrow{r}\right)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}-\int V_{xc}\left[{\rho }_0\right]{\rho }_0+E_{core}---\left(3\right)$

公式(3)中φ_i表示Kohn-Sham的轨道函数，代表参考密度ρ₀有效的KS哈密顿量， V_xc和E_ex两项分别表示交换相关势和交换相关能，E_core表示原子核之间的斥力能量。其中电 荷密度起伏δρ＝ρ-ρ₀，其在二级项中表示为

以固定冲击速度在材料内传播的波阵面前后的能量关系式为其 中E是能量，P是在该冲击方向上的负的应力张量的对角分量，V是体积。下标0表示冲击之前 的量，无下标表示冲击后的量。用于约束压力的方程为P-P₀＝U⁷ρ₀(1-ρ₀/ρ)。其中，U是冲击 速度，ρ是密度。

以上这些实施例应理解为仅用于说明本发明而不用于限制本发明的保护范围。在 阅读了本发明的记载的内容之后，技术人员可以对本发明作各种改动或修改，这些等效变 化和修饰同样落入本发明权利要求所限定的范围。