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

萃取塔數學模型及過程強化的若干研究進展

2020-06-29 04:10:04譚博仁李龍祥王勇齊濤
化工進展 2020年6期
關鍵詞:模型

譚博仁,李龍祥,王勇,齊濤

(1 中國科學院過程工程研究所,濕法冶金清潔生產技術國家工程實驗室,北京100190;2 中國科學院大學,北京100190;3 河北工程大學,河北邯鄲056000)

溶劑萃取是一種重要的化工分離技術,萃取設備主要包括:混合澄清槽、萃取塔和離心萃取器等。萃取塔具有結構簡單、密封好、生產能力大等優點,在石油、化工、生物、醫藥和環境工程等諸多領域被廣泛應用[1-2]。萃取塔根據內構件和外加能量形式不同可以分為填料萃取塔、篩板萃取塔、脈沖板環萃取塔、轉盤萃取塔、振動篩板塔、Kühni萃取塔等。

近年來,針對已有各種萃取塔在水力學、軸向擴散與相間傳質等方面開展了系統、深入的實驗研究和數學模型研究,大量可靠的經驗數學模型得以建立。隨著計算模擬技術的發展,計算流體力學(CFD)因全面、直觀等優點在萃取塔模擬中逐漸得到應用,相關模擬也從單一體系向多相復雜體系發展。在萃取過程強化方面也有一系列新技術產生,包括內構件改進、外加能量方式耦合、加入納米磁性顆粒等。本文作者將從萃取塔經驗模型、流體力學模擬以及萃取過程強化三方面對近年來萃取塔研究進展進行闡述。

1 萃取塔模型

1.1 水力學模型

萃取塔的水力學性能包括分散相持液量、特征速度、液滴直徑等,這對萃取塔的指導生產操作有重要的意義。

(1)Karr往復篩板萃取塔 Stella[3]、Smith[4]等對Karr往復板篩板萃取塔進行了系統研究,發現Karr萃取塔的水力學性能與塔徑無關,但受塔板的潤濕性影響較大,可分別采用Kumar 與Hartland[5-6]的模型對Karr 萃取塔的分散相持液量與液滴直徑進行良好的預測。Christo 等[7]研究了不同材料的塔板對持液量的影響,結果表明水相為連續相時聚四氟乙烯板的特征速度最低,而有機相為連續相時,不銹鋼板的特征速度最低。因此,將潤濕特性參數引入到Kumar 與Hartland[6]的模型中,得到了較好的預測效果。

(2)脈沖篩板萃取塔 Yadav等[8]采用50年來不同脈沖篩板萃取塔的水力學實驗數據,分析了不同數學模型的準確性。結果表明:Tribess 與Brunello[9]的液泛通量預測模型、Sreenivasulu等[10]的液滴直徑預測模型、Venkatnarsaiah與Verma[11]的持液量預測模型適用性最強。Sarkar 等[12]則在Venkatnarsaiah與Verma[11]模型的基礎上對參數進行了調整,提高了預測準確性。Sen 等[13]對脈沖篩板萃取塔與脈沖板環萃取塔的水力學性能進行了比較,研究發現脈沖篩板萃取塔持液量受脈沖強度影響更大、持液量更低。

(3)脈沖板環萃取塔 Grabin 等[14]以水相為分散相,采用聚四氟乙烯、尼龍和不銹鋼材質研究了內構件潤濕性對脈沖板環萃取塔水力學性能的影響(圖1)。在高脈沖強度下,特征速度隨板材料疏水性的增加而增大;低脈沖強度下,分散相持液量受內構件潤濕性能影響較大,不銹鋼材質內構件的分散相持液量最大。

圖1 不同內構件潤濕性對液滴流動行為的影響[15]

Torab-Mostaedi 等[16]研究發現,表面張力對脈沖板環萃取塔的分散相持液量有顯著影響,在不同操作區間建立了相應的分散相持液量預測模型,Sarkar等[12]在此模型的基礎上,根據直徑為50.8mm與76.0mm 的脈沖板環萃取塔分散相持液量數據,對模型參數進行了修正。Liu 等[17]研究了水相或有機相為連續相時,帶中心柱的脈沖板環萃取塔水力學性能,并采用Kumar 與Hartland[6,18]分散相持液量、液滴直徑預測模型在不同操作區間建立了對應的水力學模型。Wang等[19-20]研究了有脈沖與無脈沖條件下的板環萃取塔的水力學性能,分別建立了特征速度、分散相持液量以及Sauter液滴直徑的預測模型。

