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

近岸浪-風暴潮耦合模型在天津沿海的應用

2017-01-13 05:01:05李雪王智峰武雙全董勝賈婧張曉爽吳昊
海洋通報 2016年6期
關鍵詞:模型

李雪,王智峰,武雙全,董勝,賈婧,張曉爽,吳昊

(1.國家海洋局數字海洋科學技術重點實驗室,天津 300171;2.中國海洋大學 工程學院,山東 青島 266100;3.國家海洋信息中心,天津 300171;4.國家海洋局海洋環境信息保障技術重點實驗室,天津 300171)

近岸浪-風暴潮耦合模型在天津沿海的應用

李雪1,2,王智峰2,武雙全3,董勝2,賈婧2,張曉爽4,吳昊2

(1.國家海洋局數字海洋科學技術重點實驗室,天津 300171;2.中國海洋大學 工程學院,山東 青島 266100;3.國家海洋信息中心,天津 300171;4.國家海洋局海洋環境信息保障技術重點實驗室,天津 300171)

為了精細化描述天津沿海臺風天氣下近岸浪和風暴潮特征,基于非結構三角網格,建立近岸浪與風暴潮的耦合模型,其中臺風風場采用藤田臺風模型,近岸浪采用SWAN波浪模型,風暴潮采用ADCIRC模型。通過對幾次典型臺風暴潮數值模擬的驗證,耦合模型對風速、有效波高和增水的計算結果與實際觀測資料符合性均較好,能夠很好地反映臺風過程中天津沿海近岸浪和風暴潮特征,可以為天津的防災減災工作提供科學依據。

臺風浪;風暴潮;SWAN;ADCIRC

隨著我國海洋經濟及沿海地區人口的快速發展,建立和完善沿海地區的海洋數值預報系統顯得尤為重要,對海洋防災減災、保護國家和人民的生命財產安全具有重要的理論指導和現實意義。天津地處中國北方黃金海岸的中部,是山東半島與遼東半島的交匯點,華北、西北廣大地區的重要出??冢c世界上180多個國家、地區的400多個港口保持著貿易往來,發展潛力巨大。然而,天津地勢低洼,臺風暴潮災害頻發,是沿海遭受災害嚴重的地區之一,其造成的經濟損失從20世紀90年代后已有明顯增加的趨勢(葉鳳娟等,2012)。因此,在天津海域建立完整的海洋數值預報系統十分必要,也是天津防災減災的重要研究課題(傅賜福等,2013)。

臺風是熱帶風暴潮的主要驅動力,臺風的路徑、強弱、影響半徑等直接關系到風暴潮災害的強度。臺風產生的大風,往往在海面上掀起巨浪,進一步升高水位。如果臺風登陸時刻,風暴潮與天文大潮高潮位相遇,則危害性更大(馮士筰,1982)。因此,波浪和風暴潮通過不同的機制相互影響和制約,只有全面考慮這兩者之間的相互作用,才能準確地模擬出真實的水動力和波浪條件,為其他研究提供基礎。

在實際海岸環境條件下,風暴潮與近岸波浪演化之間復雜的相互作用,主要表現在兩個方面:一方面是波浪傳播對水動力產生影響,如波浪傳播變形所引起的增減水和波生近岸流 (Xie et al,1999,2001;鄭立松等,2010);另一方面則是水動力變化對波浪產生影響,如水位和流速、流向的變化引起波浪變形(白志剛等,2012)。這種相互作用在強風過程中更為明顯。

關于風暴潮對波浪影響的研究,主要通過采用考慮流速和潮位變化的波浪傳播變形方程來體現,但以往大部分的研究采用的是恒定流場或波浪場,沒有考慮實際海洋環境中波浪、水動力均隨時間和空間變化而產生的耦合作用(夏波等,2006,2012),或是僅針對某一臺風發生過程對天津沿海的影響進行了研究,并未形成能夠綜合體現近岸浪和風暴潮對天津沿海水動力環境影響的模擬系統。

本文將采用第三代波浪模擬的SWAN模型和近岸風暴潮ADCIRC模型,建立精細化的近岸浪-風暴潮耦合模型,并利用該系統對天津沿海的臺風浪和風暴潮進行模擬計算,從而為天津的防災減災工作提供科學依據。

1 近岸浪-風暴潮耦合模型的建立

1.1藤田臺風模型

