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

級聯式氣-液旋流分離器流動特性數值研究

2024-01-20 10:40:24孫治謙孫銘澤王振波
石油學報(石油加工) 2024年1期

耿 坤,孫治謙,李 騰,孫銘澤,王振波

(中國石油大學(華東) 新能源學院,山東 青島 266580)

隨著大部分天然氣井達到高含水開采階段,天然氣品質受含水率的影響顯著增加[1],實現高效的除液過程能夠大大降低天然氣生產和運輸成本[2]。氣-液旋流分離器用于天然氣除液時,具有系統簡單、運行成本低、分離效果好等優勢[3-4]。盡管憑借諸多優勢,氣-液旋流分離器在天然氣除液領域廣泛應用,但仍需要緊湊的分離器結構來提高分離設備的空間利用率[5-6]。

前人在開發新型旋流分離器結構方面已做了大量研究[7-8]。Gao等[9]設計了一種用于井下油-水分離的級聯式水力旋流器,兩級入口均采用螺旋線式,但應用中仍存在易堵塞等缺陷。鄭春峰等[10]開發了一種用于井下的三級氣-液旋流分離器,研究填補了級聯式旋流分離器結構在氣-液分離領域的空白,為不同流型下氣-液混合物分離提供了有效方法,但該裝置內部結構復雜,加工制造困難。由此可見,國內外對于氣-液分離領域的級聯式旋流分離器相關研究成果較少,且目前開發的級聯式氣-液旋流分離器結構復雜,應用存在一定的局限性[11-13]。本研究中基于傳統單級雙蝸殼式氣-液旋流分離器結構開發的級聯式氣-液旋流分離器[14-15],在能夠實現兩級分離過程的同時,兼具結構簡單、便于加工等優勢。

在級聯式氣-液旋流分離器工程應用前需要了解其內流場分布特性[16-18]。然而,由于氣-液旋流分離器內為多相強旋流場,受湍流問題發展的限制,目前對于氣-液旋流分離器內部流動機理及細節尚無法精確描述,因此大量學者選擇通過計算流體動力學(CFD)的方法對旋流分離器內流動特性進行預測[19-20]。旋流分離器內的旋轉流是產生分離作用的主要載體[21],局部渦流和短路流是影響分離性能的重要因素[22],尤其是短路流,通常會對旋流分離器性能產生負面影響[23-24]。Dong等[25]通過空氣年齡模型(Air age model)計算短路流區域;Zhao等[26]研究了溢流管壁厚對短路流及軸向速度波帶的影響。盡管這些研究在一定程度上涉及了旋流分離器內流場的流動行為,但仍不能很好地識別各種流型的分布區域及位置信息,對短路流的量化也尚無統一的方法。

本研究基于旋流分離器常用的雷諾應力模型(RSM)進行流動仿真[27-28],詳細分析了級聯式氣-液旋流分離器上下行流和內外旋流的分布尺度,探究了局部渦流與短路流的產生位置及特征,提出了一種短路流流量的理論計算方法,為級聯式氣-液旋流分離器流場特性分析提供了新思路。

1 級聯式氣-液旋流分離器幾何模型及實驗裝置

級聯式氣-液旋流分離器結構如圖1所示。入口為對稱布置的雙蝸殼結構[29],一級旋流管上方增設了二級旋流管。級聯式氣-液旋流分離器筒體直徑(Dh)為0.10 m,詳細結構參數尺寸如表1所示。

圖1 級聯式氣-液旋流分離器結構示意圖Fig.1 Structure diagram of cascade gas-liquid cyclone separator

表1 級聯式氣-液旋流分離器幾何尺寸Table 1 Geometry size of cascade gas-liquid cyclone separator

