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

多股進料攪拌釜停留時間分布的信息熵

2015-03-19 01:57:47王宇良鄭???/span>周業豐黃正梁王靖岱陽永榮
浙江大學學報(工學版) 2015年3期
關鍵詞:實驗

王宇良,鄭??。軜I豐,黃正梁,王靖岱,陽永榮

(浙江大學 化學工程聯合國家重點實驗室,化學工程與生物工程學系,浙江 杭州310027)

目前關于不同裝置內流動與混合情況的研究多是針對單股進料和單股出料系統,然而實際化學反應器、混合器、燃燒室和精餾塔等常具有多股進料或多股出料[1].此外,生理系統[2]以及濕地的水流匯合和分流過程[3]中也存在類似情況.

對于具有多股進料或出料的系統,裝置內的混合質量尤為重要.以具有多股進料的高壓聚乙烯釜式反應器為例,如果混合不良,局部引發劑濃度過高會導致產品質量下降,嚴重時甚至發生爆聚,造成停車.該反應器長徑比很大,內部結構復雜,并且目前缺乏對其中流體流動、混合情況的充分認識,這都限制了反應器的優化與放大.因此,對具有多股進料或出料的復雜系統的內部流動混合情況進行研究具有重要的實際意義.

停留時間分布可以描述系統內的流動特性,是微觀混合過程在宏觀上的表現.Nauman[4]綜述了停留時間理論的發展歷程、相關研究成果以及在不同領域的應用情況,認為停留時間理論雖然成熟但是并非停滯不前.Leclerc等[5]擴展了停留時間分布理論,通過由基本反應器相互連接的復雜網絡來理解停留時間分布,并開發了軟件包.該軟件可以用于模型的選擇和模型參數的檢驗,多股進料和出料過程可以通過卷積和最優化進行模擬.

停留時間分布常采用示蹤響應法進行測定,由實驗數據處理得到停留時間分布密度函數E(t).為了比較不同條件下測得的停留時間分布,通常是比較其統計特征值,如數學期望、方差.將這些特征值代入多釜串聯模型和軸向混合模型求取模型參數,以此反映流體的流動混合情況.但是,這些模型都是根據單股進料和單股出料的系統推導得出,應用于具有多股進料或出料的系統時存在困難.具有多股進料或出料的裝置,往往具有復雜的結構,不同的進料口與出料口的組合,得到的停留時間分布可能存在巨大差異.在某些組合下,可能導致裝置內流體存在滯留區,由停留時間分布曲線得到的無因次方差甚至會大于1[4].可見,將傳統的方差分析應用于多股進料或出料的系統時,可能得到片面甚至錯誤的結論.

信息熵(information entropy)可以用于描述混合效果,因為大多數的混合現象都可以利用概率進行解釋.信息熵由信息論之父Shannon[6]于1948年首次提出,故又稱Shannon熵(香農熵).信息熵是信息源紊亂程度的量度,反映系統運動狀態的隨機性.近年來,不少研究者將信息熵概念應用到物理、化學、生物及工程等領域中,取得了豐碩的成果[7-15].

在化學工程領域,信息熵已經被成功應用于鼓泡塔、流化床和氣力輸送等過程的復雜流動信號的分析,表明信息熵能較好地反映流體的流動特征.Manish等[11]采用信息熵理論分析了下流式鼓泡塔的示蹤劑濃度變化數據,以一個無因次參數描述裝置內混合情況,并將該參數與傳質效率進行關聯.Gui等[12]采用離散元耦合大渦模擬方法(discrete element method-large eddy simulation,DEMLES)對鼓泡流化床進行計算流體力學(computational fluid dynamics,CFD)模擬,并對顆粒的交叉徑向分布函數進行信息熵分析,研究發現平均信息熵可以很好地表征不同操作條件下的宏觀混合情況.Nedeltchev等[13]基于信息熵理論和信號劃分區域分析,提出了多相流反應器(鼓泡塔,流化床,噴射床)中不同流型辨識的新方法.Rakoczy等[14]通過停留時間分布實驗研究了旋轉磁場內一種混合器的混合特征,采用信息熵理論對停留時間分布數據進行處理,提出信息混合能力的概念用于描述在旋轉磁場下的連續流動系統的混合質量.鐘文琪等[15]基于壓力脈動信號的信息熵分析,建立了信息熵與噴動流化床中流型之間的聯系.

