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

異種鋼焊接接頭蠕變過程的有限元模擬

2015-10-29 02:27:54張伯奇蔡志鵬李克儉李軼非潘際鑾
中國機械工程 2015年2期
關鍵詞:焊縫有限元結構

張伯奇 蔡志鵬 李克儉 李軼非 潘際鑾

清華大學,北京,100084

異種鋼焊接接頭蠕變過程的有限元模擬

張伯奇蔡志鵬李克儉李軼非潘際鑾

清華大學,北京,100084

異種鋼焊接接頭常應用于高溫下工作的大型構件,其蠕變特性和持久性能對結構的安全性具有重要意義。常規的蠕變試驗和高溫持久試驗僅適用于均勻材料,對非均勻的異種鋼接頭存在局限性。為了解決該問題,針對某種構件上采用的含有過渡層的馬氏體鋼與珠光體鋼異種金屬接頭進行了有限元建模。分別利用K-R蠕變損傷模型和改進的θ Projection模型進行模擬計算,并且對比了非均勻結構與均勻的焊縫、過渡層和母材的蠕變特性。對比發現,非均勻結構中抗蠕變性能較弱的過渡層部分的蠕變速率要大于均勻的過渡層材料的蠕變速率,這說明材料的不均勻性會帶來附加的蠕變損傷。因此,僅通過均勻材料的蠕變試驗來預測非均勻結構的壽命是不可靠的,應采用蠕變試驗與有限元計算相結合的方法進行非均勻結構的壽命預測。

異種鋼焊接;非均勻材料;蠕變模擬;有限元模擬

0 引言

異種鋼焊接目前被廣泛應用于核電、火電、石油化工等裝備制造領域。異種鋼焊接技術可以充分發揮不同材料的優勢,降低成本,并提高結構設計的靈活性,是大型裝備制造中的關鍵技術之一[1-2]。

異種鋼焊接接頭大多工作在高溫環境下,其蠕變特性和持久性能對結構的服役安全性至關重要。

對于異種鋼焊接接頭,常規的蠕變和持久試驗存在難以克服的局限性。首先,蠕變和持久試驗均是針對均勻材料而設計的,難以反映異種材料之間的相互作用;然后,蠕變和持久試驗均采用單軸拉伸的加載方式,而非均勻結構必然存在復雜的多軸應力狀態[3-4],與試驗假設并不相符;最后,蠕變和持久試驗周期長、成本高,嚴重制約了整個產品的設計和生產。因此,采用有限元計算方法模擬異種鋼接頭的蠕變過程并預測其壽命是行之有效的措施。

對耐熱鋼蠕變過程的模擬已有很多研究。Kachanov[5]提出了經典的K-R蠕變損傷模型;文獻[6-7]將K-R蠕變模型推廣到多軸應力條件下;劉學等[8]擬合了P91鋼的蠕變損傷模型,并進行了多軸應力條件下的有限元計算。還有很多學者從組織劣化方面研究了蠕變損傷模型的物理本質[9]。除了上述對均勻材料蠕變模型的研究以外,還有一些學者對焊接接頭非均勻結構的蠕變行為進行了研究。張建強等[4,10]對T91鋼和HR3C鋼焊接接頭的蠕變失效進行了有限元模擬,但只采用了簡單的冪律蠕變模型,而未對接頭的損傷情況進行表征;Sunil-Goyal等[11]利用有限元方法對2.25Cr1Mo鋼熱影響區的IV型裂紋進行了預測,同樣只采用了冪律蠕變模型,而未反映加速蠕變階段應力應變的劇烈變化??梢?,對焊接接頭,特別是對異種鋼接頭這種復雜的非均勻結構進行蠕變模擬還有待進一步研究。

本文采用K-R蠕變損傷模型和θ Projection模型分別研究含有2.25Cr1Mo過渡層的馬氏體鋼與珠光體鋼焊接接頭的蠕變行為,以期探索一種合理的異種鋼接頭壽命預測方法。

1 異種鋼接頭結構與蠕變試驗