圖1右圖為級聯式氣-液旋流分離器工作過程示意圖,氣、液兩相由雙蝸殼入口進入分離器后旋轉向下運動,大部分液滴在離心力的作用下運動到一級旋流管內壁聚集形成液膜,沿壁面向下流動進入錐段,匯集至液相出口處。氣、液兩相流到達旋流管底部后折返向上,進入排氣芯管,在排氣芯管頂部擋板及兩側開口的作用下進入二級旋流管的環形空間,此時,一級旋流管內未實現分離的細小液滴以及隨短路流逃逸至二級旋流管的部分液滴,在環形空間中進行第二次離心分離,進一步分離后的液相聚集至二級旋流管底部,沿降液管向下流進一級旋流管內,最終匯集至液相出口,而氣相則在壓力作用下向上運動,由頂部氣相出口排出。

為測試級聯式氣-液旋流分離器的分離性能,驗證模擬結果的準確性,設計了如圖2所示的實驗流程。其中,霧化系統由空氣壓縮機、旋渦泵、水箱和扇型霧化噴嘴組成;引風系統由旋渦氣泵組成;分離系統由緩沖腔、級聯式氣-液旋流分離器和水箱組成;測量系統由渦街流量計、畢托管、液體轉子流量計、U型管壓差計、溫度計和壓力表等組成。實驗物料采用空氣和水,常溫實驗條件下其物性與模擬物料一致。實驗表明該旋流器工作過程穩定,具有較高的研究和應用價值[30]。

1—Liquid tank;2—Vortex pump;3—Liquid pressure gauges;4—Rotor flowmeter;5—Air compressor;6—Air regulator filter valve; 7—Atomising nozzles;8—Buffer chamber;9—Cascade gas-liquid cyclone separator;10—U-shaped differential pressure gauge; 11—Thermometer;12—Pitot tube;13—Vortex flowmeter;14—Vortex air pump圖2 級聯式氣-液旋流分離器實驗裝置流程圖Fig.2 Experimental setup flow chart of cascade gas-liquid cyclone separator

2 級聯式氣-液旋流分離器數值模擬方法

2.1 數值模型控制方程

由于氣-液旋流分離器中流動表現出強湍流特性,本研究中湍流計算選用雷諾應力(RSM)模型,該模型對于復雜各向異性湍流的渦旋、旋轉及應變率變化等現象都能實現準確預測[27]。將空氣看作不可壓縮流動介質時,其雷諾平均應力(RANS)方程如式(1)和式(2)所示。

(1)

(2)

(3)

式(3)中,等式左側分別為應力的當地時間導數和對流運輸項,等式右側分別為:

在分離性能研究時采用離散顆粒模型(DPM)模擬液滴運動,運動方程符合牛頓第二定律,如式(4)所示。

(4)

式(4)中,等式右側分別為阻力項、重力與浮力的合力項。

與固體顆粒不同,液滴在運動過程中通常會伴隨著碰撞聚并和破碎現象,為準確計算該現象,本研究在DPM模型的基礎上增加了斯托克斯碰撞聚并模型和泰勒類比破碎(TAB)模型,依據振蕩液滴的臨界變形情況衡量液滴破碎過程,其控制方程可表示為式(5)。

(5)

式(5)中,Cb、Cd、CF和Ck為無量綱常數,由前人[31]大量實驗及理論推導可得取值分別為1/3、1/2、8和5。

2.2 數值模擬條件

本研究中空氣作為連續相,黏度為1.79×10-5Pa·s,密度為1.225 kg/m3;水作為離散相,黏度為0.001 Pa·s,密度為998.6 kg/m3。綜合考慮物料的性質,設置為速度入口,兩相均勻進料速度為1~7 m/s,入口氣量為63.36~443.52 m3/h,入口液量為78 mL/min,液滴直徑為0.1~10 μm;氣相出口兩相流動充分發展,設置為收斂性較好的自由出流,液滴到達此處逃逸;液相出口在實際工作時為液封狀態,無氣流溢出,可看作壁面。其他壁面均設置為無滑移的標準壁面函數,液滴到達壁面被捕捉。本研究在瞬時狀態下進行模擬,時間步長為0.001 s,通過壓力-速度耦合的SIMPLEC算法進行雙精度計算,壓力求解方法為PRESTO!,其他變量求解方法為QUICK,與Chen等[27]采用的求解方法相同。