風場作為引發臺風、風暴潮的主要因子,對臺風暴潮的預報和研究起到至關重要的作用。風場決定波浪場,表面風場對近岸區流場的影響更大,通常近岸區的風增水和風生流現象尤為明顯。迄今,有關臺風暴潮的研究中通常采用各種經驗的或半經驗半理論的臺風模型,其方法主要有三大類:第一類是圓對稱風場模型(高橋,1939;藤田,1952;Jelesnianski,1965;Holland,1980);第二類為改進的非對稱風場模型,如非對稱結構風場模型(章家彬等,1986),橢圓形對稱的風壓模型(陳孔沫,1994),特征等壓線風場模型(朱首賢等,2003);第三類為氣象數值預報模型,如MM5,WRF。

本文采用便于描述計算的藤田模型進行風場計算。風場由氣壓梯度和臺風移動導致的風速合成而成,其中梯度風公式為:

式中,f為科氏力系數,ρα為空氣密度,r是離開臺風中心的距離,P是離開臺風中心距離為r處的氣壓;而臺風移動引起的環境風場由Ueno公式給出:

式中,Vtx和Vty表示臺風移動速度的x和y分量,R為臺風最大風速半徑,因此合成風場為:

其中,Wx和Wy表示合成風速的x和y分量,C1和C2為經驗系數,θ為考慮大氣邊界層影響之后梯度風的偏角(稱為入流角),P0是臺風中心氣壓,P∞是距離臺風中心無限遠處的氣壓。

1.2 SWAN波浪模型

波浪模擬采用第三代波浪模型的SWAN(Simulating Waves Nearshore)模型,該模型通過采用波作用量平衡方程,在源項中計入能量輸入和損耗項(底摩擦、破碎、白浪)、波與波之間非線性相互作用等(Booij et al,1999),從而比較全面合理地描述風浪生成及其在近岸區的演化過程。

在直角坐標系下,波作用量平衡方程(Ris et al,1999)可表示為:

式中:N=E(σ,θ,x,y,t)/σ為波作用密度譜;x,y為地理坐標,t為時間,θ為傳播方向,σ為相對頻率,cxy、cσ、cθ分別為波能量在地理空間和譜空間上的傳播速度。方程左側第一項為N隨時間的變化率;第二項為N在地理空間方向上的傳輸,第三項是由于流場和水深所引起的N在波浪相對頻率σ空間的變化;第四項為N以傳播速度Cθ在譜分布空間θ(譜分量中垂直于波峰線的方向)上的傳播;方程右端項的S為以譜密度表示的源匯項,具體描述為:

式中,分項Sin為風能輸入,Snl3、Snl4分別為三波和四波相互作用的非線性波能傳播,Sds,w為白帽耗散,Sds,b為底摩擦耗散,Sds,br為水深變淺引起的波能破碎耗散。

由于近岸區的波浪場對近岸區流場的改變主要體現于波浪淺水變形產生的波浪輻射應力,SWAN模型中波浪輻射應力張量的計算公式為:

式中,ρ0為水密度;g為重力加速度;n為群速度與相速度之比。

1.3 ADCIRC風暴潮模型

本文采用的風暴潮模型 ADCIRC(Advanced Circulation Model),是基于有限元方法的跨尺度(大洋、近岸、河口)水動力計算模型(Luettich et al,1990)。ADCIRC模型采用非結構網格,可以使模型在水深變化劇烈、岸線復雜的地方有較高的分辨率,而在地形變化緩慢的地方分辨率相對低一些,這樣既可以滿足計算的要求,又可以節省計算時間。同時該模型在空間上采用有限元方法進行求解,時間上采用有限差分方法,計算效率高且穩定。

在直角坐標系下,模型的控制方程為:

式中,ζ為從平均海平面算起的水位高度,H=ζ+h為總水深;U和V分別表示坐標軸x和y方向的垂向平均流速,f為科氏力系數,Ps為表面大氣壓力,η表示牛頓潮勢作用,α為有效彈性系數;τs,wind和τs,wave表示表面風應力項和波浪輻射應力項,τb、D、M分別代表底部切應力項、擴散項以及側向應力項。

為了避免或減小Galerkin有限元離散所帶來的振蕩、不守恒等數值問題,ADCIRC模型采用對短波具有阻尼作用的通用波動連續性方程GWCE(Generalized Wave Continuity Equation)來代替原有的連續性方程(Blain et al,2004),其優越性體現于在沒有對流加速度的情況下,對長波的計算更加精確;且方程對水位和速度的求解是自然解耦的,即先求得水位,后解得流速。

其中

