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

基于E/CRC磨損模型的離心泵壁面磨損特性研究

2021-05-07 09:55:08賴芬王鳳鳴朱相源常沛然李國君
哈爾濱工程大學學報 2021年5期
關鍵詞:模型

賴芬, 王鳳鳴, 朱相源, 常沛然, 李國君

(西安交通大學 熱流科學與工程教育部重點實驗室,陜西 西安 710049)

泵是把機械能轉化為流體壓力能的通用機械,其中離心泵使用最廣泛,特別是在水利水電工程中。但環境惡化引起的泥沙流失導致河流中含有泥沙,例如,中國的黃河泥沙含量為37.8 kg/m3[1],美國的科羅拉多河泥沙含量為10 kg/m3[2]。水利水電工程中的離心泵一般按清水工況設計,但其內部的實際流動為液固兩相流動。固體顆粒的存在不僅會降低離心泵的效率,還會對離心泵壁面造成侵蝕磨損。磨損將導致頻繁地維修和更換過流部件,造成輸運系統停機,輸運費用大幅度增加。Wilson等[3]指出大型礦場每小時的停機成本大約為10萬美元。因此,研究黃河和長江流域使用的含沙水離心泵內的液固兩相流動及其磨損特性具有重要意義。

但固體顆粒引起的壁面磨損是個復雜過程,最終的磨損形態受很多因素的影響,包括顆粒特性、流場特性、壁面特性等。國內外學者對壁面磨損的影響因素進行了大量的試驗研究。例如,Tian等[4]通過科里奧利磨損試驗臺測試了不同顆粒粒徑下壁面的磨損系數,結果表明顆粒形狀、粒徑分布等顆粒特性會對磨損系數產生重要的影響。陳紅生等[5]通過離心泵液固兩相流水力試驗,發現造成局部磨損的重要原因是葉輪出口附近的射流-尾流結構。Wiedenroth[6]采用超聲設備測試了4種葉輪葉片的磨損量,發現當顆粒撞擊角較大時,硬度高的壁面受到的磨損較嚴重。

隨著計算機技術的迅猛發展,數值模擬方法已成為研究各種物理現象的重要手段。目前關于離心泵液固兩相流的數值模擬及磨損預測,磨損模型主要選用Finnie[7]和Tabakoff[8]模型,研究對象主要針對離心泵的某一過流部件[1,9-11]。例如,Noon等[9]應用Finnie磨損模型預測了石灰漿輸送泵蝸殼的磨損形態,結果表明蝸殼的磨損率隨撞擊速度、質量濃度、顆粒粒徑的增大而增大。劉娟等[10]應用Finnie磨損模型對低體積分數的離散相顆粒在離心泵中的運動規律及葉輪壁面的磨損特性進行了數值研究,發現液固相密度差距越大,固體顆粒的運動跟隨性越差,固體顆粒與過流表面發生碰撞的幾率增大,葉輪壁面的磨損強度增加。黃先北等[1]采用Tabakoff 磨損模型對不同泥沙及不同入口工況下離心泵葉輪壁面的磨損規律進行了探索,發現顆粒粒徑會顯著影響葉輪壁面的磨損形態和位置,顆粒在離心泵入口分布越均勻,壁面磨損越分散,磨損位置軸對稱性越明顯。何創新[11]采用Tabakoff 磨損模型對單級雙吸中開式離心泵葉輪壁面的磨損特性進行了分析,結果表明葉輪入口前的密封體導葉能有效地抑制葉片的磨損,合理的葉片型線設計可以顯著降低葉片壓力側的磨損。

然而,大部分磨損模型的建立和推導都是基于氣固兩相流試驗的。磨損模型能否準確預測液固兩相流下的壁面磨損特性需進一步驗證。Zhang等[12]研究表明侵蝕磨損研究中心(E/CRC)提出的磨損模型不僅適用于氣固兩相流動,而且在液固兩相流動中也能獲得準確的結果。Peng等[13]指出E/CRC磨損模型與其他磨損模型相比,在CFD軟件中更易執行,模型中考慮了顆粒硬度和形狀因素,得到的數值結果與試驗結果更接近。但現有研究中未見采用E/CRC磨損模型預測離心泵壁面磨損特性的研究。因此,本文針對黃河和長江流域使用的含沙水離心泵展開了研究,建立了基于E/CRC磨損模型的離心泵壁面磨損特性預測數值方法,對各個區域的磨損形態進行了預測并對比了不同區域的最大和平均磨損率,分析了磨損率變化規律并預測了最大磨損率發生位置,另外還討論了顆粒粒徑及濃度對離心泵葉輪磨損特性的影響。

