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

基于流固耦合的離心泵啟動過程瞬態(tài)葉片動應(yīng)力特性

2016-08-04 06:16:38袁建平夏水晶宗偉偉周幫倫付燕霞
振動與沖擊 2016年12期
關(guān)鍵詞:效應(yīng)變形

袁建平, 夏水晶, 宗偉偉, 周幫倫, 付燕霞

(江蘇大學 流體機械工程技術(shù)研究中心,江蘇 鎮(zhèn)江 212013)

基于流固耦合的離心泵啟動過程瞬態(tài)葉片動應(yīng)力特性

袁建平, 夏水晶, 宗偉偉, 周幫倫, 付燕霞

(江蘇大學 流體機械工程技術(shù)研究中心,江蘇 鎮(zhèn)江212013)

為研究離心泵啟動過程的葉片的應(yīng)力和變形情況,針對IS65-50-160的離心泵啟動過程瞬態(tài)內(nèi)部流場和結(jié)構(gòu)場進行了雙向流固耦合聯(lián)合求解。其中流場計算基于RANS方程與SST湍流模型;結(jié)構(gòu)場計算基于彈性體結(jié)構(gòu)動力學方程;對葉輪葉片在雙向流固耦合作用下的變形和應(yīng)力分布進行了計算,獲得了離心泵啟動過程中瞬時效應(yīng)對葉片應(yīng)力和應(yīng)變的影響規(guī)律。研究結(jié)果表明:葉片的最大等效應(yīng)力和應(yīng)變呈振蕩上升趨勢,振動強度先減小后增大; 在離心泵啟動過程中,葉片進口與后蓋板的交界處出現(xiàn)應(yīng)力集中,葉片的變形量從后蓋板到前蓋板呈遞增的趨勢;葉片的最大等效應(yīng)力和應(yīng)變量都大于穩(wěn)態(tài)工況下的最大等效應(yīng)力和應(yīng)變量。

離心泵;流固耦合;應(yīng)力分布;應(yīng)變;數(shù)值模擬

離心泵的啟動過程內(nèi)部瞬態(tài)流動非常復(fù)雜,葉輪轉(zhuǎn)速由零迅速上升為額定值,特征雷諾數(shù)從零迅速上升至上百萬,內(nèi)部流動從層流狀態(tài)迅速變成復(fù)雜的湍流流動,泵內(nèi)部流場十分不穩(wěn)定,容易出現(xiàn)二次流、動靜干涉等物理現(xiàn)象,使泵內(nèi)產(chǎn)生強烈的振動和噪聲,甚至損壞泵機組設(shè)備[1-6]。目前關(guān)于離心泵啟動階段內(nèi)部流場瞬態(tài)特性還缺乏系統(tǒng)的研究,人們一般采用準穩(wěn)態(tài)假設(shè)方法代替瞬態(tài)過程。

流固耦合不僅考慮了流體載荷對固體變形的作用,而且還能實現(xiàn)固體變形對流動結(jié)構(gòu)的影響。近年來,眾多學者對流固耦合應(yīng)用進行了系統(tǒng)的研究,其中Benra等[7]分別采用流固耦合的兩種分析方法對無堵塞單葉片離心泵轉(zhuǎn)子部件進行了研究。Langthjem等[8]對二維離心泵進行了數(shù)值計算,得出了葉輪與流體間的相互作用是引起離心泵發(fā)出噪聲的重要因素。Guadaqn等[9]對一新型泵進行了流固耦合計算,通過試驗測試驗證了計算結(jié)果的準確性。王洋等[10]對離心泵進行了流固耦合計算,對不同工況下離心泵葉輪最大等效應(yīng)力和最大總變形量進行了對比分析。

總之,目前國內(nèi)外學者關(guān)于離心泵的研究主要集中在穩(wěn)態(tài)內(nèi)流特性等方面,而關(guān)于離心泵啟動瞬態(tài)特性的研究還不夠完善。流固耦合方法主要應(yīng)用于離心泵流動誘導(dǎo)振動的分析,在離心泵瞬態(tài)特性研究方面很少。

