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

一種識(shí)別聲源噪聲輻射區(qū)域的方法

2017-05-04 05:49:27蘇俊博朱海潮蘇常偉
船舶力學(xué) 2017年2期
關(guān)鍵詞:模態(tài)區(qū)域

蘇俊博,朱海潮,蘇常偉

(海軍工程大學(xué) 船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室,武漢 430033)

一種識(shí)別聲源噪聲輻射區(qū)域的方法

蘇俊博,朱海潮,蘇常偉

(海軍工程大學(xué) 船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室,武漢 430033)

針對(duì)大型結(jié)構(gòu)體的聲輻射問題,提出了一種利用聲輻射模態(tài)識(shí)別聲源表面的噪聲輻射區(qū)域的方法。分析了各階聲輻射模態(tài)對(duì)聲源的輻射聲功率的貢獻(xiàn)量,找出了對(duì)遠(yuǎn)場(chǎng)輻射聲功率貢獻(xiàn)最大的幾個(gè)主要的聲輻射模態(tài),然后利用這幾個(gè)主要的聲輻射模態(tài)重建聲源表面法向振速,通過聲源表面法向振速的重建結(jié)果實(shí)現(xiàn)了聲源表面噪聲輻射區(qū)域的識(shí)別。通過對(duì)平板聲源在幾種不同頻率下的噪聲輻射區(qū)域的仿真分析驗(yàn)證了文中方法的正確性。該文方法對(duì)于確定特定頻率下聲源表面的噪聲輻射區(qū)域,從而進(jìn)一步進(jìn)行輻射噪聲控制具有積極的意義。

噪聲輻射區(qū)域;聲輻射模態(tài);輻射聲功率;聲場(chǎng)重建;減振降噪

0 引 言

為降低飛機(jī)、船舶等大型結(jié)構(gòu)體的噪聲,準(zhǔn)確識(shí)別結(jié)構(gòu)體表面的噪聲輻射區(qū)域具有很重要的意義。近年來近場(chǎng)聲全息技術(shù)在噪聲源的識(shí)別與定位方面取得了很大的成功[1-2],然而降低結(jié)構(gòu)體的輻射噪聲與減小噪聲源的振動(dòng)之間并沒有必然的直接聯(lián)系,減振而不降噪的情況經(jīng)常出現(xiàn)。這是因?yàn)槁曉幢砻嬲駝?dòng)在波數(shù)域中可分為超聲速波數(shù)成分和亞聲速波數(shù)成分[2],而僅有部分的波數(shù)成分能夠?qū)h(yuǎn)場(chǎng)的輻射噪聲產(chǎn)生影響,以無限大平板為例,對(duì)遠(yuǎn)場(chǎng)輻射有貢獻(xiàn)的僅為波數(shù)域中的超聲速波數(shù)成分。因此聲源表面振動(dòng)最大的區(qū)域的噪聲輻射能力并不一定最強(qiáng)。對(duì)于降低遠(yuǎn)場(chǎng)輻射噪聲而言,準(zhǔn)確識(shí)別聲源結(jié)構(gòu)表面的噪聲輻射區(qū)域比噪聲源的識(shí)別與定位更有意義。Williams[3-4]利用超聲速聲強(qiáng)來識(shí)別聲源結(jié)構(gòu)表面的噪聲輻射區(qū)域,Magelhaes和Tenenbaum[5]將超聲速聲強(qiáng)技術(shù)擴(kuò)展到任意結(jié)構(gòu)的聲源。由于超聲速聲強(qiáng)技術(shù)是基于輻射圓內(nèi)部的波數(shù)成分積分得到,而對(duì)有限體積的結(jié)構(gòu),輻射圓外部的一些波數(shù)成分仍然能夠向遠(yuǎn)場(chǎng)輻射噪聲,因此超聲強(qiáng)技術(shù)的計(jì)算結(jié)果并不是足夠準(zhǔn)確。聲輻射模態(tài)理論[6-7]的提出為解決這一問題提供了新的途徑。Marburg[8]利用聲輻射模態(tài)理論構(gòu)造了一個(gè)聲源表面輻射聲功率貢獻(xiàn)因子,用來評(píng)價(jià)聲源結(jié)構(gòu)表面各部分對(duì)遠(yuǎn)場(chǎng)輻射噪聲的貢獻(xiàn)。本文基于聲輻射模態(tài)理論,提出了一種不同于Marburg計(jì)算表面聲功率輻射貢獻(xiàn)因子的方法,首先找出了對(duì)遠(yuǎn)場(chǎng)輻射聲功率貢獻(xiàn)最大的幾個(gè)主要的聲輻射模態(tài),然后利用這幾個(gè)主要的聲輻射模態(tài)重建聲源表面法向振速,通過聲源表面法向振速的重建結(jié)果即可實(shí)現(xiàn)聲源表面噪聲輻射區(qū)域的識(shí)別。該方法相對(duì)于Marburg方法的優(yōu)點(diǎn)在于實(shí)現(xiàn)噪聲輻射區(qū)域識(shí)別的同時(shí),找出了對(duì)遠(yuǎn)場(chǎng)輻射聲功率貢獻(xiàn)最大的幾個(gè)主要的聲輻射模態(tài),有利于進(jìn)一步利用聲輻射模態(tài)理論對(duì)聲源的輻射噪聲進(jìn)行控制。

