Method for calculating elementary reaction rate constants in hexanitrohexaazaisowurtzitane decomposition networks

Abstract

The application discloses a method for calculating the elementary reaction rate constant in a hexanitrohexaazaisowurtzitane decomposition network. The method realizes the calculation of the micro elementary reaction rate constant by identifying the generation and annihilation of single molecules through molecular dynamics simulation and using full lifetime statistics. On this basis, the application also provides a curve optimization method based on Arrhenius linear transformation and monotone decreasing error screening, and a method for constructing a kinetic model containing 11 kinds of core intermediates and 15 elementary reactions by layering according to the carbon skeleton length and tracking the main path. The method and model of the application overcome the defects of the traditional method, such as the inability to obtain the micro elementary reaction rate, the large scatter noise and the calculation disaster caused by the non-reduced dimension of the model, and provide a data source with high fidelity for the macro numerical simulation of the flow field of the energetic material detonation reaction, and have high engineering application and commercial protection value.

Description

Technical Field

[0001] This invention belongs to the field of detonation reaction analysis and calculation of energetic materials, specifically involving a method for calculating the elementary reaction rate constant in the decomposition network of hexanitrohexaazaisowulzane, as well as a curve optimization fitting method and a kinetic model construction method based on the elementary reaction rate constant. Background Technology

[0002] Hexanitrohexaazaisowulzane (CL-20, chemical formula C6H6N) 12 O 12 CL-20 is the fourth generation of energetic materials after TNT and RDX. It possesses a cage-like molecular structure, a detonation velocity of 9500 m / s, and an explosive yield three times that of TNT, making it the highest energy density conventional single-element explosive to date. As one of the highest-energy, high-density explosives known with the broadest application prospects, its thermal degradation and thermal shock evolution mechanisms under extreme conditions of high temperature, high pressure, and detonation are a core frontier in the field of energetic materials research. Accurately obtaining the reaction rate constant of CL-20 during its decomposition process is a prerequisite for conducting macroscopic numerical simulations of the detonation reaction flow field and predicting the explosive's sensitivity and thermal stability.