式中,Qx=UH,Qy=VH表示流量分量,τ0為優化相位傳播特性的數學參數。臺風施加于海面的切應力用以下公式計算:

式中,拖曳系數Cd=10-3×(0.75+0.067),如果Cd>0.003,則取Cd=0.003。

波浪輻射應力τs,wave由公式(19)和公式(20)計算得到。通過對SWAN模型計算所得輻射應力張量的計算,τs,wave將波浪淺水變形代入風暴潮計算模型中,以體現波浪對風暴潮的影響(Hedges et al,2004;Liu et al,2008,2009;Dietrich et al,2012)。同時,近岸區的水深和流速隨水位的漲落而不斷變化,ADCIRC計算所得各節點水深和流速的時間序列作為SWAN模型的輸入條件,可模擬出流場對波浪場的影響,從而實現天津近岸波浪-風暴潮的耦合計算(Sebastian et al,2014;Li et al,2014)。

2 耦合模型的驗證

2.1耦合模型參數設置

為了更好地模擬強風過程中風浪及流場的變化,減小邊界條件誤差對計算結果的影響,計算區域應該選擇得足夠大。但是,波浪與流場之間的非線性相互作用在近岸地區最為明顯,要較為準確地描述這種相互作用,需要有足夠精確的地形測量數據和采取足夠小的空間和時間步長。因此,為了快速、準確地模擬風浪之間的相互作用,本文使用由NOAA(National Oceanic and Atmospheric Administration)提供的全球陸地海洋GEODAS高程數據ETOPO1和渤海近海海圖提取的水深數據拼接而成的水深數據,水深分布見圖1;計算區域在包括本文所研究區域的基礎上,擴展至整個渤海和部分黃海海域;模型采用三角網格,并在天津海域進行加密處理,如圖2所示;計算區域網格數為133 136個,節點數為68 094個。

近岸波浪模式的計算時間步長為20 min;方向步長為10°;頻率計算范圍為0.05~1.0 Hz;底摩阻公式采用認為底摩擦系數是底床粗糙度與實際波浪條件函數的Madsen公式進行計算,粗糙度取為0.05;白帽采用Komen公式;其他參數均取默認值進行計算。

風暴潮模型的邊界條件設為初始條件為t=0時,ζ=U=V=0;海岸邊界條件為邊界的法向速度為0;潮汐開邊界由M2,S2,K1,O1共4個分潮驅動計算。為滿足CFL條件為要求,時間步長取為30 s;底摩擦系數采用適用于大洋與近海的線性與二次律混合形式;海面風應力與風速呈二次平方律關系,風拖曳系數采用Garratt公式(Luettich et al,2000),模型通過設定最小水深判斷網格的干濕狀態。

圖1 計算區域水深分布圖

圖2 計算區域局部加密網格圖

2.2典型臺風過程

將上述過程建立的近岸浪-風暴潮耦合模型應用于天津海域,為了驗證耦合系統的可靠性,現針對4次典型臺風過程進行分析計算。這4次臺風過程分別為:1972年“7203”號臺風Rita、1985年“8509”號臺風Mimie、2005年“0509”號臺風Matsa和2012年“1201”號臺風Damrey,其對應的臺風發生過程時間及塘沽站實測信息如表1所示,臺風移動路徑見圖3。

為檢驗計算所用風場的精度,將由藤田臺風模型計算所得塘沽站的計算風速、風向與實測風速、風向對比如下。從圖中可以看出,風場計算模式對各次臺風暴潮過程中最大風速的計算是較為準確的,對臺風風場過程的描述也是基本成功的,可以作為耦合模型的輸入條件。

表1 臺風發生過程及塘沽站實測值

圖3 臺風移動路徑圖

2.3近岸浪計算結果驗證

圖5分別給出了4次臺風過程中塘沽站的波高模擬值與實際觀測值隨時間變化的對比圖,其中紅色實心點表示實測值,藍色虛線表示模型不進行耦合的計算值,黑色實線表示耦合模型計算值。

圖4 (a) 7203號臺風塘沽站風速、風向對比圖

圖4 (b) 8509號臺風塘沽站風速、風向對比圖

圖4 (c) 0509號臺風塘沽站風速、風向對比圖

圖4 (d) 1210號臺風塘沽站風速對比圖

