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

噴水推進泵壓力脈動特性數值計算及分析

2014-08-07 12:17:01張明宇王永生靳栓寶魏應三付建
西安交通大學學報 2014年11期

張明宇,王永生,靳栓寶,魏應三,付建

(海軍工程大學動力工程學院, 430033, 武漢)

噴水推進泵壓力脈動特性數值計算及分析

張明宇,王永生,靳栓寶,魏應三,付建

(海軍工程大學動力工程學院, 430033, 武漢)

針對噴水推進器裝船后不均勻來流對壓力脈動特性的影響,以某巡邏艇噴水推進混流泵為研究對象,基于RANS方程和SST湍流模型,通過流體動力學軟件CFX穩態計算,進行了巡邏艇航速數值預報,所得計算值與試航值誤差為1.8%,從而驗證了計算流體動力數值計算的可信性。采用分離渦模擬方法,對敞水泵和裝船泵進行了三維非定常數值模擬,計算分析了葉輪進出口、葉輪內部、導葉內部及噴口5個截面和葉輪葉頂間隙處的壓力脈動,并對不均勻來流帶來的差別進行了研究。結果表明:在敞水泵和船后泵的葉輪出口、導葉內部,水流距葉輪越遠,壓力脈動影響越小,壓力脈動頻率取決于葉輪轉動頻率,壓力脈動幅值沿輪轂到輪緣逐漸增大,船后泵壓力脈動幅值整體高于敞水泵;對于均勻來流,敞水泵旋轉域葉輪室的壓力脈動頻率主要受導葉的影響,船后泵則受軸頻的影響,二者壓力脈動幅值在葉頂間隙處均從葉頂沿導邊到隨邊逐漸增大;對于敞水泵,流道出口壓力脈動頻率主要受葉頻控制,對于船后泵,壓力脈動頻率為軸頻。

噴水推進泵;計算流體動力學;壓力脈動;分離渦模擬

隨著艦船綜合性能的提高,噴水推進在高性能艦船、快艇上的應用越來越廣泛。作為噴水推進器核心部件的噴水推進泵(簡稱噴泵)的設計便顯得尤為重要,壓力脈動特性研究可在保證噴泵水力效率的情況下有效縮減過激的泵體振動及其引起的局部空化,更是降噪研究的基礎和前提。出于艦船的隱身性,研究噴泵內部非定常流場壓力脈動特性具有重要的意義。

轉動葉輪和靜止導葉之間的相對運動、葉片導邊的水流沖擊、葉片隨邊的脫流、局部空化等都會引起噴泵出現壓力脈動[1]。本文的研究對象噴泵為混流泵,以往泵內壓力脈動的研究通常是試驗研究,其周期長、成本高,對于所發現的非定常特性缺乏量化指標。隨著計算流體動力學(CFD)的逐步完善,國內外學者通過數值模擬的方法對泵內三維非定常流場進行了探索。Solis等通過改變葉輪形狀和徑向尺寸、采用SST湍流模型探索了減小壓力脈動對離心泵的影響[2]。Cheah等采用k-ε雙方程湍流模型對蝸殼內的脈動特性進行了模擬[3]。Zobeiri等對泵內轉子-定子的動靜干涉進行了數值模擬和分析[4]。Furukawa等基于葉表面奇點法對離心葉輪進行了無黏流場分析[5]。王春林等用大渦模擬研究了高比轉數混流泵非定常流場的壓力脈動特性[6]。施衛東等采用k-ε雙方程湍流模型研究了采樣頻率和時間對軸流泵壓力脈動特性的影響[7]。常書平基于RANS方程和SST湍流模型對噴水推進器內非定常流場進行了數值模擬[8]。靳栓寶對比分析了時間步長取葉輪旋轉3°、2°、1°所需時間對監測點壓力的影響[9]。

