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

基于雙矩形海灣模型的泰國灣潮汐潮流研究*

2013-09-20 05:42:36方國洪崔欣梅
海洋科學進展 2013年4期
關鍵詞:區域模型

吳 頔,方國洪*,滕 飛,崔欣梅

(1.國家海洋局 第一海洋研究所,山東 青島266061;2.海洋環境科學和數值模擬國家海洋局重點實驗室,山東 青島266061)

1920年,Taylor將等深矩形海灣里的潮波,分解為入射和反射的Kelvin波以及一族Poincaré波[1](稱為Taylor問題[2]),此后眾多學者在此基礎上進行研究和推廣。Godin[3-4]將幾個“Taylor模型”組合并考慮開邊界強迫。陳宗鏞[5]、方國洪和王仁樹[6]、Rienecker和 Teubner[7]以及趙進平和侍茂崇[8]在經典 Taylor模型的基礎上引入摩擦項,考慮了能量耗散問題。1973年,Brown[9]成功地采用了Defant的配置法[10],克服了Kelvin波不完全反射時求解困難的問題,同時也使配置法成為拓展Taylor問題的有力工具,為后續的許多研究者沿用。但上面的研究只適用于半無限長海區。1991年,Fang等[11]將方國洪和王仁樹的解[6]加以拓展,考慮了灣口處的強迫邊界條件,給出適用于有限區域、包含摩擦效應的Taylor問題通解,進一步推廣了 Taylor解析模型適應的范圍。之后 Xia等[12],Carbaijal[13],Rizal[14-15],Molen等[16],Kagan等[17],Roos和Schuttelaats[18-19]均針對不同海區采用有限區域的Taylor問題解進行研究。夏綜萬[20]對Taylor問題作出系統論述。

這些研究,除Rizal[15]外,都主要關注于中高緯度海域,所得潮波系統在北半球均為逆時針旋轉的系統。位于9°N左右的泰國灣,其科氏參量遠小于M2潮波角頻率,存在一個順時針旋轉的M2潮波系統,特別值得予以理論上的解釋。我們以Fang等[11]的模型為基礎,建立了一個雙矩形的海灣模型,以研究泰國灣的潮波問題。該模型不但考慮了摩擦耗散和開邊界強迫,其雙矩形結構囊括泰國灣中的曼谷灣區域,使得模擬結果更接近真實。

1 模型的設置與求解

1.1 模型設置

新建立的雙矩形海灣模型由2個不等深矩形海灣組成(圖1),左右2個矩形區域分別稱為區域一和區域二。每個區域獨立滿足如下假設:1)摩擦力與流速的一次方成比例;2)科氏參量在海灣內為一常數;3)不計作用于海水質量上的引潮力;4)直線表示壁面,其法向速度為0;開邊界(圖1中x=0處虛線)存在外海輸入潮波強迫。

圖1 雙矩形海灣模型Fig.1 A sketch map of the two-rectangular-gulf model

假設區域一和區域二的水深均為常數,記為h1和h2;2個區域長度分別為L1和L2;寬度分別為B1和B2。區域一的范圍在x方向從0到x1(x1=L1),在y方向從0到y3(y3=B1);區域二的范圍在x方向從x1到x2(x2=L1+L2),在y方向從y1到y2(y2=y1+B2)。

1.2 控制方程和邊界條件

取潮波運動方程和連續方程:

其中分別為x,y方向潮流流速為從海平面算的潮位高度;f=2Ωsinφ為科氏參數,其中Ω(7.29×10-5s-1)為地轉角速率,φ為緯度;γ,h和g分別為摩擦系數、海區水深和重力加速度。將2個區域的潮位高度及流速分別表示為(j=1時,0≤x≤x1;j=2時,x1≤x≤x2),且滿足式式(1)~式(3)。邊界條件:

其中)為x=0處水位,即開邊界強迫。銜接條件:

式(7)保證了2個區域的潮波在連接處水位、體積輸運及能通量的連續性。

設和eiσt成比例,即存在使得:

則式(1)~式(3)可改寫:

其中,

又由式(9)~式(11)計算得出:

其中,

kj=σ/cj為Kelvin波的波數為重力長波的相速度。

相應的邊界條件式(4)~式(6)改寫:

其中,Z0和θ0為開邊界處水位的振幅和遲角。

銜接條件:

1.3 方程解法

1.3.1 方程的通解

根據式(14)和(15)的條件,可以給出滿足方程的四族Kelvin波式(20)~式(23)和四族Poincaré波式(24)~式(27):

其中,。sj,n取實部為正的值。將以上代入式(13)可以推出:

式(20)~式(27)中aj,bj,κj,n和λj,n為待定常數,需通過下列關系計算。將式(20)~式(27)代入式(16)~式(19),則壁面條件式(16)和(17)改寫:

開邊界條件式(18)改寫:

銜接條件式(19)改寫:

1.3.2 配置法

采用配置法求方程組的近似解。先將式(24)和(26)所示的Poincaré波族截斷至N1階,在x=0,0≤y≤y3和x=x1,0≤y≤y3兩條線段上均勻取2 N1+2個點。設線段x=x1,y1≤y≤y2上取了J=N2+1個點,則在線段x=x2,y1≤y≤y2上對應取N2+1個點。將這2 N1+N2+3個點坐標代入邊界條件式(28)~(32)中可得2 N1+2 N2+4個方程,不難解出aj,bj,kj,nj和λj,nj(nj=1,2,…,Nj;j=1,2)。

2 雙矩形模型在泰國灣中的應用

泰國灣介于中南半島和馬來半島之間(圖2)。地理范圍為99°10′~105°00′E,6°00′~13°30′N,灣口向東南開,西北端深入的湄南河口部分稱為曼谷灣。

由于泰國灣處于低緯度地區,受弱科氏加速度的影響,其潮波系統相對復雜,尤其是頻率較大的半日潮。以M2分潮為例,其泰國灣灣口的順指針旋轉的無潮點有悖于北半球的一般規律,一直受到相關學者的關注,而其泰國灣灣頂無潮點的存在性及位置,直至今日各學者給出的結果仍不盡相同。

在1944年Dietrich就給出基本不反映泰國灣潮波主要結構的同潮圖。后來 Defant[10]、Ye和 Robinson[21]、俞慕耕[22]給出的同潮圖,雖結果有所改善,但相互之間的差異十分顯著。此后,丁文蘭[23]、Fang[24]、Yanagi和 Toshiyuki[25]、Fang等[26]、Zu等[27]均給出了泰國灣的同潮圖。其中 M2分潮在泰國灣灣口的無潮點基本確定,但灣頂無潮點的存在及位置依舊不盡相同。上述研究方法均局限于資料分析和數值模擬,未給出理論解釋。

此外,由于曼谷灣的存在,使得單一矩形區域的Taylor模型不足以刻畫泰國灣的潮汐分布。借助新建立的雙矩形海灣模型來模擬該區域K1和M2分潮的潮波系統,并針對M2分潮的影響因素,包括曼谷灣和開邊界條件,進行了相關實驗研究。

圖2 泰國灣示意圖Fig.2 A sketch map of the Gulf of Thailand

2.1 模型設置

取L1=580 mm,B1=400 km,h1=40 m,L2=110 km,B2=100 km,y1=220 km,h2=13 m。這些數值使模型等效于泰國灣,其中區域一代表泰國灣主體,區域二代表曼古灣。并取N1=19,N2=5,在線段x=0,0≤y≤y3和線段x=x1,0≤y≤y3上均勻各取20個配置點,在線段x=x2,y1≤y≤y2,上均勻取5個點(圖3)。對 K1和 M2分潮,σ分別取7.286 7×10-5s-1、1.405 2×10-4s-1。μj則仿照Fang等[11]中所給量級取0.1。根據泰國灣及曼谷灣所處緯度,取科氏參數f1=0.22×10-4s-1;f2=0.33×10-4s-1。開邊界條件的調和常數采用Fang等[26]數值模擬結果(采用東8時區遲角)。

