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

分數(shù)階單自由度線性振子雙側(cè)對稱碰撞振動分析

2021-09-27 07:05:04劉汝逾牛江川申永軍楊紹普
振動與沖擊 2021年16期
關(guān)鍵詞:振動系統(tǒng)研究

劉汝逾,牛江川,2,申永軍,2,楊紹普,2

(1.石家莊鐵道大學(xué) 機械工程學(xué)院,石家莊 050043;2.石家莊鐵道大學(xué) 省部共建交通工程結(jié)構(gòu)力學(xué)行為與系統(tǒng)安全國家重點實驗室,石家莊 050043)

碰撞振動系統(tǒng)是一種非光滑系統(tǒng),存在許多光滑系統(tǒng)不具備的復(fù)雜現(xiàn)象。由于碰撞振動的復(fù)雜性,吸引了許多學(xué)者對其進行研究。碰撞振動廣泛存在于機械設(shè)備中,例如刨煤機刨頭接觸煤炭時會發(fā)生碰撞[1]和輪船航行中會與楔形物發(fā)生碰撞[2]。呂小紅等[3]通過數(shù)值仿真研究了一類具有雙側(cè)剛性約束的單自由度碰撞振動系統(tǒng),分析了系統(tǒng)由周期運動向混沌轉(zhuǎn)遷的分岔過程。張繼業(yè)等[4]研究了單自由度自治系統(tǒng)的碰撞振動,得到了碰撞振動周期運動存在的條件及周期運動穩(wěn)定的充分性判據(jù)。羅冠煒等[5]研究了一類單自由度碰撞振動系統(tǒng)的動力學(xué)行為,描述了系統(tǒng)周期運動的規(guī)律,分析了映射奇異性對系統(tǒng)局部分岔和全局分岔的影響。馮進鈐等[6-7]以典型的Duffing單邊碰撞系統(tǒng)為例,研究了系統(tǒng)的顫振分岔,并對系統(tǒng)中的混沌鞍進行了細致的分析。李得洋等[8]利用數(shù)值仿真方法分析了單自由度含對稱彈性約束碰撞的振動系統(tǒng),得到了在兩參數(shù)平面內(nèi)周期運動的分布及轉(zhuǎn)遷規(guī)律。樂源等[9]討論了一類單自由度雙面碰撞振子的對稱型周期n-2運動以及非對稱型周期n-2運動,并通過分析對稱系統(tǒng)的Poincaré映射的對稱性證明了對稱型周期運動只能發(fā)生音叉分岔。Andreaus等[10]建立了含有雙側(cè)間隙和保險杠的單自由度雙側(cè)碰撞模型,分析了諧波激勵強度、間隙寬度、保險杠剛度和偽共振頻率等最重要的參數(shù)對單自由度系統(tǒng)響應(yīng)的影響。Kumar等[11]研究了具有剛性勢壘的Duffing-van der Pol振蕩器在白噪聲激勵下的分岔問題,利用非光滑變量變換將非光滑碰撞振動系統(tǒng)映射到連續(xù)相平面上,并考慮了兩側(cè)剛性約束不對稱布置的情況。Reboucas等[12]使用逐點映射和標(biāo)準平均與非平滑變換相結(jié)合,分析了具有恢復(fù)系數(shù)的單自由度模型的振動沖擊響應(yīng)。Xie等[13]通過隨機平均和Mellin變換的方法研究了非線性振動系統(tǒng)在高斯白噪聲激勵下的近似瞬態(tài)響應(yīng)。Li[14]針對沖擊振動系統(tǒng)的不連續(xù)性,引入非光滑變換克服不連續(xù)性,然后將隨機平均技術(shù)應(yīng)用于變換后的隨機變質(zhì)量系統(tǒng),通過求解Fokker-Planck-Kolmogorov方程,可以得到近似的解析響應(yīng),以Van der Pol變質(zhì)量振動沖擊系統(tǒng)為例,說明了該方法的有效性。Rigaud等[15]提出了由齒輪嚙合的振動沖擊引起齒輪嘎嘎聲并進行了試驗研究,設(shè)計了一個特定的試驗裝置分析正齒輪在輸入速度波動下的非線性動力學(xué)行為。Reboucas等[16]采用諧波線性化、平均化和數(shù)值模擬的方法,對單自由度振動沖擊振蕩器的頻率響應(yīng)進行了分析,并考慮了3種不同的沖擊力模型。Chen等[17]針對高斯白噪聲激勵下單自由度振動沖擊系統(tǒng)響應(yīng)的封閉形式平穩(wěn)概率密度函數(shù),提出了一種新的求解方法,采用Zhuravlev-Ivanov變換將單側(cè)勢壘的振動沖擊振蕩器轉(zhuǎn)換為無勢壘的振蕩器。Nguyen等[18]采用數(shù)學(xué)建模和動力學(xué)分析的方法,研究了慣性質(zhì)量和激勵頻率對自成移動系統(tǒng)的Duffing沖擊振動裝置的耦合影響。Yan等[19]研究了單側(cè)和雙側(cè)軟約束振動沖擊膠囊系統(tǒng)在質(zhì)量比、剛度比、接觸間隙、外激勵幅值和頻率等控制參數(shù)變化下的動力學(xué)特性,得到單側(cè)約束膠囊的動力學(xué)行為主要是周期性的,當(dāng)左側(cè)約束剛度增大時,具有兩側(cè)約束的其他兩個膠囊的動力學(xué)響應(yīng)變得復(fù)雜的結(jié)論。Xu等[20]研究了具有單邊非零偏移屏障的非線性振動系統(tǒng)的多值響應(yīng)和動力學(xué)特性,發(fā)現(xiàn)在某些情況下,沖擊系統(tǒng)可能具有兩個或4個穩(wěn)態(tài)解。然后討論了穩(wěn)態(tài)解的分類,結(jié)果表明非平凡的穩(wěn)態(tài)解可能會因鞍結(jié)分岔和Hopf分岔而失去穩(wěn)定性。Xu等[21]給出了一般非線性振動沖擊振子的Melnikov方法,以Duffing振動沖擊振蕩器為例說明了該方法的應(yīng)用,并通過相圖、Poincaré截面和分岔圖對所得結(jié)果進行了驗證。

