A Parameter Optimization Calculation Method for the Trinomial Ignition and Growth Model of Energetic Materials

Abstract

The present invention discloses a method for optimizing the calculation of parameters of a trinomial ignition growth model of energetic materials. By conducting parameter calibration tests on the trinomial ignition growth model of energetic materials, a calculation model for shock initiation of energetic materials is established. The parameters to be optimized, the parameter value ranges, and the fitness function of the trinomial ignition growth model are defined, the parameters of the trinomial ignition growth model are initialized, and the fitness function value is calculated. Based on the hybrid particle swarm algorithm, the initial parameters, and their fitness function values, the parameters of the trinomial ignition growth model are subjected to crossover and mutation changes, and iterative calculations are performed until the calculation results of the parameters of the trinomial ignition growth model coincide with the test results, and the parameters of the trinomial ignition growth model that coincide with the test results are quickly calibrated. This method is automatic, fast, and efficient, and is suitable for the rapid calibration of equation of state parameters in nonlinear finite element calculations.

Description

Technical Field

[0001] The present invention relates to the technical field of energetic materials, and particularly relates to an optimization calculation method for parameters of a three-term ignition growth model of energetic materials. Background Art

[0002] In 1985, Tarver et al. proposed a three-term ignition growth model. The three-term ignition growth model is one of the explosive detonation reaction models, which can well describe the explosive detonation reaction process and has become the main calculation model for studying the shock initiation of various explosives. The finite element calculation results of the explosive shock initiation process described by the three-term ignition growth model should be maximally consistent with the experimental results. There are numerous parameters in the three-term ignition growth model. To obtain model parameters that are consistent with the experimental results, usually, relying on manual parameter adjustment is cumbersome and accidental, consuming time and effort.

[0003] Aiming at the existing problems, there is an urgent need to provide an optimization calculation method for parameters of a three-term ignition growth model of energetic materials. Summary of the Invention

[0004] The present invention provides an optimization calculation method for parameters of a three-term ignition growth model of energetic materials. This calculation method has the advantages of being fast and accurate, can avoid a large amount of manual calculation, improve the efficiency of obtaining the parameters of the three-term ignition growth model, and the finite element calculation results obtained by using this calculation method for the parameters of the three-term ignition growth model are consistent with the experimental results.

[0005] The present invention adopts the following specific technical solutions:

[0006] An optimization calculation method for parameters of a three-term ignition growth model of energetic materials, the calculation method comprising the following steps:

[0007] Step 1, conduct a calibration test on the parameters of the three-term ignition growth model;

[0008] Step 2, establish a calculation model for shock-initiated energetic materials;

[0009] Step 3, in the three-term ignition growth model, for the JWL equation of state and the reaction rate equation of the unreacted energetic materials, select the parameters to be optimized, define the parameter value range, and at the same time determine the fitness function based on the experimental results;

[0010] Step 4, initialize the parameters of the three-term ignition growth model and calculate the fitness function value;

[0011] Step 5, based on the initial parameters and their fitness function values, through a hybrid particle swarm algorithm with crossover and genetic mutation, perform crossover and mutation changes on the parameters of the three-term ignition growth model;

[0012] Step 6: Iterative calculation is performed until the calculation results of the parameters of the trinomial ignition growth model match the test results.

[0013] Furthermore, in Step 1, a test device for calibrating the parameters of the trinomial ignition growth model is used to conduct a test for calibrating the parameters of the trinomial ignition growth model. The test device for calibrating the parameters of the trinomial ignition growth model includes a detonator, a plane wave generator, a loading energetic material, a polytetrafluoroethylene partition, a tested energetic material, a witness plate, a uniform magnetic field, a trigger probe, a pulse forming network, a single U-shaped electromagnetic particle velocity meter, a combined electromagnetic particle velocity meter, and an oscilloscope.

[0014] Among them, the tested energetic material is composed of two parts of explosive columns that are fitted through a wedge-shaped surface. The combined electromagnetic particle velocity meter is installed on the wedge-shaped surface of the explosive column, and the explosive column is placed in a uniform magnetic field.

[0015] During the test, the detonator detonates the plane wave generator to generate a plane wave. After being loaded by the loading energetic material and attenuated by the polytetrafluoroethylene partition, the shock wave detonates the tested energetic material.

[0016] When the plane wave generated by the plane wave generator reaches the trigger probe, the trigger probe is turned on, and the oscilloscope is triggered through the pulse forming network to record the signal. When the detonation wave reaches the positions of the single U-shaped electromagnetic particle velocity meter and the combined electromagnetic particle velocity meter in sequence, the single U-shaped electromagnetic particle velocity meter and the combined electromagnetic particle velocity meter move together with the detonation products behind the wave. The horizontal working section lengths of the single U-shaped electromagnetic particle velocity meter and the combined electromagnetic particle velocity meter cut the magnetic induction lines of the uniform magnetic field to generate an induced electromotive force ε. When the magnetic field strength B and the horizontal working section length L of the combined electromagnetic particle velocity meter remain unchanged, the induced electromotive force ε generated by the movement of the electromagnetic particle velocity meter is measured through the oscilloscope. According to Faraday's law of electromagnetic induction:

[0017] ε = B × L × v;

[0018] The movement speed v of the sensitive unit in the electromagnetic particle velocity meter can be calculated.

[0019] By changing the thickness of the polytetrafluoroethylene partition, the pressure when the shock wave enters the tested energetic material is adjusted, and the changing trend of the particle velocity inside the initiator charge with position under different incident pressure conditions is obtained.

[0020] In the optimization calculation of the parameters of the trinomial ignition growth model, the parameters of the trinomial ignition growth model are adjusted, and the process of the shock wave with different intensities entering the initiator charge is calculated to obtain the changing trend of the internal particle velocity. Until it matches the changing trend of the internal particle velocity in the test, the parameters of the trinomial ignition growth model of the initiator charge are obtained.

