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

粗糙表面微凸體模型構(gòu)建與接觸特性分析

2023-01-31 07:47:18何本帥王晶晶張錦華蔡安江
振動與沖擊 2023年2期
關(guān)鍵詞:變形模型

李 玲, 何本帥, 王晶晶, 張錦華, 蔡安江

(西安建筑科技大學(xué) 機(jī)電工程學(xué)院,西安 710055)

機(jī)械結(jié)合面之間的接觸在微觀層面上是兩個粗糙表面間的接觸問題。相互接觸的粗糙表面在機(jī)械系統(tǒng)運(yùn)行過程中起著傳遞運(yùn)動、負(fù)載和能量的重要角色,其接觸行為直接關(guān)系到機(jī)械零件之間的連接性能,對機(jī)械系統(tǒng)的靜態(tài)和動態(tài)行為都有著重要的影響[1-4]。因此,準(zhǔn)確地構(gòu)建粗糙表面間的接觸模型,開展粗糙表面的接觸特性研究,對于現(xiàn)實機(jī)械工程問題的解決具有重要意義。

為了研究粗糙表面的接觸特性,很多學(xué)者在接觸力學(xué)領(lǐng)域不斷創(chuàng)新。Greenwood等[5]基于統(tǒng)計學(xué)理論,將粗糙表面上的微凸體高度分布看成隨機(jī)函數(shù),最先提出粗糙表面的彈性接觸模型,即GW模型。后來,許多學(xué)者開始拓展GW模型,并提出了許多全新的模型。例如:Chang等[6]基于微凸體在塑性變形階段的體積守恒原理,提出了一種微凸體塑性接觸模型,并將兩者擬合獲得了彈塑性接觸特性,利用統(tǒng)計學(xué)理論獲得了整個粗糙表面的接觸模型,簡稱CEB模型;Zhao等[7]為了研究微凸體從彈性變形階段向塑性變形階段轉(zhuǎn)變的力學(xué)過程,基于接觸力學(xué)理論,利用數(shù)學(xué)多項式將微凸體的彈性行為和完全塑性行為緊密聯(lián)系起來,獲得了微凸體在三種變形階段的接觸特性曲線,簡稱ZMC模型;Kogut等[8]通過有限元方法深入研究了單個球體與剛性平面的接觸行為,并使用冪函數(shù)公式擬合獲得了彈性、彈塑性和完全塑性狀態(tài)下接觸載荷與變形量之間的解析表達(dá)式,簡稱KE模型;Jackson等[9]利用有限元的方法,研究了微凸體在變形過程中材料屬性和幾何形狀的影響,簡稱JG模型。國內(nèi)學(xué)者[10-12]對結(jié)合面的力學(xué)模型也進(jìn)行了研究和分析,促進(jìn)了粗糙表面接觸特性的研究和發(fā)展。

上述學(xué)者對于粗糙表面的接觸特性研究重點放在了接觸力學(xué)部分,所采用的微凸體輪廓假設(shè)與GW模型相同,均為半球體模型,并未涉及微凸體模型的重新建立。但國外學(xué)者[13-15]對粗糙表面的測量研究發(fā)現(xiàn),正弦形接觸比半球形接觸能更好地模擬粗糙表面的微凸體接觸,并將粗糙表面轉(zhuǎn)換成傅里葉級數(shù)或魏爾斯特拉斯曲線來考慮多尺度的微凸體,重新建立粗糙表面接觸模型。同時,國內(nèi)學(xué)者安琪等[16]針對磨削表面的微凸體輪廓,提出一種半周期余弦曲線回轉(zhuǎn)體模擬磨削表面微凸體的接觸模型;Shang等[17]提出一種采用拋物線回轉(zhuǎn)體等效微凸體的方法,推導(dǎo)出粗糙表面法向剛度與彈性接觸之間的理論關(guān)系。隨著這類模型廣泛應(yīng)用于粗糙表面的接觸特性研究,如何更加合理地優(yōu)化粗糙表面的接觸模型具有重要意義。

