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

基于ANSYS的輪軌接觸疲勞分析

2015-06-07 11:22:58曹世豪文良華江曉禹
關(guān)鍵詞:裂紋有限元模型

曹世豪,李 煦,文良華,江曉禹

(1.西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室,四川 成都 610031;2.西南交通大學(xué) 力學(xué)與工程學(xué)院,四川 成都 610031)

?

基于ANSYS的輪軌接觸疲勞分析

曹世豪1,李 煦2,文良華2,江曉禹2

(1.西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室,四川 成都 610031;2.西南交通大學(xué) 力學(xué)與工程學(xué)院,四川 成都 610031)

通過分析ANSYS求解斷裂參數(shù)K因子的原理位移外推法,并比較ANSYS仿真計(jì)算與斷裂力學(xué)的解析解結(jié)果,表明ANSYS計(jì)算K因子的精確度可以達(dá)到99.9%。對鋼軌表面存在裂紋的輪軌接觸疲勞問題進(jìn)行研究,在不同裂紋角度下,獲得不同位置的裂紋尖端應(yīng)力強(qiáng)度因子。結(jié)果表明:隨著裂紋角度的增加,應(yīng)力強(qiáng)度因子KI增加而KII減小;在裂紋角度比較小時(shí),裂紋以滑開型破壞為主,隨著裂紋角度的增加,其破壞形式向張開型破壞轉(zhuǎn)變;疲勞裂紋最危險(xiǎn)的位置發(fā)生在接觸斑邊緣位置,在鋼軌養(yǎng)護(hù)時(shí),應(yīng)選用固體潤滑劑。

機(jī)車工程;位移外推法;輪軌接觸;疲勞裂紋;角度

輪軌接觸疲勞是指輪軌接觸過程中,在接觸區(qū),由于車輪對鋼軌的循環(huán)力作用,使得鋼軌表面或次表面形成微裂紋,隨后微裂紋擴(kuò)展,導(dǎo)致鋼軌表面大塊剝離,甚至發(fā)生斷裂[1]。世界各國對輪軌接觸疲勞進(jìn)行了大量的分析研究,對輪軌滾動接觸疲勞的產(chǎn)生機(jī)理說法不一。

M. C. Dubourg等[2-3]認(rèn)為,過載或者重載是造成這些部位裂紋形成的原因,而進(jìn)一步的研究表明剪切應(yīng)力在裂紋的形成上起到了很大作用。江曉禹等[4-6]認(rèn)為,對于有表面微觀粗糙度的輪軌接觸區(qū),鋼軌表面存在嚴(yán)重的壓應(yīng)力和較大的殘余拉應(yīng)力,這些殘余拉應(yīng)力可能是造成鋼軌表面微觀裂紋形成或擴(kuò)展的重要因素;并且降雨或油等液態(tài)介質(zhì)粘附在粗糙表面的凹陷部位,大大加快了鋼軌的疲勞破壞速度。J.W.Ringberg[7]認(rèn)為,對鋼軌表面短裂紋進(jìn)行分析時(shí)必須考慮材料的塑性。J.Seo等[8]認(rèn)為,對于鋼軌表面長裂紋,其破壞以滑開型為主,在全滑動運(yùn)行狀態(tài)下其擴(kuò)展速率隨著摩擦系數(shù)的增加而增大。

隨著鐵路客貨運(yùn)量的增大和列車速度的提高,輪軌滾動接觸疲勞破壞變得越來越嚴(yán)重,尤其是高速重載線路,直接危害行車安全。例如,2000年10月17日,英國的一列高速列車從倫敦的King’s Cross開往Leeds的途中,由于曲線外側(cè)鋼軌的斷裂造成4人死亡70人受傷的重大事故[9]。為減少鋼軌的疲勞破壞,最有效的辦法就是提高鋼軌的抗疲勞性能和對鋼軌斷面進(jìn)行打磨;通過鋼軌打磨可以有效控制和消除疲勞裂紋,但是過度頻繁地對鋼軌打磨會增加鋼軌的維修費(fèi)用[10]。因此,對含表面裂紋的鋼軌斷裂疲勞性能進(jìn)行研究,能夠更加清晰的認(rèn)識鋼軌的疲勞斷裂機(jī)理,進(jìn)而改進(jìn)鋼軌的抗疲勞性能,延長鋼軌的使用壽命;這樣既能減少經(jīng)營成本,又能降低運(yùn)行風(fēng)險(xiǎn)。

