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

移動質量簡支梁耦合時變系統建模與實驗設計*

2015-03-13 02:24:35馬志賽周思達
振動、測試與診斷 2015年5期
關鍵詞:模態實驗質量

馬志賽, 劉 莉, 周思達, 楊 武

(北京理工大學飛行器動力學與控制教育部重點實驗室 北京,100081)

?

移動質量簡支梁耦合時變系統建模與實驗設計*

馬志賽, 劉 莉, 周思達, 楊 武

(北京理工大學飛行器動力學與控制教育部重點實驗室 北京,100081)

建立了移動質量簡支梁耦合時變系統的動力學模型,通過數值仿真分析了移動質量速度及加速度對耦合時變系統模態參數的影響,得到移動質量誘導產生的附加阻尼。設計并搭建移動質量簡支梁實驗系統,通過參考實驗得到實驗系統的初始阻尼,并分別采用頻域和時域模態參數辨識方法對質量塊不同移動速度下的實驗系統進行辨識。結果表明,所建立動力學模型能夠對移動質量問題進行準確描述,實驗系統可為時變結構動力學分析的理論研究提供實驗支持,特別是對時變結構模態參數辨識方法進行實驗驗證。

移動質量; 簡支梁; 時變系統; 實驗系統設計

引 言

實際工程結構的特性都是隨時間變化的,特別是在航空航天、機械和橋梁等應用領域結構的時變特性不能被忽略。例如,火箭發射過程中,燃料快速消耗,火箭系統質量特性不斷變化;飛行器超高速飛行時,氣動加熱會引起材料剛度和阻尼隨時間變化;航天器太陽帆板的展開會引起結構質量與剛度特性發生變化;火車通過鐵路橋時,車橋組合系統結構特性會隨著火車的移動不斷變化[1-2]。為更好地對各類時變結構進行分析與設計,有必要對時變結構的動力學分析方法展開研究,特別是對能夠獲取時變結構真實條件下整體動力學特性的實驗模態分析[3]方法進行研究。Poulimenos等[4]通過一個質量塊在簡支梁上移動的實驗系統驗證了所提出的基于時間序列的辨識方法。Spiridonakos等[5]通過一個可伸縮的機械結構對所提出的基于振動響應的故障診斷方法進行了驗證。Liu等[6]通過一個變長度懸臂梁實驗系統驗證了所提出的基于子空間的時變結構模態參數辨識方法。續秀忠等[7]搭建了變質量懸臂梁實驗系統,通過改變懸臂梁上漏斗內水的質量實現了懸臂梁質量隨時間的變化。充分考慮實驗系統的物理意義、時變特征、結構特性隨時間變化形式及其是否精確可控等因素,筆者以移動質量-梁耦合時變系統為對象開展動力學特性研究與實驗設計。

很多工程結構都會承受移動的質量,如車-橋系統、彈-炮系統、重物-橋吊系統及滑塊-導軌系統等,這些系統均可簡化為移動質量-梁模型。移動質量-梁耦合系統本質上為時變系統,在已公開的文獻中,將移動質量作用下梁結構的動力學響應問題歸納為3類[8-9]:移動載荷、移動質量和移動振子問題。移動載荷問題不考慮移動質量的慣性,王文亮等[10]給出了移動載荷作用下簡支梁動力學響應的求解方法。移動質量問題考慮移動質量的慣性但認為其剛度無限大,Michaltsos等[11-13]分別討論了移動質量慣性、向心力、科氏力、轉動慣性以及加速度對簡支梁動力學響應的影響。移動振子問題同時考慮了移動質量的慣性及其有限剛度,Biondi等[8]建立了考慮移動振子慣性情況下移動振子-梁耦合系統的動力學模型。

借鑒Biondi給出的建模方法,在僅忽略移動質量轉動慣性的前提下,筆者建立了任意激勵下移動質量-梁耦合時變系統的動力學模型,精度更高,適用范圍更廣。通過數值仿真分析了移動質量速度及加速度對耦合時變系統模態參數的影響,得到移動質量誘導產生的附加阻尼。設計并搭建移動質量簡支梁實驗系統,通過參考實驗得到實驗系統的初始阻尼,并分別采用頻域和時域模態參數辨識方法對質量塊不同移動速度下的實驗系統進行了辨識。

1 動力學模型

