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

基于故障可診斷性的齒輪箱傳感器優化布置

2021-02-26 10:39:50彭珍瑞
振動與沖擊 2021年4期
關鍵詞:模態故障評價

彭珍瑞, 劉 臻

(蘭州交通大學 機電工程學院,蘭州 730070)

齒輪箱作為調節轉速和傳遞扭矩的旋轉機械設備,其是否正常運行將直接決定機器能否正常工作。因此,有必要通過傳感器來測得數據對齒輪箱進行狀態監測與故障診斷,以便排除故障隱患,傳感器優化布置就變得尤為重要。另外,若要達到對齒輪箱故障診斷的目的,前提是結構具有故障可診斷性。因此,尋找一組傳感器集合,使其檢測的數據能夠達到期望的故障最大可診斷性就變得很有必要[1]。

近年來,國內外學者深入研究了傳感器優化布置問題,但對于齒輪箱這類復雜機械設備傳感器優化布置方面研究甚少。Kammer[2]提出有效獨立法,以Fisher信息矩陣的行列式為優化目標,從全部自由度中逐步迭代刪除對目標模態線性獨立性貢獻較小的自由度,保留對目標模態線性獨立性貢獻較大的自由度,直到得到較優的傳感器布置方案。模態動能法的主要思想是在模態動能較大的位置處布置傳感器,以提高結構模態測試時的信噪比,便于信號采集[3]。基于以上兩種方法,很多人進行了改進,劉偉等[4]提出了有效獨立平均加速度幅值法和有效獨立模態動能法兩種傳感器優化布置方法,在考慮模態線性獨立的同時能具有較高平均動態響應和模態動能。詹杰子等[5]提出了有效獨立-改進模態應變能法,最終布置的傳感器位置在具有模態獨立性前提下,可以提高測點組合的抗噪性。還有人將智能算法引入到傳感器優化布置中,以模態動能、模態信息矩陣的行列式、奇異值比值等作為優化目標,實現測點的優化[6-9]。文獻[10]中僅安裝了一些自身感興趣變量的傳感器,對于如何布置傳感器具體安裝數量及位置,以保證實現齒輪箱相關故障診斷方面沒有針對性的研究。王桂蘭等[11]從風機齒輪箱的動力學模型出發,將結構分析方法應用到齒輪箱傳感器的優化配置中,實現傳感器數量最少,識別、隔離可能出現的故障能力最大的配置目標。

以上傳感器優化布置研究大多都集中于齒輪箱的模態觀測性,而對齒輪箱可能發生的故障信號是否能被識別與分離研究較少,因此本文以齒輪箱有限元分析為基礎,從實時測得齒輪箱運行狀態信號出發,能更好的量化評價每一種故障的可檢測性及與其他故障的可分離性。將故障可診斷性應用于齒輪箱傳感器優化布置中,并將模態矩陣的奇異值比值、故障可診斷性、平均加速度幅值構成綜合評價指標,對傳感器數量和位置均進行優化,確定最合理的傳感器布置方案,使得最終布置的傳感器既反映齒輪箱的模態信息,還能對可能出現的故障實現識別與分離。通過ZDH10型齒輪箱故障診斷實驗臺實驗,利用相關系數對實驗結果進行驗證,證明了所提方法的可行性。

1 基本理論

1.1 計算模型

對于一個具有任意自由度的線性系統,其運動方程為

(1)

式中:M,C,K分別為質量矩陣、阻尼矩陣、剛度矩陣;f(t)為力向量;x(t)為響應向量。根據模態疊加原理,響應向量可以表示為

x(t)=Φq(t)

(2)

式中:Φ為模態振型矩陣;q(t)為模態坐標向量。

1.2 有效獨立平均加速度幅值法

