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

多方向耦合振動連軋機再生顫振建模及應用

2017-05-19 03:43:21劉曉潺郜志英曾令強
中南大學學報(自然科學版) 2017年3期
關鍵詞:振動模型

劉曉潺,臧 勇,郜志英,曾令強

多方向耦合振動連軋機再生顫振建模及應用

劉曉潺,臧 勇,郜志英,曾令強

(北京科技大學 機械工程學院,北京,100083)

為了研究連軋機的再生顫振機理,首先綜合考慮軋制過程中的軋輥彈性壓扁和軋件加工硬化效應的影響以及軋機工作輥的振動所造成的軋制輥縫在垂直和水平方向的變化,建立新的動態軋制過程模型,將該模型與多方向耦合振動的軋機結構模型結合,建立單機架軋機振動模型,而后對機架間的張力耦合和軋件厚波動在機架間的有時滯傳遞進行表征,建立新的多方向振動耦合的連軋機再生顫振系統模型。以所建立的模型為基礎,根據勞斯判據和積分值檢驗法分別對典型五機架冷連軋機每個單機架系統和多機架耦合的再生顫振系統進行穩定性分析。研究結果表明:所建立的再生顫振模型正確。

動態軋制過程模型;連軋機;再生顫振;時滯

軋機振動嚴重制約了軋制速度的提高,導致軋制產品質量下降,嚴重時甚至造成斷帶或者軋制設備被損壞,是目前國內外研究的熱點之一。軋機的振動模態多種多樣,振動機理復雜。軋機振動的形式主要包括軋機機架和輥系的垂直振動、傳動系統的扭轉振動和輥系的水平振動等[1?2]。軋機的自激振動機理有軋制過程和軋機結構之間的相互作用,由軋制張力波動所造成的負阻尼效應,軋機不同振動模態之間的耦合,以及相鄰機架之間的耦合作用所導致的再生顫振等[2]。其中軋機的再生顫振機理最復雜。軋機的再生顫振是指上一時段本機架軋機的振動再次引起該軋機的振動。軋機振動會在機架之間擴散,1個機架振動會帶動整個連軋機產生振動。在對寶鋼2030冷連軋機組的振動測試中,曾發現第四機架最先開始振動而后帶動其他機架一起振動的現象[1]。HU等[3?5]建立了連軋機再生顫振模型,對連軋機的再生顫振機理進行了解釋,表明機架間的張力耦合和軋件厚度波動在機架間的有時滯傳遞是導致連軋機再生顫振的重要原因。KIMURA等[6]建立了僅考慮垂直方向振動的五自由度的五機架冷連軋機再生顫振模型,其研究結果表明第五機架存在1個最優的摩擦因數使得連軋機系統穩定性達到最佳。NIROOMAND等[7]引入波的傳播理論對張力波動在機架間的傳遞進行了表征,通過仿真和實驗證明該方法更能準確地描述張力波動在機架間的傳遞。但是,HU等[8]在建立連軋機再生顫振模型時所用的動態軋制過程模型忽略了軋制過程中軋件的加工硬化和軋輥的彈性壓扁效應,而ZHAO等[4]建立的連軋機再生顫振模型只考慮了軋機在垂直方向的振動,沒有考慮軋機在水平和扭轉方向的振動。在冷軋過程中軋件的加工硬化和軋輥的彈性壓扁效應不可忽略,軋機不同方向振動之間也存在明顯的耦合作用[9]。為此,本文作者首先綜合考慮軋件的加工硬化和軋輥的彈性以及軋輥振動引起的軋制輥縫在垂直和水平方向的變化,建立了1個新的動態軋制過程模型,然后將該動態軋制過程模型和多方向耦合振動的軋機結構模型耦合,建立單機架軋機振動模型,在單機架軋機振動模型的基礎上,對機架間的耦合因素進行表征,建立新的多方向耦合振動的連軋機再生顫振模型,最后以典型冷連軋機參數為例,對模型進行驗證,以便為后續對連軋機振動機理的研究提供參考。

1 動態軋制過程模型的建立

1.1 軋制變形區幾何特征

