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

單晶高溫疲勞損傷參量的選取與壽命建模

2016-12-06 07:07:33荊甫雷王榮橋胡殿印蔣康河
航空學報 2016年9期
關鍵詞:模型

荊甫雷*,王榮橋,胡殿印,蔣康河

1.中航空天發動機研究院有限公司 基礎與應用研究中心,北京 101304 2.北京航空航天大學 能源與動力工程學院,北京 100083

單晶高溫疲勞損傷參量的選取與壽命建模

荊甫雷1,*,王榮橋2,胡殿印2,蔣康河2

1.中航空天發動機研究院有限公司 基礎與應用研究中心,北京 101304 2.北京航空航天大學 能源與動力工程學院,北京 100083

高溫疲勞損傷是引起單晶渦輪葉片破壞的主要因素之一。利用不同試驗條件下DD6標準試件的低周疲勞和蠕變-疲勞試驗結果,結合基于滑移系的黏塑性應力-應變分析,分別研究了晶體取向、應變范圍、平均應變以及保載時間等對單晶高溫疲勞損傷的影響機制。進而采用滑移剪應變最大的滑移系作為臨界滑移系,選取臨界滑移系上的最大Schmid應力、最大滑移剪應變率、循環Schmid應力比以及滑移剪應變范圍等細觀參量作為損傷參量,建立了一種新的基于臨界平面的循環損傷累積(CDA)模型。結果表明,該模型對于DD6高溫疲勞壽命預測精度基本在3倍分散帶內。

鎳基單晶高溫合金;高溫疲勞;損傷;壽命預測;臨界平面;滑移系

鎳基單晶高溫合金具有優異的耐高溫、抗氧化、抗蠕變以及抗疲勞性能,是目前制造先進航空發動機渦輪轉子葉片的主要材料[1]。單晶渦輪葉片在服役過程中承受著高溫環境和交變應力的聯合作用,高溫疲勞損傷是引起其破壞的主要因素之一[2-3]。已有的 試 驗 結 果[4-5]表 明,單 晶 的 高 溫疲勞壽命與晶體取向、試驗溫度、應變范圍、平均應變和保載時間等因素密切相關。

高溫疲勞損傷的晶體取向相關性是單晶與各向同性合金的最主要差異。文獻[6-8]認為彈性模量的方向性變化是引起上述差異的主要原因,并將晶體彈性模量及其取向參數引入到各向同性壽命模型中,建立了一系列唯象的單晶高溫疲勞壽命模型。由于直接采用宏觀應力、應變和保載時間等試驗測量值作為損傷參量,而不涉及對單晶細觀變形過程的描述,此類模型難以反映單晶疲勞裂紋通常沿滑移面萌生的細觀特征[5],也無法慮及本構模型的計算誤差對壽命預測的影響。

近年來,用于多軸疲勞壽命預測的臨界平面法(Critical Plane Approach)[9]開始被引入到單晶高溫疲勞壽命建模中[10-13]。滑移是單晶主要的非彈性變形機制,高溫疲勞損傷主要在特定的滑移面上累積并引起裂紋萌生和擴展,最容易破壞的平面是發生滑移的特定滑移面[5]。采用特定滑移面作為臨界平面,并以該面上線彈性應力-應變分析得到的滑移系細觀參量(如滑移系上的剪應力/應變、正應力/應變等)作為損傷參量進行壽命建模,在解決單晶高溫疲勞損傷晶體取向相關性的同時,還可以很好地反映其破壞的細觀特征。需要注意的是,上述模型中普遍缺少表征晶體非彈性變形的相關項(如滑移剪應變),因而無法準確體現非彈性變形對壽命的影響,目前主要用于單晶高周疲勞壽命預測,而對于低周疲勞和蠕變-疲勞壽命預測精度較低[14]。

