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

水下連發超空泡射彈的流動與阻力特性研究

2020-11-04 01:44:36施紅輝周東輝賈會霞
空氣動力學學報 2020年4期
關鍵詞:模型

施紅輝,周東輝,周 棟,賈會霞

(浙江理工大學 機械與自動控制學院,杭州 310018)

0 引 言

水下連發射彈的超空泡流動問題,來源于水下槍炮的連續發射、超空泡射彈的飽和攻擊等海戰背景,具有很強的軍事應用價值。超空泡射彈是一種新型水下高速武器[1],相比傳統的水下作戰武器,它顯著地提高了自身的運動速度和射程。超空泡射彈的一個重要用途是作為防御武器裝備在水面艦船或者潛艇上[2],通過連續發射高速且密集的射彈形成彈幕保護艦船或者潛艇免受魚雷的威脅。當超空泡武器串聯運行時,這就需要考慮2個甚至多個超空泡的相互作用的問題,由于連發超空泡射彈周圍流場相互影響、干擾,引起射彈間流場的變化,與水下單發射彈的超空泡有較大的區別。它既包含單發射彈超空泡的特性,又有自身特殊的復雜性。因此,研究水下連發射彈的超空泡演化規律與阻力特性在學術研究和工程應用上都是必要的,這方面的研究還很少見。

俄羅斯和烏克蘭的科學家在射彈空泡特性和水動力學方面進行了大量的基礎研究[3-7]。Truscott[8]利用高速攝像技術研究了彈體以小攻角入水時彈體表面誘導產生的超空泡流動態特性,并討論了長徑比與彈體頭型對入水時產生超空泡的形態影響。Yamashita[9]等對水下超高速的球形與細長體彈體進行了研究,發現超空泡的存在大大降低了水下航行體的阻力。Weiland 和Vlachos[10]通過實驗研究了圓柱形射彈高速入水空泡的產生和發展,并指出射彈發射時驅動氣體的泄漏促進了超空泡的快速形成。Lee等[11]利用能量守恒原理對射彈高速入水空泡生成、發展過程及閉合特性進行了研究,建立了高速入水條件下的入水空泡動力學模型。國內學者對水下航行體空化問題也進行了一些研究。易文俊[12]等利用Fluent軟件對水下高速航行體自然超空化流動進行了數值模擬,獲得了航行體的超空泡形態變化特性。金大橋[13]等基于均質平衡多相流理論,利用Fluent軟件對水下射彈自然超空泡減阻特性進行了數值模擬,研究了空化數對水下射彈空泡閉合部位和阻力系數的影響,重點分析了水下射彈結構參數對自然超空泡減阻特性的影響,得到了空化器直徑、模型長細比和不同尾部形狀對水下射彈超空泡減阻特性影響的規律。郭子濤和張偉等[14-15]對幾種柱形彈體水平入水進行了實驗和理論研究,建立了關于空泡擴展的理論模型,獲得了彈體阻力系數和平頭、球頭以及幾種卵形彈體頭型系數之間的關系。魏應杰等[16]對超空泡水下航行體的控制問題進行研究,建立了超空泡航行體非線性動力學模型,設計了非線性切換控制策略。施連會等[17]對超空泡航行體的動力穩定性問題進行了數值研究,計算了超空泡航行體的動力不穩定區域和相應的臨界頻率曲線。魯傳敬課題組[18-19]基于雷諾平均Navier-Stokes方程和混合多相流模型,對通氣空泡與自然空泡共同作用下的三維軸對稱的通氣空泡流動進行了數值模擬,分析了通氣量對水動力學的特性影響,給出適當的通氣量利于超空泡穩定的結論。施紅輝[20-21]課題組利用自行研制的水平超空泡發射系統對水下單發射彈進行了大量試驗,觀察到了多種工況下的超空泡形態,研究了水下射彈的水動力學特性以及水深和彈體長徑比對超空泡彈體阻力系數及空泡形狀影響。

