An Optimization Method for Determining Parameters of Condensed Explosive Reaction Rate Model

Abstract

The present invention relates to an optimization method for determining the parameters of a condensed state explosive reaction rate model, belonging to the field of energetic materials. The method comprises: obtaining and preprocessing a condensed state explosive shock detonation experimental data set; establishing a condensed state explosive reaction rate model and embedding it into nonlinear dynamic analysis software, establishing a calculation model within the software; establishing a time cost objective function and a physical profile cost objective function; performing finite element calculations on the reaction rate model parameters using the software and the calculation model to obtain numerical simulation results, comparing the numerical simulation results with the shock detonation experimental data to obtain the values ​​of the time cost objective function and the physical profile cost objective function; and finding the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold based on a genetic algorithm and a NASA II algorithm. The method achieves automatic optimization of the reaction rate model parameters by obtaining the corresponding reaction rate model parameters when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold.

Description

Technical Field

[0001] The invention relates to the technical field of energetic materials, in particular to an optimization method for determining parameters of a condensed explosive reaction rate model. Background Art

[0002] Condensed matter explosives are the primary charge in mainstream equipment, and high-fidelity reaction flow simulation of their shock detonation behavior is crucial for improving equipment safety. The reaction rate model, based on the ZND model theory, simulates the shock detonation behavior of explosives using the equation of state for the unreacted explosive and the equation of state for the detonation products. This model is the mainstream model for studying shock detonation of condensed matter explosives. After obtaining the equations of state for the unreacted explosive and the detonation products, the reaction rate model parameters must be compared with a one-dimensional shock detonation experimental dataset.

[0003] However, reaction rate model parameters often have no direct correlation with shock initiation experimental datasets. Currently, the determination of these parameters relies on manual trial and error, which not only consumes a lot of time and effort but also has a high degree of subjectivity. It is impossible to objectively evaluate whether the parameters have reached the optimal level, and it is impossible to objectively compare the performance of different types of reaction rate models on the same set of shock initiation experimental datasets.

[0004] Therefore, it is necessary to provide a method for automatically determining the parameters of the reaction rate model of condensed explosives, so as to improve the efficiency and objectivity of determining the parameters of the reaction rate model of condensed explosives. Summary of the Invention

[0005] In view of the above analysis, the present invention aims to provide an optimization method for determining the reaction rate model parameters of condensed state explosives, so as to solve the problem that the current determination of the reaction rate model parameters of condensed state explosives relies on manual labor, which is not only time-consuming and labor-intensive but also highly subjective.

[0006] The present invention provides an optimization method for determining parameters of a condensed explosive reaction rate model, the method comprising the following steps:

[0007] Acquiring a condensed explosive shock initiation experimental data set, and preprocessing the shock initiation experimental data set;

[0008] Establishing a condensed explosive reaction rate model and embedding it into nonlinear dynamic analysis software, and establishing a calculation model consistent with the impact detonation experimental conditions in the nonlinear dynamic analysis software;

[0009] Establish a time cost objective function and a physical profile cost objective function; perform finite element calculation on the reaction rate model parameters using the nonlinear dynamic analysis software and the calculation model to obtain numerical simulation results, and obtain the values ​​of the time cost objective function and the physical profile cost objective function by comparing the numerical simulation results with the corresponding data in the impact initiation experiment data set; based on the genetic algorithm and the NASAⅡ algorithm, find the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold value as the reaction rate model parameters, thereby optimizing the reaction rate model parameters.

[0010] Furthermore, the preprocessing of the shock initiation experiment data set includes:

[0011] Filtering and cleaning the shock initiation experiment data set;

[0012] Based on the cleaned shock initiation experimental data set, characteristic reference points at different Lagrangian positions are extracted;

[0013] The characteristic reference points represent data points in the shock initiation experiment data set corresponding to each shock initiation experiment that need to be compared with the numerical simulation results of the calculation model; the characteristic reference points include physical flow field information at different Lagrangian positions;

[0014] The physical flow field information includes data on the impact detonation pressure changing with time and data on the velocity changing with time.

[0015] Furthermore, the expression of the time cost objective function is:

[0016] , (2)

[0017] in, Indicates the total number of shock initiation experiments that need to be compared, represents the total number of Lagrangian positions that need to be compared in one experiment, represents the wave arrival time at the jth Lagrangian position of the i-th experiment in the numerical simulation results obtained by finite element calculation, represents the wave arrival time at the jth Lagrangian position of the i-th experiment in the shock initiation experiment dataset.

[0018] Furthermore, the expression of the physical contour cost objective function is:

[0019] , (3)

[0020] in, Indicates the total number of shock initiation experiments that need to be compared, represents the total number of Lagrangian positions that need to be compared in one experiment, represents the total number of characteristic reference points at each Lagrangian position, It represents the physical flow field information of the kth characteristic reference point at the jth Lagrangian position of the i-th experiment in the numerical simulation results obtained by finite element calculation. Represents the physical flow field information of the kth characteristic reference point at the jth Lagrangian position of the i-th experiment in the shock initiation experiment dataset.

[0021] Furthermore, the method of finding the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold value as the reaction rate model parameter based on the genetic algorithm and the NASA II algorithm includes:

[0022] Step S31, initializing configuration parameters, generating genetic individuals based on the configuration parameters to construct an initial population; the genetic individuals contain the reaction rate model parameter information;

[0023] Step S32: Calculate the values ​​of the time cost objective function and the physical outline cost objective function corresponding to each individual in the population; select excellent individuals to form the parent population through fast non-dominated sorting and crowding calculation based on the values ​​of the time cost objective function and the physical outline cost objective function corresponding to each individual;

[0024] Step S33, determine whether the termination condition is met. If so, output the optimal solution when the values ​​of the time cost objective function and the physical contour cost objective function corresponding to the individuals in the current population are both less than the threshold value as the reaction rate model parameter; otherwise, adaptively expand the value range of the parameters contained in the individuals in the parent population as needed, and cross and mutate each individual in the parent population according to the crossover probability and mutation probability to obtain a new genetic individual, obtain the offspring population based on the new genetic individual, merge the offspring population and the parent population as the population of the next iteration, and repeat steps S32 to S33.

[0025] Furthermore, the selecting excellent individuals to form the parent population by fast non-dominated sorting and crowding calculation includes:

