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

爆炸成型彈丸貫穿靶板靶后破片空間分布模型

2022-08-06 05:04:36邢柏陽侯云輝張東江楊趙兵徐偉華朱桂利
國防科技大學學報 2022年4期
關鍵詞:模型

邢柏陽,郭 銳,侯云輝,駱 強,張東江,楊趙兵,徐偉華,朱桂利,崔 浩

(1. 中國航天科技集團有限公司 第七研究院 第七設計部, 四川 成都 610100;2. 南京理工大學 機械工程學院, 江蘇 南京 210094; 3. 西安現代控制技術研究所, 陜西 西安 710065)

深空撞擊載荷是一種通過引爆自身攜帶火藥,將爆炸成型彈丸(explosively formed projectile, EFP)以超高速(2 km/s以上)射向被探測小天體表面,制造人造撞擊坑,暴露其表面風化層以下物質,為其他探測載荷開展探測或采樣提供作業環境的高效式動能探測載荷,對探測小天體內部物質成分和結構特性具有重要意義[1]。但是,在進行深空撞擊載荷的地面試驗時,只能分析反濺破片云的分布情況,如圖1所示[2]。

圖1 深空撞擊載荷的地面試驗[2]Fig.1 Land experiment of deep space impact load[2]

Corvonato等[3]運用Bernoulli雙紐線描述了包括反濺破片云在內的完整破片云,發現反濺破片云與靶后破片云的分布存在某種特定的聯系。另外,由于鋼板的結構較簡單,在進行深空撞擊載荷相關試驗前,有必要先建立EFP貫穿靶板靶后破片空間分布模型。

Piekutowski[4-9]對破片云的結構進行了首次的系統定義,提出了圓柱狀彈丸正撞擊薄板所產生的破片云的模型,此模型使用圓柱直徑、一個徑向速度和四個軸向速度表示破片云的運動。Sch?fer[10]提出了球形彈丸正撞擊薄板產生破片云的模型,并給出了中心破片速度、薄板破片的質心速度、橢球長半軸的速度和橢球短半軸的速度。蔣建偉等[11-13]通過試驗和數值仿真建立了工程計算模型,此模型可以很好地預測EFP產生的靶后破片的質量和速度的分布特性。黃炫寧等[14]運用AUTODYN軟件中光滑粒子動力學(smoothed particle hydrodynamics, SPH)算法,并結合量綱分析及正交設計理論,建立了EFP正侵徹靶后破片云形狀的數學描述模型。在分析靶后破片空間分布時,目前的文獻通常是建立靶后破片云的形狀方程或分析靶后破片云上個別重要點的速度等,還沒有文獻分析靶后破片云中所有破片的速度、質量、動能、數量與其空間位置間的關系,更沒有文獻建立靶后破片來源(靶板和侵徹體)、初始條件(靶板厚度、侵徹體著靶速度)與上述關系間的聯系。因此,有必要在考慮EFP變截面特點的基礎上,建立EFP貫穿靶板形成的靶后破片的空間分布模型。首先,驗證數值仿真方法以及數值仿真結果;然后,分析靶后破片云中所有破片的速度、質量、動能、數量與其空間位置間的關系;最后,用一系列的數值仿真結果驗證上述關系。

1 數值仿真模型

1.1 數值仿真方法驗證

Dalzell指出,當使用AUTODYN-3D對EFP貫穿鋼靶板形成靶后破片進行數值仿真時,SPH方法優于傳統的Euler方法和Lagrange方法[15]。因此本文使用AUTODYN-3D中的SPH方法對EFP貫穿靶板形成靶后破片的過程進行數值仿真。本小節將利用文獻[11-13]中試驗對本文的數值仿真方法進行驗證,文獻[11-13]中的EFP戰斗部均相同,如圖2所示。EFP戰斗部的裝藥直徑為56 mm,藥型罩為壁厚3 mm的球缺形罩,裝藥為長徑比0.86的壓裝JH-2,殼體厚度為2.52 mm。

(a) 殼體 (b) 主裝藥(a) Shell (b) Main charge

(c) 藥型罩 (d) 壓緊環(c) Liner (d) Press-on ring圖2 EFP戰斗部[11-13]Fig.2 The EFP warhead [11-13]