本文對某噴泵在敞水、船后兩種情況下進行了脈動特性分析,探索了船后不均勻來流對脈動的影響。本文設置了6組監測點,其中在葉輪葉頂間隙內設置了4個監測點,以分析葉頂間隙處壓力脈動特征,而在葉輪中間徑向設置了6個監測點,旨在為噴泵的噪聲控制提供有價值的參考。

1 計算模型和監控點分布

研究對象為比轉數等于445的混流泵,葉輪葉片數為4,導葉數為7,泵直徑為265mm,葉頂間隙為0.3 mm。為方便分析,將流道出口(泵進口)從輪緣到輪轂的監測點依次命名為D1,D2,…,D6,同理葉輪監測點依次為R1,R2,…,R6,葉輪出口監測點依次為RS1,RS2,…,RS6,導葉體監測點依次為S1,S2,…,S6,噴口監測點依次為N1,N2,…,N6,葉頂間隙從葉輪進口到出口依次為JX1,JX2,…,JX4,如圖1所示。

+:監測點

2 數值模擬

本文采用工程上應用廣泛的雷諾時均法進行CFD計算。首先選取SST湍流模型來封閉控制方程,對裝船泵內部流場進行穩態數值模擬。SST模型在近壁面區域調用k-ε模型模擬,收斂性好,而在湍流充分發展區域調用k-ω模型模擬,計算效率高[10]。在驗證數值計算的可信性之后,采用分離渦(DES)法[11]對敞水泵和船后泵進行三維非定常計算,求解其壓力脈動特性。

2.1 網格情況

采用結構化網格進行網格劃分:中葉輪網格數為106萬,導葉體網格數為102萬,由turbogrid劃分;噴口網格數為24萬,流道和計算域網格數為286萬,由ICEM中劃分;計算域的長、寬、高分別為泵直徑的30、10、8倍,網格總數為518萬。網格劃分情況如圖2所示。

(a)泵和噴口 (b)流道和計算域

網格邊界層均為0.1 mm,后處理Y+最大值為92,滿足計算要求。

2.2 域的設置和邊界條件

計算域包括噴口、導葉、葉輪和流道,其中僅有葉輪為旋轉域。設進口為速度進口,開口為壓力開口,相對靜壓值為0,出口為壓力出口,相對壓力值為0。除了3個交界面以外,其余部分均設置為壁面,葉輪葉片與輪轂為相對靜止壁面,其他為絕對靜止壁面。邊界條件如圖3所示。

圖3 邊界條件示意圖

2.3 初始條件及時間步長

泵的設計工況為2 400 r/min,穩態步長為30/(nπ)=0.004 s,n為轉速。瞬態計算以穩態結果為初始值。采樣頻率越高,所得數據的時間分辨率越高,采樣頻率過低,脈動幅值可能被低估[7]。結合以上情況,本文時間步長為每旋轉1°所需的時間,為1/(6n)=69.4 μs,共用時0.15s。

3 裝船泵穩態驗證

根據廠方提供船體阻力曲線,對不同轉速下的噴水推進器提供的凈推力進行了計算,并將凈推力曲線與船體阻力曲線疊加后通過兩曲線的交點來預報航速,結果如表1和圖4所示。

表1 2 400 r/min下噴水推進器凈推力

圖4 航速的數值模擬結果

由圖4可知,在交點處凈推力與船體阻力相等,此時泵工況為2 400 r/min,來流速度為78 km/h,泵功率為247.3 kW。由于廠商船體試航值為76.6 km/h,所以模擬誤差為1.8%,說明本研究所用數學模型和計算方法基本可信。

4 敞水泵與船后泵流道出口不均勻度

流道出口不均勻度是導致敞水泵和船后泵壓力脈動的主要原因。進水流道出口速度不均勻度[12]

(a)敞水泵 (b)船后泵

5 敞水泵壓力脈動特性計算與分析

