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

淺水中質(zhì)量源造波方法

2017-09-22 09:56:04田正林孫昭晨梁書(shū)秀
水道港口 2017年4期
關(guān)鍵詞:質(zhì)量

田正林,孫昭晨,梁書(shū)秀

(大連理工大學(xué) 海岸及近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,大連 116024)

淺水中質(zhì)量源造波方法

田正林,孫昭晨,梁書(shū)秀

(大連理工大學(xué) 海岸及近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,大連 116024)

針對(duì)質(zhì)量源造波過(guò)程中源區(qū)上方水體高度難以確定的問(wèn)題提出了改進(jìn)方法,將源區(qū)高度值考慮到源強(qiáng)方程中,避免了造波過(guò)程中源區(qū)上方水體高度值多次試算帶來(lái)的不便。該方法彌補(bǔ)了質(zhì)量源造波方法在水深小于1.5 m的淺水中所造波浪與理論波浪之間差距較大的缺陷,尤其是水平網(wǎng)格步長(zhǎng)與源區(qū)高度值相比小于10的情況,推廣了質(zhì)量源造波方法在淺水中的應(yīng)用。文中驗(yàn)證了波高對(duì)所造波浪的影響以及時(shí)間步長(zhǎng)對(duì)淺水中所造波浪的影響,得出結(jié)論:當(dāng)波高與水深的比值小于0.2時(shí)所造波浪波高與理論值的擬合較好;當(dāng)時(shí)間步長(zhǎng)小于0.005 s時(shí)所造波浪無(wú)明顯的上下漂浮現(xiàn)象,波形穩(wěn)定,符合理論波形要求。

質(zhì)量源造波;淺水造波;源區(qū)高度;改進(jìn)方法

數(shù)值波浪水槽是港口海岸工程數(shù)值模擬的基礎(chǔ),造波系統(tǒng)與消波系統(tǒng)是數(shù)值波浪水槽的重要組成部分,造波系統(tǒng)造波能力及消波系統(tǒng)消波能力又是數(shù)值波浪水槽可用性的關(guān)鍵。大多數(shù)值模擬采用加長(zhǎng)計(jì)算域或是波動(dòng)沒(méi)有達(dá)到準(zhǔn)穩(wěn)定狀態(tài)時(shí)結(jié)束模擬避免二次反射波的影響[1]。為了得到長(zhǎng)時(shí)間穩(wěn)定的數(shù)值模擬,對(duì)數(shù)值水槽的要求是一能夠降低二次反射波,二所造目標(biāo)波浪不受反射波的影響,流體域內(nèi)部造波方法能夠很好地解決該問(wèn)題。Larsen和Dancy[2]在水域內(nèi)部沿一定長(zhǎng)度的直線段上通過(guò)擾動(dòng)水面生成目標(biāo)波浪,避免了二次反射波的影響。

目前,流體內(nèi)域造波主要有動(dòng)量源造波和質(zhì)量源造波。動(dòng)量源造波是將目標(biāo)波浪作為人工源項(xiàng)加到動(dòng)量方程中生成所需波浪場(chǎng),質(zhì)量方程保持不變,自由水面以下區(qū)域均屬造波源區(qū),造波過(guò)程中在整個(gè)造波源區(qū)內(nèi)施加動(dòng)量源方程;質(zhì)量源造波是將目標(biāo)波浪作為人工源項(xiàng)加到質(zhì)量方程中生成所需波浪場(chǎng),造波過(guò)程中并不是自由水面以下區(qū)域均屬造波源區(qū),源區(qū)頂部與自由水面之間存在一定高度無(wú)質(zhì)量源的水體,正因?yàn)榇瞬糠炙w的存在使得源區(qū)高度的確定變得繁瑣。Lin和Liu[3]通過(guò)試算得出質(zhì)量源造波源區(qū)中心到靜水面之間的距離為13~12水深時(shí)所造波浪比較理想。周玲玲,丁全林,蘭慶琳等[4]指出質(zhì)量源距離水面13水深時(shí)的模擬結(jié)果較為精確,并將此處定義為質(zhì)量源的標(biāo)準(zhǔn)位置。余勇飛[5]在深水?dāng)?shù)值水槽造波時(shí)質(zhì)量源區(qū)上方水體高度采用一倍波高時(shí)所造波浪精度較高,與理論波浪匹配較好。李宏偉[6]在造孤立波時(shí)質(zhì)量源區(qū)上方水體高度采用半倍波高時(shí),所造孤立波波形與理論波形匹配較好。Chen和Hsiao[7]在FLOW-3D中開(kāi)發(fā)了水體內(nèi)域造波質(zhì)量源區(qū)域的設(shè)計(jì)方法,詳細(xì)描述了質(zhì)量源位置的確定公式,為質(zhì)量源造波方法的應(yīng)用提供了便利。