本文針對多股進料復雜流動系統,以攪拌釜為例,采用脈沖示蹤實驗,考察不同示蹤劑進料位置和不同操作條件下的停留時間分布,并結合信息熵理論方法對數據進行處理與分析,提出信息熵混合度的概念用于表征流動與混合情況,以期揭示攪拌釜裝置內的流動特性.此外,通過實驗方法驗證多股進料條件下總停留時間分布的疊加規則.

1 實驗裝置及方法

1.1 實驗裝置

實驗裝置如圖1所示.采用高為2 150 mm、內徑為210 mm的有機玻璃攪拌槽作為冷模實驗裝置,該裝置有4個進料口,1個底部出料口,裝置相關參數見表1.表1中,di為筒體內徑,ds為攪拌軸直徑;h為筒體高度;hb為中部軸承套高度;h1、h2、h3及h4為4個進料口(下文中以inlet1、inlet2、intel3及inlet4表示)的高度.攪拌軸上共裝有32層攪拌槳,其中29層為直徑190 mm的雙葉斜槳,槳面與運動方向的傾斜角為10°.此外在裝置頂部與中部分別裝有直葉渦輪,裝置底部安裝一個類似錨式的特殊槳.在中部直葉渦輪上方有一軸承套,其中軸承起到固定攪拌軸的作用,軸承套和攪拌軸的截面積占據了裝置截面積的75%左右.

圖1 具有4股進料和單股出料攪拌釜的停留時間分布實驗示意圖Fig.1 Schematic diagram of RTD experiments in stirred tank with four inlets and one outlet

表1 攪拌釜實驗裝置參數Tab.1 Parameters of experimental apparatus of stirred tank mm__

采用自來水作為實驗介質,攪拌轉速n通過變頻器進行控制.裝置內液面高度為2.0 m,液體體積為67.0 L.工業裝置在穩定運行時,inlet1的流量為其他進料口的一半,故取inlet1流量為0.5 m3/h,取inlet2、inlet3及inlet4流量為1.0 m3/h.

實驗以飽和KCl溶液作為示蹤劑,采用脈沖法測量停留時間t分布,使用電導率儀測定裝置出口處液體的電導率.t=0時,通過進料口旁的支路以注射器脈沖注入40 mL示蹤劑,同時開始記錄出口處的電導率c(t).實驗所用自來水的電導率約為170μs/cm,為了消除自來水電導率對實驗結果的影響,對c(t)進行歸一化處理,得到無因次的電導率c′(t):

式中:cmin、cmax分別為單次實驗過程中電導率的最小值與最大值.由c′(t)可以計算得到相應的停留時間分布密度函數:

1.2 分析方法

1.2.1 方差分析 采用方差分析[16]求E(t)的數學期望,得到平均停留時間

求E(t)對ˉt的二次矩,可得方差

為了便于比較不同條件下的停留時間分布,定義無因次時間θ:

則無因次停留時間分布密度函數E(θ):

已知對于平推流有

1.2.2 信息熵分析 通常,一個隨機變量的信息熵根據其概率分布定義.設X是一個離散的隨機變量,共有N種取值,其取值為χi(i=1,2,…,N)的概率p i=p(χi),則有

χi事件所對應的自信息量定義為式中:對數底a通常取2、e或10,對應信息熵的單位分別為bit、nat、det,三者可以互相換算.下文的推導中會涉及e,故本研究的信息熵中對數底均取e.

信息熵H定義為信源的加權平均信息量,即

規定當p(χi)=0時,-p(χi)loga p(χi)=0.確定概率事件的信息熵最小,完全隨機事件的信息熵最大.

概率密度為p(χ)的連續隨機變量有

式中:Δt為時間步長.

在進行停留時間分布實驗時,出口處示蹤劑濃度降到0的時刻為tm,則式(16)中N≥m,m表示離散隨機變量,有m種取值.

將信息熵理論應用于停留時間分布的分析是基于示蹤劑到達出口所需時間的不確定性.為了便于比較,可對信息熵進行歸一化處理[14,17],定義一個無因次參數M——信息熵混合度,用以表示流動情況接近平推流或者全混流的程度:

式中:Hmin為信息熵的最小值,Hmax為信息熵的最大值.對于平推流,將式(9)代入式(15),可得

對于全混流,將式(10)代入式(15)可得

對于m個進料口、單個出料口的閉式系統,應用全混流模型,可知

式中:V為裝置容積,q為進料總體積,等于各股進料流量之和.

將式(16)、(18)~(20)代入式(17),可得對于多股進料、單股出料的閉式系統,信息熵混合度為

