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

懸臂式離心泵流固耦合特性研究

2016-09-16 01:20:10牟介剛陳瑩谷云慶鄭水華錢亨
哈爾濱工程大學學報 2016年8期
關鍵詞:變形

牟介剛,陳瑩,谷云慶,鄭水華,錢亨

(浙江工業大學 機械工程學院,浙江 杭州 310014)

?

懸臂式離心泵流固耦合特性研究

牟介剛,陳瑩,谷云慶,鄭水華,錢亨

(浙江工業大學 機械工程學院,浙江 杭州 310014)

針對懸臂式離心泵流固耦合問題,采用單向耦合方法,在不同懸臂比、不同工況下,分析葉輪應力、應變的變化規律,獲取離心泵在運轉過程中最容易發生強度破壞和剛度破壞的位置,并在有預應力和無預應力下,對不同懸臂比離心泵的轉子動力學特性進行研究。結果表明:在一定工況下,懸臂比越小,葉輪的變形越小,葉輪后蓋板中間流道處最有可能發生剛度破壞,葉片出口邊緣與前、后蓋板交接處最有可能發生強度破壞;懸臂比越小,離心泵轉子系統的固有頻率越大,流固耦合作用會降低轉子系統的前三階固有頻率。

離心泵;懸臂比;流固耦合;應力應變;轉子動力學;固有頻率;數值模擬

網絡出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20160624.1127.006.html

懸臂式離心泵作為離心泵應用最廣泛的一種結構形式,其轉子泵軸的兩個支撐軸承均位于泵軸的一端,葉輪則安裝在泵軸的另一端,處于自由懸臂狀態。這種特殊結構,使得其工作狀態參數(轉速、外載荷)和實際結構參數(懸臂長度、支撐間距等)對轉子動力特性的綜合影響變得更為復雜。離心泵在高速旋轉過程中,蝸殼內流場狀態不斷變化,勢必造成葉輪所受載荷隨之變化而產生流體激勵力,繼而引起懸臂式離心泵轉子的振動。同時周期性的液力載荷又會引起葉輪及泵軸的動態變形,從而進一步影響流體的分布,直接導致離心泵運行效率低,甚至引發故障[1]。

懸臂式離心泵內部流動為復雜的三維非定常湍流流動,常伴有渦流、回流、汽蝕、水力振動等現象[2]。隨著計算流體力學的發展,運用CFD研究離心泵內部流場已成為一種主流手段,在內流場與壓力脈動的計算已經取得了一些成果[3-4]。Kelder等[5]通過勢流理論和實驗兩方面研究了離心泵內外流場的特性。Kaewnai等[6]利用CFD研究不同湍流模型對離心泵外特性的影響。當前,在旋轉機械中流固耦合問題已成為研究熱點之一,如羅永要等[7]對混流式水輪機動載荷作用下的應力特性進行研究,結果表明動應力過大導致葉片產生裂紋。Young[8]對彈性的復合海洋螺旋槳進行研究,指出葉片的變形改變液體壓力分布、空化程度,最終降低螺旋槳的效率。Brennen等[9]通過理論和試驗分析了懸臂式離心泵流體誘導轉子產生動應力,嘗試在前蓋板上增加溝槽或凸條以降低入口處渦流率,進而增加轉子系統的穩定性,但效果并不明顯。吳賢芳等[10]對離心泵關死點內部流場研究表明,流固耦合作用對揚程的影響較小,對壓力脈動影響較大。Zhang等[11]基于流固耦合方法研究了不同葉片形狀的瞬態振動性能。劉厚林等[12]對比分析單、雙向耦合的葉輪靜應力強度,指出單向耦合即可滿足葉輪靜應力分析。

泵軸結構參數是離心泵運行穩定性的重要因素,針對懸臂式離心泵的特殊結構,其懸臂比的設計對離心泵動力學特性有重要影響,但目前泵軸結構對轉子動力學影響研究鮮有報道。采用單向流固耦合計算方法研究懸臂式離心泵靜力學特性和轉子動力學特性,分析懸臂比對離心泵穩定性的影響。

1 參數建模及數值計算方法

1.1離心泵建模

離心泵流固耦合特性研究過程中,以IS 80-50-250型號離心泵為研究對象,其結構示意圖如圖1(a)。

(a)單級單吸懸臂式離心泵

(b)懸臂比圖1 離心泵結構示意圖Fig.1 Structure diagram of centrifugal pump

