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

惰性芯體物性對階梯形空腔裝藥結構爆轟驅動的影響*

2017-04-05 03:58:56李傳增王樹山宋述忠
爆炸與沖擊 2017年2期

李傳增王樹山宋述忠

(1.北京理工大學爆炸科學與技術國家重點實驗室,北京100081; 2.北京奧信化工科技發展有限責任公司,北京100044)

惰性芯體物性對階梯形空腔裝藥結構爆轟驅動的影響*

李傳增1,2,王樹山1,宋述忠2

(1.北京理工大學爆炸科學與技術國家重點實驗室,北京100081; 2.北京奧信化工科技發展有限責任公司,北京100044)

為研究惰性芯體物性對軸向階梯形空腔裝藥結構爆轟驅動的影響規律,分別加工并裝配了空腔內填充LY12鋁、尼龍1011和無填充物三種技術狀態的戰斗部并進行了靜爆實驗,采用脈沖X射線成像測試技術獲得了破片平均速度和飛散場形態;應用非線性動力學計算軟件LS-DYNA,采用ALE算法進行了數值計算,分析了三種惰性芯體物性對破片場飛散形態、沖擊波壓力和破片初速的影響規律。結果表明:惰性芯體材料的沖擊阻抗較主裝藥的沖擊阻抗越大,對沖擊波壓力的影響越明顯,作用于殼體表面的初始壓力越大,主裝藥推動破片做功的能力越強,破片的速度越大,且惰性芯體物性的影響隨著芯體半徑的增大而增大。

梯形空腔裝藥;爆轟驅動;惰性芯體

爆轟產物驅動是一個產物氣體與固體板、殼運動的耦合問題,物體運動必須與產物流場同時求解。一些學者針對帶有空腔的圓柱形裝藥結構爆轟驅動金屬殼體的規律進行了研究,給出不同的破片初速計算公式。如R.M.Lioyd[1]在帶有剛性芯體裝藥的Gurney修正式基礎上,對空腔的徑向稀疏效應進行了簡單修正,給出了帶有空腔的圓柱形裝藥破片初速的計算公式;李世才[2]采用經典的Gurney公式形式,建立虛擬的載荷系數并通過實驗標定的方法給出了計算圓筒形裝藥破片初速的計算公式。S.Wang等[3]和李傳增[4]建立了瞬時爆轟產物雙向驅動模型,在Gurney公式的基礎上,推導了計算圓柱形空腔裝藥破片初速的工程計算模型。姬聰生[5]和王樹山等[6]提出了一種采用軸向階梯形空腔裝藥結構的動能桿戰斗部,通過多層多段帶空腔的異形結構裝藥控制桿條具有不同的飛散速度,在戰斗部軸截面上形成不同的環形分布場,并采用數值計算和實驗相結合的方法對其破片飛散場特征和速度分布規律進行了研究。然而,在實際的應用中,為保證在發射沖擊加載時的裝藥強度,空腔內往往填充不同材料的惰性芯體,而惰性芯體材料的物理和力學特性對該裝藥結構的爆轟驅動將產生影響。因此,本文結合一定工程背景,分別對空腔內填充LY12鋁、尼龍1011和無填充物3種技術狀態下的預制破片戰斗部進行靜爆實驗和數值計算,從機理上探索惰性芯體物性對階梯形空腔裝藥結構爆轟驅動的影響規律。

1 靜爆實驗

1.1 實驗戰斗部機構及參數

實驗戰斗部結構如圖1所示,主裝藥為注裝B炸藥,采用軸向階梯形空腔裝藥結構;擋圈、內襯和外殼的材料均為2 mm厚的LY12鋁;端蓋和墊圈材料為尼龍1011;預制破片材料為高強度軸承鋼,尺寸為?6 mm×6 mm,共960枚,徑向2層軸向4段;加工并裝配了3發無填充物(編號為1、2、3)、3發填充尼龍1011(編號為4、5、6)、3發填充LY12鋁材料(編號為7、8、9)的戰斗部進行靜爆實驗,起爆方式為頂部中心起爆。

1.2 實驗原理與方法

如圖2所示,將受試戰斗部懸吊于600 mm×600 mm的方形鋼質防護桶內,并保證射線管、戰斗部和防護桶的軸線在同一條直線上,且與底片夾相垂直。爆炸點距射線管的距離為2 050 mm,底片夾與爆心相距820 mm,并在每次實驗后進行校準。立靶及校準完成后,采用制式電雷管起爆戰斗部,同時利用電離漆包線的方法向操作臺傳出同步觸發信號,預先設定延遲時間為150μs時,啟動脈沖450形閃光X射線系統。

