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

基于SPH-FEM方法的半球形聚能裝藥破甲特性研究

2016-08-04 07:00:07張之凡王龍侃張阿漫
振動與沖擊 2016年14期

張之凡, 李 兵, 王龍侃, 張阿漫

(1.哈爾濱工程大學 船舶工程學院,哈爾濱 150001; 2.中國人民解放軍91439部隊,遼寧 大連 116041)

?

基于SPH-FEM方法的半球形聚能裝藥破甲特性研究

張之凡1, 李兵2, 王龍侃1, 張阿漫1

(1.哈爾濱工程大學 船舶工程學院,哈爾濱150001; 2.中國人民解放軍91439部隊,遼寧 大連116041)

摘要:為了研究桿式射流的形成及破甲過程,基于SPH-FEM方法建立半球形聚能裝藥模型,對半球形聚能裝藥起爆后形成金屬射流及射流穿透雙殼的過程進行仿真模擬。通過對射流形成過程、速度衰減規律以及殼體破口形成過程進行分析,得到以下結論:在半球形聚能裝藥的爆轟作用下,藥型罩被高壓沖擊變形繼而形成了高速的金屬射流;射流頭部最大速度可達約4 174 m/s;射流在擊穿雙殼的過程中會發生斷裂形成射流斷裂塊,第一層殼在被擊穿過程中,經歷了沖塞、凹陷等過程,第二層殼直接被射流穿透,之后不斷地被射流斷裂塊擊穿。整個計算旨在對半球形聚能炸藥的工程設計提供參考。

關鍵詞:半球形聚能裝藥;SPH-FEM方法;射流;破甲

聚能裝藥[1-5]從起爆到形成高速射流,繼而進行破甲的過程是一個極其復雜的物理過程,爆轟波從藥型罩頂部掠至尾部時,將藥型罩以很大的速度向軸向擠壓,此時藥型罩在極大的壓力作用下發生塑性流動加速,最后匯成一股具有高速及高能量的金屬射流[6-7],隨后金屬射流擊穿殼板,在殼板中形成高溫、高壓、高應變率區域,實現侵徹破甲。半球形聚能裝藥能夠形成桿式射流,射流頭部速度可達3~5 km/s,不同于其它射流,桿式射流對炸高的敏感程度較低,但是藥型罩利用率較高,破甲能力較高,可以用于攻擊新型防護裝甲、武裝直升機和大型水面艦艇等目標,所以本文對半球形聚能裝藥形成桿式射流和其擊穿雙殼結構過程進行研究。隨著聚能裝藥被廣泛應用于軍事以及工程領域,國內外很多專家學者都對聚能射流的形成以及射流對結構的毀傷進行了深入的研究,Molinari[8]采用有限元方法(Finite Eliment Method,FEM)對聚能射流的形成、斷裂及穿透靶板過程進行數值模擬;Liu 等[9-10]基于SPH (Smoothed Particle Hydrodynamics) 方法對線性聚能裝藥射流的形成以及穿透效果進行研究;溫萬治等[11]基于MOCL (Markon Cell Line) 分界面跟蹤算法,采用二維多流體網格法,模擬了錐形罩聚能裝藥侵徹鋼板的全過程;張先鋒等[12]利用有限元軟件對三種典型聚能射流的成型及侵徹過程進行了數值模擬;李磊等[13]采用SPH算法實現了錐形罩聚能裝藥射流形成過程的三維數值模擬。然而,以上數值模擬的都是錐形聚能裝藥的射流形成過程,半球形聚能裝藥形成的射流寬度較大,對結構的毀傷更大,破甲效果更好,并且SPH方法[6, 14-17]的無網格特性能克服結構大變形引起的網格畸變等困難,這種特性使得它非常適宜對聚能射流形成過程進行模擬分析。因此,本文結合SPH方法和FEM 方法各自的優越性,應用SPH-FEM方法將半球形聚能裝藥的射流形成過程以及射流擊穿雙殼結構的過程進行數值仿真,首先采用SPH方法建立半球形聚能裝藥模型,模擬金屬射流的形成過程,并對射流速度進行分析。隨后,在此基礎上,建立SPH-FEM聚能射流擊穿雙殼結構模型,進一步研究金屬射流對雙殼結構的毀傷過程,分析射流頭部速度的衰減規律以及兩層殼體破口的形成,旨在為相關的工程應用提供參考。