1 液固兩相流及磨損預測數學模型

本文采用的流動介質為黃河、長江流域的含沙水,根據長江夏季泥沙濃度[14]及黃河多年平均泥沙濃度[15],泥沙濃度大約為32 kg/m3,體積分數小于3%,泵內流動為低濃度液固流。文獻[16]指出歐拉-拉格朗日方法的使用條件是離散相體積分數小于10%~12%。因此,本文符合歐拉-拉格朗日方法的使用條件,為了提高計算精確性,本文選用雙向耦合的歐拉-拉格朗日方法求解。計算時忽略顆粒與顆粒間的相互作用力,并假定液固相之間沒有質量和能量交換。計算時將液體視為連續相,液體流場通過在歐拉坐標系中求解雷諾時均方程獲得;將固體顆粒視為離散相,固體顆粒運動通過在拉格朗日坐標系中求解顆粒軌跡方程獲得。液體湍流脈動引起的顆粒擴散采用隨機游走模型。顆粒撞擊壁面后發生的動量變化由Grant 和Tabakoff(G&T)[8]提出的顆粒碰撞反彈模型模擬,磨損率由E/CRC 提出的磨損模型進行計算,詳細的數學模型如下。

1.1 連續相控制方程

1) 質量守恒方程:

(1)

2) 動量守恒方程:

(2)

(3)

式中:μt為湍流粘度;μt由標準k-ε模型確定[17-18];k為湍動能;δij為“Kronecker delta”符號。

Fi為動量交換源項分量,Fi從以下方程獲得:

(4)

1.2 離散相控制方程

文獻[16]指出顆粒相動量源項主要來自拖曳力、重力、虛擬質量力、壓力梯度力、熱泳力、布朗力、Saffman 升力和Magnus升力產生的動量。但熱泳力主要是由溫度梯度引起的,離心泵中溫度變化很小,因此,本文忽略了熱泳力;布朗力和Saffman 升力只有處理亞微觀粒子時才考慮,Magnus升力只有處理旋轉粒子時才考慮,本文所考慮的粒子不是亞微觀粒子,也不旋轉,因此,本文忽略布朗力、Saffman 升力和Magnus升力,只考慮拖曳力、重力、虛擬質量力和壓力梯度力產生的動量源項。

(5)

(6)

式中:dp為顆粒直徑;μ為液體的分子粘度;Cd為拖曳系數,定義為:

(7)

式中:a1、a2、a3為由Morsi等[19]給定的常數,Re為相對雷諾數,定義為:

(8)

Fv為用于加速顆粒周邊流體的虛擬質量力矢量,定義為:

(9)

Fp為液體中的壓力梯度力矢量,當固體顆粒通過高壓區時,它對顆粒軌跡有著重要的影響,定義為:

(10)

1.3 隨機游走模型

顆粒軌道方程中的流體速度采用瞬時速度來考慮顆粒的湍流擴散,并計算足夠多的代表性顆粒的軌跡來考慮湍流對顆粒的隨機性影響。計算時考慮顆粒與流體的離散渦之間的相互作用,流體的脈動速度假定為時間的分段常量函數,在渦團的特征生存時間內脈動速度保持為常量,滿足高斯概率密度分布函數的隨機脈動速度u′、v′、w′為:

(11)

(12)

(13)

(14)

式中ζ為服從正態分布的隨機數。

渦團的特征生存時間定義為常量:

(15)

式中:TL為流體的拉格朗日積分時間尺度;k為湍動能;ε為湍動能耗散率。

1.4 G&T顆粒碰撞反彈模型

Peng等[13]對比了2種顆粒碰撞反彈模型和5種磨損模型預測的結果,結果表明G&T顆粒碰撞反彈模型結合E/CRC磨損模型預測的結果與試驗結果最接近,因此,本文選用G&T顆粒碰撞反彈模型。模型中,顆粒撞擊壁面后,其垂直于壁面切線方向的動量變化率為法向恢復系數Vn2/Vn1,其平行于壁面切線方向的動量變化率為切向恢復系數Vt2/Vt1,分別定義為:

