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

空化熱效應對上游泵送機械密封潤滑性能的影響

2016-10-25 05:47:50陳匯龍王彬任坤騰李同趙斌娟
化工學報 2016年10期
關鍵詞:模型

陳匯龍,王彬,任坤騰,李同,趙斌娟

?

空化熱效應對上游泵送機械密封潤滑性能的影響

陳匯龍,王彬,任坤騰,李同,趙斌娟

(江蘇大學能源與動力工程學院,江蘇鎮江 212013)

機械密封端面空化現象是影響機械密封潤滑性能的重要因素。采用計算流體動力學方法,基于Antoine公式,建立了考慮空化熱效應的計算流體動力學模型,并與常用的僅考慮端面液膜粘溫效應的模型進行對比,分析了空化熱效應對密封性能的影響。結果表明:在低轉速下,空化熱效應的影響可以忽略,但在高轉速下空化熱效應使上游泵送機械密封的高壓區形成能力減弱,泵送量降低,開啟力降低;在高轉速工作條件下分析密封失效機理時,除了考慮黏溫效應之外,還要考慮空化熱效應的影響;空化熱效應使螺旋槽槽區局部溫度比僅考慮黏溫特性時稍高;考慮空化熱效應時,動環端面空化發生程度最嚴重,由動環端面沿膜厚方向至靜環槽底,空化區域越來越小,槽底空化區域最小,此規律與僅考慮黏溫關系時相反。

機械密封;動壓潤滑;空化;熱力學性質;計算流體動力學;模型;穩定性

引 言

非接觸式液體機械密封端面間液膜的空化現象對液膜的承載力具有非常重要的影響,尤其在高轉速情況下,受流體變速運動和內摩擦溫升的共同作用,液膜空化及其對密封性能的影響更加突出,因此深入研究這種影響的機理和規律性具有重要的理論和實際意義。

密封端面間液膜的空化問題最早是由Findley[1]提出的。Jakobsson等[2-3]提出了一種自適應空化邊界條件來求解Reynolds方程,即JFO空化邊界條件,由于其符合質量守恒定律,經過Elord算法[4]的改進,使該邊界條件應用的可行性大大增加。在動力潤滑理論領域,求解基于雷諾方程并考慮空化問題的控制方程大多使用上述邊界條件和算法或其演化形式[5-9]。

對于液膜空化的研究,離不開液膜壓力和溫度。Pascovici等[10]認為等黏度下得到的密封間隙溫度高于變黏度下的溫度,基于等黏度的熱流體動力模型已不適用。Lebeck[11]認為熱效應對密封性能產生的影響和其他因素同等重要。這是由于熱效應改變了潤滑介質黏度,導致密封環變形,以及有可能使液膜汽化,從而導致泄漏量增加或密封失效[12]。由于介質黏度對溫度的敏感性,很多學者在建立密封潤滑的熱流體動力潤滑模型時,考慮了黏溫關系的影響。而空化時的飽和汽化壓力同樣也對溫度非常敏感,但在運用JFO邊界條件求解雷諾方程時,通常設定汽化壓力為大氣壓力或一定值[5-9],尚未有研究者考慮熱效應對飽和汽化壓力的影響問題。

本文以螺旋槽上游泵送機械密封為研究對象,以水為介質,在考慮黏溫關系的基礎上,利用克勞修斯-克拉貝隆方程基于實驗數據擬合得到的汽化壓力隨溫度的依變關系,即Antoine公式[13],建立考慮空化熱效應的計算流體動力學模型,研究不同轉速下空化熱效應對密封潤滑性能的影響。

1 計算模型

1.1 物理模型

本文研究的螺旋槽上游泵送機械密封結構如圖1所示,靜環材料為硬質合金,動環材料為碳石墨。在靜環端面開設螺旋槽,槽型線為對數螺旋線,如圖2所示。由于g個螺旋槽沿周向均布,因此取1/g液膜作為計算區域,如圖3所示。本文采用笛卡兒坐標系,設靜環端面為-平面,圓心為坐標原點,沿膜厚指向槽底方向為軸正方向。螺旋槽幾何參數和工況參數見表1。

