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

靜壓下球形空腔吸聲覆蓋層的建模與性能分析

2016-05-04 05:54:24何世平易少強
船舶力學 2016年7期
關鍵詞:有限元變形模型

張 沖,何世平,易少強

(1.海軍工程大學 動力工程學院,武漢 430033;2.海軍航空兵學院91899部隊,遼寧 葫蘆島 125001)

靜壓下球形空腔吸聲覆蓋層的建模與性能分析

張 沖1,2,何世平1,易少強1

(1.海軍工程大學 動力工程學院,武漢 430033;2.海軍航空兵學院91899部隊,遼寧 葫蘆島 125001)

基于聲波垂直入射下的二維理論,利用COMSOL軟件建立球形空腔吸聲覆蓋層單元模型,得到了在空腔受壓變形基礎上直接進行聲學性能計算的耦合模型,避免了以往根據變形量重新建模帶來的誤差,并與靜壓下穿孔率的理論公式進行對比,驗證了模型的有效性;利用所建模型分析了靜壓下球形空腔的吸聲性能,并比較了內部氣壓對空腔變形和吸聲系數的影響。結果表明:隨著靜壓的增大,峰值上移頻帶展寬,更顯著的是吸聲效果總體變差;內部氣壓使得空腔變形量減少、低頻吸聲效果有所提高。

吸聲覆蓋層;球形空腔;靜壓;內部氣壓;吸聲性能;COMSOL

0 引 言

吸聲覆蓋層技術是水下減振降噪技術中繼低噪聲螺旋槳技術和浮筏技術以后的又一大技術進步[1-2],對其開展理論研究有助于更好地了解吸聲覆蓋層的聲學特性,對工程設計具有重要的指導意義。目前對含空腔的吸聲覆蓋層的研究主要有解析方法、數值方法等[3-5],這些研究大都是在常壓條件下進行的,沒有考慮外部靜水壓力的影響。因此在理論計算和有限元仿真中,所得到的結果與實際測試結果存在很大差別。實際使用的吸聲覆蓋層是在幾百米深的水下工作的,其必然要受到靜水壓力的作用,聲學效果也會隨之發生改變。靜壓下吸聲覆蓋層吸聲性能發生改變的原因主要有兩個方面:一是材料的動態力學性能會隨著壓力而變化;二是吸聲覆蓋層受到靜水壓力時,空腔會發生形變。對于材料力學性能的研究,國內外許多學者針對類橡膠的粘彈性材料進行了很多測試和分析。呂林海[6]利用動態黏彈譜儀對高分子聚合物力學性能的頻變特性進行了測試,并基于得出的動態力學參數計算了圓柱形空腔吸聲覆蓋層的聲學性能;黃修長等[7]利用聲管測試的方法對材料的參數進行了測試,分析了各參數對諧振型吸聲結構的影響;鄒明松等[8]采用傳遞矩陣法建立了靜壓下吸聲覆蓋層的聲阻抗求解方法。而對于吸聲覆蓋層空腔變形的研究還處于起步階段,姜聞文等[9]利用ANSYS軟件對幾種橡膠結構的吸聲覆蓋層在靜壓下的吸聲性能進行了研究,比較了不同結構的空腔對吸聲性能的影響;Panigrahi等[10]基于有限元法對不同組合空腔在靜壓下的聲學性能進行了計算,并將結果與解析解、實驗作了對比分析;陶猛[11]推導出了吸聲覆蓋層聲學特能的簡化計算方法,計算了不同靜壓下吸聲覆蓋層的聲學性能。

當前的文獻研究都是基于有限元法,首先計算靜壓下吸聲覆蓋層單元的幾何變形,然后根據變形量重新建立規則形狀的吸聲覆蓋層單元模型。由于空腔在靜壓下的變形是無規則、無規律可循的,即使得到靜壓下的變形量,二次建模時也很難與實際的空腔變形相吻合,這對后續的聲學性能分析帶來一定的誤差;另一方面,關于靜壓下吸聲覆蓋層的所有研究都沒有考慮空腔內部氣壓的影響。本文鑒于以上兩方面的不足,利用COMSOL軟件建立了球形空腔的固體力學和聲—固耦合模型,然后在空腔受壓發生形變的模型上直接劃分網格,計算不同壓力下吸聲覆蓋層的吸聲性能,避免了二次建模;并進一步分析了靜壓下空腔內部氣壓對吸聲性能的影響。

