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

干氣密封環磨合過程摩擦振動信號混沌特性分析

2021-06-10 01:27:46陳金林丁雪興唐莉萍
工程科學與技術 2021年3期
關鍵詞:關聯振動信號

陳金林,丁雪興,唐莉萍

(蘭州理工大學 石油化工學院,甘肅 蘭州 730050)

近年來,國內外對干氣密封環端面間摩擦磨損研究逐漸重視,干氣密封環的摩擦磨損通常出現在啟停階段,但實際上由于加工制造、裝配誤差和工作環境的影響,在運行階段也會出現一定程度的磨損[1],同時伴隨著溫升、劃傷等現象的產生,這些因素會進行不斷累積,影響密封端面的性能。磨合是干氣密封環在使用初期界面所必須經歷的磨損階段,掌握磨合過程的運行狀態對提高端面的磨合質量有重要意義。摩擦振動信號是摩擦系統的重要輸出參數,在干氣密封環的摩擦磨損過程中是必然存在的,并且蘊含了大量反映系統狀態的信息。摩擦振動信號通過傳感器和數據采集系統獲取,采集過程不影響機器的正常運行,可實現無損檢測。對于干氣密封環摩擦磨損特性的研究主要有理論建模和試驗研究兩個方面。因摩擦表面具有分形特征[2],孫寶財[3]、陳金林[4]等采用分形理論建立了干氣密封環端面摩擦剛度模型。Ding等[5]對干氣密封環DLC薄膜織構表面摩擦學性能進行了試驗研究。Jiang等[6]同時考慮良好的密封性能和可能的優良耐磨性能,進行了一種新的干氣密封環表面結構設計。但利用摩擦振動信號反映密封環端面的磨合過程迄今尚未見報道。

混沌理論能夠良好反映系統的非線性特征,與時頻域分析法相比,有著更好的直觀性,被廣泛用于摩擦學問題的研究[7]。朱華等[7]發現摩擦力信號具有混沌特性。之后,Ionita[8]研究一種發現系統混沌敏感性的方法,將確定性混沌稱為裝配和部件失效過程的一般建模技術,并得出混沌行為出現失穩狀態時系統會發生故障的結論。Takuji等[9]分析了一個具有干摩擦的受迫機械動力系統,該系統能產生混沌粘滑振動。Ding等[10]在環盤式摩擦磨損試驗機上進行了摩擦試驗,對摩擦噪聲進行了混沌分析。結果表明,摩擦噪聲是混沌的。Liu等[11]在球形盤上測試儀上進行了磨合磨損試驗,利用混沌吸引子分析了切向摩擦振動和法向摩擦振動的變化。Sun等[12]在銷–盤磨損試驗機上進行了磨合磨損試驗,發現從吸引子的軌跡狀態中可以分辨出摩擦副從磨合向穩定磨損階段的轉變。不同材料、不同摩擦方式,其摩擦特性往往有所不同,而摩擦振動信號、摩擦噪聲等可能受到機械振動、機械噪聲等的干擾,直接采用混沌理論對信號進行處理分析,難以得到真實摩擦特性規律。

作者利用EEMD方法提取干氣密封環磨合過程中的摩擦振動信號,經該方法降噪處理后的摩擦振動信號能更好地消除因機械振動對信號的干擾,數據更加真實。基于混沌理論構造動力系統的摩擦振動信號吸引子相圖,計算其特征量關聯維數。研究磨合狀態與混沌特征量之間的關系,可為干氣密封環磨合狀態的監測和識別提供理論基礎。

1 試 驗

1.1 試驗設備與摩擦副

試驗采用MMW–1立式萬能摩擦磨損試驗機進行干氣密封環摩擦磨損試驗,試驗裝置如圖1(a)所示。摩擦振動信號的采集采用YSV2303S型三軸加速度傳感器(頻率范圍:1~7 000 Hz,靈敏度:100 mV/g)。其中,加速度傳感器固定在下副盤上方,緊貼下試件,采樣頻率為64 kHz,每0.1 s采集6 400個點。此外,試驗還用到了丙酮清洗機、烘干機。

