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

三維水沙模型在永定河河勢分析中的應用

2014-08-08 01:00:17韓鳳霞何國建方紅衛
水利水電科技進展 2014年6期
關鍵詞:模型

韓鳳霞,何國建,方紅衛

(1.北京市水利規劃設計研究院,北京 100048;2.清華大學水利水電工程系,北京 100084;3.清華大學水沙科學與水利水電工程國家重點實驗室,北京 100084)

三維水沙模型在永定河河勢分析中的應用

韓鳳霞1,何國建2,3,方紅衛2,3

(1.北京市水利規劃設計研究院,北京 100048;2.清華大學水利水電工程系,北京 100084;3.清華大學水沙科學與水利水電工程國家重點實驗室,北京 100084)

為了檢驗永定河現狀治導工程能否有效控制河勢,采用CREST三維水沙數值模擬系統模擬分析河勢走向、水沙濃度以及河床變形。根據永定河盧梁段河流泥沙和河床變化特征,引入合理床面懸移質平衡輸沙濃度、推移質平衡輸沙率以及平面分配系數,建立永定河盧梁段三維水沙模型并進行河床變形驗證;模擬分析了河道流速場分布、懸沙濃度分布及河床變形規律,以及現狀治導工程在控導河勢方面發揮的作用,得到了河道環流、丁壩頭繞流等復雜水流規律,并對游蕩型河道流勢、流向和頂沖點位置以及洪水主槽沖刷、邊灘丁壩保護區淤積情況進行了預測。分析結果表明,河道治導工程在控制河勢方面依然起到主導作用,部分丁壩壩頂受漫流沖刷破壞,局部河道需要添加少量短丁壩以鞏固堤防安全。

三維水沙模型;平衡輸沙率;河勢;河床變形;治導工程;永定河

絕大多數天然河道流動具有三維特性,尤其工程近區泥沙運動及物質輸移具有強烈的三維特性,因此,僅僅考慮一維和二維數學模型已不能滿足工程實踐要求。隨著計算機及網絡技術的發展、相關數學模型理論的完善以及工程實踐要求的提高,三維數學模型在水利、泥沙和環境領域的應用更加廣泛,也逐步走向成熟。

泥沙三維數學模型起步于20世紀80年代,其研究精度在很大程度上依賴于河床表面條件的選取,90年代以后,三維水沙模型得到了迅速發展。Shimizu等[1]在假設垂線方向流速為對數分布的基礎上建立了簡化三維模型,用于彎道水槽動床試驗計算中,并將結果與二維模型計算值進行了比較分析;Olsen等[2-3]模擬了圓柱繞流的最大沖刷坑問題,利用k-ε紊流模型求解N-S方程,得到床面剪切力,然后根據床面剪切力與床面上泥沙濃度關系計算泥沙輸移擴散方程,并利用該方法模擬了沉沙池床面變化情況;Wang等[4]發展了一個三維泥沙數學模型,該模型采用有限元離散,計算過程采用靜壓假定,模型可以模擬沖積河流中水流、泥沙運動和河床地形變化,特別適用于模擬橋墩附近局部沖刷問題;Gessler等[5]利用美國工程兵團的三維模型CH3DSED計算了密西西比河下游彎道河流泥沙輸運過程,預測了深挖航道工程對河道沖淤的影響;Li等[6]利用考慮了浮力作用的k-ε紊流模型計算了泥沙在不同水流條件下的沉積形態;Wu等[7]在Zhu[8]開發的水動力學模型基礎上,添加泥沙輸運和河床變形計算,建立了一個河道水沙運動三維數學模型,該模型使用了基于二維泊松方程求解自由水面,通過求解懸沙輸運方程計算懸沙濃度,推移質輸運采用非平衡輸沙模式計算方法,通過泥沙連續方程計算出河床變形。

