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

基于MIKE21模型的小清河入海口感潮河段水動力模擬及洪水風險分析

2025-06-28 00:00:00胡文雁周肆訪徐征和潘維艷叢鑫李昱岐
濟南大學學報(自然科學版) 2025年3期
關鍵詞:模型

Hydrodynamic Simulation and Flood Risk Analysis of Tide-sensitive Section of Xiaoqing River Estuary Basedon MIKE 21 Model

HU Wenyan 1 , ZHOU Sifang 2 , XU Zhenghe', PAN Weiyan 1 ,CONG Xin ?1 ,LI Yuqi 1 (1.School of Water Conservancy and Environment,Universityof Jinan,Jinan 25OO22,Shandong,China; 2.Shandong Hydrology and Water Environment Technology Co.,Ltd.,Jinan 25Oo14,Shandong,China)

Abstract:Toresearchthehydraulicchangecharacteristicsofriver estuaryandtheriskcaused bydiferent extreme hydrological events,atwo-dimensional hydrodynamicnumerical modelof thetide-sensitive sectionof Xiaoqing River estuary was establishedbyusing MIKE21model.Combinedwiththeflood protectionstandardof theriversection,fourcombinationsof floodand tide encounter weresetupto simulateand analyzethe characteristicsof the water flowchanges of the Xiaoqing River estuaryunder diferentfloodand tidecombinationscenarios,andanalyze theriskof the estuaryunder diferentscenariosaccording tothesimulationresults.Theresultsshowthatwhen therunoff conditions arethesame,the lowerthetidallevelinthelowerreachesoftide-sensitivesection is,thegreater thedecrease in water levelalong he courseis,andthe fastertheflowvelocityis.When thetidallevelsinthe lowerreachesarethesame,thegreater he runoff effect in theupperreachesof tide-sensitivesectionis,the greater the variationrangeofthewaterlevel along the course is,andthefastertheflowvelocityis.Thewatersurface transversedropofthefourcombinationsoffloodtides is gradually increasing inthebendsandtheriskoffloodtidecombinations1and2is higher,andallofthemreach the super-standard flood after the model simulation.

Keywords: tide-sensitive section; hydrodynamics; flood-tide combination;numerical simulation; Xiaoqing River

入海口是一個受陸地和海洋相互作用的特殊地帶,具有生態平衡、水質凈化、水文調節、航運港口等功能。近年來,隨著臺風“泰利”“木蘭”等極端水文事件頻發,造成大量人員傷亡、房屋倒塌,不僅阻礙沿海地區經濟發展,而且還對當地自然生態環境造成難以恢復的損失。此外,洪水導致水庫溢洪,增加水資源供需矛盾,造成水資源管理困難。相關數據顯示,風暴潮、海浪和赤潮災害為我國主要海洋災害,其中風暴潮過程造成的直接經濟損失最為嚴重。2022年我國沿海共發生風暴潮13次,其中5次造成災害,直接經濟損失23.7億元。廣西壯族自治區由臺風影響引起的風暴潮次數最多,為3次。“221003”號風暴潮對山東省東營、濰坊兩地海岸防護、鹽田、水產養殖的經濟損失占風暴潮總經濟損失的 47.5%[1] ,因此加強對沿海地區入海口區域水文變化研究,增強抵御風暴能力極為重要[2-4] O

現階段,物理模型和數值模型是進行水文模擬的主要方法,相較于物理模型,數值模型具有成本低、精度高等優點,常被用于解決較為復雜的水力學問題,常用的模型如Fluent、ECOM、Delft3D、MIKE等,其中MIKE模型覆蓋面廣,操作靈便,被廣泛應用于各項水流運動模擬中[5]。曹宇航等[6]基于MIKE11水動力模型對區域河網水位及洪水過程進行模擬預測。金玲等[7]以遼寧省大連市復州河下游為研究對象,構建一維水動力模型,探究設計洪水演進過程及風險指標。郭鳳清等[8]運用MIKE21模型中FLOWMODELFM模塊,對/江蓄滯洪區水量、流速進行動態變化模擬。王天澤等搭建了北京市某科學城二維水動力模型,探究洪水進程及淹沒面積在現狀與規劃條件下的變化,為城市規劃提供大尺度、高精度模擬研究。

