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

基于J-C模型的GH907高溫合金動態(tài)本構(gòu)關(guān)系及失效關(guān)系

2021-10-28 07:14:30董澤民劉璐璐徐凱龍趙振華
機(jī)械工程材料 2021年10期
關(guān)鍵詞:模型

董澤民,陳 偉,劉璐璐,徐凱龍,趙振華

(南京航空航天大學(xué)1.能源與動力學(xué)院,航空發(fā)動機(jī)熱環(huán)境與熱結(jié)構(gòu)工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室,2.機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室,南京 210016)

0 引 言

GH907高溫合金是一種以鐵-鎳-鈷為基的低膨脹高溫合金,在650 ℃以下具有很高的強(qiáng)度、低的膨脹系數(shù)、良好的熱疲勞性能以及幾乎恒定不變的彈性模量,因而在制造航空發(fā)動機(jī)機(jī)匣、隔熱環(huán)等環(huán)形件上應(yīng)用廣泛[1]。航空發(fā)動機(jī)在運(yùn)行過程中,其機(jī)匣容易受到葉片故障產(chǎn)生的高動能碎片的沖擊,給飛行安全帶來隱患。而在研制階段,發(fā)動機(jī)機(jī)匣包容試驗(yàn)成本高,耗時長,因此數(shù)值模擬成為非常重要的研究方法。準(zhǔn)確描述機(jī)匣材料在高溫、高應(yīng)變速率下的動態(tài)本構(gòu)關(guān)系以及損傷失效關(guān)系對于航空發(fā)動機(jī)機(jī)匣的包容數(shù)值模擬至關(guān)重要。目前,在眾多常用的描述金屬材料的動態(tài)本構(gòu)模型及失效模型中,Johnson-Cook(J-C)模型[2-3]因形式簡潔、物理意義明確并且可以通過分離變量的方式標(biāo)定參數(shù)而被學(xué)者們廣泛采用。

國內(nèi)外學(xué)者已經(jīng)對航空發(fā)動機(jī)常用高溫合金的動態(tài)力學(xué)性能進(jìn)行了廣泛的研究。劉曉等[4]對GH4169合金進(jìn)行了J-C本構(gòu)模型和失效模型的擬合,對本構(gòu)模型中的對數(shù)應(yīng)變速率項(xiàng)系數(shù)進(jìn)行了線性修正,并通過動態(tài)壓縮試驗(yàn)和仿真驗(yàn)證了模型的有效性。LIU等[5]采用2套不同參數(shù)的J-C模型描述了不同溫度范圍GH4169合金的本構(gòu)關(guān)系,并通過數(shù)值仿真模擬了彈丸沖擊靶板過程,結(jié)果表明仿真預(yù)測的彈道極限及失效模式與試驗(yàn)結(jié)果具有較好的一致性。KORKMAZ等[6]通過對Nimonic 80A高溫合金進(jìn)行不同溫度、不同應(yīng)變速率的拉伸試驗(yàn)獲取了J-C模型參數(shù),并通過數(shù)值仿真模擬動態(tài)拉伸試驗(yàn)驗(yàn)證了參數(shù)的準(zhǔn)確性。HAQ等[7]采用J-C模型描述Inconel-718合金靶板的本構(gòu)和失效行為,利用數(shù)值仿真方法探究了不同頂角的圓錐形彈體對靶板破壞模式和彈道極限的影響。WANG等[8]提出修正的J-C本構(gòu)模型和失效模型來表征發(fā)動機(jī)包容環(huán)常用GH3536高溫合金的動態(tài)力學(xué)行為,并利用數(shù)值仿真對該材料蜂窩結(jié)構(gòu)在不同溫度、沖擊速度和沖擊角度下的抗沖擊性能進(jìn)行了分析,結(jié)果表明修正的模型比原始J-C模型具有更準(zhǔn)確的表征能力。UGODILINWA等[9]研究了新型航空高溫合金Haynes 282在3種不同熱處理?xiàng)l件下的準(zhǔn)靜態(tài)和高應(yīng)變速率壓縮變形行為,建立了Arrhenius和改進(jìn)的J-C本構(gòu)模型來描述材料在高應(yīng)變速率和高溫下的動態(tài)力學(xué)行為。BORA等[10]和WANG等[11]分別在J-C失效模型中引入不同表達(dá)式的Lode應(yīng)變參數(shù)來表征材料的動態(tài)失效行為。