[0026] Perform a fast non-dominated sort on all individuals in the population and stratify all individuals according to the Pareto frontier;

[0027] For the initial population, individuals are selected from the stratified individuals through a tournament selection algorithm to form a parent population;

[0028] For non-initial populations, individuals corresponding to the Pareto frontiers of the first n layers are selected from the stratified individuals so that their sum is greater than or equal to the number of individuals in the initial population; if the sum of the individuals corresponding to the Pareto frontiers of the first n layers is equal to the number of individuals in the initial population, the individuals corresponding to the Pareto frontiers of the first n layers are used as the parent population; otherwise, the crowding degree of all individuals corresponding to the Pareto frontier of the nth layer is calculated and sorted from large to small, and the first m individuals are selected so that the sum of the individuals corresponding to the Pareto frontier of the first n-1 layers plus m is equal to the number of individuals in the initial population, and the individuals corresponding to the Pareto frontier of the first n-1 layers and the m individuals selected from the Pareto frontier of the nth layer are used as the parent population.

[0029] Furthermore, the fast non-dominated sorting refers to stratifying all individuals in the population by non-dominated relationships, including the following steps:

[0030] Step S301: Pair all individuals in the population and calculate the non-dominated relationships between them. Obtain the set of individuals dominated by each individual and the set of individuals dominating the individual. Select individuals whose set of individuals dominating the individual is 0 to form the first-level Pareto frontier. These individuals are the best solutions in the population.

[0031] Step S302: Determine whether the number of individuals remaining in the population is 0. If so, the stratification is completed. Otherwise, all individuals remaining in the population are paired and the non-dominated relationship between them is calculated to obtain the set of individuals dominated by each individual and the set of individuals that dominate the individual. The individuals whose set of individuals that dominate the individual is 0 are selected to form the next layer of Pareto frontier, and step S302 is repeated.

[0032] Furthermore, the values ​​of the time cost objective function corresponding to the individuals in the set are sorted from small to large. For the individuals corresponding to the maximum and minimum values ​​of the time cost objective function, the congestion corresponding to their time cost objective function is set to infinity; for other individuals, the congestion corresponding to their time cost objective function is calculated using the following formula:

[0033] , (4)

[0034] in, is the number of individuals in the set, is the maximum value of the time cost objective function corresponding to the individual in the set, is the minimum value of the time cost objective function corresponding to the individual in the set, Indicates the The value of the time cost objective function corresponding to each individual;

[0035] Sort the values ​​of the physical contour cost objective function corresponding to the individuals in the set from small to large. For the individuals corresponding to the maximum and minimum values ​​of the physical contour cost objective function, set the crowding degree corresponding to their physical contour cost objective function to infinity. For other individuals, calculate the crowding degree corresponding to their physical contour cost objective function using the following formula:

[0036] , (5)

[0037] in, is the maximum value of the physical contour cost objective function corresponding to the individual in the set, is the minimum value of the physical contour cost objective function corresponding to the individual in the set, Indicates the The value of the physical contour cost objective function corresponding to each individual;

[0038] For each individual, the congestion degree corresponding to its time cost objective function and physical contour cost objective function is accumulated to obtain the congestion degree of the individual.

[0039] Furthermore, the adaptive expansion of the value range of the parameter included in the individuals in the parent population as needed includes:

[0040] Determine whether the value range of each parameter contained in the individuals in the parent population needs to be adaptively expanded; the adaptive expansion includes upper boundary adaptive expansion and lower boundary adaptive expansion;

[0041] If the parameter needs to be adaptively expanded on the upper boundary, the upper boundary of the parameter is multiplied by its upper boundary scaling factor as its upper boundary;

[0042] If the parameter needs to be adaptively expanded at the lower boundary, the lower boundary of the parameter is multiplied by its lower boundary scaling factor as its lower boundary.

[0043] Furthermore, in the nonlinear dynamic analysis software, the same material properties, geometric model, contact conditions, load conditions and boundary conditions as those in the shock initiation experiment are set, thereby establishing a calculation model consistent with the shock initiation experiment working conditions.

[0044] Compared with the prior art, the present invention can achieve at least one of the following beneficial effects:

[0045] 1. The present invention establishes a time cost objective function and a physical profile cost objective function, and searches for the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold value based on a genetic algorithm and a NASA II algorithm as the reaction rate model parameters, thereby realizing automatic optimization of the reaction rate model parameters. This solves the problem that the current determination of the reaction rate model parameters of condensed explosives relies on manual labor, which is not only time-consuming and labor-intensive but also highly subjective.

[0046] 2. The present invention adaptively expands the value range of the parameters contained in the individuals in the parent population as needed, thereby improving the calculation efficiency while obtaining a reasonable parameter value range.

[0047] 3. The present invention is applicable to the optimization of various reaction rate model parameters and has particular versatility. It can also compare the adaptability of different reaction rate models on the same shock initiation experimental data set, thereby evaluating the reaction rate models.

[0048] In the present invention, the above-mentioned technical solutions can be combined with each other to achieve more preferred combinations. Other features and advantages of the present invention will be described in the following description, and some advantages will become apparent from the description or be learned through practice of the present invention. The objectives and other advantages of the present invention can be realized and obtained through the contents particularly pointed out in the description and drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] The accompanying drawings are only used for the purpose of illustrating specific embodiments and are not to be considered as limiting the present invention. Throughout the drawings, the same reference symbols denote the same components.

[0050] Figure 1 A flow chart of an optimization method for determining reaction rate model parameters of condensed explosives according to an embodiment of the present invention;

[0051] Figure 2 This is a one-dimensional Lagrangian analysis experimental test device for impact initiation according to an embodiment of the present invention;

[0052] Figure 3 The physical flow field information at different Lagrangian positions in the shock initiation experiment data set corresponding to the shock initiation experiment of the first embodiment of the present invention;

[0053] Figure 4 The optimization results of the embodiment of the present invention based on the method of the present invention;

[0054] Figure 5 This embodiment of the present invention is based on the 75% optimization result of James R. Gambino. DETAILED DESCRIPTION

[0055] The preferred embodiments of the present invention will be described in detail below in conjunction with the accompanying drawings, wherein the accompanying drawings constitute a part of this application and are used together with the embodiments of the present invention to illustrate the principles of the present invention, and are not used to limit the scope of the present invention.

