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

周期矩形輪廓聲學散射體的散射性質預測算法及其應用研究

2017-04-21 00:50:36王海濤曾向陽杜博凱劉延善
振動與沖擊 2017年7期
關鍵詞:性質

王海濤, 曾向陽, 杜博凱, 劉延善

(西北工業大學 航海學院,西安 710072)

周期矩形輪廓聲學散射體的散射性質預測算法及其應用研究

王海濤, 曾向陽, 杜博凱, 劉延善

(西北工業大學 航海學院,西安 710072)

具有周期排列形式的散射體是建筑聲學中一種常見的聲學擴散結構,針對各種室內環境中最為常見的周期矩形輪廓散射體,基于光柵方程發展了一種散射性質預測算法,并利用此方法對矩形輪廓散射體的散射性質進行了預測及分析。先詳細介紹預測算法的推導過程,此算法僅利用散射體的幾何參數及入射波參數即可計算各個階次反射波的能量,從而對散射性質進行預測。通過對不同算例散射系數的計算及對比證明了預測算法具有良好的正確性及計算效率,并利用預測算法對周期矩形輪廓散射體的散射性質進行分析。分析表明其具有缺級散射、對稱散射性質,這些性質不但可以使周期矩形輪廓散射體在室內聲學擴散問題中得到更加合理的應用,還使其有潛力在其他領域,如噪聲控制、空間濾波等問題中發揮重要的作用。

周期散射體;矩形輪廓;散射系數;缺級散射;對稱散

聲散射是室內聲學中的一項重要研究內容。BERANEK[1]提出了室內聲學中響度、混響感、親切感、溫暖感和環繞感五個音質因素,并提出相應的客觀音質參數,其中每項參數都與聲場散射程度密切相關。在1996年BERANEK[2]對音質評價指標的補充中,與散射相關的音質指標同樣占據大多數。因此,為了增加聲場的擴散程度從而提高音質,需要按使用需求在室內安裝聲學散射體[3]。具有周期排列形式的矩形輪廓散射體是一種最為常見的界面聲散射結構,基于其設計制造簡單、美觀等特點,這種結構多被建筑設計人員應用于音樂廳、錄音室、劇院等對聲場散射要求較高的場所[4]。類似于微觀尺度下周期結構中振動波傳播的帶隙特性[5],周期形式的散射體也具有不同于一般聲散射結構的獨特性質,它既可以通過增加散射、抑制回聲以提高室內音質水平[6],也有可能因為散射指向性導致頻率響應發生變化,引起聲染色效應,對室內音質產生負面影響[7]。因此,對周期矩形輪廓散射體散射性質進行準確預測及估計對于改善室內音質具有十分重要的意義。

定義為非鏡面反射聲能與全部反射聲能比值的散射系數是描述結構聲學散射能力的基本參數,其計算方法一直是一個研究熱點。在1945年,RAYLEIGH[8]就已經通過平面波分解的方法對正弦型輪廓的周期散射體進行了散射系數計算方面的研究。此后,不斷有解析方法被提出,WIRGIN等[9-10]對此進行了詳細的總結。不過,由于實際中的散射體形狀通常多變,這些要求散射體形狀規則的解析方法實用性并不強。隨著有限元邊界元等數值計算方法的興起,在散射系數計算領域又出現了Boundary Element Method (BEM)[11],Finite Element-Plane Wave Decomposition (FE-PWD)[12-13]兩種方法。BEM結合了邊界元法與MOMMERTZ[14]提出的散射系數計算理論,適用于尺寸有限大的周期結構,而FE-PWD最早出現于吸聲系數的數值計算之中,結合了有限元法與平面波分解理論,適用于尺寸無限大的周期結構。另外,隨著新型數值計算方法無網格法在聲場數值仿真中的應用[15],以網格法為基礎的散射系數計算方法也獲得研究[16]。上述方法雖然有效,但由于以離散模型為基礎,存在著計算量過大、計算效率低下、上限計算頻率低等問題。因此有必要發展快速、準確的散射系數估計方法。

