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

粉末藥型罩聚能射流形成過程中溫度分布及影響分析

2016-10-10 08:06:40呂愿宏王續躍寶圖雅李曉杰王連吉
工程爆破 2016年4期
關鍵詞:測量分析

呂愿宏, 王續躍, 寶圖雅, 李曉杰, 王連吉

(1. 內蒙動力機械研究所, 呼和浩特 010010;2. 大連理工大學 機械工程學院精密與特種加工教育部重點實驗室, 遼寧大連 116024;3. 大連理工大學 工業裝備結構分析國家重點實驗室, 遼寧大連 1160244. 內蒙古航天紅崗機械有限公司, 呼和浩特 010010)

?

粉末藥型罩聚能射流形成過程中溫度分布及影響分析

呂愿宏1, 王續躍2, 寶圖雅4, 李曉杰3, 王連吉2

(1. 內蒙動力機械研究所, 呼和浩特 010010;2. 大連理工大學 機械工程學院精密與特種加工教育部重點實驗室, 遼寧大連 116024;3. 大連理工大學 工業裝備結構分析國家重點實驗室, 遼寧大連 1160244. 內蒙古航天紅崗機械有限公司, 呼和浩特 010010)

根據石油射孔彈的實際結構和幾何尺寸,運用ANSYS/LS-DYNA-2D非線性動力學有限元分析軟件,采用瞬態非線性的熱耦合計算方式,對孔隙率為11.53%的銅藥型罩射流形成過程中典型瞬態的溫度場進行描述和分析;對多孔藥型罩聚能射流的最高溫度-時間曲線進行研究,并對射流自身是否產生熔化進行判斷;對比了多孔藥型罩與密實藥型罩聚能射流軸線和外表面溫度。結果表明:聚能射流軸線溫度高,由軸線向外表面逐漸降低,最高溫度先增大,11μs增到最大1 743K后減小,最后幾乎不變,約為1 378K,多孔藥型罩比密實藥型罩聚能射流的溫度高,延伸性能和穩定性能更好。

石油射孔彈; 數值模擬; 溫度場; 多孔藥型罩; 密實藥型罩; 聚能射流

1 引言

金屬聚能射流的溫度越高,延展性能越好,射流斷裂時間越長,這決定了聚能射流良好的侵徹性能。但目前還沒有一種可以直接測量聚能射流溫度的方法,大都采用估算的方法。劉迎彬等〔1〕主要從多孔藥型罩受到爆轟波加載及卸載作用、多孔藥型罩在壓合時的塑性變形以及聚能射流拉伸過程中的塑性變形等三個方面對聚能射流引起的溫升進行了理論估算,并與密實藥型罩聚能射流的溫升進行比較。陳昊〔2〕用LS-DYNA數值模擬了球缺形紫銅藥型罩射流形成時的溫度變化情況。李裕春等〔3〕利用LS-DYNA-2D對線型聚能射流的形成過程進行了數值模擬,分析了線型聚能射流形成過程中溫度的分布特性。馬上等〔4〕提出自適應分裂質點方案,避免了物質點法中的數值斷裂,更加準確地仿真了射流形成過程,并對射流形成過程中溫度的分布情況進行了分析。

目前,大多是對某一時刻密實藥型罩聚能射流溫度的數值計算進行研究,對其是否產生熔化并未給出可靠判斷。通常,所用的粉末旋壓免燒結工藝壓制的藥型罩會帶有孔隙,為了研究多孔藥型罩聚能射流在自由運動過程中的穩定性,對其溫升進行數值計算是非常重要的。本文選用多孔材料的物態方程,對多孔藥型罩聚能射流的物態進行判斷,并與密實藥型罩聚能射流的溫度進行對比分析。

2 基本假設及多孔度

彈殼、炸藥和藥型罩都是均勻連續介質,計算采用頂部中心點起爆,裝藥結構為嚴格軸對稱結構,采用cm-g-μs單位制。粉末藥型罩采用多孔材料狀態方程進行計算,定義初始多孔度α0為〔5〕:

(1)

式中:ρ0為密實材料在常態下的初始密度;ρ00為多孔材料在常態下的初始密度;V0=1/ρ0;V00=1/ρ00;孔隙率δ=1-1/α0。

3 計算模型

3.1幾何模型

采用LS-DYNA-2D建立數值模型,該模型由外殼、炸藥和金屬藥型罩三部分組成。以普遍使用的某型號石油射孔彈為例〔6〕,裝藥高度為40mm,藥型罩底部外口徑為40mm,罩錐角為60o,罩壁厚為1.5mm,裝藥為28g的PBX9010。射孔彈的整個結構具有軸對稱性,為節約計算資源,將采用1/2有限元計算模型,截面結構和尺寸如圖1所示。

圖1 射孔彈的幾何結構和幾何尺寸Fig.1 Geometric structure and dimension of perforating charge

3.2有限元模型

模型使用SOLID162二維實體單元進行網格劃分,彈殼和炸藥、炸藥和藥型罩之間的接觸使用CONTACT_2D_AUTOMATIC_SURFACE_TO_SURFACE接觸算法〔7〕,藥型罩自身的接觸使用CONTACT_2D_AUTOMATIC_SINGLE_SURFACE接觸算法。由于藥型罩在炸藥作用下形成聚能射流的過程中存在大變形、大應變,使用Lagrange算法會造成單元嚴重畸變,因此需要使用自適應網格技術〔8〕。計算時間為20μs,不刪除炸藥和彈殼的作用。彈殼、炸藥和藥型罩均使用映射劃分網格〔9〕。本模型共生成2 908個節點,2 652個單元,其有限元模型如圖2所示。

圖2 射孔彈的有限元模型Fig.2 Finite element model of perforating charge

3.3材料模型及參數

(2)

式中:P為等熵壓力;E為內能;V為爆轟產物的相對體積;A,B,R1,R2,ω為待定常數。

炸藥PBX9010的計算參數分別為:密度ρ0=1.787 g/cm3,爆速D=0.84 cm/μs,爆壓PJ=0.34 Mbar,A=5.814 Mbar,B=6.8×10-2Mbar,R1=4.1,R2=1.0,ω=0.35,E0=0.09 Mbar。

藥型罩材料為密實和多孔紫銅,均統一采用Mie-Gruneisen狀態方程來描述〔7〕:

(3)

pH(η)= A1η+A2η2+A3η3

(4)

在以上方程中,可以用α控制孔隙度。當α=1時,為密實藥型罩,密實銅的狀態方程參數是〔10〕:A1=1.4057,A2=2.4818,A3=3.5961。當α=1.13時,為孔隙率11.53%多孔藥型罩,筆者在文〔10〕中,采用多孔材料模型對銅在α=1.13孔隙度時的狀態方程進行修正,并擬合得到:A1=2.5921,A2=-4.229,A3=19.5435。

對于外殼材料碳鋼采用Gruneisen狀態方程來描述:

(5)

式中:c為聲速;S1,S2和S3為D -μ曲線斜率的系數;γ0為Gruneisen系數;a為對γ0的一階修正;μ=ρ/ρ0-1。

鋼外殼的物態方程參數為〔11〕:c=0.4569cm/μs,S1=1.49,S2=S3=0,γ0=2.17,a=0.46。

4 模擬結果與分析

為了對聚能射流形成過程進行熱力耦合分析,需要在原來的K文件中加入熱材料關鍵字與控制關鍵字。對于LS -DYNA而言,熱(溫度)、變形和應力三者彼此相關且同時產生,DYNA提供了三者的耦合計算。因為聚能射流形成的過程瞬間完成,高應變率下的變形往往是絕熱過程,所以本文忽略導熱的影響,模型的初始溫度設定為293 K。分析類型為熱力耦合分析,整個模擬過程中采用瞬態非線性的熱耦合計算方式。

