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

三輥非對稱滾彎成形的數值仿真及試驗驗證

2016-09-13 06:58:43胡捷飛許光輝李曙生
中國機械工程 2016年15期
關鍵詞:有限元模型

王 艷 胡捷飛 許光輝 李曙生 熊 巍

1.上海理工大學,上海,200093  2.鎮江市鍛壓機床廠,鎮江,2120003.泰州職業技術學院,泰州,225300

?

三輥非對稱滾彎成形的數值仿真及試驗驗證

王艷1,2胡捷飛1許光輝2李曙生3熊巍1

1.上海理工大學,上海,2000932.鎮江市鍛壓機床廠,鎮江,2120003.泰州職業技術學院,泰州,225300

基于有限元法對三輥非對稱滾彎成形工藝進行了研究,對比分析了三輥非對稱滾彎成形和三輥對稱滾彎成形過程中,變形區應力場、板材上表面的塑形應變場及卷制力的變化規律。仿真結果表明:側輥位移進給量相同的工況下,三輥非對稱滾彎成形的卷制力大于三輥對稱滾彎成形的卷制力;三輥非對稱滾彎變形區的縱向應力和徑向應力均大于三輥對稱滾彎成形的縱向應力和徑向應力;三輥非對稱滾彎成形板材壓彎段的成形質量高于三輥對稱滾彎成形的成形質量。最后,經三輥非對稱滾彎試驗驗證,有限元模型的成形誤差為6.8%,有較高的精度。

三輥非對稱滾彎;應力場;卷制力;成形精度

0 引言

滾彎成形工藝是一種無切削、高效率、高質量的生產工藝,相比傳統的板材加工工藝, 它具有高效、省材、低成本、精度高等優點,因此,滾彎成形工藝在造船、石油化工、鍋爐、壓力容器等行業已經得到廣泛的應用[1]。三輥滾彎成形工藝是工業應用最為廣泛的一種滾彎成形工藝,根據加工時三個軋輥間的相對位置分布,它又可分為對稱滾彎成形和非對稱滾彎成形兩種。與三輥對稱滾彎成形工藝相比,三輥非對稱滾彎成形工藝一次上料即可完成筒體卷制加工,并能有效減少成形筒體的剩余直邊,可進一步提高生產效率。

滾彎成形工藝涉及軋輥與板材間的摩擦、接觸以及板材的回彈,是一個復雜的成形過程[2]。板材滾彎成形工藝是設計卷板裝備、制定彎卷工藝、提高卷板精度和實現彎卷自動控制的理論基礎,具有重要研究價值和現實意義[3]。目前研究者對于三輥非對稱滾彎成形工藝的研究主要是分析工藝參數對三輥非對稱滾彎成形的影響和三輥非對稱滾彎工藝模型的建立等方面。Tran等[4]分析了側輥位移進給量、板材尺寸等工藝參數對三輥非對稱滾彎成形的影響,并研究了滾彎過程中板材表面塑性應變的演變規律。Feng等[5]針對板材的三輥非對稱滾彎工藝建立了側輥位移進給量與板材最終成形半徑的滾彎工藝模型,并結合試驗對工藝模型進行了驗證。Zhang等[6]在考慮板材外形尺寸、工作輥直徑、板材與工作輥摩擦等因素的基礎上,對中厚板彈塑性三輥非對稱壓彎工藝過程進行了工藝模型的建立與試驗驗證。但是研究者對于三輥非對稱滾彎成形過程中,板材變形區內部應力應變狀態、板材的回彈、塑性變形、卷制力等方面卻少有研究,也未能進一步結合三輥對稱滾彎成形工藝對三輥非對稱滾彎成形工藝的特點進行對比分析。

