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

基于接觸模型和峭度最優(yōu)Laplace小波的滾動軸承量化診斷?

2012-12-26 09:08:22朱永生張優(yōu)云
振動工程學報 2012年6期
關鍵詞:振動模型

袁 幸,朱永生,張優(yōu)云,洪 軍

(1.西安交通大學潤滑理論及軸承研究所,陜西 西安 710049;2.西安交通大學機械制造系統(tǒng)工程國家重點實驗室,陜西西安 710049)

引 言

飛機、電力機車、高檔機床等重大裝備的核心技術問題之一就是軸承,它對整個裝備制造業(yè)的發(fā)展水平有著舉足輕重的作用,其運行狀態(tài)往往直接影響整臺機器的精度、可靠性及壽命[1]。 Sawalhi N將最小熵盲反卷積和譜峭度應用于滾動軸承故障檢測中[2]。Karthik K應用互信息進行了故障特征選取[3]。譚繼勇針對滾動軸承沖擊信號的檢測問題,構造了隨機共振自適應檢測算法,提高了信噪比,避免了漏峰現(xiàn)象[4]。文獻[5]提出了基于經(jīng)驗模態(tài)分解和奇異值差分譜的軸承故障診斷方法。從已有研究結果來看,以先進信號分析為基礎的滾動軸承診斷技術研究較多,人們大都針對損傷信號處理與分類,研究相關的算法,卻很少關注服役軸承具體的損傷尺寸。對于滾動軸承這樣的基礎部件,應用廣泛但規(guī)格種類繁多,工況差別大,不同損傷程度的樣本獲取困難。以機車為例,輪對軸承的故障樣本少,一個輪對檢修流水線在3個月時間檢測的1 200多個輪對軸承中僅有2例故障軸承[6]。再者,單純的基于信號處理的傳統(tǒng)方法較少考慮軸承系統(tǒng)的具體結構特點,側重于對故障存在性的檢測,不易描述損傷演化過程,因而難以準確評價軸承可以持續(xù)運行的安全余量。

本文系統(tǒng)地提出根據(jù)正、反問題進行損傷量化診斷的原理及方法。建立多體接觸振動模型來描述滾動軸承空間運動特性,涵蓋軸承界面復雜相互作用及摩擦學系統(tǒng)行為,以接觸表面局部區(qū)域尺寸變化來刻畫損傷發(fā)展趨勢,實現(xiàn)復雜工況下軸承的損傷建模。通過精細劃分損傷區(qū)域由正問題求解標準模式數(shù)據(jù)庫,將反問題計算轉(zhuǎn)化為多維優(yōu)化問題,運用灰色關聯(lián)度辨識軸承損傷位置和具體尺寸。

1 滾動軸承損傷接觸建模

1.1 接觸剛度系數(shù)

由 Hertzian理論,點接觸負荷-變形關系為[7]

式中 W為彈性趨近量,Eschmann給出了兩固體點接觸彈性變形表達式

式中 Δ,kb,_是由材料屬性決定的 Hertzian系數(shù);d分別為彈性模量、泊松比與接觸表面曲率和。點接觸的接觸力可表示為

因此軸承接觸剛度系數(shù)為

∑ d的計算方法由 Harris給出,詳見文獻 [7],一旦軸承設計尺寸和材料確定,K可由式(2)~(4)求出。

1.2 接觸關系的坐標表達

圖1是滾動軸承坐標示意圖,第i個鋼球-套圈接觸變形Wi是內(nèi)圈在X,Y方向位移(xs,ys),鋼球位置角θi和游隙c的函數(shù),如下式

設(xb,yb)為鋼球坐標,由于振動傳遞作用,考慮到鋼球自身的振動,局部接觸變形變?yōu)?/p>

θi表示第i個鋼球的位置角

式中 kc為軸承公轉(zhuǎn)速度,N為鋼球個數(shù)。設轉(zhuǎn)速為k,則,其中Db為鋼球直徑,Dp為軸承節(jié)圓直徑。