已有文獻(xiàn)質(zhì)量源區(qū)上方水體高度大多采用試算法確定,只給出經(jīng)驗(yàn)高度區(qū)間,并且試算法的成功快慢取決于個(gè)人經(jīng)驗(yàn),對(duì)于經(jīng)驗(yàn)不足者尤顯耗時(shí)。筆者造波過(guò)程中發(fā)現(xiàn)質(zhì)量源區(qū)上方水體高度為一倍波高時(shí)在水深大于1.5 m的深水情況下所造波浪極好,但在水深小于1.5 m的淺水中所造波浪并不理想,尤其是在水平網(wǎng)格步長(zhǎng)與源區(qū)高度值相比小于10時(shí)所造波浪波高與目標(biāo)波浪波高相差較大,必須不斷地加大或減小源區(qū)上方水體高度以減小差距。而通常情況下需要造多種波高的目標(biāo)波浪,這樣每次變化波高時(shí)都需要不斷調(diào)整源區(qū)上方水體高度,既費(fèi)時(shí)也是沒(méi)有必要的。故本文針對(duì)質(zhì)量源造波源區(qū)上方水體高度難以確定的問(wèn)題,提出改進(jìn)方法,避免了多次試算帶來(lái)的麻煩,進(jìn)一步完善了質(zhì)量源造波方法。

1 數(shù)值模型

基于FLUENT求解器,求解雷諾平均N-S方程,采用雷諾應(yīng)力粘性模型封閉方程組。采用有限體積法離散方程,PISO壓力速度耦合算法求解離散方程,對(duì)流項(xiàng)采用三階MUSCL格式,用改進(jìn)的QUICK方法跟蹤自由表面。邊界條件為左右側(cè)為對(duì)稱邊界條件,底部為墻邊界條件,上部為壓力進(jìn)口邊界條件。以下僅就文中的控制方程、造波方法及消波方法進(jìn)行描述。

1.1控制方程

張量形式的連續(xù)性方程和動(dòng)量方程如下

(1)

(2)

對(duì)于二維問(wèn)題i,j=1,2;ui,uj為速度;p為壓強(qiáng);ρ為密度;g為重力加速度。粘性應(yīng)力張量τij為

(3)

μ為分子運(yùn)動(dòng)粘性。

1.2造波方法

數(shù)值波浪水槽的質(zhì)量源設(shè)置于計(jì)算域內(nèi),質(zhì)量源區(qū)距離左側(cè)消波區(qū)至少一倍波長(zhǎng)的距離可以很好地消除波浪二次反射(Lin & Liu[1]),本文造波源區(qū)距離左側(cè)消波區(qū)2L(L為波長(zhǎng))。

造波源區(qū)內(nèi)質(zhì)量守恒方程為

(4)

s(xi,t)為質(zhì)量源強(qiáng)方程,其大小由設(shè)計(jì)的入射波波面方程決定,造波源區(qū)內(nèi)非零,造波源區(qū)外為零。理論上,當(dāng)質(zhì)量源寬度遠(yuǎn)小于波長(zhǎng)時(shí),可作點(diǎn)源處理,假定質(zhì)量源引起的質(zhì)量增減對(duì)應(yīng)于入射波波面的變化,并且源區(qū)內(nèi)質(zhì)點(diǎn)的水平速度等于設(shè)計(jì)波浪的水平速度,因?yàn)橘|(zhì)量源可以看做均勻分布的點(diǎn)源,所以在造波源域Ω內(nèi)每個(gè)網(wǎng)格中應(yīng)滿足下列關(guān)系[5-6,8-9]

s(xi,t)ΔxΔyΔt=2u(xi,t)ΔyΔt

(5)

化簡(jiǎn)式(5)得下式

(6)