[0056] A specific embodiment of the present invention discloses an optimization method for determining the parameters of a condensed explosive reaction rate model. Figure 1 As shown, the method includes the following steps:

[0057] Step S1, obtaining a condensed explosive shock initiation experimental data set, and preprocessing the shock initiation experimental data set;

[0058] Step S2: establishing a condensed explosive reaction rate model and embedding it into nonlinear dynamic analysis software, and establishing a calculation model consistent with the impact detonation experimental conditions in the nonlinear dynamic analysis software;

[0059] Step S3, establishing a time cost objective function and a physical profile cost objective function; performing finite element calculation on the reaction rate model parameters using the nonlinear dynamic analysis software and the calculation model to obtain numerical simulation results, and obtaining values ​​of the time cost objective function and the physical profile cost objective function by comparing the numerical simulation results with corresponding data in the impact initiation experiment data set; based on the genetic algorithm and the NASAⅡ algorithm, finding the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold value as the reaction rate model parameters, thereby optimizing the reaction rate model parameters.

[0060] Specifically, in step S1, obtaining a condensed explosive shock initiation experimental data set includes:

[0061] Conduct shock detonation experiments on condensed explosives;

[0062] Obtain the experimental data set of shock initiation of condensed explosives.

[0063] Specifically, a certain condensed matter explosive is selected, and a one-dimensional Lagrangian analysis experiment of shock detonation is carried out on the condensed matter explosive using a testing device to obtain physical flow field information of the condensed matter explosive at different Lagrangian positions under different shock intensities, so as to construct the shock detonation experimental data set; the physical flow field information includes data on the change of shock detonation pressure over time and data on the change of velocity over time, which can be expressed as a pressure-time curve and a velocity-time curve.

[0064] It should be noted that the physical flow field information of condensed matter explosives at different Lagrangian positions under a single shock intensity represents the data from a single shock initiation experiment; data from multiple shock initiation experiments constitute the shock initiation experimental dataset. One-dimensional Lagrangian analysis experiments on shock initiation of condensed matter explosives can be conducted under different initial states, including different initial densities and temperatures.

[0065] For example, based on the manganese copper piezoresistive pressure measurement technology, the explosive plane lens detonation loading technology and the air and partition attenuation technology, a one-dimensional Lagrangian analysis experimental test device for the shock initiation of condensed explosives is constructed, such as Figure 2 As shown. During the experiment, a detonator detonates a planar lens to generate a shock wave. A trigger probe is turned on, causing a pulsed constant current source to begin supplying power to the manganese-copper piezoresistive sensor. The shock wave is then waveform-adjusted by the condensed explosive planar lens to form a planar detonation wave, which detonates the loading charge. The planar detonation wave is then attenuated by the air gap and aluminum partition before loading the condensed explosive under test. An oscilloscope is used to record the time-varying voltage of the manganese-copper piezoresistive sensor at four different Lagrangian positions within the condensed explosive under test. Combined with the pre-calibrated piezoresistance relationship of the manganese-copper piezoresistive sensor, data on the time-varying impact initiation pressure at these four different Lagrangian positions is obtained. By varying the height of the air gap ring or the thickness of the aluminum partition, the loading pressure acting on the explosive under test, i.e., the impact intensity, can be controlled. This allows data on the time-varying impact initiation pressure at different Lagrangian positions of the condensed explosive under different impact intensities to be obtained, thereby constructing the shock initiation experimental data set.

[0066] Furthermore, the preprocessing of the shock initiation experiment data set includes:

[0067] Filtering and cleaning the shock initiation experiment data set;

[0068] Based on the cleaned shock initiation experimental dataset, characteristic reference points at different Lagrangian positions are extracted.

[0069] Specifically, the shock detonation experimental dataset was filtered using Python's Savitzky-Golay function to eliminate noise. The window size, moving step size, and polynomial order were set, and the window was slid across the shock detonation experimental dataset using the moving step size. A polynomial fit was performed on the data within the window to obtain the fitted data, thereby completing the filtering. The filtered shock detonation experimental dataset was cleaned by removing abnormal data.

[0070] Specifically, based on the cleaned shock initiation experiment data set, characteristic reference points at different Lagrangian positions in the shock initiation experiment data set corresponding to each shock initiation experiment are extracted. The characteristic reference points represent the data points in the shock initiation experiment data set corresponding to each shock initiation experiment that need to be compared with the numerical simulation results of the calculation model; the characteristic reference points include physical flow field information at different Lagrangian positions. For each Lagrangian position, the first characteristic reference point is the step value when the shock wave / detonation wave reaches the Lagrangian position, and subsequent characteristic reference points are selected according to the pressure-time curve or velocity-time curve corresponding to the physical flow field information; the total number of characteristic reference points at each Lagrangian position should be no less than 20.

[0071] It should be noted that the effective characteristic reference points are selected according to the actual experimental conditions. The selected characteristic reference points should conform to the physical process of shock initiation, and the abnormal points on the corresponding pressure-time curve or velocity-time curve should not be selected. For the experimental results at a Lagrangian position of a shock initiation experiment, generally at the wave arrival time and In the time period between, feature reference points are selected at equal time intervals.

[0072] For example, the shock initiation experiment data set corresponding to a shock initiation experiment includes the physical flow field information at five Lagrangian positions: 0 mm, 3 mm, 6 mm, 9 mm, and 14 mm, such as Figure 3 As shown in Figure 2, the solid line represents the experimental data at each Lagrangian position after cleaning, and the points on the solid line represent the characteristic reference points at each Lagrangian position.

[0073] It should be noted that for different shock initiation experimental methods, the above preprocessing steps can be added or deleted according to the characteristics of the experimental data set and the actual situation to complete the preprocessing. For example, since the above example uses a manganese copper piezoresistive sensor as the pressure sensor, Figure 3 It can be seen that at the three Lagrangian positions of 6mm, 9mm, and 14mm, a V-shaped pressure change occurs after the peak pressure. This is caused by the sensor packaging, so the data corresponding to the V-shaped pressure change should be eliminated.

[0074] As you can understand, the Lagrangian position refers to tracking the position of a specific fluid particle in a flow field. By marking and tracking specific particles in the flow field, the position of each particle can be expressed as a function of time. Physical quantities in the flow field, such as pressure, are functions of particle position and time. By placing pressure sensors or velocity sensors at different locations in each shock detonation experiment, the pressure or velocity of the shock / detonation wave at different Lagrangian positions can be measured, thereby studying the physical flow field characteristics at different Lagrangian positions.