文獻[11-13]中EFP的X光照片如圖3(a)下半部分所示,本文所建立的EFP數值仿真模型如圖3 (a)上半部分所示。圖3 (b)中的紅色部分為靶板的固定約束,靶板和EFP的粒子直徑均為0.5 mm。

(a) EFP形狀 (b) EFP和靶板(a) Shape of EFP (b) EFP and target圖3 數值仿真模型Fig.3 Numerical simulation model

文獻[11-13]中EFP正侵徹靶板,靶板厚度為20 mm,EFP著靶速度為2 120 m/s,靶板長、寬均為200 mm,本文對此試驗進行復現,數值仿真得到的靶后破片云如圖4上半部分所示,試驗所拍的X光照片如圖4下半部分所示。

對比圖4(a)和圖4(b)的上下部分可以發現,數值仿真復現的結果與試驗結果非常吻合,這表明本文的數值仿真方法是具有可靠性的。

(a)80μs時刻的靶后破片云 (b)95μs時刻的靶后破片云(a)Behind-armor debris clouds at 80μs (b)Behind-armor debris clouds at 95μs圖4 靶后破片云比較Fig.4 Comparison of behind-armor debris clouds

對比95 μs時刻數值仿真復現的靶后破片云徑向膨脹速度與試驗結果,如表1所示。

表1 靶后破片云徑向和頭部膨脹速度Tab.1 Radial and head expansion velocity of the behind-armor debris clouds

由表1可以發現,數值仿真復現的結果與試驗結果非常吻合,偏差不超過2%,這表明本文的數值仿真方法是具有可靠性的。

1.2 數值仿真結果驗證

為了驗證本文數值仿真結果的可信性,針對某典型EFP進行了貫穿靶板的試驗。EFP的激光高速攝影如圖5上半部分所示,所建立的數值仿真模型如圖5下半部分所示。通過幕布墻上已知尺寸的標尺可以換算得到某典型EFP最大長度為98.0 mm,最大半徑為28.5 mm。

圖5 試驗結果和數值仿真模型Fig.5 Experiment result and simulation model

數值仿真模型中靶板的材料為均質裝甲鋼,狀態方程假設為線性;EFP的材料為銅,狀態方程假設為Gruneisen,參數取自文獻[16-17]。在靶板厚度介于30~70 mm且EFP著靶速度介于1 650~1 860 m/s的條件下,當EFP著靶后經過0.5 ms,侵徹過程已經徹底完成,此時靶后破片的速度、質量、動能、數量以及空間位置已經穩定,不再變化,故取EFP著靶后0.5 ms時刻的數值仿真結果進行統計分析。對比分析SPH粒子尺寸為0.3 mm、0.4 mm以及0.5 mm時的數值仿真結果發現,SPH粒子尺寸為0.5 mm時的結果具備足夠的精度并可以節約大量的時間成本。

EFP正侵徹60 mm厚的靶板,測得EFP著靶速度約為1 650 m/s,著靶后0.5 ms時刻的靶后破片云激光高速攝影如圖6 (a)所示,對應時刻的數值仿真結果如圖6 (b)所示。

(a) 激光高速攝影 (b) 數值仿真(a) Laser high speed photography (b) Numerical simulation圖6 EFP貫穿靶板形成的靶后破片云Fig.6 Behind-armor debris cloud formed by the perforation of EFP through the target

由圖6可以發現,數值仿真得到的靶后破片云最前端的銅-鋼粘結體是由EFP殘體和沖塞塊頂部共同組成的,這與激光高速攝影得到的圖像是相吻合的。另外,二者的靶后破片云輪廓也是接近的。激光高速攝影中靶后破片云的上下不對稱,主要是由于EFP的著靶距離很大(100 m),在飛行過程中難免會存在一定的俯仰偏航,使其在著靶時并非嚴格意義上的正侵徹。盡管EFP在著靶前存在的著靶角很小,但卻使得數值仿真的結果與試驗所得的激光高速攝影照片有一定差異。但是總體而言,數值仿真中的靶后破片云與試驗所得的靶后破片云是基本吻合的,這表明本文的數值仿真結果是具有可信性的。

對比試驗中靶后破片云最前端的銅-鋼粘結體的軸向速度與數值仿真結果,如表2所示。

表2 銅-鋼粘結體的軸向速度Tab.2 Velocity of the steel-copper compound

