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

基于剖面分割的大型風力發電機葉片三維覆冰增長特性分析

2025-11-12 00:00:00范才進牛唯卓然陳沛龍劉磊毛先胤黃歡楊旗
重慶大學學報 2025年10期

中圖分類號:TM85 文獻標志碼:A 文章編號:1000-582X(2025)10-001-10

Analysis of three-dimensional icing characteristics of large wind turbine blades based on profile segmentation

FAN Caijin',NIU Wei2,ZHUO Ran',CHEN Peilong2,LIU Lei’,MAO Xianyin2, HUANG Huan2,YANG Qi2 (1.Electric Power Research Institute,China Southern Power Grid, Guangzhou 510663,P.R.China; 2.Key Laboratory of Ice Prevention amp; Disater Reducing of China Southern Power Grid Co.,Ltd.,Electric Power Research Institute of Guizhou Power Grid Co.,Ltd., Guiyang 55oo02,P.R.China)

Abstract: Ice accretion is a crucial factor affcting the safe and stable operation of wind turbines.Developing a numerical model for simulating ice formation on wind turbine blades is essential for predicting icing phenomena. Although the finite element method iscurrentlythe most widelyusedapproach,itiscomputationally intensiveand inefficient for large-scale applications. This study focuses on the blades of a 300kW wind turbine, employing a profilesegmentation method to investigate water droplet impact,freezing,and iceaccretion morphology changes on blade surfaces.A multiphase flow simulation model for air and liquidon the blade surface is developed,and formulas for local and overall collsion and freezing coefficients are derived.This approach enables characterization of overall water droplet impact and freezing behavior with reduced computational load.Results reveal that the water droplet collsion coefficient decreases gradually from the blade tip toward the root,with reductions exceeding 80% in the maximum values of β1 and a1 at approximately 0.5R. Maximum water droplet capture occurs near the blade tip ( 0.8R to 0.9R ),while significant ice accretion predominantly occurs between 0.5R and R .The overflow effect of the water film results in low freezing coeficients in the droplet colision zone buthighervalues in the overflowregion.Furthermore,closer tothe bladetip,ice growth exhibits greateriterative shape changes,and reduced linearity in the ice accretion rate.

Keywords: wind turbine blade; icing; water droplet collision; icing accretion morphology

近年來,風力發電已成為全球可再生能源的重要組成部分之一。截至2022年底,全球風力發電裝機容量已經達到 1000GW 左右,從發展趨勢來看,風力發電技術正在不斷改進和創新,例如增加風機高度、提高風機效率、降低成本以及增強可靠安全性等[13]。根據中國風能行業協會的數據,覆冰是風力發電面臨的主要問題之一,每年因此造成很大損失。在2018年的一次覆冰天氣中,中國多個省份的風力發電機組停機率高達 70% 以上,導致全國風電發電量減少了約2.5億 kW?h ,損失達1.7億元人民幣[45]。

為研究風力發電機葉片覆冰的防治方法,國內外學者對風機葉片覆冰形成機制及影響特性展開了大量研究。風力發電機葉片覆冰主要集中在前緣部分,主要由空氣中過冷卻水滴碰撞凍結而成。黎芷毓等在數值模擬 300kW 大型風力發電機葉片霧淞覆冰時,采用拉格朗日法計算了風機葉片表面的水滴碰撞系數。Shu等則基于FLUENT軟件采用歐拉法實現了對風機葉片表面水滴碰撞系數的求解,并通過人工覆冰試驗驗證了仿真計算結果,結果表明:在不同風速、溫度和水滴中值直徑條件下風機葉片表面覆冰范圍、形態差異明顯。

