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

浸沒式撞擊流反應器流場渦特性的數值研究

2022-09-13 07:58:20張建偉高偉峰董鑫馮穎
化工學報 2022年8期

張建偉,高偉峰,董鑫,馮穎

(沈陽化工大學機械與動力工程學院,遼寧 沈陽 110142)

引 言

撞擊流技術因其兩股流體高速對置撞擊,形成劇烈湍動的撞擊區,為強化傳質傳熱提供了良好條件[1],廣泛地應用在石油化工[2]、煤化工[3]等過程工業。撞擊流反應器內流體劇烈湍動,流場內存在大量的渦。渦作為一種典型的湍流結構[4-8],眾多學者已對撞擊流反應器的渦特性進行研究。渦通過有旋運動卷吸周圍流體,使其與渦核內部流體接觸,以此達到混合目的。此外渦的演變有效抑制了流場內“死區”的形成,改善了流體的混合效率。Schwertfirm 等[9]研究受限撞擊流反應器發現,渦是影響管內流動的主要因素。杜柯江等[10]對小型受限撞擊流反應器研究發現,在撞擊面上周期性地生成旋渦,撞擊區域的偏斜振蕩對渦存在影響。Gao 等[11]采用PIV 測試系統,研究了受限撞擊流反應器內撞擊駐點穩定性,發現Re和射流速度比與撞擊駐點的穩定性有關。撞擊駐點的穩定性直接影響反應器內的渦脫落與運動形式。張經緯[12]認為撞擊面兩側的渦生成和脫落是由撞擊駐點振蕩引起的。B?hm等[13]運用PIV 技術發現了對置射流燃燒器剪切層渦的脫落頻率為1~2 kHz。Sultan 等[14]確定了T 型撞擊流反應器內渦向吞噬流轉變的Reynolds數和幾何參數。Yi等[15]發現Y 型撞擊流反應器彎曲的噴嘴引起的反向渦有助于提高混合效率。溫謙等[16]通過TRPIV 技術發現淹沒射流中有序的相干渦結構在上游逐漸發展。Judd等[17]通過高分辨率紅外傳感器(IR)發現在低Re下湍流射流的撞擊區域的中心發出卷曲渦。在撞擊流反應器內由于流體在空間上存在速度梯度,流場因此產生了渦結構,其中大尺度的渦改變了流型。流場內流型的變化直接影響著流體的混合效果,通過研究流場渦特性以及渦的運動機理,為提高流體混合效率奠定基礎。

目前,研究者對撞擊流反應器內流場渦特性開展了相關理論研究。許鑫磊等[18]使用大渦模擬(LES)方法研究在層流條件下T 型撞擊流反應器內流體流動形式,揭示了非穩態吞噬流機理。Wu等[19]引用有限時間Lyapunov 指數(FTLE)場識別了湍流中的渦結構。Li 等[20]通過LES 和k-ε模型瞬態模擬對并行式湍流撞擊流進行研究,發現匯合點附近的渦旋相互作用引起重建的流場。Hrisheekesh等[21]采用格子Boltzmann方法確定了渦脫落的基本特征,得到對稱渦脫落強度隨Re的增加而增大的結論。郭棟等[22]采用Realizablek-ε模型研究氣液式撞擊流反應器粒子停留時間,并得出粒子停留時間的最佳工況。Zhang 等[23]通過LES 發現在波狀圓柱的尾流區存在多尺度和肋狀渦現象,并用渦通量揭示圓柱的邊界渦旋的產生機理。Zhao等[24]通過LES 研究了在有渦發生器的擾動下傾斜射流中的渦演過程。目前,前人多數改變撞擊流反應器的幾何參數和工藝條件來研究流場內部旋渦形成的規律。但是,關于渦的形成規律、演化形式和作用機理的研究較少,影響撞擊流反應器渦量場和渦特性的因素也有待深入研究。