圖1 滾動軸承示意圖Fig.1 Schematic diagram of rolling bearing

設(xo,yo)為外圈處的運動坐標,接觸變形同樣可表示為

1.3 動態(tài)接觸建模

球軸承接觸阻尼隨平均負荷變化,其值很小,Kramer通過頻域測試方法給出了阻尼近似計算公式[8]

故xb方向包含阻尼的動態(tài)接觸力可寫成

1.4 拖動力模型

點線接觸中兩個滾動元件之間的相對滑動速度引起的潤滑劑與滾動元件之間的摩擦力被稱為拖動力。拖動特性是決定軸承元件動態(tài)特性的關鍵因素之一,是滾動軸承動力學分析中不可缺少的基本參數(shù)。引入摩擦系數(shù)_,yb和xb方向拖動力的函數(shù)關系為

摩擦系數(shù)_由彈流潤滑模型決定,設Λ為油膜厚度與表面粗糙度的比值,則[7]

式中 _bd受潤滑參數(shù)及溫度的影響,在本文使用的潤滑模型中,摩擦系數(shù)_bd是滑 /滾率s的函數(shù)。

通過以上分析,對深溝球軸承的平面運動,考慮到剛度、阻尼、拖動力和滑 /滾率的動態(tài)接觸力表達式為

1.5 運動方程式

滾動軸承的振動存在著振動源和振動傳遞路徑的交互,如圖2所示,系統(tǒng)響應最終為傳感器接收。

圖2 滾動軸承系統(tǒng)振動模型Fig.2 Vibration model of rolling bearing system

軸承表面損傷形態(tài)決定了激振力的頻譜,界面動態(tài)接觸決定了振動系統(tǒng)的傳遞特性,最終的振動頻譜由上述二者共同決定。正問題分析中,外圈質(zhì)量和剛度將激勵系統(tǒng)較低的固有頻率。考慮到軸承的高頻諧振頻帶,設置mp,cp,kp為附加自由度來模擬軸承-軸承座結合部的約束以激勵高頻響應。將傳遞路徑的質(zhì)量、剛度、阻尼及幾何位置作為系統(tǒng)參數(shù),根據(jù)牛頓定律可建立內(nèi)圈-鋼球-外圈接觸振動微分方程組如下

式中ms,mb,mo分別為內(nèi)圈與軸的質(zhì)量,鋼球質(zhì)量,外圈質(zhì)量;ko,co為外圈剛度與阻尼;fx,fy為徑向載荷與不平衡力的總和。由于工程環(huán)境、安裝條件等因素的影響,使得軸承-軸承座結合部的材料特性、幾何尺寸和邊界條件具有微小的偏差。為了與實際工況等效,mp,cp,kp的具體數(shù)值可通過健康信號的頻譜加以修正。模型參數(shù)導致的誤差做了誤差置零處理[9],獲得數(shù)值計算與實驗結果一致的數(shù)值模型。

1.6 表面損傷建模

滾動軸承由疲勞誘發(fā)的表面損傷主要包括裂紋、凹點、剝落等。最早研究軸承損傷振動模型的是McFadden和 Smith,該模型通過軸承運轉(zhuǎn)過程中鋼球經(jīng)過故障位置產(chǎn)生的一系列沖擊與沖擊衰減函數(shù)的卷積來模擬振動[10]。Choudhury A提出外加沖擊力的形式模擬損傷所激勵的動態(tài)效應,建立了滾動軸承表面損傷的隨機振動數(shù)學模型[11]。Behzad M等通過基于模型的振動響應與傳遞函數(shù)相結合的方法較好地模擬了軸承的振動行為[12]。

