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

鋁粉塵云團爆轟溫壓效應的數值模擬

2018-03-01 01:10:22昝文濤洪滔董賀飛
兵工學報 2018年1期
關鍵詞:區域

昝文濤, 洪滔, 董賀飛

(1.北京理工大學 機電學院, 北京 100081; 2.北京應用物理與計算數學研究所, 北京 100094)

0 引言

溫壓武器作為一種新型武器近些年來得到了廣泛的關注,其利用高溫、高壓可以造成大范圍殺傷效果,鋁粉塵由于其質量輕、含能高的特點被廣泛用于溫壓武器固體藥劑中。而近些年來鋁粉塵在工業領域造成的爆炸事故也越來越多,造成了巨大的人員傷亡及經濟損失,因此研究鋁粉塵云團的爆轟溫壓效應有著重要的意義。

裴明敬等[1]通過實驗研究了含鋁溫壓燃料的能量釋放效率及爆炸沖擊波的形成及擴散過程,發現鋁粉對沖擊波能量貢獻很大,能量釋放率較高;鄭波等[2]通過高速分析系統觀察了溫壓炸藥的拋灑過程,分析了溫壓炸藥拋灑規律;李秀麗等[3]通過外場實驗研究了云爆劑的爆炸沖擊波參數,發現裝藥量30 kg的實驗彈爆炸火球直徑達到17.4 m,并存在兩個正壓作用區;洪滔等[4-6]研究了鋁顆粒的點火機制,分析了鋁顆粒的點火溫度及懸浮鋁粉塵爆轟波參數,并模擬了爆轟波管中的鋁粉塵爆轟過程。

由于溫壓武器毀傷過程包含拋灑、起爆等復雜過程,目前對于溫壓燃料殺傷效果的研究主要以實驗為主,對于其毀傷效應的數值模擬基本以商業軟件模擬含鋁溫壓炸藥爆炸為主。

為此,本文中采用時空守恒元解元(CE/SE)算法求解兩相流方程組,編制程序模擬了半徑3 m的均勻鋁粉塵云團在空氣中爆炸后的沖擊波傳播過程及產生的高溫火球特性,分析了鋁粉塵云團爆轟的溫壓效應。

1 模型

模擬過程采用了兩相流模型,假設顆粒為球形,初始直徑都相同,單個顆粒的溫度都是均勻的,忽略粒子間的作用,粒子化學反應產生的能量假定都被氣體吸收,氣體的組分都是均勻混合的,忽略粒子與氣體間的輻射作用,固相和氣相都滿足各自守恒方程[7-9]。

氣相守恒方程:

(1)

(2)

(3)

(4)

式中:下角標g、a分別代表氣體和鋁顆粒;ρ為密度;u為速度;e為內能;p為壓力;φ為氣體和固體顆粒體積分數,φg+φa=1;qa為單位質量的鋁粒子反應能;Q為氣體與粒子間的熱傳導,

(5)

λ為氣體導熱系數,Nu為Nusselt數。

固相守恒方程:

(6)

(7)

(8)

(9)

(10)

式中:I為單位體積內顆粒的質量變化率;n為單位體積內粒子數;F為氣體對粒子的壓力,

(11)

(12)

Cd為拖曳系數,

(13)

Re為雷諾數。

氣體組分方程:

(14)

式中:j為氣體組分O、N、Al;y為組分濃度。

氣體狀態方程:

(15)

式中:w為分子量;R為氣體常數;T為氣體溫度;k為氣體中組分數量。

對于鋁顆粒,采用如(16)式和(17)式所示的燃燒模型[11],由于鋁顆粒表面被Al2O3包裹,鋁顆粒溫度達到Al的熔點時鋁顆粒膨脹發生破裂從而點火。

(16)

(17)

式中:r為粒子半徑;Ti為點火溫度;d0為粒子初始直徑;Φ為氣體中氧氣的摩爾份額;m=1.75.

本文中所有物理量均采用國際單位制,模擬采用了CE/SE算法[12-14]求解歐拉方程組,對于兩相流方程組的源項,使用4階龍格- 庫塔方法求解。

2 數值模擬結果

2.1 程序驗證

對于理想氣體,氣體常數1.4,分別計算了網格數量為400、300、200、100和32的算例,計算結果如圖1所示。

圖1 t′=0.2時Sod激波管壓力Fig.1 Pressures of Sod shock tube for t′= 0.2

