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

旋流噴嘴內超臨界流體閃蒸過程的數值模擬

2016-08-11 09:47:32馮留海王江云孫中衛
石油學報(石油加工) 2016年4期

馮留海, 王江云, 趙 凡,3, 孫中衛,4, 王 娟, 毛 羽

(1.中國石油大學 重質油國家重點實驗室, 北京 102249; 2.北京低碳清潔能源研究所, 北京 102209;3.蘭州蘭石能源裝備工程研究院有限公司,甘肅 蘭州 730314; 4.新奧科技發展有限公司, 河北 廊坊 065001)

?

旋流噴嘴內超臨界流體閃蒸過程的數值模擬

馮留海1,2, 王江云1, 趙凡1,3, 孫中衛1,4, 王娟1, 毛羽1

(1.中國石油大學 重質油國家重點實驗室, 北京 102249; 2.北京低碳清潔能源研究所, 北京 102209;3.蘭州蘭石能源裝備工程研究院有限公司,甘肅 蘭州 730314; 4.新奧科技發展有限公司, 河北 廊坊 065001)

摘要:旋流噴嘴內超臨界流體中瀝青溶質的體積分數分布對顆粒成形有重要影響。根據減壓相變傳質傳熱理論開發了閃蒸相變模型,采用自定義函數(UDF)的方式植入到CFD軟件Fluent中。將閃蒸相變模型耦合多相流混合模型用于研究旋流噴嘴內超臨界流體的閃蒸相變過程,分析旋流噴嘴內壓力、速度、溫度和各相濃度分布,以預測旋流噴嘴對顆粒成形的影響。結果表明,旋流噴嘴內三相介質分層流動,從而實現戊烷溶劑與瀝青溶質的預分離,有利于形成粒徑較小且密實的瀝青顆粒。

關鍵詞:旋流噴嘴; 數值模擬; 閃蒸相變; 非平衡熱力學

超臨界流體溶劑脫瀝青技術是重質油梯級分離工藝中的重要組成部分。利用超臨界戊烷溶劑可以選擇性去除渣油中的瀝青質、稠環化合物和重金屬等雜質,分離得到加工性能較好的脫瀝青油和高軟化點的脫油瀝青[1-2]。戊烷-瀝青超臨界流體從噴嘴內閃蒸噴出,并快速膨脹造粒[3],在噴霧造粒塔內實現重組分造粒與分離。閃蒸相變噴嘴是決定造粒質量的核心部件。研究超臨界瀝青造粒中噴嘴內戊烷的閃蒸相變及流動過程,對優化噴嘴結構、改進工藝流程具有十分重要的意義[4-5]。

閃蒸相變過程廣泛存在于工業生產中。對閃蒸相變過程的研究結果[6-10]表明,閃蒸過程受溫度、壓力、汽化潛熱等因素影響較大。人們大多從實驗的角度研究閃蒸相變的影響因素,但是由于閃蒸相變過程的復雜性和工程問題的多樣性,擬合得到的閃蒸相變經驗公式具有局限性,難以全面揭示戊烷閃蒸相變機理,并且本研究對象為高壓狀態下易燃易爆的戊烷-瀝青超臨界流體,又給實驗測量帶來了極大的困難。

近年來,計算流體力學為復雜流動過程的研究提供了新的途徑,數值模擬可以通過迭代計算控制方程組得到流場的詳細信息。建立相變模型研究相變過程日趨增多。基于小時間尺度、較低環境溫度下相變過程滿足熱力學平衡過程的假設,前人成功建立了各種空穴模型[11-12]。但是超臨界流體閃蒸相變過程熱傳輸傳遞的能量只有部分用于溶劑氣化,屬于典型的非平衡熱力學過程。聶永廣等[13]研究了不同噴嘴結構對戊烷氣化率和瀝青顆粒成形的影響。結果表明,在噴嘴入口段后加入漸縮段,能夠抑制戊烷溶劑氣化速率,有利于形成粒徑較小且密實的瀝青顆粒。為進一步提高瀝青密度,需要進行多次降壓,但降壓以后瀝青相在非超臨界戊烷溶劑中的溶解度會迅速降低并大量析出。所以,上述射流噴嘴會導致兩個問題,即過長的噴嘴會增大加工難度和瀝青大量析出很容易堵塞噴嘴,給生產帶來安全隱患。為避免上述不足,筆者所在課題組提出了一種可控粒徑瀝青殘渣噴霧造粒噴嘴,擬通過旋流過程預先分離戊烷溶劑和瀝青溶質,來提高造粒質量。由于對旋流噴嘴的機理研究仍不充分,筆者根據閃蒸相變原理,將閃蒸相變過程視為壓力突變產生的沸騰過程,從而將沸騰模型擴展為閃蒸相變模型,并與Fluent計算軟件中的多相流混合模型耦合,還考慮了輻射對溫度場的影響,模擬了氣、液戊烷與瀝青的三相閃蒸相變過程,以期為旋流噴嘴的設計提供理論依據。