本文采用LES方法對浸沒式撞擊流反應器渦特性進行數值模擬研究。討論不同進口速度(v0)、噴嘴間距(L/D)對流場渦量及渦能量的影響,結合Q判據分析渦的演化過程。根據流線圖分析渦系的作用范圍,判定流體的流動形式。然后考察不同工況下流場內的平均渦量和渦能量分布規律,得到流場的最優工況。本文的研究成果可為撞擊流反應器的渦的特性提供理論依據,豐富撞擊流反應器流場的內涵,為研究渦與混合關系提供理論基礎和指導思路。

1 計算模型和數值方法

1.1 幾何模型

反應器模型如圖1 所示,通常撞擊流反應器由噴嘴、溢流口、排水口、殼體構成。撞擊流反應器由厚度為5 mm的有機玻璃制成。本文研究中,撞擊流反應器的殼體內直徑D1=130 mm,高度H=500 mm,兩側噴嘴居中對稱放置,噴嘴直徑為D=10 mm,噴嘴間距為L,溢流口直徑為D2=30 mm,其中心距反應器頂部為H2=30 mm,排水口直徑D3=30 mm。

圖1 撞擊流反應器幾何模型Fig.1 Geometric model of impinging stream reactor

1.2 控制方程

Huang 等[25]討論了高溫浮力射流的流場演化和渦結構特征,通過大渦模擬方法發現初始溫度對渦環的影響。Zhang 等[26]通過大渦模擬研究了圓柱繞流的剪切層不穩定性,討論了下游尾流剪切層中Kelvin Helmholtz 不穩定性的特征。眾多學者運用大渦模擬方法研究了渦的演化規律和成因,表明大渦模擬方法能很好地模擬流場渦的特性。本文采用大渦模擬方法研究了撞擊流反應器內的流場旋渦特性。

質量守恒方程:

動量方程:

式中,τij為亞格子尺度應力,對于亞格子應力常采用渦黏性假設確定,其定義為:

式中,CS[27]為模型常數,CS=0.21。

1.3 介質和邊界條件

對模擬方法進行設置,壓力-速度耦合采用SIMPLE 方法,其特點是較易收斂。壓力采用二階格式,動量方程采用有界中心差分格式,時間步長設置為0.01 s,Courant 數設置為1.0。介質為常溫常壓的水,密度為0.998×103kg/m3,黏度為2.98×10-3Pa·s。采用壓力出口邊界條件,并假設為充分發展,壓力為101.325 kPa。進口速度(v0)分別取0.5,0.8,1.0,1.3,1.77,2.0 m/s;無量綱噴嘴間距(L/D)分別取1, 3, 5。采用標準壁面函數進行處理,壁面設置無滑移邊界條件。因實驗時撞擊流反應器外壁為有機玻璃,壁面粗糙度系數(roughness constant)設置為默認值(0.23)。

2 模型的驗證

2.1 網格獨立性檢驗

利用Mesh 模塊,對模型進行網格劃分,噴嘴附近區域進行加密處理,以便更好地捕捉到噴嘴附近的渦結構。按上述方法繪制了305603、498602、716182 和1134135 個網格,采用相同的計算方法比較了4 組網格的徑向速度,如圖2 所示。可以看出當網格數量為716182個時,能準確地反映該位置的徑向速度分布規律,為了保證計算的精確性和節約成本,在網格數為716182個條件下模擬計算。

圖2 網格無關性檢驗Fig.2 Grid independence test

2.2 模擬可靠性檢驗

為驗證模擬的準確性,采用TR-PIV 測試系統對撞擊流反應器的速度場進行測定,實驗裝置如圖3 所示。該TR-PIV 測試系統由丹麥DANTEC 公司生產,主要包括連續激光器、同步器、CCD 攝像機及圖像處理系統四部分。雙腔Nd:YAG 連續激光器(RayPower2000)發射輸出功率為4 W、波長為532 nm 的片狀激光束。CCD 攝像機是FlowsenseEo 攝像機。采集頻率為15 Hz,連續采集300幅瞬態速度場圖像。密度為1.05~1.15 g/cm3和直徑為10~15 μm的空心玻璃球(HGS)用作示蹤顆粒,由于示蹤粒子的尺寸和密度較小,對流場的擾動可忽略不計。拍攝時攝像機采用雙幀模式。

