• <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
    欧美另类亚洲清纯唯美| 色哟哟·www| 日韩高清综合在线| av黄色大香蕉| 精品不卡国产一区二区三区| 久久婷婷人人爽人人干人人爱| 日本a在线网址| 久久人人爽人人片av| a级毛色黄片| 青春草视频在线免费观看| 中文亚洲av片在线观看爽| 欧美日韩在线观看h| 噜噜噜噜噜久久久久久91| 高清毛片免费观看视频网站| 国产aⅴ精品一区二区三区波| 不卡视频在线观看欧美| 麻豆久久精品国产亚洲av| 免费av毛片视频| 日韩av在线大香蕉| 亚洲国产精品国产精品| 你懂的网址亚洲精品在线观看 | 男女啪啪激烈高潮av片| 中文字幕熟女人妻在线| 国产真实乱freesex| 老师上课跳d突然被开到最大视频| 亚洲av五月六月丁香网| 亚洲美女视频黄频| 国产成人影院久久av| 久久精品国产亚洲av天美| 国产成人freesex在线 | 啦啦啦观看免费观看视频高清| 免费av不卡在线播放| 欧美bdsm另类| 一个人观看的视频www高清免费观看| 一级av片app| 国产伦一二天堂av在线观看| 国产三级在线视频| 91在线精品国自产拍蜜月| 美女高潮的动态| 午夜日韩欧美国产| 欧美又色又爽又黄视频| 久久精品综合一区二区三区| 久久99热6这里只有精品| 一级毛片aaaaaa免费看小| 日本熟妇午夜| 91午夜精品亚洲一区二区三区| 国产欧美日韩精品亚洲av| 久久久色成人| 国产欧美日韩精品一区二区| 午夜福利在线在线| 久久6这里有精品| 91精品国产九色| 内射极品少妇av片p| 村上凉子中文字幕在线| 69av精品久久久久久| 特大巨黑吊av在线直播| 欧美日韩综合久久久久久| 精品一区二区免费观看| 亚洲va在线va天堂va国产| 乱码一卡2卡4卡精品| 国产老妇女一区| 日日撸夜夜添| 国产真实乱freesex| 国产伦精品一区二区三区四那| 两个人的视频大全免费| 嫩草影院精品99| 麻豆乱淫一区二区| 国内精品一区二区在线观看| 日日撸夜夜添| 日本欧美国产在线视频| 18+在线观看网站| 人人妻人人澡欧美一区二区| 一区二区三区四区激情视频 | 两性午夜刺激爽爽歪歪视频在线观看| 国产精品av视频在线免费观看| 欧美成人a在线观看| 国产精品永久免费网站| 色哟哟哟哟哟哟| 床上黄色一级片| 国产精品美女特级片免费视频播放器| 熟妇人妻久久中文字幕3abv| 婷婷色综合大香蕉| 插阴视频在线观看视频| 真人做人爱边吃奶动态| 婷婷色综合大香蕉| 国产亚洲精品久久久com| 国产色婷婷99| av中文乱码字幕在线| 综合色丁香网| 免费不卡的大黄色大毛片视频在线观看 | 一卡2卡三卡四卡精品乱码亚洲| 天天一区二区日本电影三级| 久久中文看片网| 男人和女人高潮做爰伦理| 精品国内亚洲2022精品成人| 午夜福利在线在线| 亚洲精品影视一区二区三区av| 在线播放国产精品三级| 国内少妇人妻偷人精品xxx网站| 男女边吃奶边做爰视频| 一级av片app| 国产黄片美女视频| 十八禁网站免费在线| 乱系列少妇在线播放| 国产亚洲精品综合一区在线观看| 国内揄拍国产精品人妻在线| 国模一区二区三区四区视频| 99久久精品一区二区三区| 身体一侧抽搐| 观看美女的网站| 色av中文字幕| 亚洲成人久久性| 国产一区二区激情短视频| 午夜激情福利司机影院| 黄色日韩在线| 久久久久久久午夜电影| 99热这里只有精品一区| 亚洲成人久久性| 舔av片在线| 最新在线观看一区二区三区| 国产精品1区2区在线观看.| 我的女老师完整版在线观看| 国产美女午夜福利| 看十八女毛片水多多多| 国产一区二区激情短视频| 国产精品一二三区在线看| 少妇熟女aⅴ在线视频| 精品国产三级普通话版| 国产人妻一区二区三区在| 少妇的逼水好多| 色在线成人网| 久久国内精品自在自线图片| 日韩在线高清观看一区二区三区| 婷婷六月久久综合丁香| av在线播放精品| 啦啦啦啦在线视频资源| 三级经典国产精品| eeuss影院久久| 天天躁日日操中文字幕| 麻豆一二三区av精品| 无遮挡黄片免费观看| 午夜福利在线观看吧| 国产精品一区二区免费欧美| 22中文网久久字幕| av视频在线观看入口| 国产av一区在线观看免费| 悠悠久久av| 成人漫画全彩无遮挡| 精品99又大又爽又粗少妇毛片| 熟女电影av网| 超碰av人人做人人爽久久| 99热这里只有精品一区| 国产高清视频在线播放一区| 少妇猛男粗大的猛烈进出视频 | 国内精品美女久久久久久| 久久久久久大精品| 级片在线观看| 日日啪夜夜撸| 91久久精品国产一区二区成人| 尤物成人国产欧美一区二区三区| 人妻夜夜爽99麻豆av| 最近视频中文字幕2019在线8| 欧美中文日本在线观看视频| 不卡一级毛片| 69av精品久久久久久| 亚洲在线自拍视频| 精品久久久久久久久亚洲| 国产精华一区二区三区| 欧美日本视频| 国产成年人精品一区二区| 少妇人妻精品综合一区二区 | 人妻夜夜爽99麻豆av| 偷拍熟女少妇极品色| 高清毛片免费看| 成人毛片a级毛片在线播放| 久久综合国产亚洲精品| 天堂网av新在线| 欧美日本亚洲视频在线播放| 亚洲欧美日韩高清专用| 白带黄色成豆腐渣| 亚洲内射少妇av| 免费看日本二区| 插逼视频在线观看| 成人午夜高清在线视频| 日韩一区二区视频免费看| 美女大奶头视频| 国产亚洲av嫩草精品影院| 亚洲成av人片在线播放无| 日韩精品有码人妻一区| 久久精品夜色国产| 免费无遮挡裸体视频| 变态另类丝袜制服| 亚洲图色成人| 亚洲精品一区av在线观看| 国产成人a∨麻豆精品| 免费人成在线观看视频色| 97人妻精品一区二区三区麻豆| 一本久久中文字幕| 黄色一级大片看看| 亚洲精品国产成人久久av| 亚洲人成网站在线播| 悠悠久久av| 亚洲aⅴ乱码一区二区在线播放| 成人美女网站在线观看视频| 亚洲欧美日韩无卡精品| 成人性生交大片免费视频hd| 插逼视频在线观看| 免费看美女性在线毛片视频| 久久久久性生活片| av在线亚洲专区| 久久99热这里只有精品18| 久久精品夜色国产| 熟女电影av网| 少妇裸体淫交视频免费看高清| 日本一本二区三区精品| 久久精品久久久久久噜噜老黄 | 高清毛片免费观看视频网站| 欧美激情久久久久久爽电影| 乱人视频在线观看| 美女被艹到高潮喷水动态| 最新中文字幕久久久久| 夜夜爽天天搞| 99热这里只有精品一区| 成人特级黄色片久久久久久久| 俄罗斯特黄特色一大片| 成人欧美大片| 亚洲美女视频黄频| aaaaa片日本免费| 成人亚洲欧美一区二区av| 欧美性猛交黑人性爽| 国产乱人偷精品视频| 亚洲精品一区av在线观看| 国产精华一区二区三区| av福利片在线观看| 日日干狠狠操夜夜爽| 日日摸夜夜添夜夜添小说| 精品国产三级普通话版| 男人狂女人下面高潮的视频| 午夜福利在线在线| 乱系列少妇在线播放| 波野结衣二区三区在线| 国产高清有码在线观看视频| 亚洲无线观看免费| 久久精品综合一区二区三区| 啦啦啦观看免费观看视频高清| 男人舔女人下体高潮全视频| 成人鲁丝片一二三区免费| 久久精品夜夜夜夜夜久久蜜豆| 春色校园在线视频观看| 日韩精品中文字幕看吧| 欧美一区二区亚洲| 日韩欧美在线乱码| 床上黄色一级片| 久久精品国产鲁丝片午夜精品| 久久久午夜欧美精品| 亚洲精品国产av成人精品 | 久久6这里有精品| 真人做人爱边吃奶动态| 天堂√8在线中文| 亚洲久久久久久中文字幕| 男人舔奶头视频| 性色avwww在线观看| 日本熟妇午夜| 我要搜黄色片| 天堂av国产一区二区熟女人妻| 精品一区二区免费观看| 内地一区二区视频在线| 国产一区二区在线观看日韩| a级毛色黄片| 我的老师免费观看完整版| 精品久久久久久久末码| 国产精品一区二区三区四区免费观看 | 国产精品永久免费网站| 日韩精品有码人妻一区| 老司机福利观看| 神马国产精品三级电影在线观看| 成人av在线播放网站| 亚洲成人av在线免费| 狠狠狠狠99中文字幕| 久久久久久久久中文| 免费电影在线观看免费观看| 成人亚洲精品av一区二区| 精品熟女少妇av免费看| 精品久久国产蜜桃| 老司机午夜福利在线观看视频| 成人高潮视频无遮挡免费网站| 久久精品久久久久久噜噜老黄 | 国产私拍福利视频在线观看| 国产精品永久免费网站| 亚洲美女视频黄频| 国产精品久久电影中文字幕| 国产三级在线视频| 国产精品一二三区在线看| 老熟妇仑乱视频hdxx| 亚洲精品日韩在线中文字幕 | 蜜桃久久精品国产亚洲av| 人人妻,人人澡人人爽秒播| 亚洲国产色片| 色哟哟哟哟哟哟| 国产精品久久久久久av不卡| 国产在线精品亚洲第一网站| 国产乱人视频| 亚洲av熟女| 国产久久久一区二区三区| 欧美日韩一区二区视频在线观看视频在线 | 少妇熟女欧美另类| 欧美不卡视频在线免费观看| 色综合亚洲欧美另类图片| 成人无遮挡网站| 欧美最黄视频在线播放免费| 搡老岳熟女国产| videossex国产| 亚洲丝袜综合中文字幕| 丰满人妻一区二区三区视频av| 99热这里只有精品一区| 国内揄拍国产精品人妻在线| 蜜桃久久精品国产亚洲av| 国产 一区精品| 午夜激情福利司机影院| 日日啪夜夜撸| 麻豆乱淫一区二区| 国产精品福利在线免费观看| 免费在线观看成人毛片| 波多野结衣高清无吗| 男插女下体视频免费在线播放| 看免费成人av毛片| 成人午夜高清在线视频| 欧美不卡视频在线免费观看| 男女视频在线观看网站免费| 亚洲国产精品合色在线| 99国产精品一区二区蜜桃av| 亚洲色图av天堂| 最近中文字幕高清免费大全6| 国产午夜精品论理片| 露出奶头的视频| 国产激情偷乱视频一区二区| 亚洲无线在线观看| 18+在线观看网站| 天堂√8在线中文| 精品一区二区三区av网在线观看| 99热网站在线观看| 国产精品一二三区在线看| 国产一级毛片七仙女欲春2| 午夜福利在线观看吧| 亚洲av中文字字幕乱码综合| 国产成人91sexporn| 久久精品国产亚洲av天美| 最新中文字幕久久久久| 国产久久久一区二区三区| 日韩精品有码人妻一区| 国产乱人偷精品视频| 久久婷婷人人爽人人干人人爱| 精品久久久久久成人av| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲五月天丁香| 久久久久国产网址| 男女视频在线观看网站免费| 在线a可以看的网站| 久久午夜福利片| 亚洲av二区三区四区| 卡戴珊不雅视频在线播放| 国产精品一区二区三区四区免费观看 | 99热精品在线国产| 久久久久久久午夜电影| 日日摸夜夜添夜夜添av毛片| 午夜亚洲福利在线播放| 日韩国内少妇激情av| 欧美日韩乱码在线| 国内精品一区二区在线观看| 六月丁香七月| 99视频精品全部免费 在线| 性色avwww在线观看| 欧美区成人在线视频| 国产成人影院久久av| 在线天堂最新版资源| 午夜老司机福利剧场| 中文亚洲av片在线观看爽| 男人的好看免费观看在线视频| 在线观看66精品国产| 亚洲婷婷狠狠爱综合网| 免费观看在线日韩| 最近在线观看免费完整版| 麻豆久久精品国产亚洲av| 精品一区二区三区视频在线| 日韩av不卡免费在线播放| 午夜激情福利司机影院| 国产真实伦视频高清在线观看| 久久精品国产99精品国产亚洲性色| 色在线成人网| 亚洲久久久久久中文字幕| 热99在线观看视频| 精品久久久久久久久久久久久| 香蕉av资源在线| 国产色爽女视频免费观看| 亚洲18禁久久av| 国内精品久久久久精免费| 波多野结衣高清作品| 精品久久久噜噜| 欧美精品国产亚洲| 亚洲人成网站在线观看播放| 美女xxoo啪啪120秒动态图| 精品少妇黑人巨大在线播放 | 淫秽高清视频在线观看| 国产免费男女视频| 免费人成视频x8x8入口观看| 全区人妻精品视频| 能在线免费观看的黄片| 此物有八面人人有两片| 欧美丝袜亚洲另类| 欧美日韩精品成人综合77777| 精品99又大又爽又粗少妇毛片| 久久婷婷人人爽人人干人人爱| 午夜久久久久精精品| 国产精品一及| 免费观看精品视频网站| 国产精品爽爽va在线观看网站| 国产成人一区二区在线| 免费搜索国产男女视频| 夜夜看夜夜爽夜夜摸| 好男人在线观看高清免费视频| 日韩人妻高清精品专区| 免费看日本二区| 国产高清激情床上av| 又爽又黄a免费视频| 亚洲一级一片aⅴ在线观看| 精品久久久噜噜| 欧美一区二区国产精品久久精品| 欧美激情国产日韩精品一区| 日韩av在线大香蕉| 亚洲第一区二区三区不卡| 麻豆国产av国片精品| 久久欧美精品欧美久久欧美| 老司机午夜福利在线观看视频| 久久亚洲精品不卡| 国产高清三级在线| 露出奶头的视频| 两性午夜刺激爽爽歪歪视频在线观看| 欧美成人一区二区免费高清观看| 国产女主播在线喷水免费视频网站 | 白带黄色成豆腐渣| 国内精品宾馆在线| 神马国产精品三级电影在线观看| 女的被弄到高潮叫床怎么办| 丰满人妻一区二区三区视频av| 免费人成在线观看视频色| 久久久精品欧美日韩精品| 国产成人freesex在线 | 亚洲精品在线观看二区| 晚上一个人看的免费电影| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲国产色片| 午夜福利18| 99热这里只有精品一区| 五月玫瑰六月丁香| 日韩欧美在线乱码| 免费搜索国产男女视频| 我的老师免费观看完整版| 欧美色欧美亚洲另类二区| 婷婷色综合大香蕉| 丝袜喷水一区| a级毛片a级免费在线| 成人高潮视频无遮挡免费网站| 国产欧美日韩精品亚洲av| 亚洲最大成人av| 久久久久久久亚洲中文字幕| 亚洲精品色激情综合| 十八禁国产超污无遮挡网站| 国产精品人妻久久久影院| 乱码一卡2卡4卡精品| 尤物成人国产欧美一区二区三区| 91av网一区二区| 91精品国产九色| 亚洲精品成人久久久久久| 性欧美人与动物交配| 六月丁香七月| 亚洲经典国产精华液单| 午夜福利在线观看吧| 亚洲成人av在线免费| 日韩欧美精品免费久久| 精品熟女少妇av免费看| 亚洲最大成人手机在线| 亚洲第一区二区三区不卡| 男人和女人高潮做爰伦理| 成人av一区二区三区在线看| 日韩亚洲欧美综合| 国产黄片美女视频| 免费观看人在逋| 国产午夜精品论理片| 18禁在线播放成人免费| 一个人看视频在线观看www免费| 亚洲av第一区精品v没综合| 亚洲欧美日韩东京热| 最近中文字幕高清免费大全6| 国语自产精品视频在线第100页| 中国美女看黄片| 99视频精品全部免费 在线| 男人和女人高潮做爰伦理| 插逼视频在线观看| 男插女下体视频免费在线播放| 波野结衣二区三区在线| 夜夜爽天天搞| 又黄又爽又免费观看的视频| 美女xxoo啪啪120秒动态图| 少妇高潮的动态图| 男人和女人高潮做爰伦理| 人妻夜夜爽99麻豆av| 美女内射精品一级片tv| 亚洲欧美日韩东京热| 校园人妻丝袜中文字幕| 国产成人freesex在线 | 国产免费男女视频| 国产精华一区二区三区| 99国产极品粉嫩在线观看| 日本三级黄在线观看| 有码 亚洲区| 黄色一级大片看看| 国产人妻一区二区三区在| 搡老岳熟女国产| 国产高清三级在线| 在线观看美女被高潮喷水网站| 国产亚洲精品久久久久久毛片| 麻豆乱淫一区二区| 久久精品国产清高在天天线| av专区在线播放| 男女边吃奶边做爰视频| 91在线观看av| 18禁在线无遮挡免费观看视频 | aaaaa片日本免费| 自拍偷自拍亚洲精品老妇| 91精品国产九色| 免费高清视频大片| 欧美激情久久久久久爽电影| 简卡轻食公司| 午夜a级毛片| a级一级毛片免费在线观看| 少妇猛男粗大的猛烈进出视频 | 露出奶头的视频| 国产高清不卡午夜福利| 麻豆精品久久久久久蜜桃| www.色视频.com| 亚洲成av人片在线播放无| 午夜老司机福利剧场| 人人妻人人澡欧美一区二区| 丝袜喷水一区| 最近中文字幕高清免费大全6| 国产免费一级a男人的天堂| 特大巨黑吊av在线直播| 少妇的逼好多水| 看免费成人av毛片| 久99久视频精品免费| 热99在线观看视频| 18+在线观看网站| 成人午夜高清在线视频| ponron亚洲| 国产久久久一区二区三区| 综合色av麻豆| 中文亚洲av片在线观看爽| 成人鲁丝片一二三区免费| 国产精品一区二区三区四区久久| 久久久久国产精品人妻aⅴ院| 亚洲国产色片| 草草在线视频免费看| 日韩欧美在线乱码| 俺也久久电影网| 国产午夜精品久久久久久一区二区三区 | 伦理电影大哥的女人| 国产人妻一区二区三区在| 男女做爰动态图高潮gif福利片| 国产在线男女| 午夜老司机福利剧场| 两个人的视频大全免费| 麻豆乱淫一区二区| 久久亚洲精品不卡| 国产精品电影一区二区三区| av卡一久久| 最好的美女福利视频网| 午夜精品一区二区三区免费看| 日本爱情动作片www.在线观看 | 如何舔出高潮| 国产精品日韩av在线免费观看| 一级黄色大片毛片| 国产69精品久久久久777片| 国产亚洲91精品色在线| 白带黄色成豆腐渣| 国产日本99.免费观看| 精品乱码久久久久久99久播| 少妇熟女aⅴ在线视频| 一进一出抽搐动态| 国产精品一及| 少妇熟女aⅴ在线视频| 狂野欧美白嫩少妇大欣赏| 美女黄网站色视频| 日韩,欧美,国产一区二区三区 | 亚洲最大成人手机在线| 国产不卡一卡二| 久久久欧美国产精品| 给我免费播放毛片高清在线观看| 国产不卡一卡二| 在线播放无遮挡| 好男人在线观看高清免费视频| 国产高清激情床上av| 国产午夜精品论理片| 在线免费十八禁| 男女做爰动态图高潮gif福利片| 国产成人一区二区在线| 别揉我奶头~嗯~啊~动态视频| 国产成人91sexporn| 国产成人freesex在线 | 搡女人真爽免费视频火全软件 | 女生性感内裤真人,穿戴方法视频| 国产精品1区2区在线观看.| 成人鲁丝片一二三区免费| 男女边吃奶边做爰视频| 国产精品久久久久久久电影| 午夜福利成人在线免费观看|