1 旋流噴嘴幾何模型及網格劃分

本研究主要考察超臨界流體旋流噴嘴噴孔段流動分布和噴口處閃蒸相變過程。旋流噴嘴由入口段、旋流構件、噴孔段組成,結構如圖1所示。計算時,假設經過旋流構件進入噴孔段的超臨界流體達到穩定旋轉流動狀態,可以將三維旋流噴嘴結構簡化成二維模型,從而節省計算耗時。旋流噴嘴水平放置,以噴孔段入口圓心為坐標原點,采用Gambit建模軟件對旋流噴嘴和柱狀空間進行完全結構化網格劃分;在管壁附近加密網格,以考慮邊界層對模擬結果的影響,網格節點數為10421。旋流噴嘴閃蒸相變幾何模型如圖2所示,并通過Fluent軟件模擬其流場分布和閃蒸相變過程。

圖1 旋流噴嘴結構示意圖Fig.1 Schematic of swirl nozzle structure

圖2 旋流噴嘴幾何模型示意圖Fig.2 Schematic of geometric model for the swirl nozzleD1=5 mm; D2=2 mm; D3=150 mm;D4=40 mm; L1=12 mm; L2=30 mm

2 閃蒸相變數學模型和邊界條件

2.1閃蒸相變數學模型

閃蒸相變是一個復雜的多相流動傳質傳熱可壓縮流動過程[14-15]。由于難以確定閃蒸過程中氣泡直徑和成核速率等參數,本研究擬采用閃蒸模型耦合單場方法對閃蒸過程進行數值模擬[16-17]。混合模型假設介質為連續可相互穿插流體,并且允許相間存在小的滑移速度,其基本控制方程組包括質量守恒方程式(1),k相滑移速度定義式(2),k相質量方程式(3),動量守恒方程式(4),能量守恒方程式(5)。

(1)

vdr,k,j=vk,j-vm,j

(2)

(3)

(4)

(5)

式(1)中的混合物平均密度、平均速度分別由式(6)~(7)計算,式(5)中Ek的表達式見式(8)。

ρm=αlρl+αgρg

(6)

(7)

(8)

通過自定義函數(UDF)植入三相流連續方程(質量傳遞及守恒)、動量方程(三相動量傳遞及守恒)和能量方程(通過對流、導熱、輻射及相變引起的能量轉移過程),建立戊烷-瀝青體系閃蒸過程計算方法。在超臨界戊烷-瀝青體系的閃蒸相變過程中,瀝青相只發生動量和勢量交換,戊烷相除了動量和熱量交換外還發生相變,在液相戊烷和氣相戊烷間還存在質量交換,故還需計算戊烷相變造成的質量(Sm)和能量(SE)傳遞以封閉方程組,計算式為式(9)和式(10)。

(9)

(10)

2.2計算條件

模擬對象為氣、液戊烷和瀝青的三相閃蒸相變過程。模擬時,噴嘴入口施加壓力入口邊界條件,入口壓力和溫度分別為5 MPa和453.15 K,瀝青相體積分數為0.751;假設出口處流動已經局部單向化,施加壓力出口邊界條件,壁面處采用無滑移邊界條件。湍流模型采用RNGk-ε模型和標準壁面函數,考慮到輻射對溫度場的影響而添加了輻射模型。為保證計算穩定性和迭代收斂速率,壓力速度耦合算法采用針對非穩態可壓縮流動建立的PISO算法,對流項離散采用二階迎風格式,壓力離散格式考慮到流場具有旋轉和高曲率的性質,采用PRESTO算法。

3 結果與討論

3.1超臨界流體閃蒸相變模型驗證