2.3 Q準則

級聯式氣-液旋流分離器中存在著復雜的渦流現象,對于渦結構的分析也是流動特性研究的重要內容。本研究采用Q準則對渦結構進行識別,該識別方法應用廣泛,計算效率高[31-32]。

氣相流動過程中的速度梯度張量(如式(6)所示)可分解為兩部分:應變速率張量(S)和渦量張量(Ω)。

(6)

則Q值可定義為:

(7)

通過Q準則對渦結構進行識別時,可根據Q值大小較精確地識別出應變率或渦量占據主導的區域,當Q<0時,流場中黏性應力占據主導;當Q>0時,流場中渦量占據主導,可判斷為級聯式氣-液旋流分離器中的渦旋區域。

3 級聯式氣-液旋流分離器數值模擬準確性驗證

3.1 網格無關性驗證

合理的網格模型能夠有效降低模擬誤差。對級聯式氣-液旋流分離器流域進行多面體網格劃分,在結構狹縫處進行局部網格加密,以保證網格質量要求。圖3為級聯式氣-液旋流分離器網格模型。圖4為網格無關性驗證。在(8~30)×104間選取6種數量的網格模型進行計算,如圖4所示,網格數量達到123015后,級聯式氣-液旋流分離器內切向速度峰值基本不變;達到220135后,中心處切向速度值也無明顯變化。此外,計算結果表明,網格數量達到220135后,級聯式氣-液旋流分離器計算壓降保持在9300 Pa基本不變。由于數值計算效率會隨網格數量增加而降低,因此本研究選用的網格數量為220135,在高效計算的同時,能夠獲得準確可靠的流動特性及分離性能研究結果。

圖3 級聯式氣-液旋流分離器網格模型Fig.3 Grid model of cascade gas-liquid cyclone separator

圖4 網格無關性驗證Fig.4 Grid independence verification

3.2 數值模型驗證

為驗證數值模型的準確性,采用本研究的數值模型對文獻[33]中的旋流分離器結構進行模擬計算,計算所得的切向速度和軸向速度模擬指與文獻[33]中所得的切向速度和軸向速度實驗值吻合度較高(如圖5所示),因此本研究中所采用的RSM模型能夠準確地計算出旋流分離器的流動特性,這與Deng等[28]和Gao等[32]的研究結論相似。

圖5 RSM數值模型驗證Fig.5 RSM numerical model validation(a)Tangential velocity (vt);(b)Axial velocity (va)

為進一步驗證數值模型計算模擬結果的準確性,將本研究的數值模型模擬結果與實驗結果進行對比。實驗物料采用常溫空氣和水,入口氣量為63.36~443.52 m3/h,入口液量為78 mL/min,液滴直徑范圍為0.1~10.0 μm,數值模擬條件與實驗條件相同,對比結果如表2所示。由表2可以看出,實驗的與模擬的壓降及效率相對誤差均小于20%,進一步表明本研究對級聯式氣-液旋流分離器流動特性和分離性能的預測準確性較高。

4 級聯式氣-液旋流分離器流動分離特性

級聯式氣-液旋流分離器內流場分布情況是影響其分離性能的根本原因,以入口氣速為4 m/s、入口氣量為253.44 m3/h、液量為78 mL/min條件為例,深入探究級聯式氣-液旋流分離器流動特性。

4.1 旋轉流流動特性

