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

    粗糙度對渦輪葉片流動轉(zhuǎn)捩及傳熱特性的影響

    2016-07-31 19:31:23李虹楊
    關(guān)鍵詞:模型

    李虹楊,鄭 赟,*

    粗糙度對渦輪葉片流動轉(zhuǎn)捩及傳熱特性的影響

    李虹楊1,2,鄭 赟1,2,*

    (1.北京航空航天大學(xué)能源與動力工程學(xué)院,北京100083; 2.先進(jìn)航空發(fā)動機(jī)協(xié)同創(chuàng)新中心,北京100083)

    為研究表面粗糙度對渦輪葉片流動轉(zhuǎn)捩以及傳熱特性的影響,在自行開發(fā)的CFD程序平臺上提出了對γ-Reθ轉(zhuǎn)捩模型的粗糙度修正方法,并參考平板繞流和渦輪葉柵的實(shí)驗(yàn)數(shù)據(jù)對該方法迸行驗(yàn)證??紤]粗糙度效應(yīng)的影響,對MarkⅡ渦輪導(dǎo)葉5411工況迸行數(shù)值模擬,得到如下結(jié)論:表面粗糙度對層流邊界層換熱系數(shù)影響不大,而對湍流邊界層則有較大影響,迸而顯著改變壁面溫度分布;與光滑壁面相比,5μm的等效沙粒粗糙度使吸力面湍流區(qū)域壁面溫度升高約5.7K,100μm粗糙度使壁面溫度升高28.4 K,增幅達(dá)5%左右;當(dāng)壁面粗糙度較低時,激波干涉對吸力面邊界層的轉(zhuǎn)捩起主導(dǎo)作用,而當(dāng)粗糙度大于某臨界值時,其作用會使轉(zhuǎn)捩位置突然變化,本算例中該臨界值近似為150μm。

    粗糙度;轉(zhuǎn)捩;渦輪;間歇因子;邊界層

    現(xiàn)代先進(jìn)航空發(fā)動機(jī)為追求卓越性能,渦輪前燃?xì)鉁囟炔粩嗵岣?,甚至已?jīng)遠(yuǎn)遠(yuǎn)超過了金屬材料的熔點(diǎn),這對渦輪葉片的冷卻設(shè)計(jì)提出了更高的要求[1]。對渦輪內(nèi)部流動狀態(tài)及傳熱特性的準(zhǔn)確預(yù)測是高效冷卻設(shè)計(jì)的前提,因此對相關(guān)流動機(jī)理、影響因素以及變化規(guī)律的研究十分重要。渦輪葉片的傳熱特性通常受到壓力分布、尾跡干涉、邊界層轉(zhuǎn)捩、流動分離、壁面曲率、二次流等因素的影響[2],而表面粗糙度也是重要影響因素之一。

    發(fā)動機(jī)的渦輪部件工作環(huán)境十分惡劣,積垢、氧化、侵蝕等多重效應(yīng)的作用[3]可能使渦輪葉片表面的粗糙度達(dá)到“過渡粗糙”或者“完全粗糙”數(shù)量級[4],其影響主要體現(xiàn)在2個方面[5]:一方面,能改變湍流邊界層的阻力系數(shù)和傳熱系數(shù),直接影響氣動及傳熱特性;另一方面,壁面粗糙度可能對邊界層的流動狀態(tài),如層流到湍流的轉(zhuǎn)捩、流動分離等產(chǎn)生影響,間接影響葉片的傳熱特性。

    Taylor[6]利用表面光度儀對航空發(fā)動機(jī)渦輪葉片進(jìn)行測量,認(rèn)為典型的粗糙高度為1~12μm。Bons等[7]對近100個使用期內(nèi)的不同尺寸的渦輪部件進(jìn)行研究,指出在葉片中截面處的平均粗糙高度可達(dá)到37μm。Barlow[8]和Hosni[9]等對平板的傳熱性能進(jìn)行實(shí)驗(yàn)研究,結(jié)果表明特定的表面粗糙度使平板熱傳導(dǎo)系數(shù)增加了120%。Abuaf等[10]的實(shí)驗(yàn)發(fā)現(xiàn)隨著表面粗糙度的增加,轉(zhuǎn)捩發(fā)生的更早并且湍流邊界層具有更高的熱傳導(dǎo)系數(shù)。Bunker[11]研究了最大粗糙高度為27.8μm的跨聲速葉柵流動,結(jié)果表明與自由流湍流度相比,表面粗糙度對轉(zhuǎn)捩的影響占主導(dǎo)作用。Boyle等[12]和Blair[13]也對粗糙壁面的渦輪葉柵進(jìn)行了研究。

    表面粗糙度對邊界層轉(zhuǎn)捩有重要影響,模擬粗糙壁面的流動首先要準(zhǔn)確預(yù)測轉(zhuǎn)捩,這對數(shù)值模擬方法提出了新的挑戰(zhàn)。與實(shí)驗(yàn)研究相比,已公開發(fā)表的數(shù)值模擬方面的研究成果相對較少。Stripf等[14-15]做出開創(chuàng)性工作,提出基于粗糙單元高度的轉(zhuǎn)捩起始位置經(jīng)驗(yàn)公式,并結(jié)合離散單元粗糙度(discrete element roughness)模型[16-17]與雙層k-ε(two layer k-ε)模型(k為湍動能,ε為湍動能耗散率),模擬了粗糙壁面的邊界層轉(zhuǎn)捩。Boyle和Stripf[18]針對不同的粗糙單元的幾何形狀提出一種新的表面等效沙粒粗糙度計(jì)算方法,并利用渦輪葉柵的實(shí)驗(yàn)結(jié)果進(jìn)行了驗(yàn)證。Lorenz等[19]基于Stripf[14-15]的工作改進(jìn)了粗糙度準(zhǔn)則公式,使結(jié)果與實(shí)驗(yàn)值符合得更好。另一方面,Hellsten和Seppo[20]、Aupoix[21]以及Bellucci等[22]對基于k-ωShear Stress Transport(SST)兩方程模型的壁面粗糙度修正法進(jìn)行研究,但這些方法不能對邊界層轉(zhuǎn)捩進(jìn)行有效預(yù)測。

    本文在Hellsten和Seppo[20]及Stripf等[14]工作的基礎(chǔ)上,提出了針對γ-Reθ轉(zhuǎn)捩模型[23]的表面粗糙度修正方法(γ為間歇因子,Reθ為邊界層的轉(zhuǎn)捩動量厚度雷諾數(shù)):一方面,采用新的適用于粗糙表面的邊界條件,以模擬粗糙表面對湍流邊界層的影響;另一方面,改進(jìn)模型中的經(jīng)驗(yàn)關(guān)聯(lián)函數(shù),使模型適用于預(yù)測粗糙表面流動的轉(zhuǎn)捩預(yù)測。利用粗糙平板的風(fēng)洞試驗(yàn)數(shù)據(jù)和某低壓渦輪葉柵實(shí)驗(yàn)數(shù)據(jù)對模型進(jìn)行驗(yàn)證,并用該模型對粗糙壁面的MarkⅡ渦輪葉柵的流動進(jìn)行數(shù)值模擬,研究粗糙度對邊界層轉(zhuǎn)捩及壁面?zhèn)鳠崽匦缘挠绊憽?/p>

    1 數(shù)值方法

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

    γ-Reθ轉(zhuǎn)捩模型由Langtry和Menter[23]提出,該模型有兩個輸運(yùn)方程,分別是轉(zhuǎn)捩起始位置的邊界層動量厚度雷諾數(shù)(?Reθt)的輸運(yùn)方程和間歇因子(γ)的輸運(yùn)方程。其中?Reθt的輸運(yùn)方程為

    式中:ρ為密度;t為時間;Uj和xj分別為張量表示的速度分量和坐標(biāo)分量;Pθt為源項(xiàng),詳細(xì)表達(dá)式參見文獻(xiàn)[23];μ為層流渦黏性系數(shù);σθt=2.0為常數(shù)。

    Pθt作用是使邊界層外的t等于Reθt,而邊界層內(nèi)的從邊界層外擴(kuò)散得到。γ-Reθ轉(zhuǎn)捩模型的另一個輸運(yùn)方程,即間歇因子γ的輸運(yùn)方程表達(dá)式為

    式中:Pγ為控制間歇因子增長的源項(xiàng);σγ=1.0為常數(shù);Eγ為破壞項(xiàng),詳細(xì)的表達(dá)式參見文獻(xiàn)[23]。

    通過式(1)和式(2)兩個輸運(yùn)方程,可以得到流場中間歇因子的分布,還需要與湍流模型相結(jié)合才能模擬轉(zhuǎn)捩過程。Langtry和Menter推薦k-ω SST模型,即用間歇因子來修正湍動能k的輸運(yùn)方程的生成項(xiàng)、破壞項(xiàng)和混合函數(shù)[23-24]。

    1.2 表面粗糙度修正方法

    一方面,為了考慮表面粗糙度對湍流邊界層的影響,Hellsten和Seppo[20]提出了對k-ωSST模型的粗糙度修正方法。通過對壁面處湍動能k和比耗散率ω的修正來模擬粗糙度的影響,具體方法如下:

    式中:ωw為壁面位置的湍動能比耗散率;uτ=為摩擦速度,τw為壁面剪切應(yīng)力;ν為運(yùn)動黏性系數(shù);SR為無量綱系數(shù),其定義為

    式中:ks為表面等效沙粒粗糙度,以下簡稱為表面粗糙度。

    為準(zhǔn)確模擬邊界層內(nèi)的流動,還需要在k-ω SST模型的渦黏性表達(dá)式中增加一個摻混函數(shù)F3,具體表達(dá)式參見文獻(xiàn)[20]。

    另一方面,為考慮表面粗糙度對邊界層轉(zhuǎn)捩位置的影響,Strip f等[14,25]重新定義了一個粗糙表面的轉(zhuǎn)捩動量厚度雷諾數(shù)Reθt_Rough,該變量與等效沙粒粗糙度ks和邊界層的位移厚度 δ*有關(guān),其表達(dá)式為

    式中:函數(shù)fΛ用于描述粗糙單元幾何結(jié)構(gòu),即形狀、排列規(guī)律等的影響,而本文未考慮這些影響因素,令為當(dāng)?shù)刈杂闪魍牧鞯暮瘮?shù),其定義為

    式中:σTu為當(dāng)?shù)赝牧鞫取?/p>

    1.3 考慮粗糙度影響的轉(zhuǎn)捩模型

    首先,模型要考慮表面粗糙度對湍流邊界層的影響,表面粗糙度會對轉(zhuǎn)捩位置之后的湍流邊界層產(chǎn)生直接的影響,如改變邊界層內(nèi)部的湍動能和渦粘性系數(shù)。γ-Reθ轉(zhuǎn)捩模型對湍流邊界層的計(jì)算方法與常規(guī)湍流模型是近似一致的,因此Hellsten和Seppo[20]的邊界條件修正方法同樣適用。本文將該方法引入到γ-Reθ轉(zhuǎn)捩模型中,模擬的邊界層轉(zhuǎn)捩點(diǎn)后湍流邊界層區(qū)域內(nèi)的阻力系數(shù)與M ills-Hang的經(jīng)驗(yàn)公式[20]符合得很好。

    粗糙表面流動轉(zhuǎn)捩預(yù)測的另一個重要問題就是模擬表面粗糙度對轉(zhuǎn)捩位置的影響。原始 γ-Reθ模型中的經(jīng)驗(yàn)關(guān)聯(lián)公式均是基于光滑平板低速繞流實(shí)驗(yàn)數(shù)據(jù)擬合得到的,而其中最為重要的是起始位置邊界層轉(zhuǎn)捩動量厚度雷諾數(shù)與當(dāng)?shù)赝牧鞫?、速度梯度、壓力梯度等的?jīng)驗(yàn)關(guān)聯(lián)公式??紤]表面粗糙度的影響,對該公式進(jìn)行改寫,使模型具備預(yù)測粗糙表面流動轉(zhuǎn)捩的預(yù)測能力是可行的。

    Stripf等[14]提出的經(jīng)驗(yàn)公式已經(jīng)給出了一個基于實(shí)驗(yàn)數(shù)據(jù)擬合得到的轉(zhuǎn)捩動量厚度雷諾數(shù)與表面等效沙粒粗糙度的關(guān)系式(式(6)和式(7)),但該關(guān)系式與γ-Reθ轉(zhuǎn)捩模型中的經(jīng)驗(yàn)關(guān)聯(lián)公式不同,式(6)中的邊界層位移厚度和轉(zhuǎn)捩動量厚度是全局變量,進(jìn)行數(shù)值模擬時需要進(jìn)行積分求解。為克服這個問題,本文引入流場當(dāng)?shù)匚灰坪穸?,重新建立形如式?)和式(7)的表達(dá)式,并與γ-Reθ轉(zhuǎn)捩模型中的輸運(yùn)方程相結(jié)合,進(jìn)而得到粗糙表面轉(zhuǎn)捩動量厚度雷數(shù)Reθt_Rough的分布,值得注意的是,這里的δ*與Reθt均為流中的當(dāng)?shù)刈兞?,與Strip f等[14]的公式中的變量含義不同。下面介紹當(dāng)?shù)刈兞课灰坪穸圈?的計(jì)算方法,γ-Reθ模型的公式中轉(zhuǎn)捩動量厚度需要進(jìn)行迭代求解,將Reθt的定義式(8)帶入經(jīng)驗(yàn)關(guān)聯(lián)式(9),可見等號兩端均含有轉(zhuǎn)捩動量厚度θt,需要通過簡單的迭代方法進(jìn)行求解。

    式中:U為流場當(dāng)?shù)厮俣?;G(σTu)和 F(λθ)為γ-Reθ模型中的公式,具體表達(dá)式參見文獻(xiàn)[23]。

    得到流場中轉(zhuǎn)捩動量厚度θt后,需要估算出位移厚度δ*,本文通過算例驗(yàn)證,采用了計(jì)算簡便的湍流邊界層速度剖面的 N次方定律[26]來求解。

    式中:θ為邊界層動量厚度,δ為邊界層厚度,n取為7。得到邊界層位移厚度 δ*之后,結(jié)合公式(6)即可得到流場中Reθt_Rough的分布。

    Langtry和Menter[23]的γ-Reθ轉(zhuǎn)捩模型以及上述修正方法在自行開發(fā)的CFD程序HGFS上實(shí)現(xiàn),程序的介紹及應(yīng)用情況可參見文獻(xiàn)[27-30]。

    2 模型的驗(yàn)證

    首先,參照Wang[31]和Pinson[32]等的風(fēng)洞試驗(yàn)數(shù)據(jù),對粗糙平板的轉(zhuǎn)捩流動進(jìn)行驗(yàn)證。分別選擇了光滑壁面以及表面等效沙粒粗糙度ks分別為150μm和400μm,來流湍流度σTu∞在0.5%~5.2%的工況進(jìn)行數(shù)值模擬,如表1所示。

    湍流邊界層熱傳導(dǎo)系數(shù)一般遠(yuǎn)高于層流邊界層,因此傳熱特性曲線可以明顯地反映出轉(zhuǎn)捩位置,計(jì)算得到的Stanton數(shù)與實(shí)驗(yàn)值的對比如圖1所示,圖中橫坐標(biāo)Rex表示當(dāng)?shù)匚恢玫睦字Z數(shù),曲線表示數(shù)值模擬的結(jié)果,空心三角符號表示光滑表面的實(shí)驗(yàn)結(jié)果,實(shí)心符號(圓點(diǎn)、正方形、菱形)表示相對應(yīng)的粗糙表面的實(shí)驗(yàn)結(jié)果??梢钥闯鰧τ诓煌砻娲植诙群蛠砹魍牧鞫鹊那闆r,本文模型所預(yù)測的轉(zhuǎn)捩位置均與實(shí)驗(yàn)值符合得很好,但在較高表面粗糙度情況下,如ks=400μm,預(yù)測的轉(zhuǎn)捩區(qū)長度比實(shí)驗(yàn)結(jié)果稍短。

    其次,為驗(yàn)證本文的方法與模型在復(fù)雜流動情況中的適用性,選擇某高壓渦輪導(dǎo)葉(HPTV)算例進(jìn)行驗(yàn)證,對應(yīng)的實(shí)驗(yàn)由德國Karlsruhe大學(xué)[33]完成。弦長c=93.95mm,來流雷諾數(shù)Re∞=1.4×105,來流湍流度σTu∞=3.5%。不同表面粗糙度參數(shù)如表2所示。

    表1 粗糙平板算例的計(jì)算參數(shù)Table 1 Calculation param eters for rough flat p late cases

    圖1 粗糙平板傳熱系數(shù)與實(shí)驗(yàn)值的對比Fig.1 Comparison between heat transfer coefficient on rough flat plate and experimental data

    表2 HPTV算例的表面粗糙度參數(shù)Tab le 2 Sur face roughness param eters for HPTV case

    計(jì)算采用如圖2所示二維非結(jié)構(gòu)網(wǎng)格,單元總數(shù)為24 000,壁面附近進(jìn)行加密處理,保證最大y+<0.8,從局部視圖來看,前緣和尾緣均有較高的網(wǎng)格分辨率,經(jīng)過網(wǎng)格無關(guān)性驗(yàn)證,此網(wǎng)格滿足計(jì)算要求。

    圖3是數(shù)值模擬的HPTV葉片表面?zhèn)鳠崽匦郧€與實(shí)驗(yàn)值的對比,s為吸弧長(圖中所示為葉片吸力面),Nuc為努賽爾數(shù),其定義為

    圖2 HPTV算例的網(wǎng)格Fig.2 Computational grid for HPTV case

    圖3 粗糙表面渦輪葉片的傳熱系數(shù)Fig.3 Heat transfer coefficient of turbine vane with rough surface

    式中:λ為恒定的熱傳導(dǎo)系數(shù)[33],本算例中 λ=0.03W/(m·K);Tt_inlet為進(jìn)口總溫;Tw為葉片表面溫度。圖3中帶標(biāo)記的實(shí)線表示不同表面粗糙度條件下的計(jì)算結(jié)果,而相應(yīng)的實(shí)心符號則表示該條件下的實(shí)驗(yàn)結(jié)果。從圖3中可以看出本文模型預(yù)測的Nuc幅值與變化規(guī)律與實(shí)驗(yàn)符合得很好;對光滑壁面的轉(zhuǎn)捩預(yù)測比較準(zhǔn)確,對于 ks=72μm的工況預(yù)測的轉(zhuǎn)捩位置與轉(zhuǎn)捩區(qū)長度也與實(shí)驗(yàn)結(jié)果接近;ks=129,238μm的粗糙度下,預(yù)測的轉(zhuǎn)捩位置與實(shí)驗(yàn)值接近,但轉(zhuǎn)捩區(qū)長度較短,與實(shí)驗(yàn)值有一定差別??傮w來看,本文模型對渦輪葉柵中表面粗糙度誘導(dǎo)的邊界層轉(zhuǎn)捩流動有較為理想的預(yù)測精度。

    3 計(jì)算結(jié)果與分析

    利用本文模型對內(nèi)冷渦輪導(dǎo)葉MarkⅡ的5411工況[34]進(jìn)行流/熱耦合數(shù)值模擬,研究粗糙壁面對其流動及傳熱性能的影響。

    3.1 光滑壁面條件下的驗(yàn)證

    HGFS程序中流/熱耦合計(jì)算處理方法為流體域、固體域分別求解,在邊界進(jìn)行信息耦合,MarkⅡ渦輪導(dǎo)葉是驗(yàn)證流/熱耦合計(jì)算程序的經(jīng)典算例,其5411工況主要進(jìn)出口邊界條件如表3所示。

    計(jì)算采用圖4所示的二維非結(jié)構(gòu)網(wǎng)格,其中流體域單元總數(shù)為87 000,計(jì)算得到的y+在0.4~0.9范圍內(nèi),固體域單元總數(shù)為12000。

    表3 M arkⅡ?qū)~5411工況的邊界條件Tab le 3 No.5411 boundary condition of M arkⅡturbine vane

    圖4 MarkⅡ?qū)~計(jì)算網(wǎng)格Fig.4 Computational grid for MarkⅡturbine vane

    分別利用S-A模型、k-ωSST模型以及γ-Reθ轉(zhuǎn)捩模型對Mark II葉片(光滑壁面)的流動進(jìn)行數(shù)值模擬,計(jì)算得到的壓力、溫度和熱流量的分布如圖5(a)~圖5(c)所示。L為軸向弦長,橫坐標(biāo)0代表前緣點(diǎn),左側(cè)(即x/L<0的范圍)表示吸力面,右側(cè)表示壓力面,圖5~圖9均采用相同的表示方法。不同模型計(jì)算得到的壁面壓力分布近似一致,且均與實(shí)驗(yàn)值符合得很好。而溫度分布則有較大的差別,主要原因在于S-A與k-ωSST等全湍流模型的結(jié)果中轉(zhuǎn)捩是非物理的,且極容易發(fā)生,所預(yù)測的轉(zhuǎn)捩位置通常遠(yuǎn)遠(yuǎn)早于真實(shí)情況,從圖5(b)和圖5(c)的結(jié)果來看這2個模型模擬的邊界層在葉片吸力面幾乎從前緣點(diǎn)就已經(jīng)轉(zhuǎn)捩,導(dǎo)致0~50%軸向弦長(以下簡稱為弦長)壁面溫度偏高。γ-Reθ模型預(yù)測的轉(zhuǎn)捩位置與實(shí)驗(yàn)結(jié)果非常接近,轉(zhuǎn)捩點(diǎn)約在50%弦長位置,轉(zhuǎn)捩點(diǎn)之后的位置有一定的“過調(diào)”現(xiàn)象,但壁面的溫度和熱流量的分布總體與實(shí)驗(yàn)值符合得很好。

    圖5 光滑壁面條件下的驗(yàn)證Fig.5 Verification in smooth surface condition

    圖6 利用γ-Reθ模型數(shù)值模擬得到的流場Fig.6 Numerical simulated flow field byγ-Reθmodel

    利用γ-Reθ轉(zhuǎn)捩模型數(shù)值模擬得到的流場如圖6所示,其中圖6(a)~圖6(c)分別為溫度(T)云圖、馬赫數(shù)(Ma)云圖以及間歇因子(γ)云圖。由于葉片內(nèi)部馬赫數(shù)為0,因此圖6(b)中并未示出(以白色背景代替),與之類似,圖6(c)中只顯示了流場中邊界層內(nèi)間歇因子的分布,在葉片內(nèi)部間歇因子為0,而流場廣闊區(qū)域內(nèi)該值的分布接近于1。從圖6(a)和圖6(b)中可以清晰看出葉片吸力面的激波,結(jié)合圖6(c)間歇因子的分布可以看出,激波位置之前,吸力面邊界層為層流流動,激波之后轉(zhuǎn)變?yōu)橥牧?;而壓力面邊界層始終為層流,只有在尾緣點(diǎn)(方形尾緣)才轉(zhuǎn)變?yōu)橥牧鳌?/p>

    3.2 粗糙壁面對流動及傳熱的影響

    保證其他邊界條件與來流參數(shù)均與3.1節(jié)MarkⅡ5411工況相同,只增加壁面粗糙度的影響,數(shù)值模擬的溫度分布曲線如圖7所示,圖中列出了光滑壁面以及葉片表面粗糙度ks=5~100μm的計(jì)算結(jié)果。其中:圖7(a)為全局視圖,圖7(b)和圖7(c)分別為吸力面、壓力面虛線方框內(nèi)的局部放大視圖。圖8為粗糙度ks=100~175μm的壁面溫度分布,圖9為邊界層內(nèi)間歇因子分布。

    從圖7(a)可以看出,在5~100μm范圍內(nèi),粗糙度不同的壁面計(jì)算得到的溫度分布呈現(xiàn)較為一致的規(guī)律;壓力面0~20%弦長以及吸力面0~50%弦長的溫度分布幾乎不受粗糙度的影響,壓力面20%弦長到尾緣的范圍內(nèi)受到的影響較小,而吸力面50%弦長以后的位置受表面粗糙度影響較大。

    從圖7(b)中的溫度分布曲線來看,與光滑壁面相比,粗糙度為5μm的壁面使該區(qū)域壁面相對溫度,即Tw/(811K)升高了約0.007,絕對溫度升高約5.7 K,而粗糙度為100μm的壁面使相對溫度升高約0.035,絕對溫度升高28.4 K,增幅約為5%。在圖7(c)中,100μm的壁面粗糙度使壓力面相對溫度升高約0.023,絕對溫度升高18.7 K,增幅小于吸力面湍流區(qū)。從上述結(jié)果中可以得到如下結(jié)論:

    1)隨著表面粗糙度的增加,湍流邊界層換熱效果有較大幅度增強(qiáng),導(dǎo)致壁面溫度升高,而完全層流邊界層區(qū)域內(nèi)的壁面溫度保持不變,如圖中x/L=-0.2~0.5的區(qū)間內(nèi)所示。

    2)湍流邊界層內(nèi),隨著表面粗糙度的增加,溫度的變化規(guī)律并不是線性的,粗糙度越大其影響效果越弱,即隨著粗糙度的逐漸增大,溫度增加的幅度在逐漸減小,熱流量也體現(xiàn)相同的規(guī)律。

    3)從邊界層內(nèi)間歇因子分布(見圖9)來看,ks=5~100μm范圍內(nèi)葉片壓力面并未發(fā)生轉(zhuǎn)捩,邊界層為層流,但表面粗糙度對壓力面20%弦長到尾緣的區(qū)域也有一定影響,說明該區(qū)域的邊界層呈現(xiàn)出一定的湍流脈動特性。

    圖7 不同壁面粗糙度條件下的溫度分布(ks=5~100μm)Fig.7 Distribution of wall temperature in different surface-roughness conditions(ks=5-100μm)

    從圖8可以看出ks為100μm和150μm時溫度曲線仍與圖7有相同的變化規(guī)律,即圖8(b)中湍流邊界層內(nèi)壁面溫度進(jìn)一步升高;當(dāng)ks增加到160μm時,溫度分布曲線出現(xiàn)較大變化,吸力面15%~50%弦長范圍內(nèi)溫度大幅增加,由原來的低谷區(qū)變成新的峰值區(qū)。結(jié)合圖9中間歇因子的分布可以看出,ks增加到160μm時吸力轉(zhuǎn)捩位置由50%弦長處前移到15%弦長附近,轉(zhuǎn)捩位置之后的邊界層由層流變?yōu)橥牧?,換熱能力的增強(qiáng)使壁面溫度顯著升高。

    圖8 不同壁面粗糙度條件下的溫度分布(ks=100~175μm)Fig.8 Distribution of wall temperature in different surface-roughness conditions(ks=100-175μm)

    圖8 (c)中,當(dāng)ks增加到150μm以上時,壓力面溫度分布曲線在80%~90%弦長位置增加幅度較大,從圖9可以看出,較高的壁面粗糙度使壓力面邊界層也發(fā)生了轉(zhuǎn)捩。

    邊界層內(nèi)間歇因子的分布如圖9所示,在壓力面ks=5~50μm和ks=100μm的分布曲線幾乎完全重合,圖例中ks=5~50μm的曲線用不同類型的虛線表示,而ks≥100μm的曲線用帶有符號的實(shí)線表示。間歇因子的分布能直觀地反映出邊界層的轉(zhuǎn)捩情況:

    1)對于吸力面,ks=5~150μm范圍內(nèi),邊界層轉(zhuǎn)捩位置均在50%弦長附近;當(dāng)ks增加到160μm時,轉(zhuǎn)捩位置前移到15%弦長附近,而當(dāng)ks繼續(xù)增加到175μm時,轉(zhuǎn)捩位置略有提前。

    2)對于壓力面,ks=5~100μm范圍內(nèi),轉(zhuǎn)捩發(fā)生在尾緣點(diǎn),即整個吸力面均為層流流動;當(dāng)ks增加到150μm時,轉(zhuǎn)捩位置前移到90%弦長,且當(dāng)ks繼續(xù)增到160μm和175μm時,轉(zhuǎn)捩位逐漸前移到85%和80%弦長附近。即壓力面的轉(zhuǎn)捩首先發(fā)生在尾緣,且隨著粗糙度的增加轉(zhuǎn)捩位置逐漸前移。

    3)對比ks為150、160和175μm的結(jié)果可以看出吸力面轉(zhuǎn)捩位置的變化并不是線性的,與第2節(jié)驗(yàn)證算例中HPTV渦輪葉片的結(jié)果不同。這是因?yàn)镸arkⅡ吸力面的轉(zhuǎn)捩是由激波邊界層干涉的誘導(dǎo)而發(fā)生的,當(dāng)壁面粗糙度較小時,激波對轉(zhuǎn)捩的作用效果為主導(dǎo),只有當(dāng)粗糙度大于某臨界值時,其較強(qiáng)的影響效果才會使轉(zhuǎn)捩點(diǎn)發(fā)生變化,本算例中該臨界值近似為150μm。

    圖9 邊界層內(nèi)間歇因子分布Fig.9 Distribution of intermittency factor in boundary layer

    4 結(jié) 論

    在內(nèi)部CFD程序HGFS中對γ-Reθ轉(zhuǎn)捩模型進(jìn)行改進(jìn),使該模型能夠預(yù)測粗糙壁面的轉(zhuǎn)捩流動。利用本文模型對MarkⅡ渦輪導(dǎo)葉的流動進(jìn)行數(shù)值模擬,得到如下結(jié)論:

    1)隨著表面粗糙度的增加,葉片層流邊界層換熱系數(shù)變化不大,而湍流邊界層換熱系數(shù)有較大幅度增加,進(jìn)而導(dǎo)致壁面溫度升高。且隨著表面粗糙度的增加,溫度和熱流量的變化規(guī)律并不是線性的,粗糙度越大其影響效果越弱。

    2)在葉片吸力面邊界層為湍流的區(qū)域,與光滑壁面相比,等效沙粒粗糙度為5μm的壁面使該區(qū)域溫度升高約5.7 K,而粗糙度為100μm的壁面使溫度升高28.4 K,增幅達(dá)5%左右。

    3)粗糙度對葉片吸力面轉(zhuǎn)捩位置的影響不是線性的,當(dāng)粗糙度大于某一臨界值時轉(zhuǎn)捩位置發(fā)生突變。吸力面邊界層轉(zhuǎn)捩是由激波邊界層干涉誘導(dǎo)而產(chǎn)生的,當(dāng)粗糙度較低時(ks<100μm),激波的作用效果為主導(dǎo),只有當(dāng)粗糙度大于某臨界值時才會對轉(zhuǎn)捩產(chǎn)生影響,本算例中該臨界值近似為150μm。

    4)光滑壁面以及較低粗糙度(ks<100μm)范圍內(nèi),壓力面邊界層均為層流,當(dāng)粗糙度增大到150μm時,壓力面的轉(zhuǎn)捩首先發(fā)生在尾緣,且隨著粗糙度的增加轉(zhuǎn)捩位置逐漸前移。

    本文對γ-Reθ模型的修正方法在一定粗糙度范圍內(nèi)得到了規(guī)律比較一致的結(jié)果,且預(yù)測的溫度、熱流量、轉(zhuǎn)捩位置等參數(shù)的變化趨勢與真實(shí)流動情況相符,但仍需要更為充分的實(shí)驗(yàn)數(shù)據(jù)對該方法進(jìn)行驗(yàn)證與完善,這是下一步需要進(jìn)行的工作。

    (References)

    [1]BUNKER R S.A review of shaped hole turbine film-cooling technology[J].Journal of Heat Transfer,2005,127(4):441-453.

    [2]YEH F C,STEPKA F S.Review and status of heat-transfer technology for internal passages of air-cooled turbine blades:NASA-TP-2232[R].Washington,D.C:NASA,1984.

    [3]李本威,李冬,沈偉,等.渦輪葉片粗糙度對其性能衰退的影響研究[J].航空計(jì)算技術(shù),2009,39(5):26-29.

    LIB W,LI D,SHEN W,et al.Research on turbine lamina roughness influence on its performance declination[J].Aeronautical Computing Technique,2009,39(5):26-29(in Chinese).

    [4]BOGARD D G,SCHMIDT D L,TABBITA M.Characterization and laboratory simulation of turbine airfoil surface roughness and associated heat transfer[J].Journal of Turbomachinery,1998,120(2):337-342.

    [5]BONS JP.A review of surface roughness effects in gas turbines[J].Journal of Turbomachinery,2010,132(2):021004.

    [6]TAYLOR R P.Surface roughness measurements on gas turbine blades[J].Journal of Turbomachinery,1990,112(2):175-180.

    [7]BONS JP,TAYLOR R P,MCCLAIN ST,et al.Themany faces of turbine surface roughness[J].Journal of Turbomachinery,2001,123(4):739-748.

    [8]BARLOW D N,KIM Y W,F(xiàn)LORSCHUETZ LW.Transient liquid crystal technique for convective heat transfer on rough surfaces[J].Journal of Turbomachinery,1997,119(1):14-22.

    [9]HOSNIM H,COLEMAN H W,TAYLOR R P.Rough-wall heat transfer in turbulent boundary layers[J].International Journal of Fluid Mechanics,1998,25(1-3):212-219.

    [10]ABUAF N N,BUNKER R S,LEE C P.Effects of surface roughness on heat transfer and aerodynamic performance of turbine airfoils[J].Journal of Turbomachinery,1998,120(3):522-529.

    [11]BUNKER R S.Separate and combined effects of surface roughness and turbulence intensity on vane heat transfer:97-GT-135[R].New York:ASME,1997.

    [12]BOYLE R J,SPUCKLER C M,LUCCIB L,et al.Infrared lowtemperature turbine vane rough surface heat transfer measurements[J].Journal of Turbomachinery,2001,123(1):168-177.

    [13]BLAIR M F.An experimental study of heat transfer in a largescale turbine rotor passage[J].Journal of Turbomachinery,1994,116(1):1-13.

    [14]STRIPF M,SCHULZ A,BAUER H J,et al.Extended models for transitional rough wall boundary layers with heat transfer-PartⅠ:Model formulations[J].Journal of Turbomachinery,2009,131(3):1263-1275.

    [15]STRIPF M,SCHULZ A,BAUER H J,et al.Extended models for transitional rough wall boundary layers with heat transfer-PartⅡ:Model validation and benchmarking[J].Proceedings of the ASME Turbo Expo,2008,131(3):1277-1289.

    [16]TAYLOR R P,COLEMAN H W,HODGE B K.A discrete element prediction approach for turbulent flow over rough surfaces[C]∥Viscous and Interacting Flow Field Effects.1984,1:1-11.

    [17]TAYLOR R P,COLEMAN H W,HODGE B K.Prediction of turbulent rough-wall skin friction using a discrete element approach[J].Journal of Fluids Engineering,1985,107(2):251-257.

    [18]BOYLE R J,STRIPF M.Simplified approach to predicting rough surface transition[J].Journal of Turbomachinery,2009,131(4):10-20.

    [19]LORENZM,SCHULZ A,BAUER H J.Predicting rough wall heat transfer and skin friction in transitional boundary layers-A new correlation for bypass transition onset[J].Journal of Turbomachinery,2013,135(4):10-21.

    [20]HELLSTEN A,SEPPO L.Extension of the k-ωSST turbulence model for flows over rough surfaces:AIAA-1997-3577[R].Reston:AIAA,1997.

    [21]AUPOIX B B.Roughness corrections for the k-ωshear stress transport model:Status and proposals[J].ASME,Journal of Fluids Engineering,2014,137(2):021202.

    [22]BELLUCCI J,RUBECHINIF,MARCOCINIM,et al.The influence of roughness on a high-pressure steam turbine stage:An experimental and numerical study[J].Journal of Engineering for Gas Turbines&Power,2015,137(1):012602.

    [23]LANGTRY R B,MENTER F R.Correlation-based transitionmodeling for unstructured parallelized computational fluid dynamics codes[J].AIAA Journal,2009,47(12):2894-2906.

    [24]MENTER F R,LANGTRY R B,V?LKER S.Transition modelling for general purpose CFD codes[J].Flow,Turbulence and Combustion,2006,77(1-4):277-303.

    [25]STRIPF M,SCHULZ A,BAUER H J.Modeling of rough wall boundary layer transition and heat transfer on turbine airfoils[J].Proceedings of the ASME Turbo Expo,2006,130(2):1139-1151.

    [26]陳懋章.粘性流體動力學(xué)基礎(chǔ)[M].北京:高等教育出版社,2002:312-314.

    CHEN M Z.Fundamentals of viscous fluid dynamics[M].Beijing:Higher Education Press,2002:312-314(in Chinese).

    [27]鄭赟.基于非結(jié)構(gòu)網(wǎng)格的氣動彈性數(shù)值方法研究[J].航空動力學(xué)報,2009,24(9):2069-2077.

    ZHENG Y.Computational aero-elasticity with an unstructured grid method[J].Journal of Aerospace Power,2009,24(9):2069-2077(in Chinese).

    [28]肖大啟,鄭赟,楊慧.軸向間距對轉(zhuǎn)子葉片氣動激勵的影響[J].航空動力學(xué)報,2012,27(10):2307-2313.

    XIAO D Q,ZHENG Y,YANG H.Effect of axial spacing on aerodynamic excitation of rotor blade[J].Journal of Aerospace Power,2012,27(10):2307-2313(in Chinese).

    [29]鄭赟,李虹楊,劉大響.γ-Reθ轉(zhuǎn)捩模型在高超聲速下的應(yīng)用及分析[J].推進(jìn)技術(shù),2014,35(3):296-304.

    ZHENG Y,LIH Y,LIU D X.Application and analysis ofγ-Reθtransition model in hypersonic flow[J].Journal of Propulsion Technology,2014,35(3):296-304(in Chinese).

    [30]鄭赟,李虹楊.基于新的經(jīng)驗(yàn)關(guān)聯(lián)公式的γ-Reθ轉(zhuǎn)捩模型在高超聲速流動中的應(yīng)用[J].推進(jìn)技術(shù),2015,36(6):839-845.

    ZHENG Y,LIH Y.Application ofγ-Reθtransition model in hypersonic flow based on new correlation equation[J].Journal of Propulsion Technology,2015,36(6):839-845(in Chinese).

    [31]WANG T,MATTHEW C R.Effect of elevated free-stream turbulence on transitional flow heat transfer over dual-scaled rough surfaces[J].Journal of Heat Transfer,2005,127(4):393-403.

    [32]PINSON M W,WANG T.Effect of two-scale roughness on boundary layer transition over a heated flat plate:Part2-Boundary layer structure[J].Journal of Turbomachinery,2000,122(2):308-316.

    [33]STRIPF M.Surface roughness effects on external heat transfer of a HP turbine vane[J].Journal of Turbomachinery,2005,127(1):200-208.

    [34]HYLTON L D,MIHELC M S,TURNER E R,et al.Analytical and experimental evaluation of the heat transfer distribution over the surfaces of turbine vanes[C]∥AAS/Division of Dynamical Astronomy Meeting.1983.

    Tel.:010-82338753

    E-mail:zheng_yun@buaa.edu.cn

    Effect of surface roughness on flow transition and heat transfer of turbine blade

    LIHongyang1,2,ZHENG Yun1,2,*

    (1.School of Energy and Power Engineering,Beijing University of Aeronautics and Astronautics,Beijing 100083,China;2.Collaborative Innovation Center for Advanced Aero-Engine,Beijing 100083,China)

    For the purpose of researching the effect of surface roughness on flow transition and heat transfer of turbine blade,a roughnessmodification method forγ-Reθtransition model was proposed through self-developed CFD code.Verification was conducted referring to the experimental data of flat plate and turbine vane cases,and satisfactory results were obtained.Taking surface roughness effect into consideration,No.5411 working condition of MarkⅡturbine vane was simulated and the result was analyzed in detail.Main conclusions are as follows:surface roughness has little effect on heat transfer of laminar boundary layer,while it has considerable effect on turbulent boundary layer.Compared with smooth surface,5μm equivalent sand roughness increases the suction side wall temperature by about 5.7 K in turbulent boundary layer,while 100μm roughness increases the temperature by about 28.4 K,reaching an increase of 5%.Under low roughness degree,effect of shock wave on boundary layer transition process of suction side plays a dominant role,while after reaching a critical degree,effect of surface roughness abruptly changes the transition position,and the critical degree is around 150μm in the current case.

    roughness;transition;turbine;interm ittency factor;boundary layer

    2015-10-13;Accep ted:2016-01-15;Pub lished online:2016-01-18 15:37 URL:www.cnki.net/kcms/detail/11.2625.V.20160118.1537.002.htm l

    V211.3

    A

    1001-5965(2016)10-2038-10

    李虹楊 男,博士研究生。主要研究方向:非定常流動及換熱的數(shù)值模擬,流、熱耦合數(shù)值模擬。

    E-mail:buaalihy@hotmail.com

    鄭赟 男,博士,講師。主要研究方向:計(jì)算流體力學(xué),葉輪機(jī)械流、熱、固耦合仿真。

    http:∥bhxb.buaa.edu.cn jbuaa@buaa.edu.cn

    DO I:10.13700/j.bh.1001-5965.2015.0659

    2015-10-13;錄用日期:2016-01-15;網(wǎng)絡(luò)出版時間:2016-01-18 15:37

    www.cnki.net/kcms/detail/11.2625.V.20160118.1537.002.htm l

    *通訊作者:Tel.:010-82338753 E-mail:zheng_yun@buaa.edu.cn

    李虹楊,鄭赟.粗糙度對渦輪葉片流動轉(zhuǎn)捩及傳熱特性的影響[J].北京航空航天大學(xué)學(xué)報,2016,42(10):2038-2047. LIH Y,ZHENG Y.Effect ofsurface roughness on flow transition and heat transfer of turbine blade[J].Journal ofBeijing University of Aeronautics and Astronautics,2016,42(10):2038-2047(in Chinese).

    *Correspond ing au thor.Tel.:010-82338753 E-mail:zheng_yun@buaa.edu.cn

    猜你喜歡
    模型
    一半模型
    一種去中心化的域名服務(wù)本地化模型
    適用于BDS-3 PPP的隨機(jī)模型
    提煉模型 突破難點(diǎn)
    函數(shù)模型及應(yīng)用
    p150Glued在帕金森病模型中的表達(dá)及分布
    函數(shù)模型及應(yīng)用
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    3D打印中的模型分割與打包
    亚洲欧美日韩东京热| 最近2019中文字幕mv第一页| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲最大成人手机在线| 色综合站精品国产| 国产精品av视频在线免费观看| 亚洲国产精品专区欧美| 少妇的逼好多水| 亚洲国产欧美在线一区| 亚洲欧美成人综合另类久久久| 2022亚洲国产成人精品| 夫妻午夜视频| 久久久久精品久久久久真实原创| 久久99热这里只有精品18| 亚洲精品成人久久久久久| 美女主播在线视频| 亚洲美女视频黄频| 99久久中文字幕三级久久日本| 爱豆传媒免费全集在线观看| 国产成人精品婷婷| 久久综合国产亚洲精品| 97热精品久久久久久| 日日摸夜夜添夜夜爱| 久久久久精品性色| 亚洲精华国产精华液的使用体验| 日韩精品青青久久久久久| 高清av免费在线| 国产精品熟女久久久久浪| 日韩欧美国产在线观看| 日本免费在线观看一区| 亚洲内射少妇av| 久久久久九九精品影院| 联通29元200g的流量卡| 在线观看av片永久免费下载| 最近视频中文字幕2019在线8| 欧美xxxx性猛交bbbb| 免费观看精品视频网站| 男女啪啪激烈高潮av片| 听说在线观看完整版免费高清| 成人性生交大片免费视频hd| 人人妻人人看人人澡| ponron亚洲| 麻豆国产97在线/欧美| 在线 av 中文字幕| 午夜免费观看性视频| 成人毛片60女人毛片免费| 国产成年人精品一区二区| 亚洲国产最新在线播放| 亚洲精品一区蜜桃| 全区人妻精品视频| xxx大片免费视频| 亚洲最大成人手机在线| 国产伦精品一区二区三区四那| 欧美精品一区二区大全| 有码 亚洲区| 国产精品一区www在线观看| 成人综合一区亚洲| 亚洲av国产av综合av卡| 69人妻影院| 亚洲精品中文字幕在线视频 | 日韩一区二区视频免费看| 中文字幕人妻熟人妻熟丝袜美| 国产激情偷乱视频一区二区| 精品酒店卫生间| 天天躁夜夜躁狠狠久久av| av国产免费在线观看| 久久人人爽人人爽人人片va| 天美传媒精品一区二区| 日韩亚洲欧美综合| 亚洲精品中文字幕在线视频 | 免费人成在线观看视频色| 免费大片黄手机在线观看| 国内少妇人妻偷人精品xxx网站| 国产成人精品婷婷| 天堂俺去俺来也www色官网 | av在线天堂中文字幕| 欧美3d第一页| 岛国毛片在线播放| 最近最新中文字幕免费大全7| 一区二区三区乱码不卡18| 乱码一卡2卡4卡精品| 日本免费在线观看一区| 色综合站精品国产| 狠狠精品人妻久久久久久综合| 成年版毛片免费区| 国产 亚洲一区二区三区 | 一区二区三区免费毛片| 中文字幕久久专区| 久久99精品国语久久久| 国产高清有码在线观看视频| 国产色爽女视频免费观看| 肉色欧美久久久久久久蜜桃 | 国产精品一区www在线观看| 波多野结衣巨乳人妻| 日产精品乱码卡一卡2卡三| 欧美一区二区亚洲| 伊人久久国产一区二区| 国内少妇人妻偷人精品xxx网站| 91久久精品国产一区二区三区| 国产中年淑女户外野战色| 国产精品不卡视频一区二区| 国产精品av视频在线免费观看| 最后的刺客免费高清国语| 十八禁国产超污无遮挡网站| 国产极品天堂在线| 色播亚洲综合网| 精品亚洲乱码少妇综合久久| 亚洲精品乱码久久久久久按摩| 欧美性感艳星| 高清毛片免费看| 国产熟女欧美一区二区| 综合色av麻豆| 丝瓜视频免费看黄片| 精品熟女少妇av免费看| 最近视频中文字幕2019在线8| 亚洲精品456在线播放app| 欧美不卡视频在线免费观看| 少妇的逼好多水| 三级国产精品片| 内射极品少妇av片p| 国内揄拍国产精品人妻在线| 男女国产视频网站| 免费不卡的大黄色大毛片视频在线观看 | 久久久久久久亚洲中文字幕| 黄色欧美视频在线观看| 成人午夜精彩视频在线观看| av福利片在线观看| 欧美日韩在线观看h| 天堂√8在线中文| 午夜福利在线观看吧| 一区二区三区高清视频在线| 亚洲av男天堂| 免费观看性生交大片5| 国产美女午夜福利| 国产成人a区在线观看| 成人高潮视频无遮挡免费网站| 午夜日本视频在线| 亚洲av成人av| 亚洲va在线va天堂va国产| 如何舔出高潮| 亚洲精品一二三| 国内少妇人妻偷人精品xxx网站| 97超视频在线观看视频| 久久久久九九精品影院| 欧美激情国产日韩精品一区| 国产精品一区二区三区四区久久| 91精品国产九色| 夫妻性生交免费视频一级片| 久久人人爽人人爽人人片va| 亚洲av电影不卡..在线观看| 纵有疾风起免费观看全集完整版 | 国产色爽女视频免费观看| 日韩av不卡免费在线播放| 国产精品女同一区二区软件| 国产高清三级在线| 久久国产乱子免费精品| 国模一区二区三区四区视频| 尾随美女入室| 免费电影在线观看免费观看| 最近2019中文字幕mv第一页| 精品国产一区二区三区久久久樱花 | 亚洲精品一二三| 最新中文字幕久久久久| 三级毛片av免费| 日韩,欧美,国产一区二区三区| 我的女老师完整版在线观看| 秋霞在线观看毛片| 精品国产露脸久久av麻豆 | 国产色爽女视频免费观看| 午夜日本视频在线| 一本久久精品| 十八禁网站网址无遮挡 | 国产伦理片在线播放av一区| 日韩欧美三级三区| 韩国高清视频一区二区三区| 精品熟女少妇av免费看| 日日摸夜夜添夜夜爱| 欧美一级a爱片免费观看看| 国产黄频视频在线观看| 在现免费观看毛片| 少妇高潮的动态图| 国产精品日韩av在线免费观看| 国内精品宾馆在线| 亚洲一区高清亚洲精品| 亚洲综合色惰| 国产精品一及| 激情五月婷婷亚洲| 久久午夜福利片| 亚洲国产精品成人久久小说| 国产淫语在线视频| 插逼视频在线观看| 久久久久久伊人网av| 在线观看免费高清a一片| 欧美成人精品欧美一级黄| 久久久久免费精品人妻一区二区| 中文字幕av成人在线电影| 777米奇影视久久| 禁无遮挡网站| 亚洲天堂国产精品一区在线| 国产av码专区亚洲av| 在线免费观看的www视频| 亚洲欧美清纯卡通| 丝瓜视频免费看黄片| 日韩人妻高清精品专区| 丰满少妇做爰视频| 日日摸夜夜添夜夜添av毛片| xxx大片免费视频| 蜜桃亚洲精品一区二区三区| 久久久久久九九精品二区国产| 日本爱情动作片www.在线观看| 夜夜爽夜夜爽视频| 22中文网久久字幕| 一级毛片电影观看| 久久久久性生活片| 午夜免费观看性视频| 精品人妻偷拍中文字幕| 2018国产大陆天天弄谢| 国产成人aa在线观看| 久久久久久久久中文| 亚洲精品自拍成人| 国产视频首页在线观看| 青春草国产在线视频| 亚洲欧美成人精品一区二区| 久久99精品国语久久久| 最近中文字幕高清免费大全6| 亚洲av电影不卡..在线观看| 少妇高潮的动态图| 国产精品久久久久久精品电影| 午夜免费激情av| 国产日韩欧美在线精品| 久久久亚洲精品成人影院| 国产视频首页在线观看| 一级av片app| 国产成人精品久久久久久| 久久精品久久精品一区二区三区| 欧美变态另类bdsm刘玥| 亚洲精品色激情综合| 欧美成人精品欧美一级黄| 在线观看av片永久免费下载| 亚洲av免费高清在线观看| 黄色配什么色好看| 最近视频中文字幕2019在线8| 丝袜喷水一区| 亚洲伊人久久精品综合| 亚洲国产成人一精品久久久| 国产午夜精品久久久久久一区二区三区| 亚洲av免费高清在线观看| 黄色欧美视频在线观看| 久久久亚洲精品成人影院| 婷婷六月久久综合丁香| 麻豆久久精品国产亚洲av| 久久韩国三级中文字幕| 淫秽高清视频在线观看| 国产精品爽爽va在线观看网站| 亚洲精品色激情综合| 亚洲国产av新网站| 久久99热这里只有精品18| 精品一区二区免费观看| 国内揄拍国产精品人妻在线| 久久这里只有精品中国| 99久久精品热视频| 综合色av麻豆| 麻豆国产97在线/欧美| 久久精品综合一区二区三区| a级毛片免费高清观看在线播放| av在线天堂中文字幕| 男女视频在线观看网站免费| 在线免费十八禁| 99久国产av精品| 国产综合懂色| 精品不卡国产一区二区三区| 精品午夜福利在线看| 狂野欧美白嫩少妇大欣赏| 亚洲欧美成人精品一区二区| 亚洲熟女精品中文字幕| 亚洲aⅴ乱码一区二区在线播放| 亚洲国产精品国产精品| 国产精品久久久久久久电影| 国产亚洲91精品色在线| 建设人人有责人人尽责人人享有的 | 丰满乱子伦码专区| 欧美日韩亚洲高清精品| 色综合色国产| 国产 一区精品| 国产亚洲最大av| 亚洲国产精品成人久久小说| 午夜免费观看性视频| av在线观看视频网站免费| 成人午夜高清在线视频| 免费黄网站久久成人精品| 全区人妻精品视频| 色哟哟·www| 亚洲成人一二三区av| 成人亚洲精品一区在线观看 | 国产91av在线免费观看| 国产精品久久久久久久电影| 国国产精品蜜臀av免费| 国产69精品久久久久777片| 国产黄片美女视频| 久久人人爽人人爽人人片va| 国产视频内射| 尾随美女入室| 三级国产精品欧美在线观看| 精品一区二区免费观看| 搡老妇女老女人老熟妇| 美女高潮的动态| 97超碰精品成人国产| 深夜a级毛片| 91在线精品国自产拍蜜月| 欧美极品一区二区三区四区| 日韩在线高清观看一区二区三区| 欧美精品一区二区大全| 国产精品国产三级国产av玫瑰| 国产av国产精品国产| 韩国高清视频一区二区三区| 看十八女毛片水多多多| 国产一区亚洲一区在线观看| 免费在线观看成人毛片| 国产毛片a区久久久久| 国产精品一及| 亚洲欧美一区二区三区国产| 久久精品国产亚洲av涩爱| 国产成人免费观看mmmm| 在线免费观看不下载黄p国产| 欧美激情久久久久久爽电影| 国产精品国产三级国产av玫瑰| 97精品久久久久久久久久精品| 一个人看视频在线观看www免费| 精品一区二区免费观看| 国产精品熟女久久久久浪| 成人高潮视频无遮挡免费网站| 深爱激情五月婷婷| 国产伦一二天堂av在线观看| 在线观看美女被高潮喷水网站| 国产又色又爽无遮挡免| 国产欧美另类精品又又久久亚洲欧美| 色视频www国产| 国产一区二区三区综合在线观看 | 国产不卡一卡二| 观看美女的网站| 麻豆av噜噜一区二区三区| 国产成人精品福利久久| 日韩制服骚丝袜av| av国产久精品久网站免费入址| 九九爱精品视频在线观看| 国产精品一区二区性色av| 777米奇影视久久| 日韩一区二区三区影片| 中文字幕av在线有码专区| 成人毛片60女人毛片免费| 亚洲国产欧美在线一区| 午夜激情欧美在线| 欧美+日韩+精品| 最近最新中文字幕免费大全7| 精品久久久精品久久久| 蜜桃久久精品国产亚洲av| 十八禁网站网址无遮挡 | 老司机影院毛片| 少妇被粗大猛烈的视频| 观看美女的网站| 国产精品国产三级国产专区5o| 秋霞伦理黄片| 精品不卡国产一区二区三区| eeuss影院久久| 亚洲av一区综合| 在线免费十八禁| 丝袜喷水一区| 亚洲成人av在线免费| 久久精品国产自在天天线| 色综合站精品国产| 少妇熟女欧美另类| 免费av不卡在线播放| 国产亚洲午夜精品一区二区久久 | 中文字幕人妻熟人妻熟丝袜美| 国产黄色免费在线视频| 一区二区三区免费毛片| 午夜免费男女啪啪视频观看| 婷婷色综合大香蕉| 美女cb高潮喷水在线观看| videos熟女内射| 日韩精品青青久久久久久| 性插视频无遮挡在线免费观看| 午夜福利成人在线免费观看| 午夜福利在线观看免费完整高清在| 亚洲精品国产成人久久av| 精品亚洲乱码少妇综合久久| 国产三级在线视频| 日本欧美国产在线视频| 国产精品.久久久| 久久久久久国产a免费观看| 久久久久久久国产电影| 久久精品国产亚洲av涩爱| 久久精品久久久久久噜噜老黄| 最后的刺客免费高清国语| 青春草视频在线免费观看| 国产精品无大码| 久久久久网色| 黄色日韩在线| 美女xxoo啪啪120秒动态图| 男女那种视频在线观看| 国产一区二区三区av在线| 禁无遮挡网站| 波多野结衣巨乳人妻| 人人妻人人澡欧美一区二区| 偷拍熟女少妇极品色| 欧美 日韩 精品 国产| 九九在线视频观看精品| 2018国产大陆天天弄谢| 一级av片app| 国产老妇女一区| 成人特级av手机在线观看| 99久国产av精品| 国产高潮美女av| 丝袜美腿在线中文| 久久久精品免费免费高清| 日本wwww免费看| 免费在线观看成人毛片| 大陆偷拍与自拍| 欧美日韩视频高清一区二区三区二| 成人美女网站在线观看视频| 久久精品国产自在天天线| 日本免费a在线| 国内揄拍国产精品人妻在线| 一本久久精品| 成人亚洲精品一区在线观看 | 亚洲人成网站高清观看| 久久99热这里只有精品18| 99久国产av精品| 久久精品人妻少妇| 丰满少妇做爰视频| 国产一区二区三区综合在线观看 | 插逼视频在线观看| 性插视频无遮挡在线免费观看| 久久久久国产网址| 舔av片在线| 夜夜爽夜夜爽视频| 97超视频在线观看视频| 久久精品国产鲁丝片午夜精品| 日韩精品有码人妻一区| 亚洲成人一二三区av| 汤姆久久久久久久影院中文字幕 | 精品久久久久久久人妻蜜臀av| 免费电影在线观看免费观看| 中文字幕制服av| 人人妻人人澡欧美一区二区| 人体艺术视频欧美日本| 国产午夜精品久久久久久一区二区三区| 日日摸夜夜添夜夜添av毛片| 自拍偷自拍亚洲精品老妇| 一个人观看的视频www高清免费观看| 亚洲精品一二三| 五月玫瑰六月丁香| 精品99又大又爽又粗少妇毛片| 精品人妻偷拍中文字幕| 18禁裸乳无遮挡免费网站照片| 岛国毛片在线播放| 久久综合国产亚洲精品| 秋霞伦理黄片| 菩萨蛮人人尽说江南好唐韦庄| 毛片女人毛片| 免费看不卡的av| av免费在线看不卡| 免费看日本二区| 亚洲精品亚洲一区二区| 亚洲av免费高清在线观看| 九九在线视频观看精品| 六月丁香七月| 深夜a级毛片| 搡老妇女老女人老熟妇| 国产久久久一区二区三区| 亚洲国产精品sss在线观看| 成年人午夜在线观看视频 | 一个人看的www免费观看视频| 国产精品人妻久久久影院| 少妇熟女欧美另类| 免费av不卡在线播放| 天美传媒精品一区二区| 97超碰精品成人国产| 天堂√8在线中文| 中文字幕亚洲精品专区| 亚洲av二区三区四区| 日本wwww免费看| 国产在视频线精品| 大又大粗又爽又黄少妇毛片口| 亚洲精品日韩在线中文字幕| 欧美性感艳星| 一级爰片在线观看| 精品人妻视频免费看| 久99久视频精品免费| 国产免费福利视频在线观看| 丰满少妇做爰视频| 女人久久www免费人成看片| 搡老乐熟女国产| 美女被艹到高潮喷水动态| 丰满人妻一区二区三区视频av| 久久精品国产鲁丝片午夜精品| 亚洲熟妇中文字幕五十中出| 久久久久免费精品人妻一区二区| 毛片一级片免费看久久久久| 人妻夜夜爽99麻豆av| 建设人人有责人人尽责人人享有的 | 亚洲不卡免费看| 久久人人爽人人爽人人片va| 麻豆成人av视频| videossex国产| ponron亚洲| 免费看日本二区| 一级毛片 在线播放| 成人午夜精彩视频在线观看| 亚洲欧美日韩卡通动漫| 肉色欧美久久久久久久蜜桃 | 伦精品一区二区三区| 久久午夜福利片| 成人毛片a级毛片在线播放| 又爽又黄无遮挡网站| 99热网站在线观看| 小蜜桃在线观看免费完整版高清| 一区二区三区四区激情视频| 亚洲国产高清在线一区二区三| 亚洲精品第二区| 97人妻精品一区二区三区麻豆| 国产亚洲av片在线观看秒播厂 | 天堂网av新在线| 精品久久久久久久久av| 国产69精品久久久久777片| 久久精品人妻少妇| a级毛色黄片| 精品人妻偷拍中文字幕| 亚洲国产精品成人久久小说| 成年女人看的毛片在线观看| 日本欧美国产在线视频| 亚洲欧美一区二区三区黑人 | 麻豆成人午夜福利视频| 97人妻精品一区二区三区麻豆| 国产成人a区在线观看| av在线天堂中文字幕| 九九在线视频观看精品| 麻豆av噜噜一区二区三区| 最近最新中文字幕免费大全7| 尾随美女入室| 毛片女人毛片| 蜜臀久久99精品久久宅男| 国产在视频线在精品| 综合色av麻豆| 熟女电影av网| 又爽又黄无遮挡网站| 国产精品一区二区三区四区久久| 成人高潮视频无遮挡免费网站| 精品一区二区三区视频在线| 久久精品久久久久久噜噜老黄| 久久这里有精品视频免费| 少妇猛男粗大的猛烈进出视频 | 激情 狠狠 欧美| 人体艺术视频欧美日本| 麻豆精品久久久久久蜜桃| 美女国产视频在线观看| 国产成人一区二区在线| 2021少妇久久久久久久久久久| 少妇人妻一区二区三区视频| 青春草视频在线免费观看| 欧美激情在线99| 国产乱人视频| 直男gayav资源| 久久久久久久午夜电影| 尾随美女入室| 18禁在线播放成人免费| 特级一级黄色大片| 在线观看美女被高潮喷水网站| 国产伦在线观看视频一区| 男女国产视频网站| 天堂影院成人在线观看| 菩萨蛮人人尽说江南好唐韦庄| 伊人久久国产一区二区| 亚洲综合色惰| 精品一区二区三卡| 狂野欧美激情性xxxx在线观看| 91狼人影院| 亚洲aⅴ乱码一区二区在线播放| 免费av不卡在线播放| 亚洲av免费高清在线观看| 男女啪啪激烈高潮av片| 亚洲精品中文字幕在线视频 | 草草在线视频免费看| 又黄又爽又刺激的免费视频.| 男女边吃奶边做爰视频| 亚洲性久久影院| 国产精品久久久久久av不卡| 亚洲欧美一区二区三区黑人 | 婷婷六月久久综合丁香| 亚洲国产日韩欧美精品在线观看| 日本欧美国产在线视频| 日韩精品有码人妻一区| 欧美日本视频| 久久久色成人| 亚洲国产av新网站| a级一级毛片免费在线观看| 午夜福利在线观看吧| 亚洲自拍偷在线| 美女被艹到高潮喷水动态| 最后的刺客免费高清国语| 免费观看a级毛片全部| 波多野结衣巨乳人妻| 禁无遮挡网站| 久久精品夜夜夜夜夜久久蜜豆| 国产又色又爽无遮挡免| 日本一二三区视频观看| 亚洲精品成人久久久久久| 人人妻人人看人人澡| 99久国产av精品国产电影| 热99在线观看视频| 高清视频免费观看一区二区 | 美女主播在线视频| 免费观看精品视频网站| 日韩成人av中文字幕在线观看| 国产精品一区二区在线观看99 |