采用1.3.2節中方法解出相關參數出aj,bj,kj,nj和λj,nj(nj=1,2,…,Nj;j=1,2),代入式(20)~式(27)即可解出對應分潮的ζj,uj,vj。根據式(8)引入時間因子eiσt之后取實部可得:

圖3 配置點示意圖Fig.3 Setting up the Collocation points

為研究海灣模型的潮汐分布,利用式(33)導出相應分潮遲角θj及振幅Zj關系式:

2.2 潮汐分布

根據上述方法得到泰國灣-曼谷灣理論解(圖4a和圖5a);使用文獻[26]數據重繪同潮圖(圖5a和圖5b)。其中K1分潮結構較為簡單,圖4a和圖4b中潮汐分布格局一致。在泰國灣中存在一逆時針旋轉的無潮點,其位置偏向灣口,并向西岸偏移。振幅則隨著向灣頂的深入而不斷加大,最大振幅出現在曼谷灣灣頂。

圖4 K1分潮同潮圖Fig.4 Co-tidal chart of tidal constituent K1

圖5 M2分潮同潮圖Fig.5 Co-tidal chart of tidal constituent M2

M2分潮頻率大,波長僅為K1分潮的一半,故潮波結構比K1分潮復雜的多。圖5a中顯示,M2分潮振幅在泰國灣(除曼谷灣部分)較小,最大振幅不超過0.3 m,在曼谷灣中向灣頂方向不斷增大;遲角分布則較復雜,出現2個無潮點。其中一個無潮點位于泰國灣灣頂、曼谷灣口外,呈逆時針旋轉。另一個無潮點則靠近泰國灣灣口,偏向東岸,呈順時針旋轉,該無潮點不遵循北半球潮波逆時針旋轉的一般規律。對比圖5a和圖5b,M2分潮振幅分布格局一致。

為分析K1和M2分潮的構成,分別給出這2個分潮在區域一中的Kelvin波和Poincaré波分布(圖6和圖7)。由圖6中可以看出K1分潮潮波由于其頻率小、波長長,雖然處于較低緯度的泰國灣區域,其組成仍然以Kelvin波為主,并表現出明顯的Kelvin波性質。Kelvin波受北半球向右的科氏加速度影響,入射波東岸強西岸弱,反射波東岸弱西岸弱。合成結果構成了一個逆時針旋轉的潮波系統。由于模型考慮了摩擦,潮波在傳播過程中產生能量耗散,使得反射波的能量小于入射波,總的疊加效果東岸強西岸弱,因此無潮點偏向西岸。

圖6 K1分潮在區域一中的潮波分解Fig.6 Decomposition of K1 in area I

圖7 M2分潮在區域一中的潮波分解Fig.7 Decomposition of M2 in area I

由M2分潮Kelvin波和Poincaré波的組成(圖7)可見,M2分潮由于頻率高,在低緯度區慣性頻率f比潮波頻率σ小很多,故在泰國灣區域Kelvin波較弱,且由于能量耗散,入射波遠遠強于反射波,沒有疊加形成無潮點。而Poincaré波在泰國灣灣口較強,且存在一個逆時針旋轉的無潮點。Kelvin波和Poincaré波疊加后在灣頂形成一逆時針旋轉的無潮點,灣口形成一順時針旋轉的無潮點。

考慮到Poincaré波主要是受邊界條件和銜接條件控制,我們認為除低科氏加速度及地形因素影響外,曼谷灣是影響泰國灣灣頂形成逆時針旋轉的無潮點的重要因素;開邊界條件是泰國灣灣口形成順時針旋轉的無潮點的首要因素。