式中:xi為“源強(qiáng)”的水平位置;Δx為網(wǎng)格的水平尺度;Δy為網(wǎng)格的垂向尺度;右端系數(shù)2表示質(zhì)量源產(chǎn)生傳播方向相反的兩列波。

二維二階Stocks 波的水平速度方程如下

(7)

將式(7)代入(6)得源強(qiáng)方程如下

(8)

式中:H為波高;k為波數(shù);ω為頻率;D為水深。

以下通過(guò)對(duì)造波源區(qū)域源強(qiáng)方程的改進(jìn),可以得到更適合淺水區(qū)域的造波方法。造波源區(qū)簡(jiǎn)圖如圖1所示,造波源區(qū)底部為墻邊界并無(wú)水體向下噴出,考慮單位時(shí)間內(nèi)造波源區(qū)其他3個(gè)方向噴出的水體質(zhì)量守恒可推得

s(t)·Δx·Hs=2u·Hs+Δx·v

(9)

式中:u為造波源區(qū)流體質(zhì)點(diǎn)水平速度;v為造波源區(qū)流體質(zhì)點(diǎn)垂向速度;Hs為造波源區(qū)高度?;?jiǎn)(9)可得

s(t)=2uΔx+vHs

(10)

圖1 源造波區(qū)簡(jiǎn)圖Fig.1 Sketch of the source wave-maker zone

質(zhì)量源造波假定源區(qū)內(nèi)各點(diǎn)產(chǎn)生的水平速度等于波生成對(duì)應(yīng)點(diǎn)產(chǎn)生的水平速度,造波源區(qū)上方水體依靠水體自重形成的恢復(fù)力產(chǎn)生波浪,因此源區(qū)上方至自由水面高度的確定比較關(guān)鍵。仔細(xì)研究方程(10)可以發(fā)現(xiàn):造波過(guò)程中u和v的量級(jí)相當(dāng),由于Δx一般小于波長(zhǎng)的5%,非常小,而Hs一般大于23倍的水深,對(duì)于水深大于1.5 m的深水Δx相比Hs更為小量,至少小一個(gè)量級(jí),式中的第二項(xiàng)對(duì)源強(qiáng)的貢獻(xiàn)較小可以忽略,但對(duì)于水深小于1.5 m的淺水,第二項(xiàng)會(huì)隨著水深的減小而增大,其對(duì)造波源的貢獻(xiàn)也會(huì)逐漸增大,如果第二項(xiàng)繼續(xù)忽略會(huì)影響所造波浪的波高。數(shù)值波浪水槽造波時(shí)如果水深小于1.5 m,Hs小于1.0 m,Δx又小于0.1 m,可以看出此時(shí)第一項(xiàng)與第二項(xiàng)相比在10以內(nèi),不足一個(gè)量級(jí),如果水深繼續(xù)減小,Hs繼續(xù)減小,式中第二項(xiàng)對(duì)源強(qiáng)的貢獻(xiàn)會(huì)更大,不應(yīng)忽略。根據(jù)Broren.M.,Larsen.J.[10]描述,源造波的原理類似于“水下爆炸”原理,那么源爆炸瞬間質(zhì)點(diǎn)各個(gè)方向的速度大小應(yīng)該相等。所以質(zhì)量源向外噴出的流體速度各方向應(yīng)該大小相等,因此u=v,最終可得源強(qiáng)方程為

(11)

周玲玲,丁全林,蘭慶琳等[4]指出為避免數(shù)值振蕩,質(zhì)量源初始時(shí)刻應(yīng)設(shè)置為0,二階Stocks波令初始相位為-π2,并在計(jì)算歷程開(kāi)始時(shí)刻設(shè)置軟啟動(dòng)緩慢施加質(zhì)量源驅(qū)動(dòng)波面變化。本文采用此種做法,最終施加的源強(qiáng)方程如下

圖2 源區(qū)右側(cè)數(shù)值水槽示意圖Fig.2 Right side schematic diagram of numerical wave tank

(12)

式中:N為緩變控制參數(shù);t0為緩變時(shí)間。

1.3消浪層

為了吸收二次反射波和不需要的源造波,數(shù)值波浪水槽左右端應(yīng)布置至少一倍波長(zhǎng)的海綿吸波層耗散波能。本文采用的吸收方程是由Niels and David[11]提出的橢圓方程(13),x0表示消波層起始點(diǎn)x坐標(biāo);xs表示消波層長(zhǎng)度,方程從海綿層起點(diǎn)1逐漸減小到海綿層終點(diǎn)0。

