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

橢球面鏡聲波反射聚焦數(shù)值模擬研究

2016-07-26 08:15:26陳正武曾新吾
振動與沖擊 2016年2期

張 軍, 陳 鵬, 陳正武, 趙 云, 曾新吾

(1.中國空氣動力研究與發(fā)展中心 氣動噪聲控制重點實驗室,四川 綿陽 621000;2.國防科學(xué)技術(shù)大學(xué) 光電科學(xué)與工程學(xué)院,長沙 410073)

?

橢球面鏡聲波反射聚焦數(shù)值模擬研究

張軍1, 陳鵬1, 陳正武1, 趙云2, 曾新吾2

(1.中國空氣動力研究與發(fā)展中心 氣動噪聲控制重點實驗室,四川綿陽621000;2.國防科學(xué)技術(shù)大學(xué) 光電科學(xué)與工程學(xué)院,長沙410073)

摘要:據(jù)Hamilton線性理論解,給出沿橢球面鏡軸線的電弧放電等離子聲源(AD-PSS)反射波計算結(jié)果,分析由鏡面引起的相位變化對中心波、邊緣波及尾波傳播影響。在聚焦前區(qū)中心波壓力為正,邊緣波、尾波壓力為負;聚焦后區(qū)則相反。在線性條件下,反射波壓力峰值出現(xiàn)于橢球面鏡幾何焦點,越過焦點后反射波壓力幅值迅速衰減。利用有限元軟件COMSOL對橢球面鏡聲反射過程進行數(shù)值模擬,揭示反射波傳播演化過程及聲場分布規(guī)律。據(jù)KZK(Khokhlov Zabolotskava Kuznetsov)方程及等效聲源法,分析非線性效應(yīng)對聲傳播過程影響。研究表明,非線性效應(yīng)將使橢球面鏡實際焦點位置偏離幾何焦點,即正壓實際焦點出現(xiàn)在幾何焦點后,負壓實際焦點出現(xiàn)在幾何焦點前;隨非線性系數(shù)增加正壓實際焦點后移,負壓實際焦點前移。

關(guān)鍵詞:等離子體聲源;反射聲波;橢球面鏡

我國大功率強聲技術(shù)發(fā)展較快,并廣泛應(yīng)用。如,醫(yī)學(xué)的體外沖擊波碎石、高強度聚焦超聲手術(shù),工業(yè)的強聲粉碎、無損檢測及遠距離廣播等。尤特金[1]從理論上證明高壓放電在流體介質(zhì)中會產(chǎn)生較強聲波,且為寬頻帶,大部分能量集中于可聽聲頻段。再如自然界中雷雨天氣時,在劇烈閃電同時伴隨的雷聲。

Hamilton等[2]對人工條件的強聲產(chǎn)生方式進行探索。Raymond等[3]在實驗室中成功制造出自然界雷聲。用儲能電容儲存高壓電,將電容與一對高強度銅鎢合金電極組成串聯(lián)電路,需產(chǎn)生強聲時使電路導(dǎo)通,立刻有上萬伏高壓電瞬間被加載至電極產(chǎn)生劇烈放電,電極間介質(zhì)隨之被擊穿生成高溫高壓等離子體,并迅速膨脹、擠壓周圍介質(zhì)產(chǎn)生較強聲波。由于該聲波可控(如重復(fù)頻率、脈沖寬度等),實際上形成新的高強度聲源。按聲波產(chǎn)生方式命名,該聲源被稱為等離子體聲源(Plasma Sound Source,PSS)。

金明劍[4]曾對PSS的電聲特性進行研究。黃逸帆等[5]成功研制出電暈放電(Corona Discharge)方式的等離子體聲源(CD-PSS),王一博等[6]成功研制出電弧放電(Arc Discharge)方式的等離子體聲源(AD-PSS)。而CD-PSS的單脈沖聲能量較低,但聲波產(chǎn)生時刻一致性較好。為實現(xiàn)聲能聚焦,可將放電電極組成陣列,并利用相控方法控制波束。AD-PSS的單脈沖聲能量較高,但聲波產(chǎn)生時刻具有較大隨機性。張軍等[7]針對AD-PSS提出用橢球面鏡反射方式實現(xiàn)聲能量聚焦。由于橢球面鏡幾何特性較特殊,即從第一焦點發(fā)出的球面波在橢球面反射后將在第二焦點處聚焦,從而使一部分原該擴散到其它方向的聲能量聚集到目標(biāo)方向。