[0021] Furthermore, Step 2 specifically includes:

[0022] Using a non - linear finite - element calculation software, to reduce the computational load, based on the symmetry of the test device for calibrating the parameters of the three - term ignition growth model, a two - dimensional axisymmetric calculation model of shock - initiated energetic materials is established. From top to bottom, it is respectively a plane - wave generator, a loaded energetic material, a polytetrafluoroethylene separator, the energetic material to be tested, and a witness plate. In the calculation software LSDYNA, the upper end of the plane - wave generator is set as the initiation line. In the calculation, the initiation line initiates the plane - wave generator, generating a plane wave. After being stably propagated through the loaded energetic material and attenuated by the polytetrafluoroethylene separator, it initiates the energetic material to be tested.

[0023] The three - term ignition growth reaction rate model is adopted to describe the energetic material to be tested. This model consists of the JWL equation of state for unreacted energetic materials, the JWL equation of state for detonation products, and the reaction rate equation.

[0024] The JWL equations of state for unreacted energetic materials and detonation products are respectively:

[0025]

[0026]

[0027] In the above formula, P E is the initial pressure of the energetic material, P p is the product pressure of the energetic material, V E is the initial specific volume of the energetic material, V p is the product specific volume of the energetic material, C v is the heat capacity, T0 is the initial temperature of the energetic material, T p is the product temperature of the energetic material, and A, B, R1, R2, and ω are undetermined parameters.

[0028] Among them, the reaction rate equation is:

[0029]

[0030] In the above formula, λ is the reactivity of the energetic material, t is the time, ρ is the density, ρ0 is the initial density, P is the pressure, and I, G1, G2, a, b, x, c, d, y, e, g, z are constants; among them, a is the critical compressibility, which is used to define the ignition limit. When the compressibility is less than a, the energetic material does not ignite and does not detonate. Or rather, when the shock wave is strong enough to compress the energetic material to a certain degree, ignition can occur; the pressure exponent y of the combustion term is 1, the burnup order b of the ignition burnup is 2 / 3, and the burnup order c of the combustion term is 2 / 3, indicating inward spherical particle combustion; the parameters I and x control the number of ignition hot spots, and the ignition term is a function of the shock wave intensity and pressure duration; G1 and d control the reaction growth of the hot spots in the early stage after ignition, and G2 and z determine the reaction rate under high pressure; the first term in the reaction rate equation is called the ignition term, the second term is the growth term, and the third term is the rapid reaction term; during the calculation process, the maximum and minimum values of the reactivity λ in the reaction rate equation must be set to control the start and shutdown of different reaction terms; when λ > F igmax , the ignition term is taken as zero; when λ > F G1max ,, the combustion term is taken as zero; when λ < F G2min ,, the rapid reaction term is taken as zero.

[0031] Furthermore, step three specifically includes:

[0032] When calibrating the JWL equation of state parameters of the unreacted energetic material, let the parameters F of the ignition growth model igmax = 0, F G1max = 0 and F G2min = 1, turning off all reaction terms of the energetic material and treating it as an inert material. The design variables are defined as R1, R2, R3, R5, R6 in the JWL equation of state of the unreacted energetic material, and the fitness function is defined as the error based on the takeoff speed and takeoff slope at the highest point and the test results; as the parameters to be optimized, based on the value range of the parameters to be optimized for common energetic materials, the parameter value range is defined;

[0033] When calibrating the reaction rate equation, turn on all reaction terms of the energetic material. The design variables are defined as F MXIG 、F MXGR 、F MNGR 、G1, G2, and the fitness function is defined as the error between the calculated results of the takeoff speed and takeoff slope at the highest point and the test results, serving as the fitness function; when the value of the fitness function is less than 15%, the subsequent iterative calculation based on the hybrid particle swarm optimization algorithm stops.

[0034] Furthermore, step four specifically includes:

[0035] Within the value range of the parameters to be optimized, customize the division of the parameter value intervals, and use the optimal Latin hypercube sampling method and orthogonal experiment to generate a set of parameter combination schemes as the initial values of the parameters of the trinomial ignition growth model; in the hybrid particle swarm optimization algorithm, each particle carries a set of parameter combination schemes, and define the number of particle swarms of the hybrid particle swarm optimization algorithm according to the number of the initial parameter combination schemes.

[0036] Run and compile the Python script, call and start the nonlinear finite element calculation software LSDYNA, use the parametric design language to automatically establish the model, mesh generation, generate calculation files, and solve, and calculate the fitness function values of each particle in the particle swarm.

[0037] Furthermore, step five specifically includes:

[0038] Evaluate the influence of the parameters to be optimized on the fitness function values: for each value interval of the parameters to be optimized, calculate the average value and standard deviation of the corresponding fitness function values, and evaluate whether this value interval of the parameters to be optimized is conducive to the stability and reduction of the fitness function values; for the value interval of the parameters to be optimized with small change range and reduction of the fitness function values, narrow the value interval of this parameter to be optimized; for the value interval of the parameters to be optimized with small change range and increase of the fitness function values, exclude the original value of this value interval of the parameters to be optimized and narrow the value interval; for the value interval of the parameters to be optimized with large change range, do not make adjustments temporarily; feedback the redefined interval to the parameter crossover and mutation.

[0039] Parameter crossover change means that based on the evaluation results of the influence of the parameters to be optimized on the fitness function values, if the value of a certain parameter in the parameter group to be optimized is not conducive to the reduction of the fitness function values, then cross the value of this parameter to be optimized with the value of this parameter in the global optimal combination.

[0040] Parameter mutation change means that based on the evaluation results of the influence of the parameters to be optimized on the fitness function values, if the value of any parameter in the parameter group to be optimized is conducive to the reduction of the fitness function values, then mutate the value of this parameter to be optimized within the current position interval of this parameter to be optimized.

