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

非均勻載荷下厚壁圓筒穩(wěn)態(tài)蠕變應(yīng)力的計(jì)算

2015-03-07 00:34:01王建梅苗克軍徐俊良李璞
關(guān)鍵詞:有限元方法

王建梅,苗克軍,徐俊良,李璞

?

非均勻載荷下厚壁圓筒穩(wěn)態(tài)蠕變應(yīng)力的計(jì)算

王建梅,苗克軍,徐俊良,李璞

厚壁圓筒是工程中的重要結(jié)構(gòu),廣泛應(yīng)用于化工、核電、軍工、冶金等領(lǐng)域[1],長期在高溫(0.3Tm以上,Tm為熔點(diǎn))和復(fù)雜載荷下工作時(shí)需要考慮蠕變過程[2]。關(guān)于在均勻載荷下穩(wěn)態(tài)蠕變應(yīng)力的求解問題,Naumenko等分析了蠕變對厚壁圓筒平面應(yīng)力的影響[3],Boyle基于von Mises屈服準(zhǔn)則,得到了厚壁圓筒的等效應(yīng)力解析解[4]。然而,在工程中經(jīng)常遇到的是非均勻載荷的情況[5]。關(guān)于厚壁圓筒在非均勻載荷下的彈性問題的應(yīng)力求解,吳慶良等采用復(fù)變函數(shù)法計(jì)算了外表面承受非均勻載荷作用時(shí)厚壁圓筒的平面應(yīng)力[6],鄭俊德等利用應(yīng)力函數(shù)法獲得了厚壁圓筒的應(yīng)力解析解[7]。上述研究雖然給出了一些厚壁圓筒平面應(yīng)力的解析求解方法,但均未考慮在蠕變條件下受非均勻載荷厚壁圓筒的應(yīng)力求解問題,且國內(nèi)尚無相關(guān)文獻(xiàn)介紹這方面的研究。

本文以厚壁圓筒為研究對象,利用多維應(yīng)力下穩(wěn)態(tài)蠕變的基礎(chǔ)理論與厚壁圓筒平面應(yīng)力的計(jì)算方法,推導(dǎo)出一類在外邊界固定、內(nèi)邊界受非均勻載荷情況下的穩(wěn)態(tài)蠕變應(yīng)力計(jì)算公式,提出了一種多維應(yīng)力蠕變的計(jì)算方法,應(yīng)用所得公式分析了工程中油膜軸承襯套的蠕變應(yīng)力,并以蠕變試驗(yàn)為基礎(chǔ),模擬了襯套的蠕變過程。

1 理論計(jì)算

1.1 算法推導(dǎo)

為計(jì)算承受非均勻內(nèi)壓的厚壁圓筒的穩(wěn)態(tài)蠕變應(yīng)力,本文主要考慮第二階段的蠕變,即穩(wěn)態(tài)蠕變。在穩(wěn)態(tài)蠕變階段,不考慮溫度對蠕變應(yīng)力的影響時(shí),彈性力學(xué)中的平衡方程、協(xié)調(diào)方程、幾何方程、邊界方程等仍然適用[8]。同時(shí),假設(shè)圓筒處于平面應(yīng)力狀態(tài),即σz=0,對厚壁圓筒的中間截面進(jìn)行應(yīng)力分析。

設(shè)定一個(gè)內(nèi)半徑為a、外半徑為b的厚壁圓筒,其外邊界(r=b)固定,內(nèi)邊界(r=a)受到非均勻徑向壓力p(θ)的作用,計(jì)算模型如圖1所示。

圖1 厚壁圓筒計(jì)算模型示意圖

(1)

式中:am為與壓力相關(guān)的系數(shù);m和n為整數(shù),且m≤n。

邊界條件為

(2)

式中:ur和vθ為位移分量。

根據(jù)幾何形狀和邊界載荷選取應(yīng)力函數(shù)[9]

(3)

式中:Am、Bm、Cm和Dm為未知參數(shù)。

