999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

三維交錯沙波上的紊流特性數值模擬

2017-07-05 11:38:32何立群陳孝兵
水利水電科技進展 2017年4期
關鍵詞:結構模型

何立群,陳孝兵,陳 力,趙 堅

(1.河海大學水利水電學院,江蘇 南京 210098;2.河海大學水文水資源與水利工程科學國家重點實驗室,江蘇 南京 210098)

?

三維交錯沙波上的紊流特性數值模擬

何立群1,陳孝兵1,陳 力2,趙 堅1

(1.河海大學水利水電學院,江蘇 南京 210098;2.河海大學水文水資源與水利工程科學國家重點實驗室,江蘇 南京 210098)

為了揭示復雜河床沙波上的紊流結構特性,利用Fluent軟件,基于雷諾平均N-S方程及k-ω紊流模型構建了描述三維沙波上水流運動規律的數學模型,并采用試驗數據對模型進行了驗證。結果表明:所采用的數值模擬模型能夠良好的刻畫沙波上的紊流特征,沙波的三維性是影響水流紊流結構的重要因素;相比于縱向相位差φ1,三維交錯沙波的橫向相位差φ2對剪切應力與床面阻力的影響更加顯著,隨著φ2增大,沙波逐漸交錯,剪切應力與床面阻力均減小。

三維沙波;k-ω模型;紊流;剪切應力;床面阻力

沙波紊流結構的研究已取得了許多有價值的成果。毛野等[4]在循環水槽內利用粒子圖像速度場儀將流動與電腦圖像相結合,展示了沙波床面上紊流擬序結構;Noguchi 等[5]在運動沙波床面上測量了不同時間、不同床面位置的水流瞬時流速,認為波峰處的時均流速較明渠變化不大,波谷處由于反向流的作用,時均流速變化較大;陳孝兵等[6]利用NaCl示蹤試驗研究了不同床面形態驅動下的地表水-地下水水交換過程,認為近壁面流場分布是影響地下水流場的主要因素;黃華東等[7]以天然的階梯-深潭結構為研究對象,采用k-ω紊流模擬方法探討了不同流量條件下河床底部壓力結構特性及安全性質;Stoesser 等[8-9]揭示出二維沙波表面的再附著點附近,流線彎曲并且向外運動會導致泡漩渦發展成發夾渦,他們將之歸因于二次流的非穩定性;Venditti等[10-11]則統計了不同床面條件下水流分離區長度,得到分離長度大約為4~6倍波高的結論;Lefebvre等[12]則較為系統地討論了不同二維沙波床面上水流分離現象與剪切應力之間的響應關系;唐立模等[13]對有關明渠紊流與床面形態的相互作用研究進行了分析和總結。

從已有的研究來看,無論是數值模擬還是物理模型試驗,絕大多數研究集中在二維條件下,僅有極少數的研究涉及三維沙波情況,對于三維沙波條件下的紊流結構及動床阻力尚缺乏深入的研究。實際上,天然河流中沙波往往呈現為復雜交錯的三維結構,且床面的三維性將引發水流的次生流,這在某些情況下可能會加強或減弱床面上的紊流強度。有限的研究已揭示出三維床面上的紊流結構明顯有別于二維情況。Parsons等[14]通過分析Parana River底部沙波上的流速分布,討論了不同沙波形態對水流的影響,結果表明,三維沙波與二維沙波相比,流速變化較大,其研究結果揭示了真實河流中不同沙波形態之間的紊流結構異同;Omidyeganeh等[15]利用大渦模擬方法研究了三維床面上流場結構,并分析統計了不同床面形態間紊流結構的差異,結果表明三維沙波流場中存在比二維流場更多的渦結構,三維床面增加了床面的水流拖曳力;Nezu 等[16]研究了沙波背水面的三維渦結構,發現分離渦的平均周期隨雷諾數的增加而減小,泡漩渦的平均周期隨雷諾數的增加而增加,泡漩渦的渦強度更大;Maddux等[2-3]通過在室內水槽中構建三維沙波來討論二維、三維沙波形態的不同對流場分布以及沙波表面壓力分布的影響,其試驗研究方法為更進一步研究復雜三維床面上的水沙過程提供了寶貴的經驗;最近,Chen等[17]在Maddux等[2-3]物理模型基礎上,構建了三維沙波床面條件地形下數值計算模型,探討了不同流動狀態下床面阻力的變化規律以及床面壓力分布狀態和大小對來流條件的響應,并進一步分析了地表紊流對淺層地下水運動的影響,這一研究豐富了Maddux等[2-3]的物理模型研究成果。