此外,FENSAP軟件也常用于風機葉片表面覆冰增長的數值模擬,并可得出水滴碰撞系數分布。Jin等9基于FENSAP對風機葉片翼型S826和S832進行了覆冰模擬,并通過風洞試驗進行驗證,發現在相同覆冰環境下,與S826相比,S832型風機葉片的冰形更為復雜,2種翼型氣動性能差異與覆冰厚度和類型相關。Ibrahim等[研究了旋轉風力發電機葉片上結冰的比例和相似性,使用FENSAPICE軟件進行數值CFD(computational fluid dynamics)結冰模擬,并測試了風機葉片翼型縮放方法對覆冰結果的影響。

為提高風機覆冰數值模擬的效率,解決覆冰計算中的一些特殊問題,許多研究者還編寫了結冰計算軟件。Wang等[針對風機葉片提出了一種改進的多重結冰數值計算模型,該模型整合了流場計算、水滴碰撞凍結過程,可用于解決風機葉片外邊流場的周期性變化對覆冰過程造成的影響。Son等[12則編寫了名為WISE的風機葉片三維結冰計算程序,并通過FENSAP驗證了程序計算結果。

為進一步揭示覆冰對風機功率輸出特性的影響,舒立春等[3-15]在人工氣候室和自然環境條件下研究了覆冰環境參數對小型風力發電機功率輸出特性的影響;胡琴等利用了輕黏土和膩子試驗模擬不同形態覆冰對其功率的影響特性,結果顯示:相對于流線冰,角狀冰會帶來更多的氣動載荷損失。李瀚濤等[1則考慮了覆冰粗糙度對風機氣動性能的影響,發現覆冰葉片氣動性能下降主要由葉片氣動外形變化導致,而受覆冰粗糙度的影響較小。

風力發電機工況特殊,當前對葉片覆冰的防治手段仍不完善,而基于環境條件對風機葉片覆冰增長速率進行仿真模擬是一項重要手段。應用較多的是利用CFD流體力學仿真計算軟件模擬風機覆冰涉及到的流場、水滴碰撞和凍結,但該過程較為復雜,且風力發電機尺寸大、網格數量多,對計算機性能和計算時間有一定要求,一定程度上制約了該技術在風機覆冰預測、預警和防治上的應用。

筆者以 300kW 風機葉片為研究對象,采用剖面分割法來研究風機葉片不同位置水滴碰撞和凍結特性,并將計算結果和其他研究者的試驗觀測結果進行對比,在規避大量有限元網格繪制和迭代計算量的同時獲得風機葉片表面較為準確的水滴碰撞、凍結、覆冰形態變化規律,為風力發電機葉片覆冰預測、預警和防治提供數據和方法參考。

1風力發電機葉片覆冰數值計算模型

如圖1所示,風力發電機葉片覆冰主要是由空氣中隨氣流運動的過冷卻水滴碰撞凍結形成的,而風機葉片外部空氣中的過冷卻水滴隨氣流運動形成氣-液二相流,水滴運動特性主要取決于氣流對水滴的作用力。

圖1風機葉片翼型外部流場劃分

Fig.1 Divisionof external flow field ofwind turbineblade surface icing

為獲得水滴在葉片表面的碰撞特性,采用歐拉法求解風力發電機葉片外部氣-液二相流分布(利用FLUENT軟件),該分布中氣流運動滿足質量、動量和能量守恒的非線性微分方程為

式中: ρa 是氣流密度, kg/m3 ua 是氣流速度, m/s;T 表示空氣靜態溫度, K;σij 是應力張量, Pa;Ea 和 Ha 是總能量和焓;t為時間, s;g 為重力加速度, m/s2;ka 熱傳導系數, W/(m?K);νi 為二相流速度, m/s;τij 為黏性應力張量, N/m2 水滴和氣流一樣被視為連續相,其運動的質量平衡和動量方程為

式中: uw 是水滴的速度, m/s;a 是空氣中水滴的體積分數; ρw 是氣流密度, kg/m3;Cp 是液滴阻力系數; Red 是水滴相對氣流運動的雷諾數; k 是運動質點慣性的Stokes數。

式中: μ 為空氣的運動黏度, m2/s;Rd 為水滴半徑, m 。

通過氣-液二相流求解后可獲得水滴在風機葉片表面的體積分數 α ,進而可通過式(7)獲得葉片表面的水滴局部碰撞系數 β? 和整體碰撞率 a1 分別為

式中: L 是二維葉片輪廓長度, m;H 代表葉片垂直于氣流方向的最大厚度,即葉片在氣流方向上投影的最大高度, m;n 為葉片表面水滴碰撞單位的法向向量; a 和 V 分別為無窮遠處空氣中水滴體積分數和風速;dL為沿著葉片二維輪廓的長度微元;水滴碰撞風機葉片后,可通過熱平衡和質量平衡方程求解其凍結比例。如圖2所示,覆冰在風機葉片表面的凍結過程主要獲得熱量為水滴碰撞動能和凍結過程釋放潛熱,損失的熱量包括熱輻射、蒸發或升華消耗熱量,對流熱量和加熱過冷卻水滴至凍結溫度所需熱量等,熱平衡方程為

Qc+Qe+Qw+Qr=Qm+Q

圖2風機葉片表面覆冰熱平衡示意圖

結合凍結質量平衡,碰撞水滴質量 M0 為凍結冰層質量 Mice 和未凍結水量 Mw 之和,可求解水滴局部凍結系數 β3 為[18]

式中: Sj 為風機葉片上目標單元面積, m2;hc 為對流換熱系數, W/(m2?K);Tw 為水滴凍結溫度, K;χ 為蒸發系數;

e(T) 為溫度 T 時凍結表面的飽和水汽壓, kPa;w 為空氣中液態水含量, kg/m3 : V 為葉片相對氣流運動速度;

Cw 為水的比熱容, J/(kg?°C);ε 為發射率, ε=0.95:σ?R 為斯蒂芬-玻耳茲曼常量; Lf 為冰的融化潛熱, J/kg 。

2 風力發電機葉片水滴碰撞特性

以 300kW 三葉片風力發電機葉片為研究對象,該風力發電機基本參數如表1所示。風力發電機旋轉過程中葉片各位置和氣-液二相流的相對速率很大,若將風機葉片剖切為多個截面(見圖3),剖切間距需要根據可接受的相鄰截面最大水滴碰撞系數誤差決定。根據風力發電機 V-n 運行曲線(見圖4)可計算各截面所處位置合成速度,進而按式(7)(8)計算水滴在葉片表面的局部碰撞系數 β? 和整體碰撞系數 a1

表1 300kW 風力發電機基本參數

圖3風機葉片剖切示意圖

Fig.3Slices of thewind turbine blade

圖4風力發電機運行曲線 Fig.4 Windturbineoperationcurve

如圖5所示,由于不同截面形態、大小及所處位置運行線速度不同,各個截面的水滴局部碰撞系數的分布差異明顯。對于局部碰撞系數 β? 值而言,各個截面最大水滴碰撞系數均位于截面的駐點處,沿著兩邊逐漸減小,截面位置越靠近葉尖,水滴碰撞系數 β1 最大值越大,在旋轉半徑為 14.44m 截面A處,駐點 β? 最大為0.8,而在旋轉半徑為 4m 的截面G處 β? 最大僅0.03。對于碰撞范圍而言,從截面A到截面G,水滴碰撞范圍先逐步增大,而后又逐漸縮小。由此可知:風機葉片靠近輪轂的位置由于截面尺寸較大、線速度較小則不易捕獲過冷卻水滴,而靠近葉尖位置往往更易捕獲空氣中的過冷卻水滴而形成覆冰。若將對每個截面水滴局部碰撞系數進行積分處理,可得到各截面的水滴整體碰撞系數 a1 ,如圖6所示,圖中 R 為葉片截面所處位置離葉根的長度。

圖5風機葉片各截面局部水滴碰撞系數 β1

Fig.5Localwaterdropletcollisionefficiencies β? ofblade slides

a?1 沿著整個葉片的分布主表現出以下特征: 1)a1 值均在葉尖處最大,沿著葉片向輪轂方向逐漸降低,但 a1 值的零點并不在輪轂位置 ;2)a1 值的零點位置主要受到水滴粒徑大小的影響,當水滴中值體積直徑(medianvolumediameter,MVD)較大時,該位置靠近輪轂,否則靠近葉尖。根據 a1 的變化規律,風機葉片的覆冰區域主要集中在 0.5R~R 范圍內。例如在MVD為 20μm 時,在2種工況下( V=5m/s 和 8m/s 風機葉片的 7~ 14.44m 位置才會有覆冰產生,其余位置無覆冰。

β? 和 a?1 反映了葉片局部或整個截面在迎風側投影范圍內水滴碰撞捕獲的比例,為進一步對比葉片各位置水滴捕獲量的差異,假設空氣中液態水含量 w=1g/m3 ,對單位時間1s內葉展方向 1m 厚各個截面葉片表面水滴質量捕獲速率 M 進行計算,得到結果如圖7所示。

圖6風機葉片各截面整體水滴碰撞系數 a1

圖7風機葉片各截面( 1m 厚)整體水滴質量捕獲速率 M

Fig.7The overall water droplet mass capture rate M ofblade slidesatathickness of 1m

不同于水滴碰撞系數,水滴質量捕獲速率M反映了葉片整體對過冷卻水滴的捕獲速率。根據2種運行條件下的計算結果, M 值沿著葉展方向從輪轂向葉尖呈現先增加后減小的規律,這表明:盡管葉尖處的水滴碰撞系數最大,但相同時間內葉尖不再是獲得過冷卻水滴最多的位置,而靠近葉尖的 R=12~13m 處捕獲的水滴量更多。該結論可為風機葉片電熱防冰線路的優化布置提供參考。

3風力發電機葉片覆冰模擬及增長規律

3.1風機葉片表面的水滴凍結特性

風機葉片表面水滴的碰撞系數(局部 β1 或整體 a1) 可用于計算覆冰水的獲得量,通過式(9)~(10)的求解可得到該覆冰水量的凍結比例,即凍結系數 a3 或 β3 ,由此可通過Makkonen模型獲得覆冰形態和速率dM/dt。

式中,ds為風機葉片表面覆冰計算針對目標單元的面積。

風機葉片表面凍結系數 β3 分布如圖8~9所示,圖中橫坐標 s 表示沿著翼型截面的距離。對比水滴碰撞區域可知,水滴實際凍結區域(即生成覆冰區域)要大于水滴碰撞區域,這主要是由于水滴捕獲后產生水膜,水膜在風力作用下向葉片兩側溢流造成的,越靠近外溢流邊界的區域,其水滴凍結系數 β3 越接近于1,表明水膜在溢流過程中被逐漸凍結,最終被完全凍結為冰層。而在水滴碰撞區域內, β3 最大值和 β1 最大值位置并不重合,也側面證明了當空氣中液態水含量較大時,葉片翼型駐點附件捕獲的水滴大部分并未直接凍結,而是以水膜形式流向兩側。

圖8風機葉片截面A和截面D處水滴凍結系數分布

圖9風機葉片截面A和截面 D 處 β1 和 β3 的對比結果 Fig.9 Comparison resultsof β1 and β3 (section A and section D)

此外,凍結系數 β3 值受到風速 ∣V? 溫度 T 和液態水含量w的影響,當T降低,w減小,水滴碰撞區域內的 β3 值會迅速增大,溢流區域會相應減小,當風速增大時, β3 值的變化則取決于對流換熱效率的增大程度。

3.2風機葉片覆冰形態和質量變化

由于風機葉片各個截面上水滴碰撞系數和凍結系數分布不均,相同時間內形成的覆冰形態和質量也差別巨大。設定迭代步長為 dt=15min,90min 后各截面覆冰形態如圖10所示,葉片整體三維形態與文獻[20-21]試驗觀測結果一致。對于葉尖處截面A,覆冰開始時,冰層分布較寬,包含前緣和截面下表面,隨著厚度的增加,覆冰逐漸集中于A截面前緣,覆冰整體形態粗糙且不均勻,極大改變了葉片截面翼型結構。隨著截面向葉根方向移動,覆冰增長速率降低,水滴碰撞凍結位置的冰層分布逐漸均勻化,截面G~I位置基本無覆冰。

圖10風機葉片各截面覆冰形態

Fig.10 Comparison of icing shapes on blade slides

在設定條件下各個截面覆冰質量隨覆冰時間增長情況如圖11所示。各個截面覆冰質量 Ma 基本隨時間線性增長,越靠近葉根的截面其 Ma-t 曲線線性度越好,反之越差。例如對截面A ,tlt;45min 時其覆冰速率基本保存不變, t=45min 覆冰速率突然減小,而后又增大,這主要是由于葉尖位置覆冰形態變化較大造成水滴碰撞捕獲量改變形成的。而對于截面E~F等,其本身截面尺寸較大,水滴碰撞凍結區域相對較小,冰層對其流場影響不大,因此覆冰速率小且穩定。

圖11風機葉片各截面覆冰質量(各截面均寬

Fig.11 Comparison of icingmasson blade slides(wideof

4結論

1)針對風力發電機葉片覆冰過程水滴碰撞、凍結物理過程,基于氣液二相運動和熱平衡方程建立了水滴碰撞和凍結數值計算模型。推導了風機葉片截面水滴整體碰撞系數計算式。

2)對于類似 300kW 大型風力發電機葉片可采用剖切法對葉片不同段二維截面水滴碰撞特性進研究,可大幅降低計算量。一般運行條件下風機葉片水滴碰撞系數由葉尖向葉根方向逐漸減小, β1 最大值和 α1 到葉中( 0.5R 處)均可降低 80% 以上。

3)風機葉片覆冰水滴捕獲量最大位置位于靠近葉尖的 0.8R~0.9R 處,而一般覆冰均主要分別在 0.5R~R 范圍內,葉片防冰需重點關注該位置內的覆冰增長情況,對于其他位置可不用采取防冰措施。

4)由于水膜溢流作用,風機葉片截面的水滴凍結系數呈現水滴碰撞區域小而溢流區域大的特點(接近于1),從而進一步擴大了覆冰初期葉片表面的覆冰增長范圍。

5)越靠近風機葉片葉尖,覆冰增長過程冰形迭代變化程度越大,覆冰速率線性度也越差,反之,對越靠近葉根方向的截面,冰層對截面形態的影響較小,覆冰速率也更為穩定。