到目前為止,提出的大多數(shù)為滾動軸承損傷的數(shù)字信號仿真與模擬沖擊模型,未能有效地將損傷融入到滾動軸承自身的幾何學和動力學方程中。本文通過找出滾動軸承在運轉(zhuǎn)過程中鋼球通過損傷區(qū)的運動變化規(guī)律,在宏觀上建模并與接觸振動相結合為分析滾動軸承疲勞累積損傷提供一種可用的方法,損傷模型如圖 3所示。

損傷長度為L,深度為h,當鋼球運行至損傷區(qū)域時,游隙值將隨之增加h,使得W=W+h,瞬間通過長度L后,游隙值恢復為W。當鋼球處于損傷區(qū)域時,導致接觸面間的 Hertzian接觸力突然降低或變?yōu)榱?相當于滾動軸承的非線性時變系統(tǒng)的激勵發(fā)生突變使得振動能量出現(xiàn)波動,損傷程度可用L×h或L/h來定量衡量。

圖3 表面損傷實物與幾何模型Fig.3 Physical and geometric model of surface damage

2 數(shù)值計算與分析

式(15)所代表的狀態(tài)方程組非線性強,難以解析,本文采用數(shù)值方法獲取系統(tǒng)響應。在MATLAB中編制了相應程序,計算方法為定步長 Runge-Kutta法。算例采用美國凱斯西儲大學(CWRU)軸承研究中心的深溝球軸承 (6205-2RS SKF),轉(zhuǎn)速1 797 r/min,采樣頻率 12 000 Hz,CWRU提供了軸承詳細設計參數(shù)[13]。軸承受到來自轉(zhuǎn)子不平衡激勵所產(chǎn)生的強迫振動,其振動頻率為轉(zhuǎn)子的旋轉(zhuǎn)頻率。軸承也將產(chǎn)生由于總剛度連續(xù)周期變化而形成的VC(Varying compliance)振動及其諧波組合[14,15]。這是因為,在徑向載荷的作用下,各鋼球的受力情況是不一樣的,且隨著鋼球上的某一點的運動位置不同受力情況亦不一樣,隨著鋼球相對于徑向載荷作用線的移動,軸承剛度以數(shù)倍于鋼球沿靜止套圈轉(zhuǎn)動的頻率呈周期性變化。這種振動為參數(shù)激勵的通過振動,VC振動頻率與外圈損傷頻率數(shù)值上相等屬于軸承自身本體振動成分,VC=kc×nb=106.7 Hz。圖4的功率譜清晰地顯示了VC振動及其高次諧波,初步表明文中計算具有一定的可靠性。

圖4 正常狀態(tài)波形及功率譜圖Fig.4 Waveform and spectrum of the normal bearing

圖5~7分別是損傷信號波形及功率譜圖,內(nèi)、外圈和鋼球損傷頻率依次為fi=161.3 Hz,fo=106.7 Hz,fb=140.4 Hz。與圖 4比較,損傷時波形具有明顯的周期沖擊特征,這是游隙瞬間變化引起的,波形重復出現(xiàn)的時間間隔為損傷出現(xiàn)的周期0.006 2,0.009 4,0.007 1 s。

由于接觸振動模型全面考慮了振動的傳遞路徑,誘發(fā)軸承各階固有頻率,使得高頻諧振頻率與損傷頻率調(diào)制,在高頻處譜圖出現(xiàn)以損傷頻率為主的邊頻帶,幅值遠離載波中心程逐漸減小的趨勢。由于內(nèi)圈損傷位置隨軸旋轉(zhuǎn)不斷改變,損傷引起的沖擊響應幅值必然受到轉(zhuǎn)頻調(diào)制,會在損傷頻率各次諧波處出現(xiàn)邊頻譜線,造成內(nèi)圈損傷功率譜圖更加復雜,如圖5所示。試驗信號與仿真結果取得較好的一致性,也說明了本文建模方法的有效性。

圖5 內(nèi)圈損傷狀態(tài)波形及功率譜圖Fig.5 Waveform and spectrum of the bearing with inner race damage

圖6 外圈損傷狀態(tài)波形及功率譜圖Fig.6 Waveform and spectrum of the bearing with outer racedamage

