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

含夾雜物軸承鋼中裂紋的萌生與擴(kuò)展

2022-06-08 05:06:24羅敏汪久根馮毅雄馮照和
軸承 2022年2期
關(guān)鍵詞:裂紋有限元模型

羅敏,汪久根,馮毅雄,馮照和

(1.浙江大學(xué) 機(jī)械工程學(xué)院,杭州 310027;2.杭州汽輪機(jī)股份有限公司,杭州 310022)

0 引言

軸承作為機(jī)械旋轉(zhuǎn)體的支承結(jié)構(gòu),是工業(yè)機(jī)械中應(yīng)用廣泛的部件,研究滾動軸承的裂紋萌生、擴(kuò)展方式及疲勞壽命對保證其可靠性至關(guān)重要。在軸承滾動接觸表面光滑且潤滑條件適合的情況下,軸承次表層微裂紋擴(kuò)展引發(fā)的層裂在滾動接觸疲勞失效中占主導(dǎo)地位。次表層的裂紋萌生往往始于材料缺陷,特別是夾雜物引發(fā)的應(yīng)力集中現(xiàn)象促使裂紋萌生與擴(kuò)展。因此,圍繞軸承次表層分布的夾雜物展開研究,建立軸承滾子與滾道接觸區(qū)仿真模型進(jìn)行模擬具有一定的工程意義。

分析應(yīng)力對接觸疲勞的貢獻(xiàn),經(jīng)歷了靜態(tài)應(yīng)力、最大剪應(yīng)力、最大正交剪應(yīng)力與Mises應(yīng)力的發(fā)展過程,目前常用Mises應(yīng)力分析滾動接觸疲勞失效。裂紋的萌生與擴(kuò)展過程十分復(fù)雜且具有隨機(jī)性,需用統(tǒng)計分析方法研究。軸承滾動體與內(nèi)圈滾道的滾動接觸過程中發(fā)生疲勞損傷累積,會在次表層萌生微裂紋,多個微裂紋的聯(lián)合形成了向表面擴(kuò)展的較長裂紋,裂紋到達(dá)表面產(chǎn)生剝落。文獻(xiàn)[1]綜述了非金屬夾雜物對滾動軸承疲勞壽命的影響,介紹了麻點(diǎn)和剝落2種疲勞失效類型的理論分析方法和試驗(yàn)測試技術(shù)。

關(guān)于接觸疲勞的研究已有80多年的歷史,提出了多種方法分析夾雜物對疲勞過程的影響:文獻(xiàn)[2]提出用最大能量釋放判據(jù)來分析夾雜物邊界上的裂紋萌生過程,建立了彈塑性多軸疲勞裂紋萌生有限元模型,對滾動接觸疲勞裂紋的萌生壽命進(jìn)行預(yù)測;文獻(xiàn)[3]用Eshelby張量分析夾雜物的影響;文獻(xiàn)[4]用格林函數(shù)分析夾雜物對材料力學(xué)性能的影響。

近年來,文獻(xiàn)[5-15]較系統(tǒng)地分析了夾雜物對接觸疲勞的影響:文獻(xiàn)[5-6]用計算幾何的方法分析了材料的拓?fù)潆S機(jī)性和材料性能的隨機(jī)性對韋布爾分布斜率的影響,分布斜率1.29~3.36時,疲勞壽命與最大赫茲接觸應(yīng)力的9.35次方成反比;文獻(xiàn)[7-8]分析了晶粒間的疲勞損傷累積,認(rèn)為晶粒的邊界對剝落失效與疲勞壽命的離散性有顯著的影響;文獻(xiàn)[9]建立三維有限元模型分析微觀結(jié)構(gòu)對滾動接觸疲勞的影響,這與球軸承的接觸問題一致;文獻(xiàn)[10-11]用彈塑性有限元模型結(jié)合碳彌散模型分析了白層的產(chǎn)生及其方向,模擬了馬氏體的退化,認(rèn)為在應(yīng)力集中區(qū)域中萌生裂紋;文獻(xiàn)[12-13]模擬分析了夾雜物周圍蝴蝶狀裂紋的形成、萌生與擴(kuò)展過程,與試驗(yàn)結(jié)果相吻合;文獻(xiàn)[14] 提出一種結(jié)合損傷力學(xué)的內(nèi)聚力模型,用于模擬多晶軸承鋼在滾動接觸循環(huán)加載下的裂紋萌生與擴(kuò)展,模型中的損傷在晶界由內(nèi)聚力元素受損斷裂引發(fā),最終獲得的層裂方式和疲勞壽命分布與公開文獻(xiàn)中的結(jié)果相符;文獻(xiàn)[15]建立了基于連續(xù)損傷力學(xué)的有限元模型,研究深溝球軸承的滾動接觸疲勞壽命,建立了相關(guān)壽命方程。這些工作為揭示滾動疲勞機(jī)理提供了新的認(rèn)識,然而,針對軸承鋼疲勞壽命的準(zhǔn)確計算仍需進(jìn)一步研究。