總體而言,當前對于沙波上紊流結構的研究主要著眼于二維條件[4-8,17]或某一固定三維沙波情況,對沙波形態不斷變化的情況下,河流紊流結構變化的規律性研究還亟待加強。本文針對復雜的床面沙波形態,基于數學模型,分析三維沙波規律性變化條件下相應紊流結構特征,探討三維沙波變化對流場結構、床面阻力及剪切應力分布的影響規律。

1 模擬方法

1.1 控制方程

Patel等[19-20]對比了相同模型條件下,二維水槽試驗與數值模擬中沙波上方流速分布、漩渦大小、阻力系數等表征參數,結果表明k-ω紊流模型對具有顯著渦流結構的流動現象具有較好的適用性,能較好地刻畫沙波上方紊流結構。

對于三維穩態不可壓縮流,雷諾平均N-S方程為

=0

(1)

(2)

(3)

雷諾應力與紊動能k、比率耗散系數ω有關,其計算式為

(4)

其中νt=k/ω

式中δij為克羅里克符號。

紊流模型的k方程和ω方程分別為

(5)

(6)

其中μt=ρvtα=5/9β=3/40

β*=9/100σk=σω=1/2

1.2 模型說明

復雜的沙波河床地形結構會體現在沙波形態特征及其交錯分布狀況上。已有研究表明,天然河流交錯沙波可以采用兩組簡諧波疊加生成,通過控制順水流方向簡諧波的縱向相位差φ1和垂直于水流方向的簡諧波橫向相位差φ2衍生出多種沙波結構[18]。圖1(a)(b)(c)給出了3種典型的交錯分布沙波形態。

圖1 部分典型沙波形態及φ1、φ2變化對沙波影響示意圖

本文模擬中,φ1、φ2的變化范圍分別為90°~270°和0°~180°,以30°為公差遞增,不同的相位組合形成不同的沙波結構,共49個沙波模型,基本覆蓋了沙波交錯發展的各個階段。圖1(d)(e)描述了這兩組相位差交替變化下,相應的床面結構變化情況,其中A、C與B、D分別為縱向 、橫向波長中點處截面。不難看出,對比A、C截面,當φ1變化時,沙波基本形態保持一致,只改變相對高差;對比B、D截面,當φ2變化時,沙波由相互平行逐漸過渡到相互交錯,床面基本形態發生改變;可見,床面的沙波交錯式分布狀況取決于橫向相位差φ2。為了突出交錯分布的床面對水流結構的影響,在模型中保持沙波單元橫縱波長比以及波高為常量進行數值模擬計算,這一假定基本符合天然河流中沙波狀態。

1.3 計算網格與邊界條件

由于沙波在空間上具有一定的周期性和對稱性,為了優化計算時間,對模型計算范圍作如下考慮:對于每一種沙波地形,在縱向、橫向上均截取一個周期作為計算域;在垂向上,選擇遠大于波高的水深值以保證波峰處截面上弗勞德數Fr?1,水流為亞臨界流條件,此時底部沙波對水流自由表面影響可以忽略[11]。圖2為典型沙波模型及其邊界條件,模型尺寸為0.76 m×0.5 m×0.5 m;邊界條件設置為:水流左右兩側面為對稱邊界,底部為不透水邊界,沿著水流方向的上下游邊界為周期性邊界,上下游邊界之間的壓力差為水流流動的唯一驅動力,水體頂面為光滑對稱邊界。利用GAMBIT軟件進行網格劃分,沙波表面設置34層邊界層網格,第一層高0.1 mm,增長速率1.08倍,邊界層網格厚度約為1.7 cm,邊界層網格之上網格最小高度1 mm,最大高度2 mm(圖2(a))。整個模型的單元數約為120萬,經過網格獨立性檢驗后,在此數量級的單元數下,解的精度不再因網格的進一步加密而有顯著變化,滿足計算要求,能較好地刻畫沙波上方壓力分布、背水漩渦、次生流等紊流結構特性(圖2(b))。

圖2 模型網格、邊界條件設置及典型流線分布示意圖

2 模型驗證