離心泵參數為:流量Q=50 m3/h,揚程H=80 m,轉速n=2 900 r/min,進口直徑D1=80 mm;葉輪為閉式葉輪,葉片數為5,葉輪外徑D2=252 mm;葉輪材料選為HT200,密度ρ=7 200 kg/m3,彈性模量E=2.07×1011Pa,泊松比γ=0.3;泵軸材料選為45,密度ρ=7 890 kg/m3,彈性模量E=1.1×1011Pa,泊松比γ=0.28。離心泵懸臂比為泵軸伸出端與軸承支撐間距之比,如圖1(b)所示,則離心泵的懸臂比λ為

(1)

式中:l1為泵軸伸出端長度,l2為軸承支撐間距。

1.2網格劃分

離心泵數值計算模型流體域包括進水管、葉輪水體、蝸殼水體和出水管4部分,采用網格劃分軟件ICEM對流體域進行網格劃分。因葉輪和蝸殼水體結構復雜,造成離心泵非結構化網格單元數量多,進而帶來數值計算量大、耗時長的問題,故采用拓撲映射法對流體域劃分結構化六面體網格,使網格生成速度快、質量好、數據結構簡單[13]。進水管與出水管部分,因其結構簡單,對液體流入面和流出面進行O型劃分,即得到高質量、正交性良好的結構化網格;葉輪水體部分,結構雖相對復雜,但呈現循環對稱的特點,因此將葉輪平均分割成5個部分,只需建立其中一個周期的塊并做出符合實際流體流向的六面體網格;蝸殼水體部分,其結構復雜部位集中在隔舌處,可將隔舌部位的六面體網格流向分為3個方向,分別是出水口方向、第一斷面方向和第八斷面方向,確保網格流向和實際流體流向相同,同時在隔舌處進行網格加密。在流體壁面的邊界層節點應用linear排列方式,加密邊界層使流固耦合交換信息誤差減小。分別用4組不同尺度的網格對設計工況進行數值模擬,對比模擬揚程值來進行網格無關性驗證,如表1所示不同網格數量對應的模擬揚程,從表中可看出計算的揚程值與試驗值間的誤差較小,均在2%內。因此在滿足計算精度的前提下考慮到節省時間,本文最終選用網格數量為110萬左右。固體域包括固體葉輪和泵軸2部分,采用ANSYS自帶網格劃分軟件對其進行非結構化網格劃分,網格數為188 763,節點數為297 564。計算模型網格劃分情況如圖2所示。

表1 不同網格數量的揚程模擬值

1.3計算方程

懸臂式離心泵在數值計算中,內部流體可視為等溫不可壓縮流體,內部流動近似為常密度流動。采用三維、穩態、不可壓的連續性方程及RNGk-ε湍流方程進行求解[14]。其中連續方程為

(2)

動量方程為

(3)

式中:ρ為密度,ui為平均相對速度,xi為各坐標分量,t為時間。

懸臂式離心泵內部流線彎曲程度大,采用RNGk-ε湍流模型能很好地處理內部流動中旋轉和旋流情況,在RNGk-ε湍流模型中湍動能k和湍流耗散率ε的方程為

(4)

(5)

圖2 網格模型Fig.2 Grid model

1.4邊界條件設置

研究過程中分別選取0.6Q、0.8Q、1.0Q、1.2Q、1.4Q五種工況,采用ANSYS CFX中進行數值模擬,進出口分別設置為質量流量進口、自由流出口。葉輪采用旋轉坐標系,蝸殼、進出水管采用靜坐標系;定常計算時,動靜交接面采用凍結轉子??刂品匠炭臻g離散使用基于有限元體積法,對流項選擇高分辨率,收斂目標設置為10-4。采用標準壁面函數處理近壁面,固體壁面設為無滑移。輸送介質為常溫常壓條件下的清水。

1.5單向耦合過程

基于單向耦合方法,對離心泵λ分別為0.8、0.9、1.0、1.1、1.2、1.3等6種懸臂比下的葉輪應變應力分布情況進行分析。耦合過程在ANSYS Workbench中實現,為將離心泵內部流體力準確得傳遞到固體葉輪表面,需定義流固耦合面,將流體定常計算結果以壓力載荷的形式施加到流固耦合面上,定義葉輪水體與固體所有接觸面為流固耦合面,即前、后蓋板面以及葉片表面[15];并對轉子結構施加慣性載荷(重力載荷和離心載荷);假定軸承支承為剛性支承,其安裝處添加圓柱約束,軸向和徑向固定,切向自由。