[0075] Specifically, in step S2, the condensed explosive reaction rate model is divided into pressure-based, temperature-based, and entropy-based models, which are used to characterize the rate at which unreacted condensed explosive is converted into detonation products. The nonlinear dynamic analysis software, which uses the finite element method to simulate and analyze complex dynamic problems such as detonation, embeds the reaction rate model into the nonlinear dynamic analysis software so that the software can use the reaction rate model for simulation.

[0076] It should be noted that, during implementation, a condensed matter explosive reaction flow model is usually established to describe the entire process of shock ignition and detonation growth of the condensed matter explosive, including an equation of state, a mixing law, and a reaction rate model. The equation of state includes an equation of state for the unreacted condensed matter explosive and an equation of state for the detonation products, which are used to describe the physical state changes of the unreacted condensed matter explosive and the detonation products, such as the relationship between pressure, temperature, and density. The mixing law is the equilibrium condition of the equation of state and is used to solve the state of the mixture of the unreacted condensed matter explosive and the detonation products. The reaction flow model is embedded in nonlinear dynamic analysis software, enabling the software to more accurately simulate the shock ignition and detonation growth process. Since the parameters of the equation of state can be obtained from other relevant experiments, and the mixing law does not have parameters, this application only considers the optimization of the reaction rate model parameters and does not elaborate on the equation of state and the mixing law.

[0077] For example, the Ignition and Growth model is a pressure-dependent reaction rate model included with the nonlinear dynamic analysis software LS-DYNA. Its equations of state for both the unreacted explosive and the detonation product use the JWL equation of state, and the mixing rule uses the pressure and temperature equilibrium conditions. Its reaction rate model can be expressed as:

[0078] , (1)

[0079] in, Indicates the reaction progress of condensed explosives, Indicates pressure, represents the reaction rate coefficient in the ignition stage, represents the reaction rate coefficient in the growth phase, represents the reaction rate coefficient of the completion stage, 、 、 、 、 、 、 、 、 is the parameter, represents the initial density of the condensed explosive, Indicates the current density of the condensed explosive.

[0080] It should be noted that commercial nonlinear dynamic analysis software includes built-in reaction rate models. This application is applicable to both built-in reaction rate models in nonlinear dynamic analysis software and self-developed reaction rate models. When using a self-developed reaction rate model, such as an entropy-based reaction rate model, it needs to be embedded in the nonlinear dynamic analysis software.

[0081] Furthermore, in the nonlinear dynamic analysis software, the same material properties, geometric model, contact conditions, load conditions and boundary conditions as those in the shock initiation experiment are set, thereby establishing a calculation model consistent with the shock initiation experiment working conditions.

[0082] Specifically, material properties include setting parameters related to the simulation model's material (such as density, elastic modulus, and yield strength) to ensure that the material corresponding to the simulation model is consistent with the condensed explosive used in the shock initiation experiment. The geometric model includes setting the shape, size, and layout of the simulation model's components to ensure that the corresponding geometric model is consistent with the shock initiation experiment. Contact conditions include setting the contact type (such as friction contact or frictionless contact) and contact parameters (such as friction coefficient and contact stiffness) between the simulation model's components to ensure that the contact conditions of the simulation model are consistent with the shock initiation experiment. Load conditions include setting the time history and amplitude of the load applied in the simulation model to be the same as the actual load in the shock initiation experiment. Boundary conditions include setting the boundary constraints in the simulation model to be the same as the actual boundary constraints in the shock initiation experiment.

[0083] For example, a nonlinear dynamic analysis software LS-DYNA is used to create a k-file in LS-DYNA that is consistent with the shock detonation test conditions as a calculation model. The following information is defined in the k-file:

[0084] Material properties: Set material-related parameters of the simulation model, such as density, elastic modulus, and yield strength;

[0085] Geometric model: Set up the same geometric model as the shock detonation experiment, including the shape, size and layout of each component;

[0086] Contact conditions: Set the contact type and contact parameters between the components to be the same as those in the impact initiation experiment;

[0087] Load conditions and boundary conditions: Set the same load loading and boundary constraint conditions as those in the impact initiation experiment;

[0088] Output settings: Set the output data type and frequency.

[0089] It should be noted that consistency with the shock initiation test conditions means that when using nonlinear dynamic analysis software for numerical simulation analysis, the conditions and parameters set must be exactly the same as those in the actual shock initiation experiment. This ensures that the numerical simulation results accurately reflect the physical phenomena in the experiment, thereby improving the reliability of the numerical simulation results. Different nonlinear dynamic analysis software requires different calculation models, and the calculation model is established based on the nonlinear dynamic analysis software actually selected.

[0090] Specifically, in step S3, the expression of the time cost objective function is:

[0091] , (2)

[0092] in, Indicates the total number of shock initiation experiments that need to be compared, represents the total number of Lagrangian positions that need to be compared in one experiment, represents the wave arrival time at the jth Lagrangian position of the i-th experiment in the numerical simulation results obtained by finite element calculation, represents the wave arrival time at the jth Lagrangian position of the i-th experiment in the shock initiation experiment dataset.

[0093] It can be understood that for each experiment, the wave arrival time at the first Lagrangian position should be used as the reference time for subsequent Lagrangian positions.

[0094] For example, Figure 3 As shown in the shock initiation experimental data, the pressure suddenly changes when the wave arrives. The pressure-time curve at a certain Lagrangian position is traversed to find the moment when the mutation amplitude is greater than half of the loading pressure. This moment is the wave arrival time at the jth Lagrangian position of the i-th experiment in the numerical simulation results. .

[0095] It should be noted that for the wave arrival time at the first Lagrangian position of the i-th experiment, the numerical simulation results are and shock initiation experimental data are set to the same moment, so the value of its time cost objective function is determined by the numerical simulation results at the subsequent Lagrangian position and the arrival time of the shock initiation experimental data wave.

[0096] Furthermore, the expression of the physical contour cost objective function is:

[0097] , (3)

