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

新型濕法除塵系統內氣液兩相流動的數值模擬

2020-09-23 09:30:28許浩潔王軍鋒王東保張偉姚江
化工進展 2020年9期

許浩潔,王軍鋒,王東保,張偉,姚江

(江蘇大學能源與動力工程學院,江蘇鎮江212013)

現代工業生產中,粉塵不僅會造成環境污染,威脅作業人員健康[1],甚至還會帶來安全隱患[2-3],因此,對工業粉塵的有效治理十分必要。目前,常見的工業除塵方式按降塵原理可分為過濾式、沉降式、噴淋式以及靜電式。其中傳統過濾式干法除塵不僅壓降大、能耗高,同時在管道內還會形成粉塵集聚,難以滿足工業生產要求,而以噴淋式為主的濕法除塵技術因具有能耗低、結構簡單、適用范圍廣及可同時去除多種有害氣體等優點而廣泛應用于煤礦開采、化工脫硫、金屬拋光等現代工業生產過程中[4]。同時已有研究表明,擋板繞流[5]結構廣泛存在于高爐煉鐵[6]、流體換熱[7-8]、旋風除塵及濕法脫硫等工業領域,其流動狀態發生改變的同時常伴有剪切、分離、回流等現象,其中渦旋的產生可一定程度上強化化學反應、傳熱傳質及多相混合,對過程強化有著重要影響。為克服現有技術的不足,前期工作主要結合濕式除塵及擋板繞流技術優勢,同時運用多相流理論知識,設計了一種帶有擋板結構的拋光粉塵濕法處理系統[9],實現了含塵氣流與液面、液膜及液滴等多種形式間的相互作用。

濕法除塵作為典型的氣、液、固三相流動,其降塵機理宏觀上表現為液滴與顆粒物間的相互混合,其微觀降塵機理主要包括慣性碰撞、直接攔截、布朗擴散及重力沉降等,捕集過程涉及十分復雜的多相流動問題,通過傳統實驗手段研究難度較大。近年來,隨著計算流體力學(CFD)技術的快速發展[10],數值模擬已逐漸發展成為低成本、高效率的流體力學問題預測方法[11-13],且較傳統實驗方法具有節約人力物力、良好的可重復性及可獲得實驗難以測量的物理量等優點。其與實驗方法相互補充,在科學研究和工程技術中發揮著越來越重要的作用,已成為不可或缺的研究手段。

目前,針對濕法除塵技術的數值模擬研究多集中于除塵效率影響因素的探究[14],包括液滴速度、液滴粒徑、顆粒物粒徑等。其中劉曉燕等[15]對液滴速度及尺寸影響捕集效率過程進行了系統分析,發現液滴對顆粒的捕集效率隨顆粒物粒徑及液滴速度的增大不斷提高,而隨著液滴粒徑的增大則略有降低。Ding等[16]綜合分析了液滴尺寸、液滴速度及流場結構對去除效率的影響,其數值模擬結果表明,當液滴粒徑為30μm,且液滴與顆粒間相對速度為1m/s 時,其相互作用時間得到延長,除塵效率最高。王翱等[17]基于脫硫塔內環境,建立了綜合考慮慣性、攔截、擴散、熱泳和擴散泳作用的單液滴捕集顆粒模型,分析了溫度、液滴直徑及顆粒物粒徑對捕集過程及效率的影響規律。同時已有研究表明[18],除塵裝置流場分布情況可影響噴淋液滴及顆粒物運動特性,改變液滴停留時間及逃逸率,進而決定整體除塵效率。其中,魯軻等[19]通過增設導流裝置改善了礦用濕式除塵器內部的風流分布,使其脫水效率從93.2%提高到了99.3%。周山明等[20]對噴淋脫硫塔結構進行優化,使得流場分布更加合理;李彩亭等[21-22]通過在濕法除塵設備內部增設斜板及傘罩裝置,有效地延長了液滴的停留時間,進而強化了氣液相間的傳熱傳質;吳振元等[23]則從液滴直徑、液滴初速度、空塔氣速、噴淋密度等角度研究了噴淋吸收塔內液滴的運動及分布規律。