從圖1中可以看出,網格數量200時可以很好地模擬Sod問題,與理論解符合較好,網格數量繼續減小到32時也能模擬Sod激波管問題,但是與理論解的相比,波陣面寬度變寬了。

模擬了鋁粉塵在管道內的爆轟波傳播過程,鋁粉塵密度為0.304 kg/m3,顆粒半徑為1.7 μm,管道直徑為15.2 cm,從左端起爆,左端閉合,右端開口,上下為固壁,模擬得到的爆轟波速度D=1 630 m/s,詳見文獻[7], Tulis等[15]由實驗中得到的鋁粉塵爆速為1 650 m/s.

2.2 計算條件

本文數值模擬了半徑為3.0 m的鋁粉塵云團均勻分布在空氣中。鋁粉塵密度為0.304 kg/m3,與分布范圍內空氣中氧氣當量比為1,顆粒半徑為2 μm,鋁密度為2 700.000 kg/m3. 空氣初始密度為1.210 kg/m3,氣體壓力為0.10 MPa,溫度為298 K. 起爆在中心位置起爆,起爆區為2×2個網格,起爆條件φg=1,ρg=2.200 kg/m3,ug=1 400 m/s,vg=1 400 m/s,Tg=3 200 K. 由于模型中心對稱,為減小計算量,采用計算模型1/4區域,左側與下側采用全反射對稱邊界,上側與右側采用透射邊界。如圖2所示,左下角區域為半徑3.0 m鋁粉塵云團,計算模型尺寸為25 m×25 m,網格數量800×800. 計算中采用的參數[16]:λ=0.1 J/(m·s·K),qa=ηQa,Qa=31.5 kJ/g,η=80%.

圖2 模擬區域尺寸模型Fig.2 Geometry model of simulation area

2.3 數值模擬結果

2.3.1 不同時刻爆炸流場物理量分布

圖3為t=2.13 ms時爆炸流場物理量分布。圖3(a)為爆炸流場壓力分布,此時沖擊波到達3.0 m位置,到達鋁粉塵云團邊緣,由于起爆區域為正方形,因此向上和向右傳播的沖擊波陣面接近平面波,最大壓力為2.10 MPa左右,中間部分沖擊波陣面為圓弧形,壓力為1.75 MPa,波后流場壓力逐漸下降。圖3(b)為爆炸流場密度分布,向上和向右傳播的沖擊波陣面密度較高,約為2.900 kg/m3,中間部分沖擊波波陣面密度為2.600 kg/m3,波后流場密度逐漸下降,中心位置密度較低,為1.000 kg/m3. 圖3(c)為流場溫度分布,可以看出流場中沖擊波陣面處溫度迅速上升,波后流場中大部分區域溫度達到3 800 K. 圖3(d)是流場中鋁粉塵體積分數分布,此時鋁粉塵由于激波作用堆積在波陣面附近,沖擊波陣面上側與右側處體積分數約為5.2×10-5,與初始時刻濃度比為46.2%,右上部體積分數約為2.0×10-5,與初始時刻濃度比為17.8%,中心處鋁粉塵基本反應完畢剩余較少。

圖4為t=3.49 ms時爆轟流場物理量分布,此時沖擊波到達4.7 m位置。圖4(a)是沖擊波流場壓力分布,此時向上和向右傳播的沖擊波陣面壓力下降到0.90 MPa,中間部分沖擊波陣面壓力為0.75 MPa,波后壓力緩慢下降,而后小幅上升達到穩定,中心區域壓力穩定在0.75 MPa. 圖4(b)為爆轟流場密度分布,向上和向右傳播的沖擊波陣面處密度為3.800 kg/m3,中間部分沖擊波陣面處密度為3.200 kg/m3,波后流場密度下降至0.800 kg/m3,中心區域密度約1.300 kg/m3. 圖4(c)為流場溫度分布,沖擊波附近溫度逐漸上升到3 800 K,對比圖4(b)發現此部分區域同時為低密度區域,將此部分高溫低密度區域定義為火球區域,將2 500 K處設為火球邊界,可以看出此時火球與沖擊波開始發生分離。圖4(d)是流場中鋁粉塵體積分數分布,極少量鋁粉塵在沖擊波作用下運動到4.0 m處,可以看出此時鋁粉塵基本反應完畢,只有在沖擊波后上側與右側部分區域存在少量鋁粉塵,濃度與初始時刻比例約0.14%. 此時鋁粉塵速度為480 m/s,沖擊波速度為900 m/s,由于沖擊波波速較鋁粉塵顆粒速度快,因此此時沖擊波到達4.7 m位置而鋁粉塵顆粒只到達4.0 m處。當沖擊波離開云團時,沖擊波陣面與粉塵云團分離,反應區的能量無法使沖擊波保持增長,沖擊波壓力降下降。

