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

洪水對水庫水溫分層結構的影響

2013-03-15 05:39:11美,逄勇,2,王華,鮑
水資源保護 2013年5期
關鍵詞:結構

陶 美,逄 勇,2,王 華,鮑 琨

(1.河海大學環(huán)境學院,江蘇 南京 210098;2.河海大學淺水湖泊綜合治理與資源開發(fā)教育部重點實驗室,江蘇 南京 210098;3.江蘇省交通規(guī)劃設計院股份有限公司,江蘇 南京 210005)

水庫建成蓄水后,受以年為周期的水溫、入流水溫、氣象規(guī)律性變化的影響,水庫在沿水深方向上呈現(xiàn)出有規(guī)律的水溫分層,水庫水溫分層可能直接導致庫區(qū)內的水質分層和生態(tài)分層,并且對農(nóng)田灌溉、工業(yè)供水、生活用水、下游河流的水質和生態(tài)平衡等方面都產(chǎn)生重要影響,因此對水庫水溫分層規(guī)律及其影響因子的研究具有重大意義[1-4]。

目前,我國關于水庫水溫分層結構的研究主要包括兩方面:①定性研究,采用徑流-庫容比數(shù)法(即參數(shù)α-β法)[5]判定水庫水溫結構,該法在山口巖水庫、劉家峽水庫的水溫研究中有著重要的應用[6-7];②定量研究,主要采用經(jīng)驗法和數(shù)值計算法分析水庫庫區(qū)的水溫分布。經(jīng)驗法主要有朱伯芳[8]、張大發(fā)[9]的水溫一維垂向經(jīng)驗公式法,預測水庫壩前垂向水溫具體分布,模型計算主要有垂向一維水溫模型,其代表有WRE模型、MIT模型以及“湖溫一號”模型[10],同時剖面二維的水庫水溫模型已有一些研究成果,并在實際工程中有一定的應用[3,11],三維耦合浮力流k-ε模型對漫灣水庫進行了應用[12],MIKE3水溫模塊在托什河干流梯級水庫水溫研究也得到了應用[13]。

在前人研究的基礎上,筆者定性分析水庫水溫分層結構的影響因子,定量研究各個影響因子的具體影響,以烏江梯級水庫中具有典型庫容和水溫分層結構的洪家渡水庫、東風水庫、烏江渡水庫為研究對象,研究不同歷時洪水量對不同水庫水溫分層結構的具體影響,同時根據(jù)模型計算和數(shù)據(jù)結果分析,完善β值的定義,提出具有應用價值的通用表,以期對全國其他流域梯級水庫水溫研究有借鑒作用。

1 水溫數(shù)學模型介紹

水庫水溫數(shù)值計算擬采用目前美國環(huán)保局推薦使用的環(huán)境流體動力學模型EFDC(environmental fluid dynamic code)三維水動力模型。

1.1 模型基本方程[14]

a. 動量方程。

?t(mHu)+?x(myHuu)+?y(mxHvu)+

?z(mwu)-(mf+v?xmy-u?ymx)Hv=

-myH?x(gξ+p)-my(?xH-z?xH)?zp+

?z(mH-1Av?zu)+Qu

(1)

?t(mHv)+?x(myHuv)+?y(mxHvv)+

?z(mwv)-(mf+v?xmy-u?ymx)Hu=

-mxH?y(gξ+p)-mx(?yh-z?xH)?zp+

?z(mH-1Av?zv)+Qv

(2)

?t(mξ)=-gH(ρ-ρ0)ρ0-1

(3)

b. 連續(xù)方程。

?t(mξ)+?x(myHu)+?y(mxHv)+?z(mw)=0

(4)

(5)

c. 狀態(tài)方程。

ρ=ρ(p,Sa,T)

(6)

d. 溫度輸移方程。

?t(mHT)+?x(myHuT)+?y(mxHvT)+

?z(mwT)=?z(mH-1Kv?zT)+QT

(7)