由表2可以發現,數值仿真所得的靶后破片云最前端的銅-鋼粘結體的軸向速度與試驗結果相吻合,偏差不超過3 %,這再次表明本文的數值仿真結果是具有可信性的。

將數值仿真得到的鋼破片質量與試驗結果以及通過文獻[18-19]中的理論分析得到的結果作對比,如表3所示。在試驗中為了兼顧回收效率以及回收精度,采用了2 mm×2 mm孔徑的篩網,可以保障質量不小于0.5 g的靶后破片都能被獲取,故僅統計不小于0.5 g的靶后破片。同時,由于試驗中銅破片質量為174.90 g,與數值仿真結果765.33 g以及理論分析結果719.80 g差距較大,故不考慮銅破片。

表3 靶板產生的靶后破片質量比較Tab.3 Comparison of mass of behind-armor debris from target

由表3可以發現,數值仿真結果與試驗統計結果和理論分析吻合,表明本文的數值仿真結果是具有可信性的。

綜上所述,本文的數值仿真結果是具有可信性的,可以利用該數值仿真模型做進一步的分析與研究。

2 靶后破片空間分布模型

后文中的下標t表示靶板產生的靶后破片,下標p表示EFP產生的靶后破片。H0re=H0/D表示相對靶板厚度,H0表示靶板厚度,D=125 mm表示EFP戰斗部裝藥直徑。v0re=v0/vint表示相對著靶速度,v0表示EFP著靶速度,vint=1 950 m/s表示EFP成型速度。EFP正侵徹靶板,“不同相對靶板厚度”表示相對著靶速度為0.846 2,相對靶板厚度分別為0.240 0、0.320 0、0.400 0、0.480 0和0.560 0;“不同相對著靶速度”表示相對著靶速度分別為0.846 2、0.861 5、0.892 3、0.923 1和0.953 8,相對靶板厚度為0.320 0。前文驗證了數值仿真結果的可信性,為了分析更大的質量范圍,故統計質量不小于0.1 g的靶后破片。

2.1 最大相對速度和相對初始條件的關系

vremax=vmax/vint表示靶后破片云中破片的最大相對速度,vmax表示靶后破片云中破片的最大速度。最大相對速度和相對初始條件的關系可以通過一個平面加以表示,如圖7所示。圖7及后文的圖中,Target表示靶板產生的靶后破片,EFP表示EFP產生的靶后破片,后綴Fit表示擬平面。

(a) 靶板產生的靶后破片(a) Behind-armor debris from target

(b) EFP產生的靶后破片(b) Behind-armor debris from EFP圖7 最大相對速度和相對初始條件的關系Fig.7 Relationship between the relative maximum velocity and the relative initial condition

圖7的擬合平面表達式如式(1)~(2)所示。

vremaxt=-1.233H0re+1.195v0re+0.056 87

(1)

vremaxp=-1.246H0re+1.228v0re+0.043 04

(2)

式(1)~(2)的成立條件均為0.240 0≤H0re≤0.560 0,0.846 2≤v0re≤ 0.953 8,確定系數均為0.98,表明擬合效果良好。由式(1)~(2)可以發現,最大相對速度隨相對靶板厚度的減小和相對著靶速度的增加而增加。

(3)

2.2 相對最遠空間位置和相對初始條件的關系

xremax=xmax/D表示靶后破片云中破片的相對最遠空間位置,xmax表示靶后破片云中破片的最遠空間位置,由于本文僅涉及EFP正侵徹靶板,因此該空間位置通常為靶后破片云頂端。相對最遠空間位置和相對初始條件的關系可以通過一個平面加以表示,如圖8所示。

(a) 靶板產生的靶后破片(a) Behind-armor debris from target

(b) EFP產生的靶后破片(b) Behind-armor debris from EFP圖8 相對最遠空間位置和相對初始條件的關系Fig.8 Relationships between the relative farthest spatial position and the relative initial condition

圖8的擬合平面表達式如式(4)~(5)所示。

xremaxt=-8.592H0re+9.017v0re+0.075 03

(4)

xremaxp=-8.158H0re+9.204v0re-0.324 0

(5)