1 聲源噪聲輻射區(qū)域的識(shí)別

將聲源表面離散成一系列等面積的輻射單元,其輻射聲功率可表示為[9]

其中:R=(s/ 2) Re[Z],Z為輻射聲阻抗,s為輻射單元面積,v為聲源表面速度向量。將R進(jìn)行特征值分解:

其中:Φ為N×N維矩陣,Φ的列向量φi( i=1,2,…N )即為聲輻射模態(tài)向量。Λ=diag[λ1,λ2,…,λi,…]是以特征值為對(duì)角線元素的對(duì)角矩陣,特征值在對(duì)角線上按由大到小的順序依次排列。由于特征值與對(duì)應(yīng)的聲輻射模態(tài)的輻射效率成正比,因此聲輻射模態(tài)向量在聲輻射模態(tài)矩陣Φ中按照輻射效率的大小由低到高依次排列。為了更好地說明聲輻射模態(tài)的這個(gè)性質(zhì),選擇一正方形簡(jiǎn)支鋼板,對(duì)其聲輻射模態(tài)在波數(shù)域的波數(shù)譜進(jìn)行分析。假設(shè)正方形鋼板邊長(zhǎng)為0.5m,分析頻率為380 Hz,則平板前十階聲輻射模態(tài)對(duì)應(yīng)的波數(shù)譜如圖1所示。

圖1 平板前十階聲輻射模態(tài)的波數(shù)譜Fig.1 The wave number spectrum of the front ten acoustic radiationmodes

圖1中聲輻射模態(tài)波數(shù)譜圖中的圓圈代表輻射圓,從圖中可以看出,在此分析頻率下,輻射圓內(nèi)部包含主要波數(shù)成分的聲輻射模態(tài)主要是前六階,并且輻射圓內(nèi)包含的主要波數(shù)成分隨著階數(shù)的增大逐漸減少。這解釋了聲輻射模態(tài)隨著階數(shù)的增大輻射效率逐漸降低的原因。但是僅以輻射圓內(nèi)部是否包含主要波數(shù)成分來判斷該階聲輻射模態(tài)是否對(duì)遠(yuǎn)場(chǎng)輻射有影響不夠準(zhǔn)確,因?yàn)楦鶕?jù)邊角模態(tài)輻射理論,對(duì)于有限體積的聲源,輻射圓外部的少數(shù)波數(shù)成分仍然能向遠(yuǎn)場(chǎng)輻射能量。本文中以輻射聲功率來確定向遠(yuǎn)場(chǎng)輻射噪聲的主要聲輻射模態(tài),具體方法將在下節(jié)論述。假定向遠(yuǎn)場(chǎng)輻射噪聲的主要的聲輻射模態(tài)有M階,利用這M階聲輻射模態(tài)重建聲源表面的法向振速分布,則重建的法向振速分布就是向遠(yuǎn)場(chǎng)輻射噪聲的振動(dòng)形式,聲源表面重建法向振速幅值較大的區(qū)域代表了向遠(yuǎn)場(chǎng)輻射噪聲的區(qū)域。具體重建過程如下。

聲源結(jié)構(gòu)表面的法向速度向量可用聲輻射模態(tài)表示成如下形式