小清河發源于山東省市南部山區,自西向東匯集20個縣(市、區)的來水,于山東省壽光市羊口鎮東注入萊州灣。小清河是魯中地區重要的排水河道,也是山東省唯一的省級直管河道,在全省國民經濟中具有至關重要的地位[10]。近些年,隨著全球氣候變暖,小清河流域水文循環加劇,尤其是東部小清河入海口地區遭遇極端水文事件的頻率相較于其他地區大。相關調研顯示,2018年臺風“溫比亞”期間,小清河下游發生了嚴重的洪澇災害,降雨僅1d就幾乎達到了大暴雨量級。2019年臺風“利奇馬”期間,小清河全流域再次發生洪水,尤其在小清河入海口干流岔河、石村水文站均發生了有記錄以來最大洪峰體積流量 642m3/s 。小清河流域河網復雜,其入海口感潮河段常遭受上游洪水與下游潮位雙重影響。國內許多學者對小清河感潮河段都進行了研究。彭慧等[3]通過歷年觀測資料,利用組合風險率模型確定了小清河感潮河段最優洪潮遭遇組合。慕金波等[1]利用一維水質模型預測了工程建設對小清河感潮河段的水質影響。金玉休等[12]基于二維溫排水數學模型探究小清河下游感潮河段溫排水擴散與地形地貌的相互關系。李嘉[13]基于Delft3D模型系統模擬預測小清河感潮河口抗生素的傳輸與分布。上述研究主要針對小清河感潮河段的洪潮組合聯合概率與水質演化走勢,運用MIKE21模型對小清河入海口感潮河段區域開展水動力與風險的研究鮮有涉及,因此,本文中采用MIKE21模型建立小清河入海口感潮河段二維水動力數值模型,模擬不同洪潮組合情景下感潮河段的水流變化特征并開展風險分析,以期為小清河入海口水文風險應對提供科學依據。

研究區概況

