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

不同聲子晶體模型的軌道結構振動帶隙對比分析

2022-03-18 00:46:34梁玉雄馮青松陸建飛雷曉燕
振動與沖擊 2022年5期
關鍵詞:振動理論結構

梁玉雄, 馮青松, 陸建飛, 楊 舟, 雷曉燕

(1.華東交通大學 鐵路環境振動與噪聲教育部工程研究中心,南昌 330013;2.江蘇大學 土木工程與力學學院,江蘇 鎮江 212013)

隨著列車運營速度的不斷提高和鐵路運營里程的快速增長,軌道結構引發的振動與噪聲問題日益突出[1-3],因此軌道結構的動力分析一直被國內外學者所關注[4-8]。鐵路中的軌道結構由于其扣件軌枕的周期性布置特點,且鋪設長度通常較長可視為無限長,可將其作為周期性結構對其進行動力特性分析。Wu等[9]利用多層梁模型研究了有限離散支承軌道在高頻范圍的橫向振動,得到的加速度響應與實測數據吻合較好。雷曉燕等[10-14]將鋼軌考慮為歐拉梁,采用傅里葉變換法建立了軌道結構連續彈性單層梁、雙層梁、三層梁車軌耦合模型,進行了軌道結構振動響應研究,并對軌道臨界速度與軌道強振動進行了研究。Grassie等[15]研究了軌道結構的高頻垂向振動。Thompson[16-17]將軌道結構視為周期結構,利用有限元法討論了軌道結構50~5 000 Hz之間的振動特性。

Thompson還計算了周期支承鋼軌的振動衰減率,發現在“pinned-pinned”頻率附近會產生衰減區域。馬龍祥等[18]將鋼軌視為無限長的歐拉梁,利用周期結構頻域內響應的性質和疊加原理,對移動諧振荷載下軌道結構的振動進行了研究。Zhang等[19]采用2.5D有限元方法建立了離散支承軌道模型,研究了軌道結構的振動響應及軌枕墊彈簧剛度對其的影響,并與現場實測結果進行了對比。Sheng等[20-21]通過基于波數域的周期支承結構在諧波荷載作用下動態響應分析,得到了周期軌道結構振動波傳播常數和諧振特性。

近年來隨著聲子晶體理論的興起[22-26],許多學者基于聲子晶體理論研究了彈性波在周期性軌道結構中傳播特性,聲子晶體是指由兩種或兩種以上介質組成的具有彈性波帶隙特性的周期性復合材料或結構,當彈性波在聲子晶體中傳播時,受內部周期結構的作用,某些頻率范圍內的彈性波不能傳播,相應的頻率范圍稱為帶隙(也稱禁帶);而其他頻率范圍內的彈性波可以傳播,相應的頻率范圍稱為通帶。Wang等[27-28]利用聲子晶體理論通過歐拉梁模型分析了不考慮阻尼時有砟軌道結構帶隙行為和形成機制,并通過鐵木辛柯梁模型和現場測試研究了不考慮阻尼時有序和隨機無序高鐵無砟軌道中的波傳播問題。孟鐸[29]利用聲子晶體理論分析了周期性軌道結構振動帶隙特性,并對動力吸振器進行了研究。易強等[30]建立周期性軌道結構模型并將鋼軌考慮為鐵木辛柯梁,結合聲子晶體理論研究了彈性波在周期性軌道結構中傳播特性,得出在帶隙范圍內彈性波在軌道結構中無法自由傳播,且外界激勵也無法向系統輸入能量,而通帶范圍內可進行能量的輸入與傳播。馮青松等[31]通過平面波展開法將鋼軌等效為鐵木辛柯梁,對周期軌道結構垂向振動帶隙特性進行了研究,與有限元計算得到的傳輸特性進行了對比分析,得到了周期性結構中彎曲波的傳遞特性,并分析了軌道結構主要參數對帶隙特性的影響。