[0003] Traditional methods for obtaining the pyrolysis rate constant of energetic materials mainly rely on macroscopic experiments (such as thermogravimetric analysis (TGA) and differential scanning calorimetry (DSC). However, these experimental methods cannot capture ultrafast processes on the order of microseconds or even nanoseconds under extreme detonation conditions. They can only fit the "overall reaction parameters" under black-box conditions and cannot reveal the true process of microscopic elementary reactions. In recent years, molecular dynamics (MD) simulations (especially reaction molecular dynamics, such as ReaxFF-MD) have been widely used in the study of the microscopic decomposition of energetic materials. However, existing chemical evolution kinetic analysis methods have common shortcomings: they can usually only roughly count the evolution of macroscopic products (such as N2 and H2O). For the densely packed, transient macromolecules or free radical intermediates in the CL-20 decomposition reaction network, existing algorithms cannot accurately identify their entire life cycle, let alone directly calculate the microscopic rate constant of each elementary reaction in the dynamic simulation trajectory where the topology of the reaction network is changing frequently.

[0004] Furthermore, when extracting rate constants using conventional methods (such as directly fitting quantitative smoothing to the MD trajectory), a large number of "noise points" and "abnormal data points" are easily generated due to thermal fluctuations and random particle collisions. Current technologies lack an effective curve optimization mechanism based on physicochemical laws (such as Arrhenius's law) for adaptive error filtering and bad point removal of scattered microscopic reaction data. Simultaneously, since the CL-20 decomposition network contains hundreds or thousands of intermediates and elementary reactions, directly substituting the entire network into macroscopic flow field calculations would lead to a severe "numerical computation disaster." Currently, there is no systematic process in the industry that can integrate "high-fidelity calculation of microscopic elementary reaction rates," "curve error removal optimization," and "main path hierarchical dimensionality reduction model construction." Summary of the Invention

[0005] To address the shortcomings of existing technologies, such as the inability to obtain microscopic elementary reaction rates, high noise in scatter plot data, and computational disasters caused by the lack of dimensionality reduction in the model, this invention provides a method for calculating elementary reaction rate constants in the decomposition network of hexanitrohexaazaisowulzane. This invention achieves a closed-loop model across the entire network, from high-fidelity microscopic elementary constants to a macroscopic dimensionality-reduced optimized model, by introducing molecular lifetime statistics, Arrhenius monotonically decreasing error screening, and a hierarchical main path kinetics model.

[0006] To achieve the above objective, the present invention provides a method for calculating the rate constant of elementary reactions in a hexanitrohexaazaisowulzane decomposition network, comprising the following steps:

[0007] (1) The generation and annihilation times of each molecule during the decomposition reaction of hexanitrohexaazaisowurtzite at different temperatures were obtained by molecular dynamics simulation;

[0008] (2) For any molecule B in the decomposition reaction network, obtain the average lifetime τ of all molecules B with the same molecular formula at the same temperature. B ;

[0009] (3) Obtain the conditional probability ω of molecule B with the same molecular formula undergoing annihilation through each reaction pathway at the same temperature. i ;

[0010] (4) Based on the average lifespan τ B and the conditional probability ω i Calculate the elementary reaction rate constant k for the i-th decomposition reaction in which molecule B is a reactant, according to the following formula. c,i : .

[0011] In one implementation, obtaining the generation and annihilation times of each molecule during the decomposition reaction of hexanitrohexaazaisowulzane at different temperatures through molecular dynamics simulation specifically includes: assigning a fixed label to each atom; counting the product molecules twice at set time intervals during the molecular dynamics simulation and recording the set of atomic labels that make up each molecule; if the set of atomic labels that make up a molecule at the previous moment no longer forms a single molecule at the next moment but forms a union of multiple molecular label sets, then the decomposition reaction of a molecule is identified, and the annihilation time of that molecule and the generation time of its product molecules are recorded.

[0012] In one implementation, the conditional probability ω i Calculate using the following formula: , of which M B M represents the total number of chemical reactions in which molecule B participates as a reactant. i This indicates the number of times the i-th reaction occurs, in which molecule B participates as a reactant.

[0013] In one implementation, the average lifetime τ B Calculate using the following formula: , where N B α represents the quantity of molecule B produced during the reaction. j β represents the generation time of the j-th molecule B. j Let be the annihilation time of the j-th molecule B.

[0014] Those skilled in the art will understand that the pyrolysis or detonation of hexanitrohexaazaisowulzane (CL-20) under extreme conditions constitutes a highly complex multi-step cascade chemical reaction network. The initial CL-20 molecule does not transform into the final molecular fragmentation products (such as HCN and NO2) in one step. Instead, it first initiates the formation of a series of macromolecular intermediates through the breaking of specific chemical bonds (such as the N-NO2 bond). Subsequently, these intermediates undergo dozens or even hundreds of cascade elementary reactions to gradually degrade.

[0015] Based on the aforementioned multi-step cascade characteristics, the "molecule B" defined by molecular lifetime statistics in the method described in this invention is, theoretically, a universal and arbitrary designation covering any specific molecule within the entire decomposition reaction network. At the same given macroscopic temperature, due to intrinsic differences in molecular structure, chemical bond energy, and thermodynamic stability, different types of intermediate product molecules in the reaction network exhibit varying average lifetimes τ in molecular dynamics simulations. B and the conditional probability ω of evolving towards different downstream paths. i They are all different and have their own unique characteristics.

[0016] Therefore, by running the above method, it is possible to independently and accurately calculate the elementary reaction rate constant k of each specific molecule in the reaction network (specifically, molecule B includes both the initial reactant molecule of hexanitrohexaazaisowulzane itself and various intermediate product molecules or free radical fragments in the decomposition reaction network) participating in the i-th reaction pathway at a specific temperature, using "any molecule B" as the statistical unit. c,i .

[0017] In another aspect of the present invention, a method for optimizing the rate constant curve of the decomposition reaction of hexanitrohexaazaisowulzane is provided, wherein the rate constant is the rate constant of each elementary reaction calculated according to the method described above for calculating the rate constant of elementary reactions in the decomposition reaction network of hexanitrohexaazaisowulzane, and the optimization method includes the following steps:

[0018] (1) The relationship between the rate constants of each elementary reaction and temperature is transformed into lnk by variable substitution based on the Arrhenius equation. c Find the set of data points that have a linear relationship with 1 / T, and exhaustively search for all monotonically decreasing subsets of data after variable substitution;

[0019] (2) Perform linear fitting on each monotonically decreasing subset of data, calculate the relative error of each data point in the subset, and select the subset of data where the relative error of all data points is less than 2%;

[0020] (3) Select the subset with the largest number of data points from the selected subset of data as the optimized rate constant-temperature dataset, and fit the final reaction rate parameters accordingly.

[0021] In one implementation, the search for all monotonically decreasing subsets of data employs an exhaustive approach, traversing all possible combinations of subsets.

[0022] In another aspect of the present invention, a method for constructing an optimized reaction kinetic model for hexanitrohexaazaisowulzane is provided, comprising the following steps:

[0023] (1) The decomposition reaction network of hexanitrohexaazaisowurtzite was obtained through molecular dynamics simulation;

[0024] (2) The intermediate products in the reaction network are classified based on the carbon skeleton length, and the reaction network is layered according to the number of carbon atoms 6, 4, 3, and 2.

[0025] (3) Trace the main reaction pathway from hexanitrohexaazaisowurtzite to the final molecular cleavage products HCN and NO2 along the direction with the largest reaction flux, and retain the elementary reactions and corresponding intermediate products on the main pathway;

[0026] (4) Using the above analysis and calculation methods, calculate the elementary reaction rate constants of each elementary reaction on the main path at different temperatures;

[0027] (5) The elementary reaction rate constants calculated in step (4) are optimized using the above optimization method, and finally an optimized reaction kinetic model consisting of 11 intermediate products and 15 elementary reactions is constructed.

[0028] In one embodiment, the 11 intermediate products include: C6H6O 10 N 11 C6H6O8N 10 The 15 elementary reactions include: C4H4O8N8, C4H4O6N7, C4H4O4N6, C3H3O2N4, C3H3N3, C2H2O4N4, C2H2O2N3, HCN, and NO2; 12 O 12 =C6H6O 10 N 11 +NO2, C6H6N 12 O 12 =C4H4O8N8+C2H2O4N4、C6H6O 10 N 11 =C6H6O8N 10 +NO2, C6H6O 10 N 11 =C4H4O6N7+C2H2O4N4、C6H6O8N 10 =C4H4O4N6+C2H2O4N4, C4H4O8N8=C4H4O6N7+NO2, C4H4O8N8=2C2H2O4N4, C4H4O6N7=C4H 4O4N6+NO2, C4H4O6N7=C2H2O4N4+C2H2O2N3, C4H4O4N6=C3H3O2N4+NO2+HCN, C3H3O2N4= NO2+C3H3N3, C3H3N3=3HCN, C4H4O4N6= 2C2H2O2N3, C2H2O4N4= C2H2O2N3+NO2 and C2H2O2N3=2HCN+NO2.

[0029] In one implementation, the reaction pathway framework of the model is C6→C4→C3 / C2→HCN+NO2, wherein the carbon chain length gradually decreases.

[0030] In another aspect of the invention, the application of the CL-20 thermal decomposition reaction kinetic model constructed by the above-described method for calculating the elementary reaction rate constant in the CL-20 decomposition reaction network, or the method for optimizing the CL-20 decomposition reaction rate constant curve, or the method for constructing the CL-20 decomposition reaction kinetic model, to the macroscopic numerical simulation of the CL-20 explosive detonation reaction flow field is provided.

[0031] Beneficial effects

[0032] Compared with the prior art, the present invention has the following significant advantages:

[0033] 1. High-fidelity in-situ microscopic calculation of multi-step cascade reaction networks: This invention, through "atomic tag tracking" and "lifetime statistics", for the first time extracts the average lifetime and conditional probability of any individual molecule (including initial reactants and transient intermediates) in situ from the dynamically changing CL-20 decomposition MD simulation trajectory, thereby calculating the rate constant of independent elementary reactions with high fidelity.

[0034] 2. Adaptive removal of data noise based on physical and chemical laws: The curve optimization method of this invention cleverly utilizes Arrhenius linearization replacement, combined with exhaustive search of monotonically decreasing subsets and the dual checkpoint of "2% relative error", which can automatically and intelligently remove "false reactions" or "collision bad points" from MD scatter data with high noise caused by thermal fluctuations, ensuring the physical authenticity of the fitted parameters.

[0035] 3. Successfully achieved dimensionality reduction of complex networks, eliminating the curse of numerical values: The model construction method innovatively hierarchically layers the CL-20 macro-network based on carbon chain length (C6 to C2), and performs dimensionality reduction by segmenting along the main path with the highest reaction flux, ultimately reducing the dimensionality of the 10... 3 ~10 4 The network is condensed into a skeletal model of "11 core intermediates and 15 core elementary reactions". It not only fully preserves the key reaction chain of CL-20 decomposition, but also eliminates the potential "numerical curse" for numerical simulation of macroscopic detonation reaction flow fields (such as high-dimensional calculations of fluid dynamics), greatly improving simulation efficiency and accuracy. Attached Figure Description

[0036] To more clearly illustrate the technical solutions in this invention or related technologies, the accompanying drawings used in describing this invention or related technologies will be briefly introduced below. Obviously, the drawings in the following description only relate to some embodiments of the invention. For those skilled in the art, other drawings can be obtained based on these drawings without any inventive effort.

[0037] Figure 1 This is an optimized pathway diagram for the CL-20 decomposition reaction of the present invention;

[0038] Figure 2 The simulation data for the reaction C4H4O8N8=C4H4O6N7+NO2 at 2000 K are shown, where (a) is the relationship between the number of reaction occurrences and time, and (b) is the change in the number of reactant C4H4O8N8 molecules over time.

[0039] Figure 3 The histogram of molecular lifetime distribution of C4H4O8N8 molecule at 2000 K;

[0040] Figure 4 The original and optimized data are for the rate constant of the reaction HNO=H+NO;

[0041] Figure 5 The original and optimized data for the rate constant of the reaction C2H2O3N2=C2H2O2N+NO;

[0042] Figure 6 The bar chart shows the final product mass distribution obtained by simulation at 3000 K under isothermal conditions based on the original detailed model (2284 intermediate products, 6295 elementary reactions);

[0043] Figure 7 The bar chart shows the final product mass distribution obtained by simulating the simplified model of this invention (11 intermediate products and 15 elementary reactions) under isothermal conditions of 3000 K. Detailed Implementation

[0044] To better illustrate the technical means and effects of the present invention, the present invention is further described below in conjunction with non-limiting embodiments. These embodiments (including descriptions mentioned in the embodiments) are intended to illustrate implementation methods of the present invention and are not intended to limit the scope of any claims. According to the present invention, those skilled in the art will understand that many changes can be made to the specific embodiments disclosed without departing from the spirit and scope of the present invention, and the same or similar results can still be obtained.

[0045] In one embodiment of the present invention, an analytical calculation method for the rate constant of the CL-20 decomposition reaction is provided, specifically including the following technical solutions:

[0046] 1. Statistical methods for intermediate products and reaction equations in the decomposition reaction of CL-20

[0047] This invention simulates the reaction process of CL-20 at different temperatures using molecular dynamics methods, outputting atomic coordinates and determining whether chemical bonds are formed based on interatomic distances. Based on atomic coordinates, this invention establishes a method for identifying and statistically analyzing intermediate products and reaction equations in the CL-20 system. The method for identifying intermediate products involves finding the smallest set of atoms such that no atom within the set bonds with atoms outside the set. Building upon the identification of intermediate products, this invention proposes a method for identifying chemical reaction equations: assigning each atom a fixed label composed of numbers and letters. During molecular dynamics simulations, atomic motion causes positional changes, but the labels remain unchanged. Product molecules are counted twice at extremely short intervals (0.01 picoseconds), simultaneously counting the atomic labels that make up each molecule. If the set of atomic labels that made up a molecule at the previous time step no longer forms a single molecule at the next time step, but is instead the union of multiple molecular label sets, then a decomposition reaction of a molecule is identified. This method can identify all chemical reactions occurring at high temperatures in CL-20.

[0048] 2. Derivation and Calculation of the Rate Constant Calculation Method for CL-20 Decomposition Reaction Based on Molecular Lifetime

[0049] Under high-temperature conditions, the decomposition intermediates of CL-20 explosive undergo continuous decomposition and synthesis reactions. During chemical reactions, reactant molecules disappear and product molecules are generated. Based on this, the time required for a molecule to go from formation to annihilation can be defined as the molecular lifetime, and the average molecular lifetime of all molecules with the same molecular formula can be defined as the average lifetime of that molecule.

[0050] Now, let's define an extremely short time period, denoted as δt, and let P be the probability that molecule B disappears within this time period. B The number of chemical reactions that lead to the annihilation of molecule B is denoted as N, where the first to N-1 are decomposition reactions, and the last is the set of all synthesis reactions, with corresponding probabilities of occurrence in time interval δt as P1, P2, P3, ..., P N Since different chemical decompositions have different bond-breaking sites, leading to the breaking of different chemical bonds through molecular thermal vibrations, the probabilities of different decomposition reactions occurring are independent. Synthesis reactions are completed through collisions between molecule B and other molecules; this event is independent of the thermal vibrations of molecule B itself, therefore the probability of synthesis reactions occurring is independent of all decomposition reactions. Therefore, the relationship between the total probability and the sub-probabilities of each reaction path is:

[0051] (1)

[0052] Among them, P BLet P represent the probability that molecule B transforms into other molecules through a chemical reaction within a very short time interval δt. Molecule B can be transformed into different molecules through different chemical reaction pathways within this time interval; the probability of transformation through the i-th reaction pathway is denoted as P. i .

[0053] When the time interval δt is sufficiently short, the subprobabilities of each reaction path are much less than 1, therefore a first-order approximation can be applied. The formula for calculating the total probability under the first-order approximation is:

[0054] (2)

[0055] Meanwhile, the probability of a molecule annihilating through the i-th path is defined as ω. i The calculation formula is:

[0056] (3)

[0057] Where, ω i This represents the probability that molecule B will annihilate through the i-th reaction pathway after a chemical reaction has already occurred.

[0058] The lifetime of molecule B is denoted as τ. B The formula for calculating molecular lifetime using probability is:

[0059] (4)

[0060] The rate constant for the i-th decomposition reaction is denoted as k. c,i The calculation method is as follows:

[0061] (5)

[0062] The molecular lifetime and the probability of the i-th reaction occurring during molecular annihilation in the formula for calculating the reaction rate constant can both be obtained through statistical tracking of intermediate product molecules and their atomic tag sets. Using this method, large-scale statistical analysis can yield the rate constants for all decomposition reactions of CL-20 at different temperatures.

[0063] The analytical calculation method based on the rate constant of the decomposition reaction of hexanitrohexaazaisowoodsane mainly includes the following two aspects:

[0064] 1. Mean lifespan (τ) B ) calculation

[0065] According to the method proposed in this invention, a reaction rate parameter library can be constructed by statistically screening the product molecules and reaction equations throughout the entire detonation process. Applying the statistical method for intermediate products and reaction equations of the CL-20 decomposition reaction established in this invention, the generation and annihilation times of each product molecule in the entire reaction process can be determined by atomic tags and the number of reaction equations, as well as the average lifetime (τ) of molecule B. B The calculation method is as follows:

[0066] (6)

[0067] Where, N B α represents the quantity of molecule B produced through a chemical reaction at a given temperature. j β represents the generation time of the j-th molecule B. j Let be the annihilation time of the j-th molecule B.

[0068] 2. Probability (ω) i ) calculation