式中:u、v、w分別是邊界擬合正交曲線坐標x、y、z方向上的速度分量;t為時間;mx、my和分別為水平坐標變換因子;m為度量張量行列式的平方根,m=mxmy;Av為垂向紊動黏滯系數(shù);Kv為垂向紊動擴散系數(shù);f為科里奧利系數(shù);g為重力加速度;ρ為混合密度;ρ0是參考密度;H為總水深;h為未擾動的z坐標原點以下的水深;p為壓力;Sa為鹽度;T為溫度;ξ為自由勢能;Qu、Qv分別為動量在x、y方向的源匯項;QT是溫度的源匯項。

聯(lián)立各動量方程和連續(xù)方程即可解出u、v、w、ρ、p、Sa、T和ξ。

1.2 方程的離散求解

EFDC模型采用二階精度的空間有限差分格式求解控制方程,變量布置采用交錯網(wǎng)格。模型的時間積分采用具有二階精度的三層有限差分格式,采用內外模分裂方式將物理過程分解為內模和外模。外模求解采用半隱格式,利用預處理共軛梯度法同時求解二維水位場。內模求解垂向擴散項采用隱格式。動量方程內模主要是求解應力和速度的垂向分布。輸移方程中的平流項模型采用Smolarkiewicz多維正定對流輸移格式。

2 水溫數(shù)學模型的建立與率定

2.1 洪家渡水庫模型的建立

洪家渡水庫數(shù)值計算范圍為壩址處至回水段20 km,計算時段根據(jù)2006年水庫水溫實測資料選取,選擇水溫分層具有季節(jié)規(guī)律的4月、7月、10月、12月。

2.1.1初始條件的確定

a. 初始水位的確定。初始水位采用2006年洪家渡水庫壩前實測月均水位資料(表1)。

表1 2006年洪家渡水庫壩前實測月均水位和下泄水溫

b. 初始水溫條件的確定。根據(jù)水溫實測資料垂向分布情況,將水庫水溫沿水深方向分為20層,初始水溫根據(jù)2006年洪家渡水庫壩前斷面垂向水溫(表2)給定,取值為模擬月份前一個月的水溫實測值,用線性內插法設定每層的初始水溫。

2.1.2邊界條件的確定

a. 流量、水位邊界條件的確定。洪家渡水電站共設3個開邊界,入流邊界為凹水河、白甫河和木白河,出流邊界為壩前出流斷面。2006年1月1日—12月31日庫區(qū)實測日平均入、出庫流量見圖1和圖2。

表2 2006年洪家渡水庫壩前斷面垂向水溫分布

圖1 2006年洪家渡庫區(qū)實測入庫流量過程線

圖2 2006年洪家渡庫區(qū)實測出庫流量過程線

圖3 洪家渡水庫水溫率定結果

b. 水溫邊界條件的確定。主流的入庫水溫采用2006年洪家渡壩前斷面實測值(表2),支流水溫數(shù)據(jù)較小,水深較淺,河道特征明顯,垂向斷面水溫混合均勻,接近天然水溫。出流水溫根據(jù)2006年洪家渡水庫下泄水溫實測值給定(表1)。

c. 大氣邊界條件的確定。大氣邊界條件中大氣壓、空氣溫度、空氣濕度、降水量根據(jù)畢節(jié)地區(qū)氣象實測資料給定,蒸發(fā)量、太陽輻射量及云量由模型自行運算。

2.2 模型率定結果

洪家渡水庫2006年4月、7月、10月、12月壩前斷面垂向水溫率定結果見圖3。

由圖3得水溫計算值與實測值誤差為0.5~2.0℃,模型率定效果較好,可以用來模擬計算洪家渡庫區(qū)的水溫分布情況,對洪家渡水庫進行水溫率定的同時也對烏江流域的東風水庫和烏江渡水庫作了率定驗證。各水庫水溫參數(shù)率定結果見表3。

表3 各水庫水溫參數(shù)率定結果

3 不同來水量對烏江梯級水庫水溫分層結構的影響

3.1 水庫水溫判定方法

