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

基于修正BW模型的戰(zhàn)斗部爆炸破片形成機理研究

2017-08-30 12:22:29吳衛(wèi)國杜志鵬張春輝
振動與沖擊 2017年15期
關(guān)鍵詞:模型

李 營, 吳衛(wèi)國, 杜志鵬, 張 瑋, 張春輝, 張 磊

(1.武漢理工大學(xué) 交通學(xué)院,武漢 430063;2.海軍裝備研究院,北京 100161)

基于修正BW模型的戰(zhàn)斗部爆炸破片形成機理研究

李 營1,2, 吳衛(wèi)國1, 杜志鵬2, 張 瑋2, 張春輝2, 張 磊1

(1.武漢理工大學(xué) 交通學(xué)院,武漢 430063;2.海軍裝備研究院,北京 100161)

膨脹環(huán)是研究戰(zhàn)斗部爆炸破片形成機理的重要手段。通過引入修正的B-W內(nèi)聚力斷裂模型ABAQUS子程序,開展了膨脹環(huán)碎裂過程的數(shù)值仿真計算,并將整個膨脹碎裂過程劃分為四個階段,并討論了膨脹速度對破片質(zhì)量分布和破片數(shù)量的影響,進一步研究了厚度對碎裂模式的影響。得到以下結(jié)論:①膨脹碎裂過程可以分為整體塑性、穩(wěn)定頸縮、局部頸縮發(fā)展和最終碎裂形成四個階段;②破片數(shù)量隨著初始膨脹速度增加而增大,且質(zhì)量分布服從Rayleigh分布;③徑厚比的變化改變了變形后期的應(yīng)力狀態(tài),薄殼結(jié)構(gòu)的應(yīng)力以受拉為主,而厚殼結(jié)構(gòu)后期出現(xiàn)剪切應(yīng)力狀態(tài)從而使得斷裂模式由拉伸徑縮失效逐漸過渡到剪切失效。

反艦導(dǎo)彈;膨脹環(huán);爆炸破片;Rayleigh分布;應(yīng)力三軸度

反艦導(dǎo)彈成為現(xiàn)代軍艦的重要威脅。其爆炸產(chǎn)生的高速破片和沖擊波是兩個主要的載荷,能有效毀傷艙室結(jié)構(gòu),并對艦用設(shè)備和艦員構(gòu)成重要威脅。研究反艦導(dǎo)彈戰(zhàn)斗部爆炸破片形成機理對于準(zhǔn)確掌握破片特性,提高艦船防護水平具有重要意義。膨脹環(huán)是研究材料碎裂特性和戰(zhàn)斗部爆炸破片形成機理的有效手段[1]。采用數(shù)值方法輔助開展膨脹環(huán)碎裂行為的研究是目前常用的技術(shù)手段。數(shù)值仿真水平的提高主要通過2個途徑實現(xiàn),分別為建立合理的本構(gòu)關(guān)系更為準(zhǔn)確地表征材料的力學(xué)行為,另外一種為提高數(shù)值計算方法的速度和精度。有關(guān)研究表明[2-3],本構(gòu)關(guān)系表征不準(zhǔn)帶來的誤差遠遠大于計算方法帶來的誤差。

斷裂準(zhǔn)則研究是材料斷裂行為研究的基礎(chǔ),已經(jīng)經(jīng)歷了若干代發(fā)展。早期一直采用等效塑性應(yīng)變方法。Gurson等[4]通過微觀空穴形成和擴展的機理分析提出斷裂應(yīng)變與應(yīng)力三軸度相關(guān)。Johnon等[5]提出目前廣泛應(yīng)用的Johnson-cook斷裂準(zhǔn)則,綜合考慮了應(yīng)力三軸度效應(yīng)、應(yīng)變率效應(yīng)和溫度效應(yīng),得到學(xué)術(shù)界的廣泛關(guān)注,并集成到了DYNA、ABAQUS等商業(yè)軟件中。Bao等[6]對低應(yīng)力三軸度區(qū)間進行了重大改進,提出了B-W準(zhǔn)則,能較為明確地區(qū)分延性損傷、剪切失效模式,受到國外學(xué)者的廣泛重視[7-8]。