由于分數(shù)階微積分可以提供一個強有力的工具對工程中的黏彈性等行為進行建模,近年來關(guān)于分數(shù)階微積分的研究受到越來越多的關(guān)注。廖少鍇等[22]研究了非線性分數(shù)微積分求解的單步數(shù)值積分算法,討論了振子自由振動及強迫振動下參數(shù)變化對振子非線性特性的影響。楊建華等[23]利用諧波平衡法得到了含分數(shù)階導(dǎo)數(shù)阻尼的一類線性系統(tǒng)在簡諧激勵下系統(tǒng)響應(yīng)的近似解。黃燦等[24]利用分數(shù)階導(dǎo)數(shù)算子研究了線性分數(shù)階振動系統(tǒng)在諧波激勵下的穩(wěn)態(tài)響應(yīng),并討論了分數(shù)階導(dǎo)數(shù)項對剛度和阻尼的影響。銀花等[25]基于精細積分方法,提出了一種新的具有分數(shù)階導(dǎo)數(shù)型本構(gòu)關(guān)系的黏彈性結(jié)構(gòu)動力響應(yīng)的數(shù)值計算方法。王學(xué)彬[26]給出了兩種常見分數(shù)階導(dǎo)數(shù)即Riemann-Liouille分數(shù)階導(dǎo)數(shù)和Caputo分數(shù)階導(dǎo)數(shù)的拉普拉斯變換公式,并給出了具體實例。閆啟方等[27]運用黏彈性理論、分數(shù)階導(dǎo)數(shù)的性質(zhì),研究了分數(shù)階Kelvin模型描述的黏彈性材料的松弛特性、蠕變特性和動態(tài)力學(xué)行為。劉艷芹[28]將分數(shù)階微分算子引入到黏彈性介質(zhì)中的阻尼振動中建立了分數(shù)階非線性振動方程,并利用Adomian分解方法借助Mathematica軟件求解了該類分數(shù)階阻尼振動方程的近似解,研究了振子運動與方程中分數(shù)階導(dǎo)數(shù)的關(guān)系。牛江川等[29]研究了基于速度反饋分數(shù)階PID(proportional integral derivative)控制的單自由度線性振子的自由振動,利用平均法得到了系統(tǒng)的近似解析解。Yang等[30]研究了高斯白噪聲激勵下兩類分數(shù)階導(dǎo)數(shù)沖擊振蕩器的模型,利用非光滑變換和隨機平均方法得到解析近似解,并對其隨機分岔進行了討論。