上述文獻對軌道結構的振動響應研究較多,對軌道結構中振動波的特性研究較少。而在采用聲子晶體理論對軌道結構振動波特性的研究中,尚未對歐拉梁理論模型及鐵木辛柯梁理論模型的主要區別和在不同軌道結構類型中的適用性進行深入研究。鑒于此,本文利用軌道結構的周期性支承特點引入布洛赫周期性邊界條件,按是否考慮軌枕及道砟影響分別建立單層歐拉梁、單層鐵木辛柯梁,雙層歐拉梁、雙層鐵木辛柯梁四種聲子晶體分析模型,計算了不考慮阻尼和考慮阻尼兩種情況下軌道結構自由振動的特征波能帶結構,對不同聲子晶體理論分析模型的適用性及針對性,以及軌枕和道砟、阻尼對帶隙的影響進行了分析,通過有砟軌道和無砟軌道兩種典型軌道結構進行現場力錘激振測試,驗證各理論模型分析結果。

1 不同聲子晶體分析模型及其傳遞矩陣推導

1.1 不同鋼軌振動分析模型簡介

本文將軌道結構按是否考慮軌枕及道砟影響,將軌下支承分別簡化為單層支承梁模型和雙層支承梁模型,支承形式為彈性點支承,再分別按歐拉梁理論和鐵木辛柯梁理論進行分析。

單層支承梁理論分析模型中軌道結構簡化為由鋼軌、扣件組成的無限周期結構,如圖1所示。扣件簡化為彈簧單元,忽略軌枕及道砟的影響。

圖1 單層支承梁理論分析模型

雙層支承梁理論分析模型中軌道結構簡化為由鋼軌、扣件、軌枕和道砟組成的無限周期結構,如圖2所示。扣件與道砟簡化為彈簧單元,軌枕簡化為質量塊單元。軌道結構參數如表1所示。

圖2 雙層支承梁理論分析模型

表1 軌道結構參數

考慮扣件及道砟的阻尼時,扣件剛度和道砟剛度采用復剛度形式[32]

扣件豎向剛度kcrf=krf(1+iηp)

(1)

扣件垂向剛度kcrm=krm(1+iηp)

(2)

道砟豎向剛度kcsf=ksf(1+iηs)

(3)

道砟縱向剛度kcsm=ksm(1+iηs)

(4)

式中,ηp=0.20,ηs=0.068分別為扣件墊板阻尼損耗因子和道砟阻尼損耗因子。

1.2 傳遞矩陣推導過程

1.2.1 等直歐拉梁彎曲振動傳遞矩陣

根據達朗貝爾原理,等直歐拉梁的彎曲振動微分方程為

(5)

式中:υ(x,t)是梁的豎向位移;E為梁的模量;A和I分別為梁的截面面積和慣性矩;ρ為梁的線密度。

解得微分方程的解為

Y=c1chβx+c2shβx+c3cosβx+c4sinβx

(6)

其中c1=(L1+L3)/2,c2=(L2+L4)/2,c3=(L1-L3)/2,c4=(L2-L4)/2

令A=(chβx+cosβx)/2,B=(shβx+sinβx)/2,C=(chβx-cosβx)/2,D=(shβx-sinβx)/2。

由材料力學中撓度、轉角、彎矩和剪力的關系

(7)

將式(6)代入式(7)得出θ、M、Q的一般表達式,當x=0時,令Y=Y0,M=M0,Q=Q0,求得Y、θ、M、Q的一般表達式,將第0跨元胞左右側狀態向量其寫成矩陣形式為

(8)

令Tb(x)=

(9)

式(9)為歐拉梁彎曲振動的傳遞矩陣。

1.2.2 等直鐵木辛柯梁彎曲振動傳遞矩陣

根據梁的Timoshenko理論,梁的彎曲振動微分方程如下

(10)

(11)