2.3 潮流分布

為了研究泰國灣的潮流分布情況,根據2.1節中結果分別計算了K1和M2分潮的潮流橢圓要素(半長軸、半短軸、長軸傾角及遲角),并給出相應分潮的潮流橢圓(圖8和圖9)。

圖8 K1分潮的潮流橢圓圖(陰影表示順時針旋轉)Fig.8 K1 current ellipses(shaded part indicates clockwise rotation)

圖9 M2分潮的潮流橢圓圖(陰影表示順時針旋轉)Fig.9 M2 current ellipses(shaded part indicates clockwise rotation)

K1分潮潮流在區域一中主要呈逆時針旋轉,在區域二中主要呈順時針旋轉,在2個區域內均以往復流為主。泰國灣灣口東側和曼谷灣灣口東側出現較大流速。在泰國灣腹部,長軸平行于x軸,在灣頂則逐漸垂直于x軸,在與Taylor[1]給出的等深矩形海灣十分相似。

M2分潮潮流在泰國灣和曼谷灣的灣口附近都呈順時針旋轉,在2個海灣的灣頂附近則都呈現逆時針旋轉,與Fang等的理論分析[11]一致。

2.4 潮能通量

潮能通量強度又叫做能通量密度,表示單位時間內通過單位寬度的垂直斷面的潮波能量,其計算公式:

式中(Φx,Φy)為能通量密度在(x,y)方向的分量;T 為潮波周期;ρ為海水密度,取1 025 kg/m3;Z 和θ分別為潮汐的振幅和遲角;(U,V)為(x,y)方向潮流分量的振幅,(ξ,η)為相應的遲角。

圖10 K1分潮能通量密度分布圖Fig.10 The tidal energy flux density vectors of constituent K1

圖11 M2分潮能通量密度分布圖Fig.11 The tidal energy flux density vectors of constituent M2

根據式(37)及前文結果,計算了K1和M2分潮的潮能通量密度,并繪制潮能通量密度分布圖以刻畫K1和M2分潮能通量流向(圖10和圖11)。圖10中K1分潮能量從東岸輸入,并沿東岸向灣頂傳輸,自東向西掃過灣頂,并沿西岸輸出,在傳輸過程中,能量逐漸減弱。這體現出Kelvin波在北半球傳播的一般規律。圖11中M2分潮能量較小,正如2.2節所述,受開邊界及弱科氏加速度影響,其能量從東岸輸入,但直接流向西岸,且沿西岸傳輸,緩慢右轉,轉向東岸后,沿東岸傳輸直至灣頂,后沿灣頂邊界自東向西傳輸,在曼谷灣處分為2支,一部分流入曼谷灣,另一部分繼續沿灣頂傳輸。

3 敏感性實驗

曼谷灣對M2分潮在泰國灣灣頂形成逆時針旋轉的無潮點起重要作用;開邊界條件是M2分潮在泰國灣灣口形成順時針旋轉的無潮點的必要條件。為了驗證上述觀點,我們設計了如下2組實驗:1)曼谷灣實驗:刪除代表曼谷灣的區域二,其他條件同2.1節;改變區域二位置,其他條件同2.1節。2)開邊界實驗:開邊界條件采用等振幅邊界,其他條件同2.1節;開邊界條件采用等遲角邊界,其他條件同2.1節。

3.1 曼谷灣存在和位置的敏感性實驗

針對M2分潮,改變曼谷灣的存在及位置,控制其他條件與2.1節中一致,以探討曼谷灣對泰國灣M2分潮的影響。

方案1 只保留區域一,即建立一個單一矩形海灣模型(圖12),其他條件保持與2.1節相同,以此排除曼谷灣的影響。實際計算中在模型兩端均勻各取20個配置點(圖12a),沿用Fang等的解法[11],計算該矩形海灣模型的潮汐分布,得到M2分潮的同潮圖如圖12b所示。