而針對于擋板繞流的研究則多集中于強化傳熱領域[24-27],較少涉及氣力輸送、濕法除塵等常見的工業多相流動過程。因此,一種帶有擋板結構的新型濕法除塵設備有待設計與開發,以進一步分析完善擋板繞流在強化多相流動過程中的作用機制。對該新型濕法除塵系統實物裝置的前期試驗研究結果表明,擋板及噴淋液滴的加入可使系統除塵效率由30%提升至82%[28],其關鍵就在于擋板結構強化了噴淋液滴與含塵氣流間的相互作用。基于該論斷,本文綜合濕法除塵與擋板繞流技術優勢,針對所設計的帶有擋板結構的噴霧降塵裝置,通過數值模擬方法分析了不同擋板安放角下的氣相旋流場特性,同時探究了旋流場作用下噴淋液滴相運動特性隨入口風速、擋板安放角及粒徑尺寸的變化規律,并擬合推導了噴淋液滴逃逸率預測公式,為濕法除塵器的設計及性能提升提供了理論參考。

1 模型建立

1.1 物理模型

帶有擋板結構的噴霧降塵系統示意圖如圖1所示,主要由變頻風機、濾網、弓形擋板及噴頭組成,其中噴頭噴淋方向為斜向下與含塵氣流呈逆流接觸。系統運行時,右側含塵氣流在風送流場作用下進入左側集塵區(對應圖1 紅框區域),由下而上依次流經液滴噴淋區、擋板旋流區及濾網過濾區,在集塵區內實現顆粒物與液面、液膜及液滴等多種形式間的相互作用,最終從上方出口排出。其中旋流區域內液滴與顆粒物間的相互作用是去除PM10、PM2.5等細小顆粒物的關鍵,因此,為了深入探究該區域內噴淋液滴及顆粒物的運動情況,選取集塵區作為計算區域建立物理模型,其中兩塊弓形擋板分別于集塵區兩側壁面交替設置,定義擋板弦與所在壁面的夾角為安放角(θ),并設置安放角大小為可調控,本文中安放角共設置有30°、45°、60°、90°及120°五種布置方案,如圖2所示,詳細結構參數見表1。

圖1 一種帶有擋板結構的噴霧降塵系統示意圖

圖2 擋板安放角布置方案

表1 集塵區計算模型結構參數

集塵區三維計算模型如圖3 所示(240mm×400mm×942mm),與實驗模型為1∶1 比例。利用ICEM 軟件對計算區域進行六面體結構化網格劃分,并對近壁處及擋板附近進行網格加密,同時通過Z=0.4m 處中心截面豎直方向速度分量對網格數量進行了無關性驗證,以安放角θ=60°為例,結果如圖4 所示,可以發現隨著網格數量增加到675608 后,進一步增加網格數量至1537830,速度值變化量不超過5%。因此,綜合考慮計算精度與計算速度,確定最終網格數量為675608。

圖3 計算區域三維模型建立及網格劃分

圖4 網格無關性驗證結果

1.2 數學模型

因實際工況下壓力差及溫度差均較小,因此模擬中將連續相氣體和離散相液滴均視為不可壓縮流體。本文著重探究氣、液兩相流動問題,初步將含塵氣流作純凈空氣進行處理。經試算,當入口風速為2m/s時,其雷諾數Re已達到105數量級,故選用k-ε湍流模型模擬氣相湍流運動,同時,噴淋液滴體積對于整個裝置內的氣相體積分數遠小于10%,且碰撞形式多表現為離散相與壁面間的相互作用,不考慮液滴間的相互碰撞,因此采用離散相模型(DPM)進行噴淋液滴相的行為描述。

1.2.1 連續相數學模型

不可壓縮定常流動的質量守恒方程(連續性方程)見式(1)。

動量守恒方程見式(2)。

式中,ρ 為流體密度;t 為時間;ui、uj分別為i、j方向上的速度分量;xi、xj分別為i、j方向上的位移分量;p為靜壓;τij表示作用在垂直于i軸平面上沿j 方向的切應力(即應力張量);fi為作用在單位質量流體微團上的體積力(包括重力和外部體積力)。

1.2.2 離散相數學模型

僅考慮曳力和重力作用的離散相運動方程見式(3)。

式中,up為離散相運動速度;u 為連續相流體運動速度;fD(u-up)為單位質量曳力;ρp為離散相密度;g為重力加速度;ρ為流體密度。

1.3 計算方法及邊界條件

