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

基于神經網絡建模的機床滑動結合面動態特性參數識別

2018-04-24 08:07:25朱堅民周亞南何丹丹鄭洲洋
振動與沖擊 2018年7期
關鍵詞:模態有限元優化

朱堅民, 周亞南, 何丹丹, 鄭洲洋

(上海理工大學 機械工程學院,上海 200093)

機床中存在各類結合面,大量的結合面對機床結構的整體特性有著重要影響[1]。結合面的存在會降低機床結構的局部剛度,直接影響其機械性能。結合面的接觸剛度和接觸阻尼是機床剛度和阻尼的重要組成部分,眾多研究表明:機床靜剛度中30%~50%決定于結合部的剛度特性[2];機床出現的振動問題有60%以上來源于結合部;機床的阻尼90%以上來源于結合部的接觸阻尼[3]。目前,國內外對固定結合面[4]和滾動結合面[5]動態特性的理論和實驗研究較多,而對滑動結合面的研究相對較少。滑動結合面作為整機系統中重要的結合面之一,其結合面動態特性參數的準確識別對提高機床的設計水平和加工質量具有重要意義。

針對機床滑動結合面的建模及其動態特性參數的識別問題,國內外學者主要從三個方面進行研究:理論解析法、實驗測試法、理論與試驗相結合的方法。在理論解析法研究方面,Edward等[6]利用有限元理論分析了滑動導軌結合部的靜態特性及其接觸剛度的產生機理。王禹林等[7]采用吉村允孝法確定機床滑動導軌結合面的動態特性參數,并利用ANSYS建立其有限元模型。在試驗測試法研究方面,郭成龍等[8]通過模態實驗法識別了貼塑導軌滑動結合面的動態特性參數。伍良生等[9]對自行設計的滑動導軌實驗臺進行了掃頻振動測試,獲得了結合面的單位面積剛度和阻尼。在理論與實驗相結合的研究方面,Lee等[10]通過有限元模型與實驗測量機床部件位移得到的頻響函數,建立優化模型進而識別出導軌結合部的剛度矩陣,識別誤差<7%。王立華等[11]基于ANSYS軟件的優化設計模塊結合模態實驗對銑床關鍵結合面的特性參數進行了識別,識別誤差<8%。孫明楠等[12]將模態實驗數據與有限元分析模型相結合,通過擬牛頓BFGS算法識別了機床導軌結合部動態特性參數,識別誤差<9%。Deng等[13]在MATLAB與ANSYS相互集成環境下結合模態試驗建立優化任務,求得導軌結合面的動態特性參數,識別誤差<5%。理論解析與實驗測試相結合的方法充分融合了前兩種方法的優點,用實驗結果不斷修正理論模型,建模精度較高。但基于修正的理論模型進行參數優化識別時,由于每一步優化迭代都要進行費時的理論計算,導致參數識別效率較低。因此,有學者選用擬合精度和效率較高的擬合模型近似代替理論模型。汪中厚等[14]將響應面法與實驗模態測試相結合識別了滑動結合面的動態特性參數,識別誤差<5%。響應面法在很大程度上解決了每次優化迭代中需調用耗時的理論模型計算的問題,但基于響應面法等建模方法所建立的結合面擬合模型精度不高,影響了結合面參數的識別精度。

針對上述問題,以自行設計的滑動結合面實驗臺為研究對象,提出了利用神經網絡建立其工作臺-床身滑動導軌結合面的擬合模型,并結合布谷鳥優化算法對滑動結合面動態特性參數進行優化識別的方法,獲得了較高的參數識別精度。

1 基本原理

1.1 滑動結合面的建模方法

機床滑動結合面在動態力的作用下,表現出既有彈性又有阻尼,可采用一系列彈簧-阻尼單元組成的等效動力學模型近似表征工作臺和床身間滑動導軌的結合面。通過合理確定結合點的數目、每個結合點的自由度以及每個自由度的等效剛度和等效阻尼參數,可模擬滑動結合面不同條件和狀態下的動態特性。本文建立的機床滑動結合面模型,如圖1所示。圖1中,結合面上的彈簧-阻尼單元均勻分布,每個結合點上三個方向的等效剛度為kx、ky、kz,三個方向的等效阻尼為cx、cy、cz。利用有限元軟件建立機床滑動結合面模型時,可采用ANSYS提供的Matrix27單元,該單元特性類似彈簧-阻尼單元,通過定義其剛度和阻尼參數來模擬滑動結合面的接觸特性,即在兩個結合面相對應的結合點位置處生成節點,用Matrix27剛度和阻尼單元將兩個節點連接。依據Matrix27單元的特性和并聯彈簧-阻尼單元的等效原理,進而建立機床滑動結合面的剛度和阻尼參數模型。

