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

雙發回轉體水下齊射流體動力特性數值仿真

2021-11-08 06:51:20瑤,山,
水下無人系統學報 2021年5期
關鍵詞:模型

施 瑤, 高 山, 潘 光

雙發回轉體水下齊射流體動力特性數值仿真

施 瑤, 高 山, 潘 光

(1. 西北工業大學 航海學院, 陜西 西安, 710072; 2. 西北工業大學 無人水下運載技術重點實驗室, 陜西 西安, 710072)

為研究雙發回轉體水下齊射過程流體動力演化特性, 基于均質多相流理論、標準RNG-模型、Singhal空化模型以及重疊網格技術, 建立三維雙發回轉體水下齊射模型, 開展不同發射速度下雙發回轉體非定常空泡與運動姿態演變過程數值仿真, 分析了典型工況下雙發回轉體齊射過程流場結構演變、運動特性以及齊射速度對其運動特性的影響。結果表明: 航行前期空泡發展至最大, 隨著航行器向自由液面運動,空泡從其末端由下至上逐漸脫落, 并發生了潰滅現象; 由于齊射過程中流動干擾區域的存在, 雙發回轉體肩部空泡形態演變過程由不對稱演變為對稱狀, 從而導致回轉體質心先向內側偏移, 隨后不斷向外側偏移;隨著齊射速度的增大, 回轉體出筒時刻空泡長度不斷增大, 其質心由內側向外側偏移的交點向后推遲, 同時偏轉角度不斷減小。此研究對水下齊射工程技術具有一定的借鑒意義。

雙發回轉體; 水下齊射; 空泡; 回射流; 水動力特性

0 引言

齊射是指依托水下移動發射平臺將兩發或多發航行器以一定發射空間間距在極短時間內連續彈射出水的過程。單筒多航行器水下齊射技術因具有儲彈量多、攻擊扇面大、易于實現飽和攻擊等優勢, 可極大提高航行器突防概率, 因此受到各軍事強國的高度關注。然而, 在水下齊射過程中, 航行器肩部低壓區會出現空化, 并伴隨著生長、脫落以及潰滅等復雜非定常多相流現象, 加之齊射過程中, 航行器之間會產生較為嚴重的流動干擾現象, 導致其兩側非定常空化區域不對稱演變, 直接影響其運動姿態的穩定性。

水下發射過程常常涉及大尺度空泡群非定常演化, 早期研究的單空泡相關理論成果難以直接應用[1-3]。因此, 國內外相關學者針對回轉體水下發射過程非定常空泡特性與運動特性開展了相關研究。Dyment等[4]基于雷諾平均方程(Reynolds averaged Navier-Stokes, RANS)方法和流體體積(volume of fluid, VOF)多相流模型對回轉體出水過程的尾空泡進行數值仿真, 仿真結果與相關實驗吻合度良好; 王一偉等[5-6]針對回轉體水下發射過程中空泡生成、發展、脫落和潰滅等開展了系統研究, 給出了將回射流運動時間與回轉體運動時間比值作為空泡穩定性的重要判據; 魏英杰等[7]將質量輸運空化模型與混合介質RANS方法相結合, 給出了回轉體阻力系數和空化數之間的關系, 分析了肩空泡對回轉體流體動力特性的影響。權曉波等[8]采用MIXTURE多相流模型與動網格技術對回轉體水下發射過程進行了計算, 獲得了尾空泡生成演化的周期性特征, 同時發現了發射水深增加導致尾空泡壓力增大的現象。

回轉體水下齊射涉及多體間流動干擾的水動力問題, 由于此概念較新, 相關文獻較少, 但可以借鑒多體在無固壁干擾條件下相互干擾的相關研究成果[9-10]。何春濤[11]對2個并聯圓柱體入水進行了數值仿真, 分析了空泡內部輪廓的演變過程, 給出了雙體入水過程中流動干擾消失的臨界發射間距。宋武超等[12]基于勢流理論和非線性假設, 引入二維軸對稱入水空泡計算模型和影響函數, 給出了雙體并聯入水過程空泡的三維演化特性。Xu等[13]研究了雙發回轉體以不同時序出水過程, 發現當反向旋轉渦對出現時, 對次發回轉體的運動姿態產生較大的影響。

