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

水庫泄流洪水波在宜昌—大通河段傳播規律研究

2025-02-19 00:00:00仇紅亞劉亞新李妍清張冬冬熊豐
人民長江 2025年1期

摘要:

在當前水利工程和防洪管理中,水庫泄流洪水波傳播問題備受關注。以長江宜昌—大通干流河段為例,構建了一維水動力學模型,模擬了不同工況條件下洪水波的傳播特性,并對比分析了沿程洪水演進規律。研究結果表明:① 斷波洪水的形成可在不同程度上縮短下游河段的峰現時間,隨著斷波流量的增加,以運動波傳播的河段距離明顯縮短,洪水波形變指數增加;② 隨著洪峰流量持續時間的延長,同樣會縮短以運動波特性傳播的洪水波演進的距離,并增加洪水波形變指數,但同一站點或是里程上的峰現時間會有推后,洪峰歷時延長;③ 多波束斷波條件下,計算波速普遍較高,隨著沿程距離的增加波形變指數遞減,但多波束傳播條件下形變指數沿程衰減過程明顯放緩;④ 當有支流入匯時,沿程站點洪峰流量變化均表現出明顯的非一致性,高流量歷時明顯延長。研究成果對于揭示水庫泄流洪水波傳播規律具有重要意義,可為區域防洪管理提供參考。

關" 鍵" 詞:

水庫泄流; 洪水波; 傳播特性; 形變指數; 水動力學模型; 長江

中圖法分類號: P332.2

文獻標志碼: A

DOI:10.16232/j.cnki.1001-4179.2025.01.004

收稿日期:2024-07-24;接受日期:2024-09-11

基金項目:

國家重點研發計劃項目(2022YFC3002703);中國長江三峽集團有限公司科研項目(0704221);國家自然科學基金項目(52309004,U2340205)

作者簡介:

仇紅亞,男,工程師,博士,主要從事水文水資源研究。E-mail:qiu_hongya@ctg.com.cn

Editorial Office of Yangtze River. This is an open access article under the CC BY-NC-ND 4.0 license.

文章編號:1001-4179(2025) 01-0022-09

引用本文:

仇紅亞,劉亞新,李妍清,等.

水庫泄流洪水波在宜昌—大通河段傳播規律研究

[J].人民長江,2025,56(1):22-30.

0" 引 言

水庫泄流洪水波傳播是當前水利工程和防洪管理中備受關注的問題[1-3]。研究人員針對該領域展開了一系列研究工作,其中數值模擬成為研究的主要方法之一[4-5]。近年來,數學模型在水庫泄流洪水波傳播領域的應用研究取得了多項進展。例如,基于復雜的水動力模型和地形數據,研究者能夠更精確地預測洪水波的擴散和淹沒范圍,從而提前進行預警和采取應對措施[6]。同時,部分模型還考慮了水庫泄流對水質的影響,如泥沙輸移和水質變化,這對水資源保護和環境管理尤為重要。

隨著計算能力和數據處理技術的提升,數學模型不斷發展,從簡化模型向復雜的耦合系統模型過渡,增強了對不同調度條件下水庫泄洪過程的理解和預測能力[7-8]。通過建立數學模型,研究人員可以在不同尺度和復雜情景下進行模擬和預測,為防洪管理決策提供科學依據[9-10]。此外,數學模型方法還能夠充分發揮計算機技術的優勢,實現洪水波傳播過程的快速、大規模計算,并且在一定程度上彌補了野外實驗難以覆蓋各種情況的局限性[11]。因此,數學模型方法在水庫泄流洪水波傳播規律研究中具有重要意義,為解決實際工程和管理問題提供了強有力的支持。

