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

    基于RANS/LES方法的超聲速底部流場數(shù)值模擬

    2017-11-23 05:57:05張露李杰
    航空學(xué)報 2017年1期
    關(guān)鍵詞:超聲速湍流流場

    張露,李杰

    基于RANS/LES方法的超聲速底部流場數(shù)值模擬

    張露,李杰*

    西北工業(yè)大學(xué) 航空學(xué)院,西安 710072

    分別采用基于兩方程k-ω剪切應(yīng)力輸運(yùn)(SST)湍流模型的延遲DES(DDES)、更改的DDES(MDDES)和改進(jìn)的DDES(IDDES)方法,并引入可壓縮修正,結(jié)合三階MUSCL-Roe和五階 WENO-Roe兩種空間離散格式,針對超聲速底部的復(fù)雜流動現(xiàn)象,開展了數(shù)值模擬研究。計算結(jié)果表明本文方法能夠捕捉到超聲速底部流動中豐富的湍流結(jié)構(gòu),通過分析計算結(jié)果對超聲速底部的流動機(jī)理有了進(jìn)一步的認(rèn)識,為下一步的超聲速底部流動減阻改進(jìn)和雷諾平均Navier-Stokes/大渦模擬(RANS/LES)方法在非定常高可壓縮性流動中的應(yīng)用提供了參考。通過對比分析不同空間離散格式的計算結(jié)果研究了數(shù)值耗散對計算的影響,五階WENO-Roe格式的計算結(jié)果與實驗結(jié)果吻合良好;對不同RANS/LES混合方法的計算結(jié)果進(jìn)行了對比分析,結(jié)果表明IDDES方法在近壁區(qū)的表現(xiàn)優(yōu)于DDES和MDDES方法。

    RANS/LES混合方法;超聲速底部流動;壓縮性效應(yīng);湍流模型;數(shù)值模擬

    炮彈或?qū)椀燃?xì)長體類飛行器在超聲速飛行中會受到波阻、底阻和摩阻三方面阻力的影響,而其中底阻所占的比重最大,要改進(jìn)飛行品質(zhì),準(zhǔn)確計算飛行軌跡,增大射程,就要對超聲速底部流動進(jìn)行全面準(zhǔn)確的分析研究,從而采取有針對性的措施減小底阻的不利影響。

    近年來研究者們已經(jīng)利用實驗手段對超聲速底部流動開展了一些有價值的工作。Herrin和Dutton[1]通過風(fēng)洞試驗對馬赫數(shù)為2.5的超聲速底部流場近尾跡區(qū)域復(fù)雜流動結(jié)構(gòu)進(jìn)行了分析研究;Bourdon和Dutton[2]利用平面激光成像技術(shù)對超聲速底部流場進(jìn)行研究,獲得了近尾跡區(qū)域瞬時流動的影像數(shù)據(jù)。然而風(fēng)洞試驗中對模型進(jìn)行固定的支架和相關(guān)測量設(shè)備都會對結(jié)果產(chǎn)生影響,雖然可以依靠經(jīng)驗方法對試驗結(jié)果進(jìn)行修正,但目前仍然缺乏對這些經(jīng)驗方法可靠性的驗證研究。試驗技術(shù)的局限性以及計算機(jī)技術(shù)的快速發(fā)展使得數(shù)值計算成為一條研究超聲速底部流場的便捷途徑。然而超聲速底部流場內(nèi)所包含的復(fù)雜流動現(xiàn)象對數(shù)值計算方法提出了新的挑戰(zhàn)。工程計算中常用的雷諾平均Navier-Stokes(RANS)方法對附著流和小分離流動有較好的計算能力,但由于其忽略了小尺度脈動,難以準(zhǔn)確預(yù)測超聲速底部流場出現(xiàn)的大尺度渦結(jié)構(gòu)及底部的非定常脈動現(xiàn)象。大渦模擬(LES)方法雖然可以較好地預(yù)測流動分離現(xiàn)象,但對于高雷諾數(shù)流動需要極高的網(wǎng)格分辨率和非常小的時間步長才能完全捕捉到邊界層內(nèi)擬序結(jié)構(gòu)的演化[3],過高的計算成本阻礙了LES方法的廣泛應(yīng)用。

    當(dāng)前有限的計算資源和研究者們對高精度數(shù)值計算結(jié)果的需求促使RANS/LES混合方法在近年來得到了快速的發(fā)展。研究者們對RANS和LES之間不同的混合方式開展了大量的工作,在眾多混合方法中,DES類混合方法展現(xiàn)出了廣泛應(yīng)用的潛力。針對不同尺度湍流運(yùn)動的特點,Spalart等在1997年基于Spalart-Allmara(S-A)模型提出了脫體渦模擬(DES)方法[4],也被稱為DES97方法,這一方法也成為后來眾多的DES類衍生方法的雛形。進(jìn)一步的研究表明DES方法會導(dǎo)致出現(xiàn)模型應(yīng)力耗散(MSD)現(xiàn)象[5],從而引起 網(wǎng) 格 誘 導(dǎo) 分 離 (GIS)[6]和 對 數(shù) 層 不 匹 配(LLM)[7]等問題。為了解決 MSD引起的GIS問題,Spalart提出了DES97的改進(jìn)版本——延遲DES(DDES)方法[5]。Shur等[8]在此基礎(chǔ)上又提出了改進(jìn)的DDES方法(IDDES),其目的是進(jìn)一步解決模型應(yīng)力耗散引起的對數(shù)層不匹配的問題。Xiao[9-12]和 Huang[13]等采用基于兩方程湍流模型的DDES和IDDES方法對串列圓柱[9-10]、跨聲速激波震蕩[13]、基本起落架[11]和空腔流動[12]等復(fù)雜湍流問題開展了廣泛的計算研究。

    各種RANS/LES方法也開始被國內(nèi)外學(xué)者們應(yīng)用到超聲速底部流場的研究中。Forsythe等[14]采用了基于S-A和剪切應(yīng)力輸運(yùn)(SST)湍流模型的DES方法,并引入了可壓縮修正,分別采用結(jié)構(gòu)和非結(jié)構(gòu)計算網(wǎng)格,結(jié)果表明DES方法的計算結(jié)果與實驗吻合良好。Baurle等[15]采用一種分區(qū)域的RANS/LES混合方法,其將SST模型和一方程SGS模型結(jié)合來預(yù)測超聲速底部流動。Kawai和Fujii[16]采用BL和Smagrinsky模型結(jié)合的RANS/LES方法,通過與LES、MILES和RANS計算結(jié)果的比較證實該方法在計算精度和效率上都有較好的表現(xiàn)。Simon等[17]分別采用了基于S-A湍流模型的分區(qū)域和全局DES方法,通過計算對不同方法對超聲速底部流場的預(yù)測能力進(jìn)行了對比分析。薛幫猛和楊永[18]采用了基于兩方程湍流模型的DES方法,計算結(jié)果表明,相對于RANS方法,DES方法可以更好地模擬超聲速底部流場中分離渦的發(fā)展。肖志祥和符松[19]采用了基于兩方程SST模型和弱非線性修正的k-ω模型的DES和DDES方法,引入壓縮性修正,結(jié)果表明壓縮性修正可以更好地描述流場特征,DES類混合方法能夠捕捉到豐富的流動現(xiàn)象和非定常特性。高瑞澤[20]和陳琦[21]等構(gòu)造了各自的 BL-Smagorinsky混合模型,相關(guān)計算表明這些方法能夠得到優(yōu)于RANS的模擬結(jié)果。

    本文采用基于兩方程k-ωSST湍流模型的多種RANS/LES混合方法并考慮可壓縮修正,對超聲速底部流動開展了數(shù)值模擬研究,結(jié)合不同空間離散格式,分析了數(shù)值耗散效應(yīng)對流場計算結(jié)果的影響,并對不同RANS/LES混合方法在超聲速底部流動中的表現(xiàn)進(jìn)行了對比分析,為進(jìn)一步研究超聲速底部流動減阻改進(jìn)措施和RANS/LES方法在非定常高可壓縮性流動中的應(yīng)用提供了參考。

    1 數(shù)值方法及湍流模型

    本文工作采用中心有限體積法,時間離散采用二階全隱式LU-SGS-τTS算法,空間離散采用基于五階WENO插值的Roe格式。超聲速底部流場會出現(xiàn)膨脹激波等復(fù)雜的流動問題,采用Roe格式可能出現(xiàn)非物理解或影響計算的穩(wěn)定性,為了正確求解流場,本文中采用 Harten和Hyman[22]提出的熵修正方法。

    對于超聲速底部流動的數(shù)值模擬,F(xiàn)orsythe等[14]通過比較一方程S-A模型和兩方程SST模型的計算結(jié)果,得出SST模型的表現(xiàn)明顯好于S-A模型,在本文計算中采用SST模型作為RANS/LES混合方法的基礎(chǔ)湍流模型。借鑒DES97中混合長度的思想,Travin等[7]構(gòu)建了基于SST模型的DES方法。通過引入過渡函數(shù)避免LES模式在邊界層內(nèi)被錯誤的激活,Menter和Kuntz[6]提出了基于SST模型的DDES方法。

    1.1 基于k-ωSST模型的RANS/LES混合方法

    原始的SST模型[23]對近壁區(qū)流動采用 Wilcox提出的k-ω模型,而在邊界層邊緣及自由剪切層則采用k-ε模型,兩者通過混合函數(shù)進(jìn)行過渡。采用Suzen和Hoffmann[24]提出的可壓縮修正,這種修正方法在模型中的k-ε部分加入了可壓縮耗散項和壓力擴(kuò)張項,引入可壓縮修正和混合長度l后,湍動能輸運(yùn)方程可寫為

    式中:ρ為密度;k為湍動能;σk=0.5;μ為分子運(yùn)同的混合方法取不同的值,對于引入可壓縮修正后的SST模型,有

    式中:常系數(shù)α1=1.0,α2=0.4,α3=0.2。

    混合函數(shù)

    1.1.1 DDES的構(gòu)造

    不同于 Menter等[6]提出的DDES方法,本文所采用的DDES方法中的混合長度是遵循Sparlat[5]基于 S-A 湍流模型提出的 DDES方法的形式而來的,Gritskevich等[25]對基于SST模型的DDES方法的一些參數(shù)重新進(jìn)行了標(biāo)定。當(dāng)采用DDES方法時混合長度取

    式中:CDES=0.78F1+0.61(1-F1);亞格子尺度Δmax=max(Δx,Δy,Δz),Δx、Δy、Δz分別為x、y、z3個方向上的網(wǎng)格步長,延遲函數(shù)為

    其中:νt為湍流黏性;κ為Karman常數(shù),取值為0.41;S為應(yīng)變變化率張量;Ω為旋度張量。

    1.1.2 MDDES的構(gòu)造

    Vatsa和Lockard[26]在應(yīng)用DDES方法對雙圓柱算例進(jìn)行模擬時,圓柱上游會出現(xiàn)較強(qiáng)的渦黏性,為了保持DDES方法在近壁面表現(xiàn),同時避免這種非物理的解,利用延遲函數(shù)fd構(gòu)造了一個新的混合函數(shù)來對k方程的產(chǎn)生項進(jìn)行干預(yù),這就是更改的DDES(MDDES)方法。當(dāng)fd>CMMD時(常數(shù)CMMD=0.975),修改后的k方程產(chǎn)生項為

    1.1.3 IDDES的構(gòu)造

    當(dāng)采用IDDES方法時,混合長度l取值為

    式中:亞格子長度Δ=min[Cwmax(d,Δmax),Δmax],Cw=0.15。lSST、CDES和Δmax的定義與式(4)中一致。經(jīng)驗混合函數(shù)珟fd的定義為

    當(dāng)fe=0時,IDDES以DDES模式運(yùn)行;當(dāng)fe>0,珟fd=fb時,IDDES方法在近壁區(qū)以WMLES模式運(yùn)行。函數(shù)fe的具體構(gòu)成為

    式中:各模型常數(shù)Cdt1=20,Cdt2=3,Cl=5.0,Ct=1.87。

    1.2 空間離散格式

    對無黏通量項離散采用基于迎風(fēng)型通量差分裂的Roe格式。在三維貼體正交坐標(biāo)系下穿過網(wǎng)格單元界面的無黏通量可表示為

    采用Roe格式對Navier-Stokes方程空間項進(jìn)行離散求解時,首先求得網(wǎng)格單元體交界面兩側(cè)變量的值qL和qR。對網(wǎng)格交界面左右兩側(cè)的變量是通過差值得到的。利用基于平均的網(wǎng)格中心點處的值來插值得到網(wǎng)格交界面處的變量值。因此,對交界面處變量插值的精度就決定了整個空間離散格式的精度。為了研究格式耗散對計算結(jié)果的影響,本文分別采用基于MUSCL插值的MUSCL-Roe三階格式和基于改進(jìn)的 WENO方法對交界面兩邊狀態(tài)參量進(jìn)行重構(gòu)的 WENORoe五階格式。

    由MUSCL插值方法[27]得到的三階迎風(fēng)差分格式為

    以一維標(biāo)量模型方程為例,對五階 WENO守恒差分格式進(jìn)行介紹,五階精度的邊界外推值形式為

    ωk為非線性加權(quán)因子,其定義為

    2 計算網(wǎng)格

    數(shù)值模擬是基于Herrin和Dutton[1]的實驗進(jìn)行的,入口自由來流馬赫數(shù)為2.46,入口壓力p=31 415Pa,溫度為T=145K,實驗圓柱直徑D=63.5mm,基于圓柱直徑的雷諾數(shù)Re=2.858×106。

    入口速度分布與實驗觀測值一致,圓柱表面為無滑移邊界,出口為插值邊界,遠(yuǎn)場為無反射邊界條件?;趫A柱直徑D和聲速的無量綱時間步長為Δt=0.001,對應(yīng)的物理時間步長為2.63×10-4s,每個時間步內(nèi)取15個子迭代。待流場穩(wěn)定后取最后3×104時間步的計算結(jié)果進(jìn)行統(tǒng)計平均。

    采用多塊結(jié)構(gòu)化網(wǎng)格技術(shù)生成了計算網(wǎng)格,整個計算域尺寸為4.5D×9D(圓柱區(qū)域底部半徑×圓柱高度,D為圓柱直徑),入口距圓柱底部為4D,底部下游取5D。對于圓柱底部附近及下游關(guān)注區(qū)域內(nèi)的網(wǎng)格進(jìn)行了加密,壁面首層網(wǎng)格滿足y+≤1,計算網(wǎng)格如圖1和圖2所示。為了研究網(wǎng)格尺度的影響,劃分了3套計算網(wǎng)格,均采用相同的網(wǎng)格拓?fù)浣Y(jié)構(gòu)和相同的壁面首層網(wǎng)格尺度,網(wǎng)格總數(shù)分別為6×106(Mesh-C),1.5×107(Mesh-M),2.1×107(Mesh-F)。采用 DDES方法配合三階MUSCL-Roe空間格式分別對3套網(wǎng)格求解。

    這里僅對不同密度網(wǎng)格下超聲速底部流場中心軸線上的軸向時均速度U/U∞分布進(jìn)行對比(如圖3所示),細(xì)節(jié)的流場分析見3.1節(jié)。可以看出,3套網(wǎng)格計算得到的最大回流速度均比實驗值要大。粗網(wǎng)格計算得到的最大回流速度明顯超過預(yù)測,流動再附點位于x/R=2.21,與實驗值的相對誤差為17.23%;中等密度網(wǎng)格和細(xì)網(wǎng)格計算得到的最大回流速度和流動再附位置基本一致,再附點位于x/R=2.38,與實驗值的相對誤差為10.86%。從計算精度和計算效率的角度出發(fā),本文后面的工作均選用中等密度網(wǎng)格(Mesh-M)。

    3 計算結(jié)果分析

    3.1 流 場

    為了深入研究超聲速底部的流動機(jī)理,首先對數(shù)值計算得到的流場進(jìn)行分析,計算采用IDDES方法并配合五階WENO空間格式。

    圖4所示為采用IDDES方法計算得到的z=0截面的瞬時密度梯度云圖,從圖中可以看出超聲速底部流場的整個發(fā)展過程。由于圓柱底部的幾何轉(zhuǎn)折使得湍流邊界層在此處形成分離,同時形成了較強(qiáng)的膨脹波(Ⅰ)。在靠近圓柱底部形成了一個回流區(qū)域(Ⅱ),可以看出許多小尺度的渦結(jié)構(gòu)在該區(qū)域內(nèi)運(yùn)動,在底部拐角處,流動的轉(zhuǎn)折角度和膨脹波的強(qiáng)度都由回流區(qū)域的大小決定。自由剪切層(Ⅲ)將回流區(qū)與外部的高速流動區(qū)域分割開來。在膨脹波的作用下,剪切層轉(zhuǎn)向中心軸線方向。在中心軸線上存在一個軸向速度為0的點即為流動再附點,剪切層經(jīng)過再附區(qū)域(Ⅳ)后,流動轉(zhuǎn)而沿著中心軸線發(fā)展,同時形成了再壓縮激波(Ⅴ)。在再附區(qū)域內(nèi),一部分動量不足的氣流由于逆壓梯度的影響被推到回流區(qū)內(nèi),另一部分來流則向著下游尾跡區(qū)繼續(xù)發(fā)展,從圖中可以看出在回流區(qū)內(nèi)小尺度的渦結(jié)構(gòu)較為集中,而在尾跡區(qū)內(nèi)則多是較大的湍流結(jié)構(gòu),這一點在圖5中得到了進(jìn)一步的印證。

    圖5 為IDDES方法計算得到的超聲速底部流場的瞬時湍流結(jié)構(gòu)圖,圖中展示的是依據(jù)Q準(zhǔn)則[30]的等值面圖,可以看出在圓柱底部流動具有較強(qiáng)的非定常特性,底部區(qū)域呈現(xiàn)出豐富的湍流渦結(jié)構(gòu),既包含大的發(fā)卡渦同時也包含了許多小尺度的渦結(jié)構(gòu)。也表明本文所采用的方法對小尺度渦具有較強(qiáng)的捕捉能力。

    3.2 空間離散格式的影響

    采用基于SST湍流模型的DDES方法對超聲速底部流動進(jìn)行數(shù)值模擬,空間離散格式分別采用三階 MUSCL-Roe格式和五階 WENO-Roe格式,對計算結(jié)果進(jìn)行對比分析,以研究兩種空間格式對數(shù)值結(jié)果的影響。

    圖6對比了三階 MUSCL-Roe格式和五階WENO-Roe格式在底部上游1mm位置的邊界層速度型。在邊界層內(nèi),DDES方法中的RANS模式有效,可以看出兩種格式對未分離前的邊界層的預(yù)測結(jié)果沒有差異,均與實驗觀測值吻合良好。

    圖7 對比了兩種空間格式得到的底部時均壓力系數(shù)Cp分布,從圖中可以看出格式耗散對表表面壓力分布的影響,實驗中底部表明的壓力系數(shù)分布沿徑向較為平坦,而三階MUSCL-Roe格式得到表面壓力系數(shù)在徑向變化較大;五階WENO-Roe格式得到的表面壓力系數(shù)在徑向變化則相對較小,與實驗中的Cp徑向分布趨勢一致,Cp≈-0.09,與實驗觀測值-0.102的相對誤差為11.76%,這明顯優(yōu)于三階MUSCL-Roe格式。

    對比兩種空間格式得到的中心軸線上的速度分布(如圖8所示),可以看出五階格式得到的最大回流速度和再附點位置均與實驗吻合良好,而三階MUSCL-Roe格式得到的回流區(qū)域較之實驗要小,預(yù)測的再附位置更靠近上游,最大回流速度也比實驗值要大。

    在底部下游區(qū)域取x/R=0.629 9,x/R=1.259 8,和x/R=1.889 8這3個位置的徑向和軸向速度分布進(jìn)行對比,如圖9和圖10所示。從圖中可以看出,由于較大的數(shù)值耗散,三階MUSCL-Roe格式得到的自由剪切層的寬度明顯小于五階WENO-Roe格式,自由剪切層位置較之實驗更靠近中心軸線;而五階 WENO-Roe格式預(yù)測的自由剪切層的寬度和位置與實驗吻合良好。計算結(jié)果表明,對于超聲速底部流場出現(xiàn)的自由剪切層等間斷現(xiàn)象,五階 WENO-Roe格式比三階MUSCL-Roe格式具有更強(qiáng)的處理能力。

    圖11展示了兩種空間格式計算得到的截面瞬時渦量云圖,由于存在較大的耗散,三階MUSCL-Roe格式只捕捉到了底部下游區(qū)域較大的流動結(jié)構(gòu),與之相比,五階 WENO-Roe格式能夠捕捉到更多的小尺度結(jié)構(gòu)。

    3.3 RANS/LES混合方法的影響

    分別采用基于SST湍流模型的DDES,MDDES和IDDES方法對超聲速底部流動進(jìn)行數(shù)值模擬,空間離散格式采用五階WENO-Roe格式,對計算結(jié)果進(jìn)行對比分析,以研究不同RANS/LES方法的表現(xiàn)。

    圖12比較了3種方法在底部前1mm計算得到的邊界層速度分布。對比3種混合方法的計算結(jié)果,DDES和IDDES方法計算得到的速度分布基本一致,這是由于IDDES方法在底部的上游附體區(qū)域,fd函數(shù)中的(1-fdt)分支被激活,即模型處于DDES模式。而MDDES方法計算得到的邊界層速度型與DDES和IDDES方法具有較大差異,MDDES所預(yù)測的邊界層高度較之DDES和IDDES方法要薄。與實驗值相比,DDES和IDDES計算得到的邊界層高度與實驗結(jié)果吻合較好;MDDES計算得到的邊界層厚度遠(yuǎn)小于實驗觀測值,其速度分布與實驗結(jié)果存在較大差異。

    從底部壓力系數(shù)(圖13)的對比中可以看出,3種方法計算得到的底部壓力分布均較為平坦,這一趨勢與實驗結(jié)果是一致的,較之其他兩種方法,MDDES預(yù)測值與實驗值差異最大,DDES計算得到的底部壓力系數(shù)分布略小于實驗值,ID-DES方法在這里處于WMLES模式,其預(yù)測值與實驗值吻合良好。結(jié)合圖12中邊界層速度型的計算結(jié)果,IDDES方法在近壁區(qū)的表現(xiàn)比DDES和MDDES方法要好,這是由于IDDES將RANS/LES切換向邊界層內(nèi)壓迫使其在近壁面區(qū)域的表現(xiàn)更接近實驗值。

    圖14是底部下游沿中心軸線上的軸向時均速度分布對比。與實驗相比,DDES和IDDES方法所預(yù)測的回流區(qū)最大速度位置和再附點位置都與實驗吻合良好。MDDES方法所預(yù)測的回流區(qū)域比實驗要大,最大速度位置與再附點和實驗相比更靠近下游。在再附點之后,DDES和IDDES方法的軸向速度分布與實驗存在一定差距,兩種方法在x/R=4.0附近均開始偏離實驗值,這可能是隨著與圓柱底部的距離不斷增加,下游尾跡區(qū)的網(wǎng)格尺度逐漸增大所引起的。

    圖15 和圖16分別比較了圓柱底部下游區(qū)域3個不同位置的軸向和徑向時均速度。在x/R=0.629 9和x/R=1.259 8位置,3種方法計算得到的剪切層均位于r/R=0.8~1.0區(qū)域內(nèi),均與實驗結(jié)果吻合良好。在剪切層向里的回流區(qū)內(nèi)(r/R<0.8),DDES和IDDES的計算結(jié)果較為接近,也更靠近實驗值,相對而言MDDES的計算結(jié)果較實驗有一定的差異。進(jìn)一步向下游發(fā)展,在x/R=1.889 8位置,此處靠近再附區(qū)域,DDES和IDDES方法得到的剪切層高度位于r/R=0.5附近,與實驗值吻合良好,而MDDES方法的預(yù)測值比實驗要高??傮w來看,DDES和IDDES所預(yù)測的不同位置軸向和徑向速度分布更接近實驗值。

    圖17 比較了3種不同混合方法計算得到的瞬時渦量云圖。從圖中可以看出,3種方法均能夠捕捉到底部流場的小尺度渦結(jié)構(gòu),相比之下IDDES所預(yù)測下游流場渦結(jié)構(gòu)更加細(xì)膩。

    4 結(jié) 論

    考慮可壓縮修正,分別采用基于兩方程k-ω SST湍流模型的DDES、MDDES和IDDES方法對超聲速底部流場進(jìn)行數(shù)值模擬研究。計算結(jié)果表明RANS/LES方法能夠捕捉到超聲速底部流場中豐富的渦結(jié)構(gòu),深入認(rèn)識超聲速底部流場的流動機(jī)理。采用三階MUSCL-Roe和五階 WENO-Roe空間離散格式研究了格式耗散對流場計算的影響,并對不同RANS/LES混合方法的表現(xiàn)進(jìn)行了分析。

    1)對比不同空間離散格式的影響,兩種格式都能夠預(yù)測未分離前的邊界層;由于數(shù)值耗散的影響,五階 WENO-Roe格式對圓柱底部表面壓力分布、回流區(qū)最大速度、分離再附位置、自由剪切層寬度和位置的預(yù)測均比三階MUSCL-Roe格式要好;在底部下游區(qū)域,五階 WENO-Roe格式能夠捕捉到更多細(xì)小的渦結(jié)構(gòu)。

    2)對比不同RANS/LES方法的影響,DDES和IDDES方法能夠準(zhǔn)確預(yù)測邊界層的速度分布、回流區(qū)最大速度及其位置和流動再附點位置;IDDES方法在近壁面區(qū)域有較好的表現(xiàn),能夠準(zhǔn)確預(yù)測超聲速底部的表面壓力分布;DDES和IDDES方法對不同位置速度分布的預(yù)測結(jié)果比MDDES方法要好;IDDES方法預(yù)測的下游流場渦結(jié)構(gòu)更加細(xì)膩。

    [1] HERRIN J L,DUTTON J C.Supersonic base flow experiments in the near wake of a cylindrical afterbody[J].AIAA Journal,1994,32(1):77-83.

    [2] BOURDON C J,DUTTON J C.Planar visualizations of large-scale turbulent structures in axisymmetric supersonic separated flows[J].Physics of Fluids,1999,11(1):201-213.

    [3] CHOI H,MOIN P.Grid-point requirements for large eddy simulation:Chapman’s estimates revisited[J].Physics of Fluids,2012,24(1):011702.

    [4] SPALART P R,JOU W H,STRELETS M,et al.Comments on the feasibility of LES for wings,and on a hybrid RANS/LES approach[C]/Proceeding of the First AFOSR International Conference on DNS/LES.Columbus:Greyden,1997:137-147.

    [5] SPALART P R,DECK S,SHUR M,et al.A new version of detached-eddy simulation,resistant to ambiguous grid densities[J].Theoretical and Computational Fluid Dynamics,2006,20(3):181-195.

    [6] MENTER F R,KUNTZ M.Adaptation of eddy-viscosity turbulence models to unsteady separated flow behind vehicles[C]/The Aerodynamics of Heavy Vehicles:Trucks,Buses and Trains.Berlin Heidelberg:Springer-Verlag,2004:339-352.

    [7] TRAVIN A,SHUR M,STRELETS M,et al.Physical and numerical upgrades in the detached-eddy simulation of complex turbulent flows[C]/Advances in LES of Complex Flows.Norwell:Kluwer Academic,2006:239-254.

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

    [9] XIAO Z,LIU J,HUANG J,et al.Numerical dissipation effects on massive separation around tandem cylinders[J].AIAA Journal,2012,50(5):1119-1136.

    [10] XIAO Z,LUO K.Improved delayed detached-eddy simulation of massive separation around triple cylinders[J].Acta Mechanica Sinica,2015,31(6):799-816.

    [11] XIAO Z,LIU J,LUO K,et al.Investigation of flows around a rudimentary landing gear with advanced detached-eddy-simulation approaches[J].AIAA Journal,2013,51(1):107-125.

    [12] XIAO L,XIAO Z,DUAN Z,et al.Improved-delayed-detached-eddy simulation of cavity-induced transition in hypersonic boundary layer[J].International Journal of Heat and Fluid Flow,2015,51:138-150.

    [13] HUANG J,XIAO Z,LIU J,et al.Simulation of shock wave buffet and its suppression on an OAT15Asupercritical airfoil by IDDES[J].Science China-Physics Mechanics& Astronomy,2012,55(2):260-271.

    [14] FORSYTHE J R,HOFFMANN K A,CUMMINGS R M,et al.Detached-eddy simulation with compressibility corrections applied to a supersonic axisymmetric base flow[J].Journal of Fluids Engineering,2002,124(4):911-923.

    [15] BAURLE R A,TAM C J,EDWARDS J R,et al.Hybrid simulation approach for cavity flows:blending,algrithm and boundary treatment issues[J].AIAA Journal,2003,41(8):1463-1480.

    [16] KAWAI S,F(xiàn)UJII K.Computational study of supersonic base flow using hybrid turbulence methodology[J].AIAA Journal,2005,43(6):1256-1275.

    [17] SIMON F,DECK S,GUILLEN P,et al.Reynolds-averaged navier-stokes/large-eddy simulations of supersonic base flow[J].AIAA Journal,1994,44(11):2578-2590.

    [18] 薛幫猛,楊永.基于兩方程湍流模型的DES方法在超音速圓柱底部流動計算中的應(yīng)用[J].西北工業(yè)大學(xué)學(xué)報,2006,24(5):544-547.XUE B M,YANG Y.Technical details in applying DES method to computing supersonic cylinder-base flow[J].Journal of Northwestern Polytechnical University,2006,24(5):544-547(in Chinese).

    [19] 肖志祥,符松.用RANS/LES混合方法研究超聲速底部流動[J].計算物理,2009,26(2):221-230.XIAO Z X,F(xiàn)U S.Study on supersonic base flow using RANS/LES methods[J].Chinese Journal of Computational Physics,2009,26(2):221-230(in Chinese).

    [20] 高瑞澤,閻超.LES/RANS混合方法對超聲速底部流動的應(yīng)用[J].北京航空航天大學(xué)學(xué)報,2011,37(9):1095-1099.GAO R Z,YAN C.LES/RANS hybrid method for supersonic axisymmetric base flow[J].Journal of University of Aeronautics and Astronautics,2011,37(9):1095-1099(in Chinese).

    [21] 陳琦,司芳芳,陳堅強(qiáng),等.RANS/LES在超聲速突起物繞流 中 的 應(yīng) 用 研 究 [J].航 空 學(xué) 報,2013,34(7):1531-1537.CHEN Q,SI F F,CHEN J Q,et al.Study of protuberances in supersonic flow with RANS/LES method[J].Acta Aeronautica et Astronautica Sinica,2013,34(7):1531-1537(in Chinese).

    [22] HARTEN A,HYMAN J M.Self-adjusting grid methods for one-dimensional hyperbolic conservation laws[J].Journal of Computational Physics,1983,50(2):235-369.

    [23] MENTER F R.Two-equation eddy viscosity turbulence models for engineering applications[J].AIAA Journal,1994,32(8):1598-1605.

    [24] SUZEN Y B,HOFFMANN K A.Investigation of supersonic jet exhaust flow by one-and two-equation turbulence models:AIAA-1998-0322[R].Reston:AIAA,1998.

    [25] GRITSKEVICH M S,GARBARUK A V,SCHUTZE J,et al.Development of DDES and IDDES formulations for the k-ω shear stress transport model[J].Flow,Turbulence and Combustion,2012,88(3):431-449.

    [26] VATSA V N,LOCKARD D P.Assessment of hybrid rans/les turbulence model for aeroacoustics applications:AIAA-2010-4011[R].Reston:AIAA,2010.

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

    [28] BORGES R,CARMONA M,COSTA B,et al.An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[J].Journal of Computational Physics,2008,227(6):3191-3211.

    [29] JIANG G S,SHU C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics,1996,126(1):202-228.

    [30] JEONG J,HUSSAIN F.On the identification of a vortex[J].Journal of Fluid Mechanics,1995,285(2):69-94.

    Numerical simulations of supersonic base flow field based on RANS/LES approaches

    ZHANG Lu,LI Jie*
    School of Aeronautics,Northwestern Polytechnical University,Xi’an 710072,China

    Numerical investigation on supersonic base flow is performed using several Reynolds averaged Navier-Stokes(RANS)and large eddy simulation(LES)hybrid methods based on k-ωshear stress transport(SST)model with compressibility correction,including delayed detached eddy simulation(DDES),modified-DDES(MDDES)and improved-DDES(IDDES)approaches.Third-order MUSCL-Roe and fifth-order WENO-Roe spatial scheme are applied in the investigation.The numerical results show plenty of small scale turbulence structure in supersonic base flow.The complex flow physics are comprehensively understood,which provides references for the base aerodynamic drag reduction and the application of RANS/LES hybrid methods to unsteady highly compressible flows in future research.Numerical dissipation effects of two spatial schemes are investigated.Computational results show that fifth-order WENO-Roe scheme is more validated than third-order MUSCL-Roe scheme when compared with experimental data.Furthermore,comparative analysis of the computational results with several RANS/LES hybrid methods is conducted.The results show that IDDES approach has better performance in regions near the wall than DDES and MDDES approaches.

    RANS/LES hybrid methods;supersonic base flow;compressibility effects;turbulence model;numerical simula-tions

    2016-01-22;Revised:2016-02-18;Accepted:2016-05-12;Published online:2016-05-12 12:48

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

    s:National Basic Research Program of China(2015CB755800);National Natural Science Foundation of China(11172240);Aeronautical Science Foundation of China(2014ZA53002)

    V211.3

    A

    1000-6893(2017)01-120102-12

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

    10.7527/S1000-6893.2016.0145

    2016-01-22;退修日期:2016-02-18;錄用日期:2016-05-12;網(wǎng)絡(luò)出版時間:2016-05-12 12:48

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

    國家重點基礎(chǔ)研究發(fā)展計劃 (2015CB755800);國家自然科學(xué)基金 (11172240);航空科學(xué)基金 (2014ZA53002)

    *通訊作者 .E-mail:lijieruihao@163.com

    張露,李杰.基于RANS/LES方法的超聲速底部流場數(shù)值模擬[J].航空學(xué)報,2017,38(1):120102.ZHANG L,LI J.Numerical simulations of supersonic base flow field based on RANS/LES approaches[J].Acta Aeronautica et Astronautica Sinica,2017,38(1):120102.

    (責(zé)任編輯:李明敏)

    *Corresponding author.E-mail:lijieruihao@163.com

    猜你喜歡
    超聲速湍流流場
    高超聲速出版工程
    高超聲速飛行器
    大型空冷汽輪發(fā)電機(jī)轉(zhuǎn)子三維流場計算
    重氣瞬時泄漏擴(kuò)散的湍流模型驗證
    轉(zhuǎn)杯紡排雜區(qū)流場與排雜性能
    超聲速旅行
    基于HYCOM的斯里蘭卡南部海域溫、鹽、流場統(tǒng)計分析
    基于瞬態(tài)流場計算的滑動軸承靜平衡位置求解
    高超聲速大博弈
    太空探索(2014年5期)2014-07-12 09:53:28
    “青春期”湍流中的智慧引渡(三)
    欧美成人一区二区免费高清观看| 80岁老熟妇乱子伦牲交| 亚洲国产日韩欧美精品在线观看| 国产女主播在线喷水免费视频网站 | 一个人免费在线观看电影| 狂野欧美白嫩少妇大欣赏| 日韩欧美国产在线观看| 全区人妻精品视频| 最新中文字幕久久久久| 69人妻影院| 亚洲成人久久爱视频| 亚洲在线自拍视频| av.在线天堂| 日本wwww免费看| 综合色丁香网| 中文字幕亚洲精品专区| 女人久久www免费人成看片| 亚洲欧美日韩无卡精品| 激情 狠狠 欧美| 人妻少妇偷人精品九色| 亚洲av中文字字幕乱码综合| 日韩av在线免费看完整版不卡| 国产精品一区二区三区四区久久| 精品一区二区三区视频在线| 国产伦精品一区二区三区四那| 欧美xxⅹ黑人| 男女边吃奶边做爰视频| 日韩 亚洲 欧美在线| 精品国内亚洲2022精品成人| 久久热精品热| 综合色av麻豆| 三级国产精品片| 日韩欧美一区视频在线观看 | 亚洲精品第二区| 秋霞伦理黄片| 亚洲精品aⅴ在线观看| 黑人高潮一二区| 午夜免费激情av| 亚洲国产精品国产精品| 亚洲电影在线观看av| 精品一区二区免费观看| 禁无遮挡网站| 亚洲在线观看片| 99视频精品全部免费 在线| 草草在线视频免费看| 免费观看无遮挡的男女| 免费黄网站久久成人精品| 成人高潮视频无遮挡免费网站| 精品人妻视频免费看| 偷拍熟女少妇极品色| 三级国产精品片| 国产亚洲精品av在线| 神马国产精品三级电影在线观看| 久久99热这里只频精品6学生| 久久精品久久久久久久性| 尾随美女入室| 少妇丰满av| 国产av国产精品国产| 婷婷色av中文字幕| 白带黄色成豆腐渣| 黑人高潮一二区| 五月天丁香电影| 男女那种视频在线观看| 永久免费av网站大全| xxx大片免费视频| 大香蕉97超碰在线| 国产在线男女| 国产精品无大码| 少妇裸体淫交视频免费看高清| 国产69精品久久久久777片| 国产伦在线观看视频一区| av播播在线观看一区| 搡老乐熟女国产| 久久精品国产亚洲网站| 亚洲色图av天堂| 精品国产三级普通话版| 国内精品宾馆在线| 成人午夜精彩视频在线观看| 亚洲激情五月婷婷啪啪| 成人漫画全彩无遮挡| 建设人人有责人人尽责人人享有的 | 日本av手机在线免费观看| 国产精品伦人一区二区| 欧美日韩综合久久久久久| 我的女老师完整版在线观看| 久久精品夜色国产| 国产伦在线观看视频一区| 97人妻精品一区二区三区麻豆| 一二三四中文在线观看免费高清| 观看免费一级毛片| 日韩一本色道免费dvd| 国国产精品蜜臀av免费| 特大巨黑吊av在线直播| 日本免费在线观看一区| 91久久精品国产一区二区三区| 久久精品国产亚洲网站| 欧美不卡视频在线免费观看| 亚洲av福利一区| 成人性生交大片免费视频hd| 不卡视频在线观看欧美| 日本wwww免费看| 亚洲va在线va天堂va国产| 在线观看人妻少妇| 男女边摸边吃奶| 我要看日韩黄色一级片| 亚洲一区高清亚洲精品| 国产精品一区二区三区四区免费观看| 国产成年人精品一区二区| 成年人午夜在线观看视频 | 丰满人妻一区二区三区视频av| 日本欧美国产在线视频| 波多野结衣巨乳人妻| 国产爱豆传媒在线观看| 最近中文字幕高清免费大全6| 晚上一个人看的免费电影| av天堂中文字幕网| 日韩av免费高清视频| 在线免费十八禁| 国产精品人妻久久久久久| 久久99蜜桃精品久久| 国产免费视频播放在线视频 | 国产成人福利小说| 极品教师在线视频| 国产精品一区二区性色av| 日本猛色少妇xxxxx猛交久久| 熟女人妻精品中文字幕| 精品久久久久久成人av| 少妇高潮的动态图| 亚洲欧洲国产日韩| 99九九线精品视频在线观看视频| 亚洲图色成人| 国产免费福利视频在线观看| 免费大片18禁| 国产色爽女视频免费观看| 午夜福利网站1000一区二区三区| 少妇裸体淫交视频免费看高清| 丝瓜视频免费看黄片| 免费观看在线日韩| 淫秽高清视频在线观看| 国产在线一区二区三区精| 亚洲电影在线观看av| 欧美另类一区| 狠狠精品人妻久久久久久综合| 色网站视频免费| 国产男人的电影天堂91| 免费黄色在线免费观看| av天堂中文字幕网| 我要看日韩黄色一级片| 日韩欧美国产在线观看| 91aial.com中文字幕在线观看| 一二三四中文在线观看免费高清| 午夜免费观看性视频| 2021天堂中文幕一二区在线观| 久久6这里有精品| 尤物成人国产欧美一区二区三区| 午夜福利在线在线| 18禁裸乳无遮挡免费网站照片| 亚洲三级黄色毛片| 五月伊人婷婷丁香| 简卡轻食公司| 黄色日韩在线| 亚洲最大成人av| 亚洲国产成人一精品久久久| 国内少妇人妻偷人精品xxx网站| 免费看美女性在线毛片视频| 亚洲天堂国产精品一区在线| 亚洲真实伦在线观看| 亚洲国产精品成人综合色| 免费黄频网站在线观看国产| 国产三级在线视频| 在现免费观看毛片| 日本熟妇午夜| 国产精品人妻久久久影院| 亚洲精品第二区| 看十八女毛片水多多多| 我要看日韩黄色一级片| 国产成人精品一,二区| 神马国产精品三级电影在线观看| videossex国产| a级一级毛片免费在线观看| 少妇被粗大猛烈的视频| 欧美xxxx性猛交bbbb| 日韩,欧美,国产一区二区三区| 国产精品一及| 亚洲一区高清亚洲精品| 亚洲第一区二区三区不卡| 中国国产av一级| 欧美成人午夜免费资源| 91久久精品电影网| 精品久久久久久久末码| av在线亚洲专区| www.色视频.com| 国产精品美女特级片免费视频播放器| 国产精品一区二区在线观看99 | 国产在线一区二区三区精| 国产永久视频网站| 欧美成人一区二区免费高清观看| 熟妇人妻久久中文字幕3abv| 一级毛片 在线播放| 99热全是精品| 亚洲成色77777| 激情五月婷婷亚洲| 亚洲图色成人| 美女大奶头视频| 免费观看的影片在线观看| 亚洲性久久影院| 99久久中文字幕三级久久日本| 青春草国产在线视频| 欧美最新免费一区二区三区| 校园人妻丝袜中文字幕| 青春草视频在线免费观看| 精品人妻视频免费看| 毛片女人毛片| 国产精品.久久久| 亚洲国产欧美人成| 亚洲精品一区蜜桃| 免费高清在线观看视频在线观看| 亚洲精品日本国产第一区| 亚洲精华国产精华液的使用体验| xxx大片免费视频| 国产精品久久久久久久久免| 免费观看av网站的网址| 卡戴珊不雅视频在线播放| 日韩电影二区| 我的女老师完整版在线观看| 晚上一个人看的免费电影| 午夜福利网站1000一区二区三区| 91久久精品国产一区二区成人| 国产淫语在线视频| 免费电影在线观看免费观看| 蜜桃亚洲精品一区二区三区| 亚洲精品亚洲一区二区| 99久国产av精品国产电影| 亚洲无线观看免费| 久久久国产一区二区| 18禁动态无遮挡网站| 亚洲欧美一区二区三区黑人 | 亚洲av.av天堂| 一级毛片黄色毛片免费观看视频| 日韩中字成人| 少妇猛男粗大的猛烈进出视频 | 91精品伊人久久大香线蕉| 国产亚洲av嫩草精品影院| 国产精品美女特级片免费视频播放器| 极品教师在线视频| 黄色欧美视频在线观看| 亚洲,欧美,日韩| 精品国产露脸久久av麻豆 | 亚洲自偷自拍三级| 久久久精品欧美日韩精品| 卡戴珊不雅视频在线播放| 能在线免费观看的黄片| 国产精品蜜桃在线观看| 欧美丝袜亚洲另类| 精品久久久久久久人妻蜜臀av| 黄色一级大片看看| 身体一侧抽搐| 亚洲精品日韩av片在线观看| av免费在线看不卡| 国产精品精品国产色婷婷| 欧美日韩国产mv在线观看视频 | 国产极品天堂在线| 久久久久精品久久久久真实原创| 久久人人爽人人爽人人片va| 一夜夜www| 在线免费十八禁| 别揉我奶头 嗯啊视频| 久久久欧美国产精品| 国产男女超爽视频在线观看| 国产亚洲5aaaaa淫片| 日日撸夜夜添| 欧美性猛交╳xxx乱大交人| 婷婷色综合www| 蜜桃亚洲精品一区二区三区| 成年女人在线观看亚洲视频 | 99久久九九国产精品国产免费| 国产永久视频网站| 美女国产视频在线观看| 青春草亚洲视频在线观看| 18禁裸乳无遮挡免费网站照片| 亚洲成人一二三区av| 日韩视频在线欧美| 亚洲经典国产精华液单| 久久久午夜欧美精品| 日韩三级伦理在线观看| 欧美激情在线99| 欧美xxxx性猛交bbbb| 久久久a久久爽久久v久久| 不卡视频在线观看欧美| 国产乱人视频| 亚洲av.av天堂| 日本免费a在线| 能在线免费观看的黄片| 尤物成人国产欧美一区二区三区| 99re6热这里在线精品视频| 22中文网久久字幕| 爱豆传媒免费全集在线观看| 在线免费十八禁| 成人美女网站在线观看视频| 三级国产精品欧美在线观看| 一级毛片电影观看| 国产精品蜜桃在线观看| 日本免费在线观看一区| 欧美性感艳星| 久久久久久久久久久丰满| av在线蜜桃| 深爱激情五月婷婷| 男女那种视频在线观看| 亚洲乱码一区二区免费版| 一区二区三区乱码不卡18| 伊人久久精品亚洲午夜| 欧美日韩在线观看h| 春色校园在线视频观看| 天堂av国产一区二区熟女人妻| 少妇丰满av| 亚洲欧美中文字幕日韩二区| 九草在线视频观看| 国产av不卡久久| 亚洲精品456在线播放app| 夜夜看夜夜爽夜夜摸| 欧美性猛交╳xxx乱大交人| 国产在视频线在精品| 国产精品国产三级专区第一集| 美女内射精品一级片tv| 久久这里有精品视频免费| eeuss影院久久| 国产人妻一区二区三区在| 精品酒店卫生间| 精品久久国产蜜桃| 高清毛片免费看| 中文字幕av成人在线电影| 国产精品久久久久久av不卡| 午夜视频国产福利| 免费播放大片免费观看视频在线观看| 亚洲欧美一区二区三区国产| 久久久久久久久久久免费av| 国产精品爽爽va在线观看网站| 免费黄网站久久成人精品| 欧美三级亚洲精品| 免费大片18禁| 国产成人a区在线观看| 男女边摸边吃奶| 少妇裸体淫交视频免费看高清| 国产精品蜜桃在线观看| 一级毛片黄色毛片免费观看视频| 国产成人福利小说| 国产v大片淫在线免费观看| 人人妻人人澡人人爽人人夜夜 | 在现免费观看毛片| 一个人看视频在线观看www免费| 国产成年人精品一区二区| 久久久久性生活片| 最近中文字幕2019免费版| 丰满乱子伦码专区| 99久久人妻综合| 日韩一本色道免费dvd| 六月丁香七月| 精品久久久久久久人妻蜜臀av| 草草在线视频免费看| 日韩 亚洲 欧美在线| 国产精品1区2区在线观看.| 亚洲欧洲日产国产| 尤物成人国产欧美一区二区三区| 日日啪夜夜撸| 日本一二三区视频观看| 中文资源天堂在线| 精品国产露脸久久av麻豆 | 欧美激情国产日韩精品一区| 校园人妻丝袜中文字幕| 日韩国内少妇激情av| 青青草视频在线视频观看| 国产精品不卡视频一区二区| 欧美日韩视频高清一区二区三区二| 国产午夜精品论理片| 偷拍熟女少妇极品色| 亚洲,欧美,日韩| 国产黄色小视频在线观看| 久久久欧美国产精品| 国产精品爽爽va在线观看网站| 免费看不卡的av| 淫秽高清视频在线观看| 欧美+日韩+精品| 天美传媒精品一区二区| 成年版毛片免费区| 亚洲va在线va天堂va国产| 高清午夜精品一区二区三区| 婷婷色综合www| 午夜福利在线观看免费完整高清在| 久久国内精品自在自线图片| 麻豆久久精品国产亚洲av| 久久久色成人| 22中文网久久字幕| 久久久久性生活片| 婷婷色麻豆天堂久久| 国产乱来视频区| 精品人妻偷拍中文字幕| 韩国av在线不卡| 亚洲av国产av综合av卡| 99热这里只有是精品在线观看| 只有这里有精品99| 三级国产精品片| 国产一区二区三区av在线| 超碰97精品在线观看| 五月伊人婷婷丁香| 又爽又黄a免费视频| 国产精品国产三级国产av玫瑰| 十八禁国产超污无遮挡网站| 夫妻午夜视频| 国产色爽女视频免费观看| 久久人人爽人人爽人人片va| 亚洲欧美成人精品一区二区| 两个人视频免费观看高清| 欧美激情久久久久久爽电影| 午夜福利在线观看吧| 亚洲电影在线观看av| 精品一区二区三区视频在线| 男女边摸边吃奶| 天堂俺去俺来也www色官网 | 国产精品美女特级片免费视频播放器| 自拍偷自拍亚洲精品老妇| 麻豆精品久久久久久蜜桃| 91av网一区二区| 亚洲av不卡在线观看| 日韩不卡一区二区三区视频在线| 99热这里只有是精品50| 综合色av麻豆| 欧美激情久久久久久爽电影| 最近手机中文字幕大全| 边亲边吃奶的免费视频| 麻豆av噜噜一区二区三区| 看非洲黑人一级黄片| 尤物成人国产欧美一区二区三区| 亚洲av成人av| 91狼人影院| 亚洲欧美日韩卡通动漫| 免费无遮挡裸体视频| 毛片一级片免费看久久久久| 欧美日韩一区二区视频在线观看视频在线 | 国国产精品蜜臀av免费| 熟妇人妻久久中文字幕3abv| videossex国产| av福利片在线观看| 国产精品嫩草影院av在线观看| 国产高潮美女av| 一级黄片播放器| 亚洲国产精品sss在线观看| 中文字幕免费在线视频6| 亚洲国产av新网站| 亚洲国产成人一精品久久久| 亚州av有码| 欧美一级a爱片免费观看看| 丰满乱子伦码专区| 日韩国内少妇激情av| 午夜福利在线观看吧| 三级国产精品欧美在线观看| 成年免费大片在线观看| 久久99热这里只有精品18| 天堂影院成人在线观看| 国产黄a三级三级三级人| 欧美高清成人免费视频www| 视频中文字幕在线观看| 国产精品.久久久| av免费在线看不卡| 欧美 日韩 精品 国产| 久久久久久伊人网av| 亚洲丝袜综合中文字幕| 免费大片18禁| 人体艺术视频欧美日本| 乱系列少妇在线播放| 日产精品乱码卡一卡2卡三| 亚洲av电影不卡..在线观看| 五月玫瑰六月丁香| 国产午夜精品论理片| 性插视频无遮挡在线免费观看| 欧美日韩一区二区视频在线观看视频在线 | 欧美激情在线99| av.在线天堂| 99热6这里只有精品| 黄片wwwwww| 99热这里只有是精品50| 18禁在线播放成人免费| 最近中文字幕2019免费版| 久99久视频精品免费| 街头女战士在线观看网站| 国产在线一区二区三区精| 日韩av免费高清视频| 内射极品少妇av片p| 性色avwww在线观看| 天天躁夜夜躁狠狠久久av| 超碰97精品在线观看| 国产成人免费观看mmmm| 欧美高清成人免费视频www| 在线观看免费高清a一片| 三级国产精品欧美在线观看| 一级av片app| 99热网站在线观看| 亚洲av中文字字幕乱码综合| 色播亚洲综合网| 国产女主播在线喷水免费视频网站 | 精品一区二区三区人妻视频| 女人被狂操c到高潮| 一级毛片 在线播放| 欧美变态另类bdsm刘玥| 国产精品久久久久久av不卡| 国内少妇人妻偷人精品xxx网站| 99久国产av精品| 亚洲一级一片aⅴ在线观看| 午夜免费男女啪啪视频观看| 亚洲国产精品成人综合色| 成人毛片60女人毛片免费| 一级av片app| 久久久久久久久久黄片| av在线亚洲专区| 国语对白做爰xxxⅹ性视频网站| 男人舔奶头视频| 蜜桃久久精品国产亚洲av| 国产亚洲91精品色在线| 亚洲最大成人中文| 国产老妇伦熟女老妇高清| 亚洲成人精品中文字幕电影| 男女那种视频在线观看| 国产综合精华液| 夜夜看夜夜爽夜夜摸| 天堂俺去俺来也www色官网 | 免费黄网站久久成人精品| 国产黄a三级三级三级人| 九草在线视频观看| 国产老妇伦熟女老妇高清| 亚洲成人精品中文字幕电影| 国产成人a区在线观看| 欧美日韩视频高清一区二区三区二| 夜夜看夜夜爽夜夜摸| 麻豆av噜噜一区二区三区| 日韩大片免费观看网站| av在线亚洲专区| 高清欧美精品videossex| 街头女战士在线观看网站| 最近最新中文字幕大全电影3| 狂野欧美激情性xxxx在线观看| 久久久久久久久久人人人人人人| 免费在线观看成人毛片| 丰满人妻一区二区三区视频av| 色5月婷婷丁香| 一区二区三区乱码不卡18| 国产三级在线视频| 一级爰片在线观看| av在线播放精品| 午夜福利网站1000一区二区三区| 亚洲精品国产av蜜桃| 夫妻午夜视频| 免费av不卡在线播放| 秋霞在线观看毛片| 国产 亚洲一区二区三区 | 中文字幕av在线有码专区| 91精品伊人久久大香线蕉| 男人舔奶头视频| 日韩欧美国产在线观看| 国产综合精华液| 三级国产精品欧美在线观看| 国产91av在线免费观看| 亚洲精品国产成人久久av| 精品国产露脸久久av麻豆 | 日韩av在线免费看完整版不卡| 久久久久久久久久久丰满| 听说在线观看完整版免费高清| 肉色欧美久久久久久久蜜桃 | 国产一区二区三区综合在线观看 | 国产在视频线精品| 日韩欧美精品v在线| 国产日韩欧美在线精品| 欧美zozozo另类| 26uuu在线亚洲综合色| 亚洲丝袜综合中文字幕| 大片免费播放器 马上看| 免费高清在线观看视频在线观看| 免费黄频网站在线观看国产| 一夜夜www| ponron亚洲| 欧美成人午夜免费资源| 乱人视频在线观看| 寂寞人妻少妇视频99o| 中文字幕免费在线视频6| 精品一区在线观看国产| av.在线天堂| 成人亚洲精品av一区二区| 亚洲欧洲日产国产| 欧美区成人在线视频| 人妻系列 视频| 日韩欧美三级三区| 免费看光身美女| 国产麻豆成人av免费视频| 国产 亚洲一区二区三区 | 国产伦在线观看视频一区| 中文乱码字字幕精品一区二区三区 | 国产黄色视频一区二区在线观看| 男女下面进入的视频免费午夜| 久久久久久久久中文| 国产亚洲精品久久久com| 哪个播放器可以免费观看大片| 国产色爽女视频免费观看| 91久久精品国产一区二区三区| www.av在线官网国产| 国产人妻一区二区三区在| 高清日韩中文字幕在线| 亚洲人成网站在线观看播放| 日韩av在线大香蕉| 我的老师免费观看完整版| 熟妇人妻不卡中文字幕| 搡老乐熟女国产| av在线蜜桃| 又爽又黄a免费视频| 国产高清有码在线观看视频| 亚洲精品乱码久久久久久按摩| 日韩欧美国产在线观看| 两个人视频免费观看高清|