本文采用雙向流固耦合方法對離心泵啟動過程內(nèi)部流場和結(jié)構(gòu)場進行聯(lián)合求解,獲得了離心泵啟動階段葉片最大等效應(yīng)力和最大變形量隨時間的變化規(guī)律,并對葉輪進行了強度校核,為確保離心泵安全可靠的運行提供理論依據(jù)。

1實驗

為了保證離心泵瞬態(tài)計算的準確性,同時確定模擬的邊界條件,本文對模型泵進行了試驗研究。試驗裝置見圖1。離心泵啟動瞬態(tài)性能試驗裝置是一個閉式系統(tǒng),試驗臺主要由動力驅(qū)動裝置、試驗泵裝置和管路系統(tǒng)三部分構(gòu)成。

圖1 水泵試驗裝置圖Fig.1 The test rig of pumps

測試過程前,先對閥門開度進行調(diào)節(jié),盡量使得在管路特性完全相同的情況下進行對比試驗,并確保離心泵最終在設(shè)計工況下運行。采用變頻器設(shè)置離心泵的啟動時間為1 s,運行的穩(wěn)定轉(zhuǎn)速為2 900 r/min,閥門的開度調(diào)試到額定流量對應(yīng)的開度。

在1 s啟動時間過程中分別記錄了轉(zhuǎn)速,流量和進出口壓力隨時間的變化曲線,瞬態(tài)結(jié)果見圖2。從圖2可知,流量滯后于轉(zhuǎn)速和揚程到達最大值,隨著啟動時間的增長,流量到達最大流量的滯后時間也逐漸加長。從流量時域曲線可以看出,啟動過程中流量曲線變化趨勢近似于三次曲線,可分為三個階段:第一階段流量從零開始緩慢上升;第二階段流量近似于直線快速上升;第三階段流量變化趨勢與第一階段相似,最終趨于直線穩(wěn)定狀態(tài)。

圖2 1 s啟動時間下離心泵啟動過程瞬態(tài)外特性曲線Fig.2 The external characteristic curves of the centrifugal pump during its at t=1 s

瞬態(tài)計算,一般采用準穩(wěn)態(tài)假設(shè)方法代替瞬態(tài)過程。然而對于離心泵快速啟動過程,這種方法忽略了離心泵啟動過程中流體加速度和葉輪旋轉(zhuǎn)加速度對內(nèi)部流場的影響,瞬態(tài)特性與穩(wěn)態(tài)特性存在的誤差較大。本文將1 s啟動時間下離心泵啟動過程瞬態(tài)的流量曲線作為數(shù)值模擬的邊界條件(見圖2),以使得數(shù)值模擬結(jié)果較為準確。

2數(shù)值模擬方法

2.1流場數(shù)值模擬方法

流場計算包括葉輪和壓水室流道流體區(qū)域。計算時葉輪區(qū)域的流場采用旋轉(zhuǎn)坐標系下的控制方程求解,其他區(qū)域采用靜止坐標系控制方程求解。根據(jù)湍流的雷諾方程理論,連續(xù)性方程和動量方程的張量形式表達式為

(1)

(2)

對于離心泵內(nèi)部的流場計算,為了封閉流場求解方程,必須引用湍流模型。SST(Shear Strain Transport)模型的相關(guān)文獻已經(jīng)很多[11],這里不再贅述。

2.2結(jié)構(gòu)動力方程

離心泵在運行過程中,葉輪會受到內(nèi)部流體的反作用力,由于流體壓力分布是變化的,因此葉輪會隨按時間變化的載荷作用產(chǎn)生響應(yīng),而且這種響應(yīng)受結(jié)構(gòu)慣性力和阻尼作用比較顯著。根據(jù)哈密爾頓原理,彈性體的結(jié)構(gòu)動力學方程為

(3)

2.3雙向流固耦合的求解過程