(4)Kühni 萃取塔 Buchbender 等[21]將Kühni 萃取塔每個艙室分為三部分(葉輪與下部塔板的下部區、葉輪所在的中部區、葉輪與頂部塔板的上部區)對液滴的平均停留時間進行研究。結果表明:減少底部區可使萃取塔的高度降低20%而不會對液滴的停留時間有顯著影響。Rode等[22]利用現有的Sauter液滴直徑模型、液滴終端速度模型、液滴特征速度模型對液泛區域進行了成功預測。Gomes等[23]研究發現,含有溶質的液滴表面張力更小,破碎頻率更高,使液滴更小、停留時間更長、分散相持液量增加。Asadollahzadeh等[24-25]首先通過不同萃取體系、不同傳質方向進行研究,建立了Kühni萃取塔特征速度與滑移速度的預測模型,然后分別采用最大熵法、Gamma 法、逆高斯和Weibull 法對Kühni萃取塔的液滴直徑分布進行了預測,結果表明最大熵法預測效果最為準確。Arab等[26]發現由分散相到連續相傳質會使持液量增大,而連續相到分散相則相反,這可能是受Marangoni效應的影響。

對于其他攪拌萃取塔, Hemmati[27]、Asadollahzadeh[28]等分別采用三種不同的萃取體系對篩板轉盤塔、不對稱轉盤萃取塔的水力學進行研究,建立了分散相持液量與液泛速度的預測模型。Asadollahzadeh 等[29]建立了Kühni 萃取塔、Oldshue-Rushton萃取塔分散相持液量的統一預測公式。

1.2 軸向擴散和相間傳質

Stella等[3]對Karr往復板篩板萃取塔的傳質性能進行了研究,結果表明,傳質方向由連續相到分散相的傳質系數高于從分散相到連續相,這是因為由連續相到分散相的傳質過程中,液滴尺寸更小,提高了傳質面積;同時,基于Harikrishnan 等[30]的關聯式提出了傳質系數預測模型。Smith等[4]通過對模型中參數的修正,使該模型可以準確預測堿性體系,同時提出了Karr 往復板篩板萃取塔的放大準則。

Yadav 等[8]使用傳質模型對近50 年來不同脈沖篩板萃取塔的實驗數據進行了預測驗證,結果表明模型均不能達到很好的預測效果。Torab-Mostaedi等[31]通過三種萃取體系的數據,對Johnson等[32]的傳質模型中的增強因子進行了修正,使預測結果更加準確。Bahmanyar 等[33]利用Newman[34]的傳質模型得到有效擴散系數,并將有效擴散系數與Reynolds數相關聯,建立了有效擴散系數的預測模型,通過引入有效擴散系數更準確地預測了濃度隨塔高的變化規律。

Charton 等[35]對脈沖板環萃取塔的軸向擴散系數進行了深入研究,基于Buratti[36]所建立的模型,擴大了軸向擴散系數模型的預測范圍(脈沖強度<0.075m/s,0.17<板間距/塔徑<1.33),對于板間距/塔徑<1 的小型萃取塔,其預測誤差由143%降到了21%,對于工業型萃取塔,其預測誤差由95%降到了27%。Liu 等[17,37]通過50mm、70mm 與90mm內徑的萃取塔對此模型進行了驗證,預測效果良好,同時,在Kumar 與Hartland[38]模型的基礎上,提出了新的模型,實現了對不同塔徑板環萃取塔傳質系數的準確預測。Jahya 等[39]結合孔隙率提出了脈沖板環萃取塔的幾何放大方法,如式(1)~式(3)所示。

在實際的工業生產中,無脈沖的板環萃取塔也得到了廣泛的應用,Wang等[40-41]研究了板環萃取塔在有脈沖與無脈沖條件下的軸向擴散與傳質性能,并分別建立了統一的預測模型,對于軸向擴散與傳質系數的預測精度分別在20%與22%以下。

Safari 等[42]研究發現,采用更精確的軸向擴散模型計算的傳質系數低于采用平推流模型計算結果的20%,表明了軸向擴散對萃取塔性能影響的重要性。同時,發現液滴直徑的分布是影響填料脈沖萃取塔軸向擴散的重要原因,分別提出了分散相、連續相的軸向擴散系數預測模型與傳質系數預測模型。