綜上可知, 關于水下垂直發射過程中非定常空泡演化與出水彈道特性的研究主要集中在單發航行器, 而針對航行器水下齊射過程, 應重點關注流動干擾條件下航行器肩部空泡非定常演化與出水彈道特性, 但目前針對該問題的研究有限。因此, 文中通過開展三維雙發回轉體水下齊射過程水動力特性數值仿真, 獲取非定常空泡脫落及潰滅過程中回轉體運動特性的演變規律, 此研究對后續水下齊射工程技術具有一定的借鑒意義。

1 數值計算方法

1.1 控制方程

描述回轉體水下垂直發射氣液多相流動的基本控制方程包括連續性方程、動量方程和能量方程, 其基本形式如下。

連續性方程

動量方程

能量方程

1.2 湍流方程

采用標準RNG-模型, 通過修正湍流黏度, 并考慮了平均流動中的旋轉和旋流流動情況, 能夠更好地處理高應變率以及流線彎曲程度較大的流動。

其中, 湍流黏性

1.3 空化模型

采用Singhal空化模型, 其計算公式如下

上式包含的假設有:

1) 不考慮韋伯數的變化, 將其視為經驗系數;

2) 速度與相變率呈線性關系;

當地壓力小于空化壓力, 即<p時, 空化率為

式中, 經驗參數取值為C=0.02,C=0.02。

1.4 VOF多相流模型

VOF多相流模型是一種在固定Euler網格下的表面跟蹤方法, 它適用于在能夠解決混合物各相之間界面的數值網格上模擬多種互不相容流體的流動。文中求解的氣-液-汽3相流動問題中, 液相為主相, 其余2相為次相。如圖1所示, 若液相體積分數為, 氣相體積分數為, 汽相體積分數為, 則混合相體積分數為1---。

1.5 重疊網格技術

重疊網格的節點分為洞內點、計算點和插值點3種。洞內點不參與流場計算, 計算點參與流體計算, 插值點進行流場信息的傳遞, 3種網格節點在重疊網格分布如圖2所示。

圖1 VOF模型原理示意圖

圖2 重疊網格示意圖

重疊網格技術的基本思想是利用子域網格在重疊區域進行插值處理實現流場信息的實時傳遞, 解決了傳統貼體網格重組過程帶來的網格畸變等問題。重疊網格的實現流程包括網格生成、網格裝配、數值計算、網格更新直至求解完成全過程, 具體的實現流程如圖3所示。

2 數值計算模型

2.1 幾何模型與網格劃分

首先建立雙發回轉體水下齊射模型, 其中模型直徑=15 mm, 長細比/=10,為回轉體長度, 采用半球頭型。考慮到文中研究雙發回轉體水下齊射流動干擾過程, 為了控制計算量和提高計算效率, 截取1/2雙發回轉體水下齊射流場計算域開展數值計算, 如圖4所示。

圖5為回轉體模型網格劃分, 為了更好捕捉回轉體周圍流場細節演變, 對發射筒口、筒口上方流域以及自由液面進行了局部加密, 其中第1層網格高度為0.045 mm, 所對應的+為40; 為更好地節約計算資源, 在遠場設置較大網格尺度。

圖5 模型網格劃分

2.2 數值算法與網格無關性驗證

為了驗證數值仿真方法的有效性, 基于上述數值計算方法的建立及邊界條件設置, 文中自主設計了回轉體水下發射試驗裝置, 并對單發回轉體水下垂直發射過程開展實驗與數值計算對比分析。圖6所示為回轉體水下發射實驗平臺示意圖, 主要由發射系統、控制系統, 高速攝像系統以及防護回收系統組成。發射裝置由儲氣瓶、減速電機、水下導軌、發射筒及各種線路等組成, 通過控制系統調節儲氣瓶氣壓大小來使航行器達到所需的出筒速度; 高速攝像系統由控制電腦、Phantom型號高速攝像機及相應線路等組成。調節高速攝像機的白平衡、分辨率(640×1 024)、幀率(3000幀/s)及曝光時間等參數, 以拍到清晰的畫面; 實驗中選用2盞1 000 W和2盞500 W新聞燈作為背景光源, 以保證攝像視野清晰。防護回收系統由泡沫防護板、起吊裝置等組成。實驗中采用泡沫防護板來抵消航行器出水慣性載荷的沖擊, 既保證了實驗人員的安全性, 同時又防止航行器頭型損壞。此外, 仿真和實驗縮比模型為1:1, 初始條件設置一致。

