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

散射體旋轉(zhuǎn)角對(duì)二維聲子晶體帶隙結(jié)構(gòu)影響分析

2014-06-15 17:06:07杜敬濤賀彥博馮浩成

杜敬濤,賀彥博,馮浩成

(哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院,黑龍江哈爾濱150001)

散射體旋轉(zhuǎn)角對(duì)二維聲子晶體帶隙結(jié)構(gòu)影響分析

杜敬濤,賀彥博,馮浩成

(哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院,黑龍江哈爾濱150001)

為探討散射體旋轉(zhuǎn)作用對(duì)于彈性波在聲子晶體中的傳播特性產(chǎn)生的影響,采用平面波展開(kāi)法將散射體與晶格軸線夾角因素引入結(jié)構(gòu)函數(shù),研究了在不同填充率下,伴隨著旋轉(zhuǎn)角的變化,正方形、橢圓形、正六邊形散射體禁帶出現(xiàn)的位置及寬度變化情況。研究結(jié)果表明:散射體旋轉(zhuǎn)角度對(duì)二維聲子晶體帶隙結(jié)構(gòu)具有顯著影響,并且散射體旋轉(zhuǎn)角對(duì)正六邊形、正方形散射體的影響明顯大于對(duì)于橢圓形散射體的影響;然而,無(wú)論散射體取為何種形狀,旋轉(zhuǎn)角度對(duì)聲子晶體產(chǎn)生禁帶的位置(禁帶中間頻率)影響并不明顯。

聲子晶體;帶隙;禁帶;旋轉(zhuǎn)角;平面波展開(kāi)法;散射體

聲子晶體是指具有彈性波帶隙的周期性復(fù)合材料或結(jié)構(gòu),作為一種新型聲學(xué)功能帶隙材料,聲子晶體的研究不僅具有重要的理論價(jià)值,同時(shí)具有廣闊的應(yīng)用前景,可以為減振降噪提供一種全新的技術(shù)思路[1-2]。為此,聲子晶體已經(jīng)成為國(guó)際學(xué)術(shù)界的研究熱點(diǎn)。盡管對(duì)于彈性波在層狀介質(zhì)中的傳播特性的研究已有近70年的歷史,但是聲子晶體概念的提出及相關(guān)理論的研究時(shí)間并不長(zhǎng)。1992年,Sigalas等從理論上證實(shí)了周期結(jié)構(gòu)中彈性波禁帶現(xiàn)象的存在[3],隨后,Kushwaha等首次提出聲子晶體概念[4]。1995年,Martinez-Sala等從實(shí)驗(yàn)角度對(duì)彈性波禁帶特性進(jìn)行了觀測(cè)[5]。從此,聲子晶體的研究引起了全世界學(xué)者的關(guān)注。

就目前研究狀況而言,對(duì)于散射體在晶格中整齊排列理想情況的研究已經(jīng)獲得了較好的進(jìn)展。無(wú)論是二維或三維的情況,只有在一定條件下才能產(chǎn)生彈性波或聲帶波間隙。在進(jìn)行二組元二維聲子晶體的研究過(guò)程中,散射體排列方式多采用正方形[6]、六邊形[7]、氯化鑰排列[8]等。散射體截面積大部分采用圓形,少數(shù)則采用正方形。但以上研究都存在不足,即一旦晶體制備完成,在實(shí)際使用的過(guò)程中無(wú)法根據(jù)實(shí)際需要調(diào)整帶隙位置及寬度。

以上不足得到了中外學(xué)者的關(guān)注并提出了各自的解決方案。具體如K.Bertodi等提出了施加外應(yīng)力的方案,從而拓寬了聲子晶體帶隙[9-11]。此外,顏琳等通過(guò)散射體旋轉(zhuǎn)作用,從而成功地探究了聲子晶體帶隙寬度與旋轉(zhuǎn)角度的關(guān)系[12]。與此類似,吳福根也對(duì)于液體—液體環(huán)境下,正方形散射體偏轉(zhuǎn)角度對(duì)于帶隙結(jié)構(gòu)產(chǎn)生的影響進(jìn)行了研究[13]。但是他們的研究成果僅限于正方形散射體,對(duì)于其他形狀的散射體并未涉及。