級聯式氣-液旋流分離器內的旋轉流是產生離心力的主要流型,旋轉流呈邊壁為下行流、中心為上行流的分布特點。上、下行流的分布特點通常根據級聯式氣-液旋流分離器內氣流的軸向速度判斷,軸向速度為0的位置稱作零速包絡面。圖6(a)為一級旋流區內軸向速度分布曲線。在一級旋流區內,不同高度處的軸向速度分布特點相似,零速包絡面出現在R=0.55rh位置,該位置沿徑向向外為下行流,沿徑向向內為上行流。受級聯式氣-液旋流分離器邊壁與中心處壓差作用的影響,氣相在向下流動的過程中會不斷向中心匯集,因此越向下氣相流量越少,上、下行流的軸向速度均逐漸減小。圖6(b)為二級旋流區內軸向速度分布曲線。二級旋流區內軸向速度分布與一級旋流區存在較大差異,分離區底部(z=6.5Dh)受降液縫隙處溢出氣流的影響,軸向速度較大,其他位置軸向速度變化較小。

表2 級聯式氣-液旋流分離器模擬值與實驗值對比Table 2 Comparison between simulated and experimental values of cascade gas-liquid cyclone separator

圖6 不同z坐標位置下軸向速度(va)分布Fig.6 Axial velocity (va)distribution at different z-coordinates(a)Primary cyclone;(b)Secondary cyclone

沿排氣芯管向上,級聯式氣-液旋流分離器中心處軸向速度不斷降低,當旋流分離器中心線處軸向壓力梯度為正時,存在軸向速度停滯現象[34]。圖7為中心線處壓力及壓力梯度沿z坐標位置的變化。在一級旋流區內壓力梯度為負值,軸向速度不發生停滯,而排氣芯管內壓力梯度為正值,故存在明顯停滯現象。

另一種常見的描述旋流分離器流動的方法是從渦的角度進行的,級聯式氣-液旋流分離器內旋轉流場呈現蘭金組合渦的分布特點,該特點與旋流分離器內氣流的切向速度分布密切相關,前人研究[35-36]表明,旋流分離器內切向速度與徑向位置存在如式(8)

圖7 級聯式氣-液旋流分離器中心線處壓力梯度分布Fig.7 Pressure gradient distribution at centerline of cascade gas-liquid cyclone separator

所示的關系。

vt·Rn=C

(8)

當n>0時,vt與R呈負相關,流體內部幾乎無摩擦損失,為準自由渦區域;當n<0時,vt與R呈正相關,流體轉動近似剛體特性,為準強制渦區域。

圖8為級聯式氣-液旋流分離器內不同z坐標位置處切向速度分布曲線。圖8(a)為一級旋流區內切向速度分布曲線。一級旋流區內切向速度隨軸向高度的增加略有上升,不同z坐標位置處切向速度均呈M狀分布,切向速度峰值處(R=0.35rh)即為準強制渦與準自由渦分界處,n<0的準強制渦區域即為內旋流區或渦核區,n>0的準自由渦區域即為外旋流區。圖8(b)為二級旋流區內切向速度分布曲線。受流動阻力的影響,排氣芯管內切向速度隨高度的增加而不斷降低。由于二級旋流區分離空間受排氣芯管限制,氣流沿邊壁旋轉向下的同時伴隨著排氣芯管外壁附近的上行流運動,因此無明顯的組合渦分布特點。

圖8 不同z坐標位置下切向速度(vt)分布Fig.8 Tangential velocity (vt)distribution at different z-coordinates(a)Primary cyclone;(b)Secondary cyclone

從渦的角度分析旋轉流流動時可以借助Q準則進行渦結構識別,根據Q準則定義可知,Q值能夠度量某一局部位置處旋轉速率相對于應變速率的超額量[37],故可依據Q值為0的位置判斷蘭金組合渦分布范圍[32]。圖9為一級旋流區內不同z坐標位置處Q值分布曲線。圖9中由Q值判斷的準強制渦和準自由渦分界位置與圖8中由切向速度峰值判斷的準強制渦和準自由渦分界位置一致,均在R=0.35rh處。

圖9 不同z坐標位置下Q值分布Fig.9 Q-value distribution at different z-coordinates

4.2 短路流流動特性

4.2.1 局部渦流

