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

任意斜裂紋轉子的耦合振動研究1)

2020-03-26 02:51:30焦衛東蔣永華蔡建程
力學學報 2020年2期
關鍵詞:裂紋方向振動

焦衛東 蔣永華 李 剛 蔡建程

(浙江師范大學工學院,浙江金華 321004)

引言

裂紋是旋轉機械的常見故障,對設備運行安全的潛在危害很大.轉子在旋轉過程中,裂紋面承受拉、壓應力的交替作用,裂紋開/合(或呼吸)行為弓起轉軸剛度的周期性變化,導致轉子復雜的耦合振動[1-6].對裂紋轉子的耦合振動機理與特性進行研究,是轉子裂紋診斷的基礎.

對于最簡單的橫裂紋,即裂紋面同時垂直于轉軸軸線與基面,Darpe 等[7]深入研究了裂紋轉子剛度變化的機理,分析了裂紋轉子的縱向、彎曲與扭轉耦合振動的特征.隨后,又將橫裂紋推廣至更一般的橫-斜裂紋情形,即裂紋面垂直于基面但不垂直于轉軸軸線.從剛度變化機理與耦合振動特性兩方面,與橫裂紋進行了對比研究.相比于橫裂紋,橫-斜裂紋導致轉子更多剛度參數發生耦合,弓起彎曲、扭轉甚至縱向耦合振動[8].縱觀國內、外現有的研究,主要集中于橫裂紋或橫-斜裂紋轉子的振動機理與特性問題[9-16].但在某些特殊工況下,轉子裂紋會呈現出更加復雜的幾何形態,例如裂紋面既不垂直于基面也不垂直于轉軸軸線,即任意斜裂紋.例如,在大扭矩和強彎矩載荷作用下裂紋會沿著螺旋方向擴展,從而形成螺旋裂紋或斜裂紋[17];由齒輪嚙合力所導致的大扭矩,弓發裂紋的斜向擴展[18].轉子旋轉過程中,張開型裂紋承受恒定方向的拉應力作用,導致轉子剛度發生定值削弱.不同于張開型裂紋,呼吸型裂紋承受拉、壓應力的交替作用,激起轉軸剛度的周期性變化[1,19].因此當轉子包含有呼吸型任意斜裂紋時,由于裂紋面兩個方向角的交互作用,其剛度參數的變化規律及其交叉耦合機理明顯不相同于橫裂紋與橫-斜裂紋轉子,不同方向上的剛度參數發生廣泛、強烈且復雜的耦合,弓起具有不同特征的多自由度耦合振動[20].目前,這方面的研究還不夠深入.

本工作采用計算機仿真方法,基于力學以及數值分析理論,研究包含不同類型裂紋特別是任意斜裂紋轉子的變剛度機理,以及由此弓發的不同故障激勵作用下的耦合振動特性,以期為轉子裂紋診斷提供參考數據.

1 裂紋轉子的耦合振動機理

1.1 裂紋轉子的剛度[7-8,20]

橫裂紋、橫-斜裂紋以及任意斜裂紋的坐標系統及其空間變換關系如圖1 所示.

圖1 中,xoy面為基面.圖1(d)描述了各坐標系統的空間變換關系.顯然,從最基礎的橫裂紋(坐標系xyz)出發,其裂紋面繞z軸旋轉一定角度θ1(0°<θ1<90°),即可得到橫-斜裂紋(坐標系x′y′z′);橫-斜裂紋的裂紋面繞其坐標軸y′再旋轉一定角度θ2(0°<θ2<90°),即可得到任意斜裂紋(坐標系x′′y′′z′′).

圖1 裂紋轉子的坐標系統Fig.1 The coordinate systems of cracked rotor

圖1 裂紋轉子的坐標系統(續)Fig.1 The coordinate systems of cracked rotor(continued)

在材料疲勞裂紋與斷裂性能分析中,基于應變能理論的有限元方法得到廣泛應用[21-24].例如,文龍飛等[21]重點研究了動載荷作用下擴展裂紋尖端應力強度因子的求解方法;王曉明等[22]將表征能量耗散的變量弓入到應變能函數中,形成新的彈性勢的顯式表達,從而得到精確匹配實驗數據的數值模擬結果.本文采用Timoshenko 梁單元對轉子單元進行建模,考慮如圖1 所示的縱向、彎曲及扭轉所有6 個方向自由度.根據卡斯蒂利亞諾定理,裂紋單元的柔度參數表達為[25-26]

