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

盤繞式伸展臂非線性振動力學分析

2018-05-28 02:54:40樊鵬玄陳務軍張祎貝
振動與沖擊 2018年10期
關鍵詞:模態振動結構

樊鵬玄,陳務軍,張祎貝,趙 兵

(上海交通大學 空間結構研究中心,上海 200240)

盤繞式伸展臂作為支撐構件廣泛應用于太陽翼、太陽帆、空間站等航天器[1-2]。具有構造簡捷、比剛度大、收納率高、成本低等優點,是可展結構中高效、可靠的典范。盤繞臂在航天器中的重要作用決定其振動分析與控制是航天工程中的要點。目前針對盤繞臂的研究多集中于結構穩定性與屈曲荷載[3-4]、展開與收納的動靜力學特性[5-6]。振動分析方面,McEachen等[7-9]計算了盤繞臂線性振動頻率并進行了試驗驗證。戈東明、陳務軍等[10-12]進行了盤繞臂展開狀態的線性模態分析。韓建斌等[13]用縮小彈性模量的方法描述三角框受壓屈曲狀態,采用分段線性模型描述加勁索松弛前后的盤繞臂振動力學性能,對加勁索松弛前后的盤繞臂分別進行了模態分析。

盤繞臂屬柔性可展開結構體系,結構剛度受加勁索預應力影響較大。加勁索中的預應力對結構具有加載作用,引發P-Δ效應,要準確獲取預應力對盤繞臂自振頻率、模態的影響須考慮此幾何非線性的影響;航天器在太空中受到沖擊荷載以及盤繞臂展開結束時,盤繞臂將發生加勁索高頻交替張弛的大幅度振動,故考慮大幅振動下狀態非線性因素進行振動分析與控制也是盤繞臂設計的基礎。

以往的報道尚未對預應力加載作用引發P-Δ效應的影響以及大幅振動下加勁索高頻交替張弛引發的狀態非線性進行研究,文章將針對兩個問題分別進行計算。首先采用幾何精確理論及結構屈曲分析理論推導加勁索的預應力取值范圍,然后進行振動分析。小幅振動時,采用增量有限元法考慮預應力加載作用引發幾何非線性對結構模態和頻率的影響,得出小幅振動下使結構自振頻率最大的預應力取值,并繪制設計用“頻率-預應力-直徑”參考面,在參考面內可對加勁索的預應力、直徑等參數進行主動設計;大幅振動時,利用顯式動力學積分解決加勁索狀態非線性中剛度矩陣奇異的問題并獲取振動時程響應,采用自編制的動力時程信號處理程序,提取大幅振動的頻率、模態,并揭示大幅振動時加勁索使盤繞臂振幅迅速衰減的作用機理。可作為大幅非線性振動控制的參考依據。

1 盤繞式伸展臂非線性振動分析方法

以振動中加勁索是否松弛為判定依據,將盤繞臂振動問題為兩類:一是加勁索不發生松弛的小幅振動,二是加勁索發生高頻張弛切換的大幅振動。

1.1 小幅非線性振動

在加勁索預應力的加載作用下盤繞臂形成靜力平衡體系,其振動以此為初始狀態。隨預應力增大,盤繞臂變形增大,預應力加載產生的P-Δ效應增強。可基于非線性有限元法,將P-Δ效應作為荷載剛度納入計算,計算流程如下:

(1) 建立剛度矩陣和質量矩陣

加勁索剛度組成[14]

KC=KLC+KNLC

(1)

縱桿及橫桿(采用梁單元)剛度組成[15-16]

KB=KEB-KPB

(2)

式中:KEB是彈性剛度陣(與EI和l有關);KPB是荷載剛度陣(與作用在單元兩端的力P和l有關)。

依據有限元法剛度矩陣集成規則,集成得整體剛度

(3)

將各個單元質量集中到節點,根據節點編號集成得到集中質量矩陣M。

(2) 建立特征值方程并迭代求解

不考慮系統阻尼,可建立廣義特征值方程式(4),采用子空間迭代提取考慮P-Δ效應的頻率及模態

Kφ-ω2Mφ=0

(4)

式中:φ為固有振型;ω為固有頻率

如式(3),KNLC使盤繞臂剛度增大,KPB使其剛度減小,故存在一個預應力值使盤繞臂自振頻率取極值。取不同預應力值代入式(4),求解結構頻率及模態,可找到使結構頻率最大的預應力。

1.2 大幅非線性振動