[0098] in, Indicates the total number of shock initiation experiments that need to be compared, represents the total number of Lagrangian positions that need to be compared in one experiment, represents the total number of characteristic reference points at each Lagrangian position, It represents the physical flow field information of the kth characteristic reference point at the jth Lagrangian position of the i-th experiment in the numerical simulation results obtained by finite element calculation. Represents the physical flow field information of the kth characteristic reference point at the jth Lagrangian position of the i-th experiment in the shock initiation experiment dataset.

[0099] For example, it is assumed that the physical flow field information is the data of the impact detonation pressure changing with time, which is expressed as a pressure-time curve, and L characteristic reference points are set at the jth Lagrangian position of the i-th experiment. For the jth Lagrangian position of the i-th experiment: the wave arrival time at this Lagrangian position has been determined when calculating the time cost objective function. The curve corresponding to the numerical simulation result is translated to make its wave arrival time consistent with the wave arrival time of the impact detonation experiment dataset. According to the jth Lagrangian position in the impact detonation experiment dataset, the wave arrival time at the jth Lagrangian position has been determined when calculating the time cost objective function. The time of the characteristic reference point is selected from the numerical simulation results after translation, and the pressure corresponding to the data point is .

[0100] Furthermore, using the nonlinear dynamic analysis software used in step S2 and the established calculation model consistent with the shock detonation test conditions, finite element calculations are performed on the condensed explosive reaction rate model embedded in the nonlinear dynamic analysis software to obtain numerical simulation results. By comparing the numerical simulation results obtained by the finite element calculations with the corresponding data in the preprocessed shock detonation test data set using formulas (2) and (3), the values ​​of the time cost objective function and the physical profile cost objective function are obtained.

[0101] Furthermore, the method of finding the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold value as the reaction rate model parameter based on the genetic algorithm and the NASA II algorithm includes:

[0102] Step S31, initializing configuration parameters, generating genetic individuals based on the configuration parameters to construct an initial population; the genetic individuals contain the reaction rate model parameter information;

[0103] Step S32: Calculate the values ​​of the time cost objective function and the physical outline cost objective function corresponding to each individual in the population; select excellent individuals to form the parent population through fast non-dominated sorting and crowding calculation based on the values ​​of the time cost objective function and the physical outline cost objective function corresponding to each individual;

[0104] Step S33, determine whether the termination condition is met. If so, output the optimal solution when the values ​​of the time cost objective function and the physical contour cost objective function corresponding to the individuals in the current population are both less than the threshold value as the reaction rate model parameter; otherwise, adaptively expand the value range of the parameters contained in the individuals in the parent population as needed, and cross and mutate each individual in the parent population according to the crossover probability and mutation probability to obtain a new genetic individual, obtain the offspring population based on the new genetic individual, merge the offspring population and the parent population as the population of the next iteration, and repeat steps S32 to S33.

[0105] Specifically, in step S31, the configuration parameters include the number of individuals in the population, the parameters contained in each individual, the initial boundary of each parameter, the boundary scaling factor of each parameter, the crossover probability, the mutation probability, the mutation intensity, and the maximum number of iterations. Based on the parameters contained in each individual and the initial boundary of each parameter, the genetic individual is generated by encoding, which is expressed as ,in, is the number of parameters included for each individual, is the value of the first parameter of the i-th individual, is the value of the nth parameter of the ith individual. The initial population is constructed by selecting values ​​within the initial bounds of each parameter to generate genetic individuals.

[0106] Furthermore, in order to improve the stability of the calculation results, the following three modes are included when constructing the initial population:

[0107] In the normal mode, a random value is selected within the initial bounds of each parameter to generate each individual in the initial population;

[0108] Selection mode, randomly selects values ​​within the initial boundary of each parameter to generate some individuals in the initial population; the remaining individuals are generated by giving the value of each parameter;

[0109] Restart mode, taking an existing population that meets the initial bounds of each parameter as the initial population.

[0110] It should be noted that the number of parameters for each individual represents the number of genes for each individual, and the reaction rate model parameters that need to be optimized are taken as the parameters contained in each individual; the boundaries of each parameter include an upper boundary and a lower boundary, the upper boundary represents the maximum value of the value range of each parameter, and the lower boundary represents the minimum value of the value range of each parameter; accordingly, the boundary scaling factor of each parameter includes an upper boundary scaling factor and a lower boundary scaling factor; the initial boundary of each parameter is set based on experience.

[0111] Specifically, in step S32, when calculating the values ​​of the time cost objective function and the physical profile cost objective function corresponding to an individual in the population, the parameters corresponding to the individual are substituted into the reaction rate model. Finite element calculations are performed on the reaction rate model using the nonlinear dynamic analysis software and the calculation model to obtain numerical simulation results. The numerical simulation results obtained by the finite element calculations are compared with the corresponding data in the preprocessed impact initiation experimental data set using formulas (2) and (3) to obtain the values ​​of the time cost objective function and the physical profile cost objective function corresponding to the individual.

[0112] Furthermore, when performing finite element calculations on each individual, parallel computing is used to improve computational efficiency.

[0113] For example, suppose the number of individuals in the population is , the maximum number of individuals that the computing resources support for simultaneous calculation is , then there exists a minimum Make , so at most The finite element calculation of all individuals in the population can be completed with only one finite element calculation, thereby greatly improving the calculation efficiency.

[0114] It should be noted that since the time cost objective function and the physical contour cost objective function require finite element calculations to obtain their results, some unreasonable reaction rate model parameters can cause the finite element calculation to stagnate. To avoid wasting computing resources, they must be terminated promptly. Therefore, when the finite element calculation duration of an individual exceeds a time threshold, the calculation of that individual is terminated, marked, and excluded. For example, the individual's corresponding time cost objective function and physical contour cost objective function are both assigned a value of 1000 to exclude it.

[0115] Furthermore, the selecting excellent individuals to form the parent population by fast non-dominated sorting and crowding calculation includes:

[0116] Perform a fast non-dominated sort on all individuals in the population and stratify all individuals according to the Pareto frontier;

[0117] For the initial population, individuals are selected from the stratified individuals through a tournament selection algorithm to form a parent population;