Vn2/Vn1=0.993-1.76β+1.56β2-0.49β3

(16)

Vt2/Vt1=0.988-1.66β+2.11β2-0.67β3

(17)

式中:Vn、Vt分別表示顆粒撞擊速度的垂直分量和切線分量, m/s;下標1、2分別表示撞擊前和撞擊后;β為顆粒撞擊角, rad。

1.5 E/CRC磨損模型

顆粒撞擊壁面后,對壁面造成的磨損與壁面材料、顆粒特性、撞擊角等因素相關。與其他磨損模型相比,E/CRC磨損模型考慮了顆粒硬度和形狀因素,得到的數值結果與試驗結果更接近。E/CRC磨損模型中磨損率計算方程為[12]:

(18)

(19)

式中:ER為磨損率(每單位面積的質量損失),kg/m2;C和n分別取值2.17×10-7,2.41;BH為布氏硬度;FS為顆粒形狀系數,對于球形顆粒取值為0.2;Vp為顆粒速度,m/s;β為顆粒撞擊角, rad;Ai值見表1。

表1 Ai值Table 1 Values of Ai

1.6 數學模型驗證

由于文獻中關于離心泵壁面磨損率的試驗數據較少,因此,本文采用直徑為50 mm、曲率半徑為76.9 mm的 90°彎管試驗對數學模型進行驗證。雖然離心泵涉及旋轉部件,整體結構與彎管不一致,運動條件也有差別,但由磨損模型可知,磨損率主要與顆粒形狀、硬度、速度、撞擊角和局部壁面特性有關,與是否涉及旋轉部件和設備整體結構無關。運動條件差別將導致磨損模型中的撞擊角和撞擊速度差別,也就是說磨損模型中考慮了運動條件的差別。因此,彎管試驗是可以驗證數值模擬方法的。

彎管壁面磨損率的數值結果與Zeng等[20]采用電化學方法測量獲得的試驗結果對比如圖1所示。由圖1可知,沿著彎管曲率角逆時針方向,磨損率逐漸增加。彎管磨損率數值結果與試驗結果吻合良好,最大磨損率與最大磨損位置預測準確。因此,該數學模型可以準確地預測液固兩相流中壁面的磨損率與最大磨損位置。

圖1 彎管磨損率數值結果與試驗結果對比Fig.1 Comparison of elbow erosion rate between numerical results and experimental data

2 模型泵及數值設定

2.1 模型泵及數值計算域

本文研究對象為單級單吸離心泵,主要幾何參數如表2所示,其設計參數為:流量25 m3/h,揚程15 m,轉速2 500 r/min,比轉速135。數值計算域包括進口延伸段、葉輪、無葉擴壓器、蝸殼及出口延伸段,其中進口延伸段長度為葉輪進口直徑的3倍,出口延伸段長度為蝸殼出口直徑的3倍。

表2 模型泵主要幾何參數Table 2 Main geometric parameters of model pump

2.2 網格劃分

為了生成高質量的網格,采用ICEM-CFD里的結構化六面體網格對數值計算域進行劃分。為了準確地捕捉近壁面湍流,在近壁面采用加密處理,并調整第1層網格高度,將壁面y+值控制在30附近,最終的網格劃分結果如圖2所示。

圖2 計算域網格劃分Fig.2 Grid generation of computational domain

2.3 數值設定

數值計算時選取25°清水為連續相,球形沙顆粒(SiO2)為離散相。采用滑移網格技術模擬轉子與定子間的相對運動,設置進口延伸段和擴散段與葉輪的交界面為滑移界面,葉輪計算域設在旋轉坐標系,其余計算域設在靜止坐標系。進口邊界條件設為流體速度,出口邊界條件設為自由出流。球形沙顆粒從進口處垂直于進口面恒定釋放,速度與流體速度一致。非穩態計算時間步長設為0.000 1 s,即葉輪每旋轉1.5°求解一次,每個時間步長迭代次數設為20次。

2.4 網格無關性驗證

為了驗證網格無關性,選取7組結構化六面體網格對數值計算域進行劃分,并計算模型泵在設計工況下壁面的平均磨損率。平均磨損率隨網格數變化趨勢如圖3所示。由圖3可知,隨著網格數增加,平均磨損率呈減小趨勢,但當網格數增加至286萬后,平均磨損率的變化很小,在5%范圍內,因此,本文選取模型泵的計算網格數為286萬。

