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

不規(guī)則波作用下斜坡堤越浪的數(shù)值模擬

2018-03-15 11:32:57李東洋張慶河焦方騫
水道港口 2018年1期
關(guān)鍵詞:模型

李東洋,張慶河,焦方騫

(天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室,天津 300072)

斜坡堤是保護(hù)沿海地區(qū)免受潮、浪襲擊的重要水工建筑物,由于斜坡堤迎浪面較緩,波浪爬高往往較大,很容易發(fā)生越浪,因此準(zhǔn)確模擬斜坡堤越浪過程并確定越浪量,對(duì)于斜坡堤的設(shè)計(jì)具有重要意義。

早期對(duì)斜坡堤越浪的研究主要以物理模型試驗(yàn)為主,王聰?shù)萚1]在總結(jié)前人物理模型試驗(yàn)的基礎(chǔ)上,研究了不規(guī)則波作用下?lián)趵藟Ω叨鹊葘?duì)越浪量的影響。近年來,隨著計(jì)算機(jī)技術(shù)和并行計(jì)算的發(fā)展,利用數(shù)值模型研究斜坡堤越浪成為港口海岸工程學(xué)科關(guān)注的熱點(diǎn)。追蹤自由表面運(yùn)動(dòng)的斜坡堤越浪數(shù)值模型主要分為兩大類,一類基于N-S方程,如Van Gant等[2]較早基于二維不可壓縮的N-S方程研究了波浪與寬肩臺(tái)可滲斜坡堤相互作用的水動(dòng)力特性,Losada等[3]基于體積平均的RANS方程和VOF方法建立了COBRAS(Cornell Breaking Waves and Structure)模型,模擬分析了不規(guī)則波在可滲斜坡堤上的越浪過程。關(guān)大瑋[4]基于CFD軟件Flow-3D建立三維數(shù)值波浪水槽,模擬研究了規(guī)則波和不規(guī)則波作用下不可滲光滑斜坡堤的越浪過程。閆科諦等[5]基于OpenFOAM軟件建立數(shù)值波浪水槽,通過網(wǎng)格細(xì)化,刻畫出柵欄板凹槽的形狀,進(jìn)而模擬了規(guī)則波作用下非可滲堤心柵欄板護(hù)面的斜坡堤越浪過程。Xiang等[6]等通過將Fluidity和Y3D耦合的方法,模擬研究了規(guī)則波作用下扭王塊體護(hù)面斜坡堤周圍的水動(dòng)力特性。另外一類模型為無網(wǎng)格方法,如葉曉文[7]基于不可壓縮SPH數(shù)學(xué)模型研究了規(guī)則波沿可滲斜坡堤的爬坡和越浪過程, Ren[8]基于弱可壓的光滑粒子流體力學(xué)方法(WCSPH)模擬了波浪與可滲潛堤的相互作用,并分析了堤體內(nèi)部的水動(dòng)力特性。溫鴻杰[9]基于改進(jìn)的SPH方法,模擬研究了規(guī)則波在出水堤上的破碎和越浪過程。

總體上來看,目前關(guān)于斜坡堤越浪數(shù)值模擬的研究還主要集中于規(guī)則波,對(duì)護(hù)面塊體進(jìn)行三維全尺度模擬且考慮護(hù)面塊體墊層以下可滲特性的工作也比較少見,這使得利用數(shù)值模型解決斜坡堤越浪等實(shí)際工程問題仍受到一定的制約。因此,利用數(shù)值模型展現(xiàn)真實(shí)海況下斜坡堤越浪過程,需要進(jìn)一步探索。考慮到SPH類模型在壓力脈動(dòng)等方面還存在有待解決的問題,本文將利用Jesus[10]和Lara等[11]基于OpenFoam開源軟件開發(fā)的ihFoam求解器,通過OpenFoam網(wǎng)格劃分工具snappyHexMesh精細(xì)刻畫網(wǎng)格的方法生成扭王字塊,模擬不規(guī)則波作用下扭王塊體護(hù)面斜坡堤越浪過程,以使斜坡堤越浪的數(shù)值模擬在解決工程問題時(shí)達(dá)到一定的實(shí)用性,為斜坡堤的設(shè)計(jì)提供依據(jù)。

1 控制方程與邊界條件

1.1 控制方程

