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

水下加筋圓柱殼聲輻射特性分析的改進精細傳遞矩陣法

2017-11-30 05:49:47龐福振王獻忠姚熊亮
振動與沖擊 2017年22期
關鍵詞:振動結構

龐福振, 吳 闖, 王獻忠, 姚熊亮

(1. 哈爾濱工程大學 船舶工程學院,哈爾濱 150001; 2. 海軍裝備研究院,北京 100161;3. 上海交通大學 機械系統(tǒng)與振動國家重點實驗室,上海 200240)

水下加筋圓柱殼聲輻射特性分析的改進精細傳遞矩陣法

龐福振1,2, 吳 闖3, 王獻忠1, 姚熊亮1

(1. 哈爾濱工程大學 船舶工程學院,哈爾濱 150001; 2. 海軍裝備研究院,北京 100161;3. 上海交通大學 機械系統(tǒng)與振動國家重點實驗室,上海 200240)

基于Flügge殼體振動理論,將改進精細傳遞矩陣法應用于水下加筋柱殼聲輻射問題的求解;通過應用加法定理和增量存儲,改進非齊次項求解方法,構造了加筋圓柱殼的輻射噪聲改進精細傳遞矩陣法。分析了邊界條件、結構損耗因子、流體介質以及殼體厚度對加筋圓柱結構聲輻射的影響。研究表明,當邊界條件從自由、簡支到剛固,結構對外輻射聲壓逐漸降低;隨著結構損耗因子的增大,結構對外輻射聲壓大體上逐漸降低。

非齊次項; 精細積分傳遞矩陣; 加筋圓柱殼; 振動; 聲輻射

在實際工程中,特別是在現(xiàn)代國防軍事領域,加筋圓柱殼是眾多武器裝備的基本簡化模型,如:魚雷、導彈、潛艇等。因此對于加筋圓柱殼的聲振耦合機理分析,一直以來是國內外學者的研究熱點。在數(shù)值方面,F(xiàn)EM、BEM方法、SEA方法等方法得到了廣泛應用[1]。數(shù)值方法理論上來說能夠處理任意復雜的結構,但其求解精度受限于計算頻段,建模周期長,而且數(shù)值方法難以進行機理分析,這些均制約了有限元技術的發(fā)展。而傳統(tǒng)的解析方法只能給出簡單結構的解析解。因此半解析半數(shù)值法逐漸成為進行加筋圓柱殼聲輻射分析的一種有效方法。其中,傳遞矩陣法由于其對求解鏈狀結構動力問題的巨大優(yōu)勢,而受到國內外學者青睞。

Tottenham等[2]首先給出了一種求解圓柱殼結構自由振動的傳遞函數(shù)分析方法。后來Irie等[3]進行進一步推廣應用,形成了較為完善的傳遞矩陣方法,并將其成功應用到軸向具有離散彈簧支承性殼體結構的自由振動問題上,這些研究均是圍繞結構的自由振動特性進行展開研究的,并沒有涉及到殼體結構振動響應和水下聲輻射的問題。蔡顯新[4]以旋轉殼結構為計算模型,借鑒了Irie等的思想,將旋轉殼的狀態(tài)向量沿周向以級數(shù)形式展開,推導了旋轉殼的一階振動微分方程,給出了一種求解旋轉殼自由振動的半解析解,然而該方法精度有限,要求離散的分段較細。曹雷等運用Riccati傳遞矩陣法分析了水下有限長環(huán)肋圓柱殼的聲輻射性能,但是在計算一階非齊次矩陣微分方程的非齊次項時,采用了分段插值理論對非齊次項進行多項式逼近,難免會出現(xiàn)數(shù)值不準確,精度損失的問題。

為此,針對以上的不足,提出一種精細積分傳遞矩陣法。該方法主要用來解決結構振動響應和水下聲輻射的問題,對于文獻[5],在計算一階非齊次矩陣微分方程的非齊次項時,采用分段插值理論對非齊次項多項式逼近的方法所帶來數(shù)值不準確,精度損失的問題,本文采用了一種精度更高的處理方法,即譚述君等[6]提出的非齊次項精細積分法。該方法將柱殼微段劃分為精細積分步,基于積分步內的Taylor級數(shù)展開和加法定理,循環(huán)處理得到柱殼微段內非齊次項的高精度結果,數(shù)值計算時精度更高更準確。

1 理論推導

1.1 圓柱殼的場傳遞矩陣