關(guān)于分數(shù)階黏彈性材料的應(yīng)用和碰撞振動系統(tǒng)的研究成果已經(jīng)非常豐富,但是關(guān)于附加黏彈性材料碰撞振動系統(tǒng)的研究卻很少。將分數(shù)階微積分引入到碰撞振動系統(tǒng)的研究,可以更深入地理解附加黏彈性材料的碰撞振動行為。本文研究了系統(tǒng)關(guān)鍵參數(shù)對系統(tǒng)周期運動的影響,可以為附加黏彈性材料的傳動齒輪、鑄造振動篩和壓力成型機等機械設(shè)備存在的非線性問題提供借鑒,避免這些系統(tǒng)在混沌狀態(tài)下工作,提高機械設(shè)備的動態(tài)性能,實現(xiàn)對復(fù)雜設(shè)備的振動和噪聲的有效控制,為碰撞振動系統(tǒng)的優(yōu)化設(shè)計和改造提供一定的依據(jù)。

本文利用分數(shù)階微分描述黏彈性材料的黏彈性,然后基于分數(shù)階單自由度線性雙側(cè)剛性碰撞模型,研究雙側(cè)對稱碰撞振動系統(tǒng)在簡諧激勵下的穩(wěn)定性和分岔行為。對于碰撞振動的穩(wěn)態(tài)解,利用平均法得到分數(shù)階線性系統(tǒng)的等效剛度和等效阻尼;對于碰撞振動的瞬態(tài)解,利用迭代法得到更精確的瞬態(tài)固有頻率。在此基礎(chǔ)上,分析了雙側(cè)對稱碰撞振動系統(tǒng)的近似解析解。根據(jù)近似解析解,分析了對稱n-1-1周期運動的存在條件,并利用Poincaré映射研究了n-1-1周期運動的穩(wěn)定性。進一步分析了當(dāng)外界激勵頻率、分數(shù)階階次和間隙變化時該系統(tǒng)的分岔行為。

1 近似解析解

分數(shù)階單自由度雙側(cè)碰撞模型,如圖1所示。其中包含剛度系數(shù)為k的線性彈簧、阻尼系數(shù)為c的線性阻尼和系數(shù)為K1(K1>0)的黏彈性材料。振子m受到周期性的激勵力Fcos(ωt+θ)作用,其中振幅為F、激勵頻率為ω、初始相位為θ。

圖1 含間隙的單自由度雙側(cè)碰撞物理模型Fig.1 Single-degree-of-freedom bilateral vibro-impact physical model with gaps

當(dāng)沒有碰撞時的數(shù)學(xué)模型可以表示為

(1)

式中,Dp[x(t)]為x(t)關(guān)于t的p階導(dǎo)數(shù)(0≤p≤1)。這里我們采用Caputo[31]定義

(2)

式中,Γ(z)為Gamma函數(shù),滿足Γ(z+1)=zΓ(z)。

當(dāng)碰撞發(fā)生時滿足以下關(guān)系

(3)

由于存在阻尼和分數(shù)階導(dǎo)數(shù),因此在諧波激勵下系統(tǒng)解需要分別考慮穩(wěn)態(tài)解和瞬態(tài)解。首先假設(shè)系統(tǒng)運動無碰撞,然后研究其穩(wěn)態(tài)解。應(yīng)用平均法得到兩個參數(shù)Ceq(p)和Keq(p),它們表示新的阻尼系數(shù)和新的剛度系數(shù)[32-33]。

Ceq(p)=c+K1ωp-1sin(pπ/2)

(4)

Keq(p)=k+K1ωpcos(pπ/2)

(5)

然后得到式(1)穩(wěn)態(tài)解的一階近似等價整數(shù)階微分方程

(6)

此外,研究式(1)的瞬態(tài)解則無需考慮外部激勵,將運動微分方程重寫為

(7)

類似地,系統(tǒng)自由振動的新阻尼系數(shù)和新剛度系數(shù)可以表示為

(8)

(9)

從而得到式(1)瞬態(tài)解的一階近似等價整數(shù)階方程

(10)

(11)

(12)

(13)

(14)

(15)

(16)

系統(tǒng)的速度可以由式(15)和式(16)求導(dǎo)得到。

2 數(shù)值解驗證