因實際工作中石墨材料硬度較差,常用的“軟碰硬(C–SiC)”密封副不能滿足工業需求,故試驗選用“硬碰硬(SiC–SiC)”摩擦副且在靜環表面鍍上一層DLC薄膜。試驗的上試件(動環)為SiC螺旋槽密封環,下試件(靜環)為DLC薄膜密封環。制作工序如下:首先,將動環加工至特定形狀尺寸;然后,將試件表面進行打磨拋光,使得其表面粗糙度約為0.212 μm;最后,用激光打標機進行刻槽,槽深10 μm,槽數16,螺旋角16°且周向均勻分布。將靜環用超聲波清洗10 min,烘干后采用磁控濺射法進行DLC鍍膜,鍍膜厚度為3 μm。試件結構如圖1(b)、(c)所示。

圖1 測試裝置及試件Fig. 1 Testing instruments and specimens

1.2 試驗方法與試驗工況

實際生產條件下,密封端面間的接觸比壓為0.1~0.8 MPa,平均線速度范圍為0.6~1.6 m/s[13],因此選擇在載荷150、450 N,轉速500 r/min(即端面比壓0.167、0.503 MPa,平均線速度1.36 m/s)的工況下進行試驗。試驗前對試驗機進行調試,設定運行時間600 s,上下試件用乙醇和丙酮超聲清洗后干燥數分鐘。依次設置轉速、載荷及終止時間,測量振動信號、摩擦扭矩、溫度等參數,摩擦系數可根據摩擦扭矩換算得出。試驗停止后讀取試驗結果并記錄,最后將試件進行清洗烘干。

1.3 摩擦振動信號提取

一般來說,振動信號是較易于獲取的參數,但實際測量的信號包含著整個系統所產生的各種干擾信號,在此基礎上的數據分析無法反映密封面狀態的變化,需要在初始數據之上進行摩擦振動信號的有效提取。

提取摩擦振動特征信號一般采用集合經驗模態分解(EEMD)方法[14]。該方法是基于經驗模態分解(EMD)法的一種改進算法[15]。通過EMD分解后,特征信號在不同分辨率下充分展現出來,得到一系列本征模態分量(IMF),包含了原始振動信號中不同時間尺度的局部特征。但由于各個IMF分量間會存在較大的模式混疊現象。因此,EEMD方法提出在此基礎上將原始信號與特定幅值系數的高斯白噪聲進行疊加運算,之后再不斷地進行EMD分解,獲得若干個模式混疊程度小或者無混疊的IMF分量。摩擦振動信號的提取則是通過選擇包含其特征的IMF分量進行重構特征信號。

EEMD的具體計算步驟為:

1)初始化EMD總體平均次數M,令i=1。

2)對xi(t)進行i次EMD分解計算:

①在原始信號x(t)中加入一個幅值系數為k的高斯白噪聲ni(t),第i次加噪后的信號xi(t):

②用EMD分解,得到一組分量IMFcj,i(j=1,2,···,I)。其中,cj,i為第i次EMD分解后得到的第j個本征模態分量。

③若i

3)對M次分解后的IMF分量進行均值計算:

在EEMD計算過程中,需要對兩個參數,即白噪聲幅值系數k和EMD計算總次數M進行賦值。由文獻[11]可知,k的確定范圍為0.01~0.50,在進行降噪的時候,k值一般取值的范圍為0.01~0.10。加入噪聲后對于結果的影響記作e,存在以下數量關系:

式中,Mav為EMD分解計算總次數的平均值。

從式(3)中可以看出,當k值越小時,e值也越小,說明分解結果精度越高。但是,當k過小時,加入的白噪聲無法引起信號局部極值點發生改變,使得計算失去意義,因此按照規定范圍取值十分重要。另外,M值的取值范圍一般定為100~300[16]。