圖3 撞擊流反應器實驗示意圖Fig.3 Schematic diagram of impinging stream reactor experiment

在測量反應器內速度場時,先將示蹤粒子與實驗流體混合,配制成示蹤粒子液體。示蹤粒子液體由蠕動泵和經過離心泵的實驗流體一同通入反應器。待流場穩定后(粒子充分布滿整個反應器),進行標定實驗,保證測試區域的流場質點均可以被捕捉并確保了速度測試的準確性。

對撞擊流反應器進行TR-PIV 實驗測量,實驗進口速度控制在1.77 m/s,得出撞擊流反應器的徑向速度圖,如圖4 所示。圖中數值計算出的徑向速度與實驗測得徑向速度較好吻合,說明采用LES 模型能夠準確計算出撞擊流反應器內的速度場,通過模擬計算出的數據真實可靠。

圖4 實驗與模擬徑向速度對比Fig.4 Comparison of experimental and simulated radial velocity

3 結果與討論

3.1 撞擊流反應器流場速度分析

為討論撞擊流反應器內流體的運動形式,圖5給出不同軸向位置(-5≤Y/D≤5)流場的速度流線圖。撞擊流反應器內流線的變化關于反應器的圓心呈中心對稱。其原因為流體從噴嘴流出后發生撞擊形成撞擊面,撞擊面再次與反應器壁碰撞形成徑向回流,在徑向回流的作用下,反應器內形成了以反應器中心為對稱中心的渦。此外,在XOZ面內反應器渦系共分為4 部分,這表明在該平面軸向射流與徑向射流對流場內的流體具有分割作用,使得渦在各自渦系內演變,渦系之間演變不相互影響。在Y/D=±1 位置,流場內的“點源”“點匯”最多,這體現出流體發生撞擊之后,經過一定時間的運動與能量轉換,渦系發生改變,湍動程度變強。隨著截取面遠離噴嘴,渦的變化趨勢呈小渦向大渦轉化。根據流線圖可知,隨著流體與噴嘴遠離,渦的影響區域逐漸變大。Y/D=±5處的渦旋影響范圍最大,渦旋結構為近似對稱的Dean 渦。渦影響區變大,有利于混合。

圖5 不同軸向位置流線圖Fig.5 Streamline of different axial positions

3.2 撞擊流反應器內流場渦量分析

在撞擊流反應器內,由于兩噴嘴射出的流體高速撞擊,形成了劇烈湍動區,流體發生剪切、扭曲等有旋運動。渦量[28]是衡量流體有旋運動程度的物理量,其計算公式如下:

式中,ω為渦量;?為哈密頓算子;v為流體速度。渦量大的位置表示流體有旋運動程度高,湍動劇烈。

為探究撞擊流渦量場,將流場內的渦量歸一化處理成流場平面平均渦量(ωˉ)。分析了撞擊流反應器不同噴嘴間距(L/D)、不同進口速度(v0)對流場平均渦量(ωˉ)的影響,如圖6。隨著入口速度的增加,ωˉ呈先增大后減小的趨勢。v0>1.77 m/s 后,流場內的進口動量過高,不利于流體平動動能向轉動動能轉化,有旋運動程度降低,渦量因此減小。L/D=1,兩噴嘴間距過小,流體撞擊后得不到充分發展,渦量降低。噴嘴間距增大到L/D=5,此時噴嘴間距過大,導致流體在流場內耗散作用增強,能量耗散過多,導致發生撞擊時能量不足,有旋運動程度降低。

圖6 流場平均渦量對比Fig.6 Comparison of mean vorticity in flow field

為探究撞擊流反應器流場內的渦量分布規律,圖7 給出L/D=3、不同進口速度下,流場的徑向渦量分布規律。在反應器的徑向位置,其渦量分布呈“M”形。在撞擊駐點流體的速度最小,有旋運動程度低,這說明流體在撞擊中軸線位置發生能量轉換,因此渦量較低。隨著流體徑向運動,渦量逐漸增大,表明從噴嘴射出的兩股流體合流后動量增大,并在噴嘴出口附近達到峰值后,由于耗散作用渦量逐漸減小。增大進口速度后,導致進口動能增大,流場內的流體湍動程度增加,流體發生撞擊后在徑向位置的渦量逐漸增大。