圖7和圖8分別給出了仿真計算與試驗相圖對比和豎直方向位移對比。從圖中可以發現, 水下運動階段, 仿真結果和試驗結果數據吻合度較好, 包括尾空泡演變以及數值方向彈道。試驗中尾空泡由于發射筒口不均勻氣團效應導致回射流產生, 與仿真存在一定的誤差, 豎直方向彈道在出水階段最大有7.2%的誤差。誤差來源主要是試驗數據提取誤差以及出水噴濺造成的誤差, 可以認為此仿真方法達到了要求的精度。

圖6 水下垂直發射試驗系統示意圖

圖7 仿真計算與實驗相圖對比

圖8 豎直方向位移對比

針對上述回轉體模型, 開展了粗糙尺度網格(2.37×106)、中等尺度網格(4.05×106)與精細尺度網格(6.13×106)計算結果對比, 不同尺度網格計算下回轉體豎直方向位移演變如圖9所示。

圖9 網格無關性驗證

從圖9中可以看出粗糙尺度網格計算結果與中等、精細尺度網格差異較大, 而中等尺度網格與精細尺度網格計算結果基本一致。考慮到計算成本和效率, 故選取中等尺度網格(4.05×106)開展數值計算, 滿足重疊網格計算要求。

3 仿真結果與分析

3.1 流場結構演變

圖10為典型工況下速度為18 m/s時雙發回轉體水下齊射液相體積分數演變云圖。由圖可知,當回轉體尾端完全出筒時, 其速度達到最大值, 此時空化數最小, 回轉體肩部低壓區空泡發展最長。后期隨著回轉體不斷向自由液面運動, 速度不斷減小, 空化數增大, 空泡沿著其末端向回轉體頭部方向不斷脫落直至消失。圖11為雙發回轉體水下齊射壓力演變云圖。由圖可知, 雙發回轉體附近的壓力分布基本一致, 其中回轉體對流場壓力的影響主要集中在頭部和尾部駐點附近的高壓區, 其中尾端駐點附近高壓區形成主要原因如下: 當周圍流體沿回轉體頭部區域運動至尾端區域時, 由于慣性作用, 運動流體在失去回轉體幾何表面的約束后, 會呈一定角度向回轉體尾端中心線位置靠攏, 此角度與回轉體幾何形狀以及速度密切相關; 當兩股流體在尾端匯合時, 速度方向一致的流體分子將會相互疊加, 速度方向相反的流體分子將會發生激烈碰撞, 從而造成尾端駐點高壓區的形成。

圖10 雙發回轉體水下齊射液相體積分數演變云圖

圖11 雙發回轉體水下齊射壓力演變云圖

圖12和圖13分別給出了典型時刻回轉體空泡壓力與形態演變。從圖中可以發現, 雙發回轉體初期空泡末端出現了明顯的不對稱逆壓梯度, 導致空泡末端形態發生不對稱演變。由于空泡從其末端開始潰滅, 以致其末端近似駐點位置壓力較高, 加之回轉體之間流動干擾區域的影響, 空泡末端兩側駐點位置與空泡內部低壓區共同形成較高的逆壓梯度, 此逆壓梯度與航行器軸線之間的夾角小于90°, 形成不對稱逆壓梯度。同時在不對稱逆壓梯度的影響下, 產生不對稱回射流, 如圖14所示為8 ms流線圖。

圖12 典型時刻回轉體空泡壓力演變

圖13 空泡形態演變圖

圖14 8 ms流線圖

隨著航行器向自由液面運動, 不對稱回射流沿著空泡末端兩側不斷向其頭部方向運動, 并伴隨著空泡脫落及潰滅現象; 整個過程中, 隨著航行器運動姿態的演變, 回射流由不對稱發展成對稱狀, 最后空泡發生潰滅現象, 對回轉體運動姿態的穩定性產生了較大影響。

3.2 回轉體運動特性

定義回轉體在水中航行初始狀態為0時刻, 偏轉角逆時針為正。如圖15所示為雙發回轉體在軸方向上的速度和位移曲線, 可以發現雙發回轉體運動特性基本一致, 在運動過程中雙發回轉體質心首先向內側偏轉約0.6 mm, 隨后不斷向外側偏轉約2 mm。根據速度曲線可知, 雙發回轉體質心后續相互遠離趨勢增大。

圖15 X軸方向速度與位移曲線

