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

    橫向氣流中非牛頓液體射流直接數(shù)值模擬

    2016-12-06 07:06:57朱呈祥尤延鋮
    航空學(xué)報(bào) 2016年9期
    關(guān)鍵詞:牛頓流體氣液黏性

    朱呈祥*,尤延鋮

    廈門(mén)大學(xué) 航空航天學(xué)院,廈門(mén) 361005

    橫向氣流中非牛頓液體射流直接數(shù)值模擬

    朱呈祥*,尤延鋮

    廈門(mén)大學(xué) 航空航天學(xué)院,廈門(mén) 361005

    采用直接數(shù)值模擬研究了動(dòng)量比為6的非牛頓液體射流在橫向氣流中的破碎特征,重點(diǎn)分析了表面波發(fā)展、彎曲特性和展向擴(kuò)散等射流結(jié)構(gòu)及其非牛頓特征。在非牛頓液體射流的近噴嘴位置,其表面跡線隨時(shí)間擺動(dòng),并有逆橫向氣流方向運(yùn)動(dòng)的趨勢(shì)。射流的展向擴(kuò)散僅在噴入初期較為明顯,此后基本維持在接近35°的擴(kuò)散半角。射流的一次破碎,尤其是液絲與液滴、甚而衛(wèi)星液滴的生成過(guò)程都被精細(xì)地捕捉和描述。在近噴嘴的射流柱附近,橫向氣流流動(dòng)具有一定的圓柱擾流相似性,在其他區(qū)域則表現(xiàn)為極其復(fù)雜的紊流特征。射流的非牛頓特性主要體現(xiàn)在黏性系數(shù),不同位置的液體黏性系數(shù)相差超過(guò)20%。相較牛頓液體射流,剪切稀化的非牛頓射流具有更易破碎的特征。

    液體射流;橫向氣流;非牛頓流體;破碎;直接數(shù)值模擬

    液體射流噴入橫向氣流被廣泛應(yīng)用在工程實(shí)際中,包括發(fā)動(dòng)機(jī)的液體燃料噴注、農(nóng)業(yè)灌溉、噴墨打印等,深入理解液體射流的物理本質(zhì)對(duì)開(kāi)發(fā)高效率的此類(lèi)設(shè)備至關(guān)重要。

    目前國(guó)際上針對(duì)液體射流噴入橫向氣流已經(jīng)開(kāi)展了大量研究。在實(shí)驗(yàn)領(lǐng)域,Lee[1]和 Sallam[2]等分別對(duì)比了不同湍動(dòng)進(jìn)口條件下和非湍動(dòng)進(jìn)口時(shí)圓形液體射流在橫向氣流中的噴入特征;Birouk[3]、Wu[4]和Stenzler[5]等則進(jìn)一步整理了液體射流噴入橫向氣流后的發(fā)展規(guī)律,包括射流的破碎模式、噴射軌跡、穿透深度、破碎長(zhǎng)度等。此外,Gutmark[6]、Shapiro[7]、Megerian[8]和 Coletti[9]等針對(duì)射流在橫向氣流中的破碎特征也都開(kāi)展了相關(guān)的實(shí)驗(yàn)研究。

    在數(shù)值領(lǐng)域,Herrmann[10-11]以力平衡優(yōu)化水平集(Balanced Force Refined Level Set)方法為基礎(chǔ),重點(diǎn)對(duì)比了不同氣液兩相密度比下的射流破碎特征,Cavar和 Meyer[12]則運(yùn)用大渦模擬(LES)技術(shù)探究了射流的渦場(chǎng)結(jié)構(gòu)。與Cavar和Meyer的研究類(lèi)似,Galeazzo等[13]也主要關(guān)心渦的形成發(fā)展,但其重心偏向于氣液兩相的摻混。在針對(duì)液體射流噴入橫向氣流的參數(shù)研究中,Sau和Mahesh[14]根據(jù)霧化水平專(zhuān)門(mén)對(duì)液體射流的噴入條件進(jìn)行了優(yōu)化,Pai等[15]則將重心放在流量比和氣流韋伯?dāng)?shù)上,比較了它們對(duì)射流特征的影響,Muldoon和Acharya[16]則從力學(xué)角度出發(fā),分析了不同簡(jiǎn)諧力下橫向氣流中的射流結(jié)構(gòu)。此外,Margason[17]、Aalburg[18]和 Mahesh[19]等還分別對(duì)目前該領(lǐng)域的國(guó)際前沿工作進(jìn)行了整理總結(jié)。

    前文提到的所有研究都只考慮了牛頓流體。在非牛頓流體方面,Wong等[20]重點(diǎn)分析了4種不同的射流破碎類(lèi)別,并比較了最終的破碎長(zhǎng)度和液滴尺寸。Clasen等[21]則對(duì)一類(lèi)黏彈性流體開(kāi)展了實(shí)驗(yàn)研究,并著重分析了它的稀化特性。此外,Yarin[22]還專(zhuān)門(mén)針對(duì)非牛頓流體的力學(xué)和流變特性進(jìn)行了總結(jié)。然而到目前為止,國(guó)際上還沒(méi)有針對(duì)非牛頓流體開(kāi)展過(guò)射流在橫向氣流中的研究。但在工程實(shí)際問(wèn)題中,非牛頓流體其實(shí)廣泛存在,以液體燃料為例,通常將其簡(jiǎn)單處理為牛頓流體,但事實(shí)上在添加了穩(wěn)定劑和抗氧化劑后,液體燃料是具有弱的非牛頓特性的。而在液體火箭發(fā)動(dòng)機(jī)中,甚至需要采用高黏度的剪切稀化非牛頓流體作為燃料,以應(yīng)對(duì)它的儲(chǔ)存問(wèn)題。這些將直接導(dǎo)致其破碎體系相較傳統(tǒng)煤油更為復(fù)雜,并最終誘發(fā)不同的霧化效果與燃燒效率。

    因此,本文將以非牛頓流體為研究對(duì)象,開(kāi)展直接數(shù)值模擬(Direct Numerical Simulation,DNS)研究,重點(diǎn)分析液體射流在近噴嘴附近的流動(dòng)特征。介紹采用的數(shù)值方法與設(shè)置,重點(diǎn)分析液體射流的三維結(jié)構(gòu)、彎曲和擴(kuò)散特性、時(shí)域變化、渦場(chǎng)特征等,討論其非牛頓黏性變化,并與牛頓射流進(jìn)行對(duì)比。

    1 數(shù)值方法

    本文采用的數(shù)值工具為自主開(kāi)發(fā)的DNS程序Free Surface 3D(FS3D),該程序求解的是三Δ維不可壓Navier-Stokes方程組:式中:u為速度矢量;ρ為密度;p為壓力;t為時(shí)間;k為外部作用力;T為氣液兩相分界面處的表面張力。FS3D程序采用流體體積(Volume of Fluid,VOF)[23]方法捕捉氣液兩相分界面,該方法定義變量f表征單元格內(nèi)的液體體積分?jǐn)?shù),即

    f滿足以下守恒關(guān)系:

    為了精確描述氣液兩相分界面,F(xiàn)S3D程序還運(yùn)用分段線性界面計(jì)算(Piecewise Linear Interface Calculation,PLIC)[24]方法進(jìn)行了分界面重構(gòu)。此外,模擬射流進(jìn)口的湍流度對(duì)工程應(yīng)用同樣至關(guān)重要,因此FS3D程序還搭建了一個(gè)基于Kornev和Hassel[25]的來(lái)流發(fā)生器。相應(yīng)的數(shù)值方法均已在文獻(xiàn)[26-27]中進(jìn)行了氣液兩相液滴和瑞利破碎射流的實(shí)驗(yàn)驗(yàn)證,也說(shuō)明了方法的可靠與準(zhǔn)確性。

    對(duì)于非牛頓黏性,本文采用以下的冪律函數(shù)進(jìn)行模擬:

    式中:μ0為零剪切時(shí)的動(dòng)力黏度為剪切率;K和n為取決于流體和環(huán)境的模型常數(shù)。該冪律模型已被 Motzigemba[28]以及 Focke等[29]驗(yàn)證過(guò)。同時(shí),文獻(xiàn)[30]也就液態(tài)射流的黏性系數(shù)開(kāi)展了基于Schroeder等[31]的實(shí)驗(yàn)數(shù)據(jù)的驗(yàn)證。

    本文模擬的橫向氣流為空氣,液體是20%質(zhì)量分?jǐn)?shù)的剪切稀化PVP溶液。該液體的Deborah數(shù) De和Elasto-capillary數(shù)Ec都在10-8量級(jí),遠(yuǎn)低于黏彈性流體的極限值0.35和2.35,因此是典型的冪律流體。

    表1給出了液體與氣體的相關(guān)物性參數(shù)以及射流直徑D,下標(biāo)l和c分別代表液體與橫向氣流。在橫向氣流射流噴注研究中,氣液兩相的動(dòng)量比往往是決定射流流態(tài)的一個(gè)關(guān)鍵參數(shù)[32],本文定義氣液兩相動(dòng)量比q為

    此外,表2還給出了基于式(8)的氣液兩相韋伯?dāng)?shù)We與雷諾數(shù)Re:

    表1 計(jì)算條件設(shè)置Table 1 Computational condition setup

    表2 無(wú)量綱參數(shù)Table 2 Dimensionless parameters

    本文將著重分析射流在啟動(dòng)過(guò)程中近噴嘴附近的流動(dòng)特征。這主要有兩個(gè)方面的考慮,一是對(duì)于持續(xù)噴注,射流初期的流動(dòng)特征對(duì)最終的流動(dòng)形態(tài)影響顯著,這對(duì)透徹理解射流破碎過(guò)程非常重要;另一方面,射流啟動(dòng)過(guò)程也是間歇噴注的主要流動(dòng)特征,因此分析啟動(dòng)過(guò)程是研究間歇噴注的關(guān)鍵。為了研究射流的橫向噴注,本文將采用如圖1所示的矩形計(jì)算域,左側(cè)邊界面為氣流入口,以恒定速度uc吹入橫向來(lái)流,液體從xOy面內(nèi)的圓形噴嘴噴入,橫向氣流沿x正方向從左側(cè)的yOz面流入。下邊界為無(wú)滑移壁面,靠近左側(cè)邊界面開(kāi)有一個(gè)圓形射流入口,液體以帽狀速度型讀入并沿y軸正向噴入,其余邊界設(shè)置為自由出流(von Neumann)條件。對(duì)于左側(cè)的速度入口邊界,直接將來(lái)流速度型賦給虛網(wǎng)格,不考慮液體在噴嘴內(nèi)由于近壁無(wú)滑移而形成的邊界層,在工程中這對(duì)應(yīng)短直噴嘴。該模擬過(guò)程相當(dāng)于先單獨(dú)計(jì)算噴嘴流動(dòng),將其出口速度型提取出來(lái),再賦值給本計(jì)算域直接作為邊界條件讀入。該手段由于無(wú)需布置過(guò)多網(wǎng)格在噴嘴近壁邊界層內(nèi),不會(huì)導(dǎo)致計(jì)算資源浪費(fèi),因此在DNS中已被廣泛采用。為了能夠求解湍流中Kolmogorov長(zhǎng)度[33]以滿足DNS對(duì)網(wǎng)格精度的要求,本文在計(jì)算域10D×10D×5D內(nèi)布置了512×512×256的結(jié)構(gòu)化網(wǎng)格,因此最小網(wǎng)格尺度僅為4μm。

    圖1 矩形計(jì)算域示意圖Fig.1 Schematic of rectangular computational domain

    2 射流結(jié)構(gòu)

    圖2給出了非牛頓液體噴入橫向氣流后近噴嘴處的流動(dòng)結(jié)構(gòu)。圖中:t*=t ul/R,t為物理時(shí)間,ul為射流速度,R為射流半徑,射流表面顏色代表無(wú)量綱速度(為當(dāng)?shù)厮俣瘸砸后w射流速度)。射流沿著y軸正向噴入,橫向氣流沿x軸正向噴入。受橫向氣流的影響,液態(tài)射流沿x軸正向彎曲,并沿z軸展向擴(kuò)散。隨著射流深度的增加,液體被逐漸剪切成薄層,并在一定范圍內(nèi)與橫向氣流的流動(dòng)方向幾乎平行。然而這種平行流動(dòng)的趨勢(shì)在射流前部并不存在,到射流頭部液體甚而逆著橫向氣流的方向彎曲,形成了一個(gè)大的類(lèi)似bag類(lèi)型的破碎,并產(chǎn)生了若干分裂的液絲與液滴。而這些液絲與液滴由于分裂位置不同也表現(xiàn)出了不同的速度分布,靠近噴嘴出口的速度低,靠近頭部前緣的速度高。

    圖2 t*=9.2時(shí)刻的橫向氣流射流結(jié)構(gòu)Fig.2 Flow structure of jet in crossflow(t*=9.2)

    雖然 Hermann[10]也曾在相似的動(dòng)量比 (q=6.6)下開(kāi)展過(guò)橫向射流數(shù)值計(jì)算,但本文發(fā)現(xiàn)的上述流動(dòng)特征與其存在顯著區(qū)別。這主要有兩個(gè)原因:一是無(wú)量綱參數(shù)的差別,二是本文選用了非牛頓剪切稀化流體。在無(wú)量綱參數(shù)方面,盡管Her-mann教授與本文選用的q相似,但他采用的Re和We卻很高。Pai等[15]曾指出,即使采用相同的q,液體韋伯?dāng)?shù)Wel的不同同樣會(huì)產(chǎn)生不同的射流結(jié)構(gòu)。在非牛頓流體和牛頓流體方面,液體的零剪切黏性系數(shù)較高將削弱射流的破碎;在射流內(nèi)不同剪切率的位置,液體由于黏性系數(shù)不同也會(huì)表現(xiàn)出不同的破碎形式,射流的破碎也將最先發(fā)生在剪切率較高的區(qū)域。在以上兩個(gè)原因共同作用下,本文研究的液體射流形成了前文提及的類(lèi)似bag類(lèi)型的流動(dòng)與破碎特征。

    為了更好地觀察橫向氣流中非牛頓射流的結(jié)構(gòu),圖3分別給出了液體的前視與俯視圖,圖中紅色箭頭代表射流方向,藍(lán)色箭頭代表橫向氣流方向,黑色實(shí)線為射流的彎曲軌跡,并與 Wu[4](紅色實(shí)線)和Stenzler[5]等(藍(lán)色實(shí)線)進(jìn)行對(duì)比。從圖3(a)可以更明顯地看出射流的彎曲特性,但與 Wu[4]和Stenzler[5]等擬合的彎曲曲線相比,本文研究的射流彎曲軌跡明顯偏低。Wu在擬合曲線的過(guò)程中主要考慮了橫向射流的力學(xué)特性,在多種牛頓流體(水、乙醇、丙三醇等)和寬參數(shù)范圍(如動(dòng)量比3.38<q<185)基礎(chǔ)上推導(dǎo)射流的彎曲軌跡。他指出,射流的彎曲軌跡主要由兩個(gè)因素決定,一個(gè)是動(dòng)量比q,另一個(gè)就是液體的黏性系數(shù)。在他們的研究中,液體最高的黏性系數(shù)約為水的3.66倍,這明顯低于本文研究的25倍水黏性系數(shù)。Wu在文獻(xiàn)[4]中也對(duì)比了不同黏性系數(shù)流體的橫向射流,他強(qiáng)調(diào),高黏性系數(shù)液體射流會(huì)比低黏性系數(shù)液體射流的彎曲程度更高,這正與本文得到的射流結(jié)構(gòu)是相吻合的。當(dāng)然,由于本文采用的是剪切稀化非牛頓流體,而Wu所研究的均為牛頓流體,因此Wu所提到的高黏性效應(yīng)會(huì)被一定程度削弱。此外,文中僅考慮射流啟動(dòng)過(guò)程的彎曲特性,這與 Wu[4]和Stenzler[5]等的持續(xù)噴注也存在區(qū)別。

    圖3 射流結(jié)構(gòu)的前視圖和俯視圖Fig.3 Front view and top view of jet structure

    從圖3還可以發(fā)現(xiàn),射流不僅沿x軸方向彎曲,而且在y軸方向也存在彎曲,以軸向位置x/D=2D處為例,液體向下彎曲(y軸負(fù)方向)的高度Hben達(dá)到了0.875D。在圖3(b)中,射流沿z方向的展向擴(kuò)散現(xiàn)象也非常明顯。隨著軸向位置的增加,液體的寬度Wjet也逐漸增加,x/D=2.9D處射流的寬度達(dá)到了3.5D。在該圖中,前文提到的射流頭部逆橫向氣流現(xiàn)象也更為顯著,射流逆向區(qū)域的長(zhǎng)度Lrev(定義為下游未破碎區(qū)到逆向頭部前緣的x軸的距離)達(dá)到了1.875D。結(jié)合圖3還可以發(fā)現(xiàn),射流破碎主要發(fā)生在薄層區(qū)(也稱(chēng)為表面破碎),這是由于該區(qū)域內(nèi)液體被最大程度地剪切,因此黏性系數(shù)也最低。射流在該時(shí)刻形成的液滴尺寸覆蓋了D/8~D/40之間,其中圓形液滴的尺寸主要集中在D/10~D/15。射流分裂的位置不同,各小液滴的運(yùn)動(dòng)速度也不盡相同,最快的液滴速度甚至可以超過(guò)最慢液滴的4倍。

    圖4為射流的總面積S*與總質(zhì)量m*隨時(shí)間的變化關(guān)系,這里的面積與質(zhì)量均為無(wú)量綱參數(shù),無(wú)量綱關(guān)系為

    式中:射流的表面積總體上隨時(shí)間增加,具體可劃分為3個(gè)主要區(qū)間。t*<t*1=3時(shí),由于射流僅輕微波動(dòng),表面積緩慢增加;t*>t*1時(shí),射流開(kāi)始出現(xiàn)彎曲并被剪切成薄層,因此表面積快速增加,到t*2=8時(shí)表面積已經(jīng)達(dá)到了5;當(dāng)t*>t*2時(shí),面積增加開(kāi)始放緩,此時(shí)液絲與液滴的生成是面積增加的主要因素,但這種面積的增加相較液體薄層引起的面積增加是少量的。從圖4(b)可以發(fā)現(xiàn),當(dāng)t*>t*3=8.6時(shí)液體質(zhì)量出現(xiàn)突降,這是由于部分液體已經(jīng)流出了計(jì)算域,而在此之前,液體的質(zhì)量幾乎恒定。

    圖4 射流的總面積與總質(zhì)量隨時(shí)間的變化Fig.4 Temporal evolution of total surface area and total mass of jet

    3 特征分析

    3.1 彎曲特性

    液體射流在橫向氣流的噴注研究中,彎曲特性是其最重要的特征之一。圖5給出了射流在7個(gè)典型截面內(nèi)的形狀,各橫截面法向均與射流表面平行,內(nèi)部紅色區(qū)域?yàn)橐后w。在噴嘴出口位置,射流呈圓形,但到截面b位置時(shí),迎風(fēng)面呈平面狀,這種現(xiàn)象在以往的研究中都未被發(fā)現(xiàn)過(guò)。本文將其歸因于液體的高黏性系數(shù),而與非牛頓流體無(wú)關(guān)。為了證明這一點(diǎn),本文在3.5節(jié)也開(kāi)展了高黏性系數(shù)的牛頓流體計(jì)算,發(fā)現(xiàn)其同樣存在平面狀迎風(fēng)面。在以往的實(shí)驗(yàn)研究中,如文獻(xiàn)[4-5],受實(shí)驗(yàn)設(shè)置的限制,很難觀測(cè)到射流的瞬時(shí)橫截面形狀;而在數(shù)值研究中,如文獻(xiàn)[10,18],由于液體的高黏性系數(shù)會(huì)引起氣液兩相分界面處物性參數(shù)的急劇跳動(dòng),因此也很難做到該類(lèi)數(shù)值模擬。當(dāng)然,這種前平后凸的橫截面形狀仍有待通過(guò)以后的實(shí)驗(yàn)去進(jìn)一步驗(yàn)證。

    在下游截面c,由于邊緣受到的剪切較中心區(qū)域高,因此射流的迎風(fēng)面形狀逐漸彎曲,這種特征在截面d內(nèi)更為明顯。在截面e~g內(nèi),射流被逐漸拉成薄層,而且由于射流破碎,薄層的寬度逐漸減小。值得注意的是,在截面f與g內(nèi),受橫向氣流作用影響,液體薄層的邊緣向下彎曲,極大地影響非牛頓液體的黏性,這將在3.5節(jié)中進(jìn)一步討論。

    在研究射流的彎曲特性時(shí),除了單一時(shí)刻(t*=9.2)橫截面內(nèi)的射流形狀,本文還對(duì)比了近噴嘴處的射流彎曲軌跡,如圖6所示。在t*=3時(shí),射流向下游略彎曲,受橫向氣流影響,在射流前緣形成了一個(gè)小的頭部。到t*=5時(shí),該頭部被拉長(zhǎng)并沿x軸方向彎曲,形成了一個(gè)表面凹坑。該凹坑顯然對(duì)橫向氣流更為敏感,因此在t*=7時(shí)液體被進(jìn)一步拉伸剪切。在t*=7和t*=9.2之間,射流頭部被快速拉長(zhǎng),并形成了獨(dú)特的“甩尾”現(xiàn)象,該過(guò)程極為迅速,從t*=7~9.2,射流的軸向長(zhǎng)度增加了一倍,即從2.7D~5.7D。從圖6中還可以發(fā)現(xiàn),射流頭部受氣流作用已出現(xiàn)破碎。

    圖5 t*=9.2時(shí)刻的射流橫截面形狀Fig.5 Slice view of cross-section shape of jet at t* =9.2

    圖6 對(duì)稱(chēng)切面內(nèi)的射流隨時(shí)間變化Fig.6 Temporal evolution of a jet in longitudinal slice through jet center

    3.2 展向擴(kuò)散特性

    除了具有彎曲特性,液體在橫向氣流中的噴注也隨時(shí)間展向擴(kuò)散,圖7給出了射流在4個(gè)不同時(shí)刻迎風(fēng)面的左視圖,虛線代表射流中心線,黑實(shí)線與α共同表征射流擴(kuò)散角。此處的擴(kuò)散角α定義為射流中心線與最左側(cè)未破碎極限之間的夾角。與t*=3時(shí)相比,t*=5對(duì)應(yīng)的展向面積明顯增寬,并在t*=7時(shí)進(jìn)一步擴(kuò)散,此時(shí)射流兩側(cè)邊緣由于強(qiáng)氣流剪切已經(jīng)出現(xiàn)破碎現(xiàn)象,氣流剪切強(qiáng)意味著液體的黏性系數(shù)低,因此相較其他位置更易破碎(該破碎過(guò)程將在3.3節(jié)詳細(xì)解釋?zhuān)T趖*=9.2時(shí),射流出現(xiàn)了高度破碎,形成若干小液滴,但擴(kuò)散角并未出現(xiàn)太大變化。

    圖7 射流的左視結(jié)構(gòu)圖Fig.7 Left view of jet structure

    圖8是射流擴(kuò)散角的時(shí)域變化,可以看到t*=3~7時(shí),擴(kuò)散角快速增加,t*=3時(shí)擴(kuò)散角僅為8.7°,到t*=5時(shí)已增長(zhǎng)3倍多,擴(kuò)散角為27.1°,在t*=7時(shí)更是達(dá)到34.4°,但在此之后擴(kuò)散角的增速明顯放緩,t*=9.2時(shí)仍維持在34.7°。從圖8可以推斷,液體在橫向氣流中噴注所具有的擴(kuò)散特性在噴入初期最為顯著,隨著時(shí)間增加擴(kuò)散將逐漸趨于

    圖8 射流擴(kuò)散角隨時(shí)間的變化Fig.8 Temporal evolution of spreading angle

    穩(wěn)定,這與文獻(xiàn)[10]中在相近動(dòng)量比q=6.6條件下發(fā)現(xiàn)的規(guī)律是類(lèi)似的。

    3.3 破碎特征

    對(duì)于液體射流,液絲與液滴的生成無(wú)可避免。文獻(xiàn)[4]將橫向氣流中的射流破碎分為兩類(lèi):圓柱破碎(Column Breakup)和表面破碎(Surface Breakup),本文將重點(diǎn)關(guān)注如圖9所示的邊緣破碎,屬于表面破碎類(lèi)別。

    圖9 液絲與液滴的形成過(guò)程Fig.9 Formation of ligaments and droplets

    在t*=7.6時(shí),受表面張力作用,細(xì)長(zhǎng)液絲的表面形成了頸部(圖9(a)藍(lán)色箭頭處)。而表面張力與曲率半徑成反比,即曲率半徑越小受力越大,因此頸部受到的表面張力是最大的。隨著時(shí)間推進(jìn),頸部將不斷收縮并在t*=7.9時(shí)斷裂,與此同時(shí),液絲上形成了兩個(gè)新的頸部(黑色箭頭處)。相似的現(xiàn)象在圖9(d)~圖9(f)中紅色與紫色箭頭處同樣可以看到。長(zhǎng)液絲在橫向氣流的作用下變得極不穩(wěn)定,頸部的形成以及最終的液絲破碎主要都是受瑞利不穩(wěn)定性影響。從圖9中的綠色箭頭還可以發(fā)現(xiàn),這些斷裂后的小體積液體在表面張力的作用下有不斷收縮成圓球形狀的趨勢(shì)。此外,圖9(b)~圖9(e)中黑色虛線還圈出了一個(gè)細(xì)致的液滴碰撞過(guò)程(最終同樣趨于圓球形)。

    3.4 渦特性

    渦量分布通常用來(lái)描述當(dāng)?shù)氐男魈匦?。圖10為本文計(jì)算得到的t*=9.2時(shí)的瞬時(shí)渦量場(chǎng),白色箭頭代表橫向氣流方向,渦量值ω根據(jù)射流直徑D以及流速u(mài)l進(jìn)行了無(wú)量綱處理,下標(biāo)x、y、z分別代表3個(gè)方向,即

    在射流的迎風(fēng)面,氣流中的高渦量主要集中在液體表面,并且旋流層的厚度沿著射流表面增加,近似于邊界層的發(fā)展。在x/D=2.5位置,受射流彎曲影響,高渦量區(qū)將脫離液體表面,并在液體與高渦量區(qū)之間形成了一個(gè)特殊的未擾動(dòng)層,同時(shí),高渦量區(qū)在下游再貼附液體表面時(shí)形成了一個(gè)軸向漩渦(如圖10(d)所示)。

    在射流的背風(fēng)面,存在一個(gè)非常大的高渦量尾跡區(qū),如圖10(e)所示,其渦量值范圍較圖10(a)大3倍。由于近噴嘴附近的射流近似為圓形,因此其流場(chǎng)與圓柱擾流類(lèi)似。在圓柱擾流研究中,文獻(xiàn)[34]通過(guò)對(duì)比不同來(lái)流雷諾數(shù)下單個(gè)圓柱擾流的尾跡渦發(fā)展及其阻力分布時(shí)就明確指出,當(dāng)來(lái)流Rec超過(guò)40時(shí),其擾流尾跡將出現(xiàn)不穩(wěn)定,而本文研究中橫向氣流的Rec遠(yuǎn)超該臨界值,達(dá)到了1 621,故出現(xiàn)極為復(fù)雜的尾跡區(qū)。

    在射流頭部附近區(qū)域M,橫向氣流流過(guò)液體時(shí)在下游形成了對(duì)漩渦(CVP),盡管與圖10(e)同為橫向擾流,但流動(dòng)形態(tài)顯然更加穩(wěn)定,這有以下3個(gè)主要原因:①圖10(b)中的液體具有x軸方向移動(dòng)速度,因此氣液兩相的相對(duì)速度小,Rec因而也較小;②液體的尺寸較噴嘴處的射流直徑小,這將進(jìn)一步降低Rec;③從前緣滯止點(diǎn)開(kāi)始沿液體表面發(fā)展的氣流邊界層在脫體前未完全發(fā)展,因此流動(dòng)也更穩(wěn)定。在圖10(c)中,由于氣液兩相的相對(duì)速度以及液體尺寸均較小,滿足渦的穩(wěn)定發(fā)展條件,因此還形成了類(lèi)似卡門(mén)渦街的流場(chǎng)結(jié)構(gòu)。

    圖10 t*=9.2時(shí)刻的瞬時(shí)渦量場(chǎng)Fig.10 Instantaneous vorticity field at t*=9.2

    3.5 非牛頓特征

    3.5.1 黏性系數(shù)

    圖11給出的是射流的黏性系數(shù)(無(wú)量綱,除以μl)沿y軸的平均變化規(guī)律??梢园l(fā)現(xiàn),在y/D<3時(shí),液體的黏性系數(shù)不斷下降,并在y/D=3~5之間維持在一個(gè)較低值,之后又開(kāi)始增加。通過(guò)分析流場(chǎng)可知,隨著y的增加,射流變得越來(lái)越?。▓D10(a)),液體被不斷剪切導(dǎo)致其黏性系數(shù)持續(xù)下降,在噴嘴出口處液體的黏性系數(shù)為0.95,到y(tǒng)/D=3時(shí)降至0.80左右。在y/D=3~5之間,由于液體的厚度僅微弱變化,因此黏性系數(shù)也基本維持在0.80附近。y/D>5后的黏性系數(shù)增加主要是受射流頭部較厚、內(nèi)部黏性系數(shù)較高的影響。

    圖11 t*=9.2時(shí)刻液體的黏性分布沿y軸的變化Fig.11 Spatial distribution of viscosity in liquid phase along ydirection at t*=9.2

    在圖2中可以觀察到,y/D=3~5之間的射流破碎最為顯著,其實(shí)這可以從圖11中該處液體的較低黏性系數(shù)(比μ0低20%)得到解釋?zhuān)驗(yàn)轲ば韵禂?shù)越低,液體越容易破碎。從定量角度看,低黏性系數(shù)代表著較大的當(dāng)?shù)乩字Z數(shù),因而此處的液體更為不穩(wěn)定。

    3.5.2 與牛頓液體射流對(duì)比

    圖12 t*=6.7時(shí)刻牛頓射流與非牛頓射流對(duì)比Fig.12 Comparison of jet structure from Newtonian and non-Newtonian at t* =6.7

    為了更清晰地展示非牛頓射流與牛頓射流的區(qū)別,本文也開(kāi)展了牛頓射流的數(shù)值模擬,并進(jìn)行對(duì)比,如圖12所示。在t*=6.7時(shí)刻,兩者的總體射流結(jié)構(gòu)是類(lèi)似的,但存在幾處主要區(qū)別:①射流邊緣的液體破碎,對(duì)非牛頓流體而言,由于邊緣處剪切率大,液體的黏性系數(shù)低,因此分裂成液絲與液滴,但在牛頓射流中,液體才剛剛出現(xiàn)破裂;②非牛頓射流的頭部形成一個(gè)液體薄層,而在牛頓射流中并不存在;③二者的展向擴(kuò)散特性不同。牛頓射流的最寬處為2.9D,而非牛頓射流的最寬處達(dá)到3.3D,兩者的擴(kuò)散角相差也超過(guò)了20%。因此即使液體只有弱的非牛頓特性,橫向氣流中的射流特征也會(huì)出現(xiàn)顯著差別。

    4 結(jié) 論

    1)受橫向氣流作用,射流噴出后存在彎曲與展向擴(kuò)散,并呈現(xiàn)bag類(lèi)型破碎。

    2)射流在彎曲過(guò)程中迎風(fēng)面的液體表面趨平,下游液體受剪切作用被逐漸拉成薄層。

    3)液體的展向擴(kuò)散在噴注初期最為顯著,后趨于穩(wěn)定,其擴(kuò)散角最終恒定在約35°。

    4)破碎形成的小液滴最高速與最低速相差4倍,圓形液滴直徑范圍為D/10~D/15。

    5)迎風(fēng)面的渦量沿液體表面發(fā)展,而背風(fēng)面氣流形成復(fù)雜尾跡區(qū),近似于圓柱繞流。

    6)近噴嘴處的射流液體,其高/低黏性系數(shù)相差超過(guò)20%,而在破碎區(qū)液體的黏性系數(shù)最低。

    7)與牛頓流體相比,剪切稀化非牛頓射流在頭部形成單獨(dú)的液體薄層且更易破碎,其展向擴(kuò)散角也更大。

    致 謝

    本文部分工作是在德國(guó)斯圖加特大學(xué)完成,因此特別感謝Bernhard Weigand教授和Moritz Ertl博士的幫助與討論,也要感謝斯圖加特高性能計(jì)算中心對(duì)本工作的大力支持。同時(shí),本文作者也要感謝廈門(mén)大學(xué)校長(zhǎng)基金對(duì)該課題的資助。

    [1] LEE K,AALBURG C,DIEZ F J,et al.Primary breakup of turbulent round liquid jets in uniform crossflows[J].AIAA Journal,2007,45(8):1907-1916.

    [2] SALLAM K A,AALBURG C,F(xiàn)AETH G M.Breakup of round nonturbulent liquid jets in gaseous crossflow[J].AIAA Journal,2004,42(12):2529-2540.

    [3] BIROUK M, NYANTEKYI-KWAKYE B,POPPLE-WELL N.Effect of nozzle geometry on breakup length and trajectory of liquid jet in subsonic crossflow[J].Atomization and Sprays,2011,21(10):847-865.

    [4] WU P,KIRKENDALL K A,F(xiàn)ULLER R P.Breakup process of liquid jets in subsonic crossflows[J].Journal of Propulsion and Power,1997,13(1):64-73.

    [5] STENZLER J N,LEE J G,SANTAVICCA D A.Penetration of liquid jets in cross-flow[J].Atomization and Sprays,2006,16(8):887-906.

    [6] GUTMARK E J,IBRAHIM I M,MURUGAPPAN S.Circular and noncircular subsonic jets in cross flow[J].Physics of Fluids,2008,20(7):075110.

    [7] SHAPIRO S R,KING J M,CLOSKEY R T M.Optimization of controlled jets in crossflow[J].AIAA Journal,2006,44(6):1292-1298.

    [8] MEGERIAN S,DAVITIAN J,ALVES L S B.Transverse-jet shear-layer instabilities.Part 1.Experimental studies[J].Journal of Fluid Mechanics,2007,593:93-129.

    [9] COLETTI F,ELKINS C J,EATON J K.An inclined jet in crossflow under the effect of streamwise pressure gradients[J].Experiments in Fluids,2013,54:1589.

    [10] HERRMANN M.Detailed numerical simulations of the primary atomization of a turbulent liquid jet in crossflow[J].Journal of Engineering for Gas Turbine and Power,2010,132(6):061506.

    [11] HERRMANN M.The influence of density ratio on the primary atomization of a turbulent liquid jet in crossflow[J].Proceedings of Combustion Institute,2011,33(2):2079-2088.

    [12] CAVAR D,MEYER K E.LES of turbulent jet in cross flow:Part 2.POD analysis and identification of coherent structures[J].International Journal of Heat and Fluid Flow,2012,36:35-46.

    [13] GALEAZZO F C C,DONNERT G,CARDENAS C,et al.Computational modeling of turbulent mixing in a jet in crossflow[J].International Journal of Heat and Fluid Flow,2013,41:55-65.

    [14] SAU R,MAHESH K.Optimization of pulsed jets in crossflow[J].Journal of Fluid Mechanics,2010,653:365-390.

    [15] PAI M G,DESJARDINS O,PITSCH H.Detailed simulations of primary breakup of turbulent liquid jets in crossflow[R].Stanford,CA:Center for Turbulence Research,2008:451-466.

    [16] MULDOON F,ACHARYA S.Direct numerical simulation of pulsed jets in crossflow[J].Computers & Fluids,2010,39(10):1745-1773.

    [17] MARGASON R J.Fifty years of jet in crossflow research:AGARD CP[R].Paris:AGARD,1993,534:1-41.

    [18] AALBURG C,VAN LEER B,F(xiàn)AETH G M,et al.Properties of nonturbulent round liquid jets in uniform gaseous cross flows[J].Atomization and Sprays,2005,15(3):271-294.

    [19] MAHESH K.The interaction of jets with crossflow[J].Annual Review of Fluid Mechanics,2013,45:379-407.

    [20] WONG D C Y,SIMMONS M J H,DECENT S P,et al.Break-up dynamics and drop size distributions created from spiraling liquid jets[J].International Journal of Multiphase Flow,2004,30:499-520.

    [21] CLASEN C,EGGERS J,F(xiàn)ONTELOS M A,et al.The beads-on-string structure of viscoelastic threads[J].Journal of Fluid Mechanics,2006,556:283-308.

    [22] YARIN A L.Free liquid jets and films:Hydrodynamics and rheology[M].New York:Wiley,1993.

    [23] HIRT C W,NICHOLS B D.Volume of fluid (VOF)method for the dynamics of free boundaries[J].Journal of Computational Physics,1981,39:201-225.

    [24] RIDER W J,KOTHE D B.Reconstructing volume tracking[J].Journal of Computational Physics,1998,141:112-152.

    [25] KORNEV N,HASSEL E.Synthesis of homogeneous anisotropic divergence free turbulent fields with prescribed second-order statistics by vortex dipoles[J].Physics of Fluids,2007,19(6):068101.

    [26] RIDER W J,KOTHE D B.Reconstructing volume tracking[J].Journal of Computational Physics,1998,141:112-152.

    [27] GOMAA H,KUMAR S,HUBER C,et al.Numerical comparison of 3Djet breakup using a compression scheme and an interface reconstruction based VOF-code[C]/24th ILASS-Europe,2011.

    [28] MOTZIGEMBA M,ROTH N,BOTHE D,et al.The effect of non-Newtonian flow behavior on binary droplet collisions: VOF-simulation and experimental analysis[C]/Proceedings of ILASS-Europe,2002.

    [29] FOCKE C,BOTHE D.Computational analysis of binary collisions of shear thinning droplets[J].Journal of Non-Newtonian Fluid Mechanics,2011,166:799-810.

    [30] ZHU C,ERTL M,WEIGAND B.Numerical investigation on the primary breakup of an inelastic non-Newtonian liquid jet with inflow turbulence[J].Physics of Fluids,2013,25(8):083102.

    [31] SCHROEDER J,LEDERER M L,GAUKEL V,et al.Effect of atomizer geometry and rheological properties on effervescent atomization of aqueous polyvinylphrrolidone solution[C]/24th ILASS-Europe,2011.

    [32] YOU Y,LUEDEKE H,HANNEMANN K.On the flow physics of a low momentum flux ratio jet in a supersonicturbulent crossflow[J].Europhysics Letters,2012,97(2):24001.

    [33] BATCHELOR G K.The theory of homogeneous turbulence[M].Cambridge:Cambridge University Press,1953.

    [34] MUNSON B R,YOUNG D F,OKIISHI T H.Fundamentals of fluid mechanics[M].New York:John Wiley&Sons,2006.

    Direct numerical simulation of a non-Newtonian liquid jet in crossflow

    ZHU Chengxiang*,YOU Yancheng
    School of Aerospace Engineering,Xiamen University,Xiamen 361005,China

    A direct numerical simulation study of a non-Newtonian liquid jet in crossflow is carried out with a moderate momentum flux ratio 6.The emphasis of this paper mainly focuses on the flow structure of the jet,including surface behaviors,bending phenomena,spreading features and non-Newtonian characteristics.Deep into the near-field region,it can be observed that the trajectory of the jet oscillates with time and has a tendency to move backward in the reverse direction of the crossflow.The spreading angle increases only at the start of the injection but keeps nearly constant afterwards at 35°.Further insight into the flow physics is obtained by visualizing the primary breakup of the jet,especially the formation of ligaments and droplets,as well as satellite droplets.In the near-field region,the flow feature is similar to that of a circular cylinder,while showing complex turbulent behavior in other regions.The specific non-Newtonian characteristics of the fluid are observed by analyzing the shear thinning viscosity,which varies over 20%spatially inside the liquid.Compared to Newtonian fluids,the current shear thinning non-Newtonian liquid jet shows a stronger breakup feature.

    liquid jet;crossflow;non-Newtonian fluid;breakup;direct numerical simulation

    2015-09-15;Revised:2015-11-03;Accepted:2016-01-12;Published online:2016-01-31 12:57

    URL:www.cnki.net/kcms/detail/11.1929.V.20160131.1257.012.html

    s:National Natural Science Foundation of China(91441128;51276151)

    V231.2

    A

    1000-6893(2016)09-2659-10

    10.7527/S1000-6893.2016.0004

    2015-09-15;退修日期:2015-11-03;錄用日期:2016-01-12;網(wǎng)絡(luò)出版時(shí)間:2016-01-31 12:57

    www.cnki.net/kcms/detail/11.1929.V.20160131.1257.012.html

    國(guó)家自然科學(xué)基金(91441128;51276151)

    *通訊作者.Tel.:0592-2186849 E-mail:chengxiang.zhu@xmu.edu.cn

    朱呈祥,尤延鋮.橫向氣流中非牛頓液體射流直接數(shù)值模擬[J].航空學(xué)報(bào),2016,37(9):26592-668.ZHU C X,YOU Y C.Direct numerical simulation of a non-Newtonian liquid jet in crossflow[J].Acta Aeronautica et Astronautica Sinica,2016,37(9):26592-688.

    朱呈祥 男,博士,講師,碩士生導(dǎo)師。主要研究方向:氣液兩相流,非牛頓流體力學(xué),氣體動(dòng)力學(xué)。Tel:0592-2186849

    E-mail:chengxiang.zhu@xmu.edu.cn尤延鋮 男,博士,教授,博士生導(dǎo)師。主要研究方向:高超聲速流動(dòng),計(jì)算流體力學(xué)。Tel:0592-2186849

    E-mail:yancheng.you@xmu.edu.cn

    *Corresponding author.Tel.:0592-2186849 E-mail:chengxiang.zhu@xmu.edu.cn

    猜你喜歡
    牛頓流體氣液黏性
    微重力下兩相控溫型儲(chǔ)液器內(nèi)氣液界面仿真分析
    非牛頓流體
    氣液分離罐液位計(jì)接管泄漏分析
    什么是非牛頓流體
    少兒科技(2019年3期)2019-09-10 07:22:44
    富硒產(chǎn)業(yè)需要強(qiáng)化“黏性”——安康能否玩轉(zhuǎn)“硒+”
    如何運(yùn)用播音主持技巧增強(qiáng)受眾黏性
    區(qū)別牛頓流體和非牛頓流體
    玩油灰黏性物成網(wǎng)紅
    首款XGEL非牛頓流體“高樂(lè)高”系列水溶肥問(wèn)世
    CO2 驅(qū)低液量高氣液比井下氣錨模擬與優(yōu)化
    老熟女久久久| 婷婷色综合www| 狂野欧美激情性xxxx在线观看| 亚洲精品乱码久久久久久按摩| 国产熟女欧美一区二区| 亚洲国产日韩一区二区| 欧美97在线视频| 国产成人精品一,二区| 成人综合一区亚洲| 久久亚洲国产成人精品v| 91精品一卡2卡3卡4卡| 国产精品99久久久久久久久| 日韩欧美一区视频在线观看 | 91午夜精品亚洲一区二区三区| 22中文网久久字幕| 亚洲av中文av极速乱| 亚洲天堂av无毛| 亚洲精品日韩av片在线观看| 亚洲综合色惰| 国产男女内射视频| 久久ye,这里只有精品| 视频中文字幕在线观看| 91在线精品国自产拍蜜月| 啦啦啦中文免费视频观看日本| 日本爱情动作片www.在线观看| 日韩欧美精品免费久久| 久久国内精品自在自线图片| 人妻少妇偷人精品九色| 最近2019中文字幕mv第一页| 啦啦啦啦在线视频资源| 亚洲成人一二三区av| 热re99久久精品国产66热6| 天堂俺去俺来也www色官网| 韩国高清视频一区二区三区| 妹子高潮喷水视频| 亚洲av综合色区一区| 色视频www国产| 免费看日本二区| 18禁动态无遮挡网站| 成年av动漫网址| 国产日韩欧美亚洲二区| 观看免费一级毛片| 观看免费一级毛片| 91精品一卡2卡3卡4卡| 波野结衣二区三区在线| 两个人的视频大全免费| 深夜a级毛片| 亚洲伊人久久精品综合| 尾随美女入室| 80岁老熟妇乱子伦牲交| 岛国毛片在线播放| 成人免费观看视频高清| 两个人的视频大全免费| 自拍偷自拍亚洲精品老妇| 视频区图区小说| 精品一区在线观看国产| 一二三四中文在线观看免费高清| 免费大片18禁| av国产久精品久网站免费入址| 在线 av 中文字幕| 一级av片app| 精品亚洲成a人片在线观看| 老女人水多毛片| 免费黄频网站在线观看国产| 国产精品国产av在线观看| 人妻 亚洲 视频| 五月开心婷婷网| 国产黄频视频在线观看| 高清av免费在线| 国产高清有码在线观看视频| 亚洲精品国产av成人精品| 啦啦啦啦在线视频资源| 精品人妻一区二区三区麻豆| 亚洲久久久国产精品| 只有这里有精品99| a级毛片免费高清观看在线播放| 欧美日韩国产mv在线观看视频| 亚洲人与动物交配视频| 男女无遮挡免费网站观看| 成人免费观看视频高清| 亚洲精品久久久久久婷婷小说| 亚洲自偷自拍三级| 黑人猛操日本美女一级片| 热re99久久国产66热| 99热全是精品| 精品人妻一区二区三区麻豆| 国产一区二区三区综合在线观看 | 国产黄色免费在线视频| 日韩,欧美,国产一区二区三区| 亚洲内射少妇av| 97在线视频观看| 日日爽夜夜爽网站| 97精品久久久久久久久久精品| 国产精品蜜桃在线观看| 免费黄频网站在线观看国产| 日日摸夜夜添夜夜添av毛片| 大话2 男鬼变身卡| 久久国产乱子免费精品| 色哟哟·www| 午夜福利视频精品| 国产片特级美女逼逼视频| 亚洲av综合色区一区| 成人综合一区亚洲| 日本黄色日本黄色录像| av.在线天堂| 人妻夜夜爽99麻豆av| 日本午夜av视频| 国产亚洲欧美精品永久| 亚洲情色 制服丝袜| 亚洲精品日本国产第一区| 亚洲图色成人| 亚洲三级黄色毛片| 国产男女超爽视频在线观看| 久久99一区二区三区| 国产黄片美女视频| 国产精品一区二区三区四区免费观看| av网站免费在线观看视频| 亚洲精品,欧美精品| 日韩欧美 国产精品| 国产 一区精品| 婷婷色麻豆天堂久久| 日韩伦理黄色片| 噜噜噜噜噜久久久久久91| 能在线免费看毛片的网站| 妹子高潮喷水视频| 午夜影院在线不卡| 久久精品国产亚洲av涩爱| 纵有疾风起免费观看全集完整版| 日本色播在线视频| 色哟哟·www| 99久久综合免费| 免费看日本二区| 夜夜爽夜夜爽视频| 亚洲中文av在线| 日产精品乱码卡一卡2卡三| 九九久久精品国产亚洲av麻豆| 国产精品欧美亚洲77777| 一级毛片久久久久久久久女| 免费观看的影片在线观看| 人妻系列 视频| 日日啪夜夜爽| 精品久久久噜噜| 国产探花极品一区二区| 午夜福利在线观看免费完整高清在| 26uuu在线亚洲综合色| 91久久精品国产一区二区三区| 亚洲熟女精品中文字幕| 伦理电影免费视频| 国产精品无大码| 国产亚洲欧美精品永久| 日韩av在线免费看完整版不卡| 色94色欧美一区二区| 亚洲欧美精品自产自拍| 欧美精品一区二区免费开放| 亚洲国产欧美在线一区| 日本爱情动作片www.在线观看| 色婷婷久久久亚洲欧美| 亚洲图色成人| 看免费成人av毛片| 免费观看av网站的网址| 国产 一区精品| 又粗又硬又长又爽又黄的视频| 搡老乐熟女国产| 男人和女人高潮做爰伦理| 晚上一个人看的免费电影| 熟女av电影| 男人舔奶头视频| 五月玫瑰六月丁香| 啦啦啦视频在线资源免费观看| 99热全是精品| 久久精品熟女亚洲av麻豆精品| 人人妻人人澡人人爽人人夜夜| 91午夜精品亚洲一区二区三区| 亚洲国产精品专区欧美| 日本爱情动作片www.在线观看| 精品人妻熟女av久视频| 丝袜喷水一区| 伦理电影大哥的女人| 亚洲精品视频女| 国产男人的电影天堂91| 国产成人免费无遮挡视频| 水蜜桃什么品种好| 夜夜爽夜夜爽视频| 国产男女超爽视频在线观看| 欧美精品一区二区免费开放| 一本—道久久a久久精品蜜桃钙片| 国产av国产精品国产| 免费黄色在线免费观看| 精品亚洲成国产av| 国产在线免费精品| 男人舔奶头视频| 国产精品免费大片| 亚洲va在线va天堂va国产| 国产在线一区二区三区精| 亚洲在久久综合| 亚洲精品乱码久久久久久按摩| 日韩av免费高清视频| 热re99久久国产66热| 亚洲精品中文字幕在线视频 | 99久久精品一区二区三区| 91久久精品电影网| 18禁裸乳无遮挡动漫免费视频| 国产精品久久久久久av不卡| 国产亚洲一区二区精品| 国产淫语在线视频| 男女免费视频国产| 国产成人精品婷婷| 日韩亚洲欧美综合| 在线观看一区二区三区激情| 免费av不卡在线播放| 亚洲怡红院男人天堂| 最后的刺客免费高清国语| 亚洲精品一区蜜桃| 99国产精品免费福利视频| 国产白丝娇喘喷水9色精品| 国产高清三级在线| 91久久精品国产一区二区成人| 在线天堂最新版资源| 婷婷色av中文字幕| 国产精品人妻久久久影院| 久久久午夜欧美精品| 一级毛片黄色毛片免费观看视频| 日本黄色日本黄色录像| 亚洲欧美精品自产自拍| 自拍偷自拍亚洲精品老妇| 欧美国产精品一级二级三级 | 久久人妻熟女aⅴ| 国产一区亚洲一区在线观看| 大香蕉97超碰在线| 三上悠亚av全集在线观看 | 久久99热这里只频精品6学生| 嘟嘟电影网在线观看| 欧美高清成人免费视频www| 亚洲av欧美aⅴ国产| 日本与韩国留学比较| 欧美国产精品一级二级三级 | 美女大奶头黄色视频| 少妇的逼好多水| 午夜免费鲁丝| 亚洲精品中文字幕在线视频 | 日日撸夜夜添| 最近手机中文字幕大全| 亚洲国产精品国产精品| av福利片在线观看| 18禁裸乳无遮挡动漫免费视频| 中文字幕精品免费在线观看视频 | 欧美日韩视频高清一区二区三区二| 成人二区视频| 国产毛片在线视频| 在线观看免费日韩欧美大片 | av在线播放精品| 欧美成人精品欧美一级黄| 国产一区二区三区综合在线观看 | 欧美精品高潮呻吟av久久| 三级经典国产精品| 3wmmmm亚洲av在线观看| 精品少妇久久久久久888优播| 韩国高清视频一区二区三区| 国产欧美日韩精品一区二区| 国产亚洲一区二区精品| 国产成人aa在线观看| 精品一品国产午夜福利视频| 永久网站在线| 国产成人精品婷婷| 亚洲国产色片| 精品一品国产午夜福利视频| 亚洲av欧美aⅴ国产| av一本久久久久| 麻豆乱淫一区二区| 97超碰精品成人国产| 黄片无遮挡物在线观看| 免费观看性生交大片5| 久久久久久久大尺度免费视频| 91成人精品电影| 插逼视频在线观看| 一本色道久久久久久精品综合| 我要看黄色一级片免费的| 免费人成在线观看视频色| 国产免费一级a男人的天堂| 中文精品一卡2卡3卡4更新| 下体分泌物呈黄色| 久热这里只有精品99| 国产成人一区二区在线| 日本色播在线视频| 精品卡一卡二卡四卡免费| 中文字幕久久专区| 人妻夜夜爽99麻豆av| 伦理电影免费视频| 国产亚洲一区二区精品| 不卡视频在线观看欧美| 国内揄拍国产精品人妻在线| av国产久精品久网站免费入址| 十八禁高潮呻吟视频 | 久久午夜福利片| 亚洲av成人精品一区久久| 边亲边吃奶的免费视频| 国产精品欧美亚洲77777| 日韩一区二区三区影片| 欧美激情国产日韩精品一区| 国产熟女午夜一区二区三区 | 大香蕉97超碰在线| 99久久精品一区二区三区| 蜜臀久久99精品久久宅男| 日韩成人av中文字幕在线观看| 香蕉精品网在线| 国产精品.久久久| 国产成人91sexporn| 春色校园在线视频观看| 亚洲第一av免费看| 亚洲精品日韩在线中文字幕| 丰满迷人的少妇在线观看| 一级av片app| 亚洲成人av在线免费| 女人精品久久久久毛片| 国产精品熟女久久久久浪| 免费看不卡的av| 少妇被粗大猛烈的视频| 午夜老司机福利剧场| av线在线观看网站| 永久网站在线| 丰满乱子伦码专区| 最近的中文字幕免费完整| 久久久久久久久久成人| 成人国产麻豆网| 日韩免费高清中文字幕av| 久久国产精品男人的天堂亚洲 | 成人二区视频| 亚洲经典国产精华液单| 亚洲伊人久久精品综合| 大香蕉97超碰在线| 赤兔流量卡办理| 精品久久久久久久久亚洲| 国语对白做爰xxxⅹ性视频网站| 中文字幕久久专区| 黄片无遮挡物在线观看| 日韩成人av中文字幕在线观看| 又爽又黄a免费视频| 欧美三级亚洲精品| 黑人猛操日本美女一级片| 婷婷色av中文字幕| 大陆偷拍与自拍| 日韩av免费高清视频| 成人二区视频| 波野结衣二区三区在线| 日韩三级伦理在线观看| 国产乱来视频区| 不卡视频在线观看欧美| 亚洲国产精品专区欧美| 丝袜脚勾引网站| 欧美日本中文国产一区发布| 亚洲怡红院男人天堂| 国语对白做爰xxxⅹ性视频网站| 亚洲人成网站在线观看播放| 看非洲黑人一级黄片| 一本久久精品| 99久久精品一区二区三区| 成人黄色视频免费在线看| 久久午夜福利片| 欧美老熟妇乱子伦牲交| 欧美日韩在线观看h| 成人无遮挡网站| 国产伦精品一区二区三区四那| 精品久久久久久久久亚洲| 热99国产精品久久久久久7| 高清视频免费观看一区二区| 婷婷色综合大香蕉| av不卡在线播放| 尾随美女入室| 国产一区二区在线观看av| 春色校园在线视频观看| 在线播放无遮挡| 日本wwww免费看| 国产老妇伦熟女老妇高清| 日韩免费高清中文字幕av| 午夜免费男女啪啪视频观看| 午夜av观看不卡| 最后的刺客免费高清国语| 欧美变态另类bdsm刘玥| 亚洲天堂av无毛| 婷婷色综合www| 免费看不卡的av| 精品久久久久久电影网| 久久精品夜色国产| 伦精品一区二区三区| 亚洲一级一片aⅴ在线观看| 亚洲人与动物交配视频| 久久久国产精品麻豆| 国产精品免费大片| 99热这里只有精品一区| 男女免费视频国产| 久久99热这里只频精品6学生| 久久久久久久亚洲中文字幕| 亚洲av欧美aⅴ国产| 国产精品久久久久久av不卡| 国产精品麻豆人妻色哟哟久久| 精华霜和精华液先用哪个| 99久国产av精品国产电影| 99re6热这里在线精品视频| 精品人妻熟女av久视频| 成人午夜精彩视频在线观看| 日韩成人伦理影院| 国产成人精品无人区| 最近中文字幕高清免费大全6| 精品亚洲成a人片在线观看| 亚洲av中文av极速乱| 亚洲综合精品二区| a级片在线免费高清观看视频| 亚洲综合色惰| 在线免费观看不下载黄p国产| 国产成人freesex在线| 在线观看三级黄色| 女的被弄到高潮叫床怎么办| 午夜福利网站1000一区二区三区| av天堂久久9| 国产精品伦人一区二区| 国产免费一级a男人的天堂| 日日摸夜夜添夜夜爱| 日韩精品免费视频一区二区三区 | 五月天丁香电影| 亚洲国产日韩一区二区| 黑人猛操日本美女一级片| 亚洲av在线观看美女高潮| 最新的欧美精品一区二区| 一本久久精品| 国产欧美亚洲国产| 97在线视频观看| 久久精品国产a三级三级三级| 国产精品偷伦视频观看了| 中文乱码字字幕精品一区二区三区| 欧美激情国产日韩精品一区| 亚洲国产欧美在线一区| 男人狂女人下面高潮的视频| 男女免费视频国产| 少妇猛男粗大的猛烈进出视频| 成年女人在线观看亚洲视频| 26uuu在线亚洲综合色| 麻豆乱淫一区二区| 人妻夜夜爽99麻豆av| 99久久综合免费| 久热这里只有精品99| 这个男人来自地球电影免费观看 | 精品久久久久久久久亚洲| 国产美女午夜福利| 国产精品一区二区三区四区免费观看| 久久久精品94久久精品| 观看免费一级毛片| 国产欧美日韩综合在线一区二区 | 亚洲精品一二三| 国产精品三级大全| 日本91视频免费播放| 69精品国产乱码久久久| 51国产日韩欧美| 国产成人精品久久久久久| 日日摸夜夜添夜夜爱| 91精品一卡2卡3卡4卡| 亚洲综合色惰| 国产精品伦人一区二区| 91午夜精品亚洲一区二区三区| 亚洲婷婷狠狠爱综合网| 国产毛片在线视频| 一本—道久久a久久精品蜜桃钙片| 热re99久久精品国产66热6| 22中文网久久字幕| 亚洲精品国产av成人精品| 伦理电影大哥的女人| av天堂久久9| 国产淫语在线视频| 国产精品人妻久久久影院| 久久人人爽av亚洲精品天堂| 亚洲一区二区三区欧美精品| 一本大道久久a久久精品| 女人久久www免费人成看片| 久久午夜综合久久蜜桃| 老女人水多毛片| 亚洲欧洲国产日韩| 另类精品久久| 在线看a的网站| 大片免费播放器 马上看| 亚洲性久久影院| 多毛熟女@视频| 国产免费又黄又爽又色| 日韩精品有码人妻一区| 啦啦啦啦在线视频资源| a级毛片在线看网站| 大香蕉97超碰在线| 三级经典国产精品| 精品少妇久久久久久888优播| 午夜视频国产福利| 免费不卡的大黄色大毛片视频在线观看| 国产深夜福利视频在线观看| 日本与韩国留学比较| 亚洲欧美一区二区三区国产| 日本黄大片高清| 美女内射精品一级片tv| 伊人久久精品亚洲午夜| 伦理电影免费视频| 亚洲第一av免费看| 色婷婷av一区二区三区视频| 丝袜在线中文字幕| 日日啪夜夜撸| 欧美日韩av久久| 午夜福利网站1000一区二区三区| 国语对白做爰xxxⅹ性视频网站| 伦理电影大哥的女人| 国产高清有码在线观看视频| 日韩欧美一区视频在线观看 | 亚洲精品456在线播放app| 一区二区三区精品91| 一级黄片播放器| 亚洲av日韩在线播放| 大香蕉97超碰在线| 黑丝袜美女国产一区| 少妇猛男粗大的猛烈进出视频| 不卡视频在线观看欧美| 综合色丁香网| 亚洲国产日韩一区二区| 日本爱情动作片www.在线观看| 高清视频免费观看一区二区| 成人漫画全彩无遮挡| 午夜视频国产福利| av一本久久久久| 欧美高清成人免费视频www| 男女国产视频网站| 国产精品国产av在线观看| av一本久久久久| 午夜视频国产福利| 国产精品一区二区在线不卡| 国产精品不卡视频一区二区| 我要看日韩黄色一级片| 在线播放无遮挡| 丰满迷人的少妇在线观看| 亚洲无线观看免费| 成人漫画全彩无遮挡| 久久亚洲国产成人精品v| 我要看日韩黄色一级片| 久久av网站| 综合色丁香网| 日韩,欧美,国产一区二区三区| 久久97久久精品| 青春草亚洲视频在线观看| 国产精品不卡视频一区二区| 青青草视频在线视频观看| 国产成人精品婷婷| videos熟女内射| 亚洲自偷自拍三级| 国产精品国产av在线观看| 各种免费的搞黄视频| 成人午夜精彩视频在线观看| 免费大片18禁| 欧美最新免费一区二区三区| 毛片一级片免费看久久久久| 黄色日韩在线| 亚洲国产精品成人久久小说| 亚洲精品日韩在线中文字幕| 丰满饥渴人妻一区二区三| av免费观看日本| av黄色大香蕉| 欧美精品国产亚洲| 久久久欧美国产精品| 人人妻人人添人人爽欧美一区卜| 亚洲中文av在线| 日本欧美国产在线视频| 亚洲av中文av极速乱| 国产91av在线免费观看| 最近最新中文字幕免费大全7| 国产伦精品一区二区三区四那| 亚洲情色 制服丝袜| 看十八女毛片水多多多| 99久久精品国产国产毛片| 欧美xxⅹ黑人| 久久精品久久久久久久性| 国产一区二区在线观看日韩| 欧美精品高潮呻吟av久久| 亚洲精品国产成人久久av| 亚洲国产精品999| 在线观看免费高清a一片| 亚洲第一区二区三区不卡| 亚洲精品第二区| av福利片在线| 久久国产精品男人的天堂亚洲 | 中文在线观看免费www的网站| 夫妻午夜视频| 国产在线视频一区二区| 国产av国产精品国产| 99视频精品全部免费 在线| 男人和女人高潮做爰伦理| 这个男人来自地球电影免费观看 | 水蜜桃什么品种好| 一区二区三区精品91| 自线自在国产av| 在线观看三级黄色| 另类亚洲欧美激情| 日韩成人av中文字幕在线观看| 欧美日本中文国产一区发布| 成人国产麻豆网| 精品卡一卡二卡四卡免费| 国产欧美亚洲国产| 大片免费播放器 马上看| 免费av不卡在线播放| 夜夜骑夜夜射夜夜干| 高清黄色对白视频在线免费看 | 免费大片黄手机在线观看| 毛片一级片免费看久久久久| 欧美日韩综合久久久久久| 欧美 亚洲 国产 日韩一| 制服丝袜香蕉在线| 国产黄色视频一区二区在线观看| 日本与韩国留学比较| 亚洲精品国产av蜜桃| 18禁裸乳无遮挡动漫免费视频| 久久这里有精品视频免费| 久久影院123| 精品酒店卫生间| 精品久久久精品久久久| 午夜福利在线观看免费完整高清在|