目前,有關(guān)GH907合金的研究主要集中在表面涂層的抗氧化性能[12]、耐腐蝕性能[13]、熱處理工藝和顯微組織[14]等方面,關(guān)于該合金在不同溫度、不同應(yīng)變速率下的動態(tài)力學(xué)行為研究并不多見。豐建朋等[15]在分析3種Arrhenius型方程的基礎(chǔ)上,以Zener-Hollomon參數(shù)為主要變量,并綜合考慮應(yīng)變、應(yīng)變速率和溫度的影響,建立了GH907合金在應(yīng)變速率10-2~10 s-1范圍內(nèi)的本構(gòu)關(guān)系;但是該應(yīng)變速率范圍較小,且未考慮失效模型的影響。作者以GH907合金為研究對象,使用材料試驗(yàn)機(jī)在常溫下進(jìn)行光滑圓棒和缺口試樣的準(zhǔn)靜態(tài)拉伸試驗(yàn),并使用分離式霍普金森拉桿試驗(yàn)裝置進(jìn)行了應(yīng)變速率在1 000~3 000 s-1的動態(tài)拉伸試驗(yàn),還使用分離式霍普金森壓桿(SHPB)試驗(yàn)裝置進(jìn)行了溫度在20~400 ℃、應(yīng)變速率在600~3 000 s-1的動態(tài)壓縮試驗(yàn),獲得了用于描述GH907合金在較寬溫度和應(yīng)變速率范圍內(nèi)的力學(xué)性能的J-C模型參數(shù),并通過SHPB仿真驗(yàn)證本構(gòu)模型參數(shù)的有效性,擬為打靶和機(jī)匣包容數(shù)值模擬提供依據(jù)。

1 試驗(yàn)材料與方法

試驗(yàn)材料為GH907高溫合金,其化學(xué)成分[1]如表1所示。

表1 GH907高溫合金的化學(xué)成分

在100KN電子萬能材料試驗(yàn)機(jī)上進(jìn)行室溫和高溫準(zhǔn)靜態(tài)拉伸試驗(yàn)。室溫準(zhǔn)靜態(tài)拉伸用光滑圓棒試樣和缺口試樣的形狀和尺寸分別見圖1(a)和圖1(b),為了獲得較寬應(yīng)力三軸度范圍內(nèi)合金的力學(xué)性能,設(shè)計了4種缺口尺寸,缺口半徑R分別為1.0,1.5,3.0,6.0 mm。通過速率控制的方式加載,光滑試樣的名義加載應(yīng)變速率分別為0.000 1,0.001,0.01 s-1,缺口試樣的名義加載應(yīng)變速率為0.001 s-1。高溫準(zhǔn)靜態(tài)拉伸用光滑圓棒試樣的形狀和尺寸見圖1(c),通過管式加熱爐對試樣及其轉(zhuǎn)接段加熱,外置程序控制爐溫及保溫時間,溫度分別為80,160,240,320 ℃,保溫時間為1 h,應(yīng)變速率為10-3s-1。利用直徑20 mm的分離式霍普金森拉桿試驗(yàn)裝置進(jìn)行室溫動態(tài)拉伸試驗(yàn),動態(tài)拉伸試樣的形狀和尺寸見圖1(d),夾持端與拉桿進(jìn)行連接,并通過應(yīng)變片、示波器等裝置采集試驗(yàn)中的波形信號,試驗(yàn)時的應(yīng)變速率分別為1 700 s-1,2 100 s-1,2 700 s-1。通過壓桿直徑14.5 mm的SHPB試驗(yàn)裝置進(jìn)行室溫(20 ℃)和高溫動態(tài)壓縮試驗(yàn),撞擊桿長0.4 m,動態(tài)壓縮試樣的形狀和尺寸如圖1(e)所示。高溫動態(tài)壓縮試驗(yàn)中試樣的溫度通過熱電偶測定加熱爐爐膛中心溫度代替,試驗(yàn)溫度分別為200,400 ℃,通過應(yīng)變片、動態(tài)應(yīng)變儀等裝置采集試驗(yàn)中的波形信號,應(yīng)變速率分別為600 s-1,1 000 s-1,1 400 s-1,1 700 s-1,1 800 s-1,2 000 s-1。