σr、σθ和τrθ為應(yīng)力分量,表達(dá)式為

(4)

將式(3)代入式(4)得

(5)

對于多維穩(wěn)態(tài)蠕變應(yīng)力分析,將等效應(yīng)力和等效蠕變分別定義為多維應(yīng)力下的應(yīng)力分量和蠕變應(yīng)變分量,以諾頓穩(wěn)態(tài)蠕變公式為基礎(chǔ),得到如下穩(wěn)態(tài)蠕變的平面應(yīng)力應(yīng)變關(guān)系[10]

(6)

將式(5)代入式(6)得

(7)

幾何方程[9]為

(8)

將式(7)代入式(8),并對前兩項(xiàng)積分,得

(9)

式中:f1(θ)和f2(r)-∫f1(θ)dθ是積分常量。

將式(9)代入式(8)中的第3個(gè)式子,根據(jù)位移單值條件得

(10)

將式(10)代入式(9),并將邊界條件代入式(5)和式(9),得

(11)

式中

Am={3a2+m[-5a2mb2+3a2b2mm-3b2m+2(1+

m)]am}/{2(1+m)[15a2+4mb2+15a2b2+4m+

2a2mb2m+2(9a2+8b2)+9a2mb2m(a2-b2)2m2]}

(12)

Bm={a2+mb2[15a2+2mm-9a2b2m(-m+m2)+

b2m+2(16+9m2)]am}/{2(-m+m2)[15a2+4mb2+

15a2b2+4m+2a2mb2m+2(9a2+8b2)+

9a2mb2m(a2-b2)2m2]}

(13)

Cm={3a2+mb2m[5b2+2m+3a2mb2(1-m)+

3a2m+2m]am}/{2(-1+m)[15a2+4mb2+15a2b2+4m+

2a2mb2m+2(9a2+8b2)+9a2mb2m(a2-b2)2m2]}

(14)

Dm={a2+mb2+2m[-15a2b2mm-9a2m+2(m+

m2)+a2mb2(16+9m2)]am}/{2(m+m2)·

[15a2+4mb2+15a2b2+4m+2a2mb2m+2(9a2+

8b2)+9a2mb2m(a2-b2)2m2]}

(15)

將式(12)~式(15)代入式(5),可得非均勻內(nèi)壓下厚壁圓筒的穩(wěn)態(tài)蠕變應(yīng)力分量

(16)

1.2 非均勻載荷下厚壁圓筒的蠕變應(yīng)力分析

當(dāng)機(jī)械零部件的溫度超過材料的蠕變溫度時(shí),均應(yīng)按照蠕變設(shè)計(jì)準(zhǔn)則來確定設(shè)計(jì)用的許用應(yīng)力,本文采用第四強(qiáng)度理論[10-11]進(jìn)行判斷。設(shè)計(jì)中應(yīng)使零部件的最大工作應(yīng)力滿足如下關(guān)系式

(17)

式中:[σ]c為蠕變許用應(yīng)力,由不同零部件的工作要求決定。厚壁圓筒的最大蠕變應(yīng)力應(yīng)該滿足式(17)。

由式(16)可以看出,各應(yīng)力分量由余弦函數(shù)疊加得到,考慮到每一疊加項(xiàng)的蠕變分布規(guī)律類似,且疊加不改變?nèi)渥儜?yīng)力的分布規(guī)律,故只取其疊加項(xiàng)中的第一項(xiàng)進(jìn)行分析。

考慮到上述計(jì)算模型中影響穩(wěn)態(tài)蠕變應(yīng)力的主要幾何參數(shù)有厚壁圓筒的外半徑b和內(nèi)半徑a以及任意一點(diǎn)處的半徑r,因此,設(shè)定K為厚壁圓筒外半徑與內(nèi)半徑之比:K=b/a。均勻選擇厚壁圓筒的10個(gè)位置(a≤r≤b),得到其在不同K值下的蠕變應(yīng)力,如圖2所示。