1.2 神經網絡模型的建立

為了對圖1中滑動結合面的等效剛度參數和等效阻尼參數進行優化識別,本文采用2個神經網絡模型分別建立剛度參數(kx、ky、kz)與機床前四階理論固有頻率(f1、f2、f3、f4)、阻尼參數(cx、cy、cz)與機床前四階理論阻尼比(ζ1、ζ2、ζ3、ζ4)之間的模型。神經網絡模型的訓練數據來源于對包含圖1中滑動結合面的機床整機有限元模型進行理論模態分析的結果,神經網絡模型的建模與訓練原理,如圖2所示。首先,為了使神經網絡模型能得到有效充分的訓練,需要通過有效的試驗設計方法,以獲得全面反映機床滑動結合面動態特性的樣本點。本文采用優化拉丁超立方試驗設計方法確定了n組剛度參數樣本點(kxi、kyi、kzi,i=1,2,…,n)和n組阻尼參數樣本點(cxi、cyi、czi,i=1,2,…,n),這些樣本點在采樣空間中分布均勻。其次,建立包含圖1中滑動結合面剛度(阻尼)參數的機床整機有限元模型,根據每一組樣本點對機床整機有限元模型進行理論模態分析,獲得機床前四階理論固有頻率(f1i、f2i、f3i、f4i,i=1,2,…,n)和理論阻尼比(ζ1i、ζ2i、ζ3i、ζ4i,i=1,2,…,n)。實際計算時可由ISIGHT集成ANSYS進行相應的有限元分析和計算,最后依次得到n組剛度參數樣本點與機床前四階理論固有頻率、n組阻尼參數樣本點與機床前四階理論阻尼比之間的對應關系,并分別被用于對圖2中的剛度參數神經網絡模型、阻尼參數神經網絡模型進行訓練。

圖1 機床工作臺-床身滑動結合面的等效模型

Fig.1 Equivalent model of sliding joints of machine table-bed

1.3 動態特性參數的優化識別

剛度參數優化過程為:確定優化變量為剛度參數(kx,ky,kz),其上下限及變化范圍可通過理論分析和參數靈敏度分析綜合后確定,所構建的剛度參數優化目標函數和約束如式(1)所示

(1)

圖2 剛度(阻尼)參數神經網絡模型的建模與訓練

為了保證尋優的效率,本文采用布谷鳥搜索算法[15]對式(1)的優化問題進行求解,并由此確定圖1中剛度參數kx,ky,kz的最優解。阻尼參數(cx,cy,cz)的優化過程與此類似,所建立的優化目標函數和約束條件如式(2)所示

(2)

2 應用實例

本文以自行設計制造的機床滑動結合面實驗臺為研究對象,該機床的主要組成部件有床身、工作臺、砂輪箱、滑座、立柱等,機床照片如圖4所示。本文對機床的工作臺-床身滑動導軌結合面進行建模,并識別其動態特性參數。

2.1 理論模態分析

考慮圖4機床的實際約束情況,以及本文重點研究機床工作臺-床身滑動導軌結合面,在對機床進行建模時,其立柱與床身結合處使用剛性連接,床身底面建立全約束,砂輪箱-滑座及其他部件間結合面暫不考慮按固結處理。工作臺的V形滑動導軌有兩個大小相等的滑動結合面,長為1 000 mm,寬為38.77 mm;矩形滑動導軌有一個結合面,長為1 000 mm,寬為60 mm。對這三個滑動結合面按圖1所示分別建模時,將其劃分為7部分,在每兩部分的連接中心處分布一對彈簧阻尼單元,即每個結合面通過均布的6對彈簧阻尼單元連接上下兩個結構。本文所建立的機床整機有限元模型,如圖5所示。在進行模態分析時,采用ANSYS中的Block Lanczos法求取機床整機的前四階固有頻率,采用QR Damped法求取機床阻尼比。

2.2 建立神經網絡模型

2.2.1 訓練樣本與網絡初始化