以轉速為500 r/min、載荷為450 N工況條件下的干氣密封光面環z方向(切向)上的原始振動數據為例,進行分析計算(本文皆取的切向振動信號)。用MATLAB對EEMD算法進行編程,對其中相關參數進行賦值,這里的白噪聲幅值系數、總體平均次數分別取0.1、100。圖2為干氣密封環在600 s摩擦振動信號的EEMD分解結果,通過分解可以得到8個IMF分量和一個殘差R,IMF1~IMF8分量頻率由高到低變化。由于摩擦振動信號頻率高、振幅小[17],所以,選擇前兩個分量IMF1、IMF2合成摩擦振動信號特征信號。圖3為120、480、600 s時間下摩擦振動信號降噪前后的對比。從圖3可以看出:原始數據時域波形雜亂無序,沒有規律可循;提取出的摩擦振動信號幅值減小,呈現更明顯的規律變化[17]。

2 混沌分析理論

混沌理論主要是通過計算某一單變量時間序列來研究整個系統的混沌動力學行為,一般用吸引子和混沌參數進行詳細描述。混沌理論相關的研究主要包括3個方面:1)重構相空間;2)混沌吸引子和混沌參數的計算分析;3)混沌特性識別。

2.1 相空間重構及混沌參數的計算

對一個單變量時間序x1,x2,x3, ···xn,依據Takens嵌入定理[16],在相空間中構造一個矩陣:

式中,m為嵌入維數,τ為延遲時間,n為時間序列點數,N=n–(m–1)τ。

對于一個單變量混沌時間序列{x(i)},可以先求取其自相關函數A(τ),再畫出自相關函數A(τ)圖像和時間τ之間的函數圖像,當函數值開始下降至初始值的1–1/e時,此時的τ即為所要求的最佳延遲時間。

圖2 試驗600 s時刻摩擦振動信號EEMD分解結果Fig. 2 Results of EEMD decomposition of friction vibration signal at 600 s

圖3 振動信號時域波形圖Fig. 3 Time-domain waveform of vibration signal

當確定嵌入維數時,若m的取值過小,那么空間相點將會發生重疊,造成一部分信息的丟失;若取值太大,則會延長計算時間,所以應選擇最佳嵌入維數。由于C–C法算法較為復雜,過程繁瑣,程序運行較慢,所以用飽和關聯維數法進行計算。這種方法主要是算出各個嵌入維數下的關聯維數、最佳嵌入維數,這種方法也叫做飽和關聯維數法。

關聯維數在一定程度上反映了相空間中系統的運動狀態,若關聯維數大,說明系統相關聯的點數越多,相空間中狀態點越密集,混沌吸引子收斂程度高。關聯維數也象征著系統的有序度,通過對關聯維數的計算可以判斷混沌信號在小尺度下的精細復雜程度,從而表征吸引子的形態。

關聯維數一般采用G–P算法進行計算[18]。在式(4)中構造的矩陣當中任意兩矢量間的距離為:

對于任意給出的一個正數ε滿足rij≤ε的所有矢量,稱為關聯矢量,將其數目記為N1(ε);將rij>ε的數目記為N2(ε),則總距離數目為N(ε)=N1(ε)+N2(ε)。把距離不大于ε的點對在所有點對中所占比例記為C(ε),即C(ε)=N1(ε)/N(ε),C(ε)稱為密度相關函數。C(ε)有如下形式:

則關聯維數的定義為:

式中,θ為Heaviside單位函數,當x≤0時取值為0,當x>0時取值為1。

由上述計算可以看出來,C(ε)在ε→0上與變量ε存在關系limC(ε)→εD。D為關聯維數,若ε取值合適,則可以完整刻畫出吸引子的自相似結構。在計算時,可以在lgC(ε)–lg ε的圖中選擇線性較好的區域進行擬合,直線的斜率便為關聯維數D。

