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

        ?

        模擬仿真技術(shù)在磨損預(yù)測中的應(yīng)用現(xiàn)狀

        2013-04-17 11:44:48熊壘林有希
        機械制造與自動化 2013年5期
        關(guān)鍵詞:模型研究

        熊壘,林有希

        (福州大學(xué)機械工程及自動化學(xué)院,福建福州 350108)

        0 引言

        磨損作為摩擦學(xué)的一個重要方面,許多學(xué)者已經(jīng)進行了非常多的理論和實驗研究[1-2]。目前關(guān)于磨損的研究手段主要以實驗為主,但是由于摩擦學(xué)系統(tǒng)本身龐大而復(fù)雜,其中涉及的不可控制變量太多,相關(guān)實驗需要耗費大量的人力物力,而且實驗的重復(fù)性和可比性較弱,這給定量研究磨損行為造成了很大困難[3],所以有部分學(xué)者開始嘗試運用模擬仿真的手段來進行磨損的研究[4]。模擬仿真的優(yōu)點在于其變量易于控制,針對一種摩擦磨損現(xiàn)象,可以刻意的選擇相關(guān)參數(shù)和約束,從而避免其他因素的干擾。同時模擬仿真的重復(fù)性強,可比性好,更有利于有關(guān)規(guī)律和機理的深入研究[5]。本文綜述了模擬仿真技術(shù)在磨損研究中的應(yīng)用現(xiàn)狀,并就當(dāng)前存在的問題提出了建議。

        1 兩種主要磨損的產(chǎn)生機理

        1.1 粘著磨損

        當(dāng)聚合物在金屬表面滑動時產(chǎn)生粘著,粘著磨損的重要特征是軟的聚合物向硬質(zhì)金屬表面的轉(zhuǎn)移:開始時聚合物轉(zhuǎn)移到金屬上,接著聚合物本體與附著在金屬表面的聚合物轉(zhuǎn)移膜之間產(chǎn)生粘著作用,然后在剪切力的作用下,粘著點處的聚合物從本體上撕扯下來并積累在金屬表面,逐漸形成轉(zhuǎn)移膜,隨著轉(zhuǎn)移膜不斷變厚,同時在循環(huán)載荷的反復(fù)作用下,轉(zhuǎn)移膜破裂,一部分碎片脫離摩擦表面,一部分加入其他的磨損形式中,這樣就形成材料逐漸損耗的過程[6]。轉(zhuǎn)移膜的生成改變了摩擦接觸的類型,從金屬與聚合物的摩擦接觸變成聚合物與聚合物的摩擦接觸;同時也改變了金屬表面及形貌。一般來說較軟的物質(zhì)易在較硬的表面形成粘著磨損,如聚合物/金屬摩擦副;物理性質(zhì)相近或者相同者也易形成粘著磨損,如金屬/金屬摩擦副[7]。

        1.2 磨粒磨損

        外界硬質(zhì)顆粒或者對磨表面上的硬凸起物或粗糙峰在摩擦過程中引起表面材料脫落的現(xiàn)象稱為磨粒磨損。磨粒磨損很少是一種單一的磨損機制引起的,而經(jīng)常是多種磨損機制綜合作用的結(jié)果。而且磨損機制通常隨著工況條件的變化而變化。磨粒磨損過程中材料的去除機制主要分為以下兩種[8]:

        a)由塑性變形機制引起的去除過程

        1)犁溝效應(yīng)。在載荷作用下磨粒被壓入材料表面,同時向前運動,受到擠壓的材料產(chǎn)生塑性流動,向兩側(cè)隆起。這種過程不會直接導(dǎo)致材料的去除,但是在磨粒的反復(fù)作用下,材料產(chǎn)生多次變形,最后便會脫落,形成磨損。

        2)微觀切削。在載荷作用下,磨粒在材料表面發(fā)生類似刨削一樣的過程,這種過程會直接導(dǎo)致材料的去除。

        b)由斷裂機制引起的去除過程

        這種磨損機制主要針對脆性材料而言,比如陶瓷、玻璃等。磨粒在壓力作用下會與接觸表面形成極高的接觸應(yīng)力,以致于超過脆性材料的極限,從而發(fā)生脆性斷裂。

        2 磨粒磨損的模擬仿真

        根據(jù)磨粒磨損的概念可知,表面粗糙度對磨粒磨損的影響較大。由于表面形貌隨機性強、形成過程復(fù)雜、影響因素多,所以不可能通過實時監(jiān)測的手段來觀察表面形貌的變化,現(xiàn)在對表面形貌的研究多采用模擬仿真的手段[9]。趙永春[10]采用離散元原理直觀的對表面粗糙度進行了模擬,結(jié)果表明在剛性光滑平面與理想粗糙表面相互作用后,粗糙度變小,表面結(jié)構(gòu)也發(fā)生很大變化,隨著兩個接觸表面法向接近量的增大,粗糙表面脫落的顆粒數(shù)增加,磨損加劇。Zhang[11-12]基于隨機過程理論提出了描述二體磨損過程中表面形貌變化的磨損動態(tài)過程理論,通過假設(shè)剛性表面為絕對剛體,建立了軟表面與剛性表面摩擦?xí)r,受微觀切削而引起的表面粗糙度變化和磨損量的預(yù)測模型。劉峰壁[13]考慮到實際過程中硬表面也會受到磨損,于是在此模型基礎(chǔ)上進行修正,建立了兩表面均受到磨損的預(yù)測模型。馮偉等[14]利用表面輪廓曲線高度近似服從高斯分布這一統(tǒng)計學(xué)特點,分別用matlab和ansys建立了表面形貌的三維仿真模型,并用赫茲接觸模型計算了不同載荷下接觸應(yīng)力在表面的分布情況。

        吳國清等[15]運用Monte Carlo方法隨機創(chuàng)建磨粒形狀,采用打靶法隨機生成磨粒表面,采用網(wǎng)格剖分法記憶表面形貌,在三維空間內(nèi)模擬了多個磨粒同時參與磨損的隨機動態(tài)磨損過程,并通過實驗驗證了模型的有效性。岑巖宏[16]利用三維有限元法動態(tài)模擬剛性球形磨粒在彈-塑性平板上的壓入-滑動過程,模擬了不同壓入深度下表面的磨損劃傷行為,結(jié)合實驗數(shù)據(jù)建立了適合碳鋼系列材料塑性變形脊形貌與磨粒壓入深度和寬度的曲線關(guān)系。Jacobson[17]建立了兩維磨損面內(nèi)僅限于純切削過程的統(tǒng)計學(xué)模型,在模型中考慮了多個磨粒同時作用的影響,從而使它更接近于實際情況。文中利用該模型研究了磨粒尺寸、載荷、工件表面硬度對磨損率的影響,預(yù)測了接觸的磨粒數(shù)量和磨損表面形貌,并與其他研究試驗結(jié)果進行對比,驗證了模型的有效性。Podra[18]建立錐形接觸時考慮表面形貌的磨粒磨損模型,用有限元和解析方法對磨損進行了計算,但兩種方法計算結(jié)果相差較大,分析認(rèn)為產(chǎn)生差異的原因是在解析方法中沒有考慮彈性變形的影響。M.Barge[19]等著重研究了磨粒磨損中的犁溝現(xiàn)象,建立了較硬的圓頭形磨粒在較軟表面的犁溝切削模型,考察了載荷、材料和磨粒的硬度等因素對表面應(yīng)力分布、殘余應(yīng)力和磨損量的影響,并與Bowden和Tabor的犁溝切削模型進行了對比。Anders S?derberg[20]通 過 對 磨 損 過 程 的 簡 化,基 于Archard磨損模型和顯式歐拉積分求解方法,建立了制動片在磨合期時的三維有限元模型,分析了磨合過程中接觸應(yīng)力的分布和磨損變化,發(fā)現(xiàn)在磨合初期,表面粗糙度較大,應(yīng)力分布非常不均勻,隨著磨合的進行,應(yīng)力分布趨于均勻,磨損程度和應(yīng)力分布有關(guān),即接觸應(yīng)力大的區(qū)域磨損較為嚴(yán)重。

        3 粘著磨損的模擬仿真

        粘著磨損中較有代表性的是Holm模型和Archard模型[21],相比Holm模型,Archard模型更接近于粘著磨損的真實情況,因此大部分粘著磨損都是在此模型基礎(chǔ)上加以修正而成。

        Sung-San Cho[22]等從分子作用力的角度出發(fā),建立了半球形彈性體與剛性表面的粘著接觸有限元模型,計算了載荷、接觸半徑、脫粘力等相關(guān)參數(shù)的關(guān)系,結(jié)果表明較小的表面粗糙度能有效減小靜摩擦力和粘著摩擦力,但是會加劇粘著磨損。X.Yin和K.Komvopoulos[23]建立了粗糙表面粘著效應(yīng)的粘著磨損模型,以表面微凸體是否發(fā)生完全塑性變形為磨損發(fā)生準(zhǔn)則,將接觸面積、磨損率、磨損系數(shù)等目標(biāo)參數(shù)表示為法向載荷、表面粗糙度、材料的彈塑性變形參數(shù)、表面能、材料相容度等因素的函數(shù),并采用典型的陶瓷/陶瓷、陶瓷/金屬、金屬/金屬摩擦副的材料參數(shù)對模型進行計算,得到了相關(guān)影響因素與目標(biāo)參數(shù)之間的曲線關(guān)系,該模型比較完整的包含了粘著磨損過程中的各影響因素,對磨損的建模過程具有一定的借鑒意義。

        綜上所述,數(shù)值仿真方法作為新的磨損研究手段和方法己引起摩擦學(xué)界的重視,并將發(fā)揮越來越重要的作用。目前關(guān)于磨損仿真研究的主要不足在于磨損仿真模型預(yù)測結(jié)果同實際情況相差較大;建立仿真模型時的簡化處理導(dǎo)致模型精確度不夠高;利用Arhcard等經(jīng)典磨損模型進行仿真計算時缺乏關(guān)于材料磨損率規(guī)律的準(zhǔn)確數(shù)據(jù)[24]。磨損仿真結(jié)果的可靠性和實用性在相當(dāng)大程度上取決于磨損模型的正確性和精確性。為了更好地發(fā)揮數(shù)值仿真方法在摩擦磨損研究中的作用,應(yīng)在深化磨損機理研究、摩擦學(xué)系統(tǒng)分析方法及磨損表面微觀分析的基礎(chǔ)上建立更加完善和切合實際的磨損模型,并通過大量實驗研究確定關(guān)于磨損率規(guī)律的準(zhǔn)確數(shù)據(jù)。

        4 溫度-應(yīng)力場仿真

        在載荷作用下,摩擦接觸區(qū)域會產(chǎn)生大量的摩擦熱,摩擦熱嚴(yán)重影響材料的物理力學(xué)性能和化學(xué)性能,對磨損形式的改變有著直接影響。例如制動器在制動過程中溫度場分布不均,這將引起兩接觸體的熱彈性變形,從而改變了接觸面之間的壓力分布。反過來,接觸界面之間的壓力分布不均勻,使制動器溫度場分布更加不均勻,形成了一種惡性的循環(huán),這種現(xiàn)象被稱作熱彈性不穩(wěn)定(TEI)[25]。對樹脂基制動片/金屬摩擦副而言,摩擦表面溫度分布不均勻性會導(dǎo)致制動片一系列物化性能的變化,如熱彈性不穩(wěn)定性、熱分解、熱膨脹等。由于摩擦過程的特殊性,摩擦表面的溫度和應(yīng)力大小、變化規(guī)律、分布情況等都不可能實時測量,這為進一步的研究工作造成困難,模擬仿真技術(shù)使其成為可能。

        4.1 溫度場

        摩擦副表面溫度的計算模型普遍采用Block溫度模型,即半無限體表面承受單一集中熱源,并提出了閃溫的概念。接著Jaeger發(fā)展了這一理論,闡述了矩形形狀的移動熱源作用在半無限體表面上的數(shù)學(xué)模型[26]。Ling[27]運用表面形貌的隨機模型來估算滑動表面的瞬時溫升,Wang S[28]則對具有分形特征的粗糙表面的滑動摩擦局部溫升進行了研究,提出了界面溫度分布的分形理論。Tian和Kennedy[29]的研究指出,實際的滑動接觸都是多點離散接觸且滑動副物體尺寸都是有限的,在接觸微凸體處除了有局部閃現(xiàn)溫升作用外,還有名義表面溫升的作用。這要求考慮物體的有限形狀,滑動接觸的重復(fù)特性以及物體間總的“熱量平衡”。

        由上述可知,制動摩擦溫度場問題是典型的非線性問題,滑動摩擦表面的溫升由名義表面溫升和局部表面溫升(閃現(xiàn)溫度)組成。求解名義表面溫升一般都假設(shè)熱流分配系數(shù)為一常數(shù),然后按照理論分析(或某個實驗數(shù)據(jù))得出一個確定的熱流分配系數(shù)。實際上,制動過程中摩擦表面上產(chǎn)生的熱量在零件間的分配是變化的,與物體的尺寸大小,材料的熱物理性能及制動時間有關(guān),并隨材料的導(dǎo)熱性與接觸表面處的溫度梯度之比值而變化[30]。

        在制動摩擦名義表面溫升計算的熱模型中,名義接觸表面摩擦熱流邊界條件的確定是建立和求解的關(guān)鍵。K.Fridrich[31]建立了單個鋼微凸體在PEEK或碳纖維-PEEK復(fù)合材料上滑動的模型,該模型假設(shè)單位面積上所產(chǎn)生的摩擦熱取決于接觸應(yīng)力p、摩擦系數(shù)μ和滑動速度v,然后通過計算接觸點處的單位熱源強度來計算接觸溫度。Kennedy和Ling[32]最早運用數(shù)值模擬的方法研究盤式制動器的熱彈性不穩(wěn)定現(xiàn)象,但他的模型簡化為軸對稱模型,忽略了移動熱源的問題。

        4.2 應(yīng)力場

        摩擦過程中溫度場和應(yīng)力場是相互耦合、彼此影響的。M.Tirovic和A.J.Day[33]在詳細(xì)地研究了盤式制動器的壓力分布情況后,認(rèn)為摩擦制動器的摩擦襯片與對偶件表面之間產(chǎn)生的摩擦熱并不是均勻地分布在滑動表面上的,而是取決于局部應(yīng)力的分布情況;而由于溫度效應(yīng)造成熱膨脹和表面熱裂紋等結(jié)果又會反過來影響應(yīng)力分布。唐旭晟在求解溫度-應(yīng)力耦合場的做法是:首先假設(shè)初始階段應(yīng)力是均勻分布的,求出溫度場的解,然后將溫度作為施加載荷,反過來求解溫度引起的熱應(yīng)力分布。Hasan Sofuoglu等[34]分別求解了考慮溫度和不考慮溫度影響下的表面應(yīng)力分布,發(fā)現(xiàn)在溫度影響下表面應(yīng)力更大,而較大的表面應(yīng)力主要分布在第一粗糙峰上,這是因為在第一次法相接觸時,主要是這些較“突出”的微凸體承受載荷,當(dāng)施加位移時,摩擦熱在這些微凸體上造成溫度累積效應(yīng),從而產(chǎn)生熱彈塑性變形,所以表面應(yīng)力要更大一些。Gong和Komvopoulos[35]也分別計算了只考慮彈性接觸/彈塑性接觸和溫度影響下的熱彈塑性接觸的應(yīng)力大小,結(jié)果與 Hasan Sofuoglu 一致。Ji-Hoon Choi[36]建立了三維有限元盤式制動器模型,模擬了反復(fù)制動情況下盤/片的熱彈性接觸問題,深入分析了溫度、應(yīng)力及材料屬性之間的內(nèi)在關(guān)系。發(fā)現(xiàn)在初始階段應(yīng)力分布較均勻,隨著制動的進行應(yīng)力沿徑向方向發(fā)生明顯變化,局部出現(xiàn)應(yīng)力集中,應(yīng)力的變化主要是由于溫度升高導(dǎo)致表面變形而造成的,應(yīng)力分布和大小隨著溫度變化而變化,即所謂的熱應(yīng)力不穩(wěn)定性(TEI)。在研究材料屬性與應(yīng)力分布關(guān)系時,發(fā)現(xiàn)熱膨脹系數(shù)和彈性模量對應(yīng)力分布有很大影響,在仿真中將彈性模量減半后,應(yīng)力分布較之前更加均勻,這說明較軟的材料有利于增加實際接觸面積,使應(yīng)力分布更加均勻。相比各向同性材料,各向異性材料的熱彈性穩(wěn)定性更好,溫度和應(yīng)力的分布也更加均勻。

        5 結(jié)語

        如前所述,由于模擬仿真技術(shù)更易于變量的控制,與試驗研究相比,具有可操作性、可重復(fù)性、可比性等優(yōu)點,因而越來越受到摩擦磨損領(lǐng)域的研究學(xué)者們的重視。但是也存在一些問題,如模擬結(jié)果與實驗數(shù)據(jù)不吻合,建模時的簡化處理影響模型精度,個別物理量難以通過實驗得到驗證等。為了更好的使用模擬仿真技術(shù)進行研究,筆者覺得應(yīng)該從以下方面入手:首先,深化磨損機理的研究。只有真正弄清楚磨損機理,建立完善的磨損理論,才能從根本上指導(dǎo)模擬仿真的建模過程,才能保證模型的科學(xué)性和精確性。其次,盡量在模型中考慮更多的參量,過度的簡化會使仿真精度下降;再次,加強實驗驗證的環(huán)節(jié),針對所建立的模型設(shè)計專門的實驗,只有通過實驗驗證才能最終保證模型的正確性。

        [1]溫詩鑄.摩擦學(xué)原理[M].北京:清華大學(xué)出版社,1990.

        [2]黃平.摩擦學(xué)教程[M].北京:高等教育出版社,2008.

        [3]江親瑜.零件磨損過程及壽命預(yù)測的數(shù)值仿真[J].潤滑與密封,1997,29(6):29-30.

        [4]Yan L,XuJ J.Research on the simulation method of wear process[J].Tribology,1999,19(1):50-55.

        [5]趙春華.基于摩擦學(xué)與動力學(xué)的摩擦學(xué)系統(tǒng)狀態(tài)描述方法研究[D].武漢,武漢理工大學(xué),2005.

        [6]B.布尚.摩擦學(xué)導(dǎo)論[M].北京:機械工業(yè)出版社,2006.

        [7]蘇立志,李志章.耐溫聚合物基復(fù)合材料的摩擦學(xué)性能與磨損機理研究進展[J].材料科學(xué)與工程,2000(70):125-130.

        [8]劉家浚.材料磨損機理及其耐磨性[M].北京:清華大學(xué)出版社,1993.

        [9]閆小青,謝志龍,樊保圣,等.基于表面粗糙度影響的摩擦聲發(fā)射特性研究[J].Tribology,2011(31):587-591.

        [10]趙永春,劉小君,王偉,等.基于離散元法的表面形貌相互作用研究[J].合肥工業(yè)大學(xué)學(xué)報,2010(33):487-490.

        [11]Zhang B,Xie Y B.Two body micro-cutting wear model,Part I:Two-dimensional roughness model[J].Wea,1989,129:37-48.

        [12]Zhang B,Xie Y B.Two body micro-cutting wear model,Part II:Three-dimensional roughness model[J].Wear,1989,129:49-58.

        [13]Liu F B,Li X E,Xie Y B.A new model for predicting roughness of rubbing surface in three bodyabrasive wear[J].Chinese Journal of Mechanical Engineering,2000,19(1):9-12.

        [14]馮偉,嚴(yán)新平,周新聰.銷盤滑動磨損試驗的仿真建模研究[J].中國機械工程,2005(16):2141-2144.

        [15]吳國清,張曉峰,方亮,等.兩體磨料磨損的三維動態(tài)模擬[J].Tribology,2000(20):360-364.

        [16]岑啟宏,孫錕,方亮,等.二體磨料磨損犁溝及脊的三維有限元動態(tài)模擬[J].Tribology,2004(24):249-253.

        [17]Jacobson S,Wallen P,Hogmark S.Fundamental aspects of abrasive wear studied by a new numerical simulation model[J].Wear,2008,123:207-223.

        [18] Podra P,Andersson S.Finite element analysis of wear simulation of a conical spinning contact considering surface topography[J].Wear,1999,224:13-21.

        [19]M.Barge,G.Kermouche,P.Gilles,J.M.Bergheau.Experimental and numerical study of the ploughing part of abrasive wear[J].Wear,2003(255):30-37.

        [20]Anders Sderberg,S ren Andersson.Simulation of wear and contact pressure distribution at the pad-to-rotor interface in a disc brake using general purpose nite element analysis software[J].Wear,2009(267):2243-2251.

        [21]Valentin L.Popov.Contact Mechansics and Friction physical Principles and Applications[M].北京:清華大學(xué)出版社,2011.

        [22]Sung-San Cho,Seungho Park.Finite element modeling of adhesive contact using molecular potential[J].Tribology,2004(37):763-769.

        [23]X.Yin,K.Komvopoulos.An adhesive wear model of fractal surfaces in normal contact[J].Solid and Structure,2010(47):912-921.

        [24]汪選國,嚴(yán)新平,李濤生,等.磨損數(shù)值仿真技術(shù)的研究進展[J].摩擦學(xué)學(xué)報,2004(24):188-192.

        [25]高誠輝,黃健萌,林謝昭,等.盤式制動器摩擦磨損熱動力學(xué)研究進展[J].中國工程機械學(xué)報,2006(4):83-88.

        [26] O.O.Evtushenko,E.H.Ivanyk,N.V.Horbachova.Analytic methods for thermal calculation of brakes(review)[J].Materials Science,2000,36(6):857-862.

        [27]F.Ling.On temperature transients at sliding interface[J].ASME Journal of Lubrication Technology,1969,7:397-405.

        [28] Wang S.A fractal theory of the interfacial temperature distribution in the slow sliding Regime[J].ASME J of Trio.,1994,116-812.

        [29]Xuefeng Tian,F(xiàn)rancis E Kennedy Jr.Contact surface temperature models for finite bodies in dry and boundary lubricated sliding[J].ASME Journal of Tribology,1993,115:411-418.

        [30]馬保吉,朱均.摩擦制動器接觸表面溫度計算模型[J].西安工業(yè)學(xué)院學(xué)報,1999,19(1):35-39.

        [31]K.Friedrich,J.Flock,K.Varadi and Z.Neder.Numerical and finite element contact and thermal analysis of real composite-steel surfaces in sliding contact[J].Wear,1999,225-229:368-379.

        [32]Kennedy,F(xiàn).E.,Ling,F(xiàn).F.A thermal,Thermoelastic,and Wear Simulation of a High-Emergy SlidingContact Problem[J].ASME J.Lubr.Technol,1974,l97,497-508.

        [33]J.Halling.Principles of Tribology[M].北京:機械工業(yè)出版社,1975.

        [34]Hasan Sofuoglu,Alaettin Ozer.Thermomechanical analysis of elastoplastic medium in sliding contact with fractal surface[M].Tribology,2008(41):783-796.

        [35] Gong Z-Q,Komvopoulos K.Thermo mechanical analysis of semi-infinite solid in sliding contact with a fractal surface.ASME J Tribol,2005;127:331 – 42.

        [36]Ji-Hoon Choi,In Lee.Finite element analysis of transient thermoelastic behaviors in disk brakes[J].Wear,2004(257):47-58.

        猜你喜歡
        模型研究
        一半模型
        FMS與YBT相關(guān)性的實證研究
        2020年國內(nèi)翻譯研究述評
        遼代千人邑研究述論
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        視錯覺在平面設(shè)計中的應(yīng)用與研究
        科技傳播(2019年22期)2020-01-14 03:06:54
        EMA伺服控制系統(tǒng)研究
        新版C-NCAP側(cè)面碰撞假人損傷研究
        3D打印中的模型分割與打包
        妃光莉中文字幕一区二区| 在线观看国产一区亚洲bd| 国产自产拍精品视频免费看| 91亚洲免费在线观看视频| 欧美丰满老熟妇aaaa片| 亚洲18色成人网站www| 国产日韩欧美911在线观看| 国产精品自产拍av在线| 伊人久久大香线蕉av不变影院| 在线视频观看一区二区| 久久综合香蕉国产蜜臀av| 成年午夜无码av片在线观看| 国产av一区二区三区丝袜| 日本女同视频一区二区三区| 国产精品久久久爽爽爽麻豆色哟哟| 天天躁狠狠躁狠狠躁夜夜躁| 亚洲日产无码中文字幕| 性色国产成人久久久精品二区三区 | 成年女人午夜特黄特色毛片免| 久久精品女人av一区二区| 中文字幕人妻少妇引诱隔壁| 国产精品露脸张开双腿| 无码日日模日日碰夜夜爽| 国内精品久久久久久久久齐齐| 亚洲精品美女久久久久99| 99精品国产综合久久麻豆| 一本色道久久88综合日韩精品| 久久精品免费一区二区喷潮| 国产精品久久中文字幕亚洲| 美腿丝袜诱惑一区二区| 人妻少妇精品专区性色av| 亚洲国产成人资源在线桃色| 亚洲精品国产亚洲av| 啦啦啦中文在线观看日本| 日韩激情小视频| 日本久久一级二级三级| 亚洲国产精品成人久久久| 久久精品国产亚洲av高清漫画| 最新福利姬在线视频国产观看| 美女视频一区二区三区在线| 亚洲男人av天堂午夜在|