結合臺風路徑圖和風速圖可以看出,7203號臺風于7月27日左右登陸渤海灣,SWAN模型很好的模擬出了27日的波高峰值,峰值誤差僅為0.6cm,但在峰值過后模擬值下落得比觀測值快,其原因可能在于對應時間的風速模擬值衰減得較快。8509號臺風雖然未直接襲擊天津海域,但其觀測波高在19日高達近2 m,模型模擬出的波高峰值與觀測值相差1.5 cm,相對誤差為0.72%,峰值模擬十分準確。1210號臺風于2012年8月2日登陸江蘇與山東省,而后于8月3日12時轉入天津海域,其威力原本在陸地登陸期間有所減弱,但塘沽站觀測得最大波高仍高達3.17 m,模擬得波高峰值為3.41 m,相對誤差為7.27%,波高時變過程的相關系數為0.95,模擬結果良好。而0509號臺風無塘沽站的波高觀測數據,故本文采用有實測數據的山東石島站進行替代驗證,由圖可以看出耦合模型對此次臺風過程的波高變化趨勢模擬較為準確。

圖5 (a) 7203號臺風塘沽站波高模擬值與觀測值對比圖

圖5 (b) 8509號臺風塘沽站波高模擬值與觀測值對比圖

圖5 (c) 1210號臺風塘沽站波高模擬值與觀測值對比圖

圖6 0509號臺風石島站波高模擬值與觀測值對比圖

由于塘沽測站只進行波高的測量,而無波向的實測數據,故無法進行波向的結果對比。由以上4次臺風過程波高觀測值與計算值的對比分析可以看出,本耦合模型的風浪模塊能夠有效模擬出天津海域近岸風浪的峰值及時變過程,可以作為后續防災減災工作的依據。

2.4風暴潮增水計算結果驗證

為了驗證耦合模型對風暴潮增水的模擬效果,圖7(a)、(b)、(c)、(d)分別給出了4次臺風過程模擬的增減水與實際觀測值時間序列的對比圖,其中紅色實心點表示實測值,藍色虛線為風暴潮模型不進行耦合的計算值,黑色實線表示耦合模型計算結果。此外,由于渤海是我國最北部海域,除在夏季遭受臺風暴潮影響外,冬半年的氣象活動也十分劇烈,溫帶氣旋、寒潮等天氣過程所伴隨的大風也是天津海域風暴潮的誘因之一,其影響也需受到重視。故在此對2009年4月的溫帶氣旋的增減水過程進行了模擬,模型所用驅動條件為歐洲中期天氣預報中心ECMWF(European Centre for Medium-Range Weather Forecasts)的再分析數據。圖8為模擬增減水過程與實測值的對比。

圖7、圖8的對比結果表明,耦合模型對增減水過程模擬結果良好,5次風暴潮過程的峰值相對誤差均不超過15%,相關系數均在0.85以上。

結合波高與水位的對比圖可以看出,在波高出現峰值的時刻,水位也相應得以升高。且8509號臺風并未直面登陸天津海域,其增水卻高達1.5 m,有風暴潮與天文大潮高潮位相遇的可能性,這也進一步驗證近岸水動力與波浪之間復雜的相互作用與影響。而從波高及增水的未耦合值與耦合值的對比可以看出,耦合模型在一定程度上模擬出了這種相互作用,使得模型的計算值更為接近實測值。故耦合模型的波浪模塊對風暴潮波高的模擬,以及風暴潮模塊對風暴潮增水的模擬更為準確,本文建立的近岸浪-風暴潮耦合模型是科學且可靠的。

圖7 (a) 7203號臺風塘沽站增水模擬值與觀測值對比圖

圖7 (b) 8509號臺風塘沽站增水模擬值與觀測值對比圖

圖7 (c) 0509號臺風塘沽站增水模擬值與觀測值對比圖

圖7 (d) 1210號臺風塘沽站增水模擬值與觀測值對比圖

圖8 2009年4月溫帶氣旋風暴潮塘沽站增水模擬值與觀測值對比圖

3 結論

考慮到風暴潮與近岸浪之間的相互作用機制,本文構建了近岸浪-風暴潮耦合模型,并對發生在天津海域的3次典型臺風過程進行了模擬計算。其中臺風風場采用藤田臺風模型,近岸浪模型采用SWAN模型,風暴潮模型采用ADCIRC模型。通過近岸浪和風暴潮的驗證,可以得到以下結論:

(1)耦合模型能夠較好地模擬風浪的成長演變過程以及風暴潮對風浪傳播的影響,采用非結構化的三角形網格,能有效模擬出近岸風浪的變化情況;

(2)耦合模型能夠較好的模擬風暴潮過程中的實際增減水過程以及風浪對水位的影響;

