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

基于行波法的多板耦合結構振動功率流研究

2017-03-09 07:56:29漆瓊芳
振動與沖擊 2017年3期
關鍵詞:振動結構

漆瓊芳, 喻 敏, 陳 攀

(1.武漢理工大學 交通學院,武漢 430063; 2.中國艦船研究設計中心,武漢 430064)

基于行波法的多板耦合結構振動功率流研究

漆瓊芳1, 喻 敏1, 陳 攀2

(1.武漢理工大學 交通學院,武漢 430063; 2.中國艦船研究設計中心,武漢 430064)

通過引入局部坐標,離散波幅系數矩陣,組裝狀態向量矩陣,建立了離散板的行波法通用表達式,該表達式的提出可以將行波法的研究對象由簡單結構推廣到多板耦合結構。自編行波法MATLAB程序計算結構的穩態響應,將行波法半解析結果與有限元數值結果進行對比,驗證了提出的行波法通用表達式計算多板耦合結構的有效性、高效性。應用行波法分別求解經典薄板理論和Mindlin板理論的振動控制方程,計算耦合結構的主動功率流和被動功率流,結果表明:在全頻范圍內,面內剪切和旋轉慣量對主動功率流和被動功率流有很大影響,在進行功率流主動控制時應盡可能采用厚板理論。損耗因子的增加可以有效減少共振頻率范圍附近的被動功率流峰值,但對其他頻率范圍內的被動功率流峰值影響不大。

行波法;多板耦合結構;主動功率流;被動功率流

目前行波法主要用于研究半無限尺寸結構的振動波傳遞特性,如CREMER[1],計方[2-3],車馳東[4],黃修長[5]等研究了曲梁、平板、線形連接、角形連接、“T”型連接和“十”型連接板的振動傳遞特性,但半無限板不考慮彈性波在邊界條件處的散射和反射,而實際的船舶結構的板為有限尺寸,需要考慮邊界條件對彈性波傳遞的影響,本文行波法研究對象為考慮邊界條件的有限尺寸板。

行波法應用在有限板方面,KESSISSOGLOU[6]應用行波法研究有限尺寸L型板的振動功率流,分析面內波與彎曲波在L型接頭下的轉換。LIU[7-9]應用行波法分析有限L型梁、L型板在簡諧點力激勵下的主動功率流和被動功率流,研究有限尺寸耦合連接板結構的功率流主動控制。薛開等[10]采用改進傅里葉算法,應用半解析解求出了中厚矩形平板的振動功率流。趙芝梅等[11]應用行波法分析力激勵和彎矩激勵對有限L型板振動功率流的影響,探討激勵角度、激勵位置等參數對高頻彎曲波和低頻彎曲波傳遞特性的影響。焦映厚等[12]建立方鋼與L型板的耦合運動模型,用行波法分析阻振方鋼的尺寸和位置對振動功率流的影響。因此,行波法對有限尺寸板的研究對象是任意夾角耦合板和L型板,尚未涉及到T型和十字型板結構。實際船舶結構艙壁與船殼縱橫耦合連接,形成T型、十型等不同的接頭形式,因此研究多板耦合結構的能量傳遞特性具有一定工程意義。同時,文獻對于有限尺寸板振動特性研究大多是基于泊松-克西霍夫(Poisson-Kirchhoff)經典薄板理論的振動控制方程,但在進行高頻動力響應分析時必須考慮面內剪切變形的影響,因此基于厚板理論進行振動傳遞特性分析有其必要性。

本文應用行波法求解基于Mindlin厚板理論的振動控制方程,推導T型板的穩態響應公式,通過三步處理方式:①引入局部坐標;②離散波幅系數矩陣;③組裝狀態向量矩陣。提出了行波法的通用表達式,將行波法的研究對象擴展到十字型板等連接處多板耦合結構。自編行波法MATLAB程序計算耦合結構的穩態響應,將行波法半解析解結果與有限元數值結果進行對比,驗證了行波法通用表達式計算多板耦合結構的有效性、高效性及便捷性。對比厚板和經典薄板的振動控制方程行波解,探討面內剪切、轉動慣量、耦合損耗因子對主動功率流和被動功率流的影響。

1 求解Mindlin板振動控制方程