從超空泡問題的提出至今,國內外已經有了許多研究成果,但它們基本上都是針對單個水下物體或射彈的。在水下發射導彈或射彈時,武器能否連續發射直接關系到裝備的應用和作戰效果,這就需要考慮2個甚至多個超空泡的相互作用的問題,而這方面的研究還很少見。何春濤等[22]對在重力加速下低速兩射彈串聯傾斜入水問題進行了初步試驗,發現當兩射彈串聯間距小到一定程度時,后面的射彈完全進入到前面射彈形成的空泡內部。施紅輝等[23]對連發射彈以100 m/s的速度在水下做水平勻速運動誘導的超空泡流場進行了數值模擬,研究發現連發射彈經過初生、發展和相互作用后,能夠形成一個包裹兩射彈的超空泡,超空泡的尾部產生高速回射流。

針對水下連發射彈自然減速運動過程中的超空泡演化規律及運動特性,本文以Flunet流體仿真軟件為平臺,利用動網格技術分別建立了水下兩連發、三連發射彈的二維數值計算模型,通過數值計算分別獲得了兩連發射彈和三連發射彈在水下自由運動過程中的超空泡演化規律、連發射彈頭部前方的壓力分布、連發射彈的位移變化曲線、連發射彈的速度衰減曲線、連發射彈的阻力特性曲線,本文的研究結果率先揭示了水下兩連發和三連發射彈誘導的超空泡流動的典型特征及相關力學參數。

1 控制方程及數值方法

1.1 控制方程

數值計算采用流體體積函數(VOF)多相流模型描述水和水蒸汽構成的多相流動系統。VOF多相流模型將水、水蒸汽兩相當做單一介質的混合流動系統,各相共享一套動量方程,通過計算得到單元內各相的體積分數確定流動系統中各相的分布。

混合相的連續性方程和動量方程分別為:

式中:t為時間;ui、uj分別為笛卡爾坐標系中的速度分量;ρm為混合相密度;p 為流場的壓力;μm為混合相動力黏度;SM為附加的源相。

ρm,μm的表達式分別為:

式中:αv為水蒸汽相的體積分數,ρv、ρl 分別為水蒸汽和水的密度;μv、μl分別為水蒸汽和水的動力黏度。

數值計算采用RNG k-ε 湍流模型[24],由于空化現象是流線強烈彎曲導致的,該模型能夠更好地處理高彎曲流線和高應變率的流動。湍動能k 和耗散率ε 的控制方程分別為:

式中:μeff=μm+μt,μt 為湍流黏度;k 為湍動能,ε 為耗散率,αk、αε分別為k 和ε 的負向效應的普朗特數;Gk為平均速度梯度產生的湍動能;C1ε、C2ε為湍流動能耗散率的經驗常數。

1.2 空化模型

空化模型選擇Schnerr-Sauer空化模型[25],該空化模型具有較高收斂速度和計算穩定性。水蒸汽相方程的一般形式為:

式中:vv為水蒸汽相的速度矢量;Re、Rc分別為蒸發速率和冷凝速率;?B為氣核的半徑;Pv為水的飽和蒸汽壓力。

1.3 數值方法和邊界條件