文獻(xiàn)[16]分析了非金屬夾雜物深度、徑向載荷、表面摩阻力和滾道曲率半徑對表面裂紋萌生的作用,徑向載荷和表面摩阻力對裂紋萌生有顯著的影響,而滾道曲率半徑對裂紋擴(kuò)展的方向沒有影響。文獻(xiàn)[17]試驗(yàn)分析了裂紋的萌生與擴(kuò)展階段,認(rèn)為在早期就萌生了裂紋,并且主要是裂紋擴(kuò)展階段決定了疲勞壽命的長短。滾動接觸疲勞不僅是滾動軸承的重要問題,對于列車的輪軌接觸、摩擦無級變速器的壽命也是重要問題。文獻(xiàn)[18] 研究了溫度、摩阻力系數(shù)、材料的晶粒大小和接觸應(yīng)力對軌道接觸疲勞的影響。目前,有限元分析已經(jīng)是研究滾動接觸疲勞的重要方法。

本文基于內(nèi)聚力模型的損傷起始準(zhǔn)則和損傷演化規(guī)律,利用VUMAT子程序結(jié)合連續(xù)損傷力學(xué)構(gòu)造了新的損傷演化方式,實(shí)現(xiàn)循環(huán)加載下的損傷累積,擬建立基于內(nèi)聚力模型的疲勞損傷累積失效模型,模擬含夾雜物軸承次表層裂紋從萌生至擴(kuò)展到接觸面的動態(tài)過程,以更好地分析、預(yù)測含夾雜物軸承鋼滾動接觸時的裂紋萌生與擴(kuò)展,并且分析載荷和摩擦因數(shù)對裂紋萌生與擴(kuò)展過程的影響。

1 基于內(nèi)聚力模型的裂紋萌生與擴(kuò)展模擬

1.1 內(nèi)聚力單元的損傷起始準(zhǔn)則

內(nèi)聚力模型遵循的斷裂準(zhǔn)則為牽引-分離定理,通過應(yīng)力-張開位移曲線定義單元失效方式。內(nèi)聚力模型的雙線性本構(gòu)關(guān)系如圖1所示,圖中Tmax為內(nèi)聚力單元進(jìn)入損傷階段前的最大牽引力,K0為單元剛度,GIC為材料斷裂能,δin為損傷起始點(diǎn),δfail為失效斷裂點(diǎn)。牽引力為內(nèi)聚力單元受到外部應(yīng)力F作用后反饋的單元內(nèi)力,用來維持單元的完整性。在起裂階段,隨著牽引位移δ的不斷增加,牽引力隨之增加,達(dá)到損傷起始點(diǎn)δin后,材料進(jìn)入損傷階段,材料剛度退化,對內(nèi)聚力單元的拉力增大,反饋的牽引力反而減小,直到牽引位移達(dá)到δfail,牽引力減小到0,單元斷裂失效。

圖1 內(nèi)聚力模型的雙線性本構(gòu)關(guān)系

1.2 基于連續(xù)損傷力學(xué)的疲勞損傷累積過程

損傷演化模型用于描述內(nèi)聚力單元進(jìn)入損傷階段后材料不斷退化,內(nèi)部牽引力不斷下降至0的過程。為模擬高周循環(huán)載荷作用下的疲勞損傷累積失效過程,在材料本構(gòu)應(yīng)力應(yīng)變關(guān)系中引入一個基于連續(xù)損傷力學(xué)的變量D,定義為材料的損傷變量,用于降低材料的剛度。初始牽引位移為0時材料沒有損傷,此時初始總損傷值D為0;當(dāng)材料完全損傷時,總損傷值D為1。