我國現(xiàn)行的水庫環(huán)境影響評價中普遍采用經(jīng)驗公式方法——徑流-庫容比數(shù)法來判定水庫的水溫結構。徑流-庫容比數(shù)法又稱為庫水替換次數(shù)法[5]。其判別指標為:

當α<10時,為穩(wěn)定分層型;當α>20時,為混合型;當10<α<20時,為過渡型。

對于分層型水庫,當β>1.0時,洪水對水溫結構有影響,為臨時混合型;當β<0.5時,洪水對水溫結構無影響;當0.5<β<1.0,洪水對水溫結構有一定影響,但未破壞水溫的分層結構。

3.2 烏江梯級水庫水溫結構

3.2.1洪家渡水庫

洪家渡水庫壩址多年平均流量為149 m3/s,電站水庫總庫容為49.47億m3,正常蓄水位為1 140 m,相應庫容為44.97億m3,調節(jié)庫容為33.60億m3,死水位為1 076 m,死庫容為11.36億m3,水庫回水長度為85 km。

目前我國普遍采用參數(shù)α-β法判別水庫水溫結構,洪家渡壩址多年平均徑流量為46.99億m3,水庫水溫結構的α值為0.95。根據(jù)判別式對α值進行分析,洪家渡水電站水庫屬于穩(wěn)定分層型水溫結構,水庫會產(chǎn)生水溫分層現(xiàn)象。

洪家渡庫區(qū)大洪水過程多為復峰,一般陡漲緩落,洪峰主要是由1 d暴雨造成,一次洪水以3 d洪量為主。根據(jù)實測洪水峰、量系列進行計算,洪家渡水電站設計洪水(五百年一遇)洪量為6 550 m3/s,3 d洪量為11.22億m3,水庫總庫容為49.47億m3,β值為0.227,洪水對洪家渡水庫水溫的分布結構沒有影響。5 d洪量14.94億m3,7 d洪量17.48億m3,β值分別為0.302、0.353,洪水對洪家渡水庫水溫的分布結構依然沒有影響。

3.2.2東風水庫

東風水電站壩址多年平均流量為345 m3/s。電站水庫總庫容為10.16億m3,調節(jié)庫容為4.91億m3,死庫容為3.73億m3。正常蓄水位為970 m,死水位為936 m。電站大壩設計洪水標準為百年一遇(P=1%),設計洪水位為974.08 m,校核標準為千年一遇(P=0.1%),校核洪水位為977.53 m。依據(jù)前述的α和β判別法,對東風水庫水溫結構判別如下:東風壩址多年平均徑流量108.80億m3,水庫總庫容為10.25億m3,α值為10.61。根據(jù)判別式對α值進行分析可知,東風水電站水庫屬于過渡型水溫結構,水庫會產(chǎn)生水溫分層現(xiàn)象,但不穩(wěn)定。

從歷年特大暴雨來看,3 d暴雨多集中在一天,因此流域內洪峰主要是由1 d暴雨造的,一次洪水過程主要集中在3 d洪水過程線內,3 d洪量為7.76億m3,水庫總庫容為10.25億m3,β=0.757,0.5≤β≤1.0,表明洪水對東風水庫水溫的分布結構的影響不明顯;5 d洪量為10.61億m3,7 d洪量為12.95億m3,β值分別為1.035、1.263,β>1.0,表明洪水對東風水庫水溫的分布結構有影響。

3.2.3烏江渡水庫

烏江渡壩址多年平均流量為483 m3/s,電站水庫總庫容為23.00億m3,調節(jié)庫容為9.28億m3,死庫容為12.12億m3,正常蓄水位為760 m。電站大壩設計洪水標準為五百年一遇(P=0.2%),設計洪水位760.30 m,校核標準為五千年一遇(P=0.02%),校核洪水位為762.80 m。烏江渡壩址多年平均徑流量為158億m3,水庫總庫容為21.40億m3,α值為7.383。根據(jù)判別式對α值進行分析,烏江渡水庫屬于過渡型水溫結構,水庫會產(chǎn)生水溫分層現(xiàn)象,但不穩(wěn)定。