圖7 鋼球損傷狀態(tài)波形及功率譜圖Fig.7 Waveform and spectrum of the bearing with rolling element damage

3 故障量化特征提取

3.1 故障演變過程描述

將L×h設為輕度、中度、嚴重 3個階段,分別取值為(0.1,0.1),(0.6,0.8),(1.2,1.5)mm。圖8中(a)~(c)為 3個階段的外圈損傷頻譜,可見譜圖結構差異明顯。輕微損傷時損傷頻率被諧波分量所干擾,低頻調(diào)制邊頻能量較弱,隨著損傷程度加深低頻調(diào)制譜線逐漸清晰,并且譜圖能量分布變化顯著。圖8(d)為損傷發(fā)展時歸一化振動烈度變化曲線,損傷加劇沖擊越明顯,振動烈度值越大;沖擊越不明顯,甚至被VC及其諧波頻率所淹沒時,振動烈度越小。不難看出,通過剖分L×h的數(shù)值可以有效區(qū)分振動響應。以L×h的大小量化損傷區(qū)域等效于局部不同柔度的植入,從而誘發(fā)振動響應變化,形成模式的多樣性。

圖8 外圈不同損傷程度的頻譜及振動烈度值Fig.8 Spectrum and vibration intensity of the bearing with outer race damagefor different fault severities

3.2 峭度最優(yōu)Laplace小波特征提取

機械診斷屬于模式識別問題,為了增強特征向量的敏感性,達到精確提取軸承損傷時沖擊響應信號的目的,本文研究了一種峭度最優(yōu) Laplace小波。R Lind等人構造的 Laplace小波是一種復指數(shù)小波,解析表達式為[16]

式中 參數(shù)矢量γ=[a,kn]決定了小波的性能;成員變量a,kn為模態(tài)動力學參數(shù);A用來歸一化小波函數(shù)。峭度對沖擊信息非常敏感,又是區(qū)分非高斯分布的指標,高的峭度值意味著信號含有豐富的沖擊成分。基于波形匹配的思想,本文用峭度最大化準則定量描述參數(shù)矢量γ。在一定范圍內(nèi)變化a和kn,選擇使得峭度最大的a,kn值為參數(shù),經(jīng)計算獲得最優(yōu)值a=0.85,kn=16,這就是峭度最優(yōu) Laplace小波變換。本質(zhì)上,基于峭度最優(yōu)的 Laplace小波變換也就是控制小波濾波器的通帶帶寬,參數(shù)優(yōu)化相當于調(diào)整品質(zhì)因子Q,最終結果使小波沖擊特征更加顯著,得到與信號最相似的基元函數(shù)。

時域峭度、脈沖等 5個指標已被證明可較好區(qū)分多重狀態(tài)[17]。僅依靠時域信息對軸承損傷尺寸的精確識別還不夠理想,因為信息含量太少,頻譜結構變化得不到體現(xiàn)。因此,結合頻域無量綱指標組建標準模式,頻域指標定義為平均頻率、頻譜偏斜度、頻譜峭度、均方根比[18]。由于工況變動的隨機性和采集信息的噪聲影響,量化診斷具有一定的不確定性。模糊分析是解決識別過程中不確定性問題的有力工具。以 [0,1]之間的模糊隸屬度表示特征向量有助于提升特征敏感度,優(yōu)化特征空間布局,增強信息可靠性。根據(jù)經(jīng)驗可選擇升半嶺型函數(shù)表示隸屬度[17]。至此,構造了量化診斷標準模式向量

以外圈損傷尺寸 0.001,0.002,0.003 mm為例,根據(jù)模型計算特征向量,結果如圖9所示。可見,經(jīng)峭度最優(yōu) Laplace小波提取的 9個無量綱指標對差值為 0.001 mm的損傷尺寸實現(xiàn)了良好的分類,也說明動力學響應參數(shù)可以較好映射局部損傷尺寸,為量化診斷奠定了基礎。