1 ANSYS計(jì)算K因子原理及驗(yàn)證

1.1 ANSYS計(jì)算K因子的原理

大型通用有限元軟件ANSYS基于位移外推法[11-12]進(jìn)行應(yīng)力強(qiáng)度因子的計(jì)算,其實(shí)現(xiàn)方法有兩種:一種是在裂紋尖端建立1/4節(jié)點(diǎn)的奇異性單元,由ANSYS直接求出應(yīng)力強(qiáng)度因子值;另一種是先計(jì)算出裂紋尖端附近區(qū)域節(jié)點(diǎn)的位移,再線性擬合,這種方法不一定需要考慮奇異性,只需要單元網(wǎng)格足夠精細(xì)即可。式(1)[13]詳細(xì)的描述了裂紋尖端附近節(jié)點(diǎn)的位移與應(yīng)力強(qiáng)度因子之間的關(guān)系,其參考坐標(biāo)系如圖1。

圖1 二維模型的裂紋坐標(biāo)系

(1)

式中:r,θ為裂紋尖端處的極坐標(biāo);u,v分別為x,y方向的位移分量;G為材料的剪切模量;K為裂紋的應(yīng)力強(qiáng)度因子;κ是與材料泊松比有關(guān)的系數(shù),關(guān)系如式(2):

(2)

當(dāng)θ=±180o時(shí),式(1)可簡化為:

(3)

當(dāng)r→0時(shí),裂紋尖端的KⅠ,KⅡ?yàn)椋?/p>

(4)

1.2 ANSYS計(jì)算K因子的實(shí)例驗(yàn)證

本次驗(yàn)證實(shí)驗(yàn)?zāi)P蜑楹砻媪鸭y的有限寬板,長為0.8 mm,寬為0.4 mm,裂紋長度為0.1 mm,為表面垂直裂紋,在上邊界中間位置,如圖2。其中材料的彈性模量為210 GPa,泊松比為0.3,兩邊施加均勻分布拉力,大小為10 MPa。

圖2 含表面裂紋有限寬板的有限元模型

表1為不同方法計(jì)算的應(yīng)力強(qiáng)度因子。由表1可知,通過位移外推法求K因子時(shí),越接近裂紋尖端,其值越大;當(dāng)r→0時(shí),即為最終所求的K因子,大小為0.265 5 MPa·m0.5;而當(dāng)節(jié)點(diǎn)離裂紋尖端比較遠(yuǎn)的時(shí)候就會出現(xiàn)比較明顯的誤差,其誤差最大可以達(dá)到5%;因此求K因子時(shí),離裂紋尖端越近越好。ANSYS直接計(jì)算的KI的大小為0.265 7 MPa·m0.5,而斷裂力學(xué)解析解的大小為0.265 9 MPa·m0.5,兩者的相似度達(dá)到99.9%,可以說ANSYS在線彈性范圍內(nèi)計(jì)算應(yīng)力強(qiáng)度因子的精確度非常高。

表1 不同方法計(jì)算的應(yīng)力強(qiáng)度因子

(續(xù)表1)

節(jié)點(diǎn)與裂尖的距離r/μm對應(yīng)節(jié)點(diǎn)張開位移Δv/nmKI/(MPa·m0.5)KI=ΔvGκ+12πrKI=limr→0ΔvGκ+12πrANSYS直接計(jì)算斷裂力學(xué)解析解8.7511.860.263013.0014.410.262217.2516.550.261425.7520.100.25980.26550.26570.265930.0021.640.258240.0024.820.256770.0032.430.254890.0036.850.2658

2 輪軌接觸疲勞分析

