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

    Predicting lean blow-off of bluffbody stabilized flames based on Damk?hler number

    2019-02-27 08:59:22ZhonghoWANGBinHUAimingDENGJunhuZHANGQingjunZHAO
    CHINESE JOURNAL OF AERONAUTICS 2019年2期

    Zhongho WANG,Bin HU,Aiming DENG,Junhu ZHANG,Qingjun ZHAO,c,*

    aInstitute of Engineering Thermophysics,Chinese Academy of Sciences,Beijing 100190,China

    bUniversity of Chinese Academy of Sciences,Beijing 100049,China

    cKey Laboratory of Light-duty Gas-turbine,Chinese Academy of Sciences,Beijing 100190,China

    Abstract Lean Blow-Off(LBO)prediction is important to propulsion system design.In this paper,a hybrid method combining numerical simulation and Da(Damk?hler)model is proposed based on bluffbody stabilized flames.In the simulated reacting flow,Practical Reaction Zone(PRZ)is built based on OH radical concentration,and it is considered to be the critical zone that controls LBO.Da number is obtained based on the volume-averaged parameters of PRZ.The flow time scale(τf)indicates the residence time of the fresh mixture flowing through the PRZ.It is obtained based on the characteristic length and volume-averaged axial velocity of the PRZ.The chemical time scale(τc)indicates the shortest time needed to trigger the reaction of the mixture.It is obtained by commercial software CHEMKIN through monitoring the transient variation of the reactor temperature.The result shows that the average Da number under different LBO conditions is 1.135(the Da number under each LBO condition ranges from 0.673 to 1.351).This indicates that the flow time scale and chemical time scale are comparable.The combustion is in a critical state where LBO is easy to occur.With the increase of the fuel mass flow rate(the global fuel/air ratio increases from 0.004761 to 0.01095),τfincreases from 0.001268 s to 0.007249 s,and τcdecreases from 0.00124 s to 0.00089 s.Accordingly,Da number increases from 1.023 to 8.145,which shows that the combustion becomes more stable.The above results show that the method proposed in the present study can properly predict the LBO limits of combustors,which provides important technical supports for combustor design and optimization.

    KEYWORDS Bluffbody stabilized flames;Damk?hler number;Hybrid method;LBO;Time scales

    1.Introduction

    Bluffbody stabilization of flames is commonly used in aeroengine combustors,especially in augmentors.At inlets of augmentors,the oxygen concentration and pressure are relatively low(generally,the total pressure of augmentor inlet is 0.10-0.40 MPa,and the oxygen mass fraction is 12-17%),which goes against flame stabilization.Therefore,Lean Blow-Off(LBO)performance of combustors needs to be estimated during design stage,and then a high-accuracy LBO prediction method is needed to provide enough technical supports for combustor optimization.

    Early LBO prediction methods are mainly semi-empirical models which were developed from experiments and theoretical derivation.These models could be classified into two categories:PSR (Perfectly Stirred Reactor)model1and CT(Characteristic Time)model.2These two models were both proposed in the study of bluffbody stabilized flames.Longwell made a homogeneity hypothesis to the recirculation zone and proposed the conception of PSR due to the drastic mixing process within it.1,3Flows of air and fuel enter the reactor and are instantaneously and perfectly mixed with the burning gases within the reaction volume.The gases in reactor are uniform in composition and temperature.LBO occurs when the heat supplied by the reactor is insufficient to raise the fresh mixture up to its ignition temperature.3Following Longwell et al.'s study,LBO prediction models were improved by Lefebvre and Rizk,and the models contain more factors affecting LBO,such as inlet flow velocity,reactor volume,pressure,temperature4,scale and shape of bluffbody5,6,fuel type7,etc.CT model explains LBO in a different angle of view.It was firstly proposed by Zukoski and Marbe.2They deemed that main reaction does not occur in recirculation zone.Hot products in recirculation zone only play the role of ignition source.Fresh mixture contacts with hot products at the edge of recirculation zone.If the fresh mixture is ignited by hot products within the residence time, flames can be stabilized.Otherwise,LBO occurs.Mellor et al.developed time models8,9and regarded the shear layer between the free stream and hot recirculation zone as the key region which is important to flame stabilization process.On the basis of this,they proposed that if the shear layer residence time of fuel is smaller than the sum of the fuel evaporation time and ignition delay time,LBO will occur.9

    With the development of computational technique,high spatial/temporal accuracy numerical accuracy simulation has provided researchers with powerful tools in flames dynamic analysis close to LBO.In the study of flame dynamics near LBO,LES-CMC (Large Eddy Simulation-Conditional Moment Closure)model was employed by Tyliszczak et al.10LES-PDF(Probability Density Function)model by Black and Smith11and Gokulakrishnan et al.12LES-EBU(Eddy Break-Up)model by Kim et al.13LES-EDC(Eddy Dissipation Concept)model by Gokulakrishnan et al.14LES-LEM(Linear Eddy Mixing)model by Menon et al.15and LES-LTEM(Lagrangian Transport Element Method)model16by Erickson et al.When the combustion was close to LBO,similar flame behavior was detected in these simulation works.Local extinctions were observed in the above simulation works when flames were close to LBO.Lieuwen et al.summarized the relevant works and divided flame behavior into three steps based on their characteristics17:A)local extinction along the flame sheet;B)large scale wake disruption;C)the final blowout associated with wake cooling and shrinking.

    LES is well-suited for simulating flame behavior before LBO,and many details of the flow field can be reproduced.There are also some cases where LES is used for LBO predictions.Near LBO flame behavior obtained in the simulations is consistent with the experiment.However,when quantified LBO limits need to be obtained,relevant models should be improved further.LES-CMC was used for LBO prediction of two-phase flames by Tyliszczak et al.10It was found that by adjusting a model constant in the subgrid model,a reasonable LBO limit could be obtained.LES-CMC model was also employed by Zhang and Mastorakos in LBO prediction of the Sandia D-F flames.18The predicted LBO axial bulk velocity is about 60%higher than the experimental value.LES-PDF method was used by Gokulakrishnan et al.for the prediction of LBO limits of bluffbody stabilized propane flames.12Ultimate quantified LBO limits were not given in this work,but it can be known from the simulation that the LBO equivalence ratio is below 0.60.LES-PDF method was employed by Black and Smith for LBO prediction of an aero low emission fuel injector.11The fluid unsteadiness near the LBO condition was captured,but the localized fuel-rich pockets affected the accuracy of LBO prediction.There are few LES works predicting LBO under various operating conditions and combustor configurations.Even though accurate LBO prediction can be achieved under a specific condition,whether it could be used for other conditions is still questionable.It is still hard for LES to satisfy engineering needs in LBO prediction.

    At present,the primary problem impeding the application of high spatial/temporal numerical simulation in practical engineering is the huge computational cost,especially for combustion simulation,which couples chemical reaction in turbulent flow.According to Luo et al.'s work,typically,the thickness of reaction zone ofn-heptane flame is on the order of 2×10-5m,which restricts the grid size.Enormous calculation and memory resources are needed to resolve completely all flow structures.19As Menon have pointed out,LES can be very expensive and even though it can resolve the large scale structures accurately,the dynamics at the small scales still require modeling.It must be ensured that the high cost and uncertainties in the sub-grid models are not counter productive.20

    The combination of semi-empirical LBO prediction models and numerical simulation compromises the universality and computational cost,and it has been explored by some researchers.Bin Hu et al.improved Lefebvre's LBO prediction model by introducing Flame Volume(FV)concept to the model.FV is related with flow fields which are obtained by numerical simulation.FV model21,FIA (Fuel Iterative Approximation)22model,and CCFV(Characteristic Curve of Flame Volume)23model were developed successively on the basis of FV concept and improved Lefebvre's model.The prediction uncertainties are reduced from±30%to lower than±10%in average.Ahmed also adopted the flame volume conception and developed a LBO prediction method.24Zheng et al.deemed that CRZ(Central Recirculation Zone)is the critical zone controlling LBO.25The conception of Feature-Section was proposed and defined as the surfaces where axial velocity equals zero.Feature-parameter Γ was calculated based on the fuel/air ratio and the changing rate of average temperature of Feature-Sections.The inflection point of Γ-fuel/air ratio curve indicates the occurrence of LBO.Kiel and Roach26,and Knaus27,28et al.predicted LBO of bluffbody stabilized flames by combiningDamodel with numerical simulation.In their approaches,different definitions of time scales contained inDamodel were proposed.From the aspect of turbulent mixing, flow time scale was defined as the ratio of turbulent kinetic energy(k)and dissipation rate(ε).26,27The chemical time scales were defined based on fuel consumption rate of each cell26and the time needed for OH radical mole fraction reaches maxima27.From the aspect of local extinction induced by high strain rate,the flow time scale was defined as the inverse of local strain rate.The chemical time scale was defined as the inverse of extinction strain rate obtained by CHEMKIN.28When the lobes whereDa<1 merge together downstream of the bluffbody,LBO is considered to occur.

    In this paper,a hybrid LBO prediction method is proposed by combining numerical simulation andDamodel based on the Practical Reaction Zone of bluffbody stabilized flames.This method is a reasonable compromise between semiempirical models and numerical simulation.LBO limits can be accurately predicted with lower computational cost,which may be a feasible approach for LBO limits prediction in practical engineering.

    2.General idea and theory

    The LBO prediction method proposed in this study combinesDamodel and numerical simulations.PRZ is characterized from simulated reacting flow fields based on OH radical concentration,and it is considered as the critical zone controlling LBO.Danumber is defined as the ratio of flow time scale and chemical time scale,expressed asDa=τflow/τchem.Factors affecting LBO can be embodied concisely by the two time scales inDamodel.All parameters related withDanumber are obtained based on volume-averaged parameters of PRZ.

    2.1.Practical reaction zone(PRZ)

    The schematic diagram of PRZ in combustor is shown in Fig.1.The pink area indicates PRZ.Intense reactions that may affect LBO all occur there.Recirculation brings hot products upstream,but only weak reactions really occur in recirculation zone.Flames are anchored on the trailing edges of the V-gutter,and develop downstream.

    Fig.1 Schematic diagram of PRZ in bluffbody stabilized flames.

    PRZ is considered as the critical zone controlling LBO in the present study.There are many definitions of PRZ.Nonpremixed combustion stabilized by bluffbody contains many factors,such as turbulent flow,chemical reaction,heat and mass transfer,fuel injection and atomization.The definition of critical zones should take all the factors into consideration.At present,the definitions of critical zone can be divided into three categories:A)Considering the recirculation zone as PRZ.It is true that recirculation zone is important to flame stabilization.However,in the large part of recirculation zone,reactant has already been consumed and only weak reaction exists.B)PRZ defined based on fuel concentration.The main problem of this approach is that once the reaction is started,the fuel is consumed and initial fuel concentration field is destroyed.Some researchers obtain the fuel concentration field by simulating non-reacting flow,but the effects of fuel evaporation and gas expansion are not taken into consideration.C)PRZ defined by temperature.It is true that the primary characteristic of combustion is high temperature,but high temperature may also exist in products where the reaction has already been ended.

    According to the above discussion,in the present study,PRZ is defined based on OH radicals.The region where OH radical concentration is larger than a certain value is considered to be PRZ.The generation and consumption of OH radicals are always accompanied with intense heat absorbing and releasing.If OH radicals are not produced in a chemical reaction,the reaction cannot provide enough radicals and heat needed for stable combustion.29The existence of OH radicals shows whether effective reactions exist,and the concentration of OH radicals indicates the intensity of the reactions.The specific value of OH radical concentration on the PRZ boundary is determined by comparing the PRZ obtained in simulated flow field and the experiment.30

    2.2.Time scales

    There are various definitions of the time scales contained inDamodel.Systematic and comprehensive summary of the definition ofDanumber can be referred in Husain31and Cavaliere et al.'s work.32Shanbhogue et al.have pointed out thatDamodels are suitable for the prediction of local extinction,but do not fundamentally describe the ultimate LBO.17The innovation point of the current work is making it feasible to predict LBO byDamodels through the definition of PRZ.According to Shanbhogue et al.'s work,17different definitions ofDanumber describe flame behaviors from different angles.The definition ofDanumber can be divided mainly into 3 categories.The first definition is based on turbulent flow.The combustion is assumed to propagate from one Kolmogorov-scale vortex to a neighbouring one.32The flow time scale is defined as the time needed for flame propagation,and the chemical time scale is defined as the lifetime of the large eddies of the turbulence.Local extinction occurs when flow time scale is larger than chemical time scale.The second definition ofDanumber is based on strain rate.The flow time scale is typically defined as the reciprocal of the local strain rate,and the chemical time scale is defined as the reciprocal of the extinction strain rate.28Local extinction occurs whenDanumber is smaller than 1.The third definition ofDanumber is based on reaction zone.The flow time scale is typically defined as the residence time of the fresh mixture in the reaction zone,and the chemical time is generally used to indicate the LBO critical flow time.The first two definitions ofDanumber are suitable for predicting local extinction.However,when they are applied to LBO prediction,a proper criterion is lacked,which describes to what extent the extinctions may lead to LBO.The third definition ofDanumber is widely used in LBO prediction and has the potential to obtain accurate quantitative LBO limits.The pre-diction accuracy relies on the proper definition of the reaction zone and time scales.It will be very efficient to obtain the information needed in the semi-empirical models by numerical simulation.The target for the current study is the judgement of LBO,instead of the reproduction of the transient process of LBO,and thereby,RANS can satisfy the needs in this work.The detailed introductions of the two time scales will be presented below.

    2.2.1.Flow time scale(τf)

    The flow time scale in the present study is defined as the ratio of characteristic length scale(L)and characteristic velocity(v).Lis defined as the ratio of PRZ volume(VPRZ)and a characteristic area(A),expressed asL=VPRZ/A.In the present work,the characteristic area is suggested to be the product of section area(Ac)of the combustor and blockage ratio(BR),expressed asA=Ac×BR.vis the volume-averaged axial velocity of flows in PRZ.Flow time scale indicates the residence time of fresh mixture flowing through the PRZ.

    2.2.2.Chemical time scale(τc)

    Chemical time scale is defined as the shortest time needed to trigger the reaction in PRZ.It is obtained by the PSR module in CHEMKIN.

    In this paper,a reduced 16 species-23 steps chemical kinetic mechanism33is introduced into CHEMKIN.Temperatures,pressures and equivalence ratios needed in CHEMKIN are obtained by calculating the volume-averaged parameters of PRZ in simulated flow field.In CHEMKIN,the initial residence time is set long enough to establish a stable combustion,and then it is reduced slowly until a sudden drop of temperature appears,which indicates that the reaction cannot be triggered.This residence time is considered to be the minimum time needed to trigger the reaction and written as τc.

    An example of τcacquisition in CHEMKIN is shown in Fig.2 (the operating condition of the combustorisMainlet=0.184,Tinlet=546.15 K,FAR=0.00399).It is shown that there is pulse-rise of reactor temperature at the beginning of the ignition process both in Fig.2(a)and(b).This is caused by the initial high temperature(1706 K)set in PSR reactor.The initial temperature is assigned as the volumeaveraged temperature of PRZ since the fresh mixture entering the PRZ is ignited by hot gas when a stable combustion is established.The pulse-rise of reactor temperature does not affect the final calculation of chemical time scale.

    In Fig.2(a),after the reaction has been triggered,the reactor temperature drops to about 1940 K and maintains stable.The result indicates that the residence time of 0.00103 s is enough for reaction to be triggered.In Fig.2(b),when the residence time is reduced to 0.00102 s,the reactor temperature drops rapidly to the inlet temperature(546.15 K).The residence time of 0.00102 s is too short for the reaction to be triggered,and then the chemical time scale under this condition is obtained as 0.00103 s.

    3.Test rig

    The LBO experiment was conducted by Wang and Jin.30For the clarity of the work presented here,the test rig is concisely introduced.The schematic of test rig is shown in Fig.3.The combustor is a 2D rectangular combustor.Incoming air is heated by 3 BK-I heaters.A bluffbody is installed 240 mm downstream of the combustor inlet.A fuel manifold is installed 13 mm upstream of the bluffbody.Fresh mixture is ignited by high-energy ignition plug.The post-test section contains temperature and pressure probes and an exhaust nozzle.The distance between the high energy ignition plug and bluffbody is 80 mm.The outlet temperature is measured by Pt-Rh_30/Pt-Rh_6 thermocouples,and the mass flow rate of the fuel is measured by DMF-1-2 mass flow meter.

    There are 3 fuel injecting holes located symmetrically with spanwise length of 55 mm.The diameter of each injecting hole is 0.6 mm,and the fuel is injected in the reverse direction of the air flow.The fuel adapted in the experiment is Chinese RP-3 aviation kerosene.At about 960 mm downstream of the bluffbody,there are 3 temperature measuring points distributed uniformly along the vertical direction with an interval of 30 mm.The static pressure of incoming air flow is maintained 1 atm,the inlet temperature is 490-620 K,and theManumber of the incoming air is 0.17-0.25.For more details,please refer to Ref.30.

    The LBO experimental data are listed in Table 1.In the present work,the LBO data are selected and categorized based on the sameMainletorTinletto analyze the effects of single parameters on LBO.The LBO data categorizations are shown in Fig.4.In group 1 and group 2,Tinletis almost constant,and the maximum difference ofTinletis about 8 K.The average inlet temperature(Tinlet,average)of group 1 and group 2 are 495.9 K and 547.15 K,respectively.In group 3 and group 4,Mainletis almost constant,and the maximum difference ofMainletis about0.011.The average in letManumber(Mainlet,average)of group 3 and group 4 are 0.217 and 0.205,respectively.

    Fig.2 Temperature-time curve obtained by PSR module in CHEMKIN.(Mainlet=0.184,Tinlet=546.15 K,FAR=0.00399).

    Fig.3 Test rig.

    Table 1 Experimental LBO data.

    4.Numerical simulation

    4.1.Turbulence model

    Realizablek-ε model is adapted to simulate turbulent flows.34The model satisfies certain mathematical constraints on the Reynolds stresses to be consistent with the physics of turbulent flows,and could more accurately predict the spreading rate of both planar and round jets.It is also likely to provide superior performance for flows involving rotation and recirculation.

    Realizablek-ε model is different from RNGk-ε model and standardk-ε model mainly in turbulent viscosity(μt)calculation and transport equation for the dissipation rate(ε).Inkε models,μt=ρCμk2/ε.In RNGk-ε model and standardk-ε model,Cμis constant,but in realizablek-ε,Cμis a function of the mean strain and rotation rates,the angular velocity of the system rotation,and the turbulence fields(kand ε),as shown in Eq.(1)-(4).

    whereA0andAsare empirical constants related with strain rate,and

    and

    The modeled transport equations forkand ε in the realizablek-ε model are Eqs.(5)and(6).

    and

    Fig.4 LBO data grouping.

    where

    ρ denotes density,Gkrepresents the generation of turbulence kinetic energy due to the mean velocity gradients,Gbis the generation of turbulence kinetic energy due to buoyancy,YMdenotes the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate,σkand σεare the turbulent Prandtl numbers forkand ε,C1εandC3εare constants,SkandSεare user-defined source terms.The production term in the ε equation(the second term on the right hand of Eq.(6)does not involve the production ofk.It is believed that the present form better represents the spectral energy transfer.The destruction term(the next to last term on the right-hand side of Eq.(6)does not have any singularity.This is contrasted with traditionalk-ε models.

    4.2.Combustion model

    Laminar flamelet model is employed in the simulation of combustion.35The flamelet model views the turbulent flame as an ensemble of thin,laminar,locally one-dimensional flamelet structures within the flow field.The governing equations can be simplified to one dimension along the axis of the fuel and oxidizer jets,where complex chemistry calculations can be affordably performed.Mixture fraction theory is applied to non-premixed combustion simulation.A PDF table is generated based on flamelet model.36Temperature and species fraction are functions of mean mixture fraction(f)and mixture fraction variance.fcan be written in term of the atomic mass fraction as

    whereZiis the elemental mass fraction for elementi.Sincefis a conserved variable,local equivalence ratio can be conveniently calculated without considering the conversion of fuel to products.

    Turbulence-chemistry interaction is modelled by assumed PDF model.As only one type of fuel is included in the present study,β-shaped PDF closure model is adopted.

    According to Claramunt et al.'s work,the prediction accuracy of the majority of radicals through flamelet model is acceptable.37Ravikanti et al.employed laminar flamelet model in simulating the flow field of a bluffbody combustor,and the prediction of OH radicals fitted well with the experiment results.38The laminar flamelet model also performed well in Wang et al.'s work in separating the reaction zone based on OH radicals in a primary combustor.The reaction zone obtained in simulation fitted well with the photo taken during the experiment.39

    Kundu et al.'s 16 species-23 steps kinetic mechanism is employed in combustion simulation.33The chemical mechanism has been validated by the authors and other researchers.The chemical mechanism was employed by Zhang et al.in the study of radical recombination in a hydrocarbon-fueled scramjet nozzle.40The mechanism was validated by comparing the calculated ignition delay time with the experimental data.The mechanism was also employed by Siow and Yang in 3D analysis of Jet-A combustion.41The calculated CO,CO2and NO emissions fitted well with the experimental data.The mechanism has also been validated by Wang et al.39,42in the calculation of the axial temperature distribution of a swirl stabilized combustor.The calculated value fitted well with the measured one.

    4.3.Discrete model

    Euler-Lagrange approach method43is employed in the calculation of discrete phase.The dispersed phase is presumed to occupy a small volume fraction,and thus the effects of droplet collisions and break-up are neglected.The dispersed phase can exchange momentum,mass,and energy with the fluid phase.The droplet trajectories are computed individually at specified intervals during the fluid phase calculation.Spherical drag law44is employed for drag coefficient.The injection type is assumed as slim solid cone,with 5°cone angle.The droplet size distribution is assumed uniform with an average diameter of 60 μm.

    4.4.Model combustor and boundary conditions

    The computational domain of the combustor is shown in Fig.5.The size of the whole combustor and locations of the fuel injecting holes in numerical simulation are the same with those in the experiment.

    Second-order upwind spatial discretion is employed in all cases.Standard wall functions are adopted.The hydraulic diameter of the inlet and outlet of the model combustor is 159.375 mm.Kundu et al.'s 16 species-23 steps kinetic mechanism is employed in combustion simulation.33The inlet conditions in numerical simulation are the same with those in experiment.30

    Fig.5 Model combustor and boundary conditions.

    4.5.Grid generation

    The grids of the flow field are generated by commercial software GAMBIT.As shown in Fig.6,structured grids are generated all over the flow field.The total number of the grids is about 2.65 M.The work of grid dependence check has been done by Wang and Jin in their works.30The grid dependence check was performed with 1.5 M and 3.0 M cells.The results show that 1.5 M cells were enough for the computations.In the current study,2.65 M cells are used to obtain smoother data distributions.

    4.6.Validation of numerical methods

    The numerical methods employed in this work are validated based on experimental temperature data measured downstream of the bluffbody.

    Total temperature is measured under two operating conditions atx=1200 mm(Case 1:Tinlet=537 K,Mainlet=0.16,FAR=0.0048;Case2:Tinlet=523 K,Mainlet=0.155,FAR=0.010).Distributions of simulated and experimental temperature alongydirection atx=1200 mm under the two operating conditions are shown in Fig.7.

    The highest temperature appears at the centerline of the combustor,and the temperature profile alongydirection is similar with that in experiments.The error between simulation and experiment is 7.34%in average.In the present study,the accuracy of simulation is acceptable.During the experiment,the liquid fuel accumulated on the upwind surfaces of the bluffbody and formed a liquid film.When the accumulated fuel was blown downstream,the fuel film would undergo the atomization process again.This process cannot be reproduced in simulation,and thus deviation between calculated and measured values is inevitable.

    Fig.6 Grids of computational zone.

    5.Results and discussion

    5.1.Effects of FAR on volume of PRZ

    In the present study,the specific value of OH radical concentration on PRZ boundary is determined by comparing the PRZ of simulated reacting flow field and the flame pictures captured by high-speed camera in the experiment30.Currently,it is still difficult to quantitatively compare the OH concentration obtained by PLIF and numerical simulation in the entire flow field.Qualitative comparison is employed in the work.Khalil et al.pointed out that the distribution of OH radicals could be used to characterize PRZ,but the OH radicals occupy a large portion of the combustor for non-premixed mode,and it is not reasonable to define all the zone where OH exists as PRZ.45St?hr et al.also pointed out that the medium and low levels of OH may represent burned gas whose OH concentration has decayed while it was transported away from the PRZ.29In Cavaliere et al.'s work,32the behavior of PRZ was also reflected by high-OH zone rather than all the zone where OH radicals exist.Hence,a specific value of OH concentration on PRZ boundary should be determined in the current study.For the present work,the boundary value is obtained by comparing the OH concentration field with the photos of the PRZ,30and the value of 0.01 mole/m3is considered to be proper to describe the flame behavior.Fig.8 shows PRZ obtained by simulation and experiments.By comparing Fig.8(a)and(b),it can be concluded that the size and shape of PRZ in simulated flow field fit well with the photograph(inlet condition:Tinlet=537 K,Mainlet=0.16,FAR=0.0048).As the fuel is injected in the reverse direction of the incoming air flow,the fuel/air mixing is extremely nonuniform.Much fuel is gathered at the surface of the bluffbody.Under the effect of aerodynamic force,fuel drops migrate to the trailing edges of the bluffbody,and then flow into the shear layer downstream.

    Fig.7 Simulated and experimental temperature along y direction at x=1200 mm.

    Fig.8 PRZ obtained by simulation and experiment.

    Fig.9 3D description of PRZ under different FAR.

    Fig.10 Variation of VPRZ,LBOwith FARLBO.

    Fig.9 shows the 3D description of PRZ under different FAR.The axial length of PRZ increases rapidly with the increase of PRZ.Because PRZ is limited by combustor in radial direction,with the increase of fuel supply,the cross section area of PRZ changes little,and the increase of volume of Practical Reaction Zone(VPRZ)is mainly embodied in the increase of its axial length.As shown in Fig.9(e),when FAR is high,PRZ even extends out of the combustor.

    Fig.10 shows the variation ofVPRZ,LBOwith FARLBO.Fig.10 shows thatVPRZ,LBOincreases with the increase of FARLBO.Because the fuel mass flow rates set in numerical simulation are obtained from experiments,larger FARLBOwill cause larger fuel flow rate at LBO,and then enlargeVPRZ,LBO.The data points ofVPRZare fitted by a straight line.It should be noted in Fig.10 that with the continuative decrease of FAR along the fitting line,a minimum FARLBOcan be obtained as 0.0037 whenVPRZ,LBOinfinitely approximates 0.This value may be the minimum FARLBOof the current combustor style.Similar results are also obtained as FARLBO=0.002523and FARLBO=0.00342by Bin Hu et al.,which is a little lower than that in the paper due to the better fuel/air mixing in swirling flows.Owing to the lack of sufficient experimental data,the minimum FARLBOis only a conjecture and cannot be verified in the present work.It may be used as an estimated value for engineers in combustor designs.

    5.2.Effects of Mainleton time scales of PRZ

    Effect ofMainleton the Volume of PRZ(VPRZ,LBO)under different LBO conditions is shown in Fig.11.Fig.11 shows thatVPRZ,LBOincreases with the increase ofMainlet.On one hand,fuel supply increases due to the deteriorated LBO performance caused by the increase ofMainlet.On the other hand,whenMainletincreases,the turbulence exchange rate in the shear layer downstream of the bluffbody increases,and then the distribution of fuel is more uniform,which leads to more complete combustion.Hence,VPRZ,LBOincreases with the increase ofMainlet..

    Effect ofMainleton the volume-averaged equivalence ratio of PRZ under different LBO conditions(ΦPRZ,LBO)is shown in Fig.12.In the paper,ΦPRZ,LBOis obtained by Eq.(11),whereris the air/fuel ratio on a mass basis at stoichiometric ratio.

    Fig.11 Effect of Mainleton VPRZ,LBO.

    Fig.12 Effects of Mainleton ΦPRZ,LBO.

    Fig.12 shows that ΦPRZ,LBOdecreases with the increase ofMainlet..Generally speaking,the global fuel/air ratio of LBO will increase with the increase ofMainlet,and then ΦPRZ,LBOshould also increase since PRZ is a local region of entire flow field.The results shown in Fig.12 go against the authors' conjecture,which may be attributed to the non-uniform fuel distribution.Fig.13 shows the mixture fraction distribution within PRZ under differentMainlet.In Fig.13,the white area closed to the upwind surface of bluffbody is the extreme-fuel-rich region where the reaction cannot be triggered,and it is caused by the gather of fuel at the upwind surface of bluffbody.Fig.13 shows that with the increase of FAR,the white area grows and extends deeper into PRZ,which causes that the high-mixture fraction-zone in PRZ decreases.In addition,the low-mixture fraction-zone increases gradually with the increase ofMainlet.Hence,ΦPRZ,LBOdecreases with the increase ofMainlet.

    Effect ofMainletonTPRZ,LBOis shown in Fig.14.Fig.14 shows thatTPRZ,LBOincreases with the increase ofMainlet.It is found that the variation of temperature shown in Fig.14 is not consistent with the variation of equivalence ratio in Fig.12.According to Fig.12,TPRZ,LBOshould increase for the increase of fuel concentration in PRZ.This may be attributed to the enhanced turbulence mixing in PRZ.The detailed information can be illustrated in Fig.15.Fig.15 shows the variation of high temperature zone and PRZ with the increase ofMainlet.The specific value of 1800 K is assigned to qualitatively characterize the high temperature zone in PRZ.PRZ is colored green,and high temperature zone is colored red.The fraction of PRZ overlapped by high temperature zone under each condition is given below the corresponding picture.From Fig.15,it can be seen that the percentage of PRZ overlapped by high temperature zone increases(from 9.18%to 23.95%)with the increase ofMainlet(ranging from 0.189 to 0.230).Therefore,a conclusion can be drawn that the increase ofMainletactually enhances the stirring strength in PRZ,and leads to the increase ofTPRZ,LBO.

    Fig.13 Mixture fraction contour within PRZ under different Mainlet.

    Fig.14 Effects of Mainleton TPRZ,LBO.

    Fig.15 3D description of PRZ and high temperature zone under different Mainlet.

    Effects ofMainleton volume-averaged pressure of PRZ under different LBO conditions(PPRZ,LBO)are shown in Fig.16.Corresponding total pressure loss coefficient(σ)is marked near each data point.Fig.16 shows thatPPRZ,LBOslightly decreases with the increase ofMainlet,and σ increases for the increase ofMainlet.With constant inlet static pressure(1 atm),the increase of total pressure loss will surely lead to the decrease of static pressure downstream of the bluffbody,and then lead to the decrease ofPPRZ,LBO.It can also be seen in Fig.16 thatPPRZ,LBOdoes not change much.In group 1,whenMainletincreases from 0.189 to 0.230,PPRZ,LBOonly decreases by 1.90%.In group 2,whenMainletincreases from 0.184 to 0.228,PPRZ,LBOonly decreases by 1.85%.

    Effects ofMainleton flow time scale under different LBO conditions(τf,LBO)are shown in Fig.17.Fig.17 shows that τf,LBOincreases with the increase ofMainlet.From Fig.9 it can be observed that the variation ofVPRZ,LBOis mainly embodied in axial direction.With the increase ofMainlet,the axial length of PRZ increases.In addition,Mainlethas little effect on volume-averaged axial velocity of PRZ due to the partial recirculation in PRZ.Therefore,τf,LBOincreases with the increase ofMainlet.

    Effects ofMainleton chemical time scale under different LBO conditions(τc,LBO)are shown in Fig.18.Fig.18 shows that τc,LBOincreases with the increase ofMainlet.It is shown in Fig.12 and Fig.14 that ΦPRZ,LBOdecreases with the increase ofMainletandTPRZ,LBOincreases with the increase ofMainlet.τc,LBOis mainly affected by these two parameters which have different effects on τc,LBO.It can be seen in Fig.14 that the variation ofTPRZ,LBOis less than 60 K,andTPRZ,LBOis higher than 1670 K.The maximum relative increment ofTPRZ,LBOin group 1 and group 2 is 2.07%.Considering that the variation ofTPRZ,LBOis relatively small,τc,LBOis mainly affected by ΦPRZ,LBOin the paper.Hence τc,LBOincreases with the increase ofMainlet.

    5.3.Effects of Tinleton time scales of PRZ

    Effects ofTinletonVPRZ,LBOare shown in Fig.19.Fig.19 shows thatVPRZ,LBOdecreases with the increase ofTinlet.LBO limit is extended for the increase ofTinlet.Therefore,fuel concentration in the combustor decreases due to the reduction ofmf(fuel mass flow rate).Hence,VPRZ,LBOdecreases with the increase ofTinletfor the shrinkage of PRZ.

    Effects ofTinleton ΦPRZ,LBOare shown in Fig.20.ΦPRZ,LBOincreases with the increase ofTinlet.In general,the increase ofTinletextends the LBO limit,and then the global equivalence ratio in the combustor should decrease.However,it is shown in Fig.20 that the effects ofTinleton ΦPRZ,LBOincrease with the increase ofTinlet.The reason for this can be illustrated clearly in Fig.21.With the increase ofTinlet,the fuel/air ratio gradually decreases as expected,and then the low-mixture fraction-zone decreases.Meanwhile,the extreme-fuel richarea close to the surface of the bluffbody decreases because of the better fuel evaporation and more complete combustion.Hence,ΦPRZ,LBOincreases with the increase ofTinlet.

    Fig.16 Effects of Mainleton PPRZ,LBO.

    Fig.17 Effects of Mainleton τf,LBO.

    Fig.18 Effects of Mainleton τc,LBO.

    Effects ofTinletonTPRZ,LBOare shown in Fig.22.TPRZ,LBOincreases with the increase ofTinlet.The increase ofTinletis good for complete evaporation and combustion of fuel.The reaction process can also be promoted by the increase ofTinlet.Compared with the volume high temperature zone,with the increase ofTinlet,the volume of Practical Reaction Zone decreases faster since it is directly affected by the fuel concentration in the combustor(Fig.23).The high temperature zone is characterized by the specific temperature 1800 K,and it can be seen that almost all the high temperature zone is coated by PRZ.Volume fraction of high temperature zone in PRZ increases from 21.33%to 21.53%.If the boundary of the high temperature zone is adjusted to 1600 K,volume fraction of high temperature zone in PRZ increases from 67.33%to 70.66%,the increasing trend becomes more obvious,and henceTPRZ,LBOincreases with the increase ofTinlet.

    Fig.19 Effects of Tinleton VPRZ,LBO.

    Fig.20 Effects of Tinleton ΦPRZ,LBO.

    Fig.21 Mixture fraction contour within PRZ under different Tinlet.

    Effects ofTinletonPPRZ,LBOare shown in Fig.24.Similar to Fig.16,the value of σ near the data points is the total pressure loss coefficient.It is shown in Fig.24 thatPPRZ,LBOincreases with the increase ofTinlet,and σ decreases with the increase ofTinlet.The global temperature rise in the combustor decreases for the decrease of fuel supply,and thus the pressure loss caused by combustion decreases,which leads to the reduction of the total pressure loss σ.With constant inlet static pressure(1 atm),PPRZ,LBOincreases slightly with the increase ofTinlet.It should be noted thatPPRZ,LBOis not changed much.In group 3,whenTinletincreases from 492.15 K to 611.15 K,PPRZ,LBOonly increases by 0.58%.In group 4,whenTinletincreases from 500.15 K to 560.15 K,PPRZ,LBOonly increases by 0.39%.

    Effects ofTinleton τf,LBOare shown in Fig.25.Fig.25 shows that τf,LBOdecreases with the increase ofTinlet.According to Fig.19,the increase ofTinletwill decrease the volume of PRZ.The effects of the gas expansion caused by the increase ofTinleton the volume-averaged axial velocity(v)are little due to the partial recirculation contained in PRZ.Therefore,τf,LBOdecreases with the increase ofTinlet.

    Effects ofTinleton τc,LBOare shown in Fig.26.Fig.26 shows that τc,LBOdecreases with the increase ofTinlet.According to Fig.20 and Fig.22,ΦPRZ,LBOandTPRZ,LBOboth increase with the increase ofTinlet.The reaction is thus easier to be triggered,and τc,LBOis reduced.

    Fig.22 Effects of Tinleton TPRZ,LBO.

    Fig.23 3D description of PRZ and high temperature zone under different Tinlet.

    5.4.Da numbers comparison

    Danumbers under all 16 LBO conditions are shown in Fig.27.The average value of theDanumbers is 1.135,the maximum is 1.351 and the minimum is 0.673.The results show that the flow time scale and the chemical time scale are comparable,indicating a critical state of LBO.In this state, flames are extremely unstable.If there is an abnormal perturbation in the incoming flow upstream,LBO is very likely to happen.The reason for the fluctuation ofDanumbers may be attributed to the difference between heterogeneous scalar distribution in PRZ and homogeneous hypothesis in the PSR module of CHEMKIN.

    The variation ofDanumber with different FAR underTinlet=497.15 andMainlet=0.18869 is shown in Fig.28.Fig.28 shows thatDanumber increases in an exponential way with the increase of FAR.The increase of FAR directly increases the fuel concentration in the combustor,and PRZ is expanded accordingly.Hence,τfincreases with the increase of FAR.In addition,the increase of fuel concentration also leads to an intensive reaction,which reduces the chemical time substantially.For the reasons discussed above,Danumber increases rapidly with the increase of FAR.

    The work presented here provides combustor designers with a practical tool that links theDamodel and simulated reacting flow together.Relevant work will be continued and focus more on the following aspects:1)The minimum FARLBOis observed in the present study.The result may be related to the property of the fuel and the configuration of the combustor.In the future,the minimum FARLBOof different combustors may be studied.2)The fuel distribution in the combustor is non-uniform,but chemical time scales are obtained by the zero-dimensional reactor PSR in CHEMKIN.The effects of heterogeneity will be considered by building a reactor net within PRZ in future works.

    Fig.24 Effects of Tinleton PPRZ,LBO.

    Fig.25 Effects of Tinleton τf,LBO.

    Fig.26 Effects of Tinleton τc,LBO.

    Fig.27 Da numbers under 16 LBO conditions.

    Fig.28 Relationship between Da and FAR.

    6.Conclusions

    In the present study,a hybrid LBO prediction method combining numerical simulation andDamodel is proposed based on bluffbody stabilized flames.Some important conclusions are obtained as follows:

    (1)Practical Reaction Zone(PRZ)is built based on concentration of OH radicals in the simulated flow field of the combustor.Flow time scales(τf)are calculated based onVPRZand the volume-averaged axial velocity of PRZ.Chemical time scales(τc)are calculated based on volume-averaged temperature,pressure and equivalence ratio of PRZ.

    (2)τfincreases with the increase ofMainlet,but decreases with the increase ofTinlet.τcincreases with the increase ofMainlet,but decreases with the increase ofTinlet.τfis mainly determined byVPRZ,and τcis determined byTPRZand ΦPRZ.

    (3)The average value ofDanumber under LBO conditions is 1.135,andDanumber under each LBO condition ranges from 0.673 to 1.351.Under LBO conditions,the flow time scale and chemical time scale are comparable,which indicates a critical state of combustion.

    (4)Danumber increases in an exponential way when fuel/air ratio of inlet flow increases,showing that flames become stable when the fuel supply increases.

    Acknowledgements

    This work has been carried out with the supports of National Key Research and Development Program of China(No.2016YFB0901402)and National Natural Science Foundation of China(No.51476170).

    亚洲成av片中文字幕在线观看| 美女福利国产在线| 不卡av一区二区三区| 好男人电影高清在线观看| av在线老鸭窝| 国产91精品成人一区二区三区 | 亚洲全国av大片| 亚洲精品第二区| 啦啦啦视频在线资源免费观看| 国产精品 欧美亚洲| 欧美精品av麻豆av| 精品久久久久久电影网| 久久这里只有精品19| 女人精品久久久久毛片| 日韩大片免费观看网站| 久久久久视频综合| 男男h啪啪无遮挡| 手机成人av网站| 搡老岳熟女国产| 精品人妻1区二区| 国产欧美日韩精品亚洲av| 亚洲人成电影观看| 黄色a级毛片大全视频| 亚洲精品一卡2卡三卡4卡5卡 | av网站免费在线观看视频| 亚洲av国产av综合av卡| 精品国产一区二区三区久久久樱花| 动漫黄色视频在线观看| 国产精品影院久久| 久久久精品94久久精品| 欧美久久黑人一区二区| 99九九在线精品视频| 国产成人欧美| 国产精品国产av在线观看| 国产精品自产拍在线观看55亚洲 | 人妻久久中文字幕网| 老司机午夜福利在线观看视频 | 成人黄色视频免费在线看| 欧美黄色淫秽网站| 亚洲 欧美一区二区三区| 亚洲专区国产一区二区| 国产高清视频在线播放一区 | 国产精品亚洲av一区麻豆| 亚洲精品国产一区二区精华液| 高清av免费在线| 啦啦啦在线免费观看视频4| 我要看黄色一级片免费的| 亚洲欧美清纯卡通| 国产免费一区二区三区四区乱码| 欧美激情久久久久久爽电影 | 日本vs欧美在线观看视频| 9色porny在线观看| 久久久久久久大尺度免费视频| 久久亚洲精品不卡| 一级片免费观看大全| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲中文日韩欧美视频| 在线观看人妻少妇| 国产男人的电影天堂91| 日本av手机在线免费观看| 亚洲熟女毛片儿| 久久亚洲国产成人精品v| 亚洲一区二区三区欧美精品| 久久国产精品人妻蜜桃| 亚洲全国av大片| 亚洲中文av在线| 精品国产乱码久久久久久小说| 一级毛片电影观看| 999精品在线视频| 久久久久精品人妻al黑| 视频在线观看一区二区三区| 亚洲国产中文字幕在线视频| 国产精品熟女久久久久浪| 色老头精品视频在线观看| 日韩 欧美 亚洲 中文字幕| 国产成人精品无人区| 大片电影免费在线观看免费| 9热在线视频观看99| 精品一品国产午夜福利视频| 国产色视频综合| 看免费av毛片| 国产亚洲av片在线观看秒播厂| 欧美日韩中文字幕国产精品一区二区三区 | 免费在线观看完整版高清| 免费人妻精品一区二区三区视频| 中国国产av一级| 动漫黄色视频在线观看| 成人亚洲精品一区在线观看| 老司机福利观看| 欧美日韩av久久| 午夜福利乱码中文字幕| 国产精品成人在线| 精品高清国产在线一区| 亚洲伊人色综图| 亚洲性夜色夜夜综合| 女人精品久久久久毛片| 狂野欧美激情性bbbbbb| 欧美激情极品国产一区二区三区| 操美女的视频在线观看| 欧美 亚洲 国产 日韩一| 男人添女人高潮全过程视频| 一区二区三区精品91| 久久久久久久久免费视频了| 亚洲少妇的诱惑av| 午夜福利在线观看吧| 亚洲欧洲精品一区二区精品久久久| 精品卡一卡二卡四卡免费| 国产亚洲精品第一综合不卡| 亚洲精品自拍成人| 日日爽夜夜爽网站| 青春草视频在线免费观看| 午夜福利在线观看吧| 久久综合国产亚洲精品| 成人国产一区最新在线观看| 永久免费av网站大全| av一本久久久久| 啦啦啦啦在线视频资源| 99国产综合亚洲精品| 精品欧美一区二区三区在线| 动漫黄色视频在线观看| 免费观看人在逋| 欧美日韩亚洲国产一区二区在线观看 | 人成视频在线观看免费观看| 高潮久久久久久久久久久不卡| 国产精品久久久久久精品古装| 亚洲精品一区蜜桃| av国产精品久久久久影院| √禁漫天堂资源中文www| 热re99久久国产66热| 免费日韩欧美在线观看| 下体分泌物呈黄色| www.999成人在线观看| 18禁黄网站禁片午夜丰满| 国产在视频线精品| 国产av一区二区精品久久| 亚洲精品日韩在线中文字幕| 精品国产乱码久久久久久小说| 成年人黄色毛片网站| 色视频在线一区二区三区| 国产高清视频在线播放一区 | 人人妻人人添人人爽欧美一区卜| 国产成人a∨麻豆精品| 无限看片的www在线观看| 亚洲精品在线美女| av欧美777| 大香蕉久久成人网| 精品久久久精品久久久| 精品一品国产午夜福利视频| 久久久精品国产亚洲av高清涩受| 天天躁日日躁夜夜躁夜夜| 啦啦啦免费观看视频1| 国产成人精品无人区| 制服人妻中文乱码| 国产欧美日韩一区二区三区在线| 久久精品人人爽人人爽视色| 777米奇影视久久| 国产视频一区二区在线看| 啦啦啦 在线观看视频| 18禁裸乳无遮挡动漫免费视频| 色婷婷久久久亚洲欧美| 精品人妻熟女毛片av久久网站| 中文字幕人妻熟女乱码| 老司机在亚洲福利影院| 亚洲人成77777在线视频| 成人三级做爰电影| 丝瓜视频免费看黄片| 国产成人精品无人区| 黄片播放在线免费| 美女中出高潮动态图| 欧美精品一区二区免费开放| 亚洲欧美一区二区三区久久| 欧美黑人欧美精品刺激| 一本色道久久久久久精品综合| 成年人午夜在线观看视频| 国产不卡av网站在线观看| 日韩有码中文字幕| 侵犯人妻中文字幕一二三四区| 亚洲av美国av| 人妻久久中文字幕网| 一区二区三区激情视频| 国产精品一区二区在线观看99| 日韩中文字幕视频在线看片| 亚洲国产成人一精品久久久| 亚洲av日韩在线播放| 久久人人爽av亚洲精品天堂| 精品国产一区二区三区四区第35| 亚洲五月婷婷丁香| 岛国在线观看网站| 色94色欧美一区二区| 一级片'在线观看视频| 在线观看www视频免费| 美女脱内裤让男人舔精品视频| www.精华液| 国产极品粉嫩免费观看在线| 亚洲午夜精品一区,二区,三区| 国产又色又爽无遮挡免| e午夜精品久久久久久久| 免费在线观看视频国产中文字幕亚洲 | 欧美日韩国产mv在线观看视频| 久久国产精品男人的天堂亚洲| 亚洲精品自拍成人| 在线观看舔阴道视频| 青春草视频在线免费观看| 69精品国产乱码久久久| 岛国在线观看网站| 19禁男女啪啪无遮挡网站| 亚洲avbb在线观看| 精品亚洲乱码少妇综合久久| 欧美日韩视频精品一区| 老司机亚洲免费影院| 在线精品无人区一区二区三| 夫妻午夜视频| 午夜影院在线不卡| 国产麻豆69| 在线十欧美十亚洲十日本专区| 精品一品国产午夜福利视频| av片东京热男人的天堂| 老汉色∧v一级毛片| 狠狠狠狠99中文字幕| 亚洲情色 制服丝袜| 新久久久久国产一级毛片| 亚洲激情五月婷婷啪啪| 一本综合久久免费| 捣出白浆h1v1| 中文字幕最新亚洲高清| 在线观看免费视频网站a站| 久久精品久久久久久噜噜老黄| 99国产精品一区二区三区| 国产一区二区三区综合在线观看| 熟女少妇亚洲综合色aaa.| 少妇粗大呻吟视频| 亚洲人成电影观看| 日韩 亚洲 欧美在线| 国产精品自产拍在线观看55亚洲 | 日韩一卡2卡3卡4卡2021年| 日韩 亚洲 欧美在线| 黑丝袜美女国产一区| 亚洲情色 制服丝袜| 国产xxxxx性猛交| 亚洲人成电影免费在线| 妹子高潮喷水视频| a级毛片在线看网站| 91九色精品人成在线观看| 18禁国产床啪视频网站| 亚洲欧洲日产国产| 午夜福利视频在线观看免费| 日本wwww免费看| 一二三四在线观看免费中文在| 亚洲精品美女久久av网站| 自拍欧美九色日韩亚洲蝌蚪91| 啦啦啦啦在线视频资源| 麻豆av在线久日| 欧美日韩国产mv在线观看视频| 国产精品国产av在线观看| 最新在线观看一区二区三区| 欧美+亚洲+日韩+国产| 久久久久国内视频| a级毛片黄视频| 国产一区二区三区在线臀色熟女 | 日韩,欧美,国产一区二区三区| 亚洲精品乱久久久久久| www.自偷自拍.com| 99香蕉大伊视频| 久久国产精品影院| 久久久久久久久久久久大奶| www日本在线高清视频| 大香蕉久久网| 国产99久久九九免费精品| 国产av一区二区精品久久| 十分钟在线观看高清视频www| 国产日韩欧美视频二区| 日韩人妻精品一区2区三区| 丰满少妇做爰视频| 五月开心婷婷网| 中国国产av一级| 男女无遮挡免费网站观看| 老司机午夜十八禁免费视频| 亚洲精品国产av蜜桃| 一区二区三区精品91| 美女国产高潮福利片在线看| 精品免费久久久久久久清纯 | 亚洲欧洲精品一区二区精品久久久| 精品亚洲乱码少妇综合久久| 久久精品久久久久久噜噜老黄| 国产精品一区二区精品视频观看| 日韩视频一区二区在线观看| 脱女人内裤的视频| 99精品欧美一区二区三区四区| 可以免费在线观看a视频的电影网站| 999精品在线视频| 午夜福利在线观看吧| 免费日韩欧美在线观看| 亚洲欧美精品自产自拍| 高清欧美精品videossex| 岛国毛片在线播放| 啦啦啦在线免费观看视频4| 欧美激情 高清一区二区三区| 亚洲免费av在线视频| 在线亚洲精品国产二区图片欧美| 亚洲人成电影观看| 国产日韩欧美视频二区| 亚洲av欧美aⅴ国产| 亚洲专区字幕在线| 久久人妻福利社区极品人妻图片| 日本av手机在线免费观看| 两个人免费观看高清视频| 大码成人一级视频| 久久久国产欧美日韩av| 日本vs欧美在线观看视频| 日韩中文字幕欧美一区二区| 纵有疾风起免费观看全集完整版| 色播在线永久视频| 日韩三级视频一区二区三区| 我的亚洲天堂| 男女床上黄色一级片免费看| 50天的宝宝边吃奶边哭怎么回事| av视频免费观看在线观看| 久久久久久亚洲精品国产蜜桃av| 国产亚洲一区二区精品| 18在线观看网站| 老司机在亚洲福利影院| 色精品久久人妻99蜜桃| av免费在线观看网站| 狠狠狠狠99中文字幕| 亚洲欧美清纯卡通| 亚洲欧美精品自产自拍| 久久久久网色| 男女下面插进去视频免费观看| 国产免费一区二区三区四区乱码| 亚洲精品粉嫩美女一区| 日日摸夜夜添夜夜添小说| 99国产精品免费福利视频| 久久午夜综合久久蜜桃| 人妻一区二区av| 男女国产视频网站| 91麻豆精品激情在线观看国产 | 免费在线观看完整版高清| 男女高潮啪啪啪动态图| 1024香蕉在线观看| 欧美久久黑人一区二区| 亚洲成av片中文字幕在线观看| 国产亚洲欧美精品永久| 国产淫语在线视频| 国产成人啪精品午夜网站| 人人妻,人人澡人人爽秒播| 亚洲激情五月婷婷啪啪| 欧美少妇被猛烈插入视频| 亚洲综合色网址| 欧美日韩精品网址| 日本精品一区二区三区蜜桃| 国产国语露脸激情在线看| netflix在线观看网站| 啦啦啦免费观看视频1| 正在播放国产对白刺激| 亚洲,欧美精品.| 正在播放国产对白刺激| 日韩精品免费视频一区二区三区| 99热国产这里只有精品6| 不卡一级毛片| 日韩有码中文字幕| 纯流量卡能插随身wifi吗| 久久香蕉激情| 嫩草影视91久久| 操出白浆在线播放| 欧美+亚洲+日韩+国产| 操出白浆在线播放| 极品少妇高潮喷水抽搐| 一本色道久久久久久精品综合| 不卡一级毛片| 在线观看免费午夜福利视频| 亚洲全国av大片| 女警被强在线播放| 99国产极品粉嫩在线观看| 日本一区二区免费在线视频| 亚洲av成人不卡在线观看播放网 | 精品少妇一区二区三区视频日本电影| 国产有黄有色有爽视频| 一本大道久久a久久精品| 国产亚洲精品久久久久5区| 人妻人人澡人人爽人人| 亚洲国产精品一区二区三区在线| 女人爽到高潮嗷嗷叫在线视频| 欧美日韩中文字幕国产精品一区二区三区 | 欧美人与性动交α欧美精品济南到| 秋霞在线观看毛片| 国产成人影院久久av| 亚洲三区欧美一区| 亚洲色图 男人天堂 中文字幕| 国产精品免费视频内射| 欧美日韩成人在线一区二区| 老汉色∧v一级毛片| 国产精品久久久久久人妻精品电影 | 久久国产亚洲av麻豆专区| 老熟女久久久| 这个男人来自地球电影免费观看| 黄色毛片三级朝国网站| av在线老鸭窝| 99精品久久久久人妻精品| 丝袜脚勾引网站| 亚洲视频免费观看视频| 精品国产乱码久久久久久小说| 热99国产精品久久久久久7| 久久人妻福利社区极品人妻图片| 咕卡用的链子| 在线看a的网站| 亚洲国产精品成人久久小说| 亚洲精品久久久久久婷婷小说| 一二三四在线观看免费中文在| 美女大奶头黄色视频| 久久国产精品影院| 国产av又大| 久久久久久久大尺度免费视频| 中国美女看黄片| 国产老妇伦熟女老妇高清| 亚洲av电影在线进入| 黄片播放在线免费| 搡老岳熟女国产| 99久久精品国产亚洲精品| 亚洲国产成人一精品久久久| 久久精品亚洲av国产电影网| 国产老妇伦熟女老妇高清| 黄色视频不卡| 亚洲情色 制服丝袜| 国产精品香港三级国产av潘金莲| 男人操女人黄网站| 777久久人妻少妇嫩草av网站| 国产一区二区三区在线臀色熟女 | 精品国产乱码久久久久久男人| 精品高清国产在线一区| 精品免费久久久久久久清纯 | 国产成人免费无遮挡视频| 国产一区二区激情短视频 | 亚洲一码二码三码区别大吗| 亚洲五月色婷婷综合| 1024视频免费在线观看| 国产精品久久久久久精品古装| 欧美日韩亚洲国产一区二区在线观看 | 亚洲一卡2卡3卡4卡5卡精品中文| 国产无遮挡羞羞视频在线观看| 中亚洲国语对白在线视频| 亚洲欧美日韩另类电影网站| 精品国产乱码久久久久久小说| 国产一级毛片在线| 在线观看免费午夜福利视频| 国产成人系列免费观看| 精品一区在线观看国产| 岛国毛片在线播放| 一区二区三区乱码不卡18| 免费女性裸体啪啪无遮挡网站| 中文字幕精品免费在线观看视频| 中文字幕人妻丝袜制服| 下体分泌物呈黄色| 制服诱惑二区| 黑人巨大精品欧美一区二区蜜桃| 三级毛片av免费| 成人国产一区最新在线观看| 91国产中文字幕| 国产老妇伦熟女老妇高清| 国产免费现黄频在线看| 精品国内亚洲2022精品成人 | 国产成+人综合+亚洲专区| 久久精品国产亚洲av高清一级| 国产精品久久久久成人av| 亚洲国产欧美一区二区综合| 在线观看www视频免费| 人人妻人人澡人人看| 欧美精品高潮呻吟av久久| 夫妻午夜视频| 99精品欧美一区二区三区四区| 一级毛片女人18水好多| 久久女婷五月综合色啪小说| www日本在线高清视频| 十八禁人妻一区二区| 黄色 视频免费看| 亚洲av男天堂| 两个人免费观看高清视频| 久久热在线av| 国产一区二区三区在线臀色熟女 | 亚洲人成电影观看| www.自偷自拍.com| 亚洲avbb在线观看| 久久久久精品人妻al黑| 欧美 亚洲 国产 日韩一| 母亲3免费完整高清在线观看| 欧美精品亚洲一区二区| 别揉我奶头~嗯~啊~动态视频 | 婷婷丁香在线五月| 欧美午夜高清在线| 亚洲国产欧美网| 男女免费视频国产| 精品国产超薄肉色丝袜足j| 日韩制服骚丝袜av| 青春草视频在线免费观看| 我的亚洲天堂| 精品国产一区二区三区久久久樱花| 亚洲自偷自拍图片 自拍| 国产精品久久久久久人妻精品电影 | 国产视频一区二区在线看| 69精品国产乱码久久久| 午夜影院在线不卡| 9色porny在线观看| 一个人免费在线观看的高清视频 | 欧美激情极品国产一区二区三区| 91精品伊人久久大香线蕉| 狠狠狠狠99中文字幕| 久久久精品区二区三区| 亚洲av男天堂| 欧美老熟妇乱子伦牲交| 国产成人a∨麻豆精品| 欧美在线一区亚洲| 国产无遮挡羞羞视频在线观看| 操出白浆在线播放| 亚洲精品久久成人aⅴ小说| 精品一区二区三区av网在线观看 | 欧美精品人与动牲交sv欧美| 美女视频免费永久观看网站| 最新在线观看一区二区三区| 欧美老熟妇乱子伦牲交| 最近最新免费中文字幕在线| 亚洲精品国产区一区二| 一本色道久久久久久精品综合| 男人添女人高潮全过程视频| 50天的宝宝边吃奶边哭怎么回事| 午夜福利免费观看在线| 亚洲自偷自拍图片 自拍| 成人手机av| videos熟女内射| 久久久久精品国产欧美久久久 | 国产人伦9x9x在线观看| 成在线人永久免费视频| 国产精品久久久久久精品电影小说| 嫁个100分男人电影在线观看| 亚洲 国产 在线| 一级片'在线观看视频| 亚洲精品美女久久av网站| 一本大道久久a久久精品| 国产精品偷伦视频观看了| 亚洲 国产 在线| 2018国产大陆天天弄谢| 欧美日本中文国产一区发布| 99国产精品一区二区三区| 国产一区二区三区av在线| 一个人免费看片子| 人人澡人人妻人| 欧美xxⅹ黑人| 悠悠久久av| 大片电影免费在线观看免费| 欧美激情极品国产一区二区三区| 又黄又粗又硬又大视频| 国产黄频视频在线观看| 亚洲欧美色中文字幕在线| 丰满迷人的少妇在线观看| 欧美激情高清一区二区三区| 在线观看www视频免费| 一本—道久久a久久精品蜜桃钙片| 另类亚洲欧美激情| 真人做人爱边吃奶动态| 交换朋友夫妻互换小说| 狠狠精品人妻久久久久久综合| 99久久精品国产亚洲精品| 精品乱码久久久久久99久播| 我要看黄色一级片免费的| 免费观看a级毛片全部| 国产在线观看jvid| 91老司机精品| 纵有疾风起免费观看全集完整版| 激情视频va一区二区三区| av又黄又爽大尺度在线免费看| 另类亚洲欧美激情| 在线观看舔阴道视频| 欧美亚洲 丝袜 人妻 在线| 99国产极品粉嫩在线观看| 青春草亚洲视频在线观看| 国产男女超爽视频在线观看| 天堂俺去俺来也www色官网| 日本撒尿小便嘘嘘汇集6| 亚洲成人国产一区在线观看| 不卡av一区二区三区| 99精品久久久久人妻精品| 久久精品国产亚洲av香蕉五月 | 成人18禁高潮啪啪吃奶动态图| 国产一区二区激情短视频 | 成年人黄色毛片网站| 国产xxxxx性猛交| 国产成人a∨麻豆精品| 色精品久久人妻99蜜桃| 精品亚洲成a人片在线观看| 制服人妻中文乱码| 亚洲av日韩在线播放| 亚洲精品国产av蜜桃| 男女国产视频网站| 国产亚洲精品久久久久5区| 飞空精品影院首页| 国产亚洲午夜精品一区二区久久| 岛国毛片在线播放| 黑人巨大精品欧美一区二区mp4| 别揉我奶头~嗯~啊~动态视频 | 亚洲av电影在线观看一区二区三区| 一个人免费在线观看的高清视频 | 久久中文字幕一级| 丰满饥渴人妻一区二区三| 亚洲国产欧美在线一区| 国产欧美日韩精品亚洲av| 一级a爱视频在线免费观看| 日韩一卡2卡3卡4卡2021年| 91精品三级在线观看| 男人爽女人下面视频在线观看| 欧美日韩视频精品一区| 精品国产一区二区三区四区第35| 国产av国产精品国产| 1024视频免费在线观看| 精品亚洲成a人片在线观看| 一区二区三区四区激情视频| av超薄肉色丝袜交足视频| 国产野战对白在线观看|