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

        ?

        堰塞壩漫頂潰決過(guò)程數(shù)值模擬

        2012-04-17 09:30:22鐘啟明陳生水趙聯(lián)楨

        鐘啟明,陳生水,趙聯(lián)楨,任 強(qiáng),曹 偉

        (1.南京水利科學(xué)研究院巖土工程研究所,江蘇南京 210029;

        2.南京水利科學(xué)研究院水文水資源與水利工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,江蘇南京 210029;

        3.南京林業(yè)大學(xué)土木工程學(xué)院,江蘇南京 210037)

        地震、降雨等自然災(zāi)害往往引發(fā)大規(guī)?;潞蜕襟w崩塌,使受災(zāi)區(qū)域河流阻塞,形成天然堰塞壩。堰塞壩的工作條件、壩體幾何特征以及壩體物質(zhì)組成和內(nèi)部結(jié)構(gòu)都與人工土石壩存在明顯差別。首先在工作條件上,因沒(méi)有溢洪道和其他泄洪設(shè)施,上游持續(xù)來(lái)水使壩前水位快速上漲,堰塞壩很容易發(fā)生漫頂潰決;其次從壩體形態(tài)來(lái)看,堰塞壩壩頂一般凹凸不平,存在天然凹槽,且體型多呈不規(guī)則形狀,沿河流方向堆積范圍明顯超過(guò)狹窄河谷的寬度;再者從壩體結(jié)構(gòu)來(lái)講,因壩體材料沒(méi)有經(jīng)過(guò)人工篩選和施工碾壓,大部分堰塞壩結(jié)構(gòu)疏松,不均勻性強(qiáng),壩體材料多為礫石、粗顆粒,粗細(xì)混雜,粒徑級(jí)配范圍變化大[1-2]。調(diào)查表明[3],90%的堰塞壩都在其形成1a內(nèi)潰決。一旦發(fā)生潰決,所產(chǎn)生的洪水或泥石流將對(duì)下游產(chǎn)生嚴(yán)重危害。因此,有必要研究堰塞壩漫頂潰決機(jī)理及潰決過(guò)程,為正確評(píng)估堰塞壩潰決致災(zāi)后果、科學(xué)制定堰塞壩潰決防洪應(yīng)急預(yù)案提供技術(shù)支撐。

        1 堰塞壩潰決機(jī)理和潰決過(guò)程分析

        目前關(guān)于堰塞壩潰決機(jī)理和潰決過(guò)程的研究工作主要參考人工土石壩的成果,對(duì)于如何合理模擬寬級(jí)配堰塞壩材料在高強(qiáng)度水流侵蝕條件下的非平衡沖刷特性的研究工作很少涉及。王光謙等[4]提出在堰塞壩潰壩模型中考慮潰口水力沖刷和重力坍塌擴(kuò)展,但是該模型使用的沖蝕公式是在恒定均勻流條件下建立的一般平衡沖蝕公式;黃金池[5]利用水槽試驗(yàn)結(jié)果,將高強(qiáng)度水流非平衡沖刷公式引入到潰口擴(kuò)展模型,并將潰口垂直下切、橫向擴(kuò)展、壩坡溯源沖刷3種過(guò)程結(jié)合起來(lái)建立堰塞壩潰壩數(shù)學(xué)模型,但是該模型的沖蝕公式主要基于室內(nèi)平面水槽試驗(yàn),泥沙粒徑范圍為0.18~13.00mm,無(wú)法合理反映堰塞壩邊坡上礫石、塊石等大粒徑滑坡堆積物的起動(dòng)問(wèn)題,且沖蝕公式中的參數(shù)不易獲取。

        要正確模擬堰塞壩漫頂潰決過(guò)程,需要從堰塞壩材料粒徑范圍寬的特點(diǎn)入手,同時(shí)考慮壩體密實(shí)度、壩體材料強(qiáng)度、下游坡度、漫頂水流流速的影響。筆者基于上述考慮,根據(jù)堰塞壩壩體材料的寬級(jí)配特性,建立描述堰塞壩漫頂潰決潰口發(fā)展規(guī)律與流量過(guò)程的數(shù)值模型,并利用該模型對(duì)唐家山堰塞壩的潰決過(guò)程進(jìn)行數(shù)值模擬,驗(yàn)證了該模型和計(jì)算方法的合理性。

        2 堰塞壩漫頂潰決過(guò)程數(shù)值模型

        2.1 寬級(jí)配壩體顆粒臨界起動(dòng)流速

        堰塞壩漫頂潰決過(guò)程中水流的沖蝕過(guò)程以及非均勻巖土顆粒的起動(dòng)具有非恒定的特性。假設(shè)水流作用下忽略顆粒間的黏性特征,首先根據(jù)堰塞壩壩體土石料顆粒的受力情況求得顆粒起動(dòng)的臨界流速(圖1)。

        圖1 土體顆粒在壩坡上的受力示意圖Fig.1 Forces acting on a soil particle along dam slope

        對(duì)土體代表顆粒1而言,考慮水流作用下巖土體顆粒所受的力一般有重力W、水流拖曳力Fd[6]、上舉力Fl[6],其表達(dá)式分別為

        式中:ρs——土顆粒的密度;ρw——水的密度;d50——土顆粒的平均粒徑;g——重力加速度;Cd——拖曳力系數(shù),一般取0.4[7];v——水流流速;Cl——上舉力系數(shù),一般取0.1[7]。

        需要指出的是,堰塞壩體土石料級(jí)配寬泛,最大粒徑和最小粒徑相差很大。為此,引入與水流方向垂直的附加作用力R(圖1)來(lái)考慮粗顆粒對(duì)細(xì)顆粒的阻攔、遮蔽作用以及細(xì)顆粒對(duì)粗顆粒的包圍、填實(shí)作用。由于目前對(duì)R的研究很少,尚無(wú)精確理論公式,可近似假定R與顆粒的平均剪力成比例[8-11],即

        其中

        式中:R——附加作用力;φ——比例系數(shù)與顆粒面積系數(shù)的乘積;τs——不均勻顆粒的平均剪力;Km——無(wú)因次系數(shù);M——緊密系數(shù),代表顆粒組成的密實(shí)程度,反映顆粒組成的不均勻程度和細(xì)顆粒填充孔隙的程度[9];Cu——不均勻系數(shù)。

        由式(4)可得

        其中

        式中K為系數(shù),可以根據(jù)不均勻顆粒起動(dòng)流速試驗(yàn)資料及天然河流實(shí)測(cè)資料確定,在0.785~1.727范圍內(nèi)變化,此處取1.3[9]。

        如圖1所示,土顆粒起動(dòng)時(shí)受到的摩擦力可表示為

        式中:Ff——土顆粒受到的摩擦力;φ——土顆粒間的內(nèi)摩擦角;θ——壩坡坡角;c——土體的黏聚力。通過(guò)受力分析可知,土顆粒1起動(dòng)的臨界條件為

        將式(1)(2)(3)(5)(6)代入式(7)可得土顆粒在壩坡上的臨界起動(dòng)流速:

        2.2 臨界起動(dòng)流速公式驗(yàn)證

        唐家山堰塞壩潰決過(guò)程中,現(xiàn)場(chǎng)測(cè)量得到壩體顆粒的起動(dòng)流速為2.4m/s[12],因此,可利用該起動(dòng)流速來(lái)驗(yàn)證本文提出的臨界起動(dòng)流速公式的合理性。根據(jù)相關(guān)資料[13],計(jì)算參數(shù)取值為:d50=0.03m,c=60kPa,M=0.75,θ=23°,φ=45°,g=9.8m/s2,ρs=2.6×103kg/m3。將上述參數(shù)代入式(8),可計(jì)算得出臨界起動(dòng)流速vc=2.12m/s(如不考慮附加力作用,臨界起動(dòng)流速為1.35m/s),與實(shí)測(cè)結(jié)果基本一致,驗(yàn)證了寬級(jí)配顆粒起動(dòng)流速的合理性。

        2.3 壩頂溢流沖刷速率

        當(dāng)土體顆粒在水流作用下起動(dòng)后,在潰壩水流的作用下,壩頂潰口和下游壩坡發(fā)生沖蝕,由于壩體材料粒徑級(jí)配范圍變化大,因此在分析不同土體陡水槽沖蝕試驗(yàn)結(jié)果[14-16]的基礎(chǔ)上,選擇d90與d30為代表粒徑,筆者建議了一個(gè)計(jì)算堰塞壩土石料沖蝕率的經(jīng)驗(yàn)公式:

        其中

        式中:Qs——水流沖蝕率;d90——小于某粒徑土石料質(zhì)量分?jǐn)?shù)為90%所對(duì)應(yīng)的顆粒粒徑;d30——小于某粒徑土石料質(zhì)量分?jǐn)?shù)為30%所對(duì)應(yīng)的顆粒粒徑;B——潰口寬度;v*——摩阻流速;vb——潰口底流速;ˉv——水流平均流速;H——堰塞湖水位高程;Hc——潰口底部高程;J——水力梯度;N——潰口處糙率;Qb——潰口流量。

        當(dāng)漫頂潰壩發(fā)生后,水流沿著初始潰口沖蝕下游壩坡,潰口流量可采用以下寬頂堰公式計(jì)算:

        式中m為流量系數(shù),此處可取0.5[17]。

        2.4 沖蝕公式驗(yàn)證

        建議的壩體材料沖蝕公式可反映筑壩材料顆粒級(jí)配、密實(shí)度、強(qiáng)度、滲透通道傾角和壩坡、摩阻流速及水流速度對(duì)沖蝕量的影響,公式中各物理量意義明確且易于測(cè)定。為了驗(yàn)證所建議公式的合理性,運(yùn)用該公式對(duì)4種不同土體的單寬沖蝕量進(jìn)行計(jì)算分析,計(jì)算結(jié)果與陡水槽沖蝕試驗(yàn)結(jié)果[14-16]的對(duì)比見(jiàn)圖2。從圖2可以看出,計(jì)算結(jié)果位于試驗(yàn)結(jié)果的變化區(qū)域內(nèi),這表明從工程應(yīng)用角度來(lái)看,所建議的壩體材料沖蝕公式是可以接受的。

        圖2 水流流速與單寬沖蝕率的關(guān)系Fig.2 Relationship between flow velocity and erosion rate

        2.5 潰口發(fā)展

        土石壩漫頂破壞潰口發(fā)展?fàn)顩r見(jiàn)圖3。壩體初始潰口形狀并非一直保持不變,而是隨漫頂水流流速的增大和壩體沖蝕的加劇不斷加深加寬,時(shí)間段增量Δti內(nèi)水流下切深度增量為

        時(shí)間段Δt內(nèi)水流下切深度增量為

        式中j為Δt時(shí)段內(nèi)的時(shí)間步長(zhǎng)數(shù)。

        如果不考慮潰口邊坡失穩(wěn)和坍塌引起的潰口橫向擴(kuò)展,潰口底部的沖蝕速率應(yīng)和潰口邊坡的沖蝕速率大體相等,因此假定潰口的深度和寬度以同樣的速率發(fā)展[18],則水流對(duì)壩體潰口兩側(cè)的直接沖刷形成的潰口寬度增量ΔB可表達(dá)為

        時(shí)間段Δt內(nèi)水庫(kù)水位變化量為

        式中:Qin——入庫(kù)流量;Sa——庫(kù)水位為H時(shí)的水庫(kù)面積。

        潰口受到水流的連續(xù)沖蝕發(fā)生垂向下切和橫向擴(kuò)展,邊坡也隨水流沖蝕變得越來(lái)越陡,當(dāng)垂向下切深度達(dá)到臨界深度時(shí),潰口邊坡發(fā)生間歇性失穩(wěn)坍塌(圖3)。臨界深度可采用極限平衡方法導(dǎo)出[18]:

        式中βk為潰口邊坡臨界坡角。

        3 數(shù)值計(jì)算方法

        可通過(guò)迭代計(jì)算來(lái)對(duì)上述潰壩數(shù)學(xué)模型進(jìn)行求解,以模擬堰塞壩漫頂破壞潰口發(fā)展過(guò)程,得出潰口流量過(guò)程線。為證明所建議模型的有效性和合理性,針對(duì)堰塞壩漫頂引起的潰壩問(wèn)題,進(jìn)行潰壩洪水過(guò)程計(jì)算,計(jì)算實(shí)例為唐家山堰塞壩,具體求解過(guò)程見(jiàn)圖4,圖中Hbm為潰口最終底高程,Bu為潰口最終頂寬,Bm為潰口最終底寬。

        圖3 土石壩漫頂破壞潰口發(fā)展示意圖Fig.3 Breach development of earth-rocKdam due to overtopping failure

        4 模型驗(yàn)證

        “5.12”汶川大地震堵江形成的堰塞壩中以唐家山堰塞壩規(guī)模最大,對(duì)下游威脅也最大。壩體位于北川縣城上游4.6km涪江支流通口河峽谷中,其堰塞壩體積約為2037萬(wàn)m3,按我國(guó)堰塞湖壩體等級(jí)劃分,唐家山堰塞壩屬于極高危險(xiǎn)、潰決損失極嚴(yán)重的Ⅰ級(jí)風(fēng)險(xiǎn)堰塞壩[19]。

        堰塞體的基本形狀為寬頂堰,上游坡度18°~22°,下游坡度35°~40°。順河向堰頂寬度為150m左右,堰塞體除表部少量為碎石土外,主體為基巖滑坡形成的巨石、孤塊石,其中下部巖性為弱風(fēng)化~微風(fēng)化,仍保持原地層層序關(guān)系,破裂巖體結(jié)構(gòu)密實(shí),強(qiáng)度高。壩體最大壩高124.4m,順河向底長(zhǎng)803.4m,頂長(zhǎng)300.0m左右,寬611.8m,體積約為2.037×107m3,壩底高程為669.5 m,上游坡度約為20°,下游平均坡比為1∶1.28[19]。經(jīng)過(guò)應(yīng)急搶險(xiǎn),截至2008年6月1日,堰塞壩頂開(kāi)挖出1條梯形泄流槽,斷面形狀為梯形,兩側(cè)邊坡為1∶1.5。泄流渠堰頂高程為740.4m,底寬8 m,深13m,總長(zhǎng)695m,上游平緩段縱坡為0.6%,下游陡坡段縱坡分別為24%和16%[19]。

        由于泄流槽經(jīng)2次爆破,6月10日才開(kāi)始泄流,所以選取6月10日0:00為起始計(jì)算時(shí)間。根據(jù)上述基本條件,計(jì)算時(shí)間步長(zhǎng)取Δt=0.02h,計(jì)算時(shí)間設(shè)定為20 h。其他計(jì)算參數(shù)如下:堰塞體高度為124.4m,堰塞體橫河向頂寬為611.8m,堰塞體順河向長(zhǎng)度為300.0m,下游平均坡比為1∶1.28,初始潰口寬度為8.0m,初始潰口深度為13m,上游入庫(kù)流量為0m3/s,堰塞體平均粒徑為0.03m,粒徑比d90/d30為30,堰塞體孔隙率為0.4,堰塞體材料平均密度為2.6×103kg/m3,內(nèi)摩擦角為45°,堰塞體下游平均坡角為23°,堰塞體黏聚力為60kPa,緊密系數(shù)為0.75,流量系數(shù)為0.5。

        計(jì)算得出唐家山堰塞壩潰壩洪峰流量出現(xiàn)在壩體潰決后11.8h,洪峰流量為6137.4m3/s(實(shí)測(cè)最大值為6500.0m3/s[19]),見(jiàn)圖5。圖6為潰口發(fā)展過(guò)程,潰口頂寬202.0m,底寬120.7 m,深度22.0m(實(shí)測(cè)潰口呈上寬下窄的“倒梯形”,其開(kāi)口寬145.0~225.0m,底寬100.0~145.0m,坡高10.0~60.0m[19])。由圖5~6可知,計(jì)算結(jié)果與實(shí)測(cè)資料接近,表明筆者所建議的堰塞壩潰壩數(shù)值模型和數(shù)值計(jì)算方法是合理的。

        圖4 堰塞壩漫頂潰決數(shù)值模擬流程Fig.4 Flow chart of numerical simulation of barrier dam breach due to overtopping failure

        圖5 潰口流量過(guò)程對(duì)比Fig.5 Comparison of calculated and observed discharge hydrographs at breach

        圖6 潰口發(fā)展過(guò)程Fig.6 Breach development

        5 結(jié) 語(yǔ)

        根據(jù)堰塞壩壩體土石料的寬級(jí)配特性,建立描述堰塞壩漫頂潰決過(guò)程中潰口發(fā)展規(guī)律與流量過(guò)程的數(shù)值模型。該模型引入與水流方向垂直的附加作用力來(lái)考慮粗顆粒對(duì)細(xì)顆粒的阻攔、遮蔽作用以及細(xì)顆粒對(duì)粗顆粒的包圍、填實(shí)作用?;诓煌馏w陡水槽沖蝕試驗(yàn)結(jié)果,針對(duì)寬級(jí)配堰塞壩材料在高強(qiáng)度水流侵蝕條件下的沖蝕特性,提出了一個(gè)能反映壩體材料顆粒級(jí)配、密實(shí)度、強(qiáng)度、壩體坡度、水流速度及摩阻流速對(duì)沖蝕量影響的土體沖蝕公式,公式中各物理量意義明確且易于測(cè)定。利用該模型對(duì)唐家山堰塞壩的潰壩過(guò)程進(jìn)行模擬計(jì)算,得出的潰口發(fā)展規(guī)律與潰壩洪水流量過(guò)程與實(shí)測(cè)結(jié)果接近,驗(yàn)證了該模型和計(jì)算方法的合理性,從而為正確評(píng)估堰塞壩潰決致災(zāi)后果、科學(xué)制定堰塞壩潰決防洪應(yīng)急預(yù)案、減輕因堰塞壩潰決造成的損失提供了一種新的有效手段。

        [1]任強(qiáng),陳生水,鐘啟明,等.堰塞壩的形成機(jī)理與潰決風(fēng)險(xiǎn)[J].水利水電科技進(jìn)展,2011,31(5):30-34.(REN Qiang,CHEN Shengshui,ZHONG Qiming,et al.Formation mechanism and breaching failure risKof barrier dams[J].Advances in Science and Technology of Water Resources,2011,31(5):30-34.(in Chinese))

        [2]鄧明楓,陳寧生,胡桂勝,等.松散及弱固結(jié)堰塞體潰壩形式與流量過(guò)程[J].水利水電科技進(jìn)展,2011,31(1):11-14.(DENGMingfeng,CHEN Ningsheng,HUGuisheng,et al.Outburst patterns and discharge process for weakly consolidated and loose barrier dams[J].Advances in Science and Technology of Water Resources,2011,31(1):11-14.(in Chinese))

        [3]COATA JE,SCHUSTER RL.The formation and failure of natural dams[J].Geological Society of America Bulletin,1988,100(1):1054-1068.

        [4]王光謙,鐘德鈺,張紅武,等.汶川地震唐家山堰塞湖泄流過(guò)程的數(shù)值模擬[J].科學(xué)通報(bào),2008,53(24):3127-3133.(WANG Guangqian,ZHONG Deyu,ZHANGHongwu,et al.Numerical simulation of dischargeprocess for Tangjiashanbarrier lake due to Wenchuan earthquake[J].Chinese Science Bulletin,2008,53(24):3127-3133.(in Chinese))

        [5]黃金池.堰塞壩漫頂潰口流量變化過(guò)程的數(shù)值模擬[J].水利學(xué)報(bào),2008,39(10):1235-1240.(HUANG Jinchi.Numerical modeling of flow through breach of landslide dams[J].Journal of Hydraulic Engineering,2008,39(10):1235-1240.(in Chinese))

        [6]錢(qián)寧,萬(wàn)兆惠.泥沙運(yùn)動(dòng)力學(xué)[M].北京:科學(xué)出版社,2003.

        [7]韓其為,何明民.泥沙起動(dòng)規(guī)律及起動(dòng)流速[M].北京:科學(xué)出版社,1999.

        [8]李榮,王迎春.非均勻沙起動(dòng)規(guī)律研究[J].泥沙研究,1999(2):27-32.(LI Rong,WANG Yingchun.Study on laws of threshold motion for non-uniform sediment[J].Journal of Sediment Research,1999(2):27-32.(in Chinese))

        [9]秦榮昱,王崇浩.河流推移質(zhì)運(yùn)動(dòng)理論及應(yīng)用[M].北京:中國(guó)鐵道出版社,1996.

        [10]王協(xié)康,敖汝莊.泥沙起動(dòng)條件及機(jī)理的非線性研究[J].長(zhǎng)江科學(xué)院院報(bào),1999,16(4):39-41.(WANG Xiekang,AO Ruzhuang.Nonlinear study on sediment incipient motion condition and its mechanism[J].Journal of Yangtze River Scientific Research Institute,1999,16(4):39-41.(in Chinese))

        [11]楊具瑞,方鐸,何文社,等.非均勻沙起動(dòng)規(guī)律研究[J].水利學(xué)報(bào),2002,33(10):82-86.(YANG Jurui,FANG Duo,HE Wenshe,et al.Study on laws of incipientmotionfor non-uniform sediment[J].Journal of Hydraulic Engineering,2002,33(10):82-86.(in Chinese))

        [12]LIU Ning,CHEN Zuyu,ZHANG Jianxin,et al.Draining the Tangjiashan barrier lake[J].Journal of Hydraulic Engineering,ASCE,2010,136(11):914-923.

        [13]胡卸文,黃潤(rùn)秋,施裕兵,等.唐家山滑坡堵江機(jī)制及堰塞壩潰壩模式分析[J].巖石力學(xué)與工程學(xué)報(bào),2009,28(1):181-189.(HUXiewen,HUANG Runqiu,SHI Yubing,et al.Analysis of blocking river mechanism of Tangjiashan landslide an dambreaking mode of its barrier dam[J].Chinese Journal of RocKMechanics and Engineering,2009,28(1):181-189.(in Chinese))

        [14]MEYER-PETER E,MULLER R.Formulas for bed load transport[C]//Proceedings of Second Meeting International Association for Hydraulic Research.Stockholm:IAHR,1948:39-64.

        [15]SMARTG M.Sediment transport formula for steep channels[J].Journal of Hydraulic Engineering,ASCE,1984,110(3):267-276.

        [16]范家驊,陳裕泰,金德春,等.懸移質(zhì)挾沙能力水槽試驗(yàn)研究[J].水利水運(yùn)工程學(xué)報(bào),2011(1):1-16.(FAN Jiahua,CHEN Yutai,JIN Dechun,et al.Experimental studies on carrying capacity of suspended load[J].Hydro-Science and Engineering,2011(1):1-16.(in Chinese))

        [17]李家星,趙振興.水力學(xué)[M].南京:河海大學(xué)出版社,2001.

        [18]SINGH V P.Dam breach modeling technology[M].Dordrecht,Netherland:Kluwer Academic Publisher,1996.

        [19]劉寧,張建新,林偉,等.汶川地震唐家山堰塞引流除險(xiǎn)工程及潰壩洪水演進(jìn)過(guò)程[J].中國(guó)科學(xué)E輯:技術(shù)科學(xué),2009,39(8):1359-1366.(LIU Ning,ZHANG Jianxin,LIN Wei,et al.Drainage reinforcement project and dam breaKflood process of Tangjiashanbarrier due to Wenchuan earthquake[J].Science in China E:Technological Sciences,2009,39(8):1359-1366.(in Chinese))

        东京热东京道日韩av| 亚洲日本va午夜在线电影| 亚洲av成人一区二区三区av| 免费国产黄片视频在线观看| 女同重口味一区二区在线| 亚洲男人天堂黄色av| 久久久久人妻精品一区蜜桃| 国产成年无码V片在线| 国产一区二区三区四区色| 欧美人与物videos另类xxxxx | 欧美亚洲国产精品久久久久| 毛片精品一区二区二区三区| 久久理论片午夜琪琪电影网| 日韩电影一区二区三区| 欧美精品v欧洲高清| 国产伦理一区二区久久精品| av无码电影一区二区三区| 欧美疯狂性xxxxxbbbbb| 国产成人亚洲精品2020| 青青草手机视频免费在线播放 | 国产一品二品三区在线观看| 国产无套粉嫩白浆在线观看| 亚洲的天堂av无码| 偷拍自拍一区二区三区| 狠色人妻丝袜中文字幕| 国产精品多p对白交换绿帽| 亚洲精品国产字幕久久vr| 亚洲在线一区二区三区| 色综合久久中文字幕综合网| 久久人人爽人人爽人人av| 精品人妻丰满久久久a| 久久综合亚洲鲁鲁五月天| 久久婷婷五月国产色综合| 亚洲av无码成人精品区天堂| 情色视频在线观看一区二区三区| 亚洲一区二区三区综合免费在线| 97久久草草超级碰碰碰| 日中文字幕在线| 在线视频精品少白免费观看| 国产精品伦理久久一区| 文字幕精品一区二区三区老狼|