國內在三維水沙模型研究上也做了大量工作。Fang等[9-10]考慮了非飽和懸移質計算模式和全沙模式,并在永定河河床變形模擬中對平衡輸沙和非平衡輸沙模式進行了比較[10];方紅衛[11]研究了天然水沙條件下變密度紊流在一般曲線坐標系中的基本方程和邊界條件,并計算了浮運沉井施工期三維沖刷問題;Fang等[12]利用有限分析方法求解非正交曲線網格上水流泥沙輸運方程,模擬了三峽工程壩區1976年以來水庫沖淤規律;陸永軍等[13]根據紊流隨機理論,利用壁函數模擬了葛洲壩水利樞紐修建前后水沙變化以及三峽工程壩區泥沙沖淤規律;夏云峰[14]建立了基于k-ε模式的三維水沙模型,模擬了天然河道長江南京大勝關至西壩分汊河段水流流速分布、泥沙濃度分布及汊道沖淤情況等;假冬冬等[15]建立了考慮黏性河岸變形的三維水沙模型,將崩岸模型與局部網格可動技術相結合,重現了彎道實驗中河道擺動過程;吳修廣等[16]基于FVCOM模型模擬了杭州灣大潮期間泥沙輸運過程,對涉水建筑物附近三維水沙影響進行了分析評估。

三維水沙模型是完整模擬流體運動以及泥沙輸運和河床變形的主要工具,但由于目前水沙輸移理論體系尚不完善,還缺乏普遍適用的統一標準,不同研究者建立的三維水沙模型各有特點,數值求解方法各不相同。這些模型在一定條件下針對工程問題都可以取得滿足要求的結果,具有一定的理論價值和工程應用價值,同時又都存在一定的局限性。因此,對已有模型進行分析,不斷補充完善,提高精度,擴大應用范圍是發展三維水沙模型研究的有效途徑。本文是在前人工作基礎上,針對研究區域水沙特征,選取合適的三維水沙模型及離散方法,對永定河盧溝橋至梁各莊段水沙輸移及河床變形進行模擬,分析百年一遇洪水條件下,水面變化、流勢分布以及河床變化情況,為永定河河勢分析提供技術支撐。

1 研究區域特征

永定河是對北京市防洪威脅最大的一條河,為季節性河流,夏季洪水陡漲陡落,含沙量高,尤其中游河道具有沖刷和強烈堆積并存的游蕩性。永定河在北京境內的盧溝橋至梁各莊段河道(以下簡稱盧梁段)長約60km,河道兩岸堤距220~2600m,寬窄變化較大,河床粒徑約2mm,懸沙中值粒徑約為0.05mm。由于兩岸堤防及河床多為砂(沙)質,沖淤不穩,中泓游蕩,險工往往上提下錯,洪水期間險情迭出,歷來是北京市防汛重點之一。

為了治理永定河洪災,1959—1961年根據盧梁段河道特點以及“三固一束”規劃治理原則(固定險工,以改善和解決永定河防洪問題;固定河槽,以保障行洪順暢;固定灘地,以防止灘地顯沒無常;束窄河道,以使河槽逐漸刷深)實施了治導工程。為了控制河勢和保護堤防,沿河道兩側修建了大量丁壩和丁壩群。圖1為該河段中金門閘至閻家鋪段全長約13km的地形圖,兩岸共有丁壩26座(圖中長短實線為丁壩)。由圖1可見,該段河道寬窄不一,深槽和淺灘相互交錯,丁壩群與單個長丁壩相互呼應,地形極為復雜,水沙輸移具有強烈的三維特性。

圖1 永定河計算區域地形和斷面布設

2 模型原理及主要參數

本研究水沙模擬計算采用清華大學研制的CREST三維水沙數值模擬系統。三維水沙模型關鍵在于選取合適的泥沙輸移模型及邊界條件,水動力學模塊采用標準N-S方程耦合k-ε紊流模型求解,在同位網格下利用SIMPLE算法求解。

2.1 懸移質輸移

研究河段內懸沙濃度較小,不足以影響水流密度,因而在懸移質控制方程中不考慮密度變化。懸移質泥沙濃度控制方程如下:

式中:σc為紊動Schmidt數;uj為水流三維方向上的流速;ωs為泥沙顆粒沉降速度。

水面處懸移質通量為零,其濃度cn由下一層網格點濃度cn-1通過指數關系求得。床面上泥沙濃度作為求解的一個邊界條件,是正確求解懸沙濃度場的重要因素;已知床面上一層單元網格點含沙量c2,得到床面附近含沙量cb:

式中:cb*為床面附近的懸移質平衡輸沙濃度;b為床面泥沙交換層厚度,Einstein認為b=2D50[17],本文在河道計算中取為0.02倍水深,即b=0.02h。

式中:f(η)為含沙量沿垂線分布函數;ηa和ηh分別為距離床面和水面相對位置。

2.2 推移質輸移

清水沖刷河床粗化過程中,由于進口斷面沒有來沙,水流就要從河床中補充泥沙,沖刷河床,水流和河床之間通過泥沙交換趨向平衡,使得下游斷面輸沙率要大于上游斷面輸沙率,這就是推移質不平衡輸移的表現。將推移質不平衡輸沙方程推廣到平面二維情況:

式中:qb和qb*分別為推移質輸沙率和平衡輸沙率;qbx和qby為qb在x、y方向的分量;αbx和αby分別為垂線平均流速和河床切應力的方向系數;Ls為不平衡輸沙長度,取為

對推移質平衡輸沙率qb*,各研究學者在不同理論基礎上提出了大量計算公式,具有代表性的有Meyer-Peter公式、Bagnold公式、Einstein公式、Yalin公式、Engelund公式以及Achers-White公式等;錢寧等[20]對各家計算公式進行了分析統計和說明。根據公式適用范圍以及數值計算的可操作性,本文采用van Rijn提出的公式[18]:

式中:d50為中值粒徑;u*c為床面臨界剪切速率,可由Shields曲線獲得。

2.3 河床變形

根據質量守恒原理,某一點河床變形等于該點處泥沙進出通量差:

其中δb為床面粗糙層高度。

3 模型驗證及結果分析

治理永定河盧梁段這樣的游蕩型河道,主要任務就是分析河勢變化情況,并采取適當措施進行控制;只有河勢得到了一定的控制,主流才能比較穩定,游蕩程度才會有所減弱。最近幾十年,永定河河道上游來水來沙條件發生了較大變化,為了檢驗現狀丁壩能否有效控制河勢,并分析判斷是否需要新設丁壩鞏固堤防安全,根據前面建立的三維水沙數學模型(以下簡稱本文模型),對該河段河勢走向、水沙濃度以及河床變形進行模擬。由于自永定河治導工程實施以來,永定河基本未再發生較大洪水,無實際資料可供驗證;本研究利用北京市水利科學研究所完成的動床河勢物理模型試驗對模型進行率定驗證[22],該物理模型試驗平面比尺為500,垂直比尺為60,采用電木粉作為模型沙,密度為1.45t/m3,中值粒徑為0.1mm。對應物理模型試驗條件,選用標準洪水洪峰流量2500m3/s、來水懸移質濃度3kg/m3作為邊界條件。

3.1 河床變形驗證

圖2為沿程4個斷面本文模型沖淤計算值與北京市水利科學研究所物理模型試驗值對比,可見計算值與試驗值擬合較好,表明本文模型能夠反映河床變形的實際情況。

3.2 流速場分布

該河段河勢的控制、險工位置的固定,除了有賴于護岸工程外,還與灘地能否得到保護有關。長短丁壩所形成的“凸出”和“凹入”型險工分別具有以下特點:凸出型險工外形比較顯著地凸出河中,有著顯著的挑流作用,能將主流挑向對岸預定部位;如果靠溜長度較長,則水流能順著險工岸線流動,出水方向比較固定,還能起到導流作用。凹入型險工外形為凹入的弧形,有如彎道凹岸;這類險工經常靠溜,且靠溜范圍較長,變化小,水流經過險工后,出水方向頗為穩定,能順著凹岸弧線平順地流出,既能迎托水流,又能導引水流,在控制河勢方面比平順型和凸出型相對都好。

通過本文模型計算得到該區域丁壩控制下的流速和流線分布如圖3所示,可以看出,長丁壩壩身較長而突入河中,挑托主流離開堤岸,以掩護下游堤岸不受沖刷;短壩略呈下挑式,迎托水流,消減水勢,減輕急流對堤岸沖刷。長丁壩起到控制河勢作用,主流被控制在河道中間,險工處短丁壩群則起到保護堤岸的作用。圖4為局部河道流場圖,可見長丁壩控制河勢處流速較大,特別是丁壩頭部繞流流速大;丁壩群附近、兩個丁壩之間水流形成漩渦。

