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

基于FVCOM的象山港海域潮汐潮流與溫鹽結構特征數值模擬

2014-06-19 06:21:34韓松林梁書秀孫昭晨
水道港口 2014年5期

韓松林,梁書秀,孫昭晨

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

基于FVCOM的象山港海域潮汐潮流與溫鹽結構特征數值模擬

韓松林,梁書秀,孫昭晨

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

摘要:基于有限體積法的FVCOM模型,建立了象山港海域的三維潮汐潮流和溫鹽數值模型,計算中考慮了潮流、風、太陽輻射和徑流因素的影響。模擬結果與2009年的監測資料進行了對比驗證,結果表明建立的模型可以模擬該海域的水流運動和溫鹽分布特征。通過對數值結果分析得到了該海域的同潮圖、潮流橢圓圖、潮流性質和溫鹽分布等。結果表明,象山港的潮汐屬于非正規半日淺海潮;M2分潮流橢圓長軸從口門到灣頂逐漸減小,其走向與岸線的方向基本一致;狹灣內呈現往復流特征而口門外開闊海域呈旋轉流特征。灣口和灣頂部有著顯著的溫度差和鹽度差,海水溫度由灣口向灣頂部逐漸增大,鹽度分布則正好相反。狹灣內距離灣口不同位置的橫向溫度、鹽度垂向分布結構特征各不相同。

關鍵詞:FVCOM;潮流;溫鹽;數值模擬;象山港

象山港處于浙江北部沿海,北面緊靠杭州灣,南臨三門灣,東側為舟山群島,縱深60 km,口門寬約20 km,港內寬3~8 km,是一個典型的狹長形半封閉海灣[1]。港內自然環境優越,水產資源豐富,具有良好的港口資源、濕地資源和海洋生物等資源,是寧波市發展海洋經濟的重要區域。近20 a來,港內臨港工農業、水產養殖業以及旅游業等發展迅速,極大地提升了區域經濟發展水平。但是隨著海洋經濟的迅猛發展,導致沿岸的工農業開發加劇、污水的過度排放和養殖業的不合理發展等,使象山港的海洋環境和生態系統受到了嚴重的威脅。灣內水質富營養化嚴重,赤潮時有發生。潮汐潮流運動是象山港的主要動力過程,是污染物、營養鹽等輸運的主要動力因素;溫鹽分布是海洋生態系統中重要的環境因子,是影響浮游動物、浮游植物生長生存的重要影響因素。因此建立象山港海域的潮汐潮流和溫鹽模型,研究其分布特征,對海洋資源的開發和生態系統的保護以及進一步對生態系統的研究等具有重要的意義。

FVCOM模型是由UMASSD?WHOI聯合開發的一個基于非結構網格的、有限體積的、三維原始方程的海洋模型[2]。模型水平方向上采用非結構化的三角形網格,可對地形復雜的區域局部加密,以更好的模擬復雜的岸線以及島嶼;垂向采用σ坐標,有助于處理不規則的海底地形。在潮間帶,模型采用干濕網格技術考慮了潮灘對潮流的影響。因此,利用FVCOM模型非常有利于模擬具有復雜岸線和地形的象山港海域。對于象山港的潮流潮汐和溫鹽特征已有學者進行了一些研究,董禮先等[3]建立了潮波運動數值模型,研究了象山港內影響M4分潮的控制因子和機理;曹穎等[4]基于FVCOM模擬了溫排水的擴散輸運過程,朱軍政等[5]應用該模型模擬了潮流鹽度的時空分布;關于象山港海域水動力和溫鹽結構觀測資料分析的研究較少,董禮先等[6]利用1981~1990年的實測水文資料象山港內的鹽度分布和環流結構。這些研究使我們對象山港的水動力及溫鹽特征有一定的了解,但是對于象山港海域考慮多種驅動因素下的水動力特征和溫鹽結構分布展開的分析研究還相對較少,而水動力和溫鹽的模擬是建立生態系統模型的基礎,因此有必要對其進行進一步研究。