暫不考慮裝船泵進水流道壓力脈動的影響,對敞水泵各截面進行了壓力脈動特性分析。對敞水泵進行三維非定常模擬時,以直管取代裝船流道,直管同樣采用結構化網格進行離散,進口處設置為流量進口,該流量值從船后泵在2 400 r/min、78 km/h工況中的流道出口提取,為634.9 kg/s。其他設置及監測點與船后泵完全相同,同樣是先計算穩態,再以穩態為初始條件計算瞬態。敞水泵模型及邊界設置如圖6所示。

圖6 敞水泵模型及邊界設置

5.1 敞水泵壓力脈動幅值分析

各截面監測點壓力脈動幅值如圖7所示。在敞水泵中來流均勻,葉輪內部、葉輪出口、導葉內部以葉輪葉片擾動為主,從輪緣到輪轂壓力脈動幅值逐漸減弱,壓力幅值從大到小依次為葉輪出口、流道出口、葉輪中、導葉內部、噴口。噴口處水流經導葉整流后距葉輪最遠,其受葉輪影響較小;葉頂間隙的4個監測點從葉頂沿導邊到隨邊的壓力幅值受葉片影響逐漸變大。

圖7 敞水泵監測點壓力脈動幅值直方圖

5.2 敞水泵監測點時域圖分析

為方便顯示,本文間隔性地選取一些點進行分析。如圖8a,由流道出口監測點1、2、4顯示,流道出口處監測一圈出現了4個明顯的小周期。這是由于來流均勻、脈動與葉頻相符、周期性受葉片影響所致。脈動幅值由輪緣到輪轂逐漸減小,說明靠近輪緣處受葉片影響較大。如圖8b,由于葉輪最外緣受葉片影響最大,所以通過對葉輪出口與導葉內部最外緣監測點進行分析發現,在一個周期內兩者也均有4個明顯的小周期,其幅值顯示,葉片旋轉頻率對壓力脈動的影響在導葉區不如葉輪區[6]。本文研究表明,在葉輪葉頂間隙處,從泵進口沿軸向,水流受葉輪影響逐漸增強,因此選擇JX4和葉輪最外緣監測點R1來顯示葉輪區域脈動情況。如圖8c,監測一圈在葉輪區域有7個明顯的小周期,與導葉數相同。這是因為葉輪為旋轉域,水流與導葉沖擊的反作用形成了葉輪域壓力脈動的時域特性,與物理現象相符。如圖8d,在噴口處,除監測點6例外,脈動壓力值從輪緣到輪轂呈下降趨勢,這是噴口中心的旋渦引起的,雖然距離葉輪較遠,但整體波峰差值偏小,也有明顯的4個周期。

(a)流道出口處監測點1、2、4

(b)導葉內部和葉輪出口最外緣處監測點

(c)葉頂間隙和葉輪內部監測點

(d)噴口處監測點

5.3 敞水泵壓力脈動頻域分析

根據三維瞬態輸出結果,采用傅里葉變換(FFT)對各監測點進行脈動頻譜分析[13]。各組監測點的壓力脈動頻域圖如圖9所示。為方便顯示,取各組監測點1、2、4、6,Fn為軸頻的倍數。葉輪葉片數為4,葉輪對水流的影響頻率及其諧波為軸頻的4倍(葉頻);導葉葉片數為7,導葉對水流的影響頻率及其諧波為軸頻的7倍。由圖9可知,在流道出口、葉輪出口、導葉內部,壓力脈動主頻率均為葉輪葉片頻率,即軸頻的4倍,且壓力脈動幅值從輪緣到輪轂依次減小,而在導葉中各監測點壓力脈動幅值波動范圍較小,下降趨勢減緩,這是導葉整流的結果。在葉輪旋轉域的葉頂間隙和葉輪中,每個旋轉周期內7片導葉對水流的反作用被切割成7份,壓力脈動主頻率為軸頻的7倍。葉頂間隙處從葉輪進口到出口沿軸向壓力脈動幅值增加,這是受葉片影響累積的結果,葉輪中沿徑向從輪緣到輪轂壓力脈動幅值逐漸減小。在噴口處最大壓力脈動頻域在1.4 kPa以下,由于導葉的消失、錐形輪轂引起的擾流等原因,所以規律性不明顯。