除了散射系數計算外,周期形式散射體對入射聲波的指向性散射也是一項研究熱點。與光學領域中的光柵類似,周期散射體可以將入射波按照不同的階次反射向不同的方向,這種指向性散射在某些情況下會引起室內空間中的聲場分布發生變化,導致原始聲音被賦予外加的音色特點,使聽者產生主觀聽感上的厭惡情緒,嚴重影響聽音效果[17]。因此,正確地預測周期散射體指向性對于降低散射體負面影響并提高室內音質水平十分重要。

本文針對室內環境中最為常見的周期矩形輪廓散射體,基于光柵方程所表示的周期結構柵格效應發展了一種散射性質快速預測算法,可根據散射體的幾何參數求得各個階次反射波的能量,此方法可對周期矩形輪廓散射體的散射系數、指向性散射性質進行快速計算及預測。論文首先詳細介紹了快速預測算法的推導過程;然后利用無規入射散射系數及方向入射散射系數的對比計算分析了快速計算方法的正確性與適用性;最后利用快速計算方法對周期矩形輪廓散射體的指向性散射性質進行了預測分析,分析表明此種散射體具有缺級散射、對稱散射的性質。

1 散射性質預測算法

假設一束具有單位幅值的平面波入射到具有周期矩形輪廓的散射體上(見圖1),散射體的周期長度為L,凹槽長度為l,高度為H,入射波的聲壓可表示為

pi(x,y)=e-jk(sin θix-cos θiy)

(1)

式中:θi為入射角;k=2π/λ 為波數;λ為波長。與光柵類似,周期型的聲學散射體同樣具有可用光柵方程所表示的柵格效應[18-19],在二維平面內,散射波的階次及方向可根據光柵方程求得

(2)

式中:θs為散射角;λ為入射波的波;n為散射波的階次;n=0即代表鏡面反射波。

圖1 周期矩形輪廓散射體及入射波示意圖Fig. 1 Schematic of the periodic diffuser with rectangular profile and the incident wave

為了計算反射波的能量,現對階次為n、散射角為θs的反射波進行分析。由于矩形輪廓的散射體具有凹槽結構,上述反射波應來源于兩部分:①被散射體頂面所反射的入射波,此處假設其傳播路徑為Path1;②來源于進入凹槽內部,經過多重反射之后的入射波,其傳播路徑為Path2。以圖1中所示的散射體頂面為參考線,沿兩個路徑的聲波經過反射后所傳播的距離分別為

(3)

則兩個路徑上的聲波在參考線處的相位差為

(4)

假設散射體表面為剛性結構,則忽略能量在反射面處的損失,根據相位差可得到周期散射體的反射系數為

(5)

由于式(5)為周期函數,可將其表示為傅里葉級數形式

(6)

式中,αn為傅里葉系數,可表示為

(7)

其中的反射系數表示反射聲壓與入射聲壓之比,根據式(7),可得到在參考線處的反射聲壓為

(8)

從式(8)的形式可以看出,反射波的聲壓由一系列聲壓幅值為αn的聲波構成,因此,只要求得αn,即可得到對應n階反射波的能量。通過聯立式(5)及式(7),求得αn的值為

(9)

式中,τ=l/L,表示凹槽長度與周期長度的比值。最終,經過推導可得到各階次反射波的能量

(10)

式(10)給出了各階次反射波的能量,通過式(10)及式(2)所示的光柵方程,即可對周期矩形輪廓散射體的各項散射性質進行預測。

2 算例驗證

為了驗證本文所發展預測算法的正確性,論文利用此方法對周期矩形輪廓散射體的散射系數進行了計算,并與實驗測量結果及經典邊界元法結果進行了比較。由于本文所發展的算法最終給出的是各階次反射波的能量,因此需要利用下式來計算散射系數

(11)