參考文獻

[1]王月普.風力發電現狀與發展趨勢分析[J].電力設備管理,2020(11):21-22.

WangYP.Analysisofcurrentsituationanddevelopment trendofwind powergeneration[J].ElectricPowerEquipment Management, 2020(11): 21-22. (in Chinese)

[2]王孟夏,周生遠,楊明,等.計及海底電纜熱特性的可接納海上風電裝機容量評估方法[J].電力系統自動化,2021,45(6):195-202.

9 Wang MX,Zhou SY,Yang M,etal.Asessment methodforacceptable installedcapacityofoffshore windfarmsconsidering thermal characteristicsofsubmarine cables[J].AutomationofElectric Power Systems,2021,45(6):195-202.(inChinese)

[3]胡琴,楊秀余,梅冰笑,等.風力發電機葉片臨界防冰與融冰功率密度分析[J].中國電機工程學報,2015,35(9): 4997- 5002. Hu Q,Yang XY,MeiBX,etal.Analysisof treshold powerdensityforanti-icing and de-icing of wind turbine blade[J]. Proceedings of the CSEE,2015,35(19): 4997-5002.(in Chinese)

[4]ZhangYW,GuoWF,LiY,etal.Anexperimentalstudyoficingdistributiononasymmetricalaifoilforwindturbineblade in the offshore environmental condition[J]. Ocean Engineering,2023,273: 113960.