(a)流道出口

(b)葉頂間隙

(c)葉輪內部

(d)葉輪出口

(e)導葉內部

(f)噴口處

6 船后泵壓力脈動特性計算與分析

斜底船船后泵計算模型與邊界條件見圖3。以2 400 r/min、78 km/h工況進行瞬態計算,待穩定后,取最后3圈輸出結果平均值來分析脈動幅值走向,對各截面時域特性、頻域特性進行直觀展示,探究其規律及物理意義。

6.1 船后泵壓力脈動幅值分析

各監測點壓力脈動幅值如圖10所示,其中葉頂間隙只有4個監測點。

圖10 船后泵監測點脈動壓力幅值直方圖

由圖10可知,流道出口處因距離葉片導邊尚有一段距離,水流經過了流道彎管和軸之后發生紊亂,但依然呈現規律性,即在流道出口處,從輪緣到輪轂壓力脈動幅值逐漸增大,與敞水泵的情況相反。葉輪內部、葉輪出口、導葉內部以葉輪葉片擾動為主,從輪緣到輪轂脈動壓力幅值逐漸減小。整體幅值從大到小依次為葉輪出口、葉輪、導葉內部。噴口處距葉輪最遠,受葉輪影響較小,水流經導葉整流后的壓力幅值僅有微弱變化;葉頂間隙的4個監測點,從葉頂沿導邊到隨邊的累積壓力脈動幅值受葉片影響逐漸增大;船后泵整體壓力脈動幅值走向與敞水泵一致。

6.2 監測點壓力脈動時域特性分析

由于流道出口距離葉片導邊尚有一段距離,受葉輪影響微弱,流道彎管和軸引起的擾流掩蓋了葉輪對來流的反作用,所以在監測一圈后各監測點壓力脈動變化沒有規律性可尋,這與敞水泵均勻來流下呈現4個小周期存在明顯差異,且在同一監測點船后泵壓力脈動幅值明顯大于敞水泵。該結果同樣表現在旋轉域的葉輪中,由于葉輪設置為旋轉域,葉輪中和葉頂間隙處監測一圈為一個周期,一圈內流道和軸引起的擾流掩蓋了葉輪的影響,所以一圈監測點壓力脈動變化沒有規律性。葉輪出口和導葉內部最外緣監測點壓力脈動時域圖如圖11所示。

圖11 葉輪出口和導葉內部最外緣監測點壓力脈動時域圖

由圖11可見,在葉輪出口與導葉中,葉輪引起的擾流起主導作用,監測一圈出現了4個明顯的小周期,與葉頻相符,相比葉輪出口處監測點RS1,導葉中4個小周期的波峰相差較小,說明葉片旋轉頻率對壓力脈動的影響在導葉區不如葉輪區,該結果與敞水泵一致。噴口監測點壓力脈動時域圖如圖12所示。

圖12 噴口監測點壓力脈動時域圖

由圖12可見,噴口監測點壓力脈動呈現出4個小周期,但此時受葉片的影響很小,波峰差值很小,且壓力峰值從輪緣到輪轂依次遞減。結合以上葉輪出口和導葉分析可知,距葉輪越遠,受葉輪影響越小。

6.3 監測點壓力脈動頻域特性分析

