在线观看国产区-在线观看国产欧美-在线观看国产免费高清不卡-在线观看国产久青草-久久国产精品久久久久久-久久国产精品久久久

美章網 資料文庫 材料本構模型參數的獲得方法范文

材料本構模型參數的獲得方法范文

本站小編為你精心準備了材料本構模型參數的獲得方法參考范文,愿這些范文能點燃您思維的火花,激發您的寫作靈感。歡迎深入閱讀并收藏。

材料本構模型參數的獲得方法

《中國有色金屬學報》2015年第十二期

摘要:

為了獲得材料在高應變率下的本構模型參數,提出一種基于激光沖擊強化實驗和有限元模擬相結合的材料本構模型參數獲得方法。該方法首先預估材料在高應變率下的本構模型參數,然后進行有限元模擬,以模擬結果與實驗結果對比的誤差落在一定范圍內為優化目標,對預估的本構模型參數進行修正,并最終得到其具體值。基于此方法,得出TiAl合金在高應變率下的Hugoniot彈性極限(HEL)為1.46GPa。

關鍵詞:

鈦鋁合金;本構模型;高應變率;激光沖擊強化;有限元模擬

激光沖擊強化是一種新型的表面強化技術,該技術是利用強激光束輻照涂覆在金屬表面的吸收層,使其氣化并等離子化形成高壓沖擊波并向材料內部傳播,引起材料高應變率(>1×106s1)動態響應,使材料表層產生塑性變形并形成高殘余壓應力,從而提高材料表面硬度、抗磨損、抗腐蝕和抗疲勞等力學性能[15]。目前,關于激光沖擊強化中沖擊波的產生與發展,沖擊波作用下材料的高應變率動態響應,材料顯微組織變化規律等許多問題的了解仍不是很清楚[6]。因此,許多學者采用實驗與有限元模擬相結合的方法來研究激光沖擊強化對材料表面性能的影響,但是由于受高應變率的影響,材料的本構模型參數明顯不同于準靜態情況下材料的本構模型參數(如對于大多數材料,動態屈服強度隨應變率的增大而增大),為有限元模擬帶來了困難,因而,確定材料在高應變率下的本構模型參數就顯得相當重要。為了研究材料在高應變率下的本構模型,研究者大多使用膨脹環、分離式Hopkinson桿和Taylor實驗等方法來獲得材料的力學性能數據,然后,通過擬合和外推得到其本構模型參數[79]。HOGGATT等[9]運用爆炸膨脹環技術,對許多工程材料的本構關系進行了測試,趙峰等[10]運用分離式Hopkinson桿對AZ91鑄造鎂合金在不同應變率下的動態本構關系進行了研究,REN等[11]運用Taylor桿沖擊實驗對Ti-6Al-4V合金在高應變率加載條件下的動態斷裂行為進行了研究。但是,這些實驗方法所獲得的應變率范圍大都在1×102~1×105s1范圍內,無法滿足激光沖擊強化模擬的要求。要想獲得大于1×105s1應變率下的材料動態力學數據,就必須采用沖擊、炸藥爆轟或輻射等方法[8],這些實驗方法不僅非常危險,而且實驗裝置復雜、成本較高,對材料本構模型參數的確定帶來了很大的困難。為了解決這個問題,本文作者在對TiAl合金進行激光沖擊強化實驗和有限元模擬的基礎上,提出了一種獲得材料在高應變率下本構模型參數的方法,并利用該方法,得到了TiAl合金在高應變率下的理想彈塑性本構模型參數。

1本構模型參數的獲得方法

在本構模型參數獲得方法中,首先預估模型參數,然后進行有限元模擬,以模擬結果與實驗結果對比的誤差落在一定范圍內為優化目標,對預估的材料本構模型參數進行修正,最終得到其具體值。在激光沖擊強化過程中,高壓沖擊波會引起材料的高應變率動態響應,因此,材料的本構模型需考慮高應變率的影響,高應變率下材料的本構模型主要有Johnson-cook模型,Zerilli-Armstrong本構模型,理想彈塑性模型(屈服強度為在高應變率下的屈服強度)等。Johnson-cook模型綜合考慮了應變、應變率和溫度等的影響,可以較好地反映材料的真實本構關系,但是,模型參數復雜,影響因素較多,多通過測量不同應變速率、不同溫度條件下的多組應力應變曲線,利用數值擬合的方法確定。DING等[12]使用理想彈塑性模型對35CD450HRC鋼進行了模擬研究,模擬結果與實驗結果誤差在7%左右。花銀群等[13]也使用該模型對TC4進行了研究,結果與實驗數據對比誤差在2%左右,由此可以看出,在利用理想彈塑性模型對激光沖擊強化進行模擬時,可以得到較好的結果,并且該模型參數簡單,不考慮應變硬化等因素,只需確定材料在高應變率時的Hugoniot彈性極限(HEL),當材料中應力波壓力小于HEL時,材料不發生塑性應變;當應力波壓力大于HEL時,將發生永久的塑性變形[6]。