圖1 密封結構示意圖

圖2 靜環端面螺旋槽結構

圖3 液膜計算域

表1 密封面幾何參數與工況參數

1.2 控制方程及數值求解

傳統的雷諾方程是求解流體動壓潤滑問題的基本方程,但雷諾方程是由N-S方程簡化推導而來的(如在液膜厚度方向忽略了壓力及黏度變化的影響等),由于流場各參數之間的相互耦合作用,有可能使流場流動細節分析不能夠準確把握,本文使用計算流體力學有限體積法,利用SIMPLEC算法對N-S方程、能量方程、氣體輸運方程、Reynolds黏溫方程和Antoine公式進行求解,主要方程如下。

1.2.1 N-S方程 液膜流場采用如下N-S方程計算獲得

1.2.2 氣體輸運方程 在空化過程中,液相和氣相的質量傳輸由如下氣體輸運方程控制

式中,e和c分別代表空化過程中液相和氣相之間的質量傳輸,可由描述單個氣泡生長過程的Rayleigh-Plesset方程[14]推導得到,表述如下。

式中,v為飽和汽化壓力。上述方程表明,當流場壓力小于汽化壓力時,氣泡產生,空化發生,反之,氣泡潰滅,空化消失。

1.2.3 能量方程

式中,S相為熱源項,表示液膜中由于密封環轉動產生的黏性摩擦熱。

1.2.4 Reynolds 黏溫方程[15]表征流體黏度隨溫度變化的關系為

1.2.5 Antoine公式 水的飽和汽化壓力和溫度的擬合曲線關系如圖4所示,由圖可見水的飽和汽化壓力對溫度非常敏感。本文利用克勞修斯-克拉貝隆方程,基于實驗數據擬合得到飽和壓力隨溫度的依變關系,即Antoine公式,如式(6)所示。

由上述方程及相應的邊界條件建立考慮空化熱效應的計算流體動力學模型。其中黏溫方程和表征空化熱效應的Antoine公式利用Fluent中的自定義函數(UDF)功能進行編程、編譯,并嵌入到相應的接口。結合SIMPLEC算法,整個計算模型的求解程序如圖5所示。

圖5 計算流體動力學模型求解程序

求解上述方程,需使用離散格式建立離散方程。為了保證方程求解的精度及速度,動量方程和能量方程使用二階迎風格式,氣體輸運方程使用QUICK格式。

1.3 邊界條件

由于螺旋槽上游泵送機械密封端面間微尺度液膜沿端面軸對稱分布,取圖3中周期液膜為計算單元,其周期性邊界條件為

空化邊界條件為:當壓力小于空化壓力時,空化產生,空化泡內壓力為飽和汽化壓力值,空化泡壁壓力梯度為0,即

為簡化計算及保證一定的計算精度,假設液膜與動、靜環端面的對流傳熱系數一致,并采用文獻[16]的方法計算;忽略慣性力和熱輻射的影響,不考慮端面微變形對溫度、流場的影響。

1.4 對比模型

為了研究空化熱效應對上游泵送機械密封潤滑性能的影響,本文選取以下兩種模型進行對比分析。模型1:只考慮黏溫效應,即耦合Reynolds黏溫方程,為一般熱效應分析法;模型2:考慮黏溫效應和空化熱效應,即耦合Reynolds黏溫方程和Antoine公式。

1.5 模型有效性檢驗

運用Gambit對單周期液膜進行網格劃分,采用Pave方法生成四面形面網格,然后采用Cooper方法拉伸生成六面體體網格。為了保證模擬結果的可靠性,運用 5種網格劃分方案進行網格劃分,并針對工況為轉速1000 r·min-1、進口壓力0.2 MPa進行了兩種對比模型的網格無關性檢驗,如表2所示。