1 物理模型

1.1 有限元物理模型的建立

具有球形空腔周期性分布的吸聲覆蓋層結構如圖1所示,在粘彈性體內有成正三角形排列的相互平行的球形空腔。把整個吸聲覆蓋層看成是由各個以球形空腔為中心的正六棱柱體組成,由于結構對稱、重復,只需要研究其中一個棱柱體中的波傳播即可。然而要建立一個六面棱柱體的理論模型,是一件相當困難的事情。為了計算簡便,本文通過幾何轉換用圓柱體近似代替六棱柱體[12],每個吸聲覆蓋層單元就可以簡化為一段有限長粘彈性圓柱管。吸聲覆蓋層中波的傳播與損耗問題可以模型化為粘彈性圓柱管中波的傳播與損耗,這為建立理論模型提供了極大的方便。

圖1 球形空腔周期性分布的吸聲覆蓋層結構Fig.1 The structure of the absorbing layer of the periodic distribution of the spherical cavities

三維周期邊界條件下的有限元模型雖然能完整地仿真吸聲覆蓋層的聲學特性,但由于三維模型的體單元節點數過多會造成單次運算效率較低,難以滿足多次迭代運算需求。考慮到吸聲覆蓋層單元具有軸對稱性,因此在建模時將體單元模型(三維周期邊界模型)簡化為面單元模型(二維軸對稱模型),如圖2所示。

圖2左側邊界處、上下兩端均超出有限元的虛線是模型的旋轉中心軸,有限元模型僅為完整單元的一個旋轉截面。灰色區域W、V分別代表入射端水介質、橡膠覆蓋層基體,空白區域C代表球形空腔。吸聲覆蓋層模型的邊界1為固定邊界;三維模型的周期性邊界條件轉化為二維軸對稱模型的法向量位移為0的邊界(即邊界2、3);邊界4、5、6為模型的軸對稱邊界;在有限元模型中通過設置聲—固耦合邊界(即邊界7)來處理水介質—橡膠基體的耦合作用;通過設置吸收邊界9解決無限水介質區域,即在W外端面法向方向添加平面波輻射條件,使有限厚度的有限元能代表無限厚的水介質;邊界8為球形空腔的自由邊界。

聲波從上端面沿旋轉軸方向垂直入射模型單元,會產生聲波的反射、折射、輻射、吸收等作用。通過積分入射聲波、吸收聲波的能量可求得吸聲覆蓋層在聲波垂直入射條件下的吸聲特性。

圖2 吸聲覆蓋層的單元模型Fig.2 The unit cell of anechoic coatings model

1.2 模型的有效性驗證

圓柱形空腔的吸聲覆蓋層有限元模型比較成熟,陶猛[13]、譚洪波等[14]已多次驗證了其有效性,但對球形空腔單元受壓后直接計算其吸聲性能的模型還沒驗證過其可行性。Gaunaurd[15]理論分析了受壓前后球形空腔吸聲覆蓋層的穿孔率與靜水壓力的關系,

式中:μ是材料的剪切模量,P代表靜水壓力,Φ代表穿孔率。

本文正是利用(1)式穿孔率的理論模型與COMSOL有限元模型作對比,以驗證模型的可行性。在如圖2所示的模型中,穿孔率可以表示為

式中:Vc(P)、Vv分別表示靜壓P下球形空腔的體積、橡膠基體的體積。

要證明模型的有效性,只要驗證(1)式、(2)式中的Φ是否相同即可。但在靜壓P下Vc(P)是變化的,要想計算Vc(P),必須借助于Stokes公式(3)將曲面積分轉化為曲線積分

在圖示的二維軸對稱模型中,主要是對空腔內表面進行曲線積分,

式中:Asphere代表球形空腔的曲線積分,z代表球形空腔內表面每一點的z向坐標,nz代表z坐標的法向矢量,lsphere代表球形空腔的曲線。

得到空腔內表面的曲線積分后,再通過COMSOL中的“計算旋轉幾何的積分”選項就可以將求得的Asphere轉化為球形空腔受壓后的體積,即Vc(P)。在計算時橡膠基體的泊松比在0.499左右,為幾乎不可壓縮材料,故Vv是定值。

在驗證計算時,采用的橡膠材料的參數是隨頻率的改變而發生變化的,參數表達式根據實驗數據擬合得到[16],橡膠材料的楊氏模量(單位:Pa)、損耗因子和泊松比的表達式為:

由于橡膠的楊氏模量為復楊氏模量,只要知道復楊氏模量的實部E和虛部η就可以求得復楊氏模量,故可表示為:

以上各式中,f代表頻率,i代表虛數。

COMSOL有限元與解析公式通過計算得到的結果如圖3所示。從圖中可以看出,有限元與解析解的穿孔率變化曲線吻合較好。至于在靜壓較大時,有限元結果與解析公式的穿孔率存在一定的差別,主要是由于解析公式在推導的過程中,假設計算單元中含有均勻分布的小氣泡,而且氣泡在受壓變形后仍是規則的球形,故理論推導的假設與空腔受壓后呈橢球形的實際變形存在誤差。

圖3 有限元與理論結果比較圖Fig.3 Comparison of FEM and theoretical results

2 不同壓力下吸聲覆蓋層的形變及吸聲性能

圖4 不同靜壓下單元形變及應力圖Fig.4 Deformation and stress of unit under different hydrostatic pressure

與目前靜壓下吸聲覆蓋層有限元研究不同,本文在利用COMSOL建模時,選用“固體力學”、“聲—固耦合”和“移動網格”三大模塊相結合,主要好處是單元模型在計算靜壓變形后,可以依靠“移動網格”直接將計算的變形應用于聲—固耦合計算中,不需要根據變形量二次建模,有效地提高了計算精度。如圖4所示為1 MPa、3 MPa和5 MPa下的靜壓變形,吸聲覆蓋層的聲學特性的計算就是在此變形基礎上直接進行的。從圖中可以看出,吸聲覆蓋層隨著壓力的增大,球形空腔的形變量和覆蓋層厚度的壓縮量越大,且應力能越大。

圖5 不同靜壓下的吸聲系數曲線Fig.5 Absorption coefficient curves under different hydrostatic pressure

圖5比較了吸聲覆蓋層在不同靜壓下的吸聲系數,可以看出,隨著靜壓的增大,吸聲系數總體呈下降趨勢,低頻吸聲性能變差,且峰值向高頻移動,頻帶增寬100 Hz左右。從共振理論分析看,由于空腔體積和覆蓋層厚度被壓縮,等效阻抗失效,反射增強,使得共振峰上移頻帶展寬,而且低頻吸聲效果隨之下降。目前許多學者研究的圓柱形空腔,雖然表明在靜壓下吸聲峰值發生移動,但在高頻段靜壓越大吸聲系數也越大[9],即壓力大的吸聲系數曲線要穿透壓力小的曲線有上揚的走向。而本文研究球形空腔卻發現,吸聲系數曲線隨著壓力的增大總體呈下降趨勢(如在5 MPa時,吸聲系數比常壓,即0 MPa時,總體下降高達近30%),這一點與先前學者研究的圓柱形空腔存在顯著的差異。

球形空腔與圓柱形空腔雖然結構不同,但吸聲系數走向不會發生太大的不同,趨勢應該一致[17]。所以COMSOL計算結果與其他學者研究的結論存在差別的主要原因是:由于單元模型在靜壓下的變形呈不規則橢球形(圖4所示),如若用計算的變形量重新建立模型,不可能與實際變形相對應;本文在計算聲學特性時,單元采用的是受壓變形后的實際模型,避免了以往根據變形量重新建模的誤差,而這種不規則的形變造成了對聲波吸收能力的減弱、對聲波反射能力的增強,故計算的吸聲系數總體下降。同時,從工程設計方面表明,若僅僅要求常壓(即無靜壓作用)下設計的吸聲覆蓋層滿足工程需要,那么在水下使用時會達不到想要的吸聲效果,因為在水下不僅吸聲頻帶發生變化,吸聲性能也會明顯地減弱。

3 靜壓下空腔內部氣壓的影響

無論圓柱形空腔還是球形空腔,內部不是真空,都是充滿空氣的。吸聲覆蓋層在靜壓作用下被壓縮,空腔體積會縮小,內部氣體壓力會隨之增大。那么靜壓下,空腔內部氣壓會對吸聲覆蓋層的形變和吸聲性能產生怎樣的影響,至今沒有相關文獻對此進行過研究。本文在利用COMSOL有限元建模時,把內部氣壓的作用考慮在內,分析其對單元形變和吸聲系數的影響。