模型驗證基于Maddux等[2-3]的水槽試驗,試驗采用的三維沙波尺寸為1.6 m×0.9 m,平均波高為0.04 m。選取其中T3工況的試驗結果進行模型驗證,該試驗工況的平均水深為0.561 m,平均流速為0.261 m/s。通過設置相同的計算模型尺寸,采用上文所描述的邊界條件設置和網格加密要求,通過人工調節進出口壓力差,使水流平均流速接近水槽試驗平均流速。通過對比相同截面上流速分布,來判斷模型計算的合理性。模擬結果如圖3所示,在沙波上方數值計算結果略大于試驗結果,數值模型計算中能得到更大的背水漩渦。

圖3 兩個不同截面流場試驗及模擬流速結果對比

兩個截面的流速均方根誤差RMSE分別為0.031 m/s與0.028 m/s,遠小于試驗平均流速,得到的結果與試驗結果比較吻合,采用的數值計算方法可靠。

圖4 不同雷諾數與沙波形態時渦量等值線

圖5 xz剖面流線分布及x方向流速云圖(單位:m/s)

3 結果分析與討論

3.1 流場結構

通過分析床上渦結構、典型截面上的流場來探討三維沙波上的流場結構;通過分析渦量等值面的分布情況及形態差異,可以非常直觀的獲得水流在近壁面處流動情況。定義雷諾數Re=Uh/ν,式中U為平均流速、h為沙波波高、ν為運動黏滯系數,圖4給出了不同雷諾數和沙波形態下,渦結構等值面圖。由圖4可知,渦量等值面呈現兩種形態:第一種為分布在波峰與波谷之間的管狀結構,說明水流具有明顯的漩渦結構;第二種為分離區附近的類平面結構,該區域流場結構變化劇烈,渦量變化明顯。在沙波波峰上方流場結構受沙波形態影響較小,流線方向主要朝向水流方向,此處水流渦量較小;而在波峰以下近壁面區域,流場受沙波形態直接影響,剪切應力分布不均,水流呈現為繞流、漩渦等復雜的形態(圖2(b))。

對比圖4(a)(b)(c)可知,在相同沙波形態條件下,隨著雷諾數增大,渦結構形態分布基本一致。剪切層在波峰處生成,并逐漸向下游擴散[18],因此波峰后渦結構較為復雜多變。同時,由于沙波的阻擋,部分水流沿著波峰線流動,與背水漩渦相結合,形成螺旋前進的三維繞流結構;部分水流則直接越過波峰線,形成水沙分離等紊流現象。對比圖4(d)(e)(f)可知,隨著φ2增大,連續的兩個沙波由平行逐漸過渡為相互交錯,由于沙波的阻擋,管狀形態的渦量等值面減少,而由于沙波形態復雜度增大,近壁面處水流結構更為復雜,第二種形態的渦量等值面增多。

圖5為相同雷諾數條件下,不同沙波形態中x方向流速分布云圖以及流線分布。統計所有截面中分離區長度L,結果表明,L與波高h的比值范圍為3.75~5.5。Noguchi等[5]實測了沙波形成各個階段的水流特性,分析了沙波形態的影響,認為隨著沙波的發展,分離長度增大,趨近于5.5 倍波高,本文的分析結果與其接近,略有差異的原因在于模擬過程中,模型尺寸、水流強度都不盡相同。

當φ2由0°增加至180°時,沙波逐漸交錯,近壁面流態復雜,水流相互碰撞,逆流向水流受到限制,流速減小。因此,相比于圖5(a)(b),圖5(c)(d)整體流速分布均勻且x負方向流速更小。

圖6為不同沙波條件下,yz剖面流線分布及y方向流速云圖。在沙波上方,存在y、z方向流速矢量,這是由三維沙波導致的三維流動結構造成的。由于沙波阻力與主流區水體剪切力的共同作用,在波峰上方生成順時針渦體,同時在波谷上方生成逆時針渦體。在波峰與波谷連接處,即為兩種渦體交界處。

圖6 yz剖面流線分布及y方向流速云圖(單位:m/s)

對于圖6(a)(b),順時針渦體遠大于逆時針渦體,占據整個截面的75%。原因在于波峰上方水體受沙波阻力影響較小,能量耗散較小,故流速矢量較大,而波谷上方水體則相反。對比圖6(c)(d),隨著φ2增大,沙波相互交錯導致阻力減小,從而導致波谷上方逆時針渦體增大。由于兩種渦體的相互作用,導致沙波上方水體y、z方向流速變化較小。

由圖5、圖6可知,沙波的三維性促成了流場中垂直、平行于流動方向的渦結構;可以預見,這種三維紊流結構將會影響床沙質縱向、橫向的推移,沙波交錯性與上覆水體流場分布具有很強的互動性。

