亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        VMD的時(shí)頻特性及其在滾動(dòng)軸承故障診斷中的應(yīng)用

        2018-03-23 05:46:51陳保家劉浩濤邱光銀汪新波
        關(guān)鍵詞:故障診斷模態(tài)故障

        陳保家 聶 凱 劉浩濤 邱光銀 汪新波 林 云

        (1. 三峽大學(xué) 水電機(jī)械設(shè)備設(shè)計(jì)與維護(hù)湖北省重點(diǎn)實(shí)驗(yàn)室, 湖北 宜昌 443002;2. 三峽大學(xué) 機(jī)械與動(dòng)力學(xué)院, 湖北 宜昌 443002)

        滾動(dòng)軸承振動(dòng)信號(hào)中的非固有頻率沖擊成分代表故障征兆,對于早期故障,這種非固有頻率的沖擊往往很微弱,而且軸承發(fā)生故障時(shí),其振動(dòng)信號(hào)具有非線性、非平穩(wěn)特性[1],傳統(tǒng)的傅里葉變換以信號(hào)可表示為不同頻率的平穩(wěn)正弦波的線性迭加作為前提條件,但缺乏對信號(hào)局部時(shí)頻變化的精確描述.所以需要用時(shí)頻域分析方法對信號(hào)進(jìn)行診斷.

        經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)方法可將非線性、非平穩(wěn)信號(hào)自適應(yīng)分解成不同頻段的本征模態(tài)分量(Intrinsic Mode Function,IMF)進(jìn)而轉(zhuǎn)化成平穩(wěn)信號(hào)進(jìn)行分析[2-3],在滾動(dòng)軸承故障特征提取上得到了成功的應(yīng)用.Kedadouche等提出了一種基于EMD和Teager能量算子[4]結(jié)合最小熵反褶積的齒輪故障診斷方法,高強(qiáng)等利用EMD對滾動(dòng)軸承故障信號(hào)進(jìn)行后利用Hilbert變換進(jìn)行包絡(luò)解調(diào)得到包絡(luò)譜,成功地提取滾動(dòng)軸承內(nèi)外圈故障特征頻率[5].但是由于EMD存在一些明顯的缺點(diǎn),EMD分解的分量之間容易出現(xiàn)模態(tài)混疊以及EMD存在著明顯的端點(diǎn)效應(yīng),這大大制約著EMD在軸承故障診斷中的應(yīng)用[6].小波變換(Wavelet Transform)是一種時(shí)間窗和頻率窗都可改變的時(shí)頻局部化分析方法,它具有窗口大小固定但其形狀可改變的特性.但由于選取合適的小波基函數(shù)難度非常大,同時(shí)小波分析必須人為選擇分解層數(shù),缺乏對信號(hào)處理的自適應(yīng)性[7],從而顯得有很大的局限性.

        變分模態(tài)分解(Variational Mode Decomposition,VMD)是Konstantin等[8]于2014年提出的一種自適應(yīng)信號(hào)分解新方法,該方法在分解過程中通過循環(huán)迭代求取約束變分問題的最優(yōu)解來確定分解得到的IMF的頻率中心及帶寬,從而有效地分離出信號(hào)的各個(gè)頻率成分.與EMD和小波相比,VMD具有收斂快、魯棒性高的特點(diǎn),并在故障診斷領(lǐng)域得到了應(yīng)用,通過VMD分解得到的各個(gè)IMF進(jìn)行解調(diào)處理可有效提取故障特征,目前常用的方法有能量算子解調(diào)、希爾伯特包絡(luò)解調(diào)等.

        本文通過仿真信號(hào)對比分析了VMD、EMD和小波分解的噪聲魯棒性和抑制模態(tài)混疊能力,并將其應(yīng)用于滾動(dòng)軸承故障診斷實(shí)例當(dāng)中,結(jié)果顯示VMD在故障信號(hào)處理方面有優(yōu)越性.

        1 信號(hào)分解的仿真分析

        VMD的分解過程是一種自適應(yīng)求解約束性變分問題最優(yōu)解的過程,通過循環(huán)迭代確定不同的頻率中心及帶寬從而將信號(hào)分解成不同的IMF分量.詳細(xì)步驟可參考文獻(xiàn)[8-9,11,12].

        1.1 仿真信號(hào)對比分析

        為了研究VMD、EMD和小波變換的分析效果,采用以上3種方法分別對仿真信號(hào)進(jìn)行分解,仿真信號(hào)x(t)由3個(gè)簡諧信號(hào)加隨機(jī)噪聲組成,表達(dá)式為

        x(t)=f1(t)+f2(t)+f3(t)+n(t)(1)

        圖1 含噪聲信號(hào)的組成成分

        選用db10小波對仿真信號(hào)進(jìn)行3層小波分解,分解所得到的時(shí)域波形及其頻譜如圖2所示.從圖中可看出,經(jīng)過小波分解所得的各個(gè)分量中,沒有出現(xiàn)頻率為20 Hz的諧波分量,且分量中都存在一定的噪聲干擾.圖3為EMD分解結(jié)果,其中4個(gè)分量中已經(jīng)分離出了3個(gè)對應(yīng)的諧波分量,EMD分量中第4個(gè)分量和第3個(gè)分量對應(yīng)的是20 Hz和50 Hz的諧波分量,含有少量的噪聲成分,第1個(gè)分量對應(yīng)的是288 Hz的諧波分量,但是噪聲干擾嚴(yán)重,第2個(gè)分量是原始信號(hào)中沒有的偽分量,對比圖1中的原始分量波形,可以看出EMD分解的信號(hào)分量還存在不同程度的失真.總體來看,通過小波分解和EMD不僅很難辨識(shí)出原始信號(hào)中的諧波分量,而且噪聲魯棒性差.

        圖2 信號(hào)的小波分解

        圖3 信號(hào)的EMD分解

        同樣,利用VMD對x(t)進(jìn)行分析,設(shè)定模態(tài)數(shù)k=4,所得到的時(shí)域波形及其頻譜如圖4所示,VMD的前3個(gè)IMF分量將原始信號(hào)的3個(gè)諧波分量從低頻到高頻都一一分解出來了,其中第1個(gè)IMF分量和第2個(gè)IMF分量與20 Hz和50 Hz的成分幾乎相同,頻譜非常干凈無噪聲,第3個(gè)IMF分量含有少量的噪聲,第4個(gè)IMF分量為被分離出來的純噪聲.總的來看,VMD能夠?qū)⒃夹盘?hào)的各組成成分與噪聲成分進(jìn)行有效分離,具有很好的噪聲魯棒性,結(jié)合各IMF分量的時(shí)域圖和FFT譜,可以很容易推斷原始信號(hào)x(t)的組成成分.

        圖4 信號(hào)的VMD分解

        1.2 抗模態(tài)混疊性能研究

        為了研究VMD的抗模態(tài)混疊性能.本節(jié)選用簡諧信號(hào)與間斷信號(hào)進(jìn)行合成,用來模擬復(fù)雜信號(hào)中包含間斷信號(hào)的情況.圖5為信號(hào)各組成部分的時(shí)域圖,合成信號(hào)y(t)的表達(dá)式為

        圖5 間斷信號(hào)的各組成部分

        其中,f1(t)、f2(t)與式(7)中的f1(t)、f2(t)相同,諧波分量f3(t)是間斷信號(hào),其表達(dá)式為

        (3)

        選用db10小波對y(t)進(jìn)行3層小波分解,其結(jié)果如圖6所示.第3分量和第4分量沒有分解出頻率為20 Hz和50 Hz的諧波,而且第2分量和第1分量有一定的模式混淆,并不能有效地分解出仿真信號(hào)的各個(gè)成分,由此可以看出小波在分解間斷信號(hào)容易受到干擾.再利用EMD對y(t)進(jìn)行分析,分解得到7個(gè)分量和一個(gè)殘余信號(hào).為圖文簡潔,取前4個(gè)分量如圖7所示.

        圖6 信號(hào)的小波分解

        圖7 信號(hào)的EMD分解

        對比圖5和圖7可以看到,合成信號(hào)中頻率為50 Hz和頻率為250 Hz的成分被分解到IMF1中,而在IMF2中又出現(xiàn)了頻率為50 Hz的中頻成分,同樣頻率為20 Hz的低頻成分也分別被分解到IMF2和IMF3當(dāng)中,這說明經(jīng)驗(yàn)?zāi)B(tài)分解存在嚴(yán)重的模態(tài)混疊,不能將間斷信號(hào)準(zhǔn)確分離出來.同樣利用VMD對y(t)進(jìn)行分析,設(shè)定模態(tài)數(shù)k為3,結(jié)果如圖8所示,從圖中可以看出VMD方法準(zhǔn)確地將3個(gè)分量分離出來,特別是對于間斷分量f3(t)的分解,僅僅只是在信號(hào)突變處出現(xiàn)了比較輕微的誤差.通過對比可以發(fā)現(xiàn),VMD比小波變換和EMD具有更好的抗模態(tài)混疊效果.

        圖8 信號(hào)的VMD分解

        2 VMD滾動(dòng)軸承故障診斷流程

        通過以上分析,可以看出VMD方法摒棄了信號(hào)遞歸分解方式的束縛,采用一種非遞歸的、自適應(yīng)搜尋策略對信號(hào)進(jìn)行分解,因此相對于EMD和小波分解具有更好的噪聲魯棒性[9]和抗模態(tài)混疊效果,適合于從較低信噪比的軸承故障信號(hào)中將有用的特征信號(hào)分離出來,實(shí)現(xiàn)軸承的故障診斷.峭度指標(biāo)是無量綱參數(shù),由于它與軸承轉(zhuǎn)速、尺寸、載荷等無關(guān),對沖擊信號(hào)特別敏感,當(dāng)故障信號(hào)含有沖擊信號(hào)時(shí),峭度值會(huì)成明顯的變化.本文采用VMD、Teager能量算子解調(diào)和峭度指標(biāo)分量優(yōu)選方法相結(jié)合,對滾動(dòng)軸承故障信號(hào)進(jìn)行分析,提取有效的故障特征,由于文獻(xiàn)的篇幅限制,在此不再對Teager能量算子解調(diào)的基本原理進(jìn)行展開敘述,詳細(xì)的基本原理可以參考文獻(xiàn)[4].其診斷步驟如下:

        1)獲取振動(dòng)信號(hào),初始化模態(tài)數(shù),令k=2.

        2)對信號(hào)進(jìn)行VMD分解,觀察每個(gè)模態(tài)分量的中心頻率.

        3)判斷中心頻率是否相近,如果相近,則確定模態(tài)數(shù)k=k-1,否則以模態(tài)數(shù)k=k+1進(jìn)行步驟2),直到模態(tài)數(shù)設(shè)置合理.

        4)通過峭度準(zhǔn)則[10]對IMF分量進(jìn)行篩選,選取出峭度值最大的模態(tài),將其做為故障特征最敏感分量.

        5)利用Teager能量算子解調(diào)[11]的方法對步驟4)中優(yōu)選出的分量進(jìn)行包絡(luò)解調(diào)分析,求取分量包絡(luò)譜.

        6)將包絡(luò)譜與故障特征頻率進(jìn)行匹配,判斷故障類型.

        3 診斷實(shí)例

        3.1 滾動(dòng)軸承故障數(shù)據(jù)獲取

        為了進(jìn)一步驗(yàn)證該方法對滾動(dòng)軸承故障診斷的有效性和實(shí)用性,文中選取美國Case Western Reserve University電氣工程實(shí)驗(yàn)室軸承故障實(shí)驗(yàn)臺(tái)采集的軸承故障振動(dòng)信號(hào)進(jìn)行研究分析.

        試驗(yàn)臺(tái)模擬滾動(dòng)軸承4種狀態(tài):內(nèi)圈故障、外圈故障、滾動(dòng)體故障以及正常狀態(tài).本論文選取風(fēng)機(jī)端的數(shù)據(jù)作為研究對象,軸承型號(hào)6205-2RS JEM,采樣頻率12 kHz,軸承轉(zhuǎn)速r=1 797 r/min,轉(zhuǎn)頻為29.95 Hz.通過滾動(dòng)軸承不同部件故障的特征頻率計(jì)算公式可得:外圈故障頻率fo=107.34 Hz,內(nèi)圈故障頻率fi=162.18 Hz,滾動(dòng)體故障頻率fb=119.42 Hz.

        3.2 內(nèi)圈故障分析

        滾動(dòng)軸承發(fā)生內(nèi)圈故障時(shí),滾動(dòng)體經(jīng)過故障點(diǎn)處會(huì)引起沖擊振動(dòng),內(nèi)圈轉(zhuǎn)動(dòng)時(shí)將產(chǎn)生周期性沖擊[12],內(nèi)圈故障信號(hào)的時(shí)域和頻域波形圖如圖9所示.先采用db10小波對內(nèi)圈故障信號(hào)做4層分解,分解結(jié)果如圖10所示.

        圖9 內(nèi)圈故障的時(shí)頻圖

        圖10 信號(hào)的小波分解

        從小波分量中并不能看出明顯的故障特征,再根據(jù)峭度準(zhǔn)側(cè)選取峭度系數(shù)最大的分量為d1,對其進(jìn)行Teager能量算子解調(diào),其包絡(luò)譜如圖11所示.從圖中可以看到內(nèi)圈故障特征頻率162 Hz和2倍頻成分以及轉(zhuǎn)頻成分.接著采用EMD方法對內(nèi)圈故障信號(hào)進(jìn)行分解,共得到10個(gè)分量和一個(gè)殘余分量,為了描述方便,給出了前6個(gè)IMF分量時(shí)域波形,如圖12所示,同樣峭度系數(shù)最大原則選取出最優(yōu)分量為IMF1,對其進(jìn)行Teager能量算子解調(diào),其包絡(luò)譜如圖13所示,同樣可以觀察到內(nèi)圈故障特征頻率和2倍頻成分,相比小波來說,非故障特征頻率處的幅值較小,故障特征提取能力稍強(qiáng),但是轉(zhuǎn)頻部分卻沒有凸顯出來.

        圖11 小波分量能量算子解調(diào)

        圖12 信號(hào)的EMD分解

        圖13 EMD分量能量算子解調(diào)包絡(luò)

        最后采用本文的VMD方法對內(nèi)圈故障信號(hào)進(jìn)行分析,根據(jù)VMD算法提取信號(hào)特征最優(yōu)原則,盡可能避免模態(tài)混疊同時(shí)保留信號(hào)的全部特征信息的情況下做了模態(tài)分量個(gè)數(shù)的對比分析[13],當(dāng)k取8的時(shí)候,分解得到的模態(tài)分量頻帶沒有出現(xiàn)混疊,且分解得到的各分量頻譜的幅值稍微偏高,特征信息頻率的周圍含有相對較少的噪聲成分,因此在頻域中特征信息部分更為凸顯,由此可選擇k=8,其分解結(jié)果如圖14所示.利用峭度系數(shù)最大原則選出的最優(yōu)分量為IMF4,對其進(jìn)行Teager能量算子解調(diào),得到包絡(luò)譜如圖15所示,圖中可明顯觀察到軸承內(nèi)圈故障特征頻率162 Hz相近頻率,且對比度很高,并且凸顯出了軸承30 Hz的轉(zhuǎn)頻成分,有效地識(shí)別了故障特征頻率,和內(nèi)圈故障特性頻率吻合,相比小波和EMD對噪聲干擾性更小,但是這樣并不能完全體現(xiàn)VMD的優(yōu)勢,接下來將在故障信號(hào)更微弱的滾動(dòng)體實(shí)驗(yàn)里做分析,因?yàn)闇y得滾動(dòng)體振動(dòng)信號(hào)的傳遞路徑復(fù)雜,且失效形式種類多樣,因此故障特征更難被識(shí)別.

        圖14 信號(hào)的VMD分解

        圖15 VMD分量能量算子解調(diào)包絡(luò)

        3.3 滾動(dòng)體故障分析

        當(dāng)滾動(dòng)體表面出現(xiàn)局部損傷時(shí),滾動(dòng)體在深溝球軸承的內(nèi)、外圈滾道之間除了存在滾動(dòng)狀態(tài)之外,還可能存在滑動(dòng)狀態(tài),這種滑動(dòng)將會(huì)掩蓋調(diào)制故障特征,同時(shí)滾動(dòng)體的損傷點(diǎn)在內(nèi)、外圈并不是一直保持接觸,因此滾動(dòng)體故障的振動(dòng)信號(hào)表現(xiàn)出更復(fù)雜的成分,故障特征表現(xiàn)得更加微弱,這也導(dǎo)致了目前諸多故障診斷方法都無法有效提取出深溝球軸承滾動(dòng)體的故障特征頻率[14],為了凸顯本文方法的有效性,將對滾動(dòng)體故障信號(hào)進(jìn)行分析.

        圖16為滾動(dòng)體故障振動(dòng)信號(hào)的時(shí)頻譜,對滾動(dòng)體故障信號(hào)進(jìn)行小波分解、EMD分解及VMD分解,為圖文簡潔,在這里不再展示信號(hào)分解的分量圖.采用db10小波對內(nèi)圈故障信號(hào)做4層分解,通過峭度準(zhǔn)則選取峭度值最大的d4分量進(jìn)行Teager能量算子解調(diào),得到如圖17小波分量的能量算子解調(diào)包絡(luò).

        圖16 滾動(dòng)體故障的時(shí)頻圖

        圖17 小波分量能量算子解調(diào)包絡(luò)

        從圖中可看出包絡(luò)中包含很多無關(guān)的頻率成分凸顯出來,完全無法提取出故障特征頻率,同樣采用EMD對滾動(dòng)體故障信號(hào)做分解得到10個(gè)分量和1個(gè)殘余分量,選取峭度系數(shù)最大IMF4分量進(jìn)行Teager能量算子解調(diào),得到如圖18的包絡(luò)圖,包絡(luò)中依然包含很多無關(guān)的頻率,故障特征頻率無法提?。?/p>

        圖18 EMD分量能量算子解調(diào)包絡(luò)

        最后采用本文的VMD方法能對滾動(dòng)體故障信號(hào)進(jìn)行分析,根據(jù)同上文的VMD算法提取信號(hào)特征最優(yōu)原則,選擇k=7進(jìn)行7層分解,選取峭度系數(shù)最大的IMF6分量進(jìn)行Teager能量算子解調(diào),得到如圖19所示的包絡(luò)圖,可以準(zhǔn)確分解出滾動(dòng)體的故障頻率120 Hz,準(zhǔn)確識(shí)別滾動(dòng)體故障特征頻率,但是小波、EMD以及VMD都沒有分解出滾動(dòng)體故障信號(hào)的基頻,存在一定的不足.但是VMD可以準(zhǔn)確分解出滾動(dòng)體的故障頻率,而使用小波和EMD分解的故障頻率都被淹沒在噪聲里,沒有有效地提取出滾動(dòng)體故障特征頻率,由此可見,VMD相比EMD和小波來說,可分解出更復(fù)雜更微弱的故障信號(hào),對故障特性頻率的提取有更大的優(yōu)勢.

        圖19 VMD分量能量算子解調(diào)包絡(luò)

        4 結(jié) 論

        滾動(dòng)軸承早期故障的振動(dòng)信號(hào)受到強(qiáng)噪聲干擾以及信號(hào)不規(guī)則性的影響,這些因素導(dǎo)致了其故障特征信息難以提取,針對此問題本文提出了一種基于VMD和Teager能量算子結(jié)合的故障診斷方法,并應(yīng)用于滾動(dòng)軸承的故障特征提取中,通過實(shí)驗(yàn)數(shù)據(jù)得出此方法的有效性及以下結(jié)論:

        1)仿真信號(hào)分析結(jié)果表明,本方法能夠在復(fù)雜的信號(hào)里有更有效的分解能力,且具備更強(qiáng)的抑制噪聲干擾和抗模態(tài)混疊的能力.

        2)本文提出的基于VMD的診斷方法對滾動(dòng)軸承早期故障信號(hào)中的微弱特征信息能有效提取,準(zhǔn)確判別了故障特征類型,相比傳統(tǒng)的直接包絡(luò)譜分析方法、EMD、小波等方法,分解效果更好,故障特征信息提取更明顯.

        [1] 鄧飛躍.滾動(dòng)軸承故障特征提取與診斷方法研究[D].北京:華北電力大學(xué),2016.

        [2] Huang N E,Shen Z,Long S R,et al.The Empirical Mode Decomposition and the Hilbert Spectrum for Non-linear and Non-Stationary Time Series Analysis[J].Proceedings of the Royal Society A Mathematical Physi-cal&Engineering Sciences,1998,554(1971):903-995.

        [3] 杜秋華,楊曙年.經(jīng)驗(yàn)?zāi)J椒纸獾母倪M(jìn)及其對球軸承缺陷的診斷[J].振動(dòng)、測試與診斷,2007,27(1):67-70.

        [4] 李 輝,鄭海起,楊紹普.基于EMD和Teager能量算子的軸承故障診斷研究[J]. 振動(dòng)與沖擊,2008,27(10):15-17,22,188.

        [5] 高 強(qiáng),杜小山,范 虹,等.滾動(dòng)軸承故障的EMD診斷方法研究[J]. 振動(dòng)工程學(xué)報(bào),2007,20(1):15-18.

        [6] 蘇玉香,劉志剛,李科亮,等.一種改善EMD端點(diǎn)效應(yīng)的新方法及其在諧波分析中的應(yīng)用[J].電工電能新技術(shù),2008,27(2):33-37.

        [7] 葉 昊,王桂增,方崇智.小波變換在故障診斷中的應(yīng)用[J].自動(dòng)化學(xué)報(bào),1997,23(6):736-740.

        [8] Dragomiretskiy K,osso D.Variational Mode Decomposition[J].IEEE Tran on Signal Processing,2014,2(3):531-544.

        [9] 唐貴基,王曉龍.變分模態(tài)分解方法及其在滾動(dòng)軸承早期故障診斷中的應(yīng)用[J].振動(dòng)工程學(xué)報(bào),2016,29(4):638-648.

        [10] 蒲子璽,殷 紅,張 楠,等. 基于峭度準(zhǔn)則VMD及平穩(wěn)小波的軸承故障診斷[J]. 機(jī)械設(shè)計(jì)與研究,2017,33(1):67-71.

        [11] 馬增強(qiáng),李亞超,劉 政,等. 基于變分模態(tài)分解和Teager能量算子的滾動(dòng)軸承故障特征提取[J]. 振動(dòng)與沖擊,2016,35(13):134-139.

        [12] 王曉龍,唐貴基.基于變分模態(tài)分解和1.5維譜的軸承早期故障診斷方法[J].電力自動(dòng)化設(shè)備,2016,36(7):125-130.

        [13] 鄭小霞,周國旺,任浩翰,等.基于變分模態(tài)分解的風(fēng)機(jī)滾動(dòng)軸承早期故障診斷[J].軸承,2016,41(7):48-53.

        [14] 劉尚坤,唐貴基.改進(jìn)的VMD方法及其在轉(zhuǎn)子故障診斷中的應(yīng)用[J]. 動(dòng)力工程學(xué)報(bào),2016,36(6):448-453.

        猜你喜歡
        故障診斷模態(tài)故障
        故障一點(diǎn)通
        奔馳R320車ABS、ESP故障燈異常點(diǎn)亮
        國內(nèi)多模態(tài)教學(xué)研究回顧與展望
        因果圖定性分析法及其在故障診斷中的應(yīng)用
        故障一點(diǎn)通
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
        江淮車故障3例
        由單個(gè)模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
        基于LCD和排列熵的滾動(dòng)軸承故障診斷
        基于WPD-HHT的滾動(dòng)軸承故障診斷
        内射中出无码护士在线| 亚洲国产精品国自产拍性色| 亚洲国产精品久久精品| 免费看黄色电影| av无码天一区二区一三区| 国产精品成人黄色大片| 免费人成视频网站在线不卡| 国产精品自在拍在线拍| 久久国产精品二区99| 亚洲精品熟女乱色一区| 亚洲av区,一区二区三区色婷婷 | 亚洲国产精品久久久久久无码| 无码av免费精品一区二区三区| 毛片无遮挡高清免费久久| 蜜桃人妻午夜精品一区二区三区| 尹人香蕉久久99天天拍| 国产色秀视频在线播放| 日本高清色惰www在线视频| 亚洲中文字幕在线第六区| 亚洲av日韩av女同同性| 一道久在线无码加勒比| 狠狠亚洲超碰狼人久久老人 | 亚洲精品一区二区三区麻豆| 一本无码人妻在中文字幕免费| 国产精品三级在线专区1| 中文字幕人乱码中文字幕乱码在线| 精品乱人伦一区二区三区| 亚洲狠狠婷婷综合久久| 97久久成人国产精品免费| 一区二区三区国产黄色| 国产亚洲2021成人乱码| 免费人成在线观看播放国产| 日本高清一区二区三区不卡| 亚洲av色欲色欲www| aaaaa级少妇高潮大片免费看| 亚洲a∨好看av高清在线观看| 日本视频一区二区三区在线观看 | 拍摄av现场失控高潮数次| 精品黄色av一区二区三区| 精品一区中文字幕在线观看| 人妻夜夜爽天天爽一区 |