本文針對上述問題,首先應用有限元軟件ABAQUS,分別建立三輥對稱滾彎成形和三輥非對稱滾彎成形的二維動態有限元模型。其次通過對比分析對稱滾彎和非對稱滾彎成形過程中,變形區應力的分布特點、卷制力的變化規律以及板材上表面塑性應變場的分布規律,對三輥非對稱滾彎成形工藝的機理及特點進行研究。最后結合三輥非對稱滾彎成形試驗,對有限元模型的成形精度進行驗證。

1 三輥非對稱滾彎工藝模型

三輥非對稱滾彎時,通過調整側輥相對于上輥的位置來控制板材滾彎成形后的成形半徑。在上輥和下輥對變形區板材的夾持作用下,三輥非對稱滾彎變形區近似關于上輥與下輥輥心連線對稱,并假設板材在變形區內成形曲率一致。基于上述假設,根據側輥相對位置與成形后板材的幾何關系,即可推導出側輥位移進給量與板材成形半徑的數值關系。如圖1所示,非對稱滾彎時,側輥沿傾角方向的位移進給量為d,板材最終成形半徑為Rf,此時側輥輥心相對坐標為(x0,y0),距上輥輥心距離為a。圖1中,A表示側輥傾斜角與上下輥輥心連線的交點;ψ為側輥傾角。

圖1 側輥與板材成形半徑的幾何關系示意圖

假設非對稱滾彎成形時,塑性變形后的板材與直邊部分始終相切,則根據圖1所示的幾何關系,推導出側輥輥心相對坐標(xs,ys)與板材回彈前成形半徑R的關系式[7]:

(1)

式中,Rb為下輥半徑;θ為側輥輥心與上輥輥心連線與垂線的夾角;T為板厚。

三輥非對稱滾彎成形時,側輥輥心初始位置相對坐標(x0,y0)及側輥輥心相對位置坐標(xs,ys)滿足下式:

(2)

側輥的位移進給量d為

(3)

板材回彈后的成形半徑才是最終成形半徑Rf。回彈與板材的彈性模量、屈服極限、板材厚度、板材自重、板與軋輥間的摩擦等因素有關,因此,難以建立高精度的回彈模型。工程上一般采用的板材回彈模型[8],經變換后為

(4)

式中,K0、K1分別為相對強化系數和形狀系數(矩形截面取值為1.5);E為彈性模量;σs為屈服極限。

聯立式(1)~式(4)即可確定側輥位移進給量d與非對稱滾彎后板材成形半徑Rf間的數值關系。

2 有限元模型的建立

2.1幾何模型

本文采用ABAQUS/Explicit模塊的顯式中心差分算法對三輥非對稱滾彎成形和三輥對稱滾彎成形過程進行動態仿真。分別建立三輥非對稱滾彎成形和三輥對稱滾彎成形的二維有限元模型,如圖2所示。軋輥在滾彎過程中的彈性變形量很小,將其設定為剛體。板材厚度為30 mm,選用可變形的平面應力應變單元CPS4R,并沿厚度方向將板材網格劃分為6層。有限元模型中軋輥直徑參數見表1、表2。滾彎加工時上輥為主動輥,轉速為0.3 rad/s,下輥及側輥為從動輥。三輥非對稱滾彎成形加工時,側輥沿30°傾角傾斜向上進給;三輥對稱滾彎成形時,下輥也沿30°傾角傾斜向上進給,如圖2所示。

表1 三輥非對稱有限元模型軋輥直徑

表2 三輥對稱有限元模型軋輥直徑

2.2接觸定義

(a)非對稱滾彎有限元模型

(b)對稱滾彎有限元模型圖2 滾彎成形有限元模型

滾彎成形涉及復雜的接觸過程,對于板材與軋輥間的接觸,采用surface to surface接觸,并將剛體軋輥表面定義為主面,板材表面定義為從面。采用罰函數摩擦模型求解軋輥與板材間的接觸,其中三輥非對稱滾彎模型中,為減小側輥對板材的摩擦阻力,保證板材的順利送料[9],將上下輥與板材間的摩擦因數設為0.15[10],將側輥與板材間的摩擦因數設為0.1;在對稱滾彎模型中,將上輥與板材間的摩擦因數分別設為0.15,而下輥與板材間的摩擦因數設為0.1。