大幅振動中,加勁索張緊時盤繞臂是幾何不變結構體系。但松弛時,成為幾何可變體系,剛度矩陣奇異,無法基于式(4)進行特征值迭代。文章采用動力學響應直接積分,獲取盤繞臂自由振動的時程響應,再進行時程信號處理來提取大幅振動參數。

(5)

(6)

求解式及得到時程響應,利用快速傅氏變換(Fast Fourier Transform, FFT)及隨機子空間法[19](SSI/data)編制時程信號處理程序,用于提取自振頻率及模態。后文分析指出,基于時域的SSI/data法雖可同時得到自振頻率及模態,但振動信號中夾雜大量加勁索局部振動帶來的干擾,無法確定可靠的振動參數。利用基于頻域的FFT識別法處理盤繞臂自由端的振動信號,能提取出更準確的頻率值。結合兩種方法提取大幅非線性振動模態及頻率。圖1給出計算流程。

圖1 盤繞式伸展臂大幅度非線性振動計算Fig.1 Nonlinear large amplitude vibration analysis

2 加勁索預應力取值范圍

加勁索預應力對盤繞臂有加載作用(圖2),展開狀態下,設每根加勁索預應力σ0為,加勁索截面積為A,則加勁索作用力F=σ0A。在加勁索預應力作用下縱桿應保持平直,三角框處于受壓的微彎曲狀態。

圖2中BC節段受到4個加勁索力,根據力的矢量和原理將索力投影到縱桿軸線方向及三角框平面內,作用在縱桿上的軸壓力

F1=2Fcosβ

(7)

該力作用下縱桿不發生屈曲

(8)

(9)

式中:t為盤繞臂節距;EIl為縱桿截面抗彎剛度。中間節段縱桿等效于兩端固支桿,端部節段等效于一端固支一端鉸支。

圖2 加勁索作用力Fig.2 The force from the cable

加勁索作用在三角框上的力轉化為Ft,指向形心(圖3)

(10)

圖3 三角框大撓度后屈曲Fig.3 Post-buckling of triangle battern

三角框受壓發生變形,基于幾何精確理論及勢能駐值原理將平衡方程建立在變形后的位形上。總勢能見式(11),平衡方程見式(12)。

(11)

(12)

式中:b是三角框邊長;EIt為橫桿截面抗彎剛度;s為沿橫桿的弧坐標;φ為橫桿變形后切線與原始位置軸線方向的夾角。

求解非線性微分方程式(12)可得到三角框受壓彎曲后變形。以φ近似替代sinφ,可求解出小變形時三角框受壓臨界屈曲荷載

(13)

基于式(7)~式(9)、式(13)可求出加勁索預應力應滿足的條件

(14)

3 盤繞臂振動力學特性計算

3.1 盤繞式伸展臂建模

非線性有限元分析及顯式動力學積分運算過程可直接利用ABAQUS軟件進行,采用以往的盤繞臂設計參數(見表1[20])建立盤繞臂計算模型,盤繞半徑r為150 mm,總長1 884 mm,分為9節(圖4)。

圖4 盤繞式伸展臂整體模型Fig.4 Numerical model of coilable mast

加勁索用多段鉸接桿單元建模,并設置只能承受軸向拉力。

表1 模型參數Tab.1 Material properties

3.2 小幅度非線性振動計算

根據式計算出預應力范圍是:8.60~33.32 MPa,計算時對加勁索施加預應力進行幾何非線性靜力分析以獲取加載后的平衡態,基于該狀態集成彈性剛度與荷載剛度陣,采用子空間迭代法,求解自振模態及頻率。圖5給出預應力從1 MPa到35 MPa的盤繞臂頻譜(增加1 MPa計算1次),可分為3段,各段特性:

A段:加勁索預應力小于等于8 MPa,三角框尚未屈曲,加勁索是剛度薄弱項。圖6給出5 MPa時結構1階模態,整體表現為彎曲。據前文1.1部分的分析,此段預應力增大使KC的增大超過KB的減小,故預應力增大,盤繞臂一階頻率上升。

圖5 預應力與自振頻率曲線Fig.5 Relationship between prestress and frequency

圖6 預應力5 MPa下盤繞臂1階模態Fig.6 First-order mode of coilable mast of 5 MPa prestress

B段:加勁索預應力大于等于9 MPa,三角框進入后屈曲階段。此段KC的增大與KB的減小持平,盤繞臂整體剛度保持水平。圖7給出預應力15 MPa時盤繞臂的1階模態,整體表現為彎曲,但三角框成為盤繞臂剛度弱項。

