Nonequilibrium kinetics effects in Richtmyer–Meshkov instability and reshock processes

Yiming Shan ,Aiguo Xu ,Lifeng Wang and Yudong Zhang

1 National Key Laboratory of Computational Physics,Institute of Applied Physics and Computational Mathematics,Beijing 100088,China

2 National Key Laboratory of Shock Wave and Detonation Physics,Institute of Fluid Physics,China Academy of Engineering Physics,Mianyang 621900,China

3 Center for Applied Physics and Technology,MOE Key Center for High Energy Density Physics Simulations,College of Engineering,Peking University,Beijing 100871,China

4 State Key Laboratory of Explosion Science and Technology,Beijing Institute of Technology,Beijing 100081,China

5 School of Mechanics and Safety Engineering,Zhengzhou University,Zhengzhou 450001,China

Abstract Kinetic effects in the inertial confinement fusion ignition process are far from clear.In this work,we study the Richtmyer–Meshkov instability and reshock processes by using a two-fluid discrete Boltzmann method.The work begins by interpreting the experiment conducted by Collins and Jacobs (2002,J.Fluid Mech. 464,113–136).It shows that the shock wave causes substances in close proximity to the substance interface to deviate more significantly from their thermodynamic equilibrium state.The thermodynamic non-equilibrium (TNE) quantities exhibit complex but inspiring kinetic effects in the shock process and behind the shock front.The kinetic effects are detected by two sets of TNE quantities.The frist set includesandwhich correspond to the intensities of the non-organized momentum Flux (NOMF),Non-Organized Energy Flux (NOEF),the flux of NOMF and the flux of NOEF.All four TNE measures abruptly increase in the shock process.The second set of TNE quantities includesand,which denote the entropy production rates due to NOMF,NOEF and their summation,respectively.The mixing zone is the primary contributor to ,while the flow field region outside of the mixing zone is the primary contributor to.Additionally,each substance exhibits different behaviors in terms of entropy production rate,and the lighter fluid has a higher entropy production rate than the heavier fluid.

Keywords: Richtmyer–Meshkov instability,discrete Boltzmann method,nonequilibrium kinetics effects

1.Introduction

The Richtmyer–Meshkov instability [1,2] (RMI) is a fundamental hydrodynamic instability and is widely present in fields such as astrophysical applications [3] and inertial confinement fusion [4–7] (ICF).The interpretation and understanding of the physics of RMI are significant and necessary.In particular,RMI (and reshock) phenomena in ICF significantly affect the implosion performance,which has attracted much attention [8].In addition,an increasing number of studies in recent years have shown that in ICF,due to the presence of high energy density particles and extreme and violent changes in physical quantities,there are strong nonequilibrium kinetic effects,which are nonnegligible and have an impact on the implosion process.In other fields,nonequilibrium behavioral characteristics induced by flow behaviors such as shock waves are also receiving increasing attention [9–23].

Rinderknecht et al [24] noted that the systematic abnormalities in the National Ignition Facility (NIF) implosion dataset,including inferred missing energy in the hohlraum,drive asymmetry in near-vacuum hohlraums,low areal density and high burn-averaged ion temperatures compared with mainline models,and low ratios between the DD-neutron and DT-neutron yields and inferred Ti,suggest that kinetic physics may play a role.Shan et al [25] discussed that the kinetic shocks cause an anomalous large energy spread of the DD neutron signal and anomalous scaling of the neutron yield with the thickness of the CD layers,which are not explicable using hydrodynamic mechanics,and these findings were supported for the first time by experiments and simulations.Cai et al [26] presented a new framework that combines kinetics and hydrodynamics to simulate shock waves and hydrodynamic instabilities in plasmas with high density.This hybrid approach preserves the kinetic effects of ions and their impacts on shock wave propagation,plasma interpenetration,and hydrodynamic instability.For two-fluid simulation,Zhang et al [27] completed nonequilibrium molecular dynamics simulations of shock release in polystyrene (CH) under experimental conditions.Due to its lighter mass than carbon,hydrogen can stream out of the bulk of the CH foil,as observed in these simulations.This released hydrogen generates low-density plasmas in front of the inflight shell,which is in agreement with the experimental findings.Zhang et al noted that for ICF simulations,the kinetic effects of species separation are absent in single-fluid radiation-hydrodynamics codes.Although some studies have begun to pay attention to nonequilibrium kinetic effects,they are still very inadequate due to the lack of research methods,and the study of kinetic effects in ICF is still in its infancy.