3.3 懸沙濃度分布

根據懸沙輸運對流擴散方程,即可進一步計算出懸沙濃度的三維分布規律。本文使用的模型模式已經成功用于多種復雜區域的泥沙輸移及河床變形計算,如三峽水庫、德國易北河等(具體見文獻[9,12]等),從河床變形的數據對比也可看出本文模型適用于永定河的水沙計算。圖5分別為沿垂向深度0.1H、0.4H、0.7H(H為總水深)和水體表面4層的懸沙濃度(質量濃度)平面分布。

圖2 河床變形計算值與試驗值對比

圖3 流量2500m3/s時河道流速和流線分布

圖4 局部河道流場

圖5 不同深度懸沙濃度平面分布

從圖5可以看出,在垂向分布上,由于重力和濃度梯度擴散的共同作用,懸沙濃度從水體表層到水體底層濃度逐漸增加,符合懸移質泥沙輸移一般規律。同時從灘、槽部位泥沙分布來看,主槽是泥沙輸移主要區域,如斷面2,其灘面泥沙濃度不到主槽泥沙濃度一半,主槽中由于河勢影響還伴有螺旋流輸沙現象。

從同水深處平面懸沙濃度分布可見:無論是在長丁壩還是短丁壩群挑流處,由于丁壩控制主流區域流速較大,易發生沖刷,使水體含沙量較高,而在受長丁壩保護的邊灘及短丁壩群內有較多漩渦,泥沙易發生淤積,水體中含沙量較小,這也是河流洪水期典型的刷槽淤灘自然現象。

3.4 河床變形

圖6為洪水期間整個區域河床沖淤分布圖。由圖6可見:河床沖刷主要發生在兩處,一處是河道主槽區域,除了在下游局部河段外,沿程主槽均發生沖刷,平均沖刷接近50cm;第二處為丁壩頭部,該處由于挑流作用導致局部流速較大,長丁壩頭處沖深均在1m以上。因此永定河發生2 500m3/s的洪水(相當于百年一遇)情況下,在做好堤防防護的同時,還應注意丁壩頭的保護。從圖6還可見:受丁壩保護的邊灘區域泥沙落淤,淤積厚度一般在30cm以內;這些淤積物保護了堤岸使其免受沖刷,但淤積在灘面的泥沙一般難以重新被沖入主槽;洪水后期如果主槽也發生淤積,則將減少河道過流斷面面積,對防洪不利。

圖6 河床沖淤分布

由上述水流流速分布、泥沙濃度分布等計算結果分析可知,目前該河段丁壩及丁壩群在控制河勢方面仍然起到了主導作用,與過去相比河勢規律沒有發生大變化,長丁壩和短丁壩群仍然是主流靠溜處,丁壩依然是保護堤防和導引水流的主要水工建筑物。但在某些局部河段,由于洪水來流能量較大,靠近邊岸處仍會發生沖刷現象,例如圖6中a、b兩處,需要添加少量短丁壩,以起到保護該處堤防的作用。

4 結論

a.在永定河河勢變化分析上,獲得了河道橫流、斜流以及平面環流、三維環流、丁壩頭繞流等復雜水流規律,預測分析了游蕩型河道流勢、水流流向和頂沖點位置等。

b.在泥沙輸運上,得到了平面及垂向懸移質濃度的分布規律,結合河床沖淤特征,再現了洪水主槽沖刷、邊灘丁壩保護區淤積的現象。

c.模型模擬的流速場和流勢等分析結果顯示,該段河道治導工程在控制河勢方面依然起到主導作用,部分丁壩壩頂受漫流沖刷破壞,局部河道需要添加少量短丁壩以鞏固堤防安全。

[1]SHIMIZU Y,YAMAGUCHI H,ITAKURAT.Threedimensional computation of flow and bed deformation[J].Journal of Hydraulic Engineering,1990,116(9):1090-1109.

[2]OLSEN N R B,KJELLESVIG H M.Three-dimensional numerical flow modeling for estimation of maximum local scour depth[J].Journal of Hydraulic Research,1998,36(4):579-590.

[3]OLSEN N R B,KJELLESVIG H M.Three-dimensional numerical modeling of bed changes in a sand trap[J].Journal of Hydraulic Research,1999,37(2):189-198.