級聯式氣-液旋流分離器內二次渦流的存在往往會對短路流的位置和流量變化產生較大影響,部分二次渦流會伴隨主流運動穩定存在,通過Q值分布結合流線分布可以較為準確地識別級聯式氣-液旋流分離器中局部渦流情況[38]。

圖10為級聯式氣-液旋流分離器內Q值分布云圖及局部渦流流線圖。圖10(a)為擋板上方尾跡渦流,2個小渦對稱分布于擋板邊緣兩側,距離擋板較近。圖10(b)為擋板下方氣流折返形成的局部渦流,渦尺寸較大。圖10(c)為二級旋流區底部環形渦流,該位置受降液縫隙溢出氣流的影響,產生了明顯的局部循環流動。圖10(d)為頂灰環,與傳統旋流分離器相似,級聯式氣-液旋流分離器在一級旋流區頂部也存在環形氣流死區,形成頂灰環區域,但與傳統旋流分離器相比,級聯式氣-液旋流分離器內頂灰環區域相對較小。圖10(e)為排氣芯管與降液管下方縱向渦流和短路流,該縱向渦流是由于“節流效應”產生的,受降液縫隙處短路氣流的劇烈影響,氣流還在降液管與排氣芯管間的位置產生了一個小尺寸的局部湍流渦。圖10(f)為底流口處尾跡渦流,是底流口處氣流向上折返流動產生的。

圖10 級聯式氣-液旋流分離器內局部渦流分布Fig.10 Local vortex distribution of cascade gas-liquid cyclone separator

4.2.2 短路流

在級聯式氣-液旋流分離器內的各種流動形態中,短路流會對旋流分離器的分離性能產生最直接的影響[25]。受局部縱向渦流(如圖10(e)所示)的影響,通過排氣芯管底部單一截面軸向速度積分確定短路流流量會存在較大的誤差[24],因此筆者提出了一種新的短路流流量計算方法。

級聯式氣-液旋流分離器中短路流流量由qs1和qs2兩部分組成(如圖11所示)。一級旋流區中外旋流在向下流動時會不斷有氣相在壓力的作用下匯入內旋流中,而發生短路流的位置則是突然有大量流體進入排氣芯管內匯入內旋流中,因此,排氣芯管區域內(R<0.5Do)氣相流量變化可反映短路流產生位置及短路流流量。以排氣芯管底部截面P2為參照面,P2下方存在某一截面P1,為發生短路流的邊界位置;P2上方則存在某一截面P3,為短路流完全匯入排氣芯管后的位置。P1和P3的截面流量差值即為從排氣芯管底部逃逸的短路流流量qs1;qs2為從降液管底部逃逸的短路流流量,可用降液縫隙任意截面P4的流量表示。qs1和qs2的計算式如式(9)和式(10)所示。

(9)

(10)

圖11 級聯式氣-液旋流分離器短路流示意圖Fig.11 Short-circuit flow of cascade gas-liquid cyclone separator

為確定短路流發生位置(即P1位置),計算了在不同入口氣速下向心流流量(qc)隨z坐標位置的變化,如圖12所示。由圖12可以看出:在一級旋流區底部向心流流量較為穩定,這證實了一級旋流區中氣相會不斷由外旋流向內旋流匯集;而在一級旋流區上部的位置,受縱向渦流的影響,向心流流量有所增加;而在縱向渦流上方,向心流流量突然增大,該突增位置即為短路流發生的邊界位置,即P1位于z=4Dh處;不同入口氣速下縱向渦流及短路流產生的位置基本相同。

圖12 向心流流量(qc)隨z坐標位置變化曲線Fig.12 Variation curve of centripetal flow rate (qc) with z-coordinate

為進一步確定截面P1、P2位置,本研究計算了R<0.5Do范圍內截面流量qP隨著z坐標位置的變化,以入口氣速為4 m/s時為例,如圖13所示。由圖13可以看出:截面流量突然增大的位置為P1所在位置,即P1位于z=4Dh處,與圖12中確定的P1位置相同,P1截面流量qP1=80.46 m3/h;截面流量達到穩定后的位置為P3所在位置,即P3位于z=4.3Dh處,P3截面流量qP3=97.63 m3/h,qP1與qP3之差即為從排氣芯管底部逃逸的短路流流量qs1=17.17 m3/h。按照式(8)計算P4截面流量,可得到從降液管底部逃逸的短路流流量qs2=155.91 m3/h。