實際船舶結構的艙壁與船體外板通過連接線耦合,形成典型的多板耦合結構,建立如圖1所示的物理力學模型。在結構的聲學條件改變處將結構離散,即在點激勵處、線耦合處將板離散,離散板的連接線編號分別為0,1,2,…,n,n是離散板數目,離散板編號為板ij(i,j=0,1,2,3,…,n)每個離散板均有一個對偶右手坐標系,如離散板01左端坐標x01y01z01,右端坐標x10y10z10,板01的右端坐標與板12的左端坐標原點重合,y方向一致。每個離散板的位移向量是u,v,w,力向量Nxx、Nxy、Vxy,且分別與每個離散板局部坐標x、y、z一致。

圖1 多板線耦合連接示意圖Fig.1 Schematic diagram of multi-plate with one coupling line

根據Mindlin厚板理論[7],每個離散板的彎曲振動方程為:

(1)

每個離散板的面內振動方程為:

(2)

(3)

式中:D=Eh3/[12(1-μ2)]是板的彎曲剛度;E是板的彈性模量;h是板的厚度;μ為材料泊松比;κ是剪切因子;G是剪切模量;f(x,y,t)是外部點激勵。

任意離散板ij的位移和轉角表達式

(σ2-1)(λ3Ceλ3x+λ4Deλ4x)-

ky(Eeλ5x+Feλ6x)]sin(kyy)eiωt

(5)

(σ2-1)ky(Ceλ3x+Deλ4x)-

(λ5Eeλ5x+λ6Feλ6x)]cos(kyy)eiωt

(6)

任意離散板ij的平面內位移表達式:

kyIeλ9x+kyJeλ10x)sin(kyy)eiωt

(7)

λ9Ieλ9x+λ10Jeλ10x) cos(kyy)eiωt

(8)

式中:w是傅里葉變換后的橫向(z向)位移;u是傅里葉變換后的縱向(x向)位移;v是傅里葉變換子y向位移;φx是繞y軸轉角;φy是繞x軸轉。ky=mπ/Ly是y方向彈性波數,m是模態數;波數

Mindlin厚板的力與位移滿足關系式:

(9)

(10)

(11)

(12)

(13)

(14)

式中:Mxx是彎矩;Mxy是扭矩;Vx是平面外剪力;Nxx是平面內軸力;Nxy是平面內剪力;Γxy是平行于板的剪應變。

2 提出行波法通用表達式

不同于文獻中常規處理方式:將每個離散板的波幅系數和狀態向量表達式分別寫出。本文的三步處理方式:建立每個離散板的局部坐標,組裝離散板的彎曲波的波幅系數矩陣和面內波波幅系數矩陣,分離離散板的位移向量和力向量。處理后提出了行波法通用表達式,將每個離散板ij的狀態向量寫成如下統一形式:

設任意離散板ij的廣義力向量

設離散板ij的廣義位移向量

設離散板ij的彎曲波和面內波的波幅系數矩陣

離散板耦合連接線i的某點(x0,y0)受點力作用,根據級數疊加,點力可以等效為線力

離散板ij在某一模態下的廣義位移向量

(15)

設離散板ij在某一模態下的廣義力向量

(16)

eλ6xeλ7xeλ8xeλ9xeλ10x}

中厚板ij的剛度系數矩陣:

Eh/(1-μ2)Eh/2/(1+μ)}

建立行波法通用表達式后,結構的邊界條件可以一般化寫出:

若板ij左端有點力作用:

(17)

板ij、板jk、板jl、板jm四板線耦合連接,板ij與板jk夾角β1,板ij與板jl夾角β2,板ij與板jm夾角β3:

耦合線連接處的位移連續性條件:

耦合線連接處的力平衡條件

兩個離散板連接夾角為β,坐標系下的位置旋轉矩陣為Tβ

根據實際耦合結構的型式,組裝每個離散板的平衡方程矩陣,形成總體結構的平衡方程矩陣

αX=F

(18)

式中:α是平衡方程矩陣;X是波幅系數矩陣;F是外載荷矩陣。通過求逆即可求出波幅系數矩陣,從而得出離散板的狀態向量,根據Mindlin板理論的功率流表達式[6],Pin是輸入功率流:

(19)

離散板x截面的功率流Ix是:

(20)

功率流的實部是主動功率流Ixa:

(21)

功率流的虛部是被動功率流Ixr:

(22)

3 程序有效性驗證