式(4)~(5)的成立條件為0.240 0≤H0re≤0.560 0,0.846 2≤v0re≤0.953 8,確定系數分別為1.00、0.99。由式(4)~(5)可以發現,相對最遠空間位置隨相對靶板厚度的減小和相對著靶速度的增加而增加。

2.3 相對速度和相對空間位置的關系

vre=v/vmax表示靶后破片云中破片的相對速度,v表示靶后破片云中破片的速度;xre=x/xmax表示靶后破片云中破片的相對空間位置,x表示靶后破片云中破片的空間位置。相對速度和相對空間位置的關系如圖9所示。

從圖9(a)~(b)可以發現,vre與xre的關系呈線性,但不同H0re條件下的斜率不同。從圖9(c)~(d)可以發現,vre與xre的關系呈線性,不同v0re條件下的斜率無明顯差異。另外,可以認為最遠空間位置處的靶后破片擁有最大的速度,因此擬合線經過(1,1)點。采用式(6)~(8)對相對速度和相對空間位置的關系進行描述。

vre=λxre+(1-λ)

(6)

(7)

(8)

(a) 不同靶板厚度條件下靶板產生的靶后破片(a) Behind-armor debris from target with different target thickness

(b) 不同靶板厚度條件下EFP產生的靶后破片(b) Behind-armor debris from EFP with different target thickness

(c) 不同著靶速度條件下靶板產生的靶后破片(c) Behind-armor debris from target with different EFP speed

(d) 不同著靶速度條件下EFP產生的靶后破片(d) Behind-armor debris from EFP with different EFP speed圖9 相對速度和相對空間位置的關系Fig.9 Relationships between the relative velocity and the relative spatial position

式(6)~(8)的成立條件為0.240 0≤H0re≤0.560 0,0.846 2≤v0re≤0.953 8。式(7)~(8)確定系數依次為0.98、0.99。從式(7)~(8)可以發現,隨著H0re的增加,λt先減小后增加,λp一直增加,并且始終有λt>0、λp>0以及vre始終隨著xre的增加而增加。

2.4 相對質量和相對空間位置的關系

mre=m/mto表示靶后破片云中破片的相對質量,m表示靶后破片云中空間位置小于x的破片的質量,mto表示靶后破片云中破片的總質量。相對質量和相對空間位置的關系如圖10所示。

(a) 不同靶板厚度條件下靶板產生的靶后破片(a) Behind-armor debris from target with different target thickness

(b) 不同靶板厚度條件下EFP產生的靶后破片(b) Behind-armor debris from EFP with different target thickness

(c) 不同著靶速度條件下靶板產生的靶后破片(c) Behind-armor debris from target with different EFP speed

(d) 不同著靶速度條件下EFP產生的靶后破片(d) Behind-armor debris from EFP with different EFP speed圖10 相對質量和相對空間位置的關系Fig.10 Relationships between the relative mass and the relative spatial position

從圖10(a)~(b)可以發現,mre與xre的關系呈冪函數,但不同H0re條件下對應的冪函數關系明顯不同。從圖10(c)~(d)可以發現,mre與xre的關系呈冪函數,不同v0re條件下對應的冪函數關系無明顯差異。另外,由橫縱坐標的定義可知,各個工況下的擬合線均經過(1,1)點。采用式(9)~(11)對相對質量和相對空間位置的關系進行描述。

(9)

(10)

(11)

式(9)~(11)的成立條件為0.240 0≤H0re≤0.560 0,0.846 2≤v0re≤0.953 8。式(10)~(11)確定系數分別為0.84、0.99。從式(10)~(11)可以發現,隨著H0re的增加,βt一直增加,βp一直減小,并且始終有βt>0,βp>0。

2.5 相對數量和相對空間位置的關系

Nre=N/Nto表示靶后破片云中破片的相對數量,N表示靶后破片云中空間位置小于x的破片的相對數量,Nto表示靶后破片云中破片的總數量。相對數量和相對空間位置的關系如圖11所示。

(a) 不同靶板厚度條件下靶板產生的靶后破片(a) Behind-armor debris from target with different target thickness

(b) 不同靶板厚度條件下EFP產生的靶后破片(b) Behind-armor debris from EFP with different target thicknes

(c) 不同著靶速度條件下靶板產生的靶后破片(c) Behind-armor debris from target with different EFP speed