本質(zhì)上講,橢球面鏡屬于曲面聲障板。O’Neil[8]最早利用線性聲學(xué)近似理論研究簡諧波經(jīng)凹圓面反射問題,導(dǎo)出沿軸線及焦平面的反射波壓力、質(zhì)點速度及聲強解析表達式,并發(fā)現(xiàn)凹圓面焦點處聲壓幅值與圓面半徑成正比,與聲波長成反比。Lucas等[9]從線性聲波方程出發(fā),通過Bessel-Fourier變換及拋物近似獲得沿凹圓面軸線及焦平面的反射聲壓表達式,提出諸如線性聚焦增益、焦區(qū)寬度及聚焦相位改變等重要概念。二者理論模型均針對反射面上具有均勻振速分布(如壓電晶體貼片)的聚焦聲源,結(jié)果不能直接用于點源的聲反射聚焦,僅給出頻域線性模型而無法給出反射聲場形成過程,且不能考慮聲傳播的非線性效應(yīng)。Hamilton[10]利用Kirchhoff聲衍射公式推導(dǎo)出沿橢球面鏡軸線的時域反射聲場理論解,認為反射聲波由不同形態(tài)波組成,并給出如方波、N波等的特殊波形傳播演化過程。

本文利用Hamilton時域理論解給出橢球面鏡反射聲場計算結(jié)果,分析鏡面相位調(diào)制對AD-PSS反射波傳播影響;利用多物理場有限元軟件COMSOL對AD-PSS反射波傳播過程進行數(shù)值模擬,驗證理論計算的3種反射波結(jié)構(gòu),揭示時域二維反射聲場演化規(guī)律;為分析實驗中觀察到的實際焦點位置偏移現(xiàn)象[11],利用KZK方程研究非線性傳播效應(yīng)的影響。

1理論模型

在AD-PSS聲波反射系統(tǒng)組成結(jié)構(gòu)中,放電電極位于橢球面鏡第一焦點Z1處。當(dāng)控制開關(guān)導(dǎo)通后,放電電極之間會產(chǎn)生一高強度球面波聲脈沖,據(jù)橢球面幾何性質(zhì),聲波經(jīng)橢球面反射后在第二焦點Z2處聚焦,見圖1。

圖1 橢球面鏡聲波反射聚焦示意圖Fig.1 Schematic of sound wave focusing using an ellipsoidal mirror

1.1線性理論

據(jù)衍射理論,封閉空間內(nèi)任意點聲壓可通過封閉面的聲壓確定。在AD-PSS反射系統(tǒng)中,認為封閉面由橢球面鏡與從鏡面開口到無窮遠處的假想半無限曲面構(gòu)成,通過Kirchhoff公式可獲得封閉面內(nèi)任意點反射聲壓表達式[12],即

(1)

式中:p為反射波聲壓;pg為按幾何聲學(xué)給出的封閉面聲壓邊界條件;c0為介質(zhì)中小擾動聲速;R為封閉面一點到觀察點位置矢量;n為垂直封閉面外法線矢量;s為封閉面。

為簡化,式(1)略去物理量時空變量。據(jù)索墨菲輻射條件可忽略從無窮遠處假想封閉面發(fā)出的反射波對聲場貢獻。則式(1)積分只在橢球面鏡S0上進行。式(1)中[·]表示物理量經(jīng)R/c0時間延遲,如f(t)=f(t-R/c0)。

(2)

式中:p0為初始聲壓幅值;f為聲波源函數(shù);a為橢球面鏡長半軸距離。

