999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

風力機啟停機過程中的振動特征分析

2019-05-29 11:18:54趙艷潘嘉寧王振宇林葵庚蔣建群
湖南大學學報·自然科學版 2019年3期
關鍵詞:風速振動區域

趙艷 潘嘉寧 王振宇 林葵庚 蔣建群

摘? ?要:風力機運行中由于受到環境荷載和葉輪轉動的激勵作用,塔架產生振動,影響風力機壽命.為了分析風力機在啟停機過程中的振動規律和動態特性,對1.5 MW風力機塔頂的振動進行了長期監測,并采用基于數據驅動的隨機子空間法識別了結構一階自振頻率和阻尼比.對監測結果進行統計分析,發現該風力機在并網轉速附近運行的時段較多,當風力機啟動經過并網轉速時,有明顯的共振現象,而風力機停機經過并網轉速時,共振現象則不明顯.模態識別的結果表明,一階自振頻率隨運行工況不同而微幅變化,風力機啟動經過并網轉速時的阻尼比較小,風力機停機經過并網轉速時的阻尼比相對較大.此外,基于Sommerfeld效應對風力機啟停機過程中的不同振動現象進行了解釋.本文成果對于掌握同類型風力機的振動特性和動力參數優化具有借鑒意義.

關鍵詞:風力機;啟停機過程;振動監測;共振;Sommerfeld效應

中圖分類號:TK83;TV314? ? ? ? ? ? ? ? ? ? ?文獻標志碼:A

Abstract:Wind turbines vibrate under the excitation of environmental loads and rotor movement in normal operating conditions, and vibrations of the tower affect the service life of such system. In order to analyze the vibration and dynamic characteristics of the wind turbine during startup and shutdown process, vibration of a 1.5 MW turbine tower was monitored for a long term. First-order natural frequencies and damping ratios of the wind turbine were identified utilizing a data-driven stochastic subspace identification method. Monitoring data revealed that the wind turbine usually operated near the rotation speed of grid connection. In addition, obvious resonance at the turbine startup was noted as it surpassed the rotation speed of grid connection, but was not apparent at the turbine shutdown. Modal identification results revealed that the first-order natural frequencies varied slightly under different operating conditions, while the damping ratios were found to be smaller when the wind turbine started up, and its rotation speed surpassed that of grid connection, which were found to be greater at turbine shutdown. In addition, the characteristics of different vibration phenomena in the turbine during startup and shutdown processes were explained with the Sommerfeld effect. Our findings contribute to the field of vibration analysis and allow for more efficient optimization of dynamic parameters in wind turbines of similar design.

Key words:wind turbine;startup and shutdown process;vibration monitoring;resonance;Sommerfeld effect

近年來,風電產業迅猛發展,風力機的額定功率、機艙重量和塔架高度不斷增大[1].當風力機塔架的自振頻率與環境荷載或葉輪旋轉的諧波頻率相重合時,會產生共振,導致振動幅值增大,影響風力機的正常運行,并產生較大的疲勞損傷.為避免出現共振,通常設計塔架的基頻避開環境荷載頻率及葉輪的旋轉頻率(1f)、過槳頻率(3f)[2].

由于振幅與阻尼比成反比,正確地估計風力機結構的阻尼比對于掌握風力機的振動特性非常重要[3].陸上風力機的阻尼比可近似看成是結構阻尼和氣動阻尼的線性疊加[4].結構阻尼取決于材料類型和結構連接構造,鋼結構的阻尼比一般為0.5%到1.5%[5].氣動阻尼是風力機運行過程中,在空氣流動和葉輪振蕩的相互作用中產生的,它受風速、葉片轉速、幾何條件和來流參數等影響,對于不同的風力機類型和運行狀態會有所不同,目前對氣動阻尼的取值還缺乏深入的認識[6].