[4]WANGSY,JIAY.Computationalmodellingand hydroscience research:advanceinhydro-scienceand engineering[C]//Proceedingsof2ndInternational Conference on Hydro-Science and Engineering.Beijing:Tsinghua University Press,1995:2147-2157.

[5]GESSLER D,HALL B,SPASOJEVIC M,et al.Application of 3D mobile bed,hydrodynamic[J].Journal of Hydraulic Engineering,1999,125(7):737-749.

[6]LI C W,MA F X.3D numerical simulation of deposition patterns and sand disposal in flowing water[J].Journal of Hydraulic Engineering,2001,127(3):209-218.

[7]WU W M,RODI W,WENKA T.3-D calculations of flow and sediment transport in rivers[J].Journal of Hydraulic Engineering,2000,126(1):4-15.

[8]ZHU J.An introduction and guide to the computer program FAST3D[R].Karlsruhe,Germany:Unversityof Karlsruhe,1992.

[9]FANG H W,RODI W.Three-dimensional calculations of flow and suspended sediment transport in the neighborhood of the dam for the Three Gorges Project(TGP)reservoir in the Yangtze River[J].Journal of Hydraulic Research,2003,41(4):379-394.

[10]FANG H W,HE G J,WANG L X.Influence of vertical resolution and nonequilibrium model on three-dimensional calculations of flow and sediment transport[J].Journal of Hydraulic Engineering,ASCE,2010,136(2):122-128.

[11]方紅衛.浮運沉井施工期三維沖刷的數值模擬[J].中國科學:A輯,1995,25(1):107-112.(FANG Hongwei.Numerical simulationofthree-dimensionalerosionin floating caisson construction[J].Science in China:Ser A,1995,25(1):107-112.(in Chinese))

[12]FANG H W,WANG G Q.Three-dimensional mathematical model of suspended-Sediment transport[J].Journal of Hydraulic Engineering,2000,126(8):578-592.

[13]陸永軍,竇國仁,韓龍喜,等.三維紊流懸沙數學模型及應用[J].中國科學:E輯,2004,34(3):311-328.(LU Yongjun,DOUGuoren,HANLongxi,etal.Threedimensional mathematical model of suspended load and its application[J].Science in China:Ser E,2004,34(3):311-328.(in Chinese))

[14]夏云峰.感潮河道三維水流泥沙數值模型研究與應用[D].南京:河海大學,2002.

[15]假冬冬,邵學軍,王虹,等.考慮河岸變形的三維水沙數值模擬研究[J].水科學進展,2009(3):311-317.(JIA Dongdong,SHAOXuejun,WANGHong,etal.3D mathematical modeling for fluvial processes considering bank erosion[J].Advances in Water Science,2009(3):311-317.(in Chinese))

[16]吳修廣,劉光生,程文龍.基于FVCOM的杭州灣三維泥沙數值模擬[J].水利水運工程學報,2011(4):86-96.(WU Xiuguang,LIU Guangsheng,CHENG Wenlong.3D numerical simulation of sediment in the Hangzhou Bay based on FVCOM[J].Hydro-Science and Engineering,2011(4):86-96.(in Chinese))

[17]錢寧,萬兆惠.泥沙運動力學[M].北京:科學出版社,1983.

[18]van RIJN L C.Sediment transport,part I:bed load transport[J].Journal of Hydraulic Engineering,1984,110(10):1431-1456.

[19]王興奎,邵學軍,李丹勛.河流動力學基礎[M].北京:中國水利水電出版社,2002.

[20]錢寧,萬兆惠.泥沙運動力學[M].北京:科學出版社,1983.

[21]FANG H W.Three dimensional calculations of flow and bed load transport in the Elbe River(Report No.763)[R].Karlsruhe,Germany:University of Karlsruhe,2000.

[22]永定河盧溝橋下游河勢動床模型試驗[R].北京:北京市水利科學研究所,1991.

Apptication of three-dimensional flow-sediment model in flood control planning of Yongding River

