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

基于模態綜合法與動態凝聚的滾珠絲杠進給系統狀態空間建模與實驗驗證

2017-07-18 11:49:22張為民仇煒諫
振動與沖擊 2017年13期
關鍵詞:模態模型系統

楊 勇, 張為民, 仇煒諫, 陳 盛

(1. 江蘇科技大學 機械工程學院, 江蘇 鎮江 212003; 2. 同濟大學 機械與能源工程學院, 上海 201804;3. 南京郵電大學 自動化學院, 南京 210023)

基于模態綜合法與動態凝聚的滾珠絲杠進給系統狀態空間建模與實驗驗證

楊 勇1, 張為民2, 仇煒諫1, 陳 盛3

(1. 江蘇科技大學 機械工程學院, 江蘇 鎮江 212003; 2. 同濟大學 機械與能源工程學院, 上海 201804;3. 南京郵電大學 自動化學院, 南京 210023)

以往研究中為便于機械傳動環節與控制系統的集成建模,常將滾珠絲杠進給系統簡化為簡單機械模型,其在動態性能精確求解等方面存在不足。針對這一問題,充分考慮滾珠絲杠進給系統的結構柔性以及絲杠螺母副的耦合特性,構建滾珠絲杠進給系統多柔體動力學模型;在此基礎上,綜合模態綜合法和動態凝聚法,對進給系統有效簡化動力學信息準確提取,加快模型降階收斂和凝聚自由度響應求解速度,推導得到滾珠絲杠進給系統的狀態空間縮聚模型。基于內置傳感器信號的動態性能試驗測試,對模型的正確性進行驗證,其一階固有頻率計算誤差為5.54%,實驗和仿真計算得到的輸入扭矩與輸出轉速間的頻響曲線在一階頻率處極點坐標分別為32.83 dB和35.16 dB。該研究為實現滾珠絲杠進給系統與伺服控制系統的高效高精度集成建模提供了重要的理論基礎和方法指導。

滾珠絲杠進給系統; 狀態空間; 模態綜合; 動態凝聚

機床進給驅動系統中,其機械傳遞環節對整體進給系統動態特性有著關非常重要的影響,采用不同方法對滾珠絲杠進給系統建模并對其動態性能進行研究,一直受到國內外眾多學者的關注。

目前滾珠絲杠進給系統的建模方法大致可以分為五類。① 滾珠絲杠進給系統的解析法建模[1-4]:采用解析法運用結構動力學相關原理構建滾珠絲杠進給系統動力學模型,對其動力學性能預測評價。其特點是:多用作理論分析或定性分析,在精確動力學性能預測評價方面存在不足。② 滾珠絲杠進給系統的傳遞矩陣法建模[5]:采用傳遞矩陣法,以矩陣的形式推導系統數學模型,求解該矩陣能夠快速方便地對系統動力學性能進行研究分析。其特點是:不受系統部件組成形式的限制,多用作定性分析,在精確動力學性能預測評價方面存在不足。③ 滾珠絲杠進給系統的集中質量法建模[6-7]:采用集中質量法構建滾珠絲杠進給系統動力學模型,方便快速性的對系統穩定性影響因子及變化規律做出判定,利于系統定性分析與振動頻率估計。有效減少模型的自由度,降低傳遞函數階數,便于與伺服控制系統實現集成建模,但過于簡化的研究方法與建模手段,抑制了部分對系統性能造成較大影響的重要因素如系統結構柔性、結合部剛度等,同時此方法忽略了滾珠絲杠的細長桿特征。④ 滾珠絲杠進給系統的混合建模法[8-11]:采用混合建模法對滾珠絲杠進給系統進行動力學建模及動力學性能分析,重要零部件采用分布參數建模,其余部件采用集中質量建模。該建模方法相對于集中質量法有所改進,但其仍未能充分考慮系統結構柔性對其動力學性能的影響。⑤ 滾珠絲杠進給系統的有限元法建模[12]:采用有限元技術對滾珠絲杠進給系統進行動力學建模及動力學性能研究,可以方便準確的描述機床靜態特征之外,也可以實現動態時變特性問題的有效求解。但系統整體動力學模型較為復雜,自由度數目眾多,模型比較龐大,不利于與伺服控制集成。

綜上,為便于機械傳動環節與控制系統的集成,常將滾珠絲杠進給系統簡化為一個簡單的機械模型,但其在結構動力學性能準確研究預測、動力學響應精確求解等方面存在不足;而滾珠絲杠進給系統的有限元方法建模等,雖然可以實現其動態特性與動力學響應的精確求解,但是模型復雜、自由度龐大,無法實現與伺服控制系統的集成建模。