烏江渡水庫3 d洪量為20.50億m3,5 d洪量為29.70億m3,7 d洪量為36.90億m3,相應的β值分別為0.96、1.39、1.72,由此可見大洪水對烏江渡水庫水溫結構影響非常顯著。

3.3 分析方法

分析文獻[3-4]中水庫壩前斷面水溫計算的方法,再結合水溫計算公式的基本原理,可得出汛期影響水庫水溫分層結構的因子主要有:一次洪量的大小、一次洪量歷時、水庫庫容等。文獻[3-4]中在計算β值時都未明確一次洪水量的歷時,這在判定洪水對水庫水溫分層結構的影響時不完全準確,對某一水庫,當一次洪水量相同時,計算得到的β值是相同的,即汛期洪水對該水庫水溫分層結構的影響是相同的,但當這一次洪水量經(jīng)歷的時間不一樣時,其相同的一次洪水量對水庫壩前水溫分層結構的影響程度是不同的。為了具體分析上述3個因子對水庫水溫分層結構的影響,筆者對烏江流域的洪家渡、東風、烏江渡3個水庫數(shù)值進行計算,總結分析3個因子對水庫水溫分層的不同影響,同時給出3個水庫在不同的β值條件下壩前斷面水溫分布情況,補充完善β值的定義,為烏江流域其他水庫在汛期受洪水影響后水庫水溫分層情況的研究作參考和借鑒。

根據(jù)洪家渡水庫3 d、5 d、7 d洪量值,采用內插法確定洪家渡水庫汛期受洪水影響程度時一次洪水量分別為7.64億m3、14.64億m3、20.68億m3,根據(jù)洪家渡水文站歷年洪水資料,利用同頻率放大法得到在3 d、5 d、7 d洪量為上述3組數(shù)值時的洪水過程線,以洪水過程線的流量值作為模型的邊界條件,利用EFDC模型計算在此洪水過程后水庫壩前斷面的水溫分布情況,從而說明一次洪水過程在不同歷時情況下水庫不同的壩前斷面水溫分布,根據(jù)庫表與庫底的最大溫差值判定洪水對水庫水溫結構的影響程度,在通過計算β值來說明洪水影響程度的基礎上,再通過計算數(shù)據(jù)更加具體地說明水庫在汛期經(jīng)歷不同歷時的一次洪水后水庫壩前斷面的水溫變化情況。

圖4 洪家渡水庫壩前斷面水溫分布

圖5 東風水庫壩前斷面水溫分布

圖6 烏江渡水庫壩前斷面水溫分布

單個水庫的計算說明洪水歷時、洪水量的大小對水庫水溫分層結構的影響,無法說明水庫水溫分層結構的另外一個影響因子——水庫庫容,因此需要計算不同水庫在汛期一次洪水量條件下的壩前斷面水溫分層情況,從不同水庫的計算結果可以分析得到在洪水歷時相同的條件下不同庫容水庫在汛期一次洪水情況下水溫分層結構。東風水庫和烏江渡水庫類比洪家渡水庫進行計算,東風水庫一次洪水量定為4.8億m3、7.3億m3、18.9億m3,烏江渡水庫一次洪水量定為14.2億m3、29.5億m3、41.2億m3。同樣利用EFDC模型進行計算,得到不同洪水過程后水庫壩前斷面水溫分布情況。

3.4 結果及分析

利用EFDC模型計算洪家渡水庫、東風水庫、烏江渡水庫汛期經(jīng)歷不同歷時的一次洪水后壩前斷面水溫分布情況,每個水庫計算相同的一次洪水量不同歷時過程后壩前水溫分布情況,計算結果見圖4~圖6。

