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

    A selected review of vortex identification methods with applications *

    2018-10-27 09:09:58YuningZhang張宇寧XuQiu裘勖FeipengChen陳飛鵬KaihuaLiu劉凱華YuningZhang張宇寧XiangruiDong董祥瑞ChaoqunLiu
    關(guān)鍵詞:凱華祥瑞

    Yu-ningZhang(張宇寧),XuQiu(裘勖),Fei-pengChen(陳飛鵬),Kai-huaLiu,(劉凱華),Yu-ningZhang(張宇寧),Xiang-ruiDong(董祥瑞),ChaoqunLiu

    1. Key Laboratory of Condition Monitoring a ndControlfo rPowe rPlant Equipm ent(Min istryofEdu cation),School of Energy, Power andMechanicalEngineering,NorthChinaElectricPowerUniversity,Beijing 102206 , China

    2. BeijingSo uth-to-North Wat erDiversionTu anchengLak eMan agementO ffi ce,Beijing100195 ,China

    3. College ofMechanicalandTransportationEngineering,ChinaUniversityofPetroleum-Beijing,Beijing 102249, China

    4. Beijing Key Laboratory of ProcessFluidFiltrationandSeparation,ChinaUniversityofPetroleum-Beijing,Beijing 102249, China

    5. National Key Laboratory of Transient Physics, Nanjing University of Science and Technology, Nanjing 210094,China

    6. Department of Mathematics, University of Texas at Arlington, Arlington, Texas, USA

    Abstract: In the present review, recent progress on the vortex identification methods are introduced with a focus on the newly proposed omega method (Ω method).The advantages of Ω method are summarized with many illustrating examples.Furthermore, comparing with other existing methods (e.g., Q criterion and 2λ criterion), one of the characteristics of Ω method is its independence on the chosen threshold values for vortex identifications.The important parameters involved for the practical applications of Ω method are further discussed in detail together with the physical interpretation of the Ω and some suggestions of the future work.Other emerging topics (e.g., Lagrangian coherent structure and Rortex) are also introduced with comments.

    Key words: Vortex identification, turbulence generation, Ω method, vorticity concentration, Galilean invariant, Rortex

    Authors? biography

    Yu-ning Zhang (1983-), Associate Professor of School of Energy, Power and Mechanical Engineering,North China Electric Power University, Beijing,China.He obtained Bachelor of Engineering in Power Engineering and Engineering Thermodynamics (2006)and Master of Engineering in Fluid Machinery and Engineering (2008) both in Tsinghua University,Beijing, China and Ph.D.in Fluid Dynamics (2013),The University of Warwick, UK.His primary research interests include vortex identification, cavitation,high-speed photography and hydraulic machineries.He has published over 40 papers on high-profile international journals as leading/corresponding author.He has given over ten invited talks on leading conferences.He also serves as the Chairman of International Cavitation Forum 2016, an Associate Editor of IET Renewable Power Generation, member of scientific/ organizing committee of many well-established international conferences.

    Chaoqun Liu (1945-), Distinguished Professor of Department of Mathematics, University of Texas at Arlington, Arlington, Texas, USA.He also servers as the Director of Center for Numerical Simulation and Modeling.He obtained Bachelor of Science in Fluid Mechanics (1967) and Master of Science in Computational Fluid Dynamics (1981) both in Tsinghua University, Beijing, China and Ph.D.in Applied Mathematics (1989), University of Colorado at Denver, USA.His fields of interest include vortex identification, numerical analysis, multigrid, multilevel adaptive methods, computational fluid dynamics,direct numerical simulation, large eddy simulation,high-order numerical scheme, numerical combustion and software development.He has published over 200 journal/conference papers together with eight books.He served as the Principal Lecturer of many short courses relating with the turbulence generation and sustenance in boundary layer held both in USA and China.He was also the Chairman of the first (1997)and the third (2001) AFOSR International Conference on Direct Numerical Simulation and Large Eddy Simulation.

    Introduction

    Vortex is a quite common phenomenon in the natural world (e.g., vortex induced by the flying hawkmoth[1]) and our daily life.Furthermore, vortex is widely involved in a wide range of flows in the practical engineering applications including ultrasonic biomedical treatment[2-3], airplane wing aerodynamics[4], shock-wave induced flow[5], natural convection[6], wind turbine dynamics[7-9], jetting phenomenon[10].Specifically, vortex is also strongly related with various kinds of hydrodynamics including cavitation[11], pump performance[12-13], and tip clearance vortex[14-15].For example, vortex related cavitation phenomenon include bubble-particle interactions in the silt-laden flow[16-20], wave propagation in bubbly flow[21-22], collapsing dynamics[23-26], acoustic bubble dynamics (e.g., chaotic oscillations[27], nonlinear oscillator[28-29], damping[30]and stability[31-32]), and cavitation inception[33-34].Furthermore, many industrial-level flows e.g., reversible pump turbine for largescale energy storage, also suffer from many vortex induced phenomenon including rotating stall[35],vibrations[36], pressure fluctuations[37-41], tip clearance vortex[42-43], hysteresis[44-45], hump stability[46-47].Hence, methods for vortex identification are a critical issue in fluid mechanics especially hydrodynamics for revealing the underlying vortical structures with further improvement of the theory.

    In the literature, this topic has been intensively addressed by many researchers[48-68]with various kinds of vortex identification methods proposed (e.g.,Qcriterion[56]and2λcriterion[58]).For a recent review, readers are referred our previous paper with a clear classification of the existing methods with applications examples[48].Recently, Liu and his collaborators proposed several new vortex identification methods[49-55]with plenty of applications in the literature[69-70].In this review, recent progress on the vortex identification methods are selected with a focus on the newly proposed omega method.A brief introduction of each section of the present review will be given as follows.In Section 1, several typical vortex identification methods are briefly introduced with some comments.In Section 2, the omega method(Ωmethod) is introduced with a summary of its advantages based on comparisons with other existing methods (e.g.,Qcriterion[56]and2λcriterion[58]).In Section 3, several demonstrating applications of theΩmethod are given and interpreted.In Section 4,several possible future works are given and discussed in detail with a physical interpretation of theΩand other emerging topics (e.g., Lagrangian coherent structure and Rortex).In Section 5, the concluding remarks of the present review are given together with the available code for theΩvortex identification.

    1.An overview of typical vortex identification methods

    In this section, several well-employed vortex identification methods will be briefly introduced.For the convenience of the further discussions, the characteristic equation of velocity gradient tensor is given as

    HereP,Q,Rare three invariants of the velocity gradient tensor.Furthermore, the velocity gradient(?V) could be decomposed into two parts as follows

    whereAis the symmetric part,Bis the antisymmetric part.

    1.1 Q criterion

    Qcriterion is directly derived based on the second invariantQof the velocity gradient tensor given in Eq.(1) with the following expression[56]

    1.2 λ2 criterion

    The2λcriterion defines the vortex of the fluid flow as the local pressure minimum of the certain fluid domain[58].If one defines three real eigenvalues ofA2+B2asλ1,λ2andλ3with the following order (λ1≥λ2≥λ3), theλ2criterion corresponds to the existence of two negative eigenvalues (i.e.,λ2<0)[58].

    1.3 Δ criterion

    Based on the critical point theory, Chong et al.[57]defined the vortex region with

    1.4 Comparisons between different existing methods

    In this section, a brief note on the comparisons between different existing methods will be given.The aforementioned methods all require to adjust the threshold values in order to deliver a proper result.Hence, the common limitation of the aforementioned three methods is their dependence on the above thresholds.Specifically, different thresholds may lead to different vortex structures for a given flow case.However, a rigorous criterion for the determination of this threshold is still absent.

    Figure 1 shows the vortex identification results by the2λcriterion (a) andQcriterion (b) for the Taylor-G?rtler vortical structures.As shown in Fig.1,there are no essential difference among the two for this case if a proper threshold is chosen.Similar findings were also reported by the Dubief and Delcayre[71](referring to Fig.5 of the Ref.[71]).

    Fig.1 Comparison of the vortex identification by the 2λ criterion (a) and Q criterion ((b), identified with-2Q).Velocity vectors are also shown for the comparions.This figure was adapted from Figs.5(a), 5(c) of Fraňa et al.[72]

    2.Ω method

    2.1 Description of Ω method

    The definition ofΩcan be described as follows[49]:

    whereaandbare the square of Frobenius norm ofAandB(referring to Eq.(1)).For the practical calculation,Ωis further modified as follows[49]

    Fig.2 An illustrating example of the influences of the Q values on the vortex identification.Here, physical identification of the T-S wave is selected for the discussions.This figure was adapted from the Fig.3 of Liu et al.[49]

    where ε is a small positive number employed to avoid the division by zero.For detailed determination of parameter ε, readers are referred to Section 2.2.

    Different with other existing vortex identification methods, the Ω method has a normalized value ranging from 0 to 1 since both a and b are not negative.Meanwhile, new Ω method has a clear physical meaning through the definition of the vortices corresponding to the vorticity overtaking the deformation.

    Fig.3 Comparison for the vortex identification of the T-S wave.The chosen values of the Q criterion (Q =10-3) and 2λ criterion 32 (λ = -10-3 ) are optimized based on the author?s experience.This figure was adapted from the Fig.5 of Liu et al.[49]

    2.2 Advantages of Ω method

    The advantages of Ω method could be summarized as follows.Here, vortex identifications of the T-S wave and the micro wave generator are selected as examples for discussions.

    Figure 2 shows an illustrating example of the influences of the Q values (Q=10-3, 10-4and 10-2respectively) on the vortex identification.One can find that the choice of the Q values will strongly affect the identified vortex.However, it is quite difficult to determine an appropriate threshold value for the vortex identification based on the Q criterion.Similar dilemma also happens for other identification methods (e.g.,2λ criterion).

    Fig.4 Comparison for the vortex identification of micro wave generator.The chosen values of the Q criterion (Q =10-3 ) and λ2 criterion (λ = -10-3) are optimized based on the author?s experience.This figure was adapted from the Fig.9 of Liu et al.[49]

    For the Ω method, the identified vortex is universal.Figures 3, 4 show the comparisons between the λ2/Q criterions and the Ω method for the vortex identification of the T-S wave and micro wave generator respectively.The chosen values of the Q criterionand λ2criterionare optimized based on the author?s experience.It could be found that without the subjective adjustment of the values (e.g., the process involved for the choice of a proper threshold during the usage of theλ2|Qcriterions), theΩmethod could successfully capture the primary characteristics of the vortex structure.As shown in Figs 5, 6, the results of vortex identification is quite insensitive to the chosen values ofΩ.

    Fig.5 The influences of the values of Ω on the Λ vortex identification.This figure was adapted from the Fig.6 of Liu et al.[49]

    Fig.6 The influences of the values of Ω on the vortex identification during the boundary layer transition.This figure was adapted from the Fig.7 of Liu et al.[49]

    Furthermore, during the practical usage of theQcriterion, for the large values ofQ, only strong vortex will be captured with the weak vortex skipped while for the small values ofQ, there are too many weak vortex, leading to the smear of the vortex structure.However, for theΩmethod, both the strong and weak vortex could be captured simultaneously (as shown in Figs.3, 4).

    For the completeness, the influences of the parameterεin Eq.(8) on theΩmethod will be further discussed.In our previous paper, a proper determination method of parameterεwas proposed such as[50]

    Fig.7 Influences of values of ε on the vortex identification of the boundary layer transition.This figure was adapted from the Fig.1 of Dong et al.[50]

    Fig.8 Vortical structure visualization using the proposed determination method in different time steps during the boundary layer transition.This figure was adapted from the Fig.3 of Dong et al.[50]

    Figure 7 shows the tests of theεdetermination method for the vortex identification att=8.2Tduring the boundary layer transition on a flat plate based on a high order direct numerical simulation(DNS).Figure 7(b) shows the vortex identification by employing Eq.(9) for the determination ofε.For the comparisons, smaller and larger values ofεare also chosen to visualize the vortex structures as shown in Figs.7(a), 7(c).In Fig.7(a), too many clouds is shown without a clear demonstration of the vortex structure while in Fig.7(c), some significant vortex structures will not be identified.Hence, the Eq.(9) is a proper choice for the determination ofε.Figure 8 further shows the identified vortical structure of the boundary layer transition by theΩmethod in different time steps.As shown in Fig.8, many typical vortical structures and its time-domain evolutions could be well identified.

    3.Applications of Ω method

    3.1 Vortex shedding in the wake flow of moving bodies

    The vortex shedding dynamics induced by the activities of the human bodies plays an important role on the determination of the contaminant dispersions(e.g., bacteria) in indoor environments especially hospital inpatient wards, airline cabin and certain factories.Hence, the identification of the vortex patterns generated by the movement of the human bodies is essential to fulfil the above targets.

    Recently, Tao et al.[69]applied theΩmethod for the identification of vortices in the wake flow of two typical types of moving manikins (representing the human bodies).For the completeness, two scaled manikins with different shapes (slim and larger ones respectively) were carefully manufactured using 3-D printing technology.Meanwhile, an experimental study using smoke visualization system was also conducted with the results compared with the omega identification ones.

    Figure 9 shows the comparisons of the vortex structures of the wake flows induced by the two kinds of the moving manikins in different time steps.Figure 9(a) shows the slim one and Fig.9(b) shows the larger one.Both the smoke visualization (upper subplots)and the vortex identification through theΩmethod(lower subplots) are all shown.Colors in the figure represent the velocity magnitude.Generally speaking,the vortex identification through theΩmethod agrees well with the experimental observation especially in terms of the separation angle at the manikin heads and the primary patterns of the vortex structures.Furthermore,Ωmethod could also reveal the sheared wake and the entrained flow behind the moving manikins.

    3.2 Drifting snow in the turbulent atmosphere boundary layer

    In the turbulent atmosphere boundary layer with high latitudes, drifting snow frequently occurs as an important type of aeolian multiphase flow.The study of the drifting snow is of great importance from the viewpoints of the glaciological and the hydrological aspects.For example, the drifting snow could lead to a re-distribution of the snow layer, which is critical for the evaluation of possible natural geological hazards.Furthermore, through the dispersion and the transportation of the snow particles, the mass of the ice sheets will be affected, leading to the temporal or spatial imbalance.

    Recently, Huang and Wang[70]established a 3-D drifting snow computational model based on the large eddy simulation with a large amount of snow particle trajectories tracked.For the vortex identification,Ωmethod is employed to further reveal the vortex structures in the fully developed turbulent boundary layer.Figure 10 shows the primary vortex structures in the turbulent boundary layer with snow particles.As shown in Fig.10(b), regions with high concentrations of the snow streamers are of great values of vortical vorticity.Hence, the formation of the snow streamer is closely related with the local fluid vortex flow.This is quite different with the sand streamer, for which the strong surface wind is of great importance for the sand distributions.

    Fig.10 A demonstration of the vortex identification with Ω method in fully developed turbulent boundary layer with snow particles.(a) The vortex structures in the turbulent boundary layer.(b) The distributions of vortical vorticities larger than 2 000 s-2 (green patches) together with the snow particles (black points).This figure was adapted from the Fig.8 of Huang and Wang[70]

    3.3 Vortices in the reversible pump turbine

    Reversible pump turbine are widely employed in the pump hydro energy storage power plant.However,currently, the flow-induced severe instability are greatly affecting the daily normal operations of the units, leading to many unplanned stops of the power plants.Because of the frequent shifts among various kinds of working modes, there are a large amount of vortices shown in the flow passage of the reversible pump turbine.Hence, the vortex identification in the reversible pump turbines is important to reveal the inside vortex flow structures with further improvement of the design.

    Fig.11 The vortices identified by Ω method in the reversible pump turbine at 21° guide vane opening.This figure was adapted from Fig.3 of Zhang et al.[38]

    Generally speaking, five operational modes(including the turbine mode, the runaway mode, the turbine brake mode, the zero-flow-rate mode and the reverse pump mode) are usually involved during the daily operations of the turbine.Hence those five modes will be selected for the analysis of vortices using theΩmethod as shown in Zhang et al.[38].Figure 11 shows the primary characteristics and the distributions of vortices inside the turbine during five different working modes.As shown in Fig.11, for the turbine mode, limited vortices are shown in the impeller while for the reverse pump mode, a large amount of vortices occupies the majority of the turbine.Figure 12 further shows a quantitative comparison among five different working modes in terms of the vortices area.Among them, the zero-flow-rate mode is of the largest vortices area and the reverse pump mode is the second one.For other working modes, the difference is marginal.According to Fig.11,Ωvortex identification method could capture both weak (e.g., Fig.11(d)) and strong vortices (e.g.,Fig.11(e)) in the whole passage of the reversible pump turbines.

    Fig.12 A quantitative statistics of vortices areas identified by the Ω method for five different operational conditions of the reversible pump turbine at 21° guide vane opening.This figure was adapted from Fig.4 of Zhang et al.[38]

    4.Future works

    4.1 Determination of the parameter in Ω method and the interpretation of its physical meanings

    In theΩmethod, the determination of the parameterεis quite essential because it will affect the results of the vortex identification.As shown in Eq.(8), the initial intention for the addition of the small positive parameterεis to avoid the division by zero.Such treatment only involves the numerical details.Here, some emerging ideas will be given and discussed to further explore the physical meanings ofΩmethod together with the parameterε.

    For the convenience, theΩvalue for the case without the involvement of the parameter is defined as 0Ωwith the following expressions

    The addition of the parameterεcould be fulfilled through revisions of two parts (theaandbrespectively).For the convenience, the addition of theεwill be further discussed in terms of non-dimensional forms.The first method involves the revisions of terma(deformation term) such as

    with

    Here Δais the non-dimensional parameter relating with deformation term.

    If one draw an analogy between theΩdefinition and the concentration of a mass (e.g., the concentration of sugar in the water), many similarity could be revealed.For the readers? convenience, a new term “vorticity concentration (VC)” is newly proposed.The pure deformation can be defined as flows with VC=0 while the rigidly rotational flow is defined as flows with VC=1.According to theΩcriteria, the vorticity overtakes the deformation when the “vortex” is formed[49].It means that the vortex is formed in flows when VC>0.50.Hence,with the aid of the concept of VC, the physical meanings of the vortices could be clearly shown.Alternatively speaking, the addition of termaεin the Eq.(11) refers to the decrease of the vorticity concentration, leading to the decrease of theΩvalues.Quantitatively speaking, the above formulas could be revised as:

    or

    Here Δcis the non-dimensional parameter reflecting the variations of theΩdue to the addition of the parameterεain Eq.(11).

    Fig.13 The variations of the Δc and Ω with the parameter Δa

    Figure 13 shows the variations of the ΔcandΩwith the parameter Δa.In this example,0Ωis set to be 0.52 for the convenience of discussions.With the increase of Δa, Δcincreases dramatically while theΩdecreases from the 0.52 to 0.32.Especially,for Δalarger than 10-2, the decrease rate ofΩis quite significant (with a quick dilution of the vorticity).

    Alternatively, the other method for the addition ofεinvolves the revisions of the vorticity termbas follows:

    On the contrary of the first method, the revisions of termb(bεin the Eq.(15)) corresponds to the increment of the vorticity concentration.Different with the previous case, the revisions of termbinvolves the change of both numerator and denominator.Figure 14 shows the variations of the ΔcandΩwith the parameterΔb.With the increase of Δb,both ΔcandΩincrease (corresponding to the vorticity concentration enhancement).In the future,more examples relating with the parameter determination ofεwill be given.

    Fig.14 The variations of the Δc and Ω with the parameter Δb

    Hence, as a summary, for different methods of addition ofε(aεin Eq.(11) orbεin Eq.(15)), it could avoid the possible numerical errors but with some other effects, which should be noticed during the practical applications.With the revisions of the deformation term (εa), the weak vortex will be skipped and could not be clearly identified.On the contrary, with the revisions of the vorticity term (εb),some structures will be identified as the “fake” vortex,leading to the incorrectness.

    4.2 Lagrangian coherent structures

    Generally speaking, the trajectories of most typical fluid particles involved in the natural or industrial flows is highly dependent on the initial conditions.This feature causes great difficulties for the development of computational fluid models (e.g.,in terms of fluid motion predictions and comparisons)and also the experimental observations.However,different with the uncertainties of the single trajectory,Lagrangian coherent structures (LCSs) usually exists with a robust pattern and frames[65-68], which could be employed to represent the typical fluid structures.

    Recently, LCSs has made significant progress with experimental confirmations, coherent eddy transport ocean models and short-term forecasting potentials.For a comprehensive review, readers are referred to Haller[68].In the future, for industrial-level vortex flows (e.g., the aforementioned turbines), a 3-D LCSs is highly desirable and the computational efficiency is of great concern for those cases.

    4.3 Rortex

    In order to identify the vortical structures more precisely, Liu et al.[51,54-55]proposed a new vortex vector (named as “Rortex”) with a rigorous mathematical proven of its existence and an accurate definition.Comparing with the vortex, a clear relationship between the Rortex and the rigid fluid rotation is established.For the verifications of the effectiveness of the Rortex, several demonstrating examples were given including both the 2-D flow (e.g., Couette flow,rigid rotational flow) and the complex 3-D flow (e.g.,boundary layer transition on a flat plate).Moreover,many Rortex based vortex structures could be also utilized for the vortex identification e.g., Rortex vector field, Rortex lines and Rortex tubes.

    4.4 Differences between vortex identification and visualization

    There are great differences between the vortex identification and the vortex visualization.For the vortex identification, the main purpose is to identify the correct vortex structures for the various kinds of flows.However, for the vortex visualization, the primary target is to visualize the already identified vortex with a well-established view.Hence, the difference between theΩmethod and the other methods(e.g., theQcriterion and2λcriterion) is quite significant for the vortex identification.OnlyΩmethod could identify the correct vortex structures because it can identify both the strong and weak vortex.For other methods, the identified vortex is strongly dependent on the given values of a certain criterion.For large values, only strong vortex is identified with many important vortex structure omitted.For small values, many weak vortex is identified with too much noises, leading to the useless of the identified results.It looks like an unachievable task to strike a balance between the two trends.The omega method provides a solution to this longstanding dilemma for the vortex identification.In the future, the identification of the weak vortex is an emerging topic (e.g., with 0.50<Ω< 0.52) in order to reveal some unobserved vortex structures (e.g.,various kinds of vortex during the turbulence structure generation and the time/space domain evolution).

    5.Conclusions

    In the present review, recent progress on the vortex identification methods has been critically selected with a focus on the newly proposedΩmethod.The advan- tages ofΩmethod are demonstrated through comparisons with other existing methods (e.g.,Qcriterion and2λcriterion).The important parameters involved for the practical applications ofΩmethod are further discussed in detail with many suggestions of the future work.The following concluding remarks could be given:

    (1)Ωmethod overcomes the usage of the threshold values, leading to its independence of the user?s subjective experience.

    (2) The involved parameterεin theΩmethod is important for vortex identification and is quite worthy for the further investigation.

    (3) More application examples ofΩmethod is quite necessary for extending its applications into broad ranges of flows.

    TheΩcode (written both in Tecplot and Fortran as a supplementary material) is available through the following website:

    https://link.springer.com/article/10.1007%2Fs11433-016-0022-6

    Acknowledgements

    This work was supported by the Science Foundation of China University of Petroleum, Beijing(Grant No.2462016YJRC003).We are very glad to take this opportunity to appreciate the fruitful discussions with Prof.Lian-di Zhou.

    猜你喜歡
    凱華祥瑞
    《祥瑞》
    瑞意搖滾
    睿士(2024年5期)2024-05-24 17:37:36
    不忘初心 牢記使命 推進樂凱華光跨越發(fā)展
    樂凱華光印刷科技有限公司 公司簡介
    金豬慶祥瑞
    青年歌聲(2019年2期)2019-12-09 12:59:29
    Analysis of the vortices in the inner flow of reversible pump turbine with the new omega vortex identification method *
    祥瑞
    寶藏(2017年10期)2018-01-03 01:53:10
    浣溪沙·新鄰居
    揚子江詩刊(2017年3期)2017-11-13 20:07:52
    Numerical simulation of hydraulic force on the impeller of reversible pump turbines in generating mode*
    麒麟祥瑞(七絕)
    寶藏(2017年4期)2017-05-17 03:33:49
    高清毛片免费观看视频网站| 久久久久精品国产欧美久久久| 亚洲欧美一区二区三区黑人| 性色av乱码一区二区三区2| 99精品欧美一区二区三区四区| 欧美日本中文国产一区发布| www.www免费av| 大码成人一级视频| 国产亚洲精品综合一区在线观看 | 十分钟在线观看高清视频www| 久久久久久久久免费视频了| 九色亚洲精品在线播放| 男女做爰动态图高潮gif福利片 | 久久热在线av| 老司机在亚洲福利影院| 欧美在线黄色| 欧美在线一区亚洲| 美女扒开内裤让男人捅视频| 亚洲av五月六月丁香网| 91成人精品电影| 亚洲精品粉嫩美女一区| 久久久国产成人免费| 丰满人妻熟妇乱又伦精品不卡| 女人爽到高潮嗷嗷叫在线视频| 免费在线观看影片大全网站| 国产一区二区三区综合在线观看| 国产av又大| 美女午夜性视频免费| 亚洲,欧美精品.| 两个人视频免费观看高清| 精品无人区乱码1区二区| 夜夜躁狠狠躁天天躁| 精品福利观看| 久久婷婷人人爽人人干人人爱 | 一边摸一边做爽爽视频免费| 亚洲视频免费观看视频| 一本久久中文字幕| 久久久久久久久免费视频了| 黄色 视频免费看| 亚洲第一青青草原| 91国产中文字幕| 丰满人妻熟妇乱又伦精品不卡| 色婷婷久久久亚洲欧美| 国产三级在线视频| 男女之事视频高清在线观看| 欧美一级毛片孕妇| 亚洲色图av天堂| 国产又色又爽无遮挡免费看| 中亚洲国语对白在线视频| 久久久久国产一级毛片高清牌| 韩国av一区二区三区四区| 97碰自拍视频| 成在线人永久免费视频| 女人精品久久久久毛片| 757午夜福利合集在线观看| 色在线成人网| 日日干狠狠操夜夜爽| 国产亚洲欧美98| 99国产精品99久久久久| 亚洲精品在线观看二区| 香蕉丝袜av| 精品无人区乱码1区二区| 波多野结衣巨乳人妻| 免费在线观看影片大全网站| 国产精品久久久人人做人人爽| 国产精品精品国产色婷婷| 女警被强在线播放| 欧美乱码精品一区二区三区| 国产真人三级小视频在线观看| 国产亚洲精品久久久久5区| av片东京热男人的天堂| 两性午夜刺激爽爽歪歪视频在线观看 | 久久精品成人免费网站| 女性被躁到高潮视频| 制服诱惑二区| 午夜影院日韩av| 天堂影院成人在线观看| 国产av精品麻豆| 国产亚洲精品一区二区www| 久久久久久久精品吃奶| 18禁国产床啪视频网站| 老熟妇仑乱视频hdxx| 国产伦一二天堂av在线观看| 成人永久免费在线观看视频| 免费在线观看视频国产中文字幕亚洲| 久久久久国产精品人妻aⅴ院| 曰老女人黄片| 中文亚洲av片在线观看爽| 国产亚洲av高清不卡| 欧美乱色亚洲激情| 搞女人的毛片| av视频免费观看在线观看| 久久精品91无色码中文字幕| 我的亚洲天堂| 国产精品久久久久久亚洲av鲁大| 黄频高清免费视频| 亚洲精品国产色婷婷电影| 麻豆成人av在线观看| 精品国产美女av久久久久小说| 人成视频在线观看免费观看| 最近最新中文字幕大全电影3 | 不卡一级毛片| 国产三级黄色录像| e午夜精品久久久久久久| 美国免费a级毛片| 最好的美女福利视频网| 久久精品91无色码中文字幕| 精品国内亚洲2022精品成人| 嫩草影视91久久| av在线天堂中文字幕| 国产单亲对白刺激| 老司机在亚洲福利影院| 在线天堂中文资源库| 午夜老司机福利片| 国产亚洲精品综合一区在线观看 | 午夜两性在线视频| 国产一区二区激情短视频| 十八禁人妻一区二区| 十分钟在线观看高清视频www| 老司机午夜十八禁免费视频| 又紧又爽又黄一区二区| 欧美乱码精品一区二区三区| 亚洲国产日韩欧美精品在线观看 | 久久久久久久精品吃奶| 久9热在线精品视频| 一卡2卡三卡四卡精品乱码亚洲| 国产精品一区二区在线不卡| 久久 成人 亚洲| 99久久99久久久精品蜜桃| 91九色精品人成在线观看| 757午夜福利合集在线观看| 国产高清视频在线播放一区| 美女高潮喷水抽搐中文字幕| 欧美激情高清一区二区三区| 国产成人影院久久av| 欧美在线一区亚洲| 久久午夜亚洲精品久久| 夜夜夜夜夜久久久久| 国产单亲对白刺激| 久久国产亚洲av麻豆专区| 亚洲五月色婷婷综合| 男女午夜视频在线观看| 成人国产一区最新在线观看| 亚洲欧美日韩无卡精品| 亚洲,欧美精品.| 久久中文字幕一级| 久久午夜亚洲精品久久| 日韩欧美在线二视频| 两性午夜刺激爽爽歪歪视频在线观看 | 久久精品影院6| 老鸭窝网址在线观看| 人人妻人人澡人人看| 午夜免费鲁丝| 一边摸一边抽搐一进一出视频| 每晚都被弄得嗷嗷叫到高潮| 两个人免费观看高清视频| 天堂动漫精品| 免费搜索国产男女视频| 99久久国产精品久久久| www国产在线视频色| 两个人看的免费小视频| 一级毛片高清免费大全| 国产99白浆流出| 黄色女人牲交| 欧美黄色淫秽网站| 日本vs欧美在线观看视频| 窝窝影院91人妻| 久久久久久免费高清国产稀缺| 日韩 欧美 亚洲 中文字幕| 国产精品二区激情视频| 午夜福利,免费看| 国产精品久久久久久人妻精品电影| 亚洲五月婷婷丁香| 免费无遮挡裸体视频| 亚洲欧美精品综合久久99| 日本免费a在线| 亚洲激情在线av| 日本五十路高清| 久热爱精品视频在线9| 国产在线精品亚洲第一网站| 亚洲七黄色美女视频| 国产一区在线观看成人免费| 岛国在线观看网站| 99久久99久久久精品蜜桃| 亚洲人成伊人成综合网2020| 欧美成人免费av一区二区三区| 欧美最黄视频在线播放免费| 欧美日韩精品网址| 亚洲中文字幕一区二区三区有码在线看 | 中文字幕高清在线视频| 在线国产一区二区在线| 淫妇啪啪啪对白视频| 淫妇啪啪啪对白视频| 又黄又爽又免费观看的视频| 在线观看www视频免费| 久久香蕉精品热| 欧美不卡视频在线免费观看 | 此物有八面人人有两片| 免费在线观看视频国产中文字幕亚洲| 亚洲国产精品久久男人天堂| 波多野结衣一区麻豆| 国产亚洲精品久久久久5区| 长腿黑丝高跟| 嫩草影院精品99| 欧美在线一区亚洲| 最新在线观看一区二区三区| 亚洲国产精品久久男人天堂| 欧美日本视频| 波多野结衣av一区二区av| 精品国产超薄肉色丝袜足j| 天天躁夜夜躁狠狠躁躁| 久久久久久久午夜电影| 午夜精品久久久久久毛片777| 久久人妻av系列| 日韩中文字幕欧美一区二区| 亚洲精品av麻豆狂野| 九色亚洲精品在线播放| 午夜免费成人在线视频| 午夜福利,免费看| 黄频高清免费视频| 精品国产一区二区久久| 天堂影院成人在线观看| 9热在线视频观看99| 亚洲少妇的诱惑av| 亚洲中文av在线| 在线观看日韩欧美| 国产日韩一区二区三区精品不卡| 国产精品99久久99久久久不卡| 一级a爱片免费观看的视频| 欧美日韩亚洲综合一区二区三区_| 精品国产一区二区久久| 看片在线看免费视频| 女人被躁到高潮嗷嗷叫费观| 亚洲久久久国产精品| tocl精华| 免费无遮挡裸体视频| 国产黄a三级三级三级人| 99久久综合精品五月天人人| 日本五十路高清| 久久性视频一级片| 国产精品一区二区精品视频观看| 中文字幕另类日韩欧美亚洲嫩草| 国产1区2区3区精品| av视频免费观看在线观看| aaaaa片日本免费| 久久国产精品影院| 性色av乱码一区二区三区2| 国产成人av激情在线播放| 在线视频色国产色| 热99re8久久精品国产| 国产精品98久久久久久宅男小说| 极品教师在线免费播放| 91老司机精品| 搡老岳熟女国产| 99精品久久久久人妻精品| 757午夜福利合集在线观看| 日韩三级视频一区二区三区| 久久午夜综合久久蜜桃| 亚洲第一av免费看| 性少妇av在线| 午夜精品国产一区二区电影| 青草久久国产| 一级黄色大片毛片| 满18在线观看网站| 国产伦一二天堂av在线观看| videosex国产| 黄色a级毛片大全视频| av欧美777| 美女 人体艺术 gogo| 亚洲人成网站在线播放欧美日韩| 久久久国产精品麻豆| 久久伊人香网站| 亚洲 国产 在线| 久久久久国内视频| 少妇熟女aⅴ在线视频| 国产三级黄色录像| 人妻丰满熟妇av一区二区三区| 久久狼人影院| 国产亚洲欧美98| 欧美一区二区精品小视频在线| 黄色视频,在线免费观看| 午夜久久久在线观看| av天堂在线播放| 91国产中文字幕| 久久精品亚洲熟妇少妇任你| 成人免费观看视频高清| 日日摸夜夜添夜夜添小说| 精品国产一区二区三区四区第35| 最新在线观看一区二区三区| 精品人妻1区二区| 国产精华一区二区三区| 精品少妇一区二区三区视频日本电影| 久久久久久久久久久久大奶| 欧美绝顶高潮抽搐喷水| 国产三级黄色录像| 看免费av毛片| 亚洲,欧美精品.| 午夜久久久在线观看| 亚洲少妇的诱惑av| 久久久精品国产亚洲av高清涩受| 亚洲 国产 在线| 午夜福利视频1000在线观看 | 青草久久国产| 欧美成人午夜精品| 久久天躁狠狠躁夜夜2o2o| 国产免费男女视频| 久久精品人人爽人人爽视色| 日韩大尺度精品在线看网址 | 91精品国产国语对白视频| 国产私拍福利视频在线观看| 色综合欧美亚洲国产小说| 亚洲全国av大片| 大型av网站在线播放| 免费一级毛片在线播放高清视频 | 亚洲av片天天在线观看| 久久久水蜜桃国产精品网| 久久久久精品国产欧美久久久| 午夜福利18| 国产精华一区二区三区| 在线观看66精品国产| 日韩欧美免费精品| 好看av亚洲va欧美ⅴa在| 琪琪午夜伦伦电影理论片6080| 亚洲天堂国产精品一区在线| 国产高清激情床上av| a在线观看视频网站| 国产成人免费无遮挡视频| 亚洲欧美日韩另类电影网站| 亚洲性夜色夜夜综合| 嫁个100分男人电影在线观看| 日韩欧美在线二视频| 侵犯人妻中文字幕一二三四区| 日本精品一区二区三区蜜桃| 精品国产一区二区久久| 国产精品乱码一区二三区的特点 | 69精品国产乱码久久久| 国产精品亚洲美女久久久| 中文字幕人成人乱码亚洲影| 亚洲五月天丁香| 一级毛片女人18水好多| 每晚都被弄得嗷嗷叫到高潮| 女性生殖器流出的白浆| 精品日产1卡2卡| 欧美成人午夜精品| 一二三四社区在线视频社区8| 视频在线观看一区二区三区| 久久热在线av| 纯流量卡能插随身wifi吗| 91成年电影在线观看| 首页视频小说图片口味搜索| 午夜激情av网站| 丰满的人妻完整版| 午夜亚洲福利在线播放| 色综合欧美亚洲国产小说| 欧美乱色亚洲激情| 婷婷精品国产亚洲av在线| 亚洲aⅴ乱码一区二区在线播放 | 男女之事视频高清在线观看| 久久午夜综合久久蜜桃| 亚洲国产欧美日韩在线播放| 国产精品九九99| 亚洲成a人片在线一区二区| 午夜两性在线视频| 亚洲中文av在线| 丰满人妻熟妇乱又伦精品不卡| 999久久久国产精品视频| 色播亚洲综合网| 国产成+人综合+亚洲专区| 精品久久久久久久毛片微露脸| 中文字幕av电影在线播放| 99在线人妻在线中文字幕| 久久热在线av| www.999成人在线观看| 欧美一级a爱片免费观看看 | 国产精品久久久久久精品电影 | 777久久人妻少妇嫩草av网站| 51午夜福利影视在线观看| 一本综合久久免费| 亚洲av第一区精品v没综合| 国产免费男女视频| 亚洲av熟女| 天天一区二区日本电影三级 | 亚洲精品中文字幕在线视频| 麻豆国产av国片精品| 午夜视频精品福利| 99精品久久久久人妻精品| 老司机靠b影院| 久9热在线精品视频| 色哟哟哟哟哟哟| 国产熟女午夜一区二区三区| 午夜福利高清视频| 日本一区二区免费在线视频| 久久人人97超碰香蕉20202| 一级作爱视频免费观看| 少妇粗大呻吟视频| 美女 人体艺术 gogo| 久久久久亚洲av毛片大全| 精品一区二区三区四区五区乱码| 国产精品久久久久久亚洲av鲁大| 夜夜夜夜夜久久久久| av天堂在线播放| 国产麻豆69| ponron亚洲| 真人一进一出gif抽搐免费| 国产一区二区三区视频了| 亚洲一卡2卡3卡4卡5卡精品中文| 嫩草影视91久久| 桃色一区二区三区在线观看| 午夜福利18| 精品欧美一区二区三区在线| 欧美在线黄色| 一卡2卡三卡四卡精品乱码亚洲| 亚洲第一av免费看| 亚洲欧美精品综合久久99| 中文字幕人成人乱码亚洲影| 丰满人妻熟妇乱又伦精品不卡| 国产91精品成人一区二区三区| 成人18禁在线播放| 欧美老熟妇乱子伦牲交| 亚洲在线自拍视频| 深夜精品福利| 夜夜躁狠狠躁天天躁| 啦啦啦 在线观看视频| 国产99久久九九免费精品| av有码第一页| 好看av亚洲va欧美ⅴa在| 久久精品国产亚洲av高清一级| 亚洲一区中文字幕在线| 久久久久国内视频| 一区二区三区高清视频在线| aaaaa片日本免费| 色在线成人网| 国产精品永久免费网站| 国产精品免费一区二区三区在线| 国产成人av激情在线播放| 久久久久久国产a免费观看| 丁香六月欧美| 亚洲狠狠婷婷综合久久图片| 久久精品亚洲熟妇少妇任你| 色播亚洲综合网| 777久久人妻少妇嫩草av网站| 一级片免费观看大全| 又紧又爽又黄一区二区| 中文字幕人妻熟女乱码| 最好的美女福利视频网| 91麻豆av在线| 91麻豆精品激情在线观看国产| 黑人欧美特级aaaaaa片| 欧洲精品卡2卡3卡4卡5卡区| 国产精品免费视频内射| 欧美乱妇无乱码| 国产欧美日韩一区二区三区在线| 桃色一区二区三区在线观看| 国产主播在线观看一区二区| 欧美久久黑人一区二区| 国产精品日韩av在线免费观看 | 色精品久久人妻99蜜桃| 精品国产亚洲在线| 一进一出抽搐动态| 丁香欧美五月| 正在播放国产对白刺激| 久久精品影院6| 色播在线永久视频| 51午夜福利影视在线观看| 国产精品99久久99久久久不卡| 国产亚洲av嫩草精品影院| 一进一出抽搐gif免费好疼| 国产麻豆成人av免费视频| 乱人伦中国视频| 国产精品亚洲av一区麻豆| 亚洲国产精品999在线| 久久久久国产一级毛片高清牌| 国产不卡一卡二| 欧美中文综合在线视频| 男女床上黄色一级片免费看| 亚洲av片天天在线观看| 亚洲国产欧美一区二区综合| 久久精品91蜜桃| 真人一进一出gif抽搐免费| 日韩免费av在线播放| 国内精品久久久久久久电影| 亚洲美女黄片视频| 熟女少妇亚洲综合色aaa.| 欧美午夜高清在线| 免费一级毛片在线播放高清视频 | 激情视频va一区二区三区| 在线观看免费午夜福利视频| 亚洲欧美日韩无卡精品| 国产成人啪精品午夜网站| 亚洲av美国av| 色综合婷婷激情| 满18在线观看网站| 激情视频va一区二区三区| 精品熟女少妇八av免费久了| 国产av在哪里看| 少妇的丰满在线观看| 一a级毛片在线观看| 亚洲精品国产色婷婷电影| 老汉色av国产亚洲站长工具| 成人永久免费在线观看视频| 久久精品成人免费网站| 亚洲全国av大片| 欧美丝袜亚洲另类 | 国产av一区在线观看免费| 亚洲精品美女久久久久99蜜臀| 99久久99久久久精品蜜桃| 精品免费久久久久久久清纯| 午夜两性在线视频| 国产国语露脸激情在线看| 村上凉子中文字幕在线| √禁漫天堂资源中文www| 天堂动漫精品| 久久精品亚洲熟妇少妇任你| 久久性视频一级片| 日韩精品青青久久久久久| 亚洲第一欧美日韩一区二区三区| 两个人免费观看高清视频| 国产成人av激情在线播放| 亚洲无线在线观看| 9热在线视频观看99| 真人一进一出gif抽搐免费| 国产男靠女视频免费网站| 法律面前人人平等表现在哪些方面| 在线免费观看的www视频| 窝窝影院91人妻| 制服丝袜大香蕉在线| 变态另类成人亚洲欧美熟女 | 久久狼人影院| 亚洲在线自拍视频| 黄色视频,在线免费观看| 亚洲色图av天堂| tocl精华| 久久人人精品亚洲av| 人妻久久中文字幕网| 中文字幕人妻丝袜一区二区| 超碰成人久久| 如日韩欧美国产精品一区二区三区| 最好的美女福利视频网| 欧美国产日韩亚洲一区| 亚洲精品美女久久av网站| 两性午夜刺激爽爽歪歪视频在线观看 | 成人三级黄色视频| 国产精品一区二区免费欧美| 欧美日韩亚洲综合一区二区三区_| 久久国产精品影院| 亚洲成a人片在线一区二区| 丰满的人妻完整版| 亚洲午夜理论影院| 嫁个100分男人电影在线观看| 久久久久亚洲av毛片大全| 色尼玛亚洲综合影院| 国产亚洲欧美在线一区二区| 激情在线观看视频在线高清| 国产国语露脸激情在线看| 午夜a级毛片| 一级毛片精品| 久久天躁狠狠躁夜夜2o2o| 日韩精品青青久久久久久| 欧美国产日韩亚洲一区| 大陆偷拍与自拍| 一边摸一边做爽爽视频免费| 99久久精品国产亚洲精品| 婷婷六月久久综合丁香| 久久 成人 亚洲| 色精品久久人妻99蜜桃| 国产精品影院久久| 黑人欧美特级aaaaaa片| 性色av乱码一区二区三区2| 免费在线观看影片大全网站| 欧美成人免费av一区二区三区| 黄色视频不卡| 老司机在亚洲福利影院| 亚洲人成伊人成综合网2020| 精品国产美女av久久久久小说| 精品国产国语对白av| 国产亚洲av嫩草精品影院| 亚洲欧美日韩无卡精品| 国内精品久久久久精免费| 日韩av在线大香蕉| 亚洲 国产 在线| av天堂久久9| 欧美激情高清一区二区三区| 亚洲中文字幕一区二区三区有码在线看 | 丝袜在线中文字幕| 老司机福利观看| av有码第一页| 午夜两性在线视频| 亚洲成国产人片在线观看| 日韩三级视频一区二区三区| 97人妻天天添夜夜摸| 亚洲人成伊人成综合网2020| 亚洲男人天堂网一区| 亚洲va日本ⅴa欧美va伊人久久| 久99久视频精品免费| 日韩欧美一区视频在线观看| 男男h啪啪无遮挡| 一区二区三区精品91| 丁香欧美五月| 黄色毛片三级朝国网站| 18禁裸乳无遮挡免费网站照片 | 色综合亚洲欧美另类图片| 亚洲国产日韩欧美精品在线观看 | 久久 成人 亚洲| 午夜福利在线观看吧| 国产亚洲精品第一综合不卡| 亚洲 欧美一区二区三区| 一个人观看的视频www高清免费观看 | 成人免费观看视频高清| 黑人巨大精品欧美一区二区mp4| 久久精品国产99精品国产亚洲性色 | 国产精华一区二区三区| 亚洲第一av免费看| 国产1区2区3区精品| 一区二区三区精品91| 国产一区在线观看成人免费|