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

環境參數對風電齒輪箱傳動系統疲勞損傷的影響

2024-05-23 03:24:07伍源朱才朝譚建軍宋朝省張會陽
重慶大學學報 2024年3期

伍源 朱才朝 譚建軍 宋朝省 張會陽

收稿日期:2022-01-23

網絡出版日期:2022-05-07

基金項目:重慶市自然科學基金資助項目(cstc2020jcyj-msxmX0710);中央高校基本科研業務費資助項目(2020CDJ-LHSS-008,2021CDJCGJ008)。

作者簡介:伍源(1997—),男,碩士研究生,主要研究方向為風電機組疲勞損傷分析,(E-mail)2284391610@qq.com。

通信作者:朱才朝,男,教授,博士生導師,(E-mail)cczhu@cqu.edu.cn。

摘要:在風電機組全壽命周期內,長期風速概率分布會使風電齒輪箱傳動系統動載荷出現隨機特性,影響其疲勞損傷預估精度。筆者提出了一種考慮長期風速概率分布特征的風電齒輪箱傳動系統疲勞損傷預估方法,通過建立大功率海上風電機組OpenFAST-SIMPACK聯合仿真模型,計算不同平均風速與湍流強度組合工況下的風電齒輪箱傳動系統齒輪短期疲勞損傷,進而采用代理模型技術重構“平均風速、湍流強度-短期疲勞損傷”映射關系,預測齒輪長期疲勞損傷。研究結果表明:風電齒輪箱傳動系統低速級太陽輪容易發生接觸疲勞失效;在額定風速以下,低速級太陽輪短期疲勞損傷與平均風速呈正相關,在額定風速附近,平均風速與湍流強度的隨機特性均會增大其長期疲勞損傷不確定性,增大其疲勞失效風險。

關鍵詞:風電齒輪箱;傳動系統;疲勞損傷;隨機特性

中圖分類號:TH113 ?????????文獻標志碼:A ??????????文章編號:1000-582X(2024)03-132-13

大功率海上風電機組是目前最有效地開發海上風能資源的重大海洋工程裝備之一。風電齒輪箱是風電機組中傳遞力與運動的關鍵傳動裝置。據統計,海上風電齒輪箱故障率比陸上高5%以上[1]。這是由于海面粗糙度比陸地小,海表風速大、季節性變化明顯,因而由風剪切效應造成的作用在風電齒輪箱的氣動載荷波動幅值大、隨機性強,造成海上風電齒輪箱疲勞失效間隔時間短[2]。同時,海上風電裝備運維成本比陸上高10%以上[1],對海上風電齒輪箱可靠性設計提出了更高要求。因此,開展考慮長期風速概率分布特征的風電齒輪箱傳動系統疲勞損傷隨機特性分析對指導海上風電齒輪箱設計、提高其可靠性具有重要意義。

近年來,國內外學者對風電齒輪箱傳動系統疲勞損傷開展了大量研究。Dong等[3]分別建立風電機組整機動力學模型和齒輪箱動力學模型,利用解耦方法計算了不同風速作用下風電機組全局氣動載荷,分析了不同氣動載荷作用下風電齒輪箱傳動系統齒輪接觸應力概率分布特征;向東等[4]建立隨機風載作用下風電齒輪箱傳動系統短期疲勞損傷計算模型,分析了隨機風速下各齒輪的短期彎曲疲勞與接觸疲勞損傷;Nejad等[5-7]針對風電齒輪箱傳動系統疲勞損傷分析開展了大量研究,包括考慮長期風速概率分布的10 MW與5 MW風電齒輪箱傳動系統齒輪長期疲勞損傷計算模型、考慮隨機波浪作用下4種浮式平臺對風電齒輪箱傳動系統長期疲勞損傷影響等;熊中杰[8]建立考慮時變嚙合剛度和阻尼的齒輪箱動力學模型,分析不同風速、不同載荷類型下的風電齒輪箱疲勞損傷規律。文獻[3-8]報道了有效預測風電齒輪箱傳動系統長期疲勞損傷的方法,通過等間距離散平均風速,計算各風速工況下的齒輪短期疲勞損傷后加權求和得到長期疲勞損傷,但離散間距取值主要依靠設計者經驗,離散后所需仿真的環境參數組合工況數量多,計算耗時。

