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

    圓盤結(jié)構(gòu)下旋轉(zhuǎn)爆震波的不穩(wěn)定傳播特性

    2018-03-15 09:50:48夏鎮(zhèn)娟馬虎卓長飛周長省
    航空學報 2018年2期
    關(guān)鍵詞:爆震前導激波

    夏鎮(zhèn)娟,馬虎,卓長飛,周長省

    南京理工大學 機械工程學院,南京 210094

    爆震波是一種超聲速燃燒波,是激波與火焰面緊密耦合的聯(lián)合體,氣流跨過它后熱力學狀態(tài)(如壓力、溫度等)會急劇增加。爆震燃燒接近于等容燃燒,具有能量釋放率快、熱力循環(huán)效率高等優(yōu)點,與等壓燃燒相比,基于爆震燃燒的推進系統(tǒng)具有更高的熱效率,在航空航天領(lǐng)域有廣闊的應(yīng)用前景。

    旋轉(zhuǎn)爆震發(fā)動機(Rotating Detonation Engine,RDE)是基于爆震燃燒的一種新型動力推進裝置,一道或多道旋轉(zhuǎn)爆震波(Rotating Detonation Wave, RDW)在燃燒室內(nèi)連續(xù)傳播,爆震產(chǎn)物從另一端排出并產(chǎn)生連續(xù)推力。RDE的燃燒室構(gòu)型主要有3種形式:同軸圓環(huán)形、無內(nèi)柱圓筒形以及圓盤形結(jié)構(gòu)。目前實驗及數(shù)值模擬的研究工作主要集中在同軸圓環(huán)形結(jié)構(gòu),針對另外兩種結(jié)構(gòu)也有部分研究。

    20世紀60年代,Voitsekhovskii等[1]在圓盤形發(fā)動機上首次進行了旋轉(zhuǎn)爆震的實驗研究,得到了連續(xù)傳播的旋轉(zhuǎn)爆震波,并使用轉(zhuǎn)鼓式相機拍攝到旋轉(zhuǎn)爆震波的流場結(jié)構(gòu),開創(chuàng)了RDE實驗研究的先河。之后,眾多國家和機構(gòu)[1]也相繼開展了大量的實驗及數(shù)值研究。實驗方面,Bykovskii等[2-3]實現(xiàn)了多種氣態(tài)及液態(tài)燃料的旋轉(zhuǎn)爆震燃燒,并采用速度補償法記錄了旋轉(zhuǎn)爆震波的傳播過程,發(fā)現(xiàn)了兩種不同結(jié)構(gòu)的同向傳播模態(tài),以及爆震波的對撞傳播和軸向脈沖傳播模態(tài)。Anand等[4-5]進行了H2/Air在環(huán)形燃燒室內(nèi)的實驗研究,通過改變反應(yīng)物的質(zhì)量流率及發(fā)動機的尺寸,得到了4種不穩(wěn)定傳播模態(tài):非周期性無序傳播模態(tài)(重復熄滅、再起爆過程)、低頻正弦振蕩模態(tài)、多波轉(zhuǎn)變模態(tài)以及軸向脈沖傳播模態(tài)。Rankin等[6]使用OH*化學發(fā)光技術(shù)進行RDE的可視化實驗研究,觀測到燃燒室中爆震波的形狀和尺寸、其后的斜激波,以及新鮮填充氣體與燃燒產(chǎn)物之間的緩燃面等。Zhang等[7]在無內(nèi)柱空心圓筒燃燒室中進行了H2/Air旋轉(zhuǎn)爆震的實驗研究,得到了不同Laval噴管收縮比下爆震波的傳播模態(tài),并與切向不穩(wěn)定燃燒模態(tài)進行了對比。數(shù)值模擬方面,國內(nèi)外學者進行了大量的二維及三維數(shù)值研究,得到了旋轉(zhuǎn)爆震波的精細波系結(jié)構(gòu)及自持機理[8-9],Manabu等[10]揭示了旋轉(zhuǎn)爆震波“燃燒波+激波”的組合特性。Eude等[11]模擬了增加環(huán)形燃燒室寬度后的三維流場結(jié)構(gòu),發(fā)現(xiàn)當寬度變大后,燃燒室徑向尺寸的影響會越來越明顯,且在燃燒室徑向也觀察到了清晰的橫向爆震波結(jié)構(gòu)。Frolov等[12]開展了考慮黏性和推進劑分開噴注的三維數(shù)值模擬,發(fā)現(xiàn)了燃燒室中多個強度不同的爆震波頭,與實驗結(jié)果吻合。Tang等[13]對無內(nèi)柱的圓筒形發(fā)動機進行了三維數(shù)值模擬,得到了穩(wěn)定傳播的多波頭旋轉(zhuǎn)爆震波。Li等[14]數(shù)值研究了彎管中爆震波的傳播特性,得到了兩種傳播模態(tài):一種為解耦—再起爆模態(tài),隨著曲率半徑的減小,靠近內(nèi)圓附近的爆震波發(fā)生解耦,流場中的橫波傳播到內(nèi)圓附近,重新起爆爆震波;另一種為曲率半徑增大到一定值后,爆震波以穩(wěn)定模態(tài)傳播。

    Bykovskii等[15]在圓盤形發(fā)動機上進行了煤粉/空氣的兩相爆震實驗研究,實驗成功起爆并得到了脈沖爆震波和旋轉(zhuǎn)爆震波,并采用側(cè)壁面開窗的方式成功觀察到流場內(nèi)的爆震波結(jié)構(gòu)。Nakagami等[16]設(shè)計了兩側(cè)(或單側(cè))為石英玻璃壁面的圓盤形發(fā)動機,通過高速攝影及紋影技術(shù)觀測燃燒室全流場結(jié)構(gòu),便于研究爆震波的傳播特性。Ishiyama等[17-18]將這種圓盤形發(fā)動機與渦輪及壓氣機組合,設(shè)計了旋轉(zhuǎn)爆震渦輪組合發(fā)動機,并進行了冷流及燃燒實驗??梢暬瘜嶒灥膱A盤形RDE燃燒室示意圖如圖1所示,燃料與氧化劑從燃燒室外圓噴注進入燃燒室,產(chǎn)物從內(nèi)圓噴出,流道收斂,流通面積漸縮,這與傳統(tǒng)的圓環(huán)形燃燒室的等直通道有所差異。

    圖1 可視化實驗的燃燒室示意圖[16]Fig.1 Diagram of combustor used in visualization experiment[16]

    基于上述研究背景,忽略圓盤形燃燒室寬度的影響及出口流動的相互作用,將燃燒室簡化為二維環(huán)形收斂模型,噴注入口為計算域外圓,出口為計算域內(nèi)圓。由于計算域為漸縮通道,不同于傳統(tǒng)的RDE沿母線展開的二維矩形計算域,旋轉(zhuǎn)爆震波的傳播特性必然有所不同,有必要進行相關(guān)研究工作。本文開展了二維環(huán)形計算域的數(shù)值研究,得到兩種典型的旋轉(zhuǎn)爆震波傳播模態(tài),并對旋轉(zhuǎn)爆震波的不同傳播模態(tài)進行對比分析,重點研究不穩(wěn)定傳播模態(tài)的流場特征,以及爆震波參數(shù)、出口流場參數(shù)和出口增壓比的變化規(guī)律。

    1 數(shù)值方法

    1.1 求解方法

    針對2H2+O2+3.76N2的反應(yīng)混合物,本文進行二維圓盤結(jié)構(gòu)下旋轉(zhuǎn)爆震波傳播特性的數(shù)值模擬。忽略黏性、熱傳導和擴散等輸運效應(yīng),采用有限體積法離散二維化學非平衡流Euler控制方程組,采用重構(gòu)-推進方法計算無黏通量,對流項采用Steger-Warming矢通量分裂格式進行離散,使用時間算子分裂算法處理化學非平衡流中的剛性問題,化學反應(yīng)機理采用7組分8步基元反應(yīng),采用有限速率化學反應(yīng)模型,使用Arrhenius公式計算反應(yīng)速率常數(shù)。本文所用數(shù)值方法的詳細介紹及其驗證見文獻[19-20],這里不再贅述。

    1.2 計算模型與邊界條件

    本文所用計算模型為二維環(huán)形收斂計算域,如圖2所示,其中,Dout為計算域外徑,Din為計算域內(nèi)徑。

    圖2 二維環(huán)形計算域Fig.2 Two-dimensional annular computational domain

    計算域外圓邊界為噴注入口邊界,邊界上網(wǎng)格單元的流動情況由該處的壓力p決定,內(nèi)圓邊界為出口邊界。入口邊界條件為

    1)p≥p0,p0為預(yù)混氣體的噴注總壓,流動處于阻塞狀態(tài),噴注速度為0,預(yù)混氣體不能有效噴注進入燃燒室。

    2)p0>p>pcr,噴注速度按等熵膨脹公式計算:

    (1)

    式中:T0=300 K為噴注總溫;γ為反應(yīng)混合物的比熱比;R為混合物的氣體常數(shù)。pcr為臨界噴注壓力,其表達式為

    (2)

    3)p≤pcr,流動處于壅塞狀態(tài),噴注速度為

    (3)

    計算域的出口邊界為

    1) 出口為超聲速,所有守恒變量由內(nèi)部流場外推得到。

    2) 出口為亞聲速,邊界處壓力等于外界反壓,其他守恒變量由內(nèi)部流場參數(shù)外推得到,本文算例的外界反壓皆為0.1 MPa。

    點火區(qū)域為一扇形高溫高壓區(qū)域,如圖2所示,點火區(qū)壓力為2.5 MPa,溫度為3 000 K。預(yù)混氣體采用分段填充[21]的方式,第1個循環(huán)周期內(nèi),在點火段左側(cè)的扇形區(qū)域(如圖2所示)填充惰性氣體(N2),其他區(qū)域按入口邊界條件填充可燃預(yù)混氣體,以保證形成沿一個方向傳播的旋轉(zhuǎn)爆震波。

    本文的計算旨在研究二維環(huán)形計算域內(nèi)旋轉(zhuǎn)爆震波的不穩(wěn)定傳播特性及其參數(shù)的變化規(guī)律,不考慮旋轉(zhuǎn)爆震波的精細結(jié)構(gòu),綜合計算資源的限制,最終選擇1/3 mm的網(wǎng)格尺度為計算尺寸。網(wǎng)格無關(guān)性驗證已在之前的研究工作中完成[22]。

    2 結(jié)果與討論

    通過改變二維環(huán)形計算域的幾何尺寸,本文進行了圓盤結(jié)構(gòu)下的二維數(shù)值模擬工作,得到了旋轉(zhuǎn)爆震波的兩種典型傳播模態(tài):穩(wěn)定傳播模態(tài)和不穩(wěn)定傳播模態(tài)。其中,穩(wěn)定傳播模態(tài)在之前的研究工作中已進行了詳細的研究及分析[22],在此不再贅述,僅作為對比參考。本文著重研究不穩(wěn)定傳播模態(tài)下旋轉(zhuǎn)爆震波的流場特征、爆震波參數(shù)變化等。初步研究發(fā)現(xiàn),環(huán)形計算域的曲率半徑、寬度(B)以及噴注總壓對旋轉(zhuǎn)爆震波的傳播模態(tài)都有影響,計算工況及結(jié)果如表1所示。本文將選取典型工況2進行詳細分析。

    由之前的研究可得,工況1下,旋轉(zhuǎn)爆震波在二維環(huán)形計算域內(nèi)以穩(wěn)定模態(tài)傳播,傳播過程中,旋轉(zhuǎn)爆震波的波系結(jié)構(gòu)及爆震波強度皆能達到相對穩(wěn)定的狀態(tài),如圖3所示,其中t為計算時間。在本文中,工況1僅僅作為對照工況,而不再進行詳細分析。

    表1 計算工況及結(jié)果Table 1 Calculation cases and results

    圖3 RDW穩(wěn)定傳播模態(tài)的流場壓力云圖 (工況1,t:1.258~1.295 ms)Fig.3 Contours of pressure of RDW flow field in stable propagation mode (Case 1,t: 1.258-1.295 ms)

    與工況1相比,工況2中保持計算域的曲率半徑(內(nèi)圓半徑)不變,減小計算域的寬度。旋轉(zhuǎn)爆震波仍能成功起爆并沿周向傳播,但傳播過程并不穩(wěn)定,在傳播過程中重復出現(xiàn)“解耦—再起爆”的現(xiàn)象,下面將對這種不穩(wěn)定傳播模態(tài)進行詳細分析。

    2.1 不穩(wěn)定傳播流場云圖

    工況2下,旋轉(zhuǎn)爆震波起爆成功,但爆震波在傳播過程中出現(xiàn)了“解耦—再起爆”的現(xiàn)象。圖4、圖5為不穩(wěn)定傳播模態(tài)下不同時刻的流場云圖,圖4表示旋轉(zhuǎn)爆震波解耦的過程,圖5為爆震波重新起爆的過程,其中圖4(a)、圖5(a)為壓力云圖,圖4(b)、圖5(b)為溫度(T)云圖,圖4(c)為某一時刻H2質(zhì)量分數(shù)(CH2)的局部放大圖,圖5 中HS(Hot Spot)為局部熱點。圖6為不同時刻壓力云圖的局部放大圖,圖中:DW(Detonation Wave)為爆震波,RW(Reflected Wave)為反射激波,DLW(Decoupling Leading shock Wave)為解耦前導激波,TDW(Transverse Detonation Wave)為橫向爆震波。

    由圖4可得,旋轉(zhuǎn)爆震波在二維環(huán)形計算域內(nèi)傳播時,靠近內(nèi)圓出口位置,爆震波的前導激波與火焰鋒面逐漸發(fā)生分離,這主要是由于內(nèi)圓附近擴張幾何曲面的發(fā)散作用削弱了爆震波的強度,當燃燒產(chǎn)生的能量不能維持前導激波的傳播時,火焰鋒面與前導激波面逐漸分離,爆震波發(fā)生解耦。由溫度分布圖4(b)可得,工況2在圖4中的解耦開始時刻約為0.243 ms,且靠近出口的解耦區(qū)域會逐漸往流場內(nèi)部擴散,影響環(huán)形流場內(nèi)爆震波的傳播,但爆震波的解耦范圍是有限的,在該計算工況下,最大解耦區(qū)域的徑向尺寸占整個環(huán)形計算域?qū)挾鹊?/3左右。圖4(c)為t=0.257 5 ms時刻H2質(zhì)量分數(shù)云圖的局部放大圖,由圖4(b)、圖4(c)可以看出,在爆震波解耦區(qū)域,前導激波與其后的火焰鋒面之間存在一小塊未燃氣體區(qū)域(Unburned Zone,UZ)。

    圖4 RDW不穩(wěn)定傳播模態(tài)的流場云圖(解耦,工況2,t:0.233 33~0.262 49 ms)Fig.4 Contours of RDW flow field in unstable propagation mode (decoupling,Case 2,t: 0.233 33-0.262 49 ms)

    圖5 RDW不穩(wěn)定傳播模態(tài)的流場云圖(再起爆,工況2,t:0.265 66~0.308 87 ms)Fig.5 Contours of RDW flow field in unstable propagation mode (re-initiation,Case 2,t: 0.265 66-0.308 87 ms)

    圖6 不同時刻壓力云圖的局部放大圖Fig.6 Detail views of pressure contours at different moments

    由圖4(a)的壓力云圖可得,旋轉(zhuǎn)爆震波解耦之前,前導激波與其在外壁面的反射激波之間的夾角較大且維持在一定值范圍之內(nèi),這與穩(wěn)定傳播工況的壓力流場云圖一致,如圖3所示。但隨著出口附近爆震波的解耦,解耦處前導激波的傳播速度降低,逐漸落后于未解耦區(qū),而爆震波的解耦對反射激波的傳播影響較小,反射激波的傳播速度變化不大。因此,隨著時間的推移,前導激波與反射激波的夾角逐漸變小,兩者距離縮短,如圖6(a)所示,到圖4(a)壓力云圖的最后一個時刻,反射激波有趕上前導激波的趨勢,圖6(b)為該時刻壓力云圖的放大圖,由圖6(b)可得,爆震波的解耦區(qū)域約占整個環(huán)域?qū)挾鹊?/3,反射激波受解耦區(qū)域的影響,波形發(fā)生變化,但傳播速度仍比解耦的前導激波快,有追趕上前導激波的趨勢。

    結(jié)合圖4及圖6(b)的分析,前導激波與反射激波的距離逐漸縮短,最終反射激波趕上前導激波,兩道激波發(fā)生碰撞,如圖6(c)所示,碰撞點處的壓力、溫度驟然上升,遠高于流場其他位置,流場中出現(xiàn)局部熱點HS。局部熱點在未燃預(yù)混氣體中以半圓形波陣面快速膨脹傳播,如圖6(d)所示。由圖4(b)、圖4(c)可知,由于爆震波的解耦,流場出口附近存在一小塊未燃區(qū)域,重新起爆的半圓形爆震波在流場中傳播時,遇到該未燃氣體區(qū)域并迅速形成了一道沿徑向傳播的橫向爆震波TDW,如圖6(d)中實線圓圈所示。由于該區(qū)域的未燃氣體經(jīng)過前導激波壓縮,初始壓力較高,因此,在該區(qū)域形成的沿徑向傳播的爆震波的強度遠高于沿周向傳播的旋轉(zhuǎn)爆震波。

    經(jīng)過約5.8 μs,局部熱點形成的半圓形爆震波擴展到環(huán)形區(qū)域的整個徑向,在流場中形成了近似的平面爆震波,如圖6(e)、圖6(f)所示,右側(cè)的橫向爆震波傳播至內(nèi)圓附近時,由于內(nèi)圓邊界為出口無反射邊界條件,TDW直接傳出計算域而未發(fā)生反射。平面爆震波沿環(huán)形計算域傳播一段時間后,由于計算域幾何條件的影響,逐漸恢復成解耦之前的曲面爆震波形狀,其波后同樣形成了反射激波,如圖6(g)所示,這與解耦之前的流場極為相似。但重新起爆的爆震波并沒有穩(wěn)定傳播下去,而是在接下來旋轉(zhuǎn)爆震波沿周向傳播的整個過程中(該算例下旋轉(zhuǎn)爆震波傳播了12個周期)一直重復循環(huán)“解耦—再起爆”的過程,這由圖5后面時刻的流場圖可以看出,重新起爆的爆震波剛恢復解耦前的波形,在靠近出口區(qū)域又出現(xiàn)了解耦現(xiàn)象。

    2.2 流場參數(shù)變化

    為對比研究旋轉(zhuǎn)爆震波在不穩(wěn)定傳播模態(tài)與穩(wěn)定傳播模態(tài)下流場參數(shù)的差異,該節(jié)研究了不同徑向位置處旋轉(zhuǎn)爆震波的壓力、溫度、傳播速度等的變化,并與同一初始條件下的理論值進行對比,分析爆震波的不穩(wěn)定傳播對爆震參數(shù)的影響。

    2.2.1 爆震波壓力

    由之前的研究工作可得,在穩(wěn)定傳播模態(tài)下,爆震波流場內(nèi)各監(jiān)測點處的壓力值呈現(xiàn)周期性變化趨勢,且壓力峰值相對穩(wěn)定。而對于不穩(wěn)定傳播模態(tài),流場內(nèi)的壓力變化與爆震波的解耦再起爆時刻及位置有關(guān),圖7對比研究了兩種傳播模態(tài)下不同徑向位置處壓力隨時間的變化情況,其中各監(jiān)測點的布置如圖2所示,D為不同監(jiān)測點處的直徑。

    由圖7(a)可得,穩(wěn)定傳播模態(tài)下,同一監(jiān)測點處各周期之間旋轉(zhuǎn)爆震波的峰值壓力差別很小,說明爆震波在該工況下能實現(xiàn)穩(wěn)定傳播,但不同徑向位置的爆震波壓力值不同,這主要是由環(huán)形計算域的幾何條件所致,外圓壁面的收斂匯聚作用增強了爆震波強度,而內(nèi)圓出口附近的發(fā)散作用削弱了爆震波強度,使得越靠近內(nèi)圓出口,爆震波的強度及峰值壓力越小。

    由圖7(b)可得,不穩(wěn)定傳播模態(tài)下,旋轉(zhuǎn)爆震波不僅在不同徑向位置處壓力峰值不同,而且在同一監(jiān)測點的不同傳播周期內(nèi),爆震波的峰值壓力也存在較大差異,且越靠近計算域的出口,這種差異越大。這說明在工況2條件下,旋轉(zhuǎn)爆震波強度隨時間呈現(xiàn)強弱變化規(guī)律。這與2.1節(jié)的壓力流場云圖相對應(yīng),不穩(wěn)定傳播模態(tài)下旋轉(zhuǎn)爆震波重復出現(xiàn)“解耦—再起爆”的循環(huán)過程。如圖7(b)所示,虛線圓圈所示即為爆震波發(fā)生解耦時的壓力峰值,由于解耦的前導激波缺少化學反應(yīng)能量的支持,激波強度下降,壓力隨之下降。而其后一個壓力峰值卻遠大于流場中其他位置處的壓力,如圖中實線圓圈所示,說明此時流場中存在局部熱點形成的過驅(qū)爆震波,之后過驅(qū)爆震波逐漸衰減為旋轉(zhuǎn)爆震波。由圖中還可以發(fā)現(xiàn),越靠近外圓,不同周期的壓力峰值差異越小,這也與壓力云圖相對應(yīng),靠近外圓的爆震波強度較高,足以抵抗稀疏波的削弱,并未發(fā)生解耦,但仍然受到解耦部分流場的影響,爆震波壓力有小幅度波動。

    圖7 不同直徑處壓力隨時間的變化Fig.7 Variation of pressure with time in different diameters

    2.2.2 周向壓力、溫度分布

    為研究爆震波解耦和重新起爆時流場內(nèi)壓力、溫度等參數(shù)的變化情況,提取同一徑向位置、不同時刻的周向數(shù)據(jù),對比不同時刻流場中壓力、溫度的周向分布,詳細分析不穩(wěn)定傳播模態(tài)下爆震波解耦以及再起爆的過程。

    由前面的流場云圖分析已得,工況2條件下爆震波的解耦區(qū)域約占環(huán)形計算域?qū)挾鹊?/3,因此,選取直徑70 mm處監(jiān)測點的周向數(shù)據(jù)來研究爆震波的解耦—再起爆過程。

    圖8 不同時刻的壓力、溫度曲線(監(jiān)測點:D=70 mm)Fig.8 Pressure and temperature curves at different moments (monitoring point: D=70 mm)

    圖8選取了3個不同時刻周向的壓力、溫度數(shù)據(jù),分別對應(yīng)于爆震波未解耦、解耦以及再起爆階段。

    由圖8(a)可知,爆震波解耦前,前導激波與燃燒波緊密耦合,爆震波的壓力峰值為2.12 MPa,溫度峰值為2 171 K,用CEA(Chemical Equilibrium with Applications)計算理論壓力值為2.01 MPa,溫度為2 920.6 K,爆震速度為1 975.6 m/s。其中,計算初始條件取波前的壓力和溫度值(分別為82.782 kPa和190.65 K)。計算壓力值與理論值吻合較好,但溫度相差較大,這可能受計算域曲率的影響。圖8(b)為爆震波解耦時刻的曲線圖,從圖8(b)及其局部放大圖可得,溫度曲線上升的周向位置明顯落后于壓力曲線,說明此時的前導激波與燃燒波已經(jīng)分離,爆震波解耦,且該時刻下的壓力峰值僅為0.92 MPa,遠低于CEA理論計算值,另一個波峰代表反射激波,由曲線圖可得,反射激波的壓力峰值要高于解耦的前導激波,說明反射激波受爆震波解耦的影響較小,此時反射激波的強度高于解耦后的前導激波強度,這也為反射激波追趕上前導激波,重新起始爆震波提供條件。圖8(c)表示壓力和溫度的上升位置吻合良好,說明前導激波與火焰鋒面重新耦合在一起,且由圖8(c)可得,此時的壓力峰值高達3.50 MPa,高于理論計算值,這說明此處為局部熱點形成的過驅(qū)爆震波,過驅(qū)爆震波逐漸衰減后形成旋轉(zhuǎn)爆震波。

    2.2.3 傳播速度

    圖9表示工況2條件下,爆震波解耦過程中,不同徑向位置處前導激波(Leading Shock wave, LS)及反射激波RW的傳播速度隨時間的變化規(guī)律,其中時間軸與圖4流場云圖的時間一致。

    圖9 不同徑向位置前導激波與反射激波的傳播速度Fig.9 Propagation velocities of leading shock wave and reflected wave in different radial positions

    由圖9可得,爆震波解耦前,隨著直徑D的增大,前導激波與反射激波的傳播速度逐漸增大且保持一致,在直徑90 mm附近,爆震波傳播速度與理論值1 975.6 m/s接近,靠近內(nèi)圓出口的速度受幾何條件的限制,傳播速度較低。未解耦流場的速度分布與穩(wěn)定傳播工況下流場內(nèi)的速度分布相同,保證了爆震波與反射激波的相對位置不變,如圖3及圖4中未解耦時刻的流場云圖所示。對于不穩(wěn)定傳播模態(tài),隨著時間的推移,靠近出口的前導激波速度首先開始下降,隨之擴散到流場內(nèi)部。隨著直徑的增加,前導激波速度依次開始下降,這與圖4中壓力云圖的解耦順序一致,爆震波從出口位置開始解耦,并逐漸往流場內(nèi)部擴散。但由圖9可得,反射激波的傳播速度受解耦的影響較小,速度變化不大,且大于同一位置前導激波的傳播速度。因此,后面的反射激波逐漸追趕前導激波,最終發(fā)生碰撞促使流場中局部熱點的形成。

    由圖9發(fā)現(xiàn),外圓附近的前導激波的傳播速度幾乎沒有變化,兩道激波的傳播速度保持一致,這是因為外圓的收斂匯聚作用增強了爆震波強度,使之足以抵抗稀疏波的削弱而未發(fā)生解耦。

    2.3 出口參數(shù)變化

    由前所述,不穩(wěn)定傳播模態(tài)下,爆震波的“解耦—再起爆”過程對流場參數(shù)有很大的影響,該節(jié)重點研究爆震波的不穩(wěn)定傳播對出口速度、馬赫數(shù)以及增壓比的影響。

    2.3.1 出口速度及馬赫數(shù)

    圖10表示不同時刻下,計算域出口的徑向分速度沿周向的分布。其中,t=0.238 2、0.257 5、0.271 5、0.292 3 ms分別為從爆震波的穩(wěn)定段、解耦段、再起爆段和恢復穩(wěn)定段選取的時刻,以此研究同一工況條件下旋轉(zhuǎn)爆震波的不穩(wěn)定傳播對出口速度的影響。

    由圖10可得,在爆震波傳播的不同階段,出口徑向速度沿周向的分布大致相同,徑向速度的最大值在爆震波鋒面附近取得,且不同時刻的速度峰值變化不大。說明爆震波的“解耦—再起爆”過程對徑向速度分量影響很小。

    除了徑向速度的分布,本節(jié)還分析了出口平面的各速度分量及馬赫數(shù)的平均值隨時間的變化。圖11為兩個連續(xù)的“解耦—再起爆”周期內(nèi)出口速度分量及馬赫數(shù)平均值的變化曲線,其中,平均值為各參數(shù)在出口平面上沿周向的平均。

    由圖11可得,爆震波的不穩(wěn)定傳播對出口速度分量的平均值影響很小,在兩個“解耦—再起爆”周期內(nèi),爆震波的解耦及再起爆過程會引起出口速度的輕微波動,但基本保持不變。出口平均馬赫數(shù)的波動值較速度波動有所增加,但總體變化仍舊很小。由此說明,爆震波的不穩(wěn)定傳播對燃燒室出口的速度及馬赫數(shù)影響不大。

    圖10 不同時刻下出口徑向速度分量沿周向的分布Fig.10 Distribution of radial velocity component along circumference at exit for different moments

    圖11 出口馬赫數(shù)和速度分量隨時間的變化Fig.11 Variation of velocity components and Mach number at exit with time

    2.3.2 出口增壓比

    (4)

    式中:pout_total,t為t時刻燃燒室出口總壓;ρout為燃燒室出口密度;ds為出口平面上網(wǎng)格的周向長度;l為燃燒室出口的周長。

    (5)

    圖12 出口增壓比隨時間的變化Fig.12 Variation of outlet pressure amplification ratio with time

    圖12(a)為穩(wěn)定傳播模態(tài)下(工況1)燃燒室出口增壓比隨時間的變化曲線,由圖中曲線可得,隨著時間的推移,出口增壓比的波動逐漸減小,最后趨于穩(wěn)定,穩(wěn)定后增壓比處于動態(tài)平衡狀態(tài),呈現(xiàn)周期性的振蕩,但振幅很小。該工況條件下,燃燒室出口增壓比的平均值大約在1.5左右。

    圖12(b)為不穩(wěn)定傳播模態(tài)(工況2)下燃燒室出口增壓比隨時間的變化曲線,由圖可得,在開始時刻,出口增壓比變化較為紊亂,隨著時間的推移,增壓比并沒有趨于穩(wěn)定,而是呈周期性波動,波動的幅值較高且不穩(wěn)定。周期性波動的第1個峰值時刻(t=0.268 5 ms,圖12(b)中實線圓圈所示)對應(yīng)于爆震波第1個“解耦—再起爆”周期內(nèi)重新起爆階段的某一時刻,如圖13所示,此時重新起始的爆震波強度較高,屬于過驅(qū)爆震,且沿徑向傳播的橫向爆震波(圖6(d)中TDW)剛好傳至出口位置,因此,爆震波附近的出口總壓值迅速上升,增壓比隨之迅速躍升至峰值。隨著時間的推移,沿徑向傳播的TDW從出口傳出,且重新起始的過驅(qū)爆震波逐漸衰減為旋轉(zhuǎn)爆震波,出口總壓下降,增壓比迅速衰減并維持一段較為穩(wěn)定的值(如圖12(b)局部放大圖所示),直到下一個循環(huán)周期的爆震波再起爆階段,增壓比再次上升。

    圖13 工況2中壓力云圖的局部放大圖 (t=0.268 5 ms)Fig.13 Detail view of pressure contours in Case 2 (t=0.268 5 ms)

    由此可得,出口增壓比躍升的主要原因是:再起爆時,流場中的局部熱點遇到未燃氣體區(qū)域形成強度較高的沿徑向傳播的橫向爆震波TDW,TDW傳播至出口位置時引起出口總壓的急劇上升。由于計算域的內(nèi)圓邊界條件為無反射的出口邊界,TDW傳播至出口后直接傳出并沒有發(fā)生反射。因此,TDW從出口傳出的下一時刻,爆震波附近總壓迅速下降,增壓比相應(yīng)降低(見圖12(b))。重新起始的周向過驅(qū)爆震波也在增壓比的變化過程中起了一定作用。

    由上述分析可得,出口增壓比的循環(huán)周期應(yīng)與爆震波“解耦—再起爆”的循環(huán)周期保持一致。經(jīng)過計算,圖12(b)所示燃燒室出口增壓比的循環(huán)頻率如圖14所示。開始階段,頻率較低,隨著時間的推移,頻率逐漸增加,并在幾個循環(huán)過后趨于動態(tài)穩(wěn)定,頻率均值約為21.6 kHz,頻率值在均值附近上下波動,但波動幅值較小。另外,由解耦云圖可得,爆震波“解耦—再起爆”的循環(huán)頻率也有相似的變化,開始時頻率較低,約為17.4 kHz左右,隨著時間的推移,頻率逐漸上升,最終頻率值穩(wěn)定在21.6 kHz附近并上下波動,這也驗證了開始的推論:燃燒室出口增壓比的循環(huán)周期與爆震波“解耦—再起爆”的循環(huán)周期保持一致,在該算例下,兩者的循環(huán)頻率皆為21.6 kHz。對圖12(b)中動態(tài)穩(wěn)定后的出口增壓比取時間平均,得到不穩(wěn)定傳播工況2下出口增壓比的平均值約為1.5,與穩(wěn)定傳播工況1的出口增壓比數(shù)值相近。

    圖14 出口增壓比的頻率隨時間的變化Fig.14 Variation of frequency of pressure amplification ratio at exit with time

    3 結(jié) 論

    通過對二維圓盤結(jié)構(gòu)下旋轉(zhuǎn)爆震波傳播模態(tài)的數(shù)值計算,得到如下結(jié)論:

    1) 得到了旋轉(zhuǎn)爆震波不穩(wěn)定傳播模態(tài)的流場特點:不穩(wěn)定傳播模態(tài)下,旋轉(zhuǎn)爆震波重復進行“解耦—再起爆”過程,爆震波的壓力、溫度以及傳播速度皆隨爆震波的解耦及再起爆發(fā)生變化。

    2) 不穩(wěn)定傳播模態(tài)下,內(nèi)圓擴張曲面的發(fā)散作用,使前導激波與火焰鋒面分離,出口附近的爆震波首先解耦,解耦區(qū)域逐漸往流場內(nèi)部擴張;外圓附近的爆震波因幾何收斂作用增強而未發(fā)生解耦。爆震波解耦導致前導激波的傳播速度降低,但反射激波受解耦的影響較小,傳播速度大于前導激波,反射激波追趕上前導激波并發(fā)生碰撞,產(chǎn)生局部熱點,重新起始爆震波。

    3) 爆震波的不穩(wěn)定傳播對出口速度分量的平均值影響很小,爆震波的解耦及再起爆過程會引起出口速度的輕微波動,但基本保持不變。出口平均馬赫數(shù)的波動值較速度波動有所增加,但總體變化仍舊很小。出口增壓比受不穩(wěn)定傳播模態(tài)的影響較大,隨時間呈周期性變化,變化幅值較高且不穩(wěn)定。出口增壓比的循環(huán)周期與爆震波“解耦—再起爆”的循環(huán)周期保持一致,穩(wěn)定后兩者的循環(huán)頻率值皆為21.6 kHz。

    [1] VOITSEKHOVSKII B V, MITROFANOV V V, TOPCHIYAN M E. Structure of the detonation front in gases[J]. Combustion, Explosion, and Shock Waves, 1969, 5(3): 385-395 (in Russian).

    [2] BYKOVSKII F A, ZHDAN S A, VEDERNIKOV E F. Continuous spin detonations[J]. Journal of Propulsion and Power, 2006, 22 (6): 1204-1216.

    [3] BYKOVSKII F A, ZHDAN S A, VEDERNIKOV E F. Continuous detonation of syngas-air mixtures in an annular flow-type cylindrical combustor[C]∥24th International Colloquium on the Dynamics of Explosions and Reactive Systems, 2013: 1-4.

    [4] ANAND V, GEORGE A S, DRISCOLL R, et al. Characterization of instabilities in a rotating detonation combustor[J]. International Journal of Hydrogen Energy, 2015, 40(46): 16649-16659.

    [5] ANAND V, GEORGE A S, DRISCOLL R, et al. Longitudinal pulsed detonation instability in a rotating detonation combustor[J]. Experimental Thermal and Fluid Science, 2016, 77: 212-225.

    [6] RANKIN B A, RICHARDSON D R, CASWELL A W, et al. Chemiluminescence imaging of an optically accessible non-premixed rotating detonation engine[J]. Combustion and Flame, 2017, 176 (1): 12-22.

    [7] ZHANG H L, LIU W D, LIU S J. Experimental investigations on H2/air rotating detonation wave in the hollow chamber with Laval nozzle[J]. International Journal of Hydrogen Energy, 2017, 42(5): 3363-3370.

    [8] 劉世杰, 林志勇, 孫明波, 等. 旋轉(zhuǎn)爆震波發(fā)動機二維數(shù)值模擬[J]. 推進技術(shù), 2010, 31(5): 634-640.

    LIU S J, LIN Z Y, SUN M B, et al. Two-dimensional numerical simulation of rotating detonation wave engine [J]. Journal of Propulsion Technology, 2010, 31(5): 634-640 (in Chinese).

    [9] 劉世杰, 覃慧, 林志勇, 等. 連續(xù)旋轉(zhuǎn)爆震波細致結(jié)構(gòu)及自持機理[J]. 推進技術(shù), 2011, 32(3): 431-436.

    LIU S J, QIN H, LIN Z Y, et al. Detailed structure and propagating mechanism research on continuous rotating detonation wave[J]. Journal of Propulsion Technology, 2011, 32(3): 431-436 (in Chinese).

    [10] MANABU H, FUJIWARA T, PIOTR W. Fundamentals of rotating detonations[J]. Shock Waves, 2009, 19(1): 1-10.

    [11] EUDE Y, DAVIDENKO D M, G?KALP I, et al. Use of the adaptive mesh refinement for 3D simulations of a CDWRE (Continuous Detonation Wave Rocket Engine): AIAA-2011-2236[R]. Reston, VA: AIAA, 2011.

    [12] FROLOV S M, DUBROVSKII A V, IVANOV V S. Three-dimensional numerical simulation of the operation of a rotating-detonation chamber with separate supply of fuel and oxidizer[J]. Russian Journal of Physical Chemistry B, 2013, 7 (1): 35-43.

    [13] TANG X M, WANG J P, SHAO Y T. Three-dimensional numerical investigations of the rotating detonation engine with a hollow combustor[J]. Combustion and Flame, 2015, 162(4): 997-1008.

    [14] LI J, NING J G, ZHAO H, et al. Numerical investigation on the propagation mechanism of steady cellular detonations in curved channels[J]. Chinese Physics Letters, 2015, 32(4): 144-147.

    [15] BYKOVSKII F A, ZHDAN S A, VEDERNIKOV E F, et al. Detonation of a coal-air mixture with addition of hydrogen in plane-radial vortex chambers[J]. Combustion, Explosion, and Shock Waves, 2011, 47(4): 473-482.

    [16] NAKAGAMI S, MATSUOKA K, KASAHARA J, et al. Visualization of rotating detonation waves in a plane combustor with a cylindrical wall injector: AIAA-2015-0878[R]. Reston, VA: AIAA, 2015.

    [17] ISHIYAMA C, MIYAZAKI K, NAKAGAMI S, et al. Experimental study of research of centrifugal-compressor-radial-turbine type rotating detonation engine: AIAA-2016-5103[R]. Reston, VA: AIAA, 2016.

    [18] HIGASHI J, ISHIYAMA C, NAKAGAMI S, et al. Experimental study of disk-shaped rotating detonation turbine engine: AIAA-2017-1286[R]. Reston, VA: AIAA, 2017.

    [19] 卓長飛, 武曉松, 封鋒. 底部排氣彈三維湍流燃燒的數(shù)值模擬[J].固體火箭技術(shù), 2013, 36(6): 720-726.

    ZHUO C F, WU X S, FENG F. Numerical simulation of three-dimensional turbulent combustion of the base bleed projectile[J]. Journal of Solid Rocket Technology, 2013, 36(6): 720-726 (in Chinese).

    [20] 卓長飛, 武曉松, 封峰. 超聲速流動中底部排氣減阻的數(shù)值研究[J].兵工學報, 2014, 35(1): 18-26.

    ZHUO C F, WU X S, FENG F. Numerical research on drag reduction of base bleed in supersonic flow[J]. Acta Armamentarii, 2014, 35(1): 18-26 (in Chinese).

    [21] 馬虎, 武曉松, 王棟, 等. 旋轉(zhuǎn)爆震發(fā)動機數(shù)值研究[J]. 推進技術(shù), 2012, 33(5): 820-825.

    MA H, WU X S, WANG D, et al. Numerical investigation for rotating detonation engine[J]. Journal of Propulsion Technology, 2012, 33(5): 820-825 (in Chinese).

    [22] 夏鎮(zhèn)娟, 武曉松, 馬虎, 等. 圓盤結(jié)構(gòu)下旋轉(zhuǎn)爆震波的二維數(shù)值研究[J]. 推進技術(shù), 2017, 38(6): 1409-1418.

    XIA Z J, WU X S, MA H, et al. Two-dimensional numerical simulation of rotating detonation wave in plane-radial structure[J]. Journal of Propulsion Technology, 2017, 38(6): 1409-1418 (in Chinese).

    猜你喜歡
    爆震前導激波
    雷克薩斯車系爆震控制基理介紹
    一種基于聚類分析的二維激波模式識別算法
    航空學報(2020年8期)2020-09-10 03:25:34
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    基于“三思而行”的數(shù)學章前導學課設(shè)計——以《數(shù)的開方》(導學課)為例
    肺爆震傷治療的研究進展
    斜激波入射V形鈍前緣溢流口激波干擾研究
    一種S模式ADS-B前導脈沖檢測方法
    適于可壓縮多尺度流動的緊致型激波捕捉格式
    長距離爆震式點火槍設(shè)計
    焊接(2015年9期)2015-07-18 11:03:52
    3.0T磁敏感加權(quán)成像在兔顱腦爆震傷早期出血灶檢測及預(yù)后評估中的應(yīng)用
    婷婷色麻豆天堂久久| 麻豆国产97在线/欧美| 欧美日韩精品成人综合77777| 3wmmmm亚洲av在线观看| 九九久久精品国产亚洲av麻豆| 精品一区二区免费观看| 久久这里只有精品中国| 九九久久精品国产亚洲av麻豆| 一级片'在线观看视频| 我要看日韩黄色一级片| 日日摸夜夜添夜夜爱| 日韩欧美精品免费久久| 午夜爱爱视频在线播放| 婷婷色av中文字幕| 亚洲图色成人| 少妇人妻精品综合一区二区| 久久久久久久大尺度免费视频| 极品少妇高潮喷水抽搐| 成人美女网站在线观看视频| 国产精品不卡视频一区二区| 人妻制服诱惑在线中文字幕| 大陆偷拍与自拍| 男的添女的下面高潮视频| 国产精品无大码| 久久久成人免费电影| 亚洲内射少妇av| 观看美女的网站| 亚洲综合精品二区| 久久久久精品性色| 久热久热在线精品观看| 丝袜喷水一区| 欧美日本视频| 免费av不卡在线播放| av网站免费在线观看视频 | 亚洲国产精品专区欧美| 一边亲一边摸免费视频| 一个人免费在线观看电影| 亚洲国产最新在线播放| 国产黄色视频一区二区在线观看| 久久久久久久大尺度免费视频| 亚洲在线自拍视频| 中文资源天堂在线| 三级男女做爰猛烈吃奶摸视频| 精品人妻偷拍中文字幕| 国产 一区精品| 美女黄网站色视频| 久久久精品欧美日韩精品| 大片免费播放器 马上看| 成人性生交大片免费视频hd| 亚洲av不卡在线观看| 夜夜爽夜夜爽视频| 久久精品国产亚洲av涩爱| av.在线天堂| 91aial.com中文字幕在线观看| 97超视频在线观看视频| 嫩草影院新地址| 婷婷色综合www| 色视频www国产| 中文在线观看免费www的网站| 国产亚洲精品av在线| 国产精品1区2区在线观看.| 亚洲精品视频女| 国产伦精品一区二区三区视频9| 丰满乱子伦码专区| 好男人视频免费观看在线| 久久久久精品性色| 18禁动态无遮挡网站| 久久久久久久久久久丰满| 精品一区二区三卡| 看免费成人av毛片| 边亲边吃奶的免费视频| 免费播放大片免费观看视频在线观看| 国产精品三级大全| 国产精品女同一区二区软件| av卡一久久| 一级毛片我不卡| 在现免费观看毛片| 插阴视频在线观看视频| 一二三四中文在线观看免费高清| 亚洲内射少妇av| 精品午夜福利在线看| 国产 一区精品| 国产不卡一卡二| 亚洲熟女精品中文字幕| 亚洲欧美精品专区久久| 亚洲av免费在线观看| 精品久久久久久久久av| 成人特级av手机在线观看| 国产亚洲午夜精品一区二区久久 | 国产极品天堂在线| 成人国产麻豆网| 91av网一区二区| 在线观看av片永久免费下载| 国产精品.久久久| 80岁老熟妇乱子伦牲交| 色尼玛亚洲综合影院| 亚洲精品中文字幕在线视频 | 亚洲高清免费不卡视频| 精品久久久久久久末码| 亚洲内射少妇av| 日本与韩国留学比较| 日韩 亚洲 欧美在线| 在线免费观看的www视频| 亚洲熟妇中文字幕五十中出| 亚洲av福利一区| 久久精品综合一区二区三区| 免费高清在线观看视频在线观看| 永久网站在线| 我的老师免费观看完整版| 天堂中文最新版在线下载 | 亚洲av免费在线观看| 亚洲国产欧美人成| 如何舔出高潮| 国产高清不卡午夜福利| 亚洲av男天堂| 亚洲成人久久爱视频| av线在线观看网站| 国产在线一区二区三区精| 免费无遮挡裸体视频| 国内少妇人妻偷人精品xxx网站| 高清日韩中文字幕在线| 久久久色成人| 日日摸夜夜添夜夜添av毛片| 91av网一区二区| 久久99热这里只频精品6学生| 99热6这里只有精品| 91精品一卡2卡3卡4卡| 亚洲精品日韩在线中文字幕| 一级毛片久久久久久久久女| 国产在线一区二区三区精| 最近中文字幕2019免费版| 亚洲av免费在线观看| 最近中文字幕高清免费大全6| 青春草国产在线视频| 欧美一级a爱片免费观看看| 精品一区二区三区人妻视频| 国产成人91sexporn| 观看美女的网站| 91精品国产九色| 国产欧美日韩精品一区二区| 免费人成在线观看视频色| av线在线观看网站| 国产成人精品久久久久久| 99久久精品热视频| 丰满人妻一区二区三区视频av| 久久精品久久久久久噜噜老黄| 99久久精品热视频| 国产黄色小视频在线观看| 午夜福利视频精品| 久久精品夜色国产| 成人综合一区亚洲| 人人妻人人澡人人爽人人夜夜 | 国产女主播在线喷水免费视频网站 | 亚洲va在线va天堂va国产| 亚洲四区av| 91精品国产九色| 啦啦啦韩国在线观看视频| 欧美97在线视频| 女人久久www免费人成看片| 日日摸夜夜添夜夜添av毛片| 精品久久久久久成人av| 国产精品蜜桃在线观看| 蜜桃久久精品国产亚洲av| 九九爱精品视频在线观看| 成人欧美大片| 色综合色国产| 婷婷色麻豆天堂久久| 久久精品久久久久久噜噜老黄| 九色成人免费人妻av| 大片免费播放器 马上看| 99久国产av精品国产电影| 人妻少妇偷人精品九色| 极品教师在线视频| 大又大粗又爽又黄少妇毛片口| 三级经典国产精品| 婷婷色麻豆天堂久久| 欧美精品国产亚洲| 丝袜美腿在线中文| 精品一区二区三区视频在线| 毛片女人毛片| 午夜爱爱视频在线播放| 亚洲一区高清亚洲精品| 80岁老熟妇乱子伦牲交| 高清午夜精品一区二区三区| 高清午夜精品一区二区三区| 亚洲精品一二三| 久久久久久久久久人人人人人人| 亚洲电影在线观看av| 国产成人精品久久久久久| 亚洲伊人久久精品综合| 国产男人的电影天堂91| 男人爽女人下面视频在线观看| 欧美 日韩 精品 国产| 国产日韩欧美在线精品| 色5月婷婷丁香| 欧美精品一区二区大全| 哪个播放器可以免费观看大片| 免费观看av网站的网址| 三级男女做爰猛烈吃奶摸视频| 国产精品国产三级专区第一集| 婷婷色麻豆天堂久久| 美女大奶头视频| 男女边吃奶边做爰视频| 嫩草影院精品99| 亚洲av电影在线观看一区二区三区 | 人妻少妇偷人精品九色| 亚洲美女搞黄在线观看| 精品一区二区三卡| 精品熟女少妇av免费看| 精华霜和精华液先用哪个| 午夜福利成人在线免费观看| 欧美变态另类bdsm刘玥| 中文欧美无线码| 午夜福利在线观看吧| www.色视频.com| 老司机影院成人| 人妻制服诱惑在线中文字幕| 97人妻精品一区二区三区麻豆| 视频中文字幕在线观看| 日日啪夜夜爽| 欧美bdsm另类| 少妇熟女欧美另类| 亚洲电影在线观看av| 五月玫瑰六月丁香| 午夜精品国产一区二区电影 | 99re6热这里在线精品视频| 七月丁香在线播放| 成人二区视频| 一区二区三区高清视频在线| 亚洲国产精品成人综合色| 精品久久久久久久人妻蜜臀av| 最新中文字幕久久久久| 人妻系列 视频| 国产有黄有色有爽视频| 天堂√8在线中文| 18禁动态无遮挡网站| 国产精品久久久久久久电影| 精品一区在线观看国产| 青春草视频在线免费观看| 五月伊人婷婷丁香| 国产精品久久久久久精品电影小说 | 麻豆国产97在线/欧美| 亚洲国产最新在线播放| 人妻系列 视频| 国产在线男女| 亚洲欧美日韩东京热| 91午夜精品亚洲一区二区三区| 国产成人精品婷婷| 午夜精品在线福利| av线在线观看网站| 精品熟女少妇av免费看| 99热网站在线观看| 午夜福利视频1000在线观看| 国产精品伦人一区二区| 国产精品不卡视频一区二区| av播播在线观看一区| 一级二级三级毛片免费看| 三级国产精品片| 国产午夜精品论理片| 在线观看一区二区三区| 国模一区二区三区四区视频| 国产 一区 欧美 日韩| 亚洲第一区二区三区不卡| 久久99精品国语久久久| 免费大片18禁| 国产亚洲一区二区精品| 天天一区二区日本电影三级| 国产午夜精品久久久久久一区二区三区| 国产综合精华液| 91精品国产九色| 久久久久久久久久人人人人人人| 22中文网久久字幕| 狂野欧美激情性xxxx在线观看| 国产精品女同一区二区软件| 又爽又黄无遮挡网站| 色5月婷婷丁香| 少妇人妻一区二区三区视频| 天堂√8在线中文| 免费黄色在线免费观看| 日韩欧美一区视频在线观看 | 男女视频在线观看网站免费| 国产伦精品一区二区三区四那| 国产成人91sexporn| 欧美丝袜亚洲另类| 中文天堂在线官网| 久久久久网色| 国产有黄有色有爽视频| 欧美zozozo另类| 亚洲人成网站在线观看播放| 我的老师免费观看完整版| 性插视频无遮挡在线免费观看| 91狼人影院| 一边亲一边摸免费视频| 中文在线观看免费www的网站| 亚洲熟女精品中文字幕| 日日干狠狠操夜夜爽| 久久久久久久久中文| 美女大奶头视频| 国产麻豆成人av免费视频| 国产精品一区二区三区四区免费观看| 观看美女的网站| 波多野结衣巨乳人妻| 日韩欧美三级三区| 久久精品国产鲁丝片午夜精品| 国产精品国产三级专区第一集| 麻豆成人午夜福利视频| 欧美日韩综合久久久久久| 亚洲av电影在线观看一区二区三区 | 亚洲四区av| 免费无遮挡裸体视频| 免费看美女性在线毛片视频| 91狼人影院| 亚洲欧洲日产国产| 联通29元200g的流量卡| 亚洲乱码一区二区免费版| 久99久视频精品免费| 一个人看的www免费观看视频| 精品久久久久久久久av| 高清在线视频一区二区三区| 久久久久国产网址| 日本免费a在线| 波多野结衣巨乳人妻| 国产国拍精品亚洲av在线观看| 乱码一卡2卡4卡精品| 看非洲黑人一级黄片| 亚洲精品第二区| 国产爱豆传媒在线观看| 视频中文字幕在线观看| 国产在线男女| 丰满乱子伦码专区| 欧美高清成人免费视频www| 国产国拍精品亚洲av在线观看| 亚洲不卡免费看| 街头女战士在线观看网站| 久久精品夜色国产| 婷婷色麻豆天堂久久| 亚洲aⅴ乱码一区二区在线播放| 色5月婷婷丁香| 亚洲最大成人av| 久久精品熟女亚洲av麻豆精品 | 欧美xxxx性猛交bbbb| 夜夜看夜夜爽夜夜摸| 国产国拍精品亚洲av在线观看| 免费黄频网站在线观看国产| 亚洲av男天堂| 国产大屁股一区二区在线视频| 久久这里只有精品中国| 亚洲图色成人| 国产精品一及| 最近最新中文字幕大全电影3| 人体艺术视频欧美日本| 午夜激情久久久久久久| 日韩亚洲欧美综合| 黄色欧美视频在线观看| 晚上一个人看的免费电影| 黄色一级大片看看| 亚洲av成人av| 人妻系列 视频| 1000部很黄的大片| 成人毛片60女人毛片免费| 久久久精品欧美日韩精品| 亚洲国产日韩欧美精品在线观看| 18+在线观看网站| 国产欧美另类精品又又久久亚洲欧美| 亚洲精品456在线播放app| h日本视频在线播放| 日本午夜av视频| 乱系列少妇在线播放| 青春草视频在线免费观看| av国产免费在线观看| 一级毛片我不卡| 亚洲国产日韩欧美精品在线观看| 日韩精品有码人妻一区| 免费观看的影片在线观看| 亚洲精华国产精华液的使用体验| av在线观看视频网站免费| 久久99热这里只频精品6学生| 女人被狂操c到高潮| 亚洲无线观看免费| av免费观看日本| 成年人午夜在线观看视频 | 看十八女毛片水多多多| www.色视频.com| 99热6这里只有精品| 欧美三级亚洲精品| 日本av手机在线免费观看| av黄色大香蕉| 男女下面进入的视频免费午夜| 麻豆成人午夜福利视频| 伦精品一区二区三区| 3wmmmm亚洲av在线观看| 亚洲综合精品二区| 特级一级黄色大片| 少妇的逼好多水| 特大巨黑吊av在线直播| 国产av码专区亚洲av| 欧美日韩视频高清一区二区三区二| 日韩视频在线欧美| 精品久久久久久久久亚洲| 日韩亚洲欧美综合| 少妇高潮的动态图| 国产精品女同一区二区软件| 欧美 日韩 精品 国产| 三级国产精品欧美在线观看| 久久久久久国产a免费观看| 免费黄网站久久成人精品| 亚洲激情五月婷婷啪啪| 嫩草影院新地址| 日韩一区二区视频免费看| 亚洲精品乱码久久久久久按摩| 少妇裸体淫交视频免费看高清| 一级黄片播放器| 久久久久久久国产电影| 人体艺术视频欧美日本| 亚洲熟女精品中文字幕| 亚洲精品影视一区二区三区av| 国产中年淑女户外野战色| 韩国高清视频一区二区三区| 亚洲性久久影院| 18禁在线播放成人免费| 九九爱精品视频在线观看| 99视频精品全部免费 在线| 99久久精品热视频| 欧美成人精品欧美一级黄| 久久久亚洲精品成人影院| 亚洲高清免费不卡视频| 精品国内亚洲2022精品成人| 99久久中文字幕三级久久日本| 久久精品夜色国产| 身体一侧抽搐| 天天躁日日操中文字幕| 国产成人精品一,二区| 亚洲一级一片aⅴ在线观看| 丝瓜视频免费看黄片| 亚洲自拍偷在线| 最近的中文字幕免费完整| 国产午夜精品论理片| 久久久精品欧美日韩精品| 国产 一区 欧美 日韩| 日韩中字成人| 美女高潮的动态| 男的添女的下面高潮视频| 国语对白做爰xxxⅹ性视频网站| 欧美高清成人免费视频www| av免费在线看不卡| 成人鲁丝片一二三区免费| 国产男女超爽视频在线观看| 精品久久国产蜜桃| 中文在线观看免费www的网站| 日韩av在线免费看完整版不卡| 精品人妻熟女av久视频| 日韩欧美精品v在线| 国产精品一二三区在线看| 午夜福利在线在线| 国产一级毛片在线| 日韩大片免费观看网站| av又黄又爽大尺度在线免费看| 午夜爱爱视频在线播放| 久久精品人妻少妇| 一级av片app| 成人国产麻豆网| 亚洲在线自拍视频| 日韩av免费高清视频| 久久久久免费精品人妻一区二区| 青青草视频在线视频观看| 人妻制服诱惑在线中文字幕| 大香蕉97超碰在线| 亚洲精品久久久久久婷婷小说| 欧美成人精品欧美一级黄| 日韩制服骚丝袜av| 久久久精品免费免费高清| 99九九线精品视频在线观看视频| 国产白丝娇喘喷水9色精品| 精品久久久精品久久久| 一区二区三区高清视频在线| 能在线免费观看的黄片| 欧美日韩综合久久久久久| 国产精品久久久久久av不卡| 91久久精品电影网| 男人狂女人下面高潮的视频| 99久久精品热视频| 成人欧美大片| 日本免费a在线| 人体艺术视频欧美日本| 狂野欧美白嫩少妇大欣赏| 免费观看a级毛片全部| 亚洲成人一二三区av| 国产黄频视频在线观看| 99久久人妻综合| 国产成年人精品一区二区| 国产乱人偷精品视频| 亚洲内射少妇av| 人妻一区二区av| 一区二区三区免费毛片| 看十八女毛片水多多多| 国语对白做爰xxxⅹ性视频网站| 亚洲激情五月婷婷啪啪| 欧美极品一区二区三区四区| 99久久中文字幕三级久久日本| 国产精品麻豆人妻色哟哟久久 | 在线免费十八禁| 日韩av不卡免费在线播放| av女优亚洲男人天堂| 2021天堂中文幕一二区在线观| 久久精品国产亚洲av天美| 日韩制服骚丝袜av| 久久久久久久久中文| av在线天堂中文字幕| 国产黄片美女视频| 亚洲欧美中文字幕日韩二区| 91午夜精品亚洲一区二区三区| 久久久久久久国产电影| 成人国产麻豆网| 秋霞伦理黄片| 女人十人毛片免费观看3o分钟| 免费观看av网站的网址| 日韩av在线免费看完整版不卡| 晚上一个人看的免费电影| xxx大片免费视频| 亚洲无线观看免费| 春色校园在线视频观看| 亚洲av中文av极速乱| 中文字幕免费在线视频6| 亚洲国产av新网站| 国产欧美日韩精品一区二区| 青春草亚洲视频在线观看| 熟妇人妻久久中文字幕3abv| 午夜老司机福利剧场| 亚洲精品久久久久久婷婷小说| 久久人人爽人人片av| 美女主播在线视频| 国产日韩欧美在线精品| 久久精品久久精品一区二区三区| 国产男女超爽视频在线观看| 好男人在线观看高清免费视频| 亚洲精品日韩av片在线观看| 少妇人妻精品综合一区二区| 欧美三级亚洲精品| 久热久热在线精品观看| 欧美bdsm另类| 80岁老熟妇乱子伦牲交| 毛片一级片免费看久久久久| 欧美三级亚洲精品| 久久草成人影院| 国产精品久久久久久精品电影| 亚洲欧美成人精品一区二区| 久久97久久精品| 国产成人精品久久久久久| 久99久视频精品免费| 街头女战士在线观看网站| 51国产日韩欧美| 99久久精品热视频| 久久亚洲国产成人精品v| 日韩欧美精品免费久久| 日产精品乱码卡一卡2卡三| 国产乱人偷精品视频| 夫妻午夜视频| 亚洲成色77777| 日韩一区二区三区影片| 日韩欧美精品v在线| 国产单亲对白刺激| 日韩av免费高清视频| 欧美+日韩+精品| 亚洲四区av| 欧美日韩精品成人综合77777| 熟妇人妻不卡中文字幕| 一级毛片电影观看| 最近中文字幕2019免费版| 国产精品国产三级国产专区5o| 婷婷六月久久综合丁香| 欧美高清成人免费视频www| 最后的刺客免费高清国语| 中文字幕av成人在线电影| 国产久久久一区二区三区| 色哟哟·www| 亚洲av二区三区四区| 久久精品国产亚洲网站| 亚洲欧美成人综合另类久久久| 青春草视频在线免费观看| 国产av不卡久久| 免费av观看视频| 国产高清有码在线观看视频| 伊人久久精品亚洲午夜| 日韩av在线大香蕉| 国产免费福利视频在线观看| 能在线免费观看的黄片| 亚洲成人精品中文字幕电影| 女人被狂操c到高潮| 精品久久久久久电影网| 黄色欧美视频在线观看| 搡老乐熟女国产| 不卡视频在线观看欧美| 一区二区三区免费毛片| 人妻夜夜爽99麻豆av| 国产中年淑女户外野战色| 只有这里有精品99| 亚洲国产成人一精品久久久| www.av在线官网国产| 美女cb高潮喷水在线观看| 免费大片黄手机在线观看| 国产av国产精品国产| 十八禁网站网址无遮挡 | 免费看日本二区| 一级毛片 在线播放| 美女被艹到高潮喷水动态| 黑人高潮一二区| 日韩成人伦理影院| 久久97久久精品| 全区人妻精品视频| 精品人妻熟女av久视频| 久久久午夜欧美精品| 亚洲av不卡在线观看| 最新中文字幕久久久久| 久久久a久久爽久久v久久| 亚洲一级一片aⅴ在线观看|