如圖1為移動質量-梁耦合系統,其振動微分方程為

Q(x,t)+P(x,t)

(1)

圖1 移動質量-梁耦合系統Fig.1 Coupled moving-mass and beam system

其中:L,m,c和EI分別為梁的長度、單位長度質量、黏性阻尼系數和彎曲剛度;y(x,t)為梁上x處在t時刻的撓度;Q(x,t)為橫向載荷;P(x,t)為移動質量對梁的等效作用力。

當移動質量的軸向尺寸d遠小于梁的長度L時,忽略其轉動慣性[12],則移動質量對梁的等效作用力可寫為

χ(s(t))δ(x-s(t))

(2)

其中:M0為移動質量的質量;g為重力加速度;s(t)為移動質量的位移函數;δ(x-s(t))為Diracdelta函數。

χ(s(t))為窗函數,定義為

(3)

對于梁的撓度有

(4)

式(4)中各項依次表示牽連加速度、加速度法向分量、科氏加速度和向心加速度。

將式(2)和式(4)代入式(1)可得

(5)

其中:符號′表示對x的微分;符號·表示對t的微分。

設梁的撓度函數具有如下形式

(6)

其中:φi(x)(i=1,2,…,N)為梁的第i階模態振型函數。

梁的撓度函數滿足如下正交關系

(7)

將式(6)代入式(5)可得

(8)

將式(8)兩端同時乘以φj(x),并沿[0,L]積分,利用模態振型函數的正交性可得

(9)

進一步,移動質量-梁耦合系統振動微分方程的矩陣形式可寫為

K(t){q(t)}={F(t)}

(10)

其中

M(t)=diag{Mi}+M0diag{φi(s)}Φ(s)

(11)

其中:Φ(s)為梁的模態振型函數矩陣;Φ′(s)和Φ″(s)分別為Φ(s)關于移動質量位移s的一階和二階導數矩陣。

由式(10)可知,移動質量-梁耦合系統為時變系統。

需要注意的是,上述推導過程中對梁的邊界條件未加任何限制,只要梁的模態振型函數已知就可以應用該動力學模型進行分析,且在移動質量對梁作用力的等效過程中,僅忽略了移動質量的轉動慣性,建模精度更高,適用范圍更廣。對于兩端簡支梁,模態振型函數為

(12)

2 數值仿真

移動質量-梁耦合時變系統的質量、阻尼和剛度矩陣均隨移動質量的位置變化而變化,而移動質量的速度會影響阻尼和剛度矩陣,加速度只對剛度矩陣有影響。因此,通過數值仿真可分析移動質量運動形式對耦合時變系統模態參數的影響。

對于移動質量簡支梁耦合時變系統,梁的長度為L=2 m,單位長度質量為m=4.71 kg/m,彎曲剛度為EI=1 050 N·m2、粘性阻尼系數假設為c=0。移動質量塊的質量為M0=4.865 8 kg,運動形式為勻變速,即s(t)=v0t+at2/2,其中,v0為初速度,a為加速度,重力加速度取為g=9.8 m/s2。

2.1 速度對模態參數的影響

在加速度a=0,初速度分別為v01=0.05 m/s,v02=0.10 m/s和v03=0.20 m/s的情況下對質量塊從簡支梁一端運動到另一端的過程進行仿真,仿真時間分別為40,20和10 s。三種運動形式下,移動質量速度對耦合時變系統前4階模態參數的影響如圖2所示(橫坐標表示質量塊在簡支梁上的位置)。

由圖2可知, 移動質量的速度對耦合時變系統固有頻率的影響很小,對模態阻尼比影響較大,且影響程度隨移動質量速度的增大而增大。由于簡支梁兩端邊界條件的對稱性,移動質量勻速通過簡支梁的過程中,耦合時變系統模態參數的變化規律也具有對稱性。

2.2 加速度對模態參數的影響

在初速度v0=0.05 m/s,加速度分別為a1=0,a2=0.005 m/s2和a3=0.03 m/s2的情況下對質量塊從簡支梁一端運動到另一端的過程進行仿真,仿真時間分別為40,20和10s,在3種運動形式下,移動質量加速度對耦合時變系統前4階模態參數的影響如圖3所示(橫坐標表示質量塊在簡支梁上的位置)。

圖2 移動質量速度對模態參數的影響Fig.2 Modal parameters influenced by moving-mass velocity