2有限元模型構建

要想獲得正確的本構模型參數,建立合理的有限元模型是關鍵,在激光沖擊強化有限元模擬中,主要涉及到沖擊波壓力的計算與加載、材料本構模型的選擇、網格劃分與施加邊界條件等問題。由于激光沖擊強化過程中沖擊波壓力的加載時間非常短,在沖擊波傳播過程中,材料內部會發生各種應力波的反射和相互作用,因而,材料對沖擊的響應需要持續相當長的時間才能達到穩定,因此,為得到穩定的殘余應力場,有限元分析需要采用顯式動態分析與隱式靜態分析相結合的方法。顯式動態分析采用顯式動力有限元算法,分析瞬態沖擊下產生的沖擊波在材料中的傳播以及相互作用過程,得到材料動態響應。然后進行隱式靜態分析,釋放顯式動態分析后的材料內部彈性應變,計算材料內部的平衡狀態,得到穩定的殘余應力場[14]。

2.1沖擊波壓力波形計算與加載激光誘導的沖擊波在空間和時域的分布將對模擬結果產生極其重要的影響。在空間分布上,由于實驗中使用的激光脈沖能量在空間平面內近似均勻分布,因此,可認為激光誘導的沖擊波壓力在整個光斑內均勻分布。而在沖擊波的時域分布上,通常使用FABBRO等[15]提出的模型。該模型認為在水約束狀態下,激光產生的沖擊波壓力大致趨于一個六次多項式分布,沖擊波壓力脈沖的半峰值脈寬大致為激光脈寬的2~3倍,沖擊波的峰值壓力如式(2)所示。

2.2材料本構模型參數預估由于TiAl合金在高應變率下的HEL未知,因此,需要首先對其預估,在一維應變下,材料的HEL與動態屈服強度的關系為。因此,只需預估出材料的動態屈服強度,即可求出HEL。根據Campbell和Ferguson的研究,對于大多數材料,動態屈服強度隨應變率的增大而增大[8],圖3所示為普通低碳鋼屈服強度與應變率的關系[17]。從圖3可以看出:在不同應變率下,其屈服強度明顯不同,且隨著應變率的增大,屈服強度也越來越大。TiAl合金在準靜態情況下的屈服強度為600MPa[18],隨著應變率的增大,其屈服強度應該越來越大。因此,在本模擬中,可以首先設置TiAl合金在106s1應變率下的動態屈服強度為700MPa,比準靜態情況下高出17%。在模擬中使用的其他材料參數為ρ=3800kg/m3,E=165GPa,υ=0.32。

2.3有限元網格模型及邊界條件的確定在有限元模擬中沖擊應力波會在模型邊界處發生反射,反射波匯聚到中心會對結果造成一定影響,因此,為了防止應力波反射所造成的影響,需要使用半無限大三維實體模型,即在模型的側面和底面施加無反射邊界條件,讓應力波透射。模擬時,為了與實驗結果有更好的對比,所使用的模型尺寸與實驗所用的試件完全一致,即直徑為15mm、高為5mm的圓柱體,在圓柱體正中心沖擊,由于沖擊壓力與模型均關于中心對稱,為提高計算效率,僅建立1/4模型進行分析計算,并在對稱面上施加對稱邊界條件。最終的有限元模型如圖4所示。

3TiAl合金HEL的確定