基于散射體旋轉(zhuǎn)角對(duì)于聲子晶體禁帶的產(chǎn)生具有一定影響的事實(shí)將對(duì)于實(shí)際過(guò)程中人為調(diào)整聲子禁帶的位置及寬度具有極大的指導(dǎo)意義,本文在之前的基礎(chǔ)上進(jìn)行了更為深入的研究,包括正方形散射體六邊形散射體以及橢圓形散射體。

1 理論模型

依據(jù)彈性動(dòng)力學(xué)理論,在材料的連續(xù)性均勻性、完全線彈性、小變形、無(wú)初始應(yīng)力等基本假設(shè)下,分別建立描述質(zhì)點(diǎn)力、位移以及應(yīng)力應(yīng)變之間的三類基本方程[14],首先是運(yùn)動(dòng)微分方程:

式中:σ是應(yīng)力,ρ為質(zhì)量密度,u表示質(zhì)點(diǎn)的振動(dòng)位移。

描述應(yīng)變位移關(guān)系的幾何方程為

式中:ε表示應(yīng)變。描述應(yīng)力—應(yīng)變關(guān)系的物理方程為

此處采用了張量理論中的愛(ài)因斯坦規(guī)則。在求解時(shí),先將幾何方程(2)代入物理方程(3),在得到用位移表示的應(yīng)力分量后,將其代入運(yùn)動(dòng)方程(1),通過(guò)整理后可以得到以位移為未知函數(shù)的運(yùn)動(dòng)方程:

式中:u(r,t)是位移矢量,λ(r)與μ(r)是拉梅常數(shù),ρ(r)是介質(zhì)的質(zhì)量密度,它們與縱波波速及橫波波速之間的關(guān)系為

縱波波速:橫波波速:

在研究彈性波在聲子晶體內(nèi)產(chǎn)生的帶隙作用時(shí),平面波展開(kāi)法是常用的方法[14]。其基本思想是:由于聲子晶體具有結(jié)構(gòu)周期性,故而可以將彈性常數(shù)、密度等參量按照傅里葉級(jí)數(shù)展開(kāi),并與Bloch定理相結(jié)合,將彈性波波動(dòng)方程在倒格矢空間內(nèi)以平面波疊加的方式展開(kāi),即將波動(dòng)方程轉(zhuǎn)換為本征值求解,從而得到聲子晶體的能帶結(jié)構(gòu)。由于結(jié)構(gòu)的周期性,拉梅參數(shù)λ、μ和密度ρ都是空間r=(x,y)的周期性函數(shù),從而將各材料參數(shù)進(jìn)行傅里葉級(jí)數(shù)展開(kāi):(此處使用函數(shù)g統(tǒng)一來(lái)表示各材料參量)通過(guò)相關(guān)推導(dǎo),可將傅里葉展開(kāi)系數(shù)g(G)表示為

式中,G為倒格矢,對(duì)于二維晶格,倒格矢定義為

P(G)為散射體的幾何參數(shù),其定義為

由于散射體形狀多種多樣,產(chǎn)生的結(jié)構(gòu)函數(shù)也互不相同,其顯現(xiàn)出的能帶特性也有許多差別,本文重點(diǎn)討論正方形晶格條件下幾種具有代表性的散射體能帶特性,即:圓形、橢圓形、正方形及正六邊形。

圖1 不同形狀散射體示意圖Fig.1 Different shapes of scatters

對(duì)于圖1所示的圓形、橢圓形、正方形、正六邊形散射體,其形狀不同主要體現(xiàn)在結(jié)構(gòu)函數(shù)P(G)的差異,采用不同的散射體形狀,其結(jié)構(gòu)函數(shù)可表示為[15]

1)圓形:

式中:J1(·)為一階第一類Bessel函數(shù)。

2)橢圓形:

式中:Gx'=Gxaov/bov,Gy'=Gy,aov為長(zhǎng)半軸的長(zhǎng)度,bov為短半軸的長(zhǎng)度。

3)正方形:

4)正六邊形:

如前所述,在實(shí)際中由于加工誤差及環(huán)境因素的存在,散射體的排列方式可能并不是整齊劃一,而是存在一定旋轉(zhuǎn)角,經(jīng)轉(zhuǎn)化可以得出以下關(guān)系式[16]:

位移場(chǎng)由于結(jié)構(gòu)的周期性,故而可以采用Bloch定理進(jìn)行分解,產(chǎn)生結(jié)果為

式中:k為限制在第一Brillouin區(qū)內(nèi)的Bloch波矢,uk(r)是與各材料參數(shù)具有相同周期的函數(shù),同樣將其展成傅里葉級(jí)數(shù)的形式[17-18]:

將以上理論推導(dǎo)的結(jié)果采用特征矩陣的形式,可以表示為[17]

式中:Q和R分別為剛度矩陣和質(zhì)量矩陣。通過(guò)求解這個(gè)標(biāo)準(zhǔn)的矩陣特征問(wèn)題,并繪制特征頻率Ω=ωa/2πCt與波數(shù)關(guān)系曲線,即可得到聲子晶體的帶隙結(jié)構(gòu)。

2 數(shù)值結(jié)果與討論

本節(jié)采用MATLAB對(duì)上述散射體旋轉(zhuǎn)的二維聲子晶體帶隙結(jié)構(gòu)進(jìn)行仿真計(jì)算。以鋁圓柱體按照正方形晶格周期排列在環(huán)氧樹(shù)脂基體中構(gòu)成的二維聲子晶體作為分析算例,相關(guān)模型參數(shù)如表1所示。

取晶格常數(shù)a=0.01 m,為保證精度,取平面波數(shù)為625個(gè),為方便問(wèn)題的說(shuō)明,本文采用帶隙寬度ΔΩ=Ωmax-Ωmin及帶隙中心頻率Ωm=ΔΩ/2描述及對(duì)比不同形狀散射體在正方形晶格內(nèi)的帶隙特性。

表1 仿真模型材料參數(shù)Table 1 Material parameters of simulation model

2.1 圓形散射體

首先,對(duì)圓形散射體采用表1中所示的相關(guān)數(shù)據(jù),應(yīng)用所建立的模型得到描述聲子晶體帶隙結(jié)構(gòu)的特征值方程,求解特征值方程,得到如圖2所示的能帶結(jié)構(gòu)圖。

圖2 圓形散射體二維聲子晶體能帶結(jié)構(gòu)Fig.2 Band structure of 2-D phononic crystals with circle scatterers

本文所得到帶隙曲線與文獻(xiàn)[14]中對(duì)同樣模型參數(shù)的計(jì)算結(jié)果基本一致,從而驗(yàn)證了本文理論建模和程序編制的正確性。根據(jù)圖2可以明顯地看出,經(jīng)過(guò)測(cè)算,XY模態(tài)下最大正則化禁帶寬度ΔΩmax為0.472 9,Z模態(tài)下最大正則化禁帶寬度ΔΩmax為0.375 8。由于旋轉(zhuǎn)角度對(duì)于圓形散射體不產(chǎn)生影響,故而不予探討。

2.2 正方形散射體

對(duì)于正方形散射體,因存在結(jié)構(gòu)的對(duì)稱性,只用針對(duì)旋轉(zhuǎn)角度θ由0向45°變化過(guò)程進(jìn)行分析。首先對(duì)倒格子矢量依照式(14)進(jìn)行轉(zhuǎn)換,然后選取表1中的材料參數(shù),依照所建立的理論模型形成特征值矩陣并求解特征值方程,得到聲子晶體帶隙結(jié)構(gòu)圖,記錄禁帶寬度及出現(xiàn)位置。