研究區位于小清河入海口的石村水文站至羊角溝潮位站(地理位置為東經 118°2545′′-118°5208′′. 北緯 ,見圖1。該區域屬“一日兩潮”的不規則半日潮河段,受到上游洪水和下游潮汐頂托的共同影響,又稱感潮河段,河段長約50km 。小清河入海口屬華北暖溫帶半濕潤季風型大陸性氣候,年內四季分明,溫差變化較大,年平均降水量約為 600mm ,降雨主要集中在6—9月份,降雨量約占全年降水量的 72.9% 。主要支流包括淄河、新塌河以及老塌河。研究區地勢平緩,處于平原地區,河床橫比降為 1/6000~1/8000 0

2 數據來源及模型構建

2.1 數據來源

地形數據來自地理空間數據云,氣象數據來自中國氣象數據網,水文資料主要源于《淮河流域水文資料》。由于氣旋、切變線、鋒面等溫帶天氣系統,研究區較大暴雨洪水多發生在7、8月份,因此,本文中選用石村水文站及羊角溝潮位站2022 年7—8月份逐時平均降雨、蒸發、風速風向、流量及水位數據開展模型模擬。

2.2 模擬范圍

模擬范圍為小清河下游的石村水文站與羊角溝潮位站之間。王道閘是位于石村水文站與羊角溝潮位站之間的擋潮閘,其上游范圍的潮汐作用較弱。

圖1小清河入海口感潮河段概況
(資料來源:國家測繪地理信息局標準地圖服務網站,經ArcGIS10.8軟件處理。)

為了重點考慮不同洪潮對人海口風險的影響,模型重點輸出范圍選為王道閘與羊角溝潮位站之間,兩站點之間干流長度約為 29.5km ,流域面積為

2.3 網格劃分及參數設置

利用ArcGIS10.8軟件提取岸線及河道點位的經、緯度作為岸線文件land.xyz,參考小清河防洪綜合治理工程的相關報告及2019年5月實測河道大斷面數據,生成用于河道地形文件的水深文件water.xyz。用SMS11.2軟件將岸線文件進行平滑處理及非結構化三角形網格劃分,共生成19655個網格和11487個節點。轉換為Mesh文件后,導入水深文件,插值后生成河道地形文件。研究區網格劃分及地形插值結果如圖2所示。

圖2小清河入海口感潮河段網格劃分及地形插值示意圖

為了確保模型模擬結果的準確性,建立水動力模型時參考小清河防洪綜合治理工程的相關報告中給出的河道設計糙率值,建立空間變化的曼寧數文件。根據實測氣象數據,建立隨時間、空間變化的降水文件和風場文件。模型采用時間步長及時間步數分別為1800s和156。初始條件為dfsu文件,邊界條件為石村水文站逐小時流量時間序列。考慮到不同支流入流對小清河入海口河段水流影響的不同,將各支流的實測流量資料以源匯項形式加入模型中。其他參數均采用模型默認值。

根據收集到的水文資料,選擇2022年7月30日至8月6日對模型進行率定和驗證。采用Pearson-IⅢ型適線法分析率定與驗證期的場次洪水頻率,洪水頻率曲線見圖3(a)。根據經驗頻率點的分布,并考慮參數在地區上的協調,選用 Cv=0.15,Cs/Cv=4.33 ( Cv 為變差系數, Cs 為偏態系數)進行試線,經計算,王道閘處率定與驗證時期洪水重現期為100、50a的設計洪水的體積流量分別為 255.96,265.53m3/s 。

率定與驗證時期的實測潮位曲線見圖3(b)。由圖可見,研究區實測最高潮位為 2.56m ,最低潮位為 0.62m ,平均潮位為 1.45m ,平均漲潮歷時5.73h ,平均落潮歷時 5.93h 。

標題

為了評價模型模擬結果,采用相關系數、納什系數以及相對誤差分析模擬值與實測值之間的相關性和吻合程度。

2.4水動力率定及驗證結果

圖4所示為二維水動力模型在王道閘處水位率定和驗證的結果。由圖4(a)可以看出,率定期模擬水位與實測水位相差不大且變化趨勢吻合較好,相關系數為0.86,納什系數為0.66,相對誤差為11.87% ,說明模擬效果較好。在率定好的水動力模型基礎上,用2022年8月2—6日水位數據驗證模型,結果如圖4(b)所示。由圖可以看出,驗證期王道閘水位的模擬值和實測值吻合較好,相關系數為0.89,納什系數為0.70,相對誤差為 12.24% ,均在可接受范圍內。總體來說,構建的模型模擬效果較好,合理可行,可用于不同情景下小清河入海口水流變化特征的模擬預測

3 情景設置及結果分析

3.1 情景設置

本文中擬直接引用文獻報道的采用定性與定量相結合方法,獲得小清河感潮河段洪潮遭遇的組合概率結果。彭慧等通過Copula函數,首先定性分析小清河流域洪水與潮位遭遇的可行方案集,然后定量計算組合方案集的風險概率,最后得出最優洪潮組合為重現期為50a的洪水遭遇重現期為20a的潮位和重現期為50a的潮位遭遇重現期為20a的洪水,與陳翠英等[14]研究得出的小清河流域潮位遭遇組合基本一致。本文中在設置情景模擬方案時,上游采用重現期為50、20a的設計洪水過程線,下游使用羊角溝潮位站重現期為50、20a的設計潮位過程線。考慮到小清河入海口處各支流非城區段防洪標準為重現期為 20a ,因此各支流均采用重現期為20a設計洪水過程線,且洪峰與上游洪水過程洪峰同時發生。小清河入海口段洪水和潮位遭遇的4種組合方案見表1。選用小清河“180821”號洪水發生期的典型洪水過程線和典型潮位過程線,同倍比放大后得到設計洪峰流量和設計汛期最高潮位。本文中將從不同洪潮組合情景下的縱向流態分布及河道彎道處橫向流態分布2個方面進行分析,進而分析判斷感潮河段是否存在漫堤現象,探討入海口的洪水風險。

標題
表1小清河入海口感潮河段洪潮組合方案

3.2 不同洪潮組合情景下的縱向流態分布

采用MIKE21模型,按照進口斷面洪峰時期與出口斷面下游潮位最高時期去描述縱向流態在不同洪潮組合下的沿程變化。小清河入海口感潮河段洪峰與下游潮位最高時期沿程水位和流速如圖5所示。從洪峰時期水位圖可知,水位最高的為洪潮組合1,水位最低的為洪潮組合4。對比洪潮組合1、2,水位分別沿程下降 64.66% 、 69.62% ,二者水位差呈增加趨勢,水位差變化幅度從開始的 1.39% 增至 15.23% ,說明在徑流相同時,下游潮位越低,水位沿程下降幅度越大。對比洪潮組合1、3,在距王道閘 29.03km 處,水位均為 2.053m 時,組合1、3的水位下降幅度分別為 64.26% 、 59.38% ,說明在下游潮位相同時,上游徑流作用越大,水位沿程下降幅度也就越大。上述模擬結果表明,水位沿程變化幅度隨著徑流作用越大或下游潮位越低時越大,與圖5(b)中下游潮位最高時期的水位變化結果一致。

由于在下游潮位達到最高時會受到下游潮水位頂托影響,入海口泄洪受限,因此,相比洪峰時期,水位在潮位最高時期的下降幅度明顯較小。對于洪峰時期流速變化,如圖5(c)所示,流速最快的是洪潮組合2,平均流速為 1.82m/s ,流速最慢的是洪潮組合3,平均流速為 1.56m/s 。通過比較洪潮組合1、2可知,在徑流相同時,下游潮位越大,導致泄洪能力受限,造成流速越慢;通過比較洪潮組合1、3可知,當下游潮位相同時,徑流作用越大,流速越快。上述模擬結果表明,在小清河入海口感潮河段,徑流作用越大或下游潮位越低,流速越大,與圖5(d)中下游潮位最高時期的流速變化結果一致。

3.3 不同洪潮組合情景下彎道處橫向流態分布

河道彎道的水流流態通常比直道的更復雜[15]針對研究區水流流向自東北方向至向東方向處的彎道,不同洪潮組合情景下彎道處水位分布與流矢圖如圖6所示。由圖可知,在4種洪潮組合情景下,彎道處凹岸的水位均高于凸岸的,形成了水面橫比降。流矢圖中箭頭疏密程度表明流速的快慢,箭頭越密集表明流速越快;反之,流速越慢。由圖6可以看出,在進入彎道前水體中間水流的流速快于兩岸處水流的,在進入彎道后水流流速整體有所變慢。

圖5小清河入海口感潮河段洪峰與下游潮位最高時期沿程水面線、流速
(a)洪潮組合1(b)洪潮組合2(c)洪潮組合3

為了探究水面橫比降在彎道的沿程變化,在彎道處設置了3個斷面,分別是水流進入彎道處斷面(斷面1)、彎道頂部斷面(斷面3)以及前2個斷面中間處斷面(斷面2)。不同洪潮組合中3個斷面水面橫比降的模擬結果如圖7所示。由圖可以看出,從彎道入口到彎頂斷面附近,水面橫比降均呈現增大的趨勢。從斷面1到斷面3,洪潮組合2、4的水面橫比降增長率分別為 69.82% 、 70.04% ,增長率相差不大。對比洪潮組合1、3,在同為重現期50a設計潮位時,洪潮組合1在3個斷面的斷面橫比降分別為 0.111%00.0.147%0.0.188% ,洪潮組合3各斷面水面橫比降分別為 0. 084%0. 0. 111%0. 0. 142% 表明洪潮組合1在3個斷面處的水面橫比降均大于洪潮組合3的。總體來看,在下游潮位相同時,上游流量越大,彎道處水面橫比降越大,但水面橫比降增長率相差不大。

圖7小清河入海口感潮河段不同洪潮組合情景下各斷面水面橫比降

4 不同洪潮組合情景下入海口風險分析

根據《水利水電工程等級劃分及洪水標準》(SL252—2017)[16]和小清河干流治理工程歷次規劃確定防洪標準,小清河入海口感潮河段將洪水等級分為一般洪水、現狀標準內洪水以及超標準洪水,3個洪水等級對應的警戒狀態遞增,劃分界限為除澇水位和防洪水位,小清河防洪綜合治理工程的相關報告中感潮河段10個斷面樁號的防洪、除澇水位如表2所示,10個斷面樁點的位置見圖1。對比感潮河段洪峰時期與最高潮位時期的水位模擬可知,洪峰時期水位整體最高,風險最大,因此將圖5(a)中洪峰時期在10個斷面處水位模擬結果(表3)與各斷面防洪、除澇水位(表2)進行對比,得到各斷面洪水等級,如表4所示。

在洪潮組合1情景下,王道閘至引黃濟青倒虹河段3個斷面的模擬水位分別為 5.75、5.64、5.45m 均高于對應位置處的防洪水位 5.36,5.18,5.01m 因此該河段屬于超標準洪水。老塌河口至羊角溝潮位站河段7個斷面的模擬水位分別為4.67、4.47、4.33、4.33、3.63、2.78、2.04m,均在相應位置的除澇水位與防洪水位之間,因此該河段屬于現狀標準內洪水。綜上,洪潮組合1情景下小清河入海口感潮河段的洪水風險最大,對防洪最不利。洪潮組合2在王道閘至廣饒壽光分界河段與洪潮組合1的結論一致,但在羊角溝潮位站處的模擬水位為1.72m ,低于除澇水位 2.03m ,該斷面為一般洪水,因此,與洪潮組合1相比,洪潮組合2的洪水風險次之。洪潮組合3各斷面均為現狀標準內洪水,因此河段現有防洪標準可抵御洪潮組合3帶來的洪水風險。洪潮組合4洪水等級既有現狀標準內洪水,也有一般洪水,是最可能發生的情景概率,但洪水風險較小。

表2小清河入海口感潮河段各斷面除澇、防洪水位
表3小清河入海口感潮河段洪峰時期各斷面模擬水位
表4小清河入海口感潮河段不同洪潮組合情景下各斷面洪水等級
注: ① “√”表示斷面所屬洪水等級。

綜合各河段的模擬結果發現,王道閘至引黃濟青倒虹河段洪水風險相對較高,其中在洪潮組合1、2情景下,洪水等級達到超標準洪水,因此在小清河入海口發生洪潮遭遇時,要重點關注王道閘至引黃濟青倒虹河段干流的洪水走勢,重點監測雨水情變化,加強對該河段干流以及淄河口干流的堤防巡查,詳細了解當地區域洪水特點,因勢利導,綜合運用工程和非工程措施,及時采取應急防護和搶險舉措,把水文災害風險降到最低。

由不同洪潮組合情景下彎道處橫向流態特征分析可知,當水流進人河道彎道后,凸岸水位低于凹岸水位,并且隨著水流前進,凸岸與凹岸之間的水位差呈增大趨勢。根據孫東坡等[17]對河流彎道橫向流動和輸沙研究結果,橫斷面內形成的環流與縱向水流相結合形成螺旋狀流動,進而在彎道處呈現“凹沖凸淤\"現象。除了水位變化(水面橫比降)會影響斷面水面線之外,彎道處床面橫比降也會對斷面水面線產生一定影響[18],此外,進入彎道后水流產生的螺旋流將對河堤產生破壞,影響船舶正常行駛。因此,在應對水文災害的同時,要在出現潰堤險情前及時預警,提前采取應對措施。重點關注彎道處河槽水位等水文要素的變化,并且在整治河道彎道時,充分考慮水面超高特點和床面橫比降,保證船舶正常航行。

5結論

本文中采用MIKE21模型建立小清河入海口感潮河段二維水動力模型,分析不同洪潮組合情景下的水流變化特征與洪水風險,得到以下主要結論:

1)在4種洪潮組合情景下,當徑流條件相同且下游潮位越低時,或當下游潮位相同且上游徑流作用越大時,水位沿程下降幅度越大,流速越快,流速波動幅度越小。

