• <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)計分析
    亚洲精品一二三| 最近中文字幕高清免费大全6| 久久精品久久久久久噜噜老黄| 欧美日韩综合久久久久久| 久久国产精品男人的天堂亚洲 | 亚洲第一av免费看| 永久免费av网站大全| 黄色配什么色好看| 国产av一区二区精品久久| 天堂俺去俺来也www色官网| 国产成人免费无遮挡视频| 亚洲精品乱久久久久久| a级毛片在线看网站| 99国产精品免费福利视频| 国产亚洲av片在线观看秒播厂| 一二三四在线观看免费中文在 | 久久午夜综合久久蜜桃| 精品一区在线观看国产| 久久国产亚洲av麻豆专区| videossex国产| 青青草视频在线视频观看| 亚洲精品国产色婷婷电影| 欧美精品av麻豆av| 日本av手机在线免费观看| 亚洲精品国产色婷婷电影| 国产高清国产精品国产三级| 有码 亚洲区| 欧美日韩精品成人综合77777| 亚洲精品久久久久久婷婷小说| 交换朋友夫妻互换小说| 肉色欧美久久久久久久蜜桃| 男女无遮挡免费网站观看| 亚洲精品一区蜜桃| 黄色 视频免费看| 午夜福利影视在线免费观看| 高清不卡的av网站| 精品少妇内射三级| 日本欧美国产在线视频| 国产精品久久久久久精品古装| 久久精品国产a三级三级三级| 国产高清三级在线| 国产精品一区www在线观看| 夜夜骑夜夜射夜夜干| 久久精品久久精品一区二区三区| 在线观看国产h片| 国产成人91sexporn| 极品人妻少妇av视频| 欧美精品人与动牲交sv欧美| 交换朋友夫妻互换小说| 欧美精品人与动牲交sv欧美| videosex国产| 久久久久久久久久久免费av| 国产成人精品在线电影| 国产一区二区三区综合在线观看 | 肉色欧美久久久久久久蜜桃| 久久久国产精品麻豆| 在线天堂中文资源库| 国产熟女午夜一区二区三区| 久久99热6这里只有精品| 久久97久久精品| 只有这里有精品99| 伦精品一区二区三区| 精品亚洲乱码少妇综合久久| 日韩成人av中文字幕在线观看| 男的添女的下面高潮视频| 制服丝袜香蕉在线| 国产亚洲精品第一综合不卡 | 国产毛片在线视频| 麻豆精品久久久久久蜜桃| 18禁在线无遮挡免费观看视频| 少妇的逼好多水| 免费观看在线日韩| 精品人妻在线不人妻| 高清在线视频一区二区三区| 卡戴珊不雅视频在线播放| av片东京热男人的天堂| 成人手机av| 亚洲国产看品久久| 国内精品宾馆在线| 亚洲欧洲日产国产| 欧美精品一区二区免费开放| 69精品国产乱码久久久| 一级,二级,三级黄色视频| 男女下面插进去视频免费观看 | 久久青草综合色| 在线观看国产h片| 91久久精品国产一区二区三区| 精品午夜福利在线看| 十分钟在线观看高清视频www| 亚洲熟女精品中文字幕| 国产精品国产三级专区第一集| 亚洲四区av| 欧美精品av麻豆av| 黑人欧美特级aaaaaa片| 国产男女内射视频| 亚洲综合色网址| 国产精品国产三级国产av玫瑰| 少妇的逼好多水| 国产激情久久老熟女| 伦理电影免费视频| 一级片'在线观看视频| 日韩大片免费观看网站| 七月丁香在线播放| 日韩中字成人| 2022亚洲国产成人精品| 交换朋友夫妻互换小说| 满18在线观看网站| 国产黄频视频在线观看| 日本午夜av视频| 婷婷色综合大香蕉| 中国三级夫妇交换| 日本与韩国留学比较| 精品人妻偷拍中文字幕| 欧美+日韩+精品| 日韩精品免费视频一区二区三区 | kizo精华| 亚洲,欧美精品.| 中国美白少妇内射xxxbb| 日韩欧美一区视频在线观看| 久久久精品区二区三区| 亚洲成人手机| 在线观看免费高清a一片| 亚洲精品久久午夜乱码| 久久ye,这里只有精品| 一本色道久久久久久精品综合| 亚洲av欧美aⅴ国产| 国产成人精品无人区| 九九在线视频观看精品| 满18在线观看网站| 母亲3免费完整高清在线观看 | 午夜精品国产一区二区电影| 水蜜桃什么品种好| 成人无遮挡网站| 久久综合国产亚洲精品| 精品午夜福利在线看| 18禁裸乳无遮挡动漫免费视频| 99香蕉大伊视频| 天天躁夜夜躁狠狠久久av| 最新的欧美精品一区二区| 亚洲精品自拍成人| 只有这里有精品99| 亚洲精品乱码久久久久久按摩| 黑丝袜美女国产一区| 久久久欧美国产精品| 精品亚洲成国产av| 成人综合一区亚洲| 香蕉精品网在线| 亚洲欧美日韩另类电影网站| 午夜老司机福利剧场| 亚洲国产色片| 久久国内精品自在自线图片| 天堂8中文在线网| 狂野欧美激情性xxxx在线观看| 人人妻人人澡人人看| 中国三级夫妇交换| 免费播放大片免费观看视频在线观看| 日韩av在线免费看完整版不卡| 国产女主播在线喷水免费视频网站| 伊人亚洲综合成人网| 一区二区三区精品91| 国产精品成人在线| 少妇人妻久久综合中文| 国产欧美日韩综合在线一区二区| 免费观看a级毛片全部| 黄片无遮挡物在线观看| freevideosex欧美| 免费观看av网站的网址| 在线精品无人区一区二区三| 免费大片黄手机在线观看| 欧美激情国产日韩精品一区| 午夜免费观看性视频| 亚洲精华国产精华液的使用体验| 不卡视频在线观看欧美| 深夜精品福利| 日韩成人伦理影院| 国产探花极品一区二区| 99国产精品免费福利视频| 久久国产精品大桥未久av| 岛国毛片在线播放| 亚洲国产精品999| 欧美激情国产日韩精品一区| 亚洲精品一二三| 少妇的逼水好多| 18+在线观看网站| 岛国毛片在线播放| 大陆偷拍与自拍| 亚洲美女搞黄在线观看| 在现免费观看毛片| 亚洲欧美成人精品一区二区| 国产xxxxx性猛交| 欧美3d第一页| 热re99久久精品国产66热6| 亚洲欧美日韩另类电影网站| 亚洲五月色婷婷综合| 午夜影院在线不卡| 色94色欧美一区二区| 这个男人来自地球电影免费观看 | 欧美性感艳星| 人人澡人人妻人| 国产成人免费无遮挡视频| av电影中文网址| 国产精品嫩草影院av在线观看| 精品少妇内射三级| 日日爽夜夜爽网站| 多毛熟女@视频| 色吧在线观看| 国产成人91sexporn| 蜜臀久久99精品久久宅男| 九九在线视频观看精品| www.av在线官网国产| kizo精华| 伦精品一区二区三区| 国产成人精品无人区| 亚洲伊人久久精品综合| 国产极品粉嫩免费观看在线| 成人国产麻豆网| 国产一区二区在线观看av| 一本色道久久久久久精品综合| 91久久精品国产一区二区三区| 国产成人免费无遮挡视频| 美女国产高潮福利片在线看| 国产欧美日韩一区二区三区在线| 91成人精品电影| 国产一区二区激情短视频 | 国产成人免费无遮挡视频| 2021少妇久久久久久久久久久| 五月天丁香电影| a级毛色黄片| 免费人成在线观看视频色| 日韩在线高清观看一区二区三区| 久久精品国产亚洲av涩爱| videosex国产| 最近最新中文字幕大全免费视频 | 男女边摸边吃奶| 免费观看a级毛片全部| 大码成人一级视频| 日本爱情动作片www.在线观看| 欧美日韩av久久| 少妇精品久久久久久久| 欧美97在线视频| 肉色欧美久久久久久久蜜桃| 亚洲欧美一区二区三区黑人 | 久久久久久久久久人人人人人人| 日韩av不卡免费在线播放| av片东京热男人的天堂| 亚洲,一卡二卡三卡| 亚洲五月色婷婷综合| 久久精品国产亚洲av天美| 国产精品国产av在线观看| 亚洲性久久影院| 美女脱内裤让男人舔精品视频| 精品人妻偷拍中文字幕| 最新中文字幕久久久久| 亚洲欧美成人综合另类久久久| 女人精品久久久久毛片| 国产激情久久老熟女| 久久韩国三级中文字幕| 欧美 亚洲 国产 日韩一| 中文字幕免费在线视频6| 国产日韩欧美在线精品| 中国美白少妇内射xxxbb| 精品国产一区二区三区久久久樱花| 最近最新中文字幕大全免费视频 | 久久免费观看电影| 亚洲精品,欧美精品| 少妇被粗大猛烈的视频| 国产精品欧美亚洲77777| 亚洲色图综合在线观看| 精品久久久精品久久久| 人妻少妇偷人精品九色| 精品人妻一区二区三区麻豆| 精品少妇黑人巨大在线播放| 最近中文字幕2019免费版| 国产有黄有色有爽视频| 久久99精品国语久久久| 美女国产视频在线观看| 亚洲欧美成人精品一区二区| 日韩欧美一区视频在线观看| 久久av网站| 熟女电影av网| 五月伊人婷婷丁香| 国产一区二区在线观看av| 国产精品人妻久久久久久| 深夜精品福利| 2022亚洲国产成人精品| 爱豆传媒免费全集在线观看| 亚洲国产精品一区三区| 咕卡用的链子| 最近的中文字幕免费完整| 日本91视频免费播放| 深夜精品福利| 插逼视频在线观看| 永久免费av网站大全| av免费观看日本| 国产精品欧美亚洲77777| 亚洲国产最新在线播放| 国产成人欧美| 欧美97在线视频| av有码第一页| 日本黄大片高清| 国产精品国产三级国产av玫瑰| 久久99热这里只频精品6学生| 精品久久久精品久久久| 天堂俺去俺来也www色官网| 国产免费又黄又爽又色| 成人免费观看视频高清| 国产精品 国内视频| 国产成人91sexporn| 久久久久久久精品精品| 51国产日韩欧美| 久久婷婷青草| 欧美日本中文国产一区发布| 在线天堂中文资源库| 国产精品 国内视频| 女人精品久久久久毛片| 国产免费视频播放在线视频| 建设人人有责人人尽责人人享有的| 午夜福利乱码中文字幕| 大香蕉久久网| 日本色播在线视频| 丝袜美足系列| 久久精品国产亚洲av天美| 精品国产一区二区三区久久久樱花| 免费在线观看黄色视频的| 亚洲色图综合在线观看| 欧美成人午夜免费资源| 老司机影院毛片| 亚洲国产av新网站| 夜夜骑夜夜射夜夜干| 精品一区二区三卡| 大话2 男鬼变身卡| 久久人人爽人人爽人人片va| 99久久精品国产国产毛片| 热99久久久久精品小说推荐| 免费在线观看黄色视频的| 亚洲婷婷狠狠爱综合网| 又黄又爽又刺激的免费视频.| 男的添女的下面高潮视频| 黑人欧美特级aaaaaa片| 男人舔女人的私密视频| 人体艺术视频欧美日本| 国产成人精品婷婷| 日本欧美国产在线视频| 两个人免费观看高清视频| 欧美亚洲 丝袜 人妻 在线| www.熟女人妻精品国产 | 亚洲在久久综合| 亚洲色图 男人天堂 中文字幕 | 日本欧美国产在线视频| 各种免费的搞黄视频| 一二三四在线观看免费中文在 | 国产亚洲av片在线观看秒播厂| 男女边摸边吃奶| 男女边摸边吃奶| 久久精品久久久久久噜噜老黄| 美女xxoo啪啪120秒动态图| 人人妻人人添人人爽欧美一区卜| 欧美人与性动交α欧美精品济南到 | av一本久久久久| 制服诱惑二区| 国产极品粉嫩免费观看在线| 国产视频首页在线观看| 午夜免费鲁丝| 高清欧美精品videossex| 老司机影院成人| 一级黄片播放器| 蜜桃国产av成人99| 天堂俺去俺来也www色官网| 日韩视频在线欧美| 纯流量卡能插随身wifi吗| 黄色 视频免费看| 999精品在线视频| 日日啪夜夜爽| 久久av网站| 欧美日韩一区二区视频在线观看视频在线| 大陆偷拍与自拍| 精品人妻偷拍中文字幕| 成人毛片a级毛片在线播放| 国产不卡av网站在线观看| 日韩欧美一区视频在线观看| 国产av一区二区精品久久| 精品一品国产午夜福利视频| 午夜福利视频在线观看免费| 在线观看免费视频网站a站| 久久久国产欧美日韩av| 18禁国产床啪视频网站| 精品酒店卫生间| 天天躁夜夜躁狠狠久久av| 最新中文字幕久久久久| 日韩不卡一区二区三区视频在线| 国产一区二区三区av在线| 熟女电影av网| 90打野战视频偷拍视频| 永久免费av网站大全| 国产av精品麻豆| 18+在线观看网站| 成人漫画全彩无遮挡| 9热在线视频观看99| 夜夜爽夜夜爽视频| 国产成人精品在线电影| 亚洲色图 男人天堂 中文字幕 | 老熟女久久久| 美女国产高潮福利片在线看| 日韩中文字幕视频在线看片| www.色视频.com| 精品一区二区免费观看| 精品久久久久久电影网| 亚洲av国产av综合av卡| 亚洲av国产av综合av卡| 欧美日韩视频精品一区| 国产精品一国产av| 最新中文字幕久久久久| 又粗又硬又长又爽又黄的视频| 青青草视频在线视频观看| 超色免费av| 寂寞人妻少妇视频99o| 午夜福利网站1000一区二区三区| 天堂中文最新版在线下载| 久久这里只有精品19| 美女福利国产在线| 黑人欧美特级aaaaaa片| 久久综合国产亚洲精品| 桃花免费在线播放| 桃花免费在线播放| 日韩av免费高清视频| 欧美少妇被猛烈插入视频| 久久人人爽av亚洲精品天堂| 青春草亚洲视频在线观看| 天美传媒精品一区二区| 99九九在线精品视频| 丝袜美足系列| 王馨瑶露胸无遮挡在线观看| 国产成人aa在线观看| 国产一级毛片在线| 国产日韩欧美亚洲二区| 国产成人精品一,二区| 在线观看免费视频网站a站| 国产探花极品一区二区| 日本黄大片高清| 亚洲精品久久久久久婷婷小说| 亚洲熟女精品中文字幕| 91国产中文字幕| 22中文网久久字幕| 亚洲欧洲国产日韩| 亚洲一级一片aⅴ在线观看| 十分钟在线观看高清视频www| 日韩中文字幕视频在线看片| 亚洲天堂av无毛| 久久精品国产鲁丝片午夜精品| 免费不卡的大黄色大毛片视频在线观看| 国产白丝娇喘喷水9色精品| 久热久热在线精品观看| 国产精品三级大全| 国产精品成人在线| 亚洲成色77777| 美女中出高潮动态图| 18禁国产床啪视频网站| 欧美日韩一区二区视频在线观看视频在线| 亚洲熟女精品中文字幕| 日韩精品有码人妻一区| 精品一区在线观看国产| 最近手机中文字幕大全| 青春草国产在线视频| 80岁老熟妇乱子伦牲交| 熟女电影av网| 搡女人真爽免费视频火全软件| 中文字幕最新亚洲高清| www日本在线高清视频| 亚洲av免费高清在线观看| 国产无遮挡羞羞视频在线观看| 欧美 日韩 精品 国产| 国产成人午夜福利电影在线观看| 永久免费av网站大全| 我的女老师完整版在线观看| 一本久久精品| 乱码一卡2卡4卡精品| 精品国产一区二区三区久久久樱花| 亚洲天堂av无毛| 成人国产av品久久久| 亚洲一区二区三区欧美精品| 亚洲色图 男人天堂 中文字幕 | 免费在线观看黄色视频的| 中文字幕最新亚洲高清| 黄网站色视频无遮挡免费观看| 成年女人在线观看亚洲视频| 亚洲av日韩在线播放| 中文精品一卡2卡3卡4更新| 人妻人人澡人人爽人人| 女人被躁到高潮嗷嗷叫费观| 国产一区二区激情短视频 | 国产av码专区亚洲av| freevideosex欧美| 亚洲av国产av综合av卡| 如日韩欧美国产精品一区二区三区| 亚洲内射少妇av| 色婷婷久久久亚洲欧美| 日韩视频在线欧美| 欧美精品一区二区大全| 多毛熟女@视频| 人妻系列 视频| 午夜福利乱码中文字幕| 午夜福利网站1000一区二区三区| 久久综合国产亚洲精品| 飞空精品影院首页| 国产成人精品一,二区| 少妇高潮的动态图| 成人亚洲欧美一区二区av| 天天操日日干夜夜撸| 免费高清在线观看日韩| 国产永久视频网站| 亚洲欧洲国产日韩| 久久久久视频综合| 青青草视频在线视频观看| 又黄又粗又硬又大视频| 97精品久久久久久久久久精品| 中文字幕制服av| 国产69精品久久久久777片| 亚洲精品中文字幕在线视频| 日本91视频免费播放| videos熟女内射| 国产免费又黄又爽又色| 色5月婷婷丁香| 在线看a的网站| 亚洲熟女精品中文字幕| 草草在线视频免费看| 国产视频首页在线观看| 国产成人91sexporn| 一二三四中文在线观看免费高清| 亚洲美女搞黄在线观看| 亚洲经典国产精华液单| 一区二区三区乱码不卡18| 久久久精品区二区三区| 欧美成人精品欧美一级黄| 国产精品人妻久久久久久| 看免费成人av毛片| 久久久久久久久久久久大奶| 日韩av免费高清视频| 久久久久网色| 18禁国产床啪视频网站| 只有这里有精品99| 色视频在线一区二区三区| 香蕉精品网在线| 成人午夜精彩视频在线观看| 国产男女超爽视频在线观看| 菩萨蛮人人尽说江南好唐韦庄| 成人国产av品久久久| 国产精品.久久久| 超色免费av| 亚洲av成人精品一二三区| 在线观看国产h片| 在线观看美女被高潮喷水网站| 欧美日韩亚洲高清精品| 性高湖久久久久久久久免费观看| 国产麻豆69| 国产精品嫩草影院av在线观看| 免费人成在线观看视频色| 满18在线观看网站| 高清av免费在线| 亚洲国产成人一精品久久久| 国产一区二区三区综合在线观看 | 国产欧美日韩综合在线一区二区| 免费观看无遮挡的男女| 老熟女久久久| 你懂的网址亚洲精品在线观看| 亚洲欧洲精品一区二区精品久久久 | av一本久久久久| 免费人成在线观看视频色| 在线观看三级黄色| 国产片特级美女逼逼视频| 日韩在线高清观看一区二区三区| 日韩欧美一区视频在线观看| 亚洲av国产av综合av卡| 亚洲精品久久午夜乱码| 久久精品国产a三级三级三级| 精品视频人人做人人爽| 亚洲美女黄色视频免费看| 亚洲国产最新在线播放| 亚洲精品成人av观看孕妇| 精品国产一区二区三区四区第35| 欧美日韩综合久久久久久| 亚洲精品国产色婷婷电影| 中文字幕亚洲精品专区| 久久99一区二区三区| 久久青草综合色| 久久久国产欧美日韩av| 久久99精品国语久久久| 国产免费一区二区三区四区乱码| 国产淫语在线视频| 亚洲精品自拍成人| 免费观看性生交大片5| 乱人伦中国视频| 国产免费现黄频在线看| 少妇的逼好多水| 另类精品久久| videosex国产| 黑人猛操日本美女一级片| 久久久久久人妻| 男的添女的下面高潮视频| 啦啦啦中文免费视频观看日本| 青春草视频在线免费观看| av福利片在线| 国产 精品1| 免费av不卡在线播放| 免费观看a级毛片全部| 9191精品国产免费久久| 热99久久久久精品小说推荐| 精品国产一区二区久久| 性色avwww在线观看| 人人妻人人爽人人添夜夜欢视频| 亚洲综合色网址| 热99国产精品久久久久久7| 熟女av电影| 国产激情久久老熟女| 十八禁高潮呻吟视频| 欧美日韩成人在线一区二区|