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

通氣超空泡多相流場(chǎng)數(shù)值仿真方法

2013-09-19 03:13:26周景軍董春鵬尹韶平項(xiàng)慶睿
關(guān)鍵詞:模型

周景軍, 董春鵬, 尹韶平, 項(xiàng)慶睿

?

通氣超空泡多相流場(chǎng)數(shù)值仿真方法

周景軍, 董春鵬, 尹韶平, 項(xiàng)慶睿

(中國(guó)船舶重工集團(tuán)公司 第705研究所, 陜西 西安, 710075)

通氣超空泡流動(dòng)涉及多相流動(dòng)、湍流、相變及可壓縮等流體力學(xué)難點(diǎn)問(wèn)題, 流動(dòng)機(jī)理非常復(fù)雜。其中多相流模型是通氣超空泡數(shù)值仿真研究工作的重點(diǎn), 將嚴(yán)重影響通氣超空泡數(shù)值仿真結(jié)果的精度。本文有針對(duì)性地對(duì)比了目前廣泛采用的均質(zhì)平衡流模型和歐拉雙流體模型, 結(jié)合作者所在課題組多年來(lái)在水洞試驗(yàn)和數(shù)值仿真方面的研究成果, 從空泡形態(tài)和流體動(dòng)力兩方面分析了歐拉雙流體模型在預(yù)測(cè)通氣超空泡方面的優(yōu)勢(shì)。隨著研究的進(jìn)一步深入, 通氣超空泡數(shù)值仿真方法有望成為超空泡減阻技術(shù)的重要研究手段, 可以為工程設(shè)計(jì)提供參考。

通氣超空泡; 多相流; 數(shù)值仿真方法; 均值平衡流模型; 歐拉雙流體模型

0 引言

對(duì)于通氣超空泡流動(dòng)的研究, 最早可以追溯到上個(gè)世紀(jì)40年代。1946年, Reichardt首次提出通過(guò)人工通氣的方法可以在較低的速度或環(huán)境壓力下形成超空泡, 從而使超空泡試驗(yàn)可以在低速水洞中進(jìn)行, 極大地推動(dòng)了通氣超空泡的發(fā)展。

超空泡流動(dòng)由于涉及到高雷諾數(shù)、多相流、相變及可壓縮性等復(fù)雜流動(dòng)問(wèn)題, 研究和測(cè)量手段非常有限。試驗(yàn)方法無(wú)疑是研究通氣超空泡最為有效的方法, 但是作為超空泡武器主要研究手段的模型試驗(yàn)方法, 受限于試驗(yàn)設(shè)備和測(cè)試條件, 尚無(wú)法開展高速、小空泡數(shù)、大尺度的空泡水洞試驗(yàn)。大尺度模型的有動(dòng)力水下發(fā)射試驗(yàn)雖然可以提高試驗(yàn)速度、降低空泡數(shù), 但研制經(jīng)費(fèi)大、周期長(zhǎng), 通過(guò)少量條次的試驗(yàn)結(jié)果, 難于總結(jié)科學(xué)規(guī)律, 特別是無(wú)法揭示復(fù)雜流動(dòng)現(xiàn)象的結(jié)構(gòu)和機(jī)理。水洞試驗(yàn)仍然是目前最主要的研究手段, 但是由于來(lái)流速度和工作段尺度的限制, 水洞試驗(yàn)主要以機(jī)理性研究為主, 其結(jié)果如何應(yīng)用于大尺度實(shí)航試驗(yàn)國(guó)內(nèi)還未有定論。理論方面最具代表性的有基于勢(shì)流理論速度勢(shì)的面元法、以及由前蘇聯(lián)著名學(xué)者Logvinovich提出的獨(dú)立擴(kuò)張?jiān)?。面元法可以快速得到空泡形態(tài), 但是由于空泡尾部閉合方式需要進(jìn)行假設(shè), 同時(shí)無(wú)法考慮尾部泄氣模式, 主要用于自然空化的仿真??张莳?dú)立擴(kuò)張?jiān)砜梢钥焖佾@得非定??张菪螒B(tài), 對(duì)于研究空泡穩(wěn)定性具有重要作用, 但是對(duì)于航行體流體動(dòng)力尤其是尾部滑行力和尾翼上的流體動(dòng)力, 目前主要還是采用經(jīng)驗(yàn)公式獲得, 該方法嚴(yán)重依賴于空泡形態(tài)和航行體相對(duì)位置決定的浸水深度, 可以用于彈道仿真, 但工程應(yīng)用需要的前提是空泡形態(tài)穩(wěn)定可控、可測(cè), 難度很大。隨著計(jì)算機(jī)能力的迅速提高和計(jì)算流體力學(xué)(Compu- tational Fluid Dynamics, CFD)的快速發(fā)展, CFD技術(shù)已經(jīng)成為航空、航天以及船舶領(lǐng)域有力的研究手段, 可以高效地獲得流場(chǎng)細(xì)節(jié), 給出流體動(dòng)力參數(shù)。對(duì)于通氣超空泡流動(dòng)而言, 由于試驗(yàn)技術(shù)的復(fù)雜性, 求解通氣超空泡多相流場(chǎng)的CFD技術(shù)對(duì)于通氣超空泡減阻技術(shù)的發(fā)展將起到至關(guān)重要的作用, 一方面可以結(jié)合水洞試驗(yàn)驗(yàn)證模型精度, 深入研究通氣超空泡生成機(jī)理; 另一方面數(shù)值仿真方法不受來(lái)流速度、模型尺寸等的限制, 研究對(duì)象不受限制。