航行器水下航行階段, 局部當地空化數在非定常空泡的演變過程中, 可以表示為

航行器垂直向自由液面運動過程中, 由于空化數與速度的二次方成反比, 如圖15所示, 隨著航行器速度減小, 空化數急劇增大, 導致航行器末端空泡發生脫落及潰滅現象。基于此, 針對齊射過程航行器肩部空化現象展開數值分析。

圖16和17分別為回轉體在軸方向上的受力以及沿軸方向的偏轉力矩曲線。從圖中可知, 在水中航行前期, 左側回轉體所受合力為負, 指向其內側; 右側回轉體所受合力為正, 也指向其內側。而左側回轉體偏轉力矩為正, 沿逆時針方向; 右側回轉體偏轉力矩為負, 沿順時針方向, 故而回轉體此過程中所受的合力作用點在質心以下區域。此外, 雙發回轉體受到向外側的合力逐漸減小, 導致其橫向速度增大趨勢有所減小。但由于此時偏轉力矩不斷增大, 回轉體向外側偏轉的角度不斷增大, 伴隨著回轉體出水, 頭部迎流作用消失, 偏轉力矩消失, 偏轉角度此時達到最大, 如圖18所示。

圖16 X軸方向受力曲線

圖17 Y軸方向力矩曲線

圖18 偏轉角速度和偏轉角曲線

結合回轉體空泡形態具體說明, 圖19為典型時刻回轉體空泡形態與受力圖, 由于空泡的不對稱性, 導致回轉體質心偏上區域外側沾濕面積較大, 此時均布力所產生的均布力矩M如圖所示, 加之合力作用于質心以下, 回轉體在此力矩的作用下向外側偏轉, 這與回轉體偏轉角演變曲線一致。隨著雙發回轉體都不斷向外側偏轉, 外側區域演變為背流側, 內側區域演變為迎流側。因此, 相比回轉體背流側區域, 其迎流側空泡脫落速度較快, 導致空泡兩側由不對稱向對稱狀轉變, 直至發生潰滅現象, 如圖20所示。

3.3 發射速度對運動特性影響分析

圖21與圖22分別為不同發射速度(14、16、18 m/s)下回轉體水下齊射彈道及偏轉角曲線。由圖21可知, 隨著發射速度的增加, 雙發回轉體質心向外側偏移點不斷推遲, 偏轉角度不斷減小。隨著回轉體尾端出筒速度的增加, 空化數減小, 空泡發展更長。因此, 齊射速度18 m/s時空泡長度最大, 其空泡潰滅時間相比齊射速度14m/s和16m/s推遲, 導致水中航行階段橫向受力方向由內側向外側轉變交點推遲。

圖19 雙發回轉體空泡演化與受力圖

圖20 不對稱空泡脫落示意圖

圖21 雙發回轉體水下齊射彈道曲線

圖22 雙發回轉體偏轉角度曲線

此外, 隨著發射速度的增加, 回轉體在相同水深條件下的航行時間越短, 加之偏轉力矩方向不變, 導致偏轉角度不斷減小, 如圖22所示。因此基于上述分析, 可通過流動控制手段, 如被動排氣, 來改變航行器肩部附近壓力分布特性, 或采用抗空化頭型來抑制空化現象的產生, 從而提高航行器運動姿態的穩定性。

4 結論

基于均質多相流理論、標準RNG-模型、Singhal空化模型以及重疊網格技術, 建立了三維雙發回轉體水下齊射數值計算模型, 模擬了不同齊射速度下回轉體運動特性演變過程, 獲得以下結論。

1) 回轉體出筒時刻, 其發射速度達到最大值, 空化數最小, 導致空泡長度最大; 隨著回轉體不斷向自由液面運動, 其肩部空泡在回射流的作用下從其末端由下至上逐漸脫落, 并發生潰滅現象。

2) 水中航行階段, 由于齊射過程狹長流動干擾區域的存在, 雙發回轉體肩部空泡形態從不對稱狀演變為對稱狀, 導致雙發回轉體橫向受力在航行初期方向發生變化, 偏轉力矩方向不變, 從而其質心先向內偏移隨后向外偏移。

3) 回轉體齊射速度越大, 出筒時刻空泡長度越長, 其質心由內側向外側偏移的交點越向后推遲; 而回轉體所受偏轉力矩方向不變, 導致偏轉角度不斷減小。