圖7 徑向渦量分布Fig.7 Radial vorticity magnitude distribution

為討論噴嘴間距對流場渦量分布規律的影響,圖8 給出v0=1.77 m/s,不同L/D下流場Y/D=1 位置渦量分布規律。根據圖8分析可知,L/D=3時軸向渦量分布最為理想,并在撞擊區域位置存在峰值。同時流體在反應器邊緣的作用下,發生強烈的扭轉因此渦量得以升高。當L/D=1 時,流體發生撞擊后將大部分動能轉化為靜壓能,速度梯度隨之減小,因此渦量較低。噴嘴間距過大(L/D=5),流體在進入反應器后發生能量耗散,導致流體撞擊過程動量不完全轉化,速度降低進而速度梯度減小,因此渦量較低。

圖8 Y/D=1渦量分布Fig.8 Vorticity magnitude distribution at Y/D=1

為直觀描述流場內渦量作用區域,探究v0=1.77 m/s,L/D=3 工況下,撞擊流反應器的渦量三維等值面,如圖9 所示。由圖9 可知,流體經過撞擊后遠離撞擊駐點,渦量逐漸減小,駐點周圍流體由軸向流動變為徑向流動。隨著渦量等值面選取值降低,發現在徑向回流區域的流體存在有旋運動。這表明,徑向回流也是一種有旋運動的流體,但渦量值較小。根據圖9(a)~(c)分析ω>100 s-1等值面接近于圓形,這表明,流體撞擊之后以有旋運動的形式向四周擴散。當ω=75 s-1等值面邊緣出現了波動和彎曲,此時流體的有旋運動減弱,容易受撞擊駐點振蕩和撞擊面偏轉的影響。

圖9 撞擊流反應器渦量等值面Fig.9 Vorticity magnitude of iso-surface of impinging stream reactor

為探究撞擊流反應器的最優工況下,渦量在不同位置的分布規律,分析了XOY面渦量分布曲線,如圖10 所示,由于撞擊流反應器結構關于噴嘴對稱,研究Y軸正方向(0~200 mm)。在撞擊面兩側渦量存在峰值,且在Y=10 mm位置渦量最大,沿Y軸向上渦量峰值逐漸減小,渦體攜帶的能量隨著Y坐標值的增加而持續衰落,說明流體遠離駐點有旋運動程度持續減弱。

圖10 不同位置渦量分布曲線Fig.10 Vorticity magnitude distribution curve at different positions

3.3 流場內渦能量分析

在撞擊流反應器內由于流體流動形式復雜,為準確描述渦存在位置,采用Hunt 等[29]提出的Q判據來確定流場內渦的位置。同時,高助威等[30]從能量的角度理解Q可以表示單位時間、單位空間時均流渦具有的能量。其公式為:

式中,‖· ‖為張量的二范數;Ω為反對稱渦張量;S為對稱應變率張量。區域內流體的反對稱渦張量大于對稱應變率張量,即Q>0,表明該位置存在渦;反之不存在。此外,Q值越高證明該區域的流體攜帶的能量越大,有利于混合。

為探究不同進口速度(v0)和噴嘴間距(L/D)對流場渦能量的影響,討論不同工況條件下流場的平均渦能量(),如圖11 所示。根據流場平均渦能量()分析,撞擊流反應器最優工況條件為v0=1.77 m/s、L/D=3,流場內有足夠的進口動量和撞擊空間,流體經過撞擊后,渦得到了充分的發展。噴嘴間距過小,流體的動能撞擊后多轉化為撞擊駐點的靜壓能,致使渦能量不高。過大的噴嘴間距使流體在撞擊之前動能有較大的耗散導致渦能量較低。

