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

    高超聲速進氣道強制轉(zhuǎn)捩流動的大渦模擬

    2016-11-14 00:41:17趙曉慧鄧小兵毛枚良楊偉趙慧勇
    航空學報 2016年8期
    關(guān)鍵詞:大渦進氣道邊界層

    趙曉慧, 鄧小兵,*, 毛枚良, 楊偉, 趙慧勇

    1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室, 綿陽 621000 2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 綿陽 621000 3.中國空氣動力研究與發(fā)展中心 高超聲速沖壓發(fā)動機重點實驗室, 綿陽 621000

    ?

    高超聲速進氣道強制轉(zhuǎn)捩流動的大渦模擬

    趙曉慧1, 鄧小兵1,*, 毛枚良2, 楊偉1, 趙慧勇3

    1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室, 綿陽621000 2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 綿陽621000 3.中國空氣動力研究與發(fā)展中心 高超聲速沖壓發(fā)動機重點實驗室, 綿陽621000

    吸氣式高超聲速飛行器常在進氣道邊界層內(nèi)布置粗糙顆?;驕u流發(fā)生器強制流動轉(zhuǎn)捩為湍流以確保發(fā)動機正常啟動。為了清晰認識強制轉(zhuǎn)捩過程,采用隱式大渦模擬方法,對強制轉(zhuǎn)捩問題開展了數(shù)值模擬研究。針對鉆石形和斜坡形渦流發(fā)生器,計算得到渦流發(fā)生器誘導(dǎo)的流動結(jié)構(gòu),顯示出強制轉(zhuǎn)捩流動由渦流發(fā)生器產(chǎn)生的反向旋轉(zhuǎn)流向渦對主導(dǎo)。擾動沿流向增長和發(fā)展,導(dǎo)致流向渦對以偶模式或奇模式失穩(wěn),偶模式失穩(wěn)產(chǎn)生對稱形式的渦對破碎,而奇模式失穩(wěn)則導(dǎo)致非對稱(彎曲)形式的渦對破碎。流向渦對破碎后產(chǎn)生一系列發(fā)卡渦并最終促使邊界層轉(zhuǎn)捩為湍流。最后就計算網(wǎng)格和數(shù)值耗散對隱式大渦模擬結(jié)果的影響以及計算的收斂性進行了討論。

    高超聲速; 進氣道; 強制轉(zhuǎn)捩; 渦流發(fā)生器; 大渦模擬

    吸氣式高超聲速飛行器通常巡航在高度比較高,而密度則相對較低的大氣環(huán)境中,此時飛行雷諾數(shù)相對較小,自然轉(zhuǎn)捩發(fā)生的位置通常超出了進氣道長度。為了滿足發(fā)動機設(shè)計對邊界層湍流流態(tài)的要求,需要在進氣道采用強制轉(zhuǎn)捩措施。2004年,Hyper-X (X-43A)采用渦流發(fā)生器在飛行試驗中成功實現(xiàn)了強制轉(zhuǎn)捩[1],測試結(jié)果顯示,對于馬赫數(shù)7和10的飛行條件,飛行器前體需要采用粗糙帶來保證進氣道上游的全湍流狀態(tài)。

    渦流發(fā)生器強制轉(zhuǎn)捩的作用機理圖像尚未完全清晰,最早的一些研究者試圖用Tollmien-Schlichting (T-S)波解釋粗糙顆粒強制轉(zhuǎn)捩機理[2],Reshotko和Tumin[3]則認為瞬時增長在粗糙顆粒強制轉(zhuǎn)捩中起了重要作用。數(shù)值計算結(jié)果顯示強制轉(zhuǎn)捩過程與渦流發(fā)生器產(chǎn)生的流向渦對的發(fā)展有關(guān),這與G?rtler渦的二次失穩(wěn)相似,在文獻[4]中這種粗糙顆粒引起的流向反向旋轉(zhuǎn)的渦對也被稱作G?rtler渦。Li等[5-6]在研究低速和高速流動G?rtler渦時指出其二次失穩(wěn)主要有奇模式(Odd mode)和偶模式(Even mode)兩種,分別導(dǎo)致了彎曲(Sinuous)和對稱(Varicose)兩種不同形式的流向渦對破碎,并且其不穩(wěn)定模式與G?rtler渦波長有關(guān)。Tullio等[7]利用線性穩(wěn)定性理論以及求解Navier-Stokes方程研究超聲速平板利用顆粒強制轉(zhuǎn)捩時也指出,對稱和彎曲是兩種最容易出現(xiàn)的不穩(wěn)定形態(tài)。

    高超聲速邊界層的渦流發(fā)生器強制轉(zhuǎn)捩過程受渦流發(fā)生器的形狀、尺寸、分布和環(huán)境噪聲等諸多因素的影響。Berry等[8-10]對五種形狀的渦流發(fā)生器作了試驗比較,發(fā)現(xiàn)鉆石形和斜坡形渦流發(fā)生器具有較好的轉(zhuǎn)捩觸發(fā)效果。趙慧勇等[11]在中國空氣動力研究與發(fā)展中心的FL-31(直徑0.5 m)超聲速風洞中,開展了鉆石形和斜坡形渦流發(fā)生器強制轉(zhuǎn)捩研究,給出了鉆石形渦流發(fā)生器觸發(fā)轉(zhuǎn)捩的有效高度。Danehy等[12]采用平面激光誘導(dǎo)熒光(Plane Laser-Induced Fluorescence, PLIF) 技術(shù)研究了平板表面半球凸起物對邊界層的影響,岡敦殿等[13]采用基于納米示蹤的平面激光散射技術(shù)(Nano-tracer Planar Laser Scattering, NPLS)進行了超聲速平板上直立圓臺凸起物繞流流場研究,張慶虎[14]采用NPLS技術(shù)對三角楔渦流發(fā)生器產(chǎn)生的流動結(jié)構(gòu)給過較為細致的展示,這些研究都給出邊界層中突起物尾跡某一平面上的精細結(jié)構(gòu),對研究突起物尾跡的發(fā)展有很大幫助。Borg和Schneider[15]在靜聲風洞中開展了X-51前體強制轉(zhuǎn)捩研究,指出來流噪聲會使轉(zhuǎn)捩位置提前;鄧小兵[16]采用隱式大渦模擬(Implicit Large Eddy Simulation, ILES)方法,通過在計算來流中加噪聲信號來模擬試驗條件,發(fā)現(xiàn)相同的現(xiàn)象。這表明,要通過試驗準確預(yù)測轉(zhuǎn)捩位置需采用靜聲風洞以消除風洞噪聲的影響。

    相比于試驗研究,數(shù)值計算可以獲得更為全面的流場信息,而且不受風洞噪聲的影響。研究強制轉(zhuǎn)捩問題采用的數(shù)值模擬方法主要有直接數(shù)值模擬(Direct Numerical Simulation, DNS),大渦模擬(Large Eddy Simulation, LES)和雷諾平均Navier-Stokes(RANS)方程方法。其中,DNS方法最準確,Marxen[17]、Iyer[18]和Duan[19-20]等利用DNS對突起物強制轉(zhuǎn)捩流動進行研究,獲得了清晰的流動結(jié)構(gòu)和轉(zhuǎn)捩圖像。DNS方法的缺點是所需的計算資源過多,在當前的計算機硬件條件下還難以推廣應(yīng)用。涂國華等[21]曾經(jīng)利用RANS方法對強制轉(zhuǎn)捩流動進行模擬,結(jié)果表明RANS方法難以預(yù)測轉(zhuǎn)捩位置,且對轉(zhuǎn)捩后壓縮面上熱流分布的預(yù)測也不理想。周玲等[22-23]改進了k-ω-γ轉(zhuǎn)捩模式,并對高超聲速飛行器進氣道前體邊界層強制轉(zhuǎn)捩進行數(shù)值計算,指出改進模型對轉(zhuǎn)捩位置有較好的預(yù)測能力。然而采用RANS方法難以準確捕捉強制轉(zhuǎn)捩的流動結(jié)構(gòu)。LES方法能夠以遠低于DNS的計算開銷給出較為準確的非定常流動結(jié)構(gòu)。Sayadi和Moin[24]利用大渦模擬方法對馬赫數(shù)0.2可壓縮平板邊界層的H-type和K-type轉(zhuǎn)捩進行了模擬,指出渦黏性在層流和轉(zhuǎn)捩區(qū)的存在會阻礙擾動的發(fā)展,不利于轉(zhuǎn)捩預(yù)測。ILES不加入亞格子模型,在層流區(qū)不會引入亞格子耗散,能夠更準確地模擬轉(zhuǎn)捩區(qū)的擾動增長。Orlik和Fedioun[25]采用基于5階WENO格式的隱式大渦模擬方法計算了射流控制的高超聲速邊界層強制轉(zhuǎn)捩流動,結(jié)果顯示ILES能夠捕捉到豐富的流動細節(jié)。本文重點關(guān)注渦流發(fā)生器誘導(dǎo)強制轉(zhuǎn)捩過程的模擬,因而采用ILES是合適的。

    通常認為,大渦模擬方法需要采用高階精度格式,鄧小兵和金玲[26]通過對湍流數(shù)值模擬中誤差傳播和發(fā)展的動力學行為的分析,認為二階精度為基礎(chǔ)的數(shù)值方法可以應(yīng)用于湍流大渦模擬。本文采用二階有限體積隱式大渦模擬方法,對鉆石形和斜坡形渦流發(fā)生器激發(fā)的高超聲速進氣道邊界層轉(zhuǎn)捩現(xiàn)象開展數(shù)值模擬研究,計算得到了流向渦失穩(wěn)的兩種基本模式:偶模式與奇模式,同時給出了脈動量沿流向的增長規(guī)律,最后對數(shù)值耗散和網(wǎng)格分辨率共同作用下ILES方法的收斂性作了探討。

    1 研究方法

    1.1隱式大渦模擬方法

    采用基于二階有限體積的ILES方法研究強制轉(zhuǎn)捩問題。由于ILES采用數(shù)值黏性模擬亞格子耗散,數(shù)值耗散的控制就十分重要。對流項空間離散格式的數(shù)值耗散對大渦模擬結(jié)果有明顯的影響,過大的數(shù)值耗散往往導(dǎo)致能譜過度耗散,抑制不穩(wěn)定波的初始增長,所以數(shù)值耗散的控制對轉(zhuǎn)捩位置的準確預(yù)測至關(guān)重要。采用基于MUSCL( Monotone Upstream-centered Scheme for Conservation Laws)重構(gòu)的ROE格式,在光滑區(qū),該重構(gòu)表達式退化為

    (1)

    式中:qL,i+1/2和qR,i+1/2分別為i和i+1單元界面變量左側(cè)和右側(cè)重構(gòu)值;當κ=1/3時為常用的三階格式,此時格式對ILES是過耗散的,κ=1時格式變?yōu)槎A中心格式,是無耗散的,通過調(diào)整κ可以將數(shù)值耗散控制在合理水平。高超聲速流動中存在強激波,計算中對差商Δ-、Δ+采用min-mod限制器保證計算的魯棒性[27]。黏性項離散采用二階中心格式。非定常時間推進采用雙時間步方法[28],其中物理時間推進采用二階歐拉后差格式,隱式時間推進采用LU-SGS方法[29]。

    1.2算例驗證

    Duan等[30]開展了來流馬赫數(shù)Ma=2.99的超聲速湍流平板邊界層DNS研究,這里用此算例驗證所采用軟件和方法的有效性。

    算例來流密度ρ∞=0.089 1 kg/m3,來流溫度T∞=218.2 K。計算域尺寸為:流向69 mm,展向35 mm,法向110 mm。計算網(wǎng)格單元量約900萬,其中流向dx+=12,展向dy+=4.5,法向第一層網(wǎng)格dz+=0.3。壁面邊界條件為無滑移等溫壁,壁面溫度Tw=567.3 K。計算入口邊界層厚度為7.2 mm,入口邊界采用時間平均剖面疊加瞬時脈動的方法,其中時間平均剖面采用RANS計算的結(jié)果,瞬時脈動采用RRM(Rescaling-Recycling Method)[31-32]將下游x=52 mm位置處的瞬時脈動按照邊界層相似關(guān)系重新標定后疊加到入口邊界。

    物理模型采用基于改進延遲脫體渦模擬(Improved Delayed Detached-Eddy Simulation, IDDES)[33]方法的壁面?;鬁u模擬(Wall-Modeled LES, WMLES),其中湍流雷諾應(yīng)力模型采用一方程Spalart-Allmaras(S-A)模型[34],MUSCL重構(gòu)耗散參數(shù)κ=0.8。

    圖1為計算得到的邊界層平均速度剖面(時間平均加展向平均)。模擬結(jié)果與理論曲線以及文獻[30]DNS模擬結(jié)果符合較好,這表明本文所用的軟件和方法是可靠的。

    圖1 超聲速平板邊界層平均速度剖面計算結(jié)果Fig.1 Computational results of averaged velocity profile in supersonic flat boundary layer

    2 模型與網(wǎng)格說明

    計算采用的進氣道模型是在趙慧勇等[11]試驗?zāi)P突A(chǔ)上簡化得到的,如圖2所示。進氣道為分段壓縮結(jié)構(gòu),三個壓縮面均與水平方向有一定夾角。水平方向從前緣指向下游定義為x正方向,豎直方向定義為y向,渦流發(fā)生器(Vortex Generators, VG)凸起方向大致為y的負方向,展向定義為z方向。z向取了試驗?zāi)P椭芯€附近的一小段,x向長度約300 mm。渦流發(fā)生器采用了鉆石形和斜坡形兩種,后緣大致放在x=88 mm的位置(當?shù)剡吔鐚雍穸燃s為1.2 mm)。計算域展向尺寸以及模擬渦流發(fā)生器的個數(shù)和高度見表1。

    圖2 進氣道與渦流發(fā)生器模型及其網(wǎng)格(G1)Fig.2 Models and grids of inlet and vortex generators (G1)

    表1 計算模型

    Table 1 Computational models

    CaseVGshapeComputationalwidth/mmVGnumberVGheight/mm1Diamond6.021.02Ramp5.730.5

    考慮計算量,本文對鉆石形的渦流發(fā)生器只模擬了兩顆,對斜坡形的模擬了三顆。展向采用周期邊界條件。為了更好地分辨關(guān)鍵區(qū)域的流動結(jié)構(gòu),加密了渦流發(fā)生器附近和尾跡等區(qū)域的網(wǎng)格,并采用G1和G2粗細兩套網(wǎng)格計算考察了網(wǎng)格分辨率對計算結(jié)果的影響,這兩套網(wǎng)格在包含主要流動結(jié)構(gòu)的區(qū)域配置見表2。其中dx+為基于壁面單位的流向網(wǎng)格尺度,dz+為展向,dy+為法向,dy2+為法向第一層網(wǎng)格尺度。按照文獻[35]給出的LES網(wǎng)格分辨率準則,表2給出的粗網(wǎng)格G1已經(jīng)滿足了LES模擬要求。

    表2 計算網(wǎng)格

    3 數(shù)值計算和分析

    3.1計算條件

    根據(jù)趙慧勇等的試驗[11],選用馬赫數(shù)6試驗中的兩組試驗參數(shù)作為計算工況。來流馬赫數(shù)Ma=5.96,迎角α=1°,來流壓力p∞,來流溫度T∞以及第一、第二、第三壓縮面的物面溫度Tw1、Tw2、Tw3見表3。

    表3 計算條件

    3.2渦流發(fā)生器產(chǎn)生的渦結(jié)構(gòu)

    利用大渦模擬,對兩個工況進行模擬,獲取了非定常流場結(jié)構(gòu)。渦流發(fā)生器周圍的渦結(jié)構(gòu)如圖3 所示。其中Q值定義為

    (2)

    式中:Ω為基于特征速度為來流聲速和特征長度為1 mm的無量綱渦強;S為變形率張量。

    圖3 渦流發(fā)生器附近的流線和流場Q等值面Fig.3 Streamlines and Q iso-surfaces of flow around vortex generators

    根據(jù)工況1的瞬時流場可以看出,在鉆石形渦流發(fā)生器前緣根部產(chǎn)生了一個“馬蹄”渦,同時在后緣產(chǎn)生一對反向旋轉(zhuǎn)的流向渦。后緣產(chǎn)生的這對渦,其旋轉(zhuǎn)方向與外側(cè)的“馬蹄”渦反向,而位置比外側(cè)的“馬蹄”渦要遠離壁面一些,其強度也高于“馬蹄”渦。工況2的瞬時流場顯示,斜坡形的渦流發(fā)生器產(chǎn)生的渦結(jié)構(gòu)相對簡單,其主要結(jié)構(gòu)是渦流發(fā)生器側(cè)緣產(chǎn)生的一對反向旋轉(zhuǎn)的流向渦,這對渦的結(jié)構(gòu)、轉(zhuǎn)向與工況1中鉆石形渦流發(fā)生器后緣產(chǎn)生的流向渦相似。由于渦流發(fā)生器形狀、高度和分布密度不同,這里兩個工況中,鉆石形渦流發(fā)生器產(chǎn)生的流向渦更強,斜坡形渦流發(fā)生器產(chǎn)生的渦對更密集。

    3.3強制轉(zhuǎn)捩過程分析

    圖4顯示出鉆石形渦流發(fā)生器下游的主導(dǎo)結(jié)構(gòu)是尾緣產(chǎn)生的流向渦對。這對流向渦相互抬升,并在第二壓縮面上破碎并脫落出一系列的“發(fā)卡”渦,“發(fā)卡”渦很大程度上加劇了流動內(nèi)層和外層的流體交換,破壞了邊界層的穩(wěn)定性,促成轉(zhuǎn)捩。圖5顯示出斜坡形渦流發(fā)生器主要流向渦對結(jié)構(gòu)的發(fā)展,與鉆石形渦流發(fā)生器強制轉(zhuǎn)捩相似,斜坡形渦流發(fā)生器產(chǎn)生的流向渦對在第二壓縮面破碎脫落產(chǎn)生“發(fā)卡”渦,促成轉(zhuǎn)捩,明顯有別于工況1的是,工況2中流動結(jié)構(gòu)有展向的搖擺和彎曲,搖擺的頻率與“發(fā)卡”渦脫落的頻率是相關(guān)的。計算結(jié)果顯示,渦流發(fā)生器促使流動轉(zhuǎn)捩過程是其產(chǎn)生的流向渦對主導(dǎo)的,流向渦對結(jié)構(gòu)失穩(wěn)、破碎產(chǎn)生一系列“發(fā)卡”渦,使得邊界層流動轉(zhuǎn)捩為湍流。

    圖4 鉆石形渦流發(fā)生器強制轉(zhuǎn)捩瞬時流場Q等值面圖(G1, κ=0.9)Fig.4 Q iso-surfaces of forced transition instantaneous flow by diamond shaped vortex generators (G1, κ=0.9)

    圖5 斜坡形渦流發(fā)生器強制轉(zhuǎn)捩瞬時流場Q等值面圖(G2, κ=0.9)Fig.5 Q iso-surfaces of forced transition instantaneous flow by ramp shaped vortex generators (G2, κ=0.9)

    圖6 沿流向渦物理量脈動的均方根Fig.6 Root-mean-square of variable fluctuations along a stream-wise vortex

    兩種渦流發(fā)生器強制轉(zhuǎn)捩過程都與流向渦對結(jié)構(gòu)的失穩(wěn)有關(guān)。圖7為瞬時流場中壓力脈動p′沿流向渦對結(jié)構(gòu)兩個渦心的分布。圖7(a)顯示工況1瞬時流場中壓力脈動相位相同,表明工況1流向渦對失穩(wěn)由偶模式主導(dǎo),這種模式導(dǎo)致流向渦對對稱破碎;圖7(b)顯示工況2瞬時流場中壓力脈動相位相反,表明工況2流向渦對失穩(wěn)由奇模式主導(dǎo),這種模式導(dǎo)致流向渦對彎曲破碎。

    圖7 瞬時流場沿流向渦對兩個渦心的壓力脈動分布Fig.7 Pressure fluctuation of instantaneous field distribution along two vortex centers of counter-rotating vortices

    渦流發(fā)生器的幾何參數(shù)和分布密度會影響流向渦對結(jié)構(gòu)不穩(wěn)定增長的主導(dǎo)模式、頻率、破碎位置等,這些參數(shù)的具體影響有待深入研究。Li等[5-6]在研究G?rtler渦二次失穩(wěn)時指出,G?rtler渦波長(展向)較大時二次失穩(wěn)由偶模式主導(dǎo),波長較小時更容易出現(xiàn)奇模式。文中工況1渦對之間的距離較遠,對應(yīng)的波長為3 mm,工況2中渦對之間的距離較近,波長為1.9 mm,文中兩個工況在失穩(wěn)模式和波長的對應(yīng)關(guān)系上與Li等對于G?rtler渦二次失穩(wěn)的結(jié)論相符合。

    3.4計算收斂性討論

    對于本文研究的渦流發(fā)生器強制轉(zhuǎn)捩流動,流向渦對結(jié)構(gòu)破碎的位置基本決定了轉(zhuǎn)捩位置,準確模擬渦破碎位置取決于對擾動增長過程的準確模擬,而準確模擬擾動增長的關(guān)鍵則是避免數(shù)值耗散對其產(chǎn)生明顯的抑制。因此,網(wǎng)格密度和數(shù)值方法所引入耗散的大小,是準確計算轉(zhuǎn)捩位置的關(guān)鍵。本節(jié)通過對數(shù)值耗散和網(wǎng)格密度的調(diào)節(jié),考察了本文所采用的數(shù)值方法的收斂性。

    圖8和圖9給出了在G1和G2兩套網(wǎng)格上計算結(jié)果隨網(wǎng)格和耗散參數(shù)κ的變化,其中σp為壓力均方根。本文以流向渦對破碎位置作為計算收斂的判斷依據(jù)。渦破碎位置不僅受網(wǎng)格的影響,同樣受到數(shù)值格式耗散的影響。大體趨勢是網(wǎng)格越密、數(shù)值格式耗散越低,渦破碎位置越靠前,這是由于數(shù)值耗散能夠抑制擾動的增長。隨著格式耗散的降低和網(wǎng)格的加密,渦對破碎位置前移并趨于不變。如前所述,本文采用的粗網(wǎng)格G1滿足LES模擬的需求[35],如果一定的格式耗散(取決于κ值)下,渦破碎位置不再隨網(wǎng)格加密變化,則可以認為數(shù)值方法是合適的。

    從圖8和圖9中流向渦對破碎位置以及流場中脈動壓力的增長看,在格式耗散較小(κ較大)的情況下,計算更容易達到網(wǎng)格收斂;κ=0.9時,G1和G2兩套網(wǎng)格下計算得到的渦破碎位置幾乎一致,可認為實現(xiàn)了網(wǎng)格收斂;對于文中兩個工況,在G1網(wǎng)格下,采用κ=0.9計算基本得到了收斂的流向渦對破碎位置。

    圖8 第二壓縮面上流場Q等值面圖 Fig.8 Q iso-surfaces of flow on the second compress surface

    圖9 沿流向渦壓力脈動的均方根Fig.9 Root-mean-square of pressure fluctuations along a stream-wise vortex

    4 結(jié) 論

    1) 渦流發(fā)生器強制轉(zhuǎn)捩過程的物理圖像可以描述為:渦流發(fā)生器在其下游產(chǎn)生一對流向渦,這對流向渦失穩(wěn)產(chǎn)生一系列“發(fā)卡”渦并逐漸演化直至破碎,最終使邊界層轉(zhuǎn)捩為湍流。

    2) 渦流發(fā)生器下游流向渦的失穩(wěn)包括兩種模式:偶模式和奇模式,前者產(chǎn)生對稱形式的渦破碎,后者產(chǎn)生彎曲形式的渦對破碎。本文開展的數(shù)值模擬研究同時捕捉到了這兩種失穩(wěn)模式。

    3) 采用隱式大渦模擬方法計算強制轉(zhuǎn)捩問題時,數(shù)值耗散對計算結(jié)果有明顯影響,文中對網(wǎng)格和數(shù)值耗散影響作了測試:在較小的格式耗散情況下,計算更容易達到網(wǎng)格收斂;隨著網(wǎng)格加密、耗散變小,計算得到的流向渦對破碎前移并趨于不變。

    [1]BERRY S A, DARYABEIGI K, WURSTER K. Boundary-layer transition on X-43A[J]. Journal of Spacecraft and Rockets, 2010, 47(6): 922-934.

    [2]KLEBANOFF P S, TIDSTROM K D. Mechanism by which a two-dimensional roughness element induces boundary-layer transition[J]. Physics of Fluids, 1972, 15(7): 1173-1188.

    [3]RESHOTKO E, TUMIN A. Role of transient growth in roughness-induced transition[J]. AIAA Journal, 2004, 42(4): 766-770.

    [4]SESCU A, PENDYALAY R K, THOMPSON D. On the growth of G?rtler vortices excited by distributed roughness elements: AIAA-2014-2885[R]. Reston: AIAA, 2014.

    [5]LI F, MALIK M R. Fundamental and subharmonic secondary instabilities of G?rtler vortices[J]. Journal of Fluid Mechanics, 1995, 297: 77-100.

    [6]LI F, CHOUDHARI M, CHANG C L, et al. Development and breakdown of G?rtler vortices in high speed boundary layers: AIAA-2010-705[R]. Reston: AIAA, 2010.

    [7]TULLIO N D, PAREDES P, SANDHAM N D, et al. Laminar-turbulent transition induced by a discrete roughness element in a supersonic boundary layer[J]. Journal of Fluid Mechanics, 2013, 735: 613-646.

    [8]BERRY S A, HORVATH T J. Discrete-roughness transition for hypersonic flight vehicles[J]. Journal of Spacecraft and Rockets, 2008, 45(2): 216-227.

    [9]BERRY S A, AUSLENDER A H, DILLEY A D, et al. Hypersonic boundary layer trip development for Hyper-X[J]. Journal of Spacecraft and Rockets, 2001, 38(6): 853-864.

    [10]BERRY S A, NOWAK R J, HORVATH T J. Boundary layer control for hypersonic airbreathing vehicles: AIAA-2004-2246[R]. Reston: AIAA, 2004.

    [11]趙慧勇, 周瑜, 倪鴻禮, 等. 高超聲速進氣道邊界層強迫轉(zhuǎn)捩試驗[J]. 實驗流體力學, 2012, 26(1): 1-6.

    ZHAO H Y, ZHOU Y, NI H L, et al. Test of forced boundary-layer transition on hypersonic inlet[J]. Journal of Experiments in Fluid Mechanics, 2012, 26(1): 1-6 (in Chinese).

    [12]DANEHY P M, BATHEL B, IVEY C, et al. NO PLIF study of hypersonic transition over a discrete hemispherical roughness element: AIAA-2009-0394[R]. Reston: AIAA, 2009.

    [13]岡敦殿, 易仕和, 陳植, 等. 超聲速平板上直立圓臺突起物繞流流場研究[C]//中國力學大會-2013. 北京: 中國力學學會, 2013.

    GANG D D, YI S H, CHEN Z, et al. Flow field investigations of supersonic flow over a circular protuberance mounted on a flat plate[C]//CCTAM-2013. Beijing: The Chinese Society of Theoretical and Applied Mechanics, 2013 (in Chinese).

    [14]張慶虎. 超聲速流動分離及其控制的試驗研究[D]. 長沙: 國防科學技術(shù)大學, 2013.

    ZHANG Q H. Experimental investigation of supersonic flow separation and its micro-ramp control[D]. Changsha: National University of Defense Technology, 2013 (in Chinese).

    [15]BORG M P, SCHNEIDER S P. Effect of freestream noise on roughness-induced transition for the X-51A forebody[J]. Journal of Spacecraft and Rockets, 2008, 45(6): 1106-1116.

    [16]鄧小兵. 來流噪聲對高超聲速進氣道強制轉(zhuǎn)捩的影響[C]//中國力學大會-2013. 北京: 中國力學學會, 2013.

    DENG X B. Freestream noise impact on forced transition of hypersonic boundary layer[C]//CCTAM-2013. Beijing: The Chinese Society of Theoretical and Applied Mechanics, 2013 (in Chinese).

    [17]MARXEN O, IACCARINO G, SHAQFEH E S G. Numerical simulations of hypersonic boundary-layer instability with localized roughness: AIAA-2011-0567[R]. Reston: AIAA, 2011.

    [18]IYER P S, MUPPIDI S, MAHESH K. Roughness-induced transition in high speed flows: AIAA-2011-0566[R]. Reston: AIAA, 2011.

    [19]DUAN Z W, XIAO Z X, FU S. Direct numerical simulation of hypersonic transition induced by ramp roughness elements: AIAA-2014-2037[R]. Reston: AIAA, 2014.

    [20]DUAN Z W, XIAO Z X. Direct numerical simulation of geometrical parameter effects on the hypersonic ramp-induced transition: AIAA-2014-2495[R]. Reston: AIAA, 2014.

    [21]涂國華, 燕振國, 趙曉慧, 等. SA和SST湍流模型對高超聲速邊界層強制轉(zhuǎn)捩的適應(yīng)性[J]. 航空學報, 2015, 36(5): 1471-1479.

    TU G H, YAN Z G, ZHAO X H, et al. SA and SST turbulence models for hypersonic forced boundary layer transition[J]. Acta Aeronautica et Astronautia Sinica, 2015, 36(5): 1471-1479 (in Chinese).

    [22]周玲, 閆超, 孔維萱. 高超聲速飛行器前體邊界層強制轉(zhuǎn)捩數(shù)值模擬[J]. 航空學報, 2014, 35(6): 1487-1495.

    ZHOU L, YAN C, KONG W X. Numerical simulation of forced boundary layer transition on hypersonic vehicle forebody[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(6): 1487-1495 (in Chinese).

    [23]周玲, 閆超, 郝子輝, 等. 轉(zhuǎn)捩模式與轉(zhuǎn)捩準則預(yù)測高超聲速邊界層流動[J]. 航空學報, 2016, 37(4): 1092-1102.

    ZHOU L, YAN C, HAO Z H, et al. Transition model and transition criteria for hypersonic boundary layer flow[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(4): 1092-1102 (in Chinese).

    [24]SAYADI T, MOIN P. Large eddy simulation of controlled transition to turbulence[J]. Physics of Fluids, 2012, 24(11): 923-938.

    [25]ORLIK E, FEDIOUN I. Hypersonic boundary-layer transition forced by wall injection: A numerical study[J]. Journal of Spacecraft and Rockets, 2014, 51(4): 1306-1318.

    [26]鄧小兵, 金玲. 讓誤差飛一會兒——湍流大渦模擬中的動態(tài)自適應(yīng)迎風方法[C]//中國力學大會-2013. 北京: 中國力學學會, 2013.

    DENG X B, JIN L. Let the error fly for a while—A dynamic adaptive upwind method for large eddy simulations of turbulent flow[C]//CCTAM-2013. Beijing: The Chinese Society of Theoretical and Applied Mechanics, 2013 (in Chinese).

    [27]朱自強, 吳子牛, 李津, 等. 應(yīng)用計算流體力學[M]. 北京: 北京航空航天大學出版社, 1998: 54-88.

    ZHU Z Q, WU Z N, LI J, et al. Application of computational fluid dynamics[M]. Beijing: Beihang University Press, 1998: 54-88 (in Chinese).

    [28]JAMESON A. Time dependent calculations using multigrid with application to unsteady flows past airfoils and wings: AIAA-1991-1596[R]. Reston: AIAA, 1991.

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

    [30]DUAN L, BEEKMAN I, MARTN M P. Direct numerical simulation of hypersonic turbulent boundary layers with varying freestream Mach number: AIAA-2010-0353[R]. Reston: AIAA, 2010.

    [31]LUND T S, WU X, SQUIRES K D. Generation of turbulent inflow data for spatially-developing boundary layer simulations[J]. Journal of Computational Physics, 1998, 140(2): 233-258.

    [32]EDWARDS J R, CHOI J, BOLES J A. Hybrid LES/RANS simulation of a mach 5 compression-corner interaction: AIAA-2008-0718[R]. Reston: AIAA, 2008.

    [33]SHUR M L, SPALART P R, STRELETS M K, et al. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities[J]. International Journal of Heat and Fluid Flow, 2008, 29(6): 1638-1649

    [34]SPALART P R, ALLMARAS S R. A one-equation turbulence model for aerodynamic flows: AIAA-1992-0439[R]. Reston: AIAA, 1992.

    [35]GEORGIADIS N J, RIZZETTA D P, FUREBY C. Large-eddy simulation: current capabilities, recommended practices, and future research[J]. AIAA Journal, 2010, 48(8): 1772-1784.

    趙曉慧男, 博士研究生, 助理研究員。主要研究方向: 湍流數(shù)值模擬, 高超聲速空氣動力學。

    Tel.: 0816-2463219

    E-mail: xhzhao@skla.cardc.cn

    鄧小兵男, 博士, 副研究員。主要研究方向: 計算流體力學, 湍流數(shù)值模擬。

    Tel.: 0816-2463219

    E-mail: xbdeng@skla.cardc.cn

    毛枚良男, 博士, 研究員, 博士生導(dǎo)師。主要研究方向: 計算流體力學, 高階精度格式, 高超聲速空氣動力學。

    Tel.: 0816-2463296

    E-mail: mm1219@163.com

    Large eddy simulation for forced transition flow at hypersonic inlet

    ZHAO Xiaohui1, DENG Xiaobing1,*, MAO Meiliang2, YANG Wei1, ZHAO Huiyong3

    1. State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center,Mianyang621000, China 2. Computational Aerodynamics Institute, China Aerodynamics Research and Development Center,Mianyang621000, China 3. Science and Technology on Scramjet Laboratory, China Aerodynamics Research and Development Center,Mianyang 621000, China

    Forced transition is commonly used for a hypersonic air breathing vehicle to ensure its scramjet’s normal start by placing roughness elements or vortex generators in the inlet boundary layer. To get a clear image of forced transition process, an implicit large eddy simulation method is used for laminar-turbulent forced transition flows in the boundary layer of a hypersonic inlet configuration. Main structures are obtained for the forced transition flows induced by the vortex generators shaped of diamonds and ramps, which show that the flows are dominated by stream-wise counter-rotating vortices. Fluctuations grow in the stream-wise direction and cause instabilities. Two fundamental modes of the instabilities are found in the simulations, the even mode and the odd one. The even mode leads to a varicose way of breaking down of the counter-rotating vortices, while the odd one leads to a sinuous way. Strings of hairpin vortices are generated after the breaking down and finally cause the transition. A discussion is carried out on the computation convergence together with the influence of the grids and numerical dissipation.

    hypersonic; inlet; forced transition; vortex generator; large eddy simulation

    2016-02-16; Revised: 2016-02-17; Accepted: 2016-06-21; Published online: 2016-06-2715:34

    National Natural Science Foundation of China (11402289)

    . Tel.: 0816-2463219E-mail: xbdeng@skla.cardc.cn

    2016-02-16; 退修日期: 2016-02-17; 錄用日期: 2016-06-21;

    時間: 2016-06-2715:34

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

    國家自然科學基金 (11402289)

    .Tel.: 0816-2463219E-mail: xbdeng@skla.cardc.cn

    10.7527/S1000-6893.2016.0200

    V211.3; O355

    A

    1000-6893(2016)08-2445-09

    引用格式: 趙曉慧, 鄧小兵, 毛枚良, 等. 高超聲速進氣道強制轉(zhuǎn)捩流動的大渦模擬[J]. 航空學報, 2016, 37(8): 2445-2453. ZHAO X H, DENG X B, MAO M L, et al. Large eddy simulation for forced transition flow at hypersonic inlet[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2445-2453.

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

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

    猜你喜歡
    大渦進氣道邊界層
    基于AVL-Fire的某1.5L發(fā)動機進氣道優(yōu)化設(shè)計
    基于輔助進氣門的進氣道/發(fā)動機一體化控制
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    基于壁面射流的下?lián)舯┝鞣欠€(wěn)態(tài)風場大渦模擬
    軸流風機葉尖泄漏流動的大渦模擬
    基于大渦模擬的旋風分離器錐體結(jié)構(gòu)影響研究
    一類具有邊界層性質(zhì)的二次奇攝動邊值問題
    The coupling characteristics of supersonic dual inlets for missile①
    非特征邊界的MHD方程的邊界層
    某柴油機進氣道數(shù)值模擬及試驗研究
    汽車零部件(2014年2期)2014-03-11 17:46:30
    男女午夜视频在线观看| 最黄视频免费看| 国产麻豆69| 天堂8中文在线网| 在线 av 中文字幕| 国产真人三级小视频在线观看| 国产午夜精品久久久久久| 亚洲精品自拍成人| 久久热在线av| 一区二区三区激情视频| 男女床上黄色一级片免费看| 两性夫妻黄色片| 在线播放国产精品三级| 国产精品免费视频内射| 宅男免费午夜| 最新的欧美精品一区二区| 亚洲avbb在线观看| 曰老女人黄片| 少妇被粗大的猛进出69影院| 中亚洲国语对白在线视频| 777米奇影视久久| 日韩免费高清中文字幕av| 91字幕亚洲| a级片在线免费高清观看视频| 久久久水蜜桃国产精品网| 美女高潮到喷水免费观看| 一个人免费在线观看的高清视频| 麻豆国产av国片精品| 欧美久久黑人一区二区| 每晚都被弄得嗷嗷叫到高潮| 性高湖久久久久久久久免费观看| 日本撒尿小便嘘嘘汇集6| 黄色毛片三级朝国网站| 乱人伦中国视频| 欧美日韩中文字幕国产精品一区二区三区 | kizo精华| 香蕉国产在线看| 欧美久久黑人一区二区| 欧美日韩亚洲高清精品| 午夜福利免费观看在线| 一区二区av电影网| 久久精品亚洲精品国产色婷小说| 久久国产精品男人的天堂亚洲| 国产亚洲精品久久久久5区| 成人国产av品久久久| 亚洲精品国产精品久久久不卡| 制服诱惑二区| 久久国产亚洲av麻豆专区| 每晚都被弄得嗷嗷叫到高潮| 久久国产精品大桥未久av| 精品国产一区二区久久| 久久精品国产综合久久久| 亚洲伊人色综图| 亚洲av电影在线进入| 亚洲国产av新网站| 欧美性长视频在线观看| 丁香欧美五月| av超薄肉色丝袜交足视频| 免费久久久久久久精品成人欧美视频| 侵犯人妻中文字幕一二三四区| 久热这里只有精品99| 一级,二级,三级黄色视频| 黑人猛操日本美女一级片| 亚洲欧美精品综合一区二区三区| 亚洲久久久国产精品| 国产色视频综合| 国产精品一区二区免费欧美| 中文字幕制服av| 国产亚洲欧美在线一区二区| 最新在线观看一区二区三区| 亚洲va日本ⅴa欧美va伊人久久| 丝袜喷水一区| av超薄肉色丝袜交足视频| 熟女少妇亚洲综合色aaa.| 国产精品熟女久久久久浪| 丁香欧美五月| 久久久久国内视频| 国产99久久九九免费精品| 国产一区二区三区视频了| 99精品欧美一区二区三区四区| 国产黄频视频在线观看| 脱女人内裤的视频| av天堂在线播放| 一个人免费在线观看的高清视频| 欧美一级毛片孕妇| 国产精品99久久99久久久不卡| 国产成人系列免费观看| cao死你这个sao货| 久久人妻av系列| 国产精品亚洲一级av第二区| 捣出白浆h1v1| 国产成人免费观看mmmm| 国产亚洲精品一区二区www | 夜夜爽天天搞| 精品一区二区三区av网在线观看 | 1024视频免费在线观看| 午夜免费鲁丝| 好男人电影高清在线观看| 日韩熟女老妇一区二区性免费视频| 久久香蕉激情| 香蕉久久夜色| 久久久精品免费免费高清| 亚洲av国产av综合av卡| 成人18禁高潮啪啪吃奶动态图| 热re99久久精品国产66热6| www.999成人在线观看| 99香蕉大伊视频| 成人精品一区二区免费| 亚洲国产毛片av蜜桃av| 中国美女看黄片| 久久久久精品国产欧美久久久| 在线观看免费视频日本深夜| 国产三级黄色录像| 十分钟在线观看高清视频www| 国产精品久久久久成人av| 国产日韩欧美在线精品| 在线观看免费日韩欧美大片| 丝袜美足系列| 欧美成人免费av一区二区三区 | 不卡一级毛片| 精品国产一区二区三区久久久樱花| 欧美黄色淫秽网站| 亚洲欧美日韩另类电影网站| 在线天堂中文资源库| 操美女的视频在线观看| 少妇 在线观看| 国产一卡二卡三卡精品| 亚洲伊人久久精品综合| 一级a爱视频在线免费观看| 亚洲男人天堂网一区| 精品国产乱码久久久久久小说| 国产精品电影一区二区三区 | 国内毛片毛片毛片毛片毛片| 亚洲欧洲精品一区二区精品久久久| 一级a爱视频在线免费观看| 久久精品熟女亚洲av麻豆精品| 国产不卡av网站在线观看| 水蜜桃什么品种好| 新久久久久国产一级毛片| 欧美黑人欧美精品刺激| 汤姆久久久久久久影院中文字幕| 国产一区二区在线观看av| 大片免费播放器 马上看| 欧美激情久久久久久爽电影 | 伊人久久大香线蕉亚洲五| 老汉色av国产亚洲站长工具| 麻豆av在线久日| 国产区一区二久久| 人人澡人人妻人| 一区二区三区乱码不卡18| 大片免费播放器 马上看| 日本wwww免费看| 国产在线精品亚洲第一网站| 欧美人与性动交α欧美软件| 美女视频免费永久观看网站| a级片在线免费高清观看视频| 交换朋友夫妻互换小说| 99久久国产精品久久久| 国产高清videossex| 国产无遮挡羞羞视频在线观看| 成人亚洲精品一区在线观看| av天堂久久9| 亚洲专区中文字幕在线| 18在线观看网站| 久久国产亚洲av麻豆专区| 亚洲国产毛片av蜜桃av| 美女主播在线视频| 久热爱精品视频在线9| 国产欧美亚洲国产| 高清在线国产一区| 久久久国产欧美日韩av| 男人操女人黄网站| 曰老女人黄片| 成在线人永久免费视频| 国产精品一区二区免费欧美| 午夜福利欧美成人| bbb黄色大片| 在线 av 中文字幕| 人人妻人人澡人人爽人人夜夜| 久久久国产一区二区| 国产无遮挡羞羞视频在线观看| 香蕉丝袜av| 久久久久精品国产欧美久久久| 男人舔女人的私密视频| 欧美日韩成人在线一区二区| 欧美国产精品va在线观看不卡| 亚洲第一青青草原| 午夜激情久久久久久久| 悠悠久久av| 国产亚洲精品久久久久5区| 搡老乐熟女国产| 亚洲av成人不卡在线观看播放网| a在线观看视频网站| 99久久99久久久精品蜜桃| 欧美黑人精品巨大| 久久久久久亚洲精品国产蜜桃av| 日本黄色视频三级网站网址 | 日韩大码丰满熟妇| 91九色精品人成在线观看| 他把我摸到了高潮在线观看 | 久热这里只有精品99| 国产欧美日韩综合在线一区二区| 又紧又爽又黄一区二区| 中文字幕另类日韩欧美亚洲嫩草| 丁香六月欧美| 另类亚洲欧美激情| 色婷婷av一区二区三区视频| 蜜桃国产av成人99| 最近最新免费中文字幕在线| 高清欧美精品videossex| 欧美老熟妇乱子伦牲交| 人妻久久中文字幕网| 婷婷成人精品国产| 亚洲av国产av综合av卡| av电影中文网址| 精品乱码久久久久久99久播| 美女扒开内裤让男人捅视频| 天天躁夜夜躁狠狠躁躁| 在线观看免费日韩欧美大片| 老司机在亚洲福利影院| 国产亚洲精品久久久久5区| 久久久久久久久免费视频了| 久久久精品免费免费高清| 日本wwww免费看| 精品少妇一区二区三区视频日本电影| 亚洲av日韩在线播放| 黄色毛片三级朝国网站| 国产伦理片在线播放av一区| 国产激情久久老熟女| 18禁观看日本| 亚洲熟女精品中文字幕| 一级a爱视频在线免费观看| 天天添夜夜摸| 50天的宝宝边吃奶边哭怎么回事| 欧美变态另类bdsm刘玥| 成人三级做爰电影| 日本五十路高清| 在线观看免费日韩欧美大片| 欧美日韩av久久| 90打野战视频偷拍视频| 国产av国产精品国产| 精品视频人人做人人爽| 在线av久久热| 精品亚洲乱码少妇综合久久| 日韩有码中文字幕| 另类精品久久| 人成视频在线观看免费观看| 日韩免费高清中文字幕av| 91精品三级在线观看| 欧美日韩亚洲国产一区二区在线观看 | 久久精品亚洲熟妇少妇任你| 国产男女超爽视频在线观看| 亚洲精品av麻豆狂野| 80岁老熟妇乱子伦牲交| 2018国产大陆天天弄谢| 色在线成人网| 免费不卡黄色视频| 日韩制服丝袜自拍偷拍| 在线观看免费午夜福利视频| 另类精品久久| 超色免费av| 在线十欧美十亚洲十日本专区| 女性生殖器流出的白浆| 看免费av毛片| 国产三级黄色录像| 欧美大码av| 搡老岳熟女国产| 久9热在线精品视频| 亚洲一区中文字幕在线| 久久人妻熟女aⅴ| 国产av国产精品国产| 欧美午夜高清在线| 中文字幕色久视频| 国产亚洲精品第一综合不卡| 国产成人一区二区三区免费视频网站| 老汉色∧v一级毛片| 母亲3免费完整高清在线观看| 丝袜美足系列| 日韩熟女老妇一区二区性免费视频| 精品第一国产精品| avwww免费| 国产97色在线日韩免费| 精品亚洲乱码少妇综合久久| 午夜老司机福利片| 久久人人97超碰香蕉20202| 成人手机av| 老熟妇乱子伦视频在线观看| 在线看a的网站| 亚洲午夜理论影院| 色94色欧美一区二区| 久久久久久久久免费视频了| a级片在线免费高清观看视频| 亚洲精品国产一区二区精华液| 国产日韩欧美亚洲二区| 日韩有码中文字幕| 狠狠精品人妻久久久久久综合| 十八禁人妻一区二区| 手机成人av网站| 中文字幕另类日韩欧美亚洲嫩草| 成人永久免费在线观看视频 | 不卡av一区二区三区| 精品久久久久久久毛片微露脸| 50天的宝宝边吃奶边哭怎么回事| 男女床上黄色一级片免费看| 天堂8中文在线网| 久久毛片免费看一区二区三区| 精品久久蜜臀av无| 亚洲欧美色中文字幕在线| 菩萨蛮人人尽说江南好唐韦庄| 悠悠久久av| 久久性视频一级片| 激情视频va一区二区三区| 母亲3免费完整高清在线观看| 亚洲va日本ⅴa欧美va伊人久久| 黑丝袜美女国产一区| 少妇猛男粗大的猛烈进出视频| 亚洲精品美女久久久久99蜜臀| 午夜成年电影在线免费观看| 久久亚洲精品不卡| 国产精品欧美亚洲77777| 午夜视频精品福利| 亚洲一码二码三码区别大吗| 欧美国产精品va在线观看不卡| 国产欧美日韩精品亚洲av| 日本撒尿小便嘘嘘汇集6| 蜜桃国产av成人99| 99香蕉大伊视频| 国产精品久久久久久精品古装| 99riav亚洲国产免费| 亚洲精品一卡2卡三卡4卡5卡| 久久久久久久精品吃奶| 激情在线观看视频在线高清 | 大片电影免费在线观看免费| 久久精品熟女亚洲av麻豆精品| 中文字幕色久视频| 精品国产一区二区三区四区第35| 国产亚洲欧美精品永久| 久久狼人影院| 欧美国产精品一级二级三级| h视频一区二区三区| 亚洲综合色网址| 亚洲美女黄片视频| 亚洲第一欧美日韩一区二区三区 | 日韩有码中文字幕| 伊人久久大香线蕉亚洲五| 99热国产这里只有精品6| 欧美亚洲日本最大视频资源| 中文字幕人妻熟女乱码| 亚洲全国av大片| 69精品国产乱码久久久| 狠狠精品人妻久久久久久综合| 欧美人与性动交α欧美精品济南到| 激情视频va一区二区三区| 亚洲人成电影观看| 999精品在线视频| 欧美日韩视频精品一区| 欧美亚洲日本最大视频资源| 在线观看免费午夜福利视频| av线在线观看网站| 最新在线观看一区二区三区| 日本欧美视频一区| 多毛熟女@视频| 人妻 亚洲 视频| svipshipincom国产片| 精品一品国产午夜福利视频| 伦理电影免费视频| 亚洲av美国av| 怎么达到女性高潮| 性高湖久久久久久久久免费观看| 欧美日韩国产mv在线观看视频| 午夜福利在线免费观看网站| 又黄又粗又硬又大视频| 交换朋友夫妻互换小说| 亚洲五月色婷婷综合| 久久免费观看电影| 激情视频va一区二区三区| 国产色视频综合| 亚洲精品国产区一区二| 日韩人妻精品一区2区三区| 99九九在线精品视频| 涩涩av久久男人的天堂| 亚洲国产欧美日韩在线播放| 精品久久久久久久毛片微露脸| 久久久国产一区二区| 亚洲全国av大片| 电影成人av| 亚洲专区字幕在线| 露出奶头的视频| 国产成人av激情在线播放| 人人妻,人人澡人人爽秒播| 亚洲一区二区三区欧美精品| 看免费av毛片| 一区福利在线观看| 性少妇av在线| 久久久精品区二区三区| 国产伦人伦偷精品视频| 国产成人av激情在线播放| 精品国产一区二区三区四区第35| 国产成人av教育| 在线观看www视频免费| 国产av精品麻豆| 精品久久蜜臀av无| 亚洲人成伊人成综合网2020| 又紧又爽又黄一区二区| 宅男免费午夜| 少妇的丰满在线观看| 男人舔女人的私密视频| 久久九九热精品免费| 侵犯人妻中文字幕一二三四区| 亚洲国产欧美一区二区综合| 757午夜福利合集在线观看| 国产午夜精品久久久久久| 精品欧美一区二区三区在线| 99国产极品粉嫩在线观看| 香蕉国产在线看| 欧美日韩av久久| 国产亚洲精品久久久久5区| 亚洲自偷自拍图片 自拍| 视频区欧美日本亚洲| 熟女少妇亚洲综合色aaa.| 日本五十路高清| 99精品久久久久人妻精品| 少妇粗大呻吟视频| 久9热在线精品视频| 亚洲中文字幕日韩| 国内毛片毛片毛片毛片毛片| 国产高清国产精品国产三级| 丁香六月天网| 国产av又大| 欧美日韩精品网址| 久久国产精品大桥未久av| 色视频在线一区二区三区| 久久国产精品男人的天堂亚洲| 午夜老司机福利片| 国产一区二区三区视频了| 中亚洲国语对白在线视频| xxxhd国产人妻xxx| 一进一出好大好爽视频| 亚洲综合色网址| 日本五十路高清| 中文字幕人妻熟女乱码| 亚洲欧美激情在线| 交换朋友夫妻互换小说| 国产av一区二区精品久久| 国产片内射在线| 99国产精品免费福利视频| 一区二区三区国产精品乱码| 久久精品国产亚洲av高清一级| 一二三四在线观看免费中文在| 91国产中文字幕| 自线自在国产av| 一本—道久久a久久精品蜜桃钙片| 精品一品国产午夜福利视频| 精品国产乱子伦一区二区三区| 成人亚洲精品一区在线观看| 国产精品国产高清国产av | cao死你这个sao货| 亚洲色图av天堂| 国产色视频综合| 日韩制服丝袜自拍偷拍| 高清av免费在线| 欧美午夜高清在线| 国产成人免费无遮挡视频| 国产亚洲精品久久久久5区| 国产精品1区2区在线观看. | 老司机午夜福利在线观看视频 | 自拍欧美九色日韩亚洲蝌蚪91| 国产不卡av网站在线观看| 日日爽夜夜爽网站| 性色av乱码一区二区三区2| 国产精品1区2区在线观看. | 熟女少妇亚洲综合色aaa.| 亚洲,欧美精品.| 91精品国产国语对白视频| 91av网站免费观看| 亚洲欧美激情在线| 极品人妻少妇av视频| 国产真人三级小视频在线观看| 久久影院123| 国产一区二区激情短视频| 中文字幕精品免费在线观看视频| 高清毛片免费观看视频网站 | 亚洲午夜理论影院| 午夜福利在线观看吧| 悠悠久久av| 女人被狂操c到高潮| 一进一出好大好爽视频| 中文字幕av在线有码专区| 成人午夜高清在线视频| 人妻夜夜爽99麻豆av| 欧美3d第一页| 亚洲,欧美精品.| netflix在线观看网站| 啪啪无遮挡十八禁网站| 国产高潮美女av| 成年女人毛片免费观看观看9| 午夜免费观看网址| netflix在线观看网站| a级毛片在线看网站| 老汉色∧v一级毛片| 女人被狂操c到高潮| 国产精品一及| 欧美日本视频| 亚洲精品久久国产高清桃花| 欧美又色又爽又黄视频| 久久久久精品国产欧美久久久| 国产亚洲欧美在线一区二区| 精品福利观看| 国产麻豆成人av免费视频| 在线观看一区二区三区| 国产欧美日韩精品一区二区| 91久久精品国产一区二区成人 | 18禁黄网站禁片免费观看直播| 国产aⅴ精品一区二区三区波| 久久久久久久午夜电影| 国产乱人伦免费视频| 亚洲国产色片| 色噜噜av男人的天堂激情| 精品国内亚洲2022精品成人| 亚洲七黄色美女视频| 日本黄大片高清| 少妇的丰满在线观看| 亚洲精品色激情综合| 少妇的逼水好多| 免费人成视频x8x8入口观看| 999久久久国产精品视频| 熟女少妇亚洲综合色aaa.| 变态另类丝袜制服| 老熟妇乱子伦视频在线观看| 操出白浆在线播放| 一区福利在线观看| 一区二区三区国产精品乱码| 欧美成人一区二区免费高清观看 | 日本在线视频免费播放| a级毛片a级免费在线| 变态另类成人亚洲欧美熟女| 免费在线观看成人毛片| 18禁国产床啪视频网站| 亚洲专区字幕在线| 精品欧美国产一区二区三| 美女免费视频网站| 欧美最黄视频在线播放免费| 欧美zozozo另类| 国内少妇人妻偷人精品xxx网站 | 国产一区二区在线av高清观看| 国产精品国产高清国产av| 日韩有码中文字幕| 久久婷婷人人爽人人干人人爱| 看黄色毛片网站| 这个男人来自地球电影免费观看| 最好的美女福利视频网| 99热这里只有是精品50| 狠狠狠狠99中文字幕| 午夜视频精品福利| 亚洲国产精品久久男人天堂| 国产一区二区三区在线臀色熟女| 熟女电影av网| 成年女人毛片免费观看观看9| 长腿黑丝高跟| 国语自产精品视频在线第100页| 美女 人体艺术 gogo| 九九久久精品国产亚洲av麻豆 | 变态另类丝袜制服| 国语自产精品视频在线第100页| 99热6这里只有精品| 国产精品一及| 两性午夜刺激爽爽歪歪视频在线观看| 午夜精品在线福利| 日韩 欧美 亚洲 中文字幕| 国产欧美日韩精品一区二区| 色精品久久人妻99蜜桃| 国产男靠女视频免费网站| 久久久久久久久免费视频了| 岛国在线免费视频观看| 亚洲专区中文字幕在线| 午夜久久久久精精品| 观看美女的网站| 午夜视频精品福利| www.熟女人妻精品国产| av在线天堂中文字幕| 国产99白浆流出| 亚洲avbb在线观看| 天天添夜夜摸| 色哟哟哟哟哟哟| 日本三级黄在线观看| 俄罗斯特黄特色一大片| 我的老师免费观看完整版| www.精华液| 欧美性猛交╳xxx乱大交人| 不卡一级毛片| 法律面前人人平等表现在哪些方面| 舔av片在线| 色哟哟哟哟哟哟| svipshipincom国产片| 国内毛片毛片毛片毛片毛片| 每晚都被弄得嗷嗷叫到高潮| www.熟女人妻精品国产| 男女之事视频高清在线观看| 欧美一区二区精品小视频在线| 精品久久久久久久久久久久久| 日本免费a在线| 18禁观看日本| xxx96com| 啪啪无遮挡十八禁网站| 国产三级在线视频| 成人欧美大片| 嫁个100分男人电影在线观看| 国产一区二区在线观看日韩 | 亚洲专区国产一区二区| 国产成人欧美在线观看| 国产一区二区激情短视频| 99视频精品全部免费 在线 | 一二三四社区在线视频社区8| 国产精品99久久99久久久不卡| 欧美3d第一页|