輪軌接觸疲勞分析包含龐大的單元數(shù)量,并且接觸問題存在高度的邊界非線性[14],同時(shí)還需考慮材料的塑性,因此該問題的求解需要大量的迭代過程,極不容易收斂。為了減少計(jì)算量節(jié)省計(jì)算時(shí)間,本次模擬將輪軌接觸疲勞分析分為輪軌的接觸分析和鋼軌的疲勞分析兩部分。首先通過輪軌接觸分析計(jì)算出輪軌間的接觸應(yīng)力,隨后在進(jìn)行鋼軌疲勞分析的時(shí)候?qū)⑶懊嬗?jì)算的連續(xù)分布接觸應(yīng)力以節(jié)點(diǎn)力的形式施加在模型上,以此力的作用等效車輪對鋼軌的作用;車輪在鋼軌上的滾動效果通過荷載在模型上的位置移動來實(shí)現(xiàn),如圖3。

圖3 輪軌接觸疲勞模型

2.1 有限元模型

根據(jù)我國鐵路主要干線采用鋼軌類型的實(shí)際狀況,本次分析的鋼軌類型以60 kg/m鋼軌為基準(zhǔn);模型為二維含表面裂紋的鋼軌模型,且為平面應(yīng)變問題;有限元模型高為176 mm,長為600 mm,裂紋長度為50 μm;鋼軌的材料取U71Mn-鋼,其力學(xué)性能如表2[15-16];單元類型為二維8節(jié)點(diǎn)奇異性單元PLANE183,裂紋附近和荷載作用區(qū)域的單元尺寸為10 μm,鋼軌兩邊和底邊單元尺寸為10 mm;在不考慮軌枕的影響的情況下,鋼軌下端采取全約束,整體有限元模型以及裂紋尖端局部有限元模型如圖4、圖5;整個(gè)模型(b)的單元數(shù)量為39 984個(gè),節(jié)點(diǎn)數(shù)量為123 993個(gè)。

表2 U71Mn-鋼的力學(xué)性能

圖4 含表面裂紋的鋼軌有限元模型

圖5 不同角度裂紋尖端有限元模型

2.2 結(jié)果分析

15 t軸重作用下,車輪經(jīng)過鋼軌表面裂紋的整個(gè)過程中, 30°,60°,90°這三種裂紋的應(yīng)力強(qiáng)度因子KI,KII的結(jié)果如圖7。其中,x表示接觸斑中心與裂紋之間的距離。

圖6 不同角度裂紋的KI,KII變化

由圖6(a)可知, 30°,60°,90°裂紋的應(yīng)力強(qiáng)度因子KI的最大值分別為13.3,31.2,35.8 MPa·m0.5,發(fā)生在接觸斑邊緣與裂紋之間的距離為0.4,0.5,0.5 mm;可以看出,隨著裂紋角度的增加,應(yīng)力強(qiáng)度因子KI增加,其原因是隨著裂紋角度的增加,使得裂紋尖端與鋼軌表面的距離增加,變相的增加了I型裂紋的長度;從KI最大值的位置可以看出,裂紋發(fā)生張開型破壞的最危險(xiǎn)位置基本都發(fā)生在接觸斑邊緣位置。

由圖6(b)可知, 30°,60°,90°裂紋的應(yīng)力強(qiáng)度因子KII的最大值分別為26.4,17.2,7.84 MPa·m0.5,發(fā)生在接觸斑邊緣與裂紋之間的距離為0.9,1.5,0.9;可以看出,隨著裂紋角度的增加,應(yīng)力強(qiáng)度因子KII的變化趨勢與KI剛好相反,呈減小趨勢;從KII最大值的位置可以看出,裂紋發(fā)生滑開型破壞的最危險(xiǎn)位置也都發(fā)生在接觸斑邊緣位置,但是相對于張開型破壞位置,距離接觸斑邊緣的距離要遠(yuǎn)一些。

綜合圖6可知在接觸斑壓在裂紋上時(shí),無論垂直裂紋90°還是斜裂紋30°,60°,KI幾乎為0,而KII不為0,也就是說在此階段裂紋是否擴(kuò)展主要由KII的大小來確定。斜裂紋的KII要比垂直裂紋的大得多,在接觸斑壓在裂紋上時(shí),垂直裂紋的KII值小于材料的裂紋擴(kuò)展門檻值ΔKth(2.2 MPa·m0.5),此時(shí)裂紋不擴(kuò)展,而斜裂紋的KII值都已經(jīng)大于門檻值ΔKth,也就是說此階段裂紋不停的在擴(kuò)展。