[0041] Compare the fitness function value of each particle in the particle swarm with its own historical best fitness function values, and at the same time compare it with the global best fitness function values of all generations. If the fitness function value of the particle is better than its own historical best fitness function values, then modify the particle's own best fitness function value to the current particle's fitness function value; if the opposite is true, then perform the above crossover and mutation operations.

[0042] Determine whether the current global optimal fitness function value is better than the global optimal fitness function value. If the current global optimal fitness function value is better than the global optimal fitness value, then modify the global optimal fitness value to the current global optimal fitness function value; if otherwise, perform the above-mentioned crossover and mutation operations;

[0043] Determine whether the global optimal fitness function value meets the requirement that the calculation threshold is less than 15%. If it meets, output the optimal parameters of the three-term ignition growth model. If it does not meet, perform the above-mentioned crossover and mutation operations;

[0044] Furthermore, step six specifically includes:

[0045] Update the parameters of the three-term ignition growth model, iteratively calculate the fitness function, compare the fitness function values, modify the global optimal fitness function value and the local optimal fitness function value until the global optimal fitness function value meets the requirement that the calculation threshold is less than 15%, and output the optimal parameters of the three-term ignition growth model.

[0046] Beneficial effects:

[0047] The method for optimizing and calculating the parameters of the three-term ignition growth model of energetic materials according to the present invention conducts parameter calibration experiments on the three-term ignition growth model, establishes a calculation model for shock initiation of energetic materials, defines the parameters to be optimized, the parameter value ranges, and the fitness function values of the three-term ignition growth model; initializes the parameters of the three-term ignition growth model and calculates the fitness function value; based on the initial parameters and their fitness function values, performs crossover and mutation changes on the parameters of the three-term ignition growth model through a hybrid particle swarm algorithm, and iteratively calculates until the calculation results of the parameters of the three-term ignition growth model match the experimental results, quickly calibrating the parameters of the three-term ignition growth model that match the experimental results. This method has the characteristics of being automatic, fast, and efficient, and is applicable to the rapid calibration of equation of state parameters in nonlinear finite element calculations.

[0048] The present invention proposes to apply a hybrid particle swarm algorithm based on genetic variation to the calibration problem of the three-term ignition growth parameters, thereby quickly and accurately obtaining the parameters of the three-term ignition growth model, improving the calculation calibration efficiency in the process of calculating shock initiation of energetic materials, and being applicable to quickly obtaining parameters that accurately describe the shock initiation process of energetic materials. Brief description of the drawings

[0049] Figure 1 It is a flowchart of the method for optimizing and calculating the parameters of the three-term ignition growth model of energetic materials according to the present invention;

[0050] Figure 2 It is a schematic diagram of the principle of the experimental device for calibrating the parameters of the three-term ignition growth model adopted by the present invention;

[0051] Figure 3It is a specific flowchart of the parameter optimization calculation method for the trinomial ignition growth model of energetic materials in the present invention;

[0052] Figure 4 It is a schematic diagram of the calculation model for shock initiation of energetic materials in step two of the present invention;

[0053] Figure 5 It is a comparison diagram of the calculated value and experimental value of the particle velocity of unreacted energetic materials in the present invention;

[0054] Figure 6 It is a comparison diagram of the calculated value and experimental value of the particle velocity of reacted energetic materials in the present invention.

[0055] Among them, 1 - detonator, 2 - plane wave generator, 3 - loaded energetic material, 4 - polytetrafluoroethylene separator, 5 - energetic material to be measured, 6 - witness plate, 7 - uniform magnetic field, 8 - trigger probe, 9 - pulse forming network, 10 - single U - type electromagnetic particle velocity meter, 11 - combined electromagnetic particle velocity meter, 12 - oscilloscope. Specific embodiments

[0056] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0057] The present invention provides a parameter optimization calculation method for the trinomial ignition growth model of energetic materials. The trinomial ignition growth model includes the JWL equation of state of unreacted energetic materials and the reaction rate equation. The JWL equation of state of unreacted energetic materials in the trinomial ignition growth model is described in detail in Embodiment 1, and the reaction rate equation in the trinomial ignition growth model is described in detail in Embodiment 2.

[0058] Embodiment 1

[0059] As Figure 1 and Figure 3 shown, this embodiment provides a parameter optimization calculation method for the trinomial ignition growth model of energetic materials. This calculation method includes the following steps:

[0060] Step 1 S1, conduct parameter calibration experiments on the JWL equation of state of unreacted energetic materials in the trinomial ignition growth model, specifically including:

[0061] As Figure 2As shown in the structure, the experimental device for calibrating the parameters of the trinomial ignition growth model of energetic materials includes a detonator 1, a planar wave generator 2, a loaded energetic material 3, a polytetrafluoroethylene separator 4, a tested energetic material 5, a witness plate 6, a uniform magnetic field 7, a trigger probe 8, a pulse forming network 9, a single U-shaped electromagnetic particle velocity meter 10, a combined electromagnetic particle velocity meter 11, and an oscilloscope 12.

[0062] Among them, the tested energetic material 5 is composed of two parts of explosive columns that are fitted through wedge-shaped surfaces. The combined electromagnetic particle velocity meter is installed on the wedge-shaped surface of the explosive column, and the explosive column is placed in a uniform magnetic field. During the experiment, the detonator 1 detonates the planar wave generator 2 to generate a planar wave. After being loaded by the loaded energetic material 3 and attenuated by the polytetrafluoroethylene separator 4, the shock wave detonates the tested energetic material 5.

