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

基于波傳播法的橢圓柱殼自由振動特性研究

2017-06-19 19:35:26張冠軍李天勻繆宇躍
振動與沖擊 2017年12期
關鍵詞:模態振動

張冠軍, 朱 翔, 李天勻, 繆宇躍

(華中科技大學 船舶與海洋工程學院, 武漢 430074)

基于波傳播法的橢圓柱殼自由振動特性研究

張冠軍, 朱 翔, 李天勻, 繆宇躍

(華中科技大學 船舶與海洋工程學院, 武漢 430074)

基于Flügge殼體理論推導出橢圓柱殼的自由振動方程。由于橢圓柱殼周向曲率的變化,造成振動方程在周向模態階數域內不解耦,采用波傳播法將殼體位移以雙Fourier級數形式展開,橢圓截面曲率半徑以單Fourier級數形式展開,通過級數變換,將變系數的偏微分方程組轉換為關于周向模態階數相互耦合的有限階常系數線性方程組,通過求解耦合振動方程得到橢圓柱殼的固有頻率。隨后對影響殼體固有頻率的主要參數進行分析,得到橢圓柱殼對稱和反對稱模態頻率隨橢圓度、殼長比等參數的變化規律。

橢圓柱殼;自由振動;固有頻率;橢圓度參數;反對稱模態

隨著加工技術的進步,橢圓柱殼在潛器、航空航天、管道輸送等領域的應用日益增多。相比于圓柱殼結構,橢圓柱殼在潛器、航空器倉容及設備布置上具有更大優勢,橢圓形管道相對于方形管道斷面濕周面積小,密閉性更好,相對于圓形管道占用高度空間小。另一方面,即使是針對原設計的圓柱殼結構,由于制造工藝、加工方法或焊接變形,乃至巨大的深水壓力作用等諸多因素的影響,其橫截面也無法保證為一個完美的圓形,不可避免的存在各種幾何缺陷,可能產生不可忽略的橢圓度等偏差[1-2]。因此,研究橢圓柱殼的自由振動特性具有重要意義。

目前在結構工程領域對于圓柱殼體的振動特性研究較為全面,其相關文獻也較多,而對于橢圓柱殼振動方面的研究還相對較少。由于橢圓柱殼周向曲率的改變導致其振動方程是一個變系數的偏微分方程組,因此無法像圓柱殼那樣可以得到解析解,只能借助于各種數值解法進行求解,而其方程又頗為繁瑣,給求解帶來了諸多不便。Klosner[3-4]最早運用攝動解法對橢圓柱殼自由振動特性進行了研究,但其方法僅適用于橢圓度較小的情況,應用范圍受到限制。隨后Culberson和Boyd等[5-6]利用部分簡化后的Love理論和Donnel理論研究了橢圓柱殼的自由振動,但簡化使得求解殼體低頻頻率及計算較長的殼體頻率時,其結果存在較大誤差。Elsbernd等[7]采用變分法對具有初始應力的橢圓柱殼進行了研究,其方程形式頗為繁瑣復雜,應用起來并不方便,且其只對特定的橢圓度進行了計算,沒有深入討論固有頻率與橢圓度之間的關系。朱建雄等[8]使用了一種特殊的坐標變換將復雜截面形狀的柱殼映射成圓柱殼,通過對變化后的圓柱殼進行相關的振動分析可以近似求得非圓截面柱殼在不同邊界條件下的頻率及振型,但該方法中結果的準確性依賴于映射形狀函數節點數的選取,而節點數越多則公式推導及計算越復雜,且該方法涉及大量導數求積及高斯數值積分,使用上存在一定局限性。最近,Khalifa[9]基于Flügge殼體理論,采用傳遞矩陣法研究了支撐在彈性地基上的正交各向異性橢圓柱殼的自由振動特性,討論了殼體及地基相關參數對橢圓柱殼振動特性的影響。Ahmed等[10]同樣采用傳遞矩陣法研究了周向變厚度橢圓柱殼的自由振動特性。但傳遞矩陣法受單元傳遞矩陣的計算精度、傳遞矩陣連乘過程中的累積誤差以及計算機的舍入(或截斷)誤差影響,狀態向量在傳遞過程中會產生一定的精度損失,對于指數矩陣還存在收斂穩定性問題,且計算效率較低[11]。最近張盛等用有限元軟件對包括具有初始橢圓度在內的幾種含有幾何制造誤差的圓柱殼的聲振特性進行了仿真分析和計算。Tornabene等[12]采用一種和有限元法類似的廣義積分法求解了復合材料橢圓柱殼的自由振動特性。Namazi等[13]采用有限元法與實驗法研究了復合材料橢球殼的自由振動特性,并將有限元法結果與實驗結果進行對比,結果吻合較好,同時研究了材料參數及結構參數對橢球殼自由振動的影響。