行波法是基于解析解,是研究振動波傳播機理的一種有效方法,但由于解析解公式推導的復雜性,限制了該方法在復雜結構上的應用。目前,行波法的應用對象是梁、平板、L型板、箱形結構、圓柱殼等簡單結構。本文通過推導T型板的行波解表達式,通過引入坐標轉換矩陣,離散波幅系數矩陣,組裝狀態向量矩陣,推導出了多板耦合結構的動力響應公式,建立了行波解通用的表達式,將行波法的研究對象由簡單結構擴展到多板耦合結構,并自編MATLAB程序進行穩態響應計算,具體流程圖和算法如圖2。

圖2 行波法程序流程圖Fig.2 Flow chart of traveling wave method

3.1T型板有效性驗證

模型為點激勵力下的T型板,根據聲學不連續條件,該模型可以劃分為4個離散板,離散板的板長均為0.4 m,板寬均為0.6 m,板的自由邊界處均為簡單支撐,彈性模量E=2.0×1011Pa,阻尼損耗因子η=0.01,泊松比μ=0.3,復彈性模量取Ep=E(1+jη),密度ρ=7 800 kg/m3,板橫向載荷F=F0eiωt,幅值F0=1 N。通過軟件ABAQUS計算出的有限元數值結果驗證MATLAB行波法程序的正確性,有限元三維模型的網格尺寸為0.02 m×0.02 m,行波法程序是基于Poisson-Kirchhoff薄板理論。

該行波法程序是基于彈性波的模態疊加,應進行收斂性驗證,圖3顯示了激勵點在不同模態疊加數下的橫向位移曲線。可知模態疊加數分別為4和8時,激勵點處橫向位移曲線相差較大,模態疊加數為12時與模態疊加數為8時的曲線較為接近,近似重合,在計算頻率上限3 000 Hz處的橫向位移小于設定的容差1.0×10-10,程序判定收斂。在收斂性域度范圍內,本模型的模態疊加數保守取為20。

圖3 行波法的收斂性驗證Fig.3 Convergence verification of traveling wave method

圖4顯示了10 mm的T型板在在簡諧點激勵下激勵點的橫向響應,可知在計算頻率0~1 500 Hz內,行波法MATLAB程序計算結果與有限元數值計算結果吻合較好,計算頻率大于1 500 Hz時,兩者的位移響應峰值會有偏移,這與由于有限元方法的高頻不穩定有關。這說明,基于行波法通用表達式的MATLAB程序可以有效準確的計算T型板的振動響應。

圖4 T型板激勵點橫向位移Fig.4 Z-direction displacement of T-shaped plate at excited point

3.2 十字型板有效性驗證

點激勵力的十字型可以劃分為5個離散板,離散板尺寸、材料屬性和邊界條件與3.1節一致。圖4顯示了十字型板繞y軸轉角的頻響曲線,在計算頻率0~1 200 Hz范圍內,頻響曲線吻合較好,驗證了行波法通用表達式及自編行波法程序的有效性。

圖5 十字型板激勵點繞y軸轉角Fig.5 Y-axis angle of cross-shaped plate at excited point

3.3 本文行波法算法的高效性驗證

一般采用有限元法進行船舶結構振動計算,并且有限元法已有成熟的商業軟件ABAQUS、ANSYS、NASTRAN/PATRAN等,但有限元數值計算方法有一定局限性,隨著計算頻率的增加,網格密度大于彎曲波長度的1/6,有限元網格密度必須隨之增加,這增加存儲內存,增加了有限元計算時間。且有限元法具有高頻不穩定性,進行高頻動力響應分析時會影響計算精度和準確度。此外,有限元是一種數值解法,對振動波的機理研究方面沒有優勢。行波法通用表達式的提出使行波法求解復雜結構成為可能,而且行波法代表波的疊加,具有明確的物理意義,是振動傳遞特性的機理性研究的重要方法。為驗證基于行波法的通用表達式的MATLAB程序的計算高效性,對船舶典型結構進行算例分析,MATLAB程序的行波法計算和ABAQUS有限元計算均在同一臺筆記本(內存4 GB,雙核2.5 GHz)上進行,結構模型的每個離散板的尺寸均為0.6 m×0.6 m,網格尺寸為0.02 mm ×0.02 mm,物理屬性與3.1一致,計算頻率上下限為1~1 000 Hz,步長為2 Hz,模型邊界都是簡支,選擇Poisson-Kirchhoff薄板理論模型。表1列出了兩種方法的計算時間,對比可知行波法耗費的時間比有限元法要少。

表1 行波法與有限元法的計算時間