圖1 實驗戰斗部結構圖和照片(單位:mm)Fig.1 Structure diagram and photo of testing warhead(unit:mm)

圖2 實驗場地布置Fig.2 Testing site

1.3 實驗結果

為方便描述問題,本文中將惰性芯體按半徑由小到大依次記為A(起爆端)、B、C、D四段。

1.3.1 破片飛散場

2號戰斗部由于射線強度調高只得到了靜止像而沒有得到破片動態飛散像,其余8發戰斗部的破片飛散場特征如圖3所示。

圖3 1 5 0μs時各實驗戰斗部的脈沖X射線成像照片Fig.3 X-ray images at 150μs after initiating

三種技術狀態戰斗部破片飛散場平均尺寸分別為218.4、236.4和249.1 mm,內部空洞平均尺寸分別為97.5、101.1和105.7 mm;4~6號戰斗部與1~3號戰斗部相比,破片分布場密集性有所降低,徑向雙層破片已基本錯開,A段破片的堆積程度明顯降低,在底片上能夠清晰地看到5層破片的分布情況,7~9號戰斗部的破片場進一步展開,D段徑向雙層破片已全部錯開,A段破片的堆積程度較前兩種技術狀態進一步降低,呈現幾個一簇的現象,雖然其破片場展開的較為充分,但由于B、C兩段的速度差相對較小,在底片上只能較為清晰地分辨D段的4層破片。4號和7號均沒有拍攝戰斗部靜止像,其中4號中間部位出現了煙霧,對比分析5號和6號,可判斷出尼龍1011材料在瞬間高溫、高壓的爆炸沖擊載荷作用下迅速發生了塑性變形、破碎、相變甚至氣化現象。而7號芯體的端部外形保持得相對較為完整,表明鋁棒在爆炸加載過程中沒有發生較大的塑性變形,這與實驗中回收的鋁棒情況相吻合。

1.3.2 破片速度測試

由于稀疏波和端蓋效應的影響以及破片之間的摩擦、碰撞和擠壓等原因使破片發生了翻轉,很難分清底片上各破片所屬的段數,但D段破片的時間歷程是清晰的。因此本文中采用數字圖像處理技術對各技術狀態戰斗部的D段各層破片的平均速度和最大速度進行了測量。

(1)測試方法

采用脈沖X射線成像技術進行破片速度參數測量的一般原理為:在關聯時間差Δt=t2-t1范圍內,獲取兩幅不同時刻破片場的分布透視圖像,通過分析這兩幅圖像可以得到任意一個破片相對于某一標記時刻的距離,從而得到該破片在Δt時間內所飛行的距離d,進而得到該破片的大致速度。為減少測量誤差,本文中通過采集已知尺寸參照物的像素數目,計算出實測單位與像素所對應的轉換關系,其轉換關系式為:d/i=l/n,其中d為兩點間距離,i為兩點間像素個數,l為參照物實際尺寸,n為參照物的像素個數。

(2)測試結果

計算時,由于芯體材料LY12鋁棒在爆炸加載過程中外形保持較完整,所以7號實驗中選芯體直徑作為參照尺寸;4號實驗中選取整個底片作為參照尺寸;其余實驗中均選取戰斗部直徑作為參照尺寸。D段4層破片平均速度的測試結果見表1~3。

表1 無填充材料時D段破片速度(1~3號)Table 1 Fragment velocity of warhead 1~3(no filling)

表2 芯體材料為尼龍1011時D段破片速度(4~6號)Table 2 Fragment velocity of warhead 4~6(nylon 1011)

表3 芯體材料為LY12鋁時D段破片速度(7~9號)Table 3 Fragment velocity of warhead 7~9(LY12 Al)

由表1~3可知,三種技術狀態下D段破片的平均速度由小到大依次為:空腔、填充尼龍1011、填充LY12鋁。其中,當惰性芯體材料為尼龍1011和LY12鋁時,D段破片的平均速度較無填充物條件下分別增加了約8.6%和19.4%。該結果表明:帶惰性芯體的軸向階梯形裝藥結構能夠有效控制各段破片之間的速度梯度,且惰性芯體材料對破片速度分布的影響較為顯著。

2 數值計算

2.1 有限元建模方法