本文研究的模型為浸沒在無限流場介質中的有限長加筋圓柱殼,其中假設圓柱殼兩端為無限大剛性障板,模型如圖1所示。

圖1 加筋圓柱殼模型示意圖Fig.1 Model of stiffened cylindrical shell

殼體變形通常采用基于直線假設的薄殼理論進行描述,為了得到精準的數(shù)值,本章采用相對比較精確的Flügge殼體理論[7],柱殼的微元進行受力分析,得到力的平衡方程。由于本文是基于動力學理論建立的方程,很多項中都包含時間項,為了方便書寫和推導,在下文的書寫中將略去動響應的時間項e-iωt, 圓柱殼坐標系(γ,φ,x)以及位移正方向,如圖2所示。

圖2 圓柱殼坐標系Fig.2 Cylindrical shell coordinates system

根據(jù)Flügge殼體理論,殼體力平衡方程為

(1)

(2)

(3)

(4)

通過消除8個未知量,保留8個未知量,將各量無量綱化并沿周向展成三角函數(shù)形式,經(jīng)過化簡得到一階矩陣微分方程

(5)

1.2 環(huán)肋的點傳遞矩陣

由于環(huán)肋的作用,圓柱殼的狀態(tài)向量在環(huán)肋處發(fā)生改變,利用環(huán)肋與殼體連接處的連續(xù)條件和環(huán)肋的運動控制方程(面內的拉伸和彎曲振動,面外的彎曲和扭轉振動)很容易得到環(huán)肋的點傳遞矩陣Rk。圓柱殼上某一環(huán)肋位置處ξk處,由于環(huán)肋的存在,環(huán)肋左端、右端的殼體中兩個面內力和兩個面外力發(fā)生變化,狀態(tài)向量滿足式(6)

(6)

式中,Rk為8×8的點傳遞矩陣,矩陣中各非零元素為Rii=1,i=1,2,…,8。

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

式中:Rb、I1、I2、IP、J分別為環(huán)肋的中和軸半徑、對縱向對稱軸的慣性矩、對徑向對稱軸的慣性矩、極慣性矩以及扭轉常數(shù);A,ρ,G分別為環(huán)肋截面積、密度、剪切模量;E為復合彈性模量;E=E0(1+iη),E0為彈性模量;η為損耗因子。Rb=R+e;e為偏心距,內肋取負號,外肋取正號。

1.3 非齊次項動力響應的處理問題

對于線性非齊次微分方程式(5),其一般解為

(15)

式中:r(τ)={F(ξ)}-{p(ξ)}右端第二項為非齊次項引起的狀態(tài)響應;對上式各項的求解是解決聲輻射問題的第一步,也是非常重要的一步。eUΔξ的求解可參考文獻[8]的指數(shù)矩陣精細積分法。因此,只要解出非齊次項就可以得到狀態(tài)向量Z(ξ), 采用精細積方法處理非齊次項。

外流場用Helmholtz方程表示為

(16)

聲壓滿足無窮遠處邊界條件和端面邊界條件為

(17)

(18)

Helmholtz方程可采用分離變量的辦法,考慮式(17)、式(18)等邊界條件,則周圍流體域的輻射聲壓可表示為

(19)

(20)

式中,e為已知項。

α=ikml

(21)

pp=ecos(kmlξk+1)

(22)

qq=esin(kmlξk+1)

(23)

則式(20)可改寫為

(24)

由加法定理可得

(25)

同理可得

Φ0(2Δξ)=(Φ0(Δξ))2

(26)

φ0(2Δξ)=(φ0(Δξ))2

(27)

當基本區(qū)段Δξ劃分為精細區(qū)段τ(τ=Δξ/2M)時,M建議取20。E(τ)、Φ0(τ)、φ0(τ)可采用Taylor級數(shù)展開有限項進行近似

(28)

(29)

(30)

(31)

(32)

根據(jù)式(25)可得

(33)

令E1為E(τ)的實部

(34)

令E2為E(τ)的虛部

(35)

(36)

(37)

(38)

(39)

(40)

至此,已利用精細區(qū)段的加法定理求出非齊次項的精細結果。同時可將聲壓作用下的式(15)寫成矩陣的形式

(41)

集中力對圓柱殼的作用與環(huán)肋類似,僅僅改變了集中力處得狀態(tài)向量,故集中力左右兩端滿足式(42)

(42)

式中:Fk=RK-1[0 0 0 0 0fn0 0]T;IFk為集中力的點傳遞矩陣。