行波法通用表達式有以下特點:①求出的是基于行波法和模態疊加法的半解析解,對船體板構件進行振動計算時,只需要改變結構的材料參數、截面屬性、邊界條件、局部坐標、激勵屬性等就可以對結構的振動響應進行快速計算,并給出結構的位移響應曲線和功率流曲線。②與有限元法需要密集的網格單元不同,此程序對一個離散板只建立一個單元,只占據一個單元的存儲空間,存儲矩陣小,具有計算快速性的優點。若計算模型改變,只需在界面里改變模型物理參數即可對模型進行再次計算,相比有限元法減少了三維建模和網格劃分的工作量。

4 算例分析

4.1 Kirchhoff板理論和Mindlin板理論對比

圖6 T型板厚為3 mm時的橫向位移Fig.6 Z-direction displacement of 3 mm T-shaped plate

圖7 T型板厚為10 mm時的橫向位移Fig.7 Z-direction displacement of 10 mm T-shaped plate

圖8 T型板厚為17 mm時的橫向位移Fig.8 Z-direction displacement of 17 mm T-shaped plate

當T型板對邊為簡單支撐,其余邊自由時。圖6、圖7、圖8分別顯示板厚為3 mm、13 m、17 mm時,兩種板理論計算的點激勵處板的橫向位移,對比可知,在500 Hz頻率以下,兩種板理論的振動控制方程推出的行波解吻合較好,大于該頻率范圍則會有差異,且隨著板厚增加,差別越大,但總體變化趨勢不大,在相位上會有錯開。這是因為,當計算頻率一定時,板厚越小,板面內的剪切慣量和轉動慣量對振動影響越小。當板厚一定時,隨著計算頻率增大,面內剪切和轉動慣量對振動響應影響越大。板厚越小,位移響應曲線在共振頻率下的峰值和低谷越密集,在計算頻域范圍內模態數越多。

圖9 T型板厚為10 mm時的主動功率流Fig.9 Active power flow of T-shaped plate with thickness in 10 mm

圖10 T型板厚為10 mm時的被動功率流Fig.10 Reactive power flow of 10 mm T-shaped plate

T型板厚10 mm時,參考功率級10-12W,圖9、圖10分別給出了基于厚、薄板兩種理論的振動控制方程的主動功率流和被動功率流行波解。對比可知,在全頻率范圍內,兩種板理論計算的主動功率流和被動功率流均有差異,在頻率1 000 Hz以上,差別明顯,因此在進行振動主動控制時[9],若以薄板理論的控制方程為行波法基礎,以主動功率流為優化目標,則會造成很大的誤差,達不到主動控制效果,因此在進行主動控制時應在全頻域范圍內盡可能用厚板理論進行修正。

4.2 損耗因子對功率流的影響

T型板厚為10 mm時,應用Poisson-Kirchhoff薄板理論,模型尺寸與屬性與3.1節一致,計算頻率范圍取100~1 000 Hz。探討損耗因子變化對主動功率流和被動功率流的影響。

圖11 板損耗因子變化時的主動功率Fig.11 Active power flow with different loss factor

圖12 板損耗因子變化時的被動功率流Fig.12 Reactive power flow with different loss factor

圖11顯示了損耗因子變化時,T型板的主動功率流,對比可知,無損耗因子時的主動功率流比有損耗因子時的主動功率流小得多,損耗因子存在與否對主動功率流的影響顯著。當損耗因子不存在時,在計算頻率范圍內,主動功率流在-50 dB上下范圍內。隨著損耗因子的增加,主動功率流會增加,損耗因子沒有改變主動功率流的共振峰對應的頻率。這說明主動功率流與結構阻尼密切相關。當損耗因子為零,結構阻尼不能使振動能量損耗,主動功率流很小,可以被忽略。當損耗因子非零,振動能量可以通過振動傳遞,主動功率流不為零。圖12可知,損耗因子的變化對被動功率流曲線影響不大,且損耗因子的增加可以減少被動功率流在共振頻率附近的峰值,對其他頻率范圍內的被動功率流影響不大。

5 結 論

應用行波法分別求解基于Poisson-Kirchhoff薄板和Mindlin厚板理論的振動控制方程,推導T型板耦合結構的穩態響應公式,通過三步處理方式:引入局部坐標,離散波幅系數矩陣,組裝狀態向量矩陣,建立行波法通用表達式,自編行波法MATLAB程序計算結構的穩態響應,將行波法半解析解結果與有限元數值結果進行對比驗證,根據厚板理論和薄板理論分析耦合結構的主動功率流和被動功率流,總結如下:

(1)提出的通用表達式和通用計算程序可以有效地計算多塊板耦合結構的振動響應,并且計算效率較高,該思路將行波法的研究對象擴展到十字型板等連接處多塊板耦合結構。

(2)在低頻范圍內,面內剪切和旋轉慣量對穩態響應影響不大。在全頻范圍內,面內剪切和旋轉慣量對主動功率流和被動功率流有很大影響,因此在進行功率流主動控制是要適當對薄板理論進行修正。

(3)損耗因子是否存在對主動功率流影響顯著,無損耗因子時的主動功率流比有損耗因子時的主動功率流小得多。損耗因子對被動功率流曲線影響不大,損耗因子的增加可以有效減少共振頻率范圍附近的被動功率流峰值,但對其他頻率范圍內的被動功率流峰值影響不大。

[ 1 ] CREMER L, HECKL M, PETERSSON B A T. Structure-borne sound [J]. Physics Today, 1975, (28): 81

[ 2 ] JI Fang, YAO Xiongliang, YE Xi.The influence of blocking mass parameters on the vibration isolation performance of a power cabin[J]. Journal of Marine Science and Application, 2011, 10(1):25-32.

[ 3 ] 計方, 路曉東, 姚熊亮. 船體結構黏彈性夾層阻抑振動波傳遞特性研究[J]. 應用基礎與工程科學學報, 2012, 20(3): 464-471. JI Fang, LU Xiaodong, YAO Xiongliang.Research on hull structure elastic interlayer impending vibration wave propagation[J]. Journal of Basic Science and Engineering, 2012, 20(3): 464-471.

[ 4 ] 車馳東, 陳端石. 多轉角結構中轉角處附加的阻振質量對結構波傳遞的抑制[J]. 船舶力學, 2010, 14(9): 1052-1064. CHE Chidong, CHEN Duanshi.Structure-borne sound attenuation in a multi-corner structure with attached blocking masses[J]. Jounal of Ship Mechanics, 2010, 14(9): 1052-1064.

[ 5 ] 黃修長,徐時吟,章振華,等.利用波動法研究曲梁結構中的波形轉換和能量傳遞[J].振動與沖擊,2013,31(8):38-46. HUANG Xiuchang, XU Shiyin, ZHANG Zhenghua, et al. Application of wave approach in wave mode conversion and energy transmission of curved beams [J].Journal of Vibration and Shock,2013,31(8):38-46.

[ 6 ] KESSISSOGLOU N J. Power transmission in L-shaped plates including flexural and in-plane vibration[J]. The Journal of the Acoustical Society of America, 2004, 115(3): 1157-1169.

[ 7 ] LIU C C, LI F M, LIANG T W, et al. The wave and vibratory power transmission in a finite L-shaped Mindlin plate with two simply supported opposite edges[J]. Acta Mechanica Sinica, 2011, 27(5): 785-795.

[ 8 ] LIU C C, LI F, HUANG W. Active vibration control of finite L-shaped beam with travelling wave approach[J]. Acta Mechanica Solida Sinica, 2010, 23(5): 377-385.

[ 9 ] LIU C C, LI F M, FANG B, et al. Active control of power flow transmission in finite connected plate[J]. Journal of Sound and Vibration, 2010, 329(20): 4124-4135.

[10] 趙芝梅, 盛美萍. 激勵特性對L型板振動功率流的影響[J].兵工學報, 2013,34(8):986-993. ZHAO Zhimei, SHENG Meiping.Effect of driving force characteristics on the power flow of L-shaped plate[J]. Acta Armamentarii,2013,34(8):986-993.

[11] 薛開, 王久法, 李秋紅, 等. 中厚矩形板的振動功率流特性分析[J]. 振動與沖擊, 2013,32(21):78-181. XUE Kai, WANG Jiufa, LI Qiuhong, et al.Vibration power flow analysis of moderately thick rectangular plates[J].Journal of Vibration and Shock,2013,32(21):78-181.

[12] 焦映厚, 侯守武, 陳照波, 等. 采用行波法研究有限尺寸加肋L型板結構的振動特性[J]. 振動工程學報, 2013,26(6):871-878. JIAO Yinghou, HOU Shouwu, CHEN Zhaobo, et al. Research on vibration properties of finite rib-stiffened L-shaped plateusing wave method[J].Journal of Vibration Engineering, 2013,26(6):871-878.