針對這一問題,本文充分考慮滾珠絲杠進給系統結構柔性以及絲杠螺母副的耦合特性等,首先構建滾珠絲杠進給系統多柔體動力學模型;在此基礎上,綜合模態綜合法和動態凝聚法,對進給系統有效簡化動力學信息簡化提取,建立基于多柔體模型的滾珠絲杠進給系統狀態空間模型,在對進給系統動力學響應精確求解的同時,便于與伺服控制系統集成建模。

1 滾珠絲杠進給系統多柔體動力學建模

在滾珠絲杠進給系統多柔體動力學建模過程中,首先對其主要功能結構部件進行劃分和分析,確定組成該進給系統的床身、立柱、電機軸、聯軸器、絲杠等各部分的結構簡化、單元體選擇等,在此基礎上,確定對應于各個部分的剛度、材料、密度、彈性模量等參數。為了保證分析計算精度,其立柱、床身、滑枕等部件采用修正的三維實體單元建模,其一方面有效克服單元剪切自鎖問題,另一方面也確保了單元的求解精度。對于電機軸、絲桿軸等采用考慮彎曲剪切綜合變形的粱單元進行建模。導軌滑塊采用具有橫向剛度和縱向剛度的彈簧單元進行建模。橫絲杠支撐軸承采用具有徑向剛度和軸向剛度的彈簧單元進行建模。聯軸器等可以采用扭轉彈簧進行建模。實際動力學建模過程中,聯軸器扭轉剛度為2.7×104N·m/rad,軸承的軸向剛度為4.6×109N/m,徑向剛度為2.5×108N/m,導軌滑塊的豎直剛度和切向剛度分別為1.88×109N/m、1.25×109N/m。立柱、滑枕、床身等實體單元材料密度7 300 kg/m3,彈性模量為140 GPa,泊松比為0.26,絲杠軸等的材料密度為7 800 kg/m3,彈性模量為206 GPa,泊松比為0.3。

考慮絲杠螺母副內耦合特性,依據文獻[13] 分析及推導中滾珠絲杠副其他剛度與軸向剛度的物理映射關系,建立絲杠螺母副動力學模型,滾珠絲杠副型號為Steinmeyer 3426/30.60[14],其基本參數為:絲桿的公稱直徑60 mm、滾珠直徑9 mm、接觸壓力角45°、滾珠工作圈數8圈、絲杠有效長度105 mm。由其技術參數手冊可知,滾珠絲杠副的軸向剛度為k1=1.65×109N·m-1,據前述其他剛度與軸向剛度的物理映射關系可得徑向、扭轉、彎曲及耦合剛度分別為:k2=k3=8.35×108N·m-1,k4=0.94×104N·m·rad-1,k5=k6=1.73×106N·m·rad-1,kc=3.94×106N·rad-1。最后得到的滾珠絲杠進給系統多柔體結構動力學模型如圖1所示。

圖1 滾珠絲杠進給系統多柔體結構動力學模型

2 基于狀態空間法的滾珠絲杠進給系統建模

在前述滾珠絲杠多柔體動力學建模的基礎上,將滾珠絲杠進給系統的結構動力學模型抽象為結構裝配體模型,可以得到無約束狀態下滾珠絲杠進給系統的結構動力學模型并寫為

(1)

滾珠絲杠進給系統的子結構A、子結構B及聯接結構J都應該滿足結構動力學普遍方程,即

(2)

式中,i=A、B或J。

(3)

(4)

(5)

(6)

令y=xi1=xi,將其寫成矩陣形式可得

(7)

聯立式(6)和式(7),可以得到滾珠絲杠進給系統子結構A、子結構B及聯接結構J的狀態空間方程的基本形式

(8)

圖2 滾珠絲杠進給系統的子結構(聯接結構)信號傳遞框圖

根據式(1)與式(8),可以得到滾珠絲杠進給系統的狀態空間模型為

(9)

實際上,方程式(9)中的空間坐標含有非獨立自由度,因此,需要建立其廣義自由度坐標下的滾珠絲杠進給系統的狀態空間模型。將其子結構與聯接結構的獨立自由度與非獨立自由度劃分,并通過廣義坐標變換矩陣進行變換,可以得到

(10)

(11)