目前,已經(jīng)有一些學(xué)者采用數(shù)值仿真方法開展膨脹環(huán)碎裂特性的研究,但由于斷裂準(zhǔn)則的選取相對簡單,再現(xiàn)復(fù)雜碎裂機理和碎裂模式方面仍有待改進。Rusinek等[9]利用ABAQUS對膨脹環(huán)碎裂過程進行了仿真,但假設(shè)應(yīng)力應(yīng)變曲線一旦出現(xiàn)平臺就發(fā)生斷裂,與實際情況不符。Zhang等[10]認為等效塑性應(yīng)變?yōu)?時材料破壞,并采用單元刪除方法表征材料失效,此種方法未能反應(yīng)材料破壞能的破壞機理。陳磊等[12-13]考慮了材料斷裂能,并基于Johnson-cook斷裂準(zhǔn)則利用ABAQUS對韌性膨脹環(huán)的膨脹碎裂進行了,得到了破片的平均尺寸,對Grady-Kipp準(zhǔn)則進行了驗證,但未考慮低應(yīng)力三軸度區(qū)間材料破壞機理的變化對碎裂特性的影響,使得破壞模式的預(yù)測不夠準(zhǔn)確。尚未見到使用B-W準(zhǔn)則研究膨脹環(huán)碎裂特性的研究,該準(zhǔn)則對膨脹碎裂模式影響的研究有待開展。

本文將修正的B-W內(nèi)聚力斷裂模型通過材料子程序方法集成進ABAQUS, 并以此為基礎(chǔ)開展了膨脹碎裂的數(shù)值仿真研究,討論了膨脹速度對頸縮演化、破片分布的影響,研究了膨脹環(huán)徑厚對碎裂模式的影響。

1 爆炸破片形成基本理論

Mott早在1947年就給出了膨脹環(huán)破片的基本理論,示意圖見圖1。圓環(huán)在沖擊瞬間膨脹半徑為r、速度為V0,則此時的應(yīng)變率為

(1)

圖1 戰(zhàn)斗部碎裂示意圖

假設(shè)有一初始斷口,則該處形成自由界面,并形成卸載波。卸載波向兩側(cè)傳播,并形成一段無應(yīng)力區(qū)。該卸載波的波速為

(2)

Mott給出平均破片長度的預(yù)測公式,如下

x0=(2PF/ργ)1/2r/V0

(3)

式中:PF為相應(yīng)應(yīng)變率下的流動斷裂應(yīng)力;ρ為材料密度;r為膨脹環(huán)半徑;v為膨脹速度。

2 低碳鋼碎裂破壞的力學(xué)特性

2.1 塑性流動準(zhǔn)則

Johnson-cook強度模型[14]是一種黏熱塑性本構(gòu)關(guān)系,這種模型能較好地描述金屬材料的加工硬化效應(yīng)、應(yīng)變率效應(yīng)和溫度軟化效應(yīng)。Johnson-cook強度模型的形式為

(4)

塑性溫度升高為

(5)

式中:ρ,c分別為材料的密度和比熱;β為Taylor-quinney常數(shù),取為0.9。

2.2 內(nèi)聚力損傷準(zhǔn)則

BW模型在JC模型的基礎(chǔ)上對應(yīng)力三軸度區(qū)間進行了修正,將應(yīng)力三軸度對斷裂應(yīng)變的影響進行分段處理。具體公式如下:

(6)

式中:D01、D02、D03為材料參數(shù),通過不同材料試樣獲得。

作者前期工作中,開展了不同溫度、不同應(yīng)變率、多種應(yīng)力狀態(tài)下的材料斷裂試驗,對應(yīng)力率效應(yīng)和溫度效應(yīng)進行修正,得到了采用修正的B-W內(nèi)聚力失穩(wěn)斷裂準(zhǔn)則,具體參數(shù)參見文獻[15]。用等效塑性應(yīng)變作為損傷因子D的開動準(zhǔn)則。損傷開動的臨界應(yīng)變εd是一個與應(yīng)力三軸度、應(yīng)變率及溫度相關(guān)的量,可以用下式表示

(7)

式中,D01、D02、D03與式(6)相同,D1、D2、D3、D4、D5也為材料的斷裂參數(shù),通過一系列材料試驗獲得。

一旦損傷開動,假定損傷發(fā)展隨單元內(nèi)部塑性變形的發(fā)展而線性增加,即

(8)

式中:σ0為損傷開動時刻的單元應(yīng)力,Gc是損傷開動到材料完全失效(D=1)時所需要的斷裂能量;up為損傷開動后單元的塑性位移,由下式給出

up=Lmesh(Δεp)

(9)