為了提高傳動系統長期疲勞損傷的計算效率,部分學者將多項式響應面(polynomial response surface,PRS)、克里金(kriging,KRG)、人工神經網絡(artificial neural network,ANN)和徑向基函數(radial basis function,RBF)等代理模型方法引入到復雜系統動力學分析中[9-12]。Juan等[13]建立了環境參數與風電機組等效疲勞載荷和發電量之間的PRS代理模型,并開展環境參數全局靈敏度分析;Li等[14]考慮風浪等環境工況相關性,利用KRG代理模型預測漂浮式風電機組系泊纜、塔架底座和頂部橫截面的長期疲勞損傷;Zhang等[15]利用ANN和KRG等代理模型計算動態環境載荷下浮式風電機組系泊線的疲勞失效概率;Wilkie等[16]利用高斯過程回歸代理模型預測了海上風電機組疲勞損傷、疲勞可靠性等性能,分析風和波浪等環境參數對性能指標的影響。文獻[13-16]報道了高效的風電機組疲勞損傷分析的計算方法,但忽略了平均風速與湍流強度概率分布對風電齒輪箱傳動系統疲勞損傷隨機特性的影響。

筆者考慮隨機風速與風電機組拓撲結構,建立大功率海上風電機組OpenFAST-SIMPACK聯合仿真模型;考慮長期風速概率分布特征,利用Copula函數構建平均風速-湍流強度聯合概率分布,建立平均風速-湍流強度關鍵工況集合,計算對應的風電齒輪箱傳動系統齒輪短期疲勞損傷;采用代理模型技術重構“平均風速、湍流強度-疲勞損傷”映射關系,預測齒輪長期疲勞損傷。

1海上風電機組OpenFAST-SIMPACK聯合仿真建模

1.1風電機組結構及運行原理

圖1所示為某型風電機組齒輪箱傳動系統結構及原理。海上風電機組主要由風輪、主軸、齒輪箱、發電機、塔筒及塔架等組成。在隨機風速作用下,葉片將風能轉變為機械能,并通過主軸驅動齒輪箱轉動。齒輪箱主要包括低速級、中間級和高速級三級傳動結構,其中低速級和中間級為行星斜齒輪傳動,高速級為平行軸斜齒輪傳動。發電機通過聯軸器與高速軸相連。海上風電機組基本設計參數如表1所示,風電齒輪箱主要設計參數如表2所示。

1.2OpenFAST-SIMPACK聯合仿真模型

圖2所示為大功率海上風電機組OpenFAST-SIMPACK聯合仿真模型,主要包括2類不同層級的子模型,即風電機組整機全局耦合模型與風電齒輪箱傳動系統動力學模型。首先,利用OpenFAST[17]建立大功率海上風電機組整機全局耦合模型,控制策略采用變速-變槳控制[18];為了提高計算效率,風電齒輪箱傳動系統簡化為傳動比。然后,利用SIMPACK[19]建立風電齒輪箱傳動系統動力學模型,其中齒輪副采用切片理論建模、傳動軸與行星架簡化為剛體,軸承利用6×6剛度矩陣模擬;發電機和聯軸器采用集中質量單元建模[20]

為了實現OpenFAST-SIMPACK聯合仿真,利用Matlab動態修改OpenFAST風電機組整機全局耦合模型工況.inp文件,計算任意給定平均風速-湍流強度環境參數下的風電機組輪轂處6自由度氣動載荷;利用Matlab將氣動載荷文件格式轉為SIMPACK風電齒輪箱傳動系統動力學模型載荷.afs文件,并同時調用其宏命令.sjs文件進行仿真,計算風電齒輪箱傳動系統齒輪副動態嚙合力。