圖3為超臨界流體閃蒸相變過程流場內速度和溫度分布云圖。由圖3可見,超臨界流體經噴嘴噴出,在出口處壓力驟降到常壓,從而導致超臨界流體在噴嘴出口處速度激增,并發生劇烈的閃蒸相變。壓力驟降使得內能部分轉化成動能,流體在噴嘴出口達最大速度,并且劇烈的閃蒸相變過程導致旋流噴嘴出口處速度場出現明顯的非均勻特性,如圖3(a)所示。壓力驟降導致戊烷在噴嘴出口處吸熱氣化,溫度沿射流方向逐漸降低,最終趨于穩定,形成一個狹長的高溫帶,如圖3(b)所示。模擬得到的出口處溫度為312.5 K。

圖3 超臨界流體閃蒸相變過程流場內 速度(v)和溫度(T)分布Fig.3 Velocity (v) and temperature (T) distributions during flash evaporation and phase change (a) v; (b) T

戊烷由5 MPa下的飽和溫度453.5 K降至常壓下的飽和溫度309.5 K,溶劑的氣化質量分數由式(11)計算。

(11)

由式(11)計算得到的戊烷氣化質量分數為87.69%,而數值模擬計算得到的氣化質量分數為84.07%。數值模擬結果與理論計算結果基本吻合,證明了閃蒸相變模型的可靠性。

3.2旋流噴嘴內閃蒸相變流場分析

圖4為旋流噴嘴內的切向速度、壓力和溫度分布。由圖4(a)看到,流體經旋流部件在噴孔段產生了旋轉流動,具有較大的切向速度。由圖4(b)看到,由于離心力的作用,噴嘴內壓力沿徑向逐漸增加,即中心區壓力低、邊壁區壓力高;沿軸向流體壓力逐漸降至常壓。由圖4(c)看到,旋流噴嘴內溫度分布比較均勻,沒有明顯變化。

圖4 旋流噴嘴內切向速度(vt)、壓力(p)、溫度(T)分布Fig.4 Tangential velocity(vt), pressure(p) and temperature(T) distribution in swirl nozzle (a) vt; (b) p; (c) T

圖5為旋流噴嘴內氣、液戊烷和瀝青體積分數分布云圖。由圖5可見,在離心力作用下,噴孔段三相介質形成分層流動:內層是部分汽化的氣相戊烷溶劑(見圖5(a));中間為尚未汽化的液相戊烷溶劑(見圖5(b));外側為密度最大的瀝青相(見圖5(c))。

圖5 旋流噴嘴內氣、液戊烷和瀝青體積分數分布Fig.5 Volume fraction distributions of vapor and liquid pentane and asphalt in swirl nozzle (a) Vapor pentane; (b) Liquid pentane; (c) Asphalt

圖6為旋流噴嘴噴孔段10 mm、20 mm、30 mm處的切向速度和壓力分布。由圖6(a)看到,旋流噴嘴內產生了較強的旋轉流動,切向速度沿徑向呈中心準剛性渦、外部準自由渦分布,切向速度沿徑向逐漸增大,在壁面附近流體受黏性力作用導致切向速度急劇減小;噴嘴內摩擦導致的沿程損失導致切向速度沿軸向逐漸減小。由圖6(b)看到,噴孔段不同位置處壓力分布基本一致,由于離心力的作用,壓力分布沿徑向逐漸增大,即中心區壓力低而邊壁壓力高,在噴嘴出口處內外壓差為1.4 MPa,中心壓力均低于1 MPa。

圖6 旋流噴嘴噴孔段不同位置處的切向速度(vt)和壓力(p)分布Fig.6 Tangential velocity(vt) and pressure(p) distributions in different positions of orifice section (a) vt; (b) p

圖7為旋流噴嘴噴孔段10 mm、20 mm、30 mm 處的氣、液戊烷和瀝青體積分數分布。由圖7(a)可見,噴嘴中心處壓力沿軸向劇烈降低,液態戊烷氣化率大大增加,在x=20 mm處液態戊烷已基本轉換完全,并在噴嘴中心處產生了“氣芯”現象,在離心力作用下,密度較小的介質分布在噴孔段中心區,氣相溶劑體積分數沿徑向逐漸減小,在邊壁附近接近0;噴孔段壓力沿軸向逐漸減小,導致氣相溶劑體積分數逐漸增大,在噴口處基本接近于1。由圖7(b)可見,尚未氣化的液相溶劑主要分布在氣相溶劑外側。由圖7(c)可見,在離心力作用下,密度最大的瀝青相分布在近壁區,體積分數基本與氣相溶劑體積分數分布相反。綜上所述,旋流噴嘴內流體介質分層流動,實現了氣相溶劑與瀝青相的預分離。