[0063] When the planar wave generated by the planar wave generator 2 reaches the trigger probe 8, the trigger probe 8 is conducted. Through the pulse forming network 9, the oscilloscope 12 is triggered to record the signal. When the detonation wave successively reaches the positions where the single U-shaped electromagnetic particle velocity meter 10 and the combined electromagnetic particle velocity meter 11 are located, the single U-shaped electromagnetic particle velocity meter 10 and the combined electromagnetic particle velocity meter 11 move together with the post-wave detonation products. The working section lengths of the horizontal parts of the single U-shaped electromagnetic particle velocity meter 10 and the combined electromagnetic particle velocity meter 11 cut the magnetic induction lines of the uniform magnetic field to generate an induced electromotive force ε. Under the condition that the magnetic field strength B and the working section length L of the horizontal part of the combined electromagnetic particle velocity meter remain unchanged, the induced electromotive force ε generated by the movement of the electromagnetic particle velocity meter is measured by the oscilloscope 12. According to Faraday's law of electromagnetic induction:

[0064] ε = B × L × v;

[0065] the moving speed v of the sensitive unit in the electromagnetic particle velocity meter can be calculated;

[0066] By changing the thickness of the polytetrafluoroethylene separator 4, the pressure when the shock wave enters the tested energetic material is adjusted, and the changing trend of the particle velocity inside the initiator charge with position under different incident pressure conditions is obtained;

[0067] In the optimization calculation of the parameters of the trinomial ignition growth model, the parameters of the trinomial ignition growth model are adjusted, and the process of different-intensity shock waves entering the initiator charge is calculated to obtain the changing trend of the internal particle velocity; until it coincides with the changing trend of the internal particle velocity in the experiment, the parameters of the trinomial ignition growth reaction rate equation of the initiator charge are obtained.

[0068] Step S2, establish a calculation model for shock-wave initiation of energetic materials, specifically including:

[0069] Using the non - linear finite - element calculation software LSDYNA, in order to reduce the computational amount, according to the symmetry of the test device for calibrating the parameters of the three - term ignition growth model, a two - dimensional axisymmetric calculation model of shock - initiated energetic materials is established. As Figure 4 shown, from top to bottom are the plane - wave generator 2, the loaded energetic material 3, the polytetrafluoroethylene separator 4, the tested energetic material 5, and the witness plate 6. In the calculation software LSDYNA, the upper end of the plane - wave generator 2 is set as the initiation line. In the calculation, the initiation line initiates the plane - wave generator 2, generating a plane - wave, which propagates stably through the loaded energetic material 3, and after attenuation by the polytetrafluoroethylene separator 4, initiates the tested energetic material 5;

[0070] The three - term ignition growth reaction rate model is adopted to describe the tested energetic material 5. This model consists of the JWL equation of state for unreacted energetic materials, the JWL equation of state for detonation products, and the reaction rate equation;

[0071] The JWL equations of state for unreacted energetic materials and detonation products are respectively:

[0072]

[0073]

[0074] In the above formula, P E is the initial pressure of the energetic material, P p is the product pressure of the energetic material, V E is the initial specific volume of the energetic material, V p is the product specific volume of the energetic material, C v is the heat capacity, T0 is the initial temperature of the energetic material, T p is the product temperature of the energetic material, and A, B, R1, R2, and ω are undetermined parameters;

[0075] Among them, the three - term ignition growth reaction rate equation is:

[0076]

[0077] In the above formula, λ is the reactivity of the energetic material, t is time, ρ is density, ρ0 is the initial density, P is pressure, and I, G1, G2, a, b, x, c, d, y, e, g, z are constants; among them, a is the critical compressibility used to define the ignition limit. When the compressibility is less than a, the energetic material does not ignite and does not detonate. Or rather, when the shock wave is strong enough to compress the energetic material to a certain degree, ignition can occur; the pressure exponent y of the combustion term is 1, the burn order b of the ignition burnup is 2 / 3, and the burn order c of the combustion term is 2 / 3, indicating inward spherical particle combustion; the parameters I and x control the number of ignition hot spots, and the ignition term is a function of the shock wave intensity and pressure duration; G1 and d control the reaction growth of the hot spots in the early stage after ignition, and G2 and z determine the reaction rate under high pressure; the first term in the reaction rate equation is called the ignition term, the second term is the growth term, and the third term is the rapid reaction term; during the calculation process, the maximum and minimum values of the reactivity λ in the reaction rate equation must be set to control the start and shutdown of different reaction terms; when λ > F igmax , the ignition term is taken as zero; when λ > F G1max , the combustion term is taken as zero; when λ < F G2min , the rapid reaction term is taken as zero.

[0078] Step S3: Select the parameters to be optimized in the JWL equation of state for the unreacted energetic material, define the parameter value ranges, and at the same time determine the fitness function based on the test results, specifically including:

[0079] When calibrating the parameters of the JWL equation of state for the unreacted energetic material, let the parameters F igmax = 0, F G1max = 0 and F G2min = 1, turning off all reaction terms of the energetic material and treating it as an inert material; select R1, R2, R3, R5, R6 in the JWL equation of state for the unreacted energetic material as the parameters to be optimized, define the parameter value ranges based on the value ranges of the parameters to be optimized for common energetic materials, respectively within the ranges of 8000 - 10000, -0.01 - -0.06, 0.00001 - 0.00003, 10 - 15, 1.0 - 2.5, and compare the error between the calculation results and the test results based on the highest takeoff speed and takeoff slope as the fitness function;

[0080] Step S4: Initialize the parameters of the JWL equation of state for the unreacted energetic material and calculate the fitness function value, specifically including:

[0081] Within the parameter value range of the parameter to be optimized, customize the division of the parameter value interval, and use the optimal Latin hypercube sampling method and orthogonal experiment to generate a set of parameter combination schemes as the initial values of the JWL equation of state parameters for unreacted energetic materials. In the hybrid particle swarm optimization algorithm, each particle carries a set of parameter combination schemes, and the number of particle swarms of the hybrid particle swarm optimization algorithm is defined according to the number of initial parameter combination schemes.

[0082] Run and compile the Python script, call and start the nonlinear finite element calculation software LSDYNA, use the parametric design language to automatically establish the model, mesh generation, generate calculation files, and solve, and calculate the fitness function values of each particle in the particle swarm.