LS-DYNA中的多物質ALE流固耦合算法在計算爆轟驅動問題時,炸藥、空氣等流體材料在ALE單元中流動,由于網格固定,不存在單元劇烈畸變的問題,所以能夠很好地描述裝藥對破片的驅動過程[7]。為能夠完整地描述破片場的形成過程,建立全戰斗部模型。預制破片、殼體、內襯和端蓋均采用Lagrange單元進行劃分。主裝藥、空腔或惰性芯體以及空氣域采用ALE算法進行描述,并通過共節點控制物質邊界。整個有限元模型共劃分為734 290個單元,分為970個部分(PART)。為保證主裝藥爆轟產物充分驅動殼體和破片,空氣域應足夠大,本模型中取為主裝藥的3.3倍,并在空氣域邊界施加非反射邊界條件,以模擬無限空氣域。為使該戰斗部的破片場達到預定分布狀態,設計時應充分考慮稀疏波效應的影響,使最小空腔段對應的破片向外充分飛散,而最大空腔段對應的破片應在稀疏波效應的作用下向內拉伸,以減小破片場中間的空洞,所以延長最小空腔段裝藥,并用非金屬墊圈來填補破片留下的空間。有限元計算模型和流固耦合計算模型見圖4,其中空腔技術狀態按填充空氣進行處理。

圖4 有限元模型圖和流固耦合計算模型Fig.4 Numerical simulation drawing and fluid-structure interaction model

2.2 接觸算法

破片式戰斗部起爆后,內襯、殼體、破片、端蓋之間將發生復雜的接觸碰撞,且破片之間的相互碰撞摩擦等作用是不可忽略的。進行數值計算之初曾嘗試使用侵蝕面面接觸算法(ERODING_SURFACE_ TO_SURFACE)進行計算,但出現了破片之間的網格穿透的現象,使計算結果不可信。時黨勇等[8]提出采用基于Pinball搜索方式的侵蝕自動單面接觸算法(ERODING-SINGLE-SURFACE)同時對多個PART進行控制,讓程序自動追蹤各構件可能發生碰撞的區域,且可以有效模擬破片之間的相互摩擦、碰撞與翻滾效應。本文中為避免網格穿透現象的發生并取得和實際相符的計算結果,綜合采用了以上2種算法,在全面設置自動單面接觸算法的同時,也設置了徑向2層破片間的侵蝕面面接觸算法,取得了較好的計算效果。同時,為使上述接觸算法更加魯棒,預制破片之間保留了0.002 mm的微小間距,并定義了材料之間的靜、動摩擦因數。

2.3 材料本構模型及參數

計算時,采用統一的cm-μs-g單位制。炸藥采用MAT-HIGH-EXPLOSIVE-BURN高能炸藥燃燒材料模型和JWL狀態方程進行描述。Jones-Wilkens-Lee(JWL)狀態方程能夠精確地描述在爆炸驅動過程中爆轟產物的壓力、體積和能量特性,其表達式為

式中:p為壓力;V為相對體積;e為內能;A、B、R1、R2及ω為常數。

注裝B炸藥的本構模型和狀態方程參數見表4,其中,ρ為密度;E0為單位體積初始內能;D為爆速;pCJ為爆轟壓力。

空氣采用MAT_NULL材料模型和線性多項式狀態方程:

式中:C0~C6為多項式方程系數;其中,C0=-1.00×10-6,C1=C2=C3=C6=0,C4=C5=0.4;μ=ρ/ρ0 -1,ρ、ρ0分別為密度和初始密度。

外殼體、內襯、墊圈、端蓋和破片均采用MAT_PLASTIC_KENEMATIC隨動硬化模型描述。該模型的表達式為

式中:σy為屈服極限;σ0為初始屈服應力為有效塑性應變;Ep是塑性硬化模量,Ep=EtanE/(EEtan);C、p、β、Etan、E為輸入參數,β為硬化參數,E為彈性模量,Etan為切線模量。各材料具體參數見表5,其中:ν為泊松比,εf為失效應變。

表4 主裝藥材料模型參數Table 4 Parameters of B explosives

2.4 計算結果與分析

采用上述有限元模型和計算算法分別對3種技術狀態的戰斗部進行數值計算,計算時間為150μs。

2.4.1 破片場飛散形態

