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

    An integrated model for predicting the flame propagation in crimped ribbon flame arresters☆

    2018-06-29 09:15:36ZhengWangBingSunQingshanHuangFuhuaJiang
    關(guān)鍵詞:創(chuàng)城市容市政府

    Zheng Wang ,Bing Sun ,Qingshan Huang ,Fuhua Jiang *

    1 Qingdao Institute of Bioenergy and Bioprocess Technology,Chinese Academy of Sciences,Qingdao 266101,China

    2 State Key Laboratory of Safety and Control for Chemicals,SINOPEC Research Institute of Safety Engineering,Qingdao 266000,China

    3 Key Laboratory of Green Process and Engineering,Institute of Process Engineering,Chinese Academy of Sciences,Beijing 100190,China

    1.Introduction

    Crimped ribbon flame arresters are usually installed in the pipelines of the flammable and combustible gas or on the top of the storage tank[1,2].They allow the gas to pass through and prevent the passage of the burning flame.Improper design of the flame arresters may give rise to flameproof failure,which can result in a fire disaster and even lead to a great loss of lives and properties.The design method of numerical simulation is promising;however,its reliability and accuracy are greatly influenced by the mathematical model.An accurate mathematical model should be set up to help designing the crimped ribbon flame arresters.

    Many experimental studies[3-5]have been conducted on the flame arresters,but numerical studies are less.To the best of our knowledge,only Wanget al.[6]did a pioneering experimental and numerical study to investigate flame propagation and quenching in the microchannel of crimped ribbon flame arresters.There are few studies on the microchannel of flame arrester;however,many studies[7-9]have been done on the micro-combustor.Since both microchannels of flame arresters and micro-combustors are micro-size space involving flame propagation,some important conclusions obtained in the microcombustors can be learned and extended for studying the microchannel in flame arresters.The laminar model was usually assumed in the predictions of micro-combustor[10-12].Wanget al.[6]also used the laminar model to predict the flame characteristics in the microchannel of crimped ribbon flame arrester.The fuel enters the micro-combustor with tens of centimeters per second,but the fluid velocity at the entrance of the microchannel of flame arresters is more than tens of meters per second.That is to say,the fluid velocity at the inlet of the microchannel in the flame arresters is more than hundreds of times greater than that of the micro-combustor.Additionally,because the triangle microchannel yields the highest flame turbulence intensity[13],mixing between the components is enhanced.Moreover,violent chemical reaction flow is usually in a chaotic flow pattern[14,15].Thus,whether the laminar model is suitable to predict the flame propagation in the microchannel of the flame arresters or not,especially for the triangle mircochannel of crimped ribbon flame arresters,is worth discussing,and more studies on the flow regime of the flame propagation in the microchannel of the flame arresters need to be carried out.

    Since there are violent combustion reactions in the crimped ribbon flame arresters,the chemical kinetics mechanism of flammable mixtures is very important to establish the accurate mathematical model.Propane has been widely used in applications such as gas turbines and water heating due to advantages of cost savings,widely available and convenient storage[16-19].Owing to its importance,oxidation of propane has therefore been intensively investigated by experiments since 1980s[20-24].Scholars have proposed a variety of kinetics mechanisms for the system of propane-air,and some numerical studies[25,26]containing the chemical kineticsmechanisms of propane-air mixtures have been conducted in recent years.For example,Wanget al.[6]utilized a single step mechanism for propane-air to study the relationship between the flameproof velocity and the arrester structure;Mansouriet al.[27]employed a simplified chemical kinetics mechanism of propane-air,which consists of 3 chemical reactions and 5 species,to investigate the effect of hydrogen addition and equivalence ratio on characteristics of premixed propane-air flame in a burner;Konakovet al.[28]used a one step global mechanism to investigate the flame shape and its long term stability;Luoet al.[29]adopted a simplified propane oxidation reaction kinetics to investigate the transition from deflagration to detonation in an explosive vessel with three barriers.Although great achievements have been made,few comparative studies of the sensitivity of the chemical kinetics mechanisms on the accuracy of numerical prediction have been taken yet.To the best of the authors' knowledge,only Gutkowski[30]compared the performance between the one step and two steps propane-air chemical kinetics mechanisms on the flame behaviour under quenching conditions for the flame propagating of propane-air mixtures in square cross-sections channels.Therefore,further studies of the sensitivity of chemical kinetics mechanisms on the prediction of flame quenching in the flame arresters should be undertaken.Since it is difficult to determine whether the flame is extinguished or not in the microchannel with the numerical method,scholar gave a very wide temperature range[6]as the judging criterion of quenching.Since it is not suitable for industrial applications,a scientific criterion for the flame quenching in numerical simulation should be established.

    In the present work,an integrated mathematical model for the prediction of the reactive flow under high temperature in the microchannel of the crimped ribbon flame arresters isestablished,and the sensitivities of four chemical kinetics mechanisms of propane-air on the accuracy of flame propagation characteristics are investigated.An accurate chemical kinetics mechanism is identified and recommended for predicting the violent combustion reaction in the microchannel of the crimped ribbon flame arresters.The fluid flow regime of the flame propagation in the microchannel of the crimped ribbon flame arresters under the deflagration and detonation condition is determined,and a judging criterion of the flame quenching for the numerical method is proposed.The details of the flame propagation in the triangle microchannel including the distributions of temperature and species mass fraction are captured.The deep understanding of propane-air combustion and extinction characteristics in the microchannel is also obtained.It is indicated that the integrated mathematical model validated in this work can be used as a reliable tool for the modeling,designing and selecting of the industrial crimped ribbon flame arresters.

    2.Physical Structure and Mathematical Model

    2.1.Physical structure

    A crimped ribbon flame arrester,as described in Fig.1(a),consists of an arrester housing,an arrester element and other associated fittings.When the flame passes through the arrester,it penetrates into the microchannels of the arrester element which consisted of many isosceles triangular microchannels.Meanwhile,the flame is split into small flamelets.Since the microchannel could efficiently absorb heat,the flame could be quenched successfully.Hence,the microchannelis the critical part to quench the flame.The microchannel,as presented in Fig.1(b),was investigated as the physical structure in this work.

    Two microchannels(physical structures A and B)with different structure parameters were extracted from the literature[6]to verify the mathematical model in this work,and the basic parameters of the microchannel were illustrated in Table 1.

    Table 1 The basic parameters of physical structures A and B

    The arrester housing has an expansion ratio β [β =D/d,wheredis the diameter of the pipeline,andDis the diameter of the arrester element as shown in Fig.1(a)].The flame enters the pipeline with the velocityV.Since the mass is conserved,the flame velocity at the entrance of the microchannel,uin,can be calculated asV/β2.

    2.2.Mathematical model

    For the purpose of reducing the computational cost and get an efficient and reliable solution,steady mathematical models are taken here[31].

    Fig.1.The sketch of the crimped ribbon flame arrester:(a)crimped ribbon flame arrester;(b)the microchannel.

    2.2.1.Governing equations

    The conservations of continuity,momentum,and energy are kept in control volume,and these mathematical models can be expressed with the following equations.

    (1)Continuity equation

    The gas density is calculated using the ideal gas law.

    (2)Momentum equation

    wherepis the static pressure,ρg is the gravitational body force;is the stress tensor,and it is given by

    where μ is the laminar viscosity;μtis the turbulent viscosity;Iis the unit tensor.For laminar simulation,ρk=0 and μt=0.

    A standardk-ε equations is used here to solve the turbulent flow[32].

    where

    In the above equations,μ is the molecular viscosity.Gkrepresents the generation of turbulent kinetic energy due to the mean velocity gradients.C1ε=1.44,C2ε=1.92,σk=1.0,σε=1.3 andCμ=0.09.

    (3)Energy equation

    wherekeffis the effective conductivity.keff=k+kt,wherektis turbulent thermal conductivity andkt=Pr=0.85.Eis the total energy,tand it is defined as

    whereh issensible enthalpy andh=∑jYjhj.For individual species,

    whereTref=298.15 K.

    The gas viscosity,specific heat and thermal conductivity are calculated as a mass fraction-weighted average of all species.The specific heat of each species is calculated using a piecewise polynomial fit of temperature.

    2.2.2.Modeling of chemical reaction

    The mixing and transport of chemical species were modeled by solving species transport conservation equations which described the convection,diffusion,and reaction sources for each component species.The local mass fraction of each speciesYi,is predicted through the solution of a convection-diffusion equation for theith species.The conservation equation takes the following form:

    whereRiis the net production rate of speciesiby chemical reaction.

    In Eq.(10),Jiis the diffusion flux of speciesi,which arises due to gradients of concentration and temperature.The dilute approximation(also called Fick's law)is utilized to model mass diffusion due to concentration gradients,in which the diffusion flux can be written as,

    whereDi,mis the mass diffusion coefficient or speciesiin the mixture.DT,iis the thermal diffusion coefficient,andSct=0.7.

    Though detailed chemical kinetic calculations are computationally expensive,in order to predict precisely the effect of turbulent fluctuations on the chemical reaction,detailed Arrhenius chemical kinetics are incorporated in the turbulent flame by the Finite Rate/Eddy Dissipation Model(FR/EDM).

    For the finite rate model,Ri,ris given by

    The rate constant for reactionr,kf,ris computed using the Arrhenius expression:

    The EDM is based on the hypothesis that the chemical reaction is fast in relation to the transport processes of the turbulent fluid flow.The reaction rate is assumed to be proportional to a mixing time defined by the turbulent kinetic energy(k)and the specific dissipation rate(ε).Based on the contribution of Magnussen and Hjertager[33],the net production rate of speciesidue to reactionr,Ri,r,is given by the smaller of the two expressions below

    whereAis an empirical constant equal to 4.0;Bis an empirical constant equal to 0.5.

    In Eqs.(15)and(16),the chemical reaction rate is governed by the large-eddy mixing time scale,k/ε,proposed in the eddy-breakup model of Spalding[34].The finite-rate/eddy-dissipation model,in which both the Arrhenius(Eq.(13)and eddy-dissipation(Eqs.(15)and(16))reaction rates are calculated,were employed.The netreaction rate is taken as the minimum of these two rates.For the laminar flow,the Arrhenius rate is taken.

    Four different mechanisms(Jones and Lindstedt[23],Konakovet al.[28],Luoet al.[29],Mansouriet al.[27]),as presented in Tables 2-5,respectively,were applied to simulate propane-air combustion in the microchannel of crimped ribbon flame arresters.

    Table 2 Oxidation kinetics of propane in Jones and Lindstedt[23]

    Table 3 Oxidation kinetics of propane in Konakov et al.[28]

    Table 4 Oxidation kinetics of propane in Luo et al.[29]

    Table 5 Oxidation kinetics of propane in Mansouri et al.[27].

    2.3.Boundary conditions

    The boundary conditions are very important for getting correct simulation results of fluid flow[35,36].On the wall,no slip boundary condition is applied and the standard wall functions are adopted for all species at the solid wall[37-39],and a zero-gradient(zero- flux)boundary condition is set for all species.The temperature of the wall is isothermal(Twall=300 K).The boundary conditions at the inlet are set as a fixed inlet velocity,which is normal to the boundary of the inlet.The flame velocity at the entrance of the microchannel,uin,isV/β2.The inlet conditions for the turbulent kinetic energy and the dissipation rate are estimated by[31,40]

    whereDHis the hydraulic diameter of the inlet.

    The inlet temperature of the species was chosen to be 2300 K and the mass fractions of the reactants(C3H8and O2)were 0.064 and 0.232,respectively.

    The outlet of the microchannel is treated as a fully developed boundary,and all the variables are subject to

    Second order upwind scheme is used to discretize convective terms of all the governing equations,and the traditional SIMPLE(Semi-Implicit Method for Pressure Linked Equations)algorithm is adopted to deal with the pressure-velocity coupling.A 3D segregated solver with the under-relaxation method is employed to solve the conservation equations.The simulation convergence is judged by the residuals of the whole governing equations and the convergence criteria for the residuals are set to be 1.0×10-6for all the equations.

    3.Results and Discussion

    The experimental data in the contribution of Wanget al.[27]were used to verify the mathematical model.First,grid independence was conducted to ensure the results independent on the grid size.The flow pattern of the flame propagation in the microchannel was determined.Second,the mathematical model was verified and the most accurate chemical kinetics mechanism was screened.Third,temperature and species mass fraction profiles predicted by different chemical kinetics mechanisms for the physical structure A were compared.Finally,the temperature and species mass fraction profiles predicted by the most accurate chemical kinetics mechanism for the physical structure B were presented.

    3.1.Grid independence

    A fine grid size is essential for computational accuracy during the simulation,and correspondingly the computational cost is sharply increased.Therefore,a balance should be taken between the accuracy and the computational time.The grid independence study was performed for both physical structures A and B.Set the physical structure A as an example to illustrate this process.The predicted static temperature along the axial direction of the microchannel at the centerline is plotted in Fig.2.Preliminary computations with two folds finer grids were also conducted.As presented in Fig.2,there is no significant difference between the results obtained by these two sets of grids.Accordingly,the mesh intervals of 0.04 mm and 0.1 mm in the cross-section and the axial direction respectively are employed here to ensure higher accuracy as well as to reduce the computational time.

    Fig.2.Study of grid independency.

    3.2.Flow regime determination

    Fig.3.μt/ μat the center surface:(a)the physical structure A;(b)the physical structure B.

    The flow regime of the flame propagation in the microchannel is worth discussing.Because the inlet velocity is much smaller than that of the microchannel in the crimped ribbon flame arresters,some scholars utilized a laminar model to predict the flame characteristics in the micro-combustors.Since the inlet velocity of the flame entering the microchannel is up to tens of meters per second and there is violent combustion in the microchannel of the crimped ribbon flame arresters,the fluid flow regime may be turbulent.When the inlet velocity is equal to the flameproof velocity,the ratio of μt/ μ in the center surface of physical structures A and B are depicted in Fig.3.It is shown that the value ofμt/ μ is higher than 1 in the vast majority of zones.It is indicated that the flow pattern is turbulent in the microchannel of the crimped ribbon flame arresters.

    In addition,when the inlet velocity is equal to the flameproof velocity,profiles of the predicted outlet temperature at the centerline in the axial direction of the physical structure A with the turbulent and laminar model are presented in Fig.4.It can be seen that the temperature predicted by the laminar model is much higher than that predicted by the turbulent one.The temperature profile predicted by the turbulent model is much closer to the experimental results.

    Fig.4.Predicted outlet temperature profiles at the axial centerline of the physical structure A.

    For the above two reasons,the flow regime of the flame propagation in the microchannel under the deflagration and detonation conditions is inferred to be turbulent.

    3.3.Comparison with the experimental results

    Three operating conditions(V<Vf,V=Vf,V>Vf, whereVandVfrepresent the flame velocity in the inlet of the pipeline and the flameproof velocity,respectively)were utilized to verify the mathematical model established here.The experimental data ofVandVfare extracted from the literature[6].

    In this work,the auto-ignition temperature is adopted as the judging criterion,which is recommended as the quenching criterion by Zhouet al.[41]in the experiments.It is noteworthy that the auto-ignition temperature is the minimum temperature required to ignite the flammable mixture without a spark or flame being present.In other words,if the temperature of the hot products is higher than the auto-ignition temperature,the unburned gas behind the flame arrester can be ignited and thus a flame reemerges,which means the flame is failed to be quenched in the flame arresters.In contrast,if the outlet temperature is lower than the auto-ignition temperature,the flame could be quenched successfully.However,the temperature is uneven in the outlet,if one part or even one point of the outlet temperature is higher than the auto-ignition temperature,the flammable mixture can also be reignited,even though the temperature of all the other parts in the outlet are lower than the auto-ignition temperature.Therefore,the highest temperature value at the outlet surface is chosen here for comparing with the auto-ignition temperature of propane to verify the reliability of the mathematical model.

    Fig.5.Predicted outlet temperatures for the physical structure A.

    Fig.6.Predicted outlet temperatures for the physical structure B.

    Figs.5 and 6 present the predicted outlet temperatures with different inlet velocity and chemical kinetics mechanisms for physical structures A and B,respectively.As illustrated in Fig.5,for the deflagration flame,whenV<Vf, the predicted outlet temperature is lower than the auto-ignition temperature,and vice versa.For the detonation flame,the predicted results of this integrated model show the same performance,as depicted in Fig.6.When the flame velocity is equal to the flameproof velocity[6],the predicted outlet temperature using the four kinetic models aforementioned is approximately equal to the auto-ignition temperature of propane-air with the deviation range of 0.11%-12.87%and 1.23%-8.92%for physical structures A and B,respectively.It is shown that the kinetics mechanisms are all good for industrial design if a turbulent model is applied.However,different kinetics mechanisms bring about different degrees of deviations.Therefore,with the aim of getting the reliable predicted results,it is essential to find the most accurate one to be integrated in the model of fluid flow.The computed values of deviations are depicted in Table 6.

    相豐委員的問題很直接:“近幾年市政府及有關(guān)部門在違建治理方面力度大、手段硬,效果明顯,但個別區(qū)域違建現(xiàn)象仍然存在。創(chuàng)城成功之后,一些區(qū)域占道經(jīng)營、店外經(jīng)營現(xiàn)象有所反彈。對這些違法行為,條例都有明確規(guī)定,處罰措施也很具體。請問執(zhí)法部門如何加大執(zhí)法監(jiān)管力度,加強違建治理,規(guī)范經(jīng)營行為,保持市容市貌的整潔有序?”

    Table 6 Deviations between the predicted outlet temperature and the auto-ignition temperature of propane-air

    The deviation of prediction may be caused by the following reasons.First,the turbulent model is not accurate enough,and more studies are demanded.Second,there is a deviation between the inlet velocity adopted in the numerical simulation and the actual value in experiments.Since the inlet velocity of the microchannel is too difficult to measure due to the tiny dimension of the microchannel in the experiments,it is postulated as a uniform value ofV/β2in the numerical study.However,in the real system the velocity at the inlet of each microchannel has a slight difference from the postulated value because the velocity is uneven resulted from the influence of the physical structure of the arrester housing.As a consequence,these differences induce a deviation between the predicted outlet temperature and the autoignition temperature of propane-air.

    Fig.7.Temperature distribution in the physical structure A with different kinetics:(a)Jones and Lindstedt[23];(b)Konakov et al.[28];(c)Luo et al.[29];(d)Mansouri et al.[27].

    Although there are some deviations between the predicted values and the experimental values,the deviations are all within the range for the engineering applications,and it is indicated that the integrated model established here is reliable.The deviations presented in Table 6 are different for these four chemical kinetics mechanisms.It should be noted that though the deviation is small in the relative proportion,the absolute value is much higher.For example,the maximum deviation for the physical structure A predicted by the mechanism of Jones and Lindstedt[23]is-12.87%,and it means the predicted outlet temperature is 93.05 K lower than the auto-ignition temperature of propane-air.This indicates the chemical kinetics mechanism integrated in the mathematical model has great influence on the accurate prediction of the flame propagation in the microchannel.Since the outlet temperature predicted by the chemical kinetics mechanism of Mansouriet al.[27]among the four kinetics at the flameproof velocity has the relatively small deviation with the auto-ignition temperature of propaneair for both physical structures A and B,it is identified as the most accurate one to be integrated with the fluid flow model to predict the flame propagation in the microchannel.Consequently,it can be applied for the further industrial applications.

    3.4.Results for the physical structure A

    The integrated mathematical model could predict the temperature distribution and species mass fractions distribution accurately for both physical structures A and B,and the results indicate the integrated model is robust.Since the effect of chemical kinetics mechanisms on the distributions of temperature and the species mass fraction in two physical structures(A and B)is similar,only the results in the physical structure A are compared.In addition,to facilitate the comparison of the species mass fraction distribution predicted by different chemical kinetics mechanisms,only the common components(C3H8,H2O,CO2and H2O)in the four mechanisms are revealed.Some specific species(CO,H2)in the chemical kinetics mechanisms of Jones and Lindstedt[23],Luoet al.[29]and Mansouriet al.[27]are not presented.

    3.4.1.Effect of chemical kinetics mechanism on temperature

    In order to show vividly the simulation results,the center surface at the axial direction in the physical structure A is extracted.

    Fig.7 shows the temperature distribution predicted by different chemical kinetics mechanisms in the center surface of the physical structure A when the inlet velocity is equal to the flameproof velocity.Fig.7(a)-(d)represents the predicted results by chemical kinetics mechanisms from Jones and Lindstedt[23],Konakovet al.[28],Luoet al.[29]and Mansouriet al.[27],respectively.It is shown that the temperature rises first in all simulations because the flammable gas releases a significant amount of heat through combustion.Then it drops greatly due to heat loss from the cold wall.However,different chemical kinetics mechanisms give different temperature distributions.Since there are 6 steps of reactions and 2 steps of them are reverse reactions that absorbed the heat,the temperature distribution predicted by the chemical kinetics mechanism from Luoet al.[29]is more uniform than the others.In contrast,there is only 1 step reverse reaction in the model from Mansouriet al.[27],and no reverse reaction in the other two chemical kinetics mechanisms and hence a relatively sharp change of temperature is observed.

    3.4.2.Effect of chemical kinetics mechanism on species mass fraction distribution

    The comparison of species mass fraction distributions predicted by different chemical kinetics mechanisms for the physical structure A at the flameproof velocity is presented in Figs.8-11,where(a)-(d)represent the chemical kinetics mechanism from Jones and Lindstedt[23],Konakovet al.[28],Luoet al.[29]and Mansouriet al.[27],respectively.

    Fig.8.Mass fraction of C3H8 for the physical structure A:(a)Jones and Lindstedt[23];(b)Konakov et al.[28];(c)Luo et al.[29];(d)Mansouri et al.[27].

    Fig.9.Mass fraction of O2 for the physical structure A:(a)Jones and Lindstedt[23];(b)Konakov et al.[28];(c)Luo et al.[29];(d)Mansouri et al.[27].

    As revealed in Fig.8,the mass fraction of C3H8is the highest at the inlet of the microchannel,and then it decreases drastically due to consumption of propane as the chemical reaction goes on.Finally,mass fraction of C3H8tends to be stable because the stop of the reactions in the microchannel of the crimped ribbon flame arresters when the flame is quenched.This is consistent with the law of the propane combustion.

    3.4.2.2.Reactant of O2.The mass fractions for reactant of O2predicted by different chemical kinetics mechanisms are presented in Fig.9.

    As shown in Fig.9,the mass fraction of O2is the highest at the inlet of the microchannel,and then it decreases due to the reaction of propane oxidation.The mass fraction of O2is different along the axial direction due to the different consumption rates of O2.For the Konakovet al.[28],O2is consumed quickly due to the fact that C3H8reacts with O2to produce CO2and H2O.However,for the chemical kinetics mechanisms from Luoet al.[29]and Mansouriet al.[27],O2is consumed smoothly because C3H8is oxidized to CO at first,and then to CO2.The mass fraction of O2predicted by the chemical kinetics mechanism from Jones and Lindstedt[23]is bigger than the other three,as a result of the reason is that CO,which needs less O2,is the main product.

    3.4.2.3.Product of CO2.The mass fractions for product of CO2predicted by different chemical kinetics mechanisms are revealed in Fig.10.

    As presented in Fig.10,there is no product of CO2at the inlet of the microchannel.As the reaction goes on,the mass fraction of CO2increases with different rates.On one side,for the chemical kinetics mechanism of Mansouriet al.[27],there is little CO2in the first half position of the microchannel,and CO2is produced in the second half of the microchannel.The reason is that C3H8is oxidized to be CO at first,and then CO is oxidized to CO2finally.For Konakovet al.[28],the mass fraction of CO2tends to remain a high value at the distance not far from the entrance of the microchannel.On the other side,for the chemical kinetics mechanisms from Luoet al.[29]and Jones and Lindstedt[23],the mass fraction of CO2increases first and then declines since the CO2decomposes to produce CO and O2.In particular,the mass fraction of CO2predicted by the chemical kinetics mechanism from Jones and Lindstedt[23]is much smaller than the other three,because CO is the main product and CO2is the by-product.

    3.4.2.4.Product of H2O.The mass fractions for product of H2O predicted by different chemical kinetics mechanisms are illustrated in Fig.11.

    As presented in Fig.11,there is no H2O at the inlet of the microchannel.As the reaction goes on,the mass fraction of H2O increases noticeably.Finally,mass fraction of H2O tends to keep unchanged.

    The fundamental data depicted in Figs.8-11 are valuable since they are too difficult to get in the experiments due to the tiny dimensions of the microchannel.As a result,the understanding of propane-air combustion in the microchannel is deepened in this work.

    Fig.10.Mass fraction of CO2 for the physical structure A:(a)Jones and Lindstedt[23];(b)Konakov et al.[28];(c)Luo et al.[29];(d)Mansouri et al.[27].

    3.4.2.5.Influence of chemical kinetics mechanisms on the final position of reaction.Distributions of species mass fraction in the centerline of physical structure A predicted by four chemical kinetics mechanisms are plotted in Fig.12.For Mansouriet al.[27]and Konakovet al.[28],CO2is the main product.On the contrast,CO is the main product predicted by the chemical kinetics mechanism from Jones and Lindstedt[23].Apparently,for the chemical kinetics mechanism from Luoet al.[29],the final mass fraction of CO and CO2is approximately equal.Additionally,species mass fractions of the other products are different.Besides,if a chemical reaction is in progress in the microchannel,the concentrations of the reactants will reduce and the product concentration will increase[30].When the mass fractions of species reduce to be a constant,the reaction is complete.From Fig.12,we can infer that different chemical kinetics mechanisms bring about different final position of reaction.The final positions of reaction are 4 mm,4 mm,8 mm and 10 mm for the mechanisms from Jones and Lindstedt[23],Konakovet al.[28],Luoet al.[29]and Mansouriet al.[27],respectively.

    3.5.Results for the physical structure B

    Because the chemical kinetics mechanism of Mansouriet al.[27]among the four chemical kinetics models is the most accurate one to predict the flame propagation in the microchannels,the temperature and species mass fraction profiles at the center surface of the physical structure B with the chemical kinetics mechanism of Mansouriet al.[27]are revealed in Fig.13.The temperature rises first since propane combustion releases heat,and then declines enormously by the wall cooling effect.The propane is oxidized to the intermediate product of CO in the first half length,and finally CO is oxidized to produce the final product of CO2in the second half.

    It can be inferred that the prediction with the integrated mathematical model established in this work agrees well with the experimental results qualitatively and quantitatively and the verified mathematical model can be applied for further design of the crimped ribbon flame arresters.

    4.Conclusions

    In this work,a reliable and accurate mathematical model including steady fluid flow at elevated temperatures,violent combustion reaction and heat transfer was set up.The turbulent model was used to predict the flame propagation in the microchannels of the crimped ribbon flame arresters under the deflagration and detonation conditions.Four chemical kinetics mechanisms on the prediction accuracy of the flame characteristics in the microchannel were investigated.

    The following conclusions can be obtained in this work.First,the pattern of fluid flow in the microchannel of the crimped ribbon flame arresters with flame propagation is con firmed to be turbulent under the deflagration and detonation conditions.Second,the auto-ignition temperature of flammable gas instead of a wide range of temperatures could be taken as the judging criterion to determine whether the flame is quenched or not.Third,chemical kinetics mechanisms have a significant impact on the accuracy of prediction,and the kinetics mechanism of Mansouriet al.[27]among the four ones is the most accurate.The temperature and species mass fraction distributions in microchannels,which are too difficult to measure by the experimental method due to the tiny dimension of microchannels,are revealed.In other words,the fundamental insights into chemical reactions and heat loss are well portrayed and the effect of the cold wall for quenching the deflagration and detonation flames in the microchannels is also determined.Fourth,the final position of reaction in the microchannel could be predicted with the integrated model,but different chemical kinetics mechanisms give different final positions of reaction in the microchannel.

    Fig.11.Mass fraction of H2O for the physical structure A:(a)Jones and Lindstedt[23];(b)Konakov et al.[28];(c)Luo et al.[29];(d)Mansouri et al.[27].

    Fig.12.Species mass fraction profiles with different mechanisms for the physical structure A:(a)Jones and Lindstedt[23];(b)Konakov et al.[28];(c)Luo et al.[29];(d)Mansouri et al.[27].

    Fig.13.Temperature and species mass fraction profiles predicted by the chemical kinetics mechanisms of Mansouri et al.[27]for the physical structure B:(a)temperature;(b)C3H8;(c)O2;(d)CO;(e)CO2;(f)H2O.

    The obtained results demonstrate the capability of numerical approach in predicting the behavior of real accidental explosions.The integrated mathematical model verified in this work can be used as an effective tool for modeling,designing and choosing such type of crimped ribbon flame arresters with the propane-air medium in the future.

    Nomenclature

    Arpre-exponential factor,consistent units

    athe bottom side length of the isosceles triangle,mm

    Cj,rmolar concentration of speciesjin reactionr,kmol·m-3

    Cp,jthe specific heat capacity of speciesj,J·(kg·K)-1

    Dthe diameter of the arrester element,mm

    dthe diameter of the pipeline,mm

    Eractivation energy for the reaction,J·kmol-1

    hthe height of the isosceles triangle,mm

    hsensible enthalpy,h=∑jYjhj,J

    hjthe sensible enthalpy of speciesj,J

    Jjthe diffusion flux of speciesj,kg·(m-2·s)-1

    kturbulent kinetic energy,m2·s3

    kf,rforward rate constant for reactionr,kg·(m-3·s)-1

    Mw,ithe molecular weight of speciesi

    Mw,jthe molecular weight of speciesj

    Mw,Rthe molecular weight of reactantR

    NRNumber of reactions

    pthe static pressure,Pa

    Runiversal gas constant,J·(kmol-1·K)-1

    Rithe net source of chemical speciesi,consistent units

    Ri,rthe Arrhenius molarrate of creation/destruction of speciesiin reactionr,consistent units

    Shthe heat of chemical reaction,J

    uinvelocity at the entrance of the microchannel,m·s-1

    uoverall velocity vector,m·s-1

    Vvelocity at the pipeline,m·s-1

    Vfflameproof velocity,m·s-1

    Yjthe mass fraction of speciesj

    YPthe mass fraction of any product speciesP

    YRthe mass fraction of a particular reactantR

    α the deflecting angle between the narrow channel and axial direction,(°)

    β expansion ratio of flame arresters,dimensionless

    γ the basic angle of the isosceles triangle,(°)

    ε turbulent dissipation rate,m2·s3

    η′j,rrate exponent for reactant speciesjin reactionr

    η″j,rrate exponent for product speciesjin reactionr

    ν′i,rstoichiometric coefficient for reactantiin reactionr

    v″i,rstoichiometric coefficient for productiin reactionr

    v″j,rstoichiometric coefficient for productjin reactionr

    ν′R,rstoichiometric coefficient for reactantRin reactionr

    ρ density,kg·m-3

    the stress tensor,Pa

    φ general variable

    Subscripts

    eff effective

    f flameproof

    H hydraulic

    in inlet

    ref. reference

    tturbulent

    Acknowledgments

    Our gratitude goes to the Supercomputing Center of USTC(University of Science and Technology of China).

    [1]D.Lietze,Crimped metal ribbon flame arrestors for the protection of gas measurement systems,J.Loss Prev.Process Ind.15(1)(2002)29-35.

    [2]A.A.Pekalski,J.F.Zevenbergen,S.M.Lemkowitz,H.J.Pasman,A review of explosion prevention and protection systems suitable as ultimate layer of protection in chemical process installations,Process Saf.Environ.Prot.83(1)(2005)1-17.

    [3]P.Bauer,Experimental investigation on flame and detonation quenching:applicability of static flame arresters,J.Loss Prev.Process Ind.18(2)(2005)63-68.

    [4]C.Kersten,H.F?rster,Investigation of deflagrations and detonations in pipes and flame arresters by high-speed framing,J.Loss Prev.Process Ind.17(1)(2004)43-50.

    [5]Y.Okawa,C.Youn,T.Kagawa,A study of the characteristics of flow rate and extinction in a flame arrester with radial slit structure,J.Loss Prev.Process Ind.25(2)(2012)242-249.

    [6]L.Wang,H.Ma,Z.Shen,The quenching of propane deflagrations by crimped ribbon flame arrestors,J.Loss Prev.Process Ind.43(2016)567-574.

    [7]A.C.Benim,S.Iqbal,W.Meier,F.Joos,A.Wiedermann,Numerical investigation of turbulent swirling flames with validation in a gas turbine model combustor,Appl.Therm.Eng.110(2017)202-212.

    [8]M.Choi,Y.Sung,M.Won,Y.Park,M.Kim,G.Choi,D.Kim,Effect of fuel distribution on turbulence and combustion characteristics of a micro gas turbine combustor,J.Ind.Eng.Chem.48(2017)24-35.

    [9]J.Wan,W.Yang,A.Fan,Y.Liu,H.Yao,W.Liu,Y.Du,D.Zhao,A numerical investigation on combustion characteristics of H2/air mixture in a micro-combustor with wall cavities,Int.J.Hydrog.Energy39(15)(2014)8138-8146.

    [10]V.Giovannoni,R.N.Sharma,R.R.Raine,Premixed combustion of methane-air mixture stabilized over porous medium:A 2D numerical study,Chem.Eng.Sci.152(2016)591-605.

    [11]X.Li,J.Zhang,H.Yang,L.Jiang,X.Wang,D.Zhao,Combustion characteristics of nonpremixed methane micro-jet flame in co flow air and thermal interaction between flame and micro tube,Appl.Therm.Eng.112(2017)296-303.

    [12]S.Raimondeau,D.Norton,D.G.Vlachos,R.I.Masel,Modeling of high-temperature microburners,Proc.Combust.Inst.29(1)(2002)901-907.

    [13]M.Yu,K.Zheng,T.Chu,Gas explosion flame propagation over various hollowsquare obstacles,J.Nat.Gas Sci.Eng.30(2016)221-227.

    [14]J.Miao,C.W.Leung,C.S.Cheung,Z.Huang,W.Jin,Effect of H2addition on OH distribution of LPG/Air circumferential inverse diffusion flame,Int.J.Hydrog.Energy41(22)(2016)9653-9663.

    [15]C.M.R.Vendra,J.X.Wen,V.H.Y.Tam,Numerical simulation of turbulent flame-wall quenching using a coherent flame model,J.Loss Prev.Process Ind.26(2)(2013)363-368.

    [16]M.Bin Shams,E.M.Elkanzi,Z.Ramadhan,S.Rahma,M.Khamis,Gas turbine inlet air cooling system for enhancing propane recovery in a gas plant:Theoretical and cost analyses,J.Nat.Gas Sci.Eng.43(2017)22-32.

    [17]I.Dincer,C.Zamfirescu,A review of novel energy options for clean rail applications,J.Nat.Gas Sci.Eng.28(2016)461-478.

    [18]R.K.Abdrakhmanov,B.F.Boyarshinov,S.Y.Fedorov,Investigation of the local parameters of a cellular propane/butane/air flame,Int.J.Heat Mass Transf.109(2017)1172-1180.

    [19]B.Mohammadreza,T.Sadegh,F.L.Morteza,Experimental study on the effects of mixture flow rate,equivalence ratio,oxygen enhancement,and geometrical parameters on propane-air premixed flame dynamics in non-adiabatic meso-scale reactors,Energy121(2017)657-675.

    [20]M.Cathonnet,Chemical kinetic modeling of combustion from 1969 to 2019,Combust.Sci.Technol.98(4-6)(1994)265-279.

    [21]J.Simmie,Detailed chemical kinetic models for the combustion of hydrocarbon fuels,Prog.Energy Combust.Sci.29(6)(2003)599-634.

    [22]M.Cord,B.Husson,H.Lizardo,C.Juan,O.Herbinet,P.A.Glaude,R.Fournet,B.Sirjean,B.L.Frédérique,R.L.Manuel,Z.Wang,Study of the low temperature oxidation of propane,J.Phys.Chem.A116(50)(2012)12214-12,228.

    [23]W.Jones,R.Lindstedt,Global reaction schemes for hydrocarbon combustion,Combust.Flame73(3)(1988)233-249.

    [24]Y.Jiang,R.Qiu,Areduced mechanism for flame inhibition by phosphorus-containing compounds based on level of importance analysis,Chin.J.Chem.Eng.18(5)(2010)711-720.

    [25]M.Bahmani,J.Shariati,A.N.Rouzbahani,Simulation and optimization of an industrial gas condensate stabilization unit to modify LPG and NGL production with minimizing CO2emission to the environment,Chin.J.Chem.Eng.25(3)(2017)338-346.

    [26]N.Zhang,T.Qiu,B.Chen,CFD simulation of propane cracking tube using detailed radical kinetic mechanism,Chin.J.Chem.Eng.21(12)(2013)1319-1331.

    [27]Z.Mansouri,M.Aouissi,T.Boushaki,Numerical computations of premixed propane flame in a swirl-stabilized burner:Effects of hydrogen enrichment,swirl number and equivalence ratio on flame characteristics,Int.J.Hydrog.Energy41(22)(2016)9664-9678.

    [28]S.A.Konakov,S.V.Dzyubanenko,V.V.Krzhizhanovskaya,Computer simulation ?approach in development of propane-air combustor microreactor,Procedia Comput.Sci.101(2016)76-85.

    [29]C.Luo,J.Zanganeh,B.Moghtaderi,A 3D numerical study on the effects of obstacles on flame propagation in a cylindrical explosion vessel connected to a vented tube,J.Loss Prev.Process Ind.44(2016)53-61.

    [30]A.Gutkowski,Numerical analysis of effect of ignition methods on flame behavior during passing through a sudden contraction near the quenching conditions,Appl.Therm.Eng.54(1)(2013)202-211.

    [31]Q.Huang,W.Zhang,C.Yang,Modeling transport phenomena and reactions in a pilot slurry airlift loop reactor for direct coal liquefaction,Chem.Eng.Sci.135(2015)441-451.

    [32]B.Lauder,D.Spalding,Lectures in Mathematical Models of Turbulence,Academic Press,London,UK,1972.

    [33]B.F.Magnussen,B.H.Hjertager,On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion,Symp.Combust.16(1)(1977)719-729.

    [34]D.Spalding,Mixing and chemical reaction in steady con fined turbulent flames,Symp.Combust.13(1)(1971)649-657.

    [35]Q.Huang,T.Liu,J.Yang,L.Yao,L.Gao,Evaluation of radiative transfer using the finite volume method in cylindrical photoreactors,Chem.Eng.Sci.66(17)(2011)3930-3940.

    [36]Q.Huang,L.Yao,T.Liu,J.Yang,Simulation of the light evolution in an annular photobioreactor for the cultivation ofPorphyridium cruentum,Chem.Eng.Sci.84(2012)718-726.

    [37]H.K.Versteeg,W.Malalasekera,An Introcuction to Computional Fluid Dynamics:The Finite Volume Method,Wiley,New York,1995.

    [38]Q.Huang,C.Yang,G.Z.Yu,Z.Mao,3-D simulations of an internal airlift loop reactor using a steady two-fluid model,Chem.Eng.Technol.30(7)(2007)870-879.

    [39]Q.Huang,C.Yang,G.Yu,Z.S.Mao,Sensitivity study on modeling an internal airlift loop reactor using a steady 2D two-fluid model,Chem.Eng.Technol.31(12)(2008)1790-1798.

    [40]Q.Huang,C.Yang,G.Yu,Z.Mao,CFD simulation of hydrodynamics and mass transfer in an internal airlift loop reactor using a steady two-fluid model,Chem.Eng.Sci.65(20)(2010)5527-5536.

    [41]K.Zhou,Z.Li,Z.Zhou,The quenching of deflagration by crimped ribbon flame arresters,J.Chin.Univ.Sci.Technol.27(4)(1997)449-454.

    猜你喜歡
    創(chuàng)城市容市政府
    市政府召開常務(wù)會議
    市政府召開常務(wù)會議
    開發(fā)建設(shè)導(dǎo)則于城市開發(fā)建設(shè)實踐——以桃浦智創(chuàng)城為例
    市政府召開常務(wù)會議
    “創(chuàng)城”進行時
    少先隊活動(2021年6期)2021-07-22 08:44:16
    市政府召開常務(wù)會議
    一片狼藉
    心系老區(qū)環(huán)衛(wèi) 扮靚水城市容——記沁縣環(huán)衛(wèi)中心黨支部書記、主任王建軍
    我家與“創(chuàng)城”
    我家與“創(chuàng)城”
    成人亚洲精品一区在线观看| 亚洲狠狠婷婷综合久久图片| 久久天躁狠狠躁夜夜2o2o| av网站免费在线观看视频| 欧美日韩福利视频一区二区| 国内久久婷婷六月综合欲色啪| 欧美乱码精品一区二区三区| 亚洲va日本ⅴa欧美va伊人久久| 中文字幕精品免费在线观看视频| 在线观看舔阴道视频| 在线观看舔阴道视频| 国内精品久久久久精免费| 国产熟女午夜一区二区三区| 夜夜爽天天搞| 制服丝袜大香蕉在线| 久久人妻福利社区极品人妻图片| 国产精品久久视频播放| 欧美国产日韩亚洲一区| 麻豆成人av在线观看| 12—13女人毛片做爰片一| 久久天堂一区二区三区四区| 男人的好看免费观看在线视频 | 1024视频免费在线观看| 国产精品一区二区三区四区久久 | 久久中文字幕人妻熟女| 亚洲成人国产一区在线观看| 99精品欧美一区二区三区四区| 国产精品久久久久久人妻精品电影| 美女 人体艺术 gogo| 人人澡人人妻人| 精品不卡国产一区二区三区| 亚洲,欧美精品.| 亚洲五月天丁香| 人妻久久中文字幕网| 熟妇人妻久久中文字幕3abv| 正在播放国产对白刺激| 亚洲欧美激情在线| 欧美在线黄色| 电影成人av| x7x7x7水蜜桃| 老司机深夜福利视频在线观看| 日韩欧美国产一区二区入口| 视频区欧美日本亚洲| 18禁美女被吸乳视频| 狠狠狠狠99中文字幕| 国产亚洲欧美在线一区二区| 久久久久久国产a免费观看| 亚洲狠狠婷婷综合久久图片| av视频在线观看入口| 欧美日韩中文字幕国产精品一区二区三区 | 久久人人97超碰香蕉20202| 亚洲av成人不卡在线观看播放网| 巨乳人妻的诱惑在线观看| 久久狼人影院| 91精品三级在线观看| 91av网站免费观看| 又大又爽又粗| xxx96com| 大型黄色视频在线免费观看| 女生性感内裤真人,穿戴方法视频| 午夜精品在线福利| 国产激情欧美一区二区| 亚洲少妇的诱惑av| 51午夜福利影视在线观看| 中文字幕色久视频| 免费在线观看视频国产中文字幕亚洲| 一级a爱视频在线免费观看| 女人被躁到高潮嗷嗷叫费观| 天天躁夜夜躁狠狠躁躁| 亚洲电影在线观看av| 亚洲一区二区三区色噜噜| 日本黄色视频三级网站网址| 老汉色∧v一级毛片| 亚洲精品中文字幕在线视频| 伊人久久大香线蕉亚洲五| 国产xxxxx性猛交| 国内毛片毛片毛片毛片毛片| 十分钟在线观看高清视频www| 精品国产国语对白av| av天堂在线播放| 亚洲狠狠婷婷综合久久图片| 美女扒开内裤让男人捅视频| 夜夜看夜夜爽夜夜摸| 人人妻人人爽人人添夜夜欢视频| 人人妻人人澡人人看| 亚洲中文日韩欧美视频| 天天躁狠狠躁夜夜躁狠狠躁| 日本a在线网址| 精品欧美一区二区三区在线| 亚洲色图av天堂| 人人妻人人澡欧美一区二区 | а√天堂www在线а√下载| 欧美成狂野欧美在线观看| 精品国产一区二区久久| 久久精品91无色码中文字幕| 免费搜索国产男女视频| 久久天躁狠狠躁夜夜2o2o| 黄片大片在线免费观看| 老司机午夜十八禁免费视频| 欧美另类亚洲清纯唯美| 脱女人内裤的视频| 久久人人97超碰香蕉20202| 90打野战视频偷拍视频| 亚洲一区高清亚洲精品| 精品不卡国产一区二区三区| 日本精品一区二区三区蜜桃| 日本撒尿小便嘘嘘汇集6| 女人爽到高潮嗷嗷叫在线视频| 村上凉子中文字幕在线| 男女下面插进去视频免费观看| 美女扒开内裤让男人捅视频| 女性被躁到高潮视频| 国产一区二区三区视频了| 久久久国产欧美日韩av| 99久久综合精品五月天人人| 免费高清视频大片| 免费女性裸体啪啪无遮挡网站| 十分钟在线观看高清视频www| 国产免费男女视频| 久久狼人影院| 如日韩欧美国产精品一区二区三区| 我的亚洲天堂| 男女午夜视频在线观看| 少妇裸体淫交视频免费看高清 | 丰满人妻熟妇乱又伦精品不卡| 日韩欧美免费精品| 一级片免费观看大全| 99香蕉大伊视频| 国产精品电影一区二区三区| 久久人人精品亚洲av| 久久人妻熟女aⅴ| 欧洲精品卡2卡3卡4卡5卡区| 亚洲中文字幕日韩| 女生性感内裤真人,穿戴方法视频| 美女 人体艺术 gogo| 国产免费男女视频| 中文字幕人妻丝袜一区二区| 久久久久久亚洲精品国产蜜桃av| 少妇熟女aⅴ在线视频| 中文字幕人成人乱码亚洲影| 成年人黄色毛片网站| 中文字幕av电影在线播放| 国产av在哪里看| 久久人妻熟女aⅴ| 国产成人啪精品午夜网站| 欧美成人午夜精品| 亚洲一卡2卡3卡4卡5卡精品中文| 亚洲午夜理论影院| 国产精品久久久人人做人人爽| 给我免费播放毛片高清在线观看| netflix在线观看网站| 国产精品久久久av美女十八| 久久午夜亚洲精品久久| 99久久99久久久精品蜜桃| 国产又爽黄色视频| 欧美乱妇无乱码| 丝袜美足系列| 在线永久观看黄色视频| 别揉我奶头~嗯~啊~动态视频| 91老司机精品| 涩涩av久久男人的天堂| 一本久久中文字幕| 色综合婷婷激情| 亚洲国产看品久久| 操出白浆在线播放| 欧美不卡视频在线免费观看 | 97人妻精品一区二区三区麻豆 | 国产欧美日韩一区二区精品| 少妇熟女aⅴ在线视频| 91麻豆精品激情在线观看国产| 国产欧美日韩一区二区精品| 美女国产高潮福利片在线看| 国产精品,欧美在线| 在线观看免费午夜福利视频| 亚洲欧美精品综合一区二区三区| 天堂动漫精品| 美女扒开内裤让男人捅视频| www.自偷自拍.com| 午夜福利在线观看吧| 国产99久久九九免费精品| 国产亚洲精品第一综合不卡| 精品欧美国产一区二区三| 午夜日韩欧美国产| av在线播放免费不卡| 男人的好看免费观看在线视频 | 久久精品人人爽人人爽视色| tocl精华| 欧美一级a爱片免费观看看 | 美女扒开内裤让男人捅视频| 欧美不卡视频在线免费观看 | 国产99久久九九免费精品| av免费在线观看网站| 午夜福利18| 亚洲熟妇熟女久久| 一级作爱视频免费观看| 亚洲第一欧美日韩一区二区三区| 法律面前人人平等表现在哪些方面| 欧美黑人欧美精品刺激| 久久久久久亚洲精品国产蜜桃av| 欧美人与性动交α欧美精品济南到| 国产国语露脸激情在线看| 狂野欧美激情性xxxx| 久久久久久亚洲精品国产蜜桃av| 国产精品 欧美亚洲| 黄色毛片三级朝国网站| 国产成+人综合+亚洲专区| 在线观看免费视频日本深夜| 每晚都被弄得嗷嗷叫到高潮| 亚洲男人天堂网一区| 在线观看午夜福利视频| 成人亚洲精品av一区二区| 成人永久免费在线观看视频| 巨乳人妻的诱惑在线观看| 欧美在线一区亚洲| 午夜久久久久精精品| 国产99白浆流出| 中文字幕人妻丝袜一区二区| 男人操女人黄网站| 免费观看精品视频网站| 12—13女人毛片做爰片一| 久久精品人人爽人人爽视色| 国产成人欧美在线观看| 岛国在线观看网站| 亚洲av片天天在线观看| 涩涩av久久男人的天堂| 人成视频在线观看免费观看| 免费久久久久久久精品成人欧美视频| 不卡av一区二区三区| 亚洲中文字幕一区二区三区有码在线看 | 精品一区二区三区av网在线观看| 少妇熟女aⅴ在线视频| 日韩精品青青久久久久久| 在线观看免费视频网站a站| 中文亚洲av片在线观看爽| 丁香欧美五月| 性色av乱码一区二区三区2| 精品国产一区二区久久| 免费在线观看视频国产中文字幕亚洲| 日本欧美视频一区| а√天堂www在线а√下载| 欧美日韩一级在线毛片| 国产精品久久电影中文字幕| 久久精品亚洲熟妇少妇任你| 村上凉子中文字幕在线| 免费看美女性在线毛片视频| 桃色一区二区三区在线观看| 看片在线看免费视频| 每晚都被弄得嗷嗷叫到高潮| 人人妻人人爽人人添夜夜欢视频| 欧美另类亚洲清纯唯美| 亚洲全国av大片| 色精品久久人妻99蜜桃| 免费在线观看视频国产中文字幕亚洲| 久久久久久国产a免费观看| 一卡2卡三卡四卡精品乱码亚洲| 不卡一级毛片| 免费看十八禁软件| 精品卡一卡二卡四卡免费| 1024香蕉在线观看| www.精华液| 亚洲伊人色综图| 高清毛片免费观看视频网站| 欧美不卡视频在线免费观看 | 国产高清视频在线播放一区| 18禁裸乳无遮挡免费网站照片 | 免费无遮挡裸体视频| 女性生殖器流出的白浆| 激情视频va一区二区三区| 久久人妻福利社区极品人妻图片| 久久国产亚洲av麻豆专区| 久久狼人影院| 精品一区二区三区四区五区乱码| 亚洲熟女毛片儿| 中文字幕另类日韩欧美亚洲嫩草| 久久久久久久久中文| 亚洲国产看品久久| 99riav亚洲国产免费| 久久国产精品人妻蜜桃| 757午夜福利合集在线观看| 中文亚洲av片在线观看爽| 亚洲精品在线观看二区| 男女之事视频高清在线观看| svipshipincom国产片| 一本大道久久a久久精品| 亚洲国产毛片av蜜桃av| a级毛片在线看网站| 亚洲第一电影网av| 午夜福利高清视频| 久久人妻福利社区极品人妻图片| 国产av一区在线观看免费| 亚洲人成网站在线播放欧美日韩| 国产精品免费一区二区三区在线| 热99re8久久精品国产| 亚洲国产欧美网| 村上凉子中文字幕在线| 丝袜人妻中文字幕| 国产精品一区二区精品视频观看| 大型av网站在线播放| 91在线观看av| 久久久久久久午夜电影| 国产精品野战在线观看| 亚洲av成人av| 国产在线精品亚洲第一网站| 男女做爰动态图高潮gif福利片 | 精品人妻1区二区| 真人一进一出gif抽搐免费| 亚洲专区国产一区二区| av超薄肉色丝袜交足视频| 国产伦一二天堂av在线观看| 俄罗斯特黄特色一大片| 国产亚洲精品综合一区在线观看 | 999久久久精品免费观看国产| 九色亚洲精品在线播放| 亚洲最大成人中文| 日韩中文字幕欧美一区二区| 免费av毛片视频| 亚洲五月天丁香| 女人精品久久久久毛片| 在线观看舔阴道视频| 老司机靠b影院| 亚洲熟妇熟女久久| av超薄肉色丝袜交足视频| 久久午夜综合久久蜜桃| 午夜福利高清视频| 免费高清在线观看日韩| av片东京热男人的天堂| 岛国在线观看网站| 午夜福利一区二区在线看| 国产欧美日韩一区二区精品| 香蕉久久夜色| 国产主播在线观看一区二区| av网站免费在线观看视频| svipshipincom国产片| 欧美乱码精品一区二区三区| www.自偷自拍.com| 久久香蕉国产精品| 操出白浆在线播放| 好看av亚洲va欧美ⅴa在| 97碰自拍视频| 黄色视频不卡| 国产麻豆成人av免费视频| 9191精品国产免费久久| 99久久久亚洲精品蜜臀av| 色综合站精品国产| 免费在线观看亚洲国产| 欧美日韩瑟瑟在线播放| 欧美黑人欧美精品刺激| 免费不卡黄色视频| 少妇熟女aⅴ在线视频| 亚洲欧美日韩无卡精品| 18美女黄网站色大片免费观看| 人人澡人人妻人| 不卡一级毛片| 欧美日韩中文字幕国产精品一区二区三区 | 自线自在国产av| 亚洲熟妇中文字幕五十中出| 亚洲精品久久国产高清桃花| 午夜久久久在线观看| 在线观看www视频免费| 在线天堂中文资源库| 少妇的丰满在线观看| 夜夜看夜夜爽夜夜摸| 叶爱在线成人免费视频播放| 亚洲国产中文字幕在线视频| 又紧又爽又黄一区二区| 一进一出好大好爽视频| 免费看十八禁软件| 午夜免费鲁丝| 午夜福利高清视频| 欧美激情高清一区二区三区| 久久影院123| 侵犯人妻中文字幕一二三四区| 日韩欧美三级三区| 久99久视频精品免费| 亚洲va日本ⅴa欧美va伊人久久| 啦啦啦免费观看视频1| 国产色视频综合| 午夜日韩欧美国产| 亚洲片人在线观看| 可以在线观看的亚洲视频| 激情在线观看视频在线高清| 99热只有精品国产| 三级毛片av免费| 无限看片的www在线观看| 国产av在哪里看| 香蕉丝袜av| 亚洲av电影不卡..在线观看| 欧美成人午夜精品| 91精品三级在线观看| e午夜精品久久久久久久| 国产乱人伦免费视频| 91麻豆精品激情在线观看国产| 麻豆国产av国片精品| 美国免费a级毛片| 免费av毛片视频| 午夜老司机福利片| 久热爱精品视频在线9| 一个人免费在线观看的高清视频| 丰满人妻熟妇乱又伦精品不卡| 欧美亚洲日本最大视频资源| 琪琪午夜伦伦电影理论片6080| av福利片在线| 一本大道久久a久久精品| 成人三级黄色视频| 亚洲av第一区精品v没综合| svipshipincom国产片| 老司机福利观看| 一边摸一边抽搐一进一出视频| 久久中文看片网| 首页视频小说图片口味搜索| 亚洲中文字幕日韩| 日本a在线网址| 啦啦啦观看免费观看视频高清 | 男人操女人黄网站| 热re99久久国产66热| 国产又色又爽无遮挡免费看| 日韩大码丰满熟妇| 在线国产一区二区在线| 久久狼人影院| 在线观看免费午夜福利视频| 老司机在亚洲福利影院| 岛国视频午夜一区免费看| 成人亚洲精品av一区二区| 亚洲成人久久性| 色播亚洲综合网| 黑人欧美特级aaaaaa片| 韩国精品一区二区三区| 亚洲人成电影观看| 亚洲中文日韩欧美视频| 日韩av在线大香蕉| 午夜福利影视在线免费观看| 久久中文字幕一级| 精品高清国产在线一区| 亚洲久久久国产精品| 久久精品国产99精品国产亚洲性色 | 欧洲精品卡2卡3卡4卡5卡区| 妹子高潮喷水视频| 黄频高清免费视频| 久久午夜综合久久蜜桃| 亚洲欧美日韩无卡精品| 给我免费播放毛片高清在线观看| 亚洲最大成人中文| 国产一区二区三区综合在线观看| 日韩欧美在线二视频| 国产亚洲精品第一综合不卡| 午夜精品久久久久久毛片777| 侵犯人妻中文字幕一二三四区| 精品国产乱码久久久久久男人| 亚洲 欧美一区二区三区| 欧美午夜高清在线| 大码成人一级视频| 神马国产精品三级电影在线观看 | 啦啦啦韩国在线观看视频| 亚洲av成人不卡在线观看播放网| 精品免费久久久久久久清纯| 亚洲中文日韩欧美视频| 十八禁网站免费在线| 非洲黑人性xxxx精品又粗又长| 亚洲欧美激情在线| 久久人妻熟女aⅴ| 看片在线看免费视频| 美女高潮喷水抽搐中文字幕| 日本三级黄在线观看| 在线视频色国产色| 欧美成狂野欧美在线观看| 亚洲专区国产一区二区| 久久人妻福利社区极品人妻图片| 日日夜夜操网爽| 99热只有精品国产| 国产亚洲欧美精品永久| 一区二区三区高清视频在线| 黑丝袜美女国产一区| 夜夜夜夜夜久久久久| 在线观看66精品国产| 亚洲性夜色夜夜综合| 精品第一国产精品| 免费看a级黄色片| 久9热在线精品视频| 久热爱精品视频在线9| 国产伦人伦偷精品视频| 久久久久久久久中文| 日韩三级视频一区二区三区| xxx96com| 国产精品久久久久久精品电影 | 久久久久国产一级毛片高清牌| 国产精品影院久久| 欧美日韩一级在线毛片| 久久久久久免费高清国产稀缺| 免费搜索国产男女视频| 国产1区2区3区精品| 男女下面进入的视频免费午夜 | 免费在线观看完整版高清| 国产精品亚洲av一区麻豆| 黄网站色视频无遮挡免费观看| 亚洲国产欧美一区二区综合| 国产免费男女视频| 天天添夜夜摸| 国产成人精品久久二区二区免费| 亚洲av第一区精品v没综合| 久久香蕉精品热| 午夜老司机福利片| 此物有八面人人有两片| 国产麻豆成人av免费视频| 欧美一区二区精品小视频在线| 亚洲精华国产精华精| 亚洲午夜理论影院| 欧美黄色片欧美黄色片| 首页视频小说图片口味搜索| 色哟哟哟哟哟哟| 日本 欧美在线| 女性被躁到高潮视频| 99国产综合亚洲精品| 精品国产乱子伦一区二区三区| 国产精品久久久人人做人人爽| 91麻豆av在线| 久久国产精品男人的天堂亚洲| 可以免费在线观看a视频的电影网站| 天天躁狠狠躁夜夜躁狠狠躁| 黄色片一级片一级黄色片| 欧美一区二区精品小视频在线| 禁无遮挡网站| 在线十欧美十亚洲十日本专区| 19禁男女啪啪无遮挡网站| 欧美大码av| 国产在线观看jvid| 美女 人体艺术 gogo| 12—13女人毛片做爰片一| 青草久久国产| 国产亚洲精品久久久久久毛片| 夜夜爽天天搞| 亚洲一区中文字幕在线| 一个人免费在线观看的高清视频| 国产成人免费无遮挡视频| 日韩三级视频一区二区三区| 国产精品野战在线观看| 久久久久久久午夜电影| xxx96com| 欧美激情 高清一区二区三区| 日本三级黄在线观看| 无遮挡黄片免费观看| 亚洲成a人片在线一区二区| 成人国产一区最新在线观看| 午夜福利视频1000在线观看 | 99在线人妻在线中文字幕| 欧美一级a爱片免费观看看 | 久久国产亚洲av麻豆专区| 长腿黑丝高跟| 欧美乱码精品一区二区三区| 丝袜人妻中文字幕| 午夜免费成人在线视频| 亚洲天堂国产精品一区在线| 叶爱在线成人免费视频播放| 亚洲电影在线观看av| 9热在线视频观看99| 久久精品成人免费网站| 神马国产精品三级电影在线观看 | 亚洲一码二码三码区别大吗| 九色亚洲精品在线播放| 麻豆av在线久日| 在线观看66精品国产| 18美女黄网站色大片免费观看| 国产乱人伦免费视频| 国产伦人伦偷精品视频| 99国产精品一区二区三区| 亚洲九九香蕉| 咕卡用的链子| 91成年电影在线观看| 啦啦啦韩国在线观看视频| 美女高潮喷水抽搐中文字幕| 夜夜躁狠狠躁天天躁| АⅤ资源中文在线天堂| 日日夜夜操网爽| 亚洲av美国av| 欧美一级a爱片免费观看看 | 精品一区二区三区av网在线观看| 日本免费a在线| 999精品在线视频| 亚洲欧美精品综合久久99| 国产一区二区激情短视频| 久热爱精品视频在线9| 操美女的视频在线观看| 男女床上黄色一级片免费看| 女性被躁到高潮视频| 黄色片一级片一级黄色片| 欧美午夜高清在线| 欧美绝顶高潮抽搐喷水| 99精品久久久久人妻精品| 日韩欧美一区视频在线观看| 黑人操中国人逼视频| 成年人黄色毛片网站| 国产又色又爽无遮挡免费看| 中文字幕高清在线视频| 一区二区三区精品91| 啦啦啦免费观看视频1| 91精品国产国语对白视频| av网站免费在线观看视频| 神马国产精品三级电影在线观看 | 人妻久久中文字幕网| 视频区欧美日本亚洲| 亚洲中文字幕日韩| 精品日产1卡2卡| 久久精品影院6| 黄色成人免费大全| 久久久国产成人精品二区| 啦啦啦 在线观看视频| 一级黄色大片毛片| 中文字幕色久视频| 琪琪午夜伦伦电影理论片6080| 国产精品亚洲一级av第二区| 久久久久国产一级毛片高清牌| 琪琪午夜伦伦电影理论片6080| 一级黄色大片毛片| 一进一出抽搐动态| 国产在线观看jvid| 国产aⅴ精品一区二区三区波|