3.1激光沖擊強化實驗及結果實驗時所用的試件照片如圖5所示,材料為TiAl合金,名義成分為Ti-45.5Al-2Cr-2Nb-0.15B。在由中科院沈陽自動化研究所搭建的型號為SIA-LSP-1的激光沖擊強化設備上進行,所用的高能脈沖激光器型號為Extra20,脈沖能量為9J,激光波長為1064nm,脈寬為20ns,采用直徑為3mm圓形脈沖在中心沖擊一次。在激光沖擊之前,先將試件在HY4050型豪克能應力消除設備(頻率為40kHz,輸出振幅為50μm)上處理5min,以消除表面殘余應力,再用酒精清洗并用氮氣吹干,然后,在試件表面粘貼一層100um厚的黑膠帶作為吸收層,并用厚度為2mm的流動水膜作為約束層。然后,在試件中心進行單次激光沖擊,沖擊之后使用Proto-LXRD型X射線應力測試儀測量激光沖擊區域表層與深度方向的殘余應力,測量方法采用傾斜固定ψ法,其衍射條件列于表1,深度方向殘余應力采用0.5mol/L的Na2SO4+H2SO4溶液進行逐層電解腐蝕測量,X射線輻照區域直徑2mm,即測量的是激光沖擊區域中心2mm內的平均應力[19]。由于模型與沖擊壓力的對稱性,x方向與y方向殘余壓力基本相等,分布規律也一致,因此,以下討論中均使用y方向的殘余應力進行表示。測量結果如圖6所示,數據誤差為±20MPa,從圖6可以看出,殘余壓應力最大值出現在次表層,最大值為350MPa,之后,隨著深度的增加,殘余壓應力逐漸降低,在0.65mm左右減小為0。

3.2誤差范圍的確定為了獲得正確的材料本構模型參數,必須給定合理的誤差范圍,若誤差范圍太大,則得到的結果不可靠;誤差太小,由于壓力模型及實驗測量結果等因素的誤差,可能根本得不到符合要求的結果。在激光沖擊強化有限元模擬中,產生誤差的來源很多,當高功率的激光脈沖輻照金屬表面后,涂覆在金屬表面的吸收層氣化并等離子化形成等離子體,該等離子體繼續吸收激光能量,并受到約束層的約束形成沖擊波,作用在材料表面,因此,該過程是一個非常復雜的物理過程,要對其進行真實模擬是非常困難的。在有限元模擬中,將該過程簡化為壓力脈沖直接作用在材料表面,從而該過程就引入了誤差,具體表現在壓力峰值大小、空間分布和時域分布等方面。另外,采用的材料本構模型是理想彈塑性模型,該模型本身就是對真實本構模型的一種簡化,因此,該模型存在原理上的誤差。在有限元求解時,網格劃分的大小,所采用的材料阻尼等也會對結果產生影響;最后,本實驗中采用逐層剝離法測量深度方向的殘余應力,由于在剝離部分材料后,材料的表面狀態發生了改變,會引起殘余應力的重新分布,因此,也會引入適量的誤差。由于各個誤差產生的復雜性,對每個誤差進行精確估計是非常困難的,且各個誤差之間的相互影響規律也很難確定,因此,若通過計算來確定總體誤差范圍是非常困難的。而從其他學者的有限元模擬結果與實驗結果對比分析中可以發現,當采用相似的實驗條件與模擬條件時,殘余壓應力最大值誤差與殘余壓應力影響深度誤差總是在一定的范圍內變動,具有一定規律,因此,可以這樣認為,該范圍即是綜合了各種誤差之后總的誤差范圍,可以對該誤差范圍加以分析,從而得到本例中所需確定的誤差范圍。圖7所示為從文獻[1213,2023]中得到的不同材料激光沖擊強化實驗結果與模擬結果對比的誤差,其中殘余壓應力最大值誤差(Errorofthemaximumcompressiveresidualstress,EMCRS),殘余壓應力影響深度誤差(Errorofaffecteddepthofcompressiveresidualstress,EADCRS)的負值表示模擬值比實驗值小,正值表示比實驗值大。從圖7可以看出,殘余壓應力最大值與實驗值之間的誤差主要集中在±10%之內,并且正負分布比較均勻,因此,本文作者將殘余壓應力最大值誤差范圍設定為±10%。殘余壓應力深度誤差變化范圍比較大,從為37.5%~14%,但誤差仍集中在20%~14%之間,因此,將殘余壓應力影響深度誤差范圍設定為20%~14%。

