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

定子多風路空冷汽輪發電機流場模擬湍流模型驗證與分析

2018-05-14 13:31:17路義萍葛亞軍王浩然王芳
電機與控制學報 2018年5期
關鍵詞:風速模型

路義萍 葛亞軍 王浩然 王芳

摘 要:研究一種湍流模型模擬定子多風路大型空冷汽輪發電機中的具有多種流動形勢轉變特征的復雜流動。以一臺三進四出定子多風路的空冷汽輪發電機為原型,建立了包括風扇、定子、轉子、氣隙在內的通風系統實驗臺,采用熱線風速儀測量了1 000 r/min穩態工況下出風區中多個定子風溝出口截面的風速;然后,建立了實驗臺物理模型并采用三維建模軟件進行了數學建模,在相同工況下進行數值計算,同時將測量結果與采用3種湍流模型得到的計算結果進行對比。結果表明:采用RNG k-ε模型時的風速計算結果與實驗結果誤差最小,RNG k-ε模型更適合模擬定子多風路的大型空冷汽輪發電機整機內部的流場。最后,分析了采用RNG k-ε湍流輸運模型時實驗臺內部的流場分布特征。所得結論可為大型空冷汽輪發電機通風冷卻設計提供理論參考。

關鍵詞:定子多風路空冷汽輪發電機;實驗研究;數值模擬;流場;湍流模型選取

中圖分類號:TM 311, TK 121

文獻標志碼:A

文章編號:1007-449X(2018)05-0063-07

Abstract:A suitable turbulence model was researched for the simulation of the complex flow characteristics in multiple ventilation large air cooled turbo generator. A three inlets and four outlets of multiple ventilation stator ducts of turbogenerator with air coolant was comprised as the prototype, including fan, stator, rotor, air gap and ventilation system experimental bench. The values of wind speed of multiple stators duct outlet zone were measured under 1 000 r/min steadystate condition using the hotwire anemometer. Then, physical model was established by using threedimensional establishing software and was used for numerical calculation under the same condition. The measurement results are compared with the calculation results of standard k-ε, realizable k-ε and the RNG k-ε turbulence models. The results show that calculation of RNG k-ε model has the minimum error with experimental results, which is more suitable for the simulation of the distribution of the flow field in a large air cooled turbo generator. Finally, experimental bench interior flow field distribution characteristics of RNG k-ε model adopted were analyzed. The conclusions provides a theoretical reference for the design of ventilation cooling of Large Aircooled Turbo generator.

Keywords:multiple ventilation stator ducts of aircooled turbo generators; experiment; numerical simulation;flow field;selection of turbulence model

0 引 言

目前,火力發電中利用熱效率較高的燃氣輪機的復合型發電技術正在全球不斷擴大應用,要求與之配套的空冷汽輪發電機也達到高功率、大容量化。由于汽輪發電機冷卻系統以空冷系統的設備最為簡單,易于維護,安全性好等優點,近年來對空冷電機的需求增加迅速,機組的單機容量越來越大,電磁負荷增加顯著,熱負荷較高,繞組的溫升較大,當溫度升高到一定程度將危及到絕緣材料及發電機的安全運行,因而對汽輪發電機的冷卻研究越來越重要[1]。