//HAN Fengxia1,HE Guojian2,3,FANG Hongwei2,3(1.Beijing Institute of Water,Beijing 100048,China;2.Department of Hydraulic Engineering,Tsinghua University,Beijing 100084,China;3.State Key Laboratory of Hydroscience and Engineering,Tsinghua University,Beijing 100084,China)

In order to examine whether the training structures could effectively control the river regime or not,a threedimensional(3D)model CREST is applied to analyze the river regime trend,sediment transportation,and river bed deformation of Yongding River.According to the characteristics of sediment and riverbed variation of Lu-Liang section of Yongding River,a 3D numerical model is set up with a proper equilibrium concentration of suspended load,the bedload transport rate,and plane partition coefficient.The model is calibrated by bed deformation.Then,the rule of training structures played in controlling the river regime is studied by simulating the velocity field,the sediment concentration distribution and bed deformation.As a result,complex flow laws of the vortex and the spur dike flow are obtained.Additionally,the flow tendency,flow direction,scouring point location,erosion in main channel due to flood,and deposition owing to spur dike protection are predicted.The 3D numerical model is applied to simulate the river regime variation,sediment transport,and bed deformation of Yongding River.The results show that long spur dikes control the river regime,and some of them are suffered from overflow and local scour on the head.For levee safety,short groins need to be added to some local river banks where the erosion often happened.

three-dimensional flow-sediment model;equilibrium transport rate;river regime;riverbed variation;training structure;Yongding River

TV143+.4

:A

:1006-7647(2014)06-0056-06

10.3880/j.issn.1006-7647.2014.06.012

2013-10-12 編輯:熊水斌)

韓鳳霞(1967—),女,寧夏石嘴山人,高級工程師,主要從事水利規劃研究。E-mail:fengxia_han@sina.com

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 色AV色 综合网站| 亚洲国模精品一区| 日本午夜在线视频| 国产91丝袜在线播放动漫 | 亚洲永久色| 麻豆精品国产自产在线| 国产91九色在线播放| 久久人午夜亚洲精品无码区| 国产黄视频网站| 亚洲国产精品不卡在线| 色综合手机在线| 国产00高中生在线播放| 伦伦影院精品一区| 亚洲娇小与黑人巨大交| 国产综合色在线视频播放线视 | 丁香婷婷久久| 永久成人无码激情视频免费| 四虎精品黑人视频| 特级毛片免费视频| 亚洲熟女中文字幕男人总站| 日韩av资源在线| 国产乱子伦无码精品小说| 免费无遮挡AV| 亚洲成在人线av品善网好看| 国产午夜看片| 久久伊人操| 久久精品无码国产一区二区三区 | 欧洲在线免费视频| 四虎成人免费毛片| 久久免费看片| 操操操综合网| 国产人成在线视频| 精品一区二区无码av| 黄色成年视频| 污视频日本| 中文字幕亚洲另类天堂| 伊人久久久久久久| 最新国产高清在线| 久久无码av一区二区三区| 亚洲国产欧美国产综合久久 | 日韩精品无码免费一区二区三区 | 久久久久国产精品免费免费不卡| 亚洲精品无码av中文字幕| 久久精品aⅴ无码中文字幕| 国产在线拍偷自揄观看视频网站| 国产亚洲精品自在线| 一级爆乳无码av| 一级毛片免费高清视频| 天天综合网站| 亚洲精品在线影院| 四虎成人免费毛片| 国产白浆一区二区三区视频在线| 欧美色视频网站| 色综合久久无码网| 精品无码专区亚洲| 国产精品免费p区| 久草视频精品| 欧美日韩国产在线观看一区二区三区 | 99久久亚洲综合精品TS| 亚洲精品视频免费看| 亚洲网综合| 亚洲69视频| 久久久久人妻一区精品色奶水| 亚洲无码电影| 免费无码在线观看| 五月天福利视频| 99久久精品久久久久久婷婷| 免费国产无遮挡又黄又爽| 99热这里只有免费国产精品 | 色天堂无毒不卡| 九色视频线上播放| 国产成人毛片| 无码人妻免费| 亚洲成在人线av品善网好看| 亚洲天堂网在线观看视频| 国产精品美女网站| 精品国产福利在线| 精品久久香蕉国产线看观看gif| 国产激情影院| 亚洲另类第一页| 色哟哟色院91精品网站| 在线va视频|