開展實測資料分析是計算風力機結構自振頻率和阻尼比的重要研究手段.Shirzadeh等[5]采用運行模態分析法識別了海上風力機在停機狀態下的一階自振頻率和阻尼比,并與數值模擬結果進行比較.Damgaard等[7]評價了不同土壤條件下海上單樁風力機超速停機狀態下的一階自振頻率和阻尼比.Versteijlen等[8]對位于Burbo Banks風場的海上風力機進行了12次“超速停機”試驗,研究了地基阻尼對振動的影響.Hackell等[9]基于數據驅動的隨機子空間法識別了5 MW海上風力機的自振頻率、阻尼比和振型,對結構在不同運行狀態下的動態特性進行了評價.上述研究多是基于實測資料對結構的模態參數的總體規律進行分析,并驗證和改進數值模擬方法,沒有深入分析風力機在啟停機過程中的共振機制和振動響應特點.

本文對1.5 MW風力機塔頂的振動響應進行了長期監測,采用基于數據驅動的隨機子空間方法,計算了風力機結構的一階自振頻率和阻尼比.根據風力機啟停機過程中的振動幅值及阻尼比變化,詳細分析和揭示了共振特點,并采用Sommerfeld效應對這種共振機制進行了合理的解釋.

1? ?運行模態分析方法

大型結構體系在環境激勵下的模態參數識別,是一個比較熱門的研究課題[10-13].許多分析方法應運而生,時域方法有:隨機減量法[14]、自然激勵技術(Natural Excitation Technique,NExT)[15]和隨機子空間法(Stochastic subspace identification,SSI)[16];頻域方法有:峰值拾取(PeakPicking,PP)[17]和頻域分解法(Frequency Domain Decomposition,FDD)[18].SSI方法由于識別精度高,無需像傳統時域方法那樣需要前處理得到自由衰減曲線,因此在土木工程界得到越來越多的應用.它主要包括基于協方差驅動的SSI方法和基于數據驅動的SSI方法兩種.后者理論推導清晰,計算量相對較小,因此在大型復雜結構中應用廣泛.本文采用基于數據驅動的SSI方法,從環境激勵下的振動響應中提取模態參數.

風力機結構的環境激勵(輸入激勵)是難以測量的隨機激勵,在對風力機結構的振動響應進行現場監測時,采集到的數據都是在離散的時間點上的,且必然存在著噪聲干擾,因此采用線性的離散狀態方程模型:

式中:xk為第kΔt時刻系統的狀態向量,xk = x(kΔt)(其中Δt為采樣時間間隔);yk為輸出向量;A為空間矩陣;C為輸出矩陣;wk、vk分別為過程噪聲、測量噪聲,這兩種噪聲均為不可測量的噪聲.

假設wk和vk為零均值的白噪聲且互不相關,其協方差矩陣為:

式中:E表示數學期望;wp、vp均為協方差矩陣;Q、R、S均為噪聲序列;δpq為Kroneckerdelta函數;p、q為任意的兩個時間點.

基于數據驅動的SSI方法的具體運算流程如圖1所示.首先,將風力機結構的振動響應(輸出數據)組成Hankel矩陣,對Hankel矩陣進行正交三角分解(Quadrature Rectangle Decomposition,QRD),得到輸出響應向輸入激勵的投影矩陣,然后對投影矩陣進行奇異值分解(Singular Value Decomposition,SVD),得到卡爾曼濾波狀態序列,進而得到風力機結構的狀態向量和輸出向量,應用最小二乘原理估計出結構的系統矩陣,得到風力機結構的系統矩陣及模態參數.

圖1中,Y0.2i-1的第1個下標表示Hankel矩陣左上角元素的時間序數,第2個下標表示左下角元素的時間序數;yi為時刻i測點的輸出向量;j為離散的采集點數,一般假定j→+∞,但是在實際測量中j不可能是無限大的,所以采集時間應盡量長一些;Yp、Yf分別表示“過去”、“將來”的輸出;Qi為投影矩陣;U、V為正交矩陣;S為奇異值組成的對角矩陣;?祝i、Xi分別為可觀矩陣、卡爾曼濾波狀態序列.