當軋輥發生垂直和水平方向的振動時,軋制變形區動態幾何關系如圖1所示,其中:eσ為軋件入口張應力;dσ為軋件出口張應力;ve為軋件的入口速度;vr為軋輥的線速度;vd為軋件的出口速度;he為入口軋件厚度;hn為中性點處軋件厚度;hd為出口軋件厚度;hc為出口輥縫高度;xc為軋件水平出口位置,以軋件中心線為X軸,上下軋輥的中心連線為Y軸;φ為軋制接觸弧上任意一點與軋輥中心的連線和Y軸的夾角。將軋制接觸弧看作拋物線,沿接觸弧任意水平位置的軋件高度如下:

圖1 動態軋制輥縫Fig. 1 Dynamic roll gap

圖2 軋制變形區內的體積微元Fig. 2 Volume element in roll gap

將該體積微元截面視為梯形,其體積為

體積V對時間求偏導,忽略二階項,可得

將式(4)代入式(2),忽略二階項,可得

對式(1)兩邊求關于t的偏導,可得

對式(1)兩邊求關于x的偏導,可得

將式(6)和(7)代入式(5),以入口位置的軋件速度ve和入口位置軋件厚度he為初始條件,求解得到沿軋制區任意水平位置x的水平軋制速度:

將式(1)代入式(8),即得修正后的秒流量相等方程:

此外,由于軋制壓力的作用,軋制過程中存在軋輥的彈性壓扁,引入Hitchcock的軋輥彈性壓扁模型,變形后的軋輥等效半徑可以表示為

式中:Rw為軋輥原始半徑;λ為軋輥材料泊松比;p為接觸區單位軋件寬度的軋制壓力;r為壓下量;Ew為軋輥材料的彈性模量。根據軋輥半徑和軋制壓力之間的關系,進行收斂運算,就可以得到動態軋制力下的軋輥壓扁等效半徑。

1.2 軋件的出入口位置和速度

在中性點xn處,軋輥線速度近似等于軋制水平速度,根據式(8)得

對上式進行整理,得到軋件的入口速度ve:

同樣,可以得到軋件出口速度vd:

根據式(1),得到軋件的入口位置xe: d

φ為軋件出口位置和軋輥中心的連線與Y軸的夾角,因為該角度較小,滿足以下關系:

對式(15)進行整理,得到vd并將其代入式(13),化簡可得:

因為xe- xc遠大于xd- xc,對式(16)進行化簡可得出口位置為

1.3 軋制力和軋制扭矩計算

在均勻變形假設條件下,在變形區取任意微元段,軋制入口區和出口區的應力關系如圖3所示,變形區軋制力平衡方程如式(18)所示,其中出口區用下側的正號,入口區用上側的負號。

忽略二階量dσx· d h,并據tanφ= dh/(2d x),可得:

在冷軋過程中軋件的加工硬化不可忽略,軋件變形抗力隨軋件應變的變化如下:

圖3 入口區和出口區的軋制壓力分布切片法分析Fig. 3 Slab analysis on a volume element of entry region and exit region

其中:σ0,A和n為與材料相關的系數,可以通過實驗測定,對于連軋機h0是指第1機架的入口軋件厚度。將塑性變形條件式σx= kf- p及庫侖摩擦模型τs=μ p代入式(19),可得

引入軋件的加工硬化后,增加了式(21)的求解難度。從入口區到出口區,h持續減小,kf持續增加。因此假設hkf是一個常數,忽略第3項。但若直接將第3項忽略,則無法考慮到加工硬化對軋制力的影響,因此,在忽略第3項后,方程兩端同時除以hkf,得到式(22)。根據文獻[10],忽略第3項后,軋件厚度越小,根據式(22)得到的單位軋制壓力分布與根據卡爾曼方程得到的軋制壓力分布越接近。因此,對于冷軋薄板帶材時,該假設成立。

方程(22)的邊界條件為

可以得到入口區單位軋制壓力pe和出口區單位軋制壓力pd為:

式中:

在中性點處,pe= pd,由此得到軋制區中性點位置:

沿軋件和軋輥的接觸弧積分,得到垂直方向的軋制力Fv、軋制扭矩M和水平方向的軋制力Fh:

因為μtanφ? 1,對式(30)進行化簡得到:

1.4 軋制力能參數增量模型

軋制過程中的各種參數波動或者外界擾動會造成軋制力波動,軋制力的波動是造成軋機振動的直接原因。為了簡化計算,對軋制力能參數的變化量做線性化處理。將軋制力能參數對軋制工藝參數進行泰勒級數展開,并取其一階小量。軋件的入口速度ve、出口速度vd、垂直方向軋制力Fv、水平方向軋制力Fh以及軋制扭矩M隨輥縫波動量的一階泰勒展開式可表示為