本文基于粗糙表面形貌測量試驗,提出一種采用二次函數(shù)回轉(zhuǎn)體等效微凸體的辦法,建立微凸體接觸半徑與接觸變形的解析關(guān)系。然后根據(jù)幾何模型,重新推導(dǎo)出單個微凸體在彈性、彈塑性和完全塑性三種變形階段的接觸表達(dá)式,并應(yīng)用接觸力學(xué)理論和概率統(tǒng)計方法,建立粗糙表面的微觀接觸模型。最后針對模型的有效性進(jìn)行了驗證和分析,并揭示了塑性指數(shù)、微凸體尺寸參數(shù)對接觸特性的影響規(guī)律。

1 粗糙表面模型構(gòu)建

1.1 粗糙表面形貌測量試驗

本文使用Talysurf非接觸式光學(xué)輪廓儀測量表面形貌,獲取不同加工工藝下的二維表面信號。利用頻譜分析法確定表面成分有效信息的頻段范圍,并把不同頻率成分的信息分解到互不重疊的頻段上,獲得試塊表面的粗糙度以及二維輪廓曲線數(shù)據(jù)。圖1為表面形貌測量試驗。圖1(a)為合金鋼長方體試塊,上表面為待測表面,待測表面區(qū)域為20 mm×20 mm,材料為20CrMo。在前立面上作有標(biāo)記,為了便于說明,將打標(biāo)立面與上表面的交邊記為參考線。圖1(b)為測量位置,每個表面分3次測量,如圖1(b)中所標(biāo)的1,2和3等三條線段,測量時由下向上,測試長度選取3 mm左右,最小采樣間距為0.25 μm。

圖1 粗糙表面形貌測量試驗Fig.1 Rough surface morphology measurement experiment

試驗共選取三種不同加工工藝的待測表面,分別為:①Ra=0.85 μm,磨削+熱處理加工表面;②Ra=1.302 μm,磨削加工表面;③Ra=1.602 μm,銑削加工表面。將獲得的表面二維輪廓曲線數(shù)據(jù)進(jìn)行整理,每個測試表面選取一組數(shù)據(jù),使用MATLAB軟件仿真,獲得不同粗糙表面的二維輪廓,如圖2所示。

圖2 不同粗糙表面的二維輪廓Fig.2 Two-dimensional profile of different rough surfaces

由圖2可知,不同加工工藝下的粗糙表面在微觀層面上均是粗糙不平的,隨機(jī)分布著許多大小不一的凹谷和凸峰,其二維輪廓曲線差異較大。從局部放大圖可以看出,雖然各表面的粗糙度和加工工藝不同,但單個微凸體的二維輪廓曲線,均呈現(xiàn)出頂部較窄,底部較寬的特征,與二次函數(shù)曲線具有很高的相似度。

1.2 單個微凸體的輪廓模型

單個微凸體的二維輪廓特征與二次函數(shù)曲線具有很高的相似度,據(jù)此微凸體的二維輪廓可用二次函數(shù)曲線進(jìn)行假設(shè),如圖3所示,l為微凸體的直徑,h為微凸體的高度。

圖3 微凸體二維輪廓模型假設(shè)Fig.3 Assumptions of the 2D profile model of asperities

根據(jù)圖3,可定義粗糙表面上單個微凸體的二維幾何尺寸為

(1)

式中,x的范圍為微凸體正負(fù)直徑的一半。參數(shù)l和h的獲取方法為:首先,分別對粗糙表面微凸體輪廓的凸峰值與凹谷值進(jìn)行標(biāo)記;然后,對所標(biāo)記的相鄰峰值點的間隔進(jìn)行求和,得到相鄰峰值點間距的平均值的一半作為微凸體的直徑l;對所標(biāo)記的相鄰峰谷值做差處理,進(jìn)而得到各峰谷差值的平均值作為微凸體的高度h。