為了解決異種鋼母材回火參數的矛盾,本文研究的異種鋼接頭采用了帶有過渡層的結構。首先在馬氏體鋼母材一側用埋弧焊堆焊若干層2.25Cr1Mo過渡層,焊后進行(670~690 ℃)×12 h回火,并將過渡層加工成坡口形狀。然后用TIG打底焊連接帶有過渡層的馬氏體鋼與珠光體鋼,并用2.25Cr1Mo焊絲埋弧焊填充。最后進行(650~660 ℃)×20 h回火。焊接接頭的結構如圖1所示。圖2所示為經過硝酸酒精腐蝕后的焊接接頭的照片。焊接接頭的母材和焊絲的主要成分見表1。

圖1 異種鋼焊接接頭結構示意圖

圖2 焊接接頭照片

為了確定接頭持久性能的薄弱環節,對整個接頭進行了高溫持久試驗。分別采用提高溫度和增大應力的加速試驗方法進行試驗,試驗溫度在500~620 ℃之間,試驗應力在80~180 MPa之間。試樣橫跨接頭的各個區域,通過斷裂位置即可判斷其薄弱部位。試驗發現,所有試樣均在過渡層內斷裂,如圖3所示。

表1 焊接接頭各部分的主要成分(質量分數) %

圖3 接頭高溫持久試樣金相照片

可見過渡層是接頭持久性能的薄弱環節。因此,在有限元計算中著重研究過渡層的蠕變行為,模型包括馬氏體鋼母材、過渡層和焊縫三部分,而省略了珠光體鋼一側的建模。

過渡層是多層結構,由于對高合金的馬氏體鋼母材逐層稀釋,故各層過渡層的成分有所不同。但是每層過渡層厚度都很小,無法單獨制取蠕變試樣,只能忽略各過渡層之間在蠕變特性上的差別。過渡層材料的不均勻性將通過彈性模量的差別來體現。

在600 ℃、不同應力下分別對馬氏體鋼母材、過渡層和焊縫進行蠕變試驗,得到其蠕變曲線。蠕變試樣直徑為10 mm,小于焊縫和過渡層的寬度,試樣標距為55 mm。以馬氏體鋼為例,測得其蠕變速率與時間的關系如圖4所示。

圖4 600 ℃下馬氏體鋼母材的蠕變試驗數據

2 蠕變模型的擬合

2.1K-R蠕變損傷模型的擬合

蠕變過程一般分為蠕變硬化、穩態蠕變和加速蠕變3個階段。其中,前2個階段可用冪律模型很好地描述,而第3階段通常采用K-R蠕變損傷模型來描述。

K-R蠕變損傷模型應用廣泛,它不僅可以計算加速蠕變過程中構件的應力應變情況,還能預測構件的損傷部位。K-R蠕變損傷模型可用2個方程來表述,即與損傷變量耦合的蠕變速率方程和損傷演化方程,其形式如下:

(1)

(2)

式中,ε為應變;σ為應力;D為損傷變量(初值為0,單調增大,變為1時表示完全損傷);A、B、n、p、φ為常數,由材料性質決定。

在應力恒定的條件下,對式(2)積分:

(3)

式中,t為時間。

則有

(4)

將式(4)代入式(1),得

(5)

將蠕變硬化和穩態蠕變階段的應變率表示為冪律形式,得到完整的應變率方程:

(6)

式中,q、C、m為系數。

利用蠕變試驗獲得的數據,采用非線性擬合的方法對式(6)進行擬合,分別得到馬氏體鋼母材、過渡層和焊縫的蠕變模型常數。該蠕變模型可以獲得良好的擬合度,以馬氏體鋼母材在600 ℃、100 MPa下的蠕變速率曲線為例,其擬合情況如圖5所示。

圖5 采用K-R蠕變損傷模型擬合的馬氏體鋼蠕變速率曲線

2.2θ Projection模型的擬合與改進

θ Projection模型將蠕變過程視為蠕變硬化和加速蠕變2個階段,分別采用指數函數來描述,其表達式為

ε=θ1(1-e-θ2t)+θ3(eθ4t-1)

(7)

lgθi=ai+biσ+ciT+di(σ T)i=1,2,3,4

(8)

式中,ai、bi、ci、di為系數。

式(7)為蠕變曲線方程,即應變ε與時間t之間的關系。其中,θ1、θ2、θ3、θ4為參數,它們都可表示為溫度T和應力σ的函數,如式(8)所示。若溫度不變,式(8)僅為應力的函數,則最少只需測得2條不同應力下的蠕變曲線即可得到全部常數。