圖7 旋流噴嘴噴孔段不同位置處的氣、 液戊烷和瀝青體積分數分布Fig.7 Volume fraction distributions of vapor and liquid pentane and asphalt in different positions of orifice section (a) Vapor pentane; (b) Liquid pentane; (c) Asphalt

圖8為旋流噴嘴噴孔段10 mm、20 mm、30 mm 處的溫度分布。從圖8可以看出,由于溶劑在中心區部分氣化吸收了一定的熱量,導致中心區溫度降低了3~5 K,流體溫度沿軸向無明顯變化,噴嘴出口處瀝青相溫度仍大于其軟化點溫度(430~450 K),可以防止瀝青凝固而堵塞噴嘴。

圖8 旋流噴嘴噴孔段不同位置處的溫度分布Fig.8 Temperature distribution in different positions of orifice section

通過對旋流噴嘴流動相變的機理分析,建立了閃蒸相變模型,并對超臨界流體的閃蒸相變過程進行了研究。結果表明,旋流噴嘴內三相介質實現分層流動,弱化了氣相溶劑體積劇烈膨脹帶來的膨化作用,理論上能夠增大瀝青顆粒的表觀密度,對旋流噴嘴的優化改進和固體瀝青顆粒的成型具有較好的指導意義。

4 結 論

(1)閃蒸相變屬于復雜的非熱力學平衡過程,根據減壓相變過程的傳質傳熱原理建立了閃蒸相變模型,并與Fluent計算軟件中的混合模型耦合,計算旋流噴嘴內的相變過程。模擬得到的溫度和氣化率結果與理論計算吻合較好,驗證了利用模型研究閃蒸相變的可行性。

(2)戊烷溶劑在旋流噴嘴出口處大量氣化,導致形成多孔、蓬松的絮狀瀝青。旋流噴嘴內三相分層流動,實現了戊烷溶劑與瀝青相的預先分離,有利于形成粒徑較小且密實的瀝青顆粒。此外,還可以通過改變操作參數和結構參數來保證瀝青造粒質量。

(3)旋流噴嘴噴孔段長度較短,析出的瀝青相具有較大的軸向速度,而且噴嘴出口處溫度大于瀝青相軟化點溫度,這些都保證了所設計的旋流噴嘴不容易發生堵塞危險。

符號說明:

C——質量分數,%;

cp——比熱容,J/(kg·K);

D——直徑,mm;

E——比總能量,J/kg;

g——重力加速度,m2/s2;

h——比焓,J/kg;

keff——傳熱系數,W/(m2·K);

L——長度,mm;

p——壓力,MPa;

r——徑向位置,mm;

R——管徑,mm;

SE——能量傳遞源項,J;

Sm——質量傳遞源項,kg;

T——溫度,K;

T1,T2——溫度積分上、下限,K;

t——時間,s;

v——速度,m/s;

vdr——滑移速度,m/s;

vm——流體平均速度,m/s;

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

w——單位質量,kg;

x——模型方程通用坐標,mm;

α——控制方程各相體積分數;

γ——時間迭代松弛因子;

ρ——流體密度,kg/m3;

ρm——流體平均密度,kg/m3;

φ——體積分數;

下標

g, l——氣相,液相;

i,j——空間坐標;

k——氣相或液相。

參考文獻

[1] 徐春明,趙鎖奇,盧春喜,等. 重質油梯級分離新工藝的工程基礎研究[J].化工學報, 2010, 61(9): 2393-2400.(XU Chunming, ZHAO Suoqi, LU Chunxi, et al. Engineering basics of heavy oil deep stage separating process[J].Journal of Chemical Engineering of Chinese Universities, 2010, 61(9): 2393-2400.)

[2] 孫顯峰,孫學文,趙鎖奇,等.超臨界溶劑脫瀝青操作參數對遼河稠油減壓渣油脫油瀝青的影響[J].石油煉制與化工,2010, 41(2): 30-34.(SUN Xianfeng, SUN Xuewen, ZHAO Suoqi, et al. Influence of supercritical solvent deasphalting operation parameters on the de-oiled asphalt of Liaohe heavy crude vacuum residuum[J].Petroleum Processing and Petrochemicals, 2010, 41(2): 30-34.)

[3] 孫顯峰,孫學文,趙鎖奇,等.PGSS法用于脫油瀝青顆粒的制備[J].化工學報,2010, 24(2): 290-296.(SUN Xianfeng, SUN Xuewen, ZHAO Suoqi, et al. Preparation of de-oiled asphalt particles by PGSS process[J].Journal of Chemical Engineering of Chinese Universities, 2010, 24(2): 290-296.)