圖2 近似解析解和數(shù)值解擬合Fig.2 Approximate analytical solution and numerical solution fitting

在不考慮碰撞的情況下,將數(shù)值解與解析解進行比較。結(jié)果表明,近似等效整數(shù)階系統(tǒng)的穩(wěn)態(tài)解具有較高的精度,而瞬態(tài)解的精度較低。另外,分數(shù)階系數(shù)K1越大,瞬態(tài)解的擬合誤差越大。從碰撞振動系統(tǒng)的特征可以看出,碰撞振動系統(tǒng)的每個部分都是瞬態(tài)運動。但是,當(dāng)K1相對較小時,誤差要小得多。考慮到振動的影響,我們只需改變上述參數(shù)K1即可分析系統(tǒng)的擬合誤差。當(dāng)K1=0.5,K1=1.0和K1=1.5時,系統(tǒng)數(shù)值解和近似解析解的局部放大擬合圖,如圖3~圖5所示。隨著K1的增大,瞬態(tài)解的擬合誤差將會增大。

圖3 K1=0.5的瞬態(tài)解的擬合誤差Fig.3 Fitting error with K1=0.5

圖4 K1=1.0的瞬態(tài)解的擬合誤差Fig.4 Fitting error with K1=1.0

圖5 K1=1.5的瞬態(tài)解的擬合誤差Fig.5 Fitting error with K1=1.5

3 n-1-1周期運動的存在性分析

在合適的參數(shù)下,系統(tǒng)可以表現(xiàn)為n-p-q周期運動。其中:n為外界激勵周期數(shù);p和q分別為振子與剛性約束W1和W2碰撞的次數(shù)。撞擊面被視為Poincaré截面,并表示為

(17)

假設(shè)振子與剛性約束W1之間的第一次碰撞發(fā)生在時刻t=t1,則振子移動到剛性約束W2處并在時刻t=t2與W2發(fā)生碰撞,具體的流程如圖6所示。根據(jù)圖6所示的周期性n-1-1運動的特性,系統(tǒng)必須滿足

圖6 對稱周期運動相圖Fig.6 Phase portrait of symmetric periodic motion

(18)

圖的對稱周期運動相圖Fig.7 Phase portrait of symmetric periodic motion with

4 n-1-1周期運動的Poincaré映射及其穩(wěn)定性

本章利用Poincaré映射研究n-1-1周期運動的穩(wěn)定性。根據(jù)式(1),等效平面系統(tǒng)可以表示為

(19)

把式(19)的二維系統(tǒng)改成三維系統(tǒng),可以得到

(20)

同時將Poincaré橫截面定義為

(21)

根據(jù)雙側(cè)對稱碰撞的特征,碰撞振子的運動軌跡在一個周期內(nèi)分為4個階段:①碰撞振子離開剛性約束W1并移動到剛性約束W2;②碰撞振子與剛性約束W2發(fā)生碰撞;③碰撞振子離開剛性約束W2并移動到剛性約束W1;④碰撞振子與剛性約束W1發(fā)生碰撞。

因此,可以將這4個映射表示為

(22)

根據(jù)式(13)和式(14)并令τ(0)=θ-φ,可以得到Poincaré映射為

(23)

假設(shè)DP1,DP2,DP3和DP4分別對應(yīng)4個Poincaré映射的連續(xù)階段的線性化矩陣,則可以有

(24)

當(dāng)碰撞振子與剛性約束碰撞時,只改變速度而不改變位移,且速度恢復(fù)系數(shù)為r,可以得到

(25)

則與P1光滑矢量場相關(guān)的局部流可表示為

(26)

(27)

映射關(guān)系可以建立為

(28)

(29)

為了進行以下計算,可以將隱式函數(shù)式(29)表示為

(30)

(31)

因此,可以得到DP1為

DP1=

(32)

DP3=

(33)

然后可以得到碰撞振動系統(tǒng)的線性化矩陣為

DP=DP4·DP3·DP2·DP1

(34)

應(yīng)用DP的特征值可以判斷碰撞振動系統(tǒng)周期運動的穩(wěn)定性。當(dāng)兩個特征值都位于單位圓內(nèi)時,系統(tǒng)的周期運動是穩(wěn)定的。只要一個特征值不處于單位元內(nèi),周期運動在不動點就不穩(wěn)定。根據(jù)第3章中的參數(shù),可得兩個特征值為λ1=0.313 020和λ2=0.424 067,均位于單位元內(nèi)。可以得出結(jié)論:其n-1-1周期運動在不動點處是穩(wěn)定的,也可以從圖7中得到驗證。

