• <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è)計
    亚洲va日本ⅴa欧美va伊人久久| av在线播放免费不卡| 黄色女人牲交| 最近最新中文字幕大全免费视频| cao死你这个sao货| 久久天躁狠狠躁夜夜2o2o| cao死你这个sao货| a在线观看视频网站| 久久久国产成人免费| 国产精品秋霞免费鲁丝片| 如日韩欧美国产精品一区二区三区| 国产精品久久久久久人妻精品电影| 人成视频在线观看免费观看| 如日韩欧美国产精品一区二区三区| 午夜成年电影在线免费观看| 亚洲最大成人中文| 久久狼人影院| 中亚洲国语对白在线视频| 亚洲一区高清亚洲精品| 国产伦一二天堂av在线观看| 亚洲 欧美 日韩 在线 免费| 一本综合久久免费| 夜夜爽天天搞| 国产精品日韩av在线免费观看 | 国产精品日韩av在线免费观看 | 中文字幕高清在线视频| 岛国在线观看网站| 精品国产一区二区三区四区第35| 高清黄色对白视频在线免费看| 国产亚洲欧美98| 亚洲av片天天在线观看| 女人高潮潮喷娇喘18禁视频| 亚洲精品中文字幕一二三四区| 性少妇av在线| 欧美精品啪啪一区二区三区| 婷婷六月久久综合丁香| 禁无遮挡网站| 脱女人内裤的视频| 涩涩av久久男人的天堂| а√天堂www在线а√下载| 十分钟在线观看高清视频www| 真人做人爱边吃奶动态| 亚洲自拍偷在线| 国产成人精品在线电影| 亚洲成国产人片在线观看| 淫妇啪啪啪对白视频| 国产精品 欧美亚洲| 黄色丝袜av网址大全| 看黄色毛片网站| 免费无遮挡裸体视频| 国产野战对白在线观看| 欧美日本中文国产一区发布| 制服诱惑二区| 国产激情久久老熟女| 午夜福利影视在线免费观看| 久久久久久人人人人人| 亚洲无线在线观看| 中亚洲国语对白在线视频| 婷婷丁香在线五月| 欧美最黄视频在线播放免费| 一进一出抽搐gif免费好疼| 亚洲第一av免费看| 成人亚洲精品一区在线观看| 一个人观看的视频www高清免费观看 | 99久久99久久久精品蜜桃| 久久热在线av| 在线观看免费日韩欧美大片| 淫秽高清视频在线观看| 丁香欧美五月| 亚洲美女黄片视频| 日本精品一区二区三区蜜桃| 国产精品日韩av在线免费观看 | 欧美成人性av电影在线观看| 搡老妇女老女人老熟妇| 国产一区二区激情短视频| 在线观看一区二区三区| 美女高潮喷水抽搐中文字幕| 男人的好看免费观看在线视频 | tocl精华| 9色porny在线观看| 精品人妻1区二区| 一二三四社区在线视频社区8| 中亚洲国语对白在线视频| 一个人观看的视频www高清免费观看 | 好男人电影高清在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 精品无人区乱码1区二区| 99久久国产精品久久久| 桃色一区二区三区在线观看| 国产亚洲欧美98| 亚洲国产精品久久男人天堂| 国产亚洲av高清不卡| 久久国产精品男人的天堂亚洲| 性色av乱码一区二区三区2| 国产成人av教育| 人人澡人人妻人| 777久久人妻少妇嫩草av网站| 色综合亚洲欧美另类图片| 免费女性裸体啪啪无遮挡网站| 一进一出抽搐动态| 久久久久国产精品人妻aⅴ院| 日韩欧美国产一区二区入口| 国产一区二区在线av高清观看| 非洲黑人性xxxx精品又粗又长| 脱女人内裤的视频| 国产高清视频在线播放一区| or卡值多少钱| 欧美丝袜亚洲另类 | 欧美+亚洲+日韩+国产| 成人亚洲精品一区在线观看| 亚洲成a人片在线一区二区| 一级毛片精品| 少妇 在线观看| 亚洲av美国av| 后天国语完整版免费观看| 免费一级毛片在线播放高清视频 | 嫁个100分男人电影在线观看| 九色国产91popny在线| 欧美国产精品va在线观看不卡| 国产精品久久久av美女十八| 狠狠狠狠99中文字幕| 久久草成人影院| 精品卡一卡二卡四卡免费| 99久久久亚洲精品蜜臀av| 视频区欧美日本亚洲| 男人操女人黄网站| 免费一级毛片在线播放高清视频 | 亚洲精华国产精华精| 欧美国产日韩亚洲一区| 国语自产精品视频在线第100页| 久久国产精品男人的天堂亚洲| 久久精品91无色码中文字幕| 制服人妻中文乱码| 好男人电影高清在线观看| 黄色视频不卡| 中文字幕色久视频| 一本大道久久a久久精品| 一区福利在线观看| 久久久久久久久久久久大奶| 久久 成人 亚洲| 老司机靠b影院| 天天添夜夜摸| 精品国产一区二区久久| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲国产精品久久男人天堂| 久久婷婷成人综合色麻豆| 午夜福利,免费看| 国产欧美日韩综合在线一区二区| 久久久久久久久中文| 成人三级黄色视频| 久久人妻熟女aⅴ| 精品欧美一区二区三区在线| 精品高清国产在线一区| 美国免费a级毛片| 熟妇人妻久久中文字幕3abv| 久久久久久免费高清国产稀缺| 国产视频一区二区在线看| 不卡一级毛片| 免费看a级黄色片| 午夜福利18| 少妇粗大呻吟视频| 久久久国产欧美日韩av| 久久天堂一区二区三区四区| 国产精品自产拍在线观看55亚洲| 亚洲精品久久国产高清桃花| 亚洲国产欧美一区二区综合| 女人被狂操c到高潮| av电影中文网址| 少妇粗大呻吟视频| 宅男免费午夜| 国产精品久久久人人做人人爽| 国产一区二区激情短视频| av超薄肉色丝袜交足视频| 禁无遮挡网站| 精品日产1卡2卡| 高清在线国产一区| 高清毛片免费观看视频网站| 美女高潮喷水抽搐中文字幕| 欧美日韩福利视频一区二区| 91九色精品人成在线观看| 国产av在哪里看| 久久人妻福利社区极品人妻图片| 久久影院123| 黑人巨大精品欧美一区二区蜜桃| 亚洲av第一区精品v没综合| 中文字幕最新亚洲高清| 亚洲精品国产区一区二| 亚洲性夜色夜夜综合| 国产97色在线日韩免费| 久久国产精品影院| 亚洲熟妇中文字幕五十中出| 欧美乱妇无乱码| 国产在线精品亚洲第一网站| 999久久久国产精品视频| 国产成人精品久久二区二区91| 亚洲男人的天堂狠狠| 国产亚洲精品久久久久久毛片| 日韩欧美三级三区| 国产欧美日韩综合在线一区二区| 亚洲欧美日韩另类电影网站| 国产三级黄色录像| 久久久久精品国产欧美久久久| 久久人人97超碰香蕉20202| 午夜福利免费观看在线| 欧美中文综合在线视频| 亚洲精品国产区一区二| a级毛片在线看网站| 午夜成年电影在线免费观看| 午夜福利18| 国产成人精品无人区| 欧美一级a爱片免费观看看 | 熟妇人妻久久中文字幕3abv| 午夜免费观看网址| 久久精品国产99精品国产亚洲性色 | 波多野结衣高清无吗| 国产精品爽爽va在线观看网站 | 午夜久久久在线观看| 欧美日韩黄片免| 日本精品一区二区三区蜜桃| 国产亚洲欧美精品永久| 精品国产亚洲在线| 久久精品人人爽人人爽视色| 色av中文字幕| 妹子高潮喷水视频| 制服丝袜大香蕉在线| 欧美中文综合在线视频| 纯流量卡能插随身wifi吗| 国产成人精品无人区| 多毛熟女@视频| 国产亚洲精品久久久久久毛片| 国产精品精品国产色婷婷| 婷婷丁香在线五月| 国产99白浆流出| 制服丝袜大香蕉在线| 午夜久久久在线观看| 亚洲精品在线美女| 99在线人妻在线中文字幕| 精品无人区乱码1区二区| 中文字幕精品免费在线观看视频| 天天躁狠狠躁夜夜躁狠狠躁| av在线天堂中文字幕| 欧洲精品卡2卡3卡4卡5卡区| 国产午夜福利久久久久久| 日日干狠狠操夜夜爽| 久久性视频一级片| 亚洲人成伊人成综合网2020| 亚洲精品国产一区二区精华液| 日韩大尺度精品在线看网址 | 99re在线观看精品视频| 在线免费观看的www视频| 亚洲五月婷婷丁香| 亚洲一卡2卡3卡4卡5卡精品中文| 欧美一级毛片孕妇| 咕卡用的链子| 亚洲一区二区三区色噜噜| 男女午夜视频在线观看| 久久中文字幕人妻熟女| 亚洲一卡2卡3卡4卡5卡精品中文| 午夜免费鲁丝| 一边摸一边做爽爽视频免费| 欧美在线黄色| 国产一区二区三区综合在线观看| 最近最新中文字幕大全免费视频| 99香蕉大伊视频| 亚洲欧洲精品一区二区精品久久久| 国产亚洲精品av在线| 美国免费a级毛片| 精品人妻1区二区| 黄频高清免费视频| 夜夜躁狠狠躁天天躁| 韩国av一区二区三区四区| 国产精品久久久av美女十八| 这个男人来自地球电影免费观看| av免费在线观看网站| 免费看美女性在线毛片视频| 精品一区二区三区av网在线观看| 免费在线观看影片大全网站| 午夜精品在线福利| 免费久久久久久久精品成人欧美视频| 久久 成人 亚洲| 黄色毛片三级朝国网站| 女性被躁到高潮视频| 首页视频小说图片口味搜索| 欧美一级毛片孕妇| 黑人巨大精品欧美一区二区mp4| 精品国产乱码久久久久久男人| 午夜成年电影在线免费观看| 亚洲av成人不卡在线观看播放网| 91国产中文字幕| 日日爽夜夜爽网站| 国产成人系列免费观看| 女同久久另类99精品国产91| or卡值多少钱| 日韩大尺度精品在线看网址 | 在线免费观看的www视频| 男人操女人黄网站| 国产成人免费无遮挡视频| 美女午夜性视频免费| 一边摸一边抽搐一进一小说| 一区二区三区国产精品乱码| av视频免费观看在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 咕卡用的链子| 亚洲,欧美精品.| 成人亚洲精品av一区二区| 757午夜福利合集在线观看| 精品一品国产午夜福利视频| 亚洲一卡2卡3卡4卡5卡精品中文| 在线观看舔阴道视频| 精品人妻在线不人妻| 精品国产一区二区久久| 亚洲av五月六月丁香网| 美女高潮到喷水免费观看| 成人国语在线视频| 亚洲五月色婷婷综合| 日韩av在线大香蕉| 免费观看精品视频网站| 国产在线观看jvid| 99精品久久久久人妻精品| 国产精品电影一区二区三区| 亚洲,欧美精品.| 性少妇av在线| 久久香蕉精品热| АⅤ资源中文在线天堂| 成熟少妇高潮喷水视频| 免费少妇av软件| 国产精华一区二区三区| 国产免费男女视频| 亚洲成av片中文字幕在线观看| 麻豆国产av国片精品| 在线观看www视频免费| videosex国产| 午夜久久久在线观看| 身体一侧抽搐| 日本免费a在线| 亚洲av成人一区二区三| 一区二区日韩欧美中文字幕| 黄色片一级片一级黄色片| 精品久久久精品久久久| 国产精品二区激情视频| 国产国语露脸激情在线看| 亚洲国产日韩欧美精品在线观看 | 国产精品99久久99久久久不卡| 欧美老熟妇乱子伦牲交| 久久久久久亚洲精品国产蜜桃av| 精品不卡国产一区二区三区| 91成年电影在线观看| 高清在线国产一区| 人人妻人人澡欧美一区二区 | 国产精品98久久久久久宅男小说| www.999成人在线观看| 欧美色欧美亚洲另类二区 | 亚洲自拍偷在线| 亚洲激情在线av| 大码成人一级视频| 欧美亚洲日本最大视频资源| 男女下面插进去视频免费观看| 亚洲专区字幕在线| 国产私拍福利视频在线观看| 欧美国产日韩亚洲一区| 99国产极品粉嫩在线观看| 啪啪无遮挡十八禁网站| 黄片播放在线免费| 88av欧美| 美女扒开内裤让男人捅视频| 亚洲av成人不卡在线观看播放网| 国产亚洲欧美98| 亚洲熟女毛片儿| 亚洲九九香蕉| 中文字幕久久专区| 精品一品国产午夜福利视频| 日本在线视频免费播放| 国产亚洲精品第一综合不卡| 精品熟女少妇八av免费久了| 午夜激情av网站| 精品人妻在线不人妻| 桃色一区二区三区在线观看| 不卡一级毛片| 亚洲aⅴ乱码一区二区在线播放 | 99久久综合精品五月天人人| 午夜精品国产一区二区电影| 一级,二级,三级黄色视频| 又紧又爽又黄一区二区| 制服丝袜大香蕉在线| 人人澡人人妻人| 欧美av亚洲av综合av国产av| 女人被狂操c到高潮| 午夜视频精品福利| 18美女黄网站色大片免费观看| 怎么达到女性高潮| 久久久久久大精品| 正在播放国产对白刺激| 亚洲国产毛片av蜜桃av| 国产精华一区二区三区| 纯流量卡能插随身wifi吗| 电影成人av| 久久久久久久久中文| 欧美激情 高清一区二区三区| 色尼玛亚洲综合影院| 久久精品国产综合久久久| 在线永久观看黄色视频| 97超级碰碰碰精品色视频在线观看| 久久国产精品影院| 国产成人精品无人区| 9色porny在线观看| 黄色视频不卡| www.熟女人妻精品国产| 色在线成人网| 国产精品美女特级片免费视频播放器 | 超碰成人久久| videosex国产| 欧美久久黑人一区二区| 老司机靠b影院| 欧美日韩一级在线毛片| 国产亚洲精品综合一区在线观看 | 不卡av一区二区三区| 妹子高潮喷水视频| 无遮挡黄片免费观看| 欧美黑人欧美精品刺激| 三级毛片av免费| 午夜两性在线视频| 久久中文字幕一级| 啦啦啦韩国在线观看视频| 亚洲片人在线观看| 欧美乱妇无乱码| 色婷婷久久久亚洲欧美| 夜夜看夜夜爽夜夜摸| 热re99久久国产66热| 少妇 在线观看| 好男人电影高清在线观看| 亚洲性夜色夜夜综合| √禁漫天堂资源中文www| 欧美日韩黄片免| 午夜视频精品福利| 制服人妻中文乱码| 亚洲国产欧美网| 巨乳人妻的诱惑在线观看| 国产精品免费一区二区三区在线| 国产一卡二卡三卡精品| 国产单亲对白刺激| 亚洲成人免费电影在线观看| 精品乱码久久久久久99久播| 人妻久久中文字幕网| 不卡一级毛片| 亚洲精品久久国产高清桃花| 欧美性长视频在线观看| 极品教师在线免费播放| 18禁黄网站禁片午夜丰满| 少妇熟女aⅴ在线视频| 天堂√8在线中文| 熟女少妇亚洲综合色aaa.| 男女午夜视频在线观看| 国产亚洲av嫩草精品影院| 999精品在线视频| 国产高清视频在线播放一区| 淫秽高清视频在线观看| 久久香蕉激情| 国产精品一区二区免费欧美| 亚洲中文字幕一区二区三区有码在线看 | 美女高潮到喷水免费观看| 国产精品久久久久久人妻精品电影| 精品无人区乱码1区二区| 一区二区三区国产精品乱码| 在线观看66精品国产| 免费av毛片视频| av网站免费在线观看视频| 亚洲欧美日韩高清在线视频| 夜夜看夜夜爽夜夜摸| 亚洲无线在线观看| 老司机午夜福利在线观看视频| 亚洲电影在线观看av| 国产午夜精品久久久久久| 九色国产91popny在线| 欧美黑人精品巨大| 欧美另类亚洲清纯唯美| 精品久久久久久,| 操出白浆在线播放| 曰老女人黄片| 大陆偷拍与自拍| 日韩大尺度精品在线看网址 | 成人三级做爰电影| 巨乳人妻的诱惑在线观看| 操美女的视频在线观看| 午夜免费成人在线视频| 免费在线观看亚洲国产| 日韩欧美国产一区二区入口| 欧美成狂野欧美在线观看| 九色亚洲精品在线播放| 18禁美女被吸乳视频| 午夜激情av网站| 精品免费久久久久久久清纯| 亚洲一区二区三区不卡视频| 欧美日韩亚洲综合一区二区三区_| 午夜激情av网站| a在线观看视频网站| 亚洲五月色婷婷综合| 日本一区二区免费在线视频| 久久人人97超碰香蕉20202| 成年人黄色毛片网站| 亚洲熟女毛片儿| 热99re8久久精品国产| 91字幕亚洲| av福利片在线| 大型黄色视频在线免费观看| 国产精品美女特级片免费视频播放器 | 日韩精品中文字幕看吧| 搡老妇女老女人老熟妇| 在线十欧美十亚洲十日本专区| 精品免费久久久久久久清纯| 老熟妇乱子伦视频在线观看| tocl精华| 欧美亚洲日本最大视频资源| 亚洲色图 男人天堂 中文字幕| 不卡一级毛片| 真人一进一出gif抽搐免费| 久久久久久免费高清国产稀缺| tocl精华| 悠悠久久av| 久久久国产成人精品二区| 人人妻人人澡欧美一区二区 | 久久久久久免费高清国产稀缺| 大香蕉久久成人网| 激情视频va一区二区三区| 久9热在线精品视频| 日本黄色视频三级网站网址| 麻豆av在线久日| 中文亚洲av片在线观看爽| 精品欧美国产一区二区三| 老司机靠b影院| 欧美日韩一级在线毛片| 午夜福利高清视频| 无人区码免费观看不卡| 日韩三级视频一区二区三区| 高潮久久久久久久久久久不卡| 91精品三级在线观看| 亚洲欧美日韩高清在线视频| 日韩三级视频一区二区三区| 成人18禁在线播放| 午夜视频精品福利| 他把我摸到了高潮在线观看| 精品一区二区三区四区五区乱码| 亚洲一区中文字幕在线| 久久中文看片网| 国产av又大| 精品卡一卡二卡四卡免费| 在线观看www视频免费| 露出奶头的视频| 淫妇啪啪啪对白视频| 亚洲全国av大片| 日韩欧美免费精品| 亚洲精品国产色婷婷电影| 亚洲美女黄片视频| 午夜福利视频1000在线观看 | 黄频高清免费视频| 一本综合久久免费| 亚洲国产欧美一区二区综合| 天天躁夜夜躁狠狠躁躁| 色尼玛亚洲综合影院| videosex国产| 少妇 在线观看| 很黄的视频免费| 久久精品国产清高在天天线| 亚洲国产日韩欧美精品在线观看 | 人妻丰满熟妇av一区二区三区| 国产高清有码在线观看视频 | 在线国产一区二区在线| 中国美女看黄片| 女同久久另类99精品国产91| 99在线视频只有这里精品首页| 麻豆久久精品国产亚洲av| 免费在线观看黄色视频的| 国产精品久久久久久亚洲av鲁大| 性欧美人与动物交配| av欧美777| 在线播放国产精品三级| 中文字幕人妻丝袜一区二区| 我的亚洲天堂| 两个人视频免费观看高清| 亚洲专区国产一区二区| 午夜福利视频1000在线观看 | 国产不卡一卡二| 免费在线观看日本一区| 露出奶头的视频| 又黄又爽又免费观看的视频| 亚洲成av人片免费观看| 亚洲国产欧美日韩在线播放| 国产精品美女特级片免费视频播放器 | 久久国产精品影院| 日韩高清综合在线| 一区二区三区高清视频在线| 国产麻豆69| 男人舔女人下体高潮全视频| 成人18禁在线播放| 国产一卡二卡三卡精品| 欧美性长视频在线观看| 亚洲人成电影观看| 亚洲 国产 在线| 久久久久久亚洲精品国产蜜桃av| 一二三四社区在线视频社区8| 国产精品国产高清国产av| 国产麻豆成人av免费视频| 亚洲视频免费观看视频| 国产精品乱码一区二三区的特点 | 久久久国产欧美日韩av| 窝窝影院91人妻| 国产精品久久视频播放| 在线视频色国产色| 99国产精品99久久久久| 禁无遮挡网站| 亚洲国产看品久久| 国产成人精品久久二区二区免费| 一本久久中文字幕| 国产xxxxx性猛交| 精品卡一卡二卡四卡免费| tocl精华|