式中:υ(x,t)是梁的豎向位移;θ(x,t)是梁截面的彎曲轉角;G為梁的剪切模量;A和I分別為梁的截面面積和慣性矩,κ為鐵木辛柯梁的剪切系數。

令υ(x,t)=Y(x)eiωt,θ(x,t)=Θ(x)eiωt,ω為圓頻率,由式(10)和式(11)得到

(12)

解式(12)微分方程,求得方程的解具有以下兩種形式

Y(x)=C1coshk1x+C2sinhk1x+C3cosk2x+

(13)

Y(x)=C1cosk1x+C2sink1x+C3cosk2x+

(14)

其中C1,C2,C3,C4為待定系數。

(15)

(16)

將Timoshenko梁彎曲振動時的狀態向量寫為

(17)

式(17)中,

U(ω)=

(18)

進而得到

(19)

Tb(x)=U(ω)A(x)A-1(0)U(ω)-1

(20)

式中,Tb(x)為Timoshenko梁垂向彎曲振動的傳遞矩陣。

1.2.3 等直梁單層支承點處傳遞矩陣

圖3為單層支承梁模型第0跨元胞鋼軌扣件支承點處受力圖。

圖3 單層支承梁扣件節點處受力圖

根據節點處位移連續條件及力平衡方程可以得到,扣件節點右側狀態向量Ψr(0+)和左側狀態向量Ψr(0-)具有以下關系,

Ψr(0+)=TΨr(0-)

(21)

式中,

qr(x,ω)={ν(x,ω),θ(x,ω),μ(x,ω)}T

fr(x,ω)={Q(x,ω),M(x,ω),N(x,ω)}T

(22)

式(22)為鋼軌在單層支承扣件節點處傳遞矩陣,其中krf,krr,krn分別為扣件的豎向剛度,彎曲剛度、縱向剛度。

1.2.4 等直梁雙層支承點處傳遞矩陣

圖4為雙層支承梁模型第0跨元胞支承點處受力圖,第0跨元胞中包含了扣件間距內的鋼軌、扣件彈簧、軌枕質量塊、道砟彈簧。

圖4 雙層支承梁扣件節點處受力圖

根據胡克定律,扣件力的表達式為

Qt(0)=-krn[ur(0)-us(0)]

Mt(0)=-krr[θr(0)-θs(0)]

Nt(0)=-krf[vr(0)-vs(0)]

(23)

式中:ur,us分別為鋼軌和軌枕的豎向位移;θr,θs分別為鋼軌和軌枕的彎曲轉角;vr,vs分別為鋼軌和軌枕的豎向位移。

寫成矩陣的形式為

fr(0+)=fr(0-)+E1qr(0)+E2qs(0)

(24)

其中

fr(0+)={Qr(0),Mr(0),Nr(0)}T

qr(0)={νr(0),θr(0),μr(0)}T

qs(0)={νs(0),θs(0),μs(0)}T

由扣件節點力的平衡方程和鋼軌的動力學方程可得

(25)

式中:Kr={krf,krr,krn}T;Ks={ksf,ksr,ksn}T;ksf,ksr,ksn分別為道砟的豎向剛度,彎曲剛度、縱向剛度。

(26)

并將式(26)寫成矩陣形式可得,

qs=E3qr

(27)

其中

E3=

(28)

將式(27)代入式(24)可得

fr(0+)=fr(0-)+(E1+E2E3)qr(0)

(29)

根據式(29)及節點處的位移連續條件,扣件節點右側狀態向量Ψr(0+)和左側狀態向量Ψr(0-)具有以下關系

Ψr(0+)=TΨr(0-)

(30)

(31)

式(31)為鋼軌在雙層支承扣件節點處傳遞矩陣,其中

(32)

(33)

(34)

1.2.5 等直梁縱向自由振動傳遞矩陣

等直梁的縱向振動微分方程為

(35)