5 分岔分析與討論

激勵頻率ω對分數(shù)階碰撞振動系統(tǒng)的動力學(xué)行為有很大的影響。首先取第2章的系統(tǒng)參數(shù),僅更改K1=1.5。首先選取激勵頻率ω作為分岔變量,隨著激勵頻率的變化,根據(jù)數(shù)值解可以得到系統(tǒng)在碰撞面上的分岔圖,如圖8所示。可以看出,當(dāng)激勵頻率ω變化時,系統(tǒng)發(fā)生不連續(xù)擦邊分岔。當(dāng)ω為3.000~1.000時,系統(tǒng)表現(xiàn)為1-1-1周期運動、1-2-2周期運動、1-3-3周期運動、1-4-4周期運動,然后是黏滑振動。

圖8 外部激勵頻率ω的分岔圖Fig.8 Bifurcation diagram with ω

當(dāng)ω=2.682和ω=2.683時的相圖,如圖9所示。結(jié)果表明,系統(tǒng)從1-1-1周期運動經(jīng)過不連續(xù)擦邊分岔為1-2-2周期運動,而且在分岔前后都是穩(wěn)定的周期運動。ω=2.682和ω=2.683的位移時間歷程,如圖10和圖11所示。當(dāng)ω=1.780和ω=1.420時同樣也發(fā)生不連續(xù)擦邊分岔。

圖9 不連續(xù)擦邊分岔的相圖Fig.9 Phase portrait of discontinuous boundary bifurcation

圖10 ω=2.682的位移時間歷程Fig.10 Displacement time history of ω=2.682

圖11 ω=2.683的位移時間歷程Fig.11 Displacement time history of ω=2.683

圖12 分數(shù)階階次p的分岔圖Fig.12 Bifurcation diagram of fractional order p

當(dāng)p介于0.390~1.000時,系統(tǒng)表現(xiàn)為對稱的1-1-1周期運動。例如:當(dāng)p=0.800時,系統(tǒng)表現(xiàn)為對稱的1-1-1周期運動,其相圖如圖14所示。當(dāng)p=0.300時,系統(tǒng)表現(xiàn)為非對稱的1-1-1周期運動,其相圖如圖13所示。當(dāng)系統(tǒng)參數(shù)改變時,系統(tǒng)會發(fā)生擦邊分岔。在上述參數(shù)給定的情況下,分岔的分數(shù)階臨界值為p=0.276,如圖15所示。

圖13 p=0.800的對稱周期運動相圖Fig.13 Phase portrait of symmetric periodic motion with p=0.800

圖14 p=0.300的非對稱周期運動相圖Fig.14 Phase portrait of asymmetric periodic motion with p=0.300

圖15 p=0.276的擦邊分岔相圖Fig.15 Phase portrait of grazing bifurcation with p=0.276

當(dāng)p<0.276時,系統(tǒng)從非對稱周期1-1-1運動轉(zhuǎn)變?yōu)閿M周期運動。隨著p減小,系統(tǒng)表現(xiàn)為混沌運動。當(dāng)p=0.100時,可以得到系統(tǒng)的混沌吸引子,如圖16所示。

圖16 p=0.100的混沌運動Fig.16 Chaotic motion with p=0.100

圖17 關(guān)于碰撞間隙Δ的分岔圖Fig.17 Bifurcation diagram of collision gap Δ

圖18 Δ=0.850的對稱周期運動相圖Fig.18 Phase portrait of symmetric periodic motion with Δ=0.850

圖19 θ=0,Δ=0.850的非對稱周期運動相圖Fig.19 Phase portrait of asymmetric periodic motion with θ=0,Δ=0.850

圖20 θ=π,Δ=0.850的非對稱周期運動相圖Fig.20 Phase portrait of asymmetric periodic motion with θ=π,Δ=0.850

圖21 Δ=0.630的擬周期運動相圖Fig.21 Phase portrait of quasiperiodic motion with Δ=0.630

6 結(jié) 論