The RMI has been extensively studied [28–31] in theory[32–36],numerical simulations [37–50],and experiments[51–66] for general conditions.The physical comprehension of RMI is becoming clearer.However,for general fluid cases,there are also significant nonequilibrium kinetic effects,and this part is missing for the majority of existing studies.For the RMI and reshock processes,in addition to the nonequilibrium effect due to material flow and mixing,shock wave effects are also present,which generate extremely strong nonequilibrium kinetic effects,and the two kinds of nonequilibrium effects interact with each other.Unfortunately,such effects cannot be or are not easily observed and described by traditional macroscopic fluid modeling methods,so studies of nonequilibrium effects present in RMI are scarce.Therefore,it is important to develop corresponding kinetic description methods and to detect and describe the kinetic effects therein.

The numerical experimental study of complex flow includes three main steps: (i) physical modeling,(ii) discrete format selection/design,and (iii) numerical experiments and complex field analysis.The discrete Boltzmann method(DBM) [67–71] is a modeling and analysis method based on the discrete Boltzmann equation and developed from the lattice Boltzmann method [67,72].It focuses on steps (i) and (iii),and leaves the best discrete format design in step (ii) to computational mathematicians.As a modeling method,with increasing the degree of non-equilibrium/non-continuity,the DBM uses more kinetic moments of the distribution function f to capture the main feature of the system.As an analysis method,it provides various measures based on the non-conserved kinetic moments of (f -feq) to detect,describe and exhibit complex Thermodynamic Nonequilibrium (TNE)behavior,where feqis the corresponding equilibrium distribution function.The DBM is a further development of the statistical physical phase space description method.With the introduction of phase space based on the non-conserved kinetic moments of (f -feq),we can further introduce the concept of behavior feature vectors by taking a group of behavior feature quantities as elements.For example,since any definition of the nonequilibrium strength depends on the research perspective,the description of the nonequilibrium behavior of the system obtained by using the nonequilibrium strength from a single perspective is incomplete,often losing much information.Thus,it is helpful to introduce the nonequilibrium strength vector,each element of which describes the degree/extent/intensity of system deviation from equilibrium from its own perspective,to describe the nonequilibrium state of the system from multiple perspectives [68].

Since there are infinite nonconservative moments of the distribution functions,each independent component of any nonconservative moment can provide a corresponding perspective of the nonequilibrium strength.In addition,the time or space change rate of any physical quantity,such as the density,temperature,flow rate,and pressure,can provide a corresponding perspective of the nonequilibrium strength.Other parameters,such as the Knudsen numbers corresponding to different interfaces and the relaxation time of any nonconservative moment,can provide the nonequilibrium strength from the corresponding perspective.Therefore,the nonequilibrium strength vector constructed based on the nonequilibrium strengths from several different perspectives is still only one perspective for describing the nonequilibrium behavior of the system.However,compared with the description based on the nonequilibrium strength from a single perspective,the nonequilibrium strength vector provides a much more accurate description.It should also be mentioned that what DBM includes is a series of physical constraints on the discrete format,instead of the specific discrete format itself.The best discretization scheme is an open research subject of computational mathematics.

The physical function of DBM is roughly equivalent to a set of extended hydrodynamic equations (EHE),where‘extended’ means that,besides the evolution equations of conserved kinetic moments of the distribution function f,evolution equations of some most relevant nonconserved kinetic moments are included.The necessity of the extended part of the EHE increases with increasing the Knudsen (Kn)number.Therefore,being able to recover the hydrodynamic equations in the corresponding level is only one part of the physical function of DBM.To which order of TNE the DBM should include depends on the specific situation.The simplest one that meets the needs is naturally the first choice.

DBM has brought about a series of new insights in several fields,such as for multiphase flow [69,73,74],rarefied gas flow [71,75],combustion and detonation [76–78],and hydrodynamic instabilities [79–81].In contrast to conventional CFD methods,the DBM is not restricted by continuous and near-equilibrium assumptions and can describe cases of high Kn numbers with strong TNE and noncontinuity effects.When designing numerical examples in the literature[70],the description method of nonequilibrium strength vectors from different perspectives was adopted.