[1] Plesset M S, Chapman R B. Collapse of an Initially Spherical Vapour Cavity in the Neighbourhood of a Solid Boundary[J]. Journal of Fluid Mechanics, 1971, 47(2): 283-290.

[2] Plesset M S, Prosperetti A. Bubble Dynamics and Cavitation[J]. Annual Review of Fluid Mechanics, 1977, 9(1): 145-185.

[3] Brennen C E. Cavitation and Bubble Dynamics[M]. Oxford, UK: Oxford University Press, 1995.

[4] Dyment A, Flodrops J P, Paquet J B, et al. Gaseous Cavity at the Base of an Underwater Projectile[J]. Aerospace ence & Technology, 1998, 2(8): 489-504.

[5] 王一偉, 黃晨光, 吳小翠, 等. 航行體水下垂直發射空泡脫落條件研究[J]. 工程力學, 2015, 32(11): 33-39.

Wang Yi-wei, Huang Chen-guang, Wu Xiao-cui, et al. Investigation of Cavities Shedding Condition on Underwater Vehicles in the Vertical Launch Process[J]. Engineering Mechanics, 2015, 32(4): 544-550.

[6] 王一偉, 黃晨光, 杜特專, 等. 航行體垂直出水載荷與空泡潰滅機理分析[J]. 力學學報, 2012, 44(1): 39-48.

Wang Yi-wei, Huang Chen-guang, Du Te-zhuan, et al. Mechanism Analysis about Cavitation collapse Load of Underwater Vehicles in a Vertical Launching Process[J]. Chinese Journal of Theoretical and Applied Mechanics. 2012, 44(1): 39-48.

[7] 魏英杰, 閔景新, 王聰, 等. 潛射導彈垂直發射過程空化特性研究[J]. 工程力學, 2009, 26(7): 251-256.

Wei Ying-jie, Min Jing-xin, Wang Cong, et al. Research on Cavitation of Vertical Launch Submarine Missile[J]. Engineering Mechanics, 2009, 26(7): 251-256.

[8] 權曉波, 燕國軍, 李巖, 等. 水下航行體垂直發射尾空泡生成演化過程三維數值研究[J]. 船舶力學, 2014(7): 739-745.

Quan Xiao-bo, Yan Guo-jun, Li Yan, et al. Three-dimensional Numerical Study on the Evolution Process of Tail Bubble of Underwater Vehicle Vertical Launching[J]. Journal of Ship Mechanics, 2014(7): 739-745.

[9] 程麗, 張亮, 吳德銘, 等. 無升力雙體水動力干擾計算[J]. 哈爾濱工程大學學報, 2005, 26(1): 1-6.

Cheng Li, Zhang Liang, Wu De-ming, et al. Hydrodynamic Interactions between Two Underwater Non-lifting Bodies[J]. Journal of Harbin Engineering University, 2005, 26(1): 1-6.

[10] 金大橋, 王聰, 魏英杰, 等. 水下軸向串列雙圓柱體帶空泡繞流研究[J]. 哈爾濱工程大學學報, 2010, 31(10): 1329-1334.

Jin Da-qiao, Wang Cong, Wei Ying-jie, et al. Cavitating Flow Study of an Underwater Two Axial Circular Cylinder in Tandem Arrangement[J]. Journal of Harbin Engineering University, 2010, 31(10): 1329-1334.

[11] 何春濤. 典型運動體入水過程多相流動特性研究[D]. 哈爾濱: 哈爾濱工業大學, 2012.

[12] 宋武超, 魏英杰, 路麗睿, 等. 基于勢流理論的回轉體并聯入水雙空泡演化動力學研究[J]. 物理學報, 2018, 67(22): 240-256.

Song Wu-chao, Wei Ying-jie, Lu Li-rui, et al. Dynamic Characteristics of Parallel Water-entry Cavity Based on Potential Flow Theory[J]. Acta Physica Sinica, 2018, 67(22): 240-256.

[13] Xu H, Wei Y, Wang C, et al. On Wake Vortex Encounter of Axial-symmetric Projectiles Launched Successively Underwater[J]. Ocean Engineering, 2019, 189: 106382.1- 106382.11.

Numerical Simulation of Hydrodynamic Characteristics of Double-Revolving Bodies in Underwater Salvo

SHI Yao, GAO Shan, PAN Guang