隨著長江流域梯級水資源的開發,特別是三峽水庫蓄水運用以來,長江水文情勢發生了顯著的變化[12-13]。水庫進行應急搶險或是緊急防洪調度這一過程經常會引起泄流洪水在短時間內發生較大變化,甚至形成急變洪水波,而急變洪水波大多伴隨著斷波的形成[14-15]。程海云等[16]分析了三峽水庫運行后荊江河段洪水傳播時間特性,并進一步揭示了斷波洪水與沙市站水位流量關系的響應規律。鑒于目前已有研究成果僅限于長江中下游荊江河段,泄流洪水波在螺山站以下河段傳播特性有待進一步研究,因此有必要深入研究洪水波傳播時間和傳播距離以及削減過程等動力特征。本文以長江宜昌至大通干流河段為例,通過構建一維水動力學模型,模擬不同工況下水庫泄流洪水波的傳播特性,并分析其演進規律,以為區域防洪管理提供參考。

1" 一維水動力學模型構建

1.1" 研究區域

本次研究以長江宜昌至大通干流河段為研究對象。該河段長約1 183 km,其中宜昌至湖口為中游河段,沿程有清江、洞庭湖、漢江、鄱陽湖及其他支流入匯;湖口至大通河段為長江下游的一部分,長約228 km,該河段基本不受潮汐影響。宜昌至大通干流河段所處地理位置及重要站點分布如圖1所示。

1.2" 模型建立與參數率定

斷波在向下游傳播的過程中,會受到地形等多種因素的影響而發生變形,這導致水面比降在整個過程中變得緩慢,斷面呈現出漸變流的特性,因此可以使用完整的動力波模型來進行近似模擬。

本文采用河網水流數學模型進行建模以及分析計算。該模型基于圣維南原理,由連續方程(1)和動量方程(2)組成一維水動力學控制方程組,基本形式如下[17]:

B(η/t)+Q/x=0

(1)

Qt+gAηx+(Q2/A)x+gn2QQAR4/3=0

(2)

式中:B為斷面水面寬度,m;

η為斷面中心水位,m;

t為時間,s;

Q為斷面流量,m3/s;

x為距離,m;

g為重力加速度,m/s2;

A為過水斷面面積,m2;

n為河道綜合糙率;

R為斷面水力半徑,m。

采用2020年實測斷面資料對研究區域建模,在宜昌—大通干流河段設置927個斷面,斷面間距為0.24~7.89 km。采用宜昌站流量過程作為模型計算上邊界,下邊界采用大通站水位-流量關系曲線,區間入匯的包括清江、洞庭湖、漢江和鄱陽湖等水系,模型中采用支流入匯流量過程。結合地形數據,選取2016~2018年實測數據進行參數率定,2019~2022年實測數據進行驗證。圖2~3為各站流量與水位率定與驗證結果,率定驗證完成后各個河段糙率見表1。由圖2~3及表1可知,各站點流量、水位過程的驗證結果良好,模型參數合理,可用于宜昌—大通河段洪水特性分析計算[18]。

2" 宜昌—大通河段洪水波傳播特性分析

2.1" 波速計算

根據圣維南動力方程中各項的對比關系和數量級

大小,可以忽略某些次要項,將洪水波分為5類[19]:① 若忽略慣性項和附加比降,則稱為運動波;② 若僅忽略慣性項,則稱為擴散波;③ 若忽略摩阻項和河底比降(重力)項,則稱為慣性波;④ 若忽略局地慣性項,則稱為對流運動波;⑤ 若所有項均不能忽略,則稱為動力波。在分析河道洪水波傳播規律時,通常把運動波和擴散波統稱為運動波,將慣性波和動力波合并稱為動力波。

運動波波速可用式(3)計算,僅向下游傳播:

Ck=φ·V(3)

式中:Ck是運動波速;

φ是波速系數,φ=1+AVdVdA,顯然,φ是大于1的常數,隨著斷面面積和流速變化而變化,通常取值范圍為1≤φ≤2;V是研究河段斷面平均流速。

動力波的波速計算如式(4):

Cd=V+gh

(4)

式中:Cd是動力波波速;

g是重力加速度;