式(11)給出了某個入射方向上的方向散射系數,在求解無規入射散射系數時,需要以二維散射體底邊中心為球心,在散射體上方布置一個半球面,在半球面上均勻取點代表不同入射方向的聲源。當入射方向不在二維散射體所在的xoy平面內時,可通過如圖2所示的頻率偏移方式將其轉化到平面內,即

Sd(θi′,φ′,k)=Sd(θi′,0,kcosφ′)

(12)

在求得所有入射方向的散射系數之后,可按照Paris方程求得無規入射散射系數

(13)

圖2 入射波的入射方向不在xoy平面內時的轉換方式Fig. 2 Converting of the incident wave lies outside the xoy plane

為了全面驗證本文方法的正確性,此處對兩個周期矩形輪廓散射體的無規入射散射系數進行了計算,第一個算例散射體的周期長度及凹槽長度為L=0.06 m,l=0.03 m,高度為H=0.03 m。此處將計算結果與在國際標準:ISO 17497-1下的實驗測量結果[20]進行了比較,對比如圖3所示。

圖3 本文方法與實驗測量結果之間無規入射散射系數的對比Fig. 3 Comparisons of the random-incidence scattering coefficients between the method in this paper and measurement

從圖3可知,利用本文方法所得到的無規入射散射系數與實驗測量結果非常接近,兩者之間的差值在所有頻率上均<0.1,所有頻率上的平均差值為0.020,此算例可初步證明本文方法的正確性。

第二個算例中,周期矩形輪廓散射體的尺寸為L=0.2 m,l=0.1 m,H=0. 02 m,此處計算了散射體的無規入射散射系數及方向入射散射系數,并將計算結果與正確性已得到證明的3D BEM法的計算結果[21]進行了對比,對比如圖4、圖5所示。圖5中的灰度表示不同入射方向上的散射系數值,圓的徑向表示入射方向的俯仰角,圓周方向表示入射方向的方位角。

圖4 本文方法與BEM法之間無規入射散射系數的對比Fig. 4 Comparisons of the random-incidence scattering coefficients between the method in this paper and BEM

圖5 本文方法與BEM法之間方向散射系數的對比Fig. 5 Comparisons of the directional scattering coefficients between the method in this paper and BEM

從圖4所示的無規入射散射系數計算結果對比來看,與第一個算例一樣,兩者之間的差值均<0.1;在所有頻率上,兩者平均差值為0.037 6,體現了本文方法與傳統方法的結果十分接近。

從圖5所示的方向入射散射系數結果來看,首先兩種方法在圖案清晰度上有著明顯區別,BEM方法在圖案清晰度上較為模糊,而本文方法則比較清晰,這由兩種方法不同理論所導致:本文方法在計算反射波能量時,首先會對散射方向進行預判斷,對于某些入射方向有可能僅存在鏡面反射波,因此其散射系數會直接判定為0;而BEM在計算時會對每一個入射方向均進行計算,從而導致其計算結果均>0。從方向入射散射系數的分布來看,兩種方法結果十分接近,在主要散射范圍及散射系數數值方面均具有較高的一致性,當頻率為100 0 Hz時,只有較小范圍內的入射波會發生明顯的散射現象,隨著頻率的升高,發生散射的入射波的方向范圍也逐漸增大。

上述兩個算例證明本文所發展的散射性質預測算法具有良好的正確性,可對周期矩形輪廓散射體的無規入射散射系數及方向入射散射系數進行準確計算。

3 周期矩形輪廓散射體散射性質

本文利用前述算法對周期矩形輪廓散射體的散射性質進行了分析,分析表明,與其他形式的散射體相比,矩形輪廓的散射體具有一些與其幾何因素密切相關的散射性質,分別包括缺級散射性質、振蕩散射性質以及對稱散射性質,以下為具體分析。

3.1 缺級散射