綜上得到動態軋制過程模型如下:

2 連軋機再生顫振模型的建立

2.1 單機架軋機振動模型的建立

軋機振動是軋制過程和軋機結構相互作用的結果,選擇合適的軋機結構模型與動態軋制過程模型相互耦合才能建立合理的軋機振動系統模型。將軋機機架簡化為如圖4所示的彈簧質量塊模型。其中:m1為上工作輥的等效質量;k1為由工作輥偏移距產生的上支承輥系對上工作輥的水平等效剛度;c1為由工作輥偏移距產生的上支承輥系對上工作輥的水平阻尼;m2為集中于上工作輥中心的上部機架和上部輥系的等效質量;k2為機架上立柱和橫梁的等效剛度;c2為機架立柱及上橫梁的阻尼;Jm為上工作輥的等效轉動慣量;k3為上工作輥及其傳動系統的等效扭轉剛度;c3為上主傳動系統的阻尼系數[9]。

單機架軋機的振動模型如下:

圖4 軋機結構模型Fig. 4 Mill stand model

式中:w為軋件寬度。忽略軋件寬展,在軋制過程中

軋件寬度不變。

2.2 多機架軋機再生顫振模型的建立

連軋機相鄰機架間的關系表征如圖5所示,相鄰

機架通過軋件張力和厚度波動在機架間有時滯傳遞耦

合。機架間的張力波動與速度波動滿足以下關系[4]:

圖5 連軋機相鄰機架間的關系表征Fig. 5 Relationship between adjacent stands

其中:E為軋件材料的彈性模量;i為機架號;iiL,1-為第1-i機架與第i機架之間的距離;1,+iiL為第i機架與第i+1機架之間的距離。

在連軋過程中,上游機架振動所產生的軋件厚度波動在機架之間傳播,經歷過一個時間滯后傳入下游機架,即第i+1機架軋件入口厚度和第i機架出口厚度之間滿足以下關系:

其中:,1iiτ+為第i機架和第i+1機架之間的時滯時間,

ve,i+1為第i+1機架的軋件入口速度。

以單機架軋機振動模型為基礎,根據相鄰機架間的關系,可以得到相鄰的第i機架和第i+1機架再生顫振模型如下:

根據式(38)可知:每增架1個機架,系統增加7個自由度,據此可以建立任意機架耦合在一起的多方向耦合振動連軋機再生顫振模型。

3 連軋機再生顫振系統模型的驗證

3.1 單機架軋機振動系統穩定性分析

以某鋼廠五機架冷連軋機為例,對所建立的軋制過程模型、單機架軋機振動模型以及連軋機再生顫振模型進行驗證。該五機架冷連軋機每個機架的軋制工藝參數如表1所示。每個機架結構相似,因此,每個機架的結構參數設置相同,m1=4 800 kg,m2=117 900 kg,Jm=138 kg,k1=2.9×108N/m,k2=1.8×1010N/m,k3=24 MN/m,c1=3.9×105,c2=2.4×107,c3=5.9×103,w=1 m, L=4.5 m。

表1 五機架冷連軋機軋制工藝參數Table1 Process parameters of the five stand tandem cold rolling m ills

根據前面所建立的軋制過程模型,對每一機架的軋制力進行計算,將其與實測軋制力進行對比,并進行誤差分析,結果如表2所示。從表2可見:計算軋制力和實測軋制力之間的誤差較小,證明了本文所建立的軋制過程模型的正確性。根據勞斯判據,分別對每一個單機架進行穩定性分析,求解每一個單機架的臨界速度。對比分析5個機架的臨界速度,由于第1和第2機架的壓下量較大,軋件變形抗力較小,其臨界速度明顯比后續3個機架的小;后續3個機架的臨界速度比較接近,但是第4和第機架的軋制工作速度比較大,因此,第4和第5機架穩定性較差,在實際生產中,也是第4和第5機架振動較劇烈。根據秒流量相等原則,當第5機架的軋制速度為29.36 m/s時,第4機架的軋制速度為21.58 m/s,這表明在第5機架達到其臨界速度時,第4機架還遠未達到其臨界速度,因此,第5機架穩定性最差。

表2 計算軋制力和實測軋制力的對比以及每個單機架的臨界速度Table2 Comparison of calculated rolling force with measured rolling force and critical velocity for each stand