Torab-Mostaedi 等[43]研究了非對稱轉盤萃取塔的傳質性能,結果表明,傳質系數受轉速的影響較大,而受連續相流速的影響較小;在傳質由連續相到分散相的過程中,分散相流速的影響對傳質系數影響較小,而相反的傳質方向影響則較大;結合Sherwood數、Reynolds 數和分散相持液量,建立了非對稱轉盤萃取塔的傳質系數預測模型。Asadollahzadeh 等[44]對該模型參數進行修正,成功地預測了不同傳質方向下Kühni 萃取塔的傳質系數,同時發現降低萃取體系的表面張力可有效地提高傳質性能。

2 萃取塔流體力學模擬

2.1 CFD模擬

隨著計算機技術迅速發展,計算流體力學模擬方法逐漸應用于萃取設備的流體力學預測中,為萃取塔設備的設計、放大提供了新的方法。

Bardin-Monnier[45]、Modes[46]等分別對脈沖板環萃取塔與轉盤塔中單液滴的流動行為進行了模擬,研究發現由于液滴-壁面的相互作用的影響,液滴在一個板間距內的停留時間呈現出多峰分布。Sarkar 等[47]對脈沖板環萃取塔中單相流場與軸向擴散進行模擬,結果表明隨著板間距增加和孔隙率減小,軸向擴散與塔內壓降逐漸減小,而隨著板環的外徑和內徑之比的增加,軸向擴散與塔內壓降逐漸增大。Tang等[48]采用CFD對三種不同結構的脈沖篩板萃取塔內的單相流場進行了數值模擬,得到了比標準篩板結構傳質性能更好的分散-聚并型篩板結構。

Duan[49]、Bujalski[50]等結合粒子圖像測速技術(PIV)分別對篩板萃取塔與脈沖板環萃取塔內單相流進行了模擬,結果表明實測流場與CFD 模擬結果吻合較好,Duan 等[49]發現篩板之間存在較大的逆流區,通過改變降液管結構可優化連續相的流動狀態。Amokrane等[15]則結合PIV對轉盤萃取塔內的單相流動進行模擬,比較了不同湍流模型對預測效果的影響(圖2),結果表明,如果是完全湍流,標準k-ε(或者RNG k-ε)湍流模型與一階迎風格式相結合可達到足夠的預測精度;當需考局部湍流時,應優先采用RNG k-ε紊流模型和乘方格式。

CFD 模擬研究也逐漸應用于萃取塔的液-液兩相流動體系。Yadav[51]、Din[52]、Saini[53]等分別基于雙歐拉模型,標準k-ε模型通過FLUENT 軟件模擬成功地預測了兩相逆流流動的脈沖篩板塔與脈沖板環萃取塔的持液量。Saini等[53]研究發現,液滴直徑對板環塔內持液量影響較小,但會影響分散相在塔內的分布,脈沖頻率對分散相持液量影響較大,而振幅影響較小。Retieb 等[54]利用ASTRID 軟件對脈沖板環萃取塔的持液量進行了成功的預測,誤差在12%以內。

圖3 CFD模擬值與PIV測量值的結果比較[55]

研究者們還采用CFD 對塔內的流場情況進行了研究。Drumm等[55]對轉盤塔內單相以及兩相的流體力學進行模擬,并利用PIV 對流場進行了測量(如圖3),模擬結果與測量結果吻合較好。Sarkar等[56-57]通過模擬考察了操作條件、幾何尺寸對脈沖板環萃取塔持液量的影響,同時引入運輸方程對軸向擴散進行了模擬,結果表明,持液量與軸向擴散的預測值和實驗值之間的誤差分別在8%與5%以內,具有良好的一致性。Yi等[58]對液滴的曳力系數進行修正,利用CFD 對無脈沖的板環萃取塔內軸向擴散進行了預測。考慮工業萃取塔內流體流動的非對稱性,Yu[59]、Laitinen[60]等分別對脈沖板環萃取塔與Kühni塔的流體力學進行了三維模擬,結果表明較小的塔徑不存在明顯的徑向分量,可用二維模擬代替。

2.2 CFD-PBM模擬

分散相中液滴尺寸的分布對兩相流動行為的準確預測起關鍵性作用,群體平衡模型(PBM)能夠描述離散相實體的分布特性以及引起分布變化的離散相微觀行為,對雙流體模型的湍流應力、相間作用力和相間傳質至關重要[61],將PBM 和CFD 相結合建立CFD-PBM耦合模型,可提高預測的準確性。

