亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        減壓轉(zhuǎn)油線內(nèi)氣液兩相流動的計算流體力學(xué)模擬

        2017-04-11 02:50:24張呂鴻趙宗賢肖曉明
        化學(xué)工業(yè)與工程 2017年4期
        關(guān)鍵詞:渣油液滴氣相

        張呂鴻,趙宗賢,肖曉明

        (天津大學(xué)化工學(xué)院,天津 300072)

        減壓轉(zhuǎn)油線是原油減壓蒸餾單元中連接減壓加熱爐和減壓蒸餾塔的重要管線,其主要職責(zé)是將高溫的常壓渣油在較低壓降下輸送至減壓塔[1-2]。

        減壓轉(zhuǎn)油線內(nèi)輸送的渣油為氣液兩相碳氫混合物。在傳輸過程中,渣油在流動方向上受到阻力作用,壓力持續(xù)下降。在這個過程中,渣油液滴中可揮發(fā)組分逐漸汽化,使渣油的汽化率逐漸提高。由于汽化率高低直接影響渣油的拔出率,因此提高渣油氣化率是減壓轉(zhuǎn)油線設(shè)計好壞的重要考察因素。此外,為了保證減壓塔各餾出分的質(zhì)量,應(yīng)當(dāng)使減壓爐維持在較低溫度下操作。這就要求降低轉(zhuǎn)油線的溫度降控制在可接受范圍內(nèi)。

        為了要設(shè)計出優(yōu)秀的減壓轉(zhuǎn)油線,對轉(zhuǎn)油線內(nèi)部流場分布情況的了解是極重要的一環(huán)。目前對減壓轉(zhuǎn)油線內(nèi)流場的研究較少,且研究方法主要分成2類,一是一維閃蒸模型[3-6],二是采用三維計算流體力學(xué)分析[3,7]。對于一維閃蒸模型,首先是將管道抽象為一系列相互連接的單元。在每個單元內(nèi),假設(shè)渣油氣液相均勻混合,且溫度相等[2-4]。每個單元內(nèi)氣液相的交換則通過閃蒸計算獲得。流場的壓降通過兩相流壓降經(jīng)驗關(guān)聯(lián)式進行計算[8-11]。一維模型計算效率較高,但無法獲得管內(nèi)流場的具體情況。對于三維計算流體力學(xué)分析,目前研究主要采用VOF(volume of fluid)模型[3]和DPM(discrete phase model)模型[7]。由于液相渣油在占有較大質(zhì)量分數(shù)的同時其體積分數(shù)極小(<1%),采用VOF模型就需要極精密的網(wǎng)格。為了處理數(shù)量巨大的網(wǎng)格,需要采用分段計算的方式[3],這就導(dǎo)致了管道連接處的失真。DPM模型忽略渣油體積,僅考慮液滴與連續(xù)相的動量質(zhì)量傳遞[7]。目前的研究簡化了渣油蒸發(fā)過程,并將渣油蒸發(fā)速率定一個恒定的傳質(zhì)速率。而復(fù)雜碳氫混合物的蒸發(fā)過程是無法用一個常數(shù)準確模擬的[13-14]。

        本研究綜合前人研究,采用了DPM模型對減壓轉(zhuǎn)油線內(nèi)流場進行模擬。首先采用虛擬組分法對氣液相渣油物理性質(zhì)進行計算。而后為了精確模擬液滴表面?zhèn)鳠崤c傳質(zhì)的過程,考慮了液滴表面溫度對傳質(zhì)速率的影響。最后,分析了管道內(nèi)氣液相壓力、溫度、產(chǎn)物組成等因素,以期為減壓轉(zhuǎn)油線的設(shè)計提供指導(dǎo)。

        1 研究方法

        首先介紹了氣液相渣油物理性質(zhì)的描述方法,然后給出多相流模型的控制方程和計算方法。

        1.1 氣液相渣油物理性質(zhì)描述

        采用虛擬組分法描述渣油組成。虛擬組分法將渣油分為一系列窄組分,每個組分當(dāng)做渣油的一種組成。本研究依照組分的沸點對渣油進行虛擬組分分割:將高于入口溫度和低于出口溫度的組分各劃分為1種組分,而沸點在管道溫度范圍內(nèi)的組分劃分為4個組分。虛擬組分的平均沸點、平均相對分子質(zhì)量、API重度和臨界性質(zhì)由PRO/II軟件進行計算。理想氣體焓以及蒸發(fā)熱采用Cavett方法[15]計算。渣油液滴的密度采用API[16]經(jīng)驗關(guān)聯(lián)式。飽和蒸汽壓通過估算不同溫度和壓力下的氣液平衡常數(shù)獲得。由于管內(nèi)的壓力較低,假設(shè)氣相渣油為理想混合物。渣油液滴的表面張力使用Twu[17]經(jīng)驗關(guān)聯(lián)式計算。最后,氣相與液相渣油的黏度采用Letsou-Stiel[18]和Dean-Stiel方程[19]進行計算。

        1.2 連續(xù)相方程

        流體流動與傳熱過程的計算基于3個基本的物理定律,即質(zhì)量守恒、動量守恒和能量守恒。以下將具體介紹連續(xù)相中運用到的理論方程。

        1.2.1控制方程

        將質(zhì)量守恒定律運用于計算流體力學(xué)中,可以表述為單位時間內(nèi)流入流體微元中的凈質(zhì)量,等于該時間內(nèi)流體微元質(zhì)量的凈增量,具體如式(1)所示[20]。

        (1)

        (2)

        其中,P代表靜壓;而ρgi為i方向上重力的體積力;Fi為其他模型引入的體積力源項;τij為應(yīng)力張量。多組分的能量守恒方程如式(3)所示:

        (3)

        1.2.2湍流模型

        湍流模型中最常見的標準k-ε模型,可以運用于絕大多數(shù)的工業(yè)湍流計算當(dāng)中,但是標準k-ε模型各向同性的假設(shè),將在高速流動的近壁區(qū)域產(chǎn)生較大誤差。因此本研究采用Realizablek-ε模型[21]對減壓轉(zhuǎn)油線內(nèi)氣相流動進行模擬。

        1.3 離散相方程

        1.3.1渣油液滴的控制方程

        渣油液滴的受力方程如式(4)所示[21]。

        (4)

        其中FD是指單位質(zhì)量渣油液滴所受到的拖曳力。

        渣油液滴的模擬還包括液滴與氣相的傳質(zhì)過程。首先假設(shè)液滴內(nèi)物質(zhì)均勻混合,這是因為液滴內(nèi)傳質(zhì)較氣相快得多[22]。由于液滴內(nèi)傳熱速率較氣相要慢,因此認為液滴表面與內(nèi)部存在溫差。本研究采用了考慮表面溫度的傳質(zhì)模型,其中液滴傳質(zhì)過程采用式(5)[21]。

        (5)

        1.3.2渣油蒸發(fā)過程迭代求解

        渣油液滴計算的初始化采用收斂的連續(xù)相流場解。

        本研究耦合求解物質(zhì)與能量的傳輸過程,需要對每1個液滴表面進行傳質(zhì)過程迭代。由于減壓轉(zhuǎn)油線內(nèi)壓力較低,故假設(shè)氣相組分為理想混合物。渣油液滴表面的組分濃度與氣相組分關(guān)聯(lián)采用拉烏爾定律(Raoult's law),并寫成式(6)形式。

        pi,v=xi,vP=xi,lPsat,i

        (6)

        pi,v是氣相中組分i的分壓,而xi,v和xi,l是組分i在氣相和液滴中的組分摩爾分數(shù)。P是氣相的總壓,Psat,i是組分i在當(dāng)前溫度下的飽和蒸汽壓。

        液滴表面的溫度通過計算液滴表面與其附近物質(zhì)與能量平衡獲得。其中,熱傳導(dǎo)包含2個部分,一是液滴周邊的氣體混合物傳導(dǎo)至液滴表面的熱量(qo),二是液滴內(nèi)部與表面間的熱傳導(dǎo)(qi)。液滴表面附近的熱平衡方程如式(7)。

        (7)

        L指的是液滴表面渣油在表面溫度Ts下的汽化潛熱。液滴內(nèi)部熱傳導(dǎo)的計算考慮了液滴內(nèi)的對流傳熱過程。液滴內(nèi)部導(dǎo)熱系數(shù)hi,eff可以由熱傳導(dǎo)系數(shù)hi,eff和非穩(wěn)態(tài)等量熱邊界層厚度δe來計算,且表達式如式(8)。

        (8)

        本研究中,外部流場的有效熱傳導(dǎo)系數(shù)采用Reitz等[13,20]推導(dǎo)出的導(dǎo)熱方程,并考慮了內(nèi)部擴散和Stefan流動對液滴導(dǎo)熱的影響。該方程具體如式(9)。

        qo=ho,eff(Tsur-Ts)=

        (9)

        (10)

        (11)

        BM=(yFs-yFsur)/(1-yFs)

        (12)

        并且,需要設(shè)定2個溫度。其一為蒸發(fā)溫度Tvap。當(dāng)溫度低于Tvap時,由于蒸發(fā)速度緩慢,假設(shè)蒸發(fā)未發(fā)生;當(dāng)溫度高于Tvap時,則開始運用上述的渣油蒸發(fā)計算。

        最后,熱傳導(dǎo)系數(shù)的計算需要與液滴表面的傳質(zhì)過程相耦合。由于方程的高度非線性,具體求解需要一個迭代的過程。迭代首先使用時間步初期渣油液滴和液滴所在網(wǎng)格內(nèi)氣相的數(shù)據(jù)進行計算,并獲得方程(9)左右兩側(cè)的差值。之后將該差值乘上松弛因子作為迭代殘差,并以源項方式導(dǎo)入FLUENT中求解。

        1.3.3邊界條件設(shè)定

        關(guān)于離散相模型邊界條件的選擇,轉(zhuǎn)油線的入口和出口處設(shè)定為逃逸(escape)邊界條件。到達此處的液滴不在返回計算域,并脫離流場計算。其他管道內(nèi)壁邊界條件采用壁面膜(Wall-Film)模型。當(dāng)液滴與壁面接觸時,將根據(jù)接觸角與相對速度,自動判斷液滴是濺起或是附著在壁面上。液滴入口采用surface條件,渣油液滴通過surface進入計算域。

        2 實例分析

        選取3個工業(yè)上使用的減壓轉(zhuǎn)油線(外形如圖1)進行計算流體力學(xué)分析。管道內(nèi)入口和出口處的溫度和壓力為工廠中實地測量[1,7]。其中,壓力的測量采用U型管壓差計和精密真空表。溫度的測量選用熱電偶。管道的結(jié)構(gòu)參數(shù)主要包括入口與出口的管徑、管道長度、彎管式樣以及特殊部件(如褲狀三通等)。對于所分析的減壓轉(zhuǎn)油線,其結(jié)構(gòu)參數(shù)如表1所示。

        2.1 網(wǎng)格劃分

        運用ICEM 14.5軟件對幾何模型進行網(wǎng)格劃分。由于網(wǎng)格的規(guī)律性較強,故全部采用結(jié)構(gòu)化的六面體網(wǎng)格,使之在較少網(wǎng)格數(shù)下獲得較高質(zhì)量的數(shù)值解。如圖2所示,網(wǎng)格繪有良好邊界層。各個減壓轉(zhuǎn)油線幾何體的網(wǎng)格數(shù)量均控制在200萬以內(nèi)(A轉(zhuǎn)油線147萬,B轉(zhuǎn)油線172萬,C轉(zhuǎn)油線185萬)。其中,較小的網(wǎng)格主要存在于管道壁面附近以捕捉邊界層。而較大網(wǎng)格主要位于轉(zhuǎn)油線大直徑管道末端,因為這部分流場的流速較小。

        圖1 3個工業(yè)中減壓轉(zhuǎn)油線例子Fig.1 Case study of three transfer lines

        表1 減壓轉(zhuǎn)油線結(jié)構(gòu)參數(shù)

        圖2 結(jié)構(gòu)化的六面體網(wǎng)格Fig.2 Structural hex mesh

        關(guān)于網(wǎng)格的無關(guān)性驗證,由于減壓轉(zhuǎn)油線結(jié)構(gòu)較簡單,以較少網(wǎng)格即可獲得高精確的解。而渣油液滴的計算需要較精密網(wǎng)格。為了節(jié)約計算時間,采用自適應(yīng)網(wǎng)格技術(shù),即在流場中選擇速度梯度較大區(qū)域細化網(wǎng)格。

        2.2 計算結(jié)果與討論

        2.2.1轉(zhuǎn)油線內(nèi)壓力分布

        壓力分布模擬結(jié)果的相對誤差展現(xiàn)在表2中,可以看出多相流模型的壓力模擬較為精確,且相對誤差均在5%左右。圖3為減壓轉(zhuǎn)油線內(nèi)壓力沿管徑方向的分布情況。可以看出小管徑管道內(nèi)壓力變化明顯,而在接近出口的大直徑管道內(nèi)壓力變化較小。

        圖3 轉(zhuǎn)油線內(nèi)壓力分布曲線Fig.3 Pressure distributions in transfer lines

        此外還可以看出,各管道內(nèi)壓力分布的主要區(qū)別在于異徑管道的接頭處。A管道在合流處壓力陡降,B管道壓力出現(xiàn)2次下降,C管道內(nèi)壓力的變化則較為平緩。其主要原因是這3種管道分別采取了不同的合流方式。對于A型管,它的入口管道直插入大直徑管道。而C型管道2次在合流處采用褲狀三通,壓力波動明顯較小。B管道壓降變化陡峭程度介于在A和C管之間??梢钥闯鲈跍p壓轉(zhuǎn)油線內(nèi)采用褲狀三通和逐級擴徑的方式,可以使管內(nèi)流場較為平緩。

        2.2.2減壓轉(zhuǎn)油線內(nèi)溫度分布

        溫度模擬結(jié)果與真實值相對誤差見表2。管內(nèi)溫度分布圖如圖4。

        表2 模擬結(jié)果與真實數(shù)據(jù)的相對誤差

        圖4 轉(zhuǎn)油線內(nèi)溫度分布曲線Fig.4 Temperature distributions in transfer lines

        由圖4可以看出溫降主要集中在轉(zhuǎn)油線較小直徑的管道中。此外,溫度變化曲線的趨勢和壓力變化是很相似的,并在管道擴徑處流場溫度發(fā)生突降。這主要是管內(nèi)氣體進入大直徑管道時發(fā)生突然膨脹而導(dǎo)致的。同樣,A轉(zhuǎn)油線由于管道突然擴張,溫度變化相對較為劇烈。而B和C管道內(nèi)由于2次擴徑,溫度逐漸降低,變化過程相對平緩。

        2.2.3減壓轉(zhuǎn)油線內(nèi)速度分布

        轉(zhuǎn)油線內(nèi)流場的速度分布如圖5所示。

        由圖5可見,管道內(nèi)流速在管道合流處劇烈升高,之后再急劇下降。流速劇烈升高的主要原因是管道在該處擴徑,管內(nèi)流體具體膨脹。之后大直徑管道吸收了氣體膨脹的體積,使流速再次回落。此外還可以看出,采用逐級擴徑以及褲裝三通的C轉(zhuǎn)油線速度峰值最小。而采用直插式擴徑的A轉(zhuǎn)油線速度極值最大。

        圖5 轉(zhuǎn)油線內(nèi)流速分布曲線Fig.5 Velocity distributions in transfer lines

        2.2.4渣油液滴的汽化

        本研究僅監(jiān)控入口與出口處渣油液滴中各組分的含量,并將結(jié)果顯示于表3中。渣油各個組分以渣油的平均沸點進行命名,溫度單位為攝氏度。其中A轉(zhuǎn)油線的入口溫度為394.7 ℃,B轉(zhuǎn)油線為392.3 ℃,C轉(zhuǎn)油線為389.0 ℃。對于出口溫度,A轉(zhuǎn)油線為380.2 ℃,B轉(zhuǎn)油線為380.9 ℃,C轉(zhuǎn)油線為374.0 ℃。

        表3 各減壓轉(zhuǎn)油線入口與出口處渣油液滴各組分質(zhì)量分數(shù)

        由表3可以看出,轉(zhuǎn)油線內(nèi)對于平均沸點高于入口的渣油組分,蒸發(fā)的效果不明顯。主要蒸發(fā)的渣油為輕組分,且平均沸點低于轉(zhuǎn)油線出口溫度的渣油組分基本蒸發(fā)完畢。除了最重的組分,其余組分也依次部分蒸發(fā),且在液滴內(nèi)所占質(zhì)量分數(shù)有所下降。

        2.2.5液滴碰撞壁面情況

        根據(jù)減壓轉(zhuǎn)油線內(nèi)渣油液滴運動的軌跡圖,可以看出幾個液滴碰撞壁面的區(qū)域。分別為小直徑管道的彎管部位和合流處大直徑管道的壁面,具體位置如圖6所示。由于減壓轉(zhuǎn)油線的振動問題也是管道設(shè)計重要的考慮因素,這些渣油液滴作為激振源,將導(dǎo)致該部分管路的劇烈振動。因此施工過程中需給予額外固定。

        圖6 渣油液滴與管道碰撞的情況Fig.6 Collisions between residue droplets and walls

        3 結(jié)論與展望

        通過建立計算流體力學(xué)模型,對減壓轉(zhuǎn)油線內(nèi)伴隨氣液相變的兩相流進行模擬與分析,旨在對轉(zhuǎn)油線結(jié)構(gòu)優(yōu)化做出一定指導(dǎo)。獲取了3套工業(yè)中減壓轉(zhuǎn)油線的結(jié)構(gòu)參數(shù)、渣油物性和操作條件。從壓力、溫度、流速、渣油液化以及液滴軌跡這些角度對轉(zhuǎn)油線內(nèi)流場分布進行了理論研究,并得出以下結(jié)論:

        1)多相流模型雖然計算量大,但可以模擬出準確趨勢,且在模擬結(jié)果上逼近工業(yè)值;2)減壓轉(zhuǎn)油線入口管道合流方式的不同,直接影響管道內(nèi)的壓力、溫度等的分布;3)管道內(nèi)汽化的渣油主要是平均沸點低于管道入口溫度的較輕渣油組分;4)入口處小直徑管道應(yīng)兩次采用褲狀三通匯入大直徑管道,以使流場變化平穩(wěn),并減小管道的振動;5)管道內(nèi)液滴隨流場運動,且大多碰撞了彎管中曲率半徑較大的外側(cè)處管壁,需要對這部分管件加固,以預(yù)防管道的振動。

        參考文獻:

        [1]李哲. 減壓轉(zhuǎn)油線的設(shè)計[J]. 當(dāng)代化工, 2005, 34(6): 389-391

        Li Zhe. The design of vacuum transfer lines[J]. Contemporary Chemical Industry, 2005, 34(6): 389-391(in Chinese)

        [2]陳建民, 孫玉玲, 姜斌, 等. 典型工業(yè)減壓轉(zhuǎn)油線結(jié)構(gòu)和運行參數(shù)分析[J]. 化工進展, 2011, 30(7): 1 450-1 455

        Chen Jianmin, Sun Yuling, Jiang Bin,etal. Analysis of structure and operating parameters of typical industrial vacuum transfer line[J]. Chemical Industry and Engineering Progress, 2011, 30(7): 1 450-1 455(in Chinese)

        [3]秦婭. 減壓轉(zhuǎn)油線氣液兩相流動特性模擬及結(jié)構(gòu)優(yōu)化研究[D]. 天津: 天津大學(xué), 2009

        Qin Ya. A study on modeling of vapor liquid two-phase flow in the transfer line and structure optimization[D]. Tianjin: Tianjin University, 2009(in Chinese)

        [4]秦婭, 李鑫鋼, 姜斌, 等. 閃蒸二相流模型在減壓轉(zhuǎn)油線改造中的應(yīng)用[J]. 化學(xué)工程, 2010, 38(3): 26-29

        Qin Ya, Li Xingang, Jiang Bin,etal. Application of two-phase flashing flow model to retrofit of a vacuum transfer line[J]. Chemical Engineering(China), 2010, 38(3): 26-29(in Chinese)

        [5]Qin Y, Han X, Wang H,etal. Modeling of two-phase flashing flow of multicomponent mixtures in large diameter pipes[J]. Chemical Engineering and Technology, 2008, 31(11): 1 676-1 684

        [6]肖家治, 李發(fā)永, 王萬里, 等. 轉(zhuǎn)油線工藝設(shè)計方法的研究[J]. 石油大學(xué)學(xué)報: 自然科學(xué)版, 1996, 20(2): 90-92

        Xiao Jiazhi, Li Fayong, Wang Wanli,etal. Investigation on the design method of transfer lines[J]. Journal of the University of Petroleum, China, 1996, 20(2): 90-92(in Chinese)

        [7]張慶軍, 減壓轉(zhuǎn)油線的數(shù)值模擬研究[D]. 天津: 天津大學(xué), 2008

        [8]Kohoutek J, Zachoval J, Odstrcil M,etal. Solving practical industrial problems in two-phase multicomponent mixture flow-critical velocity [J]. Heat Transfer Engineering, 2001, 22: 32-24

        [9]Beggs D, Brill P. A study of two-phase flow in inclined pipes[J]. Journal of Petroleum Technology, 1973, 25: 607-617

        [10]Dukler E, Wicks M, Cleveland G. Frictional pressure drop in two-phase flow: B an approach through similarity analysis[J]. AIChE Journal, 1964, 10(1): 44-51

        [11]Mukherjee H, Brill P. Empirical equations to predict flow patterns in two-phase inclined flow[J]. International Journal of Multiphase Flow, 1985, 11(3): 299-315

        [12]Hagedorn R, Brown E. Experimental study of pressure gradients occurring during continuous two-phase flow in small-diameter vertical conduits[J]. Journal of Petroleum and Technology, 1965, 17: 475-485

        [13]Youngchul R, Rolf R. A vaporization model for discrete multi-component fuel sprays[J]. International Journal of Multiphase Flow, 2009, 35: 101-117

        [14]Sergei S. Advanced models of fuel droplet heating and evaporation[J]. Progress in Energy and Combustion Science, 2006, 32: 162-214

        [15]Simulation Science Inc. Pro/II reference manual[M]. USA: CA, 2008

        [16]Daubert E, Danner P. API technical data book——Petroleum refining[M]. 6th edition. Washington: American Petroleum Institute (API), 1997

        [17]Twu H. An internally consistent correlation for predicting the critical properties and molecular weights of petroleum and coal-tar liquids[J]. Fluid Phase Equilibrium, 1984, 16: 137-151

        [18]Letsou A,Stiel I. Viscosity of saturated nonpolar liquids at elevated pressures[J]. AIChE Journal, 1973, 19: 409-412

        [19]Dean G, Stiel S. The viscosity of nonpolar gas mixtures at moderate and high pressures[J]. AIChE Journal, 1965, 11: 526-533

        [20]Versteeg K, Malalasekera W. An introduction to computational fluid dynamics: The finite volume method[M]. Edinburgh: Pearson Education Limited, 2007

        [21]ANSYS Inc. Fluent User′s Guide 14.5[M]. Canonsburg,PA: ANSYS Inc., 2012

        [22]Ra Y, Reitz D. The application of a multi-component vaporization model to gasoline direct injection engines[J]. International Journal of Engine Research, 2003, 4: 193-218

        猜你喜歡
        渣油液滴氣相
        用油漿渣油調(diào)合10號建筑瀝青的研究
        石油瀝青(2023年5期)2023-12-08 08:35:04
        基于分子結(jié)構(gòu)的渣油沸騰床加氫轉(zhuǎn)化特點研究
        氣相過渡金屬鈦-碳鏈團簇的研究
        液滴間相互碰撞融合與破碎的實驗研究
        噴淋液滴在空氣環(huán)境下的運動特性
        塔河渣油重溶劑脫瀝青深度分離研究
        石油瀝青(2019年4期)2019-09-02 01:41:56
        中國石化石油化工科學(xué)研究院開發(fā)RHT-200系列渣油加氫催化劑
        新型釩基催化劑催化降解氣相二噁英
        預(yù)縮聚反應(yīng)器氣相管“鼓泡”的成因探討
        氣相防銹技術(shù)在電器設(shè)備防腐中的應(yīng)用
        少妇勾引视频网站在线观看| 黑人巨大av在线播放无码| 粗一硬一长一进一爽一a级| 禁止免费无码网站| 一区二区三区国产天堂| 日本免费视频| 亚洲精品在线视频一区二区| 精品国产乱码久久久久久婷婷| 午夜一区欧美二区高清三区| 亚洲av无码一区二区乱子仑| 北岛玲亚洲一区二区三区| 亚洲一区毛片在线观看| 人妻少妇精品无码专区二区 | 日本不卡在线视频二区三区| 久久久久亚洲av无码专区首| 国产精品三级在线观看无码| 妺妺窝人体色www聚色窝韩国| 国产成人av一区二区三| 国产狂喷水潮免费网站www| 国产成人无码一区二区三区在线| 青青国产成人久久91| 精品人妻一区二区三区不卡毛片| 亚洲国产精品综合久久网络 | 女同性恋亚洲一区二区| 国产流白浆视频在线观看| 色噜噜狠狠狠综合曰曰曰| 天天做天天躁天天躁| 美女裸体无遮挡免费视频国产| 国产一品二品三区在线观看| 一区二区三区中文字幕| 国产精品黄网站免费观看| 久久99精品免费国产| 一二三四区中文字幕在线| 少妇高潮尖叫黑人激情在线| 亚洲熟女av超清一区二区三区| 亚洲av乱码国产精品观| 日韩人妻无码一区二区三区久久| 国产黄页网站在线观看免费视频| 国产毛片三区二区一区| 麻豆91蜜桃传媒在线观看| 99精品免费久久久久久久久日本 |