式中:u(x,t)是梁的豎向位移;E為梁的模量;A為梁的截面面積;ρ為梁的線密度。

令μ(x,t)=U(x)sin(ωt+φ)代入式(35)得

(36)

式(36)解微分方程的解為

U=c1cosαx+c2sinαx

(37)

由材料力學

(38)

令x=0時,U=U0,N=N0,作為已知條件,求得式(37)中的系數

(39)

將式(39)中系數代入式(37),將第0跨元胞左右側狀態向量其寫成矩陣形式為

(40)

(41)

式(41)為等直梁縱向振動的傳遞矩陣。

1.3 聲子晶體分析模型振動特性求解方程

對于無縫線路鋼軌這一聲子晶體結構,由布洛赫(Bloch)定理,以及左右端點狀態向量間的傳遞矩陣,可得:

|TlUl-e-iκlI|=0

(42)

式中:I為單位矩陣;Tl為扣件節點處傳遞矩陣;Ul為振動傳遞矩陣;κ為對于給定頻率ω由式(42)求得的鋼軌結構中傳播的特征波的復波數,其實部表示特征波相位的改變,虛部表示特征波的衰減;l為周期長度,即相鄰軌枕扣件間的距離。

2 不同聲子晶體分析模型分析結果

2.1 無阻尼狀況軌道垂向彎曲振動能帶結構分析

圖5為無阻尼聲子晶體分析模型計算得到的軌道結構垂向彎曲振動能帶結構圖。

(a) 第一種特征波能帶結構(實部)

從計算結果可知,在單層梁模型中,歐拉梁理論與鐵木辛柯梁理論在低頻范圍(0~250 Hz),能帶結構圖的實部和虛部二者無明顯區別,隨著頻率的提高,鐵木辛柯理論分析得到的第二種特征波帶隙的第2階帶隙起止位置(1 013~1 028 Hz)較歐拉梁理論分析結果(1 342~1 352 Hz)有所不同,增加了一條新的高頻帶隙(2 718~2 723 Hz)。

在雙層支承梁模型中,采用歐拉梁理論與鐵木辛柯梁理論在低頻范圍(0~250 Hz)有明顯不同,且隨著頻率的提高二者計算結果差異越來越大。鐵木辛柯梁理論分析得到的第一種特征波出現了0~135 Hz的完全帶隙(波矢實部為零、波矢虛部非零)、136~190 Hz的不完全帶隙(波矢實部非零、波矢虛部非零)、191~3 500 Hz的完全帶隙,而在歐拉梁理論分析得到結果中0~134 Hz為不完全帶隙,134~3 500 Hz則為完全帶隙;鐵木辛柯理論分析得到的第二種特征波帶隙的第2階帶隙位置(1 013~1 028 Hz)較歐拉梁理論分析結果(1 342~1 352 Hz)有所不同,增加了一條高頻帶隙(2 718~2 723 Hz)。

這表明單層梁模型在低頻(0~250 Hz)范圍,基于聲子晶體理論進行無縫線路鋼軌減振控制時,采用歐拉梁理論與鐵木辛柯梁理論并無顯著差別,但在中高頻(250 Hz以上)范圍減振控制時,二者有顯著不同。而在雙層梁模型中,在0~3 500 Hz范圍基于聲子晶體理論進行無縫線路鋼軌減振控制時,采用歐拉梁理論與鐵木辛柯梁理論在0~250 Hz內和1 000 Hz以上均有顯著差別。

2.2 有阻尼狀況軌道垂向彎曲振動能帶結構分析

圖6為考慮阻尼時單層鐵木辛柯梁理論模型和雙層鐵木辛柯梁理論模型,計算得到的鋼軌能帶結構圖。