采用θ Projection模型擬合的馬氏體鋼母材100 MPa下的蠕變曲線和蠕變速率曲線如圖6所示。

(a)蠕變曲線

(b)蠕變速率曲線圖6 θ Projection模型對馬氏體鋼蠕變數據的擬合結果

由圖6可見θ Projection模型的擬合效果并不理想,主要原因是蠕變硬化階段的蠕變速率與試驗數據差別過大,而且最小蠕變速率也遠低于實際值。由圖6b可見,蠕變硬化階段的蠕變速率在對數坐標下近似呈直線,應符合冪函數形式,因此蠕變曲線在硬化階段也應符合冪函數形式??紤]到θ Projection模型對加速蠕變階段的擬合度較好,因此保留θ Projection模型的加速蠕變項,而將硬化階段改為冪函數形式。改進后的θ Projection模型表達式為

ε=θ1σθ2tθ3+θ4(eθ5t-1)

(9)

其中,θ1、θ2、θ3為材料常數,θ4、θ5為溫度和應力的函數,表達式見式(8)。用改進后的θ Projection模型擬合馬氏體鋼母材在600 ℃、100 MPa下的蠕變曲線和蠕變速率曲線,并與改進前的擬合結果進行對比,如圖7所示。

(a)蠕變曲線

(b)蠕變速率曲線圖7 改進前后的θ Projection模型對蠕變數據擬合結果對比

由圖7可見,改進后的θ Projection模型在蠕變硬化階段具有良好的擬合精度,而且最小蠕變速率與實際值很接近,蠕變加速階段蠕變速率的陡升特性也與實際相符,只是加速階段開始時間比試驗結果略晚。

與K-R蠕變損傷模型相比,改進后的θ Projection蠕變模型不僅同樣獲得了良好的擬合度,而且保持了簡潔的表達式形式,對有限元計算比較有利。但是K-R蠕變損傷模型可以預測構件發生蠕變損傷的部位和程度,仍具有一定優勢。下面分別采用以上兩種模型對焊接接頭的蠕變行為進行模擬。

3 蠕變過程的有限元模擬

3.1材料彈性模量的確定

除上文中擬合得到的蠕變模型外,進行彈性范圍內的蠕變模擬還需要知道材料的彈性模量和泊松比。泊松比通常取0.3,而彈性模量的測試存在困難。過渡層成分梯度很大,尺寸又很小,常規拉伸試驗只能反映過渡層的平均彈性模量,而無法體現材料的不均勻性。因此,本文采用通過顯微硬度換算彈性模量的方法。Bao等[12]研究了硬度測量過程中的外力做功、能量耗散與材料彈性恢復,提出了彈性模量與硬度之間的理論關系:

(10)

其中,H為硬度,F為常數,Rs為恢復阻力。對于Berkovich或Vickers壓頭,F取0.6647。Er為接觸模量,其表達式為

(11)

其中,E、μ分別為被測材料的彈性模量和泊松比,Ei、μi分別為壓頭的彈性模量和泊松比。式(10)中的Rs反映的是材料在加載和卸載過程中的能量耗散,其表達式為

(12)

式中,pm為加載時的最大載荷;hs為壓痕邊緣線在載荷方向上的彈性恢復位移。

由于本文研究的焊接接頭各部位均屬于CrMoV系耐熱鋼,具有相似的彈塑性特征,可以假設其具有相同的恢復阻力。因此,根據式(10),接觸模量與硬度的平方根成正比。利用對整個接頭的顯微硬度測試結果,以及對均勻母材的彈性模量測試結果,可以換算出不同部位的接觸模量,再通過式(11)推算出各部位的彈性模量,如圖8所示。

圖8 接頭硬度分布與彈性模量的計算

圖8中包含馬氏體鋼母材、焊縫和3層過渡層的顯微硬度測試數據。根據硬度的變化趨勢,將接頭離散成31層,每層根據其平均硬度換算彈性模量。這樣可以完整地表現出各層過渡層熔合區的軟化層、熱影響區的淬硬區、正火區、回火軟化區等區域的不均勻性。

3.2焊接接頭的有限元建模

