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

    Reverse motion characteristics of water-vapor mixture in supercavitating flow around a hydrofoil*

    2016-12-06 08:15:22XiangbinLI李向賓NanLI李楠GuoyuWANG王國(guó)玉MindiZHANG張敏弟
    關(guān)鍵詞:李楠王國(guó)

    Xiang-bin LI (李向賓), Nan LI (李楠), Guo-yu WANG (王國(guó)玉), Min-di ZHANG (張敏弟)

    1. School of Nuclear Science and Engineering, Beijing Key Laboratory of Passive Safety Technology for Nuclear Energy, North China Electric Power University, Beijing 102206, China, E-mail: lixiangbin@ncepu.edu.cn

    2. School of Mechanics and Vehicle Engineering, Beijing Institute of Technology, Beijing 100081, China

    Reverse motion characteristics of water-vapor mixture in supercavitating flow around a hydrofoil*

    Xiang-bin LI (李向賓)1, Nan LI (李楠)1, Guo-yu WANG (王國(guó)玉)2, Min-di ZHANG (張敏弟)2

    1. School of Nuclear Science and Engineering, Beijing Key Laboratory of Passive Safety Technology for Nuclear Energy, North China Electric Power University, Beijing 102206, China, E-mail: lixiangbin@ncepu.edu.cn

    2. School of Mechanics and Vehicle Engineering, Beijing Institute of Technology, Beijing 100081, China

    The supercavitation has attracted a growing interest because of its potential for high-speed vehicle maneuvering and drag reduction. To better understand the reverse flow characteristics of a water-vapor mixture in supercavitating flows around a hydrofoil,a numerical simulation is conducted using a unified supercavitation model, which combines a modified RNG -kε turbulence model and a cavitation one. By comparing the related experimental results, the reverse motion of the water-vapor mixture is found in the cavitation area in all supercavitation stages. The inverse pressure gradient leads to reverse pressure fluctuations in the cavity,followed by the reverse motion of the water-vapor two-phase interface. Compared with the water-vapor mixture area at the back of the cavity, the pressure in the vapor area is inversely and slowly reduced,a higher-pressure gradient occurs near the cavity boundary.

    supercavitation, water-vapor interface, reverse motion, pressure distribution, numerical simulation

    Introduction

    The cavitation flows are very complicated because of the large density variations in the flow field, the frequent pressure fluctuations, and the violent watervapor interface mass transfer[1]. They can be divided into four stages based on the special cavitation number: the inception cavitation, the sheet cavitation, the cloud cavitation, and the supercavitation. As a fully developed pattern, the supercavitation comes into being when the cavity covers the entire hydrofoil suction surface and gradually extends downstream[2]. The pressure inside a supercavitation area is often very low, where a relatively stable cavity forms. Thus, a free-stream method based on the potential flow theory was developed to predict the supercavitation dynamics, providing a linear and slender-body theory for two-dimensional supercavitating flows. However, in this method, it is assumed that the pressure distribution in a supercavitation area is uniform and equal to the local vaporization pressure, and that the cavity interface is streamlined and smooth.

    In recent years, the supercavitation has attracted growing interests because of its potential for highspeed vehicle maneuvering and drag reduction[3-8]. Since 2000, numerous advanced measurement methods, equipment, and computation technologies were adopted to study supercavitating flows. Wosnik et al.[9]measured the ventilated supercavitation in turbulent bubbly wakes with a high speed camera and the particle image velocimetry, Shafaghat et al.[10]used a socalled “non-dominated sorting genetic algorithm” to carry out the shape optimization of two-dimensional cavitators in supercavitating flows, Rashidi et al. studied the ventilated supercavitation by using two methods: the VOF method and the Young's algorithm[11]. As for natural supercavitating flows around a hydrofoil, Wang et al.[2]observed that the water and vapor bubbles filled the cavitation region in the supercavitation stage. Further investigations indicate that the supercavitating flows are identified in the following three features: (1) the fluctuation cavity with periodic vortex shedding, (2) the vapor and watervapor mixture coexisting inside the cavity with a tur-bulent wake, and (3) a cavity largely filled with vapor having a two-phase tail and distinct phase boundaries in the wake region. Substantial unsteady features exist although the cavity boundary is visually smooth and the water-vapor mixture moves reversely over time[12]. In addition, the velocity in the supercavitating region has a special distribution even in the core of the cavitation region, and the lower velocity area is reduced as the cavitation number decreases.

    To further understand the reverse motion mechanism of the water-vapor mixture in supercavitating flows around a hydrofoil, numerical simulations are carried out in this study, focusing on the unsteady characteristics, including the transient cavitation structures and the pressure distribution. The commercial software CFX is selected and the computational model is based on a unified supercavitation one, which combines a modified RNG -kε turbulence model and a cavitation one.

    Fig.1 Schematic diagram of the hydrofoil

    Fig.2 Sketch of the foil position in the test section

    Fig.3 Grid distribution around the foil

    1. Supercavitation model

    1.1Hydrofoil, grid, boundary and initial conditions

    A supercavitation hydrofoil investigated by Tulin[13], as shown in Fig.1, is adopted in the present study. Figure 2 gives the hydrofoil position in the test section, of 10c in length, 2.7c in height, and 1c in width (where c is the chord length of the hydrofoil, more details can be found in Ref.[12]). A two-dimenonal computational domain is chosen to correctly follow the geometry of the test section. An orthogonal mesh with 62 720 cells is generated as shown in Fig.3. Alternative mesh systems, with 41 540 and 85 580 cells, are also chosen to evaluate the mesh convergence of the computation. To ensure more accurate results, a finer mesh is used around the hydrofoil and the cavitation region as the related experiments indicate. The corresponding evaluation was carried out in the Ref.[14], and the results show that a similar pressure distribution is obtained under three kinds of mesh conditions. With both accuracy and economy in mind, the mesh of medium size, 62 720 cells, is selected for the following simulation.

    In this study, the foil lift angle θ is fixed at o 15,the reference velocity U is fixed at 10 m/s, the Reynolds number and the cavitation number are defined as follows:

    where ν is the kinematic viscosity of the liquid,vp represents the Phase-change threshold pressure of the vapor, p∞is the reference static pressure, andlρ is the density of the liquid.

    In addition, the numerical simulation conditions are as follows: the inflow velocity, the volume fractions and the turbulence quantities are specified at the inlet boundary, and the cross-sectional-averaged static pressure is imposed as the reference pressure at the outlet. A no-slip boundary condition is used at both the upper and lower walls, and the no-slip wall condition with a wall function is provided at the foil surface.

    1.2Governing equations

    A unified supercavitation model is based on the governing equations of both a turbulence model and a mass transport one. First, based on the Favre average form of the mass conservation equation and the Navier-Stokes equation for the momentum conservation, the RNG -kε turbulence model is as follows:

    Fig.4 Instantaneous cavity structures of experiments and CFD simulations

    where ρ is the density of the liquid-vapor mixture, k is the turbulent kinetic energy, ε is the turbulent dissipation rate, u is the velocity component,kp is the source term representing the production of k, the empirical constants σk=σε=0.7179, Cε1=0.085, Cε2=1.680, μ and μtrepresent the laminar viscosity and the turbulent viscosity, respectively. The turbulent viscosity is defined as follows

    where the empirical constant Cμ=0.09. It is found that the standard RNG -kε model will result in shorter cavity lengths as compared with previously reported experimental data[15]. Accordingly, a modified RNG -kε turbulence closure model is adopted in this study. Instead of ()=fρρ in the original model,the function is defined as follows

    wherevρ is the density of the liquid and the vapor,vα is the vapor volume fraction, and n is the empirical constant. The effect of the vapor phase is added.

    Secondly, considering that the cavitation process is dominated by the thermodynamics and the kinetics of the phase-change dynamics in the flow field, a mass transport model (i.e., the cavitation model) is provided based on the vapor volume fraction, just as shown in Eq.(7). The source terms m˙+and m˙-represent the evaporation and the condensation of the phases, respectively[16].

    Fig.5 Instantaneous distribution of vapor volume fraction

    Here, Γ is the diffusion coefficient, the empirical constants=5×10-4, F=50 and F=0.01, thee c vapor bubble radius RB=10-6m , p is the pressure in the vapor bubble. Several experimental studies have shown the significant effect of turbulence on cavitating flows. The local turbulent pressure fluctuationshould be considered in the definition of the vapor pressure, thus a modified model is adopted as follows:

    Fig.6 Instantaneous flow visualization of experiments and CFD simulations

    whereturbp is the local turbulent fluctuating pressure andsatp is the saturated vapor pressure. This approach is simple and yields better outcomes[17].

    The abovementioned modifications on the turbulence and cavitation models are added into CFX software by means of self-programming. Based on the unified supercavitation model, the standard RNG -kε model and the modified one are used to simulate and compare the data with the experimental ones. Just as shown in Fig.4, the cavity length simulated by the modified model (with =2n) agrees well with the experimental one. So, the empirical constant n is selected as =2n[18].

    2. Results and discussions

    Based on the above numerical model, some typical supercavitating flows around a hydrofoil are simulated in different supercavitation stages. The detailed flow structures and the pressure distribution, especially the reverse motion mechanism, are analyzed as follows.

    2.1Flow structures

    Figure 5 shows a typical cavity structure contour,where the void fraction in a cavity position is related with the vapor volume fraction. The high content vapor(the vapor volume fraction above 0.8, which is called the pure vapor-phase area) is full in the front cavity,whereas the water-vapor mixture phase occupies the rear one. The section between the pure vapor-phase area and the water-vapor mixture area is defined as the water-vapor two-phase interface.

    Figure 6 shows the transient cavitation structures,where three typical cavitation numbers such as 0.74,0.67, and 0.54 are selected. The flows with σ=0.74 corresponds to the first supercavitating stage: the fluctuation cavity with periodic vortex shedding, the flows with σ=0.54indicates the second supercavitating stage: the vapor and the water-vapor mixture coexist inside the cavity with a turbulent wake, the flows with σ=0.67 belongs to the middle developing stage[12]. The calculated results agree well with the experimental data[16]. Figure 7 gives the corresponding relative cavity length of both the experimental results and the numerical ones for different cavitation numbers. The relative cavity length, l*, is defined as the ratio of the actual cavity length to that of the hydrofoil chord. With a deviation of no more than 11%, these data also confirm the above conclusion. Here, the actual cavity length is defined as the horizontal one beginning from the foil's leading point (point A, as shown in Fig.1) to the cavity closure point. For the numerical length, it can be measured from the CFX chart easily. For the experimental one, it is obtained by means of a home-made software with an error bar of approximately 3%-10%. That error is mainly resulted from the uncertainty of judgment on the cavity closure position.

    Fig.7 Relative cavity length at different cavitation numbers

    Fig.8 Reverse movement cycle of the two-phase interface at different cavitation numbers

    Moreover, a reverse motion of the water-vapor mixture phase occurs within the cavity in all stages. The results show that the frequent exchange between the water phase and the vapor phase still exists inside the cavity, despite the fact that the cavity is already enclosed and that the cavitation contour line is relatively stable. The motion of the interface between the vapor phase and the water-vapor mixture is seen as the variation of the two-phase volume fraction inside the cavity. The cavity length is smaller with σ=0.74and σ=0.67. The water-vapor two-phase interface reverses to the hydrofoil end from the back of the cavity with a relatively short period, whereas the process lasts up to 92 ms until σ=0.54. For convenience of comparison, one reverse motion period is defined as the period when the two-phase interface starts to move from the cavity back to the hydrofoil leading point. Figure 8 shows thereverse movement cycle of the two-phaseinterfacefordifferentcavitationnumbers.The reverse movement period obviously increases with the decrease of the cavitation number (shown as absolute values), and a sharp change occurs between the cavitation numbers 0.61 and 0.57. For various cavity lengths with different cavitation numbers, a relative backward motion cycle, T*, is defined as the original period multiplied by the cavity length at σ=0.74,and then divided by the cavity length at the current cavitation number, which is also given in Fig.8 (shown as the relative value). It can be shown that the relative adverse motion period sees a linear decrease with the decreasing cavitation number. This trend shows that when the cavitation number continuously drops, the volume exchange and the mass transfer between the water and the vapor phases inside the cavity become increasingly frequent, resulting in more even supercavitating flows.

    Fig.9 Transient pressure distributions

    Fig.10 Schematic diagram of the special section location

    2.2Pressure distribution

    Fig.11 Instantaneous pressure distributions at the level sections

    The dynamic pressure is obtained to understand the reverse motion mechanism of the two-phase interface in the cavity. Figure 9 shows the transient pressure distributions. A low pressure area appears behind the hydrofoil suction surface in all stages (shown in black color), and extends downstream with the reduction of the cavitation number. Unlike the mainstream area with high pressure (shown in gay color), thepressure fluctuation near the wake of the low pressure region and the vortex entrainment make the local vortex to have an instantaneous variation, which behaves as the periodic shedding of the upper and lower vortices (corresponding to the cavitation structures in Fig.6). For example, as shown in Fig.9 with σ=0.74, at =0mst, the upper vortex near the cavity tail occupies a dominant position, which rotates in a clockwise direction (shown in red arrow), at =4mst, the lower vortex is much stronger, which rotates in a counterclockwise direction (indicated with arrows), then a periodic vortex shedding forms. The similar behavior can be seen in the case of σ=0.67and σ=0.54. As a result of adverse pressure gradient effects, some of the upper and lower vortices near the cavity enclosing point begin to shed downstream, whereas others move reversely toward the cavity front, resulting in a dynamic distribution of the water-vapor mixture phase in the cavity.

    To further reveal the pressure distribution characteristics in the supercavitation region, Fig.10 gives the position of three relative level sections, wherein the hydrofoil leading vertex is chosen as the coordinate origin, and three ordinate points are fixed: the hydrofoil leading vertex with =0y, the middle point of the suction surface with y=0.069c, and the trailing edge vertex with y=0.222c. Figure 11 gives the transient pressure variation in the corresponding positions,where the reference vaporization pressure is 3 168 Pa. Apparently, in all cases: (1) Corresponding to the cavity tail, one position exists at which the pressure dramatically decreases to a level below the vaporization pressure and reaches a minimum value, then back around the vaporization pressure. The same pressure variation trend occurs at the sections of L1 and L2,whereas the opposite occurs at the L3 section. Obviously, it corresponds to the alternate shedding of the upper and lower vortices formed at the cavity tail,where the change of the vortex-core position results in the variation of the lowest pressure location. (2) Violent pressure fluctuation occurs at some place inside the cavity, the local pressure can be seen to suddenly increase and exceeds the vaporization pressure, then rapidly decreases to a certain level near the latter. The fluctuation surface moves reversely along with time. A comparison with the above vapor volume fraction distribution demonstrates that the pressure fluctuation position corresponds exactly the water-vapor twophase interface location, and the reverse movement of the pressure fluctuation position results in the adverse motion of the two-phase interface. The pressure suddenly increases because a pure vapor region is formed in front of the two-phase interface, where the liquid is completely vaporized. (3) Compared with the watervapor mixture area at the back of the cavity, the pressure in the vapor area is inversely and slowly reduced,having a minimum value near the suction side. Moreover, the pressure distribution is more uniform in the vapor area because the frequent volume and mass transfers in the water-vapor mixture area results in a continuous pressure fluctuation. (4) At the vertical direction in the cavity, the pressure gradually increases from the L1 section to the L3 section, i.e., the upper pressure is always greater than the lower pressure in the cavity.

    Fig.12 Pressure distributions on the vertical sections

    To better display the pressure variation perpendicular to the flow direction, the relative instantaneous pressure distribution at the vertical section of =1.4xc is obtained based on the coordinates mentioned above, as shown in Fig.12. Except at σ=0.74 (in the wake region), the pressure in the mainstream area gradually decreases with the decreasing cavitation number, while behaving relatively stable in the cavitation region. Moreover, a great pressure gradient exists in the cavity boundary. Figure 13 shows the transient pressure distributions on several vertical sections for various cavitation numbers. The cavity length is shorter at =σ 0.74, hence the latter two sections (=1.5xc and =x 2.0c) are located in the wake region, where the local pressure dramatically decreases. However, the same characteristics can be found in all stages: (1) the minimum pressure position is obtained within the shear layer, beginning from the foil leading point (=0)x. A strong shear layer forms from the hydrofoil leading vertex because of its special structure. The continuous vortex motion in the shear layer leads to a local pressure reduction. (2) At every vertical section, the pressure monotonously increases from the shear layer to the upper part of the cavity. (3) A large pressure gradient exists near the upper and lower boundaries of the cavity, which indicates that the cavitation boundary is located at the position where the pressure sharply increases. The cavitation region expands with the decreasing cavitation number, whereas the pressure gradient at the boundary is gradually reduced.

    Fig.13 Instantaneous pressure distributions on the vertical sections

    3. Conclusions

    The flow characteristics in different supercavitating regimes around a hydrofoil are simulated through a related numerical model. The reverse motion mecha-nism of the interface between the water and vapor phases is analyzed based on the cavitation structure and the pressure distribution. The following is a summary of the main findings:

    (1) A reverse motion of the water-vapor mixture phase occurs within the cavity in all supercavitating stages. The relative adverse motion period sees a linear decrease with the reduction of the cavitation number. This result shows that when the cavitation number drops continuously, the volume exchange and the mass transfer between the water and vapor phases inside the cavity become more frequent, resulting in a more even supercavitating flow.

    (2) The effect of the adverse pressure gradient induces the adverse pressure fluctuation inside the cavity. Corresponding to the two-phase interface position inside the cavity, a violent pressure fluctuation occurs. This fluctuation can be observed from the local pressure that suddenly increases and exceeds the vaporization pressure, then rapidly decreases to a level near the latter, the fluctuation surface moves reversely along with time. In short, the variation in the pressure fluctuation position results in a reverse motion of the two-phase interface in the cavity.

    (3) Compared with the water-vapor mixture area at the back of the cavity, the pressure in the vapor area becomes inversely and slowly reduced, with a minimum value near the suction side. Unlike the mainstream area, a higher pressure gradient occurs near the cavity boundary, and the minimum pressure position is related to the shear layer beginning from the foil leading point.

    In this paper, we mainly focus on the qualitative analyses of the flow structure and the pressure distribution in the cavity. The adopted parameters including the inlet velocity, the chord length and the lift angle of the foil are all constants. In a future research, the systematic simulation will be carried out under different conditions, and some experimental means will be considered to directly obtain the pressure distribution, the quantitative variation will further reveal the mechanism of the reverse flow in the supercavitation region.

    Acknowledgements

    The sponsoring of the China scholarship fund,the Beijing Key Laboratory of Passive Safety Technology for Nuclear Energy is gratefully appreciated.

    References

    [1]WANG G., CAO S. and IKOHAGI T. Study on behaviors of cavitating vortices in flows around a hollow-jet valve[C]. Third International Conference on Fluid Machinery and Fluid Engineering. Beijing, China, 2000, 451-455.

    [2]WANG G., SENOCAK I. and SHYY W. et al. Dynamics of attached turbulent cavitating flows[J]. Progress in Aerospace Sciences, 2010, 37(6): 551-581.

    [3]WEILAND C., VLACHOS P. P. Time-scale for critical growth of partial and supercavitation development over impulsively translating projectile[J]. International Journal of Multiphase Flow, 2012, 38(1): 73-86.

    [4]KIM Seon-Hong, KIM Nakwan. Hydrodynamics and modeling of a ventilated supercavitating body in transition phase[J]. Journal of Hydrodynamics, 2015, 27(5): 763-772.

    [5]PAN Zhan-cheng, LU Chuan-jing and CHEN Ying. Study on the characteristic of a forced pitching supercavitating vehicle[J]. Journal of Shanghai Jiaotong University(Science), 2015, 20(6): 713-720.

    [6]LIU Tao-tao, WANG Guo-yu and DUAN Lei. Assessment of a modified turbulence model based experiment results for ventilated supercavity[J]. Transactions of Beijing Institute of Technology, 2016, 36(3): 247-251(in Chinese).

    [7]KARIM M. M., AHMMED M. S. Numerical study of periodic cavitating flow around NACA0012 hydrofoil[J]. Ocean Engineering, 2012, 55(15): 81-87.

    [8]ZHI S. Numerical investigations of supercavitation around blunt bodies of submarine shape[J]. Applied Mathematical Modelling, 2013, 37(20-21): 8836-8845.

    [9]WOSNIK M., ARNDT R. E. A. Measurements in high void-fraction bubbly wakes created by ventilated supercavitation[J]. Journal of Fluids Engineering, 2013,135(1): 011304.

    [10] SHAFAGHAT R., HOSSEINALIPOUR S. M. and NOURI N. M. et al. Shape optimization of two-dimensional cavitators in supercavitating flows, using NSGA II algorithm[J]. Applied Ocean Research, 2008, 30(4): 305-310.

    [11] RASHIDIAI.,PASANDIDEH-FARDA M. and PASSANDIDEH-FARD M. et al. Numerical and experimental study of a ventilated supercavitating vehicle[J]. Journal of Fluids Engineering, 2014, 136(10): 101301.

    [12] LI X., WANG G. and ZHANG M. et al. Structures of supercavitating multiphase flows[J]. International Journal of Thermal Sciences, 2008, 47(10): 1263-1275.

    [13] TULIN M. P. Supercavitating propellers[C]. Proceedings of Fourth ONR Symposium on Naval Hydrodynamics. Rhode-Saint-Genese, Belgium, 2001, 239-286.

    [14] LI Xiang-bin, WANG Guo-yu and YU Zhi-yi et al. Multiphase fluid dynamics and transport processes of low capillary number cavitating flows[J]. Acta Mechanica Sinica, 2009, 25(2): 161-172.

    [15] OLIVER C. D., FORTES P. R. and REBOUD J. L. Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation[J]. Journal of Fluids Engineering, 2003, 125(1): 38-45.

    [16] SCHNERR G. H., SAUE J. Physical and numerical modelling of unsteady cavitation dynamics[C]. The 4th International Cconference on Mmultiphase Flow. New Orleans, USA, 2001.

    [17] SINGHAL A. K., ATHAVALE M. M. and LI H. Y. et al. Mathematical basis and validation of the full cavitation model[J]. Journal of Fluids Engineering, 2002, 124(3): 617-624.

    [18] LI Xiang-bin, WANG Guo-yu and ZHANG Bo et al. Evaluation of the RNG -kε model on the numerical simulations of supercavitating flows around a hydrofoil[J]. Chinese Journal of Hydrodynamics, 2008, 23(2): 181-188(in Chinese).

    (May 21, 2014, Revised August 18, 2014)

    * Project supported by the National Natural Science Foundation of China (Grant No. 50679001).

    Biography: Xiang-bin LI (1975-), Male, Ph. D., Lecturer

    猜你喜歡
    李楠王國(guó)
    在研究的路上鐫刻生命的印記
    Taking Robotics, AI, IoT to the World
    What Makes You Tired
    一滴水中的王國(guó)
    趣味(語文)(2020年5期)2020-11-16 01:34:54
    地下王國(guó)
    她的2000億打工王國(guó)
    逃離鼠王國(guó)
    ON A MULTI-DELAY LOTKA-VOLTERRA PREDATOR-PREY MODEL WITH FEEDBACK CONTROLS AND PREY DIFFUSION?
    建立新王國(guó)
    NBA特刊(2018年21期)2018-11-24 02:47:48
    Modeling of thermodynamics of ice and water in seasonal ice-covered reservoir *
    天堂影院成人在线观看| 国产视频一区二区在线看| 久久九九热精品免费| 免费看美女性在线毛片视频| 九色成人免费人妻av| 欧美黑人欧美精品刺激| 草草在线视频免费看| 亚洲美女视频黄频| 亚洲精品一区av在线观看| 午夜免费成人在线视频| 亚洲国产高清在线一区二区三| www日本在线高清视频| 欧美日本视频| 美女被艹到高潮喷水动态| 亚洲在线自拍视频| 一本精品99久久精品77| 欧美国产日韩亚洲一区| 欧美另类亚洲清纯唯美| 亚洲欧美精品综合久久99| 久久久久国产一级毛片高清牌| 成人亚洲精品av一区二区| 精品久久蜜臀av无| 国产单亲对白刺激| 日本黄色片子视频| 99久久精品国产亚洲精品| 女警被强在线播放| 欧美日韩福利视频一区二区| 久久午夜亚洲精品久久| 久久午夜亚洲精品久久| 香蕉国产在线看| 99在线视频只有这里精品首页| 久久午夜亚洲精品久久| 好男人在线观看高清免费视频| 高清毛片免费观看视频网站| 欧美成人性av电影在线观看| 黄频高清免费视频| 麻豆国产97在线/欧美| 午夜福利欧美成人| 最新中文字幕久久久久 | 日韩欧美国产一区二区入口| 午夜福利高清视频| 精品电影一区二区在线| 最好的美女福利视频网| 欧美在线黄色| 一区福利在线观看| 久久精品aⅴ一区二区三区四区| 老司机午夜福利在线观看视频| 又爽又黄无遮挡网站| 91麻豆av在线| 国产v大片淫在线免费观看| 欧美另类亚洲清纯唯美| 日本在线视频免费播放| 精品99又大又爽又粗少妇毛片 | 亚洲天堂国产精品一区在线| 黄色丝袜av网址大全| 欧洲精品卡2卡3卡4卡5卡区| 亚洲国产精品久久男人天堂| svipshipincom国产片| 久久久久精品国产欧美久久久| 淫妇啪啪啪对白视频| 国产1区2区3区精品| 757午夜福利合集在线观看| 国产精品亚洲av一区麻豆| 在线免费观看不下载黄p国产 | 国产又黄又爽又无遮挡在线| 一级毛片女人18水好多| 搡老熟女国产l中国老女人| 白带黄色成豆腐渣| 中文亚洲av片在线观看爽| 久久久久久人人人人人| 亚洲国产精品999在线| 久久婷婷人人爽人人干人人爱| 日本精品一区二区三区蜜桃| 亚洲av电影在线进入| cao死你这个sao货| 日本一二三区视频观看| 久久久久性生活片| 黄色丝袜av网址大全| 脱女人内裤的视频| 99国产精品99久久久久| 一二三四在线观看免费中文在| 成人国产一区最新在线观看| 日本与韩国留学比较| 久久亚洲精品不卡| 欧美黄色片欧美黄色片| 叶爱在线成人免费视频播放| 这个男人来自地球电影免费观看| 1024手机看黄色片| 久久草成人影院| 法律面前人人平等表现在哪些方面| x7x7x7水蜜桃| 一边摸一边抽搐一进一小说| 国模一区二区三区四区视频 | h日本视频在线播放| 亚洲色图av天堂| 日韩中文字幕欧美一区二区| 久久精品人妻少妇| 亚洲成人中文字幕在线播放| 午夜a级毛片| 国产高清视频在线观看网站| 国产av不卡久久| av在线天堂中文字幕| 最近最新中文字幕大全免费视频| tocl精华| 熟女少妇亚洲综合色aaa.| 亚洲在线观看片| 99久久99久久久精品蜜桃| 激情在线观看视频在线高清| 中文亚洲av片在线观看爽| 久久欧美精品欧美久久欧美| 日韩 欧美 亚洲 中文字幕| 精品乱码久久久久久99久播| 女生性感内裤真人,穿戴方法视频| 黄频高清免费视频| www.www免费av| 三级毛片av免费| 熟女少妇亚洲综合色aaa.| 高潮久久久久久久久久久不卡| 亚洲一区二区三区不卡视频| 婷婷精品国产亚洲av| 久久久久久人人人人人| 一a级毛片在线观看| 一二三四社区在线视频社区8| 91av网站免费观看| 高清毛片免费观看视频网站| 免费在线观看成人毛片| 最新在线观看一区二区三区| 成熟少妇高潮喷水视频| 中文字幕久久专区| 日本熟妇午夜| 成人av一区二区三区在线看| 2021天堂中文幕一二区在线观| 国产一区二区在线观看日韩 | 亚洲专区中文字幕在线| 国产一区二区在线观看日韩 | 日韩欧美国产一区二区入口| 97人妻精品一区二区三区麻豆| 国产欧美日韩一区二区三| 久久精品综合一区二区三区| 欧美在线一区亚洲| 在线观看日韩欧美| 国产激情久久老熟女| 久久久久久久精品吃奶| 99久久无色码亚洲精品果冻| 一进一出抽搐动态| 婷婷丁香在线五月| 日日夜夜操网爽| 白带黄色成豆腐渣| 伊人久久大香线蕉亚洲五| 啦啦啦免费观看视频1| www.999成人在线观看| 香蕉丝袜av| 国产淫片久久久久久久久 | 日本黄色视频三级网站网址| 亚洲中文字幕一区二区三区有码在线看 | xxx96com| 久久久久久久久中文| av福利片在线观看| 三级国产精品欧美在线观看 | 午夜激情欧美在线| 色尼玛亚洲综合影院| 一本综合久久免费| 999精品在线视频| 他把我摸到了高潮在线观看| 偷拍熟女少妇极品色| 日韩有码中文字幕| 天堂影院成人在线观看| netflix在线观看网站| 女同久久另类99精品国产91| 嫩草影视91久久| 久久婷婷人人爽人人干人人爱| 亚洲性夜色夜夜综合| 午夜免费激情av| 色av中文字幕| 成人一区二区视频在线观看| 国产乱人伦免费视频| 悠悠久久av| 国产精品久久久人人做人人爽| 老汉色av国产亚洲站长工具| 欧美xxxx黑人xx丫x性爽| 在线免费观看不下载黄p国产 | 成年免费大片在线观看| 亚洲av电影不卡..在线观看| 天堂网av新在线| 免费看光身美女| 婷婷精品国产亚洲av在线| 97超视频在线观看视频| 欧美日韩乱码在线| 在线观看日韩欧美| 亚洲国产中文字幕在线视频| 欧美绝顶高潮抽搐喷水| 欧美日韩黄片免| 国产91精品成人一区二区三区| 久久精品91无色码中文字幕| 亚洲男人的天堂狠狠| 国产伦在线观看视频一区| 成人特级黄色片久久久久久久| 午夜福利18| 亚洲精品美女久久久久99蜜臀| 亚洲九九香蕉| 国产美女午夜福利| 精品国内亚洲2022精品成人| 久久国产精品影院| 精品久久久久久,| 国产精品一区二区三区四区久久| 最近最新免费中文字幕在线| 亚洲精品456在线播放app | 欧美精品啪啪一区二区三区| 日韩人妻高清精品专区| 国产av不卡久久| 国产精华一区二区三区| 午夜成年电影在线免费观看| 久久精品国产综合久久久| 伊人久久大香线蕉亚洲五| 国产午夜精品久久久久久| 丝袜人妻中文字幕| 成人三级做爰电影| 亚洲va日本ⅴa欧美va伊人久久| 免费看a级黄色片| 亚洲精品国产精品久久久不卡| 嫁个100分男人电影在线观看| 亚洲欧美精品综合久久99| 在线观看舔阴道视频| 日本五十路高清| 叶爱在线成人免费视频播放| 久久久久性生活片| 99国产精品99久久久久| 18禁观看日本| 久久精品国产清高在天天线| 国产精品久久视频播放| 亚洲av电影不卡..在线观看| 亚洲最大成人中文| 无遮挡黄片免费观看| 18禁美女被吸乳视频| 国产97色在线日韩免费| 久久久久久久久中文| 午夜福利在线在线| 成人特级av手机在线观看| 天堂√8在线中文| 国产亚洲av嫩草精品影院| 精品国产乱子伦一区二区三区| 成年女人看的毛片在线观看| 一进一出好大好爽视频| 最近最新免费中文字幕在线| 国内久久婷婷六月综合欲色啪| 老熟妇乱子伦视频在线观看| 国产又色又爽无遮挡免费看| 91在线精品国自产拍蜜月 | 最新在线观看一区二区三区| 9191精品国产免费久久| 老熟妇仑乱视频hdxx| 黄色成人免费大全| 搡老岳熟女国产| 亚洲精品中文字幕一二三四区| 精品久久蜜臀av无| 欧美绝顶高潮抽搐喷水| 亚洲国产精品sss在线观看| 午夜福利免费观看在线| 成人av在线播放网站| 精品国内亚洲2022精品成人| 免费在线观看影片大全网站| www.999成人在线观看| 国产精华一区二区三区| 夜夜躁狠狠躁天天躁| 麻豆久久精品国产亚洲av| 久久精品影院6| tocl精华| 日日摸夜夜添夜夜添小说| 91字幕亚洲| 亚洲av电影在线进入| av片东京热男人的天堂| 久久精品国产综合久久久| 日韩欧美精品v在线| 亚洲欧美一区二区三区黑人| 特大巨黑吊av在线直播| 日韩精品中文字幕看吧| 亚洲欧美日韩无卡精品| 我要搜黄色片| 色av中文字幕| 一二三四在线观看免费中文在| 国产精品久久久久久精品电影| 男女视频在线观看网站免费| 亚洲avbb在线观看| 在线观看66精品国产| 哪里可以看免费的av片| 久久久国产成人免费| xxx96com| 两个人的视频大全免费| 日韩欧美三级三区| a在线观看视频网站| 欧美日韩亚洲国产一区二区在线观看| 免费无遮挡裸体视频| 波多野结衣高清作品| 国产成人aa在线观看| 99国产精品99久久久久| 国产三级中文精品| 99国产综合亚洲精品| 很黄的视频免费| 国产亚洲精品综合一区在线观看| 国产精品综合久久久久久久免费| 久久精品影院6| 操出白浆在线播放| 中文字幕久久专区| 国产高潮美女av| 国产午夜精品久久久久久| 久久这里只有精品19| 在线观看舔阴道视频| 美女大奶头视频| 免费搜索国产男女视频| 亚洲天堂国产精品一区在线| 精品一区二区三区视频在线观看免费| 床上黄色一级片| 麻豆成人午夜福利视频| 中文字幕高清在线视频| 精品人妻1区二区| cao死你这个sao货| 日本一本二区三区精品| 99久久久亚洲精品蜜臀av| 天堂√8在线中文| 三级毛片av免费| 国产精品国产高清国产av| 五月伊人婷婷丁香| 99久久无色码亚洲精品果冻| www.熟女人妻精品国产| 免费在线观看影片大全网站| 久久久久久国产a免费观看| 无人区码免费观看不卡| 国产熟女xx| 久久亚洲精品不卡| 欧美乱妇无乱码| 国产1区2区3区精品| 在线视频色国产色| a级毛片在线看网站| 免费在线观看日本一区| 两个人的视频大全免费| 午夜福利免费观看在线| 美女扒开内裤让男人捅视频| 精品国产亚洲在线| 亚洲狠狠婷婷综合久久图片| 亚洲欧美精品综合久久99| 激情在线观看视频在线高清| 久久这里只有精品中国| 久久九九热精品免费| 亚洲av片天天在线观看| 在线观看一区二区三区| 这个男人来自地球电影免费观看| 国产精品电影一区二区三区| 国产成人系列免费观看| 欧美xxxx黑人xx丫x性爽| 国产精品野战在线观看| 亚洲自拍偷在线| 国模一区二区三区四区视频 | 少妇丰满av| 午夜福利视频1000在线观看| 中文字幕高清在线视频| 最新中文字幕久久久久 | 午夜精品一区二区三区免费看| 村上凉子中文字幕在线| 久久久久九九精品影院| 久久久久免费精品人妻一区二区| 国产高清激情床上av| 麻豆成人午夜福利视频| 亚洲国产日韩欧美精品在线观看 | 国产精品久久久人人做人人爽| 欧美日韩一级在线毛片| 色老头精品视频在线观看| 99在线人妻在线中文字幕| 法律面前人人平等表现在哪些方面| 国产亚洲精品久久久久久毛片| 国产精品女同一区二区软件 | 国产精品99久久99久久久不卡| 床上黄色一级片| 女人被狂操c到高潮| 国产激情欧美一区二区| 日韩免费av在线播放| 精品国产乱码久久久久久男人| 99久久综合精品五月天人人| 国产熟女xx| 人人妻,人人澡人人爽秒播| 精品国产乱码久久久久久男人| 日韩有码中文字幕| 欧美乱码精品一区二区三区| 草草在线视频免费看| 毛片女人毛片| 国产伦在线观看视频一区| 免费看日本二区| 国产精品98久久久久久宅男小说| 国产精品一及| 无人区码免费观看不卡| 亚洲国产精品合色在线| 中文字幕人妻丝袜一区二区| e午夜精品久久久久久久| 国产亚洲av高清不卡| 听说在线观看完整版免费高清| 欧美中文综合在线视频| 免费在线观看影片大全网站| 精品国产乱码久久久久久男人| 欧美中文综合在线视频| 两个人的视频大全免费| 一夜夜www| 精品一区二区三区av网在线观看| 婷婷精品国产亚洲av在线| 亚洲国产精品成人综合色| 怎么达到女性高潮| 搡老岳熟女国产| 两性夫妻黄色片| 在线十欧美十亚洲十日本专区| 亚洲av美国av| 夜夜夜夜夜久久久久| 在线a可以看的网站| 亚洲av日韩精品久久久久久密| 欧美中文综合在线视频| 亚洲专区中文字幕在线| 51午夜福利影视在线观看| 九色国产91popny在线| 欧美日本亚洲视频在线播放| 性色avwww在线观看| 看黄色毛片网站| 九九久久精品国产亚洲av麻豆 | 亚洲av成人精品一区久久| 在线视频色国产色| 国产欧美日韩一区二区精品| 欧美日韩国产亚洲二区| 一进一出好大好爽视频| 国产成人一区二区三区免费视频网站| 曰老女人黄片| 色综合站精品国产| 欧美成人免费av一区二区三区| 精品国产亚洲在线| 国产伦人伦偷精品视频| 校园春色视频在线观看| 一本久久中文字幕| 国产精品 欧美亚洲| 国产爱豆传媒在线观看| 国产精品久久视频播放| 国产成人精品久久二区二区免费| 亚洲美女视频黄频| 国内精品久久久久久久电影| 亚洲人成网站在线播放欧美日韩| 日本a在线网址| 国产高清三级在线| 巨乳人妻的诱惑在线观看| 给我免费播放毛片高清在线观看| 久久精品亚洲精品国产色婷小说| 99在线人妻在线中文字幕| 国产精品永久免费网站| 嫩草影视91久久| 成人无遮挡网站| 国产成+人综合+亚洲专区| 97超视频在线观看视频| 欧美极品一区二区三区四区| 免费在线观看影片大全网站| 免费搜索国产男女视频| 中文资源天堂在线| 国产蜜桃级精品一区二区三区| 欧美黄色淫秽网站| 免费在线观看影片大全网站| 高清毛片免费观看视频网站| 搡老熟女国产l中国老女人| 久久久水蜜桃国产精品网| 丰满人妻熟妇乱又伦精品不卡| 日本与韩国留学比较| 成人av在线播放网站| 亚洲国产高清在线一区二区三| 国产成人精品无人区| 91在线观看av| 在线免费观看的www视频| 国产精品 欧美亚洲| 神马国产精品三级电影在线观看| 毛片女人毛片| 亚洲人成网站在线播放欧美日韩| 国产黄色小视频在线观看| 亚洲av成人一区二区三| 99热只有精品国产| 人妻久久中文字幕网| 日本五十路高清| 国产成人aa在线观看| 国产一级毛片七仙女欲春2| 国模一区二区三区四区视频 | 亚洲av日韩精品久久久久久密| 久久午夜综合久久蜜桃| 少妇人妻一区二区三区视频| 国产成人aa在线观看| 女同久久另类99精品国产91| 国产黄色小视频在线观看| 麻豆一二三区av精品| 亚洲精品美女久久久久99蜜臀| 亚洲精品乱码久久久v下载方式 | 免费在线观看影片大全网站| www.www免费av| 一个人看的www免费观看视频| 嫩草影院精品99| 99久久精品一区二区三区| 亚洲av免费在线观看| av黄色大香蕉| 国产主播在线观看一区二区| 久久久国产成人精品二区| 好男人电影高清在线观看| 午夜免费激情av| 婷婷精品国产亚洲av| а√天堂www在线а√下载| 香蕉久久夜色| 国产三级黄色录像| 九色国产91popny在线| 国产不卡一卡二| 国产美女午夜福利| 老司机深夜福利视频在线观看| 国产精品一区二区三区四区久久| 亚洲精品456在线播放app | 午夜福利在线观看免费完整高清在 | 国内揄拍国产精品人妻在线| 很黄的视频免费| 一进一出抽搐gif免费好疼| 国产亚洲精品久久久com| 夜夜看夜夜爽夜夜摸| 久久亚洲真实| 少妇熟女aⅴ在线视频| 亚洲狠狠婷婷综合久久图片| 啦啦啦韩国在线观看视频| 国产91精品成人一区二区三区| 国内久久婷婷六月综合欲色啪| 日韩 欧美 亚洲 中文字幕| av片东京热男人的天堂| 一级作爱视频免费观看| 日本成人三级电影网站| 免费人成视频x8x8入口观看| 精品酒店卫生间| 91久久精品国产一区二区成人| 亚洲在久久综合| 一个人免费在线观看电影| 中文字幕av在线有码专区| 丰满乱子伦码专区| 成人一区二区视频在线观看| 高清av免费在线| 国产精品电影一区二区三区| АⅤ资源中文在线天堂| 国产又黄又爽又无遮挡在线| 高清日韩中文字幕在线| 蜜桃久久精品国产亚洲av| 日韩中字成人| 一卡2卡三卡四卡精品乱码亚洲| 免费观看a级毛片全部| 免费看av在线观看网站| 亚洲国产精品sss在线观看| 国产精品一二三区在线看| 神马国产精品三级电影在线观看| 最近的中文字幕免费完整| 久久热精品热| 亚洲人与动物交配视频| 22中文网久久字幕| 久久久久网色| 国产女主播在线喷水免费视频网站 | 免费播放大片免费观看视频在线观看 | 晚上一个人看的免费电影| av国产免费在线观看| 亚洲久久久久久中文字幕| 久久亚洲精品不卡| 直男gayav资源| 色综合亚洲欧美另类图片| 国产乱人偷精品视频| 亚洲av.av天堂| 男女那种视频在线观看| 黑人高潮一二区| 日本一本二区三区精品| 亚洲18禁久久av| 国产欧美日韩精品一区二区| 亚洲欧洲日产国产| 亚洲国产日韩欧美精品在线观看| 国产精品精品国产色婷婷| 久久久久免费精品人妻一区二区| 我要看日韩黄色一级片| 亚洲无线观看免费| 国产高清有码在线观看视频| 少妇高潮的动态图| 日韩欧美精品免费久久| 亚洲国产高清在线一区二区三| 男人的好看免费观看在线视频| 中文亚洲av片在线观看爽| 日韩高清综合在线| 身体一侧抽搐| 日韩 亚洲 欧美在线| 日产精品乱码卡一卡2卡三| 干丝袜人妻中文字幕| 欧美日本亚洲视频在线播放| 亚洲欧美精品综合久久99| 精品国产三级普通话版| 日韩欧美国产在线观看| 五月玫瑰六月丁香| 日韩一区二区视频免费看| 噜噜噜噜噜久久久久久91| 伦精品一区二区三区| 亚洲欧美日韩卡通动漫| 三级经典国产精品| videossex国产| 国产一级毛片七仙女欲春2| 熟女人妻精品中文字幕| 只有这里有精品99| 人妻少妇偷人精品九色| 亚洲精品久久久久久婷婷小说 | 国产精品不卡视频一区二区| 校园人妻丝袜中文字幕| 久久久久久久久久久免费av| 国语对白做爰xxxⅹ性视频网站| 人妻夜夜爽99麻豆av| 日韩 亚洲 欧美在线| 国产三级在线视频| 内地一区二区视频在线| 亚州av有码| 在线播放国产精品三级| 亚洲国产欧洲综合997久久,| 久久久欧美国产精品| 69人妻影院| 插逼视频在线观看| 精品人妻熟女av久视频| 自拍偷自拍亚洲精品老妇|