• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    繞多孔孔板通氣氣體與液體兩相橫射流旋渦特性分析

    2017-08-16 08:12:48劉濤濤王國(guó)玉張耐民黃彪
    兵工學(xué)報(bào) 2017年7期
    關(guān)鍵詞:旋渦空泡空化

    劉濤濤, 王國(guó)玉, 張耐民, 黃彪

    (1.北京理工大學(xué) 機(jī)械與車(chē)輛學(xué)院, 北京 100081;2.北京宇航系統(tǒng)工程研究所, 北京 100076)

    繞多孔孔板通氣氣體與液體兩相橫射流旋渦特性分析

    劉濤濤1, 王國(guó)玉1, 張耐民2, 黃彪1

    (1.北京理工大學(xué) 機(jī)械與車(chē)輛學(xué)院, 北京 100081;2.北京宇航系統(tǒng)工程研究所, 北京 100076)

    為了研究通氣氣體與液體兩相流旋渦特性,采用RNGk-ε湍流模型并結(jié)合Level Set界面捕捉方法,對(duì)繞多孔孔板氣體與液體兩相橫射流流動(dòng)特性進(jìn)行了數(shù)值模擬,并與實(shí)驗(yàn)結(jié)果進(jìn)行了比較。研究結(jié)果表明:液相橫流受到射流氣體的阻礙作用在孔口上游、形成分離鞍點(diǎn)和馬蹄渦,此分離鞍點(diǎn)隨距壁面高度的增加逐漸靠近孔心,形成分離線(xiàn);液相橫流繞過(guò)射流氣體后形成兩個(gè)較為封閉的分離旋渦,此分離旋渦隨距壁面高度的增加逐漸遠(yuǎn)離孔心。射流氣體內(nèi)部反旋轉(zhuǎn)渦對(duì)的發(fā)展過(guò)程可分為3個(gè)特征階段:流動(dòng)位于孔口附近時(shí),反旋轉(zhuǎn)渦對(duì)從壁面逐漸形成,渦核間距和高度不斷增大,影響面積不斷擴(kuò)張;隨著流動(dòng)向下游發(fā)展,反旋轉(zhuǎn)渦對(duì)影響面積不斷收縮直至消失;當(dāng)流動(dòng)發(fā)展至下游某一位置時(shí),反旋轉(zhuǎn)渦對(duì)在射流氣體頂端再次形成,隨著反旋轉(zhuǎn)渦對(duì)的不斷發(fā)展,在平板壁面誘導(dǎo)出2次渦對(duì)。

    流體力學(xué); 氣體與液體兩相流; 橫射流; 數(shù)值模擬; 旋渦結(jié)構(gòu)

    0 引言

    主動(dòng)通氣作為一種有效且易于實(shí)現(xiàn)的流場(chǎng)不穩(wěn)定性調(diào)控措施,廣泛應(yīng)用于各類(lèi)流體機(jī)械以及高速水下和水面航行體中。由于通氣氣體與液體(簡(jiǎn)稱(chēng)氣液)兩相流是一種復(fù)雜的多相湍流[1],往往伴隨著氣、液之間的界面產(chǎn)生、運(yùn)動(dòng)等復(fù)雜物理過(guò)程,將產(chǎn)生復(fù)雜的渦旋結(jié)構(gòu)[2]。

    針對(duì)通氣空化產(chǎn)生的復(fù)雜渦旋結(jié)構(gòu),人們開(kāi)展了諸多實(shí)驗(yàn)與數(shù)值計(jì)算研究。Semenenko[3]根據(jù)實(shí)驗(yàn)結(jié)果和理論研究指出渦環(huán)泄氣是由流動(dòng)分離引起的復(fù)雜旋渦結(jié)構(gòu)產(chǎn)生的,空泡尾流區(qū)充滿(mǎn)了泡沫狀的水氣混合物。Kunz等[4]采用數(shù)值模擬對(duì)軸對(duì)稱(chēng)航行體在一定攻角下的穩(wěn)態(tài)和非穩(wěn)態(tài)自然空化及不含相變的通氣空泡進(jìn)行了研究,對(duì)比分析了航行體表面張力、空泡形態(tài)和阻力,描述了宏觀流場(chǎng)的旋渦結(jié)構(gòu)。文獻(xiàn)[5-7]對(duì)繞回轉(zhuǎn)體通氣空化進(jìn)行了一系列實(shí)驗(yàn)研究,通過(guò)Time-Resolved 粒子圖象測(cè)速(PIV) 技術(shù)對(duì)尾流區(qū)域的旋渦結(jié)構(gòu)進(jìn)行了分析。M?kiharju等[8]通過(guò)X射線(xiàn)實(shí)驗(yàn)和數(shù)值模擬對(duì)平板通氣局部空化進(jìn)行了細(xì)致研究,分析了高雷諾數(shù)下通氣空化特征,研究了通氣空化中的傅汝德數(shù)、雷諾數(shù)和韋伯?dāng)?shù)的影響,通過(guò)對(duì)空泡尾部閉合區(qū)旋渦結(jié)構(gòu)的觀測(cè),發(fā)現(xiàn)閉合區(qū)域上游的湍流邊界層導(dǎo)致壓力波動(dòng)作用于空氣- 水的交界面并使其發(fā)生分離。Ji等[9-10]提出了“三分量”空化模型來(lái)捕捉空泡的發(fā)展,研究發(fā)現(xiàn)不可凝結(jié)氣體抑制了自然空化下回射流的推進(jìn)引起的大尺度旋渦結(jié)構(gòu),并且隨著通氣速率的增加,自然空泡會(huì)較為顯著地被抑制。Wang等[11]對(duì)繞回轉(zhuǎn)體通氣空化空泡脫落機(jī)制進(jìn)行了實(shí)驗(yàn)與數(shù)值計(jì)算研究,發(fā)現(xiàn)通氣和回射流會(huì)在空泡內(nèi)部上游和下游分別形成一個(gè)主渦,同時(shí)在二者之間誘導(dǎo)出一個(gè)2次渦,認(rèn)為2次渦的運(yùn)動(dòng)和能量是空泡發(fā)生斷裂脫落的關(guān)鍵因素。于嫻嫻等[12]對(duì)回轉(zhuǎn)體通氣云狀空化進(jìn)行了研究,發(fā)現(xiàn)通氣和回射流相互作用下形成的旋渦結(jié)構(gòu)會(huì)造成通氣云狀空化發(fā)展過(guò)程中空泡的部分脫落。時(shí)素果等[13]對(duì)繞空化器通氣空化流動(dòng)進(jìn)行了數(shù)值計(jì)算,指出濾波器模型(FBM)能更加準(zhǔn)確地捕捉通氣空泡的渦旋結(jié)構(gòu)。段磊等[14]對(duì)繞回轉(zhuǎn)體渦環(huán)泄氣方式下通氣空化非定常流動(dòng)特性進(jìn)行了研究,結(jié)果表明空泡閉合位置的高壓與空泡區(qū)域的低壓形成較大逆壓梯度,使空泡區(qū)域出現(xiàn)流動(dòng)分離,進(jìn)而在空泡區(qū)域產(chǎn)生復(fù)雜的旋渦結(jié)構(gòu),此旋渦結(jié)構(gòu)與主流相互作用、引起了空泡斷裂。胡曉等[15]比較了大渦模擬(LES)和RNGk-ε湍流模型對(duì)通氣空泡尾部氣體泄漏方式和空泡外形的影響,發(fā)現(xiàn)LES的瞬態(tài)計(jì)算結(jié)果更符合通氣空泡的特性。王復(fù)峰等[16]采用實(shí)驗(yàn)和數(shù)值模擬相結(jié)合的方法對(duì)繞帶不同尺度空化器的通氣空化流場(chǎng)進(jìn)行了研究,發(fā)現(xiàn)大空化器后部流動(dòng)分離明顯,存在更加復(fù)雜的旋渦運(yùn)動(dòng)。劉濤濤等[17]應(yīng)用一種分域湍流模型對(duì)繞回轉(zhuǎn)體通氣超空化流動(dòng)進(jìn)行了研究,指出基于密度修正的模型可以較好地體現(xiàn)前端空泡的可壓縮性,F(xiàn)BM模型可以捕捉尾部的多尺度空泡渦團(tuán)結(jié)構(gòu)。

    當(dāng)引入的主動(dòng)射流氣體以一定的角度從孔狀或有限縫槽結(jié)構(gòu)等較小過(guò)流斷面通入水中時(shí),與外部水流體相互作用,呈現(xiàn)出橫射流特性(JICF)。關(guān)于橫射流研究目前已經(jīng)取得了很多有價(jià)值的成果,主要集中在單相流體介質(zhì),對(duì)多相流則涉及較少。Margason[18]將橫射流流場(chǎng)總結(jié)為3個(gè)主要特征,使其成為一個(gè)理論體系:1)流體從射流孔射出后,在橫流的推動(dòng)下射流向下偏轉(zhuǎn),同時(shí)橫流從射流兩側(cè)繞過(guò)射流,在橫流的剪切作用下,射流形成了一對(duì)反旋轉(zhuǎn)渦對(duì)(CVP),該旋轉(zhuǎn)渦對(duì)主宰著射流孔附近的流動(dòng);2)由于受到射流的阻滯作用,橫流在射流孔前端會(huì)形成分離點(diǎn)和馬蹄渦,馬蹄渦尺度遠(yuǎn)小于CVP;3)繞過(guò)射流后,橫向主流會(huì)在射流孔下游形成非定常的尾跡,尾跡渦的強(qiáng)度是3種渦結(jié)構(gòu)中最弱的。Fric等[19-20]在3渦結(jié)構(gòu)模型的基礎(chǔ)上參照自由射流的特征提出了4渦模型,將射流孔內(nèi)部溢出的剪切層渦環(huán)作為第4個(gè)渦系。Andreopouos等[21]通過(guò)實(shí)驗(yàn)首次發(fā)現(xiàn)在下游某處CVP下方還存在一對(duì)轉(zhuǎn)向與之相反的渦對(duì),稱(chēng)之為2次渦對(duì)。對(duì)于2次渦對(duì)的起源目前還存在爭(zhēng)議。Morton等[22]認(rèn)為2次渦對(duì)應(yīng)該為馬蹄渦的分支,Yuan等[23]通過(guò)LES反駁了這一說(shuō)法,他們將該渦命名為壁面尾跡渦結(jié)構(gòu)。Hale等[24]通過(guò)流油實(shí)驗(yàn)和數(shù)值模擬研究了多孔射流附近的流動(dòng),認(rèn)為2次渦對(duì)是由進(jìn)入射流下方低壓區(qū)的橫流下洗造成的。Yao等[25]對(duì)3孔布置的橫流中射流進(jìn)行了直接數(shù)值模擬,其研究表明,當(dāng)孔間距為1倍孔徑時(shí),中間射流孔的CVP尺度受到兩側(cè)射流的影響,CVP被明顯抑制。Roger等[26]對(duì)并排孔進(jìn)行了LES數(shù)值研究,發(fā)現(xiàn)相鄰兩孔間距是影響射流下游流場(chǎng)渦結(jié)構(gòu)的重要參數(shù),當(dāng)間距較小時(shí),兩股射流的摻混得到增強(qiáng),射流孔下游的CVP尺度增大。吳海玲等[27]對(duì)二維橫向射流進(jìn)行數(shù)值模擬,比較了標(biāo)準(zhǔn)k-ε模型、RNGk-ε模型、Realizablek-ε模型等不同湍流模型對(duì)射流流動(dòng)與傳熱特性預(yù)報(bào)的準(zhǔn)確性,發(fā)現(xiàn)RNGk-ε模型、Realizablek-ε模型對(duì)流場(chǎng)及壁面對(duì)流換熱特性的預(yù)測(cè)優(yōu)于標(biāo)準(zhǔn)k-ε模型。趙馬杰等[28]采用LES方法研究了高雷諾數(shù)下的橫側(cè)射流,指出JICF進(jìn)場(chǎng)迎風(fēng)渦是由于射流出口上游剪切層Kelvin-Helmholtz不穩(wěn)定性引起的。

    為了進(jìn)一步研究通氣空化產(chǎn)生的復(fù)雜旋渦特性,本文利用ANSYS CFX商業(yè)軟件,采用均質(zhì)多相流模型和RNGk-ε湍流模型,并采用Level Set界面捕捉方法對(duì)繞多孔平板流場(chǎng)進(jìn)行數(shù)值模擬,從橫射流角度分析了流場(chǎng)的渦旋結(jié)構(gòu)。

    1 數(shù)值方法

    1.1 基本控制方程

    采用均質(zhì)平衡流模型,則Farve平均的Navier-Stokes方程為

    (1)

    (2)

    式中:ρ為流體密度;ui、uj為速度分量;p為壓強(qiáng);μ和μt分別為層流和湍流黏性系數(shù)。

    1.2 湍流模型

    標(biāo)準(zhǔn)k-ε模型是典型的雷諾時(shí)均化(RANS)湍流模型,它把渦黏系數(shù)和湍動(dòng)能及湍動(dòng)能耗散聯(lián)系在一起,該模型由Launder等于1972年提出[29],對(duì)于均相平衡流動(dòng)的數(shù)值計(jì)算,標(biāo)準(zhǔn)k-ε模型的控制方程為

    Pt-ρmε,

    (3)

    (4)

    式中:Pt為湍動(dòng)能生成項(xiàng);ρm、μm分別為混合密度和混合湍流黏性;μt定義為

    (5)

    Cε1、Cε2、σε、σk、Cμ分別為模型常數(shù)。

    RNGk-ε模型由Yakhot等[30]采用“重整化群”的數(shù)學(xué)方法推導(dǎo)得出,與標(biāo)準(zhǔn)k-ε模型的不同之處主要是在ε方程中增加了R項(xiàng),即

    (6)

    式中:R為流場(chǎng)變化度,定義為

    (7)

    經(jīng)化簡(jiǎn)得到RNGk-ε模型的ε方程為

    (8)

    相應(yīng)的模型系數(shù)取值由理論分析得出,具體結(jié)果[31]列于表1中。在η<η0區(qū)域,RNGk-ε模型中的系數(shù)Cε1小于標(biāo)準(zhǔn)k-ε模型中的系數(shù)Cε1;在η>η0區(qū)域,RNGk-ε模型中的系數(shù)Cε1大于標(biāo)準(zhǔn)k-ε模型中的系數(shù)Cε1,因此與標(biāo)準(zhǔn)k-ε模型相比,在高應(yīng)變率及流線(xiàn)彎曲程度較大的流場(chǎng)中,RNGk-ε模型將產(chǎn)生較小的湍流黏性系數(shù)。

    表1 各湍流模型ε方程中Sε的表達(dá)式及模型常數(shù)

    1.3LevelSet方法

    在固定的歐拉坐標(biāo)系中,含有相界面兩相介質(zhì)的流動(dòng)可用Navier-Stokes方程[32-33]描述為

    (9)

    (10)

    ρm(x)=ρa(bǔ)+(ρl-ρa(bǔ))Hε(φ(x)),

    (11)

    μm(x)=μa+(μl-μa)Hε(φ(x)),

    (12)

    式中:ρa(bǔ)、μa分別為氣相密度和氣相湍流黏性;ρl、μl分別為液相密度和液相湍流黏性.

    Level Set方法的主要思想是將界面定義為一個(gè)函數(shù)的零等值面(線(xiàn)),即φ(x,t)=0. 令φ以適當(dāng)?shù)乃俣纫苿?dòng),使其零等值的面就是物質(zhì)界面。

    Level Set函數(shù)φ定義為

    (13)

    Heaviside函數(shù)Hε被定義為

    Hε(d)=

    (14)

    式中:ε1是一個(gè)小量規(guī)整參數(shù),恒為正。Heaviside函數(shù)也可以作為區(qū)分計(jì)算區(qū)域內(nèi)介質(zhì)種類(lèi)的方法和指標(biāo)。

    借助于Level set函數(shù)φ,相界面的內(nèi)在集合特性參數(shù)可被確定為

    法向向量

    (15)

    界面曲率

    (16)

    因而(9)式中的表面張力項(xiàng)可表示為

    (17)

    這樣,(9)式就可以像求解單相流體的Navier-Stokes方程一樣方便地求解了。

    1.4 計(jì)算邊界條件及設(shè)置

    圖1給出了平板幾何參數(shù),平板全長(zhǎng)300 mm,寬70 mm,其表面設(shè)置4個(gè)通氣孔,直徑d=2.6 mm,間距(孔心之間距離)L=6.6d,位于距工作段前端60 mm處,工作段前后各有一段光滑過(guò)渡的導(dǎo)流段。圖2為計(jì)算區(qū)域及邊界設(shè)置條件。計(jì)算域總長(zhǎng)3 500 mm,高度190 mm,邊界條件設(shè)置為:速度入口,壓力出口,速度、壓力和含氣率在流動(dòng)方向的導(dǎo)數(shù)為常數(shù),通氣孔為質(zhì)量流量入口,固壁為絕熱、無(wú)滑移壁面條件。實(shí)驗(yàn)和數(shù)值計(jì)算中采用的來(lái)流速度v=5.1 m/s,通氣體積流量系數(shù)Qs=m/(4ρa(bǔ)vd2)=1.12,其中m為總通氣質(zhì)量流量(kg/s),ρa(bǔ)為實(shí)驗(yàn)工況下氣體密度(單位:kg/m3,實(shí)驗(yàn)工況pe=4 atm,Te=25 ℃). 數(shù)值計(jì)算中時(shí)間步長(zhǎng)設(shè)定為0.000 5 s,總計(jì)算時(shí)間為0.1 s.

    圖1 計(jì)算模型幾何參數(shù)Fig.1 Geometrical parameters of plate model

    圖2 計(jì)算區(qū)域及邊界設(shè)置條件Fig.2 Computational domain and boundary conditions

    圖3給出了平板工作段及孔口附近的網(wǎng)格圖。為了更準(zhǔn)確地捕捉平板表面及孔口附近的旋渦結(jié)構(gòu),在近壁處及通氣孔附近進(jìn)行網(wǎng)格加密。近壁面y+值在30~150之間,滿(mǎn)足壁面函數(shù)要求。

    圖3 平板工作段及孔口附近的網(wǎng)格圖Fig.3 Computational grids around plate

    分別取網(wǎng)格數(shù)為50萬(wàn)、100萬(wàn)、150萬(wàn)、180萬(wàn)和200萬(wàn)的網(wǎng)格進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證,其中為保證網(wǎng)格精度使計(jì)算準(zhǔn)確,網(wǎng)格選擇的標(biāo)準(zhǔn)為所有網(wǎng)格質(zhì)量均在0.7以上,所得到的阻力系數(shù)分布如圖4所示。

    圖4 網(wǎng)格無(wú)關(guān)性驗(yàn)證Fig.4 Verification of grid independence

    由圖4可知,當(dāng)網(wǎng)格數(shù)達(dá)到150~200萬(wàn)時(shí),可以發(fā)現(xiàn)所得到的阻力系數(shù)變化不大,因此在節(jié)省計(jì)算時(shí)間和保證計(jì)算結(jié)果精度的前提下,將網(wǎng)格數(shù)取為180萬(wàn)。

    2 結(jié)果與討論

    2.1 數(shù)值與實(shí)驗(yàn)比較

    圖5給出了特定時(shí)刻數(shù)值模擬計(jì)算得到的空泡形態(tài)與實(shí)驗(yàn)結(jié)果的對(duì)比。從圖5中可以看出,本文采用的數(shù)值模擬方法計(jì)算得到的空泡形態(tài)發(fā)展過(guò)程與實(shí)驗(yàn)結(jié)果具有較好的一致性,即空泡緊貼通氣孔后部形成連續(xù)條狀,隨著流動(dòng)向下游發(fā)展,空泡在通氣孔后聚集并不斷向周向膨脹、發(fā)生交匯,水氣交界面清晰。由于實(shí)驗(yàn)中無(wú)法完全確保氣體持續(xù)均勻通入,實(shí)驗(yàn)觀測(cè)到的氣泡非定常特性與數(shù)值計(jì)算存在較大差異。為了進(jìn)一步說(shuō)明數(shù)值計(jì)算方法的可行性,圖6給出了實(shí)驗(yàn)觀測(cè)的氣泡外形與數(shù)值計(jì)算的對(duì)比,其中實(shí)驗(yàn)觀測(cè)結(jié)果為多次結(jié)果的平均值,可以看出,數(shù)值計(jì)算得到的氣泡外形與實(shí)驗(yàn)結(jié)果吻合較好。本文中采用的模型孔間距較大,各孔間流動(dòng)發(fā)展規(guī)律基本一致,因此在接下來(lái)的分析中主要選取中心處單孔的流動(dòng)特性進(jìn)行分析。

    圖5 實(shí)驗(yàn)與數(shù)值模擬得到的氣泡形態(tài)對(duì)比Fig.5 Comparison of simulated and experimental cavity shapes

    圖6 實(shí)驗(yàn)與數(shù)值模擬得到的氣泡外形對(duì)比Fig.6 Comparison of simulated and experimental cavity geometric shapes

    2.2 孔口附近旋渦結(jié)構(gòu)

    圖7和圖8分別給出了孔口附近Oxy平面內(nèi)渦量分布和相應(yīng)的流線(xiàn)分布(以z/d=-3.3為圓心)。從圖7和圖8中可以看出,射流氣體在液相橫流作用下逐漸轉(zhuǎn)變方向,由與橫流正交變?yōu)榕c橫流流動(dòng)方向平行。由于受到射流的阻滯作用,橫流在孔口上游x/d=0.583處形成分離鞍點(diǎn)和馬蹄渦,此分離鞍點(diǎn)隨距壁面高度的增加逐漸靠近孔心,形成分離線(xiàn)。橫流繞過(guò)射流后,在孔口下游形成尾跡渦,此尾跡渦主要來(lái)源于壁面邊界層。

    圖7 Oxy平面渦量云圖(z/d=-3.3)Fig.7 Contour of spanwise vorticity on vertical plane Oxy (z/d=-3.3)

    圖8 Oxy平面流線(xiàn)分布(z/d=-3.3)Fig.8 Streamline patterns on vertical plane Oxy(z/d=-3.3)

    為進(jìn)一步分析孔口附近的旋渦結(jié)構(gòu),圖9給出了距壁面不同高度Oxz平面內(nèi)的流線(xiàn)分布(y/d=0,y/d=0.07,y/d=0.15,y/d=0.17),其中藍(lán)色線(xiàn)條為圓孔輪廓。由圖9可以看出,當(dāng)y/d<0.17時(shí),流線(xiàn)均存在3個(gè)分離點(diǎn)、2個(gè)螺旋節(jié)點(diǎn)、1個(gè)臨界節(jié)點(diǎn)和3條分離線(xiàn)。其中分離點(diǎn)S1位于孔口上游,與馬蹄渦相關(guān)聯(lián),由其伸出的兩條分離線(xiàn)L1和L2指出了馬蹄渦的跡線(xiàn),該分離線(xiàn)繞過(guò)氣流后向外擴(kuò)張發(fā)展,逐漸遠(yuǎn)離孔中心線(xiàn)。通氣孔下游的分離點(diǎn)S3則對(duì)應(yīng)于定常流動(dòng)條件下剛性圓柱體繞流下游的滯止點(diǎn),位于孔口下游,與該分離點(diǎn)直接連接的兩個(gè)螺旋節(jié)點(diǎn)N2和N3,在孔口下游邊沿形成兩個(gè)分離旋渦,該分離旋渦主要由射流氣體背面逆壓梯度造成,這一現(xiàn)象與文獻(xiàn)[34]中基本一致。分離點(diǎn)S2位于分離旋渦下游,橫流繞過(guò)射流后于分離點(diǎn)S2處發(fā)生匯合并沿著流動(dòng)的對(duì)稱(chēng)線(xiàn)向下游發(fā)展。圖10給出了分離點(diǎn)S1、S2、S3以及分離旋渦渦核位置L隨平板高度的變化趨勢(shì)。結(jié)合圖9和圖10可以看出,隨著y/d的增大,分離點(diǎn)S1逐漸靠近孔心,S2逐漸遠(yuǎn)離孔心,當(dāng)y/d增大到0.17時(shí),分離點(diǎn)S1消失,S2位置保持相對(duì)穩(wěn)定。當(dāng)y/d<0.06時(shí),分離點(diǎn)S3與分離旋渦渦核位置基本重合,隨著y/d的增大,分離點(diǎn)S3與分離旋渦渦核逐漸遠(yuǎn)離孔心,且兩者間的距離逐漸增大,同時(shí)分離旋渦逐漸凸顯,影響范圍不斷擴(kuò)大。

    圖9 Oxz平面流線(xiàn)分布Fig.9 Streamline patterns on a vertical plane Oxz

    圖10 滯止點(diǎn)和分離旋渦渦核偏移孔心距離與y/d的關(guān)系Fig.10 Relation among stagnation points, counter-rotating vortex cores and y/d

    2.3 射流氣體內(nèi)部旋渦結(jié)構(gòu)

    為了進(jìn)一步對(duì)射流氣體內(nèi)部的旋渦結(jié)構(gòu)進(jìn)行分析,圖11給出了不同x/d位置處Oyz平面內(nèi)流線(xiàn)分布和氣相體積分?jǐn)?shù)云圖,可以看出在射流內(nèi)部形成了以射流孔中心線(xiàn)為對(duì)稱(chēng)軸的一對(duì)渦量相等、旋轉(zhuǎn)方向相反的渦對(duì),即CVP. 圖12給出了CVP渦核位置隨流向距離的變化趨勢(shì),其中,橫坐標(biāo)表示沿流向距孔心的距離,縱坐標(biāo)分別表示CVP渦核偏離孔中心線(xiàn)的距離E和距平板壁面高度H.

    圖11 不同流向截面位置流線(xiàn)分布Fig.11 Streamline patterns on different vertical planes Oyz

    圖12 CVP渦核位置隨流向距離變化曲線(xiàn)Fig.12 Relation between CVP cores and x/d

    結(jié)合圖11和圖12可以看出,CVP隨流向的發(fā)展過(guò)程可以劃分為以下3個(gè)特征階段:

    1)初步發(fā)展階段(x/d=0~10):當(dāng)x/d=0時(shí),CVP尚未形成,孔心處流線(xiàn)平行于y軸,隨著向孔邊界的移動(dòng),流線(xiàn)逐漸偏離y軸,趨于平行x軸,氣體流線(xiàn)整體呈現(xiàn)發(fā)散狀。隨著流動(dòng)向下游發(fā)展,在x/d=0.5處CVP開(kāi)始形成,此位置與分離旋渦的分離點(diǎn)S3位置基本一致,說(shuō)明CVP的形成可能與橫流繞流分離旋渦有關(guān)。由于壁面的存在,此時(shí)尚未形成較為完整的渦對(duì)。當(dāng)x/d=1時(shí),射流氣體內(nèi)部已經(jīng)形成較為完整的CVP,此階段渦核偏離距離和高度呈現(xiàn)快速增長(zhǎng)趨勢(shì)。從x/d=1到x/d=5,由于黏性的作用,CVP流向斷面的影響面積不斷擴(kuò)張,渦核偏離距離和高度呈現(xiàn)緩慢增長(zhǎng)趨勢(shì)。隨著流動(dòng)進(jìn)一步向下游發(fā)展,CVP影響面積變化不顯著,渦核偏離距離持續(xù)增大,渦核高度發(fā)展趨勢(shì)逐漸趨于平緩,在x/d=10處達(dá)到峰值。

    2)衰弱階段(x/d=10~35):這一階段,CVP流向斷面的影響面積不斷收縮,渦核偏離距離持續(xù)增大,在x/d=25處達(dá)到最大,對(duì)應(yīng)的渦核高度逐漸降低。當(dāng)流動(dòng)發(fā)展到x/d=30時(shí),CVP基本消失,相應(yīng)的渦核偏離距離和高度趨于0.

    3)2次發(fā)展階段(x/d=35~85):這一階段,CVP再次產(chǎn)生并不斷發(fā)展,其影響面積遠(yuǎn)遠(yuǎn)大于初次CVP的影響面積。當(dāng)x/d=36時(shí),射流氣體頂端流線(xiàn)發(fā)生扭曲,但仍未產(chǎn)生CVP. 隨著流動(dòng)向下游發(fā)展至x/d=37時(shí),在射流氣體內(nèi)部再次形成較為完整的CVP,且影響面積較小。結(jié)合圖11可以看出,剛形成的CVP位于射流氣體頂端孔中心線(xiàn)附近,渦核間距較小,渦核高度較高,且左渦核高度要高于右渦核。從x/d=37到x/d=75,CVP影響面積不斷擴(kuò)張,渦核間距在x/d=45前呈現(xiàn)快速增長(zhǎng),躍過(guò)x/d=45后,渦核間距持續(xù)增大但增長(zhǎng)速度減緩。渦核高度變化趨勢(shì)與渦核間距變化趨勢(shì)存在顯著差異,在此階段渦核高度持續(xù)下降至x/d=60,躍過(guò)x/d=60后,渦核高度保持相對(duì)平穩(wěn)趨勢(shì)。從圖10(c)可以看出,在x/d=60處,CVP在平板壁面誘導(dǎo)出2次渦對(duì),這一現(xiàn)象與Andreopoulos等[35]觀測(cè)的實(shí)驗(yàn)現(xiàn)象基本一致。隨著流動(dòng)向下游發(fā)展,2次渦對(duì)被卷吸入CVP后消失, CVP影響范圍擴(kuò)大。在x/d=85位置處,2次渦對(duì)再次產(chǎn)生。由于2次渦對(duì)的周期性產(chǎn)生- 消失,使渦核相對(duì)保持在某一高度。

    3 結(jié)論

    本文基于均相流模型,采用RNGk-ε湍流模型并耦合對(duì)Level Set界面捕捉方法對(duì)繞多孔平板流場(chǎng)進(jìn)行數(shù)值模擬,分析了通氣氣液兩相橫射流旋渦特性,主要結(jié)論如下:

    1)液相橫流受到射流氣體的阻礙作用,在孔口上游形成分離鞍點(diǎn)和馬蹄渦,此分離鞍點(diǎn)隨距壁面高度的增加逐漸靠近孔心,形成分離線(xiàn)。

    2)由于射流氣體與液相橫流相互作用,橫流在射流背面出現(xiàn)繞流,形成兩個(gè)較為封閉的分離旋渦,此分離旋渦隨距壁面高度的增加逐漸遠(yuǎn)離孔心。

    3)射流氣體內(nèi)部CVP的發(fā)展過(guò)程可分為3個(gè)特征階段:初步發(fā)展階段、衰弱階段和2次發(fā)展階段。流動(dòng)位于孔口附近時(shí),CVP從壁面逐漸形成,渦核間距和高度不斷增大,影響面積不斷擴(kuò)張;隨著流動(dòng)向下游發(fā)展,CVP影響面積不斷收縮直至消失;當(dāng)流動(dòng)發(fā)展至下游某一位置時(shí),CVP在射流氣體頂端再次形成,隨著CVP的不斷發(fā)展,在平板壁面誘導(dǎo)出2次渦對(duì)。

    References)

    [1] Brennen C E. Fundamentals of multiphase flow[M]. Cambridge, UK:Cambridge University Press, 2005.

    [2] Cantro F, Crespo A, Manuel F, et al. Equlibrium of ventilated cavities in tip vortices[J]. Journal of Fluids Engineering, 1997, 119(4):759-767.

    [3] Semenenko V N. Artificial supercavitation physics and calculation, ADP012080[R]. Kiev, Ukraine: Institute of Hydeomechanics of Ukrainian Academy of Science, 2001.

    [4] Kunz R F, Boger D A, Chyczewski T S, et al. Multiphase CFD analysis of natural and ventilated cavitation about submerged bodies[C]∥3th ASME/JSME Joint Fluids Engineering Conference. San Francisco, CA, US: ASME, 1999.

    [5] Wosnik M, Schauer T J, Arndt R E A. Experimental study of ventilated supercavitating vehicle[C]∥5th International Symposium on Cavitation. Osaka, Japan: Osaka University, 2003.

    [6] Schauer T. An experimental study of a ventilated supercavity vehicle[D]. Minneapolis, MN, US: University of Minnesota, 2003.

    [7] Wosnik M, Schauer T, Arndt R E A. Measurements in high void-fraction bubbly wakes created by ventilated supercavitation[J]. Journal of Fluids Engineering, 2013,135(1):531-538.

    [8] M?kiharju S A. The dynamics of ventilated partial cavities over a wide range of Reynolds numbers and quantitative 2D X-ray densitometry for multiphase flow[D]. MI, US: University of Michigan, 2012.

    [9] Ji B, Luo X W, Zhang Y, et al. A three-component model suitable for natural and ventilated cavitation[J]. Chinese Physics Letter, 2010, 27(9): 38-43.

    [10] Ji B, Luo X W, Peng X X, et al. Numerical investigation of the ventilated cavitating flow around an underwater vehicle based on a three-component cavitation model[J]. Journal of Hydrodynamics, 2010, 22(6): 753-759.

    [11] Wang Y W, Huang C G, Wei Y P, et al. Ventilated partial cavitation shedding caused by the interaction of re-entry and air injection[C]∥Proceedings of the 8th International Symposium on Cavitation. Singapore: Nanyang Technological University, 2012.

    [12] 于嫻嫻, 王一偉, 黃晨光, 等. 軸對(duì)稱(chēng)航行體通氣云狀空化非定常特征研究[J]. 船舶力學(xué), 2014, 18(5): 499-506. YU Xian-xian, WANG Yi-wei, HUANG Chen-guang, et al. Unsteady characteristics of ventilated cloud cavity around symmetrical bodies[J]. Journal of Ship Mechanics, 2014, 18(5):499-506. (in Chinese)

    [13] 時(shí)素果, 王國(guó)玉, 余志毅, 等. FBM湍流模型在非定常通氣超空化流動(dòng)計(jì)算中的評(píng)價(jià)與應(yīng)用[J]. 船舶力學(xué), 2012, 16(10):1100-1106. SHI Su-guo, WANG Guo-yu, YU Zhi-yi, et al. Evaluation of the filter based turbulence model for computation of unsteady ventilated supercavitating flows[J]. Journal of Ship Mechanics, 2012, 16(10):1100-1106. (in Chinese)

    [14] 段磊, 王國(guó)玉, 付細(xì)能. 渦環(huán)泄氣方式下通氣空化的非定常流動(dòng)特性研究[J]. 兵工學(xué)報(bào), 2014, 35(5):712-718. DUAN Lei, WANG Guo-yu, FU Xi-neng. Research of the unsteady characteristics of ventilated cavitating flows in the form of gas leakage by toroidal vortex[J]. Acta Armamentarii, 2014, 35(5): 712-718. (in Chinese)

    [15] 胡曉, 郜冶, 彭輝. 湍流模型對(duì)空泡形態(tài)影響的數(shù)值研究[J]. 計(jì)算力學(xué)學(xué)報(bào), 2015,32(1): 129-135. HU Xiao, GAO Ye, PENG Hui. Numerical investigation on influence of turbulence model on cavity shape[J]. Chinese Journal of Computational Mechanics, 2015,32(1): 129-135. (in Chinese)

    [16] 王復(fù)峰, 王國(guó)玉, 張敏弟, 等. 帶空化器回轉(zhuǎn)體空化流場(chǎng)研究[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2015, 36(7):959-964. WANG Fu-feng, WANG Guo-yu, ZHANG Min-di, et al. Study of cavitating flows around an axisymmetric body with cavitators[J]. Journal of Harbin Engineering University, 2015, 36(7):959-964. (in Chinese)

    [17] 劉濤濤, 王國(guó)玉, 段磊. 基于實(shí)驗(yàn)觀測(cè)的分域湍流模型在通氣超空化中的評(píng)價(jià)[J]. 北京理工大學(xué)學(xué)報(bào), 2016, 36(3): 247-251. LIU Tao-tao, WANG Guo-yu, DUAN Lei. Assessment of a modified turbulence model based experiment results for ventilated supercavity[J]. Transactions of Beijing Institute of Technology, 2016, 36(3): 247-251. (in Chinese)

    [18] Margason R J. Fifty years of jet in cross flow research[C]∥Proceedings of the AGARD Symposium on Computational & Experimental Assessment of Jets in Cross Flow. London, UK: Advisory Group for Aerospace Research and Development, 1993.

    [19] Fric T F. Structure in the near field of the transverse jet[D]. Pasadena, CA, US: California Institute of Technology, 1990.

    [20] Fric T F, Roshko A. Vortical structure in the wake of a transverse jet[J]. Journal of Fluid Mechanics, 1994, 279:1-47.

    [21] Andreopoulos J, Rodi W. Experimental investigation of jets in a cross flow[J]. Journal of Fluid Mechanics, 1984, 138:93-127.

    [22] Morton B R, Ibberson A. Jets deflected across flow[J]. Experimental Thermal and Fluid Science, 1996, 12(2):112-133.

    [23] Yuan L L, Street R L, Ferziger J H. Large eddy simulation of a round jet in cross-flow[J]. Journal of Fluid Mechanics, 1999, 379:71-104.

    [24] Hale C A, Plesniak M W, Ramadhyani S. Structural features and surface heat transfer associated with a row of short hole jets in cross flow[J]. International Journal of Heat and Fluid Flow, 2000, 21:542-553.

    [25] Yao Y, Maidi M. Direct numerical simulation of single and multiple square jets in cross flow[J]. Journal of Fluids Engineering, 2011, 133(3):1201.

    [26] Roger F, Gourara A, Most J M, et al. Numerical investigation on the reactive gas mixing through interaction between twin square jets side-by-side and a cross flow[J]. Journal of Chemical Engineering, 2004, 238:45-55.

    [27] 吳海玲, 陳聽(tīng)寬, 羅毓珊. 應(yīng)用不同紊流模型的二維橫向射流傳熱數(shù)值模擬研究[J]. 西安交通大學(xué)學(xué)報(bào), 2001, 35(9):903-907. WU Hai-ling, CHEN Ting-kuan, LUO Yu-shan. Numerical simulation of 2D jet to cross flow heat transfer with different turbulence models[J]. Journal of Xi’an Jiaotong University, 2001, 35(9):903-907. (in Chinese)

    [28] 趙馬杰, 曹長(zhǎng)敏, 張宏達(dá), 等. 高雷諾數(shù)湍流橫側(cè)射流的大渦模擬[J]. 推進(jìn)技術(shù), 2016, 37(5): 834-843. ZHAO Ma-jie, CAO Chang-min, ZHANG Hong-da, et al. Large eddy simulation of a high Reynolds number turbulent jet in cross flow[J]. Journal of Propulsion Technology,2016, 37(5): 834-843. (in Chinese)

    [29] Launder B E, Spalding D B. The numerical computation of turbulent flows[J]. Computational Methods in Applied Mechanics and Engineering, 1974, 3(2): 269-289.

    [30] Yakhot V, Orszag S A. Renormalization group analysis of turbulence[J]. Journal of Scientific Computing, 1986, 1(1): 3-5.

    [31] Speziale C G, Thangam S. Analysis of RNG based turbulence model for separated flows[J]. Journal of Engineering Science, 1992, 30(10): 1379-1388.

    [32] Huang B, Wang G Y, Zhang M D, et al. Level set model for simulation of cavitation flows[J]. Journal of Ship Mechanics, 2011, 15(3):207-216.

    [33] 吳欽, 王國(guó)玉, 付細(xì)能, 等. 繞平板氣液兩相流數(shù)值計(jì)算方法[J]. 北京理工大學(xué)學(xué)報(bào), 2014, 34(5):475-500. WU Qin, WANG Guo-yu, FU Xi-neng, et al. Numerical study on multiphase flows around a plane[J]. Transactions of Beijing Institute of Technology, 2014, 34(5):475-500. (in Chinese)

    [34] Kelso R M, Lim T T, Perry A E. An experimental study of rounds jets in cross flow[J]. Journal of Fluid Mechanics, 1996, 306:111-144.

    [35] Andreopoulos J, Rodi W. Experimental investigation of jets in a cross flow[J]. Experimental Thermal and Fluid Science, 1996, 12(2):112-133.

    Analysis of Vortex Dynamics of Gas-liquid Two-phase Crossflows around a Porous Plate

    LIU Tao-tao1, WANG Guo-yu1, ZHANG Nai-min2, HUANG Biao1

    (1.School of Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China; 2.Beijing Institute of Astronautical System Engineering, Beijing 100076, China)

    In order to investigate the vortex structures of multiphase flows around a porous plate, the flow characteristics of gas-liquid two-phase crossflows were numerically simulated by using the RNGk-εturbulence modeland the Level Set model. The simulated results are compared with the experimental results. The research results show that the separation point and the horseshoe vortex are formed since the crossflow is blocked by gas jet, which can be observed at the upstream of the jet hole. With the increases in the distance away from the wall, the separation point gradually gets close to the jet exit hole. The crossflow detours the jet and forms two counter-rotating vortices on the lateral edges of jet, and the vortex evolution depends upon the distance away from the wall. The counter-rotating vortex pairs (CVP) are formed in the jet region. The development process of the vortex pairs can be devided into 3 stages: in the near wake-region (i.e., close to the edge of the jet exit hole), the counter-rotating vortex pairs gradually form on the wall, and the increases in the heights of the vortex cores ae well as the distance between the cores lead to the expansion of vortex area. As the flow develops toward downstream, the area of CVP shrinks and even disappears entirely. Further downstream, the CVP appears again on the top of the jet corresponding to the formation of secondary vortex pairs in the near-wall region.

    fluid mechanics; gas-liquid two-phase flow; crossflow jet; numerical simulation; vortex structure

    2016-10-24

    國(guó)家自然科學(xué)基金項(xiàng)目(51239005)

    劉濤濤(1989—),男,博士研究生。E-mail:liutaotao_0708@126.com

    王國(guó)玉(1961—),男,教授,博士生導(dǎo)師。E-mail: wangguoyu@bit.edu.cn

    TV131.3+2

    A

    1000-1093(2017)07-1375-10

    10.3969/j.issn.1000-1093.2017.07.016

    猜你喜歡
    旋渦空泡空化
    功率超聲作用下鋼液中空化泡尺寸的演變特性
    鋼鐵釩鈦(2023年5期)2023-11-17 08:48:34
    水下航行體雙空泡相互作用數(shù)值模擬研究
    小心,旋渦來(lái)啦
    大班科學(xué)活動(dòng):神秘的旋渦
    旋渦笑臉
    山間湖
    三維扭曲水翼空化現(xiàn)象CFD模擬
    不同運(yùn)動(dòng)形式下水物相互作用空化數(shù)值模擬
    基于LPV的超空泡航行體H∞抗飽和控制
    基于CFD的對(duì)轉(zhuǎn)槳無(wú)空泡噪聲的仿真預(yù)報(bào)
    船海工程(2015年4期)2016-01-05 15:53:28
    美女cb高潮喷水在线观看| 色视频www国产| 亚洲人成网站在线播| 成人精品一区二区免费| 精品一区二区三区视频在线观看免费| 亚洲美女黄片视频| 91久久精品国产一区二区三区| 国产精品久久视频播放| 国产亚洲欧美98| 日韩一区二区视频免费看| 黄色日韩在线| 99热网站在线观看| 精品久久久久久久末码| 尤物成人国产欧美一区二区三区| 网址你懂的国产日韩在线| 国产精品久久久久久亚洲av鲁大| 变态另类丝袜制服| 给我免费播放毛片高清在线观看| 老师上课跳d突然被开到最大视频| 国产精品一区二区三区四区免费观看 | 国产精品乱码一区二三区的特点| 少妇人妻一区二区三区视频| а√天堂www在线а√下载| 在线看三级毛片| a级一级毛片免费在线观看| 国产一级毛片七仙女欲春2| 亚洲综合色惰| 美女cb高潮喷水在线观看| 亚洲天堂国产精品一区在线| 国产极品精品免费视频能看的| 精品一区二区免费观看| 俺也久久电影网| 国产91av在线免费观看| 99久久九九国产精品国产免费| 国产一区二区三区av在线 | 精品欧美国产一区二区三| 一个人免费在线观看电影| 亚洲va在线va天堂va国产| 亚洲综合色惰| 亚洲最大成人av| 一个人观看的视频www高清免费观看| 成年免费大片在线观看| av在线亚洲专区| 日本五十路高清| 国产一区二区激情短视频| 1024手机看黄色片| 亚洲av中文字字幕乱码综合| 亚洲av中文av极速乱| 有码 亚洲区| videossex国产| 午夜精品国产一区二区电影 | 欧美日韩综合久久久久久| 热99在线观看视频| 亚洲婷婷狠狠爱综合网| 国产久久久一区二区三区| 国产精品一区二区三区四区免费观看 | 国产淫片久久久久久久久| 国产综合懂色| 两个人的视频大全免费| 久久久久久大精品| 色av中文字幕| 日韩精品中文字幕看吧| 香蕉av资源在线| 欧美xxxx黑人xx丫x性爽| 国产久久久一区二区三区| 欧美最新免费一区二区三区| 国产成年人精品一区二区| 美女cb高潮喷水在线观看| 欧美一级a爱片免费观看看| 成年女人看的毛片在线观看| 国产精品一区二区三区四区久久| 狂野欧美白嫩少妇大欣赏| 日韩一区二区视频免费看| 国语自产精品视频在线第100页| 18禁在线播放成人免费| 麻豆国产97在线/欧美| 欧美成人a在线观看| 男人和女人高潮做爰伦理| 欧美另类亚洲清纯唯美| 看片在线看免费视频| 波野结衣二区三区在线| 中文亚洲av片在线观看爽| 亚洲真实伦在线观看| 久久这里只有精品中国| 国产一区二区在线av高清观看| 日韩欧美免费精品| 亚洲精品影视一区二区三区av| 日本欧美国产在线视频| 久久99热6这里只有精品| 日本在线视频免费播放| 美女被艹到高潮喷水动态| 久久鲁丝午夜福利片| 精品熟女少妇av免费看| 欧美日韩在线观看h| 国产精品免费一区二区三区在线| 男人的好看免费观看在线视频| 乱码一卡2卡4卡精品| 你懂的网址亚洲精品在线观看 | 亚洲第一区二区三区不卡| 国产精品人妻久久久影院| 在现免费观看毛片| 亚洲无线观看免费| 热99在线观看视频| 熟女人妻精品中文字幕| 免费看光身美女| 男人舔女人下体高潮全视频| 亚洲最大成人手机在线| 99久久精品一区二区三区| 亚洲激情五月婷婷啪啪| 岛国在线免费视频观看| 久久久久久久亚洲中文字幕| 22中文网久久字幕| 日日干狠狠操夜夜爽| 直男gayav资源| 国产视频一区二区在线看| 亚州av有码| 在现免费观看毛片| av在线老鸭窝| 人人妻,人人澡人人爽秒播| 免费大片18禁| 中文字幕人妻熟人妻熟丝袜美| 国产精品国产高清国产av| 亚洲自拍偷在线| 国产高清有码在线观看视频| 国产精品一区二区免费欧美| 此物有八面人人有两片| 欧美在线一区亚洲| 联通29元200g的流量卡| 少妇人妻一区二区三区视频| 亚洲精品影视一区二区三区av| 亚洲av成人精品一区久久| 国产高清视频在线播放一区| 日本爱情动作片www.在线观看 | 国产精品一区二区三区四区久久| 别揉我奶头~嗯~啊~动态视频| 亚洲人成网站在线观看播放| 深夜a级毛片| 亚洲高清免费不卡视频| 精品久久久久久久人妻蜜臀av| 无遮挡黄片免费观看| 久久6这里有精品| 欧美国产日韩亚洲一区| 寂寞人妻少妇视频99o| 久久精品综合一区二区三区| 精品99又大又爽又粗少妇毛片| 欧美三级亚洲精品| 亚洲性夜色夜夜综合| 男人狂女人下面高潮的视频| 国产免费一级a男人的天堂| 听说在线观看完整版免费高清| 男人狂女人下面高潮的视频| 亚洲一区二区三区色噜噜| 九九热线精品视视频播放| 又爽又黄无遮挡网站| 国产白丝娇喘喷水9色精品| 国产精品一及| 欧美3d第一页| 99视频精品全部免费 在线| 色哟哟·www| 久久久欧美国产精品| 亚洲激情五月婷婷啪啪| 卡戴珊不雅视频在线播放| 丝袜喷水一区| 成人三级黄色视频| 毛片女人毛片| 人人妻人人看人人澡| 少妇丰满av| 国产亚洲精品久久久com| 18+在线观看网站| 18禁黄网站禁片免费观看直播| 午夜福利在线观看吧| 青春草视频在线免费观看| 波多野结衣高清无吗| 一夜夜www| av黄色大香蕉| 人人妻人人澡欧美一区二区| 春色校园在线视频观看| 欧美人与善性xxx| 九九久久精品国产亚洲av麻豆| 国产精品嫩草影院av在线观看| 国产精品亚洲美女久久久| 草草在线视频免费看| 成人性生交大片免费视频hd| 少妇高潮的动态图| 亚州av有码| 麻豆av噜噜一区二区三区| 成年版毛片免费区| 床上黄色一级片| 久久久国产成人精品二区| av福利片在线观看| 欧美性感艳星| 观看免费一级毛片| 深夜a级毛片| 日本欧美国产在线视频| 亚洲av五月六月丁香网| 变态另类丝袜制服| 又爽又黄a免费视频| 神马国产精品三级电影在线观看| 插阴视频在线观看视频| 国产精品亚洲美女久久久| av免费在线看不卡| 色av中文字幕| 欧美一区二区精品小视频在线| 日韩欧美免费精品| 真人做人爱边吃奶动态| 久久久久久伊人网av| 国产黄色视频一区二区在线观看 | 亚洲中文日韩欧美视频| 免费电影在线观看免费观看| 国产69精品久久久久777片| 亚洲国产高清在线一区二区三| 久久精品综合一区二区三区| 亚洲人成网站高清观看| 久久精品国产亚洲av香蕉五月| 日韩 亚洲 欧美在线| 特级一级黄色大片| 亚洲内射少妇av| 欧美极品一区二区三区四区| 亚洲av中文av极速乱| 久久久久国内视频| 伦精品一区二区三区| 日日撸夜夜添| 成年av动漫网址| 国产精品一区二区性色av| 成人午夜高清在线视频| 国产亚洲精品综合一区在线观看| 久久精品影院6| 亚洲人成网站在线播| 日韩av不卡免费在线播放| 亚洲婷婷狠狠爱综合网| 成人高潮视频无遮挡免费网站| 亚洲内射少妇av| 亚洲不卡免费看| 天天躁夜夜躁狠狠久久av| 国产精品久久视频播放| 国内精品一区二区在线观看| АⅤ资源中文在线天堂| 亚洲第一电影网av| 一级a爱片免费观看的视频| 国产真实乱freesex| 国产精品国产高清国产av| 久久久久久国产a免费观看| 又爽又黄a免费视频| 欧美日本亚洲视频在线播放| 69av精品久久久久久| 内地一区二区视频在线| 联通29元200g的流量卡| 午夜福利成人在线免费观看| 国内精品一区二区在线观看| 国产三级在线视频| 欧美精品国产亚洲| 亚洲欧美中文字幕日韩二区| 国产亚洲精品综合一区在线观看| 久久国产乱子免费精品| 日本a在线网址| 女的被弄到高潮叫床怎么办| 99久久九九国产精品国产免费| 十八禁网站免费在线| 亚洲成a人片在线一区二区| 日韩精品中文字幕看吧| 97超视频在线观看视频| 亚洲成人久久爱视频| 亚洲最大成人中文| 黑人高潮一二区| 亚洲人与动物交配视频| 成人精品一区二区免费| 日产精品乱码卡一卡2卡三| 精品久久久久久久末码| 国产精品人妻久久久久久| 日韩欧美精品v在线| 一区二区三区免费毛片| 亚洲丝袜综合中文字幕| 最新中文字幕久久久久| 美女内射精品一级片tv| 精品一区二区三区av网在线观看| 亚洲精品一卡2卡三卡4卡5卡| 亚洲自偷自拍三级| av女优亚洲男人天堂| 亚洲无线观看免费| 国产精品,欧美在线| 好男人在线观看高清免费视频| 黄色欧美视频在线观看| 日韩一本色道免费dvd| 国产精品爽爽va在线观看网站| 51国产日韩欧美| 亚洲人成网站在线观看播放| 村上凉子中文字幕在线| 日日摸夜夜添夜夜爱| 中文在线观看免费www的网站| 精品一区二区三区视频在线观看免费| 久久久精品94久久精品| 99热这里只有是精品在线观看| 一个人看的www免费观看视频| 亚洲精品色激情综合| 精品人妻熟女av久视频| 尾随美女入室| 搡女人真爽免费视频火全软件 | 日本黄色片子视频| h日本视频在线播放| 在线免费观看的www视频| 能在线免费观看的黄片| 国产探花极品一区二区| 99久久精品热视频| 精品午夜福利视频在线观看一区| 神马国产精品三级电影在线观看| 热99re8久久精品国产| 精品一区二区三区人妻视频| 最近手机中文字幕大全| 国产精品久久久久久久电影| 内射极品少妇av片p| 五月玫瑰六月丁香| av在线老鸭窝| 国产一区亚洲一区在线观看| 久久久久久久久久久丰满| 国产一区二区三区av在线 | 人人妻人人看人人澡| 亚洲熟妇熟女久久| 啦啦啦观看免费观看视频高清| www.色视频.com| 久久草成人影院| 国产精品福利在线免费观看| 国产日本99.免费观看| 1000部很黄的大片| 老司机影院成人| 小说图片视频综合网站| 性色avwww在线观看| 国产精品1区2区在线观看.| 成人午夜高清在线视频| 伦理电影大哥的女人| 一级毛片久久久久久久久女| 亚洲av中文字字幕乱码综合| 成人特级av手机在线观看| 尤物成人国产欧美一区二区三区| 午夜a级毛片| 一级av片app| 婷婷精品国产亚洲av在线| 亚洲国产精品成人久久小说 | 五月玫瑰六月丁香| 免费观看精品视频网站| 亚洲精品粉嫩美女一区| 真人做人爱边吃奶动态| 亚洲最大成人av| 亚洲人成网站在线观看播放| 舔av片在线| 国产欧美日韩精品一区二区| 亚洲精品粉嫩美女一区| 看黄色毛片网站| 国产成人福利小说| 人人妻人人看人人澡| 久久久久久久午夜电影| 91在线观看av| 国产精品女同一区二区软件| 精品少妇黑人巨大在线播放 | 国产欧美日韩精品亚洲av| 国产 一区 欧美 日韩| 九九久久精品国产亚洲av麻豆| 成年av动漫网址| 最新中文字幕久久久久| 国产乱人偷精品视频| 菩萨蛮人人尽说江南好唐韦庄 | 一区二区三区免费毛片| 国产精品亚洲一级av第二区| 高清毛片免费观看视频网站| 久久精品91蜜桃| videossex国产| 亚洲成人av在线免费| 国模一区二区三区四区视频| 插逼视频在线观看| 性色avwww在线观看| 国产亚洲精品久久久久久毛片| 18禁在线无遮挡免费观看视频 | 国产黄色视频一区二区在线观看 | 亚洲国产高清在线一区二区三| 高清毛片免费观看视频网站| 精品一区二区三区人妻视频| 日韩 亚洲 欧美在线| 大型黄色视频在线免费观看| 97超级碰碰碰精品色视频在线观看| 最好的美女福利视频网| 特大巨黑吊av在线直播| 久久久成人免费电影| 成熟少妇高潮喷水视频| 在线观看一区二区三区| 亚洲天堂国产精品一区在线| 十八禁网站免费在线| 大又大粗又爽又黄少妇毛片口| 国产国拍精品亚洲av在线观看| 国产精品一区二区免费欧美| 有码 亚洲区| 精品久久久久久久人妻蜜臀av| 九九热线精品视视频播放| 亚洲欧美成人精品一区二区| a级毛片a级免费在线| 日本一本二区三区精品| 成人av一区二区三区在线看| 男人舔奶头视频| 男女之事视频高清在线观看| 天堂影院成人在线观看| 最好的美女福利视频网| 乱系列少妇在线播放| 欧美人与善性xxx| 亚州av有码| 精品久久久久久成人av| 亚洲精品亚洲一区二区| 天天躁日日操中文字幕| 99精品在免费线老司机午夜| 国产不卡一卡二| 1024手机看黄色片| 国产精品1区2区在线观看.| 激情 狠狠 欧美| 淫妇啪啪啪对白视频| 最新中文字幕久久久久| aaaaa片日本免费| 夜夜看夜夜爽夜夜摸| 卡戴珊不雅视频在线播放| 中文资源天堂在线| 免费不卡的大黄色大毛片视频在线观看 | 免费看a级黄色片| 简卡轻食公司| av在线亚洲专区| 国产成人91sexporn| 麻豆精品久久久久久蜜桃| 在线观看免费视频日本深夜| 综合色av麻豆| 最近最新中文字幕大全电影3| 午夜福利18| 能在线免费观看的黄片| 熟女电影av网| 欧美在线一区亚洲| 91在线观看av| 女生性感内裤真人,穿戴方法视频| 中文在线观看免费www的网站| 色噜噜av男人的天堂激情| 日本黄色片子视频| 91午夜精品亚洲一区二区三区| 97热精品久久久久久| 国产探花极品一区二区| a级毛片a级免费在线| 联通29元200g的流量卡| 丰满人妻一区二区三区视频av| 女的被弄到高潮叫床怎么办| 日本与韩国留学比较| 99在线视频只有这里精品首页| 波野结衣二区三区在线| 一区二区三区四区激情视频 | 久久久a久久爽久久v久久| 婷婷精品国产亚洲av在线| 亚洲最大成人手机在线| 又黄又爽又免费观看的视频| 久久久久久大精品| 可以在线观看毛片的网站| 免费不卡的大黄色大毛片视频在线观看 | 午夜精品在线福利| 观看免费一级毛片| 婷婷六月久久综合丁香| 一卡2卡三卡四卡精品乱码亚洲| 欧美日本亚洲视频在线播放| 午夜爱爱视频在线播放| 国产伦一二天堂av在线观看| 丰满人妻一区二区三区视频av| 干丝袜人妻中文字幕| 国产成人精品久久久久久| 搡老妇女老女人老熟妇| 国产精品电影一区二区三区| 菩萨蛮人人尽说江南好唐韦庄 | av.在线天堂| 男女视频在线观看网站免费| av在线亚洲专区| 国产一区二区在线av高清观看| 黄色欧美视频在线观看| 国产成年人精品一区二区| 欧美另类亚洲清纯唯美| 18禁黄网站禁片免费观看直播| 色视频www国产| 婷婷精品国产亚洲av在线| 国产精品美女特级片免费视频播放器| 看十八女毛片水多多多| 成人午夜高清在线视频| 国产探花极品一区二区| 我要搜黄色片| 欧美+日韩+精品| 欧美性感艳星| 欧美日本视频| 国产精品日韩av在线免费观看| 亚洲天堂国产精品一区在线| 日韩制服骚丝袜av| 欧洲精品卡2卡3卡4卡5卡区| 色噜噜av男人的天堂激情| 桃色一区二区三区在线观看| 精品一区二区三区视频在线观看免费| 色哟哟·www| 久久人妻av系列| 日韩精品有码人妻一区| 老司机午夜福利在线观看视频| 日本五十路高清| 精品人妻熟女av久视频| 国产精品久久久久久久电影| 日本爱情动作片www.在线观看 | 露出奶头的视频| 波多野结衣巨乳人妻| 色哟哟·www| 可以在线观看毛片的网站| 国产欧美日韩精品亚洲av| 黄色欧美视频在线观看| 天堂√8在线中文| 国产精品野战在线观看| 男女边吃奶边做爰视频| 大又大粗又爽又黄少妇毛片口| 国产黄a三级三级三级人| 天堂影院成人在线观看| 我的女老师完整版在线观看| 欧美另类亚洲清纯唯美| 男人舔女人下体高潮全视频| 91精品国产九色| 精品久久久久久久久亚洲| 亚洲精品456在线播放app| 一级a爱片免费观看的视频| 99久国产av精品| 精品久久久久久久久亚洲| 校园春色视频在线观看| 黄片wwwwww| 亚洲不卡免费看| 99热网站在线观看| 观看免费一级毛片| 免费不卡的大黄色大毛片视频在线观看 | 午夜精品在线福利| 国产亚洲精品久久久久久毛片| 一级毛片我不卡| 亚洲经典国产精华液单| 一个人观看的视频www高清免费观看| 国产精品一区二区三区四区免费观看 | 亚洲欧美精品综合久久99| 日本色播在线视频| 晚上一个人看的免费电影| 日韩精品有码人妻一区| 色尼玛亚洲综合影院| 日韩中字成人| 老司机午夜福利在线观看视频| 亚洲精品一区av在线观看| 可以在线观看的亚洲视频| 三级国产精品欧美在线观看| 久久久国产成人免费| 日韩一本色道免费dvd| 免费在线观看成人毛片| 欧美日本亚洲视频在线播放| 国产单亲对白刺激| 亚洲成人中文字幕在线播放| 国产精品久久电影中文字幕| 亚洲综合色惰| 亚洲av成人精品一区久久| 床上黄色一级片| 97超碰精品成人国产| 欧美性感艳星| 看免费成人av毛片| 亚洲精品影视一区二区三区av| 国产熟女欧美一区二区| 99久久精品国产国产毛片| 一区二区三区高清视频在线| 麻豆乱淫一区二区| 午夜影院日韩av| 午夜精品国产一区二区电影 | 成人二区视频| 亚洲国产精品sss在线观看| 久久久国产成人精品二区| 日韩中字成人| 国产乱人偷精品视频| 深夜精品福利| 精品久久久久久久人妻蜜臀av| 99九九线精品视频在线观看视频| 国产成人freesex在线 | 人人妻人人澡人人爽人人夜夜 | 国模一区二区三区四区视频| 日本-黄色视频高清免费观看| 好男人在线观看高清免费视频| 又粗又爽又猛毛片免费看| 国产日本99.免费观看| 少妇被粗大猛烈的视频| 别揉我奶头 嗯啊视频| 黑人高潮一二区| 国产在线男女| 最近中文字幕高清免费大全6| 亚洲专区国产一区二区| 国产亚洲av嫩草精品影院| 12—13女人毛片做爰片一| 久99久视频精品免费| 最好的美女福利视频网| 亚洲久久久久久中文字幕| 亚洲一级一片aⅴ在线观看| 人妻丰满熟妇av一区二区三区| 激情 狠狠 欧美| 国产白丝娇喘喷水9色精品| 男插女下体视频免费在线播放| 嫩草影视91久久| 欧美成人免费av一区二区三区| 看黄色毛片网站| 国产精品久久久久久亚洲av鲁大| 人妻久久中文字幕网| 在线观看一区二区三区| 国产在视频线在精品| 免费高清视频大片| 成人欧美大片| 69人妻影院| 欧美区成人在线视频| 非洲黑人性xxxx精品又粗又长| 97人妻精品一区二区三区麻豆| 亚洲久久久久久中文字幕| 国产av不卡久久| 国内精品一区二区在线观看| 精品久久久久久久末码| 亚洲激情五月婷婷啪啪| 午夜福利成人在线免费观看| 国产国拍精品亚洲av在线观看| 久久久久免费精品人妻一区二区| 日韩欧美精品免费久久| 国产精品综合久久久久久久免费|