本文采用非結構化的三角形網格,建立了象山港海域的三維潮流和溫鹽模型,綜合考慮潮流、風、太陽輻射和徑流等因素作用的影響,對象山港海域潮汐、潮流和溫鹽分布進行了數值模擬。通過對數值結果分析給出了主要分潮的同潮圖和潮流橢圓圖,溫度、鹽度的表、底層分布及狹灣內橫向、縱向的溫鹽垂向分布,并結合近期現場實測潮位、流速和溫鹽數據資料討論了它們的分布特征,為進一步了解象山港海域的水動力以及溫鹽分布特征提供參考,也為下一步對該海域物質輸運及生態系統研究奠定基礎。

1 計算模型

1.1模型方程

模型控制方程包括連續性方程、動量方程、溫鹽擴散方程和密度方程。σ坐標下的控制方程如下

式中:t是時間;D是總水深;ζ為水位;u,v,ω分別為σ坐標下x,y和σ方向的速度分量,σ取值從海底處的-1到海表面處的0;g為重力加速度;ρ為海水密度;ρo為水體參考密度;f為科氏參數;Fx,Fy,FT和FS分別是水平向動量、溫度和鹽度擴散項。T為海水溫度,S為海水鹽度;H?是水體吸收的太陽輻射;Km和Kh分別為垂向紊動粘性系數和熱擴散系數,模型中國采用修正的MY-2.5湍流閉合模型[7]求解。

海表面和底部速度邊界條件為

式中:(τsx,τsy)和(τbx,τby)分別為海表面風切應力和海底摩擦切應力在x和y方向上的分量。

海表面和底部溫鹽邊界條件為

式中:Qn(x,y,t)為表面凈熱通量,包括短波輻射、長波輻射、感熱和潛熱通量四部分;SW(x,y,0,t)是海表面處短波輻射通量;cp為海水比熱系數;P∧和E∧分別為降雨和蒸發量;AH為水平熱擴散系數;α是海底地形的坡度;n為垂直于坡度軸線方向。

模型采用模式分裂法求解,二維外膜中的控制方程在三角形單元積分后,通過改進的四階龍格庫塔進行求解;三維內膜的動量方程求解采用一種簡單的顯式和隱式相結合的方法,其中速度的局地變化采用一階精度的前差格式積分。對流項采用二階精度的迎風顯式格式求解,垂向擴散項采用隱式求解。具體離散求解可以參考文獻[8]。

1.2模型設置

模型計算范圍為象山港全域(29°24′~29°50′N,121°23′~122°05′E),包括象山港狹灣、佛渡水道和牛鼻山水道,如圖1所示。開邊界取為兩段,分別為灣口東北側的郭巨鎮—六橫島連線和東南側的爵溪鎮—臺門連線。象山港灣內岸線曲折,島嶼眾多,為了更好地擬合復雜的地形條件,對島嶼和岸線以及水深變化劇烈的區域進行了加密,計算網格見圖2。水平分辨率在島嶼和岸線周圍為約80 m,在灣口部為400 m左右,水平網格單元數為18 645,網格節點數為9 958,垂直向均勻分為11個σ層。開邊界水位由杭州灣大模型提供[9],該模型邊界由實測潮位提供,并在大范圍內驗證良好。模型中考慮風、徑流、潮流和密度流的共同作用,其中風、空氣溫度、相對濕度、凈短波輻射采用NCEP每隔6 h平均的再分析資料,凈熱通量采用文獻[10]中公式計算。徑流考慮了鳧溪河、顏公河和墻頭排水口等淡水的注入,流量采用文獻[11]中的徑流量數據,入海口位置見圖1。

圖1 計算區域及站位分布圖Fig.1Calculated area and location of stations

圖2 數值計算網格劃分Fig.2Gird of calculated area