典型的焊接接頭高溫持久試樣斷口附近的解剖金相照片如圖3所示。以圖3中的過渡層結構為依據在ABAQUS軟件中進行建模,采用平面應變模型,在一側施加120 MPa拉伸載荷,如圖9所示。

圖9 焊接接頭幾何模型、載荷及邊界條件

圖9中包含馬氏體鋼母材、焊縫以及交錯排列的3層過渡層。根據彈性模量的變化,在各層熔合區附近的軟化層部位細分了多層結構,并減小網格尺寸,以更好地體現材料的非均勻性,如圖10所示。

圖10 焊接接頭模型的網格劃分

3.3焊接接頭初始狀態應力應變分布

在拉伸載荷作用下,焊接接頭初始狀態的應力應變分布如圖11所示。

由圖11可見,初始狀態下硬度和彈性模量較低的過渡層區域整體應變較大,而且過渡層之間的熔合區軟化層存在應變集中。另外,各層過渡層的淬硬區存在應力集中,特別是焊道交錯排列形成的尖角處應力集中現象更加明顯。

(a)Mises等效應力

(b)橫向正應變圖11 焊接接頭初始狀態應力應變分布

3.4兩種蠕變模型的對比

經過1000 h蠕變后,兩種蠕變模型的計算結果如圖12所示。本文僅關注應變分布情況,因為蠕變失效往往是由應變決定的。

(a)蠕變損傷模型得到的橫向正應變

(b)蠕變損傷模型的損傷變量

(c)改進型θ Projection模型得到的橫向正應變圖12 經1000 h蠕變后兩種模型的計算結果(變形量放大15倍)

由圖12可見,兩種蠕變模型得到的應變分布類似,都是在較軟的過渡層部分蠕變變形較大,并出現縮頸現象,但是θ Projection模型的變形程度更大。另外,過渡層熔合區軟化層在初始狀態下存在的應變集中現象,經過蠕變后得到緩解。這是因為蠕變應變逐漸增大,比初始狀態的應變大2個數量級,因此軟化層的應變集中現象被掩蓋。

觀察圖3所示的實際斷口,發現斷裂發生在過渡層之間的熔合區,并在熔合區出現正在擴展的裂紋。由于有限元計算中沒有引入裂紋的因素,只能反映出過渡層危險性較高。實際上在蠕變的開始階段,應變集中的軟化層可能出現微裂紋,成為過渡層中的薄弱環節。因此,兩種模型對實際結構蠕變行為的預測都是合理的,但仍存在一定局限性。

以過渡層中最危險部位的應變為例,對比兩種模型的差別,如圖13所示。

圖13 兩種模型對過渡層中最危險部位應變的計算結果對比

由圖13可見,改進型θ Projection模型的蠕變加速階段出現時間較早,但蠕變率的增長比較平緩,而蠕變損傷模型的加速蠕變階段出現較晚,但蠕變率突然快速增長。結合損傷變量的分析發現,采用蠕變損傷模型對加速蠕變開始時間的預測更為精確,而改進型θ Projection模型對應變的預測比較保守,安全性較高。

3.5材料非均勻性對蠕變行為的影響

非均勻結構的蠕變行為不同于均勻材料,這是因為結構各部分的蠕變特性不同,變形量存在差異,因此相互之間會產生拘束作用,出現復雜的多軸應力狀態。為了研究材料非均勻性對結構蠕變行為的影響,將相同載荷下焊接接頭中的焊縫內最危險部位的蠕變曲線與均勻焊縫材料的蠕變曲線進行對比,如圖14所示。

圖14 接頭中的焊縫與均勻焊縫材料蠕變曲線對比

由圖14可見,兩種模型的計算結果都顯示,在焊接接頭中的焊縫的蠕變量要大于均勻的焊縫材料的蠕變量。這說明材料的非均勻性會使結構整體的持久性能比其中最弱的單一材料更弱。

4 結論

(1)K-R蠕變損傷模型和改進型θ Projection模型對馬氏體鋼與珠光體鋼異種接頭中各部分材料的蠕變數據都能獲得良好的擬合結果。

(2)在非均勻結構的蠕變模擬中,蠕變硬化階段兩種模型計算結果基本相同。改進型θ Projection模型更早進入加速蠕變階段,但應變速率增長比較平緩;K-R蠕變損傷模型加速蠕變階段開始較晚,但應變速率快速增長。K-R模型對加速蠕變階段開始時間的預測比較精確,而改進型θ Projection模型安全性更高。