[4] 劉美麗,毛羽,王娟,等.氣粒兩相流與夾套耦合傳熱的數值模擬[J].化學反應工程與工藝,2011, 27(2): 121-132.(LIU Meili, MAO Yu, WANG Juan, et al. Numerical simulation of coupled heat transfer between gas-particle flow and jackets[J].Chemical Reaction Engineering and Technology, 2011, 27(2): 121-132.)

[5] 毛羽,徐春明,趙鎖奇,等.瀝青殘渣造粒多級相變噴霧造粒噴嘴:CN,101987287[P].2010-07-07.

[6] HENRY R E, FAUSKE H K. The two-phase critical flow of one-component mixtures in nozzles, orifices, and short tubes[J].Heat Transfer, 1971, 93(2): 179-187.

[7] WALLIS G B. Critical two-phase flow[J].International Journal of Multiphase Flow, 1980, 6(1-2): 97-112.

[8] NILPUENG K, WONGWISES S. Experimental investigation of two-phase flow characteristics of HFC-134a through short-tube orifices[J].International Journal of Refrigeration, 2009, 32(5): 854-864.

[9] SAURY D, HARMAND S, SIROUX B F. Flash evaporation from a water pool: Influence of the liquid height and of the depressurization rate[J].International Journal of Thermal Sciences, 2005, 44(10): 953-965.

[10] 郭迎利,鄧煒,嚴俊杰,等.初始條件對瞬態閃蒸過程的影響[J].工程熱物理學報,2008, 29(8): 1335-1338.(GUO Yingli, DENG Wei, YAN Junjie, et al. Influence of the initial condition on pool water instantaneous flash evaporation[J].Journal of Engineering Thermophysics, 2008, 29(8): 1335-1338.)

[11] 曹東剛,何國強,潘宏亮,等.三種空穴模型在可調汽蝕文氏管數值模擬中的對比研究[J].西北工業大學學報,2013, 31(4): 596-601.(CAO Donggang, HE Guoqiang, PAN Hongliang, et al. Comparative investigation among three cavitation models for simulating venture[J].Journal of Northwestern Polytechnical University, 2013, 31(4): 596-601.)

[12] KUNZ R F, BOGER D A, STINEBRING D R, et al. A preconditioned Navier-Stokes method for two-phase flow with application to cavitation prediction[J].Computers Fluids, 2000, 29(8): 849-875.

[13] 聶永廣,毛羽,王江云,等.高壓射流中的戊烷閃蒸過程數值模擬[J].石油學報(石油加工),2012, 28(5): 814-821.(NIE Yongguang, MAO Yu, WANG Jiangyun, et al. Numerical simulation of flash evaporation in high-pressure pentane jet[J].Acta Petrolei Sinica (Petroleum Processing Section), 2012, 28(5): 814-821.)

[14] NEROORKAR K, SCHMIDT D. Model of vapor-liquid equilibrium of gasoline-ethanol blended fuels for flash boiling simulations[J].Fuel, 2011, 90(2): 655-673.

[15] SHER E, BAR-KOHANY T, RASHKOVAN A. Flash-boiling atomization[J].Progress in Energy and Combustion Science, 2008, 34(4): 417-439.

[16] PRASANTH K S, Prasad B V S S, VENKATARATHNAM G, et al. Influence of surface evaporation on stratification in liquid hydrogen tanks of different aspect ratios[J].International Journal of Hydrogen Energy, 2007, 32(12): 1954-1960.

[17] 聶永廣,毛羽,王江云,等.超臨界瀝青噴霧造粒中的戊烷閃蒸CFD模擬[J].化學工業與工程,2013, 30(4): 68-72.(NIE Yongguang, MAO Yu, WANG Jiangyun, et al. CFD simulation of supercritical pentane flash boiling in asphalt granulation[J].Chemical Industry and Engineering, 2013, 30(4): 68-72.)

收稿日期:2015-04-22

基金項目:國家重點基礎研究發展規劃“973”項目(2010CB226902)基金資助

文章編號:1001-8719(2016)04-0741-07

中圖分類號:TE65

文獻標識碼:A

doi:10.3969/j.issn.1001-8719.2016.04.012

Numerical Simulation of Flash Evaporation for Supercritical Fluid in Swirl Nozzle

FENG Liuhai1,2, WANG Jiangyun1, ZHAO Fan1,3, SUN Zhongwei1,4, WANG Juan1, MAO Yu1