2 結果分析與討論

2.1潮汐潮流特性分析

采用2009年6~7月烏沙山站連續15 d的實測潮位資料和0916、0917、0918、0919四個測流站(圖1)大潮期和小潮期連續25 h實測海流資料與模型計算結果進行了比較,以驗證模型的精度及可靠性。圖3和圖4分別給出潮位、潮流流速流向的驗證結果。從圖中可以看出,烏沙山潮位的計算值與實測值吻合程度較好;除個別時刻外,計算得到的流速和流向與實測值相對誤差較小。流速沿水深遞減,表層流速最大,底層最小,底層流速約為表層的50%。

圖3 烏沙山站模型模擬潮位與實測值對比Fig.3Comparison of observed and simulated tidal level at Wushashan

圖40917 、0919測站模型模擬結果與實測值對比Fig.4Comparison of observed and simulated tidal speed and direction at 0917 and 0919

由潮位計算結果經潮汐調和分析得到象山港M2分潮的同潮圖(圖5)。潮汐的類型可以通過不同分潮的振幅比和判斷。象山港海域F均大于0.5,除靠近外海的牛鼻山水道外,G均大于0.04,屬于非正規半日淺海潮。M2分潮在潮位中占主導地位,其次是S2分潮;狹灣內的M2分潮振幅HM2均在1.3 m以上,由灣口的1.3 m逐漸增到灣頂部的1.7 m,灣口M2分潮的遲角與灣頂只差4°;狹灣內淺水分潮作用明顯,M4分潮振幅HM4由從口門處的0.1 m增到灣頂的0.5 m。潮波自口門傳入后,由于不斷受地形及邊界的反射作用,逐漸由前進潮波轉為駐波性質;潮差由灣口向灣頂逐漸增大,港頂平均潮差可達3.7 m。受淺海分潮影響,狹灣內漲落潮歷時不對稱明顯,漲潮歷時均大于落潮歷時,越往港內漲潮歷時越長,港頂最大漲落潮時差約170 min。

利用潮流調和分析結果,根據公式計算可知象山港海域的潮流屬半日潮流區,淺海分潮明顯。象山港海域M2分潮流在潮流中占主導地位,圖6為該海域M2分潮潮流橢圓圖。M2分潮流橢圓長軸由牛鼻山水道至狹灣口門之間約為1.0~1.3 m/s,由狹灣口門至灣頂,M2分潮流逐漸減弱,其橢圓長軸在西滬港口和鐵港口附近分別為0.7 m/s和0.4 m/s左右,走向與岸線的方向基本一致。從調和分析結果可得,M2分潮流的橢圓率在狹灣口門外海域介于0.5~0.8,潮流運動形式呈旋轉流特征;狹灣內橢圓率在0~0.1,呈往復流特征,這主要是受灣內地形和邊界的制約。灣內潮流的流向與岸線基本平行,落潮流速大于漲潮流速,最大流速發生在狹灣口門附近,狹灣內的流速由灣口向灣頂逐漸減小,灣內漲、落潮流最大實測流速分別可達1.54 m/s和1.83 m/s[1]。漲潮流通過牛鼻山水道和佛渡水道傳入象山港港口門后,沿著主水道向灣內傳播,至西滬港口門處分出一支傳入港內,在缸爿山附近由于水道狹窄,水流在此比較湍急,流速可達到1 m/ s以上;在烏沙山附近,受白石山島—中央山島—銅山島一線島嶼的影響,水流在此分為南、北兩路;過島后匯合繼續向西推進,在強蛟處又分成兩支,分別進入黃墩港和鐵港。落潮時,流路與漲潮過程相反,3個支港內的水體匯入主水道后沿原路退出象山港。灣內3個支港西滬港、鐵港、黃墩港海域在低潮位時絕大部分潮灘露出,高潮位時淹沒。

2.2溫鹽結構分析