根據(jù)測量表面的二維輪廓數(shù)據(jù),隨機(jī)選取單個微凸體的數(shù)據(jù)點,通過對比二次函數(shù)模型和半球體模型的擬合效果,來驗證二次函數(shù)模型的正確性,如圖4所示。從整體來看,二次函數(shù)擬合曲線與實測數(shù)據(jù)更吻合,尤其是在微凸體底部。半球體模型無法實現(xiàn)完全擬合的原因是采用單一的曲率半徑,在重載的情況下會對單個微凸體的解析計算帶來誤差。值得注意的是,半球體模型的擬合曲線呈現(xiàn)出半橢圓形態(tài)是因為輪廓數(shù)據(jù)點的水平方向和垂直方向的放大比率不同,產(chǎn)生了畸變現(xiàn)象。本文提出的二次函數(shù)模型實現(xiàn)了對微凸體整體輪廓的擬合,能夠在載荷較大的情況下,更加精確地進(jìn)行結(jié)合面接觸特性分析。

圖4 微凸體二維輪廓擬合驗證Fig.4 Two-dimensional contour fitting verification of asperity

綜上所述,微凸體的二維輪廓采用二次函數(shù)曲線擬合。對于微凸體的三維輪廓,可通過二次函數(shù)曲線沿著z軸旋轉(zhuǎn)一周獲得的回轉(zhuǎn)體模型進(jìn)行表征。與半球體模型的建立類似,本文模型基于以下假設(shè):①將兩粗糙表面之間的接觸等效為一剛性平面與一等效粗糙表面的接觸;②微觀接觸表面的微凸體輪廓為二次函數(shù)回轉(zhuǎn)體,其高度分布符合高斯分布;③在接觸變形過程中只考慮微凸體變形,宏觀基體變形和微凸體之間的相互作用不予考慮。

2 粗糙表面理論建模

2.1 單個微凸體的微觀接觸

圖5為單個微凸體與剛性平面之間的接觸原理。圖5中:d為接觸間隙;r為接觸半徑;δ=h-d為接觸變形。

注:虛線表示微凸體的初始形狀;實線表示施加載荷發(fā)生接觸變形后的形狀。圖5 單個微凸體微觀接觸模型Fig.5 Micro-contact model of single asperity

結(jié)合圖5與式(1),可得單個微凸體的接觸半徑r與接觸變形δ的關(guān)系為

(2)

隨著接觸載荷的逐步增大,接觸變形δ和接觸半徑r也將隨之增大。當(dāng)接觸載荷逐漸增大到超過材料的屈服極限時,微凸體的變形將逐漸由彈性向彈塑性并最終向完全塑性轉(zhuǎn)化,下面將分階段對微凸體的變形過程進(jìn)行分析。

2.1.1 彈性變形階段

當(dāng)δ<δc,接觸變形較小時,微凸體處于彈性變形階段,可視為一個等效曲率半徑為R的半球體與剛性平面的彈性接觸,這與半球體模型假設(shè)類似。因此,微凸體頂端可采用HERTZ接觸理論[18]的處理辦法,求得彈性變形階段的平均接觸壓力pe的表達(dá)式為

(3)

根據(jù)Tabor[19]的研究結(jié)論,在最大HERTZ壓力等于0.6H,即pe=0.4H時微凸體開始出現(xiàn)塑性變形。對于一般的情況,開始出現(xiàn)塑性變形的條件為

(4)

式中:k為最大平均接觸壓力的系數(shù),k=0.454+0.41v,v為較軟材料的泊松比;H為較軟材料的硬度。

將式(4)代入式(3),可推導(dǎo)出材料初始塑性屈服點的臨界變形量δc為

(5)

根據(jù)二次函數(shù)回轉(zhuǎn)體模型,求得微凸體頂點輪廓線的曲率半徑re為