圖3 實(shí)際波形與理論波形的對(duì)比Fig.3 Comparisons of present wave elevation and theory wave elevation

海綿層消波方程μ(X)

(13)

數(shù)值波浪水槽源區(qū)右半側(cè)示意圖,如圖2所示。

2 數(shù)值模擬結(jié)果

本文數(shù)值波浪水槽長(zhǎng)36 m,高0.7 m,水深0.5 m,質(zhì)量源高度Hs為水深的23,△x小于波長(zhǎng)的5%。以下對(duì)所造波浪結(jié)果、波形沿程變化、波高對(duì)目標(biāo)波浪的影響以及時(shí)間步長(zhǎng)對(duì)波形的影響進(jìn)行描述。

圖4 波形沿程變化圖Fig.4 Wave elevation displacement of present model with distance from the source zone

2.1實(shí)際波形與理論波形對(duì)比

為了驗(yàn)證所造波浪的有效性,給出水深D=0.5 m,波高H=0.06 m,周期T分別為0.8 s、1.0 s、1.2 s、1.4 s時(shí)的實(shí)造波浪波形與理論波浪波形的對(duì)比(圖3)。從圖3中可以看出:不同周期時(shí)的實(shí)際波形與理論波浪波形匹配較好,雖波谷處振幅值比理論值略小但仍都可以滿足精度要求。

2.2波浪沿程變化情況

為了說(shuō)明波浪沿程變化情況,給出了水深D=0.5 m,波高H=0.06 m,周期T分別為0.8 s、1.0 s、1.2 s、1.4 s時(shí)造波源區(qū)(約5倍波長(zhǎng)處)右側(cè)波浪10倍波長(zhǎng)范圍內(nèi)的波浪沿程變化情況(圖4)。從圖4可以看出:所造波浪沿程衰減較小,滿足目標(biāo)波浪精度要求。

圖5 所造波浪波高與理論波高的對(duì)比Fig.5 Comparisons of the targeted wave height and theoretical wave height

2.3波高對(duì)目標(biāo)波浪的影響

為了檢驗(yàn)改進(jìn)的質(zhì)量源造波方法所造不同波高的波浪與理論波浪的差距,本部分計(jì)算了D=0.5 m,T=1.2 s,波高H分別為0.06 m、0.08 m、0.10 m、0.12 m造波源區(qū)(約5倍波長(zhǎng)處)右側(cè)波浪10倍波長(zhǎng)范圍內(nèi)的波浪隨時(shí)間的變化(圖5)。從圖5可以看出:隨著波高的增大,所造波浪波高與理論波高的差距在增大;波高H=0.10 m時(shí)所造波浪波高只是在波峰處略小于理論值,波谷處仍與理論值擬合較好,但是當(dāng)波高H=0.12 m時(shí)所造波浪波峰和波谷處均與理論值有偏差,所造波浪效果變差。由此可知當(dāng)波高與水深的比值超過(guò)0.2時(shí)所造波浪波高與理論值的擬合較差,所以此種造波方法適合于波高與水深的比值在0.2以下的情況。

2.4時(shí)間步長(zhǎng)的影響

圖6 波形隨時(shí)間變化圖Fig.6 Time step effects on wave surface

為了檢驗(yàn)時(shí)間步長(zhǎng)對(duì)改進(jìn)的質(zhì)量源造波方法所造波浪的影響,給出了水深D=0.5 m,周期T=1.2 s,波高H=0.06 m時(shí)所造波浪波形隨時(shí)間的變化。時(shí)間步長(zhǎng)驗(yàn)證之前需要對(duì)網(wǎng)格進(jìn)行反復(fù)調(diào)整驗(yàn)證,當(dāng)網(wǎng)格連續(xù)3次調(diào)整后的波形無(wú)較大變化后停止,進(jìn)行時(shí)間步長(zhǎng)驗(yàn)證。時(shí)間步長(zhǎng)分別取0.01 s、0.005 s、0.002 s、0.000 2 s進(jìn)行驗(yàn)證。對(duì)比所得波面由圖6可以看出:當(dāng)dt=0.01 s、dt=0.005 s時(shí),隨著時(shí)間的增長(zhǎng)波面出現(xiàn)了幅度較大的上下飄動(dòng)現(xiàn)象,而當(dāng)dt=0.002 s、dt=0.000 2 s時(shí),波面已無(wú)明顯的上下飄動(dòng)現(xiàn)象,波面較穩(wěn)定。

