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

船體橫蕩運動對轉(zhuǎn)子-浮筏氣囊耦合系統(tǒng)非線性動力學(xué)特性的影響*

2022-12-28 05:09:24王軍偉杜曉蕾
潤滑與密封 2022年12期
關(guān)鍵詞:系統(tǒng)

王軍偉 李 明 謝 旋 杜曉蕾

(西安科技大學(xué)理學(xué)院 陜西西安 710054)

近年來,隨著我國船舶制造技術(shù)水平的不斷提高,船舶的制造向著大型化、復(fù)雜化發(fā)展。然而艦船在海上作業(yè)時通常會受到風(fēng)浪的作用使得船體發(fā)生周期性的搖蕩等牽連運動,這些運動會對船舶結(jié)構(gòu)產(chǎn)生影響,尤其是船用轉(zhuǎn)子-軸承系統(tǒng),從而影響艦船的可靠性和安全性。為保證艦船能平穩(wěn)航行,研究這些運動對船用轉(zhuǎn)子系統(tǒng)的影響具有重要意義。

傳統(tǒng)對于轉(zhuǎn)子動力學(xué)的分析主要是針對基礎(chǔ)固連于地面上的轉(zhuǎn)子-軸承系統(tǒng),比如泵、大型風(fēng)機和發(fā)電機組的轉(zhuǎn)子系統(tǒng),在分析其動力學(xué)特性的時候都是假設(shè)轉(zhuǎn)子的基礎(chǔ)靜止不動而且剛度較大。目前對于考慮牽連運動的轉(zhuǎn)子動力學(xué)的研究主要集中在機載、船載和車載方面。楊蛟和曹樹謙[1]利用航空發(fā)動機雙轉(zhuǎn)子模型試驗臺分別模擬飛機做橫滾、俯仰和偏航運動時系統(tǒng)的振動情況。林富生等[2-4]考慮轉(zhuǎn)子系統(tǒng)的多種因素,如裂紋、等加速、等減速,針對機動飛行條件下的轉(zhuǎn)子系統(tǒng)的非線性動力學(xué)行為進行了討論。HOU等[5]采用拉格朗日方程推導(dǎo)了飛機旋轉(zhuǎn)時含內(nèi)間隙和赫茲接觸力的轉(zhuǎn)子-滾珠軸承系統(tǒng)的非線性動力學(xué)方程,對機動載荷作用下的系統(tǒng)響應(yīng)進行了數(shù)值分析。祝長生和陳擁軍[6-7]利用Lagrange方程建立了飛行器在任意空間作機動飛行時機載轉(zhuǎn)子系統(tǒng)的動力學(xué)模型,并討論了機動飛行對發(fā)動機轉(zhuǎn)子系統(tǒng)動力學(xué)特性的影響。李杰等人[8]考慮航空發(fā)動機雙轉(zhuǎn)子中介軸承的耦合作用及陀螺力矩的影響,建立了機動飛行條件下雙轉(zhuǎn)子-滾動軸承支承耦合系統(tǒng)動力學(xué)模型,研究了不同轉(zhuǎn)速比下轉(zhuǎn)子的動力學(xué)特性。HAN和LI[9]基于非慣性參考系,建立了垂蕩運動下船用轉(zhuǎn)子-軸承系統(tǒng)的動力學(xué)模型,利用數(shù)值方法分析了系統(tǒng)的穩(wěn)態(tài)響應(yīng)。