2.2 空間吸引子演化

通過以上參數的計算可以在m維相空間中完整刻畫出吸引子演化軌跡,但由于人類對于高維視覺的局限性,通常采用主成分分析法將m維吸引子圖形降維至3維空間中,使其容易進行觀測分析,具體方法如下。

在空間中選取3個主矢量來投影吸引子,計算m階矩陣Y:

計算得到Y的特征值為λ1,λ2,λ3,…,λm(λ1≥λ2≥λ3≥,…,≥λm≥0)。其中,λ1、λ2和λ3為矩陣Y的主特征值,將重構的矩陣X降維至3個主矢量方向,得到矩陣為:

式中,ξ1、ξ2、ξ3為矩陣的3個主矢量。矩陣α的每一行均對應一個坐標,將這些點線連接起來即可得到吸引子軌跡。

2.3 混沌特性的判別

2.3.1 主分量分析法

主分量分析法又被稱為PCA分布法[19–20]。根據相關資料[21]可以得到,混沌信號和噪聲的譜圖主分量分布圖有著顯著差異,噪聲的譜圖是一條與橫坐標相平行的直線,而混沌信號的主分量譜圖為一條斜率為負的直線。1963年,美國氣象學者洛倫茲在研究天氣預報的規律中發現在某一確定性的方程中會存在混沌行為,便根據大氣對流模型發現了洛倫茲吸引子[22]。圖4為洛倫茲吸引子y軸分量的主分量譜圖。從圖4中可以看出,洛倫茲的主分量譜圖過定點且斜率為負,存在混沌特征。因此,可以用該方法識別某一時間序列是否有混沌特征。

圖4 洛倫茲方程y方向的主分量譜圖Fig. 4 Principal component spectrum of the y direction of the lorentz equation

2.3.2 最大Lyapunov指數判別法

最大Lyapunov指數表現了系統在相空間中相鄰軌道間收縮或發散的平均指數率,是一種整體特征,值可正、可負、也可為0。若系統的最大Lyapunov指數大于0,那么該系統一定具有混沌特征,可將此作為分析時間序列是否混沌的一個判斷依據。

3 試驗分析及討論

3.1 摩擦振動信號相空間重構

對于提取出的摩擦振動信號分別利用第2節所介紹的自相關函數法和飽和關聯維數法計算相空間重構參數延遲時間和最佳嵌入維數。在轉速為500 r/min、載荷為150 N的工況參數下,以摩擦進行到第480 s時的試驗數據為例,利用MATLAB軟件編寫算法求得其自相關函數和嵌入維數。自相關函數的圖像結果如圖5(a)所示。由圖5(a)可知,當延遲時間的值為2時,自相關函數值開始小于初始值的1–1/e倍,因此可以得出,當延遲時間取為2時,為最佳延遲時間。

圖5(b)為摩擦振動信號的雙對數變化曲線,無標度區間為–5.75~–4.75,將嵌入維數從2增加至30,找到其最佳擬合曲線,得到的關聯維數結果如圖5(c)所示。當嵌入維數達到24時,關聯維數在可接受誤差范圍內變化(一般取相對誤差rD=((D2(m2)–D2(m1))/D2(m1))≤5%),因此,m=24為最佳嵌入維數,可以認為m=24時的關聯維數即為該時間序列的關聯維數值。按照該方法計算得到的摩擦試驗過程中各時間序列的最佳嵌入維數和延遲時間列于表1、2中。

圖5 延遲時間和嵌入維數Fig. 5 Time delay and optimal embedding dimension

表1 500 r/min、150 N摩擦振動信號最佳嵌入維數和延遲時間Tab. 1 Optimal embedding dimension and delay time for frictional vibration signals under 500 r/min and 150 N

表2 500 r/min、450 N摩擦振動信號最佳嵌入維數和延遲時間Tab. 2 Optimal embedding dimension and delay time for frictional vibration signals under 500 r/min and 450 N