(d) 不同著靶速度條件下EFP產生的靶后破片(d) Behind-armor debris from EFP with different EFP speed圖11 相對數量和相對空間位置的關系Fig.11 Relationships between the relative number and the relative spatial position

由圖11(a)~(b)可以發現,當v0re一定而H0re變化時,Nre隨xre的增加呈冪函數增加。由圖11(c)~(d)可以發現,當H0re一定而v0re變化時,Nre隨xre的增加呈冪函數增加。另外,由橫縱坐標的定義可知,各個工況下的擬合線均經過(1,1)點。采用式(12)對相對數量和相對空間位置的關系進行描述。

(12)

采用式(12)描述各工況時的確定系數如表4所示。從表4可以發現,確定系數均接近于1.00,表明使用式(12)作為上述9種工況下的18組數據的擬合公式是合理的。

表4 各工況下的確定系數Tab.4 Working condition and determination coefficient

3 靶后破片空間分布模型驗證

由于在試驗中極難準確地測量并統計所有靶后破片的速度(質量、動能、數量)與空間位置的關系、區分靶后破片來源(靶板或EFP),因此進行2組數值仿真,以驗證EFP貫穿靶板靶后破片空間分布模型的可靠性。EFP正侵徹靶板,第一組初始條件為H0re=0.480 0、v0re= 0.892 3,第二組初始條件為H0re=0.560 0、v0re=0.861 5。由式(1)~(2)可以得到最大相對速度,由式(4)~(5)可以得到相對最遠空間位置,其結果列于表5~6,并與相應的數值仿真結果進行對比。

表5 最大相對速度和相對最遠空間位置(第一組)Tab.5 Relative maximum velocity and relative farthest spatial position (the 1st group)

表6 最大相對速度和相對最遠空間位置(第二組)Tab.6 Relative maximum velocity and relative farthest spatial position (the 2nd group)

由表5~6可以發現,由空間分布模型得到的最大相對速度和相對最遠空間位置與相應的數值仿真的結果均是比較吻合的,表明空間分布模型是可靠的。

設Ekre表示空間位置介于xre1和xre2之間的靶后破片云中破片的相對動能,Ekre=Ek/Eke;Eke表示靶后破片云中破片的最大動能,Eke=0.5·mto·vmax2;Ek表示空間位置介于xre1和xre2之間的靶后破片云中破片的真實動能,Ek=0.5 [(mre1-mre2)mto]·[(vre1+vre2)vmax/2]2;mre1-mre2表示空間位置介于xre1和xre2之間的靶后破片云中破片的相對質量;0.5(vre1+vre2)表示空間位置介于xre1和xre2之間的靶后破片云中破片的平均速度。相對動能與相對空間位置的關系如式(13)。

(13)

第一組的相對動能與相對空間位置的關系如圖12所示,曲線由空間分布模型得到,數據點由數值仿真得到,其中,βt=1.951,λt=1.114,βp=4.537,λp=1.108,Δxre=0.1。

圖12 相對動能與相對空間位置的關系(第一組)Fig.12 Relationship between the relative kinetic energy and the relative spatial position (the 1st group)

第二組的相對動能與相對空間位置的關系如圖13所示,曲線由空間分布模型得到,數據點由數值仿真得到,其中,βt=2.093,λt=1.208,βp=2.826,λp=1.154,Δxre=0.1。

圖13 相對動能與相對空間位置的關系(第二組)Fig.13 Relationship between the relative kinetic energy and the relative spatial position (the 2nd group)

在圖12~13中,相對位置介于nΔxre和(n+1)Δxre之間的所有靶后破片的相對動能之和,以橫坐標為xre=0.05+nΔxre的點表示,其中n=0,1,2,3,4,5,6,7,8,9。由于此次數據統計中僅收集了10個點,所以曲線不夠光滑。從圖12~13中可以發現,分布模型均可以有效地預測數據點的變化趨勢,以上均表明空間分布模型是可靠的。

第一組的相對數量和相對空間位置的關系如圖14所示,曲線由空間分布模型得到,數據點由數值仿真得到。

圖14 相對數量和相對空間位置的關系(第一組)Fig.14 Relationship between the relative number and the relative spatial position (the 1st group)

第二組的相對數量和相對空間位置的關系如圖15所示,曲線由空間分布模型得到,數據點由數值仿真得到。