1基本理論

1.1SPH基本方程

考慮材料的本構關系以及人工黏度Πij,得到具有材料強度的流體動力學控制方程的SPH形式如下[6]:

式中:ρ為密度;v為速度;a、b表示坐標方向;p為壓力;e為內能;t為時間;x為粒子的位移;Wij為粒子j對粒子i產生影響的光滑函數,本文中所應用的光滑函數為分段三次樣條光滑函數[6];ε為應變率;τ為剪切應力;Πij為人工黏度[6]。

1.2物態方程

(1)爆炸產生的爆轟產物采用Jones-Wilkins-Lee狀態方程[18]:

(2)

式中:η為爆炸氣體的密度和原始爆炸物的初始密度的比值;ρ0為初始炸藥的密度;P為爆轟產物的壓力;e0為高能炸藥單位質量的內能;A、B、R1、R2為實驗數據擬合所得的系數,其具體取值見表1。

表1 TNT爆轟產物Jones-Wilkins-Lee狀態方程參數

(2) 藥型罩和鋼的的狀態方程采用固體力學的Mie-Gruneisen狀態方程[19-20]

(3)

式中:Γ為Gruneisen常數,pH為沖擊Hugoniot曲線上的點的壓力,其表達式為[6]:

pH(ρ)=

(4)

表2 金屬Mie-Gruneisen狀態方程參數

(3) 藥型罩和鋼的屈服強度采用Johnson-Cook屈服模型[21]:

(5)

表3 金屬材料本構模型參數

2聚能射流形成的數值模擬

2.1計算模型

半球形聚能裝藥結構如圖1所示,其中裝藥采用TNT炸藥,炸藥高l=0.21 m,半徑R=0.08 m,裝藥中心點起爆;藥型罩材料為紫銅,厚度為d=0.004 m。采用不均勻分布的粒子對該模型進行模擬,炸藥的粒子間距為dx1=0.001 m,藥型罩的粒子間距為dx2=0.002 5 m。

圖1 聚能裝藥計算模型Fig.1 Model of shaped charge

2.2模擬結果及分析

圖2為半球形聚能裝藥射流形成過程及不同時刻的速度云圖,初始時刻半球形藥型罩如圖中t=0時刻所示;t≈48 μs時,在爆轟壓力的作用下藥型罩的頂部被壓垮,并且頂部開始出現外翻,此時頭部速度達到約3 281 m/s;t≈54 μs時藥型罩繼續翻轉變形,由于藥型罩還受到軸向拉伸和徑向壓縮,從圖中可以看出,整個藥形罩向軸線處匯聚并出現拉伸,形成了初始射流,此時頭部速度約為3 822 m/s;隨著射流的形成,t≈66 μs時藥型罩的原內表面漸漸變成外表面,原外表面變成內表面,從圖中可以看出,藥型罩頂部已經完全被翻轉到外面成為射流的頭部,此時頭部速度達到峰值,約為4 174 m/s,與實驗得到的X光照片提供的數據(約4 000 m/s)相吻合[22];隨著金屬射流繼續向前運動,射流金屬拉長變細,頭部速度開始減小,t≈84 μs時頭部速度約為3 876 m/s,從圖中還可以看出射流尾部出現杵;t≈102 μs時尾部杵的長度不斷增大,此時射流頭部速度繼續降低,約為3 821 m/s。由此可得,藥型罩在爆轟壓力的作用下發生高度變形,同時由于在變形過程中還受到軸向拉伸和徑向壓縮,藥型罩表面激烈變形和碰撞,進而在軸線上聚合形成能量密度更高的聚能射流,且射流頭部的速度較高,尾部杵的速度較低,基本上遵循線性的速度分布,符合射流速度分布的基本規律[22]。

