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

    Level Set Method for Simulation of Cavitating Flows

    2011-06-07 07:52:18
    船舶力學(xué) 2011年3期
    關(guān)鍵詞:北京理工大學(xué)相間空化

    (School of Vehicle and Transportation Engineering,Beijing Institute of Technology,Beijing 100081,China)

    Level Set Method for Simulation of Cavitating Flows

    HUANG Biao,WANG Guo-yu,ZHANG Min-di,CHEN Wei-qing

    (School of Vehicle and Transportation Engineering,Beijing Institute of Technology,Beijing 100081,China)

    In this study,the Level set method is presented and applied to multiphase cavitating flows.Herein,this method is applied to flow solver that fully couples the mass,momentum,and transport equation;in numerical simulation,the Level set method with a compressive discretization scheme in both time and space to minimize the smearing of the free surface at the interface is used for simulation of nature and ventilated cavity.In these computations,three cases,including the steady nature cavitating flow around the axisymmetric body and the Clark-Y hydrofoil,and the ventilated cavitating flow around the disk are conducted,compared with the results obtained by simulation with the widelyused mixture model and experimental data,it is shown that,the Level set method is presented specifically to cavitating flows that decreases the interface diffusion where needed.Such an approach is a valid method that could reduce smearing at the cavity interface and also capture the sharp interface between different phases in an efficient and precise manner.

    cavitating flow;mixture model;Level set method

    Biography:HUANG Biao(1985-),male,Ph.D.student of Beijing Institute of Technology,E-mail:huangbiao@bit.edu.com.

    1 Introduction

    Cavitating and capturing the motion of the cavity interface in multi-phase fluid flows are a challenging research field of computational fluid dynamics.The motions of cavity interfaces(nature or ventilated supercavities)are important physical phenomenon in many physically interesting.A number of methods have been developed to resolve and track this interface in Computational Fluid Dynamics(CFD)simulating multifuild flows with sharp interfaces listing advantages and disadvantages each.The Level Set method and VOF model are becoming more and more popular to capture the motion of a free surface.The volume of fluid(VOF)method has been widely applied in tracking free interface since it was proposed by Hirt and Nichols[1],in the VOF method,the volume fraction function C is defined first,and the value of C for each cell corresponds to the fraction of fluid filling each cell.Here,the free interface can then be tracked using the interface reconstruction technique and advection algorithm.Important progress has been made in the VOF method by many researchers,such as the SLIC by Noh and Woodward[2],the FLAIR by Ashgriz and Poo[3],the refined reconstruction by Youngs[4],and the CIC-SAM by Ubbink and Issa[5].In parallel,the proposed Level-set formulation is based on the method for interface-fluid flow by reinitializing the interface using a signed distance from the distance.The advantage of the Level set method is the maintenance of the sharp interface which in the VOF method can only be remained by the complicated reconstruction of the interface[6].Numerical algorithms for the Level set equation have been further developed over the past decade.Sussman et al[7]applied the Level set method to incompressible fluid flows in their research;the Level set method is combined for solving the mass conservation equation and the momentum transport equation.Additionally,Shu Bitan et al[8]used the Level set method into OpenFOAM capturing the free interface in incompressible fluid flows;the terminal velocity of rising air bubble in water and the detachment of a bubble from a wall are successfully simulated.For caviating flows,the species-mass-conservation-based Level set method was valued for natural and ventilated cavitating flows.

    In present work,in order to capture the distinct interface between different phases in cavitating flows,we provide the results for the validation study of the Level set method.This method is tested on several selected problems for which there are available experimental data of nature and ventilated cavitating flows.Three different test cases for modeling the steady and transient free surface tension have been considered.Special attention is paid to the advantages of this method compared with the numerical results with the other multiphase model-mixture model.

    2 Physical and numerical model

    The set of governing equations consists of the conservative form of the Reynolds averaged Navier-Stokes equations,plusing a volume fraction transport equation to account for the cavitation dynamics,and the turbulence closure.

    2.1 Mixture model

    The Navier-Stokes equations in their conservative form governing a Newtonian fluid without body forces and heat transfers are presented below in the Cartesian coordinates:

    The mixture property,φm,can be expressed as:

    where φ can be density,viscosity,and so on.The mixture model is very simple which treats both phases symmetrically.It is appropriate for a calculation of non-disperse liquid-liquid or gas-liquid two phase flow.

    2.2 Level set method with mass transfer

    The free surface flow model is a homogeneous Eulerian-Eulerian multiphase model,in the Level set method,the interface between the two phases is represented by a continuous scalar function φ( x, t),which is set to zero on the interface,is positive on one side,and negative on the other side.This way both phases are identified,and the location of the physical interface is associated with the surface φ=0.The function φ is called the Level set function and is typically defined as the signed distance to the interface:i.e.,φ=-d( x, t ) on one side of the interface and φ=+d( x, t ) on the other,where d( x, t ) is the shortest distance from the point x to the interface[6].

    When the interface is advected by the flow,the evolution of the Level set function is given by:

    In the Level set fluid-fluid formulation,the density and viscosity are typically interpolated across the interface as follows:

    where subscripts 1 and 2 denote the values corresponding to the two different phases,respectively.Here,Hε(φ ( x ))is a smoothed Heaviside function which is defined in Ref.[7].

    When the interface between the fluids is highly curved,then the effect of the surface tension force in the multiphase simulation becomes particularly important.Surface tension acts at the interface between two fluids.Computationally,this is awkward to deal with and the role of surface tension is generally incorporated as a continuum surface force.Using the continuum surface force model,the momentum transport equation for incompressible flow is:

    where,the model is based on the continuum surface force model of Brackbill et al[10]which models the surface tension force as a volume force concentrated at the interface.σ is the surface tension coefficient,n is the interface normal vector pointing from the primary fluid to the second fluid(calculated using the volume fraction gradient),and κ is the surface curvature defined by:

    For nature cavitating flows,mass transfer is performed within a volume-fraction-based Level set method:

    where ρφis the mixture density,ujis the velocity component in Cartesian coordinates,t is the time,+is the condensation rate and˙-is the evaporation rate.In presented study,the source term˙+and˙-denote vapor generation(evaporation)and condensation rates which are depicted in Ref.[11].

    In computations,the Level set method uses a compressive discretization scheme in both time and space to minimize the smearing of the free surface at the interface.Typically,this reduces smearing at the interface to two or three layers of cells.The Level-set method is an interface capturing approach rather than an interface reconstruction approach.

    2.3 Turbulence model

    In present study,the SST model[12]is chosen for steady simulation and the DES model based on the SST formulation is chosen for transient simulation.The advantage of this combination is that the accurate prediction of turbulent boundary layers up to separation and in separated regions carries over from the SST model.In addition,the SST model supports the formulation of a zonal DES formulation[13],which is less sensitive to grid resolution restrictions than the standard DES formulation.

    Detached-Eddy Simulation model:

    All the coefficients are listed for completeness in Ref.[13],and the eddy-viscosity is defined as:

    The length scales of the model in terms of k and ω read:

    The idea behind the DES model for Strelets[13]is to switch the SST-RANS model to a LES model in regions where the turbulent length,predicted by the RANS model is larger than the local grid spacing.In this case,the length scale used in the computation of the dissipation rate in the equation for the turbulent kinetic energy is replaced by the local grid spacing.Here,Δ is based on the largest dimension of the grid cell,Δ=max( Δx,Δy,Δz).

    For wall-bounced separated flows,the above formulation results in a hybrid model that functions as the standard SST model inside the whole attached boundary layer,and as its sub grid-scale version in the near wake.

    3 Results and discussion

    3.1 Steady natural cavitation on cylindrical bodies

    Rouse and McNown[14]carried out a series of experiments on natural cavitation about ax-isymmetric configurations.For each configuration,measurements were made across a range of cavitation numbers.The hemispherical configuration is computed here,the experiments were performed at Reynolds number greater than 100 000 based on the diameter.As shown in Fig.1,the boundary conditions used in the present simulations include inflow,outflow and no-slip,standard boundary conditions for incompressible flow are applied:the velocity is imposed at the inlet and the pressure is fixed at the domain outlet.A no-slip boundary condition is used at both the upper and lower walls.Here,an appropriate grid of dimension is used,as shown in Fig.2,the O type mesh is chosen because of the rounded leading edge.Based on the diameter of the cylinder,a value of Re=1.36×105is used for simulations in computations,and the diameter of the cylinder is 0.02m.

    Fig.1 Outline of the computational domain

    Fig.2 Computational grids around body

    Fig.3(a)and(b)show flowfield snapshots of the steady cavitating flow over a hemispherical axisymmetric body,the contours of volume fraction of vapor phase are shown respectively when the cavitation number is 0.3.For all cases,the presences of the bubble manifest itself as a decrease in magnitude,flattering and length-ening of the pressure minimum along the surface.Also,bubble closure gives rise to an overshoot in pressure recovery due to the local stagnation associated with free-stream liquid flowing over the convex curvature at the aft end of the bubble.But the cavity interfaces captured with different model are noticeable different,in the downstream region of the nature cavity,the cavity interface obtaind by Level set method is sharper.In addition,in the cavitation area,the vapor volume fraction is moderately larger than the results obtained with the mixture model.

    Fig.3(c)shows the comparisons between predicted and measured surface pressure distributions in parallel.Both the models match the experimental data satisfactorily.Differences in the performance are more noticeable in the closure region,where the vapor phase condenses.Fig.3(d)shows the corresponding density distribution along the surface.As seen from density plots,the liquid phase first expands and vapor phase appears uniformly inside the cavity,then the vapor phase is compressed,in a shock like fashion,back to the liquid phase.The differences among the two models in density profiles are significant.This implies that each model generates different compressibility characteristics,although they produce very similar steadystate pressure distributions;which model produces the correct compressibility is an open question and needs further investigation.

    Fig.3 Contour of volume fraction of vapor phase and density distribution with different models(σ=0.3)(a)Contour of volume fraction of vapor phase with Level set method(b)Contour of volume fraction of vapor phase with mixture model(c)Comparisons between predicted and measured surface pressure distributions(d)Density distribution along the surface with different models

    3.2 Unsteady natural cavitation around a Clark-Y hydrofoil

    Unsteady nature cavitating flow over the Clark-Y hydrofoil is investigated using the Level set method.In such cavitating case,liquid-vapor mass transfer occurs through local fluctuations about the saturation pressure.Here,the geometry is based on the Clark-Y profiles,has a 0.07m chord,and is twisted a total of 8 degrees.The computational model setup approximates to the experiments of Wang et al[15],simulating the hydrofoil within inviscid wall.In the experiment,the tunnel was run various speeds and pressures,visual comparisons are made with the results at V∞=10m/s,and σ=0.8.

    Fig.4 Outline of the computational domain

    Fig.5 Computational grids around a hydrofoil

    Fig.6 displays comparisons of the unsteady cavity around the Clark-Y hydrofoil between the mixture model and the Level set method,the temporal evolution of the computed and experimentally observed flow structures with cloud cavitation(σ=0.8)is shown.Here,the experimental visual images are shown in Fig.6(a),the black and white images correspond to an instantaneous visualization of one cycle,and Fig.6(b)and(c),show the time sequences of flow structures predicted by different models,the blue color corresponds to liquid and the red color corresponds to a highly concentrated bubble flow.

    Fig.6 Time evolution of cloud cavitation with the numerical results and experimental visuals

    There is a good agreement between the numerical and experimental results concerning the external shape and global structure of both attached cavity and vapor cloud shedding is achieved;but the vapor/water interface and vapor volume fraction with different models are distinct.The cavity got by Level set method is more obvious with larger vapor volume fraction.And also,the Level set method scheme prediction displayed the increased cavity size.This is because,for the Level set method in present study,the anti-diffusive characteristics of this scheme predict the nature cavity with much sharper,and hence it provides crisp interface for both steady-state and transient simulations.

    3.3 Axisymmetric ventilated disk cavitor

    The case relevant to the cavitating flows is ventilated cavitating flows around an axisymmetric disk cavitator.In this case,liquid flow around the disk and a gas injected into the wake form supercavities.The geometric configuration is displayed below in Fig.7,with the disk positioned such that the flat front is perpendicular to the free-stream velocity vector of a liquid flow.Such a case,where the air gas is injected into a liquid flow to form a cavity,is referred to as artificial or ventilated cavitation,the diameter of disc Dnis 15mm.In numerical simulation,the grid distribution is shown in Fig.8.

    Fig.7 Schematic diagram and computational domain

    Fig.8 Computational grid for flow analysis

    There are two cases predicted with different multiphase models,the assumed conditions for the cases are:for case-1,U∞=8m/s,and σ=2.37,ventilated rate Q is 400L/h;for case-2,U∞=3m/s,and σ=22.56,ventilated rate Q is 600L/h.

    Solution comparisons with and without the Level set method are evaluated in Fig.9.In our computation,the gas volume fraction isosurface 0.8 is defined as ventilated cavity.In both cases,without the interface sharpening of Level set method,the predicted profiles of the ventilated cavity tend to converge on a relatively smaller cavity size to the experimental visualizations and the numerical results obtained by Level set method.Whereas,the scheme with Level set method converges a cavity size shape that is quite similar to experimental photographs.

    Fig.9 Visual comparisons of the predicted cavities with observations from experiment

    Compared with the experimental data,the predicted cavity profiles at different models are examined in Fig.10.The numerical results vividly show significant benefit in the prediction using the free surface sharpening approach,the axial cavity width(Dc)with Level set method are well predicted.As we know,mixture model allows mutual penetration and mixing between different phases,there is large numerical dissipation in vacuoles length and radius of the direction,particularly reflected in the length direction.But,the Level set method includes complete interface features,and with the surface tension model,it is easy to handle diffusion and pulsation characteristics in numerical simulations.

    4 Conclusions

    In this present study,the Level set method for nature and ventilated cavitating flows are validated,various cases including the steady cavitation around the cylindrical body,unsteady cavitating flow around the Clark-Y hydrofoil and the ventilated cavitating around a disk are presented.Compared with the results obtained by mixture model and experimental visualizations:

    (1)For nature cavitating flows,although all three models give satisfactory predictions in overall pressure distributions and unsteady cavity shapes,differences are observed in the closure region of the cavity with the mixture model and Level set method,resulting from the differences in compressibility characteristics handled by each model.The cavity interface obtaind by Level set method is sharper.In addition,in the cavitation area,the vapor volume fraction is moderately larger than the results obtained with the mixture model.

    (2)For ventilated cavitating flows,the Level set method modified forms resulted in improved numerical accuracy and the ability to better resolve the cavity shape with the same cells.The scheme with Level set method converging a ventilated cavity size shape is very similar to experimental photographs and data,compared with the results obtained by mixture model.

    Fig.10 Plots of the cavity radius versus the axial distance

    [1]Hart C W,Nichols B D.Volume of fluid(VOF)method for the dynamics of free boundaries[J].J Comput Phys,1981,39(1):201-225.

    [2]Noh W F,Woodward P R.SLIC(simple line interface calculation)[C]//In:van de Vooren A,Zandbergen P(ends).Proceedings of Fifth International Conference on Fluid Dynamics,59 of Lecture Notes in Physics,Berlin:Springer,1976:330-340.

    [3]Ashgriz N,Poo J Y.FLAIR:flux line-segment model for advection and interface reconstruction[J].J compute Phys,1991,93(2):449-468.

    [4]Youngs D L.Time-dependent multi-material flow with large fluid distortion[M].In:Morton K W,Baines M J(eds).Numerical Method for Fluid Dynamics,New York:Academic Press,1982:273-285.

    [5]Ubbink O,Issa R I.A method for capturing sharp fluid interfaces on arbitrary meshes[J].J Comput Phys,1999,153(1):26-50.

    [6]Sussman M,Puckett E G.A coupled Level set and volume-of-fluid method for computating 3d and axisymmetric incompressible two-phase flows[J].Journal of Computational Physics,2000,162:301-337.

    [7]Sussman M,Smereka P,Osher S.A Level set approach for computating solutions to incompressible two-phase flow[J].Journal of Computational Physics,1994,114:146-159.

    [8]Shu Bitan,et al.Implementation of the Level Set Method into OpenFORM for capturing the free interface in incompressible fluid flows[J/OL].Documentation of OpenFOAM.www.openfoam.com.

    [9]Zwart P.Industrial CFD applications of free surface and cavitation flows[R].VKI Lecture Series:Industrial Two-Phase Flow CFD,2005.

    [10]Brackbill J U,Kothe D B,Zemach C.A continuum method for modeling surface tension[J].Journal of Computational Physics,1992,100:335-354.

    [11]Kubota A,Kato H,Yamaguchi H.A new modeling of cavitating flows:a numerical study of unsteady cavitation on a hydrofoil section[J].J Fluid Mech.,1992,240:59-96.

    [12]Menter F R.Zonal two equation k-ω turbulence models of aerodynamics flow[R].AIAA paper 1993-2906,1993.

    [13]Strelets M.Detached eddy simulation of massively separated flows[R].AIAA-01-0879,2001.

    [14]Rouse H,McNown J S.Cavitation and pressure distribution,head forms at zero angle of yaw[R].Studies in Engineering Bulletin 32,State University of lowa,1948.

    [15]Wang G Y,Senocak I,Shyy W.Dynamics of attached turbulent cavitating flows[J].Progress in Aerospace Sciences,2001,37(6):551-581.

    [16]Gopalan S,Katz J.Flow structure and modeling issues in the closure region of attached cavitation[J].Phys.of Fluids.2000,12:895-911.

    [17]Foeth E J,Vandoorne C W H,Vantersiega T,Wieneke B.Time resolved PIV and flow visualization of 3D sheet cavitation[J].Exp.in Fluids,2006,40:503-513.

    [18]Adrian R J.Twenty years of particle image velocimetry[J].Exp.in Fluids,2005,39:159-169.

    Level set方法在空化流動計算中的應(yīng)用

    黃 彪,王國玉,張敏弟,陳維清

    (北京理工大學(xué)機(jī)械與車輛工程學(xué)院,北京 100081)

    文章將Level set方法應(yīng)用于自然空化與通氣空化流動計算中,通過耦合入相間的質(zhì)量傳輸方程,建立了一個用于描述氣、汽、液多相空化流動的數(shù)值模擬方法,分別對繞圓頭回轉(zhuǎn)體和Clark-Y型水翼的自然空化流動和繞圓盤的通氣空化流動進(jìn)行了數(shù)值模擬研究,并與實(shí)驗(yàn)結(jié)果進(jìn)行了對比,研究結(jié)果表明:相較常用的多相流模型(mixture model),Level set方法通過時間和空間上的壓縮離散方案,減小了相間界面的擴(kuò)散,準(zhǔn)確地捕捉到了相間的界面,可以有效地運(yùn)用于空化多相流動。

    空化流動;混相模型;Level set方法

    TV131.32

    A

    黃 彪(1985-),男,北京理工大學(xué)機(jī)械與車輛工程學(xué)院博士研究生;

    陳維清(1985-),男,北京理工大學(xué)機(jī)械與車輛工程學(xué)院碩士研究生。

    TV131.32

    A

    1007-7294(2011)03-0207-10

    date:2010-07-24

    Supported by National Natural Science Foundation of China(Grant No.50679001 and Grant No.50979004)

    王國玉(1961-),男,博士,北京理工大學(xué)機(jī)械與車輛工程學(xué)院教授/博士生導(dǎo)師;

    張敏弟(1970-),女,博士,北京理工大學(xué)機(jī)械與車輛工程學(xué)院副教授;

    猜你喜歡
    北京理工大學(xué)相間空化
    功率超聲作用下鋼液中空化泡尺寸的演變特性
    鋼鐵釩鈦(2023年5期)2023-11-17 08:48:34
    北京理工大學(xué)機(jī)械與車輛學(xué)院簡介
    北京理工大學(xué)通信與網(wǎng)絡(luò)實(shí)驗(yàn)室
    新型分支相間導(dǎo)流排
    變壓器相間過流保護(hù)不正確動作原因的探討
    電子制作(2017年19期)2017-02-02 07:08:58
    三維扭曲水翼空化現(xiàn)象CFD模擬
    Design of Two-wheeled Mobile Control Robot with Holographic Projection
    不同運(yùn)動形式下水物相互作用空化數(shù)值模擬
    國家航天立法研討會在北京理工大學(xué)舉行
    太空探索(2015年1期)2015-07-18 11:02:13
    輸電線路相間距離保護(hù)應(yīng)用探討
    河南科技(2014年5期)2014-02-27 14:05:44
    最近2019中文字幕mv第一页| 国产一区二区三区综合在线观看| 亚洲美女搞黄在线观看| 亚洲精品一区蜜桃| 热99久久久久精品小说推荐| 欧美日本中文国产一区发布| 国产熟女欧美一区二区| 制服丝袜香蕉在线| 午夜福利在线观看免费完整高清在| 午夜福利,免费看| h视频一区二区三区| 女人被躁到高潮嗷嗷叫费观| 各种免费的搞黄视频| 自拍欧美九色日韩亚洲蝌蚪91| av在线播放精品| 丝袜脚勾引网站| 建设人人有责人人尽责人人享有的| 亚洲精品一区蜜桃| 国产免费一区二区三区四区乱码| 亚洲精品视频女| 寂寞人妻少妇视频99o| 免费观看a级毛片全部| 日韩中文字幕视频在线看片| 欧美国产精品一级二级三级| 97在线视频观看| 成人影院久久| 日本爱情动作片www.在线观看| 国产成人免费无遮挡视频| 在现免费观看毛片| 中文字幕制服av| 精品国产一区二区三区四区第35| 国产xxxxx性猛交| av视频免费观看在线观看| 男女无遮挡免费网站观看| 精品人妻一区二区三区麻豆| 一区福利在线观看| 熟女av电影| 看十八女毛片水多多多| 亚洲成人av在线免费| 国产精品av久久久久免费| 国产成人精品久久二区二区91 | 午夜精品国产一区二区电影| 久久精品国产a三级三级三级| 春色校园在线视频观看| 亚洲一级一片aⅴ在线观看| 亚洲人成电影观看| 黄色 视频免费看| 日本wwww免费看| 欧美日韩一级在线毛片| 国产精品免费大片| 卡戴珊不雅视频在线播放| 一级毛片我不卡| 亚洲精品成人av观看孕妇| 男人爽女人下面视频在线观看| 久久久精品国产亚洲av高清涩受| 成人漫画全彩无遮挡| 日本欧美国产在线视频| 亚洲美女搞黄在线观看| 黄色 视频免费看| 欧美精品高潮呻吟av久久| 国产不卡av网站在线观看| 国产乱人偷精品视频| 菩萨蛮人人尽说江南好唐韦庄| 熟女少妇亚洲综合色aaa.| 亚洲国产精品一区三区| 国产老妇伦熟女老妇高清| 亚洲精品国产av蜜桃| 久久这里有精品视频免费| av女优亚洲男人天堂| 狠狠精品人妻久久久久久综合| 秋霞伦理黄片| 国产有黄有色有爽视频| 国产精品偷伦视频观看了| 久久这里只有精品19| 激情五月婷婷亚洲| 日本wwww免费看| 美女国产视频在线观看| 国产精品免费大片| 两个人看的免费小视频| 色哟哟·www| 亚洲色图 男人天堂 中文字幕| 亚洲精品视频女| 国产精品秋霞免费鲁丝片| 1024视频免费在线观看| 少妇被粗大猛烈的视频| 欧美精品一区二区免费开放| 少妇人妻 视频| 亚洲经典国产精华液单| 十八禁网站网址无遮挡| 国产精品.久久久| 免费女性裸体啪啪无遮挡网站| 久久久国产精品麻豆| 丝袜美腿诱惑在线| 日本-黄色视频高清免费观看| 日产精品乱码卡一卡2卡三| 丝袜人妻中文字幕| 亚洲av中文av极速乱| 一级片免费观看大全| 91aial.com中文字幕在线观看| 人人妻人人澡人人看| 成年女人毛片免费观看观看9 | av视频免费观看在线观看| 蜜桃在线观看..| 尾随美女入室| h视频一区二区三区| 久久久久精品人妻al黑| 国产精品一区二区在线观看99| freevideosex欧美| 欧美日韩一级在线毛片| 2018国产大陆天天弄谢| 91精品国产国语对白视频| 不卡av一区二区三区| 精品99又大又爽又粗少妇毛片| 深夜精品福利| 中文字幕人妻丝袜一区二区 | 国产熟女欧美一区二区| 国产精品秋霞免费鲁丝片| 国产精品久久久久久av不卡| 日韩制服丝袜自拍偷拍| 久久久国产精品麻豆| 久久精品国产自在天天线| 国产人伦9x9x在线观看 | 免费观看在线日韩| 精品久久久久久电影网| 卡戴珊不雅视频在线播放| 午夜激情av网站| 18禁动态无遮挡网站| 国产精品欧美亚洲77777| 国产精品偷伦视频观看了| 大码成人一级视频| 一区二区三区乱码不卡18| 久久久久网色| av在线老鸭窝| 国产熟女欧美一区二区| 十分钟在线观看高清视频www| 久久精品久久精品一区二区三区| 亚洲三区欧美一区| 久久这里有精品视频免费| 久久99一区二区三区| 国产av一区二区精品久久| 亚洲欧美色中文字幕在线| 精品少妇黑人巨大在线播放| 亚洲成色77777| 日产精品乱码卡一卡2卡三| 色吧在线观看| 国产免费视频播放在线视频| 在现免费观看毛片| 超碰97精品在线观看| 极品少妇高潮喷水抽搐| 国产精品不卡视频一区二区| 麻豆av在线久日| 可以免费在线观看a视频的电影网站 | 色婷婷久久久亚洲欧美| 日韩中文字幕欧美一区二区 | 国产精品三级大全| 99热国产这里只有精品6| 国产女主播在线喷水免费视频网站| 七月丁香在线播放| 91精品伊人久久大香线蕉| 黄网站色视频无遮挡免费观看| 999久久久国产精品视频| 宅男免费午夜| 街头女战士在线观看网站| 天美传媒精品一区二区| 美女视频免费永久观看网站| 91精品三级在线观看| www.av在线官网国产| 亚洲成人av在线免费| 亚洲av成人精品一二三区| 亚洲精品av麻豆狂野| 最近最新中文字幕免费大全7| 久久人人爽av亚洲精品天堂| 欧美人与性动交α欧美软件| 免费黄网站久久成人精品| 亚洲精品一二三| 大陆偷拍与自拍| 日韩不卡一区二区三区视频在线| 色视频在线一区二区三区| 亚洲五月色婷婷综合| 日韩一区二区视频免费看| 电影成人av| 国产视频首页在线观看| 久久精品夜色国产| 欧美97在线视频| 人体艺术视频欧美日本| 亚洲av国产av综合av卡| 国产毛片在线视频| 亚洲三区欧美一区| 青春草国产在线视频| 亚洲精华国产精华液的使用体验| 免费观看在线日韩| 99热网站在线观看| 色哟哟·www| 啦啦啦啦在线视频资源| 日韩,欧美,国产一区二区三区| 成人国产麻豆网| 国产精品熟女久久久久浪| 免费看不卡的av| 免费人妻精品一区二区三区视频| 9色porny在线观看| 麻豆精品久久久久久蜜桃| 国产成人一区二区在线| 国产黄频视频在线观看| av在线播放精品| 久久久久久免费高清国产稀缺| 考比视频在线观看| 久久亚洲国产成人精品v| 免费观看在线日韩| 美女脱内裤让男人舔精品视频| 最新的欧美精品一区二区| 亚洲国产欧美在线一区| 在线 av 中文字幕| 久久精品国产鲁丝片午夜精品| 中文字幕av电影在线播放| 精品卡一卡二卡四卡免费| 伊人久久国产一区二区| 秋霞伦理黄片| 国产高清不卡午夜福利| 大片电影免费在线观看免费| 亚洲五月色婷婷综合| 亚洲男人天堂网一区| 哪个播放器可以免费观看大片| 好男人视频免费观看在线| 午夜影院在线不卡| 免费观看在线日韩| 国产精品免费大片| 99热全是精品| 国产毛片在线视频| 久久久久精品久久久久真实原创| 亚洲欧洲日产国产| 亚洲五月色婷婷综合| 熟女av电影| 欧美bdsm另类| 国产高清国产精品国产三级| 哪个播放器可以免费观看大片| av网站在线播放免费| 婷婷色综合大香蕉| 亚洲色图 男人天堂 中文字幕| 午夜日韩欧美国产| 黑人巨大精品欧美一区二区蜜桃| 伦理电影大哥的女人| 国产精品无大码| 777米奇影视久久| 下体分泌物呈黄色| 久久鲁丝午夜福利片| 少妇人妻久久综合中文| 欧美精品亚洲一区二区| 亚洲成av片中文字幕在线观看 | 国产视频首页在线观看| 如何舔出高潮| 国产一区二区激情短视频 | 日本欧美国产在线视频| 国产淫语在线视频| 人成视频在线观看免费观看| 亚洲一级一片aⅴ在线观看| 少妇人妻精品综合一区二区| 少妇的逼水好多| 亚洲av在线观看美女高潮| 亚洲国产欧美网| 在线亚洲精品国产二区图片欧美| 久久鲁丝午夜福利片| 成人手机av| 久久久精品94久久精品| 女人被躁到高潮嗷嗷叫费观| 中文字幕另类日韩欧美亚洲嫩草| 日韩伦理黄色片| 高清不卡的av网站| 欧美成人午夜免费资源| av有码第一页| 少妇人妻久久综合中文| 日韩欧美精品免费久久| 一个人免费看片子| 久久青草综合色| 人妻人人澡人人爽人人| 亚洲av.av天堂| 制服人妻中文乱码| 日日摸夜夜添夜夜爱| 国产成人91sexporn| 午夜免费鲁丝| 丝瓜视频免费看黄片| 人妻一区二区av| 亚洲美女黄色视频免费看| 啦啦啦在线免费观看视频4| 亚洲精品国产色婷婷电影| 男女国产视频网站| 九色亚洲精品在线播放| 精品一品国产午夜福利视频| 国产免费福利视频在线观看| 人人澡人人妻人| 高清在线视频一区二区三区| 午夜日本视频在线| 精品国产国语对白av| 考比视频在线观看| 久久久精品国产亚洲av高清涩受| 精品卡一卡二卡四卡免费| 考比视频在线观看| 国产欧美日韩综合在线一区二区| 欧美 亚洲 国产 日韩一| 国产精品不卡视频一区二区| 99热网站在线观看| 男人爽女人下面视频在线观看| 国产成人午夜福利电影在线观看| 亚洲精品,欧美精品| 久久国产亚洲av麻豆专区| 精品亚洲成a人片在线观看| 日韩,欧美,国产一区二区三区| 国产日韩欧美在线精品| 国产不卡av网站在线观看| 只有这里有精品99| 人人澡人人妻人| 国产乱来视频区| 岛国毛片在线播放| 婷婷色综合大香蕉| 婷婷色综合www| 精品视频人人做人人爽| 丝袜人妻中文字幕| 国产成人av激情在线播放| 91成人精品电影| 超色免费av| 色婷婷久久久亚洲欧美| 久久久久精品人妻al黑| 精品一区二区免费观看| 9热在线视频观看99| 老汉色av国产亚洲站长工具| 色哟哟·www| 乱人伦中国视频| 亚洲一区二区三区欧美精品| 丝袜美足系列| 18+在线观看网站| 欧美人与性动交α欧美精品济南到 | 色哟哟·www| 亚洲一区中文字幕在线| 久久精品aⅴ一区二区三区四区 | 精品国产乱码久久久久久小说| 精品国产露脸久久av麻豆| 中国三级夫妇交换| 黄色 视频免费看| 香蕉丝袜av| 亚洲综合色惰| 人妻系列 视频| 熟女电影av网| 亚洲色图综合在线观看| 亚洲少妇的诱惑av| 中文天堂在线官网| av在线播放精品| 欧美在线黄色| 国产国语露脸激情在线看| 欧美人与性动交α欧美软件| 一区二区日韩欧美中文字幕| 边亲边吃奶的免费视频| av在线老鸭窝| 久久国产精品大桥未久av| 男女边吃奶边做爰视频| 香蕉精品网在线| 成人亚洲欧美一区二区av| 97在线视频观看| 丝袜在线中文字幕| 夫妻性生交免费视频一级片| 亚洲内射少妇av| 亚洲国产欧美日韩在线播放| 日本-黄色视频高清免费观看| 久久久久久伊人网av| 国产成人aa在线观看| 亚洲在久久综合| 一本—道久久a久久精品蜜桃钙片| 建设人人有责人人尽责人人享有的| 在线天堂最新版资源| 一级毛片 在线播放| 一级a爱视频在线免费观看| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | av视频免费观看在线观看| 精品亚洲乱码少妇综合久久| 国产有黄有色有爽视频| 久久综合国产亚洲精品| 日韩电影二区| 免费在线观看黄色视频的| 美女国产高潮福利片在线看| 人人妻人人爽人人添夜夜欢视频| 久久久国产一区二区| 亚洲精品在线美女| 欧美xxⅹ黑人| 黑人猛操日本美女一级片| 高清在线视频一区二区三区| 欧美日韩精品成人综合77777| 午夜福利视频精品| 国产成人精品一,二区| 国产精品国产三级专区第一集| 久久午夜综合久久蜜桃| 纯流量卡能插随身wifi吗| 男的添女的下面高潮视频| 男男h啪啪无遮挡| 精品一品国产午夜福利视频| 菩萨蛮人人尽说江南好唐韦庄| 人妻系列 视频| www日本在线高清视频| 亚洲精品成人av观看孕妇| kizo精华| 免费高清在线观看视频在线观看| 婷婷色综合www| 欧美+日韩+精品| 成年女人毛片免费观看观看9 | 久久久精品免费免费高清| 日产精品乱码卡一卡2卡三| 我要看黄色一级片免费的| 日韩,欧美,国产一区二区三区| 亚洲内射少妇av| 熟妇人妻不卡中文字幕| 日韩av不卡免费在线播放| 免费观看无遮挡的男女| 免费看不卡的av| 黑人猛操日本美女一级片| 亚洲国产看品久久| 亚洲伊人久久精品综合| 麻豆精品久久久久久蜜桃| 亚洲第一区二区三区不卡| 999久久久国产精品视频| 免费观看在线日韩| 一区二区三区四区激情视频| 日韩一区二区视频免费看| 亚洲国产欧美在线一区| 国产精品麻豆人妻色哟哟久久| 日韩精品有码人妻一区| 欧美精品亚洲一区二区| 人人妻人人澡人人爽人人夜夜| 免费久久久久久久精品成人欧美视频| 国产深夜福利视频在线观看| 如日韩欧美国产精品一区二区三区| 久久久久精品性色| 国产一区二区激情短视频 | 久久午夜福利片| 国产毛片在线视频| 久久久久人妻精品一区果冻| 亚洲三级黄色毛片| 国产成人91sexporn| 国产一级毛片在线| 一级,二级,三级黄色视频| a级毛片黄视频| 国产有黄有色有爽视频| 人人妻人人添人人爽欧美一区卜| 国产野战对白在线观看| 国产国语露脸激情在线看| 亚洲欧洲日产国产| 夫妻午夜视频| 丝瓜视频免费看黄片| 中文欧美无线码| a级毛片黄视频| 日韩熟女老妇一区二区性免费视频| 免费不卡的大黄色大毛片视频在线观看| 天天躁夜夜躁狠狠躁躁| 一区在线观看完整版| 99热全是精品| 国产免费福利视频在线观看| 卡戴珊不雅视频在线播放| 亚洲国产看品久久| 久久久精品国产亚洲av高清涩受| 97精品久久久久久久久久精品| 在线天堂最新版资源| 一区二区三区四区激情视频| 国产熟女午夜一区二区三区| 蜜桃在线观看..| 七月丁香在线播放| 国产在视频线精品| 精品国产乱码久久久久久男人| 99国产综合亚洲精品| 国产免费一区二区三区四区乱码| 80岁老熟妇乱子伦牲交| 久久女婷五月综合色啪小说| 日韩不卡一区二区三区视频在线| 99九九在线精品视频| 成年动漫av网址| av在线观看视频网站免费| 最近中文字幕高清免费大全6| 女的被弄到高潮叫床怎么办| 伊人亚洲综合成人网| 精品少妇内射三级| 午夜影院在线不卡| 街头女战士在线观看网站| 亚洲成人手机| 久久久国产欧美日韩av| 欧美xxⅹ黑人| 欧美日韩亚洲国产一区二区在线观看 | 国产1区2区3区精品| 新久久久久国产一级毛片| 精品国产一区二区三区久久久樱花| av网站在线播放免费| 久久久久久久大尺度免费视频| 最近中文字幕2019免费版| 国产一区二区 视频在线| 久久久a久久爽久久v久久| 日本av免费视频播放| 不卡av一区二区三区| 黄色一级大片看看| 久久久精品94久久精品| 亚洲精品一区蜜桃| av.在线天堂| 亚洲欧美一区二区三区国产| 欧美日韩成人在线一区二区| 久久人人97超碰香蕉20202| 少妇人妻久久综合中文| 成人午夜精彩视频在线观看| 国产av国产精品国产| 国产有黄有色有爽视频| 亚洲精品成人av观看孕妇| 看非洲黑人一级黄片| 欧美精品高潮呻吟av久久| 亚洲欧洲精品一区二区精品久久久 | 女人高潮潮喷娇喘18禁视频| 丰满乱子伦码专区| 国产白丝娇喘喷水9色精品| 在线 av 中文字幕| 两性夫妻黄色片| 国产爽快片一区二区三区| 777米奇影视久久| 国产成人午夜福利电影在线观看| 美女福利国产在线| 国产精品 国内视频| 亚洲国产欧美日韩在线播放| 成人国语在线视频| 久久毛片免费看一区二区三区| 在线观看免费高清a一片| 国产成人精品无人区| 妹子高潮喷水视频| 少妇的丰满在线观看| 国产午夜精品一二区理论片| 国产免费视频播放在线视频| 亚洲欧美一区二区三区黑人 | 免费高清在线观看视频在线观看| 熟妇人妻不卡中文字幕| 国产激情久久老熟女| 欧美人与善性xxx| 性色avwww在线观看| 亚洲精品日韩在线中文字幕| 亚洲美女搞黄在线观看| 国产男女内射视频| 久久久久精品人妻al黑| 在线 av 中文字幕| 少妇被粗大猛烈的视频| 波多野结衣av一区二区av| 丰满迷人的少妇在线观看| 色网站视频免费| 秋霞在线观看毛片| av免费在线看不卡| 国产成人午夜福利电影在线观看| 9191精品国产免费久久| a级毛片在线看网站| www.自偷自拍.com| 免费av中文字幕在线| 黄色配什么色好看| a 毛片基地| 欧美bdsm另类| 日韩欧美一区视频在线观看| 一二三四中文在线观看免费高清| 成人毛片60女人毛片免费| 丝袜在线中文字幕| 777米奇影视久久| 九草在线视频观看| 下体分泌物呈黄色| 99热国产这里只有精品6| 男人添女人高潮全过程视频| 亚洲 欧美一区二区三区| 久久99一区二区三区| 国产成人欧美| 啦啦啦在线观看免费高清www| 日本黄色日本黄色录像| a级毛片黄视频| 男女免费视频国产| 又黄又粗又硬又大视频| 不卡视频在线观看欧美| 我的亚洲天堂| 精品亚洲成国产av| 亚洲经典国产精华液单| 国精品久久久久久国模美| 国产精品.久久久| 亚洲欧洲精品一区二区精品久久久 | 18在线观看网站| 亚洲国产色片| 天天操日日干夜夜撸| 亚洲欧美精品综合一区二区三区 | 观看美女的网站| 中文字幕人妻丝袜一区二区 | 久久久久久久精品精品| 女的被弄到高潮叫床怎么办| 一区二区三区精品91| a级毛片黄视频| 欧美激情高清一区二区三区 | 美女国产视频在线观看| 九九爱精品视频在线观看| 免费观看a级毛片全部| 97精品久久久久久久久久精品| 在线天堂最新版资源| 精品少妇久久久久久888优播| 国产精品国产av在线观看| 男女高潮啪啪啪动态图| 18+在线观看网站| 成人漫画全彩无遮挡| 精品人妻熟女毛片av久久网站| 日韩成人av中文字幕在线观看| 黄片小视频在线播放| 精品第一国产精品| 亚洲男人天堂网一区| 青青草视频在线视频观看| 黄色配什么色好看| 亚洲av综合色区一区| 丝袜美腿诱惑在线| av国产久精品久网站免费入址| 亚洲精品视频女| 精品亚洲成a人片在线观看| 成人漫画全彩无遮挡| 成人毛片60女人毛片免费| 国产在线一区二区三区精| 国产成人av激情在线播放| 欧美最新免费一区二区三区| 亚洲三区欧美一区| 成年女人毛片免费观看观看9 | 飞空精品影院首页|