本文基于Flügge[14]的殼體理論,采用波傳播法并結合Marguerre[15]對非圓截面曲率半徑的處理方法對殼體的自由振動方程進行了推導,將變系數的偏微分方程組轉換為周向模態階數相互耦合的有限階常系數線性方程組進行求解,得到橢圓柱殼的固有頻率及振型。通過與已有文獻結果對比,證明了本文處理橢圓柱殼周向變曲率的方法以及理論建模和求解方法的有效性。其后進行了更加深入的研究,分析討論了橢圓度和殼長L/mr0比對殼體對稱和反對稱模態頻率的影響。

1 研究對象

圖1 橢圓柱殼幾何參數及坐標系Fig.1 Geometry parameters and coordinate system of elliptic cylindrical shell

現引入下列無量綱坐標:

(1)

2 基本理論推導

2.1 殼體的位移及內力基本關系

根據Flügge薄殼理論,殼體的幾何方程為:

(2)

式中:εx和εs為殼體中面內各點的線應變;γxs為剪應變;κx,κs和τ為中面內各點主曲率及扭率的改變。

殼體的物理方程為:

(3)

式中:K=Eh/(1-μ2)為殼體的拉壓剛度;D=Eh3/12(1-μ2)為殼體的彎曲剛度;Nx,Ns分別為x和s方向單位寬度上的面內力;Nsx和Nxs為平面內單位寬度上的剪切力;Mx,Ms和Mxs,Msx分別為單位寬度上的彎矩和扭矩。

殼體內的靜力平衡方程為:

(4)

式中:Qx和Qs分別為x和s方向單位寬度上的橫剪力。

將式(4)中第4式和第5式代入到第2式和第3式并聯立第1式,可得到一組新的平衡方程:

(5)

將式(2),式(3)代入式(5)中,可得到矩陣形式表達的關于位移變量的平衡方程:

(6)

其中各算子表達式如下:

2.2 真空中橢圓柱殼體的自由振動方程

為進一步求解式(6),需將Lij(i,j=1,2,3)轉變為常系數的偏微分形式。首先需要給出橢圓柱殼截面的曲率半徑r(s)的具體表達式。Marguerre首先提出了用傅里葉級數表示非圓截面曲率半徑的方法。隨后,Romano等[16]將Marguerre的方法進行了簡化,只取其傅里葉級數的前兩項來表示橢圓截面的曲率半徑,其無量綱表達式如下:

(7)

ε=3Q-36/35Q3

(8)

圓柱殼在周向有無數條對稱軸,其對稱模態和反對稱模態的振型相同,可偏轉一定角度重合。對于橢圓柱殼,只有長軸和短軸兩條對稱軸,有必要分別討論對稱與反對稱模態的變化規律(本文中所討論橢圓柱殼的對稱與反對稱模態均相對于長軸)。將殼體各位移在軸向和周向展開成雙Fourier級數的形式。對于對稱模態,其展開形式為

(9)

對于反對稱模態,其展開形式為

(10)

式中:m是軸向半波數;n是周向模態階數;Um,n,Vm,n,Wm,n為不同m,n下殼體軸向、周向和徑向位移Fourier幅值系數;km為殼體的軸向半波數,不同邊界條件下km的可根據梁函數推導出其近似值[17]。需要指出的是基于波傳播法的位移展開形式對于簡支邊界條件,長殼短殼均適用,但其它邊界條件則主要適用于中長殼。如果在殼體軸向位移函數假設形式上采用Fourier級數擴展形式,則對于短殼也可以擴展到任意邊界。

將式(7)和式(9)代入到平衡方程式(6)中,通過級數變換可將式(6)中的變系數轉變為常系數,得到三個關于Um,n,Vm,n,Wm,n沿周向模態階數n相互耦合新的方程

(11)

3 數值計算及分析

進一步考察式(11)可知,式(11)中的Um,n,Vm,n,Wm,n關于不同周向模態階數n是相互耦合的,且n均為奇數或偶數。表明式(11)中,n的奇數項與偶數項是非耦合的,則Um,n,Vm,n,Wm,n可以根據n的奇偶數性分開求解。對于任意給定的軸向半波數m,Um,n,Vm,n,Wm,n都可以按n的奇偶值分成兩組,由于橢圓度ε的影響,每一組中Um,n,Vm,n,Wm,n又會相互耦合成無窮多個線性方程,需對方程組進行截斷求解,截取的項數與解的精確度有關。