連發射彈的計算模型為軸向串聯的兩發射彈和軸向串聯的三發射彈,定義相鄰兩射彈的間距為H,射彈的物理模型是直徑D=6 mm,長度L=48 mm的圓柱體,如圖1所示,來源于參考文獻[21]的試驗模型,材質為5005鋁鎂合金,質量為3.56 g。由于數值模擬中射彈設定的速度很高,計算求得的弗勞德數Fr?30,根據重力效應的判斷準則,可以不考慮重力的影響。由于射彈為軸對稱回轉體,數值計算采用二維軸對稱模型。整個網格區域寬度是200 mm,長度是1000 mm,足夠避免邊界效應。數值模擬中采用射彈運動、流體靜止的方法研究連發射彈的自然超空泡,這種方法更接近真實的環境。通過用戶自定義函數(UDF)編寫程序嵌入到Flunet中模擬具有初動能的連發射彈自由運動,在運動過程中網格的更新和生成采用動態層法實現。分別對水下兩連發和三連發射彈進行了數值模擬,連發射彈彼此的間距H =10D=60 mm,這兩種情況的計算模型如圖2所示。圖2(a)表示兩連發射彈的計算域及邊界條件設置,圖2(b)表示三連發射彈的計算域及邊界條件設置,這兩種計算模型的計算域邊界條件劃分相同,右側邊界條件為壓力出口邊界條件,左側邊界條件為無滑移壁面條件,計算域上側設定為無滑移壁面條件,射彈設置為無滑移壁面條件,射彈從計算域右側運動到左側,射彈的初始發射速度都設定為相同的速度v0=200 m/s,之后射彈做自由運動。整個計算域采用四邊形結構網格劃分,射彈附近的網格進行加密處理。

圖1 射彈的物理模型Fig.1 Physical model of projectile

圖2 連發射彈的計算域和邊界條件設置Fig.2 Computational region and boundary conditions of successively fired projectiles

對于上述計算模型壓力與速度之間的耦合求解采用PISO 算法,壓力場的空間離散采用PRESTO!格式,各項體積率離散采用CICSAM 格式,密度和動量采用二階迎風離散格式,空間離散采用二階迎風格式,時間離散采用一階隱格式。

2 數值方法驗證

在研究水下連發射彈超空泡流動前,先進行數值方法的驗證。由于水下連發超空泡射彈試驗難度很大,并且缺乏相關的文獻數據,所以選取水下單發超空泡射彈試驗進行數值模擬驗證。采用上文所述的數值模擬方法對參考文獻[21]中的工況8進行了數值計算,工況8采用的射彈和圖1所示的射彈模型相同,質量為3.56 g,試驗中射彈的初速度大小為83.1 m/s,水深為90 mm。數值計算的整個網格區域寬度是200 mm,長度是1000 mm,右側邊界條件為壓力出口邊界條件,左側邊界條件為無滑移壁面條件,網格為四邊形網格,射彈周圍的網格進行加密處理,加密區的網格尺寸為0.2 mm,時間步長取4×10-7s。

在水下深度一定的情況下,射彈高速無推力自由運動過程中,忽略重力和浮力的影響,射彈形成的超空泡內外壓力差Δp 相同,空化數僅由射彈速度vp決定,空化數σ 的表達式[26]如下:

式中:h 為水的深度,p0為一個標準大氣壓,pc為水的飽和蒸汽壓,這里取值為3540 Pa。

超空泡狀態下圓柱體的阻力系數可以近似表達為[23]:

式中:Cd0是空化數σ=0的阻力系數[26](圓盤空化器的Cd0≈0.82)。

對于水下單發射彈無推力自由運動,由牛頓第二定律可得:

式中:m 是射彈質量,S 是射彈位移,A0是射彈的截面積。

定義速度衰減系數kv=ρlA0Cd0/2m ,通過積分方程(12),可得:

關于單發超空泡射彈的位移隨時間變化,圖3顯示了數值模擬結果、實驗結果、公式(13)的理論計算結果三者之間的對比,縱坐標表示的是無量綱位移(S/D),從圖中可以看出數值模擬的射彈位移變化趨勢與實驗和理論計算結果基本一致,在數值上吻合良好。圖4比較了數值模擬和實驗在t=2 ms時獲得的射彈水下空泡形狀,從圖中可以看出兩者的空泡形狀基本一致。經過上述的對比,有效地驗證了數值模擬方法的準確性。

圖4 數值模擬和實驗在t=2 ms時的超空泡形狀Fig.4 Comparisons of supercavity configurations between numerical simulation and experiments at t=2 ms