受壓前、后空腔內氣體壓力存在如下關系:

式中:p0、p分別代表壓縮前、后空腔內的壓強,ρ0、ρ分別代表壓縮前、后空腔內空氣的密度,A0、A分別代表壓縮前后二維球形空腔的面積,γ代表氣體常數。其中,p0=0.1 MPa是標準大氣壓力,γ=1.4。

隨著壓力的增加,就會有剩余的空氣壓力作用在空腔內壁,對空腔形變產生一定的作用。由(5)式可知,剩余的空氣壓力為:

在COMSOL模型中,將△p作為邊界載荷加載在空腔內壁上。圖6和圖7所示分別為有、無空氣壓力作用下,空腔體積的減少量相對于原空腔體積的變化率和靜壓為5 MPa時吸聲系數曲線。

圖6 相對體積變化率Fig.6 Relative volume change rate curves

圖7 吸聲系數曲線比較圖Fig.7 The comparison of absorption coefficient curves

通過相對體積比較圖可以看出,在空腔內部,有空氣壓力的變形量要比無空氣壓力的變形量小,原因是空腔受到靜壓壓縮時,內部氣體會給空腔壁沿著法向向外作用的支撐力,以阻礙空腔的壓縮變形。通過吸聲系數比較圖可以看出,在低頻段(1 600 Hz以下頻段),有空氣壓力的吸聲系數要大于無空氣壓力的吸聲系數,高頻段正好相反。這是因為無空氣壓力的孔腔被壓縮量大,吸聲覆蓋層基體上部的等效阻抗不再是漸變的,與水阻抗失配的增加,使得反射更加增強,孔腔的固有頻率也會隨之升高,低頻吸聲效果故而比有空氣壓力的要差,高頻反而要強[18]。

4 結 語

本文基于聲波垂直入射條件下吸聲覆蓋層的解析理論,利用COMSOL有限元軟件建立靜水壓力下球形空腔模型,將計算的穿孔率與現有解析公式進行了對比,驗證了模型的有效性,然后計算了不同靜壓下模型單元的變形量,并在單元形變的基礎上采用移動網格模塊直接計算聲—固耦合的吸聲系數,最后討論了空腔內部氣壓對空腔形變和吸聲性能的影響。結果表明:利用移動網格直接在單元形變基礎上計算吸聲系數時,峰值向高頻移動,頻帶變寬,而且與其他學者根據變形量二次建模得到的結果相比,吸聲效果隨著靜壓的增大整體呈下降趨勢;內部氣壓的作用可以使得空腔變形量減小,低頻吸聲性能有所提升。

[1]孟曉宇,肖國林,陳 虹.國外潛艇聲隱身技術現狀與發展綜述[J].艦船科學技術,2011,33(11):135-139. Meng Xiaoyu,Xiao Guolin,Chen Hong.Review of the present situation and development of acoustic stealth technology for submarines abroad[J].Ship Science and Technology,2011,33(11):135-139.

[2]蘇 強,王桂波,朱鵬飛,宋 楊.國外潛艇聲隱身前沿技術發展綜述[J].艦船科學技術,2014,36(1):1-9. Su Qiang,Wang Guibo,Zhu Pengfei,et al.Review of foreign submarine acoustic stealth frontier technologies development[J].Ship Science and Technology,2014,36(1):1-9.

[3]湯渭霖,何世平,范 軍.含圓柱形空腔吸聲覆蓋層的二維理論[J].聲學學報,2005,30(4):289-295. Tang Weilin,He Shiping,Fan Jun.Two-dimensional model for acoustic absorption of viscoelastic coating containing cylindrical holes[J].Acta Acustica,2005,30(4):289-295.

[4]何世平,湯渭霖,何 琳,等.變截面圓柱形空腔覆蓋層吸聲系數的二維近似解[J].船舶力學,2006,10(1):120-127. He Shiping,Tang Weilin,He Lin,et al.Analysis of acoustic characteristics of anechoic coating containing varying sectional cylindrical cavity[J].Journal of Ship Mechanics,2006,10(1):120-127.

[5]Easwaran V,Munjal M L.Analysis of reflection characteristics of a normal incidence plane wave on resonant sound absorbers:A finite element approach[J].Journal of the Acoustical Society of America,1993,3(3):1308-1318.

