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

畸形波生成、演化過程研究——波面時間過程隨位置的變化

2011-12-28 04:48:46崔成張寧川裴玉國
海洋通報 2011年4期

崔成,張寧川,裴玉國

(大連理工大學 港口海岸及近海工程國家重點實驗室,遼寧 大連 116024)

畸形波生成、演化過程研究
——波面時間過程隨位置的變化

崔成,張寧川,裴玉國

(大連理工大學 港口海岸及近海工程國家重點實驗室,遼寧 大連 116024)

基于完全非線性波浪數(shù)值模型,模擬實測畸形波的生成、演化過程,進而分析波面時間過程和畸形波參數(shù)的變化規(guī)律。研究結(jié)果表明:在畸形波形成前、后均可能出現(xiàn)連續(xù)的大波和“深谷”現(xiàn)象。在畸形波生成、演化過程中出現(xiàn)的“深谷”和畸形波具有伴生性而非獨立異常波浪現(xiàn)象。上述大波、“深谷”和畸形波的連續(xù)作用對海洋結(jié)構(gòu)物和船舶造成的潛在威脅可能遠超過“瞬態(tài)”畸形波,此點應引起高度重視。

畸形波;生成、演化過程;深谷

已有一定量關(guān)于畸形波的研究和探討,但是關(guān)于畸形波形成、演化過程的研究還相對較少。Osborne[1]等根據(jù) B-F不穩(wěn)定性使用非線性薛定諤方程模擬了畸形波的發(fā)展過程,并指出在巨大波峰前后都會形成深谷。Stansell[2]分析北海實測記錄時,給出了實測記錄中包含深谷的時間序列,其中波谷5.38 m,有效波高3.79 m,波谷與有效波高的比值1.4。Clauss[3,4]指出,在惡劣的海況條件下,不僅發(fā)生獨立的大波還會發(fā)生洶涌的大波群,并通過物理實驗和數(shù)值水槽分別模擬了上述兩種極端波浪。

現(xiàn)有理論研究和實測資料均間接表明,在畸形波形成和演化過程中,可能接連出現(xiàn)波高略小于畸形波的的連續(xù)大波和巨大深谷。由于他們一起對海洋結(jié)構(gòu)物和船舶造成巨大的潛在威脅可能遠超過“瞬態(tài)”畸形波,因此有必要對畸形波生成、演化過程進行進一步深入研究和探討。

本文采用完全非線性數(shù)值模型模擬畸形波生成、演化過程,分析畸形波和連續(xù)大波、深谷之間的伴生關(guān)系。

1 畸形波模擬的非線性數(shù)值模型

將二維不可壓縮流體的連續(xù)方程和完全非線性N-S方程作為畸形波模擬的非線性數(shù)值模型的控制方程,設(shè)定恰當?shù)臄?shù)值水槽邊界條件,基于VOF方法求解,獲得含有畸形波的波動過程。

1.1 控制方程

采用 VOF數(shù)值方法時,不可壓二維流體的連續(xù)方程可記為:

式中 u,v分別表示 x,y方向的速度分量;gx,gy分別表示x,y方向的體積加速度;p為流體壓力;ρ為流體的密度;ν為流體的運動粘滯系數(shù)。

1.2 邊界條件

1.2.1 自由表面條件 波浪自由表面可以根據(jù)體積函數(shù)分布計算得到。體積函數(shù)F滿足下面的控制方程:

1.2.2 數(shù)值水槽造波邊界 在此考慮了兩種情況:①模擬實測畸形波序列(新年波,北海畸形波,北海深谷)的生成、演化過程時,數(shù)值水槽造波邊界采用單波列疊加模型;②人工數(shù)值生成畸形波序列時,數(shù)值水槽造波邊界采用雙波列疊加模型。

單波列疊加模型的波面可描述為:

N是樣本個數(shù),?t為采樣間距,i表示組成波編號i = 0,1,2,……, M;M表示組成波個數(shù)。

雙波列疊加模型的波面可描述為:

根據(jù)造波理論,造波邊界處的水平速度 U0(t)的表達見式(7)

1.2.3 數(shù)值水槽吸收邊界 數(shù)值造波水槽遠端邊界采用的是海綿層邊界。二次反射的吸收采用了主動吸收的方法[7]。