圖3 移動質量加速度對模態參數的影響Fig.3 Modal parameters influenced by moving-mass acceleration

由圖3可知,移動質量的加速度對耦合時變系統固有頻率的影響很小,對模態阻尼比影響較大,且影響程度隨移動質量加速度的增大而增大。需要說明的是,當移動質量作勻加速運動時,其速度隨時間的增加而增大,因而,模態參數的變化會受到速度與加速度的綜合影響,且變化規律不再具有對稱性。

綜上可知,移動質量簡支梁耦合時變系統的固有頻率受移動質量運動速度及加速度的影響很小,模態阻尼比對移動質量的速度及加速度則較為敏感。在數值仿真中,系統自身的粘性阻尼系數均假設為c=0,因此,仿真得到的阻尼可認為是由移動質量誘導產生的附加阻尼。由圖2和圖3可知,附加阻尼可能為負值,因而,當系統自身阻尼較小時,附加阻尼效應可能導致系統的某一階或某幾階模態阻尼比為負值。

3 實驗系統設計

移動質量簡支梁實驗系統包括時變結構、激勵系統、力和運動傳感器、測量與分析系統、電機系統和支撐臺體6部分,其設計方案及實物如圖4所示。

圖4 移動質量簡支梁實驗系統設計方案及實物圖Fig.4 Experimental system of the moving-mass and simply supported beam

時變結構包括簡支梁和移動質量塊,通過改變質量塊的質量大小、運動形式(位置、速度和加速度)等方式來模擬不同時變特性的時變結構。簡支梁為2 000mm×60mm×10mm(長×寬×高)的鋼質梁,移動質量塊的質量變化范圍為2.080 0~4.865 8kg。

激勵系統包括激振器和功率放大器,激振器用于輸出不同類型的帶寬信號或者單頻信號,功率放大器用于對輸出信號的幅值進行調節。實驗系統中,激振器選用ModalShop2025E輕型模態激振器,功率放大器選用SmartAmpTM的2100E21-400功率放大器。

力和運動傳感器分別用來測量時變結構所受到的激勵信號和該激勵作用下結構的動力學響應信號。力傳感器選用PCBTM的288D01機械阻抗傳感器,運動傳感器選用PCBTM的333B30加速度傳感器。

測量與分析系統選用LMSTM的SCADASIII系統,可用于預設力和運動傳感器的采樣頻率和采樣時間,并對傳感器采集到的激勵信號和結構動力學響應信號進行處理。

電機系統包括電機、編碼器、運動控制器和同步繞線器,分別選用FaulhaberTM的2642W024CR直流微電機、IE2-512電磁編碼器、MCDC3003S運動控制器和塑制同步繞線器。編碼器和運動控制器用來控制電機輸出端的運動形式(位置、速度和加速度),電機輸出端直接與同步繞線器固連,以繞線的方式驅動質量塊在簡支梁上移動。電機、編碼器和運動控制器構成閉合回路,能夠對電機輸出端的運動形式進行精確控制,進而實現對質量塊運動形式(位置、速度和加速度)的精確控制。此外,通過電機系統可預設力和運動傳感器采樣開始時刻質量塊所需滿足的運動信息標準,移動過程中,當質量塊的運動信息滿足該預設運動信息標準時,測量與分析系統啟動力和運動傳感器開始采樣,采樣終止時刻則通過測量與分析系統直接預設,實現了質量塊任意運動形式下任意時間段內激勵信號和結構動力學響應信號的采集。

支撐臺體由金屬型材搭建而成,為時變結構中的簡支梁提供兩端簡支的邊界條件。

4 實驗結果

實驗前,在簡支梁下表面沿軸向均勻布置15個傳感器,自左向右依次編號為1~15,如圖4(a)所示。5號傳感器位置為激勵的施加位置,選用機械阻抗傳感器同時測量該處的激勵力和加速度響應信號;其余14個傳感器均為加速度傳感器,用以測量所在位置的加速度響應信號。

4.1 參考實驗

對于移動質量簡支梁實驗系統,把從簡支梁中點到其右側與中點相距0.8m的區間做80等分,自左向右分別將質量塊置放于等分點處,依次進行81組時不變結構動力學實驗。實驗過程中,M0=4.865 8 kg,激勵為0~256Hz的隨機激勵,激勵力和加速度響應信號的采樣頻率均為512Hz,采樣時間為4s,選用H1估計[3](平均次數:30次)得到所有響應點與參考點(5號點)之間頻響函數隨質量塊位置的變化規律,其中,5號響應點的原點頻響函數如圖5所示。