氣囊隔振器也被稱作為空氣彈簧,因其具有固有頻率低、承載力大、剛度可調(diào)、無蠕變等特點,成為大型旋轉(zhuǎn)機械中性能優(yōu)異的隔振器。研究人員對氣囊浮筏隔振系統(tǒng)進行了深入研究。徐偉等人[10]將氣囊隔振裝置應(yīng)用到船舶旋轉(zhuǎn)機械的隔振之中,發(fā)現(xiàn)氣囊隔振器能明顯地減小動力設(shè)備的振動能量。呂志強等[11]將氣囊隔振器嵌入到浮筏裝置中,通過調(diào)節(jié)氣囊中的壓力保持了浮筏的平衡性。施亮等人[12]建立了主機氣囊隔振器對中姿態(tài)響應(yīng)的線性化模型,為系統(tǒng)優(yōu)化設(shè)計提供了理論依據(jù)。ZHANG等[13]介紹了一種用于船舶推進軸系的智能浮筏系統(tǒng)(IFRS)并對IFRS的力學(xué)特性進行了建模和分析。趙興乾等[14]以某型船用推進電機隔振系統(tǒng)校核為背景,建立隔振系統(tǒng)簡化模型,分析了船舶不同姿態(tài)下軸系對中校核。ZHAO等[15-16]先將氣囊-浮筏隔振結(jié)構(gòu)嵌入到船用旋轉(zhuǎn)機械的轉(zhuǎn)子-軸承系統(tǒng)中,建立了多維耦合的氣囊-浮筏隔振力學(xué)模型,探討了系統(tǒng)的非線性動力學(xué)特性,接著建立了在基礎(chǔ)激勵作用下考慮轉(zhuǎn)子-浮筏氣囊耦合系統(tǒng)的動力學(xué)模型,分析了系統(tǒng)的動力學(xué)特性。李鵬超[17]在文獻[15-16]的基礎(chǔ)上研究了沖擊激勵下系統(tǒng)的動力學(xué)特性,并在氣囊隔振器中安裝限位器來對船用旋轉(zhuǎn)機械系統(tǒng)的沖擊響應(yīng)進行限制,最后對比分析了參數(shù)變化對系統(tǒng)非線性動力學(xué)行為的影響。

轉(zhuǎn)子-軸承系統(tǒng)作為旋轉(zhuǎn)機械的核心部件,在運轉(zhuǎn)的時候難免會產(chǎn)生振動,這些振動會對艦船產(chǎn)生諸多不利影響,比如振動產(chǎn)生的噪聲會使得船員的工作和居住環(huán)境惡化,振動會使船體結(jié)構(gòu)發(fā)生疲勞損壞,還會影響艦船的隱蔽性等等。采用浮筏-氣囊隔振裝置可以減少振動對艦船的影響,但轉(zhuǎn)子-軸承系統(tǒng)的運動情況就會變得更為復(fù)雜。本文作者著重研究了船體橫蕩運動對轉(zhuǎn)子-浮筏氣囊耦合系統(tǒng)的影響。

1 橫蕩作用下系統(tǒng)動力學(xué)模型

1.1 系統(tǒng)運動方程的建立

圖1所示為艦船搖蕩運動和船用轉(zhuǎn)子-浮筏氣囊系統(tǒng)的結(jié)構(gòu)示意圖,其中O-X0Y0Z0為艦船相對于地面的坐標(biāo)系,o1-x1y1z1為轉(zhuǎn)子相對于艦船的坐標(biāo)系,o2-x2y2z2為浮筏相對于艦船的坐標(biāo)系。為了簡化問題的分析,現(xiàn)作以下假設(shè):將轉(zhuǎn)子看作具有集中質(zhì)量的圓盤[18],質(zhì)量為m1;浮筏和軸承支座之間剛性連接,且軸承為短軸承,浮筏和軸承可以看作為一個整體,質(zhì)量為m2;氣囊可看作為在x、y方向上同時具有線性剛度和阻尼的彈簧,其中,剛度分別為kx、ky;阻尼分別為dx、dy;假定船體橫蕩的時候,轉(zhuǎn)子系統(tǒng)在z方向的運動與船體一致,即忽略系統(tǒng)z方向的運動。

圖1 船用氣囊浮筏轉(zhuǎn)子-軸承系統(tǒng)示意