高周疲勞的內(nèi)聚力模型認(rèn)為總損傷為靜態(tài)過載損傷和動態(tài)循環(huán)疲勞載荷損傷之和。單元靜態(tài)損傷值D1可采用內(nèi)聚力單元自帶的基于位移的損傷演化來表示,定義動態(tài)循環(huán)疲勞損傷值為D2,獲得總損傷值D為

D=D1+D2。

(1)

內(nèi)聚力模型自帶的基于位移的損傷演化是靜態(tài)線性的,可表示為

(2)

當(dāng)δ=δin時材料剛進(jìn)入損傷階段,D1=0,單元開始損傷;δ=δfail時單元損傷結(jié)束,D1=1,單元失效。

對于D2,由于軸承的高周疲勞中微觀晶粒的塑性應(yīng)變很小,假設(shè)這種疲勞損傷為準(zhǔn)脆性損傷,非線性方程的廣泛通用形式為

(3)

式中:N為循環(huán)次數(shù);Δσ為相關(guān)應(yīng)力的變化范圍;τs為晶界的抗應(yīng)力;σm為法向平均應(yīng)力;m為損傷定律指數(shù)。

由于本文所研究的軸承滾子與內(nèi)圈滾道接觸模型的載荷為壓力,因此單元節(jié)點(diǎn)均處于壓縮狀態(tài),主應(yīng)力對該模型主要產(chǎn)生的Ⅱ型滑開裂紋的擴(kuò)展只能起到抑制作用。可假設(shè)沿單元切向作用的剪應(yīng)力是造成疲勞損傷、裂紋萌生與擴(kuò)展的唯一應(yīng)力,此時法向平均應(yīng)力σm對剪切裂紋的形成沒有影響,由此可得動態(tài)循環(huán)疲勞損傷值D2的損傷演化方程為與切向應(yīng)力有關(guān)的函數(shù),即

(4)

式中:Δτs為作用于單元的切向應(yīng)力變化范圍。

則內(nèi)聚力單元內(nèi)部反饋的牽引力為

T′=(1-D)Tmax,

(5)

內(nèi)聚力單元的剛度為

K′=K(1-D),

(6)

式中:T′為材料損傷后內(nèi)部牽引力;K′為材料損傷后更新的單元剛度;K為單元初始剛度。

在高周疲勞循環(huán)中,為提高計算效率采用Lemaitre提出的“jump-in cycles”算法,該方法已廣泛用于基于連續(xù)損傷力學(xué)的滾動接觸疲勞模擬。假設(shè)在損傷值增量ΔD為定值的前提下,一定的循環(huán)次數(shù)ΔN內(nèi)應(yīng)力變化范圍保持不變,如圖2所示。在每個循環(huán)段ΔN結(jié)束后對內(nèi)聚力單元剛度K和總損傷值D進(jìn)行更新,直到總損傷值D達(dá)到1,判斷該內(nèi)聚力單元失效。

圖2 分段線性損傷演化過程

選擇損傷增量ΔD為0.2,疲勞累積失效模型的損傷演化步驟如下:

1)將每個內(nèi)聚力單元的初始總損傷值和初始循環(huán)次數(shù)設(shè)置為零,即

(7)

式中:j為第j個單元,j=1,2,…,n;i為第i個循環(huán)段,i=1,2,…,n;n為內(nèi)聚力單元個數(shù)。

3)計算每個內(nèi)聚力單元的動態(tài)疲勞損傷速率,即

(8)

4)選擇具有最大動態(tài)疲勞損傷速率的內(nèi)聚力單元作為當(dāng)前循環(huán)段i的臨界單元,即

(9)

5)當(dāng)前循環(huán)段i中的循環(huán)數(shù)為

(10)

6)此時循環(huán)次數(shù)為

N=N+ΔNi。

(11)

7)在下個循環(huán)段開始時,內(nèi)聚力單元的損傷值更新為

(12)

單元剛度更新為

(13)