圖2 半球形聚能裝藥射流形成過程(速度云圖)Fig.2 Shaped charge jet formation (Velocity contour)

圖3為射流頭部和尾部的速度時間曲線,圖中空心方塊為射流尾部的測點A,實心圓圈為射流頭部的測點B。從圖中可以看出,射流頭部和尾部的速度變化趨勢相近,都為先增大后減小,整個變化趨勢與上圖中的速度云圖相對應。對于射流尾部,射流末端的粒子具有最大速度值,由測點A的速度曲線可以得到,t≈42 μs時射流尾部的速度達到峰值,約為2 010 m/s,隨后速度不斷減小,最終趨于穩定。對于射流頭部,t≈66 μs時速度達到最大值,約為4 174 m/s,隨后速度減小,t≈108 μs時速度趨于平緩,最終頭部速度約3 700 m/s。由此可得,射流頭部及尾部的速度先增大后減小,頭部最大速度可以達到約4 174 m/s,尾部最大速度可以達到約2 010 m/s,與錐形罩相比,射流頭部速度較低,但是藥型罩利用率較高,射流直徑較粗,有利于進一步破甲。

圖3 射流頭部和尾部的速度時間曲線(測點A和B如左圖中標注:A點為射流尾部的測點,B為射流頭部的測點)Fig.3 Velocity-time curves of jet head and jet tail (As left figure shows, A and B are test points of jet head and jet tail)

3聚能射流擊穿雙殼的數值模擬

3.1計算模型

由于結構的對稱性,所以本文采用1/2模型進行計算,模型中所用的狀態方程及本構關系如1.2節所示,聚能裝藥采用2.1節中的SPH半球形聚能裝藥計算模型,兩層殼的長和寬皆為a=b=0.5 m,網格單元為四邊形單元,網格單元大小為0.002 5 m,對中心區域進行局部畫細,局部加密區域長為a1=0.15 m,寬為b1=0.15 m,網格單元大小為0.001 25 m。炸藥和藥型罩為光滑粒子,殼板為 Lagrange 單元網格,光滑粒子與 Lagrange 單元網格之間采用點到面侵蝕接觸算法。炸高h=0.15 m,第一層殼厚為6 mm,第二層殼厚為20 mm,兩層殼之間的垂直距離為h′=0.25 m,整個計算模型如圖4所示。

圖4 聚能裝藥作用對結構毀傷的計算模型Fig.4 Model of double-shells subjected to shaped charge

3.2模擬結果及分析

3.2.1射流擊穿雙層殼的過程

圖5為聚能射流的形成及穿透雙層殼的過程,聚能裝藥起爆后爆轟產物會推動半球形藥型罩,形成射流彈丸,t≈94 μs時,射流到達第一層殼,隨后射流完全擊穿第一層殼,由于射流速度很大,射流頭部形狀沒有太大變化,依然保持細長的形狀;t≈114 μs時,射流開始斷裂,頭部出現小塊射流;t≈136 μs時射流繼續斷裂,此時已經斷裂成了四塊,從圖中可以發現,頭部斷裂的小塊射流也在不斷被拉長;t≈168 μs,射流開始穿透第二層殼,隨后斷裂的射流陸續擊穿第二層殼,由于第二層殼的阻力作用,射流開始不斷的堆積,斷裂的射流漸漸融合;t≈240 μs時,射流的主體部分已經穿透第二層殼,由于射流斷裂塊的堆積,射流主體在擊穿第二層殼時發生了“反射”現象,有小部分粒子由于第二層殼的阻力作用而向射流的反方向飛散,從圖中還可以看出,已經穿透第二層殼的頭部射流金屬已經發生飛散,而射流的主要部分由于殼體的阻力作用發生堆積,導致殼體的破口不斷增大。整個過程中不僅發生了射流的斷裂,還發生了小部分粒子的反射以及穿透之后的飛散現象。