離心泵結構強度計算的靜力學方程[16]:

(6)

(7)

式中:K為剛度矩陣,D為彈性矩陣,B為應變矩陣,δ為位移,F為所受的力,σ為應力。

根據第四強度理論計算等效應力σe:

(8)

式中:σ1為第一主應力,σ2為第二主應力,σ3為第三主應力。

2 葉輪應變應力分析

2.1葉輪應變分析

圖3為不同懸臂比λ下葉輪的最大總變形隨流量變化曲線圖。由圖3可知,葉輪在離心泵徑向力、軸向力和扭矩共同作用下發生應變,總變形的最小值和最大值分別處于λ=0.8在Q=40m3/h和λ=1.3在Q=70m3/h下。各個工況下,最大總變形量與懸臂比的大小呈正相關,懸臂比越小,葉輪變形越小。不同懸臂比的離心泵運行時,葉輪產生的最大總變形隨工況的增大而呈現先下降后上升的變化趨勢。這是因為不同工況下離心泵徑向力、軸向力、扭矩都隨流量而變,由離心泵定常計算結果可知,軸向力隨流量增大而減小,徑向力隨流量先減小后增大,設計工況下徑向力最小。扭矩隨流量增大而增大,在三者綜合載荷作用下葉輪總變形先下降后上升。從最大總變形的變化趨勢可以看出,大懸臂比對不同工況的敏感程度更高,即改變運行工況時大懸臂比離心泵的葉輪變形梯度更大,小懸臂比的變形量較平緩。各個工況下,不同懸臂比的離心泵運轉時葉輪產生的變形量方差依次為1.93×10-4、6.3×10-5、3.3×10-5、7.55×10-5、1.82×10-4。在設計工況點變形量較為集中,區間為0.19~0.21mm;越偏離設計工況,葉輪最大變形量隨懸臂比的改變發生的變化越大。故盲目增大懸臂比會加劇離心泵不穩定運行,內部流體激勵力引起葉輪較大變形,反過來較大變形影響內部流場。為了保證葉輪運轉過程中的剛度要求,避免選擇較大的懸臂比,且避免懸臂泵在過大偏離設計工況點下運行。

圖4 為在設計工況點不同懸臂比λ下葉輪的應變分布云圖??芍?,葉輪發生了彎曲變形,變形幾乎成軸對稱,并且沿著葉輪半徑增大的方向逐漸增加。葉輪后蓋板出口邊緣中間流道處發生最大變形,此處壓力最大,說明此處最有可能發生剛度破壞。改變懸臂比對葉輪產生最大總變形有較明顯的影響。隨著懸臂比增加,葉輪小變形區域逐漸減小,變形量逐漸增大。此外,懸臂比取大值時葉輪的最大應變明顯高于懸臂比取小值,說明泵軸伸出端越短葉輪變形越小,離心泵運行過程越穩定;泵軸伸出端越長越容易引起葉輪后蓋板邊緣中間流道處的較大應變。

圖3 不同懸臂比不同載荷下葉輪的最大總變形Fig.3 The maximum total deformation of impeller under different cantilever ratios and loads

圖4 設計工況下葉輪應變分布圖Fig.4 Distribution of impeller total deformation under design condition

2.2葉輪應力分析

圖5為不同懸臂比λ下葉輪的最大等效應力隨流量變化曲線圖。由圖5可知,葉輪的最大等效應力均隨流量的增加而降低,降低速度先緩后急。最大等效應力為78.958MPa,小于葉輪材料HT200的許用應力200MPa,說明葉輪的強度足夠。在0.6Q工況下,葉輪的最大等效應力略有不同;當λ=1.0,最大等效應力最小,為77.73MPa;當λ=1.3時,最大等效應力最大,為78.958MPa,最大值與最小值相差1.2MPa。懸臂比的改變并不會引起葉輪的最大等效應力出現明顯變化,因此,對于懸臂式離心泵在小工況運轉時,應充分考慮葉輪的強度要求。

圖6 為在設計工況點不同懸臂比λ下葉輪的應力分布云圖。由圖6可知,葉輪的等效應力呈循環對稱分布,這與葉輪循環對稱結構相符(如圖中局部放大圖所示最大等效應力集中在葉片背面出口邊緣與前后蓋板相交接的位置),此處前后蓋板的壓力載荷最大。最大等效應力在59.11~59.13MPa,遠小于葉輪的許用應力,且不同懸臂比下等效應力的最大值與最小值的相對偏差在1%以內,說明離心泵懸臂比對葉輪最大等效應力值的影響很小。

