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

    Effect of non-condensable gas on a collapsing cavitation bubble near solid wall investigated by multicomponent thermal MRT-LBM?

    2021-03-11 08:32:58YuYang楊雨MingLeiShan單鳴雷QingBangHan韓慶邦andXueFenKan闞雪芬
    Chinese Physics B 2021年2期
    關(guān)鍵詞:韓慶

    Yu Yang(楊雨), Ming-Lei Shan(單鳴雷), Qing-Bang Han(韓慶邦), and Xue-Fen Kan(闞雪芬)

    1Jiangsu Key Laboratory of Power Transmission and Distribution Equipment Technology,Hohai University,Changzhou 213022,China

    2College of Computer and Information,Hohai University,Nanjing 210000,China

    Keywords: multicomponent,cavitation bubble,non-condensable gas,lattice Boltzmann method

    1. Introduction

    Cavitation is ubiquitous in the natural and industrial fields. A cavitation bubble collapsing near a solid wall will lead the extreme conditions to form locally,such as high pressure, temperature and velocity.[1]Owing to these cavitation properties,research on the potential of cavitation exploitation is currently an interesting topic. Some optimized applications of cavitation have been established, such as cavitation degradation,cavitation sterilization,and lab-on-chip.[2–4]However,cavitation treatment is difficult to optimize because the mechanisms of the cavitation phenomenon are not completely illustrated in the presently reviewed literature.

    Many numerical simulation methods have been suggested to investigate the collapsing cavitation bubbles because experiments are difficult to conduct and are deficient in capturing the internal mechanism of cavitation bubbles. Computational fluid dynamics (CFD) methods have been successfully employed for the simulation of multiphase flow.[5]However,they require complex interface-tracking or interface-capturing technology, and interactions among the constituent molecules on a microscopic scale cannot be considered. In contrast to conventional CFD methods,the lattice Boltzmann method(LBM)based on the microscopic nature and mesoscopic characteristics performs well in the numerical simulations of complex gas-liquid phase change. The main advantages, such as clear physical background,perfect parallel computing scheme, linear convection term, and extraordinary capability of dealing with complex boundaries, distinguish the LBM from other numerical methodologies. The Shan–Chen[6]pseudopotential method stands out in comparison with other multiphase lattice Boltzmann (LB) models because of its simple implementation and automatic separation of phases or components without interface-tracking or interface-capturing.

    In the last two decades, many efforts have been made to apply the pseudopotential LBM to investigating the phase change and heat transfer phenomena,such as cavitation,evaporation, and boiling.[7–9]Several models have been established by introducing thermal equations into LB multiphase models.[10–12]Generally, thermal LB models are divided into three categories: the multispeed approach,[13]the doubledistribution-function (DDF) approach,[14,15]and the hybrid approach.[16,17]Among them,the pseudopotential DDF model proposed by Zhang and Chen has been widely used in the simulation of multiphase flow with phase change.[18]Later,Gong and Chen improved the thermal LBM by adding a new form of the source term to the temperature distribution function,thus reducing the spurious currents.[19]Meanwhile,multicomponent (MCP) LB models are being developed rapidly.[20,21]Bao and Schaefer considerably improved the multicomponent multiphase model,but thermal effects have not yet been incorporated into this model.[22]Qin et al. studied the heat transfer of multicomponents by using a single temperature evolution equation.[23]Ikeda et al. proposed a thermal immiscible multicomponent LB model.[24]In their model, each component employs the respective temperature equation, and then the temperature is coupled at each node by using the densityweighted combination. Later, one addressed many aspects of the multicomponent thermal model, such as the large density ratio, spurious current and numerical instability.[25–28]This advances the LBM closer to the realistic physical phenomena in simulating multicomponent multiphase flow.

    Based on single-component DDF LBM, Yang et al.[29]used the improved thermal multiple-relaxation-time (MRT)-LBM model to simulate a collapsing cavitation bubble and discussed the effects of the position offset parameter and pressure difference on the maximum temperature of the bubble. However, cavitation bubble contains vapor and non-condensable gas in real physical phenomena. Some studies indicate that experiments with a small non-condensable gas concentration inside the bubble show its strong influence on the dynamics and thermodynamics of the bubble.[30,31]However, the internal influence of the non-condensable gas on the collapsing cavitation bubble cannot be explained at a microscopic level.Therefore,it is of great significance to study the effect of noncondensable gas inside bubble on the bubble collapsing near a solid wall by using numerical methods. To the best of the author’s knowledge, little researches have been reported on the dynamic and thermodynamic behaviors of multicomponent cavitation bubble by using LBM in the literature.

    In the present work,a developed multicomponent thermal MRT-LBM is presented to numerically study the collapsing cavitation bubble. The rest of this paper is organized as follows. The basis numerical model is presented in Section 2. In Section 3 the model verification and computation domain are described. In Section 4 the effect of the non-condensable gas inside bubble on a collapsing bubble near a wall is illustrated.In Section 5 the conclusions and discussion are presented.

    2. Multicomponent lattice Boltzmann methods

    2.1. Pseudopotential MRT-LBM

    The multicomponent evolution equation of the density distribution function with the MRT collision operator can be given by[32,33]

    and the diagonal matrix is written as

    The macroscopic Naiver–Stokes(NS)equation can be derived from Eq. (1). The mass and momentum equations of each component are as follows:[35,36]

    where ασis the volume fraction satisfying α1+α2=1. The range of α1and α2are both 0–1. According to Eq.(4),we can obtain the continuity equation for the mixture

    For the incompressible or slightly compressible fluids,the viscous stress tensor is

    The average dynamic viscosityμ =α1μ1+α2μ2.

    Using the transformation matrix, equation (1) can be transformed into

    The multicomponent moment space can be obtained from the following expression:

    mσ = M fσ=Mij,σfj,σ=(ρσ,eσ,?σ,jx,σ,qx,σ,jy,σ,qy,σ,pxx,σ,pxy,σ)T, (8)

    where Fσ=(Fx,σ,Fy,σ) refers to the total force of the fluid–fluid interaction acting on the system. Both the intramolecular force and intermolecular force exist in the multicomponent system. The intramolecule force can be expressed as

    where Gσσis the intramolecule interaction strength. In this study,the subscripts σ =1 and σ =2 represent the gas phase and liquid phase,respectively. The gas is regarded as the ideal fluid and thus G11=0. The liquid is the non-ideal fluid. The weight wi=1/3 for |ei|2=1 and the weight wi=1/12 for|ei|2=2 on D2Q9 lattice for the nearest-neighbor interactions.The interaction potential ψ2(x)is adopted to mimic the interaction force of liquid,which is given by[37]

    where a=0.4963(RTc)2/pc, b=0.1873RTc/pc, and Rg=1 is the gas constant,with Tcbeing the critical temperature,and the critical pressure.The constants a=1 and b=4 are used in the present study. For the gas component,the pseudopotential is equal to its density. The intermolecule force is determined by[39]

    The improved forcing scheme proposed by Li et al.[40]is adopted in the present model. An adjustable coefficient ε as shown in Eq. (17) is employed to adjust the pseudopotential LB model to achieve good thermodynamic consistency and improve the mechanical stability of the model. In this work,based on the investigation of the thermodynamic consistency,ε is set to be 1.86 for the liquid component. For the gas component,just set ε =0.

    2.2. Thermal MRT-LBM

    The temperature in an EOS can be derived by solving the target temperature equation[41]

    where α =k/ρcvis the thermal diffusivity, k is the thermal conductivity, ρ is the macroscopic density, cvis the specific heat at a constant volume,and φ is the source term of the temperature.

    In order to solve the target temperature equation, an improved thermal LB model based on the MRT collision operator is used in this model.[42]To simulate the temperature evolution process of gas–liquid two-phase component, the temperature distribution function of the σ-th component is given by

    In numerical implementations, ?tφσ≈ [φσ(t) ?φσ(t ?δt)]/δt.

    Unless otherwise specified,the unit adopted in this paper is the lattice unit of the LBM.The basic units of length,time,and mass are lu (lattice unit), ts (time step), and mu (mass unit),respectively. Thus,the units of density,pressure,velocity, and viscosity are expressed as mu/lu3, mu/(ts2·lu), lu/ts,and lu2/ts,respectively.

    3. Verification of spherical cavitation bubble

    3.1. Validation of Laplace law

    For the cavitation bubble case, Laplace law shows that the pressure difference between both sides of the bubble is inversely proportional to the bubble radius, which can be expressed as

    where pinand poutare the pressure inside and outside the bubble,respectively,R0is the initial bubble radius and θ is the surface tension. An approximately pure gas bubble surrounded by liquid is placed in the center of the periodic domain in a gravity-free field. The density field is initialized into[43]where ρinand ρoutare the density inside and outside the bubble,respectively. (x0,y0)is the center position of the domain.The width of the phase interface W =5. The hyperbolic tangent function tanh(x)=(e2x?1)/(e2x+1). The computational domain is a 201×201 lattice system. Initially, the density of gas(liquid)is 0.003(0.00001)and 0.000001(0.42)inside and outside the bubble, respectively. The initial radius R0is set to be 20, 25, 30, 35, and 40. The results at the different reduced temperatures are shown in Fig.1. The linear relation between the pressure difference ?p and 1/R0can be clearly observed,indicating that the simulation results accord well with those from the Laplace law.

    Fig.1. Plots of ?p versus 1/R0,showing validation of Laplace law.

    3.2. Comparison between simulation results from RP equation and Rayleigh bubble collapse

    The Rayleigh–Plesset (RP) equation can be used to describe the dynamical behavior of cavitation bubble growth and collapse. In order to verify the consistency of the LBM simulation results of the cavitation bubble collapse with the numerical solution of the RP equation, a model of a spherical cavitation bubble with radius R(t)in an infinite liquid field is established as shown in Fig.2. The bubble is assumed to be a mixture of non-condensable gas and vapor. The liquid temperature T∞is considered to be homogeneous. The pressure inside the bubble pb(t)includes the vapor pressure pvand the partial pressure of non-condensable gas pg. The pressure far from the bubble is p∞(t),which is regarded as a known input parameter to control the growth and collapse of the bubble.

    Under the above assumptions,the RP equation is written as[1]

    According to the Meyer equation

    the relationship between γ and cvcan be obtained by combining Eqs.(28)and(29)as follows:

    Fig.2. A spherical cavitation bubble in infinite liquid field.

    The specific heat at constant volume is chosen to be cv=3 in this paper. Thus, the adiabatic coefficient is γ =1.3. For this calculation comparing the numerical and RP equation solutions,the solution of the RP equation is obtained by substituting the corresponding LBM parameters into Eq.(27). The initial densities of the system are the same as those in Subsection 3.1,and the surface tension is chosen to be θ =0.01426.The pressure difference ?p= p∞?pb=0.0325 is obtained by artificially setting a density difference to trigger bubble collapse. Meantime, the single-component (SCP) model is added as the control group,that is,the bubble containing noncondensable gas is replaced by the vapor bubble in the computation domain. The initial densities of vapor and liquid are ρv=0.00306 and ρl=0.406 at 0.6Tc, respectively. The numerical simulation domain is 201×201. Periodic boundary conditions are adopted in all directions. The instantaneous radius is normalized by the initial radius R0, and the time is non-dimensionalized by the total simulation time step.

    Fig.3.Comparison among radius evolutions with time of collapsing bubble,obtained from three calculation methods.

    Figure 3 shows that the MRT-LBM simulation results accord with the results from the law of Rayleigh bubble collapse. It is verified that the multicomponent LBM method can simulate the dynamic behavior of a cavitation bubble containing non-condensable gas. Additionally, the slope of the blue dot curve is higher than that of the red dot curve. It is illustrated that the collapse velocity of the bubble containing noncondensable gas is lower than that of the vapor bubble. Thus,the non-condensable gas weakens the collapse velocity of a cavitation bubble.

    3.3. Validation of adiabatic law

    In this section,it is shown that the thermal model can be used to simulate the temperature evolution of the collapsing bubble. Because the collapse time of the cavitation bubble is nearly instantaneous (on the order of microseconds), it is generally considered that the process from cavitation bubble shrinkage to complete collapse is adiabatic. The formula for calculating the maximum temperature at the instant of bubble collapse under the adiabatic law is given by[1]

    where T∞is the ambient temperature in the liquid, R0is the initial bubble radius, and γ =1.3 is a constant representing the adiabatic coefficient of the ideal gas. The pressure difference ?p=0.0325. The evolution process of the temperature field of a collapsing cavitation bubble in an infinite domain and the corresponding central temperature curve of the bubble are shown in Fig.4. The instantaneous temperature is normalized by the critical temperature Tc, and the time is nondimensionalized by the total simulation time step.

    Fig.4. Plot of T/Tc versus t?,showing verification of adiabatic law.

    It is worth mentioning that the maximum temperature of the bubble collapse can theoretically reach infinity. This is because the theoretical model usually neglects the effects of heat dissipation and viscous dissipation.However,such a high temperature cannot be reached due to dissipation and other factors in the experiment or the LBM simulation.These methods have a common trend of temperature increase. As shown in Fig.4,the trend of the LBM numerical simulation is consistent with that of the theoretical solution.

    4. Effect of non-condensable gas in bubble on collapsing bubble near solid wall

    Due to the complexity and unpredictable nature of cavitation, the progress of revealing its thermal behavior and consequences need to be analyzed by multicomponent thermal LBM.The thermal effect of the collapsing bubble can be described more accurately by considering the action of noncondensable gas. In order to study the thermodynamic mechanism of the bubble containing non-condensable gas,a numerical model of cavitation bubble near a solid wall is established.

    4.1. Computational model

    The computational simulation domain of the cavitation bubble collapse near a solid wall is set to be a 201×201 lattice as shown in Fig.5, where R0is the initial radius of the bubble,b is the distance from the bubble center to the solid wall,and the offset parameter λ =b/R0is defined as the standoff distance between the bubble center and the lower wall. The bottom boundary is set to be a rigid plane boundary with a bounce-back boundary condition.[34]Periodic boundary conditions are adopted in the left and right direction.Additionally,an anti-bounce-back approach is adopted at the top of the computational domain by the open boundary.

    Fig.5. Computational domain.

    The initial bubble volume is constant, and thus the volume fraction of the gas and vapor are important adjustable parameters. Three cases of the initial volume fractions of the non-condensable gas and vapor in the bubble are given in Table 1. They are used to investigate the different morphologies and thermodynamics of collapsing under different volume ratios of the gas to the vapor. In this study, only the noncondensable gas and vapor volume in the bubble are considered to be variable,and the remaining initial conditions are invariant. In the simulation,the initial radius of the bubble is determined to be R0=40.The initial liquid density is ρl=0.406 at 0.6Tc. The pressure difference ?p= p∞?pb=0.0358. In our present work,[44,45]the obvious secondary collapse can be observed when the offset parameter λ =1.4 of the SCP bubble collapses near the wall. Therefore,the offset parameter is chosen to be λ =1.4 as the reference parameter in this MCP model.

    Table 1. Initial volume fraction of non-condensable gas and vapor in bubble.

    4.2. Cavitation bubble collapse near solid wall

    The results are presented, with the pressure field on the left and the temperature on the right,and due to the symmetry of the cavitation bubbles in Figs. 6–8, emphasis is placed on the temperature field.

    The cavitation bubble in case 1 can be considered as an almost vapor-filled bubble. As shown in Figs. 6(a)–6(c), the bubble deforms from a spherical bubble to an elongated bubble in a direction perpendicular to the solid wall. The radial liquid flow is retarded because of the solid wall, and a region with reduced pressure/temperature is formed near the solid wall. A conical high-pressure region is formed above the bubble due to the rebound of the liquid and the relatively high liquid velocity. The temperature inside the bubble is observed to increase gradually. The bubble collects the energy from the surrounding liquid,resulting in a lower temperature at the interface between the bubble and liquid. A micro-jet flow is formed on the top of the bubble as shown in Fig.6(d). The first collapse occurs as shown in Fig.6(e), and the jet directs toward the boundary and hits the lower bubble wall. Hence,the bubble becomes toroidal. Owing to the driving effect of the jet,the toroidal bubble rapidly collapses as shown in Fig.6(f),which is called the second collapse. Because the time interval between the first collapse and the second collapse is extremely short, the high-temperature and high-pressure shockwaves created by the second collapse hit the wall almost simultaneously with the first collapse as shown in Figs.6(g)and 6(h). The waves damage the wall and form a hotspot as shown in Figs.6(i)and 6(j). Then,the residual wave diffuses into the surrounding liquid until it vanishes.

    Figure 7 shows the simulation results for case 2. The evolution processes of the pressure and temperature fields before the first collapse are displayed in Figs. 7(a)–7(c) and are similar to those in Figs. 6(a)–6(e). After the first collapse (Figs. 7(d) and 7(e)), the high-temperature and highpressure shockwaves hit the solid wall perpendicularly. The second collapse is shown in Fig.7(f). The visible bubble wall disappears completely, and a high-pressure and hightemperature toroidal bubble is formed. Meanwhile, owing to the water hammer effect, a clear high-pressure and hightemperature domain is formed on the solid wall. The central temperature of the toroidal bubble is reduced. As shown in Fig.7(g), when the waves created by the high-pressure and high-temperature toroidal bubble propagate outward, another low-pressure toroidal region is formed in situ due to the rebound effect of the liquid.

    As shown in Figs. 7(g)–7(j), the pressure waves created by the second collapse are converged and superimposed on the symmetry axis,forming a high-pressure region,which immediately diffuses upward and downward. The temperature waves first diffuse the radiation waves from the collapse points to the surroundings. Then, the high-temperature waves reach the solid wall together with the high-pressure waves. Therefore,a hotspot region is formed on the wall,which intensifies the cavitation erosion. According to the second law of thermodynamics, energy can be transferred spontaneously from a high-temperature region to a low-temperature region. The high-temperature waves created by the bubble radiate in all directions and eventually vanish. The lower temperature region of the toroidal bubble gradually reaches the surrounding liquid temperature.

    Figure 8 shows the evolutions of pressure and temperature in case 3. The bubble is considered to be a non-condensable gas-filled bubble with little vapor. The evolution processes in Figs.8(a)–8(f)before the second collapse are similar to those in Figs. 7(a)–7(f). Except the collapse time increasing. The pressure wave and jet flow formed by the first collapse reach the bottom wall and rebound immediately(Fig.8(e)). Subsequently,a high-pressure and high-temperature toroidal bubble is formed at the second collapse as shown in Fig.8(f). The interval between the first collapse and the second collapse is longer. There is no timely liquid replenishment in the bottom region near the wall after the pressure waves has been created by the first collapse rebound (Fig.8(g)). A low-pressure region is formed on the bottom wall below the bubble.The pressure waves hit the wall driven by the liquid pressure gradient.Meanwhile, the high-temperature waves radiate outward and arrive at the solid wall together with the pressure waves,until they dissipate.

    Fig.6. Evolutions of pressure field and temperature field around collapsing bubble near solid wall(case 1).

    Fig.7. Evolutions of pressure field and temperature field around collapsing bubble near solid wall(case 2).

    Fig.8. Evolutions of pressure field and temperature field around collapsing bubble near solid wall(case 3).

    The above results indicate that the bubbles in all three cases undergo secondary collapse. However, the pressure field, temperature field evolution and collapse velocity of the bubble are different,and the damage to the wall is particularly diverse. Due to the percentage of gas volume in the bubble is the only variable in the simulation, the influencing factor for the cavitation bubble collapse is determined. Furthermore,the damage effects on the wall are analyzed in the following subsection.

    4.3. Effects of damage on wall

    The crucial roles of pressure and temperature in wall damage are shown by a comparison in Fig.10, which illustrates the time sequences of the two parameters pwand Tw,which are the pressure and the temperature against the rigid boundary,respectively. In order to explain and understand the phenomena shown in Fig.10,the bubble profiles and temperature field images of the first and second collapse are displayed in Fig.9.

    As can be seen from Fig.10(a), it is seemed that only one shock wave reaches the solid wall. However,as shown in Fig.9(a), under the same pressure difference, the bubble collapse velocity is fastest in case 1. The volume of the toroidal bubble is relatively small in Figs.9(a1)and 9(b1),so the collapse time of the toroidal bubble is relatively short. The time interval between the two collapses is essentially related to the deformation velocity of the collapsing cavitation bubble,which is also referred to as the collapse time of the toroidal bubble. Therefore,because the time interval between the first collapse and the second collapse is extremely short,the hightemperature and high-pressure shock waves created by the second collapse hit the wall almost simultaneously with the first collapse.

    Fig.9. Density field(left)and temperature field(right)at first and second collapse in three cases.

    Fig.10. Time sequences of pw and Tw on solid wall for different cases(λ =1.4).

    When the volume non-condensable in the bubble is larger,the bubble is not easily broken down by the jet flow, and the time step of the first collapse is longer. This indicates that the non-condensable gas has an impeding effect on the collapsing bubble. A larger volume of the toroidal bubble is shown in Fig.9(a3), and the corresponding collapse time increases.Although the bubble profiles in Figs.9(b2)and 9(b3)are similar, the surrounding density distributions are quite different.As the interval increases, the shockwaves created by the first collapse already rebound off the bottom wall in all directions.There is no timely liquid replenishment in the bottom region near the wall. The pressure waves created by the second collapse hit the wall,driven by the liquid pressure gradient,which causes pwto increase sharply as shown in Fig.10(c). In addition, as the collapse time increases, the temperature dissipation and exchange between the bubble and the surrounding liquid also increase. The internal temperature of the bubble decreases. Therefore,the maximum Twfor case 1 is the highest,and the maximum Twfor case 3 is the lowest in the three cases.

    Thus, the non-condensable gas in the bubble can affect the deformation velocity and the profile of the bubble, which results in variations in the thermodynamic and kinetic behavior of the bubble. When the non-condensable gas concentration increases,the distinction of the pressure and temperature on the wall after the second collapse are more obvious. In addition, it is found that the thermodynamic behavior of the bubble is not completely synchronized with the dynamic behavior due to the thermal dissipation in the collapse process that needs further considering.

    5. Discussion and conclusions

    In this study, a multicomponent thermal MRT-LBM method is adopted to investigate cavitation bubble collapse with non-condensable gas. The temperature changes of the two-component two-phase flow are well-captured by the LBM. Furthermore, it is demonstrated that the multicomponent thermal MRT-LBM method is an effective tool for simulating multicomponent multiphase flow. Some conclusions can be drawn as follows.

    (i) The numerical results obtained by the present model comply with the Laplace law and adiabatic law. In addition,the LBM simulation results of cavitation bubble collapse are consistent with the numerical solution of the RP equation for the partial pressure gas.

    (ii) A numerical model of cavitation bubble near a solid wall is established. Three cases of different initial volume fractions of non-condensable gas and vapor in the bubble are considered in the simulation to investigate the different morphologies and thermodynamic behavior of collapsing bubble.The pressure and temperature field evolution provide clear physical pictures for understanding the mechanism of bubble collapse near the wall. The results show that the noncondensable gas concentration in the bubble significantly influences the collapse consequence.

    (iii) The pressure and temperature against the solid wall during the first collapse and the second collapse are explained by analyzing the bubble profiles of each collapse in detail.When considering the non-condensable gas, the thermodynamic and kinetic behavior of the bubble vary. The distinction of the pressure and temperature on the wall after the second collapse are more obvious as the non-condensable gas concentration increases. Therefore,it is of great importance to study the effect of non-condensable gas inside bubble on cavitation.

    The influence of the non-condensable gas in the bubble can provide an insight into the engineering or industry of the cavitation applications. Future work will focus on the cavitation bubble with mixture non-condensable gas and vapor collapse near a thermal/cold wall.

    猜你喜歡
    韓慶
    Effect of porous surface layer on wave propagation in elastic cylinder immersed in fluid
    滲漏是三峽庫區(qū)砂質(zhì)土橘園氮磷流失的主要途徑*
    “丹麥小馬達”環(huán)球航線:愛跟上奮斗的節(jié)奏
    Effect of viscosity on stability and accuracy of the two-component lattice Boltzmann method with a multiple-relaxation-time collision operator investigated by the acoustic attenuation model
    林海守望者
    韓慶玲:紅手印留下“鐵腿女書記”
    韓慶芳:我站立的地方, 就是我的陣地
    當代黨員(2020年10期)2020-07-04 12:33:47
    基于Logistic人口發(fā)展模型方程的安徽省人口問題研究
    Investigation of cavitation bubble collapse near rigid boundary by lattice Boltzmann method*
    七月的最后一天,或忘記
    參花(上)(2013年9期)2013-06-10 15:40:50
    久久午夜福利片| 超碰av人人做人人爽久久| 丰满乱子伦码专区| 亚洲自拍偷在线| 欧美又色又爽又黄视频| 欧美日韩综合久久久久久 | 亚洲无线在线观看| 精品久久国产蜜桃| 99在线人妻在线中文字幕| 国内揄拍国产精品人妻在线| 国产乱人视频| 亚洲精品在线观看二区| 婷婷精品国产亚洲av在线| 成年女人看的毛片在线观看| 国产黄片美女视频| 在线观看免费视频日本深夜| 日韩精品中文字幕看吧| 日日摸夜夜添夜夜添av毛片 | 99久久成人亚洲精品观看| 夜夜夜夜夜久久久久| 精品午夜福利视频在线观看一区| 日本黄色片子视频| 久久精品国产99精品国产亚洲性色| 亚洲av免费在线观看| 久久久久亚洲av毛片大全| 此物有八面人人有两片| 久久欧美精品欧美久久欧美| 国产蜜桃级精品一区二区三区| 国产熟女xx| 国产亚洲精品av在线| 熟妇人妻久久中文字幕3abv| 国产在线男女| 国产成人a区在线观看| av在线老鸭窝| 欧美黄色淫秽网站| 欧洲精品卡2卡3卡4卡5卡区| 在线播放国产精品三级| 欧美成人免费av一区二区三区| 特大巨黑吊av在线直播| 欧美日韩亚洲国产一区二区在线观看| 校园春色视频在线观看| 欧美成人免费av一区二区三区| 日本在线视频免费播放| 国产大屁股一区二区在线视频| 内地一区二区视频在线| 日韩欧美在线二视频| 一本精品99久久精品77| 脱女人内裤的视频| 又黄又爽又刺激的免费视频.| 99热这里只有是精品50| 97超级碰碰碰精品色视频在线观看| 久久久国产成人免费| 亚洲欧美激情综合另类| 久久99热6这里只有精品| 成人亚洲精品av一区二区| 国产精品久久电影中文字幕| 直男gayav资源| 国产精品永久免费网站| 午夜久久久久精精品| 国产v大片淫在线免费观看| 亚洲av不卡在线观看| 亚洲男人的天堂狠狠| 精品久久久久久久久久久久久| 亚洲avbb在线观看| 老女人水多毛片| 99精品在免费线老司机午夜| 免费在线观看日本一区| 日日夜夜操网爽| 亚洲,欧美,日韩| 白带黄色成豆腐渣| 成年女人看的毛片在线观看| 国产精品98久久久久久宅男小说| 亚洲成人免费电影在线观看| 两个人的视频大全免费| 欧美日本亚洲视频在线播放| 国产私拍福利视频在线观看| 18禁在线播放成人免费| 舔av片在线| 丝袜美腿在线中文| 一边摸一边抽搐一进一小说| 能在线免费观看的黄片| 亚洲一区高清亚洲精品| 久久草成人影院| 久久久久久久亚洲中文字幕 | 美女高潮喷水抽搐中文字幕| 男女床上黄色一级片免费看| 亚洲av电影在线进入| 欧美性猛交╳xxx乱大交人| 两个人视频免费观看高清| 欧美成人一区二区免费高清观看| 别揉我奶头~嗯~啊~动态视频| 首页视频小说图片口味搜索| 蜜桃亚洲精品一区二区三区| 国产综合懂色| 波多野结衣巨乳人妻| netflix在线观看网站| 亚洲国产精品合色在线| 18+在线观看网站| 内射极品少妇av片p| 免费av不卡在线播放| 国产一级毛片七仙女欲春2| 村上凉子中文字幕在线| 高清日韩中文字幕在线| 国内毛片毛片毛片毛片毛片| 18禁裸乳无遮挡免费网站照片| 中文字幕精品亚洲无线码一区| 99久国产av精品| 淫妇啪啪啪对白视频| 久久久久免费精品人妻一区二区| 国产精品1区2区在线观看.| 精品国产亚洲在线| 亚洲国产精品sss在线观看| 我的女老师完整版在线观看| 亚洲性夜色夜夜综合| 久久天躁狠狠躁夜夜2o2o| 可以在线观看的亚洲视频| 一级av片app| 国产成人啪精品午夜网站| 日日夜夜操网爽| 日韩精品中文字幕看吧| 一夜夜www| 亚洲精品久久国产高清桃花| 很黄的视频免费| 婷婷丁香在线五月| 啦啦啦观看免费观看视频高清| 18+在线观看网站| 国产老妇女一区| av在线观看视频网站免费| 欧美黑人欧美精品刺激| 亚洲av第一区精品v没综合| 最后的刺客免费高清国语| 在线看三级毛片| 午夜久久久久精精品| 伊人久久精品亚洲午夜| 欧美日韩黄片免| 91av网一区二区| 他把我摸到了高潮在线观看| 黄色视频,在线免费观看| 在线观看午夜福利视频| 国产中年淑女户外野战色| 国产av一区在线观看免费| 亚洲国产日韩欧美精品在线观看| 亚洲国产高清在线一区二区三| 国产av在哪里看| 欧美精品国产亚洲| 亚洲黑人精品在线| 色精品久久人妻99蜜桃| 丰满的人妻完整版| 一个人免费在线观看电影| 97热精品久久久久久| 99热这里只有是精品在线观看 | 久久精品夜夜夜夜夜久久蜜豆| 亚洲18禁久久av| 男女之事视频高清在线观看| 好男人在线观看高清免费视频| 亚洲av不卡在线观看| 欧美最黄视频在线播放免费| 少妇高潮的动态图| 2021天堂中文幕一二区在线观| 婷婷六月久久综合丁香| 中文字幕熟女人妻在线| 亚洲国产精品成人综合色| 男女视频在线观看网站免费| 欧美日韩福利视频一区二区| 又爽又黄a免费视频| 伊人久久精品亚洲午夜| 91字幕亚洲| 黄色女人牲交| 琪琪午夜伦伦电影理论片6080| 给我免费播放毛片高清在线观看| 一级av片app| 欧美高清成人免费视频www| 国产一级毛片七仙女欲春2| 国产一区二区在线观看日韩| 成人特级av手机在线观看| 国产精华一区二区三区| 美女高潮的动态| 国产av一区在线观看免费| 91狼人影院| 精品一区二区三区人妻视频| 日韩欧美一区二区三区在线观看| 黄色视频,在线免费观看| 午夜福利视频1000在线观看| 欧洲精品卡2卡3卡4卡5卡区| 99国产精品一区二区三区| 欧美日韩黄片免| 婷婷六月久久综合丁香| 一个人免费在线观看电影| 麻豆国产av国片精品| 特大巨黑吊av在线直播| 9191精品国产免费久久| av在线蜜桃| 在线观看一区二区三区| 国产精品自产拍在线观看55亚洲| 亚洲最大成人手机在线| 国产精品伦人一区二区| 中文资源天堂在线| 在线观看av片永久免费下载| 高清日韩中文字幕在线| 国产精品久久久久久人妻精品电影| 性欧美人与动物交配| 国产在线男女| 两个人的视频大全免费| 亚州av有码| 国模一区二区三区四区视频| 两人在一起打扑克的视频| 国产 一区 欧美 日韩| 国产白丝娇喘喷水9色精品| 亚洲自拍偷在线| 欧美日韩中文字幕国产精品一区二区三区| 午夜精品在线福利| 蜜桃久久精品国产亚洲av| 伊人久久精品亚洲午夜| 久久国产乱子伦精品免费另类| 一个人免费在线观看的高清视频| 精品一区二区三区视频在线| 午夜精品一区二区三区免费看| 在线a可以看的网站| 亚洲国产精品999在线| 综合色av麻豆| 熟女人妻精品中文字幕| www.熟女人妻精品国产| 一个人看的www免费观看视频| 中文字幕人成人乱码亚洲影| 成人午夜高清在线视频| 男人狂女人下面高潮的视频| 免费一级毛片在线播放高清视频| 午夜福利在线观看免费完整高清在 | 欧美3d第一页| 熟妇人妻久久中文字幕3abv| 午夜影院日韩av| 深爱激情五月婷婷| 亚洲第一区二区三区不卡| 久久99热这里只有精品18| 国产黄色小视频在线观看| 少妇人妻精品综合一区二区 | 天堂av国产一区二区熟女人妻| 国产精品久久久久久久电影| 亚洲精品成人久久久久久| 偷拍熟女少妇极品色| 国产精品嫩草影院av在线观看 | 性插视频无遮挡在线免费观看| 亚洲avbb在线观看| 熟妇人妻久久中文字幕3abv| 最好的美女福利视频网| 琪琪午夜伦伦电影理论片6080| 国产精品美女特级片免费视频播放器| 可以在线观看毛片的网站| 亚洲成人免费电影在线观看| 欧美黄色片欧美黄色片| a在线观看视频网站| 国产aⅴ精品一区二区三区波| 国产精品三级大全| 亚洲欧美日韩高清专用| 窝窝影院91人妻| ponron亚洲| 国产欧美日韩精品亚洲av| 亚洲欧美清纯卡通| 午夜久久久久精精品| 日本成人三级电影网站| 我的女老师完整版在线观看| 男女之事视频高清在线观看| 国产真实乱freesex| 99精品在免费线老司机午夜| 日韩欧美 国产精品| 九色成人免费人妻av| 精品人妻熟女av久视频| 十八禁国产超污无遮挡网站| 精品福利观看| 欧美日韩黄片免| 久久人人爽人人爽人人片va | 嫩草影院精品99| 51午夜福利影视在线观看| 精品日产1卡2卡| 国产精品一区二区三区四区免费观看 | 日韩亚洲欧美综合| 99久久成人亚洲精品观看| 小说图片视频综合网站| 亚洲av美国av| 亚洲五月婷婷丁香| 高清在线国产一区| av福利片在线观看| 国产欧美日韩一区二区三| 国产午夜福利久久久久久| 五月玫瑰六月丁香| 日韩亚洲欧美综合| 亚洲成av人片在线播放无| 日韩欧美国产一区二区入口| 俄罗斯特黄特色一大片| 一级黄片播放器| 国产亚洲精品综合一区在线观看| 此物有八面人人有两片| 欧美一区二区亚洲| 亚洲精品456在线播放app | 欧美不卡视频在线免费观看| 午夜a级毛片| 毛片女人毛片| 亚洲成av人片在线播放无| 精品久久久久久久久亚洲 | 淫秽高清视频在线观看| 夜夜夜夜夜久久久久| 日韩国内少妇激情av| 亚洲片人在线观看| 最后的刺客免费高清国语| 国产精品久久电影中文字幕| 婷婷精品国产亚洲av| 亚洲综合色惰| 国产不卡一卡二| 成人午夜高清在线视频| 99精品在免费线老司机午夜| 亚洲精品成人久久久久久| 欧美激情在线99| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 午夜影院日韩av| 制服丝袜大香蕉在线| 成年免费大片在线观看| 特级一级黄色大片| 久久精品夜夜夜夜夜久久蜜豆| 久久久色成人| 国产激情偷乱视频一区二区| 1000部很黄的大片| 九色成人免费人妻av| 18禁黄网站禁片午夜丰满| 18+在线观看网站| 97碰自拍视频| 成人特级av手机在线观看| 又爽又黄无遮挡网站| 国产探花极品一区二区| avwww免费| 亚洲精品久久国产高清桃花| 亚洲欧美日韩高清在线视频| 国产高潮美女av| 我要看日韩黄色一级片| 女人被狂操c到高潮| 网址你懂的国产日韩在线| av欧美777| 桃红色精品国产亚洲av| 欧美成人免费av一区二区三区| av天堂中文字幕网| 91麻豆精品激情在线观看国产| 日韩高清综合在线| 淫妇啪啪啪对白视频| av在线观看视频网站免费| 91久久精品国产一区二区成人| 又爽又黄a免费视频| 每晚都被弄得嗷嗷叫到高潮| 色综合站精品国产| 国内精品久久久久久久电影| 99久久成人亚洲精品观看| 免费在线观看日本一区| 成人毛片a级毛片在线播放| 深爱激情五月婷婷| 国产精品永久免费网站| 欧美+亚洲+日韩+国产| 欧美一区二区亚洲| 日韩欧美在线二视频| 成人毛片a级毛片在线播放| 精品久久久久久久久av| 欧美色欧美亚洲另类二区| 看十八女毛片水多多多| 天堂动漫精品| 精品一区二区三区视频在线| 免费av观看视频| 国产精品久久久久久精品电影| 亚洲成人久久爱视频| 国产精品一区二区三区四区免费观看 | 免费电影在线观看免费观看| 久久久国产成人免费| 757午夜福利合集在线观看| 他把我摸到了高潮在线观看| 99在线视频只有这里精品首页| 伊人久久精品亚洲午夜| av天堂在线播放| 日韩欧美精品免费久久 | 国产91精品成人一区二区三区| 99热6这里只有精品| 小说图片视频综合网站| 丰满乱子伦码专区| 中文字幕免费在线视频6| 精品一区二区三区人妻视频| 国产黄a三级三级三级人| 亚洲第一欧美日韩一区二区三区| 全区人妻精品视频| 亚洲精品粉嫩美女一区| 欧美日韩亚洲国产一区二区在线观看| 欧美区成人在线视频| 欧美性感艳星| 国产成人啪精品午夜网站| 麻豆久久精品国产亚洲av| 色播亚洲综合网| 老司机午夜福利在线观看视频| 久久香蕉精品热| 男女之事视频高清在线观看| 免费观看精品视频网站| 亚洲成人精品中文字幕电影| 国产成人影院久久av| 乱码一卡2卡4卡精品| 成人三级黄色视频| 给我免费播放毛片高清在线观看| 国产成+人综合+亚洲专区| 在线免费观看的www视频| 欧美黄色淫秽网站| 老熟妇乱子伦视频在线观看| 亚洲内射少妇av| 色尼玛亚洲综合影院| 好男人电影高清在线观看| 午夜福利在线观看吧| 国产亚洲精品av在线| 久久精品国产清高在天天线| 国产精品伦人一区二区| 久久国产精品影院| 欧美性猛交╳xxx乱大交人| 成人美女网站在线观看视频| 五月伊人婷婷丁香| 日韩 亚洲 欧美在线| 国产午夜福利久久久久久| 亚洲av第一区精品v没综合| 国产色婷婷99| 日本 av在线| 亚洲国产色片| 又黄又爽又免费观看的视频| www.色视频.com| 99久久精品国产亚洲精品| 亚洲一区二区三区色噜噜| 最近视频中文字幕2019在线8| 亚洲人成网站高清观看| 18禁黄网站禁片免费观看直播| 欧美色视频一区免费| 免费看光身美女| 舔av片在线| 在线播放国产精品三级| 午夜影院日韩av| 国产精华一区二区三区| 日本熟妇午夜| 色综合站精品国产| av在线老鸭窝| 97人妻精品一区二区三区麻豆| 国产欧美日韩一区二区精品| 国语自产精品视频在线第100页| 亚洲最大成人av| 亚洲av成人精品一区久久| 国产免费av片在线观看野外av| 日本黄色片子视频| 日韩人妻高清精品专区| 色噜噜av男人的天堂激情| 色综合站精品国产| 男女那种视频在线观看| 午夜福利18| 简卡轻食公司| 国产一区二区三区视频了| 国产精品久久视频播放| 熟妇人妻久久中文字幕3abv| 色哟哟·www| 一卡2卡三卡四卡精品乱码亚洲| 亚洲成人久久性| 国模一区二区三区四区视频| 色5月婷婷丁香| 亚洲18禁久久av| 日韩欧美在线二视频| 亚洲自拍偷在线| 久久精品影院6| 久久人妻av系列| 国产欧美日韩一区二区三| 热99re8久久精品国产| 黄色一级大片看看| 久久久久久大精品| 国产 一区 欧美 日韩| 村上凉子中文字幕在线| 18禁裸乳无遮挡免费网站照片| 搞女人的毛片| 亚洲无线在线观看| a级一级毛片免费在线观看| 亚洲18禁久久av| 国内精品久久久久精免费| 久久人人爽人人爽人人片va | 精品人妻1区二区| 国产视频内射| 国产精品久久久久久亚洲av鲁大| 日韩有码中文字幕| 欧美激情久久久久久爽电影| 国产精品精品国产色婷婷| 日日摸夜夜添夜夜添av毛片 | a在线观看视频网站| 国产视频一区二区在线看| 日韩高清综合在线| 亚洲成人免费电影在线观看| 亚洲精品影视一区二区三区av| 亚洲经典国产精华液单 | 国产精品久久久久久久电影| 欧美xxxx黑人xx丫x性爽| 国产精品人妻久久久久久| 亚洲激情在线av| 在线播放国产精品三级| 无遮挡黄片免费观看| 精品免费久久久久久久清纯| 亚洲,欧美,日韩| 超碰av人人做人人爽久久| 色5月婷婷丁香| 国产精品女同一区二区软件 | 欧美最新免费一区二区三区 | 久久久精品大字幕| 国产成年人精品一区二区| 国内精品久久久久精免费| 长腿黑丝高跟| 免费观看精品视频网站| 久久精品影院6| 美女被艹到高潮喷水动态| 精品久久国产蜜桃| a级一级毛片免费在线观看| 女生性感内裤真人,穿戴方法视频| 免费一级毛片在线播放高清视频| 亚洲熟妇熟女久久| 日日摸夜夜添夜夜添小说| 琪琪午夜伦伦电影理论片6080| 天堂av国产一区二区熟女人妻| 国产精品日韩av在线免费观看| 看黄色毛片网站| 国产一区二区三区在线臀色熟女| 亚洲无线观看免费| 少妇裸体淫交视频免费看高清| 亚洲欧美日韩东京热| 一夜夜www| 黄色一级大片看看| 午夜福利在线观看吧| 亚洲精品在线美女| 1024手机看黄色片| 亚洲专区中文字幕在线| 中文字幕精品亚洲无线码一区| 自拍偷自拍亚洲精品老妇| 国产单亲对白刺激| avwww免费| 午夜久久久久精精品| 亚洲熟妇中文字幕五十中出| 小蜜桃在线观看免费完整版高清| 深爱激情五月婷婷| 人妻制服诱惑在线中文字幕| 成人无遮挡网站| 亚洲一区二区三区色噜噜| 热99re8久久精品国产| 日韩成人在线观看一区二区三区| 校园春色视频在线观看| 久久久成人免费电影| av黄色大香蕉| 久久精品影院6| 亚洲午夜理论影院| 51午夜福利影视在线观看| 国产黄片美女视频| 国产亚洲精品久久久com| 一区二区三区高清视频在线| 很黄的视频免费| 精品人妻1区二区| 久久草成人影院| 色综合亚洲欧美另类图片| 国产综合懂色| 麻豆国产97在线/欧美| 亚洲成av人片在线播放无| 69人妻影院| 99视频精品全部免费 在线| 欧美日韩福利视频一区二区| 村上凉子中文字幕在线| 亚洲国产精品久久男人天堂| 国产精品一区二区三区四区久久| 国产在线男女| 精品免费久久久久久久清纯| 免费av不卡在线播放| 男人的好看免费观看在线视频| 日韩欧美 国产精品| xxxwww97欧美| 欧美在线黄色| 国产 一区 欧美 日韩| 国产黄a三级三级三级人| 欧美日本视频| 99国产精品一区二区蜜桃av| av福利片在线观看| 亚洲av免费高清在线观看| 波多野结衣高清无吗| 如何舔出高潮| 美女黄网站色视频| 啦啦啦韩国在线观看视频| 国内毛片毛片毛片毛片毛片| 一个人看的www免费观看视频| 久久久色成人| 国产亚洲精品av在线| 欧美午夜高清在线| 啦啦啦观看免费观看视频高清| 国内毛片毛片毛片毛片毛片| 又紧又爽又黄一区二区| av女优亚洲男人天堂| 一区二区三区高清视频在线| 国产精品嫩草影院av在线观看 | 蜜桃久久精品国产亚洲av| 少妇人妻一区二区三区视频| 亚洲第一电影网av| 色哟哟·www| 国产免费男女视频| 日本一二三区视频观看| 此物有八面人人有两片| 亚洲人成网站在线播| 国产av不卡久久| 特级一级黄色大片| 男插女下体视频免费在线播放| 日韩欧美国产一区二区入口| 国产久久久一区二区三区| 国产一区二区亚洲精品在线观看| 国产成+人综合+亚洲专区| www.www免费av| 亚洲av成人不卡在线观看播放网| 精品日产1卡2卡| 国产高清视频在线观看网站| 一本综合久久免费| 国产成人影院久久av| 午夜免费成人在线视频| 亚洲第一电影网av| a级毛片a级免费在线| 亚洲经典国产精华液单 |