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

循環(huán)射流混合槽壓力波動時間序列差異性分析

2016-09-16 01:20:36孟輝波趙耀禹言芳宋明遠吳劍華
哈爾濱工程大學學報 2016年8期
關(guān)鍵詞:符號

孟輝波, 趙耀, 禹言芳, 宋明遠 , 吳劍華

(1. 沈陽化工大學 遼寧省高效化工混合技術(shù)重點實驗室, 遼寧 沈陽 110142; 2. 天津大學 化工學院, 天津 300072)

?

循環(huán)射流混合槽壓力波動時間序列差異性分析

孟輝波1,2, 趙耀1, 禹言芳1,2, 宋明遠1, 吳劍華1,2

(1. 沈陽化工大學 遼寧省高效化工混合技術(shù)重點實驗室, 遼寧 沈陽 110142; 2. 天津大學 化工學院, 天津 300072)

為了研究循環(huán)射流混合槽內(nèi)壓力波動信號的前向和后向符號時間序列的動力學特征的差異性,利用等概率原則對PFS(pressure fluctuation signals)時間序列進行符號化轉(zhuǎn)換,通過修正的Shannon熵選取最佳符號集大小和子序列長度。利用子序列編碼圖、時間不可逆轉(zhuǎn)性、秩次圖和秩次距離等參數(shù)對PFS時間序列進行STSA(symbolic time series analysis)研究。研究結(jié)果表明:修正的Shannon熵最小時確定優(yōu)化符號化參數(shù)為n=2和L=10。前向與后向時間序列的子序列編碼分布相似但其頻數(shù)不等,表明CJT(circulating jet tank)內(nèi)湍流流動呈現(xiàn)多尺度混沌確定性特征。PFS的時間不可逆性值隨著周向角的增加不斷增加,隨z/H的增大呈現(xiàn)先降低再上升最后降低的趨勢,隨雷諾數(shù)的增加呈現(xiàn)"W"型分布。時間不可逆性值與秩次距離對PFS前后向序列動力學特征的差異性判斷相互吻合。

循環(huán)射流混合槽;壓力波動;符號時間序列分析;時間不可逆性;秩次距離;湍流;混沌

網(wǎng)絡出版地址:http://www.cnki.net/kcms/detail/23.1390.U.20160624.1127.022.html

循環(huán)射流混合槽作為一種化工過程強化的單元操作設備,有效地解決了聚氯乙烯(Polyvinyl chloride,PVC)傳統(tǒng)生產(chǎn)槽底部密封泄漏以及高能耗難題,有效改善槽內(nèi)混合效果,增加有效生產(chǎn)時間,提高生產(chǎn)效率,廣泛應用于精細化工生產(chǎn)中的聚合反應過程[1-2]。壓力波動信號(pressure fluctuation signals, PFS)作為湍流流動的動力學顯著指標,它隱含著系統(tǒng)內(nèi)部復雜的流體動力學信息[3]。本課題組前期對循環(huán)射流混合槽(circulating jet tank,CJT)內(nèi)瞬態(tài)脈動壓力的動力學特性的研究取得一定進展[4-10],但是對PFS時間序列的動力學特征的相似性或差異性研究未見相關(guān)文獻報道。

符號時間序列分析(symbolic time series analysis,STSA)作為一種由混沌時間序列分析、信息理論和符號動力學理論相結(jié)合的數(shù)據(jù)分析工具,具有對噪聲不敏感的性質(zhì),有效降低了噪聲對流動固有特征的影響,因此被廣泛應用于化工、機械、經(jīng)濟學、天文學、地理學等不同領(lǐng)域[11-15]。

因此,本文利用等概率區(qū)間法將不同軸向、周向和徑向位置的PFS轉(zhuǎn)化成符號時間序列,分別分析符號時間序列的直方圖、時間不可逆性Tfb、秩次距離的分布特性,從而為進一步提取CJT內(nèi)壓力波動信號固有特征提供一種思路。