圖5 聚能射流擊穿雙殼的過程Fig.5 Penetration of shaped charge jet onto double-shells

3.2.2射流頭部速度衰減規律

圖6為射流頭部的速度時間曲線,從曲線中可以看出,速度的趨勢為先增大后減小,之后又增大最終趨于穩定,最大速度可以達到約4 174 m/s。隨著金屬射流的形成,頭部速度不斷增大,t≈66 μs時速度達到最大值,隨后速度開始不斷減??;t≈94 μs時射流開始擊穿第一層殼,從曲線可以明顯看出,此時曲線切線斜率增大,這正是由于第一層殼的阻力作用;t≈168 μs時射流開始擊穿第二層殼,隨后殼的阻力作用以及射流斷裂引起速度迅速減??;t≈182 μs時頭部斷裂部分擊穿第二層殼,頭部速度達到最小值,約為1 871 m/s,之后,射流頭部完全穿透第二層殼,速度開始增大,t≈182 μs后速度趨于穩定,約為2 200 m/s。由此看出,第二層殼的阻力作用對速度衰減的影響較大,其次,射流斷裂也對速度減小有一定的影響。

圖6 射流頭部速度-時間曲線Fig.6 Velocity-time curve of jet head

3.2.3殼體破口形成過程分析

圖7為第一層殼Von Mises Stress圖,從圖中可以看出,聚能裝藥在起爆后約108 μs時形成的射流已經作用于殼體,此時殼體處于沖塞凹陷階段,出現大的塑性變形,此時最大應力約達到1.4 GPa;隨著聚能載荷的進一步作用,殼體由于塑性變形而變?。籺≈120 μs時表面開始出現破口,此時殼體已經處于凹陷圓盤化階段,破口的直徑約為0.021 4 m,將其與炸藥直徑進行無量綱化,表示為C,此時C≈0.133;隨后殼體被射流完全擊穿,殼體中心出現圓形破孔。

圖8為第二層殼Von Mises Stress圖,從圖中可以看出,t≈178 μs時射流頭部開始作用于第二層殼,圖中初始破口的直徑較小,這是由于射流在擊穿第一層殼之后出現了斷裂現象,射流頭部斷裂部分首先對第二層殼作用;隨后在后續射流斷裂塊的作用下,殼體直接被擊穿,t≈200 μs時C≈0.064 5;隨著射流斷裂塊陸續穿透殼體,射流不斷的堆積,斷裂的射流漸漸融合,導致破口不斷增大,且破口形狀不再規則,t≈254 μs時C≈0.136。由此可得,第一層殼在被擊穿過程中,經歷了沖塞、凹陷等過程,第二層殼直接被射流擊穿,之后不斷的被射流斷裂塊擊穿。

圖7 第一層殼Von Mises Stress圖Fig.7 Von Mises Stress contour of the first layer of shells

圖8 第二層殼Von Mises Stress圖Fig.8 Von Mises Stress contour of the second layer of shells

