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

玻璃微珠液滴碰撞分離過程數(shù)值研究

2024-01-05 00:15:20李慧玲胡曉磊余子寒謝能剛

李慧玲, 胡曉磊, 余子寒, 謝能剛

(安徽工業(yè)大學(xué) 機械工程學(xué)院, 安徽 馬鞍山 243032)

0 引 言

玻璃微珠是近年來發(fā)展起來的一種用途廣泛、性能特殊的新型無機非金屬材料.玻璃微珠具有強度高、流動性好、光學(xué)性能優(yōu)異等顯著特點,因此被廣泛應(yīng)用于軍事、航空、建材及其他高科技領(lǐng)域作光反射材料、凈化拋光材料、填充材料、載體材料等.隨著科學(xué)技術(shù)的不斷發(fā)展和用戶對產(chǎn)品性能要求的不斷提高,玻璃微珠的直徑越來越小,范圍可達(dá)到數(shù)微米到數(shù)百微米[1].玻璃微珠的生產(chǎn)方法分為物理法和化學(xué)法,采用物理法制備玻璃微珠的技術(shù)發(fā)展較為成熟.生產(chǎn)中常用的固相玻璃粉末法就是一種物理法,該方法常將玻璃粉料和燃料噴入球化爐內(nèi),使玻璃粉末顆粒充分受熱熔化,獲得熔融的玻璃微珠液滴,在此過程中玻璃微珠液滴之間大多發(fā)生碰撞分離現(xiàn)象,但由于爐內(nèi)溫度過高,該微觀運動過程難以捕捉[2].

為準(zhǔn)確預(yù)測玻璃微珠液滴碰撞過程,需要對玻璃微珠液滴碰撞過程進(jìn)行數(shù)值模擬.液滴碰撞過程是一個涉及多尺度、非線性問題的復(fù)雜流動過程,對于一些時間和空間的微觀信息,實驗難以捕捉.與實驗研究相比,數(shù)值模擬在流場及參數(shù)設(shè)置方面具有更大柔性.在此方面,夏盛勇等[3]采用流體體積(VOF)方法對直徑為10~200 μm的三氧化二鋁液滴碰撞過程進(jìn)行了數(shù)值模擬,獲得了反彈、大變形后聚合和自反分離等結(jié)果類型,研究了不同Weber數(shù)下的結(jié)果,并得到了它們之間的臨界Weber數(shù).Fostiropoulos等[4]研究了水乳化重油燃料爆炸液滴的碰撞行為,用3個VOF傳輸方程分離共存水液體和水蒸氣流體界面,分析模擬結(jié)果,得到液滴破碎只發(fā)生在燃料-空氣界面的一部分,類似于噴氣的特征.尹強等[5]將VOF方法用于研究循環(huán)冷卻系統(tǒng)中海水液滴碰撞過程,得到了聚合和自反分離兩種結(jié)果類型及二者的臨界Weber數(shù),分析了流場結(jié)構(gòu)、液滴直徑和海水濃度對碰撞結(jié)果的影響.

Amani等[6]使用守恒的水平集(level set)方法模擬了所有體系中正面和偏離中心的二元液滴碰撞,介紹了一種新的薄層穩(wěn)定方法來數(shù)值解析薄膜,并對每個案例進(jìn)行了深入的能量分析,涵蓋了廣泛的碰撞范圍,為碰撞過程提供了更深入的見解.Ashna等[7]采用LBM方法研究了液滴的蒸發(fā)過程,引入了有效表面張力測量蒸發(fā)速率對碰撞結(jié)果的影響,得到了其隨Stefan數(shù)或Weber數(shù)增加而減少,隨時間的增加而增加的規(guī)律.Liu等[8]通過將實驗可視化與數(shù)值模擬相結(jié)合,研究了雙乳化液滴的水動力剪驅(qū)動二元碰撞,發(fā)現(xiàn)了過渡和反轉(zhuǎn)這兩種典型碰撞運動,并得到了它們之間的決定影響因素.Zhang等[9]用分子動力學(xué)的方法研究了尺寸對二元碰撞納米液滴動力學(xué)的影響,對不同直徑液滴的拉伸因子、能量耗散和碰撞結(jié)果進(jìn)行分析,得到了納米尺度上原子相互作用引起的動能耗散對液滴分離的發(fā)生或其他方面有顯著影響.