式(2)右端的3項分別稱為中心波、邊緣波及尾波。空間上此3種波分別對應(yīng)從橢球面鏡頂點、開口邊緣、頂點與開口邊緣間剩余表面發(fā)出的反射波。由觀察點與橢球面鏡的空間位置變化引起的相位改變會使3種波呈不同傳播特性。式(2)中待定求解變量為

(3)

式中:t1,t2選取為保證式(2)積分時間具有一致性,即積分上限大于下限。

在焦點Z2上,由于中心波、邊緣波到達時間相同(t1=t2),式(2)右端第3項尾波項積分為0,反射波聲壓表達式為

(4)

1.2非線性理論

各向均勻無旋流體介質(zhì)中非線性聲束傳播可用KZK方程描述,即

(5)

式(5)左端為聲波傳播項,右端各項為影響波傳播的各種效應(yīng),即衍射效應(yīng)、吸收效應(yīng)及非線性效應(yīng)。

為便于計算,引入坐標(biāo)變換將KZK方程無量綱化

σ=z/df,ρ=r/af,τ=ω0t′,P=p/p0

(6)

式中:df為橢球面鏡出口面到第二焦點距離;af為出口面圓半徑;ω0為聲波脈沖中心頻率。

KZK方程可無量綱化為

(7)

式中:G為線性聚焦增益;A為無量綱吸收系數(shù);N為無量綱非線性系數(shù);即

(8)

式(7)可用時域有限差分(FDTD)及算子分離方法(Operator Split Method)求解[13]。FDTD方法難以處理由橢圓特征曲線所致復(fù)雜邊界,為獲得聲波反射問題定解,利用等效聲源方法將邊界條件簡化。

將從橢球面鏡上發(fā)出的二次反射波等效為從其出口面發(fā)出的波,見圖2。區(qū)別在于出口面的反射波幅值、相位相對鏡面已發(fā)生改變。利用等效方法將不規(guī)則橢球面邊界轉(zhuǎn)化成簡單的圓形活塞聲源邊界,即

(9)

由階躍函數(shù)帶來的數(shù)值邊界會使反射波計算結(jié)果出現(xiàn)數(shù)值震蕩。為減小計算誤差,可在計算邊界處設(shè)置吸收條件或增大計算區(qū)域。

圖2 等效聲源方法示意圖Fig.2 Schematic of Equivalent source method

2數(shù)值模擬

算例的參數(shù)為:初始強聲脈沖波形為三角波,脈沖持時T=4 μs;橢球面鏡長半軸a=13.8 cm,短半軸b=7.75 cm,截斷后橢球面深度d=12.41 cm。傳播介質(zhì)為空氣時取小擾動聲速c0=340 m/s,密度ρ0=1.2 kg/m3。傳播介質(zhì)為水時取c0=1 500 m/s,密度ρ0=1 000 kg/m3。

2.1AD-PSS反射波沿橢球面鏡軸線傳播

取傳播介質(zhì)為空氣,在橢球面鏡軸線設(shè)置6個觀察點,即σ=linspace(0.4, 1.4, 6)(linspace為Matlab語言中線性分布函數(shù))。利用式(2)進行計算,可得各觀察點反射聲波結(jié)果。

在衍射效應(yīng)作用下AD-PSS反射波中出現(xiàn)中心波C,邊緣波E及尾波W等波形,分別對應(yīng)從橢球面鏡頂點、邊緣及頂點與邊緣間曲面發(fā)出的波,見圖3。

圖3 AD-PSS反射波沿橢球面鏡軸線傳播Fig.3 The evolution of AD-PSS reflection wave along thesymmetric axis of an ellipsoidal mirror