(a)a≤r≤(5a+3b)/8

(b)(a+b)/2≤r≤b圖2 不同K值下厚壁圓筒的蠕變應(yīng)力分布

由圖2a可知,當(dāng)r≤(5a+3b)/8時(shí),半徑越小,則蠕變應(yīng)力隨著K值增大呈現(xiàn)出的先增大后減小的趨勢越明顯。由圖2b可知:隨著K值增大,蠕變應(yīng)力在半徑r=(a+b)/2處近似線性遞減;當(dāng)r>(a+b)/2時(shí),半徑越大,則蠕變應(yīng)力隨著K值的增大而減小的速率變得越緩慢。當(dāng)K值一定時(shí),隨著半徑的增大,厚壁圓筒的蠕變應(yīng)力逐漸降低,說明最大蠕變應(yīng)力出現(xiàn)在內(nèi)邊界上;厚壁圓筒的最大蠕變應(yīng)力隨著K值的增加呈現(xiàn)出先增大后減小的規(guī)律,在K=1.58時(shí)蠕變應(yīng)力最大。

2 油膜軸承襯套的蠕變應(yīng)力計(jì)算

油膜軸承最薄弱的零件就是襯套中的巴氏合金,蠕變是其失效的一種形式[12]。油膜軸承襯套在工作過程中長期承受油膜壓力與高溫的作用,會(huì)產(chǎn)生微觀蠕變,導(dǎo)致對油膜軸承的潤滑性能和壽命產(chǎn)生不可忽略的負(fù)面影響,因此,軸承中巴氏合金的蠕變不可忽略[13]。工程中油膜軸承襯套的外邊界是固定的,內(nèi)邊界承載區(qū)120°范圍內(nèi)承受徑向非均勻油膜壓力的作用。選擇某工況下襯套內(nèi)邊界軸向承受油膜壓力最大的中間截面,截面內(nèi)周向油膜壓力可用余弦函數(shù)表示

(18)

式中:poil(θ)為油膜軸承襯套內(nèi)壁所受的徑向油膜壓力,MPa;θ為與水平方向的夾角;bm為油膜壓力相關(guān)系數(shù),其值見表1。襯套內(nèi)半徑為110 mm,外半徑為128 mm,K=1.16。

表1 油膜壓力相關(guān)系數(shù)

將油膜軸承襯套的相關(guān)參數(shù)代入式(16)與式(17)中,得到其承載區(qū)的蠕變應(yīng)力分布,如圖3所示。

圖3 油膜軸承襯套的蠕變應(yīng)力分布

由圖3可知:在周向,襯套的蠕變應(yīng)力隨著周向節(jié)點(diǎn)坐標(biāo)的增加呈現(xiàn)先增大后減小的形態(tài)分布;在徑向,襯套的蠕變應(yīng)力隨著半徑的增大而減小。這表明,襯套所受載荷越大的區(qū)域,蠕變應(yīng)力也越大。

3 蠕變試驗(yàn)

3.1 試驗(yàn)過程

采用有限元分析軟件ANSYS求解蠕變時(shí),所需的蠕變參數(shù)由單軸蠕變拉伸試驗(yàn)確定[14]。為求得襯套內(nèi)壁材料的蠕變參數(shù),針對襯套內(nèi)壁材料的性能要求,在原有WDW-E100D型電子式萬能試驗(yàn)機(jī)的基礎(chǔ)上,設(shè)計(jì)了一套滿足國家標(biāo)準(zhǔn)[15]要求的蠕變試驗(yàn)設(shè)備,增加了應(yīng)變測量、加熱、保溫和控溫的功能。試驗(yàn)所用試件如圖4所示,試驗(yàn)設(shè)備如圖5所示。選擇13 MPa-60 ℃、18 MPa-60 ℃和20 MPa-60 ℃ 3種工況進(jìn)行蠕變試驗(yàn),根據(jù)蠕變試驗(yàn)數(shù)據(jù)繪制3種工況下的蠕變速率曲線,如圖6所示,根據(jù)蠕變試驗(yàn)結(jié)果計(jì)算的穩(wěn)態(tài)蠕變速率值見表2。