當n取有限項p時,可以得到3p個線形方程(m>0),并寫成矩陣形式:

[M(q)-Ω2I]X(q)=0(q=1,2)

(12)

式中:X(1)[Um,1Um,3…Um,2p-1?Vm,1Vm,3…Vm,2p-1?Wm,1Wm,3…Wm,2p-1],對應于n取奇數;X(2)[Um,0Um,2…Um,2p-2?Vm,2Vm,4…Vm,2p-2?Wm,0Wm,2…Wm,2p-2],對應于n取偶數。

M(q)(q=1,2)為3p階的實數方陣,由式(11)中的系數循環迭代而成,I為單位矩陣,0表示零向量。

為滿足式(12)中的X(q)有非零解,則要求:

(13)

求解式(13)可得橢圓柱殼的固有頻率。

3.1 收斂性分析

這里選取周向模態階數n=10進行分析,利用QR分解法求解式(13),得到真空中橢圓柱殼的自振頻率。計算的模型參數:殼厚比h/r0=1/100,泊松比μ=0.3,楊氏模量E=2.1×1011Pa,材料密度ρs=7 800 kg/m3,殼體兩端簡支。在后文關于不同橢圓度橢圓柱殼的對比分析及橢圓柱殼與圓柱殼的對比中,殼體的平均半徑r0是相等的,即殼體的截面周長是相等的。橢圓柱殼固有頻率隨截斷項數變化的計算結果,如圖 2 所示。

圖2 無量綱固有頻率隨截斷項數的變化(n=1)Fig.2 Dimensionless natural frequency changes with the truncated number(n=1)

從圖2可知,隨著截斷項數p的增加,橢圓柱殼固有頻率最終會趨于穩定,即表明計算結果趨于收斂。當殼體橢圓度較小時,計算結果趨于收斂所需的截斷項數相對較少,當橢圓度較大時,則需更高的截斷項數以保證結果收斂。隨著截斷項數的增加,計算結果的精度也會更高,但過多的截斷項數會導致計算效率的降低,在n=10,L/mr0≤20范圍內,取截斷項數p=25可以保證足夠的精度。通過計算發現,在其它周向模態階數(包括奇數階數和偶數階數)下的收斂性情況都也都大體一致,因此在本文后續的計算中取p=25。

3.2 數值方法及驗證

為了驗證本文計算方法及程序的可靠性,同基于變分原理進行求解橢圓柱殼固有頻率的文獻[7]結果進行了對比,如表1所示。從表1可知,兩者結果吻合很好,證明了本方法的準確性與可靠性。

表1 橢圓柱殼的無量綱固有頻率Ω對比驗證 (h/r0=0.002,μ=0.3,ε=0.5,E=2.1×1011,兩端簡支)Tab.1 Dimensionless natural frequency verification and comparison of elliptic cylindrical shell (h/r0=0.002,μ=0.3,ε=0.5,E=2.1×1011, Simply supported-Simply supported)

3.3 橢圓柱殼模態振型

圓柱殼有無數條對稱軸,且對稱和反對稱振型可

通過偏轉重合,而橢圓柱殼只有長軸和短軸兩條對稱軸,其對稱振型和反對稱振型不同,因此,有必要對橢圓柱殼的對稱振型和反對稱振型進行研究。

圖3給出了不同周向模態階數下橢圓柱殼的橫截面振型圖,可以看出除周向模態階數n=0以外,其他周向模態階數均存在振型明顯不同的對稱模態和反對稱模態。

圖3 橢圓柱殼模態振型圖Fig.3 Mode shape of elliptic cylindrical shell

3.4 橢圓柱殼固有頻率隨殼體參數的變化

3.4.1 殼體固有頻率隨橢圓度的變化

橢圓度是衡量橢圓扁平程度的主要物理參數,同時殼體長度的不同也將影響殼體的固有頻率,因此,研究橢圓度及殼體長度對橢圓柱殼固有頻率的影響很有必要。圖4給出了幾種不同殼長比L/mr0下各周向模態階數n對應的橢圓殼的無量綱固有頻率Ω隨橢圓度ε的變化曲線。模型參數:ρs=7 800 kg/m3,E=2.1×1011Pa,μ=0.3,h/r0=1/100,殼體兩端簡支。