圖13 截面流量(qP)隨z坐標位置變化曲線Fig.13 Variation curve of cross-sectional flow rate (qP)with z-coordinate

在計算qs1、qs2基礎上得到級聯式氣-液旋流分離器中短路流流量占入口流量的比值,如圖14所示。由圖14可以看出:當入口氣速由1 m/s增加至7 m/s時,從排氣芯管底部逃逸的短路流流量qs1由4.93 m3/h升高至33.63 m3/h,從降液管底部逃逸的短路流流量qs2由35.66 m3/h升高至278.10 m3/h,短路流流量占入口流量的比值也由64.07%增至70.29%。

圖14 短路流流量及其占比隨入口氣速(vin)的變化曲線Fig.14 Variation curves of short-circuit flow rate and percentage with inlet air velocity (vin)

4.3 分離性能

分離效率是衡量級聯式氣-液旋流分離器分離性能的直接指標,由氣相出口處液滴逃逸情況體現。圖15為不同入口氣速下分離效率隨液滴直徑的變化曲線。傳統旋流分離器分離效率通常會隨著短路流流量增加而有所降低[23]。而本研究中的級聯式氣-液旋流分離器分離效率隨入口氣速增加不斷增大,并未因短路流流量占比增大而出現明顯的降低現象。一方面是由于隨著入口氣速的增加,主導分離的離心力不斷增大;另一方面是由于級聯式氣-液旋流分離器中增設的二級旋流區能夠有效分離短路流攜帶的液滴,緩解“跑粗”現象對分離過程的負面影響。不同入口氣速下分離效率的變化規律具有相似性,液滴直徑小于0.2 μm時存在明顯的“魚鉤效應”,極小的液滴在流動過程中因更易發生團聚而使其分離效率略高。當液滴直徑大于0.2 μm時,分離效率隨液滴直徑的增加不斷提高,直徑6 μm以上的液滴分離效率可達100%。

圖15 級聯式氣-液旋流分離器分離效率(η)變化曲線Fig.15 Variation curves of separation efficiency (η)of cascade gas-liquid cyclone separator

5 結 論

(1)Q準則能夠準確識別級聯式氣-液旋流分離器內蘭金組合渦分布情況,Q值為0處即為準自由渦與準強制渦的分界位置。新型級聯式氣-液旋流分離器中一級旋流區內上下行流分界線位于R=0.55rh處,內、外旋流分界線位于R=0.35rh處;二級旋流區流動不呈現明顯的蘭金組合渦分布特點。

(2)根據匯入排氣芯管區域內氣相流量的變化能夠有效判斷短路流的發生位置及短路流流量,氣-液旋流分離器中短路流流量通常較高,當入口氣速為7 m/s時,級聯式氣-液旋流分離器中一級旋流區內短路流流量高達入口流量的70.28%。

(3)級聯式氣-液旋流分離器內存在大量的局部二次渦流,但二級旋流區能夠有效分離短路流攜帶的液滴,提高旋流器分離效率,兩級分離后級聯式氣-液旋流分離器對直徑6 μm以上的液滴能夠實現完全分離。

符號說明:

A——截面面積,m2;

a,b——排氣芯管開口寬度和高度,m;

CD——阻力系數;

Dg——氣相出口直徑,m;

Dh——筒體直徑,m;

Dl——液相出口直徑,m;

Do——排氣芯管直徑,m;

Dt——降液管直徑,m;

Dp——液滴直徑,μm;

ep——壓降相對誤差,%;

eη——分離效率相對誤差,%;

h——雙蝸殼入口高度,m;

ht——降液管長度,m;