(6)

根據(jù)HERTZ接觸理論進(jìn)行推導(dǎo),將半球體模型的等效曲率半徑R用二次函數(shù)回轉(zhuǎn)體模型的頂點輪廓線曲率半徑re代替,進(jìn)而彈性變形階段的接觸面積ae的表達(dá)式為

(7)

結(jié)合式(3)和式(7),整理可得彈性變形階段接觸載荷fe的表達(dá)式為

(8)

由式(8)可得,彈性變形階段接觸剛度ke的表達(dá)式為

(9)

2.1.2 完全塑性變形階段

當(dāng)δ>δp時,微凸體將發(fā)生完全塑性變形。δp為材料在完全塑性屈服臨界點對應(yīng)的變形量,δp=tδc。根據(jù)文獻(xiàn)[20-21]的研究,選擇比率t=110。在這個階段,平均接觸壓力pp為較軟材料的硬度[22],是一個定值,即

pp=H

(10)

在回轉(zhuǎn)體模型接觸區(qū)域近似圓形的情況下,可根據(jù)Johnson[23]的研究將完全塑性變形階段的接觸面積ap表示為

(11)

式中,rp為完全塑性變形階段微凸體的接觸半徑。

結(jié)合式(10)和式(11),完全塑性變形階段接觸載荷fp的表達(dá)式為

(12)

由式(12)可得,完全塑性變形階段接觸剛度kp的表達(dá)式為

(13)

2.1.3 彈塑性變形階段

當(dāng)δc<δ<δp時,微凸體既發(fā)生彈性變形,也發(fā)生塑性變形,平均接觸壓力與接觸變形的關(guān)系將變得極為復(fù)雜,無法進(jìn)行準(zhǔn)確的定量描述。采用改進(jìn)的HERMIT多項式插值方法[24]進(jìn)行過渡處理,增強(qiáng)平均接觸壓力在彈塑性變形始末的連續(xù)性和光滑性,其表達(dá)式為

(14)

根據(jù)文獻(xiàn)[25]的方法,彈塑性變形階段的接觸面積aep可通過構(gòu)建樣板函數(shù),實現(xiàn)其在彈塑性變形始末的連續(xù)、光滑和單調(diào)性要求,表達(dá)式為

(15)

式中,C1和C2分別為二次項系數(shù)和三次項系數(shù),且同時滿足以下條件

(16)

將式(18)代入式(17)中,可得

(17)

根據(jù)式(14)和式(15)推導(dǎo)出彈塑性變形階段接觸載荷fep的表達(dá)式為

(18)

由式(18)可得,彈塑性變形階段接觸剛度kep的表達(dá)式為

(19)

2.2 粗糙表面之間的微觀接觸

兩個粗糙表面在接觸過程中,有如下定義:名義接觸面積An上如果有N個微凸體,則所有接觸的微凸體數(shù)目可以表示為

(20)

式中:d為給定的接觸表面平均距離;η為微凸體的面積密度;An為名義接觸面積;Ф(z)為微凸體高度分布的概率密度函數(shù),其表達(dá)式為[26]

(21)

(22)

2.2.1 彈性變形階段

(23)

(24)

(25)

2.2.2 彈塑性變形階段

(26)

(27)

(28)

2.2.3 完全塑性變形階段

(29)

(30)

(31)

整個粗糙表面的接觸面積、接觸載荷和接觸剛度分別是所有接觸的微凸體在三個變形階段的累加和。且為了確保所得到的結(jié)果具有廣泛適用性,分別用An,EAn和σ/EAn對上述接觸參數(shù)進(jìn)行無量綱處理,獲得無量綱接觸面積A*、無量綱接觸載荷F*和無量綱接觸剛度K*的表達(dá)式,即

(32)

(33)

(34)

3 模型驗證