圖1 試樣的形狀和尺寸Fig.1 Geometry and dimensions of specimens: (a) room temperature quasi-static tension smooth specimen; (b) room temperature notched quasi-static tension specimen; (c) high temperature quasi-static tension smooth specimen; (d) dynamic tension specimen and (e) dynamic compression specimen

2 試驗(yàn)結(jié)果與討論

2.1 拉伸性能

由圖2可以看出,GH907合金具有明顯的屈服平臺,是一種典型的韌性金屬材料。在準(zhǔn)靜態(tài)拉伸試驗(yàn)中,盡管應(yīng)變速率不同,但試驗(yàn)結(jié)果仍具有較好的重復(fù)性。

圖2 光滑試樣的準(zhǔn)靜態(tài)拉伸真應(yīng)力-真應(yīng)變曲線Fig.2 True stress-true strain curves of smooth specimens during quasi-static tension

由圖3可以看出:隨著缺口半徑的增大,缺口試樣發(fā)生破壞時的載荷整體上呈現(xiàn)出減小的趨勢,并且發(fā)生破壞時的位移幾乎相同,遠(yuǎn)小于光滑圓棒試樣破壞時對應(yīng)的位移。這是因?yàn)閷τ贕H907合金這種塑性材料來說,缺口帶來的三向應(yīng)力狀態(tài)以及應(yīng)力集中現(xiàn)象,約束了材料內(nèi)部的塑性變形。缺口試樣的屈服強(qiáng)度和抗拉強(qiáng)度都比光滑試樣明顯提高,同時達(dá)到破壞時的變形量也更小,脆性有所增強(qiáng)。

圖3 光滑試樣和不同尺寸缺口試樣的準(zhǔn)靜態(tài)拉伸載荷-位移曲線(應(yīng)變速率0.001 s-1)Fig.3 Load-displacement curves of smooth and different-notch-sized specimens during quasi-static tension (0.001 s-1 strain rate)

由圖4可以看出,在變形溫度80~320 ℃范圍內(nèi),GH907合金準(zhǔn)靜態(tài)拉伸載荷-位移曲線彈性段的斜率基本一致,進(jìn)入強(qiáng)化段呈現(xiàn)出明顯的軟化趨勢,且隨著溫度的升高,合金達(dá)到峰值載荷所需位移有所增大。

圖4 不同溫度下光滑試樣的準(zhǔn)靜態(tài)拉伸載荷-位移曲線(應(yīng)變速率0.001 s-1)Fig.4 Load-displacement curves of smooth specimens during quasi-static tension at different temperatures (0.001 s-1 strain rate)

由圖5可以看出,在較高應(yīng)變速率(1 7002 700 s-1)的條件下,GH907合金沒有明顯的屈服段,同時屈服應(yīng)力和流動應(yīng)力都表現(xiàn)出明顯的應(yīng)變速率強(qiáng)化效應(yīng),高應(yīng)變速率下達(dá)到抗拉強(qiáng)度的應(yīng)變遠(yuǎn)小于準(zhǔn)靜態(tài)條件下的。