由圖8可知, 30°裂紋應(yīng)力強(qiáng)度因子峰值KII max是KI max的1.98倍,60°裂紋的KII max是KI max的0.55倍,90°裂紋的KII max是KI max的0.22倍。鋼軌表面的疲勞裂紋屬于張開型和滑開型同時(shí)存在的復(fù)合型裂紋,在裂紋角度比較小的時(shí)候,裂紋以滑開型破壞為主,隨著裂紋角度的增加,其破壞形式向張開型破壞為主轉(zhuǎn)變。

圖7 不同角度裂紋的Kmax變化

總之,車輪在趨向裂紋的時(shí)候,裂紋開始緩慢張開,此時(shí)若鋼軌表面存在水或油等液態(tài)介質(zhì)時(shí),它們會隨著裂紋的張開而進(jìn)入裂紋;由于裂紋最大張開位置發(fā)生在接觸斑壓在裂紋一小段的地方,隨后裂紋面在閉合擠壓的過程中,裂紋里的液體無法被排出;由于液體的不可壓縮性,會使裂紋表面受到垂直壓力阻止裂紋閉合,加速裂紋擴(kuò)展[17]。另一方面,鋼軌養(yǎng)護(hù)中,為了減少鋼軌磨損,通常會采取涂抹潤滑劑的辦法,降低輪軌間的摩擦系數(shù),以達(dá)到減小磨損的目的;此時(shí)應(yīng)避免使用液體潤滑劑,優(yōu)先采用固體潤滑劑。

3 結(jié) 論

1)計(jì)算裂紋尖端的K因子時(shí),所選路徑距離裂紋尖端越近,結(jié)果精確度越高;且ANSYS計(jì)算K因子的精確度可以達(dá)到99%以上。

2)裂紋發(fā)生疲勞破壞的最危險(xiǎn)位置基本都發(fā)生在接觸斑邊緣位置,在鋼軌養(yǎng)護(hù)時(shí),優(yōu)先使用固體潤滑劑。

3)鋼軌表面疲勞裂紋屬于張開型和滑開型同時(shí)存在的復(fù)合型裂紋;隨著裂紋角度的增加,應(yīng)力強(qiáng)度因子KI增加,而KII卻減小;在裂紋角度比較小的時(shí)候,裂紋以滑開型破壞為主,隨著裂紋角度的增加,其破壞形式向張開型破壞轉(zhuǎn)變。

[1] 陳顏堂,劉冬雨,方鴻生.鋼軌鋼的滾動接觸疲勞[J].鋼鐵研究學(xué)報(bào),2000,12(5):50-54. Chen Yantang,Liu Dongyu,Fang Hongsheng.Rolling contact fatigue of rail steel [J].Journal of Iron and Steel Research,2000,12(5):50-54.

[2] Dubourg M C,Villechaise B.Analysis of multiple fatigues cracks-part I theory [J].ASME,Journal of Tribology,1992,144:455-461.

[3] Dubourg M C,Villechaise B,Godet M.Analysis of multiple fatigue cracks-part II results [J].ASME,Journal of Tribology,1992,144:462-468.

[4] 江曉禹,金學(xué)松.輪軌接觸表面有液態(tài)介質(zhì)時(shí)的接觸問題研究[J].工程力學(xué),2005,22(2):27-32. Jiang Xiaoyu,Jin Xuesong.The analysis of wheel/rail contact surface with liquid [J].Engineering Mechanics,2005,22(2):27-32.

[5] 江曉禹,金學(xué)松.考慮表面微觀粗糙度的輪軌接觸彈塑性分析[J].西南交通大學(xué)學(xué)報(bào),2001,36(6):588-590. Jiang Xiaoyu,Jin Xuesong.Elastic-Plastic analysis of contact problems of wheel/rail with surface micro-roughness [J].Journal of Southwest Jiaotong University,2001,36(6):588-590.

[6] 江曉禹,金學(xué)松.輪軌間的液態(tài)介質(zhì)和表面微觀粗糙度對接觸表面疲勞損傷的影響[J].機(jī)械工程學(xué)報(bào),2004,40(8):18-23. Jiang Xiaoyu,Jin Xuesong.Influence of liquid and micro-roughness on the fatigue of wheel-rail contact [J].Journal of Mechanical Strength,2004,40(8):18-23.