早期,文獻[2]通過改變副槽出風孔的直徑及采用傾斜副槽風道,利用電風速計測量靜止的徑向風溝出風孔的風速,研究副槽的通風均勻性;文獻[3]在總結靜態實驗及電機冷卻事故的基礎上,建立了旋轉動態模型;文獻[4]對東方電機有限公司研發的220 MW大型空冷汽輪發電機的結構特點、通風及運行測試結果進行了系統全面的研究;文獻[5]對220 MW定子多風路大型空冷汽輪發電機的定子與氣隙,采用計算流體動力學(CFD)方法,在電機總風量不變的情況下,研究了改變定子、氣隙風量分配百分比對定子各個徑向通風溝內的流體速度、溫度變化的影響;文獻[6]以包括定轉子及端部的一體化物理模型為研究對象,采用與文獻[5]相同的湍流Standard k-ε兩方程等模型,反演多風路電機中的流場特征;文獻[7]研究了3種湍流k-ε兩方程模型對空冷汽輪發電機中轉子熱流場的影響;文獻[8-9]采用CFD方法分別研究了大型空冷汽輪發電機定子部分、端部的流動及傳熱特性;文獻[10]采用實驗研究方法研究了定子向氣隙射流對定轉子間環形空間中流動的影響;文獻[11]建立了轉速相對較低的風力發電機中圓盤形的定轉子系統的通風實驗臺,研究冷卻空氣經定子射流并沖擊轉子表面流動特征及表面局部傳熱性能;文獻[12]介紹了定子多風路電機實驗臺研發考慮的相關因素及強度計算。綜上可見,國內外針對空冷汽輪發電機的冷卻研究多數采用數值模擬與實驗研究相結合的方法。與此同時,其他類型電機也采用同類方法研究[13-15]。

空冷汽輪發電機在冷卻過程中,其內部流場處于湍流狀態,任一湍流模型都不是普遍適用的,不同通風特征流場應該選擇的湍流模型應該有所不同。現有的湍流模型很多,針對水力旋轉機械及通常的空冷汽輪發電機轉子,一些學者也提出了湍流模型應用的選擇[16-17]。對于定子三進四出多風路的大型空冷汽輪發電機,氣隙內部存在正、逆向流動,旋轉射流,風扇與轉子內流體高速旋轉流,存在多種流動形態轉變的復雜流場,尚未見轉子旋轉狀態下流場實驗與相應的多種湍流模型計算比較的文獻。鑒于此,本文以一臺定子多風路的空冷汽輪發電機為原型,與生產廠商一起研發設計并制造了包括風扇、定子、轉子、氣隙在內的通風系統實驗臺。首先,基于ANSYS Fluent軟件平臺,在網格獨立條件下,采用Standard k-ε、Realizable k-ε和RNG k-ε 3種湍流模型時數值計算電機實驗臺內部在1 000 r/min工況下的流場分布。然后,采用熱線風速儀測量電機,定子各風溝出口處的風速大小,并與3種湍流模型的計算結果進行對比,分析正確湍流模型下的流場特征,旨在探討不同湍流模型在定子多風路大型空冷汽輪發電機內部復雜流場模擬中的適用性。

1 物理模型的建立

1.1 實驗臺的搭建及通風系統

本文的電機通風實驗臺是以某大型空冷汽輪發電機風路為原型,經流動的相似性分析、強度計算、簡化等設計過程,最終制造而成。考慮到空冷汽輪發電機為兩側中心對稱供風結構,實驗臺僅考慮本體部分通風,不考慮端部通風及整機熱問題,采用整機的一半風路,即只在電機實驗臺的一端設置風扇,另一端為整機的中心對稱面且采用封閉結構,具體結構包括風扇、定子軸徑向環板與配風筒、定子塊隔成的風溝、轉軸上的筋板隔成的8個水平槽(替代副槽)與轉子環上的徑向風孔組成的轉子風路、氣隙等,如圖1所示。其中:擋風罩前端面為實驗臺的空氣入口,轉軸上的風扇同時為8個定子配風筒、8個轉子水平槽、定轉子間氣隙供風;沿周向布置的8根定子配風筒及環板使定子通風成為整機三進四出、本實驗臺中為兩進兩出風路,進風區頂部采用鐵皮封閉,出風區為敞開形式,見圖2,出口所在機殼頂面設置為實驗臺空氣出口。搭建完成的實驗臺三維剖面結構和實體模型如圖1、圖2所示。