[0118] For non-initial populations, individuals corresponding to the Pareto frontiers of the first n layers are selected from the stratified individuals so that their sum is greater than or equal to the number of individuals in the initial population; if the sum of the individuals corresponding to the Pareto frontiers of the first n layers is equal to the number of individuals in the initial population, the individuals corresponding to the Pareto frontiers of the first n layers are used as the parent population; otherwise, the crowding degree of all individuals corresponding to the Pareto frontier of the nth layer is calculated and sorted from large to small, and the first m individuals are selected so that the sum of the individuals corresponding to the Pareto frontier of the first n-1 layers plus m is equal to the number of individuals in the initial population, and the individuals corresponding to the Pareto frontier of the first n-1 layers and the m individuals selected from the Pareto frontier of the nth layer are used as the parent population.

[0119] It should be noted that it is necessary to ensure that the number of individuals in the parent population is the same as the number of individuals in the initial population.

[0120] Specifically, the Pareto front is divided into several layers, each of which is a set of individuals in a non-dominated relationship. The first-layer Pareto front is the set of individuals that are not dominated by any other individual. When i is greater than 1, the i-th layer Pareto front is the set of individuals that are not dominated by any individual other than in layers 1 to (i-1). If individual a is not inferior to individual b on all objective functions and is superior to individual b on at least one objective function, then individual a dominates individual b; otherwise, individual a does not dominate individual b. If individual a does not dominate individual b, and individual b does not dominate individual a, then individuals a and b are in a non-dominated relationship.

[0121] For example, assuming that the values ​​of the time cost objective function and the physical profile cost objective function of individual a and individual b are (Ta, Pa) and (Tb, Pb) respectively, individual a dominates individual b only in the following two cases: Case 1, Ta ≥ Tb and Pa > Pb; Case 2, Ta > Tb and Pa ≥ Pb.

[0122] Specifically, the fast non-dominated sorting refers to stratifying all individuals in the population by non-dominated relationships, including the following steps:

[0123] Step S301: Pair all individuals in the population and calculate the non-dominated relationships between them. Obtain the set of individuals dominated by each individual and the set of individuals dominating the individual. Select individuals whose dominating set is empty to form the first-level Pareto frontier, which is the best solution set in the population.

[0124] Step S302: Determine whether the number of individuals remaining in the population is 0. If so, the stratification is completed. Otherwise, all individuals remaining in the population are paired and the non-dominated relationship between them is calculated to obtain the set of individuals dominated by each individual and the set of individuals dominating the individual. The individuals whose set of individuals dominating the individual is empty are selected to form the next layer of Pareto frontier, and step S302 is repeated.

[0125] Exemplarily, assuming that the number of individuals in the population is N, all individuals in the population are paired up to calculate the non-dominated relationship between them. For individual k, the set of individuals sk dominated by k and the set of individuals nk dominating k are obtained, and all individuals with ni=empty set (i=1,...,N) constitute the first-level Pareto frontier F1; all remaining individuals in the population are paired up to calculate the non-dominated relationship between them, and the sk and nk corresponding to individual k are updated. All individuals corresponding to ni=0 (i=1,...,N-F1) constitute the second-level Pareto frontier F2; this cycle is repeated until the number of remaining individuals in the population is 0, at which point all N individuals are stratified into n-level Pareto frontiers F1 to Fn.

[0126] Specifically, the crowding degree is used to measure the distribution density of individuals in the target space. For each individual, the crowding degree is the sum of the distances of its adjacent individuals on each objective function; the greater the crowding degree of the individual, the sparser the solutions around the individual and the better the diversity. The values ​​of the time cost objective function corresponding to the individuals in the set are sorted from small to large. For the individuals corresponding to the maximum and minimum values ​​of the time cost objective function, the crowding degree corresponding to their time cost objective function is set to infinity; for other individuals, the crowding degree corresponding to their time cost objective function is calculated using the following formula:

[0127] , (4)

[0128] in, is the number of individuals in the set, is the maximum value of the time cost objective function corresponding to the individual in the set, is the minimum value of the time cost objective function corresponding to the individual in the set, Indicates the The value of the time cost objective function corresponding to each individual;

[0129] Sort the values ​​of the physical contour cost objective function corresponding to the individuals in the set from small to large. For the individuals corresponding to the maximum and minimum values ​​of the physical contour cost objective function, set the crowding degree corresponding to their physical contour cost objective function to infinity. For other individuals, calculate the crowding degree corresponding to their physical contour cost objective function using the following formula:

[0130] , (5)

[0131] in, is the maximum value of the physical contour cost objective function corresponding to the individual in the set, is the minimum value of the physical contour cost objective function corresponding to the individual in the set, Indicates the The value of the physical contour cost objective function corresponding to each individual;

[0132] For each individual, the congestion degree corresponding to its time cost objective function and physical contour cost objective function is accumulated to obtain the congestion degree of the individual.

[0133] For example, assume that the number of individuals in the initial population is N, and for the non-initial population Q, its number of individuals is 2N. Perform a fast non-dominated sort on the population Q to obtain the p-layer Pareto front F1 to Fp, and select the individuals corresponding to the first n layers of the Pareto front F1 to Fn so that their individual sum is ≥ N; if the sum of the individuals corresponding to F1 to Fn is equal to N, then the individuals corresponding to F1 to Fn are used as the parent population; otherwise, calculate the crowding degree of all individuals corresponding to Fn and sort them from large to small, select the first m individuals from them, so that the sum of the individuals corresponding to F1 to Fn-1 plus m is equal to N, and use the individuals corresponding to F1 to Fn-1 and the m individuals selected from Fn as the parent population.

[0134] Specifically, in step S33, the termination condition includes that there is an individual in the current population whose corresponding time cost objective function and physical contour cost objective function have values ​​less than a threshold or the maximum number of iterations is reached.

[0135] Furthermore, the adaptive expansion of the value range of the parameter included in the individuals in the parent population as needed includes:

[0136] Determine whether the value range of each parameter contained in the individuals in the parent population needs to be adaptively expanded; the adaptive expansion includes upper boundary adaptive expansion and lower boundary adaptive expansion;

[0137] If the parameter needs to be adaptively expanded on the upper boundary, the upper boundary of the parameter is multiplied by its upper boundary scaling factor as its upper boundary;

[0138] If the parameter needs to be adaptively expanded at the lower boundary, the lower boundary of the parameter is multiplied by its lower boundary scaling factor as its lower boundary.

[0139] Specifically, for each parameter contained in the individuals in the parent population: check whether it falls on the upper boundary or lower boundary of the parameter. If it falls on the upper boundary of the parameter, it is determined that the parameter needs to be adaptively expanded to the upper boundary; if it falls on the lower boundary of the parameter, it is determined that the parameter needs to be adaptively expanded to the lower boundary.

