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

    基于VOF方法的大水滴袋狀破碎的仿真研究

    2018-10-08 12:06:40張文英常士楠雷夢(mèng)龍
    關(guān)鍵詞:水環(huán)無量不穩(wěn)定性

    張文英, 常士楠, 蔣 斌, 雷夢(mèng)龍

    (北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 100191)

    0 引 言

    凍雨的形成是由于大量存在于高空云頂?shù)谋ЯW踊蛞簯B(tài)水滴在暖層中凝結(jié)下落,穿過較薄的冷層,由于高層大氣缺少凝結(jié)核而未能發(fā)生相變,仍以液態(tài)水的形式存在。以凍雨或凍毛細(xì)雨形式存在高空中的大水滴(d>50 μm),我們稱之為過冷大水滴(Supercooled Large Droplets,SLD)。

    由于大水滴在空中較之地面的數(shù)密度和雨強(qiáng)貢獻(xiàn)率較大,它對(duì)飛機(jī)表面的撞擊結(jié)冰比常規(guī)尺寸的水滴危害更大[1]。研究常規(guī)尺寸的水滴可以忽略在運(yùn)動(dòng)中的形變和破碎,假設(shè)其始終為球形,且不會(huì)發(fā)生飛濺等現(xiàn)象。而大水滴由于直徑較大,水滴的Oh數(shù)(Ohnesorge數(shù))較小(通常遠(yuǎn)小于0.1),意味著水滴表面張力和內(nèi)部粘性力的量級(jí)小于所受的氣動(dòng)力,這表征著在水滴撞擊機(jī)翼表面前就會(huì)發(fā)生水滴的變形和破碎,并且在撞擊飛機(jī)迎風(fēng)表面后會(huì)發(fā)生飛濺,生成尺寸更小的子水滴,或?qū)⒃谖丛O(shè)防護(hù)的較大區(qū)域上發(fā)生二次撞擊,出現(xiàn)不可控的結(jié)冰,危害飛機(jī)的安全飛行。由此可見過冷大水滴導(dǎo)致的飛機(jī)結(jié)冰現(xiàn)象較之常規(guī)水滴更為復(fù)雜,引起了各國(guó)研究人員的重視[2-14]。

    圍繞運(yùn)動(dòng)液滴的破碎與變形,國(guó)內(nèi)外的研究人員傾注了大量的心血,并取得了一些重要的研究成果。國(guó)外對(duì)于大水滴的實(shí)驗(yàn)研究起步較早,Hinze[2]早在上世紀(jì)中期就開展了在高液-氣密度比和高Re數(shù)條件下液滴破碎的研究,提出了在一定Oh數(shù)條件下,破碎模式隨著We數(shù)改變而變化。Krzeczkowski[3]在Hinze的基礎(chǔ)上進(jìn)行了繼續(xù)研究,確定了破碎模式和We數(shù)的對(duì)應(yīng)關(guān)系,并且分析了典型的袋狀和剪切破碎的形態(tài)特征。Chou和Faeth[4]做了在高液-氣密度比和低Oh數(shù)(Oh?0.1)條件下的液滴破碎的瞬態(tài)特性研究,提出了袋裝破碎時(shí)間與水滴特征時(shí)間相關(guān)的結(jié)論。Rimbert和Castanet[5]認(rèn)為可以用R-T不穩(wěn)定性來解釋袋狀破碎的機(jī)理,或者是由于液滴內(nèi)部的流動(dòng)導(dǎo)致。Guildenbecher[6]對(duì)液滴破碎的實(shí)驗(yàn)文獻(xiàn)做了綜述,并且提出液滴袋狀破碎的機(jī)理是其內(nèi)部液體由兩極向赤道位置的流動(dòng)。

    盡管在實(shí)驗(yàn)文章中,研究人員們已經(jīng)給出相對(duì)可靠的破碎模式和We數(shù)對(duì)應(yīng)的范圍,但是由于測(cè)量技術(shù)的限制,難以得到流場(chǎng)中的壓力分布和液滴內(nèi)部的流動(dòng),因此水滴變形和破碎的機(jī)理尚未明確。伴隨著計(jì)算機(jī)技術(shù)和圖形學(xué)的發(fā)展,研究人員針對(duì)液滴的變形和破碎在數(shù)值仿真方面也取得了一定的成就。Zaleski[7]進(jìn)行了大液滴在Re=1000條件下的數(shù)值仿真,追蹤了相界面的細(xì)節(jié)變化和液滴變形處氣流的渦流結(jié)構(gòu)。Desjardins[8]、Lebas[9]、Sander和Weigand[10]各自做了液滴在不同We數(shù)條件下變形破碎的仿真研究,但是由于研究的重點(diǎn)不同,他們更關(guān)注的是直接數(shù)值仿真或者是大渦模擬的算法應(yīng)用,在防冰以及大水滴的破碎等方面應(yīng)用略顯不足。

    國(guó)內(nèi)的相關(guān)研究主要以仿真為主。西北工業(yè)大學(xué)的胡劍平[11]以及上海交通大學(xué)的羅輝[12]等人側(cè)重于應(yīng)用現(xiàn)有的水滴動(dòng)力學(xué)模型去仿真研究SLD對(duì)飛機(jī)部件結(jié)冰的影響。他們的研究結(jié)果證明了,考慮大水滴動(dòng)力學(xué)特性的仿真結(jié)果相較于只考慮常規(guī)尺寸水滴更貼近實(shí)驗(yàn)結(jié)果,充分證明了研究大水滴運(yùn)動(dòng)變形的必要性。然而僅僅從宏觀角度對(duì)水滴運(yùn)動(dòng)行為學(xué)進(jìn)行研究,不能從機(jī)理上解釋水滴變形破碎和飛濺的現(xiàn)象,因此,針對(duì)單個(gè)水滴的運(yùn)動(dòng)和變形破碎進(jìn)行研究是很有意義的。在此方面,北京航空航天大學(xué)的常士楠課題組[13-14]做了關(guān)于水滴袋狀、復(fù)合破碎過程中變形特性研究,對(duì)相關(guān)的理論模型進(jìn)行了改進(jìn),實(shí)現(xiàn)了水滴變形量和破碎時(shí)間的預(yù)測(cè)。

    本文應(yīng)用VOF方法對(duì)不同韋伯?dāng)?shù)下的水滴進(jìn)行數(shù)值模擬,關(guān)注兩相界面變化,得到與實(shí)驗(yàn)一致的變化規(guī)律,應(yīng)用R-T不穩(wěn)定性解釋袋狀破碎現(xiàn)象;統(tǒng)計(jì)了水滴的變形、破碎時(shí)間和運(yùn)動(dòng)速度等瞬態(tài)特性,與前人的實(shí)驗(yàn)數(shù)據(jù)和理論模型加以對(duì)比,揭示其特性與原理。

    1 數(shù)值方法與物理模型

    在兩相流的數(shù)值仿真方向,主要有兩大類方法:界面追蹤法(Interface tracking Method)和界面捕獲法(Interface capturing Method),它們分別是基于拉格朗日和歐拉方法。前者將兩相流其中一相視為大量微小粒子的集合,然而它的局限性在于不能自動(dòng)修改拓?fù)浣Y(jié)構(gòu),并且很難實(shí)現(xiàn)并行運(yùn)算。而基于歐拉方法的界面捕獲法,包括本文采用的VOF(Volume of Fluid)方法和LSM(Level Set Method)方法。VOF方法追蹤每個(gè)網(wǎng)格內(nèi)液體的體積分?jǐn)?shù),因此有著較好的質(zhì)量守恒特性;但是,在兩相界面陡峭的條件下,LSM方法能夠更好地捕捉界面。但是LSM方法守恒性較差,隨著數(shù)值迭代的積累,在界面變化時(shí)液體體積可能會(huì)丟失,并且它的計(jì)算開銷較之VOF方法更大。綜上所述,本文選擇了VOF方法進(jìn)行仿真。

    本文的研究對(duì)象為相對(duì)于靜止空氣有較大的初速度的大水滴,由于與環(huán)境溫差小并且液滴尺寸相對(duì)較小,在研究過程中可以忽略傳熱傳質(zhì)現(xiàn)象,湍流影響和重力作用[15],因此水滴運(yùn)動(dòng)變形的流場(chǎng)可以視為二維軸對(duì)稱不可壓縮的層流流場(chǎng),求解相應(yīng)的N-S方程和連續(xù)方程。動(dòng)量方程的控制方程,即考慮表面張力項(xiàng)的N-S方程為:

    式中F為液體體積分?jǐn)?shù),定義如下:

    由F定義可知,式(1)中

    ρ(F)=ρlF+ρg(1-F)(3)

    ρl、ρg分別是液相和氣相的密度。

    同理可得

    μ(F)=μlF+μg(1-F)(4)

    式中μl、μg分別是液相和氣相的動(dòng)力黏度。u為速度矢量,t是時(shí)間,應(yīng)變率張量由D=(u+uT)/2給出。式(1)中最后一項(xiàng)即為表面張力項(xiàng),σ為表面張力系數(shù),而κ為氣液界面曲率,δs為表面狄拉克δ函數(shù)。

    數(shù)值模擬的計(jì)算模型示意圖如圖1所示,流場(chǎng)為1.5×40 mm的矩形,x軸為軸對(duì)稱條件,其他三個(gè)邊界均為常壓壓力出口條件。初始條件為流場(chǎng)空氣靜止且為常溫常壓(20 ℃,一個(gè)大氣壓),一個(gè)直徑為0.62 mm的初速度為u0的球形水滴出現(xiàn)在流場(chǎng)中,并從距入口1.5 mm處開始運(yùn)動(dòng)。

    圖1 水滴在靜止空氣中運(yùn)動(dòng)破碎示意圖Fig.1 Schematic of a high-speed droplet movement and breakup in still air

    影響水滴變形、破碎的參數(shù)主要有水滴初始直徑d0,水滴-氣流相對(duì)速度u,水滴-氣流表面張力σ,水滴、氣流的密度ρl、ρg,動(dòng)力黏度μl、μg等。進(jìn)行無量綱化后得到無量綱參數(shù)表達(dá)式如下:

    We=ρgu2d0/σ(5)

    Re=ρgud0/μg(6)

    采用上述無量綱參數(shù),可列出算例條件見表1。

    由于VOF模型對(duì)于網(wǎng)格質(zhì)量和密度的依賴性較強(qiáng),本文選用了三套網(wǎng)格進(jìn)行獨(dú)立性驗(yàn)證。在水滴運(yùn)動(dòng)變形的區(qū)域網(wǎng)格尺寸分別是25 μm,10 μm和5 μm,網(wǎng)格數(shù)目分別是25萬,60萬和120萬。驗(yàn)算表明,當(dāng)網(wǎng)格尺寸較大時(shí),液相界面較粗糙,不能準(zhǔn)確把握變化迅速的水滴變形破碎的特征。但是當(dāng)網(wǎng)格過密時(shí),增加了大量的計(jì)算時(shí)間且計(jì)算結(jié)果沒有明顯的改善,綜合比較之后,選擇了最小尺寸為10 μm的網(wǎng)格進(jìn)行計(jì)算。

    表1 算例初始條件設(shè)置(Oh=0.00473)Table 1 Summary of the simulation initial conditions (Oh=0.00473)

    本文對(duì)壓力耦合方程的求解采用了半隱方法,即SIMPLEC算法,在文獻(xiàn)[19]中對(duì)二維層流流動(dòng)的計(jì)算結(jié)果表明,SIMPLEC算法的收斂特性遠(yuǎn)優(yōu)于SIMPLE算法。壓強(qiáng)計(jì)算采用了PRESTO格式,其他物理量采用QUICK格式離散。

    2 計(jì)算結(jié)果與分析

    2.1 大水滴典型袋狀破碎形態(tài)特征

    在不同的We數(shù)下,液滴的破碎呈現(xiàn)出不同的形態(tài)特征,Guildenbecher[6]對(duì)這些破碎模式做了總結(jié)性的分類,它們分別是振動(dòng)破碎、袋狀破碎、復(fù)合模式破碎、剪切破碎和破壞性破碎??梢杂杀?看出,這些破碎類型之間并不存在We數(shù)的明確界限,不同的研究人員往往給出了不同的結(jié)果。這主要是因?yàn)槟壳皩?duì)破碎模型的定義只是定性而不是定量的,并且由于實(shí)驗(yàn)條件的不同(激波管、風(fēng)洞或者噴嘴),都對(duì)水滴破碎產(chǎn)生了較大的影響。

    對(duì)于表1所列出的速度范圍,即流體相對(duì)速度較低時(shí),大部分?jǐn)?shù)值分析選擇了層流流動(dòng)進(jìn)行計(jì)算,此時(shí)對(duì)水滴破碎的分析往往基于不穩(wěn)定性理論。對(duì)這一部分的分析往往考慮以下因素:表面張力作用(R-P不穩(wěn)定性),水滴被環(huán)境空氣加速/減速(R-T不穩(wěn)定性),以及在液滴邊緣表現(xiàn)出的剪切不穩(wěn)定性(Kelvin-Helmholtz不穩(wěn)定性,簡(jiǎn)稱K-H不穩(wěn)定性)[20]。最后一項(xiàng)常見于剪切破碎模式下的分析,本文暫不考慮。在對(duì)于非紊流袋狀破碎的研究中,R-T不穩(wěn)定性被廣泛應(yīng)用于解釋實(shí)驗(yàn)結(jié)果[21-23],本文中也將R-T不穩(wěn)定性應(yīng)用于解釋袋狀破碎初期的現(xiàn)象。

    表2 不同的液滴破碎模式We數(shù)范圍劃分(Oh數(shù)遠(yuǎn)小于0.1)Table 2 Different breakup regimes as a function of We(Oh is far less than 0.1)

    因此,本文針對(duì)不同研究人員確定的袋狀破碎區(qū)間,選擇了We數(shù)為14~22的水滴,作為典型袋狀破碎算例進(jìn)行形態(tài)特征的分析和研究。如圖2所示,可以觀察到水滴在靜止空氣中以中等速度運(yùn)動(dòng),并逐漸發(fā)生變形破碎的過程。從中可以明顯地觀察到破碎的典型特征和發(fā)生的時(shí)間,在接下來的章節(jié)中將對(duì)這些特征和各類重要的瞬態(tài)性能進(jìn)行分析和研究。 在水滴破碎過程中通常使用特征時(shí)間tc來對(duì)各種瞬態(tài)性能進(jìn)行無量綱化分析,tc的表達(dá)式如下:

    圖2 We=16水滴隨時(shí)間變形破碎過程(水滴運(yùn)動(dòng)方向從左至右)各圖對(duì)應(yīng)無量綱時(shí)間為:t/tc=a) 0, b) 0.33, c) 0.55, d) 0.77, e) 1.22, f) 1.77, g) 2.12, h) 2.19.此處tc=452 μs.Fig.2 Numerical simulation for droplet moving from left to right with We=16 at t/tc=a) 0, b) 0.33, c) 0.55, d) 0.77, e) 1.22, f) 1.77, g) 2.12, h) 2.19. (tc=452 μs)

    根據(jù)仿真中觀察到的水滴變形和破碎的形態(tài)特征,本文首次提出,將袋狀破碎分為五個(gè)特征階段,具體如下:

    1) 梯形變形階段

    圖2(a-c)的過程,水滴由球形開始發(fā)生形變。由于水滴運(yùn)動(dòng)方向從左至右,在迎風(fēng)側(cè)的變化較為迅速,由圖3(b)可見,水滴迎風(fēng)側(cè)靠近駐點(diǎn)的位置由于該處的高壓被擠壓變平。而與之相對(duì)應(yīng)的,水滴在背風(fēng)側(cè)的變化較為緩慢,仍基本保持球體不變,這與早期的實(shí)驗(yàn)所觀察到的現(xiàn)象一致[4]。此后由于慣性力和氣動(dòng)阻力的共同作用,水滴整體逐漸被擠壓變平,在無量綱時(shí)間t/tc=0.55時(shí),可以觀察到明顯的梯形結(jié)構(gòu)。

    2) 延伸發(fā)展階段

    圖2(c-d)是水滴縱向延伸并且厚度均勻的階段。由圖3(c)可見,在運(yùn)動(dòng)前方駐點(diǎn)處和液體內(nèi)部壓力作用下,水滴在x方向受力逐漸被擠壓。又由于水滴邊緣處空氣的速度增大,靜壓減小,水滴在y方向不斷延展,最終形成了一個(gè)厚度均勻的圓餅形水滴。這個(gè)由球體逐漸變?yōu)楸馄降倪^程在無量綱時(shí)間t/tc=0.77時(shí),可以顯著觀察到。

    3) 水環(huán)形成階段

    圖2(d-f)是水環(huán)結(jié)構(gòu)出現(xiàn)的階段。水滴由厚度均勻的餅狀向中心較薄邊緣較厚的袋狀結(jié)構(gòu)發(fā)展,中部的液體向邊緣移動(dòng),水環(huán)不斷變厚并且作為支撐結(jié)構(gòu)承力,在t/tc=1.77時(shí)可以明顯觀察到水環(huán)的形成。而水滴中心部位反速度方向延展,形成袋狀薄膜。這主要是由于袋狀凹陷處,氣流在駐點(diǎn)處的高壓區(qū)擠壓液體在水滴內(nèi)部由赤道位置向極點(diǎn)位置運(yùn)動(dòng)導(dǎo)致。并且由圖3(e)可以看出,液滴的表面出現(xiàn)了壓力波動(dòng),然而這些擾動(dòng)還處在表面,由于袋狀薄膜的厚度未達(dá)到臨界值,這些微小的擾動(dòng)沒有穿透液滴。

    4) 袋狀破碎階段

    圖2(f-g)是袋狀結(jié)構(gòu)不斷向后延展變薄,最終破碎斷裂與水環(huán)分離的過程。這個(gè)過程可以用R-T不穩(wěn)定性來解釋。R-T不穩(wěn)定現(xiàn)象可以簡(jiǎn)單表述為:當(dāng)兩種不同密度流體的分界面做加速運(yùn)動(dòng)且加速度方向垂直于界面時(shí),此界面是否穩(wěn)定取決于加速度的方向是否從大密度流體指向小密度流體。由于受到氣動(dòng)阻力,并且水滴所受加速度方向由輕質(zhì)流體(空氣)指向重質(zhì)流體(水),就會(huì)出現(xiàn)R-T不穩(wěn)定現(xiàn)象,表現(xiàn)為袋狀薄膜上出現(xiàn)壓力的擾動(dòng),見圖3(f-g)。

    這個(gè)過程可以用R-T不穩(wěn)定性來解釋,如文獻(xiàn)中所述[21],假定傳統(tǒng)的二維不穩(wěn)定性理論可以擴(kuò)展至三維條件,根據(jù)R-T理論,可以證得在液滴表面出現(xiàn)特定波長(zhǎng)的擾動(dòng)時(shí),擾動(dòng)振幅的增長(zhǎng)速度最快,這一特定波長(zhǎng)見式(9):

    其中f是盤狀液滴在靜止空氣中的加速度,其表達(dá)式如式(10)所示:

    這種特定波長(zhǎng)的擾動(dòng)迅速在液滴的迎風(fēng)面上形成并發(fā)展,擾動(dòng)的幅度也逐漸增大,形成孔洞最終穿透液膜并迅速擴(kuò)大.在無量綱時(shí)間t/tc=2.12時(shí),可以觀察到膜狀結(jié)構(gòu)破裂,與水環(huán)分離。

    5) 子液滴回縮階段

    圖2(g-h)是袋狀薄膜破碎成更小的子液滴的過程。如圖2(h)所示,可以觀察到中心薄膜上和與水環(huán)相接處,均出現(xiàn)了正弦波的形狀。

    按照實(shí)驗(yàn)觀測(cè)和R-P不穩(wěn)定性的假設(shè)[20],受到恒定加速度作用的黏性流體在運(yùn)動(dòng)過程中會(huì)由于不穩(wěn)定性誘導(dǎo)而破碎。即對(duì)于特定形狀的流體,存在一個(gè)特征頻率(fc)的擾動(dòng),這個(gè)擾動(dòng)的周期(λc)一定而成長(zhǎng)率最大,在運(yùn)動(dòng)過程中占據(jù)統(tǒng)治地位,導(dǎo)致流體發(fā)生破碎。

    在這一假設(shè)條件下對(duì)實(shí)驗(yàn)和仿真結(jié)果進(jìn)行分析:Jain[18]在液滴破碎實(shí)驗(yàn)中觀察到子液滴的粒徑呈雙峰分布,分別是水環(huán)破碎和袋狀薄膜破碎產(chǎn)生的,即袋狀薄膜破碎后產(chǎn)生的子液滴粒徑均一。此外,如圖2(g)所示,本文仿真結(jié)果中同樣觀察到,臨界破碎時(shí)薄膜表面出現(xiàn)的正弦波是等振幅的,與前文實(shí)驗(yàn)觀察結(jié)論結(jié)合,可以由理論計(jì)算預(yù)測(cè)袋狀破碎后子水滴的粒徑分布。

    在無量綱時(shí)間t/tc=2.19時(shí),可見到斷裂后的正弦周期狀液面不斷回縮,在表面張力的作用下呈表面積最小化的趨勢(shì),形成圓形液滴。此外,雖然袋狀薄膜破碎,但外圍的水環(huán)仍保持了原有形狀不變。

    圖3 We=16水滴袋裝破碎壓力云圖(水滴運(yùn)動(dòng)方向從左至右)(tc =452 μs)Fig.3 The pressure nephogram of water drop (We=16) during bag breakup(tc =452 μs)

    2.2 大水滴典型袋狀破碎瞬態(tài)參數(shù)分析

    2.2.1 水滴形變及破碎時(shí)間分析

    水滴袋狀破碎過程中的形變和破碎所需的時(shí)間是值得關(guān)注的重要參數(shù)。在前人的工作中,Gordon[16]建立了簡(jiǎn)易的圓柱體水滴模型用于預(yù)測(cè)水滴破碎時(shí)間,模型假設(shè)液體為剛性,氣動(dòng)力將水滴從柱體空間內(nèi)推出即為水滴破碎;此后又有Rourke&Amsden提出基于彈簧振子系統(tǒng)的TAB模型[17];本課題組創(chuàng)立了將變形過程等效為橢球形水滴的BTB模型[14],認(rèn)為在變形過程中,當(dāng)半個(gè)水滴的質(zhì)心在直徑方向上發(fā)生一定的相對(duì)位移后破碎。BTB模型引入修正系數(shù)Cm,將相對(duì)變形量D=d/d0表示成We數(shù)的函數(shù),如式(11)所示。

    (D-1+D5-2D-4)=80CmWe(0

    本文對(duì)袋狀破碎進(jìn)行仿真研究,并且將仿真結(jié)果與基于BTB模型的計(jì)算結(jié)果和前人的實(shí)驗(yàn)結(jié)果[3-4,18]進(jìn)行比對(duì),如圖4所示。其中包括來自Krzeczkowski(We=13.5)、Chou &Faeth (We=15)和Mohit Jain(We=20)。圖4中給出了不同We數(shù)下袋狀破碎水滴隨著時(shí)間形變的曲線(實(shí)驗(yàn)為系列點(diǎn)),并以箭頭標(biāo)出了相應(yīng)的水滴破碎點(diǎn)。其中y軸坐標(biāo)為水滴無量綱形變量,即水滴迎風(fēng)直徑d與初始直徑d0的比值(在水滴破碎后,迎風(fēng)直徑視為水環(huán)的外徑);x軸坐標(biāo)為無量綱時(shí)間,即實(shí)際時(shí)間t與特征時(shí)間tc之比。圖中水滴We數(shù)為16~22,Oh均小于0.005。

    圖4 水滴相對(duì)變形量d/d0與無量綱時(shí)間t/tc曲線圖Fig.4 Drop cross-stream diameter as a function of time during bag breakup

    首先需要明確的是,水滴破碎變形過程主要受到以下幾種力的作用:氣動(dòng)力,表面張力和黏性力。其中氣動(dòng)力是發(fā)生變形的主因,表面張力和水滴內(nèi)部的黏性力是減緩和阻止水滴變形的作用力。式(5)中給出的常用參數(shù)正是對(duì)這三個(gè)作用力進(jìn)行了無量綱化:We數(shù)代表氣動(dòng)力比表面張力,Oh數(shù)代表了水滴內(nèi)部黏性力比表面張力。在本文關(guān)注的袋狀破碎大水滴領(lǐng)域,水滴Oh數(shù)遠(yuǎn)小于0.01,We>10,即黏性力遠(yuǎn)小于水滴表面張力的量級(jí),而氣動(dòng)力又遠(yuǎn)大于表面張力,即在水滴運(yùn)動(dòng)的短暫時(shí)間間隔內(nèi)動(dòng)量來不及傳遞,這意味著水滴將先變形而不是改變速度。因此本文分析的重點(diǎn)放在了氣動(dòng)力和水滴形變上,兼顧表面張力的影響。

    由圖中三組曲線(仿真組、實(shí)驗(yàn)組和擬合組)可以看出,在每組中曲線中,隨著We數(shù)的增加,水滴的相對(duì)變形量也不斷增大。此外,在仿真組和實(shí)驗(yàn)組的曲線中可以觀察到,水滴破碎后曲線斜率明顯增大,即水滴殘余的水環(huán)迅速擴(kuò)大。聯(lián)系圖3(g,h)的壓力云圖可以分析得出,水滴袋狀薄膜向后凹陷處,存在一個(gè)迎風(fēng)駐點(diǎn)處的高壓區(qū),而水滴邊緣處由于氣流速度增大而靜壓減小,導(dǎo)致液滴破碎后外側(cè)水環(huán)在y方向受力,有一個(gè)側(cè)向的加速度,快速膨脹。但隨著水滴的破碎,這個(gè)壓差迅速消失,水環(huán)保持一個(gè)勻速的膨脹速度直到破碎。本文之所以關(guān)注此處水環(huán)的變化,是因?yàn)樵谄淦扑楹笮纬傻乃坞m然數(shù)量較少,但卻是子水滴群中體積分?jǐn)?shù)最大的部分,在子水滴群中占統(tǒng)治地位。

    對(duì)于水滴袋狀破碎過程中的變形,在Chou和Faeth的實(shí)驗(yàn)文獻(xiàn)中,不僅給出了實(shí)驗(yàn)數(shù)據(jù),還提出了相關(guān)的經(jīng)驗(yàn)公式,指出在水滴變形初期過程,相對(duì)變形量與無量綱時(shí)間呈線性正相關(guān),此后為二次函數(shù)關(guān)系,具體公式如式(12)所示:

    在對(duì)本文仿真數(shù)據(jù)進(jìn)行擬合后,發(fā)現(xiàn)在We=16的條件下,結(jié)果與經(jīng)驗(yàn)公式吻合很好,如式(13),且線性段相關(guān)系數(shù)R2=0.995:

    但隨著We數(shù)的增長(zhǎng),線性段曲線斜率增大,逐漸偏離經(jīng)驗(yàn)公式,分析數(shù)據(jù)后認(rèn)為這一經(jīng)驗(yàn)公式更適用于We較低(We=16)的袋裝破碎。并且從圖中數(shù)組曲線可以看出,對(duì)于相對(duì)形變仿真的結(jié)果在全范圍內(nèi)都與試驗(yàn)數(shù)據(jù)吻合較好;而BTB模型主要關(guān)注的是水滴在袋狀薄膜臨界破碎時(shí)期變形量和破碎時(shí)間,因此在破碎后期與實(shí)驗(yàn)數(shù)據(jù)貼合更好。

    在水滴破碎時(shí)間的預(yù)測(cè)上,BTB模型效果很好:雖然比起實(shí)驗(yàn)結(jié)果略微偏大,但最大相對(duì)誤差小于15%,并且在We=20的條件下與實(shí)驗(yàn)結(jié)果(圖5)極為貼近。而仿真數(shù)據(jù)則普遍偏小,對(duì)這一結(jié)果進(jìn)行分析,可以參照水滴破碎試驗(yàn)結(jié)果。考慮到實(shí)際實(shí)驗(yàn)中水滴臨界破碎時(shí),不斷拉伸的袋狀薄膜厚度不斷減小,曲率半徑隨之減小,表面張力不斷增加,與氣動(dòng)力達(dá)到平衡,保持形態(tài)的穩(wěn)定不發(fā)生破碎。直到變形過程中,液膜曲率半徑與水環(huán)半徑一致不再變化,氣動(dòng)力的擾動(dòng)克服了表面張力,薄膜才發(fā)生破碎。但在仿真中,所使用的網(wǎng)格步長(zhǎng)為10 μm,大于了臨界破碎的水膜厚度,在計(jì)算中無法實(shí)現(xiàn)水膜的急劇變薄的過程,導(dǎo)致較之實(shí)驗(yàn)過早地發(fā)生了破碎。對(duì)此應(yīng)做的改進(jìn)是結(jié)合計(jì)算經(jīng)濟(jì)性,在運(yùn)算中使用自適應(yīng)網(wǎng)格,在薄膜處進(jìn)一步細(xì)化網(wǎng)格,才能更好地表現(xiàn)出水膜破碎的完整形態(tài)。

    圖5 水滴在We=20條件下的袋狀破碎試驗(yàn)照片[18](氣流方向從右至左, tc=110 μs)Fig.5 Experimental observations[18] for We=20(tc=110 μs)

    此外從水滴破碎點(diǎn)中可以看出,隨著We數(shù)的增長(zhǎng),三組數(shù)據(jù)都顯示出一致的趨勢(shì):水滴破碎的時(shí)間減小,而臨界破碎的相對(duì)變形量明顯增加。

    2.2.2 水滴速度與阻力系數(shù)

    在研究袋狀破碎過程中大水滴的速度和受力時(shí),常常要考慮形變帶來的影響。顯而易見,隨著水滴在不同階段中形狀的改變和迎風(fēng)直徑的增大,受到的阻力比起初始球形有著明顯的增大,在對(duì)阻力和速度進(jìn)行分析時(shí),常常會(huì)引入阻力系數(shù)CD,將形變與受力以經(jīng)驗(yàn)系數(shù)的形式聯(lián)系起來,CD的表達(dá)式如式(14)所示,其中u是水滴質(zhì)心的瞬時(shí)速度,a是計(jì)算得到的加速度。

    圖6給出了不同條件下水滴主體無量綱速度變化量(u0-u)/u0和無量綱時(shí)間t/tc的關(guān)系,以及阻力系數(shù)CD曲線。其中除本文仿真結(jié)果外,還對(duì)照Chou和Faeth[4]的實(shí)驗(yàn)結(jié)果(We=15)以及Jain[18]的三維數(shù)值實(shí)驗(yàn)結(jié)果,驗(yàn)證仿真正確性。

    (a)

    (b)

    從圖6(a)中可以觀察到,液滴質(zhì)心無量綱速度在運(yùn)動(dòng)初期和末期變化平緩,在t/tc=0.6~1.7之間,發(fā)生較為劇烈的變化,即在水滴呈盤形的階段,速度的相對(duì)改變量約為7%。如圖6(b)所示,阻力系數(shù)CD隨著水滴的扁平化,也就是水滴迎風(fēng)直徑的增長(zhǎng)而不斷增大,并且在迎風(fēng)直徑最大時(shí)達(dá)到峰值(CD~1.8),此后隨著袋狀薄膜開始形成,迎風(fēng)直徑不再顯著變化,而水滴的背風(fēng)面形成了光滑的過渡,水滴前后壓差減小,CD顯著降低。

    此外,在起始階段球形液滴的阻力系數(shù)約為0.58,與理論值和實(shí)驗(yàn)數(shù)據(jù)吻合較好。但在峰值處,根據(jù)Chou & Faeth[4]的阻力系數(shù)擬合曲線,袋狀破碎的阻力系數(shù)最大值約為1.6,而仿真的結(jié)果(CD,max~1.8)比起實(shí)驗(yàn)結(jié)果偏大,但與不可壓縮流(Re<2000)中盤形的阻力系數(shù)理論值1.7相近[24]。

    3 結(jié) 論

    本文應(yīng)用VOF方法,對(duì)大水滴典型袋狀破碎過程中的特性進(jìn)行了探究,得到的主要結(jié)論如下:

    1) 袋狀破碎分為五個(gè)特征階段:梯形變形階段、延伸發(fā)展階段、水環(huán)形成階段、袋狀破碎階段和子液滴形成階段;

    2) 水膜的破碎可以用R-T不穩(wěn)定性來解釋,子液滴的粒徑可以應(yīng)用R-P不穩(wěn)定性進(jìn)行預(yù)測(cè);

    3) 水滴破碎后殘余的水環(huán)迅速擴(kuò)大。這是由于水滴袋狀薄膜凹陷處,即迎風(fēng)駐點(diǎn)處存在一個(gè)高壓區(qū),水環(huán)在y向受力導(dǎo)致的;

    4) 水滴變形初期過程,相對(duì)變形量與無量綱時(shí)間呈線性正相關(guān),此后為二次函數(shù)關(guān)系的關(guān)系式,與經(jīng)驗(yàn)公式吻合很好;

    5) 水滴破碎時(shí)間隨著We數(shù)的增長(zhǎng)減小,而臨界破碎的相對(duì)變形量明顯增加。

    6) 阻力系數(shù)CD隨著水滴迎風(fēng)直徑的增長(zhǎng)而不斷增大,最大值約為1.8,此后顯著降低。

    猜你喜歡
    水環(huán)無量不穩(wěn)定性
    烏雷:無量之物
    參觀水環(huán)中心
    劉少白
    藝術(shù)品(2020年8期)2020-10-29 02:50:02
    可壓縮Navier-Stokes方程平面Couette-Poiseuille流的線性不穩(wěn)定性
    論書絕句·評(píng)謝無量(1884—1964)
    炳靈寺第70 窟無量壽經(jīng)變辨識(shí)
    西藏研究(2017年3期)2017-09-05 09:45:07
    增強(qiáng)型體外反搏聯(lián)合中醫(yī)辯證治療不穩(wěn)定性心絞痛療效觀察
    前列地爾治療不穩(wěn)定性心絞痛療效觀察
    制何首烏中二苯乙烯苷對(duì)光和熱的不穩(wěn)定性
    中成藥(2014年11期)2014-02-28 22:29:49
    水環(huán)壓縮機(jī)組在乙炔生產(chǎn)中的應(yīng)用
    国产精品偷伦视频观看了| 91成年电影在线观看| 99国产极品粉嫩在线观看| 午夜免费成人在线视频| 青春草视频在线免费观看| 欧美97在线视频| 在线观看人妻少妇| 另类亚洲欧美激情| 日韩,欧美,国产一区二区三区| 岛国毛片在线播放| 叶爱在线成人免费视频播放| 欧美精品人与动牲交sv欧美| 久久久久国内视频| 夜夜夜夜夜久久久久| 老司机影院毛片| 91av网站免费观看| 欧美精品人与动牲交sv欧美| 精品国产乱码久久久久久小说| 日韩一卡2卡3卡4卡2021年| 亚洲国产av影院在线观看| 操美女的视频在线观看| 人人澡人人妻人| 国产成人精品久久二区二区免费| 97人妻天天添夜夜摸| 国产日韩欧美视频二区| 亚洲全国av大片| 999久久久国产精品视频| 亚洲成人免费av在线播放| 久久亚洲精品不卡| 国产精品.久久久| 亚洲九九香蕉| 制服诱惑二区| 久久久久国产一级毛片高清牌| 狂野欧美激情性xxxx| 欧美中文综合在线视频| 色播在线永久视频| 九色亚洲精品在线播放| 精品人妻一区二区三区麻豆| 欧美老熟妇乱子伦牲交| 国产不卡av网站在线观看| 亚洲精品国产精品久久久不卡| 中文字幕人妻丝袜制服| 久9热在线精品视频| 亚洲国产av影院在线观看| 人妻久久中文字幕网| 欧美大码av| 欧美另类亚洲清纯唯美| 亚洲精品久久午夜乱码| 黄色视频在线播放观看不卡| 麻豆乱淫一区二区| 国产亚洲精品第一综合不卡| 热99久久久久精品小说推荐| 黄频高清免费视频| 午夜精品久久久久久毛片777| 亚洲av电影在线观看一区二区三区| 国产片内射在线| 搡老熟女国产l中国老女人| 亚洲av成人不卡在线观看播放网 | 亚洲国产av影院在线观看| 99国产精品99久久久久| 国产野战对白在线观看| 大型av网站在线播放| 欧美一级毛片孕妇| 九色亚洲精品在线播放| 在线精品无人区一区二区三| 如日韩欧美国产精品一区二区三区| 久久综合国产亚洲精品| 久久久久久久大尺度免费视频| 欧美日本中文国产一区发布| 精品久久蜜臀av无| 一本久久精品| 精品一品国产午夜福利视频| 欧美日韩成人在线一区二区| 窝窝影院91人妻| 国产高清视频在线播放一区 | 亚洲欧美清纯卡通| 国产一区有黄有色的免费视频| 国产免费福利视频在线观看| av网站在线播放免费| 亚洲va日本ⅴa欧美va伊人久久 | 首页视频小说图片口味搜索| 国产成+人综合+亚洲专区| 丝瓜视频免费看黄片| 在线天堂中文资源库| 美女主播在线视频| av电影中文网址| 一个人免费看片子| 欧美精品一区二区大全| av天堂久久9| 欧美日韩黄片免| 久久久久精品人妻al黑| 日本vs欧美在线观看视频| 夜夜骑夜夜射夜夜干| 亚洲国产毛片av蜜桃av| 男女边摸边吃奶| 国产片内射在线| 建设人人有责人人尽责人人享有的| 18禁国产床啪视频网站| 欧美日韩福利视频一区二区| 亚洲精品国产av成人精品| 亚洲国产毛片av蜜桃av| 69av精品久久久久久 | 老熟妇乱子伦视频在线观看 | 成年女人毛片免费观看观看9 | 国产激情久久老熟女| a 毛片基地| 午夜福利一区二区在线看| 欧美日韩成人在线一区二区| 99久久综合免费| 97在线人人人人妻| 国产福利在线免费观看视频| 久久影院123| 999久久久精品免费观看国产| 精品少妇黑人巨大在线播放| 一进一出抽搐动态| 欧美精品啪啪一区二区三区 | 一本—道久久a久久精品蜜桃钙片| 久久 成人 亚洲| 亚洲成人免费电影在线观看| 飞空精品影院首页| 亚洲精品国产一区二区精华液| 人人妻人人添人人爽欧美一区卜| 水蜜桃什么品种好| 巨乳人妻的诱惑在线观看| 一二三四社区在线视频社区8| 日韩大码丰满熟妇| 成人国产av品久久久| 制服诱惑二区| 丰满迷人的少妇在线观看| 一级a爱视频在线免费观看| h视频一区二区三区| 国产精品久久久久久人妻精品电影 | 亚洲国产毛片av蜜桃av| 亚洲av电影在线进入| tocl精华| 50天的宝宝边吃奶边哭怎么回事| 一边摸一边抽搐一进一出视频| 日本欧美视频一区| 肉色欧美久久久久久久蜜桃| 午夜老司机福利片| 久久青草综合色| 亚洲欧美激情在线| 另类亚洲欧美激情| 国产精品偷伦视频观看了| 久久久国产一区二区| 美国免费a级毛片| 丰满饥渴人妻一区二区三| 在线永久观看黄色视频| 久久久久国内视频| 免费av中文字幕在线| 中文字幕精品免费在线观看视频| 精品少妇一区二区三区视频日本电影| 日本91视频免费播放| av有码第一页| 久久性视频一级片| 国产精品亚洲av一区麻豆| 久久午夜综合久久蜜桃| 成年人午夜在线观看视频| 亚洲av成人不卡在线观看播放网 | 中文字幕最新亚洲高清| 亚洲va日本ⅴa欧美va伊人久久 | 十八禁网站免费在线| 久久久久久久久久久久大奶| 色婷婷av一区二区三区视频| 国精品久久久久久国模美| 在线观看www视频免费| 两个人免费观看高清视频| 国产有黄有色有爽视频| www.精华液| 美女视频免费永久观看网站| av网站免费在线观看视频| 黑人欧美特级aaaaaa片| 欧美av亚洲av综合av国产av| 嫁个100分男人电影在线观看| 中文字幕高清在线视频| 久久久欧美国产精品| 国产国语露脸激情在线看| av不卡在线播放| 亚洲第一青青草原| 狂野欧美激情性xxxx| 侵犯人妻中文字幕一二三四区| 国产在线视频一区二区| 国产亚洲欧美精品永久| 午夜精品国产一区二区电影| 性色av乱码一区二区三区2| 国产片内射在线| 欧美日韩亚洲高清精品| 国产福利在线免费观看视频| 脱女人内裤的视频| 国产精品国产av在线观看| 99久久国产精品久久久| 国产亚洲av高清不卡| 五月开心婷婷网| 亚洲成国产人片在线观看| 精品国产国语对白av| 丝袜喷水一区| 69精品国产乱码久久久| 王馨瑶露胸无遮挡在线观看| 亚洲精品一二三| 热99re8久久精品国产| 亚洲欧美清纯卡通| 亚洲成av片中文字幕在线观看| 三级毛片av免费| 亚洲va日本ⅴa欧美va伊人久久 | 欧美另类一区| 精品人妻1区二区| 18在线观看网站| 窝窝影院91人妻| 日韩电影二区| 亚洲伊人久久精品综合| 每晚都被弄得嗷嗷叫到高潮| 免费看十八禁软件| 欧美 亚洲 国产 日韩一| 国产精品二区激情视频| 亚洲成国产人片在线观看| 久久久精品区二区三区| 人成视频在线观看免费观看| 在线观看舔阴道视频| 另类亚洲欧美激情| 亚洲精品乱久久久久久| 久久久久久久国产电影| 超色免费av| 精品久久久精品久久久| 亚洲精品美女久久av网站| 大片电影免费在线观看免费| 久久精品国产a三级三级三级| 亚洲精华国产精华精| 在线观看免费高清a一片| 午夜福利乱码中文字幕| 女人精品久久久久毛片| 亚洲 国产 在线| 纵有疾风起免费观看全集完整版| 9191精品国产免费久久| 亚洲精品国产精品久久久不卡| 成人三级做爰电影| 狠狠精品人妻久久久久久综合| 日韩人妻精品一区2区三区| 精品高清国产在线一区| 午夜福利一区二区在线看| 日本撒尿小便嘘嘘汇集6| 亚洲五月色婷婷综合| 成在线人永久免费视频| 国产免费视频播放在线视频| 亚洲男人天堂网一区| 国产精品亚洲av一区麻豆| 老汉色∧v一级毛片| 久久精品国产亚洲av高清一级| 看免费av毛片| 久久精品aⅴ一区二区三区四区| 日韩制服骚丝袜av| 老司机靠b影院| 69av精品久久久久久 | 精品亚洲成国产av| 欧美日韩福利视频一区二区| 国产精品1区2区在线观看. | 精品视频人人做人人爽| 女人精品久久久久毛片| 亚洲伊人色综图| 一区二区三区精品91| 在线 av 中文字幕| 叶爱在线成人免费视频播放| 大片免费播放器 马上看| 欧美午夜高清在线| 成人三级做爰电影| 国产成人系列免费观看| 别揉我奶头~嗯~啊~动态视频 | 日韩大片免费观看网站| 亚洲 国产 在线| 一区二区三区精品91| 国产老妇伦熟女老妇高清| 女性生殖器流出的白浆| 青青草视频在线视频观看| 国产真人三级小视频在线观看| 热99re8久久精品国产| 性色av一级| 男女下面插进去视频免费观看| 午夜日韩欧美国产| 国产欧美亚洲国产| 亚洲情色 制服丝袜| 在线观看免费视频网站a站| 国产成人av教育| 黄网站色视频无遮挡免费观看| 满18在线观看网站| 亚洲专区字幕在线| 国产一区二区三区在线臀色熟女 | 日韩欧美免费精品| 国产99久久九九免费精品| 国产欧美日韩一区二区三 | 亚洲精品一区蜜桃| 视频区图区小说| 久久久久久免费高清国产稀缺| 老鸭窝网址在线观看| 人妻人人澡人人爽人人| 久久女婷五月综合色啪小说| 中文字幕av电影在线播放| 国产伦理片在线播放av一区| 正在播放国产对白刺激| 欧美午夜高清在线| 熟女少妇亚洲综合色aaa.| 在线观看人妻少妇| 女人精品久久久久毛片| 国产亚洲精品第一综合不卡| 亚洲av男天堂| videos熟女内射| 99精品久久久久人妻精品| 亚洲中文字幕日韩| 亚洲精华国产精华精| 热99国产精品久久久久久7| 欧美日韩中文字幕国产精品一区二区三区 | 国产成人欧美在线观看 | 亚洲精品国产区一区二| 欧美成狂野欧美在线观看| 亚洲精品成人av观看孕妇| 51午夜福利影视在线观看| 国产成人精品久久二区二区免费| 90打野战视频偷拍视频| 日日爽夜夜爽网站| 国产在线视频一区二区| 人人澡人人妻人| 欧美午夜高清在线| 国产高清国产精品国产三级| 亚洲成人国产一区在线观看| 国产片内射在线| av欧美777| 国产男人的电影天堂91| 搡老岳熟女国产| 亚洲伊人久久精品综合| 新久久久久国产一级毛片| 18禁黄网站禁片午夜丰满| 少妇粗大呻吟视频| 久久久国产成人免费| 麻豆乱淫一区二区| 在线观看免费日韩欧美大片| 另类精品久久| 免费观看a级毛片全部| 亚洲欧美激情在线| 欧美日韩成人在线一区二区| 性色av一级| 亚洲人成电影观看| 在线观看www视频免费| 久久久久国内视频| 激情视频va一区二区三区| 一级,二级,三级黄色视频| 男人爽女人下面视频在线观看| 日韩欧美免费精品| 免费日韩欧美在线观看| 欧美精品一区二区大全| 久久久久国产精品人妻一区二区| 18在线观看网站| 亚洲欧美日韩另类电影网站| 国产黄色免费在线视频| 亚洲精品成人av观看孕妇| 麻豆av在线久日| 一级a爱视频在线免费观看| 电影成人av| 丰满迷人的少妇在线观看| 宅男免费午夜| 精品第一国产精品| 国产一区二区三区综合在线观看| 久久 成人 亚洲| 国产一区二区激情短视频 | 最近最新中文字幕大全免费视频| 99国产精品一区二区三区| 亚洲九九香蕉| 老司机影院成人| 性色av乱码一区二区三区2| 高清欧美精品videossex| 黄色视频,在线免费观看| 欧美日韩国产mv在线观看视频| 久久久国产精品麻豆| 久久 成人 亚洲| 日韩熟女老妇一区二区性免费视频| 免费一级毛片在线播放高清视频 | 黄色片一级片一级黄色片| 亚洲 国产 在线| 精品国产一区二区久久| 99国产极品粉嫩在线观看| 亚洲伊人久久精品综合| 国产在线观看jvid| 欧美日韩亚洲国产一区二区在线观看 | 精品第一国产精品| 精品人妻熟女毛片av久久网站| 免费在线观看日本一区| 午夜福利免费观看在线| 又黄又粗又硬又大视频| 精品福利永久在线观看| 97精品久久久久久久久久精品| 久久精品成人免费网站| 久久久精品区二区三区| 精品国产乱码久久久久久小说| 91字幕亚洲| 亚洲国产av新网站| 国产成人系列免费观看| 国产亚洲午夜精品一区二区久久| 久久中文看片网| 淫妇啪啪啪对白视频 | 高清在线国产一区| 日本精品一区二区三区蜜桃| 99国产精品99久久久久| 国产精品av久久久久免费| 美国免费a级毛片| 亚洲精品中文字幕一二三四区 | 国产一区有黄有色的免费视频| 欧美一级毛片孕妇| 国产精品久久久人人做人人爽| 男女高潮啪啪啪动态图| 99re6热这里在线精品视频| 国产成人av激情在线播放| 在线精品无人区一区二区三| 亚洲中文日韩欧美视频| 天堂中文最新版在线下载| 国产精品熟女久久久久浪| 国产精品久久久久久人妻精品电影 | 妹子高潮喷水视频| 99香蕉大伊视频| www.自偷自拍.com| 亚洲第一欧美日韩一区二区三区 | 中文字幕人妻丝袜一区二区| 啦啦啦视频在线资源免费观看| 亚洲人成77777在线视频| 无限看片的www在线观看| 男女下面插进去视频免费观看| 一二三四在线观看免费中文在| 久久久久久久久久久久大奶| 亚洲精品在线美女| 国产精品秋霞免费鲁丝片| 91成人精品电影| 中文字幕色久视频| 在线观看免费视频网站a站| 成人影院久久| 十八禁网站网址无遮挡| 在线观看舔阴道视频| 亚洲七黄色美女视频| 丝袜脚勾引网站| 亚洲av国产av综合av卡| 中文字幕制服av| av片东京热男人的天堂| 亚洲国产看品久久| 黑人欧美特级aaaaaa片| 亚洲av电影在线进入| 免费在线观看视频国产中文字幕亚洲 | 亚洲av成人一区二区三| 天堂中文最新版在线下载| 国产精品一区二区精品视频观看| 动漫黄色视频在线观看| 波多野结衣一区麻豆| 韩国精品一区二区三区| 日韩中文字幕视频在线看片| 91麻豆精品激情在线观看国产 | 亚洲一码二码三码区别大吗| 秋霞在线观看毛片| 国产极品粉嫩免费观看在线| 免费黄频网站在线观看国产| 夫妻午夜视频| 捣出白浆h1v1| 咕卡用的链子| av网站免费在线观看视频| 性色av一级| 午夜老司机福利片| 国产淫语在线视频| 久久热在线av| 午夜福利免费观看在线| 亚洲午夜精品一区,二区,三区| 久久久久国内视频| 欧美+亚洲+日韩+国产| 成人国语在线视频| 1024香蕉在线观看| 精品高清国产在线一区| 永久免费av网站大全| 欧美成人午夜精品| 免费在线观看黄色视频的| 少妇精品久久久久久久| 91大片在线观看| 性高湖久久久久久久久免费观看| 五月天丁香电影| 丝袜人妻中文字幕| 午夜精品国产一区二区电影| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲成人国产一区在线观看| 午夜日韩欧美国产| 老鸭窝网址在线观看| 国产一级毛片在线| 狂野欧美激情性xxxx| 国产精品国产av在线观看| 一级毛片精品| 亚洲一码二码三码区别大吗| 欧美亚洲日本最大视频资源| 一本—道久久a久久精品蜜桃钙片| 在线观看免费视频网站a站| 黄色怎么调成土黄色| 香蕉丝袜av| 狂野欧美激情性bbbbbb| 18禁黄网站禁片午夜丰满| 美女午夜性视频免费| 成人18禁高潮啪啪吃奶动态图| 国产免费现黄频在线看| a级毛片黄视频| 91国产中文字幕| 韩国精品一区二区三区| svipshipincom国产片| www.999成人在线观看| 亚洲精品国产av成人精品| 又大又爽又粗| 日韩制服丝袜自拍偷拍| 免费久久久久久久精品成人欧美视频| 亚洲国产日韩一区二区| 欧美激情久久久久久爽电影 | 日日爽夜夜爽网站| 永久免费av网站大全| 国产成人欧美| 欧美一级毛片孕妇| 亚洲精品美女久久久久99蜜臀| 又黄又粗又硬又大视频| 欧美成狂野欧美在线观看| 亚洲成人免费av在线播放| 欧美大码av| 国产精品自产拍在线观看55亚洲 | 免费av中文字幕在线| 女人被躁到高潮嗷嗷叫费观| 另类亚洲欧美激情| 亚洲精品乱久久久久久| 亚洲少妇的诱惑av| 国产麻豆69| 国产精品麻豆人妻色哟哟久久| av天堂久久9| 99精国产麻豆久久婷婷| 侵犯人妻中文字幕一二三四区| 操出白浆在线播放| 国产淫语在线视频| 真人做人爱边吃奶动态| 亚洲激情五月婷婷啪啪| 少妇粗大呻吟视频| 国产精品免费大片| 大香蕉久久网| av免费在线观看网站| 三级毛片av免费| 国产免费av片在线观看野外av| 国产精品 欧美亚洲| 女性生殖器流出的白浆| 久热爱精品视频在线9| 国产极品粉嫩免费观看在线| 精品人妻在线不人妻| bbb黄色大片| 精品一品国产午夜福利视频| 中文字幕高清在线视频| 亚洲精品日韩在线中文字幕| av天堂在线播放| 搡老乐熟女国产| 黑丝袜美女国产一区| 亚洲三区欧美一区| 女人精品久久久久毛片| 亚洲久久久国产精品| 超碰成人久久| 两性夫妻黄色片| 亚洲第一av免费看| 一级毛片女人18水好多| 最近最新中文字幕大全免费视频| 狠狠精品人妻久久久久久综合| 欧美国产精品一级二级三级| 91成年电影在线观看| 亚洲成人免费电影在线观看| 中文字幕色久视频| 亚洲精品粉嫩美女一区| 一级a爱视频在线免费观看| 性少妇av在线| 欧美老熟妇乱子伦牲交| svipshipincom国产片| 99精国产麻豆久久婷婷| 国产无遮挡羞羞视频在线观看| 777米奇影视久久| 五月开心婷婷网| 黑人巨大精品欧美一区二区mp4| 国产精品秋霞免费鲁丝片| 亚洲成人国产一区在线观看| 欧美国产精品va在线观看不卡| 精品熟女少妇八av免费久了| 国产精品久久久人人做人人爽| 精品久久久精品久久久| 亚洲熟女精品中文字幕| 一区二区日韩欧美中文字幕| 国产成人精品久久二区二区免费| 免费人妻精品一区二区三区视频| tube8黄色片| 亚洲精品av麻豆狂野| 美女大奶头黄色视频| 嫩草影视91久久| 久久热在线av| cao死你这个sao货| 精品人妻1区二区| 高潮久久久久久久久久久不卡| 亚洲国产精品一区二区三区在线| 欧美日韩亚洲高清精品| av一本久久久久| 久久人妻熟女aⅴ| 亚洲全国av大片| 久久影院123| 99国产精品一区二区蜜桃av | 午夜福利影视在线免费观看| 欧美精品一区二区免费开放| 巨乳人妻的诱惑在线观看| 91精品伊人久久大香线蕉| 免费看十八禁软件| 欧美精品高潮呻吟av久久| 日日夜夜操网爽| 这个男人来自地球电影免费观看| 国产在线观看jvid| 精品高清国产在线一区| 国产精品成人在线| 午夜福利乱码中文字幕| 一本久久精品| 操出白浆在线播放| 精品乱码久久久久久99久播| 手机成人av网站| 国产97色在线日韩免费| 亚洲精品国产精品久久久不卡|