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

    耦合界面張力的三維流體界面不穩(wěn)定性的格子Boltzmann 模擬*

    2022-03-04 02:10:04馬聰劉斌梁宏
    物理學(xué)報(bào) 2022年4期
    關(guān)鍵詞:不穩(wěn)定性雷諾數(shù)表面張力

    馬聰 劉斌 梁宏

    (杭州電子科技大學(xué)理學(xué)院,杭州 310018)

    采用介觀格子Boltzmann 方法模擬界面張力作用下三維流體界面的Rayleigh-Taylor (RT)不穩(wěn)定性的增長過程,主要分析表面張力對(duì)流體界面動(dòng)力學(xué)行為及尖釘和氣泡后期增長的影響機(jī)制.首先發(fā)現(xiàn)三維RT 不穩(wěn)定性的發(fā)生存在臨界表面張力(σc),其值隨著流體Atwood 數(shù)的增大而增大,且數(shù)值預(yù)測(cè)值與理論分析結(jié)果 σc (ρh-ρl)g/k2 一致.另外,隨著表面張力的增大,不穩(wěn)定性演化過程中界面卷吸程度和結(jié)構(gòu)復(fù)雜性逐漸減弱,系統(tǒng)中界面破裂形成離散液滴的數(shù)目也顯著減少.相界面的后期動(dòng)力學(xué)行為也從非對(duì)稱發(fā)展轉(zhuǎn)向始終保持關(guān)于中軸線對(duì)稱.尖釘與氣泡振幅在表面張力較小時(shí)對(duì)其變化不顯著,當(dāng)表面張力增大到一定值后,可以有效地抑制尖釘與氣泡振幅的增長.進(jìn)一步發(fā)現(xiàn),高雷諾數(shù)三維RT 不穩(wěn)定性在不同表面張力下均經(jīng)歷4 個(gè)不同的發(fā)展階段:線性階段、飽和速度階段、重加速和混沌混合階段.尖釘與氣泡在飽和速度階段以近似恒定的速度增長,其漸進(jìn)速度的值與修正的勢(shì)流理論模型結(jié)果一致.受非線性Kelvin-Helmholtz 旋渦的剪切作用,尖釘與氣泡隨后的增長被加速,導(dǎo)致在重加速階段的演化速度超過勢(shì)流模型的解析解.重加速階段不能持續(xù)發(fā)展下去,尖釘與氣泡在不穩(wěn)定性后期的增長速度會(huì)隨時(shí)間上下波動(dòng),這表明不穩(wěn)定性的演化進(jìn)入了混沌混合階段.通過數(shù)值分析,證實(shí)了三維RT 不穩(wěn)定性在后期的混沌混合階段具有二次增長的規(guī)律,并且尖釘與氣泡增長率總體上隨著表面張力的增大而逐漸減少.

    1 引言

    Rayleigh-Taylor (RT)不穩(wěn)定性是低密度流體加速高密度流體時(shí)兩者密度交界面的不穩(wěn)定現(xiàn)象,廣泛存在于自然現(xiàn)象(如卷云的形成、超新星爆發(fā)、地下鹽丘和火山島的形成)和工程應(yīng)用(如超燃沖壓發(fā)動(dòng)機(jī)、慣性約束核聚變、氣象學(xué)、海洋運(yùn)動(dòng)學(xué))中,也是多相流體界面不穩(wěn)定性、湍流混合等基礎(chǔ)問題的理論基礎(chǔ)[1,2].因此,研究RT 不穩(wěn)定性的發(fā)展與演化規(guī)律對(duì)于理論研究和工程實(shí)踐都具有重要的科學(xué)價(jià)值和現(xiàn)實(shí)意義.著名學(xué)者Rayleigh[3]和Taylor[4]最早描述了RT 不穩(wěn)定性現(xiàn)象,提出了著名的線性穩(wěn)定性理論:交界面處初始擾動(dòng)隨時(shí)間以指數(shù)形式增長.后來,Lewis[5]實(shí)驗(yàn)證實(shí)了線性穩(wěn)定性理論可以有效地描述擾動(dòng)振幅的增長達(dá)到 0 .4λ(λ是初始擾動(dòng)的波長),并發(fā)現(xiàn)不穩(wěn)定性隨后進(jìn)入非線性增長階段.2001 年,Waddell等[6]實(shí)驗(yàn)研究了二維單模RT 不穩(wěn)定性問題,觀察到緊接著線性階段,尖釘與氣泡振幅在非線性階段的平均增長速度接近于常數(shù).Goncharov[7]理論分析單模RT 不穩(wěn)定性的非線性增長,并給出尖釘與氣泡在非線性階段以恒定速度增長的勢(shì)流模型:

    其中,ub是氣泡的速度;us是尖釘?shù)乃俣?g是 重力加速度;At是Atwood 數(shù);k是波數(shù);Cg對(duì)于三維流動(dòng)取 1,而對(duì)于二維情形則取3.受Goncharov工作[7]的啟發(fā),Sohn[8]進(jìn)一步分析了流體黏性和表面張力對(duì)氣泡漸進(jìn)速度的影響,提出了包含這些因素的修正勢(shì)流模型,三維情形下可以表述為

    其中ν是流體黏性,ρh表示重流體的密度.Glimm等[9]基于前追蹤方法模擬了二維單模RT 不穩(wěn)定性問題,首次發(fā)現(xiàn)氣泡與尖釘?shù)难莼俣冉?jīng)過恒定速度增長后會(huì)出現(xiàn)再加速行為.RT 不穩(wěn)定性的再加速現(xiàn)象在后續(xù)的實(shí)驗(yàn)[10]和數(shù)值模擬[11]中被進(jìn)一步證實(shí),并被歸因于氣泡與尖釘界面處形成的Kelvin-Helmholtz 渦旋的作用.Bian 等[12]采用直接數(shù)值模擬方法研究了不同雷諾數(shù)和阿特伍德數(shù)下渦量對(duì)三維單模RT 不穩(wěn)定性的再加速增長階段的影響,并確定了再加速過程與渦量之間有很強(qiáng)的相關(guān)性.2012 年,Ramaprabhu 等[13]數(shù)值研究了混相流體的三維單模RT 不穩(wěn)定性的后期增長規(guī)律,首次發(fā)現(xiàn)不穩(wěn)定性在高雷諾數(shù)時(shí)經(jīng)歷4 個(gè)不同的發(fā)展階段,包括線性增長、飽和速度增長、再加速和混沌混合階段.同期,Wei 和Livescu[14]基于直接數(shù)值模擬方法高精度模擬了低阿特伍德數(shù)流體的二維單模RT 不穩(wěn)定性的后期增長過程,同樣觀察到高雷諾數(shù)的單模RT 不穩(wěn)定性經(jīng)歷上述4 個(gè)發(fā)展階段,并且報(bào)道氣泡振幅在后期的混沌混合階段具有平均二次增長的規(guī)律.Hu 等[15]模擬了不同雷諾數(shù)下的二維單模RT 不穩(wěn)定性問題,發(fā)現(xiàn)在中等雷諾數(shù)條件下,氣泡與尖釘在后期階段會(huì)出現(xiàn)反復(fù)的加速與減速現(xiàn)象.李德梅等[16]近期利用離散Boltzmann 方法研究了可壓縮流體的二維單模RT 不穩(wěn)定性的非平衡效應(yīng).另外,Liang 等[17-19]針對(duì)非混相流體的單模RT 不穩(wěn)定性后期演化過程也開展了一系列介尺度數(shù)值研究,分析了廣泛的雷諾數(shù)和阿特伍德數(shù)對(duì)不穩(wěn)定性的后期增長階段和相界面動(dòng)力學(xué)行為的影響.

    上述研究豐富了人們對(duì)單模RT 不穩(wěn)定性演化機(jī)制的認(rèn)識(shí),但均未考慮兩相流體間界面張力對(duì)不穩(wěn)定性增長的影響.已有研究表明,表面張力會(huì)顯著影響多相流體界面輸運(yùn)現(xiàn)象及流體界面不穩(wěn)定性行為,特別是RT 不穩(wěn)定性在界面張力作用下可以顯示出毛細(xì)波、收縮、破裂等獨(dú)特的界面動(dòng)力學(xué)行為.鑒于此,一些學(xué)者嘗試探索表面張力對(duì)單模RT 不穩(wěn)定性增長的影響,如Daly[20]采用有限差分方法模擬了界面張力對(duì)單模RT 不穩(wěn)定性早期增長的影響,發(fā)現(xiàn)增大表面張力可以減小線性增長率;Zhang 等[21]利用多相流格子Boltzmann 方法模擬了表面張力作用下二維單模RT 不穩(wěn)定性的演化過程,發(fā)現(xiàn)存在表面張力作用可以抑制氣泡與尖釘?shù)脑鲩L;Young 和Ham[22]模擬了表面張力作用下單模與多模RT 不穩(wěn)定性的演化機(jī)制,發(fā)現(xiàn)增大界面張力可以抑制混合層的增長;Matsuoka[23]采用邊界積分法模擬了耦合表面張力的可壓縮流體的二維單模RT 不穩(wěn)定性的演化過程,發(fā)現(xiàn)存在表面張力有利于誘導(dǎo)界面發(fā)生夾斷行為以及抑制不穩(wěn)定性湍流現(xiàn)象的發(fā)生;Sohn[24]與夏同軍等[25]分析了表面張力對(duì)單模RT 不穩(wěn)定性的飽和速度增長階段的影響,給出了包含表面張力效應(yīng)的氣泡漸進(jìn)速度的解析表達(dá)式;黃皓偉等[26]數(shù)值研究了表面張力對(duì)高雷諾數(shù)的二維單模RT 不穩(wěn)定性后期增長的影響,發(fā)現(xiàn)增大表面張力可以減弱演化過程中相界面結(jié)構(gòu)的復(fù)雜程度以及有效抑制后期階段相界面破裂行為,并進(jìn)一步給出了不同表面張力下氣泡與尖釘?shù)暮笃谠鲩L率.

    綜上所述,對(duì)耦合界面張力的單模RT 不穩(wěn)定性的研究已取得了一些進(jìn)展,但相關(guān)研究還不夠深入.絕大多工作[20-25]局限于不穩(wěn)定性發(fā)展的中前期階段,所考慮的雷諾數(shù)較小,且所研究的物理工況均為二維情形[20-26],而對(duì)表面張力作用下三維單模RT 不穩(wěn)定性的演化后期及氣泡與尖釘增長的描述尚未報(bào)道.本文將基于介觀格子Boltzmann方法系統(tǒng)研究高雷諾數(shù)的三維單模RT 不穩(wěn)定性的后期演化機(jī)制,主要分析界面張力對(duì)相界面動(dòng)力學(xué)行為及氣泡與尖釘后期增長的影響.

    2 數(shù)值模型與問題描述

    格子Boltzmann 方法是基于氣體動(dòng)理學(xué)理論的介觀數(shù)值方法,具有清晰的物理背景和粒子演化特性,并行計(jì)算效率高,可以從底層方便地描述流體內(nèi)部與環(huán)境間的相互作用,因而在模擬多相多組分等復(fù)雜流體輸運(yùn)問題上具有很大的優(yōu)勢(shì)[27].根據(jù)流體間作用力的不同物理背景,現(xiàn)有的多相流格子Boltzmann 模型主要可以分為4 類:顏色模型、偽勢(shì)模型、自由能模型和相場模型.不同于前3 類模型,相場格子Boltzmann 模型需要求解明確的界面追蹤方程,相界面求解精度高,適用于模擬界面大拓?fù)渥兓亩嘞嗔鲉栴}[28,29].因此,本文采用相場格子Boltzmann 模型研究三維非混相Rayleigh-Taylor 不穩(wěn)定性的后期增長規(guī)律.該模型利用兩套粒子分布函數(shù)fi和gi,并采用多松弛碰撞算子以提高模型處理高雷諾數(shù)流動(dòng)的數(shù)值穩(wěn)定性.雙分布函數(shù)的多松弛格子Boltzmann 方法的演化方程可以表示為如下的統(tǒng)一形式[30]:

    其中,fi(x,t) 為序參數(shù)的粒子分布函數(shù),gi(x,t) 為刻畫流場的粒子分布函數(shù),(x,t) 和(x,t) 則為粒子分布函數(shù)所對(duì)應(yīng)的平衡態(tài)分布函數(shù),Mf和Mg表示碰撞矩陣,Sf和Sg為松弛矩陣,F(xiàn)i(x,t)和Gi(x,t) 分別表示源項(xiàng)和外力分布函數(shù).為了恢復(fù)正確的相場方程和不可壓流體動(dòng)力學(xué)方程,平衡態(tài)分布函數(shù)(x,t) 和(x,t) 分別構(gòu)造為[30]

    其中,η是與遷移率相關(guān)的調(diào)節(jié)參數(shù);cs為聲速;ci和ωi分別為速度空間的離散速度和權(quán)系數(shù),取值依賴于選取的格子模型.針對(duì)三維的Cahn-Hilliard 方程,格子Boltzmann 方法至少需要7 個(gè)離散速度才能滿足恢復(fù)宏觀控制方程所約束的矩條件[31],因此,為了提高計(jì)算效率,在序參數(shù)分布函數(shù)的演化方程(3a)中采用高效的D3Q7 格子模型,其對(duì)應(yīng)的權(quán)系數(shù)ωi為ω01/4,ω1-61/8,c2/4,離散速度ci定義為

    另外,針對(duì)三維流體動(dòng)力學(xué)方程,在流場分布函數(shù)的演化方程(3b)中采用通用的D3Q15 格子模型,其對(duì)應(yīng)的權(quán)系數(shù)為ω02/9,ω1-61/9,ω7-141/72,c2/3,且離散速度ci設(shè)定為

    三維多松弛格子Boltzmann 模型的其他元素中碰撞矩陣Mf和Mg的選取可參考文獻(xiàn)[31,32].為了恢復(fù)正確的相界面追蹤的Cahn-Hilliard 方程,Liang 等[30]在演化方程(3a)中引入關(guān)于時(shí)間導(dǎo)數(shù)的局部源項(xiàng)以提高相界面求解的精度:

    另一方面,本文考慮了外力引入格子Boltzmann方法所產(chǎn)生的離散效應(yīng),并采用了 He 等[33]提出的完全消除離散誤差的外力格式:

    其中,G為流體系統(tǒng)所受的外力;Fa是為了恢復(fù)正確的動(dòng)量方程而引入的額外界面力,其定義為表示兩相流體間的表面張力;μ為相場模型中的化學(xué)勢(shì),

    其中β,k為與界面張力(σ)和界面厚度(D)相關(guān)的模型參數(shù),.通過計(jì)算粒子分布函數(shù)的零階矩和一階矩,可以得到宏觀統(tǒng)計(jì)量[30]

    另外,流體密度ρ可由序參數(shù)φ的線性插值計(jì)算,

    其中,ρl為輕流體的密度.在實(shí)際模擬中,需要采用合適的差分格式離散格子Boltzmann 模型的導(dǎo)數(shù)項(xiàng),本文采用顯式的歐拉格式計(jì)算時(shí)間導(dǎo)數(shù)項(xiàng):

    同時(shí)利用各向同性的二階差分格式計(jì)算空間梯度和拉普拉斯算子:

    其中χ表示任意一個(gè)變量.

    為了研究流體不穩(wěn)定性的后期演化規(guī)律,本文考慮的物理工況為z軸方向足夠長的長方體通道.如圖 1 所示,Lx,Ly,Lz分別表示管道的長度、寬度與高度,且Lx×Ly×LzW×W×16W.初始時(shí)刻,重流體位于管道的上方而輕流體置于管道的下方,并在兩相流體界面z8W處施加一個(gè)微小的擾動(dòng):

    圖1 三維單模RT 不穩(wěn)定性問題的示意圖Fig.1.Schematic of three-dimensional single-mode RT instability problem.

    根據(jù)相場理論,序參量的初始分布設(shè)定為如下的雙曲正切函數(shù):

    其值在體相區(qū)取0 和1,而在界面區(qū)域從0 連續(xù)變化至1.為了表征RT 不穩(wěn)定性的演化特性,需要引入兩個(gè)重要的無量綱參數(shù),即雷諾數(shù)(Reynolds number,Re)和阿特伍德數(shù)(Atwood number,At),分別定義為[33]

    其中,ν表示流體黏性,重流體的密度ρh給定為 1,而輕流體的密度ρl可以根據(jù)阿特伍德數(shù)的值來確定.為了描述重力效應(yīng),對(duì)管道中輕重流體均施加一個(gè)豎直方向上的浮力:

    本文著重分析表面張力對(duì)高雷諾數(shù)RT 不穩(wěn)定性的后期演化特性的影響,除特別聲明,流動(dòng)的雷諾數(shù)和Atwood數(shù)分別固定為Re5000,At0.1,其他物理參數(shù)給定為W100,D4,松弛矩陣Sf和Sg分別設(shè)定為

    3 數(shù)值結(jié)果和討論

    通過數(shù)值計(jì)算發(fā)現(xiàn),當(dāng)表面張力較大時(shí),RT 不穩(wěn)定性的初始擾動(dòng)會(huì)發(fā)生振蕩直至相界面變?yōu)榉€(wěn)定的水平平面,即RT 不穩(wěn)定性現(xiàn)象將不會(huì)發(fā)生.這表明RT 不穩(wěn)定性的發(fā)生存在一個(gè)臨界表面張力,當(dāng)表面張力大于該臨界值時(shí),RT 不穩(wěn)定性行為將不會(huì)出現(xiàn),而當(dāng)表面張力小于該臨界值時(shí),RT 不穩(wěn)定性現(xiàn)象將會(huì)發(fā)生.本文統(tǒng)計(jì)了不同流體Atwood 數(shù)下三維RT 不穩(wěn)定性發(fā)生的臨界表面張力,結(jié)果如圖 2 所示.可以發(fā)現(xiàn),臨界表面張力隨著Atwood 數(shù)的增大而呈現(xiàn)遞增規(guī)律.進(jìn)一步理論分析了RT 不穩(wěn)定性現(xiàn)象發(fā)生的臨界表面張力.根據(jù)線性穩(wěn)定性理論[6,10],RT 不穩(wěn)定性的擾動(dòng)振幅在初始階段的增長具有指數(shù)形式,即ha1eγt+a2e-γt,其中h表示擾動(dòng)振幅,t是演化時(shí)間,a1和a2為擬合系數(shù),γ是線性增長因子且可表示為[6]

    要使RT 不穩(wěn)定性能夠發(fā)生,線性增長因子中被開方數(shù)的值必須大于 0,從而可以推導(dǎo)出臨界表面張力的解析表達(dá)式為

    其中波數(shù)k2π/W.圖2 進(jìn)一步給出了不同Atwood 數(shù)下臨界表面張力的理論值,可以發(fā)現(xiàn)理論值基本落在數(shù)值模擬所預(yù)測(cè)的范圍內(nèi),即數(shù)值模擬結(jié)果與理論分析結(jié)果一致.在接下來的研究中,只考慮RT 不穩(wěn)定性能夠發(fā)生的情形,即σ<σc.

    圖2 不同Atwood 數(shù)下三維RT 不穩(wěn)定性的臨界表面張力Fig.2.Critical surface tensions of three-dimensional RT instability at various Atwood numbers.

    圖3 給出了不同界面張力作用下三維單模RT 不穩(wěn)定性的相界面的演化過程.可以看出,三維RT 不穩(wěn)定性在演化初期顯示出相似的界面動(dòng)力學(xué)特征,重流體向下運(yùn)動(dòng)而輕流體向上升起,即輕流體和重流體之間相互滲透,且滲透深度隨著時(shí)間的演化而逐漸增加,從而輕流體上升形成了氣泡,而重流體下落形成了尖釘.緊接著,尖釘繼續(xù)往下運(yùn)動(dòng),氣泡繼續(xù)上升.與此同時(shí),兩流體剪切作用產(chǎn)生的Kelvin-Helmholtz 不穩(wěn)定性開始發(fā)展,其強(qiáng)度也隨時(shí)間逐漸增強(qiáng),并開始影響相界面的動(dòng)力學(xué)行為.在Kelvin-Helmholtz 不穩(wěn)定性的作用下,尖釘向上卷起,形成了經(jīng)典的蘑菇結(jié)構(gòu)(見圖3(a)—(c)中t8 時(shí)刻與圖3(d)的t10 時(shí)刻).同時(shí)也注意到尖釘?shù)木砥鸱仍诒砻鎻埩^小時(shí)差異不大,但繼續(xù)增大表面張力會(huì)減弱尖釘卷起程度,并且卷起發(fā)生的時(shí)刻也會(huì)隨之推遲.接下來觀察到三維單模RT 不穩(wěn)定性在不同的表面張力作用下,表現(xiàn)出的顯著不同的界面動(dòng)力學(xué)行為.當(dāng)表面張力較小(σ10-6或 1 0-5)時(shí),尖釘隨著時(shí)間繼續(xù)往下運(yùn)動(dòng),在流體間剪切力的作用下,卷起表面變得不光滑,出現(xiàn)了凹凸不平的結(jié)構(gòu),并且在卷起尾端形成了4 個(gè)尺寸相同的類似于“卷發(fā)”的界面結(jié)構(gòu)(見t10 時(shí)刻).隨后,不穩(wěn)定性系統(tǒng)的非線性效應(yīng)越來越劇烈,在高強(qiáng)度的Kelvin-Helmholtz 不穩(wěn)定性的剪切作用下,產(chǎn)生了許多不同尺度的旋渦,進(jìn)而使得流體界面變得異常不穩(wěn)定,演化后期相界面發(fā)生了多層次卷起、劇烈變形、混沌破裂,形成了非常復(fù)雜的拓?fù)浣Y(jié)構(gòu),并且使得流體系統(tǒng)中產(chǎn)生了大量的游離的離散液滴.另外可以發(fā)現(xiàn),在低界面張力的三維RT 不穩(wěn)定性的增長后期,相界面的演化圖案關(guān)于中心軸失去了對(duì)稱性.相界面的非對(duì)稱性發(fā)展在高雷諾數(shù)下混相流體的三維單模RT 不穩(wěn)定性現(xiàn)象[13]中同樣被觀察到.當(dāng)表面張力增至σ2×10-4,流體間的剪切作用減弱,流體界面卷起表面相對(duì)比較光滑,未出現(xiàn)凹凸不平的結(jié)構(gòu),形成的“卷發(fā)”長度也隨之減小.此后,非線性Kelvin-Helmholtz 不穩(wěn)定性的強(qiáng)度隨時(shí)間逐漸增強(qiáng),致使相界面在多個(gè)位置發(fā)生卷起甚至出現(xiàn)破裂行為,最終在演化后期也形成了非常復(fù)雜的拓?fù)浣Y(jié)構(gòu).但不同于上述表面張力較低時(shí)的情形,高雷諾數(shù)三維單模RT 不穩(wěn)定性在較大表面張力的作用下,相界面演化圖案始終保持著關(guān)于中軸線對(duì)稱的特性,并且整個(gè)演化過程中界面破裂的現(xiàn)象減少,系統(tǒng)中形成離散液滴的數(shù)量也相應(yīng)地減少.特別地,將表面張力增大至σ5×10-4,在卷起尾端未觀察到明顯的“卷發(fā)”結(jié)構(gòu),雖然在演化后期仍然觀察到界面的多層次卷起行為,但界面卷起程度和界面結(jié)構(gòu)的復(fù)雜性顯著減弱,系統(tǒng)中生成的離散液滴數(shù)目也顯著減少.另外,在整個(gè)不穩(wěn)定性的演化過程中,相界面始終維持關(guān)于中軸線的對(duì)稱性.

    圖3 表面張力對(duì)高雷諾數(shù)三維單模RT 不穩(wěn)定性相界面演化的影響,R e=5000,A t=0.1 (a) σ=1×10-6 ;(b)σ=1×10-5 ;(c) σ=2×10-4 ;(d)σ=5×10-4Fig.3.Effect of the surface tension on the evolution of phase interface in three-dimensional single-mode RT instability with high Reynolds number,R e=5000,A t=0.1 :(a) σ=1×10-6 ;(b) σ=1×10-5 ;(c) σ=2×10-4 ;(d) σ=5×10-4 .

    為了更細(xì)致地顯示表面張力對(duì)界面動(dòng)力學(xué)行為的影響,給出了上述界面張力作用下三維RT 不穩(wěn)定性中對(duì)角平面 (xy) 的界面演化圖案,結(jié)果如圖4 所示.可以看出,在最初階段,流體界面在不同界面張力作用下均表現(xiàn)出相似的行為,重流體在重力作用下向下運(yùn)動(dòng)形成尖釘,而密度較小的流體向上升起形成了氣泡,并受Kelvin-Helmholtz不穩(wěn)定性的影響,界面在兩層位置發(fā)生卷吸現(xiàn)象,從而形成了兩對(duì)旋轉(zhuǎn)方向相反的旋渦.界面卷起的程度會(huì)隨著表面張力的增大而逐漸減弱,卷起發(fā)生的時(shí)間也會(huì)隨之而推遲.界面兩層卷起行為是三維單模RT 不穩(wěn)定性獨(dú)有的特征,在二維單模RT 不穩(wěn)定性的演化過程中并未出現(xiàn).接下來,對(duì)于低表面張力(σ10-6或 1 0-5)情形,形成的兩對(duì)旋渦隨著時(shí)間演化而逐漸增大,導(dǎo)致在各自卷起尾端處分別形成了一對(duì)二級(jí)旋渦.在多個(gè)旋渦相互作用下,界面在多個(gè)位置發(fā)生卷起行為,形成了一些不同尺度的旋渦結(jié)構(gòu).最終,在高流體界面剪切作用下,中軸線上界面的多處位置發(fā)生卷吸、變形、混沌的破裂,形成了許多離散的小液滴.另外,進(jìn)一步可以發(fā)現(xiàn),對(duì)于低表面張力情形,流體界面在演化后期難以維持關(guān)于中軸線的對(duì)稱性.當(dāng)表面張力較大(σ2×10-4)時(shí),兩對(duì)一級(jí)旋渦隨時(shí)間繼續(xù)增長,伴隨著卷起部分的長度越來越長,并逐漸開始收縮,與中軸線附近的流體界面發(fā)生接觸,相比低表面張力情形,在界面卷起的尾端未出現(xiàn)二級(jí)旋渦的界面結(jié)構(gòu).隨后,非線性Kelvin-Helmholtz 不穩(wěn)定性的影響逐漸加強(qiáng),導(dǎo)致界面多個(gè)位置發(fā)生卷起變形,形成較復(fù)雜的結(jié)構(gòu),但相界面在整個(gè)演化過程中始終保持中軸線的對(duì)稱性.當(dāng)界面張力繼續(xù)增大時(shí),流體界面變得相對(duì)簡單與光滑,界面卷起的長度和程度也隨之減小,未觀察到劇烈拓?fù)渥兓膹?fù)雜界面結(jié)構(gòu).

    圖4 表面張力對(duì)高雷諾數(shù)三維單模RT 不穩(wěn)定性的對(duì)角平面(x=y 平面)的相界面演化的影響,R e=5000,A t=0.1 (a)σ=1×10-6 ;(b) σ=1×10-5 ;(c) σ=2×10-4 ;(d)σ=5×10-4Fig.4.Effect of the surface tension on the evolution of phase interface in the diagonal plane of three-dimensional single-mode RT instability with high Reynolds number,R e=5000,A t=0.1 :(a) σ=1×10-6 ;(b) σ=1×10-5 ;(c) σ=2×10-4 ;(d)σ=5×10-4.

    上述研究分析了表面張力對(duì)三維RT 不穩(wěn)定性的相界面動(dòng)力學(xué)的影響,而尖釘與氣泡振幅及增長速度是描述RT 不穩(wěn)定性演化特征的幾個(gè)重要物理量.為了進(jìn)一步表征界面張力效應(yīng),定量統(tǒng)計(jì)了廣泛表面張力條件下尖釘與氣泡振幅、增長速度隨時(shí)間的演化規(guī)律.尖釘與氣泡的振幅定義為尖釘與氣泡的實(shí)時(shí)位置與初始位置在z方向的間距,因此在初始時(shí)刻,尖釘與氣泡振幅均為零.圖 5 給出了尖釘與氣泡振幅在不同表面張力下隨時(shí)間變化的演化曲線.可以看出,尖釘與氣泡振幅在不同表面張力下隨時(shí)間演化不斷增長,并在較小表面張力時(shí),同一時(shí)刻可以獲得更大的尖釘與氣泡的振幅.具體而言,當(dāng)表面張力在較小范圍內(nèi),尖釘與氣泡振幅的演化曲線幾乎相互重合,即界面張力對(duì)不穩(wěn)定性增長的影響不再顯著.在非混相不可壓流體的RT 不穩(wěn)定性的演化中,尖釘與氣泡的動(dòng)力學(xué)特征理論上由單位質(zhì)量的重力、黏性耗散力、界面張力之間的競爭決定[34].針對(duì)表面張力極小的情形,表面張力遠(yuǎn)遠(yuǎn)小于系統(tǒng)中重力與耗散力,因此對(duì)不穩(wěn)定性中尖釘與氣泡發(fā)展的影響可以忽略.當(dāng)表面張力增大到一定值后,它可以有效地抑制尖釘與氣泡振幅的增長,即尖釘與氣泡振幅隨著表面張力的增長而顯著減小.進(jìn)一步模擬了不同Atwood 數(shù)下界面張力對(duì)三維單模RT 不穩(wěn)定性增長的影響,同樣發(fā)現(xiàn)表面張力較小時(shí)對(duì)三維不穩(wěn)定性中尖釘與氣泡的增長影響不大,繼續(xù)增大表面張力可以有效抑制不穩(wěn)定性的發(fā)展.圖 6 描述了不同界面張力下尖釘與氣泡隨時(shí)間演化的增長速度.根據(jù)速度演化曲線,發(fā)現(xiàn)高雷諾數(shù)的三維單模RT 不穩(wěn)定性在不同界面張力下的增長仍然可以劃分為4 個(gè)階段:線性增長階段、飽和速度階段、重加速階段和混沌混合階段.當(dāng)表面張力充分小時(shí),氣泡和尖釘?shù)脑鲩L速度曲線在線性階段相互重合,繼續(xù)增大表面張力可以抑制尖釘與氣泡的線性增長,這符合線性穩(wěn)定性理論的分析結(jié)果:線性增長因子如(23)式所示,在表面張力較小時(shí)差別不大,隨著表面張力的增大而逐漸減小.緊接著,尖釘與氣泡均以近似恒定的速度增長,這表明三維RT 不穩(wěn)定性的演化進(jìn)入了飽和速度增長階段.經(jīng)典的勢(shì)流理論模型[7]預(yù)測(cè)了理想流體的單模RT 不穩(wěn)定性中尖釘與氣泡在飽和速度階段的漸進(jìn)速度.隨后,Sohn[24]基于上述勢(shì)流模型分析了流體黏性和界面張力對(duì)氣泡飽和速度增長的影響,給出了氣泡漸進(jìn)速度的理論解,如(2)式所示.受Goncharov 勢(shì)能模型[7]的啟發(fā),尖釘在飽和速度增長階段的漸進(jìn)速度可以表示為

    圖5 表面張力對(duì)隨時(shí)間演化的尖釘(左)和氣泡(右)振幅的影響Fig.5.Influence of surface tension on the time evolutions of spike (left) and bubble (right) amplitudes.

    圖6 表面張力對(duì)隨時(shí)間演化的尖釘(左)和氣泡(右)增長速度的影響.黑色和藍(lán)色虛線分別表示修正的勢(shì)能模型所預(yù)測(cè)尖釘與氣泡漸進(jìn)速度在 σ=1×10-5 σ=5×10-4 時(shí)的解析解Fig.6.Influence of surface tension on the time evolutions of spike (left) and bubble (right) growth velocities.The black and blue dotted lines represent the analytical solutions of the spike and bubble asymptotic velocities from the modified potential flow model at σ=1×10-5 and σ=5×10-4 .

    根據(jù)Sohn[24]改進(jìn)的勢(shì)流模型,可以發(fā)現(xiàn)氣泡與尖釘在界面張力σ<1×10-5時(shí)飽和速度的解析解差別極小,因此在圖6 中僅給出σ1×10-5和σ5×10-4時(shí)尖釘與氣泡漸進(jìn)速度的理論值.可以發(fā)現(xiàn),本文通過格子Boltzmann 方法所預(yù)測(cè)的尖釘與氣泡漸進(jìn)速度與勢(shì)流理論模型的結(jié)果一致,進(jìn)一步可以發(fā)現(xiàn),尖釘與氣泡在表面張力較小時(shí)更早地進(jìn)入飽和速度階段,并且該階段的持續(xù)時(shí)間隨著表面張力的增大而增長.接下來,受系統(tǒng)中非線性Kelvin-Helmholtz 旋渦的剪切作用,使得尖釘和氣泡的增長速度超過勢(shì)流模型的理論值,這預(yù)示著不穩(wěn)定性的演化進(jìn)入了重加速階段.重加速階段是由遠(yuǎn)離初始中心平面的輕流體與尖釘?shù)慕缑嫣幃a(chǎn)生的渦旋向氣泡尖端傳播,導(dǎo)致氣泡重新加速,其發(fā)生必須滿足兩個(gè)條件:(a)渦旋需要在豎直方向上比氣泡更快地移動(dòng);(b)渦旋的結(jié)構(gòu)必須要保持足夠長時(shí)間.結(jié)合相界面動(dòng)力學(xué)圖案,可以發(fā)現(xiàn)對(duì)所有表面張力情形,高雷諾數(shù)三維RT 不穩(wěn)定性均可以產(chǎn)生保留時(shí)間足夠長的旋渦,導(dǎo)致出現(xiàn)重加速階段,但由于旋渦的強(qiáng)度隨著界面張力的增大而減小,表面張力較小的RT 不穩(wěn)定性可以更早地進(jìn)入重加速階段.重加速階段不能持續(xù)地發(fā)展下去,從圖 6 可以發(fā)現(xiàn),對(duì)所有表面張力情形,尖釘與氣泡在演化后期的增長速度隨時(shí)間上下波動(dòng),其值反復(fù)地加速與減速,這表明RT 不穩(wěn)定性的發(fā)展最終進(jìn)入了混沌混合階段.已有的高精度直接數(shù)值模擬研究顯示單模RT 不穩(wěn)定性在高雷諾數(shù)條件下會(huì)出現(xiàn)后期的混沌混合增長階段,在該階段尖釘與氣泡的振幅呈現(xiàn)平均二次增長的規(guī)律,即hs,bαs,bgAtt2,其中αs,b表示尖釘與氣泡的后期增長率.αs,b反映了不穩(wěn)定性在混沌混合階段的發(fā)展規(guī)律,因此是RT 不穩(wěn)定性后期研究中最為重要的物理統(tǒng)計(jì)量.為了計(jì)算不穩(wěn)定性湍流混合階段的增長率,一些學(xué)者基于振幅與演化時(shí)間的二次關(guān)系提出多種不同的測(cè)量方法,Liang 等[19]近期在二維單模RT 不穩(wěn)定性的后期研究中詳細(xì)地比較了這些測(cè)量方法,指出 Olson 等[35]提出的計(jì)算增長率的方法具有良好的收斂性與一致性,并被用于本文測(cè)量不同表面張力下三維單模RT 不穩(wěn)定性的后期增長率.Olson 等[35]對(duì)演化后期的振幅與時(shí)間的二次關(guān)系式兩邊開方得到,因此通過與演化時(shí)間t的線性擬合可以獲得擬合直線的斜率,進(jìn)而通過計(jì)算可以得到尖釘與氣泡的后期增長率αs,b.圖 7 給出了不同表面張力條件下尖釘與氣泡振幅的開方與的關(guān)系曲線以及在后期相應(yīng)的線性擬合直線.從圖 7 可以發(fā)現(xiàn),兩者在不穩(wěn)定性的演化后期確實(shí)成線性關(guān)系,這也驗(yàn)證了三維單模RT 不穩(wěn)定性在后期混沌混合階段具有二次增長的規(guī)律.進(jìn)一步通過線性擬合,可以獲得尖釘增長率在表面張力σ1×10-7,1×10-6,1×10-5,2×10-4,5×10-4時(shí)分別為0.1293,0.1304,0.1348,0.1106,0.1069,而對(duì)應(yīng)的氣泡增長率依次為0.1286,0.1268,0.1275,0.0976,0.0691.可以發(fā)現(xiàn),對(duì)于相同的界面張力,尖釘?shù)暮笃谠鲩L率大于氣泡的增長率,即尖釘比氣泡運(yùn)動(dòng)得更快.另外,當(dāng)界面張力較小時(shí),界面張力對(duì)尖釘與氣泡增長率的影響不再顯著,而當(dāng)界面張力增加至一定值時(shí),尖釘與氣泡增長率則隨著界面張力的繼續(xù)增大而減小,即界面張力會(huì)抑制尖釘與氣泡的后期增長.

    圖7 不同界面張力下三維RT 不穩(wěn)定性中尖釘(左)與氣泡(右)振幅的開方 與演化時(shí)間 的關(guān)系曲線,實(shí)線表示在后期混沌混合階段的線性擬合結(jié)果Fig.7.Relationsbetweenthesquare roots of the spike and bubble amplitudes and the time in the three-dimension-al RT instability with different surface tensions.The solid lines represent the linear fitting results at the late-time chaotic mixing stage.

    4 結(jié)論

    本文采用多相流的介觀格子Boltzmann 方法模擬了三維非混相、不可壓流體的RT 不穩(wěn)定性問題,系統(tǒng)分析表面張力對(duì)不穩(wěn)定性中界面動(dòng)態(tài)過程和尖釘與氣泡定量增長的影響.通過數(shù)值模擬發(fā)現(xiàn),三維RT 不穩(wěn)定性的發(fā)生存在著一個(gè)臨界表面張力,當(dāng)表面張力小于該臨界值時(shí),RT 不穩(wěn)定性現(xiàn)象才會(huì)發(fā)生.根據(jù)線性穩(wěn)定性理論建立了臨界表面張力的解析表達(dá)式,發(fā)現(xiàn)臨界表面張力隨著流體Atwood 數(shù)的增大而增大,且數(shù)值預(yù)測(cè)的臨界值與理論分析結(jié)果一致.另外發(fā)現(xiàn),界面張力對(duì)流體界面動(dòng)力學(xué)起著穩(wěn)定因子的作用,即增大表面張力可以抑制界面的卷吸與破裂行為,相界面的后期演化也從混沌的非對(duì)稱發(fā)展轉(zhuǎn)化為關(guān)于中軸線的對(duì)稱發(fā)展.還定量統(tǒng)計(jì)了表面張力對(duì)尖釘與氣泡振幅、演化速度和無量綱加速度的影響.在表面張力較小時(shí),界面張力對(duì)尖釘與氣泡增長的影響不再顯著,而繼續(xù)增大界面張力則可以抑制尖釘與氣泡的增長.根據(jù)速度演化曲線,高雷諾的三維RT 不穩(wěn)定性在所有界面張力情形下均呈現(xiàn)線性增長、飽和速度增長、重加速、混沌混合4 個(gè)不同的發(fā)展階段.尖釘與氣泡在飽和速度階段的漸進(jìn)速度與包含界面張力效應(yīng)的勢(shì)流理論模型相符合.進(jìn)一步發(fā)現(xiàn),隨著界面張力的增大,由于系統(tǒng)中形成旋渦的數(shù)目與強(qiáng)度逐漸減弱,導(dǎo)致不穩(wěn)定性進(jìn)入飽和速度增長和重加速階段的時(shí)間隨之推遲.最后,通過尖釘與氣泡振幅的開方與時(shí)間的線性擬合,從數(shù)值上證明了三維單模RT 不穩(wěn)定性在后期的混沌混合階段符合二次增長的規(guī)律,并且發(fā)現(xiàn)尖釘與氣泡后期增長率與表面張力呈現(xiàn)出遞減關(guān)系.

    猜你喜歡
    不穩(wěn)定性雷諾數(shù)表面張力
    可壓縮Navier-Stokes方程平面Couette-Poiseuille流的線性不穩(wěn)定性
    基于Transition SST模型的高雷諾數(shù)圓柱繞流數(shù)值研究
    神奇的表面張力
    小布老虎(2016年4期)2016-12-01 05:46:08
    MgO-B2O3-SiO2三元體系熔渣表面張力計(jì)算
    上海金屬(2016年2期)2016-11-23 05:34:45
    失穩(wěn)初期的低雷諾數(shù)圓柱繞流POD-Galerkin 建模方法研究
    基于轉(zhuǎn)捩模型的低雷諾數(shù)翼型優(yōu)化設(shè)計(jì)研究
    增強(qiáng)型體外反搏聯(lián)合中醫(yī)辯證治療不穩(wěn)定性心絞痛療效觀察
    民機(jī)高速風(fēng)洞試驗(yàn)的阻力雷諾數(shù)效應(yīng)修正
    前列地爾治療不穩(wěn)定性心絞痛療效觀察
    CaF2-CaO-Al2O3-MgO-SiO2渣系表面張力計(jì)算模型
    上海金屬(2014年3期)2014-12-19 13:09:06
    亚洲三区欧美一区| 啦啦啦免费观看视频1| 久久久久久亚洲精品国产蜜桃av| 久热这里只有精品99| 一级作爱视频免费观看| 黄频高清免费视频| 老司机靠b影院| 国产欧美日韩一区二区精品| www.www免费av| 一级片免费观看大全| 人成视频在线观看免费观看| 亚洲一码二码三码区别大吗| 亚洲男人天堂网一区| 男人舔女人的私密视频| 在线av久久热| 日本 av在线| 成人一区二区视频在线观看| 一进一出抽搐动态| 久久亚洲精品不卡| 啦啦啦免费观看视频1| 国内精品久久久久精免费| 亚洲欧洲精品一区二区精品久久久| 99精品久久久久人妻精品| 俺也久久电影网| 好男人在线观看高清免费视频 | 曰老女人黄片| 亚洲av五月六月丁香网| 亚洲一区中文字幕在线| 美女国产高潮福利片在线看| 老司机在亚洲福利影院| 妹子高潮喷水视频| 18禁观看日本| 国产一区在线观看成人免费| 黑丝袜美女国产一区| 亚洲成av人片免费观看| 两个人看的免费小视频| 一区二区三区激情视频| 久久中文看片网| 亚洲七黄色美女视频| 97超级碰碰碰精品色视频在线观看| 中文字幕久久专区| 老司机靠b影院| 男女做爰动态图高潮gif福利片| 精品久久久久久久毛片微露脸| 亚洲黑人精品在线| 两个人视频免费观看高清| 十八禁网站免费在线| 国产精品自产拍在线观看55亚洲| 免费看美女性在线毛片视频| 亚洲精品中文字幕一二三四区| 18禁美女被吸乳视频| 99国产综合亚洲精品| 精品国产乱码久久久久久男人| 男女之事视频高清在线观看| 18美女黄网站色大片免费观看| 麻豆成人午夜福利视频| 中出人妻视频一区二区| 亚洲天堂国产精品一区在线| 免费观看精品视频网站| 色播在线永久视频| 超碰成人久久| 在线观看www视频免费| 亚洲精品一区av在线观看| 在线免费观看的www视频| 啪啪无遮挡十八禁网站| 欧美久久黑人一区二区| 国产精品综合久久久久久久免费| 久久久精品欧美日韩精品| 午夜激情福利司机影院| 国产精品国产高清国产av| av欧美777| 老鸭窝网址在线观看| www.精华液| 精品久久久久久久人妻蜜臀av| 国产99久久九九免费精品| 欧美日韩福利视频一区二区| 免费高清视频大片| 欧美一区二区精品小视频在线| 免费一级毛片在线播放高清视频| 亚洲成a人片在线一区二区| 亚洲欧美一区二区三区黑人| 国产精品自产拍在线观看55亚洲| 色播在线永久视频| 一级作爱视频免费观看| 国产熟女午夜一区二区三区| 在线观看日韩欧美| 精品午夜福利视频在线观看一区| 最新美女视频免费是黄的| 亚洲精品av麻豆狂野| 亚洲欧美日韩高清在线视频| 9191精品国产免费久久| 成人av一区二区三区在线看| 欧美午夜高清在线| av欧美777| 国产熟女xx| 少妇的丰满在线观看| 亚洲一码二码三码区别大吗| 欧美午夜高清在线| av欧美777| 国产一级毛片七仙女欲春2 | 久久久精品国产亚洲av高清涩受| 中文字幕最新亚洲高清| 天天一区二区日本电影三级| 日韩精品中文字幕看吧| 97碰自拍视频| 777久久人妻少妇嫩草av网站| 欧美日韩黄片免| 亚洲成人久久性| 好男人电影高清在线观看| 国产一区在线观看成人免费| 天堂√8在线中文| 最近最新中文字幕大全免费视频| 午夜精品久久久久久毛片777| 一卡2卡三卡四卡精品乱码亚洲| 久久久久久久久免费视频了| 人人妻人人看人人澡| 国产精品免费视频内射| 亚洲精品一区av在线观看| 国产精品,欧美在线| 日本免费a在线| 久9热在线精品视频| 美女扒开内裤让男人捅视频| 精品卡一卡二卡四卡免费| 国产99白浆流出| 国产一区在线观看成人免费| 1024视频免费在线观看| 日本一本二区三区精品| 久久久久亚洲av毛片大全| 制服丝袜大香蕉在线| 色av中文字幕| 亚洲中文字幕一区二区三区有码在线看 | 99国产精品99久久久久| 午夜福利欧美成人| 麻豆成人午夜福利视频| 两性午夜刺激爽爽歪歪视频在线观看 | 草草在线视频免费看| www.www免费av| 一级毛片女人18水好多| 日本五十路高清| 久久精品91蜜桃| 99久久国产精品久久久| 可以在线观看的亚洲视频| 美女扒开内裤让男人捅视频| 亚洲午夜理论影院| 一本大道久久a久久精品| 久久精品91蜜桃| 十八禁网站免费在线| 午夜激情福利司机影院| 美女国产高潮福利片在线看| 老汉色av国产亚洲站长工具| 亚洲一区高清亚洲精品| 熟女电影av网| 国产精品,欧美在线| 免费观看人在逋| 老鸭窝网址在线观看| 男人的好看免费观看在线视频 | 午夜福利一区二区在线看| 亚洲熟妇中文字幕五十中出| 久久亚洲真实| 两人在一起打扑克的视频| 人妻丰满熟妇av一区二区三区| 国产精品1区2区在线观看.| 天天躁狠狠躁夜夜躁狠狠躁| 久久久久久国产a免费观看| 1024视频免费在线观看| 可以免费在线观看a视频的电影网站| 看片在线看免费视频| 一边摸一边做爽爽视频免费| 丁香欧美五月| 午夜福利免费观看在线| 1024手机看黄色片| 级片在线观看| 亚洲一区中文字幕在线| 中文字幕人成人乱码亚洲影| 欧美成人一区二区免费高清观看 | 国产精品1区2区在线观看.| 悠悠久久av| 一区福利在线观看| 中文字幕人妻丝袜一区二区| 99国产综合亚洲精品| 国产精品美女特级片免费视频播放器 | 女同久久另类99精品国产91| 国产一区在线观看成人免费| 亚洲九九香蕉| 无人区码免费观看不卡| 亚洲人成网站高清观看| av在线播放免费不卡| 在线天堂中文资源库| 18禁裸乳无遮挡免费网站照片 | 国产午夜福利久久久久久| a级毛片a级免费在线| 成年女人毛片免费观看观看9| 精品不卡国产一区二区三区| 精华霜和精华液先用哪个| 免费电影在线观看免费观看| 日韩视频一区二区在线观看| 69av精品久久久久久| 亚洲五月天丁香| 亚洲一区高清亚洲精品| 色哟哟哟哟哟哟| 日韩欧美免费精品| 国产99久久九九免费精品| 波多野结衣av一区二区av| 成人午夜高清在线视频 | 色综合站精品国产| 国产精品永久免费网站| 女人高潮潮喷娇喘18禁视频| 午夜福利视频1000在线观看| 天天一区二区日本电影三级| 国产精品日韩av在线免费观看| 欧美另类亚洲清纯唯美| 欧美+亚洲+日韩+国产| 满18在线观看网站| 人人妻人人看人人澡| 国产高清激情床上av| 1024视频免费在线观看| 亚洲av电影不卡..在线观看| 欧美最黄视频在线播放免费| 午夜日韩欧美国产| 国产精品亚洲一级av第二区| 久久精品人妻少妇| 99国产精品一区二区蜜桃av| 成人欧美大片| 成人三级黄色视频| 久久午夜亚洲精品久久| 亚洲天堂国产精品一区在线| 午夜免费鲁丝| 一级作爱视频免费观看| 亚洲国产中文字幕在线视频| 亚洲专区字幕在线| 久久精品人妻少妇| www.999成人在线观看| 久久久久九九精品影院| 久久久久精品国产欧美久久久| 午夜激情av网站| 免费看a级黄色片| 国产黄a三级三级三级人| 无遮挡黄片免费观看| 免费在线观看影片大全网站| 亚洲国产欧洲综合997久久, | 国产激情偷乱视频一区二区| 黄色 视频免费看| 国产欧美日韩一区二区三| 国产精品久久视频播放| 亚洲 欧美一区二区三区| 不卡一级毛片| 88av欧美| 久99久视频精品免费| 女性被躁到高潮视频| АⅤ资源中文在线天堂| 正在播放国产对白刺激| 午夜a级毛片| 国产成人av教育| 国语自产精品视频在线第100页| 久久久久久国产a免费观看| 国产成人啪精品午夜网站| 亚洲全国av大片| 一级片免费观看大全| 国产成人精品久久二区二区免费| 欧美黄色片欧美黄色片| 国产久久久一区二区三区| 久久久久亚洲av毛片大全| 久久精品成人免费网站| 国产精品乱码一区二三区的特点| 在线天堂中文资源库| 99热这里只有精品一区 | 久久久久免费精品人妻一区二区 | 欧美国产日韩亚洲一区| 久久国产乱子伦精品免费另类| 嫁个100分男人电影在线观看| 老司机深夜福利视频在线观看| 亚洲熟妇中文字幕五十中出| 精品熟女少妇八av免费久了| 久久久久九九精品影院| 狠狠狠狠99中文字幕| 麻豆国产av国片精品| 久久精品国产99精品国产亚洲性色| 欧美大码av| 欧美激情久久久久久爽电影| 精华霜和精华液先用哪个| 精品免费久久久久久久清纯| 精品久久蜜臀av无| 色婷婷久久久亚洲欧美| 亚洲片人在线观看| 中国美女看黄片| 国产激情偷乱视频一区二区| 亚洲av成人av| 91国产中文字幕| av福利片在线| 欧美中文综合在线视频| 97超级碰碰碰精品色视频在线观看| 精品久久久久久久人妻蜜臀av| 999精品在线视频| 亚洲熟妇中文字幕五十中出| 中文字幕人妻熟女乱码| 久久久久久久久中文| 国产精品久久久久久亚洲av鲁大| 琪琪午夜伦伦电影理论片6080| a在线观看视频网站| 91麻豆精品激情在线观看国产| 99riav亚洲国产免费| 校园春色视频在线观看| 亚洲国产日韩欧美精品在线观看 | 老汉色∧v一级毛片| 90打野战视频偷拍视频| 欧美中文日本在线观看视频| ponron亚洲| 视频在线观看一区二区三区| 在线十欧美十亚洲十日本专区| 变态另类丝袜制服| 亚洲人成网站高清观看| 免费电影在线观看免费观看| 黄色片一级片一级黄色片| 日韩欧美国产在线观看| 亚洲精品国产一区二区精华液| 国产亚洲精品一区二区www| 日本五十路高清| 亚洲狠狠婷婷综合久久图片| 亚洲成国产人片在线观看| 波多野结衣巨乳人妻| 婷婷精品国产亚洲av在线| 色综合站精品国产| 欧美黄色淫秽网站| 欧美精品亚洲一区二区| 亚洲男人天堂网一区| 中文字幕精品亚洲无线码一区 | 久久人人精品亚洲av| 在线观看免费午夜福利视频| 中文亚洲av片在线观看爽| 欧美成人性av电影在线观看| 久久精品国产综合久久久| 波多野结衣高清作品| 午夜成年电影在线免费观看| avwww免费| 老司机午夜十八禁免费视频| 好男人电影高清在线观看| av福利片在线| 久久久久久九九精品二区国产 | 两个人看的免费小视频| 午夜福利18| 91国产中文字幕| 一级毛片精品| av天堂在线播放| 国产成人精品久久二区二区91| 最新美女视频免费是黄的| 黄网站色视频无遮挡免费观看| 听说在线观看完整版免费高清| www国产在线视频色| 国产精品国产高清国产av| 欧美午夜高清在线| 中文字幕人妻熟女乱码| 老司机深夜福利视频在线观看| 宅男免费午夜| 99久久精品国产亚洲精品| 久热爱精品视频在线9| 中文字幕人妻丝袜一区二区| 一夜夜www| 久久精品91无色码中文字幕| 我的亚洲天堂| 成人国产一区最新在线观看| 人成视频在线观看免费观看| 成人精品一区二区免费| 亚洲精品久久成人aⅴ小说| 在线看三级毛片| 精品少妇一区二区三区视频日本电影| 视频区欧美日本亚洲| 午夜久久久久精精品| 啦啦啦免费观看视频1| 热99re8久久精品国产| 亚洲人成网站在线播放欧美日韩| 嫩草影院精品99| 高潮久久久久久久久久久不卡| 啦啦啦免费观看视频1| 日韩有码中文字幕| 欧美精品亚洲一区二区| 听说在线观看完整版免费高清| 最近最新中文字幕大全电影3 | 黄色女人牲交| 精品无人区乱码1区二区| 久热爱精品视频在线9| 淫秽高清视频在线观看| 国产人伦9x9x在线观看| 一级黄色大片毛片| 欧美激情久久久久久爽电影| 欧美激情 高清一区二区三区| x7x7x7水蜜桃| 亚洲国产日韩欧美精品在线观看 | 久久国产乱子伦精品免费另类| 99国产综合亚洲精品| 免费在线观看日本一区| 长腿黑丝高跟| 久久 成人 亚洲| 国产一区二区三区视频了| 成人手机av| 欧美性长视频在线观看| 男女视频在线观看网站免费 | 国产成+人综合+亚洲专区| 国内毛片毛片毛片毛片毛片| 国产精品电影一区二区三区| 波多野结衣高清作品| 成人18禁高潮啪啪吃奶动态图| 久久久久久九九精品二区国产 | 免费女性裸体啪啪无遮挡网站| 亚洲人成伊人成综合网2020| 久久国产精品影院| 88av欧美| 亚洲人成伊人成综合网2020| 一级毛片精品| 亚洲av成人不卡在线观看播放网| 女性生殖器流出的白浆| 97超级碰碰碰精品色视频在线观看| 99精品在免费线老司机午夜| 国产亚洲精品一区二区www| 精品卡一卡二卡四卡免费| 欧美日本亚洲视频在线播放| 久久天堂一区二区三区四区| a在线观看视频网站| 成人三级做爰电影| 一本大道久久a久久精品| 天天躁狠狠躁夜夜躁狠狠躁| 免费看十八禁软件| 久久热在线av| 成人av一区二区三区在线看| 成人亚洲精品一区在线观看| e午夜精品久久久久久久| 韩国av一区二区三区四区| 国产成人系列免费观看| 久久久久国内视频| 日日摸夜夜添夜夜添小说| 亚洲精品国产区一区二| 久久久久久大精品| 91大片在线观看| 国产在线观看jvid| 亚洲国产中文字幕在线视频| 啦啦啦免费观看视频1| 丰满人妻熟妇乱又伦精品不卡| 99久久精品国产亚洲精品| 久久中文字幕一级| 天天躁夜夜躁狠狠躁躁| 中文字幕另类日韩欧美亚洲嫩草| 男女那种视频在线观看| 亚洲国产精品sss在线观看| 亚洲成a人片在线一区二区| 男人舔奶头视频| 老司机午夜十八禁免费视频| 国产高清视频在线播放一区| 国产伦人伦偷精品视频| 欧美大码av| 巨乳人妻的诱惑在线观看| 男人操女人黄网站| 欧美不卡视频在线免费观看 | 国产片内射在线| 国产亚洲精品av在线| 黄片小视频在线播放| www.熟女人妻精品国产| 神马国产精品三级电影在线观看 | 日日干狠狠操夜夜爽| 精品国产乱子伦一区二区三区| 久久中文看片网| 99热6这里只有精品| 国产一区二区三区视频了| 午夜精品在线福利| 精品福利观看| av在线天堂中文字幕| 欧美久久黑人一区二区| 国产av在哪里看| 中亚洲国语对白在线视频| 看免费av毛片| 1024香蕉在线观看| 国产又黄又爽又无遮挡在线| www.自偷自拍.com| 亚洲国产精品成人综合色| 国产av一区在线观看免费| 村上凉子中文字幕在线| 久久香蕉激情| 手机成人av网站| 国产亚洲欧美在线一区二区| 99久久国产精品久久久| 亚洲欧美激情综合另类| 欧美国产日韩亚洲一区| 好男人电影高清在线观看| 久久久久九九精品影院| 国产伦一二天堂av在线观看| 欧美日韩乱码在线| 亚洲成人免费电影在线观看| 视频在线观看一区二区三区| 嫩草影视91久久| 熟女少妇亚洲综合色aaa.| 亚洲国产欧洲综合997久久, | 久久精品aⅴ一区二区三区四区| 亚洲在线自拍视频| 欧美午夜高清在线| 亚洲av成人一区二区三| 最近最新免费中文字幕在线| 亚洲成人免费电影在线观看| xxxwww97欧美| 男女做爰动态图高潮gif福利片| 国内久久婷婷六月综合欲色啪| 国产一区二区在线av高清观看| 亚洲自偷自拍图片 自拍| 男女之事视频高清在线观看| 亚洲专区字幕在线| 亚洲中文av在线| 亚洲熟妇熟女久久| 亚洲精品久久国产高清桃花| 日韩大尺度精品在线看网址| 男人舔女人下体高潮全视频| 熟女少妇亚洲综合色aaa.| 国产亚洲精品第一综合不卡| e午夜精品久久久久久久| 亚洲三区欧美一区| 淫秽高清视频在线观看| 精品一区二区三区av网在线观看| 亚洲专区中文字幕在线| 色播亚洲综合网| 国产激情偷乱视频一区二区| 国产亚洲欧美在线一区二区| 国产乱人伦免费视频| 亚洲精品久久国产高清桃花| 亚洲五月色婷婷综合| 免费女性裸体啪啪无遮挡网站| 亚洲av电影不卡..在线观看| 久久午夜综合久久蜜桃| 国产99白浆流出| 69av精品久久久久久| 亚洲精品一卡2卡三卡4卡5卡| 欧美日韩一级在线毛片| 亚洲精品在线美女| 最近在线观看免费完整版| 久久久久久久久免费视频了| 日日爽夜夜爽网站| 中文字幕人成人乱码亚洲影| 亚洲人成电影免费在线| 黄色成人免费大全| 老鸭窝网址在线观看| 欧美乱码精品一区二区三区| 精品国产美女av久久久久小说| 精品一区二区三区四区五区乱码| 欧美日韩精品网址| 亚洲一码二码三码区别大吗| 国产免费男女视频| 满18在线观看网站| 欧美一区二区精品小视频在线| 亚洲一卡2卡3卡4卡5卡精品中文| 此物有八面人人有两片| 色尼玛亚洲综合影院| 久久亚洲精品不卡| 色综合婷婷激情| 成人亚洲精品一区在线观看| 香蕉久久夜色| 亚洲成av人片免费观看| 国产高清视频在线播放一区| 日韩欧美在线二视频| 他把我摸到了高潮在线观看| 搡老熟女国产l中国老女人| 日韩欧美免费精品| 色哟哟哟哟哟哟| 三级毛片av免费| 婷婷六月久久综合丁香| 成人三级做爰电影| 99国产极品粉嫩在线观看| 村上凉子中文字幕在线| 欧美zozozo另类| 久久国产亚洲av麻豆专区| 欧美激情 高清一区二区三区| 欧洲精品卡2卡3卡4卡5卡区| 中文亚洲av片在线观看爽| 男人操女人黄网站| 一二三四在线观看免费中文在| 色播亚洲综合网| 国产1区2区3区精品| 国产一区二区在线av高清观看| 国语自产精品视频在线第100页| 在线观看免费日韩欧美大片| 美女午夜性视频免费| 国产高清videossex| av在线天堂中文字幕| 国产精品国产高清国产av| 国产精品久久久久久人妻精品电影| 少妇 在线观看| 色播亚洲综合网| 亚洲av片天天在线观看| 亚洲第一青青草原| www.精华液| 久久亚洲精品不卡| www.熟女人妻精品国产| 久久欧美精品欧美久久欧美| 久久久国产成人精品二区| 欧美在线黄色| 亚洲成国产人片在线观看| 日韩欧美国产一区二区入口| 桃红色精品国产亚洲av| 亚洲va日本ⅴa欧美va伊人久久| 成人午夜高清在线视频 | 黄色 视频免费看| 禁无遮挡网站| 久久久久亚洲av毛片大全| 老司机福利观看| 精品久久久久久久毛片微露脸| 久久精品亚洲精品国产色婷小说| 亚洲黑人精品在线| 夜夜看夜夜爽夜夜摸| 午夜福利成人在线免费观看| 亚洲人成77777在线视频| 欧美性长视频在线观看| 热99re8久久精品国产| 一a级毛片在线观看| 黄色毛片三级朝国网站| 不卡av一区二区三区| 久久九九热精品免费| 欧美黑人精品巨大| 中文字幕av电影在线播放| 一本久久中文字幕| 一a级毛片在线观看| 中文字幕久久专区|