(3)本文構建的近岸浪-風暴潮耦合模型能夠為臺風暴潮災害及其預警報提供技術支撐;為政府和國家管理部門制定沿海地區發展規劃提供有力的支撐,為制定防災減災應急預案、補償、救濟措施和政策提供參考依據。

Blain C A,Massey T,2004.Application of a coupled discontinuouscontinuous Galerkin finite element shallow water model to coastal ocean dynamics.Ocean Modelling,10(3-4):283-315.

Booij N,Ris R C,Holthuijsen L H,1999.A third-generation wave model for coastal regions:Part I.Model description and validation.Journal of Geophysical Research,104(C4):7649-7666.

Dietrich J C,Tanaka S,Westerink J J,et al,2012.Performance of the unstructured-mesh SWAN+ADCIRC model in computing hurricane waves and surge.Journal of Scientific Computing,52(2):468-497.

Hedges T S,Mase H,2004.Modified Hunt's equation incorporating wave setup.Waterway,Port,Coastal and Ocean Engineering,130(3): 109-113.

Li X B,Sun X Y,Liu Y,et al,2014.Statistical and causes analysis of storm surges along Tianjin coast during the past 20 years.Marine Science Bulletin,16(1):15-24.

Liu B,Liu H,Guan C,2009.Numerical simulation of tropical cyclone intensity using an air-sea-wave coupled prediction system.Adv. Geosci,18:19-43.

Liu H,Peng M,2008.The effect of wave-current interactions on the storm surge and inundation in Charleston Harbor Hurricane Hugo 1989. Ocean Modelling,20(3):252-269.