“——”對稱模態頻率,“……”反對稱模態頻率圖4 固有頻率隨橢圓度的變化關系 (m=1,兩端簡支)Fig.4 The relationship between the natural frequency with ellipticity (m=1,simply supported-simply supported)

由圖4可知,橢圓柱殼的固有頻率相對于圓柱殼有明顯的不同,對應于同一周向模態階數n,橢圓柱殼的對稱和反對稱模態頻率值不相等,兩頻率的差值隨橢圓度的增大而變大。當n=1,3,5等奇數時,橢圓柱殼對稱模態頻率大于反對稱模態頻率,當n=2,4等偶數時,反對稱模態頻率大于對稱模態頻率,呈現出相反的變化規律。殼體殼長比L/mr0較小時,對應于各周向模態階數,對稱和反對稱模態頻率均差異較大。當L/mr0較大時,除了n=1以外,其他周向模態階數對應的對稱和反對稱模態頻率差別不大。當L/mr0=20時,只有n=1對應的對稱與反對稱模態頻率差異比較明顯,其它周向模態階數n對應的對稱與反對稱模態隨頻率逐漸趨于重合。

殼長比L/mr0不同,橢圓柱殼基頻所對應的周向模態階數n也不同。當L/mr0>10時,橢圓柱殼的基頻對應于周向模態階數n=2,L/mr0不斷減小時,則橢圓柱殼基頻對應的周向模態階數n不斷增大。其它周向模態階數n越接近于基頻所對應的周向模態階數n,則其對應的頻率也越低,周向模態階數n=0所對應的頻率一般較高。

在實際工程中,對于柱殼結構,一般都按照理想圓柱殼模型采取減振降噪措施,當存在橢圓度時,殼體固有頻率將發生改變。但若能控制殼體的橢圓度,使橢圓柱殼和理想圓柱殼的固有頻率差值在一定范圍內(如<5%),則可近似按照理想圓柱殼模型采取減振降噪措施。圖5 給出了L/mr0=3和L/mr0=20時,不同周向模態階數對應的殼體相對頻率比值Ωε/Ωc隨ε的變化曲線,其中Ωε表示不同橢圓度下橢圓柱殼無量綱固有頻率,Ωc表示半徑為r0的圓柱殼的無量綱固有頻率,即ε=0時的無量綱頻率。

從圖5可知,不同周向模態階數n對應的殼體相對頻率比值隨橢圓度的變化是不同的,殼長比不同時其變化規律也不相同,但若能將殼體橢圓度控制在ε≤0.2的范圍內,各周向模態階數對應的固有頻率相對比值差值均<5%,則可以近似將橢圓柱殼視為理想圓柱殼開展減振降噪工作。

圖5 殼體相對頻率比值隨橢圓度的變化(m=1,兩端簡支)Fig.5 Relative frequency ratio varies with ellipticity (m=1,simply supported-simply supported)

3.4.2 固有頻率隨殼長比的變化

圖6給出了在給定橢圓度ε的情況下,各周向模態階數n對應的橢圓殼體的固有頻率Ω和殼長比L/r0的變化曲線。

圖6 橢圓殼體固有頻率和殼長比之間的關系(m=1,兩端簡支)Fig.6 The relationship between the natural frequency with length ratio (m=1, simply supported-simply supported)

從圖6可知:當橢圓度一定時,橢圓柱殼各周向模態階數n對應的固有頻率隨殼長比的變化都較為相似,在殼長比較小時,無量綱頻率下降較快,當殼長比逐漸增大時,其下降的速率逐漸減小并趨于平緩。表明在給定殼體半徑r0以及軸向半波數m情況下,當殼體長度L較短時,殼體固有頻率對長度的變化較敏感,長度的改變會引起殼體固有頻率的較大變化,但當殼體較長時,殼體長度的改變對結構固有頻率的影響較小,與圓柱殼表現出類似的規律。文中只給出了對稱模態頻率的變化規律,反對稱模態與之類似。

4 結 論

本文基于Flügge的薄殼理論推導出真空中橢圓柱殼的自由振動方程,采用波傳播法將殼體位移以雙Fourier級數形式展開,周向曲率半徑以單Fourier級數形式展開,通過級數變換將變系數的偏微分方程組轉換為關于周向模態階數相互耦合的有限階常系數線性方程組,并由此求解出橢圓柱殼在不同參數下的自振頻率,得出如下結論:

(1)橢圓柱殼關于長軸和短軸對稱,且對稱和反對稱模態頻率及振型不同。

(2)對應于同一周向模態階數n,橢圓柱殼的對稱和反對稱模態頻率值不相等,兩頻率的差值隨橢圓度的增大而變大。殼長比較小時,對應于各周向模態階數,對稱和反對稱模態頻率均差異較大;殼長比較大時,除了n=1以外,其他周向模態階數對應的對稱和反對稱模態頻率差別不大。

(3)殼體固有頻率隨著殼體長度的增加而減小,且減小的速率隨著殼體長度的增加而降低。

(4)在實際工程中,若能將橢圓度控制在ε≤2.2的范圍內,則可以近似將橢圓柱殼視為理想圓柱殼展開減振降噪工作。

[1] 龔有根, 賀玲鳳. 含有初始凹陷缺陷圓柱殼穩定承載能力的實驗研究與數值計算[J]. 實驗力學, 2010, 25(1): 73-80.

GONG Yougen, HE Lingfeng. Experimental study and numerical calculation of stability and load-carrying capacity of cylindrical shell with initial dent[J]. Journal of Experimental Mechanics, 2010, 25(1): 73-80.

[2] 張盛, 金翔, 周樺. 加肋圓柱殼制造誤差對聲學性能的影響研究[J]. 中國艦船研究, 2011, 6(4): 43-50.

ZHANG Sheng, JIN Xiang, ZHOU Hua. Influence of construction error on sound radiation for ring-stiffened cylindrical shell[J]. Chinese Journal of Ship Research, 2011, 6(4): 43-50.

[3] KLOSNER J M, POHLE F V. Natural frequencies of an infinitely long noncircular cylindrical shell[R]. PIBAL Rept, 476, 1958.

[4] KLOSNER J M. Frequencies of an infinitely long noncircular cylindrical shell[R]. PIBAL Rept. 552, 1959.

[5] CULBERSON L D, BOYD D E. Free vibrations of freely supported oval cylinders[J]. AIAA Journal, 1971, 9(8): 1474-1480.

[6] CARL E K, BOYD D E. Free vibrations of noncircular cylindrical shell segments[J]. AIAA Journal, 1971, 9(2): 239-244.

[7] ELSBERND G F, LEISSA A W. The vibrations of non-circular cylindrical shells with initial stresses[J]. Journal of Sound and Vibration, 1973, 29(3): 309-329.

[8] 朱建雄,曹志遠,李國豪. 非圓柱殼在各種邊界條件下的自由振動分析[J]. 力學學報. 1992, 24(2): 171-179.

ZHU Jianxiong, CAO Zhiyuan, LI Guohao. Free vibration of noncircular cylindrical shells with arbitrary boundary conditions[J]. Journal of Theoretical and Applied Mechanics, 1992, 24(2): 171-179.

[9] KHALIFA M. Effects of non-uniform Winkler foundation and non-homogeneity on the free vibration of an orthotropic elliptical cylindrical shell[J]. European Journal of Mechanics-A/Solids, 2015, 49: 570-581.

[10] AHMED, KHALIFA M. Simplified equations and solutions for the free vibration of an orthotropic oval cylindrical shell with variable thickness[J]. Mathematical Methods in the Applied Sciences, 2011, 34(14):1789-1800.

[11] 曹雷, 馬運義, 黃玉盈. 環肋加強變厚度圓柱殼的自由振動[J]. 華中科技大學學報(城市科學版), 2007, 24(2): 63-66.

CAO Lei, MA Yunyi, HUANG Yuying. Free vibration of ring-stiffened circular cylindrical shell with variable thickness[J]. Journal of Huazhong University of Science and Technology(Urban Science Edition), 2007, 24(2): 63-66.

[12] TORNABENE F, FANTUZZI N, BACCIOCCHI M, et al. Free vibrations of composite oval and elliptic cylinders by the generalized differential quadrature method[J]. Thin-Walled Structures, 2015, 97: 114-129.

[13] NAMAZI M M, AGHANAJAFI C. The analyze of vibration of composite elliptical shell[J]. Australian Journal of Basic and Applied Sciences, 2010(7):1542-1554.

[14] FLüGGE W. Stress in shells[M]. 2nd ed. Berlin and New York: Spring-Verlag, 1973.

[15] MARGUERRE K. Stability of the cylindrical shells of variable curvature[R]. NASA TM 1302, 1951.

[16] ROMANO F, KEMPNER J. Stresses in short noncircular cylindrical shells under lateral pressure[J]. Journal of Applied Mechanics, 1962, 29(4): 669-674.