本文所發展的算法給出了周期矩形輪廓散射體對入射波的反射方向及各階次反射波的能量,從式(10)所示的反射波能量表達式可以預測這種形式的散射體當尺寸具有特殊比例時會表現出與光柵類似的缺級散射性質。式(10)中:En為鏡面反射波以外其他階次反射波的能量,其中n為反射波階次,τ為凹槽長度與周期長度的比值。當n與τ共同決定cos 2nπτ的值為1時,會使此階反射波能量為0,從而出現缺級散射現象。例如,當τ=1/2,即凹槽長度為周期長度一半時,偶數階次的反射波能量會為0(n=±2,±4,…),從而發生偶數階次缺級;同理,當τ=1/3時,3的倍數階次n=±3,±6,…的反射波會發生缺級現象。另外,根據反射聲波能量表達式可知,對應于n=±1的反射波不會發生缺級散射現象,因為τ取值范圍在0~1,在此情況下n=±1會使cos 2nπτ取值位于-1~1,因此En不會為0。

為了對比證明周期矩形輪廓的缺級散射性質,本文利用無網格法[22]實現了BEM邊界積分方程理論,并對兩個三維周期矩形輪廓散射體的散射聲場進行了計算,其中一個散射體L=0.2 m,l=0.1 m,滿足τ=1/2的關系,另一個散射體L= 0.2 m,l=0.13 m,兩個散射體的高度均為H=0.04 m,且周期數均為15。為了清楚的顯示反射波的方向,接收點設置在與圖2中一個相同的半圓面上,兩個示例性頻率上的計算結果如圖6所示。

根據式(2),散射體周期長度為0.2 m,在320 0 Hz入射角為45°的情況下,應有四束反射波,分別為45°(n=0),10°(n=-1),-21°(n=-2)以及-62°(n=-3)。從圖6可知,凹槽長度為l=0.1 m的散射體其n=-2階的反射波消失,而l=0.13 m的散射體則具有所有階次的反射波,這證明了由本文方法所預測的周期矩形輪廓散射體的缺級散射性質。4 400 Hz的情況同樣證明了此結論,在垂直入射的情況下,對于具有倍數比例關系的第一個散射體,其n=±2階次的反射波消失,而第二個散射體則存在所有階次的反射波。

另外,圖6所示的反射波的能量也與本文算法的結果一致。例如,對于第一個散射體,在4 400 Hz時按照式(10)的計算,其歸一化鏡面反射波能量應為0.994 7,占據了絕大部分的反射能量,而圖4顯示了相同的情況,鏡面反射波的能量要遠遠大于其他階次反射波的能量,這也進一步證明了本文所發展的算法具有良好的正確性。

圖6 表示缺級散射的兩種不同尺寸散射體的聲場Fig. 6 Sound fields of two diffusers with different dimensions to illustrate the missing order scattering

3.2 對稱散射

式(10)所示的反射波能量表達式顯示周期矩形輪廓散射體還具有對稱散射性質。E0表達式中的2τ(1-τ)是一個關于τ=0.5的對稱函數,而En表達式中的cos 2nπτ同樣是一個關于τ=0.5的對稱函數。因此,當兩個散射體的凹槽長度與周期長度之比關于τ=0.5對稱時,它們會具有相同的散射系數。為了證明此預測,此處對多個散射體的無規入射散射系數進行了計算,這些散射體的周期長度及高度均為L=0.2 m,H=0.04 m,凹槽長度分別為l=0.1 m,0.2 m,0.3 m,…,0.9 m。計算結果如圖7所示。