[0069] Molecule B may have multiple decomposition and synthesis reaction pathways. Let M be the total number of reactions in which molecule B participates. B The number of times the i-th reaction occurs is denoted as M. i Therefore, the conditional probability of the i-th reaction (ω) i )for:

[0070] (7)

[0071] Among them, M B M represents the total number of chemical reactions in which molecule B participates as a reactant. i This indicates the number of times the i-th reaction occurs, in which molecule B participates as a reactant.

[0072] According to the method proposed in this invention, the reaction rate constant can be calculated from the molecular lifetime and the probability of the reaction occurring.

[0073] In another embodiment of the present invention, a method for optimizing the rate constant curve of the CL-20 decomposition reaction is provided, specifically including the following technical solutions:

[0074] Molecular lifetime methods can be used to obtain the rate constants of thermal decomposition reactions at different temperatures, forming rate constant-temperature relationship curves. However, the molecular dynamics simulation system for reaction force fields is limited (1×10⁻⁶). 5 (Atomic scale), the rate constant-temperature curve obtained by this method has fluctuation errors and bad points (step data points caused by fluctuation errors), which need to be eliminated before parameter fitting.

[0075] This invention proposes to eliminate bad pixels using an exhaustive approach. First, variable substitution is performed based on the Arrhenius equation:

[0076] (8)

[0077] Where k c Let T be the rate constant, R be the gas constant (8.314 J / mol / K), A be the pre-exponential factor, and E be the gas constant. a is the activation energy of the reaction. The latter two are the reaction rate constants that need to be calculated using fitting.

[0078] Taking the logarithm of both sides of the Arrhenius equation:

[0079] (9)

[0080] Therefore, ln k c The ideal relationship with 1 / T is a linear decreasing relationship. The rate constant-temperature curve obtained from molecular lifetime calculations is converted to ln k. c Given a data point set of 1 / T, search for all monotonically decreasing subsets of data. With an original set of 30 data points, up to four million monotonically decreasing subsets of data can be found.

[0081] A linear fit is performed on each subset, and the relative error at each point is calculated after the fit is complete. The formula for calculating the relative error at the i-th point is:

[0082] (10)

[0083] Among them, A and E a To fit the pre-exponential factor and activation energy of the reaction obtained from this subset of data, {ln k c,i} represents a subset of the set consisting of reaction rate constants at different temperatures.

[0084] A subset of data points with relative errors less than 2% was selected, and then the dataset with the longest data chain within this subset was chosen as the final optimized rate constant-temperature dataset. This dataset was then used to fit the final pre-exponential factor A and activation energy E of the thermal decomposition reaction. a .