其中:c為輻射模態(tài)展開系數(shù)。假設(shè)將輻射體表面離散成N個(gè)輻射單元,并計(jì)算得到了N階聲輻射模態(tài)。由于(3)式具有良好的收斂性[10],用少量的聲輻射模態(tài)就可近似地表示出振速向量,因此表面速度可用截?cái)嗟穆曒椛淠B(tài)表示為

其中:矩陣Φ(N×N1)表示只取了前N1階聲輻射模態(tài)。當(dāng)輻射體表面速度向量中有N2個(gè)元素是已知的,可以得到由這N2個(gè)元素組成的方程組

為保證(5)式的適定性,輻射體表面速度向量中已知元素的數(shù)量大于聲輻射模態(tài)的數(shù)量即可,而不要求知道輻射體表面全部的法向速度。由(5)式可以求得聲輻射模態(tài)系數(shù)

其中:(Φ′)+(N1×N2)是Φ′(N2×N1)的偽逆。

利用M階主要的聲輻射模態(tài)可由(7)式對(duì)聲源表面法向振速中向遠(yuǎn)場(chǎng)輻射噪聲的部分

其中:Φ0為由M階主要的聲輻射模態(tài)組成的矩陣,c0為這M階聲輻射模態(tài)對(duì)應(yīng)的聲輻射模態(tài)展開系數(shù)。

2 仿真驗(yàn)證

本節(jié)通過對(duì)一帶有無限大障板的平板的數(shù)值仿真驗(yàn)證本文方法的有效性,并在本節(jié)中提出了尋找對(duì)遠(yuǎn)場(chǎng)輻射噪聲有貢獻(xiàn)的聲輻射模態(tài)的方法。本節(jié)中選取一個(gè)受單點(diǎn)簡(jiǎn)諧激勵(lì)的四周為無限大障板的簡(jiǎn)支鋼質(zhì)平板為研究對(duì)象。平板的長(zhǎng)度、寬度和厚度為L(zhǎng)x×Ly×h=0.5m×0.5m× 0.008 m,楊氏模量E取2×1011Pa,υ為泊松比取0.28,材料密度ρ0=7 800 kg/m3,空氣中聲速c0=343 m/s。用公式(8)計(jì)算該簡(jiǎn)支板的模態(tài)頻率,可知該簡(jiǎn)支板的(1,1)、(2,2)、(3,3)階模態(tài)頻率分別為153 Hz,612 Hz,1 378 Hz。

圖2 1 378 Hz激勵(lì)下平板表面法向振速幅值Fig.2 The amplitude of velocity of the plate at the exciting frequency of 1 378 Hz

當(dāng)以頻率為1 378 Hz、幅值為10 N的簡(jiǎn)諧力激勵(lì)平板,平板將以(3,3)階模態(tài)振動(dòng)。將平板表面均勻分為32×32個(gè)面積單元,其表面法向振速幅值如圖2所示。本節(jié)將對(duì)以(3,3)階模態(tài)振動(dòng)的該簡(jiǎn)支平板在不同分析頻率下的噪聲輻射區(qū)域進(jìn)行分析。

2.1 確定輻射聲功率的主要聲輻射模態(tài)

能否找出對(duì)遠(yuǎn)場(chǎng)輻射噪聲有主要貢獻(xiàn)的聲輻射模態(tài),對(duì)聲源表面的噪聲輻射區(qū)域的識(shí)別結(jié)果的準(zhǔn)確與否有著至關(guān)重要的影響。聲輻射模態(tài)對(duì)遠(yuǎn)場(chǎng)輻射噪聲的貢獻(xiàn)可以用該階聲輻射模態(tài)的輻射聲功率表示,而聲源的輻射聲功率可以表示成如下的形式[7]

其中:λi為第i階聲輻射模態(tài)對(duì)應(yīng)的特征值,ci為聲輻射模態(tài)展開系數(shù)。用前M階聲輻射模態(tài)的輻射聲功率來估計(jì)聲源的輻射聲功率,則聲源的輻射聲功率的估計(jì)誤差可由下式表示:

對(duì)于上述由簡(jiǎn)諧力激勵(lì)的簡(jiǎn)支平板而言,分析頻率為300 Hz時(shí),該簡(jiǎn)支平板的輻射聲功率的估計(jì)誤差與聲輻射模態(tài)的截止階數(shù)之間的關(guān)系如圖3所示。從圖3中輻射聲功率的估計(jì)誤差與聲輻射模態(tài)的截止階數(shù)之間的關(guān)系曲線中可以看出,輻射聲功率的估計(jì)誤差隨聲輻射模態(tài)截止階數(shù)的增大迅速減小,因此只有前面少量的聲輻射模態(tài)對(duì)遠(yuǎn)場(chǎng)的輻射噪聲產(chǎn)生影響。更需要引起重視的是,在圖3中的局部放大圖中可以看到在2-5階,6-10階,輻射聲功率的估計(jì)誤差并沒有隨著截止階數(shù)的增大而降低,這說明在平板此種振動(dòng)形式下3-5階、7-10階聲輻射模態(tài)對(duì)輻射聲功率基本沒有貢獻(xiàn),不是輻射聲功率的主要聲輻射模態(tài),這些模態(tài)的存在對(duì)聲源表面噪聲輻射區(qū)域的識(shí)別造成干擾,所以在重建聲源表面的法向振速時(shí)要將這些非主要聲輻射模態(tài)去除。第i階聲輻射模態(tài)的輻射聲功率可由下式計(jì)算得到:

圖3 輻射聲功率估計(jì)誤差Fig.3 The estimation error of radiated sound power

圖4 輻射聲功率估計(jì)誤差(重新排序)Fig.4 The estimation error of radiated sound power (rearranged acoustic radiationmodes)

將聲輻射模態(tài)按輻射聲功率由大到小重新進(jìn)行排序,并再次計(jì)算輻射聲功率的估計(jì)誤差,得到圖4。從圖4中可以看出,經(jīng)過重新排序后,輻射聲功率估計(jì)誤差在聲輻射模態(tài)截止階數(shù)為2時(shí)就下降到了1.5‰,因此輻射聲功率的主要模態(tài)為前兩階,此時(shí)的前兩階聲輻射模態(tài)即為未重新排序以前的第一階和第六階聲輻射模態(tài),利用此兩階聲輻射模態(tài)重建平板表面法向振速,即可實(shí)現(xiàn)對(duì)平板表面噪聲輻射區(qū)域的識(shí)別。

表1 各分析頻率下平板的主要聲輻射模態(tài)Tab.1 The dom inant acoustic radiation modes at the analyzing frequencies

當(dāng)分析頻率為500 Hz、700 Hz、900 Hz、1 100 Hz、1 300 Hz和1 500 Hz時(shí),經(jīng)分析得到的輻射聲功率的主要聲輻射模態(tài)如表1所示。

2.2 平板聲源的噪聲輻射區(qū)域分析

由2.1節(jié)得到的各分析頻率下向遠(yuǎn)場(chǎng)輻射噪聲的主要聲輻射模態(tài),利用(7)式可完成對(duì)平板表面法向振速的重建,從而利用此重建結(jié)果可實(shí)現(xiàn)各分析頻率下平板表面噪聲輻射區(qū)域的識(shí)別。為驗(yàn)證本文方法的準(zhǔn)確性,利用Marburg的方法做了同樣的分析,并將本文得到的結(jié)果與Marburg的方法得到的結(jié)果作比對(duì),如圖5、圖6所示。圖5、圖6中第一列是平板表面法向振速幅值,第二列是Marburg的方法得到的平板表面各區(qū)域?qū)椛渎暪β实呢暙I(xiàn)因子,第三列是本文方法重建得到的法向振速的幅值,圖5中由上到下每一行分別是分析頻率為300 Hz、500 Hz、700 Hz和900 Hz時(shí)的仿真結(jié)果,圖6中由上到下每一行分別是分析頻率為1 100 Hz、1 300 Hz和1 500 Hz時(shí)的仿真結(jié)果。

從圖5、圖6中可以觀察到當(dāng)分析頻率在300-900 Hz時(shí),平板表面的噪聲輻射區(qū)域集中在四個(gè)角上,隨著分析頻率的繼續(xù)增大,平板表面噪聲輻射區(qū)域也越來越大,這樣的分析結(jié)果與邊角輻射模態(tài)理論一致。此外,通過對(duì)本文提出的方法得到的識(shí)別結(jié)果與Marburg的方法得到的結(jié)果進(jìn)行比較可知二者基本一致,驗(yàn)證了本文方法對(duì)噪聲輻射區(qū)域識(shí)別結(jié)果的準(zhǔn)確性。由于本文中的方法在識(shí)別噪聲輻射區(qū)域的過程中,首先找出了向遠(yuǎn)場(chǎng)輻射聲功率的主要聲輻射模態(tài),這一結(jié)果不僅有利于聲源表面噪聲輻射區(qū)域的識(shí)別,而且對(duì)利用聲輻射模態(tài)進(jìn)行聲輻射主動(dòng)控制具有很重要的意義,可有效減少控制的通道數(shù)。