本文通過將滲流運(yùn)動(dòng)的Darcy-Forchheimer方程引入到兩相流求解器interFoam中,生成全新的考慮孔隙介質(zhì)流的ihFoam求解器,采用主動(dòng)吸收式速度入口造波,通過求解體積平均的RANS方程,模擬水體和空氣兩相不可壓縮流體的運(yùn)動(dòng)。紊流模型選用RNGk-ε模型,采用VOF法捕捉自由表面。模型控制方程如下

(1)

(2)

(3)

(4)

式中:D50為組成孔隙介質(zhì)的中值粒徑;KC=umT/(nD50),其中um為最大震蕩速度;T為震蕩周期。α和β分別為Darcy-Forchheimer線性項(xiàng)系數(shù)和非線性項(xiàng)系數(shù),2011年Lara等[13]提出了α和β關(guān)于孔隙率和中值粒徑的經(jīng)驗(yàn)公式如下

(5)

當(dāng)n=1時(shí),即純水體情況時(shí),式(2)轉(zhuǎn)化為標(biāo)準(zhǔn)的雷諾時(shí)均方程。

圖1 數(shù)值水槽和邊界條件示意圖Fig.1 Numerical wave flume and boundary condition diagram

1.2 邊界條件

數(shù)值水槽和邊界條件示意圖如圖1所示,左側(cè)入口設(shè)置為主動(dòng)吸收造波邊界。右側(cè)出口設(shè)置為輻射消波邊界。數(shù)值水槽前后側(cè)邊界設(shè)置為可滑移固邊界,法向速度為0。底面及護(hù)面塊體設(shè)置為固壁邊界,數(shù)值邊界條件設(shè)置為無滑移邊界條件,壁面處的流速為0。

2 扭王塊體護(hù)面斜坡堤越浪模擬

為了驗(yàn)證數(shù)值模型的合理性,我們收集已有物理模型試驗(yàn)數(shù)據(jù)與數(shù)值模型結(jié)果進(jìn)行對(duì)比分析。

2.1 物理試驗(yàn)?zāi)P驮O(shè)置

圖2 物理模型試驗(yàn)尺寸示意圖[14]Fig.2 Physical model size diagram[14]

陳偉秋等[14]通過物模試驗(yàn)研究了不規(guī)則波作用下扭王塊體護(hù)面斜坡堤的平均越浪量。試驗(yàn)是在長(zhǎng)175 m、寬1.2 m、高1.5 m的波浪水槽中進(jìn)行的,試驗(yàn)?zāi)P偷木唧w尺寸如圖2所示,斜坡堤左右兩側(cè)坡度均為1:1.5,堤頂高度為0.5 m,斜坡堤堤底寬為2.11 m,兩端各有坡度為1:2的拋石棱體作為護(hù)底,海側(cè)斜坡鋪設(shè)一層扭王字塊,單個(gè)扭王字塊的高度是6 cm。試驗(yàn)水深0.4 m,入射波為不規(guī)則波,采用JONSWAP譜進(jìn)行模擬,有效周期為2 s,有效波高分別為6.6 cm、7.1 cm、7.6 cm。

表1 多孔介質(zhì)物理參數(shù)表Tab.1 Physical parameter of porous media

模擬過程中除扭王塊體采用全尺度模擬外,斜坡堤堤體其他組成部分均作為孔隙介質(zhì)處理,各層物理參數(shù)見表1。

2.2 三維數(shù)值水槽不規(guī)則波驗(yàn)證

試驗(yàn)中采用JONSWAP譜進(jìn)行模擬,JONSWAP譜的表達(dá)形式為

(6)

(7)

(8)

(9)

式中:S(f)為波浪譜密度函數(shù);Hs為有效波高;fp為譜峰周期;f為波浪頻率;γ為譜峰因子,取3.3;σ是波浪密度的標(biāo)準(zhǔn)差,當(dāng)f≤fp時(shí),σ=0.07,當(dāng)f>fp時(shí),σ=0.09。

3-a 波面歷時(shí)圖 3-b 計(jì)算波譜與理論波譜對(duì)比圖圖3 Hs=7.1 cm情況下數(shù)值造波結(jié)果Fig.3 Hs=7.1 cm Simulation results of free surface