電機通風實驗臺中,空氣由定子機座前端擋風罩所圍成的進口處進入,經軸流風扇吸入后共分為3個風路,一路沿著軸向直接進入實驗臺定子和轉子之間的氣隙;一路沿軸向進入轉子環端部,再分別進入8個轉子水平槽道,經轉子環上的徑向出風口進入氣隙;一路沿軸向分別經4個長和4個短的定子配風筒進入兩定子進風區,此外,經定子鐵心徑向風溝進入氣隙;3路空氣在氣隙內混合后經過Ⅰ與Ⅱ兩個定子出風區的出口面流出,見圖1。

2 數值計算

2.1 基本假設

針對于電機通風實驗臺內空氣流動,涉及到流動狀態等物理參數,為了選擇正確的控制方程來對流場進行描述,對計算域做出以下基本假設:1)實驗臺內部空氣流場參數不隨時間變化,僅研究穩態流動;2)忽略重力對實驗臺內部空氣流動的影響;3)電機內流體流動的雷諾數很大(Re>2 300),屬于湍流,采用湍流模型對電機內流場進行求解;4)電機內流體流速遠小于聲速,即馬赫數很小,故將流體作為不可壓縮流體處理。

2.2 數學模型

實驗臺在穩定運行時,轉子和風扇相對定子進行高速轉動,轉子內的空氣跟隨主軸旋轉,風扇周圍的空氣跟隨扇葉做周向和軸向的運動。在計算流體動力學分析中采用多重參考系,認為轉子內部及風扇周圍的空氣相對其他部分做固定轉速旋轉,采用多重參考系進行計算。在三維湍流流場求解過程中采用質量、動量守恒方程,其通用控制方程[18]為

2.3 邊界條件與求解

由于實驗臺進出口直接和周圍環境空氣相通,所以入出口均為壓力邊界條件,表壓力均為0 Pa;在轉子內表面及風扇扇葉和輪轂的表面這些與旋轉固體相接觸的表面設置旋轉壁面邊界條件。近壁面處采用標準壁面函數,近壁面Y+大于30。計算域內部網格節點的離散方程組采用隱式、分離求解,其中速度與壓力耦合方程采用SIMPLE算法求解,殘差取1103。經多次改變網格數量和質量,獲得網格獨立收斂解。

3 實驗研究及結果分析

3.1 實驗測量

本文的實驗測量采用熱線風速儀測量電機通風實驗臺在1 000 r/min穩定工況下定子風溝出口所在表面的空氣速度,新購置熱線風速儀不需標定。考慮實驗測量的方便性以及被工字鋼條隔開的每個定子風溝出風口截面圓周方向的對稱性,本文選取以下位置進行測量:軸向測量位置如圖2所示,從非傳動端看,沿軸向對各測量的定子出風區Ⅰ與Ⅱ中的風溝進行標號,測量1-11號定子風溝表面風速,見圖1、2中標號;周向風溝出口截面測量位置如圖3所示,從非傳動端看的左側定子風溝周向測量1-9號位置的風速,周向和軸向兩個標號確定出風溝的出口面的唯一位置。然后,在每一個風溝截面上沿旋轉方向,分別在風溝圓弧形截面沿周向中間弧線上先標定測量1/5、1/2與4/5位置共3點位置,最后測量風速。

測量發現,周向出風截面上的風速是不均勻的,風速大小與旋轉方向有關,被工字鋼條隔開的每個小風溝截面上,迎風面風速大于背風面;存在軸向拉桿的位置,如圖3中的周向標號為3的風溝出風面,風速總體偏低。實驗時,每個測點測量多次,直到讀數穩定為止。

3.2 實驗結果的不確定度分析

對于一個有價值的測量結果必須進行實驗不確定度評價。不確定度的范圍反映了實驗數據與真值之間的靠近程度。不確定度愈小,實驗數據與真實值越靠近,測量的質量越高。本文實驗所測的數據的不確定度為直接測量的不確定度。