由圖3看出,① 時間上。在焦前區(qū)(σ<1)觀察點上3種波到達順序不同。中心波先到,尾波次之,邊緣波最后到。隨觀察點靠近焦點(σ=1),邊緣波相對中心波延遲時間縮短。時間軸上相當(dāng)于中心波向右傳播,邊緣波向左傳播。焦后區(qū)(σ>1),邊緣波越過中心波先到觀察點,因橢圓曲線幾何特點使3種波到觀察點的聲程不同。② 反射波相位。在焦前區(qū)中心波壓力為正,尾波、邊緣波壓力為負。表明除從橢球面鏡頂點發(fā)出的反射波外,初始壓力為正的聲波脈沖在鏡面其余部分反射后會經(jīng)歷180°相位改變成為負壓。在焦后區(qū)中心波滯后于尾波、邊緣波,并成為負壓的主要貢獻者。按射線理論聲波在剛性面反射時不存在相移。而由波動理論觀點,Kirchhoff公式由單極子源(時間偏導(dǎo)數(shù)項)及偶極子源(空間偏導(dǎo)數(shù)項)組成,反射波相位改變表明從鏡面發(fā)出的二次輻射波中偶極子源占主導(dǎo)地位。③ 反射聲波幅值。靠近焦點過程中,中心波、邊緣波及尾波幅值均逐漸增加,此為橢球面鏡幾何聚焦特性的體現(xiàn)。與拋物面鏡不同,橢球面鏡邊緣反射波幅值更低(在反射前球面波傳播路徑更長)。而在焦后區(qū),聲壓幅值更低的邊緣波(E)出現(xiàn)在聲壓幅值更高的尾波(W)前。實際上,在非線性效應(yīng)作用下,高壓區(qū)較低壓區(qū)傳播速度更快,邊緣波、尾波最終會疊加。

2.2聲波反射聚焦過程數(shù)值模擬

為更清楚認識AD-PSS反射波傳播過程,利用有限元軟件COMSOL進行數(shù)值模擬。在COMSOL中用參數(shù)化曲線建立橢球面鏡幾何模型,將鏡面設(shè)為剛性邊界。在橢球面鏡第一焦點設(shè)置高斯脈沖點源,計算域內(nèi)材料選為空氣,整個計算區(qū)域用三角形網(wǎng)格剖分。一個聲波長內(nèi)用6個計算網(wǎng)格,網(wǎng)格大小與時間步長之間滿足CFL條件。數(shù)值模擬結(jié)果見圖4。由圖4看出,從AD-PSS放電電極發(fā)出的聲波呈球面波形式擴散。一部分聲波前向傳播而形成直達波;另部分聲波受鏡面反射成反射波,反射波含3部分,即中心波、邊緣波及尾波(圖4(b))。邊緣波、尾波相位與中心波相反。在向橢球面鏡第二焦點傳播中反射波波陣面逐漸匯聚,聲波壓力峰值不斷升高;過第二焦點后反射波波陣面發(fā)散較快。因此,為獲得最高聲能量利用率,應(yīng)將目標(biāo)置于橢球面鏡第二焦點附近。

圖4 COMSOL軟件模擬AD-PSS反射波傳播過程Fig.4 The simulation results of the evolution of AD-PSS reflection wave using software COMSOL

2.3非線性效應(yīng)對反射波實際焦點位置影響

對AD-PSS強聲脈沖,反射波在橢球面鏡第二焦點壓力較高。適當(dāng)調(diào)整放電參數(shù)可使第二焦點附近反射波峰值壓力達幾十MPa。故應(yīng)考慮聲波傳播的非線性效應(yīng)。

圖5 KZK方程計算的AD-PSS反射波Fig.5 AD-PSS reflection waves given by KZK equation

由于空氣介質(zhì)非線性系數(shù)較小,為考察非線性效應(yīng)對聲波傳播影響,選傳播介質(zhì)為水,其非線性系數(shù)β= 3.5,吸收系數(shù)α=0.014 dB/m (MHz)。通過量綱歸一化關(guān)系,獲得線性聚焦增益G=24.5,吸收系數(shù)A=2.6E-3,非線性系數(shù)N=1.4。利用KZK方程進行數(shù)值計算,所得結(jié)果見圖5。其中實線表示KZK方程解,虛線表示Hamilton線性理論解。由圖5看出,焦前區(qū),通過KZK方程計算的邊緣波到達時間更晚(E2>E1),因KZK方程對聲衍射算子用拋物近似,其效果相當(dāng)于對邊緣波傳播的實際距離采取菲涅耳近似(即開方近似)。據(jù)非線性理論,聲波實際傳播速度為小擾動聲速c0與局部粒子速度u的疊加,故聲波高壓區(qū)傳播速度較低壓區(qū)快。隨傳播距離增加非線性效應(yīng)使中心波到達時間提前。焦后區(qū)非線性效應(yīng)使尾波能趕上邊緣波,從而避免線性條件下幅值低的正壓出現(xiàn)在幅值高的正壓之前。

