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

    平面激波沖擊并排液滴的三維數(shù)值模擬研究1)

    2022-12-18 06:11:16吳潤龍李祝軍丁航
    力學(xué)學(xué)報 2022年11期
    關(guān)鍵詞:激波液滴氣流

    吳潤龍 李祝軍 丁航

    (中國科學(xué)技術(shù)大學(xué)近代力學(xué)系,合肥 230026)

    引言

    在高速氣流沖擊下,懸浮空氣中的液滴會發(fā)生持續(xù)變形并最終可能破碎.這一非常復(fù)雜的界面演化和流動耦合問題廣泛存在于工程應(yīng)用和武器設(shè)計之中,比如液態(tài)化學(xué)武器的高速拋灑[1]、超聲速飛行時飛行器表面的雨滴侵蝕防護(hù)[2]、火箭引擎中液態(tài)燃料的摻混過程[3]等.因此,激波沖擊下液滴破碎的過程對工業(yè)生產(chǎn)實踐具有非常重要的指導(dǎo)意義.另一方面,由于液滴在激波沖擊下的氣動變形和破碎過程存在大密度比、多介質(zhì)和復(fù)雜波系結(jié)構(gòu),甚至出現(xiàn)界面拓?fù)鋷缀巫兓拖嘧兊葟?fù)雜流動問題,給研究帶來了極大的挑戰(zhàn).

    關(guān)于液滴氣動破碎問題的研究已有廣泛的實驗研究,在液滴演化形態(tài)特征、液滴變形與破碎的物理機(jī)制、破碎特征時間[4]以及破碎后形成的子液滴尺寸分布等方面得到了豐富的研究成果.Hinze[3]采用韋伯?dāng)?shù)(W e=ρU2D/σ)來分析液滴在不同流動條件下的液滴氣動破碎機(jī)制,其中 ρ 和U為來流的密度和特征速度,D和 σ 為液滴的直徑和表面張力系數(shù).韋伯?dāng)?shù)表征了來流的氣動慣性作用與液滴表面張力作用的相對大小.研究發(fā)現(xiàn),隨著韋伯?dāng)?shù)的增加,液滴會出現(xiàn)不同的破碎模式.Pilch等[5]以及Guildenbecher 等[6]給出了5 種不同破碎模式的劃分,包括振蕩破碎[7-8](vibrational)、袋狀破碎[9-10](bag)、袋蕊破碎[11](bag-and-stamen)、剪切破碎[12-13](shear)和毀滅性破碎[14-16](catastrophic).在較小韋伯?dāng)?shù)(<100)下,液滴發(fā)生震蕩破碎、袋狀破碎、或袋蕊破碎.在較大韋伯?dāng)?shù)下(>100)液滴則發(fā)生剪切破碎,并生成大量小液滴[17].韋伯?dāng)?shù)>350 以上的毀滅性破碎模態(tài)是否存在目前尚有爭議[18].

    前人針對液滴破碎過程提出了不同的理論模型來進(jìn)行解釋.Nicholls 等[12]提出了“剪切剝離”(shearstripping)模型,認(rèn)為在高速氣流沖擊下,氣?液兩相間會形成邊界層,邊界層的剪切作用使液滴在迎風(fēng)面赤道附近不斷剝離出子液滴并最終導(dǎo)致液滴整體破碎.Liu 等[19]提出了“薄層細(xì)化”(sheet-thinning)模型,他們認(rèn)為是高速氣流使液滴迎風(fēng)面赤道處液體向流向彎曲形成薄片,液體薄片因毛細(xì)不穩(wěn)定性而形成拉絲,并最終斷裂生成大量子液滴.Theofanous 等[18,20]在實驗中采用激光誘導(dǎo)熒光技術(shù)研究了高韋伯?dāng)?shù)下液滴的氣動破碎過程,提出液滴破碎主要由兩種機(jī)制所主導(dǎo),分別為“瑞利-泰勒穿透”(Rayleigh-Taylor (RT) piercing)和“剪切誘導(dǎo)夾帶”(shear-induced entrainment).前者是由RT 不穩(wěn)定性(RT instability)與表面張力起主導(dǎo)作用導(dǎo)致液滴破碎,相應(yīng)韋伯?dāng)?shù)范圍是10 1000.而當(dāng)韋伯?dāng)?shù)處100~1000 之間時,則是兩種機(jī)制相互競爭.

    由于受到實驗條件的限制,實驗研究對于液滴形態(tài)時空演化的流動細(xì)節(jié)難以清楚刻畫,特別是對高馬赫數(shù)、高韋伯?dāng)?shù)下液滴發(fā)生剪切破碎的研究.此時界面波的尺度較小且實驗觀測受到液霧遮擋的影響,因此對界面演化和液滴內(nèi)部流場的測量極其困難.近年來,人們開始使用直接數(shù)值模擬的方法來研究液滴破碎問題,以解決以上問題.Sembian 等[21]通過實驗并結(jié)合數(shù)值模擬研究了水柱在激波沖擊下內(nèi)部波系的演化過程.他們在實驗中觀察到水柱內(nèi)生成了空化氣泡,結(jié)合數(shù)值結(jié)果分析后指出: 膨脹波在液滴內(nèi)匯聚產(chǎn)生低壓區(qū),導(dǎo)致液相空化并形成氣泡.Meng 等[22]使用擴(kuò)散界面方法數(shù)值研究了三維液滴的剪切破碎模式,描述了剪切破碎機(jī)制下液滴的變形和破碎特征,發(fā)現(xiàn)界面的RT 不穩(wěn)定性會使回流區(qū)紊亂.Liu 等[23]使用五方程模型模擬了二維和三維情況下液滴的剪切破碎過程,并根據(jù)數(shù)值結(jié)果將液滴變形破碎分成3 個不同的階段,分別是不穩(wěn)定性生成階段、液滴扁平化階段以及子液滴剝離階段.Dorschner 等[24]使用擴(kuò)散界面方法數(shù)值模擬了在超聲速流場中帶黏性和表面張力的水滴的變形破碎過程,著重分析了迎風(fēng)面赤道液絲的演化過程,并進(jìn)一步研究了影響液絲循環(huán)斷裂過程的因素.

    實際問題中流場往往存在眾多液滴,并且液滴間可能會發(fā)生相互作用.Yoshida 等[25]通過實驗研究了在入射激波馬赫數(shù)為1.4 在流向上相距為8.8 mm 的兩個直徑都為1.4 mm 的球形液滴的演化過程,他們發(fā)現(xiàn)前面的液滴位移要顯著大于后面的液滴.Igra 等[26]使用全息干涉儀可視化了在激波沖擊下兩個串聯(lián)柱狀水滴在破碎過程各個階段的相互作用.他們發(fā)現(xiàn)前面液滴的變形、位移以及加速度與單個液滴相似,但前面水柱的存在使得后面液滴的變形、位移以及阻力系數(shù)都要小于前面液滴,并且后面液滴達(dá)到最大變形的時間要晚于前面液滴.Chen 等[27]使用基于五方程模型的多相流求解器分別模擬了激波沖擊兩個水平排列水柱以及兩個垂直排列水柱的演化過程,發(fā)現(xiàn)液滴間相互作用與激波在液滴間的演化密切相關(guān).總體來說,相對于單液滴氣動破碎研究,目前多液滴氣動耦合變形和破碎研究還相對較少,特別是三維多液滴之間的相互作用機(jī)制及其對氣動變形帶來的影響有待進(jìn)一步研究.這些研究將有助于高速氣流沖擊下液滴云兩相流模型中相間作用建模.

    本文采用基于三維切割網(wǎng)格算法的二階守恒型清晰界面方法,研究兩個并排球形甘油液滴在激波沖擊下的動力學(xué)演化過程.由于考慮高速氣流沖擊作用,流動建模中忽略了液滴的表面張力和黏性作用,因此流動控制方程為三維無黏可壓縮歐拉方程.研究中著重探究激波沖擊后的復(fù)雜波系結(jié)構(gòu)生成以及并排液滴耦合作用所帶來的界面非對稱演化.此外,還研究了兩液滴的間距參數(shù)對通道側(cè)界面形態(tài)以及液滴形態(tài)演化的影響.

    1 物理問題描述及數(shù)值方法介紹

    1.1 問題介紹

    如圖1 所示,考慮的物理問題是兩個大小相同、并排分布的甘油液滴受到空氣中平面激波沖擊后的動力學(xué)演化過程.初始液滴的半徑固定為R=0.895 mm,液滴間距離為L.入射激波馬赫數(shù)為Ms=2.67,初始環(huán)境壓力為一個標(biāo)準(zhǔn)大氣壓(101 325 Pa).甘油密度為1263 kg/m3,波前空氣密度為1.25 kg/m3.計算中采用了拉伸結(jié)構(gòu)網(wǎng)格(416 × 316 × 316),其中在液滴附近區(qū)域采用均勻密網(wǎng)格,其網(wǎng)格解析度為每個初始半徑內(nèi)分布50 個網(wǎng)格,在遠(yuǎn)離界面的計算區(qū)域則使用較粗的網(wǎng)格解析度.

    1.2 控制方程

    流體運(yùn)動的控制方程為三維歐拉方程

    其中守恒量Q和通量F分別為

    其中,ρ是密度,u=(u,v,w)T是流體速度,p是壓強(qiáng),E=e+|u|2/2 是比總能,e為比內(nèi)能,I為單位矩陣.

    本文使用剛性氣體狀態(tài)方程[28]來近似液體的本構(gòu)關(guān)系,以便封閉控制方程.具體表達(dá)式如下

    其中,γ是比熱比,P∞是和流體熱力學(xué)性質(zhì)相關(guān)的一個常數(shù).對于本文所涉及的流體,空氣取γ=1.4,P∞=0,甘油取γ=1.504 8,P∞=299 MPa[29].

    1.3 數(shù)值方法介紹

    本文使用實驗室已有的三維可壓縮守恒清晰界面數(shù)值方法[30]對激波沖擊并排液滴問題進(jìn)行數(shù)值模擬,以下對該數(shù)值方法進(jìn)行簡單介紹.

    該方法使用水平集(level set)方法[31]來捕捉兩相界面的位置及其演化,其中界面演化的水平集方程為

    式中,φ (x,t)為符號距離函數(shù),其取值為零的等值面表示界面位置,vφ為界面速度.基于水平集函數(shù)所表征的界面位置,采用切割網(wǎng)格方法對背景結(jié)構(gòu)網(wǎng)格進(jìn)行切割和組裝,在界面附近生成非結(jié)構(gòu)網(wǎng)格,并保證重構(gòu)界面與非結(jié)構(gòu)網(wǎng)格單元面時時重合;在遠(yuǎn)離界面處仍然采用背景網(wǎng)格.關(guān)于切割網(wǎng)格方法的具體細(xì)節(jié)可以參考文獻(xiàn)[30,32-33].

    使用任意拉格朗日?歐拉框架(arbitrary Lagrangian-Eulerian)對歐拉方程進(jìn)行有限體積離散.在網(wǎng)格單元面上的通量計算中,單相的黎曼問題求解使用AUSM+-up 方法[34],兩相的黎曼問題求解則使用精確黎曼求解器[29].對于水平集方程中的時間項,使用2 階龍格庫塔法進(jìn)行離散,距離符號函數(shù)的空間梯度項使用5 階WENO 格式[35]進(jìn)行離散.考慮到在兩相界面發(fā)生拓?fù)鋷缀巫兓虺霈F(xiàn)長時間計算的情況下,φ 場會偏離符號距離函數(shù).為了保持符號距離函數(shù)性質(zhì),在本文中求解了符號距離函數(shù)的重新再初始化方程[36]

    式中,τ為虛擬時間,φ0是重新初始化前(τ=0)的 φ場,s(φ0)是符號距離函數(shù).對于空間項 ? φ,這里使用二階ENO 格式[37-38]來離散,而時間項離散則使用顯式的二階龍格?庫塔法.

    關(guān)于數(shù)值程序的驗證,Shen 等[30]已將馬赫數(shù)為2.67 的激波沖擊單個甘油液滴后界面初期演化的數(shù)值結(jié)果與Theofanous 等[39]的實驗結(jié)果進(jìn)行了對比,證實了兩者在液滴的整體形態(tài)與局部界面特征上基本吻合.由于本工作采用了基本一致的物理參數(shù),并使用同一套三維數(shù)值程序?qū)げ_擊并排甘油液滴進(jìn)行數(shù)值模擬研究,因此可以認(rèn)為數(shù)值結(jié)果可靠.

    2 結(jié)果分析

    2.1 沖擊初期波系演化

    對于激波沖擊并排液滴動力學(xué)問題,液滴外流場分為兩個區(qū)域,一個是由兩液滴相鄰區(qū)域構(gòu)成的通道區(qū)(如圖1 所示),另一個則是不相鄰區(qū)域構(gòu)成的非通道區(qū).當(dāng)入射激波IS 接觸液滴后,激波后續(xù)演化在通道側(cè)和非通道側(cè)是非對稱的.圖2 給出了從入射激波接觸液滴后前5 μs 時間內(nèi)波系演化情況.在1 μs 時,入射激波IS 在界面上發(fā)生了反射和透射,反射和透射的波系類型由兩相流體的聲阻抗Z=ρc決定,其中c為聲速.具體地,透射到甘油液滴中的聲波強(qiáng)度可以由空氣和甘油之間的透射系數(shù)決定

    其中,下標(biāo)a表示空氣,d表示甘油液滴.根據(jù)上游駐點(diǎn)處空氣和甘油液滴的狀態(tài),可以估算出Id/Ia=0.003,這說明來流空氣中的能量只有極少量的部分透射進(jìn)液滴,而絕大部分能量則反射到氣相中,在液滴上游形成了弓形反射激波RS.而透射進(jìn)液滴的波系為透射壓縮波TCW.由于壓縮波的傳播速度為當(dāng)?shù)芈曀?且甘油液滴內(nèi)的波前聲速為cd=597 m/s,小于空氣中入射激波IS 的運(yùn)動速度vIS=Mscg=899 m/s,因此在圖2(b)中TCW 的形狀是向左凸的.在赤道界面(X=0)附近,入射激波的反射模式由于入射激波與界面夾角的變化,變成了馬赫反射,并形成了馬赫桿MS.和非通道側(cè)不同的是,在通道側(cè),入射激波IS 連接了兩個液滴的馬赫桿MS.在2 μs 時,通道側(cè)出現(xiàn)了新的波系結(jié)構(gòu).由于兩個液滴各自的弓形反射激波RS 向上游傳播,在通道側(cè)兩道RS 發(fā)生了相交,進(jìn)而產(chǎn)生了兩道新的反射激波RS2.由激波動力學(xué)可以知道,在兩道RS 和兩道RS2 形成的交叉點(diǎn)處必然會形成馬赫桿,從而讓四道波陣面前后的壓力、速度、密度這些流場變量得到匹配.因此,在3 μs 時可以看到,在通道上游出現(xiàn)了馬赫桿MS2,并且MS2 連接了RS 和RS2,形成了兩個對稱分布的三波點(diǎn).而在通道側(cè)下游,當(dāng)MS 跨過赤道(X=0)向下游傳播時,受彎曲界面上產(chǎn)生的稀疏波的作用發(fā)生彎曲.與非通道側(cè)不同的是,連接馬赫桿MS 的入射激波IS 在通道側(cè)也發(fā)生了彎曲,并向下游略微凸起.這是因為通道的存在導(dǎo)致入射激波IS 在通道側(cè)相比非通道側(cè)距離界面更近,從而也更容易受到赤道后曲界面處的稀疏波作用,進(jìn)而波陣面發(fā)生彎曲.此外,隨著液滴兩側(cè)彎曲馬赫桿MS 向下游的不斷運(yùn)動,最終它們在液滴后駐點(diǎn)處相遇,形成反射繞射激波RDS.從2~3 μs 時間內(nèi)波系演化過程中還可以發(fā)現(xiàn),隨著MS 不斷向下游傳播,入射激波最終與兩道馬赫桿MS 發(fā)生了合并形成了融合激波CS.在4 μs 時,如圖 2(e)中所示,反射繞射激波RDS 與融合激波CS 的相交點(diǎn)會進(jìn)一步形成馬赫桿MS3.

    圖2 沖擊初期波系演化過程 (L/R=1)Fig.2 The evolution of wave structures at the early stage of impact (L/R =1)

    在5 μs 時,從圖2(f)中可以看到,通道側(cè)上游的馬赫桿MS2 波陣面較為平整,而兩道反射激波RS2 的強(qiáng)度明顯減弱.相比通過非通道側(cè)RS 的彎曲波陣面,上游來流通過MS2 的波陣面后的流動改變是截然不同的.通過RS 會導(dǎo)致氣流向偏離液滴的方向偏轉(zhuǎn),而通過MS2 的平直波陣面的氣流方向則不發(fā)生改變,這也必然會誘導(dǎo)通道側(cè)與非通道側(cè)的界面形態(tài)發(fā)生持續(xù)的非對稱演化.在液滴下游RDS 反向繞射時,通道側(cè)的部分與背風(fēng)面的分離激波融合形成了一道強(qiáng)度更強(qiáng)的分離激波,并在兩道激波連接處形成一個折點(diǎn)TP,而非通道側(cè)的部分則未發(fā)生融合.此外,需要指出的是: 在平面激波沖擊早期(0~ 5 μs),由于液滴的變形可以忽略不計,而且入射激波透射進(jìn)入液滴內(nèi)的能量極其微弱,因此平面激波沖擊并排液滴后液滴外部波系的演化過程與激波沖擊并排剛性小球的波系演化過程基本上是完全一致的.

    通道的存在使得液滴兩側(cè)的流動波系出現(xiàn)了非對稱結(jié)構(gòu),進(jìn)而也導(dǎo)致了壓力場的非對稱變化.圖3顯示了與波系演化相對應(yīng)的前5 μs 的壓力場變化過程.這里采用壓力系數(shù)

    來表征壓力場的變化,其中p為當(dāng)?shù)貕毫?p0為初始環(huán)境壓力,pd為初始入射激波的波后氣流壓力.圖3顯示,在2 μs 時兩個并排液滴的弓形反射激波在通道下游側(cè)相交,形成了由兩道反射激波RS2 圍成的高壓區(qū).但由于受到下游氣流膨脹的影響,該區(qū)域壓強(qiáng)相對于初始液滴赤道附近的壓強(qiáng)要低,且隨著時間發(fā)展持續(xù)降低.可以看到由于弓形反射激波相交形成的高壓在3 μs 時已經(jīng)逐漸減弱,在4 μs 時基本和非通道側(cè)壓力區(qū)別不大.相對而言,在通道側(cè)液滴上游的高壓區(qū)則仍然得到了保持.從圖3(e)和圖3(f)中 (即4~ 5 μs) 壓力云圖演化可以看到,在液滴間通道的下游形成了一個非常明顯的低壓區(qū).該低壓區(qū)是由兩個液滴間喉道后氣流的快速膨脹所導(dǎo)致的.而且,由于該低壓區(qū)同時受到兩個液滴背風(fēng)面稀疏波的共同作用,其壓力要比液滴非通道側(cè)對應(yīng)位置的壓力更低.

    圖3 并排液滴與激波相互作用早期的壓力云圖 (L/R=1)Fig.3 The evolution of wave structures at the early stage of impact (L/R=1)

    總體來說,在高速氣流沖擊作用下,由于液滴和周圍空氣的密度差異大,早期液滴界面沒有呈現(xiàn)明顯的變形,但是并排液滴的幾何效應(yīng)對波系結(jié)構(gòu)的發(fā)展和壓力場的建立起到了顯著作用.

    2.2 沖擊中期界面形態(tài)演化

    圖4 和圖5 給出了并排雙液滴在平面激波沖擊下的三維形態(tài)學(xué)演化過程及X-Y平面(Z=0)上的壓力及馬赫數(shù)分布.其中圖4 顯示了3 μs,7 μs 和10 μs的液滴側(cè)面和正面視圖.雖然激波及波后來流在3 μs 時都已經(jīng)繞過了液滴到達(dá)了液滴下游(見圖2(d)),但是液滴整體界面形態(tài)并未發(fā)生明顯變化,如圖4(a1)和圖4(a2)所示.雖然在7 μs 時入射激波已遠(yuǎn)離液滴,但在周圍氣流的持續(xù)作用下,液滴界面形態(tài)開始出現(xiàn)輕微的變化.此外,在平面激波沖擊下由于并排液滴的反射與繞射波系的相互作用,壓力和馬赫數(shù)分布對于通道側(cè)和非通道側(cè)是不對稱的.如圖4(b1)和圖4(c1)所示,在兩液滴的通道側(cè)下游,馬赫數(shù)明顯大于非通道側(cè)的馬赫數(shù).圖4(b1)和圖4(b2)中,在液滴迎風(fēng)面出現(xiàn)了一圈圈漣漪狀的表面波;同時,在液滴背面以后駐點(diǎn)為中心出現(xiàn)了一圈凹陷狀態(tài)的界面變形,但這一變形對于通道側(cè)和非通道側(cè)明顯是不對稱的.具體來說,非通道側(cè)凹陷的區(qū)域更大,而通道側(cè)凹陷的區(qū)域更小;這表明通道側(cè)與非通道側(cè)的不對稱流動對界面形態(tài)演化產(chǎn)生了影響.圖4(c1)和圖4(c2)顯示在10 μs 時,液滴的迎風(fēng)面的表面波進(jìn)一步增長,在迎風(fēng)面中間出現(xiàn)了一道凹槽結(jié)構(gòu);該現(xiàn)象與單液滴實驗中觀察到的情形[39]是類似的.此外,液滴背風(fēng)面的界面褶皺線出現(xiàn)了一定的不對稱性.靠近通道側(cè)的液滴側(cè)面褶皺線是與來流方向基本上是垂直的,而在非通道側(cè)的褶皺線與來流方向存在一個較為明顯的夾角.

    圖4 并排液滴的三維形態(tài)演化及X-Y 平面(Z=0)上的壓力與馬赫數(shù)分布,標(biāo)號為1 的列為液滴側(cè)面,標(biāo)號為2 的列為液滴正面(L/R=1)Fig.4 Three-dimensional morphological evolution of side-by-side droplets,pressure and Mach number contours in the X-Y plane (Z=0).The column with label 1 corresponds to the side of the droplet,and the column with the label 2 corresponds to the front of the droplet (L/R=1)

    圖5 顯示了13 μs,16 μs 和19 μs 的液滴側(cè)面和正面視圖.如圖5(a1)和圖5(a2)所示,兩個并排液滴由于流向上的壓差作用,整體逐漸被壓扁,特別是液滴后駐點(diǎn)附近的界面變得更為扁平.隨著時間的發(fā)展,如圖5(b1)和圖5(c1)所示,通道側(cè)與非通道側(cè)在液滴界面演化的差別越來越顯著.從側(cè)面視圖可以看出,由于界面KH 不穩(wěn)定性和空氣動力學(xué)的共同作用,在初始液滴赤道附近出現(xiàn)山脊?fàn)盥∑?其中非通道側(cè)的隆起部分相較于通道側(cè)更為明顯.

    圖5 激波沖擊并排液滴(L/R=1)后界面三維形態(tài)演化及X-Y 平面(Z=0)上的壓力與馬赫數(shù)分布.其中液滴側(cè)面視圖標(biāo)號為 1,液滴正面標(biāo)號為 2Fig.5 Snapshots of side-by-side droplets (L/R=1) after being impacted by a planar shock,pressure and Mach number contours in the X-Y plane(Z=0).The column with label 1 shows the side view of the droplet,and the column with the label 2 shows the front view

    為了定量衡量液滴間相互作用對液滴形態(tài)演化的影響,圖6 給出了并排雙液滴(L/R=1)以及單個液滴的幾何參數(shù)隨時間變化的對比圖.并排雙液滴的幾何參數(shù)采用在垂直來流方向上的通道側(cè)半徑和非通道側(cè)半徑(如圖5(c1)所示),這里半徑指的是流動方向界面最大截面中液滴一側(cè)界面到中軸線的最大距離.從圖6 中并排液滴非通道側(cè)半徑與單液滴半徑的對比,可以發(fā)現(xiàn)兩者半徑隨時間變化的規(guī)律基本是一樣的.這是由于液滴間相互作用所導(dǎo)致的流動變化主要發(fā)生在液滴通道一側(cè),對流動中期液滴非通道側(cè)的流場幾乎沒有改變.因此,到此時為止液滴非通道側(cè)的半徑改變與單個液滴相比幾乎是一樣的.另一方面,通道側(cè)半徑的增加比非通道側(cè)更慢.這是因為在通道側(cè)由于兩個液滴的弓形反射激波的相交,導(dǎo)致通道側(cè)相比非通道側(cè)出現(xiàn)了更高的壓強(qiáng)場(這一點(diǎn)在上一節(jié)波系演化分析中已提及),從而抑制了通道側(cè)半徑的增長.對于不同液滴間距下液滴通道側(cè)的半徑演化情況,下面將進(jìn)一步進(jìn)行探究.

    圖6 并排液滴兩側(cè)半徑(L/R=1)及單液滴半徑隨時間的變化Fig.6 The radius evolution of side-by-side droplets (L/R=1) and single droplet

    在本文中,由于計算資源的限制,本文僅考慮直至沖擊中期液滴界面的演化過程,沒有進(jìn)一步研究沖擊后期的液滴破碎過程.另一方面,由于液滴中期的界面形態(tài)及流場狀態(tài)對液滴后續(xù)氣動破碎的過程至關(guān)重要,因此可以預(yù)見,通道側(cè)與非通道側(cè)的流場與界面差異必然導(dǎo)致液滴兩側(cè)破碎的時間與剝離出子液滴的快慢出現(xiàn)不同.

    2.3 通道間距對液滴相互作用的影響

    圖7 給出了液滴距離與初始半徑比(L/R)分別為1,0.5 以及0.2 在不同時刻X-Y平面(Z=0)上的液滴輪廓圖.這里液滴與激波相互作用問題的流動參數(shù)與之前數(shù)值模擬是一致的,僅僅改變了并排液滴之間的距離.通過觀察圖7 可以發(fā)現(xiàn),在不同通道間距下,液滴通道側(cè)與非通道側(cè)的界面變形都是不對稱的.此外,通道間距的改變對非通道側(cè)的界面形態(tài)的演化影響很小,例如非通道側(cè)的界面特征結(jié)構(gòu)的相對位置基本不變,而且界面變形的大小也基本相同.但對于通道側(cè)來說,液滴間初始距離的改變對液滴界面形態(tài)的影響是非常明顯的,特別是赤道附近以及背風(fēng)面界面形態(tài)的演化.首先,當(dāng)非通道側(cè)赤道處界面出現(xiàn)波動時,通道側(cè)的波動則更為明顯.其次,當(dāng)并排液滴距離減小時,同一時刻通道側(cè)的界面波動幅值隨之發(fā)生顯著增加.以L/R=0.2 的問題為例(如圖7(c)所示),可以觀察到液滴赤道(X=0)界面出現(xiàn)了非常明顯的凹陷.這是由于隨著通道初始距離變小,上游來流在經(jīng)過通道時被阻塞的程度增大,進(jìn)而導(dǎo)致通道內(nèi)氣流壓力升高,從而擠壓通道側(cè)赤道附近界面向兩側(cè)運(yùn)動.

    圖7 不同液滴間距下X-Y 平面(Z=0)液滴的截面圖Fig.7 The cross-sectional profile of droplets in the X-Y plane (Z=0)with different droplet gaps

    另一方面,液滴間通道下游出口處的氣流膨脹,會誘導(dǎo)該處并排液滴之間的界面相互靠近.圖8 給出了在L/R=0.5 和0.2 通道間距下并排液滴在Z=0的X-Y截面內(nèi)的壓力云圖和界面形狀.可以看到在兩個通道間距下均發(fā)生通道下游出口處界面靠近乃至閉合的情形,其中L/R=0.2 構(gòu)型大約在16 μs 發(fā)生局部界面閉合.界面閉合所導(dǎo)致的氣流阻塞顯著升高了通道內(nèi)的壓力.此外在圖8(b1)~圖8(c2)可以看到在液滴內(nèi)部出現(xiàn)了高壓區(qū).其形成機(jī)制類似于激波沖擊下單液滴內(nèi)部的高壓形成[33],都是由于液滴內(nèi)部的壓縮波與彎曲界面作用而匯聚形成的.不同之處在于,并排液滴的界面不對稱性導(dǎo)致了高壓區(qū)在液滴內(nèi)部的位置和形狀都發(fā)生了一定的變化.圖9給出了兩液滴之間通道中心點(diǎn)(即圖7中的坐標(biāo)原點(diǎn))的壓力隨時間的變化.早期(0~5 μs)通道中心點(diǎn)的壓力主要受到初始激波IS 沖擊的作用先增大,并在弓形反射激波RS 的作用下減弱.隨后,通道出口的逐漸減小,逐步阻塞通道內(nèi)氣體流動,使得通道內(nèi)的壓力會逐漸上升.當(dāng)局部界面閉合發(fā)生時,通道內(nèi)的壓力達(dá)到其峰值.可以發(fā)現(xiàn)通道內(nèi)壓力峰值的出現(xiàn)與通道間距大小相關(guān)聯(lián),其中最小間距(L/R=0.2)的峰值最大,且峰值出現(xiàn)的時間也最早.由于界面閉合處承受著巨大的流向壓差,導(dǎo)致其后續(xù)發(fā)生斷裂破碎,大大減弱了對氣流的阻塞作用,從而使得通道內(nèi)的壓力逐漸降低.

    圖8 不同通道間距下,并排液滴的壓力云圖.標(biāo)號為2 的列為L/R=0.2Fig.8 Pressure contours of side-by-side droplets under different channel spacing.The column labeled 1 corresponds to L/R=0.5 and the column labeled 2 corresponds to L/R=0.2

    圖9 不同L/R 下通道中心(X=Y=Z=0)處壓力值隨時間變化Fig.9 Time variation of pressure at the channel center (X=Y=Z=0)with different L/R

    圖10 和圖11 分別給出了在L/R=0.2 構(gòu)型下在19 μs 和22 μs 兩個時刻的液滴界面三維形態(tài).可以看到在19 μs 時,兩個液滴在赤道后方已經(jīng)通過凸起的界面連接在一起,形成了一個類似楔形面的結(jié)構(gòu).在22 μs 時,這個結(jié)構(gòu)就已經(jīng)發(fā)生了斷裂破碎,生了許多脫離的條狀液絲和液滴,并且脫離的液絲與兩個液滴原來連接位置的連線方向大致平行.這些現(xiàn)象初步展現(xiàn)了并排液滴間的相互作用,特別是通道側(cè)與非通道側(cè)在氣動破碎過程上的差異性.通過以上流動分析,通道內(nèi)氣流和兩側(cè)界面演化的耦合是這些復(fù)雜液滴相互作用出現(xiàn)的直接原因.

    圖10 t=19 μs 時刻的液滴的三維界面形態(tài)(L/R=0.2)Fig.10 Snapshot of the two droplets with L/R=0.2 at t=19 μs

    圖11 t=22 μs 時刻的液滴的三維界面形態(tài)(L/R=0.2)Fig.11 Snapshot of the two droplets with L/R=0.2 at t=22 μs

    3 結(jié)論

    本文采用三維守恒清晰界面數(shù)值方法數(shù)值研究了平面激波沖擊下并排液滴的動力學(xué)過程并分析了演化初期液滴外波系結(jié)構(gòu)演化過程.發(fā)現(xiàn)沖擊初期(0~ 5 μs)并排液滴之間通道側(cè)與非通道側(cè)的波系結(jié)構(gòu)發(fā)展截然不同.在液滴通道側(cè)由于兩個液滴的弓形反射激波相交形成了新的反射激波以及馬赫桿,而且入射激波在通道側(cè)連接了并排液滴馬赫桿,因此與非通道側(cè)反射激波形成的彎曲波陣面區(qū)別明顯.沖擊初期雖然入射激波已掃過并排液滴,但由于慣性效應(yīng),液滴形狀仍然基本保持不變.給出了沖擊中期(7~ 16 μs)并排液滴三維界面形態(tài)的演化圖像,研究了入射激波遠(yuǎn)離并排液滴時的界面形態(tài)演化規(guī)律,揭示了通道下游出口處新的流動現(xiàn)象,包括氣流膨脹導(dǎo)致的界面逐漸接近至閉合及后續(xù)氣流阻塞導(dǎo)致的界面破碎等.最后,研究了并排液滴間距對液滴相互作用的影響,發(fā)現(xiàn)液滴間距大小對通道內(nèi)壓力峰值的出現(xiàn)影響顯著,如壓力峰值的出現(xiàn)與通道下游出口閉合是同步發(fā)生的.另外,更小液滴間距不僅帶來更大壓力峰值,也使得峰值出現(xiàn)的時間更早.這些液滴間耦合作用機(jī)制的發(fā)現(xiàn)對準(zhǔn)確建立高速氣流中稠密兩相流計算的相間作用模型有一定的借鑒意義.

    猜你喜歡
    激波液滴氣流
    氣流的威力
    一種基于聚類分析的二維激波模式識別算法
    基于HIFiRE-2超燃發(fā)動機(jī)內(nèi)流道的激波邊界層干擾分析
    液滴間相互碰撞融合與破碎的實驗研究
    噴淋液滴在空氣環(huán)境下的運(yùn)動特性
    斜激波入射V形鈍前緣溢流口激波干擾研究
    適于可壓縮多尺度流動的緊致型激波捕捉格式
    固體運(yùn)載火箭變軌發(fā)動機(jī)噴管氣流分離研究
    飛片下的空氣形成的“超強(qiáng)高速氣流刀”
    基于停留時間分布的氣流床氣化爐通用網(wǎng)絡(luò)模型
    国产私拍福利视频在线观看| 午夜精品久久久久久毛片777| 亚洲第一电影网av| 小蜜桃在线观看免费完整版高清| .国产精品久久| 免费看a级黄色片| 国产精品日韩av在线免费观看| 俺也久久电影网| 成人亚洲精品av一区二区| 欧美日本视频| 国内少妇人妻偷人精品xxx网站| 国产精品,欧美在线| 精品久久久久久久末码| 国产精品久久久久久久久免| 人人妻人人看人人澡| 欧美日韩黄片免| 亚洲综合色惰| 99久久无色码亚洲精品果冻| 久久草成人影院| 国产亚洲精品av在线| 欧美一区二区精品小视频在线| 亚洲av免费在线观看| 午夜福利在线观看吧| 精品99又大又爽又粗少妇毛片 | 成人午夜高清在线视频| 啦啦啦啦在线视频资源| 欧美高清成人免费视频www| 国内少妇人妻偷人精品xxx网站| 国产亚洲91精品色在线| 在线观看午夜福利视频| 欧美国产日韩亚洲一区| h日本视频在线播放| av国产免费在线观看| 极品教师在线免费播放| 狂野欧美白嫩少妇大欣赏| 窝窝影院91人妻| 国产高清视频在线播放一区| 哪里可以看免费的av片| 国内精品宾馆在线| 老师上课跳d突然被开到最大视频| 亚洲va日本ⅴa欧美va伊人久久| 亚洲自拍偷在线| 久久久精品欧美日韩精品| 欧美xxxx性猛交bbbb| 窝窝影院91人妻| 国产高潮美女av| 十八禁网站免费在线| 国内精品宾馆在线| 欧美一级a爱片免费观看看| 白带黄色成豆腐渣| av在线蜜桃| 久久精品国产自在天天线| 舔av片在线| 在线免费观看不下载黄p国产 | 成人午夜高清在线视频| 精品午夜福利视频在线观看一区| 熟妇人妻久久中文字幕3abv| 日本五十路高清| 乱人视频在线观看| 桃色一区二区三区在线观看| av天堂中文字幕网| 色综合亚洲欧美另类图片| 久久午夜福利片| 精品日产1卡2卡| 久久久久久久久中文| 俄罗斯特黄特色一大片| 两个人视频免费观看高清| 桃红色精品国产亚洲av| 成熟少妇高潮喷水视频| 一个人看视频在线观看www免费| av.在线天堂| 国内揄拍国产精品人妻在线| 精品久久久久久久久av| 深爱激情五月婷婷| 亚洲av熟女| 男插女下体视频免费在线播放| 亚洲成a人片在线一区二区| 日本在线视频免费播放| 免费观看精品视频网站| 亚洲午夜理论影院| av在线蜜桃| 他把我摸到了高潮在线观看| 少妇裸体淫交视频免费看高清| 成人性生交大片免费视频hd| 久久人妻av系列| 老司机深夜福利视频在线观看| 大型黄色视频在线免费观看| 国产91精品成人一区二区三区| 国产精品精品国产色婷婷| 两个人视频免费观看高清| 一级黄片播放器| 中文字幕久久专区| 又紧又爽又黄一区二区| 自拍偷自拍亚洲精品老妇| 十八禁网站免费在线| 男女视频在线观看网站免费| 亚洲四区av| 国产亚洲欧美98| 大又大粗又爽又黄少妇毛片口| 日本一本二区三区精品| 精品日产1卡2卡| 九色国产91popny在线| 国产伦在线观看视频一区| 女同久久另类99精品国产91| 日韩欧美国产一区二区入口| 97人妻精品一区二区三区麻豆| 97超视频在线观看视频| 亚洲内射少妇av| av在线天堂中文字幕| 天天躁日日操中文字幕| 久久精品国产亚洲网站| 又黄又爽又刺激的免费视频.| 久久久久国产精品人妻aⅴ院| 99久久久亚洲精品蜜臀av| 日韩欧美精品v在线| 99久久中文字幕三级久久日本| 欧美性猛交黑人性爽| 久久6这里有精品| 韩国av在线不卡| 亚洲久久久久久中文字幕| 国产av一区在线观看免费| 搡老熟女国产l中国老女人| 在线播放国产精品三级| 级片在线观看| 亚洲三级黄色毛片| 中文资源天堂在线| 国产精品人妻久久久久久| 真实男女啪啪啪动态图| 免费av观看视频| 变态另类成人亚洲欧美熟女| 97碰自拍视频| 他把我摸到了高潮在线观看| 久久中文看片网| 国产黄a三级三级三级人| av在线天堂中文字幕| 亚洲欧美日韩高清专用| 小蜜桃在线观看免费完整版高清| 深夜a级毛片| 男女视频在线观看网站免费| 精品一区二区三区av网在线观看| 男女啪啪激烈高潮av片| 校园人妻丝袜中文字幕| 日本三级黄在线观看| 夜夜看夜夜爽夜夜摸| 国产精品三级大全| 久久精品国产99精品国产亚洲性色| 网址你懂的国产日韩在线| 精品乱码久久久久久99久播| 人妻少妇偷人精品九色| 中文资源天堂在线| 99久久精品热视频| 精品一区二区三区视频在线| 中国美白少妇内射xxxbb| eeuss影院久久| 一本久久中文字幕| 国产伦精品一区二区三区四那| 国产三级在线视频| 99久久无色码亚洲精品果冻| 亚洲av成人精品一区久久| 男女啪啪激烈高潮av片| av视频在线观看入口| 久久国产精品人妻蜜桃| 成人国产一区最新在线观看| 天美传媒精品一区二区| 久久精品国产亚洲av天美| 午夜福利在线观看免费完整高清在 | 国产精品伦人一区二区| 小说图片视频综合网站| 色哟哟·www| 最新在线观看一区二区三区| 国产精品1区2区在线观看.| 免费不卡的大黄色大毛片视频在线观看 | 午夜激情欧美在线| 国内精品一区二区在线观看| 亚洲一级一片aⅴ在线观看| 18禁黄网站禁片午夜丰满| 日韩国内少妇激情av| 国产精品一区二区免费欧美| 国产精品综合久久久久久久免费| 国产黄色小视频在线观看| 高清在线国产一区| 成人国产一区最新在线观看| а√天堂www在线а√下载| 国产av一区在线观看免费| 午夜福利高清视频| 女人十人毛片免费观看3o分钟| 一本精品99久久精品77| 国产爱豆传媒在线观看| 制服丝袜大香蕉在线| 乱码一卡2卡4卡精品| 欧美+亚洲+日韩+国产| 91久久精品国产一区二区成人| 亚洲在线自拍视频| 国产 一区精品| 桃红色精品国产亚洲av| 搡老熟女国产l中国老女人| 91av网一区二区| 欧美激情在线99| 又爽又黄无遮挡网站| 成人综合一区亚洲| 白带黄色成豆腐渣| 成人三级黄色视频| 狂野欧美激情性xxxx在线观看| 熟女人妻精品中文字幕| a级一级毛片免费在线观看| 精品一区二区三区人妻视频| 亚洲内射少妇av| bbb黄色大片| 国产精品电影一区二区三区| 12—13女人毛片做爰片一| 欧美黑人欧美精品刺激| 久久国内精品自在自线图片| 亚洲人成网站高清观看| 九九在线视频观看精品| 国产免费av片在线观看野外av| 亚洲aⅴ乱码一区二区在线播放| 精品久久久久久成人av| 亚洲成人免费电影在线观看| 亚洲狠狠婷婷综合久久图片| 在线播放国产精品三级| 欧美又色又爽又黄视频| 久久午夜福利片| 久久久久久久久中文| 成人国产麻豆网| 热99re8久久精品国产| 黄色日韩在线| 国国产精品蜜臀av免费| 99热这里只有是精品在线观看| 成人永久免费在线观看视频| 婷婷丁香在线五月| 人妻丰满熟妇av一区二区三区| av.在线天堂| 动漫黄色视频在线观看| 国产精品爽爽va在线观看网站| av女优亚洲男人天堂| 色视频www国产| 中国美白少妇内射xxxbb| 日韩中字成人| 色精品久久人妻99蜜桃| 一本精品99久久精品77| 亚洲一区二区三区色噜噜| 久久中文看片网| 一进一出抽搐动态| 国产免费男女视频| 日本 av在线| 免费观看的影片在线观看| 日韩欧美在线二视频| 国产亚洲精品av在线| 一a级毛片在线观看| 亚洲五月天丁香| 精品日产1卡2卡| 国产精品1区2区在线观看.| www.色视频.com| 国产精华一区二区三区| 亚洲欧美日韩高清在线视频| 国产高清视频在线播放一区| a级毛片免费高清观看在线播放| 欧美日韩亚洲国产一区二区在线观看| av.在线天堂| 久久精品影院6| 可以在线观看的亚洲视频| 99热这里只有是精品50| 欧美三级亚洲精品| 久久人人爽人人爽人人片va| 乱系列少妇在线播放| 网址你懂的国产日韩在线| 欧美日本亚洲视频在线播放| 国产一区二区三区在线臀色熟女| 麻豆av噜噜一区二区三区| 欧美日韩国产亚洲二区| 少妇高潮的动态图| 九色成人免费人妻av| 最近在线观看免费完整版| 尾随美女入室| 国产乱人视频| 免费看美女性在线毛片视频| 成人无遮挡网站| 日韩在线高清观看一区二区三区 | a级毛片a级免费在线| 欧美3d第一页| 看免费成人av毛片| 少妇丰满av| 少妇猛男粗大的猛烈进出视频 | 神马国产精品三级电影在线观看| 免费观看人在逋| 亚洲欧美精品综合久久99| 欧美日韩瑟瑟在线播放| 色综合亚洲欧美另类图片| 久久久色成人| 欧美一级a爱片免费观看看| www日本黄色视频网| 极品教师在线免费播放| 99久国产av精品| 国产高清不卡午夜福利| 波多野结衣高清作品| 国产亚洲欧美98| 人妻夜夜爽99麻豆av| 偷拍熟女少妇极品色| 精品久久久久久成人av| 色5月婷婷丁香| 免费看日本二区| 亚洲av中文字字幕乱码综合| 2021天堂中文幕一二区在线观| 我要看日韩黄色一级片| 欧洲精品卡2卡3卡4卡5卡区| 欧美日韩黄片免| 成年女人永久免费观看视频| 夜夜夜夜夜久久久久| 不卡视频在线观看欧美| 少妇人妻一区二区三区视频| 国产高清有码在线观看视频| 国产伦精品一区二区三区四那| 久久这里只有精品中国| 日本一二三区视频观看| 成人精品一区二区免费| 精品人妻偷拍中文字幕| 亚洲国产精品成人综合色| 亚洲精品乱码久久久v下载方式| 内射极品少妇av片p| 女生性感内裤真人,穿戴方法视频| 亚洲人成网站在线播| 亚洲成人免费电影在线观看| 欧美性猛交╳xxx乱大交人| 老司机福利观看| 女同久久另类99精品国产91| 亚洲狠狠婷婷综合久久图片| av黄色大香蕉| 欧美日韩黄片免| 亚洲av一区综合| 中文字幕av在线有码专区| 91在线观看av| 99热网站在线观看| 国产亚洲欧美98| 精品福利观看| 午夜福利欧美成人| 99热这里只有是精品在线观看| 亚洲午夜理论影院| 久久精品夜夜夜夜夜久久蜜豆| 国产麻豆成人av免费视频| 国产伦在线观看视频一区| 热99在线观看视频| 亚洲最大成人手机在线| 成人美女网站在线观看视频| 亚洲乱码一区二区免费版| 狠狠狠狠99中文字幕| 亚洲欧美日韩高清在线视频| 少妇人妻一区二区三区视频| 午夜福利视频1000在线观看| 国产成人一区二区在线| 日本精品一区二区三区蜜桃| 麻豆av噜噜一区二区三区| 成人国产一区最新在线观看| 日韩 亚洲 欧美在线| 老熟妇仑乱视频hdxx| 精品一区二区三区视频在线| 1000部很黄的大片| 亚洲国产色片| 人妻夜夜爽99麻豆av| 精品欧美国产一区二区三| av.在线天堂| 午夜福利视频1000在线观看| 久久香蕉精品热| 国产成人影院久久av| 小说图片视频综合网站| 国产成人影院久久av| 黄色一级大片看看| 久久久久久久午夜电影| 国产精品久久久久久亚洲av鲁大| 亚洲熟妇熟女久久| 免费看日本二区| 国产极品精品免费视频能看的| 国产精品,欧美在线| 久久久久九九精品影院| 国产激情偷乱视频一区二区| 亚洲真实伦在线观看| 深夜精品福利| 桃红色精品国产亚洲av| 国产精品一区二区免费欧美| 久久久久久九九精品二区国产| 男女边吃奶边做爰视频| 国产视频内射| av福利片在线观看| 露出奶头的视频| 亚洲精品乱码久久久v下载方式| 两个人的视频大全免费| 久久6这里有精品| 免费一级毛片在线播放高清视频| 久久精品国产亚洲av香蕉五月| 99热6这里只有精品| av.在线天堂| 久久6这里有精品| 国产精品不卡视频一区二区| 丰满人妻一区二区三区视频av| x7x7x7水蜜桃| 国产午夜精品论理片| 国产精品久久久久久亚洲av鲁大| 亚洲国产欧洲综合997久久,| 免费搜索国产男女视频| 亚洲人成伊人成综合网2020| 日日摸夜夜添夜夜添av毛片 | 91久久精品国产一区二区三区| 老师上课跳d突然被开到最大视频| 午夜福利高清视频| 嫩草影院入口| 欧美+日韩+精品| 国产精品人妻久久久影院| 成人亚洲精品av一区二区| 日韩精品青青久久久久久| 男插女下体视频免费在线播放| 日韩,欧美,国产一区二区三区 | 久久香蕉精品热| 成年人黄色毛片网站| 国产女主播在线喷水免费视频网站 | 免费观看的影片在线观看| 国产成人a区在线观看| 亚洲成人精品中文字幕电影| 国产精华一区二区三区| 欧美+日韩+精品| 亚洲自偷自拍三级| 51国产日韩欧美| 免费av观看视频| www.色视频.com| 美女高潮的动态| 亚洲色图av天堂| 国产黄a三级三级三级人| 很黄的视频免费| 伊人久久精品亚洲午夜| www.色视频.com| 精品人妻熟女av久视频| 国产伦在线观看视频一区| 久久精品国产亚洲av香蕉五月| 午夜福利视频1000在线观看| 天天一区二区日本电影三级| 有码 亚洲区| 欧美成人免费av一区二区三区| 日本成人三级电影网站| 国模一区二区三区四区视频| 亚洲精华国产精华液的使用体验 | 亚洲aⅴ乱码一区二区在线播放| 国产久久久一区二区三区| 免费无遮挡裸体视频| 成年女人永久免费观看视频| 国产69精品久久久久777片| 婷婷亚洲欧美| 久久人人爽人人爽人人片va| 亚洲av电影不卡..在线观看| 久久热精品热| 给我免费播放毛片高清在线观看| 欧美激情国产日韩精品一区| 麻豆国产97在线/欧美| 天堂影院成人在线观看| 无人区码免费观看不卡| 亚洲无线观看免费| 国产精品三级大全| 日本一本二区三区精品| 丰满的人妻完整版| 如何舔出高潮| 国产主播在线观看一区二区| 成人综合一区亚洲| 偷拍熟女少妇极品色| 老熟妇乱子伦视频在线观看| 日韩亚洲欧美综合| 国产午夜精品久久久久久一区二区三区 | 狠狠狠狠99中文字幕| 国产精品亚洲一级av第二区| 国产精品一区www在线观看 | 观看免费一级毛片| 亚洲最大成人av| 国产免费男女视频| 亚洲va在线va天堂va国产| 12—13女人毛片做爰片一| 久久久久性生活片| 搡女人真爽免费视频火全软件 | 内射极品少妇av片p| 欧美成人一区二区免费高清观看| 亚洲成人久久性| 美女高潮的动态| 黄色女人牲交| 欧美国产日韩亚洲一区| 一区福利在线观看| 亚洲精品乱码久久久v下载方式| 午夜a级毛片| 中文字幕人妻熟人妻熟丝袜美| АⅤ资源中文在线天堂| 不卡视频在线观看欧美| 禁无遮挡网站| 亚洲乱码一区二区免费版| 国产精品永久免费网站| 有码 亚洲区| 天天一区二区日本电影三级| 亚洲四区av| 午夜影院日韩av| 99久久成人亚洲精品观看| 国产精品三级大全| 日韩精品中文字幕看吧| 日本一二三区视频观看| 午夜影院日韩av| 精品久久久久久久久久久久久| 变态另类丝袜制服| 国产精品乱码一区二三区的特点| 国产精华一区二区三区| 亚洲美女黄片视频| 老女人水多毛片| 久久久精品欧美日韩精品| 国产精品98久久久久久宅男小说| 一进一出抽搐动态| 久久这里只有精品中国| 真人做人爱边吃奶动态| 欧美黑人巨大hd| 22中文网久久字幕| 欧美成人免费av一区二区三区| 亚洲自偷自拍三级| 两性午夜刺激爽爽歪歪视频在线观看| 日韩欧美精品免费久久| 人人妻,人人澡人人爽秒播| 欧美国产日韩亚洲一区| 91久久精品国产一区二区三区| 看十八女毛片水多多多| 18禁黄网站禁片午夜丰满| 久久精品久久久久久噜噜老黄 | 日本a在线网址| 国产精品久久久久久久电影| 亚洲欧美激情综合另类| 久久婷婷人人爽人人干人人爱| av黄色大香蕉| 亚洲18禁久久av| 午夜视频国产福利| 欧美成人a在线观看| 日韩人妻高清精品专区| 色av中文字幕| 精品久久国产蜜桃| 国产黄片美女视频| 久久九九热精品免费| 免费在线观看影片大全网站| 无人区码免费观看不卡| 日韩欧美一区二区三区在线观看| 午夜视频国产福利| 免费大片18禁| 白带黄色成豆腐渣| 久久午夜亚洲精品久久| 97人妻精品一区二区三区麻豆| 亚洲在线自拍视频| 狠狠狠狠99中文字幕| 悠悠久久av| 欧美极品一区二区三区四区| 国产一区二区激情短视频| 国产乱人伦免费视频| 无人区码免费观看不卡| 日日撸夜夜添| 伦精品一区二区三区| 久久九九热精品免费| 美女大奶头视频| 免费在线观看日本一区| av视频在线观看入口| 一进一出抽搐gif免费好疼| 亚洲男人的天堂狠狠| 精品久久国产蜜桃| 免费av不卡在线播放| 亚洲欧美激情综合另类| 国产精品国产三级国产av玫瑰| 欧美+日韩+精品| 一本一本综合久久| 黄片wwwwww| 麻豆一二三区av精品| 欧美xxxx黑人xx丫x性爽| 免费看a级黄色片| 精品久久久久久久久久免费视频| 观看美女的网站| 三级国产精品欧美在线观看| 精品99又大又爽又粗少妇毛片 | 亚洲av五月六月丁香网| 51国产日韩欧美| 日日撸夜夜添| 内射极品少妇av片p| 亚洲国产色片| 久久99热这里只有精品18| 日日摸夜夜添夜夜添av毛片 | 国产三级在线视频| 淫妇啪啪啪对白视频| 欧美成人性av电影在线观看| 国产爱豆传媒在线观看| 又爽又黄a免费视频| 少妇裸体淫交视频免费看高清| av在线蜜桃| 亚洲国产精品久久男人天堂| av在线天堂中文字幕| 不卡一级毛片| 九九爱精品视频在线观看| 乱人视频在线观看| 色哟哟哟哟哟哟| 亚洲精品乱码久久久v下载方式| 国产精品自产拍在线观看55亚洲| 久久精品91蜜桃| 一卡2卡三卡四卡精品乱码亚洲| 欧美一区二区国产精品久久精品| 久久香蕉精品热| 如何舔出高潮| 91久久精品国产一区二区三区| 日韩中文字幕欧美一区二区| 精品日产1卡2卡| 亚洲最大成人av| 日本黄大片高清| 欧美国产日韩亚洲一区| 制服丝袜大香蕉在线| 午夜福利高清视频| 国产又黄又爽又无遮挡在线| 亚洲欧美日韩卡通动漫| 午夜福利在线在线| 日本一本二区三区精品| 免费搜索国产男女视频| 欧美日本亚洲视频在线播放| 国产精品电影一区二区三区| 日韩高清综合在线| 日本色播在线视频| 久久久久久久久中文|