3 結(jié) 論

本文提出了一種識(shí)別聲源表面噪聲輻射區(qū)域的方法,該方法利用對(duì)遠(yuǎn)場(chǎng)輻射聲功率有貢獻(xiàn)的聲輻射模態(tài)重建聲源表面振動(dòng),識(shí)別出了聲源結(jié)構(gòu)表面向遠(yuǎn)場(chǎng)輻射噪聲的區(qū)域。前面的研究表明:(1)各階聲輻射模態(tài)之間能夠?qū)崿F(xiàn)輻射聲功率解耦,聲源的輻射聲功率能夠用各階聲輻射模態(tài)的輻射聲功率的線性疊加表示,并且存在對(duì)遠(yuǎn)場(chǎng)輻射聲功率貢獻(xiàn)較大的幾個(gè)主要的聲輻射模態(tài),利用這幾個(gè)主要的聲輻射模態(tài)重建聲源表面法向振速即可實(shí)現(xiàn)對(duì)聲源噪聲輻射區(qū)域的識(shí)別;(2)聲源結(jié)構(gòu)表面振動(dòng)最大的區(qū)域并不一定是對(duì)遠(yuǎn)場(chǎng)輻射聲功率貢獻(xiàn)最大的區(qū)域,當(dāng)選定的分析頻率遠(yuǎn)小于結(jié)構(gòu)聲源的激勵(lì)頻率時(shí),聲源表面的邊角部分對(duì)輻射聲功率貢獻(xiàn)最大。本文中是以平板聲源為研究對(duì)象,但本文的研究方法和結(jié)論同樣適用于具有其他幾何外形的聲源。

[1]Maynard JD,Williams E G,Lee Y.Near field acoustic holography I.Theory of generalized holography and the developmentof NAH[J].JAcoust.Soc.Am.,1985,78(4):1395-1413.

[2]Williams Earl G.Fourier acoustics:Sound radiation and nearfield acoustical holography[M].London:Academic Press, 1999:28-35.

[3]Williams E.Supersonic acoustic intensity[J].JAcoust.Soc.Am.,1995,97(1):121-127.

[4]Williams E.Supersonic acoustic intensity on planar sources[J].JAcoust.Soc.Am.,1998,104(5):2845-2850.

[5]Magelhaes M B S,Tenenbaum R B.Supersonic acoustic intensityfor arbitrarily shaped sources[J].Acta Acust.Acust., 2006,92:189-201.

[6]Borgiotti G V.The power radiated by a vibrating body in an acoustic fluid and its determination from boundary measurements[J].JAcoust.Soc.Am.,1990,88(4):1884-1893.

[7]Cunefare K A.The design sensitivity and control of acoustic power radiated by three-dimensional structures[D].Ph.D.thesis,The Pennsylvania State University,1990.

[8]Steffen Marburg.Surface contributions to radiated sound power[J].JAcoust.Soc.Am.,2013,133(6):3700-3705.

[9]Elliott S J,Johnson M E.Radiation modes and the active control of sound power[J].JAcoust.Soc.Am.,1993,94(4): 2194-2204.

[10]聶永發(fā),朱海潮.利用源強(qiáng)密度聲輻射模態(tài)重建聲場(chǎng)[J].物理學(xué)報(bào),2014,63(10):104303. Nie YongFa,Zhu Haichao.Acoustic field reconstruction using source strength density acoustic radiation modes[J].Acta Phys.Sin.,2014,63(10):104303.

A method to identify the noise radiation regions on the surface of acoustic source

SU Jun-bo,ZHU Hai-chao,SU Chang-wei
(National Key Laboratory on Ship Vibration&Noise,Naval University of Engineering,Wuhan 430033,China)