由圖3可以看出,伴隨著填充率的提高,禁帶寬度呈現(xiàn)出先擴(kuò)大后縮小的趨勢(shì),且在填充率很小時(shí)(f=0.1時(shí))產(chǎn)生禁帶寬度極小,且在該填充率下,當(dāng)旋轉(zhuǎn)角度超過(guò)25°時(shí),禁帶已經(jīng)消失。在填充率小于或等于0.4時(shí),伴隨著旋轉(zhuǎn)角度的增大,禁帶寬度縮小,但是速度逐步減緩。當(dāng)填充率大于0.4時(shí),伴隨旋轉(zhuǎn)角度的增加,禁帶寬度下降明顯,且填充率越大,下降速率越快。

圖3 不同散射體填充率條件下,正方形散射體禁帶寬度與偏轉(zhuǎn)角度關(guān)系圖Fig.3 Relationship between rotation angle and gap width of square scatterers under various filling fractions

圖4 不同散射體填充率條件下,正方形散射體禁帶中間頻率與偏轉(zhuǎn)角度關(guān)系圖Fig.4 Relationship between rotation angle and midgap of hexagon scatterers under various filling fractions

由圖4可以看出,伴隨著填充率的提高,禁帶中間頻率呈現(xiàn)逐步上升的趨勢(shì),即產(chǎn)生最大禁帶的位置逐漸上升。當(dāng)填充率等于0.1,結(jié)合圖3可知,在前幾階特征頻率下并沒(méi)有禁帶產(chǎn)生,在第8、9階特征頻率下有極窄禁帶產(chǎn)生,且當(dāng)旋轉(zhuǎn)角度超過(guò)25°時(shí),禁帶基本消失。

圖5為當(dāng)填充率f=0.4時(shí),旋轉(zhuǎn)角θ分別為0及45°時(shí)Z模態(tài)下帶隙結(jié)構(gòu)圖。

圖5 正方形散射體Z模態(tài)禁帶結(jié)構(gòu)圖Fig.5 Z mode of band gap of rectangular scatterers

2.3 正六邊形散射體

由于正六邊形散射體的對(duì)稱性,取旋轉(zhuǎn)角度θ范圍在0~30°變化,分析過(guò)程同2.2節(jié),得到聲子晶體帶隙結(jié)構(gòu)圖,記錄禁帶寬度及位置。如圖6及圖7所示禁帶寬度及禁帶中間頻率都會(huì)產(chǎn)生顯著變化。

圖6 不同散射體填充率條件下,正六邊形散射體禁帶寬度與偏轉(zhuǎn)角度關(guān)系圖Fig.6 Relationship between gap width of hexagon scatterers and rotation angle under various filling fractions

如圖6所示,與正方形散射體的情況類似,伴隨著填充率f的提高,禁帶寬度總體呈現(xiàn)出先擴(kuò)大再縮小的變化趨勢(shì)。此外,當(dāng)填充率過(guò)低或填充率過(guò)高時(shí),禁帶寬度較小,禁帶寬度的變化也較小。伴隨著填充率的上升,禁帶寬度的變化也更加劇烈,當(dāng)f>0.3時(shí),禁帶寬度均呈現(xiàn)出先擴(kuò)大后縮小的趨勢(shì),且當(dāng)旋轉(zhuǎn)角θ=15°時(shí),禁帶寬度均達(dá)到最大值。

如圖7所示,伴隨著填充率f的提高,禁帶中間頻率總體呈現(xiàn)上升趨勢(shì),與正方形的情況相類似,伴隨著旋轉(zhuǎn)角度的增加,禁帶中間頻率較為固定,即禁帶出現(xiàn)的位置較為固定。在實(shí)驗(yàn)過(guò)程中,也發(fā)現(xiàn),無(wú)論填充率取到怎樣的值,在旋轉(zhuǎn)角θ=15°時(shí),禁帶出現(xiàn)的位置均有下降。