圖12 泰國灣的矩形海灣模型Fig.12 A single-rectangular-gulf model for the Gulf of Thailand

對比圖12b與圖5a結果可見,M2分潮在泰國灣仍然存在一順時針旋轉的無潮點,但其位置向灣內深入,而灣頂逆時針旋轉的無潮點完全消失。

方案2 令y1=160 km;其他條件同2.1節。

方案3 令y1=100 km;其他條件同2.1節。

方案2和3的結果分別示于圖13。對比圖13與圖5a結果可知,隨著曼谷灣向東平移,M2分潮在泰國灣灣口順時針旋轉的無潮點位置變化不大,灣頂的無潮點向灣頂方向發生退化,但可以看出這個無潮點隨曼谷灣位置東移向東北移動。

綜合上述3個方案,可以看出曼谷灣是影響M2分潮在泰國灣灣頂形成逆時針旋轉的無潮點的重要因素,其存在促使M2分潮在泰國灣灣頂形成無潮點,且曼谷灣的位置決定了該無潮點的位置。

圖13 曼谷灣實驗的M2分潮同潮圖Fig.13 The co-tidal chart of M 2 resulted from the sensitive experiments in the Gulf of Bangkok

3.2 開邊界條件的敏感性實驗

著眼于研究泰國灣M2分潮對開邊界的響應,具體方案:

方案1 開邊界振幅均取0.2 m,其他條件同2.1節;

方案2 開邊界遲角均取0°,其他條件同2.1節。

依舊采用2.1節中的方法,分別計算方案1、方案2模型的潮汐分布,并給出相應同潮圖(圖14)。

圖14 開邊界實驗的M2分潮同潮圖Fig.14 The co-tidal chart of M2 resulted from the sensitive experiments under the open boundary condition

與圖5a結果比較可見,當開邊界振幅為常數時,其振幅結構變化不大,灣口的無潮點退化至東岸,但依然保持順時針旋轉的結構。當開邊界遲角給定均勻值時,灣口順時針的無潮點完全消失,潮波進入海灣后差不多以前進波的形式向灣頂推進,并在曼谷灣灣口外形成一個逆時針旋轉的無潮點。對比2組實驗結果,不難看出,改變開邊界條件不影響灣頂形成逆時針的無潮點;而M2分潮在泰國灣灣口順時針旋轉的無潮點,主要依賴于開邊界遲角自東向西增加的梯度條件,即潮波傳入泰國灣的方向;而灣口的振幅變化對該無潮點形成僅起到加強作用。遲角自東向西的增加則是南海M2分潮潮波自東北向西南傳播的結果。

5 結 語

本文利用雙矩形海灣模型研究了泰國灣的潮波問題。研究表明:

1)K1分潮潮波在泰國灣中主要表現Kelvin波的性質;

2)M2分潮在泰國灣中Kelvin波的性質較弱,灣口和灣頂主要受Poincaré波控制;

3)曼谷灣的存在是M2分潮在泰國灣灣頂形成逆時針旋轉無潮點的主要原因,其位置決定了該無潮點的位置;

4)開邊界遲角自東向西增加的梯度條件是M2分潮在泰國灣灣口順時針無潮點形成的必要條件。

[1] TAYLOR G I.Tidal oscillations in gulfs and rectangular basins[J].Proc.Lond.Math.Soc.,1922,s2-20(1):148-181.

[2] HENDERSHOTT M C,SPERANZA A.Co-oscillating tides in long,narrow bays:the Taylor problem revisited[J].Deep Sea Research and Oceanographic Abstracts,1971,18(10):959-980.

[3] GODIN G.Some remarks on the tidal motion in a narrow rectangular sea of constant depth[J].Deep Sea Research and Oceanographic Abstracts,1965,12(4):461-468.

