• <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)計分析
    99九九线精品视频在线观看视频| 熟女人妻精品中文字幕| 日韩人妻高清精品专区| a级一级毛片免费在线观看| 婷婷色综合www| 菩萨蛮人人尽说江南好唐韦庄| 网址你懂的国产日韩在线| 国产大屁股一区二区在线视频| 99久久精品国产国产毛片| 久久久久久久久大av| 国产熟女欧美一区二区| 狂野欧美白嫩少妇大欣赏| a级毛片免费高清观看在线播放| 91aial.com中文字幕在线观看| 欧美高清成人免费视频www| 国语对白做爰xxxⅹ性视频网站| 最黄视频免费看| 深夜a级毛片| 亚洲国产精品专区欧美| av在线观看视频网站免费| 黄色日韩在线| 国产免费又黄又爽又色| 黄色视频在线播放观看不卡| 偷拍熟女少妇极品色| 搡女人真爽免费视频火全软件| 精品熟女少妇av免费看| 久久久久国产网址| 精品国产露脸久久av麻豆| 亚洲激情五月婷婷啪啪| 国产精品99久久久久久久久| 三级国产精品片| 高清不卡的av网站| 国产精品一区二区在线不卡| 亚洲怡红院男人天堂| 亚洲精品亚洲一区二区| 三级经典国产精品| 91精品一卡2卡3卡4卡| av线在线观看网站| 日韩一区二区三区影片| 久久久色成人| 午夜福利视频精品| 午夜免费鲁丝| 爱豆传媒免费全集在线观看| 精品人妻视频免费看| 寂寞人妻少妇视频99o| 成人高潮视频无遮挡免费网站| 国产淫语在线视频| 亚洲av在线观看美女高潮| a级毛色黄片| 久久久久精品性色| 18禁在线无遮挡免费观看视频| 激情 狠狠 欧美| 爱豆传媒免费全集在线观看| 国产伦精品一区二区三区视频9| 国产黄频视频在线观看| 人妻夜夜爽99麻豆av| 不卡视频在线观看欧美| 丝袜喷水一区| 日韩亚洲欧美综合| 国产在线免费精品| 日韩亚洲欧美综合| 成人亚洲欧美一区二区av| 黄片wwwwww| 亚洲成人av在线免费| 人妻 亚洲 视频| 人人妻人人添人人爽欧美一区卜 | 99精国产麻豆久久婷婷| 亚洲精品一二三| 婷婷色av中文字幕| 丝袜脚勾引网站| 能在线免费看毛片的网站| 久久影院123| 国产精品免费大片| 最近最新中文字幕免费大全7| 亚洲精品,欧美精品| 高清黄色对白视频在线免费看 | 日韩亚洲欧美综合| 一二三四中文在线观看免费高清| 激情五月婷婷亚洲| 国产精品免费大片| 亚州av有码| 亚洲av.av天堂| 欧美精品一区二区免费开放| 日韩人妻高清精品专区| 欧美日韩综合久久久久久| 黄色配什么色好看| 亚洲在久久综合| 中国美白少妇内射xxxbb| 国产 一区 欧美 日韩| 丝袜喷水一区| 亚洲图色成人| 中文字幕亚洲精品专区| 少妇 在线观看| 国产亚洲精品久久久com| 亚洲aⅴ乱码一区二区在线播放| 国产精品欧美亚洲77777| 自拍欧美九色日韩亚洲蝌蚪91 | 亚洲精品自拍成人| 日韩大片免费观看网站| 国产成人免费观看mmmm| 亚洲欧美日韩无卡精品| 国产av精品麻豆| 新久久久久国产一级毛片| 一本久久精品| 欧美亚洲 丝袜 人妻 在线| 亚洲在久久综合| av女优亚洲男人天堂| 18禁在线无遮挡免费观看视频| 国产精品女同一区二区软件| 欧美少妇被猛烈插入视频| 中国三级夫妇交换| 少妇人妻久久综合中文| 亚洲精品日韩av片在线观看| 日产精品乱码卡一卡2卡三| 身体一侧抽搐| 国产高清不卡午夜福利| 岛国毛片在线播放| 成人一区二区视频在线观看| 美女视频免费永久观看网站| 国产中年淑女户外野战色| 亚洲欧美日韩卡通动漫| av又黄又爽大尺度在线免费看| 日韩精品有码人妻一区| 国产精品国产三级专区第一集| 精品人妻偷拍中文字幕| 天堂中文最新版在线下载| 岛国毛片在线播放| 欧美日本视频| 国产淫语在线视频| 老女人水多毛片| 久久久久久人妻| 熟妇人妻不卡中文字幕| 美女视频免费永久观看网站| 三级国产精品片| 在线免费十八禁| 国产免费一区二区三区四区乱码| 一级二级三级毛片免费看| 夫妻午夜视频| 国产成人精品婷婷| 精品国产三级普通话版| 2022亚洲国产成人精品| 中文天堂在线官网| 国产亚洲5aaaaa淫片| 日本与韩国留学比较| 亚洲美女黄色视频免费看| 91aial.com中文字幕在线观看| 久久精品人妻少妇| av在线观看视频网站免费| 2022亚洲国产成人精品| 哪个播放器可以免费观看大片| 2018国产大陆天天弄谢| 91久久精品国产一区二区三区| 国产亚洲最大av| 日韩制服骚丝袜av| 久久久久久久久大av| 黄片无遮挡物在线观看| 国产精品久久久久久久电影| 黑丝袜美女国产一区| 老司机影院毛片| 国产高清不卡午夜福利| 三级经典国产精品| www.av在线官网国产| 国产精品99久久久久久久久| 国产成人一区二区在线| 内射极品少妇av片p| 97超碰精品成人国产| 国产午夜精品一二区理论片| 夜夜爽夜夜爽视频| 人人妻人人看人人澡| 边亲边吃奶的免费视频| 黄片无遮挡物在线观看| 亚洲精品久久久久久婷婷小说| 亚洲成人中文字幕在线播放| 欧美日韩精品成人综合77777| 精品久久久久久久久亚洲| 少妇丰满av| 久久精品国产自在天天线| 99久久精品热视频| 五月玫瑰六月丁香| 又粗又硬又长又爽又黄的视频| 国产亚洲91精品色在线| 欧美另类一区| 激情五月婷婷亚洲| 青春草亚洲视频在线观看| 少妇裸体淫交视频免费看高清| 视频中文字幕在线观看| 校园人妻丝袜中文字幕| 成人漫画全彩无遮挡| 搡老乐熟女国产| 中国国产av一级| 亚洲色图综合在线观看| 国精品久久久久久国模美| 91久久精品国产一区二区成人| 亚洲色图综合在线观看| 亚洲av电影在线观看一区二区三区| 久久影院123| 天堂8中文在线网| 成人18禁高潮啪啪吃奶动态图 | 在线观看国产h片| 久久久国产一区二区| 日韩强制内射视频| 青青草视频在线视频观看| 一区二区三区四区激情视频| 卡戴珊不雅视频在线播放| 美女内射精品一级片tv| 亚洲va在线va天堂va国产| 色婷婷久久久亚洲欧美| 亚洲欧美一区二区三区国产| 99热6这里只有精品| 99热6这里只有精品| a 毛片基地| 熟女人妻精品中文字幕| 水蜜桃什么品种好| 97精品久久久久久久久久精品| 国产片特级美女逼逼视频| 水蜜桃什么品种好| 亚洲av综合色区一区| 中文字幕亚洲精品专区| 高清毛片免费看| 精品99又大又爽又粗少妇毛片| 免费黄网站久久成人精品| 国产深夜福利视频在线观看| 精品亚洲成国产av| 国产成人一区二区在线| 亚洲成人av在线免费| 午夜福利高清视频| 久久久午夜欧美精品| 亚洲伊人久久精品综合| 亚洲精品国产色婷婷电影| av卡一久久| 国产淫片久久久久久久久| 国产v大片淫在线免费观看| 国产亚洲av片在线观看秒播厂| 色5月婷婷丁香| .国产精品久久| 欧美精品人与动牲交sv欧美| 黄片无遮挡物在线观看| 99热网站在线观看| 久久影院123| 国产亚洲精品久久久com| 美女视频免费永久观看网站| 亚洲av成人精品一区久久| 国产一区二区三区综合在线观看 | 国产av精品麻豆| 亚洲精品自拍成人| 色婷婷av一区二区三区视频| 亚洲精品国产色婷婷电影| 成人无遮挡网站| 国产精品国产三级国产专区5o| 久久久久久久久久久丰满| 国产精品人妻久久久久久| 国产免费一级a男人的天堂| 菩萨蛮人人尽说江南好唐韦庄| 亚洲色图综合在线观看| 乱系列少妇在线播放| 国产精品伦人一区二区| 老女人水多毛片| 一级a做视频免费观看| 日日摸夜夜添夜夜爱| 好男人视频免费观看在线| 日韩 亚洲 欧美在线| 久久久久久人妻| kizo精华| 欧美xxⅹ黑人| 少妇丰满av| 国国产精品蜜臀av免费| av卡一久久| 少妇 在线观看| 亚洲人成网站高清观看| 亚洲精华国产精华液的使用体验| 成人亚洲欧美一区二区av| 免费高清在线观看视频在线观看| 激情 狠狠 欧美| 日本黄色片子视频| 看免费成人av毛片| 欧美xxⅹ黑人| 26uuu在线亚洲综合色| 国产高清国产精品国产三级 | 美女中出高潮动态图| 亚洲欧美成人精品一区二区| 日本免费在线观看一区| 亚洲av欧美aⅴ国产| 草草在线视频免费看| 国产男女内射视频| 国产精品精品国产色婷婷| 亚洲成人手机| 久久精品夜色国产| av在线老鸭窝| 十八禁网站网址无遮挡 | 最后的刺客免费高清国语| 日韩国内少妇激情av| 青春草亚洲视频在线观看| 狠狠精品人妻久久久久久综合| 免费观看无遮挡的男女| 一本一本综合久久| 国产精品伦人一区二区| 人妻 亚洲 视频| 日本免费在线观看一区| 国产精品秋霞免费鲁丝片| 国产高清有码在线观看视频| 亚洲欧美日韩另类电影网站 | 欧美日韩综合久久久久久| 哪个播放器可以免费观看大片| 国产视频首页在线观看| 日韩亚洲欧美综合| 久久久久久久久久久丰满| 大香蕉久久网| 自拍偷自拍亚洲精品老妇| 高清av免费在线| 我要看黄色一级片免费的| 精品99又大又爽又粗少妇毛片| av在线app专区| 尤物成人国产欧美一区二区三区| av线在线观看网站| 亚洲美女搞黄在线观看| 中文字幕免费在线视频6| 极品教师在线视频| 干丝袜人妻中文字幕| 国产精品一二三区在线看| 国产淫片久久久久久久久| 久久国内精品自在自线图片| 2021少妇久久久久久久久久久| 丰满少妇做爰视频| 麻豆成人午夜福利视频| 日韩亚洲欧美综合| 99久国产av精品国产电影| 高清欧美精品videossex| 我要看黄色一级片免费的| 亚洲欧美清纯卡通| 一区二区三区乱码不卡18| videossex国产| 毛片女人毛片| 在线 av 中文字幕| 熟女人妻精品中文字幕| 在线 av 中文字幕| 精品人妻一区二区三区麻豆| 黄色一级大片看看| 亚洲欧美精品专区久久| 伦理电影免费视频| 黄色视频在线播放观看不卡| 国产日韩欧美在线精品| 国产男女超爽视频在线观看| 99热国产这里只有精品6| 国产高清有码在线观看视频| 日韩一区二区三区影片| 成人午夜精彩视频在线观看| 欧美日韩一区二区视频在线观看视频在线| 少妇丰满av| 亚洲色图综合在线观看| 热re99久久精品国产66热6| 在线播放无遮挡| 日日啪夜夜撸| 一边亲一边摸免费视频| 高清欧美精品videossex| 精品亚洲乱码少妇综合久久| 人妻一区二区av| 免费人妻精品一区二区三区视频| a级毛片免费高清观看在线播放| 免费观看的影片在线观看| 蜜桃久久精品国产亚洲av| 色综合色国产| 久久精品久久久久久噜噜老黄| 啦啦啦中文免费视频观看日本| 欧美xxxx性猛交bbbb| 亚洲国产精品国产精品| 人人妻人人澡人人爽人人夜夜| 日本午夜av视频| 有码 亚洲区| 大片免费播放器 马上看| 老熟女久久久| 熟女电影av网| 啦啦啦在线观看免费高清www| 在线免费十八禁| 久久久精品免费免费高清| 亚洲电影在线观看av| 看非洲黑人一级黄片| 人人妻人人添人人爽欧美一区卜 | av卡一久久| 亚洲av男天堂| 91久久精品国产一区二区成人| 18禁动态无遮挡网站| 男女无遮挡免费网站观看| 交换朋友夫妻互换小说| 有码 亚洲区| 热99国产精品久久久久久7| 老熟女久久久| 国产白丝娇喘喷水9色精品| 国产伦精品一区二区三区视频9| 美女内射精品一级片tv| 国产高清三级在线| 夜夜看夜夜爽夜夜摸| 日韩大片免费观看网站| 能在线免费看毛片的网站| videos熟女内射| 人体艺术视频欧美日本| 免费播放大片免费观看视频在线观看| 亚洲国产精品一区三区| 午夜激情福利司机影院| 欧美一级a爱片免费观看看| 亚洲第一区二区三区不卡| 精品一区二区免费观看| 人人妻人人添人人爽欧美一区卜 | 久久久久久久久久久丰满| 最近中文字幕2019免费版| 亚洲美女视频黄频| 18禁在线无遮挡免费观看视频| 精品人妻一区二区三区麻豆| 国产男女超爽视频在线观看| 国内揄拍国产精品人妻在线| 亚洲精品久久午夜乱码| 成人二区视频| 大又大粗又爽又黄少妇毛片口| 欧美日韩综合久久久久久| av免费在线看不卡| 中文字幕精品免费在线观看视频 | 51国产日韩欧美| 亚洲丝袜综合中文字幕| 国产乱人视频| 国产精品三级大全| 成人二区视频| 黄色日韩在线| 亚洲精品色激情综合| 久久久久网色| 高清视频免费观看一区二区| 日本色播在线视频| 美女内射精品一级片tv| 老女人水多毛片| 中文字幕av成人在线电影| 在线观看美女被高潮喷水网站| 久久久久久久久久久丰满| 在线看a的网站| 在线观看国产h片| 免费看日本二区| 一级毛片黄色毛片免费观看视频| 日本猛色少妇xxxxx猛交久久| 建设人人有责人人尽责人人享有的 | 一本一本综合久久| 日韩在线高清观看一区二区三区| 日本色播在线视频| 高清欧美精品videossex| 日本欧美视频一区| 亚洲真实伦在线观看| 亚洲,一卡二卡三卡| 亚洲人成网站高清观看| 亚洲真实伦在线观看| 日本vs欧美在线观看视频 | 亚洲欧美一区二区三区国产| 成人午夜精彩视频在线观看| 亚洲成人中文字幕在线播放| av国产久精品久网站免费入址| 在线观看免费视频网站a站| 欧美日韩视频精品一区| 午夜老司机福利剧场| 国产在线男女| 最新中文字幕久久久久| 热99国产精品久久久久久7| 成人午夜精彩视频在线观看| 国产一区二区三区av在线| 视频中文字幕在线观看| 国产 精品1| 中国美白少妇内射xxxbb| 亚洲美女视频黄频| 青春草国产在线视频| 黄色日韩在线| 一级二级三级毛片免费看| 晚上一个人看的免费电影| 久久国产精品大桥未久av | 自拍欧美九色日韩亚洲蝌蚪91 | 一本久久精品| 亚洲丝袜综合中文字幕| 深夜a级毛片| 99久久综合免费| 久久国产精品男人的天堂亚洲 | 伦理电影大哥的女人| 精品久久久久久电影网| 亚洲国产精品成人久久小说| 在线看a的网站| 欧美日韩视频精品一区| 大香蕉97超碰在线| 十八禁网站网址无遮挡 | 日本与韩国留学比较| 毛片女人毛片| 日韩,欧美,国产一区二区三区| 18禁裸乳无遮挡动漫免费视频| 人妻系列 视频| 最近中文字幕2019免费版| 欧美另类一区| 99久久人妻综合| 亚洲av欧美aⅴ国产| 日日摸夜夜添夜夜爱| 男人添女人高潮全过程视频| 精品久久久久久久久亚洲| 丝袜脚勾引网站| 亚洲经典国产精华液单| 激情五月婷婷亚洲| 尤物成人国产欧美一区二区三区| 男女边摸边吃奶| 97热精品久久久久久| 80岁老熟妇乱子伦牲交| 亚洲熟女精品中文字幕| 久久国内精品自在自线图片| a级一级毛片免费在线观看| 日本黄大片高清| 99久久中文字幕三级久久日本| 日本爱情动作片www.在线观看| 久久人人爽人人片av| 日韩电影二区| 永久网站在线| 18禁裸乳无遮挡动漫免费视频| 亚洲精品久久午夜乱码| 建设人人有责人人尽责人人享有的 | 国产精品一区二区在线观看99| 日韩人妻高清精品专区| 亚洲一级一片aⅴ在线观看| 亚洲aⅴ乱码一区二区在线播放| 国产片特级美女逼逼视频| 成年av动漫网址| 91精品一卡2卡3卡4卡| 国产乱人视频| 精品久久久久久久久av| 在线 av 中文字幕| 大又大粗又爽又黄少妇毛片口| 国产男女超爽视频在线观看| 亚洲欧美日韩无卡精品| 夫妻午夜视频| 国产亚洲最大av| 欧美日韩视频高清一区二区三区二| 久久精品久久精品一区二区三区| 欧美一区二区亚洲| 中文精品一卡2卡3卡4更新| 男女国产视频网站| 国产伦精品一区二区三区视频9| 国产老妇伦熟女老妇高清| 国产一区二区三区综合在线观看 | av国产免费在线观看| 国产精品久久久久成人av| 国产成人一区二区在线| 亚洲精品色激情综合| 中文字幕制服av| 日韩人妻高清精品专区| 久久午夜福利片| 91狼人影院| 欧美少妇被猛烈插入视频| 久久精品国产亚洲av天美| 日韩av在线免费看完整版不卡| 国语对白做爰xxxⅹ性视频网站| 高清视频免费观看一区二区| 国产亚洲91精品色在线| 内地一区二区视频在线| 国产一区二区三区av在线| 97精品久久久久久久久久精品| 波野结衣二区三区在线| 伦理电影免费视频| 国产伦精品一区二区三区四那| 日韩中字成人| 国语对白做爰xxxⅹ性视频网站| 久久毛片免费看一区二区三区| 免费不卡的大黄色大毛片视频在线观看| 嘟嘟电影网在线观看| 纵有疾风起免费观看全集完整版| 91精品国产国语对白视频| 蜜桃在线观看..| 亚洲国产欧美在线一区| 欧美成人精品欧美一级黄| 国内少妇人妻偷人精品xxx网站| 六月丁香七月| 色5月婷婷丁香| 黑人猛操日本美女一级片| 又黄又爽又刺激的免费视频.| 日本色播在线视频| 狂野欧美激情性bbbbbb| 国产熟女欧美一区二区| 久久久久久久精品精品| 午夜视频国产福利| 亚洲精品乱久久久久久| 精品亚洲乱码少妇综合久久| 五月伊人婷婷丁香| 99热这里只有是精品在线观看| 你懂的网址亚洲精品在线观看| 在线观看免费高清a一片| 久久久久国产网址| 亚洲欧美日韩卡通动漫| 国产伦理片在线播放av一区| 亚洲欧洲日产国产| 99热这里只有是精品50| 亚洲av成人精品一二三区| 性色av一级| 2022亚洲国产成人精品| 51国产日韩欧美| 大话2 男鬼变身卡| 亚洲四区av| av.在线天堂| 能在线免费看毛片的网站| 大香蕉久久网| 精品少妇久久久久久888优播| 十八禁网站网址无遮挡 | 久久久久精品久久久久真实原创| 人妻系列 视频| 亚洲av中文av极速乱| 国产老妇伦熟女老妇高清| 亚洲av日韩在线播放| 好男人视频免费观看在线| 九草在线视频观看| 国产黄片视频在线免费观看| 人妻一区二区av| 国产综合精华液| 国产伦理片在线播放av一区| 日韩强制内射视频| 国产成人一区二区在线| 人人妻人人添人人爽欧美一区卜 | 欧美日韩综合久久久久久| 久久 成人 亚洲| 国产在线视频一区二区| 啦啦啦视频在线资源免费观看| 亚洲最大成人中文|