有效獨立平均加速度幅值法(effective independent-average acceleration amplitude,EI-AAA)是針對有效獨立法(effective independence,EI)所選取的測點沒有較大的能量,從驅動點頻響函數的近似表達式出發,推導出的一種傳感器優化布置方法。為了使齒輪箱的初選測點集合不但具有盡可能高的平均加速度幅值,而且還要盡可能的保證測試模態具有較好的線性獨立性。利用有效獨立-平均加速度幅值法,即利用式(3)進行迭代,不斷刪除對目標模態獨立性貢獻較小且具有較低加速度響應的節點,即刪除Hii接近于0的節點,保留對目標模態的線性獨立性貢獻較大且具有較大加速度響應的節點,即保留Hii接近于1的節點,直到達到所要求的傳感器數目為止。

H=diag(Φ(ΦTΦ)-1Φ)·diag(ΦΦT)

(3)

1.3 K均值聚類算法

K均值聚類算法是一種迭代型聚類算法。本文中隨著對齒輪箱有限元網格劃分的細化,使得某些自由度在重要模態中的振型值相近,這些自由度在動力載荷下的響應也近似相等,這就會導致信息冗余。因此,采用歐氏距離式(4)作為相似性指標,將齒輪箱s個節點d階振型的數據歸為K類,進而可以有效地避免信息冗余問題[12]。

(4)

式中:K為聚類中心個數;s為節點數;Φi為第i個節點對應的振型;φk為第k個聚類中心振型。

1.4 核密度估計

核密度估計[13]是一種用于估計概率密度函數的非參數方法,本文用核密度對齒輪箱狀態的頻域信號進行密度函數估計。對于一組頻譜數據{c1,c2,c3,…,cn/2},其服從密度函數f(c)的分布,則f(c)的核密度估計為

(5)

式中,w為帶寬;Kw(u)為核函數。本文選應用最廣泛的高斯核函數為核函數

(6)

那么,f(c)的概率密度函數就可以寫為

(7)

1.5 相關函數信息融合

相關函數信息融合是利用實測的兩振動信號的相關函數來確定權值[14]。當傳感器布置數量在兩個及兩個以上時,為了獲得更全面的信息,需要將兩個傳感器采集的信息融合為一個。主要步驟如下:

步驟1利用傳感器測得m個齒輪箱的運行狀態信息X1(n),X2(n),…,Xm(n)。

步驟2利用互相關函數表達式式(8)對測得的每個信號分別與其余信號進行相關,相關程度越大,其對應的權值就越大。

步驟3采用相關信號的能量式(9)、式(10)來表征相關程度,即能量越大,相關程度越大。

步驟4通過式(11)確定各信號權重。由式(12)確定融合信號。

離散性隨機振動信號的互相關函數表達式為

(8)

式中:N為信號長度;RXiXj直接反映兩信號Xi(n)和Xj(n)之間的相關性。則振動信號的能量為

(9)

信號xi(n)與其他信號的總相關能量為

(10)

根據權值與相關函數的能量成正比原則得到

ρ1∶ρ2∶…∶ρm=E1∶E2∶…∶Em
ρ1+ρ2+…+ρm=1

(11)

由信號權重得到融合后的信號

X(n)=ρ1X1(n)+ρ2X2(n)+…+ρmXm(n)

(12)

1.6 故障可診斷性

故障可診斷性包括故障可檢測性和故障可分離性[15-16]。本文中,故障可檢測性FD(fi)指齒輪箱發生某種故障與未發生故障的概率密度函數的最小K-L散度值,故障可分離性FI(fi,fj)指齒輪箱發生某種故障與齒輪箱發生另外一種故障的概率密度函數的最小K-L散度值。K-L散度的計算公式為

(13)

式中,pi(c)和pj(c)為對兩個不同故障的概率分布。故障可檢測性FD(fi)可表示為

FD(fi)=min[K(pi‖pNF)]

(14)

式中:PNF為齒輪箱未發生故障的概率密度函數;pi為齒輪箱發生某種故障的概率密度函數。當FD(fi)越大,表明故障可檢測性越強;當FD(fi)=0時,表示故障不能被檢測。故障可分離性FI(fi,fj)可表示為

FI(fi,fj)=min[K(pi‖pj)]