圖9 外圈不同損傷尺寸的特征向量Fig.9 Feature vectors of the bearing with outer race damage for different fault size

4 局部損傷的定量檢測

4.1 定量識別的灰色關聯(lián)度分析

量化診斷正問題實際上是通過對含有不同損傷位置和尺寸的滾動軸承進行動力學求解,獲得損傷特征并轉(zhuǎn)化為模糊向量,由計算得到的離散值構建損傷數(shù)據(jù)庫。本文將反問題轉(zhuǎn)化為軸承的實際振動輸出與接觸振動數(shù)值模型的數(shù)值輸出間差異最小化為優(yōu)化目標的優(yōu)化問題。利用關聯(lián)尋優(yōu)求解出與輸入值最相似的樣本點,進而測出損傷位置及損傷長度L,深度h。正問題計算得到標準模式數(shù)據(jù)庫

式中i為模式類別,k=9表示特征向量個數(shù)。令待檢模式為?j(k)={xj(1),xj(2),… ,xj(9)},其中j為模式類別,k的意義同上。記Si對Sj的關聯(lián)系數(shù)為Yij(k),Yij(k)由下式確定[19]

式中 Δij=|Si(k)-Sj(k)|稱為絕對差,d為分辨系數(shù),量化診斷時可取d=0.5,則Si與Sj的關聯(lián)度為

關聯(lián)度從物理意義講是對灰色過程的白化處理,從數(shù)學角度講是空間幾何曲線的比較。反問題辨識損傷尺寸可用如下數(shù)學關系式描述

正、反問題進行滾動軸承量化診斷的流程如圖10。

圖10 量化診斷流程Fig.10 Flow chart for the quantitative identification

4.2 案例驗證

案例同樣源自 CWRU軸承中心,采用電火花機床分別在測試軸承的內(nèi)、外圈和鋼球上加工深度為11 mils(0.279 4 mm)的盲孔,盲孔的直徑為 7 mils(0.177 8 mm),14 mils(0.355 6 mm),21 mils(0.533 4 mm),這樣就可以通過實驗數(shù)據(jù)完成量化診斷。測試軸承分別在 4種不同載荷(載荷分別為0~ 3馬力)以及相對應的 1 797,1 772,1 750,1 730 r/min四種轉(zhuǎn)速下測試。振動信號由16通道DAT數(shù)據(jù)記錄儀采集,采樣頻率為12 k Hz。 CWRU未提供第4種工況下的外圈數(shù)據(jù),因此本文用內(nèi)、外圈和鋼球在前 3種工況下的樣本各 20個按照圖 10的流程進行量化診斷。

正問題計算區(qū)間定為L=0.1~1 mm,h=0.1~1 mm,搜索步長0.001 mm。表1給出了任意一種工況下h的反問題分析結果,表 2為 3種工況下L的分析結果。計算出損傷尺寸預測值的相對誤差,結果表明,損傷深度預測相對誤差不超過 4%,損傷長度的相對誤差不超過 7%。

表1 量化診斷結果(h/mm)Tab.1 Results for the quantitative identification(h/mm)

表 2 量化診斷結果(L/mm)Tab.2 Results f or thequantitative identif ication(L/mm)

5 結 論

(1)滾動軸承的振動響應是軸承內(nèi)部幾何關系、運動關系在外部載荷作用下的綜合體現(xiàn)。結合軸承界面復合約束作用及摩擦學行為,建立了計及傳遞路徑的滾動軸承動態(tài)接觸振動模型,獲取了軸承不同損傷狀態(tài)下的本體振動信號。該模型可以敏銳地捕捉到軸承損傷時振動信號的沖擊信息且涵蓋了頻譜的中、高頻成分,理論結果與案例數(shù)據(jù)分析相吻合。