目前, 通氣超空泡數(shù)值仿真主要還是延用上世紀(jì)90年代發(fā)展起來(lái)的用于模擬自然空化的均質(zhì)平衡流模型, 國(guó)內(nèi)大量的研究結(jié)果表明, 均質(zhì)平衡流模型在預(yù)測(cè)通氣超空泡方面盡管具有較高的計(jì)算效率, 但在計(jì)算精度方面尤其是通氣量仿真方面和試驗(yàn)結(jié)果相比誤差較大。

均質(zhì)平衡流模型預(yù)測(cè)通氣超空泡的結(jié)果國(guó)外公開成果相對(duì)國(guó)內(nèi)較少, 最具代表性的是美國(guó)賓州大學(xué)Kunz領(lǐng)導(dǎo)的課題組, 他們開發(fā)了一套完整的空泡流數(shù)值計(jì)算仿真程序, 并且可以實(shí)現(xiàn)和航行體6自由度方程進(jìn)行耦合求解。計(jì)算結(jié)果如圖1所示。

因?yàn)闆](méi)有試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比, 計(jì)算精度無(wú)法評(píng)價(jià)。但是從空泡形態(tài)來(lái)看, 均質(zhì)平衡流模型的計(jì)算結(jié)果有明顯的特點(diǎn), 即空泡尾部無(wú)法閉合且一直延伸到計(jì)算流場(chǎng)尾部。

針對(duì)上述問(wèn)題, 本文基于商業(yè)軟件CFX, 首次采用歐拉雙流體模型進(jìn)行通氣超空泡流場(chǎng)的仿真工作, 通過(guò)多組水洞試驗(yàn)驗(yàn)證, 證明所采用的歐拉雙流體模型在預(yù)測(cè)通氣超空泡流動(dòng)方面相對(duì)均質(zhì)平衡流模型更能給出合理的流場(chǎng)細(xì)節(jié), 對(duì)工程設(shè)計(jì)具有一定的參考價(jià)值。

1 空泡流動(dòng)特點(diǎn)

空化問(wèn)題中, 不論是自然空化還是通氣空化首先涉及到的是多相流動(dòng)問(wèn)題。對(duì)于自然空化而言, 空泡的生成是由局部靜壓驅(qū)動(dòng), 當(dāng)局部壓力低于飽和蒸汽壓時(shí), 產(chǎn)生自然空化, 自然空化同時(shí)還受到水中氣核的影響。對(duì)于通氣空化而言, 其過(guò)程與自然空化過(guò)程相反, 通入的不可凝結(jié)氣體排開周圍的水, 使局部壓力升高, 導(dǎo)致空化數(shù)降低生成空泡。自然空化尤其是云狀空化更加符合均質(zhì)平衡流動(dòng)假設(shè)條件, 目前大量研究成果均證明該方法在預(yù)測(cè)自然空化方面具有較高的計(jì)算精度。但對(duì)于大尺度的通氣超空泡而言, 其內(nèi)部為通入的不可凝結(jié)氣體, 空泡外部為充滿水的外流場(chǎng), 整個(gè)流場(chǎng)更加像是分層流動(dòng), 只有在空泡邊界和空泡尾部閉合區(qū)域2種流體才存在一定程度的混合, 大尺度的通氣超空泡已很難滿足均質(zhì)平衡流動(dòng)的假設(shè), 歐拉雙流體模型分別對(duì)每相的動(dòng)量方程進(jìn)行求解, 更為符合通氣超空泡實(shí)際流型, 但是兩相界面上的作用模型對(duì)空泡流動(dòng)仿真精度有較大影響, 是研究的重點(diǎn), 總之, 兩相是否求解一套動(dòng)量方程是2種多相流模型的主要區(qū)別。