圖7 預應力15 MPa下盤繞臂1階模態Fig.7 1st mode of coilable mast in 15 MPa prestress

C段:隨預應力增加,P-Δ效應更加顯著,此時KC的增大小于KB的減小,整體剛度下降,自振頻率降低。33 MPa力下結構的1階頻率降低為0(表 2),相應模態(圖8)為縱桿屈曲屈曲模態。這與式(14)的理論推導恰好吻合。

圖8 預應力33 MPa下盤繞臂1階模態Fig.8 1st mode of coilable mast in 33 MPa prestress

表2 結構固有頻率Tab.2 Natural vibration frequency Hz

計算不同直徑加勁索、不同預應力下盤繞臂自振頻率,繪制圖9,可得到設計用“頻率-預應力-直徑”包絡面(以O為頂點的扇形),設計時應使加勁索直徑、預應力落在面中,保證結構剛度最大。

圖9 不同直徑加勁索的盤繞臂頻率曲線Fig.9 Frequency curves of different cable diameters

3.3 大幅度非線性振動計算

3.3.1 動力學響應時程積分

如圖10,固定端板O點,在自由端施加Y向位移δ,使盤繞臂發生變形U。施加位移δ時,靠近固定端的加勁索先松弛,隨δ增大,加勁索松弛的節數增多(表3)。此外,預應力小的加勁索更容易出現松弛;當δ大于或等于5 mm時,9節加勁索全部松弛。

圖10 振動前狀態Fig.10 Deformation before vibration

令δ等于5 mm,使盤繞臂自由振動并進行動力學響應時程積分,繪制自由端Y向位移U-t曲線于圖11,盤繞臂的振幅從5 mm開始衰減,最后穩定在1.4 mm左右。

表3 不同δ下松弛的節數Tab.3 Number of relaxation steps with various δ

圖11 盤繞式伸展臂自由端動力學響應Fig.11 Dynamic response of free-tip

3.3.2 頻率提取及模態識別

取各節三角框端部節點的位移時程數據用于計算結構模態,結構變形主要在y,z方向,x向變形很小,為提高計算效率,取各節點y,z向位移響應進行模態識別。將3.3.1中仿真得到的數據代入基于MATLAB平臺編寫的SSI/data程序中計算,采用頻率、阻尼及振型三重穩定判定準則來判定第j階模態是否穩定,設置容差見式(15)。

(15)

式中:fi,j、ξi,j分別表示選取系統階次為i時,識別的第j階模態的頻率和阻尼比,MACj表示第j階模態的置信準則(Modal Assurance Criterion,MAC),定義見式(16)。

(16)

式中:ψi,j表示當選取系統階次為i時,識別的第j階模態的振型。

取預應力為10 MPa的進行計算,將識別的系統階次定為2 940~3 000。取前240階模態(頻率<40 Hz)繪制穩定圖12。圖中模態十分密集,大部分是加勁索局部模態,無法從圖中直接確定主結構模態。故借助FFT方法得到頻響曲線來判定所需提取的結構模態。

基于圖11的位移時程信號,利用FFT處理得到盤繞臂自振頻率見表4(表中“小幅振動”指前文3.2部分計算結果)。可看出12.5 MPa預應力下的盤繞臂在小幅振動時頻率最大,同時在大幅非線性振動下頻率也最大。由此得出整體剛度越大,抵抗大幅度振動并維持剛度不降低的能力越強。

圖12 穩定圖Fig.12 Stability diagram

表4 結構自振頻率Tab.4 Natural vibration frequency Hz

振幅衰減區的自振頻率為22.73 Hz,在穩定圖14中,22.67 Hz附近穩定點較集中,以此為衰減區模態,振型為圖14。該階模態包含加勁索交替張弛的非線性現象,振型比較復雜同法可確定振幅穩定區結構頻率為24.97 Hz,振型為結構的一階彎曲變形。

圖13 22.73 Hz附近的穩定圖Fig.13 Stability diagram near 22.73 Hz

圖14 22.67 Hz對應的振型Fig.14 Vibration mode of 22.67 Hz

3.3.3 結果討論

(1) 能量守恒驗證

顯式動力學積分采用中心差分來近似積分,易出現累積誤差。通過驗證盤繞臂振動中的能量守恒來評價計算是否正確[21]。結構大幅非線性振動中的能量關系應滿足式(17)

ETOTAL=EK+EI

(17)

式中:EK為動能;EI為內能;ES為應變能。