2? ?工程監測概況

本文監測的1.5 MW風力機位于某沿海風電場.風力機的切入和切出風速分別為3 m/s和22 m/s(10 min平均風速),額定風速為11 m/s.當10 min平均風速超過22 m/s時,葉輪將停止工作,以避免風荷載過大引起的潛在的損壞.風力機葉輪的并網轉速為9 r/min,額定轉速為17.3 r/min,轉速范圍為0~17.3 r/min±10%.在正常運行狀態下,機艙繞塔架軸線隨著風向自動旋轉,以對準風向.

風電場的數據采集與監控系統(SCADA)遠程連接風力機與主控機房,記錄了監測期間的風速、風向、葉輪轉速及葉輪槳距角的變化,采樣頻率為1/7 Hz.在塔筒高66 m處、東偏南80°方位上安裝雙軸速度計采集風力機的振動響應,采樣頻率為16 Hz,如圖2所示.

2.1? ?運行狀態分析

圖3為監測期間風速的出現概率,平均風速為5.46 m/s,風速的變化范圍為0~23.52 m/s.圖4為風速與風向之間的關系,當風速較小時,在各個方向上出現的頻次差別較小,較大的風速主要出現在方位角-10°~20°和200°~240°附近.

圖5為監測期間不同葉輪轉速的出現概率.在9.5~10.5 r/min附近的出現概率很高,約為30.42%,額定轉速17.3 r/min的概率為22.23%.

2.2? ?允許的頻帶范圍

風力機系統中,當葉片旋轉通過塔架時,會引起共振現象,其激勵頻率主要為葉輪的旋轉頻率f和過槳頻率.對于一個有3個葉片的風力機,過槳頻率為3f及其倍數:

式中:f為葉輪轉速對應的旋轉頻率.如果葉片的過槳頻率f3n接近塔架基頻時,風力機結構的振幅將顯著增大,產生共振[4].DNV規范規定,風力機結構設計時,風力機基頻應該在1f、3f頻率帶的基礎上,預留出±10%的安全度,以避免發生共振[2]:

式中:f3n為轉子正常工作時對應的旋轉頻率1f和過槳頻率3f;f1,n為結構的第一階自振頻率.本文風力機的工作轉速范圍為9 r/min(并網轉速)~17.3 r/min(額定轉速),因此,葉輪轉動通過塔架產生的旋轉頻率1f為0.15~0.288 Hz,過槳頻率3f為0.45~0.865 Hz.而預留10%安全度的1f、3f頻率帶分別為0.135~0.316 8 Hz、0.405~0.951 5 Hz,如圖6所示.

3? ?啟停機過程中的振動特征

基于SSI算法,識別出監測的風力機的基頻為0.425 Hz,雖然避開了3f頻率帶,但沒有預留10%的安全度,與并網轉速(9 r/min)工況下的過槳頻率3f很接近.圖7為并網轉速附近,風力機葉輪轉速和風速之間的關系,可見在2~5 m/s風速區間上,隨著風速的增大,葉輪轉速并沒有隨之增大,而是近乎被“抑制”在10 r/min附近.該風電場風速在2~5 m/s區間的比例很高,達到36.16%,風力機在啟停機過程中都不可避免地要通過并網轉速,因此,非常有必要掌握此時風力機的振動特性.

3.1? ?啟動工況