[0140] For example, assume that each individual contains n parameters, [A1, A2, ..., A n ]; the initial upper bound of each parameter is [M1,M2,…,M n ], the initial lower bound of each parameter is [L1,L2,…,L n ]; the upper boundary scaling factor corresponding to each parameter is [J1,J2,…,J n ], the lower boundary scaling factor corresponding to each parameter is [K1, K2,…, K n ]. In the first round of iteration, if the first parameter A1 falls below its upper boundary M1, it is determined that the parameter needs to be adaptively expanded to the upper boundary, and As its upper bound; if the nth parameter A n There exists a point that falls on its lower boundary L n If the parameter needs to be adaptively expanded to the lower boundary, as its lower boundary.

[0141] Furthermore, two individuals in the parent population are selected as parent individuals based on the crossover probability, and the parameter values ​​contained in the parent individuals are updated using a random multi-point crossover method to generate new genetic individuals; the parameter values ​​contained in the individuals in the parent population are randomly changed based on the mutation probability and mutation intensity to mutate and generate new genetic individuals; the new genetic individuals generated by crossover and mutation are combined to obtain the offspring population.

[0142] It should be noted that if a new genetic individual generated by crossover or mutation already exists, it will be excluded. It is necessary to ensure that the number of individuals in the offspring population is the same as the number of individuals in the initial population. When the termination condition is met, all individuals in the first-level Pareto front are output. The parameters contained in these individuals are the optimal parameter set for the reaction rate model, thereby optimizing the reaction rate model parameters.

[0143] Performance evaluation:

[0144] In order to evaluate the performance of the method of the present invention, the ignition and growth model was used as the reaction rate model for LX-07 explosive, and the I, , G1, c, d, y, G2, e, g and z are optimized. The optimization results based on the method of the present invention are as follows Figure 4 As shown, the optimization results based on James R. Gambino 75% are as follows Figure 5 The parameter optimization results of the two are shown in Table 1:

[0145] Table 1 Parameter optimization results of the present invention and James R. Gambino 75%

[0146] ;

[0147] contrast Figure 4 and Figure 5 As can be seen, at a loading pressure of 5 GPa, the optimization results of the present invention's method correspond to smaller values ​​for the time-cost objective function. At a loading pressure of 2.25 GPa, the optimization results of the present invention are more consistent with experimental trends at the 20 mm and 25 mm Lagrangian positions. At a loading pressure of 1.8 GPa, the reaction rates of the optimization results of the present invention's method are slower at the 30 mm and 35 mm Lagrangian positions, and their von Neumann peaks are more consistent with experimental results. Furthermore, as shown in Table 1, the I value of the optimization results of the present invention's method is much larger than that of James R. Gambino's 75% optimization result. For HMX-based LX-07 explosives, a larger I value is more consistent with common sense, and the G2 value is generally larger than the G1 value. Therefore, the optimization results of the present invention are superior to James R. Gambino's 75% optimization result.

[0148] Compared with the prior art, the optimization method for determining the reaction rate model parameters of condensed explosives provided by the present invention has the following beneficial effects:

[0149] 1. The present invention establishes a time cost objective function and a physical profile cost objective function, and searches for the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold value based on a genetic algorithm and a NASA II algorithm as the reaction rate model parameters, thereby realizing automatic optimization of the reaction rate model parameters. This solves the problem that the current determination of the reaction rate model parameters of condensed explosives relies on manual labor, which is not only time-consuming and labor-intensive but also highly subjective.

[0150] 2. The present invention adaptively expands the value range of the parameters contained in the individuals in the parent population as needed, thereby improving the calculation efficiency while obtaining a reasonable parameter value range.

[0151] 3. The present invention is applicable to the optimization of various reaction rate model parameters and has particular versatility. It can also compare the adaptability of different reaction rate models on the same shock initiation experimental data set, thereby evaluating the reaction rate models.

[0152] Those skilled in the art will appreciate that all or part of the process steps of the above-described embodiments can be implemented by instructing related hardware through a computer program, and the program can be stored in a computer-readable storage medium, such as a magnetic disk, an optical disk, a read-only memory, or a random access memory.

[0153] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by any technician familiar with this technical field within the technical scope disclosed by the present invention should be covered by the scope of protection of the present invention.

Claims

1. An optimization method for determining parameters of a condensed explosive reaction rate model, characterized in that: The method comprises the following steps: Acquiring a condensed explosive shock initiation experimental data set, and preprocessing the shock initiation experimental data set; Establishing a condensed explosive reaction rate model and embedding it into nonlinear dynamic analysis software, and establishing a calculation model consistent with the impact detonation experimental conditions in the nonlinear dynamic analysis software; Establish the time cost objective function and the physical contour cost objective function, where: The expression of the time cost objective function is: ， （2） Where N represents the total number of shock initiation experiments that need to be compared, and M represents the total number of Lagrangian positions that need to be compared in one experiment. represents the wave arrival time at the jth Lagrangian position of the i-th experiment in the numerical simulation results obtained by finite element calculation, represents the wave arrival time at the jth Lagrangian position of the i-th experiment in the shock initiation experiment data set, represents the numerical simulation result of the wave arrival time at the first Lagrangian position of the i-th experiment, represents the shock initiation experimental data of the wave arrival time at the first Lagrangian position of the i-th experiment; The expression of the physical contour cost objective function is: ， （3） Where N represents the total number of shock initiation experiments that need to be compared, M represents the total number of Lagrangian positions that need to be compared in one experiment, and L represents the total number of characteristic reference points at each Lagrangian position. It represents the physical flow field information of the kth characteristic reference point at the jth Lagrangian position of the i-th experiment in the numerical simulation results obtained by finite element calculation. Represents the physical flow field information of the kth characteristic reference point at the jth Lagrangian position of the i-th experiment in the shock initiation experiment data set; Finite element calculation is performed on the reaction rate model parameters using the nonlinear dynamic analysis software and calculation model to obtain numerical simulation results, and the values ​​of the time cost objective function and the physical profile cost objective function are obtained by comparing the numerical simulation results with the corresponding data in the impact initiation experimental data set; based on the genetic algorithm and the NASAⅡ algorithm, the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold is found as the reaction rate model parameter, thereby optimizing the reaction rate model parameters.