2)4種洪潮組合在彎道處均會形成水面橫比降,上游流量越大,各斷面在彎道處水面橫比降越大,因此,小清河入海口感潮河段要重點監測彎道凹岸水位,加強對彎道堤岸的巡查,改善彎道水流流態,把對船舶通航穩定性影響的降到最小。

3)洪潮組合1、2的洪水為超標準洪水,要重點關注這2種組合情景下存在的危險性。對王道閘至各支流洪水趨勢開展重點監測,加強堤防巡查,采取有效措施應對水文災害。

參考文獻:

[1]佚名.2022年中國海洋災害公報(摘登)[N].中國自然資源報,2023(5).

[2] 宋曉猛,張建云,占車生,等.氣候變化和人類活動對水文循環影響研究進展[J].水利學報,2013,44(7):779.

[3] 彭慧,劉璐,姜鈞耀,等.感潮河段上游洪水與河口潮位遭遇研究[J].人民黃河,2021,43(8):44.

[4]朱日清,姜尚堃,楊恒偉.山東省水旱災害防御技術支撐工作經驗做法[J].山東水利,2023(4):12.

[5]穆聰,李家科,鄧朝顯,等.MIKE 模型在城市及流域水文-環境模擬中的應用進展[J].水資源與水工程學報,2019,30(2): 71.