圖3 t=2.13 ms時刻的流場壓力、密度、溫度和鋁粉塵體積分數Fig.3 Pressure, density and temperature of flow field, and volume fraction of aluminum particles for t=2.13 ms

圖4 t=3.49 ms時刻的流場壓力、密度、溫度和鋁粉塵體積分數Fig.4 Pressure, density and temperature of flow field, and volume fraction of aluminum particles for t=3.49 ms

圖5為t=10.99 ms時刻流場的壓力、密度和溫度分布圖,可以看出此時沖擊波到達10.0 m處。圖5(a)為流場壓力分布圖,沖擊波陣面壓力衰減到0.36 MPa,波后區域壓力逐漸下降,中心區域壓力約0.17 MPa. 圖5(b)為流場密度分布圖,沖擊波陣面處密度約為2.700 kg/m3,沖擊波后密度逐漸下降,中心區域密度約為0.400 kg/m3,此區域半徑約為6.5 m,對比圖5(c)可以看出此部分區域為火球區域。圖5(c)為流場溫度分布圖,沖擊波后流場溫度上升至800 K左右,在火球處溫度迅速上升,火球中心區域溫度在3 800K左右,火球與沖擊波產生明顯分離。圖5(d)為流場中O2體積分數φO分布,可以看出火球內部中心區域中O2體積分數只有0.02左右,火球內部靠外區域中O2體積分數稍高為0.03,火球外O2體積分數為0.23. 在鋁粉塵爆炸過程中會產生氣態鋁,而在火球擴散過程中氣態鋁與O2發生反應釋放能量,維持火球中心的高溫。

圖5 t=10.99 ms時刻的流場壓力、密度、溫度和氧氣體積分數Fig.5 Pressure, density and temperature of flow field,and volume fraction of O2 for t=10.99 ms

圖6 t=29.32 ms時刻的流場壓力、密度、溫度和氧氣體積分數Fig.6 Pressure, density and temperature of flow field, and volume fraction of O2 for t=29.32 ms

圖6為t=29.32 ms時刻流場物理量分布。圖6(a)為流場壓力分布圖,此時沖擊波傳播至20.0 m處,沖擊波陣面壓力衰減至0.23 MPa,沖擊波后流場壓力逐漸下降,中心區域壓力為0.10 MPa. 圖6(b)為流場密度分布圖,沖擊波陣面密度下降至2.100 kg/m3,波后流場密度逐漸下降,到達中心火球區域時密度約為0.200 kg/m3. 圖6(c)為流場溫度分布圖,沖擊波后流場中溫度約800 K,此時火球傳播至8.0 m位置,火球附近流場溫度由800 K迅速上升至3 000 K,由于火球進一步擴散,火球內部溫度下降到3 500 K以上。圖6(d)為流場中O2體積分數分布,流場中O2含量在火球處急劇下降,從火球外0.23下降至0.02左右,可以看出火球區域為一個低壓、低密度、低氧、高溫區域。

2.3.2 流場中不同位置的物理量隨時間的變化曲線

圖7 數據采集點及不同時刻流場近場壓力和鋁粉塵體積分數曲線Fig.7 Data collection points, and pressure in near field and volume fraction of aluminum distribution lines at different times

選取對角線上的點作為研究對象,圖7(a)為數據采集點示意圖,圖7(b)和圖7(c)為不同時刻近場壓力曲線以及鋁粉塵體積分數曲線。起爆后,初始沖擊波向外傳播,鋁粉塵被點火發生反應,沖擊波壓力在3.0 m處粉塵云團邊界位置時上升達到最大值1.75 MPa,隨沖擊波繼續向外傳播,云團內鋁粉塵被氣體加速向外擴散并消耗減少,波后壓力逐漸下降。隨時間增加,如圖7(c)所示,鋁粉塵云團可加速擴散至3.5 m處,極少量粉塵顆粒可以到達4.0 m處反應完畢。

圖8 不同時刻的流場壓力、密度和溫度曲線Fig.8 Pressure, density and temperature distributions at different times