圖7結果顯示,τ=0.1與τ=0.9的散射系數曲線完全重合,τ=0.2與τ=0.8、τ=0.3與τ=0.7、τ=0.4與τ=0.6的情況同樣如此,證明周期矩形輪廓散射體具有對稱散射性質。從散射系數的結果來看,當τ=0.5時,散射體具有最大的散射系數,隨著τ值偏離0.5,散射系數數值有減小的趨勢。另外,從圖7所示的無規入射散射系數計算結果來看,周期矩形輪廓散射體表現出了明顯的振蕩散射性質,即散射系數的曲線具有類周期性的變化規律,且這種周期性隨著頻率的升高逐漸減弱,散射系數最終收斂于某一固定值。這是由于反射波能量E0與En的表達式中均含有參數cosΔφ。cosΔφ是一個入射波波長相關的周期性函數,當入射波波長發生變化時,cosΔφ所發生的周期性變化也會導致各階次反射波的能量隨之發生類似的周期性變化。這說明矩形輪廓散射體的散射系數對頻率較為敏感,較小的頻率變化也有可能對散射系數產生明顯的影響。

圖7 周期矩形輪廓散射體對稱散射示意圖Fig. 7 Symmetrical scattering of the periodic diffuser with rectangular profile

4 結 論

針對室內空間中常見的周期矩形輪廓散射體,本文從研究方法及散射性質兩方面展開了研究。論文首先基于光柵方程發展了一種散射性質預測算法,此算法僅利用二維散射體的幾何參數及入射波的相關參數即可計算經過散射體反射的聲波能量,而無需傳統方法中耗時的數值仿真,論文通過對不同尺寸散射體散射系數的計算及對比證明了算法具有良好的正確性。論文還利用所推導的算法對周期矩形輪廓散射體的散射性質進行了分析,分析表明矩形輪廓散射體具有多種特殊的散射性質:

(1)缺級散射,當散射體的凹槽長度與周期長度符合一定比例關系時會導致某些階次的反射波消失。

(2)對稱散射,當兩個周期矩形輪廓散射體的凹槽長度與周期長度之比τ關于τ=0.5對稱時,它們會具有相同的散射系數,并且當一個散射體的凹槽長度與周期長度之比為0.5時,散射系數具有最大值。

基于這些性質的分析,周期矩形輪廓散射體不但可以在室內聲學中的獲得更加合理的應用,還有潛力在噪聲控制、空間濾波等其他領域發揮重要的作用。

[ 1 ] BERANEK L L. Music, acoustics & architecture[M]. New York: Wiley, 1962.

[ 2 ] BERANEK L L. How they sound: concert and opera halls[M]. New York: Acoustical Society of American, 1996.

[ 3 ] 王季卿. 聲場擴散與廳堂音質[J]. 聲學學報, 2001, 26(5): 417-421. WANG Jiqing. Sound diffusion and auditorium acoustics[J]. Acta Acustica, 2001, 26(5): 417-421.

[ 4 ] 樂意, 趙其昌, 沈勇, 等. 大型廳堂的建筑聲學設計方法研究[J]. 南京大學學報(自然科學版), 2011, 47(2): 208-217. LE Yi, ZHAO Qichang, SHEN Yong, et al. Study on the acoustic design of large auditoriums[J]. Journal of Nanjing University(Natural Sciences), 2011, 47(2): 208-217.

[ 5 ] 陳榮, 吳天行. 非對稱周期結構中耦合波的傳播特性[J]. 振動與沖擊, 2015, 34(1): 68-73. CHEN Rong, WU Tianxing. Coupled wave propagation in asymmetric periodic structures[J]. Journal of Vibration and Shock, 2015, 34(1): 68-73.

[ 6 ] 蔣國榮. 室內聲場模擬中的界面聲散射[J]. 聲學技術, 2009, 28(6): 697-700. JIANG Guorong. Sound scattering in room acoustic modeling[J]. Technical Acoustics, 2009, 28(6): 697- 700.

[ 7 ] 劉海生, 龔農斌. 室內聲學中散射研究進展[J]. 應用聲學, 2005, 24(2): 126-132. LIU Haisheng, GONG Nongbin. Progress in diffusers researches in room acoustics[J]. Applied Acoustics, 2005, 24(2): 126-132.

[ 8 ] RAYLEIGH J W. The theory of sound[M]. New York: Dover, 1945.