3.2 雷諾剪切應力

雷諾剪切應力為

(7)

圖7(a)(b)分別為不同φ1、φ2時,雷諾數與流場中最大剪切應力關系曲線。由圖7可看出,剪切應力τ與雷諾數之間存在冪函數關系τ=aReb(式中a、b為系數,圖7(a)中a=1.35×10-11,b=2.035),不同φ值對剪切應力影響不同。當φ2保持不變,φ1由90°增加到270°,Re=13 000時,(τmax-τmin)/τmin=15.7%。由于φ1變化并沒有改變沙波形態的起伏度,當沙波波高一致時,波峰處的時均流速變化不大,故當雷諾數相同時,雷諾應力變化較小。

圖7 雷諾數與剪切應力關系

圖7(b)為沙波形態按φ2變化時,剪切應力變化規律,其中虛線的經驗系數a=3.02×10-12、b=2.13,實線的經驗系數與圖7(a)一致。由圖7可知,當雷諾數一致時,φ2增大,沙波逐漸交錯發展,剪切應力減小,在曲線中表現為b≈2,但a隨著φ2增大而減小。根據上文分析,沙波的交錯發展導致上方流速分布均勻,故水流之間相互作用減小。

無論床面是否交錯,最大雷諾剪切應力均分布在沙波背水坡上。水流越過沙波頂點后,沙波對流場的直接影響不復存在[4],由于沒有沙波的作用,背水坡處流場結構很不穩定。如圖4、圖5所示,在背水坡處渦結構較為明顯,流線復雜多變,說明流速變化劇烈,因此剪切應力較大。

3.3 沙波床面阻力

床面阻力大小直接影響河道泄流能力、輸沙能力,對水位預報和河床變形預報十分重要。床面形態隨水流強度變化而變化,而床面阻力又隨床面形態的變化而變化。目前,對于沙波交錯性所帶來的床面阻力分布方面的研究極少。

圖8給出了不同床面形態條件下總形狀阻力與摩擦阻力分布趨勢。顯然,φ1變化對沙波的總形狀阻力與摩擦阻力影響較小;而當沙波逐漸交錯,即φ2增加時,摩擦阻力與總形狀阻力均明顯減小。總形狀阻力與摩擦阻力的極大值與極小值分別出現在φ2=0°、180°時。這是由于沙波逐漸交錯使得縱向截面上的沙波起伏度降低(參見圖1(e)),對水流阻力干擾逐漸減小。

圖8 不同床面形態下阻力分布

圖8中的總形狀阻力和摩擦阻力趨勢面的最佳擬合公式分別為

Ffd=a/(1+bφ1+cφ2)

(9)

(10)

若取Re=6 000,式(9)系數為a=0.15,b=0.005,c=0.015,復相關系數R2=0.973;方程(10)系數為a=0.07,b=-6.82,c=-6.65,復相關系數R2=0.968。

4 結 語

基于雷諾平均N-S方程及k-ω紊流模型構建了描述三維沙波上水流運動規律的數學模型,試驗結果表明,所建立的模型能準確刻畫沙波上方水流運動特性。對模型交錯性變化與近壁面流場結構、應力分布之間關系進行對比分析,結果表明:三維沙波的交錯性是影響紊流結構的重要因素;在相同水動力條件下,流場最大剪切應力、床面阻力對φ1敏感性較低,φ2增大將使剪切應力與動床阻力減小。

本文僅著重對三維交錯性沙波對水流影響做了分析討論,而實際上沙波波高、波長等因素也是影響水沙運動的重要因素,今后需要進一步探討更為復雜床面形態變化與近壁面紊流結構、泥沙輸移之間的互動關系。

[ 1 ] PARSONS D R,BEST J.Bedforms: views and new perspectives from the third international workshop on marine and river dune dynamics (MARID3)[J].Earth Surface Processes & Landforms,2013,38(3): 319-329.

[ 2 ] MADDUX T B,NELSON J M,MCLEAN S R.Turbulent flow over three-dimensional dunes: 1.free surface and flow response[J].Journal of Geophysical Research: Earth Surface,2003,108(F1): 6009(1-10).

[ 3 ] MADDUX T B,MCLEAN S R,NELSON J M.Turbulent flow over three-dimensional dunes: 2.fluid and bed stresses[J].Journal of Geophysical Research: Earth Surface,2003,108(F1): 6010(1-11).