圖11 不同工況下流場平均渦能量Fig.11 -Q of flow field under different working conditions

3.4 撞擊流反應器流場渦結構

為明確流場內的渦核區的存在范圍與形狀,討論了v0=1.77 m/s,L/D=3 工況下,流場不同Q值的等值線渦核分布,結果如圖12所示。可以看出撞擊流反應器內渦主要分布在撞擊區域和出口區域。結合圖12(a)分析,此時Q取值較小(0.01 s-2),因此其分布范圍也更廣,主要存在于反應器的撞擊區域和反應器上半部分,這是由于上部存在出口,流場內的流體可以得到充分的發展,極大地促進渦的形成。同時,在反應器的出口位置存在渦環,表明流體從反應器的出口流出時,部分流體的流向與出口位置的流體速度垂直,進而在出口處形成渦環。隨著Q值的增加渦核的分布范圍逐漸變小,如圖12(b)所示,Q=0.04 s-2時整個撞擊區域的渦核結構凸顯,撞擊區域的渦包括馬蹄渦與肋狀渦。馬蹄渦主要存在于撞擊駐點較遠的撞擊面邊緣,而靠近撞擊區的中心渦多以肋狀渦為主,這表明流場內馬蹄渦的存在區域內流場的湍動程度相對于肋狀渦較低。由于撞擊流反應器撞擊區域的渦多以肋狀渦的形式存在,其湍動程度高、演變時間短,有利于物料的傳遞;此外,馬蹄渦存在于反應器壁附近,由于馬蹄渦的尺度較大,可將沉積的物料卷吸湍動區域,強化傳質傳熱。

圖12 撞擊區域的渦核分布Fig.12 Vortex core distribution in impact region

3.5 撞擊流反應器流場渦演化過程

研究了進口速度v0=1.0 m/s下,流場XOZ面的Q值分布規律。根據渦在反應器內的位置和運動規律,將其演化過程總結為5 個階段,分別為:撞擊階段、徑向運動階段、徑向分離階段、周向運動階段、周向破裂階段,其變化過程如圖13所示。

如圖13(a)、(b)所示,發生撞擊時,流體首先呈渦的形式發生撞擊,并且渦核呈圓狀。進口速度為1.0 m/s,噴嘴間距為30 mm,在0.03 s 后兩股流體發生撞擊,如圖13(b)所示。兩股流體撞擊后,撞擊面壓力變大,渦體受壓力作用發生變化,形狀由圓形變成橢圓形,如圖13(c)所示,此時渦核完成撞擊。

圖13 XOZ面渦演化過程示意圖Fig.13 Schematic diagram of the vortex evolution of XOZ surface process

圖13(c)~(k)所示為徑向運動階段,在撞擊面剪切力的作用下使渦發生分離,分離結束如圖13(d)所示。分離之后渦對繼續沿撞擊面位置向反應器邊緣運動,如圖13(e)所示,但隨著渦體在徑向運動,其渦體之間發生演化,該現象稱為徑向融合。

如圖13(e)~(g)所示,同向運動的一對渦發生變化,不再以同等能量進行運動,而是一個渦逐漸變強,另一個逐漸變弱,這表明流體發生撞擊之后,撞擊駐點偏移已經存在因此會存在該現象。當出現徑向融合時,如圖13(e)所示,撞擊中心的渦能量逐漸降低,說明撞擊駐點振蕩會導致其附近的渦對能量降低,徑向融合結束后,如圖13(g)所示,駐點附近渦能量降到最低,此時徑向運動的渦能量最高。隨后在撞擊面邊緣擺動的作用下融合的渦體發生分離,繼續向壁面運動。

當渦結束徑向運動之后開始進行徑向分離階段[圖13(l)、(m)],兩對渦到達反應器邊緣,此刻渦對并未立即分裂,而是能量有所增大,增加至如圖13(l)所示,才發生分裂過程,這表明渦不是在撞擊壁面后立即發生演變,而是當渦體的能量增加到一定數值后才發生分裂,該數值稱為渦脫落強度的無量綱環流強度[21](Γω,max),如圖13(m) 所示,渦分離基本結束。