[0083] Step S5: Based on the initial parameters and their fitness function values, evaluate the influence of the parameters to be optimized on the fitness function values, and use the hybrid particle swarm optimization algorithm with cross genetic mutation to perform parameter crossover and mutation changes on the JWL equation of state parameters of unreacted energetic materials, specifically including:

[0084] Evaluate the influence of the parameters to be optimized on the fitness function values; for each value interval of the parameters to be optimized, calculate the average value and standard deviation of the corresponding fitness function values, and evaluate whether this value interval of the parameter to be optimized is conducive to the stability and reduction of the fitness function values; for the value interval of the parameter to be optimized with a small change amplitude and a decreasing fitness function value, narrow the value interval of this parameter to be optimized; for the value interval of the parameter to be optimized with a small change amplitude and an increasing fitness function value, exclude the original value of this value interval of the parameter to be optimized and narrow the value interval; for the value interval of the parameter to be optimized with a large change amplitude, do not make adjustments temporarily; feedback the redefined interval to the parameter crossover and mutation;

[0085] Parameter crossover change means that based on the evaluation result of the influence of the parameter to be optimized on the fitness function value, if the value of a certain parameter in the parameter group to be optimized is not conducive to the reduction of the fitness function value, then cross the value of this parameter to be optimized with the value of this parameter in the global optimal combination.

[0086] Parameter mutation change means that based on the evaluation result of the influence of the parameter to be optimized on the fitness function value, if the value of any parameter in the parameter group to be optimized is conducive to the reduction of the fitness function value, then mutate the value of this parameter to be optimized within the current position interval of this parameter to be optimized.

[0087] Compare the fitness function value of each particle in the particle swarm with its own historical optimal fitness function value, and at the same time compare it with the global optimal fitness function value of all generations; if the fitness function value of the particle is better than its own historical optimal fitness function value, then modify the particle's own optimal fitness function value to the current particle's fitness function value; otherwise, perform the above crossover and mutation operations.

[0088] Determine whether the current global optimal fitness function value is better than the global optimal fitness function value; if the current global optimal fitness function value is better than the global optimal fitness function value, then modify the global optimal fitness value to the current global optimal fitness function value; if otherwise, perform the above crossover and mutation operations;

[0089] Determine whether the global optimal fitness function value meets the requirement that the calculation threshold is less than 15%; if it meets, output the optimal parameters of the three-term ignition growth model; if it does not meet, perform the above crossover and mutation operations;

[0090] Step S6, iterative calculation until the calculation result of the JWL equation of state parameters of the unreacted energetic material coincides with the experimental result, specifically including:

[0091] Update the JWL equation of state of the unreacted energetic material, iteratively calculate the fitness function value, compare the fitness function values, modify the global optimal fitness function value and the local optimal fitness function value until the global optimal fitness function value meets the requirement that the calculation threshold is less than 15%, and output the optimal parameters of the JWL equation of state of the unreacted energetic material. Figure 5 The calculated value of the particle velocity of the unreacted energetic material is compared with the experimental value.

[0092] The above model parameter optimization calculation method, through carrying out the parameter calibration test of the three-term ignition growth model, establishing the calculation model of the shock initiation of energetic materials, defining the parameters to be optimized, the parameter value range, and the fitness function value of the JWL equation of state of the unreacted energetic material, initializing the parameters of the three-term ignition growth model, calculating the fitness function, based on the hybrid particle swarm algorithm and the initial parameters and their fitness function values, performing the crossover and mutation changes of the parameters of the unreacted JWL equation of state, and iterative calculation until the calculation result of the parameters of the unreacted JWL equation of state coincides with the experimental result, quickly calibrating to obtain the parameters of the unreacted JWL equation of state that coincide with the experimental result. This method has the characteristics of being automatic, fast, and efficient, and is applicable to the rapid calibration of the equation of state parameters in nonlinear finite element calculations.

[0093] Since the hybrid particle swarm algorithm based on genetic variation is applied to the calibration problem of the three-term ignition growth parameters, the parameters of the unreacted JWL equation of state can be obtained quickly and accurately, improving the calculation calibration efficiency of the process of calculating the shock initiation of energetic materials, and is applicable to quickly obtaining the parameters that accurately describe the shock initiation process of energetic materials.

[0094] Example Two

[0095] This example provides a method for optimizing the calculation of the parameters of the three-term ignition growth model of energetic materials, and this calculation method includes the following steps:

[0096] The first step is to conduct a calibration test for the reaction rate equation parameters of the three-term ignition growth model:

[0097] As Figure 2 shown in the structure, the calibration test device for the three-term ignition growth model parameters of energetic materials includes a detonator 1, a plane wave generator 2, a loaded energetic material 3, a polytetrafluoroethylene partition 4, a measured energetic material 5, a witness plate 6, a uniform magnetic field 7, a trigger probe 8, a pulse forming network 9, a single U-shaped electromagnetic particle velocity gauge 10, a combined electromagnetic particle velocity gauge 11, and an oscilloscope 12.

[0098] Among them, the measured energetic material 5 is composed of two parts of explosive columns fitted through a wedge-shaped surface. The combined electromagnetic particle velocity gauge 11 is installed on the wedge-shaped surface of the explosive column, and the explosive column is placed in the uniform magnetic field.

[0099] During the test, the detonator 1 detonates the plane wave generator 2 to generate a plane wave. After being loaded by the loaded energetic material 3 and attenuated by the polytetrafluoroethylene partition 4, the shock wave detonates the measured energetic material 5.