2. The optimization method for determining the reaction rate model parameters of condensed explosive according to claim 1, wherein: The preprocessing of the shock initiation experiment data set includes: Filtering and cleaning the shock initiation experiment data set; Based on the cleaned shock initiation experimental data set, characteristic reference points at different Lagrangian positions are extracted; The characteristic reference points represent data points in the shock initiation experiment data set corresponding to each shock initiation experiment that need to be compared with the numerical simulation results of the calculation model; the characteristic reference points include physical flow field information at different Lagrangian positions; The physical flow field information includes data on the impact detonation pressure changing with time and data on the velocity changing with time.

3. The optimization method for determining the reaction rate model parameters of condensed explosive according to claim 1, characterized in that: The method of finding the optimal solution when the values ​​of the time cost objective function and the physical profile cost objective function are both less than a threshold value as the reaction rate model parameter based on the genetic algorithm and the NASA II algorithm includes: Step S31, initializing configuration parameters, generating genetic individuals based on the configuration parameters to construct an initial population; the genetic individuals contain the reaction rate model parameter information; Step S32: Calculate the values ​​of the time cost objective function and the physical outline cost objective function corresponding to each individual in the population; select excellent individuals to form the parent population through fast non-dominated sorting and crowding calculation based on the values ​​of the time cost objective function and the physical outline cost objective function corresponding to each individual; Step S33, determine whether the termination condition is met. If so, output the optimal solution when the values ​​of the time cost objective function and the physical contour cost objective function corresponding to the individuals in the current population are both less than the threshold value as the reaction rate model parameter; otherwise, adaptively expand the value range of the parameters contained in the individuals in the parent population as needed, and cross and mutate each individual in the parent population according to the crossover probability and mutation probability to obtain a new genetic individual, obtain the offspring population based on the new genetic individual, merge the offspring population and the parent population as the population of the next iteration, and repeat steps S32 to S33.

4. The optimization method for determining the reaction rate model parameters of condensed explosive according to claim 3, characterized in that: The method of selecting excellent individuals to form the parent population by fast non-dominated sorting and crowding calculation includes: Perform a fast non-dominated sort on all individuals in the population and stratify all individuals according to the Pareto frontier; For the initial population, individuals are selected from the stratified individuals through a tournament selection algorithm to form a parent population; For non-initial populations, individuals corresponding to the Pareto frontiers of the first n layers are selected from the stratified individuals so that their sum is greater than or equal to the number of individuals in the initial population; if the sum of the individuals corresponding to the Pareto frontiers of the first n layers is equal to the number of individuals in the initial population, the individuals corresponding to the Pareto frontiers of the first n layers are used as the parent population; otherwise, the crowding degree of all individuals corresponding to the Pareto frontier of the nth layer is calculated and sorted from large to small, and the first m individuals are selected so that the sum of the individuals corresponding to the Pareto frontier of the first n-1 layers plus m is equal to the number of individuals in the initial population, and the individuals corresponding to the Pareto frontier of the first n-1 layers and the m individuals selected from the Pareto frontier of the nth layer are used as the parent population.

5. The optimization method for determining the reaction rate model parameters of condensed explosive according to claim 4, characterized in that: The fast non-dominated sorting method is to stratify all individuals in the population by non-dominated relationships, including the following steps: Step S301: Pair all individuals in the population and calculate the non-dominated relationships between them. Obtain the set of individuals dominated by each individual and the set of individuals dominating the individual. Select individuals whose set of individuals dominating the individual is 0 to form the first-level Pareto frontier. These individuals are the best solutions in the population. Step S302: Determine whether the number of individuals remaining in the population is 0. If so, the stratification is completed. Otherwise, all individuals remaining in the population are paired and the non-dominated relationship between them is calculated to obtain the set of individuals dominated by each individual and the set of individuals that dominate the individual. The individuals whose set of individuals that dominate the individual is 0 are selected to form the next layer of Pareto frontier, and step S302 is repeated.

6. The optimization method for determining the reaction rate model parameters of condensed explosives according to claim 5, characterized in that: Sort the values ​​of the time cost objective function corresponding to the individuals in the set from small to large. For the individuals corresponding to the maximum and minimum values ​​of the time cost objective function, set the congestion corresponding to their time cost objective function to infinity; for other individuals, calculate the congestion corresponding to their time cost objective function using the following formula: ， （4） Where N is the number of individuals in the set, T max is the maximum value of the time cost objective function corresponding to the individual in the set, T min It is the minimum value of the time cost objective function corresponding to the individuals in the set, and T(k+1) represents the value of the time cost objective function corresponding to the k+1th individual; Sort the values ​​of the physical contour cost objective function corresponding to the individuals in the set from small to large. For the individuals corresponding to the maximum and minimum values ​​of the physical contour cost objective function, set the crowding degree corresponding to their physical contour cost objective function to infinity. For other individuals, calculate the crowding degree corresponding to their physical contour cost objective function using the following formula: ， （5） Among them, P max is the maximum value of the physical profile cost objective function corresponding to the individual in the set, P min It is the minimum value of the physical contour cost objective function corresponding to the individuals in the set, and P(k+1) represents the value of the physical contour cost objective function corresponding to the k+1th individual; For each individual, the congestion degree corresponding to its time cost objective function and physical contour cost objective function is accumulated to obtain the congestion degree of the individual.

7. The optimization method for determining the reaction rate model parameters of condensed explosives according to claim 3, characterized in that: Adaptively expanding the value range of the parameters contained in the individuals in the parent population as needed includes: Determine whether the value range of each parameter contained in the individuals in the parent population needs to be adaptively expanded; the adaptive expansion includes upper boundary adaptive expansion and lower boundary adaptive expansion; If the parameter needs to be adaptively expanded on the upper boundary, the upper boundary of the parameter is multiplied by its upper boundary scaling factor as its upper boundary; If the parameter needs to be adaptively expanded at the lower boundary, the lower boundary of the parameter is multiplied by its lower boundary scaling factor as its lower boundary.

8. The optimization method for determining the reaction rate model parameters of condensed explosives according to claim 1, characterized in that: In the nonlinear dynamic analysis software, the same material properties, geometric model, contact conditions, load conditions and boundary conditions as those in the shock initiation experiment are set, so as to establish a calculation model consistent with the shock initiation experiment working conditions.

Images