為獲取神經網絡的訓練樣本數據,首先需確定剛度參數、阻尼參數的范圍。根據吉村允孝理論[16],估算出機床工作臺-床身滑動結合面的剛度為3.87 kN/mm,阻尼為2.58 N·s/mm。根據已有相關文獻的研究結果及對結合面上各彈簧-阻尼單元的剛度和阻尼參數進行靈敏度分析后,確定本文研究的滑動結合面剛度變化范圍為(10~104)N/mm,阻尼參數的變化范圍為(0.5~2.8)N·s/mm。兼顧參數優化迭代所需時間的限制,最終確定剛度參數的范圍為(0.5~6.0)kN/mm,阻尼參數的范圍為(1.0~2.8)N·s/mm。

利用ISIGHT的試驗設計模塊中優化拉丁超立方的方法確定120組設計變量樣本點(kxi、kyi、kzi,i=1,2,…,120)、(cxi、cyi、czi,i=1,2,…,120),然后通過ISIGHT集成ANSYS進行相應的有限元分析和計算,樣本點的計算采用APDL語言,得到圖5模型的前四階理論固有頻率(f1i、f2i、f3i、f4i,i=1,2,…,120)和理論阻尼比(ζ1i、ζ2i、ζ3i、ζ4i,i=1,2,…,120)。該方法無需人工干預,可自動完成神經網絡訓練樣本數據的采集。

圖3 動態特性參數優化識別流程圖

圖4 機床滑動結合面實驗臺

圖5 機床的有限元模型

根據圖2的建模原理,采用兩個3層BP神經網絡模型。剛度參數神經網絡模型的輸入為滑動結合面剛度參數kx、ky、kz,輸出為圖5模型的前四階固有頻率f1、f2、f3、f4;阻尼參數神經網絡模型的輸入為滑動結合面阻尼參數cx、cy、cz,輸出為圖5模型的前四階阻尼比ζ1、ζ2、ζ3、ζ4。神經網絡隱含層節點個數根據經驗公式采用6個神經元,隱含層的激勵函數采用sigmoid函數,輸出層的激勵函數采用purelin函數。為提高BP神經網絡的訓練速度與精度,網絡訓練前,采用遺傳算法(GA)對BP神經網絡的初始權值和閾值進行優化[17],遺傳算法的種群個體數和遺傳代數的設定可通過多次反復調整確定,種群個體數設為50,遺傳代數設為250。交叉概率、變異概率的確定可參照已有類似文獻的研究結果,確定交叉概率取0.85,變異概率取0.05。采用實數編碼方法,通過MATLAB的遺傳算法工具箱進行求解,將得到的最優解作為BP網絡的初始權值和閾值。

2.2.2 網絡訓練與模型檢驗

利用上述120組樣本數據分別對圖2中的2個神經網絡進行訓練,網絡訓練前對樣本進行歸一化處理,將輸入-輸出參數轉化為[0,1]區間的值。對網絡進行訓練的過程是對一個非線性映射逼近的過程,也是權值和閾值不斷調整優化的過程,為了說明本文采用GA優化網絡初始權值和閾值的必要性,本文比較了GA-BP和BP神經網絡模型的訓練誤差曲線,如圖6所示。從圖6可知:GA-BP網絡在迭代1 000步達到設定誤差,低于BP網絡的3 210步,收斂速度明顯提升。在網絡輸入變量的變化范圍內,隨機另選非訓練樣本點的其他數值,對BP及GA-BP剛度參數神經網絡模型的精度進行測試,測試結果,如表1所示。從表1可知,GA-BP神經網絡模型預測值與有限元計算值更為接近。對BP和GA-BP模型在相同的目標精度下的訓練精度進行對比結果,如表2所示。綜上可知,GA-BP模型的訓練精度更高,且GA-BP神經網絡模型具有更強的泛化能力。

(a) BP神經網絡模型

(b) GA-BP神經網絡模型

2.3 實驗模態分析

對機床工作臺-床身滑動導軌進行模態試驗,確定機床的實驗固有頻率、阻尼比、模態振型。模態實驗的工作原理,如圖7所示。采用NI9234采集模塊的NIcDAQ-9172數據采集系統和ModalView模態分析軟件,通過Kistler SN2075427型激振力錘在工作臺一側施加激振力,分布在工作臺上的BK4525B型三向加速度傳感器測量工作臺振動信號。實驗測試現場,如圖8所示。