[6]呂林梅,溫激鴻,趙宏剛,溫熙森.黏彈性材料的動態力學特性及其對覆蓋層吸聲性能的影響[J].物理學報,2014,63 (15):1-6. Lü Linmei,Wen Jihong,Zhao Honggang,et al.Dynamical mechanical property of viscoelastic materials and its effect on acoustic absorption of coating[J].Acta Phys.Sin.,2014,63(15):1-6.

[7]黃修長,朱蓓麗,胡 碰,華宏星.靜水壓力下橡膠動態力學參數的聲管測量方法[J].上海交通大學學報,2013,47 (10):1503-1508. Huang Xiuchang,Zhu Beili,Hu Peng,et al.Measurement of dynamic properties of rubber under hydrostatic pressure by water-filled acoustic tube[J].Journal of Shanghai Jiaotong University,2013,47(10):1503-1508.

[8]鄒明松,吳文偉,余曉麗,廖彬彬.靜水壓下聲學覆蓋層聲阻抗研究[J].艦船科學技術,2013,35(3):57-60. Zou Mingsong,Wu Wenwei,Yu Xiaoli,et al.Calculation of acoustic coating’s impedance under hydrostatic pressure[J]. Ship Science and Technology,2013,35(3):57-60.

[9]姜聞文,陳光冶,朱 彥.靜水壓變化下橡膠結構吸聲性能的計算與分析[J].噪聲與振動控制,2006(5):55-57. Jiang Wenwen,Chen Guangye,Zhu Yan.Computation and analysis of sound absorption performance of rubber structures under variable hydraulic pressure[J].Noise and Vibration Control,2006(5):55-57.

[10]Panigrahi S N,Jog C S,Munjal M L.Multi-focus design of underwater noise control linings based on finite element analysis[J].Applied Acoustics,2008,69(12):1141-1153.

[11]陶 猛,卓琳凱.靜水壓力下吸聲覆蓋層的聲學性能分析[J].上海交通大學學報,2011,45(9):1340-1344. Tao Meng,Zhuo Linkai.Effect of hydrostatic pressure on acoustic performance of sound absorption coating[J].Journal of Shanghai Jiaotong University,2011,45(9):1340-1344.

[12]何世平.粘彈性圓柱管中波的傳播研究及吸聲覆蓋層聲學特性研究[D].上海:上海交通大學,2004. He Shiping.Research of wave propagation in viscoelastic tube and analysis of acoustic characteristics of anechoic coating [D].Shanghai:Shanghai Jiaotong University,2004.

[13]陶 猛,卓琳凱.基于ANSYS的吸聲覆蓋層聲學性能計算與分析[J].振動與沖擊,2011,30(1):87-90. Tao Meng,Zhuo Linkai.Simulation and analysis for acoustic performance of a sound absorption coating using ANSYS software[J].Journal of Vibration and Shock,2011,30(1):87-90.

[14]譚紅波,趙 洪,徐海亭.有限元法分析空腔周期分布粘彈性層的聲特性[J].聲學學報,2003,28(3):277-282. Tan Hongbo,Zhao Hong,Xu Haiting.Sound characteristics of the viscoelastic layer containing periodic cavities by the finite element method[J].Acta Acustica,2003,28(3):277-282.

[15]Gaunaurd G,Callen E,Barlow J.Pressure effects on the dynamic effective properties of resonating perforated elastomers [J].J Acoust.Soc.Am.,1984,76(1):173-177.

[16]胡 碰.靜水壓力下聲學覆蓋層聲學性能模塊化方法研究[D].上海:上海交通大學,2008. Hu Peng.Research of acoustic coating’s properities under hydrostatic pressure by modular-FEM method[D].Shanghai: Shanghai Jiaotong University,2008.

[17]姚熊亮,劉文賀,劉慶杰,張 妍.水深與腔型對隔聲去耦瓦吸聲系數的影響[J].哈爾濱工程大學學報,2007,28(6): 605-610. Yao Xiongliang,Liu Wenhe,Liu Qingjie,et al.Influence of depth and cavity shape on absorption coefficients of soundisolating decoupled tiles[J].Journal of Harbin Engineering University,2007,28(6):605-610.

[18]羅 忠,朱 錫,林志駝,王衛忠.水下吸聲覆蓋層結構及吸聲機理研究進展[J].艦船科學技術,2009,31(8):23-30. Luo Zhong,Zhu Xi,Lin Zhituo,et al.A review of underwater anechoic coating structure and absorption theories[J].Ship Science and Technology,2009,31(8):23-30.