圖8為風力機啟動過程中的運行狀態和振動速度時程圖.圖中1區域對應的是停機狀態,機艙位置處風速v低于切入風速(3 m/s),葉輪轉速Ω低于并網轉速(9 r/min),此時風力機只是空轉運行.2區域為并網轉速狀態,風速v在切入風速(3 m/s)附近,此時若沒有更大的風能驅動葉輪加速轉動,葉輪轉速Ω長時間處于9~10 r/min.3區域中,葉輪轉速Ω隨風速v增大而增大.可以觀察到,2區域(并網轉速狀態)的風速v和轉速Ω雖然都小于3區域(正常運行狀態),但2區域的振動速度V(最大值為0.175 m/s)卻大于3區域(最大值為0.1 m/s).這是因為,2區域在并網轉速(9 r/min)附近,對應的葉輪過槳頻率3f為0.45 Hz,風力機塔架的一階自振頻率為0.425 Hz,兩者很接近,容易發生共振.在3區域,平均風速大于6 m/s,葉輪轉速大于10 r/min,此時的過槳頻率3f大于塔架自振頻率,不會發生共振.

從圖8(c)中可以看到,葉輪轉速小于2 r/min時,風力機為空轉運行,此時對應的振動速度很小;風力機啟動過程中,轉速從2 r/min快速增大8 r/min,達到并網轉速附近,對應于1區域到2區域的過渡過程;在8 r/min增大到10 r/min過程中,振動速度快速增大;在轉速大于10 r/min后,振動速度有所降低.

3.2? ?停機工況

圖9為風力機停機過程中的運行狀態和振動速度時程圖.圖中4區域為正常運行狀態,葉輪轉速Ω隨風速v減小而減小.隨著風速v進一步減小,風力機進入并網轉速狀態,并準備脫網(5區域).當風速v小于切入風速后,風力機停機(6區域).

可以觀察到,5區域與圖8中2區域啟動工況的轉速差不多,但是比2區域的振動速度V要小,4、5區域中的振動速度相差不多.從圖9(c)中也可以看到,隨著轉速從額定轉速減小,與脫網經過并網轉速時的振動速度峰值相差不大;當風力機為停機狀態時,振動速度很小.

可見,風力機啟動和停機過程中的振動特點是不同的,此外,對其他時段啟停機過程的振動也進行了分析,可以得到類似的結論.

4? ?振動影響機制分析

4.1? ?運行模態識別

采用隨機子空間方法對實測振動數據進行自振頻率和阻尼比識別.啟停機工況下的一階自振頻率和阻尼比,如圖10和11所示.這里的阻尼比ξtotal,主要包括結構阻尼比ξstruct和氣動阻尼比ξaero兩部分[4]:

從圖10和圖11可以看出,該風力機的一階自振頻率隨運行工況的變化有微幅變化.停機狀態(1和6區域)對應的一階自振頻率最小,1區域均值為0.425 Hz,6區域為0.423 Hz,而其他區域的自振頻率都略有增大,這是由于葉輪轉動時,在重力和離心力的作用下導致葉輪的剛度增加,產生離心剛化效應,從而影響葉輪自身的動力特性[19].

識別得到的阻尼比離散性較大,但是從總體趨勢上還是有一些規律:2和5區域同為并網轉速狀態,2區域的阻尼比較小,大多低于2%,圖8中2區域的振動速度較大;而5區域的阻尼比較大,大都大于2%,圖9中5區域的振動速度較小,這也與文獻[3]得到的結論基本一致,阻尼比是影響結構振幅的重要因素之一.1和6區域(停機狀態)對應的阻尼比分別為0.548%和0.237%,此時葉輪轉速很低,為靜止或者空轉狀態,氣動阻尼可以忽略,可近似看做結構阻尼比.其他啟停機過程的阻尼比也可以得到類似的結論.

4.2? ?Sommerfeld效應分析

風力機在啟停機過程中不同的振動現象還可以從Sommerfeld效應得到進一步解釋[20-21].

通常,把不受系統響應影響的能量源稱為理想的能量源,隨著系統動態響應而發生變化的能量源為非理想能量源.Arnold Sommerfeld在研究一個由不平衡轉子和柔性支撐組成的振動系統時,首次觀察到這種非理想現象[21],該能量源的一部分能量被用于增加該系統的振動而不是增加轉子轉速,如圖12所示.