表1固有頻率計算值與預測值絕對誤差

Tab.1Absoluteerrorbetweenthecalculationandpredictedvaluesofnaturalfrequencies

固有頻率有限元計算值/HzBP預測值/Hz絕對誤差GA?BP預測值/Hz絕對誤差f141.06241.1580.09641.0700.00839.45239.323-0.12939.442-0.01036.79836.9340.13636.8090.01140.38940.287-0.10240.381-0.00834.57334.7140.14134.5840.011f258.42057.4530.03358.4220.00258.30958.272-0.03758.306-0.00356.41056.4550.04556.4120.00258.52358.485-0.03858.522-0.00157.82957.8830.05457.8330.004f3109.783109.8640.081109.7880.005109.471109.365-0.106109.463-0.008106.095106.1800.085106.1010.006109.676109.581-0.095109.671-0.005107.559107.6920.133107.5680.009f4126.827126.9400.113126.8340.007125.836125.699-0.137125.828-0.008123.203123.3310.128123.2100.007126.552126.446-0.106126.546-0.006121.667121.8600.193121.6760.009

表2 BP和GA-BP模型訓練性能對比

圖7 模態試驗工作原理

在Modal View軟件中建立工作臺的模型,如圖9所示。在該模型上設置了37個測點,使用單點激振多點拾振的錘擊法對結合面進行模態測試,床身通過地腳螺栓固定,螺栓的預緊力矩很大,可以認為床身與地面的連接方式為剛性連接。實驗模態分析結果,如表3所示。實驗振型,如圖10所示。

圖8 實驗測試現場

2.4 參數識別

由于剛度參數影響結合部結構的固有頻率、模態振型和振動幅值,阻尼參數對結構的固有頻率和模態振型影響較小,而對結構振動幅值影響較大。因此本文首先對結合面的剛度參數進行識別,然后基于已識別的剛度參數,對結合面的阻尼參數進行識別。

圖9 試驗模態測試模型及測點布置

模態階次i固有頻率fi/Hz阻尼比ζi/%136.894.23258.066.223108.295.024123.575.17

(a)第1階振型

(b)第2階振型

(c)第3階振型

(d)第4階振型

對剛度參數進行識別時,定義設計變量kx、ky、kz的變化范圍為(0.5~6.0)kN/mm,根據式(1)所示的目標函數及約束條件,采用布谷鳥算法進行尋優求解,優化識別過程的迭代收斂曲線,如圖11所示。識別的剛度參數,如表4所示。

圖11 剛度參數識別收斂曲線

阻尼參數識別時將已經識別的滑動結合面的剛度參數作為已知量,定義設計變量cx、cy、cz的變化范圍為(1.0~2.8)N·s/mm,根據式(2)所示的目標函數及約束條件,采用布谷鳥算法進行尋優求解,優化識別的迭代收斂曲線,如圖12所示。識別的阻尼參數值,如表5所示。

表4 剛度參數的識別結果

圖12 阻尼參數識別收斂曲線

Fig.12 The convergence curve of damping parameter identification

2.5 參數識別誤差的檢驗

根據上述優化所得的結合面參數,對圖5有限元模型中工作臺-床身滑動導軌結合面的剛度及阻尼參數進行設置,進行理論模態分析,并與實驗模態分析所得的固有頻率和阻尼比進行對比,如表6、表7所示。

表5 阻尼參數的識別結果

表6 實驗固有頻率與理論固有頻率的對比

由表6可知,理論模態固有頻率和實驗模態固有頻率的誤差<3.5%。由表7可知,理論模態分析的阻尼比和實驗測試值誤差<3%。剛度及阻尼參數的識別誤差均低于已有文獻的研究,驗證了本文提出的剛度、阻尼參數識別的正確性。考慮文中主要研究工作臺滑動導軌,截取工作臺理論模態振型,如圖13所示。與圖10試驗模態振型對比可得,理論振型與實驗振型符合較好。

表7 實驗阻尼比與理論阻尼比的對比

(a)第1階振型

(b)第2階振型

(c)第3階振型

(d)第4階振型

3 結 論

本文基于有限元分析數據,建立了滑動結合面動態特性參數的神經網絡模型,并將該模型作為其有限元分析的替代模型,結合機床實驗模態分析結果,優化識別出滑動結合面的剛度與阻尼參數。