1.4 加筋圓柱殼聲輻射的求解

1.4.1 集中力作用下的狀態(tài)向量

(43)

1.4.2 聲壓作用下的狀態(tài)向量

(44)

1.4.3 聲輻射的求解

根據(jù)線性疊加原理,對應于任意階周向波數(shù)n下的徑向位移滿足

(45)

但是pmn是未知的,故必先求解pmn才能得到wn(x)。

根據(jù)連續(xù)性條件,在流體與結構的接觸面上,流體徑向速度必須等于結構的徑向速度

(46)

由于結構界面處上任意一點均滿足連續(xù)條件,因而可在結構長度方向進行取M個點,其中M既要大于波數(shù)m,又要滿足M=2πλ-1。取M個點代入式(46)中,并將式(19)、式(45)代入式(46)中,則可轉化為求解各波數(shù)對應的聲壓系數(shù)的方程

(47)

可通過Moore-Penrose廣義求逆法求解。具體過程為pmn

(48)

(49)

顯然上式中只有一個未知的列向量{pm1}n,所以可以求解。

[U]n{p}n={Q}n

(50)

將[U]進行奇異值分解,可得

[U]=[A][D][V]T

(51)

對矩陣[U]進行Moore-Penrose求逆,可得

[U]-1=[V][D]-1[A]T

(52)

將其代入式(24),可得

(53)

將求得的系數(shù)矩陣{p}代入式(19),可得到流場中的輻射聲壓。由水下輻射聲壓可求解對應的輻射聲壓級:Lp=20 lg(p/p0),參考基準聲壓值為1 μPa。圓柱殼上任意一點的徑向位移為

(50)

結構表面徑向均方振速為

(51)

速度級為

速度級的基準為

2 PITMM的驗證與數(shù)值計算

2.1 驗證模型簡介

采用文獻[9-10]和文獻[5]中的數(shù)據(jù)對本文方法的準確性進行驗證。對比結果見表1、表2。

文獻[9]中參數(shù)為:材料E0=2.06×1011Pa; 泊松比μ=0.3; 密度ρ=7 800 kg/m3; 損耗因子取η=0.01; 密度ρ0=1 000 kg/m3; 聲波傳播速度取c0=1 500 m/s; 柱殼幾何參數(shù)長L=0.6 m; 半徑R=0.2 m; 厚度h=0.003 m; 柱殼的內肋為0.002 m×0.030 m; 環(huán)肋數(shù)量9根。徑向激勵力作用在殼體(L/2,0,R)處,幅值為1 N。 測量點為(L/2,θ,2 m)。

文獻[10]中參數(shù)為:材料、流體介質參數(shù)同文獻[9]。柱殼幾何參數(shù)長L=0.6 m; 半徑R=0.175 m; 厚度h=0.002 m; 柱殼的外肋0.002 m×0.025 m; 環(huán)肋數(shù)量5根。 徑向激勵力作用在殼體(L/2,0,R)處, 幅值為1 N。測量點(L/2,θ,1.175 m), 如圖3所示。

圖3 驗證模型示意圖Fig.3 Model schematic diagram

2.2 計算結果的對比

(1) 算例1與公開文獻數(shù)據(jù)對比。

(2) 算例2與公開文獻數(shù)據(jù)對比。

文獻[9]采用能量法,用Hamilton原理和Green函數(shù)方法研究環(huán)肋柱殼的聲輻射特性;文獻[10]基于Flügge殼體理論,僅考慮徑向激振力,計及殼體受環(huán)肋的反力和反力矩作用,分析加筋柱殼的聲輻射特性。通過表1、表2中試驗值與計算值對比分析可知,本文方法較文獻[5]、文獻[9]方法精度要高,與文獻[10]方法精度相當,除精度高外,本文方法最大特點是易于編程,計算效率高。

表1 簡支加筋柱輻射聲壓級對比

表2 簡支加筋柱殼輻射聲壓級對比

3 加筋柱殼結構振動噪聲特性研究

采用文獻[10]的模型參數(shù),基于非齊次項精細積分傳遞矩陣法分析加筋柱殼的一些參數(shù)對聲輻射的影響。除“3.1”節(jié)算例外,其余算例中模型邊界條件均為兩端簡支,模型自振頻率計算值如表3所示。表3中,m為縱向半波數(shù),n為周向波數(shù)。

3.1 邊界條件對加筋柱殼輻射噪聲的影響