以水下兩連發射彈為例進行網格無關性驗證,建立了3 中不同密度的網格,網格數分別為45.5 萬(case1,網格最小尺寸為0.4 mm)、53.7萬(case2,網格最小尺寸為0.2 mm)、65.8萬(case3,網格最小尺寸為0.1 mm),取時間步長為4×10-7s,分別采用上述三種網格密度對發射速度為200 m/s的兩連發射彈入水進行了數值模擬,得到了前發射彈速度隨時間的變化曲線,如圖5所示。從圖中可知,隨著網格密度的增加,case2的網格密度和case3的網格密度計算所得的前發射彈運動速度之間的差異已經小到可以忽略不計,綜合考慮計算精度和計算成本,選取了網格密度為case2的網格用于本文計算。

在不影響計算正常的前提下,選取了三種時間步長值(d t1=4×10-7s,d t2=2×10-7s,d t3=1×10-7s),針對不同的時間步長對兩連發射彈入水的計算結果進行了分析驗證,獲取了前發射彈速度隨時間的變化曲線,如圖6所示。通過對比可知,三種時間步長下前發射彈的速度變化規律一致,這說明此范圍內的時間步長對于本文的數值模擬是適用的,考慮時間成本的限制,選擇了時間步長d t1=4×10-7s進行數值模擬。

圖5 不同網格密度的計算結果Fig.5 Calculated results of different mesh densities

圖6 不同時間步長的計算結果Fig.6 Calculated results of different time steps

3 結果與討論

水下連發射彈誘導的自然超空泡的發展是一個非定常過程。水下兩連發和三連發射彈的超空泡演化過程如圖7所示,圖7中紅色代表水相,藍色代表水蒸汽相。圖7(a)表示的是水下兩連發射彈超空泡演化過程,對其進行分析,可以將兩連發射彈超空泡的發展過程分成3個階段:第1階段,射彈1和射彈2各自形成獨立的超空泡,射彈被包裹在空泡內部,只有射彈頭部沾濕,射彈受到的阻力主要為壓差阻力。第2階段,兩射彈運動過程中,在射彈1超空泡流場的影響下,射彈2的頭部逐漸靠近并最終進入射彈1超空泡的尾部,兩發射彈的超空泡發生部分融合,隨后第二個超空泡的大部分被分離(圖中對應t=1.0~1.2 ms),分離后的空泡在周圍水壓的作用下慢慢變小并最終發生潰滅,射彈2進入射彈1的超空泡內部運動后并逐漸追上射彈1。t=1.5 ms時,射彈2已經完全在射彈1超空泡的內部運動。第3階段,射彈2與射彈1發生追尾碰撞,碰撞后兩者共同運動,文獻[22]在開展的兩個平頭圓柱體串聯低速入水試驗中有觀察此種現象,這也證明了數值模擬結果的準確性。t=2.6 ms時,射彈2追上射彈1并發生追尾,在兩射彈發生追尾后,超空泡壁面上出現一個波動現象,射彈1頭部流動分離點附近的空泡壁面先出現波谷(圖中t=2.6 ms時,①處的放大圖中橢圓標記處),之后發生追尾處的空泡壁面出現波峰(圖中t=2.7 ms時,②處的放大圖中橢圓標記處),推測此現象產生的原因是由兩射彈發生追尾碰撞產生的瞬態擾動引起的;該波動現象一直存在,隨著時間的推移,存在波動的空泡壁面出現頸縮,空泡從此位置逐漸開始分裂。圖7(b)給出了水下三連發射彈超空泡演化過程,其發展過程可以分為4個階段:第1階段,三發射彈各自形成獨立的超空泡。第2階段,射彈3進入射彈2的超空泡,射彈2進入射彈1形成的超空泡。第3階段,射彈3從射彈2脫離的超空泡內部穿出,穿出時其頭部與水接觸再次生成空泡,生成的空泡與射彈2脫離的空泡相互融合;第4階段,射彈2先追趕上射彈1并發生追尾碰撞,碰撞后兩者共同運動,隨后射彈3追趕上射彈2并發生追尾碰撞,碰撞后三者共同運動。t=2.9 ms時,射彈2和射彈1發生追尾碰撞,射彈3 部分進入射彈1 的超空泡內。t=4.9 ms時,射彈3與射彈2發生追尾碰撞。射彈發生碰撞后同樣在超空泡壁面上出現擾動現象,與兩連發射彈的情形一樣,這里不再詳細敘述。隨著時間的推移,三連發射彈存在擾動的空泡壁面出現頸縮,并且空泡從此位置開始分裂。