起爆后,爆轟波以球面波的形式從A段起爆端開始向前傳播,且最先與內襯最頂端接觸,隨后向兩側擴展,內襯開始向外膨脹變形,并發生破碎。由于延長了頂層破片處的裝藥,所以稀疏波效應的影響在一定程度上有所減小,抑制了頂層破片向內飛散的情況。此后,在爆轟波巨大壓力下,內襯不斷變形并擠壓里層破片,里層破片推動外層破片向外膨脹,A段破片率先向外擴張,其他三段依次展開,由于戰斗部的軸向階梯形空腔裝藥結構和爆轟波傳播至各段的時間差,戰斗部初始呈現圓臺狀鼓包。在端蓋處,受端蓋的側向膨脹和稀疏波效應的影響,端部壓力迅速衰減,從而減小了對兩端破片的推動作用,造成戰斗部兩端的破片速度減小,并向戰斗部軸線方向收縮。150μs時,三種技術狀態戰斗部的破片飛散場如圖5所示。圖5(a)和5(b)中后三段破片分布場外部整體輪廓基本呈線性,差別不明顯。而相對于圖5(a)和5(b),圖5(c)中分布場外部整體輪廓明顯具有一定的弧度,且從軸向觀察,D段破片展開后明顯較前兩者稀疏,說明其B、C兩段環形破片場的展開半徑相差不大,速度梯度較小。上述現象與實驗結果吻合較好。

圖5 1 5 0μs時三種技術狀態戰斗部的破片飛散場計算結果Fig.5 Numerically simulated fragment status at 150μs

2.4.2 惰性芯體物性對沖擊波壓力的影響

主裝藥在一端起爆時,爆轟波以球面波的形式向前傳播,由于惰性芯體的存在,爆轟波的傳播過程將受到一定的擾動,從而導致作用于殼體內表面上的沖擊波壓力發生變化而影響破片速度的大小。本文中選取與殼體內表面相鄰的空氣域單元作為觀測點,對沖擊波壓力進行測量,按惰性芯體半徑從小到大的順序將觀測點依次記為A、B、C、D,各觀測點的壓力-時間歷程曲線和沖擊波壓力峰值分別見圖6和表6。

圖6 各觀測點壓力-時間歷程曲線Fig.6 Pressure-time history curve

表7 各觀測點沖擊波壓力峰值Table 7 Shock wave peak pressure

以無填充物技術狀態戰斗部觀測點的沖擊波壓力峰值為參照值,空腔填充LY12鋁和尼龍1011技術狀態戰斗部的沖擊波壓力峰值分別為無填充物的0.99、1.28、2.03、2.75倍和0.98、1.04、1.15、1.19倍。從爆轟物理學的角度分析,當爆轟產物向空腔內部飛散時,在惰性芯體中產生爆炸沖擊波,同時在爆轟產物中產生了反射沖擊波或反射稀疏波,這種反射波的性質取決于主裝藥和惰性芯體的物理特性,如果炸藥的沖擊阻抗ρ0D(ρ0為炸藥的初始密度,D為爆速)比惰性芯體的沖擊阻抗ρm0Dm(ρm0為惰性芯體的原始密度,Dm為沖擊波速度)小,反射時界面上的壓力P將高于爆轟波的C-J壓力PH,則產物中的反射波為沖擊波。反之,如果炸藥的沖擊阻抗大于惰性芯體的沖擊阻抗,界面上的壓力將低于爆轟波的C-J壓力,則傳入產物中的反射波為稀疏波。如果兩者的沖擊阻抗相等,則界面處不發生反射現象,入射波強度則不變地傳入惰性芯體中去。B炸藥的密度約為LY12鋁材料密度的0.6倍,是尼龍1011材料密度的1.61倍。當B炸藥的密度為1.68 g/cm3時,其爆速D為7840 m/s,爆轟波C-J面壓力約為25.82 GPa,根據泰特狀態方程,其爆炸沖擊波在鋁介質中的傳播速度約為7 168 m/s,分界面壓力約為28 GPa,即鋁材料中爆炸沖擊波的初始壓力比炸藥的爆轟波C-J壓力要高,并且隨著材料的沖擊阻抗或動力學剛度的增大而增大。因此,主裝藥的沖擊阻抗大于尼龍1011材料的沖擊阻抗,而小于LY12鋁材料的沖擊阻抗??梢?在爆轟波作用下,LY12鋁材料在該過程中向爆轟產物中反射的為沖擊波,增加了作用在殼體表面上的壓力。相反,當惰性芯體材料為尼龍1011和空腔(空氣)時向爆轟產物中反射的均為稀疏波,造成沖擊波壓力下降,而空腔(空氣)中的的稀疏波效應更強。可見,惰性芯體物性對沖擊波壓力的影響與各芯體材料的沖擊阻抗相關,且其沖擊阻抗較主裝藥的沖擊阻抗越大,對沖擊波壓力的影響越明顯,作用于殼體表面的初始壓力越大。