圖5 常溫下光滑試樣的準(zhǔn)靜態(tài)和動態(tài)拉伸真應(yīng)力-真應(yīng)變曲線Fig.5 True stress-true strain curves of smooth specimens during quasi-static tension and dynamic tension at room temperature

2.2 壓縮性能

由圖6可以看出,GH907合金具有明顯的應(yīng)變硬化效應(yīng),但是不同加載應(yīng)變速率下,硬化曲線幾乎平行,表明合金在壓縮應(yīng)變速率600~2 000 s-1范圍內(nèi)沒有明顯的應(yīng)變速率強(qiáng)化效應(yīng)。這歸因于在高的加載速率下,材料內(nèi)部的熱量來不及散發(fā)引起的熱軟化效應(yīng)與應(yīng)變速率強(qiáng)化效應(yīng)的相互抵消。

圖6 常溫下試樣的動態(tài)壓縮真應(yīng)力-真應(yīng)變曲線Fig.6 True stress-true strain curves of specimens during dynamic compression at room temperature

由圖7可以看出,隨溫度從室溫升到400 ℃,GH907合金硬化階段的整體真應(yīng)力-真應(yīng)變曲線有明顯的下移趨勢,表現(xiàn)出了材料的溫度軟化效應(yīng),同時卸載應(yīng)變也隨著溫度的升高呈現(xiàn)下降趨勢。

圖7 不同溫度下試樣的動態(tài)壓縮真應(yīng)力-真應(yīng)變曲線Fig.7 True stress-true strain curves of specimens during dynamic compression at different temperatures

3 本構(gòu)、失效模型的建立及驗(yàn)證

3.1 本構(gòu)模型的建立

J-C本構(gòu)模型[2]的表達(dá)式為

(1)

式(1)右側(cè)3項(xiàng)依次表征材料的硬化效應(yīng)、應(yīng)變速率效應(yīng)及溫度效應(yīng)。參數(shù)A,B,n可以通過室溫下光滑圓棒試樣的準(zhǔn)靜態(tài)拉伸試驗(yàn)獲得。在參考應(yīng)變速率和室溫下,式(1)的后2項(xiàng)均等于1,則式(1)可改寫為

(2)

參數(shù)A的值等于材料的初始屈服強(qiáng)度即塑性應(yīng)變?yōu)?時的應(yīng)力值,簡單計算可得A=666.2 MPa;參數(shù)B,n由參數(shù)A以及拉伸試驗(yàn)測得的應(yīng)力、應(yīng)變數(shù)據(jù)根據(jù)式(2)擬合獲得,并且考慮到引伸計在頸縮段劇烈變化時對應(yīng)變的測定并不準(zhǔn)確,參考文獻(xiàn)[14]的處理方法,僅采用頸縮前的數(shù)據(jù)進(jìn)行擬合,最終獲得B=2 365.4 MPa,n=0.974。

參數(shù)C可以通過常溫動態(tài)壓縮試驗(yàn)獲得。在常溫動態(tài)壓縮試驗(yàn)條件下,式(1)右側(cè)第3項(xiàng)等于1,則式(1)可以改寫為

(3)

按照式(3),繪制應(yīng)力項(xiàng)隨對數(shù)無量綱應(yīng)變速率的變化曲線,固定應(yīng)變的值即可獲得相應(yīng)的C值,C值為曲線的斜率。計算得到C=0.005。參數(shù)m的值可以通過高溫動態(tài)壓縮試驗(yàn)獲得。在固定的應(yīng)變速率下,式(1)可改寫為

(4)

按照式(4),繪制不同試驗(yàn)溫度下的對數(shù)應(yīng)力項(xiàng)隨對數(shù)無量綱溫度的變化曲線,固定應(yīng)變的值即可獲得相應(yīng)的m值,m值為曲線的斜率。計算得到m=2.0。