為驗證所建模型的有效性,基于圖2中粗糙表面(Ra=1.302 μm)的實測數(shù)據(jù),材料彈性模量E1=E2=211 GPa,泊松比v1=v2=0.28,材料硬度H=1.96 GPa,微凸體的直徑l=16.790 μm,微凸體的高度h=1.082 μm。選取表1中[27]的塑性指數(shù)ψ=2.0,根據(jù)式(6)求得σ=0.037 μm,η=4.459×1010m-2,進(jìn)行模型驗證。

表1 不同塑性指數(shù)的表面參數(shù)Tab.1 Surface parameters for different plasticity indices

3.1 單個微凸體接觸模型對比驗證

圖6為單個微凸體接觸面積a、接觸載荷f分別與變形量的關(guān)系,比較了本文模型與CEB模型、ZMC模型和KE模型的預(yù)測結(jié)果。如圖6所示,本文模型的曲線均呈現(xiàn)連續(xù)光滑地變化趨勢,且與對比模型的趨勢是一致的。隨著變形量的增加,本文模型逐漸接近ZMC模型和KE模型,最終在完全塑性階段相互重合。其中,CEB模型和KE模型的曲線都發(fā)生了跳躍的現(xiàn)象,這是由于CEB模型忽略了微凸體的彈塑性變形階段,平均接觸壓力從彈性階段的2/3kH增大到塑性階段的kH;KE模型在曲線擬合過程中忽略了兩個屈服臨界點處曲線的平滑性。此外,本文模型和ZMC模型均采用樣板函數(shù),所以曲線是連續(xù)光滑變化的。值得注意的是,上述模型的微凸體接觸半徑均隨著接觸變形而變化的,但本文基于二次函數(shù)回轉(zhuǎn)體的微凸體建模實現(xiàn)了對微凸體整體輪廓的擬合,能夠在載荷較大的情況下,更加精確地進(jìn)行結(jié)合面接觸特性分析。

圖6 不同微凸體接觸模型對比Fig.6 Comparison of different asperity contact models

3.2 粗糙表面接觸模型對比驗證

圖7為不同粗糙表面接觸模型的無量綱接觸面積A*、無量綱接觸載荷F*分別與無量綱平均距離d*的關(guān)系。圖7(a)和圖7(b)顯示,所有模型預(yù)測的接觸曲線均隨著平均距離的增大而逐漸減小為零。在平均距離較大時,粗糙表面的微凸體大部分處于彈性變形階段,四種模型的差距很小。當(dāng)平均距離為零時,四種模型的差距最大。從小到大來看,接觸面積依次為0.064,0.068,0.075和0.111,接觸載荷依次為6.42×10-4,8.05×10-4,8.39×10-4和1.13×10-3。本文模型的接觸曲線均與ZMC模型和KE模型較為接近,而CEB模型的曲線明顯偏大。這一現(xiàn)象的原因是本文模型、ZMC模型和KE模型都考慮了微凸體彈塑性階段,但與其他模型的差異在于本文模型采用二次函數(shù)回轉(zhuǎn)體模擬微凸體輪廓,建立了更加直觀的微凸體接觸半徑與接觸變形的解析關(guān)系。結(jié)合3.1節(jié)和3.2節(jié)的對比驗證,可進(jìn)一步說明了本文模型的有效性。

圖7 不同粗糙表面接觸模型對比Fig.7 Comparison of different rough surface contact models

4 討論與分析

4.1 塑性指數(shù)對接觸面積的影響

根據(jù)表1的表面參數(shù),繪制無量綱接觸面積A*與無量綱接觸載荷F*的關(guān)系,如圖8所示。不同塑性指數(shù)對應(yīng)的表面參數(shù)分別為Ψ1=0.7,β1=0.033 9;Ψ2=1.5,β2=0.046 7;Ψ1=2.5,β1=0.060 1。

圖8 不同塑性指數(shù)下接觸面積與接觸載荷的關(guān)系Fig.8 Relationship between contact area and contact load under different plasticity indices