將破口大小與炸藥直徑進行無量綱化,表示為C,兩層殼的破口時間曲線如圖9所示,虛線代表第一層殼,從圖中可以看出,t≈120 μs時第一層殼出現破口,初始破口與炸藥直徑無量綱化之后C≈0.133,隨著射流載荷的作用,破口越來越大;t≈128 μs時射流完全擊穿第一層殼,破口處射流直徑基本保持不變,此時C≈0.164,之后破口基本穩定;t≈254 μs時C≈0.174。實線代表第二層殼,從圖中可以看出,t≈178 μs時第二層殼出現破口,初始破口與炸藥直徑無量綱化之后C≈0.054 9,初始破口小于第一層殼的初始破口,這是因為第二層殼首先被射流斷裂塊擊穿;t≈202 μs時后續的射流斷裂塊作用于殼體,導致破口進一步增大;t≈236 μs射流主體開始擊穿第二層殼,且由于射流不斷的堆積,斷裂的射流漸漸融合,導致破口再一次增大;最后破口大小趨于穩定,t≈254 μs時C≈0.136。由此可得,金屬射流在擊穿雙殼結構的過程中,由于射流斷裂導致第二層殼的初始破口小于第一層殼的初始破口,隨著射流的進一步作用,破口進一步增大,最終第一層殼的破口約為第二層殼破口的1.3倍。

圖9 第一層殼和第二層殼的破口-時間曲線(圖中C指將破口大小與炸藥直徑的無量綱化)Fig.9 Crevasse-time curves

4結論

本文首先采用SPH方法建立半球形聚能裝藥模型,對金屬射流的形成進行分析,得到以下結論:藥型罩在高壓的爆轟壓力作用下發生高速變形,并由于受到軸向拉伸和徑向壓縮,進在軸線上聚合形成聚能射流;射流頭部最大速度可以達到約4 174 m/s,尾部最大速度可以達到約2 010 m/s,射流頭部的速度較高,尾部杵的速度較低,基本上遵循線性的速度分布,符合射流速度分布的基本規律。在此基礎上,建立SPH-FEM半球形聚能裝藥擊穿雙層殼模型,對射流速度以及毀傷過程進行分析,得到以下結論:

(1) 射流頭部速度先增大后減小,在擊穿第二層殼以后速度有所增大但最終趨于穩定,頭部最大速度可以達到約4 174 m/s,第二層殼的阻力作用對速度衰減的影響較大;

(2) 第一層殼在被擊穿過程中,經歷了沖塞、凹陷等過程,第二層殼直接被射流擊穿,之后不斷的被射流斷裂塊擊穿;

(3) 第一層殼的初始破口大于第二層殼的初始破口,約為2.4倍,這是由于第二層殼先被射流斷裂塊擊穿,隨著射流的進一步作用,破口進一步增大,最終第一層殼的破口約為第二層殼破口的1.3倍。

參 考 文 獻

[1] 黃正祥, 張先鋒, 陳惠武,等. 藥型罩錐角對聚能桿式侵徹體成型的影響[J]. 南京理工大學學報, 2005, 29(6): 645-657.

HUANG Zheng-xiang, ZHANG Xian-feng, CHEN Hui-wu, et al.Influence on formed mechanism of jetting projectile charge by liner angle[J]. Journal of Nanjing University of Science and Technology, 2005, 29(6): 645-657.

[2] 曹麗娜. 聚能射流和破甲過程數值模擬方法的研究[D]. 長春: 長春工業大學, 2010.

[3] Perez E, Fauquignon D,Chanteret P. Fundamental studies of shaped charge mechanisms[C]// Proc 3rd Intern Symp on Ballistics. Karlsruhe,Germany, 1977.

[4] 初文華. 處理非連續問題的三維SPH算法及其在沖擊動力學問題中的應用[D].哈爾濱:哈爾濱工程大學,2013.

[5] 賈鑫, 黃正祥, 祖旭東,等. 聚能裝藥垂直侵徹橡膠復合裝甲的變形研究[J]. 工程力學, 2013, 30(2): 451-457.

JIA Xin, HUANG Zheng-xiang, ZU Xu-dong, et al. Research on deformation of rubber composite armor against shaped charge vertical penetration[J]. Engineering Mechanics, 2013, 30(2): 451-457.

[6] Liu G R, Liu M B. Smoothed particle hydrodynamics: a meshfree particle method[D]. Singapore:World Scientific Publishing, 2003.