[0085] In another embodiment of the present invention, a method for constructing an optimization model for the decomposition reaction of CL-20 is also provided, specifically including the following technical solutions:

[0086] This invention proposes an analytical method for analyzing the decomposition reaction mechanism of CL-20 based on its carbon architecture, optimizing the reaction pathway by statistically analyzing changes in molecular length based on carbon chains. The final optimized reaction pathway framework is: C6→C4→C3 / C2→HCN+NO2. Figure 1 As shown, Figure 1 This diagram shows the optimized decomposition reaction pathway of CL-20 according to the present invention. It illustrates the process of CL-20 gradually decomposing from its initial state into HCN and NO2 through various intermediate products, as well as the small molecule products generated at each step of the reaction. It visually demonstrates the decreasing carbon skeleton decomposition rule of C6→C4→C3 / C2→HCN+NO2. The transformation relationships of each substance in the diagram correspond one-to-one with the 15 elementary reactions recorded in Table 1.

[0087] CL-20 molecular dynamics simulations yielded 2284 intermediate products and 6295 corresponding elementary reaction equations. The optimized model's reaction pathway contained 11 intermediate products and 15 elementary reactions. The rate equations for the elementary reactions were Arrhenius equations, with rate parameters including the pre-exponential factor A and the activation energy E. a Two items. The reaction equations and rate parameters obtained by this invention are shown in Table 1, and the specific optimized reaction pathway is shown in the attached figure. After thermal decomposition, the intermediate product formed by CL-20 is a mixture of HCN and NO2, and the subsequent combustion process can be described by a general reaction model.

[0088] Table 1. Optimized CL-20 decomposition reaction rate parameters

[0089]

[0090] The optimized model of the CL-20 decomposition reaction consists of 11 intermediate products and 15 elementary reactions, and can be used for macroscopic numerical simulation of the flow field in explosive detonation reactions. Based on the reaction model, a multiphase flow macroscopic numerical simulation program was developed using C++. The simulation algorithm for solid explosives uses the material point method, and the simulation algorithm for gaseous products uses the implicit Eulerian method. The numerical simulation employs an adaptive time step with a time step interval of 10. -6 Up to 10 -8 s, grid resolution is 10 -3 m.

[0091] In multiphase flow simulations, each product has a different thermodynamic state, and the chemical reaction generates different material fields, forming a reaction flow field. According to the energy and momentum conservation relationships of molecular collisions, the distribution of the heat of reaction after the reaction is completed is inversely proportional to the reciprocal of the molecular mass of the product.