圖3 網格無關性研究設計工況Fig.3 Grid independence study

3 計算結果與分析

為了探索離心泵磨損初期各個區域的磨損形態和磨損率以及其受顆粒粒徑及濃度的影響,本文對設計工況、不同顆粒粒徑工況、不同顆粒濃度工況下的模型泵進行了研究。由于磨損初期壁面磨損較輕,壁面變化引起的通流型線變化較小,對磨損發生的情況和磨損條件的影響較小,因此,本文忽略通流型線的變化對磨損的影響。

文獻[2]指出黃河砂中值粒徑為30 μm,多砂河常見粒徑為90 μm。雖然實際泵流動的兩相流顆粒的直徑不會是某一均勻直徑,但采用平均粒徑計算的結果與采用一定顆粒直徑范圍計算的結果差別很小。為了分析粒徑對磨損形態和磨損率的影響,本文設計工況下選取平均粒徑60 μm,選取顆粒濃度32 kg/m3;不同顆粒粒徑工況時保持其他參數不變,顆粒直徑設為30、90、120 μm;不同顆粒濃度工況時保持其他參數不變,顆粒濃度設為19.5、44.5、57 kg/m3;具體計算工況如表3所示。

表3 計算工況Table 3 Calculation conditions

3.1 設計工況下離心泵壁面的磨損特性

為了探索離心泵壁面的磨損特性,對設計工況下模型泵各個區域的磨損形態進行了預測并對比了不同區域的最大和平均磨損率,分析了葉片與后蓋板交界處磨損率隨葉片曲率角的變化規律并預測了最大磨損率發生位置。

圖4為模型泵各個區域磨損形態的時間演化。由圖4可知,與其他過流部件相比,蝸殼和出口延伸段的磨損較輕,葉輪的磨損較嚴重,尤其是葉片前緣及葉片壓力側尾緣附近區域。與葉輪前蓋板相比,葉輪后蓋板磨損區域面積更大。當磨損時間由0.5 s增至2.0 s,蝸殼及出口延伸段的磨損率變化較小,葉輪及擴壓段的磨損率及磨損區域面積增加明顯,尤其是葉片壓力側尾緣附近區域。葉片前緣的磨損區域面積增加不明顯。當磨損時間為2 s時,前蓋板靠近葉片壓力側尾緣區域、后蓋板中部區域、后蓋板外圍區域均遭受了嚴重的磨損。

圖5為不同區域的最大和平均磨損率對比圖。由圖5可知,最大磨損率最大值發生在葉片和后蓋板區域,為7.9×10-4kg/m2;平均磨損率最大值發生在后蓋板區域,為1.4×10-6kg/m2;最大和平均磨損率最小值均發生在進口延伸段,分別為1.0×10-8和 2.6×10-9kg/m2;最大和平均磨損率分布規律不一致,例如,前蓋板的最大磨損率低于葉片,但平均磨損率高于葉片;同一區域的最大和平均磨損率差距巨大,例如,葉片的最大磨損率是平均磨損率的2 079倍。由以上分析可以推測,葉輪后蓋板是磨損最嚴重的區域,最大磨損率發生在葉片與后蓋板交界處。

圖6為葉片與后蓋板交界處磨損率隨葉片曲率角的變化規律。由圖6可知,葉片壓力側與吸力側的磨損率隨葉片曲率角的變化規律不同,葉片吸力側的磨損率曲線較平緩,葉片壓力側的磨損率曲線波動較大,葉片前緣吸力側的磨損率大于壓力側,葉片尾緣壓力側的磨損率大于吸力側;但無論是葉片壓力側還是吸力側,磨損率均是在葉片前緣取得最大值,最大磨損率發生在葉片曲率角59.8°附近,為7.9×10-4kg/m2。

3.2 顆粒粒徑對葉輪磨損特性的影響