重復(fù)步驟2—7,直到損傷值D達(dá)到臨界值1判斷單元失效。當(dāng)內(nèi)聚力單元失效時,單元內(nèi)部分離,認(rèn)為微裂紋在該失效位置處萌生,此時的循環(huán)次數(shù)Nin即裂紋初次萌生階段。該循環(huán)過程持續(xù)進(jìn)行,將有新的單元失效,可以模擬裂紋的擴(kuò)展行為。直到當(dāng)前失效單元存在z坐標(biāo)為0的節(jié)點(diǎn),即微裂紋擴(kuò)展至軸承滾道表面時,認(rèn)為將在該位置發(fā)生剝落現(xiàn)象,此時的循環(huán)次數(shù)Nfail即軸承總壽命。

2 有限元模型

2.1 軸承滾動接觸有限元模型

滾子與內(nèi)圈的滾動接觸中,可將圓柱滾子軸承的內(nèi)圈表面擴(kuò)展成無限半空間,將滾子與滾道之間的法向接觸等效模擬為赫茲接觸,得到分布載荷為

(14)

式中:x為接觸點(diǎn)的橫坐標(biāo);a為接觸半寬;pmax為最大接觸應(yīng)力。

采用ABAQUS軟件建立NU308型內(nèi)圈無擋邊圓柱滾子軸承滾子與滾道接觸區(qū)仿真模型,在模型中距離表面0.5a深度處設(shè)定直徑為20 μm的圓形夾雜物,模型參數(shù)見表1。

表1 仿真模型參數(shù)

2.2 模型載荷、邊界及材料參數(shù)

研究滾動接觸疲勞損傷累積導(dǎo)致的裂紋萌生與擴(kuò)展過程,需要模擬滾子在內(nèi)滾道上循環(huán)滾動的運(yùn)動過程。如圖3所示,通過使表面壓力分布載荷在接觸面上x軸坐標(biāo)為-2.5a~2.5a區(qū)域內(nèi),以離散的步驟單向從右向左移動進(jìn)行模擬,施加載荷的有效區(qū)域?yàn)?2a~2a,將一次循環(huán)分解為多個離散載荷。移動載荷的法向分量為

圖3 滾動接觸模擬載荷示意圖

(15)

式中:xc為載荷中心的橫坐標(biāo)。

移動載荷的切向分量方向?yàn)閺淖笙蛴遥杀硎緸榕c法向分量相關(guān)的函數(shù),即

t(x)=fμP(x),

(16)

式中:fμ為摩擦因數(shù)。

將模型表面-2a~2a區(qū)域分割為10塊,第1個分析步為第1塊區(qū)域施加載荷,載荷分布方式為自定義場,以此進(jìn)行每個分析步的載荷設(shè)置,并設(shè)置上個分析步的載荷不傳遞到下個分析步;以同樣的方式在每個分析步中設(shè)置表面切向載荷。邊界約束方式為模型底部完全固定,模型左右邊界約束x軸方向自由度,在模型所有單元間插入厚度為0的內(nèi)聚力單元。

內(nèi)聚力單元的最大牽引力Tmax與鋼材屈服強(qiáng)度相關(guān),其彈性模量、剪切模量、密度與鋼基體材料相同。計算獲得相關(guān)參數(shù),并結(jié)合GCr15的S-N曲線獲得(4)式中參數(shù)τs和m的值,見表2。獲得滾動軸承赫茲接觸內(nèi)聚力有限元模型,如圖4所示。

表2 內(nèi)聚力單元參數(shù)

圖4 滾動接觸內(nèi)聚力模型

2.3 損傷累積有限元分析流程

ABAQUS自帶的內(nèi)聚力模型在受載損傷過程中只考慮了靜態(tài)過載損傷,需要通過子程序VUMAT將基于連續(xù)損傷力學(xué)的動態(tài)循環(huán)疲勞損傷與靜態(tài)過載損傷結(jié)合,并聯(lián)合ABAQUS自帶的python腳本端口,提取每個循環(huán)段加載后ABAQUS程序計算結(jié)果ODB文件中每個單元的總損傷值、應(yīng)力、剛度作為下個循環(huán)段的初始值,生成新的inp文件并提交ABAQUS進(jìn)行下個循環(huán)段的計算。損傷值等于1時認(rèn)為單元失效,子程序中根據(jù)損傷值更改狀態(tài)變量并傳入主程序進(jìn)行單元刪除,繼續(xù)判斷失效單元的節(jié)點(diǎn)z軸坐標(biāo)是否等于0,不是則進(jìn)入下個循環(huán),是則結(jié)束循環(huán),輸出完整的疲勞損傷過程,分析流程如圖5所示。