[4] GODIN G.The M2tide in the Labrador Sea,Davis Strait and Baffin Bay[J].Deep Sea Research and Oceanographic Abstracts,1965,12(4):469-477.

[5] CHEN Z Y.A tidal mode of rectangular shallow basins[J].Oceanologia et Limnologia Sinica,1965,7(2):85-93.陳宗鏞.長方形淺水海灣的一種潮波模式[J].海洋與湖沼,1965,7(2):85-93.

[6] FANG G H,WANG R S.Tides and tidal streams in gulfs[J].Oceanologia et Limnologia Sinica,1966,8(1):60-77.方國洪,王仁樹.海灣的潮汐與潮流[J].海洋與湖沼,1966,8(1):60-77.

[7] RIENECKER M M,TEUBNER M D.A note on frictional effects in Taylor's problem[J].Journal of Marine Research,1980,38(2):183-191.

[8] ZHAO J P,SHI M C.The tidal wave reflection problem and the influence of friction in Semi-enclosed rectangular basin[J].Acta Oceanologica Sinica,1988,10(3):259-269.趙進平,侍茂崇.半封閉矩形海灣中潮波反射問題及摩擦的影響[J].海洋學報,1988,10(3):259-269.

[9] BROWN P J.Kelvin wave reflection in a semi-infinite canal[J].Journal of Marine Research,1973,31(1):1-10.

[10] DEFANT A.Physical Oceanography[M].New York:Pergamon Press,1961:590.

[11] FANG Z,YE A L,FANG G H.Solutions of tidal motions in a semi-closed rectangular gulf with open boundary condition specified[M]∥Tidal Hydrodynamics,New York:John Wiley &Sons Inc.,1991:153-168.

[12] XIA Z W,CARBAJAL N,SUNDERMANN J.Tidal current amphidromic system in semi-enclosed basins[J].Continental Shelf Research,1995,15(2-3):219-240.

[13] CARBAJAL N.Two applications of Taylor's problem solution for finite rectangular semi-enclosed basins[J].Continental Shelf Research,1997,17(7):803-817.

[14] RIZAL S.The role of non-linear terms in the shallow water equation with the application in three-dimensional tidal model of the Malacca Strait and Taylor's Problem in low geographical latitude[J].Continental Shelf Research,2000,20(15):1965-1991.

[15] RIZAL S.Taylor's problem-influences on the spatial distribution of real and virtual amphidromes[J].Continental Shelf Research,2002,22(15):2147-2158.

[16] MOLEN J,GERRITS J,SWART H E.Modelling the morphodynamics of a tidal shelf sea[J].Continental Shelf Research,2004,24(4-5):483-507.

[17] KAGAN B A,ROMANENKOV D A.Effect of hydrodynamic properties of the sea bottom on the tidal dynamics in a rectangular basin[J].Atmospheric and Oceanic Physics,2006,42(6):777-784.

[18] ROOS P C,SCHUTTELAATS H M.Horizontally viscous effects in a tidal basin:extending Taylor's problem[J].Journal of Fluid Mechanics,2009,640:421-439.

[19] ROOS P C,SCHUTTELAATS H M.Influence of topography on tide propagation and amplification in semi-enclosed basins[J].Ocean Dynamics,2011,61(1):21-38.

[20] XIA Z W.Taylor problem of tidal wave[M].Beijing:Ocean Press,2011.夏綜萬.海洋潮波的Taylor問題[M].北京:海洋出版社,2011.

[21] YE A L,ROBINSON I S.Tidal dynamics in the South China Sea[J].Geophysical Journal of the Royal Astronomical Society,1983,72(3):691-707.

[22] YU M G.Preliminary discussion on tidal characteristic of South China Sea[J].Acta Oceanologica Sinica,1984,6(3):293-300.俞慕耕.南海潮汐特征的初步探討[J].海洋學報,1984,6(3):293-300.