[ 9 ] WIRGIN A. Reflection from a corrugated surface[J]. Journal of the Acoustical Society of America, 1980, 68(2): 692-699.

[10] CHESNEAUX J M, WIRGIN A. Reflection from a corrugated surface revisited[J]. Journal of the Acoustical Society of America, 1994, 96(2): 1116-1129.

[11] KOSAKA Y, SAKUMAY T. Numerical examination on scattering coefficients of architectural surfaces using the boundary element method[J]. Acoustical Science & Technology, 2005, 26(2): 136-144.

[12] SAKAMOTO S, MUKAI H, TACHIBANA H. Numerical study on sound absorption characteristics of resonance-type brick/block walls[J]. Journal of the Acoustical Society of Japan, 2000, 21(1): 9-15.

[13] EMBRECHTS J J, GEETERE L D, VERMEIR G, et al. Calculation of the random incidence scattering coefficients of a sine-shaped surface[J]. Acustica United with Acustica, 2006, 92(4): 593-603.

[14] MOMMERTZ E. Determination of scattering coefficients from the reflection directivity of architectural surface[J]. Applied. Acoustics, 2000, 60(2): 201-203.

[15] 李鴻秋, 陳國平, 史寶軍. 多聯通封閉空間聲場響應的基于核重構的最小二乘無網格解法[J]. 振動與沖擊, 2012, 31(8): 148-163. LI Hongqiu, CHEN Guoping, SHI Baojun. Acoustic response in multi-domain based on least-square point collocation method and reproducing kernel particle method[J]. Journal of Vibration and Shock, 2012, 31(8): 148-163.

[16] 王海濤, 曾向陽. 周期結構聲散射系數的無網格數值計算方法[J]. 計算物理, 2013, 30(2): 229-236. WANG Haitao, ZENG Xiangyang. Meshless models for numerical calculation of scattering coefficients of periodic structures[J]. Chinese Journal of Computational Physics, 2013, 30(2): 229-236.

[17] CHOI Y J. Effects of periodic type diffusers on classroom acoustics[J]. Applied Acoustics, 2013, 74(5): 694-707.

[18] HOLFORD R L. Scattering of sound waves at the ocean surface: a diffraction theory[J]. Journal of the Acoust Society of America, 1981, 70(4): 1103-1115.

[19] HOLFORD R L. Scattering of sound waves at a periodic, pressure-release surface: an exact solution[J] Journal of the Acoust Society of America, 1981, 70(4): 1116-1128.

[20] CHOI Y J, JEONG D U, KIM J Y. Some issues in measurement of the random-incidence scattering coefficients in a reverberation room [J].Acta Acustica United with Acustica,2008,94(5):769-773.

[21] LEE H, SAKUMA T. Numerical characterization of acoustic scattering coefficients of one-dimensional periodic surfaces[J]. Applied Acoustics, 2015, 88: 129-136.

[22] 張雄, 劉巖. 無網格法[M]. 北京: 清華大學出版社, 2004.

Predicting scattering properties of a periodic-type diffuser with a rectangular profile

WANG Haitao, ZENG Xiangyang, DU Bokai, LIU Yanshan

(School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072,China)

The periodic-type diffuser is a kind of important acoustical structures in architectures. It is usually beneficial to improve diffusion and prevent echoes. However, sometimes it also leads to side effects of coloration because of its periodically arranged structure. Here, a grating effect-based method was developed to predict scattering properties of a periodic-type diffuser with a rectangular profile most widely used in room acoustics. Firstly, the derivation of the grating effect-based method was introduced in detail. This method was capable of calculating the scattering energy of reflection waves using the geometrical parameters of the diffuser without time-consuming numerical simulations. Then, the calculation and comparison of scattering coefficients in examples showed that the proposed method has a good accuracy to predict scattering properties. Lastly, the scattering properties of the periodic-type diffuser with a rectangular profile were analyzed using the grating effect-based method. The results demonstrated that the periodic-type diffuser with a rectangular profile has properties of missing order scattering and symmetric scattering; based on these special properties, the periodic-type diffuser with a rectangular profile can not only be more reasonably used in room acoustics, but also play an important role in other engineering fields, such as, spatial filtering, spectrum control, and noise control.