圖15 振動過程中能量守恒驗證Fig.15 Verification of energy conservation in vibration

能量守恒定律要求振動中總能量等于常數。由于體系處于彈性范圍,還應有EI=ES。圖15顯示在振動中總能量最大值57.39 J,最小值56.97 J,平均值57.26 J,最大偏差0.51%,基本守恒;內能和應變能也基本相等。滿足能量守恒定律;振動穩定后,盤繞臂總能量在動能和應變能間相互轉換,保持動態平衡。

(2) 時程信號處理法與特征值迭代法對比

采用子空間迭代法進行特征值迭代考察的是結構剛度矩陣的固有屬性,當結構尺寸、材料參數、加勁索的預應力確定后,特征值迭代所求解出的結果是唯一確定、靜態不變的,表征結構小幅振動特性。

時程信號處理不僅可以考察不同振幅、不同預應力下加勁索發生不同程度松弛后的振動參數,也可考察同一次大幅振動中,由于振幅衰減使加勁索松弛節數越來越少,剛度特性發生改變的動態過程。

當然,采用時程信號處理同樣可計算小幅振動參數,取預應力10 MPa、12.5 MPa、15 MPa,采用時程信號處理法求解得到自振頻率分別為28.36 Hz、28.57 Hz、28.40 Hz,與特征值迭代的結果(表4第4列)相對誤差分別為1.97%、1.98%、1.70%,兩種法計算結果的統一性再次驗證了計算的正確性。

(3) 加勁索作用機理

文中計算均未考慮阻尼,圖11顯示大幅振動有振幅衰減階段,而圖16顯示小幅振動無振幅衰減段。為揭示加勁索的 “穩定索”功能,令δ等于5 mm,計算無索盤繞臂自振響應并進行對比(圖17)。結果顯示隨時間推移,振幅未出現衰減。這表明只有存在加勁索,才會發生振幅衰減現象。

圖16 盤繞式伸展臂小幅振動位移時程Fig.16 Displacement response of small amplitude vibration

圖17 無加勁索伸展臂振動位移時程Fig.17 Displacement response of non-cable mast

振動中加勁索受拉張緊時,縱桿對其做正功,縱桿變形能轉變成加勁索變形能,加勁索松弛后,加勁索的變形能又轉變成加勁索局部自振的動能。從圖15的能量轉化可看出,如此往復5個周期,體系能量最終轉變為整體結構中縱橫桿和勁索的應變能、動能,以及加勁索局部自振動的動能,并處于動態平衡。

(4) 試驗驗證

制作盤繞臂模型進行模態測試來驗證文章的計算方法見圖18。制作的盤繞臂參數見表1,由于加工條件限制,減小加勁索直徑為0.3 mm。

圖18 盤繞臂模態試驗裝置Fig.18 Modal test facility of coilable mast

采用掃描激光測振儀測試盤繞臂的自振頻率及模態。在盤繞臂不同位置施加振動激勵,使其分別發生彎曲、扭轉、彎曲+扭轉三種形式的振動。在盤繞臂一個面的各層節點上共布置20各信號測點,對盤繞臂進行振動掃頻,掃頻范圍為0~25 Hz。得到不同振動形式下頻譜響應曲線峰值點與相應模態見表5。綜合三種振動模式的測試結果,得到盤繞臂的各階固有頻率及對應模態見表 6(表 5中陰影部分)。

試驗中施加到加勁索中的預應力約為20 MPa,盤繞臂的最大振幅為0.4 mm,加勁索未發生松弛,故計算中采用小幅非線性振動分析方法。由于結點采用鑄鐵制作,相對盤繞臂整體結構有較大質量,故在計算中考慮節點質量,并以集中質量點形式添加到計算模型中,每個節點質量為28.7 g。計算結果表明,節點集中質量的存在使盤繞臂自振頻率有較大降低。

從表6可看出計算結果的模態基本與試驗結果吻合,但自振頻率偏大約3%。這是由于加工制作的結點存在間隙,導致實際結構剛度偏小;同時計算中未考慮結構阻尼也會導致計算結果偏大。總體來看計算值與實驗值較為接近,計算方法是可靠的。

表5 盤繞臂振動響應曲線峰值點及對應模態Tab.5 The peak points and modals of vibration response curve of coilable mast Hz

表6 自振頻率及模態試驗與計算結果對比Tab.6 Comparison of the experiments and simulation in natural frequencies and modal Hz

4 結 論