式中:MΓ、CΓ、KΓ、QΓ分別為滾珠絲杠進給系統廣義坐標下的質量、阻尼、剛度、載荷矩陣,且

(12)

據式(6)、式(7)、式(11)可得廣義坐標獨立自由度下的滾珠絲杠進給系統狀態空間方程

(13)

基于廣義坐標獨立自由度的滾珠絲杠進給系統狀態空間建模示意圖見圖3。

圖3 基于廣義坐標獨立自由度的滾珠絲杠進給系統狀態空間模型

Fig.3 Ball screw feed system state space model based on generalized coordinate independent freedom

3 基于模態綜合法與動態凝聚的滾珠絲杠進給系統狀態空間建模

基于上述廣義坐標獨立自由度下的滾珠絲杠進給系統狀態空間方程,可以方便的實現與伺服控制系統的集成,但是由于滾珠絲杠進給系統有限元模型的復雜性,勢必會造成滾珠絲杠進給系統狀態空間方程模型自由度數量的龐大性,因此必須對上述模型進行縮減,實現滾珠絲杠進給系統狀態空間信息的準確提取與可靠簡化,從而確保滾珠絲杠進給系統多柔體動力學模型與伺服控制模型的高效率高精度集成。

在前述基礎上,采用動態凝聚技術,在保證計算精度的同時,加快模型降階收斂和凝聚自由度響應求解速度,推導得到縮聚后的滾珠絲杠進給系統的狀態空間模型,以便于與伺服控制高效高精度集成。

滾珠絲杠進給系統的廣義特征方程

KΦ-MΦΩ=0

(14)

式中,Φ、Ω分別為前n階固有振型矩陣和固有頻率矩陣。

物理坐標和廣義模態坐標之間的廣義坐標變換,如式(15)所示

xn×1=Φn×npn×1

(15)

式中,p為模態坐標。由式(15)可以看出,系統物理坐標運動均可由前n個模態坐標主振動分量疊加而成,由前述分析可知,滾珠絲杠進給系統節點空間自由度較多,因此其振型矩陣Φn×n亦較大,從而造成上述坐標變換過程的計算較繁瑣,且耗時較長。

基于模態綜合的基本原理,系統物理坐標下的運動向量xn×n可以近似用前m階(m≤n)以模態主坐標p為權重系數的振型常量Φ疊加表達,表達成矩陣形式為

(16)

將式(16)代入滾珠絲杠進給系統動力學方程式(11),可得

(17)

通過模態阻尼技術以模態阻尼系數[15-17]引入系統阻尼因素影響,得到系統阻尼矩陣CΓ,如式(18)所示

CΓ=2CL(μTKΓμ)1/2=2CLΛ

(18)

式中:CL=diag[cL,j],cL,j為模態阻尼比;Λ為模態剛度矩陣其可由前述多柔體建模及其結構動態特性與模態分析求解得到。模態阻尼系數可以通過實驗模態測試得到,有關系統模態阻尼比的測試計算方法研究已經較多[18-19],且已經比較成熟,而且目前這些測試計算方法已經被成功集成于常用的振動模態測試儀器(比較知名的振動測試儀品牌有LMS(LMS International)、B&K(Brüel & Kjaer)等)中,在實際的模態測試實驗及工程實踐中得到廣泛應用[20]。綜合以往相關方面研究,借助對機床樣機的實驗測試,最終確定模態阻尼比取值為0.05。

利用式(19)可以實現以廣義模態坐標形式表達的物理坐標下的邊界約束及平衡條件

Hpm×1≈0

(19)

式中,H為廣義模態坐標下的邊界約束及平衡系數矩陣。

按照獨立度自由度s和非獨立自由度r對上式進行分塊,并最終可以得到

(20)

將式(20)代入滾珠絲杠進給系統動力學方程式(17),則可寫出類似于式(13)的縮聚后的滾珠絲杠進給系統狀態空間方程,系統動力學響應求解可以通過式(16)、式(20)的坐標變換進行求解,如圖4所示。

圖4 基于滾珠絲杠進給系統狀態空間模型的動力學響應求解示意圖

4 實驗驗證