渦結束徑向分離之后開始進行周向運動,該階段的典型特征是渦對在徑向回流的作用下,向噴嘴與反應器夾角間運動,能量進一步提升,如圖13(m)、(n)所示,此時渦的運動達到極限。

周向破裂階段,渦體受反應器的空間限制發生渦核破裂,如圖13(o)所示。之后由于徑向回流的持續作用,渦系向撞擊駐點擴散。在整個演化過程中渦的徑向運動和周向運動的時間最長,兩者共占整個演化過程的34.1%。

在撞擊流反應器內徑向運動、徑向分離和周向運動為渦典型的演化過程,因此圖14~圖16 對上述過程持續時間進行分析。在徑向運動和徑向分離過程中,增大進口速度可有效縮短渦演化時間,表明流場進口速度增大后,徑向射流所攜帶的能量高,使渦徑向演化時間縮短。L/D=3 渦的徑向分離時間最短,縮短徑向分離時間可使渦劇烈湍動,提高混合效率。周向運動過程中在L/D=3,v0=1.77 m/s的條件下,周向運動時間較長,渦在周向運動過程中,通過渦的卷積和運輸作用使沉積在反應器邊緣的物料再次到達湍動區,使物料重新參與到反應過程中。因此根據各階段的演化時間分析,在L/D=3,v0=1.77 m/s 條件下渦的徑向分離時間最短,周向運動時間較長,可確定其為流體混合的最佳工況。

圖14 XOY面渦徑向運動時間Fig.14 Radial movement time of vortex under XOY surface

圖15 XOY面渦徑向分離時間Fig.15 Radial separation time of vortex under XOY surface

圖16 XOY面渦周向運動時間Fig.16 Peripheral motion time of vortex under XOY surface

3.6 撞擊流反應器徑向射流渦周期性分析

流體在渦的作用下劇烈擾動,使徑向射流與壁面撞擊后在其兩側的渦交替出現。采用Q判據分析了不同噴嘴間距(L/D)下撞擊流反應器XOZ面徑向射流兩側渦的演化規律,如圖17 所示。如圖17(a)紅圈標出,在撞擊面的兩側存在一對渦結構,由于兩噴嘴對稱撞擊,撞擊后形成的徑向射流方向與反應器壁面剛好垂直,流體在撞擊壁面之后改變原有運動方向,從徑向射流轉為周向流動,產生一對渦結構。經過0.04 s 之后在撞擊面的右側渦消失,僅在左側存在一個渦,這表明由于撞擊駐點的振蕩導致撞擊面發生偏移,撞擊面流體多數向左側流動,左側存在渦結構,同理可解釋圖17(d)。如圖17(d)~(e)所示,再經過0.06 s 之后撞擊面兩側再次出現一對渦,此刻撞擊面的流體速度與反應器壁面再次垂直。整個過程渦的演化周期(tc)為0.18 s。

圖17 渦周期性演化示意圖Fig.17 Schematic diagram of vortex periodic evolution

為清晰闡述L/D對tc的影響,表1給出v0=1.77 m/s 條件下的L/D與tc的關系。分析表1 可知,演化周期(tc)隨著噴嘴間距(L/D)的增加呈先增大后減小趨勢,在L/D=3 時,tc存在最大值。由于渦的演化過程伴隨著渦能量的變化,渦的交變周期越長,渦在反應器邊緣存在的時間越長,越有利于提高反應器內流體的混合效率。此外,周期性演變的渦,有效地改善流場“死區”問題。

表1 撞擊流反應器不同L/D和tc關系Table 1 Comparison of L/D and tc of impinging stream reactor

4 結 論

本文通過大渦模擬(LES)方法研究了撞擊流反應器內流場渦運動規律,利用渦量(ω)和Q值定量分析了渦的有旋運動和能量,得到結論如下。