從圖12可以看到,當能量供給不斷增加(例如:風速逐漸增大),振幅由a開始增大到b,接近基頻(共振轉速Ω),這時,隨著能量供給進一步增加,轉子轉速增大不多(增加到Ω1),但振幅卻從b增加到較大的值c.之后,如果進一步增加能量供給,轉速增大到更高的值Ω2,更加遠離共振轉速Ω,振幅c降低到一個較低的值d.這種a→b→c→d的發展過程對應于風力機的啟動過程.當能量供給從大到小時(例如:風速逐漸減小),也可以觀察到類似的現象,但過渡路徑是d→e→b→a,對應于風力機的停機過程[21].這也解釋了圖8(c)和圖9(c)中塔頂振動速度與葉輪轉速之間的關系,在并網轉速附近,當轉速從小到大變化時,振動速度有突然增大再減小的現象,但是轉速從大到小時,振動速度變化不大.

此外,從圖12還可以看出,外部能量的大小決定是否有可能快速通過共振區域.當風速達到風力機的切入風速,但增加緩慢時,隨著葉輪過槳頻率接近塔架自振頻率,結構發生共振會消耗風能,葉輪不可能進一步加速轉動.如圖7,隨著風速從2 m/s增大到5 m/s,葉輪轉速基本沒有增大,約為10 r/min左右.若風速快速增大,有足夠的能量輸入,則能夠驅動葉輪達到更高的轉速從而穿越共振區[22].

5? ?結? ?論

由于風力機葉輪在啟停機過程中,不可避免地要經過并網轉速,風力機在接近并網轉速時,過槳頻率3f與結構的一階自振頻率很接近,會引起瞬時共振.為了了解此時風力機的振動特性,本文依據大量實測數據,采用運行模態分析方法,分析了風力機在啟停機過程中的振動和模態參數特征,并通過Sommerfeld效應對其進行解釋.主要結論如下:

1)對風力機塔頂的振動進行了長期監測,發現該風力機在并網轉速附近運轉的時段較多,當風力機啟動經過并網轉速時塔頂振動明顯加強,風力機停機經過并網轉速時,振幅變化不明顯.

2)基于數據驅動的SSI方法識別了風力機在啟停機過程中一階自振頻率和阻尼比.風力機的一階自振頻率隨運行工況不同而微幅變化,停機狀態的一階自振頻率最小;葉輪轉動時,由于離心剛化效應,自振頻率略有增大.風力機啟動經過并網轉速時,阻尼比較小,對應的振動速度較大;風力機停機經過并網轉速時,阻尼比較大,對應的振動速度較小;在相同的葉輪轉速下,振幅與阻尼比成反比.

3)Sommerfeld效應可以解釋風力機在啟停機過程中不同的振動現象.在風力機啟動過程中,隨著葉輪轉速增加接近共振頻率,風力機結構的共振會消耗風能,葉輪轉速增加緩慢,但風力機結構振幅增大.當有足夠的風能輸入,才能驅動葉輪達到更高的轉速從而穿越共振.

4)在進行風力機塔架設計時,應開展塔架基頻優化,避開1f 和3f頻率帶10%以上,并采取合適的控制策略,使風力機快速穿越共振區,以減少可能出現的共振現象.

參考文獻

[1]? ? 柯世堂,曹九發,王瓏,等. 風力機塔架-葉片耦合模型風致響應時域分析[J]. 湖南大學學報(自然科學版),2014,41(4):87—93.

KE S T,CAO J F,WANG L,et al. Time-domain analysis of the wind-induced responses of the coupled model of wind turbine tower-blade coupled system[J]. Journal of Hunan University(Natural Sciences),2014,41(4):87—93.(In Chinese)