由此可見,若發(fā)生相同的等效塑性應(yīng)變,大單元尺寸的損傷將超過小尺寸的損傷。事實上,當(dāng)材料發(fā)生一裂紋演化為主的破壞時,局部塑性應(yīng)變往往集中在小單元內(nèi),公式描述的損傷發(fā)展模型將材料的抗損傷發(fā)展通過材料的斷裂能表征,既可以描述材料的分離過程,又可以減少數(shù)值計算對網(wǎng)格尺寸的依賴。

(10)

一旦損傷達到1,材料失去承載力,利用單元消除技術(shù)將失效單元消除。

2.3 材料模型驗證

采用ABAQUS進行數(shù)值仿真計算,通過實驗和數(shù)值仿真的比對可以看出:改進的BW模型較好再現(xiàn)了Taylor桿高速撞擊后較為復(fù)雜的花瓣狀破壞模式[16]。說明改進的BW模型能較好預(yù)測低碳鋼在高應(yīng)變率、復(fù)雜應(yīng)力狀態(tài)下的損傷失效行為,為進一步開展膨脹環(huán)碎裂計算奠定了基礎(chǔ)。

定義應(yīng)力張量的6個分量,更新各時間步中的應(yīng)力應(yīng)變參量,通過FORTRAN編譯器編寫VUMAT材料程序代碼,定義應(yīng)力三軸度、應(yīng)變率和溫度對材料動態(tài)損傷特性的影響。在ABAQUS開展計算時,調(diào)用子程序開展動態(tài)過程計算,并最終得到材料動態(tài)碎裂過程及終點效應(yīng)。

A/MPaB/MPanCmGc/kND01249.28890.7460.0580.9425-6.743D02D03D1D2D3D4D50.0451.3250.2961.184-1.4650.0058.07

3 數(shù)值計算

3.1 數(shù)值模型

利用ABAQUS建立了膨脹環(huán)的三維仿真模型,圓環(huán)內(nèi)徑r=20 mm、外徑R=21 mm,高度h=1 mm,如圖3所示。共劃分了10 000個C3D8R單元,其中膨脹環(huán)厚度方向上劃分5個單元。ABAQUS中總體坐標(biāo)系為笛卡爾坐標(biāo),無法直接加載徑向速度,通過Matlab編程實現(xiàn)初始膨脹速度V0。材料參數(shù)設(shè)置見表1。

圖3 計算模型示意圖

3.2 膨脹基本過程

以初速度100 m/s的膨脹環(huán)的膨脹過程為例,整個膨脹過程分為四個典型階段,如圖3所示。(Ⅰ)整體塑性階段(0~18.03 μs):隨著膨脹環(huán)膨脹,整體均勻變形,并發(fā)生產(chǎn)生較大塑性,塑性溫升可達到204 K。(Ⅱ)穩(wěn)定頸縮階段(18.03~25.2 μs):局部區(qū)域塑性變形急劇增加,并形成一系列頸縮;非頸縮區(qū)域的塑性變形停止,溫度不再升高(圖5中B點)。(Ⅲ)局部頸縮發(fā)展階段(25.2~29.4 μs):隨著時間的推移,部分頸縮區(qū)域快速發(fā)展,溫度升高進一步加快(如圖5中C點);一部分頸縮則由于相鄰快速發(fā)展頸縮區(qū)形成的卸載波,使得塑性變形受到抑制(Arrested necking),并最終成為破片中的殘余頸縮(圖5中A點);(Ⅳ)最終碎裂階段(29.4 μs以后),局部斷口形成,碎裂形成階段。

圖4 破片形成的四個階段

圖5為不同速度下頸縮發(fā)展變化的示意圖。當(dāng)膨脹速度為400 m/s時,單位無量綱周長(該點沿周長方向的距離/膨脹環(huán)半徑) 上的頸縮為7個,其中有2個頸縮位置后期持續(xù)發(fā)展,并最終形成斷口。當(dāng)局部頸縮發(fā)展越快,對周圍頸縮區(qū)域的抑制作用越強。

圖5 局部區(qū)域的塑性溫升

當(dāng)膨脹速度達到100 m/s時,單位無量綱周長上的頸縮為13個,其中5個在后期形成斷口,斷口形成的自由界面產(chǎn)生mott卸載波,其余頸縮受其作用,應(yīng)力卸載,塑性變形不再增加。

(a) V0=400 m/s

(b) V0=1 000 m/s

3.3 初始膨脹速度對碎裂特性的影響