3.2 摩擦振動信號混沌特性識別

在轉速為500 r/min,載荷為150、450 N的工況條件下,對摩擦振動信號時間序列進行主成分分析。圖6(a)為在480 s時,兩個工況下摩擦振動時間序列的主分量譜圖。從圖6(a)可以看出,兩個工況下摩擦振動信號的主成分譜圖都可近似看成一段斜率為負的直線,這說明在干氣密封環的摩擦過程中摩擦振動信號存在混沌特征。

圖6(b)為整個磨合過程的摩擦振動信號最大Lyapunov指數數值。從圖6(b)可以看出,在整個磨合試驗過程當中,各個工況下摩擦振動信號的最大Lyapunov指數皆大于0,這說明隨著磨合過程的進行,摩擦振動時間序列始終保持著混沌特征,可以利用混沌理論研究。

圖6 主分量譜圖和最大Lyapunov指數Fig. 6 Principal component and maximum Lyapunov index spectrum

3.3 摩擦振動混沌吸引子演變規律

在摩擦副的運動過程中,摩擦系數的變化規律可以很好地反映摩擦副的摩擦狀態。首先,分析摩擦系數的變化趨勢;然后,用摩擦系數曲線的變化為參照對象分析摩擦振動信號的混沌特征。

圖7為干氣密封環在不同工況下的摩擦系數隨摩擦時間的變化情況。

圖7 不同工況下摩擦系數隨時間的變化Fig. 7 Variation of friction coefficients with time under different working conditions

由圖7可知,隨著時間的增加,摩擦系數皆呈現減小的趨勢,隨后保持穩定。當轉速為500 r/min、載荷為150 N時,摩擦系數在400 s左右附近達到平穩,數值為0.10;當轉速為500 r/min、載荷為450 N時,摩擦系數在260 s保持平穩,值保持在0.055微幅波動。這是因為在摩擦的初期,類金剛石薄膜表面在局部接觸點處由于高應力的作用,發生黏著磨損,進而引起摩擦系數值較大;在之后的摩擦過程中,由于黏著效應產生一定量磨屑,這些磨屑在摩擦過程中逐漸向類金剛石薄膜表面轉移,形成一種轉移膜,該轉移膜具有低的剪切強度,使得其摩擦系數減小,并且這種轉移膜使得后期的整個滑動過程中都具有較低的摩擦系數。在同種轉速下,隨著載荷的增加,摩擦系數也逐漸降低。這是因為隨著界面壓力的升高,接觸點數目和尺寸都進而增加,加快了表面DLC薄膜石墨化進程,降低了摩擦系數。

根據表1、2中求得的最佳嵌入維數和延遲時間計算出摩擦振動信號吸引子,其在500 r/min、450 N工況下演化過程如圖8所示。從圖8(a)~(d)可以看出,在干氣密封環摩擦60~240 s時間內,吸引子相軌跡向內收斂,其體積和軌跡曲率半徑較大,在此階段內發生不均勻磨損,處于自適應過程,這時摩擦表面只有少數粗糙微凸體發生接觸摩擦,實際接觸面積小,使得表面微凸體發生破壞、塑性變形等,摩擦系統較為無序,表現為相軌跡不穩定,呈發散狀。在圖8(e)~(j)中,即時間300~600 s時間內,吸引子相軌跡在較小的特定范圍內周而復始,軌跡呈現聚集狀態,代表此時間區間內系統達到穩定。

圖8 500 r/min、450 N工況下摩擦振動信號混沌吸引子演化Fig. 8 Chaotic attractor evolution of friction vibration signal under 500 r/min, 450 N working condition

