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

    S形進氣道錘激波的動態(tài)特性

    2018-03-15 09:50:30劉庭申馮立好王晉軍
    航空學報 2018年2期
    關(guān)鍵詞:進氣道馬赫數(shù)激波

    劉庭申,馮立好,王晉軍

    北京航空航天大學 航空科學與工程學院,北京 100083

    進氣道的錘激波載荷由發(fā)動機強喘振引起。發(fā)動機喘振主要由壓氣機葉片的部分或全部失速所產(chǎn)生,可導致進入發(fā)動機的流量突然減少,同時伴隨著發(fā)動機進口壓力的急劇升高。這種壓力的突變會產(chǎn)生向上游傳播的壓力波,一般稱為錘激波。錘激波的峰值壓力遠遠大于定常狀態(tài)的壓力值,通常能達到定常值的兩倍以上。因此盡管錘激波現(xiàn)象并不常見,錘激波載荷卻是飛機進氣道結(jié)構(gòu)強度設(shè)計的決定性因素之一[1]。如果進氣道的結(jié)構(gòu)強度無法承受錘激波載荷的沖擊,發(fā)動機出現(xiàn)喘振時可能會導致進氣道的破壞,進而威脅飛機的飛行安全。如果進氣道結(jié)構(gòu)強度設(shè)計過于保守,則會使得進氣道重量太大影響到飛機的性能。所以,對進氣道錘激波現(xiàn)象進行研究,獲得較為準確的進氣道錘激波載荷具有非常重要的價值。

    國外已經(jīng)進行了一些關(guān)于進氣道錘激波現(xiàn)象的研究。一些研究者采用試驗方法探究了飛機進氣道的錘激波現(xiàn)象[1-6],通過各種方法誘使發(fā)動機喘振,進而進行數(shù)據(jù)的測量。Marshall[2]依據(jù)大量的試驗測量數(shù)據(jù)總結(jié)了發(fā)動機喘振時壓力變化過程以及喘振超壓的估算方法,可以在只做較少次數(shù)的發(fā)動機喘振試驗的情況下估算多種不同狀態(tài)的發(fā)動機喘振壓升情況,為進氣道的錘激波研究提供了基礎(chǔ)。Nugent和Holzman[5]對F-111A飛機進行飛行試驗,獲得了多種飛行狀態(tài)下發(fā)生喘振時發(fā)動機和進氣道壓力隨時間變化情況以及錘激波的傳播速度等數(shù)據(jù)。在數(shù)值計算方面,早期的研究主要針對二維管道開展[7-8]。Hindash等[7]通過對二維管道的錘激波模擬分析了進氣道的受載情況并與試驗數(shù)據(jù)對比驗證了計算的準確性,同時文章進一步分析了進氣道出口不均勻壓升產(chǎn)生的錘激波的傳播特點以及不同進氣道幾何形狀對錘激波傳播過程的影響。

    此后一些學者陸續(xù)進行了三維數(shù)值模擬研究,并且應用到了多種戰(zhàn)斗機進氣道錘激波現(xiàn)象的研究中[9-12]。Goble等[9]進行了F-22飛機進氣道的錘激波三維數(shù)值模擬研究,他們采用在進氣道出口設(shè)定隨時間變化的壓力作為邊界條件來產(chǎn)生錘激波,計算了多種飛行條件下的錘激波情況,總結(jié)了進氣道所受的壓力分布規(guī)律,并發(fā)現(xiàn)了進氣道的支路和抽吸孔排氣可以降低進氣道所受的錘激波峰值壓力。類似的錘激波模擬方法應用到了瑞典JAS-39飛機的進氣道錘激波研究中[10],作者討論了動態(tài)計算時間步長的選取,以及加快計算速度的方法,探討了不同計算模型在計算速度和準確性方面的優(yōu)劣。根據(jù)計算結(jié)果獲得了進氣道壁面的錘激波載荷數(shù)據(jù),并為飛機的進氣道強度設(shè)計提供了依據(jù)。數(shù)值方法的應用可以較為準確地獲得進氣道所受錘激波壓力的分布情況,為進氣道結(jié)構(gòu)強度設(shè)計提供更豐富的信息。Gridley等[13]進一步提出了將不同飛行狀態(tài)的錘激波載荷和各狀態(tài)出現(xiàn)的概率相結(jié)合的方法,允許在小概率情況下出現(xiàn)錘激波載荷超出進氣道結(jié)構(gòu)設(shè)計強度,從而進一步降低進氣道的結(jié)構(gòu)重量。

    國內(nèi)關(guān)于進氣道錘激波載荷的相關(guān)文獻較少。主要由朱宇和沈天榮[14]歸納總結(jié)了國外錘激波載荷的研究方法和成果,而關(guān)于試驗和數(shù)值計算的研究成果很少。目前,國內(nèi)主要采用工程估算的方法評估發(fā)動機喘振超壓比,以此作為進氣道錘激波載荷。與此同時,如果能輔以數(shù)值模擬計算錘激波沿管道運動過程中載荷變化情況,可以更加有效地為進氣道的結(jié)構(gòu)設(shè)計提供幫助。

    雖然國外已經(jīng)對錘激波現(xiàn)象進行了一些研究,但是大部分只提供了進氣道壁面壓力分布,較少關(guān)注錘激波過程中進氣道內(nèi)部的流動結(jié)構(gòu)變化,尤其是內(nèi)部復雜的三維流動特性。本文對進氣道錘激波進行研究,探討適用于進氣道錘激波計算的三維非定常數(shù)值模擬方法,獲得錘激波傳播過程中進氣道壁面壓力變化情況,并且對進氣道內(nèi)部的三維流動結(jié)構(gòu)和動態(tài)演化過程進行分析總結(jié)。

    1 模型和計算方法

    本文采用的計算模型為M2129型S形進氣道。這是一種圓形截面進氣道。它由AGARD(Advisory Group for Aerospace Research and Development)設(shè)計,并進行了風洞試驗[15]。試驗在來流馬赫數(shù)為0.21下進行,包含了兩種不同的流量條件,兩種流量下喉道的馬赫數(shù)分別為0.412和0.794。試驗所測得的進氣道壁面上的壓力分布數(shù)據(jù)可以獲得,經(jīng)常被用于數(shù)值模擬結(jié)果的驗證對比[16-18]。

    M2129型S形進氣道的剖面幾何形狀見圖1。該進氣道模型包括外部壁面、筆直的入口段和出口段,以及中間的S段。為了接近常見飛行器的進氣道尺寸,本文中的進氣道模型尺寸相比于原試驗中的模型擴大了4.65倍。其入口段的直徑為0.6 m,出口段直徑約為0.7 m,入口段長度為0.6 m,中間S段的中軸線長度為2.184 m。由于模型是對稱的,僅選取了一半的結(jié)構(gòu)進行計算。為了記錄進氣道壁面壓力隨時間的動態(tài)變化,選取了上下壁面各11個點進行記錄,即P1~P11和S1~S11。22個點的編號和位置也在圖1中標出,第1個點同進氣道入口的距離和后面每相鄰兩個點沿x方向的距離均為0.3 m。

    圖1 M2129型S形進氣道幾何模型Fig.1 M2129 S-shaped inlet geometry model

    圖2 進氣道壁面網(wǎng)格劃分Fig.2 Grid distribution on inlet wall

    整個計算域使用結(jié)構(gòu)網(wǎng)格進行劃分,進氣道壁面上網(wǎng)格如圖2所示。為了精確模擬進氣道內(nèi)部的流動,進氣道內(nèi)部流場采用了較密的網(wǎng)格。近壁面的網(wǎng)格也進行了加密來模擬壁面流動。網(wǎng)格數(shù)量一共約200萬。

    本文采用非定常雷諾平均Navier-Stokes方程進行數(shù)值計算,采用耦合求解器求解。外流場設(shè)定為壓力遠場邊界條件,給定來流的馬赫數(shù)和靜壓,在進氣道出口給定反壓。首先采用定常計算的方法獲得初值,再在進氣道出口給定壓力隨時間變化的函數(shù)來模擬發(fā)動機喘振的壓力變化,進而進行錘激波的動態(tài)計算。這種錘激波的數(shù)值模擬方法將在后文中詳細介紹。動態(tài)計算的時間步長為2.5×10-5s。計算所選用的湍流模型為k-ε模型。這種湍流模型已經(jīng)由Nichols[16]在M2129型進氣道模型上進行了靜態(tài)的計算,并與試驗結(jié)果對比驗證了正確性。

    本文為了驗證計算的準確性,將風洞試驗數(shù)據(jù)和計算所獲得的結(jié)果進行了對比。試驗數(shù)據(jù)來自文獻[15]。因為試驗的低流量狀態(tài)和后面用于錘激波計算的流量狀況比較接近,選取該狀態(tài)的試驗數(shù)據(jù)進行對比驗證。

    試驗條件是來流馬赫數(shù)為0.21,進氣道喉道處的馬赫數(shù)為0.412,雷諾數(shù)為1.158×106。圖3為試驗和計算獲得的進氣道上下壁面的壓力分布結(jié)果,圖中橫坐標為各點x坐標與該方向總長度L的比值,縱坐標為靜壓p與來流總壓pt的比值。上壁面壓力的計算結(jié)果與試驗數(shù)據(jù)非常吻合,下壁面計算結(jié)果數(shù)值略高于試驗數(shù)據(jù),但是誤差在3%以內(nèi)。表明本文所用方法可以對M2129型進氣道內(nèi)部的流動情況進行準確的模擬。

    圖3 進氣道上下壁面壓力計算與試驗結(jié)果對比Fig.3 Comparison between numerical and experimental results of inlet upper and lower walls

    為進一步驗證非定常計算的準確性,對激波管問題進行了數(shù)值模擬,并與理論解進行了比較,計算所用激波管長度為1 m,均分為等長的高壓段和低壓段,之間由薄片隔開。高壓段壓力為1 Pa,密度為1 kg/m3。低壓段壓力為0.1 Pa,密度為0.125 kg/m3,計算從兩段氣體之間的薄片突然撤離開始。圖4為時間t=0.2 s時激波管內(nèi)的壓力p、密度ρ、速度v分布情況的計算結(jié)果和理論值的對比。結(jié)果表明本文所用數(shù)值方法可以準確地模擬激波傳播的非定常過程。

    圖4 激波管內(nèi)壓力、密度和速度的數(shù)值與理論結(jié)果對比Fig.4 Comparison of pressure,density and velocity between numerical and theoretical results of shock tube

    2 錘激波波形設(shè)置

    通常所采用的計算方法是直接給定喘振時的壓力變化波形作為邊界條件[7-11],進一步進行錘激波的計算,本文采用同樣的方法。壓力波形的確定依賴于發(fā)動機喘振的實際測量數(shù)據(jù),不同發(fā)動機發(fā)生喘振時可能會有不同的壓力變化波形。波形的差異會影響錘激波的計算結(jié)果,因此應當盡量準確地模擬實際的壓力變化情況。

    通過對喘振時發(fā)動機壓力波形的試驗測量發(fā)現(xiàn)不同的發(fā)動機型號以及不同的喘振狀況產(chǎn)生的壓力變化波形會有所差異。但是其隨時間變化規(guī)律通常為迅速達到峰值之后緩慢地下降[7]。因此本文選用圖5所示三角形波形進行錘激波的數(shù)值模擬,其超壓比為2.1,即峰值壓力為定常值的2.1倍。出口壓力一開始保持恒定,在t=0.01 s時瞬間達到峰值。之后壓力呈線性下降,經(jīng)過0.1 s,即在t=0.11 s時壓力恢復到初始值。在此之后,壓力保持不變,直到計算結(jié)束。這種波形的壓力曲線如圖5所示,圖中縱坐標為各時刻壓力同定常時壓力pss的比值。

    圖5 錘激波壓力波形Fig.5 Pressure waveform of hammershock

    3 壓力變化與流動結(jié)構(gòu)

    本節(jié)進行錘激波計算的狀態(tài)為來流馬赫數(shù)Ma∞=0.65,來流靜壓p∞=101 325 Pa,靜溫T∞=288.15 K,進氣道定常時的換算流量Wc=60 kg/s(即進氣道出口馬赫數(shù)Mae=0.44,入口馬赫數(shù)Ma0=0.63)。進氣道上下壁面各點靜壓隨時間變化記錄在圖6中。各點所對應的位置已經(jīng)在圖1中顯示,圖中將各點的靜壓值除以來流靜壓p∞從而進行無量綱處理。圖6(a)為上壁面各點的結(jié)果,圖6(b)為下壁面各點的結(jié)果。從圖中可以明顯地看出沿著進氣道入口到出口各點壓力曲線的變化規(guī)律。

    圖6 進氣道上下壁面壓力隨時間變化曲線Fig.6 Time history of pressure on inlet upper and lower walls

    圖6(a)中沿實線箭頭方向依編號次序逐個經(jīng)過P1~P11的曲線,說明越靠近進氣道出口的點壓力開始升高的時間越早,反映了錘激波沿著進氣道逐漸向上游傳播的過程。觀察圖6(a)中壓力曲線在達到峰值之后的下降段,從P1~P11各點的曲線也沿著虛線箭頭方向排布,表明越靠近進氣道出口的點,壓力達到峰值之后下降得越慢。由于該種變化規(guī)律,使得曲線的形態(tài)呈現(xiàn)出越靠近進氣道入口的點的曲線顯得更尖,而越接近進氣道出口的點的曲線更寬。關(guān)于下壁面的點S1~S11,圖中箭頭也標示出了相似的變化規(guī)律。根據(jù)最靠近進氣道入口的點P1和S1的壓力達到峰值的時間,可計算出錘激波在進氣道內(nèi)移動的平均速度為306 m/s,馬赫數(shù)約為0.9。而錘激波相對于上游來流的平均傳播速度為481 m/s,馬赫數(shù)達1.42。

    圖7統(tǒng)計了上下壁面各11個點的峰值壓力和定常時的壓力值??梢钥吹娇傮w上靠近入口位置的載荷略大于出口附近載荷,最大峰值壓力位于上壁面靠近入口的彎折處。上下壁面各點的峰值壓力均達到來流靜壓的2.3~2.65倍左右,說明錘激波能夠誘導產(chǎn)生比定常狀態(tài)更高的載荷。

    圖7 進氣道上下壁面峰值壓力和定常狀態(tài)壓力Fig.7 Peak pressure and steady pressure on the inlet upper and lower walls

    圖8顯示了錘激波經(jīng)過進氣道時壁面上的壓力同來流靜壓比值(p/p∞)的分布云圖;圖9為上下壁面沿程壓力分布曲線,從圖中可以看到錘激波沿進氣道向上游傳播時壓力變化的過程。圖8(a)為t=0.01 s時的壓力云圖,此時進氣道出口壓力剛產(chǎn)生突變,在進氣道壁面上還未產(chǎn)生明顯的影響。在t=0.015 s時,即圖8(b)所示,錘激波剛經(jīng)過進氣道的第1個彎折處,此時錘激波位置的下壁面受到的壓力略高于上壁面。在t=0.017 5 s和t=0.02 s兩個時刻,即圖8(c)和圖8(d)所示,錘激波正經(jīng)過進氣道第2個彎道,此時可以明顯看出錘激波處上壁面的壓力較高。由圖8所示的壓力云圖可以發(fā)現(xiàn),錘激波向上游傳播的過程中進氣道受到的最大壓力的區(qū)域位于進氣道入口附近彎折處的上壁面,這一區(qū)域所受峰值壓力達到來流靜壓的2.65倍左右。在錘激波傳出進氣道之后,進氣道出現(xiàn)了中部壓力較高,入口和出口附近壓力較低的特點,如圖8(e)所示。在這一中部高壓區(qū)逐漸后移并消散之后,進氣道便呈現(xiàn)出越靠近入口壓力越低的特點。

    圖8 進氣道壁面壓力云圖Fig.8 Pressure contours on inlet wall

    圖9 進氣道上下壁面沿程壓力分布Fig.9 Pressure distribution along inlet upper and lower walls

    對于圖9所示壓力沿程分布曲線。在t=0.01 s時剛剛產(chǎn)生錘激波,因而可以看到相應時刻曲線在出口處(x=3.4 m)壓力突增,而其他部位的壓力分布與穩(wěn)態(tài)時相同。在t=0.02 s時,錘激波傳至進氣道入口附近??梢钥吹?,在激波面下游,上下壁面壓力沿進氣道軸向以相近的梯度逐漸下降,與所設(shè)定的錘激波波形的形態(tài)類似。錘激波傳出進氣道之后,壁面壓力整體上開始逐漸下降。在t=0.03 s時,進氣道壁面壓力呈現(xiàn)入口處低、出口處高的特點,并且上下壁面沿軸向壓力以近似的梯度增高。此時入口處壓力已經(jīng)降至接近當?shù)胤€(wěn)態(tài)時壓力。而在之后時刻,入口附近的壓力會進一步下降直至低于當?shù)胤€(wěn)態(tài)時的壓力,典型的如圖中展示的t=0.08 s時刻的曲線。在該時刻,進氣道前段所受壓力低于其穩(wěn)態(tài)時壓力,而進氣道后段所受的壓力仍高于穩(wěn)態(tài)壓力。值得注意的是此時進氣道在靠近入口彎道處的下壁面區(qū)域。由于進氣道壁面在此處突然開始彎曲,在穩(wěn)態(tài)時,該處便是一個明顯的低壓區(qū)域。但在t=0.08 s時刻,該區(qū)域的壓力進一步降低,同時其壓力梯度也更大。這會使得進氣道在該處受到較大的吸力載荷以及較大的剪切應力。

    不同時刻進氣道對稱面上壓力云圖以及流線圖顯示在圖10中,圖中壓力均除以來流靜壓以進行無量綱處理。圖10(a)~圖10(d)顯示了錘激波在進氣道內(nèi)的傳播過程,激波面的法向基本與進氣道中軸線方向一致。位于錘激波上游的壓力和流線未受到影響,依舊與錘激波發(fā)生前一樣,但是在錘激波下游,則可以觀察到回流的產(chǎn)生。值得注意的是,圖10(c)錘激波經(jīng)過進氣道靠近入口的彎折處時,由于進氣道的彎曲,導致錘激波下游回流的流線向上壁面傾斜,錘激波處上壁面出現(xiàn)明顯的壓力集中,使得此處成為進氣道內(nèi)部峰值壓力最高的部位。在錘激波傳播過程中進氣道下游產(chǎn)生回流,而上游氣流仍從進氣道入口流入,因而在進氣道內(nèi)部堆積大量流體,形成高壓區(qū)域。圖10(d)中錘激波即將從進氣道傳出,激波面后方出現(xiàn)的一段高壓區(qū)域,使得流動出現(xiàn)分離,一部分繼續(xù)向上游流動,另一部分則流向下游,并在進氣道中部誘導產(chǎn)生復雜的螺旋結(jié)構(gòu)。在錘激波離開進氣道后,聚集在進氣道中部的高壓區(qū)導致氣流同時向進氣道入口和出口流出。由于在進氣道出口所設(shè)定的波形持續(xù)時間較長,因而入口附近壓力下降更快,因此這一高壓區(qū)域逐漸向后移動。圖10(e)中高壓區(qū)位于進氣道中部,流動在中部分離分別向入口和出口流出,而圖10(f)由于高壓區(qū)的后移,流動分離區(qū)域也移動到了進氣道的后段。中部高壓區(qū)持續(xù)了較短時間,在圖10(f)之后不久便完全消散,此時距離錘激波的產(chǎn)生時間僅過了0.02 s,而此時進氣道出口的壓力還處于較高的狀態(tài),進而導致進氣道再次出現(xiàn)回流,回流從進氣道入口流出,在進氣道外部產(chǎn)生旋渦,如圖10(g)所示。此時呈現(xiàn)出越靠近進氣道入口,壓力下降越快的特點?;亓鞒掷m(xù)了較長的時間,并隨著進氣道出口壓力的下降而逐漸減弱。在進氣道出口壓力恢復到初始壓力時,進氣道內(nèi)流動尚未恢復正常,而是直到0.195 s以后,即距錘激波發(fā)生0.185 s之后,流動恢復到初始方向并逐漸穩(wěn)定,如圖10(h)所示。

    圖10 不同時刻進氣道對稱面流線圖和壓力云圖Fig.10 Streamlines and pressure contours on symmetry plane of inlet at different time

    圖11 不同時刻進氣道內(nèi)部三維流線圖 (左:側(cè)視;右:俯視)Fig.11 3D streamlines in inlet duct at different time (left:side view; right:top view)

    為了進一步分析進氣道內(nèi)部的三維流動特征,圖11給出了不同時刻進氣道內(nèi)部的三維流線圖,圖中每一時刻均由兩種不同視角流線圖組成。圖中流線顏色表示速度的大小,單位為m/s。t=0.01 s時的流線圖如圖11(a)所示,此時進氣道出口剛產(chǎn)生壓力突變,內(nèi)部流動還未受到影響。在t=0.015 s時(圖11(b)),錘激波傳播至進氣道中部,可以看到錘激波下游的回流速度遠小于錘激波上游氣流的流速。從俯視圖中可以看到此時進氣道內(nèi)部流動沒有出現(xiàn)明顯的側(cè)向運動,流線基本與進氣道中軸線方向一致。t=0.02 s時(圖11(c)),位于進氣道中心的旋渦使得流線旋轉(zhuǎn)并互相纏繞著向進氣道兩側(cè)移動。同時值得注意的是,在錘激波下游的回流中,處于進氣道中心旋渦處的氣流流速非常低(小于5 m/s),而旋渦周邊靠近壁面的氣流流速相對較高,尤其是下壁面附近,流速為錘激波下游回流中最高。t=0.025 s時(圖11(d)),位于進氣道中部的旋渦可以在圖10(e)中觀察到其在對稱面的二維投影。而在進氣道后段的上壁面附近還可以觀察到另外兩個對稱的旋渦。t=0.03 s時(圖11(e)),氣流從高壓區(qū)同時向進氣道入口和出口流出,在進氣道內(nèi)越遠離高壓區(qū)的地方流速越快。圖11(f)為t=0.08 s中部高壓區(qū)消失后出現(xiàn)的回流,此時回流的流速遠高于錘激波經(jīng)過進氣道時其下游回流的流速,流線基本與進氣道中軸線平行,沒有出現(xiàn)明顯的三維流動特性。

    4 進氣道換算流量的影響

    進一步選用3種不同的流量條件進行錘激波過程的數(shù)值模擬。3種計算條件下穩(wěn)態(tài)時的換算流量分別為30、45、60 kg/s(對應的進氣道出口馬赫數(shù)分別為Mae=0.21,0.31,0.44)。通過改變進氣道出口的反壓來達到預計的流量條件,而來流條件以及非定常計算時錘激波波形均與第3節(jié)所述相同。

    為了對比換算流量變化的影響,圖12詳細列出了上下壁面各點在3種不同流量條件下的峰值壓力??梢钥闯?,隨著換算流量(或進氣道出口馬赫數(shù))的減小,進氣道各處所受的峰值壓力逐漸增大。

    圖12 不同流量條件時上下壁面各點的峰值壓力Fig.12 Peak pressure of each point on upper and lower walls for different flow rate conditions

    圖13 不同流量條件時點P1處壓力隨時間變化曲線Fig.13 Time history of pressure at probe point P1 for different flow rate conditions

    圖13給出了點P1在3種流量條件下壓力隨時間的變化過程,其中點P1的位置已經(jīng)在圖1中標出??梢钥闯鲈诓煌瑩Q算流量時,除了壓力的數(shù)值呈現(xiàn)差異,壓力開始變化的時間也有所不同,說明錘激波的傳播速度不同。依據(jù)進氣道入口附近點達到峰值的時間,在Wc=30 kg/s(出口馬赫數(shù)Mae=0.21)時,錘激波傳播速度達391 m/s(馬赫數(shù)Ma=1.15),而在Wc=60 kg/s(出口馬赫數(shù)Mae=0.44)時,錘激波相對進氣道的平均傳播速度為306 m/s(Ma=0.9)。

    5 超壓比的影響

    錘激波波形的超壓比(Over Pressure Ratio, OPR)表示波形壓力變化峰值與初始壓力的比值。為分析超壓比的差異對進氣道所受壓力的影響,將錘激波波形變更為峰值壓力達3倍初始壓力,而壓力突變?yōu)榉逯档臅r間以及壓力恢復到初始值的時間均與圖5所示的原波形相同。除了波形的改變外,計算時其他條件均保持相同。

    兩種超壓比條件下上下壁面峰值壓力如圖14所示。表明增大波形的超壓比顯著增大了進氣道壁面所受的峰值壓力。超壓比由2.1增大到3使得設(shè)定的波形的峰值壓力增大了約1.43倍。而進氣道各點的峰值壓力的增大比例為1.34~1.42倍,略低于波形中峰值壓力的增大比例。

    圖15給出了進氣道上壁面點P1的壓力變化曲線,可以看出在超壓比為3的條件下,點P1達到峰值的時間更早,說明錘激波在進氣道內(nèi)傳播速度更快。在超壓比為3時,錘激波在進氣道內(nèi)傳播的平均速度為379 m/s(Ma=1.11)。

    圖14 不同超壓比時上下壁面各點的峰值壓力Fig.14 Peak pressure of each point on upper and lower walls for different Over Pressure Ratios(OPRs)

    圖15 不同超壓比時點P1處壓力隨時間變化曲線Fig.15 Time history of pressure at probe point P1 for different over pressure ratios

    6 中心線幾何形狀的影響

    進氣道的幾何形狀會對錘激波壓力的分布情況產(chǎn)生影響。對于S形進氣道幾何造型方面的研究較多,其中Lee和Boedicker[19]提出的關(guān)于中心線形狀的3種曲線方程得到了非常廣泛的應用。本文選取其中兩條曲線與原曲線進行對比,探討其錘激波壓力分布的差異。3條曲線如圖16所示,曲線1為原M2129型進氣道中心線,曲線2和3來自文獻[19],圖中橫縱軸分別表示中心線各點x坐標和z坐標與對應方向總長度的比值。對于前半段(Δx<0.5),相比于曲線1,曲線2初始時曲率較小,方向變化緩慢,之后曲率迅速增大,而曲線3則一開始曲率最大,方向變化最快,之后變得最為平緩。在后半段(Δx>0.5),曲線2和3變化規(guī)律相反。計算時除了中心線形狀外,其他參數(shù)的設(shè)定完全相同。

    圖16 進氣道中心線幾何形狀Fig.16 Geometry shape of inlet centerlines

    圖17 不同中心線時上下壁面各點的峰值壓力Fig.17 Peak pressure of each point on upper and lower walls with different centerlines

    圖17給出了3種不同中心線形狀的進氣道上下壁面各點在錘激波演化過程中的峰值壓力。對比可以看到峰值壓力的變化在進氣道前半段的上壁面最為明顯,最大峰值壓力的位置和大小均有所不同。曲線2開始最平緩之后曲率迅速變大,相應的進氣道靠近入口位置峰值壓力較小,最大峰值所處位置更靠后。曲線3開始曲率很大,之后平緩,相應的進氣道模型受到最大峰值壓力的地方更靠前。而進氣道后半段各點的峰值壓力受幾何形狀的影響較小。通過對比可以看出進氣道中心線形狀的變化對其錘激波壓力分布具有重要的影響。但是峰值壓力還受到其他因素的影響,如所處進氣道的位置、進氣道截面積變化等。因而各處峰值壓力的變化并不與曲率分布完全對應。

    7 結(jié) 論

    本文選取M2129型S形進氣道模型進行了錘激波動態(tài)過程的數(shù)值模擬。利用在進氣道出口設(shè)定壓力波形的方法模擬喘振時發(fā)動機的壓力變化。詳細分析了來流馬赫數(shù)為0.65,進氣道初始換算流量為60 kg/s的條件下錘激波經(jīng)過進氣道時壁面的壓力變化和內(nèi)部的三維流動情況。并進一步分析了不同換算流量、超壓比和進氣道中心線幾何形狀對進氣道所受錘激波壓力的影響。結(jié)果表明:

    1) 在馬赫數(shù)為0.65,進氣道初始換算流量為60 kg/s(出口馬赫數(shù)Mae=0.44)的條件下,錘激波以相對進氣道大約306 m/s(Ma=0.9)的速度向上游傳播,相對激波上游來流的速度約481 m/s(Ma=1.42)。在錘激波傳播過程中伴隨三維復雜流動結(jié)構(gòu)和回流的產(chǎn)生。通過記錄進氣道上下壁面各點的壓力隨時間變化情況,發(fā)現(xiàn)進氣道上下壁面所受峰值動態(tài)壓力達到了來流靜壓的2.3~2.65倍。錘激波經(jīng)過進氣道彎道處時,彎道外側(cè)的壁面壓力會大于內(nèi)側(cè),尤其靠近入口的彎折處,上壁面所受的峰值壓力為整個進氣道所受錘激波壓力最高的部位。在錘激波傳出進氣道后,進氣道內(nèi)部流動不會立刻恢復到初始流動,而是經(jīng)歷了由進氣道中部同時向入口和出口流出以及隨后變?yōu)橥耆亓鞯攘鲃訝顟B(tài),經(jīng)過了大約0.185 s后才恢復為正常狀態(tài)。

    2) 通過對不同換算流量條件下錘激波演化過程的對比,發(fā)現(xiàn)較低流量時進氣道所受峰值壓力更大,錘激波相對進氣道傳播速度更快。對比不同超壓比時錘激波演化過程,發(fā)現(xiàn)較高超壓比時進氣道受到峰值壓力更大,錘激波傳播速度更快。

    3) 對具有不同中心線幾何形狀的進氣道進行錘激波計算,表明了中心線形狀對錘激波壓力的大小和分布有重要影響。不同中心線形狀時,S形進氣道前半段的彎道所受的錘激波壓力產(chǎn)生了明顯的改變,曲率較大的地方易產(chǎn)生大的錘激波峰值壓力,而后半段受到的影響較小。

    [1] YOUNG L C, BEAULIEU W D. Review of hammershock pressures in aircraft inlets[J]. Journal of Aircraft, 1975, 12(4): 210-216.

    [2] MARSHALL F L. Prediction of inlet duct overpressures resulting from engine surge[J]. Journal of Aircraft, 1973, 10(5): 274-278.

    [3] KURKOV A P, SOEDER R H, MOSS J E. Investigation of the stall hammershock at the engine inlet[J]. Journal of Aircraft, 1975, 12(4): 198-204.

    [4] BELLMAN D, HUGHES D. The flight investigation of pressure phenomena in the air intake of an F-111A airplane: AIAA-1969-0488[R]. Reston, VA: AIAA, 1969.

    [5] NUGENT J, HOLZMAN J K. Flight-measured inlet pressure transients accompanying engine compressor surges on the F-111A airplane: NASA TN D-7696[R]. Washington, D.C.: NASA, 1974.

    [6] EVANS P J, TRUAX P. YF-16 air induction system design loads associated with engine surge[J]. Journal of Aircraft, 1975, 12(4):205-209.

    [7] HINDASH I, BUSH R, COSNER R. Computational modeling of inlet hammershock wave generation: AIAA-1990-2005[R]. Reston, VA: AIAA, 1990.

    [8] PAYNTER G, MAYER D, TJONNELAND E. Flow stability issues in supersonic inlet flow analyses: AIAA-1993-0290[R]. Reston, VA: AIAA, 1993.

    [9] GOBLE B D, KING S, TERRY J, et al. Inlet hammershock analysis using a 3-D unsteady Euler/Navier-Stokes code: AIAA-1996-2547[R]. Reston, VA: AIAA, 1996.

    [10] YTTERSTRCM A, AXELSON E. Hammershock calculations in the air intake of JAS 39 GFUPEN using dual timestepping: AIAA-1999-3113[R]. Reston, VA: AIAA, 1999.

    [11] BECKER J. Dynamic hammershock effects on the air intake design of supersonic aircraft[M]∥Structures under Shock and Impact Ⅲ. Southampton: WIT Press, 1994: 419-426.

    [12] HSIEH T, WARDLAW A B, COLLINS P, et al. Numerical investigation of unsteady inlet flowfields[J]. AIAA Journal, 1987, 25(1): 75-81.

    [13] GRIDLEY M, SYLVESTER T, TRUAX P. Impact of a probabilistic approach on inlet hammershock design loads: AIAA-1999-2114[R]. Reston, VA: AIAA, 1999.

    [14] 朱宇, 沈天榮. 飛機進氣道錘擊波載荷評估方法研究[J]. 航空發(fā)動機, 2015, 41(3): 6-11.

    ZHU Y, SHEN T R. Evaluation approach of hammer shock loading for aircraft inlet[J]. Aeroengine, 2015, 41(3): 6-11 (in Chinese).

    [15] ANDERSON B H, REDDY D R, KAPOOR K. Study on computing separating flows within a diffusion inlet S-duct[J]. Journal of Propulsion & Power, 1994, 10(5): 661-667.

    [16] NICHOLS R. Calculation of the flow in a circular S-duct inlet: AIAA-1991-0174 [R]. Reston, VA: AIAA, 1991.

    [17] STANLEY R, MOHLER J. Wind-US flow calculations for the M2129 S-duct using structured and unstructured grids: AIAA-2004-0525[R]. Reston, VA: AIAA, 2004.

    [18] DEBIASI M, HERBERG M R, YAN Z, et al. Control of flow separation in S-ducts via flow injection and suction: AIAA-2008-0074[R]. Reston, VA: AIAA, 2008.

    [19] LEE C, BOEDICKER C. Subsonic diffuser design and performance for advanced fighter aircraft: AIAA-1985-3073[R]. Reston, VA: AIAA, 1985.

    猜你喜歡
    進氣道馬赫數(shù)激波
    高馬赫數(shù)激波作用下單模界面的Richtmyer-Meshkov不穩(wěn)定性數(shù)值模擬
    爆炸與沖擊(2024年7期)2024-11-01 00:00:00
    一維非等熵可壓縮微極流體的低馬赫數(shù)極限
    基于AVL-Fire的某1.5L發(fā)動機進氣道優(yōu)化設(shè)計
    基于輔助進氣門的進氣道/發(fā)動機一體化控制
    載荷分布對可控擴散葉型性能的影響
    一種基于聚類分析的二維激波模式識別算法
    航空學報(2020年8期)2020-09-10 03:25:34
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    斜激波入射V形鈍前緣溢流口激波干擾研究
    適于可壓縮多尺度流動的緊致型激波捕捉格式
    The coupling characteristics of supersonic dual inlets for missile①
    丁香欧美五月| 岛国在线观看网站| 一进一出抽搐动态| 亚洲五月天丁香| 欧美zozozo另类| 久久精品综合一区二区三区| 97超视频在线观看视频| 久久婷婷人人爽人人干人人爱| 午夜影院日韩av| 久久人人精品亚洲av| 中文亚洲av片在线观看爽| 一个人看的www免费观看视频| 超碰成人久久| 两性午夜刺激爽爽歪歪视频在线观看| 好看av亚洲va欧美ⅴa在| 国产精品亚洲美女久久久| 91在线观看av| 亚洲av成人一区二区三| 国产激情欧美一区二区| 两个人看的免费小视频| 精品久久久久久久久久久久久| 精品久久久久久成人av| 成熟少妇高潮喷水视频| 欧美另类亚洲清纯唯美| 法律面前人人平等表现在哪些方面| 亚洲av电影不卡..在线观看| 免费看日本二区| 久久久久性生活片| 国内精品久久久久久久电影| 成年免费大片在线观看| 美女 人体艺术 gogo| 在线看三级毛片| 欧美色欧美亚洲另类二区| 久久久久性生活片| 在线看三级毛片| 美女高潮的动态| 久久久久国产一级毛片高清牌| www.精华液| 99国产综合亚洲精品| 亚洲七黄色美女视频| 亚洲性夜色夜夜综合| 国产精品日韩av在线免费观看| 国产三级黄色录像| 久久精品影院6| 极品教师在线免费播放| 久久天躁狠狠躁夜夜2o2o| av女优亚洲男人天堂 | 最近最新免费中文字幕在线| 国产成人精品无人区| 亚洲欧美一区二区三区黑人| 日韩欧美免费精品| 舔av片在线| 亚洲精品中文字幕一二三四区| 母亲3免费完整高清在线观看| 免费观看的影片在线观看| 欧美色视频一区免费| x7x7x7水蜜桃| 成人三级做爰电影| 国产高清videossex| 欧美日本亚洲视频在线播放| 久久精品国产综合久久久| 两性午夜刺激爽爽歪歪视频在线观看| 黄色日韩在线| 99视频精品全部免费 在线 | 久久天躁狠狠躁夜夜2o2o| 两性夫妻黄色片| 中文字幕精品亚洲无线码一区| 两人在一起打扑克的视频| 国产精品 欧美亚洲| 欧美大码av| 国产激情偷乱视频一区二区| 又黄又爽又免费观看的视频| 亚洲欧美日韩东京热| 香蕉国产在线看| 最近最新中文字幕大全免费视频| 曰老女人黄片| 国产精品电影一区二区三区| 亚洲专区字幕在线| 色老头精品视频在线观看| 亚洲欧美日韩东京热| 国产午夜精品论理片| 激情在线观看视频在线高清| 在线观看免费视频日本深夜| 好男人电影高清在线观看| 国内揄拍国产精品人妻在线| 亚洲av美国av| 露出奶头的视频| 亚洲欧美日韩卡通动漫| 亚洲五月天丁香| 99riav亚洲国产免费| 俺也久久电影网| 最近最新中文字幕大全免费视频| 一个人免费在线观看的高清视频| 日韩 欧美 亚洲 中文字幕| 国产伦精品一区二区三区视频9 | 久久人妻av系列| 99久久成人亚洲精品观看| 欧美国产日韩亚洲一区| 99久久精品一区二区三区| 欧美不卡视频在线免费观看| 999久久久精品免费观看国产| 成年女人看的毛片在线观看| 亚洲成人久久性| 亚洲精品一区av在线观看| av国产免费在线观看| 国产精品免费一区二区三区在线| 亚洲精品美女久久av网站| 日本免费a在线| 中国美女看黄片| 国产真实乱freesex| av片东京热男人的天堂| 婷婷丁香在线五月| 两个人看的免费小视频| 欧美乱色亚洲激情| 亚洲欧美精品综合一区二区三区| 两个人的视频大全免费| www日本黄色视频网| 国产高清视频在线观看网站| 香蕉丝袜av| 桃红色精品国产亚洲av| 国产成人aa在线观看| 国产69精品久久久久777片 | 国产一区二区在线观看日韩 | 一个人看的www免费观看视频| 国产爱豆传媒在线观看| 精华霜和精华液先用哪个| 18禁国产床啪视频网站| 国产午夜精品久久久久久| 日韩欧美免费精品| 精品一区二区三区视频在线 | 色视频www国产| 国产综合懂色| 99国产精品99久久久久| av福利片在线观看| 免费看a级黄色片| 午夜免费观看网址| 听说在线观看完整版免费高清| 这个男人来自地球电影免费观看| 亚洲成av人片免费观看| 成人鲁丝片一二三区免费| 一本精品99久久精品77| 欧美高清成人免费视频www| 亚洲电影在线观看av| 美女大奶头视频| 精品国产三级普通话版| 国产精品九九99| 国产成人影院久久av| 夜夜夜夜夜久久久久| 国产探花在线观看一区二区| 欧美日韩国产亚洲二区| 99视频精品全部免费 在线 | 99久久综合精品五月天人人| 国产高清三级在线| 一级毛片精品| 免费无遮挡裸体视频| 日韩欧美 国产精品| av欧美777| 亚洲av电影在线进入| 日韩欧美免费精品| 午夜激情欧美在线| 黑人欧美特级aaaaaa片| 国产aⅴ精品一区二区三区波| 亚洲一区二区三区不卡视频| 小说图片视频综合网站| 村上凉子中文字幕在线| 欧美黄色片欧美黄色片| 日韩av在线大香蕉| 亚洲人与动物交配视频| 亚洲欧洲精品一区二区精品久久久| 又黄又粗又硬又大视频| 老汉色∧v一级毛片| 丰满的人妻完整版| 亚洲九九香蕉| 国产淫片久久久久久久久 | 人妻久久中文字幕网| 狠狠狠狠99中文字幕| 91久久精品国产一区二区成人 | 欧美色视频一区免费| 两人在一起打扑克的视频| 真人做人爱边吃奶动态| 国产毛片a区久久久久| 色哟哟哟哟哟哟| 人人妻,人人澡人人爽秒播| avwww免费| 搡老妇女老女人老熟妇| 99国产精品一区二区蜜桃av| 国产成人一区二区三区免费视频网站| 成人18禁在线播放| 午夜日韩欧美国产| 午夜福利在线在线| 久久精品综合一区二区三区| 女警被强在线播放| 免费观看精品视频网站| 国产亚洲精品久久久久久毛片| 久久国产乱子伦精品免费另类| 美女高潮喷水抽搐中文字幕| 国产成人精品久久二区二区91| 欧美色欧美亚洲另类二区| 亚洲 欧美一区二区三区| 变态另类成人亚洲欧美熟女| 成人国产一区最新在线观看| 在线视频色国产色| 久久精品综合一区二区三区| 青草久久国产| 熟女少妇亚洲综合色aaa.| 亚洲av片天天在线观看| 久久精品91蜜桃| 国产精品久久久久久精品电影| 18禁美女被吸乳视频| 非洲黑人性xxxx精品又粗又长| 高清在线国产一区| 在线十欧美十亚洲十日本专区| 久久久久亚洲av毛片大全| 免费电影在线观看免费观看| 亚洲精华国产精华精| 又爽又黄无遮挡网站| 国产精品亚洲美女久久久| 午夜免费成人在线视频| 九九热线精品视视频播放| 最近在线观看免费完整版| 级片在线观看| 一二三四在线观看免费中文在| 色噜噜av男人的天堂激情| 久久中文字幕人妻熟女| 国产精品美女特级片免费视频播放器 | 在线免费观看不下载黄p国产 | 国产午夜福利久久久久久| 在线永久观看黄色视频| 国产成人精品久久二区二区91| 欧美3d第一页| 国产免费男女视频| 人人妻人人澡欧美一区二区| 亚洲国产欧美人成| 黄色片一级片一级黄色片| 日韩欧美在线乱码| 国产精品综合久久久久久久免费| 九色成人免费人妻av| 男人的好看免费观看在线视频| 国产人伦9x9x在线观看| 两性午夜刺激爽爽歪歪视频在线观看| 日韩欧美国产一区二区入口| 长腿黑丝高跟| 久久午夜亚洲精品久久| 男女那种视频在线观看| 国产极品精品免费视频能看的| 国产成年人精品一区二区| 18禁裸乳无遮挡免费网站照片| 色尼玛亚洲综合影院| 亚洲午夜精品一区,二区,三区| 国产三级在线视频| 最近视频中文字幕2019在线8| 色视频www国产| 亚洲一区二区三区不卡视频| 亚洲av成人不卡在线观看播放网| 国产精品日韩av在线免费观看| 国产精品美女特级片免费视频播放器 | 国产精品1区2区在线观看.| 色综合亚洲欧美另类图片| a在线观看视频网站| 久久久久精品国产欧美久久久| 亚洲国产高清在线一区二区三| 国产真人三级小视频在线观看| 国产伦人伦偷精品视频| 狂野欧美激情性xxxx| 亚洲自偷自拍图片 自拍| 美女大奶头视频| 狠狠狠狠99中文字幕| 99热这里只有精品一区 | 最好的美女福利视频网| 国产精品自产拍在线观看55亚洲| 免费在线观看成人毛片| 无限看片的www在线观看| 午夜福利在线观看免费完整高清在 | 欧美在线黄色| 黄色视频,在线免费观看| 国产成人一区二区三区免费视频网站| 一级毛片精品| 成人国产一区最新在线观看| 亚洲av电影在线进入| 夜夜夜夜夜久久久久| 91在线观看av| 男人舔女人下体高潮全视频| 国产精品久久视频播放| 午夜久久久久精精品| 免费在线观看影片大全网站| 男女做爰动态图高潮gif福利片| 欧美三级亚洲精品| 波多野结衣高清无吗| 久久中文字幕人妻熟女| 99精品久久久久人妻精品| 国产av麻豆久久久久久久| 最新在线观看一区二区三区| 美女被艹到高潮喷水动态| 男人舔奶头视频| 久久中文字幕人妻熟女| 欧美极品一区二区三区四区| 国产精品99久久99久久久不卡| 他把我摸到了高潮在线观看| 午夜免费激情av| 小蜜桃在线观看免费完整版高清| 女人高潮潮喷娇喘18禁视频| 中文字幕熟女人妻在线| 久久久久性生活片| 看黄色毛片网站| 久久久久久九九精品二区国产| 国产午夜精品久久久久久| 亚洲人成伊人成综合网2020| 在线十欧美十亚洲十日本专区| 他把我摸到了高潮在线观看| 国产欧美日韩一区二区精品| 亚洲天堂国产精品一区在线| 亚洲 欧美 日韩 在线 免费| 久久精品国产清高在天天线| 亚洲色图 男人天堂 中文字幕| 麻豆一二三区av精品| 少妇人妻一区二区三区视频| 在线观看美女被高潮喷水网站 | 无限看片的www在线观看| 在线观看一区二区三区| 又大又爽又粗| 日韩大尺度精品在线看网址| 日韩欧美国产一区二区入口| 制服丝袜大香蕉在线| 色尼玛亚洲综合影院| 1000部很黄的大片| 免费大片18禁| 亚洲精品乱码久久久v下载方式 | 动漫黄色视频在线观看| 午夜福利免费观看在线| 亚洲欧美精品综合久久99| 亚洲精品一卡2卡三卡4卡5卡| 色吧在线观看| 99久久国产精品久久久| 亚洲国产精品合色在线| 美女高潮喷水抽搐中文字幕| 无人区码免费观看不卡| 国产精品99久久久久久久久| 97人妻精品一区二区三区麻豆| 这个男人来自地球电影免费观看| 99热精品在线国产| 欧美日本亚洲视频在线播放| 99国产精品一区二区蜜桃av| 国产 一区 欧美 日韩| 大型黄色视频在线免费观看| 免费大片18禁| 黄色丝袜av网址大全| 两个人的视频大全免费| 国产精品av久久久久免费| 精品久久蜜臀av无| 在线观看66精品国产| 国产午夜精品久久久久久| 久久久久国产一级毛片高清牌| 在线观看66精品国产| 一二三四社区在线视频社区8| 色视频www国产| 又紧又爽又黄一区二区| 久久久久久大精品| 熟女少妇亚洲综合色aaa.| 日本黄色视频三级网站网址| 99精品欧美一区二区三区四区| 热99re8久久精品国产| 激情在线观看视频在线高清| 亚洲精品美女久久av网站| 又爽又黄无遮挡网站| 99视频精品全部免费 在线 | 欧美国产日韩亚洲一区| 99国产精品一区二区三区| 午夜精品一区二区三区免费看| 久久久国产成人免费| 国产男靠女视频免费网站| 国产亚洲精品av在线| 99久久成人亚洲精品观看| 欧美国产日韩亚洲一区| 亚洲午夜精品一区,二区,三区| 日韩欧美三级三区| 三级毛片av免费| 免费av不卡在线播放| 床上黄色一级片| 午夜福利成人在线免费观看| 成人特级黄色片久久久久久久| 90打野战视频偷拍视频| 国产精品 欧美亚洲| 亚洲成人久久性| 欧美黄色片欧美黄色片| 最近最新免费中文字幕在线| 看免费av毛片| 免费看a级黄色片| 又大又爽又粗| 国产激情欧美一区二区| 成人三级黄色视频| 高清毛片免费观看视频网站| 欧美日本视频| 91麻豆精品激情在线观看国产| 国产日本99.免费观看| 麻豆成人av在线观看| a级毛片a级免费在线| 欧美国产日韩亚洲一区| 97碰自拍视频| 亚洲欧美日韩卡通动漫| 窝窝影院91人妻| www.www免费av| 99精品久久久久人妻精品| 757午夜福利合集在线观看| 综合色av麻豆| 亚洲人成网站高清观看| 国产一区二区在线观看日韩 | 亚洲18禁久久av| 久久久久国产一级毛片高清牌| 12—13女人毛片做爰片一| 日日干狠狠操夜夜爽| 一级毛片高清免费大全| 精品国产美女av久久久久小说| 无人区码免费观看不卡| 网址你懂的国产日韩在线| 亚洲精品美女久久久久99蜜臀| 老司机午夜福利在线观看视频| 夜夜看夜夜爽夜夜摸| 偷拍熟女少妇极品色| 99热这里只有是精品50| 五月伊人婷婷丁香| 香蕉国产在线看| 色老头精品视频在线观看| 成人特级av手机在线观看| 蜜桃久久精品国产亚洲av| av天堂在线播放| 亚洲五月天丁香| 国产伦在线观看视频一区| 国产成年人精品一区二区| 熟女少妇亚洲综合色aaa.| www.自偷自拍.com| a在线观看视频网站| 亚洲国产精品合色在线| 国产精品爽爽va在线观看网站| 两个人看的免费小视频| 99久久综合精品五月天人人| 热99re8久久精品国产| 一个人看的www免费观看视频| 国产一区二区在线av高清观看| 日韩欧美在线二视频| 欧美av亚洲av综合av国产av| 真人一进一出gif抽搐免费| 一二三四在线观看免费中文在| 舔av片在线| 精品国产乱子伦一区二区三区| 91av网一区二区| 午夜激情欧美在线| 亚洲成人免费电影在线观看| 色播亚洲综合网| 91在线观看av| 国产午夜精品久久久久久| 欧美大码av| 久久精品夜夜夜夜夜久久蜜豆| 国产69精品久久久久777片 | 亚洲欧美日韩高清专用| 麻豆成人av在线观看| 国产伦人伦偷精品视频| 麻豆一二三区av精品| 88av欧美| 欧美大码av| 国产精品av视频在线免费观看| 欧美丝袜亚洲另类 | 久久国产乱子伦精品免费另类| 成人性生交大片免费视频hd| 午夜两性在线视频| 狠狠狠狠99中文字幕| 欧美日本视频| 亚洲精品国产精品久久久不卡| 天堂影院成人在线观看| 日本 av在线| 特大巨黑吊av在线直播| 国产三级中文精品| 午夜精品久久久久久毛片777| 欧美日韩综合久久久久久 | 精品久久久久久久毛片微露脸| 99久久成人亚洲精品观看| 97人妻精品一区二区三区麻豆| 操出白浆在线播放| 熟女电影av网| 免费搜索国产男女视频| 搞女人的毛片| 亚洲性夜色夜夜综合| 看免费av毛片| 精品乱码久久久久久99久播| 午夜福利成人在线免费观看| 日本与韩国留学比较| 亚洲人与动物交配视频| 人妻夜夜爽99麻豆av| 黄片大片在线免费观看| 亚洲欧美精品综合一区二区三区| 日日干狠狠操夜夜爽| 熟妇人妻久久中文字幕3abv| 国产黄a三级三级三级人| 淫妇啪啪啪对白视频| 1024香蕉在线观看| 午夜福利成人在线免费观看| 国产av在哪里看| 久久这里只有精品19| 国产高潮美女av| 一本一本综合久久| 亚洲国产高清在线一区二区三| 最近最新免费中文字幕在线| h日本视频在线播放| 一个人观看的视频www高清免费观看 | 精品国内亚洲2022精品成人| 国产91精品成人一区二区三区| 99久久成人亚洲精品观看| 国产欧美日韩一区二区三| 欧美日韩亚洲国产一区二区在线观看| 婷婷精品国产亚洲av在线| 免费看十八禁软件| 婷婷亚洲欧美| 亚洲av第一区精品v没综合| 国产伦一二天堂av在线观看| a级毛片a级免费在线| 香蕉久久夜色| 极品教师在线免费播放| 9191精品国产免费久久| 国产精品九九99| 亚洲人成网站高清观看| 深夜精品福利| 国内精品久久久久精免费| 男女床上黄色一级片免费看| 久久国产精品人妻蜜桃| 精品乱码久久久久久99久播| 真实男女啪啪啪动态图| 免费av毛片视频| 久久午夜综合久久蜜桃| 国产精品久久久人人做人人爽| 亚洲欧美精品综合久久99| 美女免费视频网站| 中文字幕人妻丝袜一区二区| 夜夜躁狠狠躁天天躁| www.999成人在线观看| 免费在线观看日本一区| 久久久久久久久久黄片| 两个人看的免费小视频| 长腿黑丝高跟| 麻豆一二三区av精品| 听说在线观看完整版免费高清| 国产亚洲精品一区二区www| 一个人免费在线观看电影 | 亚洲九九香蕉| 国产免费男女视频| 欧美在线一区亚洲| 国产精品一区二区精品视频观看| 给我免费播放毛片高清在线观看| 午夜视频精品福利| 中文字幕精品亚洲无线码一区| 亚洲成人精品中文字幕电影| 人妻久久中文字幕网| 草草在线视频免费看| 亚洲天堂国产精品一区在线| 十八禁网站免费在线| 国产1区2区3区精品| 天天添夜夜摸| 亚洲最大成人中文| 91久久精品国产一区二区成人 | 在线国产一区二区在线| 亚洲性夜色夜夜综合| 国产精品乱码一区二三区的特点| 久久久精品欧美日韩精品| 99国产精品一区二区三区| a级毛片在线看网站| 久久香蕉精品热| bbb黄色大片| 日韩欧美精品v在线| 人妻久久中文字幕网| 国产精品日韩av在线免费观看| 亚洲欧美精品综合一区二区三区| xxxwww97欧美| 黄片小视频在线播放| 九色国产91popny在线| 国产伦人伦偷精品视频| 国产高潮美女av| 色尼玛亚洲综合影院| 亚洲九九香蕉| 1000部很黄的大片| 日本三级黄在线观看| 听说在线观看完整版免费高清| 亚洲精品在线美女| 熟女少妇亚洲综合色aaa.| a级毛片a级免费在线| 亚洲国产高清在线一区二区三| 中出人妻视频一区二区| 久久久久久国产a免费观看| 欧美日韩综合久久久久久 | 此物有八面人人有两片| 国产aⅴ精品一区二区三区波| 性色av乱码一区二区三区2| 欧美一区二区精品小视频在线| 黄色丝袜av网址大全| 一个人免费在线观看电影 | 99久久综合精品五月天人人| 琪琪午夜伦伦电影理论片6080| 在线观看免费视频日本深夜| 久久草成人影院| 欧美一区二区国产精品久久精品| 亚洲人成网站在线播放欧美日韩| 亚洲第一欧美日韩一区二区三区| 国产三级黄色录像| 网址你懂的国产日韩在线| 亚洲成a人片在线一区二区| www.999成人在线观看| 成人高潮视频无遮挡免费网站| 成人无遮挡网站| 欧美一区二区精品小视频在线| 精品久久久久久成人av| 亚洲性夜色夜夜综合| 精品国产亚洲在线| 亚洲成a人片在线一区二区| 国产蜜桃级精品一区二区三区| 国产激情偷乱视频一区二区| 男人的好看免费观看在线视频| 男人舔女人下体高潮全视频| 国产熟女xx|