[23] DING W L.Distribution of tides and tidal currents in the South China Sea[J].Oceanologia et Limnologia Sinica,1986,17(6):468-480.丁文蘭.南海潮汐和潮流的分布特征[J].海洋與湖沼,1986,17(6):468-480.

[24] FANG G H.Tide and tidal current charts for the marginal seas adjacent to China[J].Chinese Journal of Oceanology and Limnology,1986,4(1):1-16.

[25] YANAGI T,TOSHIYUKI T.Clockwise phase propagation of semi-diurnal tides in the Gulf of Thailand[J].Journal of Oceanography,1998,54(2):143-150.

[26] FANG G H,KWOK Y K,YU K J,et al.Numerical simulation of principal tidal constituents in the South China Sea,Gulf of Tonkin and Gulf of Thailand[J].Continental Shelf Research,1999,19(7):845-869.

[27] ZU T T,GAN J P,EROFEEVA S Y.Numerical study of the tide and tidal dynamics in the South China Sea[J].Deep Sea Research Part I:Oceanographic Research Papers,2008,55(2):137-154.

猜你喜歡
區域模型
一半模型
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
分割區域
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
關于四色猜想
分區域
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
主站蜘蛛池模板: 99久久精品免费看国产免费软件| 亚洲成av人无码综合在线观看| 中文字幕自拍偷拍| 亚洲伊人天堂| 欧美福利在线观看| 欧美一区二区三区不卡免费| 免费国产在线精品一区 | 中文字幕在线不卡视频| 国产网站免费观看| 久久精品国产精品一区二区| 欧美精品伊人久久| 国产va在线| 欧美精品v日韩精品v国产精品| 亚洲天堂自拍| 就去吻亚洲精品国产欧美| 91精品国产无线乱码在线| 亚洲swag精品自拍一区| 欧美亚洲日韩中文| 成人永久免费A∨一级在线播放| 狠狠色狠狠综合久久| 免费在线a视频| 成人国产三级在线播放| 在线观看国产网址你懂的| 国产精品亚洲一区二区在线观看| 97久久免费视频| 国产成人无码播放| 超薄丝袜足j国产在线视频| 色网在线视频| 亚洲第一区欧美国产综合| 欧美一区日韩一区中文字幕页| 婷婷中文在线| 国产精品漂亮美女在线观看| 欧美日韩中文字幕在线| 永久免费精品视频| 精品国产毛片| 国产又大又粗又猛又爽的视频| 日韩在线永久免费播放| 亚洲日韩精品欧美中文字幕 | 好紧太爽了视频免费无码| 亚洲乱伦视频| 精品无码视频在线观看| 久久人搡人人玩人妻精品一| 国产精品v欧美| 亚洲欧洲综合| 成人毛片免费在线观看| 这里只有精品国产| 国产日韩久久久久无码精品| 亚洲美女操| 免费xxxxx在线观看网站| 精品久久久久久久久久久| 草逼视频国产| 亚洲第一黄色网| 国产欧美在线视频免费| 欧美一级在线看| 国产精品久久久久久久久| 国产免费黄| 国产一级做美女做受视频| 正在播放久久| 99九九成人免费视频精品| 精品久久久久成人码免费动漫| 亚洲国产成人久久精品软件| 亚洲av成人无码网站在线观看| 亚洲制服丝袜第一页| 亚洲欧美成人网| 亚洲欧美激情小说另类| 91福利国产成人精品导航| 黄色在线不卡| 亚洲国产综合精品一区| 99re这里只有国产中文精品国产精品| 91色在线视频| 久久永久免费人妻精品| 日本精品一在线观看视频| 欧美精品啪啪| 欧美日韩91| 日韩欧美网址| 国产精鲁鲁网在线视频| 毛片免费观看视频| 高潮爽到爆的喷水女主播视频| 亚洲欧美另类中文字幕| 亚洲综合天堂网| 欧美成人二区| 思思99思思久久最新精品|