為方便與敞水泵對比,仍取各組監測點1、2、4、6進行分析。各組監測點壓力脈動頻域特性如圖13所示。由圖13可見,與裸泵相比,監測一圈內,流道出口以流道彎管和軸引起的擾流為主,流道和軸對水流的影響頻率為軸頻,而敞水泵則因為來流均勻,脈動頻率主要為葉頻。由圖13b、13c可見,葉輪旋轉域頻率及其諧波為軸頻,這是不均勻來流掩蓋了導葉對水流反作用的結果,與敞水泵有明顯的不同。在葉輪出口、導葉內部和噴口,水流脈動主要受葉輪的影響,葉輪對水流的影響頻率及其諧波為葉頻,與敞水泵相似,只是因為受到不均勻來流影響。船后泵的壓力脈動幅值上大于敞水泵。

(a)流道出口

(b)葉頂間隙

(c)葉輪內部

(d)葉輪出口

(e)導葉內部

(f)噴口

7 結 論

采用分離渦方法對噴泵在敞水和船后的內部壓力脈動進行了數值計算和對比,得出以下結論。

(1)流道出口距離葉片導邊尚有一段距離,船后泵中水流經過流道彎管和軸之后產生紊亂,但依然具有規律性,即在流道出口,從輪緣到輪轂的壓力脈動幅值逐漸增大,該結果與敞水泵相反。

(2)葉輪出口和導葉內部的壓力脈動主頻以葉頻為主,壓力脈動幅值從輪緣到輪轂逐漸減小。

(3)旋轉域葉輪在均勻來流的敞水泵中的壓力脈動頻率以導葉對水流反作用頻率為主,即為7倍軸頻。壓力脈動幅值從輪緣到輪轂逐漸減小,減小幅度減緩;船后泵中,來流紊亂,壓力脈動頻率為軸頻。葉頂間隙,從葉頂沿導邊至隨邊累積壓力脈動幅值受葉片影響逐漸增大。

(4)噴口處,敞水泵和船后泵的壓力脈動幅值都很微弱,從殼體到軸心呈下降趨勢,監測點6突然增大,這是輪轂消失引起的旋渦所致。噴口處壓力脈動頻率均為葉頻。

(5)對比分析顯示,壓力脈動主要受葉輪的影響,從葉輪出口到噴口各監測點的壓力脈動幅值逐漸減小,說明隨著流體遠離葉輪,葉輪對壓力脈動的影響逐漸減小。

(6)相同的監測位置,船后泵脈動壓力幅值總體比敞水泵大。

[1] 施衛東, 鄒萍萍, 張德勝, 等.高比轉速斜流泵內部非定常壓力脈動特性 [J].農業工程學報, 2011, 27(4): 147-152.

SHI Weidong, ZOU Pingping, ZHANG Desheng, et al.Unsteady flow pressure fluctuation of high-specific-speed mixed-flow pump [J].Transactions of the CSAE, 2011, 27(4): 147-152.

[2] MOISES S, BAKIR S, KHELLADI F, et al.Pressure fluctuations reduction in centrifugal pumps: Influence of impeller geometry and radial gap [J/CD].ISRN Mechanical Engineering, 2011: ID 479594.

[3] CHEAH K W, LEE T S, WINOTO S H.Unsteady fluid flow study in a centrifugal pump by CFD method [C]∥4th International Symposium on Fluid Machinery and Fluid Mechanics.Beijing, China: Tsinghua University Press, 2009: 66-71.

[4] ZOBEIRI A, KUENY J L, FARHAT M, et al.Pump-turbine rotor-stator in generating mode pressure fluctuation in distributor channel [C]∥The 23th IAHR Symposium on Hydraulic Machinery and Systems.Yokohama, Japan: IAHR, 2006: 17-21.

[5] FURUKAWA A, TAKAHARA H.Pressure fluctuation in a vaned diffuser downstream from a centrifugal pump impeller [J].International Journal of Rotating Machinery, 2003, 9(4): 285-292.

[6] 王春林, 賈飛, 吳志旺, 等.高比轉數混流泵非定常流場壓力脈動特性 [J].排灌機械工程學報, 2013, 31(2): 103-108.