h是研究斷面平均水深。動力波可沿上下游兩個特征方向傳播,同流量條件下動力波波速一般大于運動波波速。大多數河道水流運動屬于運動波范疇,但在流量發生急劇變化時波型也會發生相應的轉換。

2.2" 洪水波形變指數

洪水波傳播過程中,波形不斷發生變化,對于洪水過境的每一個斷面都會出現一個漲水過程和退水過程,漲水過程的流量變率坡度如下:

dQdt=Qp-Qitp-ti

(5)

據此,引入洪水波形變指數[20]:

FI=(Qp-Qi)/Qi(tp-ti)/Δt

(6)

式中:Qp是洪峰流量;

Qi是洪水起漲流量;

tp是洪峰出現時的時間;

ti是洪水起漲時的時間;

Δt是時間差,此處為1 h。FI越接近0,洪水過程越平穩,數值越大洪水波動越劇烈。

2.3" 研究方案和計算工況

(1) 不同斷波流量計算方案。

根據三峽水庫和葛洲壩樞紐出庫流量的系列統計數據,計算模型上邊界宜昌站恒定流流量取15 000 m3/s,假定在5 h內流量分別驟增至 2 0000,25 000,30 000,35 000 m3/s,即出現斷波流量分別為5 000,10 000,15 000,20 000 m3/s,增加后的流量持續5 d后,在5 h內再驟減至15 000 m3/s,圖4(a)為4組斷波流量出現過程,其中支流清江、洞庭湖、漢江和鄱陽湖來水分別按照無入流和近10 a汛期平均流量兩個方案計算。

(2) 不同洪峰持續時間計算方案。

假設宜昌站以15 000 m3/s的恒定流起漲流量開始計算,假定在5 h內流量驟增至 25 000 m3/s,增加后流量再分別持續1,2,3,4,5 d后,再在5 h內驟減至15 000 m3/s,共計5 組不同歷時的斷波洪水過程(圖4(b)),其中支流清江、洞庭湖、漢江和鄱陽湖來水分別按照無入流和近10 a汛期平均流量兩個方案計算。

(3) 多波束泄流影響計算方案。

假設宜昌站以15 000 m3/s的恒定流起漲流量開始計算,假定在5 h內流量驟增至 25 000 m3/s,流量持續2 d后,再施放一組10 000 m3/s的斷波洪水過程,兩組洪峰持續3 d后在5 h內驟減至25 000 m3/s,繼續持續2 d后再在5 h 內驟減至15 000 m3/s,考察這2組斷波洪水的傳播過程(圖4(c)),其中支流清江、洞庭湖、漢江和鄱陽湖來水分別按照無入流和近10 a汛期平均流量兩個方案計算。

2.4" 結果與分析

2.4.1" 不同斷波流量的影響

洪峰流量大小通常對下游防洪形勢有著重要影響。對比不同斷波流量條件下洪水波的傳播規律(圖5),發現同一流量條件下,沙市站以上河段洪峰流量基本無削減,僅持續時間存在不同程度減少。當斷波流量為5 000 m3/s時,枝城和沙市站洪峰持續時間分別從5 d降至73.85 h和23.80 h;隨著距離的增加洪峰坦化作用越明顯,螺山、漢口、九江和大通站洪峰分別削減了2.88%,7.95%,11.25%和12.88%。然而,隨著斷波流量的增加,洪峰沿程坦化過程加劇,當斷波流量為15 000 m3/s時,螺山、漢口、九江和大通站洪峰分別削減了3.55%,8.25%,14.96%和19.73%。

當斷波流量增加,下游各站點的峰現時間均出現不同程度提前,以漢口站為例,當斷波流量為5 000 m3/s時,峰現時間出現在164.10 h,當斷波流量分別為10 000,15 000,20 000 m3/s時,峰現時間分別為148.20,146.10,144.95 h。此外,當有支流入匯時,洪水波傳播過程會更加復雜,支流的入匯導致洪峰流量在觀測站點出現了不同程度的增加,其中洞庭湖和鄱陽湖的入匯導致螺山和大通站的洪峰流量不同程度增加,洪峰持續時間更長。從表2可見,隨著斷波流量的增加,同一站點洪水波形變指數急劇增加,表明洪水波動越劇烈;隨著沿程距離的增加洪水波形變指數遞減,表明波峰在不斷的削減,但高流量的斷波洪水在同一站點仍存在較大的形變指數。

