• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Molecular Dynamics Simulation of Shock Response of CL-20 Cocrystals Containing Void Defects

    2024-02-29 08:23:12ChnglinLiWeiYngQingGnYjunWngLinLingWenoZhngShungfeiZhuChnggenFeng
    Defence Technology 2024年1期

    Chnglin Li , Wei Yng , Qing Gn ,*, Yjun Wng , Lin Ling , Weno Zhng ,Shungfei Zhu , Chnggen Feng

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

    b School of Materials Science & Engineering, Beijing Institute of Technology, Beijing 100081, China

    Keywords: CL-20 co-crystals Molecular dynamics simulation Reactive force field Impact response Hot spot Void defect

    ABSTRACT To investigate the effect of void defects on the shock response of hexanitrohexaazaisowurtzitane(CL-20)co-crystals, shock responses of CL-20 co-crystals with energetic materials ligands trinitrotoluene (TNT),1,3-dinitrobenzene (DNB), solvents ligands dimethyl carbonate (DMC) and gamma-butyrolactone (GBL)with void were simulated,using molecular dynamics method and reactive force field.It is found that the CL-20 co-crystals with void defects will form hot spots when impacted, significantly affecting the decomposition of molecules around the void.The degree of molecular fragmentation is relatively low under the reflection velocity of 2 km/s, and the main reactions are the formation of dimer and the shedding of nitro groups.The existence of voids reduces the safety of CL-20 co-crystals, which induced the sensitivity of energetic co-crystals CL-20/TNT and CL-20/DNB to increase more significantly.Detonation has occurred under the reflection velocity of 4 km/s,energetic co-crystals are easier to polymerize than solvent co-crystals,and are not obviously affected by voids.The results show that the energy of the wave decreases after sweeping over the void, which reduces the chemical reaction frequency downstream of the void and affects the detonation performance, especially the solvent co-crystals.

    1.Introduction

    CL-20,as one of the most important energetic materials applied,was limited application due to its disadvantages of high sensitivity[1] and low mechanical properties [2].The co-crystals formed by CL-20 and other compounds can greatly reduce the sensitivity [3],improve the safety of materials and effectively expand the scope of application on the basis of retaining high energy,which provides a new idea for the preparation of high-energy and low-sensitivity energetic materials.The design of co-crystals is an effective method to improve its properties from the molecular level [4],which could also balance its high energy and insensitive properties[5],better than the mixed explosive[6].At present,a variety of CL-20 co-crystals with energetic materials(TNT,DNB,benzotrifuroxan(BTF) [7-10], etc.) have been prepared, showing excellent performance in high energy and insensitive.In addition,the co-crystals of CL-20 prepared with solvents(DMC[11],GBL[5],etc.)also benefit from reducing its sensitivity.To evaluate the property of CL-20 cocrystals, its shock response characteristics gained much attention recently,of which molecular dynamics simulation[12]method and reactive force field [13] with low gradients correction (ReaxFF-lg)[14]show an advantage.The initial parameters of ReaxFF lg reactive force field are obtained by training the physical and chemical properties of nitromethane,TATB,RDX,and PETN,and are generally applicable to them.However, as a nitro compound, CL-20 has similar physical and chemical properties to the above compounds,so ReaxFF-lg reactive force field has also been used by many researchers to explore the microscopic properties of CL-20 [15-17],and has drawn a conclusion consistent with the experimental results.

    Recent shock simulations focus on the relationship between shock velocity and the response of CL-20 co-crystals[18].Liu et al.[15] studied the influence of shock velocity on the response of CL-20/TNT co-crystal, CL-20 would be polymerized to form a dimer and no decomposition occurred under impact at low velocity(0.8-1.0 km/s) while under shock velocity of 2-4 km/s, the cocrystals structure decomposition caused by the impact is dominant,and the amount of dimer decreases gradually.Zhang et al.[16]simulated the shock decomposition process of CL-20/BTF, the results showed that the shock wave sensitivity of CL-20 at different directions has slight differences.According to Liu et al.[17],CL-20/HMX co-crystals undergo periods of induction, fast compression,slow compression,and expansion under shock velocity of 4-10 km/s,and different wavelet front velocities were observed in CL-20 and HMX, respectively.However, During the preparation of energetic materials, defects of different kinds, shapes, and sizes often occur.The above studies consider little about the effect of defects on the shock response of CL-20 co-crystal, which would come from the preparation and crystallization processes, such as voids, doping,and dislocations [19-21].Although Hang et al.[22-24] revealed the adverse effects of defects on the sensitivity of CL-20/NQ,CL-20/TNT, and CL-20/DNB co-crystals by static equilibrium model, the shock response of CL-20 co-crystals with defects have not been reported.

    Fig.1.(a) Molecular structures of CL-20/TNT; (b) CL-20/DNB; (c) CL-20/DMC; (d) CL-20/GBL co-crystals molecular structures and perfect/defective co-crystals supercells.

    Table 1Cell parameters of built CL-20 perfect and void defect co-crystals.

    In this paper,in order to study the impact of defects on the shock response of CL-20 co-crystals,and the effect of the ligands,the void collapse,initial chemical reaction,and hot spot formation of CL-20/TNT, CL-20/DNB, CL-20/DMC and CL-20/GBL co-crystals were studied by molecular dynamics simulation and ReaxFF-lg reactive force field,the influence of defects on shock sensitivity of CL-20 cocrystals was discussed from the molecular level, aims to provide a reference for the shock response mechanism of CL-20 co-crystals.

    2.Methods

    2.1.Model establishment and calculation

    The supercells of perfect and void-containing energetic cocrystals CL-20/TNT [7] and CL-20/DNB [9], solvent co-crystals CL-20/DMC [11] and CL-20/GBL [5] are established based on crystal data of Cambridge crystal database, shown in Fig.1.The specific sizes and atomic numbers of each supercell are shown in Table 1.All simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [25].Periodic boundary conditions are applied to these supercells in thea,b, andcdirections.The conjugate gradient algorithm [26] is used to minimize lattice energy, Maxwell Boltzmann random distribution [27]is used to assign the initial velocity of the atom,and the Berendsen thermostat[28]is used to relax the supercell for 15 ps at 298 K.The relaxation structure is obtained by 15 ps relaxation at 298 K and 0 GPa with 0.1 fs timestep using the NPT ensemble and Nos'e-Hoover thermostat and barostat [29,30], and the atomic nearest neighbor table is generated every ten timesteps.The comparison of the size and density of the perfect supercell after relaxation with the experimental value is shown in Table 2.It can be seen that the error between the simulation value of the lattice constant and the experimental value is small, which shows that the ReaxFF-lg force field can accurately describe the physical properties of CL-20 cocrystals.

    Set directioncas a nonperiodic boundary condition and set up a reflective wall [31] to realize the shock induced simulation.Supercells whose void radii are 0 and 20 ? are respectively reflected along directioncat a velocity of 2 km/s and 4 km/s.The total simulation time lasts 5 ps,and 0.1 fs timestep were used during the simulation.The spatial position of the atom, molecular type information and position are output every 50 fs, the spatial position of the atom is output every 0.04 ps,the bond level is calculated every 1 fs, and the average bond level is calculated every 5 fs.Derive bond-level data,and determine the breakage and formation cut-off of interatomic bonds based on the bond-level table between C,H,O,and N[32].Then observe the phenomenon of void collapse and hotspot generation during shock process, count the number of molecules generated and the frequency of chemical reactions, and compare the results of perfect and void-containing co-crystals.

    Table 2Cell parameters and density values (ρ/(g·cm-3)) of CL-20 Perfect co-crystals cells after relaxation.

    2.2.Data processing method

    (1) Wavefront analysis method of shock model

    In order to study the wavefront moving state after the shock,the shock data at different times are analyzed to investigate the propagation process of the shock wave in the model.The wavefront velocity under each impacted state is obtained by linear fitting.The position determined by the average value of the two velocities in front of and behind the wavefront is taken as the wavefront position.The shock wave position-time curve is drawn, and the wave velocity is approved by the shock Hugoniot formula.

    According to the average velocity of each particle along the impact direction after the shock wave sweep, determine the particle velocity behind the wave of the relatively static reference system,and then combine it with the initial velocity of the impact model to obtain the particle velocity behind the wave relative to the impact model.

    μpis the particle velocity behind the wave, μp0is the particle velocity behind the wave of the relatively static reference system,andVs0is the initial shock velocity of the system, in km/s.In this paper,the reflective wall impact method is used to give the model a velocityVs0in the opposite direction to the impact direction,so the wave velocityDof the wavefront relative to the impact model is as follows:

    whereD0is the wavefront velocity of the stationary reference system, in km/s.The above equation and the impact Hugoniot equationPn=ρ0upD were combined as

    wherePnis the wavefront pressure,unit GPa;ρ0is the density of the system before impact, g/cm3; D0is the wavefront velocity of the relative static reference system,and Vs0is the initial shock velocity of the system, in km/s.Since the average velocity of each impact model atom along the impact direction after wavefront scanning is close to 0 [34],μp0is regarded as 0.

    (2) Analysis of two-dimensional temperature distribution

    In order to investigate the morphology of high-temperature hot spots caused by void collapse,a two-dimensional temperature field diagram was calculated.The regional temperature analysis is based on the energy equipartition theorem[33],which connects the total translational kinetic energy and total energy of atoms in the region,and converts them to obtain the regional temperature.The average velocity of the centroid was subtracted from the velocity of the atom to obtain the relative velocity of the atom, the relative translational kinetic energy was calculated, then the energy equalization theorem was applied to calculate the temperature in the region.The average value of the velocity of the atom in the region swept by the shock model wavefront is close to 0 [34], the atomic translational kinetic energy in this region is about the relative translational kinetic energy, and the temperature of the area is not affected by the wavefront is regarded as 298 K.The influence region of the wavefront determines the temperature based on the average translational kinetic energy, which is determined by the energy equalization theorem.

    whereCkeis the translational kinetic energy in kJ/mol;NAis the Avogadro constant with the value of 6.02×1023;kis the Boltzmann constant with the value of 1.38×10-23J/K;Tis the temperature of the region, unit K.According to the formula, the conversion relationship between the translational kinetic energy of atoms and the temperature in the region can be obtained as follows:

    Each model region is divided into 100 × 40, calculate the average temperature of each region, thus drawing the twodimensional temperature field diagram of void collapse.

    (3) Analysis of atomic trajectories and molecular species

    In the process of molecular dynamics simulation, the positions of all atoms in the system follow the equations of motion determined by the forcefield.The physical quantity of each atom is output at a preset timestep of atomic information.The Open Visualization Tool (OVITO) [35] is used to process the trajectory data, and the atomic number, atomic type, charge, atomic coordinates, atomic velocity, kinetic energy, potential energy, and other physical quantities are output every 0.1 ps.In addition, in order to avoid accidental chemical bonds, the cutoff of average bond order [32] is created to determine whether there is a bond between atoms, thus identify the species of molecular fragments under the condition of maintaining the integrity of molecules.

    In order to investigate the chemical reaction rules during each shock, the main reaction products and their distribution were analyzed.The products include undecomposed CL-20 molecules,co-crystal ligand molecules, and chemical species such as N2,NO2,NO3, HCN, H2, CO2,and H2O.In the simulation process, the types and positions of molecular fragments are derived every 0.1 ps,and the relative number and distribution of the main groups are analyzed over time, so as to obtain the curve of the number of molecules over time and the distribution of each molecule in the supercell at a specific time.The relative contentCmof particles is defined as follows:

    wherenis the number of molecules in the system andn0is the number of molecules in the unreacted system.In the initial unreacted state,Cmof CL-20 molecule equals to 1.

    (4) Chemical elementary reaction analysis

    In order to investigate the chemical reaction path in the two cases of shock compression and wavefront reverse stretching,here uses the chemical reaction frequency statistical method to analyze.This method obtains the chemical reaction through the change of the bonding relationship between adjacent timesteps (that is, it is regarded as one step completed within the step interval),and does not discuss the reaction process within the step interval.An example of elementary reaction discrimination is shown in Fig.2.

    Fig.2.Examples of organic substitution reaction.

    As shown in Fig.2, each group is numbered at past and next timesteps, obtain the affiliation of each atomic group to the molecule.Analyze the changes of molecular groups at adjacent timesteps,tracing the new affiliated molecules of atoms after reactions,and traverse the molecular structure containing the newly assigned atoms,finally integrate the results obtained from the traversal into the elementary reaction formula.In addition, this article uses CL-20/DNB as the benchmark to normalize the reaction frequency obtained by this method,which is convenient for the comparison of different crystal orientation supercells.

    3.Results and discussion

    3.1.Analysis of shock wavefront displacement

    The wavefront position of different perfect impact models with time is shown in Fig.3.Determine the wavefront moving speed of each impact speed according to Eq.(2),and calculate the wavefront pressure using the impact Hugoniot formula, which is compared with the simulation results in this paper, as shown in Table 3.The corresponding column analysis diagram is shown in Fig.4.

    It can be seen from Fig.3 that the wavefront velocity of CL-20/DMC and CL-20/GBL solvent co-crystals under 2 km/s impact is significantly lower than that of CL-20/TNT and CL-20/DNB energetic co-crystals,while the difference is small under 4 km/s impact,so it can be inferred that energetic co-crystals is more prone to detonation.

    The error between the simulation result of pressure and the calculation result of impact Hugoniot formula is less than 5%,which proves that ReaxFF-lg reactive force field can accurately express the impact response of CL-20 co-crystals.The wavefront pressure at 2 km/s shock is about 30 GPa, which is less than the CL-20 cocrystals detonation pressure, indicating that these two velocity shocks are not enough to cause CL-20 detonation.However,under the impact of 4 km/s velocity, the pressure of each co-crystal’s wavefront has been greater than the detonation pressure.It is believed that CL-20 has been detonated at this time, and the detonation pressure of energetic co-crystals is the highest.

    Fig.3.Wavefront velocity of different models.

    Table 3Comparison between simulated value and calculateda value of shock wave velocity pressure.

    Fig.4.Wavefront velocity histogram of different models.

    3.2.Hotspot generation process

    Taking CL-20/GBL co-crystals as an example, the twodimensional temperature field diagram of shock induced supercells with void defect is drawn based on the energy equipartition theorem, and analyzed the void collapse and hot spot generation process at different timesteps, as shown in Fig.5.

    Fig.5.Hot spot evolution of CL-20/GBL co-crystals with defects shocked by(a)2 km/s and (b) 4 km/s.

    It can be seen from Fig.5 that in the CL-20/GBL co-crystals containing voids,impact of the shock wave intensifies the thermal motion of the atoms upstream of the cavity, accelerates the collision with the downstream atoms in the cavity to generate hot spots,and no mass accumulation occurs when the cavity is collapsed.Under low velocity of 2 km/s, the hot spot is slightly flat, and the recovery speed of the depression is slow after the wavefront sweeps through the cavity.Under high velocity of 4 km/s, the hot spot range became large and tapered,indicating that the increase in shock velocity can extend the size of the hot spot in the impact direction, which is consistent with the high-speed impact simulation results of RDX with voids[36].The central temperatures of hot spots formed by void collapse of four kinds of CL-20 cocrystals are slightly different, under the impact of 2 km/s, the maximum temperatures and pressures of energetic co-crystals and solvent co-crystals are both about 4050 K and 45.0 GPa under the impact of 4 km/s, the maximum temperatures of energetic cocrystals and solvent co-crystals are about 12,000 K, 95.8 GPa and 10,500 K, 92.2 GPa, respectively.In energetic co-crystals, the hot spot temperature of CL-20/TNT is slightly higher than 12,000 K,which is supposed to be related to the different chemical energy release of different ligand molecules.Under the same impact conditions, the void in energetic co-crystals delays the shockwave longer,and there is more energy exchange between the wavefront and the hot spot, which corresponds to the higher hot spot temperature of the energetic co-crystals above.In addition,the residual energy is less after the wavefront sweeps through the void, which makes it more difficult for downstream molecules to react than upstream.

    3.3.Quantity and distribution of main species

    This chapter investigates the relative content, number, and molecules fragmentation degree of species in the perfect and voidcontaining CL-20 co-crystals under the reflect velocity of 2 km/s and 4 km/s, as shown in Figs.6 and 7.The dotted line represents perfect supercells, and the solid line represents void-containing supercells.

    Fig.6.(a)Changes of the numbers of CL-20;(b)Ligand;(c)Molecules fragmentation degree;(d)NO2 in co-crystals under impact of 2 km/s(solid line is void containing co-crystals,dotted line is perfect co-crystals).

    Under the reflect velocity of 2 km/s, the number of CL-20 molecules and co-crystals ligands in the co-crystals containing defects decreases continuously.Molecule’s fragmentation degree almost remains unchanged before the cavity collapsed (0-1.76 ps), and then increases step by step.The corresponding curve in the perfect model changes smoothly, as shown in Fig.6.Nitrogen containing small species such as N2and NO2are generated by shock, and the number of small species such as H2O, CO2and HCN is very small,similar to the results of CL-20/TNT co-crystals simulation [37].Before the void collapsed in each supercell, the reaction speed of CL-20 and ligand is the same.However, after the void-containing supercell undergoes the void collapse stage, the reaction speed around the hot spot accelerates.In the void-containing supercells,the number of molecules involved in the reaction is significantly greater than that of the perfect supercells.Due to the compound effect of void collapse and reverse stretching, molecular fragmentation degree also shows an increasing trend, while the perfect model remains relatively flat in the above parameters.After the wavefront was swept,no CL-20 molecules in the supercell had been recovered to complete molecules.For perfect CL-20/DMC and voidcontaining CL-20/DMC, the reaction intensity and the molecular fragmentation degree change relatively smoothly, which is presumed to be caused by the stronger C-H…O hydrogen bond between molecules [11].

    The molecular reaction amount of energetic co-crystals components is much higher than that of solvent co-crystals.It can be seen that the shock response properties between energetic cocrystals and solvent co-crystals are quite different.The introduction of voids has a greater impact on the number of CL-20 and its ligand molecules,molecular fragmentation degree,and the number of products in the energetic co-crystals than solvent co-crystals.The molecular reaction of each component in the void-containing supercells is more intense than that of the perfect supercells,indicating that the hot spots caused by voids will significantly increase the sensitivity of the CL-20 co-crystals, especially the energetic co-crystals.The molecular fragmentation and the number of small species generated in CL-20 co-crystals are lower than those of pure CL-20 [38], indicating that co-crystals modification can significantly reduce sensitivity [5].At this shock velocity, CL-20 molecules and co-crystal ligands react significantly only when hot spots are generated,and there is no significant change in molecular fragmentation degree, with no detonation reaction occurring.

    Under the reflect velocity of 4 km/s, the molecular reaction amount, the number of product species and the molecular fragmentation of these supercells are significantly higher than those under 2 km/s impact, as shown in Fig.7.Before the void collapse(0-1.04 ps), the number of molecules in perfect and defective cocrystals is similar, but the reaction rate of atoms near the hot spot is much higher than that of other atoms after void collapse,resulting in a more dramatic change in the number of components of the void-containing co-crystals than in the perfect co-crystals.When the wavefront sweeps through the whole co-crystals, most of the CL-20 molecules react completely,on the contrary,the ligand molecules react incompletely, which is consistent with the simulation results of defect-containing CL-20/DNB [39].Due to the difference in the density of nitro functional groups, the reaction volume of solvent ligands is significantly lower than that of energetic ligands, and the amount of N2and NO2generated in the solvent co-crystals is also lower, which proves that the detonation performance and sensitivity of solvent co-crystals are lower than that of energetic co-crystals.The molecular number of CL-20/DMC supercell is the most smoothly under the shock of this velocity,preliminarily judged that the sensitivity reduction effect of CL-20/DMC co-crystals is the best among the four types of co-crystals,but the detonation performance is weak.In energetic co-crystals,CL-20/DNB’s fragmentation degree and the amount of N2species are slightly less than those in CL-20/TNT, and the number of residual CL-20 molecules and in CL-20/DNB is larger.It can be seen that the sensitivity of CL-20/DNB of the two co-crystals is lower,which is consistent with the conclusion of CL-20/DNB preparation and test results[39].

    Table 4Molecular distribution of co-crystals ligands under the velocity of 4 km/s.

    Table 4 (continued)

    The results show that the existence of voids has a great influence on the number of co-crystals ligand molecules.It is speculated that the molecular will dimerize first as the wavefront scanning under low-velocity reflects.However,the CL-20 reaction amount and the ligand molecular reaction amount in the void-containing supercells are relatively less than perfect supercells, speculated that the existence of voids causes the wavefront energy to concentrate on the cavity collapse, thus reducing the energy of wavefront and slows down the chemical reaction frequency downstream of the cavity.In order to verify this conjecture, in four kinds of CL-20 co-crystal shock models containing voids, the distribution of co-crystal molecules with wavefront positions upstream of the void,downstream of the void,and swept through the model were plotted,as shown in Table 4, corresponding to the number of co-crystals ligand molecules-time change curve in Fig.5(b).

    3.4.Chemical elementary reaction frequency

    The chemical elementary reaction analysis method was used to investigate the influence of voids on the chemical reaction of shock induced CL-20 co-crystals, and the shock response process was divided into two stages: cavity collapse and hot spot fading.The high-frequency elementary reactions and frequency of perfect and void-containing CL-20 co-crystals under shock velocity of 2 km/s and 4 km/s were investigated, as shown in Tables 5 and 6.

    The main reactions under the shock of 2 km/s include adsorption of CL-20 with ligands, nitro shedding and dimerization reactions,as shown in Table 5.In the cavity collapse stage(0.64-1.76 ps), the dimerization reaction between CL-20 and ligand occurred mainly,and the frequency of CL-20/GBL was the highest,but the net reaction frequency was low,indicating that the wavefront energy is close to the activation energy of the reaction.The solvent cocrystals CL-20/DMC and CL-20/GBL have a higher net frequency of denitration reactions, but the energetic ligands also undergo nitro shedding, which leads to similar NO2production in each supercell during the cavity collapse stage, corresponding to the result of the aforementioned NO2quantity change curve (see Fig.6(d)).In the hot spot fading stage (1.76-5.00 ps), various dimerization reactions were more frequent,CL-20/GBL still has the highest frequency, corresponding to the analysis of the tightest bond between the CL-20 and the ligand in chapter 3.2,and have the better desensitization effect.The number of small carboncontaining molecules generated is very small.It is speculated that CL-20 molecules have a weak ring-opening reaction under the impact of 2 km/s, and the dimerization reaction frequency is greater than the corresponding reverse reaction,meaning that the velocity cannot trigger the detonation of each CL-20 co-crystal[16].

    It is found that the chemical reaction of perfect co-crystals is generated during the wavefront scanning and sparse stretching[40], while the hot spot formed in the void-containing co-crystals,which increases the reaction rate around the hot spot,resulting in a significant difference in reaction frequency.Besides the frequency of dimerization reaction between CL-20 and ligand is higher in perfect co-crystals, indicating that there is little correlation with this reaction and the void.The generation reaction frequency of small species in the void-containing supercell is higher.It is speculated that this type of reaction is mainly caused by abnormally high-speed atoms impacting the hot spots downstream of the cavity.In the hot spot fading stage, the net frequency of nitro shedding of CL-20 in void-contained co-crystals is much greater than that in perfect co-crystals, of which CL-20/DNB shows the highest frequency, which is consistent with the above conclusion,the molecular reaction intensity of CL-20 and the production of N2and NO2in the void-contained co-crystals are higher than those in perfect co-crystals (seen in Fig.6).

    It can be seen from Table 6 that under the shock of 4 km/s, thepriority of nitro group shed from CL-20 molecules in the cavity collapse stage(0.48-1.04 ps)is much higher than that from ligand TNT and DNB molecules.In CL-20/DMC, the ligand lost its methyl group and self-decomposed, while the other ligands did not undergo high-frequency self-decomposition reaction.At this shock velocity, there is still a dimerization reaction between CL-20 molecules and ligand molecules, but with the exception of CL-20/TNT,the generation frequency of each dimer is lower than 2 km/s,indicating that the priority of the polymerization reaction is reduced,and the molecules tend to decompose.

    Table 5Reaction frequency of CL-20 co-crystals unit under shock velocity of 2 km/s.

    In the hot spot fading stage (1.04-4.00 ps), high-frequency reactions involve small species,generating H2O,HN2O,H2O2et al.It is considered that the detonation of each co-crystal has occurred under the impact of this velocity.The main reaction frequency of the perfect co-crystal is higher than that of the void-containing cocrystal, the gap between the two in the solvent co-crystal is large,and it is consistent with the above conclusion that the effect of void defects on the detonation performance of solvent co-crystal is higher than that of energetic co-crystal.The molecules in the energetic co-crystals CL-20/TNT and CL-20/DNB undergo polymer reactions and stack to form super-large carbon-containing groups.CL-20/TNT is the most prone to polymer reactions, while solvent co-crystals are more difficult for this reaction,and the existence of voids has little effect on the size of the super-large molecular clusters.The results of the two-dimensional temperature field diagram show that the void in CL-20/TNT delays the shockwave longer and has the most energy exchange.The reaction frequency statistics table shows that the small species reaction frequency of CL-20/TNT is greater than that of the other co-crystals.It is speculated that both are factors that affect the hot spot formation temperature.

    4.Discussion

    Here in two energetic eutectic CL-20/TNT, CL-20/DNB and two solvent co-crystals CL-20/DMC, CL-20/GBL perfect and defectcontaining model impact response are studied, the effect of void defects on the formation of hot spots,the type,and quantity of mainproducts and the reaction frequency of species were analyzed.Nomura et al.[36] simulated the RDX impact with defects and showed that the hot spot spreads in a cone shape when the shock velocity is high.The research in this paper shows that the hot spot of void-containing CL-20 co-crystals would be elongated along the shock direction,which is similar to Nomura’s conclusion.Liu et al.[15]simulated the piston impact of CL-20/TNT energetic co-crystal and showed that there is no detonation at 2 km/s impact, CL-20 polymerization and decomposition reactions exist simultaneously;under 4 km/s impact,the direct decomposition of perfect and void-containing co-crystals dominates, which is consistent with the product analysis and reaction frequency results in this research.It is speculated that CL-20 co-crystals with solvent components cannot produce super-large molecular clusters under the high-speed impact, and defects have little effect on the size of super-large molecular clusters generated.

    Table 6Reaction frequency of CL-20 co-crystals unit under shock velocity of 4 km/s.

    Here in the statistical results of the shock response frequency of each co-crystals model (see Table 6) show that CL-20 solvent cocrystal is extremely difficult to generate super-large molecular clusters at the initial stage of impact, and it is difficult to cause aggregation.Wang et al.[39] speculated that the CL-20/DNB cocrystal has lower sensitivity than CL-20/TNT in their co-crystal preparation experiments, here compares the number of small species products (see Figs.3 and 4) to verify this conclusion.It is also found that the effect of void defects on CL-20/TNT is greater than CL-20/DNB.When detonation occurs under the high-speed impact (4 km/s), there is little difference in the number of species and reaction frequency among these types of co-crystals, but the introduction of voids would significantly reduce the frequency of species generation of each co-crystal when the hot spot faded,and reduce the detonation performance of each CL-20 co-crystal, and more significant impact on the solvent co-crystals are observed.

    5.Conclusions

    ReaxFF-lg force field was performed to analyze the formation of shock-induced hot spots, the number of molecules, the main products,and the reaction frequency of CL-20/TNT,CL-20/DNB,CL-20/DMC, and CL-20/GBL co-crystals.The conclusions are listed as follows.

    (1) The hot spot evolution process shows that as the shock velocity increases,the size of the hot spot becomes longer along the shock direction,and the shape of the hot spot becomes a cone at 4 km/s.The different types of co-crystal ligands have a significant impact on the hot spot temperature.Under the same conditions, the void in the energetic co-crystal has a greater retardation effect on the wavefront,and there is more energy exchange between the wavefront and the void,resulting in a high hot spot temperature.

    (2) Under the reflect velocity of 2 km/s,the main reaction of CL-20 co-crystal is the synthesis of dimers and the shedding of nitro groups, without the formation of a large number of species and no detonation.Compared with perfect crystal,the existence of voids aggravates the chemical reactions,reduces the safety of CL-20 co-crystals,and the sensitivity of void-containing energetic co-crystal increases more obviously.

    (3) The sensitivity of CL-20/DNB in energetic co-crystals is lower than that of CL-20/TNT under the impact of 4 km/s, indicating that the sensitivity reduction effect is better.Energetic co-crystal is more prone to polymerization reaction than solvent co-crystals and would unaffected by voids, while solvent co-crystal is less sensitive and harder to agglomerate.Under the shock of this velocity, both the perfect and the hole-containing CL-20 co-crystals detonate, and the wavefront energy is reduced after sweeping through the cavity,reducing the frequency of chemical reactions downstream of the cavity and weakening the detonation performance.In addition, the effect of void defects on the detonation performance of CL-20 solvent co-crystal is higher than that of CL-20 energetic co-crystal.

    Declaration of competing interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    Acknowledgements

    This project was supported by the National Natural Science Foundation of China (22275018), and the Project of State Key Laboratory of Explosion Science and Technology (Beijing Institute of Technology) (Grant No.QNKT20-04).

    久久精品夜夜夜夜夜久久蜜豆| 日本免费一区二区三区高清不卡| 亚洲精品一区av在线观看| 久久精品影院6| 国模一区二区三区四区视频| 国产精品永久免费网站| 变态另类成人亚洲欧美熟女| 久久久成人免费电影| 老鸭窝网址在线观看| 热99re8久久精品国产| 人妻丰满熟妇av一区二区三区| 国产精品1区2区在线观看.| 亚洲精品一区av在线观看| 亚洲18禁久久av| 国产成人福利小说| 亚洲精品色激情综合| 一区福利在线观看| 国产精品香港三级国产av潘金莲| 国产精品香港三级国产av潘金莲| 亚洲18禁久久av| 我要搜黄色片| 亚洲五月婷婷丁香| 女警被强在线播放| 欧美3d第一页| 亚洲男人的天堂狠狠| 欧美日韩亚洲国产一区二区在线观看| 国产精品 欧美亚洲| 亚洲在线观看片| 国产激情偷乱视频一区二区| 99国产精品一区二区三区| 麻豆久久精品国产亚洲av| ponron亚洲| 一进一出好大好爽视频| 国产高清视频在线观看网站| 一二三四社区在线视频社区8| 搡老妇女老女人老熟妇| 欧美性感艳星| 欧美一区二区亚洲| 久久午夜亚洲精品久久| av专区在线播放| 国产免费一级a男人的天堂| 色吧在线观看| 精品一区二区三区视频在线 | 在线免费观看的www视频| 午夜老司机福利剧场| 色尼玛亚洲综合影院| 精品福利观看| 高清日韩中文字幕在线| 精品一区二区三区视频在线 | 免费观看人在逋| 国产亚洲欧美98| 国产熟女xx| 亚洲精品一区av在线观看| 听说在线观看完整版免费高清| 非洲黑人性xxxx精品又粗又长| 小蜜桃在线观看免费完整版高清| 五月伊人婷婷丁香| 成人无遮挡网站| 成人18禁在线播放| av福利片在线观看| 亚洲国产欧洲综合997久久,| bbb黄色大片| 真人做人爱边吃奶动态| 亚洲国产精品sss在线观看| 亚洲av一区综合| 欧美中文日本在线观看视频| 亚洲专区国产一区二区| 亚洲av成人av| 两人在一起打扑克的视频| 91九色精品人成在线观看| 亚洲真实伦在线观看| 亚洲av电影在线进入| 精品无人区乱码1区二区| 亚洲第一欧美日韩一区二区三区| 蜜桃久久精品国产亚洲av| 亚洲avbb在线观看| 精品久久久久久久久久免费视频| av专区在线播放| 女同久久另类99精品国产91| 久久中文看片网| 久久性视频一级片| 午夜免费男女啪啪视频观看 | 狂野欧美激情性xxxx| 国产aⅴ精品一区二区三区波| 国产综合懂色| 国产又黄又爽又无遮挡在线| 真人做人爱边吃奶动态| 亚洲无线观看免费| 男女做爰动态图高潮gif福利片| 欧美xxxx黑人xx丫x性爽| 尤物成人国产欧美一区二区三区| 91久久精品电影网| 制服丝袜大香蕉在线| 精品午夜福利视频在线观看一区| 国产精品久久久久久久电影 | 国产熟女xx| 在线视频色国产色| 女人十人毛片免费观看3o分钟| 欧洲精品卡2卡3卡4卡5卡区| 精品一区二区三区人妻视频| 叶爱在线成人免费视频播放| 99久国产av精品| 激情在线观看视频在线高清| 日韩欧美国产一区二区入口| 91麻豆精品激情在线观看国产| 免费电影在线观看免费观看| 国产免费av片在线观看野外av| 日本精品一区二区三区蜜桃| 两个人视频免费观看高清| 超碰av人人做人人爽久久 | 99久久综合精品五月天人人| 成人亚洲精品av一区二区| 国产激情欧美一区二区| 美女大奶头视频| 亚洲精品影视一区二区三区av| 中文字幕人妻熟人妻熟丝袜美 | 又爽又黄无遮挡网站| 欧美日韩瑟瑟在线播放| 久久久久久久久久黄片| 国产精品美女特级片免费视频播放器| 成人高潮视频无遮挡免费网站| 亚洲片人在线观看| 午夜精品久久久久久毛片777| 真人做人爱边吃奶动态| 亚洲片人在线观看| 一夜夜www| 久久久久九九精品影院| 很黄的视频免费| 午夜精品一区二区三区免费看| 两性午夜刺激爽爽歪歪视频在线观看| 久久精品91蜜桃| 波多野结衣高清无吗| 18禁国产床啪视频网站| 国产一区二区亚洲精品在线观看| 99国产综合亚洲精品| 极品教师在线免费播放| 久久久久国内视频| 51国产日韩欧美| 久久久久久久亚洲中文字幕 | 日本 av在线| 国产精华一区二区三区| 天堂av国产一区二区熟女人妻| 日本与韩国留学比较| 欧美激情在线99| 黑人欧美特级aaaaaa片| 国产男靠女视频免费网站| 69人妻影院| 久久久久精品国产欧美久久久| 免费在线观看影片大全网站| 亚洲aⅴ乱码一区二区在线播放| 国产在视频线在精品| 免费在线观看影片大全网站| 偷拍熟女少妇极品色| 亚洲国产精品久久男人天堂| 日韩免费av在线播放| 国产伦人伦偷精品视频| 国产高清三级在线| 国产一区在线观看成人免费| 成人av在线播放网站| 搞女人的毛片| 长腿黑丝高跟| 欧美日韩亚洲国产一区二区在线观看| 午夜福利高清视频| 啪啪无遮挡十八禁网站| 精品一区二区三区av网在线观看| 亚洲国产精品999在线| 久久精品91无色码中文字幕| 在线播放无遮挡| 婷婷六月久久综合丁香| 亚洲欧美激情综合另类| 夜夜看夜夜爽夜夜摸| 久久久精品欧美日韩精品| 性欧美人与动物交配| 亚洲国产日韩欧美精品在线观看 | 成年女人永久免费观看视频| 成人无遮挡网站| 99热这里只有精品一区| 久久国产精品影院| 天天添夜夜摸| 99久国产av精品| 国产精品电影一区二区三区| 久久久久国产精品人妻aⅴ院| av中文乱码字幕在线| 高清日韩中文字幕在线| 欧美一区二区国产精品久久精品| 最近视频中文字幕2019在线8| 欧美乱妇无乱码| 久久久国产成人精品二区| 成人精品一区二区免费| 免费看a级黄色片| a级毛片a级免费在线| 很黄的视频免费| 一本综合久久免费| 此物有八面人人有两片| 综合色av麻豆| 性欧美人与动物交配| 激情在线观看视频在线高清| 久久精品国产亚洲av香蕉五月| 青草久久国产| 看免费av毛片| 国产99白浆流出| 国产高清视频在线观看网站| 日日干狠狠操夜夜爽| 美女 人体艺术 gogo| 久久久久久久亚洲中文字幕 | 琪琪午夜伦伦电影理论片6080| 欧美极品一区二区三区四区| av天堂在线播放| 亚洲精华国产精华精| 日韩大尺度精品在线看网址| 亚洲中文字幕日韩| 神马国产精品三级电影在线观看| 他把我摸到了高潮在线观看| 男女视频在线观看网站免费| 免费看日本二区| 免费av观看视频| 一区二区三区激情视频| 美女大奶头视频| 女生性感内裤真人,穿戴方法视频| 一本久久中文字幕| 久久久久精品国产欧美久久久| 成人国产综合亚洲| 日韩 欧美 亚洲 中文字幕| tocl精华| 国产精品98久久久久久宅男小说| 一个人免费在线观看电影| 亚洲在线自拍视频| 啦啦啦免费观看视频1| 午夜精品久久久久久毛片777| 激情在线观看视频在线高清| 欧美一区二区亚洲| 丰满人妻熟妇乱又伦精品不卡| 国产毛片a区久久久久| 人人妻人人澡欧美一区二区| 亚洲熟妇中文字幕五十中出| 男插女下体视频免费在线播放| 国内久久婷婷六月综合欲色啪| 搡老熟女国产l中国老女人| 成人av一区二区三区在线看| 午夜福利欧美成人| 搡老妇女老女人老熟妇| 国产v大片淫在线免费观看| 99久久九九国产精品国产免费| 日本黄色视频三级网站网址| 国产麻豆成人av免费视频| 国产av麻豆久久久久久久| 国产高清视频在线播放一区| 色综合站精品国产| 久9热在线精品视频| 中文字幕高清在线视频| 国产精品综合久久久久久久免费| 老汉色av国产亚洲站长工具| 国产精品嫩草影院av在线观看 | 亚洲va日本ⅴa欧美va伊人久久| 97碰自拍视频| 变态另类丝袜制服| 美女 人体艺术 gogo| 欧美丝袜亚洲另类 | 国产精品三级大全| 人妻久久中文字幕网| 亚洲成人久久性| 搡老妇女老女人老熟妇| 久久久精品大字幕| av国产免费在线观看| 久久久色成人| 两性午夜刺激爽爽歪歪视频在线观看| 九九久久精品国产亚洲av麻豆| 最近在线观看免费完整版| 国产精品野战在线观看| 日本与韩国留学比较| 免费看十八禁软件| 国产探花在线观看一区二区| 亚洲 国产 在线| 深爱激情五月婷婷| 欧美一级毛片孕妇| 国产精品综合久久久久久久免费| 日韩欧美免费精品| 男女那种视频在线观看| 久久欧美精品欧美久久欧美| 老汉色av国产亚洲站长工具| 成人亚洲精品av一区二区| 国产亚洲精品久久久久久毛片| 丰满乱子伦码专区| 成人18禁在线播放| 人人妻,人人澡人人爽秒播| 亚洲aⅴ乱码一区二区在线播放| 两个人看的免费小视频| 久久香蕉精品热| 国产黄a三级三级三级人| 亚洲成av人片免费观看| 亚洲成av人片在线播放无| 欧美av亚洲av综合av国产av| 免费人成视频x8x8入口观看| 天堂影院成人在线观看| 成年女人看的毛片在线观看| 香蕉丝袜av| 国产久久久一区二区三区| 精品国产超薄肉色丝袜足j| 精品午夜福利视频在线观看一区| 九九在线视频观看精品| 免费看美女性在线毛片视频| 亚洲人与动物交配视频| 制服人妻中文乱码| 人妻夜夜爽99麻豆av| 久久久久久大精品| 亚洲av电影不卡..在线观看| 少妇人妻一区二区三区视频| 女人被狂操c到高潮| 91字幕亚洲| 男人和女人高潮做爰伦理| 波野结衣二区三区在线 | 高潮久久久久久久久久久不卡| 亚洲美女视频黄频| 99热这里只有精品一区| 在线视频色国产色| 男人和女人高潮做爰伦理| xxx96com| 黄片大片在线免费观看| 99久久精品国产亚洲精品| 亚洲国产日韩欧美精品在线观看 | 蜜桃久久精品国产亚洲av| 亚洲国产精品999在线| 老司机午夜十八禁免费视频| av中文乱码字幕在线| 国产麻豆成人av免费视频| 人妻夜夜爽99麻豆av| 禁无遮挡网站| 欧美三级亚洲精品| 午夜福利视频1000在线观看| 亚洲国产高清在线一区二区三| 久久天躁狠狠躁夜夜2o2o| 手机成人av网站| 一夜夜www| 国产成+人综合+亚洲专区| 午夜福利免费观看在线| 日韩欧美在线乱码| 听说在线观看完整版免费高清| 日韩大尺度精品在线看网址| 亚洲精品粉嫩美女一区| 一级毛片高清免费大全| 欧美性猛交黑人性爽| 国产精品电影一区二区三区| 久久久久久久午夜电影| 午夜日韩欧美国产| 欧美一区二区亚洲| 国产精品乱码一区二三区的特点| 69av精品久久久久久| 欧美日韩瑟瑟在线播放| 夜夜爽天天搞| 欧美最新免费一区二区三区 | 99精品久久久久人妻精品| 不卡一级毛片| 午夜精品久久久久久毛片777| av在线蜜桃| 给我免费播放毛片高清在线观看| 色哟哟哟哟哟哟| 成年人黄色毛片网站| 久久香蕉精品热| 99riav亚洲国产免费| 黄色成人免费大全| 桃红色精品国产亚洲av| 色哟哟哟哟哟哟| 真人一进一出gif抽搐免费| 精品国产亚洲在线| 99精品久久久久人妻精品| 国产高潮美女av| 精品熟女少妇八av免费久了| 久久欧美精品欧美久久欧美| 老司机午夜十八禁免费视频| 国产真实伦视频高清在线观看 | 国产精品美女特级片免费视频播放器| 嫩草影视91久久| bbb黄色大片| 欧美+亚洲+日韩+国产| 少妇人妻精品综合一区二区 | 精品国产美女av久久久久小说| 精品人妻一区二区三区麻豆 | 欧美成人免费av一区二区三区| 精品人妻1区二区| e午夜精品久久久久久久| 亚洲美女黄片视频| 欧美日韩亚洲国产一区二区在线观看| 国产精品久久久久久精品电影| 日本在线视频免费播放| 欧美黄色片欧美黄色片| 99精品在免费线老司机午夜| or卡值多少钱| 两性午夜刺激爽爽歪歪视频在线观看| 99久久无色码亚洲精品果冻| 小说图片视频综合网站| 精品欧美国产一区二区三| 每晚都被弄得嗷嗷叫到高潮| 国内精品久久久久久久电影| 亚洲自拍偷在线| 99久久综合精品五月天人人| 99久久精品热视频| 国产色爽女视频免费观看| 99热只有精品国产| 看片在线看免费视频| 精品久久久久久久毛片微露脸| 欧美日本视频| 日本a在线网址| 欧美激情久久久久久爽电影| 成人无遮挡网站| 亚洲 欧美 日韩 在线 免费| 可以在线观看的亚洲视频| 亚洲精品亚洲一区二区| 欧美丝袜亚洲另类 | 波多野结衣高清无吗| 国产视频内射| 欧美bdsm另类| 1024手机看黄色片| 婷婷精品国产亚洲av在线| 夜夜夜夜夜久久久久| 丁香欧美五月| 丰满人妻熟妇乱又伦精品不卡| 宅男免费午夜| 国产精品电影一区二区三区| 极品教师在线免费播放| 老熟妇乱子伦视频在线观看| 午夜福利视频1000在线观看| 亚洲精品在线美女| 非洲黑人性xxxx精品又粗又长| 日本与韩国留学比较| 99久久综合精品五月天人人| 高清日韩中文字幕在线| 国产欧美日韩精品一区二区| 香蕉丝袜av| 亚洲av五月六月丁香网| 青草久久国产| 日韩欧美国产一区二区入口| 欧美成人免费av一区二区三区| 久久精品国产亚洲av香蕉五月| 久久人妻av系列| 免费高清视频大片| 韩国av一区二区三区四区| a级毛片a级免费在线| 欧美乱码精品一区二区三区| 久久久国产精品麻豆| 亚洲欧美日韩高清在线视频| 国产精品av视频在线免费观看| 99热6这里只有精品| 老熟妇仑乱视频hdxx| www.www免费av| 天天一区二区日本电影三级| 国内毛片毛片毛片毛片毛片| 黄色丝袜av网址大全| 国产av在哪里看| 欧美色视频一区免费| 亚洲色图av天堂| 国产私拍福利视频在线观看| 国内精品久久久久精免费| 国产精品爽爽va在线观看网站| 伊人久久精品亚洲午夜| 午夜免费激情av| 12—13女人毛片做爰片一| 欧美色欧美亚洲另类二区| 啦啦啦免费观看视频1| 级片在线观看| 丰满乱子伦码专区| 国产精品一区二区三区四区免费观看 | x7x7x7水蜜桃| 欧美bdsm另类| 国产色婷婷99| 国产精品精品国产色婷婷| 亚洲av五月六月丁香网| 好男人在线观看高清免费视频| av天堂在线播放| 亚洲人成网站在线播| 日韩欧美三级三区| 国产亚洲精品av在线| 免费人成视频x8x8入口观看| 老汉色∧v一级毛片| 国产视频内射| 精品国产美女av久久久久小说| 美女被艹到高潮喷水动态| 制服丝袜大香蕉在线| 亚洲专区国产一区二区| 女警被强在线播放| 精华霜和精华液先用哪个| 国产 一区 欧美 日韩| 成人特级av手机在线观看| 伊人久久精品亚洲午夜| 国产亚洲精品久久久com| 日韩精品青青久久久久久| 亚洲第一欧美日韩一区二区三区| 亚洲,欧美精品.| 国产欧美日韩精品一区二区| 国产一区二区亚洲精品在线观看| 真人做人爱边吃奶动态| 90打野战视频偷拍视频| 看黄色毛片网站| 琪琪午夜伦伦电影理论片6080| 国产激情偷乱视频一区二区| 久久精品影院6| 99热这里只有是精品50| 亚洲精品色激情综合| 两个人看的免费小视频| 免费电影在线观看免费观看| av黄色大香蕉| 国产高清激情床上av| 法律面前人人平等表现在哪些方面| 色综合亚洲欧美另类图片| 国产av麻豆久久久久久久| 日本a在线网址| 午夜两性在线视频| xxxwww97欧美| 国内揄拍国产精品人妻在线| 国产成人啪精品午夜网站| 美女免费视频网站| 少妇的丰满在线观看| 少妇人妻一区二区三区视频| 少妇的丰满在线观看| av国产免费在线观看| 老司机福利观看| 午夜福利18| 午夜福利成人在线免费观看| 国产国拍精品亚洲av在线观看 | 国产高清视频在线观看网站| 一本综合久久免费| 国产亚洲精品久久久com| 日本与韩国留学比较| 亚洲天堂国产精品一区在线| a级毛片a级免费在线| 国产精品1区2区在线观看.| 一进一出抽搐动态| 九九久久精品国产亚洲av麻豆| 真人做人爱边吃奶动态| 欧美黑人巨大hd| 一进一出抽搐gif免费好疼| 婷婷丁香在线五月| 九色成人免费人妻av| 色老头精品视频在线观看| 久9热在线精品视频| 国产一区在线观看成人免费| 成人一区二区视频在线观看| 日韩av在线大香蕉| 国产一级毛片七仙女欲春2| 久久久色成人| 成人国产一区最新在线观看| 可以在线观看的亚洲视频| av福利片在线观看| 国产99白浆流出| 中亚洲国语对白在线视频| 婷婷精品国产亚洲av| 一个人观看的视频www高清免费观看| 婷婷亚洲欧美| 69av精品久久久久久| 黄色日韩在线| 十八禁人妻一区二区| 国产一级毛片七仙女欲春2| 亚洲成av人片免费观看| 不卡一级毛片| 1024手机看黄色片| 嫩草影视91久久| 国产亚洲欧美98| 亚洲精品成人久久久久久| 最新在线观看一区二区三区| 波多野结衣高清无吗| 欧美色欧美亚洲另类二区| 无遮挡黄片免费观看| 精品日产1卡2卡| 亚洲激情在线av| 免费观看精品视频网站| 19禁男女啪啪无遮挡网站| 午夜日韩欧美国产| 搡老妇女老女人老熟妇| 亚洲专区国产一区二区| 给我免费播放毛片高清在线观看| 999久久久精品免费观看国产| 亚洲国产色片| 国产精品影院久久| 黑人欧美特级aaaaaa片| 天美传媒精品一区二区| 麻豆久久精品国产亚洲av| 亚洲男人的天堂狠狠| 久久久久久大精品| 日本成人三级电影网站| 国产单亲对白刺激| 亚洲精品粉嫩美女一区| 精品福利观看| 夜夜夜夜夜久久久久| 国产欧美日韩一区二区三| 国产精品亚洲美女久久久| 精品一区二区三区视频在线 | 亚洲成av人片免费观看| 一a级毛片在线观看| 小蜜桃在线观看免费完整版高清| 99riav亚洲国产免费| 丁香欧美五月| 波多野结衣巨乳人妻| 高清毛片免费观看视频网站| 国产乱人视频| 91麻豆av在线| 九九久久精品国产亚洲av麻豆| 免费观看精品视频网站| 黄色视频,在线免费观看| 国产精品一区二区三区四区免费观看 | 99国产精品一区二区蜜桃av| 99久久久亚洲精品蜜臀av| 国产精品久久电影中文字幕| 欧美性猛交╳xxx乱大交人| 熟妇人妻久久中文字幕3abv| 夜夜爽天天搞| 久久精品91蜜桃| 18禁黄网站禁片午夜丰满| 亚洲国产精品sss在线观看| 亚洲国产精品久久男人天堂| 午夜福利免费观看在线| 国产午夜福利久久久久久| 蜜桃久久精品国产亚洲av| 99国产综合亚洲精品| 一区二区三区高清视频在线| 高清在线国产一区| 男女之事视频高清在线观看|