以自行研制的機床滑動結合面實驗臺的工作臺-床身滑動導軌結合面為例,進行了建模、實驗、參數識別等分析,驗證了本文方法的有效性,參數識別誤差優于已有研究的結果。

本文方法適用于滑動結合面的同時,也適用于機床滾動導軌,螺栓,軸承,絲杠等結合面的動態特性參數識別,可為機床的整機理論分析提供準確的結合面參數。

[1] 趙宏林,丁慶新,曾鳴,等. 機床結合部特性的理論解析及應用[J]. 機械工程學報,2008, 44(12): 208-214.

ZHAO Honglin, DING Qingxin, ZENG Ming, et al. Theoretic analysis on and application of behaviors of machine tool joints[J]. Chinese Journal of Mechanical Engineering, 2008, 44(12): 208-214.

[2] BURDEKIN M, BACK N, COWLEY A. Analysis of the local deformation in machine joints[J]. ARCHIVE: Journal of Mechanical Engineering Science, 1979, 21: 25-32.

[3] VAFAEI S, RAHNEJAT H, AINI R. Vibration monitoring of high speed spindles using spectral analysis techniques[J]. International Journal of Machine Tools & Manufacture, 2002, 42(11): 1223-1234.

[4] 李玲,蔡安江,蔡力鋼. 栓接結合部動態特性參數辨識新方法[J]. 振動與沖擊,2014, 33(14): 15-19.

LI Ling, CAI Anjiang, CAI Ligang. New method to identify dynamic characteristics of bolted joints[J]. Journal of Vibration and Shock, 2014, 33(14): 15-19.

[5] 劉耀,黃玉美. 機床滾珠導軌中圓柱面-球面結合面靜特性分析及試驗研究[J].機械工程學報,2013, 49(21): 25-30.

LIU Yao, HUANG Yumei. Theoretical analysis and experimental study on static characteristics of the cylindrical-spherical joint surfaces of linear ball guide on machine tool[J]. Chinese Journal of Mechanical Engineering, 2013, 49(21): 25-30.

[6] CHLEBUS E, DYBALA B. Modeling and calculation of properties of sliding guideway[J]. International Journal of Machine Tools and Manufacture, 1999, 39(12): 1823-1839.

[7] 王禹林,吳曉楓,馮虎田. 基于結合面的大型螺紋磨床整機靜動態特性優化[J]. 振動與沖擊,2012, 31(20): 147-152.

WANG Yulin, WU Xiaofeng, FENG Hutian. Static and dynamic characteristics optimization for a whole large-sized thread grinder based on joint surface[J]. Journal of Vibration and Shock, 2012, 31(20): 147-152.

[8] 郭成龍,袁軍堂,汪振華,等.貼塑導軌結合面法向動態特性參數的識別及分析[J].振動與沖擊,2013, 32(3): 95-98.

GUO Chenglong, YUAN Juntang, WANG Zhenhua, et al. Identification and analysis for dynamic characteristic parameters of a plastic-coated guide’s joint[J]. Journal of Vibration and Shock, 2013, 32(3): 95-98.

[9] 伍良生,王澤林,屈重年,等.滑動結合面動態特性測試系統設計[J].機械設計與制造,2012(7): 9-11.

WU Liangsheng, WANG Zelin, QU Chongnian, et al. Design of dynamic characteristic testing system for sliding joint[J]. Machinery Design & Manufacture, 2012(7): 9-11.

[10] LEE W J, KIM S I. Joint stiffness identification of an ultra-precision machine for machining large-surface micro-feature[J]. International Journal of Precision Engineering and Manufacturing, 2009, 10(5): 115-121.

[11] 王立華,羅建平,劉泓濱,等.銑床關鍵結合面動態特性研究[J].振動與沖擊,2008, 27(8): 125-129.

WANG Lihua, LUO Jianping, LIU Hongbin, et al. Research on dynamic characteristics of key machine joint surfaces of the numerically controled milling machine[J]. Journal of Vibration and Shock, 2008, 27(8): 125-129.

[12] 孫明楠,米良,干凈,等.數控機床導軌結合部動態特性參數優化識別方法研究[J].四川大學學報(工程科學版),2012, 44(3): 217-223.