(3)兩種模型的計算結果都表明,非均勻結構整體的持久性能要比其中最弱的單一材料更弱。因此,用各部分材料分別做蠕變試驗來預測結構整體的蠕變特性是不可靠的,應采用蠕變試驗與有限元計算相結合的方法進行壽命預測。

[1]馬琰, 陳蘭.異種金屬焊接技術的應用實踐[J].中國新技術新產品, 2011(19): 159.

Ma Yan, Chen Lan.Applications of Dissimilar Joint Welding Technologies[J].China New Technologies and Products, 2011(19): 159.

[2]Lundin C D.Dissimilar Metal Welds-transition Joints Literature Review[J].Welding Journal, 1982, 61(S2): 58-63.

[3]史春元, 范英, 陳字剛,等.模擬異種鋼接頭蠕變力學特性的數值分析[J].大連鐵道學院學報,1994, 15(4): 84-88.

Shi Chunyuan, Fan Ying, Chen Zigang, et al. Numerical Analysis of Creep Mechanical Characteristics of Simulated Dissimilar Metal Welded Joints[J].Journal of Dalian Railway Institute, 1994, 15(4): 84-88.

[4]張建強, 郭嘉琳, 張國棟, 等.馬氏體/珠光體異種耐熱鋼焊接接頭蠕變失效數值模擬[J].焊接學報,2009,30(2):87-90.

Zhang Jianqiang, Guo Jialin, Zhang Guodong, et al. Numerical Simulation on Creep Failure of Dissimilar Welded Joint between Martensitic Heat-resistant Steel and Pearlitic Heat-resistant Steel[J].Transactions of the China Welding Institution, 2009, 30(2): 87-90.

[5]Kachanov L M.On Time to Rupture in Creep Condition (in Russian)[J].Izviiestia Akademii Nauk SSSR, Otdelenie Tekhnicheskikh Nauk, 1958,8(1):26-31.[6]Hayhurst D R. Creep Rupture under Multi-axial States of Stress[J].J. Mech. Phys. Solids, 1972, 20(6): 381-390.

[7]Johnson A E. Complex-stress Creep of Metals[J].Metallurg. Rev., 1960, 5(20): 447-506.

[8]劉學, 段向兵, 閆平, 等. P91鋼高溫蠕變的數值研究[J]. 電力建設, 2009, 30(5): 52-55.

Liu Xue, Duan Xiangbing, Yan Pin,et al.Numerical Study of P91 Steel High Temperature Creep[J].Electric Power Construction, 2009, 30(5): 52-55.

[9]張俊善.材料的高溫變形與斷裂[M].北京: 科學出版社, 2007.

[10]張建強, 張國棟, 姚兵印, 等.T91/HR3C異種耐熱鋼焊接接頭蠕變失效有限元模擬[J].焊接學報, 2013, 34(7): 68-72.

Zhang Jianqiang, Zhang Guodong, Yao Bingyin, et al. Numerical Simulation on Creep Failure of T91/HR3C Heat-resistant Steels Welded Joint[J].Transactions of the China Welding Institution, 2013, 34(7): 68-72.

[11]Sunil-Goyal, Laha K, ChandravathiK S, et al.Prediction of Type IV Cracking Behaviour of 2.25Cr-1Mo Steel Weld Joint Based on Finite Element Analysis[J].Transaction of the Indian Institute of Metals, 2010, 63(2/3): 461-466.

[12]Bao Y W, Wang W, Zhou Y C.Investigation of the Relationship between Elastic Modulus and Hardness Based on Depth-sensing Indentation Measurements[J].Acta Materialia, 2004, 52(18): 5397-5404.

(編輯陳勇)

Finite Element Analysis of Creep Behavior of Dissimilar Steel Welded Joint

Zhang BoqiCai ZhipengLi KejianLi YifeiPan Jiluan

Tsinghua University,Beijing,100084