圖4 蠕變試驗(yàn)試件

圖5 蠕變試驗(yàn)設(shè)備

載荷/MPa131820穩(wěn)態(tài)蠕變速率/10-10s-10 9722 8134 444

圖6 蠕變速率曲線

3.2 諾頓穩(wěn)態(tài)蠕變本構(gòu)方程

利用ANSYS計(jì)算穩(wěn)態(tài)蠕變的諾頓公式為

(19)

不考慮溫度變化對蠕變應(yīng)力的影響,即c3=0,則式(19)變?yōu)?/p>

(20)

對式(20)兩邊取對數(shù),得

(21)

圖7 等效蠕變速率與等效蠕變應(yīng)力的關(guān)系

4 有限元模擬

根據(jù)油膜軸承襯套的實(shí)際尺寸,采用ANSYS建立二維模型。為提高求解精度,應(yīng)用隱式蠕變算法,選擇單元求解能力強(qiáng)且支持隱式蠕變算法的PLANE183單元,將諾頓蠕變系數(shù)輸入ANSYS中。考慮襯套的實(shí)際受載情況,外邊界施加周向與徑向位移約束,內(nèi)邊界承載區(qū)120°范圍內(nèi)施加徑向載荷。ANSYS求解結(jié)果如圖8所示。

圖8 承載區(qū)蠕變應(yīng)力云圖

由圖8可知:襯套周向的應(yīng)力隨著載荷的增大而增大;在襯套承受的最小油膜壓力處,蠕變應(yīng)力接近于0,表明在油膜壓力為0的地方,襯套基本不存在蠕變;在承受最大油膜壓力處,蠕變應(yīng)力最大,襯套蠕變最明顯;蠕變應(yīng)力隨著半徑的增大而減小,表明襯套在非承載區(qū)的蠕變不明顯。模擬結(jié)果與第2節(jié)中本文方法所得結(jié)果一致。

選定r為110、120 mm,θ為60°、90°分別對應(yīng)的4個(gè)位置,采用本文方法與有限元方法對蠕變應(yīng)力進(jìn)行計(jì)算,結(jié)果對比如圖9所示。

(a)r=110 mm

(b)r=120 mm

(c)θ=60°

(d)θ=90°圖9 本文方法與有限元方法計(jì)算的蠕變應(yīng)力對比

由圖9a和9b可知,本文方法與有限元方法計(jì)算的周向蠕變應(yīng)力基本相同;由圖9c和9d可知,本文方法與有限元方法計(jì)算的徑向蠕變應(yīng)力吻合較好。有限元方法較本文方法計(jì)算的結(jié)果略大,隨著半徑的增加,本文方法與有限元方法計(jì)算結(jié)果的誤差逐漸減小。本文方法與有限元方法計(jì)算結(jié)果的最大誤差在r=110 mm、θ=90°的位置,在此位置本文方法求得的蠕變應(yīng)力為0.090 88 MPa,有限元方法求得的蠕變應(yīng)力為0.094 43 MPa,二者的相對誤差僅為3.9%,表明本文方法能夠滿足工程的實(shí)際需求。

5 結(jié) 論

(1)本文推導(dǎo)出了非均勻載荷下厚壁圓筒的穩(wěn)態(tài)蠕變應(yīng)力計(jì)算公式,計(jì)算了油膜軸承襯套在承受油膜壓力時(shí)的蠕變應(yīng)力;提出了一種新的多維蠕變應(yīng)力求解方法,可以為多維蠕變應(yīng)力解析求解提供理論依據(jù)。