As mentioned above,kinetic effects and two-fluid numerical simulations are essential and realistic but lacking in fields such as ICF.In the present paper,we report on singlemode RMI and reshock processes numerically simulated with a two-fluid DBM [76].In the two-fluid DBM,each fluid substance has its own independent distribution function and control equation,and the model can provide nonequilibrium kinetics information for each substance.The numerical simulation conditions are the same as in [82],and the DBM results are in good agreement with the experimental results.Then,we investigate the nonequilibrium kinetic effects,molecular mixing fraction and entropy production rate in RMI and reshock and obtain some new insights.Investigating various nonequilibrium behaviors and kinetic effects in single-mode RMI and reshock problems by using the two-fluid DBM is expected to help improve the physical cognition of hydrodynamic instabilities.The remainder of the article is organized as follows.The two-fluid DBM is presented in section 2.Section 3 presents the numerical simulation results obtained by the DBM and the discussion.The conclusions are presented in section 4.

2.The two-fluid discrete Boltzmann method

The two-fluid discrete Boltzmann-BGK equation can be written as [76]:

where σ is the index of the substance,with σ=A or B in the two-fluid case.i is the index of the discrete velocity,and α is the direction of the spatial coordinate system.andare the discrete distribution function and the discrete equilibrium distribution function of substance A or B,respectively.viαis the discrete molecular velocity,rαis the spatial coordinate,and t is time.τσis the relaxation time.The number density,mass density and velocity of substance σ are:

and the fluid exhibits overall number density n,mass density ρ and velocity u values are:

The substance σ temperature Tσand the mixture temperature T are:

Here,we use the one-step (relaxation collision) model in this paper to simulate the flow field evolution and the interactions of different substance particles,and

cα=vα-uαis the fluctuation of the molecular velocity in the α direction,and c2=cαcα.R is the perfect gas constant.In the one-step relaxation collision model of the two-fluid DBM,the interactions between the two substance particles are expressed in the ‘equilibrium’ distribution functions of each substance fσ,eq: the density,velocity and temperature used to calculate the equilibrium distribution function of each substance are related to both substances.In detail,the density used to calculate the equilibrium distribution function of each substance is the substance A or B particle density ρσ,and the velocity and temperature are the mixture average velocity u and mixture average temperature T.In this way,the interactions between two substance particles are naturally contained in the two-fluid DBM.

It should be noted that the principle of DBM velocity space discretization is that the values of moments used to describe the system must remain unchanged after discretization.The detection and description of nonequilibrium behavior are the core features of the DBM.The nonconservative moments ofare used to characterize the degree of nonequilibrium in the flow field.We define:

where the summation is for the grid nodes.

Furthermore,a kind of non-dimensionalized and normalized nonequilibrium intensity quantities can be defined as below:

These new defined non-dimensionalized and normalized nonequilibrium quantities represent their respective proportions in the total nonequilibrium intensity Dsum.This new approach to define nonequilibrium quantities can also be easily generalized to other nonequilibrium quantities.The two quantities,and the normalized nonequilibrium quantitiesinvestigate how the system deviates from equilibrium in different perspectives.

All these perspectives are reasonable and have complementary effects.For convenience of description,we can construct a nonequilibrium degree vector,each of its components represents the nonequilibrium degree from its own perspective.Two examples of nonequilibrium degree vector are as below:

The entropy production rate is contributed by two parts:viscous stressand heat fulx

The total entropy production rateis:

3.Results and discussion

3.1.Nonequilibrium kinetic analysis of RMI and reshock

In this paper,we use the DBM to reproduce the experimental results with the same simulation parameters as those mentioned in the literature [82] and analyze the nonequilibrium kinetic effects for this example.

As a model construction and complex field analysis method,what the DBM presents are the basic constraints on the discrete formats.The DBM itself does not give specific discrete formats.The specific discrete formats should be chosen according to the specific problem to be simulated.The discrete velocity model used here is as follows:

where k=0,1,…,N and i=1,2,…,M,with N=4 and M=12.vkis the value of the discrete velocity.The discrete velocity model diagram is shown in figure 1.