本文基于內置傳感器信號的動態性能試驗測試方法[21-23]對滾珠絲杠進給系統的狀態空間模型進行驗證。基于內置傳感器的進給系統動態性能測試采用集成安裝于機床進給系統內的編碼器、伺服電動機電流(轉矩)傳感器等本體信息,由于其為機床控制系統的反饋信號,是進給系統動力學響應特性的直接反映,包含了進給系統運行過程中的輸入輸出信號,是進給系統動態性能測試的理想信號,且由于無需采用額外的外置傳感器,可以有效降低測試成本。滾珠絲杠進給系統動態性能測試原理如圖5所示。本文所研究機床采用西門子數控系統,西門子公司已經成功開發并集成了基于內置傳感器的伺服控制響應測試軟件與模塊,并在實際中得到廣泛應用[24-25],滾珠絲杠進給系統動態性能測試系統則基于西門子840D系統面板控制單元PCU50以及操作界面軟件HMI Advanced實現。

圖5 基于內置傳感器的滾珠絲杠進給驅動系統動態性能測試原理

基于滾珠絲杠進給系統狀態空間模型的實際轉速與輸入扭矩之間的頻響曲線見圖6,其實驗測試與仿真計算比較見表1。圖表綜合分析可以看出,幅值對應極點位置與前述進給系統固有頻率數值基本吻合,由于采用模態綜合法等對滾珠絲杠進給系統模型進行縮聚,從而會導致仿真計算與實驗測試間存在誤差,但前兩階仿真計算固有頻率為47.96 Hz、90.74 Hz,與實測固有頻率結果45.44 Hz、85.94 Hz誤差較小,一階固有頻率誤差為5.54%。在機械固有頻率處實測結果與仿真計算結果相差均約在3 dB左右,仿真計算和實驗測試結果吻合度較好,從而驗證了滾珠絲杠進給系統狀態空間建模的正確性。

(a)

(b)

表1 實驗測試與仿真計算比較

5 結 論

(1) 本文充分考慮滾珠絲杠進給系統結構柔性以及進給系統內部耦合特性,首先構建滾珠絲杠進給系統多柔體動力學模型。在此基礎上,采用模態綜合法和動態凝聚法,加快模型降階收斂和凝聚自由度響應求解速度,對進給系統有效簡化動力學信息準確提取,建立滾珠絲杠進給系統的狀態空間縮聚模型,在確保動力學響應求解精度的同時,便于與伺服控制系統集成。

(2) 基于內置傳感器信號的動態性能試驗測試方法對滾珠絲杠進給系統的狀態空間模型進行了驗證,其一階固有頻率計算誤差為5.54%,實驗和仿真得到的輸入扭矩與輸出轉速間的頻響曲線在一階頻率處極點坐標分別為32.83 dB和35.16 dB,仿真計算和實驗測試結果吻合度較好。

(3) 本文研究為實現滾珠絲杠進給系統與伺服控制系統的高效高精度集成建模提供了重要的理論基礎和方法指導。

[1] 張會端,譚慶昌,李慶華. 機床傳動絲杠的動力分析[J]. 農業機械學報, 2009,40(9): 220-226.

ZHANG Huiduan, TAN Qingchang, LI Qinghua. Dynamic analysis of the machine drive screw[J].Transactions of the Chinese Society for Agricultural Machinery, 2009,40(9): 220-226.

[2] GU U C, CHENG C C. Vibration analysis of a high-speed spindle under the action of a moving mass[J]. Journal of Sound and Vibration, 2004, 278(4/5): 1131-1146.

[3] 張會端,孫俊嶺. 彈性支承條件下傳動絲杠的橫向振動分析[J]. 長春大學學報, 2011,21(2): 16-20.

ZHANG Huiduan, SUN Junling. The lateral vibration analysis of the drive screw under the elastic supports[J].Journal of Changchun University, 2011,21(2): 16-20.

[4] 王永強,張承瑞. 滾珠絲杠進給系統仿真建模[J]. 振動與沖擊, 2013,32(3): 46-49.

WANG Yongqiang, ZHANG Chengrui. Simulation modeling of a ball screw feed drive system[J]. Journal of Vibration and Shock, 2013,32(3): 46-49.

[5] 吳文鏡,劉強. 機床動力學建模的拓展傳遞矩陣法[J]. 機械工程學報, 2010,46(21): 69-75.

WU Wenjing, LIU Qiang.Transfer matrix method for dynamic modeling of machine tools[J]. Journal of Mechanical Engineering, 2010,46(21): 69-75.

[6] CHEN J S, HUANG Y K, CHENG C C. Mechanical model and contouring analysis of high-speed ball-screw drive systems with compliance effect[J]. International Journal of Advanced Manufacturing Technology, 2004, 24(3/4): 241-250.