圖5 損傷累積有限元分析流程

3 仿真結(jié)果與討論

3.1 裂紋的萌生與擴(kuò)展

有限元模型的仿真結(jié)果如圖6所示。全過程包括裂紋的萌生、擴(kuò)展、進(jìn)一步擴(kuò)展以及擴(kuò)展至表面。

圖6 疲勞裂紋的萌生與擴(kuò)展(pmax=2.0 GPa)

隨模型受載循環(huán)次數(shù)的增加,循環(huán)次數(shù)Nin為3.23×106時萌生了第1條裂紋,所處位置為夾雜物周圍。隨著循環(huán)過程繼續(xù),多個內(nèi)聚力單元失效,初始微裂紋圍繞夾雜物逐漸擴(kuò)展且萌生了新的裂紋分叉。裂紋逐漸擴(kuò)展為更長的裂紋并向表面延伸,循環(huán)次數(shù)Nfail為17.92×106時裂紋擴(kuò)展到軸承滾道表面,預(yù)計將會在表面形成橢圓形剝落。裂紋擴(kuò)展時期的循環(huán)次數(shù)占總循環(huán)次數(shù)的81.97%。這與文獻(xiàn)[11,14]的軸承滾動接觸仿真模型研究結(jié)論一致,與文獻(xiàn)[17]的試驗(yàn)觀測結(jié)果也基本相符,認(rèn)為在接觸應(yīng)力較低的高周循環(huán)中,滾動接觸的疲勞裂紋擴(kuò)展階段平均壽命占總壽命的比例較大。

3.2 載荷對裂紋萌生與擴(kuò)展的影響

考慮載荷對裂紋萌生與擴(kuò)展的影響,將載荷增大到2.5 GPa,仿真結(jié)果如圖7所示。隨著模型受載循環(huán)次數(shù)的增加,循環(huán)次數(shù)Nin為1.06×106時萌生了第1條裂紋,所處位置為夾雜物周圍;循環(huán)次數(shù)Nfail為2.32×106時裂紋擴(kuò)展到軸承滾道表面。裂紋擴(kuò)展時期的循環(huán)次數(shù)占總循環(huán)次數(shù)的54.31%。相比載荷為2.0 GPa的情況,載荷增大為2.5 GPa后,裂紋萌生位置與擴(kuò)展方式基本不變,但裂紋擴(kuò)展壽命占總壽命的百分比明顯減小。

圖7 更高載荷下疲勞裂紋的萌生與擴(kuò)展(pmax=2.5 GPa)

3.3 摩擦因數(shù)對裂紋萌生與擴(kuò)展的影響

考慮摩擦因數(shù)對裂紋萌生與擴(kuò)展的影響,將摩擦因數(shù)增大到0.15,仿真結(jié)果如圖8所示。循環(huán)次數(shù)Nin為3.02×106時萌生了第1條裂紋,所處位置相比摩擦因數(shù)為0.05時相對表面距離明顯減小,位于夾雜物頂部;循環(huán)次數(shù)Nfail為12.74×106時裂紋擴(kuò)展到軸承滾道表面。裂紋擴(kuò)展時期的循環(huán)次數(shù)占總循環(huán)次數(shù)的76.30%。相比摩擦因數(shù)為0.05情況,摩擦因數(shù)增大到0.15后裂紋萌生位置距離表面更近,且裂紋擴(kuò)展到表面時夾雜物右側(cè)的裂紋明顯更多,這是由于接觸區(qū)摩擦力指向右。裂紋萌生壽命和擴(kuò)展壽命都相對減小且裂紋擴(kuò)展壽命占總壽命的百分比略有減小。

圖8 更大摩擦因數(shù)下疲勞裂紋的萌生與擴(kuò)展(fμ=0.15)

4 結(jié)論