另一方面,Levkovitch[15]和 Tinga[16]等從單晶滑移面上微裂紋擴展過程出發,采用與Majumdar損傷率模型相類似的形式,以各個滑移系上的Schmid應力與滑移剪應變率作為損傷參量,提出了單晶損傷率模型。其中,滑移剪應變率需要借助于單晶黏塑性本構模型進行獲取,這確保了該模型能夠考慮本構模型的計算誤差對壽命預測的影響;而基于滑移系的壽命建模則可以很好地描述單晶高溫疲勞損傷的晶體取向相關性。從基于滑移系的黏塑性本構模型方程[15-16]可以看出,滑移剪應變率與Schmid應力直接相關,特別是在背應力達到飽和時,二者呈指數函數關系,即此時損傷只與Schmid應力這一個參量相關,因此對于具有明顯非彈性變形和非對稱循環的情況,單晶損傷率模型的壽命預測誤差可能較大[17]。

在上述研究的基礎上,文獻[18]在循環損傷累積(Cyclic Damage Accumulation,CDA)模型[8]中引入臨界平面的概念,并采用單晶黏塑性本構模型計算得到的最大Schmid應力、最大滑移剪應變率、滑移剪應變范圍以及宏觀應變比、拉伸/壓縮保載頻率等作為損傷參量,所得壽命模型預測精度較高。但該模型中仍包含宏觀測量值,未能準確描述平均應變和保載時間等對單晶高溫疲勞損傷的影響機制,導致其適用范圍有限,無法滿足單晶渦輪葉片壽命評估的要求。

鑒于此,本文以國產第二代鎳基單晶高溫合金DD6為對象,分別研究晶體取向、應變范圍、平均應變以及保載時間等對單晶高溫疲勞損傷的影響機制;在此基礎上,結合已有模型的優點,提出一種完全以黏塑性應力應變分析得到的滑移系細觀參量作為損傷參量、能夠考慮本構模型計算誤差并反映破壞細觀特征及損傷晶體取向相關性的單晶高溫疲勞壽命模型,進而完成材料常數的獲取與驗證。

1 DD6高溫疲勞試驗數據

目前已有的DD6標準試件高溫疲勞試驗數據如表1所示[4-5,17]。其 中,760 ℃ 下 的試 驗 數 據比較全面,包含不同的晶體取向、應變范圍、應變比以及拉伸(壓縮)保載時間等多種情況。因此本文首先采用760℃下的DD6高溫疲勞試驗數據,進行損傷參量的選取,并建立單晶高溫疲勞壽命模型;然后利用該模型對980℃下DD6在[001]取向的低周疲勞和蠕變-疲勞進行壽命預測,確定其在不同溫度下的壽命預測精度;最后補充不同加載速率下的DD6標準試件高溫疲勞試驗,對該模型進行了初步驗證。

表1 DD6標準試件高溫疲勞試驗數據Table 1 High temperature fatigue test data of DD6standard specimens

2 損傷參量的選取

2.1 晶體取向和應變范圍對損傷的影響

利用 DD6在760 ℃下[001]/[011]/[111]取向、應變比為-1、無保載時間的低周疲勞試驗數據,研究晶體取向和應變范圍對單晶高溫疲勞損傷的影響機制。采用基于滑移系的Walker黏塑性本構模型[18-20],進行上述試驗條件下的應力-應變分析,得到滑移系細觀參量隨時間的變化,如圖1所示。其中,Oct表示八面體滑移系,Cube表示六面體滑移系。以基本達到穩定的第2個循環作為計算循環,確定各個滑移系上的典型損傷參量,如單晶損傷率模型所采用的最大Schmid應力和最大滑移剪應變率[16]、與單晶非彈性變形量直接相關的滑移剪應變范圍[18]等。上述損傷參量與對稱循環低周疲勞壽命間的關系如圖2所示。