本文中熱線風速儀儀器誤差為+3%。采用文獻[19]中的方法算得風速測量結果的不確定度為+8.3%,即風速測量結果將比真實值大8.3%。

3.3 實驗結果與數值計算結果的比較

將實驗方法測得的軸向1-11號定子出風區風溝出口各個周向小弧形截面所有點測得的風速分別求得關于軸向每條風溝完整出風面的幾何平均風速值,作為軸向1-11號定子出風區風溝出口截面上平均風速,然后,再利用Fluent軟件采用前述3種湍流模型計算出相同試驗條件下模型試驗臺中的定子兩出風區1-8及9-11號風溝出口截面的風速平均值,兩種結果比較見圖4。

本文實驗臺中的流場是風扇與轉子組成的串聯旋轉坐標系下的流動疊加固定坐標系下的定子與氣隙多風路復雜流動,進行實驗并數值模擬研究是非常必要的。從圖4的對比中可以看出,除了兩端位置,3種湍流模型計算結果僅RNG k-ε模型計算結果與實驗測量結果總體趨勢基本一致,測得的風道靠近風扇端11號風溝風速最大,預測值與實驗測量值趨勢不一致,實驗結果表明,靠近風扇端的Ⅱ風區,沿流動方向11、10與9號風溝風量逐漸減小,采用CFD數值模擬時,風量先增大,然后減小,邊界上11號風溝計算值偏低,在Ⅰ風區中,8條風溝中的風速曲線呈下凹的拋物線形式,中部的3~4號風溝風速較大,數值計算得到的1號風溝風速計算值偏高;一般而言,CFD數值計算時邊界位置計算準確性較差;實驗測量風速總是大于采用各種湍流模型計算值,測量儀器說明書中介紹儀器本身存在+3%的系統誤差,測量不確定度為+8.3%;在兩風區中,RNG k-ε模型計算結果和實驗都更接近,Realizable k-ε模型計算結果的風速最小;與實驗結果相比,在Ⅱ風區中,刨除計算起始邊界11號風溝,3種湍流模型的計算誤差均較小,而對于風速變化曲線下凹時,Standard k-ε、Realizable k-ε兩模型的計算出的風速極值及位置誤差大,中間風溝風速本應很大,模擬出的數值偏低,導致曲線形式不正確,RNG k-ε模型準確預測出了極值位置,風速分布兩邊低中間高這種下凹的曲線形式是正確的,總體上RNG k-ε湍流模型計算結果誤差最大處為14.7%,誤差原因包括假擴散誤差、離散誤差、舍入誤差等。此外,如考慮實驗結果的不確定度為+8.3%,CFD數值模擬風速也是比風速真實數值低,原因一方面是湍流模型本身是對湍流客觀規律的近似表達,另一方面,定子多風路大型電機風路本身涉及多種流動形式的急劇轉變,特別是在定子風溝兩端邊界位置、氣隙中位置。

因此,采用CFD軟件中的兩方程模型研究定子多風路大型電機內部溫度場時,一般而言,數值計算出來的定子峰值溫度相應的比真實的測量溫度要高,依此峰值溫度數據作為依據,判斷定子額定工況運行是否超溫,能更安全些。通過上述實驗研究,對于定子多風路的大型空冷汽輪發電機整機三維湍流流場數值模擬,推薦采用RNG k-ε湍流輸運模型。

4 數值計算結果及分析

4.1 壓力分析