圖8(a)為不同時刻遠場的流場壓力曲線,可以看出在小于4.0 m時,鋁粉塵反應釋放能量,流場壓力較高,4.0 m以后沖擊波壓力迅速下降。到達28.0 m位置時,沖擊波陣面壓力下降至0.19 MPa. 圖9(a)為沖擊波和火球隨時間傳播距離s曲線,圖9(b)為沖擊波和火球隨時間傳播速度曲線,從圖中可以看出,沖擊波傳播速度隨時間急速下降,45 ms時刻沖擊波傳播速度下降至500 m/s,此時沖擊波傳播至28.0 m處。根據沖擊波對暴露人員的損傷程度[17],沖擊波超壓達到0.07 MPa時就會對人體產生嚴重殺傷,達到0.1 MPa時會造成50%致死。根據圖8(a),沖擊波在24.5 m處時沖擊波陣面壓力下降至0.20 MPa,沖擊波在28.0 m處時波陣面壓力為0.19 MPa,因此可以認為半徑24.5 m范圍內致人死亡,24.5~28.0 m范圍內對人體造成嚴重損傷。

圖9 不同時刻沖擊波與火球的傳播距離和速度Fig.9 Propagation distances and velocities of shock wave and fireball at different times

圖8(b)和圖8(c)為不同時刻流場密度與溫度曲線。起爆后,沖擊波陣面處密度隨沖擊波傳播逐漸上升,而后開始下降;沖擊波波后流場密度逐漸下降,到達火球區域時密度產生急速下降,可以看出火球區域密度逐漸穩定在0.150 kg/m3左右。從不同時刻溫度曲線可以看出隨時間增大,火球擴散速度降低,從圖8(b)可以看出,33 ms以后火球傳播速度基本接近為0 m/s,到45 ms時,火球半徑為10.0 m,火球中心溫度在3 500 K以上,在半徑10.0 m范圍內以高溫及超壓殺傷為主。

圖10 t=43.87 ms時刻沿對角線處流場參數分布Fig.10 Parameters of flow field at the points of diagonal for t=43.87 ms

圖10為43.87 ms時刻沿對角線區域的流場參數曲線。圖10(a)為密度分布曲線,此時沖擊波到達27.0 m位置,沖擊波陣面上密度為1.900 kg/m3,波后流場密度下降,在火球區域密度急速下降至0.120 kg/m3. 圖10(b)為壓力分布曲線,沖擊波陣面處壓力降為0.19 MPa,波后流場壓力下降。在沖擊波后壓縮波,在實驗中也觀察到了此現象[3]。在對圖11中流場壓力分布隨時間變化曲線分析后得出,在沖擊波與火球發生分離時,波后壓力下降,而火球內由于鋁粉塵及氣態鋁反應放熱溫度較高,壓力較高,因此沖擊波后壓力下降后會產生小幅躍升,隨沖擊波傳播距離增大,火球內部鋁反應完畢,火球內中心區域壓力下降,因此在火球邊緣處形成一個壓縮波,而火球傳播速度隨時間增大而急劇下降,造成壓縮波與火球發生脫離,在沖擊波后形成二次波,在火球邊緣附近壓力有小幅躍升而后下降。圖10(c)為流場氣體速度分布曲線,沖擊波陣面處氣體速度165 m/s,波后流場中氣體速度下降,在壓縮波處氣體速度為70 m/s,在火球邊緣處氣體速度由10 m/s躍升至30 m/s,而后氣體速度下降。圖10(d)為流場中溫度曲線,沖擊波陣面后氣體溫度在350 K左右,到達火球處時氣體溫度迅速上升,火球中心溫度達到3 500 K以上。

圖11 不同時刻的局部流場壓力曲線Fig.11 Pressure distributions of local flow field at different times

3 結論

本文數值模擬了當量比為1、密度為0.304 kg/m3、顆粒半徑為2 μm的鋁粉塵在空氣中均勻分布形成的半徑3.0 m的云團,在被中心起爆后對周圍環境的溫壓效應,得出以下結論:

1)初始沖擊波在粉塵云團中傳播并引起鋁粉塵點火發生化學反應形成爆轟波,波陣面壓力在沖擊波到達云團邊界3.0 m處時達到最大值1.75 MPa. 隨沖擊波傳播的距離增加,沖擊波壓力下降, 這是由于沖擊波與云團分離,粉塵中的能量無法支持沖擊波造成的。

2) 鋁粉塵反應形成的火球隨時間可向外擴散到達10.0 m位置,火球內部為密度0.120 kg/m3、溫度3 500 K以上的低密度、低氧區域。

3) 鋁云團爆轟在半徑10.0 m區域內對周邊毀傷以高溫及沖擊波超壓為主,在半徑24.5 m區域內沖擊波超壓可致人死亡,24.5~28.0 m范圍內可對人體造成嚴重殺傷。