2.3材料模型

仿真試驗選用的鋼板材料為Q235B鋼板,材料屬性[6]如下:屈服應力σs為235MPa;彈性模量E為210GPa;泊松比μ為0.3;密度ρ為7800kg/m3;強化系數K為1900。假設板材材料為各向同性,且不考慮板材自重對滾彎成形的影響。

滾彎成形過程中板材發生彈塑性變形,必須在ABAQUS中設置板材材料的本構關系。因雙線性強化模型能夠較好地擬合材料的真實本構關系[11],本文采用雙線性強化模型定義板材的本構關系,以模擬板材在滾彎成形過程中發生的塑性變形。本構方程為

(5)

式中,σ為應力;ε為應變;εe為彈性應變極限。

3 仿真結果分析

3.1三輥非對稱滾彎成形工藝過程

三輥非對稱滾彎仿真模擬過程如圖3所示,與三輥對稱滾彎成形過程相同[11],可將三輥非對稱滾彎成形過程分為三個階段:①壓彎階段。此時側輥沿傾角斜向進給,對板材端部進行壓彎。②過渡滾彎階段。上輥轉動對板材端部進行滾彎。③穩態滾彎階段。上輥繼續轉動對板材進行連續的滾彎。過渡滾彎階段是壓彎階段到穩態滾彎階段的一個過渡階段,歷時較短,因此,可將滾彎成形后的板材沿長度方向主要分為壓彎段和穩態滾彎段,而穩態滾彎段筒體的成形半徑即為筒體的整體半徑。

(a)壓彎階段

(b)過渡滾彎階段

(c)穩態滾彎階段圖3 非對稱滾彎模擬過程

3.2變形區應力場演變分析

圖4所示為滾彎成形過程中變形區厚度方向截面應力分布曲線。由圖4a,在側輥位移進給量相同的工況下,三輥非對稱滾彎成形時,變形區內縱向應力和徑向應力均大于三輥對稱滾彎成形的縱向應力和徑向應力。對于三輥對稱滾彎成形,其變形區內截面縱向應力遠大于徑向應力,因此,變形區內徑向應力對三輥對稱滾彎成形的影響可忽略。但是由于上下輥的夾持作用,三輥非對稱滾彎變形區內徑向應力要明顯大于對稱滾彎成形的徑向應力,因此,必須考慮徑向應力對非對稱滾彎成形的影響。板材的回彈是一個彈性應變能釋放的過程,滾彎成形后板材的殘余應力會隨著彈性應變能釋放量的增加而減小。圖4b所示為回彈后板材厚度方向截面應力分布曲線。由圖4b,三輥非對稱滾彎成形后,板材截面應力要小于對稱滾彎成形的應力,即非對稱滾彎成形后板材的殘余應力比對稱滾彎成形的殘余應力小。在滾彎變形區內,板材受力狀態為上壓下拉,這種受力狀態會增大板材回彈量[12]。非對稱滾彎成形時,變形區截面縱向應力大于對稱滾彎成形的縱向應力,因此,其回彈量也會相應增大,而隨著板材回彈量的增大,殘余應力也會隨之減小。

(a)變形區截面應力分布曲線

(b)回彈后殘余應力分布曲線圖4 截面厚度方向上應力演變

3.3板材周向應變場分布分析

滾彎成形后,板材上表面的塑性應變場分布體現了滾彎成形質量。圖5所示為滾彎成形后板材上表面的塑性應變場分布情況。由圖5可知,三輥對稱滾彎成形與三輥非對稱滾彎成形后,在壓彎段,板材均存在較大的殘余塑性應變。對稱滾彎成形后,壓彎段與穩態滾彎段的塑性應變差為0.013,非對稱滾彎成形的兩段間的塑性應變差為0.003。板材純彎曲時上表面應變與成形半徑的關系如下:

εr=1-r/rm

(6)

式中,εr為板材上表面應變;r為上表面成形半徑;rm為板材中心層成形半徑。

圖5 板材上表面塑性應變場分布

由式(6)可知,壓彎段與穩態滾彎段間的塑性應變差必然導致壓彎段與穩態滾彎段的成形半徑不一致,而兩段成形半徑的一致性是滾彎成形質量的重要體現。因此,三輥非對稱滾彎成形后,壓彎段與穩態滾彎段間成形半徑的一致性要高于對稱滾彎成形兩段成形半徑的一致性,即三輥非對稱滾彎成形后板材壓彎段成形質量較高。隨著滾彎過程由壓彎階段過渡至穩態滾彎階段,板材上表面的塑性應變波動也逐漸趨于穩定。在穩態滾彎階段,與三輥非對稱滾彎成形相比,對稱滾彎成形后板材上表面的塑性應變分布波動較小,分布更加均勻,說明對稱滾彎成形后板材穩態滾彎段的成形均勻度較高。經計算,在側輥位移進給量為50 mm的工況下,穩態滾彎時,三輥非對稱滾彎成形與對稱滾彎成形,上輥與板材間的接觸弧長分別為35 mm和20 mm。與板材軋制過程相似,滾彎成形過程中,接觸弧長的增大會增大接觸力,加劇軋輥的彈性變形,降低成形質量。而較小的接觸弧長有益于提高成形的精度,并減小卷制力。因此,在側輥位移進給量為50 mm的工況下,三輥非對稱滾彎成形會增大滾彎變形區板材與軋輥間的接觸面積,增大卷制力并加劇板材上表面的塑形應變波動,進而會降低滾彎成形質量。

3.4滾彎卷制力變化分析

滾彎成形過程中,側輥與板材的接觸力即為滾彎成形的卷制力。側輥位移進給量為50 mm時,滾彎加工過程中卷制力隨時間的變化曲線如圖6所示。由圖6,在壓彎階段,三輥非對稱滾彎與對三輥對稱滾彎的卷制力相近,并在壓彎階段結束時達到峰值142.5 kN。但隨著滾彎過程由壓彎階段過渡至穩態滾彎階段后,非對稱滾彎與對稱滾彎的卷制力均有所減小。在穩態滾彎階段,三輥非對稱滾彎成形所需的卷制力約為134.3 kN,而對稱滾彎成形的卷制力約為118.2 kN。在三輥對稱滾彎成形的壓彎階段,變形區近似關于上下輥輥心連線對稱,但隨著滾彎成形過程由壓彎階段向穩態滾彎階段過渡,滾彎變形區也發生了轉移[6]。而三輥非對稱滾彎成形過程中,上輥和下輥對板材的夾持,使得非對稱滾彎變形區的轉移量小于三輥對稱滾彎成形的轉移量,因此,三輥非對稱滾彎由壓彎向穩態滾彎階段過渡后,其卷制力的變化幅度會小于對稱滾彎成形時卷制力的變化幅度。

圖6 滾彎成形卷制力變化曲線

圖7所示為穩態滾彎階段的卷制力隨側輥位移進給量的變化曲線。由圖7,在側輥位移進給量相同的工況下,穩態滾彎階段,三輥非對稱滾彎成形的卷制力要大于三輥對稱滾彎成形的卷制力,且隨側輥位移進給量的增大,三輥非對稱滾彎成形與三輥對稱滾彎成形的卷制力的差值也逐漸增大。在側輥位移進給量從30 mm逐漸增加至70 mm的過程中,三輥非對稱滾彎成形與三輥對稱滾彎成形的卷制力的差值也逐漸從13.1 kN增加至27.2 kN。

圖7 不同側輥位移進給量時卷制力變化曲線

4 試驗驗證