通過圖2可以看出,在雙對數坐標系中,最大Schmid應力、最大滑移剪應變率以及滑移剪應變范圍與低周疲勞壽命近似呈線性關系。需要注意的是,在應變范圍較大、壽命較低(≤102)的情況下,材料進入屈服階段,此時最大Schmid應力和最大滑移剪應變率基本不變,無法描述不同應變范圍下低周疲勞壽命的差異;而滑移剪應變范圍對于壽命的變化則比較敏感。因此,在壽命模型中引入滑移剪應變范圍作為損傷參量可以有效提高低周疲勞的壽命預測精度。

圖1 不同取向滑移系細觀參量隨時間的變化Fig.1 Microscopic parameters on slip systems vs time in different orientations

圖2 典型損傷參量與壽命的關系Fig.2 Relationship between typical damage parameters and life

此外,在[001]取向受載時,只有八面體滑移系上存在非零的Schmid應力、滑移剪應變等細觀參量,因此其損傷完全由八面體滑移系貢獻;在[011]取向受載時,八面體滑移系和六面體滑移系上均存在細觀參量,且八面體滑移系上的細觀參量大于六面體滑移系,因此其損傷主要由八面體滑移系貢獻;而在[111]取向受載時,六面體滑移系上的細觀參量明顯大于八面體滑移系,因此其損傷主要由六面體滑移系貢獻。

綜上所述,單晶高溫疲勞損傷的晶體取向相關性主要體現在損傷最大的滑移系上,而應變范圍對單晶高溫疲勞損傷的影響主要體現在最大Schmid應力、最大滑移剪應變率和滑移剪應變范圍等損傷參量的差異上。

2.2 平均應變對損傷的影響

利用DD6在760℃下、[001]取向、應變比為0.05、無保載時間的低周疲勞試驗數據,研究平均應變對單晶高溫疲勞損傷的影響機制。通過黏塑性應力-應變分析得到不同平均應變下滑移系細觀參量隨時間的變化,如圖3所示。

通過[001]取向對稱循環(見圖1)與非對稱循環(見圖3)的滑移系細觀參量對比可以看出,二者的差異主要源于滑移系上循環Schmid應力比(最小Schmid應力/最大Schmid應力)的明顯不同:對稱循環Schmid應力比近似為宏觀應變比-1,如圖1(a)所示;而非對稱循環Schmid應力比隨著平均應變的減小而顯著增大,并逐漸趨近于宏觀應變比0.05,如圖3(a)所示。因此,循環Schmid應力比可能是表征平均應變對單晶高溫疲勞損傷影響的有效參量。

圖3 不同平均應變下滑移系細觀參量隨時間的變化Fig.3 Microscopic parameters on slip systems vs time under different average strains

2.3 保載時間對損傷的影響

利用DD6在760℃下、[001]取向、應變比為-1、不同拉伸(壓縮)保載時間(60(0)、0(60)、30(30))的蠕變-疲勞試驗數據,研究保載時間對單晶高溫疲勞損傷的影響機制。通過黏塑性應力-應變分析得到相同應變范圍、不同保載時間下滑移系細觀參量隨時間的變化,如圖4所示。通過[001]取向無保載時間(圖1)與不同保載時間(圖4)的滑移系細觀參量對比可以看出,二者的最大Schmid應力、最大滑移剪應變率、循環Schmid應力比等損傷參量基本相同。保載時間對于疲勞壽命的影響主要體現在滑移剪應變范圍上:保載時間的存在導致了滑移剪應變范圍的明顯增加,并且對于不同的拉伸(壓縮)保載情況,滑移剪應變范圍也存在一定的差異。因此,滑移剪應變范圍可能是表征保載時間對單晶高溫疲勞損傷影響的有效參量。

圖4 不同拉伸(壓縮)保載時間下滑移系細觀參量隨時間的變化Fig.4 Microscopic parameters on slip systems vs time under different tension(compression)dwell time

3 壽命建模