本文對空氣中懸浮鋁粉塵云團的爆轟進行了數值模擬,分析了自由場中鋁粉塵云團爆轟對周邊環境的殺傷規律,對模擬溫壓彈的殺傷及鋁粉塵爆炸安全事故對周邊環境的殺傷及防護提高了認識。

)

[1] 裴明敬, 毛根旺, 胡華權, 等, 含鋁溫壓燃料性能研究[J]. 含能材料, 2007, 15(5):442-446.

PEI Ming-jing, MAO Gen-wang, HU Hua-quan, et al. Characteristic of the thermobaric explosive contained aluminum powders[J]. Chinese Journal of Energetic Materials, 2007, 15(5):442-446.(in Chinese)

[2] 鄭波, 陳力, 丁雁生, 等, 溫壓炸藥爆炸拋灑的運動規律[J]. 爆轟與沖擊, 2008, 28(5): 433-437.

ZHENG Bo, CHEN Li, DING Yan-sheng, et al. Dispersal process of explosion production of thermobaric explosive[J]. Explosion and Shock Waves, 2008, 28(5):433-437.(in Chinese)

[3] 李秀麗, 惠君明, 王伯良. 云爆劑爆炸/沖擊波參數研究[J]. 含能材料, 2008, 16(4):410-414.

LI Xiu-li, HUI Jun-ming, WANG Bo-liang. Blast/shock wave parameters of single-event FAE[J]. Chinese Journal of Energetic Materials, 2008, 16(4) :410-414.(in Chinese)

[4] 洪滔, 秦承森. 鋁顆粒激波點火機制初探[J]. 爆轟與沖擊, 2003, 23(4):295-299.

HONG Tao, QIN Cheng-sen. Mechanism of shock wave ignition of aluminum particle[J]. Explosion and Shock Waves, 2003, 23(4):295-299.(in Chinese)

[5] 洪滔, 秦承森. 懸浮鋁粉塵爆轟波參數[J]. 含能材料, 2004, 12(3):129-133.

HONG Tao, QIN Cheng-sen. Parameters of detonation in suspended aluminum dust[J]. Chinese Journal of Energetic Materials, 2004, 12(3):129-133.(in Chinese)

[6] 洪滔, 秦承森. 爆轟波管中鋁粉塵爆轟的數值模擬[J]. 爆轟與沖擊, 2004, 24(3):193-200.

HONG Tao, QIN Cheng-sen. Numerical simulation of dust detonation of aluminum powder in explosive tubes[J]. Explosion and Shock Waves, 2004, 24(3):193-200.(in Chinese)

[7] 昝文濤, 洪滔, 董賀飛. 基于CE/SE方法關于RDX-AL懸浮粉塵在空氣中的兩相爆轟的數值模擬[J]. 爆炸與沖擊, 2016, 36(5):603-610.

ZAN Wen-tao, HONG Tao, DONG He-fei. Numerical simulation of two phase detonation of suspending RDX-AL dust in air with CE/SE method[J]. Explosion and Shock Waves, 2016, 36(5):603-610.(in Chinese)

[8] 昝文濤, 洪滔, 董賀飛. 帶管道連接的空間中懸浮鋁粉塵爆轟波傳播數值模擬[J]. 含能材料, 2017, 25(6): 508-514.

ZAN Wen-tao, HONG Tao, DONG He-fei. Numerical simulation of detonation of suspending aluminum dust in air in the universal spaces connected by channel[J]. Chinese Journal of Energetic Materials, 2017, 25(6):508-514.(in Chinese)

[9] Zan W T, Hong T, Dong H F. Simulation of double-front detonation of suspended mixed RDX and aluminum dust in air[J]. Chinese Physics Letter, 2017, 34(7): 074701.

[10] Steinberg T A, Wilson D B, Benz F. The combustion phase of burning particle[J]. Combustion and Flame, 1992, 91(2):200-208.

[11] Price E W. Combustion of metalized propellants[C]∥Progress in Astronautics and Aeronautics: Fundamenals of Solid-Propellant Combustion. NY, US:AIAA, 1984:479-513.

[12] Chang S C. The method of space-time conservation element and solution element-a new approach for solving the Navier-Stokes and Euler equations[J]. Journal of Computational Physics, 1995,119(2): 295-324.

[13] Zhang D L, Wang J T, Wang G. High-order CE/SE method and applications[J]. Chinese Journal of Computational Physics, 2009, 26(2):211-220.