圖5 原點頻響函數隨質量塊位置的變化規律Fig.5 Driving-point frequency response function varying with mass position

圖6 固有頻率隨質量塊位置的變化規律Fig.6 Natural frequencies varying with mass position

采用最小二乘復頻域法[3](leastsquarescomplexfrequencydomain,簡稱LSCFD)辨識得到實驗系統前4階固有頻率隨質量塊位置的變化規律如圖6所示(圓圈表示辨識結果, 實線表示只考慮實驗系統質量矩陣受移動質量位移影響而忽略其速度及加速度影響的無阻尼固有頻率仿真結果),前4階模態阻尼比隨質量塊位置的變化規律如圖7所示。

由圖6可知,動力學模型與實驗系統在固有頻率上能夠較好地吻合, 驗證了動力學模型對于移動質量問題的適用性。圖7中阻尼的辨識結果可作為質量塊通過簡支梁過程中耦合時變系統的初始阻尼,綜合考慮初始阻尼與附加阻尼可得到質量塊通過簡支梁過程中實驗系統的真實阻尼。

圖7 模態阻尼比隨質量塊位置的變化規律Fig.7 Modal damping ratio varying with mass position

4.2 時變結構動力學實驗

為對數值仿真結果進行驗證,研究移動質量速度對系統模態參數的影響,進行如下三組時變結構動力學實驗:質量塊的運動形式依次為0.05,0.10和0.20m/s的勻速運動, 對應的采樣時間依次為16,8和4s,每組實驗均進行50次。實驗過程中M0=4.865 8 kg,力和運動傳感器的采樣開始時刻為質量塊通過簡支梁中點的時刻,采樣頻率為512Hz。

4.2.1 頻域辨識結果

由于結構特性隨時間的變化,時變結構的響應信號具有非平穩性,傳統的基于傅里葉變換的譜分析已無法反映非平穩信號頻率分量的時間變換特征,可通過時頻分析[14-15]的方法獲取響應信號隨時間變化的頻譜。采用平滑偽魏格納-維爾分布(smoothedpseudoWigner-Villedistribution[16],簡稱SPWVD)對時變結構的加速度響應信號進行時頻分析,估計出移動質量三種不同速度下各響應點處加速度響應信號時間相關自功率譜的平均值,如圖8所示。由圖8(a~c)可知,加速度響應信號的時間相關自功率譜中的“脊線”能夠清晰地指示出實驗系統固有頻率的變化趨勢。根據加速度響應信號的時間相關自功率譜, 采用基于時頻分布系數的加權最小二乘方法[17]對實驗系統的固有頻率進行辨識,辨識結果如圖8(d~f)所示(圓圈表示辨識結果,實線表示仿真結果)。比較圖8(d~f)中辨識結果與仿真結果可得,所選用頻域方法能夠準確地辨識得到實驗系統的前4階固有頻率。

4.2.2 時域辨識結果

對于移動質量簡支梁實驗系統,如將阻尼比高于8%[3]或負阻尼的極點作為數學極點直接濾掉將導致固有頻率的辨識結果失真。考慮到移動質量誘導產生的附加阻尼,文中將-1%<ξ<5%作為阻尼比估計值的篩選區間,將阻尼比位于區間內的極點作為物理極點加以保留,而將阻尼比位于區間外的極點作為數學極點直接濾掉。

采用泛函序列-多分量時變自回歸滑動平均(functionalseriesvectortime-dependentautoregressivemovingaverage,簡稱FS-VTARMA)[18]方法對實驗系統的固有頻率進行辨識。對于每組實驗,分別基于50次實驗數據進行辨識,將滿足篩選條件的固有頻率辨識結果進行統計,刻畫出實驗系統固有頻率的變化趨勢,如圖8(g~i)所示(圓圈表示辨識結果,虛線表示仿真結果)。比較圖8(g~i)中辨識結果與仿真結果可得,所選用時域方法無法辨識出移動質量簡支梁實驗系統的第1階固有頻率,且辨識精度依賴于采集數據的長度,采集數據越長,辨識精度越高。

5 結 論