The third-order Runge–Kutta scheme is used for temporal integration.The fifth-order weighted essentially nonoscillatory (WENO) scheme is used for calculating the spatial derivatives.The flow field parameters are set as follows:space steps dx=dy=5 × 10-5,time step dt=3 × 10-6,and grid number Nx × Ny=128 × 4755.The viscosity of the flow field is found to be small by evaluating the experimental and numerical simulation conditions mentioned in the [60,82].In DBM,viscosity coefficient μ=τP,where P is the pressure.In order to match the experimental conditions,τ is chosen as small as possible,provided that numerical stability is satisfied.In this paper,τ=2 × 10-6.Here,we define the substance below the interface as A and the substance above the interface as B,and ρA＜ρB.Atwood number At=(ρB-ρA)/(ρB+ρA)=0.6.The specific heat ratio of substances A and B is γ=1.25.The shock wave Mach number Ma=1.21 and the shock conditions satisfy the Rankine–Hugoniot relations [83].The interface initial perturbation isy=A0cos(kx),where perturbation wavenumber k=2π/λ=981.75,in which λ=0.0064 is the perturbation wavelength and amplitude A0=0.0002.The shock is initiated at y=0.001 05 in the lighter fluid and propagates in the positive direction of the y-axis.The light-heavy fluid initial interface is located at y=0.0032.Periodic boundary conditions are used for the left and right boundaries.The upper boundary has a reflection boundary condition,and the lower boundary has an inflow boundary condition.The parameters are selected considering grid independence,and figure 2 displays the results of the two-dimensional flow field simulation.

Figure 2(a) shows the initial flow field diagram,and the shock wave propagates along the positive y-axis.After passing through the material interface,RMI develops,and the shock front continues to propagate in the y-axis positive direction,which corresponds to figure 2(b).When the shock front reaches the upper boundary,a reflected shock wave is produced and propagates in the negative y-axis direction,as shown in figure 2(c).The reflected shock then reshocks the interface,and since this reflected shock wave is incident from the heavier fluid to the lighter fluid,the RMI interface is reversed (corresponding to figure 2(d)).As a result of the reshock,a transmitted shock propagates in the negative y-axis direction,while a reflected shock starts to propagate upward,and the reversed RMI interface continues to evolve,as shown in figure 2(e).Since the DBM results are all dimensionless,we recover the magnitude of the amplitude A of RMI and time t for comparison with the results in [60,82],and the details on recovering the magnitude in the DBM can be found in [84].Figure 3 shows the comparison between the experimental results mentioned in [60,82] with the DBM simulation results after recovering the magnitude.The DBM results are in good agreement with the experimental results,which indicates that the DBM including only the first order TNE,the time step,space step,discrete velocity model and other parameter settings in the DBM program meet the needs of the current study.

As the shock front continues to propagate along the positive y-axis,it interacts with the upper boundary at some point,generating a reflected wave propagating in the negative y-axis direction,and the reflected wave interacts with the substance interface again,i.e.reshock.For reshock,we define the mixing degree Mix:

where FAand FBare the mass fractions of substances A and B,respectively.Therefore,the range of Mixis 0 ≤Mix≤1.Mix=1 means complete mixing,and Mix=0 means no mixing.The change in the mixing degree for the two substances in the mixing zone before and after the reshock is depicted in figure 7.The mixing degree at the substance interface is 1,which represents the most fully mixed positions,and it decreases from the interface outward.From figures 7 (b)–(f),the shock front starts to interact with the interface,it can be seen that the range of the mixing zone is initially compressed and then expanded,and the RMI interface reverses due to reshock.

In the reshock process,we calculate the evolution of the different orders of the nonequilibrium intensityand(see fgiures 8–11),where (a) to (f) in these figures correspond to (a)–(f) in figure 7,respectively.Forfigure 8(a) displays that before the shock front reaches the mixing zone,there are nonequilibrium effects in the mixing zone because the RMI is developing.When the shock front reaches the mixing zone,we can observe in figure 8(b)that an interrupted interface is caused by the strong nonequilibrium effects of the shock wave,whose values are much greater than the nonequilibrium intensity in the mixing zone.In figure 8(c),the reshock wave reaches the interface,is reflected and transmits within the mixing zone,and the shock front deforms For figures 8 (d)–(f),since the reshock wave has already passed through the mixing zone,what propagates in the mixing zone at this time is the transmitted shock generated by the interaction between the shock front and substance interface.The nonequilibrium situation in the mixing zone is complicated,but the intensity is generally decreasing.