[ 4 ] 毛野,張志軍,袁新明,等.沙波附近紊流擬序結構特性初步研究[J].河海大學學報(自然科學版),2002,30(5): 56-61.(MAO Ye,ZHANG Zhijun,YUAN Xinming,et al.Characteristics of turbulent coherent structures over sastrugi in open channels[J].Journal of Hohai University (Natural Sciences),2002,30(5): 56-61.(in Chinese))

[ 5 ] NOGUCHI K,NEZU I,SANJOU M.Turbulence structure and fluid-particle interaction in sediment-laden flows over developing sand dunes[J].Environmental Fluid Mechanics,2008,8(5/6): 569-578.

[ 6 ] 陳孝兵,趙堅,李英玉,等.床面形態驅動下潛流交換試驗[J].水科學進展,2014,25(6): 835-841.(CHEN Xiaobing,ZHAO Jian,LI Yingyu,et al.Experimental study of bedform-driven hyporheic exchange[J].Advances in Water Science,2014,25(6): 835-841.(in Chinese))

[ 7 ] 黃華東,漆力健,余國安,等.階梯-深潭結構潭底壓力特征研究[J].長江科學院院報,2014,31(8): 50-54.(HUANG Huadong,QI Lijian,YU Guoan,et al.Characteristics of pressure at the bottom of step-pool[J].Journal of Yangtze River Scientific Research Institute,2014,31(8): 50-54.(in Chinese))

[ 8 ] STOESSER T,BRAUN C,VILLABA G,et al.Turbulencestructures in flow over two dimensional dunes[J].Journal of Hydraulic Engineering,2008,134(1): 42-55.

[ 9 ] 李佳佳,李志偉,張長寬,等.基于DANS方程的粗糙床面明渠水力特性研究[J].河海大學學報(自然科學版),2014,42(3): 217-222.(LI Jiajia,LI Zhiwei,ZHANG Changkuan,et al.Research on hydraulic characteristics of open-channel flow over rough bed based on DANS equations[J].Journal of Hohai University (Natural Sciences),2014,42(3): 217-222.(in Chinese))

[10] VENDITTI J G.Turbulent flow and drag over fixed two-and three-dimensional dunes[J].Journal of Geophysical Research: Earth Surface,2007,112(4): 1-21.

[11] LOPEZ F,FERNANDEZ R,BEST J.Turbulence and Coherent flow structures associated with bed form amalgamation: an experimental study of the ripple-dune transition[C]//Building Partnerships.New York: ASCE,2000: 1-10.

[12] LEFEBVRE A,PAARLBERG A J,WINTER C.Flow separation and shear stress over angle-of-repose bed forms: a numerical investigation[J].Water Resources Research,2014,50(2): 986-1005.

[13] 唐立模,孫會東,劉全帥.明渠紊流與床面形態相互作用研究進展[J].水利水電科技進展,2015,35(2): 77-84.(TANG Limo,SUN Huidong,LIU Quanshuai.Research development of the interaction between turbulence structure and bedforms in open channel[J].Advances in Science and Technology of Water Resources,2015,35(2): 77-84.(in Chinese))

[14] PARSONS D.R,BEST J.L,ORFEO O,et al.Morphology and flow fields of three-dimensional dunes,Rio Paran,Argentina: results from simultaneous multibeam echo sounding and acoustic Doppler current profiling[J].Journal of Geophysical Research: Earth Surface,2005,110(4): 1-9.

[15] OMIDYEGANEH M,PIOMELLI U.Large-eddy simulation of three-dimensional dunes in a steady,unidirectional flow:part 1 turbulence statistics[J].Journal of Fluid Mechanics,2013,721(4): 454-483.

[16] NEZU I,ADOTA A.Three dimensional structure of space time correlation on coherent vortices generated behind dune crest.[J].Journal of Hydraulic Researches,1999,37(1): 945-958.

[17] CHEN XIAOBING,CARDENAS M B,CHEN LI.Three-dimensional versus two-dimensional bed form-induced hyporheic exchange[J].Water Resources Research,2015,51(4): 2923-2936.

[18] RUBIN D M.Cross-bedding,bedforms,and paleocurrents[M].Tulsa:Society of Economic Paleontologists and Mineralogists,1987.

[19] PATEL V C,YOON J Y.Application of turbulence models to separated flow over rough surfaces[J].Journal of Fluids Engineering,1995,117(2): 234-241.