(2)根據(jù)蠕變試驗(yàn)數(shù)據(jù)得到了諾頓穩(wěn)態(tài)蠕變本構(gòu)方程,并利用所得參數(shù)模擬了襯套的蠕變應(yīng)力。本文方法的計(jì)算結(jié)果與模擬結(jié)果吻合,證實(shí)了本文方法的精確性,精度符合工程實(shí)際需求。

(3)得到了不同外內(nèi)徑比K下厚壁圓筒的蠕變應(yīng)力及其分布,證實(shí)最大蠕變應(yīng)力發(fā)生在厚壁圓筒內(nèi)邊界最大受力處,厚壁圓筒的最大蠕變應(yīng)力隨著壁厚的增加呈現(xiàn)出先增大后減小的規(guī)律,且當(dāng)K=1.58時(shí)蠕變應(yīng)力最大。這一結(jié)論可望為工程中厚壁圓筒的結(jié)構(gòu)設(shè)計(jì)提供理論參考。

[1] LIANG Yaping, WANG Huizhen, REN Xingmin. Analytical solution for spatially axisymmetric problem of thick-walled cylinder subjected to different linearly varying pressures along the axis and uniform pressures at two ends [J]. Science in China: Series G Physics, Mechanics and Astronomy, 2008, 51(1): 98-104.

[2] ALTENBACH H, GORASH Y, NAUMENKO K. Steady-state creep of a pressurized thick cylinder in both the linear and the power law ranges [J]. Acta Mechanica, 2008, 195(1/2/3/4): 263-274.

[3] NAUMENKO K, ALTENBACH H, GORASH Y. Creep analysis with a stress range dependent constitutive model [J]. Archive of Applied Mechanics, 2009, 79(6/7): 619-630.

[4] BOYLE J T. The creep behavior of simple structures with a stress range-dependent constitutive model [J]. Archive of Applied Mechanics, 2012, 82(4): 495-514.

[5] WU Qingliang, Lü Aizhong, GAO Yongtao, et al. Stress analytical solution for plane problem of a double-layered thick-walled cylinder subjected to a type of non-uniform distributed pressure [J]. Journal of Central South University, 2014, 21(5): 2074-2082.

[6] 吳慶良, 呂愛鐘. 一類非均布荷載作用下厚壁圓筒平面問題的應(yīng)力解析解 [J]. 工程力學(xué), 2011, 28(6): 6-10, 18. WU Qingliang, Lü Aizhong. Stress analytical solution for plane problem of a thick-walled cylinder subjected to a type of non-uniform distributed pressures [J]. Engineering Mechanics, 2011, 28(6): 6-10, 18.

[7] 鄭俊德, 張艷秋, 王文軍, 等. 非均勻載荷下套管強(qiáng)度的計(jì)算 [J]. 石油學(xué)報(bào), 1998, 19(1): 119-123. ZHENG Junde, ZHANG Yanqiu, WANG Wenjun, et al. Calculation of casing strength under nonuniform loading [J]. Acta Petrolei Sinica, 1998, 19(1): 119-123.

[8] 穆霞英. 蠕變力學(xué) [M]. 西安: 西安交通大學(xué)出版社, 1989: 60-83.

[9] 徐芝綸. 彈性力學(xué) [M]. 北京: 高等教育出版社, 2006: 55-63.

[10]平修二. 金屬材料的高溫強(qiáng)度理論設(shè)計(jì) [M]. 郭廷瑋,李安定, 徐介平, 譯. 北京: 科學(xué)出版社, 1983: 259-265.

[11]機(jī)械設(shè)計(jì)手冊編委會(huì). 機(jī)械設(shè)計(jì)手冊: 第5卷 [M]. 3版. 北京: 機(jī)械工業(yè)出版社, 2004: 31-33.

[12]SIMMONS T E, CONTARINI A, GIANNI N. 油膜軸承變形和壓力分析 [J]. 鋼鐵, 2009, 44(3): 93-96. SIMMONS T E, CONTARINI A, GIANNI N, et al. Deflection and pressure analysis of oil film bearings [J]. Iron and Steel, 2009, 44(3): 93-96.