[5]胡琴,王歡,邱剛,等.風力發電機葉片覆冰量化分析及其應用[J].電工技術學報,2022,37(21):5607-5616. Hu Q,Wang H,QiuG,etal. Quantitative analysis of wind turbine blade icingand itsapplication[J].TransactionsofChina Electrotechnical Society,2022,37(21): 5607-5616.(in Chinese)

[6]VillapandoF,ReggioM,Iinca A.Predictionof iceaccretionandanti-icing heating poweron wind turbine bladesusing standard commercial software[J].Energy,2016,114:1041-1052.

[7]黎芷毓,蔣興良,韓興波,等.風力發電機葉片的霧淞覆冰數值模擬[J].哈爾濱工業大學學報,2022,54(3):155-162. LiZY,JiangXL,HanXB,etal.Numericalsimulationonimeicingofwindturbineblades[J].JouralofHarbinIstituteof Technology,2022,54(3): 155-162. (in Chinese)

[8]ShuLC,LiangJ,HuQ,etal.Studyonsmallwind turbie icinganditsperformance[J].ColdRegions Scienceand Technology, 2017,134: 11-19.

[9]JinJY,VirkMS.Experimentalstudyof iceaccretiononS826amp;S832windturbineblade profiles[J].ColdRegions Science and Technology,2020,169:102913.