雙向流固耦合的求解過程中數(shù)據(jù)交換是雙向的,即將流體分析結(jié)果傳遞給固體結(jié)構(gòu)分析,固體結(jié)構(gòu)分析的結(jié)果又反向傳遞給流體分析。

本文應(yīng)用CFX軟件對流場進行非定常數(shù)值計算,網(wǎng)格變形采用軟件提供的動網(wǎng)格技術(shù),運用ANSYS軟件對結(jié)構(gòu)的瞬態(tài)動力學進行分析,采用MFS(Multi-Field Solver)功能實現(xiàn)流場數(shù)據(jù)與結(jié)構(gòu)場數(shù)據(jù)的實施交換。求解過程見圖3。

圖3中的Δt為計算過程中的時間步長,流固耦合的計算要求固體域設(shè)置的時間步長必須與流體域設(shè)置的步長相一致,所以在設(shè)置求解時間步長時,需要對流場和結(jié)構(gòu)場進行綜合考慮。

圖3 雙向流固耦合過程Fig.3 FSI simulations with two-way coupling

3計算模型、網(wǎng)格劃分和邊界條件

3.1計算模型

本文選用IS65-50-160低比轉(zhuǎn)速離心泵為研究對象,主要設(shè)計參數(shù)如下:流量Qd=25 m3/h,揚程H=32 m,轉(zhuǎn)速n=2 900 r/min,比轉(zhuǎn)速ns=65.5。模型泵過流部件的主要幾何參數(shù)見表1

表1 過流部件主要設(shè)計參數(shù)

3.2網(wǎng)格劃分

流體計算域包括進水段、前泵腔、葉輪、蝸殼、后泵腔和出水段,結(jié)構(gòu)區(qū)域包括葉輪和泵軸。利用ICEM軟件對模型泵計算區(qū)域進行網(wǎng)格生成,考慮到蝸殼和葉輪流道內(nèi)流動的復(fù)雜性,流體區(qū)域采用自適應(yīng)性比較強的非結(jié)構(gòu)化四面體網(wǎng)格實現(xiàn)復(fù)雜結(jié)構(gòu)網(wǎng)格劃分。為了準確的進行數(shù)值計算,對葉輪進口邊位置、蝸殼隔舌壁面附近等計算部位進行局部加密,固體區(qū)域采用六面體網(wǎng)格為主的網(wǎng)格劃分方式。其中,流體區(qū)域網(wǎng)格總數(shù)為1 472 548,固體區(qū)域有限元的網(wǎng)格節(jié)點數(shù)為57 877,單元數(shù)為30 908,見圖4。泵軸與葉輪的材料分別為45鋼和HT200,葉輪的特性參數(shù)分別為密度ρ=7 800 kg/m3,泊松比μ=0.25,彈性模量Ε=122 GPa。

圖4 流場與結(jié)構(gòu)場計算模型Fig.4 The computational flow model of the pump

計算總時間設(shè)置為1 s,時間步長設(shè)置為0.000 8 s。邊界條件采用總壓進口,質(zhì)量流量出口,靜止區(qū)域蝸殼壁面設(shè)置為無滑移壁面;出口邊界條件和葉輪旋轉(zhuǎn)的轉(zhuǎn)速采用CFX自帶的CEL語言加載實驗得到的流量曲線來加以控制;湍流模型采用SST模型,空間離散為二階精度。

4結(jié)果分析

4.1耦合作用下的穩(wěn)定狀態(tài)分析

離心泵在穩(wěn)定工況下運行時,不僅結(jié)構(gòu)體對流體有制約的作用,而且泵內(nèi)部復(fù)雜的非定常流動產(chǎn)生的載荷也會反作用于結(jié)構(gòu)體,產(chǎn)生動應(yīng)力分布。動應(yīng)力的存在對結(jié)構(gòu)的安全產(chǎn)生威脅,會引起結(jié)構(gòu)疲勞破壞,因此,本節(jié)基于流固耦合作用主要探討了穩(wěn)定工況下離心泵葉輪應(yīng)力分布和葉片的變形情況。