Vibration power flow of multi-plate coupled structures based on traveling wave approach

QI Qiongfang1, YU Min1, CHEN Pan2

(1. School of Transportation, Wuhan University of Technology, Wuhan 430063, China;2. China Ship Development and Design Center, Wuhan 430064, China)

The general expressions of the traveling wave approach for a discrete plate were established, by introducing local coordinates, discrete wave amplitude coefficient matrix, and assemble state vector matrix. These general expressions led to that the study object of the traveling wave approach was extended from simple plates to multi-plate coupled structures. Structural steady vibration responses calculated with the self-compiled program of the traveling wave approach were compared with the numerical results with the finite element method. This traveling wave approach semi-analytical solution was verified to be effective and convenient. Based on the classical thin plate theory and Mindlin plate theory, respectively, the traveling wave method was used to compute active power flow and passive power flow of multi-plate coupled structures. The calculated results showed that within the whole frequency range, shear and rotary inertia inplane have a great influence on active power flow and passive power flow; when conducting active power flow control, using the thick plate theory is recomended as far as possible; increase in loss factor can effectively reduce the peak value of passive power flow near the frequency range of resonances, but this influence is not large within the other frequency range.

traveling wave approach; multi-plate coupled structure; active power flow; passive power flow; reactive power flow

航空科學基金(20142365002)

2015-10-09 修改稿收到日期:2016-01-03

漆瓊芳 女,碩士生,1989年11月生

喻敏 女,博士,副教授,1978年11月生

TP533

A

10.13465/j.cnki.jvs.2017.03.025

猜你喜歡
振動結構
振動的思考
科學大眾(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
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 国产成年女人特黄特色大片免费| 亚洲第一区精品日韩在线播放| 波多野结衣一区二区三区88| 国产福利微拍精品一区二区| 亚洲乱码视频| 亚洲欧美精品在线| 国产美女91视频| 无码区日韩专区免费系列| 精品国产美女福到在线直播| 国产一区二区三区精品久久呦| 99精品影院| 国产免费a级片| 中国一级特黄大片在线观看| 国产成人凹凸视频在线| 99国产在线视频| 伊大人香蕉久久网欧美| 久久久久久尹人网香蕉| 国产精品自在自线免费观看| 久久精品女人天堂aaa| 一区二区理伦视频| 欧美成a人片在线观看| 人人爽人人爽人人片| 国产高清在线精品一区二区三区| 老司国产精品视频91| 日韩毛片免费视频| 亚洲欧美不卡中文字幕| 精品第一国产综合精品Aⅴ| 中文字幕第4页| 一级毛片免费播放视频| 国产午夜人做人免费视频中文 | 国产成a人片在线播放| 国产精品太粉嫩高中在线观看 | 国产黄色爱视频| 久久大香伊蕉在人线观看热2| 青青草国产一区二区三区| 福利姬国产精品一区在线| 久久久久亚洲Av片无码观看| 国产aaaaa一级毛片| 伊人网址在线| 亚洲一区二区无码视频| 全部免费毛片免费播放| 污网站在线观看视频| 欧美一级在线看| 四虎永久在线精品影院| 国产亚洲视频免费播放| 国产99视频精品免费观看9e| 久久不卡精品| 国产内射一区亚洲| 欧美一区二区人人喊爽| 视频一本大道香蕉久在线播放| 一本色道久久88| 亚洲Aⅴ无码专区在线观看q| 国产一级特黄aa级特黄裸毛片 | 成人国产免费| 国产视频入口| 蜜臀av性久久久久蜜臀aⅴ麻豆 | 国产一级小视频| 国产情侣一区| 国产波多野结衣中文在线播放| 亚洲a级毛片| 99精品国产电影| 亚洲男人的天堂网| 欧美日韩一区二区在线播放| 亚洲va欧美va国产综合下载| 久久精品女人天堂aaa| 直接黄91麻豆网站| 小蝌蚪亚洲精品国产| 欧美成人第一页| 国产精品无码在线看| 国产成人亚洲精品色欲AV| 亚洲精品在线91| 亚洲欧美日韩中文字幕在线| 国产成人亚洲综合A∨在线播放| 亚洲an第二区国产精品| 91精品专区国产盗摄| 9999在线视频| 白浆免费视频国产精品视频| 国产AV毛片| 中文字幕乱妇无码AV在线| 91精品久久久久久无码人妻| 亚洲系列中文字幕一区二区| 欧美国产日韩在线|