給定系統初始量yc0=1×10?6,對vr5=29.36 m/s時第5機架系統的時域和頻域進行仿真分析,結果分別如圖6和圖7所示。從圖6和圖7可見:第5機架在其臨界速度處做周期震蕩,垂直系統的振動頻率為166 Hz,扭轉系統的振動頻率為166 Hz和17.7 Hz,水平系統的振動頻率為166 Hz和44.3 Hz,但是扭轉系統和水平系統的振動以166.1 Hz為主,17.7 Hz和44.3 Hz的振動比較微弱,這表明第5機架自激振動以3倍頻為主。

圖6 vr5=29.36 m/s時第5機架的時域響應Fig. 6 Time domain response for the 5th stand at vr5of 29.36 m/s

圖7 vr5=29.36 m/s時第5機架的頻域響應Fig. 7 Frequency domain response for the 5th stand at vr5of 29.36 m/s

3.2 多機架軋機再生顫振系統穩定性分析

第4和第5機架穩定性較差,因此,建立第4和第5架機架耦合的2機架連軋機再生顫振模型。該模型是一個線性時滯微分系統,其特征方程是超越方程,具有無窮多個特征根,適用于無時滯系統的穩定性判定準則失效,需要采用適用于線性時滯微分系統的穩定性判定準則對其進行穩定性分析[11?13]。采用積分值檢驗法對連軋機再生顫振模型的穩定性進行分析,求解得到2機架軋機再生顫振模型的臨界速度為20.59 m/s(這里的臨界速度是指第5機架的軋輥線速度)。可見機架之間的耦合作用大幅度減小了系統的臨界速度。

給定初始量yc04=1×10?6,對vr5=20.59 m/s時2機架軋機再生顫振模型的時域和頻域特性進行仿真分析,結果如圖8所示。從圖8可以看出:在考慮時滯因素影響后,軋機振動波形形成了一連串的“小葫蘆”波形。這種“小葫蘆”波形是相鄰機架之間的相互作用造成的,與文獻[1]中的時域測試波形相同;由于第5機架的穩定性比第4機架的穩定性差,第5機架工作輥振動比第4機架更加劇烈;當vr5=20.59 m/s時,第4和第5機架之間的時滯時間為0.284 5 s。從圖8可以看出:由于第4和第5機架之間的張力耦合,第4機架的振動會直接引起第5機架的振動,但是張力耦合所導致的第5機架振動比較小,經過0.284 5 s后,第4機架振動所導致的軋件出口厚度波動傳遞到第5機架,加劇第5機架的振動,使得第5機架出現第1個“小葫蘆”波形,第5機架的第1個小葫蘆波形恰好落后于第4機架的第1個小葫蘆波形1個時滯時間。從圖9可以看出:在考慮了時滯因素的影響后,系統的振動頻率更復雜,但是第4和第5機架的振動仍然集中在3倍頻振動頻率附近,但是由于軋制速度剛達到其臨界速度,扭轉和水平方向的振動并不劇烈,振動頻率尚未集中;隨著軋機速度逐漸增大,整個機架將以4倍頻振動頻率發散。

以上對軋機振動模型的仿真分析和現場測試結果的比較,驗證了本文所建立的動態軋制過程、單機架軋機振動模型和多機架軋機再生顫振模型的有效性,對連軋機相鄰機架之間的耦合關系進行了研究,為后續深入研究連軋機再生顫振機理及其穩定性提供了理論基礎。

圖8 vr5=20.59 m/s多機架再生顫振模型的時域響應Fig. 8 Time domain response for multi-stand model at vr5of 20.59 m/s

圖9 vr5=20.59 m/s多機架再生顫振模型的頻域響應Fig. 9 Frequency domain response for multi-stand model at vr5of 20.59 m/s

4 結論

1) 建立了綜合考慮軋制過程中軋件的加工硬化和軋輥彈性壓扁效應以及軋制輥縫在垂直和水平方向波動的動態軋制過程模型,使用該模型對某鋼廠五機架連軋機每個機架的軋制力進行了計算,將計算結果與測試結果進行誤差分析,證明了該模型的有效性。

2) 將建立的動態軋制過程模型與多方向耦合振動的軋機結構模型結合,建立了多方向耦合振動的單機架軋機振動模型。基于該模型,使用勞斯判據對某五機架連軋機的每個機架穩定性進行分析,并求解了其臨界軋制速度。分析結果表明,該五機架連軋機的第4和第5機架的較容易發生振動,第5機架振動最容易發散。

