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

基于生死單元法的雙輥鑄軋過程熱-力耦合數值模擬

2015-10-29 02:30:35黃華貴劉文文杜鳳山
中國機械工程 2015年11期
關鍵詞:模型

黃華貴 劉文文 王 巍 杜鳳山

1.國家冷軋板帶裝備及工藝工程技術研究中心,秦皇島,0660042.燕山大學,秦皇島,066004

基于生死單元法的雙輥鑄軋過程熱-力耦合數值模擬

黃華貴1,2劉文文1,2王巍2杜鳳山1,2

1.國家冷軋板帶裝備及工藝工程技術研究中心,秦皇島,0660042.燕山大學,秦皇島,066004

以二輥φ160 mm×150 mm立式鋁帶實驗鑄軋機為對象,基于MSC.Marc商用有限元軟件及其二次開發接口,引入純鋁固-液兩相材料本構模型和界面壓力熱阻數學模型,建立了雙輥鑄軋過程熱-力耦合非線性有限元模型。采用生死單元法模擬鋁液的連續澆入,解決了輥套和鑄軋區鋁帶材間的連續耦合傳熱問題。通過數值模擬,給出了鑄軋速度、澆鑄溫度、熔池高度等因素對KISS點位置和輥面溫度時間歷程的影響規律,并對典型工況溫度模擬結果進行了實驗驗證。

雙輥鑄軋;溫度場;有限元法;生死單元法;數值模擬

0 引言

雙輥薄帶鑄軋工藝是一種集金屬快速凝固和熱軋于一體的復合工藝,具有短流程和低能耗等優點,屬于冶金及材料研究領域的前沿技術[1]。近年來,國內外學者圍繞鑄軋區金屬熔體及帶坯與鑄軋輥間的換熱行為、鑄軋帶材組織性能控制以及鑄軋工藝,開展了大量的理論與實驗研究工作[2-7]。經過近十年的快速發展和應用實踐,該技術已逐漸被應用于鋁合金、低碳鋼以及不銹鋼和硅鋼等難變形金屬帶材的工業生產[8]。作為雙輥鑄軋工藝的重要基礎,鑄軋區內固-液相組成及其對軋制力和耦合傳熱行為的影響一直是國內外學者研究的熱點,也是鑄軋輥冷卻結構的設計和鑄軋工藝優化的重要依據。

數值模擬技術作為鑄軋過程仿真的重要手段[9],可實現鑄軋過程熱-力耦合建模和輥體熱平衡的動態仿真分析,其核心技術包括:金屬熔體固-液混合相(或半固態)變形抗力模型的建立、鑄軋區帶坯與鑄軋輥表面換熱系數的計算、熔體連續澆入模型的實現等。在實際鑄軋生產過程中,鑄軋區內金屬形態主要由液相、固相和固-液兩相混合體組成,且需要軋制多圈后輥體溫度場才會達到平衡狀態。一般認為,金屬熔體具有牛頓流體的流變性能,固相金屬的塑性流變性能目前已有較為成熟的理論模型,而固-液兩相混合體鑄軋變形過程的流變應力曲線則需通過半固態金屬的熱壓縮實驗測試獲得[10-11]。湛利華等[5,12]在Gleeble上開展了鋁合金連續鑄軋過程流變行為物理模擬研究,并將其成功應用于鑄軋工藝熱-力耦合分析。鑄軋區金屬熔體、固相帶坯與輥面之間換熱行為方面,國內外學者主要采用接觸熱阻測試實驗,得到不同溫度、接觸壓力和表面粗糙度等條件下的接觸換熱系數,并應用于鑄軋輥溫度場模擬[13-15]。為獲得鑄軋輥輥體的平衡態溫度場,對鑄軋輥與帶坯進行耦合傳熱有限元建模時,持續多圈的鑄軋過程模擬需要大量的網格作為支持,這給模型求解帶來較大困難。因此,解決金屬熔體連續澆入建模問題,對實現鑄軋過程熱-力耦合和輥體熱平衡的動態仿真分析具有重要意義。

本文以燕山大學二輥φ160 mm×150 mm鑄軋機為對象,以鑄軋區流變本構模型和接觸換熱數學模型為基礎,通過MSC.Marc二次開發子程序接口,采用生死單元法建立了純鋁雙輥鑄軋過程的熱-力耦合有限元模擬模型,解決了金屬熔體的連續澆鑄建模問題。目前,生死單元法在帶鋼卷取[16]、焊接[17]和工程施工[18]等工藝研究方面已有較多應用。