Li等[10]采用CLSVOF方法的數(shù)值模型,研究了空心液滴對干燥平面沖擊的動力學(xué)和傳熱問題,比較了高空化氣泡壓力的空心液滴與連續(xù)致密液滴的動力學(xué)和傳熱特性,得到了液滴在這兩種條件下的碰撞規(guī)律.Qian等[11]利用CLSVOF方法和自適應(yīng)網(wǎng)格細(xì)化,對Newton液體和聚合物液體的碰撞行為進(jìn)行了數(shù)值研究.對液滴的變形過程、內(nèi)部流場和能量演化進(jìn)行了詳細(xì)分析,得到了二者在液滴碰撞過程中的不同之處,主要分析了聚合物液滴碰撞的一些變化特征.由上述文獻(xiàn)可知,國內(nèi)外對液滴碰撞機理的研究已然十分成熟.然而,截至目前,國內(nèi)外尚未見公開發(fā)表針對玻璃微珠液滴碰撞物理規(guī)律及模型的研究.

鑒于對熔融玻璃微珠液滴碰撞過程實驗觀測的困難性,本文采用CLSVOF方法[12]開展兩個相同尺寸玻璃微珠液滴對心碰撞的數(shù)值模擬.首先闡述了液滴碰撞的基本模型以及相關(guān)參數(shù)計算公式,隨后建立了玻璃微珠碰撞過程數(shù)值模型和計算方法,采用建立的數(shù)值方法對Qian等[13]的實驗開展數(shù)值驗證,最后研究了玻璃微珠液滴對心碰撞分離的物理規(guī)律.

1 液滴碰撞基本模型

圖1給出了兩個相同尺寸液滴對心碰撞模型示意圖,其中D0是液滴直徑,v是液滴速度,L為兩液滴的初始間距.計算區(qū)域為5D0×2.5D0.由于模型具有軸對稱特性,本文采用二維軸對稱計算條件,計算域底部為軸對稱邊界條件,其他三邊均為壓力出口,計算區(qū)域如圖1所示.

圖1 物理模型Fig. 1 The physical model

2 數(shù)值計算方法

液滴碰撞是兩相流問題,目前,兩相界面的捕獲方法主要有:VOF法[14]、水平集[15]以及 CLSVOF方法.VOF法難于計算相界面曲率,會導(dǎo)致解的振蕩以及界面的陡峭變化被抹平.水平集法在計算時不可避免地引入數(shù)值耗散和舍入誤差,造成封閉界面內(nèi)流體質(zhì)量損失,從而導(dǎo)致質(zhì)量不守恒.VOF法和水平集法具有各自的優(yōu)勢和劣勢,但是這兩種方法的優(yōu)、劣勢正好可以互補.因此,Sussman等[16]提出使用VOF法和水平集法相耦合的方法,即CLSVOF方法進(jìn)行相界面追蹤,該方法能準(zhǔn)確捕獲相界面的同時,保證質(zhì)量守恒.

對于流體力學(xué)問題,涉及到的基本控制方程如下所示:

?·v=0,

(1)

(2)

式中v為速度矢量,m/s2;p為流體壓力,Pa;F為表面張力源項;g為重力加速度,m/s2;μ為流體動力黏度.CLSVOF方法在捕獲兩相界面時涉及到的水平集函數(shù)為

(3)

式中Ω1,Ω2分別是液相和氣相,Γ為兩相界面,L為區(qū)域內(nèi)點到界面的距離.

通過引入Heaviside函數(shù)來使界面的密度和黏度光滑過渡:

(4)

式中a為1.5倍的最小網(wǎng)格尺寸.