1 符號時間序列基礎理論

1.1符號時間序列分析

采用等概率法將PFS時間序列P(iΔt) (i=1, 2, …,N)轉(zhuǎn)化為符號序列S(iΔt)[11]:

(1)

式中:N為信號長度,n為符號集大小,F(xiàn)i為分位數(shù)閾值。

在符號序列S(iΔt)中任意連續(xù)選取長度為L的符號串,構(gòu)成子序列。子序列在S(iΔt)中出現(xiàn)的次數(shù),稱為該子序列頻數(shù)k。確定符號集大小n和子序列長度L后劃分為N-L+1個子序列,子序列再經(jīng)過十進制編碼轉(zhuǎn)換為編碼序列X。Finney等[14]定義修正的Shannon熵:

(2)

式中:Nobs為S(iΔt)中出現(xiàn)頻數(shù)不為0的子序列數(shù)量,Nobs≤K。K為子序列十進制編碼總數(shù),K=nL。j為子序列編碼,pj,L是長度為L第j個子序列編碼出現(xiàn)的頻率。對于完全隨機序列,所有子序列頻數(shù)相同,H(L)=1。對于完全確定過程產(chǎn)生的序列,只出現(xiàn)某個子序列,則H(L)=0。

1.2時間不可逆轉(zhuǎn)性Tfb

按照正向和反向?qū)⒃紩r間序列數(shù)據(jù)轉(zhuǎn)換得到兩個不同的符號序列,分別求出兩個序列中子序列出現(xiàn)頻率,再進行歐氏范數(shù)計算,便可得到時間不可逆轉(zhuǎn)性Tfb:

(3)

式中:pf,j和pb,j分別表示前向與后向序列中第j個子序列出現(xiàn)頻率。Tfb值越大說明前向和后向符號序列差別越大,時間不可逆性程度越大,表明系統(tǒng)越復雜,不確定性越大。反之,如果系統(tǒng)呈現(xiàn)出良好規(guī)律性變化,那么Tfb趨于0。

1.3符號序列秩次圖與秩次距離

同一個子序列在分別降序排序的兩個不同時間序列的秩次為R1,j和R2,j,可以組成一組對應數(shù)據(jù)對(R1,j,R2,j),將所有子序列秩次點在坐標圖中描繪出來,得到兩個時間序列對應符號序列秩次圖。兩個時間序列若相似,秩次點將分布在靠近圖中對角線位置上。而秩次距離[15]是按照同一子序列編碼在不同時間序列中出現(xiàn)秩次來計算兩個序列間距離:

(4)

式中:x、y分別為兩個時間序列,Xm為符號序列中第m個編碼,Rx(Xm)和Ry(Xm)為Xm分別在x、y中秩次,px(Xm)和py(Xm)為Xm對應編碼子序列分別在x、y中出現(xiàn)頻率。DL的范圍為0~1,其中DL越大表明兩個時間序列間相似性越小,差異越大;反之,DL越小表明兩個時間序列相似性較高。當DL為0時,表明兩個時間序列完全一致。

2 CJT內(nèi)瞬態(tài)PFS實驗

CJT內(nèi)瞬態(tài)PFS時間序列測量裝置示意圖參見文獻[1-2],實驗工質(zhì)為蒸餾水,循環(huán)系統(tǒng)動力源采用德國Wilo(MHI802)不銹鋼多級離心泵,用不銹鋼金屬管浮子流量計(LZD-50/Y10/RR1/ESK型)控制實驗循環(huán)流量Q在1~9 m3/h范圍內(nèi),相應的Re為

(5)

式中:ρ是實驗工質(zhì)的密度,kg/m3;dj是射流孔直徑,m;μ是流體粘度,Pa·s;N0是射流孔總個數(shù)。

Re變化范圍為3 660~32 940,可知CJT內(nèi)流體流動為湍流流動[16]。