1 雙輥鑄軋熱-力耦合有限元建模

1.1模型簡化

二輥φ160 mm×150 mm立式實驗鑄軋機如圖1所示,設計最大軋制力50 kN(5 t),最大鑄軋力矩1200 N·m,主電機功率為3 kW,具有振動鑄軋功能[19],以改善鑄軋帶材的內部質量。為便于有限元建模,進行如下假設:①鑄軋輥為剛性輥,鋁液入口速度和溫度恒定,帶坯與輥面不發生相對滑動;②鑄軋區內熔體、鑄坯與鑄軋輥沿寬度方向傳熱均勻,軋制力軸向分布均勻,鑄軋成形過程簡化為平面應變問題;③輥套與輥內冷卻水對流換熱系數為定值;④鑄軋輥與熔體間的換熱系數取常數,輥面與帶坯間的換熱系數與接觸壓力、溫度有關。

1.壓下裝置 2.平衡彈簧 3.鑄軋輥 4.鑄軋帶材5.澆鑄系統 6.鑄軋輥 7.振動機構圖1 二輥鑄軋機結構簡圖

為減小計算規模,根據傳熱和變形的對稱性,取二分之一結構進行建模,在MSC.Marc軟件平臺上建立純鋁雙輥鑄軋過程熱-力耦合仿真模型,如圖2所示。以鑄軋厚度為3 mm的鋁帶產品為例進行建模分析,輥套厚30 mm,取鑄軋區出口中心為坐標原點o,高度方向為y軸,輥面A點(圖2)為驗證實驗溫度檢測點。圖2中,ya為溶池液面高度,yk為單元“殺死”的臨界坐標,h2為輥套和冷卻水之間的對流換熱系數。

圖2 鑄軋過程熱-力耦合模型

1.2傳熱邊界條件

(1)鑄軋帶坯與輥套之間接觸換熱邊界表達為

qf=h1(tZP-tGTW)

(1)

其中,qf為軋輥與帶坯間的熱流密度;tZP、tGTW分別為帶坯表面溫度和輥套表面溫度;h1為帶坯與輥套之間的接觸換熱系數,與界面接觸壓力有關,文獻[14]通過界面熱阻測試實驗,給出了其回歸公式:

(2)

(2)輥套與內側冷卻水間的對流換熱邊界表達為

qw=h2(tGTN-tS)

(3)

式中,qw為輥套與冷卻水之間的熱流密度;tGTN和tS為輥套內表面溫度和冷卻水溫度。

模型中輥體、冷卻水初始溫度和環境溫度均為30 ℃。h2取14kW/(m2·K)。

(3)輥套與空氣間的對流換熱和輻射換熱等價為綜合換熱系數,取25W/(m2·K)。

(4)鑄軋帶坯中心對稱面取絕熱邊界條件。

1.3材料本構模型

工業純鋁的凝固溫度區間為614~659 ℃,凝固潛熱為3.9414J/g,其熱物性參數如表1所示[20]。

表1 純鋁熱物性參數

考慮固-液兩相混合體和固相純鋁的熱變形流變特性差異,引用文獻[12]中純鋁的材料本構模型:

(4)

材料本構模型通過MSC.Marc二次開發接口子程序urpflo.f嵌入有限元模型[21]。鑄軋輥材料選用42CrMo,相關熱物性參數直接從MSC.Marc材料數據庫中讀取。

1.4生死單元法及連續澆鑄建模

由于鑄軋輥表面與鑄坯間的傳熱具有間斷周期性特點,生產實踐表明,鑄軋輥達到熱平衡至少需要工作10圈以上,鑄軋帶材數十米。若采用傳統帶-輥耦合傳熱建模方法,則入口處需設置數量巨大的“金屬熔體”網格,給模型求解帶來困難。

如圖2所示,本文利用MSC.Marc中的生死單元法(二次開發接口子程序uactive.f),在鑄軋入口澆鑄區上方預先劃分足夠數量的網格,并將節點坐標y大于ya(即熔池液面高度)的單元全部“殺死”。隨著鑄軋輥的轉動,鑄軋區內金屬網格會將“死單元”帶入鑄軋區,當其節點坐標y小于ya時,通過子程序將新進入鑄軋區的“死單元”重新“激活”。同樣,為避免鑄軋出口鋁帶長度過大而導致模型網格數量過多,在出口方向以yk作為單元“殺死”的臨界坐標,單元節點坐標y小于yk時將被“殺死”。因此,整個模擬過程中僅節點坐標在yk和ya之間的單元保持“激活”狀態,參與模型求解的單元數量大大減少,解決了連續澆鑄模擬和快速仿真問題。