單層歐拉梁模型和雙層歐拉梁模型,計算得到的軌道能帶結構圖對比情況與鐵木辛柯梁類同。由圖6可知,單層鐵木辛柯梁模型和雙層鐵木辛柯梁模型的兩種特征波均在190~500 Hz范圍由通帶變為不完全帶隙(波矢實部非零、波矢虛部非零),低頻帶隙截止頻率至500 Hz受阻尼影響較大。此外特征波波矢虛部在低頻帶隙截止頻率處反映出來的衰減峰值受阻尼影響較大,在雙層鐵木辛柯梁模型中尤為顯著,阻尼主要影響129~500 Hz內的帶隙。

考慮鋼軌扣件阻尼和道砟阻尼后,特征波的帶隙由無阻尼狀況下的完全帶隙改變為不完全帶隙,無阻尼狀況時的計算值為零的虛波矢和實波矢,在有阻尼狀況時則出現非零值,即通帶在阻尼的影響下也會有微小的衰減,禁帶的頻帶寬度會有微小的展寬,禁帶的中心位置受阻尼的影響可忽略不計。

2.3 軌道縱向振動能帶結構分析

圖7為有阻尼和無阻尼時,單層梁理論模型和雙層梁理論模型計算得到的鋼軌縱向振動能帶結構。無阻尼模型中均出現一條不完全帶隙,單層梁模型為0~70 Hz,雙層模型為0~77 Hz。從圖7(b)中能帶結構圖虛部的衰減幅值可以看出,有阻尼單層梁模型和雙層梁模型的主要衰減區域仍然分別位于0~70 Hz內和0~77 Hz內,但單層梁理論模型特征波的帶隙由0~70 Hz內的完全帶隙擴展為0~479 Hz內的不完全帶隙,雙層理論模型特征波的帶隙擴展為0~15 Hz的完全帶隙和16~480 Hz的不完全帶隙。

(a) 縱向振動特征波能帶結構(實部)

3 現場測試驗證結果

為驗證本文中不同理論分析模型的適用性,選取了昌贛高鐵吉安段和京九鐵路蓮塘段典型的無砟軌道和有砟軌道段落,采用LMS數據采集系統分別進行力錘激振測試驗證,如圖8所示。測試時選擇20跨鋼軌,Pi為第i跨跨中測點位置,在Pi處布置垂向加速度傳感器。通過力錘在不同Pi處對軌道激振時測量各測點的響應,并以錘擊點單位錘擊力的跨中響應作為分析軌道結構彈性波傳播特性的參考值,通過分析振動傳遞系數可以得到軌道結構中波的傳播特性,與理論計算得到的帶隙進行對比。

因測試設備和傳感器測試量程受限,測試僅得到了0~1 500 Hz范圍的有效測試數據,通過LMS Test.lab多通道分析儀,將力錘輸入力的頻譜和同步的振動響應采用FRF(頻率響應函數)方法和PSD(功率譜密度譜方法)進行處理,對頻響函數的振幅、傳輸特性與能帶結構理論計算結果進行了對比分析。

圖9(a)為實測得到的有砟軌道垂向振動頻響實部曲線與理論計算帶隙對比分析圖,實測結果中0~135 Hz內頻響曲線有先迅速衰減后又迅速增強,并非全頻段完全衰減,有阻尼雙層鐵木辛柯梁理論分析得到的0~135 Hz范圍內為不完全帶隙,振動波實部不為零,虛部非零衰減,這與振動波傳播時在帶隙1起始頻率處迅速衰減,在帶隙1終止頻率處迅速增強的規律基本一致,因此測試結果與考慮阻尼的0~135 Hz范圍的理論帶隙1基本吻合。實測頻響曲線在帶隙2起始頻率(136 Hz)處衰減,衰減持續至帶隙2終止頻率(190 Hz)時迅速增加的規律,與有阻尼雙層鐵木辛柯梁分析得到的理論帶隙2基本吻合,證明了扣件阻尼和道砟阻尼使得無阻尼模型的理論帶隙由完全帶隙改變為不完全帶隙,帶隙的頻帶寬度發生了展寬現象。