Sen 等[62]對Kumar-Hartland 阻力模型的參數進行了優化,采用CFD-PBM 耦合模擬使分散相持液量預測誤差減小到了5.6%;Hlawitschka等[63]基于雙膜理論,并采用CFD-PBM 耦合對轉盤塔的傳質過程進行模擬,較好地預測了溶質在塔內的濃度分布。Amokrane等[64]首先通過攪拌槽反應器中測量的液滴粒徑分布數據,對群體平衡方程(PBE)的破碎模型參數與進行了擬合,再通過CFD模擬方法,選取k-ε湍流模型研究了脈沖板環萃取塔內的單相流動,并用PIV 實驗進行了驗證,最后利用CFDPBM耦合模型,對脈沖板環萃取塔內Sauter液滴直徑進行了準確預測。Attarakih 等[65]提出了一種基于簡化二元群體平衡模型模的流體力學與傳質層次化模擬的方法,該方法利用一維CFD 模型和的二元種群平衡模型相耦合,在顆粒群平衡實驗室(PPBLAB)軟件中進行實現。

液滴破碎實驗數據的缺乏是群體平衡模型應用的局限性之一,Fang等[66]利用高速相機研究了脈沖板環萃取塔內親水性內構件對液滴破碎聚并行為的影響,提出了新的破碎頻率、聚并頻率與破碎子液滴尺寸分布模型。Korb 等[67]用二(2-乙基己基)磷酸與硫酸鋅體系,在內徑為32mm 的Kühni 萃取塔中進行了單液滴破碎實驗,提出了新的滴破碎概率和子液滴尺寸分布模型。Zhou等[68]測量了脈沖板環萃取塔內的液滴破碎行為,提出了新的破碎頻率、子液滴尺寸分布經驗模型。Yu 等[69]利用Zhou 等[68]提出的新的破碎模型,對脈沖板環萃取塔內的兩相流動進行模擬,研究了液滴粒徑范圍和分散相入口初始液滴粒徑分布的影響,模擬結果與實驗數據吻合較好。

2.3 其他模擬方法

基于簡化的群體平衡模型,液液萃取模塊(LLECMOD)、ReDrop 程序等相關模擬軟件因計算負荷低、速度快等優點,在不同類型萃取塔的模擬中有一定的應用。

Jaradat 等[70-71]利用LLECMOD 模擬了不同脈沖條件、流比下,填料萃取塔和篩板萃取塔的穩定性。結果表明:脈沖強度對高界面張力體系的影響顯著,而流量的變化對分散相持液量、液滴直徑和低界面張力體系的溶劑濃度剖面有很大影響。

Hlawitschka等[72-73]通過激光誘導熒光(LIF)與PIV分別測量兩相分布與速度場,同時使用LLECMOD與FLUENT 對Kühni 塔內的水力學進行了模擬,結果表明,LLECMOD 使用一維計算的預測準確度高;FLUENT對持液量的預測值稍低于實驗值,并需要較高的計算負荷,但可給出液滴、兩相分布、湍動能損耗和兩相流速等詳細信息。Buchbender等[74]采用單液滴實驗測量了Kühni 萃取塔內分散相持液量、Sauter平均直徑等水力學參數,并基于蒙特卡洛法,利用ReDrop 程序中的群體平衡模型對其水力學性能進行計算,與實驗數據吻合良好。

3 萃取塔的過程強化

3.1 新型萃取塔內構件

3.1.1 復合陶瓷內構件脈沖萃取塔

在鹽湖提鋰和稀土元素的分離過程中,氯化物對不銹鋼萃取塔內構件有很強的腐蝕性,Yi等[75-76]分別用陶瓷篩板與新型復合陶瓷內構件對篩板萃取塔(圖4)進行了改進。

圖4 復合陶瓷脈沖萃取塔[76]

研究表明,與陶瓷篩板內構件相比,新型復合陶瓷內構件的持液量升高50%,液滴直徑減小40%,軸向擴散系數降低50%,傳質單元高度降低50%,最優操作條件下傳質單元高度可低至0.2m,采用脈沖篩板萃取塔模型可以對陶瓷篩板與復合陶瓷內構件的水力學性能、軸向擴散以及傳質性能進行滿意預測。