表2 網格無關性檢驗

表中為單周期液膜網格單元數量;1和2分別為兩模型的單周期液膜開啟力;1和2分別為前后兩次開啟力的相對誤差。隨著網格數量的增加,兩模型的開啟力趨于穩定,當相對誤差小于0.5%時,認為網格數量的變化對計算結果的影響可以忽略,因此綜合考慮數值計算精度和時間成本,選用網格數為479037的模型進行計算。

李京浩[17]采用開放式的密封結構,同時利用透明材料制作其中一密封環,實驗觀察斜直線槽密封間隙中的水介質流動狀態,觀測出如圖6(a)所示的空化區域(圖中白色發亮部分)。圖6(b)為采用文獻[17]的參數,利用建立的空化熱效應模型對斜直線槽機械密封進行計算的結果。實驗觀測和數值計算結果表明,空化位置和形狀特征基本一致。

圖6 空化熱效應模型適用性驗證

2 計算結果及討論

2.1 空化熱效應對開啟力的影響

為便于開啟力的對比分析,定義相對開啟力增量為

式中,r為模型1的開啟力,ra為模型2的開啟力。

圖7為兩種模型的液膜開啟力隨轉速變化曲線圖。由圖7(a)可見,轉速由小增大過程中,液膜開啟力并沒有持續上升,而是存在極大值,對于本文研究對象及工況而言,極大值出現在轉速約為5000 r·min-1處。這一變化規律主要與液體黏溫特性有關,低轉速時,液膜內摩擦熱量較小,黏度變化較小,流體的動壓效應和泵送效應隨轉速的增強而增強,開啟力呈上升狀態;轉速超過5000 r·min-1之后,液膜內摩擦熱量增大,溫度升高,且隨著轉速的升高而明顯升高,黏度隨之下降,削弱了動壓效應和泵送效應,開啟力隨轉速逐步降低。由圖7(a)、(b)還可以看出,低轉速時,空化熱效應并不明顯,在常規的5000 r·min-1以下,空化熱效應與開啟力幾乎無關;轉速高于5000 r·min-1后,受空化熱效應的影響,開啟力呈現降低現象,且隨著轉速的增大開啟力的下降量增大,即值增大。這說明低轉速時,液膜空化程度低甚至未發生空化,隨轉速的提升空化不斷加劇,空化熱效應導致的開啟力負增量越來越明顯。

圖7 兩模型的開啟力對比

為了深入分析產生上述現象的原因及空化熱效應對密封潤滑性能的影響機理,本文將選取3000 r·min-1和13000 r·min-1兩個轉速對兩種模型的壓力場、溫度場、泵送量及空化區域進行對比分析。

2.2 空化熱效應對壓力分布的影響

圖8為兩種模型分別在轉速為3000 r·min-1和13000 r·min-1下的液膜厚度方向不同位置(-1.5 μm,0 μm,1.5 μm,3 μm,下同)的壓力分布云圖。兩模型的高壓區域均出現在槽根附近,且液膜軸方向壓力分布規律基本一致。由圖8(a)、(b)可知,在轉速為3000 r·min-1時,最大壓力max分別為4.78 MPa和4.77 MPa,大小基本相等,且兩者高壓區域的大小和變化梯度也基本一致;而在轉速為13000 r·min-1時,由圖8 (c)、(d)可知,模型2的最大壓力值max4.73 MPa,而模型1的最大壓力值為max5.12 MPa,模型2槽根高壓區最大壓力值和變化梯度明顯比模型1小,導致模型2的液膜開啟力為707.8 N,而模型1的液膜開啟力為811.7 N,如圖7(a)所示。

圖8 轉速為分別為3000 r·min-1和13000 r·min-1時兩模型壓力分布云圖