2考慮長期風速分布的風電齒輪箱傳動系統疲勞損傷預估

2.1長期風速概率分布

根據IEC61400-1標準[21],將10 min時序風速的平均值作為平均風速,其概率密度函數通常服從威布爾分布

2.2關鍵環境工況

為了從中高效地抽選權重占比較高的平均風速-湍流強度環境參數組合工況,結合拉丁超立方抽樣法[27]與最大差異法[28]兩者優點,建立圖3所示關鍵環境工況分析流程,以保證在抽樣樣本數量較少的條件下,同時保持各抽樣維度均勻性與差異最大性,具體步驟如下。

2.3齒輪長期疲勞損傷

3結果討論與分析

根據中國某海上風電場歷年風速統計數據,建立平均風速-湍流強度聯合概率密度函數,并選取關鍵環境工況,分別計算各關鍵環境工況下風電齒輪箱傳動系統齒輪副動態嚙合力(各工況仿真時間100 s),并計算齒輪短期和長期疲勞損傷;最后分析環境參數對齒輪短期和長期疲勞損傷的影響規律。

3.1環境參數對風電齒輪箱傳動系統齒輪短期疲勞損傷影響

圖6與圖7分別為平均風速概率密度函數和湍流強度概率密度函數,其中平均風速威布爾分布函數為,湍流強度伽馬分布函數為。

圖8所示為平均風速與湍流強度之間的相關性。從圖中可以看出,平均風速與湍流強度彼此之間存在明顯的“上尾相關性”,符合Gumbel Copula函數特征,因此,根據式(3)和(4)可得平均風速-湍流強度聯合概率密度函數,如圖9所示。

圖10所示為平均風速、湍流強度分別對低速級太陽輪接觸與彎曲應力均值與方差的影響。從圖10(a)中可以看出,當平均風速低于額定風速時,低速級太陽輪接觸與彎曲應力均值隨平均風速增加而增大;當平均風速高于額定風速時,接觸與彎曲應力均值基本穩定;低速級太陽輪接觸應力均值遠大于彎曲應力均值。從圖10(b)中可以看出,接觸應力方差在10~40 MPa波動,彎曲應力方差在0~25 MPa波動,總體上,低速級太陽輪接觸應力方差大于彎曲應力方差。

圖11所示為在額定風速(平均風速9.58 m/s,湍流強度0.11)作用下風電齒輪箱傳動系統各級齒輪短期接觸與彎曲疲勞損傷。可以看出,低速級太陽輪的短期接觸疲勞損傷最大,同時各級齒輪短期接觸疲勞損傷均大于彎曲疲勞損傷。

圖12所示為平均風速、湍流強度對低速級太陽輪短期接觸與彎曲疲勞損傷的影響。從圖12(a)中可知,當平均風速低于額定風速時,低速級太陽輪短期接觸與彎曲疲勞損傷隨平均風速增加而增大;當平均風速高于額定風速時,平均風速對其影響較小。其主要原因是當風速低于額定風速時,風電機組以追蹤最佳風能利用系數運行,載荷變化明顯;當風速高于額定風速時,由于變槳系統限制了風電機組功率超發,載荷較為穩定。從圖12(b)中可以看出,當平均風速低于8.0 m/s時,低速級太陽輪短期接觸與彎曲疲勞損傷主要受平均風速影響,湍流強度影響較小;當平均風速在8.0~9.5 m/s時,其損傷會隨著湍流強度的增加而增大;當平均風速高于9.5 m/s時,由于控制策略作用,湍流強度對其損傷影響較小。

3.2環境參數對風電齒輪箱傳動系統齒輪長期疲勞損傷影響

為了驗證本文方法計算結果的精度與效率,將本文方法與常規方法計算的齒輪長期疲勞損傷進行對比分析,如式(15)所示。在常規計算方法中,在平均風速3~25 m/s、湍流強度0~0.25的區間內,以平均風速間隔0.75 m/s、湍流強度間隔0.05選取115個環境參數組合工況,直接根據概率累加計算風電齒輪箱傳動系統齒輪長期疲勞損傷。