圖5 Tenova脈沖萃取塔與內構件[77]

3.1.2 Tenova脈沖萃取塔

Li 等[77]對板環塔的內構件進行了改進(Tenova脈沖萃取塔),如圖5 所示,板的邊緣與環的內部為鋸齒狀波紋,以強化液滴在內構件邊緣的破碎,增加比表面積、提高傳質系數。

實驗結果表明,與脈沖板環萃取塔相比,Tenova 脈沖萃取塔的持液量更低、液滴直徑更小、軸向擴散系數更小[78-79]。將該萃取塔分別應用于萃取反應速率較快的體系(硫酸-叔胺-異癸醇,一級動力學反應速率為1.4s-1)和萃取反應速率相對較慢的體系(硫酸銅-1-(2-羥基-5-壬苯基)乙酮肟,一級動力學反應速率為2.1×10-2s-1),結果表明對于萃取反應速率較快的體系,Tenova脈沖萃取塔傳質效果與傳統板環塔一致;對于萃取速率較慢的體系,Tenova脈沖萃取塔可以更好地消除兩相流量對傳質效果的影響,增大傳質系數[80]。

3.1.3 新型填料萃取塔柱

張敏卿等[81]通過對板波結構紋填料的不同齒度的網孔(網孔3.5×5.9,網孔1.7×3.4,網孔1.6×2.6)研究發現,網孔越小,切割出的分散相液滴的直徑越小,分散相液滴群總表面積越大,傳質效率越高,與鮑爾環相比傳質單元高度低15%到45%。

為了解決化工溶劑脫瀝青,潤滑油精制,含酚廢水處理等高運行通量的分離過程,清華大學化工系在VKB 型規整填料的基礎上又開發了SBW 新型的舌形板波紋規整填料,在FG 系列蜂窩格柵規整填料的基礎上開發了新型DFG導向格柵規整填料。蔡衛濱等[82-83]分別對SBW 舌形板波紋規整填料與DFG導向格柵規整填料(圖6)萃取塔的水力學性能和傳質性能進行了系統研究,結果表明SBW 舌形板波紋規整填料萃取塔的傳質單元高度為0.6~1.4m,優于舌形格柵及蜂窩格柵等規整填料,液泛通量為90~150m3·m-2·h-1;DFG導向格柵規整填料萃取塔的傳質單元高度為2.2~3.1m,液泛通量為123~154m3·m-2·h-1。

圖6 2種新型填料的照片

在化工分離領域,磁場[84]、電磁場[85-86]強化具有清潔、高效等優點,是最具發展潛力的過程強化技術之一。Amani 等[87]考察了磁性納米顆粒和磁場對傳質效果的影響。在規整填料中,用四個相同的鋅鐵氧磁芯體棒產生磁場(圖7),在有機相(分散相)中加入質量分數為0.001%~0.005%的尖晶石型錳鐵氧體納米粒子。結果表明在磁場的影響下,磁性納米粒子布朗運動加強,使得顆粒與周圍媒介相對速度增大,提高了甲苯-乙酸水體系的傳質效果。

3.2 新型萃取塔型式

3.2.1 攪拌-脈沖萃取塔

近些年,Kockmann 團隊[88]對不同類型萃取塔的效率進行了總結,將攪拌萃取塔與脈沖萃取塔相結合,設計了內徑為15mm 的攪拌-脈沖萃取塔(APC),如圖8所示。葉輪通過高速攪拌增強液滴的破碎,由底部產生的脈沖使得分散相和連續相可以順利地逆流通過板間隙,液滴在塔板處聚結并進入下一個單元艙室,保證液滴的破碎聚并區域不斷更新。

圖7 填料塔內四極磁場圖[87]

圖8 攪拌-脈沖萃取塔[89]

通過實驗分析,在高轉速下(>900r/min)該萃取塔內分散相持液量接近50%,每米傳質單元數高達17~25級。該萃取塔已被成功應用于芳香酸、3,5-二硝基苯甲酰-(R,S)-亮氨酸等手性物質的分離與5-羥甲基糠醛的富集,收率高、分離效果好[89-90]。