[20] YOON J Y,PATEL V C.Numerical model of turbulent flow over sand dune[J].Journal of Hydraulic Engineering,1996,122(1): 10-18.

Numerical simulation of turbulent characteristics over interlaced three-dimensional sand waves

HE Liqun1, CHEN Xiaobing1, CHEN Li2, ZHAO Jian1

(1.College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China; 2.State key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China)

In order to investigate the structural characteristics of turbulent flow over sand waves on a complicated riverbed, a numerical model for simulating flow over three-dimensional sand waves was established using the Fluent software based on the Reynolds-averaged Navier-Stokes (RANS) equations and thek-ωturbulence model. The model was verified using experimental data. The results show that the established model can depict the turbulent characteristics over sand waves well. The three-dimensional properties of the sand waves are important factors influencing the turbulent structure of flow. The shear stress and bed form resistance are much more sensitive to the transversal phase differenceφ2than to the longitudinal phase differenceφ1. With the increase ofφ2, the sand waves become interlaced, and the shear stress and bed form resistance both decreased.

three-dimensional sand waves;k-ωmodel; turbulent flow; shear stress; bed form resistance

國家自然科學基金(41401014,51279045);中央高校基本科研業務費專項(2016B03414)

何立群(1979—),男,碩士研究生,主要從事潛流交換研究。E-mail:wwwhlq1@qq.com

陳孝兵(1985—),男,副教授,博士,主要從事潛流交換研究。E-mail:x.chen@hhu.edu.cn

10.3880/j.issn.1006-7647.2017.04.004

TV131

A

1006-7647(2017)04-0019-06

2016-09-27 編輯:鄭孝宇)

猜你喜歡
結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 亚洲精品片911| 亚洲成在线观看| 日韩a在线观看免费观看| 亚洲中文无码h在线观看| 国产在线观看高清不卡| 国产小视频在线高清播放| 久草视频中文| 国产aaaaa一级毛片| 亚洲一区二区三区麻豆| 人妻丝袜无码视频| 国产第一福利影院| 51国产偷自视频区视频手机观看| 午夜视频免费一区二区在线看| 国产亚洲欧美在线中文bt天堂| 91精品国产综合久久不国产大片| 一本大道无码日韩精品影视| 丝袜亚洲综合| 无码专区国产精品第一页| 亚洲综合专区| 国产又大又粗又猛又爽的视频| 国产精品毛片一区| 亚洲高清国产拍精品26u| 亚洲天堂网在线播放| 国产91在线免费视频| 国产精品流白浆在线观看| 国产美女丝袜高潮| 日韩专区欧美| 亚洲,国产,日韩,综合一区 | 一区二区无码在线视频| 免费观看无遮挡www的小视频| 国产综合网站| 99久久国产综合精品2023| 黄色网在线免费观看| 男女性色大片免费网站| 麻豆精品在线播放| 国产精品永久免费嫩草研究院| 亚洲天堂网在线观看视频| 国产 日韩 欧美 第二页| 欧美天堂在线| 久久亚洲综合伊人| 亚洲成人黄色在线观看| 欧美日韩成人在线观看| 婷婷激情五月网| 亚洲性影院| 美女裸体18禁网站| 国产黄视频网站| 狠狠色婷婷丁香综合久久韩国 | 欧美国产中文| 小说 亚洲 无码 精品| www精品久久| 日韩精品高清自在线| 高清码无在线看| 制服丝袜在线视频香蕉| 国产精品尤物在线| 国产又色又刺激高潮免费看| 第一页亚洲| 国产美女无遮挡免费视频网站 | 性视频一区| 亚洲无码免费黄色网址| 久久久久亚洲AV成人网站软件| 欧美亚洲第一页| 黄色网在线免费观看| 免费无码又爽又刺激高| 欧美在线视频不卡第一页| 久久男人资源站| 在线看AV天堂| 日韩免费无码人妻系列| 亚洲国产天堂在线观看| m男亚洲一区中文字幕| AⅤ色综合久久天堂AV色综合| 国产激情无码一区二区APP| 久久婷婷综合色一区二区| 在线观看欧美精品二区| 色哟哟国产精品| 无码国内精品人妻少妇蜜桃视频| 亚洲伊人久久精品影院| 自拍偷拍一区| 亚洲成人在线网| 精品欧美一区二区三区久久久| 91福利免费| 一区二区三区四区在线| 亚洲AV无码乱码在线观看代蜜桃 |