為了分析比較不同湍流模型計算結果方便,圖5給出了采用Standard k-ε與RNG k-ε湍流模型時,實驗臺沿周向極角11.25°截面的空氣壓力分布圖。從圖中可以看出2種湍流模型的計算結果的壓力場分布趨勢相同,壓力數值大小有差別。相同點是壓力最大值均位于實驗臺轉子風扇扇葉頂端與擋風罩之間的間隙位置,最小值位于風扇入口吸力面處;對于風扇和轉子等旋轉部件,由于旋轉作用,在旋轉中心形成了負壓,在轉子外邊緣及以上定子區域均為正壓,因而沿徑向方向壓力呈階梯分布;3種湍流模型計算結果的最高壓力值大致相同,在570 Pa 左右。不同點是RNG k-ε模型計算結果的最低值是-711 Pa,比其他兩種模型要低310 Pa左右,這主要是由于RNG k-ε模型能夠更好地模擬這種風扇連同轉子的旋轉流動效應,在相同網格時形成了更小的負壓。由于RNG k-ε模型計算結果壓差大,進入實驗臺的總空氣量最多。

另外,定子配風筒相通的Ⅱ、Ⅰ定子進風區風溝背部壓力沿空氣流動方向逐漸降低,直到出風區出口表壓力為0;靠近風扇的定子兩個Ⅱ風區下方的氣隙位置壓力與定子Ⅰ進風區下方氣隙壓力接近,RNG k-ε模型預測出的氣隙中靜壓略高于其他模型。

4.2 速度分析

圖6僅給出了實驗臺中采用RNG k-ε湍流模型時,沿軸向極角11.25°截面的速度分布云圖,其他兩模型計算結果的分布規律與之相同,僅數值大小有差異。最大速度都在風扇扇葉頂端處,最小速度都主要分布在入出口處。水平槽道內的空氣流速隨旋轉半徑增大逐漸增加,靠近轉軸的區域速度低,靠近轉子環的區域速度高;氣隙內部的前方空氣受到風扇提供的壓力作用,克服流動阻力,部分從Ⅱ出風區流出,其他與定轉子徑向流進氣隙的空氣相遇,并分別由出風區流出;沿軸向靠近傳動端風扇處氣隙截面內的空氣流速大于非傳動端截面內的流速。風溝出口處截面平均風速分析如前述,此處不贅述。

5 結 論

本文建立了多風路的通風系統旋轉流實驗臺,在網格無關性的條件下,比較了Standard k-ε模型、Realizable k-ε模型和RNG k-ε模型等3種湍流模型在1 000 r/min的工況下實驗臺內部的流場分布特征并與實驗結果進行了對比,得到以下結論:

1)3種湍流模型僅RNG k-ε模型實驗誤差在14.7%以內,且風溝出風截面平均風速計算結果曲線與測量得到的風速分布曲線趨勢一致,滿足工程計算要求;

2)RNG k-ε模型的計算結果與實驗結果誤差最小,更適合模擬定子多風路的大型空冷汽輪發電機內部的流場分布。

參 考 文 獻:

[1] 金煦, 袁益超, 劉聿拯. 大型空冷汽輪發電機冷卻技術的現狀與分析[J]. 大電機技術, 2004, (4): 33.

JIN Xu, YUAN Yichao, LIU Yuzheng. State and analysis of large aircooled turbogenerator[J]. Large Electric Machine and Hydraulic Turbine, 2004, (4):33.

[2] 李廣德, 張偉紅. 空冷汽輪發電機的通風系統設計[J] . 大電機技術, 1998, (4): 12.

LI Guangde, ZHANG Weihong. Design of aircooling turbo generator ventilation system[J].Large Electric Machine and Hydraulic Turbine,1998(4): 12.

[3] 廖毅剛, 侯小全. 全空冷汽輪發電機通風冷卻研究[J]. 東方電氣評論, 2008(1): 1.

LIAO Yigang, HOU Xiaoquan. Study on ventilation system of full aircooling steamturbogenerator[J]. Dongfang Electric Review, 2008(1):1.

[4] 令紅兵, 胡德劍. 東方大型空冷汽輪發電機的設計開發[J]. 東方電氣評論, 2009, 23(92): 29.

LING Hongbing, HU Dejian. Large capacity aircooled turbo generator developed by DEC[J]. Dongfang Electric Review. 2009, 23(92): 29.