建立了一種結(jié)合連續(xù)損傷力學(xué)的內(nèi)聚力模型,用于模擬滾動接觸循環(huán)加載下夾雜物周圍的裂紋萌生與擴(kuò)展,得出以下結(jié)論:

1)該模型可以模擬含夾雜物軸承次表層裂紋從萌生至擴(kuò)展到接觸面的過程,裂紋萌生位置為夾雜物周圍,初始壽命Nin為3.23×106,總壽命Nfail為17.92×106,裂紋擴(kuò)展時期的循環(huán)次數(shù)占總循環(huán)次數(shù)的81.97%。

2)在載荷較低的高周循環(huán)中,裂紋首先萌生在夾雜物周圍,滾動接觸的疲勞裂紋擴(kuò)展階段平均壽命占總壽命的比例較大。在較大的接觸載荷作用下,裂紋萌生位置與擴(kuò)展方式?jīng)]有太大變化,但裂紋擴(kuò)展壽命占總壽命的百分比明顯減小。

3)增大接觸區(qū)的摩擦因數(shù),裂紋萌生位置更接近接觸區(qū)表面及夾雜物右側(cè),即摩擦力方向的裂紋數(shù)量明顯增多。裂紋萌生壽命和擴(kuò)展壽命都相對減小且裂紋擴(kuò)展壽命占總壽命的百分比略有減小。

猜你喜歡
裂紋有限元模型
一半模型
裂紋長度對焊接接頭裂紋擴(kuò)展驅(qū)動力的影響
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
3D打印中的模型分割與打包
磨削淬硬殘余應(yīng)力的有限元分析
預(yù)裂紋混凝土拉壓疲勞荷載下裂紋擴(kuò)展速率
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 亚洲国产清纯| 中文精品久久久久国产网址| 日韩欧美国产区| 最新国产网站| 久久久久无码精品国产免费| 色哟哟色院91精品网站| 九色在线观看视频| 亚洲黄色片免费看| 欧美日韩动态图| 国产成人精品一区二区秒拍1o| 亚洲第一黄片大全| 久996视频精品免费观看| 高清无码不卡视频| 欧美国产中文| 日韩精品免费一线在线观看| 色欲国产一区二区日韩欧美| 日韩国产综合精选| 亚洲国产成人麻豆精品| 国产一区二区网站| 欧美日韩另类国产| 婷婷色在线视频| 欧美va亚洲va香蕉在线| 四虎精品黑人视频| 国产v精品成人免费视频71pao| 亚洲精品你懂的| 成人va亚洲va欧美天堂| 一级毛片基地| 四虎永久在线视频| av一区二区人妻无码| 亚洲一区二区三区中文字幕5566| 伊人久久青草青青综合| 国产一在线| 亚洲三级影院| 亚卅精品无码久久毛片乌克兰| 三上悠亚在线精品二区| 无码福利视频| 在线免费无码视频| 干中文字幕| 亚洲精品图区| 亚洲精品午夜无码电影网| 精品国产www| 亚洲精品动漫| 天天干伊人| 色偷偷综合网| 免费一级毛片在线观看| 日本91在线| 欧美区日韩区| 毛片免费在线| 久草性视频| 国产精品手机视频一区二区| 人人91人人澡人人妻人人爽| 国产三级视频网站| 无码电影在线观看| 国产成人区在线观看视频| 四虎国产永久在线观看| 亚洲美女久久| 成人小视频网| 久久久久久午夜精品| 91探花国产综合在线精品| 黄色福利在线| 天天色综网| 伊人中文网| 亚洲美女AV免费一区| 欧美综合区自拍亚洲综合天堂| 久久综合丝袜日本网| 亚洲国模精品一区| 国产成人无码久久久久毛片| 日韩经典精品无码一区二区| 成人福利在线免费观看| 国产精品成| 99一级毛片| 亚洲精品无码久久久久苍井空| 国产精品3p视频| 99在线免费播放| 99热最新网址| 97久久精品人人做人人爽| 国产视频自拍一区| 亚洲 日韩 激情 无码 中出| 欧美日韩一区二区在线播放| 色偷偷男人的天堂亚洲av| 久久99久久无码毛片一区二区| 日本欧美视频在线观看|