根據宜昌—大通河段模型斷面間的河道里程和洪

峰在斷面間的過境時間差,計算得出沿程斷波傳播波

速,并與運動波和動力波的計算波速對比分析,結果如圖6所示。當無區間支流入匯,宜昌站出現斷波流量為5 000 m3/s 和10 000 m3/s時,計算波速在宜昌站以下750 km內均大于2u(u為運動波波速),但仍小于計算動力波速Cd,傳播特性符合過渡區洪水過程;750 km以下至大通河段,計算運動波速Ck基本處于運動波波速Ck范圍(1u≤Ck≤2u),符合運動波傳播特性;而當斷波流量達到15 000 m3/s以上時,宜昌站以下約150 km范圍內,計算波速在動力波速上下變化,表現出動力波傳播特性,150 km以下范圍至大通河段,計算波速大于2u,小于動力波速Cd,表現為過渡區洪水波特性。當有區間支流入匯,計算波速表現出類似變化,所不同的是當斷波流量為10 000 m3/s時,宜昌站以下河段計算波速均大于2u,已無符合運動波傳播河段。

2.4.2" 不同洪峰流量持續時間的影響

上游洪水持續時間是預測評估下游可能的洪澇災害的重要指標。不同的上游可能洪水持續時間在下游的洪水演進過程計算結果如圖7所示。在無區間入流條件下,上游洪峰持續時間越長,下游觀測站點的洪峰流量越大,且洪峰的削減過程越緩慢,例如宜昌站洪峰流量持續1 d時,枝城、沙市、螺山、漢口、九江和大通站分別出現峰值流量為24 600,22 100,16 400,15 600,15 300 m3/s和15 300 m3/s,峰現時間分別為29.55,33.75,68.85,103.30,143.00 h和192.20 h,而當洪水持續時間延長到5 d時,以上各站點的峰值流量為

2 500,25 100,24 900,23 300,21 500 m3/s和20 600 m3/s,峰現時間分別為125.35,127.10,132.25,148.20,175.05 h和214.85 h。可見,在相同的洪峰流量條件下,延長上游洪水的持續時間可導致下游洪峰的大幅提高和歷時的增加。

此外,通過對比有區間支流入匯情況可見,支流的來流不同程度增加了觀測站點洪峰流量,隨著上游洪水歷時的延長,支流匯入情況下洪水在各站點的洪水歷時大幅延長,洪峰流量增加了1.12%~127.78%。從表2可見,隨著洪峰持續時間的增加,同一站點洪水波形變指數急劇增加,表明洪水波動越劇烈;隨著沿程距離的增加洪水波形變指數遞減,表明波峰也在不斷的削減,但洪峰持續時間越長形變指數衰減越慢。

計算波速結果如圖8所示,當無區間支流入匯,計算波速在宜昌站以下均大于u,隨著洪峰持續時間的延長,宜昌站以下處于運動波傳播狀態的河段距離也隨之減少,當洪峰流量分別持續1,2,3,4 d和5 d時,計算波速分別在宜昌站以下550,670,740,810 km和880 km處進入到(u,2u)范圍內,在該范圍內的洪水波符合運動波傳播規律,其余上游河段計算波速大于2u,在絕大部分河段小于動力波速Cd,傳播特性符合過渡區洪水過程。當有區間支流入匯,計算波速表現出類似變化,所不同的是當斷波持續時間達到5 d時,宜昌站以下河段計算波速均大于2u,已無符合運動波傳播河段。

2.4.3" 多波束斷波傳播特征分析