Figures 10 and 11 show the evolution ofandin the mixing zone during the reshock.It has been discovered that the diagrams ofandandare comparable,respectively,but they are not identical.In addition,the distribution of each order of the nonequilibrium quantities in the reshock demonstrates that the nonequilibrium quantity diagrams exhibit a pattern similar to that obtained by the Schlieren method,which is commonly used in shock wave experiments,and display flow field details that cannot be obtained by direct observation of the flow field.This demonstrates that,similar to the Schlieren method,the nonequilibrium quantities help grasp the flow details within the flow fields,which is one of the advantages of the DBM.The above results depict the global nonequilibrium effect of the mixing zone during the reshock process.Due to the implementation of the two-fluid DBM in this study,each substance has its independent distribution function f and Boltzmann evolution equation,allowing information regarding the flow of each substance to be separately provided,which is unavailable using the single-fluid simulation method.To analyze the detailed behavior of each substance during the reshock process,here,we provide the nonequilibrium intensity distributions of substances A and B in the mixing zone.

Figure 12 shows the various orders of the nonequilibrium intensity of each substance in the reshock process.It comprises graphs (A) to (F),which correspond to the six moments in figure 7,similar to the previous figures.In addition,(A)–(F) also contain (a),(b),(c) and (d),which representandrespectively.Overall,we can observe that there are overlaps between the nonequilibrium intensity diagrams of substances A and B due to the diffusion effect between the two substances,as shown in figure 12.For the spike and bubble structures in figure 12 (A),which are already in a nonlinearly stable development stage,the nonequilibrium effects generated by substances A and B at the spike and bubble are not consistent.However,in general,the nonequilibrium intensityis greater within the substance that is ‘forced’ to change its flow state,as shown in thepart of figure 12 (A) (for example,since the spike develops in substance A,A is here called the substance that is ‘forced’ to change its flow state,and the opposite for the bubble).It can be seen from figure 12 (B) that although the existence of temperature differences between two substances leads to the generation of a heat flux at the macroscopic level if the two substances are viewed separately,there is a difference in the heat flux generated by each of the two substances for the one process due to their different characteristics.The macroscopic manifestation of the flow system behavior is determined by the combined response of both substances.In figure 12 (B),the reshock wave has reached the mixing zone and started to interact with the interfacial material,and the nonequilibrium effect due to the presence of the shock wave is captured in the nonequilibrium intensity images of both the A and B substances,whereas the nonequilibrium effect due to the RMI evolution within the mixing interface is not obvious in theandproflies because it is too small in comparison to the nonequilibrium effect produced by the shock front.This conclusion is identical to that for the first shock.In addition to leading to discontinuity of macroscopic physical quantities,the shock wave also leads to discontinuity of nonequilibrium quantities,as well as expansion of the nonequilibrium region and intensification of the nonequilibrium effect,as shown in theandproflies.In fgiure 12 (C),the propagation of the shock front in the mixing zone is captured.Specifically,theandproflies depict the propagation of the shock front in the mixing zone,and in the upper portion of theandproflies of substance B,three regions exhibit a high degree of nonequilibrium intensity.The distributions ofandreveal that the spike is reversing at this point as a result of the shock.The above is also illustrated in theanddistributions.Despite the fact that this conclusion can be reached from both profiles,the phenomenon is not expressed in the same way.In figures 12 (D)–(F),the reshock wave has already passed through the mixing zone.However,due to the interaction between the shock front and substance interface,the transmitted shock waves propagate within the mixing zone.In addition to the nonequilibrium effect generated by the shock wave,the nonequilibrium effect generated by the interface reversal and mixing at the substance interface is also captured,and the two kinds of nonequilibrium effects interact with each other.

Figure 1.Discrete velocity model.

Figure 2.Flow field evolution of single-mode RMI and reshock: (a) t=0,(b) t=0.06,(c) t=0.162,(d) t=0.264,and (e) t=0.3.

Figure 3.Comparison of amplitude A with time between DBM results and experimental data.

Figure 4.Different orders of the nonequilibrium intensity at the substance interface before and after the shock wave passes through the interface for the first time: (a) ∣Δ*2 ∣,(b) (c) and (d)the left and right images of each picture represent before and after the shock wave passes through the interface,respectively.

Figure 5.Nonequilibrium quantity distribution of substances A and B along the y-axis through the bubble before (a) and after (b) the shock.

Figure 6.Nonequilibrium quantity distribution of substances A and B along the y-axis through the spike before (a) and after (b) the shock.