而在低壓區域,由圖8可知,模型1在槽的周向擴散側出現了明顯的局部低壓區,而模型2沒有這一現象。由圖8(a)、(c)可知,模型1的低壓區域增加明顯,這是由于在高轉速時,液膜流體周向速度加快,螺旋槽臺階擴散段相應的壓力下降區域增大。為了更好地顯示和分析低壓區的壓力分布情況,將轉速13000 r·min-1時兩模型的低壓區壓力分布進行細化,如圖9所示。由圖可知,兩圖的低壓區域大小基本一致,但模型1由于沒有考慮空化熱效應的影響,飽和汽化壓力為溫度293 K時的飽和汽化壓力2.3 kPa(絕對壓力)。而模型2低壓區域壓力明顯比模型1的高,而且低壓區域壓力也不再是一定值。這是由于密封轉動生成黏性摩擦熱導致液膜溫度升高,產生空化熱效應,水的飽和汽化壓力隨溫度的升高而明顯升高,即出現空化時相應的空泡內壓力也升高。由于兩者都考慮了空化的作用,因此飽和汽化壓力表征了液膜間隙各處溫度能出現的最小壓力值,從而出現了兩模型低壓區域分布不同的情況。

圖9 轉速為13000 r·min-1時兩模型低壓區壓力分布云圖

由上述分析可知,由于膜厚很小,壓力分布在膜厚方向基本沒有變化。雖然兩模型的低壓區壓力分布不同,模型2因考慮空化熱效應低壓區壓力增大,但最大的壓力僅為34 kPa,對開啟力的貢獻很有限,主要還是槽根處高壓區的壓力分布決定液膜開啟力的大小。在高轉速情況下,由于空化熱效應,導致模型2的高壓區形成能力減弱,從而使液膜的開啟力明顯減小。

2.3 空化熱效應對泵送量的影響

由圖10可知,在低轉速時,兩模型的泵送量都隨著轉速的增加呈先上升后逐漸趨于平緩直至下降的規律,這是因為隨著轉速的增加,液膜溫度升高,黏度降低,到達一定轉速后剪切力開始減小,泵送量減小。

圖10 兩模型在不同轉速下的泵送量

當轉速在10005000 r·min-1之間時,兩模型的泵送量基本相等;而當轉速繼續增加時,考慮空化熱效應的模型2的泵送量小于模型1,且隨轉速的增大其差值越來越大;通過泵送量變化規律與開啟力的變化趨勢[圖7(a)]的比較表明,當轉速高于5000 r·min-1后,雖然泵送量在增大,但動壓效應因溫度的升高和黏度的下降而減弱,開啟力呈下降趨勢;同時說明空化的熱效應在高轉速時使泵送效應減弱。

2.4 空化熱效應對液膜溫度的影響

由圖 11可知,兩模型的動環端面平均溫度都隨轉速的增加而增加,這是間隙液膜黏性耗散熱增加所致。

圖11 兩模型的動環端面平均溫度

圖12為兩模型液膜厚度方向不同位置溫度分布云圖。由圖可知,螺旋槽區域及附近溫度低于外徑處的溫度,這是因為一方面螺旋槽區域液膜厚度較大且內徑處線速度較小,導致速度梯度減小,黏性耗散熱減小;另一方面,螺旋槽的泵送作用使對流傳熱增強,導致液膜溫度較低且趨向平均,這與文獻[18]的分析一致。兩模型各個剖面的溫度分布基本一致,模型2僅在螺旋槽槽區局部溫度比模型1稍高。

圖12 轉速為3000、13000 r·min-1時兩模型的溫度分布云圖

2.5 空化熱效應對空化區域的影響

由圖 13可知,兩種模型動環端面的空化泡平均體積分數均隨著轉速的增加而增加,說明空化熱效應并沒有改變因液膜流速增大而導致臺階擴散通道壓力下降和空化產生的規律。同時,圖13還表明,動環端面的空化程度與空化熱效應有關,隨著轉速的增大,考慮空化熱效應時的動環端面空化程度增幅更大。