式中:M為以信息熵表示的混合程度.數值介于0~1.0,其中0表示平推流,1.0表示全混流.顯然,M越大,則E(t)越接近于全混流狀態,混合質量越好.

2 結果與討論

當裝置有多股進料時,不但要掌握總停留時間分布,還要掌握各個進料口對應的停留時間分布.后者往往對裝置的設計與優化具有重要的指導意義.

2.1 4股進料單股示蹤時停留時間分布的方差分析

為了研究不同位置進料所對應的停留時間分布,保持進口1、進口2、進口3和進口4的流量不變,只在其中一個進料口處脈沖注入飽和KCl溶液作為示蹤劑.實驗考察了示蹤劑進料位置(進口1~進口4)、攪拌轉速n(300~700 r/min)對停留時間分布的影響,結果如圖2所示.

根據1.2.1節的方差分析方法,求得平均停留時間與無因次方差σ2θ,結果如表2所示.

表2 不同攪拌轉速和示蹤劑進料位置下的平均停留時間與無因次方差___Tab.2 Influences of agitation speed and tracer inlet position on mean residence time and dimensionless variance___

由表2可以發現,從進口1到進口4,越靠近裝置出口,則平均停留時間越短,說明示蹤劑在軸向的返混較小.尤其對進口4,其平均停留時間遠低于其他示蹤劑進料點,說明示蹤劑主要分布在裝置下部,裝置上部相當于滯留區,示蹤劑未能在整個裝置內實現良好混合.然而從進口1到進口4,無因次方差卻逐漸增大,表示越接近全混流狀態,顯然這與實際情況矛盾.由此可見,無因次方差并不適用于這類多股進料條件下的停留時間分布數據的分析,不能很好地反映不同進料位置、不同操作條件下的混合情況.下面將采用信息熵的方法來對停留時間分布數據進行處理和分析.

圖2 不同攪拌轉速和示蹤劑進料位置下的停留時間分布實驗結果Fig.2 RTD experimental results under different agitation speeds and tracer inlet positions

2.2 4股進料單股示蹤時停留時間分布的信息熵分析

M與攪拌雷諾數及示蹤劑進料位置的關系如圖3所示.從圖3可以發現,M與攪拌雷諾數Re及示蹤劑進料位置都有密切關系,但是兩者的影響效果有所不同.相同的示蹤劑進料位置下,隨著Re增大,液體湍流程度加強,示蹤劑在裝置內分布更加均勻,M隨之增大.在相同的攪拌雷諾數下,示蹤劑進料位置離裝置出口越遠,則示蹤劑在整個裝置內的分布越均勻,M越大.示蹤劑進料位置為進口1與進口2的M很接近,這是由于進口1與進口2均位于裝置頂部,與出料口之間的距離很接近,如圖1所示.

圖3 攪拌雷諾數和示蹤劑進料位置對信息熵混合度的影響Fig.3 Influences of Reynolds number and tracer inlet position on M

圖3中進口4的M值遠小于其他進料點,造成這種差異主要有2方面原因.一方面,裝置出料口位于底部,且進料、出料流量不變,則液體流動的平均速度一定且方向往下,導致位于裝置頂部的進口1~進口2注入的示蹤劑的平均停留時間較長.而進口4距離裝置出口較近,平均停留時間較短.另一方面,裝置中部的軸承套和攪拌軸的截面積占據了裝置截面積的75%左右,軸承套起到分區擋板的作用,因此在進口4注入的示蹤劑難以返混至裝置上部區域,上部區域相當于滯留區.總之,進口1和進口2中注入的示蹤劑可以隨主體流動通過軸承套與裝置內壁的間隙進入下部區域,在整個裝置內實現較均勻的分散.進口4注入的示蹤劑主要停留在下部區域,較快地離開裝置,難以到達裝置上部區域.

此外,存在一個有趣的現象:當Re較小時,進口3的M明顯小于進口2的M;而當Re足夠大時,兩者的M值非常接近.攪拌槽內存在3種擴散作用:主體對流、漩渦擴散以及分子擴散.inlet3位于裝置中部位置,且裝置內液體流動的平均速度往下.因此,攪拌雷諾數較小時,流體湍動程度較低,示蹤劑主要隨主體流動往下擴散,較快地進入下部區域,未能在上部區域充分擴散.示蹤劑在整個裝置內分布不夠均勻,因此M較小.而當攪拌雷諾數增大時,3種擴散作用均有所增強,軸向返混加劇,示蹤劑能同時在上下兩區域實現較好地擴散.