綜合考慮晶體取向、應變范圍、平均應變以及保載時間等因素對單晶高溫疲勞損傷的影響,在單晶損傷率模型[16]中引入臨界平面[10-14]的概念,以滑移剪應變最大的滑移系作為臨界滑移系,其所在的滑移面作為首先發生疲勞破壞的臨界平面;采用與CDA模型[8,18]相類似的形式,建立單晶高溫疲勞壽命與臨界滑移系上細觀損傷參量的函數關系,得到一種基于臨界平面的CDA模型,其為式中:N為單晶高溫疲勞壽命;α為臨界滑移系的類型(Oct或Cube)和分別為臨界滑移系上在宏觀應力最大和最小時所對應的Schmid應力為臨界滑移系上的循環Schmid應力比和Δγα分別為臨界滑移系上的最大滑移剪應變率和滑移剪應變范圍;Aα、mα、nα、zα和aα為與溫度相關的材料常數。

對式(1)兩邊取對數后,分別針對臨界滑移系為八面體滑移系([001]取向和[011]取向)以及臨界滑移系為六面體滑移系([111]取向)的情況,進行多元線性回歸分析,確定其材料常數如表2所示。由于目前DD6在980℃下只有[001]取向的試驗數據,其臨界滑移系為八面體滑移系,本文只給出了在該溫度下的八面體滑移系材料常數。

基于臨界平面的CDA模型對760℃和980℃下DD6標準試件高溫疲勞壽命的預測結果如圖5和圖6所示,單晶損傷率模型[16]的預測結果如圖7所示。可以看出,基于臨界平面的CDA模型計算的壽命基本落在試驗壽命的3倍分散帶內,其精度顯著高于單晶損傷率模型。

表2 基于臨界平面的CDA模型材料常數Table 2 Material constants of CDA model based on critical plane

圖5 760℃下的疲勞壽命預測(基于臨界平面的CDA模型)Fig.5 Fatigue life prediction at 760 ℃ (CDA model based on critical plane)

圖6 980℃下[001]取向疲勞壽命預測(基于臨界平面的CDA模型)Fig.6 Fatigue life prediction in orientation [001]at 980℃(CDA model based on critical plane)

圖7 760℃下疲勞壽命預測(損傷率模型[16])Fig.7 Fatigue life prediction at 760℃ (Damage rate model[16])

4 模型驗證

為了進一步驗證前面建立的基于臨界平面的CDA模型,補充了DD6標準試件在不同加載速率下的高溫疲勞試驗。試驗溫度分別為760℃和980℃,采用機械應變幅控制,載荷譜為三角波,周期為3s,每個狀態進行2次試驗,如表3所示。

利用基于臨界平面的CDA模型進行上述試驗條件下的高溫疲勞壽命預測,結果如圖8所示。

表3 壽命模型驗證試驗數據Table 3 Test data for verification of life model

圖8 基于臨界平面的CDA模型驗證Fig.8 Verification of CDA model based on critical plane

可以看出,不同加載速率下的壽命預測精度均在3倍分散帶內,表明該模型可以比較準確地描述單晶高溫疲勞損傷行為。

5 結 論

1)單晶高溫疲勞損傷所呈現出的晶體取向相關性主要體現在損傷最大的臨界滑移系上,對于[001]取向和[011]取向為八面體滑移系,而[111]取向為六面體滑移系;應變范圍對單晶高溫疲勞損傷的影響主要體現在最大Schmid應力、最大滑移剪應變率、滑移剪應變范圍等細觀參量的差異上;而平均應變和保載時間對單晶高溫疲勞損傷的影響則可以分別通過循環Schmid應力比和滑移剪應變范圍來表征。