圖13 動環端面平均氣泡體積分數

圖14為轉速為3000 r·min-1和13000 r·min-1時液膜厚度方向不同位置的空化區域分布云圖。由圖可知,兩模型的液膜空化區域均發生在槽臺階擴散低壓區,且隨著轉速的增大空化區域增大。但兩模型的空化程度及分布規律不同,模型1的空化程度沿軸向槽底方向逐漸加劇,而模型2正好相反,也就是模型1的槽底為空化程度最高部位,而模型2的動環表面槽區位置為空化程度最高位置,槽底空化程度非常低。由于空化導致氣泡增加,混合黏度和能量耗散降低,因此,上述兩模型空化區域的不同導致模型2僅在螺旋槽凹槽處溫度比模型1稍高一些。由前面綜合分析已知,圖9處模型2低壓區壓力在膜厚方向上基本沒有變化,但低壓區的溫度在槽深方向上是逐漸降低的,根據飽和汽化壓力和溫度的對應關系,溫度降低,飽和汽化壓力降低,但槽底附近的壓力大小及分布和動環面上基本相等,導致壓力并沒有降低到對應溫度的飽和汽化壓力,從而使空化區域減小。這些現象說明空化熱效應對空化程度及分布有顯著的影響,其影響結果與溫度分布、壓力分布等因素密切相關。

圖 14 轉速為3000和13000 r·min-1時兩模型的空化區域分布云圖

3 結 論

(1)在低轉速下,空化熱效應對密封性能的影響可以忽略,但在高轉速下空化熱效應使上游泵送機械密封的高壓區形成能力減弱,泵送量降低,開啟力降低。在高轉速工作條件下分析密封失效機理時,除了考慮黏溫效應之外,還要考慮空化熱效應的影響。

(2)空化熱效應對密封端面間液膜溫度分布的影響小,相比于僅考慮黏溫效應的模型,螺旋槽槽區局部溫度稍有提升。

(3)與僅考慮黏溫效應的模型相比,考慮空化熱效應后,得到的空化分布和程度明顯不同。僅考慮黏溫特性時,在液膜厚度方向空化區域大小基本相等,但靜環槽底的氣泡體積分數最大,空化程度最嚴重。而綜合考慮黏溫特性和空化熱效應時,動環端面空化程度最嚴重,由動環端面沿膜厚方向至靜環槽底,空化區域越來越小,槽底空化區域最小。

符 號 說 明

cp ——水的比熱容,J·kg-1·K-1 Fcond——冷凝系數 Fr——模型1的開啟力,N Fra——模型2的開啟力,N Fvap——蒸發系數 F1,F2——分別為模型1、模型2的單周期液膜開啟力,N G——相對開啟力增量,% N——單周期液膜網格數 n——相數 Re,Rc——分別為氣泡產生、潰滅源項 rB——氣泡半徑,mm T——水溫,℃ v——空泡相 vm——質量平均速度,m·s-1 α——氣相體積分數 αnuc——氣核體積分數 β——黏溫系數,取0.03℃-1 δ1,δ2——分別為模型1、模型2前后兩次的開啟力相對誤差,% η——溫度T時的黏度,kg·m-1·s-1 η0——溫度T0時的黏度,kg·m-1·s-1 λ——水的熱導率 μm——混合黏性系數 ρm——混合物密度,kg·m-3

References

[1] FINDLEY J A. Cavitation in mechanical face seals [J]. Journal of Lubrication Technology, 1968, 90 (2): 356-364.

[2] JAKOBSSON B, FLOBERG L. The finite journal bearing, considering vaporization [J]. Transactions of Chalmers University of Technology, 1957, 190: 1-116.

[3] OLSSON K O. Cavitation in dynamically loaded bearings [J]. Transaction of Chalmers University of Technology, 1965, 308: 1-60.