[14] Wang J T, Zhang D L, Liu K X. A Eulerian approach based on CE/SE method for 2D multimaterial elastic-plastic flows[J]. Chinese Journal of Computational Physics, 2007, 24(4):395-401.

[15] Tulis A J, Selman J R. Detonation tube studies of aluminum particles dispersed in air[J]. Symposium (International) on Combustion, 1982, 19(1):655-663.

[16] 洪滔, 秦承森, 林文洲. 懸浮RDX炸藥和鋁顆粒混合粉塵爆轟的數值模擬[J]. 爆炸與沖擊, 2009, 29(5):468-473.

HONG Tao, QIN Cheng-sen, LIN Wen-zhou. Numerical simulation of detonation in suspended mixed RDX and aluminum dust[J]. Explosion and Shock Waves, 2009, 29(5):468-473.(in Chinese)

[17] GJB 5212—2004 云爆彈定型實驗規程[S]. 北京: 中國人民解放軍總裝備部, 2004.

GJB 5212—2004 Finalizing test procedures for fuel-air-explosive ammunition[S]. Beijing: General Armament Department of PLA, 2004.(in Chinese)

猜你喜歡
區域
分割區域
探尋區域創新的密碼
科學(2020年5期)2020-11-26 08:19:22
基于BM3D的復雜紋理區域圖像去噪
軟件(2020年3期)2020-04-20 01:45:18
小區域、大發展
商周刊(2018年15期)2018-07-27 01:41:20
論“戎”的活動區域
敦煌學輯刊(2018年1期)2018-07-09 05:46:42
區域發展篇
區域經濟
關于四色猜想
分區域
公司治理與技術創新:分區域比較
主站蜘蛛池模板: 欧美国产日韩在线| 欧美性精品不卡在线观看| 亚洲日本精品一区二区| 亚洲第一区在线| 国产在线精品美女观看| 永久在线精品免费视频观看| 青青草久久伊人| 一本大道AV人久久综合| 国产精品爽爽va在线无码观看| 99精品国产自在现线观看| 伊人久综合| 超碰91免费人妻| a级毛片毛片免费观看久潮| 黄色网在线免费观看| 欧美亚洲国产日韩电影在线| 亚洲 成人国产| 国产精品网曝门免费视频| 成人午夜网址| 好久久免费视频高清| 国产午夜精品鲁丝片| 亚洲欧美h| 亚洲成人网在线观看| 亚洲精品在线影院| 婷婷午夜影院| 免费Aⅴ片在线观看蜜芽Tⅴ| 婷婷五月在线| 亚洲男人天堂2020| 澳门av无码| 成色7777精品在线| 国产超碰在线观看| 97免费在线观看视频| 五月天香蕉视频国产亚| 国产真实乱子伦精品视手机观看 | 国产在线观看成人91| 国产chinese男男gay视频网| 久久久91人妻无码精品蜜桃HD| 国产精品思思热在线| 亚洲va精品中文字幕| 日韩在线影院| 精品国产成人av免费| 久久久精品国产亚洲AV日韩| 亚洲国产av无码综合原创国产| 国产在线小视频| 欧美不卡视频在线| 亚洲成AV人手机在线观看网站| 久久国产乱子伦视频无卡顿| 欧美一级专区免费大片| 国产精品视频a| 亚洲综合激情另类专区| 另类综合视频| 久久综合丝袜日本网| 99免费视频观看| 久久精品娱乐亚洲领先| 国产乱人乱偷精品视频a人人澡| 久久亚洲AⅤ无码精品午夜麻豆| 国产综合无码一区二区色蜜蜜| 亚洲AV人人澡人人双人| AⅤ色综合久久天堂AV色综合| 92午夜福利影院一区二区三区| 国产不卡网| 激情成人综合网| 91精品小视频| 啪啪啪亚洲无码| 国产夜色视频| 国产特一级毛片| 国产无码制服丝袜| 欧美激情第一区| 亚洲成人在线网| 亚洲国产91人成在线| 97国产成人无码精品久久久| 国产精品欧美亚洲韩国日本不卡| 亚洲国产欧美国产综合久久| 波多野结衣一区二区三区四区视频| 波多野结衣AV无码久久一区| 精品国产欧美精品v| 亚洲一区二区日韩欧美gif| 日本在线亚洲| 中文无码影院| 国产成人一区在线播放| 精品人妻无码中字系列| 国产视频a| 玖玖精品在线|