[7] 鄭平泰, 楊濤, 秦子增. 聚能射流形成過程的理論建模與分析[J]. 國防科技大學學報, 2006, 28(3): 28-32.

ZHENG Ping-tai, YANG Tao, QIN Zi-zeng. Theoretical modeling and analysis of the formation process of shaped charge jet[J]. Journal of National University of Defense Technology,2006, 28(3): 28-32.

[8] Molinari J F. Finite element simulation of shaped charges[J]. Finite Elements in Analysis and Design,2002,38:921-936.

[9] Feng D L, Liu M B, Li H Q,et al. Smoothed particle hydrodynamics modeling of linear shaped charge with jet formation and penetration effects[J]. Computers & Fluids, 2013, 86: 77-85.

[10] Liu M B, Liu G R, Zong Z,et al. Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology[J]. Computers & Fluids, 2003, 32: 305-322.

[11] 溫萬治, 恢壽榕, 趙衡陽,等. 聚能裝藥侵徹鋼板全過程的數值模擬[J]. 爆炸與沖擊, 2001, 21(2): 126-130.

WEN Wan-zhi, HUI Shou-rong, ZHAO Heng-yang, et al. Numerical simulation for penetration of a steel slab by a shaped charge[J]. Explosion and Shock Waves, 2001, 21(2): 126-130.

[12] 張先鋒, 陳惠武. 三種典型聚能射流侵徹靶板數值模擬[J]. 系統仿真學報, 2007,19(19): 4399-4410.

ZHANG Xian-feng, CHEN Hui-wu. Computional study of three typical shaped charge jets[J]. Journal of System Simulation, 2007,19(19): 4399-4410.

[13] 李磊, 沈兆武, 李學嶺,等. SPH方法在聚能裝藥射流三維數值模擬中的應用[J]. 爆炸與沖擊,2012,32(3):316-322.

LI Lei, SHEN Zhao-wu, LI Xue-ling, et al. Application of SPH method to numerical simulation of shaped charge jet[J]. Explosion and Shock Waves, 2012, 32(3): 316-322.

[14] Zhang A M, Yang W S, Yao X L. Numerical simulation of underwater contact explosion[J]. Applied Ocean Research, 2012, 34:10-20.

[15] Zhang A M, Yang W S, Huang C, et al. Numerical simulation of column charge underwater explosion based on SPH and BEM combination[J]. Computers & Fluids, 2013, 71:169-178.

[16] Zhang A M, Ming F R, Wang S P. Coupled SPHS-BEM method for transient fluid-structure interaction and applications in underwater impacts[J]. Applied Ocean Research, 2013, 43: 223-233.

[17] Zhang Z F, Ming F R, Zhang A M.Damage characteristics of coated cylindrical shells subjected to underwater contact explosion[J]. Shock and Vibration, 2014(1):1-15.

[18] Dobratz B M. LLNL Explosive Handbook. UCRL-52997[M]. Livermore, CA: Lawrence Livermore National Laboratory, 1981.

[19] Shin Y S, Lee M, Lam K Y,et al. Modeling mitigation effects of water shield on shock waves [J]. Shock and Vibration, 1998, 5:225-234.

[20] Libersky L D, Randles P W, Carney T C,et al. High strain Lagrangian hydrodynamics: a three-dimensional SPH code for dynamic material response[J]. Journal of Computational Physics, 1993, 109: 67-75.

[21] Johnson G R, Cook W H. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[C]//Proceedings of the 7th International Symposium on Ballistics.Hauge, Netherlands,1983.

[22] 廖海平. 聚能侵徹體對雙層反應裝甲的沖擊起爆 [D]. 南京:南京理工大學, 2003.

基金項目:國家自然科學基金(U1430236;51479041;51279038)

收稿日期:2015-03-30修改稿收到日期:2015-07-29

通信作者張阿漫 男,教授,博士生導師,1981年3月生

中圖分類號:O385

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2016.14.011