筒件的成形半徑是滾彎變形區應力、應變及卷制力對滾彎成形的影響的直接體現,也是衡量滾彎成形質量的重要標準。本文通過三輥非對稱試驗,分析有限元數值仿真后筒件的成形半徑精度,即可判斷出有限元模型的準確性和可靠性。采用W1130×2500型三輥非對稱卷板機,對尺寸為30 mm×1500 mm×7500 mm的Q235B鋼板進行三輥非對稱滾彎加工試驗。驅動右輥沿側輥傾角方向傾斜向上進給50 mm對板材進行三輥非對稱滾彎加工。圖8為三輥非對稱卷板機滾彎加工后的筒件與有限元數值仿真后筒件的對比圖。

(a)試驗滾彎加工后筒件

(b)仿真試驗滾彎加工后筒件圖8 滾彎試驗與仿真結果對比

將滾彎成形后的筒件沿周向分為五段,采用數字半徑測量儀,對各段筒件的中截面內徑值進行測量。同時也將仿真成形后的筒件沿周向分為五段,并應用MATLAB數值擬合對各段板材成形半徑進行測量。將試驗實測成形半徑值與仿真計算成形半徑值列于表3。由表3可知,經三輥非對稱滾彎成形后的筒件平均半徑為1197.5 mm,在相同的加工工況下,有限元仿真得到的筒件平均半徑為1115.1 mm,有限元模型成形半徑誤差為6.8%。整體上分析,有限元模型成形半徑誤差不超過7.0%,有較高的精度,驗證了有限元模型的準確性和可靠性。

表3 滾彎試驗值、仿真模擬值對比

誤差產生的原因主要是有限元模型采用了二維建模,假設板材材料為各向同性,忽略板材自重及軋輥彈性變形對滾彎成形的影響。并且有限元模型采用了雙線性強化模型,相對于材料的真實本構模型存在一定誤差。同時有限元模型難以真實模擬出滾彎成形過程中板材的回彈,也會對有限元模型的精度產生影響,但其整體誤差不超過7.0%,有著較高的精度和可靠性,可滿足工程應用要求。

5 結論

(1)在側輥位移進給量相同的工況下,穩態滾彎階段,三輥非對稱滾彎成形變形區內,板材厚度方向截面的縱向應力和徑向應力均大于三輥對稱滾彎成形的縱向應力和徑向應力。

(2)在側輥位移進給量相同的工況下,穩態滾彎階段,三輥非對稱滾彎成形的卷制力大于三輥對稱滾彎成形的卷制力,且隨側輥位移進給量的增大,三輥非對稱滾彎成形與三輥對稱滾彎成形的卷制力的差值也逐漸增大。

(3)三輥非對稱滾彎成形后,板材壓彎段與穩態滾彎段的等效塑性應變差值小于三輥對稱滾彎成形的等效塑性應變差值,因而三輥非對稱滾彎成形后板材壓彎段成形質量較高。

(4)三輥非對滾彎試驗結果表明,有限元模型的成形半徑誤差為6.8%,有著較高的精度和可靠性。

[1]刑偉榮. 卷板機的現狀與發展[J].鍛壓裝備與制造技術, 2010, 45(2): 10-16.

Xing Rongwei. Status and Development of the Bending Machine for Plate[J]. China Metal Forming Equipment & Manufacturing Technology, 2010, 45(2): 10-16.

[2]楊建國, 方洪淵, 萬鑫, 等. 大尺寸板材滾彎過程應力場特點分析[J].塑性工程學報, 2008, 16(5):684-687.

Yang Jianguo, Fang Hongyuan, Wan Xin, et al. Character of Stress in Large Dimensional PlateduringThree-roll Bending[J]. Journal of Plasticity Engineering, 2008, 16(5):684-687.

[3]仉志強, 宋建麗, 付建華, 等. 板材彎卷成形工藝的研究現狀[J].塑性工程學報, 2014, 21(1): 1-6.