前期研究結(jié)果表明PFS的最大波動頻率為500 Hz[4],因此本實驗設定采樣頻率fs為1 kHz?;贑JT的軸對稱結(jié)構(gòu),選取整個CJT的1/4作為測量區(qū)域。如圖1所示沿壁面周向和軸向設置20個測量孔,通過調(diào)節(jié)傳感器標尺,可使每個測量孔內(nèi)沿徑向距離r/R在0.685~1.0范圍內(nèi)均勻布置10個測量點。

圖1 CJT內(nèi)瞬態(tài)壓力波動信號測量點布置圖Fig.1 Locations of measurement points for instantaneous

3 瞬態(tài)PFS符號差異性分析與討論

3.1符號集大小n和子序列長度L的選擇

對應特定時間序列當H(L)值最小時,對應n與L就是最佳參數(shù)。圖2為PFS時間序列修正的Shannon熵隨符號集大小n和子序列長度L的變化規(guī)律。從圖2可看出,不同符號集下修正的Shannon熵曲線隨著L的增加呈現(xiàn)先降低再上升的趨勢。但是,隨著L的進一步增加,系統(tǒng)的隨機成分占據(jù)主導地位,H(L)值趨近于1。通過計算不同軸向、周向和徑向位置的序列的H(L)值,發(fā)現(xiàn)當符號集大小n=2且L為10~11時得到修正的Shannon熵值比其他情況下都要小。因此,選取符號集n=2,作為本文的最佳符號集大小。

圖2 修正的Shannon熵與符號化參數(shù)的關(guān)系Fig.2 Modified Shannon entropy versus symbolization parameters

子序列長度L在8~12范圍內(nèi)時H(L)值與其真實最小值之間誤差E的變化規(guī)律如圖3所示。從圖3可以看出,隨著L的增加修正的Shannon熵的誤差先減小再增大,當L=8時修正的Shannon熵的誤差范圍為0.28~2.09%;當L為9、10、11時,最大誤差分別為0.91%、0.17%、0.39%;當L=12時誤差范圍為0.29~0.89%。因此選取L=10為最佳的子序列長度。

圖3 不同L的修正Shannon熵的誤差曲線Fig.3 Modified Shannon entropies under different L of error curves

3.2壓力波動信號的直方圖

圖4為不同周向位置處的PFS前向、后向時間序列子序列編碼分布圖。從圖可以看出CJT內(nèi)的PFS前向與后向時間序列子序列編碼頻數(shù)分布相似。由圖4還能看出,所有直方圖中子序列編碼頻數(shù)大小不等,說明CJT內(nèi)復雜的動力學結(jié)構(gòu)時間序列不是隨機的,具有一定的確定性。在θ為π/6和π/4處,時間序列中子序列編碼頻數(shù)主要集中在編碼0和1 023,兩者頻數(shù)高于3 000,說明此處PFS波動奇異性較高。在周向位置θ=π/3處,頻數(shù)相對較大且超過200次的編碼個數(shù)為13個;在置θ=5π/12處,頻數(shù)相對較大且超過200次的編碼個數(shù)為21個。由此看出π/3和5π/12處,時間序列中子序列出現(xiàn)的頻數(shù)分布隨機性增強。結(jié)合前期PIV實驗及LES瞬態(tài)模擬分析發(fā)現(xiàn)[7-9],水平浸沒射流與周圍流體相互摻混卷吸誘導形成多尺度流體動力學結(jié)構(gòu),并且漩渦尺度隨著射流長度即θ值的增大而增大。從而側(cè)面驗證壓力波動多尺度結(jié)構(gòu)的存在[6]。

3.3PFS時間不可逆轉(zhuǎn)性Tfb