隨著初始膨脹速度的增加,根據(jù)公式(1)膨脹環(huán)的應(yīng)變率呈線性增加,材料臨界斷裂應(yīng)變發(fā)生增大(公式(6))。同時,從公式(2)可以看出,Mott卸載波的傳播速度降低,均可對最終碎裂形態(tài)有影響。

圖7為初始膨脹速度為400 m/s、600 m/s、800 m/s和1 000 m/s的膨脹環(huán)碎裂最終形態(tài)。可以看出隨著初始膨脹速度的增大,破片數(shù)量逐漸增多,碎片體積逐漸減小。表2可以看出,相同膨脹速度的膨脹環(huán),最終破片數(shù)量接近但不完全相同,每個速度計算5組,統(tǒng)計得到該速度下的破片總數(shù),分別為26、106、147和180個。擬合一次函數(shù)直線如圖8所示。

圖7 不同速度的破片形態(tài)

圖8 破片數(shù)量與初速度的關(guān)系

初速度環(huán)1破片數(shù)環(huán)2破片數(shù)環(huán)3破片數(shù)環(huán)4破片數(shù)環(huán)5破片數(shù)破片樣本數(shù)40064565266002221212022106800263129293214710003535333839180

圖9為統(tǒng)計不同初速度膨脹環(huán)最終破片的質(zhì)量分布特性。采用n=2的Weibull分布,即Rayleigh分布可以可以較好地表示其破片分布特性[12],即

相應(yīng)地,Rayleigh概率密度函數(shù)為

(a) 破片質(zhì)量/g

(b) 破片質(zhì)量/g

(c) 破片質(zhì)量/g

(d) 破片質(zhì)量/g

3.4 厚度對碎裂機理的影響

有關(guān)研究表明,薄殼戰(zhàn)斗部的破壞模式主要以拉伸失效為主,但也有研究實驗中出現(xiàn)剪切為主的失效模式。

開展了徑厚比α=R/(R-r)分別為0.05、0.33和0.4,初始膨脹速度均為500 m/s的膨脹碎裂過程模擬。如圖10,三種工況分別代表三種典型破壞模式:① 延性頸縮破壞,該破壞模式下,斷口頸縮,主要發(fā)生延性破壞;② 復(fù)合損傷,材料發(fā)生延性擴孔破壞和剪切損傷破均存在的復(fù)合損傷;③ 剪切損傷,主要機理為形成局部剪切帶或絕熱剪切帶,剪切帶與徑向方向呈一定角度。

如圖11所示,隨著α的不同,膨脹環(huán)在變形過程中應(yīng)力狀態(tài)發(fā)生較大改變。在變形初期,膨脹環(huán)與單軸拉伸應(yīng)力狀態(tài)類似,應(yīng)力三軸度一直維持在1/3。變形后期,徑厚比較小(α=0.05)的環(huán)應(yīng)力三軸度逐漸超過1/3,呈現(xiàn)多軸拉伸趨勢;徑厚比較大的環(huán)應(yīng)力三軸度逐漸開始下降,并低于0。

應(yīng)力三軸度與破壞模式的關(guān)系如圖12所示。已有研究表[1],σ*<-1/3時,材料受純壓作用,材料不發(fā)生損傷;-1/3<σ*<0時,材料受力從純壓向純剪過度,為壓剪聯(lián)合作用階段,主要破壞機理為剪切破壞;0<σ*<1/3時,材料受力從純剪切向拉伸過度,為拉剪聯(lián)合作用階段,破壞模式為延性擴孔與剪切混合模式;σ*>1/3時,損傷機理為延性損傷。在區(qū)間(1)中,已經(jīng)從實驗和理論上證明,金屬材料不會發(fā)生損傷和破壞;在區(qū)間(2)、(3)中,金屬材料主要發(fā)生剪切型損傷破壞;在區(qū)間(4)中,金屬材料主要發(fā)生延性損傷破壞。

圖12 應(yīng)力三軸度對損傷模式轉(zhuǎn)變的影響

上述分析可以看出,徑厚比的不同主要改變了膨脹環(huán)變形后期的應(yīng)力狀態(tài),進一步改變了微觀損傷斷裂機理。在中等厚度的戰(zhàn)斗部中會出現(xiàn)拉伸和剪切混合破壞模式。

4 結(jié) 論

本文給出了基于修正B-W模型的材料碎裂本構(gòu)關(guān)系,研究了戰(zhàn)斗部爆炸破片形成機理,討論了破壞模式變化規(guī)律。主要得到以下幾點結(jié)論:

(1) 基于JC模型的流動應(yīng)力準(zhǔn)則和改進B-W模型的內(nèi)聚力斷裂準(zhǔn)則可以有效模擬膨脹環(huán)碎裂過程。

(2) 膨脹環(huán)膨脹碎裂過程可以分為整體塑性階段、穩(wěn)定頸縮階段、局部頸縮發(fā)展階段和最終碎裂形成階段;局部頸縮發(fā)展形成的Mott卸載波對鄰近區(qū)域的應(yīng)力卸載作用明顯,抑制了鄰近頸縮的進一步發(fā)展。

(3) 隨著初始膨脹速度的增大,破片數(shù)量呈線性增加,且破片質(zhì)量分布服從Rayleigh分布。

(4) 徑厚比能顯著影響碎裂機理,隨著徑厚比的增大,碎裂模式逐漸從延性頸縮失穩(wěn)破壞,過渡到延性破壞與剪切破壞交替的混合破壞模式,并最終過渡到剪切帶破壞;破壞模式改變的主要原因是膨脹環(huán)變形后期應(yīng)力狀態(tài)的差異。

[1] MOTT N F. Fragmentation of shell cases[J]. Proceedings of the Royal Society of London Series A Mathematical and Physical Sciences, 1947, 189: 300-308.

[2] SCHEFFLER D R, ZUKA J A. Practical aspects of numerical simulation of dynamic events: material interfaces[J]. International Journal of Impact Engineering, 2000, 24(8): 821-842.

[3] ZUKAS J A, SCHEFFLER D R. Practical aspects of numerical simulations of dynamic events: effects of meshing[J]. International Journal of Impact Engineering, 2000, 24(9): 925-945.

[4] GURSON A L. Plastic flow and fracture behavior of ductile meterials incorporating void nucleation,growth and interaction[D]. Brown University, 1975.

[5] JOHNSON G R, COOK W H. Fracture Characteristics of three metals subjected to various strains, strain rates, temperatures and pressures[J]. Engineering Fracture Mechanics, 1985, 21: 31-48.

[6] BAO Y, WIERZBICKI T. On fracture locus in the equivalent strain and stress triaxiality space[J]. International Journal of Mechanical Sciences, 2004, 46(1): 81-98.

[7] BARSOUM I, FALESKOG J. Rupture mechanisms in combined tension and shear—Experiments[J]. International Journal of Solids and Structures, 2007, 44: 1768-1786.

[8] BARSOUM I, FALESKOG J F. Rupture mechanisms in combined tension and shear-Micromechanics[J]. International Journal of Solids and Structures, 2007, 44: 5481-5498.

[9] RUSINEK A, ZAERA R. Finite element simulation of steel ring fragmentation under radial expansion[J]. International Journal of Impact Engineering, 2007, 34(4): 799-822.

[10] ZHANG H, RAVI-CHANDAR K. On the dynamics of necking and fragmentation—II.Effect of material properties, geometrical constraints and absolute size[J]. Int J Fract, 2008, 150: 3-36.

[11] 陳磊,周風(fēng)華,湯鐵鋼. 韌性金屬圓環(huán)高速膨脹過程的有限元模擬[J]. 力學(xué)學(xué)報, 2011, 43(5): 861-870.

CHEN Lei, ZHOU Fenghua, TANG Tiegang. Finite element simulations of the high velocity expansion and fragmentation of ductile metallic rings[J]. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(5): 861-870.

[12] 鄭宇軒. 韌性材料的動態(tài)碎裂特性研究[D]. 合肥:中國科技大學(xué), 2013.

[13] 鄭宇軒,陳磊,胡時勝,等. 韌性材料沖擊拉伸碎裂中的碎片尺寸分布規(guī)律[J]. 力學(xué)學(xué)報, 2013, 45(4): 580-587.