[2]? ? DNVGL-ST-0126? Support structures for wind turbines[S]. Norway:Det Norske Veritas,2016:35—36.

[3]? ? DEVRIENDT C,JORDAENS P J,SITTER G D,et al. Damping estimation of an offshore wind turbine on a monopile foundation[J]. Iet Renewable Power Generation,2013,7(4):401—412.

[4]? ? HU W H,SEBASTIAN T,ROLF G R,et al. Vibration-based structural health monitoring of a wind turbine system. Part I:Resonance phenomenon[J]. Engineering Structures,2015,89:260—272.

[5]? ? SHIRZADEH R,DEVRIENDT C,BIDAKHVIDI M A,et al. Experimental and computational damping estimation of an offshore wind turbine on a monopile foundation[J]. Journal of Wind Engineering & Industrial Aerodynamics,2013,120:96—106.

[6]? ? TEMPEL J V D. Design of support structures for offshore wind turbines[D]. Holland:Offshore Engineering and Wind Energy Sections,Technische Universiteit Delft,2006:51—58.

[7]? ? DAMGAARD M,IBSEN L B,ANDERSEN L V,et al. Cross-wind modal properties of offshore wind turbines identified by full scale testing[J]. Journal of Wind Engineering and Industrial Aerodynamics,2013,116(5):94—108.

[8]? ? VERSTEIJLEN W G,METRIKINE A,HOVING J S,et al. Estimation of the vibration decrement of an offshore wind turbine support structure caused by its interaction with soil[J]. Knowledge Technology & Policy,2011,27(30):290—300.

[9]? ? HACKELL M W,ROLFES R. Monitoring a 5 MW offshore wind energy converter condition parameters and triangulation based extraction of modal parameters[J]. Mechanical Systems & Signal Processing,2013,40(1):322—343.

[10]? 劉佩,連鵬宇,張茉顏,等. 基于環境振動的某設置防震縫結構的動力特性[J]. 湖南大學學報(自然科學版),2017,44(1):95—101.

LIU P,LIAN P Y,ZHANG M Y,et al. Dynamic characteristics of a building with seismic joints based on ambient vibration[J]. Journal of Hunan University(Natural Sciences),2017,44(1):95—101. (In Chinese)

[11]? BROWNJOHN J M W,MAGALHAES F,CAETANO E,et al. Ambient vibration re-testing and operational modal analysis of the Humber Bridge[J]. Engineering Structures,2010,32(8):2003—2018.

[12]? SHI W,SHAN J,LU X. Modal identification of Shanghai World Financial Center both from free and ambient vibration response[J]. Engineering Structures,2012,36(4):14—26.

[13]? AMERI N,GRAPPASONNI C,COPPOTELLI G,et al. Ground vibration tests of a helicopter structure using OMA techniques[J]. Mechanical Systems and Signal Processing,2013,35(1/2):35—51.

[14]? IBRAHIM S R. Random decrement technique for modal identification of structures[J]. Journal of Space craft and Rockets,2012,14(11):696—700.

[15]? JAMES G H I,CARNE T H,LAUFFER J P. The natural excitation technique (NExT) for modal parameter extraction from operating wind turbine[J]. Nasa Sti/recon Technical Report N,1995,93(4):260—277.

[16]? PEETERS B,ROECK G D. Reference-based stochastic subspace identification for output-only modal analysis[J]. Mechanical Systems and Signal Processing,1999,13(6):855—878.

[17]? 任偉新.環境振動系統識別方法的比較分析[J]. 福州大學學報(自然科學版),2001,29(6):80—86.

REN W X. Comparison of system identification methods using ambient vibration measurements[J]. Journal of Fuzhou University(Natural Science),2001,29(6):80—86. (In Chinese)

[18]? BRINCKER R,ZHANG L M,ANDERSON P. Modal identification from ambient response using frequency domain decomposition[C]// Proceedings of the 18th International Modal Analysis Conference. USA:San Antonio,2000:625—630.