2.4.3 破片飛散速度

由于單枚破片在飛散的過程中受到相鄰破片之間的碰撞、摩擦以及翻轉等因素的影響,各破片之間的速度存在一定的差異,且現有的計算方法和測試結果所得到的均為破片的平均速度值。因此,本文中對各技術狀態下的1/4有限元模型域所包含的240枚破片的速度進行了統計并求其平均值,統計結果見表7。由于戰斗部結構的對稱性,該統計結果能夠充分反映模擬戰斗部的破片速度分布。表7中的數據表明,各段破片之間存在著明顯的速度梯度,且破片平均速度隨惰性芯體半徑的增加而減小,這與美國海軍實驗場關于帶軸向空腔的圓柱形戰斗部(空腔條件下)的實爆實驗結論相符[9]。A段與D段破片位于戰斗部兩端,由于端蓋的側向膨脹,端部壓力迅速衰減,作用在兩段破片的推動作用減小,因此兩端的破片速度被明顯拉低。

惰性芯體物性對沖擊波壓力的影響規律的分析結果顯示,在戰斗部結構不變的情況下,惰性芯體物性相對主裝藥的沖擊阻抗是造成各段觀測點沖擊波壓力不同的主要原因。LY12鋁材料的沖擊阻抗大于主裝藥的沖擊阻抗,在爆轟波作用下,其向爆轟產物內部反射沖擊波,這在一定程度上增強了作用于外殼體表面上的壓力,進一步推動了破片向外飛散。相反,當惰性芯體材料為尼龍1011材料和空氣時,則向爆轟產物內反射稀疏波,造成作用于殼體表面的壓力下降,弱化了對殼體的推動作用。同時,根據瞬時爆轟的假設,惰性芯體的存在破壞了爆轟產物流場的流動平衡,惰性芯體的沖擊阻抗越大,在爆炸沖擊載荷作用下獲得速度所需的壓強越大,從而促使更多的爆轟產物對外做功推動破片飛散。因此,惰性芯體的沖擊阻抗越大,主裝藥推動破片做功的能力越強,破片的速度越大。

表7 各段破片的平均速度Table 6 Average velocity of fragment

3 結 論

(1)靜爆實驗結果表明:帶惰性芯體的軸向階梯形裝藥結構能夠有效地控制各段破片之間的速度梯度,且芯體物性對戰斗部破片速度的影響較為顯著,當填充芯體材料為尼龍1011和LY12鋁時,破片初速比空腔條件下分別提高了約8.6%和19.4%。

(2)應用非線性動力學計算軟件LS-DYNA,采用ALE算法,分別建立了實驗戰斗部的有限元模型,計算結果表明:破片飛散場形態與實驗結果吻合良好;惰性芯體物性對軸向階梯形裝藥結構爆轟驅動的影響與各芯體材料的沖擊阻抗相關,其沖擊阻抗較主裝藥的沖擊阻抗越大,對沖擊波壓力的影響越明顯,作用于殼體表面的初始壓力越大,主裝藥推動破片做功的能力越強,破片的速度越大,且芯體物性的影響隨著芯體半徑的增大而增大。

[1] Lioyd R M.Conventional warhead system physics and engineering design[M].USA:American Institute of Aeronauts and Astronautics,1998.

[2] 李世才.HQ-9戰斗部的總體設計[D].北京:北京理工大學,1992.

[3] Wang S,Li C.KE-Rod initial velocity of hollow cylindrical charge[J].Defence Science Journal,2011,61(1):25-29.

[4] 李傳增.中口徑艦炮前向攔截反導彈藥技術研究[D].北京:北京理工大學,2011.

[5] 姬聰生.動能桿戰斗部爆轟驅動研究[D].北京:北京理工大學,2004.

[6] 王樹山,李傳增.前向攔截形反導彈藥系統及其實現方法:中國,ZL 201010051895.3[P].2014-6-10.

[7] 時黨勇,李裕春,張勝民.基于ANSYS/LS-DYNA 8.1進行顯式動力分析[M].北京:清華大學出版社,2005:264-282.

[8] 時黨勇,張慶明,夏長富.多層預制破片戰斗部數值仿真方法及起爆方式影響[J].解放軍理工大學學報(自然科學版),2009,10(6):553-558. Shi Dangyong,Zhang Qingming,Xia Changfu.Numerical simulation method and different initiation modes for prefabricated multilayerfragment warhead[J].Journal of PLA University of Science and Technology(Natural Science Edition),2009,10(6):553-558.