SUN Mingnan, MI Liang, GAN Jing, et al. An optimum identification method of dynamic characteristic parameters of guideway joint on a NC machine tool[J]. Journal of Sichuan University(Engineering Science Edition), 2012, 44(3): 217-223.

[13] DENG Congying, YIN Guofu, FANG Hui, et al. Dynamic characteristics optimization for a whole vertical machining center based on the configuration of joint stiffness[J]. International Journal of Advanced Manufacturing Technology, 2015, 76(5/6/7/8):1225-1242.

[14] 汪中厚,袁可可,李剛. 基于響應面法的滑動結合面動態特性參數優化識別[J]. 中國機械工程,2016, 27(5): 622-626.

WANG Zhonghou, YUAN Keke, LI Gang. Optimization identification for dynamic characteristic parameters of sliding joints based on response surface methodology[J]. China Mechanical Engineering, 2016, 27(5): 622-626.

[15] YANG X, DEB S. Cuckoo search via lévy flights[C]//World Congress on Nature & Biologically Inspired Computing. Piscataway:IEEE Publications, 2009: 210-214.

[16] 廖伯瑜,周新民,尹志宏. 現代機械動力學及其工程應用[M]. 北京:機械工業出版社,2004.

[17] 李松,劉力軍,解永樂.遺傳算法優化BP神經網絡的短時交通流混沌預測[J]. 控制與決策,2011, 26(10): 1581-1585.

LI Song, LIU Lijun, XIE Yongle. Chaotic prediction for short-term traffic flow of optimized BP neural network based on genetic algorithm[J]. Control and Decision, 2011, 26(10): 1581-1585.

猜你喜歡
模態有限元優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
磨削淬硬殘余應力的有限元分析
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 国产精品免费福利久久播放| 国产在线一二三区| 国产精品网址在线观看你懂的| 2020精品极品国产色在线观看 | 久久黄色小视频| 欧美在线黄| 浮力影院国产第一页| 不卡视频国产| 九九精品在线观看| 日本黄色a视频| 国产一区二区三区日韩精品| 久久精品人人做人人爽97| 看av免费毛片手机播放| 精品亚洲欧美中文字幕在线看| 欧美日韩一区二区三区在线视频| 四虎亚洲精品| 91视频99| 国产在线一区二区视频| 亚洲国产精品成人久久综合影院| 国产精品毛片一区| 538精品在线观看| 草草线在成年免费视频2| 日韩精品专区免费无码aⅴ| 青青草原国产免费av观看| 免费观看国产小粉嫩喷水| 中文字幕日韩视频欧美一区| 久久青青草原亚洲av无码| 免费人成在线观看成人片| 成人字幕网视频在线观看| 免费国产好深啊好涨好硬视频| 日韩欧美国产精品| 久久99国产综合精品女同| 久久国产精品影院| 国产高清自拍视频| 四虎成人免费毛片| 日本国产精品一区久久久| 日韩免费无码人妻系列| 国产黄网永久免费| 免费欧美一级| 91精品啪在线观看国产60岁| 日韩国产欧美精品在线| 国产国模一区二区三区四区| 日韩午夜福利在线观看| 欧美精品一区在线看| 国产精品无码AV片在线观看播放| 97精品国产高清久久久久蜜芽| 亚洲欧美国产视频| 亚洲日韩在线满18点击进入| 黄网站欧美内射| 精品福利视频网| 亚洲日韩精品无码专区| 国产杨幂丝袜av在线播放| 亚洲黄色成人| 天堂岛国av无码免费无禁网站 | 在线国产毛片手机小视频| 日韩欧美网址| 国产剧情无码视频在线观看| 国产欧美日韩另类| 国产综合网站| 毛片手机在线看| 国产日韩欧美视频| 亚洲美女一区| 国产日韩欧美视频| jizz在线免费播放| 影音先锋亚洲无码| 国产精品免费p区| 九色视频最新网址| 国产91熟女高潮一区二区| 亚洲男人的天堂久久香蕉| av免费在线观看美女叉开腿| 国产裸舞福利在线视频合集| 丰满少妇αⅴ无码区| 在线观看91香蕉国产免费| 免费不卡在线观看av| 蜜桃视频一区二区三区| 无码人妻热线精品视频| 91亚洲精品第一| 国产美女丝袜高潮| 呦视频在线一区二区三区| 亚洲区视频在线观看| 亚洲系列中文字幕一区二区| 九色视频在线免费观看|