[7] Ringberg J W.Shear made growth of short surface-breaking RCF cracks [J].Wear,2005,258:955-963.

[8] Seo J, Kown S,Jun H,et al.Fatigue crack growth behavior of surface crack in rails [J].Procedia Engineering,2010,15(2):865-872.

[9] 史密斯.鋼軌滾動接觸疲勞的進(jìn)一步研究[J].中國鐵道科學(xué),2002,23(3):43-46. Smith R A.Rolling contact fatigue of rails:what remains to be done [J].China Railway Science,2002,23(3):43-46.

[10] 金學(xué)松,杜星,郭俊,等.鋼軌打磨技術(shù)研究進(jìn)展[J].西南交通大學(xué)學(xué)報(bào),2010,45(1):1-11. Jin Xuesong,Du Xing,Guo Jun,et al.State of arts of research on rail grinding [J].Journal of Southwest Jiaotong University,2010,45(1):1-11.

[11] 趙海濤,戰(zhàn)玉寶,楊永騰.基于ANSYS的應(yīng)力強(qiáng)度因子計(jì)算[J].煤礦機(jī)械,2007,28(2):22-23. Zhao Haitao,Zhan Yubao,Yang Yongteng.Calculation of stress intensity factors using ANSYS [J].Coal Mine Machinery,2007,28(2):22-23.

[12] 趙偉,向陽開.三點(diǎn)彎曲梁裂縫應(yīng)力強(qiáng)度因子有限元分析[J].重慶交通大學(xué)學(xué)報(bào):自然科學(xué)版,2007,26(增刊1):1-3. Zhao Wei,Xiang Yangkai.Research on stress intensity factor of the concrete based on FEM [J].Journal of Chongqing Jiaotong University:Natural Science,2007,27(Sup1):1-3.

[13] Alegre J M,Cuesta I I.Some aspects about the crack growth FEM simuiations under mixed-mode loading[J].International Journal of Fatigue,2010,32(4):1090-1095

[14] 魏龍海,陳春光,王明年,等.三維離散元模型及計(jì)算參數(shù)選取研究[J].重慶交通大學(xué)學(xué)報(bào):自然科學(xué)版,2008,27(4):618-623. Wei Longhai,Chen Chunguang,Wang Mingnian,et al.Study on three-dimensional discrete element method and parameter adoption [J].Journal of Chongqing Jiaotong University:Natural Science,2008,27(4):618-623.

[15] 王文健,劉啟躍.PD3和U71Mn鋼軌疲勞裂紋擴(kuò)展速率研究[J].機(jī)械強(qiáng)度,2007,29(6):1026-1029. Wang Wenjian,Liu Qiyue.Study on fatigue crack growth rate of PD3 and U71Mn rail steel [J].Journal of Mechanical Strength,2007,29(6):1026-1029.

[16] 周小林,向延念,陳秀芳.U71Mn50 kg/m普通碳素鋼軌鋼疲勞裂紋擴(kuò)展速率實(shí)驗(yàn)研究[J].中國鐵道科學(xué),2004,25(3):86-90. Zhou Xiaolin,Xiang Yannian,Chen Xiufang.Experimental study on fatigue crack growth rate of U71Mn50 kg/m ordinary carbon rail [J].China Railway Science,2004,25(3):86-90.

[17] 金學(xué)松,杜星,郭俊,等.鋼軌打磨技術(shù)研究進(jìn)展[J].西南交通大學(xué)學(xué)報(bào),2010,45(1):1-11. Jin Xuesong,Du Xing,Guo Jun,et al.State of arts of research on rail grinding[J].Journal of Southwest Jiaotong University,2010,45(1):1-11.

Fatigue Analysis of Wheel/Rail Contact Based on ANSYS

Cao Shihao1, Li Xu2, Wen Lianghua2, Jiang Xiaoyu2

(1.MOE Key Laboratory of High-Speed Railway Engineering,Southwest Jiaotong University, Chengdu 610031, Sichuan,China;2.School of Mechanics & Engineering, Southwest Jiaotong University, Chengdu 610031, Sichuan, China)