periodic type diffuser; rectangular profile; scattering coefficient; missing order scattering; symmetrical scattering

國家自然科學基金(11374241);中央高校基本科研業務費專項資金資助(3102015BJ(II)JJZ06)

2015-11-26 修改稿收到日期: 2016-02-14

王海濤 男,博士,講師,1986年生

曾向陽 男,博士,教授,1974年生

TU112

A

10.13465/j.cnki.jvs.2017.07.012

猜你喜歡
性質
含有絕對值的不等式的性質及其應用
MP弱Core逆的性質和應用
弱CM環的性質
一類非線性隨機微分方程的統計性質
數學雜志(2021年6期)2021-11-24 11:12:00
隨機變量的分布列性質的應用
一類多重循環群的剩余有限性質
完全平方數的性質及其應用
中等數學(2020年6期)2020-09-21 09:32:38
三角函數系性質的推廣及其在定積分中的應用
性質(H)及其攝動
九點圓的性質和應用
中等數學(2019年6期)2019-08-30 03:41:46
主站蜘蛛池模板: 亚洲一区黄色| 欧美国产视频| 久热re国产手机在线观看| 亚洲高清免费在线观看| 伦精品一区二区三区视频| 亚洲天堂福利视频| 91po国产在线精品免费观看| 欧美成人午夜影院| 无码福利视频| 久久毛片免费基地| 欧美日韩国产一级| 青草视频免费在线观看| 亚洲丝袜中文字幕| 亚洲一区二区三区国产精华液| 永久免费精品视频| av无码一区二区三区在线| 亚洲精品综合一二三区在线| 女人18毛片水真多国产| 精品无码一区二区在线观看| 欧美日韩国产在线播放| 精品综合久久久久久97| 亚洲福利网址| 国产福利大秀91| 日本免费精品| 自拍欧美亚洲| 久久人妻xunleige无码| 欧美一级专区免费大片| 国产人碰人摸人爱免费视频| 婷婷午夜天| 国产性爱网站| 婷五月综合| 999国内精品久久免费视频| 国产成人精品一区二区免费看京| 日韩中文欧美| 一本大道无码日韩精品影视| 在线免费无码视频| 国产91九色在线播放| 欧美一级夜夜爽| 美女亚洲一区| 亚洲人成电影在线播放| 婷婷99视频精品全部在线观看| 99久久国产自偷自偷免费一区| 日本尹人综合香蕉在线观看| 国产成人一区二区| 国产精品手机视频| 香蕉网久久| 亚洲男女天堂| 伊人欧美在线| 97国产精品视频自在拍| 国产精品久久久久鬼色| 国产免费高清无需播放器| 在线观看国产网址你懂的| 国产成人亚洲欧美激情| 午夜福利网址| 麻豆a级片| 亚洲男人在线| 日本免费新一区视频| 久久久久亚洲AV成人网站软件| 精品国产Av电影无码久久久| 国产91特黄特色A级毛片| 欧美一级高清免费a| 九色在线视频导航91| 一级香蕉视频在线观看| 欧美日本中文| 最新亚洲人成网站在线观看| 蜜臀AV在线播放| 99这里只有精品在线| 伊人激情久久综合中文字幕| 99视频在线看| 在线日本国产成人免费的| 亚洲国产精品不卡在线| 内射人妻无码色AV天堂| 亚洲中文字幕国产av| 国产一级在线观看www色 | 亚洲一区网站| 久久精品这里只有精99品| 久久综合亚洲色一区二区三区| 色综合久久88色综合天天提莫| 国产97区一区二区三区无码| 精品国产欧美精品v| 一级看片免费视频| 91麻豆国产精品91久久久|