圖13所示為本文方法與常規方法計算的風電齒輪箱傳動系統齒輪長期疲勞損傷計算結果與誤差對比。從圖中可以看出,本文方法與常規方法計算的結果誤差小于8%。整個仿真分析在一臺CPU型號為i7-4790K,主頻為4.00 GHz的臺式計算機上完成,其中常規方法計算耗時約72 h,而本文方法計算僅耗時約20?h,效率提高260%。

圖14所示為平均風速、湍流強度對低速級太陽輪長期接觸與彎曲疲勞損傷隨機特性的影響。其中,黃色陰影部分表示當平均風速確定,湍流強度概率分布如圖14(a)時,或當湍流強度確定,平均風速概率分布如圖14(b)時,95%置信區間內的齒輪長期疲勞損傷。從圖14(a)中可以看出,低速級太陽輪長期接觸與彎曲疲勞損傷隨著平均風速的增加而先增后減,同時由于湍流強度的隨機性,尤其當平均風速位于10~11 m/s時齒輪長期疲勞損傷出現明顯的隨機不確定性。從圖14(b)中可以看出,低速級太陽輪長期接觸與彎曲疲勞損傷隨湍流強度的增加也呈現出先增后減的趨勢;當湍流強度在0.14附近時齒輪長期疲勞損傷最大,同時其隨機不確定性也越明顯,疲勞失效風險增大。

4結??論

筆者考慮長期風速概率分布特征,建立了大功率海上風電機組OpenFAST-SIMPACK聯合仿真模型,通過關鍵環境工況分析,計算對應的風電齒輪箱傳動系統齒輪短期疲勞損傷,基于代理模型預測齒輪長期疲勞損傷,得出結論如下:

1) 通過關鍵環境工況選取與基于代理模型的風電齒輪箱傳動系統齒輪長期疲勞損傷計算,相對于常規方法,其計算效率可以提高260%,且誤差小于8%。

2) 風電齒輪箱傳動系統各級齒輪接觸疲勞損傷均高于彎曲疲勞損傷,其中低速級太陽輪接觸疲勞損傷最大,易發生疲勞失效。

3) 在額定風速以下時,低速級太陽輪短期疲勞損傷與平均風速呈正相關,而在額定風速附近時,其主要受湍流強度影響;平均風速與湍流強度的隨機特性會增大低速級太陽輪長期疲勞損傷不確定性。

參考文獻

[1]??李垚, 朱才朝, 陶友傳, 等. 風電機組可靠性研究現狀與發展趨勢[J]. 中國機械工程, 2017, 28(9): 1125-1133.

Li Y, Zhu C C, Tao Y C, et al. Research status and development tendency of wind turbine reliability[J]. China Mechanical Engineering, 2017, 28(9): 1125-1133. (in Chinese)

[2]??王磊. 海上風電機組系統動力學建模及仿真分析研究[D]. 重慶: 重慶大學, 2011.

Wang L. Study on systematic dynamic model and simulation for offshore wind turbine[D]. Chongqing: Chongqing University, 2011. (in Chinese)

[3]??Dong W B, Xing Y H, Moan T, et al. Time domain-based gear contact fatigue analysis of a wind turbine drivetrain under dynamic conditions[J]. International Journal of Fatigue, 2013, 48: 133-146.

[4]??向東, 蔣李, 沈銀華, 等. 風電齒輪箱在隨機風載下的疲勞損傷計算模型[J]. 振動與沖擊, 2018, 37(11): 115-123.

Xiang D, Jiang L, Shen Y H, et al. Fatigue damage calculation model for wind turbine gearboxes under random wind loads[J]. Journal of Vibration and Shock, 2018, 37(11): 115-123. (in Chinese)

[5]??Nejad A R, Gao Z, Moan T. On long-term fatigue damage and reliability analysis of gears under wind loads in offshore wind turbine drivetrains[J]. International Journal of Fatigue, 2014, 61: 116-128.