[17] ZHANG X M, LIU G R, LAM K Y. Vibration analysis of thin cylindrical shells using wave propagation approach[J]. Journal of Sound and Vibration, 2001, 239(3): 397-403.

Free vibration characteristics of an elliptic cylindrical shell based on the wave propagation method

ZHANG Guanjun, ZHU Xiang, LI Tianyun, MIAO Yuyue

(School of Naval Architecture and Ocean Engineering, Huazhong University of Science and Technology, Wuhan 430074, China)

The free vibration equations of an elliptic cylindrical shell were derived based on the Flügge shell theory. Vibration equations about the circumferential wave number are not decoupled due to the varied circumferential curvature. The shell’s displacements were expanded in double Fourier series in the wave propagation method and the circumferential curvature was expanded in single Fourier series. The partial differential equations with variable coefficients were converted into a set of linear equations which couple with each other about circumferential wave numbers. The natural frequencies of the elliptic cylindrical shell were obtained by solving the coupled equations. The influences of main parameters of the elliptic cylindrical shell, such as ellipticity parameter and shell length ratio, on the vibration characteristics were discussed. The symmetric and anti-symmetric modes of the elliptic cylindrical shell were both considered.

elliptic cylindrical shell; free vibration; natural frequency; ellipticity parameters; anti-symmetric mode

國家自然科學基金(51379083;51479079;51579109);高等學校博士學科點專項科研基金(20120142110051)

2015-10-27 修改稿收到日期: 2016-04-15

張冠軍 男,博士生,1989年生

李天勻 男,博士,教授,1969年生

O327;U663

A

10.13465/j.cnki.jvs.2017.12.031

猜你喜歡
模態振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
This “Singing Highway”plays music
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 57pao国产成视频免费播放| 欧美国产日产一区二区| 国产91丝袜在线播放动漫| 综合色区亚洲熟妇在线| 精品人妻无码中字系列| 欧美精品成人| 国外欧美一区另类中文字幕| 亚洲日韩国产精品无码专区| 99爱在线| 国产又粗又爽视频| 免费国产高清精品一区在线| 国产精品福利在线观看无码卡| 熟妇人妻无乱码中文字幕真矢织江 | 国产精品自拍露脸视频| 人妻21p大胆| 国产91视频观看| 久久亚洲AⅤ无码精品午夜麻豆| 热思思久久免费视频| 狼友视频国产精品首页| 国产丝袜第一页| 欧美啪啪网| 久青草国产高清在线视频| 欧美一区福利| 国产成人精品免费视频大全五级| 欧美乱妇高清无乱码免费| 美女无遮挡拍拍拍免费视频| 国产精品久久久久久久久| 亚洲国产一区在线观看| 色久综合在线| 亚洲一级毛片在线观| 亚洲欧洲免费视频| 国产高清在线精品一区二区三区| 26uuu国产精品视频| 国产午夜一级淫片| 久久香蕉国产线看精品| 色婷婷成人| 一本色道久久88| 国产精品色婷婷在线观看| 无码av免费不卡在线观看| 国产一二视频| 亚洲人成网线在线播放va| 在线不卡免费视频| 黄色片中文字幕| 老司机精品久久| 国产91蝌蚪窝| 国产大片喷水在线在线视频| 亚洲国产精品久久久久秋霞影院| 亚国产欧美在线人成| 国产女人在线| www.99在线观看| 国产精品无码翘臀在线看纯欲| 欧美第九页| 婷婷色婷婷| 婷婷成人综合| 国产精品一区二区国产主播| 色爽网免费视频| 露脸国产精品自产在线播| 成人午夜视频网站| 精品日韩亚洲欧美高清a| 草逼视频国产| 熟女日韩精品2区| 午夜激情婷婷| 干中文字幕| 亚洲综合片| 国产菊爆视频在线观看| 久热中文字幕在线| 国产丝袜第一页| 欧类av怡春院| 日韩一级毛一欧美一国产| 亚洲侵犯无码网址在线观看| 精品伊人久久久大香线蕉欧美| 激情六月丁香婷婷| 九九九精品成人免费视频7| 日本在线国产| 精品欧美一区二区三区久久久| 免费播放毛片| 日本在线亚洲| 国产v欧美v日韩v综合精品| 久久精品中文字幕免费| 欧美日韩午夜视频在线观看| 三上悠亚精品二区在线观看| 在线综合亚洲欧美网站|