圖5為CJT內(nèi)PFS符號化時間序列的Tfb變化曲線。圖5(a)揭示了Tfb值隨著周向角度增大不斷增加。結(jié)合圖1所示射流結(jié)構(gòu)發(fā)現(xiàn):當θ=π/6時,測量點靠近射流起始段,此區(qū)域射流軸中心未受湍動摻混影響流體流動比較穩(wěn)定;當θ=5π/12時測量點位于射流主體沖擊圓筒內(nèi)壁形成壁面射流與循環(huán)混合主體區(qū)內(nèi)截面二次流[8-10]相互耦合摻混,湍流PFS脈動比較激烈,Tfb達到最大值即PFS時間序列的時間不可逆性增強。不同軸向位置下Tfb分布如圖5(b)所示,隨z/H的減小CJT內(nèi)的流體受多股水平浸沒射流的相互卷吸合并耦合作用,導致槽內(nèi)流體波動性和不穩(wěn)定性增強,其中z/H=0.45處Tfb值約為z/H=0.65處1倍。徑向位置r/R在0.685~1處的Tfb的變化規(guī)律如圖5(c)所示。從圖中看出,Tfb隨r/R的增大呈現(xiàn)先上升后降低的近似對稱分布趨勢。在徑向位置r/R=0.860處的測量點位于射流軸線附近,此區(qū)域內(nèi)流體保留射流出口速度能力最大,流體微團的壓力脈動特性最強,致使Tfb值達到最大值;Tfb不完全對稱分布是由于射流中心線兩側(cè)(左側(cè)為一維平面擋板,右側(cè)為二維圓柱曲面)的邊界條件差異引起;近壁區(qū)的Tfb數(shù)值較小主要是因為在壁面附近的邊界層消弱流體微團的脈動,導致波動穩(wěn)定性增強。如圖5(d)所示,Tfb值隨著雷諾數(shù)的增加呈現(xiàn)“W”型趨勢,變化范圍為0.006~0.017;由圖還可看出在Re=10 980時CJT內(nèi)PFS的波動特性穩(wěn)定;然而,當Re≥29 332時擋板與圓筒壁面形成固體壁面對射流中心高速流體卷吸周圍低速流體能力的限制作用逐步加強,致使Tfb值趨于穩(wěn)定。上述研究表明,CJT內(nèi)流體流動混合呈現(xiàn)整體混沌確定性結(jié)構(gòu),與前期研究結(jié)果相吻合[1,6]。

圖4 不同周向位置下PFS前向和后向時間序列的直方圖Fig.4 Forward and backward sequence histograms of PFS time series at different circumferential positions

圖5 時間不可逆轉(zhuǎn)性Tfb在不同雷諾數(shù)及位置下的變化曲線Fig.5 Temporal irreversibility Tfb under different Re and positions

3.4符號序列秩次圖與秩次距離

圖6為不同周向位置處PFS的前向和后向時間序列的秩次圖,圖中對角線方向秩次點代表前后向時間序列的共同特征,坐標系中最左下部分即秩次點(256,256)以內(nèi)點代表CJT內(nèi)湍流射流混合引起的壓力波動主要變化特征,圖中右上部分秩次點代表CJT內(nèi)湍流射流混合引起的壓力波動突變特征;圖中平行于坐標軸的直線表明前向和后向子序列的頻數(shù)出現(xiàn)大量相同的情況,與圖4頻數(shù)高峰分布相吻合。由于CJT內(nèi)θ在π/6和π/4處測量點均位于射流初始階段,致使PFS秩次點(64,64)以內(nèi)分布的點比π/3和5π/12更靠近對角線,造成θ=π/6和π/4時各自的PFS的前向和后向的動力學特征相似性較高。如圖6(d)所示,當θ=5π/12時后向序列的主要變化與前向序列具有較強的差異性。同時隨著射流沿程的增加,射流沖擊圓柱壁面,射流軸線左側(cè)部分與循環(huán)混合主體區(qū)耦合摻混,而射流軸線右側(cè)部分形成壁面射流二次流,上述兩種不同流動結(jié)構(gòu)[17-19]誘導流體微團壓力波動的奇異性逐漸增多。