ZHENG Yuxuan, CHEN Lei, HU Shisheng, et al. Characteristics of fragment size distribution of ductile materials fragmentized under high strain rate tension[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(4): 580-587.

[14] JOHNSON G R, COOK W H. A constitutive model and data for metals subjected to large strains, high strain rates and high temperature[C]. Proceedings of the seventh international symposium on ballistics, Netherland, 1983.

[15] 李營. 液艙防爆炸破片侵徹作用機理研究[D]. 武漢:武漢理工大學(xué), 2014.

[16] 郭子濤. 彈體入水特性及不同介質(zhì)中金屬靶的抗侵徹性能研究[D]. 哈爾濱:哈爾濱工業(yè)大學(xué), 2012.

Blast fragments forming mechanism of a warhead based on modified BW model

LI Ying1,2, WU Weiguo1, DU Zhipeng2, ZHANG Wei2, ZHANG Chunhui2, ZHANG Lei1

(1. School of Transportation, Wuhan University of Technology, Wuhan, 430063, China;2. Naval Academy of Armament, Beijing 100161, China)

Expanding ring is a useful means to study a warhead’s blast fragments forming mechanism. Here, through introducing a modified BW cohesive force fracture model’s ABAQUS subroutine, the fragmentation process of an expanding ring was simulated numerically. The expanding fragmentation process was divided into four stages. The effects of expanding velocity on fragments’ mass distribution and fragments number were discussed. The effects of thickness on the fragmentation mode were also studied. It was shown that the fragmentation process is divided into 4 stages including overall plastic deformation, steady necking, local necking developing and fragmentation; the fragments number increases with increase in the initial expanding velocity and the mass distribution of fragments obeys Rayleigh distribution; the variation of the radius-thickness ratio changes stress status of end deformation stage, the stress of thin shell structures is mainly tensile, while the shear stress status appears in the end stage of thick shell structures, so their fracture modes transfer gradually from tension necking failure to shear failure.

anti-ship missile; expanding ring; blast fragmentation; Rayleigh distribution; stress tri-axiality

國防基礎(chǔ)研究項目(B1420133057);國家自然科學(xué)基金(51509196);非線性力學(xué)國家重點實驗室開放基金(LNM201505);中央高校基本科研業(yè)務(wù)費專項資金(2014-YB-20)。

2016-03-10 修改稿收到日期:2016-06-14

李營 男,博士生,1988年生

O389;TJ410

A

10.13465/j.cnki.jvs.2017.15.018

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 成年A级毛片| 亚洲一区二区在线无码 | 99在线视频精品| 国产精品入口麻豆| 91精品啪在线观看国产60岁 | 黄色网在线| 国产杨幂丝袜av在线播放| 91在线一9|永久视频在线| 亚洲伊人天堂| 欧美成人免费一区在线播放| 亚洲国产成人精品一二区| 99视频在线免费| 欧美激情福利| 国产成人做受免费视频| 免费在线成人网| 91成人免费观看| 伊人91在线| 日本成人福利视频| 欧美成人看片一区二区三区| 免费在线看黄网址| 欧美五月婷婷| 在线免费无码视频| 久久免费看片| 99在线观看视频免费| 亚洲欧美极品| 国产一级在线观看www色 | 国产精品专区第1页| 亚洲精品在线91| 久久人人97超碰人人澡爱香蕉 | 欧美日韩国产在线播放| 91区国产福利在线观看午夜 | 广东一级毛片| 综合色88| 国产情侣一区二区三区| 亚洲天堂在线免费| 日韩中文无码av超清| 国产成人8x视频一区二区| 在线精品视频成人网| 97超爽成人免费视频在线播放| 六月婷婷激情综合| 欧美午夜视频在线| 老色鬼久久亚洲AV综合| 伊人蕉久影院| 999福利激情视频| 精品成人免费自拍视频| 欧美精品影院| 51国产偷自视频区视频手机观看 | 亚洲国产成人久久77| 免费一级大毛片a一观看不卡| 国产一区二区网站| 色欲国产一区二区日韩欧美| 亚洲欧洲一区二区三区| 亚洲精品不卡午夜精品| 国产亚洲视频中文字幕视频| 国内自拍久第一页| 国产香蕉一区二区在线网站| 综合色88| 精品人妻无码区在线视频| 最新日本中文字幕| 亚洲综合极品香蕉久久网| 日本三区视频| 在线中文字幕网| 欧美综合区自拍亚洲综合天堂| 亚洲看片网| 欧美不卡视频一区发布| 日本成人在线不卡视频| 欧美无专区| 国产打屁股免费区网站| 国产乱子伦手机在线| 国产在线精品人成导航| 国产正在播放| 日韩在线成年视频人网站观看| 91成人精品视频| 成人夜夜嗨| 亚洲视频一区| 亚洲全网成人资源在线观看| 国产成人综合日韩精品无码不卡| 一本色道久久88亚洲综合| 91精品日韩人妻无码久久| 国产又粗又猛又爽视频| 日本一区二区不卡视频| 伊人久久大香线蕉成人综合网|