3.3模擬結果與實驗結果的對比在進行有限元模擬時,首先設定TiAl合金的動態屈服強度為700MPa進行模擬。由于實驗中X射線應力測試儀的輻照區域直徑為2mm,測量的是中心2mm內的平均應力,所以模擬結果也必須以2mm內的平均應力為基準進行對比分析。因此,在模擬沖擊的中心區域選取如圖8所示直徑為2mm的圓,在其上選取均勻分布的9個點作為測量點,對其殘余壓應力求平均值,作為有限元模擬的結果。最終得到的殘余壓應力分布如圖9(a)所示,其中FEA表示有限元模擬結果,Experiment表示實驗結果。模擬得到的殘余壓應力最大值為383MPa,比實驗值350MPa大9.4%;在誤差范圍內,但殘余壓應力影響深度達到0.95mm,比實驗值大46%,嚴重超出了誤差范圍,并且在深度方向上模擬得到的值均比實驗值大,只有個別數據能夠落在實驗所測數據的誤差范圍之內,說明預估的材料動態屈服強度比實際值小,使材料過早進入屈服狀態,塑性變形加大,造成殘余壓應力變大。改變材料的動態屈服強度為850MPa,再次進行模擬,其結果如圖9(b)所示,殘余壓應力最大值為300MPa,影響深度為0.68mm,與實驗值的誤差分別為14.2%和4.6%,殘余壓應力最大值超出誤差范圍,并且壓應力在深度方向上均比實驗值小,且大部分數據都落在實驗數據的誤差范圍之外,說明所選的動態屈服強度大于實際值,使塑性變形比實際值要小,造成殘余壓應力減小。因此,TiAl合金的動態屈服強度應在700~850MPa之間,取為中間值775MPa進行模擬,所得結果如圖9(c)所示,殘余壓應力最大值為326MPa,比實驗值小6.8%;影響深度為0.7mm,比實驗值大7.7%,均在誤差范圍內。并且在0.15~0.7mm內,模擬值與實驗值吻合的非常好,說明TiAl合金的動態屈服強度為775MPa是可信的,由此得到的HEL為1.46GPa。但實驗所得的殘余壓應力最大值出現在次表層,而模擬所得的最大值出現在表層,殘余壓應力最大值出現的位置不相符,造成這種現象的原因主要有兩方面:一方面,由于有限元模擬僅僅考慮了激光沖擊強化的力學效應,而忽略了其熱效應。在實驗中,吸收層在激光輻照下會迅速氣化并等離子化,在該過程中,等離子化的氣體溫度高達數千攝氏度,若吸收層太薄,則等離子化的氣體會對試件表面產生熱影響,使材料表面溫度快速升高,當沖擊完畢后,冷卻水會快速帶走大部分熱量,使材料快速冷卻收縮并在材料表層形成拉應力,這將抵消掉材料表面部分殘余壓應力,使表層殘余壓應力減小;另一方面,當沖擊波壓力為2倍的HEL時,材料表層塑性變形將達到飽和,表層殘余壓應力達到最大值,繼續增大沖擊波壓力,表層殘余壓應力將會減小,殘余壓應力最大值將會向材料內部移動,造成材料表層殘余壓應力的減小[6]。由于模擬中使用的沖擊波峰值壓力是2.65GPa,小于2倍的HEL,因此,可能是熱效應造成了實驗中表層殘余壓應力的減小,但具體原因仍需后續深入研究。

3.4本構模型參數獲得的方法在3.3節中得到TiAl合金的HEL為1.46GPa,由于目前尚無TiAl合金在高應變率下的材料本構模型參數,因此,無法對該結果進行直接驗證,但從圖9(c)可以看出,模擬結果與實驗結果在趨勢與數值上都是非常吻合的,證明該結果是可信的,該參數獲得方法是有效的。由于該方法對材料沒有特殊要求,只要能獲得激光沖擊強化實驗數據,且選定合理的誤差范圍,即可對任何具有單一參數的本構模型進行很好的確定。但是對于材料在高應變率下的其他本構模型,如Johnson-cook模型,本構關系如式(4)所示。

4結論

1)提出了一種基于激光沖擊強化實驗和有限元模擬相結合的材料本構模型參數獲得方法,該方法能對材料在高應變率下具有單一參數的本構模型進行很好的確定,對需要確定多參數的本構模型有一定局限性。2)建立激光沖擊強化有限元模擬的分析模型,并提出材料本構模型參數的預估方法。3)提出本構模型構建方法中誤差范圍的確定辦法,得出了TiAl合金在高應變率下的HEL為1.46GPa。

參考文獻