圖5 不同懸臂比不同載荷下葉輪的最大等效應力Fig.5 The maximum equivalent stress of impeller under different cantilever ratios and loads

圖6 設計工況下葉輪等效應力分布圖Fig.6 Distribution of impeller equivalent stress under design condition

3 轉子動力學特性分析

轉子的穩定性限制了離心泵的高效運行范圍,同時也影響了離心泵的安全可靠性。在靜力學分析基礎上,對不同懸臂比懸臂式離心泵進行有預應力模態分析,并將有預應力模態與自由模態的結果進行對比,其中預應力指的是內部流體計算的壓力載荷、轉子系統的慣性力和離心力。表2為不同懸臂比下轉子系統的前6階固有頻率。

表2 不同懸臂比下葉輪固有頻率

由表2可知,隨著λ的增大,葉輪的各階固有頻率逐漸減小,且減小幅度逐漸降低,說明泵軸結構的外形尺寸對固有頻率有顯著的影響。轉子系統前3階固有頻率在有預應力的情況下明顯低于無預應力情況,后三階固有頻率在兩種情況幾乎相等,因在無預應力時忽略了阻尼作用和附加質量,有預應力時考慮了離心泵的離心力和流體力,其中流體力相當于附著在固體葉輪表面的附加質量,在離心泵運轉過程中固體葉輪與周圍流體產生摩擦,同時內部流體也產生摩擦,引起能量消耗,增大了阻尼作用,從而降低了轉子系統的固有頻率。在有預應力、λ=1.3情況下,葉輪的第一階固有頻率最小,為140.49 Hz,遠離葉片通過頻率(zn/60=5×2 900/60=241.5 Hz),故懸臂式離心泵在運轉過程中不會發生共振現象。

在有預應力、λ=0.8情況下離心泵轉子系統前6階振型如圖7所示。由圖7可知,前6階振型表現為葉輪整體振動,葉輪較軸端振動較大。第1階振型表現為繞z軸的扭轉,振動出現節圓。第2階和第3階振型相同,但角度分布不同,分別沿y、x軸擺動。葉輪第4、5階各出現1條節徑,節徑方向互相垂直,同時發生彎曲振動和扭曲振動。第6階振型為沿z軸的傘形彎曲振動。轉子系統的第2階和第3階、第4階和第5階振型相同固有頻率相近,這種振型成對出現的現象與葉輪循環對稱的結構有關,因循環對稱結構具有重根模態。轉子系統的振動幅度隨著階數的增加逐漸增長,階數越高節徑數越多,離心泵轉子模態頻率越大。無預應力下轉子振型表現形式與有預應力時基本一致。

圖7 有預應力下轉子系統前6階振型Fig.7 The first six mode of vibration considering pre-stress

4 結論

1)懸臂比與流量對葉輪應力應變有較大的影響,當流量一定時,懸臂式離心泵葉輪的應變隨著懸臂比的增大而增大,最佳懸臂比為0.8;當懸臂比一定時,葉輪的應變隨著流量的增加呈先減小后增大的變化趨勢,葉輪應變沿著半徑增大的方向逐漸增大,最大變形發生在葉輪后蓋板中間流道處。

2)懸臂比的改變并不會引起葉輪的最大等效應力出現明顯變化,葉輪最大等效應力隨流量的增加而減小,應力集中出現在葉片出口邊緣與前、后蓋板交接處。

3)隨著懸臂比的增加,葉輪固有頻率逐漸減小,且減小幅度逐漸降低;流固耦合作用能夠降低葉輪的固有頻率,葉輪最小固有頻率為140.49 Hz,不會發生共振現象。葉輪振型表現為扭曲、擺動、彎曲振動,且第2階和第3階、第4階和第5階振型相同固有頻率相近。

[1]ZHANG Desheng, XU Yan, LIU Yang, et al. Effects on the performance and flow field of the axial flow pump hydraulic components with fluid-solid interaction[J]. Journal of food, agriculture & environment, 2013, 11(3/4): 2000-2004.

[2]CHEAH K W, LEE T S, WINOTO S H. Unsteady analysis of impeller-volute interaction in centrifugal pump[J]. International journal of fluid machinery and systems, 2011, 4(3): 349-359.

[3]LUCIUS A, BRENNER G. Unsteady CFD simulations of a pump in part load conditions using scale-adaptive simulation[J]. International journal of heat and fluid flow, 2010, 31(6): 1113-1118.