3) 機架間的耦合作用降低了系統的穩定性,減小了系統的臨界速度。考慮機架間耦合作用的連軋機再生顫振模型時域仿真結果與軋機測試結果更接近,計算得到的臨界速度也更接近軋機的實際工作狀態,證明了所建立的再生顫振模型的有效性,為后續對連軋機再生顫振機理的研究奠定了基礎。

[1] 鄒家祥, 徐樂江. 冷連軋機系統振動控制[M]. 北京: 冶金工業出版社, 1998: 1?87.

ZHOU Jiaxiang, XU Lejiang. Tandem mill vibration control[M]. Beijing: Metallurgical Industry Press, 1998: 1?87.

[2] YUN I S, W ILSON W R D, EHMANN K F. Review of chatter studies in cold rolling[J]. International Journal of Machine Tools and Manufacture, 1998, 38(12): 1499?1530.

[3] HU Peihua, ZHAO Huyue, EHMANN K F. Third-octave-mode chatter in rolling, Part 3: Stability of a multi-stand mill[J]. Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering M anufacture, 2006, 220(8): 1293?1303.

[4] ZHAO Huyue, EHMANN K F. Stability analysis of chatter in tandem rolling mills, Part 1: single- and multi-stand negative damping effect Journal of Manufacturing Science and Engineering[J]. Journal of Manufacturing Science and Engineering, 2013, 135(3): 031001.

[5] ZHAO Huyue, EHMANN K F. Stability analysis of chatter in tandem rolling mills—Part 2: the regenerative effect, Journal of Manufacturing Science and Engineering[J]. Journal of Manufacturing Science and Engineering, 2013, 135(3): 031002.

[6] KIMURA Y, SODANI Y, NISHIURA N, et al. Analysis of chaffer in tandem cold rolling m ills[J]. ISIJ International, 2003, 43(1): 77?84.

[7] NIROOMAND M R, FOROUZAN M R, SALIM I M. Theoretical and experimental analysis of chatter in tandem cold rolling mills based on wave propagation theory[J]. ISIJ International, 2015, 55(3): 637?646.

[8] HU PEIHUA, EHMANN K F. A dynamic model of the rolling process. Part I: homogeneous model[J]. International Journal of Machine Tools and Manufacture, 2000, 40(1): 1?19.

[9] 曾令強, 臧勇, 郜志英, 等. 軋機整體耦合建模問題研究[J]. 機械工程學報, 2015, 51(14): 46?53.

ZENG Lingqiang, ZANG Yong, GAO Zhiying, et al. Study on overall coup led modeling of the rolling m ill[J]. Journal of M echanical Engineering, 2015, 51(14): 46?53.

[10] ZHAO Huyue. Regenerative chatter in cold rolling[D]. Evanston: Northwestern University, Department of Mechanical Engineering, 2008: 37?60.

[11] 孫健, 陳杰. 時滯系統穩定性分析與應用[M]. 北京: 科學出版社, 2012: 1?24.

SUN Jian, CHEN Jie. Stability analysis and application of time-delay systems[M]. Beijing: Science Press, 2012: 1?24.

[12] 王在華, 胡海巖. 時滯動力系統的穩定性與分岔: 從理論走向應用[J]. 力學進展, 2013, 43(1): 3?20.

WANG Zaihua, HU Haiyan. stability and bifurcation of delayed dynam ic systems: from theory to application[J]. Advances in mechanics, 2013, 43(1): 3?20.

[13] 張冬梅, 俞立. 線性時滯系統穩定性分析綜述[J]. 控制與決策, 2008, 23(8): 841?849.

ZHANG Dongmei, YU Li. Survey on the stability analysis of linear time delay systems[J]. Control and Decision, 2008, 23(8): 841?849.

(編輯 趙俊)

M ultidirectional regenerative chatter model of tandem rolling m ills and its app lication

LIU Xiaochan, ZANG Yong, GAO Zhiying, ZENG Lingqiang

(School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China)

