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

    振動激發(fā)對激波反射的影響

    2017-11-20 03:14:49彭俊張子健周凱胡宗民姜宗林
    航空學(xué)報(bào) 2017年8期
    關(guān)鍵詞:極線馬赫馬赫數(shù)

    彭俊, 張子健, 周凱, 胡宗民,*, 姜宗林

    1.中國科學(xué)院力學(xué)研究所 高溫氣體動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100190 2.中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院, 北京 100049

    振動激發(fā)對激波反射的影響

    彭俊1,2, 張子健1,2, 周凱1,2, 胡宗民1,2,*, 姜宗林1,2

    1.中國科學(xué)院力學(xué)研究所 高溫氣體動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100190 2.中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院, 北京 100049

    定常激波反射分為規(guī)則反射和馬赫反射,在不同條件下2種反射結(jié)構(gòu)之間會相互轉(zhuǎn)變。高超聲速流動中的激波反射問題常面臨高溫氣體效應(yīng),隨著溫度逐漸升高,最先出現(xiàn)的是空氣分子振動激發(fā)。通過理論分析和定量計(jì)算,研究了振動激發(fā)對激波反射及其轉(zhuǎn)變規(guī)律的影響。首先給出考慮振動激發(fā)的空氣熱力學(xué)模型,并分析其與量熱完全氣體的差異以及對激波關(guān)系的影響;接著分析在規(guī)則反射和馬赫反射中,振動激發(fā)對激波反射流場的影響規(guī)律;最后討論振動激發(fā)對激波反射轉(zhuǎn)變2個準(zhǔn)則點(diǎn)的影響。研究結(jié)果表明,振動激發(fā)使激波極線的整體輪廓變大,且這種差異在經(jīng)過一次激波反射之后被明顯放大,會對激波反射的流場產(chǎn)生重大影響;對于激波反射的轉(zhuǎn)變,振動激發(fā)使轉(zhuǎn)變的2個準(zhǔn)則點(diǎn)都變大,且對規(guī)則反射向馬赫反射轉(zhuǎn)變的脫體準(zhǔn)則影響更大。

    激波反射; 振動激發(fā); 馬赫反射; 轉(zhuǎn)變準(zhǔn)則; 高超聲速流動

    激波廣泛存在于自然界和人類科技活動中,是氣體動力學(xué)最具特色的基本物理現(xiàn)象之一,是能夠在流動內(nèi)部誘導(dǎo)旋渦的重要物理過程,表現(xiàn)出強(qiáng)間斷與非線性的氣動物理特征。在激波傳播過程中,當(dāng)傳播條件發(fā)生改變時,包括壁面出現(xiàn)、流動介質(zhì)屬性改變、其他流場結(jié)構(gòu)存在等,激波會與其產(chǎn)生相互作用而發(fā)生反射,形成激波反射流場結(jié)構(gòu)。在不同條件下,會形成各種具有不同特征的激波反射結(jié)構(gòu),包括規(guī)則反射、馬赫反射等,出現(xiàn)反射激波、馬赫桿、滑移線、三波點(diǎn)等典型流場結(jié)構(gòu)[1]。在某些條件下,甚至?xí)纬呻p馬赫反射、弱激波反射以及一些過渡或反轉(zhuǎn)激波反射結(jié)構(gòu),出現(xiàn)雙馬赫桿、前向射流[2-3]、多超聲速區(qū)[4-5]、拐點(diǎn)[1,6]等特殊流場結(jié)構(gòu)。當(dāng)產(chǎn)生激波反射的條件發(fā)生變化時,各種不同反射結(jié)構(gòu)之間會發(fā)生相互轉(zhuǎn)變,存在各種轉(zhuǎn)變準(zhǔn)則[1,7-8]。自von Neumann[9-10]于20世紀(jì)40年代對激波反射現(xiàn)象進(jìn)行深入研究開始,經(jīng)過70多年的努力,人們提出了各種分析理論、試驗(yàn)技術(shù)和計(jì)算方法,發(fā)現(xiàn)了各種激波反射結(jié)構(gòu),得到了各種波系結(jié)構(gòu)產(chǎn)生的不同物理機(jī)制,對激波反射現(xiàn)象的認(rèn)識越來越深入并日臻完善。

    激波反射問題得到廣泛深入的研究,除了由于其本身具有的重要流動物理意義之外,更重要的是由于其具有很強(qiáng)的工程應(yīng)用背景,且在不同應(yīng)用背景之下,激波反射現(xiàn)象具有不同的利弊。在超/高超聲速飛行器機(jī)身外表面、翼面以及進(jìn)氣道內(nèi),激波在壁面上的反射,會與壁面附近的邊界層發(fā)生相互作用,在激波反射點(diǎn)附近區(qū)域,流動現(xiàn)象很復(fù)雜,可能會發(fā)生流動分離、激波非定常運(yùn)動、旋渦運(yùn)動、壓力波動、混合和湍流等問題[11-12]。一旦出現(xiàn)流動分離,產(chǎn)生了分離泡,激波波系會發(fā)生變化,進(jìn)而引起壓力載荷的變化,在流動再附點(diǎn)還會引起熱流的增加。在超燃沖壓發(fā)動機(jī)進(jìn)氣道和隔離段中,超聲速氣流進(jìn)入進(jìn)氣道后形成的定常激波在上下壁面間不斷反射,形成一系列激波串,完成來流壓縮過程[13-14]。在發(fā)動機(jī)尾噴管的超聲速射流中,斜激波在對稱軸反射形成馬赫盤以及一系列激波/膨脹波相關(guān)的多級反射的菱形結(jié)構(gòu),對發(fā)動機(jī)性能產(chǎn)生重要影響[13-15]。在預(yù)混可燃?xì)怏w中的激波反射能形成局部高溫高壓區(qū),可實(shí)現(xiàn)爆轟直接點(diǎn)火或爆燃向爆轟轉(zhuǎn)變,是爆轟現(xiàn)象研究和爆轟發(fā)動機(jī)研制的關(guān)鍵技術(shù)[16-19]。工程上利用與控制激波反射,非常依賴于能否清晰地認(rèn)識各種反射結(jié)構(gòu)的特征及其產(chǎn)生機(jī)理,以及在產(chǎn)生條件發(fā)生改變時,各種激波反射結(jié)構(gòu)之間的轉(zhuǎn)變準(zhǔn)則。

    自21世紀(jì)以來,隨著高超聲速飛行器的深入研究,激波反射的研究范圍逐漸拓寬到高超聲速領(lǐng)域。高超聲速氣流通過激波甚至多道激波壓縮,空氣巨大的動能轉(zhuǎn)變?yōu)閮?nèi)能,激波波后溫度急劇升高。空氣在如此高的溫度下將會發(fā)生分子振動激發(fā)、解離和電離等熱化學(xué)反應(yīng),表現(xiàn)出不同于量熱完全氣體的特性,使流動規(guī)律發(fā)生變化,這就是高溫氣體效應(yīng)[20-21]。在高溫氣體效應(yīng)影響下,激波反射會具有與量熱完全氣體流動不同的現(xiàn)象和特性,出現(xiàn)與量熱完全氣體理論預(yù)測不一樣的激波反射結(jié)構(gòu)和特征,同時各種反射之間的轉(zhuǎn)變準(zhǔn)則也會發(fā)生變化[22]。Tarnavskii[23]采用等效比熱比方法研究了考慮大氣真實(shí)物理特性時在高超聲速飛行器發(fā)動機(jī)進(jìn)氣道內(nèi)氣流偏轉(zhuǎn)角對不同激波反射結(jié)構(gòu)形成的影響,以及計(jì)算了在反射結(jié)構(gòu)中不同區(qū)域內(nèi)的氣流參數(shù)。高云亮等[24]通過化學(xué)反應(yīng)流的數(shù)值模擬研究發(fā)現(xiàn),強(qiáng)激波誘導(dǎo)的高溫氣體效應(yīng)使馬赫反射的三波點(diǎn)軌跡角變小、馬赫桿高度變低、楔頂附體激波傾角變小,并進(jìn)一步加劇了馬赫桿的向前突出。李季[25]建立了化學(xué)非平衡流激波關(guān)系的理論計(jì)算方法,得到了考慮化學(xué)非平衡效應(yīng)的激波極線圖,并結(jié)合數(shù)值模擬研究了化學(xué)非平衡對定常激波規(guī)則反射和馬赫反射遲滯效應(yīng)的影響。因此,研究考慮高溫氣體效應(yīng)下的激波反射問題,對進(jìn)一步發(fā)展激波反射理論、認(rèn)識新的激波反射物理規(guī)律和解決更為復(fù)雜的工程實(shí)際問題具有重要意義。

    空氣所處的溫度區(qū)間不同,高溫氣體效應(yīng)中占主導(dǎo)的熱化學(xué)反應(yīng)過程也不同,如表1(表中p為壓力,且1 atm=1.01×105Pa)所示[11],隨著溫度逐漸升高,最先出現(xiàn)的是空氣分子振動激發(fā);且在較大的溫度范圍內(nèi)(800~2 500 K),主要的熱化學(xué)反應(yīng)過程是分子振動激發(fā),而分子解離、電離等化學(xué)反應(yīng)過程基本沒發(fā)生,此時可以只考慮振動激發(fā)對流動的影響。因此,研究振動激發(fā)對激波反射的影響是研究高溫氣體效應(yīng)對激波反射影響的第1步。

    在各種復(fù)雜激波反射結(jié)構(gòu)中,定常激波反射是比較簡單和典型的類型。如圖1所示,定常激波反射分為規(guī)則反射和馬赫反射,圖中θ為入射氣流偏轉(zhuǎn)角,β為激波角,Ma1、Ma2和Ma3分別為來流、入射波后和反射波后的馬赫數(shù),Ma4為馬赫桿后的馬赫數(shù)。定常激波反射最早由Mach[26]通過試驗(yàn)研究首次發(fā)現(xiàn),并由von Neumann[9-10]進(jìn)行了深入研究。其中規(guī)則反射結(jié)構(gòu)由2道激波組成:入射激波和反射激波。而馬赫反射結(jié)構(gòu)由3道激波(入射激波、反射激波以及馬赫桿)和1條滑移線組成,這4個間斷面相交于1個三波點(diǎn)。因此,以定常激波反射的2種反射結(jié)構(gòu)為研究對象,對初步認(rèn)識激波反射現(xiàn)象中的高溫氣體效應(yīng)具有重要作用。

    本文主要以定常激波反射的2種反射結(jié)構(gòu),即規(guī)則反射和馬赫反射為研究對象,研究振動激發(fā)對激波反射流場特性和激波反射轉(zhuǎn)變準(zhǔn)則的影響。首先給出振動激發(fā)下空氣的熱力學(xué)模型,然后利用激波極線討論振動激發(fā)對激波關(guān)系和激波反射流場特性的影響,最后討論振動激發(fā)對激波反射轉(zhuǎn)變準(zhǔn)則的影響。

    表1不同溫度區(qū)間空氣所發(fā)生的熱化學(xué)反應(yīng)(p=1atm)[11]

    Table1Temperaturerangesofthermochemicalreactionprocessforair(p=1atm)[11]

    Temperaturerange/KThermochemicalreactionprocess<800Caloricallyperfectgas800-2500Vibrationexcitation2500-4000O2beginstodissociate4000-9000N2beginstodissociate;O2almostcompletelydissociated>9000N2almostcompletelydissociated;ionizationbegins

    圖1 定常激波反射Fig.1 Steady shock reflections

    1 振動激發(fā)下的空氣熱力學(xué)模型

    假設(shè)空氣是由79%的N2和21%的O2組成的,即

    (1)

    各組分的相對分子質(zhì)量為

    MN2=28,MO2=32

    (2)

    則各組分氣體的質(zhì)量分?jǐn)?shù)分別為

    (3)

    (4)

    各組分氣體的氣體常數(shù)分別為

    (5)

    (6)

    式中:R0=8.314 J/mol為普適氣體常數(shù)??諝獾臍怏w常數(shù)為

    R=mN2RN2+mO2RO2=288.3 J/(kg·K)

    (7)

    常溫下,空氣一般按量熱完全氣體來處理,其比焓h、定壓比熱cp和比熱比γ等熱力學(xué)參數(shù)滿足的關(guān)系為

    (8)

    (9)

    (10)

    式中:T為溫度。

    考慮振動激發(fā)后,空氣已不再是量熱完全氣體,比焓h中將會多出振動能項(xiàng)ev,即

    (11)

    (12)

    式中:Tve,N2和Tve,O2分別為N2和O2的振動特征溫度,取值為Tve,N2=3 340 K,Tve,O2=2 230 K[27]。此時,定壓比熱和比熱比則分別為

    (13)

    (14)

    根據(jù)量熱完全氣體的熱力學(xué)公式(式(8)~式(10))以及考慮振動激發(fā)的熱力學(xué)公式(式(11)~式(14)),2種空氣模型的差異隨溫度的變化如圖2所示。可見,由于分子振動能激發(fā),2種模型的差異主要體現(xiàn)在能量上,以及與能量相關(guān)的氣體屬性,如比熱、比熱比等,進(jìn)而產(chǎn)生流動上的差異。當(dāng)溫度在300 K以下時,2種模型幾乎一模一樣,并沒有表現(xiàn)出任何可觀的差異。隨著溫度的上升,空氣分子的振動自由度開始激發(fā),2種模型的差異逐漸變得明顯。當(dāng)溫度上升到800 K時,2種空氣模型比焓的差異δh/h為2.50%,定壓比熱的差異δcp/cp為8.67%,考慮振動激發(fā)的比熱比從量熱完全氣體的1.4降低到1.353,此時振動激發(fā)對空氣熱力學(xué)性質(zhì)的影響已不可忽略[22]。當(dāng)溫度達(dá)到2 500 K時,2種空氣模型比焓的差異為12.64%,定壓比熱的差異為25.11%,考慮振動激發(fā)的比熱比從量熱完全氣體的1.4降低到1.296,此時空氣分子振動自由度幾乎完全激發(fā),將極大地影響空氣的熱力學(xué)性質(zhì)及以空氣為介質(zhì)的流動。

    對于在溫度范圍為800~2 500 K內(nèi)只考慮振動激發(fā)而不考慮化學(xué)反應(yīng)所帶來的誤差,可以通過比較正激波波后參數(shù)在量熱完全氣體、振動激發(fā)和熱化學(xué)反應(yīng)(振動激發(fā)+化學(xué)反應(yīng))這3種空氣模型下的差異得到。設(shè)定來流溫度T1=226 K、壓力p1=0.1 atm、馬赫數(shù)Ma1=7,此時量熱完全氣體正激波波后溫度T2=2 366 K,接近所給溫度范圍的上限。振動激發(fā)的計(jì)算方法將在2.1節(jié)中介紹,熱化學(xué)反應(yīng)的計(jì)算則采用Gaseq小程序,3種空氣模型的計(jì)算結(jié)果如表2所示,表中p2與Ma2分別為正激波后的壓力與馬赫數(shù),括號內(nèi)的百分比代表振動激發(fā)或熱化學(xué)反應(yīng)的結(jié)果與量熱完全氣體相應(yīng)結(jié)果的差異。從表中數(shù)據(jù)可見,即使溫度已經(jīng)比較接近所給溫度范圍的上限,振動激發(fā)的影響仍起主導(dǎo)作用,化學(xué)反應(yīng)的影響較小。這說明,在2 500 K以下,只考慮振動激發(fā)對流動的影響是合理的。

    圖2 2種空氣模型的對比Fig.2 Comparison of two gas models

    表2 量熱完全氣體、振動激發(fā)和化學(xué)反應(yīng)3種空氣模型下正激波波后參數(shù)對比(Ma1=7)

    Table2Comparisonofparametersbehindnormalshockforthreegasmodelsofcaloricallyperfectgas,vibrationexcitationandthermochemicalreaction(Ma1=7)

    GasmodelT2/Kp2/atmMa2Caloricallyperfectgas23665.700.397Vibrationexcitation2114(-10.7%)5.89(3.3%)0.379(-4.5%)Thermochemicalreaction2090(-11.7%)5.91(3.7%)0.375(-5.5%)

    2 振動激發(fā)下的激波關(guān)系

    在激波反射問題中,無論是入射激波、反射激波或者是馬赫桿都要用到激波關(guān)系。因此,在研究振動激發(fā)對激波反射的影響之前,先來研究振動激發(fā)對激波關(guān)系的影響。在von Neumann提出的分析定常激波反射問題的雙激波理論和三激波理論中[10],限制激波反射求解的邊界條件是氣流偏轉(zhuǎn)角和靜壓,因此Kawamura和Saito[28]提出用氣流偏轉(zhuǎn)角和靜壓組成的關(guān)系曲線——激波極線來表示和分析激波反射過程。同時,激波極線也是全面表示激波關(guān)系的一個有力手段。

    2.1 理論計(jì)算方法與驗(yàn)證

    激波極線可以通過斜激波關(guān)系得到。即已知一個氣流偏轉(zhuǎn)角θ,通過斜激波關(guān)系得到激波極線上的一個點(diǎn)。但注意到,由于最大氣流偏轉(zhuǎn)角θm未知,θ的取值范圍需要進(jìn)一步確定,給該方法帶來一定困難。然而,激波角β的取值范圍一直都是0°~90°。因此,另一個可行的途徑是已知一個激波角β,通過正激波關(guān)系得到斜激波解和氣流偏轉(zhuǎn)角θ,從而對應(yīng)激波極線上的一個點(diǎn)。量熱完全氣體的激波關(guān)系存在解析解,因此很容易得到量熱完全氣體的激波極線圖。然而,考慮振動激發(fā)時,空氣的比焓是溫度的非線性函數(shù),如式(11)~式(12),正激波關(guān)系沒有解析表達(dá)式。此時,可從基本方程出發(fā)推導(dǎo)得到適用于計(jì)算機(jī)求解的牛頓迭代公式,其表達(dá)式為

    (15)

    (16)

    (17)

    (18)

    為對該理論計(jì)算方法進(jìn)行驗(yàn)證,在考慮振動激發(fā)的激波極線上取氣流偏轉(zhuǎn)角θ分別為20°、30° 和40°,將所得理論解與數(shù)值模擬的結(jié)果進(jìn)行對比,如表3所示。來流溫度T1=226 K,馬赫數(shù)Ma1=10。數(shù)值模擬求解的是考慮振動激發(fā)的二維Euler方程,空間上采用AUSMPW+格式[29-30]離散,時間上采用LUSGS隱式格式[31]推進(jìn),并采用MUSCL格式[32]對原始變量進(jìn)行重構(gòu)。

    從表3中各個氣流偏轉(zhuǎn)角θ下跨斜激波溫度比T21、壓力比p21和波后馬赫數(shù)Ma2等參數(shù)的理論解和數(shù)值模擬結(jié)果的對比發(fā)現(xiàn),兩者是非常一致的,說明以上提出的激波極線的理論計(jì)算方法是正確的。圖3為考慮振動激發(fā)時30° 斜劈的壓力云圖。從圖中清晰可見一道斜激波,通過與圖上根據(jù)激波角的理論值37.36° 畫出的一條直線對比發(fā)現(xiàn),兩者是平行的,即理論激波角與數(shù)值模擬得到的結(jié)果也是一致的。

    表3考慮振動激發(fā)時斜激波關(guān)系的理論解和CFD結(jié)果對比(T1=226K,Ma1=10)

    Table3ComparisonoftheoreticalandCFDsolutionsforobliqueshockrelationwithvibrationexcitation(T1=226K,Ma1=10)

    MethodT21p21Ma2θ=20°θ=30°θ=40°θ=20°θ=30°θ=40°θ=20°θ=30°θ=40°Theory4.367.4111.0721.6943.8271.804.443.052.02CFD4.367.4111.0721.6943.8271.804.443.052.02

    圖3 T1=226 K, Ma1=10時30° 斜劈的壓力云圖(CFD)Fig.3 Pressure contours of a 30° wedge at T1=226 K, Ma1=10 (CFD)

    2.2 對比分析

    圖4為量熱完全氣體和考慮振動激發(fā)的激波極線,激波前溫度T1=226 K,馬赫數(shù)Ma1=10。從圖4可見,激波波前狀態(tài)相同時,同一個氣流偏轉(zhuǎn)角對應(yīng)2個靜壓比和激波角,即出現(xiàn)2個激波解,其中靜圧比和激波角大的為強(qiáng)激波解,另一個為弱激波解。對比量熱完全氣體的結(jié)果,考慮振動激發(fā)后激波極線的整體輪廓變大了。對于強(qiáng)激波解,考慮振動激發(fā)的波后壓力明顯變大,激波角也變大。而對于弱激波解,考慮振動激發(fā)的波后壓力變小,激波角也變小。在激波極線中,最明顯的特征就是最大氣流偏轉(zhuǎn)角,它是斜激波與脫體激波的臨界角,與規(guī)則反射向馬赫反射的轉(zhuǎn)變密切相關(guān)。圖4的結(jié)果表明,考慮振動激發(fā)后最大氣流偏轉(zhuǎn)角與量熱完全氣體的相比明顯增大。后面的討論指出,這將對規(guī)則反射和馬赫反射之間的轉(zhuǎn)變產(chǎn)生重要影響。

    圖4 激波極線對比Fig.4 Comparison of shock polars

    考慮振動激發(fā)后,強(qiáng)激波解與弱激波解的波后壓力變化方向不同。而強(qiáng)激波解所對應(yīng)的波后流動一般為亞聲速,弱激波解的波后一般為超聲速流動。這說明,振動激發(fā)對強(qiáng)激波和弱激波后的高超聲速飛行器的氣動力將會產(chǎn)生不同方向的影響,尤其對周圍存在跨聲速流動的飛行器如航天飛機(jī)等的俯仰力矩影響更大。如圖5所示,當(dāng)鈍頭體以大攻角飛行時,背風(fēng)面壓力受振動激發(fā)的影響較迎風(fēng)面的小得多。迎風(fēng)面上,在頭部產(chǎn)生弓形激波,波后存在一個聲速面。聲速面前是亞聲速流動,對應(yīng)強(qiáng)激波解,振動激發(fā)使波后壓力變大;而聲速面后是超聲速流動,對應(yīng)弱激波解,振動激發(fā)使波后壓力變小,這就產(chǎn)生一個由振動激發(fā)導(dǎo)致的附加俯仰力矩。

    取氣流偏轉(zhuǎn)角θ=30°,求出強(qiáng)激波解和弱激波解所對應(yīng)的跨激波壓力比、溫度比、密度比(ρ21)、速度比(u21)等參數(shù),如表4所示。可見,振動激發(fā)對激波后溫度、密度和馬赫數(shù)的影響較大,超過10%。其中,振動激發(fā)對強(qiáng)激波波后密度的影響超過20%,由于激波前后密度比與激波脫體距離直接相關(guān),這說明振動激發(fā)對強(qiáng)激波的激波脫體距離也會產(chǎn)生重大影響。而振動激發(fā)對波后壓力和激波角的影響較小,在5%以內(nèi)。對于激波后的速度,振動激發(fā)對弱激波的影響較小,而對強(qiáng)激波的影響很大,達(dá)19.63%。另外,除激波后的溫度和密度之外,振動激發(fā)對激波后的壓力、速度、馬赫數(shù)和激波角的影響,在強(qiáng)激波和弱激波中影響方向是相反的,具有不同的變化趨勢。

    圖5 大攻角飛行鈍頭體流場Fig.5 Flow field of blunt body with large angle of attack

    表4 T1=226 K, Ma1=10時的斜激波關(guān)系Table 4 Oblique shock relations at T1=226 K, Ma1=10

    Gasmodelp21T21ρ21β/(°)u21Ma2WeakshockwavePerfect45.098.485.3138.520.7912.72Vibration43.847.415.9137.360.8013.05Difference/%-2.77-12.6111.30-3.011.3012.13StrongshockwavePerfect114.3920.045.7182.290.2190.49Vibration119.6017.176.9784.060.1760.44Difference/%4.55-14.3222.062.15-19.63-10.20

    3 振動激發(fā)下的激波反射

    定常激波反射分為規(guī)則反射(RR)和馬赫反射(MR),可以借助激波極線進(jìn)行分析,如圖6所示,圖中θa、θb、θc、θⅠ、θⅡ?yàn)榕c極線a、b、c、Ⅰ、Ⅱ?qū)?yīng)的氣流偏轉(zhuǎn)角。當(dāng)反射激波極線位于極線Ⅰ右側(cè)時,即θ>θⅠ時,可以形成馬赫反射,比如極線c。極線Ⅰ與縱軸以及入射激波極線共同交于一點(diǎn),對應(yīng)馬赫反射向規(guī)則反射轉(zhuǎn)變的von Neumann準(zhǔn)則。當(dāng)反射激波極線位于極線Ⅱ的左側(cè)時,即θ<θⅡ時,可以形成規(guī)則反射,比如極線a。極線Ⅱ與縱軸相切,對應(yīng)規(guī)則反射向馬赫反射轉(zhuǎn)變的脫體準(zhǔn)則。當(dāng)反射激波極線位于極線Ⅰ與極線Ⅱ之間時,即θⅠ<θ<θⅡ時,規(guī)則反射和馬赫反射都可能形成,即出現(xiàn)雙解區(qū)。

    圖6 規(guī)則反射和馬赫反射的激波極線 Fig.6 Shock polars of regular reflection and Mach reflection

    3.1 規(guī)則反射

    當(dāng)氣流偏轉(zhuǎn)角θ小于脫體準(zhǔn)則點(diǎn)θⅡ時,可以形成如圖1(a)所示的規(guī)則反射,此時反射激波極線與縱軸有2個交點(diǎn),即反射激波也存在強(qiáng)解與弱解。在實(shí)際中一般都是出現(xiàn)弱激波解,即激波反射的解一般取壓力激波極線中與縱軸相交壓力較小的點(diǎn),在速度激波極線中則是取速度較大的點(diǎn)。圖7為量熱完全氣體與考慮振動激發(fā)的壓力激波極線和速度激波極線。表5為反射激波后流場參數(shù)的對比,表中p31、T31、ρ31、u31分別對應(yīng)圖7 中③區(qū)和①區(qū)氣流的壓力比、溫度比、密度比和速度比,β2-3為反射激波與②區(qū)氣流的夾角。來流溫度T1=226 K,馬赫數(shù)Ma1=10,氣流偏轉(zhuǎn)角θ=30°。

    從圖7可見,無論是入射激波還是反射激波,考慮振動激發(fā)后,激波極線與量熱完全氣體的相比其整體輪廓都變大了。從反射激波極線的差異可見,振動激發(fā)帶來的激波極線的差異經(jīng)過一次激波反射后被放大了。即經(jīng)過一次激波反射,振動激發(fā)對波后流場的影響變得更大。結(jié)合表4和表5可得,振動激發(fā)對激波后壓力的影響,從入射激波的-2.77%被放大到反射激波的-2.92%,影響大小及其變化對壓力而言不是很明顯。然而對于激波后的速度,振動激發(fā)的影響從入射激波的1.30%被放大到反射激波的9.00%,變化及影響大小都非常明顯。另外還可算出反射激波后的溫度、密度、馬赫數(shù)和反射激波角,以及振動激發(fā)的影響大小。結(jié)果表明,從入射激波到反射激波,振動激發(fā)對激波后的溫度、密度、馬赫數(shù)和激波角等流場參數(shù)的影響都被明顯放大。

    圖7 規(guī)則反射的激波極線分析Fig.7 Shock polar analysis of regular reflection

    表5 T1=226 K, Ma1=10, θ=30° 時規(guī)則反射的激波關(guān)系Table 5 Shock relations of regular reflection at T1=226 K, Ma1=10, θ=30°

    Gasmodelp31T31ρ31β2-3/(°)u31Ma3Perfect261.3316.2116.1256.340.4891.22Vibration253.7012.7019.9748.500.5331.63Difference/%-2.92-21.6523.88-13.929.0033.61

    3.2 馬赫反射

    當(dāng)氣流偏轉(zhuǎn)角θ大于von Neumann準(zhǔn)則點(diǎn)θⅠ時,可以形成如圖1(b)所示的馬赫反射。此時在壓力激波極線圖上,反射激波極線與入射激波極線的強(qiáng)解部分相交,交點(diǎn)對應(yīng)的轉(zhuǎn)角和壓力就是馬赫桿后在三波點(diǎn)附近以及反射激波后的氣流偏轉(zhuǎn)角和壓力。而在速度激波極線圖上,相應(yīng)的氣流偏轉(zhuǎn)角對應(yīng)反射激波極線上的2個點(diǎn),其中速度較大的點(diǎn)3對應(yīng)反射激波后的狀態(tài),速度較小的點(diǎn)4對應(yīng)馬赫桿后在三波點(diǎn)附近的狀態(tài)。圖8 為量熱完全氣體與考慮振動激發(fā)的壓力激波極線和速度激波極線,表6為反射激波后滑移線兩側(cè)流場參數(shù)的對比,表中θ3為反射激波的氣流偏轉(zhuǎn)角;p31、u31、ρ31、T31對應(yīng)圖8中③區(qū)與①區(qū)的壓力比、速度比、密度比和溫度比,u41、ρ41、T41的含義與之類似;Δu34,1為u41和u31的差值,Δρ34,1和ΔT34,1的含義與之類似。來流溫度T1=226 K,馬赫數(shù)Ma1=10,氣流偏轉(zhuǎn)角θ=40°。

    圖8 馬赫反射的激波極線分析Fig.8 Shock polar analysis of Mach reflection

    可見,與規(guī)則反射一樣,從入射激波到反射激波,振動激發(fā)對流場參數(shù)的影響被明顯放大。然而馬赫桿是單激波,振動激發(fā)對馬赫桿后流場參數(shù)的影響并沒有被放大。由此,振動激發(fā)會對滑移線兩側(cè)流場參數(shù)的差異產(chǎn)生較大影響。從表6可見,振動激發(fā)對滑移線兩側(cè)速度差異的影響為11.76%,而對滑移線兩側(cè)密度差異的影響達(dá)51.81%?;凭€兩側(cè)的速度差異與KH(Kelvin-Helmholtz)不穩(wěn)定性直接相關(guān),考慮振動激發(fā)后滑移線兩側(cè)的速度差異變得更大,這說明振動激發(fā)對流動穩(wěn)定性也會產(chǎn)生一定影響。另外,振動激發(fā)對三波點(diǎn)后的氣流偏轉(zhuǎn)角產(chǎn)生了24.24%的影響。而對三波點(diǎn)后壓力的影響,與規(guī)則反射不同,考慮振動激發(fā)后壓力變大,影響大小為4.93%。

    表6 T1=226 K, Ma1=10, θ=40°時馬赫反射的激波關(guān)系Table 6 Shock relations of Mach reflection at T1=226 K, Ma1=10, θ=40°

    4 振動激發(fā)對激波反射轉(zhuǎn)變的影響

    如前所述,可以通過控制氣流偏轉(zhuǎn)角θ來實(shí)現(xiàn)定常激波反射在2種反射結(jié)構(gòu)之間轉(zhuǎn)變,且這種轉(zhuǎn)變存在遲滯現(xiàn)象。即理論上僅當(dāng)氣流偏轉(zhuǎn)角θ大于脫體準(zhǔn)則點(diǎn)θⅡ時,規(guī)則反射才向馬赫反射轉(zhuǎn)變;僅當(dāng)氣流偏轉(zhuǎn)角θ小于von Neumann準(zhǔn)則點(diǎn)θⅠ時,馬赫反射才向規(guī)則反射轉(zhuǎn)變,遲滯區(qū)為θⅠ<θ<θⅡ。

    圖4中指出,考慮振動激發(fā)后最大氣流偏轉(zhuǎn)角與量熱完全氣體的相比明顯增大。而從圖7和圖8可見,振動激發(fā)引起的最大氣流偏轉(zhuǎn)角的差異在經(jīng)過一次激波反射后被明顯放大,這將對規(guī)則反射向馬赫反射轉(zhuǎn)變的脫體準(zhǔn)則點(diǎn)的預(yù)測產(chǎn)生重大影響。圖9為激波反射轉(zhuǎn)變的2個準(zhǔn)則點(diǎn)以及對應(yīng)的激波角和氣流偏轉(zhuǎn)角隨來流馬赫數(shù)的變化。表7為2個激波反射轉(zhuǎn)變準(zhǔn)則點(diǎn)在2種氣體模型下的具體差異。來流溫度T1=226 K,馬赫數(shù)Ma1=10。

    可見,激波反射轉(zhuǎn)變的遲滯現(xiàn)象在馬赫數(shù)Ma1≈2.2時開始出現(xiàn),且隨著馬赫數(shù)的增大,遲滯區(qū)逐漸變寬??紤]振動激發(fā)后,激波反射轉(zhuǎn)變的2個準(zhǔn)則點(diǎn)θⅠ、θⅡ都相應(yīng)地變大。且脫體準(zhǔn)則點(diǎn)變大較明顯,而von Neumann準(zhǔn)則點(diǎn)變化較小,這導(dǎo)致激波反射轉(zhuǎn)變的遲滯區(qū)明顯變寬。從表7的數(shù)據(jù)來看,在馬赫數(shù)Ma1=10時,考慮振動激發(fā)后von Neumann準(zhǔn)則點(diǎn)θⅠ的差異在0.5°以內(nèi),影響大小為2.20%;而對于脫體準(zhǔn)則點(diǎn)θⅡ,振動激發(fā)的影響較大,角度差異約為3°,影響大小為9.59%。脫體準(zhǔn)則點(diǎn)θⅡ約3°的差異,對于認(rèn)識存在振動激發(fā)的高超聲速流動中的激波反射轉(zhuǎn)變現(xiàn)象,以及在試驗(yàn)中觀察和在工程中控制馬赫反射具有重大意義。

    圖9 規(guī)則反射和馬赫反射之間的轉(zhuǎn)變準(zhǔn)則Fig.9 Transition criterion between regular reflection and Mach reflection

    表7T1=226K,Ma1=10時規(guī)則反射和馬赫反射之間的臨界轉(zhuǎn)變角

    Table7CriticaltransitionanglesbetweenregularreflectionandMachreflectionatT1=226K,Ma1=10

    GasmodelθⅠ/(°)θⅡ/(°)Perfect19.5230.91Vibration19.9533.87Difference/%2.209.59

    5 結(jié) 論

    1) 研究結(jié)果表明,振動激發(fā)使激波極線的整體輪廓變大,對強(qiáng)激波和弱激波具有不同的影響規(guī)律。振動激發(fā)使強(qiáng)激波波后壓力變大,而使弱激波波后壓力變小。振動激發(fā)對激波流場的影響在經(jīng)過一次激波反射之后被明顯放大,分別對規(guī)則反射和馬赫反射的流場參數(shù)產(chǎn)生重大影響。

    2) 對于激波反射類型的轉(zhuǎn)變,振動激發(fā)使轉(zhuǎn)變的2個準(zhǔn)則點(diǎn)都變大,遲滯區(qū)變寬,且對規(guī)則反射向馬赫反射轉(zhuǎn)變的脫體準(zhǔn)則影響更大。對認(rèn)識存在振動激發(fā)的高超聲速流動中的激波反射轉(zhuǎn)變現(xiàn)象,以及在試驗(yàn)中觀察和在工程中控制馬赫反射具有重大意義。

    [1] BEN-DOR G. Shock wave reflection phenomena[M]. 2nd ed. Berlin: Springer-Verlag Press, 2007.

    [2] LI H, BEN-DOR G. Analysis of double-Mach-reflection wave configuration with convexly curved Mach stems[J]. Shock Waves, 1999, 9(5): 319-326.

    [3] VASILEY E I, BEN-DOR G, ELPERIN T, et al. The wall-jetting effect in Mach reflection: Navier Stokes simulations[J]. Journal of Fluid Mechanics, 2004, 511: 363-379.

    [4] TESDALL A M, HUNTER J K. Self-similar solutions for weak shock reflection[J]. SIAM Journal on Applied Mathematics, 2002, 63(1): 42-61.

    [5] SKEWS B W, LI G, PATON R. Experiments on Guderley Mach reflection[J]. Shock Waves, 2009, 19(2): 95-102.

    [6] MORIOKA T, SUZUKI Y, HONMA H. Radiation observation of strong shock wave reflection in air[C]//22th International Symposium on Shock Waves. London: Imperial College, 2000: 1201-1206.

    [7] BEN-DOR G. A state-of-the-knowledge review on pseudo-steady shock-wave reflections and their transition criteria[J]. Shock Waves, 2006, 15(3): 277-294.

    [8] VASILEY E I, ELPERIN T, BEN-DOR G. Analytical reconsideration of the von Neumann paradox in the reflection of a shock wave over a wedge[J]. Physics of Fluids, 2008, 20(4): 819-299.

    [9] VON NEUMANN J. Oblique reflection of shocks: Explosion research report 12[R]. Washington, D.C.: Navy Department, 1943.

    [10] VON NEUMANN J. Refraction, intersection and reflection of shock waves: NAVORD Report 203-45[R]. Washington, D.C.: Navy Department, 1943.

    [11] ANDERSON J D. Hypersonic and high-temperature gas dynamics[M]. 2nd ed. Reston, VA: AIAA, 2006: 385-386.

    [12] KIM H D, SETOGUCHI T. Shock induced boundary layer separation[C]//8th International Symposium on Experimental and Computational Aerothermodynamics of Internal Flows, 2007.

    [13] HEISER W, PRATT D, DALEY D, et al. Hypersonic airbreathing propulsion[M]. Reston, VA: AIAA, 1994.

    [14] GOONKO Y P, LATYPOV A F, MAZHUL I I, et al. Structure of flow over a hypersonic inlet with side compression wedges[J]. AIAA Journal, 2015, 41(3): 436-447.

    [15] KLEINE H, SETTLES G S. The art of shock waves and their flowfields[J]. Shock Waves, 2008, 17(5): 291-307.

    [16] JIANG Z L, TAKAYAMA K. Reflection and focusing of toroidal shock waves from coaxial annular shock tubes[J]. Computers and Fluids, 2013, 27(5-6): 553-562.

    [17] JACKSON S, GRUNTHANER M, SHEPHERD J. Wave implosion as an initiation mechanism for pulse detonation engines: AIAA-2003-4820[R]. Reston, VA: AIAA, 2003.

    [18] IZUMI K, ASO S, NISHIDA M. Experimental and computational studies focusing process of shock waves reflected from parabolic reflectors[J]. Shock Waves, 1994, 3(3): 213-222.

    [19] 滕宏輝, 張德良, 李輝煌, 等. 用環(huán)形激波聚焦實(shí)現(xiàn)爆轟波直接起爆的數(shù)值模擬[J]. 爆炸與沖擊, 2005, 25(6): 512-518.

    TENG H H, ZHANG D L, LI H H, et al. Numerical investigation of detonation direct initiation induced by toroidal shock wave focusing[J]. Explosion and Shock Waves, 2005, 25(6): 512-518 (in Chinese).

    [20] 姜宗林. 觸摸高溫氣體動力學(xué)[J]. 力學(xué)與實(shí)踐, 2006, 28(5): 1-7.

    JIANG Z L. Feeling high temperature gas dynamics[J]. Mechanics in Engineering, 2006, 28(5): 1-7 (in Chinese).

    [21] 樊菁. 高超聲速高溫氣體效應(yīng)判據(jù)[J]. 力學(xué)學(xué)報(bào), 2010, 42(4): 591-596.

    FAN J. Criteria on high-temperature gas effects around hypersonic vehicles[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(4): 591-596 (in Chinese).

    [22] 楊旸, 姜宗林, 胡宗民. 激波反射現(xiàn)象的研究進(jìn)展[J]. 力學(xué)進(jìn)展, 2012, 42(2): 141-161.

    YANG Y, JIANG Z L, HU Z M. Advances in shock wave reflection phenomena[J]. Advances in Mechanics, 2012, 42(2): 141-161 (in Chinese).

    [23] TARNAVSKII G A. Influence of flow angularities in a hypersonic ramjet diffuser of the formation of the shock-wave structure of the real gas flow[J]. Journal of Engineering Physics and Thermophysics, 2004, 77(3): 651-662.

    [24] 高云亮, 李進(jìn)平, 胡宗民, 等. 準(zhǔn)定常強(qiáng)激波馬赫反射波形的數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報(bào), 2008, 26(4): 456-461.

    GAO Y L, LI J P, HU Z M, et al. A numerical investigation on the Mach reflection patterns of quasi-steady strong shock waves[J]. Acta Aerodynamica Sinica, 2008, 26(4): 456-461 (in Chinese).

    [25] 李季. 高溫非平衡效應(yīng)下的激波干擾與激波反射[D]. 合肥: 中國科學(xué)技術(shù)大學(xué), 2015.

    LI J. On shock interactions and reflections with high temperature non-equilibrium effects[D]. Hefei: University of Science and Technology of China, 2015 (in Chinese).

    [26] MACH E. Uber den verlauf von funkenwellen in der ebene und im raume[J]. Sitzugsbr Akad Wiss Wien, 1878, 78: 819-838.

    [27] 汪志誠. 熱力學(xué)·統(tǒng)計(jì)物理[M]. 第四版. 北京: 高等教育出版社, 2008: 208.

    WANG Z C. Thermodynamics and statistical physics[M]. 4th ed. Beijing: Higher Education Press, 2008: 208 (in Chinese).

    [28] KAWAMURA R, SAITO H. Reflection of shock waves-1: Pseudo-stationary case[J]. Journal of Physical Society Japan, 1956, 11(5): 584-592.

    [29] KIM K H, KIM C, RHO O H. Methods for the accurate computations of hypersonic flows: I. AUSMPW+ scheme[J]. Journal of Computational Physics, 2001, 174(1): 81-119.

    [30] KIM K H, LEE J H, RHO O H. An improvement of AUSM schemes by introducing the pressure-based weight functions[J]. Computers & Fluids, 1998, 27(3): 311-346.

    [31] YOON S, JAMESON A. Lower-upper Symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations[J]. AIAA Journal, 2012, 26(9):1025-1026.

    [32] LEER B V. Towards the ultimate conservative difference scheme: V. A second-order sequel to Godunov’s method[J]. Journal of Computational Physics, 1979, 32(1): 101-136.

    (責(zé)任編輯: 王嬌)

    *Correspondingauthor.E-mail:huzm@imech.ac.cn

    Effectsofvibrationexcitationonshockreflections

    PENGJun1,2,ZHANGZijian1,2,ZHOUKai1,2,HUZongmin1,2,*,JIANGZonglin1,2

    1.StateKeyLaboratoryofHigh-temperatureGasDynamics,InstituteofMechanics,ChineseAcademyofSciences,Beijing100190,China2.SchoolofEngineeringScience,UniversityofChineseAcademyofSciences,Beijing100049,China

    SteadyshockreflectionsincluderegularreflectionandMachreflection,whichcantransittoeachotherundercriticalconditions.High-temperaturegaseffectisinevitableinhypersonicshockreflections.Astemperatureincreases,thevibrationexcitationofairmoleculescomesfirst.TheoreticalanalysisandquantitativecalculationareconductedtostudytheeffectsofvibrationexcitationonshockreflectionsandthetransitionsbetweenregularreflectionandMachrefection.Athermodynamicmodelforairwithvibrationexcitationispresentedandthencomparedwiththecaloricallyperfectgasmodel.Theinfluencesofvibrationexcitationonshockrelations,ontheflowfieldsinregularreflectionandMachreflection,andonthetransitioncriteriabetweenthemareanalyzed.Theresultsshowthatvibrationexcitationmayenlargetheoverallprofileoftheshockpolarascomparedwiththeshockpolarinthecaloricallyperfectgas.Inaddition,thedifferenceintheoverallpolarprofilesisamplifiedsignificantlyforthereflectedshock,andmayalterthereflectionconfiguration.Regardingtheshockreflectiontransitioncriteria,vibrationexcitationmaycauseincreasesofbothtransitionangles,i.e.,thedetachmentcriterionandvonNeumanncriterion,andtheincrementoftheformerismuchlargerthanthelatter.

    shockreflection;vibrationexcitation;Machreflection;transitioncriterion;hypersonicflow

    2016-12-15;Revised2017-01-04;Accepted2017-02-24;Publishedonline2017-03-231750

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

    s:NationalNaturalScienceFoundationofChina(11672308,11532014)

    2016-12-15;退修日期2017-01-04;錄用日期2017-02-24; < class="emphasis_bold">網(wǎng)絡(luò)出版時間

    時間:2017-03-231750

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

    國家自然科學(xué)基金 (11672308,11532014)

    .E-mailhuzm@imech.ac.cn

    彭俊, 張子健, 周凱, 等. 振動激發(fā)對激波反射的影響J. 航空學(xué)報(bào),2017,38(8):121055.PENGJ,ZHANGZJ,ZHOUK,etal.EffectsofvibrationexcitationonshockreflectionsJ.ActaAeronauticaetAstronauticaSinica,2017,38(8):121055.

    http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

    10.7527/S1000-6893.2017.121055

    V211.22

    A

    1000-6893(2017)08-121055-11

    猜你喜歡
    極線馬赫馬赫數(shù)
    東風(fēng)風(fēng)行T5馬赫版
    汽車觀察(2022年12期)2023-01-17 02:19:58
    一維非等熵可壓縮微極流體的低馬赫數(shù)極限
    破解定值有妙法,極點(diǎn)極線顯神威
    載荷分布對可控?cái)U(kuò)散葉型性能的影響
    穿越“馬赫谷”
    一道高考試題的背景簡介
    27馬赫,刺破蒼穹
    馬赫波反射中過度壓縮系數(shù)的計(jì)算
    NF-6連續(xù)式跨聲速風(fēng)洞馬赫數(shù)控制方式比較與研究
    半柔壁噴管初步實(shí)驗(yàn)研究
    毛片一级片免费看久久久久| 午夜激情福利司机影院| 高清毛片免费看| 少妇人妻精品综合一区二区| 丰满乱子伦码专区| 国产视频首页在线观看| 亚洲自偷自拍三级| 国产一区亚洲一区在线观看| 伦理电影免费视频| 男女边摸边吃奶| 亚洲精品国产av蜜桃| 欧美成人a在线观看| 久久久久久久久大av| 亚洲国产日韩一区二区| 亚洲丝袜综合中文字幕| 欧美人与善性xxx| 99视频精品全部免费 在线| 永久免费av网站大全| h日本视频在线播放| 最近的中文字幕免费完整| 精品一区在线观看国产| 国产精品人妻久久久影院| 午夜福利影视在线免费观看| 人妻少妇偷人精品九色| 亚洲精品成人av观看孕妇| 人妻系列 视频| 晚上一个人看的免费电影| 高清毛片免费看| 国产日韩欧美亚洲二区| 精品亚洲乱码少妇综合久久| 妹子高潮喷水视频| 插阴视频在线观看视频| 精品一区二区三区视频在线| 国产在视频线精品| 插阴视频在线观看视频| 六月丁香七月| av卡一久久| 夜夜骑夜夜射夜夜干| 久久国内精品自在自线图片| 性色avwww在线观看| 天堂8中文在线网| 亚洲高清免费不卡视频| 最近手机中文字幕大全| 高清欧美精品videossex| 国产亚洲精品久久久com| 久久6这里有精品| 色吧在线观看| 国产永久视频网站| 国产免费一级a男人的天堂| av.在线天堂| 免费在线观看成人毛片| 一级毛片 在线播放| www.色视频.com| 高清日韩中文字幕在线| 日韩av免费高清视频| 丰满乱子伦码专区| 国产男人的电影天堂91| 久久久久人妻精品一区果冻| 久久午夜福利片| 国产精品熟女久久久久浪| 国产精品99久久99久久久不卡 | 国产人妻一区二区三区在| 综合色丁香网| 国产高清有码在线观看视频| 啦啦啦视频在线资源免费观看| 日本av手机在线免费观看| 综合色丁香网| 一二三四中文在线观看免费高清| 亚洲成人手机| 亚洲四区av| 毛片一级片免费看久久久久| 插阴视频在线观看视频| 国产精品久久久久久精品电影小说 | 日本vs欧美在线观看视频 | 能在线免费看毛片的网站| 国产一区亚洲一区在线观看| 天天躁夜夜躁狠狠久久av| 日日摸夜夜添夜夜添av毛片| 国产色爽女视频免费观看| 亚洲欧美一区二区三区国产| av国产精品久久久久影院| 亚洲一级一片aⅴ在线观看| 男女啪啪激烈高潮av片| 欧美三级亚洲精品| 黄色配什么色好看| 色5月婷婷丁香| 精品视频人人做人人爽| 国产久久久一区二区三区| 综合色丁香网| 国产一区二区三区av在线| 三级国产精品欧美在线观看| 深夜a级毛片| 高清av免费在线| 国产成人精品婷婷| 亚洲欧洲国产日韩| 国产男人的电影天堂91| 国产高清有码在线观看视频| 久久久久网色| 如何舔出高潮| av卡一久久| av在线app专区| 久久久久久伊人网av| 久久婷婷青草| 日韩成人伦理影院| 久久久久国产网址| 国产成人精品福利久久| 99精国产麻豆久久婷婷| 一个人看的www免费观看视频| 亚洲av免费高清在线观看| 91精品伊人久久大香线蕉| 国产精品精品国产色婷婷| 美女xxoo啪啪120秒动态图| 嫩草影院入口| 啦啦啦啦在线视频资源| 偷拍熟女少妇极品色| 欧美区成人在线视频| 免费大片黄手机在线观看| 不卡视频在线观看欧美| 亚洲av中文av极速乱| 国产乱人视频| 国产一区二区在线观看日韩| 亚洲精品乱久久久久久| 哪个播放器可以免费观看大片| 国产永久视频网站| 狂野欧美激情性bbbbbb| 中文字幕精品免费在线观看视频 | 最后的刺客免费高清国语| 女人十人毛片免费观看3o分钟| 国产男女内射视频| 亚洲国产欧美人成| 一区二区三区四区激情视频| 在现免费观看毛片| 嘟嘟电影网在线观看| 国产高清三级在线| 欧美+日韩+精品| 制服丝袜香蕉在线| 一区在线观看完整版| 在线观看三级黄色| 色视频在线一区二区三区| 国产精品久久久久成人av| 蜜桃亚洲精品一区二区三区| 涩涩av久久男人的天堂| 在线看a的网站| 久久久久久人妻| 国产又色又爽无遮挡免| 最新中文字幕久久久久| 女性被躁到高潮视频| 在线观看免费视频网站a站| 午夜日本视频在线| 91精品国产国语对白视频| 18禁动态无遮挡网站| 大香蕉久久网| 久久人人爽人人爽人人片va| 熟女人妻精品中文字幕| 日韩国内少妇激情av| 日韩中字成人| 大香蕉97超碰在线| 少妇的逼水好多| 蜜臀久久99精品久久宅男| 欧美性感艳星| 又爽又黄a免费视频| 夜夜看夜夜爽夜夜摸| 午夜激情久久久久久久| 热re99久久精品国产66热6| 免费看光身美女| 精品国产乱码久久久久久小说| 99久久综合免费| 国产色爽女视频免费观看| 22中文网久久字幕| 国产av国产精品国产| 欧美高清性xxxxhd video| 欧美日韩综合久久久久久| 亚洲av.av天堂| 精品一区二区三区视频在线| 欧美激情极品国产一区二区三区 | 亚洲人成网站在线播| 精品国产露脸久久av麻豆| 在线观看三级黄色| 晚上一个人看的免费电影| 久久久久久久久久人人人人人人| 国产视频首页在线观看| 亚洲精品成人av观看孕妇| 尾随美女入室| 亚洲av二区三区四区| 国产亚洲一区二区精品| 国产91av在线免费观看| 久久久久久久久久久丰满| 黄色欧美视频在线观看| 夫妻性生交免费视频一级片| 成人免费观看视频高清| 日韩 亚洲 欧美在线| 97精品久久久久久久久久精品| 中文字幕亚洲精品专区| 久久国产精品男人的天堂亚洲 | 免费大片黄手机在线观看| 日产精品乱码卡一卡2卡三| 一边亲一边摸免费视频| 大片免费播放器 马上看| 国产精品蜜桃在线观看| 80岁老熟妇乱子伦牲交| 我的老师免费观看完整版| 最近最新中文字幕免费大全7| 久久97久久精品| 性色av一级| 人妻 亚洲 视频| 久久人妻熟女aⅴ| 美女国产视频在线观看| 精品久久久精品久久久| 91狼人影院| 婷婷色综合大香蕉| 丝袜脚勾引网站| 欧美人与善性xxx| 亚洲av.av天堂| 午夜激情福利司机影院| 黄色一级大片看看| a级毛片免费高清观看在线播放| 国产v大片淫在线免费观看| 大陆偷拍与自拍| 精品午夜福利在线看| 久久久久精品性色| 亚洲成色77777| 一级a做视频免费观看| 婷婷色av中文字幕| 国产高潮美女av| 在线观看国产h片| 成年av动漫网址| 交换朋友夫妻互换小说| videossex国产| 午夜免费观看性视频| 蜜臀久久99精品久久宅男| 美女福利国产在线 | 成人国产麻豆网| 国产精品嫩草影院av在线观看| 日日摸夜夜添夜夜添av毛片| 欧美精品国产亚洲| 久久久久久人妻| 国产高潮美女av| 亚洲,一卡二卡三卡| 亚洲人成网站高清观看| 看非洲黑人一级黄片| 男人爽女人下面视频在线观看| 日韩国内少妇激情av| 五月玫瑰六月丁香| 亚洲人成网站在线播| 午夜福利高清视频| 亚洲国产最新在线播放| 99国产精品免费福利视频| 十分钟在线观看高清视频www | 香蕉精品网在线| 久久人妻熟女aⅴ| 午夜免费鲁丝| 人人妻人人看人人澡| 色婷婷av一区二区三区视频| 国产精品蜜桃在线观看| 汤姆久久久久久久影院中文字幕| 国产精品熟女久久久久浪| 国产一区二区三区综合在线观看 | 夜夜爽夜夜爽视频| 99热国产这里只有精品6| 夜夜骑夜夜射夜夜干| 亚洲久久久国产精品| 亚洲精品乱久久久久久| 成年av动漫网址| 日韩av在线免费看完整版不卡| 欧美成人a在线观看| 国产又色又爽无遮挡免| 欧美精品国产亚洲| 中文字幕人妻熟人妻熟丝袜美| 国产黄色免费在线视频| 亚洲国产色片| 插逼视频在线观看| 91精品伊人久久大香线蕉| 久久毛片免费看一区二区三区| 美女国产视频在线观看| 久久热精品热| 激情五月婷婷亚洲| 日本欧美视频一区| www.色视频.com| 日本黄色片子视频| 日韩av不卡免费在线播放| 91久久精品国产一区二区成人| 欧美一级a爱片免费观看看| 亚洲成人中文字幕在线播放| 一级黄片播放器| 欧美日韩综合久久久久久| 美女xxoo啪啪120秒动态图| 国产日韩欧美在线精品| 七月丁香在线播放| 成年av动漫网址| 国产极品天堂在线| a 毛片基地| 久久女婷五月综合色啪小说| 亚洲欧美日韩卡通动漫| 国产精品99久久久久久久久| 欧美另类一区| 免费少妇av软件| 成人国产麻豆网| 国产亚洲一区二区精品| 午夜日本视频在线| 国产人妻一区二区三区在| 国产高清不卡午夜福利| 久久 成人 亚洲| 九九爱精品视频在线观看| 美女脱内裤让男人舔精品视频| 国产精品女同一区二区软件| 九九爱精品视频在线观看| 午夜福利在线观看免费完整高清在| 亚洲精品久久午夜乱码| 欧美老熟妇乱子伦牲交| 久久久精品94久久精品| 日韩电影二区| 啦啦啦在线观看免费高清www| 久久久久久久国产电影| 日本午夜av视频| 国产中年淑女户外野战色| 一级毛片黄色毛片免费观看视频| 大片免费播放器 马上看| 在线精品无人区一区二区三 | 亚洲国产精品999| 国产成人精品婷婷| 99热这里只有是精品50| 亚洲av男天堂| 国产在线一区二区三区精| 久久国产乱子免费精品| 欧美3d第一页| 久久99精品国语久久久| 少妇人妻精品综合一区二区| 中文字幕av成人在线电影| 久久国内精品自在自线图片| 精品国产露脸久久av麻豆| 国产av国产精品国产| 亚洲欧美成人精品一区二区| 亚洲精品国产成人久久av| 成人美女网站在线观看视频| 亚洲一级一片aⅴ在线观看| 国产爽快片一区二区三区| 街头女战士在线观看网站| 日本欧美国产在线视频| 美女脱内裤让男人舔精品视频| 女人十人毛片免费观看3o分钟| 欧美老熟妇乱子伦牲交| 我要看日韩黄色一级片| 亚洲欧美成人精品一区二区| 噜噜噜噜噜久久久久久91| 国产免费又黄又爽又色| 18禁裸乳无遮挡动漫免费视频| 亚洲不卡免费看| 欧美日韩精品成人综合77777| 青春草国产在线视频| av网站免费在线观看视频| 一区二区三区乱码不卡18| 女性生殖器流出的白浆| 一二三四中文在线观看免费高清| 免费大片黄手机在线观看| 一级毛片久久久久久久久女| 高清视频免费观看一区二区| 免费av中文字幕在线| 亚洲三级黄色毛片| 中文精品一卡2卡3卡4更新| 直男gayav资源| 精华霜和精华液先用哪个| 欧美丝袜亚洲另类| 亚洲精品日韩在线中文字幕| 亚洲欧美一区二区三区国产| 日韩av在线免费看完整版不卡| 成年人午夜在线观看视频| 成年美女黄网站色视频大全免费 | 国产精品女同一区二区软件| 麻豆成人av视频| 伦理电影免费视频| 在线播放无遮挡| a级毛片免费高清观看在线播放| 成年女人在线观看亚洲视频| 国产精品欧美亚洲77777| 成人美女网站在线观看视频| 午夜激情久久久久久久| 久久婷婷青草| 久久久精品94久久精品| 亚洲精华国产精华液的使用体验| 国产精品精品国产色婷婷| 免费看不卡的av| 97在线视频观看| 美女中出高潮动态图| 中文精品一卡2卡3卡4更新| 久久久久久久亚洲中文字幕| 国产成人精品婷婷| 一级毛片aaaaaa免费看小| 国产亚洲91精品色在线| 精品一品国产午夜福利视频| 精品久久国产蜜桃| 久久久欧美国产精品| 黑人猛操日本美女一级片| 免费人成在线观看视频色| 国产男人的电影天堂91| 国产欧美日韩精品一区二区| 精品久久国产蜜桃| 免费看不卡的av| 国产视频首页在线观看| 日日摸夜夜添夜夜添av毛片| 少妇裸体淫交视频免费看高清| 国产毛片在线视频| 天美传媒精品一区二区| 大香蕉久久网| 国模一区二区三区四区视频| 国产精品久久久久久久电影| 一区二区三区四区激情视频| 小蜜桃在线观看免费完整版高清| 激情五月婷婷亚洲| 久久av网站| 久久人人爽av亚洲精品天堂 | 久久韩国三级中文字幕| h日本视频在线播放| 大片免费播放器 马上看| 免费av不卡在线播放| 国产精品久久久久久av不卡| 网址你懂的国产日韩在线| 91精品国产国语对白视频| 国产精品.久久久| 一二三四中文在线观看免费高清| 国产男女超爽视频在线观看| 亚洲va在线va天堂va国产| 亚洲国产高清在线一区二区三| 夜夜骑夜夜射夜夜干| av.在线天堂| 色婷婷久久久亚洲欧美| 久久久久久九九精品二区国产| 国产精品99久久99久久久不卡 | 九草在线视频观看| 王馨瑶露胸无遮挡在线观看| 日本免费在线观看一区| 国产老妇伦熟女老妇高清| 男的添女的下面高潮视频| 狂野欧美激情性bbbbbb| 91aial.com中文字幕在线观看| 国产69精品久久久久777片| 中文欧美无线码| 如何舔出高潮| 少妇人妻久久综合中文| 国产色爽女视频免费观看| 亚洲av福利一区| 秋霞伦理黄片| 在线观看国产h片| 精品少妇黑人巨大在线播放| av卡一久久| 国产av码专区亚洲av| 夜夜骑夜夜射夜夜干| 日本黄大片高清| 熟女人妻精品中文字幕| 中文欧美无线码| 晚上一个人看的免费电影| 久久久a久久爽久久v久久| av又黄又爽大尺度在线免费看| 欧美一级a爱片免费观看看| 日本猛色少妇xxxxx猛交久久| 另类亚洲欧美激情| 日本av手机在线免费观看| 亚洲中文av在线| 午夜老司机福利剧场| 亚洲欧美精品自产自拍| 久久久久国产精品人妻一区二区| 久久99精品国语久久久| 中文字幕免费在线视频6| 久久久成人免费电影| 亚洲精品一区蜜桃| av又黄又爽大尺度在线免费看| 久久这里有精品视频免费| av一本久久久久| 国产毛片在线视频| 国产成人a∨麻豆精品| 亚洲成人av在线免费| 最近手机中文字幕大全| 中文字幕av成人在线电影| 天堂俺去俺来也www色官网| 少妇人妻一区二区三区视频| 欧美成人午夜免费资源| 国产欧美另类精品又又久久亚洲欧美| 80岁老熟妇乱子伦牲交| 天美传媒精品一区二区| 国产男女超爽视频在线观看| 麻豆成人午夜福利视频| av在线观看视频网站免费| 国产一级毛片在线| 免费播放大片免费观看视频在线观看| 少妇精品久久久久久久| 夫妻午夜视频| 亚洲在久久综合| 王馨瑶露胸无遮挡在线观看| 亚洲精品aⅴ在线观看| 亚洲精品第二区| 观看av在线不卡| 极品少妇高潮喷水抽搐| 国产亚洲午夜精品一区二区久久| 国产成人一区二区在线| 国产精品久久久久久av不卡| 我的老师免费观看完整版| 国产乱来视频区| 五月伊人婷婷丁香| 国产一区有黄有色的免费视频| 免费久久久久久久精品成人欧美视频 | 久久av网站| 美女福利国产在线 | 国产免费视频播放在线视频| 久久久久久久久久人人人人人人| 在线免费十八禁| 国产午夜精品久久久久久一区二区三区| 亚洲国产成人一精品久久久| 十分钟在线观看高清视频www | 亚洲久久久国产精品| 久久午夜福利片| 交换朋友夫妻互换小说| 国产精品久久久久久久久免| 亚洲av男天堂| 国产精品国产av在线观看| 校园人妻丝袜中文字幕| 欧美区成人在线视频| av女优亚洲男人天堂| 自拍偷自拍亚洲精品老妇| 亚洲精品乱码久久久v下载方式| 国产欧美另类精品又又久久亚洲欧美| 国产精品不卡视频一区二区| 美女内射精品一级片tv| 久久热精品热| 精品99又大又爽又粗少妇毛片| 各种免费的搞黄视频| 久久久欧美国产精品| 高清不卡的av网站| 成年美女黄网站色视频大全免费 | 亚洲精品国产av成人精品| 亚洲一区二区三区欧美精品| 欧美高清性xxxxhd video| 自拍欧美九色日韩亚洲蝌蚪91 | 日本与韩国留学比较| 熟妇人妻不卡中文字幕| 高清av免费在线| 免费看光身美女| 国产亚洲精品久久久com| 只有这里有精品99| 中文字幕人妻熟人妻熟丝袜美| 一级毛片电影观看| 三级经典国产精品| 人妻少妇偷人精品九色| 看非洲黑人一级黄片| 国产免费又黄又爽又色| 日本av手机在线免费观看| 亚洲电影在线观看av| 99热全是精品| 男女下面进入的视频免费午夜| 国产色婷婷99| 校园人妻丝袜中文字幕| 国产成人精品久久久久久| 性色av一级| 十分钟在线观看高清视频www | 九草在线视频观看| 国产片特级美女逼逼视频| 哪个播放器可以免费观看大片| 99热这里只有是精品50| 视频中文字幕在线观看| 亚洲第一区二区三区不卡| 热re99久久精品国产66热6| 日韩大片免费观看网站| 久久av网站| 美女cb高潮喷水在线观看| 日产精品乱码卡一卡2卡三| 日日啪夜夜爽| 中国国产av一级| h日本视频在线播放| 日韩 亚洲 欧美在线| 搡女人真爽免费视频火全软件| 亚洲丝袜综合中文字幕| 久久韩国三级中文字幕| 日韩av在线免费看完整版不卡| 丝瓜视频免费看黄片| 激情五月婷婷亚洲| 色网站视频免费| 国产视频内射| 爱豆传媒免费全集在线观看| 国产探花极品一区二区| 亚洲国产色片| 又大又黄又爽视频免费| 少妇的逼水好多| 欧美成人一区二区免费高清观看| 免费高清在线观看视频在线观看| 新久久久久国产一级毛片| 成人影院久久| 97在线视频观看| 亚洲精华国产精华液的使用体验| 久久影院123| 街头女战士在线观看网站| 欧美极品一区二区三区四区| 欧美日韩在线观看h| 久久99热6这里只有精品| 亚洲天堂av无毛| 欧美成人精品欧美一级黄| 国产午夜精品久久久久久一区二区三区| 99热这里只有是精品在线观看| av在线蜜桃| 少妇 在线观看| 国产黄频视频在线观看| 中文字幕免费在线视频6| 99久久精品热视频| 欧美成人午夜免费资源| 美女国产视频在线观看| tube8黄色片| 国产国拍精品亚洲av在线观看| 成人午夜精彩视频在线观看| 亚洲欧洲日产国产| 人妻一区二区av| 国产免费视频播放在线视频| 高清午夜精品一区二区三区| 在现免费观看毛片| 亚洲精品自拍成人| 欧美 日韩 精品 国产| 搡女人真爽免费视频火全软件| tube8黄色片| a级毛片免费高清观看在线播放|