據(jù)圖5計算結(jié)果,由線性理論所得反射波實際焦點位置與橢球面鏡幾何焦點位置一致。而在非線性條件下實際焦點位置相對幾何焦點位置發(fā)生偏離,符合實驗結(jié)果。沿橢球面鏡軸線,由線性、非線性理論計算所得反射波正壓峰值變化見圖6。由圖6看出,非線性效應(yīng)使焦點處反射波峰值正壓較線性理論解更高。據(jù)焦點處反射聲壓頻域解,焦點處壓力幅值與聲波頻率成正比,焦區(qū)-6 dB寬度與聲波頻率成反比[14]。在非線性條件下,聲波能量不斷從低頻泵到高頻,使焦點處聲壓峰值提高,但在總能量守恒前提下焦點處能量增加必以減弱其它位置能量為代價。反射波實際焦點位置相對幾何焦點偏離,可通過自折射理論解釋[15]。

圖6 反射波正壓峰值沿橢球面鏡軸線變化Fig.6 The peakpositive pressure along the axis of the mirror

圖7 橢球面鏡實際焦點位置偏移量隨非線性系數(shù)N變化Fig.7 The offsets between actual focus and geometric focus vs. nonlinear coefficients

由式(7)右端第3項看出,非線性效應(yīng)主要受非線性系數(shù)N控制。橢球面鏡焦距越長聲脈沖越短,初始聲壓幅值越高,N值越大。據(jù)KZK方程對實際焦點位置偏移量隨N值變化進行計算,結(jié)果見圖7。其中帶圓圈實線表示反射波正壓的實際焦點位置與幾何焦點位置差值(即實際焦點位置的偏移量);帶方塊實線表示反射波負壓實際焦點偏移量。在非線性效應(yīng)影響下負壓實際焦點出現(xiàn)在幾何焦點之前,正壓實際焦點出現(xiàn)在幾何焦點之后。故負壓實際焦點位置偏移量為負值,正壓實際焦點位置偏移量為正值。由圖7看出,隨非線性系數(shù)增加(如縮短聲脈沖寬度、提高聲脈沖壓力幅值或增大橢球面鏡焦距),負壓實際焦點向橢球面鏡焦前區(qū)移動,正壓實際焦點向焦后區(qū)移動。正壓或負壓實際焦點位置與幾何焦點位置的偏移量與非線性系數(shù)N近似呈線性關(guān)系。

3結(jié)論

本文針對橢球面鏡反射聚焦方式,通過數(shù)值模擬,分析反射波傳播及聲場分布規(guī)律,總結(jié)如下:

(1) 據(jù)Hamilton理論解,給出沿橢球面鏡軸線AD-PSS反射波計算結(jié)果,分析中心波、邊緣波及尾波傳播特性。焦前區(qū),中心波先于邊緣波、尾波到達,中心波壓力為正,邊緣波、尾波壓力為負。焦后區(qū)則相反。

(2) 線性條件下,反射波壓力峰值出現(xiàn)于橢球面鏡幾何焦點處。越過第二焦點后反射波壓力幅值迅速衰減。利用有限元軟件COMSOL對橢球面鏡聲反射聚焦過程進行數(shù)值模擬,揭示反射波傳播演化過程及二維聲場分布規(guī)律。