[6]??Nejad A R, Bachynski E E, Kvittem M I, et al. Stochastic dynamic load effect and fatigue damage analysis of drivetrains in land-based and TLP, spar and semi-submersible floating wind turbines[J]. Marine Structures, 2015, 42: 137-153.

[7]??Wang S S, Nejad A R, Moan T. On design, modelling, and analysis of a 10-MW medium-speed drivetrain for offshore wind turbines[J]. Wind Energy, 2020, 23(4): 1099-1117.

[8]??熊中杰. 隨機風速下風力發電機組齒輪箱疲勞斷裂壽命研究[D]. 南京: 南京理工大學, 2019.

Xiong Z J. Research on fatigue fracture life of wind turbine gearbox under random wind[D]. Nanjing: Nanjing University of Science and Technology, 2019. (in Chinese)

[9]??Heidebrecht A, MacManus D G. Surrogate model of complex non-linear data for preliminary nacelle design[J]. Aerospace Science and Technology, 2019, 84: 399-411.

[10]??Palmer K, Realff M. Metamodeling approach to optimization of steady-state flowsheet simulations: model generation[J]. Chemical Engineering Research and Design, 2002, 80(7): 760-772.

[11]??Jia Z Y, Davis E, Muzzio F J, et al. Predictive modeling for pharmaceutical processes using kriging and response surface[J]. Journal of Pharmaceutical Innovation, 2009, 4(4): 174-186.

[12]??Milovanovi? S, von Sydow L. A high order method for pricing of financial derivatives using radial basis function generated finite differences[J]. Mathematics and Computers in Simulation, 2020, 174: 205-217.

[13]??Murcia J P, Réthoré P E, Dimitrov N, et al. Uncertainty propagation through an aeroelastic wind turbine model using polynomial surrogates[J]. Renewable Energy, 2018, 119: 910-922.

[14]??Li X A, Zhang W. Probabilistic fatigue evaluation of floating wind turbine using combination of surrogate model and copula model[C]//Proceedings of the AIAA Scitech 2019 Forum, San Diego, California. Reston, Virginia: AIAA, 2019: AIAA2019-0247.

[15]??Zhao Y L, Dong S. Probabilistic fatigue surrogate model of bimodal tension process for a semi-submersible platform[J]. Ocean Engineering, 2021, 220: 108501.

[16]??Wilkie D, Galasso C. Impact of climate-change scenarios on offshore wind turbine structural performance[J]. Renewable and Sustainable Energy Reviews, 2020, 134: 110323.

[17]??Jonkman B, Jonkman J. FAST v8. 16.00 a-bjj[EB/OL]. 2016-07-16[2021-12-20]. https://openfast.readthedocs.io/en/v3.5.2/_downloads/5f2ddf006568adc9b88d8118dc3f1732/FAST8_README.pdf.

[18]??陳旭, 朱才朝, 宋朝省, 等. 緊急停機工況下風力發電機系統動態特性分析[J]. 機械工程學報, 2019, 55(5): 82-88.

Chen X, Zhu C C, Song C S, et al. Dynamic characteristics analysis of wind turbine under emergency shutdown events[J]. Journal of Mechanical Engineering, 2019, 55(5): 82-88. (in Chinese)

[19]??陳巖松, 朱才朝, 譚建軍, 等. 多工況下兆瓦級海上風電齒輪箱均載性能優化設計[J]. 重慶大學學報, 2022, 45(09): 1-14.

Chen Y S, Zhu C C, Tan J J, et al. Optimal design of load sharing performance of megawatt level offshore wind turbine gearbox under multi-operating conditions[J]. Journal of Chongqing University, 2022, 45(09): 1-14.(in Chinese)

[20]??劉華朝, 朱才朝, 柏厚義. 輪齒修形對兆瓦級風電齒輪箱NVH性能的影響[J]. 振動與沖擊, 2016, 35(24): 158-163, 188.