Figure 7.Mixing degree versus time in the mixing zone during reshock: (a) before the reshock wave reaches the mixing zone at t=0.246;(b)to (f): the reshock wave interacts with the substance interface at t=0.252,0.258,0.264,0.270 and 0.276.

Figure 9.versus time in the mixing zone during reshock;(a)–(f) in this fgiure correspond to (a)–(f) in fgiure 7,respectively.

Figure 10.versus time in the mixing zone during reshock;(a)–(f) in this figure correspond to (a)–(f) in figure 7,respectively.

Figure 11.versus time in the mixing zone during reshock;(a)–(f) in this fgiure correspond to (a)–(f) in fgiure 7,respectively.

Figure 12.Different orders of the nonequilibrium intensity of substances A and B in the reshock process.(A)–(F) correspond to the six moments referred to in figure 7.

Figure 13.Variation of different orders of the nonequilibrium intensities in the mixing zone with time t.

Figure 14.Normalized nonequilibrium intensity quantities in the mixing zone with time t.

The definition of the molecular mixing fraction/ratio[85,86] is:

The molecular mixing fractions of the mixing zone and the whole flow field are denoted by Rmand Ra,respectively.Figure 15 shows the molecular mixing fraction Rmand Raprofiles versus time t.Even for the same process,different study scopes result in distinct variations in the molecular mixing fraction,as shown in figure 15.For Rm,RMI begins to occur when the initial shock front impacts the interface,and Rmsignificantly increases.Then,as the shock front passes through the mixing zone and RMI continues to develop,Rmslowly increases as the fluid instability develops.Reshock occurs at t=0.246,and at the stage when the reshock wave interacts with the substance interface,both the decrease in the width of the mixing zone due to the shock front impacting the mixing zone and the promotion of the mixing of the substances caused by the reshock result in a rapid increase in Rm.Then,the mixing zone starts to expand as the shock front passes through the mixing zone,which corresponds to the inclusion of unmixed fluid into the mixing zone,resulting in a decrease in the degree of mixing in the mixing zone.However,as the RMI continues to develop,Rmincreases,and the rate is greater than that in the initial shock,which indicates that reshock can promote the mixing of substances in the flow field.

Figure 15.Molecular mixing fraction profiles versus time t.Rm:mixing zone;Ra: whole flow field.

Figure 16.Entropy production rate of substances A and B in the whole field.(a) of substances A and B;(b) of substances A and B;(c) entropy production rate of substances A and B.

For Ra,initially,Raincreases when RMI begins to develop after the shock front impacts the interface.Subsequently,Ragradually increases as the RMI interface evolves.When t=0.246,reshock occurs,and Radeclines,which is contrary to the variation in Rmat this time.This occurs because although the degree of mixing increases,the area of the mixing zone is reduced,leading to a decrease in Ra.However,when the shock front passes through the mixing zone,Rabegins to increase at a much faster rate than the first shock.The above shows that reshock promotes the flow field substance mixing,but the degree of mixing does not monotonically increase.

3.2.Entropy production rate

4.Conclusions and discussions

The entropy production rates are the second set of TNE measures in the current work.Regarding the entropy production rate of each substance,the following results can be obtained: for the RMI and reshock processes studied in this paper,(i) the mixing zone is the primary contributor towhile the flow field region excluding the mixing zone is the primary contributor to(ii) Since the shock wave is weaker in the reshock process,the increase in the entropy production rate during reshock is less than that observed during the initial shock.(iii) Each substance behaves differently in terms of entropy production rates,with the lighter fluid having a higher entropy production rate than the heavier fluid.

Acknowledgments

The authors thank Feng Chen,Junfeng Wu,Chao Zhang,Zheng Yan,Jiahui Song,Dejia Zhang,Jie Chen,Hanwei Li,and Yingqi Jia for helpful discussions on the DBM.The authors also thank the anonymous reviewers for their valuable comments and suggestions.This work was supported by the National Natural Science Foundation of China (under Grant Nos.12172061,11875001,11575033,and 11975053),the opening project of the State Key Laboratory of Explosion Science and Technology (Beijing Institute of Technology)(under Grant No.KFJJ23-02M),Foundation of National Key Laboratory of Shock Wave and Detonation Physics,and the Foundation of National Key Laboratory of Computational Physics.