(3) 通過KZK方程計算結(jié)果,分析非線性效應(yīng)對聲傳播影響。非線性效應(yīng)使反射波到達時間及實際焦點位置發(fā)生改變,即正壓實際焦點出現(xiàn)在幾何焦點后,負壓實際焦點出現(xiàn)在幾何焦點前。隨非線性系數(shù)增加,正壓實際焦點位置后移,負壓實際焦點位置前移。本文研究結(jié)果對認識橢球面鏡聲反射傳播規(guī)律具有參考價值。

參 考 文 獻

[1] 尤特金. 液電效應(yīng)[M]. 北京:科學(xué)出版社, 1962.

[2] Hamilton M F, Tjotta J N, Tjotta S. Nonlinear effects in the farfield of a directive sound source[J].J. Acoust Soc. Am., 1985, 78(1): 202-216.

[3] Schaefer R, Zagaeski M,Yoshikawa S. High source level sparker for navy app lications[R]. Phase I SBIR, Final Report for Naval Air Warfare Center, Contract No. N00421-03-P-1081, 2003.

[4] 金明劍. 水下等離子體聲源電特性的基礎(chǔ)性研究[D].北京:中科院電工研究所, 2004.

[5] 黃逸帆. 等離子體聲源技術(shù)[C]//2011全國第一屆水下安保會議技術(shù)學(xué)術(shù)論文集,2011:198.

[6] 王一博,曾新吾. 水中沖擊波源壓強的理論估算[C]//2011全國第一屆水下安保會議技術(shù)學(xué)術(shù)論文集, 2011:194.

[7] 張軍,曾新吾. 強聲波脈沖在水下的自反射聚焦[J]. 聲學(xué)技術(shù), 2010, 29(6):114-115.

ZHANG Jun, ZENG Xin-wu. The self-focusing of intense sound wave in water[J]. Acoustic Technology, 2010,29(6):114-115.

[8] Oneil H T. Theory of focusing radiators [J]. J.Acoust. Soc. Am., 1949, 21(5): 516-526.

[9] Lucas B G, Muir T G. The field of a focusing source[J]. J. Acoust. Soc. Am., 1982, 73(4): 1289-1296.

[10] Hamilton M F. Transient axial solution for the reflection of a spherical wave from a concave ellipsoidal mirror[J]. J. Acoust. Soc. Am., 1993, 93(3): 1256-1269.

[11] 張軍,曾新吾,陳聃,等. 水下強聲波脈沖負壓的產(chǎn)生和空化氣泡運動[J]. 物理學(xué)報, 2012, 61(18): 184302.

ZHANG Jun, ZENG Xin-wu, CHEN Dan,et al. The generation of negative pressure of intensive acoustic pulse in water and cavitation bubble dynamics[J]. Acta Physica Sinca, 2012, 61(18): 184302.

[12] 季家镕. 高等光學(xué)教程[M]. 北京:科學(xué)出版社, 2007.

[13] Lee Y S,Hamilton M F. Time-domain modeling of pulsed finite-amplitude sound beams[J].J. Acoust Soc. Am., 1994, 97(2): 907-917.

[14] 王鴻樟,錢盛友. 球面脈沖波在凹橢球面上反射的聚焦聲場[J].上海交通大學(xué)學(xué)報, 1996, 30(5): 109-114.

WANG Hong-zhang, QIAN Sheng-you. The reflecto-focused sound field of spherical pulse wave from concave ellipsoidal surface[J]. Journal of Shanghai Jiaotong University, 1996, 30(5): 109-114.

[15] 張軍.等離子體聲源強聲反射聲場特性的理論和實驗研究[D]. 長沙:國防科技大學(xué), 2013.

基金項目:國家自然科學(xué)基金資助項目(11504417)

收稿日期:2014-09-23修改稿收到日期:2014-12-26

通信作者曾新吾 男,教授,博士生導(dǎo)師, 1963年3月生

中圖分類號:O43

文獻標(biāo)志碼:A

DOI:10.13465/j.cnki.jvs.2016.02.034

Numerical study on the focusing of sound wave using a concave ellipsoidal mirror

ZHANG Jun1, CHEN Peng1, CHEN Zheng-wu1, ZHAO Yun2, ZENG Xin-wu2