(a) 實測FRF頻響實部曲線與理論帶隙

圖9(b)為實測得到的有砟軌道垂向振動傳輸系數曲線與理論計算得到的波矢虛部曲線及帶隙對比分析圖,有砟軌道實測垂向振動傳輸系數曲線與理論計算波矢虛部曲線及帶隙對比分析圖中,傳輸系數曲線顯示在0~250 Hz內彈性波傳播被抑制,只有非常小的振動能量傳遞,傳輸系數曲線因實際軌道結構中阻尼和失諧等因素影響表現為非平滑曲線,從而證實了在低頻段(0~250 Hz)內有阻尼雙層鐵木辛柯梁理論帶隙1和理論帶隙2位置基本吻合。而在其他頻段頻響虛部曲線出現了若干理論分析中未出現的窄帶隙。

在圖10(a)無砟軌道垂向振動頻響實部曲線與理論計算帶隙對比分析圖中,以及圖10(b)無砟軌道實測傳輸曲線與理論計算波矢曲線對比圖中,振動頻響實部曲線和傳輸系數曲線顯示0~135 Hz頻段(有阻尼單層支承鐵木辛柯梁理論帶隙1)內彈性波傳播明顯被抑制,在800~1 050 Hz頻段內彈性波傳播被抑制, 有阻尼單層支承鐵木辛柯梁理論帶隙2(1 013~1 028 Hz)包含在該頻段(800~1 050 Hz)內,實測傳輸曲線在理論帶隙之外的其他頻段出現了若干衰減段,即在低頻段(0~250 Hz)內的帶隙情況與單層支承鐵木辛柯梁理論帶隙1基本吻合,其他頻段的帶隙位置和寬度與理論值理論禁帶展寬后的帶隙并不完全吻合。

(a) 無砟軌道實測FRF頻響實部曲線曲線與理論帶隙

根據以上對比分析,低頻段內,無砟軌道相比有砟軌道,其實測結果與理論帶隙對應關系更為吻合,有砟軌道FRF函數中幅值較小的頻帶較寬,與理論計算結果有一定差異,這種差異的原因與實際軌道結構中材料、幾何缺陷以及制造誤差引起的周期失諧,軌枕、扣件、道砟等軌道結構參數與理論分析參數的差異,現場測試條件,測試儀器的靈敏度等多種因素有關,結合近年來學者們對失諧周期結構振動特性的研究文獻,周期結構的失諧會引起結構中傳播波或振動在節點處反射,導致波或振動能量局限于一個很小的幾何范圍內,形成局部振蕩,出現局部化現象, 使得通帶將變窄, 帶隙將變寬[33],此外隨著失諧程度的增加,結構的波動局部化程度增強,甚至使得禁帶頻率個數增加[34],因此考慮到實際軌道結構中周期失諧的復雜性,周期失諧為引起的實測結果與理帶隙的差異的主要原因之一。

4 結 論

本文通過分析四種聲子晶體分析模型鐘軌道結構自由振動的特征波能帶結構,及有砟軌道和無砟軌道兩種典型軌道結構的現場力錘激振測試,針對不同模型的主要區別及其在不同類型軌道結構中的適用性,得到以下主要結論:

(1) 無阻尼單層歐拉梁理論模型與無阻尼鐵木辛柯梁理論模型在低頻范圍(0~250 Hz)二者無顯著區別,低頻帶隙的形成與鋼軌的固有共振頻率相關,隨著頻率的提高,鐵木辛柯理論模型分析得到的第二種特征波帶隙的第2階帶隙起止位置(1 013~1 028 Hz)與帶寬均較歐拉梁理論分析結果(1 342~1 352 Hz)有所不同,且增加了一條新的高頻帶隙(2 718~2 723 Hz)。