為了進一步探索離心泵壁面磨損特性的影響因素,對不同顆粒粒徑工況下(工況1、2、3)離心泵的磨損率進行了計算。由上述分析可知,葉輪是離心泵內磨損最嚴重的過流部件,最大磨損率發生在葉片與后蓋板交界處。因此,這里分析顆粒粒徑對葉輪磨損形態及磨損率變化規律的影響,結果如圖7、圖8所示。圖7為不同顆粒粒徑工況葉輪磨損云圖變化,其中左邊為前蓋板,中間為葉片,右邊為后蓋板。由圖7可知,葉輪磨損形態受顆粒粒徑影響顯著。隨著顆粒粒徑的增大,磨損區域面積逐漸減小。這主要是由于在保持顆粒濃度不變的情況下,隨著顆粒粒徑的增大,顆粒數減小,葉輪壁面受到的撞擊次數減小。當顆粒粒徑為30 μm 時,葉輪前后蓋板中部及外圍均受到了不同程度的磨損,葉片前緣及葉片壓力側尾緣附近區域磨損嚴重;當顆粒粒徑增大至90 μm 時,前蓋板及葉片的磨損區域面積顯著減小,而后蓋板的磨損區域面積幾乎不變,但同一位置的磨損率顯著減??;當顆粒粒徑增大至120 μm 時,前蓋板的磨損區域只有葉片壓力側尾緣附近區域,后蓋板的磨損區域仍然是中部及外圍。

圖8為不同顆粒粒徑工況葉片與后蓋板交界處磨損率隨葉片曲率角的變化規律。其中,圖8(a)為葉片壓力側,葉片曲率角在-22.5°~59.8 °變化;圖8(b)為葉片吸力側,葉片曲率角在-32.5°~59.8°變化。由圖8(a)可知,在各個粒徑工況下,葉片壓力側磨損率基本是隨葉片曲率角先減小后增大,最大峰值在葉片曲率角59.8°附近,第2峰值在葉片曲率角-22.5°附近;由圖8(b)可知,在各個粒徑工況下,葉片吸力側磨損率基本是隨葉片曲率角先小幅度變化到葉片曲率角55°附近迅速上升,在葉片曲率角59.8°附近達到最大值;隨著顆粒粒徑增大,葉片壓力側及吸力側磨損率總體上呈減小趨勢。

圖4 磨損形態的時間演化Fig.4 Time evolution of erosion pattern

圖5 設計工況t=10 s時不同區域的最大和平均磨損率對比Fig.5 Comparisons of the maximum and average erosion rates for different regions under design condition at t=10 s

3.3 顆粒濃度對葉輪磨損特性的影響

為了進一步探索離心泵壁面磨損特性的影響因素,對不同顆粒濃度工況下(工況4、5、6)離心泵的磨損率進行了計算并分析顆粒濃度對葉輪磨損形態及磨損率變化規律的影響,結果如圖9、圖10所示。圖9為不同顆粒濃度工況葉輪磨損云圖變化,其中左邊為前蓋板,中間為葉片,右邊為后蓋板。由圖9可知,葉輪磨損形態受顆粒濃度影響顯著。隨著顆粒濃度的增大,同一位置的磨損率顯著增大。這主要是由于在保持顆粒粒徑不變的情況下,隨著顆粒濃度的增大,顆粒數增大,葉輪壁面受到的撞擊次數增多。但葉輪各個部分的磨損區域面積幾乎不變,這說明改變顆粒濃度,不會改變顆粒撞擊角。當顆粒濃度為19.5 kg/m3時,嚴重磨損主要集中在葉片前緣附近區域;當顆粒濃度增大至44.5 kg/m3時,同一位置的磨損率增大,葉片壓力側尾緣附近變成嚴重磨損區域;當顆粒濃度增大至57 kg/m3時,同一位置的磨損率進一步增大,葉片壓力側尾緣附近嚴重磨損區域面積增加明顯,葉片前緣附近嚴重磨損區域面積變化較小。

圖6 磨損率隨葉片曲率角的變化規律(dp=60 μm,Cm=32 kg/m3)Fig.6 Variation of the erosion rate along the blade curvature angle(dp=60 μm,Cm=32 kg/m3)

圖7 不同顆粒粒徑葉輪磨損云圖Fig.7 Impeller erosion contours for different particle diameters

圖10為不同顆粒濃度工況葉片與后蓋板交界處磨損率隨葉片曲率角的變化規律。其中,圖10(a)為葉片壓力側,葉片曲率角在-22.5°~59.8°變化;圖10(b)為葉片吸力側,葉片曲率角在-32.5°~59.8°變化。由圖10(a)可知,在各個濃度工況下,葉片壓力側磨損率的分布基本上是葉片前緣尾緣大,中間區域小,在葉片曲率角59.8°附近達到最大值,在葉片曲率角-22.5°附近達到第2峰值;由圖10(b)可知,在各個濃度工況下,葉片吸力側磨損率的分布基本上是葉片前緣大其余區域小,也是在葉片曲率角59.8°附近達到最大值;隨著顆粒濃度增大,葉片壓力側及吸力側磨損率隨葉片曲率角的變化規律較一致,同一位置上的磨損率基本上呈增大趨勢。