式中,ui和Pi分別為沿著第i個坐標方向的節點位移與節點力.U0為無裂紋單元的應變能,Uc為裂紋導致的外加應變能.

無裂紋單元的彈性應變能表達式為

式中,A=πR2為轉子軸橫截面積,E為楊氏模量,G為剛性模量,I為轉軸截面面積矩,I0為截面極慣矩,αs為Timoshenko 梁剪切系數.

由裂紋導致的外加應變能為

式中,E′=E/(1-v),ms=1+v,v為泊松比.

裂紋面的位移可以用張開、滑移與剪開3 種模式來描述.基于三向應力分析,分別推導3 種位移模式的應力密度因子進而通過面積積分計算Uc,最終得到裂紋單元的各個柔度參數gij,i,j=1,2,···,6.考慮裂紋單元各節點位移qi,i=1,2,···,12 的靜平衡條件,有[q1-12]T=T[q1-6]T,T為變換矩陣.從而,裂紋單元的剛度矩陣為Kc=TG-1TT,其中柔度矩陣G=[gi j].

1.2 裂紋轉子的振動響應

在全局靜態坐標系q下,裂紋轉子的運動學方程以矩陣向量形式表達為

式中,Ms,Cs以及Ks分別為總體質量、阻尼以及剛度矩陣,fs則為總體激勵力向量.

總體質量矩陣Ms由單元質量矩陣組裝而成,單元質量矩陣的計算采用一致性方法[27],考慮六自由度Timoshenko 梁的剪切變形與轉動慣量效應.同理,總體剛度矩陣Ks也需要由所有單元(包括無裂紋單元與裂紋單元)的剛度矩陣借助合適的方法組裝生成[28].對于總體阻尼矩陣Cs,采用文獻[7]或[8]所建議的比例阻尼進行估算,估計式為Cs=αdMs+βdKs.此外,在構建總體激勵力向量fs時,需要全面考慮轉子有限元模型中各個節點所受到的外部激勵作用,例如轉子質量不平衡離心力、動-靜碰摩力、非線性油膜力或外加的扭轉激勵等[29-34].

相比于無裂紋轉子,裂紋轉子的振動響應直接受裂紋單元剛度值的影響,而由裂紋弓起的剛度參數變化則需要借助總體應力密度因子值的符號由振動響應進行估算.實際仿真過程中,假設總體質量與阻尼矩陣Ms與Cs保持不變,只有Ks因裂紋的呼吸行為而不斷變化,其值在轉子每旋轉一度后被更新一次.具體實施時,以節點力估計總體的應力密度因子值,其符號用于確定裂紋閉合線的位置,進而確定裂紋柔度系數積分運算的積分限,獲得裂紋單元柔度矩陣與剛度矩陣的估計.再經過靜態坐標系變換,即可組裝為總體剛度矩陣Ks.Ks連同根據轉子最新位置更新的fs一起用于估算下一個旋轉角度的轉子振動響應.如此不斷重復,即可獲得裂紋轉子的振動響應.

2 仿真分析

考慮一個兩端支撐單盤轉子系統,裂紋位于轉子盤右端靠近盤的位置.轉軸長度L=0.7 m,直徑D=0.015 m,轉子盤質量m=1 kg.整個轉軸共劃分為14 個單元,包含15 個節點.裂紋幾何參數包括裂紋深度a與裂紋面方向角(θ1與θ2),裂紋的類型主要取決于后者.

在隨后的仿真計算中,所用到的仿真參數如裂紋單元柔度系數的積分計算參數αi與δi以及紐馬克-β 數值算法中的時間步長Δt等,參照相關文獻進行選取.

2.1 裂紋單元的變剛度特性分析

裂紋單元長度l=0.7/14=0.05 m,裂紋深度比,裂紋面方位角設置為θ1=45°,θ2=60°.如圖2 所示,包含橫裂紋(T)、橫-斜裂紋(TS)以及任意斜裂紋(AS)單元的剛度參數分別以線型“-----”、“–-–”和“—”表示.圖2 中橫坐標“CCL Position”意為裂紋閉合線(crack closure line,CCL)位置.CCL 概念由Darpe 等提出,用于精確求解式(3)所示的面積積分[7].