(1.StateKeyLaboratoryofHeavyOilProcessing,ChinaUniversityofPetroleum,Beijing102249,China;2.NationalInstituteofClean-and-LowCarbonEnergy,Beijing102209,China;3.LanzhouLSEnergyEquipmentEngineeringInstituteCo.Ltd.,Lanzhou730314,China;4.ENNScience&DevelopmentCo.Ltd.,Langfang065001,China)

Abstract:The volume fraction distribution of the asphalt solute in supercritical solution at swirl nozzle obviously affects the particle formation. Based on the theory of heat and mass transfer in phase change caused by abrupt pressure drop, the flash evaporation model suitable for supercritical fluid model was developed and added to Fluent by the UDF method. By using the flash evaporation model coupling with multiphase mixture model the process of the flash evaporation in the swirl nozzles was studied. A detailed analysis of the distributions of pressure, velocity, temperature and concentration of each phase inside nozzle was carried out in order to predict the important effect of swirl nozzle on pellet-forming. Simulation results showed that there was three-phases stratified flow in the swirl nozzle, so as the pre-separation of pentane solvent and asphalt solute was realized, which was conducive to form smaller and more compact particles.

Key words:swirl nozzle; numerical simulation; flash evaporation; thermal non-equilibrium

第一作者: 馮留海,男,博士研究生,從事多相流動的數值模擬與實驗方面的研究

通訊聯系人: 王江云,男,助理研究員,博士,從事多相流動與分離、腐蝕及燃燒過程的數值模擬與實驗方面的研究;Tel:010-89733293;E-mail:wangjy@cup.edu.cn;毛羽,男,教授,博士,從事多相流動及燃燒、氣固分離及液體霧化技術、化工過程裝備優化等方面的研究;Tel:010-89733293;E-mail:maoyu@cup.edu.cn

主站蜘蛛池模板: 亚洲日韩高清无码| 精品91在线| 国产成a人片在线播放| 人妻精品久久无码区| 91视频区| 91在线一9|永久视频在线| 亚洲 成人国产| 国产波多野结衣中文在线播放| 黄色网站不卡无码| 91口爆吞精国产对白第三集 | 亚洲欧美另类色图| 亚洲天堂福利视频| 国产成本人片免费a∨短片| 污网站在线观看视频| 一级一级一片免费| 成人蜜桃网| 国产美女无遮挡免费视频| 欧美精品三级在线| 永久免费无码成人网站| 欧美www在线观看| 婷婷亚洲视频| 2048国产精品原创综合在线| 福利在线一区| 91久久偷偷做嫩草影院精品| 欧洲精品视频在线观看| 热99精品视频| 97国产精品视频人人做人人爱| aaa国产一级毛片| 国产精品短篇二区| 成人福利在线视频| 国产 在线视频无码| 国产精品第页| 狼友视频一区二区三区| 国产尹人香蕉综合在线电影| 日日拍夜夜嗷嗷叫国产| 国产一区二区精品福利| 国产精品嫩草影院视频| 2022国产无码在线| 无码内射在线| 亚洲精品另类| 天天综合亚洲| 久久久久免费看成人影片| 香蕉国产精品视频| 国产香蕉一区二区在线网站| 亚洲午夜综合网| V一区无码内射国产| 国产高清在线观看| 国产91高清视频| 1024国产在线| 国产91在线|中文| 国产激爽大片高清在线观看| 亚洲精品男人天堂| 亚洲无码高清一区二区| 国产精鲁鲁网在线视频| 久久成人国产精品免费软件 | 国产精品免费p区| 成年A级毛片| 四虎亚洲国产成人久久精品| 香蕉在线视频网站| 国产导航在线| 538国产在线| 国产区免费精品视频| 欧美精品一区在线看| 凹凸国产分类在线观看| 2021天堂在线亚洲精品专区| 日韩福利视频导航| 性色在线视频精品| 91精品最新国内在线播放| 人人91人人澡人人妻人人爽| 久久国产精品77777| 成人永久免费A∨一级在线播放| 亚洲人成网站日本片| 亚洲色图欧美激情| 91在线高清视频| 国产综合精品日本亚洲777| 999精品视频在线| 99久久精品久久久久久婷婷| 色妞永久免费视频| 国产亚洲男人的天堂在线观看| 美女高潮全身流白浆福利区| 最新亚洲人成无码网站欣赏网| 久久这里只有精品23|