2)基于臨界平面的CDA模型以滑移剪應變最大的滑移系作為臨界滑移系,能夠反映單晶高溫疲勞破壞的細觀特征和損傷的晶體取向相關性;采用臨界滑移系上的最大Schmid應力、最大滑移剪應變率、循環Schmid應力比以及滑移剪應變范圍等細觀參量作為損傷參量,能夠描述應變范圍、平均應變、保載時間以及加載速率等對單晶高溫疲勞損傷的影響,并可慮及本構模型的計算誤差;該模型對于DD6高溫疲勞壽命的預測結果基本落在試驗壽命的3倍分散帶內,能夠滿足單晶渦輪葉片壽命評估的要求。

[1] 胡壯麒,劉麗榮,金濤,等.鎳基單晶高溫合金的發展[J].航空發動機,2005,31(3):1-7.HU Z Q,LIU L R,JIN T,et al.Development of the Nibased single crystal superalloys[J].Aeroengine,2005,31(3):1-7(in Chinese).

[2] 陶春虎,鐘培道,王仁智,等.航空發動機轉動部件的失效與預防[M].北京:國防工業出版社,2008:46-83.TAO C H,ZHONG P D,WANG R Z,et al.Failure analysis and prevention for rotor in aero-engine[M].Beijing:National Defense Industry Press,2008:46-83 (in Chinese).

[3] 丁智平,劉義倫,尹澤勇.鎳基單晶高溫合金蠕變-疲勞壽命評估方法進展[J].機械強度,2003,25(3):254-259.DING Z P,LIU Y L,YIN Z Y.Development on evaluating method of creep-fatigue life of single crystal nickelbased superalloys[J].Journal of Mechanical Strength,2003,25(3):254-259(in Chinese).

[4] 李影,蘇彬,吳學仁.高溫下取向對DD6單晶高溫合金低周疲勞壽命的影響[J].航空材料學報,2001,21(2):22-25.LI Y,SU B,WU X R.Orientation dependence of low cycle fatigue life of single-crystal nickel-base superalloy DD6 under high temperature[J].Journal of Aeronautical Materials,2001,21(2):22-25(in Chinese).

[5] 李影,蘇彬.DD6單晶合金高溫低周疲勞機制[J].航空動力學報,2003,18(6):732-736.LI Y,SU B.Mechanisms of low cyclic fatigue of DD6alloy at elevated temperature[J].Journal of Aerospace Power,2003,18(6):732-736(in Chinese).

[6] LI S X,ELLISON E G,SMITH D J.The influence of orientation on the elastic and low cycle fatigue properties of several single crystal nickel base superalloys[J].Journal of Strain Analysis for Engineering Design,1994,29(2):147-153.

[7] 岳珠峰,陶仙德,尹澤勇,等.一種鎳基單晶超合金高溫低周疲勞的晶體取向相關性模型[J].應用數學和力學,2000,21(4):373-381.YUE Z F,TAO X D,YIN Z Y,et al.A crystallographic model for the orientation dependence of low cyclic fatigue property of a nickel-base single crystal superalloy[J].Applied Mathematics and Mechanics,2000,21(4):373-381(in Chinese).

[8] 石多奇,楊曉光,于慧臣.一種鎳基單晶和定向結晶的疲勞壽命模型[J].航空動力學報,2010,25(8):1871-1875.SHI D Q,YANG X G,YU H C.Fatigue life prediction model for nickel-based single crystal and directionally solidified superalloy[J].Journal of Aerospace Power,2010,25(8):1871-1875(in Chinese).

[9] KAROLCZUK A,MACHA E.A review of critical plane orientations in multiaxial fatigue failure criteria of metallic materials[J].International Journal of Fracture,2005,134(1):267-304.

[10] ARAKERE N K,SWANSON G.Effect of crystal orientation on fatigue failure of single crystal nickel base turbine blade superalloys[J].Journal Engineering for Gas Turbines and Power,2002,124(1):161-175.

[11] SWANSON G,ARAKERE N K.Effect of crystal orientation on analysis of single-crystal nickel-based turbine blade superalloys:NASA TP-2000-210074[R].Washington,D.C.:NASA,2000:11-54.