圖5 象山港M2分潮同潮圖Fig.5Calculated co?tidal chart in Xiangshan Bay

圖6 象山港表層M2分潮潮流橢圓分布Fig.6Calculated surface component tidal current ellipses in Xiangshan Bay

圖7 模擬表層溫度、鹽度與實測值對比Fig.7Comparison of observed and simulated temperature and salinity

溫度和鹽度采用2009年象山港赤潮監控區水環境監測數據進行驗證,監測時間段為4~10月,監測頻率為半月一次。圖7為不同測站海水表層溫度、鹽度模擬值與實測值的對比。象山港海域水溫隨季節變化顯著,春季水溫平均值為15.6℃,夏季水溫平均值為29.1℃,最高實測溫度為32.1℃,發生在7月下旬,秋季水溫平均值為23.6℃。象山港鹽度隨月份變化明顯,4月份鹽度平均值為23.7,7月份平均鹽度為28.3,10月份鹽度平均值為23.6。7月份鹽度偏高的原因是由于夏季長江淡水向東北方向擴展,灣內鹽度受臺灣暖流等外海高鹽水控制所致[6]。8月份受降水和徑流量的影響,鹽度出現明顯降低。圖8為港內不同站位7月15日溫度、鹽度模擬值與實測值對比,可以看出從空間分布上,灣口和灣頂部有著顯著的溫度差和鹽度差。象山港海水溫度由灣口向灣頂部逐漸增大,港口和港頂的溫度差為4.4℃;與溫度分布相反,海水鹽度呈灣口部高,由灣口部到灣口部逐漸降低的趨勢,灣口和灣頂的鹽度差為3.2。總體來看,模型模擬結果與實測資料吻合較好,基本反映了象山港海域年度溫度、鹽度分布變化規律,說明本文模型設置和參數的選取是合適的。

圖8 不同站位7月15日表層溫度、鹽度對比Fig.8Comparison of observed and simulated temperature and salinity at different stations on July 15

(1)溫鹽的平面分布。圖9-a和圖9-b是計算海區夏季表、底層平均水溫分布。由圖9可以看出,象山港表層水溫在27~29℃,灣頂水溫最高,從灣頂往外水溫逐漸下降,而且等值線在口門外較灣內密,溫差較大。底層水溫分布與表層基本相似,較表層普遍低0.5℃左右。等溫線在灣內呈拋物線型,離兩岸陸地越近,溫度越高,這是因為離陸地近的地方水深較淺,水體升溫快,中間深槽內溫度較低,且在底層表現尤為明顯。

圖9 夏季7月份模擬溫度、鹽度平面分布Fig.9Distribution of temperature and salinity in July

圖9-c和圖9-d為計算海區夏季平面表、底層平均鹽度分布,鹽度平面分布的總體特征是牛鼻山水道和佛渡水道海域鹽度最高,鹽度值在30左右,從狹灣口向灣頂逐漸遞減,灣頂部鹽度為23~26。鹽度的垂向分布為表層最低,底層略高一點。鹽度等值線的主要特征為在表層向灣口凸出,而底層向灣頂凸出,這種分布特征與象山港的余流特征是一致的,象山港狹灣內的表層余流方向指向灣口而底層余流方向指向灣頂[6]。從鹽度大面分布的等值線可以看出,多數時間港內鹽度北岸較南岸稍高,這是由于潮波傳播在北岸較南岸快(圖5)。

圖10 象山港溫度和鹽度縱向斷面分布Fig.10Longitudinal distribution of temperature and salinity in Xiangshan Bay

圖11 象山港溫度和鹽度橫向斷面分布Fig.11Latitudinal distribution of temperature and salinity in Xiangshan Bay