圖7 連發射彈水下運動誘導的超空泡演化過程Fig.7 Evolution of supercavities induced by successively fired projectiles

圖8 表示的是兩連發射彈在時間t=0.7 ms時的壓力云圖以及射彈頭部中心點前4D 處的壓力分布,圖8(b)中的橫坐標為無量綱位移S1/D,其中S1表示點到射彈頭部中心點的水平位移,縱坐標表示無量綱壓力p/p0。從圖中可知,高壓區出現在前、后發射彈的頭部表面,前發射彈頭部中心點處的壓力達到85倍的大氣壓,后發射彈頭部中心點處的壓力為57倍的大氣壓,在射彈頭部前方流域壓力值迅速降低,距離平頭頂點越遠,壓力越小;低壓區出現在兩發射彈的空泡內,由于空化效應其壓力值和飽和蒸汽壓力保持一致(約為3540 Pa)。此時刻前發射彈超空泡的尾部與后發射彈的頭部相距的距離很小,從而引起兩發射彈周圍的流場互相干擾,具體表現為前發射彈超空泡尾部空泡在后發射彈頭部附近的高壓區作用下部分潰滅,后發射彈頭部中心點處的壓力及頭部前方的高壓區范圍小于前發射彈的原因是前發射彈超空泡尾部存在回射流[19],導致后發射彈相對前方流體的速度減小,而且前發射彈形成的空泡射流為汽水混合物,密度遠小于水的密度。

圖9給出了三連發射彈在時間t=0.7 ms時的壓力云圖以及射彈頭部中心點前4D 處的壓力分布。從圖中可以看出射彈1頭部駐點處的壓力最大,其大小約為85倍的大氣壓,而射彈3和射彈2頭部前方的壓力值比較接近,這說明后發射彈都會受到前一發射彈所形成的超空泡流場的影響。

圖8 t=0.7 ms,兩連發射彈的壓力分布圖Fig.8 Pressure distribution of two successively fired projectiles at t=0.7 ms

圖9 t=0.7 ms,三連發射彈的壓力分布圖Fig.9 Pressure distribution of three successively fired projectiles at t=0.7 ms

圖10 給出了連發射彈速度隨時間變化的曲線。圖10(a)表示的是兩連發射彈速度變化,可以看出射彈1和射彈2在運動初期的速度衰減趨勢相同,并且速度衰減很快,這是由于兩射彈的速度開始很高,受到水的阻力大;在0.7 ms<t<2.6 ms,射彈2的速度衰減逐漸變慢到速度保持不變,射彈1的速度繼續衰減,這是因為射彈2逐漸接近射彈1尾部空泡時,受到前發射彈尾部空泡低壓區的影響,受到的壓差阻力減小,當射彈2進入到射彈1的空泡中時,射彈2周圍的介質變成了水蒸氣相,密度遠小于水,導致射彈2受到很小的阻力,這里我們不妨標記為0,而射彈1繼續受到水的阻力作用;在時間t=2.6 ms,兩射彈的速度發生躍變,射彈1的速度突然變大,射彈2的速度突然變小,這是因為射彈1和射彈2發生了追尾碰撞。碰撞后,射彈1、射彈2的速度衰減相同。圖10(b)表示的是三連發射彈速度隨時間的變化。在射彈運動初期,射彈1、射彈2、射彈3 速度的衰減相同。在1 ms<t<2.1 ms,射彈3開始進入射彈2的超空泡中,所以此時間段射彈3的速度基本保持不變;在2.1 ms≤t<2.5 ms,射彈3從前一發射彈脫離的超空泡中穿出,受到水的阻力其速度又開始減小;2.5 ms≤t<4.9 ms時,射彈3進入了射彈1的超空泡中導致其所受阻力幾乎為0,從而速度保持不變。射彈2 在1.2 ms≤t<2.9 ms進入了射彈1的超空泡中引起它的速度保持不變,t=2.9 ms時,射彈1和射彈2發生追尾碰撞引起兩射彈的速度發生躍變;2.9 ms≤t<4.9 ms,射彈1和射彈2速度變化相同。t=4.9 ms時,射彈3與射彈2發生碰撞,射彈1、射彈2、射彈3三者的速度發生躍變;t>4.9 ms,射彈1、射彈2、射彈3速度衰減相同。

