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

    基于溫度異重流模型的連續(xù)分層流數(shù)值模擬方法研究

    2017-08-27 05:36:11牛明昌丁勇馬衛(wèi)狀韓盼盼
    船舶力學 2017年8期
    關(guān)鍵詞:尾流層流流線

    牛明昌,丁勇,馬衛(wèi)狀,韓盼盼

    (1.哈爾濱工程大學船舶工程學院,哈爾濱150001;2.上海交通大學海洋工程國家重點實驗室,上海200240)

    基于溫度異重流模型的連續(xù)分層流數(shù)值模擬方法研究

    牛明昌1,丁勇1,馬衛(wèi)狀1,韓盼盼2

    (1.哈爾濱工程大學船舶工程學院,哈爾濱150001;2.上海交通大學海洋工程國家重點實驗室,上海200240)

    文章基于溫度異重流模型,通過指定流場中溫度的分布及密度與溫度的函數(shù)關(guān)系實現(xiàn)密度的連續(xù)分層,并假設(shè)流場中的熱傳導在數(shù)值過程中可以忽略不計以確保密度分層的穩(wěn)定。分別采用RANS方法和LES方法對低傅汝德數(shù)和高傅汝德數(shù)情況下的圓柱繞流進行數(shù)值模擬,其尾流場中的lee波、葫蘆狀流線以及高傅汝德數(shù)情況下的波包與實驗值相符,低傅汝德數(shù)情況下圓柱上下游速度剖面曲線的最小值與實驗值僅存在2%的誤差,且其分布規(guī)律相同。通過改變流體的導熱系數(shù),探討了Peclet數(shù)對該方法收斂穩(wěn)定性的影響,得到了收斂的數(shù)值條件。

    分層流;溫度異重流;尾流場;大渦模擬

    0 引言

    由于太陽的熱輻射及大洋鹽熱環(huán)流等因素,海洋通常是分層的,其密度、溫度和鹽度沿垂直方向有明顯的變化。當潛體在密度分層流體中運動時,流體質(zhì)點會受到體積效應下的浮力作用,產(chǎn)生異于均勻環(huán)境的水動力學現(xiàn)象,如可以明顯觀察到流場中源生內(nèi)波的存在。連續(xù)分層作為海洋中典型的密度分層形式受到學者們格外的關(guān)注。Chomaz,Bonneton和Hopfinger[1]通過對連續(xù)分層流中的球體進行拖曳實驗,研究了近場尾流在不同傅汝德數(shù)下的形態(tài);王進[2]在連續(xù)分層流中對三種不同長徑比拖曳潛體激發(fā)的內(nèi)波進行了實驗,得到了內(nèi)波轉(zhuǎn)捩傅汝德數(shù)經(jīng)驗公式;Long[3]研究了連續(xù)分層流中不同尺度障礙物對定常繞流控制方程求解的影響;Miles[4]在Long模型的基礎(chǔ)上用漸近分析的方法研究了不同傅汝德數(shù)條件下半圓柱繞流形成的lee波;Stevenson[5-6]通過對連續(xù)分層流中的圓柱進行拖曳實驗,先后研究了尾流場中的定常內(nèi)波和流體粘性對lee波的影響,其結(jié)果與Lighthill[7]的彌散波模型相吻合;Boyer[8-9]在連續(xù)分層流水槽中系統(tǒng)地研究了圓柱繞流的流動現(xiàn)象;丁勇[10]基于多相流混合模型研究了線性分層流中的圓柱繞流現(xiàn)象。目前,國內(nèi)外對連續(xù)分層流的研究主要集中在理論和實驗方面,尤其缺少相關(guān)的數(shù)值模擬研究。

    本文通過指定流場中溫度的分布及密度與溫度的函數(shù)關(guān)系實現(xiàn)密度的連續(xù)分層,并假設(shè)數(shù)值過程中通過熱傳導改變的溫度值很小,建立了密度連續(xù)分層流的數(shù)值模擬方法。采用RANS方法和LES方法對其中的圓柱繞流現(xiàn)象進行數(shù)值模擬,探討研究了該方法獲得收斂結(jié)果的數(shù)值條件、收斂結(jié)果的穩(wěn)定性以及尾流場時間歷程的演變,將流場特征與Boyer[8]的實驗值進行對比,證實了此種方法在數(shù)值模擬密度連續(xù)分層流體時的可行性。

    1 數(shù)值模型

    1.1 溫度場假設(shè)

    為了與流體動力學控制方程中的對流、擴散項保持一致,本文中的熱對流是指由于對流引起的溫度物理量輸運,與熱力學中熱對流的意義相異,熱傳導是指由于擴散引起的的溫度物理量輸運。

    流體的密度主要由溫度所決定,為了保證密度垂向的穩(wěn)定分層,需保證溫度在垂向的分布是穩(wěn)定的,而根據(jù)熱力學第二定律,經(jīng)過足夠長時間之后,溫度在垂向的分布必然趨于均勻。因此本文假設(shè)通過熱傳導改變的溫度值在潛體尾跡演化的時間跨度內(nèi)足夠小,即熱傳導的影響可忽略不計,溫度在垂向上的分布短時間內(nèi)是穩(wěn)定不變的,如此便可實現(xiàn)穩(wěn)定的密度分層形式。此時,溫度輸運方程(能量方程)為對流占優(yōu)擴散方程,方程如下所示。

    其中:T為溫度,k為流體的導熱系數(shù),c為流體的比熱容,Γ=k/c為相對于T的廣義擴散系數(shù)。方程左邊為溫度場的局部變化率,右邊兩項分別為熱對流及熱擴散項。熱傳導可以忽略的程度在方程中體現(xiàn)為對流項及擴散項對局部場變化率影響程度的相對大小。

    對流與擴散的強度之比可以用Pe(Peclet數(shù))來衡量,其表達式如下所示:

    其中:u為特征速度,L為特征長度,α=Γ/ρ為特征擴散系數(shù)。Pe數(shù)越大,對流占優(yōu)程度越高,那么必然存在一個臨界Peclet數(shù)Pecr,當Pe>Pecr時,假設(shè)成立,流場中可以形成穩(wěn)定的密度分層。

    1.2 模擬方法

    整個流場的溫度、密度分布及計算模型示意圖如圖1所示。規(guī)定該圖的原點在圓柱圓心位置,x為流向,z為垂向,H為計算域的高度,U為遠方來流速度,g為重力加速度,ρ0和T0分別為流體域中的平均密度和平均溫度。為了與實驗保持一致,計算域的高度設(shè)置為0.2 m,圓柱的直徑設(shè)置為0.024 m,入口及出口邊界足夠遠。

    文中用的無量綱參數(shù):內(nèi)傅汝德數(shù):Fr=U/Nd;雷諾數(shù):Re=Ud/ν;無量綱時間:t′= tU/d,t在數(shù)值模擬條件下為計算時間,無量綱時間用tn′表示,在實驗條件下為繞流時間,無量綱時間用te′表示;N為浮頻率,N=(g△ρ/ρ0H)1/2;ν為流體的運動粘性系數(shù),其它參數(shù)的含義如圖1所示。

    圖1 計算模型示意圖Fig.1 Sketch of the calculation model

    溫度T的垂向分布形式如(3)式所示。當T=300 K時,ρ=1 018.58 kg/m3;當T= 300.2 K時,ρ=998 kg/m3。由此可得整個流域溫度的垂向分布函數(shù)如(4)式所示。

    此時,浮頻率的值為1,與實驗相符。值得注意的是,流體的運動粘性系數(shù)ν雖為流體的固有屬性,但是受溫度的影響很大。根據(jù)第10屆ITTC關(guān)于淡水和鹽水粘性系數(shù)的規(guī)定,在T=300 K附近,如果溫度變化5 K,粘性系數(shù)相差約0.07×10-6m2·s-1,本文流場的最大溫差為0.2 K,粘性系數(shù)在該范圍內(nèi)的變化非常小,因此可以將粘性系數(shù)視作常數(shù)。

    Boyer的實驗表明,當Fr=0.4時,lee波已經(jīng)不是尾流場的主要特征,這時尾流場中出現(xiàn)明顯的渦結(jié)構(gòu),所以規(guī)定當Fr>0.4時為高傅汝德數(shù),F(xiàn)r<0.4時為低傅汝德數(shù)。對低傅汝德數(shù)下的情形,可通過RANS方法予以模擬實現(xiàn),對高傅汝德數(shù)下的情形,則應通過LES方法實現(xiàn)。在高傅汝德數(shù)條件下,尾流場中湍流脈動成份的水動力作用會增強,而RANS方法抹去了瞬時脈動成份,LES方法則沒有這種弊端。理論研究[11-13]認為,當橫向距離大于πd時,利用大渦模擬進行計算,三維的影響就變得非常小,因此本文橫向設(shè)置為0.084 m(3.5d)。

    采用O型網(wǎng)格對流體域進行空間離散,二維及三維情況下的網(wǎng)格數(shù)量分別為30萬和150萬。對于RANS方法,壓力離散格式采用“PRESTO!”,對流項離散格式采用“Second Order Upwind”,擴散項離散格式采用中心差分格式,時間離散格式采用“First Order Upwind”。對于LES方法,壓力離散格式采用“Standard”,對流項離散格式采用“Bounded Central Differencing”,擴散項采用中心差分格式,時間離散格式采用“Second Order Implicit”。同時兩種方法都是采用PISO算法進行求解,時間步長為0.005 s。

    2 數(shù)值結(jié)果與討論

    2.1 低傅汝德數(shù)條件下的數(shù)值結(jié)果

    在低傅汝德數(shù)條件下,計算了三組工況,具體的參數(shù)如表1所示。

    工況1條件下的流線圖如圖2所示,圖中(a)為數(shù)值結(jié)果,(b)為實驗結(jié)果。該圖表明低速條件下雖然在圓柱的下游形成一對與實驗值較為相符的lee波,但是下游的流線形狀并不能與實驗值完全相符。說明在低速條件下對流項的影響較弱,此時熱擴散并不能完全忽略不計。

    工況2條件下結(jié)果穩(wěn)定時的流線圖及相同實驗條件下的流線圖如圖3所示,圖中(a)為數(shù)值結(jié)果,(b)為實驗結(jié)果。(a)圖中可以觀察到在圓柱的后方有兩對較明顯的lee波峰線,其特征與(b)圖的實驗結(jié)果相吻合。圖(a)、(b)在靠近x軸線處lee波的形狀與大小基本一致,不同的是數(shù)值結(jié)果中x軸線附近的流線上下兩側(cè)間距較大,實驗條件下x軸線附近上下兩側(cè)的流線在緊鄰圓柱及圓柱后方第一個波包后端幾近接觸。

    表1 工況表(低傅汝德數(shù))Tab.1 Conditions of calculation (low Froude number)

    圖2 流線圖對比(Fr=0.018,Re=12)Fig.2 Comparison of streamline(Fr=0.018,Re=12)

    圖3 流線圖對比(Fr=0.08,Re=60)Fig.3 Comparison of streamline(Fr=0.08,Re=60)

    工況3條件下的流線圖及實驗條件下的流線圖如圖4所示,圖中(a)為數(shù)值結(jié)果,(b)為實驗結(jié)果。圖(a)表明緊鄰圓柱的尾流場流線呈葫蘆狀,與圖(b)中Boyer實驗得到的流線圖非常吻合,這是該階段尾流場的主要特征。另外,圖(a)中還存在斜向傳播的lee波,這種斜向傳播的lee波在實驗結(jié)果中也可以觀察到,只是lee波的位置略有差異,數(shù)值模擬過程中沒有觀察到轉(zhuǎn)子的出現(xiàn)。在該工況下,數(shù)值結(jié)果與實驗結(jié)果符合較好。

    圖4 流線圖對比(Fr=0.17,Re=98.7)Fig.4 Comparison of streamline(Fr=0.17,Re=98.7)

    2.2 高傅汝德數(shù)條件下的數(shù)值結(jié)果

    在高傅汝德數(shù)條件下,計算了兩組工況,具體的參數(shù)如表2所示。

    圖5為Fr=0.88,Re=480條件下流線圖的實驗結(jié)果及數(shù)值結(jié)果,圖(a)為實驗結(jié)果,圖(b)為數(shù)值結(jié)果。該圖表明數(shù)值和實驗條件下的尾流場流線形狀大致相同,在緊鄰圓柱的后方都有一對渦出現(xiàn),在距離圓柱6倍直徑和16倍直徑的尾流場中存在著兩個波包,第一個波包下包絡(luò)著一對渦,圖(b)中第二個波包所在的位置比圖(a)中第二個波包出現(xiàn)的位置距離圓柱近,這是由于波包是一個不穩(wěn)定的尾流場特征,它在水平及垂向位置上不斷震蕩,震蕩幅度可以達到5倍的圓柱直徑。該數(shù)值方法基本模擬出了實驗工況下的結(jié)果。

    表2 工況表(高傅汝德數(shù))Tab.2 Conditions of calculation (high Froude number)

    圖5 流線圖對比(Fr=0.88,Re=480)Fig.5 Comparison of streamline(Fr=0.88,Re=480)

    圖6 渦量圖對比(Fr=1.77,Re=960)Fig.6 Comparison of vortex(Fr=1.77,Re=960)

    表3 不同導熱系數(shù)下的工況表Tab.3 Conditions of calculation at different heat conductivity

    圖7為Fr=0.018,Re=12時不同Pe數(shù)條件下圓柱繞流流線圖,圖(a)為Pe=7 212的情形,圖(b)為Pe=72 120的情形,兩幅圖中的流線明顯不同于Pe=72.12時的情形(圖2(a))。當Pe=7 212及Pe=72 120時,尾流場流線形狀與實驗結(jié)果更加相符(圖2(b))。圖8為Fr=0.08,Re=60時不同Pe數(shù)條件下圓柱繞流流線圖,圖(a)為Pe=32 054的情形,圖(b)為Pe=320 540的情形,兩幅圖也不同于Pe= 320.54的情形(圖3(a)),當Pe>32 054時,尾流場流線形狀與實驗結(jié)果相符(圖3(b))。圖9為Fr= 0.17,Re=98.7時不同Pe數(shù)條件下圓柱繞流流線圖,圖(a)為Pe=68 114的情形,圖(b)為Pe=681 140的情形,結(jié)果與圖7、圖8的規(guī)律相符。但是當傅汝德數(shù)增大到一定程度時,Pe數(shù)的值已經(jīng)很大了,改變流體導熱系數(shù)對流場的影響變得非常有限。由于連續(xù)分層流中的圓柱繞流是一個準穩(wěn)態(tài)過程,流線隨時間變化,圖7~9各圖中(a)與(b)數(shù)值時刻均相同。

    圖7 流線圖(Fr=0.018,Re=12)Fig.7 Diagram of streamline(Fr=0.018,Re=12)

    圖8 流線圖(Fr=0.08,Re=60)Fig.8 Diagram of streamline(Fr=0.08,Re=60)

    圖9 流線圖(Fr=0.17,Re=98.7)Fig.9 Diagram of streamline(Fr=0.17,Re=98.7)

    分析表明,在同樣導熱系數(shù)條件下,隨著來流速度的增加,Pe數(shù)越來越大,對流占優(yōu)程度越高。雖然低傅汝德數(shù)工況下結(jié)果與實驗值不完全相符,但在高傅汝德數(shù)工況下能夠取得與實驗一致的結(jié)果。因此如果在低速條件下的數(shù)值結(jié)果能夠與實驗值吻合,那么高速條件下該方法理論上也能夠適用。

    2.4 尾流場時間歷程的研究

    在模擬過程中,數(shù)值時間并不能與實驗時間完全相對應,但在某數(shù)值時刻,若上游速度剖面曲線與某一實驗時刻結(jié)果相對應,同時刻下游速度剖面數(shù)值結(jié)果也應與同實驗時刻的結(jié)果相對應。圖10~13為Fr=0.018,Re=12時數(shù)值及實驗條件下上下游(x=±7.5d)速度剖面對比圖,圖10、圖11為Pe= 7 212時的情形,圖12、圖13為Pe=72 120時的情形。速度曲線的極值代表該剖面處速度的最大值或最小值,在上游速度剖面曲線吻合的時刻,下游速度剖面基本吻合,四組結(jié)果中波谷值的相對誤差最大為2%。另外,在Pe=7 212及Pe=72 120兩種情形下,實驗時刻te′=62對應的數(shù)值時刻均是實驗時刻te′=21對應數(shù)值時刻的3倍左右(實驗時刻對應的倍數(shù)也為3),在時間演化上,數(shù)值結(jié)果與實驗結(jié)果相符。

    圖10 上下游速度剖面(Pe=7 212,tn′=7.88,te′=21)Fig.10 Velocity profiles of upstream and downstream(Pe=7 212,tn′=7.88,te′=21)

    圖11 上下游速度剖面(Pe=7 212,tn′=24.06,te′=62)Fig.11 Velocity profiles of upstream and downstream(Pe=7 212,tn′=24.06,te′=62)

    圖12 上下游速度剖面(Pe=72 120,tn′=6.64,te′=21)Fig.12 Velocity profiles of upstream and downstream(Pe=72 120,tn′=6.64,te′=21)

    圖13 上下游速度剖面(Pe=72 120,tn′=18.9,te′=62)Fig.13 Velocity profiles of upstream and downstream(Pe=72 120,tn′=18.9,te′=62)

    3 結(jié)論

    基于溫度異重流模型,通過指定流場中溫度的分布及密度與溫度的函數(shù)關(guān)系實現(xiàn)了密度的連續(xù)分層,并假設(shè)流場中的熱傳導在數(shù)值過程中忽略不計以確保連續(xù)分層的穩(wěn)定性。分別采用RANS方法和LES方法數(shù)值模擬了低傅汝德數(shù)和高傅汝德數(shù)工況下的圓柱繞流,其尾流場特征與實驗相符合,證實了這種方法在數(shù)值模擬密度連續(xù)分層流時的可行性。具體結(jié)論如下:

    (1)隨著傅汝德數(shù)的增大,Peclet數(shù)逐漸變大,能量方程變?yōu)閷α髡純?yōu)擴散方程,流場中的熱傳導可以忽略不計,密度分層能在較長時間內(nèi)保持穩(wěn)定,因此這種方法在數(shù)值模擬連續(xù)分層流時的穩(wěn)定性和精度都比較高。

    (2)采用RANS方法成功模擬出了低傅汝德數(shù)條件下尾流場中的lee波,緊鄰圓柱的尾流場流線呈現(xiàn)葫蘆狀,其特征與實驗值相符。圓柱上下游(x=±7.5d)速度剖面曲線出現(xiàn)了一個極小值和兩個相等的極大值,其變化規(guī)律與實驗定量一致,峰值的誤差在2%以內(nèi),且在時間歷程演變方面與實驗相符合。

    (3)采用LES方法對連續(xù)分層流中高傅汝德數(shù)條件下的圓柱繞流進行了數(shù)值模擬。當Fr=0.88時,距離圓柱6倍直徑和16倍直徑的尾流場中分別出現(xiàn)了一個波包,且其在水平及垂向位置上不斷發(fā)生震蕩。當Fr=1.77時,尾流場進入了完全湍流狀態(tài),圓柱后方出現(xiàn)了很強的渦結(jié)構(gòu),由于分層的抑制作用,渦街在垂向的寬度保持不變。

    (4)RANS方法能夠捕捉到流場中流線的變化,但對渦量的描述不夠細致,LES方法則彌補了這一弊端,能夠精確模擬出流場中的各種成分。因此,RANS方法適合于低傅汝德數(shù)下的流動,LES方法則同時適合低傅汝德數(shù)和高傅汝德數(shù)下的流動。

    [1]Chomaz J M,Bonneton P,Hopfinger E J.The structure of the near wake of a sphere moving horizontally in a stratified fluid[J].Journal of Fluid Mechanics,1993,254:1-21.

    [2]王進,尤云祥,胡天群,等.密度分層流體中不同長徑比拖曳潛體激發(fā)內(nèi)波特性試驗[J].科學通報,2012,08(57): 606-617. Wang Jin,You Yunxiang,Hu Tianqun,et al.The characteristics of internal waves excited by towed bodies with different aspect ratios in a stratified fluid[J].Chin Sci Bull,2012,08(57):606-617.

    [3]Long R R.Some aspects of the flow of stratified fluids[J].Tellus,1955,7(3):341-357.

    [4]Miles J W H.Lee waves in a stratified flow.Part 2.Semi-circular obstacle[J].Journal of Fluid Mechanics,1968,33(4): 803-814.

    [5]Stevenson T N.Some two-dimensional internal waves in a stratified fluid[J].Journal of Fluid Mechanics,1968,33(04): 715-720.

    [6]Stevenson T N,Chang W L,Laws P.Viscous effects in lee waves[J].Geophysical&Astrophysical Fluid Dynamics,1979, 13(1):141-151.

    [7]Lighthill M J.On waves generated in dispersive systems to travelling forcing effects,with applications to the dynamics of rotating fluids[M].Hyperbolic Equations and Waves,Springer,1970:124-152.

    [8]Boyer D L,Davies P A,Fernando H,et al.Linearly stratified flow past a horizontal circular cylinder[J].Philosophical Transactions of the Royal Society of London A:Mathematical,Physical and Engineering Sciences,1989,328(1601):501-528.

    [9]Xu Y,Fernando H J,Boyer D L.Turbulent wakes of stratified flow past a cylinder[J].Physics of Fluids(1994-present), 1995,7(9):2243-2255.

    [10]丁勇,韓盼盼,段菲,馬衛(wèi)狀.線性分層流中圓柱繞流數(shù)值模擬方法研究[J].哈爾濱工程大學學報,2016,9:1-5. Ding Yong,Han Panpan,Duan Fei,Ma Weizhuang.Numerical study of linearly stratified flow past a cylinder based on multiphase mixture model[J].Journal of Harbin Engineering University,2016,9:1-5.

    [11]M B.Large eddy simulation of the subcritical flow past a circular cylinder:numerical and modeling aspects[J].International Journal for Numerical Methods in Fluids,1998,28(9):1281-1302.

    [12]M B.Numerical and modeling influences on large eddy simulations for the flow past a circular cylinder[J].International Journal of Heat And Fluid Flow,1998,19(5):512-521.

    [13]Kravchenko A G M P.Numerical studies of flow over a circular cylinder at ReD=3 900[J].Physics of Fluids,2000,12(2): 403-417.

    Research on the numerical simulation methods of continuously stratified flows based on thermal density current model

    NIU Ming-chang1,DING Yong1,Ma Wei-zhuang1,Han Pan-pan2
    (1.College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China;2.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

    Based on thermal density current model,the continuous stratification of density is achieved by specifying the distribution of temperature and regarding density as a function of temperature.It is assumed that the heat conduction is negligible to ensure the stability of density stratification.The RANS method and LES method are used respectively to simulate the flows past a cylinder at low Froude numbers and high Froude numbers.The lee waves,gourd-shaped streamlines and wave packets are consistent with the experimental results.For the minimum values of velocity profile curves,there are only 2%errors between the calculated and experimental values,and the curve shapes are the same.The influences of Peclet numbers on the stability of this method are discussed by changing the heat conductivity coefficient of fluid.Finally,the numerical conditions of convergence are obtained.

    stratified fluid;thermal density current;wake field;large eddy simulation

    U661.1

    A

    10.3969/j.issn.1007-7294.2017.08.002

    1007-7294(2017)08-0941-09

    2017-03-26

    ××減震降噪工程專項計劃

    牛明昌(1993-),男,碩士研究生;丁勇(1959-),男,博士,教授,通訊作者,E-mail:dingyong@hrbeu.edu.cn。

    猜你喜歡
    尾流層流流線
    層流輥道電機IP56防護等級結(jié)構(gòu)設(shè)計
    防爆電機(2022年5期)2022-11-18 07:40:18
    摻氫對二甲醚層流燃燒特性的影響
    層流切應力誘導microRNA-101下調(diào)EZH2抑制血管新生
    幾何映射
    任意夾角交叉封閉邊界內(nèi)平面流線計算及應用
    飛機尾流的散射特性與探測技術(shù)綜述
    雷達學報(2017年6期)2017-03-26 07:53:06
    錐形流量計尾流流場分析
    水面艦船風尾流效應減弱的模擬研究
    X80鋼層流冷卻溫度場的有限元模擬
    大型綜合交通樞紐流線組織設(shè)計
    黄色女人牲交| 舔av片在线| 国产成人系列免费观看| 欧美三级亚洲精品| 欧美日韩瑟瑟在线播放| 麻豆一二三区av精品| xxx96com| 国产激情久久老熟女| 一卡2卡三卡四卡精品乱码亚洲| 麻豆国产av国片精品| 免费看十八禁软件| av视频在线观看入口| 欧美一区二区国产精品久久精品| 色噜噜av男人的天堂激情| 日本与韩国留学比较| 日韩免费av在线播放| 日韩中文字幕欧美一区二区| 欧美乱码精品一区二区三区| 国产伦精品一区二区三区四那| 99热这里只有是精品50| 丰满人妻熟妇乱又伦精品不卡| 精品99又大又爽又粗少妇毛片 | 亚洲国产高清在线一区二区三| 天堂av国产一区二区熟女人妻| 国产成人啪精品午夜网站| 国产单亲对白刺激| 亚洲精品久久国产高清桃花| 亚洲欧美一区二区三区黑人| 中文字幕熟女人妻在线| 老汉色av国产亚洲站长工具| 黄色女人牲交| 久久精品国产亚洲av香蕉五月| 18禁黄网站禁片免费观看直播| 亚洲欧美日韩无卡精品| 大型黄色视频在线免费观看| 亚洲一区二区三区色噜噜| 国产精品1区2区在线观看.| 国内精品一区二区在线观看| 久久热在线av| 亚洲欧洲精品一区二区精品久久久| 午夜成年电影在线免费观看| 日本三级黄在线观看| 欧美一级毛片孕妇| 久久精品国产99精品国产亚洲性色| 中文字幕精品亚洲无线码一区| 禁无遮挡网站| 桃红色精品国产亚洲av| 午夜免费激情av| 日本一二三区视频观看| 99久久无色码亚洲精品果冻| 欧美zozozo另类| 两个人视频免费观看高清| 99国产综合亚洲精品| 亚洲性夜色夜夜综合| 亚洲精品美女久久久久99蜜臀| 18禁国产床啪视频网站| 久久精品国产清高在天天线| 午夜a级毛片| 亚洲激情在线av| 国产私拍福利视频在线观看| 国产伦精品一区二区三区视频9 | 国产成人精品久久二区二区免费| 国产精品一区二区免费欧美| 丰满人妻熟妇乱又伦精品不卡| 亚洲精品乱码久久久v下载方式 | 97超级碰碰碰精品色视频在线观看| 久9热在线精品视频| 国模一区二区三区四区视频 | 亚洲欧美日韩东京热| 少妇熟女aⅴ在线视频| 日韩国内少妇激情av| 一个人看视频在线观看www免费 | 久久久久久久久免费视频了| 听说在线观看完整版免费高清| 麻豆久久精品国产亚洲av| 欧美日韩综合久久久久久 | 99久久成人亚洲精品观看| 精品欧美国产一区二区三| 亚洲精华国产精华精| 日韩高清综合在线| 免费观看的影片在线观看| 一本久久中文字幕| 毛片女人毛片| 精品熟女少妇八av免费久了| 熟女人妻精品中文字幕| 非洲黑人性xxxx精品又粗又长| 免费在线观看成人毛片| 欧美中文日本在线观看视频| 在线看三级毛片| 中文字幕精品亚洲无线码一区| 亚洲欧美激情综合另类| 久久精品国产99精品国产亚洲性色| 特大巨黑吊av在线直播| 精品久久久久久成人av| 亚洲天堂国产精品一区在线| 高清在线国产一区| 伦理电影免费视频| 欧美国产日韩亚洲一区| www.熟女人妻精品国产| 国产 一区 欧美 日韩| 欧美大码av| 亚洲 欧美一区二区三区| 亚洲av成人一区二区三| 欧美日韩瑟瑟在线播放| 亚洲专区国产一区二区| 国产亚洲精品久久久久久毛片| 国产亚洲精品久久久久久毛片| 亚洲中文字幕一区二区三区有码在线看 | 在线国产一区二区在线| 欧美国产日韩亚洲一区| 日本a在线网址| 国产成人精品久久二区二区91| 日本a在线网址| 久久久国产成人免费| 别揉我奶头~嗯~啊~动态视频| 欧美色欧美亚洲另类二区| 又爽又黄无遮挡网站| 精品99又大又爽又粗少妇毛片 | 亚洲欧美日韩高清专用| 亚洲国产精品成人综合色| 啪啪无遮挡十八禁网站| 18美女黄网站色大片免费观看| 久久午夜亚洲精品久久| 久久精品91无色码中文字幕| 免费观看的影片在线观看| 麻豆成人午夜福利视频| 99久久精品国产亚洲精品| 狂野欧美激情性xxxx| 久久久国产成人精品二区| 狂野欧美激情性xxxx| 国产人伦9x9x在线观看| 国产1区2区3区精品| 国产成人精品久久二区二区免费| www日本在线高清视频| 黄色成人免费大全| 亚洲va日本ⅴa欧美va伊人久久| 国产伦一二天堂av在线观看| 天天一区二区日本电影三级| 黑人操中国人逼视频| 成人国产综合亚洲| 国产黄片美女视频| 两个人的视频大全免费| 在线视频色国产色| 色在线成人网| 欧美午夜高清在线| 天堂动漫精品| 亚洲无线在线观看| 亚洲性夜色夜夜综合| 精品国产乱码久久久久久男人| 免费大片18禁| 欧美zozozo另类| 亚洲五月婷婷丁香| 夜夜爽天天搞| 午夜免费激情av| bbb黄色大片| 女人被狂操c到高潮| 久久午夜综合久久蜜桃| 可以在线观看毛片的网站| 国产伦精品一区二区三区四那| 99国产精品99久久久久| 亚洲国产精品合色在线| 夜夜看夜夜爽夜夜摸| 又粗又爽又猛毛片免费看| 男人舔奶头视频| 国产午夜福利久久久久久| 两个人的视频大全免费| 国产av在哪里看| 亚洲欧美日韩东京热| 51午夜福利影视在线观看| 午夜亚洲福利在线播放| 噜噜噜噜噜久久久久久91| 国产高清有码在线观看视频| 男插女下体视频免费在线播放| 极品教师在线免费播放| 在线a可以看的网站| 国产精品影院久久| 国产乱人视频| 国产爱豆传媒在线观看| 白带黄色成豆腐渣| 亚洲人成网站在线播放欧美日韩| 搡老岳熟女国产| 久久精品aⅴ一区二区三区四区| 亚洲精品国产精品久久久不卡| 国产黄色小视频在线观看| 最近最新中文字幕大全免费视频| 国产精品 国内视频| 欧美日韩福利视频一区二区| 欧美极品一区二区三区四区| 日日摸夜夜添夜夜添小说| 久久精品亚洲精品国产色婷小说| 一二三四在线观看免费中文在| 日韩人妻高清精品专区| 最新中文字幕久久久久 | 午夜福利免费观看在线| 岛国在线免费视频观看| 午夜精品久久久久久毛片777| 国内精品久久久久久久电影| 国产精品永久免费网站| 韩国av一区二区三区四区| 国产精品亚洲一级av第二区| 亚洲乱码一区二区免费版| 丰满的人妻完整版| 国产精品一区二区三区四区久久| 观看美女的网站| 国产亚洲精品久久久com| 免费在线观看成人毛片| 欧洲精品卡2卡3卡4卡5卡区| 一级毛片精品| 国产精品98久久久久久宅男小说| 人人妻人人澡欧美一区二区| 午夜免费观看网址| 国产av不卡久久| 午夜福利在线观看吧| xxxwww97欧美| 国产在线精品亚洲第一网站| 搡老熟女国产l中国老女人| 国语自产精品视频在线第100页| 久久久久久人人人人人| 国产主播在线观看一区二区| or卡值多少钱| 性色avwww在线观看| 久久草成人影院| 国内精品久久久久精免费| 超碰成人久久| 黄色 视频免费看| av国产免费在线观看| 日韩三级视频一区二区三区| 久久精品91无色码中文字幕| 亚洲人成伊人成综合网2020| 久久久精品大字幕| 日本熟妇午夜| 日韩欧美国产一区二区入口| 黄片小视频在线播放| 亚洲aⅴ乱码一区二区在线播放| 日韩有码中文字幕| 国产成人欧美在线观看| 国产精品电影一区二区三区| 在线观看66精品国产| 18禁黄网站禁片免费观看直播| 最近视频中文字幕2019在线8| 1024手机看黄色片| 成人18禁在线播放| 91九色精品人成在线观看| 黄色 视频免费看| 精品国产美女av久久久久小说| 欧美成人性av电影在线观看| 99久久精品国产亚洲精品| 国产探花在线观看一区二区| 曰老女人黄片| 美女免费视频网站| 淫秽高清视频在线观看| 亚洲成av人片免费观看| 国产精品久久电影中文字幕| 国产免费男女视频| 一区二区三区高清视频在线| 99热只有精品国产| 99久久无色码亚洲精品果冻| 成年女人毛片免费观看观看9| 中文字幕高清在线视频| 九九在线视频观看精品| 少妇丰满av| 日韩欧美免费精品| 午夜激情欧美在线| 成人精品一区二区免费| ponron亚洲| 国产黄a三级三级三级人| 嫁个100分男人电影在线观看| 成人亚洲精品av一区二区| 日韩欧美在线乱码| 老熟妇仑乱视频hdxx| 国产精品免费一区二区三区在线| 欧美激情在线99| e午夜精品久久久久久久| 美女被艹到高潮喷水动态| 国产午夜福利久久久久久| 1024香蕉在线观看| 色精品久久人妻99蜜桃| 一卡2卡三卡四卡精品乱码亚洲| 国产成人精品久久二区二区免费| 国产极品精品免费视频能看的| 国产av在哪里看| 免费av毛片视频| 日本免费a在线| x7x7x7水蜜桃| 久久久精品大字幕| 一区二区三区激情视频| 丰满的人妻完整版| 免费观看人在逋| 可以在线观看的亚洲视频| 高潮久久久久久久久久久不卡| 99视频精品全部免费 在线 | 三级毛片av免费| 国产成人av激情在线播放| 久久国产精品影院| 精品久久久久久,| 亚洲av电影在线进入| 亚洲人成电影免费在线| 丁香六月欧美| 亚洲成人精品中文字幕电影| 国产成人福利小说| 1024手机看黄色片| 亚洲在线自拍视频| 亚洲成人免费电影在线观看| 男女做爰动态图高潮gif福利片| 香蕉丝袜av| 日本撒尿小便嘘嘘汇集6| 一级毛片精品| 身体一侧抽搐| 哪里可以看免费的av片| 久久性视频一级片| 亚洲在线自拍视频| 99久久久亚洲精品蜜臀av| 免费高清视频大片| av女优亚洲男人天堂 | 亚洲国产精品999在线| 日本a在线网址| 少妇的丰满在线观看| 国内精品一区二区在线观看| 啦啦啦免费观看视频1| cao死你这个sao货| 成熟少妇高潮喷水视频| 色老头精品视频在线观看| 男女午夜视频在线观看| 一个人免费在线观看的高清视频| 最近在线观看免费完整版| 亚洲av成人精品一区久久| 麻豆久久精品国产亚洲av| 国产精品一区二区三区四区久久| 欧美成狂野欧美在线观看| 真实男女啪啪啪动态图| 国产亚洲av高清不卡| 国产私拍福利视频在线观看| 亚洲av电影在线进入| 国产成人影院久久av| 国产在线精品亚洲第一网站| 精品国产亚洲在线| 免费av毛片视频| 亚洲五月婷婷丁香| 亚洲av成人精品一区久久| 国产视频一区二区在线看| 五月玫瑰六月丁香| 毛片女人毛片| 真实男女啪啪啪动态图| 国产一区二区三区视频了| 国产成年人精品一区二区| 青草久久国产| 日韩大尺度精品在线看网址| 两个人的视频大全免费| 国产三级黄色录像| 亚洲色图 男人天堂 中文字幕| 看片在线看免费视频| 午夜福利18| 99国产精品一区二区蜜桃av| 成人18禁在线播放| 夜夜看夜夜爽夜夜摸| 长腿黑丝高跟| 国产成人av教育| 亚洲午夜理论影院| 亚洲国产精品999在线| 亚洲欧美一区二区三区黑人| 国产激情久久老熟女| 手机成人av网站| 国产精品久久久人人做人人爽| 国产亚洲精品久久久com| 夜夜躁狠狠躁天天躁| 午夜福利高清视频| 国产三级黄色录像| 欧美中文日本在线观看视频| 欧美日韩精品网址| 午夜免费观看网址| 麻豆成人午夜福利视频| 9191精品国产免费久久| 啦啦啦免费观看视频1| 夜夜爽天天搞| 欧美大码av| 国产高清视频在线播放一区| 嫩草影院精品99| 99久国产av精品| 欧美三级亚洲精品| 九九热线精品视视频播放| 亚洲一区高清亚洲精品| 国产精品1区2区在线观看.| 日本精品一区二区三区蜜桃| 18禁黄网站禁片午夜丰满| 无限看片的www在线观看| 日韩 欧美 亚洲 中文字幕| 十八禁人妻一区二区| 国产探花在线观看一区二区| 我的老师免费观看完整版| 免费大片18禁| 三级毛片av免费| 免费在线观看视频国产中文字幕亚洲| 女人高潮潮喷娇喘18禁视频| 狠狠狠狠99中文字幕| 最新美女视频免费是黄的| 女同久久另类99精品国产91| 中文字幕人妻丝袜一区二区| 亚洲成人久久爱视频| 一本综合久久免费| 精品国产三级普通话版| 十八禁人妻一区二区| 一级毛片精品| 成人三级黄色视频| 国模一区二区三区四区视频 | www日本黄色视频网| 十八禁人妻一区二区| 中文字幕av在线有码专区| 亚洲 国产 在线| 午夜免费观看网址| 成年人黄色毛片网站| 一级毛片精品| 精品久久久久久,| 老汉色av国产亚洲站长工具| 免费av不卡在线播放| 久久久久久大精品| 日本撒尿小便嘘嘘汇集6| 国产精品99久久99久久久不卡| 哪里可以看免费的av片| 亚洲一区高清亚洲精品| 国内精品一区二区在线观看| 国产精品久久久人人做人人爽| 最近在线观看免费完整版| 9191精品国产免费久久| 亚洲片人在线观看| 男人和女人高潮做爰伦理| 99久久99久久久精品蜜桃| 色噜噜av男人的天堂激情| 欧美色视频一区免费| 亚洲va日本ⅴa欧美va伊人久久| 日韩欧美一区二区三区在线观看| 亚洲成人久久爱视频| 日本熟妇午夜| 男女下面进入的视频免费午夜| 桃红色精品国产亚洲av| 99精品在免费线老司机午夜| 18禁黄网站禁片免费观看直播| 在线观看日韩欧美| 欧美激情久久久久久爽电影| 1024香蕉在线观看| 国内精品久久久久久久电影| 国产成+人综合+亚洲专区| 亚洲激情在线av| 在线观看日韩欧美| 国产精品国产高清国产av| 露出奶头的视频| 特大巨黑吊av在线直播| 欧美日韩黄片免| 亚洲欧美一区二区三区黑人| 丁香六月欧美| 久久精品国产亚洲av香蕉五月| 亚洲av免费在线观看| 12—13女人毛片做爰片一| 亚洲欧美一区二区三区黑人| 九色国产91popny在线| 国产精品久久久av美女十八| 国产淫片久久久久久久久 | 国产不卡一卡二| 最近在线观看免费完整版| 亚洲色图av天堂| 九九热线精品视视频播放| 国产 一区 欧美 日韩| 午夜福利18| 在线观看日韩欧美| 一级黄色大片毛片| 黄色视频,在线免费观看| 观看免费一级毛片| 一级a爱片免费观看的视频| 后天国语完整版免费观看| 亚洲欧美日韩东京热| 19禁男女啪啪无遮挡网站| 亚洲第一电影网av| 99国产极品粉嫩在线观看| 日本一二三区视频观看| 在线十欧美十亚洲十日本专区| 久久精品影院6| 在线播放国产精品三级| 最好的美女福利视频网| 一本一本综合久久| 手机成人av网站| 欧美一区二区国产精品久久精品| 中文字幕最新亚洲高清| 亚洲成人中文字幕在线播放| 18禁黄网站禁片午夜丰满| 一二三四在线观看免费中文在| 91av网站免费观看| 人妻丰满熟妇av一区二区三区| 亚洲国产精品sss在线观看| 真人一进一出gif抽搐免费| 欧美激情久久久久久爽电影| 哪里可以看免费的av片| 国产成人啪精品午夜网站| 伊人久久大香线蕉亚洲五| 性色avwww在线观看| 搞女人的毛片| 国产av不卡久久| 99热这里只有精品一区 | 亚洲中文av在线| 久久伊人香网站| 亚洲乱码一区二区免费版| 午夜精品在线福利| 亚洲最大成人中文| 国产精品九九99| 中文字幕最新亚洲高清| 午夜影院日韩av| 无人区码免费观看不卡| 精品不卡国产一区二区三区| 我要搜黄色片| 国产激情欧美一区二区| 九九久久精品国产亚洲av麻豆 | 国产伦精品一区二区三区四那| 国产精品久久久人人做人人爽| 日本黄色片子视频| 老汉色∧v一级毛片| 国产精品久久视频播放| 亚洲欧美日韩东京热| 久久香蕉国产精品| 国产又黄又爽又无遮挡在线| 中文字幕人成人乱码亚洲影| 国产成人av教育| 亚洲熟女毛片儿| 国产精品 国内视频| 欧美乱色亚洲激情| 国产亚洲欧美在线一区二区| 国产又黄又爽又无遮挡在线| 久久精品aⅴ一区二区三区四区| 亚洲欧美精品综合久久99| 99久久久亚洲精品蜜臀av| 国模一区二区三区四区视频 | 欧美性猛交黑人性爽| 国产激情欧美一区二区| 国产91精品成人一区二区三区| 欧美另类亚洲清纯唯美| 女人高潮潮喷娇喘18禁视频| 在线观看免费午夜福利视频| 久久精品人妻少妇| 国产成人av教育| 午夜福利在线在线| 国内精品久久久久久久电影| 国产av在哪里看| av片东京热男人的天堂| 国产精品综合久久久久久久免费| 久久性视频一级片| 在线观看免费视频日本深夜| 色播亚洲综合网| 国产伦在线观看视频一区| 亚洲av日韩精品久久久久久密| 欧美最黄视频在线播放免费| 亚洲,欧美精品.| www.熟女人妻精品国产| 欧美日本亚洲视频在线播放| 久久这里只有精品19| 亚洲 国产 在线| 亚洲国产色片| 国产午夜福利久久久久久| 最近最新免费中文字幕在线| 亚洲最大成人中文| 天天添夜夜摸| netflix在线观看网站| 国产午夜精品久久久久久| 国产高清三级在线| 国产精品久久久久久精品电影| 亚洲欧美精品综合久久99| 久久久国产欧美日韩av| 国产精品美女特级片免费视频播放器 | 老司机午夜福利在线观看视频| 日韩免费av在线播放| 亚洲成人久久爱视频| 久久久久久人人人人人| 最近在线观看免费完整版| 身体一侧抽搐| 最近在线观看免费完整版| 999精品在线视频| 黄片大片在线免费观看| 亚洲精品一卡2卡三卡4卡5卡| 欧美激情在线99| 国产精品99久久99久久久不卡| 国模一区二区三区四区视频 | 亚洲性夜色夜夜综合| 国产成人精品久久二区二区免费| 国产一区二区在线av高清观看| 国内精品久久久久精免费| 欧美黑人欧美精品刺激| www.精华液| av国产免费在线观看| 亚洲第一欧美日韩一区二区三区| 亚洲自拍偷在线| 两人在一起打扑克的视频| 人人妻人人澡欧美一区二区| 在线免费观看不下载黄p国产 | 国产男靠女视频免费网站| 午夜激情欧美在线| 999久久久精品免费观看国产| 午夜福利视频1000在线观看| 一二三四在线观看免费中文在| 中国美女看黄片| 午夜两性在线视频| 午夜福利视频1000在线观看| 一区二区三区激情视频| 给我免费播放毛片高清在线观看| 久久香蕉精品热| 青草久久国产| 999久久久国产精品视频| 在线免费观看不下载黄p国产 | 97碰自拍视频| 国产真实乱freesex| 男人舔女人的私密视频| 黄色片一级片一级黄色片| 婷婷精品国产亚洲av在线| 国产亚洲欧美98| 一级a爱片免费观看的视频| 女人被狂操c到高潮| 丁香欧美五月| 精品久久久久久成人av| 久久精品国产亚洲av香蕉五月| 国产精品永久免费网站|