4.1聚能射流形成過程中的溫度場

對多孔藥型罩聚能射流形成過程溫度場進行計算,設定計算時間為20 μs。通過數值模擬得到了聚能射流形成過程中不同時刻溫度場云圖。為了便于研究,計算結果隱藏外殼和炸藥,取幾個典型瞬態的溫度場進行觀察和分析,如圖3所示。

圖3 聚能射流形成過程中的溫度場Fig.3 The temperature field during the jet formation process

圖3給出了聚能射流形成過程中6個典型瞬態溫度場云圖,其中:(a)表示多孔藥型罩的初始狀態,初始溫度為293 K。(b)表示1.4 μs時爆轟波到達多孔藥型罩并使其的溫度上升到994.3 K。(c)表示在爆轟產物的作用下藥型罩向對稱軸線運動、碰撞,射流頭部已經形成。射流頭部附近沿對稱軸線呈現出最高溫度,罩頂部的溫度較高,罩底部的溫度低較。(d)表示11 μs時射流對稱軸線的最高溫度達到最大值1 743 K,射流、杵體和兩翼溫度均勻,射流溫度最高,杵體次之,兩翼最低,原因是兩翼只受到爆轟波的沖擊壓縮和爆轟產物的作用,還未在中心對稱軸線發生碰撞。(e)表示14 μs時射流內部軸線上溫度高,由中心向外溫度逐漸降低,從射流的內部軸線最高1 243.26 K降到聚能射流的外表面最低855.46 K。(f)表示隨著射流的自由運動,射流中心最高溫度區域被進一步拉伸擴大,杵體溫度開始變得不均勻。

4.2聚能射流溫度分析

由圖3可以看出,多孔藥型罩形成的聚能射流的最高溫度出現在聚能射流軸線上靠近中心位置處,對射流不同時刻的最高溫度進行測量,測量結果如圖4中的Highest temperature曲線所示。

圖4 銅的沖擊溫度與時間的關系曲線Fig.4 Impact T-t relationship curves of copper

聚能射流不同時刻最高溫度處壓力測量結果見圖5。

圖5 聚能射流最高溫度處的壓力與時間關系曲線Fig.5 P-t relationship curves at highest temperature of shaped charge jet

由圖4中的Highest temperature曲線和圖5可以看出,1 μs時炸藥爆轟波未到達藥型罩,溫度是初始設定值293 K。2 μs時爆轟波作用于藥型罩,多孔藥型罩的孔隙被壓實,粉末壓實產生能量沉積,溫度驟然上升到966 K。 2 ~11 μs,罩體發生碰撞產生射流,射流中心發生巨大的塑形切變,隨著壓力的升高,聚能射流的最高溫度隨時間不斷升高,11 μs時,溫度和壓力同時達到最大值,最高溫度是1 743 K,最大壓力是35 GPa。11 ~14 μs,射流開始進入自由狀態,壓力下降,隨著壓力的下降溫度也下降。14 ~20 μs,射流完全進入自由狀態,壓力較低,溫度幾乎不變,由于自由拉伸溫度略有升高,平均值1 378 K為銅的常壓熔點。

通常,高壓會使材料的熔點增高,簡單的高壓熔點可用Simon公式計算〔12〕:

(6)

式中:Tm0為零壓熔化溫度;Tm是pm壓力下的熔點;a、c是材料常數,由下式決定:

a=Pc=Qδ2/3[eq(1-δ-1/3)-δ2/3]

(7)

式中:Pc為沖擊熔點處的冷壓值;δ=ρ/ρ0K,對于銅,ρ0K=9.05;Q=59.7166GPa≈0.6Mbar;q=9.88854〔13〕;銅沖擊熔化點的密度ρm=12.98〔14〕。δ=12.98/9.05= 1.4343,Pc=1.3684Mbar。