[19]? FUNG E H K,YAU D T W. Effects of centrifugal stiffening on the vibration frequencies of a constrained flexible arm[J]. Journal of Sound & Vibration,1999,224(5):809—841.

[20]? SAMANTARAY A K. Steady-state dynamics of a non-ideal rotor with internal damping and gyroscopic effects[J]. Nonlinear Dynamics,2008,56(4):443—451.

[21]? SOMMERFELD A. Beitr ge zum dynamischen ausbau der festigkeitslehre[J]. Physikal Zeitschr,1902:266—286.

[22]? BRASIL R M L R F,FEITOSA L C S,BALTHAZAR J M. A nonlinear and non-ideal wind generator supporting structure[J]. Applied Mechanics & Materials,2006,5/6:433—442.

猜你喜歡
風速振動區域
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
中立型Emden-Fowler微分方程的振動性
基于GARCH的短時風速預測方法
關于四色猜想
分區域
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
考慮風速分布與日非平穩性的風速數據預處理方法研究
主站蜘蛛池模板: 丰满人妻中出白浆| 伊伊人成亚洲综合人网7777| 国产黑丝一区| 四虎在线观看视频高清无码 | 国产门事件在线| 青青草一区| 国产成年女人特黄特色毛片免| 爱色欧美亚洲综合图区| 女同国产精品一区二区| 精品欧美一区二区三区久久久| 欧美日本激情| 国产区91| 狠狠色噜噜狠狠狠狠奇米777| 亚洲AV无码乱码在线观看代蜜桃 | 国产精品人成在线播放| 亚洲一区无码在线| аⅴ资源中文在线天堂| 久久精品日日躁夜夜躁欧美| 亚洲αv毛片| 自偷自拍三级全三级视频 | 国产欧美专区在线观看| 日韩性网站| 尤物精品视频一区二区三区| 国产精品免费久久久久影院无码| 久久香蕉国产线看观看亚洲片| 日韩第一页在线| 国产高清在线观看91精品| 国产性生交xxxxx免费| 中文字幕天无码久久精品视频免费| 久久综合亚洲鲁鲁九月天| 亚洲成人精品| 欧洲极品无码一区二区三区| 色哟哟国产精品| 亚洲欧洲日韩久久狠狠爱| 人人澡人人爽欧美一区| 国产丝袜无码精品| 中文字幕第4页| 久热99这里只有精品视频6| 3p叠罗汉国产精品久久| 中文字幕久久精品波多野结| 欧美.成人.综合在线| 五月激情婷婷综合| 国产精品女主播| 国产精品自拍合集| 91无码国产视频| 国产在线自乱拍播放| 视频二区亚洲精品| 国产亚洲日韩av在线| 国产农村精品一级毛片视频| h网址在线观看| 自偷自拍三级全三级视频 | 亚洲AⅤ波多系列中文字幕| 久久婷婷国产综合尤物精品| 456亚洲人成高清在线| 69av免费视频| 成人一级免费视频| 超碰免费91| 亚洲AV成人一区国产精品| 2020国产精品视频| 色香蕉影院| 国产青榴视频在线观看网站| 国产成人精品午夜视频'| 国产丝袜第一页| 亚洲毛片一级带毛片基地| 国产一区二区丝袜高跟鞋| 中文字幕亚洲第一| 丰满的熟女一区二区三区l| 啪啪永久免费av| 亚洲成人免费看| 国产人人射| 欧美乱妇高清无乱码免费| 亚洲精品爱草草视频在线| 99青青青精品视频在线| 国产福利2021最新在线观看| 狠狠操夜夜爽| 亚洲swag精品自拍一区| 日韩av高清无码一区二区三区| 亚洲国产无码有码| 久久毛片网| 亚洲精品欧美日本中文字幕| 亚洲一欧洲中文字幕在线| 97狠狠操|