[13]王建梅, 薛亞文, 馬立新, 等. 蠕變對巴氏合金ZChSnSb11-6力學(xué)性能和顯微組織的影響 [J]. 中國有色金屬學(xué)報(bào), 2014, 24(10): 2513-2518. WANG Jianmei, XUE Yawen, MA Lixin, et al. Influence of creep on mechanical properties and microstructures of Babbitt alloy ZChSnSb11-6 [J]. The Chinese Journal of Nonferrous Metals, 2014, 24(10): 2513-2518.

[14]GLADIUS L, SHAW K M. Creep constitutive model and component lifetime estimation: the case of niobium-modified 9Cr-1Mo steel weldments [J]. Journal of Materials Engineering and Performance, 2011, 20(7): 1310-1314.

[15]中國鋼鐵工業(yè)協(xié)會(huì). GB/T 2039—2012 金屬材料單軸拉伸蠕變試驗(yàn)方法 [S]. 北京: 中國標(biāo)準(zhǔn)出版社, 2012.

(編輯 葛趙青)

(太原科技大學(xué)教育部重型機(jī)械工程技術(shù)研究中心,030024,太原)

根據(jù)多維應(yīng)力下穩(wěn)態(tài)蠕變的基礎(chǔ)理論,采用應(yīng)力函數(shù)法推導(dǎo)出一類外邊界固定、內(nèi)邊界受非均勻載荷的厚壁圓筒穩(wěn)態(tài)蠕變應(yīng)力計(jì)算公式;分析了不同外內(nèi)徑比K下厚壁圓筒的蠕變應(yīng)力和應(yīng)力分布情況,證實(shí)最大蠕變應(yīng)力發(fā)生在厚壁圓筒內(nèi)邊界最大受力處,厚壁圓筒的最大蠕變應(yīng)力隨著壁厚的增加呈現(xiàn)出先增加后減小的規(guī)律,且當(dāng)K=1.58時(shí)蠕變應(yīng)力最大。應(yīng)用所得計(jì)算公式,針對工程中油膜軸承襯套因承受油膜壓力而發(fā)生的蠕變問題進(jìn)行了蠕變應(yīng)力解析求解。通過蠕變拉伸試驗(yàn),得到了襯套內(nèi)壁材料的諾頓穩(wěn)態(tài)蠕變本構(gòu)方程,并根據(jù)蠕變系數(shù),利用有限元分析軟件ANSYS計(jì)算了襯套的蠕變應(yīng)力。最后,比較了上述2種方法得到的蠕變應(yīng)力,表明文中給出的計(jì)算方法與有限元方法所得結(jié)果吻合,驗(yàn)證了該解析算法的正確性,從而可以為厚壁圓筒的結(jié)構(gòu)設(shè)計(jì)和多維蠕變應(yīng)力解析求解提供理論依據(jù)。

厚壁圓筒;穩(wěn)態(tài)蠕變應(yīng)力;非均勻載荷;蠕變試驗(yàn);諾頓本構(gòu)方程;有限元;油膜軸承襯套

Steady-State Creep Stress Calculation of Thick-Walled Cylinder under Non-Uniform Load

WANG Jianmei,MIAO Kejun,XU Junliang,LI Pu

(Heavy Machinery Engineering Research Center of Ministry of Education, Taiyuan University of Science and Technology, Taiyuan 030024, China)

Based on the basic theory of steady-state creep under multi-dimensional stress state, a steady-state creep stress calculation formula of thick-walled cylinder with fixed outer boundary was derived by the stress-function method under non-uniform internal pressure. The creep stress and stress distribution of several kinds of thick-walled cylinder were analyzed, which showed that the maximal creep stress occurred in the area with the largest stress on the inner boundary of the thick-walled cylinders, increased first and then decreased with the wall thickness, and the largest creep stress was produced when the ratio of outer diameter to internal diameter is 1.58. According to the formula, the steady-state creep stress of oil-film bearing bushing subject to a non-uniform oil-film pressure was analyzed, and the creep process was simulated by using the finite element software in ANSYS, where the creep parameters were obtained based on creep tensile test. Calculation results of the above two methods were analyzed, and the results obtained by the proposed method were coincident well with those by finite element method. This research may provide a theoretical guidance for the structural design of thick-walled cylinders and the analytical calculation of multi-dimensional creep stress.