由于前期磨合的原因,干氣密封環摩擦副表層經受較高比壓、熱效應,建立起了穩定的彈性接觸條件,磨損率變得很小。從圖7可以看出,隨著時間的推移,該工況下的摩擦系數經歷由大到小的過程,在260 s后一段時間內趨于平穩。由摩擦系數的變化曲線(由大逐漸變小)能夠明顯看出摩擦過程的變化,混沌吸引子演化軌跡的變化規律(軌跡半徑由大變小)與其相一致,因此可以利用混沌吸引子識別密封環端面的磨合狀態。

圖9為轉速500 r/min、載荷150 N下的摩擦振動信號演化。從圖9可以看出,當載荷為150 N時,在60~360 s時間內,吸引子軌跡半徑較大,吸引子軌跡收斂于中心某一點,并沿著中心往復運動,存在逐漸收斂趨勢。在480~600 s時間里,吸引子的半徑維持穩定,發散趨勢不明顯,基本收斂于一個極小空間內。通過分析可以得到,摩擦振動信號吸引子動力學演化規律和摩擦系數的變化具有一致性,而吸引子的“收斂—穩定”變化規律和摩擦過程中的磨合、穩定磨損階段相對應,可用于干氣密封環磨合過程的識別和監測。

圖9 500 r/min、150 N工況下摩擦振動信號混沌吸引子演化Fig. 9 Chaotic attractor evolution of friction vibration signal under 500 r/min, 150 N working condition

3.4 關聯維數D變化規律

圖10為轉速為500 r/min,載荷為150、450 N時,摩擦振動信號關聯維數D隨著時間的變化。在轉速為500 r/min、載荷為450 N的工況下,在磨合階段內(60~300 s),關聯維數由3.26上升至7.81;在穩定階段(300~600 s)關聯維數值在8.1左右微幅波動。從關聯維數的變化趨勢可以得出:在磨合階段,關聯維數相對較小,說明此時相空間中相關聯的點數少,相點間的距離大,吸引子的體積較大;到了正常磨損階段,關聯維數在某一較大值保持平穩,摩擦系統維持在動平衡態,吸引子收斂穩定。

圖10 關聯維數隨時間變化Fig. 10 Correlation dimension changes with time

當轉速為500 r/min、載荷為150 N時,在磨合期到穩定期D值分別從1.69上升至6.03后,小幅度上下波動,關聯維數呈現出先上升,然后保持在較大的值的趨勢變化,這種規律同載荷為450 N時的摩擦振動信號相一致。這是因為密封環磨合階段即為摩擦動力學系統的自組織階段,所謂的“自組織”就是在某一空間中依靠自身發展所形成空間有序結構。在摩擦的初期,系統較為“粗放”,沒有精細復雜成分,隨后系統向穩定狀態緩慢過渡,在內在機制的驅使下從“粗糙”向“精致”發展,系統微小尺度的變化增加,從而關聯維數呈現上升趨勢。隨著磨合階段的結束,摩擦學系統為近似有序的隨機狀態,這個過程被稱為磨損混沌狀態,即處于吸引子附近的混沌震蕩狀態[23]。在這個階段里,密封環端面間的摩擦磨損行為維持在平穩有序的結構狀態上,溫度分布均勻,摩擦系數變化穩定,密封環磨損表面相互適應,達到某種動態平衡;關聯維數穩定在一個較大值,系統處于一種精密復雜的結構形式。

關聯維數D計算了吸引子間的復雜性程度,對吸引子運動軌跡的不均勻性及其動態結構進行了良好表征。這個參數不僅能將吸引子的變化過程進行定量化,同時也可以反映系統中磨合狀態的變化,識別監測密封端面間的磨合程度。

混沌吸引子、關聯維數、摩擦系數3個參數相互關聯,對磨合狀態的評定結論一致。混沌吸引子可定性判斷磨合情況,而關聯維數、摩擦系數則可定量判斷。關聯維數與摩擦系數呈負相關,摩擦初期,吸引子發散,吸引子軌跡曲率半徑大,關聯維數小,摩擦系數大。隨著時間推進,吸引子逐漸收斂,吸引子軌跡曲率半徑降低,關聯維數變大;摩擦進入向磨合狀態時極度收斂,軌跡曲率半徑基本不變,關聯維數穩定在某一值附近。磨合狀態以后,3個參數都趨于穩定,基本保持不變,上下微弱浮動。