Penetration characteristics of hemispherical shaped charge based on SPH-FEM method

ZHANG Zhi-fan1, LI Bing2, WANG Long-kan1, ZHANG A-man1

(1. College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China;2. The 91439th Unit of PLA, Dalian 116041, China)

Abstract:In order to investigate the formation of rod-liked jet and the penetration properties, a SPH-FEM model of hemispherical shaped charge was established to simulate the formation of metal jet and its penetration into double-shells. Through the analysing the formation process of shaped charge jet, the attenuation of velocity and the process of crevasse forming, it is shown that a metal jet with a high speed will be generated after the liner gets deformation under the detonation of shaped charge; the maximum velocity of jet head reaches about 4 174 m/s; during the penetration process of the first layer of shells by the jet, it experiences plug failure, denting, petaling etc. In addition, the jet breakup happens before the jet penetrates the second layer of shells and the initial crevasse occurs. The calculation and analysis presented above may be helpful for designing shaped charges.

Key words:hemispherical shaped charge; SPH-FEM; metal jet; penetration

第一作者 張之凡 女,博士生,1990年1月生

主站蜘蛛池模板: 亚洲无线一二三四区男男| 97在线免费视频| 国产精品白浆无码流出在线看| 亚洲人人视频| 国内自拍久第一页| 午夜无码一区二区三区| 亚洲国产精品无码AV| 九九香蕉视频| 人妻精品久久无码区| 99在线观看免费视频| 国产人碰人摸人爱免费视频| 欧美午夜在线视频| www成人国产在线观看网站| 午夜天堂视频| 国产精品久久久精品三级| 国产在线一区视频| 五月天综合婷婷| 国产精品精品视频| 91激情视频| 狠狠ⅴ日韩v欧美v天堂| 女人18毛片久久| Aⅴ无码专区在线观看| 国产内射在线观看| 久久久久亚洲AV成人人电影软件| 国产高清无码麻豆精品| a在线亚洲男人的天堂试看| 欧美激情综合一区二区| 国产无码性爱一区二区三区| 午夜国产大片免费观看| 日本午夜三级| 8090成人午夜精品| 亚洲国产精品VA在线看黑人| 亚洲精品欧美日本中文字幕| 国产精品色婷婷在线观看| 久久综合九色综合97网| 女人av社区男人的天堂| 国产99久久亚洲综合精品西瓜tv| 99热6这里只有精品| 欧美自拍另类欧美综合图区| 狠狠亚洲婷婷综合色香| h网址在线观看| 日韩欧美91| 91精品情国产情侣高潮对白蜜| 五月激激激综合网色播免费| 精品国产美女福到在线直播| 亚洲美女高潮久久久久久久| 天堂成人av| 久久婷婷五月综合色一区二区| 四虎在线观看视频高清无码 | 色噜噜在线观看| 亚洲AV成人一区国产精品| 1769国产精品免费视频| 亚洲综合狠狠| 婷婷成人综合| 91无码网站| 国产精品福利一区二区久久| 亚洲啪啪网| 九月婷婷亚洲综合在线| 欧洲欧美人成免费全部视频| 91久久大香线蕉| 婷婷丁香色| 国产av无码日韩av无码网站| 无码国产伊人| 午夜爽爽视频| 91网在线| 一本大道在线一本久道| 国产色婷婷视频在线观看| 日韩欧美中文字幕在线精品| 91成人在线观看视频| 露脸真实国语乱在线观看| 97精品国产高清久久久久蜜芽| 国产亚洲欧美另类一区二区| 日韩黄色大片免费看| 日韩国产一区二区三区无码| 亚洲国产成人自拍| 久久一日本道色综合久久| 91在线丝袜| 最新精品国偷自产在线| 欧美中文字幕第一页线路一| 国产精品白浆在线播放| 99久久精品国产自免费| 91在线无码精品秘九色APP|