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

    Nonlinear Model for Vibration Analysis of Fluid-conveying Pipes via the Incremental Harmonic Balance Method

    2011-06-22 05:07:10MENGDanCHENLiang
    船舶力學(xué) 2011年12期
    關(guān)鍵詞:理工大學(xué)增量諧波

    MENG Dan,CHEN Liang

    (1 School of Civil Engineering,Qingdao Technological University,Qingdao 266033,China;2 School of Architectural,Qingdao Technological University,Qingdao 266033,China;3 Architectural Design and Research Institute,Qingdao Technological University,Qingdao 266033,China)

    Nonlinear Model for Vibration Analysis of Fluid-conveying Pipes via the Incremental Harmonic Balance Method

    MENG Dan1,CHEN Liang2,3

    (1 School of Civil Engineering,Qingdao Technological University,Qingdao 266033,China;2 School of Architectural,Qingdao Technological University,Qingdao 266033,China;3 Architectural Design and Research Institute,Qingdao Technological University,Qingdao 266033,China)

    A new nonlinear model of a fluid-conveying pipe is presented.Considering the effect of the internal fluid and excitation,the nonlinear differential equations are derived based on Kane’s equation in combination with a Ritz method by using non-linear Lagrange strain theory.The equations of motion are linearized in the neighborhood of the equilibrium position to carry on the multimode linear stability analysis.In addition,the time histories for the displacements are obtained by using the Incremental Harmonic Balance(IHB)method.The validity of the new model is substantiated by comparing the model and results with those proposed by Pa?doussis(1998).

    nonlinear model;Kane’s method;vibration analysis;Incremental Harmonic Balance(IHB)method

    Biography:MENG Dan(1980-),male,Ph.D.,lecturer of Qingdao Techndogical University,E-mail:md1101@163.com.

    1 Introduction

    Vibrational dynamics of pipes conveying fluid has been the subject of numerous scientific investigations due to their wide applications.During the period from 1970 to 2000,based on the linear equations of motion for pipes supported at both ends proposed by Feodos’ev and Housner in 1950s,the literature on this topic literally exploded,with linear and nonlinear,twoand three-dimensional,theoretical and experimental studies,not only for the plain pipe and pipes with added spring supports,masses and dashpots,somewhere along the length of the pipe,but also for pipes supported at both ends.This was the result of gradual realization that this problem is a new paradigm in dynamics,displaying in a simple system a wide variety of the dynamical features of more complex engineering problems.Moreover,the fluid-conveying pipe system has the advantage that one can fairly easily conduct experiments to verify the theory.A thorough literature review will not be undertaken in this paper.For that,the reader is referred to Ref.[1-4].Most of the previous nonlinear equations were derived from the infinitesimal strain theory either by Newton approach[4]or by Hamilton’s principle[5].

    In the study of nonlinear systems,it is often of interest to find the time response for specific initial conditions.In the case of systems that can be represented as first order systems,excellent codes are nowadays available from netlib,such as LSODE written by Hindmarch[6]or those by Hairer et al[7].They enable the understanding of the dynamics of nonlinear systems and the classification of the types of response.The analyst,however,does not have that much of a choice when implicit second order ordinary differential equations are to be solved.The most popular method that can be used for finding periodic solutions of implicit second order ordinary differential equations is known as the Harmonic Balance(HB)method in the vibration literature[8],and as the Describing Function(DF)approach in the control and aeroelasticity literature[9].It is computationally very efficient for obtaining steady state solutions of nonlinear dynamics problems,but it is usually too restrictive,mainly because:the solution is assumed to have one dominant frequency component-or at most two-and the non-linearity has to be small;and one needs to know a great deal about the solution a priori or to carry enough terms in the solution and check the order of the coefficients of all the neglected harmonics,since otherwise one might obtain an inaccurate approximation,as shown by Nayfeh and Mook[8].

    Since the fluid-conveying pipes is the system which have beam-like shapes and undergo excitation,in this paper we will start from a new angle of view to address this problem.In order to drive equations of motion including the influence of internal fluid and excitation,the Kane’s method is adopted,and the effect of nonlinearity is taken into account with the help of the general Lagrange strain theory instead of the conventional infinitesimal strain theory.Utilizing the Rayleigh-Ritz method,the discretized equation of motion is obtained.This is an alternative approach to the existing methods utilizing Newtonian balance equations or Hamilton’s principle is thus provided.Then the equations of motion are linearized in the neighborhood of the equilibrium position,and solved by the complex damping method for multimode linear stability analysis.The results are compared with existing works[1].Furthermore,the resulting equations with continuous cubic nonlinearities are solved using the Incremental Harmonic Balance(IHB)method[10],the histories of the displacement are computed,and the results are also compared with existing works.

    2 Equations of motion

    A simply-supported slender shape pipe conveying fluid attached to a rigid base A and the configurations before and after deformation are exhibited by Fig.1.The in-plane deformation u→of the pipe can be represented by two scalar variables in 2-D space.Conventionally,Cartesian variables u1and u2,as shown in the figure,are widely utilized.In the present modelling method,a non-Cartesian variable along with one Cartesian variable is used instead of two Cartesian variables:s denotes the arc length stretch,rather than u1,and measures the Cartesian distance of a generic point in the axial direction of the undeformed configuration.O is a reference point identifying a point fixed in a floating frame in the system,while r→is a po-sition vector from O to the location of P in the undeformed body.

    Fig.1 Schematics of a simply-supported straight pipe conveying fluid:(a)model of pipe before deformation;(b)the pipe after deformation and kinematics subject to excitation;(c)the cross-section of the pipe

    In this study,Kane’s method is employed to derive the equations of motion.This method provides a direct and systematic way of deriving ordinary differential equations of motion for continua such as beams,plates and shells.With the assumptions given above,the equations of motion for pipes can be obtained as follows[11]:

    here,the expression of generalized inertia forceand generalized force Fiare

    in which ρ and ρfare the mass per unit length of the pipe and internal fluid,respectively,are the time derivatives of the generalized co-ordinates.,andare the inertia velocity and accelerations of the center point P of a differential element of the pipe and internal fluid,respectively.

    For the present modelling method,the following strain energy expression is used,in which Tais the apparent tension and M2is the bending moment[11-14].

    In the pipe’s wide applications,some of them will be employed in the water.In this case the pipe oscillates transversely to a stationary and uniform flow of free stream velocity,in which the inertia force and the lift forces are exerted on the system as the external force.The expression for the inertia force and the lift forces are exerted on the system H used in this section,reads[15]

    in which FL= ρeU2DCL/2 is the lift force for VIV,ρ′=CaρeπD2/4 is the added mass per unit length of the pipe,C′=δωsρeD2is the equivalent fluid damping.In which Cais added mass coefficient,taken as 1.0;ρeis the density of sea water,D is the outer diameter of the pipe,δ is the viscous force coefficient,taken as 0.8;ωs=2πStU/D is the Strouhal frequency,Stis the Strouhal number,taken as 0.2;U is velocity of currents and waves,CLis the vortex lift coefficient.

    Suppose that the vertical excitation which the motion of the pipe subjected to is dictated by a excitation function as

    in which Z0and ? are the acceleration amplitude and angular frequency of the harmonic vibration,respectively.Thus,substituting Eqs.(3)-(9)into Eq.(1)we can get the nonlinear equations of transverse vibration for the pipe conveying fluid[16]

    It is interesting to compare and contrast these equations with those previously published.The equations proposed in this paper are compared with those gained by Pa?doussis and his co-workers[1,17].Using the same notations of this paper,the equations of Pa?doussis and his coworkers can be rewritten as

    The nonlinear terms of this model will not be detailed out here,the specific content of nonlinear terms can be found in the literature of Pa?doussis and his co-workers easily.Neglect the external fluid and excitation,the main difference of Eq.(10)from Eq.(11)is shown in nonlinear terms,Fihkj,Sihkj,Nihkjand Dihkj,which have different form.But the linear terms are the same in both equations of motion.

    3 Multimode linear stability analyses

    The external force and the vertical excitation are not considered,introducing the dimensionless quantities

    the equation of motion can be expressed in the following dimensionless form

    The discretized equations of motion can be verified by computing the natural frequencies of the pipe.To this end,it is necessary to obtain the linearized equations of motion in the neighborhood of the equilibrium position.When the excitation and the external force are neglected,the linearized equation of Eq.(13)is rendered by

    In order to find the natural frequencies from Eq.(14),it is convenient to transform Eq.(14)into

    Assuming the solution of Eq.(16)as Y(t)=Y,the complex eigenvalues λndepend on the internal fluid velocity because of the matrix[C]and[K].The critical speeds of the internal fluid are the values of V at which the imaginary part of any of the minimum eigenvalues to be zero.The real part of the eigenvalues is related to the amount of damping exhibited by a given eigenvalue.

    Fig.2 Comparison of results of Pa?doussis and the present formulation for:(a)ζ=0.1;(b)ζ=0.8

    Pa?doussis probed into the static and dynamic instabilities of pipes conveying fluids[1].Therefore,the results for the present model will be validated with those of Pa?doussis under two different values of non-dimensional parameter ζ(ζ=0.1 and ζ=0.8)with the clampedclamped boundary condition.The numerical results are shown in Fig.2,which exhibits the relationship between the dimensionless flow velocity and the dimensionless frequencyWhereis the frequency of the structure in rad/s,whileVl.For ζ=0.1,the divergence occurs first when the flow velocity reaches 6.28,and if the flow velocity is increased further,the second mode divergence occurs at flow velocity of 8.99.Upon further increasing the flow velocity,the first and second mode frequencies increase together and join with the third mode.The first mode then decreases to zero and the second and third modes move together to the fourth one.In the case of ζ=0.8,there is a stabilized region between u=8.99 and 9.6.From both figures,it is clear that the present formulation provides excellent agreement with those of Pa?doussis[1].It must be pointed out that not only do the curves match well but also the various physical phenomena that occur in those results have been correctly predicted by the current formulation.

    4 Incremental Harmonic Balance(IHB)method

    The purpose of this section is to describe the method that can be used for finding periodic solutions of nonlinear oscillators of the type of Eq.(13).Contrary to small parameter techniques that require the nonlinear terms to be small,it is preferable to have a scheme that is applicable in a more general case.Urabe[18]was one of the first to introduce a multi-harmonic balancing procedure.He combined it with the Newton-Raphson iteration technique to treat strong nonlinearities and he referred to this as a Galerkin approximation with a characteristic function containing trigonometric functions.Similarly,the Incremental Harmonic Balance(IHB)method was developed by Lau et al[19]to take into account multiple harmonic components due to the nonlinear forces.Hence,in the usual terminology,the IHB can be referred to as a combined Incremental Ritz-Galerkin Harmonic Balance method or as a Harmonic Balance Newton-Raphson(HBNR)method.Indeed,Ferri[20]showed that the IHB and the HBNR methods are equivalent.The IHB method was successfully applied to various types of nonlinear structural systems,ranging from systems with continuous cubic nonlinearities to dry friction damper problems where the nonlinear force is discontinuous[21].From the assumption that the value of n2for 5,with βihkj=βihjkand γihkj=γihjkin mind,Eq.(13)can be rewritten as

    And now all derivatives are now with respect to τ.The first step in applying the IHB method is similar to a Newton-Raphson procedure:Starting from an initial guess of the solution of Eq.(18)is(yk0,?0,fi0),a neighbouring solution is reach through an incrementation(Δyk,Δ?,Δfi)using a Taylor series expansion.Neglecting the non-linear terms in Δyk,Δ? and Δfi,a linearized incrementation equation is obtained:

    The second step in the IHB method is the process of harmonic balance.Starting from known or estimated coefficients akmand bkm,the assumed or sought periodic solution is written as

    By using the sought periodic solution,it is straightforward to obtain

    By differentiation of Eq.(22),one obtainsandas well.Substituting Eq.(22)and,into Eq.(19)and make use of the harmonic balance method in one period yields

    It can be seen that the number of unknowns exceeds the number of equations by two,Δ? and Δfi.Hence,one needs to add a constraint.Recalling that the solution sought is for a given acceleration amplitude or of a predefined frequency,the constraint is either Δfi=0 or Δ?=0.If Δfi=0,fiis said to be the active increment and the solution for a different fiis obtained by incrementing fifrom point to point.Similarly,if Δ?=0,? is the active increment.Once the incremental coefficient ΔAkis found,the assumed solution can be updated and the process is repeated until a convergence criterion is met.This criterion may be defined as the infinite norm of the incremental vecter having to be smaller than a certain value.

    5 Results and discussion

    The dynamic responses with transverse displacement being our primary concern are computed with same parameters and material properties except that the internal fluid velocity is prescribed as,V=0ms-1and V=131.3ms-1.And the comparisons of dynamic responses are made between this new model and Pa?doussis’ model[1,4].Both of them use the same basis functions with the clamped-clamped boundary condition.In order to make up the vortex lift,the external transverse load q herein is dictated by a harmonic excitation with frequency and constant amplitude being ?fand P0,respectively,namely,q ()t=P0sin?ft.The transverse displacements at x=L/2 of the nonlinear models are plotted in Figs.3 and 4.Figures show that,in this case,although the new model and Pa?doussis’ model have different nonlinear terms,as seen in Eqs.(10)and(11),there is no difference between the dynamic responses obtained from these two models.

    Fig.3 Dynamic responses of the transverse displacement for the non-linear model when the internal fluid velocity is V=0m/s:(a)the present model;(b)Pa?doussis’ model

    Fig.4 Dynamic responses of the transverse displacement for the non-linear model when the internal fluid velocity is V=131.3m/s:(a)the present model;(b)Pa?doussis’ model

    From the Morrison wave theory,the hydrodynamic drag forces are involves two key factors:drag force and inertial mass coefficients.The drag force coefficients describe the pipes’surface roughness,which are affected by Reynolds Number.Over a time,converging onto a steady state as depicted in Fig.11.In this study,only the steady flow is taken into consideration,without considering the impact of lift force,the dynamic reaction of the surrounding water on the pipe is simplified as the equivalent fluid damping,according to Blevins[15],which is adopted as the damping force on the pipe.In the present model,the equivalent fluid damping is embodied in the second term of the equation.The dynamic responses of the transverse displacement as a function of the damping force caused by external fluid velocity are discussed and the results are presented in Fig.5.In this case,the black line represents the vibration state of the pipe without the impact of the fluid damping,which corresponds to the result of Pa?doussis’ model.But we can see from the figure,the drag forces play a significant role in reducing the maximum amplitudes of dynamic responses,as the external fluid velocity increases,as expected,the amplitude of vibration decreases accordingly.The cause lies in the fluid damping term,which changes with the external fluid velocity.

    Fig.5 Dynamic responses of the transverse displacement as a function of external fluid velocity

    When the structure is subjected to earthquake,its dynamic characteristic will change a lot.Then the dynamic responses of the simply-supported pipe with vertical excitation prescribed by Z(t)=Z0sin?t are discussed next.Some case studies are carried out to analyze the response characteristics.Figs.6-7 depict the histories of the displacement of the pipe under different conditions.The internal flow speeds from 0 to 230m/s are considered in the parametric studies.Losses due to internal flow frictions in the pipe system clearly deter the internal flow speed and are not considered in this investigation.The dynamic responses of the pipes transporting the internal fluids are illustrated in Fig.6.The results clearly exhibit that the displacement responses of the pipes with bending rigidity tend to converge onto steady state response amplitudes.And it demonstrates that as the internal flow speed increases,the response amplitude has decreased.

    Fig.6 Dynamic responses of the transverse displacement as a function of internal fluid velocity when the excitation is given by Z0=1.0m/s2,?=24πrad/s:(a)V=0m/s;(b)V=65.7m/s;(c)V=164.2m/s;(d)V=230.0m/s

    Through the introduction and analysis of a simple discrete system,the important modelling issues motion-induced stiffness variations,and their contribution to non-linear equations of motion have been investigated.Fig.7 demonstrates the responses of the transverse displacement at mid-point of the pipe with respect to combined internal fluid and overall motion ex citations,which manifests that the increasing angular frequency of overall motion reduces the response magnitude,while the impact of the internal fluid velocity is not significant,as expected.If one chooses a special set of deformation variables such that the strain energy of the structure is accurately characterized in a quadratic form in these variables,then it is possible to develop differential equations of motion which are totally linear in the deformation variables and yet still exhibit proper stiffness variations.This is the approach taken in the modeling method developed in this paper,wherein the strain energy of the 3-D beam is expressed in a quadratic form by using the stretch variable of the beam.If one chooses the Cartesian set of deformation variables,the stiffness variation due to large overall motion cannot be properly predicted by a linear model.If proper motion-induced stiffness variation is to be accommodated,at least up to second degree terms in the stretch strain should be retained.

    Fig.7 Dynamic responses of the transverse displacement as a function of angular frequency when the external fluid velocity is U=1m/s and Z0=1.0m/s2:(a)?=2π;(b)?=6π;(c)?=10π;(d)?=24π

    6 Conclusions

    In this paper,the nonlinear equations of 2-D motion of pipe conveying incompressible fluid have been derived based on the Kane’s method,they are quite similar to those obtained by Pa?doussis(1998).But in their most general form the present equations account for the possible existence of:(1)a stationary and uniform external flow of free stream velocity,and(2)vertical excitation.The equations of motion are solved by the complex damping method for the natural frequency.The Incremental Harmonic Balance method was proposed to predict the existence of periodic solutions.As suggested in the literature,it is capable of solving problems defined by implicit equations.Furthermore,the numerical results from the proposed model are compared with those from Pa?doussis’ model.

    It is still noteworthy that the approach developed here can also be used to handle arbitrary boundary conditions.Adaptation of the model for a specific set of geometric boundary conditions requires only the selection of a suitable set of basic functions.

    [1]Pa?doussis M P.Fluid-structure interactions:Slender structures and axial flow[M].Academic Press,London,1998.

    [2]Meng Dan,Guo Haiyan,Xu Sipeng.Nonlinear dynamic model of a fluid-conveying pipe undergoing overall motions[J].Appl.Math.Model,2011,35:781-796.

    [3]Meng Dan,Guo Haiyan,Xu Sipeng.Flow-induced vibration stability analysis of fluid-conveying pipes[J].J Vib.and Shock,2010,29(6):86-89.(in Chinese)

    [4]Semler C,Li G X,Pa?doussis M P.The non-linear equations of motion of pipes conveying fluid[J].J Sound Vib.,1994,169:577-599.

    [5]Folley C N,Bajaj A K.Spatial nonlinear dynamics near principal parametric resonance for a fluid-conveying cantilever pipe[J].J Fluid Struct.,2005,21:459-484.

    [6]Hindmarsh A C.Scientific Computing,North-Holland[M].ODEPACK,a systematized collection of ODE solvers,Amsterdam,1983.

    [7]Hairer E,Norsett S P,Wanner G.Solving ordinary differential equations[M].Springer-Verlag,Berlin,1993.

    [8]Nayfeh A H,Mook D T.Nonlinear Oscillations[M].Wiley-Interscience,New York,1979.

    [9]Gelb A,Vander Velde W E.Multiple-input describing functions and nonlinear system design[M].McGraw-Hill,New York,1968.

    [10]Semler C,Gentleman W C,Pa?doussis M P.Numerical solutions of second order implicit non-linear ordinary differential equations[J].J Sound Vib.,1996,2:553-574.

    [11]Kane T R,Ryan R R,Banerjee A K.Dynamics of a cantilever beam attached to a moving base[J].J Guid.Control Dynam.,1987,10(2):139-151.

    [12]Yoo H H,Ryan R R,Banerjee A K.Dynamics of flexible beams undergoing overall motion[J].J Sound Vib.,1995,181(2):261-278.

    [13]Timoshenko S P,Gere J M.Mechanics of Materials[M].Litton Educational Publishing,New York,1972.

    [14]Chucheepsakul S,Monprapussorn T,Huang T.Large strain formulation of extensible flexible marine pipes transporting fluid[J].J Fluids Struct.,2003,17:185-224.

    [15]Blevins R D.Flow-induced vibrations[M].2nd edn.Van Nostrand&Co.,New York,1990.

    [16]Meng Dan,Guo Haiyan,Xu Sipeng.A new nonlinear model for vibration analysis of fluid-conveying pipes undergoing overall motions[C].Proceedings of Second International Conference on Modelling and Simulation,2009,7:248-253.

    [17]Wadham-Gagnon M,Pa?doussis M P,Semler C.Dynamics of cantilevered pipes conveying fluid.Part 1:Nonlinear equations of three-dimensional motion[J].J Fluid Struct.,2007,23:545-567.

    [18]Urabe M.Galerkins procedure for nonlinear periodic systems[J].Archive for Rational Mechanics and Analysis,1965,20:150-152.

    [19]Lau S L,Cheung Y K,Wu S Y.A variable parameter incrementation method for dynamic instability of linear and nonlinear vibration of elastic systems[J].Journal of Applied Mechanics,1982,49:849-853.

    [20]Ferri A A.On the equivalence of the incremental harmonic balance method and the harmonic balance Newton_Raphson method[J].J Applied Mech.,1986,53:455-457.

    [21]Pierre C,Ferri A A,Dowell E H.Multi-harmonic analysis of dry friction damped systems using an incremental harmonic balance method[J].J Applied Mech.,1985,52:958-964.

    基于增量諧波平衡法的輸流管道非線性振動分析

    孟 丹1, 陳 亮2,3
    (1青島理工大學(xué) 土木工程學(xué)院,山東 青島266033;2青島理工大學(xué) 建筑學(xué)院,山東 青島266033;3青島理工大學(xué) 建筑設(shè)計(jì)研究院,山東 青島 266033)

    考慮內(nèi)、外流體的影響,利用Kane方法和拉格朗日應(yīng)變理論建立了輸流管道的二維非線性動力學(xué)模型。將動力學(xué)模型在平衡位置附近線性化,進(jìn)行輸流管道的線性穩(wěn)定性分析。采用增量諧波平衡法求解包含三次非線性項(xiàng)動力學(xué)模型的穩(wěn)態(tài)周期解。探討了不同參數(shù)影響下輸流管道的非線性動力響應(yīng)以及內(nèi)流流速的大小對管道動力特性的影響。

    非線性模型;Kane方法;振動分析;增量諧波平衡法

    O353.1

    A

    孟 丹(1980-),男,博士,青島理工大學(xué)土木工程學(xué)院講師,從事海洋工程非線性動力分析方面的研究;

    陳 亮(1980-),男,碩士生。

    O353.1

    A

    1007-7294(2011)12-1416-13

    date:2011-07-07

    Supported by National Natural Science Foundation of China(50808105);Shandong Province Natural Science Foundation(BS2009SF003)

    猜你喜歡
    理工大學(xué)增量諧波
    昆明理工大學(xué)
    提質(zhì)和增量之間的“辯證”
    “價(jià)增量減”型應(yīng)用題點(diǎn)撥
    昆明理工大學(xué)
    昆明理工大學(xué)
    浙江理工大學(xué)
    基于均衡增量近鄰查詢的位置隱私保護(hù)方法
    虛擬諧波阻抗的并網(wǎng)逆變器諧波抑制方法
    基于ELM的電力系統(tǒng)諧波阻抗估計(jì)
    基于ICA和MI的諧波源識別研究
    電測與儀表(2015年8期)2015-04-09 11:50:00
    成人三级黄色视频| 国产精品一区二区三区四区免费观看| 一级av片app| 性色avwww在线观看| 一区二区三区免费毛片| 国模一区二区三区四区视频| 男女视频在线观看网站免费| 观看美女的网站| 日韩av在线大香蕉| 久久6这里有精品| 欧美激情在线99| 久久午夜福利片| 亚洲国产精品成人久久小说 | 三级国产精品欧美在线观看| www.色视频.com| 欧美激情在线99| 老司机福利观看| 免费一级毛片在线播放高清视频| 99久久精品一区二区三区| 日韩亚洲欧美综合| 99在线视频只有这里精品首页| 精品无人区乱码1区二区| 女的被弄到高潮叫床怎么办| 成人亚洲精品av一区二区| 久久久久性生活片| 可以在线观看的亚洲视频| 久久久久性生活片| 人妻久久中文字幕网| 亚洲一区高清亚洲精品| 99久国产av精品| 别揉我奶头 嗯啊视频| 亚洲在线观看片| 少妇的逼水好多| 国产精品国产高清国产av| 白带黄色成豆腐渣| 欧美区成人在线视频| 菩萨蛮人人尽说江南好唐韦庄 | 日韩成人av中文字幕在线观看| 秋霞在线观看毛片| 日韩成人伦理影院| 欧美日韩精品成人综合77777| 深爱激情五月婷婷| 一卡2卡三卡四卡精品乱码亚洲| 国产探花极品一区二区| avwww免费| 女人十人毛片免费观看3o分钟| 久久精品国产99精品国产亚洲性色| av在线亚洲专区| 国产亚洲5aaaaa淫片| 久久久精品欧美日韩精品| 亚洲人与动物交配视频| 啦啦啦观看免费观看视频高清| 国产成人精品久久久久久| 美女黄网站色视频| 日韩 亚洲 欧美在线| 日韩视频在线欧美| 此物有八面人人有两片| 偷拍熟女少妇极品色| 麻豆精品久久久久久蜜桃| 人人妻人人看人人澡| 又粗又爽又猛毛片免费看| 成年女人永久免费观看视频| 久久人人爽人人片av| 欧美高清成人免费视频www| 久久久午夜欧美精品| 久久人人爽人人爽人人片va| 国产美女午夜福利| 国产伦理片在线播放av一区 | 麻豆国产av国片精品| 大又大粗又爽又黄少妇毛片口| 国产精品一区二区性色av| 亚洲av成人av| 欧美色欧美亚洲另类二区| 在线播放无遮挡| 国产av在哪里看| 亚洲av一区综合| 免费看日本二区| 岛国毛片在线播放| 免费不卡的大黄色大毛片视频在线观看 | 久久欧美精品欧美久久欧美| 日本免费a在线| 婷婷色av中文字幕| 一个人看的www免费观看视频| 热99在线观看视频| 18禁裸乳无遮挡免费网站照片| 国产亚洲av片在线观看秒播厂 | 成人综合一区亚洲| 日日啪夜夜撸| 日日摸夜夜添夜夜爱| 亚洲久久久久久中文字幕| 国产精品.久久久| 免费av毛片视频| 免费观看的影片在线观看| 亚洲精品乱码久久久v下载方式| 精品国内亚洲2022精品成人| 国内精品一区二区在线观看| 波多野结衣高清无吗| 男女做爰动态图高潮gif福利片| 波多野结衣高清作品| 亚洲av第一区精品v没综合| 美女cb高潮喷水在线观看| 青春草国产在线视频 | 边亲边吃奶的免费视频| 成人国产麻豆网| 日韩强制内射视频| 欧美3d第一页| 免费av不卡在线播放| 亚洲综合色惰| 久久久久网色| 国产伦在线观看视频一区| 国产一区二区在线观看日韩| 亚洲欧美清纯卡通| 欧美变态另类bdsm刘玥| 亚洲成人中文字幕在线播放| 成人无遮挡网站| 十八禁国产超污无遮挡网站| 欧美一区二区精品小视频在线| or卡值多少钱| 晚上一个人看的免费电影| 一个人看的www免费观看视频| 国产精品,欧美在线| 夜夜爽天天搞| 嘟嘟电影网在线观看| 亚洲精品乱码久久久v下载方式| 免费av观看视频| 熟女电影av网| 国产精品福利在线免费观看| 久久久久久久久大av| 久久久精品欧美日韩精品| 男插女下体视频免费在线播放| 亚洲欧美精品综合久久99| 成人一区二区视频在线观看| 人妻制服诱惑在线中文字幕| 91午夜精品亚洲一区二区三区| 亚洲美女视频黄频| 超碰av人人做人人爽久久| 国产精品久久久久久精品电影小说 | 国产高清有码在线观看视频| 久久久久网色| 精品久久久久久久久久久久久| 国产亚洲精品av在线| 成年版毛片免费区| 中文字幕av成人在线电影| 亚洲综合色惰| 黄色视频,在线免费观看| 午夜福利在线观看吧| 成年免费大片在线观看| 又黄又爽又刺激的免费视频.| 国产伦理片在线播放av一区 | av在线蜜桃| 国内揄拍国产精品人妻在线| 久久久久久久午夜电影| 一本久久精品| 日本三级黄在线观看| 少妇高潮的动态图| 国产亚洲5aaaaa淫片| 日韩av在线大香蕉| 亚洲在线观看片| 人人妻人人澡欧美一区二区| 69av精品久久久久久| 老女人水多毛片| 精品人妻偷拍中文字幕| 欧美成人a在线观看| 成人综合一区亚洲| 国产不卡一卡二| 国产精华一区二区三区| 熟妇人妻久久中文字幕3abv| 国产一区亚洲一区在线观看| 精品日产1卡2卡| 国产精品久久久久久久电影| 一级毛片久久久久久久久女| 最近最新中文字幕大全电影3| 男人舔奶头视频| 91久久精品国产一区二区成人| 久久人妻av系列| 久久欧美精品欧美久久欧美| 69人妻影院| 欧美成人精品欧美一级黄| 日韩视频在线欧美| 色5月婷婷丁香| 少妇人妻一区二区三区视频| 国产精品女同一区二区软件| 中文字幕制服av| 久久精品国产亚洲av涩爱 | 免费看日本二区| 少妇高潮的动态图| 99热网站在线观看| 日韩亚洲欧美综合| 黄色欧美视频在线观看| 1000部很黄的大片| 中国国产av一级| a级毛色黄片| 精品少妇黑人巨大在线播放 | 哪里可以看免费的av片| 欧美xxxx性猛交bbbb| 欧美+日韩+精品| 欧美日韩精品成人综合77777| 美女内射精品一级片tv| 中文在线观看免费www的网站| 亚洲经典国产精华液单| 日韩中字成人| 久久午夜福利片| 亚洲一区高清亚洲精品| 亚洲人成网站在线观看播放| 亚洲av中文av极速乱| 日韩成人av中文字幕在线观看| 亚洲av电影不卡..在线观看| 看免费成人av毛片| 黄片wwwwww| 99riav亚洲国产免费| 亚洲在线自拍视频| 男人舔奶头视频| 91久久精品国产一区二区三区| 精品无人区乱码1区二区| 亚洲国产精品合色在线| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲欧洲日产国产| 三级毛片av免费| 成人美女网站在线观看视频| 老熟妇乱子伦视频在线观看| 91麻豆精品激情在线观看国产| 一级毛片我不卡| 国产色爽女视频免费观看| 成人毛片a级毛片在线播放| 久久久精品94久久精品| 久久久国产成人精品二区| 伦理电影大哥的女人| 我的女老师完整版在线观看| 禁无遮挡网站| 亚洲三级黄色毛片| 亚洲成人久久性| 中文字幕熟女人妻在线| 国产69精品久久久久777片| 国产色爽女视频免费观看| 国产老妇女一区| 日韩制服骚丝袜av| 国模一区二区三区四区视频| 成人欧美大片| 中文欧美无线码| 精品免费久久久久久久清纯| 欧洲精品卡2卡3卡4卡5卡区| 人妻夜夜爽99麻豆av| 国产成人福利小说| 久久精品夜色国产| 亚洲av中文字字幕乱码综合| 欧美+亚洲+日韩+国产| 欧美一区二区亚洲| 91精品国产九色| 欧美激情久久久久久爽电影| 久久久午夜欧美精品| 淫秽高清视频在线观看| 亚洲一区高清亚洲精品| 亚洲国产精品合色在线| 热99在线观看视频| av视频在线观看入口| 午夜久久久久精精品| 国产精品久久久久久久电影| 少妇的逼好多水| 免费在线观看成人毛片| 黄色配什么色好看| 日韩三级伦理在线观看| 久久人人爽人人片av| 亚洲无线观看免费| avwww免费| 日日干狠狠操夜夜爽| 成年av动漫网址| av.在线天堂| a级毛片a级免费在线| 美女高潮的动态| .国产精品久久| 亚洲国产欧美在线一区| 青春草国产在线视频 | 免费观看a级毛片全部| 插阴视频在线观看视频| 日韩在线高清观看一区二区三区| 免费av毛片视频| 久久人人爽人人片av| 欧美激情国产日韩精品一区| 亚洲精品成人久久久久久| 中国国产av一级| 人妻系列 视频| 亚洲,欧美,日韩| 91aial.com中文字幕在线观看| 久久6这里有精品| 亚洲欧美日韩东京热| 国产三级在线视频| 国产高清有码在线观看视频| 国产久久久一区二区三区| 欧美三级亚洲精品| 日韩三级伦理在线观看| 精品人妻熟女av久视频| 国产av麻豆久久久久久久| 最好的美女福利视频网| 国产黄片视频在线免费观看| 精品人妻一区二区三区麻豆| 国产爱豆传媒在线观看| 亚洲七黄色美女视频| 日韩欧美国产在线观看| 久久精品人妻少妇| 国产精品一区www在线观看| 夜夜看夜夜爽夜夜摸| 夜夜爽天天搞| 国产精品.久久久| 在线观看美女被高潮喷水网站| 麻豆一二三区av精品| 亚洲av免费在线观看| 99久久人妻综合| av免费观看日本| 青春草国产在线视频 | 亚洲av第一区精品v没综合| 18+在线观看网站| 变态另类丝袜制服| 亚洲三级黄色毛片| 床上黄色一级片| 观看免费一级毛片| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲激情五月婷婷啪啪| 可以在线观看的亚洲视频| 精品一区二区三区人妻视频| 久久午夜亚洲精品久久| 亚洲人成网站在线播放欧美日韩| 18禁裸乳无遮挡免费网站照片| 亚洲色图av天堂| 亚洲av第一区精品v没综合| 日韩成人伦理影院| 免费大片18禁| 亚洲无线观看免费| 狂野欧美激情性xxxx在线观看| 99精品在免费线老司机午夜| 麻豆国产97在线/欧美| 亚洲欧美日韩高清在线视频| 99热网站在线观看| 久久久a久久爽久久v久久| www.色视频.com| 九九爱精品视频在线观看| av国产免费在线观看| 久久精品国产亚洲网站| 99国产极品粉嫩在线观看| 国产麻豆成人av免费视频| 精品一区二区三区人妻视频| 国产真实乱freesex| 国产一级毛片七仙女欲春2| 高清毛片免费观看视频网站| 高清午夜精品一区二区三区 | 日韩精品青青久久久久久| 99在线视频只有这里精品首页| 精品熟女少妇av免费看| 1024手机看黄色片| 男的添女的下面高潮视频| 国产老妇伦熟女老妇高清| 亚洲中文字幕日韩| 1024手机看黄色片| 久久热精品热| 美女 人体艺术 gogo| 91久久精品电影网| 国产成人影院久久av| 三级男女做爰猛烈吃奶摸视频| 亚洲aⅴ乱码一区二区在线播放| 成人高潮视频无遮挡免费网站| 亚洲av一区综合| 亚洲性久久影院| 深爱激情五月婷婷| 久久鲁丝午夜福利片| 在线观看美女被高潮喷水网站| 91狼人影院| 亚洲av一区综合| 欧美激情国产日韩精品一区| 久久6这里有精品| 欧美激情久久久久久爽电影| 好男人视频免费观看在线| 免费av不卡在线播放| 美女内射精品一级片tv| 在线天堂最新版资源| 久久久久九九精品影院| av在线蜜桃| 久久久国产成人免费| 两个人的视频大全免费| 在线免费观看不下载黄p国产| 男女啪啪激烈高潮av片| 成年免费大片在线观看| 国产一区二区在线观看日韩| 色尼玛亚洲综合影院| 九九久久精品国产亚洲av麻豆| 亚洲第一电影网av| 国产精品免费一区二区三区在线| 日韩大尺度精品在线看网址| 国产精品日韩av在线免费观看| 天堂av国产一区二区熟女人妻| 亚洲国产欧美人成| 亚洲性久久影院| 国产黄色小视频在线观看| 亚洲最大成人中文| 欧美xxxx黑人xx丫x性爽| 国产精品1区2区在线观看.| 天天躁日日操中文字幕| 欧美三级亚洲精品| 可以在线观看的亚洲视频| 欧美+亚洲+日韩+国产| 国产日本99.免费观看| 久久99精品国语久久久| 亚洲av不卡在线观看| 色噜噜av男人的天堂激情| 成人亚洲精品av一区二区| 精品国产三级普通话版| 精品人妻一区二区三区麻豆| 长腿黑丝高跟| 国产亚洲精品av在线| 在线观看美女被高潮喷水网站| 国产精品乱码一区二三区的特点| 18禁在线无遮挡免费观看视频| 女同久久另类99精品国产91| 99久久精品国产国产毛片| 1000部很黄的大片| a级毛色黄片| 亚洲成人中文字幕在线播放| 高清在线视频一区二区三区 | av在线天堂中文字幕| 欧美激情国产日韩精品一区| 美女内射精品一级片tv| 精品一区二区三区人妻视频| 久久人人爽人人片av| 国产精品无大码| 欧美另类亚洲清纯唯美| 97超视频在线观看视频| 精品久久久久久久人妻蜜臀av| 欧美性猛交╳xxx乱大交人| 国产成人freesex在线| 少妇猛男粗大的猛烈进出视频 | 青春草国产在线视频 | 欧美色欧美亚洲另类二区| 成人三级黄色视频| 高清在线视频一区二区三区 | 精品免费久久久久久久清纯| 三级男女做爰猛烈吃奶摸视频| 给我免费播放毛片高清在线观看| 亚洲精品自拍成人| 婷婷色av中文字幕| 99热网站在线观看| 99久久精品一区二区三区| av在线天堂中文字幕| 日韩,欧美,国产一区二区三区 | 国产色爽女视频免费观看| 长腿黑丝高跟| 成人永久免费在线观看视频| 99久久无色码亚洲精品果冻| 亚洲一区二区三区色噜噜| 国产黄色视频一区二区在线观看 | 欧美日韩乱码在线| 亚洲成人av在线免费| 国产亚洲精品久久久com| 天天躁夜夜躁狠狠久久av| 欧美色欧美亚洲另类二区| 欧美bdsm另类| 真实男女啪啪啪动态图| 九九热线精品视视频播放| 天堂av国产一区二区熟女人妻| 国产精品一区二区三区四区久久| 国产亚洲欧美98| 成人亚洲精品av一区二区| 在线观看免费视频日本深夜| 国产一区二区亚洲精品在线观看| 午夜久久久久精精品| 成人国产麻豆网| 一卡2卡三卡四卡精品乱码亚洲| 亚洲国产精品sss在线观看| 成人特级黄色片久久久久久久| 91午夜精品亚洲一区二区三区| 美女内射精品一级片tv| 免费一级毛片在线播放高清视频| 婷婷六月久久综合丁香| 美女 人体艺术 gogo| 成人亚洲欧美一区二区av| 日韩强制内射视频| 国产色爽女视频免费观看| 国产精品福利在线免费观看| 最近2019中文字幕mv第一页| 久久人人精品亚洲av| 晚上一个人看的免费电影| 亚洲美女视频黄频| 真实男女啪啪啪动态图| 国产成人福利小说| av在线亚洲专区| 91在线精品国自产拍蜜月| 国内精品一区二区在线观看| 日韩高清综合在线| 精品一区二区三区人妻视频| 韩国av在线不卡| 国产真实伦视频高清在线观看| 中出人妻视频一区二区| 精品一区二区免费观看| 中文字幕久久专区| 综合色丁香网| 黄色日韩在线| 欧美色视频一区免费| 又粗又爽又猛毛片免费看| 又爽又黄无遮挡网站| 久久精品夜夜夜夜夜久久蜜豆| 老师上课跳d突然被开到最大视频| 亚洲熟妇中文字幕五十中出| 免费av观看视频| 狠狠狠狠99中文字幕| 在线观看66精品国产| 在线观看午夜福利视频| 日韩一本色道免费dvd| 色视频www国产| 欧美bdsm另类| 成人欧美大片| 亚洲国产精品成人综合色| 两个人的视频大全免费| 女的被弄到高潮叫床怎么办| 欧美精品国产亚洲| 国产免费男女视频| 日本欧美国产在线视频| 国产 一区精品| 国产亚洲av片在线观看秒播厂 | 成人亚洲欧美一区二区av| 欧美xxxx性猛交bbbb| 成人国产麻豆网| 国产精品野战在线观看| 女人十人毛片免费观看3o分钟| 国产精品一区二区性色av| 国产单亲对白刺激| 亚洲国产欧洲综合997久久,| 两性午夜刺激爽爽歪歪视频在线观看| 男人的好看免费观看在线视频| 最好的美女福利视频网| av在线老鸭窝| 观看美女的网站| 永久网站在线| 国产一级毛片在线| АⅤ资源中文在线天堂| 岛国毛片在线播放| 国产美女午夜福利| 久久久成人免费电影| 人人妻人人看人人澡| 好男人视频免费观看在线| 国产三级中文精品| 国产在视频线在精品| 黄色视频,在线免费观看| 熟女电影av网| 性欧美人与动物交配| 日韩欧美三级三区| 日韩国内少妇激情av| 亚洲最大成人手机在线| av国产免费在线观看| 99精品在免费线老司机午夜| 免费搜索国产男女视频| 亚洲精品日韩在线中文字幕 | 日韩强制内射视频| 美女脱内裤让男人舔精品视频 | 国内少妇人妻偷人精品xxx网站| 精品日产1卡2卡| 国产精品三级大全| 变态另类成人亚洲欧美熟女| 日韩一区二区三区影片| 久久精品国产清高在天天线| 91狼人影院| 在现免费观看毛片| 少妇丰满av| 99久久精品一区二区三区| 18禁在线无遮挡免费观看视频| 亚洲精品久久久久久婷婷小说 | 国产精品人妻久久久久久| 性色avwww在线观看| 国产午夜精品久久久久久一区二区三区| 99热这里只有是精品50| 成人午夜高清在线视频| 成人毛片a级毛片在线播放| 欧洲精品卡2卡3卡4卡5卡区| 嫩草影院精品99| 久久精品91蜜桃| 成人二区视频| 亚洲精品粉嫩美女一区| 夫妻性生交免费视频一级片| 桃色一区二区三区在线观看| 亚洲国产欧洲综合997久久,| 亚洲色图av天堂| 久久久久性生活片| 国产高清视频在线观看网站| 老熟妇乱子伦视频在线观看| 免费在线观看成人毛片| 国产一区亚洲一区在线观看| 欧美不卡视频在线免费观看| 韩国av在线不卡| 日日撸夜夜添| 国产在线男女| 91在线精品国自产拍蜜月| 熟女电影av网| 国产麻豆成人av免费视频| 国产三级在线视频| 国产精品久久电影中文字幕| 91麻豆精品激情在线观看国产| 天堂av国产一区二区熟女人妻| av在线播放精品| av.在线天堂| 国产精品一二三区在线看| 久久99热这里只有精品18| 成人综合一区亚洲| 狂野欧美激情性xxxx在线观看| 少妇高潮的动态图| 久久久色成人| 日产精品乱码卡一卡2卡三| 九九爱精品视频在线观看| 亚洲自偷自拍三级| kizo精华| 夜夜爽天天搞| 人妻夜夜爽99麻豆av| 麻豆精品久久久久久蜜桃| 色噜噜av男人的天堂激情| 国产又黄又爽又无遮挡在线| 少妇的逼好多水| 日日摸夜夜添夜夜添av毛片| 久久6这里有精品| 亚洲第一区二区三区不卡| 国产麻豆成人av免费视频|