1) 在僅忽略移動質量轉動慣性的前提下,建立了任意激勵下移動質量-梁耦合時變系統的動力學模型,精度更高,適用范圍更廣。基于該動力學模型,可分析移動質量速度與加速度對耦合時變系統模態參數的影響,得到移動質量誘導產生的附加阻尼,對實際工程結構的設計具有指導意義。

2) 移動質量簡支梁實驗系統中質量塊的質量大小可調,運動形式精確可控,實現了質量塊任意運動形式下任意時間段內激勵信號和結構動力學響應信號的采集,能夠模擬不同時變特性的時變結構,物理意義明確,時變特征明顯,可用于時變結構模態參數辨識方法的實驗驗證。

圖8 加速度響應信號的時間相關自功率譜及固有頻率辨識結果Fig.8 Time-dependent auto power spectral of acceleration signals and estimated results of natural frequencies

3) 頻域和時域兩種方法辨識結果表明,移動質量簡支梁實驗系統具有良好的可觀測性,系統響應自由度的選擇十分方便,且研究該系統動態特性所需要的全部數據都是可以測量的。

4) 現有已公開的模態參數辨識方法尚不能準確地辨識出隨機激勵下該實驗系統的阻尼特性用于相關建模理論的研究,因此筆者未給出阻尼的辨識結果。然而結合所建動力學模型,提出的實驗系統可為阻尼辨識方法的進一步研究提供實驗驗證。

[1] Au F T K, Jiang R J, Cheung Y K. Parameter identification of vehicles moving on continuous bridges[J]. Journal of Sound and Vibration, 2004, 269(1-2): 91-111.

[2] Marchesiello S, Bedaoui S, Garibaldi L, et al. Time-dependent identification of a bridge-like structure with crossing loads[J]. Mechanical Systems and Signal Processing, 2009, 23(6): 2019-2028.

[3] Heylen W, Lammens S, Sas P. Modal analysis theory and testing[M]. Leuven: Katholieke Universiteit Leuven, 2007: A1.1-A3.34,B4.11.

[4] Poulimenos A G, Fassois S D. Output-only stochastic identification of a time-varying structure via functional series TARMA models[J]. Mechanical Systems and Signal Processing, 2009, 23(4): 1180-1204.

[5] Spiridonakos M D, Fassois S D. An FS-TARbased method for vibration-response-based fault diagnosis in stochastic time-varying structures: experimental application to a pick-and-place mechanism[J]. Mechanical Systems and Signal Processing, 2013, 38(1): 206-222.

[6] Liu K, Deng L. Identification of pseudo-natural frequencies of an axially moving cantilever beam using a subspace-based algorithm[J]. Mechanical Systems and Signal Processing, 2006, 20(1): 94-113.

[7] 續秀忠. 時變線性結構模態參數識別的理論及實驗研究[D]. 上海:上海交通大學, 2003.

[8] Biondi B, Muscolino G. New improved series expansion for solving the moving oscillator problem[J]. Journal of Sound and Vibration, 2005, 281(1-2): 99-117.

[9] 胡紅生, 錢蘇翔. 移動質量激勵下梁動態響應分析與試驗研究[J]. 振動、測試與診斷, 2010, 30(2): 153-157.

Hu Hongsheng, Qian Suxiang. Test and dynamic response analysis of beam under a moving mass[J]. Journal of Vibration, Measurement & Diagnosis, 2010, 30(2): 153-157. (in Chinese)

[10]王文亮, 杜作潤. 結構振動與動態子結構方法[M]. 上海: 復旦大學出版社, 1985:114-119.

[11]Michaltsos G T, Sophianopoulos D, Kounadis A N. The effect of a moving mass and other parameters on the dynamic response of a simply supported beam[J]. Journal of Sound and Vibration, 1996, 191(3): 357-362.

[12]Michaltsos G T. The influence of centripetal and coriolis forces on the dynamic response of light bridges under moving vehicles[J]. Journal of Sound and Vibration, 2001, 247(2): 261-277.

[13]Michaltsos G T. Dynamic behaviour of a single-span beam subjected to loads moving with variable speeds[J]. Journal of Sound and Vibration, 2002, 258(2): 359-372.

[14]Cohen L. Time-frequency analysis[M]. New Jersay: Prentice Hall, 1995:70-81.