由圖8可知:①隨著接觸載荷的增大,不同塑性指數(shù)對應(yīng)的接觸面積也隨之增大,呈現(xiàn)出近線性的關(guān)系,這是因為載荷增大,兩表面的接觸間隙減小,更多的微凸體受到擠壓開始接觸,增大了接觸面積;②接觸載荷越大,塑性指數(shù)對接觸面積的影響越大;同等載荷條件下,塑性指數(shù)越大,接觸面積越小。這是由于當(dāng)塑性指數(shù)較大時,對應(yīng)的表面粗糙度越大,雖然會有更多的接觸微凸體進(jìn)入塑性屈服狀態(tài),但單位面積上發(fā)生接觸的微凸體數(shù)目較少,其接觸面積也就越小。

4.2 塑性指數(shù)對接觸剛度的影響

圖9為無量綱接觸剛度K*與無量綱接觸載荷F*的關(guān)系。由圖9可知:①不同塑性指數(shù)對應(yīng)的接觸剛度均隨著接觸載荷的增加而增大,相同載荷條件下,塑性指數(shù)越大,接觸剛度越小,這是由于塑性指數(shù)越大,將會有更多的接觸微凸體進(jìn)入塑性屈服狀態(tài),表面抵抗彈性變形能力變?nèi)?,?dǎo)致接觸剛度降低;②接觸載荷越大,塑性指數(shù)對接觸剛度的影響越大。在接觸載荷F*=1×10-4時,三種塑性指數(shù)對應(yīng)的接觸剛度分別為3.09,0.48,0.11,兩兩之間的差值分別為2.61,0.37。這說明,接觸載荷較大時,不同塑性指數(shù)的接觸剛度差異較大,塑性指數(shù)是影響接觸剛度的主要因素。

圖9 不同塑性指數(shù)下接觸剛度與接觸載荷的關(guān)系Fig.9 Relationship between contact stiffness and contact load under different plasticity indices

4.3 微凸體尺寸參數(shù)對接觸剛度的影響

為研究微凸體尺寸參數(shù)對接觸剛度的影響,選取圖2中的粗糙表面(Ra2=1.302 μm),測量位置如圖1(b)所示的三條線段,獲得三組微凸體尺寸參數(shù):l1=15.273 μm,h1=1.039 μm;l2=16.790 μm,h2=1.082 μm;l3=17.985 μm,h3=1.065 μm。繪制不同微凸體尺寸參數(shù)的無量綱接觸剛度K*與無量綱接觸載荷F*的關(guān)系,如圖10所示。由圖10可知:①隨著接觸載荷的增加,接觸間隙開始減小,相互接觸的微凸體數(shù)目開始增多,接觸面積變大,接觸剛度亦隨之增大;②同等載荷下,微凸體直徑與高度的比值越大,接觸剛度越大,這是由于微凸體直徑與高度的比值越大,接觸半徑越大,從而產(chǎn)生更大的接觸面積,粗糙表面抵抗彈性變形的能力就越強(qiáng);③接觸載荷越大,微凸體尺寸參數(shù)對接觸剛度的影響較大。在接觸載荷F*=1×10-4時,三種微凸體尺寸參數(shù)對應(yīng)的接觸剛度分別為0.25,0.21,0.18,兩兩之間的差值分別為0.04,0.03。與圖10的剛度差值對比,進(jìn)一步表明塑性指數(shù)對接觸剛度的影響更大。

圖10 不同微凸體尺寸參數(shù)下接觸剛度與接觸載荷的關(guān)系Fig.10 Relationship between contact stiffness and contact load under different asperity size parameters

5 結(jié) 論

(1) 提出一種采用二次函數(shù)回轉(zhuǎn)體等效微凸體的建模方法,實現(xiàn)了對微凸體整體輪廓的擬合,建立了接觸半徑與接觸變形的解析關(guān)系,彌補(bǔ)了半球體模型中單一曲率半徑的缺陷。