Luettich R A,Westerink J J,2000.ADCIRC user manual:A (Parallel) Advanced Circulation Model for Oceanic,Coastal and Estuarine Waters.(http://www.unc.edu/ims/adcirc/documentv47/ADCIRC_title_ page.html).

Luettich R A,Westerink J J,Scheffner N W,1990.ADCIRC:An advanced three-dimensional circulation model for shelves,coasts, and estuaries,report I:Theory and methodology of ADCIRC-2DDI and ADCIRC-3DL.Washington D.C.:U.S.Army Corps of Engineers,Waterways Experiment Stations.

Ris R C,Booij N,Holthuijsen L H,1999.A third-generation wave model for coastal regions:Part II.Verification.Journal of Geophysical Research,104(C4):7649-7666.

Sebastian A,Proft J,Dietrich J C,et al,2014.Characterizing hurricane storm surge behavior in Galveston Bay using the SWAN ADCIRC model.Coastal Engineering,88:171-181.

The SWAN team,2011.SWAN user manual(version 40.85).Delft: Delft University of Technology.

Xie L,Pietrafesa L J,Wu K J,1999.A numerical study of wave-current interactions in the coastal zone. Proceedings,The 1999 Oceanographic Society Annual Meeting.Reno,Nevada.

Xie L,Wu K J,Pietrafesa L,Zhang C,2001.A numerical study of wave-current interaction through surface and bottom stresses:Winddriven circulation in the South Atlantic Bight under uniform winds. Journal of Geophysical Research,106(16):841-855.

白志剛,裴麗,陽磊,2012.潮間帶風電場水動力數值模擬.港工技術,49(3):1-4.

陳超輝,譚言科,李霖,等,2011.藤田公式改進及其在臺風預報中的應用.氣象科技,39(1):1-8.

馮士筰,1982.風暴潮導論.北京:科學出版社.

傅賜福,于福江,劉秋興,等.2013,四象限非對稱風場對風暴潮的改進研究.海洋通報,32(6):626-632.

傅賜福,于福江,王培濤,等,2013.濱海新區溫帶風暴潮災害風險評估研究.海洋學報,35(1):55-62.

夏波,張慶河,蔣昌波,2012.基于非結構網格的波流耦合數值模式研究.海洋與湖沼,44(6):1451-1456.

夏波,張慶河,楊華,2006.水動力時空變化對近岸風浪演化的影響--以渤海灣西南岸為例.海洋通報,25(5):1-8.

楊支中,沙文鈺,朱首賢,等,2005.一種新型的非對稱臺風海面氣壓場合風場模型.海洋通報,24(1):62-68.

葉鳳娟,李玉杰,牛福新,等,2012.風暴潮災害對天津濱海新區的影響分析.天津航海,123(2):38-39.

鄭立松,余錫平,2010.波浪輻射應力對風暴增水的影響研究.南水北調與水利科技,8(4):74-78.

(本文編輯:袁澤軼)

Application of the surface wave and storm surge coupled model in Tianjin coastal areas

LI Xue1,2,WANG Zhi-feng2,WU Shuang-quan3,DONG Sheng2, JIA Jing2,ZHANG Xiao-shuang4,WU Hao2
(1.Key Laboratory of Digital Ocean,State OceanicAdministration,Tianjin 300171,China;2.College of Engineering,OceanUniversity of China,Qingdao266100,China;3.National Marine Data and Information Service,Tianjin 300171,China;4.KeyLaboratoryof Marine Environmental Information Technology,National MarineData and Information Service,Tianjin300171,China)

In order to reasonably describe the coastal wave and the characteristics of storm surge in Tianjin coastal area, based on an unstructured triangular grid,this paper established a refined coastal wave-storm coupled simulation model. Fujita typhoon wind model was adopted as the wind field,and coastal wave selected SWAN wave model,and storm surge by ADCIRC model.By verifications with numerical simulations of three typical typhoon storms,the computational results of wind speed,significant wave height and storm surge were agreeable with measured data,approving that this coupled model can reflect the coastal wave and the characteristics of storm during a typhoon process in Tianjin coastal area,so as to provide a scientific basis for disaster prevention and reduction in Tianjin.

typhoon wave;storm surge;SWAN;ADCIRC

P731.22

A

1001-6932(2016)06-0657-09

10.11840/j.issn.1001-6392.2016.06.008

2015-04-24;

2015-09-14

國家海洋局數字海洋科學技術重點實驗室開放基金(KLDO201406);山東省海洋工程重點實驗室開放基金(201362045);國家自然科學基金(51279186;51479813;51509226);中央高校基本科研業務費專項(201513040)。

李雪(1990-),女,在讀博士,主要從事海洋環境數值預報研究。電子郵箱:lixueouc@gmail.com。

王智峰,博士,講師。電子郵箱:wzf1984@ouc.edu.cn。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 午夜a级毛片| 国内精自线i品一区202| 精品少妇三级亚洲| 国产91色| 欧美日韩国产在线观看一区二区三区| 另类专区亚洲| 国产成人精品三级| 制服丝袜亚洲| 国产欧美视频在线观看| 四虎影视国产精品| 欧美日韩国产成人高清视频| 欧洲熟妇精品视频| 国产精品夜夜嗨视频免费视频| 操国产美女| 精品国产91爱| 99视频只有精品| 国产欧美成人不卡视频| 国产精品 欧美激情 在线播放 | 88av在线| 久久精品娱乐亚洲领先| 欧美特级AAAAAA视频免费观看| 成人国产小视频| 强奷白丝美女在线观看| 亚洲精品成人7777在线观看| 2021国产在线视频| 一级毛片基地| 国产99在线| 天天综合网站| 国产一级片网址| 国产成人精品一区二区三在线观看| 超碰精品无码一区二区| 国产精品一区二区国产主播| 视频一区亚洲| 欧美性精品| 国产成人AV综合久久| 国产凹凸一区在线观看视频| 午夜啪啪网| 五月激情综合网| 狠狠做深爱婷婷综合一区| 国产va免费精品| 日韩区欧美区| 国产成人免费高清AⅤ| 91日本在线观看亚洲精品| 黄色网在线| 婷婷六月激情综合一区| 97se亚洲综合不卡 | a级免费视频| 国产精品久久久久久久伊一| 久久午夜影院| 少妇被粗大的猛烈进出免费视频| 亚欧美国产综合| 又粗又硬又大又爽免费视频播放| 91小视频版在线观看www| 试看120秒男女啪啪免费| 日本三区视频| 亚洲综合中文字幕国产精品欧美| 婷婷开心中文字幕| 国产无码在线调教| 天天综合天天综合| 露脸一二三区国语对白| 毛片在线播放a| 国产91精品调教在线播放| 在线永久免费观看的毛片| 欧美午夜久久| 美美女高清毛片视频免费观看| 奇米精品一区二区三区在线观看| 免费观看男人免费桶女人视频| 欧美日韩一区二区三| 国产麻豆aⅴ精品无码| 欧美 亚洲 日韩 国产| 国产女人水多毛片18| 国产亚洲成AⅤ人片在线观看| 亚洲视频免| 制服丝袜国产精品| 国产福利在线免费| 日韩麻豆小视频| 妇女自拍偷自拍亚洲精品| 久久99久久无码毛片一区二区| 中文字幕免费在线视频| 久久久亚洲色| 免费国产黄线在线观看| 一级毛片a女人刺激视频免费|