[1]PEYREP,FABBROR,MERRIENP,LIEURADEHP.Lasershockprocessingofaluminiumalloys.Applicationtohighcyclefatiguebehaviour[J].MaterialsScience&Engineering,1996,210(1):102113.

[2]喬紅超,高宇,趙吉賓,陸瑩,趙亦翔.激光沖擊強化技術的研究進展[J].中國有色金屬學報,2015,25(7):17441755.QIAOHong-chao,GAOYu,ZHAOJi-bin,LUYing,ZHAOYi-xiang.Researchprocessoflaserpeeningtechnology[J].TheChineseJournalofNonferrousMetals,2015,25(7):17441755.

[3]HUYong-xiang,LIKang-mei,QIChen-jie,YAOZQ,GRANDHIRV.Sizeeffectonindentationdepthofoxygen-freehighpuritycopperinducedbylasershockprocessing[J].TransactionsofNonferrousMetalsSocietyofChina,2012,22(S2):s573s578.

[4]張青來,鮑士喜,王榮,錢陽,張永康,李興成.激光沖擊強化對AZ31和AZ91鎂合金表面形貌和電化學腐蝕性能的影響[J].中國有色金屬學報,2014,24(10):24652473.ZHANGQing-lai,BAOShi-xi,WANGRong,QIANYang,ZHANGYong-kang,LIXing-cheng.EffectoflasershockprocessingonsurfacemorphologyandelectrochemicalcorrosionresistanceofAZ31andAZ91alloys[J].TheChineseJournalofNonferrousMetals,2014,24(10):24652473.

[5]ZHANGXing-quan,CHENLiu-san,YUXiao-liu,ZUOLi-sheng,ZHOUYu.Effectoflasershockprocessingonfatiguelifeoffastenerhole[J].TransactionsofNonferrousMetalsSocietyofChina,2014,24(4):969974.

[6]李應紅.激光沖擊強化理論與技術[M].北京:科學出版社,2013:20,133.LIYing-hong.Lasershocktheoryandtechnology[M].Beijing:SciencePress,2013:20,133.

[7]FABBROR,PEYREP,BERTHEL.Physicsandapplicationsoflaser-shockprocessing[J].JournalofLaserApplications,1998,10(6):265.

[8]MEYERSMA.材料的動力學行為[M].張慶明,劉彥,黃風雷,呂中杰,譯.北京:國防工業出版社,2006:225.MEYERSMA.Dynamicbehaviorofmaterials[M].ZHANGQing-ming,LIUYan,HUANGFeng-lei,LÜZhong-jie,transl.Beijing:NationalDefenceIndustryPress,2006:225.

[9]HOGGATTCR,RECHTRF.Stress-straindataobtainedathighratesusinganexpandingring[J].ExperimentalMechanics,1969,9(10):441448.

[10]趙峰,李玉龍,索濤,黃衛東,劉建睿.高應變率下鑄造鎂合金AZ91的動態壓縮性能及破壞機理[J].中國有色金屬學報,2009,19(7):11631168.ZHAOFeng,LIYu-long,SUOTao,HUANGWei-dong,LIUJian-rui.DynamiccompressivebehavioranddamagemechanismofcastmagnesiumalloyAZ91[J].TheChineseJournalofNonferrousMetals,2009,19(7):11631168.

[11]RENYu,TANCheng-wen,ZHANGJing,WANGFu-chi.DynamicfractureofTi-6Al-4ValloyinTaylorimpacttest[J].TransactionsofNonferrousMetalsSocietyofChina,2011,21(2):223235.

[12]DINGK,YEL.Lasershockpeeningperformanceandprocesssimulation[M].NewYork:CrcPress,2006:89,109.

[13]花銀群,蔡崢嶸,陳瑞芳,姜輝,吉光.TC4鈦合金激光搭接沖擊強化的實驗和數值模擬[J].激光技術,2010,34(5):632635,639.HUAYin-qun,CANZheng-rong,CHENRui-fang,JIANGHui,JIGuang.ExperimentandnumericalsimulationofoverlappinglasershockprocessinginTC4titaniumalloy[J].LaserTechnology,2010,34(5):632635,639.

[14]胡永祥.激光沖擊處理工藝過程數值建模與沖擊效應研究[D].上海:上海交通大學,2008.HUYong-xiang.Researchonthenumericalsimulationandimpacteffectsoflasershockprocessing[D].Shanghai:ShanghaiJiaoTongUniversity,2008.