(1. School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072, China; 2. Key Laboratory of Unmanned Underwater Vehicle Technology of Ministry of Industry and Information Technology, Xi’an 710072, China)

To study hydrodynamic characteristics in the underwater salvo process of double-revolving bodies, a three-dimensional underwater salvo model is built in this study based on the homogeneous multiphase flow theory, standard RNG-model, Singhal cavitation model, and overlapping mesh technique. Numerical simulations of the evolution of an unsteady cavity are conducted to determine themovement attitudes of double-revolving bodies at different launching velocities. The flow structure evolution, motion characteristics, and salvo velocity during the salvo process under a typical condition are analyzed. The results show that the maximum development of the cavity is observed during the early stage of water navigation. As the revolving bodies move towards the free surface, the cavity gradually sheds from its end towards its top. The structure collapses owing to a flow interference region in the salvo process, because the evolution of the cavity in the shoulder of the double-revolving bodies is from asymmetric to symmetric, causing mass center deflection of the revolving bodies from the inside to the outside. As the salvo velocity increases, the length of the cavity increases at the outlet of the tube moment, the intersection point of the deflection of the center of mass from the inside to the outside is delayed backward, and the angle of deflection decreases.

double-revolving bodies; underwater salvo; cavity; hydrodynamic characteristic

施瑤, 高山, 潘光. 雙發回轉體水下齊射流體動力特性數值仿真[J]. 水下無人系統學報, 2021, 29(5): 524-532.

TJ630.1; O351.2

A

2096-3920(2021)05-0524-09

10.11993/j.issn.2096-3920.2021.05.003

2020-11-05;

2020-12-09.

國家自然科學基金(52172324); 中央高校基本科研業務費(3102019JC006).

施 瑤(1988-), 男, 博士, 副研究員, 主要研究方向為水下發射水動力學.

(責任編輯: 許 妍)

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲精品无码久久毛片波多野吉| 欧美日韩在线第一页| 91在线免费公开视频| 国产精品视频系列专区| 久久精品中文字幕少妇| 波多野结衣一区二区三区88| 99青青青精品视频在线| 日本午夜影院| 欧美三级视频网站| 69国产精品视频免费| 亚洲天堂伊人| 欧美日韩导航| 欧美中文字幕在线二区| 日韩在线影院| 免费毛片全部不收费的| 久久精品国产亚洲麻豆| 人妻中文久热无码丝袜| 日本伊人色综合网| 日本91视频| 99久久精品国产自免费| 国产成人毛片| 国产高清毛片| 久久免费视频播放| 欧美三级日韩三级| 国产亚洲精品资源在线26u| 国产午夜精品一区二区三| 成人免费网站久久久| 亚洲综合18p| 久久精品无码专区免费| 曰AV在线无码| 不卡无码h在线观看| 高清无码一本到东京热| 亚洲天堂网2014| 国产人免费人成免费视频| 国产精品3p视频| 亚洲国产精品一区二区高清无码久久| 色综合久久久久8天国| 亚洲成a人片7777| 男人天堂亚洲天堂| 中文字幕欧美成人免费| 亚洲精品久综合蜜| 在线播放精品一区二区啪视频| 玩两个丰满老熟女久久网| 在线亚洲小视频| 成人一级黄色毛片| 久久久久久高潮白浆| 久久精品视频亚洲| 爽爽影院十八禁在线观看| 中文毛片无遮挡播放免费| 色婷婷在线播放| 国产成人1024精品| 第一页亚洲| 九九热精品免费视频| 精品视频免费在线| 亚洲大尺码专区影院| 动漫精品啪啪一区二区三区| 成人va亚洲va欧美天堂| AV网站中文| 不卡午夜视频| 久久精品无码专区免费| 欧日韩在线不卡视频| 99热国产这里只有精品无卡顿"| 欧美在线导航| 欧美成人午夜视频| 在线看片国产| 2021国产精品自产拍在线| 亚洲天堂网视频| 国产资源免费观看| 在线视频亚洲色图| 69综合网| 国产亚洲欧美在线中文bt天堂| 亚洲 欧美 偷自乱 图片 | 欧美a在线看| 五月婷婷中文字幕| 国产精品无码一区二区桃花视频| 欧美激情,国产精品| 午夜欧美理论2019理论| 日日噜噜夜夜狠狠视频| 91高清在线视频| 中文字幕在线欧美| 国产高颜值露脸在线观看| 日韩精品专区免费无码aⅴ|