基于短軸承(軸承長度遠(yuǎn)遠(yuǎn)小于其直徑)理論,考慮系統(tǒng)在船體橫蕩作用下(即在y方向的運動),依據(jù)牛頓運動定律,系統(tǒng)的運動微分方程可表示為

(1)

式中:F為非線性油膜力,N;e為偏心距,m;ω為轉(zhuǎn)子轉(zhuǎn)速,rad/s。

1.2 非線性油膜力模型

轉(zhuǎn)子與筏體一般通過滑動軸承連接,文中基于短軸承油膜力(軸承長度遠(yuǎn)遠(yuǎn)小于其直徑,即油膜壓力沿軸向的變化遠(yuǎn)遠(yuǎn)小于沿周向的變化)簡化Reynolds方程得出油膜力的表達式如下:

(2)

式中:p為油膜壓力;φ、θ分別為偏位角和周向方位角;油膜厚度h=c+ecosθ=c(1+εcosθ),ε為偏心率,ε=e/c,c為油膜間隙;e和η分別為軸承偏心量和潤滑油黏度。

圖2 滑動軸承結(jié)構(gòu)示意

利用半Sommerfeld條件對式(2)進行2次積分,同時認(rèn)為軸承兩端油膜壓力為0,即有邊界條件p|z=-L/2=p|z=L/2=0,得到非線性油膜力的徑向和周向分力Fr、Fτ表達式如下:

(3)

通過坐標(biāo)變換,將油膜力在徑向和周向分力Fr,F(xiàn)τ變換到x、y方向

(4)

1.3 運動方程量綱一化

利用油膜間隙c和轉(zhuǎn)子質(zhì)量m1對式(1)進行量綱一化,目的是消除量綱的影響,使得研究的問題更加廣泛。相關(guān)的參數(shù)表達式如表1所示。

表1 量綱一化參數(shù)表達式

(5)

式中:質(zhì)量比n=m2/m1;頻率比ν=ω/ω0。

方程(5)為多個量綱一化參數(shù)控制的二階非線性微分方程組。對于這類方程的求解,一般先對微分方程進行降階處理,得到8個一階方程組:

(6)

其中量綱一化后的非線性油膜力在徑向和周向的表達式為

(7)

油膜力轉(zhuǎn)換到x、y方向

式中:Xr=X1-X2,Yr=Y1-Y2。

2 非線性動力學(xué)特性分析

采用Runge-Kutta法對式(6)進行數(shù)值求解,方程中的相關(guān)參數(shù)取值范圍如表2所示。

表2 相關(guān)參數(shù)取值范圍

2.1 不考慮橫蕩運動時系統(tǒng)的動力學(xué)特性

為了對比研究橫蕩運動對系統(tǒng)動力學(xué)特性的影響,文中首先計算了在相同的系統(tǒng)參數(shù)下無橫蕩運動的系統(tǒng)動力學(xué)特性,結(jié)果如圖3所示。

圖3所示為量綱一轉(zhuǎn)子轉(zhuǎn)速Ω=0.4~3.6時系統(tǒng)在無橫蕩作用下穩(wěn)態(tài)響應(yīng)分岔圖及最大Lyapunov指數(shù),此時的系統(tǒng)參數(shù)為:σ=3,α=0.05,n=60,Ωxn=Ωyn=0.6,λ=0.2,ζ=0.1。在該參數(shù)下系統(tǒng)主要受到不平衡力和非線性油膜力的作用。在轉(zhuǎn)子轉(zhuǎn)速較低時,即Ω=0.4~2.23,系統(tǒng)受到不平衡力的影響較大,系統(tǒng)為穩(wěn)定的周期1運動狀態(tài);隨著轉(zhuǎn)速的繼續(xù)增大,系統(tǒng)的運動狀態(tài)出現(xiàn)分岔現(xiàn)象,由原來的周期1運動狀態(tài)經(jīng)倍周期分岔為周期2運動狀態(tài);當(dāng)轉(zhuǎn)速大于2.62時,系統(tǒng)的運動狀態(tài)又由周期2變?yōu)橹芷?運動狀態(tài);當(dāng)轉(zhuǎn)子轉(zhuǎn)速繼續(xù)增大時,系統(tǒng)發(fā)生準(zhǔn)周期分岔現(xiàn)象;轉(zhuǎn)子轉(zhuǎn)速Ω=3.32時,Lyapunov指數(shù)大于0,系統(tǒng)的運動表現(xiàn)為混沌狀態(tài)。綜合來看系統(tǒng)無橫蕩作用時的動力學(xué)行為表現(xiàn)為:周期1→周期2→周期1→擬周期→混沌。