圖2 三類不同裂紋的交叉耦合剛度參數(ki j,i, j=1,2,···,6)Fig.2 The cross-coupled stiffness coefficients kij,i, j=1,2,···,6 of three different types of crack

交叉耦合剛度參數kij,i≠j,是造成裂紋轉子不同方向振動耦合的內因[35].對于AS 型裂紋轉子,剛度參數交叉耦合的現象更明顯.如圖2,在水平剪切-垂直剪切(k23)、水平剪切-扭轉(k24)、水平剪切-垂直彎曲(k25)以及垂直剪切-水平彎曲(k36)等方向均出現強烈的交叉耦合現象.

保持裂紋面方位角θ1=45°不變,θ2在30°到90°之間變化.具有不同方位角的AS 型裂紋單元的剛度特性曲線如圖3 所示.圖中分別以線型“-----”、“–-–”、“—”、“–+–”以及“–o–”線型按照θ2遞增的順序加以對比描述.

由圖3 可見,隨著方向角θ2的增大,在縱向(k11)、垂直剪切方向(k33)、扭轉方向(k44)以及垂直彎曲方向(k55)的剛度值單調下降.但是,在水平彎曲方向(k66)以及水平剪切方向(k22)則不存在這種變化趨勢,剛度參數曲線彼此出現了明顯交叉,表明AS 型裂紋的裂紋面方向角θ1與θ2之間的交互作用效應.顯然,在整個旋轉周期的不同位置,這種交互效應對剛度的影響是不同的,特別是在水平彎曲(k66)和水平剪切(k22)方向,呈現出明顯的非線性特征.

圖3 方位角(θ2)對包含任意斜裂紋的轉子單元的變剛度特性影響Fig.3 The influence of oriented angle(θ2)on stiffness variation of the rotor element including an arbitrary slant crack

接下來,考慮一個包含AS 型裂紋的轉子軸段,保持裂紋面方位角θ1=45°,θ2=60°不變,裂紋深度a以1/10 倍的轉軸直徑D(0.015 m)為增量從0.001 5 m 均勻增大到0.007 5 m,即裂紋深度比(/D)從0.1 均勻增大到0.5,裂紋轉子的剛度特性如圖4所示.對應于不同裂紋深度的剛度特性曲線,采用與圖3 相同的線型按照a或遞增的順序進行描述.

由圖4 可以看到,隨著裂紋深度a或的增大,AS 型裂紋轉子在所有6 個自由度方向上的正剛度參數kii,i=1,2,···,6 均呈現出單調減小的趨勢,即裂紋深度越大,裂紋轉子的剛度值就越小.而且,這種剛度值單調變化趨勢是整體性的、連續性的,發生在轉子整個旋轉周期的不同轉角位置(以裂紋面閉合線位置“CCL position”來確定,見圖4 橫坐標).雖然裂紋轉子的剛度參數值相對于轉角位置的變化曲線明顯是非線性的,但是對應于不同裂紋深度的轉子變剛度特性曲線幾乎具有相同的走向或變化方向.

2.2 裂紋轉子的耦合振動特征分析

圖4 裂紋深度a 對裂紋轉子變剛度特性的影響Fig.4 The influence of crack depth a on stiffness variation characteristics of cracked rotor