2 模擬結果分析

2.1鋁帶單元“激活”與“殺死”處理

圖3是鑄軋速度為2.4m/min時,完成第一圈軋制后的輥面溫度變化云圖。從圖3中可以看出,金屬熔體入口液面始終保持不變,當鑄軋出口帶坯延伸至yk點后,單元即被“殺死”。在有限元法中,對于被“殺死”的單元只是將其單元剛度矩陣乘上一個較小的參數(如1.0×10-6),并非將其從實際模型中刪除。由于被“殺死”單元的單元載荷、質量和熱邊界條件等其他同類參數均為0,故可大大提高模型求解效率。鋁液入口單元溫度初值則可在模型中設置,不受單元生死操作的影響。

(a)時間τ=0(b)時間τ=4 s

(c)時間τ=12 s(d)時間τ=18 s圖3 鑄軋過程輥面溫度分布變化

模擬獲得了鑄軋區穩態溫度場(圖4a)和帶坯等效應力分布(圖4b),從圖4中可以看出鑄軋區內液相區、固-液兩相區、固相區和KISS點高度,以及鋁帶等效應力分布。結合鑄軋輥與鋁帶間接觸壓力模擬結果(圖5)和式(2),可計算出鑄軋區內鋁帶和鑄軋輥間的接觸換熱系數分布曲線(圖5),接觸壓力從入口到出口先增大后減小,接觸換熱系數與接觸壓力的變化趨勢相一致。

(a)溫度場(b)等效應力圖4 穩態鑄軋區(τ=150 s)

圖5 鋁帶與鑄軋輥間接觸壓力和接觸換熱系數

2.2鑄軋速度對KISS點高度和出口溫度的影響

取熔池液面高度ya=30 mm,鑄軋速度v分別為1.6 m/min、2.4 m/min、3.2 m/min,澆鑄溫度t0=700 ℃,KISS點高度和鋁帶出口溫度變化規律如圖6所示。從圖6中可以看出,當v=3.2 m/min時,鋁帶出口溫度處于固液兩相區,會發生漏鋼事故;當v=1.6 m/min時,出口溫度較低,KISS點高度過高,易造成卡殼事故。模擬結果表明,當v=2.4 m/min時,KISS點高度由開始的12.74 mm隨時間逐漸降低并穩定于8 mm(圖6b),鑄軋過程穩定。

(a)出口溫度

(b)KISS點高度(v=2.4 m/min)圖6 鑄軋速度對KISS點高度和鋁帶出口溫度的影響

2.3澆鑄溫度對KISS點高度和出口溫度的影響

(a)出口溫度

(b)KISS點高度圖7 澆鑄溫度對KISS點高度和出口溫度的影響

取熔池液面高度ya=30 mm,鑄軋速度v=2.4 m/min,澆鑄溫度t0分別為680 ℃、700 ℃、720 ℃,KISS點高度和鋁帶出口溫度的變化規律如圖7所示。從圖7中可以看出,三種工況下KISS點高度隨鑄軋圈數的變化趨勢基本相似,澆鑄溫度每降低20 ℃,鋁帶出口溫度約降低10 ℃(圖7a),相同時刻KISS點高度升高約1 mm(圖7b)。

2.4熔池高度對KISS點高度和出口溫度的影響

取澆鑄溫度t0=700 ℃,鑄軋速度v=2.4 m/min,熔池高度ya分別為20 mm、30 mm和40 mm,KISS點高度和鋁帶出口溫度變化規律如圖8所示。當熔池高度為20 mm時,鋁帶出口溫度超過625 ℃,會出現漏鋼事故(圖8a);當熔池高度由30 mm升高至40 mm時,穩態KISS點高度由8 mm提高至19.3 mm,鋁帶出口溫度由570 ℃降低至498 ℃。

2.5輥套熱平衡分析及實驗驗證