不同周向位置下秩次距離DL如圖7所示,DL值隨著θ的增加不斷增大。從圖7與圖6的對比可以看出,DL值的分布規(guī)律與秩次圖對前向與后向序列差異性做出判斷相吻合,并與不可逆參數(shù)分布規(guī)律相一致,說明符號序列秩次圖和秩次距離能夠較好地對時間序列間差異性做出定性和定量地判斷,可作為一種有效提取PFS時間序列湍動特性的統(tǒng)計分析手段。

圖6 不同周向位置下前向和后向時間序列的秩次圖Fig.6 Rank order diagrams of forward and backward sequence of PFS time series at different circumferential positions

圖7 不同周向位置下前向和后向時間序列的秩次距離DL Fig.7 Rank order distances of forward and backward sequence of PFS time series under different circumferential positions

4 結(jié)論

對循環(huán)射流混合槽內(nèi)壓力波動信號的前向和后向符號時間序列的動力學特征的差異性研究得到以下結(jié)論:

1)利用等概率原則對CJT內(nèi)PFS時間序列進行符號化轉(zhuǎn)換,基于修正的Shannon熵優(yōu)化得到最佳符號集大小n=2,子序列長度L=10。

2)時間不可逆轉(zhuǎn)性Tfb值隨著射流沿程的增大不斷增加。θ=π/6時,測量點靠近射流噴嘴起始段,由于射流中心部分流體與周圍低速流體摻混較弱致使PFS波動比較穩(wěn)定。當θ=5π/12時,射流沖擊圓筒內(nèi)壁形成壁面射流與循環(huán)混合主體區(qū)流體摻混導致其PFS時間序列的時間不可逆性增強。隨z/H的減小CJT內(nèi)的流體受更多股水平浸沒射流的相互卷吸合并等耦合作用,導致槽內(nèi)流體波動性和不穩(wěn)定性增強,其中z/H=0.45處Tfb值約為z/H=0.65處1倍。Tfb值隨著雷諾數(shù)的增加呈現(xiàn)“W”型趨勢。

3)符號序列秩次圖與秩次距離能夠較好地對PFS時間序列間差異性做出定性和定量地判斷,DL與Tfb變化趨勢在周向和軸向位置下大體相同。DL值分布在0~0.016,表明CJT內(nèi)湍流壓力波動在整體上呈現(xiàn)出確定性結(jié)構(gòu)特征。

[1]MENG Huibo, WANG Feng, YU Yanfang, et al. Multiscale entropy analysis of instantaneous pressure fluctuation in a novel jet tank[J]. Chemical engineering & technology, 2013, 36(12): 2137-2147.

[2]MENG Huibo, WANG Wei, WU Jianhua, et al. Experimental study on instantaneous pressure fluctuation time series in the novel tank agitated by multiple horizontal jets[J]. Chemical engineering research and design, 2012, 90(11): 1750-1764.

[3]BI H T. A critical review of the complex pressure fluctuation phenomenon in gas-solids fluidized beds[J]. Chemical engineering science, 2007, 62(13): 3473-3493.

[4]MENG Huibo, SONG Mingyuan, YU Yanfang, et al. Recurrence quantity analysis of the instantaneous pressure fluctuation signals in the novel tank with multi-horizontal submerged jets[J]. Chemical and biochemical engineering quarterly, 2016, 30(1): 19-31.

[5]MENG Huibo, WANG Feng, YU Yanfang, et al. Delay time correlation of pressure fluctuation signals in the novel circulating jet tank[J]. Chemical and biochemical engineering quarterly, 2013, 27(3): 251-257.

[6]MENG Huibo, WANG Yanfen, YU Yanfang, et al. Analysis of pressure fluctuations induced by multi-horizontal submerged jets in the novel jet tank[J]. The Canadian journal of chemical engineering, 2014, 92(5): 935-944.