(2)局部損傷的出現(xiàn)引起內(nèi)部游隙的變化,進而影響到軸承的振動特性。不同損傷位置和尺寸形成了不同的損傷模式;反之,運用模式向量也可反求損傷位置和具體尺寸。

(3)峭度對信號中的沖擊成分非常敏感,同時峭度又是一種無量綱指標,它不受信號絕對水平的影響。根據(jù)特征波形匹配原理,研究峭度最優(yōu) Laplace小波復合特征提取方法,全面挖掘了振動信號時、頻域動態(tài)信息。

(4)將量化診斷中的反問題轉(zhuǎn)化為優(yōu)化辨識,以實測特征向量作為輸入,運用灰色關聯(lián)度辨識具體損傷尺寸。診斷實踐表明,識別誤差不超過7%,這為方法的進一步實際應用提供了有力的依據(jù)。

[1] Robert B Randall,Jére Antoni.Rolling element bearing diagnostics—— A tutorial[J].Mechanical Systems and Signal Processing,2011,25(2):485—520.

[2] Sawalhi N,Randall R B,Endo H.The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis[J].Mechanical Systems and Signal Processing,2007,21(6):2 616—2 633.

[3] Karthik K,Nataraj C.Feature selection for fault detection in rolling element bearings using mutual information[J].Journal of Vibration and Acoustics,2011,133(6):061001—11.

[4] 譚繼勇,陳雪峰,何正嘉.沖擊信號的隨機共振自適應檢測方法 [J].機械工程學報,2010,46(23):61— 67.Tan Jiyong,Chen Xuefeng,He Zhengjia.Impact signal detection method with adaptive stochastic resonance[J].Journal of Mechanical Engineering,2010,46(23):61—67.

[5] 張超,陳建軍,徐亞蘭.基于 EMD分解和奇異值差分譜理論的軸承故障診斷方法[J].振動工程學報,2011,24(5):539—545.Zhang Chao,Chen Jianjun,Xu Yalan.A bearing fault diagnosis method based on EMD and difference spectrum theory of singular value[J].Journal of Vibration Engineering,2011,24(5):539—545.

[6] 屈梁生.機械監(jiān)測診斷中的理論與方法 [M].西安:西安交通大學出版社,2009.Qu Liangsheng.Theory and Methodology for Mechanical Condition Monitoring and Fault Diagnosis[M].Xi'an:Xi'an Jiao Tong University Press,2009.

[7] Harris T A.Rolling Bearing Analysis[M].4th ed.New York:John Wiley and Sons,2001.

[8] Kramer E.Dynamic of Rotors and Foundations[M].Berlin:Springer,1993.

[9] 陳雪峰.小波有限元理論與裂紋故障診斷的研究[D].西安:西安交通大學,2004.Chen Xuefeng.Research on the Theory of Wavelet Finite Element Method and Cracks Fault Diagnosis[D].Xi'an:Xi'an Jiao Tong University,2004.

[10]Mc Fadden P D,Smith J D.Model for the vibration produced by a single point defect in a rolling element bearing[J].Journal of Sound and Vibration,1984,96(1),69—82.

[11]Choudhury A, Tandon N.Vibration response of rolling element bearings in a rotor bearing system to a local defect under radial load[J].ASME,Journal of Tribology,2006,128(2):252—261.

[12]Behzad M,Bastami A R,Mba D.A new model for estimating vibrations generated in the defective rolling element bearings[J].ASME,Journal of Vibration and Acoustics,2011,133(4):041011—8.

[13]The Case Western Reserve University Bearing Data Center Website.Bearing data center seeded fault test data[EB/OL].http:http://www.eecs/cwru/edu/laboratory/bearing.

[14]TIWARI M,GUPTA K,PRAKASH O.Effect of radial internal clearance of a ball bearing on the dynamics of a balanced horizontal rotor[J].Journal of Sound and Vibration,2000,238(5):723—756.

[15]TIWARIM,GUPTA K,PRAKASH O.Dynamic response of an unbalanced rotor supported on ball bearings[J].Journal of Sound and Vibration,2000,238(5):757—779.