圖5為葉片的等效應(yīng)力分布,葉輪內(nèi)部等效應(yīng)力分布不均勻,最大等效應(yīng)力為466 200 Pa;在葉輪的進口邊,葉片與前蓋板的交界處出現(xiàn)了應(yīng)力集中,這可能是因為在葉輪的進口水流速度大,葉片受到了水流的沖擊作用。

圖5 葉片表面應(yīng)力分布Fig.5 The stress distribution on the blades

圖6為葉片在穩(wěn)定工況下的變形量分布,最大變形量位于葉片與前蓋板交界處靠近出口的位置,最大變形量為6.797 μm。葉片主要表現(xiàn)為流場壓力產(chǎn)生的彎曲和扭轉(zhuǎn)變形,而由于離心力產(chǎn)生的拉伸變形并不明顯,表明離心力比流體作用力影響小。

圖6 葉片變形分布Fig.6 The deformation distribution of the impeller

4.2耦合作用下的瞬時應(yīng)力分析

圖7為離心泵啟動1 s內(nèi)葉輪最大等效應(yīng)力的變化曲線。見圖所示,隨著葉輪旋轉(zhuǎn)轉(zhuǎn)速的增加,最大等效應(yīng)力呈現(xiàn)振蕩上升趨勢。在時間為0.1 s之前,最大等效應(yīng)力較波動強烈,這可能是因為啟動初始階段內(nèi)部流場較混亂,內(nèi)部二次流,動靜干涉,旋渦等復(fù)雜流動造成的;在0.2~0.35 s內(nèi),最大等效應(yīng)力波動減弱;在0.35 s以后,最大等效應(yīng)力隨著轉(zhuǎn)速的增大而逐漸增大,波動幅值也逐漸變大,當轉(zhuǎn)速到達額定轉(zhuǎn)速時,最大等效應(yīng)力達到最大值。

圖7 最大等效應(yīng)力隨時間的變化曲線Fig.7 The changing curve of maximum equivalent stress with time

圖8表明隨著葉輪旋轉(zhuǎn)轉(zhuǎn)速的提高,葉片應(yīng)力的變化趨勢是先減小后增大,這與葉輪整體應(yīng)力變化趨勢相同。由于葉輪自身結(jié)構(gòu)的不對稱性,使得葉輪在啟動加速過程中受力不均勻,葉片進口受到的等效應(yīng)力明顯大于葉片出口。葉片末端受力不大,但是越接近應(yīng)力集中點,應(yīng)力梯度越明顯。在t=0.337 s之前,在葉片與后蓋板相交處靠近進口的位置出現(xiàn)葉片應(yīng)力集中現(xiàn)象,這是因為離心泵啟動初始階段,葉輪的轉(zhuǎn)速較低和葉輪進口流體的慣性作用,使得葉片的進口與后蓋板接合處為葉片受力的支撐點,因而此處的應(yīng)力大于葉輪進口葉片與前蓋板接合處的應(yīng)力;在t=0.337 s之后,應(yīng)力集中出現(xiàn)在葉片與前蓋板相交處靠近進口的位置,應(yīng)力集中的地方受到流體作用力的大小為此時最大應(yīng)力值,即為葉輪最容易產(chǎn)生疲勞破壞的位置。

4.3耦合作用下的變形分析

圖9中可知,最大變形量的變化曲線與最大應(yīng)力變化曲線相似,在離心泵啟動起始階段0~0.1 s內(nèi),最大變形量波動較為劇烈;在0.1 s左右時,啟動過程振蕩的峰值幾乎為零;0.1 s之后,葉片最大變形量隨著轉(zhuǎn)速的增加呈現(xiàn)振蕩上升的趨勢,并且振動幅值也不斷的增大;當轉(zhuǎn)速達到額定轉(zhuǎn)速時,曲線整體趨于水平,此時振動幅值達到最大值。

圖8 葉片等效應(yīng)力分布Fig.8 The equivalent stress distribution of impeller

圖9 葉片最大變形量隨時間的變化曲線Fig.9 The changing curve of the maximum deformation of the impeller