[10]IbrahimGM,PopeK,NatererGF.Scaling formulationofmultiphaseflowanddroplettrajectories withrimeicearetionona rotating wind turbine blade[J].JournalofWind Engineering and IndustrialAerodynamics,2023,232:105247.

[11]Wang Q,YiX,LiuY,etal.Simulationandanalysisof windturbineiceaccretionunderyawconditionviaanImproved MultiShot Icing Computational Model[J].Renewable Energy,2020,162: 1854-1873.

[12]Son C,Kim T.Developmentofanicingsimulationcode forrotatingwind turbines[J].Journalof Wind Engineeringand IndustrialAerodynamics,2020,203:104239.

[13]舒立春,任曉凱,胡琴,等.環境參數對小型風力發電機葉片覆冰特性及輸出功率的影響[J].中國電機工程學報,2016, 36(21): 5873-5878. Shu LC,RenXK,HuQ,etal.Influencesofenvironmentalparametersonicingcharacteristicsandoutputpowerofsmallwind turbine[J].中國電機工程學報,2016,36(21):5873-5878.(inChinese)

[14]舒立春,李瀚濤,胡琴,等.自然環境葉片覆冰程度對風力機功率損失的影響[J].中國電機工程學報,2018,38(18):5599- 5605. Shu LC,LiHT,Hu Q,etal.Efectsof ice degreeof blades on powerloses of wind turbines atnatural environments[J]. Proceedings of the CSEE,2018,38(18): 5599-5605.(in Chinese)

