亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        渤海風(fēng)暴減水特征及其對深水航路影響的數(shù)值模擬

        2014-06-05 14:35:01傅賜福付翔吳少華于福江2董劍希
        海洋學(xué)報(bào) 2014年3期
        關(guān)鍵詞:風(fēng)暴潮潮位航路

        傅賜福,付翔,吳少華,于福江2,董劍希

        (1.國家海洋環(huán)境預(yù)報(bào)中心,北京 100081;2.國家海洋局海洋災(zāi)害預(yù)報(bào)技術(shù)研究重點(diǎn)實(shí)驗(yàn)室,北京 100081)

        渤海風(fēng)暴減水特征及其對深水航路影響的數(shù)值模擬

        傅賜福1,付翔1,吳少華1,于福江12,董劍希1

        (1.國家海洋環(huán)境預(yù)報(bào)中心,北京 100081;2.國家海洋局海洋災(zāi)害預(yù)報(bào)技術(shù)研究重點(diǎn)實(shí)驗(yàn)室,北京 100081)

        利用環(huán)渤海9個(gè)沿岸站近10 a潮位資料分析渤海海域的風(fēng)暴減水特征,結(jié)果表明:渤海年均出現(xiàn)50 cm和100 cm風(fēng)暴減水分別超過30 d和6 d,每年的9月至翌年4月份風(fēng)暴減水最為頻繁;建立了一套精細(xì)化天文潮-風(fēng)暴潮耦合模型用于渤海深水航路的潮位預(yù)報(bào),各站天文潮模擬驗(yàn)證的平均均方根誤差為18.5 cm,由此計(jì)算得到航路代表點(diǎn)的潮汐特征值并作潮汐預(yù)報(bào);后報(bào)模擬了近10 a重大風(fēng)暴減水過程,模擬與實(shí)測吻合較好,說明該耦合模型可為該航路的潮位預(yù)報(bào)提供有益參考。

        風(fēng)暴減水;深水航路;高分辨率;風(fēng)暴潮-天文潮耦合;數(shù)值模擬

        1 引言

        渤海包括遼東灣、渤海灣、萊州灣,由于其半封閉地形緣故,且位于中緯度地區(qū),一年四季均有風(fēng)暴潮發(fā)生,歷史上曾多次發(fā)生強(qiáng)風(fēng)暴潮災(zāi)害,尤其以由西風(fēng)帶天氣系統(tǒng)引起的溫帶風(fēng)暴潮最為頻繁和嚴(yán)重[1—2]。渤海屬超淺海,平均水深不足20 m,當(dāng)離岸風(fēng)長時(shí)間吹刮海面,較大風(fēng)暴減水與天文低潮等相疊加會(huì)引起的異常低潮位,將嚴(yán)重影響了大型油輪及巨型貨輪的正常航行和錨泊。

        從老鐵山水道航行至曹妃甸錨地的其中一段航路水深約25 m(見圖3b,虛線框),該航路為大型貨船進(jìn)出渤海的深水航路,是連接大連港與曹妃甸港的唯一航路,為了滿足吃水較深的大型船只對該航路的航運(yùn)要求,需對該航路及附近海域的潮位及風(fēng)暴減水進(jìn)行深入的研究。

        國內(nèi)許多學(xué)者對風(fēng)暴減水的研究主要側(cè)重近岸工程的極低水位計(jì)算,分別基于構(gòu)造溫帶天氣系統(tǒng)進(jìn)行數(shù)值模擬[3]、統(tǒng)計(jì)分析[4]、重大減水過程分析[5]等思路進(jìn)行研究。本文則主要分析減水特征并利用經(jīng)檢驗(yàn)的天文潮-風(fēng)暴潮耦合模型進(jìn)行數(shù)值模擬研究,以期對深水航路的潮位預(yù)報(bào)有益。

        2 近10年渤海風(fēng)暴減水特征

        2.1 渤海風(fēng)暴減水天氣系統(tǒng)特征

        分析環(huán)渤海各站同時(shí)出現(xiàn)明顯風(fēng)暴減水過程,發(fā)現(xiàn)造成渤海風(fēng)暴減水的天氣系統(tǒng)主要由冷空氣造成。

        (1)冷空氣配合氣旋:由于冷空氣東移南下恰遇氣旋東移,渤海風(fēng)向變化由東北—北—西北,氣旋出海前的東風(fēng)使北黃海的一部分海水流入渤海,渤海風(fēng)暴潮情況呈現(xiàn)先增后減,并且此天氣形勢下的風(fēng)暴潮量值比較大。這種減水過程的特點(diǎn)是整個(gè)渤海的風(fēng)暴減水持續(xù)時(shí)間長,且遼東灣減水值略大于渤海灣、萊州灣,而渤海海峽外減水最小。

        (2)西路冷空氣南下:受西路空氣南下影響,渤海前期為偏北風(fēng),渤海北部先開始出現(xiàn)減水,隨著高壓南下東移,渤海風(fēng)向由北—西北—西,渤海西部和南部也開始減水,因此渤海西部、西南部的渤海灣和萊州灣能出現(xiàn)較大減水,減水值甚至大于渤海北部。

        2.2 近10年渤海風(fēng)暴減水特征

        利用環(huán)渤海9個(gè)沿岸驗(yàn)潮站(表1)近10 a的實(shí)測潮位數(shù)據(jù),進(jìn)行渤海減水特征分析。增減水分離采用觀測潮位直接減去周期性天文潮部分,得到逐時(shí)增減水值。

        表1 環(huán)渤海9個(gè)驗(yàn)潮站位置及數(shù)據(jù)信息

        統(tǒng)計(jì)各站大于50 cm和大于100 cm減水的天數(shù),分別做年平均和多年平均月際變化分析,見圖1和圖2。從年平均變化來看,減水大于50 cm:渤海灣、萊州灣一年中出現(xiàn)的天數(shù)最多,為38~50 d,遼東灣約35 d,渤海海峽口約15 d;減水大于100 cm:渤海灣、萊州灣6~8 d,遼東灣為6 d,渤海海峽口為1 d。總的來說,渤海海峽內(nèi)比海峽口更容易出現(xiàn)減水且渤海海峽內(nèi)的風(fēng)暴減水量值比海峽口偏大。

        圖1 各站年均大于50 cm(a)和100 cm(b)風(fēng)暴減水的天數(shù)

        從年平均月際變化來看,5—8月出現(xiàn)減水的情況較少,9月至翌年4月整個(gè)渤海減水頻繁。以各站的平均天數(shù)來看,11月至3月渤海每月大于50 cm的減水都超過3 d,12月最多,超過6 d;大于100 cm的減水12月和3月平均每月超過1 d。這是由于冬、春季冷空氣頻繁南下,渤海多受偏北風(fēng)影響,整個(gè)渤海易出現(xiàn)減水。

        3 天文潮-風(fēng)暴潮耦合模式的建立及參數(shù)設(shè)置

        由于風(fēng)暴潮影響范圍較廣,模型計(jì)算區(qū)域必須足夠大以刻畫大規(guī)模海水的運(yùn)動(dòng),又要考慮重點(diǎn)區(qū)域復(fù)雜地形,本文采用基于非結(jié)構(gòu)三角形網(wǎng)格的ADCirc模型,該模型是目前國際上較先進(jìn)的風(fēng)暴潮數(shù)值模式,由北卡羅來納大學(xué)的Luettich和美國圣母大學(xué)的Westerink教授于1992年研制后經(jīng)不斷完善發(fā)展。ADCirc是基于有限元方法、垂向平均二維、正壓的水動(dòng)力學(xué)模式,具有計(jì)算速度快,精確性和穩(wěn)定性高的特點(diǎn)[6—7]。

        圖2 環(huán)渤海各站大于50 cm(a)和大于100 cm(b)風(fēng)暴減水的月均天數(shù)

        3.1 控制方程及主要參數(shù)設(shè)置

        ADCirc模式在計(jì)算過程通過基于垂直平均的原始連續(xù)方程和運(yùn)動(dòng)方程來求解自由表面起伏、二維流速等3個(gè)變量,在運(yùn)動(dòng)方程中,除了考慮平流項(xiàng)、科氏力項(xiàng)、風(fēng)應(yīng)力項(xiàng)和底摩擦項(xiàng)外,還考慮了潮汐和側(cè)向黏性項(xiàng)等。在球坐標(biāo)系下,連續(xù)方程和運(yùn)動(dòng)方程表示為:

        式中,λ,φ為經(jīng)度和緯度;?為從海平面起算的自由表面高度;U,V為深度平均的海水水平流速;H=?+h為海水總水深;R為地球半徑;f=2Ωsinφ為科氏力參數(shù),Ω為地球自轉(zhuǎn)角速度;g為重力加速度;Ps為海水自由表面大氣壓;ρ0為海水密度;η為牛頓引潮勢;α為地球有效彈性因子;τsλ、τsφ為自由表面應(yīng)力;τbλ、τbφ為底摩擦應(yīng)力;Dλ、Dφ為動(dòng)量方程的水平擴(kuò)散項(xiàng)。

        初始條件為:ζ=u=v=0。

        海岸邊界條件:邊界的法向速度為0。

        開邊界條件:輻射邊界條件,文中由M2、S2、K2、 N2、2N2、L2、MU2、NU2、T2、K1、O1、P1、Q1、J1、M1、OO1等16個(gè)分潮驅(qū)動(dòng)計(jì)算,該16個(gè)分潮調(diào)和常數(shù)取自全球潮汐模型NAO99。

        模式中求解所需物理變量的過程中,空間采用有限元法離散,時(shí)間采用有限差分法。將連續(xù)方程和運(yùn)動(dòng)方程通過引入空間變量數(shù)值加權(quán)參數(shù)進(jìn)行結(jié)合求解,提高了計(jì)算結(jié)果的穩(wěn)定性,時(shí)間步長取為30 s,滿足CFL條件要求。

        底摩擦力τb與深度平均流呈二次平方律關(guān)系,底摩擦系數(shù)Cf采用線性與二次律混合形式,同時(shí)適用于大洋與近海,公式如下:

        式中,Cfmin為最小底摩擦因子;Hbreak為臨界水深;參數(shù)θ用來控制混合公式接近其上下限的快慢;參數(shù)λ描述摩擦因子隨水深增加而增大的快慢。

        海面風(fēng)應(yīng)力與風(fēng)速呈二次平方律關(guān)系,風(fēng)拖曳系數(shù)采用Garratt公式,Cd=0.001(0.75+0.667→W),(Cd的值在0.002 2和0.003之間)。

        文中所用風(fēng)場由美國研究機(jī)構(gòu)共同研究的新一代中尺度預(yù)報(bào)模式和同化系統(tǒng)WRF經(jīng)過GTS等資料同化后的再分析風(fēng)場提供,區(qū)域范圍5°~45°N、105°~155°E,空間分辨率0.1°×0.1°,約10 km[8]。

        3.2 網(wǎng)格剖分

        文中所用水深數(shù)據(jù)由3部分組成:(1)國家海洋環(huán)境預(yù)報(bào)中心業(yè)務(wù)化風(fēng)暴潮預(yù)報(bào)系統(tǒng)所用水深數(shù)據(jù),空間分辨率2′,這部分?jǐn)?shù)據(jù)主要用于黃海中部以南;(2)渤海海域、黃海北部海域局部精細(xì)化水深數(shù)據(jù),空間分辨率12″,這部分?jǐn)?shù)據(jù)主要用于渤海及黃海北部海域;(3)渤海海域精細(xì)化電子海圖1∶5 000,這部分?jǐn)?shù)據(jù)主要用于訂正修改航線及附近海域。這3類數(shù)據(jù)通過訂正、融合后,插值到網(wǎng)格格點(diǎn)。

        圖3 計(jì)算區(qū)域網(wǎng)格劃分(a)、渤海高分辨率水深(m)(b)及所涉及潮位站和航路代表點(diǎn)(c)

        計(jì)算區(qū)域包括渤海、黃海和東海北部海域,黃海網(wǎng)格分辨率2~5 km,渤海網(wǎng)格分辨率0.5~2 km,該套網(wǎng)格包括92 304個(gè)三角形單元,47 167個(gè)網(wǎng)格點(diǎn)。渤海及航路的網(wǎng)格剖分、水深分布以及所涉及潮位站和航路代表點(diǎn)位置見圖3(文中水深及潮位值均換算至平均海平面)。通過多種水深數(shù)據(jù)訂正融合,插值后的網(wǎng)格水深更加準(zhǔn)確,航線沿線的水深在24.5~28 m之間,基本反映了精細(xì)化電子海圖的航路信息。

        4 風(fēng)暴潮過程模擬驗(yàn)證

        4.1 天文潮模擬驗(yàn)證

        模擬了2010年12月1日至2011年12月31日,共396 d,提取2011年1月1日至12月31日全年的潮汐逐時(shí)值作為模擬值,各站潮汐調(diào)和分析值由國家海洋環(huán)境預(yù)報(bào)中心根據(jù)前一年的調(diào)和分析計(jì)算給出。圖4為曹妃甸、大沽燈塔、龍口、老虎灘4個(gè)站2011年8月份天文潮模擬值與調(diào)和分析值對比,表2為各驗(yàn)潮站逐時(shí)天文潮模擬值與調(diào)和分析值的均方根誤差。可以看出,在渤海大部分驗(yàn)潮站的均方根誤差普遍小于20 cm,說明所建立的模式能夠很好刻畫渤海的潮汐特征。

        另外,在這重點(diǎn)闡述大沽燈塔的天文潮模擬情況。與其他沿岸站不同的是,大沽燈塔是位于航道上的潮位站,所處航道水深約18 m,本文認(rèn)為大沽燈塔模擬的結(jié)果可以很大程度上反映該航路上其他點(diǎn)的模擬質(zhì)量,大沽燈塔所處海域平均潮差較大(約3 m),均方根誤差為18.49 cm,吻合較好,可以從一定程度上說明航路各點(diǎn)模擬天文潮與實(shí)測相差較小,是可信的。表3為航路各點(diǎn)計(jì)算得到的潮汐特征值,可以看出航路自東向西的潮汐類型由正規(guī)半日潮向不正規(guī)半日潮過渡,平均潮差約55~80 cm,由東向西逐漸減小。

        圖4 曹妃甸(a)、大沽燈塔(b)、龍口(c)、老虎灘(d)2011年8月份天文潮對比

        表2 各驗(yàn)潮站逐時(shí)天文潮模擬值與調(diào)和分析值的均方根誤差

        表3 航線4個(gè)代表點(diǎn)潮汐特征值

        4.2 風(fēng)暴減水過程模擬驗(yàn)證

        2005年12月21—22日,受西路冷空氣南下的影響,渤海受7~8級西-西北風(fēng)控制,渤海沿岸各站出現(xiàn)了較大減水。圖5為該過程21日14時(shí)、22日08時(shí)東亞地面天氣圖,圖6為該過程葫蘆島、塘沽風(fēng)暴潮、潮位模擬與實(shí)測對比圖。

        2007年3月4—6日,受西北路冷空氣南下和江淮氣旋出海的共同影響,渤海海域發(fā)生了一次強(qiáng)風(fēng)暴潮過程,5日清晨開始,由于江淮氣旋已經(jīng)出海至朝鮮半島,渤海受8~9級北到西北風(fēng)控制,各沿岸站依次開始出現(xiàn)減水。圖7為該過程5日08時(shí)、5日20時(shí)東亞地面天氣圖,圖8為葫蘆島、黃驊、龍口、煙臺(tái)風(fēng)暴潮、潮位模擬與實(shí)測對比。

        從圖6和圖8可以看出,各站無論從過程風(fēng)暴潮、潮位模擬都與實(shí)測相當(dāng)吻合,說明該模型能客觀刻畫渤海風(fēng)暴潮。

        圖5 2005年12月21日14時(shí)(a)、22日08時(shí)(b)東亞地面天氣圖(單位:hPa)

        圖6 葫蘆島(a)、塘沽(b)風(fēng)暴潮,葫蘆島(c)、塘沽(d)潮位模擬與實(shí)測對比

        圖7 2007年3月5日08時(shí)(a)、5日20時(shí)(b)東亞地面天氣圖(單位:hPa)

        圖8 2007年3月葫蘆島(a)、黃驊(b)、龍口(c)、煙臺(tái)(d)風(fēng)暴潮、潮位模擬結(jié)果與實(shí)測結(jié)果對比

        圖9為上述兩次過程模擬航路各代表點(diǎn)的風(fēng)暴潮、潮位過程曲線,綜合圖6、圖8和圖9可以看出:渤海風(fēng)暴減水具有范圍廣、時(shí)間長、減水大的特征,整個(gè)渤海的海水受離岸風(fēng)作用通過渤海海峽大量地流出,因此無論在渤海海灣內(nèi)或是渤海中央?yún)^(qū)的減水量值差別不大,在海灣內(nèi)的減水略大于海盆中的減水;由于這種天氣系統(tǒng)持續(xù)時(shí)間基本都超過1 d,因此渤海的較大減水過程持續(xù)時(shí)間往往超過1 d;另外,這種天氣系統(tǒng)給渤海帶來較強(qiáng)西北風(fēng),因此渤海的最大減水量值可以達(dá)到150~250 cm,一旦減水低值與天文低潮疊加將出現(xiàn)極低潮位,給大型油輪及巨型貨輪的正常航行和錨泊造成嚴(yán)重影響。

        圖9 兩次過程航路各代表點(diǎn)模擬風(fēng)暴潮、潮位過程曲線

        5 結(jié)論

        (1)利用環(huán)渤海9個(gè)沿岸站近10年潮位資料統(tǒng)計(jì)分析渤海海域的風(fēng)暴減水特征,結(jié)果表明:受天氣系統(tǒng)和地形影響,渤海年均出現(xiàn)50 cm和100 cm風(fēng)暴減水分別超過30 d和6 d,每年的9月至翌年4月份風(fēng)暴減水最為頻繁;渤海風(fēng)暴減水具有范圍廣、時(shí)間長、減水大的特征,海灣內(nèi)與渤海中央?yún)^(qū)減水量值差別不大,海灣內(nèi)的減水略大于渤海中央?yún)^(qū)的減水。

        (2)建立了一套精細(xì)化天文潮-風(fēng)暴潮耦合模型對各站天文潮進(jìn)行模擬驗(yàn)證,各站平均均方根誤差為18.5 cm,由此計(jì)算得到航路代表點(diǎn)的潮汐特征值表明航路從東向西潮汐類型由正規(guī)半日潮向不正規(guī)半日潮過渡,平均潮差緩慢遞減。

        (3)利用該耦合模型后報(bào)模擬了近10 a較大過程引發(fā)的風(fēng)暴減水和極低潮位,各站模擬與實(shí)測吻合較好,說明耦合模型可為大型船只在渤海深水航路航行提供有益的潮位預(yù)報(bào)。

        [1]吳少華,王喜年,戴明瑞,等.渤海風(fēng)暴潮概況及溫帶風(fēng)暴潮數(shù)值模擬[J].海洋學(xué)報(bào),2002,24(3):28-34.

        [2]王喜年.海洋災(zāi)害及預(yù)報(bào)[M].北京:海洋出版社,1991:43-88.

        [3]吳少華,王喜年,于福江,等.連云港溫帶風(fēng)暴潮及可能最大溫帶風(fēng)暴潮的計(jì)算[J].海洋學(xué)報(bào),2002,24(5):8-18.

        [4]董勝,于亞群,余海靜.海岸帶風(fēng)暴潮減水的統(tǒng)計(jì)分析[J].海洋科學(xué),2004,13(4):70-74.

        [5]胡筱敏,熊學(xué)軍,云升軍,等.2007年3月遼東灣西部減水特征分析[J].海洋科學(xué)進(jìn)展,2008,26(3):277-284.

        [6]Blain C,Westerink J,Luettich R,et al.Grid convergence studies for the prediction of hurricane storm surge[J].International Journal For Numerical Methods In Fluids,1998,26:369-401.

        [7]Westerink J,Luettich R,et al.A basin-to channel scale unstructured grid hurricane stormsurge model applied to Southern Louisiana[J].Monthly Weather Review,2008,136:833-863.

        [8]趙洪,楊學(xué)聯(lián),邢建勇,等.WRF與MM5對2007年3月初強(qiáng)冷空氣數(shù)值預(yù)報(bào)結(jié)果的對比分析[J].海洋預(yù)報(bào),2007,24(2):1-8.

        Numerical simulation study on deepwater channel influenced by negative storm surge and its features in Bohai Sea

        Fu Cifu1,F(xiàn)u Xiang1,Wu Shaohua1,Yu Fujiang12,Dong Jianxi1
        (1.National Marine Environmental Forecasting Center,Beijing 100081,China;2.Key Laboratory of Research on Marine Hazards Forecasting of State Oceanic Administration,Beijing 100081,China)

        The negative storm surge features based on the last 10 a tidal data of 9 coastal tide gauges in Bohai Sea has analyzed in this paper,the results show that:average annual negative storm surge of 50 and 100 cm is more than 30 d and 6d respectively,and September to April next year is the most frequently.A tide-storm surge coupled model with high resolution has been simulated tide,the average RMSis 18.5 cm for the coastal tide gauge,and the tidal characteristic quantities of each point at deepwater channel has been calculated for tide prediction.The large negative storm surge for last 10 a has been simulated base on the coupled model and the calculated results agree well with the observations,indicates that the coupled model can provide a useful reference on the tide forecast for the deepwater channel in Bohai Sea.

        negative storm surge;deepwater channel;high resolution;surge-tidecoupled;numerical simulation

        P731.34

        A

        0253-4193(2013)03-0030-09

        2013-03-09;

        2013-05-07。

        海洋公益性行業(yè)科研專項(xiàng)經(jīng)費(fèi)項(xiàng)目(2012418017)。

        傅賜福(1983—),男,福建省泉州市人,工程師,從事風(fēng)暴潮預(yù)報(bào)預(yù)警及研究工作。E-mail:fucf@nmefc.gov.cn

        傅賜福,付翔,吳少華,等.渤海風(fēng)暴減水特征及其對深水航路影響的數(shù)值模擬[J].海洋學(xué)報(bào),2014,36(3):30—38,

        10.3969/j.issn.0253-4193.2014.03.004.

        Fu Cifu,F(xiàn)u Xiang,Wu Shaohua,et al.Numerical simulation study on deepwater channel influenced by negative storm surge and its features in Bohai Sea[J].Acta Oceanologica Sinica(in Chinese),2014,36(3):30—38,doi:10.3969/j.issn.0253-4193.2014.03.004.

        猜你喜歡
        風(fēng)暴潮潮位航路
        基于距離倒數(shù)加權(quán)的多站潮位改正方法可行性分析
        2012年“蘇拉”和“達(dá)維”雙臺(tái)風(fēng)影響的近海風(fēng)暴潮過程
        唐山市警戒潮位標(biāo)志物維護(hù)研究
        基于實(shí)時(shí)航路的PFD和ND的仿真研究
        防范未來風(fēng)暴潮災(zāi)害的綠色海堤藍(lán)圖
        科學(xué)(2020年4期)2020-11-26 08:27:00
        基于多變量LSTM神經(jīng)網(wǎng)絡(luò)模型的風(fēng)暴潮臨近預(yù)報(bào)
        多潮位站海道地形測量潮位控制方法研究
        基于改進(jìn)的OLS-RBF模型的感潮河段潮位預(yù)測研究
        應(yīng)召反潛時(shí)無人機(jī)監(jiān)聽航路的規(guī)劃
        托勒密世界地圖與新航路的開辟
        国产免费爽爽视频在线观看| 国产一区二区三区观看视频| 精品人妻夜夜爽一区二区| 青青草成人免费在线视频| 久久久中文久久久无码| 无码av天天av天天爽| 精品日产卡一卡二卡国色天香| 免费一级特黄欧美大片久久网 | 射死你天天日| 亚洲V无码一区二区三区四区观看 久久精品国产亚洲综合色 | 日韩亚洲精品中文字幕在线观看 | 亚洲男人第一无码av网站| 亚洲精品亚洲人成在线下载| 白色白色在线视频播放平台| 国产美女高潮流白浆视频| 穿着白丝啪啪的av网站| 国产精品久久久久久妇女| 国产人妻人伦精品1国产盗摄| 亚洲中文字幕无码久久2018| 亚洲AV无码日韩综合欧亚| 亚洲综合久久精品少妇av| 亚洲乱码中文字幕久久孕妇黑人| 亚洲人成网7777777国产 | 国产品精品久久久久中文| 侵犯了美丽丰满人妻中文字幕| 国产精品videossex国产高清| 国产乱子伦精品免费无码专区| 级毛片无码av| 青青草成人在线播放视频| 亚洲国产精品ⅴa在线观看| āV第三区亚洲狠狠婷婷综合久久| 日本一区二区三区在线视频观看| 色中文字幕在线观看视频| 国产成人无码区免费内射一片色欲 | 人妻无码第一区二区三区| 最近中文字幕mv在线资源| 人妻无码∧V一区二区| 天天色天天操天天日天天射| 国产日产综合| 制服丝袜视频国产一区| 亚洲精品熟女av影院|