在上游梯級水庫調度的過程中,水庫下游通常會產生系列的間斷波束,多波束的傳播增加了下游水流的非恒定特性。對比多波束條件下洪水波的傳播過程,發現間隔出現的斷波波束在沿程傳播過程中會不斷坦化融合,在宜昌站的階梯狀斷波流量過程傳播到螺山站時流量曲線已不再明顯彎折,漢口站及以下基本呈現單波束形態(圖9)。無論是漲水波還是退水波,隨著距離的增加斷波效應均逐漸減弱,洪峰減小,歷時延長,螺山、漢口、九江和大通站洪峰分別削減了2.24%,10.56%,16.85%和21.47%。當有支流入匯時,沿程站點洪峰流量變化表現出明顯的非一致性,支流水流的匯入不同程度地增加了支流入匯處下游站點的洪峰流量,且高流量歷時明顯延長。從表2可見,隨著沿程距離的增加洪水波形變指數遞減,表明波峰在不斷的削減,但多波束傳播條件下形變指數沿程衰減過程明顯放緩。

多波束斷波條件下計算波速結果如圖10所示,當無區間支流入匯,計算波速在宜昌站以下河段均大于2u,無符合運動波傳播河段。在宜昌站以下180 km內計算波速在動力波速Cd上下波動,該河段洪水運動符合動力波傳播規律,其余河段計算波速符合過渡區洪水運動規律。當有區間支流入匯,計算波速表現出類似變化。

3" 結 論

通過對長江宜昌至大通河段水庫泄流洪水波傳播的深入研究,本文構建了一維水動力學模型并進行了詳細分析。通過模擬對比不同水庫泄流洪水波與運動波和動力波的傳播規律,得出以下主要結論:

(1) 當水庫在短時間內下泄流量急劇變化時,即斷波形成后,隨著斷波流量增加,下游各站點的峰現時間不同程度提前,越靠下游峰現時間變化越大;斷波流量每增加5 000 m3/s,各站點峰現時間平均提前約7.5 h,當斷波流量增加到1 5000 m3/s時,宜昌至大通絕大部分河段洪水波運動符合過渡區洪水波特性,局部表現為動力波特性;隨著斷波流量的增加,同一站點洪水波形變指數急劇增加,并且斷波流量越大,形變指數衰減越慢。

(2) 隨著斷波流量持續時間的延長,下游各站點洪峰流量不同程度增加,越靠下游站點洪峰增加越大;同時斷波流量持續時間的延長也推遲了各站點的峰現時間,但峰現時間的影響是越靠下游站點影響越小;斷波流量持續時間越長,符合運動波傳播規律的河段距離越短;隨著洪峰持續時間的增加,洪水波的形變指數衰減速度變緩。

(3) 水庫下游產生的系列間斷波束增加了下游水流的非恒定特性,對比多波束條件下洪水波的傳播過程,發現間隔出現的斷波波束在沿程傳播過程中會不斷坦化融合,隨著距離的增加斷波效應均逐漸減弱,計算波速普遍較高,全河段無符合運動波傳播范圍,隨著沿程距離的增加洪水波形變指數遞減,但多波束傳播條件下形變指數沿程衰減過程明顯放緩。

需要指出的是,當有支流入匯時,以上3種條件下的沿程站點洪峰流量變化均表現出明顯的非一致性,支流水流的匯入不同程度地增加了支流入匯的下游站點的洪峰流量,且高流量歷時明顯延長。洪水波形變指數也隨著支流匯入的大小和沿程距離的變化而發生調整,但總的演化趨勢和無支流匯入情況一致。本研究不僅深化了對長江下游水庫泄流洪水波傳播規律的理解,也為實際工程和管理決策提供了重要的技術支持。未來的研究可以進一步結合更多實測數據和復雜水文條件,揭示斷波洪水的形成和演進機制,以應對日益復雜的洪水風險管理挑戰。

參考文獻:

[1]" MANUELA I B,PHILIPPE N.Spatial variability in Alpine Reservoir regulation:deriving reservoir operations fromstreamflow using generalized additive models[J].Hydrology and Earth System Sciences,2023,27(3):673-687.