將由上述數(shù)據(jù)處理得到的J-C本構(gòu)模型參數(shù)A=666.2 MPa,B=2 365.4 MPa,n=0.974,C=0.005,m=2.0代入式(1)即建立了GH907合金的J-C本構(gòu)模型。利用該本構(gòu)模型擬合不同溫度和不同應(yīng)變速率下動態(tài)壓縮真應(yīng)力-真應(yīng)變曲線并與試驗(yàn)曲線進(jìn)行對比。由圖8可知,不同條件下的試驗(yàn)曲線與擬合曲線的塑性段基本一致,說明擬合所得的J-C本構(gòu)模型能夠較好地表征GH907合金的塑性流變行為。

圖8 J-C本構(gòu)模型擬合真應(yīng)力-真應(yīng)變曲線與試驗(yàn)曲線對比Fig.8 Comparison between fitting true stress-true strain curves ofJ-C constitutive model and test curves

3.2 失效模型的建立

J-C失效模型[3](也稱為失效準(zhǔn)則)在沖擊數(shù)值模擬中應(yīng)用廣泛,該模型將失效應(yīng)變表述為應(yīng)力三軸度、應(yīng)變速率和溫度的乘積關(guān)系,互不耦合。該模型表達(dá)式為

(5)

σ*=σm/σ

(6)

式中:D1,D2,D3,D4,D5為材料常數(shù);σ*為應(yīng)力三軸度;σm為平均應(yīng)力;σ為Mises等效應(yīng)力。

參數(shù)D1,D2,D3可以通過參考應(yīng)變速率(0.01 s-1)和室溫(20 ℃)下光滑圓棒試樣和缺口試樣的準(zhǔn)靜態(tài)拉伸試驗(yàn)獲得。在參考應(yīng)變速率和室溫下,式(5)的后2項(xiàng)均等于1,則式(5)可改寫為

εf=D1+D2exp(D3σ*)

(7)

根據(jù)應(yīng)力三軸度計算公式,即式(6),光滑圓棒試樣準(zhǔn)靜態(tài)拉伸試驗(yàn)時的應(yīng)力三軸度為1/3。對于缺口試樣,根據(jù)BRIDGMAN[16]的研究結(jié)果,初始應(yīng)力三軸度計算公式為

(8)

式中:d0和R0分別為試樣缺口截面的直徑和缺口處的半徑。

失效應(yīng)變可以采用TENG等[17]的定義,即

(9)

式中:A0和Af分別為試樣初始截面面積和斷口截面面積。

將拉伸試驗(yàn)數(shù)據(jù)代入式(8)、式(9),即可得到失效應(yīng)變與應(yīng)力三軸度數(shù)據(jù),通過式(7)擬合后,得到失效應(yīng)變隨應(yīng)力三軸度的變化曲線,如圖9所示,擬合參數(shù)D1=0.806,D2=-2.38×10-5,D3=6.334。

圖9 失效應(yīng)變與應(yīng)力三軸度的關(guān)系Fig.9 Relationship between fracture strain and stress triaxiality

在室溫、同一應(yīng)力三軸度下,式(5)右側(cè)的最后1項(xiàng)為1,則式(5)可改寫為

(10)

通過不同應(yīng)變速率下的拉伸試驗(yàn),利用式(10)擬合失效應(yīng)變與對數(shù)無量綱應(yīng)變速率,擬合結(jié)果見圖10,獲得參數(shù)D4=-0.017。

圖10 失效應(yīng)變與對數(shù)無量綱應(yīng)變速率的關(guān)系Fig.10 Relationship between fracture strain and logarithmic nondimensional strain rate

在參考應(yīng)變速率、同一應(yīng)力三軸度下,式(5)右側(cè)的中間項(xiàng)為1,則式(5)可改寫為

εf=[D1+D2exp(D3σ*)](1+D5T*)

(11)