由于具有收斂性好、精確度高等優點,標準k-ε 模型已廣泛應用于模擬各類湍流運動[29-30]。1986年,Yakhot 和Orszag[31]在此基礎上考慮了渦流對湍流的影響進而提出了RNG k-ε 模型,該模型提高了對旋流流動的預測精度[32-33]。本文選用RNG k-ε 湍流模型,并采用不考慮顆粒間相互作用的DPM 模型對裝置內的噴淋液滴及固體顆粒等離散相進行描述,同時考慮重力效應。其中,由于噴淋液滴與含塵氣流呈逆流接觸,離散液滴相的存在會對氣相流場產生一定的擾動,因此,為提高數值計算的準確性與可信度,模擬過程中采用連續相與離散相間耦合計算方法,其具體計算步驟為先計算連續相流場到收斂,然后加入離散相顆粒,再重新對連續相進行計算,直到連續相與離散相的計算結果都不會因為繼續耦合而發生改變。對液滴相采用穩態的追蹤方式,以獲得其在擋板結構下的運動情況。

設定氣相介質為空氣,液滴介質為水。含塵氣流以一定速度從速度入口注入,為避免回流,設置出口為壓力出口,結合操作工況和模型尺寸估算,確定湍流強度為5%,水力直徑為0.16m。噴頭軸線方向為斜向下且與豎直方向呈45°,采用實心錐(solid cone)方式模擬加入噴淋液滴,設定同時加入流場的液滴粒徑相同,噴淋錐角固定為60°,噴淋初始速度為10m/s。十個噴頭沿集塵區長度方向均勻布置,其噴淋總量固定為600mL/min。邊界條件的設置過程中,對于氣相而言,壁面均采用標準無滑移壁面;由于發生撞擊后液滴會吸附于壁面形成液膜,因此設置壁面對液滴相的邊界條件為碰撞即捕捉。

1.4 試驗驗證

1.4.1 實驗方案

為了驗證數值模擬方法的可靠性,搭建了如圖5 所示的實驗平臺,集塵區由透明亞克力板制成。其中以高濃度煙霧為可視化介質,近似地通過煙霧流線來表征氣相流場流線,煙霧由風機從底部入口給入,在連續片光激光的照射下,通過相機捕捉獲得煙霧的流動過程。通過微型轉輪風速儀對中心平面上的Z方向上的速度分量進行實驗測量。預測結果整體趨勢吻合較好,其中由于測量儀器對氣相流場的干擾,導致右側近壁區內的實驗值略低于數值模擬值。通過以上對比可以看出,數值模擬結果基本反映了集塵區內的實際流動情況,驗證了該數值方法的可靠性,渦旋結構的實驗結果將在后續部分進行討論。

圖5 實驗裝置

圖6 渦內Z方向速度分量的實驗與數值模擬結果對比

1.4.2 誤差分析

本文涉及的實驗測量參數包括速度分量和渦旋結構可視化。為了保證實驗測量精度,對相同工況下中心界面上的Z軸方向速度分量進行了多次測量取均值處理,而轉輪風速儀作為侵入式測量儀器,對于流場的系統誤差難以避免。在渦旋結構可視化實驗過程中,對比發現不同Y值處截面的煙霧流動具有相似的規律,為保證與數值模擬結果的一致性,選取對裝置中心截面進行拍攝。

2 數值模擬結果及分析

2.1 渦旋結構

由于集塵區結構沿長度方向具有延伸性,故選取Y方向中心截面模擬結果對氣相流場特性進行分析。同時在數值模擬過程中發現,相較于擋板安放角,氣相場流動特征受入口風速影響較小,且測試過程中發現當入口風速為1m/s 時,含塵氣流與噴淋液滴接觸較充分。因此本文以v=1m/s 為主要工況開展氣相場的數值模擬工作。

圖7為不同擋板安放角對應的氣相流線圖,其中θ=0°表示未設置擋板工況。從圖7(a)可以觀察到,未設置擋板時大部分氣流均沿著遠離入口一側壁面快速上升至出口,僅在入口上方壁面處由于剪切作用產生了一定量的回流,整體流動狀態較為平穩;增設兩塊安放角為30°的擋板后,可以發現氣流在擋板安置點的近壁處出現了多個小尺度附壁渦旋;調整擋板安放角為45°時,可以發現相較于30°,兩塊擋板后方的渦旋尺度顯著增大,而附壁小尺度渦旋數量逐漸減少;隨著安放角進一步增大至90°,對應圖7(e),可以發現兩擋板下游渦旋充分發展,其作用范圍達到最大,渦旋流動逐漸成為集塵區域主流運動。由以上分析可知,弓形擋板的設置可誘導產生多個渦旋流動,顯著提升集塵區的氣相湍流強度,為噴淋液滴與含塵顆粒間的相互作用創造了有利條件。

