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

    基于IDDES框架的γ -Reθ轉(zhuǎn)捩模型

    2019-09-11 08:18:46易淼榮趙慧勇樂嘉陵肖保國鄭忠華
    航空學(xué)報 2019年8期
    關(guān)鍵詞:來流邊界層湍流

    易淼榮,趙慧勇,樂嘉陵,肖保國,鄭忠華

    中國空氣動力研究與發(fā)展中心 高超聲速沖壓發(fā)動機(jī)技術(shù)重點實驗室,綿陽 621000

    邊界層轉(zhuǎn)捩對高超聲速飛行器的設(shè)計有利有弊。一方面,邊界層轉(zhuǎn)捩導(dǎo)致氣動熱的急劇增加,嚴(yán)重影響熱防護(hù)系統(tǒng)的設(shè)計,同時轉(zhuǎn)捩還影響飛行器的摩阻、壓阻、飛行器的性能和氣動力控制;另一方面,在高超聲速進(jìn)氣道上實現(xiàn)邊界層強(qiáng)制轉(zhuǎn)捩,可以減少流動分離,提高流動抗逆壓的能力,便于進(jìn)氣道的起動和超燃發(fā)動機(jī)的正常工作。在長達(dá)一個多世紀(jì)的研究后,人們對于邊界層轉(zhuǎn)捩問題的研究取得了長足的進(jìn)步,發(fā)展了線性穩(wěn)定性、弱非線性穩(wěn)定性、瞬態(tài)增長、三波共振、二次失穩(wěn)等理論[1]描述邊界層轉(zhuǎn)捩的過程,發(fā)展了靜音風(fēng)洞來模擬天上真實的來流擾動條件[2]。但是,針對一個給定的較復(fù)雜構(gòu)型,人們?nèi)圆荒艽_切地給出邊界層會在什么條件下、什么位置由層流轉(zhuǎn)捩為湍流,特別是在高超聲速條件下,由于受到太多因素及各種因素之間互相耦合作用的影響,轉(zhuǎn)捩預(yù)測變得更加困難。

    目前能夠進(jìn)行邊界層轉(zhuǎn)捩研究的數(shù)值方法大致可分為4種:第1種是通過精細(xì)的空間和時間分辨率來解析邊界層發(fā)展過程中擾動的物理發(fā)展過程,從而捕獲邊界層由層流轉(zhuǎn)捩為湍流的詳細(xì)信息,主要為DNS(Direct Numerical Simulation)和LES(Large Eddy Simulation)方法,被用來研究典型狀態(tài)下的轉(zhuǎn)捩機(jī)制和詳細(xì)的轉(zhuǎn)捩過程[3-4],但是由于該方法對網(wǎng)格、格式和計算量的要求極大,很難針對高雷諾數(shù)、復(fù)雜構(gòu)型開展研究。第2種是通過線性和非線性穩(wěn)定性理論,通過對擾動增長率的研究從而獲得邊界層內(nèi)擾動波的增長直至失穩(wěn)過程,如LST(Linear Stability Theory)和PSE(Parabolized Stability Equation)方法[5-6],其有明確的物理機(jī)制,能夠模擬非穩(wěn)定波線性(甚至是弱非線性階段)的發(fā)展,但無法模擬Breakdown和湍斑生成的過程,也無法模擬Bypass轉(zhuǎn)捩。第3種是通過對大量試驗結(jié)果進(jìn)行總結(jié)分析,提取出關(guān)鍵的流場參數(shù),給定一個標(biāo)準(zhǔn),當(dāng)參數(shù)達(dá)到該標(biāo)準(zhǔn)時,則認(rèn)為轉(zhuǎn)捩發(fā)生,稱之為轉(zhuǎn)捩準(zhǔn)則,如NASP(National AeroSpace Plane)的尖前緣平板轉(zhuǎn)捩準(zhǔn)則:Reθ/Mae=305[7],其中Reθ為以動量厚度為長度尺度定義的雷諾數(shù),Mae為邊界層外緣馬赫數(shù)。該方法形式簡單,操作方便,但是通用性不好,只能針對類似的構(gòu)型在類似的條件下開展轉(zhuǎn)捩預(yù)測,而針對其他情況的預(yù)測就會出現(xiàn)偏差。

    由于粗糙顆粒明顯改變了邊界層內(nèi)流場的結(jié)構(gòu),其誘導(dǎo)的尾跡渦具有G?rtler渦的性質(zhì),這種反向旋轉(zhuǎn)的流向渦對導(dǎo)致邊界層內(nèi)的速度剖面出現(xiàn)大的變形,使得流向速度在法向和展向均存在拐點。G?rtler渦結(jié)構(gòu)在經(jīng)歷瞬態(tài)增長后,可能會與第2模態(tài)之間發(fā)生非線性耦合作用[21]或因其條帶結(jié)構(gòu)本身的不穩(wěn)定性[22]而二次失穩(wěn)并最終導(dǎo)致轉(zhuǎn)捩的發(fā)生。如果采用足夠精細(xì)的數(shù)值模擬方法,通過對渦的形成和破碎過程進(jìn)行捕捉,理論上其可以捕獲流向渦的形成、發(fā)展、失穩(wěn)到轉(zhuǎn)捩為湍流的整個過程。要捕獲這些過程,目前可以采用的方法有DNS、LES和RANS/LES混合方法,而目前最有可能大范圍工程應(yīng)用的是RANS/LES混合方法。Yoon等[23]首先將DES方法應(yīng)用于高超聲速邊界層內(nèi)單個粗糙顆粒誘導(dǎo)轉(zhuǎn)捩的算例,結(jié)果表明DES模型可以分辨大渦結(jié)構(gòu),可以用來模擬粗糙顆粒誘導(dǎo)的強(qiáng)制轉(zhuǎn)捩;肖志祥等采用IDDES方法分別研究了高超聲速下粗糙顆粒誘導(dǎo)平板強(qiáng)制轉(zhuǎn)捩[24]和凹槽誘導(dǎo)轉(zhuǎn)捩[25]的情況,給出了轉(zhuǎn)捩渦結(jié)構(gòu)的發(fā)展過程,表明采用IDDES方法對強(qiáng)制轉(zhuǎn)捩開展研究是可行的,但其在進(jìn)行模擬時,采用的是強(qiáng)制粗糙顆?;蛘甙疾凵嫌蔚牧鲃訛槿珜恿鞯那闆r,因此在流過粗糙顆?;虬疾壑?,有一個從層流到湍流的轉(zhuǎn)捩過程。IDDES方法之所以能夠模擬強(qiáng)制轉(zhuǎn)捩,是因為其能夠解析更小尺度的渦,能夠捕獲粗糙顆?;虬疾壅T導(dǎo)出來的渦結(jié)構(gòu),但由于IDDES方法仍然是一個針對湍流的計算方法,當(dāng)沒有間歇因子來控制其湍動能方程時,針對沒有明顯大渦結(jié)構(gòu)、由小擾動經(jīng)過線性增長和非線性增長之后破碎再轉(zhuǎn)捩為湍流的過程,其卻無法捕獲。而且純粹的IDDES方法也無法反映來流湍流度對轉(zhuǎn)捩的影響。根據(jù)普渡大學(xué)對于X-51A的試驗結(jié)果[26],來流湍流度對粗糙顆粒誘導(dǎo)的轉(zhuǎn)捩仍然有重要的影響,因此來流湍流度也是一個不可忽略的因素。

    針對γ-Reθ轉(zhuǎn)捩模型和IDDES方法各自的優(yōu)點和缺點,如果將γ-Reθ轉(zhuǎn)捩模型與IDDES方法結(jié)合起來,則理論上可以克服各自的缺點,從而同時滿足對自然轉(zhuǎn)捩和強(qiáng)制轉(zhuǎn)捩的模擬。實際上,早在2011年Sorensen等[27]就嘗試過將γ-Reθ轉(zhuǎn)捩模型與DES方法結(jié)合起來,并進(jìn)行了大迎角機(jī)翼的計算,隨后Alam等[28]則將DDES方法與k-kL-ω轉(zhuǎn)捩模型[29]進(jìn)行了結(jié)合,并將其應(yīng)用在低壓渦輪機(jī)葉片的轉(zhuǎn)捩預(yù)測中。喬磊等[30]也嘗試過將DDES方法與γ-Reθ轉(zhuǎn)捩模型相結(jié)合,對低速平板算例進(jìn)行了計算。以上結(jié)果均表明轉(zhuǎn)捩模型與DES類方法結(jié)合后,能夠像原始的γ-Reθ轉(zhuǎn)捩模型一樣對自然轉(zhuǎn)捩進(jìn)行模擬。但以上均為試探性研究,研究的構(gòu)型都較為簡單,都是低速領(lǐng)域的轉(zhuǎn)捩,沒有進(jìn)一步發(fā)揮基于IDDES框架的轉(zhuǎn)捩模型相對于傳統(tǒng)轉(zhuǎn)捩模型能夠捕捉更精細(xì)的流場結(jié)構(gòu)從而更精確地模擬強(qiáng)制轉(zhuǎn)捩的這一優(yōu)勢。

    本文希望通過結(jié)合γ-Reθ轉(zhuǎn)捩模型對自然轉(zhuǎn)捩模擬的優(yōu)點和IDDES方法在模擬強(qiáng)制轉(zhuǎn)捩時能分辨更精細(xì)物理結(jié)構(gòu)從而捕獲轉(zhuǎn)捩過程的優(yōu)點,將γ-Reθ轉(zhuǎn)捩模型添加到IDDES方法的框架中去,從而實現(xiàn)對自然轉(zhuǎn)捩和強(qiáng)制轉(zhuǎn)捩更精確的模擬。在沒有粗糙顆粒、凹槽等強(qiáng)制轉(zhuǎn)捩因素時,網(wǎng)格只需滿足RANS的網(wǎng)格要求即可,此時在邊界層內(nèi)RANS方法會被激活,并對轉(zhuǎn)捩進(jìn)行模擬。在有粗糙顆粒、凹槽等強(qiáng)制轉(zhuǎn)捩因素時,采用更精細(xì)的網(wǎng)格,邊界層內(nèi)大部分區(qū)域LES方法被激活,可對渦結(jié)構(gòu)的出現(xiàn)及發(fā)展過程進(jìn)行捕捉,從而實現(xiàn)對強(qiáng)制轉(zhuǎn)捩的模擬。

    1 計算方法

    本文搭建的基于IDDES框架的的γ-Reθ轉(zhuǎn)捩模型是在課題組開發(fā)的大規(guī)模并行CFD平臺AHL3D[31]上添加模塊建立的,AHL3D 采用基于網(wǎng)格平均的有限體積方法離散求解Navier-Stokes方程組,定常時間推進(jìn)采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法,非定常計算則包括顯示的龍格-庫塔法和隱式的雙時間步法。無黏通量采用重構(gòu)-推進(jìn)方法,重構(gòu)采用3階MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)方法,本文在此基礎(chǔ)上又添加了5階WENO(Weighted Essentially Non-Oscillatory)重構(gòu)方法。無黏通量的推進(jìn)采用Steger-Warming通量分裂、AUSM+(improved Advection Upstream Splitting Method)、AUSMPW+(AUSM+ by Pressure-based Weighted function)和LDFSS(Low Diffusion Flux-Splitting Scheme)4種計算格式。黏性通量采用改進(jìn)的Gauss定理計算。湍流模型中包含了本文需要用到的剪切應(yīng)力輸運(SST)雙方程模型。作者團(tuán)隊在文獻(xiàn)[20]中已經(jīng)實現(xiàn)了γ-Reθ轉(zhuǎn)捩模型的搭建和應(yīng)用,因此本文以該轉(zhuǎn)捩模型為基礎(chǔ),先實現(xiàn)IDDES方法,然后將IDDES方法與γ-Reθ轉(zhuǎn)捩模型相結(jié)合,重新對轉(zhuǎn)捩模型進(jìn)行標(biāo)定,從而實現(xiàn)本文計算方法的建立。

    1.1 γ-Reθ轉(zhuǎn)捩模型

    γ-Reθ模型詳細(xì)的計算公式可參考文獻(xiàn)[9]。其通過間歇因子γ對SST模型湍動能k方程的產(chǎn)生項和破壞項進(jìn)行控制,進(jìn)而控制整個轉(zhuǎn)捩過程。間歇因子輸運方程為

    (1)

    式中:ρ為密度;Uj為xj方向的雷諾平均速度;μ為分子黏性系數(shù);μt為湍流黏性系數(shù);σf=1.0為常數(shù);Pγ和Eγ分別為γ方程的產(chǎn)生項和破壞項,Pγ的形式為

    Pγ=Flengthca1ρS(γFonset)0.5(1-ce1γ)

    (2)

    其中:Fonset為決定轉(zhuǎn)捩開始的閾函數(shù);Flength為轉(zhuǎn)捩區(qū)域長度;其余相關(guān)參數(shù)的含義詳見文獻(xiàn)[9]。

    (3)

    式中:Reν為當(dāng)?shù)販u雷諾數(shù);Reθc為臨界動量厚度雷諾數(shù)。

    (4)

    (5)

    (6)

    (7)

    Reθt=f3(Tu)f4(λθ)

    (8)

    此外,γ-Reθ模型還考慮了分離誘導(dǎo)轉(zhuǎn)捩的作用。

    (9)

    γeff=max(γ,γsep)

    (10)

    式中:γsep為分離區(qū)的間歇因子;γeff為最終的間歇因子;Freattach為與再附區(qū)有關(guān)的參數(shù),具體形式參見文獻(xiàn)[9]。其可通過改變經(jīng)驗參數(shù)s1來調(diào)節(jié)對分離誘導(dǎo)轉(zhuǎn)捩的敏感度。最后,轉(zhuǎn)捩模型對SST模型k方程的作用方式為

    (11)

    式中:Pk和Dk分別為SST模型中k方程的產(chǎn)生項和破壞項。

    1.2 γ-Reθ轉(zhuǎn)捩模型與IDDES方法的結(jié)合

    γ-Reθ模型通過添加兩個輸運方程求解間歇因子γ,且通過γ作用于k方程的產(chǎn)生項和破壞項從而與SST湍流模型結(jié)合。在形式上,添加基于SST模型的IDDES方法則更為簡單,可以通過修改湍動能輸運方程破壞項中的長度尺度函數(shù)而實現(xiàn)[32]。SST模型中湍動能方程為

    (12)

    (13)

    長度尺度函數(shù)從DES到IDDES發(fā)展歷程的具體情況可參見文獻(xiàn)[33],本文僅給出lIDDES的計算形式。lIDDES與網(wǎng)格尺度Δ有關(guān)。

    Δ=min{max[Cwdw,Cwhmax,hmin],hmax}

    (14)

    式中:hmax=max(Δx,Δy,Δz),hmin=min(Δx,Δy,Δz),Δx、Δy、Δz分別為x、y、z這3個方向上的網(wǎng)格尺度;dw為距壁面距離;Cw為與亞格子模型無關(guān)的經(jīng)驗常數(shù),一般取Cw=0.15。

    lIDDES的形式為

    (15)

    (16)

    fdt=1-tanh[(Cdtrdt)3]

    (17)

    (18)

    fB=min{2exp(-9α2),1.0}

    (19)

    其中:α=0.25-dw/hmax。至此,lIDDES求出。

    將RANS框架換成IDDES框架后,會對轉(zhuǎn)捩模型產(chǎn)生2個方面的影響:

    2 IDDES方法的驗證

    超聲速混合層流動是典型的以大尺度渦結(jié)構(gòu)為主導(dǎo)的可壓縮自由剪切層湍流流動。本節(jié)通過Goebel的混合層試驗[34]來檢驗本文的IDDES方法對大尺度渦結(jié)構(gòu)的捕捉能力以及對平均量和脈動量的預(yù)測能力。

    Goebel的試驗構(gòu)型[34]如圖1所示,混合區(qū)總高48 mm,總長500 mm,隔板下部傾角為2.5 °,尖劈厚度為0.5 mm。混合區(qū)入口參數(shù)如表1所示。計算中為滿足混合區(qū)入口處的附面層厚度,參考文獻(xiàn)[33]中的做法,將隔板上部入口選在混合區(qū)入口前67 mm處,隔板下部入口選在混合區(qū)入口前42 mm處。

    計算求解二維非定常Navier-Stokes方程組,采用基于SST模型的IDDES方法進(jìn)行模擬。無黏通量重構(gòu)采用五階WENO方法,分裂采用AUSMPW+格式,黏性通量采用改進(jìn)的高斯定理,時間推進(jìn)采用隱式雙時間步法,物理時間步長為5×10-8s,子迭代最大步數(shù)為15,子迭代CFL數(shù)為0.8。上下壁面采用對稱邊界條件,隔板采用絕熱壁條件,出口采用外推?;旌蠀^(qū)網(wǎng)格量為1 001×261,垂直壁面第1層網(wǎng)格法向長度取0.01 mm。隔板上下的初始條件分別設(shè)為上下入口的來流。計算經(jīng)過1 ms后流場充分發(fā)展并穩(wěn)定,再計算1 ms,統(tǒng)計1~2 ms時的平均量。

    計算得到的t=1 ms時的流場結(jié)構(gòu)如圖2所示。從圖中可以看出,本文IDDES方法能夠分辨混合層中大尺度渦結(jié)構(gòu)的形成、發(fā)展和脫落過程。t=1~2 ms期間的3個不同流向位置的時間統(tǒng)計平均流向速度U無量綱化之后的分布如圖3所示,圖中h=48 mm為混合區(qū)總高度,U2=399 m/s為混合層下層入口速度,ΔU=301 m/s為上下層入口速度之差。x=50 mm時,計算得到的混合層比試驗稍薄,這可能是由于來流條件中沒有給定脈動值,造成剪切層發(fā)展比試驗值要稍慢導(dǎo)致的。而到x=100 mm和150 mm時,計算與試驗吻合得非常好,此時旋渦開始變大并脫落,混合層厚度明顯增加。

    圖1 Goebel的混合層試驗示意圖[34]Fig.1 Sketch of Goebel’s mixing layer test[34]

    表1 混合區(qū)入口來流條件Table 1 Free stream conditions of mixing region inlet

    圖2 二維超聲速混合層流場結(jié)構(gòu)Fig.2 Flow field structures of two-dimensional supersonic mixing layer

    圖4和圖5則分別為流向脈動速度均方根urms和法向脈動速度均方根vrms的分布,流向脈動速度均方根在x=50 mm時仍然比試驗值偏小,這與平均速度顯示混合層偏薄相一致。而在x=100 mm 和150 mm位置的流向速度均方根則與試驗吻合良好,表明本文計算方法能夠有效捕捉到流場大渦結(jié)構(gòu)。而法向速度均方根則在x=50 mm時與試驗吻合較好,在x=100 mm和150 mm時較試驗值要明顯偏大,計算得到的流向、法向速度脈動在同一位置基本是同等量級的,但試驗測量到的流向脈動要大于法向脈動,計算與試驗的差別可能是由于計算是采用二維模擬而導(dǎo)致的??傮w來說,計算結(jié)果與試驗結(jié)果的對比表明本文計算方法足以分辨超聲速混合層中的大渦結(jié)構(gòu),可應(yīng)用于可壓縮湍流的模擬。

    圖3 不同流向位置的時間統(tǒng)計平均流向速度分布Fig.3 Distribution of time-average streamwise velocity in different streamwise positions

    圖4 流向脈動速度均方根分布Fig.4 Root-mean-square distribution of streamwise velocity fluctuations

    圖5 法向脈動速度均方根分布Fig.5 Root-mean-square distribution of vertical velocity fluctuations

    3 基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型的建立及驗證

    γ-Reθ轉(zhuǎn)捩模型可以反映雷諾數(shù)對轉(zhuǎn)捩的影響,也可以通過式(8)反映壓力梯度、湍流度對轉(zhuǎn)捩的影響,但由于其對轉(zhuǎn)捩位置的預(yù)測受到計算方法、格式等多方面因素的影響,因此每一個框架上搭建的γ-Reθ轉(zhuǎn)捩模型都可能需要進(jìn)行重新標(biāo)定。作者團(tuán)隊基于在文獻(xiàn)[19-20]中采用γ-Reθ模型進(jìn)行計算的經(jīng)驗,得到以下幾點結(jié)論:

    1)在AHL3D框架上搭建的基于RANS的γ-Reθ轉(zhuǎn)捩模型如果直接采用原始的經(jīng)驗關(guān)系式進(jìn)行高超聲速轉(zhuǎn)捩模擬,得到的轉(zhuǎn)捩位置整體會比實際轉(zhuǎn)捩位置大大提前。

    2)通過重新標(biāo)定將轉(zhuǎn)捩位置整體推遲后,添加Forsythe[35]的可壓縮性修正可有效反映可壓縮性對轉(zhuǎn)捩的影響,使得該方法能在較寬馬赫數(shù)范圍內(nèi)對轉(zhuǎn)捩進(jìn)行模擬。

    基于以上經(jīng)驗,本文擬對基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型進(jìn)行如下修改:

    1)加入Forsythe的可壓縮性修正。

    2)經(jīng)驗耦合關(guān)系式仍然采用Langtry的原始關(guān)系式,參考文獻(xiàn)[15]中的做法,將式(3)中的2.193調(diào)整為3.5。

    3)式(9)中,常數(shù)2會致使γ在分離區(qū)直接變?yōu)?,導(dǎo)致k值的增加,從而致使本應(yīng)由LES方法處理的區(qū)域被程序自動判斷為RANS處理區(qū)域,這樣會導(dǎo)致分離區(qū)的精細(xì)流動因采用RANS方法而被抹去,實際上,當(dāng)將RANS方法改為IDDES方法后,對分離流的捕捉要更加精細(xì),因此,為了適當(dāng)弱化模型對分離誘導(dǎo)轉(zhuǎn)捩?;墓δ埽?dāng)需要用IDDES方法來分辨流場內(nèi)的精細(xì)流動時,將式(9)中的常數(shù)2改為1。

    3.1 馬赫數(shù)為3.5的尖錐轉(zhuǎn)捩

    Chen和Malik[36]在NASA蘭利研究中心的馬赫數(shù)Ma=3.5低擾動風(fēng)洞中進(jìn)行了半錐角為5°的尖錐轉(zhuǎn)捩測量,其來流單位雷諾數(shù)Re由2.81 ×107/m變化到7.8×107/m,得到了相同馬赫數(shù)條件下雷諾數(shù)對轉(zhuǎn)捩位置的定量影響規(guī)律。采用該算例不僅可以驗證模型對可壓縮流轉(zhuǎn)捩的模擬能力,也可以驗證模型是否能正確反應(yīng)雷諾數(shù)對轉(zhuǎn)捩的影響規(guī)律。

    本文選取了其試驗中的6個不同來流雷諾數(shù)進(jìn)行計算,并與試驗得到的轉(zhuǎn)捩位置進(jìn)行對比。計算設(shè)置的來流條件如表2所示,表中T0和P0分別為來流總溫和總壓。尖錐迎角為0°,壁面溫度統(tǒng)一設(shè)為絕熱壁,來流湍流度設(shè)為0.4%,來流黏性比RT=1.0。計算采用軸對稱模型進(jìn)行,方法為基于IDDES框架的γ-Reθ轉(zhuǎn)捩模型,無黏通量采用WENO五階重構(gòu)和AUSMPW+格式。模型全長為381 mm。尖錐流向和法向的網(wǎng)格量為341×201,壁面第1層網(wǎng)格距離壁面0.001 mm。

    由于試驗是通過測量溫度恢復(fù)系數(shù)來判斷流動是層流還是湍流,而真實的溫度壁面條件在計算中并不好給定,因此本文采用摩阻系數(shù)來判斷轉(zhuǎn)捩位置,且只與試驗對比轉(zhuǎn)捩位置,計算結(jié)果如圖6所示,其中Cf為摩阻系數(shù),r為溫度恢復(fù)系數(shù)。從圖中可知,隨著來流雷諾數(shù)的增加,轉(zhuǎn)捩位置由約0.28 m逐漸前移至0.1 m左右。模型可以反映可壓縮條件下雷諾數(shù)變化對轉(zhuǎn)捩的影響。

    表2 馬赫數(shù)為3.5的尖錐轉(zhuǎn)捩來流條件Table 2 Free stream condition of sharp cone transition tests with Mach number of 3.5

    圖6 不同Re下的尖錐轉(zhuǎn)捩Fig.6 Transition of sharp cone with different Re

    3.2 馬赫數(shù)為8.1的雙鍥轉(zhuǎn)捩

    Neuenhahn和Olivier[37]在德國TH2高超聲速激波風(fēng)洞中進(jìn)行了來流馬赫數(shù)為8.1的雙鍥中邊界層與激波干擾的試驗。試驗研究了壁溫和熵層對激波邊界層干擾的影響。在研究者后續(xù)針對該問題的CFD研究中發(fā)現(xiàn),該試驗中包含了邊界層轉(zhuǎn)捩現(xiàn)象,因為單純的層流或湍流均無法正確得到壓力分布和分離區(qū)大小。其中Krause[16]和Zhang[17]等均采用該算例對可壓縮條件下的γ-Reθ轉(zhuǎn)捩模型進(jìn)行了驗證。本文采用二維計算,模擬了尖前緣、前緣半經(jīng)R為0.5 mm和1.0 mm這3個不同模型在同一來流條件下的試驗情況。構(gòu)型如圖7所示,第1楔長度L1=180 mm,傾斜角為9°,第2楔長度L2=255 mm,拐角為11.5°。來流馬赫數(shù)為8.1,總壓為5.7 MPa,總溫為1 430 K,單位雷諾數(shù)為3.8×106/m。壁溫設(shè)為300 K,來流湍流度設(shè)為0.5%,來流黏性比RT設(shè)為1。計算方法同3.1節(jié),網(wǎng)格總量約為3.4萬,第1、2楔面的流向網(wǎng)格數(shù)分別為121和101,法向網(wǎng)格數(shù)為131,壁面第1層網(wǎng)格距離壁面小于0.002 mm。

    計算得到的壁面壓力(p)分布如圖8所示,從圖中可知,層流與湍流計算出來的拐角位置的分離泡大小有很大區(qū)別,而試驗值均介于層流和湍流計算結(jié)果之間,確實存在明顯的邊界層轉(zhuǎn)捩現(xiàn)象。同時,分離泡的大小隨R由0 mm增加到1.0 mm而逐漸增大,表明轉(zhuǎn)捩位置隨著前緣鈍度的增大而逐漸后移。計算結(jié)果反應(yīng)出了邊界層轉(zhuǎn)捩給分離泡大小帶來的影響,也反應(yīng)出了前緣鈍度(0~1.0 mm范圍內(nèi))對轉(zhuǎn)捩的影響規(guī)律。

    圖7 雙楔試驗構(gòu)型示意圖[37]Fig.7 Sketch of double wedge geometry[37]

    圖8 不同前緣半徑下雙楔的壁面壓力分布Fig.8 Wall pressure distribution of double wedge with different leading edge radiuses

    3.3 平板上單個粗糙顆粒誘導(dǎo)的轉(zhuǎn)捩

    對基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型進(jìn)行修改并經(jīng)過自然轉(zhuǎn)捩算例驗證后,對單個粗糙顆粒誘導(dǎo)高超聲速邊界層轉(zhuǎn)捩的情況進(jìn)行模擬,分析轉(zhuǎn)捩過程。本文選取了Tirtey等在高超聲速VKI H3風(fēng)洞中進(jìn)行的強(qiáng)制轉(zhuǎn)捩試驗[38-39]進(jìn)行模擬。Tirtey等研究了圓柱形、鉆石型、斜坡型的單個粗糙顆粒誘導(dǎo)轉(zhuǎn)捩的情況,本文選取其圓柱形粗糙顆粒在高雷諾數(shù)下的情況進(jìn)行模擬。模型如圖9所示[38],一塊長 290 mm、寬100 mm、前緣鈍度為0.5 mm的平板,在距前緣60 mm處安裝了一個粗糙顆粒,其參數(shù)為:a=4 mm,b=1.02 mm。本文只模擬了其長170 mm、寬40 mm的區(qū)域,流場區(qū)域高20 mm,如圖10所示。來流馬赫數(shù)為6,來流單位雷諾數(shù)為2.7×107/m,來流總壓為3.1 MPa,總溫為500 K,壁溫在試驗開始時均為292±2 K,試驗結(jié)束后層流區(qū)最大溫升為4 K,湍流區(qū)最大溫升為12 K,由于這里壁溫變化都不大,因此選取300 K為壁面溫度。整個流場的網(wǎng)格量約為3 500萬,平板第1層網(wǎng)格距離壁面0.002 mm,對應(yīng)的y+小于0.1。粗糙顆粒下游的流場流向、法向、展向的網(wǎng)格分別為601×151×207,對粗糙顆粒周圍的流場進(jìn)行了加密。粗糙顆粒后緣的流向和展向網(wǎng)格間距分別為0.005 mm和0.03 mm,對應(yīng)無量綱網(wǎng)格尺度Δx+<0.1,Δz+<1,尾跡區(qū)的網(wǎng)格間距隨著流向距離的增大而逐漸增加,到x=170 mm位置時,Δx=0.6 mm,Δz=0.2 mm,對應(yīng)的Δx+≈18,Δz+≈6。

    計算得到的以Q等值面表示的流場渦結(jié)構(gòu)如圖11所示(圖中T為溫度),可以看出粗糙顆粒誘導(dǎo)出的馬蹄渦和剪切層除了向下游發(fā)展外,還有明顯的往展向拓寬趨勢,這符合湍流的性質(zhì),也與試驗結(jié)果相符。圖12給出了壁面Stanton數(shù)(St)分布。可知,采用基于RANS的γ-Reθ轉(zhuǎn)捩模型得到的粗糙顆粒尾跡在x=70 mm(粗糙顆粒后緣下游10 mm處)位置時寬度約為12 mm,在x=170 mm位置的寬度約為14 mm,試驗測量得到的尾跡寬度在x=70 mm時約為9 mm,在x=170 mm時約為22 mm,而采用基于IDDES方法的γ-Reθ轉(zhuǎn)捩模型得到的尾跡寬度在x=70 mm時約為9 mm,在x=170 mm時約為20 mm。顯然本文模擬結(jié)果與試驗情況更為接近。

    圖9 平板模型示意圖[38]Fig.9 Sketch of panel model geometry[38]

    圖10 流場及網(wǎng)格拓?fù)涫疽鈭DFig.10 Sketch of topology of flow field and grid

    圖11 以溫度染色的Q等值面圖Fig.11 Iso-surfaces of Q colored by temperature

    圖12 壁面St分布Fig.12 St distribution on wall

    圖13為本文計算得到的St在尾跡區(qū)(對稱面上)及非尾跡區(qū)的分布與試驗對比,可知,計算與試驗基本吻合,St在粗糙顆粒后緣均先劇烈增大(由粗糙顆粒后緣誘導(dǎo)出的尾跡渦導(dǎo)致)并迅速衰減,再經(jīng)歷一個相對緩慢的增長過程,達(dá)到極大值后開始隨著流向距離的增加逐漸減小。若定義St達(dá)到全湍流值之前的極大值點為轉(zhuǎn)捩結(jié)束位置,則計算得到的轉(zhuǎn)捩結(jié)束位置為x≈120 mm,對應(yīng)的轉(zhuǎn)捩雷諾數(shù)約為3.24×106,而試驗得到的轉(zhuǎn)捩結(jié)束位置為x≈108 mm,對應(yīng)的轉(zhuǎn)捩雷諾數(shù)約為2.92×106。二者在轉(zhuǎn)捩位置和尾跡展向發(fā)展的差別的原因可能是本文的計算方法雖然一方面捕捉到了粗糙顆粒尾渦的形成、發(fā)展及失穩(wěn)過程,另一方面也?;俗匀晦D(zhuǎn)捩時擾動發(fā)展導(dǎo)致轉(zhuǎn)捩的過程,由于來流擾動是以湍動能k值的形式給定的,而非直接給定具體形式的脈動,因此無法反映來流擾動與粗糙顆粒之間的耦合作用,而在該試驗中,風(fēng)洞為常規(guī)噪聲風(fēng)洞,因此邊界層內(nèi)粗糙顆粒之前的擾動可能已經(jīng)足夠大并與粗糙顆粒產(chǎn)生了耦合作用,而目前本文的方法并不能反映出這種耦合作用。

    圖13 尾跡區(qū)與非尾跡區(qū)的St分布Fig.13 St distribution of wake region and out of wake region.

    圖14為計算得到的粗糙顆粒附近對稱面流線和壁面摩擦力線圖與Tirtey給出的流場結(jié)構(gòu)對比。從圖中可以看出,粗糙顆粒前緣誘導(dǎo)出了微弱的斜激波,而在顆粒前緣上游則存在較大的分離區(qū),較大的分離渦又誘導(dǎo)出二次反向分離渦,致使粗糙顆粒周圍的摩擦力線存在3條分離線和3條再附著線,分離渦向下游發(fā)展,形成了馬蹄渦結(jié)構(gòu),數(shù)值模擬對粗糙顆粒附近的流場結(jié)構(gòu)的模擬與試驗基本一致。

    圖15給出了對稱面上的流場結(jié)構(gòu),分別為x方向的密度梯度、流向渦大小及速度脈動的均方根分布,其中u∞表示來流速度。從圖中可知,在對稱面上,粗糙顆粒上部除誘導(dǎo)出了斜激波外,還誘導(dǎo)出了較強(qiáng)的剪切層,該剪切層厚度很薄,擾動量級很大,但基本維持在邊界層外緣,實際上,在粗糙顆粒下游,距離壁面較近的底層也出現(xiàn)了剪切層,底層部分的剪切層在向下游發(fā)展的同時也在向上層發(fā)展,同時厚度逐漸增厚,渦結(jié)構(gòu)最終破碎,擾動擴(kuò)散至整個邊界層,擾動在x=120 mm位置附近達(dá)到一個極大值之后,擾動幅值不再繼續(xù)增加,此時轉(zhuǎn)捩已經(jīng)完成。

    圖14 粗糙顆粒附近流場結(jié)構(gòu)Fig.14 Flow field structure around roughness element

    圖15 對稱面流場結(jié)構(gòu)Fig.15 Flow field structures on symmetry plane

    圖16、圖17則分別為瞬態(tài)流場中不同法向位置的溫度分布和不同流向位置流向渦的渦量分布(圖中Δx′為不同流向截面與粗糙顆粒中心的距離)。從圖中可知,粗糙單元后存在兩層剪切層(與對稱面上的結(jié)果相一致),上部剪切層很薄,在向下游的傳播過程中逐漸擴(kuò)散到外流并最終因耗散作用而徹底消失。而下部分的剪切層則與粗糙單元兩側(cè)的馬蹄渦一開始分別各自向下游發(fā)展并緩慢擴(kuò)散,直到三簇渦在展向上越來越靠近并最終開始相互耦合,之后較大的渦結(jié)構(gòu)開始破碎成較小的渦結(jié)構(gòu),流動變得不穩(wěn)定,開始變成湍流。在x=120 mm處,即Δx′=60b時,渦結(jié)構(gòu)已經(jīng)完全破碎,渦量在展向分布趨于均勻,表明此時轉(zhuǎn)捩已經(jīng)完成,與前面壁面熱流系數(shù)和對稱面速度脈動均方根分布體現(xiàn)的結(jié)果一致。

    圖16 壁面不同法向位置的瞬態(tài)溫度分布Fig.16 Instantaneous temperature distribution of different vertical distances from wall

    圖17 不同流向位置的流向渦渦量云圖Fig.17 Streamwise vorticity magnitude contours of different streamwise positions

    4 結(jié) 論

    1)將IDDES方法與γ-Reθ轉(zhuǎn)捩模型結(jié)合起來并發(fā)揮各自的優(yōu)勢,從而建立形成了新的基于IDDES框架的轉(zhuǎn)捩模型。

    2)通過超聲速混合層、高超聲速自然轉(zhuǎn)捩及粗糙顆粒誘導(dǎo)的強(qiáng)制轉(zhuǎn)捩等一系列算例的計算,表明新建立的方法在模擬自然轉(zhuǎn)捩時,能夠像γ-Reθ模型一樣通過基于轉(zhuǎn)捩動量厚度的經(jīng)驗關(guān)系式得到轉(zhuǎn)捩位置,而在模擬粗糙顆粒誘導(dǎo)的強(qiáng)制轉(zhuǎn)捩時,又能發(fā)揮IDDES方法能夠捕捉到粗糙顆粒在邊界層中誘導(dǎo)出的擾動的特點,能夠?qū)Υ植陬w粒誘導(dǎo)出的剪切層和馬蹄渦結(jié)構(gòu)的形成、耦合和失穩(wěn)過程進(jìn)行較為精細(xì)地刻畫和模擬,從而達(dá)到模擬強(qiáng)制轉(zhuǎn)捩的目的。

    3)計算結(jié)果表明新建立的方法在強(qiáng)制轉(zhuǎn)捩模擬上要優(yōu)于傳統(tǒng)的γ-Reθ轉(zhuǎn)捩模型。

    猜你喜歡
    來流邊界層湍流
    兩種典型來流條件下風(fēng)力機(jī)尾跡特性的數(shù)值研究
    能源工程(2022年2期)2022-05-23 13:51:48
    基于HIFiRE-2超燃發(fā)動機(jī)內(nèi)流道的激波邊界層干擾分析
    不同來流條件對溢洪道過流能力的影響
    重氣瞬時泄漏擴(kuò)散的湍流模型驗證
    一類具有邊界層性質(zhì)的二次奇攝動邊值問題
    彈發(fā)匹配驗證試驗系統(tǒng)來流快速啟動技術(shù)研究
    非特征邊界的MHD方程的邊界層
    “青春期”湍流中的智慧引渡(三)
    “青春期”湍流中的智慧引渡(二)
    弱分層湍流輸運特性的統(tǒng)計分析
    19禁男女啪啪无遮挡网站| 9191精品国产免费久久| 热re99久久精品国产66热6| 久久精品亚洲av国产电影网| 不卡一级毛片| 日本撒尿小便嘘嘘汇集6| 免费搜索国产男女视频| 乱人伦中国视频| 亚洲色图av天堂| 久久国产亚洲av麻豆专区| 久久性视频一级片| 天堂中文最新版在线下载| 深夜精品福利| 欧美 亚洲 国产 日韩一| 国产成人啪精品午夜网站| 91国产中文字幕| 亚洲性夜色夜夜综合| 亚洲精品久久成人aⅴ小说| 欧美日韩黄片免| 亚洲av成人一区二区三| 日韩欧美在线二视频| 这个男人来自地球电影免费观看| 狠狠狠狠99中文字幕| 999久久久国产精品视频| 在线观看免费高清a一片| 亚洲男人的天堂狠狠| 黄色片一级片一级黄色片| 久久精品91无色码中文字幕| 亚洲国产精品一区二区三区在线| 国产在线观看jvid| 少妇粗大呻吟视频| 欧美乱色亚洲激情| 日日摸夜夜添夜夜添小说| 级片在线观看| 狂野欧美激情性xxxx| 波多野结衣一区麻豆| 桃红色精品国产亚洲av| 久久人妻福利社区极品人妻图片| 天天躁夜夜躁狠狠躁躁| 欧美日韩精品网址| 真人做人爱边吃奶动态| 久久影院123| 多毛熟女@视频| 欧美成人性av电影在线观看| 欧美黄色片欧美黄色片| 国产真人三级小视频在线观看| 动漫黄色视频在线观看| 曰老女人黄片| 手机成人av网站| av电影中文网址| 日韩大尺度精品在线看网址 | 男女下面进入的视频免费午夜 | 亚洲成av片中文字幕在线观看| 麻豆久久精品国产亚洲av | 啦啦啦免费观看视频1| 人人妻,人人澡人人爽秒播| 亚洲中文字幕日韩| 多毛熟女@视频| 女人精品久久久久毛片| 午夜精品在线福利| 久久99一区二区三区| 久久久精品欧美日韩精品| 男女床上黄色一级片免费看| 欧美在线黄色| 国产一区二区三区视频了| 怎么达到女性高潮| www.精华液| 很黄的视频免费| 一区二区三区精品91| 久久精品国产亚洲av高清一级| 丝袜美腿诱惑在线| 不卡一级毛片| 亚洲欧美精品综合一区二区三区| 亚洲av片天天在线观看| 最新美女视频免费是黄的| 欧美一级毛片孕妇| 国产成人一区二区三区免费视频网站| 一个人观看的视频www高清免费观看 | 手机成人av网站| 国产亚洲精品综合一区在线观看 | 夫妻午夜视频| 免费不卡黄色视频| 日本vs欧美在线观看视频| 亚洲人成电影观看| 桃红色精品国产亚洲av| 欧美精品一区二区免费开放| 如日韩欧美国产精品一区二区三区| 国产免费现黄频在线看| 国产精品久久久人人做人人爽| www.自偷自拍.com| 99久久人妻综合| 涩涩av久久男人的天堂| 妹子高潮喷水视频| 一边摸一边做爽爽视频免费| 国产乱人伦免费视频| 亚洲欧美日韩另类电影网站| 亚洲一区中文字幕在线| 国产99久久九九免费精品| 成人18禁在线播放| 成人18禁在线播放| 精品一区二区三区av网在线观看| 午夜免费鲁丝| www.精华液| 极品人妻少妇av视频| 精品久久久久久电影网| 视频区图区小说| 久久性视频一级片| 精品久久久久久电影网| 97超级碰碰碰精品色视频在线观看| 精品久久蜜臀av无| √禁漫天堂资源中文www| 亚洲国产中文字幕在线视频| 一级片'在线观看视频| aaaaa片日本免费| 99香蕉大伊视频| 在线观看日韩欧美| 在线观看日韩欧美| e午夜精品久久久久久久| 国产国语露脸激情在线看| 一级片免费观看大全| 日本a在线网址| 少妇裸体淫交视频免费看高清 | 一区二区三区国产精品乱码| 黄色a级毛片大全视频| 国产精品一区二区免费欧美| 丰满迷人的少妇在线观看| 欧洲精品卡2卡3卡4卡5卡区| www国产在线视频色| 十分钟在线观看高清视频www| 亚洲一区高清亚洲精品| 欧美日韩一级在线毛片| 一边摸一边抽搐一进一小说| 亚洲黑人精品在线| 日本a在线网址| 精品国产乱码久久久久久男人| 亚洲av电影在线进入| 国产亚洲欧美在线一区二区| 中文字幕人妻丝袜一区二区| 国产精品免费一区二区三区在线| 丰满的人妻完整版| 两个人免费观看高清视频| 巨乳人妻的诱惑在线观看| 亚洲美女黄片视频| 国产av在哪里看| 国产av在哪里看| 国产男靠女视频免费网站| av天堂在线播放| 波多野结衣高清无吗| 国产深夜福利视频在线观看| 自线自在国产av| 韩国精品一区二区三区| 18禁观看日本| 自线自在国产av| 婷婷精品国产亚洲av在线| 成年版毛片免费区| 精品少妇一区二区三区视频日本电影| 国产视频一区二区在线看| 欧美国产精品va在线观看不卡| 亚洲av成人一区二区三| 成人精品一区二区免费| 国产亚洲欧美在线一区二区| 午夜福利在线免费观看网站| 免费女性裸体啪啪无遮挡网站| 欧美日韩乱码在线| a级毛片在线看网站| 日本免费a在线| 亚洲av日韩精品久久久久久密| 国产精品久久久av美女十八| 女人被躁到高潮嗷嗷叫费观| 亚洲欧美激情在线| 成人特级黄色片久久久久久久| 日本黄色视频三级网站网址| 午夜精品国产一区二区电影| 亚洲av成人不卡在线观看播放网| 亚洲伊人色综图| 无遮挡黄片免费观看| 国产又色又爽无遮挡免费看| 国产精品国产av在线观看| 18禁国产床啪视频网站| 免费在线观看完整版高清| 成人免费观看视频高清| 大型av网站在线播放| 51午夜福利影视在线观看| 欧美成人性av电影在线观看| 超色免费av| 男人操女人黄网站| 国产国语露脸激情在线看| 亚洲欧美激情综合另类| 久久久国产成人精品二区 | 不卡av一区二区三区| 成年版毛片免费区| 黄色成人免费大全| www.自偷自拍.com| 午夜影院日韩av| 亚洲第一av免费看| 黄色女人牲交| 亚洲激情在线av| 在线免费观看的www视频| 99国产精品免费福利视频| 欧美激情极品国产一区二区三区| 国产一区二区在线av高清观看| 亚洲狠狠婷婷综合久久图片| 老鸭窝网址在线观看| 午夜福利在线观看吧| 99riav亚洲国产免费| 麻豆国产av国片精品| 亚洲人成电影免费在线| 国产成人一区二区三区免费视频网站| 夜夜夜夜夜久久久久| 久久国产精品人妻蜜桃| 久久青草综合色| 侵犯人妻中文字幕一二三四区| 日韩免费av在线播放| 精品国产一区二区三区四区第35| 99在线人妻在线中文字幕| 99精品欧美一区二区三区四区| 国产三级黄色录像| 亚洲精品中文字幕一二三四区| 制服诱惑二区| 成熟少妇高潮喷水视频| 亚洲精品av麻豆狂野| 热re99久久国产66热| 亚洲国产精品合色在线| 欧美+亚洲+日韩+国产| av视频免费观看在线观看| 免费高清视频大片| 深夜精品福利| 免费高清在线观看日韩| 日韩av在线大香蕉| 两性夫妻黄色片| 自拍欧美九色日韩亚洲蝌蚪91| 国产精品一区二区免费欧美| 怎么达到女性高潮| 亚洲一区二区三区不卡视频| 亚洲专区国产一区二区| 欧美亚洲日本最大视频资源| 在线视频色国产色| 亚洲精品一二三| 免费观看人在逋| 欧美乱妇无乱码| 性少妇av在线| xxx96com| 在线十欧美十亚洲十日本专区| 999久久久精品免费观看国产| 成人av一区二区三区在线看| 亚洲精品久久成人aⅴ小说| 精品久久久久久久毛片微露脸| 黄片大片在线免费观看| 桃色一区二区三区在线观看| 国产成人欧美在线观看| 最近最新中文字幕大全电影3 | 女人爽到高潮嗷嗷叫在线视频| 国产精品亚洲一级av第二区| 久久久久亚洲av毛片大全| 亚洲精品一二三| 久久中文字幕人妻熟女| 久久九九热精品免费| 国产单亲对白刺激| 午夜免费鲁丝| 麻豆成人av在线观看| 人人妻,人人澡人人爽秒播| 亚洲,欧美精品.| 精品国产乱码久久久久久男人| 久久精品成人免费网站| 国产麻豆69| 在线播放国产精品三级| 国产精品久久久久久人妻精品电影| 午夜福利在线观看吧| 精品无人区乱码1区二区| 欧美久久黑人一区二区| 欧美不卡视频在线免费观看 | 一边摸一边做爽爽视频免费| 女人高潮潮喷娇喘18禁视频| 亚洲精品美女久久久久99蜜臀| 色精品久久人妻99蜜桃| 亚洲成国产人片在线观看| 看片在线看免费视频| 亚洲人成网站在线播放欧美日韩| 国产深夜福利视频在线观看| 亚洲,欧美精品.| 国产精品久久视频播放| 狂野欧美激情性xxxx| xxxhd国产人妻xxx| 99国产精品一区二区三区| 免费在线观看影片大全网站| 免费av毛片视频| 日韩av在线大香蕉| 国产精品亚洲av一区麻豆| 看黄色毛片网站| 亚洲 欧美 日韩 在线 免费| 欧美在线黄色| 日韩视频一区二区在线观看| 一边摸一边抽搐一进一小说| 日本vs欧美在线观看视频| 国产片内射在线| 一级a爱视频在线免费观看| 欧美在线黄色| 天堂中文最新版在线下载| 国产欧美日韩综合在线一区二区| 成人影院久久| 国产精品亚洲一级av第二区| 午夜影院日韩av| 久久久精品欧美日韩精品| 国产精品美女特级片免费视频播放器 | 色综合婷婷激情| 久热爱精品视频在线9| 国产在线观看jvid| 欧美日韩乱码在线| 国产区一区二久久| xxxhd国产人妻xxx| 热re99久久精品国产66热6| 亚洲熟妇熟女久久| 久久久国产成人免费| 91大片在线观看| 国产乱人伦免费视频| 欧美av亚洲av综合av国产av| 淫秽高清视频在线观看| 久久久国产成人免费| 人人澡人人妻人| 侵犯人妻中文字幕一二三四区| 巨乳人妻的诱惑在线观看| 黄色视频,在线免费观看| 69精品国产乱码久久久| 国产精品一区二区在线不卡| 亚洲熟女毛片儿| 久久亚洲精品不卡| 99精国产麻豆久久婷婷| 男人操女人黄网站| 国产人伦9x9x在线观看| 在线十欧美十亚洲十日本专区| 乱人伦中国视频| 久久久水蜜桃国产精品网| 男人的好看免费观看在线视频 | 美女国产高潮福利片在线看| 三上悠亚av全集在线观看| 国产蜜桃级精品一区二区三区| 18禁美女被吸乳视频| 国产极品粉嫩免费观看在线| 最近最新中文字幕大全电影3 | 亚洲伊人色综图| 19禁男女啪啪无遮挡网站| 成人手机av| 国产精品一区二区在线不卡| 欧美激情久久久久久爽电影 | 久久久久久人人人人人| 丰满的人妻完整版| 亚洲精华国产精华精| 欧美一级毛片孕妇| 亚洲专区字幕在线| 无遮挡黄片免费观看| 日本免费a在线| 精品欧美一区二区三区在线| 不卡av一区二区三区| 天天躁夜夜躁狠狠躁躁| 99精品欧美一区二区三区四区| 18禁裸乳无遮挡免费网站照片 | 欧美日韩精品网址| 国产欧美日韩一区二区三| 高清黄色对白视频在线免费看| 中文字幕色久视频| 欧美不卡视频在线免费观看 | 亚洲精品一区av在线观看| 久久中文看片网| 1024视频免费在线观看| www.熟女人妻精品国产| 亚洲五月色婷婷综合| 亚洲,欧美精品.| 一区在线观看完整版| 久久国产亚洲av麻豆专区| 国产高清视频在线播放一区| 日韩三级视频一区二区三区| 日本一区二区免费在线视频| 国产蜜桃级精品一区二区三区| 日韩av在线大香蕉| 高清在线国产一区| 丰满饥渴人妻一区二区三| 亚洲熟妇熟女久久| 免费在线观看亚洲国产| 欧美国产精品va在线观看不卡| avwww免费| 国产精品 国内视频| 久久久久精品国产欧美久久久| 天堂俺去俺来也www色官网| 亚洲久久久国产精品| 一级作爱视频免费观看| 99国产精品99久久久久| 精品人妻1区二区| 男女之事视频高清在线观看| 久久精品国产99精品国产亚洲性色 | 国产精品影院久久| 欧美日韩乱码在线| 男女高潮啪啪啪动态图| 精品国产亚洲在线| 亚洲av电影在线进入| 国产成人欧美在线观看| 日韩精品中文字幕看吧| 欧美黑人精品巨大| 桃红色精品国产亚洲av| 日韩大尺度精品在线看网址 | 午夜两性在线视频| 国产欧美日韩精品亚洲av| 午夜精品久久久久久毛片777| 热re99久久国产66热| 最好的美女福利视频网| 国产成人精品久久二区二区91| www.熟女人妻精品国产| 国产一区二区三区视频了| 亚洲av成人av| 香蕉国产在线看| 看免费av毛片| 精品国产超薄肉色丝袜足j| 他把我摸到了高潮在线观看| 黄频高清免费视频| 亚洲欧美精品综合久久99| 高清欧美精品videossex| 亚洲美女黄片视频| 欧美日韩乱码在线| 人妻久久中文字幕网| 久久国产精品影院| 欧美老熟妇乱子伦牲交| 国产国语露脸激情在线看| 啦啦啦在线免费观看视频4| 色哟哟哟哟哟哟| 老汉色av国产亚洲站长工具| tocl精华| 成人精品一区二区免费| 色在线成人网| 一边摸一边抽搐一进一小说| 亚洲精品国产色婷婷电影| 成人av一区二区三区在线看| 淫秽高清视频在线观看| 国产麻豆69| 欧美乱色亚洲激情| 久久狼人影院| 麻豆av在线久日| 国产日韩一区二区三区精品不卡| 亚洲精品中文字幕在线视频| 日韩三级视频一区二区三区| 一级作爱视频免费观看| 亚洲精品久久午夜乱码| 精品国产超薄肉色丝袜足j| 日韩 欧美 亚洲 中文字幕| 国产精品久久久久久人妻精品电影| 国产高清国产精品国产三级| 亚洲全国av大片| 91av网站免费观看| 欧美另类亚洲清纯唯美| 欧美国产精品va在线观看不卡| 黄色毛片三级朝国网站| 国产视频一区二区在线看| 国产精品一区二区三区四区久久 | 亚洲狠狠婷婷综合久久图片| 精品欧美一区二区三区在线| 国产精品亚洲一级av第二区| 视频区欧美日本亚洲| 黄色片一级片一级黄色片| 国产精品国产高清国产av| 国产精品二区激情视频| 夜夜夜夜夜久久久久| 女同久久另类99精品国产91| 久久久久久大精品| www.熟女人妻精品国产| 韩国精品一区二区三区| 老鸭窝网址在线观看| 国产成人免费无遮挡视频| 日韩 欧美 亚洲 中文字幕| 亚洲色图 男人天堂 中文字幕| 亚洲va日本ⅴa欧美va伊人久久| 人人妻人人爽人人添夜夜欢视频| 首页视频小说图片口味搜索| 亚洲精品一卡2卡三卡4卡5卡| 国产一区二区在线av高清观看| 欧美中文日本在线观看视频| 亚洲专区国产一区二区| 国产一区二区激情短视频| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲精品在线美女| 国内毛片毛片毛片毛片毛片| 美女扒开内裤让男人捅视频| 黄色 视频免费看| 51午夜福利影视在线观看| 久久久久久久久免费视频了| 国产野战对白在线观看| 日韩大码丰满熟妇| 咕卡用的链子| 黄片小视频在线播放| 两人在一起打扑克的视频| 亚洲五月色婷婷综合| 精品久久久久久久毛片微露脸| 高清在线国产一区| 亚洲伊人色综图| 国产精品久久视频播放| 99久久99久久久精品蜜桃| 亚洲精品一二三| 十分钟在线观看高清视频www| 亚洲精品国产区一区二| 亚洲精品成人av观看孕妇| 午夜福利在线免费观看网站| 久久久久久久久免费视频了| 国产高清激情床上av| 国产成人一区二区三区免费视频网站| 香蕉久久夜色| 男女午夜视频在线观看| 黄色视频不卡| 男人操女人黄网站| 国产精品99久久99久久久不卡| 国产av精品麻豆| xxxhd国产人妻xxx| 国产91精品成人一区二区三区| 中文字幕色久视频| 好男人电影高清在线观看| 丝袜人妻中文字幕| 国产精品乱码一区二三区的特点 | 久久亚洲精品不卡| 亚洲成a人片在线一区二区| 亚洲精品在线美女| 18禁美女被吸乳视频| 日日摸夜夜添夜夜添小说| 免费av毛片视频| 婷婷丁香在线五月| 国产av在哪里看| 99香蕉大伊视频| 亚洲人成网站在线播放欧美日韩| 久久99一区二区三区| 精品国产国语对白av| 免费高清视频大片| 黑人猛操日本美女一级片| 亚洲国产精品999在线| 女警被强在线播放| 亚洲 欧美 日韩 在线 免费| 1024视频免费在线观看| 国产精品乱码一区二三区的特点 | 精品乱码久久久久久99久播| 亚洲成国产人片在线观看| 日韩免费av在线播放| 成人手机av| 麻豆成人av在线观看| 亚洲一码二码三码区别大吗| 精品久久久精品久久久| 亚洲中文av在线| 高潮久久久久久久久久久不卡| 亚洲激情在线av| 国产精品 欧美亚洲| 男人舔女人下体高潮全视频| 久久精品国产99精品国产亚洲性色 | 一边摸一边抽搐一进一小说| 两性午夜刺激爽爽歪歪视频在线观看 | 757午夜福利合集在线观看| 久久中文字幕一级| 窝窝影院91人妻| 宅男免费午夜| 国内久久婷婷六月综合欲色啪| 精品国产一区二区三区四区第35| 夜夜夜夜夜久久久久| 97碰自拍视频| 日韩一卡2卡3卡4卡2021年| 99国产极品粉嫩在线观看| 夜夜夜夜夜久久久久| 91成人精品电影| 久久这里只有精品19| 三上悠亚av全集在线观看| 久久午夜亚洲精品久久| 中文亚洲av片在线观看爽| 久久久久久久精品吃奶| 色婷婷久久久亚洲欧美| 精品国产乱码久久久久久男人| 91在线观看av| 90打野战视频偷拍视频| 免费一级毛片在线播放高清视频 | 国产在线观看jvid| 日本撒尿小便嘘嘘汇集6| 亚洲 国产 在线| 99香蕉大伊视频| 亚洲一区二区三区色噜噜 | 亚洲欧美激情综合另类| 免费在线观看影片大全网站| www.www免费av| 亚洲精品一卡2卡三卡4卡5卡| 99国产精品一区二区蜜桃av| 两性午夜刺激爽爽歪歪视频在线观看 | 不卡av一区二区三区| 国内毛片毛片毛片毛片毛片| netflix在线观看网站| 久久久久久久午夜电影 | 在线永久观看黄色视频| 多毛熟女@视频| 国产91精品成人一区二区三区| 成人18禁在线播放| av欧美777| 校园春色视频在线观看| 久久狼人影院| 熟女少妇亚洲综合色aaa.| 日韩精品免费视频一区二区三区| 老司机福利观看| av欧美777| 成人18禁在线播放| 国产欧美日韩一区二区三| 级片在线观看| 99热只有精品国产| 国产精品一区二区免费欧美| 18禁国产床啪视频网站| 最新美女视频免费是黄的| 日韩欧美三级三区| 精品久久久久久久久久免费视频 | 欧美日韩乱码在线| 亚洲专区字幕在线| 男女之事视频高清在线观看| 老司机深夜福利视频在线观看| 大陆偷拍与自拍| 成人永久免费在线观看视频| 亚洲九九香蕉| 欧美乱色亚洲激情| 亚洲国产精品一区二区三区在线| 老熟妇乱子伦视频在线观看| 亚洲国产欧美一区二区综合|