[7]YU Yanfang, JIANG Xiuhui, MENG Huibo, et al. Computational simulation of mixing performance in the circulating jet mixing tank[J]. International journal of chemical reactor engineering, 2016, 14(2): 621-636.

[8]王艷芬. 新型循環(huán)射流混合槽內(nèi)瞬態(tài)流場PIV實驗研究及噴嘴形狀優(yōu)化[D]. 沈陽: 沈陽化工大學, 2014: 18-32.

WANG Yanfen. PIV experimental research of instantaneous flow characteristics and optimization of nozzle shape in novel Circulating Jet Tank[D]. Shenyang: Shenyang University of Chemical Technology, 2014: 18-32.

[9]MENG Huibo, WANG Wei, YU Yanfang, et al. Investigation of the effect of outlet structures on the jet flow characteristics in the circulating jet tank[J]. International journal of chemical reactor engineering, 2014, 12(1): 35-45.

[10]YU Yanfang, WU Jianhua, MENG Huibo. Numerical simulation process aspects of the novel static circulating jet mixer[J]. The Canadian journal of chemical engineering, 2011, 89(3): 460-468.

[11]CRUTCHFIELD J P, PACKARD N H. Symbolic dynamics of noisy chaos[J]. Physica D, 1983, 7(1/2/3): 201-223.

[12]LEHRMAN M, RECHESTER A B, WHITE R B. Symbolic analysis of chaotic signals and turbulent fluctuations[J]. Physical review letters, 1997, 78(1): 54-57.

[13]DAW C S, FINNEY C E A, KENNEL M B. Symbolic approach for measuring temporal “irreversibility”[J]. Physical review E, 2000, 62(2): 1912-1921.

[14]FINNEY C E A, GREEN J B JR, DAW C S. Symbolic time-series analysis of engine combustion measurements[R]. SAE Technical Paper 980624, 1998.

[15]SUYAL V, PRASAD A, SINGH H P. Symbolic analysis of slow solar wind data using rank order statistics[J]. Planetary and space science, 2012, 62(1): 55-60.

[16]ZUGHBI H D, RAKIB M A. Investigations of mixing in a fluid jet agitated tank[J]. Chemical engineering communications, 2002, 189(8): 1038-1056.

本文引用格式:

孟輝波, 趙耀, 禹言芳,等. 循環(huán)射流混合槽壓力波動時間序列差異性分析[J]. 哈爾濱工程大學學報, 2016, 37(8): 1157-1162.

MENG Huibo, ZHAO Yao, YU Yanfang,et al. Difference analysis of pressure fluctuation time series in the circulating jet tank[J]. Journal of Harbin Engineering University, 2016, 37(8): 1157-1162.

Difference analysis of pressure fluctuation time series in the circulating jet tank

MENG Huibo1,2, ZHAO Yao1, YU Yanfang1,2, SONG Mingyuan1, WU Jianhua1,2

(1. Liaoning Key Laboratory of Chemical Technology for Efficient Mixing, Shenyang University of Chemical Technology, Shenyang 110142, China; 2. School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China)

In this study, we investigated the differences in the flow dynamic characteristics in forward and backward time series of pressure fluctuation signals (PFS) in the circulating jet tank (CJT). We used the principle of equivalent possibility to symbolically transform the PFS time series. We also used modified Shannon entropy to choose the optimal parameters, including symbolic set size n and subsequence length L. We performed a symbolic time series analysis (STSA) for the PFS time series using a sequence coding histogram, time irreversibility, a rank order figure, and rank order distance. The results show that the symbolic set size n=2 and subsequence length L=10 achieved the minimum of modified Shannon entropy. The non-uniform frequency distributions of forward and backward sequence codes are similar, which indicates the existence of multi-scale chaotic deterministic characteristics in the dynamical system of the CJT. The time irreversibility value increases with increasing circumference angles. First it decreases, then increases, and finally decreases again with increasing z/H. It also shows a W-type distribution with increasing the Reynolds number. Our difference evaluation results by time irreversibility have good agreement with the rank order distance of the forward and backward time series of the PFS.

