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

        ?

        直墻前多向不規(guī)則波波高平面分布的數(shù)值模擬研究

        2021-11-04 01:48:32羅天翔朱良生徐潤(rùn)剛程趙德
        中國(guó)港灣建設(shè) 2021年10期
        關(guān)鍵詞:直墻入射波波高

        羅天翔,朱良生*,徐潤(rùn)剛,程趙德

        (1.華南理工大學(xué)土木與交通學(xué)院,廣東 廣州 510640;2.中交第四航務(wù)工程勘察設(shè)計(jì)院有限公司,廣東 廣州 510230)

        0 引言

        由于離岸建筑物的布置在港灣工程的重要性逐漸提高,如今工程設(shè)計(jì)不僅關(guān)注波要素對(duì)直墻的作用,對(duì)于直墻前沿水域的波高平面分布也要有全面的了解,以便能夠合理地確定直墻前沿的泊位、浮碼頭等工程布置位置。實(shí)際波浪往往是不規(guī)則且方向不單一的,多向不規(guī)則波更符合天然波浪的特性。因此研究正向、斜向入射后的多向不規(guī)則波經(jīng)歷反射、疊加作用后波浪在直墻建筑物前沿水域的平面空間分布具有重要意義。

        陳漢寶[1]使用數(shù)值模擬和物理模型探討了斜向波與直立式防波堤的相互作用,波高最大值與斜角的關(guān)系,給出了對(duì)應(yīng)的波態(tài)圖,定性分析了波浪特性;馮衛(wèi)兵等[2]利用淺水定常莆田經(jīng)驗(yàn)譜分解結(jié)果等模擬非線性隨機(jī)波在直墻前的運(yùn)動(dòng)和動(dòng)力過(guò)程,形成了波高、周期等一系列極值序列、分布曲線,但二階非線性理論解不能完全表現(xiàn)近岸波浪的非線性,對(duì)波高等要素的統(tǒng)計(jì)分布停留在直墻上的概率分布;俞聿修等[3]探討了波浪斜向性和多向性對(duì)不規(guī)則波作用于直墻堤上波浪荷載的影響,并提出了實(shí)用的計(jì)算方法;駱俊彬等[4]采用物理模型試驗(yàn)證明不規(guī)則波作用下的直墻前的波谷壓力概率分布可以用韋伯分布來(lái)描述,并與相對(duì)波高有關(guān);姚國(guó)權(quán)等[5]在從堤頭多角度入射的波況下研究了規(guī)則波與不規(guī)則波反射問(wèn)題,推導(dǎo)了規(guī)則波沿直立堤外側(cè)波高分布的計(jì)算式,沿直立堤方向橫截面的不規(guī)則波和規(guī)則波的波高分布,這里未強(qiáng)調(diào)波浪的非線性;周枝榮等[6]結(jié)合物理試驗(yàn)和相關(guān)研究成果,提出斜向波作用在矩形重力墩結(jié)構(gòu)上最大總水平波浪力簡(jiǎn)化計(jì)算公式;張文海等[7]分析了統(tǒng)一方程,并以此為基礎(chǔ)建立了從深水到淺水域都有效的數(shù)值模型。劉思等[8]首先根據(jù)給定的波要素、群性參數(shù)和方向分布參數(shù),給出了多向不規(guī)則波群的數(shù)值模擬方法,在此基礎(chǔ)上利用有限元法求解改進(jìn)的Boussinesq 方程的數(shù)值計(jì)算模型,建立了數(shù)值水池在指定位置模擬所要求群性的多向不規(guī)則波入射波浪邊界條件的計(jì)算方法。欒英妮[9]采用Mike21 BW 數(shù)學(xué)模型,并選用文獻(xiàn)[10]中的計(jì)算模型進(jìn)行計(jì)算并與其提供的試驗(yàn)結(jié)果比較,分析結(jié)果表明,Mike21 BW 數(shù)學(xué)模型與前人試驗(yàn)結(jié)果相近,適用于非線性較強(qiáng)、淺水以及有限水深情況下的計(jì)算;張娜等[11]使用Mike21 BW 數(shù)學(xué)模型對(duì)波浪在緩坡上的破碎現(xiàn)象及爬坡過(guò)程進(jìn)行了精確模擬,并將此方法成功應(yīng)用于實(shí)際工程中。

        關(guān)于波浪傳播入射后直墻前波高分布的研究,國(guó)內(nèi)外學(xué)者已對(duì)其進(jìn)行了大量的理論推算、物理模型以及數(shù)值模擬,并得到有效的計(jì)算方法和經(jīng)驗(yàn)公式。但前人的研究多考慮于直墻上波浪力、反射系數(shù),對(duì)波高分布的研究討論集中在直墻與波浪交匯處以及沿堤方向的橫切面,相比前人的研究,本文完全考慮不規(guī)則波的多向性、非線性,且對(duì)直墻前整體水域平面的波高分布以及波長(zhǎng)和腹點(diǎn)、節(jié)點(diǎn)位置、間距的關(guān)系進(jìn)行了較系統(tǒng)研究。本研究采用基于Boussinesq 方程的Mike21 BW 模型,對(duì)正向、斜向入射的多向不規(guī)則波在直墻前的反射、疊加過(guò)程進(jìn)行數(shù)值模擬,探討入射波與直墻前有效波高的平面空間分布關(guān)系,不同波況的有效波高分布以及得到腹點(diǎn)、節(jié)點(diǎn)間距和位置與入射波波長(zhǎng)的關(guān)系。

        1 試驗(yàn)研究方法

        1.1 Mike21 BW 數(shù)學(xué)模型簡(jiǎn)介

        計(jì)算采用丹麥DHI 的Mike21 中的BW 波浪模型,該模型通過(guò)求解沿垂向積分的Boussinesq方程獲得沿水深的平均流速、水位變化以及波高等物理量。模型中考慮了波浪的折射、繞射、反射、底部損耗等因素,可以很好地滿足計(jì)算內(nèi)容的要求。對(duì)方程采用ADI(交替方向法)算法,式中各項(xiàng)保證二階精度。新形式的Boussinesq 方程改進(jìn)了色散關(guān)系式,使BW 模型適合模擬從深水到淺水的傳播。

        1.2 基本方程

        式中:x、y 為水平坐標(biāo),m;t 為時(shí)間,s;ξ 為高出平均水位的水面高度,m;p、q 分別為x、y 方向的流量密度;h 為水深,m;D 為平均水深,m;c為謝才阻力系數(shù),c=M·h1/6,m0.5/s;M 為曼寧系數(shù),m1/3/s;E 為紊動(dòng)“渦黏”系數(shù);g 為重力加速度,g=9.81 m/s2。

        1.3 方向譜

        模型中,本文關(guān)于方向譜的設(shè)置包括對(duì)主波向的定義、方向譜的形式以及具體傳播參數(shù)的設(shè)置,波浪頻譜選擇Jonswap 譜:

        式中:f 為波頻;f0為譜峰頻率,s-1;γ 為譜峰尖度因子;H1/3為有效波高,m;T0為平均周期,s;β為峰型參數(shù);α 為γ 的函數(shù)。方向分布類型選擇頻率相關(guān)分布,方向分布函數(shù)G(θ)具體表現(xiàn)形式為:

        式中:s 為頻率相關(guān)形狀參數(shù);θmain為主波向,(°)。具體傳播參數(shù)設(shè)置為主波向θmain,最大偏離角度Δθmax兩項(xiàng),為了簡(jiǎn)化試驗(yàn),本文中Δθmax固定為常用值30°。

        1.4 計(jì)算實(shí)例

        1.4.1 試驗(yàn)概況及計(jì)算條件

        計(jì)算水域尺度為1 000 m×7 000 m,水深8 m。造波線距離左側(cè)邊界200 m,左、上、下側(cè)邊界分別設(shè)置海綿層,右側(cè)邊界設(shè)置為直立實(shí)體結(jié)構(gòu),具體布置情況如圖1 所示。本文數(shù)值計(jì)算多向不規(guī)則波經(jīng)直立實(shí)體結(jié)構(gòu)全反射后全水域波浪的水平分布,各計(jì)算工況條件見表1,其中主波向角度θmain指波浪的主波向與造波線法向的夾角。

        圖1 試驗(yàn)地形示意圖Fig.1 Schematic diagram of experimental topography

        表1 波浪初始入射條件Table 1 Initial wave incident conditions

        1.4.2 模型設(shè)置

        采用矩形網(wǎng)格劃分,空間步長(zhǎng)Δx=Δy=1 m,時(shí)間步長(zhǎng)Δt=0.1 s。近岸采用直立實(shí)體結(jié)構(gòu),反射系數(shù)為1,模型四周在開邊界上設(shè)置海綿層,吸收超出模型范圍的波浪,造波線后面產(chǎn)生的波浪也被置于該處的海綿層吸收。

        試驗(yàn)造波方式選擇多向不規(guī)則波,Jonswap 譜參數(shù)的譜峰尖度因子γ=3.3,σa=0.07,σb=0.09。取穩(wěn)定后200 個(gè)左右波的實(shí)時(shí)波面數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析,測(cè)點(diǎn)間距1 m,采取的測(cè)點(diǎn)與直墻的最遠(yuǎn)距離為入射波的2 倍深水波長(zhǎng),造波線設(shè)置在距離左、上以及下邊界200 m 處,依據(jù)模型的計(jì)算要求,最小周期(截?cái)嘀芷冢┰O(shè)置為0.7 s。

        消波系數(shù)按式(6)計(jì)算:

        式中:α 和r 為給定的常數(shù);Nsponge為消波線的數(shù)量。海綿層的寬度一般為波長(zhǎng)的1~2 倍。

        2 數(shù)學(xué)模型驗(yàn)證

        選取Demirbilek 和Nwogu 的物模試驗(yàn)工況對(duì)數(shù)學(xué)模型進(jìn)行對(duì)比、驗(yàn)證。模型布置見圖2。試驗(yàn)從造波端至右側(cè)斜坡邊界一共布置9 個(gè)波高儀,造波端左側(cè)布置的海綿層厚度為5 m,右側(cè)斜坡邊界的坡度為1∶12。波浪在深水區(qū)經(jīng)歷一段等水深傳播至由3 段斜坡組成的復(fù)合斜坡上,深水區(qū)的水深為0.531 m,復(fù)合斜坡平均坡度為1∶10.2。經(jīng)過(guò)斜坡后波浪在等水深的淺水段傳播,淺水區(qū)的水深為0.031 m。模擬波浪為不規(guī)則波,有效波高HS=0.075 m,譜峰周期TP=1.5 s。

        圖2 Demirbilek 和Nwogu 的試驗(yàn)布置示意圖Fig.2 Schematic diagram of Demirbilek and Nwogu experimental layout

        數(shù)值模型參數(shù):空間步長(zhǎng)Δx=0.05 m,時(shí)間步長(zhǎng)Δt=0.002 s,底摩擦系數(shù)fw=0.011,換算成謝才系數(shù)(Chezy Number)為55。在出口處布置1~2倍波長(zhǎng)的海綿層。模型模擬時(shí)長(zhǎng)為900 s,模擬100 s后波浪特征趨于穩(wěn)定,穩(wěn)定后的波浪用作統(tǒng)計(jì)分析。

        沿程波高的數(shù)模值和物模值對(duì)比見圖3,試驗(yàn)結(jié)果表明波高的數(shù)模試驗(yàn)值與物模試驗(yàn)值比較一致,說(shuō)明模型可用于不規(guī)則波在復(fù)雜地形上的傳播模擬。

        圖3 Demirbilek 和Nwogu 物模試驗(yàn)值與MIKE21 BW 模型有效波高試驗(yàn)值對(duì)比Fig.3 Comparison of the experimental values of the Demirbilek and Nwogu physical model with the experimental values of the effective wave height of the MIKE21 BW model

        3 計(jì)算結(jié)果分析

        3.1 離岸方向HS 空間分布

        3.1.1 正向入射波況

        在不規(guī)則波入射直墻后反射的情況中,也將最大波高Hmax值處視為腹點(diǎn),最小波高Hmin值處視為節(jié)點(diǎn)。本文采用兩種波浪周期求得波長(zhǎng)值進(jìn)行分析,分別為譜峰周期TP對(duì)應(yīng)波長(zhǎng)LP和平均周期對(duì)應(yīng)波長(zhǎng)L。

        圖4 為對(duì)3 種主波向?yàn)?°的工況結(jié)果進(jìn)行統(tǒng)計(jì)得到的各測(cè)點(diǎn)的HS。圖中x 軸為距直墻相對(duì)距離s/LP或s/L,s 為距直墻實(shí)際距離,m;LP為入射波當(dāng)前水深譜峰周期TP對(duì)應(yīng)的波長(zhǎng)值,m;L為入射波當(dāng)前水深平均周期對(duì)應(yīng)的波長(zhǎng)值;y 軸為該點(diǎn)有效波高HS與入射波有效波高HS0的比波高HS/HS0。

        圖4 主波向?yàn)?°時(shí)直墻前的波高空間分布曲線Fig.4 Spatial distribution curve of wave height in front of a straight wall when the main wave direction is 0°

        相對(duì)距離分別為s/L、s/LP時(shí),不同工況下的波高空間分布曲線一致,其中相對(duì)距離為s/LP時(shí),其曲線的腹點(diǎn)和節(jié)點(diǎn)的分布位置與規(guī)則波入射直墻反射后形成的駐波相同。在遠(yuǎn)離直墻的方向上,腹、節(jié)點(diǎn)處的波高值變化很大,振幅隨著離堤距離的增加變小,2 個(gè)平均波長(zhǎng)距離后腹、節(jié)點(diǎn)處的波高值變化與平均波高之比小于5%。根據(jù)波高與波浪能量存在H2-E 的關(guān)系,波浪入射直墻一定時(shí)間后,直墻前波能為2E,HS為入射波的倍,即HS/HS0≈1.4,而試驗(yàn)工況離堤距離達(dá)到一定值,波高分布曲線穩(wěn)定后腹點(diǎn)處的比波高值在1.4 左右,這也說(shuō)明了試驗(yàn)的合理性。這表明不同工況的直墻前波高空間分布取決于入射波當(dāng)前水深下的波長(zhǎng),其中譜峰周期TP對(duì)應(yīng)波長(zhǎng)LP與之關(guān)聯(lián)性最強(qiáng)。

        3.1.2 斜向入射波況

        限于篇幅,在此給出平均周期為12 s 時(shí)的不同主波向工況條件下的有效波高分布圖及其腹點(diǎn)包絡(luò)線,各入射波向波況下波高空間分布曲線如圖5 如示。

        圖5 各入射波向波況下波高空間分布曲線Fig.5 Wave height spatial distribution curve under various incident wave directions and wave conditions

        在靠近直墻的部分,主波向?yàn)?°的工況波高整體高于斜向入射的工況,且入射角度越大,比波高整體越低。當(dāng)將波浪分解成橫向、縱向傳遞,入射角度越小,正向傳遞的波浪能量越大,直墻前波浪的反射、疊加作用越大。但在遠(yuǎn)離直墻一定距離之后,由于斜向入射條件下的波高空間分布曲線整體相位被拉長(zhǎng),腹、節(jié)點(diǎn)的分布并不均勻,這是因?yàn)樵谶h(yuǎn)離直墻的位置,由于海浪的多向性,疊加作用強(qiáng)于主波向正向的工況,主波向角度大的曲線不一定低于主波向角度小的曲線,靠近直墻的腹點(diǎn)不一定比遠(yuǎn)離直墻的腹點(diǎn)高。直墻上的比波高值普遍小于2。由圖5 可知,不論何種波況下,直墻前平面波高最低點(diǎn)為離墻最近的節(jié)點(diǎn)處,且該點(diǎn)波高比之后的節(jié)點(diǎn)處波高值小很多。

        3.2 反射后腹點(diǎn)、節(jié)點(diǎn)間距Lf 與波長(zhǎng)L 的關(guān)系

        分別對(duì)腹點(diǎn)和節(jié)點(diǎn)的間距統(tǒng)計(jì)分析,Lf與L的關(guān)系尚不明確,下面通過(guò)彌散方程以及試驗(yàn)數(shù)據(jù)探討兩者之間的聯(lián)系。首先確定Lf如何統(tǒng)計(jì)分析,Lf的測(cè)量方式有兩種,一種是腹點(diǎn)的間距Lf1,第二種是節(jié)點(diǎn)的間距Lf2。將各種工況的波高空間分布曲線中,根據(jù)上述對(duì)不規(guī)則波入射直墻后的腹點(diǎn)、節(jié)點(diǎn)的判斷方式對(duì)腹點(diǎn)、節(jié)點(diǎn)的位置和波高值進(jìn)行統(tǒng)計(jì),相鄰的腹點(diǎn)、節(jié)點(diǎn)的間距即為L(zhǎng)f的計(jì)算方法。將各工況的Lf1和Lf2統(tǒng)計(jì),并引入公式計(jì)算其離散度σ。其表達(dá)式如下:

        由表2、表3 可知,腹點(diǎn)和節(jié)點(diǎn)在主波向?yàn)?°時(shí)分布均勻,它們的間距可以整體被考慮,在主波向角度變大時(shí),分布不再均勻且整體并無(wú)明顯規(guī)律。

        表2 HS 統(tǒng)計(jì)中各波況Lf1 的離散度Table 2 The dispersion degree of each wave condition Lf1 in HS statistics

        表3 HS 統(tǒng)計(jì)中各波況Lf2 的離散度Table 3 The dispersion degree of each wave condition Lf2 in HS statistics

        根據(jù)3.1 節(jié)中對(duì)主波向正向入射的波況,波高空間分布曲線與當(dāng)前水深、譜峰周期對(duì)應(yīng)波長(zhǎng)LP,再加上入射角度、疊加作用等因素影響,直墻前Lf、T 和L 的關(guān)系有待進(jìn)一步確定。受限于篇幅,下面只列出了與Lf值對(duì)應(yīng)關(guān)系最強(qiáng)的LP的計(jì)算結(jié)果和結(jié)果分析。

        3.2.1 正向入射波況

        表4 直墻前第1 個(gè)節(jié)點(diǎn)和其相鄰腹點(diǎn)中點(diǎn)比值Table 4 The ratio of the midpoint between the first node in front of the straight wall and its adjacent antinode

        表5 正向入射時(shí)各波況下的Lf /(1 2Lp)Table 5 LLp)under each wave condition at normal incidence f /(1 2主波向θmain/(°)平均周期/s 8 1.02 1.01 1.01 0 10 12

        3.2.2 斜向入射波況

        隨著主波向角度的變大,Lf值的離散度增大,波高分布曲線變得越來(lái)越不規(guī)則,很難得到準(zhǔn)確的與L 的關(guān)系式,難以計(jì)算出腹點(diǎn)或節(jié)點(diǎn)的準(zhǔn)確位置。

        由表4 可知,對(duì)于不同入射方向,第1 個(gè)節(jié)點(diǎn)與相鄰腹點(diǎn)的中點(diǎn)吻合。為了討論方便,將直墻前第1 段腹點(diǎn)間距定義為L(zhǎng)1f。即在L1f中,腹點(diǎn)和節(jié)點(diǎn)的對(duì)應(yīng)關(guān)系不變。因此可以通過(guò)入射波況確定L1f對(duì)應(yīng)腹點(diǎn)、節(jié)點(diǎn)的位置,可見對(duì)直墻前第1 段腹點(diǎn)間距和波長(zhǎng)的討論是有意義的。下面列出L1f與Lp以及不同入射角θmain的對(duì)應(yīng)關(guān)系。

        在其他入射波要素相同的情況下,隨著主波向角度θmain的變大,L1f變大。

        表6 各波況下的L1f/(1 2cos θmain Lp)Table 6 L1f/(1 2cos θmain Lp)under various wave conditions主波向θmain/(°)平均周期/s 8 10 12 1.02 1.01 1.01 11.25 0.97 1.01 1.00 22.5 0.98 1.00 1.03 33.75 0.99 0.98 1.00 45 0.98 1.00 1.02 0

        3.3 直墻前HS 橫向切面的空間分布情況

        在波浪正向入射的試驗(yàn)中,直墻前橫切面上的等波高線大多與直墻平行,由于入射波是多向不規(guī)則波,等波高線不嚴(yán)格平行十分正常;在波浪斜向入射的試驗(yàn)中,直墻前的等波高線與橫切面并不平行,波高分布呈現(xiàn)相位差。

        4 結(jié)語(yǔ)

        1)多向不規(guī)則波入射后,離岸方向有效波高空間分布情況如下:主波向正向時(shí),不同周期的條件下,堤前相同相對(duì)距離s/LP內(nèi)有效波高空間分布曲線一致,腹、節(jié)點(diǎn)分布規(guī)律與規(guī)則波入射直墻的波況一致,在遠(yuǎn)離直墻方向腹點(diǎn)對(duì)應(yīng)的比波高值逐漸變小,節(jié)點(diǎn)對(duì)應(yīng)比波高值逐漸變大;斜向入射波況中,隨著主波向偏離法向入射的角度變大,波高分布曲線不規(guī)則程度變大,腹、節(jié)點(diǎn)間距變大且不均勻;各工況在距離直墻一定距離后波高值趨于穩(wěn)定,比波高值在1.4 左右;整體波高分布最小值處為堤前第1 個(gè)節(jié)點(diǎn),且該點(diǎn)波高遠(yuǎn)小于其他節(jié)點(diǎn)處波高。

        2)多向不規(guī)則波入射直墻后,腹、節(jié)點(diǎn)間距與譜峰周期計(jì)算得到波長(zhǎng)LP有關(guān)。正向入射時(shí),腹點(diǎn)間距Lf1相較與節(jié)點(diǎn)間距Lf2離散度較小,腹點(diǎn)間距Lf1與LP關(guān)系式為斜向入射時(shí),堤前第1 段L1f與 LP以及主波向角度θmain對(duì)應(yīng),關(guān)系式為以上關(guān)系式可以依據(jù)入射波確定直墻前兩個(gè)腹點(diǎn)的位置和間距以及第1 個(gè)節(jié)點(diǎn)的位置,在工程實(shí)踐中可以考慮將泊位或浮碼頭的位置設(shè)置在第1 個(gè)節(jié)點(diǎn)位置及其位置附近。

        3)波浪正向入射的試驗(yàn)中,直墻前橫切面的波高分布基本為等波高線;波浪斜向入射的試驗(yàn)中,直墻前的等波高線基本不平行于橫切面,波高分布出現(xiàn)相位差。

        猜你喜歡
        直墻入射波波高
        基于FHDI-GNWM 數(shù)據(jù)的全球超越概率波高宏觀分布特征分析
        “直墻平底”之于碑刻·篆刻·書法的藝術(shù)價(jià)值
        SHPB入射波相似律與整形技術(shù)的試驗(yàn)與數(shù)值研究
        直墻半圓拱巷道圍巖應(yīng)力分布解析
        基于漂流浮標(biāo)的南大洋衛(wèi)星高度計(jì)有效波高研究
        深埋直墻拱形隧道穩(wěn)定性研究
        非平整港池的多向不規(guī)則波試驗(yàn)研究
        直墻半圓拱可縮性U型鋼支架卡纜臨界約束力理論研究
        瞬態(tài)激勵(lì)狀態(tài)下樁身速度以及樁身內(nèi)力計(jì)算
        對(duì)機(jī)械波半波損失現(xiàn)象的物理解釋
        電子科技(2015年11期)2015-03-06 01:32:24
        国产成人啪精品视频免费网| 久久久受www免费人成| 国产精自产拍久久久久久蜜| 久久中国国产Av秘 入口| 国产成人av一区二区三| 亚洲欧美中文日韩在线v日本| 男女裸交无遮挡啪啪激情试看| 国产在线不卡AV观看| av天堂一区二区三区精品| 色吧噜噜一区二区三区| 国内露脸少妇精品视频| 亚洲AV永久无码精品导航| 久久精品国产亚洲av一| 大地资源网在线观看免费官网 | 人妻无码久久一区二区三区免费| 免费一级a毛片在线播出| 久久九九精品国产不卡一区| 精品国产乱码久久久久久婷婷| 日本亚洲国产一区二区三区| 巨臀精品无码AV在线播放| 国产不卡视频在线观看| 日本亚洲欧美色视频在线播放| 999国产精品视频| 日韩精品一区二区三区免费观影| 国产精品国产三级国产aⅴ下载| 51久久国产露脸精品国产| 韩国日本亚洲精品视频| 中文字幕亚洲五月综合婷久狠狠| 男人扒开添女人下部免费视频| 一本一本久久a久久精品综合| 国产高跟丝袜在线诱惑| 亚洲男人天堂一区二区| 日韩成人大屁股内射喷水 | 亚洲欧美日韩综合久久| 无码伊人久久大香线蕉| 媚药丝袜美女高清一二区| 99精品国产一区二区三区a片 | 日本中文字幕有码网站| 7777奇米四色成人眼影| 精品免费久久久久国产一区| 高潮内射主播自拍一区|