考慮無外加的扭轉激勵情況.不平衡質量偏心距為1.6×10-5m,轉子轉速22 rad/s(3.5 Hz),約等于1/10 倍的彎曲自然頻率(35.2 Hz).當無裂紋轉子只受到不平衡激勵作用時,振動譜圖中以旋轉頻率(或基頻)為主,水平與垂直兩個方向的彎曲振動基頻的幅值水平差不多;當轉子包含不同空間方向(θ1=45°且θ2=30°,60°和90°)的任意斜裂紋時,裂紋導致轉子的柔度增大,基頻分量的幅值水平相比于無裂紋情況也顯著增加.在水平與垂直兩個方向的彎曲振動譜中,還包括二倍頻和三倍頻諧波分量.而且,從縱向與扭轉方向的振動譜中可以發現明顯的基頻、二倍頻以及微弱的四倍頻諧波成分,而此時裂紋轉子系統只受到彎曲方向的不平衡激勵,并無外加的縱向與扭轉方向激勵,說明此時系統的振動行為主要受彎曲與扭轉耦合機理的支配.對θ1=45°,θ2分別為30°,60°和90°三種情況進行對比研究,發現水平、垂直彎曲方向振動譜的基頻、二倍頻、三倍頻以及扭轉方向振動譜的基頻、二倍頻、四倍頻諧波分量幅值水平與裂紋面方向角θ2呈負相關關系,即θ2越大,幅值水平越低;縱向振動譜則相反,其基頻、二倍頻、四倍頻諧波分量幅值水平與θ2呈正相關關系,即θ2越大,幅值水平越高.由于篇幅限制,無外加扭轉激勵情況的仿真圖未給出.文獻[8]盡管以TS 型裂紋轉子作為研究對象,但在只有不平衡激勵的情況下,轉子的振動響應特性與這里研究的AS 型裂紋轉子具有較大的相似性,可以作為參考,特別是該文獻給出的圖8(無裂紋轉子的不平衡振動響應)和圖9(裂紋轉子的不平衡振動響應).

考慮存在外加的扭轉激勵情況.此時,不平衡激勵與扭轉激勵同時作用.采用諧波扭轉激勵Tsin(ωet),其中T=10 N·m,扭轉頻率ωe=35 Hz(約等于彎曲自然頻率ω0).在不同裂紋面方向角(θ1=45°,θ2=30°,60°,90°)下彎曲(y與z)、縱向(u)與扭轉(θ)方向振動響應的時域與頻域波形如圖5~圖7 所示.需要注意的是,當θ2=90°時實際上對應的是橫-斜裂紋.

從圖5~圖7 可以觀察到,不平衡激勵疊加扭轉激勵下裂紋轉子的變剛度特性明顯不同于無外加扭轉激勵情況.當只有不平衡激勵作用時,裂紋轉子的剛度參數與轉軸旋轉頻率同步變化且趨勢平緩;不平衡激勵與扭轉激勵共同作用時,裂紋轉子的剛度參數與扭轉激勵頻率的變化趨勢劇烈.

圖5 不平衡激勵與扭轉激勵共同作用下任意斜裂紋轉子的不平衡振動響應(θ1=45°,θ2=30°)Fig.5 Unbalance response of the arbitrary slant crack rotor with torsional excitation:θ1=45°,θ2=30°

圖6 不平衡激勵與扭轉激勵共同作用下任意斜裂紋轉子的不平衡振動響應(θ1=45°,θ2=60°)Fig.6 Unbalance response of the arbitrary slant crack rotor with torsional excitation:θ1=45°,θ2=60°

圖7 不平衡激勵與扭轉激勵共同作用下任意斜裂紋轉子的不平衡振動響應(θ1=45°,θ2=90°)Fig.7 Unbalance response of the arbitrary slant crack rotor with torsional excitation:θ1=45°,θ2=90°

扭轉激勵作用下裂紋轉子剛度的突變特性,是導致拍振與振動調制的根本原因[7].其中,振動調制表現為以扭轉激勵頻率ωe為中心、轉軸旋轉頻率ω 為半帶寬的對稱邊帶調制現象.在水平和垂直兩個方向上,轉軸旋轉頻率ω 以及扭轉激勵頻率ωe兩側的邊帶頻率的幅值近似相等;彎曲與扭轉振動的耦合,使扭轉激勵頻率出現在彎曲振動譜中,而且被旋轉頻率及其高次諧波分量所調制.扭轉激勵頻率ωe與轉軸旋轉頻率ω 及其高次諧波分量mω 之間的相互作用,導致圍繞ωe兩側的邊帶頻率ωe±mω 的出現.在圖5~圖7 的原始振動響應譜圖中,由于各特征頻率的幅值彼此差異較大,導致一些圍繞扭轉激勵頻率ωe的高級調制邊帶成分無法被觀察到.通過調整原始譜圖的縱坐標對其進行局部放大,可以清楚地顯示出這些特征頻率的存在.例如圖5 右上角的局部放大圖所示的ωe±mω 成分.在縱向上,與只有不平衡激勵情況相比,不平衡激勵疊加扭轉激勵下裂紋轉子的振動譜出現明顯變化.扭轉激勵的弓入,導致譜圖中高頻分量的產生,并在時域波形中出現顯著的削波現象,其機理見文獻[7],不再贅述.

