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

    反射激波作用下三維凹?xì)庵缑嫜莼臄?shù)值研究1)

    2021-05-31 01:34:08崔竹軒丁舉春
    力學(xué)學(xué)報(bào) 2021年5期
    關(guān)鍵詞:氣柱渦量不穩(wěn)定性

    崔竹軒 丁舉春 司 廷

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

    引言

    Richtmyer-Meshkov (RM) 不穩(wěn)定性是指一個(gè)具有初始擾動(dòng)的流體分界面受到瞬時(shí)沖擊力(比如激波)作用后,界面上的初始小擾動(dòng)隨著時(shí)間不斷增長(zhǎng),在演化后期界面附近會(huì)生成許多小尺度渦,并最終形成湍流混合的現(xiàn)象.在1960 年,Richmyer[1]提出了激波加速下流體流動(dòng)穩(wěn)定性的理論推導(dǎo)和分析預(yù)測(cè),Meshkov[2]在1969 年對(duì)這一問(wèn)題開(kāi)展了一系列激波管實(shí)驗(yàn),證實(shí)了Richtmyer 的理論預(yù)測(cè).開(kāi)展激波-界面相互作用研究不僅在流動(dòng)不穩(wěn)定性、旋渦動(dòng)力學(xué)以及湍流形成機(jī)理等方面有著重要的學(xué)術(shù)價(jià)值,而且在工程領(lǐng)域也有廣泛的應(yīng)用.例如在慣性約束核聚變[3-5]中,氘氚靶丸外層為燒蝕材料,內(nèi)層為氘氚燃料.當(dāng)高能激光照射球形靶丸,表面燒蝕層材料被高溫?zé)g后形成向中心傳播的匯聚激波.受加工精度限制,靶丸各層分界面存在初始小擾動(dòng),這些小擾動(dòng)在激波作用后發(fā)生RM 不穩(wěn)定性演化,從而加劇了燃料和外層殼體材料的混合,使得中心區(qū)壓力和溫度下降,導(dǎo)致點(diǎn)火失敗[6].在超燃沖壓發(fā)動(dòng)機(jī)中[7],燃料和氧化劑在燃燒室的滯留時(shí)間在微秒量級(jí),激波和燃料/氧化劑界面相互作用可以促進(jìn)燃料與氧化劑的混合,提高燃燒和推進(jìn)效率[8].另外,RM不穩(wěn)定性也被用來(lái)解釋超新星爆發(fā)等天文現(xiàn)象[9].

    由于RM 不穩(wěn)定性在諸多科學(xué)和工程領(lǐng)域有重要價(jià)值,自從其被發(fā)現(xiàn)以來(lái),科研工作者已經(jīng)采用理論分析、數(shù)值計(jì)算、實(shí)驗(yàn)研究等手段開(kāi)展了大量研究[10-12].作為RM 不穩(wěn)定性的經(jīng)典案例,平面激波與氣柱相互作用已被大量研究.Ding 等[13-14]利用肥皂膜技術(shù)生成三維柱形界面,研究了三維性對(duì)重氣柱和輕氣柱的影響,并將Haas 和Sturtevant[15]提出的理論模型拓展到三維情形.Ou 等[16]實(shí)驗(yàn)和數(shù)值研究了激波沖擊橢圓重氣柱的不穩(wěn)定性發(fā)展,重點(diǎn)關(guān)注了激波聚焦和射流現(xiàn)象,并分析Atwood 數(shù)(At=(ρ2-ρ1)/(ρ2+ρ1),其中ρ1和ρ2分別為激波沖擊前氣柱外和氣柱內(nèi)氣體密度) 和橢圓氣柱的長(zhǎng)寬比[17]對(duì)氣柱演化的影響.Zou 等[18]利用粒子圖像測(cè)速(PIV) 方法研究了激波沖擊不同長(zhǎng)寬比橢圓氣柱后的速度場(chǎng),分析了長(zhǎng)寬比為0.25~0.4 的橢圓氣柱,得到了隨著長(zhǎng)寬比增大,渦量和環(huán)量的最大值增大、渦心距縮短、氣柱演化變快的結(jié)論.

    激波沖擊一道擾動(dòng)界面后,界面會(huì)持續(xù)變形,同時(shí)激波也會(huì)變成一道擾動(dòng)激波,擾動(dòng)激波到達(dá)反射壁面后形成一道反射激波,二次沖擊處于演化中的流體界面.因此,考察反射激波和變形界面相互作用的不穩(wěn)定性演化至關(guān)重要[19].二次作用之前,激波和界面的形狀會(huì)影響斜壓渦量的生成和壓力分布,從而影響二次沖擊之后的不穩(wěn)定性演化.反射激波二次沖擊界面的不穩(wěn)定性發(fā)展受到反射距離、激波強(qiáng)度等諸多參數(shù)的影響,因此反射激波對(duì)界面演化造成的影響比單次激波沖擊情形要更加復(fù)雜,同時(shí),研究反射激波與界面相作用對(duì)于研究湍流混合有重要價(jià)值[20-23].

    目前,反射激波誘導(dǎo)RM 不穩(wěn)定性的研究相對(duì)較少.實(shí)驗(yàn)研究一直是研究激波作用下RM 不穩(wěn)定性發(fā)展的重要手段.Jacobs 等[24]利用振蕩產(chǎn)生三維短波長(zhǎng)擾動(dòng),在豎直激波管中研究了反射激波沖擊Air/SF6界面的不穩(wěn)定性發(fā)展,獲得了混合寬度增長(zhǎng)的定量數(shù)據(jù),發(fā)現(xiàn)反射激波作用后,不同初始條件下的混合寬度以相同增長(zhǎng)率線性增長(zhǎng).Si 等[25]進(jìn)行了不同反射距離27~125 mm 下反射激波沖擊Air/SF6氣泡的實(shí)驗(yàn)研究,發(fā)現(xiàn)了反射距離較小時(shí)會(huì)出現(xiàn)方向相反的射流以及界面振蕩現(xiàn)象,反射距離較長(zhǎng)時(shí)則會(huì)出現(xiàn)小尺度的渦結(jié)構(gòu).Mohaghar 等[26]則在傾斜激波管中,利用平面激光誘導(dǎo)熒光(PLIF)技術(shù),獲得了單模及多模界面的濃度演化圖.此外,已有的理論研究多關(guān)注適用于二維單模界面混合寬度增長(zhǎng)率的模型,例如,Brouillette 和Strurtevant[27]給出的廣義Richtmyer 模型.鑒于RM 不穩(wěn)定性的實(shí)驗(yàn)研究對(duì)測(cè)量技術(shù)要求很高且理論分析很難獲得大擾動(dòng)界面演化的解析解,數(shù)值模擬便成為一種廣泛采用的有效研究方法[28-30].數(shù)值計(jì)算方法近年來(lái)發(fā)展十分迅速.Hill 等[20]利用大渦模擬,對(duì)不同馬赫數(shù)下反射激波沖擊平直Air/SF6界面進(jìn)行模擬,得到反射激波作用后,混合層寬度增長(zhǎng)率隨激波馬赫數(shù)增大而變大的規(guī)律.Shankar 和Lele[31]通過(guò)數(shù)值模擬研究了反射激波沖擊Air/SF6氣簾界面,得到了反射距離越大,反射激波作用后界面會(huì)以更快的速度混合并產(chǎn)生小尺寸結(jié)構(gòu)的結(jié)論.

    在已開(kāi)展的反射激波沖擊二維氣柱研究中,實(shí)驗(yàn)和數(shù)值計(jì)算均取得了一定進(jìn)展.實(shí)驗(yàn)方面,Zhang等[32]利用PIV 觀測(cè)反射激波與SF6無(wú)膜重氣柱相互作用后的流場(chǎng),獲得了相應(yīng)的速度場(chǎng)和渦量場(chǎng),He 等[33]利用無(wú)膜氣柱研究了反射激波作用SF6和Ar 兩種重氣柱的演化,分析了氣柱幾何尺寸的變化規(guī)律.數(shù)值計(jì)算方面,Wang 等[34]使用VAS2D(2-Dimensional &Axisymmetric Vectorized Adaptive Solver) 程序研究了不同反射距離下Air/SF6氣柱的演化,分析了界面幾何尺寸,特征點(diǎn)以及環(huán)量的變化規(guī)律,Sha 等[35]基于大渦模擬研究了反射激波作用Air/SF6氣柱后的波系以及射流的演化.這些研究主要集中在二維重氣柱方向,三維性和Atwood 數(shù)對(duì)氣柱界面演化的影響還有待更深入研究.

    雖然反射激波作用氣柱的RM 不穩(wěn)定性已有一些研究,但反射激波作用輕氣柱以及反射激波作用三維氣柱的研究尚未見(jiàn)報(bào)道,其中涉及的物理規(guī)律和機(jī)理仍有待揭示.鑒于上述討論,本文采用數(shù)值計(jì)算方法,以三維凹?xì)庵鶠檠芯繉?duì)象,選取反射距離以及凹?xì)庵娜S性為變量,對(duì)比研究不同反射距離時(shí)二維和三維凹?xì)庵缑娴难莼^(guò)程,從而進(jìn)一步揭示反射激波與氣柱界面相互作用的RM 不穩(wěn)定性演化規(guī)律和機(jī)理.

    1 數(shù)值計(jì)算方法

    1.1 控制方程與數(shù)值程序

    激波與界面作用時(shí)間很短,黏性、重力、氣體擴(kuò)散以及熱傳導(dǎo)等因素對(duì)氣柱演化影響很小,同時(shí)忽略化學(xué)反應(yīng)過(guò)程,這類(lèi)流動(dòng)可以用可壓縮無(wú)黏歐拉方程來(lái)表示

    方程中U表示守恒變量,F,G,H分別表示x,y,z方向上的通量.可以表示為

    其中,u,v,w表示流體在x,y,z各方向上的速度,ρ 表示混合氣體的密度,ρi表示界面內(nèi)單種氣體的密度,p表示混合氣體的壓力,E表示流體單位體積所包含的總能量.為了使得方程組封閉,還需增加氣體的狀態(tài)方程

    其中,Cpi表示定壓比熱容,Cvi表示定容比熱容.

    本文采用課題組自主開(kāi)發(fā)的HOWD (high-order WENO and double-flux methods) 程序[36]進(jìn)行數(shù)值模擬,該程序使用均勻網(wǎng)格,能夠?qū)崿F(xiàn)五階WENO(weighted essentially non-oscillatory)[37]重構(gòu)與雙通量算法的耦合,可以消除界面處的非物理震蕩,在準(zhǔn)確模擬激波的同時(shí),還可以捕捉流場(chǎng)中的小渦結(jié)構(gòu).由于小尺度的渦結(jié)構(gòu)在RM 不穩(wěn)定性的研究中十分重要,對(duì)數(shù)值格式提出了更高的要求,而已有的數(shù)值結(jié)果證明,WENO 格式對(duì)這些問(wèn)題有更強(qiáng)的模擬能力[21,37],因此,該程序選用WENO 格式來(lái)離散方程.

    1.2 計(jì)算模型和程序驗(yàn)證

    在已有單次激波沖擊氣柱的實(shí)驗(yàn)中,主要采取肥皂膜技術(shù)生成氣柱[38-39].由于肥皂膜存在表面張力,氣柱形狀可以用Young-Laplace 方程描述[40],即氣柱界面上各處的主曲率之和為常數(shù).在柱坐標(biāo)系中,代入主曲率計(jì)算公式可以得到

    其中,r=表示在某一高度上氣柱的半徑,rzz為半徑關(guān)于高度的二階導(dǎo)數(shù),rz為半徑關(guān)于高度的一階導(dǎo)數(shù).存在一個(gè)特殊情況即氣柱內(nèi)外壓力差相等.此時(shí)Δp=0,此時(shí)為凹?xì)庵?氣柱界面為極小曲面[41],此時(shí)氣柱界面的表達(dá)式為

    其中,r0表示z=0 時(shí)對(duì)稱(chēng)面的氣柱半徑.二維氣柱半徑取r=17.5 mm,凹?xì)庵吔缑姘霃脚c二維氣柱相同,可計(jì)算得到r0=13.7 mm.

    計(jì)算模型如圖1(a) 所示,模擬在方形截面激波管實(shí)驗(yàn)段中有一凹?xì)庵缑?實(shí)驗(yàn)段尾端封閉,即為固壁邊界.反射距離定義為氣柱圓心至右側(cè)固壁端的距離.整個(gè)計(jì)算域?yàn)?00 mm×140 mm×20 mm.為節(jié)省計(jì)算時(shí)間,對(duì)于二維算例只計(jì)算半個(gè)氣柱,計(jì)算域?yàn)?00 mm×70 mm×20 mm;對(duì)于三維算例計(jì)算對(duì)稱(chēng)的1/4 個(gè)氣柱,計(jì)算域?yàn)?00 mm×70 mm×10 mm.

    為了驗(yàn)證網(wǎng)格收斂性,以平面激波與二維氣柱作用作為參考.為了便于對(duì)比,本文計(jì)算初始條件與已開(kāi)展的實(shí)驗(yàn)[38-39]設(shè)置相同,即激波馬赫數(shù)Ma=1.29,T0=293.15 K,P0=101 325 Pa.二維工況氣柱內(nèi)氣體為97%的N2與3%的SF6混合氣體,氣柱外則為純SF6,T0=293 K,P0=101 325 Pa,凹?xì)庵鶅?nèi)氣體為95%的N2與5%的SF6混合氣體,氣柱外則為純SF6.

    如圖1(b) 所示,分別在網(wǎng)格大小為0.8 mm,0.4 mm,0.2 mm 和0.1 mm 時(shí)進(jìn)行計(jì)算,比較t=80μs時(shí)對(duì)稱(chēng)軸上的密度分布,可見(jiàn)網(wǎng)格大小為0.2 mm 和網(wǎng)格大小為0.1 mm 的計(jì)算結(jié)果已十分接近,因此二維氣柱算例采用尺寸為0.1 mm 的均勻網(wǎng)格.由于三維凹?xì)庵憷?jì)算量較大,采用尺寸為0.2 mm 的均勻網(wǎng)格.

    圖1 計(jì)算模型和收斂性驗(yàn)證Fig.1 Computational domain and grid convergence validation

    已有研究已經(jīng)對(duì)HOWD 程序的準(zhǔn)確性進(jìn)行了驗(yàn)證[13-14,16],這里采用無(wú)反射激波作用二維輕氣柱與三維凹?xì)庵缑嫜莼臄?shù)值計(jì)算結(jié)果與Ding[36]的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比.圖2 顯示的是典型時(shí)刻的紋影結(jié)果,可以觀察到數(shù)值計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果吻合,包括激波波系和界面結(jié)構(gòu)都能一一對(duì)應(yīng).證明了計(jì)算程序的可靠性和適用性.

    圖2 無(wú)反射激波情況下二維氣柱(上)和三維凹?xì)庵?下)數(shù)值計(jì)算結(jié)果與相應(yīng)的實(shí)驗(yàn)結(jié)果對(duì)比,單位為μs.對(duì)比圖中上面是實(shí)驗(yàn)紋影結(jié)果,下面是數(shù)值紋影結(jié)果Fig.2 Validation of the numerical results in comparison with experimental ones for 2D(top)and concave cylinders(bottom),respectively(time unit:μs).The experimental schlieren images are listed on top

    2 反射激波對(duì)氣柱界面演化的影響

    2.1 二維氣柱

    圖3 展示了反射距離分別為40 mm,80 mm,120 mm 時(shí)二維輕氣柱演化的數(shù)值計(jì)算結(jié)果,圖的上半部分顯示渦量分布(其中深色為正渦量,淺色為負(fù)渦量),下半部分為紋影結(jié)果.激波第一次沖擊二維氣柱界面后,氣柱界面上生成斜壓渦量,且距離對(duì)稱(chēng)面越遠(yuǎn)的界面上的渦量越大,這是因?yàn)榧げ_擊氣柱界面過(guò)程中壓力梯度和密度梯度的夾角在逐漸變化.受到激波壓縮以及斜壓渦量的影響,上游界面變得扁平,隨后從對(duì)稱(chēng)面附近開(kāi)始,上游氣柱界面產(chǎn)生向下游的凹陷并生成射流結(jié)構(gòu),隨后向下游運(yùn)動(dòng),與下游界面距離縮短形成一個(gè)橋狀結(jié)構(gòu).下游界面在600μs附近會(huì)因旋渦誘導(dǎo)[36]產(chǎn)生負(fù)渦量.

    圖3 不同反射距離(40 mm,80 mm,120 mm)下二維輕氣柱的演化圖像,并與無(wú)反射情況的結(jié)果進(jìn)行對(duì)比.入射激波從左往右運(yùn)動(dòng),最右側(cè)是固壁.圖中數(shù)字表示時(shí)刻,單位為μsFig.3 Evolution images of a two-dimensional light gas cylinder under reshock with reflected distance of 40 mm,80 mm,120 mm,respectively,in comparison with those without reshock(time unit:μs)

    當(dāng)反射激波再次作用界面后,下游界面的負(fù)渦量強(qiáng)度會(huì)更高.以反射距離40 mm 的工況為例,如圖4 所示,320μs 時(shí)反射激波剛接觸界面,密度梯度與壓力梯度幾乎重合,因此生成的斜壓渦量強(qiáng)度較小.隨著激波與界面作用位置的變化,生成的斜壓渦量強(qiáng)度逐漸變大,在360 μs 時(shí)運(yùn)動(dòng)速度較快的透射激波再次運(yùn)動(dòng)到界面外形成前導(dǎo)激波,并隨著位置的變化密度梯度與壓力梯度夾角逐漸變小,斜壓渦量強(qiáng)度變小.總體來(lái)說(shuō),斜壓渦量的強(qiáng)度沿右側(cè)界面從對(duì)稱(chēng)面到最外側(cè)先變大再變小,這與單一的入射激波作用氣柱界面的情況完全不同.

    圖4 反射距離為40 mm 時(shí),反射激波作用變形界面過(guò)程中的精細(xì)結(jié)構(gòu)變化Fig.4 Variation of fine structures during the interaction of the reshock with distorted interface at reflected distance of 40 mm

    如圖5 所示,經(jīng)過(guò)一段時(shí)間的輸運(yùn),負(fù)渦量形成一個(gè)順時(shí)針旋轉(zhuǎn)的渦,并與第一次激波沖擊形成的渦一起誘導(dǎo)氣泡形成射流狀結(jié)構(gòu),將氣泡分割成兩個(gè)部分.反射激波的作用使得右側(cè)界面上的負(fù)渦量強(qiáng)度更大,從而導(dǎo)致界面發(fā)展更快.由于反射激波與入射激波方向相反,為了描述方便,用左側(cè)(L)和右側(cè)(R)來(lái)區(qū)分方向.反射距離為40 mm 時(shí),界面仍然處在發(fā)展前期,彎曲程度較小,因此生成的渦量強(qiáng)度較弱,導(dǎo)致兩個(gè)渦的強(qiáng)度并不均等,生成的射流狀結(jié)構(gòu)在渦2 的影響下方向逐漸偏轉(zhuǎn)至指向右側(cè),渦1 本身也向右側(cè)運(yùn)動(dòng).

    圖5 氣泡被射流分割Fig.5 Bubble is cut offby jet flow

    不同反射距離下,反射激波作用界面時(shí)的激波形狀和界面形態(tài)顯著不同,因此產(chǎn)生的渦結(jié)構(gòu)也有所區(qū)別.反射距離為80 mm 時(shí),由于氣泡彎曲程度較大導(dǎo)致生成渦量強(qiáng)度較大,使得渦1 和渦2 強(qiáng)度相當(dāng),生成的射流狀結(jié)構(gòu)將界面推向左側(cè),導(dǎo)致氣柱最左側(cè)界面向左側(cè)移動(dòng).對(duì)于反射距離為120 mm 的工況,反射激波作用時(shí)氣泡已經(jīng)分裂,因此反射激波作用后的演化更為復(fù)雜,生成了數(shù)個(gè)旋轉(zhuǎn)方向各異的小渦,由于氣泡最上部界面與激波幾乎垂直,生成的負(fù)渦量較大,導(dǎo)致氣泡1 自身快速形成渦對(duì).整體上,激波二次作用會(huì)使得氣柱在演化后期有生成3 個(gè)渦對(duì)的趨勢(shì):激波第一次作用后由于正渦量輸運(yùn)形成的渦2,反射激波作用后生成的負(fù)渦量輸運(yùn)形成的渦1,以及氣泡1 在演化后期形成的渦3.隨著反射距離的變化,這些渦的強(qiáng)度和生成時(shí)間會(huì)有所不同,導(dǎo)致后期界面的發(fā)展演化也不同,反射距離40 mm 時(shí)生成的3 個(gè)渦對(duì)均清晰可見(jiàn),反射距離為80 mm 時(shí)可以看到在1600μs 尚未完成氣泡1 到渦3 的演化(如圖5(b) 所示),而反射距離120 mm 時(shí)則可以觀察到渦3 快速形成,而渦1 難以分辨(如圖3 所示).渦結(jié)構(gòu)的演化決定了變形界面的結(jié)構(gòu)形態(tài),進(jìn)而直接改變氣柱界面的演化規(guī)律.

    2.2 三維凹?xì)庵?/h3>

    三維凹?xì)庵缑嬖谌肷浼げ胺瓷浼げㄗ饔孟碌难莼^(guò)程與二維情況有顯著區(qū)別.圖6 展示了三維凹?xì)庵莼募y影圖,對(duì)應(yīng)反射距離為80 mm.由于紋影圖是由各個(gè)高度的紋影圖俯視疊加而成,因此可以觀察到界面很厚.由紋影圖像可以辨識(shí)出由邊界面生成的外層輪廓以及由對(duì)稱(chēng)面生成的內(nèi)層輪廓.隨著時(shí)間的推移,這兩層輪廓間距開(kāi)始加大,并且在斜壓渦量的影響下形成了向前的射流狀結(jié)構(gòu).與二維不同的是,這個(gè)射流最終突破最前側(cè)界面并在對(duì)稱(chēng)面附近破壞了橋狀結(jié)構(gòu),這使得氣柱在寬度方向發(fā)展的更快.反射距離為80 mm 的工況,從1000 μs就開(kāi)始發(fā)生與固壁的相互作用,這遠(yuǎn)早于二維的情況.在反射激波沖擊后,橋狀結(jié)構(gòu)斷裂開(kāi),之前的射流部分分離并擴(kuò)散成環(huán).為了便于分析,圖7 給出了反射激波作用過(guò)程中較短時(shí)間內(nèi)氣柱的三維圖像,可以看出當(dāng)激波剛穿過(guò)界面且未能恢復(fù)成平面激波時(shí),也具有很強(qiáng)的三維性,而三維凹?xì)庵缑娴难莼訌?fù)雜.

    圖6 反射距離80 mm 時(shí)凹?xì)庵莼募y影圖,單位為μsFig.6 Sequence of schlieren frames showing the evolution of a concave light gas cylinder with reflected distance of 80 mm(time unit:μs)

    圖7 反射距離80 mm 時(shí),不同時(shí)刻反射激波作用凹?xì)庵娜S圖像,單位為μsFig.7 Three-dimensional frames during the interaction of a concave light gas cylinder with the reshock at reflected distance of 80 mm(time unit:μs)

    為了更好地分析三維凹?xì)庵诜瓷浼げㄗ饔孟碌难莼^(guò)程,將凹?xì)庵衅?分層研究各自的演化情況.如圖8 所示,選取z=0,z=10 mm 兩處沿著x–y平面切片,分別記為S1 和S2 平面,選取y=0 處沿x–z方向切片,記為S3 平面.

    圖8 三維凹?xì)庵慕孛媸疽鈭DFig.8 Schematics of the cross section of 3D concave cylinder

    圖9 給出了不同反射距離下3 個(gè)典型截面上的界面演化過(guò)程.通過(guò)對(duì)比S1 和S2 兩個(gè)截面可以發(fā)現(xiàn),射流突破橋狀結(jié)構(gòu)以及后期下游界面形成的環(huán)狀結(jié)構(gòu)都發(fā)生在S1 截面,與二維情況顯著不同.由S3截面可以看到,由于凹?xì)庵娜S性,激波第一次作用后在截面的上半部分生成了負(fù)渦量.在渦量誘導(dǎo)下,邊界面向右運(yùn)動(dòng)的速度增大,而對(duì)稱(chēng)面向右的運(yùn)動(dòng)速度則減小,并在S3 截面形成了尖釘和氣泡結(jié)構(gòu).因此可以觀察到S1 截面上處于對(duì)稱(chēng)位置處的運(yùn)動(dòng)速度要高于S2 截面上同位置的界面.而處于非對(duì)稱(chēng)x–z截面上的界面會(huì)受到兩個(gè)方向渦量(對(duì)氣柱z>0的部分來(lái)說(shuō),為向y軸正方向和z軸正方向)的疊加影響,將氣體向y=0 及z=±10 mm 方向推動(dòng),因此可以在S1 截面上觀察到橋狀結(jié)構(gòu)更短,兩個(gè)氣泡離得更近,以及600μs 后涌向邊界面的氣體在S1 截面的氣泡中形成了另一個(gè)小的氣泡,進(jìn)一步推動(dòng)S1 截面上的氣柱界面加速演化.因此,S1 截面上界面的演化速度高于S2 截面且生成了更多復(fù)雜結(jié)構(gòu).

    圖9 不同反射距離(40 mm,80 mm,120 mm)下3 個(gè)典型截面上的界面演化過(guò)程.每個(gè)圖的上半部分為渦量圖,下半部分為紋影圖,單位為μs.為便于觀察,S3 截面圖放大至S1(或S2)截面圖的1.5 倍(界面和激波位置看著略有不同)Fig.9 Evolution process of concave light gas cylinder under reshock along three typical cross section(time unit:μs).The plotting scale of the S3 cross section is 1.5 times larger than other frames in order to be observed clearly,thus the scale and the relative location of the deformed interface and shock waves shown in these frames seem different

    對(duì)于S2 截面,在演化初期凹?xì)庵缑娴难莼c二維氣柱情況基本相同,但是移動(dòng)速度受S3 截面渦量的影響減慢,形態(tài)變化更快,在600μs 時(shí)已經(jīng)形成了厚度很薄的橋狀結(jié)構(gòu),而二維工況在800 μs 后才出現(xiàn)類(lèi)似結(jié)構(gòu).更大的變化出現(xiàn)在反射激波與界面作用后.氣泡前部的渦(即圖5 中二維氣柱工況下的渦2)在演化后期強(qiáng)度減小并被氣泡所代替,通過(guò)觀察S3 截面上界面的演化過(guò)程,可以發(fā)現(xiàn)在演化后期受渦量影響,界面發(fā)展有向z=0 平面靠攏的趨勢(shì).

    反射距離對(duì)界面演化過(guò)程有重要影響.與二維情況類(lèi)似,反射激波以及反射距離的改變對(duì)界面演化的影響,主要是由于反射激波作用界面時(shí),界面和激波的形態(tài)不同所造成的,從而導(dǎo)致新生成的斜壓渦量強(qiáng)度以及分布不同.觀察S1 截面,從渦量圖中可以看到,反射激波作用后生成的負(fù)渦量在氣柱左側(cè)界面誘導(dǎo)生成了渦對(duì).由于反射距離為40 mm 時(shí)反射激波作用時(shí)間早且與固壁相互作用強(qiáng)烈,因此渦對(duì)位置以及界面整體形狀與其他工況不同.對(duì)比反射距離為80 mm 與反射距離120 mm 的工況,可以看到反射激波作用時(shí)界面的的整體形狀相似,因此反射激波作用后演化后期兩個(gè)工況界面形狀類(lèi)似,但80 mm 工況反射激波作用更早,因此環(huán)狀結(jié)構(gòu)演化更快,擴(kuò)散的更明顯.在S2 截面上同樣由于反射激波作用時(shí)界面的整體形狀相似,也滿足反射距離越短,界面擴(kuò)散越明顯的規(guī)律.

    3 定量分析

    為了更好地分析氣柱界面的演化情況,提取二維和三維凹?xì)庵淖钭髠?cè)和最右側(cè)界面位置以及高度與寬度數(shù)據(jù),分別如圖10 和圖11 所示.由于對(duì)稱(chēng)性,氣柱高度的變化可以反映遠(yuǎn)離對(duì)稱(chēng)面氣柱界面的位置變化.氣柱高度的變化趨勢(shì)總體上為:激波第一次作用后,由于激波的壓縮效應(yīng),氣柱寬度減小的同時(shí)高度快速增長(zhǎng),隨后增長(zhǎng)放緩,激波第二次作用后再次快速增長(zhǎng),隨后增長(zhǎng)趨于飽和,開(kāi)始出現(xiàn)下降趨勢(shì).演化后期高度的不斷震蕩則是由于部分氣柱界面跟隨渦對(duì)整體旋轉(zhuǎn)導(dǎo)致距離對(duì)稱(chēng)面最遠(yuǎn)處的界面不斷變化.

    圖10 二維氣柱界面位移和特征尺度隨時(shí)間的變化情況Fig.10 Variation of displacements and characteristic scales of the 2D cylindrical interface with the time

    圖11 三維凹?xì)庵煌孛嫔系奶卣鞒叨入S時(shí)間的變化情況.邊界面(a)和對(duì)稱(chēng)面(b)上氣柱高度與寬度隨時(shí)間的變化趨勢(shì)Fig.11 Variation of characteristic scales of 3D concave cylindrical interface with the time along different cross sections.The height and width of the interface along(a)the boundary cross section and(b)the symmetry cross section

    如圖10 所示,與無(wú)反射激波的工況進(jìn)行對(duì)比,反射激波作用二維氣柱界面后由于再次壓縮使氣柱的高度快速增長(zhǎng),反射激波作用變形界面后生成負(fù)渦量,氣泡1 的位置由左側(cè)逐漸順時(shí)針旋轉(zhuǎn),先從左側(cè)轉(zhuǎn)向遠(yuǎn)離對(duì)稱(chēng)軸方向引起氣柱高度的增長(zhǎng),隨后繼續(xù)旋轉(zhuǎn)并形成渦對(duì)使氣柱高度增長(zhǎng)趨于飽和.而反射距離為40 mm 的算例在1300μs 后氣柱高度突然增大,這是由于氣柱下游界面到達(dá)固壁無(wú)法繼續(xù)向前運(yùn)動(dòng),從而沿固壁向外圍擴(kuò)散.此外,在激波作用后,氣柱寬度迅速下降,隨后近似線性增長(zhǎng),激波第二次作用后寬度再次迅速下降,隨后恢復(fù)線性增長(zhǎng),最后增長(zhǎng)速度有放緩趨勢(shì).

    二維氣柱的左側(cè)界面在激波第一次作用后,位于對(duì)稱(chēng)面附近的界面在激波和渦量誘導(dǎo)的雙重作用下,幾乎以恒定速度快速向右側(cè)運(yùn)動(dòng),而由于右側(cè)界面受到激波沖擊較晚,開(kāi)始運(yùn)動(dòng)后速度又小于左側(cè)界面,此時(shí)氣柱寬度迅速減小.200μs 后離對(duì)稱(chēng)面稍遠(yuǎn)的氣泡結(jié)構(gòu)取代了對(duì)稱(chēng)面成為氣柱最左端界面,因此界面運(yùn)動(dòng)速度有了突變,運(yùn)動(dòng)速度減慢到小于右側(cè)界面運(yùn)動(dòng)速度,氣柱寬度開(kāi)始重新以線性增長(zhǎng).反射激波首先作用右側(cè)界面,雖然并沒(méi)有使右側(cè)界面的增長(zhǎng)停滯,但仍可導(dǎo)致氣柱寬度快速縮短.除了受到固壁限制的反射距離為40 mm 的工況外,其他工況右側(cè)界面均繼續(xù)以近乎勻速向右側(cè)緩慢運(yùn)動(dòng).反射激波作用左側(cè)界面后,左側(cè)界面反向向左運(yùn)動(dòng),這是由于新生成的渦對(duì)誘導(dǎo)產(chǎn)生射流的結(jié)果,氣柱寬度由此恢復(fù)增長(zhǎng),除反射距離120 mm 的算例外,其他工況在界面演化后期逐漸趨于飽和,而反射距離為120 mm時(shí)由于距離固壁較遠(yuǎn),在1600μs 時(shí)仍處于增長(zhǎng)階段.

    對(duì)于三維凹?xì)庵?邊界面上氣柱界面的運(yùn)動(dòng)速度較快,而在對(duì)稱(chēng)面上較慢.反射距離為40 mm 時(shí),由于邊界面比較貼近固壁,將很快與固壁作用,使氣柱寬度的增長(zhǎng)趨于停滯.在邊界面,由于界面運(yùn)動(dòng)速度更快,反射距離為80 mm 時(shí)也出現(xiàn)了界面與固壁作用導(dǎo)致增長(zhǎng)停滯的現(xiàn)象.而在對(duì)稱(chēng)面截面上,氣柱界面的演化更接近二維情況.

    如圖11 所示,三維凹?xì)庵缑娴奶卣鞒叽珉S時(shí)間的變化與二維情況顯著不同,主要在于三維性導(dǎo)致凹?xì)庵莼^(guò)程中存在大量的沿Z軸方向的運(yùn)動(dòng).在測(cè)量時(shí),主要考察豎直方向上兩端邊界面(最上或最下邊界面)和中心對(duì)稱(chēng)面兩個(gè)截面,且僅考慮始終在相應(yīng)截面上的界面結(jié)構(gòu).在邊界面,當(dāng)反射激波作用后,界面高度增長(zhǎng)幅度不及二維情況,并在短暫的增長(zhǎng)后迅速下降,在演化的后期才由于環(huán)狀結(jié)構(gòu)的快速擴(kuò)張以及與固壁相互作用等因素快速反彈.而在對(duì)稱(chēng)面,由于此處界面形狀演化更快而運(yùn)動(dòng)速度更慢,因此更加接近于反射距離更長(zhǎng)的二維工況.

    需要強(qiáng)調(diào)的是,由于三維凹?xì)庵谪Q直方向上出現(xiàn)了界面的負(fù)曲率效應(yīng),導(dǎo)致界面的演化與二維氣柱相比擁有更加復(fù)雜的結(jié)構(gòu),從而使對(duì)稱(chēng)面以及邊界面上的界面演化發(fā)生較大變化.這種結(jié)構(gòu)的變化與三維空間斜壓渦量的產(chǎn)生密切相關(guān),也與反射激波帶來(lái)的復(fù)雜波系結(jié)構(gòu)密不可分,不同現(xiàn)象及其演化機(jī)理仍需要進(jìn)一步的深入研究.

    4 結(jié)論

    本文通過(guò)數(shù)值模擬研究了反射激波作用下二維氣柱和三維凹?xì)庵难莼^(guò)程,考察了不同的反射距離對(duì)界面特征尺寸變化的影響規(guī)律.首先,不同反射距離下二維氣柱的演化歸因于反射激波作用于已變形氣柱時(shí)激波與氣柱的瞬時(shí)形態(tài),從而導(dǎo)致斜壓渦量的生成和分布產(chǎn)生明顯差異.通過(guò)渦量分析,發(fā)現(xiàn)反射激波作用界面時(shí),能夠在上游界面產(chǎn)生反方向的斜壓渦量,促使氣泡快速分裂成兩部分,并最終形成多個(gè)渦對(duì),以此來(lái)加速界面變形.其次,反射激波作用三維凹?xì)庵慕缑嫜莼榷S情況更加復(fù)雜,各個(gè)截面上的界面演化都存在三維性,由此產(chǎn)生的渦量分布也具有三維性.通過(guò)渦量誘導(dǎo)的界面運(yùn)動(dòng)揭示了各個(gè)截面上界面運(yùn)動(dòng)速度差異,并得到了反射激波作用后的三維環(huán)狀結(jié)構(gòu).最后,通過(guò)提取不同截面上界面幾何尺寸的演化數(shù)據(jù),對(duì)比分析了三維氣柱界面演化與二維氣柱的不同之處.本文的數(shù)值分析為進(jìn)一步開(kāi)展相關(guān)內(nèi)容的實(shí)驗(yàn)研究奠定了基礎(chǔ),也為揭示復(fù)雜界面三維性的影響機(jī)理提供了基礎(chǔ)數(shù)據(jù)和有效手段.

    猜你喜歡
    氣柱渦量不穩(wěn)定性
    預(yù)壓縮氣墊包裝系統(tǒng)靜力及動(dòng)力學(xué)特性研究
    包裝工程(2024年1期)2024-01-20 06:17:38
    激波誘導(dǎo)雙層氣柱演化的偏心效應(yīng)研究
    氣體物理(2022年2期)2022-03-31 12:49:16
    含沙空化對(duì)軸流泵內(nèi)渦量分布的影響
    磁控條件下激波沖擊三角形氣柱過(guò)程的數(shù)值研究?
    可壓縮Navier-Stokes方程平面Couette-Poiseuille流的線性不穩(wěn)定性
    自由表面渦流動(dòng)現(xiàn)象的數(shù)值模擬
    增強(qiáng)型體外反搏聯(lián)合中醫(yī)辯證治療不穩(wěn)定性心絞痛療效觀察
    航態(tài)對(duì)大型船舶甲板氣流場(chǎng)的影響
    前列地爾治療不穩(wěn)定性心絞痛療效觀察
    The application of numerical simulation of delta wing with blunt leading edge using RANS/LES hybrid method
    成人国产麻豆网| 午夜影院在线不卡| 国精品久久久久久国模美| 午夜av观看不卡| 日韩欧美一区视频在线观看| av电影中文网址| 精品午夜福利在线看| 精品亚洲成国产av| 亚洲精品日本国产第一区| 卡戴珊不雅视频在线播放| 国产精品一区二区在线观看99| 欧美在线黄色| 亚洲精品国产色婷婷电影| 免费不卡的大黄色大毛片视频在线观看| 久久精品国产亚洲av高清一级| 人妻 亚洲 视频| 一级a爱视频在线免费观看| √禁漫天堂资源中文www| 日本午夜av视频| 老司机影院毛片| 亚洲精品美女久久久久99蜜臀 | 少妇熟女欧美另类| 99九九在线精品视频| 亚洲 欧美一区二区三区| 亚洲国产欧美在线一区| 国产精品人妻久久久影院| 国产精品人妻久久久影院| 少妇人妻久久综合中文| 国产精品欧美亚洲77777| 一区二区三区四区激情视频| 99久久精品国产国产毛片| 啦啦啦啦在线视频资源| 三上悠亚av全集在线观看| 久久婷婷青草| 日韩,欧美,国产一区二区三区| 欧美激情 高清一区二区三区| 欧美精品高潮呻吟av久久| 另类亚洲欧美激情| 波野结衣二区三区在线| 精品人妻偷拍中文字幕| 最新中文字幕久久久久| 国产激情久久老熟女| 久久精品久久精品一区二区三区| 日韩免费高清中文字幕av| 亚洲av免费高清在线观看| 久久精品久久久久久久性| videos熟女内射| 99re6热这里在线精品视频| 久久免费观看电影| 中国三级夫妇交换| 午夜影院在线不卡| 999精品在线视频| 丰满乱子伦码专区| 国产成人午夜福利电影在线观看| 色哟哟·www| 黄频高清免费视频| 自线自在国产av| 波多野结衣一区麻豆| 国产精品不卡视频一区二区| 免费看av在线观看网站| 久久人人爽av亚洲精品天堂| 亚洲美女黄色视频免费看| 在线观看免费视频网站a站| 老司机影院成人| 多毛熟女@视频| 国产精品一国产av| 少妇的逼水好多| 国产精品免费视频内射| 亚洲美女视频黄频| 黄片小视频在线播放| 亚洲一区中文字幕在线| av片东京热男人的天堂| 午夜老司机福利剧场| 国产 一区精品| 精品人妻在线不人妻| 日韩欧美精品免费久久| 日本爱情动作片www.在线观看| 9色porny在线观看| 亚洲综合色网址| 亚洲激情五月婷婷啪啪| 久久久国产精品麻豆| 亚洲精品一区蜜桃| 日韩 亚洲 欧美在线| 在线观看美女被高潮喷水网站| 交换朋友夫妻互换小说| 少妇被粗大的猛进出69影院| av电影中文网址| 日韩三级伦理在线观看| 久久99一区二区三区| 亚洲视频免费观看视频| 亚洲三区欧美一区| 亚洲欧美一区二区三区久久| 视频区图区小说| 精品国产露脸久久av麻豆| 亚洲精品在线美女| 校园人妻丝袜中文字幕| 欧美日韩一级在线毛片| 制服丝袜香蕉在线| 国产一区有黄有色的免费视频| 日本av手机在线免费观看| 国产熟女午夜一区二区三区| 熟女电影av网| 国产又色又爽无遮挡免| 极品少妇高潮喷水抽搐| 777久久人妻少妇嫩草av网站| 欧美国产精品一级二级三级| 久久99精品国语久久久| 最近中文字幕2019免费版| 巨乳人妻的诱惑在线观看| 天堂俺去俺来也www色官网| 一本久久精品| 丰满少妇做爰视频| 国产亚洲欧美精品永久| 大香蕉久久成人网| 岛国毛片在线播放| 亚洲精品一区蜜桃| 日日爽夜夜爽网站| 看非洲黑人一级黄片| 尾随美女入室| 女性被躁到高潮视频| 自拍欧美九色日韩亚洲蝌蚪91| 91午夜精品亚洲一区二区三区| 国产av一区二区精品久久| 女人高潮潮喷娇喘18禁视频| 日本91视频免费播放| 午夜激情av网站| 久久午夜综合久久蜜桃| 日韩中字成人| videosex国产| 国产精品久久久av美女十八| 亚洲图色成人| 免费日韩欧美在线观看| 欧美 亚洲 国产 日韩一| 国产麻豆69| 久久久久精品人妻al黑| 超碰成人久久| 青春草亚洲视频在线观看| 成人国产麻豆网| 欧美精品一区二区大全| 欧美黄色片欧美黄色片| 国产一区有黄有色的免费视频| 亚洲欧美一区二区三区国产| www.熟女人妻精品国产| 国产精品偷伦视频观看了| 成人免费观看视频高清| tube8黄色片| 少妇 在线观看| 亚洲色图综合在线观看| 黑人猛操日本美女一级片| 精品久久久久久电影网| 超色免费av| a级毛片在线看网站| 91国产中文字幕| 日本vs欧美在线观看视频| 精品国产一区二区久久| 视频区图区小说| 精品久久久久久电影网| 日韩中文字幕视频在线看片| 大码成人一级视频| 欧美国产精品一级二级三级| 久久ye,这里只有精品| 亚洲精品av麻豆狂野| 97精品久久久久久久久久精品| 韩国av在线不卡| 欧美成人午夜精品| 久久精品久久久久久久性| 国产精品国产三级国产专区5o| 久热久热在线精品观看| 韩国av在线不卡| 丰满少妇做爰视频| 亚洲,一卡二卡三卡| 久久午夜福利片| 国产97色在线日韩免费| 亚洲熟女精品中文字幕| 亚洲一级一片aⅴ在线观看| 少妇的丰满在线观看| 中国三级夫妇交换| 久久久精品国产亚洲av高清涩受| 交换朋友夫妻互换小说| 国产高清不卡午夜福利| 国产成人精品福利久久| 国产不卡av网站在线观看| 亚洲国产精品国产精品| 少妇 在线观看| 热re99久久精品国产66热6| 一区二区三区乱码不卡18| 国产免费一区二区三区四区乱码| 中文精品一卡2卡3卡4更新| 国产一区二区三区av在线| 最近最新中文字幕大全免费视频 | www.av在线官网国产| 青春草国产在线视频| 夜夜骑夜夜射夜夜干| 久久精品久久精品一区二区三区| 中文精品一卡2卡3卡4更新| 免费播放大片免费观看视频在线观看| 黄片无遮挡物在线观看| 18在线观看网站| 婷婷色综合大香蕉| 久久久国产欧美日韩av| 亚洲欧美日韩另类电影网站| 一边亲一边摸免费视频| 搡女人真爽免费视频火全软件| 久久99蜜桃精品久久| 你懂的网址亚洲精品在线观看| 黑丝袜美女国产一区| 在线精品无人区一区二区三| 美女中出高潮动态图| 丝瓜视频免费看黄片| 一区二区三区激情视频| 好男人视频免费观看在线| 精品一区在线观看国产| 免费不卡的大黄色大毛片视频在线观看| 少妇人妻久久综合中文| 男女下面插进去视频免费观看| 国产成人精品福利久久| 老汉色av国产亚洲站长工具| 热99久久久久精品小说推荐| 人人妻人人爽人人添夜夜欢视频| 深夜精品福利| 国产乱来视频区| 最近的中文字幕免费完整| 成年人午夜在线观看视频| 久久ye,这里只有精品| 国产人伦9x9x在线观看 | 一区二区三区四区激情视频| 毛片一级片免费看久久久久| 久久国产亚洲av麻豆专区| 一级毛片 在线播放| 中文欧美无线码| 午夜福利在线观看免费完整高清在| 另类亚洲欧美激情| 午夜福利乱码中文字幕| 免费少妇av软件| 一级爰片在线观看| 波多野结衣av一区二区av| 欧美成人精品欧美一级黄| 国产精品久久久久成人av| 高清不卡的av网站| 国产精品女同一区二区软件| 国产激情久久老熟女| 在线看a的网站| 波野结衣二区三区在线| 成年人免费黄色播放视频| 精品一品国产午夜福利视频| 在线观看www视频免费| 久久久久网色| 午夜福利影视在线免费观看| 国产一区二区 视频在线| 国产精品免费大片| 少妇人妻 视频| 国产麻豆69| 国产一区亚洲一区在线观看| 精品少妇久久久久久888优播| 亚洲欧洲日产国产| 精品一区二区三区四区五区乱码 | 久久久久精品人妻al黑| 亚洲图色成人| 成人二区视频| 久久久精品区二区三区| 99热全是精品| 老司机亚洲免费影院| 亚洲综合精品二区| 色播在线永久视频| 在线观看国产h片| 国产精品一区二区在线不卡| 国产亚洲最大av| 亚洲精品国产色婷婷电影| 夜夜骑夜夜射夜夜干| √禁漫天堂资源中文www| 亚洲国产精品一区三区| 中文字幕精品免费在线观看视频| 人妻系列 视频| 久久精品国产鲁丝片午夜精品| 午夜免费观看性视频| 亚洲精华国产精华液的使用体验| 日本黄色日本黄色录像| 香蕉丝袜av| 亚洲美女视频黄频| 国产av国产精品国产| 国产精品熟女久久久久浪| 伊人久久大香线蕉亚洲五| 免费在线观看完整版高清| freevideosex欧美| 国产一区亚洲一区在线观看| 曰老女人黄片| 免费黄频网站在线观看国产| 免费大片黄手机在线观看| 欧美亚洲 丝袜 人妻 在线| 亚洲av福利一区| 亚洲精品av麻豆狂野| 天堂俺去俺来也www色官网| 国产一区有黄有色的免费视频| 国产在线免费精品| 欧美精品一区二区免费开放| 在线观看国产h片| a级毛片黄视频| 波多野结衣av一区二区av| 熟女电影av网| 亚洲内射少妇av| 春色校园在线视频观看| 日日啪夜夜爽| 老司机影院毛片| 国产有黄有色有爽视频| 制服诱惑二区| 国产精品女同一区二区软件| 久热久热在线精品观看| 亚洲视频免费观看视频| 天天操日日干夜夜撸| 美女国产视频在线观看| 中文字幕av电影在线播放| 少妇熟女欧美另类| 丝袜喷水一区| 国产成人一区二区在线| 18在线观看网站| 多毛熟女@视频| 美国免费a级毛片| 丝袜在线中文字幕| 男人舔女人的私密视频| 久久婷婷青草| 免费大片黄手机在线观看| 成人亚洲欧美一区二区av| 免费观看a级毛片全部| 成人影院久久| 国产片特级美女逼逼视频| 老熟女久久久| 岛国毛片在线播放| 如日韩欧美国产精品一区二区三区| 国产精品三级大全| 亚洲欧美中文字幕日韩二区| 国产精品麻豆人妻色哟哟久久| 天堂中文最新版在线下载| 一本—道久久a久久精品蜜桃钙片| 一区二区三区乱码不卡18| tube8黄色片| 国产成人精品久久二区二区91 | 精品国产一区二区三区四区第35| 成人亚洲欧美一区二区av| 日日撸夜夜添| 精品久久久精品久久久| 三上悠亚av全集在线观看| 精品少妇黑人巨大在线播放| 亚洲精品,欧美精品| av国产精品久久久久影院| 妹子高潮喷水视频| 亚洲在久久综合| 日韩精品有码人妻一区| 大片免费播放器 马上看| 乱人伦中国视频| 天堂中文最新版在线下载| 国产欧美日韩综合在线一区二区| 日本av手机在线免费观看| www.自偷自拍.com| 叶爱在线成人免费视频播放| 成人国产av品久久久| 亚洲国产成人一精品久久久| 亚洲中文av在线| 最新中文字幕久久久久| 咕卡用的链子| 欧美老熟妇乱子伦牲交| 99热国产这里只有精品6| 高清欧美精品videossex| 欧美精品高潮呻吟av久久| 日本免费在线观看一区| 日本wwww免费看| 久久97久久精品| 国产精品蜜桃在线观看| 久久鲁丝午夜福利片| 91国产中文字幕| 国产精品一区二区在线不卡| 一区二区三区激情视频| 欧美激情 高清一区二区三区| 日本午夜av视频| av卡一久久| 精品亚洲成a人片在线观看| 亚洲欧美一区二区三区黑人 | 黄色 视频免费看| 精品少妇黑人巨大在线播放| 天天操日日干夜夜撸| 黑人猛操日本美女一级片| 国产精品一二三区在线看| 欧美+日韩+精品| 五月伊人婷婷丁香| 午夜福利视频精品| 女人被躁到高潮嗷嗷叫费观| 岛国毛片在线播放| 欧美日韩精品成人综合77777| 国精品久久久久久国模美| 亚洲av在线观看美女高潮| 免费女性裸体啪啪无遮挡网站| 熟女少妇亚洲综合色aaa.| 青春草国产在线视频| 国产探花极品一区二区| 色婷婷av一区二区三区视频| 下体分泌物呈黄色| 91精品三级在线观看| 少妇的丰满在线观看| 国产成人aa在线观看| 国产精品.久久久| 久久久久国产一级毛片高清牌| 777久久人妻少妇嫩草av网站| 伦精品一区二区三区| 亚洲欧美成人综合另类久久久| 一本大道久久a久久精品| 日韩一区二区视频免费看| 国产精品无大码| 一本大道久久a久久精品| 在线精品无人区一区二区三| 精品久久久久久电影网| 久久精品夜色国产| 欧美最新免费一区二区三区| 色哟哟·www| 三上悠亚av全集在线观看| 满18在线观看网站| 国产成人av激情在线播放| 久久久久久久亚洲中文字幕| 视频在线观看一区二区三区| 午夜日韩欧美国产| 国产精品一区二区在线不卡| 999精品在线视频| 肉色欧美久久久久久久蜜桃| 国产精品一区二区在线不卡| 99久久人妻综合| 麻豆av在线久日| 国产成人a∨麻豆精品| 免费少妇av软件| 这个男人来自地球电影免费观看 | 国产精品三级大全| 国产无遮挡羞羞视频在线观看| 美女大奶头黄色视频| 日日爽夜夜爽网站| 午夜福利,免费看| 日韩制服骚丝袜av| a 毛片基地| 亚洲精品美女久久久久99蜜臀 | 伊人久久大香线蕉亚洲五| 亚洲国产精品一区二区三区在线| 天堂8中文在线网| 丁香六月天网| 最新的欧美精品一区二区| 亚洲成人手机| 蜜桃在线观看..| 亚洲美女视频黄频| 人妻少妇偷人精品九色| 日本黄色日本黄色录像| 免费不卡的大黄色大毛片视频在线观看| 国产 精品1| 久久久久精品人妻al黑| 亚洲精品乱久久久久久| 久久久精品区二区三区| 欧美日韩综合久久久久久| 久久久精品94久久精品| www.自偷自拍.com| 母亲3免费完整高清在线观看 | 狂野欧美激情性bbbbbb| 女性生殖器流出的白浆| 两个人看的免费小视频| 欧美av亚洲av综合av国产av | 亚洲精品日韩在线中文字幕| 亚洲婷婷狠狠爱综合网| 秋霞伦理黄片| 少妇人妻久久综合中文| 亚洲精品成人av观看孕妇| 99精国产麻豆久久婷婷| 久久午夜综合久久蜜桃| 波多野结衣av一区二区av| 欧美人与性动交α欧美精品济南到 | 最近最新中文字幕大全免费视频 | av在线app专区| 精品人妻一区二区三区麻豆| 亚洲精品久久久久久婷婷小说| av电影中文网址| 国产在线视频一区二区| 国产一区亚洲一区在线观看| 乱人伦中国视频| 日韩伦理黄色片| 侵犯人妻中文字幕一二三四区| videosex国产| 国产一区二区 视频在线| 精品亚洲成a人片在线观看| 街头女战士在线观看网站| 色婷婷久久久亚洲欧美| 女性生殖器流出的白浆| 成人影院久久| av在线老鸭窝| 欧美日韩一级在线毛片| 秋霞在线观看毛片| 国产成人精品久久久久久| 久久青草综合色| av视频免费观看在线观看| 精品一区二区三区四区五区乱码 | 欧美变态另类bdsm刘玥| 国产精品亚洲av一区麻豆 | 99热全是精品| 精品国产一区二区三区久久久樱花| 两个人看的免费小视频| 巨乳人妻的诱惑在线观看| 免费在线观看黄色视频的| 在线免费观看不下载黄p国产| 人人妻人人澡人人爽人人夜夜| 国产亚洲最大av| 91成人精品电影| 欧美人与性动交α欧美精品济南到 | 日日摸夜夜添夜夜爱| 精品国产一区二区三区四区第35| 精品国产一区二区久久| 国产日韩欧美视频二区| 成人毛片a级毛片在线播放| 婷婷色综合www| 国产精品一区二区在线观看99| 亚洲色图综合在线观看| 宅男免费午夜| 国产亚洲最大av| 国产av一区二区精品久久| 可以免费在线观看a视频的电影网站 | 日本猛色少妇xxxxx猛交久久| 久久久国产欧美日韩av| 精品国产乱码久久久久久小说| 人妻少妇偷人精品九色| 国产高清不卡午夜福利| 久久99热这里只频精品6学生| 免费看不卡的av| 在线观看三级黄色| 国产乱人偷精品视频| 婷婷色综合大香蕉| 亚洲三区欧美一区| 1024视频免费在线观看| 亚洲欧美中文字幕日韩二区| 久久精品国产亚洲av高清一级| 国产精品免费大片| 天天躁狠狠躁夜夜躁狠狠躁| 中文字幕人妻熟女乱码| 建设人人有责人人尽责人人享有的| 寂寞人妻少妇视频99o| 精品久久久精品久久久| 成年人免费黄色播放视频| 亚洲av电影在线观看一区二区三区| 秋霞伦理黄片| 又粗又硬又长又爽又黄的视频| 自线自在国产av| 日韩av在线免费看完整版不卡| 精品久久久精品久久久| 精品少妇黑人巨大在线播放| 日韩中文字幕视频在线看片| 伊人久久国产一区二区| 搡女人真爽免费视频火全软件| 亚洲精品日本国产第一区| 国产欧美亚洲国产| 久久综合国产亚洲精品| 秋霞伦理黄片| 亚洲成国产人片在线观看| 伊人久久大香线蕉亚洲五| 两性夫妻黄色片| 精品国产超薄肉色丝袜足j| 麻豆av在线久日| 久久久久精品性色| 成年美女黄网站色视频大全免费| av电影中文网址| 欧美激情高清一区二区三区 | 最近最新中文字幕免费大全7| videos熟女内射| 成人国产av品久久久| 免费高清在线观看视频在线观看| 国产老妇伦熟女老妇高清| 波野结衣二区三区在线| 99久久综合免费| videosex国产| 大片免费播放器 马上看| 国产高清不卡午夜福利| 人妻人人澡人人爽人人| 看十八女毛片水多多多| 精品少妇内射三级| 欧美日韩综合久久久久久| 天天躁日日躁夜夜躁夜夜| 国产精品av久久久久免费| 美国免费a级毛片| 国产熟女欧美一区二区| 一本大道久久a久久精品| 男女免费视频国产| 老司机亚洲免费影院| 美女大奶头黄色视频| 亚洲欧美一区二区三区国产| 亚洲精品日韩在线中文字幕| 国产精品国产三级国产专区5o| 在线观看免费日韩欧美大片| 蜜桃国产av成人99| 69精品国产乱码久久久| 女人被躁到高潮嗷嗷叫费观| 日韩精品有码人妻一区| 亚洲精品一区蜜桃| 色吧在线观看| 免费少妇av软件| 免费久久久久久久精品成人欧美视频| 久久狼人影院| 97在线视频观看| a级片在线免费高清观看视频| 久久午夜综合久久蜜桃| 欧美人与善性xxx| 日本av手机在线免费观看| 看免费成人av毛片| 午夜日韩欧美国产| 黄频高清免费视频| a 毛片基地| 免费人妻精品一区二区三区视频| 久久久久久免费高清国产稀缺| 如日韩欧美国产精品一区二区三区| 国产成人欧美| 视频区图区小说| 精品亚洲成国产av| 欧美精品国产亚洲| 精品人妻在线不人妻| 青春草亚洲视频在线观看| 中文字幕人妻丝袜一区二区 | 黑人欧美特级aaaaaa片| 99re6热这里在线精品视频|