盤繞臂的振動特性不僅決定振動控制措施,還決定加勁索的預應力取值。文章通過幾何精確理論和屈曲分析理論的推導,得出靜力條件下加勁索預應力被動滿足的式;基于非線性有限元法,文章計算并總結了P-Δ效應對盤繞臂小幅振動的影響,總結出設計用“頻率-預應力-直徑”參考面,預應力在該取值范圍內可進行主動設計以達到對結構剛度、振動特性進行主動控制的目的;基于顯式動力學積分獲取結構振動時程信號,采用自編制的時程信號處理程序來獲取大幅振動的模態、頻率,在大幅振動控制中可起到參考作用,并揭示了大幅振動時加勁索有助于振幅迅速衰減的機理。文章提出的計算方法也可作為預應力張拉結構體系非線性振動分析的參考案例。

參 考 文 獻

[1] 丁鋒.典型桿狀構架式展開機構[J].上海航天,2006(1): 35-40.

DING Feng, Typical mast-liked deployment mechanism with truss structure[J].Aerospace Shanghai, 2006(1): 35-40.

[2] MCEACHEN M, TRAUTT T.Confirmation of new analytics for ultra-light lattice column strength using a 40-m flight article[C]//50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference.4-7 May 2009, Palm Springs, California.

[3] 劉濤, 韓涵, 冀賓,等.盤繞式伸展機構非線性屈曲模式的試驗研究與數值仿真[J].南京航空航天大學學報, 2015, 47(6): 897-903.

LIU Tao, HAN Han, JI Bin, et al.Experimental investigation and numerical simulation on nonlinearlinear buckling mode of coil-able mast[J].Journal of Nanjing University of Aeronauticus & Astronautics, 2015, 47(6): 897-903.

[4] 劉濤, 金玉龍, 呂榕新,等.FASTMast伸展機構屈曲模式的理論分析與試驗[J].宇航學報, 2015(4): 404-409.

LIU Tao, JIN Yulong, Lü Rongxin, et al.Theoretical analysis and experimental investigation on the buckling mode of FASTMast deployable truss structure[J].Journal of Astronautics, 2015(4): 404-409.

[5] 薛紜, 劉昭.受圓柱面約束彈性細桿的摩擦平衡分析[J].力學季刊, 2016(1): 139-148.

XUE Yun, LIU Zhao.Analysis of equilibrium with friction of an elastic rod constrained by a cylindrical surface[J].Chinese Quarterly of Mechanics, 2016(1): 139-148.

[6] 張金龍.盤繞式空間伸展臂力學行為研究[D].上海:上海交通大學, 2014.

[7] LUBINSKI A, ALTHOUSE W S.Helical buckling of tubing sealed in packers[J].Journal of Petroleum Technology, 2013.

[8] MCEACHEN M.Validation of SAILMAST technology and modeling by ground testing of a full-scale flight article[C]//48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 4-7 January 2010, Orlando, Florida.

[9] MCEACHEN M, TRAUTT T, MURPHY D.The ST8 SAILMAST Validation Experiment[C]//46th AIAA/ASME/ASCE/AHS/ ASC Structures, Structural Dynamics&Materials Conference, 18-21 April 2005, Austin, Texas.

[10] 戈冬明, 陳務軍, 付功義, 等.盤繞式空間可展折疊無鉸伸展臂的屈曲分析理論研究[J].計算力學學報, 2007, 24(5): 615-619.

GE Dongming, CHEN Wujun, FU Gongyi, et al.Buckling analysis on hingeless coilable mast[J].Chinese Journal of Computational Mechanics, 2007, 24(5): 615-619.

[11] 戈冬明, 陳務軍, 董石麟.盤繞式空間可展伸展臂展開狀態結構動力研究[C].第四屆海峽兩岸結構與巖土工程學術研討會, 杭州: 2007.

[12] 戈冬明.盤繞式空間可展伸展臂折疊屈曲機理與結構動力分析[D].上海:上海交通大學,2007.

[13] 韓建斌, 黃海, 馬海波.考慮斜拉索松弛的盤繞式伸展臂振動模型[J].北京航空航天大學學報, 2014, 40(7): 970-977.

HAN Jianbin, HUANG Hai, MA Haibo.Vibration model of coilable mast considering slack diagonals[J].Journal of Beijing University of Aeronautics and Astronautics,2014, 40(7): 970-977.

[14] 陳政清,楊孟剛.梁桿索結構幾何非線性有限元-理論、數值實現與應用[M].北京: 人民交通出版社, 2013.