圖10所示為離心泵啟動過程中三個不同時刻葉片變形量分布。從圖10可知:在啟動過程中,葉片的總變形量先減小后增大,葉片的變形量分布不均勻;葉片的工作面與背面的變化情況基本一致;最小變形量位于葉片與后蓋板的交界處靠近出口的位置,最大變形量出現(xiàn)在葉片與前蓋板的交界處靠近出口的位置,可能是因為此處離心力和總壓相對較大,使得葉尖的變形最大;整個加速過程中,葉片的變形量從后蓋板到前蓋板呈遞增的趨勢,并且越靠近出口處,這種遞增的趨勢越明顯。

圖10 葉片總變形分布Fig.10 Total deformation of the impeller distribution

4.4葉片的強度分析

通過上文對葉片應(yīng)力分析可以得到:在離心泵啟動階段,葉片的最大等效應(yīng)力發(fā)生在葉片的進口處葉片與前蓋板相交的位置;最大等效應(yīng)力大于泵運行在設(shè)計工況下的最大應(yīng)力。因此對葉輪強度校核時,葉輪受到的最大應(yīng)力應(yīng)為啟動階段流體作用在葉輪上的最大等效應(yīng)力。

構(gòu)件強度校核與材料和所受載荷類型有關(guān),本文中葉片的材料為HT200,材料的許用應(yīng)力[σ]=80 MPa。強度校核要求:構(gòu)件的最大工作應(yīng)力應(yīng)小于或等于許用應(yīng)力狀態(tài)下構(gòu)件的強度條件。即

σmax≤[σ]

(4)

式中:σmax為構(gòu)件最大工作應(yīng)力;[σ]為構(gòu)件的許用應(yīng)力。

由此可見,葉輪在離心泵啟動過程中滿足強度要求。

5結(jié)論

本文通過雙向流固耦合方法對離心泵啟動過程瞬態(tài)的應(yīng)變和應(yīng)力進行了分析,得出了以下結(jié)論。

(1) 在離心泵啟動過程中,葉片的等效應(yīng)力和應(yīng)變量均呈振動上升趨勢,振動強度先減小后增大。

(2) 在離心泵啟動初始階段,葉片的等效應(yīng)力和應(yīng)變量變化頻率較高。葉片進口與后蓋板的交界處出現(xiàn)應(yīng)力集中,葉片的變形量從后蓋板到前蓋板呈遞增的趨勢,越靠近出口處,這種遞增的趨勢越明顯。

(3) 離心泵啟動過程中,葉片的等效應(yīng)力和應(yīng)變量的最大值都大于穩(wěn)態(tài)工況下的等效應(yīng)力和應(yīng)變量。

因此,在實際應(yīng)用過程中,應(yīng)盡量避免泵頻繁的啟動,防止由于頻繁的啟動造成葉片進口疲勞破壞或者葉片外緣產(chǎn)生較大的變形,導(dǎo)致葉片失效,影響泵正常運行。

[1] 關(guān)醒凡.泵的理論與設(shè)計[M].北京:機械工業(yè)出版社,1986.

[2] 袁壽其. 低比速離心泵理論與設(shè)計[M].北京:機械工業(yè)出版社,1997.

[3] Bathe K J,Zhang H.A flow-condition-based interpolation finite element procedure for incompressible fluid flows [J].Computers & Structures, 2002,80:1267-1287.

[4] Kohno H,Bathe K J.A nine-node quadrilateral FCBI element for incompressible fluid flows [J].International Journal for Numerical Methods in Fluids,2006,51: 673-699.

[5] Brennen C E.Hydrodynamics and Cavitation of pumps[M]. Vienna: Springer, 2008.

[6] Kato C,Yamade Y,Wang Hong, et al. Prediction of the noise from a multi-stage centrifugal pump [C]//ASME Fluilds Engineering Division Summer meeting,PART B. Houston, TX,2005.