加筋圓柱殼模型兩端邊界分別是兩端自由(F-F)、兩端簡支(S-S)、兩端剛固(C-C)。聲壓測量點為(L/2,θ,1.175 m),即柱殼長度中點處周向一系列點。模型計算頻率點為f=4 000 Hz,對比結果見圖4。

圖4 各測點不同邊界輻射聲壓級圖線Fig.4 Circumferential cylindrical shell SPL

圖4比較了兩端不同邊界條件對水下加筋柱殼輻射聲壓級的影響。從圖4可知,兩端邊界條件從自由、簡支、到剛固,隨著固定的加強,結構對外輻射聲壓逐漸降低,這是因為兩端固定越強結構產(chǎn)生振動越困難,輻射聲壓也就越小。

3.2 結構損耗因子對加筋柱殼輻射聲壓級的影響

加筋圓柱殼模型的結構損耗因子分別為η=0.01;η=0.05;η=0.1。聲壓測量點為(L/2,θ,1.175 m),即柱殼長度中點處周向一系列點。模型計算頻率點為f=4 000 Hz。對比結果見圖5。

圖5 各測點不同阻尼輻射聲壓級圖線Fig.5 Circumferential cylindrical shell SPL

圖5比較了不同結構損耗因子對水下加筋柱殼輻射聲壓級的影響。由圖5可知,當損耗因子較小時(η=0.01和η=0.05),其對加筋柱殼水下輻射聲壓影響不明顯,而當損耗因子較大時(η=0.1),其對加筋柱殼水下輻射聲壓影響顯著。從聲壓分布來看,在靠近激勵力處(θ=0°),損耗因子對結構輻射聲壓影響小;在遠離激勵力處(θ=180°),損耗因子對結構輻射聲壓影響大。

3.3 流體介質對加筋柱殼輻射聲壓級的影響

水的質量密度為1 000 kg/m3;聲在水中的傳播速度1 500 m/s;空氣的質量密度為1.293 kg/m3;聲在空氣中傳播速度為340 m/s。聲壓測量點為(L/2,θ,1.175 m),即柱殼長度中點處周向一系列點。模型計算頻率點為f=4 000 Hz。對比結果見圖6。

圖6 各測點不同介質輻射聲壓級圖線Fig.6 Circumferential cylindrical shell SPL

圖6比較了不同流體介質對加筋柱殼輻射聲壓級的影響,由圖6可知,結構在流體介質阻抗大的水中比在體介質阻抗小的空氣中對外輻射聲壓明顯要大,因為水的阻抗相對空氣要大,聲波造成水介質的壓強變化量大于空氣,故出現(xiàn)此現(xiàn)象。

3.4殼體厚度對加筋柱殼徑向均方振速級和輻射聲壓級的影響

加筋圓柱殼模型的殼體厚度分別為h=0.002 m;h=0.01 m;h=0.02 m。結構表面徑向均方振速測量點(L/2,0,R),流場中聲壓測量點(L/2,0,1.175 m),模型計算頻域0~1 000 Hz。對比結果見圖7、圖8。

圖7 不同殼體厚度加筋柱殼徑向均方振速級的對比Fig.7 Lv(db) comparison of stiffened cylindrical shell with different thickness

圖8 不同殼體厚度加筋柱殼輻射聲壓級的對比Fig.8 SPL(db) comparison of stiffened cylindrical shell with different thickness

圖7、圖8比較了不同殼體厚度對水下加筋柱殼徑向均方振速級和輻射聲壓級的影響。從圖7可知,隨著殼體厚度的增加徑向均方振速級的峰值點向后推移,峰值大小基本無變化;在非峰點區(qū),隨著殼體厚度的增加,加筋柱殼徑向均方振速逐漸減小。從圖8可知,隨著殼體厚度的增加輻射聲壓級的峰值點也向后推移;在非峰點區(qū),基本滿足隨著殼體厚度的增加,加筋柱殼輻射聲壓級逐漸減小的規(guī)律。對于峰值點向后推移,是因為隨著殼體厚度的增加,結構剛度增加,其自振頻率也隨之增加。

4 結 論