(2)溫鹽的垂向分布。圖10為7月份沿象山港縱向的溫度和鹽度垂向分布,斷面的位置如圖10-a所示。水溫沿縱向斷面在垂向上分層明顯,表層高、底層低,表底層溫度差一般在0.4~0.6℃;溫度從灣頂到灣口逐漸遞減,與平面分布一致。鹽度垂向則分層較明顯,表層、底層鹽度差為0.5~1。從圖10可以看出,在灣口和灣頂附近鹽度等值線較密,在垂向上均為表層向灣口、底層向灣頂的傾斜狀,與董禮先等[6]研究中實測數據得到的鹽度縱向剖面分布基本一致。

圖11給出狹灣內距離港頂不同距離的橫向斷面的溫度和鹽度分布。從模擬結果來看,狹灣內不同位置的溫度、鹽度垂向結構不同。B-C斷面位于狹灣外段,在4 m以上水層,水溫水平方向較為均勻,4 m以下水層,水溫在兩側高,中間深槽內低;鹽度分布表現為底層高于表層,在斷面中部表層局部鹽度較低,表、底鹽度差約為1.0,在水深6 m以下水層,靠近北岸一側存在有一股高鹽水,與文獻[6]中實測潮均鹽度分布相似。D-E斷面位于狹灣中段,水溫在4 m以下水層表現為北岸低南岸高,鹽度表現為北岸高南岸低,與平面分布一致。F-G斷面位于狹灣內段,水溫的分布與B-C斷面相似;鹽度分布表現為在水深4 m以上,兩側鹽度高,中間鹽度低,在4 m水深以下,鹽度呈現兩側低,中間深槽內高的特征。

3 結論

平面采用非結構化網格,基于有限體積法的FVCOM海洋模型建立了象山港三維潮汐、潮流和溫鹽數值模型,模擬了象山港內三維水動力過程和溫鹽的季節變化。通過將模擬結果與實測潮位、潮流和溫鹽資料的比較,驗證了模型的有效性和精度,說明模型能夠較好地反映了象山港潮汐、潮流的特征以及溫鹽的變化規律及分布特征。通過對數值模擬結果進行詳細分析,得到如下主要結論:

(1)象山港的潮汐類型屬于非正規半日淺海潮。狹灣內的M2分潮振幅由口門處的1.3 m逐漸增到灣頂部的1.7 m,淺水分潮作用明顯,M4分潮振幅由從口門處的0.1 m增到到港頂的0.5 m。狹灣內漲落潮歷時不對稱明顯,漲潮歷時均大于落潮歷時,越往港內漲潮歷時越長。

(2)由狹灣口門至灣頂,M2分潮流逐漸減弱,其橢圓長軸在灣口、西滬港口和鐵港口附近分別為1.2 m/s、0.7 m/s和0.4 m/s左右,走向與岸線的方向基本一致。M2分潮流的橢圓率在狹灣口門外海域介于0.5~0.8,狹灣內則在0~0.1,因此灣內為往復流性質,口門外開闊海域呈旋轉流特征。

(3)灣口和灣頂部有著顯著的溫度差和鹽度差。象山港海水溫度由灣口向灣頂部逐漸增大,港口和港頂的溫度差為4.4℃;與溫度分布相反,海水鹽度呈灣口部高,由灣口部到灣頂部逐漸降低的趨勢,灣口和灣頂的鹽度差為3.2。由于狹灣灣內水深地形以及寬度差別明顯,距離灣口不同位置的橫向斷面溫度、鹽度分布結構特征明顯不同。

對象山港水動力場和溫鹽場的準確計算,說明該模型能夠適用于象山港復雜地形的水動力的數值模擬,為進一步對該海域物質輸運及生態系統研究奠定基礎。

致謝:本文采用的實測潮位潮流和溫鹽資料均來源于寧波市海洋環境監測中心2009年的海洋水文調查,特致感謝。

參考文獻:

[1]中國海灣志編輯委員會.中國海灣志:第五分冊[M].北京:海洋出版社,1992.