圖3 無橫蕩作用下系統(tǒng)隨轉(zhuǎn)速變化的

2.2 橫蕩作用下轉(zhuǎn)子轉(zhuǎn)速對系統(tǒng)動力學(xué)特性的影響

圖4所示為有橫蕩作用下系統(tǒng)隨轉(zhuǎn)子轉(zhuǎn)速變化的穩(wěn)態(tài)響應(yīng)分岔圖及最大Lyapunov指數(shù),系統(tǒng)參數(shù)為:σ=3,α=0.05,n=60,Ωxn=Ωyn=0.6,λ=0.2,ζ=0.1,A=300。在該參數(shù)下,系統(tǒng)受到橫蕩慣性力、不平衡力和非線性油膜力的共同作用。當(dāng)轉(zhuǎn)速Ω=0.4~2.23時與無橫蕩作用相比,系統(tǒng)動力學(xué)特性由原來的周期1運動變?yōu)閿M周期運動;當(dāng)轉(zhuǎn)子轉(zhuǎn)速繼續(xù)增大,在Ω=2.29~2.6區(qū)間內(nèi),系統(tǒng)的運動狀態(tài)由準(zhǔn)周期分岔為兩支,但仍是擬周期運動狀態(tài);當(dāng)轉(zhuǎn)子轉(zhuǎn)速大于2.6,經(jīng)準(zhǔn)周期分岔的兩支又合為一支,系統(tǒng)處于擬周期運動狀態(tài);當(dāng)轉(zhuǎn)子轉(zhuǎn)速大于3.32時,轉(zhuǎn)子做高速旋轉(zhuǎn),此時系統(tǒng)在橫蕩慣性力、不平衡力和非線性油膜力的共同作用下,最大Lyapunov指數(shù)大于0,系統(tǒng)的運動表現(xiàn)為混沌狀態(tài);與無橫蕩作用時相比,系統(tǒng)進入混沌的轉(zhuǎn)速相差不多,說明橫蕩運動并沒有使得轉(zhuǎn)子提前或滯后進入混沌運動狀態(tài)。

圖4 有橫蕩作用時系統(tǒng)隨轉(zhuǎn)速變化的穩(wěn)

圖5所示為量綱一轉(zhuǎn)子轉(zhuǎn)速Ω=1.07時,系統(tǒng)的穩(wěn)態(tài)響應(yīng)。系統(tǒng)其他參數(shù)為:σ=3,α=0.05,n=60,λ=0.2,Ωxn=Ωyn=0.6,ζ=0.1,ν=68.12,A=300。從時域響應(yīng)圖5(a)中可知轉(zhuǎn)子的y方向相對位移變化不大;從頻譜響應(yīng)圖5(b)中可知出現(xiàn)了橫蕩頻率f0和工頻f1及橫蕩和工頻的組合頻率2(f0+f1),但是工頻f1占主要成分,說明橫蕩運動對轉(zhuǎn)子的動力學(xué)特性影響不大;如圖5(c)(d)所示,轉(zhuǎn)子的軸心軌跡被限制在橢圓區(qū)域內(nèi),龐加萊截面為一閉合的曲線;系統(tǒng)的最大Lyapunov指數(shù)為-0.058,綜合判斷此時系統(tǒng)做擬周期運動。

圖5 Ω=1.07時系統(tǒng)動力學(xué)穩(wěn)態(tài)響應(yīng)