圖15 相對數量和相對空間位置的關系(第二組)Fig.15 Relationship between the relative number and the relative spatial position (the 2nd group)

圖14中靶板和EFP的確定系數分別為0.97和0.98,圖15中靶板和EFP的確定系數分別為0.89和0.95,均表明空間分布模型是可靠的。同時,可以發現相對數量隨相對空間位置的增加而呈指數增加。

綜上所述,空間分布模型是可靠的。

4 結論

本文建立了EFP貫穿靶板靶后破片空間分布模型,考慮了EFP變截面的特點,可以在區分靶后破片來源(來自靶板或EFP)的基礎上快速預測靶后破片云中所有破片的速度(質量、動能、數量)與其空間位置的定量關系,可以為建立毀傷評估軟件系統奠定基礎,為深空撞擊載荷相關試驗的設計提供一定的指導,并得到以下幾點規律。

1)最大相對速度、相對最遠空間位置隨相對靶板厚度的減小而增加,最大相對速度、相對最遠空間位置隨EFP相對著靶速度的增加而增加;

2)相對速度總是隨相對空間位置的增加而呈線性增加,相對質量、相對數量總是隨相對空間位置的增加而呈冪函數增加;

3)通過靶后破片空間分布模型得到的最大相對速度和相對最遠空間位置與數值仿真結果偏差的絕對值的均值為4.098%。

4)通過靶后破片空間分布模型得到的相對數量和相對空間位置的關系與數值仿真結果的確定系數的均值為0.95。

致謝

埃因霍溫理工大學(Eindhoven University of Technology)的博士生劉磊在論文撰寫階段與作者們進行了富有成果的討論,謹致謝意!

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: a天堂视频| 伊人久久大香线蕉影院| 特级毛片8级毛片免费观看| 亚洲精品动漫| 天堂网国产| 午夜精品福利影院| 精品国产免费第一区二区三区日韩| 又黄又爽视频好爽视频| 成年人福利视频| 久久精品视频亚洲| 在线播放国产一区| 人人91人人澡人人妻人人爽 | 日日碰狠狠添天天爽| 97超碰精品成人国产| 久久无码av三级| 亚洲综合婷婷激情| 亚洲熟女中文字幕男人总站| 国内精品免费| 日韩区欧美国产区在线观看| 一区二区三区成人| 在线国产毛片手机小视频| 狠狠干欧美| 免费无码一区二区| 免费观看亚洲人成网站| 国产美女久久久久不卡| 激情综合网址| 99精品国产电影| 全部毛片免费看| 中文字幕在线看| 久久99精品久久久久久不卡| 在线观看欧美精品二区| 国产亚洲欧美在线中文bt天堂| 亚洲第一精品福利| 亚洲码一区二区三区| 欧美精品亚洲二区| 欧类av怡春院| 国产精品夜夜嗨视频免费视频| 亚洲成人www| 成人午夜在线播放| 欧美国产日产一区二区| 天堂亚洲网| 亚洲国产精品无码久久一线| 日本www在线视频| 欧美精品xx| 色婷婷狠狠干| 狠狠亚洲婷婷综合色香| 一本久道热中字伊人| 91福利在线观看视频| 野花国产精品入口| 久久大香伊蕉在人线观看热2| 久久99热这里只有精品免费看| 蜜芽一区二区国产精品| 亚洲中文无码h在线观看| 99精品国产自在现线观看| 亚洲成年人网| 欧美午夜在线播放| 亚洲欧美在线综合图区| 成人字幕网视频在线观看| 免费一极毛片| 亚洲AV无码一区二区三区牲色| 亚洲aⅴ天堂| 91网在线| 久久综合亚洲色一区二区三区| 美女被操91视频| 亚洲美女一区二区三区| 国产永久在线视频| 国产午夜一级毛片| 亚洲天堂精品视频| 國產尤物AV尤物在線觀看| 亚洲va精品中文字幕| 色婷婷国产精品视频| 婷婷色在线视频| 91精品情国产情侣高潮对白蜜| 国产精品极品美女自在线网站| 欧美一级高清片久久99| 人妻精品久久无码区| 色悠久久久久久久综合网伊人| a毛片免费在线观看| 日本不卡在线播放| 国产精品无码翘臀在线看纯欲| 亚洲色图综合在线| 日本不卡在线播放|