根據上述模擬結果,確定適于本實驗鑄軋機的合理工藝參數為:熔池高度30 mm、澆鑄溫度700 ℃、鑄軋速度2.4 m/min,并以以上參數為基礎進行輥套熱平衡分析,模擬獲得的輥面及輥套徑向溫度隨時間變化曲線如圖9所示,圖中,h′為測點距輥面的距離。

(a)出口溫度

(b)KISS點高度圖8 熔池高度對KISS點高度和鋁帶出口溫度的影響

(a)輥面溫度

1.h′=5 mm 2.h′=10 mm 3.h′=15 mm4.h′=20 mm 5.h′=25 mm 6.h′=30 mm(b)輥套徑向溫度圖9 輥套熱平衡模擬結果

鑄軋過程輥套內各點溫度隨時間呈周期性變化,距輥面越近,溫度波動周期性越明顯(圖9b)。軋制9圈后輥套溫度變化達到穩定狀態,輥面最高溫度為298 ℃,最低溫度為125 ℃(圖9a),輥套內壁溫度為45 ℃(圖9b)。

為驗證模型精度,利用Raytek紅外測溫儀對相同工況下(圖10a)輥面溫度進行采集。由于工作過程軋輥處于轉動狀態,無法對輥面進行定點測溫,本文通過事先在輥面中心A點(圖2)對應的輥身端部貼簽做標識,分別連續鑄軋1圈、2圈、3圈后停機對A點進行溫度測試,結果表明,實測值與有限元模擬結果基本吻合(圖10b)。

(a)實驗現場

(b)輥面溫度對比圖10 溫度模擬結果驗證

3 結束語

本文較為系統地分析了雙輥鑄軋過程熱-力耦合有限元建模的核心數學模型,并采用生死單元方法,解決了金屬熔體連續澆鑄建模問題,為鑄軋輥熱平衡分析和工藝優化提供了新途徑。以二輥φ160 mm×150 mm立式實驗鑄軋機為對象,模擬分析了鑄軋速度、熔池高度、澆鑄溫度對鑄軋區KISS點高度和帶材出口溫度的影響規律,給出了合理的鑄軋工藝規程及相應工況下的輥套熱平衡演變過程,并通過鑄軋實驗對溫度模擬結果進行了驗證。

[1]趙紅陽, 胡林, 李娜. 雙輥薄帶鑄軋技術的進展及熱點問題評述[J]. 鞍鋼技術, 2007(6): 1-5.

Zhao Hongyang, Hu Lin, Li Na. Development of Two-high Thin Strip Casting and Rolling Technology and Key Point Review[J].Angang Technology,2007(6):1-5.

[2]Zhao Hu, Li Peijie, He Liangju. Coupled Analysis of Temperature and Flow During Twin-roll Casting of Magnesium Alloy Strip[J]. Journal of Materials Processing Technology, 2011, 211(6): 1197-1202.

[3]Wang Bo, Zhang Jieyu, Li Xiangmei, et al. Simulation of Solidification Microstructure in Twin-roll Casting Strip[J]. Computational Materials Science, 2010, 49(1):S135-S139.

[4]陳守東, 陳敬超. 雙輥薄帶連鑄凝固組織模擬微觀模型的驗證[J]. 材料熱處理學報, 2012, 33(7): 154-158.

Chen Shoudong, Chen Jingchao. Validation of Micro Model for Solidification Process of Twin-roll Continuous Casting Thin Strip[J]. Transactions of Materials and Heat Treatment, 2012,33(7):154-158.

[5]湛利華, 李曉謙, 唐朝陽, 等. 鋁帶坯連續鑄軋過程熱力耦合有限元分析[J]. 中國機械工程, 2005, 16(11): 979-984, 996.

Zhan Lihua, Li Xiaoqian, Tang Zhaoyang,et al. Thermo-mechanical Coupling Analysis of Aluminum Continuous Roll Casting Processes[J]. China Mechanical Engineering,2005, 16(11): 979-984, 996.

[6]Parka C M, Choia J T.Thermal Crown Analysis of the Roll in the Strip Casting Process[J]. Journal of Materials Processing Technology, 2009, 209(8):3714-3723.

[7]朱志華, 肖文鋒, 李曉謙. 雙輥連續鑄軋過程中軋制界面接觸壓力切塊法建模與數值模擬[J]. 中國機械工程, 2005, 16(13): 1091-1095.