從圖4~圖6可以看出,對于同一個水庫,當一次洪水量相同時,洪水歷時越長,洪水對水庫水溫分層結構影響越小,以烏江渡水庫為例,在經(jīng)歷一次洪水量為41.2億m3時,若歷時3 d,水庫庫表與庫底的溫差為1.7℃;歷時7 d,溫差為3.0℃,可以看出一次洪水量歷時越長對水庫結構影響越小,該規(guī)律對其他兩水庫也同樣適用。從圖4~圖6可以看出,對于同一個水庫,當一次洪水量相同時,洪水歷時越長,洪水對水庫水溫分層結構影響越小,以烏江渡水庫為例,在經(jīng)歷一次洪水量為41.2億m3時,若歷時3 d,水庫庫表與庫底的溫差為1.7℃;歷時7 d,水庫庫表與庫氏的溫差為3.0℃,可以看出一次洪水量歷時越長對水庫水溫分層結構影響越小,該規(guī)律對其他兩水庫也同樣適用。當一次洪水量相同,洪水歷時也相同,庫容越大的水庫洪水過程對該水庫的壩前斷面水溫分層結構影響越小。以洪家渡(庫容45億m3)水庫和烏江渡水庫(庫容21億m3)為例,當洪水歷時3 d,一次洪水量為14億m3、26億m3時,洪家渡水庫庫表與庫底的溫差分別為4.8℃、2.7℃,烏江渡水庫庫表與庫底的溫差分別為3.2℃、2.2℃,即一次洪量增加12億m3時洪家渡水庫庫表與庫底的溫差變化了2.1℃,烏洪渡水庫庫表與庫底的溫差變化了1.1℃,說明洪水對洪家渡水庫水溫分層結構影響較小。

為了完善β值的定義,給β值補充洪水歷時的概念,以烏江流域3座水庫不同洪水過程的計算結果為基礎,總結得到不同庫容的水庫在不同β值時水庫壩前斷面庫表水溫與不同水深處的水溫的差值(表4~表6),即通用表,以溫度差來說明壩前斷面的水溫情況,通過β值的量化來說明汛期一次洪水對水庫水溫結構的影響。

表4 3 d洪量不同β值時不同庫容的水庫庫表與不同水深處的溫差值

表5 5 d洪量不同β值時不同庫容的水庫庫表與不同水深處的溫差值

表6 7 d洪量不同β值時不同庫容的水庫庫表與不同水深處的溫差值

從表4~表6可以看出,在3 d洪量條件下,對于庫容為10億m3的水庫,當β=0.15~1.95時,庫表與庫底相應的溫差為4.99~0.65℃,對于庫容為45億m3的水庫,當β值為0.15~1.95時,庫表與庫底相應的溫差為4.50~0.60℃。可以看出,對于同一個水庫,β值越大,庫表與庫底溫差越小,水溫分層結構越不明顯。對于庫容為10億m3的水庫,當β=0.15時,在3 d洪量條件下和7 d洪量條件下,庫表與庫底相應的溫差分別為4.99℃、7.02℃,由此可以說明,對于同一水庫,洪水歷時越長,洪水對水庫水溫分層結構影響越小,水溫分層結構越明顯,庫表與庫底溫差越大。

分析水庫數(shù)值計算結果結合其他水庫水溫歷史數(shù)據(jù),對β值中的一次洪水量補充定義,一次洪水量以歷時7 d作為參考,對于分層型水庫(包括穩(wěn)定分層型和過渡型水庫),當β>1.0時,洪水對水溫結構有影響,為臨時混合型;當β<0.5時,洪水對水溫結構無影響;當0.5<β<1.0,洪水對水溫結構有一定影響,但未破壞水溫的分層結構,仍為分層型水庫。對于其他歷時的一次洪水量β值的定義區(qū)間也適用,各個區(qū)間所說的“有影響”、“有一定影響”、“臨時混合”的程度不一樣,具體即表現(xiàn)為不同歷時的一次洪水量情況下,水庫庫表與庫底的溫差不相同。以東風水庫的模型數(shù)值計算結果闡述β的補充定義,當β=0.15<0.5時,洪量歷時3 d時,庫表與庫底的溫差為5.28℃,歷時7 d時,庫表與庫底的溫差為7.02℃,這兩個溫差都可稱為對東風水庫水溫分層結構無影響,水庫水溫結構仍然為穩(wěn)定分層型;當β=0.75>0.5時,洪量歷時3 d時,庫表與庫底的溫差為2.36℃,歷時7 d時,庫表與庫底的溫差為3.21℃,這兩個溫差都可稱為對東風水庫水溫分層結構有一定的影響,歷時3 d的洪水對水溫分層結構影響較大;當β=1.95>1時,洪量歷時3 d條件下,庫表與庫底的溫差為0.75℃,歷時7 d條件下,庫表與庫底的溫差為1.12℃,水庫水溫結構都為臨時混合型,歷時3 d混合程度較高。