針對傳統(tǒng)傳遞矩陣法以及精細傳遞矩陣法在處理非齊次項的不足,通過將非齊次方程Duhamel項精細積分法引入精細傳遞矩陣法,構造了求解水下加筋柱殼聲輻射問題的改進精細傳遞矩陣法,該方法采用積分步內Taylor級數(shù)展開和加法定理,具有計算精度高、易于編程的優(yōu)點。通過將本文方法計算結果與公開文獻、有限元等計算結果進行比較,以驗證本文方法的有效性。在此基礎上,分析了加筋柱殼邊界條件、結構損耗因子、流體介質及殼體厚度等參數(shù)的改變對其水下聲輻射特性的影響規(guī)律。其中,流體介質的改變對加筋柱殼水下聲輻射影響最大,結構在阻抗大的介質中聲輻射明顯大于阻抗小的介質;邊界固定越強,加筋柱殼向外輻射能力越低,輻射聲壓越小;損耗因子越大,吸收結構振動的能量越多,造成加筋柱殼向外輻射聲壓越小;殼體厚度對柱殼輻射聲壓的影響,基本滿足隨著殼體厚度增加輻射聲壓逐漸減小的規(guī)律。

[ 1 ] 瞿葉高, 華宏星, 孟光,等.基于區(qū)域分解的圓錐殼-圓柱殼-圓錐殼組合結構自由振動[J].振動與沖擊,2012, 31(22): 1-7.

QU Yegao, HUA Hongxing, MENG Guang, et al. A domain decomposition method for free vibration analysis of joined conical-cylindrical-conical shell [J]. Journal of Vibration and Shock, 2012, 31(22): 1-7.

[ 2 ] TOTTENHAM H, SHIMIZU K. Analysis of the free vibration of cantilever cylindrical thin elastic shells by the matrix progression method [J]. International Journal of Mechanical Sciences, 1972, 14(5): 293-310.

[ 3 ] IRIE T, YAMADA G, MURAMOTO Y. Free vibration of joined conical-cylindrical shells [J]. Journal of Sound and Vibration, 1994, 95(1): 31-39.

[ 4 ] 蔡顯新.計算旋轉薄殼自由振動的傳遞矩陣法[J].航空學報,1992,13(1): 64-70.

CAI Xianxin. A TMM the free vibration of shellls of revolution [J]. Acta Aeronautica et Estronautica Sinica, 1992, 13(1): 64-70.

[ 5 ] 曹雷,馬運義,黃玉盈,等.基于raccati傳遞矩陣法分析水下有限長環(huán)肋圓柱殼的聲輻射性能[J].振動與沖擊,2009, 28(9): 149-154.

CAO Lei, MA Yunyi, HUANG Yuying, et al. Analysis of acoustic radiation of ring-stiffened cylindrical shell with finite length in underwater based on riccati transfer matrix method [J]. Journal of Vibration and Shock, 2009, 28(9): 149-154.

[ 6 ] 譚述君, 鐘萬勰.非齊次動力方程duhamel項的精細積分[J].力學學報,2007, 39(3): 374-381.

TAN Shujun, ZHONG Wanxie. Precise integration method for duhamel terms arising from non-homogenous dynamic systems [J]. Chinese Journal of Theoretical and Applied Mechanics, 2007, 39(3): 374-381.

[ 7 ] 陳美霞, 駱東平, 楊叔子.殼間連接形式對雙層殼聲輻射性能的影響[J].振動與沖擊,2005, 24(5): 77-80.

CHEN Meixia, LUO Dongping, YANG Shuzi. Effect of connection type between inner and outer shells on characteristics of sound radiation from finite double cylindrical shell [J]. Journal of Vibration and Shock, 2005, 24(5): 77-80.

[ 8 ] 鐘萬勰.結構動力方程的精細時程積分法[J].大連理工大學學報,1994, 34(2): 131-136.

ZHONG Wanxie. On precise time-integration method for structural dynamics [J]. Journal of Dalian University of Technology, 1994, 34(2): 131-136.

[ 9 ] 陳美霞, 駱東平,曹鋼,等.環(huán)肋柱殼振動和聲輻射試驗分析[J].華中科技大學(自然科學版),2003,31(4): 102-104.

CHEN Meixia, LUO Dongping, CAO Gang, et al. Analysis of vibration and sound radiation from ring-stiffened cylindrical shell [J]. Journal of Huazhong University (Science and Technology), 2003, 31(4): 102-104.

[10] 謝官模, 駱東平.環(huán)肋柱殼在流場中聲輻射性能分析[J].中國造船,1995,131(4): 37-45.

XIE Guanmo, LUO Dongping. Acoustic radiation analysis of ring-stiffened cylindrical shells in fluids[J]. Ship Building of China, 1995, 131(4): 37-45.