2 多相流模型控制方程

2.1 均質(zhì)平衡流模型

均質(zhì)平衡流模型認(rèn)為兩相之間擁有共同的速度場(chǎng)和壓力場(chǎng), 目前采用較多的是混合物模型(mixture model)和流體體積(volume of fluid, VOF) 函數(shù)模型2種。Mixture模型認(rèn)為單元體積內(nèi)擁有統(tǒng)一的速度場(chǎng)、壓力場(chǎng), 兩相之間不存在相對(duì)速度, 流場(chǎng)密度由兩相的體積分?jǐn)?shù)和密度共同決定, 從控制方程來(lái)看屬于單流體模型, 只求解一套質(zhì)量、動(dòng)量和能量守恒方程。VOF模型采用在固定的歐拉網(wǎng)格下的表面跟蹤方法, 該方法適合于兩相或多相不相混合的流動(dòng), 此時(shí)并不關(guān)心某相的流動(dòng), 而是相之間界面的運(yùn)動(dòng)規(guī)律, 在該模型中, 只求解一套動(dòng)量方程, 全流場(chǎng)的每個(gè)計(jì)算單元內(nèi), 分別得到各流體組分所占有的體積率, 通過(guò)求解混合介質(zhì)的動(dòng)量方程和處理穿過(guò)區(qū)域的各流體的體積分率來(lái)模擬2種或多種不能混合的流體。

1) 連續(xù)性方程

(2)

2) 動(dòng)量方程

(4)

3) 體積分?jǐn)?shù)方程

4) 體積分?jǐn)?shù)守恒方程

考慮氣體的可壓縮性時(shí), 還需求解理想氣體的狀態(tài)方程。

2.2 歐拉雙流體模型

歐拉雙流體模型分別求解每相的連續(xù)性方程、動(dòng)量方程和能量方程。兩相之間的作用通過(guò)相界面進(jìn)行傳遞, 通過(guò)不同的模型來(lái)考慮。

1) 連續(xù)性方程

由于空氣的溶解度較低, 認(rèn)為水和空氣之間不發(fā)生質(zhì)量輸運(yùn), 水-氣兩相作用時(shí)的連續(xù)性方程為

空化時(shí)的連續(xù)性方程為

(9)

2) 動(dòng)量方程

水-氣兩相作用時(shí)的動(dòng)量方程為

空化時(shí)的動(dòng)量方程為

(12)

3) 體積守恒方程

(13)

如果不考慮可壓縮性以及忽略源相, 得

3 數(shù)值仿真結(jié)果與分析

2種多相流模型預(yù)測(cè)通氣超空泡的有效性將通過(guò)下面2組水洞試驗(yàn)結(jié)果進(jìn)行檢驗(yàn)。

3.1 后支撐模型仿真和試驗(yàn)?zāi)P蛯?duì)比

由圖2可知, 從空泡形態(tài)來(lái)看, 當(dāng)通氣量較小時(shí), 空泡內(nèi)部充滿氣水混合物, 兩相混合較為均勻, 2種多相流模型仿真結(jié)果差別不大, 說(shuō)明均質(zhì)平衡流假設(shè)此時(shí)成立。由圖3可以看出, 隨著通氣量的增加, 局部空泡逐漸發(fā)展成超空泡, 2種多相流模型預(yù)測(cè)結(jié)果有明顯區(qū)別。由圖4可以看出, 歐拉雙流體模型預(yù)測(cè)結(jié)果和試驗(yàn)結(jié)果非常接近, 均質(zhì)平衡流模型生成同樣長(zhǎng)的空泡需要的通氣量大約為試驗(yàn)的5倍。