圖6所示為其他參數(shù)不變,量綱一轉(zhuǎn)子轉(zhuǎn)速Ω=2.45時系統(tǒng)的動力學(xué)穩(wěn)態(tài)響應(yīng)。隨著轉(zhuǎn)子轉(zhuǎn)速的升高,轉(zhuǎn)子y方向的相對位移整體明顯增大;頻譜響應(yīng)圖6(b)與圖5(b)中相比出現(xiàn)f1/2頻率及組合頻率,而且f1/2頻率占主要成分,表明系統(tǒng)受到非線性油膜力的影響較大;轉(zhuǎn)子的軸心軌跡為“香蕉狀”,其運動范圍較圖5(c)有所增加;龐加萊截面表現(xiàn)為兩堆點集構(gòu)成的曲線,在該參數(shù)下系統(tǒng)最大Lyapunov指數(shù)為-0.022 1,表明系統(tǒng)仍然是擬周期運動。

圖6 Ω=2.45時系統(tǒng)動力學(xué)穩(wěn)態(tài)響應(yīng)

圖7所示為其他參數(shù)不變,量綱一轉(zhuǎn)子轉(zhuǎn)速Ω=3.52時系統(tǒng)的動力學(xué)穩(wěn)態(tài)響應(yīng)。時域響應(yīng)圖7(a)與圖5(a)相比位移大小變化不大,但其振幅變化比較劇烈;頻譜響應(yīng)圖7(b)與圖6(b)、圖5(b)相比,出現(xiàn)2(f0+f1)/5組合頻率,并且此組合頻率占主要成分,說明此時非線性油膜力和橫蕩慣性力共同作用使得系統(tǒng)的運動更為復(fù)雜;轉(zhuǎn)子的軸心軌跡表現(xiàn)為在橢圓區(qū)域內(nèi)的振蕩,龐加萊截面為一定區(qū)域內(nèi)的點集;系統(tǒng)的最大Lyapunov指數(shù)為0.000 10,綜合判斷系統(tǒng)在此參數(shù)下處在混沌運動狀態(tài)。

圖7 Ω=3.52時系統(tǒng)動力學(xué)穩(wěn)態(tài)響應(yīng)

2.3 橫蕩頻率對系統(tǒng)動力學(xué)特性的影響

圖8所示為量綱一橫蕩頻率ν=60時系統(tǒng)的動力學(xué)穩(wěn)態(tài)響應(yīng)。系統(tǒng)其他參數(shù)為:Ω=2.15,σ=3,α=0.05,n=60,Ωxn=Ωyn=0.6,λ=0.2,ζ=0.1,A=300。由時域響應(yīng)圖8(a)可以看出此時轉(zhuǎn)子的相對位移變化比較劇烈,從頻譜圖8(b)中可知出現(xiàn)了橫蕩運動頻率f0和工頻f1,但是相比工頻f1,橫蕩頻率f0占比較大,表明此時系統(tǒng)受到橫蕩的影響比較大;龐加萊截面表現(xiàn)為“帶狀”有規(guī)律的點,此時的系統(tǒng)處在擬周期運動狀態(tài)。

圖9所示為其他參數(shù)不變,量綱一橫蕩頻率ν=300時系統(tǒng)的動力學(xué)穩(wěn)態(tài)響應(yīng)。對比圖8(a)可以看到,在該橫蕩頻率下轉(zhuǎn)子的相對位移明顯變小,而且轉(zhuǎn)子的位移變得非常的平穩(wěn);從頻譜圖9(b)中可見,工頻f1起主要作用,說明此時的系統(tǒng)受到橫蕩作用影響很小;轉(zhuǎn)子的軸心軌跡由于偏心質(zhì)量的影響呈橢圓狀,龐加萊截面為一條封閉的曲線;綜合判斷在該參數(shù)下系統(tǒng)的運動狀態(tài)為擬周期。