(15)

式中,pi和pj為齒輪箱發生一種故障和另外一種故障的概率密度函數,當FI(fi,fj)越大時,表明故障fi和fj的可分離性越明顯,當FI(fi,fj)=0時,表明故障fi和fj能被分離。

1.7 相關系數

若齒輪箱處于大致相同的故障狀態,那么相同工況下采集到的振動信號在頻譜上會有相似的密度函數曲線,而不同故障狀態下振動信號頻譜的密度函數曲線則有較大的差異,利用相關系數衡量密度函數曲線之間的相似性或者差異性。假定采集到齒輪箱某兩種狀態頻譜的密度函數為f1(c)和f2(c)(c=1,2,3,…,n/2)則它們之間的相關系數為

(16)

(17)

(18)

相關系數ρ12的取值范圍為[-1,1],其值越大,就表明兩者相似度越高,反之,相似度就越差。

2 評價指標

目前,傳感器的布置方案有多個評價指標,如奇異值比值、平均加速度幅值等,但這些都是反映齒輪箱的模態觀測性,并不能反映齒輪箱中可能出現故障的識別與分離,因此,本文在評價傳感器布置方案時,引入了故障可診斷性。另外,如果僅僅以奇異值比值、平均加速度幅值、故障可診斷性其中一個指標進行評價時,評價目標比較單一,僅僅只能考慮到一個因素,而合理的傳感器布置,應需考慮多個因素。再加上這三者之間并不會發生沖突,因為故障可診斷性是從測點實時測得的信號出發,而奇異值比值和平均加速度幅值是從齒輪箱的振型出發,則故障可診斷性與奇異值比值和平均加速度幅值這兩者之間分別獨立,也就不會存在沖突;另外,平均加速度幅值準則是保證傳感器能夠布置在具有較大響應的節點處,而奇異值比值是盡可能的保證測點振型之間的正交性[17],也就是保證各測點振型的獨立性,并不會影響測點的加速度響應,則奇異值比值與平均加速度幅值是相互獨立的,所以三個評價指標是相互獨立的,并不會發生沖突。因此,本文將奇異值比值、故障可診斷性、平均加速度幅值三個評價準則構成綜合指標評價,評價不同布置方案,確定最優傳感器布置方案。

2.1 奇異值比值

模態矩陣的奇異值比(singular value ratio,SVR)作為判斷傳感器布置優劣程度的一項重要指標,為模態矩陣奇異值的最大值與最小值之比

(19)

式中:svd(·)為奇異值向量;Φ為所布置傳感器構成的模態振型矩陣。為統一綜合評價指標最大化,取S的倒數,具體函數轉換為

(20)

f1越大則傳感器布置效果越好。

2.2 平均加速度幅值

結構模態測試時一般都采用加速度傳感器,傳感器應該布置在具有較大的平均加速度幅值(average acceleration amplitude,AAA)的位置處,有利于數據的采集和提高測量的抗噪能力,計算公式為

Y=diag(ΦΦT)

(21)

式中,Y為一列向量,當所有元素都較大的時候,測點位置會較好。因此,建立平均加速度函數為

f2=mean(Y)

(22)

式中,mean(Y)為所布置位置的平均加速度幅值的整體平均值。

2.3 故障可診斷性

故障可診斷性是對齒輪箱可能發生故障可被檢測和與其余故障可分離的綜合體現,當每一種故障的可檢測性和可分離性都較大的時候,測點的布置位置會較好,因此將故障可診斷性的平均值f3作為目標函數。