利用式(11)擬合失效應(yīng)變與無量綱溫度的關(guān)系曲線,擬合結(jié)果見圖11,獲得參數(shù)D5=0.163。將由上述數(shù)據(jù)處理得到的參數(shù)D1=0.806,D2=-2.38×10-5,D3=6.334,D4=-0.017,D5=0.163代入式(5),即得到了GH907合金的J-C失效模型。

圖11 失效應(yīng)變與無量綱溫度的關(guān)系Fig.11 Relationship between fracture strain and nondimensional temperature

3.3 本構(gòu)模型的驗(yàn)證

通過上述分離變量法獲得模型的參數(shù)之后,通過與試驗(yàn)結(jié)果進(jìn)行比較來驗(yàn)證GH907合金本構(gòu)模型的有效性是十分必要的,其中有限元數(shù)值仿真是一種常用的方法。針對常、高溫SHPB壓縮試驗(yàn),在ABAQUS軟件中建立與之完全相同的三維幾何模型,試驗(yàn)用入射桿、透射桿長度均為1.5 m,直徑為14.5 mm;試樣尺寸為φ5.9 mm×3.9 mm。考慮到模型的對稱性,試驗(yàn)桿與試樣均只需建立1/4模型,如圖12所示,該有限元模型忽略了撞擊桿,選擇直接在入射桿的端面施加均布的應(yīng)力脈沖,脈沖的大小來自實(shí)際試驗(yàn)的采集。在對稱面施加對稱邊界條件;壓桿和試樣之間的接觸為硬接觸,光滑無摩擦;應(yīng)變值取自貼于入射桿和透射桿對應(yīng)位置的單元,與試樣接觸端面相距60 cm;對接觸區(qū)域網(wǎng)格進(jìn)行加密以獲得更加準(zhǔn)確的結(jié)果。霍普金森壓桿材料為60Si2MnA彈簧鋼,采用線彈性模型,其彈性模量為206 GPa,泊松比為0.29。GH907合金的材料參數(shù)采用前文擬合J-C本構(gòu)模型的數(shù)據(jù)。

圖12 SHPB壓縮試驗(yàn)有限元模型Fig.12 Finite element model of SHPB compression experiment

圖13給出了試驗(yàn)及仿真得到的3種條件下SHPB壓縮時應(yīng)變速率隨時間的變化曲線。由于試樣數(shù)量的限制,試驗(yàn)與仿真對比研究時未使用黃銅片加以整形來滿足常應(yīng)變速率加載。由圖13可以看出,在動態(tài)壓縮過程中,GH907合金應(yīng)變速率變化曲線的試驗(yàn)結(jié)果和仿真結(jié)果在應(yīng)變速率水平及持續(xù)時間上均有較好的一致性。由圖14可以看出:在動態(tài)壓縮過程中,GH907合金真應(yīng)力-真應(yīng)變曲線的試驗(yàn)結(jié)果與仿真結(jié)果非常接近,吻合較好。由表2可以看出,不同溫度、不同應(yīng)變速率下仿真得到的試樣幾何尺寸和最大應(yīng)力與試驗(yàn)結(jié)果相差很小,相對誤差均在2%以內(nèi),進(jìn)一步驗(yàn)證了擬合得到的J-C模型參數(shù)的有效性。

表2 試樣幾何尺寸和最大應(yīng)力仿真結(jié)果與試驗(yàn)結(jié)果的對比

圖13 動態(tài)壓縮過程中應(yīng)變速率變化曲線的仿真結(jié)果和試驗(yàn)結(jié)果對比Fig.13 Comparison between simulation and test results of strain rate change curves during dynamic compression

圖14 動態(tài)壓縮過程中真應(yīng)力-真應(yīng)變曲線的仿真結(jié)果和試驗(yàn)結(jié)果對比Fig.14 Comparison between simulation and test results of true stress-true strain curves during dynamic compression

4 結(jié) 論