Zhang Zhiqiang, Song Jianli, Fu Jianhua, et al. Research Status of Plate Bending and Roll-bending Process[J]. Journal of Plasticity Engineering, 2014, 21(1): 1-6.

[4]Tran Q H, Champliaud H, Feng Zhengkun. Ana-lysis of the Asymmetrical Roll Bending Process through Dynamic FE Simulations and Experimental Study[J].The International Journal of Advanced Manufacturing Technology, 2014, 75(5/8):1233-1244.

[5]Feng Zhengkun, Champliaud H. Modeling and Simulation of Asymmetrical Three-roll Bending Process[J]. Simulation Modeling Practice and Theory, 2011, 19(9): 1913-1917.

[6]Zhang Zhiqiang, Song Jianli, Fu Jianhua, et al. A Refined Model of Three-roller Elastoplastic Asymmetrical Pre-bending of Plate[J]. Journal of Iron and Steel Research, International, 2014, 21(3):328-334.

[7]全建輝, 崔鳳奎, 楊建璽, 等.基于ANSYS/LS-DYNA的花鍵冷滾軋成形數值模擬[J].中國機械工程, 2008, 19(4):419-421.

Quan Jianhui, Cui Fengkui, Yang Jianxi, et al. Numerical Simulation of Involute Spline Shaft’s Cold-rolling Forming Based on ANSYS/LS-DYNA[J]. China Mechanical Engineering, 2008, 19(4):419-421.

[8]李森, 陳富林, 李斌, 等. 四輥卷板機彎卷過程中側輥位移計算[J].鍛壓技術, 2011, 36(6):76-79.

Li Sen, Chen Fulin, Li Bin, et al. Calculation of Side Roller Displacement for Four-roll Bending Plate Process[J]. Forging & Stamping Technology, 2011, 36(6):76-79.

[9]龔靖平, 黎向鋒, 左敦穩, 等.小直徑開縫襯套雙軸柔性滾彎成形的三維有限元分析[J].中國機械工程, 2015, 26(8):1117-1123.

Gong Jingping, Li Xiangfeng, Zuo Dunwen, et al.Three-dimensional FEA of Manufacturing Process of Small-diameter Split Sleeve by Two-axle Bending[J]. China Mechanical Engineering, 2015, 26(8):1117-1123.

[10]俞漢清, 陳金德. 金屬塑性成形原理[M].北京:機械工業出版社, 1999.

[11]仉志強, 宋建麗, 付建華,等.中厚板三輥彈塑性壓彎建模與試驗[J].塑性工程學報, 2013, 20(5):102-106.

Zhang Zhiqiang, Song Jianli, Fu Jianhua, et al. Modelling and Testing for Three-roller Elastic-plastic Bending of Moderate-thick Plate[J]. Journal of Plasticity Engineering, 2013, 20(5):102-106.

[12]劉章光, 高海濤, 劉繼偉, 等.TA15鈦合金板材冷折彎成形的有限元模擬及實驗研究[J].熱加工工藝, 2014, 43(13):129-133.

Liu Zhangguang, Gao Haitao, Liu Jiwei, et al. Finite Element Simulation and Experimental Study on Cold Bending Forming of TA15 Titanium Alloy Sheet[J].Hot Working Technology, 2014, 43(13):129-133.

(編輯陳勇)

Numerical Simulation and Experimental Verification of Asymmetrical Three-roll Bending Process

Wang Yan1,2Hu Jiefei1Xu Guanghui2Li Shusheng3Xiong Wei1

1.University of Shanghai for Science and Technology, Shanghai, 200093 2.Forging Machine Tool Plant of Zhenjiang, Zhenjiang, Jiangsu,212000 3.Taizhou Polytechnic College, Taizhou, Jiangsu,225300