To explore the regenerative chatter mechanism in rolling, a new dynam ic rolling process model was established firstly. The new model considers the work roll flattening effect, the rolled piece working hardening effect and the roll gap variation generated in vertical and horizontal direction by the work roll vibration. A single stand vibration model was built by coupling the new dynam ic rolling process model with multi direction mill stand model. And then the interaction factors between adjacent rolling m ill stands, namely the inter-stand tension variation and strip gauge variation passed onto the next stand after a delay time, were represented in functions. Finally, a new multi direction multi stand regenerative chatter model was built. Based on the models, stabilities of each single stand vibration model and multi stand regenerative chatter model were analyzed using the Routh criterion and the integral criterion respectively. The results prove that the proposed model is correct.

dynam ic rolling process model; tandem rolling mills; regenerative chatter; time delay

TG331;TH113

A

1672?7207(2017)03?0635?09

10.11817/j.issn.1672-7207.2017.03.011

2016?03?25;

2016?06?28

國家自然科學基金資助項目(51175035);教育部博士點基金資助項目(20100006110024);北京高校青年英才計劃項目(YETP0367);中央高校基本科研業務費專項資助項目(FRF-BR-14-006A);國家高技術研究發展計劃(863計劃)項目(2013AA031302) (Project(51175035) supported by the National Natural Science Foundation of China; Project(20100006110024) supported by Ph.D Programs Foundation of M inistry of Education of China; Project(YETP0367) supported by Beijing Higher Education Young Elite Teacher Project; Project(FRF-BR-14-006A) supported by the Fundamental Research Funds for the Central Universities; Project(2013AA031302) supported by the National High Technology Research and Development Program (863 Program) of China)

臧勇,教授,博士生導師,從事機械動力學和塑性加工研究;E-mail: yzang@ustb.edu

猜你喜歡
振動模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
重要模型『一線三等角』
This “Singing Highway”plays music
重尾非線性自回歸模型自加權M-估計的漸近分布
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 欧美日一级片| 97国产在线视频| 欧美国产综合色视频| a毛片免费在线观看| 欧美精品1区| 天天综合天天综合| 国产性爱网站| 精品夜恋影院亚洲欧洲| 欧美日本在线一区二区三区| 全部免费毛片免费播放| 中文字幕人成乱码熟女免费| 成人福利在线免费观看| 天天躁夜夜躁狠狠躁躁88| 丁香婷婷在线视频| 精品视频福利| 久久久久青草大香线综合精品| 色综合久久无码网| h视频在线播放| 亚洲国产日韩欧美在线| 国产精品污视频| 成人免费网站在线观看| 亚洲aⅴ天堂| 性做久久久久久久免费看| 蝴蝶伊人久久中文娱乐网| 自拍偷拍一区| 玩两个丰满老熟女久久网| 国产一级裸网站| 久久亚洲国产视频| 粉嫩国产白浆在线观看| 视频一本大道香蕉久在线播放 | 日韩黄色大片免费看| 国产精品黄色片| 国产午夜不卡| 青青草91视频| 亚洲乱伦视频| 国产尤物jk自慰制服喷水| 尤物精品视频一区二区三区| 一区二区日韩国产精久久| 国产精品欧美日本韩免费一区二区三区不卡| 91国内外精品自在线播放| 国产免费精彩视频| 亚洲日韩图片专区第1页| 漂亮人妻被中出中文字幕久久| 日韩精品亚洲一区中文字幕| 亚洲欧美h| 国产打屁股免费区网站| 另类专区亚洲| 国产区成人精品视频| 久久免费看片| 欧洲欧美人成免费全部视频 | 日韩精品一区二区三区视频免费看| 亚洲视频四区| 国产人前露出系列视频| 国产精品入口麻豆| 亚洲中文字幕国产av| 国产麻豆91网在线看| 日韩精品一区二区三区大桥未久| 成人一区在线| 中文字幕亚洲精品2页| 国产乱人免费视频| 日韩黄色大片免费看| 无码又爽又刺激的高潮视频| 欧美亚洲一区二区三区导航| 手机看片1024久久精品你懂的| 亚洲精品天堂自在久久77| 天堂av综合网| 欧美一道本| 无码aⅴ精品一区二区三区| 日本道中文字幕久久一区| 国产成人亚洲无吗淙合青草| 国产成人精品18| 亚洲第一网站男人都懂| 国产精品蜜臀| 日本不卡视频在线| 国产精品福利在线观看无码卡| 国产人人干| 国产成人精品视频一区视频二区| 国产免费羞羞视频| 午夜啪啪福利| 国产真实自在自线免费精品| 国产高清在线观看91精品| 午夜啪啪福利|