[15]舒立春,于周,李瀚濤,等.潔凈與覆冰條件下風力發電機輸出功率計算方法及現場試驗驗證[J].電工技術學報,2023, 38(11): 3041-3051. ShuLC,Yu Z,LiHT,etal.Calculation methodof wind turbine output powerunder cleanand icing conditions Verified by field tests[J].Transactions ofChina Electrotechnical Society,2023,38(11): 3041-3051. (inChinese)

[16]胡琴,楊大川,蔣興良,等.葉片模擬冰對風力發電機功率特性影響的試驗研究[J].電工技術學報,2020,35(22):4807- 4815. HuQ,YangDC,JiangXL,etal.Experimentalstudyontheeffectofbladesimulatedicingonpowercharacteristicsof wind turbine[J]. Transactions of China Electrotechnical Society,2020,35(22): 4807-4815.(in Chinese)

[17]李瀚濤,舒立春,胡琴,等.考慮覆冰粗糙度影響的風力發電機葉片氣動性能數值仿真[J].電工技術學報,2018,33(10): 2253-2260. LiHT,ShuLC,HuQ,etal.Numericalsimulationofwindturbine blades aerodynamic performance basedoniceroughne effect[J].Transactions ofChina Electrotechnical Society,2018,33(10):2253-2260.(in Chinese)