The stress field of deformation zone, the stain field of upper plate surface and the changes of roll bending forces were analyzed during the asymmetrical and symmetrical roll bending process by the finite element(FE) simulation method. The simulation results show that the bending force during the asymmetrical roll bending process is larger than that during symmetrical roll bending under the same feed distance of the side roll. Compared to the asymmetrical roll bending process, the longitudinal and radial stress in deformation zone in the symmetrical roll bending process are smaller. The forming accuracy of the asymmetrical roll bending process of the pre-bending section is higher than that of symmetrical roll bending process.The experimental verification of a typical plates for the asymmetrical roll bending process was conducted, the result show that the forming error of FEA model is 6.8%.The FEA model has a higher forming accuracy.

asymmetrical three-roll bending; stress field; bending force; forming accuracy

2015-10-26

2015江蘇省重點研發計劃資助項目(BE2015139)

TG306

10.3969/j.issn.1004-132X.2016.15.018

王艷,女,1969年生。上海理工大學機械工程學院教授、博士。主要研究方向為磨削加工與特種加工。胡捷飛,男,1990年生。上海理工大學機械工程學院碩士研究生。許光輝,男,1975年生。鎮江市鍛壓機床廠工程師。李曙生,男,1972年生。泰州職業技術學院機電技術學院教授、博士。熊巍,男,1982年生。上海理工大學材料科學與工程學院講師。

猜你喜歡
有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 久无码久无码av无码| 日本道综合一本久久久88| 色综合天天操| 久久精品丝袜| 亚洲男人在线天堂| 欧美国产在线看| 欧美精品xx| 99草精品视频| 亚洲欧美一级一级a| 国产成人一区免费观看| 国产在线自乱拍播放| 日本91视频| 激情国产精品一区| 国产中文一区二区苍井空| 亚洲AV永久无码精品古装片| 日本精品αv中文字幕| 2020最新国产精品视频| 久久不卡国产精品无码| 欧美日韩福利| 农村乱人伦一区二区| 波多野结衣无码AV在线| 色噜噜狠狠狠综合曰曰曰| 伊人91在线| 一级毛片免费不卡在线| 热思思久久免费视频| 国产一区二区免费播放| 99精品免费欧美成人小视频| 久久久精品无码一区二区三区| 日韩a级片视频| 欧美一区二区自偷自拍视频| 妇女自拍偷自拍亚洲精品| 久久久四虎成人永久免费网站| 91色国产在线| jizz在线免费播放| 国产精品亚洲一区二区三区z| 亚洲一区波多野结衣二区三区| 91人人妻人人做人人爽男同| 亚洲精品爱草草视频在线| 91精品国产91久无码网站| 成人韩免费网站| 怡红院美国分院一区二区| 国产亚洲视频中文字幕视频| 性喷潮久久久久久久久| 精品久久久久成人码免费动漫| 国产成人久久777777| 影音先锋丝袜制服| 国产香蕉在线视频| 亚洲欧洲AV一区二区三区| 五月激情婷婷综合| 中文字幕av一区二区三区欲色| 日本不卡在线播放| 亚洲中文无码h在线观看| 999国产精品永久免费视频精品久久| 97在线视频免费观看| 国产农村妇女精品一二区| 欧美成一级| 伊人成色综合网| 国产成人无码综合亚洲日韩不卡| 国产亚洲精品自在线| 激情综合五月网| 性视频久久| 免费看黄片一区二区三区| 国产日韩久久久久无码精品| 国产精品毛片一区视频播| 欧美第二区| 色婷婷亚洲十月十月色天| 国产精品永久免费嫩草研究院| 亚洲v日韩v欧美在线观看| 麻豆精品国产自产在线| 亚洲综合精品香蕉久久网| 中文字幕 91| 国产手机在线小视频免费观看| 久久伊人色| 亚洲视频影院| 日本91在线| 国产女人18毛片水真多1| 91久久精品日日躁夜夜躁欧美| 亚洲国产成人无码AV在线影院L| 97精品伊人久久大香线蕉| 国产免费黄| 一区二区三区四区精品视频 | 国产黑人在线|