[6]曹宇航,袁文秀,李衛東,等.基于MIKE11模型的平原感潮河網城市防洪規劃研究[J].中國農村水利水電,2017(4):17.

[7]金玲,許士國,于德全.中小河流洪水風險分析中的洪水演進計算研究[J].水電能源科學,2014,32(10):48.

[8]郭鳳清,屈寒飛,曾輝,等.基于MIKE21的泄江蓄滯洪區洪水危險性快速預測[J].自然災害學報,2013,22(3):144.

[9]王天澤,王遠航,馬帥,等.基于MIKEFLOOD耦合模型的洪水淹沒風險分析:以北京市某科學城為例[J].水利水電技術(中英文),2022,53(7):1.

[10]劉曉麗.市小清河生態系統健康評價、修復與服務價值研究[D].:山東大學,2008:2.

[11]慕金波,鄭海春,董景廉.一維有限段模型在小清河河口水質預測中的應用[J].環境科學進展,1993(3):45.

[12]金玉休,朱龍海,吳建政,等.小清河感潮河段與河口區溫排水模擬及對地貌的響應[J].海洋科學,2014,38(11):74.

[13]李嘉.感潮河口抗生素的反應遷移機制及傳輸模擬[D].煙臺:中國科學院煙臺海岸帶研究所,2018:109.