1.2.4 數(shù)值水槽水底邊界 數(shù)值水槽水底為平底,邊界條件為:Uy(i, 1) = 0,Ux(i, 1) = Ux(i, 2)。其中Uy(i, 1)表示水底層網(wǎng)格垂向速度分量,Ux(i, 1)和Ux(i, 2)分別表示水底層及其上一層網(wǎng)格水平速度分量。

2 數(shù)值水槽的有效性驗證

2.1 數(shù)值波浪水槽有效性驗證

水槽長為6 000 cm,劃分1 200×30個網(wǎng)格,網(wǎng)格尺寸5 cm×3 cm,時間步長0.02 s,計算時間150 s,采集點距造波板2 002.5 cm處(滿足大于6倍譜峰周期對應的波長,目標譜選用修正 PM 譜(ITTC,1972)。表達式見式(9)~(11)。

圖1給出了數(shù)值水槽的布置圖。為了驗證本文所建立的數(shù)值模型有效性,使用該模型模擬一組條件為水深d = 50 cm;Hs= 5.00 cm和Tp= 1.50 s的不規(guī)則波列。模擬結(jié)果有效波高和譜峰周期分別為Hs=5.08 cm和Tp= 1.54 s。圖2給出了數(shù)值計算的波面時間過程及其傅里葉能量譜和目標波面及譜的比較,對比發(fā)現(xiàn)兩者計算結(jié)果基本吻合,可見本文所建立的數(shù)值水槽對不規(guī)則波浪的模擬具有一定的有效性。

圖1 數(shù)值水槽示布置(單位:cm)Fig. 1 Sketch of the computation domain (unit: cm)

圖2 數(shù)值計算波面及傅里葉譜與目標波面及傅里葉譜的對比Fig. 2 Comparison of target wave surface elevation and computed wave surface elevation (left) as well as comparison of target spectrum and computed spectrum (right)

再使用該模型模擬在d = 50 cm;Hs= 6 cm和Tp= 2.5 s條件下模擬一組包含畸形波的不規(guī)則波列。圖3給出模擬結(jié)果,該畸形波參數(shù)為:Hs= 6.05 cm,Hj=15.51 cm, ηj=10.36 cm,α1=2.52,α2= 2.08,α3= 3.90,α4= 0.67。其中Hs表示波列有效波高,Hj表示畸形波波高,ηj表示畸形波波峰,α1表示畸形波波高與有效波高的比,α2表示畸形波波高與前一個相鄰波高的比,α3表示畸形波波高與后一個相鄰波高的比,α4表示畸形波波峰與畸形波波高的比。

以上 2算例表明本模型不但可以有效地模擬常規(guī)不規(guī)則波列,還可以模擬包含畸形波的不規(guī)則波列。

3 畸形波生成、演化過程

作為具有代表性的算例,在前面提到的數(shù)值水槽中,定點模擬目標波面為天然條件下實測得到的三組波面時間過程(新年波,北海畸形波,北海深谷[2,5])和 4組人造畸形波,進而分別討論上述 7組不同畸形波波面的生成、演化過程。

在數(shù)值水槽中,指定在距離造波邊界大于6倍有效波長處生成目標波面。計算波面生成、演化過程時,考察了畸形波生成位置前后各 4~6倍有效波長范圍內(nèi),100倍平均周期時段內(nèi)波面變化。波面空間采樣間隔5 cm;時間采樣間隔0. 02 s。

圖3 模擬畸形波Fig. 3 Simulated freak wave

3.1 實測天然畸形波的生成、演化過程

3.1.1 實測新年波的生成、演化過程 表1詳細給出了實測[2]和模擬新年波相關(guān)參數(shù)值的比較。其中的T表示平均周期,L表示平均波長,T2= 2πL / (g th2πd / L),圖4給出實測和模擬的新年波波面時間過程比較。表1和圖4的結(jié)果顯示,新年波模擬結(jié)果與實測值吻合良好。

指定新年波(畸形波)在X=1 992.5 cm(距離造波板7倍有效波長)處生成。圖2分別給出了畸形波生成位置前后各4倍有效波長范圍內(nèi)、8個代表點處、采集的 30倍譜峰周期時段內(nèi)波面時間變化過程。

其中8個代表點的選取原則為該點波面時間過程含有滿足畸形波定義參數(shù)的較大波浪。將波面時間過程進行 T*= (t – tc) / Tp,η*= η / Hs無量綱和時間歸零處理,較大波高出現(xiàn)在T*=0時刻。

圖5(a) ~ (d)分別給出了畸形波生成位置前4倍有效波長范圍內(nèi)、4個代表點處、采集的30倍譜峰周期時段內(nèi)波面時間變化過程。4個代表點為X =1 492.5 cm、1 727.5 cm、1 927.5 cm和1 992.5 cm,距離畸形波的生成位置分別為:-1.7、-0.9、-0.2和0倍有效波長(-表示畸形波生成前)。在該4組序列中,分別出現(xiàn)滿足α1= 2.16的波浪(α1> 2,僅能滿足畸形波定義的第一個參數(shù),該波浪與鄰近的波浪一起表現(xiàn)為連續(xù)的大波),滿足α1= 2.18和α2=2.18 的波浪(α1> 2 和 α2> 2 滿足畸形波定義的第一和第二個參數(shù),該波浪的波谷很大,形成了深谷),滿足 α1= 2.13,α3= 2.73和 α4= 0.69的波浪(α1>2,α3>2和α4>0.65滿足畸形波定義的第一、第三和第四個參數(shù),該波浪表現(xiàn)為波峰值較大的大波)和滿足 α1= 2.21,α2= 2.18,α3= 3.66和 α4= 0.72的畸形波(α1>2,α2>2,α3>2 和 α4>0.65 滿足畸形波定義所有參數(shù))。

上述即為畸形波生成過程,可綜合描述為:畸形波生成前,將首先出現(xiàn)連續(xù)大波(大波群),該波群繼而形成較大的波谷,再生成畸形波。

圖5(e) ~ (h)分別給出了畸形波生成位置后4倍有效波長范圍內(nèi)、4個代表點處、采集的30倍譜峰周期時段內(nèi)波面時間變化過程。4個代表點為X =2 102.5 cm、2 157.5 cm、2 247.5 cm和2 397.5 cm,距離畸形波的生成位置分別為:0.4、0.6、0.9和1.4倍有效波長。在該 4組序列中,分別出現(xiàn)滿足α1=2.26,α2=2.71 和 α3= 2.44 的波浪(α1>2,α2>2和 α3>2滿足畸形波定義的第一、第二和第三個參數(shù),該波浪表現(xiàn)為獨立的大波),滿足 α1=2.11和α2=2.43 的波浪(α1>2和 α2>2滿足畸形波定義的第一和第二個參數(shù),該波浪的波谷很大,形成了深谷),滿足 α1=2.14,α3=2.18 和 α4=0.65 的波浪(α1>2,α3>2和α4>0.65滿足畸形波定義的第一,第三和第四個參數(shù),該波浪表現(xiàn)為波峰值較大的大波)和滿足α1=2.00的波浪(α1>2滿足畸形波定義的第一個參數(shù),該波浪與鄰近的波浪一起表現(xiàn)為連續(xù)的大波)。

上述即為畸形波演化過程,可綜合描述為:畸形波生成后,將首先演化為獨立的大波(不同時滿足畸形波定義的所有參數(shù)),然后出現(xiàn)較大的波谷,再演化成大波群,直至恢復到常規(guī)波浪形態(tài)。

表2給出生成、演化過程中出現(xiàn)的較大波浪的畸形波參數(shù)。本算例的畸形波生成、演化過程可簡記為:連續(xù)大波→深谷→大峰值波→畸形波→獨立大波→深谷→大峰值波→連續(xù)大波。空間范圍約3.1倍有效波長,時間長度約10倍平均周期。

表1 實測和模擬新年波相關(guān)參數(shù)值比較匯總Tab. 1 Coefficient values for the simulated and measured New Year wave records

圖4 實測和模擬的新年波時間過程Fig. 4 Measured and simulated New Year wave

圖 6 實測和模擬北海畸形波時間過程 Fig. 6 Measured and simulated North Sea freak wave

圖5 新年波生成、演化過程中不同代表點波面時間過程Fig. 5 Time-process of water-surface elevation at different representative positions

3.1.2 實測北海畸形波的生成、演化過程 圖 6給出實測和模擬的北海畸形波波面時間過程,表 3匯總給出了實測和模擬的北海畸形波參數(shù)比較。表3和圖6的結(jié)果顯示,北海畸形波模擬結(jié)果與實測值吻合良好。

指定北海畸形波在X = 1 962.5 cm(距離造波板約10倍有效波長)處生成。圖7(a)~(g)分別給出了畸形波生成位置前后各 4倍有效波長范圍內(nèi)、7個代表點處、采集的30倍譜峰周期時段內(nèi)波面時間變化過程。7個代表點為X = 1 663.5 cm、1 748.5 cm、1 930.5 cm、1 962.5 cm、2 037.5 cm、2 107.5 cm、2 231.5 cm,距離畸形波的生成位置分別為:-1.5、-1.1、-0.2、0、0.4、0.7、1.4倍有效波長(-表示畸形波生成前,+表示畸形波生成后)。表4匯總給出了北海畸形波的生成、演化過程中較大波的畸形波參數(shù)。

表2 新年波生成、演化過程中出現(xiàn)的極端波浪的相關(guān)參數(shù)值Tab. 2 Coefficient values for the extreme waves

圖7 北海畸形波生成、演化過程中不同代表點波面時間過程Fig.7 Time-process of water-surface elevation at different representative positions

表3 實測和模擬北海畸形波相關(guān)參數(shù)值比較匯總Tab. 3 Coefficient values for the simulated and measured North Sea freak wave records

表4 北海畸形波生成、演化過程中出現(xiàn)的極端波浪的相關(guān)參數(shù)值Tab. 4 Coefficient values for the extreme waves

綜合比較實測北海畸形波與新年波的生成、演化過程,兩者完全類似(與新年波的生成、演化過程的不同之處為:畸形波演化過程中不包含滿足條件 α1>2,α3>2 和 α4>0.65 的大峰值波),不贅述。空間范圍約3倍有效波長,時間長度約9倍平均周期。

3.1.3 實測北海深谷的生成、演化過程 圖8給出實測和模擬的北海深谷波面時間過程。表5匯總給出了實測和模擬的北海深谷參數(shù)比較。

表5和圖8的結(jié)果顯示,北海深谷模擬結(jié)果與實測值吻合良好。

指定北海深谷在X = 1 710.5 cm(距離造波板約9倍有效波長)處生成。圖9(a)~(g)分別給出了深谷生成位置前后各6倍有效波長范圍內(nèi)、7個代表點處、采集的 30倍譜峰周期時段內(nèi)波面時間變化過程。在X = 2 425.5 cm的位置采集到波面時間過程中包含畸形波。

7個代表點分別為X = 1 410.5 cm、1 710.5 cm、1 895.5 cm、2 425.5 cm、2 765.5 cm、2 835.5 cm、2 925.5 cm,距離畸形波生成位置分別為-4.6、-3.3、-2.4、0、 1.5、1.9、2.3倍有效波長(-表示畸形波生成前,+表示畸形波生成后)。表6匯總給出了與北海深谷相伴發(fā)生的畸形波的生成、演化過程中較大波的畸形波參數(shù)。

圖9 北海深谷伴生畸形波的生成、演化過程中不同代表點波面時間過程Fig. 9 Time-process of water-surface elevation at different representative positions

表5 實測和模擬北海深谷相關(guān)參數(shù)值比較匯總Tab 5 Coefficient values for the simulated and measured North Sea deep trough records

表6 與北海深谷伴生畸形波的生成、演化過程中出現(xiàn)的極端波浪的相關(guān)參數(shù)值Tab 6 Coefficient values for the extreme waves

綜合比較該畸形波與前兩組畸形波的生成、演化過程完全類似(與新年波的生成、演化過程的不同之處為:畸形波演化過程中不包含滿足條件α1>2,α3>2和 α4>0.65大峰值波),不贅述。空間范圍約7倍有效波長,時間長度約13倍平均周期。

上述3個算例綜合比較可以看出:

(1) 在畸形波生成過程中,首先出現(xiàn)不完全滿足畸形波全部定義的較大波浪或連續(xù)的大波;進而出現(xiàn)波谷較大的波浪(通常稱深谷),逐漸發(fā)展會形成滿足全部定義的畸形波(滿足α1>2,α2>2,α3>2 和 α4>0.65)。

(2) 畸形波演化過程中,首先會演化為獨立大波,再演化為深谷,最終分化成連續(xù)大波。

(3)深谷與畸形波生成、演化過程有伴生性,而非獨立出現(xiàn)的異常波浪。

表7匯總給出了伴隨畸形波的生成、演化過程中無量綱深谷。(α5= ηtmax/ Hs,α51表示發(fā)生在畸形波之前無量綱波谷,α52表示發(fā)生在畸形波之后無量綱波谷)。比較可見,發(fā)生在畸形波生成、演化過程中的兩次深谷最大值基本相同。

表7 深谷相關(guān)參數(shù)Tab 7 Coefficient values for the deep troughs

3.2 人造畸形波的生成、演化過程

模擬4組人造畸形波,對應的無量綱水深d / L分別為0.15,0.28,0.65,0.82。采用雙波列疊加模型做為數(shù)值水槽的造波邊界 ,4組人造畸形波造波條件分別為:(1)Hs= 5.2 cm,Tp= 1.8 s,d = 50 cm;(2)Hs= 4.2 cm,Tp= 1.0 s,d = 50 cm;(3)Hs=4.0 cm,Tp= 0.95 s,d = 80 cm;(4)Hs= 4.5 cm,Tp= 1.2 s,d = 150 cm。

在無量綱水深d / L = 0.15 條件下,指定畸形波在X=2515.5 cm(距離造波板7倍有效波長)處生成。圖7(a)~(i)分別給出了畸形波生成位置前后各4倍有效波長范圍內(nèi)、9個代表點處、采集的30倍譜峰周期時段內(nèi)波面時間變化過程。9個代表點分別為X=1 590.5 cm、1 940.5 cm、2 005.5 cm、2 515.5 cm、2 785.5 cm、3 005.5 cm、3 315.5 cm、3 510.5 cm和3 810.5 cm距離畸形波生成位置分別為-2.6、-1.6、-1.4、0、+0.8、+1.4、+2.2、+2.8和+3.6倍有效波長(-表示畸形波生成前,+表示畸形波生成后)。表8匯總給出了無量綱水深d / L = 0.15條件下人造畸形波的生成、演化過程中較大波的畸形波參數(shù)。

圖10 d·L-1 = 0.15 條件下畸形波生成、演化過程中不同代表點波面時間過程Fig. 10 Time-process of water-surface elevation at different representative positions

表8 d·L-1 = 0.15 條件下畸形波生成、演化過程中出現(xiàn)的極端波浪的相關(guān)參數(shù)值Tab 8 Coefficient values for the extreme waves

表9 深谷相關(guān)參數(shù)值Tab 9 Coefficient values for the deep troughs

該畸形波的生成、演化過程與前面3組畸形波的生成、演化過程完全類似(與新年波的生成、演化過程的不同之處為:畸形波演化過程包含滿足條件α1>2和α3>2的大波),不贅述,空間范圍約6.2倍有效波長,時間長度約 10倍平均周期。無量綱水深d / L 分別為0.28,0.65,0.82條件下模擬的畸形波生成、演化過程與前面結(jié)果完全一致,不贅述,空間范圍分別為6.2、3.4和4.6倍有效波長,時間長度分別為:11、11和13倍平均周期。表9給出4組人造畸形波生成、演化過程中無量綱深谷最大值。從表9中也可以看出,畸形波生成、演化過程中兩次無量綱深谷最大值基本相同。與模擬實測資料得到的結(jié)果一致。對比表7和表9中的無量綱波谷值,發(fā)現(xiàn)表9中的值都大于等于表7中的值。理論上講,波浪的波高越大波谷就會相應相對較大。

4 結(jié) 論

基于完全非線性波浪數(shù)值模型,模擬了畸形波的生成、演化過程,通過對傳播波列時間和空間的演變過程的模擬及對波列中異常波浪(畸形波和深谷)特征參數(shù)分析,得出如下結(jié)論:

(1)畸形波的生成、演化過程為:從連續(xù)大波(大波波群)→深谷→畸形波→獨立大波→深谷→連續(xù)大波。

(2)在本數(shù)值計算研究范圍內(nèi),畸形波生成、演化過程中出的大波、“深谷”和畸形波均具有伴生性,畸形波生成前、后兩次形成深谷,其最大值基本相同,深谷非獨立異常波浪。

(3)畸形波生成、演化過程的空間范圍為3~7倍有效波長,時間長度為9~13倍平均周期。

[1] Osborne A R, Onorato M, Serio M. The nonlinear dynamics of rogue waves and holes in deep-water gravity wave trains [J]. Physics Letters A, 2000, 275(5-6):386-393.

[2] Stansell P. Distributions of extreme wave, crest and trough heights measured in the North Sea [J]. Ocean Engineering, 2004, 32(8~9):1015-1036.

[3] Clauss G F. Dramas of the sea: episodic waves and their impact on offshore structures [J]. Applied Ocean Research, 2002, 24(3) :147-161.

[4] Clauss G F, Hennig J, Schmittner C E. Modelling Extreme Wave Sequences for the Hydrodynamic Analysis of Ships and Offshore Structures [C]. Proceedings PRADS 2004-9th International Symposium on Practical Design of Ships and Other Floating Structures, Travemünde, Germany.

[5] Walker D A G, Taylor P H, Taylor R E. The shape of large surface waves on the open sea and the Draupner New Year wave [J].Applied Ocean Research, 2004, 26 (3~4): 73-83.

[6] Kriebel D L. Efficient simulation of extreme waves in a Random Sea [C].Abstract for Rogue waves 2000 workshops.Brest, 2000,1-2.[7] Wang Y X. Numerical wave channel with absorbing wave-maker[J].China Ocean Engineering, 1995, 9(2): 79-88.

[8] Ren B, Wang Y X. Numerical simulation of random wave slamming on structures in the splash zone [J].Ocean Engineering, 2004,31(5-6): 547-560.

Research on the generation and evolution of freak wave——time-process of the water surface elevations’ variation with positions

CUI Cheng, ZHANG Ning-chuan, PEI Yu-guo

(State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China)

Simulations of generation and evolution of freak wave were carried out based on a fully-nonlinear numerical model. The process of water surface elevation and parameters of freak wave were analyzed consequently. The results showed that freak wave was not an independent abnormal phenomenon, but was very likely accompanied by the successive large waves (huge wave sequence) and deep trough. The potential threat to vessels and maritime structures caused by the combined effect of the extreme waves and deep trough might be far beyond a transient freak, which should be paid more attention.

freak wave; generation and evolution; deep trough

P731.21

A

1001-6932(2011)04-0387-10

2010-09-16;

2011-04-08

國家自然科學基金(11602402-4022-36);博士后科學基金(20080441115)。

崔成(1984-),男,遼寧海城人,博士。主要從事隨即波浪與結(jié)構(gòu)相互作用研究。電子郵箱:chengcui1984@163.com。

主站蜘蛛池模板: 婷婷综合在线观看丁香| 婷婷丁香色| 丰满人妻一区二区三区视频| 色有码无码视频| 91久久国产综合精品女同我| 久久精品中文字幕免费| 亚洲无码91视频| a天堂视频在线| 国产在线自在拍91精品黑人| 91九色最新地址| 国产91特黄特色A级毛片| 精品国产女同疯狂摩擦2| 亚洲二三区| 九九精品在线观看| 青青青视频免费一区二区| 看国产毛片| 久久这里只有精品国产99| 国产一在线| 黄色三级毛片网站| 人妻丰满熟妇啪啪| 午夜丁香婷婷| 亚洲天堂啪啪| 亚洲国产日韩欧美在线| 欧美第二区| 亚洲综合欧美在线一区在线播放| 成人午夜网址| 欧美性精品| 国产sm重味一区二区三区| 国产高清在线精品一区二区三区| 国内精品久久人妻无码大片高| 欧美精品一区在线看| 欧美综合成人| 国产又粗又猛又爽视频| 国产成人禁片在线观看| 亚洲无码37.| 亚洲国产精品VA在线看黑人| 国产三级精品三级在线观看| 亚洲欧美精品在线| 亚洲开心婷婷中文字幕| 亚洲日韩第九十九页| 欧美另类一区| 亚洲国产亚综合在线区| 91丨九色丨首页在线播放| 日韩第九页| 1769国产精品免费视频| 在线免费不卡视频| 亚洲无码在线午夜电影| 青青久在线视频免费观看| 国产乱人伦精品一区二区| 熟女日韩精品2区| 98超碰在线观看| 久久久久久久97| www.91中文字幕| 91无码网站| 国内精品久久人妻无码大片高| 国产激情无码一区二区三区免费| 精品国产aⅴ一区二区三区| 亚洲精品无码AV电影在线播放| 国产又大又粗又猛又爽的视频| 国产18在线| 亚洲天天更新| 亚洲熟妇AV日韩熟妇在线| 亚洲日韩在线满18点击进入| 国产精品无码AV片在线观看播放| 毛片网站观看| 99久久无色码中文字幕| 天天躁夜夜躁狠狠躁图片| 国产视频大全| 国产一区二区三区在线精品专区| 国产福利微拍精品一区二区| 日本午夜网站| 国产精品xxx| 国产毛片高清一级国语| 97视频在线精品国自产拍| 国产精品99在线观看| 久久久精品国产SM调教网站| 国产午夜在线观看视频| 美女免费精品高清毛片在线视| 国产色伊人| 国内精品免费| 亚州AV秘 一区二区三区| 国产在线一区视频|