[0092] Based on this decomposition reaction model and fluid-structure interaction numerical method, combined with a general combustion reaction database, a systematic macroscopic numerical simulation of the reaction process of CL-20 explosive under various typical working conditions was conducted, obtaining its reaction intensity distribution and product dynamic evolution law. The simulation conditions covered typical stimulus loading scenarios such as high-speed impact and hot-trigger ignition. The impact condition involved an impact velocity of 500–1000 m / s and an impact surface of concrete. The hot-trigger condition involved an ignition temperature of 600–1000 K and an ignition zone radius of 0.01–0.03 m. The numerical simulation reproduced the complete evolution process of the explosive from thermal decomposition, combustion, explosion, and even detonation under different forces and thermal loads, quantitatively characterizing the changes in temperature, pressure flow field, and product composition over time. Based on the above application method, the reaction intensity level of CL-20 explosive under different force and thermal loading environments can be predicted, and its chemical explosion safety in complex multi-condition environments can be systematically evaluated.

[0093] Example

[0094] Example 1: Calculation of the rate constant of elementary reactions in the CL-20 decomposition network

[0095] 1. Simulation parameter settings

[0096] Molecular dynamics simulations of an ε-CL-20 supercell (84,672 atoms) were performed using the ReaxFF-lg reactive force field. Thirty-five temperature points (600K–4000K, 100K intervals) were investigated. The ensemble was NVT, the thermobath algorithm was Nosé-Hoover, the integration step size was 0.1 fs, the trajectory was output every 0.01 ps, and the total simulation time for each temperature point was 100 ps.

[0097] 2. Atomic Tag Tracking and Molecular Recognition

[0098] Each atom is assigned a fixed label (consisting of numbers and letters). In each output frame, chemical bonds are determined based on interatomic distances (ReaxFF bond order threshold), and a depth-search algorithm is used to identify connected components, with each connected component corresponding to a molecule. The set of atom labels for each molecule is recorded.

[0099] 3. Record of molecular creation and annihilation events

[0100] Compare the lists of molecules in adjacent frames:

[0101] If a certain set of atomic tags appears in the current frame, and that set split into multiple disconnected subsets in the previous frame, then record that the molecule was generated in the current frame.

[0102] If a set of atomic tags in the current frame splits into multiple subsets, then the annihilation of that molecule in the current frame is recorded.

[0103] 4. Mean Life Calculation

[0104] Taking the reaction C4H4O8N8 = C4H4O6N7 + NO2 at 2000 K as an example, during the entire simulation process, a total of 1344 C4H4O8N8 molecules were generated and annihilated (N... B =1344). Lifetime per molecule (β) j -α j )Statistical distribution as Figure 3 As shown. Calculate the average lifespan according to formula (6) in the instruction manual:

[0105] .

[0106] 5. Statistical analysis of reaction conditional probability

[0107] The total number of times this molecule participates as a reactant in chemical reactions across all annihilation events is M. B =1344 (each annihilation event corresponds to one reaction). Among them, M is the number of times the target decomposition path (generating C4H4O6N7+NO2) is traversed. i =84 times. Calculate the conditional probability according to formula (7) in the instruction manual: .

[0108] 6. Rate constant calculation: Calculate according to formula (5) in the instruction manual:

[0109] .

[0110] 7. Comparison with traditional event statistics methods

[0111] Traditional event statistics only utilize time-step data of the reaction. In 100 ps, ​​with 10,000 time steps, this reaction occurs in only 84 time steps (see...). Figure 2 The effective data volume was 84 molecules. The calculated relative standard deviation (RSD) of the rate constant was 85.0%. However, this invention utilizes lifetime data from all 1344 molecules, with an RSD of only 18.5%, increasing the effective data volume by 16 times (>10 times) and significantly improving accuracy. The comparison results are shown in Table 2.

[0112] Table 2. Comparison of effective data volume and accuracy of different methods at 2000 K

[0113]

[0114] Example 2: Optimization of the reaction rate constant curve

[0115] 1. Input data

[0116] Taking the reactions HNO=H+NO and C2H2O3N2=C2H2O2N+NO as examples, the original rate constants {ln k} at 35 temperature points (600~4000K) were obtained using the method in Example 1. c,i}

[0117] 2. Arrhenius Transformation

[0118] According to formula (9) in the instruction manual, each data point (T, k) c The transformation is (1 / T, ln k) into (1 / T, ln k). c ), to obtain a dataset with linear relationships.

[0119] 3. Exhaustive search for monotonically decreasing subsets

[0120] Write a program (such as a Python recursive backtracking algorithm) to generate all possible subsets of the transformed dataset and retain only those subsets that satisfy the monotonically decreasing condition. For 35 points, theoretically there are at most C(35, k) combinations, but traversal can be performed efficiently through pruning.

[0121] 4. Linear Fitting and Relative Error Calculation