圖8 為安放角等于60°時,不同時刻煙霧在兩擋板間流動情況的可視化實驗結果。從圖中可以清楚地觀察到,擋板作用下煙霧在裝置內形成了明顯的旋轉流動,與數值模擬所得渦旋結構較好吻合。

圖7 不同擋板安放角對應的氣相流線圖

圖8 氣相流場的可視化實驗結果

2.2 氣相速度分布

由圖9 氣相速度分布云圖可知,未設置擋板時,由底部入口進入的氣流基本保持與入口相同寬度快速上升至出口,氣相在整個裝置內充滿效果較差,即不能夠充分利用集塵區域內空間,存在大量低速“流動死區”(對應速度云圖中藍色區域)。增設擋板后,可以發現當安放角θ 為30°、45°時,“流動死區”仍大量存在,主要集中在擋板附近區域,渦旋內流動速度較低;調整安放角為60°、90°及120°時,可以明顯地觀察到,由于擋板的剪切、分離作用,氣流在流經兩塊擋板時,由于流道收縮,流速迅速增加,隨后在剪切作用下,分別在兩擋板之間及上部擋板與出口之間形成了兩個速度分布由外圈向中心區域呈逐漸減小趨勢的渦旋流動,此時,氣體較好地充滿了整個集塵區域內部,“流動死區”情況得到明顯改善。

圖9 不同擋板安放角下的速度分布云圖

為進一步分析渦旋區域內氣相流動情況,模擬獲得了安放角θ=60°工況下的渦旋區域氣相矢量圖,如圖10 所示。可以發現氣流行至上部擋板尾部時,在與豎直壁面碰撞后產生了分離,大部分氣流繼續沿著直流通道快速上升,少部分氣流在產生撞擊后以相對較低的速度緊貼壁面向下運動,進而形成了所觀察到的旋流現象。圖11 為擋板間渦旋中心位置處速度分布,可以發現,渦旋呈現明顯的不對稱性,渦旋中心偏向下方有擋板一側,且在渦旋中心出現速度極小值,基本為零,在兩側靠近壁面處產生速度極大值。同時可以觀察到,安放角為60°和120°兩工況時,由于具有相等的直流通道寬度,因此其二者速度分布較好吻合;而當θ=90°時,由于直流通道進一步變窄,渦內速度極大值相對于60°和120°兩種工況均有所提高,且渦旋中心進一步向右側偏移。

圖10 渦旋區域氣相矢量圖

圖11 渦旋中心位置高度處的速度分布

2.3 壓力分布情況

壓降是評價除塵設備的一個重要指標,圖12為不同安放角對應的集塵區內靜壓分布云圖。從圖中可以看出,由于擋板的阻礙作用,使得原本較為均勻[對應圖12(a)無擋板工況]的壓力分布趨向于復雜化,在局部形成顯著的壓力梯度,具體表現為第一塊擋板下方區域壓力迅速升高,且由于流道迅速擴張,氣流流經擋板后在擋板后方區域產生了明顯的低壓回流區,直流通道處的流體在壓力梯度作用下不斷向該區域涌入,進而形成流體的旋轉運動。

圖12 不同擋板安放角下的靜壓分布云圖

圖13 進、出口壓降隨安放角及入口風速的變化規律

2.4 噴淋液滴運動軌跡

為探究不同工況下噴淋液滴的運動規律,本文分別對粒徑為60μm、120μm和160μm的液滴的運動軌跡進行了DPM 數值模擬。圖14 為不同粒徑噴淋液滴的運動軌跡,從圖14(a)中可以發現未設置擋板時,粒徑為120μm的噴淋液滴以10m/s的初速度進入氣相后,其運動速度迅速向當地氣相流速衰減,并在以1m/s 的速度從底側邊進入的含塵氣流作用下,除部分在慣性作用下直接撞擊到與入口相對的壁面及底部液面,大部分液滴則緊貼著遠離入口一側壁面快速上升至出口并產生逃逸。該過程中含塵氣流與噴淋液滴接觸過程較為平緩,集塵區空間利用率較低。

圖14 不同粒徑下噴淋液滴的運動軌跡