3 結(jié)論

將質(zhì)量源區(qū)高度值考慮到造波源強(qiáng)方程之中,量化源區(qū)上方水體高度,克服了水深小于1.5 m時(shí)造波過(guò)程中源區(qū)上方水體高度難以確定的困難,得到了適合于淺水的質(zhì)量源造波方法。該方法可以在水深小于1.5 m的淺水中造出比較理想的波浪,更適合于水平網(wǎng)格步長(zhǎng)相比于源區(qū)高度值小10倍的情況。文中驗(yàn)證了波高對(duì)所造波浪的影響、時(shí)間步長(zhǎng)對(duì)淺水中所造波浪的影響及波浪沿程變化,得出結(jié)論:當(dāng)波高與水深的比值小于0.2時(shí)所造波浪波高與理論波高值的擬合較好;當(dāng)時(shí)間步長(zhǎng)小于0.005 s時(shí)所造波浪無(wú)明顯的上下漂浮現(xiàn)象,波形穩(wěn)定,符合理論波形要求,波浪沿程衰減較小。

[1]Lin P, Liu L F. Discussion of "Vertical variation of the flow across the surf zone" [Coast. Eng. 45 (2002) 169-198][J]. Coastal Engineering, 2004, 50(3):161-164.

[2]Larsen J, Dancy H. Open boundaries in short wave simulations-A new approach[J]. Coastal Engineering, 1983, 7(3):285-297.

[3]Lin P, Liu L F. Internal Wave-maker for Navier-stokes Equations Models[J].Journal of Waterway, Port, Coastal and Ocean Engineering, 1999,125(4): 207-215.

[4]周玲玲,丁全林,蘭慶琳,等.基于LB方法的質(zhì)量源造波的數(shù)值波浪水槽[J].水電能源科學(xué),2015,33(3):99-103. ZHOU L L,DING Q L,LAN Q L, et al.Numerical Wave Tank with Mass Source Generating Waves Based on LB Method[J]. China Civil Engineering Journal, 2015, 33(3): 99-103.

[5]余勇飛.波浪的數(shù)值模擬及其流場(chǎng)結(jié)構(gòu)分析[D].哈爾濱:哈爾濱工業(yè)大學(xué),2013.

[6]李宏偉.造波理論與方法研究[D].哈爾濱:哈爾濱工程大學(xué),2013.

[7]Chen Y L, Hsiao S C. Generation of 3D water waves using mass source wavemaker applied to Navier-Stokes model[J]. Coastal Engineering, 2016, 109:76-95.

[8]房卓,張寧川,臧志鵬.透空式梳式防波堤的數(shù)值模擬和波浪透射系數(shù)的研究[J].水道港口,2011,32(2):86-93. FANG Z, ZHANG N C, ZANG Z P. Numerical simulations of open comb-type breakwater and research on its wave transmission coefficient[J].Journal of Waterway and Harbor,2011,32(2):86-93.

[9]FANG Z, ZHANG N C, ZANG Z P. Experimental and Numerical Study on Hydrodynamic Performance of Impermeable Comb-Type Breakwater[J]. Journal of Ship Mechanics, 2012, 16(6):632-645.

[10]Brorsen M, Larsen J. Source generation of nonlinear gravity waves with the boundary integral equation method[J]. Coastal Engineering, 1987, 11(2):93-113.

[11]Jacobsen N G, Fuhrman D R, Fredsφe J. A wave generation toolbox for the open-source CFD library: OpenFoam[J]. International Journal for Numerical Methods in Fluids, 2012, 70(9):1 073-1 088.

The method of mass source wavemaker in shallow water

TIANZheng-lin,SUNZhao-chen,LIANGShu-xiu

(StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,China)