式(2)的表面張力源項F為

F=σκ(φ)?H(φ)n,

(5)

式中,σ是液相的表面張力系數(shù),n是兩相界面的法線,κ(φ)是兩相界面的曲率.法線n和曲率κ(φ)的表達(dá)式如下:

(6)

(7)

光滑后的密度和黏度為

ρ(φ)=ρg+(ρl-ρg)H(φ),

(8)

μ(φ)=μg+(μl-μg)H(φ),

(9)

式中ρ(φ)為各相的平均密度,kg/m3;ρg為氣相的密度,kg/m3;ρl為液相的密度,kg/m3;μ(φ)為各相的平均動力黏度,N·s/m2;μg為氣相的動力黏度,N·s/m2;μl為液相的動力黏度,N·s/m2.

對于碰撞結(jié)果描述涉及到的Weber數(shù)和Reynolds數(shù)公式定義如下:

(10)

(11)

式中,ρ為液滴密度,kg/m3;d為液滴直徑,m;u為液滴速度,m/s;σ為表面張力系數(shù),N/m;μ為液滴動力黏度,N·s/m2.

采用PISO方法耦合求解速度場和壓力場,時間項采用Euler方法進(jìn)行離散,壓力項采用PRESTO方法進(jìn)行求解,動量方程采用二階迎風(fēng)格式進(jìn)行求解,水平集函數(shù)采用一階迎風(fēng)格式離散求解.

3 數(shù)值算法的驗證

為驗證數(shù)值方法是否適合計算液滴碰撞及其可靠性,本文針對文獻(xiàn)[13]正十四烷液滴在1 atm(1 atm=101 325 Pa)的常溫氮氣環(huán)境下的實驗進(jìn)行數(shù)值模擬驗證計算,計算流體物性參數(shù)見表1.

針對文獻(xiàn)[13]中圖4(h)近對心碰撞工況開展了數(shù)值計算,計算結(jié)果如圖2所示,其中實驗結(jié)果為三維拍攝圖.對比計算及實驗結(jié)果可以看出,計算最終結(jié)果是自反分離,取得了與實驗基本一致的液滴變形過程,只是在時刻上略微有所不同.

表1 正十四烷液滴的物理參數(shù)

(a) 實驗結(jié)果[13] (b) 計算結(jié)果 (a) Experimental results[13] (b) Calculation results

4 計算結(jié)果及分析

4.1 玻璃微珠液滴及氣相物性參數(shù)

玻璃微珠液滴的尺寸一般介于幾微米到幾百微米之間,本文選擇兩個尺寸相同的玻璃微珠液滴計算,直徑均為20 μm.模擬中所采用的物性參數(shù)是參考彭壽等[2]研究玻璃微珠球化時所用的參數(shù),如表2所示.

表2 玻璃微珠液滴的物理參數(shù)

4.2 玻璃微珠液滴對心碰撞數(shù)值研究

本文僅研究兩個相同尺寸玻璃微珠液滴對心碰撞分離的規(guī)律特性.在Weber數(shù)為149的工況下,液滴碰撞分離形成了多個衛(wèi)星液滴,但最終只有一個較大的衛(wèi)星液滴保持穩(wěn)定形態(tài),因此選取它來進(jìn)行玻璃微珠液滴碰撞分離的規(guī)律特性研究.通過對兩玻璃微珠液滴碰撞過程中的形態(tài)變化、 能量變化、 中心處壓力變化和關(guān)鍵時刻的速度矢量圖以及壓力云圖的分析, 研究玻璃微珠液滴碰撞后發(fā)生分離的原因.模擬結(jié)果如圖3所示.