[14]陳翠英,向征,王海軍,等.山東省小清河流域暴雨-潮位遭遇組合規律研究[J].水電能源科學,2022,40(2):10.

[15]高偉,楊中華.彎道水流特性和數值模擬方法研究進展[J].水電能源科學,2009,27(1):112.

[16]水利水電工程等級劃分及洪水標準:SL252—2017[S].北京:中國水利水電出版社,2017:3.

[17]孫東坡,朱岐武,張耀先,等.彎道環流流速與泥沙橫向輸移研究[J].水科學進展,2006(1):61.

[18]余海逖,周蘇芬,閆旭峰,等.U形彎道河床橫比降對水流運動特性影響試驗研究[C]/中國水利學會.2018學術年會論文集:第三分冊.北京:中國水利學會,2018:434.

(責任編輯:劉飚)

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 91热爆在线| 国内精自线i品一区202| 国产成人夜色91| 麻豆国产精品视频| 九九九国产| 国产主播在线观看| 欧美狠狠干| 青青青亚洲精品国产| 午夜激情婷婷| 天堂网亚洲综合在线| 亚洲国产一区在线观看| 亚洲女人在线| 中文字幕亚洲无线码一区女同| 国产日本一线在线观看免费| 亚洲欧美日韩成人在线| 成人在线第一页| 2019年国产精品自拍不卡| 亚洲国产欧洲精品路线久久| 国产精品男人的天堂| 国产成人亚洲毛片| 国产精品理论片| 五月激情婷婷综合| 午夜免费视频网站| 日本福利视频网站| 国产天天色| 精品超清无码视频在线观看| 色天堂无毒不卡| 婷婷色一二三区波多野衣| 毛片视频网址| 亚洲婷婷丁香| 国产久操视频| 久久99精品久久久久纯品| 国产精品观看视频免费完整版| 99久久国产综合精品2023| 又爽又大又光又色的午夜视频| 欧美日韩精品一区二区在线线| 国产毛片基地| 久久美女精品国产精品亚洲| 日韩欧美高清视频| 中文字幕在线看| 免费在线播放毛片| 国产黄色免费看| 精品久久人人爽人人玩人人妻| 国产伦精品一区二区三区视频优播| 欧美中文字幕在线视频| 国产精品成人一区二区| 亚洲第一黄片大全| 在线无码九区| 国产亚洲一区二区三区在线| 亚洲综合色吧| 911亚洲精品| 伊人激情综合网| 亚洲成人在线网| 一本一本大道香蕉久在线播放| 波多野结衣爽到高潮漏水大喷| 欧美高清日韩| 日韩欧美中文亚洲高清在线| 91亚洲国产视频| 99久久免费精品特色大片| 色九九视频| 欧美一道本| 少妇人妻无码首页| 欧美精品黑人粗大| 免费亚洲成人| 狠狠干欧美| 久久99蜜桃精品久久久久小说| 久久成人免费| 亚洲AⅤ波多系列中文字幕| 亚洲—日韩aV在线| 99re热精品视频中文字幕不卡| 欧美无专区| 欧美色伊人| 亚洲国产欧美自拍| 麻豆国产精品一二三在线观看| 欧美精品二区| 国产激情在线视频| 国产精品福利一区二区久久| 亚洲毛片一级带毛片基地| 国产手机在线观看| 国产一级无码不卡视频| 成人免费网站在线观看| 91精品啪在线观看国产91|