附錄A

Animprovedprecisetransfermatrixmethodfortheanalysisofacousticradiationofring-stiffenedcylindricalshells

PANG Fuzhen1,2, WU Chuang3, WANG Xianzhong1, YAO Xiongliang1

(1. College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China; 2. Naval Academy of Armament, Beijing 100161, China; 3. State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, Shanghai 200240, China)

Based on the Flügge shell theory, an improved precise transfer matrix method (IPTMM) was applied to the solution of acoustic radiation for stiffened cylindrical shells. According to the techniques of addition theorem and increment storage, calculation of inhomogeneous was improved. Then, the improved precise transfer matrix method was obtained. The effects of different end boundary condition, structural damping, fluid medium and shell thickness were presented. The results show that when the boundary conditions from freedom supported and simply supported change to rigidly fixed, external radiation acoustic pressure of structure gradually reduces. With the increasing of structural damping, external radiation acoustic pressure of structure gradually reduces.

non-homogeneous item; precise integration transfer matrix; ring-stiffened cylindrical shell; vibration; acoustic radiation

國家自然科學基金項目(51209052);黑龍江省青年科學基金資助項目(QC2011C013);上海交通大學海洋工程國家重點實驗室基金(1307);國防預研項目(4010403010103);中央高校基本科研業(yè)務費資助項目(HEUCF40117);中國博士后基金(2014M552661)

2016-05-17 修改稿收到日期: 2016-09-09

龐福振 男,博士,副教授,1980年生

吳闖 男,博士生,1990年生

U661.44

A

10.13465/j.cnki.jvs.2017.22.021

猜你喜歡
振動結構
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
This “Singing Highway”plays music
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 久久综合九色综合97婷婷| 在线无码av一区二区三区| 伊人久久大香线蕉影院| 国产丰满大乳无码免费播放| 狠狠做深爱婷婷综合一区| 日韩精品毛片人妻AV不卡| 狠狠亚洲五月天| 欧美不卡二区| 国产成年无码AⅤ片在线 | 亚洲色图在线观看| 欧美日韩国产综合视频在线观看| 亚洲第一综合天堂另类专| 999精品视频在线| 无码又爽又刺激的高潮视频| 欧美精品H在线播放| av一区二区无码在线| 欧美笫一页| 国产精品永久久久久| 亚洲无码不卡网| 国产91精品调教在线播放| 亚洲精品国产成人7777| 自拍偷拍欧美日韩| 免费女人18毛片a级毛片视频| 成人亚洲国产| 丁香综合在线| 人妻丰满熟妇αv无码| 全免费a级毛片免费看不卡| 在线观看国产精品一区| 欧美日韩综合网| hezyo加勒比一区二区三区| 国产内射在线观看| 中文字幕伦视频| 国产永久在线视频| 国产精品久久久久久影院| 美女一级毛片无遮挡内谢| 欧美人人干| 波多野结衣国产精品| 丁香亚洲综合五月天婷婷| 国产精品尤物铁牛tv| 久久久噜噜噜| 片在线无码观看| 伊人天堂网| 中文字幕人成人乱码亚洲电影| 日韩欧美视频第一区在线观看| 色综合综合网| 色噜噜狠狠狠综合曰曰曰| 欧美a级完整在线观看| 六月婷婷综合| 91在线精品麻豆欧美在线| 国产综合网站| 毛片手机在线看| 97久久超碰极品视觉盛宴| 欧美另类视频一区二区三区| 久久精品人人做人人综合试看| 97超级碰碰碰碰精品| yjizz国产在线视频网| 天天爽免费视频| 久久精品国产免费观看频道| 亚洲第一成年人网站| 又粗又大又爽又紧免费视频| 国产精品手机视频一区二区| 18禁黄无遮挡免费动漫网站| 亚洲日本中文字幕天堂网| 黄色一及毛片| 精品国产免费观看一区| 欧美成人在线免费| 88av在线看| 午夜欧美理论2019理论| 国产视频大全| 国产精品视频久| 91精品国产91欠久久久久| 久久人妻xunleige无码| 国产欧美专区在线观看| 国产理论一区| 国产一级视频在线观看网站| 免费看的一级毛片| 99ri国产在线| 日本成人在线不卡视频| 一级一级一片免费| 精品亚洲欧美中文字幕在线看| hezyo加勒比一区二区三区| 日韩成人午夜|