不規(guī)則波采用線性波疊加法生成[15],根據(jù)等分頻率法將波浪能量分割成70份。因?yàn)閿?shù)值模擬的側(cè)壁采用的是可滑移邊界條件,避免了壁面對(duì)內(nèi)部流場(chǎng)的影響,所以可以利用寬度較窄的計(jì)算水槽反演試驗(yàn)結(jié)果,具體數(shù)值水槽設(shè)置如下:長(zhǎng)15 m,高0.6 m,寬0.24 m,網(wǎng)格水平空間步長(zhǎng)取L/100,為了準(zhǔn)確追蹤自由表面,在其靜水位上下2倍有效波高范圍內(nèi)取30層網(wǎng)格,上下網(wǎng)格漸變,寬度方向上取15層網(wǎng)格。數(shù)值模擬時(shí)間為220 s,取x=12 m(即堤腳)處測(cè)點(diǎn)進(jìn)行分析,以有效波高為7.1 cm為例,圖3-a顯示了其水位歷時(shí)曲線。取后100個(gè)波,即20~220 s的實(shí)測(cè)數(shù)據(jù)進(jìn)行快速傅里葉分析,得到不規(guī)則波浪的能量譜。圖3-b顯示了計(jì)算模擬波譜與理論波譜的對(duì)比情況,由圖可知兩者比較接近。模擬有效波高與理論有效波高誤差為2.4%,滿足《波浪模型試驗(yàn)規(guī)程》[16]要求。

圖4 AutoCAD三維建模中扭王字塊及扭王字塊排放示意圖 圖5 扭王塊體和胸墻經(jīng)snappyHexMesh處理后局部網(wǎng)格剖分示意圖Fig.4 AutoCAD three-dimensional model of accropode blocks Fig.5 Accropode blocks and parapet local mesh encryption by snappyHexMesh

2.3 扭王塊體全尺度模擬的網(wǎng)格刻畫

OpenFOAM提供了強(qiáng)大的網(wǎng)格生成工具,通過自帶的網(wǎng)格處理工具snappyHexMesh將扭王字塊和胸墻轉(zhuǎn)變成貼體網(wǎng)格模擬,可以逼近真實(shí)工況下斜坡堤上的波浪爬坡和越浪。

首先參照《防波堤設(shè)計(jì)與施工規(guī)范》[17]中規(guī)定的扭王字塊的尺寸,使用AutoCAD三維建模繪制出扭王字塊的外部輪廓。周雅等[18-19]通過物理模型試驗(yàn)研究得出扭王字塊護(hù)面塊體的排放型式對(duì)越浪量的影響,因此為了簡(jiǎn)化模型,將扭王字塊規(guī)則排放在斜坡堤上,如圖4所示。圖5顯示了精細(xì)網(wǎng)格刻畫后的扭王塊體和胸墻局部網(wǎng)格形狀。由圖可以看出數(shù)模所采用的網(wǎng)格能較好描述扭王護(hù)面塊體和胸墻的外部形狀。

6-a越浪量歷時(shí)圖6-b累計(jì)越浪量圖圖7 平均越浪量計(jì)算值與試驗(yàn)值比較圖6 Hs=7.1cm情況下越浪量計(jì)算結(jié)果Fig.6Hs=7.1cmSimulationresultsofovertoppingdischargeFig.7Comparisonofnumericalaveragedischargeandexperimentalaveragedischarge

2.4 越浪模擬與試驗(yàn)結(jié)果的比較

8-a t=174.36 s 9-a t=174.36 s

圖6顯示了有效波高為7.1 cm情況下的波浪越浪量歷時(shí)圖6-a和累計(jì)越浪量圖6-b,在有效周期一定的情況下,越浪量隨有效波高的增大而明顯增大。取20~220 s(100個(gè)波)時(shí)間內(nèi)的越浪量計(jì)算平均值以得到平均越浪量。圖7顯示了平均越浪量數(shù)值模擬結(jié)果與試驗(yàn)值的比較,由圖可知斜坡堤越浪的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果吻合較好,表明利用snappyHexMesh網(wǎng)格處理工具生成貼體網(wǎng)格的方法對(duì)扭王護(hù)面塊體進(jìn)行全尺度模擬,對(duì)墊層以下塊石利用滲流模型模擬斜坡堤越浪可以獲得較好的效果。

2.5 越浪過程和堤前附近流場(chǎng)分析