[15]FABBROR,FOURNIERJ,BALLARDP,DEVAUXD,VIRMONTJ.Physicalstudyoflaser-producedplasmainconfinedgeometry[J].JournalofAppliedPhysics,1990,68(2):775.

[16]戴峰澤.基于納秒激光沖擊波效應的金屬表面形貌與性能研究[D].鎮江:江蘇大學,2014.DAIFeng-ze.Mechanismresearchofnano-secondlaserinducedshockwaveonthesurfacetopographyandperformanceofmetals[D].Zhenjiang:JiangsuUniversity,2014.

[17]柳永寧,朱金華,周惠久.普通低碳鋼的強度與溫度、應變速率的關系[J].材料科學進展,1990,4(4):285290.LIUYong-ning,ZHUJin-hua,ZHOUHui-jiu.Strengthofalowcarbonsteelatdifferenttemperaturesandstrainrates[J].MaterialsScienceProcess,1990,4(4):285290.

[18]李興華,楊紹利.鈦鋁合金制備技術現狀及新進展[J].材料導報,2011,25(7):94100.LIXing-hua,YANGShao-li.TheprogressofpreparationandprocessingtechnologyofTiAlalloy[J].MaterialsReview,2011,25(7):94100.

[19]喬紅超,趙亦翔,趙吉賓,陸瑩.激光沖擊強化對TiAl合金組織和性能的影響[J].光學精密工程,2014,22(7):17661773.QIAOHong-chao,ZHAOYi-xiang,ZHAOJi-bin,LUYing.EffectoflaserpeeningonmicrostructuresandpropertiesofTiAlalloy[J].OpticsandPrecisionEngineering,2014,22(7):17661773.

[20]VOOTHALURUR,LIUCR,CHENGGJ.Finiteelementanalysisofthevariationinresidualstressdistributioninlasershockpeeningofsteels[J].JournalofManufacturingScienceandEngineering,2012,134(6):30203023.

[21]胡永祥,姚振強,胡俊.激光沖擊強化殘余應力場的數值仿真分析[J].中國激光,2006,33(6):846851.HUYong-xiang,YAOZhen-qiang,HUJun.Numericalsimulationofresidualstressfieldforlasershockprocessing[J].ChineseJournalofLasers,2006,33(6):846851.

[22]郭乃國.激光沖擊處理40Cr鋼及其殘余應力場數值模擬[D].鎮江:江蘇大學,2007.GUONai-guo.Lasershockprocessingof40CrsteelandnumericalsimulationofresidualstressfieldinducedbyLSP[D].Zhenjiang:JiangsuUniversity,2007.

[23]武敬偉.激光沖擊Fe-Ni合金殘余應力場的數值模擬[D].鎮江:江蘇大學,2007.WUJing-wei.NumericalsimulationofresidualstressfieldinducedbylasershockprocessingonFe-Nialloy[D].Zhenjiang:JiangsuUniversity,2007.

作者:陳松林 趙吉賓 喬紅超 楊灝 單位:中國科學院 沈陽自動化研究所 中國科學院大學

主站蜘蛛池模板: 欧美系列第一页 | 一本岛一区在线观看不卡 | 伊人色在线观看 | 欧美日韩国产亚洲一区二区 | 丁香六月婷婷激情 | 亚洲精品国产成人99久久 | 国产精品久久久久一区二区三区 | 在线观看国产一区二区三区99 | 99久久精品国产一区二区 | 日本欧美一区二区三区片 | 在线观看你懂得 | 樱花草在线播放 | 羞羞视频免费网站入口 | 水蜜桃视频在线观看免费 | 91中文字字幕乱码 | 偷拍第一页| 手机视频在线 | 国产在线高清一级毛片 | 亚洲欧美视频网站 | 在线成人亚洲 | 成人性色生活影片 | 中文国产成人精品少久久 | 婷婷综合缴情亚洲五月伊 | 国产精品视频久 | 狠狠色丁香婷婷综合激情 | 欧美精品在线视频观看 | 天堂网avtt| 国色天香社区在线观看免费播放 | 国产精品免费在线播放 | 亚洲国产精品lv | 中国精品视频一区二区三区 | 欧美日韩在线视频观看 | 五月花在线视频 | 骚骚网站| 亚洲人成电影院 | 自w时看的视频 | 最新国产在线观看福利 | 自拍偷拍第1页 | 伊人久久综合 | 亚洲欧美一区二区三区 | 性视频网站在线 |