(2) 根據(jù)二次函數(shù)模型,重新推導(dǎo)出微凸體在彈性、彈塑性和完全塑性的接觸表達(dá)式,并應(yīng)用接觸力學(xué)理論和概率統(tǒng)計方法,建立了粗糙表面的微觀接觸模型。

(3) 對比分析了CEB模型、ZMC模型和KE模型與本文模型的差異。研究發(fā)現(xiàn),本文模型的預(yù)測值與ZMC模型和KE模型較為接近,而CEB模型明顯偏大;本文模型實現(xiàn)了對微凸體整體輪廓的擬合,更適用于重載情況下進(jìn)行接觸特性分析。

(4) 塑性指數(shù)是影響接觸剛度的主要因素,塑性指數(shù)越大,其接觸面積越小,抵抗變形的能力越弱;微凸體尺寸參數(shù)對接觸剛度的影響較小,微凸體直徑與高度的比值越大,接觸面積和接觸剛度越大。

猜你喜歡
變形模型
一半模型
重要模型『一線三等角』
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
“我”的變形計
變形巧算
例談拼圖與整式變形
會變形的餅
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 麻豆国产原创视频在线播放| 亚洲天堂777| 国产丝袜第一页| 欧美国产日本高清不卡| 国产呦精品一区二区三区下载| 伊人久久婷婷五月综合97色| 国产成人高清精品免费| 精品一区二区三区水蜜桃| 久久99国产综合精品1| 国产经典免费播放视频| 中文字幕在线免费看| 国产黄色视频综合| 亚洲青涩在线| 国产迷奸在线看| 国产精品白浆无码流出在线看| 色成人综合| 国产真实二区一区在线亚洲| 香蕉综合在线视频91| 免费亚洲成人| 久久网欧美| 91精品综合| 色婷婷电影网| 在线精品自拍| 在线观看国产黄色| 国产精品久久久精品三级| 99国产精品一区二区| 日韩欧美色综合| av在线人妻熟妇| 多人乱p欧美在线观看| 欧美狠狠干| 国产91av在线| 亚洲精品国产精品乱码不卞| 天堂网亚洲系列亚洲系列| 亚洲成人精品久久| 99无码熟妇丰满人妻啪啪 | 亚洲欧洲自拍拍偷午夜色无码| 亚洲欧美日韩精品专区| 亚洲欧美国产视频| 欧美激情,国产精品| 91精品国产综合久久香蕉922| 777国产精品永久免费观看| 日本伊人色综合网| 色婷婷国产精品视频| 久久精品国产999大香线焦| 东京热av无码电影一区二区| 毛片最新网址| 亚洲第一天堂无码专区| 在线观看国产黄色| 亚洲综合婷婷激情| 国产靠逼视频| 夜夜爽免费视频| 日本一区二区三区精品视频| 三上悠亚一区二区| 国产永久无码观看在线| 国产精品第一区在线观看| 成人免费午间影院在线观看| 青青操国产视频| 91激情视频| 久久久精品无码一区二区三区| 日韩黄色在线| 91久久偷偷做嫩草影院精品| 亚洲无线国产观看| 中国一级毛片免费观看| 欧美第一页在线| 国产91在线|日本| 国产精品成人啪精品视频| 国内老司机精品视频在线播出| 精品国产免费观看一区| 日韩精品一区二区三区大桥未久| 久久久久人妻一区精品色奶水| 欧美精品v| 国产乱子伦视频在线播放| 欧亚日韩Av| 亚洲AV免费一区二区三区| 国产熟睡乱子伦视频网站| 黄色片中文字幕| 999国产精品永久免费视频精品久久| 国产精品无码久久久久久| 亚洲综合一区国产精品| 狠狠躁天天躁夜夜躁婷婷| 久草视频福利在线观看| 99精品高清在线播放|