Through analyzing the principle of displacement extrapolation method to solve the fracture parameterKand comparing the results of ANSYS calculation with the results of the fracture mechanics, it was showed that the accuracy ofKcalculated by ANSYS could reach 99.9%. Then the contact fatigue problem of wheel/rail with surface crack was analyzed. From different angles of crack, the stress intensity factors of crack tips at different locations were obtained. The results show that, the stress intensity factorKIincreases with the increase of crack angle, meanwhile, theKIIdecreases. When the crack angle is small, the major damage of crack is sliding mode; with the increase of crack angle, the major damage of crack transforms sliding mode into opening mode. The most dangerous location of crack is at the edge of the contact area, and the solid lubricant should be used firstly in the curing process of the rail.

locomotive engineering; displacement extrapolation method; wheel/rail contact; fatigue crack; angle

10.3969/j.issn.1674-0696.2015.04.34

2013-11-10;

2014-05-12

國家自然科學(xué)基金重點(diǎn)項(xiàng)目(U1134202/E050303);國家自然科學(xué)基金項(xiàng)目(11472230)

曹世豪(1988—),男,河南焦作人,博士研究生,主要從事輪軌接觸疲勞方面的研究。E-mail: 531148108@qq.com。

江曉禹(1965—),男,貴州遵義人,教授,博士生導(dǎo)師,主要從事輪軌接觸疲勞,復(fù)合材料力學(xué)方面的研究。E-mail:xiaoyujiang8@sohu.com。

U260.33

A

1674-0696(2015)04-171-05

猜你喜歡
裂紋有限元模型
一半模型
裂紋長度對焊接接頭裂紋擴(kuò)展驅(qū)動力的影響
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
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
主站蜘蛛池模板: 98超碰在线观看| 九色视频一区| 国产成人亚洲精品蜜芽影院| 中文一区二区视频| 国产三级毛片| 在线观看网站国产| 国产噜噜在线视频观看| 99久久精品国产综合婷婷| 亚亚洲乱码一二三四区| 国产精品无码AV片在线观看播放| 日韩欧美一区在线观看| av在线手机播放| 国产白浆视频| 亚洲91在线精品| 国产91蝌蚪窝| 草草影院国产第一页| 亚洲二三区| 亚洲三级色| 成年人国产视频| 高清欧美性猛交XXXX黑人猛交| 亚洲中文无码av永久伊人| 天堂在线视频精品| 午夜丁香婷婷| 青青草一区| 亚洲日韩图片专区第1页| 婷婷午夜影院| 国产精品女人呻吟在线观看| 五月天久久综合国产一区二区| 国内精品伊人久久久久7777人| 国产一区二区三区精品欧美日韩| a级毛片免费在线观看| 国产XXXX做受性欧美88| 欧美综合中文字幕久久| 日本尹人综合香蕉在线观看| 日韩一二三区视频精品| 亚洲人成网站在线播放2019| 国产主播在线一区| 亚洲国产精品无码久久一线| 亚洲区欧美区| 国产肉感大码AV无码| 成人福利在线视频| 亚洲第一页在线观看| 99久久精品国产自免费| 99久久精品久久久久久婷婷| 国产高颜值露脸在线观看| 国产欧美精品专区一区二区| 欧美日韩亚洲国产主播第一区| 欧美日韩免费观看| 日本不卡在线播放| 国产凹凸视频在线观看| 蜜芽国产尤物av尤物在线看| 色婷婷丁香| 国产美女在线免费观看| 波多野结衣视频一区二区| 亚洲天堂在线免费| 亚洲一区二区三区中文字幕5566| 亚洲欧美日韩中文字幕在线一区| 精品一区二区三区波多野结衣| 欧美日韩精品一区二区在线线 | 国产亚洲欧美在线中文bt天堂| 亚洲天堂.com| 国产男女免费完整版视频| 亚洲天堂色色人体| 国产精品视频猛进猛出| 2021国产精品自产拍在线| 国产色网站| 亚洲AV电影不卡在线观看| 国产人在线成免费视频| 2021天堂在线亚洲精品专区| 日韩国产一区二区三区无码| 欧美在线伊人| 波多野结衣久久精品| 亚洲天天更新| 亚洲一区二区三区香蕉| 欧美视频免费一区二区三区| 91无码视频在线观看| 国产视频入口| 欧美日本在线播放| 日本成人精品视频| 久久精品国产亚洲麻豆| 日本人妻丰满熟妇区| 欧美精品导航|