[9] 隋樹元,王樹山.終點效應學[M].北京:國防工業出版社,2000.

Influence of inert core stuffing’s physical properties on the impact of detonation driving of scalar hollow charge

Li Chuanzeng1,2,Wang Shushan1,Song Shuzhong2
(1.State Key Laboratory of Explosion Science and Technology Beijing Institute of Technology,Beijing100081,China; 2.Beijing Auxin Chemical Technology Ltd,Beijing100044,China)

In the present work,to find out how the inert core's physical properties influence the detonation driving of the axial scalar hollow charge,we fabricated three warheads separately stuffed with air,nylon 1011 and aluminum 12,carried out a static explosion test using the pulse X-ray imaging testing technology,and obtained the average velocity of the fragment field and the characteristics of the velocity distribution.Also,using the LS-dyna software and the ALE algorithm,we carried out numerical simulation to analyze the influences of the physical properties of three kinds of inert core’s stuffings on the fragment field's shape,shockwave pressure,and initial velocity.The results show that the detonation driving is related with the shock resistance and dynamic stiffness of the stuffing:if the shock resistance of the stuffing is greater than the main charge,the impact on the shock wave pressure is more obvious;the greater the initial pressure on the shell surface,the greater the power of the main charge to drive the fragment,and the faster the fragment’s velocity;and the influence of the inert core’s physical properties become greater as the inert core or cavity radius increases.

scalar hollow charge;detonation driving;inert core

O381國標學科代碼:13035

:A

10.11883/1001-1455(2017)02-0291-08

(責任編輯 王小飛)

2015-04-12;

:2015-07-28

李傳增(1981- ),男,博士,高級工程師,lcz@auxin-tech.com.cn。

主站蜘蛛池模板: 久久婷婷人人澡人人爱91| 欧美成人aⅴ| 亚洲最大综合网| 亚洲成AV人手机在线观看网站| 欧美精品H在线播放| 天天综合网色| 尤物特级无码毛片免费| 国产午夜福利亚洲第一| 毛片基地视频| 真人免费一级毛片一区二区| 亚洲动漫h| 欧美精品另类| 亚洲天堂2014| 青青国产视频| 麻豆AV网站免费进入| 久久久亚洲国产美女国产盗摄| 在线观看av永久| 国产偷国产偷在线高清| 国产精鲁鲁网在线视频| 国产在线精品99一区不卡| 精品国产网站| 成人午夜精品一级毛片| 香蕉在线视频网站| 天天做天天爱天天爽综合区| 欧美一级片在线| 亚洲高清无在码在线无弹窗| 国产精品视频系列专区| 国产亚洲现在一区二区中文| 亚洲AV无码乱码在线观看代蜜桃| 国产成人精品优优av| 91精品日韩人妻无码久久| 91精品专区国产盗摄| 91精品国产91久无码网站| 亚洲欧美一级一级a| 亚洲婷婷六月| 国产亚洲视频中文字幕视频| 一级福利视频| 亚洲福利网址| 精品91视频| 国产精品亚洲综合久久小说| 日本免费精品| 狠狠亚洲婷婷综合色香| 欧美国产成人在线| 91精品啪在线观看国产| 午夜无码一区二区三区在线app| 国产一区二区三区视频| 国产欧美日韩在线在线不卡视频| 都市激情亚洲综合久久| 久久久久国色AV免费观看性色| 男女性色大片免费网站| 亚洲国产天堂久久综合226114| 任我操在线视频| 亚洲精品va| 九九免费观看全部免费视频| 在线播放91| 亚洲第一成年免费网站| 干中文字幕| 日本三级黄在线观看| 国产99视频精品免费视频7 | 中文字幕2区| 婷婷丁香在线观看| 成·人免费午夜无码视频在线观看 | 色老头综合网| 国产精彩视频在线观看| 国产精品男人的天堂| 91麻豆国产精品91久久久| 在线播放国产一区| 亚洲精品不卡午夜精品| 亚洲欧美自拍中文| 亚洲色图欧美视频| 亚洲欧美成人综合| 毛片网站在线播放| 91精品啪在线观看国产| 二级特黄绝大片免费视频大片| 99re热精品视频国产免费| 免费在线不卡视频| 色窝窝免费一区二区三区| 呦女亚洲一区精品| 青青操视频免费观看| 欧美日韩精品一区二区视频| 国产日韩久久久久无码精品| 狠狠色狠狠色综合久久第一次 |