4 結 論

a. 當一次洪水歷時相同,若洪水量相同,水庫的庫容越大,洪水對水溫的分層結構影響越小;若水庫庫容相同,洪水量越小,洪水對水溫的分層結構影響越小。當一次洪水量相同,若水庫庫容相同,洪水歷時越長,洪水對水溫的分層結構影響越小。

b. 通過分析水庫數(shù)值計算結果和參考相關文獻,對參數(shù)α-β判定法中的β值進行了時間上的定義,深入闡述了β值的具體含義,制定了不同洪量對烏江流域水庫水溫分層結構影響的通用表。該水溫研究方法對全國其他流域水電梯級開發(fā)水溫研究可提供參考和借鑒。

[ 1 ] 吳莉莉,王惠民,吳時強.水庫的水溫分層及其改善措施[J].水電站設計,2007,23(3):97-100.(WU Lili,WANG Huimin, WU Shiqiang. Water temperature hierarchy of the reservoir and its improving measures [J].Design of Hydroelectric Power Station,2007,23(3):97-100.(in Chinese))

[ 2 ] PETTS G E.蓄水河流對環(huán)境的影響[M]. 王兆印,曾慶華,呂秀貞,等,譯.北京:中國環(huán)境科學出版社,1988:25-27.

[ 3 ] 戴群英.水庫庫區(qū)及下游河道水溫預測研究[D].南京:河海大學,2006.

[ 4 ] 鄧云.大型深水庫的水溫預測研究[D].南京:河海大學,2003.

[ 5 ] SL 278—2002 水利水電工程水文計算規(guī)范[S].

[ 6 ] 詹曉群,陳建,胡建軍.山口巖水庫水溫計算及其對下游河道水溫影響分析[J].水資源保護,2005,21(1):29-35.(ZHAN Xiaoqun,CHEN Jian,HU Jianjun. Mountain rock reservoir water temperature calculation and the analysis of influence on the downstream river water temperature[J].Water Resources Protection,2005,21(1):29-35. (in Chinese))

[ 7 ] 周孝德,張華麗.劉家峽水庫水溫影響回顧評價研究[D].西安:西安理工大學,2008.

[ 8 ] 張大發(fā).水庫水溫分析及估算[J].水文,1984(1):19-27.(ZHANG Dafa. Analysis and estimation of water temperature of reservoir[J].Journal of China Hydrology,1984(1):19-27. (in Chinese))

[ 9 ] 朱伯芳.水庫水溫度估算[J].水利學報,1985(2):12-21.(ZHU Bofang. Estimation of water temperature of reservoir[J]. Journal of Hydraulic Engineering,1985(2):12-21.(in Chinese))

[10] 雒文生,宋星原.水環(huán)境分析及預測[M].武漢:武漢水利電力大學出版社,2000:27-29.

[11] 江春波,馬強,付清潭.二維溫度分層流的數(shù)值模擬[J].水力發(fā)電學報,2003,29(2):24-26.(JIANG Chunbo,MA Qiang,FU Qingtan. Numerical simulation of two-dimensional temperature stratified flow[J]. Journal of Hydroelectric Engineering,2003,29(2):24-26. (in Chinese))