[12] NAIK R A,DELUCA D P,SHAH D M.Critical plane fatigue modeling and characterization of single crystal nickel superalloys[J].Journal of Engineering for Gas Turbines and Power,2004,126(2):391-400.

[13] SHI D Q,HUANG J,YANG X G,et al.Effects of crystallographic orientations and dwell types on low cycle fatigue and life modeling of a SC superalloy[J].International Journal of Fatigue,2013,49(1):31-39.

[14] 劉金龍.鎳基單晶/定向凝固渦輪葉片鑄造模擬及其合金低循環疲勞行為研究[D].北京:北京航空航天大學,2011:101-130.LIU J L.Investigation on cast simulation and low cycle fatigue behavior of Ni-based single crystal and directionally solidified turbine blade[D].Beijing:Beihang University,2011:101-130(in Chinese).

[15] LEVKOVITCH V,SIEVERT R,SVEBDSEN B.Simulation of deformation and lifetime behavior of a FCC single crystal superalloy at high temperature under low-cycle fatigue loading[J].International Journal of Fatigue,2006,28(12):1791-1802.

[16] TINGA T,BREKELMANS W A M,GEERS M G D.Time-incremental creep-fatigue damage rule for single crystal Ni-base superalloys[J].Materials Science and Engineering A,2009,508(1-2):200-208.

[17] 荊甫雷.單晶渦輪葉片熱機械疲勞壽命評估方法研究[D].北京:北京航空航天大學,2013:59-75.JING F L.Research on thermo-mechanical fatigue life assessment of single crystal turbine blades[D].Beijing:Beihang University,2013:59-75(in Chinese).

[18] 王榮橋,荊甫雷,胡殿印.基于臨界平面的鎳基單晶高溫合金疲勞壽命預測模型[J].航空動力學報,2013,28(11):2587-2592.WANG R Q,JING F L,HU D Y.Fatigue life prediction method based on critical plane for nickel-based single crystal superalloys[J].Journal of Aerospace Power,2013,28(11):2587-2592(in Chinese).

[19] JORDAN E H,WALKER K P.A viscoplastic model for single crystals[J].Journal of Engineering Materials and Technology,1992,114(1):19-26.

[20] JORDAN E H,SHI S X,WALKER K P.The viscoplastic behavior of Hastelloy-X single crystal[J].International Journal of Plasticity,1993,9(1):119-139.

Damage parameter determination and life modeling for high temperature fatigue of single crystals

JING Fulei1,*,WANG Rongqiao2,HU Dianyin2,JIANG Kanghe2
1.Basic and Applied Research Center,AVIC Academy of Aeronautic Propulsion Technology,Beijing 101304,China 2.School of Energy and Power Engineering,Beihang University,Beijing 100083,China

High temperature fatigue damage is a major factor causing the failure of single crystal turbine blades.The influence mechanisms of crystal orientation,strain range,mean strain and dwell time on the high temperature fatigue damage of nickel-based single crystal superalloys are studied respectively with the results of low cycle fatigue and creep-fatigue tests on DD6standard specimens under different testing conditions and with viscoplastic stress-strain analysis based on slip systems.Furthermore,the slip system with the max slip shear strain is utilized as the critical slip system where the max Schmid stress,max slip shear strain rate,cyclic Schmid stress ratio and slip shear strain range are selected as the damage parameters,and a new cyclic damage accumulation(CDA)model based on critical plane is proposed.The results indicate that the predicted high temperature fatigue life of DD6with the proposed CDA model based on critical plane is basically within a factor three of the experimental life.

nickel-based single crystal superalloy;high temperature fatigue;damage;life prediction;critical plane;slip system

2015-08-31;Revised:2015-11-03;Accepted:2015-12-01;Published online:2015-12-08 14:22

URL:www.cnki.net/kcms/detail/11.1929.V.20151208.1422.002.html

s:National Natural Science Foundation of China(51375031);Aeronautical Science Foundation of China(2015ZBN3004)