[7] Benra F K,Dohmen H J. Comparison of pump impeller orbit curves obtained by measurement and FSI simulation [C]//ASME PVP2007.San Antonio,TX,2007.

[8] Langthjem M A.A numerical study of flow-induced noise in a two-dimensional centrifugal pump,Part I:hydrodynamics [J].Fluid and Structures,2004,19(3):349-368.

[9] Guadaqni,Gualtiero,F(xiàn)iore,et al. A fluid-structure analysis of the structure and fluid dynamic behavior of a new disposable pulsatile pump for cardiopulmonary bypass [J].American Society of Mechanical Engineers,Bioengineering Division(Publication) BED,2001,50:220-230.

[10] 王洋,王洪玉,徐小敏,等.沖壓焊接離心泵葉輪有限元計算[J].排灌機械工程學報,2011,29(3):109-113.

WANG Yang,WANG Hong-yu, XU Xiao-min,et al. Finite element computation for impeller of stamping and welding centrifugal pump[J].Journal of Drainage and Irrigation Machinery Engineering, 2011,29(3):109-113.

[11] 王國玉,霍毅,張博,等.湍流模型在軸流泵性能預(yù)測中的應(yīng)用與評價 [J].北京理工大學學報,2009,29 (4) :309-313.

WANG Guo-yu,HUO Yi, ZHANG Bo, et al. Evaluation of turbulence models for predicting the performance of an axial-flow pump[J].Transactions of Beijing Institute of Technology, 2009,29 (4) :309-313.

[12] 陳向陽,袁丹青,楊敏官,等. 基于流固耦合方法的300MWe級反應(yīng)堆主泵葉片應(yīng)力分析[J].機械工程學報,2010,46(3):111-115.

CHEN Xiang-yang,YUAN Dan-qing,YANG Min-guan, et al.Blade stress of the reactor coolant pump of 300 MW nuclear power plant in China based on fluid-solid coupling method[J].Journal of Mechanical Engineering,2010,46(3):111-115.

[13] Xu H,Tan M G,Liu H L,et al.Fluid-structure interaction study on diffuser pump with a two-way coupling method [C]//5th International Symposium on Fluid Machinery and Fluids Engineering.Jeju,Korea,2012.

[14] 鄭小波,羅興琦,鄔海軍.軸流式葉片的流固耦合振動特性分析[J].西安理工大學學報,2005,21(3):342-346.

ZHENG Xiao-bo,LUO Xing-qi,WU Hai-jun.Analysis of fluid-solid coupling dynamic characteristics for the axial flow blades[J].Journal of Xi’an University of Technology,2005,21(3):342-346.

[15] 張麗霞,張偉,潘際奎,等. 基于流固耦合理論的混流式葉片動力學分析[J].清華大學學報:自然科學版,2008,48(5):773-776.

ZHANG Li-xia,ZHANG Wei,PAN Ji-kui, et al.Dynamic analysis of a Francis turbine based on the fluid-structure interaction theorem[J].Journal of Tsinghua University Science and Technology,2008,48(5):773-776.

[16] Saito S.The transient characteristics of a pump during start up [J]. Bulletin of the JSME, 1982,201(25):372-379.

[17] 吳大轉(zhuǎn),許斌杰,李志峰,等.離心泵瞬態(tài)操作條件下內(nèi)部流動的數(shù)值模擬[J].工程熱物理學報,2009,30(5):781-783.

WU Da-zhuan,XU Bin-jie,LI Zhi-feng,et al.Numerical simulation on internal flow of centrifugal pump during transient operation[J].Journal of Engineering Thermophysics,2009,30(5):781-783.

[18] 王學. 基于ALE方法求解流固耦合問題[D].長沙:國防科學技術(shù)大學,2006.

Transient stress characteristic during centrifugal pumps start-up based on fluent-structure interaction

YUAN Jian-ping, XIA Shui-jing, ZONG Wei-wei, ZHOU Bang-lun, FU Yan-xia

(Research Center of Fluid Machinery Engineering and Technology, Jiangsu University, Zhenjiang 212013, China)