本文研究了具有兩個剛性約束的分數(shù)階線性碰撞振動系統(tǒng)在正弦激勵下的穩(wěn)定性和分岔行為。通過引入分數(shù)階系統(tǒng)的近似等價整數(shù)階系統(tǒng),得到了系統(tǒng)的近似解析解。基于近似解析解得到了雙側(cè)對稱n-1-1周期運動的存在條件,并用Poincaré映射分析了n-1-1周期運動的穩(wěn)定性。通過改變外激勵頻率、分數(shù)階階次和間隙,利用數(shù)值解對雙側(cè)對稱碰撞振動系統(tǒng)的分岔行為進行了深入的分析。結(jié)果表明,當(dāng)激勵頻率變化時,系統(tǒng)會出現(xiàn)不連續(xù)擦邊分岔,且系統(tǒng)在分岔前后始終是穩(wěn)定的周期運動。此外,在另一組參數(shù)下,隨著分數(shù)階階次的減小,出現(xiàn)了音叉分岔和擦邊分岔。另外,當(dāng)以間隙為控制參數(shù)時,隨著間隙的減小,對稱不動點變得不穩(wěn)定,而且會產(chǎn)生一對反對稱倍周期序列,將導(dǎo)致擬周期運動和混沌。該模型的計算過程和結(jié)論可為分數(shù)階碰撞振動系統(tǒng)的動力學(xué)行為研究提供參考。

猜你喜歡
振動系統(tǒng)研究
振動的思考
Smartflower POP 一體式光伏系統(tǒng)
FMS與YBT相關(guān)性的實證研究
遼代千人邑研究述論
WJ-700無人機系統(tǒng)
ZC系列無人機遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
振動與頻率
視錯覺在平面設(shè)計中的應(yīng)用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統(tǒng)研究
中立型Emden-Fowler微分方程的振動性
主站蜘蛛池模板: 国产原创第一页在线观看| 欧美日本不卡| 国产h视频在线观看视频| 欧美日韩国产在线人成app| 婷婷激情亚洲| 一级爆乳无码av| 欧美激情网址| 好久久免费视频高清| 亚洲国产精品日韩av专区| 青青青伊人色综合久久| 夜精品a一区二区三区| 嫩草在线视频| 国产精品永久不卡免费视频| 国产精品无码AⅤ在线观看播放| 国产精品视频猛进猛出| 国产精品香蕉在线观看不卡| 一区二区三区四区在线| a级毛片免费播放| 国产人人乐人人爱| 三上悠亚精品二区在线观看| 综合天天色| 亚洲国产黄色| 免费又黄又爽又猛大片午夜| 久久久精品久久久久三级| 一本一本大道香蕉久在线播放| 国产视频大全| 久久九九热视频| 国产95在线 | 欧美国产综合色视频| 日韩精品中文字幕一区三区| 日本免费福利视频| 国产成人一区| 国产成人无码播放| 找国产毛片看| 99热这里只有精品在线播放| 高清久久精品亚洲日韩Av| 精品国产福利在线| 成人一级免费视频| 国产成人久久综合777777麻豆| 国产日韩久久久久无码精品| 国产高清在线观看91精品| 亚洲黄色片免费看| 欧美啪啪网| 国产精品无码一二三视频| 日本在线免费网站| 亚洲第一av网站| 97视频免费在线观看| 日韩亚洲综合在线| 国产人人射| 99九九成人免费视频精品| 国产日本视频91| 日韩一区二区在线电影| 精品国产免费人成在线观看| 欧美日韩亚洲综合在线观看 | 亚洲区欧美区| 日韩视频免费| 亚洲男人的天堂久久精品| 国产精品va| 欧美性天天| 五月天福利视频| 一区二区三区成人| 午夜限制老子影院888| 尤物特级无码毛片免费| 72种姿势欧美久久久大黄蕉| 亚洲va在线∨a天堂va欧美va| 亚洲一区二区三区国产精品 | 国产国拍精品视频免费看| 色妺妺在线视频喷水| yjizz视频最新网站在线| 亚洲综合激情另类专区| 人妻21p大胆| 露脸国产精品自产在线播| 国产清纯在线一区二区WWW| 91亚洲视频下载| 国产爽爽视频| 伊人久久精品无码麻豆精品| 日本精品视频一区二区| 波多野衣结在线精品二区| 在线观看亚洲人成网站| 国产一区成人| 久久久久免费精品国产| 亚洲综合在线最大成人|