對比方差分析和信息熵分析的結果可知,對于多股進料的復雜流動系統,停留時間分布的無因次方差并不能準確地反映裝置內的流動與混合情況,在某些情況下甚至得到錯誤的結論.而信息熵混合度M可以準確反映流體的流動與混合特性.這是因為無因次方差的計算需要已知平均停留時間,而由停留時間分布數據計算得到的平均停留時間與裝置的結構、進出料的情況相關,往往偏離理論值(V/q).由此得到的無因次方差只能反映該曲線的性質而不能準確反映流動與混合的情況.而式(17)明確了M的意義為停留時間分布曲線的信息熵H接近全混流Hmax或平推流Hmin的程度.由M值的大小可以初步判斷流體的混合程度:M值越接近于1,則越接近全混流;M值越接近于0,則越接近平推流.

2.3 總停留時間分布的方差分析和信息熵分析

通過實驗可測得總停留時間分布,方法如下:在t=0時刻,進口1~進口4同時脈沖注入示蹤劑(10~20 m L),同時記錄裝置出口處液體的電導率變化.由于inlet1流量為其他進料口的1/2,為了保持各進料口示蹤劑脈沖濃度一致,inlet1加入的示蹤劑體積為其他進料口的1/2.實驗結果如圖4所示.對數據分別進行方差分析和信息熵分析,結果如表3所示.

從圖4可以看出,隨著攪拌轉速增大,停留時間分布曲線的峰變窄,拖尾現象加劇.由表3可知,隨著攪拌轉速的增大,由方差分析得到,平均停留時間與無因次方差均呈增大的趨勢,由信息熵分析得到的M值也呈增大的趨勢.

圖4 不同攪拌轉速下的總停留時間分布實驗結果Fig.4 Overall RTD experimental results under different agitation speeds

攪拌轉速增大,則攪拌雷諾數Re隨之增大,裝置內流體的湍流程度加強,返混更加劇烈,混合更加接近于全混流,所以平均停留時間增大,無因次方差增大.同樣,E(t)紊亂程度增大,M呈增大的趨勢,更趨近于1.

實驗裝置攪拌器共安裝多層攪拌槳,層間距低于攪拌槳直徑.旋轉攪拌槳能同時產生徑向流和軸向流,混合效率較高.且實驗條件下Re>180 000,液體處于完全湍流狀態[18].多股注入的示蹤劑能快速均勻地分散到整個裝置內.這與全混流的情況相似,故認為其接近于全混流.而由表3的計算結果可知,信息熵分析得到的M接近于1,說明攪拌釜接近于全混流,信息熵分析結果與實際情況相符合.

綜上所述,信息熵理論可應用于停留時間分布的分析,對于多股進料的復雜流動系統,與無因次方差相比,M能夠更加準確地反映流體的流動與混合特性.在比較不同進料口的流體在整個裝置中的流動與混合情況時,M值的優勢更加明顯.

2.4 總停留時間分布疊加規則的信息熵驗證

針對多股進料和出料的復雜流動系統,不少研究者通過理論推導與實驗方法得到總停留時間分布的疊加規則.Buffham等[19]最早對此進行研究,認為具有多股進料和出料系統的總停留時間分布密度函數遵守如下疊加規則:

式中:q是總體積流量,qej是第j股出料的流量,E kj是第j個出料口因第k個進料口注入的示蹤劑所產生的停留時間分布密度函數.Ritchie等[20]提出了不同的疊加規則:

式中:qkj是第k個進料口的進料中從第j個出料口流出的流量

莊震萬[1]采用多釜串聯模型,通過理論推導,建立了多股進料多股出料系統在理想、非理想流動情況下的總停留時間分布的疊加規則,結果與式(22)相同.對于本研究中的4股進料單股出料系統,根據式(22)或(23)均可得到相同的總停留時間分布疊加規則:

式中:Ek為第k個進料口注入的示蹤劑所對應的停留時間分布.

對圖2的實驗數據按式(24)進行處理,得到總停留時間ts分布的計算值.實驗測得的總停留時間分布曲線如圖4所示.將計算值與實驗值進行對比,典型曲線如圖5所示.分別計算實驗值與計算值的信息熵混合度,結果見表4.

圖5 不同攪拌轉速下總停留時間分布計算值與實驗值的對比Fig.5 Comparison of calculated values with experimental data of overall RTD under different agitation speed