圖7 不同散射體填充率條件下,正六邊形散射體禁帶中間頻率與偏轉(zhuǎn)角度關(guān)系圖Fig.7 Relationship between midgap of hexagon scatterers and rotation angle under various filling fractions

圖8 為當(dāng)填充率f=0.5時(shí),旋轉(zhuǎn)角θ分別為0及15°時(shí)Z模態(tài)下帶隙結(jié)構(gòu)圖。

圖8 正六邊形散射體Z模態(tài)禁帶結(jié)構(gòu)圖Fig.8 Z mode of band gap of hexagon scatterers

2.4 橢圓形散射體

取橢圓短半軸與長(zhǎng)半軸長(zhǎng)度之比為aov∶bov=0.8∶1。分析過(guò)程如2.2所述,得到聲子晶體帶隙結(jié)構(gòu)圖,并記錄禁帶寬度及出現(xiàn)位置。

如圖9所示,伴隨著填充率f的提升,禁帶寬度ΔΩ總體呈現(xiàn)逐步擴(kuò)大的趨勢(shì)。當(dāng)填充率f較大時(shí)(f>0.5),伴隨著旋轉(zhuǎn)角度θ的上升,禁帶寬度有著明顯上升的趨勢(shì)。但是當(dāng)填充率f<0.5時(shí),伴隨著旋轉(zhuǎn)角度的上升,禁帶寬度呈現(xiàn)略微下降的趨勢(shì)。

關(guān)于禁帶產(chǎn)生的位置,如圖10所示,伴隨著填充率f的提升,禁帶中間頻率Ωm逐漸上升,當(dāng)填充率f保持不變時(shí),禁帶中間頻率穩(wěn)定,不存在明顯的變化。旋轉(zhuǎn)角對(duì)于禁帶中間頻率(即禁帶產(chǎn)生的位置)不存在明顯影響。

圖9 不同散射體填充率條件下,橢圓形散射體禁帶寬度與散射體偏轉(zhuǎn)角度關(guān)系圖Fig.9 Relationship between gap width of elliptical scatterers and rotation angle under various filling fractions

圖10 不同散射體填充率條件下,橢圓形散射體禁帶中間頻率與散射體偏轉(zhuǎn)角關(guān)系圖Fig.10 Relationship between gap width of elliptical scatterers and rotation angle under various filling fractions

圖11 為出現(xiàn)最大禁帶寬度f(wàn)=0.7,θ分別取0及45°時(shí),其Z模態(tài)禁帶結(jié)構(gòu)圖。

圖11 橢圓形散射體Z模態(tài)禁帶結(jié)構(gòu)圖Fig.11 Z mode of band gap of elliptical scatterers

3 結(jié)論

本文得到如下結(jié)論:1)散射體的旋轉(zhuǎn)角對(duì)于禁帶寬度具有顯著影響。根據(jù)仿真結(jié)果可以看出,當(dāng)散射體形狀取為正方形時(shí),在填充率固定的情況下,伴隨著旋轉(zhuǎn)角度的上升,禁帶寬度總體呈現(xiàn)縮小趨勢(shì)。但是當(dāng)散射體形狀取為正六邊形時(shí),同樣在填充率固定的情況下,伴隨著旋轉(zhuǎn)角度的上升,禁帶寬度總體呈現(xiàn)先擴(kuò)大再縮小的趨勢(shì)。2)散射體旋轉(zhuǎn)角對(duì)于正六邊形及正方形散射體的影響明顯大于對(duì)于橢圓形散射體的影響。3)無(wú)論散射體取為何種形狀,旋轉(zhuǎn)角度對(duì)于聲子晶體產(chǎn)生禁帶的位置(禁帶中間頻率)影響并不明顯。但較為明顯的是,伴隨著填充率的上升,最大禁帶出現(xiàn)的位置都呈現(xiàn)上升趨勢(shì)。

本文研究成果對(duì)于實(shí)際情況具有較好的指導(dǎo)意義,在實(shí)際應(yīng)用中可以通過(guò)機(jī)械裝置對(duì)散射體進(jìn)行人為偏轉(zhuǎn),從而達(dá)到人為控制聲子晶體帶隙的目的。