圖8 不同顆粒粒徑下磨損率隨葉片曲率角的變化Fig.8 Variations of the erosion rates along the blade curvature angle for different particle diameters

圖9 不同顆粒濃度葉輪磨損云圖Fig.9 Impeller erosion contours for different particle concentrations

圖10 不同顆粒濃度下磨損率隨葉片曲率角的變化Fig.10 Variations of the erosion rates along the blade curvature angle for different particle concentrations

4 結論

1) 葉輪是離心泵內磨損最嚴重的過流部件,最嚴重的磨損發生在葉片前緣及葉片壓力側尾緣附近區域。

2) 進口延伸段的最大和平均磨損率小于其它區域,葉片和后蓋板的最大磨損率大于其它區域,葉片與后蓋板交界處葉片曲率角59.8°附近的最大磨損率達到最大值。

3) 前蓋板及葉片的磨損區域面積隨粒徑增大顯著減小,而后蓋板的磨損區域面積隨粒徑增大幾乎不變,同一位置的磨損率隨粒徑增大呈減小趨勢。

4) 葉輪各個部分的磨損區域面積隨濃度增大幾乎不變,同一位置的磨損率隨濃度增大呈增大趨勢,葉片壓力側及吸力側磨損率隨濃度增大變化規律較一致。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 丰满人妻久久中文字幕| 色哟哟精品无码网站在线播放视频| 四虎成人精品在永久免费| 91色国产在线| 日本一区高清| 日本道综合一本久久久88| 萌白酱国产一区二区| 国产精品55夜色66夜色| 另类重口100页在线播放| 日本一区二区不卡视频| 欧美精品v欧洲精品| 在线观看91精品国产剧情免费| 中文字幕无码制服中字| 精品亚洲国产成人AV| 无码人妻热线精品视频| 她的性爱视频| 色悠久久久久久久综合网伊人| 亚欧美国产综合| 女人毛片a级大学毛片免费| 欧美国产日韩在线| 一区二区三区成人| 欧美精品高清| 九九精品在线观看| 久久综合伊人77777| 亚洲色欲色欲www网| 91精品伊人久久大香线蕉| 国产十八禁在线观看免费| 国产精品第三页在线看| 亚洲区第一页| 久青草网站| 亚洲中文制服丝袜欧美精品| 亚洲清纯自偷自拍另类专区| 欧美日韩高清| 无码一区二区三区视频在线播放| 992Tv视频国产精品| 亚洲国产第一区二区香蕉| 国产91无毒不卡在线观看| 欧美激情视频在线观看一区| 欧美日韩动态图| 欧美第一页在线| 国产菊爆视频在线观看| 成人在线不卡视频| 亚洲侵犯无码网址在线观看| 99在线国产| 欧美国产日产一区二区| 伊人久久福利中文字幕| 91精品福利自产拍在线观看| 丁香婷婷综合激情| 久久激情影院| 国产亚洲欧美另类一区二区| 亚洲乱码在线播放| …亚洲 欧洲 另类 春色| 美女潮喷出白浆在线观看视频| 欧美午夜理伦三级在线观看| 怡红院美国分院一区二区| 中文字幕av一区二区三区欲色| 日日噜噜夜夜狠狠视频| 女人av社区男人的天堂| 日韩在线中文| 黄色国产在线| 久久精品免费国产大片| 国产区福利小视频在线观看尤物| 国产精品v欧美| 国产精品xxx| 国产成人综合在线视频| 成人字幕网视频在线观看| 免费在线成人网| 精品色综合| 亚洲妓女综合网995久久| 国产波多野结衣中文在线播放 | 看国产毛片| 日韩中文字幕亚洲无线码| 精品一区二区无码av| 亚洲水蜜桃久久综合网站| 欧美怡红院视频一区二区三区| 91精品小视频| 亚洲Av激情网五月天| 国产丝袜第一页| yjizz国产在线视频网| 国产精品亚洲а∨天堂免下载| 国产乱子伦手机在线| a毛片在线|