WANG Chunlin, JIA Fei, WU Zhiwang, et al.Pressure fluctuation of unsteady flow in high Specific speed mixed-flow pump [J].Journal of Drainage and Irrigation Machinery Engineering, 2013, 31(2): 103-108.

[7] 施衛東, 姚捷, 張德勝, 等.采樣頻率和時間對軸流泵壓力脈動特性的影響 [J].排灌機械工程學報, 2013, 31(3): 190-195.

SHI Weidong, YAO Jie, ZHANG Desheng, et al.Influence of sampling frequency and time on pressure fluctuation characteristics of axial-flow pump [J].Journal of Drainage and Irrigation Machinery Engineering, 2013, 31(3): 190-195.

[8] 常書平, 王永生, 魏應三, 等.噴水推進器內非定常壓力脈動特性 [J].江蘇大學學報: 自然科學版, 2012, 33(5): 522-527.

CHANG Shuping, WANG Yongsheng, WEI Yingsan, et al.Pressure fluctuation of unsteady flow in waterjet [J].Journal of Jiangsu University: Natural Science Edition, 2012, 33(5): 522-527.

[9] 靳栓寶, 王永生, 常書平, 等.混流泵內流場壓力脈動特性研究 [J].農業機械學報, 2013, 44(3): 64-68.

JIN Shuanbao, WANG Yongsheng, CHANG Shuping, et al.Pressure fluctuation of interior flow in mixed-flow pump [J].Transactions of the Chinese Society for Agricultural Machinery, 2013, 44(3): 64-68.

[10]靳栓寶, 王永生, 楊瓊方.基于數值試驗的噴水推進軸流泵的一體化設計 [J].中國造船, 2010, 51(3): 39-45.

JIN Shuanbao, WANG Yongsheng, YANG Qiongfang.Integration design of water jet axial flow pump based on numerical experimentation [J].Shipbuilding of China, 2010, 51(3): 39-45.

[11]劉學強, 伍貽兆, 程克用.用基于M-SST模型的DES數值模擬噴流流場 [J].力學學報, 2004, 36(4): 401-406.

LIU Xueqiang, WU Yizhao, CHENG Keyong.Computation of lateral turbulent jets using M-SST DES model [J].Acta Mechanica Sinica, 2004, 36(4): 401-406.

[12]魏應三, 王永生.噴水推進器進水流道不均勻度統一描述 [J].武漢理工大學學報, 2009, 31(8): 159-163.

WEI Yingsan, WANG Yongsheng.Unified description of out flow non-uniformity of waterjet duct [J].Journal of Wuhan University of Technology, 2009, 31(8): 159-163.

[13]黎義斌, 李仁年, 王秀勇, 等.低比轉數混流泵壓力脈動特性的數值模擬 [J].排灌機械工程學報, 2013, 31(3): 205-209.

LI Yibin, LI Rennian, WANG Xiuyong, et al.Numerical analysis of pressure fluctuation in low specific speed mixed-flow pump [J].Journal of Drainage and Irrigation Machinery Engineering, 2013, 31(3): 205-209.

(編輯 苗凌)

NumericalAnalysisforPressureFluctuationofWaterjetPump

ZHANG Mingyu,WANG Yongsheng,JIN Shuanbao,WEI Yingsan,FU Jian

(College of Marine Power Engineering, Naval University of Engineering, Wuhan 430033, China)