圖10 連發射彈速度衰減曲線Fig.10 Velocity attenuation of successively fired projectiles

圖11 顯示的是連發射彈位移隨時間變化的曲線,其中縱坐標為無量綱位移S/D,坐標原點為初始時刻射彈的重心位置,位移的正方向為射彈運動的速度方向。圖11(a)表示的是兩連發射彈的無量綱位移變化,從圖中可以看出隨著時間增加,射彈2的位移相比射彈1的位移越來越大,證明兩射彈之間的距離在逐漸縮短,在t=2.6 ms,射彈2的位移恰好比射彈1的位移大10D,說明射彈2在此時追上射彈1,之后兩射彈的位移差始終為10D,表明追尾后射彈1和射彈2共同運動。圖11(b)給出了三連發射彈無量綱位移的變化,射彈2與射彈1的位移差隨時間逐漸增加,射彈3與射彈2的位移差隨時間逐漸增加,在t=2.9 ms時,射彈2 的位移恰好比射彈1 的大10D,表明此時射彈2追上射彈1,t=4.9 ms時,射彈3的位移恰好比射彈2的大10D ,表明此時射彈3追上射彈2。

圖11 連發射彈位移隨時間變化曲線Fig.11 Displacement of successively fired projectiles with time

射彈在水中運動時,會受到摩擦阻力和壓差阻力的影響,這兩部分阻力的和為射彈受到的總的阻力。對該阻力進行無量綱化,則得到無量綱阻力系數。無量綱阻力系數[23]的定義如下:

式中:Fd為射彈阻力。

根據式(14),可得到如圖12所示的連發射彈的阻力系數變化曲線圖。圖12(a)表示的是兩連發射彈的阻力系數隨時間(射彈1和射彈2發生碰撞前)的變化,在運動的起始階段,射彈1和射彈2剛與水接觸,射彈要帶動周圍的流體質點一起運動,受到水的阻力很大,所以此階段射彈1、射彈2的阻力系數較大;隨著射彈1與射彈2的超空泡尾部逐漸接近,射彈2受到射彈1超空泡流場的影響,從而引起射彈2頭部附近的壓力降低,導致射彈2受到的壓差阻力減小,對射彈2而言,射彈1相當于起到導流的作用,因此射彈2的阻力系數逐漸降低;隨著射彈2進入到前一發射彈的空泡中運動,此時射彈2的阻力系數降低到幾乎為0,而射彈1繼續受到水的阻力,它阻力系數在某一定值附近震蕩變化。

圖12(b)表示的是三連發射彈的阻力系數隨時間(射彈1和射彈2碰撞前)的變化曲線。射彈1的阻力系數在起始時刻較大,之后逐漸減小到某一定值附近振蕩變化;射彈2的阻力系數隨時間增加逐漸減小到幾乎為零,這是由于射彈2在逐漸接近射彈1的過程中,受到射彈1的超空泡流場的影響所致;射彈3的阻力系數隨時間增加先減小到趨近于0,這是因為當射彈3進入射彈2形成的超空泡內部時所受阻力很小,在時間t=2.1 ms時射彈3從射彈2的超空泡內部穿出,其頭部與水突然接觸受到一個巨大的瞬時作用力,導致射彈3的阻力系數突然增大到一個很高的值,之后射彈3進入射彈1形成的超空泡內,射彈3的阻力系數又減小到接近于0。