(1) VOF模型???? (2) 歐拉雙流體模型

(a) 通氣量0.000 01 kg/s仿真結(jié)果

(a)VOF模型

(b)歐拉雙流體模型

圖3 通氣量0.000 1kg/s空泡形態(tài)仿真結(jié)果

Fig. 3 Simulation result of the cavity shape with air mass flow of 0.000 1kg/s

(a)試驗(yàn)結(jié)果

(b)歐拉雙流體模型仿真結(jié)果

圖4 歐拉雙流體模型和試驗(yàn)結(jié)果空泡形態(tài)對(duì)比

Fig. 4 Cavity shape comparison between Euler two- fluid model and experiment

圖5為2種多相流模型空泡內(nèi)部流動(dòng)計(jì)算結(jié)果??梢钥闯? 均質(zhì)平衡流模型計(jì)算結(jié)果內(nèi)部分離不明顯, 接近單相流動(dòng), 通入的不可凝結(jié)氣體被水迅速帶走。歐拉雙流體模型內(nèi)部氣體分離現(xiàn)象明顯, 內(nèi)部氣體流動(dòng)除了壓力之外受外部影響較小, 符合實(shí)際流動(dòng)情況, 預(yù)測(cè)結(jié)果更為合理。

(a) 均質(zhì)平衡流模型內(nèi)部流場(chǎng)速度分布

(b) 歐拉雙流體模型內(nèi)部流場(chǎng)速度分布

圖5 不同多相流模型內(nèi)部流場(chǎng)計(jì)算結(jié)果

Fig. 5 Calculation results of the inner flow field with different multiphase models

3.2 前支撐模型仿真和試驗(yàn)?zāi)P蛯?duì)比

為了進(jìn)一步說(shuō)明歐拉雙流體模型在計(jì)算通氣超空泡方面的計(jì)算精度, 對(duì)水洞試驗(yàn)前支撐試驗(yàn)進(jìn)行了仿真, 仿真結(jié)果如圖6所示, 可以看出, 歐拉雙流體模型可以很好地得到空泡由于重力效應(yīng)導(dǎo)致的上漂現(xiàn)象, 空泡尾部閉合狀態(tài)與試驗(yàn)符合較好, 生成同樣尺度的空泡, 需要的通氣量與試驗(yàn)相差在15%以內(nèi), 具有較高的計(jì)算精度。

(a) 仿真結(jié)果

試驗(yàn)結(jié)果

另外, 本次試驗(yàn)采用六分力天平同時(shí)測(cè)出了空泡發(fā)展過(guò)程中航行體尾部滑行力的變化規(guī)律。模擬結(jié)果如圖7所示, 可以看出, 歐拉雙流體模型和試驗(yàn)結(jié)果兩者變化規(guī)律符合較好, 但均質(zhì)平衡流模型只能得到尾部升力的定性變化規(guī)律。

4 結(jié)論

本文通過(guò)試驗(yàn)結(jié)果對(duì)2種多相流模型在預(yù)測(cè)通氣超空泡方面進(jìn)行了評(píng)價(jià), 具體結(jié)論如下。

1) 均質(zhì)平衡流模型在預(yù)測(cè)自然空化、通氣初生空化方面計(jì)算精度可以滿足工程需要, 計(jì)算量相對(duì)歐拉雙流體模型較小。研究結(jié)果同時(shí)表明, 對(duì)于自然空化和通氣初始階段, 空泡兩相之間由于湍流的作用混合較為均勻。

2) 在預(yù)測(cè)大尺度通氣空泡流方面, 歐拉雙流體模型在預(yù)測(cè)通氣空泡形態(tài)和流體動(dòng)力方面都具有較高的精度。均質(zhì)平衡流模型只能進(jìn)行定性預(yù)報(bào), 無(wú)法模擬出兩相分離狀態(tài), 空泡壁面和尺度難以進(jìn)行精確捕捉, 尤其是通氣量和試驗(yàn)相比誤差較大。

3)隨著通氣超空泡數(shù)值仿真方法研究的進(jìn)一步深入, 可以和水洞試驗(yàn)有效互補(bǔ), 為通氣超空泡航行體流體動(dòng)力特性研究提供重要研究手段。