[2]CHEN C S,LIU H D,Beardsley R C.An unstructured grid,finite?volume,three?dimensional primitive equations ocean model:ap?plication to coastal ocean and estuaries[J].Journal of Atmospheric and Ocean Technology,2003,20:159-186.

[3]董禮先,蘇紀蘭.象山港潮波響應和變形研究,II.象山港潮波數值研究[J].海洋學報,1999,21(2):1-8. DONG L X,SU J L.Tide response and tide wave distortion study in the Xiangshan Bay,II.Numerical modelling study in the Xiang?shan Bay[J].Acta Oceanologica Sinica,1992,21(2):1-8.

[4]曹穎,朱軍政.基于FVCOM模式的溫排水三維數值模擬研究[J].水動力學研究與進展:A輯,2009,24(4):432-439. CAO Y,ZHU J Z.Numerical simulation of 3D cooling water based on FVCOM[J].Journal of Hydrodynamics,2009,24(4):432-439.

[5]朱軍政,曹穎.FVCOM模型在象山港三維潮流鹽度計算中的應用[J].海洋環境科學,2010,29(6):899-903. ZHU J Z,CAO Y.Application of FVCOM for computation of 3D tidal flow and salinity in Xiangshan Bay[J].Marine Environmental Science,2010,29(6):899-903.

[6]董禮先,蘇紀蘭.象山港鹽度分布和水體混合,I.鹽度分布和環流結構[J].海洋與湖沼,2000,31(2):151-158. DONG L X,SU J L.Salinity distribution and mixing in Xiangshan Bay,I.Salinity distribution and circulation pattern[J].Oceanolo?gia et Limnologia Sinaca,2000,31(2):151-158.

[7]Mellor G L,Yamada T.Development of a turbulence closure model for geophysical fluid problem[J].Rev.Geophys.Space.Phys.,1982,20:851-875.

[8]Chen C S,Beardsley R C,Cowles G.An unstructured grid,finite?volume coastal ocean model FVCOM user manual[R].USA:Uni?versity of Massachusetts?Dartmouth,2006.

[9]熊偉,劉必勁,孫昭晨,等,寧波舟山近海三維潮汐潮流數值模擬[J].水道港口,2011,32(6):399-407. XIONG W,LIU B J,SUN Z C,et al.3D numerical simulation of tide and tidal currents in sea adjacent to Ningbo and Zhoushan[J].Journal of Waterway and Harbor,2011,32(6):399-407.

[10]Ahsan A K,Blumberg A F.Three?dimensional hydrothermal model of Onondaga Lake,New York[J].Journal of Hydraulic Engi?neering,ASCE,1999,125(9):912-923.

[11]陳偉,王小華.象山港鹽度鋒面及其動力成因[J].東海海洋,1998,16(1):1-9. CHEN W,WANG X H.Tidal fronts and their formation in the Xiangshan Bay[J].DongHai Marine Science,1998,16(1):1-9.

交通運輸部將加快五大國際航運中心建設

2014年9月3日發布的《國務院關于促進海運業健康發展的若干意見》指出,要提升海運業國際競爭力。引導要素和產業集聚,加快建設國際海運交易和定價中心,打造國際航運中心。據此,交通運輸部將在若干意見的框架下,推動現代航運服務業發展指導意見的出臺,加快國際航運中心建設,包括上海、天津濱海新區、大連、武漢和重慶。(殷缶,梅深)

Biography:HAN Song?lin(1986-),male,doctor student.

中圖分類號:TV 143;O 242.1

文獻標識碼:A

文章編號:1005-8443(2014)05-0481-08

收稿日期:2013-11-04;修回日期:2013-12-26

基金項目:國家海洋局海洋公益性行業科研專項(201105009);國家自然科學基金(51279028)

作者簡介:韓松林(1986-),男,河南省林州市人,博士研究生,主要從事近海水動力和物質輸運研究。

Numerical simulation of tides,tidal currents and temperature?salinity structures in Xiangshan Bay based on FVCOM