8-b t=174.56 s 9-b t=174.56 s

8-c t=174.76 s 9-c t=174.76 s圖8 最大單波越浪過程三維示意圖 圖9 堤前附近流場(chǎng)圖 Fig.8 3D process of the largest Fig.9 Variation diagram of single-wave overtopping velocity fields

圖8顯示了有效波高為7.6 cm,在174.5 s前后發(fā)生越浪時(shí)波面的三維過程,圖9顯示了Y=0.1 m,XZ剖面扭王字塊附近的流場(chǎng)變化。從圖中可以看出,越浪發(fā)生前,波浪沿著坡面向上爬升,在向上水體的帶動(dòng)下,水面抬高,波浪拍擊胸墻表面,直至越過堤頂形成越浪。在整個(gè)越浪過程中,由于扭王字塊護(hù)面的作用,堤內(nèi)的滲流流速相對(duì)很小。

從越浪過程、越浪波面和流速變化過程看,本文所建立的模型較準(zhǔn)確地反映了波浪與扭王塊護(hù)面斜坡堤的相互作用,具備了為斜坡堤設(shè)計(jì)優(yōu)化提供支持的能力。

表2模型與原型波要素

Tab.2 Wave elements under experimental model and prototype size

類別有效波高(m)有效周期(s)水深(m)實(shí)驗(yàn)室尺度0.07120.40原型尺度0.071×25=1.7752×251/2=100.40×25=10.00

2.6 原型工況的數(shù)值模擬

現(xiàn)階段對(duì)于原型尺寸下的波浪與斜坡堤相互作用的研究,大多是通過比尺縮放,進(jìn)行實(shí)驗(yàn)室物模試驗(yàn),再將試驗(yàn)值根據(jù)《波浪模型試驗(yàn)規(guī)程》換算成原型,而直接研究原型尺寸下斜坡堤越浪的工作還比較少。因此,為進(jìn)一步論證模型可否應(yīng)用于描述現(xiàn)場(chǎng)尺度條件下的斜坡堤越浪,首先依據(jù)《波浪模型試驗(yàn)規(guī)程》[16]將文獻(xiàn)[14]中的物理模型試驗(yàn)按照1∶25的比尺放大成原型工況,通過模擬該尺度下的斜坡堤越浪,并與文獻(xiàn)[14]按照重力相似計(jì)算得出的越浪量進(jìn)行比較,從而驗(yàn)證模型是否可用于模擬計(jì)算現(xiàn)場(chǎng)尺度下的斜坡堤越浪。以模型有效周期為2 s,有效波高0.071 m為例,按照重力相似準(zhǔn)則確定原型波要素如表2所示。

圖10 原型工況下累計(jì)越浪量Fig.10 The accumulated overtopping discharge under prototype size

原型模擬數(shù)值水槽長(zhǎng)×寬×高=375 m×6 m×15 m,依然按照上述方法建模。圖10顯示了原型尺度下累計(jì)越浪量的數(shù)值模擬結(jié)果,單個(gè)越浪發(fā)生時(shí)其實(shí)歷時(shí)很短,可以看成是一個(gè)瞬時(shí)過程,圖中折線每次的抬升就代表發(fā)生了一次越浪,“臺(tái)階”的高度即為單個(gè)越浪發(fā)生時(shí)的越浪量,從而可以得出,單波最大越浪量發(fā)生在950 s左右,單波最大越浪量約為1.3 m3/m。取100~1 100 s(100個(gè)波)時(shí)間內(nèi)的越浪量計(jì)算平均值得到平均越浪量,原型尺度下數(shù)值計(jì)算的越浪量為7.89×10-3m3/m/s。根據(jù)模型試驗(yàn)越浪量折算到原型的平均越浪量為6.68×10-5×251.5=8.36×10-3m3/m/s,與原型數(shù)值模擬結(jié)果較為接近。合理的模型試驗(yàn)可以反映原型結(jié)果是比較公認(rèn)的,因此,從原型直接模擬結(jié)果與模型試驗(yàn)折算到原型的越浪量比較吻合來看,利用OpenFOAM直接模擬原型工況下的斜坡堤越浪是可行的,并推薦采用原型尺寸進(jìn)行數(shù)值研究,進(jìn)而優(yōu)化確定設(shè)計(jì)階段中斜坡堤的斷面型式和頂高程等。