Liu H C, Zhu C C, Bai H Y. The effect of gear modification on the NVH characteristics of a megawatt level wind turbine gearbox[J]. Journal of Vibration and Shock, 2016, 35(24): 158-163, 188. (in Chinese)

[21]??International Electrotechnical Commission. Wind turbines -?Part 1: design requirements : IEC 61400-1:2005 [S]. IEC, 2006.

[22]??Li H X, Cho H, Sugiyama H, et al. Reliability-based design optimization of wind turbine drivetrain with integrated multibody gear dynamics simulation considering wind load uncertainty[J]. Structural and Multidisciplinary Optimization, 2017, 56(1): 183-201.

[23]??龔偉俊, 李為相, 張廣明. 基于威布爾分布的風速概率分布參數估計方法[J]. 可再生能源, 2011, 29(6): 20-23.

Gong W J, Li W X, Zhang G M. The estimation algorithm on the probabilistic distribution parameters of wind speed based on Weibull distribution[J]. Renewable Energy Resources, 2011, 29(6): 20-23. (in Chinese)

[24]??Sklar A. Fonctions de repartition an dimensions et leurs marges[J]. Publ. inst. statist. univ. Paris, 1959, 8: 229-231.

[25]??涂志斌, 黃銘楓, 樓文娟, 等. 基于Copula函數的風浪多方向極限狀態曲線[J]. 振動與沖擊, 2021, 40(14): 1-9, 46.

Tu Z B, Huang M F, Lou W J, et al. Dimensional environmental contour lines of the wind and wave based on Copula functions[J]. Journal of Vibration and Shock, 2021, 40(14): 1-9, 46. (in Chinese)

[26]??侯亞楠. Copula函數的估計及其應用[D]. 武漢: 華中科技大學, 2013.

Hou Y N. The estimation and application of copula function[D]. Wuhan: Huazhong University of Science and Technology, 2013. (in Chinese)

[27]??Khowaja K, Shcherbatyy M, H?rdle W K. Surrogate models for optimization of dynamical systems[EB/OL]. 2021: arXiv: 2101.10189. https://arxiv.org/abs/2101.10189.pdf.

[28]??Martini M, Guanche R, Armesto J A, et al. Met-ocean conditions influence on floating offshore wind farms power production[J]. Wind Energy, 2016, 19(3): 399-420.

[29]??周金宇, 謝里陽, 韓文欽, 等. 基于Nataf變換的載荷相關系統風險預測方法[J]. 機械工程學報, 2009, 45(8): 137-141.

Zhou J Y, Xie L Y, Han W Q,et al. Method for syetem risk prediction with load dependency based on nataf transformation[J]. Journal of Mechanical Engineering, 2009, 45(8): 137-141. (in Chinese)

[30]??張立波, 程浩忠, 曾平良, 等. 基于Nataf逆變換的概率潮流三點估計法[J]. 電工技術學報, 2016, 31(6): 187-194.

Zhang L B, Cheng H Z, Zeng P L, et al. A three-point estimate method for solving probabilistic load flow based on inverse Nataf transformation[J]. Transactions of China Electrotechnical Society, 2016, 31(6): 187-194. (in Chinese)

[31]??International Organization for Standardization. Calculation of load capacity of spur and helical gears -?Part 6: calculation of service life under variable load: ISO 6336-6:2006 [S]. Switzerland: International Organization for Standardization, 2006.

[32]??International Organization for Standardization. Calculation of load capacity of spur and helical gears -?Part 2: calculation of surface durability (pitting): ISO 6336-2:2006 [S]. Switzerland: International Organization for Standardization, 2006.

[33]??International Organization for Standardization. Calculation of load capacity of spur and helical gears -?Part 3: calculation of tooth bending strength: ISO 6336-3:2006 [S]. Switzerland: International Organization for Standardization, 2006.

[34]??Miner M A. Cumulative damage in fatigue[J]. Journal of Applied Mechanics, 1945, 12(3): A159-A164.