(1)分析了撞擊流反應器內撞擊區流體流動的形式和渦的作用范圍。發現流場存在軸向射流、徑向射流和徑向回流3 種流型,三者的作用使反應器內產生了4 個對等的渦系。流體靠近噴嘴時,渦的作用范圍小,能量高,遠離噴嘴后渦的范圍逐漸增大,能量降低,在Y/D=±5 位置形成了范圍最大的類Dean渦。

(2)平均渦量和平均渦能量在L/D=3,v0=1.77 m/s 工況下達到最大,流體湍動程度最劇烈,有利于混合。利用三維渦量等值線,分析了撞擊流反應器內渦量的空間分布,高渦量的流體主要存在于撞擊區靠近噴嘴位置。

(3)采用Q判據研究了流場內渦的演化過程,根據渦所在位置和運動路徑,將演化過程分為撞擊階段、徑向運動階段、徑向分離階段、周向運動階段、周向破裂階段5 個階段,其中徑向運動階段和周向運動階段占渦演化過程的34.1%。

(4)確定了XOZ面徑向射流兩側渦的演化周期,周期為0.15~0.20 s。反應器的渦結構主要為馬蹄渦和肋狀渦,肋狀渦主要存在于噴嘴附近,大尺度馬蹄渦主要存在于撞擊區域附近,更有利于混合。

主站蜘蛛池模板: 国产精品久久久久婷婷五月| 国产主播在线一区| 日本日韩欧美| 成人在线天堂| 欧美综合区自拍亚洲综合绿色 | 色欲国产一区二区日韩欧美| 美女扒开下面流白浆在线试听| 欧美日韩中文字幕二区三区| 91丝袜乱伦| 国产日韩欧美成人| 久久国语对白| 欧美成人午夜视频免看| 四虎AV麻豆| 日本不卡免费高清视频| 十八禁美女裸体网站| 欧美不卡视频在线观看| 日韩一区二区三免费高清| 亚洲va在线∨a天堂va欧美va| 亚洲精品国产综合99| 无码精品国产VA在线观看DVD| 亚洲国产日韩在线成人蜜芽| 456亚洲人成高清在线| 日韩av手机在线| 久草国产在线观看| 亚洲中文字幕在线一区播放| 国产成人艳妇AA视频在线| 亚洲成在人线av品善网好看| 在线看AV天堂| 国产99视频精品免费视频7| 久久亚洲AⅤ无码精品午夜麻豆| 国产性生交xxxxx免费| 伊人色在线视频| 国产成人av一区二区三区| 九九九久久国产精品| 国产精品偷伦在线观看| 夜夜操天天摸| 久久精品这里只有精99品| 日韩视频福利| 国产人成乱码视频免费观看| 蜜臀av性久久久久蜜臀aⅴ麻豆 | 91欧美亚洲国产五月天| 国产精品永久久久久| 露脸一二三区国语对白| 黄色网在线免费观看| 国产精品亚洲一区二区三区在线观看| 亚洲av无码牛牛影视在线二区| 国产在线日本| 在线免费观看AV| 日韩在线播放欧美字幕| www亚洲精品| 久久人搡人人玩人妻精品| 无码AV高清毛片中国一级毛片| 亚洲AV无码久久精品色欲| 国产精品不卡片视频免费观看| 伊人久久婷婷五月综合97色| 最新国产成人剧情在线播放| 好吊色妇女免费视频免费| 亚洲日韩日本中文在线| 国产乱子伦手机在线| 正在播放久久| 粗大猛烈进出高潮视频无码| 欧美亚洲日韩不卡在线在线观看| 91色国产在线| 青青草国产精品久久久久| 国产精品美人久久久久久AV| 色综合中文| 久久午夜夜伦鲁鲁片无码免费| 婷婷六月激情综合一区| 911亚洲精品| 99成人在线观看| 久久亚洲国产视频| 国产在线观看第二页| 综合色亚洲| 国产精品人成在线播放| 免费国产不卡午夜福在线观看| 99久久人妻精品免费二区| 久久中文字幕不卡一二区| 一本无码在线观看| 欧美啪啪精品| 亚洲性网站| 日韩中文无码av超清| 日韩经典精品无码一区二区|