[0122] For each subset, perform a least-squares linear fit to obtain the fitting parameters ln A and -E. a / R. Calculate the relative error of each original data point in the subset according to formula (10) in the instruction manual (i.e., the absolute value of the difference between the fitted value and the original value divided by the original value, and take the percentage). Select the subset with a relative error of ≤2% for the data points.

[0123] 5. Selecting the optimal subset and the final fit

[0124] Within the selected subset, the subset with the most valid data points is chosen. For the reaction HNO=H+NO, the optimal subset contains 14 data points (out of a total of 24, excluding 10 bad points); for the reaction C2H2O3N2=C2H2O2N+NO, the optimal subset contains 8 data points (excluding 16 bad points). Linear regression is then performed again on the optimal subset to obtain the final pre-exponential factor A and activation energy E. a The optimized curve is compared with the original data. Figure 4 and Figure 5 It can be seen that the bad points were effectively eliminated, and the fitted curve conforms to the expectations of the Arrhenius theory.

[0125] Example 3: Kinetic Model of CL-20 Decomposition Reaction

[0126] 1. Obtaining the original reaction network

[0127] Using the method described in Example 1, the entire decomposition reaction process of CL-20 at 35 temperature points was analyzed, identifying all 2284 intermediate products and 6295 elementary reactions. The Arrhenius parameters (pre-exponential factor A and activation energy E) for each elementary reaction were calculated. a ).

[0128] 2. Simplify model selection

[0129] (1) Carbon skeleton length analysis: The number of carbon atoms in each intermediate product was counted, and the reaction network was divided into layers according to the number of carbon atoms 6, 4, 3, and 2.

[0130] (2) Main path tracing: from CL-20 (molecular formula C6H6N) 12 O 12 Starting from [the point], the main path is traced along the direction with the highest reaction flux (i.e., the product of the rate constant of each branch reaction and the reactant concentration) to the final products HCN and NO2. This main path is as follows:

[0131] Step 1: CL-20 = C6H6O 10 N 11 +NO2

[0132] Step 2: CL-20 = C4H4O8N8 + C2H2O4N4

[0133] Step 3: C4H4O8N8 = C4H4O6N7 + NO2

[0134] Step 4: C4H4O6N7 = C4H4O4N6 + NO2

[0135] Step 5: C4H4O4N6 = C3H3O2N4 + NO2 + HCN

[0136] Step 6: C3H3O2N4 = NO2 + C3H3N3

[0137] Step 7: C3H3N3 = 3HCN

[0138] And the parallel path: C4H4O8N8→2C2H2O4N4→2C2H2O2N3+2NO2→4HCN+2NO2

[0139] It also includes the reaction from C4H4O6N7 directly to C2H2O4N4+ C2H2O2N3, and the subsequent decomposition of C2H2O4N4.

[0140] This main pathway conforms to the framework of gradually decreasing carbon chain length: C6→C4→C3 / C2→HCN+NO2.

[0141] (3) Retain the core reaction: retain all the elementary reactions (15 in total) on the main path and delete the other branch reactions. These 15 reactions are completely consistent with those listed in step (4) of claim 7.

[0142] (4) Parameter verification: Using the rate constant data of the corresponding reactions in the detailed reaction network (obtained from Example 1), the Arrhenius parameters (A and E) were refitted for these 15 elementary reactions. a The fitting method was least squares, and the fitting temperature range was 600~4000K. The final parameters are shown in Table 1.

[0143] 3. Constant Temperature Verification

[0144] (1) Verification conditions: The constant temperature was set at 3000K (close to the detonation temperature of CL-20), the time step was 1ns, the total number of steps was 20,000, and the corresponding chemical reaction time was 2×10 -5 s is sufficient to complete the reaction.

[0145] (2) Comparison objects:

[0146] Detailed model: 2284 intermediate products + 6295 elementary reactions (obtained from the original ReaxFF MD statistics).

[0147] Simplified model: 11 intermediate products + 15 elementary reactions (this invention).

[0148] (3) Combustion reaction model:

[0149] The combustion reaction model uses the general combustion reaction database GRI-Mech. After molecular fragmentation, the general combustion reaction model is used to simulate subsequent evolution.

[0150] (4) Result Comparison: The final product mass distributions obtained by the two models are shown in the figure below. Figure 6 and Figure 7 The mass fractions of the main products CO2 and N2 are compared in Table 3. The relative errors of CO2 and N2 between the simplified model and the detailed model are 8.06% and 2.62%, respectively, both within acceptable ranges, demonstrating the effectiveness of the simplified model.

[0151] Table 3. Comparison of main product mass fractions between simplified and detailed models at 3000 K

[0152]

[0153] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the various embodiments of the present invention.

Claims