[4] ELROD H G. A cavitation algorithm [J]. Journal of Lubrication Technology,1981, 103 (3): 350-354.

[5] LEBECK A O. Experiments and modeling of zero leakage backward pumping mechanical face seals [J].Tribology Transactions,2008, 51 (4): 389-395.

[6] DJAMA? A, BRUNETIèRE N, Tournerie B. Numerical modeling of thermohydrodynamic mechanical face seals [J]. Tribology Transactions, 2010, 53 (3): 414-425.

[7] QIU Y, KHONSARI M M. Performance analysis of full-film textured surfaces with consideration of roughness effects [J]. Journal of Tribology, 2011, 133 (2): 021704(1-10).

[8] MENG X K, BAI S X, PENG X D. Lubrication film flow control by oriented dimples for liquid lubricated mechanical seals [J]. Tribology International, 2014, 77: 132-141.

[9] 唐飛翔, 孟祥鎧, 李紀云, 等. 基于質量守恒的LaserFace液體潤滑機械密封數值分析 [J]. 化工學報, 2013, 64 (10): 3694-3700. TANG F X, MENG X K, LI J Y,. Numerical analysis of LaserFace liquid mechanical seal based on mass conservation [J]. CIESC Journal, 2013, 64 (10): 3694-3700.

[10] PASCOVICI M D, ETSION I. A thermo-hydrodynamic analysis of a mechanical face seal [J]. Journal of Tribology, 1992, 114: 639-645.

[11] LEBECK A O. Principles and Design of Mechanical Face Seals [M]. New York: Wiley-Interscience Publication, 1991.

[12] BRUNETIèRE N, MODOLO B. Heat transfer in a mechanical face seal [J]. International Journal of Thermal Sciences, 2009, 48 (4): 781-794.

[13] 中國科學院. 科學數據庫——工程化學數據庫[DB/OL]. http://www.sdb.ac.cn,2004. Chinese Academy of Sciences. Scientific Database—Chemical Engineering Database [DB/OL]. http://www.sdb.ac.cn, 2004.

[14] BRENNEN C E. Cavitation and Bubble Dynamics [M]. Oxford: Oxford University Press, 1995.

[15] 溫詩鑄, 黃平. 摩擦學原理[M]. 4版. 北京: 清華大學出版社, 2012: 9. WEN S Z, HUANG P. Principles of Tribology [M]. 4th ed. Beijing: Tsinghua University Press, 2012: 9.

[16] CHEN H L, XU C, ZUO M Z,. The thermal and mechanical deformation study of up-stream pumping mechanical seal [J]. IOP Conference Series: Materials Science and Engineering, 2015, 72: 042032.

[17] 李京浩. 機械密封空化效應的數值計算方法與實驗研究[D]. 北京: 清華大學, 2011. LI J H. Numerical computing method and experimental study for cavitation in mechanical seals [D]. Beijing: Tsinghua University, 2011.

[18] QIU Y, KHONSARI M M. Thermohydrodynamic analysis of spiral groove mechanical face seal for liquid applications [J]. Journal of Tribology, 2012, 134 (2): 021703 (1-11).

Influence of cavitation thermal effect on lubrication properties of upstream pumping mechanical seal

CHEN Huilong, WANG Bin, REN Kunteng, LI Tong, ZHAO Binjuan

(School of Energy and Power Engineering, Jiangsu University, Zhenjiang 212013, Jiangsu, China)