通過前述對氣相旋流場的分析可知,該復雜旋流場的存在不僅可提高氣相湍流強度,同時由于曳力作用還能夠有效延長液滴運行軌跡以及提高液滴與擋板及豎直壁面間的碰撞概率以減少夾帶現象。通過對比可以發現,當液滴粒徑較小(60μm)時,由于液滴跟隨性較好,大部分液滴均沿著氣相渦旋外圈高速區迅速向上運動,僅有少量液滴在流經擋板時與擋板發生碰撞捕集,出口仍存在明顯的逃逸現象。隨著液滴粒徑增大至120μm,對于單個噴淋液滴,由于質量的增加使得慣性力在顆粒的運動方程中所占比重增大,因此與60μm 粒徑液滴運動軌跡相比,可以觀察到第一塊擋板下方液滴的分布數量明顯增多。同時,由于液滴的慣性力與受到的來自逆流含塵氣流的曳力逐漸趨于平衡,導致液滴運動速度變化減緩,同時液滴軌跡也變得更加復雜,且被夾帶至第二塊折流板上方的液滴數量極少,這有助于噴淋液滴與含塵氣流間的充分混合。而當液滴粒徑繼續增大至160μm 左右,慣性力已遠大于氣流所引起的流動阻力,此時絕大多數液滴將直接沿著噴淋軸線方向撞向壁面及水池液面。

2.5 噴淋液滴逃逸率及駐留時間表現

工業生產過程中,逃逸率是衡量除塵設備性能的一個重要指標,逃逸現象不僅會影響設備除塵效率,同時還會造成二次污染。模擬過程中,定義逃逸率為從上方出口離開的液滴數量與累計噴淋液滴總量的比值。圖15 為入口風速1m/s 時,不同擋板設置方案下的液滴逃逸率情況。從圖中可以觀察到,未設置擋板時液滴逃逸現象嚴重。保持入口風速不變,增設不同安放角的擋板后,可以發現噴淋液滴的逃逸率得到了明顯改善,其中當安放角θ=60°時,改善效果最為明顯,該工況下粒徑范圍在40~160μm 之間的液滴逃逸率均維持在5%以下;調整安放角為90°,當噴淋液滴粒徑大于80μm時,已基本消除逃逸現象。而對于安放角為120°及45°工況,其分別在粒徑增大至100μm 及140μm 時,才較好地改善逃逸現象。

圖15 不同工況下噴淋液滴逃逸率的變化規律

由以上分析可知,擋板結構對噴淋液滴逃逸率的改善效果隨安放角呈現如下趨勢:60°>90°>120°>45°。分析擬合多種工況下的數值模擬結果,得到了考慮擋板安放角θ及液滴粒徑Dp的逃逸率計算公式[式(4)]。

本文所討論的駐留時間(Te)指的是所有噴淋液滴在碰撞到壁面、弓形擋板或從出口逃逸之前的在裝置內的平均運行時間。圖16 為噴淋液滴駐留時間隨擋板安放角、入口風速及液滴粒徑的變化規律,由于當θ=60°、90°時旋流現象最為明顯,因此僅對這兩種擋板安置方案下的駐留時間進行了統計分析。從圖16 中可以發現,當入口風速較低(0.5m/s)時,駐留時間隨著液滴粒徑的增大呈單調遞減趨勢,此時液滴受到來自氣相的曳力作用較小,隨著液滴粒徑的不斷增大,其質量不斷增加,質量力逐漸主導液滴運動,使得被氣流攜帶至旋流區的液滴數量迅速減少,進而導致噴淋液滴平均駐留時間減少。當入口風速增大至1m/s 時,可以發現駐留時間隨著粒徑的增大呈現先增加后減少的趨勢,在100μm 前后達到極大值,此時液滴所受慣性力與曳力相平衡,其運動速度緩慢,運動軌跡復雜;當入口風速進一步增大到2m/s 時,液滴慣性力遠大于所受曳力,使得噴淋液滴更易在下方擋板附近處即與碰面產生碰撞吸附,因而從圖16 中可以觀察到該工況下液滴駐留時間始終保持較低水平,且隨液滴粒徑與擋板安放角的調整變化較小。對于實際運行過程中液滴粒徑集中范圍內(80~140μm),可以發現1m/s 工況下相同粒徑液滴的駐留時間表現要明顯優于其他工況,且當安放角θ=60°時表現達到最優工況。

圖16 不同工況下的噴淋液滴駐留時間表現

3 結論