3 結(jié)語

本文基于OpenFOAM建立數(shù)值波浪水槽,控制方程采用了體積平均的RANS方程和包含非線性項(xiàng)的滲流方程描述斜坡堤內(nèi)部的水體滲流流動(dòng)。利用AutoCAD三維建模繪制出扭王塊體的外部輪廓,利用snappyHexMesh網(wǎng)格劃分工具對(duì)斜坡堤扭王護(hù)面塊體進(jìn)行全尺度網(wǎng)格精細(xì)劃分。采用建立的模型模擬研究了正向入射不規(guī)則波與扭王塊體護(hù)面斜坡堤的相互作用,模擬結(jié)果表明,平均越浪量的計(jì)算值與實(shí)驗(yàn)結(jié)果吻合較好。通過將實(shí)驗(yàn)室尺度放大成原型尺度,模擬研究了原型條件下的斜坡堤越浪,從計(jì)算結(jié)果可以看出,平均越浪量的計(jì)算值與模型通過重力相似準(zhǔn)則折算成原型的結(jié)果較為吻合。從實(shí)驗(yàn)室尺度和原型尺度模擬結(jié)果看,對(duì)護(hù)面塊體進(jìn)行全尺度模擬的數(shù)值波浪水槽目前已可以較為合理地描述復(fù)雜護(hù)面塊體斜坡堤的越浪過程,可以應(yīng)用于設(shè)計(jì)階段對(duì)斜坡堤斷面型式和頂高程的優(yōu)化確定。

需要指出的是,本文僅給出了波浪正向入射情況下斜坡堤的越浪模擬結(jié)果,有關(guān)護(hù)面塊體本身的穩(wěn)定性以及斜向入射波與斜坡堤的相互作用有待進(jìn)一步研究。

[1]王聰, 陳國平, 嚴(yán)士常, 等. 不規(guī)則波作用下斜坡堤越浪量試驗(yàn)研究[J]. 水運(yùn)工程, 2017(2): 49-52.

WANG C, CHEN G P, YAN S C,et al. Experimental study on wave overtopping of sloping seawall under irregular wave action[J]. Port & Waterway Engineering, 2017(2): 49-52.

[2]Van Gent M R A. Wave action on and in permeable structures[D]. The Netherlands: Delft University, 1995.

[3]Losada I J, Lara J L, Guanche R, et al. Numerical analysis of wave overtopping of rubble mound breakwaters[J]. Coastal Engineering, 2008, 55(1): 47-62.

[4]關(guān)大瑋. 規(guī)則波與不規(guī)則波的海堤越浪數(shù)值模擬[D]. 廣州:華南理工大學(xué), 2016.

[5]閆科諦, 張慶河. 柵欄板護(hù)面斜坡堤越浪數(shù)值模擬研究[J]. 中國港灣建設(shè), 2016, 36(2):16-19.

YAN K D, ZHANG Q H. Numerical simulation of overtopping on sloping dike with fence panels[J]. China Harbour Engineering, 2016, 36(2):16-19.

[6]Xiang J, Latham J P, Vire A, et al. Coupled fluidity/Y3D technology and simulation tools for numerical breakwater modelling[J]. Coastal Engineering Proceedings, 2012, 1(33):66.

[7]葉曉文. 斜坡式防波堤越浪過程的SPH模擬[D]. 大連:大連理工大學(xué), 2011.

[8]Ren B, Wen H, Dong P, et al. Improved SPH simulation of wave motions and turbulent flows through porous media[J]. Coastal Engineering, 2016, 107: 14-27.

[9]溫鴻杰. 非線性波浪與可滲結(jié)構(gòu)相互作用的SPH模型[D]. 大連:大連理工大學(xué), 2016.

[10]Jesus M D, Lara J L, Losada I J. Three-dimensional interaction of waves and porous coastal structures. Part I: Numerical model formulation[J]. Coastal Engineering, 2012, 64:57-72.

[11]Lara J L, del Jesus M, Losada I J. Three-dimensional interaction of waves and porous coastal structures. Part II: Experimental validation[J]. Coastal Engineering, 2012, 64: 26-46.

[12]Higuera P, Lara J L, Losada I J. Three-dimensional interaction of waves and porous coastal structures using OpenFOAM?. Part I: Formulation and validation[J]. Coastal Engineering, 2014, 83: 243-258.