[0100] When the plane wave generated by the plane wave generator 2 reaches the trigger probe 8, the trigger probe is conducted. Through the pulse forming network 9, the oscilloscope 12 is triggered to record the signal. When the detonation wave reaches the positions of the single U-shaped electromagnetic particle velocity gauge 10 and the combined electromagnetic particle velocity gauge 11 in sequence, the single U-shaped electromagnetic particle velocity gauge 10 and the combined electromagnetic particle velocity gauge 11 will move together with the post-wave detonation products. The horizontal working section lengths of the single U-shaped electromagnetic particle velocity gauge 10 and the combined electromagnetic particle velocity gauge 11 cut the magnetic induction lines of the uniform magnetic field to generate an induced electromotive force ε. By changing the thickness of the polytetrafluoroethylene partition, the pressure when the shock wave enters the measured energetic material is adjusted, and the changing trend of the particle velocity inside the initiator charge with position under different incident pressure conditions is obtained.

[0101] The second step is to establish a calculation model for the shock initiation of energetic materials:

[0102] Using nonlinear finite element calculation software, to reduce the calculation amount, according to the symmetry of the calibration test device for the three-term ignition growth model parameters of energetic materials, a two-dimensional axisymmetric calculation model for the shock initiation of energetic materials is established. As Figure 4 shown, from top to bottom are the plane wave generator 2, the loaded energetic material 3, the polytetrafluoroethylene partition 4, the measured energetic material 5, and the witness plate 6; in the calculation software LSDYNA, a detonation line is set at the upper end of the plane wave generator. In the calculation, the detonation line detonates the plane wave generator to generate a plane wave, which propagates stably through the loaded energetic material and detonates the measured energetic material after being attenuated by the polytetrafluoroethylene partition. Monitoring points are set on the axis of the measured energetic material at positions 0 mm, 2 mm, 4 mm, 6 mm, 8 mm, and 10 mm from the top to monitor the particle velocities at these points.

[0103] The trinomial ignition growth reaction rate model is adopted to describe the energetic material to be measured. This model consists of the JWL equation of state for the unreacted energetic material, the JWL equation of state for the detonation products, and the reaction rate equation.

[0104] The JWL equations of state for the unreacted energetic material and the detonation products are respectively:

[0105]

[0106]

[0107] In the above formula, P E is the initial pressure of the energetic material, P p is the product pressure of the energetic material, V E is the initial specific volume of the energetic material, V p is the product specific volume of the energetic material, C v is the heat capacity, T0 is the initial temperature of the energetic material, T p is the product temperature of the energetic material, and A, B, R1, R2, and ω are undetermined parameters;

[0108] Among them, the trinomial ignition growth reaction rate equation is

[0109]

[0110] In the above formula, λ is the reaction degree of the energetic material, t is the time, ρ is the density, ρ0 is the initial density, P is the pressure, and I, G1, G2, a, b, x, c, d, y, e, g, z are constants; among them, a is the critical degree of compression, which is used to limit the ignition boundary. When the degree of compression is less than a, the energetic material does not ignite and does not detonate. Or rather, when the shock wave is strong enough to make the energetic material reach a certain degree of compression, it can ignite, thus stipulating a necessary condition for the initiation of the energetic material; the pressure index y of the combustion term is 1, the ignition burnup order b is 2 / 3, and the burnup order c of the combustion term is 2 / 3, indicating inward spherical particle combustion; the parameters I and x control the number of ignition hot spots, and the ignition term is a function of the shock wave intensity and pressure duration; G1 and d control the reaction growth of the hot spots in the early stage after ignition, and G2 and z determine the reaction rate under high pressure; considering the characteristics of the above equations, the first term in the reaction rate equation is called the ignition term, the second term is the growth term, and the third term is the rapid reaction term. During the calculation process, the maximum and minimum values of the reaction degree λ in the reaction rate equation must be set to control the start and shutdown of different reaction terms; when λ > F igmax , the ignition term is taken as zero; when λ > F G1max , the combustion term is taken as zero; when λ < F G2min , the rapid reaction term is taken as zero.

[0111] In the third step, define the variables to be optimized in the reaction rate equation of energetic materials, and define the value ranges and fitness functions of the variables to be optimized:

[0112] When calibrating the reaction rate equation, turn on all reaction terms of the energetic material and select F MXIG , F MXGR , F MNGR , G1, and G2 as design variables. Based on the takeoff speed and takeoff slope at the highest point, compare the error between the calculated result and the experimental result as the fitness function.

[0113] In the fourth step, initialize the parameters of the reaction rate equation of energetic materials and calculate the fitness function value:

[0114] Within the value range of the parameters to be optimized, customize the parameter value intervals. Use the optimal Latin hypercube sampling method and orthogonal experiment to generate a set of parameter combination schemes as the initial values of the parameters of the reaction rate equation of energetic materials. In the hybrid particle swarm optimization algorithm, each particle carries a set of parameter combination schemes, and define the number of particle swarms in the hybrid particle swarm optimization algorithm according to the number of initial parameter combination schemes.

[0115] Run and compile the Python script, call and start the nonlinear finite element calculation software, use the parametric design language to automatically establish the model, mesh generation, generate calculation files, and solve, and calculate the fitness function of each particle in the particle swarm.

[0116] In the fifth step, based on the initial parameters and their fitness function values, evaluate the influence of the variables to be optimized on the fitness function value, and use the particle swarm optimization algorithm based on crossover and mutation to perform parameter crossover and mutation changes of the reaction rate equation of energetic materials, specifically including:

[0117] Evaluate the influence of the variables to be optimized on the fitness function value; for each value interval of the variables to be optimized, calculate the average value and standard deviation of the corresponding fitness function value, and evaluate whether this value interval of the variable to be optimized is conducive to the stability and reduction of the fitness function value; for the value interval of the variable to be optimized with a small change range and a decreasing fitness function value, narrow the value interval of this variable to be optimized; for the value interval of the variable to be optimized with a small change range and an increasing fitness function value, exclude the original value of this value interval of the variable to be optimized and narrow the value interval; for the value interval of the variable to be optimized with a large change range, do not make adjustments temporarily; feedback the redefined interval to the parameter crossover and mutation;

