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

    Flame behavior,shock wave,and instantaneous thermal field generated by unconfined vapor-liquid propylene oxide/air cloud detonation

    2023-07-31 13:30:14ConglingYeQingleiDuLijunLiuQiZhng
    Defence Technology 2023年7期

    Cong-ling Ye ,Qing-lei Du ,Li-jun Liu ,Qi Zhng ,*

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

    b School of Safety Science and Emergency Management, Wuhan University of Technology, Wuhan, Hubei, 430070, China

    Keywords:Numerical simulation Fuel-air explosive Fuel/air cloud Flame Shock wave Thermal field

    ABSTRACT Energy output and heating effects are essential for vapor-liquid fuel/air cloud detonation in the fuel-air explosive(FAE)applications or explosion accidents.The purpose of this study is to examine the dynamic large-size flame behavior,shock wave propagation law,and instantaneous thermal field generated by unconfined vapor-liquid propylene oxide (PO)/air cloud detonation.Based on computational fluid dynamics(CFD)and combustion theory,a numerical simulation is used to study the detonation process of a PO/air cloud produced by a double-event fuel-air explosive (DEFAE) of 2.16 kg.The large-scale flame behavior is characterized.The flame initially spreads radially and laterally in a wing shape.Subsequently,the developed flame increases with a larger aspect ratio.Moreover,the propagation laws of shock waves at different heights are discussed.The peak pressure of 1.3 m height level with a stepwise decline is obviously different from that of the ground with an amplitude of reversed'N'shape.In the vast majority of the first 6.9 m,the destructive effect of the shock wave near the ground is greater than that of the shock wave at 1.3 m height.Furthermore,the dynamic instantaneous isothermal field is demonstrated.The scaling relationship of various isotherms in the instantaneous thermal field with the flame and initial cloud is summarized.The comprehensive numerical model used in this study can be applied to determine the overpressure and temperature distribution in the entire fuel/air cloud detonation field,providing guidance for assessing the extent of damage caused by DEFAE detonation.

    1.Introduction

    Burning and explosion of fuel/air mixture clouds [1-8] are common safety accidents in the industrial production process,such as vapor cloud explosions caused by the leakage of liquefied petroleum gas [9-13],gas explosions induced by the break-up of natural gas transportation pipelines [14-17] etc.Moreover,the energy output and heating effect of vapor-liquid fuel cloud detonations have many military and engineering applications.Propylene oxide(PO) [18] is a typical liquid organic energetic fuel that is commonly used in fuel air explosives (FAE)and rocket propellants[19].Double-event FAE(DEFAE),which uses PO as the main fuel,is a new type of high-energy explosives.After the center charge explodes,the PO fuel is forced to disperse and ultimately forms a PO/air mixed cloud covering a large area.Explosive shock waves and thermal damage of the vapor-liquid PO/air cloud are important research issues.

    In recent decades,many scholars have carried out in-depth analyses and discussions of the rapid reaction mechanism of PO using various experimental methods.Saxena et al.studied the dissociation mechanism of PO,which was examined using Laser-Schlieren measurements in incident shock waves over 32-717 Torr and 1429-1936 K using 5%acetone diluted in krypton[20].Ramalingam et al.developed a detailed kinetic mechanism to describe both the oxidation reaction and pyrolysis chemistry of PO[21].Li et al.conducted an experiment in a vertical detonation tube to study the detonation parameters of fuel PO drops dispersed in air.The results showed that the relationship of the critical initiation energy with the equivalence ratio was a U-shaped curve,and the critical initiation energy of fuel-air was the least,whereas the equivalence of PO was 1.05 [22].The above studies are mainly concentrated in confined spaces such as shock tubes for studying the explosion parameters of PO because it is difficult,inefficient and costly to conduct large-scale experimental research on the explosion of fuel/air mixtures under unconstrained conditions,and the theoretical calculation is also very complicated.Therefore,there are relatively few studies on the detonation of vapor-liquid clouds formed with PO as fuel in an unconstrained space.Numerical simulation[23-25]is an economical and effective tool for starting large-scale explosive experiments that are difficult to perform in reality.

    The detonation pressure of PO/air mixtures has been extensively studied.Du et al.simulated a second explosion induced by PO aerosols,and its overpressures were much higher than those of the first explosion.The peak overpressure approached 2.6 MPa[26].Bai et al.numerically investigated the shockwave overpressure fields near the ground of multisource FAE explosions [27].Liu et al.studied the deflagration-to-detonation transition (DDT) of PO/air mixtures in an experimental tube.The pressure histories at different points along the tube and the trajectory of the pressure wave during the DDT process were obtained and discussed [28].Hasson et al.also studied the DDT process and found good agreement between the measured pressures and C-J calculations [29].Liu et al.explored detonation pressures at various concentrations in a confined cylindrical tube [30].It is known that the detonation behavior in a confined space easily forms a triple point owing to shock wave reflection at the wall.Thus,the detonation overpressure value is so large that it has no reference value for the study of the radial detonation pressure of a fuel/air cloud in an unconstrained open space.Therefore,it is necessary to compare and analyze the pressures near the ground and in the radial direction where the fuel distribution density is the highest,and discover their propagation laws.

    For the small-size flame behavior of various fuels,including PO fuel [31],experiments and numerical simulations have been conducted.Burluka et al.studied the laminar flames of three C3H6O isomers and found that most of these flames exhibited a nonlinear dependency of flame speed on the stretch rate [32].Ji et al.presented an experimental investigation on the flame radiation characteristics of a pool fire,including flame emissivity,diameter and temperature[33].The steady laminar diffusion flamelet model was applied by Lemos et al.to model the kinetics of a diffusive jet flame,where the fuel was a mixture of natural gas,CO2and H2.Turbulence effects have been incorporated in laminar flamelets through probability density functions [34].The structure of a developing fireball upon ignition of two-phase pentane release was presented[35].Knyazkov et al.also reported the experimental and model characterization of the flame structure of PO under stoichiometric and fuel-rich conditions at atmospheric pressure.The mole fractions of the species in three premixed flames stabilized on a flatflame burner were quantitatively measured using flame sampling molecular beam mass spectrometry[36].Gong et al.measured and compared the laminar flame speeds of C3oxygenated fuels (npropanol,propanal,and acetone)and hydrocarbons(propane)in a combustion bomb [37].However,there are few studies on the characterization of large-sized flame development.Most of these studies used fireballs measured using high-speed photography and infrared imaging to characterize the explosive flame [38-40].However,optical measurements are prone to interference from soot and fumes,and most of the measured temperatures are the surface temperatures of the fireball.The temperature distribution inside the flame is difficult to obtain,and the ability to distinguish the details of the temperature differences is too poor to capture the flame development images.

    The thermal damage of fuel/air cloud detonation is not only reflected in the lethal effect of intense heat conduction and convection inside the flame,but also in the damage caused by thermal radiation outside the flame[41,42].Accordingly,it is not reasonable to only use the fireball size to represent detonation temperature field.The thermal effect of explosive products on the target is an unsteady heat transfer process that is affected by many factors,and its measurement has been a difficult problem in various countries.

    In the above context,based on fluid dynamics and combustion theory,this study focuses on the detonation mechanism of the vapor-liquid PO/air cloud under high-energy ignition using the computational fluid dynamics (CFD) simulation tool.The CFD software used for the numerical simulation is Fluent.In addition,the concept of the instantaneous isotherm model of the fuel/air cloud detonation field is proposed by integrating the three heat transfer modes to reflect the dynamic characteristics of the transient temperature field quantitatively and evaluate the entire range of high-temperature fields,including thermal radiation,more reasonably.As a method for evaluating the thermal damage of the detonation field,the dynamic instantaneous isotherm model can better simulate how the isotherm coverage changes with detonation time and reasonably evaluate the consequences of thermal radiation.Theoretically,a detonation temperature field close to reality can be obtained.To explore the detonation of vapor-liquid PO/air mixture cloud and its heating effect in open space,three issues are mainly emphasized:1)the characterization of large-size flame behavior;2)the propagation laws of shock waves at different heights;and 3)the dynamic instantaneous isothermal field of cloud detonation and its scaling relationship with the flame and initial cloud.

    2.Numerical simulation

    2.1.Numerical methods

    Numerical models of the detonation process of liquid fuel are very important for the simulation of such a process.Advances in computational fluid mechanics have provided a basis for further insights into the dynamics of fuel flow and combustion theory.In this study,the continuous phase model solves compressible Navier-Stokes equations using the SIMPLE algorithm and the standardk-ε turbulent model.The discrete phase model (DPM) and discrete ordinate (DO) radiation model were adopted to simulate the dynamics of PO droplets and the radiation of the PO/air cloud detonation temperature field,respectively.

    (1) The continuity equation is

    where,tis the detonation time(s),ρ is the fluid density(kg/m3),is the fluid velocity (m/s),andSmis the mass transferred from the dispersed second phase to the continuous phase.In this study,Smis the mass from the vaporization of the PO droplets (kg).

    (2) The momentum conservation equation in an inertial reference frame is

    (3) The equation of energy conservation is

    where,the first three items on the right represent energy transport caused by heat conduction,species diffusion and viscous dissipation,respectively,Shis the source item including chemical reaction energy and radiation energy(W/m3),Eis the total energy of PO fuel(J/kg),his the sensible enthalpy(J/kg),λ is the thermal conductivity which is expressed as a function of temperature,λtis the turbulent thermal conductivity (W/m·K),is the diffusion flux ofj'species(mol/(m2·s)).

    (4) Species transport model

    To solve the conservation of chemical species,the local mass fraction of each speciesYj’was predicted through the solution of a convection-diffusion equation for thej'th species.Therefore,the species transport equation takes the following general form

    (5) Turbulence model

    For turbulence simulation,the Reynolds-averaged approach is a non-direct numerical simulation method.The Reynolds-averaged Navier-Stokes (RANS) equations govern the transport of averaged flow quantities,with the entire range of turbulence scales being modeled.Here,considering our computing resources and model applicability,the well-developed standard κ-ε model [43] is used which is a semi-empirical model based on transport equations for turbulence kinetic energy and its dissipation rate.It has been widely used to simulate complex turbulent processes by scholars[44-46].

    The turbulence kinetic energy κ and its dissipation rate ε were obtained from the following transport equations:

    The turbulent viscosity,ut,is computed by combining κ and ε as follows

    where,Gkrepresents the generation of turbulence kinetic energy owing to the mean velocity gradients,Gbis the generation of turbulence kinetic energy owing to buoyancy,andYMrepresents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate.The constants,andCμare 1.44,1.92 and 0.09,respectively.σkand σεare the turbulent Prandtl numbers for κ and ε,namely,1.0 and 1.3,respectively.is a coefficient associated with buoyancy,which is taken as 0 and 1 when the main flow is parallel and vertical to the gravity direction,respectively.

    (6) Discrete phase model

    The motion trajectories of micron-sized PO droplets can be predicted using the DPM and stochastic tracking model (STM) by integrating the force balance on the particle.As the discrete and continuous phases are interact with each other,a two-way coupling model is used to solve the two-phase equations.The exchanges of mass,heat,and momentum from the PO droplets to the continuous phase are qualitatively described in Fig.1.

    Fig.1.Mass,heat and momentum transfer between the discrete PO droplet particle and continuous gas phase.

    So,the mass,momentum and heat exchange are respectively computed as

    where,is the initial mass flow rate of the PO droplet (kg/s),is the initial mass of droplet particles(kg),CDis drag coefficient,Reis the relative Reynolds number,Ppis density of the particle(kg/m3),dpis diameter of the particle (m),upis the particle velocity (m/s),ufis the gas velocity(m/s),is the mass of the droplet particle on the cell entry(kg),mp,outis the mass of the droplet particle on the cell exit(kg),cppis the heat capacity of the droplet particle (J/kg·K),cgpis the heat capacity of the gas product species (J/kg·K),Tbpis the boiling point temperature of PO (K),Trefis the reference temperature of the droplet particle for enthalpy(K),Tp,inis the temperature of the droplet particle on cell entry(K),Tp,outis the temperature of the droplet particle on cell exit (K),Hpyris the heat of pyrolysis as volatiles are evolved(J/kg).

    (7) DO radiation model

    The DO model spans the entire range of optical thicknesses.The computational cost is moderate for typical angular discretization,and the memory requirements are modest.The DO model considers the radiative transfer equation (RTE);thus,its equation is written as

    Fig.2 illustrates the process of radiative heat transfer.

    Fig.2.Schematic diagram of radiative heat transfer.

    The model includes gas-phase species such as combustion products (H2O,CO2);thus,the absorption and emission of the gas are significant.The radiative properties of the products,such as absorption coefficients,mainly refer to Ju·s work[47],as shown in Table 1.

    Table 1 Absorption coefficients of vapor and carbon dioxide.

    (8) Chemical reaction kinetics model

    The combustion of the PO/air mixture can be simplified as a one-step reaction.The chemical reaction of gaseous PO/air with a stoichiometric volume fraction is

    The kinetic rate constant for reaction,kr,is computed using Arrhenius expression

    where,Aris the pre-exponential factor,βris the temperature exponent,Eris the activation energy for the reaction(J/kmol),Rris the universal gas constant(J/kmol·K).

    2.2.Model construction and setting of boundary conditions

    A static experiment of 2.16 kg DEFAE 1.3 m above the ground is served as a research object.The DEFAE bomb and fuel/air cloud generated in three-dimensional (3D) space are symmetric.Thus,the entire flow field is simplified as a two-dimensional (2D) rotational axisymmetric model that transforms a 3D problem into a 2D one,and its boundary conditions before diffusion are as shown in Fig.3(a),the ambient temperature is 300 K and the atmospheric pressure is 0.1 MPa.The fuel/air distribution before ignition is based on fully dispersed fuel/air cloud field at 100 ms in the dispersion process,as shown in Fig.3(b).The ignition position is right below the DEFAE warhead near the ground as shown in Fig.3(c).The initial simulation parameters of the PO are obtained from Liu·s work [46].

    Fig.3.Computational model of PO/air cloud detonation generated by 2.16 kg FAE: (a) Computational domain and mesh generation;(b) Fuel/air cloud before ignition;(c) Ignition position.

    Because the fuel disperses rapidly in the early stage owing to the central explosion,and the detonation process is mainly concentrated in the fuel/air cloud center,a local grid refinement method around the central axis is presented.As the velocity of the cloud front decreases,the expectation of the refined grid resolution in the far field is decreased.The same is true for pressure.Thus,the entire computing domain is meshed in divergent unstructured grids,and three grid sizes are adopted to facilitate the comparison,as shown in Table 2.

    Table 2 Mesh setting.

    2.3.Verification of numerical model

    2.3.1.Verification of cloud size

    Fig.4 shows that the simulated and experimental fuel/air cloud fields have a good coincidence in shape and size at 40 ms in the fuel dispersion process.At 100 ms after the center explosion,which drives fuel to disperse,the radii and heights of the PO/air clouds simulated by refined,moderate and coarse mesh types are in agree with the experimental results,as shown in Fig.5.The errors between the three simulations with different grids and the experiment were 3.16%,5.79%,and 14.2%in the radius of the PO/air cloud,and 6.15%,7.7%,and 14.2% in the height of PO/air cloud,which indicates that the finer mesh density can simulate more accurately.However,finer grids do not significantly improve accuracy.Therefore,to improve computational efficiency,a divergent unstructured moderate grid is adopted.In the fuel dispersal simulation using a moderate grid,the cloud radius approaches 4.02 m when it is detonated,as shown in Fig.5.

    Fig.4.Comparison of simulated and experimental fuel/air cloud fields at 40 ms in the dispersion process (a) Experimental fuel/air cloud field;(b) Simulated fuel/air cloud field.

    Fig.5.Comparison of simulated fuel/air cloud sizes in various mesh types and experiment·s result.

    There is a bit difference between the experimental and numerical simulation results at 0 m from the comparison of simulated and experimental fuel/air cloud fields at 40 ms in the dispersion process.It is because the experimental image taken form a long distance is a side projection;however,the simulated picture is the middle sectional view since the entire flow field is simplified as a two-dimensional (2D) rotational axisymmetric model.At the center of the cloud,a nearly fuel-free zone is formed by the central explosive drive,as shown in Fig.4(b)

    2.3.2.Verification of temperature and pressure

    After the PO/air cloud is sufficiently dispersed,its radius approaches 4.02 m.At this time which is the start time of detonation process,the cloud is detonated by the high-energy trinitrotoluene(TNT)explosive.Fig.6(a)shows the temperature field that cuts off 625 K below at 6 ms in the detonation process.The simulated results of line-1 and line-2 are in good agreement with the trend of the temperatures in the transverse and axial directions obtained from the experimental results,as shown in Fig.6(b).

    Fig.6.Verification of temperature: (a) Temperature field at 6 ms in detonation simulation;(b) Comparison of temperature between simulation and experiment.

    For pressure verification,the overpressure curves of a static explosion experiment of 5.16 kg DEFAE are used for comparison with the simulation results,as shown in Fig.7(a) and Fig.7(b).In Fig.7(c),P1,P2,P3,and P4 represent the monitoring sites which are 2 m,4 m,6 m,and 8 m from the cloud center respectively.It can be seen that the first two pressure curves of the simulated shock wave are consistent with the experimental pressures at 2 m and 4 m.The simulated peak overpressures at 6 m and 8 m are similar to the measured values,and the simulated pressure curves are slightly flat.

    Fig.7.Verification of pressure: (a) Layout of FAE device;(b) Fireball generated by fuel/air cloud detonation;(c) Comparison between numerical simulation and experimental overpressure of 5.16 kg DEFAE.

    3.Results and discussion

    3.1.Flame behavior

    As can be seen from the flame development of the PO/air mixture cloud detonation in 10 ms in Fig.8,the flame grows most rapidly and the chemical reaction rate is the fastest within 3 ms after TNT initiation.At 1.25 ms,the flame front approaches a maximum spread rate of 2710 m/s,as seen in Fig.9.Meanwhile,the dynamic process of the flame shape variation can be observed at every transient moment.The flame in the wing shape mainly extends radially on both sides of the radial direction,and the radius of the flame coverage approaches 4.95 m at 3 ms.The flame radius is the distance of the farthest edge of the flame from the central axis,which is also called flame coverage.According to the temporal position of the flame front and flame spread rate in Fig.9,in 3-5 ms,the flame has a secondary acceleration development process in which the flame acceleration is much smaller than that in the first 3 ms.The secondary acceleration of flame development is attributed to the increase in the thermal expansion rate of the detonation products,which pushes the heat flow to continue to expand outward.In addition,the flame front velocity has a slight downward trend in the first 2 ms.This is because in the first 2 ms or so,the cloud detonation flame is still in its initial development stage,and the flame front is extremely unstable.The overpressure and flame produced by the high-energy initiation are in contact with the cloud field in a hemispherical shape.When the cloud is detonated,the flame spreads outward in the radial direction.There are certain differences in the concentrations of liquid and gaseous fuel and the sizes of fuel droplet radially at different heights,which may also affect the propagation of detonation flame.The reaction of fuel droplet with large size is slower than that of small-size droplet and gaseous fuel.Although the flame speed decreases at a certain stage,the thermal expansion effect of products accelerates the flame again due to the continuous cloud detonation until the end of the cloud detonation.At 10 ms,the flame radius is close to 6 m,which is 1.49 times the initial cloud size.During the subsequent 15-47.5 ms,the flame formed by the PO/air cloud detonation is uplifted,its aspect ratio increases,and the flame changes from a macroscopic flat shape to a more complex messy shape,as shown in Fig.10.

    Fig.8.Flame development of PO/air mixture cloud detonation.

    Fig.9.Temporal position of flame front and flame spread rate in 10 ms.

    Fig.10.Flame development of PO/air mixture cloud detonation in 15-47.5 ms.

    According to Hasegawa and Sato·s model [48] which is applicable to liquid fuel explosions with fuel masses less than 10 kg and fireball temperatures of approximately 3600 K,the fireball diameter is calculated asD=5.25M0.314,whereMis the fuel mass(kg).The diameter of explosion fireball produced by 2.16 kg of the undispersed pure PO fuel is about 6.68 m,which is 0.54 times of the maximum flame diameter induced by DEFAE detonation with the same mass of PO fuel.This embodies the advantages of DEFAE.Owing to the initial dispersion effect,PO fuel is more widely distributed with air mixing in space;therefore,the flame coverage generated by detonation is larger.

    Flame development is closely related to the rate of chemical reaction.Fig.11 shows that the reaction rate of PO in the cloud field approached its peak at 1 m “r-n" means the reaction rate atnm from the cloud center on the 1.3 m height level.The combustion rate outside the cloud is still large at 4.2 m,because the disturbance and thermal expansion caused by detonation cause the fuel to move out of the cloud field.Meanwhile,in Fig.12(a),it is observed that there is a small fluctuation in the combustion rate from 4 to 4.5 m,which is because the explosive flame(as shown in Fig.12(b))sucks up trace amounts of residual gaseous PO fuel (as shown in Fig.12(c)) and oxygen in the air (as shown in Fig.12(d)),and the flame edge becomes unstable.In addition,it is also found that the combustion rates oscillate between 1 m and 3.5 m.This is because the concentration distribution of PO fuel is not uniform in the cloud.Thus,the fuel PO and oxygen do not necessarily react in accordance with the stoichiometric ratio,and the reaction rate differs because of the fuel concentration.Moreover,some of the heat absorbed by the liquid-phase PO as it transforms into gaseous PO also causes a variation in the reaction rate.

    Fig.11.Reaction rate of PO of various monitoring points in the cloud at 1.3 m height level.

    Fig.12.Fluctuation in combustion rate from 4 to 4.5 m at around 3.75 ms and its corresponding phenomenon:(a)Reaction rate of PO of various monitoring points outside the cloud at 1.3-m height level;(b) Flame generated by cloud detonation at 3.75 ms;(c) Residual gaseous fluel at 3.75 ms;(d) Oxygen distribution at 3.75 ms.

    The chemical reaction range is 1.12 times of the initial premixed PO/air cloud.Compared with many combustion and explosion experiments[49-51]in closed spaces,the range of chemical reaction fields outside the premixed fuel/air zone in this simulation experiment is much smaller,which is due to the cloud detonation induced by the high-energy initiation of TNT,and the reaction rate in the fuel concentration area is extremely fast.As shown in Fig.13,the combustion reaction time at 2 m and 3 m in the radial direction of 1.3 m horizontal height are 0.07 ms and 0.13 ms,respectively.In Fig.13,P(1.3,2)and P(1.3,3)mean the point which is 2 m and 3 m away from the cloud center on the 1.3 m height level.Only a small amount of fuel at the edge of the cloud diffused outward because most of the fuel is consumed inside the cloud.

    Fig.13.Molar concentration profiles of oxygen and products at the 1.3 m height level: (a) Molar concentration at P(1.3,2);(b) Molar concentration at P(1.3,3).

    3.2.Shock wave propagation law

    The variation in the peak overpressure of a shock wave with a specific distance near the ground is discussed by two explosion experiments of 8 kg and 100 kg in Liu·s work [52].The specific distance isR/M1/3,whereRis the radial distance from the center of the fuel/air cloud andMis the mass of the FEA.To verify the 2.16 kg FAE simulation in this study,Liu's work and the static explosion experiment of 5.16 kg FAE are used to conduct a comparative analysis,as shown in Fig.14.The trend of 2.16 kg FAE simulation result is consistent with those of the experiments.

    Fig.14.Peak pressure of FAE explosion varies with specific distance.

    As can be seen from Fig.15,at 0.75 ms,the flame and shock wave front formed by the high energy of the TNT explosion detonated the PO/air cloud,and the pressure front was no longer orderly and expanded more rapidly in the radial direction than in the axial direction.This was similar to the flame development described above.According to the C-J theory,the detonation front is regarded as an ideal strong discontinuity surface without thickness,which is the overlap of the pressure wave front and flame front.When they pass through,the original explosives are immediately transformed into detonation products and release chemical energy.PO/air cloud detonation is assumed to be a typical characteristic of one-wave and two-zone structures.During 0-1.75 ms,the shock wave propagates much faster at the 1.3 m height level than near the ground,as shown in Fig.16.This is because there was little fuel near the ground.On the other hand,the fuel PO at the 1.3 m height level forms a detonation state owing to its adequate concentration.

    Fig.16.Temporal position of shock wave front and its velocity at two radial horizontal planes.

    As can be seen from the pressure diagrams at 1 ms,1.5 ms,2 ms and 3 ms in Fig.15,although the shock wave front near the ground propagates slowly at the beginning,it subsequently generates stronger shock wave points by the superposition of the incident and reflected shock waves of the PO/air cloud detonation.Therefore,the shock wave front near the ground speeds up continuously,and eventually the radial and axial pressure wave fronts form a regular semicircle and expand outwards from 4 to 10 ms,as shown in Fig.17.The shock wave has a faster velocity than that at the 1.3 m height level from 1.8 to 6 ms,shown in Fig.16.

    Fig.17.Pressure field of PO/air mixture cloud detonation during 4-10 ms.

    Furthermore,it can be seen from the pressure curves of all monitoring points at a horizontal line of 1.3 m height in Fig.18 that the pressure within the cloud range is much larger than that outside the cloud,indicating that the release of the cloud detonation energy is highly concentrated.According to the peak pressure curve in the red dashed line in Fig.18,PO/air cloud detonation can be divided into three stages along the radial path of the shock wave propagation.The peak pressures in the first stage are 16.8-22.2 bar,which is mainly influenced by the high-energy initiation mode of TNT in the first 1 m.The PO concentration distribution within 1.5-3.5 m is the most abundant,forming a relatively stable selfsustaining detonation stage,and the peak pressures are 11.92-12.82 bar.Thereafter,the peak pressures at the edge of the cloud decayed rapidly,and the shock wave in the air outside the 4.5 m day freely,as shown by the pressure curves at the shock wave attenuation stage in Fig.18.

    Fig.18.Pressure propagation at the 1.3 m height level.

    From the pressure analysis near the ground,it can be seen that under the influence of the high-energy initiation mode within the first 1 m,the first pressure peaks at 0.5 m and 1 m are both very large,approximately 21 bar,as shown in Fig.19(a).To avoid this interference,the second peak values of the pressure curves at the two monitoring points are selected as the pressure peak values generated by cloud detonation.It can be seen that the peak pressures generated by cloud detonation near the ground show an increasing trend in the first 3 m represented by the red dotted dashed line in Fig.19(a).In particular,the pressure rising rate increases sharply from 2 to 3 m,which may be due to the stable detonation formed by sufficient PO fuel at 1.3 m height right above.Consequently,the peak pressure of 24 bar at 3 m near the ground approaches the maximum value,which is close to the theoretical CJ detonation pressure of PO.Subsequently,the shock wave attenuates with increasing distance,as shown in Fig.19(b).The peak pressure decays rapidly at the cloud edge.

    Fig.19.Pressure propagation near the ground: (a) Pressure curves near the ground;(b) Pressure attenuation outside the cloud.

    A comprehensive comparison of the peak pressures on the ground and 1.3 m height shows that the two peak pressure curves within the first 1.5 m and outside of 6.9 m affected by the TNT ignition source is close.In the range of 2.5-6.9 m,the peak pressures on the ground are always greater than that at the 1.3 m height,as shown in Fig.20.The peak pressures on the ground show a large amplitude of reversed 'N' shape within the cloud,which is caused by the combined action of high-energy initiation pressure attenuation of the TNT and the superposition of cloud detonation pressure reflection and incidence.The peak pressure curve along the horizontal line of 1.3 m height shows 12.26 bar at 3.5 m,4.65 bar at 4 m,and 1.94 bar at 4.5 m.It can be seen that the pressure attenuation at the edge of the cloud field is severe,because the fuel concentration at 4 m is already extremely thin,which is not enough to maintain the development of detonation.At 4.5 m,it transforms into a shockwave propagation in air.In conclusion,when a fuel/air cloud explodes near the ground,it is better to evaluate the damage effect of the shock wave by the pressure at the ground.

    Fig.20.Comparison of peak pressures along the 1.3 m height level and near the ground.

    3.3.Thermal field

    The dynamic instantaneous thermal field generated by the PO/air cloud detonation is shown in Fig.21.Restricted by the ground,the entire temperature field expands to the open space in a semielliptical form.In the first 2 ms,each isotherm presents a significant wave shape because PO/air mixtures with different concentrations in the cloud have different heat absorption,light shading effects,and refractive indices of light.The region with a higher fuel concentration has a better heat absorption effect,and shows the characteristics of the wave trough on the isotherm in the heat transfer direction.For example,as shown in Fig.21,during 0.75-1 ms,the whole isotherms appear to be wavy lines because the PO/air cloud absorbs the heat from the flame and radiation in both the axial and radial directions.The PO/air cloud has better heat absorption in the axial direction because the fuel PO is radially dispersed and there is more abundant fuel in the radial direction.For example,at a level of 1.3 m height,the isotherms of 400-800 K still have obvious wavy shapes at 2 ms,however,these isotherms gradually turn into a circular arc axially.

    Fig.21.Instantaneous thermal field containing heat radiation.

    The temperature in the far field is primarily affected by the thermal radiation of the flame in the explosion field.With the exhaustion of the fuel,there is no shelter in the heat transfer direction,and the isotherm is gradually smoothed.The coverage of each isothermal line also varies with flame development.From a flat ellipse at 10 ms,the isotherm coverage shape develops into a semi-circle at 40 ms,and then turns into a slender ellipse with a larger aspect ratio at 60 ms.

    It is not reasonable to simply use the fireball size to represent the detonation temperature field.Accordingly,this study proposes the use of a quantified isotherm instead of fireball size to evaluate the range of influence of the thermal field produced by fuel/air cloud detonation.By comparing the initial PO/air cloud with a radius of 4.02 m in Fig.22(a),flame with a radius of 6 m at 10 ms in Fig.22(b)and the instantaneous thermal field considering radiation at 10 ms in Fig.22(c),the real range of thermal damage is far larger than the initial cloud size and flame coverage.The area enclosed by the 400 K isotherm at 10 ms approaches 209.9 m2.In Fig.22(c),it can be observed that the entire thermal field diverges to the ambient air with the cloud as the center,presenting an ellipsoid shape with decreasing temperature.The temperature is ultra-high within 10 m,and the area enclosed by the isothermal line at 1400 K is 35.1 m2,which is 16.7% of the area enclosed by the isotherm at 400 K.To facilitate analysis and reading of data of the coverage range of isotherms in the transient thermal field at 10 ms,the size of each isothermal line and its relationship with the flame and initial PO/air cloud are obtained as shown in Fig.23 by summarizing the data of isotherms in Fig.22(c).For instance,at 10 ms,the coverage range of the 400 K isotherm reaches 22.56 m;the coverage size of the isotherm below 1800 K is always greater than the initial cloud size;the coverage radius of the isotherm of 1183 K is 6 m which is equal to the flame radius.

    Fig.22.Comparison of initial PO/air cloud size,flame coverage,and instantaneous thermal field considering radiation at 10 ms: (a) Fuel/air cloud before ignition;(b)Flame profile at 10 ms after ignition;(c) Transient temperature field on at 10 ms after ignition.

    Fig.23.Analysis on the size of isothermal line in the instantaneous thermal field at 10 ms and its relation to flame and initial fuel/air cloud.

    4.Conclusions

    Based on CFD and combustion theory,a numerical method to simulate the entire process of vapor-liquid PO/air cloud detonation generated by 2.16 kg DEFAE is presented systematically to understand the dynamic large-size flame behavior,shock wave propagation law,and instantaneous thermal field.The key findings are as follows.

    (1) The maximum coverage of detonation flame and the chemical reaction zone of PO are respectively 1.49 and 1.12 times of the initial cloud size.At the early stage of detonation,the flame spreads radially and laterally in a wing shape.Subsequently,the developed flame rises higher with a larger aspect ratio and more complex shape,which is different from the initial flat outline.

    (2) As for shock wave,the peak pressure of 1.3 m height level with a stepwise decline is significantly different from that of the ground with an amplitude of reversed'N'shape.The PO/air cloud generates a self-sustaining stable detonation with peak pressures of 11.92-12.82 bar in the well-fueled area on a level with 1.3 m height.The maximum peak pressure on the ground was 24 bar at 3 m.Furthermore,in the vast majority of the first 6.9 m,the overpressure damage of the shock waves on the ground is greater than that level with 1.3 m height.Ultimately,when the fuel/air cloud explodes near the ground,it is better to evaluate the damage effect of the shock wave on the far field by ground pressure.

    (3) Furthermore,the novel concept of the instantaneous isotherm model in the detonation field of a fuel/air cloud is proposed to more reasonably evaluate the scope of the entire dynamic thermal field,including the heat radiation,which quantitatively reveals the thermal damage range of each isotherm and its scale relationship with the initial cloud and flame.For instance,at 10 ms,the coverage of 400 K isotherm reaches 22.56 m which is 2.8 times of the initial cloud.In addition,the coverage radii of these isotherms below 1800 K are always greater than the initial cloud radius.In particular,the coverage radius of the 1183 K isotherm is 6 m which is equal to that of the flame.

    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.

    Acknowledgments

    The research presented in this paper was supported by the National Natural Science Foundation of China (Grant No.11972089).

    丝袜喷水一区| 欧美日韩亚洲国产一区二区在线观看 | 亚洲国产精品999| 男人舔女人的私密视频| 日韩 欧美 亚洲 中文字幕| 97在线人人人人妻| 国产一区二区三区综合在线观看| 国产成人影院久久av| 日韩一卡2卡3卡4卡2021年| 亚洲国产精品一区三区| 亚洲三区欧美一区| 亚洲成人国产一区在线观看 | 侵犯人妻中文字幕一二三四区| 国产野战对白在线观看| 久久精品人人爽人人爽视色| 美女主播在线视频| 国产1区2区3区精品| 另类精品久久| 国产成人一区二区在线| 亚洲熟女精品中文字幕| 啦啦啦在线免费观看视频4| 亚洲精品成人av观看孕妇| 欧美97在线视频| 少妇裸体淫交视频免费看高清 | 又粗又硬又长又爽又黄的视频| 咕卡用的链子| 亚洲图色成人| 亚洲九九香蕉| 婷婷丁香在线五月| 一本大道久久a久久精品| av线在线观看网站| 免费观看人在逋| 少妇人妻久久综合中文| 美女福利国产在线| 精品亚洲乱码少妇综合久久| 亚洲欧美精品自产自拍| 91麻豆av在线| 一边亲一边摸免费视频| 王馨瑶露胸无遮挡在线观看| 国产亚洲av片在线观看秒播厂| 亚洲成色77777| av福利片在线| 精品卡一卡二卡四卡免费| 午夜两性在线视频| 99国产精品免费福利视频| 久久99精品国语久久久| 亚洲精品久久午夜乱码| 丝瓜视频免费看黄片| 岛国毛片在线播放| 国产成人91sexporn| 校园人妻丝袜中文字幕| 啦啦啦啦在线视频资源| 岛国毛片在线播放| 操出白浆在线播放| 一本久久精品| 男人添女人高潮全过程视频| 久久天躁狠狠躁夜夜2o2o | 女人久久www免费人成看片| 国产色视频综合| 国产片内射在线| 亚洲欧美一区二区三区国产| 亚洲视频免费观看视频| 亚洲五月婷婷丁香| 国产1区2区3区精品| 欧美+亚洲+日韩+国产| 国产精品.久久久| 日韩伦理黄色片| 啦啦啦在线免费观看视频4| 三上悠亚av全集在线观看| 成年人免费黄色播放视频| 欧美精品一区二区免费开放| 9色porny在线观看| 午夜免费成人在线视频| 最近最新中文字幕大全免费视频 | 黄片小视频在线播放| 汤姆久久久久久久影院中文字幕| 日本av手机在线免费观看| 免费看av在线观看网站| 色网站视频免费| 七月丁香在线播放| 97精品久久久久久久久久精品| 人妻 亚洲 视频| 久久天躁狠狠躁夜夜2o2o | 成人亚洲欧美一区二区av| 极品人妻少妇av视频| 午夜91福利影院| 新久久久久国产一级毛片| 99精国产麻豆久久婷婷| 欧美另类一区| 少妇粗大呻吟视频| 国产有黄有色有爽视频| 亚洲欧洲日产国产| 国产精品国产三级国产专区5o| 国产成人一区二区在线| 免费在线观看黄色视频的| 丝袜美足系列| 亚洲国产看品久久| 国产精品偷伦视频观看了| 久久人人爽av亚洲精品天堂| 乱人伦中国视频| 天天添夜夜摸| 国产成人啪精品午夜网站| 日韩 亚洲 欧美在线| 热99久久久久精品小说推荐| 少妇人妻 视频| 99九九在线精品视频| 新久久久久国产一级毛片| 精品免费久久久久久久清纯 | 亚洲精品美女久久av网站| 人人妻人人澡人人看| 精品国产国语对白av| 精品亚洲成国产av| 麻豆av在线久日| 中国国产av一级| 少妇粗大呻吟视频| 高清黄色对白视频在线免费看| 男女边吃奶边做爰视频| 亚洲国产成人一精品久久久| 国产熟女欧美一区二区| av片东京热男人的天堂| 久热爱精品视频在线9| 蜜桃在线观看..| 丝袜人妻中文字幕| 久久精品国产综合久久久| 久久国产精品男人的天堂亚洲| 永久免费av网站大全| 欧美国产精品一级二级三级| 国产爽快片一区二区三区| 成人亚洲欧美一区二区av| 国产男人的电影天堂91| 免费看av在线观看网站| 国产在线视频一区二区| 久久久久国产一级毛片高清牌| a级毛片黄视频| 18禁国产床啪视频网站| 菩萨蛮人人尽说江南好唐韦庄| 亚洲欧美清纯卡通| 久久亚洲国产成人精品v| 亚洲天堂av无毛| 午夜福利视频精品| 午夜免费观看性视频| 嫁个100分男人电影在线观看 | 亚洲欧美精品自产自拍| 91老司机精品| 国产亚洲精品第一综合不卡| 久久精品成人免费网站| 青青草视频在线视频观看| 国产精品久久久久成人av| 国产精品偷伦视频观看了| 国产国语露脸激情在线看| 成人国语在线视频| 青春草亚洲视频在线观看| 侵犯人妻中文字幕一二三四区| 国产91精品成人一区二区三区 | 一级黄色大片毛片| 日韩欧美一区视频在线观看| 免费女性裸体啪啪无遮挡网站| 操美女的视频在线观看| 美女福利国产在线| 在线亚洲精品国产二区图片欧美| 久久精品国产亚洲av涩爱| 最新的欧美精品一区二区| 亚洲精品国产一区二区精华液| 老司机亚洲免费影院| 日本91视频免费播放| 久久精品人人爽人人爽视色| 精品国产乱码久久久久久小说| 午夜免费鲁丝| 欧美亚洲日本最大视频资源| 精品少妇黑人巨大在线播放| 国产精品成人在线| 青青草视频在线视频观看| 91字幕亚洲| 女性生殖器流出的白浆| 女警被强在线播放| 色播在线永久视频| 亚洲国产毛片av蜜桃av| 男女之事视频高清在线观看 | 欧美精品亚洲一区二区| 午夜福利一区二区在线看| 性色av乱码一区二区三区2| 亚洲欧洲国产日韩| 一级毛片女人18水好多 | 欧美xxⅹ黑人| 国产极品粉嫩免费观看在线| 搡老乐熟女国产| 欧美 亚洲 国产 日韩一| 天天躁夜夜躁狠狠躁躁| 尾随美女入室| 中国美女看黄片| 中文精品一卡2卡3卡4更新| 日韩免费高清中文字幕av| 18禁观看日本| 岛国毛片在线播放| 欧美亚洲 丝袜 人妻 在线| 91麻豆av在线| 天天躁日日躁夜夜躁夜夜| 丁香六月天网| 婷婷色麻豆天堂久久| 精品高清国产在线一区| 国产91精品成人一区二区三区 | 日本午夜av视频| 日本欧美国产在线视频| 国产精品一区二区在线观看99| 免费久久久久久久精品成人欧美视频| 97在线人人人人妻| av福利片在线| 久久人人爽人人片av| 90打野战视频偷拍视频| 两个人看的免费小视频| 美女中出高潮动态图| 美女主播在线视频| 亚洲成人免费电影在线观看 | 狂野欧美激情性bbbbbb| 色网站视频免费| 制服人妻中文乱码| 精品一品国产午夜福利视频| 国产精品亚洲av一区麻豆| 亚洲精品国产区一区二| 视频在线观看一区二区三区| 超碰成人久久| videosex国产| 黄片播放在线免费| 成年动漫av网址| 亚洲人成电影免费在线| 国产亚洲一区二区精品| 欧美在线一区亚洲| 欧美人与性动交α欧美软件| videosex国产| 校园人妻丝袜中文字幕| 热99久久久久精品小说推荐| 自拍欧美九色日韩亚洲蝌蚪91| 久久久久久久大尺度免费视频| 亚洲国产精品一区二区三区在线| 精品少妇内射三级| 久久鲁丝午夜福利片| 国产亚洲欧美精品永久| 亚洲国产精品国产精品| 亚洲一卡2卡3卡4卡5卡精品中文| 少妇 在线观看| 99国产综合亚洲精品| 国产日韩一区二区三区精品不卡| 精品一品国产午夜福利视频| 日韩中文字幕视频在线看片| 丝瓜视频免费看黄片| 一区二区三区精品91| 免费在线观看视频国产中文字幕亚洲 | 国产欧美亚洲国产| 人人澡人人妻人| 国产精品久久久av美女十八| 亚洲国产av影院在线观看| 亚洲精品第二区| 香蕉丝袜av| 脱女人内裤的视频| 男女床上黄色一级片免费看| 一级毛片电影观看| 巨乳人妻的诱惑在线观看| av福利片在线| 美女视频免费永久观看网站| 一区二区三区精品91| 又粗又硬又长又爽又黄的视频| 亚洲精品国产av蜜桃| 亚洲一区二区三区欧美精品| 欧美国产精品va在线观看不卡| 精品少妇黑人巨大在线播放| 久久热在线av| av视频免费观看在线观看| 国产成人91sexporn| 婷婷成人精品国产| 大香蕉久久成人网| 国产黄色视频一区二区在线观看| 久久精品国产亚洲av高清一级| 亚洲成人手机| 视频区欧美日本亚洲| 黄片播放在线免费| 免费在线观看完整版高清| 2018国产大陆天天弄谢| 久久精品国产a三级三级三级| 久久九九热精品免费| 欧美精品啪啪一区二区三区 | 日韩人妻精品一区2区三区| 美国免费a级毛片| 国产一区二区三区综合在线观看| 欧美日韩视频高清一区二区三区二| h视频一区二区三区| 少妇人妻久久综合中文| 久久精品国产a三级三级三级| 日本av免费视频播放| a 毛片基地| 香蕉丝袜av| 真人做人爱边吃奶动态| 欧美日本中文国产一区发布| 9色porny在线观看| 久久av网站| 国产精品免费大片| 国产xxxxx性猛交| av片东京热男人的天堂| 日本五十路高清| 在线观看免费午夜福利视频| 亚洲精品中文字幕在线视频| 成年人黄色毛片网站| 水蜜桃什么品种好| 天天操日日干夜夜撸| 我要看黄色一级片免费的| 国产一区二区激情短视频 | 婷婷成人精品国产| 国产精品一国产av| 国产一级毛片在线| 制服诱惑二区| 久久久精品区二区三区| 免费高清在线观看视频在线观看| 亚洲成人免费电影在线观看 | 成年av动漫网址| 制服诱惑二区| 悠悠久久av| 国产成人影院久久av| 黄色视频在线播放观看不卡| 电影成人av| 免费在线观看视频国产中文字幕亚洲 | 亚洲国产精品一区二区三区在线| 日本欧美国产在线视频| 亚洲国产精品一区二区三区在线| 老司机亚洲免费影院| 国产有黄有色有爽视频| 亚洲自偷自拍图片 自拍| 黑丝袜美女国产一区| 国产成人免费无遮挡视频| 捣出白浆h1v1| 啦啦啦在线免费观看视频4| 亚洲,欧美,日韩| 国产高清国产精品国产三级| 亚洲人成网站在线观看播放| 欧美另类一区| 亚洲视频免费观看视频| 久久热在线av| 99国产精品免费福利视频| 亚洲国产精品一区三区| 欧美黑人欧美精品刺激| 久久综合国产亚洲精品| 老汉色av国产亚洲站长工具| 人人妻人人添人人爽欧美一区卜| 人人妻人人爽人人添夜夜欢视频| 久久人人爽av亚洲精品天堂| 久久久久精品人妻al黑| 日日爽夜夜爽网站| 国产精品.久久久| 啦啦啦在线免费观看视频4| 国产欧美日韩精品亚洲av| 精品一区二区三区av网在线观看 | 国产精品国产三级国产专区5o| 日韩av在线免费看完整版不卡| 欧美日韩av久久| 亚洲av欧美aⅴ国产| a级毛片黄视频| 99国产精品一区二区蜜桃av | 如日韩欧美国产精品一区二区三区| 亚洲中文字幕日韩| 亚洲一区二区三区欧美精品| 天堂俺去俺来也www色官网| 久久久久久久久久久久大奶| 久热爱精品视频在线9| 日韩人妻精品一区2区三区| 水蜜桃什么品种好| 国产精品国产三级专区第一集| 99热全是精品| 日韩av在线免费看完整版不卡| 97精品久久久久久久久久精品| 日韩一本色道免费dvd| 狂野欧美激情性bbbbbb| 亚洲五月婷婷丁香| 每晚都被弄得嗷嗷叫到高潮| 最黄视频免费看| 中文字幕制服av| 日本一区二区免费在线视频| 黑人巨大精品欧美一区二区蜜桃| 精品久久久久久电影网| 国产精品麻豆人妻色哟哟久久| 少妇的丰满在线观看| 在线观看一区二区三区激情| 又粗又硬又长又爽又黄的视频| 久久狼人影院| 久久ye,这里只有精品| 1024香蕉在线观看| 黄色a级毛片大全视频| 国产成人一区二区在线| 黑人猛操日本美女一级片| 涩涩av久久男人的天堂| 人成视频在线观看免费观看| 欧美97在线视频| 永久免费av网站大全| 午夜免费鲁丝| 亚洲 国产 在线| 欧美日韩一级在线毛片| 国产无遮挡羞羞视频在线观看| 免费女性裸体啪啪无遮挡网站| 高潮久久久久久久久久久不卡| avwww免费| 1024香蕉在线观看| 大型av网站在线播放| 亚洲欧美清纯卡通| 亚洲成色77777| 啦啦啦视频在线资源免费观看| tube8黄色片| 国产黄色免费在线视频| 一边摸一边抽搐一进一出视频| 女性生殖器流出的白浆| 久久ye,这里只有精品| 欧美精品高潮呻吟av久久| 别揉我奶头~嗯~啊~动态视频 | 老司机在亚洲福利影院| 丝袜美腿诱惑在线| 大片免费播放器 马上看| 18禁国产床啪视频网站| 亚洲综合色网址| 韩国精品一区二区三区| 香蕉国产在线看| 一边摸一边抽搐一进一出视频| 老司机影院毛片| 精品一区二区三区av网在线观看 | 男女高潮啪啪啪动态图| 国产亚洲av片在线观看秒播厂| 99精品久久久久人妻精品| 97精品久久久久久久久久精品| 天堂中文最新版在线下载| 中文字幕高清在线视频| 亚洲第一青青草原| 2018国产大陆天天弄谢| 国产欧美日韩一区二区三 | 久久精品人人爽人人爽视色| 亚洲欧美日韩高清在线视频 | 一区二区三区乱码不卡18| 人妻一区二区av| 国产激情久久老熟女| 国产亚洲av片在线观看秒播厂| 亚洲人成电影免费在线| 日韩一区二区三区影片| 最近最新中文字幕大全免费视频 | 亚洲国产日韩一区二区| 欧美大码av| 在线观看免费视频网站a站| 亚洲成国产人片在线观看| 国产深夜福利视频在线观看| 欧美日韩视频精品一区| 国产视频一区二区在线看| 51午夜福利影视在线观看| 久久精品国产综合久久久| 岛国毛片在线播放| 精品少妇久久久久久888优播| 热99久久久久精品小说推荐| 亚洲色图综合在线观看| 大陆偷拍与自拍| 下体分泌物呈黄色| 亚洲色图 男人天堂 中文字幕| 国产深夜福利视频在线观看| 一边亲一边摸免费视频| 亚洲久久久国产精品| 精品少妇内射三级| 人成视频在线观看免费观看| 蜜桃国产av成人99| 日韩大片免费观看网站| 免费看不卡的av| 亚洲,一卡二卡三卡| 黑丝袜美女国产一区| 国产精品久久久久久精品电影小说| 久久99精品国语久久久| 乱人伦中国视频| 亚洲成国产人片在线观看| av视频免费观看在线观看| 久久精品久久精品一区二区三区| 夜夜骑夜夜射夜夜干| 免费观看av网站的网址| 国产日韩欧美视频二区| 久久这里只有精品19| 青春草视频在线免费观看| 两个人看的免费小视频| 黑丝袜美女国产一区| 黄色视频在线播放观看不卡| 观看av在线不卡| 无遮挡黄片免费观看| 久久久久精品国产欧美久久久 | cao死你这个sao货| 国产亚洲一区二区精品| 亚洲av男天堂| 1024香蕉在线观看| 好男人视频免费观看在线| 亚洲av欧美aⅴ国产| 老司机靠b影院| 免费不卡黄色视频| 久热这里只有精品99| 国产99久久九九免费精品| 中文字幕另类日韩欧美亚洲嫩草| 久久国产精品影院| 99热国产这里只有精品6| 天天躁狠狠躁夜夜躁狠狠躁| 丰满迷人的少妇在线观看| 亚洲人成77777在线视频| 成年人午夜在线观看视频| 香蕉国产在线看| 91成人精品电影| 亚洲欧洲精品一区二区精品久久久| 女人久久www免费人成看片| 国产精品久久久久成人av| 欧美亚洲 丝袜 人妻 在线| 成人国产一区最新在线观看 | 国产一区二区三区综合在线观看| 欧美 亚洲 国产 日韩一| 在线观看免费日韩欧美大片| videos熟女内射| 欧美少妇被猛烈插入视频| 美女高潮到喷水免费观看| 色婷婷久久久亚洲欧美| 麻豆国产av国片精品| 午夜影院在线不卡| 黄色怎么调成土黄色| 性色av乱码一区二区三区2| 国产精品一二三区在线看| 99久久精品国产亚洲精品| 国产精品国产av在线观看| 在线观看免费高清a一片| 亚洲av日韩在线播放| 老汉色∧v一级毛片| 精品久久蜜臀av无| 精品人妻熟女毛片av久久网站| 人妻人人澡人人爽人人| 国产黄频视频在线观看| 精品人妻一区二区三区麻豆| 亚洲av电影在线观看一区二区三区| 国产免费现黄频在线看| 97精品久久久久久久久久精品| 午夜老司机福利片| 国产欧美日韩一区二区三 | 久久国产精品影院| 老司机午夜十八禁免费视频| 91麻豆精品激情在线观看国产 | 欧美 日韩 精品 国产| 成人三级做爰电影| 九草在线视频观看| 水蜜桃什么品种好| 菩萨蛮人人尽说江南好唐韦庄| 免费看不卡的av| 咕卡用的链子| 中文欧美无线码| 午夜福利视频在线观看免费| 我的亚洲天堂| 国产成人精品久久二区二区91| 亚洲精品中文字幕在线视频| 啦啦啦在线观看免费高清www| 99re6热这里在线精品视频| 热99国产精品久久久久久7| 欧美精品一区二区免费开放| 国产又爽黄色视频| 黄色怎么调成土黄色| 久久99一区二区三区| 老司机靠b影院| 一边亲一边摸免费视频| videosex国产| 日韩一区二区三区影片| 各种免费的搞黄视频| 亚洲,欧美,日韩| 在线 av 中文字幕| 午夜老司机福利片| a级片在线免费高清观看视频| 男男h啪啪无遮挡| 亚洲av国产av综合av卡| 97在线人人人人妻| 我的亚洲天堂| 欧美在线一区亚洲| 老熟女久久久| 久久久久久久久免费视频了| 免费女性裸体啪啪无遮挡网站| 悠悠久久av| 夫妻午夜视频| 精品久久久久久久毛片微露脸 | 国产野战对白在线观看| 色视频在线一区二区三区| 国产成人精品久久二区二区91| 男人舔女人的私密视频| 国产高清国产精品国产三级| 少妇裸体淫交视频免费看高清 | 精品人妻在线不人妻| 亚洲国产最新在线播放| 亚洲色图 男人天堂 中文字幕| 精品视频人人做人人爽| 精品卡一卡二卡四卡免费| 青草久久国产| 日韩精品免费视频一区二区三区| 好男人视频免费观看在线| 欧美精品一区二区免费开放| 精品视频人人做人人爽| 99久久99久久久精品蜜桃| 亚洲美女黄色视频免费看| 日韩精品免费视频一区二区三区| 99久久99久久久精品蜜桃| 香蕉国产在线看| 亚洲精品久久成人aⅴ小说| 91九色精品人成在线观看| 一区在线观看完整版| av欧美777| 精品卡一卡二卡四卡免费| 香蕉国产在线看| 久热这里只有精品99| 国产在视频线精品| 精品熟女少妇八av免费久了| 亚洲视频免费观看视频| 老司机亚洲免费影院| 国产精品秋霞免费鲁丝片| 欧美人与善性xxx| 水蜜桃什么品种好| 日韩大片免费观看网站| 国产在线免费精品| 黄网站色视频无遮挡免费观看| 国产亚洲欧美精品永久| 成人亚洲欧美一区二区av| 免费在线观看影片大全网站 | 国产精品香港三级国产av潘金莲 |