(1) 常溫下在0~3 000 s-1應(yīng)變速率范圍內(nèi),拉伸時GH907高溫合金具有明顯的應(yīng)變速率強(qiáng)化效應(yīng),而壓縮時對應(yīng)變速率不敏感。

(2) GH907高溫合金在20~400 ℃溫度范圍內(nèi)、同一應(yīng)變速率下動態(tài)壓縮時產(chǎn)生軟化效應(yīng)。

(3) 基于不同的應(yīng)變速率和溫度下的拉伸和壓縮試驗(yàn)數(shù)據(jù),擬合得到GH907高溫合金的Johnson-Cook(J-C)本構(gòu)模型和失效模型;利用上述模型對動態(tài)壓縮過程進(jìn)行數(shù)值模擬,模擬得到的應(yīng)變速率變化曲線、真應(yīng)力-真應(yīng)變曲線與試驗(yàn)曲線吻合較好,試樣幾何尺寸和最大應(yīng)力與試驗(yàn)結(jié)果的相對誤差均在2%以內(nèi),驗(yàn)證了擬合所得本構(gòu)參數(shù)的準(zhǔn)確性。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 情侣午夜国产在线一区无码| 免费国产高清视频| 午夜无码一区二区三区在线app| 黄色片中文字幕| 中国毛片网| 国产亚洲精| 久操线在视频在线观看| 亚洲欧美色中文字幕| AV网站中文| 亚洲视频在线观看免费视频| 人妻出轨无码中文一区二区| 国产极品美女在线观看| 欧美色丁香| 亚洲Av综合日韩精品久久久| 日韩欧美高清视频| 国产精品九九视频| 丰满的熟女一区二区三区l| 国内精品久久九九国产精品 | 久一在线视频| 夜精品a一区二区三区| 97视频免费在线观看| 国产一级在线观看www色 | 永久免费精品视频| 免费人成黄页在线观看国产| 怡春院欧美一区二区三区免费| 精品国产三级在线观看| 亚洲精品色AV无码看| 国产爽歪歪免费视频在线观看| 亚洲Aⅴ无码专区在线观看q| 欧美激情,国产精品| 99草精品视频| 亚洲AV一二三区无码AV蜜桃| 色妞永久免费视频| 国产成人综合网| 国产激情国语对白普通话| 亚洲视屏在线观看| 国产精品久线在线观看| 亚洲成人播放| 中文字幕日韩视频欧美一区| 青青草国产一区二区三区| 欧美午夜在线视频| 久久精品无码国产一区二区三区| 国产日韩丝袜一二三区| 国产三级韩国三级理| 国产毛片基地| 久久久波多野结衣av一区二区| 亚洲黄色网站视频| 久久99国产精品成人欧美| 国产女人爽到高潮的免费视频 | 亚洲第一在线播放| 澳门av无码| 久久人体视频| 一级毛片不卡片免费观看| 亚洲精品无码在线播放网站| 日本精品视频一区二区| 精品国产香蕉在线播出| 国产一区二区人大臿蕉香蕉| 国产伦精品一区二区三区视频优播 | 最新无码专区超级碰碰碰| 国产亚洲日韩av在线| 久久精品中文字幕免费| 国产成人亚洲毛片| 欧美一道本| 国产大片喷水在线在线视频| 国产午夜福利亚洲第一| 99热亚洲精品6码| 国产精品综合色区在线观看| 欧美一区二区福利视频| 欧美精品黑人粗大| 欧美一级高清片欧美国产欧美| 亚洲欧美日韩天堂| 日韩欧美高清视频| 亚洲国产午夜精华无码福利| 99精品国产电影| 亚洲国产天堂在线观看| 欧美性猛交一区二区三区| 国产欧美视频综合二区| 福利在线一区| 亚洲人成人无码www| 国产91导航| 真实国产乱子伦高清| 国产乱子伦视频三区|