(8)

式中:γ是Gruneisen系數,對于銅,γ=2〔14〕,則c=1.3。

由上式可得到Simon熔化方程中的參數a和c均已知,Tm0=1 356 K。可以求解出圖5中不同沖擊壓強下的熔化溫度,即圖4中的Simon melting point熔化曲線。由圖4可見,9 ~13 μs,聚能射流最高溫度高于Simon公式計算得到的高壓熔點,高溫處局部產生熔化。13 ~20 μs,盡管聚能射流最高溫度略低于Simon公式計算得到的高壓熔點,由于高應變率下的變形往往是絕熱過程,所以聚能射流熔化部分仍然保持熔化狀態。

從圖4中的Highest temperature曲線可以看出,14 ~20 μs聚能射流的最高溫度是略有升高的,共升高了14 K,原因是熱產生的另外一個來源是畸變變形作用。現在對14 μs時射流橫截上的速度進行測量,測量截面如圖6中的1-1′所示。

圖6 聚能射流速度測量截面Fig.6 The velocity measurement cross section of shaped charge jet

設定截面中間的位置點坐標為0,左半軸取負值,右半軸取正值,其截面上的速度分布見圖7。

圖7 聚能射流截面上的點速度Fig.7 point velocity on measurement cross section of shaped charge jet

由圖7可以看出,截面中間位置速度最大,截面兩外端速度最小,最大相差125 m/s,速度梯度加劇拉伸變形,從而使溫度升高。此外,射流截面中間處速度高,外表面低,也是射流頭部呈尖錐狀的原因。

4.3多孔藥型罩和密實藥型罩聚能射流溫度比較

由圖4可以看出,14 μs后聚能射流最高溫度比較恒定,標志著聚能射流基本形成。下面對14 μs時多孔藥型罩和密實藥型罩聚能射流軸線和外表面溫度進行比較和分析。

聚能射流的軸線和外表面溫度測量位置點如圖8中的點1 ~20所示。外表面溫度測量位置點:1 ~8是杵體外表面測量位置點,9 ~15是兩翼外表面測量位置點,16 ~20是射流外表面測量位置點。

圖8 聚能射流溫度測量位置Fig.8 The temperature measurement position of shaped charge jet

14 μs時多孔藥型罩和密實藥型罩聚能射流軸線測量位置點和外表面測量位置點的溫度分布如圖9所示。

圖9 聚能射流溫度的對比圖Fig.9 The temperature comparison of inner and outside of shaped charge jet

由圖9分析可知聚能射流的溫度,見表1。

表1 聚能射流溫度

由圖9和表1可以看出,多孔藥型罩聚能射流軸線和外表面的溫度均比密實藥型罩的高,這主要是由于多孔藥型罩壓垮時孔隙吸收了大量能量,多孔藥型罩受沖擊波加載及卸載引起的。

軸線上靠近射流的溫度比靠近杵體的溫度高,最高溫度出現在射流軸線上位置點13,由于測量位置1和20也在聚能射流軸線上,故溫度較高。射流溫度最高,杵體次之,兩翼最低,原因是兩翼只受到爆轟波的沖擊壓縮和爆轟產物的作用,還未在中心對稱面上發生碰撞。射流和兩翼過渡區溫度梯度大。多孔藥型罩射流外表面平均溫度是855.46 K(582.31 ℃),密實罩是687.48 K(414.33 ℃),多孔罩比密實罩高167.98 K。銅的再結晶退火溫度是500 ~700 ℃〔15〕,多孔藥型罩聚能射流外表面達到了銅的再結晶退火溫度,延伸性能更好。

通過以上對比分析,多孔藥型罩比密實藥型罩聚能射流的溫度高,延伸性能和穩定性能更好。這對提高聚能射流侵徹能力具有重要意義,與多孔藥型罩高穿深的實驗結論〔16〕相符。

5 結論