Model and absorption performance of anechoic coating embedding sphere cavities

ZHANG Chong1,2,HE Shi-ping1,YI Shao-qiang1
(1.College of Power Engineering,Naval Univ.of Engineering,Wuhan 430033,China; 2.Unit No.91899,Naval Air Force Academy,Huludao 125001,China)

Based on the 2D theory of sound wave normally impinging on the absorption layer,the acousticsolid interaction model of anechoic coating which contains sphere cavities was made by COMSOL software and the model can calculate absorption coefficient directly instead of remodeling according to the volume of deformation under hydrostatic pressure.Perforation rate was compared between COMSOL and analytical solution in order to prove the validity of the model.The absorption property and influence of the internal air pressure were analysed based on the interaction model.The results indicate that with increasing of hydrostatic pressure the peak frequency is towards to the high frequency and becomes wider,the more difference is that effect of sound absorption becomes worse;internal air pressure makes the deformation of the cavity decrease and absorption coefficient in low frequency increase.

anechoic coating;sphere cavity;hydrostatic pressure;internal air pressure; absorption property;COMSOL

TB564

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

1007-7294(2016)07-0909-08

2016-03-19

國家自然科學基金(11374369)

張 沖(1987-),男,碩士生,E-mail:chzhang0101@foxmail.com;何世平(1972-),男,副教授,碩士生導師。

猜你喜歡
有限元變形模型
一半模型
重要模型『一線三等角』
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
重尾非線性自回歸模型自加權M-估計的漸近分布
“我”的變形計
例談拼圖與整式變形
會變形的餅
3D打印中的模型分割與打包
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 77777亚洲午夜久久多人| 91精品日韩人妻无码久久| 欧美成人免费午夜全| 9999在线视频| 国产成人精品免费av| 欧美日韩国产综合视频在线观看| 色婷婷在线播放| 欧美日韩资源| 青草精品视频| 亚洲系列中文字幕一区二区| 免费在线a视频| 草草影院国产第一页| 在线观看无码av免费不卡网站 | 亚国产欧美在线人成| 日韩乱码免费一区二区三区| 欧美日韩午夜| 国产美女视频黄a视频全免费网站| 国产日韩欧美黄色片免费观看| 69av免费视频| 久久精品无码国产一区二区三区 | 国产精品第一区在线观看| 国产打屁股免费区网站| 欧美精品啪啪一区二区三区| 伊人久久婷婷五月综合97色| 免费AV在线播放观看18禁强制| 永久免费av网站可以直接看的 | 国产精品妖精视频| 欧美色图第一页| 亚洲无线视频| 少妇露出福利视频| 久久国产精品电影| 在线观看的黄网| 国产91线观看| 国产福利不卡视频| 91久久国产综合精品女同我| 大学生久久香蕉国产线观看| 国内嫩模私拍精品视频| 无码福利日韩神码福利片| 亚洲成A人V欧美综合| 精品视频一区二区三区在线播| 中文字幕永久在线看| 国产熟睡乱子伦视频网站| 原味小视频在线www国产| 欧美天天干| 国产最爽的乱婬视频国语对白 | 亚洲妓女综合网995久久| 久久人搡人人玩人妻精品| 午夜精品久久久久久久99热下载| 国产日韩欧美精品区性色| 自拍偷拍欧美日韩| 国产成人精品男人的天堂| 国产日韩欧美视频| 亚洲精品国产乱码不卡| 欧美精品影院| 色婷婷亚洲十月十月色天| 国产精品私拍在线爆乳| 中国国语毛片免费观看视频| 亚洲人成网址| 国产成人高清精品免费| 美女裸体18禁网站| 91在线国内在线播放老师| 99久久精品国产精品亚洲 | 亚洲高清国产拍精品26u| 国产一级无码不卡视频| 色成人亚洲| 一级毛片在线免费看| 91精品国产无线乱码在线| 国产黑人在线| 大香网伊人久久综合网2020| 亚洲成a人片77777在线播放| 日本高清免费不卡视频| 一级毛片免费的| 国内99精品激情视频精品| 亚洲av日韩av制服丝袜| 免费看一级毛片波多结衣| 日本高清成本人视频一区| 手机在线免费不卡一区二| 国产情侣一区二区三区| 日本午夜网站| 国产精品99在线观看| 日本精品αv中文字幕| 色婷婷久久|