[5] 李偉力, 楊雪峰, 顧德寶. 空冷汽輪發電機冷卻氣流風量對定子內流體的影響[J]. 中國電機工程學報,2009, 29(21): 53.

LI Weili, YANG Xuefeng, GU Debao. Influence of air current flow change on fluid flow and heat transfer of aircooled turbogenerator with multipathventilation[J]. Proceedings of the CSEE, 2009, 29(21):53.

[6] 路義萍, 洪光宇, 湯璐, 等. 多風路大型空冷汽輪發電機三維流場計算. 中國電機工程學報, 2013, 33(3):133.

LU Yiping, HONG Guangyu, TANG Lu,et al.Calculation of 3D flow field of large aircooled turbogenerators with multipath ventilation[J].Proceedings of the CSEE, 2013, 33(3): 133.

[7] 路義萍, 潘慶輝, 孫雪梅. 湍流模型變化對汽輪發電機轉子熱流場影響[J]. 電機與控制學報, 2014, 21(1): 85.

LU Yiping, PAN Qinghui, SUN Xuemei, et al. Effect of turbulence models on temperature and flow field for turbogenerator rotor[J]. Electric Machines and Control, 2014, 21(1):85.

[8] 王芳, 郭瑞倩, 安志華, 等. 空冷發電機定子三維溫度場分布與試驗對比[J]. 電機與控制學報, 2013, 17(12): 47.

WANG Fang, GUO Ruiqian, AN Zhihua, et al.Air cooled generator stator temperature field distribution and experimental comparison[J]. Electric Machines and Control, 2013, 17(12): 47.

[9] HAN Jichao, LI Weili, WANG Likun, et al.Calculation and analysis of the surface heattransfercoefficient and temperature fields on the threedimensional complex end windings of a large turbo generator[J]. IEEE Transactions on Industrial Electronics, 2014, 61(10):5222.

[10] Jeng, T.M. , Tzeng, S.C., Chang, H.Flow visualization in an annulus between coaxis rotatingcylinders with a circular jet on stationary outer cylinder[J].International Communications in Heat and Mass Transfer,2012, 39 (8): 1119.

[11] PELLE J, HARMAND S. Convective heat transfer in a rotor stator system air gap with a centered natural suction of fluid[J]. Journal of Heat Transfer, 2011, 133(11):147.

[12] 路義萍, 潘慶輝. 多風路空冷汽輪發電機PIV流場測試實驗臺研發[J]. 哈爾濱理工大學學報, 2013, 19(4): 32.

LU Yiping, PAN Qinghui. Development of experimental setupof PIV flow field test for multiple ventilation ducts of turbo generator with air coolant[J]. Journal of Harbin University of Science and Technology, 2013, 19(4): 32.

[13] 韓家德, 杜鵬, 路義萍. 凸極電機定子風路變化對熱流場影響[J]. 電機與控制學報, 2016, 20(12): 59.

HAN Jiade, DU Peng, LU Yiping. Effect of statorventilation ducts changes on thermal and flow field of salient synchronous motor[J]. Electric Machines and Control, 2016,20(12): 59.

[14] 溫嘉斌, 侯健, 于喜偉. 定子通風槽鋼對中型高壓電機內溫度場的影響[J]. 電機與控制學報, 2016, 20(8): 40.

WEN Jiabin, HOU Jian, YU Xiwei. Influence of stator ventilation channel on the temperature field in the middlesize high voltage motor[J]. Electric Machines and Control, 2016, 20(8): 40.

[15] 路義萍, 張東學, 孫博,等. 某無刷勵磁機通風冷卻數值模擬研究[J]. 電機與控制學報, 2016, 20(6): 26.

LU Yiping, ZHANG Dongxue, SUN Bo, et al. Numerical simulation of ventilation cooling system for brushless exciter[J]. Electric Machines and Control, 2016, 20(6): 26.