For ununiform inflow to impeller due to hull and duct, the dynamic characteristics of waterjet pump (mixed-flow pump) in patrol craft are numerically analyzed.The steady flow is simulated by SST turbulence model based on RANS and the calculated speed of ship from CFD exceeds the trial voyage by 1.8%.The unsteady flow of bare pump and the same pump fixed on craft is simulated by CFD method for DES turbulence.The pressure fluctuation of impeller inlet, impeller blade channels, impeller outlet, guide vane channels, nozzle and tip of impeller is supervised and discussed, and the difference caused by the ununiformity is also analyzed.The results show that on the impeller outlet, guide vane channels and the nozzle of both bare pump and fixed pump, the farther the flow to the impeller the weaker the pressure pulsation effects.The pressure fluctuation frequency is mainly dominated by the impeller rotation frequency, and the pressure fluctuation amplitude gradually increases from hub to tip, however, the corresponding amplitude is lager in pump fixed on craft.In the revolver impeller, the pressure fluctuation frequency is mainly dominated by the guide vanes reaction in bare pump while by shaft frequency in pump fixed on craft.The pressure fluctuation amplitude in the tip gradually increases along the axis from the leading edge to the trailing edge.In the duct outlet, the pressure fluctuation frequency is mainly dominated by the impeller rotation frequency in bare pump while by shaft frequency in pump fixed on craft.

waterjet pump; computational fluid dynamics; pressure fluctuation; detached-eddy simulation

2014-02-18。

張明宇(1989—),男,碩士生;王永生(通信作者),男,教授。

國家自然科學基金資助項目(51309229)。

時間:2014-09-01

10.7652/xjtuxb201411009

TH313

:A

:0253-987X(2014)11-0051-07

網絡出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20140901.1009.007.html

主站蜘蛛池模板: 一区二区在线视频免费观看| 在线精品亚洲国产| 91小视频版在线观看www| 欧美有码在线| 国产精品福利在线观看无码卡| 国产色爱av资源综合区| 成年午夜精品久久精品| 美女国内精品自产拍在线播放 | 国产成人你懂的在线观看| 日本高清免费一本在线观看 | 日韩美毛片| 91九色国产porny| 国产成人成人一区二区| 国产激情无码一区二区APP | 国产在线小视频| 精品视频一区二区三区在线播| 一级福利视频| a在线观看免费| 四虎综合网| 在线观看精品国产入口| 欧洲高清无码在线| 国产成人精品第一区二区| AV不卡国产在线观看| 国产精品一区在线麻豆| 久久99国产综合精品女同| yjizz视频最新网站在线| 性色生活片在线观看| 制服丝袜国产精品| 日韩在线网址| 欧美午夜一区| 亚洲毛片一级带毛片基地| 九九九九热精品视频| 中文字幕在线免费看| 99精品热视频这里只有精品7| 伊人婷婷色香五月综合缴缴情| 亚洲av综合网| 国产精品污污在线观看网站| 中文字幕 欧美日韩| 久久无码av三级| 欧美爱爱网| www.99精品视频在线播放| 国产新AV天堂| 特级做a爰片毛片免费69| 亚洲精品国产综合99久久夜夜嗨| 亚洲精品另类| 99尹人香蕉国产免费天天拍| 国产亚洲精品资源在线26u| 都市激情亚洲综合久久| 91福利片| 国产成人精品在线1区| 谁有在线观看日韩亚洲最新视频| 欧美 亚洲 日韩 国产| 97超级碰碰碰碰精品| 国产微拍精品| 欧美在线综合视频| 国产精品无码在线看| 国产欧美日韩精品综合在线| 国产成人精品三级| 日本不卡在线视频| 中国一级毛片免费观看| 欧美精品啪啪一区二区三区| 狠狠操夜夜爽| 亚洲欧洲日产国产无码AV| 久久视精品| 亚洲精品人成网线在线| 中文字幕亚洲综久久2021| 午夜福利无码一区二区| 亚洲人成在线精品| 亚洲无码高清视频在线观看| 青青草国产免费国产| 久久永久视频| 超碰aⅴ人人做人人爽欧美| 五月天丁香婷婷综合久久| 免费99精品国产自在现线| 久久久精品无码一区二区三区| 国产sm重味一区二区三区| 一级毛片在线播放免费观看| 国产极品嫩模在线观看91| 国产特一级毛片| 久久人人97超碰人人澡爱香蕉| 国产男女免费完整版视频| 97在线免费|