circulating jet tank; pressure fluctuation; symbolic time series analysis; time irreversibility; rank order distance; turbulent flow; chaos

2015-06-10.網(wǎng)絡出版日期:2016-06-24.

國家自然科學基金項目(21476142, 21306115, 21106086); 遼寧省高等學校優(yōu)秀人才計劃(LR2015051); 遼寧省博士科研啟動基金項目(20131090); 遼寧省教育廳科研項目計劃(L2013164).

孟輝波(1981),男,副教授;

禹言芳, E-mail : taroyy@163.com.

10.11990/jheu.201506031

TQ051.7

A

1006-7043(2016)08-1157-06

猜你喜歡
符號
幸運符號
符號神通廣大
學符號,比多少
幼兒園(2021年6期)2021-07-28 07:42:14
“+”“-”符號的由來
靈魂的符號
散文詩(2017年17期)2018-01-31 02:34:20
怎樣填運算符號
變符號
倍圖的全符號點控制數(shù)
圖的有效符號邊控制數(shù)
草繩和奇怪的符號
主站蜘蛛池模板: 国产正在播放| 国产成人无码AV在线播放动漫 | 国产h视频免费观看| 9999在线视频| 国产97视频在线| 国产亚洲精品97AA片在线播放| 国产免费网址| 日韩不卡高清视频| 日韩欧美国产中文| 毛片免费在线视频| 国产熟睡乱子伦视频网站| 国产第一页屁屁影院| 最近最新中文字幕在线第一页| 日本亚洲欧美在线| 欧美天天干| 91久草视频| 九九热视频在线免费观看| 无码中文字幕乱码免费2| 激情综合网激情综合| 国产主播在线观看| 国产女人综合久久精品视| 亚洲精品午夜天堂网页| 成人在线不卡视频| 久爱午夜精品免费视频| 国产va在线观看| 亚洲AⅤ综合在线欧美一区| 色首页AV在线| 国产精品视频999| 欧美亚洲一区二区三区在线| 蜜桃臀无码内射一区二区三区 | 亚洲浓毛av| 在线国产三级| 欧美日韩精品在线播放| 热久久国产| 日韩在线中文| 又污又黄又无遮挡网站| 国产原创第一页在线观看| 99精品影院| 国产亚洲精| 青草娱乐极品免费视频| 国模私拍一区二区| 国产精品观看视频免费完整版| 国产亚洲精品91| 成人av手机在线观看| 亚洲无码在线午夜电影| 99热这里只有精品国产99| 久久综合五月| 午夜影院a级片| 日本影院一区| 91蜜芽尤物福利在线观看| 亚洲最猛黑人xxxx黑人猛交| av一区二区三区在线观看| 日本欧美一二三区色视频| 亚洲热线99精品视频| 国产91色| 性网站在线观看| 国产极品粉嫩小泬免费看| 国产成人综合久久| 免费看a级毛片| 久久窝窝国产精品午夜看片| 中国精品自拍| 国产swag在线观看| 日韩欧美在线观看| 丰满人妻一区二区三区视频| 伊在人亚洲香蕉精品播放 | 亚洲天堂免费在线视频| 国产精品第页| 日韩中文精品亚洲第三区| 88av在线看| 91探花国产综合在线精品| 亚洲黄色成人| 手机在线免费毛片| 亚洲人妖在线| 国产91线观看| 日韩欧美综合在线制服| 欧亚日韩Av| 精品国产美女福到在线直播| 成人va亚洲va欧美天堂| 精品久久久无码专区中文字幕| 欧美一级爱操视频| 日韩二区三区| 一本无码在线观看|