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

    Two-phase SPH simulation of vertical water entry of a two-dimensional structure *

    2018-10-27 09:10:36LinMa馬林XiZhao趙曦HuaLiu劉樺
    關(guān)鍵詞:馬林

    Lin Ma (馬林), Xi Zhao (趙曦), Hua Liu (劉樺)

    1. Key Laboratory of Hydrodynamics (Ministry of Education), School of Naval Architecture, Ocean and civil Engineering, Shanghai Jiao Tong University, Shanghai, 200240 China

    2. State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

    3. College of Ocean Science and Engineering, Shanghai Maritime University, Shanghai 200135, China

    Abstract: In view of the fact that the SPH model is easy to handle the flows with the free surface of large deformation, a 2-D flow induced by vertical water entry of a 2-D structure is simulated using the two-phase SPH model.The local pressure of the boundary particles is obtained by pressure of the fluid particles nearby through a modified kernel approximation.To evaluate the accuracy of the method, water entry of a 2-D symmetric wedge with fixed separation point of the free surface on the wedge surface is simulated.The pressure distribution of the wedge at the initial stage agrees well with the analytical results available.Evolution of the free surface and the air flow in the cavity induced by the water entry are obtained.A higher speed air jet is found at the neck of the cavity when the neck of the cavity becomes smaller.For the case of a horizontal cylinder entering the water with an unknown separation point of flow on the model surface, the early stage of the water entry is simulated for the rigid body with different density.Evolution of the free surface deformation of the half-buoyant cylinder and neutrally buoyant cylinder water entry is compared with the experimental data.The effects of the density ratio and Froude number on the pinch-off of the cavity are discussed.It is found that the pinch-off time remains almost constant for different density ratio and Froude number.Meanwhile, for a given Froude number, the dimensionless pinch-off depth and the location of the cylinder at the time of pinch-off increase with the density ratio.Further, for a given density ratio, these two parameters increase with the Froude number and, however, the relative cavity shape appears to be a self-similar shape when Fr≥ 8 .35.

    Key words: Smoothed particle hydrodynamics (SPH), two-phase, water entry, density ratio

    Introduction

    Due to the complicated hydrodynamic features during the entry process, such as the deformation of free surface, evolution of the cavity region, fluid impact force and motion of the falling object, the water entry of a body has been considered as one of the challenge problems in hydrodynamics.May[1]systematically investigated the cavity-running behavior of vehicles and summarized the water entry problems.Recently, using the technology of highspeed cameras, Truscott et al.[2]designed a series of experiments to study water entry problems, including hydrophilic sphere, hydrophobic sphere, spinning ball and so on.Their experiments have shown miraculous underwater phenomenon after solid objects entering water.

    Since the interactions between the fluid and structure usually lead to the large free surface deformation and wave breaking, the numerical simulation of the violent flow with the free surface belongs to one of the complicated flows for computational fluid dynamics.For example, Zhao and Faltinsen[3]used the potential flow theory to solve vertical water entry of a 2-D axisymmetric body.Wu et al.[4]used the boundary element method to compute the hydrodynamic flows associated with the water entry of a wedge.The smoothed particle hydrodynamics (SPH) method is one of the mesh-free methods.It was firstly presented to solve astrophysical problems by Lucy[5].Besides the applications of the SPH method in multi-phase flows and heat transfer, it should be mentioned that the SPH method is of great advantages in simulating the large deformation and breaking of the free surface,e.g., Liu et al.[6]and Zhang et al.[7].Oger et al.[8]developed a SPH formulation with a spatially varying resolution, based on the so-called variable smoothing length technique, to capture the violent flow with the free surface induced by water entry of a 2-D wedge.Gonget al.[9-10]implemented the two phase SPH model to simulate water entry of a wedge to study the water entry problem of a wedge entering the free surface and validated with the experimental data for the late stage of the wedge running with a cavity.

    The robust treatment of solid boundary conditions is also quite important, especially for the boundaries surrounding the moving structure.The repulsive boundary, proposed by Monaghan[11], provided a simple treatment of the solid wall to prevent particle penetrations, but it could cause the pressure oscillations in the fluid region.Cummins and Rudman[12]developed the mirroring boundary to implement SPH simulation, which could be the most accurate one but it was achieved at the expense of the CPU time, and also the numerical program could become heavy for the complex boundary configurations.Gómez-Gesteira and Dalrymple[13]implemented the so-called dynamic boundary to provide an alternative way to treat the solid boundary with economic computing time.

    The evolution of the cavity generated by water entry of a body has been considered as the hydrodynamic flows which follow the similarity of the Froude number in the scale of engineering structures.The difference between the water density and the density of the structures could also influence the cavity evolution for the case of free dropping test of the vertical water entry.To describe the free surface profiles caused by a cylinder dropped into calm water,Greenhow and Lin[14]presented a sequence of photographs for water entry of a half-buoyant and neutrally buoyant cylinder.Semenov and Wu[15]presented a general similarity solution for water entry problems of a wedge with its inner angle fixed and its sides in expansion with flow detachment, in which the speed of expansion is a free parameter.Kiara et al.[16]applied SPH method to simulate the water entry of horizontal cylinders with vertical impact velocity,either kept constant or freely falling, without and with spin, into quiescent water under the effect of gravity.They presented the effects of non-dimensional submergence time, Froude numbers, spin ratios and mass ratios on the flow and cavity forming stages.However,further study on the buoyant effect on the cavity evolution is needed.

    In this paper, a new local pressure method is proposed which could resolve complex boundary or moving boundary conveniently.Details on the new local pressure method are discussed in later sections.A 2-D wedge water entry is simulated to validate the accuracy of the new local pressure method.The pressure distribution of the wedge at the initial stage agrees well with the analytical results.Furthermore, a study on vertical water entry of a two-dimensional cylinder with a moving curvature boundary based on the new local boundary is investigated.The first case is the test example of constant speed.The SPH results are compared with experimental results in literature.The second case is the impact of a horizontal free-motion cylinder which means the fluid-structure coupling should be considered.The horizontal cylinder has only one degree of freedom, namely the vertical motion.The drag force on the cylinder is obtained by the integration of the pressure of the particles on the cylinder.Both of the half-buoyant cylinder and neutrally buoyant cylinder are simulated.The evolution of the free surface and the depth of penetration computed by the SPH model are compared with experimental and other numerical results.Then,the effects of the density ratio and Froude number on the pinch-off time, pinch-off location and the shape of the cavity are discussed.

    1.SPH formulation

    1.1 Governing equations

    The flow is assumed to be non-viscous, namely ignoring the shear stress of the fluid.Both water flow and air flow are simulated.The governing equations are the Euler equations formulated as follows:

    whereρrefers to fluid density,trefers to time,urefers to the velocity vector of fluid particles,grefers to the gravitational acceleration,prefers to pressure, D/Dtrefers to the material derivative and? is the gradient operator.

    The SPH formulation of governing equations can be written as:

    whereWrefers to an interpolating kernel function,superscriptsaandbrefer to different particles,is an artificial viscosity term.In this paper, a spline function proposed by Monaghan and Lattanzio[17]is adopted:

    wheread=2/3h, 1 0/7 πh2, 1/πh3for 1-D, 2-D, 3-D condition,R=|r′-r|h,his the smoothing length.

    Fig.1 Water (black) and air (white) particles near the interface

    Two different particles are used in the water and air flow, namely water particles and air particles.The particles 3haway from the interface, are the same species in the support domain, as shown in Fig.1.When the particle is close to the interface enough, like particleb, both water particles and air particles are involved in the support domain.For particleb, the momentum equation can be written as

    where subscriptiandjrepresent water particle and air particle,N1andN2are the number of water particle and the number of air particle in the support domain of particleb.

    For classical SPH method, the pressure of the fluid particles is calculated through an equation of state as follows

    1.2 Boundary condition

    If a functionf(x) is sufficiently smooth in its definition domain, it can be expressed as Taylor series expansion in the vicinity ofxi

    Multiplying by a smoothing functionWboth sides of the above equation and integrating over the support domain, we can obtain the following equation

    If we replace the functionf(x) as the pressure function of boundary particles in Eq.(10), then

    which can be written as particle approximation as

    Based on the new local boundary treatment,problems with complex boundary or moving boundary can be solved.In order to validate the new approach, a 2-D cylinder water entry is simulated.

    2.Validation

    2.1 Water entry of a 2-D wedge

    In this section, a 2-D symmetric wedge water entry problem is simulated.The wedge is assumed to enter water vertically at a constant speed of 1 m/s.The deadrise angle of the wedge is 15°.The sketch of the simulation is shown in Fig.2.

    Fig.2 Sketch of vertical water entry of a 2-D wedge (β is the deadrise angle, v0 =1m/s)

    Fig.3 Pressure distributions on the wedge (t=0.004s)

    The pressure distribution of the wedge at the initial stage is shown in Fig.3, whereyrefers to the distance between the lowest point of the wedge and the initial free surface.The two-phase results have been compared with the analytical results obtained by Zhao et al.[3].

    The evolution of the cavity induced by the impact is shown in Fig.4.At the first stage, a splash appears in a very short time after the water entry.As the wedge descending, the water is forced outward along the wedge?s equator.The expansion of the water is resisted by hydrostatic pressure.Then the cavity begins to shrink.Eventually, the cavity collapses until the moment of deep seal.

    Fig.4 Evolution of the free surface

    The velocity field of the air flow in the cavity is shown in Fig.5.It is found that the speed of the air in the cavity?s neck becomes higher when the cavity begins to shrink.From Fig.6, the speed of the air flowing through the neck can be obtained as

    wheredis the width of the neck,Nairrefers to the number of the air particles in the cavity below the neck.

    Fig.5 Velocity distribution at t =0.236s

    The correlation between the dimensionless speed of the air flowing through the neck and the Froude number which is defined asshown in Fig.7.It is obvious that the air maintains low speed whenFr′≤13.WhenFr′≥ 26, the dimensionless speed increases rapidly and attains 182 whenFr′= 1 00.The details can refer to Ma and Liu[18].

    Fig.6 Sketch of the air flowing through the neck

    Fig.7 The development of the di me nsionless air spe ed a t the neck (un is the speed of the airflowing through the neck,v0 is the speed of the wedge and , d is the width of the neck)

    3.Numerical results and discussion

    3.1 Water entry of a horizontal cylinder

    In this section, classical problems of a horizontal circular cylinder water entry are solved by using the SPH model with the new boundary treatment presented previously.First of all, a test case with constant speed is simulated.The results are compared with the experimental data obtained by Miao[19].Actually,water entry problem with constant speed can be seen as a moving boundary problem.The drag force of the cylinder when it penetrates the water is obtained by the integration of the pressure of the boundary particles on the cylinder.Furthermore, a test case of free motion is simulated.In this condition, unlike constant speed case, fluid-structure coupling should be considered.The motion of the cylinder is governed by its gravity and the drag force.The cylinder has only one degree of freedom, which means only the vertical motion is allowed.

    Firstly, a horizontal circular cylinder is forced into an initially calm water surface with a constant vertical velocityv0.The SPH model is compared with the water entry tests conducted by Miao[19].The initial geometry of the computational domain is shown in Fig.8.

    Fig.8 Sketch of the horizontal cylinder water entry

    Fig.9 Non-reflection boundary treatment

    As shown in Fig.8, the water depth is 1 m and the diameter of the cylinder is 0.125 m.The cylinder enters into the water surface with a constant speed of 0.64 m/s.In our simulation, the computational domain could not be as large as the experiment of Miao,therefore, the length of the tank is set to be 1.6 m.In order to remove the influence of wall effects, a sponge layer is placed in the region within a certain distancedfrom the solid boundary to absorb the sound wave as shown in Fig.9[9].We use the non- reflection boundary treatment to prevent the sound wave from disturbing pressure field.In the sponge layer, the density of the fluid particles is obtained through

    where the coefficient of wave absorptionθtakes the form

    where

    wheredis the width of the sponge layer,ypis the coordinate of fluid particles,ywis the coordinate of the boundary of the sponge layer.

    The initial average distance of the particles is 0.0025 m, and the total number of the fluid particles is 400×640.The non-dimensional impact force is defined as

    whereFyis the total vertical hydrodynamic force reduced from the water entry including the buoyancy force,ρ, the speed of the cylinder is a constantv0=0.64m/s ,Dis the diameter of the cylinder.Fyis obtained by the integration of the pressure of particles on the cylinder

    Fig.10 Time history of slamming coefficient during water entry process ( t* =2 v0 t/D)

    Att=0, the cylinder is forced into the calm water surface.The slamming coefficientCsas a function of non-dimensional time is shown in Fig.10.The maximum impact pressure occurs at the instant when the cylinder impacting the free surface.From the figure, we can see the slamming coefficient decreases as the cylinder penetrates the water before 2v0t|D= 0.4(t= 0.039s).After 2v0t|D= 0.4, the slamming coefficient increases because the wet surface increases.At the time on 2v0t|D= 3.0(t=0.30s), the cylinder is totally wetted, and the slamming coefficient tending towards stability.The contribution of the buoyancy to the slamming coefficient is defined as

    where/(gD) .For this case in the paper,Csb=4 .71.As shown in Fig.10, the slamming coefficient obtained by SPH model agrees well with the experimental results.

    Secondly, the cylinder is considered as dropped into the calm water with free motion.Two kinds of cylinder will be simulated, one is half-buoyant, and the other is neutrally buoyant.The diameter of both cylinders is 0.11 m.The SPH simulation results are compared with the experiments conducted by Greenhow and Lin[14].

    For the case of constant speed, the cylinder is a moving boundary actually, but for the free motion,fluid-structure coupling should be considered.With the new local boundary method, the pressure on the cylinder could be obtained conveniently by Eq.(12).The drag force is obtained by the integration of pressure of the particles on the cylinder as Eq.(15).

    Fig.11 Evolution of the free surface deformation for the case of a half-buoyant cylinder and initial speed 3.13 m/s.Left panel: SPH simulations, Right panel: high-speed images by Greenhow and Lin[14]

    The simulation starts from the moment when the horizontal cylinder impacts the calm free surface.In the experiment of Greenhow and Lin, the cylinder is released from a height of 0.50 m between the lowest point of the horizontal cylinder and the calm free surface.The initial velocity of the cylinder is set toThe free fall time ist0≈v0|g≈0.139s when the lowest point of the cylinder impacts the free surface.

    Figure 11 illustrates the evolution of the free surface of a half-buoyant cylinder water entry.As shown in Fig.11(a), we can see two jets form from the contact points of the cylinder and the liquid.As one can see, the SPH model can deal with the problem of large deformation well.

    Figure 12 shows the time series of the penetration depth of the half-buoyant cylinder.The results of the SPH model are compared with the experimental data and computed values by the CIP method[20].Satisfactory agreements can be obtained for the half-buoyant cylinder.As shown in the figure, at the first stage whent<0.35s, the penetration depth is linearly proportional to the time because the drag force is basically equal to the gravity.As the cylinder descends, the wet surface increases.Finally, the drag force is greater than the gravity and the acceleration of the penetration depth decreases.

    Fig.12 Time history of penetration depth during a half-buoyant cylinder water entry

    The evolution of the free surface of a neutrally buoyant cylinder is shown in Fig.13.Similar to the half-buoyant cylinder, two jets form when the cylinder impacts the free surface.As the cylinder descends, an axisymmetric cavity appears below the free surface.The open cavity adjoins the cylinder near its equator,and it has a similar size with the cylinder.As the horizontal cylinder penetrates the water, its momentum is transferred to the fluid by forcing it outward radially, as shown in Fig.13(b).The expansion of the cavity is resisted by the hydro-static pressure and then the direction of the radial flow is reversed which makes the cavity begin to shrink as shown in Fig.13(c).Eventually, the cavity collapses below the surface which is called deep-seal.When the cavity collapses, one jet forms and moves up towards the free surface as shown in Fig.13(d).

    Figure 14 shows the time history of submergence depth of the neutrally buoyant cylinder.The experimental datum tagged by a question mark att=0.34s is questionable because this datum deviates from other data.As shown in the figure, the SPH results agree well with the CIP results and experimental data.

    Fig.13 Evolution of the free surface deformation for the case of a neutrally buoyant cylinder and initial speed 3.13 m/s.Left panel: SPH simulations, Right panel: high-speed images by Greenhow and Lin[14]

    3.2 Effect of the density ratio on pinch-off

    In this section, the effect of density ratio and Froude number on pinch-off is investigated numerically.The density ratio is defined asand the Froude number is defined as

    Fig.14 Time history of penetration depth during a neutrally buoyant cylinder water entry

    Fig.15 Dimensionless pinch-off time as a function of the Froude number for different density ratios ( r= ρc|ρwater)

    The depth of the water is 1.40 m and the diameter of the cylinder is 0.11 m.The horizontal cylinder is dropped into the water at initial speed at the range of 1-5 m/s.The Froude number values 0.93 to 23.20.The density ratio is set at the range of 1-6, say from the small density to the large density of the cylinder.The related parameters for the numerical experiments are summarized in Table1.In Fig.15, we plot data for the dimensionless pinch-off time as a function of the Froude number for different density ratio.It is obvious that the pinch-off time is independent of the Froude number for any density ratio.The pinch-off time is about 250 ms after impact.Figure 16 shows the dimensionless pinch-off locationh*and the depth of the center of the cylinder at the time of pinch-offH*as a function of the Froude number for different density ratio, whereh*=h|D,H*=H|D.In Fig.16, (a), (b) illustrate thath*andH*increase as the Froude number.For the same Froude number, a heavy cylinder could generate a big cavity in the process of water entry.In the limit of 0.93≤Fr≤23.20 and 1≤r≤6,h*andH*can scale linearly withForh*,c( 1)=0.80,ForH*,c(1)=Fig.16(c)shows thath|His almost a constant 0.3 for different density ratio.In summary, both Froude number and density ratio can alter the values ofh*andH*, but neither of them can alter the values ofh|Hevidently.

    Fig.16 The relationships between the characteristic coefficient and the Froude number.(a) The dimensionless pinchoff depth.(b) The dimensionless depth of the center of cylinder at the time of pinch-off.(c) The scaled pinchoff depth h|H

    Table1 Parameter values of water entry

    Fig.17 The water entry cavity shapes for different Froude numbers a nd different density ratios.(a) Computed cavity shape for different density ratios at Fr = 23.2 .(b)Computed cavity shape for different Froude number at r = 2 .(c) Computed cavity shape scaled with the length of the cavity for different Froude number at r= 2

    Computed cavity shapes are shown in Fig.17.The origin (x,y) =(0,0) corresponds to the location of the pinch-off point.In the first figure, the correlation between cavity shape and density ratio is shown in the condition ofFr= 23.2.Forr=1, the cavity shape is small because the speed of the cylinder decays sharply after the cylinder water entry.For bigger density ratio, the effect of the gravity could provide driving force to the cylinder, and makes the cavity shape longer and wider at the moment of pinch-off.In the second figure, the cavity shapes are scaled withτv0, in whichτrepresents pinch-off time.Five cavity shapes show different Froude number in the condition ofr=2.Based on this figure, the third figure is obtained.In the figure,ymaxis the length of the cavity.It is clearly shown that, the relative cavity shape is slimmer and forms to a more streamlined structure as the Froude number increases.Notably, the cavity shape appears to be a self-similar shape whenFr≥ 8 .35.

    4.Conclusions

    Vertical water entry of a 2-D structure is simulated based on the two-phase SPH model.A new local boundary treatment method in the framework of the coupling boundary method is presented.The method is derived by a corrective kernel approximation and validated by a 2-D symmetric wedge water entry.The pressure distribution of the wedge at the initial stage agrees well with the analytical results.Based on the results of 2-D wedge water entry, a high speed air jet is found when the neck area of the cavity becomes small.

    In order to verify the capability of the new method for complex boundary, horizontal cylinder water entry problems are simulated.First of all, test case of constant speed is calculated.In this condition,the cylinder is actually treated as a moving boundary.The evolution of the free surface and the slamming coefficient during cylinder water entry obtained through SPH model are compared with the experimental data.Secondly, test cases with free motion are simulated.In this situation, fluid-structure coupling is considered.The motion of the cylinder is controlled by the drag force and the gravity and the flow field is influenced by the cylinder water entry.A half-buoyant cylinder and a neutrally buoyant cylinder are simulated by the new method.It is found that the evolution of the free surface and the time history of penetration of depth agree well with the measured data available.

    The effect of the density ratio and Froude number on the pinch-off is studied.It is shown that the dimensionless pinch-off depth and the dimensionless submerged depth of the center of cylinder at the time of pinch-off increase as the density ratio for the same ratio increases, however the ratio of the pinch-off depth and the corresponding submerged depth remains a constant.For a given Froude number, larger density ratio could make the pinch-off depth and the corresponding submerged depth greater.The effects of Froude number and density ratio on the pinch-off depth and the submerged depth are obtained through numerical experiments.For a given density ratio, the relative cavity shape becomes slimmer as the Froude number increases and appears to be a self-similar shape for the cases of greater Froude number.

    猜你喜歡
    馬林
    芬蘭女總理的“熱舞風(fēng)波”
    藥檢證明芬蘭總理沒吸毒
    馬林:我又一次回來了
    羽毛球(2022年6期)2022-06-07 18:09:40
    電能計(jì)量裝置故障分析及維修技術(shù)
    馬林: 野蠻生長 真實(shí)怒放
    羽毛球(2021年4期)2021-05-07 02:57:12
    性感過頭?
    From cashier t o world's youngest prime minister
    From cashier to world
    桑娜·馬林, 一股溫柔的新活力
    桑娜·馬林,一股溫柔的新活力
    莫愁(2020年4期)2020-03-16 02:41:04
    免费看光身美女| 亚洲中文日韩欧美视频| 又黄又粗又硬又大视频| 波多野结衣高清无吗| 久久久成人免费电影| 18禁美女被吸乳视频| 在线免费观看不下载黄p国产 | 一本一本综合久久| 夜夜躁狠狠躁天天躁| 欧美日韩精品网址| 国产成人系列免费观看| 一个人免费在线观看电影| 一进一出抽搐动态| 日日摸夜夜添夜夜添小说| 亚洲 欧美 日韩 在线 免费| 国产精品av视频在线免费观看| 少妇裸体淫交视频免费看高清| 国产成人aa在线观看| 国产成人av激情在线播放| 国产亚洲精品久久久久久毛片| 男人舔女人下体高潮全视频| 香蕉丝袜av| x7x7x7水蜜桃| 9191精品国产免费久久| 国产av在哪里看| 国产精品三级大全| 日韩亚洲欧美综合| 可以在线观看毛片的网站| 国产黄片美女视频| 日韩成人在线观看一区二区三区| or卡值多少钱| 亚洲欧美一区二区三区黑人| 99久国产av精品| 午夜免费男女啪啪视频观看 | 午夜激情欧美在线| 色综合婷婷激情| 日韩免费av在线播放| 国产私拍福利视频在线观看| 久久九九热精品免费| 黄色丝袜av网址大全| 成人国产综合亚洲| 国产一区二区在线观看日韩 | 国产精品野战在线观看| 国产高潮美女av| 中文字幕人妻丝袜一区二区| 精品欧美国产一区二区三| 最近最新中文字幕大全电影3| 波野结衣二区三区在线 | 美女高潮的动态| 18+在线观看网站| 中文字幕av成人在线电影| 久久精品影院6| 亚洲成av人片免费观看| 国产成人av激情在线播放| 国产99白浆流出| 97超视频在线观看视频| 桃红色精品国产亚洲av| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 亚洲专区国产一区二区| 搡女人真爽免费视频火全软件 | 免费搜索国产男女视频| www国产在线视频色| 国产免费男女视频| 91av网一区二区| 男女之事视频高清在线观看| 色尼玛亚洲综合影院| 亚洲av二区三区四区| 国产主播在线观看一区二区| 精品久久久久久久末码| 男女那种视频在线观看| 日本免费一区二区三区高清不卡| 亚洲性夜色夜夜综合| 免费在线观看成人毛片| 天天一区二区日本电影三级| 他把我摸到了高潮在线观看| 在线观看日韩欧美| 国产亚洲精品一区二区www| 一进一出抽搐gif免费好疼| 免费高清视频大片| 三级毛片av免费| 亚洲中文日韩欧美视频| 国产不卡一卡二| 男女视频在线观看网站免费| 黑人欧美特级aaaaaa片| 国产极品精品免费视频能看的| 亚洲av免费在线观看| 亚洲最大成人手机在线| 中文字幕高清在线视频| 久久久久亚洲av毛片大全| 听说在线观看完整版免费高清| 亚洲一区高清亚洲精品| 三级男女做爰猛烈吃奶摸视频| 丰满的人妻完整版| 十八禁人妻一区二区| 观看免费一级毛片| 99久久久亚洲精品蜜臀av| 久久精品影院6| 亚洲人与动物交配视频| 99久久九九国产精品国产免费| avwww免费| 最近最新免费中文字幕在线| 久久人妻av系列| 不卡一级毛片| 久久这里只有精品中国| 国产黄a三级三级三级人| 69av精品久久久久久| 法律面前人人平等表现在哪些方面| 国产精品99久久99久久久不卡| 久久精品国产亚洲av香蕉五月| 亚洲精品久久国产高清桃花| 国产 一区 欧美 日韩| 久久久久国产精品人妻aⅴ院| 欧美乱色亚洲激情| h日本视频在线播放| ponron亚洲| 亚洲激情在线av| 制服丝袜大香蕉在线| 国产一区二区三区在线臀色熟女| 偷拍熟女少妇极品色| 中文字幕人妻丝袜一区二区| 久久婷婷人人爽人人干人人爱| 天堂动漫精品| 亚洲国产精品成人综合色| 久久久国产成人免费| 亚洲成a人片在线一区二区| 久久久久久久亚洲中文字幕 | 日韩国内少妇激情av| 免费在线观看影片大全网站| 一区福利在线观看| 亚洲av电影在线进入| 久久精品91无色码中文字幕| 99久久精品国产亚洲精品| 国内久久婷婷六月综合欲色啪| 国产成+人综合+亚洲专区| 精品一区二区三区视频在线观看免费| 国产成人系列免费观看| 69av精品久久久久久| 成人亚洲精品av一区二区| 欧美绝顶高潮抽搐喷水| 给我免费播放毛片高清在线观看| 免费观看人在逋| 亚洲国产精品成人综合色| 在线观看免费午夜福利视频| 亚洲人成网站在线播| 国产精品野战在线观看| 内射极品少妇av片p| 国内精品美女久久久久久| eeuss影院久久| 69人妻影院| 免费在线观看成人毛片| 日本一本二区三区精品| 波多野结衣巨乳人妻| 黄色日韩在线| 亚洲精品一区av在线观看| 在线免费观看的www视频| 美女高潮喷水抽搐中文字幕| 女生性感内裤真人,穿戴方法视频| 男女视频在线观看网站免费| 国产精品久久久久久久久免 | 变态另类丝袜制服| 俄罗斯特黄特色一大片| 日韩高清综合在线| 国产日本99.免费观看| 国产色婷婷99| h日本视频在线播放| 午夜免费成人在线视频| 19禁男女啪啪无遮挡网站| 日本 欧美在线| 成年免费大片在线观看| 国产探花在线观看一区二区| 国产亚洲精品久久久久久毛片| 国产不卡一卡二| 免费看美女性在线毛片视频| 51国产日韩欧美| 两个人看的免费小视频| 欧美不卡视频在线免费观看| 久久欧美精品欧美久久欧美| www.999成人在线观看| 国产精品久久视频播放| 亚洲国产高清在线一区二区三| 99精品欧美一区二区三区四区| 国产视频内射| 欧美最新免费一区二区三区 | 久久久久久久久大av| 好男人在线观看高清免费视频| 99久久精品一区二区三区| 3wmmmm亚洲av在线观看| 禁无遮挡网站| www日本在线高清视频| 国产蜜桃级精品一区二区三区| 麻豆国产av国片精品| 又黄又粗又硬又大视频| 老汉色av国产亚洲站长工具| 欧美日韩亚洲国产一区二区在线观看| 成人国产一区最新在线观看| 美女cb高潮喷水在线观看| 中文资源天堂在线| 波多野结衣高清作品| www国产在线视频色| 成年女人看的毛片在线观看| 欧美成人性av电影在线观看| 国产伦精品一区二区三区四那| 午夜福利欧美成人| 国产美女午夜福利| 三级男女做爰猛烈吃奶摸视频| 国产精品一区二区三区四区久久| 一卡2卡三卡四卡精品乱码亚洲| 我的老师免费观看完整版| 男女之事视频高清在线观看| 99久久精品一区二区三区| 国产真实伦视频高清在线观看 | 精品久久久久久久毛片微露脸| 久久伊人香网站| 精品人妻偷拍中文字幕| www.www免费av| 亚洲一区二区三区不卡视频| 国产av不卡久久| 欧美成人免费av一区二区三区| 日韩中文字幕欧美一区二区| 一夜夜www| 国产成人a区在线观看| 99国产极品粉嫩在线观看| 欧美黑人巨大hd| 午夜免费男女啪啪视频观看 | 丰满人妻熟妇乱又伦精品不卡| 国产午夜福利久久久久久| 我要搜黄色片| 少妇的逼水好多| 婷婷六月久久综合丁香| 国产精品亚洲av一区麻豆| 动漫黄色视频在线观看| 日本免费a在线| 看免费av毛片| 一进一出好大好爽视频| 欧美成人免费av一区二区三区| 久久久久久人人人人人| 午夜老司机福利剧场| 久久精品国产综合久久久| 亚洲激情在线av| 亚洲一区高清亚洲精品| 日韩欧美国产一区二区入口| 嫁个100分男人电影在线观看| 国产午夜精品论理片| 亚洲激情在线av| 又爽又黄无遮挡网站| av在线蜜桃| 两个人视频免费观看高清| 欧美最黄视频在线播放免费| 高清日韩中文字幕在线| 国产97色在线日韩免费| 大型黄色视频在线免费观看| 亚洲av免费高清在线观看| 在线十欧美十亚洲十日本专区| 国产精品三级大全| 国产 一区 欧美 日韩| 亚洲精品成人久久久久久| av女优亚洲男人天堂| 婷婷亚洲欧美| 我要搜黄色片| 国产高清videossex| 国产精品99久久久久久久久| 日本撒尿小便嘘嘘汇集6| 成人一区二区视频在线观看| 国产爱豆传媒在线观看| 超碰av人人做人人爽久久 | 亚洲在线观看片| 成人高潮视频无遮挡免费网站| 精品人妻1区二区| 欧美一区二区亚洲| 成年女人毛片免费观看观看9| 黄色丝袜av网址大全| 国产精品一区二区三区四区免费观看 | 免费在线观看亚洲国产| 国产免费男女视频| 美女大奶头视频| 日韩欧美免费精品| 免费看日本二区| 757午夜福利合集在线观看| 亚洲美女视频黄频| 十八禁人妻一区二区| 国产在线精品亚洲第一网站| 不卡一级毛片| 长腿黑丝高跟| 少妇人妻精品综合一区二区 | netflix在线观看网站| 久久精品国产清高在天天线| 久久伊人香网站| 欧美av亚洲av综合av国产av| x7x7x7水蜜桃| 成人av一区二区三区在线看| 亚洲av第一区精品v没综合| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 精品国产三级普通话版| 热99re8久久精品国产| 亚洲一区二区三区不卡视频| 床上黄色一级片| 人妻久久中文字幕网| 免费无遮挡裸体视频| 亚洲国产欧美人成| 免费在线观看日本一区| 亚洲18禁久久av| 国产成人啪精品午夜网站| 久久国产乱子伦精品免费另类| 亚洲成人久久爱视频| 欧美三级亚洲精品| 十八禁人妻一区二区| 国产精品一区二区三区四区免费观看 | 亚洲精品乱码久久久v下载方式 | 天堂网av新在线| 国产不卡一卡二| 美女高潮的动态| 亚洲 国产 在线| 免费观看的影片在线观看| 99久久九九国产精品国产免费| 俄罗斯特黄特色一大片| 真人做人爱边吃奶动态| 国产黄a三级三级三级人| 国产高清三级在线| 久久久久久久久久黄片| 高清毛片免费观看视频网站| 久久精品亚洲精品国产色婷小说| 免费大片18禁| 国产伦一二天堂av在线观看| 听说在线观看完整版免费高清| 可以在线观看毛片的网站| 亚洲av美国av| 真实男女啪啪啪动态图| 免费电影在线观看免费观看| 国产在视频线在精品| 国产一区二区在线观看日韩 | 亚洲专区中文字幕在线| 狂野欧美白嫩少妇大欣赏| 麻豆久久精品国产亚洲av| 99riav亚洲国产免费| 老司机在亚洲福利影院| 啪啪无遮挡十八禁网站| 国产黄片美女视频| 午夜福利在线观看吧| 老熟妇乱子伦视频在线观看| 久久久久久久精品吃奶| 成人国产综合亚洲| 欧美成人a在线观看| 色老头精品视频在线观看| 18禁在线播放成人免费| 无限看片的www在线观看| 乱人视频在线观看| netflix在线观看网站| 欧美性猛交黑人性爽| 夜夜躁狠狠躁天天躁| 久久这里只有精品中国| 久久久久精品国产欧美久久久| 日韩 欧美 亚洲 中文字幕| 最新在线观看一区二区三区| 亚洲美女黄片视频| 首页视频小说图片口味搜索| 久久久久亚洲av毛片大全| 热99re8久久精品国产| 精品人妻偷拍中文字幕| 中文字幕精品亚洲无线码一区| 不卡一级毛片| 日韩有码中文字幕| 日韩精品青青久久久久久| 别揉我奶头~嗯~啊~动态视频| 午夜福利在线观看免费完整高清在 | 18禁裸乳无遮挡免费网站照片| 黄色日韩在线| 久久亚洲真实| 国产老妇女一区| 午夜福利在线观看吧| av天堂中文字幕网| 亚洲精品久久国产高清桃花| 亚洲精品在线美女| 国产真实乱freesex| 精品一区二区三区人妻视频| 麻豆成人午夜福利视频| 嫩草影视91久久| 三级毛片av免费| 日韩精品青青久久久久久| 美女 人体艺术 gogo| 嫩草影视91久久| 国产中年淑女户外野战色| 神马国产精品三级电影在线观看| 色视频www国产| 国产男靠女视频免费网站| 老司机午夜十八禁免费视频| 在线视频色国产色| 亚洲久久久久久中文字幕| 国产高清激情床上av| 丰满人妻熟妇乱又伦精品不卡| www.色视频.com| 国内久久婷婷六月综合欲色啪| 九九久久精品国产亚洲av麻豆| 亚洲精品乱码久久久v下载方式 | 国产亚洲精品综合一区在线观看| 少妇裸体淫交视频免费看高清| 搞女人的毛片| 欧美黑人巨大hd| 又紧又爽又黄一区二区| 亚洲av免费高清在线观看| 三级国产精品欧美在线观看| 日韩有码中文字幕| 男女午夜视频在线观看| 蜜桃亚洲精品一区二区三区| 国产主播在线观看一区二区| 中文字幕人妻熟人妻熟丝袜美 | 欧美日本视频| 国产伦人伦偷精品视频| 激情在线观看视频在线高清| 国模一区二区三区四区视频| 日本黄大片高清| 国产淫片久久久久久久久 | 国产一区在线观看成人免费| 国产精品久久久久久久久免 | 中文字幕人妻丝袜一区二区| 一级毛片高清免费大全| 超碰av人人做人人爽久久 | 露出奶头的视频| 91在线观看av| 亚洲七黄色美女视频| 欧美不卡视频在线免费观看| 午夜免费激情av| 欧美三级亚洲精品| 狂野欧美白嫩少妇大欣赏| 女人高潮潮喷娇喘18禁视频| 老司机福利观看| 九九久久精品国产亚洲av麻豆| 一本精品99久久精品77| 免费高清视频大片| 欧美日本视频| 色播亚洲综合网| 国产亚洲精品久久久久久毛片| 一本一本综合久久| 国产精品影院久久| 婷婷丁香在线五月| 国产精品免费一区二区三区在线| 国产av麻豆久久久久久久| 成人特级av手机在线观看| 久久久久免费精品人妻一区二区| 99国产精品一区二区蜜桃av| 欧美区成人在线视频| a在线观看视频网站| 亚洲欧美日韩东京热| 欧美性感艳星| 国产综合懂色| av国产免费在线观看| 精品日产1卡2卡| av国产免费在线观看| 女警被强在线播放| 久久精品亚洲精品国产色婷小说| 中文字幕人妻丝袜一区二区| www.www免费av| 麻豆一二三区av精品| 丰满人妻一区二区三区视频av | 国产精品乱码一区二三区的特点| 欧美一区二区国产精品久久精品| 久久国产精品人妻蜜桃| 国产高清激情床上av| 亚洲一区高清亚洲精品| 午夜免费成人在线视频| 欧美黄色淫秽网站| www.色视频.com| 啦啦啦观看免费观看视频高清| 综合色av麻豆| 757午夜福利合集在线观看| 亚洲精品久久国产高清桃花| 精品人妻偷拍中文字幕| 99久久综合精品五月天人人| 日韩欧美三级三区| 亚洲成av人片在线播放无| 久久6这里有精品| 精品不卡国产一区二区三区| 一进一出抽搐动态| 3wmmmm亚洲av在线观看| 国产成人a区在线观看| a级毛片a级免费在线| www日本在线高清视频| 99riav亚洲国产免费| 国产极品精品免费视频能看的| 夜夜爽天天搞| 欧美不卡视频在线免费观看| 久久久久久久久大av| 老熟妇仑乱视频hdxx| 18禁美女被吸乳视频| 亚洲国产精品合色在线| 久久亚洲真实| 三级男女做爰猛烈吃奶摸视频| 99riav亚洲国产免费| 在线天堂最新版资源| 毛片女人毛片| 熟女电影av网| www.www免费av| 一级黄色大片毛片| 国产成人影院久久av| av天堂中文字幕网| 日日夜夜操网爽| 丰满人妻一区二区三区视频av | 色老头精品视频在线观看| 18禁美女被吸乳视频| 99国产精品一区二区三区| 偷拍熟女少妇极品色| 国产精品99久久久久久久久| 免费大片18禁| 欧美+日韩+精品| 国产精品精品国产色婷婷| 免费在线观看日本一区| 久久中文看片网| 精品久久久久久久久久久久久| av国产免费在线观看| 欧美中文日本在线观看视频| 女人高潮潮喷娇喘18禁视频| 亚洲人成伊人成综合网2020| 欧美一区二区精品小视频在线| 久久九九热精品免费| 国产黄片美女视频| 日本a在线网址| 久久久久免费精品人妻一区二区| a级毛片a级免费在线| 看免费av毛片| 变态另类丝袜制服| 国产激情偷乱视频一区二区| 国产黄a三级三级三级人| 亚洲av中文字字幕乱码综合| 亚洲无线在线观看| 少妇的逼水好多| 久久久精品欧美日韩精品| 搡老熟女国产l中国老女人| 制服丝袜大香蕉在线| 久久久久国产精品人妻aⅴ院| 日韩av在线大香蕉| av天堂中文字幕网| 色精品久久人妻99蜜桃| 在线观看66精品国产| 亚洲五月天丁香| 小说图片视频综合网站| 欧美乱妇无乱码| 日韩免费av在线播放| 国产中年淑女户外野战色| www国产在线视频色| 悠悠久久av| 国产真实乱freesex| 身体一侧抽搐| 欧美xxxx黑人xx丫x性爽| 成年版毛片免费区| 夜夜爽天天搞| 少妇高潮的动态图| 啪啪无遮挡十八禁网站| 国产视频一区二区在线看| 18禁裸乳无遮挡免费网站照片| 天天一区二区日本电影三级| 狠狠狠狠99中文字幕| 欧美日本视频| 我要搜黄色片| 岛国在线观看网站| 欧美在线一区亚洲| 99热这里只有精品一区| 99精品在免费线老司机午夜| 亚洲熟妇熟女久久| 国产精品电影一区二区三区| 淫秽高清视频在线观看| 小蜜桃在线观看免费完整版高清| 欧美日韩中文字幕国产精品一区二区三区| 国产探花在线观看一区二区| 久99久视频精品免费| 成人性生交大片免费视频hd| 一区二区三区高清视频在线| 久久久国产成人免费| 国产精品综合久久久久久久免费| 久久久国产成人免费| 色哟哟哟哟哟哟| 女警被强在线播放| 村上凉子中文字幕在线| 校园春色视频在线观看| 香蕉av资源在线| 亚洲欧美日韩高清在线视频| 亚洲国产欧美人成| 一本综合久久免费| 欧美在线一区亚洲| 国产精品1区2区在线观看.| 亚洲精华国产精华精| 日本黄大片高清| 午夜激情福利司机影院| 日本五十路高清| 成年女人毛片免费观看观看9| 99久久精品热视频| 精品人妻一区二区三区麻豆 | 亚洲va日本ⅴa欧美va伊人久久| 免费搜索国产男女视频| 国产高潮美女av| 国内久久婷婷六月综合欲色啪| 深爱激情五月婷婷| 欧美日韩黄片免| 欧美成人免费av一区二区三区| 久久精品综合一区二区三区| 国产黄片美女视频| 色视频www国产| 亚洲内射少妇av| 婷婷精品国产亚洲av| 国产视频内射| 亚洲精品美女久久久久99蜜臀| 在线十欧美十亚洲十日本专区| 黄色片一级片一级黄色片| 国产精品久久久久久精品电影| 91九色精品人成在线观看| 最近在线观看免费完整版| 一区福利在线观看| 国产亚洲精品综合一区在线观看| 欧美乱妇无乱码| 国产高清激情床上av| 久久精品国产综合久久久| 1024手机看黄色片| 精品国产三级普通话版| 国产午夜精品论理片| 成人av在线播放网站| 热99在线观看视频| 国产成+人综合+亚洲专区|