表4 總停留時間分布的信息熵混合度Tab.4 Information entropy mixing degree of overall RTD

從圖5及表4可以發現,按疊加規則計算得到的總停留時間分布與實驗值吻合良好.通過信息熵分析可知,信息熵混合度的計算值(Mcal)與實驗值(Mexp)的相對偏差δ在2%以內,驗證了疊加規則式(18)的正確性.Mexp略小于Mcal,這是由于實驗的系統誤差,4個進料口處的示蹤劑脈沖并非完全相同,且脈沖進入裝置內的時間有所差別.inlet1示蹤劑注入點到裝置的距離較遠,故實驗時產生的峰延后.從圖5可以發現,實驗的E(t)峰值低于計算值;峰值過后,隨著t增大,兩者差距不斷縮小,直到實驗的E(t)超過計算值.可見,與計算值相比,實驗的E(t)分布較均勻,則信息熵較小,這與方差分析類似.

實驗和計算得到的M值均大于0.95,說明在整個裝置內流體流動情況接近于全混流;在實驗所采用的攪拌雷諾數范圍內,流體處于完全湍流狀態;隨著攪拌雷諾數增大,流體混合效果增強,M值有所增大.

3 結 論

(1)對于多股進料的復雜流動系統,與無因次方差相比,M能夠更加準確地反映流體的流動與混合特性.M介于0~1.0,其中0代表平推流,1.0代表全混流.

(2)M隨著攪拌雷諾數的增大而增大,隨示蹤劑進料位置與裝置出口距離的降低而減小.裝置中部的軸承套起到分區擋板的作用,上部分區物料容易進入下部分區,但下部分區物料難以返混至上部分區.

(3)對實驗測得的總停留時間分布和由疊加規則計算得到的總停留時間分布分別進行信息熵分析,驗證了疊加規則的正確性.M值均大于0.95,說明裝置內整體混合情況良好,接近于全混流.

):

[1]莊震萬.具有多股進料和出料的系統的流動模型[J].中國科學:A 輯,1980,10:12.ZHUANG Zhen-wan.Flow model of equipment with multiply inlets and oulets[J].Science in China:Series A,1980,10:12.

[2]江暉.氨芐西林與舒巴坦匹酯經口服在肉雞體內的藥物動力學及組織分布研究[D].武漢:.華中農業大學,2012:7- 20.JIANG Hui.Pharmacokinetics and tissue distribution of ampicillin and sulbactam pivoxyl following oral administration in chickens[D].Wuhan:Huazhong Agricultural University,2012:7- 20.

[3]WANG H,JAWITZ J W.Hydraulic analysis of cellnetwork treatment wetlands[J].Journal of Hydrology,2006,330(3/4):721- 734.

[4]NAUMAN E B.Residence time theory[J].Industrial and Engineering Chemistry Research,2008,47(10):3752- 3766.

[5]LECLERC J P,CLAUDEL S,LINTZ H G,et al.Theoretical interpretation of residence-time distribution measurements in industrial processes[J].Oil and Gas Science and Technology,2000,55(2):159- 169.

[6]SHANNON C E,A mathematical theory of communication [J].Bell System Technical Journal,1948(27):379- 423,623- 656.

[7]CHATZISAVVAS K C,MOUSTAKIDIS C C,PANOS C P.Information entropy,information distances,and complexity in atoms[J].Journal of Chemical Physics,2005,123(17):174111.

[8]STAHURA F L,GODDEN J W,XUE L,et al.Distin-guishing between natural products and synthetic molecules by descriptor Shannon entropy analysis and binary QSAR calculations[J].Journal of Chemical Information and Computer Sciences,2000,40(5):1245- 1252.

[9]STRAIT B J,DEWEY T G.The Shannon information entropy of protein sequences[J].Biophysical Journal,1996,71(1):148- 155.

[10]LI H,BAO Y Q,OU J P.Structural damage identification based on integration of information fusion and shannon entropy[J].Mechanical Systems and Signal Processing,2008,22(6):1427- 1440.

[11]MANISH P,MAJUMDER S K.Quality of mixing in a downflow bubble column based on information entropy theory[J].Chemical Engineering Science,2009,64(8):1798- 1805.

[12]GUI N,FAN J.Numerical study of particle mixing in bubbling fluidized beds based on fractal and entropy analysis[J].Chemical Engineering Science,2011,66(12):2788- 2797.