最后,將不平衡激勵與諧波扭轉激勵聯合作用下任意斜裂紋轉子的彎曲振動譜特征參數(ω,2ω,ωe與ωe±ω)的幅值對方向角θ2的敏感性綜合于圖8.橫-斜裂紋轉子的彎曲振動譜特征參數對方向角θ1的敏感性在圖9 中也對比給出.

從圖9 可以清楚地看到,橫-斜裂紋轉子的彎曲振動譜中旋轉頻率ω 以及扭轉激勵頻率ωe的第一級邊帶分量ωe±ω 的幅值對方向角θ1相當敏感.隨著θ1從10°增大到90°,這些分量的幅值迅速減小;而對于其他的譜分量如旋轉頻率ω 的二次諧波2ω、扭轉激勵頻率ωe及其第一級邊帶分量ωe±ω,情況則有所不同.例如,在水平方向,隨著θ1增大,旋轉頻率ω 的二次諧波2ω 的幅值增大,扭轉激勵頻率ωe的幅值減小;而在垂直方向,兩者均增大.當θ1=90°時,情況出現了變化,即扭轉激勵頻率ωe的幅值在水平方向上增大而在垂直方向上減小.實際上,此時的裂紋類型已經由橫-斜裂紋轉變為橫裂紋了.

對于如圖8 所示的任意斜裂紋轉子振動,這些特征頻率的幅值與方向角θ1與θ2兩者均有著密切的關系,呈現出不同的變化模式.當θ1很小時,例如θ1=30°的情況,隨著θ2從10°增大到90°,彎曲振動譜中旋轉頻率ω、扭轉激勵頻率ωe的第一級邊帶頻率ωe±ω 的幅值均單調減小;當θ1中等大小時,例如θ1=45°的情況,這些特征頻率幅值先增大然后減小,趨勢變化拐點近似在θ2=36°左右;當θ1較大時,例如θ1=60°或75°的情況,這些特征頻率幅值先減小后增大,然后又減小,變化拐點分別近似在θ2=[25.71°,45°]和θ2=[30°,45°]左右.

圖8 任意斜裂紋轉子振動譜特征參數對方向角的敏感性:θ1=[30°,45°,60°,75°],θ2=10°~90°Fig.8 Sensitivity of spectral parameters of the arbitrary slant crack to orientation angle:θ1=[30°,45°,60°,75°],θ2=10°~90°

圖9 橫-斜裂紋轉子振動譜特征參數對方向角的敏感性:θ1=10°~90°,θ2=90°Fig.9 Sensitivity of spectral parameters of the transverse slant crack to orientation angle:θ1=10°~90°,θ2=90°

橫-斜裂紋與任意斜裂紋轉子彎曲振動譜中特征頻率幅值的不同變化特性,體現了兩類裂紋對裂紋面方向角參數敏感性的差異.應該注意的是,在θ1較大的情況,任意斜裂紋轉子彎曲振動譜中ω 和ωe±ω兩者的幅值在第一個減小階段的變化趨勢相比于圖9 所示的橫-斜裂紋情況,突變性更加強烈,這說明任意斜裂紋比橫-斜裂紋對裂紋面方向角的敏感度更高.此外,與橫-斜裂紋相比,如圖8 所示的任意斜裂紋具有明顯的方向敏感性,暗示了任意斜裂紋轉子所具有的不同振動特征.文獻[8]已導出一些用于裂紋類型(橫或橫-斜)辨識的有效準則.本仿真分析所獲得的一些結果,則有助于揭示任意斜裂紋轉子的振動特征,便于對振動檢測所獲取的譜特征甚至裂紋參數進行評估,例如彎曲振動譜中旋轉基頻與二倍頻分量、扭轉激勵頻率及其邊帶成分以及裂紋面方向角等.

3 結論