L1——一級旋流管長度,m;

L2——二級旋流管長度,m;

n——切向速度指數;

P——截面;

p——靜壓,Pa;

p′——脈動壓力,Pa;

qc——向心流流量,m3/h;

qin——入口流量,m3/h;

qP——截面流量,m3/h;

qs1——排氣芯管底部逃逸的短路流流量,m3/h;

qs2——降液管底部逃逸的短路流流量,m3/h;

R——徑向位置,m;

rh——筒體半徑,m;

rp——初始液滴半徑,μm;

t——時間,s;

u——瞬時速度,m/s;

u′——脈動速度,m/s;

up——液滴運動速度,m/s;

va——軸向速度,m/s;

vin——入口氣速,m/s;

vt——切向速度,m/s;

w——網格數量;

x,y,z——三維坐標位置,m;

Δd——降液縫隙寬度,m;

Δp——壓降模擬值,Pa;

Δpexp——壓降實驗值,Pa;

η——分離效率模擬值,%;

ηexp——分離效率實驗值,%;

μ——氣相動力黏度,Pa·s;

μp——液滴動力黏度,Pa·s;

ρ——氣相密度,kg/m3;

ρp——液滴密度,kg/m3;

σ——液滴表面張力,N/m;

下標:

i,j,k——矢量方向。

主站蜘蛛池模板: 亚洲 日韩 激情 无码 中出| 欧美精品一二三区| 久热中文字幕在线| 久久特级毛片| 欧美另类视频一区二区三区| 欧美成a人片在线观看| 热久久综合这里只有精品电影| 亚洲高清日韩heyzo| 国产区91| 国产主播福利在线观看| 国产综合精品一区二区| 色综合五月婷婷| 成人伊人色一区二区三区| 国产玖玖视频| 免费午夜无码18禁无码影院| 國產尤物AV尤物在線觀看| 日韩成人免费网站| 免费毛片全部不收费的| 成人国产三级在线播放| 国产精品第一区| 成人av手机在线观看| 国产精品免费电影| 欧美在线视频不卡第一页| 国产成人高清在线精品| 国产综合网站| 日本一区二区三区精品国产| 4虎影视国产在线观看精品| 看av免费毛片手机播放| 永久免费无码成人网站| 国产丝袜第一页| 日韩精品亚洲一区中文字幕| 亚洲欧美精品日韩欧美| 日韩高清中文字幕| 久久久久亚洲av成人网人人软件 | 亚洲国产成人久久77| 丝袜亚洲综合| 国产精品99久久久久久董美香| 久996视频精品免费观看| 欧美不卡在线视频| 午夜日b视频| 亚洲自偷自拍另类小说| 伊人久综合| 亚洲日本在线免费观看| 欧美精品xx| 欧美在线精品一区二区三区| 成人国产三级在线播放| 动漫精品啪啪一区二区三区| 无码精品国产VA在线观看DVD| 国产精品主播| 欧美一级专区免费大片| 国产丝袜丝视频在线观看| av在线人妻熟妇| 国产一级毛片高清完整视频版| 青青草国产在线视频| 热九九精品| 久青草免费视频| 国产高清在线观看91精品| 国产在线观看99| 精品国产www| 久久亚洲中文字幕精品一区| 亚洲第一黄片大全| 中文字幕 91| 喷潮白浆直流在线播放| 久久精品午夜视频| 99久视频| 亚洲日韩精品综合在线一区二区| 真实国产乱子伦高清| 五月婷婷丁香综合| 少妇被粗大的猛烈进出免费视频| 亚洲综合色婷婷| 欧美三级不卡在线观看视频| 中文字幕日韩欧美| 亚洲美女久久| 亚洲欧洲日韩综合色天使| 欧美国产综合色视频| 国产乱子伦精品视频| 国产全黄a一级毛片| 理论片一区| 18禁影院亚洲专区| 亚洲综合极品香蕉久久网| 三区在线视频| 国产理论一区|