(1. Key Laboratory of Aerodynamic Noise Control, China Aerodynamics Research and Development Center, Mianyang 621000, China; 2. College of Opto-electronic Science and Engineering, National University of Defense Technology, Changsha 410073, China)

Abstract:Based on the Hamilton’s theoretical linear solution, the numerical results of the reflection waves of arc-discharge plasma sound source (AD-PSS) along the symmetric axis of an ellipsoidal mirror were presented. The propagation characteristics of the center wave, edge wave and wake wave were analyzed. In the front of the far focus of the mirror, the pressure of center wave is positive and those of the edge wave and wake wave are negative. Beyond the far focus, the situation is reversed. Under linear condition, the peak amplitude of the reflection wave collides with the far focus of mirror. The pressure amplitude of the reflection wave decays very quickly when going beyond the far focus. The reflection process of the sound wave was simulated using the FEM software COMSOL, and the evolution of reflection wave and the distribution of sound field were analyzed. Based on the khokhlov zabolotskava kuznetsov (KZK) equation, the influence of nonlinear coefficient on the propagation of reflection wave was discussed.

Key words:plasma sound source; reflection sound wave; ellipsoidal mirror

第一作者 張軍 男,博士,助理研究員,1983年11月生

郵箱:xwzeng@nudt.edu.cn

主站蜘蛛池模板: 亚洲欧美日韩中文字幕在线一区| 国产精品yjizz视频网一二区| 在线无码av一区二区三区| 成色7777精品在线| 国产欧美日韩免费| 高清无码一本到东京热| 免费欧美一级| 在线欧美国产| 欧美福利在线播放| 国产精品冒白浆免费视频| 成人夜夜嗨| 亚洲日韩精品欧美中文字幕| 日本日韩欧美| 国产成年女人特黄特色大片免费| 亚洲三级视频在线观看| 国产69精品久久久久妇女| 欧美成人午夜视频免看| 国产成人精品视频一区二区电影 | 中文字幕人妻无码系列第三区| 国产新AV天堂| 午夜国产不卡在线观看视频| 青草午夜精品视频在线观看| 91青青草视频| 欧美激情视频一区二区三区免费| 色一情一乱一伦一区二区三区小说| av大片在线无码免费| 亚洲综合婷婷激情| 色婷婷电影网| 国产精品久久久免费视频| 国产一区二区三区日韩精品| 国产女人18水真多毛片18精品| 日韩AV无码一区| 国产欧美精品午夜在线播放| 亚洲高清无码久久久| 亚洲欧美另类日本| 欧美中文字幕在线视频| 无码人中文字幕| 中文字幕色站| 欧美专区在线观看| 一本大道在线一本久道| 欧美日本在线| 久久天天躁狠狠躁夜夜2020一| 91在线中文| 91网址在线播放| 国产黄在线免费观看| 国产不卡一级毛片视频| 91麻豆精品国产91久久久久| 欧美在线黄| 97se亚洲综合在线天天| 91精品专区| 午夜性爽视频男人的天堂| 91视频区| 91小视频在线观看免费版高清| 在线观看亚洲国产| 国产成人久视频免费| 国产色偷丝袜婷婷无码麻豆制服| 看国产毛片| 欧美无专区| 国产精品女熟高潮视频| 日本亚洲成高清一区二区三区| 国产青榴视频在线观看网站| 蜜臀AV在线播放| 女同国产精品一区二区| 日韩天堂在线观看| 在线欧美国产| 久久久久国产一区二区| www欧美在线观看| 亚洲欧美综合在线观看| 粗大猛烈进出高潮视频无码| 亚洲码在线中文在线观看| 久久精品一卡日本电影| 国产二级毛片| 毛片基地美国正在播放亚洲 | 亚洲精品福利网站| 玖玖精品视频在线观看| 88av在线看| 久久先锋资源| 伊人天堂网| 香蕉色综合| 欧美日在线观看| 日本欧美精品| 欧美日韩中文国产va另类|