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

    激波/湍流邊界層干擾壓力脈動特性數(shù)值研究1)

    2021-11-09 06:26:08童福林段俊亦周桂宇李新亮
    力學(xué)學(xué)報 2021年7期
    關(guān)鍵詞:物面外層邊界層

    童福林 段俊亦 周桂宇 李新亮

    * (中國空氣動力研究與發(fā)展中心空氣動力學(xué)國家重點實驗室,四川綿陽 621000)

    ? (中國科學(xué)院力學(xué)研究所高溫氣體動力學(xué)國家重點實驗室,北京 100190)

    ** (中國空氣動力研究與發(fā)展中心計算空氣動力研究所,四川綿陽 621000)

    ?? (中國科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

    引言

    激波與湍流邊界層的相互作用問題廣泛存在于各類高速飛行器中.已有研究表明[1-2],在強激波作用下,干擾區(qū)內(nèi)壓力脈動峰值急劇升高,分離激波的非定常運動導(dǎo)致物面脈動壓力低頻能量的急劇增強,這將使得飛行器結(jié)構(gòu)出現(xiàn)振動疲勞問題,進而嚴重影響飛行安全.因此,深入研究激波/湍流邊界層干擾區(qū)內(nèi)壓力脈動特性對認識和理解激波與湍流邊界層的相互作用機制非常重要.

    自20 世紀70 年代以來,國內(nèi)外學(xué)者一直致力于對該問題的研究.早期研究主要以風(fēng)洞試驗為主,受限于實驗測量技術(shù),風(fēng)洞試驗研究多局限于物面壓力脈動信號,并以此獲得干擾區(qū)內(nèi)典型流動結(jié)構(gòu)的特征信息,如分離激波非定常運動特性、湍流結(jié)構(gòu)、速度場等.Settle 等[3]研究了激波強度對壓縮拐角物面壓力脈動強度的影響規(guī)律.Ringuette 等[4]研究表明,低雷諾數(shù)下,干擾區(qū)內(nèi)物面壓力脈動強度降低而間歇性增強.在分離激波低頻振蕩方面,Erengil等[5]采用條件取樣和變窗口系綜平均技術(shù)分析了物面壓力分布與分離激波非定常運動的相關(guān)性.結(jié)果表明,分離泡內(nèi)的低頻壓力脈動與激波位置密切相關(guān):隨著激波往下游移動,分離泡內(nèi)壓力升高,反之則降低.Dolling 等[6]比較了不同構(gòu)型誘導(dǎo)的激波/湍流邊界層干擾區(qū)物面壓力脈動特性,研究發(fā)現(xiàn),分離激波的低頻振蕩機制主要來源于下游流場的局部或全局脈動,其過零頻率與上游湍流邊界層的特征頻率無關(guān).此外,Selig 等[7]基于物面壓力與邊界層質(zhì)量通量的時空關(guān)聯(lián)分析,給出了干擾區(qū)內(nèi)湍流結(jié)構(gòu)的演化規(guī)律.Piponniau 等[8]采用本征正交分解和線性隨機估計方法,通過實驗測得的物面壓力脈動信息對速度場進行了重構(gòu).

    由于風(fēng)洞試驗?zāi)軌虻玫降膲毫π畔O為有限,為此國內(nèi)外學(xué)者針對干擾區(qū)內(nèi)壓力脈動特性也開展了大量的高精度數(shù)值模擬研究.特別是近些年,直接數(shù)值模擬方法(direct numerical simulation,DNS)已經(jīng)在激波/湍流邊界層干擾問題方面取得了長足進展,其中Pirozzoli 等[9-11]的工作最具代表性.Pirozzoli等[9]分析比較了干擾區(qū)內(nèi)物面壓力脈動與剪切層壓力脈動的差異,并基于聲共振機制成功解釋了分離激波的非定常運動現(xiàn)象.Bernardini 等[10]著重探討了干擾區(qū)內(nèi)物面壓力脈動特性的演化規(guī)律,如脈動強度、時空關(guān)聯(lián)、對流速度和頻譜等,發(fā)現(xiàn)物面壓力脈動對流速度明顯降低而其展向尺度則顯著增強.隨后,Volpiani 等[11]進一步考察了壁面溫度對物面壓力脈動特性的影響規(guī)律.與此同時,國內(nèi)童福林等針對物面壓力脈動特性也開展了較為詳細的直接數(shù)值模擬研究,歸納了激波強度[12]、壁面溫度[13]、膨脹效應(yīng)[14-15]等因素的影響規(guī)律.

    盡管在激波/湍流邊界層干擾區(qū)壓力脈動特性研究方面,國內(nèi)外已開展了大量的風(fēng)洞試驗和數(shù)值研究工作,但是這些研究大多關(guān)注物面壓力脈動特性及其與分離激波非定常運動的相互關(guān)系方面,而關(guān)于邊界層內(nèi)壓力脈動在干擾區(qū)內(nèi)演化機制的認識仍然是很有限的.實際上,Duan 等[16]的研究結(jié)果表明,邊界層內(nèi)壓力脈動結(jié)構(gòu)特征與物面壓力脈動完全不同,其特征頻率沿法向?qū)⒔档投卣鞒叨葘⒓眲≡龃?然而,目前學(xué)術(shù)界對邊界層內(nèi)外層大尺度壓力脈動結(jié)構(gòu)與分離激波的相互作用機制仍然缺乏充分的認識,亟待開展相關(guān)研究工作.

    本文采用直接數(shù)值模擬方法研究M∞=2.25,Reθ=3567 來流下激波角β=33.2°的入射激波/平板湍流邊界層干擾區(qū)內(nèi)壓力脈動特性.為了便于比較和驗證結(jié)果,計算參數(shù)的選取與Dupont 等[17]的實驗和Fang 等[18]的DNS 相近.計算結(jié)果與前人的DNS 結(jié)果及風(fēng)洞實驗數(shù)據(jù)進行了仔細驗證,著重探討干擾區(qū)外層壓力脈動與物面壓力脈動演化機制的差異,特別是激波干擾對脈動強度、功率譜、特征尺度及時空關(guān)聯(lián)特性等的影響規(guī)律.

    1 DNS 計算

    控制方程為曲線坐標系(τ,ξ,η,ζ)下的三維可壓縮無量綱Navier-Stokes 方程組

    其中J?1為直角坐標系(x,y,z)變換為曲線坐標系(ξ,η,ζ) 的Jacobian 矩陣,U為守恒變量,F,G和H分別為(x,y,z)方向上對應(yīng)的無黏通量,Fv,Gv和Hv分別為(x,y,z)方向上對應(yīng)的黏性通量,具體表達式參見文獻[19].本文采用的解算器是高精度有限差分軟件OPENCFD-SC,該求解器目前已在壓縮拐角[19-20]、膨脹角入射激波/湍流邊界層干擾[14-15,21]等復(fù)雜流動問題中得到廣泛應(yīng)用和驗證確認,可以保證計算結(jié)果的可靠性.

    對流項計算時,采用Martin 等[22]優(yōu)化構(gòu)造的WENO_SYMBO_LMT 格式以及Steger-Warming流通量分裂方法,這樣能夠在極大抑制激波間斷區(qū)的強數(shù)值振蕩同時也保證了對邊界層內(nèi)不同尺度湍流脈動結(jié)構(gòu)的高分辨率捕捉.黏性項采用八階中心差分格式進行離散,時間方向上推進采用三階Runge-Kutta 方法.

    圖1 為計算模型示意圖,來流方向為圖中從左往右,來流參數(shù)如下:馬赫數(shù)M∞=2.25,基于單位長度來流雷諾數(shù)Re∞=2.5×104mm?1,靜溫T∞=169.44 K.模型長度Lx×Ly×Lz=137.6 mm×12.7 mm×4.4 mm,這里L(fēng)x,Ly和Lz分別為流向、法向和展向長度,坐標系原點取為計算域入口處.采用與文獻[21]類似的方法生成入射激波,在計算域上邊界xin=82.2 mm左右兩側(cè)分別設(shè)置為波前和波后參數(shù),即:在xin左側(cè)流動參數(shù)取為自由來流參數(shù),而在xin右側(cè)則按照激波關(guān)系式Rankine-Hugoniot 給出.這里激波角β取為33.2°,圖中入射激波在壁面上的名義入射點為xs=101.6 mm.

    圖1 計算模型Fig.1 Computational model

    計算網(wǎng)格點數(shù)為Nx×Ny×Nz=3700×300×250 (流向×法向×展向).圖2 為計算網(wǎng)格示意圖,為了便于展示,流向間隔了10 個點,法向間隔了5 個點.網(wǎng)格生成采用代數(shù)解析方法,這里流向網(wǎng)格點采用三段分布:上游轉(zhuǎn)捩區(qū)0 mm

    圖2 計算網(wǎng)格示意圖Fig.2 Sketch of the computational grid

    邊界條件設(shè)置如下:計算域入口處取為相同來流下的層流解(laminar inflow),具體分布如圖3 所示.壁面取為無滑移和等溫壁,壁溫Tw=321.9 K,上邊界采用簡單無反射邊界條件.出口處采用超聲速出口無反射條件及緩沖區(qū),以便消除擾動波反射帶來的影響.計算域展向兩側(cè)采用周期性邊界條件.

    圖3 入口處層流剖面Fig.3 Laminar profile at inlet

    計算時,通過在上游壁面添加多頻吹吸擾動從而轉(zhuǎn)捩生成充分發(fā)展湍流邊界層,進而與下游入射激波產(chǎn)生相互作用.這里,法向擾動速度Vbs的具體表達式如下[18]

    DNS 計算時取lmax=10,mmax=5,擾動幅值A(chǔ)=0.2,擾動頻率?=0.628U∞/δ,這里δ表示參考點處邊界層厚度,下文類似.擾動帶起始點分別為xa=7.62 mm 和xb=20.32 mm.式中相位差φl和φm取為0 到1 的隨機數(shù).

    表1 給出了參考點xref處湍流邊界層參數(shù)與Fang 等[18]DNS 結(jié)果比較情況,表中Reδ,Reδ*,Reθ分別為基于邊界層厚度、位移厚度和動量厚度的雷諾數(shù),Cf為摩阻系數(shù).兩者較為接近,這也表明本文計算結(jié)果的準確性.為了進一步驗證本文DNS 結(jié)果的可靠性,還與以往數(shù)值模擬結(jié)果和風(fēng)洞實驗數(shù)據(jù)進行了比較,包括湍流邊界層的速度剖面、脈動強度,物面壓力頻譜以及平均壓力分布.

    表1 參考點xref 處湍流邊界層參數(shù)Table 1 Turbulent boundary layer parameters at xref

    首先,圖4~ 圖6 分別給出了上游湍流邊界層在參考點xref處的統(tǒng)計特性.圖4 為平均流向速度剖面,定義為

    圖4 平均流向速度剖面Fig.4 Mean streamwise velocity profile at xref

    式中ρ和u分別為密度和流向速度,上標 ? 表示平均量.需要特別說明的是,本文中平均指的是時間和展向平均.速度分布在黏性底層y+<10 符合線性律,在30

    圖5 湍流脈動強度Fig.5 Turbulence intensities at xref

    圖6 物面壓力功率譜Fig.6 Power spectral density for wall pressure at xref

    此外,圖7 還給出了干擾區(qū)內(nèi)物面壓力分布情況,其中物面壓力P*和流向坐標x*分別定義為

    圖7 物面壓力分布Fig.7 Distribution of wall pressure

    式中,P1是反射激波的波后壓力值,xr*為轉(zhuǎn)換后的反射激波腳對應(yīng)的平均流向位置,這里采用與Fang等[18]類似的方法,即(P1+P∞)/2 對應(yīng)的流向位置.圖中紅色虛線代表壓力分布的無黏解.可見,計算得到的壓力分布與Fang 等[18]的DNS 結(jié)果及Dupont 等[17]的實驗結(jié)果相符較好.在黏性干擾作用下,干擾區(qū)內(nèi)壓力梯度變化較為緩和且在下游逐步趨近于無黏解.

    2 脈動強度

    圖8 給出了上游湍流邊界層(turbulent boundary layer,TBL)壓力脈動強度沿法向分布情況.可見,峰值位置出現(xiàn)在近壁區(qū),約為y/δ=0.04,其量值沿物面法向急劇降低,在邊界層外緣,約為近壁峰值的28%.采用當(dāng)?shù)匚锩婕羟袘?yīng)力τw進行無量綱化后,計算得到的脈動分布與Bernardini 等[28]和Duan 等[16]的DNS 結(jié)果符合較好.

    圖8 上游湍流邊界層壓力脈動強度分布Fig.8 Distribution of pressure fluctuation intensity for the incoming TBL

    從圖9 可以看出,與上游湍流邊界層相比,干擾區(qū)內(nèi)壓力脈動強度呈現(xiàn)截然不同的分布規(guī)律.從量值上來看,脈動峰值呈現(xiàn)先急劇升高,這主要是激波的增強作用,隨后在下游再附區(qū)伴隨著擾動邊界層的恢復(fù)則呈現(xiàn)緩慢衰減的流向演化趨勢.

    圖9 干擾區(qū)壓力脈動強度云圖Fig.9 Contour of pressure fluctuation intensity in the interaction region

    為了定量描述這一現(xiàn)象,圖10 給出了物面壓力脈動強度沿流向的分布情況.計算得到脈動強度與Fang 等[20]的DNS 結(jié)果符合良好,兩者在整體往下偏移0.02 的情況下與Dupont 等[17]的實驗數(shù)據(jù)也較為吻合.這主要是由于壓力傳感器截止頻率太小導(dǎo)致實驗數(shù)據(jù)明顯偏低[17].圖中符號S和R分別表示分離點和再附點位置,通過平均摩阻的過零點確定.可見,分離區(qū)內(nèi)物面壓力脈動急劇增強,其峰值約為上游的4.3 倍,隨后在下游再附區(qū)脈動峰值緩慢降低,在x*?xr*=12 時,仍約為上游的2.2 倍.另外,從峰值位置來看,其在分離區(qū)呈現(xiàn)逐漸遠離壁面的態(tài)勢,這與激波誘導(dǎo)生成分離區(qū)剪切層相關(guān),如圖9中平均音速線所示(粉色點劃線).強脈動主要出現(xiàn)在干擾區(qū)內(nèi)壓縮波系、剪切層附近等強剪切區(qū)域.此后,在下游再附邊界層,其峰值位置仍集中在外層區(qū)域.

    圖10 物面壓力脈動強度分布Fig.10 Distribution of wall pressure fluctuation intensity

    3 功率譜

    沿外層y/δ=0.8 和物面分別選取5 個流向站位進行對比分析,其中站位1 位于上游湍流邊界層內(nèi),站位2 位于分離區(qū)內(nèi),站位3~ 5 位于下游再附邊界層內(nèi),各站位具體位置見圖9.如無特別說明,下文中外層均為y/δ=0.8 處.

    圖11 比較了上游湍流邊界層站位1 壓力脈動的預(yù)乘功率譜ωφ(ω),其中ω和φ(ω)分別為角頻率和功率譜密度.功率譜密度的計算采用Pwelch 方法和Hamming 窗函數(shù),同時對分段數(shù)據(jù)按50%重疊處理.為了便于定量比較,圖中還采用當(dāng)?shù)孛}動強度的平方進行了歸一化處理.顯而易見,外層壓力脈動的時間尺度明顯要大于物面處,前者峰值頻率約為ωδ/U∞=3.0,后者約為ωδ/U∞=12.3,這與Duan 等[13]的研究結(jié)論是一致的.這是由于外層壓力脈動主要表征為大尺度脈動結(jié)構(gòu).

    圖11 站位1 壓力脈動預(yù)乘功率譜Fig.11 Pre-multiplied power spectral density for pressure fluctuations at station 1

    圖12 分別給出了干擾區(qū)外層和物面壓力脈動預(yù)乘功率譜的變化情況.從兩者的比較來看,以下幾個方面值得特別關(guān)注.首先,相較于上游站位1,分離區(qū)站位2 處兩者都出現(xiàn)了峰值頻率向更低頻區(qū)偏移的現(xiàn)象,這與激波的低頻振蕩運動密切相關(guān).

    圖12 干擾區(qū)內(nèi)不同站位壓力脈動預(yù)乘功率譜Fig.12 Pre-multiplied power spectral density for pressure fluctuations at various streamwise locations

    如圖12(a)和圖12(b)所示,外層和物面壓力峰值頻率分別約為ωδ/U∞=0.7 和ωδ/U∞=4.4.以往研究表明[29],在分離激波低頻振蕩的作用下,分離區(qū)內(nèi)物面壓力脈動的時間尺度通常比上游湍流邊界層低1~ 2 個數(shù)量級.本文計算結(jié)果進一步驗證了該結(jié)論,同時結(jié)果還表明,外層壓力脈動也同樣符合這一規(guī)律.在功率譜的定性分布方面,顯然,外層壓力脈動的變化較物面壓力更為劇烈,這很可能是由于外層站位非常靠近分離激波,受其非定常運動的影響更為直接,而物面壓力處于分離泡內(nèi)(見圖9),影響則要弱得多.

    從再附區(qū)站位3~ 5 的結(jié)果來看,外層壓力脈動呈現(xiàn)快速恢復(fù)的態(tài)勢,且峰值頻率向高頻區(qū)偏移,逐步趨近于上游站位1,這與物面壓力脈動的演化過程則完全不同.從圖12(a)中可以清楚看到,再附邊界層內(nèi)物面壓力脈動的低頻能量占比仍相對較高,峰值頻率依然維持在較低頻區(qū),并沒有出現(xiàn)往高頻區(qū)恢復(fù)的現(xiàn)象.造成兩者該差異的原因可以認為主要有兩方面:一方面,上游激波誘導(dǎo)分離形成了剪切層,分離泡上方剪切層內(nèi)不穩(wěn)定渦結(jié)構(gòu)沿流向急劇發(fā)展,并進入到下游再附邊界層內(nèi),這些不穩(wěn)定結(jié)構(gòu)極大地抑制了物面壓力脈動結(jié)構(gòu)的恢復(fù);另一方面,外層壓力脈動法向位置相對較高,盡管在站位2 處激波干擾對其影響顯著,但是過激波后外層壓力脈動受近壁區(qū)流動結(jié)構(gòu)的影響則要小得多,主要表征為外層大尺度脈動結(jié)構(gòu)的破碎過程,更為詳細的定量分析可見下一節(jié)的兩點相關(guān)統(tǒng)計結(jié)果.

    4 兩點相關(guān)

    通過兩點相關(guān)性考察干擾區(qū)外層壓力脈動擬序結(jié)構(gòu)長度尺度的變化規(guī)律,且定量分析比較其與物面壓力脈動結(jié)構(gòu)尺度的差異.兩點相關(guān)系數(shù)Cpp(Δx,Δz)定義為

    式中p′為雷諾平均后得到的壓力脈動,Δx和Δz分別代表流向間距和展向間距,x0為特征點流向位置,yref分別取物面和外層,上橫線代表時間和展向平均.

    圖13 分別給出了上游湍流邊界層外層和物面脈動壓力兩點相關(guān)系數(shù)的分布情況.以往的零壓力梯度平板結(jié)果表明[28],小間距時脈動壓力相關(guān)性曲線可近似為圓形,這是小尺度壓力脈動結(jié)構(gòu)的各向同性所決定的;大間距時,大尺度脈動結(jié)構(gòu)的各向異性導(dǎo)致相關(guān)性曲線在展向出現(xiàn)了急劇拉伸,近似為橢圓形分布,表征了展向拉長的大尺度結(jié)構(gòu).顯然,圖13(a)中的結(jié)果也完全符合該規(guī)律.然而,對于外層壓力脈動,如圖13(b)所示,整體來看,外層脈動結(jié)構(gòu)的流向尺度和展向尺度都明顯大于物面壓力脈動.這里以相關(guān)性系數(shù)0.3 為參考確定其特征長度,外層結(jié)構(gòu)的流向和展向尺度分別約為0.7δ和1.0δ,而物面壓力脈動則分別約為0.14δ和0.2δ.此外,盡管外層相關(guān)性等值線小間距時仍趨近于圓形,但隨著間距增大,等值線在流向和展向都被急劇拉伸,展向明顯強于流向,這表明外層壓力脈動擬序結(jié)構(gòu)仍以展向大尺度結(jié)構(gòu)為主,但其特征尺度要大得多.

    圖13 站位1 壓力脈動兩點相關(guān)系數(shù)Fig.13 Two-point correlation coefficient of pressure fluctuations at station 1

    圖14(a)和圖14(b)分別給出了站位2 物面和外層壓力脈動的兩點相關(guān)等值線,其中等值線從外到內(nèi)依次為0.3,0.5,0.7 和0.9.為了便于比較說明,圖中也給出了站位1 的結(jié)果.從分布規(guī)律上來看,等值線形態(tài)基本類似于上游湍流邊界層,仍然以展向拉伸的橢圓形分布為主,但整體呈現(xiàn)向外擴張的趨勢,這表明激波干擾作用增大了外層和物面壓力脈動結(jié)構(gòu)尺度,特別是在展向.需要特別關(guān)注的是,從圖14(b)中可以清楚看到,同樣以等值線0.3 為參考,站位1 處的展向尺度約為1.0δ,而站位2 處增大到約2.0δ,而流向尺度兩者基本維持在0.7δ.這很可能是由于站位2 外層位于入射激波附近(見圖9),激波的展向結(jié)構(gòu)增強了壓力脈動的展向相關(guān)性.更為細致的定量分析還有待深入研究.

    圖14 站位2 兩點相關(guān)系數(shù)(虛線:站位1)Fig.14 Two-point correlation coefficient of pressure fluctuations at station 2(dashed lines:station 1)

    如圖15 所示,再附區(qū)站位3~ 5 外層和物面壓力脈動的兩點相關(guān)等值線與上游湍流邊界層也基本類似.小間距時,曲線接近于圓形,這表明小尺度結(jié)構(gòu)并沒有發(fā)生實質(zhì)改變,仍以各向同性為主.從大尺度結(jié)構(gòu)的演化歷程來看,兩者差異顯著.研究發(fā)現(xiàn),大間距時,圖15(a)~ 圖15(c)中等值線以向外擴張為主,這表明再附區(qū)內(nèi)物面壓力大尺度脈動結(jié)構(gòu)的特征尺度仍在增大,而圖15(d)~ 圖15(f)等值線則以向內(nèi)收縮為主,對應(yīng)的是外層壓力大尺度脈動結(jié)構(gòu)特征尺度的減小.

    圖15 站位3~ 5 兩點相關(guān)系數(shù)(上:物面;下:y/δ=0.8;虛線:站位1)Fig.15 Two-point correlation coefficient at station 3~ 5 (top:wall;down:y/δ=0.8;dashed lines:station 1)

    為了定量描述該現(xiàn)象,圖16 進一步給出了干擾區(qū)內(nèi)物面及外層壓力脈動積分尺度的分布.積分尺度可以用來表征脈動場最大尺度結(jié)構(gòu)的特征尺度,這里流向和展向積分尺度分別定義為[10]

    可見,外層壓力脈動流向和展向積分尺度遠大于物面壓力脈動,這與Duan 等[16]的研究結(jié)果是一致的.從整體變化規(guī)律上來看,外層壓力Λx和Λz在站位2 處急劇增大,隨后在站位3~ 5 則又緩慢減小,而物面壓力的Λx和Λz整體上呈現(xiàn)逐步增大的趨勢,這與圖15 中相關(guān)性分析結(jié)果是相符的.從量值變化上來看,以圖16(b)中的展向積分尺度Λz為例,對于外層壓力,在激波的增強作用下,站位1 處Λz=0.73δ增大到站位2 處Λz=1.17δ,隨后在再附區(qū)又緩慢衰減到站位5 的Λz=1.02δ;而物面壓力,從站位1 處Λz=0.26δ逐步增大到站位5 處的Λz=0.82δ,后者約為前者的3.15 倍.可以看到,在站位5 處,此時再附區(qū)內(nèi)物面壓力脈動結(jié)構(gòu)的展向和流向積分尺度已非常接近于外層的積分尺度.綜上所述,可以認為這是由于外層壓力脈動大尺度結(jié)構(gòu)在過激波后呈現(xiàn)快速破碎的過程,而物面上大尺度結(jié)構(gòu)仍在緩慢增長,這也解釋了圖12 中物面壓力脈動低頻能量的相對升高而外層壓力脈動峰值頻率往高頻區(qū)的偏移現(xiàn)象.

    圖16 壓力脈動積分尺度分布Fig.16 Distribution of integral scale for pressure fluctuations

    5時空關(guān)聯(lián)

    進一步分析干擾區(qū)外層和物面壓力脈動的時空關(guān)聯(lián)特性,從而得到向下游傳播的壓力波的對流速度.

    圖17 分別給出了站位1 處物面和外層壓力的時空關(guān)聯(lián)云圖,圖中時空關(guān)聯(lián)系數(shù)Rpp(Δx,Δt)定義為

    式中Δt為延遲時間.以往研究結(jié)果表明[28,30],由于脈動結(jié)構(gòu)往下游傳播,物面壓力時空關(guān)聯(lián)云圖呈現(xiàn)傾斜的類橢圓形分布,如圖17(a)所示,本文計算結(jié)果也完全符合該結(jié)論.相較于物面壓力,從圖17(b)中可以看到,外層壓力時空關(guān)聯(lián)云圖的橢圓率和傾斜角都明顯增大.

    圖17 站位1 壓力脈動時空關(guān)聯(lián)系數(shù)Fig.17 Space-time correlation coefficient for pressure fluctuations at station 1

    圖18 和圖19 分別給出了分離區(qū)站位2 和再附區(qū)站位3~ 5 壓力脈動時空關(guān)聯(lián)系數(shù)與站位1 (黑色虛線)的比較情況.如圖18(a)所示,可以發(fā)現(xiàn),相較于上游湍流邊界層,分離區(qū)物面脈動壓力時空關(guān)聯(lián)曲線更為聚集,同時傾斜角也急劇減小,但整體形態(tài)仍滿足橢圓形(只是等值線0.1 較為雜亂).值得注意的是,對于站位2 外層脈動壓力,顯然,圖18(b)中等值線形態(tài)變化非常劇烈,已完全不符合橢圓形分布,特別是等值線0.1 和0.3,在時間方向上可近似為直線分布.造成這個現(xiàn)象的原因很可能是由于站位2 外層處于強激波區(qū)(見圖9),在分離激波的強干擾下,此時泰勒凍結(jié)假設(shè)完全失效,脈動壓力在時空上表征為強非線性.此外,從圖19 中還可以看到,無論是再附區(qū)外層還是物面,其脈動壓力時空關(guān)聯(lián)等值線整體形態(tài)與上游站位1 都較為類似.從站位3~站位5,強相關(guān)等值線被急劇拉長壓縮,且傾斜角增大,呈現(xiàn)逐步逼近上游站位1 的演化過程.從定性比較的角度來看,外層脈動壓力時空關(guān)聯(lián)曲線的恢復(fù)速度明顯要快于物面.

    圖18 站位2 壓力脈動時空關(guān)聯(lián)系數(shù)(虛線:站位1)Fig.18 Space-time correlation coefficient for pressure fluctuations at station 2 (dashed lines:station 1)

    圖19 站位3~ 5 壓力脈動時空關(guān)聯(lián)系數(shù)(上:物面;下:y/δ=0.8;虛線:站位1)Fig.19 Space-time correlation coefficient at station 3~ 5 (top:wall;down:y/δ=0.8;dashed lines:station 1)

    壓力脈動對流速度的計算方法如下[16,31]:首先給定延遲時間Δt,隨后依據(jù)時空關(guān)聯(lián)曲線通過

    獲得流向間距rx,進而得到對流速度Uc=rx/Δt.采用上述方法計算得到了上游湍流邊界層物面和外層脈動壓力對流速度隨延遲時間的分布情況.如圖20 和圖21 所示.總體來看,外層脈動壓力對流速度要遠大于物面脈動壓力,這與Duan 等[16]的研究結(jié)果也較為符合.另外,對于外層脈動壓力,可見,其對流速度對延遲時間的依賴性較小,基本維持在Uc=0.94U∞,與當(dāng)?shù)仄骄飨蛩俣冉咏?然而,物面脈動壓力對流速度隨延遲時間則變化劇烈,在小延遲時間(小尺度脈動)時,對流速度約為Uc=0.56U∞,而在大延遲時間(大尺度脈動) 時,其增大到Uc=0.73U∞.已有研究表明[13],這是由于物面大尺度脈動壓力與外層大尺度結(jié)構(gòu)相關(guān),而小尺度脈動壓力則主要受近壁區(qū)結(jié)構(gòu)影響,故前者對流速度明顯大于后者.

    從圖20 和圖21 中還可以看到,分離區(qū)物面脈動壓力對流速度降低,約為Uc=0.2U∞,隨后下游再附區(qū)則又逐漸增大,但仍小于上游湍流邊界層.Na等[32]在強逆壓力梯度平板的DNS 結(jié)果中也發(fā)現(xiàn)了類似規(guī)律,本文的計算結(jié)果進一步證實了該現(xiàn)象在入射激波/湍流邊界層干擾區(qū)內(nèi)的適用性.從圖20(b)中可以看到,采用當(dāng)?shù)赝鈱铀俣萓e對Uc進行歸一化后,分離區(qū)和再附區(qū)各站位處物面脈動壓力的對流速度分布均吻合較好,約為Uc=0.25Ue,與Bernardini 等[10]在跨聲速激波/湍流邊界層干擾下的研究結(jié)論也較為接近,但后者為初始分離(incipient separation)工況.如圖21 所示,干擾區(qū)外層不同站位處脈動壓力對流速度分布與上游結(jié)果較為類似,其對延遲時間的依賴性依然較弱,但其量值減小到約為Uc=0.4Ue.結(jié)合之前的功率譜和兩點相關(guān)分析,可以推測外層大尺度脈動結(jié)構(gòu)特征尺度的減小可能是造成該現(xiàn)象的重要因素.

    圖20 干擾區(qū)內(nèi)各站位物面壓力脈動對流速度Fig.20 Convection velocity for wall pressure fluctuations at various stations in the interaction region

    圖21 干擾區(qū)內(nèi)各站位外層壓力脈動對流速度Fig.21 Convection velocity for pressure fluctuations in the outer layer

    6 結(jié)論

    本文采用直接數(shù)值模擬方法研究了來流馬赫數(shù)2.25 下激波角33.2°的入射激波對平板湍流邊界層壓力脈動特性的影響規(guī)律,詳細地分析了外層和物面壓力脈動場的典型統(tǒng)計特征,如脈動強度、功率譜密度、兩點相關(guān)和時空關(guān)聯(lián)特性等,得到以下結(jié)論:

    (1) 數(shù)值模擬準確捕捉到了干擾區(qū)內(nèi)脈動壓力.研究發(fā)現(xiàn),分離區(qū)外層脈動壓力以低頻能量為主,高頻能量占比較小.與物面脈動壓力不同的是,再附區(qū)外層脈動壓力峰值頻率逐漸往高頻區(qū)偏移,而物面壓力脈動能量仍維持在相對較低頻率.

    (2) 激波干擾對壓力脈動結(jié)構(gòu)有顯著的增強作用,干擾區(qū)內(nèi)展向積分尺度明顯大于流向積分尺度.積分尺度在外層和物面的演化歷程差異明顯,前者經(jīng)歷了先急劇增長后緩慢衰減的過程,而后者整體上呈現(xiàn)增長趨勢.

    (3) 激波干擾極大地降低了壓力脈動場的對流速度.時空關(guān)聯(lián)結(jié)果表明,干擾區(qū)下游外層壓力波的對流速度仍明顯高于物面,前者約為0.4Ue,而后者約為0.25Ue.

    猜你喜歡
    物面外層邊界層
    一種溶液探測傳感器
    傳感器世界(2022年4期)2022-11-24 21:23:50
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    讓吸盤掛鉤更牢固
    一種購物袋
    科技資訊(2016年6期)2016-05-14 13:09:55
    新型單面陣自由曲面光學(xué)測量方法成像特性仿真
    一類具有邊界層性質(zhì)的二次奇攝動邊值問題
    專題Ⅱ 物質(zhì)構(gòu)成的奧秘
    “人”字變身
    賽車空動優(yōu)化的秘密
    汽車之友(2014年21期)2014-11-03 17:33:45
    非特征邊界的MHD方程的邊界層
    国产黄色免费在线视频| 国产精品国产av在线观看| 国产无遮挡羞羞视频在线观看| 超碰成人久久| e午夜精品久久久久久久| 久久人妻福利社区极品人妻图片| 亚洲伊人色综图| 女性被躁到高潮视频| www.熟女人妻精品国产| 国产成人影院久久av| 国产欧美日韩一区二区三| 久久精品91蜜桃| 99在线视频只有这里精品首页| 女人高潮潮喷娇喘18禁视频| 欧美乱码精品一区二区三区| 免费av中文字幕在线| 别揉我奶头~嗯~啊~动态视频| 黄色视频,在线免费观看| 麻豆久久精品国产亚洲av | 波多野结衣av一区二区av| 身体一侧抽搐| www.自偷自拍.com| 80岁老熟妇乱子伦牲交| 日本精品一区二区三区蜜桃| 天天添夜夜摸| 电影成人av| av在线天堂中文字幕 | 美女高潮到喷水免费观看| 免费av毛片视频| 精品乱码久久久久久99久播| av天堂在线播放| 亚洲精品成人av观看孕妇| 亚洲午夜理论影院| 法律面前人人平等表现在哪些方面| 成人黄色视频免费在线看| 男女下面插进去视频免费观看| 日韩精品免费视频一区二区三区| 国产视频一区二区在线看| 国产成人精品久久二区二区免费| 91麻豆av在线| 亚洲色图 男人天堂 中文字幕| 色综合站精品国产| 久久久国产精品麻豆| 国产精品九九99| www.999成人在线观看| 999久久久国产精品视频| 老司机在亚洲福利影院| 高清欧美精品videossex| 可以在线观看毛片的网站| 色尼玛亚洲综合影院| 久久精品国产99精品国产亚洲性色 | 又紧又爽又黄一区二区| 国产成人av激情在线播放| 视频区图区小说| 嫁个100分男人电影在线观看| 午夜免费激情av| 久久精品成人免费网站| 最近最新中文字幕大全免费视频| 国产一区二区激情短视频| 亚洲av美国av| 久久精品影院6| 精品一品国产午夜福利视频| 身体一侧抽搐| 老熟妇乱子伦视频在线观看| 国产激情欧美一区二区| 日本三级黄在线观看| 久久亚洲真实| 午夜福利欧美成人| www.自偷自拍.com| 亚洲人成网站在线播放欧美日韩| 看片在线看免费视频| 亚洲欧洲精品一区二区精品久久久| 亚洲人成电影免费在线| 亚洲一区二区三区欧美精品| 中文字幕另类日韩欧美亚洲嫩草| 日本 av在线| 新久久久久国产一级毛片| 亚洲成人久久性| 国产无遮挡羞羞视频在线观看| 国产精品久久久久成人av| 久久青草综合色| 亚洲欧洲精品一区二区精品久久久| aaaaa片日本免费| 国产极品粉嫩免费观看在线| a级毛片在线看网站| 国产视频一区二区在线看| 午夜福利影视在线免费观看| 午夜91福利影院| 两人在一起打扑克的视频| 99riav亚洲国产免费| 久久人妻熟女aⅴ| 亚洲精品粉嫩美女一区| 人人妻人人爽人人添夜夜欢视频| 亚洲精品一区av在线观看| 黄网站色视频无遮挡免费观看| 看黄色毛片网站| 亚洲在线自拍视频| 老司机午夜福利在线观看视频| 成人黄色视频免费在线看| 免费搜索国产男女视频| 99国产综合亚洲精品| 19禁男女啪啪无遮挡网站| 欧美老熟妇乱子伦牲交| 亚洲自拍偷在线| 夜夜看夜夜爽夜夜摸 | 女人爽到高潮嗷嗷叫在线视频| 可以免费在线观看a视频的电影网站| 国产精华一区二区三区| 久久人妻熟女aⅴ| 天堂中文最新版在线下载| 母亲3免费完整高清在线观看| 成人亚洲精品一区在线观看| 久久人妻熟女aⅴ| 中亚洲国语对白在线视频| a级毛片在线看网站| 国产精品久久视频播放| 在线观看免费日韩欧美大片| 国产精品久久视频播放| 黄片小视频在线播放| 亚洲伊人色综图| 亚洲精品一区av在线观看| 亚洲国产精品sss在线观看 | 无遮挡黄片免费观看| 国产一区二区三区在线臀色熟女 | 亚洲va日本ⅴa欧美va伊人久久| 欧美亚洲日本最大视频资源| 啦啦啦免费观看视频1| 香蕉丝袜av| 视频区图区小说| 欧美最黄视频在线播放免费 | 欧美亚洲日本最大视频资源| 91字幕亚洲| 久热这里只有精品99| 99热只有精品国产| 欧美日韩精品网址| 欧美成人午夜精品| 婷婷精品国产亚洲av在线| 亚洲国产精品999在线| 80岁老熟妇乱子伦牲交| 99精品欧美一区二区三区四区| 中文字幕人妻丝袜一区二区| 成年人免费黄色播放视频| 欧美黑人欧美精品刺激| 99在线人妻在线中文字幕| 国产免费现黄频在线看| 亚洲av第一区精品v没综合| 色精品久久人妻99蜜桃| 国产色视频综合| 最新在线观看一区二区三区| 亚洲精品国产一区二区精华液| 欧美日韩一级在线毛片| 欧美日韩一级在线毛片| 男人的好看免费观看在线视频 | 在线永久观看黄色视频| 亚洲国产看品久久| 男女下面进入的视频免费午夜 | 一区二区日韩欧美中文字幕| 亚洲国产精品sss在线观看 | 亚洲av成人av| 91麻豆精品激情在线观看国产 | av片东京热男人的天堂| 日本一区二区免费在线视频| 女警被强在线播放| 国产国语露脸激情在线看| 国产极品粉嫩免费观看在线| 久99久视频精品免费| 视频区图区小说| 很黄的视频免费| 日本wwww免费看| 亚洲免费av在线视频| 国产区一区二久久| www国产在线视频色| 成熟少妇高潮喷水视频| 成熟少妇高潮喷水视频| 丝袜在线中文字幕| 免费av毛片视频| 91av网站免费观看| 少妇裸体淫交视频免费看高清 | 老汉色∧v一级毛片| 国产一区二区三区视频了| 美女午夜性视频免费| 欧美激情极品国产一区二区三区| 在线看a的网站| 在线观看日韩欧美| 色综合欧美亚洲国产小说| 国产精品一区二区免费欧美| 夜夜躁狠狠躁天天躁| 亚洲欧美一区二区三区久久| 国产精品乱码一区二三区的特点 | 免费日韩欧美在线观看| 亚洲午夜精品一区,二区,三区| 亚洲五月色婷婷综合| 动漫黄色视频在线观看| 少妇粗大呻吟视频| 老熟妇仑乱视频hdxx| 老司机靠b影院| 一级a爱片免费观看的视频| 在线观看一区二区三区激情| 午夜a级毛片| 免费在线观看完整版高清| 久久中文字幕人妻熟女| 日韩欧美免费精品| 黄频高清免费视频| 真人一进一出gif抽搐免费| svipshipincom国产片| 老司机深夜福利视频在线观看| 久久久久亚洲av毛片大全| 亚洲精品美女久久久久99蜜臀| 中国美女看黄片| 国产成人欧美在线观看| 亚洲色图综合在线观看| 亚洲中文日韩欧美视频| 老司机福利观看| 日日爽夜夜爽网站| 亚洲欧美激情在线| av天堂在线播放| 淫秽高清视频在线观看| 日韩大码丰满熟妇| 性少妇av在线| 久久伊人香网站| 怎么达到女性高潮| 亚洲专区国产一区二区| 黑人巨大精品欧美一区二区蜜桃| 成人国产一区最新在线观看| 国产成人系列免费观看| 亚洲精品美女久久av网站| 人成视频在线观看免费观看| 青草久久国产| 法律面前人人平等表现在哪些方面| 久久国产亚洲av麻豆专区| 夜夜看夜夜爽夜夜摸 | 亚洲色图av天堂| 精品一区二区三区视频在线观看免费 | 免费观看人在逋| 久久性视频一级片| av片东京热男人的天堂| 色尼玛亚洲综合影院| 真人做人爱边吃奶动态| 久久精品亚洲熟妇少妇任你| 精品国产超薄肉色丝袜足j| 涩涩av久久男人的天堂| 午夜福利一区二区在线看| 久久久久久久午夜电影 | 国产人伦9x9x在线观看| 国产伦一二天堂av在线观看| 一区二区三区激情视频| 村上凉子中文字幕在线| 黄片小视频在线播放| 国产精品av久久久久免费| 大香蕉久久成人网| 大码成人一级视频| 黄色怎么调成土黄色| 在线十欧美十亚洲十日本专区| 精品国产亚洲在线| 亚洲国产精品999在线| 99国产精品免费福利视频| 中文字幕另类日韩欧美亚洲嫩草| 免费看十八禁软件| 最新在线观看一区二区三区| 国产精品一区二区在线不卡| 91精品国产国语对白视频| 一二三四在线观看免费中文在| 三上悠亚av全集在线观看| 国产不卡一卡二| 91av网站免费观看| 在线观看日韩欧美| 级片在线观看| 中文字幕精品免费在线观看视频| 亚洲中文字幕日韩| 国产精品香港三级国产av潘金莲| 日韩欧美三级三区| 亚洲av熟女| 一进一出抽搐动态| 黄色丝袜av网址大全| 免费日韩欧美在线观看| 日本五十路高清| 国产成年人精品一区二区 | 男男h啪啪无遮挡| 欧美亚洲日本最大视频资源| 国产一区在线观看成人免费| 午夜福利免费观看在线| 别揉我奶头~嗯~啊~动态视频| 亚洲精品国产色婷婷电影| 亚洲欧美日韩另类电影网站| 一区福利在线观看| 久久亚洲真实| 国产精品久久电影中文字幕| 日韩免费高清中文字幕av| 老汉色∧v一级毛片| 日本 av在线| 国产精品自产拍在线观看55亚洲| 中文欧美无线码| 淫妇啪啪啪对白视频| 一级黄色大片毛片| 亚洲第一av免费看| 久久中文字幕一级| 热99re8久久精品国产| 精品久久久久久久久久免费视频 | 大型黄色视频在线免费观看| 久久天堂一区二区三区四区| 电影成人av| 一级片'在线观看视频| 久久久久久久久中文| 天天躁夜夜躁狠狠躁躁| 亚洲精品国产精品久久久不卡| 日韩人妻精品一区2区三区| 日本vs欧美在线观看视频| 丁香六月欧美| 老司机亚洲免费影院| 亚洲国产毛片av蜜桃av| 亚洲九九香蕉| 黄色片一级片一级黄色片| 国产成人av激情在线播放| 日韩中文字幕欧美一区二区| 亚洲精品一二三| 国产精品美女特级片免费视频播放器 | 精品国产美女av久久久久小说| 亚洲av日韩精品久久久久久密| 69精品国产乱码久久久| 国产成人av激情在线播放| 亚洲五月色婷婷综合| 天天影视国产精品| 免费在线观看完整版高清| 国产黄色免费在线视频| 啦啦啦在线免费观看视频4| 免费看a级黄色片| 91av网站免费观看| e午夜精品久久久久久久| 黑人欧美特级aaaaaa片| 国产亚洲av高清不卡| avwww免费| 久久99一区二区三区| 最新在线观看一区二区三区| 亚洲中文日韩欧美视频| 三上悠亚av全集在线观看| 欧美激情极品国产一区二区三区| 日韩大码丰满熟妇| 制服人妻中文乱码| 亚洲五月色婷婷综合| 免费高清视频大片| 欧美精品一区二区免费开放| 亚洲精品国产色婷婷电影| 少妇 在线观看| 丰满的人妻完整版| 免费人成视频x8x8入口观看| 欧美丝袜亚洲另类 | 高清av免费在线| 国产熟女xx| 99久久人妻综合| videosex国产| 精品熟女少妇八av免费久了| 怎么达到女性高潮| 国产一区二区三区视频了| 欧美激情 高清一区二区三区| 啪啪无遮挡十八禁网站| 校园春色视频在线观看| 久久久久国产精品人妻aⅴ院| 亚洲美女黄片视频| av在线天堂中文字幕 | 国产aⅴ精品一区二区三区波| 亚洲色图av天堂| 老汉色av国产亚洲站长工具| 精品一区二区三区四区五区乱码| 亚洲精品久久午夜乱码| 色在线成人网| 久久久久久久精品吃奶| 欧美中文综合在线视频| 欧美人与性动交α欧美软件| 久久亚洲精品不卡| 久久精品成人免费网站| 一级毛片精品| av中文乱码字幕在线| 高清黄色对白视频在线免费看| 美女扒开内裤让男人捅视频| 水蜜桃什么品种好| 国产精品国产av在线观看| 夜夜爽天天搞| 成人18禁高潮啪啪吃奶动态图| 少妇 在线观看| 99精品欧美一区二区三区四区| 真人一进一出gif抽搐免费| 久久国产亚洲av麻豆专区| 村上凉子中文字幕在线| 性色av乱码一区二区三区2| 成人亚洲精品av一区二区 | 亚洲人成77777在线视频| 久久久国产成人精品二区 | 长腿黑丝高跟| 亚洲人成伊人成综合网2020| 超碰成人久久| 精品乱码久久久久久99久播| 18禁裸乳无遮挡免费网站照片 | 欧美乱码精品一区二区三区| 欧美+亚洲+日韩+国产| 国产免费现黄频在线看| 国产精品秋霞免费鲁丝片| netflix在线观看网站| 婷婷精品国产亚洲av在线| 亚洲欧美一区二区三区黑人| 国产aⅴ精品一区二区三区波| 黑人巨大精品欧美一区二区mp4| 国产一区二区三区综合在线观看| 99久久国产精品久久久| 久久精品亚洲熟妇少妇任你| 国产免费av片在线观看野外av| 青草久久国产| 啦啦啦在线免费观看视频4| 老司机靠b影院| videosex国产| 免费不卡黄色视频| 咕卡用的链子| 在线观看免费视频网站a站| 黄片大片在线免费观看| 亚洲精品久久午夜乱码| 午夜精品国产一区二区电影| 亚洲性夜色夜夜综合| 国产成人欧美在线观看| 亚洲五月色婷婷综合| 成年人免费黄色播放视频| 亚洲免费av在线视频| 亚洲色图av天堂| 99在线人妻在线中文字幕| 亚洲免费av在线视频| 国产乱人伦免费视频| 黄片小视频在线播放| 亚洲人成伊人成综合网2020| 又黄又爽又免费观看的视频| 婷婷六月久久综合丁香| 国产av又大| 精品国产国语对白av| 久久亚洲精品不卡| 香蕉国产在线看| 精品国产国语对白av| 男人舔女人下体高潮全视频| 国产精品成人在线| 高潮久久久久久久久久久不卡| 久久人人97超碰香蕉20202| 黄频高清免费视频| 久久久久精品国产欧美久久久| 妹子高潮喷水视频| 免费高清在线观看日韩| 成熟少妇高潮喷水视频| 丰满饥渴人妻一区二区三| 国产欧美日韩一区二区三| 丝袜美腿诱惑在线| 正在播放国产对白刺激| 男人舔女人的私密视频| 日日干狠狠操夜夜爽| 精品久久久久久久久久免费视频 | 国产男靠女视频免费网站| 国产成人精品在线电影| 久久国产精品人妻蜜桃| 91成人精品电影| 男人的好看免费观看在线视频 | 免费观看精品视频网站| 搡老岳熟女国产| 精品熟女少妇八av免费久了| 亚洲国产精品合色在线| 国产一区二区激情短视频| 午夜久久久在线观看| 嫩草影视91久久| 色精品久久人妻99蜜桃| netflix在线观看网站| 亚洲成人国产一区在线观看| 在线观看日韩欧美| 手机成人av网站| 琪琪午夜伦伦电影理论片6080| 婷婷丁香在线五月| 999久久久精品免费观看国产| 欧美在线一区亚洲| 日日夜夜操网爽| 亚洲精品久久成人aⅴ小说| 一进一出抽搐动态| 国产精品一区二区精品视频观看| 国产欧美日韩精品亚洲av| 色婷婷av一区二区三区视频| 老司机在亚洲福利影院| 久久久精品欧美日韩精品| 国产av一区在线观看免费| 午夜影院日韩av| 热99re8久久精品国产| 亚洲熟女毛片儿| 日本免费a在线| 精品国产乱子伦一区二区三区| 欧美av亚洲av综合av国产av| 美女 人体艺术 gogo| 午夜两性在线视频| 在线看a的网站| 午夜日韩欧美国产| 国产蜜桃级精品一区二区三区| 精品国产一区二区久久| 亚洲欧美日韩另类电影网站| 免费女性裸体啪啪无遮挡网站| 亚洲精品久久午夜乱码| 精品人妻在线不人妻| 国产精品秋霞免费鲁丝片| 国产精品1区2区在线观看.| 亚洲性夜色夜夜综合| 久久久久久久久免费视频了| 一区福利在线观看| av片东京热男人的天堂| 老汉色∧v一级毛片| 亚洲五月天丁香| 日韩精品中文字幕看吧| 99国产综合亚洲精品| 91大片在线观看| 亚洲精品久久成人aⅴ小说| 欧美日韩一级在线毛片| 在线十欧美十亚洲十日本专区| 一二三四社区在线视频社区8| 欧美另类亚洲清纯唯美| 亚洲精品一区av在线观看| 精品熟女少妇八av免费久了| av免费在线观看网站| 国内毛片毛片毛片毛片毛片| 欧美最黄视频在线播放免费 | 中文字幕人妻丝袜一区二区| 欧美日韩亚洲国产一区二区在线观看| 久久久久久人人人人人| 这个男人来自地球电影免费观看| 美国免费a级毛片| 国产欧美日韩精品亚洲av| 久久久国产一区二区| 精品电影一区二区在线| av在线天堂中文字幕 | 亚洲国产精品合色在线| 一个人免费在线观看的高清视频| 精品国产乱子伦一区二区三区| 99re在线观看精品视频| 久久这里只有精品19| 每晚都被弄得嗷嗷叫到高潮| 免费看十八禁软件| 色婷婷久久久亚洲欧美| 国产亚洲精品久久久久5区| 在线国产一区二区在线| 免费高清视频大片| 丰满迷人的少妇在线观看| 亚洲片人在线观看| 女人被狂操c到高潮| 精品电影一区二区在线| 国产成人影院久久av| 欧美性长视频在线观看| 亚洲狠狠婷婷综合久久图片| 脱女人内裤的视频| 久久天堂一区二区三区四区| 自拍欧美九色日韩亚洲蝌蚪91| 99久久久亚洲精品蜜臀av| 日韩高清综合在线| 又黄又爽又免费观看的视频| 淫秽高清视频在线观看| 欧美大码av| 80岁老熟妇乱子伦牲交| 老司机福利观看| 天天影视国产精品| 十分钟在线观看高清视频www| 18禁观看日本| 一区在线观看完整版| 亚洲中文字幕日韩| 这个男人来自地球电影免费观看| 97人妻天天添夜夜摸| www.www免费av| 亚洲欧美日韩另类电影网站| 久久久久久久久久久久大奶| 757午夜福利合集在线观看| 亚洲九九香蕉| 9191精品国产免费久久| 国产成人精品在线电影| 精品国内亚洲2022精品成人| 日韩欧美三级三区| 91老司机精品| 国产精品一区二区三区四区久久 | 亚洲午夜理论影院| 色综合站精品国产| 丁香六月欧美| 五月开心婷婷网| 日日夜夜操网爽| 99在线人妻在线中文字幕| 黑人欧美特级aaaaaa片| 精品国产美女av久久久久小说| 黑人猛操日本美女一级片| 久久国产精品男人的天堂亚洲| 首页视频小说图片口味搜索| 亚洲自拍偷在线| www国产在线视频色| 国产在线观看jvid| 欧美成人免费av一区二区三区| 久久久国产欧美日韩av| 亚洲成国产人片在线观看| 久久人人精品亚洲av| 国产亚洲欧美98| 亚洲自拍偷在线| 成人三级做爰电影| 久久亚洲精品不卡| 人妻久久中文字幕网| 日韩 欧美 亚洲 中文字幕| 婷婷丁香在线五月| 国产又爽黄色视频| 高清av免费在线| 久久国产精品男人的天堂亚洲| 欧美丝袜亚洲另类 | 亚洲五月天丁香| 黄色成人免费大全| 久久性视频一级片| 久久久久久免费高清国产稀缺| 成人黄色视频免费在线看| 女人被躁到高潮嗷嗷叫费观| 三级毛片av免费| 国产野战对白在线观看| 大码成人一级视频| 国产蜜桃级精品一区二区三区| 中文字幕最新亚洲高清| 人人妻人人爽人人添夜夜欢视频| 国产日韩一区二区三区精品不卡| 欧美黄色片欧美黄色片| 亚洲精品久久午夜乱码|