[15] 杜家政, 王麗, 羅明智.基于屈曲分析結果的單元荷載剛度矩陣計算[C].北京力學會第18屆學術年會.北京: 2012年2月.

[16] 張年文.荷載剛度矩陣的推導和應用[J].茂名學院學報, 2004, 14(3): 45-48.

ZHANG Nianwen.Derivation and application of geometric stiffness matrix[J].Journal of Maoming College, 2004, 14(3): 45-48.

[17] 王勖成.有限單元法[M].北京: 清華大學出版社, 2003.

[18] 劉晶波,杜修力.結構動力學[M].北京: 機械工業出版社, 2005.

[19] 祁泉泉.基于振動信號的結構參數識別系統方法研究[D].北京:清華大學, 2011.

[20] 陳務軍,張淑杰.空間可展開結構體系與分析導論[M].北京: 中國宇航出版社, 2006.

[21] 莊茁,由小川,廖劍暉,等.基于ABAQUS的有限元分析和應用[M].北京: 清華大學出版社, 2009.

[22] BAB S, KHADEM S E, ABBASI A, et al.Dynamic stability and nonlinear vibration analysis of a rotor system with flexible/rigid blades[J].Mechanism and Machine Theory, 2016,105: 633-653.

[23] HOU L, CHEN Y, CAO Q, et al.Nonlinear vibration analysis of a cracked rotor-ball bearing system during flight maneuvers[J].Mechanism and Machine Theory, 2016,105: 515-528.

[24] BAB S, KHADEM S E, ABBASI A, et al.Dynamic stability and nonlinear vibration analysis of a rotor system with flexible/rigid blades[J].Mechanism and Machine Theory, 2016, 105: 633-653.

[25] SOLTANI P, SABERIAN J, BAHRAMIAN R.Nonlinear vibration analysis of single-walled carbon nanotube with shell model based on the nonlocal elasticity theory[J].Journal of Computational & Nonlinear Dynamics, 2015, 11(1): 011002.

猜你喜歡
模態振動結構
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
國內多模態教學研究回顧與展望
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 香蕉蕉亚亚洲aav综合| 91精品伊人久久大香线蕉| 高潮毛片免费观看| 高h视频在线| 国产丝袜丝视频在线观看| 亚洲91精品视频| 一本综合久久| 40岁成熟女人牲交片免费| 国产91麻豆免费观看| 亚洲制服中文字幕一区二区| 成人在线观看不卡| a亚洲天堂| 久草视频一区| 国产99在线观看| 五月婷婷中文字幕| 亚洲综合经典在线一区二区| 中文字幕在线免费看| 国产成人三级| 亚洲人成亚洲精品| 欧美日韩国产精品va| 黄色网址手机国内免费在线观看| 美女毛片在线| 最新精品国偷自产在线| 亚洲成人播放| 国产在线观看一区精品| 欧美一区精品| 亚洲综合狠狠| 国产精品9| а∨天堂一区中文字幕| 亚洲国产欧美自拍| 成人在线亚洲| 97超级碰碰碰碰精品| 久久久久青草大香线综合精品| 久久精品视频一| 亚洲综合第一页| 天天色天天综合| 114级毛片免费观看| 777午夜精品电影免费看| 欧美日韩国产成人高清视频| 亚洲日韩高清在线亚洲专区| 国产综合精品一区二区| 國產尤物AV尤物在線觀看| 一级看片免费视频| 国产人成网线在线播放va| 日本91视频| 国产精品亚洲精品爽爽| 综合色88| av一区二区三区高清久久| 精品自拍视频在线观看| 久久免费视频播放| 欧美国产综合色视频| 国产精品自拍露脸视频| 在线国产91| 国产H片无码不卡在线视频 | 国产喷水视频| 色视频久久| 欧美激情视频在线观看一区| AV在线天堂进入| 日本午夜三级| 国产97视频在线观看| 欧美日韩亚洲国产| 精品欧美一区二区三区久久久| 夜夜操天天摸| 热re99久久精品国99热| 精品国产91爱| 国产精品永久在线| 国产成人精品一区二区三区| 国国产a国产片免费麻豆| 色婷婷狠狠干| 亚亚洲乱码一二三四区| 午夜无码一区二区三区| 国产日本欧美亚洲精品视| 成人小视频网| 国产无码在线调教| 五月婷婷伊人网| 尤物精品视频一区二区三区| 国语少妇高潮| 国产91小视频在线观看| 男女性色大片免费网站| 91成人在线免费视频| 一级毛片在线直接观看| 国产a网站|