[15]Qian S. Introduction to time-frequency and wavelet transforms[M]. New Jersay: Prentice Hall, 2002: 147-198.

[16]Feng Zhipeng, Liang Ming, Chu Fulei. Recent advances in time-trequency analysis methods for machinery fault diagnosis: areview with application examples[J]. Mechanical Systems and Signal Processing, 2013, 38(1): 165-205.

[17]周思達. 線性時變結構時頻域模態參數辨識理論及實驗研究[D]. 北京:北京理工大學, 2012.

[18]Spiridonakos M D, Fassois S D. Parametric identification of a time-varying structure based on vector vibration response measurements[J]. Mechanical Systems and Signal Processing, 2009, 23(6): 2029-2048.

10.16450/j.cnki.issn.1004-6801.2015.05.018

*北京理工大學研究生科技創新專項計劃資助項目;北京理工大學基礎研究基金資助項目(20120142009)

2013-11-18;

2014-03-18

O327

馬志賽,男,1988年6月生,博士生。主要研究方向為時變系統模態參數辨識。 E-mail:mazhisai2008@163.com 通信作者簡介:劉莉,女,1964年5月生,博士、教授。主要研究方向為飛行器總體設計、飛行器結構分析與設計、飛行動力學與控制等。 E-mail:liuli@bit.edu.cn

猜你喜歡
模態實驗質量
記一次有趣的實驗
“質量”知識鞏固
質量守恒定律考什么
做個怪怪長實驗
做夢導致睡眠質量差嗎
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
質量投訴超六成
汽車觀察(2016年3期)2016-02-28 13:16:26
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
主站蜘蛛池模板: 亚洲国产中文在线二区三区免| 伊人久久婷婷| 国产精品第| 麻豆精品国产自产在线| 国产XXXX做受性欧美88| 午夜精品一区二区蜜桃| 国内精品久久人妻无码大片高| 色综合热无码热国产| 国产在线视频欧美亚综合| 国产va免费精品观看| 久久久精品无码一二三区| Jizz国产色系免费| 无码中字出轨中文人妻中文中| 人妻少妇乱子伦精品无码专区毛片| 欧美视频在线不卡| 国产福利小视频在线播放观看| 久久中文字幕不卡一二区| 日韩中文精品亚洲第三区| 91亚洲精品国产自在现线| 亚洲高清在线天堂精品| 日韩区欧美国产区在线观看| 欧美人与性动交a欧美精品| 91精品国产情侣高潮露脸| 噜噜噜久久| 亚洲午夜天堂| 国产无吗一区二区三区在线欢| 丁香六月激情婷婷| 亚洲欧美在线精品一区二区| 九色视频一区| 波多野结衣亚洲一区| 九九久久精品免费观看| 日韩国产欧美精品在线| 69综合网| 成人在线不卡视频| 国产无码高清视频不卡| 青青草国产精品久久久久| 又爽又大又黄a级毛片在线视频 | 欧美午夜在线观看| 91丝袜在线观看| 成年人福利视频| 国产亚洲日韩av在线| 久久精品中文字幕免费| 波多野结衣在线一区二区| 亚洲a级毛片| 婷婷色一二三区波多野衣| 黄色网页在线播放| 久久国产精品无码hdav| 97一区二区在线播放| 亚洲AV人人澡人人双人| 色噜噜狠狠狠综合曰曰曰| 国内精品九九久久久精品| 国产喷水视频| 欧美一级色视频| 久青草免费视频| 国产丰满大乳无码免费播放| 中文字幕亚洲第一| 熟妇无码人妻| 久久精品一品道久久精品| 欧美亚洲欧美| 国产精品浪潮Av| 日韩视频福利| 成年免费在线观看| 国产在线自在拍91精品黑人| 国产99精品视频| 欧美日韩北条麻妃一区二区| 国产免费好大好硬视频| 亚洲中文无码av永久伊人| 91在线精品免费免费播放| 亚洲精品国偷自产在线91正片| 精品免费在线视频| 国产SUV精品一区二区| 91在线日韩在线播放| 国产小视频a在线观看| 免费在线看黄网址| 久久99国产精品成人欧美| 久久精品一品道久久精品| 欧洲一区二区三区无码| 福利片91| 亚洲午夜福利精品无码不卡| 在线a视频免费观看| 国产尤物在线播放| 国产 在线视频无码|