4 結 論

1)原始數據時域波形雜亂無序,沒有規律可循,提取出的摩擦振動信號幅值減小、呈現更明顯的規律變化。因此,通過集合經驗模式分解方法(EEMD)能夠有效提取干氣密封環在磨合磨損過程中的摩擦振動信號。

2)摩擦振動信號具有混沌特性。密封環端面間在從跑合階段到穩定階段的摩擦振動吸引子變化趨勢為“收斂—穩定”;關聯維數隨時間的增加由小變大,到磨合階段,達到平穩,上下微弱浮動。這些皆與摩擦系數變化規律具有較好的一致性。因此,混沌吸引子和關聯維數可反映密封環端面間的磨合過程。

3)可利用摩擦振動信號的關聯維數分析干氣密封環在不同磨損階段的定量演化規律,從而進一步對密封端面的磨損進行監控預警。

猜你喜歡
關聯振動信號
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奇趣搭配
中立型Emden-Fowler微分方程的振動性
智趣
讀者(2017年5期)2017-02-15 18:04:18
基于LabVIEW的力加載信號采集與PID控制
主站蜘蛛池模板: 日本在线免费网站| 亚洲欧洲日韩综合| 任我操在线视频| 午夜人性色福利无码视频在线观看| 四虎精品国产AV二区| 直接黄91麻豆网站| 天天色天天综合| 2048国产精品原创综合在线| 久久精品午夜视频| 99在线国产| 国产区免费精品视频| аv天堂最新中文在线| 日韩一级毛一欧美一国产| 国国产a国产片免费麻豆| 无码aaa视频| 亚洲天堂视频网站| 亚洲A∨无码精品午夜在线观看| 国产91在线免费视频| 国国产a国产片免费麻豆| 国产一区二区三区免费观看| 美女被躁出白浆视频播放| 97视频在线精品国自产拍| 国产一区二区三区夜色| 亚洲综合九九| 国产黄网站在线观看| 亚洲色图综合在线| 国产熟睡乱子伦视频网站| 999精品免费视频| 国产一区二区精品福利| 欧美精品高清| www成人国产在线观看网站| 国产精品蜜臀| 色吊丝av中文字幕| 国产在线一区视频| 国产无遮挡裸体免费视频| 无码'专区第一页| 超清无码熟妇人妻AV在线绿巨人| 91精品免费高清在线| 波多野结衣视频网站| 毛片a级毛片免费观看免下载| 日本欧美中文字幕精品亚洲| 久久77777| 亚洲视频影院| 日韩在线成年视频人网站观看| 伊人成人在线视频| 超级碰免费视频91| 久久a级片| 久久精品视频一| 亚洲国产精品美女| 国产美女精品人人做人人爽| 国产主播在线一区| 日韩精品资源| 国产精品青青| 国产精品久久久久久久久kt| 性色在线视频精品| 日韩黄色大片免费看| 国产va在线观看免费| 亚洲人成成无码网WWW| 日韩无码视频专区| 国产成人综合久久精品下载| 色婷婷色丁香| 欧洲日本亚洲中文字幕| 毛片卡一卡二| 国产成人乱无码视频| 成人国产三级在线播放| 亚洲欧美自拍中文| 一区二区三区成人| 免费亚洲成人| 激情国产精品一区| 91免费国产在线观看尤物| 国产91小视频在线观看| 国产不卡在线看| 亚洲成A人V欧美综合天堂| 无码人妻免费| 亚洲精品欧美日韩在线| 色综合中文综合网| 国产aⅴ无码专区亚洲av综合网| 欧美精品啪啪| 免费国产小视频在线观看| 国产精品视频导航| 欧洲在线免费视频| 自慰网址在线观看|