(1)通過數值模擬得到了帶殼射孔彈多孔藥型罩聚能射流形成過程中的溫度分布:射流溫度最高,杵體次之,兩翼最低;射流內部軸線上溫度高,由中心向外溫度逐漸降低,14 μs時,溫度從射流的中心軸線最高的1 243.26 K降到射流的外表面最低的855.46 K。

(2)分析了聚能射流最高溫度隨時間的變化規律:初始溫度是293 K,隨后,聚能射流最高溫度逐漸升高,最高溫度達到了1 743 K,然后,溫度開始逐漸下降,下降到1 378 K后,溫度幾乎不變,聚能射流在之后的形成過程中局部產生熔化。

(3)多孔藥型罩比密實藥型罩聚能射流軸線測量點平均溫度高170.46 K,杵體外表面測量點平均溫度高258.24 K,兩翼外表面測量點平均溫度高184.71 K,射流外表面測量點平均溫度高167.98 K。多孔藥型罩聚能射流中心軸線的溫度和外表面的溫度均高于密實藥型罩,且14 μs時多孔藥型罩射流外表面溫度582.31 ℃達到再結晶退火溫度500 ~700 ℃,說明多孔藥型罩形成聚能射流后自由運動更加穩定。

〔1〕 劉迎彬,沈兆武,劉天生,等. 聚能粒子流溫升的理論計算[J]. 高壓物理學報,2013,27(6): 877-883.

LIU Ying-bin, SHEN Zhao-wu, LIU Tian-sheng, et al. Theoretical considerations on the temperaturerise of shaped charge particle jet[J]. Chinese Journal of High Press Physics, 2013,27(6):877-883.

〔2〕 陳昊. 聚能金屬射流形成及侵徹過程中的動態變形研究[D]. 南京:南京理工大學, 2012.

CHEN Hao. Research on dynamic deformation in process of jet formation and penetration[D]. Nanjing:Nanjing University of Science and Technology, 2012.

〔3〕 李裕春,吳騰芳,徐全軍,等. 線型聚能裝藥射流形成過程的數值模擬[J]. 解放軍理工大學學報(自然科學版), 2002,3(3):71-75.

LI Yun-chun, WU Teng-fang, XU Quan-jun, et al. Numerical simulation of linear shaped charge jet formation[J]. Journal of PAL University of Science and Technology(Science and Technology), 2002,3(3):71-75.

〔4〕 馬上,張雄. 聚能裝藥射流形成的自適應物質點法模擬[J]. 固體力學學報,2009,30(5):504-508.

MA Shang, ZHANG Xiong. Adaptive material point method for shaped charge jet formation[J]. Chinese Journal of Solid Mechanics,2009,30(5):504-508.

〔5〕 湯文輝,張若棋. 物態方程理論及計算概論[M]. 北京: 高等教育出版社,2008:247.

TANG Wen-hui, ZAHNG Ruo-qi.Introduction to theory and computation of equations of state[M]. Beijing:Higher Education Press, 2008:247.

〔6〕 白錫忠,常熹. 油氣井射孔彈及其應用[M]. 長沙: 湖南文藝出版社,1992:263.

BAI Xi-zhong, CHANG Xi. The well perforator and its application[M]. Changsha:Hunan Art Press,1992:263.

〔7〕 LS-DYNA keyword user's manual[Z]. Livermore:Livermore Software Technology Corporation, 2012:75-88.

〔8〕 張紅松,胡仁喜,康士廷. ANSYS 14.5 LS-DYNA非線性有限元分析實例指導教程[M]. 北京:機械工業出版社,2013:319-324.

ZHANG Hong-song, HU Ren-xi, KANG Shi-ting.ANSYS 14.5 LS-DYNA non-linear finite element analysis examples tutorials[M]. Beijing:Machinery Industry Press, 2013:319-324.