HAN Song?lin,LIANG Shu?xiu,SUN Zhao?chen
(State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China)

Abstract:Based on the unstructured grid,finite?volume coastal ocean model(FVCOM),the hydrodynamic and temperature?salinity numerical model were established in Xiangshan Bay.The tidal current,wind,the solar radia?tion and river discharge were considered in the model.The comparison of the simulated results with the measured data about tide,tidal current,temperature and salinity showed that the model could simulate the hydrodynamic and the distribution of temperature?salinity structures.The co?tidal chart,the component tidal current ellipses were ob?tained from the result.The results show that the tide of Xiangshan Bay is mainly irregular semidiurnal shallow tide. The major axis of M2tidal current component ellipse decreases from the mouth of the bay to the top and the direc?tion is parallel to the coastal line.The tidal current is rectilinear current in the fjord and rotary current at the outside of the bay mouth area.There are obvious temperature difference and salinity difference between the mouth and head of the bay.The temperature increases from the bay mouth to bay head and the salinity is just opposite.The vertical profile characteristic of temperature and salinity at different location in the fjord is different.

Key words:FVCOM;tidal current;temperature and salinity;numerical simulation;Xiangshan Bay

主站蜘蛛池模板: 日韩欧美国产成人| 日本欧美视频在线观看| 欧美亚洲香蕉| 久久久亚洲色| 热99精品视频| 亚洲欧洲一区二区三区| 69av在线| 色屁屁一区二区三区视频国产| 青青草原国产免费av观看| 天天综合色网| 999精品视频在线| 免费国产在线精品一区| 亚洲国产成人在线| 国产精品天干天干在线观看| 激情综合婷婷丁香五月尤物| 日韩av无码精品专区| 有专无码视频| 麻豆AV网站免费进入| 超薄丝袜足j国产在线视频| 熟女日韩精品2区| 亚洲天堂网2014| 全部免费毛片免费播放 | 无码精油按摩潮喷在线播放 | 欧美成人综合视频| 久久精品视频亚洲| 亚洲精品色AV无码看| 亚洲三级视频在线观看| 亚洲资源在线视频| 国产中文一区a级毛片视频| 国产肉感大码AV无码| 成人中文在线| 国产精品成人观看视频国产 | 中文字幕在线永久在线视频2020| 免费看美女自慰的网站| 99精品一区二区免费视频| 亚洲精品自在线拍| 亚洲手机在线| 国产日韩欧美精品区性色| 蜜臀AV在线播放| 超薄丝袜足j国产在线视频| 99r在线精品视频在线播放| 欧美精品二区| 国产日韩AV高潮在线| 亚洲人免费视频| 免费a级毛片视频| 国产无套粉嫩白浆| 日韩无码视频网站| 国产日产欧美精品| 国产精女同一区二区三区久| 国产精品无码AⅤ在线观看播放| 久久综合激情网| 91免费观看视频| 日韩福利视频导航| 国产精品女主播| 97精品国产高清久久久久蜜芽| 69国产精品视频免费| 国产精品爽爽va在线无码观看| 亚洲成人网在线播放| 国产青青操| 精品91视频| 日韩欧美国产精品| 亚洲熟女中文字幕男人总站| 99re热精品视频国产免费| 麻豆AV网站免费进入| 日韩成人在线网站| 最新国语自产精品视频在| 国产自产视频一区二区三区| 久久久久久久蜜桃| 久久青草免费91线频观看不卡| 国产精品一区二区国产主播| 久久久久久尹人网香蕉 | 国产精品粉嫩| 一区二区理伦视频| 国产大全韩国亚洲一区二区三区| 香蕉在线视频网站| av手机版在线播放| 亚洲一级毛片在线观| 久久无码av一区二区三区| 无码高潮喷水专区久久| 久久久精品无码一二三区| 成人一区在线| 香蕉eeww99国产在线观看|