圖8 ν=60時系統(tǒng)動力學(xué)響應(yīng)

圖9 ν=300時系統(tǒng)動力學(xué)響應(yīng)

2.4 橫蕩幅值對系統(tǒng)動力學(xué)特性的影響

圖10和圖11所示分別為不同橫蕩幅值下轉(zhuǎn)子x、y方向的振幅隨轉(zhuǎn)速的變化。系統(tǒng)參數(shù)為:σ=3,α=0.05,n=60,Ωxn=Ωyn=0.6,λ=0.2,ζ=0.1,ν=136.8。整體上看,在轉(zhuǎn)速低于3.0時轉(zhuǎn)子x、y方向的振幅都會隨著橫蕩幅值的增大而增大。橫蕩幅值一定時且轉(zhuǎn)子轉(zhuǎn)速在0.5~1.12區(qū)間內(nèi),轉(zhuǎn)子的振幅會隨著轉(zhuǎn)子轉(zhuǎn)速的增大而增大,且y方向的振幅要比x方向振幅大;當(dāng)轉(zhuǎn)速為1.12~2.28時,轉(zhuǎn)子的振幅呈下降趨勢,且y方向的振幅下降得比較明顯;轉(zhuǎn)速為2.28~2.6時,轉(zhuǎn)子的振幅會突然增大然后迅速減小;轉(zhuǎn)速為2.6~3.01時轉(zhuǎn)子x方向的振幅變化不明顯,但y方向的振幅呈下降趨勢;當(dāng)轉(zhuǎn)速繼續(xù)增大時,轉(zhuǎn)子的振幅呈不斷上升趨勢。

圖10 不同橫蕩幅值下轉(zhuǎn)子x方向振動幅值隨轉(zhuǎn)速變化

圖11 不同橫蕩幅值下轉(zhuǎn)子y方向振動幅值隨轉(zhuǎn)速變化

圖12所示為量綱一橫蕩幅值A(chǔ)=100時系統(tǒng)的動力學(xué)穩(wěn)態(tài)響應(yīng)。系統(tǒng)其他參數(shù)為:σ=3,α=0.05,n=60,λ=0.2,Ωxn=Ωyn=0.6,ζ=0.1,ν=136.8。可以看到,當(dāng)橫蕩幅值較小時轉(zhuǎn)子y方向的相對位移變化較小,頻譜響應(yīng)圖由橫蕩頻率f0及工頻f1組成,但是橫蕩頻率f0占比較小,說明此時橫蕩對系統(tǒng)的影響較小;軸心軌跡表現(xiàn)為在橢圓區(qū)域內(nèi)的振蕩,龐加萊截面為點集構(gòu)成的曲線,綜合判斷此時的系統(tǒng)運動狀態(tài)為擬周期。

圖12 A=100時系統(tǒng)動力學(xué)穩(wěn)態(tài)響應(yīng)

圖13所示為其他參數(shù)不變,量綱一橫蕩幅值A(chǔ)=500時系統(tǒng)的動力學(xué)穩(wěn)態(tài)響應(yīng)。此時橫蕩幅值較大,對比圖12(a)可以看到,橫蕩幅值變大時,轉(zhuǎn)子y方向的相對位移出現(xiàn)明顯的波動;頻譜響應(yīng)中橫蕩頻率f0占比較圖12(b)中明顯增大,此時轉(zhuǎn)子軸心軌跡在船體橫蕩影響下y方向的運動更加明顯;龐加萊截面仍是有規(guī)律的點集組成的線狀結(jié)構(gòu),表明此時系統(tǒng)仍處在擬周期運動狀態(tài)。

圖13 A=500時系統(tǒng)動力學(xué)穩(wěn)態(tài)響應(yīng)

3 結(jié)論