圖3為兩相同尺寸玻璃微珠液滴碰撞后自反分離的計算結(jié)果.兩液滴突破氣膜后發(fā)生聚合,0.5 μs聚合后在表面張力的作用下發(fā)生徑向伸展變形,形成片狀射流;在2.0 μs時液滴徑向拉伸到極限,之后液滴徑向收縮,并同時向軸向伸展;在6.0 μs時,液滴軸向拉伸到極限,不再延展,但徑向仍在不斷收縮,形成液橋;在11.5 μs時液滴發(fā)生分離,液橋分裂形成3個衛(wèi)星液滴,在表面張力的作用下,大衛(wèi)星液滴兩側(cè)的小衛(wèi)星液滴迅速與左右兩個大液滴融合,最后形成兩個完整的液滴.

圖3 計算結(jié)果(We=149, Re=135.52, Ur=44 m/s)Fig. 3 Calculation results (We=149, Re=135.52, Ur=44 m/s)

4.2.1 玻璃微珠液滴對心碰撞過程中的能量變化

許多研究表明在液滴碰撞過程中,存在3種能量即液滴動能Ek、表面能Es和黏性耗散能Ev在發(fā)生變化[17-18].為研究玻璃微珠液滴對心碰撞過程中的能量變化情況,本文對上述3種能量進(jìn)行計算分析.液滴碰撞過程能量守恒可表示為

Ek0+Es0=Ekn+Esn+Ev,

(12)

式中Ek0,Es0分別是液滴初始動能和初始表面能;Ekn,Esn和Ev分別是液滴碰撞過程中的動能、表面能和黏性耗散能.初始狀態(tài)下的動能和表面能計算公式如下:

(13)

Es0=σS=4πσR2,

(14)

式中ρ為液滴密度;R為液滴半徑;v0=1/(2Ur),Ur為液滴碰撞的相對速度;σ為表面張力系數(shù).

碰撞過程中的動能和表面能需要計算液滴的體積和表面積,由于模擬環(huán)境為二維平面,無法直接得到液滴的體積和表面積.為此本文使用當(dāng)量半徑Re,它是由現(xiàn)有的二維平面的液滴面積演化而來的.我們假設(shè)二維平面的液滴面積S為三維液滴的中心截面面積,并且把它等量成球的中心截面Se,則該球的半徑即為當(dāng)量半徑Re,計算公式為

(15)

液滴碰撞過程中的動能計算公式為

(16)

式中v為液滴速度.

液滴碰撞過程中的表面能計算公式為

(17)

液滴碰撞過程中的黏性耗散能計算公式為

Ev=Ek0+Es0-Ekn-Esn.

(18)

由上述公式計算出各個時刻相應(yīng)的能量,并將它們與初始總能量即Ek0+Es0相除進(jìn)行無量綱化處理,得到We=149時玻璃微珠液滴碰撞分離過程的能量變化曲線,如圖4所示.圖中Es,Ek和Ev曲線分別代表表面能、動能和黏性耗散能的變化情況.從圖中可以看出,在玻璃微珠液滴相互靠近,突破氣膜阻力,碰撞后軸向收縮,徑向拉伸的過程中,動能先迅速下降,在液滴即將拉伸到極限時,動能緩慢減小,在2.0 μs動能減小至零;由于液滴融合,表面積減小,表面能緩慢下降,在2.0 μs降到最小;由于動能和表面能減小,黏性耗散能持續(xù)增大,在2.0 μs增至最大.2.0~2.5 μs動能和表面能基本保持不變,這是因為液滴軸向繼續(xù)收縮,徑向變化不大,在2.5 μs軸向收縮到極限.2.5 μs后,液滴徑向收縮,軸向拉伸,在4.0 μs徑向收縮和軸向拉伸達(dá)到平衡,11.5 μs液滴發(fā)生分離,產(chǎn)生3個衛(wèi)星液滴.在這個過程中,動能先增大,當(dāng)徑向收縮和軸向拉伸達(dá)到平衡時,動能增至最大,隨后在表面張力的作用下,液滴徑向收縮,軸向拉伸形成液橋,在二者的持續(xù)作用下,液橋逐漸變細(xì),在12.0 μs被夾斷,在此期間,動能緩慢減少至零.液滴在2.5 μs后,由于徑向收縮、軸向拉伸,表面積增大,表面能逐漸增大,在6.0 μs液滴軸向拉伸到極限時增加到最大值.隨后由于液滴徑向不斷收縮,液橋變細(xì),導(dǎo)致表面能小幅減?。谝簶虮粖A斷,形成3個大小不一的衛(wèi)星液滴后,由于小衛(wèi)星液滴和大液滴的融合,表面能發(fā)生小幅振蕩,最后和初始表面能相當(dāng),符合VOF方法的質(zhì)量守恒特征.黏性耗散能由于動能和表面能的增加,急劇下降,在液滴收縮拉伸達(dá)到平衡后緩慢增長.