V232.4;O346.2

A

1000-6893(2016)09-2749-08

10.7527/S1000-6893.2015.0326

2015-08-31;退修日期:2015-11-03;錄用日期:2015-12-01;網絡出版時間:2015-12-08 14:22

www.cnki.net/kcms/detail/11.1929.V.20151208.1422.002.html

國家自然科學基金(51375031);航空科學基金(2015ZBN3004)

*通訊作者.Tel.:010-56680678 E-mail:jingfulei@163.com

荊甫雷,王榮橋,胡殿印,等.單晶高溫疲勞損傷參量的選取與壽命建模[J].航空學報,2016,37(9):27492-756.JING F L,WANG R Q,HU D Y,et al.Damage parameter determination and life modeling for high temperature fatigue of single crystals[J].Acta Aeronautica et Astronautica Sinica,2016,37(9):27492-756.

荊甫雷 男,博士,工程師。主要研究方向:航空發動機熱端部件結構強度。

Tel.:010-56680678

E-mail:jingfulei@163.com

*Corresponding author.Tel.:010-56680678 E-mail:jingfulei@163.com

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 免费在线国产一区二区三区精品| 日本91视频| 97亚洲色综久久精品| 伊人无码视屏| 亚洲日韩国产精品无码专区| 热热久久狠狠偷偷色男同| 国产免费观看av大片的网站| 中文无码精品a∨在线观看| 亚洲午夜福利精品无码| 日韩国产精品无码一区二区三区| 国产成人凹凸视频在线| 精品视频在线观看你懂的一区| 午夜国产理论| 中文字幕久久波多野结衣| 国产成年女人特黄特色毛片免 | 亚洲欧美另类日本| 97人妻精品专区久久久久| 午夜不卡视频| 午夜精品福利影院| 色成人亚洲| 精品无码一区二区在线观看| 亚洲一区二区约美女探花| 67194在线午夜亚洲| 91破解版在线亚洲| 久久国产精品电影| 国产97公开成人免费视频| 日韩欧美综合在线制服| 成人噜噜噜视频在线观看| 五月六月伊人狠狠丁香网| 91亚洲精选| 亚洲黄网在线| 午夜日韩久久影院| 久久久久九九精品影院| 亚洲国产成人精品无码区性色| 久久精品国产国语对白| 日本午夜三级| 男人天堂亚洲天堂| 91小视频在线播放| 国产美女精品人人做人人爽| 扒开粉嫩的小缝隙喷白浆视频| 亚洲一区二区精品无码久久久| 无码日韩精品91超碰| 综合色区亚洲熟妇在线| 九九热精品视频在线| 中文字幕亚洲专区第19页| 亚洲色成人www在线观看| 久久综合成人| 国产精品永久久久久| 免费高清毛片| 92精品国产自产在线观看| 久久一本日韩精品中文字幕屁孩| 毛片卡一卡二| 免费看久久精品99| 国产成人一区| 久久亚洲国产最新网站| 婷婷成人综合| 无码精品国产dvd在线观看9久| 国产精品视频导航| 福利国产微拍广场一区视频在线| 亚洲AV无码乱码在线观看裸奔 | 国内精品伊人久久久久7777人| 69综合网| 久久精品国产精品青草app| 亚洲人成电影在线播放| 无遮挡一级毛片呦女视频| 亚洲无线视频| 欧美亚洲欧美区| 中国国产A一级毛片| 亚洲日韩精品无码专区| 久久久久亚洲精品成人网 | 亚洲熟女中文字幕男人总站| 国产日韩欧美在线视频免费观看| 广东一级毛片| 亚洲精品成人片在线观看| 在线观看欧美国产| 高潮毛片无遮挡高清视频播放| 98精品全国免费观看视频| 搞黄网站免费观看| 国产精品深爱在线| 亚洲成人动漫在线| 99一级毛片| 欧美色99|