[2]" 張建云,盛金保,金君良,等.全國水庫大壩應急管理存在問題和對策建議[J].中國應急管理科學,2022(9):23-30.

[3]" 盧程偉.流域水庫群蓄滯洪區綜合防洪調度研究與應用[D].武漢:華中科技大學,2019.

[4]" 胡瓊方,李秋平,鄒濤,等.三峽出庫急變洪水波傳播規律及其對水文測報的影響[J].水利水電快報,2019,40(2):1-5,25.

[5]" 程海云,陳力,許銀山.斷波及其在上荊江河段傳播特性研究[J].人民長江,2016,47(21):30-34.

[6]" 盧程偉,陳莫非,張余龍,等.斷波在朱沱—三峽壩址庫區河段傳播規律分析[J].長江科學院院報,2021,38(8):14-18,24.

[7]" 張茜.佛子嶺水庫泄洪洪水演進模擬研究[D].淮南:安徽理工大學,2020.

[8]" 王麗學,林鳳偉,汪可欣,等.基于蒙特卡羅模擬的泄洪風險率計算[J].人民長江,2008,39(19):20-22.

[9]" TENG J,JAKEMAN A J,VAZE J,et al.Flood inundation modelling:a review of methods,recent advances and uncertainty analysis[J].Environmental Modelling amp; Software,2017,90:201-216.

[10]王麗萍,張驗科,紀昌明,等.模擬最大熵法及其在水庫泄洪風險計算中的應用[J].水利學報,2011,42(1):27-32.

[11]高帥領,夏軍強,董柏良,等.雨水口泄流對城市洪澇影響的數學模型[J].浙江大學學報(工學版),2022(3):590-597.

[12]邴建平,鄧鵬鑫,張冬冬,等.三峽水庫運行對鄱陽湖江湖水文情勢的影響[J].人民長江,2020,51(3):87-93.

[13]徐長江,徐高洪,戴明龍,等.三峽水庫蓄水期洞庭湖區水文情勢變化研究[J].人民長江,2019,50(2):6-12.

[14]李曉昭.基于MIKE11的三峽庫區洪水波傳播規律研究[D].武漢:華中科技大學,2018.

[15]孫立剛,李國定,陳金海,等.“東方之星”輪翻沉事件再現仿真系統研究[J].中國水運(下半月),2016,16(6):78-80,153.

[16]程海云,陳力.三峽水庫泄水波與沙市站水位流量響應關系研究[J].人民長江,2017,48(19):29-34,41.

[17]DHI.A modelling system for rivers and channels:reference manual[R].Denmark:DHI Water amp; Environment,2004.

[18]張冬冬,戴明龍,徐高洪,等.三峽水庫蓄水期洞庭湖出湖水量變化[J].水科學進展,2019,30(5):613-622.

[19]程思學.基于KdV-Burgers型方程的非破碎涌波生成理論研究[D].杭州:浙江大學,2022.

[20]DIPSIKHA D,ANUPAL B,ARUP K S.Characterization of dam-impacted flood hydrograph and its degree of severity as a potential hazard[J].Natural Hazards,2022,112(3):1989-2011.

(編輯:郭甜甜)

Propagation law of reservoir flood wave in main channel of Changjiang River from Yichang to Datong section

QIU Hongya1,LIU Yaxin1,2,3,LI Yanqing4,ZHANG Dongdong4,XIONG Feng4

(1.China Three Gorges Corporation,Yichang 443133,China;

2.China Yangtze Power Co.,Ltd.,Yichang 443002,China;

3.Hubei Provincial Key Laboratory of Wisdom Yangtze River and Hydropower Science,Yichang 443002,China;

4.Bureau of Hydrology,Changjiang Water Resources Commission,Wuhan 430010,China)

Abstract:

In current water conservancy engineering and flood management,the propagation of flood waves from reservoir discharge is a matter of great concern.Taking the Yichang-Datong section in main channel of Changjiang River as an example,a one-dimensional hydrodynamic model was constructed to simulate the propagation characteristics of flood waves under different working conditions and to compare and analyze the flood evolution patterns along the course.The research results show that:① The formation of break wave floods can shorten the peak arrival time in downstream sections to varying degrees.With the increase of break wave flow,the distance of river sections where flood waves propagate as kinematic waves significantly shortens,and the flood wave deformation index increases;② With the extension of the duration of peak flood flow,the distance of flood wave propagation with kinematic wave characteristics will also be shortened,and the flood wave deformation index will increase.However,the peak arrival time at the same site or a certain mileage will be delayed,and the duration of the peak flood will be extended;③ Under multi-beam break wave conditions,the calculated wave speed is generally higher.As the distance along the course increases,the deformation index declines,but the decay process of the deformation index along the course under multi-beam propagation conditions is significantly slowed down;④ When there are tributary inflows,the changes in peak flow at sites along the course all exhibit significant inconsistencies,and the duration of high flow is significantly extended.The research findings are of great significance for revealing the propagation laws of flood waves from reservoir discharge and can provide references for regional flood management.

Key words:

reservoir release; flood wave; propagation characteristics; deformation index; hydrodynamic model; Changjiang River

主站蜘蛛池模板: 欧美精品在线看| 亚洲精品波多野结衣| 国产精品yjizz视频网一二区| 91一级片| 狠狠综合久久| 精品国产成人三级在线观看| 久久五月天综合| 欧美三级视频网站| 久热99这里只有精品视频6| 国产黄色片在线看| 伦伦影院精品一区| 日韩高清一区 | A级毛片无码久久精品免费| 亚洲国产欧美国产综合久久 | 国产日韩欧美在线视频免费观看| 欧美在线精品一区二区三区| 久久久久免费看成人影片| 97青草最新免费精品视频| 丝袜美女被出水视频一区| 久久精品娱乐亚洲领先| 国产丝袜一区二区三区视频免下载| 国产中文一区二区苍井空| 少妇人妻无码首页| 麻豆AV网站免费进入| 91精品国产自产91精品资源| AV网站中文| 在线播放91| 国产成人高清亚洲一区久久| 国产香蕉在线| 国产91久久久久久| 国产乱人伦精品一区二区| 亚洲精品第一在线观看视频| 亚洲成人黄色在线观看| 国产真实乱子伦视频播放| 嫩草在线视频| 麻豆精选在线| 欧美日韩成人在线观看| 免费a级毛片视频| 国产成人精品三级| 男人天堂亚洲天堂| 无码国内精品人妻少妇蜜桃视频| 日本精品中文字幕在线不卡| 一区二区欧美日韩高清免费| 丝袜亚洲综合| 国产高清在线精品一区二区三区| 亚洲无线一二三四区男男| 2020极品精品国产| 久久久久亚洲精品无码网站| 成人免费视频一区| 无码电影在线观看| 国产成人区在线观看视频| 一本大道视频精品人妻| 香蕉伊思人视频| 在线观看免费黄色网址| 亚洲成a人片77777在线播放| 小说区 亚洲 自拍 另类| 在线观看国产精美视频| 国产v精品成人免费视频71pao| 免费无码AV片在线观看国产| 88av在线| 成年av福利永久免费观看| 91精品国产自产在线老师啪l| 亚洲高清在线天堂精品| 亚洲天堂在线免费| 亚洲AV人人澡人人双人| 激情综合婷婷丁香五月尤物 | 日韩第九页| 2021国产v亚洲v天堂无码| 国产一区亚洲一区| 久久a毛片| 国产JIZzJIzz视频全部免费| 婷婷色中文网| 亚洲天堂区| 99久久国产精品无码| 九色最新网址| h视频在线播放| 午夜精品一区二区蜜桃| 国产极品嫩模在线观看91| 欧美日韩一区二区三区在线视频| 欧美色视频在线| 亚洲无码视频喷水| 麻豆国产在线观看一区二区 |