[4]率志君, 張權, 陳春來, 等. 多級離心泵整機流場三維非穩態湍流壓力脈動特性分析[J]. 哈爾濱工程大學學報, 2013, 34(3): 306-311.

SHUAI Zhijun, ZHANG Quan, CHEN Chunlai, et al. Characteristic analysis of three-dimensional unsteady turbulent pressure fluctuation in whole flow field of multi-stage centrifugal pump[J]. Journal of Harbin Engineering University, 2013, 34(3): 306-311.

[5]KELDER J D H, DIJKERS R J H, VAN ESCH B P M, et al. Experimental and theoretical study of the flow in the volute of a low specific-speed pump[J]. Fluid dynamics research, 2001, 28(4): 267-280.

[6]KAEWNAI S,CHAMAOOT M,WONGWISES S. Predicting performance of radial flow type impeller of centrifugal pump using CFD[J]. Journal of mechanical science and technology, 2009, 23(6): 1620-1627.

[7]羅永要, 王正偉, 梁權偉. 混流式水輪機轉輪動載荷作用下的應力特性[J]. 清華大學學報:自然科學版, 2005, 45(2): 235-237, 257.

LUO Yongyao, WANG Zhengwei, LIANG Quanwei. Stress of francis turbine runners under fluctuant work conditions[J]. Journal of Tsinghua University: science and technology, 2005, 45(2): 235-237, 257.

[8]YOUNG Y L. Fluid-structure interaction analysis of flexible composite marine propellers[J]. Journal of fluids and structures, 2008, 24(6): 799-818.

[9]BRENNEN C E, ACOSTA A J. Fluid-induced rotordynamic forces and instabilities[J]. Structural control and health monitoring, 2006, 13(1): 10-26.

[10]吳賢芳, 談明高, 劉厚林, 等. 流固耦合作用對離心泵關死點內流的影響[J]. 應用基礎與工程科學學報, 2015, 23(1): 172-181.

WU Xianfang, TAN Minggao, LIU Houlin, et al. Effect of FSI on inner flow field in a centrifugal pump at shut off condition[J]. Journal of basic science and engineering, 2015, 23(1): 172-181.

[11]ZHANG Yu, HU Sanbao, ZHANG Yunqing, et al. Optimization and analysis of centrifugal pump considering fluid-structure interaction[J]. The scientific world journal, 2014, 2014: 131802.

[12]劉厚林, 徐歡, 吳賢芳, 等. 基于流固耦合的導葉式離心泵強度分析[J]. 振動與沖擊, 2013, 32(12): 27-30. LIU Houlin, XU Huan, WU Xianfang, et al. Strength analysis of a diffuser pump based on fluid-structure interaction[J]. Journal of vibration and shock, 2013, 32(12): 27-30.

[13]LIN Hongwei, TANG Kai, JONEJA A, et al. Generating strictly non-self-overlapping structured quadrilateral grids[J]. Computer-aided design, 2007, 39(9): 709-718.

[14]譚磊, 曹樹良, 桂紹波, 等. 帶有前置導葉離心泵空化性能的試驗及數值模擬[J]. 機械工程學報, 2010, 46(18): 177-182.

TAN Lei, CAO Shuliang, GUI Shaobo, et al. Experiment and numerical simulation of cavitation performance for centrifugal pump with inlet guide vane[J]. Journal of mechanical engineering, 2010, 46(18): 177-182.

[15]CAMPBELL R L, PATERSON E G. Fluid-structure interaction analysis of flexible turbomachinery[J]. Journal of fluids and structures, 2011, 27(8): 1376-1391.

[16]PEI Ji, YUAN Shouqi, YUAN Jianping. Dynamic stress analysis of sewage centrifugal pump impeller based on two-way coupling method[J]. Chinese journal of mechanical engineering, 2014, 27(2): 369-375.

本文引用格式:

牟介剛,陳瑩,谷云慶,等. 懸臂式離心泵流固耦合特性研究[J]. 哈爾濱工程大學學報, 2016, 37(8): 1111-1117.

MOU Jiegang, CHEN Ying, GU Yunqing, et al. Research on fluid-structure interaction characteristics of cantilever centrifugal pump[J]. Journal of Harbin Engineering University, 2016, 37(8): 1111-1117.

Research on fluid-structure interaction characteristics of cantilever centrifugal pump

MOU Jiegang, CHEN Ying, GU Yunqing, ZHENG Shuihua, QIAN Heng