通過上述分析可以得出,液滴動能是前期液滴發(fā)生碰撞融合的主要動力,在整個碰撞過程中對于液滴的收縮拉伸變化有著重要作用,在收縮拉伸過程中,動能會出現(xiàn)峰值.液滴在碰撞后徑向拉伸過程表面能減小;在軸向拉伸至分離過程中,表面能先增大后減小,恢復(fù)至初始表面能水平.液滴碰撞過程中黏性耗散能變化受動能和表面能影響,但從最終結(jié)果來看,液滴初始動能絕大部分轉(zhuǎn)化為黏性耗散能.由此可以看出,液滴碰撞過程中的分離結(jié)果受動能的起伏變化和表面能的變化的影響較大,在軸向拉伸前期,動能增大,表面能增大,為后面液滴拉伸分離提供了能量支持.

4.2.2 玻璃微珠液滴總壓變化

玻璃微珠液滴在發(fā)生碰撞的過程中,其靜壓和動壓都在發(fā)生變化,只分析其中一個都比較片面.而總壓為兩者之和,可以較為全面地反映液滴的壓力變化規(guī)律.圖5為兩液滴碰撞過程液滴平均總壓變化折線圖.從圖中可以看出, 兩液滴從零時刻開始相互靠近,壓縮氣膜,液滴總壓升高.0.5~2.0 μs為兩液滴碰撞融合,并徑向拉伸到極限的過程,液滴總壓持續(xù)下降,在2.0 μs降至最低.2.0~4.0 μs為液滴徑向收縮,軸向拉伸達(dá)到平衡的過程,液滴總壓持續(xù)增大,在4.0 μs達(dá)到平衡時增至最大.4.0~6.0 μs為液滴徑向收縮,軸向拉伸到極限的過程,液滴總壓減?。?.0 μs之后,液滴軸向不再拉伸,徑向仍在不斷收縮,中間液橋部分形態(tài)變化均勻,因此液滴總壓保持不變.在11.5 μs時,可以明顯看出液橋形態(tài)發(fā)生變化,兩端變細(xì),液滴總壓也忽然增大,這與末端夾斷機制一致[19-20].液橋斷裂之后,小衛(wèi)星液滴和大液滴發(fā)生融合,形成穩(wěn)定形態(tài)的過程中,液滴總壓降至穩(wěn)定狀態(tài).

4.2.3 玻璃微珠壓力云圖和速度矢量圖分析