In order to predict the stress and deformation during the transient state of the centrifugal pump startup period, a coupled solution of a flow field in the pump and a structural response of the blades was established using a Fluid-Structure Interaction method. The flow field prediction was based on the Reynolds-averaged N-S equations and the SST turbulence model, while the structure prediction was based on elastic structural dynamic equation. Based on the transient numerical simulation results, the Fluid-Structure Interaction method was used to solve the internal flow field and the structure field during the centrifugal pump startup period. The results showed that: the maximum equivalent stress and strain of the blade show a vibration that increases; the vibration intensity decreases first and then increases; an intense increase occurs in the impeller inlet near the hub; the blade deformation shows a trend that increases from the hub to the shroud; and the maximum equivalent stress and strain of the blades are greater than those in the steady-state conditions.

centrifugal pump; fluid-structure interaction; stress distribution; strain; numerical simulation

10.13465/j.cnki.jvs.2016.12.031

2015-04-10修改稿收到日期:2015-07-02

袁建平 男,博士,研究員,1970年12月生

TH212;TH213.3

A

猜你喜歡
效應(yīng)變形
鈾對大型溞的急性毒性效應(yīng)
懶馬效應(yīng)
場景效應(yīng)
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
“我”的變形計
變形巧算
例談拼圖與整式變形
應(yīng)變效應(yīng)及其應(yīng)用
會變形的餅
偶像效應(yīng)
主站蜘蛛池模板: 尤物成AV人片在线观看| 免费国产不卡午夜福在线观看| 午夜性爽视频男人的天堂| 亚洲精品你懂的| a级毛片免费看| 精品国产成人高清在线| 浮力影院国产第一页| av大片在线无码免费| 亚洲无码A视频在线| www.91在线播放| 国产噜噜噜| 久久一色本道亚洲| 国产成人91精品| 午夜日韩久久影院| 人人看人人鲁狠狠高清| 精品国产一二三区| 色综合久久88色综合天天提莫| 18禁黄无遮挡免费动漫网站| 色婷婷成人| 黄色一级视频欧美| 国产中文在线亚洲精品官网| 亚洲精品成人片在线播放| 黄片一区二区三区| 国内黄色精品| 久久精品欧美一区二区| 亚洲自拍另类| 国产色婷婷| 亚洲中文字幕无码mv| a级毛片免费看| 国产男女XX00免费观看| 亚洲人成在线精品| 一级香蕉视频在线观看| 一级毛片网| 新SSS无码手机在线观看| 国产微拍一区二区三区四区| AV无码国产在线看岛国岛| 国产乱子伦手机在线| 亚洲精品国产精品乱码不卞| 91成人精品视频| 蝴蝶伊人久久中文娱乐网| 色综合a怡红院怡红院首页| 国产经典三级在线| 在线无码九区| 一区二区三区四区精品视频| av在线5g无码天天| 8090成人午夜精品| 国产成人超碰无码| 精品久久久久成人码免费动漫| 亚洲成A人V欧美综合| 欧美日韩导航| 精品無碼一區在線觀看 | 免费a级毛片18以上观看精品| 久久国产精品电影| 精品少妇人妻一区二区| 亚洲综合九九| 久久精品免费国产大片| 亚洲第一色视频| 91视频区| 国产一区在线视频观看| 全裸无码专区| 亚洲黄色视频在线观看一区| 少妇精品久久久一区二区三区| 美女无遮挡免费网站| 久久美女精品国产精品亚洲| 欧美激情视频二区| 国产精品3p视频| 国产高清在线观看91精品| 三级毛片在线播放| 成人a免费α片在线视频网站| 国产精品自在线拍国产电影| 2048国产精品原创综合在线| 亚洲 欧美 偷自乱 图片 | 国产人成午夜免费看| 日本a级免费| 国产精品女主播| 国产麻豆另类AV| 午夜免费视频网站| 99热6这里只有精品| 欧美一区二区三区不卡免费| 99精品免费在线| 88av在线播放| 成年免费在线观看|