• <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).

    av线在线观看网站| 亚洲精品国产色婷婷电影| 亚洲精品中文字幕一二三四区 | 少妇粗大呻吟视频| 精品国产乱子伦一区二区三区 | 黄网站色视频无遮挡免费观看| 成人18禁高潮啪啪吃奶动态图| 日韩有码中文字幕| 亚洲精品中文字幕一二三四区 | 欧美午夜高清在线| 日韩制服丝袜自拍偷拍| 亚洲国产精品成人久久小说| 91老司机精品| 国产精品亚洲av一区麻豆| 亚洲av国产av综合av卡| 精品久久久久久久毛片微露脸 | 亚洲性夜色夜夜综合| 欧美黑人精品巨大| 黄色 视频免费看| 精品少妇黑人巨大在线播放| 91老司机精品| 亚洲免费av在线视频| 午夜福利在线观看吧| 考比视频在线观看| 欧美日韩亚洲高清精品| 一本大道久久a久久精品| 国产片内射在线| 黄色片一级片一级黄色片| 国产福利在线免费观看视频| 国产在线视频一区二区| 美女午夜性视频免费| 国产一级毛片在线| 国产精品熟女久久久久浪| 午夜免费鲁丝| 亚洲国产欧美在线一区| 久9热在线精品视频| 亚洲av日韩在线播放| 老汉色av国产亚洲站长工具| 日本欧美视频一区| 免费在线观看影片大全网站| 欧美精品啪啪一区二区三区 | 下体分泌物呈黄色| 日日摸夜夜添夜夜添小说| 正在播放国产对白刺激| av国产精品久久久久影院| 免费久久久久久久精品成人欧美视频| 一区二区三区精品91| 久久中文字幕一级| 99re6热这里在线精品视频| 国产xxxxx性猛交| 国产欧美日韩一区二区三区在线| 午夜激情av网站| 天天躁日日躁夜夜躁夜夜| 久久国产精品人妻蜜桃| 精品欧美一区二区三区在线| 精品国产乱码久久久久久小说| 电影成人av| 亚洲性夜色夜夜综合| 91精品三级在线观看| 成年动漫av网址| 国产精品影院久久| 激情视频va一区二区三区| 一边摸一边抽搐一进一出视频| 亚洲av成人一区二区三| 国产成人av激情在线播放| 国产激情久久老熟女| 桃花免费在线播放| 久久香蕉激情| 亚洲人成电影观看| 亚洲一码二码三码区别大吗| 成人免费观看视频高清| 亚洲精品一二三| 美女中出高潮动态图| 黑人巨大精品欧美一区二区蜜桃| 日韩一区二区三区影片| 免费日韩欧美在线观看| 精品高清国产在线一区| 黑人欧美特级aaaaaa片| av网站在线播放免费| 99久久精品国产亚洲精品| 一区在线观看完整版| 色精品久久人妻99蜜桃| 香蕉丝袜av| 90打野战视频偷拍视频| 国产精品二区激情视频| 欧美变态另类bdsm刘玥| 国产又爽黄色视频| 女人高潮潮喷娇喘18禁视频| 亚洲国产欧美日韩在线播放| 99热全是精品| 色老头精品视频在线观看| 亚洲欧美一区二区三区久久| 久久久久久久精品精品| 免费久久久久久久精品成人欧美视频| av片东京热男人的天堂| 国产欧美日韩综合在线一区二区| 亚洲五月婷婷丁香| 一级,二级,三级黄色视频| 狠狠婷婷综合久久久久久88av| 亚洲中文av在线| 亚洲国产日韩一区二区| 黄片小视频在线播放| 午夜免费观看性视频| 免费少妇av软件| 天天躁狠狠躁夜夜躁狠狠躁| 国产深夜福利视频在线观看| 成人av一区二区三区在线看 | 老鸭窝网址在线观看| 久久人妻熟女aⅴ| 中文字幕另类日韩欧美亚洲嫩草| 久久久国产成人免费| 超色免费av| 久久免费观看电影| 精品福利永久在线观看| 黄频高清免费视频| 中文字幕另类日韩欧美亚洲嫩草| 久久国产精品大桥未久av| 女人久久www免费人成看片| www.999成人在线观看| kizo精华| 国产精品国产av在线观看| 日韩电影二区| 久久国产精品男人的天堂亚洲| 国产亚洲精品第一综合不卡| 精品亚洲乱码少妇综合久久| 好男人电影高清在线观看| 老熟妇仑乱视频hdxx| 欧美亚洲日本最大视频资源| 啦啦啦啦在线视频资源| 99热全是精品| 国产精品国产三级国产专区5o| 少妇 在线观看| 日韩 亚洲 欧美在线| av免费在线观看网站| a 毛片基地| 老熟妇乱子伦视频在线观看 | 国产野战对白在线观看| 国产亚洲精品第一综合不卡| 在线精品无人区一区二区三| 51午夜福利影视在线观看| 欧美激情高清一区二区三区| 日韩欧美免费精品| 欧美激情高清一区二区三区| 国产伦人伦偷精品视频| 亚洲九九香蕉| 久久中文字幕一级| 男女之事视频高清在线观看| 亚洲av日韩在线播放| 亚洲av欧美aⅴ国产| av超薄肉色丝袜交足视频| 中国美女看黄片| 国产精品偷伦视频观看了| 成人黄色视频免费在线看| 久久香蕉激情| 国产成人av教育| 日本vs欧美在线观看视频| 日本撒尿小便嘘嘘汇集6| 18禁裸乳无遮挡动漫免费视频| 在线av久久热| 精品久久蜜臀av无| 2018国产大陆天天弄谢| 国产成人啪精品午夜网站| 夫妻午夜视频| av天堂久久9| 精品少妇久久久久久888优播| 日本av手机在线免费观看| 淫妇啪啪啪对白视频 | 午夜福利影视在线免费观看| 一区二区三区乱码不卡18| 王馨瑶露胸无遮挡在线观看| 午夜91福利影院| 人妻人人澡人人爽人人| 精品高清国产在线一区| 亚洲av片天天在线观看| 99九九在线精品视频| 又黄又粗又硬又大视频| 欧美 亚洲 国产 日韩一| 我要看黄色一级片免费的| 欧美精品啪啪一区二区三区 | 婷婷成人精品国产| 老司机靠b影院| 女人被躁到高潮嗷嗷叫费观| 久久久久久久大尺度免费视频| 91字幕亚洲| 新久久久久国产一级毛片| 搡老岳熟女国产| 国产真人三级小视频在线观看| 日韩 欧美 亚洲 中文字幕| av欧美777| 黄片小视频在线播放| 美女高潮到喷水免费观看| 欧美激情高清一区二区三区| 久久免费观看电影| 久久性视频一级片| 国产国语露脸激情在线看| av片东京热男人的天堂| 婷婷丁香在线五月| 搡老乐熟女国产| 日日爽夜夜爽网站| 夜夜骑夜夜射夜夜干| 亚洲专区国产一区二区| 飞空精品影院首页| 精品亚洲成国产av| 男人添女人高潮全过程视频| 亚洲精品一二三| 久久精品人人爽人人爽视色| 免费女性裸体啪啪无遮挡网站| 麻豆av在线久日| 在线观看免费午夜福利视频| 最新在线观看一区二区三区| 黑人猛操日本美女一级片| 考比视频在线观看| 在线天堂中文资源库| 99精国产麻豆久久婷婷| 交换朋友夫妻互换小说| 桃红色精品国产亚洲av| 三级毛片av免费| 美女扒开内裤让男人捅视频| 91老司机精品| 夜夜骑夜夜射夜夜干| 久久人人爽人人片av| 久久久久久久国产电影| 精品一区二区三卡| 国产又色又爽无遮挡免| 成人免费观看视频高清| 亚洲黑人精品在线| 欧美日韩视频精品一区| 18禁裸乳无遮挡动漫免费视频| 久久午夜综合久久蜜桃| 精品国产国语对白av| 97精品久久久久久久久久精品| 18禁观看日本| 成年动漫av网址| 国产亚洲欧美在线一区二区| 亚洲av日韩在线播放| 后天国语完整版免费观看| 亚洲免费av在线视频| 久久影院123| 人妻人人澡人人爽人人| 建设人人有责人人尽责人人享有的| 精品人妻1区二区| 1024香蕉在线观看| 俄罗斯特黄特色一大片| 香蕉国产在线看| 91大片在线观看| 欧美日韩精品网址| 国产精品一区二区精品视频观看| 老熟妇仑乱视频hdxx| 日韩精品免费视频一区二区三区| 捣出白浆h1v1| 69av精品久久久久久 | 国产一区二区三区av在线| 在线av久久热| 日本av手机在线免费观看| 国产老妇伦熟女老妇高清| 亚洲国产欧美日韩在线播放| 国产精品香港三级国产av潘金莲| 国产一区二区三区综合在线观看| 久久精品国产亚洲av高清一级| 久久性视频一级片| 91麻豆av在线| 日日摸夜夜添夜夜添小说| 日韩 亚洲 欧美在线| 97在线人人人人妻| 在线精品无人区一区二区三| 国产成人精品久久二区二区91| 免费在线观看黄色视频的| 欧美另类一区| 丰满少妇做爰视频| 日韩中文字幕欧美一区二区| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲国产日韩一区二区| 久久免费观看电影| 亚洲av男天堂| 亚洲 欧美一区二区三区| 黄色视频在线播放观看不卡| 嫩草影视91久久| 亚洲av美国av| 精品福利永久在线观看| 少妇人妻久久综合中文| 青青草视频在线视频观看| a级毛片黄视频| 97在线人人人人妻| www.999成人在线观看| 久久免费观看电影| 18禁裸乳无遮挡动漫免费视频| 亚洲色图综合在线观看| 97精品久久久久久久久久精品| 国产亚洲av片在线观看秒播厂| 国产一区二区 视频在线| 欧美性长视频在线观看| 丝瓜视频免费看黄片| 亚洲成av片中文字幕在线观看| 啪啪无遮挡十八禁网站| e午夜精品久久久久久久| 亚洲精品av麻豆狂野| 国产精品国产三级国产专区5o| 精品国产一区二区三区四区第35| 国产免费现黄频在线看| 嫁个100分男人电影在线观看| 亚洲综合色网址| 一区二区三区乱码不卡18| 亚洲精品国产av蜜桃| 午夜福利,免费看| 美女大奶头黄色视频| 十分钟在线观看高清视频www| 啦啦啦在线免费观看视频4| 一区二区三区激情视频| 久久亚洲精品不卡| 久久人妻熟女aⅴ| 日本vs欧美在线观看视频| 中文字幕av电影在线播放| 另类精品久久| 宅男免费午夜| 久久天躁狠狠躁夜夜2o2o| 婷婷丁香在线五月| 一区二区三区乱码不卡18| 国产深夜福利视频在线观看| 在线天堂中文资源库| 久久九九热精品免费| 欧美精品av麻豆av| a 毛片基地| 日韩电影二区| 国产免费福利视频在线观看| 国产区一区二久久| 9191精品国产免费久久| 成年av动漫网址| 91成年电影在线观看| 99久久精品国产亚洲精品| 在线十欧美十亚洲十日本专区| 在线永久观看黄色视频| 欧美日韩福利视频一区二区| 中文字幕精品免费在线观看视频| av视频免费观看在线观看| 在线看a的网站| 可以免费在线观看a视频的电影网站| 丰满少妇做爰视频| 97在线人人人人妻| 亚洲人成77777在线视频| 午夜老司机福利片| 91精品国产国语对白视频| 12—13女人毛片做爰片一| 老汉色av国产亚洲站长工具| 亚洲美女黄色视频免费看| 亚洲情色 制服丝袜| 天堂俺去俺来也www色官网| 国产不卡av网站在线观看| 亚洲精华国产精华精| 欧美中文综合在线视频| 9热在线视频观看99| 一区福利在线观看| 亚洲av欧美aⅴ国产| 亚洲avbb在线观看| 欧美日韩av久久| 日韩欧美一区视频在线观看| 男女下面插进去视频免费观看| 大片电影免费在线观看免费| 女性被躁到高潮视频| 精品一区二区三卡| 欧美精品一区二区大全| 色综合欧美亚洲国产小说| 热99国产精品久久久久久7| 日本欧美视频一区| 大型av网站在线播放| 午夜久久久在线观看| 少妇 在线观看| www.999成人在线观看| 大码成人一级视频| 欧美精品啪啪一区二区三区 | 亚洲人成77777在线视频| 菩萨蛮人人尽说江南好唐韦庄| 精品卡一卡二卡四卡免费| 极品少妇高潮喷水抽搐| 少妇裸体淫交视频免费看高清 | 999精品在线视频| 日韩三级视频一区二区三区| 一进一出抽搐动态| 国产免费一区二区三区四区乱码| 精品人妻一区二区三区麻豆| 免费女性裸体啪啪无遮挡网站| 成人三级做爰电影| 下体分泌物呈黄色| 热99久久久久精品小说推荐| 欧美亚洲 丝袜 人妻 在线| 9色porny在线观看| 成人av一区二区三区在线看 | 看免费av毛片| 亚洲精品国产av蜜桃| 99精国产麻豆久久婷婷| 少妇的丰满在线观看| 国产精品久久久久久精品电影小说| 80岁老熟妇乱子伦牲交| 欧美亚洲 丝袜 人妻 在线| 高潮久久久久久久久久久不卡| 纯流量卡能插随身wifi吗| 亚洲av电影在线进入| av又黄又爽大尺度在线免费看| 黄色视频在线播放观看不卡| 亚洲精品一二三| 美女主播在线视频| 亚洲欧美一区二区三区久久| 巨乳人妻的诱惑在线观看| 色精品久久人妻99蜜桃| 人人澡人人妻人| 亚洲专区中文字幕在线| 久久人妻福利社区极品人妻图片| 777米奇影视久久| 一二三四在线观看免费中文在| 久久av网站| 精品国产乱码久久久久久小说| 国产成人a∨麻豆精品| 天天躁狠狠躁夜夜躁狠狠躁| 免费观看人在逋| 亚洲熟女精品中文字幕| 久久国产亚洲av麻豆专区| 一边摸一边做爽爽视频免费| av欧美777| 国产一区二区三区av在线| 精品卡一卡二卡四卡免费| 精品国产一区二区三区四区第35| 亚洲精品久久久久久婷婷小说| 国产99久久九九免费精品| 男女午夜视频在线观看| 美女大奶头黄色视频| 午夜免费鲁丝| 91成人精品电影| 1024视频免费在线观看| av国产精品久久久久影院| 国产亚洲精品第一综合不卡| 一级片免费观看大全| 精品国内亚洲2022精品成人 | 免费观看人在逋| 欧美精品av麻豆av| 男人舔女人的私密视频| 男女免费视频国产| 黑人欧美特级aaaaaa片| 午夜免费成人在线视频| av超薄肉色丝袜交足视频| 亚洲 欧美一区二区三区| 无遮挡黄片免费观看| 汤姆久久久久久久影院中文字幕| 日日夜夜操网爽| 国产又爽黄色视频| 欧美精品一区二区免费开放| 欧美日韩国产mv在线观看视频| 日韩电影二区| 69精品国产乱码久久久| 一本大道久久a久久精品| 宅男免费午夜| 免费不卡黄色视频| 秋霞在线观看毛片| 亚洲人成77777在线视频| svipshipincom国产片| 12—13女人毛片做爰片一| 成年美女黄网站色视频大全免费| 亚洲av电影在线观看一区二区三区| 国产免费一区二区三区四区乱码| 老熟女久久久| a级毛片在线看网站| 色综合欧美亚洲国产小说| 黄色毛片三级朝国网站| 动漫黄色视频在线观看| 欧美午夜高清在线| 久久免费观看电影| 美女扒开内裤让男人捅视频| 我要看黄色一级片免费的| 97在线人人人人妻| 视频区欧美日本亚洲| 桃红色精品国产亚洲av| 欧美另类一区| 日韩 亚洲 欧美在线| av欧美777| 人成视频在线观看免费观看| 99久久精品国产亚洲精品| 我的亚洲天堂| 韩国高清视频一区二区三区| 日韩欧美一区二区三区在线观看 | 国产又爽黄色视频| 亚洲九九香蕉| 亚洲国产av影院在线观看| 日日摸夜夜添夜夜添小说| 老熟妇仑乱视频hdxx| 肉色欧美久久久久久久蜜桃| 大码成人一级视频| 我要看黄色一级片免费的| 人人妻,人人澡人人爽秒播| 美女高潮到喷水免费观看| 国产一区二区三区av在线| 热99国产精品久久久久久7| 国产精品二区激情视频| 日韩中文字幕视频在线看片| 丝袜在线中文字幕| 欧美精品人与动牲交sv欧美| 国产成人免费无遮挡视频| 青春草视频在线免费观看| 99精品久久久久人妻精品| 日韩 亚洲 欧美在线| 国产男女超爽视频在线观看| 久久ye,这里只有精品| 亚洲成人免费av在线播放| 日本av免费视频播放| 曰老女人黄片| 90打野战视频偷拍视频| 丰满迷人的少妇在线观看| 午夜福利视频在线观看免费| 亚洲一卡2卡3卡4卡5卡精品中文| 久久久久网色| 国产一区二区三区在线臀色熟女 | 国产片内射在线| tube8黄色片| 一区在线观看完整版| 精品福利观看| 叶爱在线成人免费视频播放| 黄色怎么调成土黄色| av欧美777| 日本撒尿小便嘘嘘汇集6| 老汉色∧v一级毛片| av超薄肉色丝袜交足视频| 午夜福利影视在线免费观看| 亚洲精品美女久久久久99蜜臀| √禁漫天堂资源中文www| 新久久久久国产一级毛片| 69精品国产乱码久久久| 久热爱精品视频在线9| 香蕉丝袜av| 丝袜美足系列| 高清视频免费观看一区二区| 午夜成年电影在线免费观看| 国产野战对白在线观看| 69精品国产乱码久久久| 国产视频一区二区在线看| 男女免费视频国产| 91九色精品人成在线观看| 捣出白浆h1v1| 中文字幕最新亚洲高清| 久久狼人影院| 亚洲免费av在线视频| 成年av动漫网址| 两性夫妻黄色片| av免费在线观看网站| 精品免费久久久久久久清纯 | 国产亚洲欧美在线一区二区| 国产日韩一区二区三区精品不卡| 精品视频人人做人人爽| 国产区一区二久久| 成人国产一区最新在线观看| 久9热在线精品视频| 建设人人有责人人尽责人人享有的| 欧美精品人与动牲交sv欧美| 又黄又粗又硬又大视频| 91成年电影在线观看| 日日摸夜夜添夜夜添小说| 精品乱码久久久久久99久播| 国产激情久久老熟女| 少妇精品久久久久久久| 日韩制服骚丝袜av| 好男人电影高清在线观看| 成人黄色视频免费在线看| 国产精品一二三区在线看| 色精品久久人妻99蜜桃| 国产老妇伦熟女老妇高清| 免费观看人在逋| 欧美日韩国产mv在线观看视频| 女性生殖器流出的白浆| 精品少妇久久久久久888优播| kizo精华| 中文字幕最新亚洲高清| 国产主播在线观看一区二区| 国产男人的电影天堂91| 久久久国产欧美日韩av| 极品人妻少妇av视频| 老熟女久久久| 亚洲欧美成人综合另类久久久| 中文字幕制服av| 国产熟女午夜一区二区三区| 国产精品久久久久久精品古装| 日韩大码丰满熟妇| 久久人人爽人人片av| 亚洲激情五月婷婷啪啪| 狂野欧美激情性bbbbbb| 伊人亚洲综合成人网| av免费在线观看网站| 久久久久久久久久久久大奶| 老司机在亚洲福利影院| 天天操日日干夜夜撸| 亚洲一区中文字幕在线| 69av精品久久久久久 | 国产在视频线精品| 美女脱内裤让男人舔精品视频| 久久精品国产综合久久久| 久久国产精品大桥未久av| 99九九在线精品视频| 国产av精品麻豆| 中文字幕精品免费在线观看视频| 欧美日韩亚洲综合一区二区三区_| 亚洲国产欧美网| 美女脱内裤让男人舔精品视频| 人妻人人澡人人爽人人| 久久人妻福利社区极品人妻图片| 亚洲精品国产av成人精品| 少妇裸体淫交视频免费看高清 | av免费在线观看网站| 1024视频免费在线观看| 精品少妇一区二区三区视频日本电影| 日韩大码丰满熟妇| 国产精品一二三区在线看| 两性夫妻黄色片| 久久久久久人人人人人| 国产欧美日韩综合在线一区二区| 亚洲av成人不卡在线观看播放网 | 亚洲国产欧美日韩在线播放| 69精品国产乱码久久久| 伊人久久大香线蕉亚洲五| 老司机午夜十八禁免费视频| 无遮挡黄片免费观看|