[1]張榮英,姜根山,王璋奇,等.聲子晶體的研究進(jìn)展及應(yīng)用前景[J].聲學(xué)技術(shù),2006,25(1):35-42.ZHANG Rongying,JIANG Genshan,WANG Zhangqi,et al.Progess of phononic crystal and the application perspectives[J].Technical Acoutics,2006,25(1):35-42.

[2]溫激鴻,韓小云,王剛,等.聲子晶體研究概述[J].功能材料,2003,34(4):1-3. WEN Jihong,HAN Xiaoyun,WANG Gang,et al.Review of phononic crystals research[J].Journal of Functional Materials,2003,34(4):1-3.

[3]SIGALAS M M,ECONOMOU E N.Elastic and acoustic wave band structure[J].Journal of Sound and Vibration,1992,158(2):377-382.

[4]KUSHWAHA M S,HALEVI P,ROUHANI B D.Acoustic band structure of periodic elastic composite[J].Phys Rev Lett,1993,71(1):2022-2025.

[5]MARTINEZ R S,SANCHO J,SANCHEZ J V,et al.Sound attenuation by sculpture[J].Nature,1995,378(1):241-248.

[6]KUSHWAHA M S,HALEVI P,MARTíNEZ G.Theory of acoustic band structure of periodic elastic composites[J].Phys Rev B,1994,49(4):2313-2321.

[7]MESEGURE F,HOLGADO M.Rayleigh-wave attenation by a semi-infinite two-dimensional elastic-band-gap crystal[J].Physical Review B,1999,59(19):12169-12172.

[8]VASSEUR J O,DRYMIER P A,CHENNI B,et al.Experimental and theoretical evidence for the existence of absolute acoustic band gaps in two-dimensional solid phononic crystals[J].Physical Review Letters,2001,86(1):3021-3024.

[9]BERTODI K,BOYCE M C.Mechanically triggered transformation of phononic band gap in periodic elastomeric structures[J],Physical Review B,2008(77):052105.

[10]MULIN T,DESCHANEL S,BERTODI K,et al.Pattern transformation triggered by deformation[J].Phsical.Review Letter,2007(99):084301.

[11]BERTODI K,BOYCE M C.Wave propagation and instabilities in monolithic and periodically structured elastomeric materials undergoing large deformation[J].Phsical Review,2008(78):184107.

[12]YAN L,ZHAO H P,WANG X Y.Influence of scatterer's tropism on the band gaps of two dimension phononic crystal[J].Technical Acoustics,2006,25(6):608-612.

[13]WU F G,LIU Z Y,LIU Y Y.Acoustic band gaps created by rotating square rods in a two dimensional lattice[J].Phys Rev E,2002,66:046628.

[14]溫熙森,溫激鴻.聲子晶體[M].北京:國(guó)防工業(yè)出版社,2009:61-68.

[15]王毅澤.磁電彈性聲帶隙材料中的彈性波傳播與局部化研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2009:18-32.WANG Yize.Elastic wave propagation and localization in magnetoelectroelastic band gap materials and structures[D].Harbin:Harbin Institute of Technology,2009:18-32.

[16]WANG R Z,WANG X H,GU B Y,et al.Effects of shapes and orientations of scatterers and lattice symmetries on the photonic band gap in two-dimensional photonic crystals[J].J Appl Phys,2001,90(9):4307-4313.

[17]趙言誠(chéng).二維聲子晶體結(jié)構(gòu)設(shè)計(jì)及其特性研究[D].哈爾濱:哈爾濱工程大學(xué),2006:15-34.ZHAO Yancheng.Two dimensional phononic crystal design and its characteristics investigations[D].Harbin:Harbin Engineering University,2006:15-34.

Effects of the rotation angle of scatterers on the band gap structure of a two-dimensional phononic crystal