1. A method for calculating the rate constant of elementary reactions in a decomposition reaction network of hexanitrohexaazaisowulzane, characterized in that, Includes the following steps: (1) The generation and annihilation times of each molecule during the decomposition reaction of hexanitrohexaazaisowoodsane at different temperatures were obtained by molecular dynamics simulation; (2) For any molecule B in the decomposition reaction network, obtain the average lifetime τ of all molecules B with the same molecular formula at the same temperature. B ; (3) Obtain the conditional probability ω of molecule B with the same molecular formula undergoing annihilation through each reaction pathway at the same temperature. i ; (4) Based on the average lifespan τ B and the conditional probability ω i Calculate the elementary reaction rate constant k for the i-th decomposition reaction in which molecule B is a reactant, according to the following formula. c,i : 。 2. The method according to claim 1, characterized in that, The generation and annihilation times of each molecule during the decomposition reaction of hexanitrohexaazaisowulzane at different temperatures, obtained through molecular dynamics simulations, specifically include: Each atom is assigned a fixed label. During the molecular dynamics simulation, the product molecules are counted twice at set time intervals, and the set of atomic labels that make up each molecule is recorded. If the set of atomic labels that make up a molecule at the previous time does not form a single molecule at the next time but forms a union of multiple molecular label sets, then the decomposition reaction of a molecule is identified, and the annihilation time of the molecule and the generation time of its product molecules are recorded.

3. The method according to claim 1, characterized in that, The conditional probability ω i Calculate using the following formula: Among them, M B M represents the total number of chemical reactions in which molecule B participates as a reactant. i This indicates the number of times the i-th reaction occurs, in which molecule B participates as a reactant.

4. The method according to claim 1, characterized in that, The average lifespan τ B Calculate using the following formula: Where, N B α represents the quantity of molecule B produced during the reaction. j β represents the generation time of the j-th molecule B. j Let be the annihilation time of the j-th molecule B.

5. A method for optimizing the rate constant curve of the decomposition reaction of hexanitrohexaazaisowulzane, characterized in that, The rate constants are the elementary reaction rate constants calculated according to the method described in claim 1, and the optimization method includes the following steps: (1) The relationship between the rate constants of each elementary reaction and temperature is transformed into lnk by variable substitution based on the Arrhenius equation. c Find the set of data points that have a linear relationship with 1 / T, and exhaustively search for all monotonically decreasing subsets of data after variable substitution; (2) Perform linear fitting on each monotonically decreasing subset of data, calculate the relative error of each data point in the subset, and select the subset of data where the relative error of all data points is less than 2%; (3) Select the subset with the largest number of data points from the selected subset of data as the optimized rate constant-temperature dataset, and fit the final reaction rate parameter accordingly.

6. The method according to claim 5, characterized in that, The search for all monotonically decreasing subsets of data employs an exhaustive approach, traversing all possible subset combinations.

7. A method for constructing a kinetic model of the decomposition reaction of hexanitrohexaazaisowulzane, characterized in that, The method includes the following steps: (1) The decomposition reaction network of hexanitrohexaazaisowurtzite was obtained through molecular dynamics simulation; (2) The intermediate products in the reaction network are classified based on the carbon skeleton length, and the reaction network is layered according to the number of carbon atoms 6, 4, 3, and 2. (3) Trace the main reaction pathway from hexanitrohexaazaisowurtzite to the final molecular cleavage products HCN and NO2 along the direction with the largest reaction flux, and retain the elementary reactions and corresponding intermediate products on the main pathway; (4) According to the method described in claim 1, calculate the elementary reaction rate constants of each elementary reaction on the main path at different temperatures; (5) The elementary reaction rate constants calculated in step (4) are optimized by the optimization method according to claim 5, and finally an optimized reaction kinetic model consisting of 11 intermediate products and 15 elementary reactions is constructed.

8. The method according to claim 7, characterized in that, The 11 intermediate products include: C6H6O 10 N 11 C6H6O8N 10 , C4H4O8N8, C4H4O6N7, C4H4O4N6, C3H3O2N4, C3H3N3, C2H2O4N4, C2H2O2N3, HCN and NO2; Comprehensive of the 15 types of bases listed: C6H6N 12 O 12 =C6H6O 10 N 11 +NO2, C6H6N 12 O 12 =C4H4O8N8+C2H2O4N4,C6H6O 10 N 11 =C6H6O8N 10 +NO2, C6H6O 10 N 11 =C4H4O6N7+C2H2O4N4,C6H6O8N 10 =C4H4O4N6+C2H2O4N4, C4H4O8N8=C4H4O6N7+NO2, C4H4O8N8=2C2H2O4N4, C4H4O6N7=C4H4O4N6+NO2, C4H4O6N7=C2H2O4N4+C2H2O2N3, C4H4O4N6=C3H3O2N4+NO2+HCN, C3H3O2N4=NO2+C3H3N3, C3H3N3=3HCN, C4H4O4N6=2C2H2O2N3, C2H2O4N4=C2H2O2N3+NO2 sum C2H2O2N3=2HCN+NO2.

9. The method according to claim 7, characterized in that, The reaction pathway framework of the model is C6→C4→C3 / C2→HCN+NO2, in which the carbon chain length gradually decreases.

10. The application of the kinetic model of the thermal decomposition reaction of hexanitrohexaazaisowurtzite constructed by the method of claim 1, or the optimization method of claim 5, or the method of claim 7, in the macroscopic numerical simulation of the detonation reaction flow field of hexanitrohexaazaisowurtzite explosive.

Images