[0118] Parameter crossover change means that based on the evaluation result of the influence of the variable to be optimized on the fitness function value, if the value of a certain parameter in the variable group to be optimized is not conducive to the reduction of the fitness function value, then cross the value of this variable to be optimized with the value of this variable to be optimized in the global optimal combination;

[0119] Parameter mutation change refers to, based on the evaluation result of the influence of the parameter to be optimized on the fitness function value, if the value of a certain parameter in the parameter group to be optimized is conducive to the reduction of the fitness function value, then mutate the value of the parameter to be optimized within the current position interval of the parameter to be optimized;

[0120] Compare the fitness function value of each particle in the particle swarm with its own historical best fitness function value, and at the same time compare it with the global best fitness function value of all generations. If the fitness function value of each particle in the particle swarm is better than its own historical best fitness function value, then modify the particle's own best fitness function value to the current particle's fitness function value; if the opposite is true, then perform the above crossover mutation operation. Judge whether the current global best fitness function value is better than the global best fitness function value. If the current global best fitness function value is better than the global best fitness value, then modify the global best fitness value to the current global best fitness function value; if the opposite is true, then perform the above crossover mutation operation. Judge whether the global best fitness function value meets the calculation threshold requirement of less than 15%. If it is equal, output the optimal parameters of the reaction rate equation; if it is not equal, then perform the above crossover mutation operation.

[0121] The sixth step is to perform iterative calculations until the calculation result of the parameters of the energetic material reaction rate equation is in the best agreement with the experimental result:

[0122] Update the parameters of the energetic material reaction rate equation, iteratively calculate the fitness function, compare the fitness function value with the particle best fitness function value and the global best fitness function value, and modify the global best fitness function value and the local best fitness function value until the global best fitness function value meets the calculation threshold requirement of less than 15%, and output the parameters of the energetic material reaction rate equation. Figure 6 Compare the calculated value of the particle velocity calibrated for the reaction rate equation with the experimental value.

[0123] Obviously, those skilled in the art can make various changes and modifications to the embodiments of the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalent technologies, then the present invention is also intended to include these changes and modifications.

Claims

1. A parameter optimization calculation method for the trinomial ignition growth model of energetic materials, characterized in that, Including the following steps: Step 1, conduct a calibration test on the parameters of the three-term ignition growth model; Step 2, establish a computational model for shock initiation of energetic materials, specifically including: Using a non-linear finite element computational software, establish a two-dimensional axisymmetric computational model for shock initiation of energetic materials. From top to bottom, it is respectively a plane wave generator, a loaded energetic material, a polytetrafluoroethylene separator, the energetic material to be measured, and a witness plate; in the computational software LSDYNA, set the upper end of the plane wave generator as the initiation line. In the calculation, the initiation line initiates the plane wave generator to generate a plane wave, which propagates stably through the loaded energetic material and attenuates through the polytetrafluoroethylene separator, and then initiates the energetic material to be measured; Use the three-term ignition growth model to describe the energetic material to be measured, which is composed of the JWL equation of state of the unreacted energetic material, the JWL equation of state of the detonation products, and the reaction rate equation; The JWL equations of state of the unreacted energetic material and the detonation products are respectively: ； ； In the above formula, P E is the initial pressure of the energetic material, P p is the product pressure of the energetic material, V E is the initial specific volume of the energetic material, V p is the product specific volume of the energetic material, C v is the heat capacity, T 0 is the initial temperature of the energetic material, T p is the product temperature of the energetic material, A and B and R 1, R 2 and ω are undetermined parameters; Among them, the reaction rate equation is: ； In the above formula, λ is the reactivity of the energetic material, t is time, ρ is density, ρ 0 is the initial density, P is pressure, I , G 1, G 2, a , b , x , c , d , y , e , g , z are constants; among them, a is the critical compressibility, which is used to define the ignition limit. When the compressibility is less than a , the energetic material does not ignite and does not detonate. Or rather, when the shock wave is strong enough to compress the energetic material to a certain degree, ignition can occur; the pressure exponent y of the combustion term is 1, the burn order b of ignition is 2 / 3, and the burn order c of the combustion term is 2 / 3, indicating inward spherical particle combustion; the parameters I and x control the number of ignition hot spots, and the ignition term is a function of the shock wave intensity and pressure duration; G 1 and d control the reaction growth at the early stage of hot spots after ignition, G 2 and z determine the reaction rate under high pressure; the first term in the reaction rate equation is called the ignition term, the second term is the growth term, and the third term is the rapid reaction term; during the calculation, the maximum and minimum values of the reactivity λ in the reaction rate equation must be set to control the start and shutdown of different reaction terms; when λ > F igmax , the ignition term is taken as zero; when λ > F G1max , the combustion term is taken as zero; when λ < F G2min , the rapid reaction term is taken as zero; Step 3, in the three-term ignition growth model, for the JWL equation of state of the unreacted energetic material and the reaction rate equation, select the parameters to be optimized, define the parameter value range, and at the same time determine the fitness function based on the test results; specifically including: When calibrating the JWL equation of state parameters of unreacted energetic materials, the parameters of the ignition growth model are set F igmax = 0, F G1max = 0 and F G2min = 1, turning off all reaction terms of the energetic material and treating it as an inert material. The design variables are defined as the A , B , R 1, R 2, ω in the JWL equation of state of the unreacted energetic material. The fitness function is defined as the error based on the takeoff velocity and takeoff slope at the highest point and the test results; When calibrating the parameters of the reaction rate equation, all reaction terms of the energetic material are turned on, and the design variables are defined as F igmax 、 F G1max 、 F G2min 、 G 1、 G 2. The fitness function is defined as the error between the calculated results of the takeoff speed and takeoff slope at the highest point and the experimental results; Step 4, initialize the parameters of the three-term ignition growth model and calculate the fitness function value; Step 5, based on the initial parameters and their fitness function values, through a hybrid particle swarm algorithm of cross genetic variation, perform parameter crossing and variation of the three-term ignition growth model; Step 6, perform iterative calculations until the calculation results of the parameters of the three-term ignition growth model match the test results.