[13]NEDELTCHEV S,SHAIKH A.A new method for identification of the main transition velocities in multiphase reactors based on information entropy theory[J].Chemical Engineering Science,2013,100:2- 14.

[14]RAKOCZY R L,KORDAS M,STORY G,et al.The characterization of the residence time distribution in a magnetic mixer by means of the information entropy[J].Chemical Engineering Science,2013,105:191- 197.

[15]鐘文琪,章名耀.基于信息熵分析的噴動流化床流動特性[J].化工學報,2005,56(12):2303- 2308.ZHONG Wen-qi,ZHANG Ming-yao.Flow characteristics in spou-t fluid bed by Shannon entropy analysis[J].Journal of Chemical Industry and Engineering,2005,56(12):2303- 2308.

[16]陳甘棠.化學反應工程[M].北京:化學工業出版社,2007:88- 109.

[17]NEDELTCHEV S,OOKAWARA S,OGAWA K.A fundamental approach to bubble column scale-up based on quality of mixedness[J].Journal of Chemical Engineering of Japan,1999,32(4):431- 439.

[18]王凱,虞軍.化工設備設計全書:攪拌設備[M].北京:化學工業出版社,2003:17- 31.

[19]BUFFHAM B A,KROPHOLLER H W.The washout curve,residence-time distribution,andFcurve in tracer kinetics[J].Mathematical Biosciences,1970,6:179- 184.

[20]RITCHIE B W,TOBGY A H.Residence time analysis in systems having many connections with their environment[J].Industrial and Engineering Chemistry Fundamentals,1978,17(4):287- 291.

猜你喜歡
實驗
我做了一項小實驗
記住“三個字”,寫好小實驗
我做了一項小實驗
我做了一項小實驗
記一次有趣的實驗
有趣的實驗
小主人報(2022年4期)2022-08-09 08:52:06
微型實驗里看“燃燒”
做個怪怪長實驗
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
主站蜘蛛池模板: 1024你懂的国产精品| 国产精品永久久久久| 欧美有码在线| 国产91高跟丝袜| 亚洲综合网在线观看| 久久成人国产精品免费软件| 六月婷婷综合| 亚洲国产精品一区二区第一页免| 日韩在线网址| 91人妻在线视频| 免费在线看黄网址| 国产精品私拍在线爆乳| 亚洲欧美另类日本| 国产欧美精品一区二区| 久久6免费视频| 国产欧美视频综合二区| 国产精品高清国产三级囯产AV| 亚洲欧美日韩色图| 日本人真淫视频一区二区三区| 免费人成网站在线观看欧美| 国产精品高清国产三级囯产AV | 亚洲一级无毛片无码在线免费视频| 午夜免费视频网站| 国产96在线 | 亚洲精品第一页不卡| 国产拍揄自揄精品视频网站| 中文字幕中文字字幕码一二区| 亚洲国产综合精品一区| 福利一区在线| 亚洲国产成人综合精品2020| 国产精品亚洲一区二区在线观看| 久久综合干| 强奷白丝美女在线观看| 欧美一级爱操视频| 色老二精品视频在线观看| 国产国语一级毛片在线视频| 91亚洲视频下载| 在线观看精品自拍视频| 国产91av在线| 青青操国产视频| 欧美一区二区啪啪| 久久久久无码精品国产免费| 精品国产成人a在线观看| 精品国产黑色丝袜高跟鞋 | 亚洲丝袜中文字幕| 国产精品手机在线观看你懂的| 国产日本视频91| 国产美女精品人人做人人爽| 亚洲成人免费看| 黄色网站在线观看无码| 国产麻豆91网在线看| 亚洲人成成无码网WWW| 国产69精品久久久久妇女| 国产毛片久久国产| 久久精品只有这里有| 91破解版在线亚洲| 一本视频精品中文字幕| 无码精品福利一区二区三区| 国产乱子伦精品视频| 国产白浆视频| 91青青在线视频| 日韩欧美综合在线制服| 日本免费精品| 在线a网站| 毛片基地视频| 精品色综合| 午夜福利在线观看成人| 在线观看欧美国产| 久久动漫精品| 成人午夜天| 亚洲午夜天堂| 一边摸一边做爽的视频17国产| 欧美啪啪一区| 五月天久久综合国产一区二区| 国产精品私拍在线爆乳| 精品丝袜美腿国产一区| 综合网天天| 日韩大乳视频中文字幕| 日本精品中文字幕在线不卡| 在线观看视频一区二区| 日本少妇又色又爽又高潮| 亚洲人成网站在线播放2019|