從上述玻璃微珠液滴能量變化和總壓變化分析可知,液滴碰撞分離過程有4個重要狀態(tài),即:2.0 μs徑向拉伸到極限,4.0 μs徑向收縮和軸向拉伸達(dá)到平衡狀態(tài),6.0 μs軸向拉伸到極限和11.5 μs液滴液橋夾斷分離.下面對這4種狀態(tài)的壓力云圖和速度矢量圖進(jìn)行分析.圖6給出的是兩相同尺寸玻璃微珠液滴發(fā)生碰撞2.0 μs時,即液滴徑向拉伸到極限的速度矢量圖(圖6(a))(紅色箭頭代表液相速度矢量,藍(lán)色箭頭代表氣相速度矢量,后同)和液滴內(nèi)部壓力云圖(圖6(b)).從液滴的速度矢量圖可以看出液滴的速度分布上下對稱,上方的液滴拉伸部分速度向上,上端部分速度向下擴(kuò)散;外部的速度在液滴周圍呈大小不一的漩渦擴(kuò)散,總體趨勢向上,但上端部分漩渦有向下壓縮液滴的趨勢.下方液滴速度分布和上方相同,但方向相反,具有很好的對稱性.液滴拉伸部分向上下兩側(cè)擴(kuò)散的速度矢量箭頭密集,但兩端的速度矢量箭頭稀疏,這是由于啞鈴形玻璃微珠液滴中間部分較細(xì),其界面曲率κ較大,在此處的表面張力較大,即式(2)中的動力源項F較大,使此處產(chǎn)生指向液滴內(nèi)部的較大速度.這也表示液滴徑向拉伸已達(dá)到極限,沒有持續(xù)拉伸的趨勢了,與此時刻之后液滴的形態(tài)變化相符.此時的液滴形態(tài)呈啞鈴形,液滴壓力云圖與之相符,且壓力是從兩端向中心逐漸減小,這和圖5中此時刻的平均總壓變化相符.

圖4 玻璃微珠液滴碰撞分離過程能量變化曲線

(a) 速度矢量圖 (b) 壓力云圖 (a) Velocity vectors (b) Pressure contours

圖7為玻璃微珠液滴徑向收縮和軸向拉伸達(dá)到平衡時刻4.0 μs的速度矢量圖和內(nèi)部壓力云圖.從圖7(a)可以看出,此時液滴的速度大多集中在左右兩端,并分別向兩端擴(kuò)散.液滴兩端的速度矢量箭頭明顯比較密集,這是由于液滴形態(tài)呈鈴鐺形狀,中間左右兩端部分突出,其界面曲率κ就較大,則在此處的表面張力,即式(2)中的動力源項F也就越大,使此處液體的速度越大.圖7(b)為玻璃微珠液滴壓力云圖,圖中顯示液滴壓力從上下兩端向中心逐漸增大,從中心向左右逐漸增大.液滴壓力和速度的分布的共同作用促成了此刻的液滴形態(tài).

圖8為玻璃微珠液滴軸向拉伸到極限時刻6.0 μs的速度矢量圖和內(nèi)部壓力云圖.從圖8(a)可以看出,此時液滴的速度分布均勻,并且從中間分別向左右兩端擴(kuò)散,具有良好的對稱性.此時的液滴形態(tài)呈柱狀,形態(tài)勻稱,界面曲率κ沒有明顯變化,所以液滴速度分布均勻.在端部分布有速度漩渦,與液滴內(nèi)部速度方向相反,阻礙液滴的軸向拉伸.圖8(b)為玻璃微珠液滴壓力云圖,兩端的壓力比中部壓力大,但兩者相差不大.速度矢量圖和壓力云圖都與液滴的形態(tài)變化特性相符.

(a) 速度矢量圖 (b) 壓力云圖 (a) Velocity vectors (b) Pressure contours

(a) 速度矢量圖 (b) 液滴內(nèi)部壓力云圖 (a) Velocity vectors (b) Pressure contours

(a) 速度矢量圖(a) Velocity vectors

(b) 液橋局部放大速度矢量圖 (c) 壓力云圖 (b) Partially amplified velocity vectors of the liquid bridge (c) Pressure contours

圖9為兩相同尺寸玻璃微珠液滴碰撞分離時刻11.5 μs的速度矢量圖和壓力云圖.從圖9(a)可以看出,此時液滴的速度大多集中在液橋部分,液滴形態(tài)呈啞鈴狀.圖9(b)為圖9(a)液橋部分速度矢量圖的放大圖,從中可以清晰地看出液橋部分的速度從兩端向中部擴(kuò)散,液橋外部末端靠近大液滴處有速度漩渦形成,速度從中心向兩端擴(kuò)散.圖9(c)為玻璃微珠液滴壓力云圖,圖中顯示液滴壓力都集中在液橋的部位,并從液橋中心向左右兩邊增大,導(dǎo)致末端壓力最大.綜合速度和壓力分析可知兩者的分布與液滴的末端夾斷機制一致,這也是導(dǎo)致液滴碰撞分離的重要原因.