[13]Lara J L, Losada I J, Maza M, et al. Breaking solitary wave evolution over a porous underwater step[J]. Coastal Engineering, 2011, 58(9):837-850.

[14]陳偉秋, 陳凌彥, 王登婷,等. 斜坡堤后坡砌石護(hù)面穩(wěn)定厚度的模型試驗(yàn)研究[J]. 水運(yùn)工程, 2016(6): 93-98.

CHEN W Q, CHEN L Y, WANG D T, et al. Experimental study on thickness of inner block revetment of sloped seawall[J]. Port & Waterway Engineering, 2016(6): 93-98.

[15]俞聿修, 柳淑學(xué). 隨機(jī)波浪及其工程應(yīng)用[M]. 大連:大連理工大學(xué)出版社, 2011.

[16]JTJ/T234-2001,波浪模型試驗(yàn)規(guī)程[S].

[17]JTS154-1-2011, 防波堤設(shè)計(jì)與施工規(guī)范[S].

[18]周雅, 林登榮, 李慶銀, 等. 不規(guī)則波作用下斜坡堤越浪量試驗(yàn)研究[J]. 水道港口, 2016(4):331-335.

ZHOU Y, LIN D R, LI Q Y, et al. Experimental research of wave overtopping on sloping dike under irregular waves[J]. Journal of Waterway and Harbor, 2016(4): 331-335.

[19]鐘雄華,陳國平,嚴(yán)士常,等.不同擺放方式扭王字塊體穩(wěn)定性研究[J].水道港口,2016(5):479-483.

ZHONG X H, CHEN G P, YAN S C,et al. A stability analysis of different stack types of accropode[J]. Journal of waterway and harbor, 2016(5):479-483.

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 精品三级在线| 波多野结衣无码AV在线| 精品成人一区二区三区电影 | 91麻豆精品国产高清在线| 久久成人免费| 一区二区三区毛片无码| 国产不卡国语在线| 国产精品自在线拍国产电影| 久久这里只有精品8| 欲色天天综合网| 国产精品短篇二区| 亚洲中文字幕久久精品无码一区| 日韩 欧美 小说 综合网 另类| 狠狠亚洲五月天| 91精品日韩人妻无码久久| 日韩高清成人| 国产人成网线在线播放va| 午夜久久影院| 欧亚日韩Av| 18禁色诱爆乳网站| 97se亚洲综合在线天天| 天天综合天天综合| 国产欧美中文字幕| 欧美无遮挡国产欧美另类| a级毛片网| 欧美一级专区免费大片| 久久精品亚洲中文字幕乱码| 乱人伦中文视频在线观看免费| 久久久久中文字幕精品视频| 手机精品福利在线观看| 强奷白丝美女在线观看| 久久久久久久蜜桃| 成人久久精品一区二区三区| 亚洲成人网在线观看| 毛片免费观看视频| 精品亚洲麻豆1区2区3区| 色综合天天娱乐综合网| 天天爽免费视频| 天天操天天噜| 欧美一级片在线| 欧美午夜网| 国产微拍精品| 国产一级裸网站| 国产理论最新国产精品视频| 色久综合在线| 亚洲一级毛片在线观| 久草视频精品| 久久中文字幕2021精品| 久久福利片| 日本午夜三级| 欧美亚洲另类在线观看| 国产呦视频免费视频在线观看 | 深爱婷婷激情网| 色135综合网| 黄色一及毛片| 午夜福利网址| 国模视频一区二区| 91高清在线视频| 国产乱子伦无码精品小说| 国产白浆视频| 欧美国产日韩在线播放| 91精品网站| 日韩欧美国产中文| 久久熟女AV| 欧美日韩在线观看一区二区三区| 乱人伦视频中文字幕在线| 久久不卡国产精品无码| 日韩小视频在线播放| 91成人在线免费视频| 看国产毛片| 激情网址在线观看| 全裸无码专区| 婷婷色狠狠干| a级高清毛片| 四虎成人免费毛片| 67194亚洲无码| 亚洲精品成人福利在线电影| 欧美国产日韩另类| 日本AⅤ精品一区二区三区日| 日韩东京热无码人妻| 欧美国产日韩另类| 国产麻豆精品在线观看|