〔9〕 郝好山,胡仁喜,康士廷. ANSYS 12.0 LS-DYNA非線性有限元分析從入門到精通[M]. 北京:機械工業出版社, 2010:102-104.

HAO Hao-shan, HU Ren-xi, KANG Shi-ting. ANSYS 12.0 LS-DYNA non-linear finite element analysisfrom entry to master[M]. Beijing: Machinery Industry Press, 2010:102-104.

〔10〕 呂愿宏,王續躍,李曉杰,等. 粉末藥型罩石油射孔彈侵徹性能的數值模擬[J]. 工程爆破,2015,21(2):1-4.

LV Yuan-hong,WANG Xu-yue,LI Xiao-jie, et al.The numerical simulation on penetrating performance of powdered liner in perforation of oil[J]. Engeering Blasting, 2015,21(2):1-4.

〔11〕 時黨勇,李裕春,張勝民. 基于ANSYS_LS-DYNA8.1進行顯式動力分析[M]. 北京:清華大學出版社,2005:117.

SHI Dang-yong, LI Yu-chun, ZHANG Sheng-min. Based on ANSYS_LS-DYNA8.1 explicit dynamic analysis[M]. Beijing: Tsinghua University Press, 2005:117.〔12〕 經福謙. 實驗物態方程導引[M]. 北京: 科學出版社, 1999:385-386.

JING Fu-qian. Experimental equation of state guidance[M]. Beijing: Science Press, 1999:385-386.

〔13〕 徐錫甲,張萬箱. 實用物態方程理論導引[M]. 北京:科學出版社,1986:535.

XU Xi-jia, ZHANG Wan-xiang. Practical equation of state theory guide[M]. Beijing: Science Press, 1986:535.〔14〕 MEYER M A. 材料的動力學行為[M]. 張慶民,劉彥,黃風雷,等,譯. 北京:國防工業出版社,2006:94.

MEYER M A. Dynamic behavior of materials[M]. ZHANG Qing-min, LIU Yan,HUANG Feng-lei, et al, translation. Beijing: National Defence Industry Press, 2006:94.

〔15〕 胡傳炘. 熱加工手冊[M]. 北京:北京工業大學出版社,2002:1151-1152.

HU Chuan-xin. Thermal processing handbook[M]. Beijing: Beijing University of Technology Press,2002:1151-1152.

〔16〕 李如江,沈兆武,劉天生. 多孔藥型罩聚能射流低炸高大穿深機理研究[J]. 含能材料,2008,16(4):424-427.

LI Ru-jiang, SHEN Zhao-wu, LIU Tian-sheng. Deep penetration mechanism of jet produced by shaped charge with porous liner at low stand off distance[J]. Chinese Journal of Energetic materials, 2008,16(4):424-427.

Temperaturedistributionandinfluenceanalysisduringshapedchargejetformationprocessofpowderliner

LVYuan-hong1,WANGXu-yue2,BAOTu-ya4,LIXiao-jie3,WANGLian-ji2

(1.MachineryDynamicalInstituteofInnerMongolia,Hohhot010010,China;2.KeyLaboratoryforPrecisionandNon-traditionalMachiningTechnologyoftheMinistryofEducation,SchoolofMechanicalEngineering,DalianUniversityofTechnology,Dalian116024,Liaoning,China;3.StateKeyLaboratoryofStructuralAnalysisforIndustrialEquipment,DalianUniversityofTechnology,Dalian116024,Liaoning,China;4.InnerMongoliaAerospaceMachineCorporation,Hohhot010010,China)