f3=mean(FD(fi)+FI(fi,f2)

(23)

2.4 綜合評價指標

為避免某個指標信息被淹沒的情況,對各個指標進行歸一化處理并求和,即可得綜合評價指標函數的表達式

(24)

式中,f為傳感器布置結果的綜合體現,如果綜合評價指標越大,說明各項評價指標都較好,傳感器布置更合理。

3 傳感器優化布置

由于目前大多數的研究都集中于齒輪箱的模態觀測性,而當齒輪箱發生某種故障時,故障是否能更容易檢測到并將它們區分開,研究較少,也就是說當傳感器布置在哪個位置處(如輸入端、輸出端等)時,既能夠檢測到這種故障,又能將這種故障與其余故障更好的辨識出來。因此,為了使安裝的傳感器獲取的信息不但可以反映齒輪箱的模態信息,還能對可能出現的故障實現識別與隔離,構建了一種基于故障可診斷性的齒輪箱傳感器優化布置方法,其流程如圖1所示。

圖1 傳感器優化布置流程圖Fig.1 Flow chart of optimal sensor placement

3.1 初選測點

先建立齒輪箱有限元模型,進行模態分析,提取齒輪箱的模態振型。由于有限元網格劃分精細,大部分自由度有相似的動力特性,會引起信息冗余。本文根據自由度在重要模態中振型值的相似性,使用K均值聚類算法按隸屬度對有限元劃分的節點進行自動分類,將動力特性相似的自由度劃分為一簇。在每一簇中再利用有效獨立-平均加速度幅值法選出精英個體(即所含振動信息豐富的自由度),再將這些精英個體組為一體,再從這些精英個體中選出具有較大平均加速度幅值和較大模態獨立性的節點作為初選測點。對齒輪箱節點聚類步驟如下:

步驟1從所提取的齒輪箱振型矩陣中隨機選取K個節點所對應的振型作為最初的聚類中心。

步驟2根據每個節點所對應的振型與這些聚類中心的歐氏距離,將節點分到距離最近的聚類中心所對應的類中。

步驟3更新聚類中心,即求得每一類中所有元素的均值,并將其作為新的聚類中心。

步驟4判斷聚類中心和目標函數的值是否發生改變,若不變,則輸出結果,若改變,則返回步驟2。

3.2 故障可診斷性判斷

首先,由3.1節初選的節點找到實體中齒輪箱所對應的位置。并在所有的初選位置上分別布置傳感器,測得齒輪箱正常的以及可能發生故障(如主動輪斷齒、從動輪磨損等)的狀態信號,其次,對測得的信號進行快速傅里葉變換得到其對應的頻譜圖,接著對其頻譜圖進行核密度估計式(7)求得概率密度函數。最后,用式(14)確定其對應故障的可檢測性;用式(15)確定故障的可分離性,用式(23)確定在每一個節點處的故障可診斷性。

由動態規劃算法的思想,選擇故障可診斷性最大的節點位置,然后確定此節點位置,再從剩余的節點中分別與其兩兩組合測得可能發生故障的信號,然后將這兩個位置的信號利用相關函數數據融合算法[13]進行信息融合,再利用核密度估計對其分別進行概率密度函數估計,確定它們的故障可診斷性,再選擇故障可診斷性最大的兩個節點位置組合,然后確定這兩個節點位置,依次確定第三個、第四個,直至所有初選測點進行信息融合并確定它們的故障可診斷性為止。

假設初選測點是A,B,C,D四個位置,有a(正常),b,c,d四種狀態。具體的故障可診斷性步驟:

步驟1分別在A,B,C,D這四個位置處,測得a,b,c,d四種狀態的信號。

步驟2對在每一個位置測得的信號進行快速傅里葉變換得到其對應的頻譜圖,利用核密度估計求得它們的概率密度函數。

步驟3求b,c,d的概率密度函數與正常信號(a)概率密度函數之間的K-L散度值確定其對應故障的可檢測性;再分別求b,c,d兩兩之間概率密度函數之間的K-L散度值,確定故障的可分離性,進而確定在每個節點處的故障可診斷性。

步驟4選擇故障可診斷性最大的位置,確定此位置并分別與其余位置組合,判斷各組的故障可診斷性。假如A最大,就確定A位置,繼續判斷A與B,A與C,A與D的故障可診斷性。當確定的位置在兩個及兩個以上時,對測得的信號進行信息融合。

步驟5返回步驟2,直到初選測點組合完成為止,也就是A,B,C,D四個位置組合完為止。

3.3 綜合指標確定布置方案

當需要確定最合理的傳感器數目時,應使用盡量少的傳感器來使得各目標函數值均較好。合理的傳感器布置應該綜合考慮各項指標,因此,以奇異值比值、平均加速度幅值和故障可診斷性構成的綜合評價指標式(24)對由3.2節確定的節點位置進行評價,最終確定傳感器最優布置位置。

4 應用實例

4.1 齒輪箱主要參數

齒輪箱試驗臺主要由伺服電機、ZDH10型齒輪箱、磁粉制動器組成,如圖2所示。以ZDH10型齒輪箱為研究對象,其具體參數:材料為HT150,其泊松比為0.25,密度為7 100 kg/m3,彈性模量為1.5×1011Pa。在齒輪箱的各零件中,齒輪故障所占的比例最大,能達到60%以上。因此本文就齒輪可能發生的故障進行故障可診斷性判斷,齒輪具體參數如表1所示。齒輪狀態如表2所示。可能發生的故障如圖3所示。選取INV982X系列的壓電式加速度傳感器。

圖2 故障診斷裝置Fig.2 Fault diagnosis setup

圖3 齒輪故障Fig.3 Gear fault

表1 齒輪參數Tab.1 Gear parameters

表2 齒輪狀態Tab.2 Gear state

4.2 齒輪箱傳感器位置初選測點

依照圖紙在Solidworks中建立ZDH10型齒輪箱的實體模型,然后導入ANAYS Workbench 17.0進行網格劃分,計算前6階模態,提取振型矩陣。將齒輪箱箱體劃分為10 207個節點,5 143個單元。首先根據各節點振型的動力相似性,利用K均值聚類將節點聚為4類,對每一類用EI-AAA法篩選出100個節點。再將每一類的100個節點集結在一起,用EI-AAA法從這400個節點選擇5個節點,作為初選測點。具體初選的測點位置如表3所示。節點對應測點位置如圖4所示。對應在故障實驗臺上的具體位置如圖2所示。

表3 初選測點位置Tab.3 Initial measuring points

圖4 初選測點位置Fig.4 Initial measuring points

4.3 故障可診斷性判斷

在實驗中,驅動電機的轉速為348 r/min,磁粉制動器的施加電流為0.1 A,采樣時間為6 s,采樣頻率5 kHZ,振動傳感器測量齒輪箱垂直方向上的振動數據,當傳感器布置在節點7460、節點6579、節點7788時,齒輪箱表面都比較平整,可以直接將傳感器貼在上方。在節點6925(輸入端)和節點9121(輸出端)處布置傳感器時,會存在傳感器在節點的左側和右側(背對伺服電機)兩種偏移情況,通過試驗發現在這兩個位置處,測得的信號基本一樣,因此就將傳感器布置在較為平整的一邊,輸入端在左側較平整,輸出端在右側較平整。

4.3.1 確定第一個位置

分別在節點7460、節點6959、節點6925、節點7788和節點9121處,測得F1,F2,F3,F4和F5五種狀態的信號,進行核密度估計,以節點7460為例,其對應的頻譜及其核密度估計圖如圖5所示,然后利用K-L散度確定F2,F3,F4,F5四種故障的可檢測性與可分離性,具體數值如表4所示。

表4 故障可診斷性量化評價Tab.4 Quantitative evaluation of fault diagnosability

圖5 五種狀態下的頻譜圖及其對應的密度函數Fig.5 Spectrum and their corresponding density functions under five states

當在節點6925處時,故障可診斷性的均值最大,確定第一個位置為:節點6925。

4.3.2 確定第二個位置

在確定第一個位置的基礎上,繼續判斷節點6925+節點7460、節點6925+節點9121、節點6925+節點7788、節點6925+節點6579的故障可診斷性。當判斷以上兩個節點位置的故障可診斷性時,分別測得各組節點F1,F2,F3,F4和F5等五種狀態的信號,并進行對應狀態的信息融合,再對融合信號的頻譜圖進行核密度估計,進而確定四種故障的可診斷性,具體數值如表5所示。

表5 兩信號融合下的故障可診斷性量化評價Tab.5 Quantitative evaluation of fault diagnosability based on two-signal fusion

當傳感器布置在節點6 925和7 460處時,故障可檢測性和故障可分離性的均值最大,因此確定前兩個位置即:節點6925+節點7460。

4.3.3 依次確定第三、第四、第五位置

在前兩個位置的基礎上,然后繼續判斷節點6925+節點7460+節點9121、節點6925+節點7460+節點7788、節點6925+節點7460+節點6579。由于篇幅限制,就沒列舉具體數值。當在節點7788+節點7460+節點6925處時,故障可檢測性和故障可分離性的均值較大,因此,確定前三個位置即:節點6925+節點7460+節點6579,接下來繼續判斷節點6925+節點7460+節點6579+節點9121、節點6925+節點7460+節點6579+節點7788。當在節點6925+節點7460+節點6579+節點7788處時,故障可檢測性和故障可分離性的均值較大,確定傳感器前四個位置即節點6925+節點7460+節點6579+節點7788;然后確定節點6925+節點7460+節點6579+節點7788+節9121的故障可診斷性。具體故障可診斷性變化圖如圖6(b)所示。

4.4 綜合評價確定方案

根據故障可診斷性得到布置不同數量的傳感器對應的不同布置結果,首先根據式(19),式(22)計算出不同布置方案對應的奇異值比值、平均加速度幅值。利用式(24)對不同的傳感器布置方案進行綜合評價,確定最優傳感器布置方案。具體評價結果如表6所示,則選擇兩個傳感器即節點6925和節點7460時齒輪箱的故障可診斷性和模態可觀測性最大。

表6 綜合評價指標Tab.6 Comprehensive evaluation index

單個評價指標結果如圖6所示,如果僅僅以奇異值比值(判斷兩測點振型的正交性,測點至少有兩個,因此從2開始)、平均加速度幅值、故障可診斷性其中一個進行評價時,奇異值比值是在布置兩個傳感器(即節點6925和7460)時,效果較好,如果存在沖突,平均加速度幅值和故障可診斷性在兩個(即節點6925和7460)測點時,效果會最差。而在兩個測點時,平均加速度幅值和故障可診斷性的值都較大,即效果較好,則這個評價指標之間不存在沖突。因此,利用綜合評價指標能更好地量化評價傳感器位置。由圖6也可知,當傳感器布置在兩個節點時,奇異值比值和故障可診斷性都最好,平均加速度幅值也僅次于最好值,也印證了綜合評價指標的合理性。

圖6 三個評價指標Fig.6 Three evaluation indexes

4.5 方法對比

當利用有效獨立平均加速度幅值法直接對齒輪箱進行傳感器優化布置時,得到的節點位置為節點9121和節點7788,故障可診斷性均值為0.187 22,遠小于在節點7460和節點6925處的值。具體故障可診斷性量化評價如表7所示。

表7 EI-AAA法故障可診斷性量化評價Tab.7 Quantitative evaluation of fault diagnosability by EI-AAA

4.6 結果驗證

為了證明在節點6925和節點7460兩個位置處,所布置傳感器的合理性,分別測得齒輪箱F1,F2,F3,F4,F5的狀態信號,進行核密度估計,并求它們在這兩個位置處兩兩之間的概率密度函數的相關系數。當他們之間的相關系數越小,即相關度越低,說明兩種故障越容易區分。相反,這兩種故障越難區分。為了證明在節點6925和節點7460比其他位置更為合理,與其他三種布置方案進行比較,相關系數評價如表8所示,當在節點6925和節點7460兩個位置處時,他們之間的相關系數最小,任意兩故障也更容易區分。從而也證明了傳感器在節點6925和節點7460兩個位置處布置的合理性。

表8 相關系數評價Tab.8 Correlation coefficient evaluation

5 結 論

本文構建了一種齒輪箱故障可診斷性的傳感器優化布置方法,以ZDH10型齒輪箱為例進行傳感器優化布置,獲得了較好的布置效果。

(1) 利用K均值聚類算法和有效獨立平均加速度幅值法從各聚類自由度中選擇出所含信息較豐富的自由度作為待選測點,能夠有效避免信息冗余。

(2) 將多個傳感器信息進行信息融合,并利用核密度估計對其對應的頻譜進行密度函數估計,再利用K-L散度進行故障可診斷性判斷,確定每一種可能發生的故障的可檢測性和可分離性。

(3) 通過使用由模態矩陣的奇異值比值、故障可診斷性、平均加速度幅值構成的綜合評價指標對傳感器數量和位置進行了優化,綜合考慮了各方面因素,可獲得最合理的傳感器布置方案。

猜你喜歡
模態故障評價
SBR改性瀝青的穩定性評價
石油瀝青(2021年4期)2021-10-14 08:50:44
故障一點通
奔馳R320車ABS、ESP故障燈異常點亮
國內多模態教學研究回顧與展望
故障一點通
基于Moodle的學習評價
基于HHT和Prony算法的電力系統低頻振蕩模態識別
江淮車故障3例
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
保加利亞轉軌20年評價
主站蜘蛛池模板: 久久精品这里只有国产中文精品| 一区二区午夜| 性做久久久久久久免费看| 四虎永久免费在线| 一级毛片免费的| 国产精品所毛片视频| 色婷婷在线影院| 女人av社区男人的天堂| 欧美福利在线观看| 1级黄色毛片| 毛片视频网址| 福利一区在线| 日韩小视频在线播放| 国产亚洲精品资源在线26u| 高清欧美性猛交XXXX黑人猛交| 亚洲日韩欧美在线观看| 欧美日韩高清在线| 91视频日本| 日本影院一区| 国产女人综合久久精品视| 五月婷婷欧美| 国产精品女熟高潮视频| 国产成人av一区二区三区| 国产另类视频| hezyo加勒比一区二区三区| 日韩在线永久免费播放| 国产丝袜无码精品| 国产日本一区二区三区| 国产Av无码精品色午夜| 精品国产一区91在线| 不卡无码网| 五月激情婷婷综合| 好紧太爽了视频免费无码| 欧美亚洲日韩中文| 在线欧美a| 永久免费AⅤ无码网站在线观看| 欧洲亚洲欧美国产日本高清| 四虎AV麻豆| 国产成人综合网| 国模视频一区二区| 亚洲天堂久久久| 91精品国产91欠久久久久| 四虎精品黑人视频| 国产情侣一区| 五月天综合网亚洲综合天堂网| 国产中文一区二区苍井空| 久久青草热| 亚洲精品午夜天堂网页| 色婷婷亚洲综合五月| 丁香婷婷在线视频| 一级一级一片免费| 超碰精品无码一区二区| 波多野结衣无码中文字幕在线观看一区二区| 无码在线激情片| 婷婷伊人五月| 色男人的天堂久久综合| 国产凹凸一区在线观看视频| 亚洲午夜福利在线| 青草视频免费在线观看| 亚洲天堂区| 日本日韩欧美| 国产Av无码精品色午夜| 免费观看成人久久网免费观看| 久久久亚洲色| 国产欧美视频在线| 欧美视频免费一区二区三区| 欧美在线伊人| 国产精品原创不卡在线| 国产精品视频系列专区 | 一级全黄毛片| 亚洲国产成人精品青青草原| 久久精品只有这里有| 午夜毛片福利| 亚洲一级毛片在线播放| 久久无码av三级| 天天综合网色中文字幕| 丰满人妻被猛烈进入无码| 久久精品欧美一区二区| 免费无遮挡AV| 久久这里只有精品免费| 日本不卡在线播放| 国产又爽又黄无遮挡免费观看 |