[16]何正嘉,訾艷陽,張西寧.現(xiàn)代信號處理及工程應用[M].西安:西安交通大學出版社 ,2007.He Zhengjia,Zi Yanyang,Zhang Xining.Modern Signal Processing and Engineering Application[M].Xi'an:Xi'an Jiao Tong University Press,2007.

[17]徐敏,虞和濟,張瑞林.設備故障診斷手冊 [M].西安:西安交通大學出版社,1998.Xu Min,Yu Heji,Zhang Ruilin.Handbook of Machinery Fault Diagnosis[M].Xi'an:Xi'an Jiao Tong University Press,1998.

[18]Chen P,Masatoshi T,Toshio T.Fault diagnosis method for machinery in unsteady operating condition by instantaneous power spectrum and genetic programming[J].Mechanical Systems and Signal Processing,2005,19(1):175—194.

[19]鄧聚龍.灰色系統(tǒng)基本方法 [M].武漢:華中科技大學出版社,2005.Deng Julong.The Primary Methods of Grey System Theory[M].Wuhan: Huazhong University of Science& Technology Press,2005.

猜你喜歡
振動模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
重要模型『一線三等角』
This “Singing Highway”plays music
重尾非線性自回歸模型自加權M-估計的漸近分布
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 97久久人人超碰国产精品| 久草性视频| 免费人成视网站在线不卡| 日本在线亚洲| 国产精彩视频在线观看| 亚洲天堂精品视频| 538精品在线观看| 亚洲国产成人麻豆精品| 久久久噜噜噜| 精品国产一区二区三区在线观看| 亚洲最新地址| 激情视频综合网| 黄色福利在线| 成年人国产网站| 爱做久久久久久| 亚洲国产成人综合精品2020| 91色在线观看| 亚洲综合日韩精品| 呦女亚洲一区精品| 毛片免费观看视频| 日韩av电影一区二区三区四区| 成人午夜视频网站| 在线另类稀缺国产呦| 国产高清又黄又嫩的免费视频网站| 一级毛片不卡片免费观看| 午夜高清国产拍精品| 拍国产真实乱人偷精品| 一级片免费网站| 99这里只有精品免费视频| 国产制服丝袜91在线| 欧美三級片黃色三級片黃色1| 日韩欧美在线观看| 国产一区二区网站| 日韩精品一区二区三区中文无码| 久久国产精品77777| 久久综合成人| 国产欧美日韩另类精彩视频| 国产亚洲视频中文字幕视频| 久久精品中文字幕免费| 亚洲人妖在线| 中文天堂在线视频| 久久亚洲国产一区二区| 手机在线国产精品| 国产成人调教在线视频| 又黄又湿又爽的视频| 久久熟女AV| 国产免费福利网站| 日韩成人在线网站| 亚洲中文无码av永久伊人| 国产成人乱码一区二区三区在线| 久久成人免费| 国产无码精品在线| 麻豆精品国产自产在线| 国产成人AV大片大片在线播放 | 婷婷伊人五月| 色播五月婷婷| 91视频99| 欧美第一页在线| 正在播放久久| 国产swag在线观看| 日韩一二三区视频精品| 91精品日韩人妻无码久久| 香蕉精品在线| 国产欧美日韩视频怡春院| 亚洲精品欧美日本中文字幕| 亚洲一欧洲中文字幕在线| 欧美精品一区二区三区中文字幕| 国内精自视频品线一二区| 国产96在线 | 制服丝袜在线视频香蕉| 精品乱码久久久久久久| 亚洲精品国产自在现线最新| 久久九九热视频| 国产精品尤物在线| 国产成人精品18| 国产国拍精品视频免费看| 91在线播放免费不卡无毒| 久操线在视频在线观看| 亚洲最大在线观看| 国产成人精品一区二区免费看京| 99re在线观看视频| 国产乱人乱偷精品视频a人人澡|