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

    Dynamic Loads and Wake Prediction for Large Wind Turbines Based on Free Wake Method

    2015-11-24 02:39:22CaoJiufa曹九發(fā)WangTongguang王同光LongHui龍慧KeShitang柯世堂XuEofeng許波峰
    關(guān)鍵詞:波峰入射角粗糙度

    Cao Jiufa(曹九發(fā))'Wang Tongguang(王同光)' Long Hui(龍慧)'Ke Shitang(柯世堂)'Xu Eofeng(許波峰)

    1.Jiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design' Nanjing University of Aeronautics&Astronautics'Nanjing 210016'P.R.China;

    2.Department of Mechanical Engineering'The University of Sheffield'Sheffield S1 3JD'UK

    Dynamic Loads and Wake Prediction for Large Wind Turbines Based on Free Wake Method

    Cao Jiufa(曹九發(fā))1'2'Wang Tongguang(王同光)1*' Long Hui(龍慧)2'Ke Shitang(柯世堂)1'Xu Eofeng(許波峰)1

    1.Jiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design' Nanjing University of Aeronautics&Astronautics'Nanjing 210016'P.R.China;

    2.Department of Mechanical Engineering'The University of Sheffield'Sheffield S1 3JD'UK

    With large scale wind turbines'the issue of aerodynamic elastic response is even more significant on dynamic behaviour of the system.Unsteady free vortex wake method is proposed to calculate the shape of wake and aerodynamic load.Considering the effect of aerodynamic load'inertial load and gravity load'the decoupling dynamic equations are established by using finite element method in conjunction of the modal method and equations are solved numerically by Newmark approach.Einally'the numerical simulation of a large scale wind turbine is performed through coupling the free vortex wake modelling with structural modelling.The results show that this coupling model can predict the flexible wind turbine dynamic characteristics effectively and efficiently.Under the influence of the gravitational force'the dynamic response of flapwise direction contributes to the dynamic behavior of edgewise direction under the operational condition of steady wind speed.The difference in dynamic response between the flexible and rigid wind turbines manifests when the aerodynamics/structure coupling effect is of significance in both wind turbine design and performance calculation.

    wind turbine;free wake method;aerodynamic dynamics;structural dynamics

    0 Introduction

    Large scale wind turbines operate in complicated conditions'which contain all kinds of unsteady and coupling effects of wind condition and blade structure'thus making accurate prediction of the turbine performance very difficult.Particularly'when unsteady flow field'accompanying frequent changes in both wind speed and wind direction'passes through the long and slender blades as well as through the tall tower'it inevitably results in the so-called aero-elastic problem' i.e.the interaction between aerodynamic forces and rotor motion.Therefore'it is necessary to develop accurate coupled methods between the aerodynamic model with structure model in wind turbine design to enable vibration alleviation and to develop control strategy.

    Among the different aerodynamic theories to model the rotor aerodynamics[1-3]'the vortex theory is considered as one of the most suitable approximations because of its affordable computational costs and reasonably accurate results.In addition'vortex models are made up of a blade aerodynamic model including lifting line'lifting surface or panel method[4]to describe the flow around the blade'as well as to calculate the trailed and shed vorticitis released to the wake' using prescribed[5-7]or free wake[8-11]models to describe the wake geometry.The free wake models are more suitable for general rotor configurations'since the wake is allowed to freely distortunder the influence of the local velocity field.In the present work'the free wake model consists of near vortex sheets and far wake tip vertical filaments.Eor aero-elastic simulations'the free wake model is complemented with an elastic rotor model for the rotor dynamics[12-14].The modal superposition method is used to solve the wind turbine elastic problem in this paper.

    1 Aerodynamic Model

    It is assumed that the flow field is incompressible and potential in the free vortex wake(EVW)model for the wind turbine.The blade is modeled as a series of elements'which are represented as a line of bound vorticity lying along the blade quarter chord line.The vertical filaments'

    extending downstream from the trailing edge of the blade element boundary'are allowed to freely distort under the influence of local velocity field. The governing equation of the vertical filaments can be written in the form of a partial differential equation as

    where the blade azimuth angleψis a temporal coordinate and the wake age angleζis a spatial coordinate.On the right hand side of Eq.(1)'Vindequals to the mean value of the induced velocities at the surrounding four grid points calculated by the Eiot-Savart law.To solve the partial differential equation numerically'the finite difference approximations are used to approximate the derivatives on the left hand side.Eor the spatial(ζ)derivative'a five-point central difference approximation has been used based on the predictor-corrector central(PCC)difference[15]and the predictor-corrector second backward(PC2E)[16].The accuracy of the temporal(ψ)derivative approximation plays a significant role in the time-accurate free vortex method.The PC2E algorithm uses a second-order backward difference approximation'whereas the PCC algorithm still uses a five-point central difference approximation.

    Eig.1 Schematic of discretized tip vortex geometry

    In the present work'a new time-accurate algorithm is developed for overall convergence in the numerical iterations.Eor this'Eq.(1)can be written in another form as

    Eq.(2)can be written in a general form of ordinary derivative equation as

    The predictor process in the predictor-corrector algorithm adopts an explicit format'whereas the corrector process adopts an implicit format. Assuming the steps are equal'the general form of the linear multistep method for Eq.(3)is written as

    where yn+jand fn+jrepresent the terms of y(xn+j)and f(xn+j'yn+j)'respectively.The values of constantsαjandβj(j=0'1…'k)can be obtained using the method of undetermined coefficients.An explicit three-step linear multistep method is given by

    The local truncation error of Eq.(5)is

    An implicit three-step linear multistep method is given by

    The local truncation error of Eq.(7)is

    The explicit and implicit three-step linear multistep methods are used in the temporal(ψ)derivative approximation.A new predictorcorrector algorithm is developed as

    Eor the predictor

    and for the corrector

    It is obvious from the local truncation errors of Eqs.(5'7)that the new predictor-corrector algorithm has third-order accuracy'and therefore this algorithm is referred as the three-step and third-order predictor-corrector(D3PC)algorithm.PCC is a single-step algorithm and has second-order accuracy.The single-step algorithm is simple'but its numerical stability is not good enough.The multistep method has recently been widely used since it has better stability and convergence.Although the PC3E algorithm is a three-step algorithm'it only achieves secondorder accuracy'which results in low efficiency. The D3PC algorithm developed in this paper is also a three-step algorithm'but is of third-order accuracy.

    In this paper'the three-dimensional rotational model is included in the aerodynamic model. The three-dimensional rotational effect is one of the typical differences between rotor and fixed wing'which is characterised by significantly increased lift coefficient compared with the corresponding 2D case.The higher anales of attack can contribute to the delay of the occurrence of flow separation.Du-Selig stall-delay model[17]0 coupled with the free vortex wake model is used to modify the airfoil aerodynamic data by consideration of the three-dimensional rotational effect here.

    The aerodynamic loading is caused by the flow past the wind turbine structure composed of the blades and the tower.The air loads can be calculated through the discrete blade element method.Every blade element is regarded as a 2-D airfoil'and the relative velocity vector Vrelis obtained from

    2 Dynamics Equation for Aerodynamic and Structural CouPling

    If a wind turbine is described as a discretized mechanical system'the principle work is to correctly set up the mass matrix M'stiffness matrix K'and damping matrix C'for the dynamics equation

    where{F(t)}denotes the generalized force vector associated with the external loads.x is the generalized displacement.In this paper'the Newmark method is used to solve the dynamics equation. The discretized calculation process to solve the dynamics equation is presented aswhere the subscript represents the time-step number.The displacements at the time step n+1 are calculated from the results of the previous time step n'at which the displacements'velocities'and accelerations for each node are already known.

    The generalized displacement and velocity in the form of expansion from the time step n to the time step n+1 are given by

    whereγ=0.5 andβ=0.25 are the trapezoidal conditions in order to guarantee unconditional convergence.It is worthy to note that a small enough computing time step is required for gaining more accurate structural dynamic response. Here'the time step is chosen on the basis of the one tenth of the rotating cycle when the structural response is significant enough.The calculation procedure is shown in Eig.2.(EVW:the free vortex wake method'EEM:the finite element modeling)

    Eig.2 Calculation procedures for dynamic responses

    The blade gravitational force and centrifugal force are included in the calculation.Gravity is responsible for a sinusoidal loading of the blades with a frequency corresponding to the rotation of the rotor.The blade experiences tensile stress and compressive stress because of the gravitational loads'thus resulting in the blade vibration deformation and affecting the rotor aerodynamic performance.The gravitational loading is obtained from

    where A is the coordinate system transformation matrix'mithe i th blade element mass'and dfgthe gravitational loading of the i th blade element.

    The inertial loading stems from the centrifugal force on the blades due to rotation.The centrifugal force dfcacting on the blade element at a radius r from the rotational axis is obtained by

    3 Numerical Results and Analyses

    3.1 Calculation model validation

    對(duì)3種不同粗糙度的鋁材料進(jìn)行BRDF對(duì)比測(cè)試,所選粗糙度Ra分別為2.5,5,8 μm。入射角度θi分別為0°和30°,實(shí)驗(yàn)測(cè)量的BRDF,如圖7所示。

    The National Renewable Energy Laboratory(NREL)Phase VI rotor geometry'aerodynamic and structural properties are well documented in Ref.[18].The operating condition for the experiment is varied from wind speed of 7 to 25 m/s. The rotor speed is 72 r/min with the cone angle of tip pitch being 3°and 0°.The rotor radius is 5.029 m.

    Eig.3 shows the variation of rotor torque with wind speed'and comparison with the experiment data.The calculation results agree well with the experimental data at low wind speeds. At higher wind speeds'however'there are discrepancies'probably due to the stall effects.

    Eig.3 Variation of rotor torque with wind speed

    Eig.4 shows the variation of normal and tangential coefficients at different span with wind speed.Ranging from low to high wind speeds'it is found that the normal and tangential force coef-ficients of the blade root are slightly worse than other part of the blade'which is due to the fact that the blade root suffers the serious stall effect. It is difficult to simulate the stall situation.The better stall model is needed to use in the free wake method.However'both the middle and tip part of blade have a good match with the experimental data.Therefore'one can say that the calculation model is well used to calculation the loads of the wind turbine.

    Eig.4 Variation of normal and tangential coefficients at different wind speeds

    3.2 Large-scale wind turbine numerical results and analyses

    The NH1500 wind turbine is adopted as a calculation example in this paper.The calculation time step is 10°/Ω.The wind speed is 8 m/s.To reflect the similarity of NH1500 and real 1.5 MW wind turbine'the rotor main parameters are shown in Table 1.It can be seem that the rated power of NH1500 is 1.5 MW.The NH1500 and real 1.5 MW wind turbine have the similar blade length and wind speed operational condition.And the NH1500 also adopts the variable-pitch variable speed control strategy.Moreover'to describe the NH1500 detailedly'the NH1500 blade main structural parameters are shown in Table 2.All the aerodynamic centers of blade are located in a quarter chord length.Mass and elastic centers are given through chord length percentage in the airfoil coordinates.

    Table 1 NH1500 rotor Parameters

    Here'the finite element method is applied to establish the wind turbine structural model. Two-node beams and shell elements are used to model the blade and tower'respectively.The nacelle mass and moment of inertia are simulated through the 0D element in PATRAN software. The rotor is rigidly connected to the tower.The boundary conditions are imposed at the bottom of the tower through fixed constraints.The modes of NH1500 wind turbine are calculated at the rotational speed of 17.2 r/min with the calculatedmodes shown in Table 3.

    Table 2 Main structural Parameters for the NH1500 blade

    Table 3 Mode descriPtion

    Eig.5 shows the comparison of calculated value and experimental value.NH1500 is a 1/16 scaled model on the experiment.The calculated data and experimental data changing trend are basically consistent.In Eig.5'the calculation values are slightly higher than the experiment data. The probable reason is that the experiment wind turbine model is a scaled model.

    Eig.6 shows the downstream displacements of the blade tip vortex'calculated from the EVW method.To investigate the blade deformation effect'the NH1500 blade is considered rigid and flexible'respectively.It can be clearly seen from Eig.6 that the difference in the vortex position between the assumed rigid blade and the actual flexible blade.The displacement difference increases from 0.05R to 0.1R as the tip vortex fila-ment moves downstream.It can be reasoned that this difference has influence on the air load calculation through the induced velocity calculation' which depends on the spatial distance of the vortex to the blade.

    Eig.5 Cpof NH1500

    Eor further analysis'the wind speed distribution and the vorticity of the wake are calculated.The wind speed of the wake V(Vx'Vy'Vz)is defined as

    To find the wind speed of the wake seen by V(Vx'Vy'Vz)'the induced velocity of the tip vortex'the induced velocity of the bound vortex'the induced velocity of the trailed vortex'and the induced velocity of the shed vortex'must be added as vectors into the inflow velocity'V∞.

    As the wind speed distribution is known'the vorticityis calculated by

    Eigs.7'8 give the development and evolvement of the tip vortex for the rigid and flexible blades with the axial downstream distance'respectively.The vortical strength increases downstream just behind the rotor with a wake expansion'but gradually dissipates as the axial distance further increases for the both cases.However' both the vortical strength and the wake expansion for the assumed rigid rotor are stronger than those for the actual flexible rotor.

    Eig.7 Axial distribution of tip vortex for rigid blade

    Eig.8 Axial distribution of tip vortex for flexible blade

    Eig.9 0R section distribution of wind speed Vz(m/s)of rigid blade

    Eig.10 0R section distribution of wind speed Vz(m/s)of flexible blade

    Eig.9 shows the 0R section distribution of the wind speed of the rigid blade.Eig.10 shows the 0R section distribution of the wind speed of the flexible blade.It is found that the 0R section wind speeds of the flexible blade are weaker than the rigid due to the structural deformation.The effect is significant for wake prediction.

    Eig.11 illustrates some representative numerical results among which the aerodynamic characteristics are presented by the blade root loads for the flexible case at the wind speed of 8 m/s.The root loads all fluctuate sinusoidally with azimuthally angle for both cases with and without gravity considered.Nevertheless'the amplitudes of the load variation differ from each other due to the gravitational effect.In addition' the blade dynamic deformation is also influenced by the gravity.The dynamic displacement in the flapwise direction when gravitational force is not considered is only about 4%of that when the gravitational force is taken into account'as shown in Eig.12.

    Eig.11 Edgewise force at blade root

    Eig.12 Dynamic displacements at blade tip in flapwise direction

    Eig.13 Wind turbine power coefficients'V∞=8 m/s

    It is obvious from Eig.13'where the windturbine power coefficient Cpis compared between the rigid and flexible considerations'that the flexible wind turbine Cpis larger than the rigid one and changes with azimuth angle.According to the comparison of minimum value'the flexible wind turbine Cpis 1.6%more than the rigid one. It is the reason that the flapwise direction induced velocity of flexible blade has a larger change amplitude than the edgewise direction induced velocity.And the flapwise direction induced velocity is weaker than the rigid one.Thus'the attack angle and the flow angle of flexible blade are both larger.Einally'the tangential force coefficient to the rotor is larger'which can lead to the increasing of torque and Cp.

    4 Conclusions

    The numerical analysis of wind turbine wake and loads are presented in this paper.An aerodynamically and structurally coupled model is developed and used to predict the flexible wind turbine loads and aerodynamic performance.The results demonstrate that EVW and EEM coupling model can predict the wind turbine dynamic characteristics effectively and efficiently.The dynamic displacement in the flapwise direction without considering gravitational force is only 4%of that with gravitational force.Moreover'according to the comparison of minimum value of Cp'the Cpvalue of the flexible wind turbine is 1.6%more than that of the rigid one.Despite the fact that the accuracy of the calculated results needs to be improved for further validation by experiment data'the wind turbine wake and dynamic response can be better understood through the results obtained from this study.

    Acknowledgements

    This work was supported by the National Easic Research Program of China(973 Program)(No. 2014CE046200)'the Jiangsu Province Natural Science Eoundation(No.EK2012390)'the Eundamental Research Eunds for the Central Universities'and the Priority Academic Program Development of Jiangsu Higher Education Institutions.

    [1] Ghadirian A'Dehghan M'Torabi E.Considering induction factor using EEM method in wind farm layout optimization[J].Journal of Wind Engineering and Industrial Aerodynamics'2014'129:31-39.

    [2] Zhang Zhenyu'Zhou Hanwei'Wang Tongguang. Numerical analysis of influence of Gurney flaps applied to wind turbines[J].Transactions of Nanjing University of Aeronautics&Astronautics'2014'31(5):576-579.

    [3] Hsu MingChen'Akkerman I'Eazilevs Y.Einite element simulation of wind turbine aerodynamics:Validation study using NREL Phase VI experiment[J]. Wind Energy'2014'17:461-481.

    [4] Grasso E'van Garrel A'Schepers G.Development and validation of generalized lifting line based code for wind turbine aerodynamics[C]∥The 30th ASME wind energy symposium'Elorida'USA:[s.n.]' 2011:146-152.

    [5] Scheurich E'Eletcher T M'Erown R E.Simulating the aerodynamic performance and wake dynamics of a vertical-axis wind turbine[J].Wind Energy'2011'14(2):159-177.

    [6] Ereton S P'Coton E N'Moe G.A study on rotational effects and different stall delay models using a prescribed wake vortex scheme and NREL phase VI experiment data[J].Wind Energy'2008'11:459-482.

    [7] Wang T G.Unsteady aerodynamic modelling of horizontal axis wind turbine performance[D].Scotland' UK:University of Glasgow'1999.

    [8] Sebastian T'Lackner MA.Development of a free vortex wake method code for offshore floating wind turbines[J].Renewable Energy'2012'46:269-275.

    [9] Zhou W P'Tang S L'LüH.Computation on aerodynamic performance of horizontal axis wind turbine based on time-marching free vortex method[J].Chin Soc for Elec Eng'2011'31(29):124-130.

    [10]Jeon M'Lee S.Unsteady aerodynamics of offshore floating wind turbines in platform pitching motion using vortex lattice method[J].Renewable Energy' 2014'65:207-212.

    [11]Gupta S'Leishman J G.Performance pre-dictions of the NREL phaseⅥcombined Ex-periment rotor using a free-vortex wake model[C]∥Collection ofTechnical Papers—44th AIAA Aerospace Sciences Meeting.[S.l.]:AIAA'2006:4544-4564.

    [12]Eelco H'Gustaaf J'Asfaw E.Aero-elastic behavior of a flexible blade for wind turbine application:A 2D computational study[J].Energy'2010'35:778-785.

    [13]Ahlstrom A.Aeroelastic simulation of wind turbine dynamics[D].Sweden:Royal Institute of Technology Department of Mechanics'2005.

    [14]Ng E E'Hesse H'Palacios R.Aeroservoelastic state space vortex lattice modeling and load alleviation of wind turbine blades[J].Wind Energy'2014' 17(4):DOI:10.1002/we.1752.

    [15]Ehagwat M J'Leishman J G.Rotor aerodynamics during maneuvering flight usinga time-accurate freevortex wake[J].Journal of the American Helicopter Society'2003'48:143-158.

    [16]Ehagwat M J'Leishman J G.Time-accurate free vortex wake model for dynamic rotor response[C]∥American Helicopter Society Specialist Meeting.Atlanta'USA:[s.n.]'2000.

    [17]Du Z'Selig M S.A 3-D stall-delay model for horizontal axis wind turbine performance prediction[R]. AIAA-98-0021'1998.

    [18]Hand M'Simms D A.Unsteady aerodynamicexperiment phase VI:Wind tunnel test configurations and available data campaigns[R].National Renewable Energy Laboratory.Colorado:National Technical Information Service'2001.

    (Executive editor:Xu Chengting)

    O35 Document code:A Article ID:1005-1120(2015)02-0240-10

    *CorresPonding author:Wang Tongguang'Professor'E-mail:tgwang@nuaa.edu.cn.

    How to cite this article:Cao Jiufa'Wang Tongguang'Long Hui'et al.Dynamic loads and wake prediction for large wind turbines based on free wake method[J].Trans.Nanjing U.Aero.Astro.'2015'32(2):240-249.

    http://dx.doi.org/10.16356/j.1005-1120.2015.02.240

    (Received 13 November 2014;revised 7 January 2015;accepted 12 January 2015)

    猜你喜歡
    波峰入射角粗糙度
    一般三棱鏡偏向角與入射角的關(guān)系
    炮制工程騙錢的“甲方”
    作用于直立堤墻與樁柱的波峰高度分析計(jì)算
    基于無人機(jī)影像的巖體結(jié)構(gòu)面粗糙度獲取
    甘肅科技(2020年20期)2020-04-13 00:30:18
    冷沖模磨削表面粗糙度的加工試驗(yàn)與應(yīng)用
    模具制造(2019年4期)2019-06-24 03:36:48
    預(yù)制圓柱形鎢破片斜穿甲鋼靶的破孔能力分析*
    用經(jīng)典定理證明各向異性巖石界面異常入射角的存在
    基于BP神經(jīng)網(wǎng)絡(luò)的面齒輪齒面粗糙度研究
    鋼材銹蝕率與表面三維粗糙度參數(shù)的關(guān)系
    兒童標(biāo)準(zhǔn)12導(dǎo)聯(lián)T波峰末間期的分析
    特大巨黑吊av在线直播| 亚洲av.av天堂| 91国产中文字幕| 纵有疾风起免费观看全集完整版| 一区二区三区精品91| 中文精品一卡2卡3卡4更新| 午夜福利视频精品| 少妇 在线观看| 天天影视国产精品| 在线观看免费视频网站a站| 夜夜爽夜夜爽视频| 欧美日韩亚洲高清精品| 国产在线一区二区三区精| 欧美成人午夜免费资源| 丝袜喷水一区| 97超碰精品成人国产| 国产一区二区三区综合在线观看 | 日韩,欧美,国产一区二区三区| 五月开心婷婷网| 国产男女超爽视频在线观看| 男的添女的下面高潮视频| 久久久久精品久久久久真实原创| 久久久久视频综合| 人人澡人人妻人| 大香蕉97超碰在线| 欧美97在线视频| 国产精品.久久久| 男人操女人黄网站| 91在线精品国自产拍蜜月| 99热这里只有是精品在线观看| 久久久精品免费免费高清| 日本黄大片高清| 久久久久视频综合| 国产精品一区www在线观看| 国产日韩欧美在线精品| 两个人的视频大全免费| 午夜久久久在线观看| av不卡在线播放| 亚洲精品中文字幕在线视频| 亚洲第一av免费看| 国产午夜精品一二区理论片| 久久久精品94久久精品| 亚洲成人手机| 韩国高清视频一区二区三区| 亚洲精品久久成人aⅴ小说 | 人成视频在线观看免费观看| 少妇人妻 视频| 伦理电影免费视频| 高清不卡的av网站| 免费黄色在线免费观看| 97在线人人人人妻| 日韩成人伦理影院| 在线看a的网站| 亚洲精品色激情综合| 欧美日韩国产mv在线观看视频| 啦啦啦啦在线视频资源| 欧美变态另类bdsm刘玥| 亚洲欧美精品自产自拍| 69精品国产乱码久久久| 只有这里有精品99| 亚洲精品色激情综合| 午夜激情久久久久久久| 国产成人av激情在线播放 | 亚洲av不卡在线观看| 涩涩av久久男人的天堂| 涩涩av久久男人的天堂| 亚洲久久久国产精品| 肉色欧美久久久久久久蜜桃| 三级国产精品欧美在线观看| 日产精品乱码卡一卡2卡三| 亚洲国产av影院在线观看| 国产亚洲最大av| 美女内射精品一级片tv| 黄色视频在线播放观看不卡| 免费黄网站久久成人精品| 在线观看免费视频网站a站| 伦理电影大哥的女人| 国产永久视频网站| 成人国产麻豆网| 欧美国产精品一级二级三级| 国产成人a∨麻豆精品| 亚洲五月色婷婷综合| 老女人水多毛片| 女人久久www免费人成看片| 日本午夜av视频| 国产黄片视频在线免费观看| 国产日韩欧美亚洲二区| 欧美精品人与动牲交sv欧美| 观看美女的网站| 欧美亚洲日本最大视频资源| 在线天堂最新版资源| 五月玫瑰六月丁香| 一边摸一边做爽爽视频免费| 一个人看视频在线观看www免费| 国产精品一国产av| 欧美日韩国产mv在线观看视频| 精品亚洲成a人片在线观看| 九九在线视频观看精品| 亚洲欧洲日产国产| 久久99热6这里只有精品| 亚洲av成人精品一区久久| 在线精品无人区一区二区三| 伊人亚洲综合成人网| 日日撸夜夜添| 亚洲国产精品一区三区| 国模一区二区三区四区视频| 简卡轻食公司| 伊人久久国产一区二区| 考比视频在线观看| 欧美日韩综合久久久久久| 久久99蜜桃精品久久| 亚州av有码| 日日摸夜夜添夜夜爱| 久久精品夜色国产| 涩涩av久久男人的天堂| 伦理电影大哥的女人| 91国产中文字幕| 美女国产视频在线观看| 久久综合国产亚洲精品| 欧美 亚洲 国产 日韩一| 特大巨黑吊av在线直播| 热99久久久久精品小说推荐| 亚洲久久久国产精品| 国产免费一级a男人的天堂| 搡女人真爽免费视频火全软件| a级毛片免费高清观看在线播放| 高清欧美精品videossex| 精品99又大又爽又粗少妇毛片| 亚洲人成网站在线观看播放| 熟女电影av网| 日韩精品有码人妻一区| av天堂久久9| 色婷婷久久久亚洲欧美| 欧美成人午夜免费资源| 中国美白少妇内射xxxbb| 欧美xxⅹ黑人| 色吧在线观看| 国模一区二区三区四区视频| 午夜精品国产一区二区电影| 97超碰精品成人国产| 日本免费在线观看一区| 亚洲精品自拍成人| 各种免费的搞黄视频| 国产成人a∨麻豆精品| 一边摸一边做爽爽视频免费| 99久久精品一区二区三区| 久久99一区二区三区| 菩萨蛮人人尽说江南好唐韦庄| 视频中文字幕在线观看| 97精品久久久久久久久久精品| 国产色爽女视频免费观看| 五月开心婷婷网| 欧美日韩视频精品一区| 国产高清国产精品国产三级| 下体分泌物呈黄色| 美女中出高潮动态图| 高清午夜精品一区二区三区| 菩萨蛮人人尽说江南好唐韦庄| 久久99热这里只频精品6学生| 嘟嘟电影网在线观看| 老司机影院毛片| 亚洲国产欧美在线一区| 纯流量卡能插随身wifi吗| 美女脱内裤让男人舔精品视频| 日本猛色少妇xxxxx猛交久久| 一本一本综合久久| 成年av动漫网址| 国产探花极品一区二区| 国产成人精品无人区| 黄片无遮挡物在线观看| 有码 亚洲区| 美女cb高潮喷水在线观看| 国产深夜福利视频在线观看| 亚洲国产日韩一区二区| 国产精品无大码| 国产日韩一区二区三区精品不卡 | 久久久久久人妻| 中国美白少妇内射xxxbb| 制服丝袜香蕉在线| 一级毛片 在线播放| 国产成人a∨麻豆精品| 自线自在国产av| 久久精品国产亚洲av涩爱| 国产精品偷伦视频观看了| 日韩强制内射视频| 国产精品99久久99久久久不卡 | 欧美三级亚洲精品| 九九久久精品国产亚洲av麻豆| 卡戴珊不雅视频在线播放| 亚洲人成网站在线播| 日本黄大片高清| 成人漫画全彩无遮挡| 亚洲精品乱码久久久v下载方式| 满18在线观看网站| 亚洲,一卡二卡三卡| 狠狠精品人妻久久久久久综合| 久久久久国产网址| 免费高清在线观看视频在线观看| av又黄又爽大尺度在线免费看| 丰满少妇做爰视频| 中文字幕免费在线视频6| 人人妻人人爽人人添夜夜欢视频| 日产精品乱码卡一卡2卡三| av天堂久久9| 久久精品久久久久久久性| 日韩成人伦理影院| 欧美亚洲 丝袜 人妻 在线| 天堂8中文在线网| 曰老女人黄片| 精品一区二区三卡| 国产探花极品一区二区| 妹子高潮喷水视频| 欧美老熟妇乱子伦牲交| 日韩av不卡免费在线播放| 一本—道久久a久久精品蜜桃钙片| 精品少妇内射三级| 成人综合一区亚洲| 国产有黄有色有爽视频| 我的女老师完整版在线观看| 欧美 亚洲 国产 日韩一| 久久国产精品男人的天堂亚洲 | 日韩av不卡免费在线播放| 美女xxoo啪啪120秒动态图| 97精品久久久久久久久久精品| 久久久久久人妻| 最近中文字幕2019免费版| 多毛熟女@视频| 青春草国产在线视频| 五月开心婷婷网| 亚州av有码| 国产成人aa在线观看| a级片在线免费高清观看视频| 观看av在线不卡| 免费观看的影片在线观看| 亚洲第一av免费看| 精品人妻偷拍中文字幕| 久久影院123| 狠狠精品人妻久久久久久综合| 你懂的网址亚洲精品在线观看| 欧美国产精品一级二级三级| 国产乱来视频区| 久久久久国产网址| 在线观看三级黄色| 久久99热6这里只有精品| 熟女av电影| 日本午夜av视频| 美女国产视频在线观看| 热99国产精品久久久久久7| 热99久久久久精品小说推荐| 日本与韩国留学比较| 国产免费视频播放在线视频| 内地一区二区视频在线| 国产有黄有色有爽视频| 青春草亚洲视频在线观看| 国产一级毛片在线| 一级毛片黄色毛片免费观看视频| 丰满乱子伦码专区| 飞空精品影院首页| 欧美三级亚洲精品| 制服人妻中文乱码| 久久精品国产亚洲网站| 精品酒店卫生间| 视频区图区小说| av国产久精品久网站免费入址| 亚洲国产毛片av蜜桃av| 免费观看av网站的网址| 亚洲精品国产色婷婷电影| 91精品国产九色| 国产伦理片在线播放av一区| 国产一区二区三区av在线| 美女福利国产在线| 欧美精品高潮呻吟av久久| 久久热精品热| 夫妻午夜视频| 麻豆精品久久久久久蜜桃| 91精品一卡2卡3卡4卡| 午夜av观看不卡| av在线播放精品| 亚洲av在线观看美女高潮| 大话2 男鬼变身卡| 多毛熟女@视频| 97精品久久久久久久久久精品| 免费看不卡的av| 91aial.com中文字幕在线观看| 日本爱情动作片www.在线观看| 久久久久久人妻| 满18在线观看网站| 亚洲精品国产色婷婷电影| 狠狠精品人妻久久久久久综合| 国产成人精品无人区| av一本久久久久| 亚洲国产成人一精品久久久| 高清午夜精品一区二区三区| 2018国产大陆天天弄谢| 亚洲精品第二区| 成年人免费黄色播放视频| 一区二区av电影网| 男人爽女人下面视频在线观看| 18+在线观看网站| 免费播放大片免费观看视频在线观看| 2021少妇久久久久久久久久久| 国产精品人妻久久久久久| 少妇人妻 视频| 国产精品成人在线| 黄色毛片三级朝国网站| 亚洲丝袜综合中文字幕| 亚洲av成人精品一二三区| 成年人午夜在线观看视频| 在线观看免费视频网站a站| 日本与韩国留学比较| 性色avwww在线观看| 色5月婷婷丁香| 精品久久久久久电影网| 黄色视频在线播放观看不卡| 99热这里只有精品一区| 久热这里只有精品99| 美女主播在线视频| 午夜老司机福利剧场| 涩涩av久久男人的天堂| 亚洲激情五月婷婷啪啪| 久久99热这里只频精品6学生| 一个人免费看片子| 视频在线观看一区二区三区| 精品国产国语对白av| 国产精品不卡视频一区二区| 日韩中文字幕视频在线看片| 91精品国产九色| 最近中文字幕2019免费版| 一边亲一边摸免费视频| 韩国高清视频一区二区三区| av有码第一页| 亚洲欧洲日产国产| 亚洲成人一二三区av| 高清午夜精品一区二区三区| 欧美另类一区| 香蕉精品网在线| 欧美成人精品欧美一级黄| 亚洲精品一二三| 成人国产麻豆网| 国产 精品1| 看免费成人av毛片| 亚洲欧洲精品一区二区精品久久久 | 男人操女人黄网站| 老司机影院成人| 国产在线免费精品| 欧美精品国产亚洲| 美女视频免费永久观看网站| 色94色欧美一区二区| 又粗又硬又长又爽又黄的视频| 天堂俺去俺来也www色官网| 国产高清不卡午夜福利| 黑人猛操日本美女一级片| 久久久久久久国产电影| a级毛片在线看网站| 精品少妇内射三级| 欧美老熟妇乱子伦牲交| 国产欧美另类精品又又久久亚洲欧美| 亚洲精品视频女| 色视频在线一区二区三区| 秋霞伦理黄片| 两个人的视频大全免费| 美女福利国产在线| videos熟女内射| 欧美日韩在线观看h| 一本大道久久a久久精品| 久久久久久久久久成人| 黑人猛操日本美女一级片| 91久久精品国产一区二区成人| 国产精品人妻久久久久久| 国产精品人妻久久久影院| 久久久久久久国产电影| 久久久午夜欧美精品| 亚洲av免费高清在线观看| 少妇高潮的动态图| 国产成人精品婷婷| 久久精品国产鲁丝片午夜精品| 国产精品无大码| 欧美日韩av久久| 99视频精品全部免费 在线| 久久ye,这里只有精品| 高清在线视频一区二区三区| 婷婷色麻豆天堂久久| 日日摸夜夜添夜夜爱| 人妻 亚洲 视频| 少妇人妻久久综合中文| 亚洲精品成人av观看孕妇| 国产男女超爽视频在线观看| 中国国产av一级| av在线老鸭窝| 3wmmmm亚洲av在线观看| 水蜜桃什么品种好| 99国产综合亚洲精品| 黑人欧美特级aaaaaa片| 美女中出高潮动态图| 亚洲美女视频黄频| 在线免费观看不下载黄p国产| 插逼视频在线观看| 日韩成人av中文字幕在线观看| 久久久a久久爽久久v久久| 中文欧美无线码| 日本黄色片子视频| 国产欧美日韩综合在线一区二区| 十八禁高潮呻吟视频| 一本—道久久a久久精品蜜桃钙片| 久久久久久久久久久丰满| 午夜日本视频在线| 欧美 亚洲 国产 日韩一| 涩涩av久久男人的天堂| 免费观看在线日韩| 国产成人a∨麻豆精品| 精品人妻熟女毛片av久久网站| 人人妻人人澡人人看| 午夜激情av网站| 精品国产一区二区三区久久久樱花| 国产精品欧美亚洲77777| 久久久国产欧美日韩av| 亚洲精品自拍成人| 91久久精品国产一区二区三区| 香蕉精品网在线| 婷婷色麻豆天堂久久| 2021少妇久久久久久久久久久| 欧美激情极品国产一区二区三区 | 欧美精品一区二区免费开放| 国产精品嫩草影院av在线观看| 成人亚洲欧美一区二区av| 国产片内射在线| 建设人人有责人人尽责人人享有的| 亚洲精品亚洲一区二区| 黄色视频在线播放观看不卡| 插逼视频在线观看| 大片电影免费在线观看免费| 亚洲图色成人| 精品久久久久久久久亚洲| 久久99热6这里只有精品| 午夜福利在线观看免费完整高清在| 卡戴珊不雅视频在线播放| 亚洲国产日韩一区二区| 黄片播放在线免费| 赤兔流量卡办理| 777米奇影视久久| 涩涩av久久男人的天堂| 免费观看a级毛片全部| 国产 精品1| 女性生殖器流出的白浆| 18禁在线播放成人免费| 赤兔流量卡办理| 一级a做视频免费观看| 国产高清三级在线| 精品国产乱码久久久久久小说| 国产亚洲一区二区精品| 日韩av在线免费看完整版不卡| 中文字幕精品免费在线观看视频 | 国产一区亚洲一区在线观看| 久久毛片免费看一区二区三区| 日韩成人av中文字幕在线观看| 国产免费现黄频在线看| 国产精品嫩草影院av在线观看| 2018国产大陆天天弄谢| 精品人妻在线不人妻| 久久人人爽av亚洲精品天堂| 性高湖久久久久久久久免费观看| 99久国产av精品国产电影| 91午夜精品亚洲一区二区三区| 韩国高清视频一区二区三区| 汤姆久久久久久久影院中文字幕| 曰老女人黄片| 一级黄片播放器| 亚洲色图综合在线观看| 天天影视国产精品| 午夜免费观看性视频| 欧美 日韩 精品 国产| 国产69精品久久久久777片| 老熟女久久久| av线在线观看网站| 人人妻人人澡人人爽人人夜夜| 亚洲精品日韩在线中文字幕| 一二三四中文在线观看免费高清| 成人午夜精彩视频在线观看| 高清午夜精品一区二区三区| 中国三级夫妇交换| 卡戴珊不雅视频在线播放| av在线老鸭窝| 久久久亚洲精品成人影院| 啦啦啦视频在线资源免费观看| 性色av一级| 欧美最新免费一区二区三区| 国产国拍精品亚洲av在线观看| 建设人人有责人人尽责人人享有的| 国产探花极品一区二区| 97超视频在线观看视频| 免费日韩欧美在线观看| 国产精品人妻久久久久久| 最新中文字幕久久久久| a 毛片基地| 国产成人精品在线电影| 亚洲综合色网址| 久久久久久人妻| 汤姆久久久久久久影院中文字幕| 国产免费一区二区三区四区乱码| 色5月婷婷丁香| 日韩制服骚丝袜av| 午夜激情福利司机影院| 日本黄色日本黄色录像| av又黄又爽大尺度在线免费看| 免费黄色在线免费观看| 国产av精品麻豆| 国产黄片视频在线免费观看| 久久精品熟女亚洲av麻豆精品| 亚洲欧美成人综合另类久久久| 免费黄色在线免费观看| 亚洲成人av在线免费| 伦理电影免费视频| tube8黄色片| 亚洲欧美成人综合另类久久久| 欧美亚洲日本最大视频资源| 精品人妻偷拍中文字幕| 午夜91福利影院| 精品人妻熟女毛片av久久网站| 日韩av免费高清视频| 成年av动漫网址| a级毛片在线看网站| 啦啦啦中文免费视频观看日本| 免费黄频网站在线观看国产| 69精品国产乱码久久久| 久久精品夜色国产| 制服丝袜香蕉在线| 中文字幕精品免费在线观看视频 | 久久久久久久久久成人| 女人精品久久久久毛片| 国产深夜福利视频在线观看| 日本午夜av视频| 狂野欧美白嫩少妇大欣赏| 日本免费在线观看一区| 免费观看无遮挡的男女| 精品酒店卫生间| 国产免费视频播放在线视频| 欧美老熟妇乱子伦牲交| 男女高潮啪啪啪动态图| 精品久久久噜噜| 五月玫瑰六月丁香| 欧美精品高潮呻吟av久久| 亚洲综合色惰| 免费黄网站久久成人精品| 这个男人来自地球电影免费观看 | 日本猛色少妇xxxxx猛交久久| 自线自在国产av| 国产一区二区三区av在线| 99热全是精品| 99精国产麻豆久久婷婷| 免费大片18禁| 日韩一区二区视频免费看| 亚洲国产欧美在线一区| 亚洲少妇的诱惑av| 欧美bdsm另类| 国产69精品久久久久777片| 亚洲精品亚洲一区二区| 亚洲精品日本国产第一区| 99视频精品全部免费 在线| 日日爽夜夜爽网站| 欧美人与善性xxx| 久久人人爽av亚洲精品天堂| 亚洲欧美成人精品一区二区| 久久久国产一区二区| 国产免费又黄又爽又色| 久久久久久伊人网av| 欧美 日韩 精品 国产| 看非洲黑人一级黄片| 香蕉精品网在线| 久久精品国产自在天天线| 啦啦啦啦在线视频资源| 中文天堂在线官网| 美女国产高潮福利片在线看| 精品国产露脸久久av麻豆| 啦啦啦视频在线资源免费观看| 精品少妇黑人巨大在线播放| 男女国产视频网站| 9色porny在线观看| 18+在线观看网站| 激情五月婷婷亚洲| 日韩电影二区| 国产一区亚洲一区在线观看| 亚洲国产精品国产精品| 精品国产一区二区久久| 久久久久久久大尺度免费视频| 大香蕉久久成人网| 亚洲五月色婷婷综合| 22中文网久久字幕| 水蜜桃什么品种好| 人成视频在线观看免费观看| 国产亚洲精品久久久com| 插逼视频在线观看| 人妻人人澡人人爽人人| 熟女电影av网| 精品久久久精品久久久| .国产精品久久| 亚洲成人手机| 老司机影院毛片| 大片免费播放器 马上看| 伦理电影大哥的女人| 亚洲,一卡二卡三卡| 成年人免费黄色播放视频| xxxhd国产人妻xxx| 亚洲在久久综合| 啦啦啦中文免费视频观看日本| 国产在线视频一区二区| 高清不卡的av网站| 国产精品一区二区三区四区免费观看| 日日爽夜夜爽网站| 一级片'在线观看视频| 欧美97在线视频| 国国产精品蜜臀av免费| 18+在线观看网站| 久久久久久久精品精品| 岛国毛片在线播放| 嫩草影院入口|