DU Jingtao,HE Yanbo,F(xiàn)ENG Haocheng
(College of Power and Energy Engineering,Harbin Engineering University,Harbin 150001,China)

In order to investigate the influence of rotation of scatterers on the propagation characteristic of elastic waves in phononic crystals,the plane wave expansion method is utilized in this paper.The factor of rotation angle between scatterers and lattice axis is introduced into the structure function to explore the change of position as well as the width of the band gap structure for square,ellipse and hexagon scatterers corresponding to the change of rotation angle and filling fractions.The results demonstrated that the rotation angle of scatterers has significant impact on the band gap structure of two dimensional phononic crystals.In the meantime,such impact on the hexagon and square scatterers is more obvious than the elliptical scatterers.However,no matter what the shape of the scatterer is,the influence of rotation angle on the band gap position(mid-gap)of phononic crystals is not obvious.

phononic cyrstals;bandgap;forbidden band;rotation angles;plane wave expansion method;scatterers

10.3969/j.issn.1006-7043.201309069

http://www.cnki.net/kcms/doi/10.3969/j.issn.1006-7043.201309069.html

O481.1

A

1006-7043(2014)11-1358-06

2013-09-22.網(wǎng)絡(luò)出版時(shí)間:2014-09-22.

國(guó)家自然科學(xué)基金資助項(xiàng)目(11202056).

杜敬濤(1981-),男,副教授,博士.

杜敬濤,E-mail:dujingtao@hrbeu.edu.cn.

主站蜘蛛池模板: 99国产在线视频| 国产午夜福利亚洲第一| 亚洲视频免费播放| 国产毛片不卡| 99在线视频免费观看| 国产男女免费视频| 99ri国产在线| 亚洲黄网在线| 国产69精品久久| 四虎永久在线精品国产免费 | 五月天久久综合国产一区二区| 日本色综合网| 国产一级无码不卡视频| 亚洲av片在线免费观看| 色综合久久无码网| 天天色综合4| 一本久道久综合久久鬼色| 日韩精品无码免费专网站| 亚洲国产精品美女| 久热99这里只有精品视频6| 久草视频精品| 日本人又色又爽的视频| 国产簧片免费在线播放| 人妻无码AⅤ中文字| 91青青视频| 精品国产免费观看一区| 欧美a在线视频| 亚洲人成网线在线播放va| 国产福利不卡视频| 亚洲精品成人片在线播放| 国产xx在线观看| 免费国产在线精品一区| 91福利在线观看视频| 永久毛片在线播| 国产成人精品18| 日韩不卡高清视频| 青青草国产免费国产| 国产综合日韩另类一区二区| 伊人AV天堂| 国产十八禁在线观看免费| 亚洲无码高清视频在线观看| 91成人免费观看| 国产理论精品| 国产无码精品在线| 国产精品成人一区二区不卡| 日韩亚洲综合在线| 精品久久人人爽人人玩人人妻| h视频在线播放| 日韩国产欧美精品在线| 美女裸体18禁网站| 日本91在线| 国产午夜一级毛片| 99久久国产自偷自偷免费一区| 在线亚洲小视频| 中文字幕乱码中文乱码51精品| 亚洲经典在线中文字幕| 日韩福利在线观看| 真实国产乱子伦高清| 亚洲综合色婷婷| 九九这里只有精品视频| 欧美国产在线看| 亚洲A∨无码精品午夜在线观看| 国产精品 欧美激情 在线播放| 91精品情国产情侣高潮对白蜜| 亚洲欧美在线看片AI| 中日无码在线观看| 亚洲精品无码人妻无码| 毛片视频网| 99热最新在线| 欧美一级视频免费| 青青草国产一区二区三区| 亚洲无限乱码| 99久久性生片| 久久综合九色综合97婷婷| 幺女国产一级毛片| www.亚洲色图.com| 超薄丝袜足j国产在线视频| 国产精品毛片一区视频播| 91成人免费观看| 国产精品欧美激情| 亚洲成aⅴ人在线观看| 亚洲综合狠狠|