[12] 陸俊卿,張小峰,強繼紅,等.水庫水溫數(shù)學模型及其應用[J].水力發(fā)電學報,2008,27(5):124-129.(LU Junqing,ZHANG Xiaofeng, QIANG Jihong,et al. Reservoir water temperature mathematical model and its application[J]. Journal of Hydroelectric Engineering,2008,27(5):124-129. (in Chinese))

[13] 周孝德,孫東遷.托什干河上游梯級水庫水溫預測研究[D].西安:西安理工大學,2008.

[14] 劉夏明,李俊清,豆小敏,等.EFDC模型在河口水環(huán)境模擬中的應用及進展[J].環(huán)境科學與技術,2011,34(6G):136-140.(LIU Xiaming, LI Junqing, DOU Xiaomin,et al. The application and advance of environmental fluid dynamics code (EFDC) in estuarine water environment[J]. Environmental Science & Technology,2011,34(6G):136-140. (in Chinese))

猜你喜歡
結構
DNA結構的發(fā)現(xiàn)
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
循環(huán)結構謹防“死循環(huán)”
論《日出》的結構
縱向結構
縱向結構
我國社會結構的重建
人間(2015年21期)2015-03-11 15:23:21
創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 亚洲成肉网| 国产美女在线观看| 99精品在线视频观看| 色老二精品视频在线观看| 高清欧美性猛交XXXX黑人猛交| 国产欧美网站| 最近最新中文字幕在线第一页| 国产第一页屁屁影院| 亚洲妓女综合网995久久| 91丨九色丨首页在线播放 | www.精品视频| 日本国产一区在线观看| 免费Aⅴ片在线观看蜜芽Tⅴ| av一区二区三区在线观看 | 天天做天天爱天天爽综合区| 亚洲欧州色色免费AV| 国产免费黄| 免费看a级毛片| 久久久久亚洲精品成人网| 中文无码影院| 四虎AV麻豆| 日本午夜精品一本在线观看| 草草线在成年免费视频2| 亚洲无线一二三四区男男| 午夜毛片福利| 欧美成人免费一区在线播放| 国产成人在线无码免费视频| 一本大道香蕉中文日本不卡高清二区| 在线一级毛片| 亚洲日本精品一区二区| 欧美特级AAAAAA视频免费观看| 白浆视频在线观看| 国产91高跟丝袜| 国产高清精品在线91| 波多野结衣一区二区三区AV| 国产福利免费视频| 毛片手机在线看| 亚洲欧洲自拍拍偷午夜色| 在线播放精品一区二区啪视频| 激情無極限的亚洲一区免费| 狠狠色婷婷丁香综合久久韩国| 精品综合久久久久久97超人| 国产精品视频猛进猛出| 国产麻豆精品手机在线观看| 天堂在线亚洲| 久久国产精品77777| 青青热久免费精品视频6| 青草精品视频| 日韩精品一区二区三区中文无码 | 日韩在线影院| 91久久青青草原精品国产| 日韩在线第三页| 久久精品国产999大香线焦| 四虎综合网| 老司机久久精品视频| 国产欧美专区在线观看| 国产主播喷水| a级毛片毛片免费观看久潮| 宅男噜噜噜66国产在线观看| 国产对白刺激真实精品91| 一本大道香蕉久中文在线播放| 日韩二区三区无| 一级做a爰片久久毛片毛片| 成人国产精品2021| 波多野结衣一区二区三区四区 | 精品自窥自偷在线看| 在线观看国产精品一区| 欧美全免费aaaaaa特黄在线| 中文成人无码国产亚洲| 在线视频精品一区| 国产手机在线小视频免费观看| 人妻熟妇日韩AV在线播放| 国产成人免费高清AⅤ| 日本三级欧美三级| 久996视频精品免费观看| 亚洲成人一区二区三区| 特黄日韩免费一区二区三区| 激情六月丁香婷婷四房播| 国产精彩视频在线观看| 国产91丝袜在线播放动漫| 三上悠亚精品二区在线观看| 久久久成年黄色视频|