選用RNG k-ε 湍流模型及DPM 離散相模型數值模擬了帶有擋板結構的新型濕法除塵系統內的氣相場分布特征及噴淋液滴運動情況,詳細分析了不同工況下噴淋液滴的逃逸率及駐留時間表現,得到以下結論。

(1)交替設置的弓形擋板在裝置內誘導產生了速度分布由外圈向中心呈逐漸減小趨勢的渦旋流動,其合理布置可不同程度改善低速“流動死區”。

(2)氣相渦旋的存在不僅改善了噴淋液滴的逃逸現象,同時還提高了其在裝置內的駐留時間。

(3)分析給出了液滴逃逸率的預測公式,且在v=1m/s、θ=60°工況下獲得了駐留時間最優表現,為實際生產過程提供了一定的借鑒與參考。

符號說明

Dp—— 噴淋液滴粒徑,μm

fD(u-up)—— 單位質量曳力,m/s2

fi—— 單位質量流體微團上的體積力,m/s2

H—— 集塵區高度,mm

L—— 集塵區寬度,mm

L1—— 初始擋板高度,mm

L2—— 擋板間距,mm

L3—— 集塵區出口段高度,mm

Li—— 集塵區入口寬度,mm

Lb—— 擋板弦長,mm

Lo—— 集塵區出口寬度,mm

p—— 靜壓,Pa

R—— 擋板曲率半徑,mm

Te—— 噴淋液滴駐留時間,s

u—— 連續相流體運動速度,m/s

ui—— i方向上的速度分量,m/s

uj—— j方向上的速度分量,m/s

up—— 離散相相對速度,m/s

v—— 氣相入口風速,m/s

xi—— i方向上的位移分量,m xj—— j方向上的位移分量,m

ηe—— 噴淋液滴逃逸率,%

θ—— 擋板安放角,(°)

ρ—— 流體密度,kg/m3

ρp—— 離散相密度,kg/m3

τij—— 作用在垂直于i 軸平面上沿j 方向的切應力(即應力張量),kg/(m?s2)

主站蜘蛛池模板: 色网站在线免费观看| 免费无码网站| 91精品国产一区| 黄网站欧美内射| 噜噜噜综合亚洲| 久久九九热视频| 亚洲日产2021三区在线| 国产精品国产主播在线观看| 亚洲精品第五页| 日本www色视频| 91小视频版在线观看www| 久久精品aⅴ无码中文字幕| 尤物精品视频一区二区三区| a在线观看免费| 极品私人尤物在线精品首页| 1级黄色毛片| 国产精品亚洲一区二区三区在线观看 | 青青青国产视频手机| 中国精品久久| 熟妇丰满人妻| 亚洲有无码中文网| 国产小视频a在线观看| 国产丝袜第一页| 麻豆精品国产自产在线| 国产99在线观看| 国产乱子伦一区二区=| 亚洲国产看片基地久久1024| 国产成人免费观看在线视频| 麻豆精选在线| 国产综合色在线视频播放线视| 在线国产综合一区二区三区| 国模沟沟一区二区三区| 欧美亚洲欧美区| 91网在线| 99久视频| 91精品日韩人妻无码久久| a级毛片在线免费| 色婷婷视频在线| 91九色最新地址| 久久黄色一级视频| 亚洲,国产,日韩,综合一区| 欧美黄网在线| 日韩在线第三页| 人与鲁专区| 久久久久久午夜精品| 欧亚日韩Av| 99久久人妻精品免费二区| 国产视频大全| 国产精品99在线观看| 亚洲三级色| 日韩二区三区| 九九香蕉视频| 欧美成一级| 亚洲91精品视频| 亚洲女同一区二区| 欧美激情二区三区| 国产区网址| 中文字幕有乳无码| 成人午夜亚洲影视在线观看| 天天综合网站| 国产哺乳奶水91在线播放| 亚洲无码高清一区二区| 亚洲黄色高清| 国产精品部在线观看| 国产亚洲成AⅤ人片在线观看| 18禁不卡免费网站| 国产在线精品99一区不卡| 一级毛片在线播放| 国产男女免费完整版视频| 国产免费人成视频网| 亚洲免费三区| 国产日本一区二区三区| 久久午夜夜伦鲁鲁片无码免费| 国产成人三级| 国模沟沟一区二区三区| 国产中文一区a级毛片视频| 无码 在线 在线| 国模沟沟一区二区三区| 国产成人成人一区二区| 71pao成人国产永久免费视频| 国产精品亚欧美一区二区三区| 人与鲁专区|