[18]韓興波,蔣興良,畢聰來,等.基于分散型旋轉圓導體的覆冰參數預測[J].電工技術學報,2019,34(5):1096-1105. HanXB,JiangXL,BiCL,etal.Predictionoficing environmentparametersbasedondecentralizedrotatingconductors[J] Transactions of China Electrotechnical Society,2019,34(5):1096-1105.(in Chinese)

[19]Makkonen L.Modeling of ice accretion on wires[J]. Journal ofApplied Meteorology,1984,23(6): 929-939.

[20]黎芷毓.風力發電機葉片霧淞覆冰數值模擬及自然覆冰驗證研究[D].重慶:重慶大學,2021. LiZY.Studyof numerical simulationof rime icingof wind turbine bladeand natural icing verification[D].Chongqing: ChongqingUniversity,2021.(in Chinese)

[21]梁健,舒立春,胡琴,等.風力機葉片雨淞覆冰的三維數值模擬及試驗研究[J].中國電機工程學報,2017,37(15): 4430- 4436,4584. Liang J,ShuLC,HuQ,etal.3-Dnumericalsimulationsandexperimentsonglaze iceaccretionof windturbineblades[J]. Proceedingsof the CSEE,2017,37(15): 4430-4436,4584.(in Chinese)

(編輯鄭潔)

主站蜘蛛池模板: 91精品情国产情侣高潮对白蜜| 国产日韩精品欧美一区灰| 亚洲无码免费黄色网址| 国产精品无码制服丝袜| 有专无码视频| 在线观看国产黄色| 国产剧情国内精品原创| 中国一级特黄视频| 亚洲 日韩 激情 无码 中出| 国产18在线播放| 成人精品视频一区二区在线| 国产精品污视频| 亚洲视频四区| 午夜毛片免费观看视频 | 91午夜福利在线观看精品| 看国产一级毛片| 国产欧美日韩另类精彩视频| 婷婷六月激情综合一区| 毛片免费网址| 国产精品hd在线播放| 国产人成乱码视频免费观看| 亚洲一级毛片在线播放| 日韩av在线直播| 专干老肥熟女视频网站| 免费观看欧美性一级| 欧洲高清无码在线| 久久伊人操| 欧美有码在线观看| 一级黄色片网| 欧美在线三级| 亚洲精品国产乱码不卡| 午夜三级在线| 天天爽免费视频| 国产成人91精品| 韩日免费小视频| 97人妻精品专区久久久久| 日韩在线永久免费播放| 欧美视频在线第一页| 国产精品女主播| 亚洲第一香蕉视频| 2021国产乱人伦在线播放| 在线观看无码av五月花| 国产精品第页| 国产特级毛片aaaaaaa高清| 亚洲天堂久久| 国产日韩欧美一区二区三区在线| 伊人中文网| 99热这里只有免费国产精品 | 日本欧美一二三区色视频| 亚洲综合第一区| 日本亚洲国产一区二区三区| 毛片a级毛片免费观看免下载| 91小视频在线| 亚洲国产一区在线观看| 色噜噜狠狠色综合网图区| 亚洲av无码片一区二区三区| 亚洲精品欧美重口| 热99精品视频| 91综合色区亚洲熟妇p| 成人一级黄色毛片| 国产第一页亚洲| 成人无码区免费视频网站蜜臀| 免费毛片全部不收费的| 中文字幕乱妇无码AV在线| 亚洲av无码人妻| 国产成人啪视频一区二区三区| 中文字幕av一区二区三区欲色| 亚洲大尺码专区影院| 欧美日韩免费观看| 欧美日韩亚洲国产主播第一区| 91青青在线视频| 毛片一区二区在线看| 国产美女91视频| 综合色婷婷| 一级成人a做片免费| 手机精品福利在线观看| 中文字幕日韩丝袜一区| 日韩高清无码免费| 国产小视频在线高清播放 | 无码一区18禁| 国产乱码精品一区二区三区中文 | 国产美女精品人人做人人爽|