Soboll 團隊[91-94]對直徑為32mm 的攪拌-脈沖萃取塔的水力學、傳質單元高度進行了研究,并與直徑為15mm 的同型攪拌-脈沖萃取塔進行比較,以研究該強化萃取塔的放大效應。結果表明:相比于直徑15mm 的萃取塔,將直徑擴大到32mm 后,液泛通量增大,每米傳質單元數有所下降,但仍是傳統攪拌萃取塔的兩倍。對于該類型萃取塔的數學模型還比較缺乏,本文作者團隊目前正在結合模型方法與CFD 模擬對該萃取塔的性質進行系統研究,為工業設計放大與應用提供基礎。

3.2.2 L-型脈沖萃取塔

Amani 等[95-96]設計了一種L-型(垂直-水平)脈沖篩板萃取塔,并研究了不同物系、脈沖強度對萃取塔中能量耗散、特征速度、軸向擴散與傳質系數的影響,建立了預測公式。Mohammadi[97]、Panahinia[98]等分別研究了萃取塔的持液量、不同傳質方向對軸向擴散與傳質系數的影響,并建立了相應的預測模型;Rafie 等[99-100]將內構件改為填料,考察了不同操作條件對軸向返混與傳質性能的影響。

4 結語

萃取塔作為溶劑萃取的一種重要單元設備被廣泛應用。目前萃取塔的研究主要是圍繞著液-液兩相的流體力學、液滴群行為、傳質過程模擬、設備強化等方面開展。近年來,針對不同類型的萃取塔,研究者們相繼建立了更完善、準確的經驗數學模型,廣泛應用于工程領域。計算機模擬方法的出現為萃取塔的研究提供了為更豐富的理論基礎,適用性較寬,可為萃取塔內復雜流場的研究提供支持,也為萃取塔的結構優化、工業放大提供可靠保證。

符號說明

da,dd,d—— 分別為圓環內徑、圓盤直徑、圓環外徑,m

hc—— 連續相表觀傳質單元高度,m

ε—— 孔隙率

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产免费羞羞视频| 久久精品这里只有国产中文精品| 性色一区| 亚洲欧洲日韩久久狠狠爱| 色综合国产| 在线无码九区| 国产va免费精品| 3p叠罗汉国产精品久久| 综合网久久| 五月丁香伊人啪啪手机免费观看| 国产又大又粗又猛又爽的视频| 国产女同自拍视频| 欧美在线一二区| 久久无码av三级| 亚洲电影天堂在线国语对白| 免费三A级毛片视频| 色天堂无毒不卡| 亚洲无码一区在线观看| 野花国产精品入口| 国产日韩精品欧美一区灰| 亚洲成人网在线播放| 福利在线不卡| V一区无码内射国产| 高清久久精品亚洲日韩Av| 尤物国产在线| 国产综合亚洲欧洲区精品无码| 精品一区二区三区中文字幕| 亚洲黄网在线| 中国毛片网| 中国一级特黄视频| 亚洲综合18p| a欧美在线| 国产色爱av资源综合区| 国产自在线播放| 国产黑丝视频在线观看| 国产理论精品| 天天摸天天操免费播放小视频| 99999久久久久久亚洲| 99久久国产精品无码| 成年女人a毛片免费视频| 98超碰在线观看| 国产欧美视频综合二区| 亚洲日韩精品欧美中文字幕| 内射人妻无码色AV天堂| 偷拍久久网| 国产亚洲欧美在线视频| 99热这里都是国产精品| 国产高清无码第一十页在线观看| 91免费片| 亚洲久悠悠色悠在线播放| 玖玖免费视频在线观看| 丁香婷婷久久| 无码中字出轨中文人妻中文中| 国产成人a在线观看视频| 国产一区二区三区免费观看| 国产91线观看| 国产成人啪视频一区二区三区| 四虎精品免费久久| 欧美精品H在线播放| 欧美成人精品欧美一级乱黄| 国产主播喷水| 四虎成人在线视频| 2020最新国产精品视频| 日韩小视频在线观看| 亚洲高清在线天堂精品| 欧美国产在线精品17p| 欧美精品另类| 久久影院一区二区h| 精品夜恋影院亚洲欧洲| 国产不卡在线看| 最近最新中文字幕在线第一页| 超薄丝袜足j国产在线视频| 国产区人妖精品人妖精品视频| 伊人久久青草青青综合| 国产资源免费观看| 婷婷伊人久久| 欧美色99| 国产精品久久自在自线观看| 中文字幕在线观| 麻豆国产精品一二三在线观看 | 香蕉在线视频网站| 国产在线观看99|