(1)當(dāng)轉(zhuǎn)子的轉(zhuǎn)速較低時,船用轉(zhuǎn)子-浮筏氣囊系統(tǒng)受到橫蕩作用時會由原來的單周期同步運動變?yōu)閿M周期運動;同時,轉(zhuǎn)子y方向的相對位移也會隨之變大。隨著轉(zhuǎn)子轉(zhuǎn)速的增加,系統(tǒng)出現(xiàn)分岔現(xiàn)象,后又變?yōu)閿M周期直至混沌。

(2)當(dāng)頻率比較小即橫蕩頻率較大時,橫蕩對系統(tǒng)的非線性動力學(xué)影響起主導(dǎo)作用;隨著頻率比的增大,轉(zhuǎn)子的振動幅值減小。

(3)當(dāng)轉(zhuǎn)子轉(zhuǎn)速一定時,橫蕩幅值增大會使得轉(zhuǎn)子y方向的振動幅值增大。

猜你喜歡
系統(tǒng)
Smartflower POP 一體式光伏系統(tǒng)
WJ-700無人機系統(tǒng)
ZC系列無人機遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
基于PowerPC+FPGA顯示系統(tǒng)
基于UG的發(fā)射箱自動化虛擬裝配系統(tǒng)開發(fā)
半沸制皂系統(tǒng)(下)
FAO系統(tǒng)特有功能分析及互聯(lián)互通探討
連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
一德系統(tǒng) 德行天下
PLC在多段調(diào)速系統(tǒng)中的應(yīng)用
主站蜘蛛池模板: 国产精品性| 久久精品视频一| 黄色在线网| 国产精品久久久久婷婷五月| 91精品国产一区自在线拍| 中文字幕人妻av一区二区| 97se亚洲综合| 欧美在线伊人| 国产污视频在线观看| 国产屁屁影院| 日韩无码视频专区| 91精品久久久无码中文字幕vr| 久久无码av三级| 毛片在线播放网址| 国产成人精品2021欧美日韩| 91精品国产自产在线老师啪l| 无遮挡国产高潮视频免费观看| 国产第一页屁屁影院| 亚洲欧美一区二区三区图片| 亚洲无码A视频在线| 国语少妇高潮| 欧美一级高清片欧美国产欧美| 日韩精品成人网页视频在线| 六月婷婷激情综合| 91精品视频网站| 67194亚洲无码| 亚洲不卡影院| 日本精品影院| 色男人的天堂久久综合| 国产福利一区在线| 国产传媒一区二区三区四区五区| 91免费观看视频| 中文字幕在线日韩91| 另类欧美日韩| 亚洲精品无码AV电影在线播放| 免费看a级毛片| 久久无码av一区二区三区| 97se亚洲| 97精品国产高清久久久久蜜芽| 青青操视频在线| 99爱视频精品免视看| 成年人久久黄色网站| 精品国产香蕉伊思人在线| 九九香蕉视频| 国产拍揄自揄精品视频网站| V一区无码内射国产| 婷婷开心中文字幕| 国内精品久久久久久久久久影视| 9啪在线视频| 高清欧美性猛交XXXX黑人猛交| 久久女人网| 午夜影院a级片| 色香蕉影院| 亚洲欧洲国产成人综合不卡| 国产精品开放后亚洲| 91原创视频在线| 91po国产在线精品免费观看| 久久公开视频| 亚洲妓女综合网995久久| 亚洲制服中文字幕一区二区| 97超级碰碰碰碰精品| 亚洲欧洲美色一区二区三区| 欧美中文字幕在线播放| 亚洲 日韩 激情 无码 中出| 亚洲欧洲日本在线| 97在线公开视频| 色播五月婷婷| 国产精品欧美激情| 亚洲无码四虎黄色网站| 一级毛片在线免费视频| 国产成人精品三级| 精品国产三级在线观看| 亚洲精品在线观看91| 内射人妻无码色AV天堂| 国产精品视频公开费视频| 无套av在线| 欧美午夜网站| 一级毛片免费播放视频| 亚洲第一页在线观看| 国产你懂得| 99re视频在线| 日日噜噜夜夜狠狠视频|