In view of the problem that the distance between the top of the source region and the still water level is difficult to determine in the process of mass source wavemaker, the improved method was proposed in this paper. In the mass source functions took the height of the source region into account, the inconvenience of many trials to get the height of the mass source region can be eliminated in the process of wavemaker. This method remedies the defect that the difference between targeted wave and theoretical wave in the shallow water whose depth less than 1.5 m, especially horizontal grid step is an order of magnitude 10 smaller than the source region height, the mass source wave-maker method is promoted in the shallow water.The effects of wave height and the time step on simulated waves were studied. The conclusions show that the fitting between made wave and theoretical wave is good when the ratio of wave height and water depth is less than 0.2;when the time step is less than 0.005 s, the floating phenomenon of the simulated waves is well controlled, waves surface is stable and conforms to the requirements of the wave theory.

mass source wavemaker; shallow-water wavemaker; height of source region; improved method

TV 139.2

:A

:1005-8443(2017)04-0325-05

2017-01-13;

:2017-04-06

國(guó)家自然科學(xué)基金資助項(xiàng)目(51279028);國(guó)家海洋局公益性行業(yè)科研專項(xiàng)(201405025-1)

田正林(1985-),男,黑龍江省哈爾濱人,博士研究生,主要從事波浪與結(jié)構(gòu)物相互作用的研究。

Biography: TIAN Zheng-lin(1985-),male, doctor student.

猜你喜歡
質(zhì)量
聚焦質(zhì)量守恒定律
“質(zhì)量”知識(shí)鞏固
“質(zhì)量”知識(shí)鞏固
質(zhì)量守恒定律考什么
做夢(mèng)導(dǎo)致睡眠質(zhì)量差嗎
焊接質(zhì)量的控制
關(guān)于質(zhì)量的快速Q(mào)&A
初中『質(zhì)量』點(diǎn)擊
質(zhì)量投訴超六成
你睡得香嗎?
民生周刊(2014年7期)2014-03-28 01:30:54
主站蜘蛛池模板: 思思99热精品在线| 国产精品久久精品| 91久久偷偷做嫩草影院电| 国产精品青青| 亚洲天堂精品在线| 国产精品久久久久久久伊一| 中文字幕66页| 91久久国产综合精品女同我| 国产美女精品一区二区| 精品国产一区二区三区在线观看| 日韩中文精品亚洲第三区| 久久久91人妻无码精品蜜桃HD | 日韩精品一区二区深田咏美| 免费国产不卡午夜福在线观看| yy6080理论大片一级久久| 国产全黄a一级毛片| 99精品一区二区免费视频| 都市激情亚洲综合久久| 99在线观看免费视频| 老司国产精品视频| 亚洲动漫h| 国产精品午夜电影| 无码AV高清毛片中国一级毛片| 亚洲天堂日韩av电影| 91人妻在线视频| 夜夜拍夜夜爽| 日韩精品成人网页视频在线| 久久亚洲国产最新网站| 国产大片黄在线观看| 午夜视频免费试看| 福利在线一区| 久久青草视频| 成人免费一级片| 国产av无码日韩av无码网站 | 国产Av无码精品色午夜| 国产高清无码麻豆精品| a在线亚洲男人的天堂试看| 精品国产免费观看一区| 一本色道久久88亚洲综合| 欧美一级高清视频在线播放| 蜜桃视频一区二区三区| 久草性视频| 亚洲天堂首页| 无码电影在线观看| 欧美在线伊人| 99久久精品久久久久久婷婷| 九九热精品视频在线| 色婷婷色丁香| 久久午夜影院| 91破解版在线亚洲| 国产在线自乱拍播放| 国模视频一区二区| 国产玖玖视频| 91小视频在线观看免费版高清| 久久精品娱乐亚洲领先| 亚洲中文字幕无码mv| 国产女人在线观看| 最新亚洲人成无码网站欣赏网| 秋霞国产在线| 国产福利免费在线观看| 97超爽成人免费视频在线播放| 国产91透明丝袜美腿在线| 青青青视频免费一区二区| 伊人蕉久影院| 女同久久精品国产99国| 国产一区成人| 新SSS无码手机在线观看| 国产偷国产偷在线高清| 国产亚洲视频免费播放| 色综合天天娱乐综合网| 国产欧美日韩91| 日韩亚洲综合在线| 免费欧美一级| 国产精品制服| 伊人精品视频免费在线| 婷婷色一区二区三区| 久久久久无码国产精品不卡| 毛片网站免费在线观看| 亚洲一区网站| 免费无码又爽又黄又刺激网站 | 亚洲浓毛av| 亚洲中文字幕97久久精品少妇|