Accordingtoactualstructureandgeometrysizeofpetroleumperforatingcharge,thetypicaltransienttemperaturefieldofcopperlinerwith11.53%porosityjetwasdescribedandanalyzedbyusingnonlineardynamicsfiniteelementsoftwareANSYS/LS-DYNA-2Dandnonlineartransientthermalcouplingcalculationmethod.Thehighesttemperature-timecurveofporouslinershapedchargejetwasstudiedandthejetitselfmeltornotwasjudged.Theaxisandtheoutersurfacetemperatureofporouslinerjetandsolidlinerjetwerealsocompared.Theresultshowedthatthetemperatureofaxiswashighest,decreasingfromaxistooutersurface,themaximumpeakoftemperatureinitiallyincreasedto1 743Kin11μs,thendecreased,finallystabilizedat1 738K.Incomparisontosolidcopperliner,thejettemperatureofporouslinerwashigher,andwithbetterstabilityandextensionperformance.

Petroleumperforatingcharge;Numericalsimulation;Temperaturefield;Porousliner;Solidliner;Shapedchargejet

1006-7051(2016)04-0016-06

2016-01-23

國家自然科學基金項目(11272081,51321004)

呂愿宏(1988-),男,碩士,研究方向為石油射孔技術。E-mail: 737052990@qq.com

TD235.22

Adoi: 10.3969/j.issn.1006-7051.2016.04.004

猜你喜歡
測量分析
隱蔽失效適航要求符合性驗證分析
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
滑動摩擦力的測量與計算
測量的樂趣
電力系統及其自動化發展趨勢分析
測量
中西醫結合治療抑郁癥100例分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 欧美亚洲一二三区| 国产剧情国内精品原创| 日本高清有码人妻| 亚洲国产日韩欧美在线| 无码人中文字幕| 婷婷六月综合网| 女同国产精品一区二区| AV在线天堂进入| 天天做天天爱天天爽综合区| 国产免费人成视频网| 久久精品一品道久久精品| 一本久道久久综合多人| 毛片大全免费观看| 国产第一页亚洲| 中文字幕66页| 国产一区二区网站| 99精品这里只有精品高清视频| 欧美在线黄| 亚洲精品不卡午夜精品| 99re视频在线| 国产杨幂丝袜av在线播放| 亚洲一区毛片| 国产杨幂丝袜av在线播放| 国产精品黑色丝袜的老师| 中文字幕无码av专区久久| 青草娱乐极品免费视频| 欧美一级高清片欧美国产欧美| 欧美一区中文字幕| 99国产在线视频| 亚洲色图另类| 国产日韩精品一区在线不卡 | 亚洲性日韩精品一区二区| 亚洲综合片| 91久久大香线蕉| 亚洲日韩精品无码专区97| 精品福利一区二区免费视频| 九色最新网址| 四虎精品国产AV二区| 四虎永久免费在线| 日韩久草视频| 国产内射一区亚洲| 国产成人亚洲综合A∨在线播放| 久久午夜夜伦鲁鲁片无码免费| 国产精品香蕉在线观看不卡| 黄色网在线| 在线日韩一区二区| 亚洲人成亚洲精品| 久久亚洲国产视频| 成人在线视频一区| 国产美女自慰在线观看| 中文字幕啪啪| 国产91小视频在线观看| 精品无码日韩国产不卡av | 99免费视频观看| 啦啦啦网站在线观看a毛片| 国产网站免费看| 久久久久久尹人网香蕉| 99人妻碰碰碰久久久久禁片| 久久综合九色综合97婷婷| 制服丝袜国产精品| 亚洲欧洲日韩久久狠狠爱| 新SSS无码手机在线观看| 日韩精品无码免费专网站| 51国产偷自视频区视频手机观看 | 久草网视频在线| 伊人久久福利中文字幕| 最新精品国偷自产在线| 免费一级无码在线网站 | 久久亚洲黄色视频| 成人国产小视频| 999精品色在线观看| 欧美国产菊爆免费观看| 亚洲区一区| 国产精品亚洲精品爽爽| 尤物精品视频一区二区三区| 在线观看亚洲天堂| 日本高清有码人妻| 麻豆国产原创视频在线播放| jizz国产在线| 在线99视频| 日韩大片免费观看视频播放| 中文字幕久久波多野结衣|