Dissimilar welded joints are often used in large components working under high temperature. The creep rupture property of dissimilar welded joint is of great significance for the safety of the structure. Conventional creep test and creep rupture test are only applicable to homogeneous materials, but not to inhomogeneous materials. In order to solve this problem, a finite element model was built according to a dissimilar welded joint used in a real component which made up of martensitic steel and pearlitic steel.By using K-R creep damage model and advanced θ Projection model respectively, it is found that the creep rupture property of inhomogeneous dissimilar joint is different from that of homogeneous base metal, weld metal and transition layer. The creep rate of transition layer in dissimilar joint is higher than the same homogeneous material, which means the heterogeneity of structure leads to additional creep damage. The conclusion can be made that it is unreliable to predict the creep life of a dissimilar steel structure only according to creep rupture tests on different parts for homogenous materials,the finite elements analysis based on creep rupture test is also necessary.

dissimilar steel welding; heterogeneous material; creep simulation; finite element analysis

2014-01-22

TG407DOI:10.3969/j.issn.1004-132X.2015.02.025

張伯奇,男,1986年生。清華大學機械工程學院博士研究生。主要研究方向為焊接冶金等。發表論文6篇,獲國家專利2項。蔡志鵬,男,1974年生。清華大學機械工程學院副研究員。李克儉,男,1989年生。清華大學機械工程學院博士研究生。李軼非,男,1989年生。清華大學機械工程學院博士研究生。潘際鑾,男,1927年生。中國科學院院士,清華大學機械工程學院教授、博士研究生導師。

猜你喜歡
焊縫有限元結構
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
基于焊縫余高對超聲波探傷的影響分析
TP347制氫轉油線焊縫裂紋返修
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
機器人在輪輞焊縫打磨工藝中的應用
論《日出》的結構
光譜分析在檢驗焊縫缺陷中的應用
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 午夜日本永久乱码免费播放片| 九九九精品视频| 国产99免费视频| 99激情网| 国产婬乱a一级毛片多女| 欧美精品1区| 国产精品女熟高潮视频| 在线播放91| 国产乱子精品一区二区在线观看| 制服丝袜国产精品| 亚洲欧洲日韩综合| 免费看美女毛片| 五月激情综合网| 亚洲人成网站18禁动漫无码| 欧洲一区二区三区无码| 亚洲 欧美 偷自乱 图片 | 欧美精品一区在线看| 8090成人午夜精品| 国产成人精品一区二区| 免费观看精品视频999| 久久狠狠色噜噜狠狠狠狠97视色 | 国产日韩精品欧美一区喷| 97国产精品视频自在拍| 亚洲视频欧美不卡| 伊人久久综在合线亚洲91| 国产97公开成人免费视频| 精品91自产拍在线| 久久黄色影院| 国产精品主播| 国产幂在线无码精品| 伊人无码视屏| 久久国产乱子伦视频无卡顿| 国产精品v欧美| 久久精品无码中文字幕| 亚洲高清在线天堂精品| 亚洲国产欧美目韩成人综合| 自拍亚洲欧美精品| 国产一级视频在线观看网站| 国内视频精品| 国产91精品最新在线播放| 成人91在线| 国产成人欧美| 四虎永久免费网站| 国产喷水视频| 伊人91在线| 国产精品不卡片视频免费观看| h视频在线播放| 亚洲最大福利网站| 欧美一区二区丝袜高跟鞋| 日韩麻豆小视频| 婷婷色一区二区三区| 99在线视频网站| 毛片网站免费在线观看| 国产免费久久精品99re丫丫一| 欧美专区在线观看| 欧美特黄一级大黄录像| 看国产一级毛片| 国产精品香蕉| 亚洲无码37.| 国产99久久亚洲综合精品西瓜tv| 国产精品毛片一区| 精品久久久久成人码免费动漫| 国产在线观看91精品亚瑟| 熟妇丰满人妻av无码区| 99久久精品国产自免费| 五月婷婷亚洲综合| 国产精品亚洲va在线观看| 亚洲男人在线| 国产成人a在线观看视频| 亚洲va在线∨a天堂va欧美va| 丰满人妻被猛烈进入无码| 欧美日本在线观看| 欧美福利在线| 亚洲精品无码AⅤ片青青在线观看| 成人日韩精品| 亚洲中文无码h在线观看 | 天堂成人av| 国产精品人莉莉成在线播放| 国产精品久久久久久久久久98| 99无码中文字幕视频| 亚洲欧美成aⅴ人在线观看| 91麻豆精品国产高清在线|