亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        壓縮性修正對γ-Reθ轉(zhuǎn)捩模型的影響研究

        2015-03-28 11:06:40夏陳超姜婷婷郭中州陳偉芳
        空氣動力學(xué)學(xué)報 2015年5期
        關(guān)鍵詞:方法模型

        夏陳超,姜婷婷,郭中州,陳偉芳

        (浙江大學(xué)航空航天學(xué)院,浙江杭州 310027)

        壓縮性修正對γ-Reθ轉(zhuǎn)捩模型的影響研究

        夏陳超,姜婷婷,郭中州,陳偉芳*

        (浙江大學(xué)航空航天學(xué)院,浙江杭州 310027)

        將基于流場當(dāng)?shù)刈兞康摩?Reθ轉(zhuǎn)捩模型加入可壓縮RANS方程計算程序,對平均流場控制方程、湍流和轉(zhuǎn)捩控制方程進行基于LU-SGS的隱式緊耦合求解。結(jié)合三種壓縮性修正方法對高超聲速平板、雙楔和圓錐流動進行了數(shù)值模擬,給出了壁面斯坦頓數(shù)、熱流值和實驗值的比較,以及湍動能和間歇因子分布云圖。數(shù)值計算結(jié)果表明,相同來流湍流度下不同壓縮性修正方法得到的轉(zhuǎn)捩位置差別較大。相比于未修正的模型,基于當(dāng)?shù)伛R赫數(shù)壓縮性修正后的轉(zhuǎn)捩位置最靠后,其對帶有分離的雙楔流動預(yù)測較為準(zhǔn)確;基于來流馬赫數(shù)的壓縮性修正僅使得轉(zhuǎn)捩位置稍有延遲;對湍動能源項的壓縮性修正可提高模型在轉(zhuǎn)捩后湍流段的熱流預(yù)測精度。

        高超聲速;RANS方程;γ-Reθ轉(zhuǎn)捩模型;緊耦合;壓縮性修正

        0 引言

        邊界層轉(zhuǎn)捩是目前流體力學(xué)領(lǐng)域的研究熱點和難點,也是航空航天工程中的關(guān)鍵問題。雖然國內(nèi)外研究者已經(jīng)做了大量卓有成效的工作,但對轉(zhuǎn)捩相關(guān)機理的認(rèn)識和理論分析仍有待提升。對高超聲速飛行器而言,流動形態(tài)的變化關(guān)系到飛行器的升阻特性、熱防護以及發(fā)動機中的流動品質(zhì)等。準(zhǔn)確預(yù)測轉(zhuǎn)捩位置、轉(zhuǎn)捩區(qū)長度以及氣動力/熱環(huán)境對飛行器設(shè)計有重要意義。

        相比于直接數(shù)值模擬(DNS)和大渦模擬(LES),基于RANS的數(shù)值模擬是當(dāng)前工程上預(yù)測湍流和轉(zhuǎn)捩較為經(jīng)濟和高效的方法,其中Huang[1]、Walters[2]和Lantry[3]等的研究具有代表性。Lantry等提出的基于經(jīng)驗關(guān)系式的γ-Reθ轉(zhuǎn)捩模型在k-ω SST兩方程湍流模型的基礎(chǔ)上添加了兩個輸運方程,一個是間歇因子γ的輸運方程,一個是動量厚度雷諾數(shù)Reθ的輸運方程。γ-Reθ模型是基于流場當(dāng)?shù)刈兞康霓D(zhuǎn)捩模型,可以和現(xiàn)代CFD程序相兼容,相比于目前常用的eN方法[4],具有更強的適用性。針對低速平板、翼型/機翼和渦輪機械的眾多研究[5-7]表明γ-Reθ模型對自然轉(zhuǎn)捩和分離轉(zhuǎn)捩具有較好的預(yù)測能力。

        一方面,雖然原始的γ-Reθ模型基于不可壓實驗發(fā)展而來,但 Krause[8]、You[9]、孔維萱[10]以及鄭赟[11]等采用γ-Reθ模型對高超聲速雙楔和尖錐等流動的數(shù)值計算表明該模型也能較準(zhǔn)確地預(yù)測高超聲速流動中的轉(zhuǎn)捩問題。另一方面,為增強該模型對高速流動問題的適用性,有必要進行適當(dāng)?shù)膲嚎s性修正。Kaynak[12]給出了一種基于來流馬赫數(shù)的壓縮性修正關(guān)系式,對超聲速平板的邊界層轉(zhuǎn)捩進行了預(yù)測;張曉東[13]則對經(jīng)驗關(guān)系式進行了基于當(dāng)?shù)伛R赫數(shù)的修正,通過對雙楔模型的數(shù)值計算,得到的壁面壓力系數(shù)與實驗值吻合較好。總體而言,該模型目前主要應(yīng)用于低速流動問題,其在高速流動情況下的壓縮性修正研究還很少。本文將γ-Reθ模型編入課題組發(fā)展的基于SST湍流模型的可壓縮RANS方程耦合求解程序,結(jié)合高超聲速平板和圓錐等典型流動問題,通過數(shù)值模擬和對比,初步探索了三種壓縮性修正方法對γ-Reθ模型預(yù)測轉(zhuǎn)捩的影響,為該模型的進一步改進提供參考。

        1 轉(zhuǎn)捩模型

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

        γ-Reθ轉(zhuǎn)捩模型包含間歇因子和轉(zhuǎn)捩動量厚度雷諾數(shù)兩個輸運方程,其中間歇因子的輸運方程為:

        其中,F(xiàn)length為轉(zhuǎn)捩區(qū)長度控制函數(shù),F(xiàn)onset為轉(zhuǎn)捩發(fā)生函數(shù),F(xiàn)turb用于抑制再層流化,相關(guān)參數(shù)詳見文獻[3]。

        1.2 壓縮性修正

        由于原始的γ-Reθ模型中轉(zhuǎn)捩動量厚度雷諾數(shù)Reθt關(guān)系式是通過不可壓平板實驗獲得,因此將該模型應(yīng)用于高速可壓縮流問題時有必要進行適當(dāng)?shù)目蓧嚎s性修正(Compressibility Correction,文中用CC表示)。本文實現(xiàn)并對比了三種壓縮性修正方法。

        (1)CC=1

        張曉東[13]基于高超聲速風(fēng)洞實驗數(shù)據(jù),給出了一種基于當(dāng)?shù)伛R赫數(shù)的壓縮性修正關(guān)系式:

        (2)CC=2

        Kaynak[12]給出了一種基于來流馬赫數(shù)的修正關(guān)系式:

        雖然文獻[12]給出了該關(guān)系式的適用范圍(來流馬赫數(shù)小于2.4),本文仍將其視為一種潛在的修正方法,評估其壓縮性修正效果。

        (3)CC=3

        Papp等人[14]在對其發(fā)展的基于SSGZ k-ε模型的工程轉(zhuǎn)捩模型進行壓縮性修正時在湍動能方程中添加了一個源項,形式如下:

        本文參考該方法,將上述源項加入SST湍流模型的湍動能方程中進行修正。其中,,湍流馬赫數(shù),M't為考慮轉(zhuǎn)捩效應(yīng)的轉(zhuǎn)捩湍流馬赫數(shù),系數(shù)λ的作用是當(dāng)轉(zhuǎn)捩湍流馬赫數(shù)小于λ時不起作用,取值0.2,α1=2.5,α2=2.0,ε=βωk。

        壓縮性修正1和2的實施是將F(Ma)乘到轉(zhuǎn)捩動量厚度雷諾數(shù),壓縮性修正3則是將Pe直接加到湍動能k方程的源項中。

        2 數(shù)值方法

        將 γ-Reθ轉(zhuǎn)捩模型加入課題組發(fā)展的可壓縮RANS程序,流動控制方程為雷諾平均的三維NS方程、SST兩方程湍流模型方程以及γ-Reθ轉(zhuǎn)捩輸運方程。計算程序基于結(jié)構(gòu)網(wǎng)格的有限體積方法,其中無黏通量采用計算效率和分辨率均較高的AUSMPW+格式[15]并采用MUSCL進行高階插值,黏性通量采用中心差分進行離散。

        為提高計算效率,時間推進采用Yoon和Jameson[16]提出的LU-SGS隱式方法,同時采用當(dāng)?shù)貢r間步長加速收斂。為克服源項剛性,對湍流和轉(zhuǎn)捩輸運方程中的負值耗散項進行隱式處理以增強矩陣對角占優(yōu)[12]:

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

        3.1 平板流動

        高超聲速平板流動算例取自Mee[17]的激波風(fēng)洞實驗結(jié)果。平板長1.5m,來流馬赫數(shù)6.2,來流靜溫690K,來流靜壓5.4kPa,單位雷諾數(shù)為2.6×106。計算網(wǎng)格為224×150(流向和法向),在壁面和前緣進行了加密,壁面第一層網(wǎng)格y+小于1(約0.5)。

        圖1和圖2為兩種來流湍流度情況下的平板壁面斯坦頓數(shù)分布,斯坦頓數(shù)定義為:

        其中,qw為熱流,ρ∞和u∞分別為來流密度和速度,h0∞為來流總焓,hw壁面焓。圖中對比了無壓縮性修正的γ-Reθ模型、三種壓縮性修正的γ-Reθ模型、完全層流狀態(tài)和完全湍流狀態(tài)(SST模型)的結(jié)果??梢钥闯?,層流和湍流的斯坦頓數(shù)與實驗值吻合較好,轉(zhuǎn)捩模型實現(xiàn)了層流到湍流的過渡,隨著來流湍流度的增大,轉(zhuǎn)捩位置前移,但相比于實驗值,轉(zhuǎn)捩區(qū)長度過長。相同湍流度條件下不同壓縮性修正方法的轉(zhuǎn)捩起始位置差異較大,方法1轉(zhuǎn)捩位置最靠后,方法2居中,相比于無壓縮性模型均延遲了轉(zhuǎn)捩位置,方法3由于沒有改變轉(zhuǎn)捩模型經(jīng)驗關(guān)系式,其轉(zhuǎn)捩位置與無壓縮性模型基本一致,但其在一定程度上降低了湍流段的熱流。

        圖3和圖4分別為湍流度3%時無壓縮性修正和采用壓縮性方法2的流場湍動能云圖,湍動能的急劇增大意味著轉(zhuǎn)捩的發(fā)生,可以看出,進行壓縮性修正后轉(zhuǎn)捩位置有了明顯的后移,即該壓縮性修正方法限制了湍動能的發(fā)展。數(shù)值計算也表明,無壓縮性修正時最大湍動能為61 611m2/s2,壓縮性修正方法2得到的最大湍動能減少到了60 257 m2/s2。

        圖1 湍流度2%的平板壁面斯坦頓數(shù)Fig.1 Stanton number on flat plate(Tu=2%)

        圖2 湍流度3%的平板壁面斯坦頓數(shù)Fig.2 Stanton number on flat plate(Tu=3%)

        圖3 無壓縮性修正的平板湍動能云圖Fig.3 Contour of turbulence kinetic energy without compressibility correction

        圖4 壓縮性方法2的平板湍動能云圖Fig.4 Contour of turbulence kinetic energy with compressibility correction method 2

        3.2 雙楔流動

        雙楔幾何外形[18]如圖5所示,其第一個和第二個坡的角度分別為 9°和 20.5°,前緣鈍化半徑0.5mm。計算馬赫數(shù)為8.1,來流單位雷諾數(shù)為3.8 ×106,來流靜壓520Pa,來流靜溫106K,壁面溫度為300K,來流湍流度0.9%。鈍頭雙楔計算網(wǎng)格如圖6所示,網(wǎng)格在壁面、頭部和拐角附近進行了加密,并保證第一層y+值小于1(約為0.3)。

        圖5 雙楔流動示意圖Fig.5 Schematic of blunted double wedge flow

        圖6 雙楔模型計算網(wǎng)格Fig.6 Grid of blunted double wedge

        圖7 雙楔模型壁面壓力分布Fig.7 Comparison of pressure coefficient of double wedge with blunt leading edge

        圖7 和圖8分別為雙楔模型壁面壓力系數(shù)和斯坦頓數(shù)分布??梢钥闯?,完全層流狀態(tài)在拐角的分離區(qū)過大,斯坦頓數(shù)過低,完全湍流狀態(tài)則未發(fā)生分離,斯坦頓數(shù)過高。轉(zhuǎn)捩模型的結(jié)果與實驗值更接近,其中,壓縮性修正方法1的分離區(qū)大小、壓力和斯坦頓數(shù)與實驗值最為吻合,壓縮性修正方法3與未修正的轉(zhuǎn)捩模型相比,拐角后的湍流段壁面斯坦頓數(shù)較低。

        圖8 雙楔模型壁面斯坦頓數(shù)分布Fig.8 Comparison of Stanton number of double wedge with blunt leading edge

        圖9 至圖11給出了不同壓縮性方法在拐角處的局部放大流線圖和間歇因子云圖,可以看出,間歇因子分布和流動的分離是相對應(yīng)的,分離起始點之前的間歇因子保持很小的值,表明流動還是層流;分離之后,間歇因子迅速增長。從壓縮性修正方法1到3,分離區(qū)逐漸減小,壓縮性減弱,這與壁面壓力和熱流的分布趨勢是一致的。

        圖9 壓縮性方法1拐角處流線和間歇因子Fig.9 Streamline and intermittency factor around corner with compressibility correction method 1

        圖10 壓縮性方法2拐角處流線和間歇因子Fig.10 Streamline and intermittency factor around corner with compressibility correction method 2

        圖11 壓縮性方法3拐角處流線和間歇因子Fig.11 Streamline and intermittency factor around corner with compressibility correction method 3

        3.3 圓錐流動

        針對高超聲速圓錐流動[19]進行數(shù)值計算,圓錐錐角為7°,頭部鈍化半徑5 mm,如圖12所示。來流馬赫數(shù)7.15,來流靜溫214 K,壁溫300 K,來流靜壓7 722 Pa,單位雷諾數(shù)為9.64×106。采用半模計算,網(wǎng)格量為106×78×29(流向、法向和周向),壁面第一層網(wǎng)格間距為1×10-6m。

        圖12 圓錐幾何外形Fig.12 Geometry of cone

        由于不易實現(xiàn)計算與實驗條件的完全一致,本文計算了多個來流湍流度的情況,圖13至圖15分別為來流湍流度0.2%、0.35%和0.5%條件下的圓錐表面熱流分布對比。可以看出,完全層流的計算結(jié)果在層流段與實驗值吻合,完全湍流的計算結(jié)果在湍流段熱流比實驗值偏高,轉(zhuǎn)捩模型的結(jié)果與實驗值更為接近,且隨著來流湍流度的提高,轉(zhuǎn)捩位置前移。相同湍流度條件下三種壓縮性修正方法的轉(zhuǎn)捩位置相差較大,方法1轉(zhuǎn)捩位置最靠后(來流湍流度為0.2%時未發(fā)生轉(zhuǎn)捩),方法2居中,方法3轉(zhuǎn)捩位置與無壓縮性修正基本相同,趨勢與3.1中的平板一致。壓縮性修正方法3在轉(zhuǎn)捩后的湍流段熱量值較低,與實驗值更接近。

        圖13 來流湍流度0.2%壁面熱流分布Fig.13 Heat flux on wall(Tu=0.2%)

        圖14 來流湍流度0.35%壁面熱流分布Fig.14 Heat flux on wall(Tu=0.35%)

        圖15 來流湍流度0.5%壁面熱流分布Fig.15 Heat flux on wall(Tu=0.5%)

        4 結(jié)論

        基于可壓縮RANS方程計算程序,引入γ-Reθ轉(zhuǎn)捩模型,采用三種壓縮性修正方法對高超聲速平板、雙楔和圓錐進行了數(shù)值模擬并與實驗結(jié)果進行了比較。研究表明,三種壓縮性修正后的轉(zhuǎn)捩模型在相同湍流度條件下的轉(zhuǎn)捩位置相差較大。基于當(dāng)?shù)伛R赫數(shù)的壓縮性修正有效地預(yù)測了帶有分離的雙楔流動,同時,相比于未修正的轉(zhuǎn)捩模型,該修正會較大程度地延遲轉(zhuǎn)捩位置?;趤砹黢R赫數(shù)的壓縮性修正僅使得轉(zhuǎn)捩位置稍有延遲。湍動能方程源項的壓縮性修正能降低轉(zhuǎn)捩后湍流段壁面熱流,使結(jié)果與實驗值更加接近,對轉(zhuǎn)捩位置沒有影響。此外還可看出,對于不同的流動特性,相同的壓縮性修正方式效果是不一樣的,即很難得到對任何流動情況都適用的模型,實際情況中,應(yīng)根據(jù)不同的流動特性進行湍流和轉(zhuǎn)捩模型的選擇。

        總體而言,要將γ-Reθ轉(zhuǎn)捩模型應(yīng)用于高速流動問題還有許多值得改進的工作。相比于壁面壓力系數(shù),高超聲速條件下壁面熱流值預(yù)測精度的提高仍是一個挑戰(zhàn)。本文中研究的壓縮性修正方法均較為簡單,將當(dāng)?shù)伛R赫數(shù)和湍動能源項修正相結(jié)合可能是一個更好的嘗試??紤]更多轉(zhuǎn)捩影響因素、適用范圍更廣的壓縮性修正方法是今后的一個探索方向。

        參考文獻:

        [1]Suzen Y B,Huang P G.An intermittency transport equation for modeling flow transition[R].AIAA 2000-0287.2000

        [2]Walters D K,Leylek J H.A new model for boundary-layer transition using a single-point RANS approach[J].Journal of Turbomachinery,2002,126(1):193-202.

        [3]Langtry R B,Menter F R.Transition modeling for general CFD applications in aeronautics[C].Reno:43rd AIAA Aerospace Sciences Meeting and Exhibit,2005.

        [4]Malik M,Balakumar P.Instability and transition in three-dimensional supersonic boundary layers[R].AIAA-92-5049.1992

        [5]Malan P,Suluksna K,Juntasaro E.Calibrating the γ-Reθtransition model for commercial CFD[R].AIAA 2009-1142.2009

        [6]Langtry R B,Menter F R.Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes[J].AIAA Journal,2009,47(12):2894-2906.

        [7]Wang Gang,Liu Yi,Wang Guangqiu,et al.Transitional flow simulation based on γ-Reθttransition model[J].Acta Aeronautica et Astronautica Sinica,2014,35(1):70-79.(in Chinese)王剛,劉毅,王光秋,等.采用γ-Reθt模型的轉(zhuǎn)捩流動計算分析[J].航空學(xué)報.2014,35(1):70-79.

        [8]Krause M,Behr M,Ballmann J.Modeling of transition effects in hypersonic intake flows using a correlation-based intermittency model[C].Dayton:15th AIAA International Space Planes and Hypersonic Systems and Technologies Conference,2008.

        [9]You Y,Luedeke H,Eggers T,Hannemann K.Application of the γ-Reθttransition model in high speed flows[C].Tours:18th AIAA/3AF International Space Planes and Hypersonic Systems and Technologies Conference,2012.

        [10]Kong Weixuan,Yan Chao,Zhao Rui.γ-Reθmodel research for high-speed boundary layer transition[J].Acta Aerodynamica Sinica,2013,31(1):120-126.(in Chinese)孔維萱,閻超,趙瑞.γ-Reθ模式應(yīng)用于高速邊界層轉(zhuǎn)捩的研究[J].空氣動力學(xué)學(xué)報.2013,31(1):120-126.

        [11]Zheng Yun,Li Hongyang,Liu Daxiang.Application and analysis of γ-Reθtransition model in hypersonic flow[J].Journal of Propulsion Technology,2014,35(3):296-304.(in Chinese)鄭赟,李虹楊,劉大響.γ-Reθ轉(zhuǎn)捩模型在高超聲速下的應(yīng)用及分析[J].推進技術(shù),2014,35(3):296-304.

        [12]Kaynak U.Supersonic boundary-layer transition prediction under the effect of compressibility using a correlation-based model[J].Proceedings of the Institution of Mechanical Engineers,Part G:Journal of Aerospace Engineering,2012,226(7):722-739.

        [13]Zhang X D,Gao Z H.Numerical discuss to complete empirical correlation in Langtry's transition model[J].Applied Mathematics and Mechanics,2010,31(5):544-552.(in Chinese)張曉東,高正紅.關(guān)于補充Langtry的轉(zhuǎn)捩模型經(jīng)驗修正式的數(shù)值探討[J].應(yīng)用數(shù)學(xué)和力學(xué).2010,31(5):544-552.

        [14]Papp J L,Kenzakowski D C,Dash S M.Extensions of a rapid engineering approach to modeling hypersonic laminar to turbulent transitional flows[R].AIAA 2005-0892.

        [15]Kim K H,Kim C,Rho O.Methods for the accurate computations of hypersonic flows:I.AUSMPW+scheme[J].Journal of Computational Physics,2001,174(1):38-80.

        [16]Yoon S,Jameson A.Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations[J].AIAA Journal,1988,26(9):1025-1026.

        [17]Mee D J.Boundary-layer transition measurements in hypervelocity flows in a shock tunnel[J].AIAA Journal,2002,40(8):1542-1548.

        [18]Neuenhahn T,Olivier H.Influence of the wall temperature and the entropy layer effects on double wedge shock boundary layer interactions[R].AIAA 2006-8136.

        [19]Maclean M,Mundy E,Wadhams T,Holden M,Johnson H,Candler G.Comparisons of transition prediction using PSE-chem to measurements for a shock tunnel environment[R].AIAA 2007-4490.

        Effects of compressibility correction on γ-Reθtransition model

        Xia Chenchao,Jiang Tingting,Guo Zhongzhou,Chen Weifang*
        (School of Aeronautics and Astronautics,Zhejiang University,Hangzhou 310027,China)

        Boundary layer transition plays a significant role in the prediction of aerodynamic and aerothermodynamics characteristics of hypersonic vehicle.The commonly used correlation based on γ-Reθtransition model is implemented into a Reynolds-Averaged Navier-Stokes solver in order to evaluate its capability in hypersonic flows.The mean flow and turbulent flow equations are solved simultaneously based on the LU-SGS method.Three compressibility correction methods are adopted to simulate the hypersonic flows around a flat plate,a double wedge and a cone body.Results of Stanton number,heat flux on walls,contours of turbulence kinetic energy and intermittency factor are presented.The results of γ-Reθtransition model show better flow essential than that of full laminar or full turbulent results.The transition onset positions obtained by different compressibility methods under the same free stream turbulence intensity varied significantly.Compressibility correction based on local Mach number delays the transition onset greatly,while the correction of source term of turbulent kinetic energy equation improves the accuracy of heat flux on turbulent section of the wall.Rational selection of compressibility correction method should rely on the type of flow.Elaborate and reliable compressibility correction method still need to be further investigated.

        hypersonic;RANS equations;γ-Reθtransition model;coupled method;compressibility correction

        V211.3

        :Adoi:10.7638/kqdlxxb-2014.0068

        0258-1825(2015)05-0603-07

        2014-07-14;

        :2014-10-10

        國家重點基礎(chǔ)研究發(fā)展計劃(2014CB340201)

        夏陳超(1989-),男,福建寧德人,博士生,流體力學(xué)專業(yè).E-mail:aeroxia@zju.edu.cn

        陳偉芳*(1970-),男,博士,教授,研究方向為高超聲速空氣動力學(xué).E-mail:chenwfnudt@163.com

        夏陳超,姜婷婷,郭中州,等.壓縮性修正對γ-Reθ轉(zhuǎn)捩模型的影響研究[J].空氣動力學(xué)學(xué)報,2015,33(5):603-609.

        10.7638/kqdlxxb-2014.0068 Xia C C,Jiang T T,Guo Z Z,et al.Effects of compressibility correction on γ-Reθtransition model[J].Acta Aerodynamica Sinica,2015,33(5):603-609.

        猜你喜歡
        方法模型
        一半模型
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        學(xué)習(xí)方法
        可能是方法不對
        3D打印中的模型分割與打包
        用對方法才能瘦
        Coco薇(2016年2期)2016-03-22 02:42:52
        FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
        四大方法 教你不再“坐以待病”!
        Coco薇(2015年1期)2015-08-13 02:47:34
        賺錢方法
        免费观看的a级毛片的网站| 亚洲女av中文字幕一区二区| 伊人久久综合无码成人网| 最近中文字幕完整版免费| 六月丁香久久| 少妇高潮紧爽免费观看| 亚洲一区二区三区偷拍厕所| 乱中年女人伦av一区二区| 日本大片免费观看完整视频| 国产AV无码无遮挡毛片| 日韩美女av一区二区| 777米奇色狠狠俺去啦| 国产女合集小岁9三部| 中文字幕麻豆一区二区| 久久免费亚洲免费视频| 亚洲精品国偷拍自产在线观看| 国产精品高潮呻吟av久久无吗| 亚洲无码美韩综合| 亚洲av毛片在线免费看| 中文字幕人妻熟在线影院| 久久永久免费视频| 黑人一区二区三区啪啪网站| 国产av精品一区二区三区久久| 精品亚洲一区二区三区在线观看 | 久久精品性无码一区二区爱爱| 亚洲一区二区三区av天堂| 午夜熟女插插xx免费视频| 无码一区二区三区老色鬼| 四虎无码精品a∨在线观看| 人妻少妇精品视频专区二区三区| 国产乱了真实在线观看| 在线亚洲综合| 久久精品亚洲成在人线av| 噜噜综合亚洲av中文无码| 精品国产av 无码一区二区三区| 两个人免费视频大全毛片| 青青草视频在线观看网| 成人午夜福利视频镇东影视 | 亚洲免费不卡av网站| 亚洲三级视频一区二区三区| 天天躁日日躁狠狠躁av|