[35]??International Organization for Standardization. Calculation of load capacity of spur and helical gears -?Part 5: strength and quality of materials: ISO 6336-5:2003 [S]. Switzerland: International Organization for Standardization, 2003.

[36]??龍騰, 劉建, Wang G G, 等. 基于計算試驗設計與代理模型的飛行器近似優化策略探討[J]. 機械工程學報, 2016, 52(14): 79-105.

Long T, Liu J, Wang G G, et al. Discuss on approximate optimization strategies using design of computer experiments and metamodels for flight vehicle design[J]. Journal of Mechanical Engineering, 2016, 52(14): 79-105. (in Chinese)

[37]??Okpokparoro S, Sriramula S. Uncertainty modeling in reliability analysis of floating wind turbine support structures[J]. Renewable Energy, 2021, 165: 88-108.

[38]??Rashki M, Azarkish H, Rostamian M, et al. Classification correction of polynomial response surface methods for accurate reliability estimation[J]. Structural Safety, 2019, 81: 101869.

[39]??Wang H, Li W B, Qian Z H, et al. Reconstruction of wind pressure fields on cooling towers by radial basis function and comparisons with other methods[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2021, 208: 104450.

[40]??Bhosekar A, Ierapetritou M. Advances in surrogate based modeling, feasibility analysis, and optimization: a review[J]. Computers & Chemical Engineering, 2018, 108: 250-267.

(編輯??呂建斌)

主站蜘蛛池模板: 亚洲人成网线在线播放va| 91网站国产| 久久综合干| a国产精品| 自拍亚洲欧美精品| 99久久精品免费看国产电影| 国产午夜看片| 国产毛片不卡| 亚洲动漫h| 日韩在线播放欧美字幕| 色网在线视频| 国产亚洲精品无码专| 人妻丰满熟妇αv无码| 天天综合网在线| 国产丝袜丝视频在线观看| 国产无码精品在线| 欧美高清视频一区二区三区| 国产欧美精品一区二区| 蜜桃视频一区| 欧美精品在线免费| 亚洲三级a| 色哟哟国产成人精品| 国产经典三级在线| 全色黄大色大片免费久久老太| 正在播放久久| 成年片色大黄全免费网站久久| 欧美在线国产| 国产一区免费在线观看| 国产Av无码精品色午夜| 久久久久亚洲精品成人网| 免费不卡在线观看av| 国产综合精品一区二区| 国产精品亚欧美一区二区三区| 国产91线观看| 中文字幕人妻无码系列第三区| 亚洲中文精品人人永久免费| 四虎成人在线视频| 欧美日一级片| 伊人大杳蕉中文无码| 日本久久久久久免费网络| 国产网站在线看| 天天综合亚洲| 午夜免费小视频| 国产系列在线| 日韩毛片免费| 国产99视频在线| 精品一区二区三区波多野结衣| 亚洲欧美在线综合图区| 亚洲国产系列| 中文字幕人成乱码熟女免费| 国产精品综合久久久 | 99久久国产综合精品女同| 婷婷五月在线| 国产精品永久不卡免费视频| 波多野结衣视频网站| 国产亚洲精品97在线观看| 欧美日韩精品综合在线一区| 在线观看无码a∨| 婷婷色丁香综合激情| 综合色在线| 激情六月丁香婷婷| 欧美精品v欧洲精品| 久久综合伊人 六十路| 亚洲第一视频区| 热99精品视频| 精品人妻无码中字系列| 欧美性久久久久| 久久国产拍爱| 亚洲第一成网站| 亚洲国产AV无码综合原创| 亚洲国产成人麻豆精品| 成人欧美在线观看| 男女男免费视频网站国产| 国产精品爽爽va在线无码观看| 中文字幕一区二区人妻电影| 中国毛片网| 亚洲精品成人片在线观看| 日韩小视频在线观看| 亚洲第一中文字幕| 在线精品欧美日韩| 91精品国产自产在线老师啪l| 青青草91视频|