通氣超空泡流動(dòng)數(shù)值仿真方法的精度主要與多相流模型有關(guān), 空泡尾部泄氣量、航行體尾部沾濕面流體動(dòng)力等還與湍流模型有直接關(guān)系。當(dāng)航行體進(jìn)入高速巡航段時(shí)還涉及自然空化水-汽-氣三相流動(dòng)問(wèn)題。總之, 通氣超空泡流動(dòng)非常復(fù)雜, 通氣超空泡流動(dòng)數(shù)值仿真工作還有大量問(wèn)題需要解決, 通過(guò)試驗(yàn)研究增加對(duì)通氣超空泡流動(dòng)機(jī)理的認(rèn)識(shí)是進(jìn)一步完善通氣超空泡流動(dòng)數(shù)值仿真方法的有效途徑。

[1] Reichardt H. The Laws If Cavitation Bubbles as Axially Symmetrical Bodies in a Flow[R]. Ministry of Aircraft Productuin Reports and Translations. Britian, 1946: 322-326.

[2] Martin W, Travis J S, Roger E A A. Experimental Study of a Ventilated Supercavitating Vehicle[C]//5th International Symposium on Cavitation, Osaka, Japan, 2003: 11-13

[3] 蔣潔明, 魯傳敬, 胡天群, 等. 軸對(duì)稱體通氣空泡的水動(dòng)力試驗(yàn)研究[J]. 力學(xué)季刊, 2004, 25(4): 450-456.Jiang Jie-ming, Lu Chuan-jing, Hu Tian-qun, et al. The Experimental Research on the Cavitating Flow Round An Asisymmetric Body with Ventilation[J]. Chinese Quarterly of Mechanics, 2004, 25(4): 450-456.

[4] 蔣增輝, 于開平, 張嘉鐘, 等. 超空泡航行體尾部流體動(dòng)力特性試驗(yàn)?zāi)P椭畏绞降倪x擇研究[J]. 機(jī)械科學(xué)與技術(shù), 2007, 26(12): 1648-1651. Jiang Zeng-hui, Yu Kai-ping, Zhang Jia-zhong, et al. Selecting the Strutting Position for Models of Experiments on Hydrodynamic Characteristics of the Aft Section of an Underwater Supercavitating Vehicle[J]. Mechanical Science and Technology for Aerospace Engineering, 2007, 26(12): 1648-1651.

[5] Logvinovich G V. Hydrodynamics of Free-Boundary Flow[R]. U. S. Department of Commerce, Washington, DC, 1972, Trans- lated from the Russian (NASA-TT-F-658).

[6] Qin Q, Song C S, Arndt R E. A Virtual Single-phase Natural Cavitation Model and Its Application to Cav2003 Hydrofoil [C]//5th International Symposium on Cavitation, Osaka, Japan, 2003.

[7] Leroux J B, Oliver C D, Astolfi J A. A Joint Experimental and Numerical Study of Mechanisms Associated to Instability of Partial Cavitation on Two-dimensional Hydrofoil[C]//5th International Symposium on Cavitaion, Osaka, Japan, 2003.

[8] Oliver C D, Astolfi J A. Numerical Prediction of the Cavitation Flow on a Two-dimensional Symmetrical Hydrofoil with a Single Fluid Model[C]//5th International Symposium on Cavitation, Osaka, Japan, 2003.

[9] Kubota A, Kato H, Yamaguchi H. A New Modelling of Cavitating Flows: a Numerical Study of Unsteady Cavitation on a Hydrofoil Section[J]. J Fluid Mech, 1992, 240: 59-96.

[10] Kunz R F, Boger D A, Stinebring D R, et al. A Proconditioned Navier-Stokes Method for Two Phase Flows with Application to Cavitation[J]. Computers and Fluids, 2000, 29: 849-875.

[11] Kunz R F, Boger D A, Chyczeweki T S, et al. Multi-phase CFD Analysis of Natural and Ventilated Cavitation about Submerged Bodies[C]//Third ASME/JSME Joint Fluids Engi- neering Conference. San Francisco, California, 1999.

[12] Jules W. Lindau, Robert F. Kunz, Jason M. Mulherin, James J. Dreyer, David R. Stinebring. Fully Coupled 6-Dof to Urans Modeling of Cavitating Flows Around a Supercavitating Vehicle[C]//Fifth International Symposium on Cavitation (CAV2003)Osaka, Japan, 2003.