Cavitation occurred at mechanical seal faces is an important factor affecting lubrication properties of the mechanical seal. A computational fluid dynamics model was established from Antoine equation with a consideration of cavitation thermal effect. The cavitation thermal effect on the sealing performance was analyzed and compared to results of commonly used model of viscosity-temperature effect of liquid film on seal faces. The results indicated that the influence of cavitation thermal effect was negligible at low rotating speed whereas weakened the capacity of forming high-pressure region in mechanical seal at high rotating speed, which reduced pumping rate and opening force. Both viscosity-temperature effect and cavitation thermal effect were needed to analyze seal failure mechanism at high rotating speed. The local temperature at the spiral groove was slightly higher from cavitation thermal effect than that from viscosity-temperature effect. With the cavitation thermal effect, degree of cavitation occurred most seriously in the rotating ring face and cavitation space became smaller with the smallest one in groove bottom from the rotating ring face to the stationary ring groove bottom, which was contrary to those by considering viscosity-temperature effect.

mechanical seal; hydrodynamic lubrication; cavitation; thermodynamic properties; computational fluid dynamics; model; stability

2016-04-08.

Prof.CHEN Huilong, huji@ujs.edu.cn

10.11949/j.issn.0438-1157.20160456

TH 117.2

A

0438—1157(2016)10—4334—10

國家自然科學基金項目(51279067)。

2016-04-08收到初稿,2016-07-12收到修改稿。

聯系人及第一作者:陳匯龍(1961—),男,博士,教授。

supported by the National Natural Science Foundation of China (51279067).

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 婷婷亚洲视频| 日韩久草视频| 啪啪永久免费av| 日韩A∨精品日韩精品无码| 国产一级无码不卡视频| 亚洲国产成人久久77| 伊人久久大香线蕉成人综合网| 国内精自视频品线一二区| 五月婷婷丁香色| 精品五夜婷香蕉国产线看观看| 亚洲无码37.| 国产视频 第一页| 中文字幕人成人乱码亚洲电影| 五月天婷婷网亚洲综合在线| 四虎精品黑人视频| 国产成人啪视频一区二区三区 | 日韩无码黄色| 国产精品自在拍首页视频8| 亚洲永久精品ww47国产| 国产精品视频导航| 91成人免费观看| 新SSS无码手机在线观看| 国产精品网拍在线| 香蕉久久国产超碰青草| 欧美成人怡春院在线激情| 中文字幕首页系列人妻| 婷婷色丁香综合激情| 这里只有精品在线| 无码国内精品人妻少妇蜜桃视频| 狠狠色噜噜狠狠狠狠奇米777| 欧美视频在线第一页| 国产精品99一区不卡| 成人在线不卡| 在线五月婷婷| 国产网站免费| 日韩国产综合精选| 亚洲国产成人麻豆精品| 国产视频资源在线观看| 欧美成人精品欧美一级乱黄| 久久99热这里只有精品免费看| 日本91视频| 色综合综合网| 成年人国产网站| 国产第一页屁屁影院| 成人字幕网视频在线观看| 精品国产成人高清在线| 尤物精品视频一区二区三区 | 免费观看男人免费桶女人视频| 九月婷婷亚洲综合在线| 亚洲成a人在线播放www| 国产精品免费p区| 欧美一级99在线观看国产| 国产精品内射视频| a天堂视频| 中文字幕色站| 国产精品成人免费视频99| 久久这里只有精品8| 亚洲AⅤ综合在线欧美一区| 粉嫩国产白浆在线观看| 无码网站免费观看| 精品99在线观看| 国产欧美亚洲精品第3页在线| 日韩精品一区二区三区免费在线观看| 91毛片网| 国产va视频| 亚洲天堂视频在线免费观看| 白浆免费视频国产精品视频| 欧美精品成人一区二区在线观看| 国产亚洲精品资源在线26u| 91免费国产在线观看尤物| 久久国产精品麻豆系列| 蜜桃视频一区二区| 久久成人18免费| 蜜臀AV在线播放| 人妻无码中文字幕第一区| 日韩无码黄色网站| lhav亚洲精品| 51国产偷自视频区视频手机观看| 在线人成精品免费视频| 四虎精品国产AV二区| 亚洲国产欧美国产综合久久| 久草美女视频|