5 結(jié) 論

本文對玻璃微珠液滴碰撞分離過程進(jìn)行了數(shù)值模擬研究,闡述了計算方法,并驗證了計算方法可靠性,主要結(jié)論如下:

1) 在碰撞過程中,液滴的初始表面能和最終表面能相當(dāng),動能則絕大部分轉(zhuǎn)化為黏性耗散能.

2) 在液滴軸向拉伸分離前期,動能增大、表面能增大為后期液滴分離供能,對玻璃微珠液滴碰撞后分離有著重要作用.

3) 通過對玻璃微珠液滴能量和平均總壓變化的研究,得到了液滴碰撞分離過程有4種狀態(tài),即徑向拉伸到極限、徑向收縮和軸向拉伸達(dá)到平衡、軸向拉伸到極限和液滴液橋夾斷分離.

4) 液滴的速度分布與液滴的形態(tài)密切相關(guān),不規(guī)則的液滴形態(tài)導(dǎo)致界面曲率κ較大,以此影響界面的表面張力,使得不規(guī)則的液體部分產(chǎn)生較大的不同指向的速度分布.

5) 末端夾斷機制是玻璃微珠液滴碰撞分離的主要原因.

主站蜘蛛池模板: 免费欧美一级| 欧美19综合中文字幕| 午夜视频免费一区二区在线看| 久久久噜噜噜久久中文字幕色伊伊 | 日本在线亚洲| 国产va免费精品观看| 国产女同自拍视频| 久久黄色视频影| 69免费在线视频| 91破解版在线亚洲| 国产不卡一级毛片视频| 亚洲精品无码人妻无码| 视频二区国产精品职场同事| 成年免费在线观看| 国产原创演绎剧情有字幕的| 国产成人乱无码视频| 2020最新国产精品视频| JIZZ亚洲国产| 九色综合视频网| 国产亚洲精品无码专| 91精品伊人久久大香线蕉| 亚洲AV无码久久天堂| 少妇高潮惨叫久久久久久| 久久久久国产一级毛片高清板| 国产女人18水真多毛片18精品| 国产精品开放后亚洲| 91av国产在线| 激情無極限的亚洲一区免费| 全部毛片免费看| 精品国产福利在线| 久久综合婷婷| 国产成人av一区二区三区| 欧美高清视频一区二区三区| 91色老久久精品偷偷蜜臀| 亚洲中文字幕手机在线第一页| 91久久大香线蕉| 亚洲成人免费在线| 亚洲欧美另类日本| www.狠狠| 亚洲免费黄色网| 成人国产精品网站在线看| 国产后式a一视频| 日韩中文精品亚洲第三区| 无码区日韩专区免费系列| 欧美精品v| 亚洲中文在线看视频一区| 国产全黄a一级毛片| 亚洲日本www| 黄色网站不卡无码| 国产成人精品男人的天堂下载 | 国产精品美女在线| 免费黄色国产视频| 国产乱人伦偷精品视频AAA| 亚洲三级成人| 在线看免费无码av天堂的| 国产永久无码观看在线| 久久久久久尹人网香蕉| 亚洲免费毛片| av天堂最新版在线| 欧美日韩91| 看国产毛片| 欧美日本中文| 国产农村1级毛片| 亚洲水蜜桃久久综合网站| 激情成人综合网| 成人免费一区二区三区| 亚洲高清无码久久久| 亚洲综合色婷婷中文字幕| 亚洲成a人片7777| jizz在线免费播放| 亚洲视频影院| 无码中文AⅤ在线观看| 中文字幕亚洲无线码一区女同| 国产乱码精品一区二区三区中文 | 97在线观看视频免费| 青青草原偷拍视频| 国产亚洲欧美另类一区二区| 玖玖免费视频在线观看| 无码电影在线观看| 欧美亚洲一区二区三区导航| 亚洲青涩在线| 国产第二十一页|