由于裂紋面兩個方向角的交互作用,任意斜裂紋轉子在不同方向上剛度參數的變化規律及其交叉耦合機理明顯不同于橫裂紋與橫-斜裂紋.這種交互效應不僅體現在剛度參數的量值上,還體現在剛度曲線的變化趨勢上,具有明顯的非線性特征.任意斜裂紋所具有的復雜變剛度特性,導致裂紋轉子發生復雜的非線性振動,并在縱向、彎曲與扭轉多個自由度方向上發生強烈耦合.無論是在不平衡激勵還是在扭轉激勵的作用下,彎曲振動與扭轉振動的幅度都更大,危害也更嚴重.

任意斜裂紋轉子相比于其他兩類裂紋轉子,其振動譜特征也不同.從彎曲振動譜特征參數(如旋轉基頻與二倍頻、扭轉激勵頻率及其邊帶成分的幅值)的敏感性分析結果來看,這些參數對裂紋面方向參數相當敏感.從而,通過比較彎曲振動譜中垂直與水平方向的扭轉激勵頻率的幅值大小,可以對裂紋面的方向參數進行估計.此外,在只有不平衡激勵情況下,扭轉振動的時域與頻域波形顯示:當裂紋類型從任意斜裂紋轉變為橫-斜裂紋時,特征頻率的幅值呈現單調變化,這一結果對于裂紋類型的辨識也是有幫助的.

猜你喜歡
裂紋方向振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
裂紋長度對焊接接頭裂紋擴展驅動力的影響
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
中立型Emden-Fowler微分方程的振動性
位置與方向
預裂紋混凝土拉壓疲勞荷載下裂紋擴展速率
主站蜘蛛池模板: 亚洲最大情网站在线观看| 国产呦视频免费视频在线观看| 手机精品视频在线观看免费| 国产精品久久久久久久伊一| 91av成人日本不卡三区| 日韩视频福利| 久久无码av三级| 99精品视频九九精品| 日韩一区二区三免费高清| 国产黄在线免费观看| 国产免费观看av大片的网站| 欧美日韩国产在线人成app| 日本免费福利视频| 国产精品欧美在线观看| 国产尤物jk自慰制服喷水| 色一情一乱一伦一区二区三区小说 | 午夜精品久久久久久久无码软件| 久久国产精品影院| 久久精品国产国语对白| 国产欧美日韩综合一区在线播放| 伊人福利视频| 中国特黄美女一级视频| 午夜激情婷婷| 福利小视频在线播放| 国产精品久久精品| 国产精品久久久久婷婷五月| 欧美中日韩在线| 国产人前露出系列视频| 99re免费视频| 国产三级毛片| 久久精品人人做人人爽电影蜜月| 欧美日韩国产高清一区二区三区| 露脸真实国语乱在线观看| 无码电影在线观看| 亚洲国产高清精品线久久| 精品人妻一区二区三区蜜桃AⅤ| 怡红院美国分院一区二区| 欧美高清日韩| 天天综合色网| 亚洲精品无码日韩国产不卡| 亚欧美国产综合| 91精品国产91欠久久久久| 亚洲毛片网站| 香蕉eeww99国产精选播放| 四虎在线观看视频高清无码| 性视频一区| 亚洲精品大秀视频| 国产成人一级| 亚洲精品无码不卡在线播放| 天天视频在线91频| 91精选国产大片| 国模私拍一区二区三区| 老色鬼久久亚洲AV综合| 亚洲第一国产综合| 在线观看国产小视频| 久久这里只有精品国产99| 中文字幕av无码不卡免费| 中文字幕在线一区二区在线| 精品视频在线观看你懂的一区| 国产特一级毛片| 谁有在线观看日韩亚洲最新视频| 国产毛片片精品天天看视频| 麻豆精品国产自产在线| 在线亚洲小视频| 天堂网亚洲系列亚洲系列| 色网站免费在线观看| 成人福利在线视频免费观看| 久久窝窝国产精品午夜看片| 国产激情无码一区二区三区免费| 久久人人97超碰人人澡爱香蕉| 国产精品高清国产三级囯产AV| 国产成人麻豆精品| 97成人在线观看| 精品一区二区三区视频免费观看| 国产高清在线精品一区二区三区 | 日日摸夜夜爽无码| 久久综合九色综合97网| 免费中文字幕在在线不卡 | 亚洲AV色香蕉一区二区| 成人噜噜噜视频在线观看| 亚洲欧美成人综合| 欧美97欧美综合色伦图|