圖12 連發射彈阻力系數隨時間變化曲線Fig.12 Drag coefficient curves of successively fired projectiles

4 結 論

1)水下連發射彈各自形成的超空泡流場互相影響,表現為前發射彈超空泡的流場能夠減小后發射彈頭部前方的高壓區范圍和駐點處壓力,后發射彈的超空泡與前發射彈的超空泡在相互作用后,后發射彈發生脫離自身形成的超空泡進入前發射彈超空泡內部的行為。

2)后發射彈進入到前發射彈的超空泡中時,后發射彈所受阻力幾乎為0,而前發射彈頭部沾濕繼續受到水的阻力作用,導致后發射彈逐漸追趕上前發射彈并發生追尾。

3)前、后發射彈發生追尾后在前發射彈的頭部流動分離點和追尾處,超空泡壁面出現擾動,存在擾動的空泡壁面會發生頸縮,空泡逐漸從此位置分裂。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 午夜在线不卡| 99久视频| 国产精品55夜色66夜色| 人妻少妇乱子伦精品无码专区毛片| 久久这里只精品国产99热8| 久草中文网| 中文毛片无遮挡播放免费| 热伊人99re久久精品最新地| 国产精品嫩草影院视频| 永久免费无码日韩视频| 免费在线国产一区二区三区精品 | 国产第八页| 亚洲欧美日韩高清综合678| 第九色区aⅴ天堂久久香| 中文字幕无码电影| 欧洲高清无码在线| 国产精品人莉莉成在线播放| 国内精品免费| 久久香蕉国产线看观看亚洲片| 高清欧美性猛交XXXX黑人猛交| 久久美女精品国产精品亚洲| 国产无遮挡猛进猛出免费软件| 狠狠躁天天躁夜夜躁婷婷| h视频在线观看网站| 亚洲美女一区| 国产精品熟女亚洲AV麻豆| 色综合天天娱乐综合网| 欧美性爱精品一区二区三区| 最新国产麻豆aⅴ精品无| 日韩AV无码免费一二三区| 91人人妻人人做人人爽男同| 波多野结衣一二三| 99色亚洲国产精品11p| 色婷婷成人| 国产精品毛片在线直播完整版| 一级毛片在线免费看| 久久亚洲国产视频| 91丝袜乱伦| 国产亚洲精品97AA片在线播放| 在线播放国产99re| 国产欧美自拍视频| 四虎成人在线视频| 人妻无码中文字幕一区二区三区| 91福利在线观看视频| 亚洲欧美精品在线| 国产91小视频在线观看| 综合色天天| 最新国产午夜精品视频成人| 亚洲人妖在线| 国产成人久视频免费| 毛片一区二区在线看| 国产网站免费观看| 精品国产电影久久九九| 成人字幕网视频在线观看| 久久香蕉国产线看观看亚洲片| 四虎永久在线视频| 一级一级特黄女人精品毛片| 亚洲无码高清视频在线观看| 91精品国产麻豆国产自产在线| 免费观看男人免费桶女人视频| 国产一级二级三级毛片| 亚洲人成影院午夜网站| 午夜三级在线| 在线观看国产小视频| 免费人成在线观看视频色| 亚洲成年人片| 亚洲bt欧美bt精品| 久久9966精品国产免费| 久久综合亚洲鲁鲁九月天| 国产精品无码AⅤ在线观看播放| 亚洲精品动漫| 国产精品欧美亚洲韩国日本不卡| 亚洲视频四区| 精品国产污污免费网站| 无码日韩精品91超碰| 日本三级欧美三级| 黄色污网站在线观看| 亚洲中文字幕97久久精品少妇| 婷婷综合色| 无套av在线| 国产香蕉97碰碰视频VA碰碰看| 亚洲另类国产欧美一区二区|