• <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)計分析
    在线观看免费日韩欧美大片| 美女xxoo啪啪120秒动态图| 美女主播在线视频| 亚洲国产看品久久| 久久久久久久大尺度免费视频| 一级黄片播放器| 一级毛片我不卡| 人妻少妇偷人精品九色| 婷婷色av中文字幕| 欧美精品高潮呻吟av久久| 极品少妇高潮喷水抽搐| 丝袜在线中文字幕| 精品99又大又爽又粗少妇毛片| 久久久久久久大尺度免费视频| 麻豆精品久久久久久蜜桃| 午夜福利影视在线免费观看| 国产精品欧美亚洲77777| 大码成人一级视频| 久久久久国产一级毛片高清牌| 亚洲内射少妇av| 熟女少妇亚洲综合色aaa.| 亚洲综合色惰| 国产亚洲午夜精品一区二区久久| 在线观看美女被高潮喷水网站| 国产日韩一区二区三区精品不卡| 亚洲av电影在线观看一区二区三区| 99久久中文字幕三级久久日本| 五月伊人婷婷丁香| 国产精品国产三级专区第一集| 哪个播放器可以免费观看大片| 日韩av在线免费看完整版不卡| 1024视频免费在线观看| 亚洲一区二区三区欧美精品| 国产精品.久久久| 女性生殖器流出的白浆| 国产免费福利视频在线观看| 一级片'在线观看视频| 又黄又粗又硬又大视频| 亚洲第一av免费看| 纵有疾风起免费观看全集完整版| 自拍欧美九色日韩亚洲蝌蚪91| 午夜福利影视在线免费观看| 国产免费又黄又爽又色| av网站免费在线观看视频| 国产精品人妻久久久影院| 日韩欧美精品免费久久| 天天躁狠狠躁夜夜躁狠狠躁| 满18在线观看网站| av又黄又爽大尺度在线免费看| av电影中文网址| 伊人久久国产一区二区| 制服丝袜香蕉在线| av又黄又爽大尺度在线免费看| 日本欧美国产在线视频| 波多野结衣一区麻豆| 看免费成人av毛片| 亚洲欧美一区二区三区黑人 | 一二三四中文在线观看免费高清| 国产精品嫩草影院av在线观看| 亚洲欧美日韩另类电影网站| 国产精品香港三级国产av潘金莲 | videosex国产| 亚洲av男天堂| 不卡av一区二区三区| av在线老鸭窝| 国产成人精品久久久久久| 各种免费的搞黄视频| 亚洲综合精品二区| 天天影视国产精品| 亚洲av日韩在线播放| av又黄又爽大尺度在线免费看| 2021少妇久久久久久久久久久| 午夜福利视频在线观看免费| 男女下面插进去视频免费观看| 久久久国产一区二区| www日本在线高清视频| 一级爰片在线观看| 伊人亚洲综合成人网| 有码 亚洲区| 曰老女人黄片| 熟女av电影| 国产精品无大码| 啦啦啦中文免费视频观看日本| 黄色 视频免费看| 欧美日本中文国产一区发布| 日韩中字成人| 国产成人a∨麻豆精品| 18禁观看日本| 黑丝袜美女国产一区| 日本vs欧美在线观看视频| 亚洲精品国产av成人精品| 欧美精品高潮呻吟av久久| 成人毛片60女人毛片免费| 有码 亚洲区| 秋霞在线观看毛片| 精品人妻熟女毛片av久久网站| 波多野结衣一区麻豆| 啦啦啦视频在线资源免费观看| 熟妇人妻不卡中文字幕| 久久国产精品大桥未久av| 一级毛片 在线播放| 夜夜骑夜夜射夜夜干| 国产色婷婷99| av在线观看视频网站免费| 精品少妇久久久久久888优播| 国产精品蜜桃在线观看| 欧美日韩综合久久久久久| 18禁裸乳无遮挡动漫免费视频| 中文乱码字字幕精品一区二区三区| 国产淫语在线视频| 男女边摸边吃奶| 国产人伦9x9x在线观看 | 精品一区二区免费观看| 久久久久精品性色| 亚洲成av片中文字幕在线观看 | 三级国产精品片| 欧美黄色片欧美黄色片| 日韩制服骚丝袜av| 亚洲国产精品成人久久小说| 国产成人一区二区在线| 人妻一区二区av| 亚洲一码二码三码区别大吗| 久久热在线av| 韩国精品一区二区三区| 国产福利在线免费观看视频| 激情五月婷婷亚洲| 三上悠亚av全集在线观看| 免费高清在线观看日韩| 国产男人的电影天堂91| 午夜福利视频在线观看免费| av网站免费在线观看视频| 久久影院123| 久久久久久免费高清国产稀缺| 九草在线视频观看| 大片电影免费在线观看免费| 另类精品久久| 美女脱内裤让男人舔精品视频| 啦啦啦在线观看免费高清www| 亚洲精品成人av观看孕妇| 精品人妻熟女毛片av久久网站| 亚洲国产av新网站| 国产精品免费视频内射| 伦理电影大哥的女人| 日韩大片免费观看网站| 国产精品国产三级国产专区5o| 看免费成人av毛片| 一级片'在线观看视频| 亚洲三区欧美一区| 欧美最新免费一区二区三区| 97在线视频观看| 亚洲av综合色区一区| 90打野战视频偷拍视频| 有码 亚洲区| 国产精品嫩草影院av在线观看| 亚洲少妇的诱惑av| 狠狠精品人妻久久久久久综合| 熟女av电影| 日本-黄色视频高清免费观看| 亚洲精品,欧美精品| 七月丁香在线播放| 欧美国产精品一级二级三级| 天美传媒精品一区二区| 久久ye,这里只有精品| 美女xxoo啪啪120秒动态图| 久久 成人 亚洲| videossex国产| 国语对白做爰xxxⅹ性视频网站| 亚洲国产精品成人久久小说| 久久人人爽人人片av| 午夜福利视频在线观看免费| 巨乳人妻的诱惑在线观看| 伦理电影大哥的女人| 中文欧美无线码| 在线观看www视频免费| 一边亲一边摸免费视频| 热re99久久国产66热| 婷婷成人精品国产| 热99久久久久精品小说推荐| 国产成人精品无人区| 国产精品亚洲av一区麻豆 | 亚洲熟女精品中文字幕| 丰满乱子伦码专区| 亚洲国产色片| 丁香六月天网| 90打野战视频偷拍视频| 免费播放大片免费观看视频在线观看| 免费久久久久久久精品成人欧美视频| 91精品伊人久久大香线蕉| 女的被弄到高潮叫床怎么办| 亚洲人成网站在线观看播放| 香蕉丝袜av| 精品国产乱码久久久久久小说| 久久久久人妻精品一区果冻| 曰老女人黄片| 亚洲国产毛片av蜜桃av| 久久精品人人爽人人爽视色| 久久久久精品久久久久真实原创| 日日爽夜夜爽网站| 国产精品蜜桃在线观看| 宅男免费午夜| 国产成人精品无人区| 国产精品人妻久久久影院| 咕卡用的链子| 狠狠精品人妻久久久久久综合| 日本黄色日本黄色录像| av免费观看日本| 校园人妻丝袜中文字幕| 少妇人妻精品综合一区二区| 巨乳人妻的诱惑在线观看| 永久网站在线| videos熟女内射| 国产免费现黄频在线看| 日韩精品有码人妻一区| 日本爱情动作片www.在线观看| 秋霞在线观看毛片| 美女中出高潮动态图| 90打野战视频偷拍视频| 街头女战士在线观看网站| 桃花免费在线播放| 国精品久久久久久国模美| 曰老女人黄片| 国产在线一区二区三区精| 国产一区二区三区综合在线观看| 久久久久久久久久久久大奶| 日韩大片免费观看网站| 久久国产亚洲av麻豆专区| 少妇 在线观看| 国产亚洲av片在线观看秒播厂| 在线精品无人区一区二区三| 亚洲国产精品一区二区三区在线| 日韩精品有码人妻一区| 一级毛片电影观看| 亚洲欧美中文字幕日韩二区| 中文字幕制服av| 一区福利在线观看| 欧美xxⅹ黑人| 极品少妇高潮喷水抽搐| 久久久精品免费免费高清| 色网站视频免费| 免费在线观看黄色视频的| 制服人妻中文乱码| 亚洲欧美一区二区三区国产| 涩涩av久久男人的天堂| 国产成人一区二区在线| 亚洲av电影在线进入| 久久久久久久久久久久大奶| 最黄视频免费看| 中文精品一卡2卡3卡4更新| 久久综合国产亚洲精品| 国产在视频线精品| 人妻系列 视频| 免费少妇av软件| 高清视频免费观看一区二区| 免费女性裸体啪啪无遮挡网站| 国产精品无大码| xxx大片免费视频| 久久久国产一区二区| 亚洲av电影在线观看一区二区三区| 欧美日韩视频高清一区二区三区二| 欧美日韩精品成人综合77777| 国产欧美日韩综合在线一区二区| 亚洲精品视频女| 美女主播在线视频| 国产精品免费大片| 日日爽夜夜爽网站| 少妇人妻精品综合一区二区| 免费观看av网站的网址| 精品少妇内射三级| 老汉色av国产亚洲站长工具| 男女国产视频网站| 精品国产一区二区三区四区第35| 国产成人91sexporn| 在线观看www视频免费| 满18在线观看网站| 亚洲精品国产av成人精品| 久久精品国产亚洲av涩爱| 有码 亚洲区| 十八禁高潮呻吟视频| 超色免费av| www日本在线高清视频| 成人免费观看视频高清| 日韩成人av中文字幕在线观看| 国产精品香港三级国产av潘金莲 | 亚洲人成电影观看| 美女国产高潮福利片在线看| 日本免费在线观看一区| 午夜福利视频精品| 岛国毛片在线播放| 久久久久久久久免费视频了| 成人午夜精彩视频在线观看| 国产精品人妻久久久影院| 18禁观看日本| 伊人亚洲综合成人网| 婷婷色综合大香蕉| 国产精品一国产av| 久久久久精品久久久久真实原创| 搡老乐熟女国产| 一边摸一边做爽爽视频免费| 亚洲欧美色中文字幕在线| 亚洲av在线观看美女高潮| av视频免费观看在线观看| 搡女人真爽免费视频火全软件| 性色av一级| 欧美精品人与动牲交sv欧美| av不卡在线播放| 国产白丝娇喘喷水9色精品| 中国国产av一级| 日韩电影二区| 在线观看三级黄色| 欧美av亚洲av综合av国产av | 91在线精品国自产拍蜜月| 亚洲视频免费观看视频| 精品一区在线观看国产| 精品视频人人做人人爽| 又粗又硬又长又爽又黄的视频| 欧美成人午夜免费资源| 日韩中字成人| 亚洲精品在线美女| 久久99精品国语久久久| 精品少妇一区二区三区视频日本电影 | 性色av一级| 久久久亚洲精品成人影院| 国语对白做爰xxxⅹ性视频网站| 欧美激情极品国产一区二区三区| 免费观看性生交大片5| a级片在线免费高清观看视频| 大片电影免费在线观看免费| 国产色婷婷99| 久久久国产一区二区| 在线观看人妻少妇| 18禁裸乳无遮挡动漫免费视频| 国产深夜福利视频在线观看| 久久久精品国产亚洲av高清涩受| 九九爱精品视频在线观看| 波多野结衣一区麻豆| 国产精品一二三区在线看| 一级,二级,三级黄色视频| 欧美 亚洲 国产 日韩一| 亚洲成人av在线免费| 一区二区三区四区激情视频| 不卡视频在线观看欧美| 综合色丁香网| 日韩,欧美,国产一区二区三区| 青春草亚洲视频在线观看| 久久婷婷青草| 97在线视频观看| 久久精品熟女亚洲av麻豆精品| 美女脱内裤让男人舔精品视频| 91精品伊人久久大香线蕉| 国产一区二区三区综合在线观看| 看免费av毛片| 国产成人精品福利久久| 久久青草综合色| 99re6热这里在线精品视频| 久久影院123| 国产成人午夜福利电影在线观看| 一级毛片电影观看| 国产成人午夜福利电影在线观看| 中国国产av一级| 18+在线观看网站| 国产男女内射视频| 亚洲精品久久成人aⅴ小说| xxx大片免费视频| 久久久久久久亚洲中文字幕| 一本色道久久久久久精品综合| 十八禁高潮呻吟视频| 国产精品秋霞免费鲁丝片| 国产精品久久久久久久久免| 精品国产乱码久久久久久小说| 美女视频免费永久观看网站| 中文字幕人妻熟女乱码| www.自偷自拍.com| 性少妇av在线| 国产精品麻豆人妻色哟哟久久| 日日撸夜夜添| 天天影视国产精品| 午夜久久久在线观看| 亚洲精品久久成人aⅴ小说| 国产精品免费大片| 免费看不卡的av| 97在线视频观看| 久久久精品免费免费高清| 色哟哟·www| 人妻 亚洲 视频| 免费av中文字幕在线| 亚洲精品久久成人aⅴ小说| 日本猛色少妇xxxxx猛交久久| 亚洲第一av免费看| 夫妻午夜视频| 欧美黄色片欧美黄色片| 啦啦啦中文免费视频观看日本| 欧美日本中文国产一区发布| 国产1区2区3区精品| 女人久久www免费人成看片| 啦啦啦在线免费观看视频4| 日本黄色日本黄色录像| av卡一久久| 在线天堂中文资源库| 黄色毛片三级朝国网站| 大陆偷拍与自拍| 亚洲精品国产av成人精品| 免费观看av网站的网址| 人妻少妇偷人精品九色| 这个男人来自地球电影免费观看 | a级片在线免费高清观看视频| 制服人妻中文乱码| 国产一区二区三区综合在线观看| 国产精品久久久久久av不卡| 你懂的网址亚洲精品在线观看| 午夜老司机福利剧场| 国产av码专区亚洲av| 欧美人与性动交α欧美软件| 纵有疾风起免费观看全集完整版| 秋霞在线观看毛片| 国产熟女午夜一区二区三区| 国产在线一区二区三区精| 国产成人免费观看mmmm| 天天影视国产精品| av一本久久久久| 黄片小视频在线播放| 色视频在线一区二区三区| 欧美日本中文国产一区发布| 久久午夜福利片| 中文字幕另类日韩欧美亚洲嫩草| 丝袜脚勾引网站| 男女啪啪激烈高潮av片| 亚洲成人一二三区av| 国产欧美日韩综合在线一区二区| videos熟女内射| 大码成人一级视频| 国产不卡av网站在线观看| 中国国产av一级| 精品人妻熟女毛片av久久网站| 国产1区2区3区精品| 欧美日韩综合久久久久久| av不卡在线播放| 欧美精品人与动牲交sv欧美| 丰满饥渴人妻一区二区三| 制服诱惑二区| 久久精品国产鲁丝片午夜精品| 91精品伊人久久大香线蕉| 日韩熟女老妇一区二区性免费视频| 嫩草影院入口| 亚洲国产精品一区二区三区在线| 一级黄片播放器| 免费观看a级毛片全部| 丝瓜视频免费看黄片| 亚洲精品乱久久久久久| 少妇被粗大的猛进出69影院| 看免费av毛片| 人体艺术视频欧美日本| 欧美精品亚洲一区二区| 一区在线观看完整版| 成年av动漫网址| 亚洲一级一片aⅴ在线观看| 国产精品香港三级国产av潘金莲 | 久久久国产一区二区| 国产不卡av网站在线观看| 久久99热这里只频精品6学生| 性高湖久久久久久久久免费观看| 99久国产av精品国产电影| 男女边摸边吃奶| 婷婷色av中文字幕| 日日撸夜夜添| 久久婷婷青草| 日本爱情动作片www.在线观看| 90打野战视频偷拍视频| 美女国产高潮福利片在线看| 黄片无遮挡物在线观看| 伊人久久大香线蕉亚洲五| 女人精品久久久久毛片| 少妇的逼水好多| 久久青草综合色| 国产精品av久久久久免费| 美女国产高潮福利片在线看| 国产一区二区三区av在线| 乱人伦中国视频| 精品午夜福利在线看| 青春草视频在线免费观看| 搡老乐熟女国产| 美女高潮到喷水免费观看| 夜夜骑夜夜射夜夜干| 亚洲激情五月婷婷啪啪| 99久久综合免费| 午夜免费男女啪啪视频观看| 超色免费av| 高清欧美精品videossex| 18+在线观看网站| 精品午夜福利在线看| 欧美在线黄色| 色哟哟·www| 亚洲视频免费观看视频| 国产精品 欧美亚洲| 久久久a久久爽久久v久久| 国产精品秋霞免费鲁丝片| 亚洲美女搞黄在线观看| 亚洲欧洲国产日韩| 18禁观看日本| 日产精品乱码卡一卡2卡三| 亚洲五月色婷婷综合| 久久久国产精品麻豆| 黄片播放在线免费| 国产探花极品一区二区| 日韩欧美精品免费久久| 这个男人来自地球电影免费观看 | 伦理电影大哥的女人| 五月天丁香电影| 国产精品无大码| 精品亚洲成国产av| 伦理电影免费视频| 成人免费观看视频高清| 人妻少妇偷人精品九色| 26uuu在线亚洲综合色| 免费少妇av软件| 欧美日韩视频高清一区二区三区二| 精品国产一区二区三区久久久樱花| 一二三四在线观看免费中文在| 婷婷色综合大香蕉| 亚洲国产看品久久| 国精品久久久久久国模美| 最近最新中文字幕大全免费视频 | 国产精品免费视频内射| 日韩不卡一区二区三区视频在线| 精品视频人人做人人爽| 日日摸夜夜添夜夜爱| 永久网站在线| 飞空精品影院首页| 一级片'在线观看视频| 国产在线视频一区二区| 亚洲国产色片| 99热国产这里只有精品6| 一级爰片在线观看| 亚洲情色 制服丝袜| 一级毛片黄色毛片免费观看视频| 国产日韩欧美视频二区| 99热网站在线观看| 免费看不卡的av| 秋霞在线观看毛片| 精品国产一区二区久久| 在线观看免费视频网站a站| 日日摸夜夜添夜夜爱| 成人毛片a级毛片在线播放| 狠狠精品人妻久久久久久综合| 日本猛色少妇xxxxx猛交久久| av有码第一页| 久久久久国产精品人妻一区二区| 精品久久久久久电影网| 考比视频在线观看| 午夜91福利影院| 老鸭窝网址在线观看| 天天躁夜夜躁狠狠躁躁| 男女国产视频网站| 99国产综合亚洲精品| 韩国av在线不卡| 成人国产av品久久久| av.在线天堂| 激情视频va一区二区三区| 日韩三级伦理在线观看| 欧美日本中文国产一区发布| a 毛片基地| 国产精品秋霞免费鲁丝片| 1024香蕉在线观看| 嫩草影院入口| 日产精品乱码卡一卡2卡三| 久久久a久久爽久久v久久| 天美传媒精品一区二区| 精品福利永久在线观看| 天天躁夜夜躁狠狠躁躁| 国产人伦9x9x在线观看 | 国产日韩欧美视频二区| 18禁裸乳无遮挡动漫免费视频| 久久精品国产亚洲av涩爱| 亚洲国产色片| 国产一区二区在线观看av| 亚洲,欧美精品.| 亚洲av免费高清在线观看| 少妇被粗大猛烈的视频| 少妇被粗大的猛进出69影院| 国语对白做爰xxxⅹ性视频网站| 如日韩欧美国产精品一区二区三区| 老熟女久久久| 久久久国产一区二区| 久久av网站| 超碰成人久久| 最新中文字幕久久久久| 亚洲第一av免费看| 久久午夜综合久久蜜桃| 女性生殖器流出的白浆| 国产精品免费视频内射| 最近最新中文字幕免费大全7| 不卡视频在线观看欧美| 国产成人aa在线观看| 国产97色在线日韩免费| 日本黄色日本黄色录像| 一二三四在线观看免费中文在| 七月丁香在线播放| 日本黄色日本黄色录像| 久久久久国产精品人妻一区二区| 日韩大片免费观看网站| 黄色一级大片看看| 最近的中文字幕免费完整| 国产精品亚洲av一区麻豆 | 亚洲欧洲国产日韩| 亚洲国产精品999| 美女国产高潮福利片在线看| 亚洲伊人久久精品综合| 国产探花极品一区二区| 97精品久久久久久久久久精品| 成人国语在线视频| 久久久久国产网址| 欧美成人午夜精品| 久久久国产一区二区| 国产av精品麻豆| 亚洲精品中文字幕在线视频| 国产精品欧美亚洲77777|