[16] 張宏飛, 曹紅松, 趙捍東, 等.數值仿真中湍流模型的選擇[J] . 彈箭與制導學報, 2006, 26(4): 242.

ZHANG Hongfei, CAO Hongsong, ZHAO Handong, et al. The Choice of turbulence model innumerical simulation[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2006, 26(4): 242.

[17] 黃彪, 王國玉. 基于k-ω SST模型的DES方法在空化流動計算中的應用[J]. 中國機械工程, 2010, 21(1): 85.

HUANG Biao, WANG Guoyu. Application of a SSTDES turbulence model for computation of cavitating flows[J]. China Mechanical Engineering,2010,21(1): 85.

[18] Fluent Inc.(2012). Fluent 14.5 User′s Manual. Fluent Inc, USA.

[19] Kline, S. J.,and F.A.McClintock: Describing uncertainties in singlesample experiments, mech. Eng.,p.3, January 1953.

(編輯:張 楠)

猜你喜歡
風速模型
一半模型
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
基于GARCH的短時風速預測方法
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
GE在中國發布2.3-116低風速智能風機
主站蜘蛛池模板: 国产精品女熟高潮视频| 亚洲伦理一区二区| 精品三级网站| 一区二区偷拍美女撒尿视频| 色视频久久| 蜜臀AVWWW国产天堂| 亚洲伦理一区二区| 香蕉在线视频网站| 久久精品一品道久久精品| 福利小视频在线播放| 亚洲AⅤ综合在线欧美一区| 在线观看网站国产| 亚洲精品综合一二三区在线| 免费一级大毛片a一观看不卡| 日韩毛片免费视频| 国产精品理论片| 国产一级特黄aa级特黄裸毛片| 国内精自视频品线一二区| 国产裸舞福利在线视频合集| 久久午夜夜伦鲁鲁片不卡| 久久国产精品波多野结衣| 色国产视频| 99re视频在线| 欧洲日本亚洲中文字幕| 黄色片中文字幕| 99视频在线免费| 亚洲永久视频| 精品视频免费在线| 精品久久久久久成人AV| 98精品全国免费观看视频| 嫩草影院在线观看精品视频| 精品无码一区二区三区电影| 国产精欧美一区二区三区| 国产小视频在线高清播放| 亚洲黄色视频在线观看一区| 极品国产一区二区三区| 任我操在线视频| 99热这里只有精品国产99| 国产欧美在线观看精品一区污| 亚洲国产日韩欧美在线| 91成人免费观看| 亚洲视频色图| 国产美女无遮挡免费视频| 呦女精品网站| 亚洲美女一级毛片| 54pao国产成人免费视频| 亚洲人成网站18禁动漫无码| 伊人久久精品亚洲午夜| 91久久偷偷做嫩草影院免费看| 欧美色亚洲| 激情亚洲天堂| 国产国产人免费视频成18| 亚洲日韩高清在线亚洲专区| 国产av剧情无码精品色午夜| 国产精品男人的天堂| 精品福利国产| h网站在线播放| 亚洲视频一区在线| 色网站在线视频| 精品综合久久久久久97| 国产一级妓女av网站| av在线手机播放| 特级做a爰片毛片免费69| 中文字幕乱码二三区免费| 国产毛片高清一级国语| 日本午夜网站| 一本无码在线观看| 大学生久久香蕉国产线观看 | a级毛片一区二区免费视频| 亚洲精品成人片在线观看| 久久国产黑丝袜视频| 国产精品lululu在线观看| 色老二精品视频在线观看| 一级一毛片a级毛片| 亚洲精品成人福利在线电影| 亚洲午夜福利精品无码| 精品国产成人高清在线| 国产一级二级在线观看| 毛片免费视频| 中文字幕无码中文字幕有码在线 | 午夜人性色福利无码视频在线观看| 91黄色在线观看|