For the acoustic radiation problems of large scale vibrating structure,amethod based on acoustic radiationmodes is proposed to identify the noise radiation regions on the structure.To identify the noise radiation regions on the surface of acoustic source,the contributions of acoustic radiationmodes to radiated sound power are analyzed,and themain acoustic radiationmodes,which have themost contribution to the radiated sound power,are selected.By the selected acoustic radiationmodes,the velocities on the surface are reconstructed.The identification of noise radiation regions is realized by the reconstructed resultof the velocities.The validity of the proposedmethod is validated by numerical simulating for the noise radiation regions of a rectangular plate at several frequencies.The proposed method has a positive significance for identifying the noise radiation regions of acoustic source and further reducing the radiation noise in the farfield.

noise radiation regions;acoustic radiationmodes;radiated sound power; acoustic field reconstruction;vibration and noise reduction

TB532 TB533

:Adoi:10.3969/j.issn.1007-7294.2017.02.014

2016-06-14

國家自然科學(xué)基金項(xiàng)目(51305452)

蘇俊博(1985-),男,博士研究生,E-mail:sujunbo0618@163.com;

朱海潮(1963-),男,教授,博士生導(dǎo)師。

1007-7294(2017)02-0237-07

猜你喜歡
模態(tài)區(qū)域
永久基本農(nóng)田集中區(qū)域“禁廢”
分割區(qū)域
車輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對(duì)比
關(guān)于四色猜想
分區(qū)域
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
高速顫振模型設(shè)計(jì)中顫振主要模態(tài)的判斷
基于嚴(yán)重區(qū)域的多PCC點(diǎn)暫降頻次估計(jì)
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
由單個(gè)模態(tài)構(gòu)造對(duì)稱簡(jiǎn)支梁的抗彎剛度
主站蜘蛛池模板: 99精品一区二区免费视频| 国产精品亚洲精品爽爽| 免费不卡视频| 国产麻豆永久视频| 无码专区在线观看| 国产日本欧美在线观看| 久青草免费视频| 国产欧美视频综合二区| 天天躁日日躁狠狠躁中文字幕| 午夜人性色福利无码视频在线观看| 波多野结衣无码视频在线观看| 国产在线视频欧美亚综合| 亚洲Aⅴ无码专区在线观看q| 狠狠亚洲五月天| 亚洲AV无码乱码在线观看代蜜桃 | 毛片三级在线观看| 日韩免费成人| 国产高清在线丝袜精品一区| 久久中文字幕不卡一二区| 欧美伦理一区| 国模视频一区二区| 久久99热这里只有精品免费看 | 成人日韩精品| 久久91精品牛牛| 国产成人精品一区二区三在线观看| 91青青草视频在线观看的| 三区在线视频| 欧美成人一级| 日本在线免费网站| 伊人久久大香线蕉aⅴ色| 欧美全免费aaaaaa特黄在线| jizz在线免费播放| 一本二本三本不卡无码| 亚洲综合在线最大成人| 伊人久久婷婷五月综合97色 | 日韩A∨精品日韩精品无码| 亚洲国产精品不卡在线 | 久久鸭综合久久国产| 亚洲第一视频网站| 国产精品久久久久久久久久久久| 欧美亚洲国产一区| 91精品国产自产在线观看| 久久精品66| 亚洲精品另类| 亚洲第一色视频| 91免费国产高清观看| av天堂最新版在线| 9啪在线视频| 91成人免费观看| 东京热高清无码精品| 国产成人乱码一区二区三区在线| 91小视频在线| 国产精品成人不卡在线观看| 在线免费亚洲无码视频| 欧美精品三级在线| 久久综合色视频| 欧美亚洲国产精品第一页| 制服丝袜亚洲| 久久精品一品道久久精品| 日韩美毛片| 在线观看网站国产| 丝袜久久剧情精品国产| 国产欧美日韩另类精彩视频| 91蝌蚪视频在线观看| 国产精品亚欧美一区二区三区| 亚洲日韩精品综合在线一区二区| 亚洲成人一区二区三区| 无码内射在线| 四虎精品免费久久| 野花国产精品入口| 黄色片中文字幕| 久久性视频| 在线看片中文字幕| 国产精品亚洲精品爽爽| 亚洲中文字幕23页在线| 少妇人妻无码首页| 亚洲男人的天堂久久香蕉网| 四虎永久免费地址在线网站| 国产真实自在自线免费精品| 日韩二区三区无| 国产精品久久久久久久久久98| 国产农村妇女精品一二区|