Zhu Zhihua, Xiao Wenfeng, Li Xiaoqian. Modeling of Rolling Pressure Using Slab-method and the Numerical Simulation during Twin-roll Continuous Roll Casting Process[J]. China Mechanical Engineering,2005, 16(13): 1091-1095.

[8]Jerry W S, Robert J, Comstock J R. Method of Continuous Casting Non-oriented Electrical Steel Strip:US, 2006015z142A1[P].2006-07-13.

[9]許志強, 孟哲儒, 杜鳳山, 等. 雙輥薄帶鑄軋數值模擬研究現狀及展望[J]. 燕山大學學報, 2014,39(2):95-101.

Xu Zhiqiang, Meng Zheru, Du Fengshan, et al. Current Situation and Prospect of Twin-roll Strip Casting Process Numerical Simulation[J]. Journal of Yanshan University, 2014,39(2):95-101.

[10]謝金樂, 劉允中, 吳匯江, 等. 半固態7050鋁合金熱壓縮變形行為[J]. 特種鑄造及有色合金, 2011, 31(9): 816-819.

Xie Jinle, Liu Yunzhong, Wu Huijiang, et al. Hot Compression Behavior of Semi-solid 7050 Aluminum Alloy[J]. Special Casting & Nonferrous Alloys, 2011, 31(9): 816-819.

[11]湛利華, 鐘掘, 李曉謙, 等. 連續鑄軋流變行為的物理模擬及其應力-應變關系的演變[J]. 中國有色金屬學報,2004,14(12):1995-2002.

Zhan Lihua, Zhong Jue, Li Xiaoqian, et al. Physical Simulation of Rheological Behavior and Stress-strain Relation Evolvement in Continuous Roll Casting Process[J]. The Chinese Journal of Nonferrous Metals, 2004,14(12): 1995-2002.

[12]湛利華. 鋁合金連續鑄軋過程流變行為研究及熱-力耦合分析[D]. 長沙:中南大學, 2005.

[13]劉菊. 固體界面接觸熱阻及熱導率測量的實驗研究[D]. 武漢:華中科技大學, 2011.

[14]張云湘. 瞬態法接觸熱阻實驗研究及其在快速鑄軋工藝參數仿真中的應用[D]. 長沙:中南大學, 2004.

[15]Wang De, Zhou Cheng, Xu Guojin, et al. Heat Transfer Behavior of Top Side-pouring Twin-roll Casting[J]. Journal of Materials Processing Technology, 2014, 214(6): 1275-1284.

[16]王曉晨, 楊荃, 劉瑞軍, 等. 基于ANSYS有限元法的熱卷箱內中間坯溫度場分析[J]. 北京科技大學學報, 2013(4):454-458.

Wang Xiaochen, Yang Quan, Liu Ruijun, et al. Temperature Field Analysis of Intermediate Slabs in the Hot Coil Box Based on ANSYS Finite Element Method[J]. Journal of University of Science and Technology Beijing, 2013(4): 454-458.

[17]趙洪運, 舒鳳遠, 張洪濤, 等. 基于生死單元的激光熔覆溫度場數值模擬[J]. 焊接學報, 2010, 31(5): 81-84.

Zhao Hongyun, Shu Fengyuan, Zhang Hongtao, et al. Numerical Simulationon Temperature Field of Laser Cladding Based on Birth-death Element Method[J]. Transactions of the China Welding Institution,2010, 31(5): 81-84.

[18]鄭江, 葛鴻鵬, 王先鐵, 等. 局部位形約束生死單元法及其在施工力學分析中的應用[J]. 建筑結構學報,2012, 33(8): 101-108.

Zheng Jiang, Ge Hongpeng, Wang Xiantie, et al. A Method of Element Birth and Death of Local Configuration Con-straint and Its Application in Construction Mechanics[J]. Journal of Building Structures,2012, 33(8): 101-108.

[19]李堯. 雙輥薄帶振動鑄軋過程仿真模擬及實驗研究[D]. 秦皇島:燕山大學, 2011.

[20]王祝堂. 鋁合金及其加工手冊[M].3版. 長沙:中南大學出版社,2005.

[21]陳火紅, 尹偉奇, 薛小香. MSC.Marc二次開發指南[M]. 北京:科學出版社, 2004.

(編輯袁興玲)

Thermal-mechanical Coupled Modelling and Numerical Simulation for Twin-roller Casting Process with Technique of Deactivate and Reactivate Element