[13] 郭建紅, 魯傳敬, 陳瑛, 等. 基于輸運(yùn)方程類空化模型的通氣空泡流數(shù)值模擬[J]. 力學(xué)季刊, 2009, 30(3): 378-384.

[14] Singhal A K, Vaidya N, Leonard A D. Multi-dimensional Simulation of Cavitating Flows Using a PDF Model for Phase Change[C]//ASME Fluids Engineering Division Summer Meeting, 1997.

(責(zé)任編輯: 陳 曦)

Numerical Simulation Method for Ventilated Supercavitating Multiphase Flow Field

ZHOU Jing-jun, DONG Chun-peng, YIN Shao-ping, XIANG Qing-rui

(The 705 Research Institute, China Shipbuilding Industry Corporation, Xi′an 710072, China)

Ventilated supercavitating flow involves such topics in fluid mechanics as multiphase flow, turbulence, phase change and compressibility, its mechanism is very complex. The multiphase flow model has attracted much more attention in the study of numerical simulation of supercavitating flow, however its accuracy in simulation is not satisfactory. In this paper, the homogeneous model, which are widely used in the world, are compared with the Euler two-fluid model by combining with the authors′ research by means of water tunnel experiments and numerical simulation. The advantages of the Euler two-fluid model in predicting ventilated supercavitation is analyzed in terms of cavity shape and hydrodynamics of a vehicle. Numerical simulation of ventilated supercavitation is expected to become an important approach of drag-reduction technology through supercavitation.

ventilated supercavitation; multiphase flow; numerical simulation method; homogeneous model; Euler two-fluid model

TJ630.1; O351.2

A

1673-1948(2013)03-0165-06

2012-10-12;

2012-11-13.

周景軍(1981-), 男, 博士, 主要研究方向?yàn)轸~雷總體技術(shù).

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产麻豆va精品视频| 91精品国产福利| 亚洲成在线观看| 成人在线综合| 亚洲中文字幕23页在线| 中文字幕乱码二三区免费| 久久精品aⅴ无码中文字幕 | 国产麻豆aⅴ精品无码| 国产欧美另类| 国产美女无遮挡免费视频| 日韩中文欧美| 国产在线观看人成激情视频| 免费99精品国产自在现线| 日本欧美中文字幕精品亚洲| 综合久久五月天| 久久精品亚洲专区| 色综合网址| 免费高清a毛片| 无码一区18禁| 一本一道波多野结衣一区二区 | 久久综合亚洲色一区二区三区| 亚洲免费黄色网| 国产一级在线播放| 国产AV毛片| 精品人妻系列无码专区久久| 少妇精品在线| 综合色亚洲| 亚洲激情99| 久久九九热视频| 欧美中文字幕第一页线路一| 色悠久久综合| 国产亚洲男人的天堂在线观看| 免费日韩在线视频| 久久国产精品麻豆系列| 18禁色诱爆乳网站| 亚洲午夜综合网| 九九香蕉视频| 青青草91视频| 久久成人18免费| 尤物特级无码毛片免费| 亚洲欧美一区二区三区图片 | 亚洲国产清纯| 人妻出轨无码中文一区二区| 国产精品性| 日韩高清欧美| 最新亚洲人成网站在线观看| 无码一区二区三区视频在线播放| 18禁黄无遮挡网站| 亚洲三级片在线看| 播五月综合| 国产欧美一区二区三区视频在线观看| 欧美19综合中文字幕| 四虎精品国产永久在线观看| 超级碰免费视频91| 99热这里只有成人精品国产| 91精品国产一区自在线拍| 国产精品七七在线播放| 国产不卡国语在线| 国产成年女人特黄特色大片免费| 亚洲a免费| 欧美综合成人| 国产无码高清视频不卡| 精品人妻无码中字系列| 国产av剧情无码精品色午夜| 欧美成人精品在线| 日韩欧美中文字幕在线韩免费| 青青热久麻豆精品视频在线观看| 日韩国产黄色网站| 免费国产高清视频| 欧美中文字幕第一页线路一| 色综合成人| 老司机久久99久久精品播放| 精品久久蜜桃| 日韩欧美国产中文| 潮喷在线无码白浆| 欧美福利在线观看| 久久婷婷六月| 亚洲色图综合在线| 日韩精品久久无码中文字幕色欲| 国产精品网曝门免费视频| 亚洲 成人国产| 国产靠逼视频|