(College of Mechanical Engineering, Zhejiang University of Technology, Hangzhou 310014, China)

To address cantilever centrifugal pump fluid-structure interaction problems, we used the one-way coupling method to analyze the stress and strain distribution of the impeller with different cantilever ratios and working conditions. We then determined the positions of the centrifugal pump that are most likely to be associated with rigidity and strength failure during operation. We also studied the prestressed and non-prestressed centrifugal pump rotor’s dynamic characteristics at different cantilever ratios. The results show that, under a certain condition, the smaller the cantilever ratio, the smaller the deformation of the impeller. In addition, the maximum total deformation usually occurs in the middle flow channel of the hub, the maximum equivalent stress usually occurs in the outlet of the blade and impeller cover junction, and the smaller the cantilever ratio, the greater the natural frequency of the centrifugal pump rotor system. Lastly, the first three natural frequencies of the rotor system are reduced by fluid-structure interaction.

centrifugal pump; cantilever ratio; fluid-structure interaction; stress and strain; rotor dynamics; natural frequency; numerical simulation

2015-06-10.網絡出版日期:2016-06-24.

國家自然科學基金項目(51476144);浙江省自然科學基金項目(LQ15E050005).

牟介剛(1963-), 男, 教授,博士生導師;

谷云慶(1982-), 男, 講師,博士.

谷云慶,E-mail:guyunqing@zjut.edu.cn.

10.11990/jheu.201506030

TH311

A

1006-7043(2016)08-1111-07

猜你喜歡
變形
變形記
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
柯西不等式的變形及應用
“變形記”教你變形
不會變形的云
“我”的變形計
會變形的折紙
童話世界(2018年14期)2018-05-29 00:48:08
變形巧算
例談拼圖與整式變形
會變形的餅
主站蜘蛛池模板: 亚洲国产系列| 亚洲国产精品日韩专区AV| 日韩大片免费观看视频播放| 亚洲一级毛片在线观播放| 亚洲精品国产日韩无码AV永久免费网| 国产粉嫩粉嫩的18在线播放91| 91精品啪在线观看国产91九色| 97视频免费在线观看| 伊人色综合久久天天| 亚洲色图欧美在线| 色综合天天操| 日韩福利在线观看| 色网站在线视频| 久久天天躁狠狠躁夜夜2020一| 国产美女丝袜高潮| 国产亚洲精品在天天在线麻豆| 中文字幕啪啪| 亚洲国产综合第一精品小说| 天天激情综合| 国产91精品久久| 亚洲人人视频| 97在线视频免费观看| 国产主播一区二区三区| 亚洲成在线观看| 国产婬乱a一级毛片多女| 亚洲区一区| 日韩精品成人在线| 成人福利免费在线观看| 亚洲va欧美va国产综合下载| 亚洲第一国产综合| 亚洲va欧美va国产综合下载| 国产日韩欧美视频| 国产精品无码一二三视频| 国产男女免费完整版视频| 精品三级网站| 扒开粉嫩的小缝隙喷白浆视频| 四虎亚洲精品| 国产jizz| 亚洲视频免费在线看| 青青草原偷拍视频| 日本一区高清| 成人无码区免费视频网站蜜臀| 日本三级欧美三级| 国产精品毛片一区| 国产一区二区三区夜色| 国产污视频在线观看| 超碰91免费人妻| 国产亚洲高清视频| 国产高清在线丝袜精品一区| 欧美另类第一页| 青青草原国产免费av观看| 在线a网站| 久久国产精品麻豆系列| 亚洲第一色视频| 亚洲国产成人精品一二区| 视频在线观看一区二区| 97人人模人人爽人人喊小说| 国产亚洲欧美在线人成aaaa| 伊人久久婷婷| 免费xxxxx在线观看网站| 亚洲精品国产首次亮相| 99久久亚洲精品影院| 国产精品女人呻吟在线观看| 一区二区三区精品视频在线观看| 久久精品免费看一| 久久精品国产在热久久2019| 国产精品亚洲一区二区三区z| 香蕉久久国产超碰青草| 国产二级毛片| 区国产精品搜索视频| 国产99视频精品免费观看9e| 91午夜福利在线观看精品| 日本在线国产| 国产香蕉国产精品偷在线观看| 99视频精品全国免费品| 一本大道无码高清| 99re精彩视频| 国产精品视频观看裸模| 精品第一国产综合精品Aⅴ| 91精品国产无线乱码在线 | 成人一区在线| 久久黄色小视频|