2. The calculation method according to claim 1, characterized in that In Step 1, use a calibration test device for the parameters of the three-term ignition growth model to conduct a calibration test on the parameters of the three-term ignition growth model. The calibration test device for the parameters of the three-term ignition growth model includes a detonator, a plane wave generator, a loaded energetic material, a polytetrafluoroethylene separator, the energetic material to be measured, a witness plate, a uniform magnetic field, a trigger probe, a pulse forming network, a single U-shaped electromagnetic particle velocity meter, a combined electromagnetic particle velocity meter, and an oscilloscope; Among them, the energetic material to be measured is composed of two parts of explosive columns that are fitted through a wedge-shaped surface. Install the combined electromagnetic particle velocity meter on the wedge-shaped surface of the explosive column and place the explosive column in a uniform magnetic field; During the test, the detonator initiates the plane wave generator to generate a plane wave, which is loaded through the loaded energetic material and attenuated through the polytetrafluoroethylene separator, and then the shock wave initiates the energetic material to be measured; When the plane wave generated by the plane wave generator reaches the trigger probe, the trigger probe is turned on, and the oscilloscope is triggered by the pulse forming network to record the signal. After the detonation wave reaches the positions where the single U-shaped electromagnetic particle velocity gauge and the combined electromagnetic particle velocity gauge are located in sequence, the single U-shaped electromagnetic particle velocity gauge and the combined electromagnetic particle velocity gauge move together with the post-wave detonation products. The horizontal working section lengths of the single U-shaped electromagnetic particle velocity gauge and the combined electromagnetic particle velocity gauge cut the magnetic induction lines of the uniform magnetic field, generating an induced electromotive force ε. When the magnetic field strength B and the horizontal working section length L of the combined electromagnetic particle velocity gauge remain unchanged, the induced electromotive force ε generated by the movement of the electromagnetic particle velocity gauge is measured by the oscilloscope. According to Faraday's law of electromagnetic induction: ε = B × L × v; the moving speed v of the sensitive unit in the electromagnetic particle velocity gauge can be calculated; By changing the thickness of the polytetrafluoroethylene partition, the pressure when the shock wave enters the measured energetic material is adjusted, and the variation trend of the particle velocity inside the initiator charge with position under different incident pressure conditions is obtained; In the optimization calculation of the parameters of the trinomial ignition growth model, the parameters of the trinomial ignition growth model are adjusted, and the process of the shock wave with different intensities entering the initiator charge is calculated to obtain the variation trend of the internal particle velocity; until it coincides with the experimental variation trend of the internal particle velocity, the parameters of the trinomial ignition growth model of the initiator charge are obtained.

3. The calculation method according to claim 1, wherein Step 4 specifically includes: Within the value range of the parameter to be optimized, the value range of the parameter is customarily divided, and a group of parameter combination schemes are generated by using the optimal Latin hypercube sampling method and the orthogonal experiment as the initial values of the parameters of the trinomial ignition growth model. In the hybrid particle swarm optimization algorithm, each particle carries a group of parameter combination schemes, and the number of the particle swarm of the hybrid particle swarm optimization algorithm is defined according to the number of the initial parameter combination schemes; Run the compiled Python script, call and start the nonlinear finite element calculation software LSDYNA, use the parametric design language to automatically establish the model, mesh generation, generate the calculation file, solve, and calculate the fitness function value of each particle in the particle swarm.

4. The calculation method according to claim 3, characterized in that, Step 5 specifically includes: Evaluate the influence of the parameter to be optimized on the fitness function value: For each value range of the parameter to be optimized, calculate the average value and standard deviation of the corresponding fitness function value, and evaluate whether this value range of the parameter to be optimized is conducive to the stability and reduction of the fitness function value; for the value range of the parameter to be optimized with a small change amplitude and a decreasing fitness function value, narrow the value range of this parameter to be optimized; for the value range of the parameter to be optimized with a small change amplitude and an increasing fitness function value, exclude the original value of this value range of the parameter to be optimized and narrow the value range; for the value range of the parameter to be optimized with a large change amplitude, do not make adjustments temporarily; feedback the redefined range to the parameter crossover mutation; Parameter crossover variation refers to, based on the evaluation result of the influence of the parameter to be optimized on the fitness function value, if the value of a certain parameter in the parameter group to be optimized is not conducive to the reduction of the fitness function value, then cross the value of this parameter to be optimized with the value of this parameter of the global optimal combination; Parameter mutation change refers to, based on the evaluation result of the influence of the parameter to be optimized on the fitness function value, if the value of any parameter in the parameter group to be optimized is beneficial to the decrease of the fitness function value, then mutate the value of the parameter to be optimized in the current position interval of the parameter to be optimized; Compare the fitness function value of each particle in the particle swarm with its own historical optimal fitness function value, and at the same time compare it with the global optimal fitness function value of all generations. If the fitness function value of the particle is better than its own historical optimal fitness function value, then modify the particle's own optimal fitness function value to the current particle's fitness function value; if the opposite is true, then perform the above crossover mutation operation; Judge whether the current global optimal fitness function value is better than the global optimal fitness function value. If the current global optimal fitness function value is better than the global optimal fitness value, then modify the global optimal fitness value to the current global optimal fitness function value; if the opposite is true, then perform the above crossover mutation operation; Judge whether the global optimal fitness function value meets the calculation threshold requirement of less than 15%. If it meets, output the optimal parameters of the three-term ignition growth model. If it does not meet, then perform the above crossover mutation operation.

5. The calculation method according to claim 4, characterized in that, Step six specifically includes: Update the parameters of the three-term ignition growth model, iteratively calculate the fitness function, compare the fitness function values, modify the global optimal fitness function value and the local optimal fitness function value until the global optimal fitness function value meets the calculation threshold requirement of less than 15%, and output the optimal parameters of the three-term ignition growth model.

Images