Huang Huagui1,2Liu Wenwen1,2Wang Wei2Du Fengshan1,2

1.National Engineering Research Center for Equipment and Technology of Cold Strip Rolling,Qinhuangdao,Hebei,066004 2.Yanshan University,Qinhuangdao,Hebei,066004

Aiming at the φ160 mm×150 mm vertical type experimental twin-rollers continuous caster, a thermal-mechanical coupled FEM model of twin-roller casting process was established with MSC.Marc and its’ secondary development interface based on the solid-liquid two phases constitutive model of pure aluminum and the thermal resistance model considering the contact pressure. The technique of deactivate and reactivate element method was used to simulate the continuous casting of aluminums liquid, and the problem of coupled heat transmission between the cast roll shell and aluminum strip in the cast rolling area was therefore resolved. From the simulation results, the influences of casting speed, casting temperature, the height of cast rolling area on the position of KISS point and the time history of temperature on the roll surface were presented. An twin-roller casting experiment was carried out, and the numerical simulation results are in good agreements with the experimental data.

twin-roller casting; temperature field; finite element method(FEM); deactivate and reactivate element method; numerical simulation

2014-06-30

國家自然科學基金資助項目(51474189);河北省自然科學基金資助項目(E2013203377)

TG233.6;TP391.9DOI:10.3969/j.issn.1004-132X.2015.11.013

黃華貴,男,1978年生。燕山大學機械工程學院教授、博士。主要研究方向為材料加工工藝及裝備。獲國家發明專利6項,發表論文30余篇。劉文文,男,1990年生。燕山大學機械工程學院碩士研究生。王巍,女,1978年生。燕山大學機械工程學院講師。杜鳳山,男,1960年生。燕山大學機械工程學院教授、博士研究生導師。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 无码日韩精品91超碰| www精品久久| 99久久精品免费看国产免费软件 | 国产高潮流白浆视频| 色亚洲激情综合精品无码视频| 在线观看亚洲精品福利片| 久久人人妻人人爽人人卡片av| av一区二区无码在线| yy6080理论大片一级久久| 免费在线色| 日韩欧美91| 人妻精品久久无码区| 国产真实自在自线免费精品| 精品国产电影久久九九| 国产91高清视频| 国产理论一区| 日韩精品无码不卡无码| 91精品免费高清在线| 99人妻碰碰碰久久久久禁片| 精品免费在线视频| 久久国产高清视频| 97超级碰碰碰碰精品| 国产在线拍偷自揄拍精品| 妇女自拍偷自拍亚洲精品| 免费一极毛片| a毛片在线| 黄色成年视频| 九九久久精品免费观看| 在线免费亚洲无码视频| 亚洲二区视频| 青草91视频免费观看| 99这里只有精品免费视频| 国产精品大白天新婚身材| 国产精品欧美亚洲韩国日本不卡| 婷婷综合色| 丝袜亚洲综合| 97视频精品全国在线观看| 免费三A级毛片视频| 青青草原偷拍视频| 欧美在线伊人| 国产在线观看91精品| 黑人巨大精品欧美一区二区区| 成人免费视频一区二区三区 | 免费毛片网站在线观看| 中文字幕va| 无码乱人伦一区二区亚洲一| 亚洲人成人伊人成综合网无码| 一区二区影院| 72种姿势欧美久久久大黄蕉| 黄色网在线免费观看| 亚洲欧洲综合| 999国内精品久久免费视频| 国产91成人| 国语少妇高潮| 久久精品国产精品一区二区| 日韩国产亚洲一区二区在线观看| 日本在线欧美在线| 日韩精品免费在线视频| 人人艹人人爽| 国产一线在线| 亚洲男人天堂久久| 久久伊伊香蕉综合精品| 国产福利在线免费| 99热这里只有精品2| 亚洲午夜福利精品无码不卡| 国产精品冒白浆免费视频| 一本大道香蕉高清久久| 精品国产美女福到在线不卡f| 国产主播在线一区| 国产超碰在线观看| 国产swag在线观看| 精品一区二区三区视频免费观看| 亚洲高清无码精品| 丰满人妻中出白浆| 国产精品久久久久久久久kt| 亚洲国产亚洲综合在线尤物| 91精品最新国内在线播放| 欧美激情,国产精品| 欧美成人在线免费| 亚洲另类国产欧美一区二区| 免费国产在线精品一区| 就去色综合|