[7] POP L C. Particularities of modeling ball screw based NC axes as finite degrees of freedom dynamic systems[J]. Buletinul Institution polytechnic Din Lasi, 2005, 5: 1-6.

[8] ALEYAASIN M, EBRAHIMI M, WHALLEY R. Vibration analysis of distributed-lumped rotor systems[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 189(2): 545-558.

[9] WHALLEY R, ABDUL-AMEER A A, EBRAHIMI K M. The axes response and resonance identification for a machine tool[J]. Mechanism and Machine Theory, 2011, 46(8): 1171-1192.

[10] WHALLEY R, EBRAHIMI M, ABDUL-AMEER A A. Hybrid modelling of machine tool axis drives[J]. International Journal of Machine Tools & Manufacture, 2005, 45(14): 1560-1576.

[11] 董亮,湯文成,劉立. 滾珠絲杠進給系統混合建模及其振動時變性分析[J]. 振動與沖擊, 2013,32(20): 196-202.

DONG Liang, TANG Wencheng, LIU Li. Hybird modeling and time-varying analysis of vibration for a ball screw drive[J].Journal of Vibration and Shock, 2013,32(20): 196-202.

[12] HOLROYD G, PISLARU C, FORD D G. Determination of stiffness and damping sensitivity for computer numerically controlled machine tool drives[J]. Proceedings of the Institution of Mechanical Engineers Part C-Journal of Mechanical Engineering Science, 2003, 217(10): 1165-1177.

[13] 楊勇,張為民,楊濤. 基于Kriging元模型的機床進給驅動系統動態特性優化[J]. 農業機械學報, 2013,44(5): 288-293.

YANG Yong, ZHANG Weimin, YANG Tao. Dynamic characteristic optimization of feed system based on Kriging metamodel[J].Transactions of the Chinese Society for Agricultural Machinery, 2013,44(5):288-293.

[14] Steinmeyer. Single nut with 4-point-contact[EB/OL]. http://www.steinmeyer.com/en/products/drive-technology.

[15] THOMSON W. Theory of vibration with applications[M]. Beijing: Tsinghua University Press, 2005.

[16] CHOPRA A K. Dynamics of structures: theory and applications to earthquake engineering[M]. Beijing: Higher Education Press, 2007.

[17] CLOUGH R W, PENZIEN J, WANG G Y. Dynamics of structures[M]. Beijing: Higher Education Press, 2006.

[18] NIEHUES K, SCHWARZ S, ZAEH M F. Reliable material damping ratio determination in machine tool structures[J]. Production Engineering, 2012, 6(4/5): 475-484.

[19] 管迪華. 模態分析技術[M]. 北京: 清華大學出版社, 1996.

[20] HUNG J, LAI Y, LIN C, et al. Modeling the machining stability of a vertical milling machine under the influence of the preloaded linear guide[J]. International Journal of Machine Tools and Manufacture, 2011, 51(9): 731-739.

[21] ANDOLFATTO L, LAVERNHE S, MAYER J. Evaluation of servo, geometric and dynamic error sources on five-axis high-speed machine tool[J]. International Journal of Machine Tools and Manufacture, 2011, 51(10): 787-796.

[22] 周玉清,梅雪松. 基于無傳感器信息的數控機床伺服進給系統性能評估研究[J]. 機械工程學報, 2012(20): 32.

ZHOU Yuqing, MEI Xuesong. Research on performance evaluation of CNC machine tools servo feed system based on non-sensor information[J]. Journal of Mechanical Engineering, 2012(20): 32.

[23] TOUNSI N, BAILEY T, ELBESTAWI M A. Identification of acceleration deceleration profile of feed drive systems in CNC machines[J]. International Journal of Machine Tools and Manufacture, 2003, 43(5): 441-451.

[24] 陳先鋒. SIEMENS數控技術應用工程師[M]. 北京:人民郵電出版社, 2010.

[25] 胡國清,張旭宇. 西門子SINUMEIRK840D sl/840Di sl數控系統應用工程師手冊[M]. 北京:國防工業出版社, 2013.

State space modeling and test verification for a ball screw feed system based on modal synthesis and dynamic condensation

YANG Yong1, ZHANG Weimin2, QIU Weijian1, CHEN Sheng3

(1. College of Mechanical Engineering, Jiangsu University of Science and Technology, Zhenjiang 212003, China;2. College of Mechanical Engineering, Tongji University, Shanghai 201804, China;3. College of Automation, Nanjing University of Posts and Telecommunication, Nanjing 210023, China)

In order to integrated modeling mechanical transmission segment and servo control system, a ball screw feed system is simplified as a simple mechanical model in previous studies. However, this model is not good for precisely solving its dynamic performance. Here, aiming at this problem, fully considering the structural flexibility and coupling characteristics of screw-nut pair, a ball screw feed system was modeled as a multi-flexible body dynamic model. Then, the dynamic information of the effectively simplified feed system was extracted correctely with the modal synthesis method and dynamic condensation one to speed up the model order reduction and convergence, and condensed DOFs’ responses solving. The state-space reduction model of the ball screw feed system was derived finally. Based on the built-in sensors, the system’s dynamic performance tests were conducted to verify the correctness of this model. It was shown that the calculation error of the model’s the first order natural frequency is 5.54%, the test and simulation results of the system’s frequency responses are 32.83 dB and 35.16 dB at the first-order natural frequency point, respectively. The study results provided an important theoretical basis and a technical guide for integrated modeling a ball screw feed system and its servo control system.

ball screw feed system; state-space; modal synthesis; dynamic condensation

高檔數控機床與基礎制造裝備科技重大專項(2011ZX04016-021);國家科技支撐計劃課題(2015BAF06B02-4);江蘇省高校自然科學研究面上資助項目(16KJB460017)

2016-03-21 修改稿收到日期:2016-05-12

楊勇 男,博士,講師,1985年12月生

張為民 男,博士,教授,1965年1月生

TH132.1

A

10.13465/j.cnki.jvs.2017.13.008

猜你喜歡
模態模型系統
一半模型
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
3D打印中的模型分割與打包
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
主站蜘蛛池模板: 亚洲免费福利视频| 99re精彩视频| 国产精品免费福利久久播放 | 无码又爽又刺激的高潮视频| 国产成人喷潮在线观看| 亚洲精品麻豆| 欧美日韩国产综合视频在线观看 | 欧美黄网站免费观看| 国产乱子伦精品视频| 亚洲成人播放| 免费看的一级毛片| 国产成人毛片| 强奷白丝美女在线观看| 国产毛片不卡| 又爽又大又光又色的午夜视频| 日韩精品成人在线| 久久综合一个色综合网| 国产在线精品香蕉麻豆| 麻豆国产精品一二三在线观看| 亚洲精品无码抽插日韩| 国产国产人在线成免费视频狼人色| 成人精品午夜福利在线播放| 波多野结衣无码中文字幕在线观看一区二区 | 亚洲天堂视频网| 国产成人精品综合| 久久综合五月婷婷| 国产在线98福利播放视频免费| 高清精品美女在线播放| 亚洲VA中文字幕| 欧美伊人色综合久久天天| 亚洲精品爱草草视频在线| 91精品国产91久无码网站| 久久久久人妻一区精品色奶水| 亚洲人成影视在线观看| 2020精品极品国产色在线观看| 久久黄色小视频| 国产SUV精品一区二区6| 色婷婷在线播放| 香蕉色综合| 日韩精品无码免费一区二区三区 | 国产乱肥老妇精品视频| 国产日韩欧美一区二区三区在线| 国产精品分类视频分类一区| 在线观看欧美国产| 国产美女91呻吟求| 高潮毛片免费观看| 真实国产乱子伦视频| 不卡视频国产| 国产免费a级片| 一级毛片免费高清视频| 狠狠色噜噜狠狠狠狠色综合久| 动漫精品中文字幕无码| 高清欧美性猛交XXXX黑人猛交 | 另类专区亚洲| 欧美日韩在线成人| 国产欧美视频综合二区| 亚洲精品动漫| 91久久偷偷做嫩草影院| 国产噜噜噜视频在线观看| 老司国产精品视频91| 日本精品中文字幕在线不卡| 国产福利小视频高清在线观看| 国产精品成人一区二区| 国产在线八区| 激情综合激情| 国产亚洲欧美日韩在线观看一区二区| 97se亚洲综合在线韩国专区福利| 久久黄色一级片| 国产精品成人观看视频国产| 91精品久久久久久无码人妻| 欧美国产在线一区| 美女免费黄网站| 午夜无码一区二区三区| 国产精品尤物在线| 在线观看无码av免费不卡网站| 91青青在线视频| 欧美黄网在线| 激情视频综合网| 久久久久青草线综合超碰| 国产成人毛片| 91久久国产成人免费观看| 免费人成在线观看成人片|