(2) 有阻尼狀況下的軌道結構聲子晶體模型彎曲特征波和縱向振動特征波的帶隙,與無阻尼狀況下相比,鋼軌的禁帶的頻帶寬度受阻尼影響會有微小的展寬,禁帶的中心位置受阻尼的影響可忽略不計。

(3) 有砟軌道鋼軌垂向振動傳輸特性實測結果與考慮阻尼時的雙層鐵木辛柯梁模型分析結果在低頻段(0~250 Hz)基本吻合。而無砟軌道的鋼軌振動傳輸特性實測結果與考慮阻尼時的單層鐵木辛柯梁模型分析結果在低頻段(0~250 Hz)均基本吻合。

(4) 采用聲子晶體帶隙理論減振控制時,因扣件和道砟共同影響鋼軌在低頻范圍存在剪切變形相關,對于有砟軌道應采用雙層鐵木辛柯梁模型。對于當無砟軌道在低頻(0~250 Hz)范圍可采用單層歐拉梁模型或單層鐵木辛柯梁模型,在高頻范圍(250 Hz以上)應采用鐵木辛柯梁理論。

猜你喜歡
振動理論結構
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
堅持理論創新
當代陜西(2022年5期)2022-04-19 12:10:18
神秘的混沌理論
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
理論創新 引領百年
相關于撓理論的Baer模
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
主站蜘蛛池模板: 永久免费精品视频| 国产精品一区不卡| 五月丁香伊人啪啪手机免费观看| 亚洲第一视频网| 97在线观看视频免费| 国产精品短篇二区| 欧美亚洲欧美区| 亚洲视频欧美不卡| 99re热精品视频国产免费| 香蕉eeww99国产精选播放| 国产午夜不卡| 欧美日韩va| 99视频免费观看| www.91在线播放| 欧美a级在线| 99国产精品一区二区| 一级毛片基地| 福利国产微拍广场一区视频在线| 亚洲伦理一区二区| 91成人免费观看| 国产精品区视频中文字幕| 国产精品污视频| 欧美日本激情| 亚洲欧美另类中文字幕| 欧洲成人免费视频| 美女被操91视频| 欧美色亚洲| 欧美国产在线看| 精品国产自在现线看久久| 特级精品毛片免费观看| 国产午夜人做人免费视频中文| 国产v欧美v日韩v综合精品| 最新无码专区超级碰碰碰| 青青操视频在线| 国产尹人香蕉综合在线电影| 亚洲,国产,日韩,综合一区| 996免费视频国产在线播放| 国产精品xxx| 欧美中文字幕一区| 亚洲AV永久无码精品古装片| 456亚洲人成高清在线| 精品夜恋影院亚洲欧洲| 伊人AV天堂| 亚洲av无码牛牛影视在线二区| 国产成人免费视频精品一区二区| 国产精品视频3p| 久青草网站| 国产欧美精品一区二区| 国产欧美在线观看精品一区污| 无码一区二区三区视频在线播放| 日韩在线网址| 亚洲Av激情网五月天| Aⅴ无码专区在线观看| 久久无码免费束人妻| 成色7777精品在线| 国产麻豆精品手机在线观看| lhav亚洲精品| 欧美视频二区| 久久精品国产电影| 就去吻亚洲精品国产欧美| 国产成人久久777777| 97在线碰| 中文无码影院| 国内自拍久第一页| 97超碰精品成人国产| 国产尹人香蕉综合在线电影| 岛国精品一区免费视频在线观看 | 狠狠v日韩v欧美v| 四虎综合网| 精品国产中文一级毛片在线看| 国内精品免费| 国产黑丝一区| 在线视频亚洲欧美| 国产呦精品一区二区三区网站| 欧美成人午夜视频| 久久精品人妻中文系列| 毛片免费视频| 波多野结衣无码中文字幕在线观看一区二区| 激情爆乳一区二区| 亚洲AV无码乱码在线观看裸奔 | 国产欧美在线观看视频| 精品国产香蕉在线播出|