thick-walled cylinder; steady-state creep stress; non-uniform load; creep test; Norton constitutive equation; finite element method; oil-film bearing bushing

2015-01-14。 作者簡介:王建梅(1972—),女,教授。 基金項(xiàng)目:國家自然科學(xué)基金資助項(xiàng)目(51205269);山西省自然科學(xué)基金資助項(xiàng)目(2012011018-2);山西省回國留學(xué)人員科研資助項(xiàng)目(2013-093)。

時(shí)間:2015-06-17

http:∥www.cnki.net/kcms/detail/61.1069.T.20150617.0902.001.html

10.7652/xjtuxb201509002

TH123

A

0253-987X(2015)09-0008-06

猜你喜歡
有限元方法
新型有機(jī)玻璃在站臺(tái)門的應(yīng)用及有限元分析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機(jī)制的探討
學(xué)習(xí)方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 无码人中文字幕| 乱人伦99久久| 中文字幕在线日韩91| 制服丝袜亚洲| 天堂在线亚洲| 亚洲第一视频网| 成人免费视频一区| 成年网址网站在线观看| 青草精品视频| 亚洲成a人片77777在线播放| 免费全部高H视频无码无遮掩| 亚洲狠狠婷婷综合久久久久| 国产97视频在线| 久久久久青草线综合超碰| 亚洲国产综合精品一区| 92精品国产自产在线观看| 午夜福利网址| 亚洲精品福利视频| 999精品视频在线| 色AV色 综合网站| 国产精品亚洲一区二区在线观看| 日本精品αv中文字幕| 国产91高清视频| 久久国产黑丝袜视频| 国产成人久久综合777777麻豆| 日本不卡视频在线| 亚洲手机在线| 中文毛片无遮挡播放免费| 国产白浆在线| 亚洲免费人成影院| 欧美人与牲动交a欧美精品| 亚洲第一综合天堂另类专| 制服丝袜在线视频香蕉| 第一页亚洲| 一级做a爰片久久免费| 美女视频黄又黄又免费高清| 色亚洲激情综合精品无码视频 | 曰韩人妻一区二区三区| 日韩无码黄色| 国产剧情无码视频在线观看| 国产乱子伦精品视频| 国产日韩精品一区在线不卡| 一级不卡毛片| 成人午夜亚洲影视在线观看| 五月婷婷综合网| 精品国产成人高清在线| 国产第一福利影院| 欧美日韩在线成人| 亚洲嫩模喷白浆| 亚洲欧洲日韩久久狠狠爱| 成人午夜精品一级毛片| 国产一区二区人大臿蕉香蕉| 国产乱码精品一区二区三区中文| 男人天堂伊人网| 中文字幕调教一区二区视频| 成人精品视频一区二区在线| 很黄的网站在线观看| 丁香六月激情综合| 日本免费a视频| 亚洲一区无码在线| 免费人成视网站在线不卡| 国产精品三区四区| 亚洲欧美国产视频| 欧美精品亚洲精品日韩专区| 亚洲日韩高清在线亚洲专区| 91福利在线看| 天天躁狠狠躁| 欧美日韩国产综合视频在线观看| 婷婷久久综合九色综合88| 91外围女在线观看| 欧美一级在线播放| 亚洲va精品中文字幕| 精品久久久久久成人AV| 国内丰满少妇猛烈精品播| 日韩一区二区三免费高清 | 亚洲成人一区二区| 亚洲91精品视频| 久久久久中文字幕精品视频| 中国毛片网| 日本黄色a视频| 日本午夜在线视频| 色婷婷电影网|