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

基于輸水損失動態變化特性的渠系優化配水模型研究

2021-08-20 07:13:42廖相成胡鐵松
水利學報 2021年7期
關鍵詞:模型

廖相成,胡鐵松

(武漢大學水資源與水電工程科學國家重點實驗室,湖北武漢 430072)

1 研究背景

渠系輸配水計劃編制的主要任務是在渠系物理特性固定條件下,依據未來農業氣象條件和灌溉可利用水量變化合理確定灌區各級渠道引水量大小、持續時間和引水頻次。因此,輸配水直接關系到灌溉服務質量,輸配水多了會導致灌區水資源浪費、地下水位上升和鹽堿化等問題;反之,輸配水少了則導致干旱缺水和作物減產。我國大型灌區骨干渠系水利用系數在2015年才達到0.597[1],大量的水資源在渠系輸配水過程中因滲漏損失浪費掉了。從這個角度講,如何進行渠系優化配水以減少渠系滲漏損失,對于我國節水型社會建設以及最嚴格水資源管理具有重要的意義。

目前渠系優化配水研究主要聚焦在渠系優化配水模型構建和模型求解算法上[2],并呈現出以下特點:對田間渠系配水的研究多于對骨干渠系配水的研究,對具體某一級渠系配水的研究多于對整體灌溉渠系配水的研究[3],對求解算法的研究多于對模型構建的研究[4](如模型目標中滲漏量最小是否合適與全面),尤其是等流量和不等流量條件下采用0-1規劃進行田間渠系輪灌優化配水編組的研究較多[5-10]。滲漏損失最小化是渠系優化配水在模型構建方面研究的主要目標之一[6,8,11],但目前各種優化配水模型中對滲漏損失特征的描述和準確量化研究還存在明顯不足。渠系優化配水是根據不同渠道控制區域內的種植結構、作物需水量和土壤質地的時空差異性來合理調配各級渠道適宜引水過程以實現滲漏損失最小化的目標。從這個意義上講,準確量化輸配水過程滲漏損失是能否實現這一目標的關鍵,也是優化配水的主要挑戰之一。部分學者[12]采用簡單的打折法確定輸配水滲漏損失,Mishra等[13]在印度普爾納哈拉分干渠輸配水研究中假設農渠和斗渠分別按照20%和15%計算;部分學者如Delgoda 等[14]、Liu 等[15]、Guo 等[16]等假設渠道中的水流始終處于穩態,將渠道每公里輸水損失看成一個常數,或者假定渠床土壤透水系數和指數具有時不變特征,不隨配水時間進程變化,采用如考斯加科夫(Kostiakov)公式等經驗公式法計算滲漏損失[17-20]。然而在實際的渠道輸配水過程中,土壤透水能力以及流量、流速等水力因素都是顯著改變的,這導致滲漏損失具有明顯的隨時間、隨流程動態變化的特征。忽視這些特征將導致滲漏損失被嚴重高估或低估,與實測結果相比,Shah等[21]和Zhang等[22]發現考斯加科夫公式的計算值高估了1.5倍和2.5倍,Akkuzu[23]發現莫里茲(Moritz)公式和戴維斯-威爾遜(Davis-Wilson)公式估算的滲漏損失遠低于實測值。

土壤水分入滲是水分在土壤水勢梯度的作用下通過土壤孔隙進入土壤并在其中運動的過程。Mishra 等[13]認為除受渠床土壤性質、渠道斷面形式與襯砌狀況、渠道水深等水力特性以及地下水埋深等多種因素的影響之外,渠床土壤濕潤狀況是渠床土壤入滲速度的重要制約因素;另一方面,渠床土壤濕潤狀況也受到渠道滲漏、地下水補給等因素的影響,是一個動態變化的量,從水量平衡來看,渠道輸配水過程中渠床土壤當前階段的濕潤狀況只與前一階段的濕潤狀況以及前一階段的地下水埋深、流量等因素有關,而當前階段的濕潤狀況一旦確定,渠床土壤濕潤狀況未來階段的發展不再受當前階段以前各階段狀況的影響,即渠床土壤濕潤狀況具有無后效性。渠床土壤濕潤狀況的動態變化及無后效性特征對渠床土壤入滲具有重要影響,但目前渠系優化配水滲漏計算經驗公式法并沒有考慮這些特征,而是將相關的模型參數處理為常數,故而產生了較大的計算誤差。因此,建立能夠描述渠床濕潤狀況動態變化及無后效性特征的渠系優化配水模型是減少渠系滲漏損失的關鍵。

本文著力描述渠系滲漏損失動態衰減特征,通過對渠床濕潤狀況無后效性的描述,提出了渠床土壤透水系數和指數動態計算方法,在此基礎上建立了一種新的考慮渠系滲漏損失動態特征的渠系優化配水模型,并在河套灌區骨干渠系優化配水中取得了良好的效果。

2 模型與方法

本文建立了一種新的考慮渠系滲漏損失動態特征的渠系優化配水模型,它包括兩個模塊,一個是以輸水損失最小化為目標的渠系優化配水模型,一個是基于前期影響含水量的渠道輸水損失動態計算模型。

2.1 基于前期影響含水量的渠道輸水損失動態計算模型(DCLM)要準確計算渠道輸水的滲漏損失,僅把水力參數作為變量是不夠的,還必須考慮其他因素的變化。渠床土壤干濕程度直接關系其透水性能,對輸配水滲漏損失有較大影響,渠床土壤濕潤度越高,水分入滲鋒面土水勢越高,水分入滲鋒面與渠床水勢梯度就越小,因此土壤入滲能力越低[24-27]。渠系輸配水過程初期或者是輪灌渠道,渠床土壤干燥,此時土壤入滲處于自由滲漏階段,渠床土壤入滲速度較大;反之,若渠道行水時間較長,處于輸配水過程中后期,渠床土壤已被充分濕潤,渠床滲漏可能處于穩滲階段,渠床土壤入滲速度較小。

渠床土壤含水率是準確表征其濕潤程度的重要指標,但存在直接測量成本高和時空變異性大的缺陷,難以在灌區這樣大尺度輸配水管理中應用。借鑒水文學前期影響雨量的概念,本文探索性提出渠床土壤前期影響含水量和相對濕潤度作為間接衡量渠床土壤干濕程度的指標,描述渠床土壤干濕程度在渠道輸水過程中的變化,并將其與常用的計算渠道輸水損失的經驗公式參數結合,從而形成了既形式簡便、又能考慮渠床土壤透水性變化且滿足渠系輸配水管理的渠道輸水滲漏損失動態計算模型。

渠床土壤前期影響含水量的計算需考慮渠道輸水的下滲水量補給、土壤水的蒸散發消耗及土壤蓄水容量上限,其計算公式為:

式中:Wa,t+1、Wa,t分別為第t+1 時段、第t時段的前期影響含水量,m3;Vt為第t個時段渠道的輸水量,m3;K為渠床土壤蓄水的消退系數;Wam為渠床土壤蓄水容量,m3。用式(1)計算Wa時,可取渠道停水較長時間后的Wa等于零,或取渠道連續大流量運行較長時間后的Wa等于Wam,由此向后逐時段推算,便可求得渠道輸水過程中各個時段的Wa。渠床土壤的干濕程度可用相對濕潤度表示:

式中:θt為第t個時段渠床土壤的相對濕潤度,反映了渠床土壤的透水性特征。考斯加科夫公式中渠床土壤透水系數A和指數m隨渠床土壤的相對濕潤度θt變化而動態變化,本文將其視為與渠床土壤相對濕潤度動態有關的動態變量A(θt)和m(θt),鑒于考斯加科夫公式在我國應用的廣泛性,提出下述基于前期影響含水量的渠道輸水損失動態計算方法:

式中:σ(t)為每公里渠道輸水損失系數,即單位渠長的流量損失率,1/km;Qn為渠道凈流量,m3/s。從式(3)可以看出,單位渠長流量損失率σ不僅與渠道凈流量Qn有關系,還與渠床土壤動態透水系數A(θt)和動態指數m(θt)密切相關,而A(θt)和m(θt)是隨渠床土壤相對濕潤度θ的變化而變化的,是θ的函數。A(θt)和m(θt)的具體函數形式依據實際問題確定,本文簡化為線性關系:

式中a、b、c、d為常數。

渠道輸水損失動態計算模型中的參數K,Wam以及函數關系A(θt)、m(θt)需要結合實測的渠道輸水資料進行率定,率定過程如下:首先根據渠道輸水的流量數據,采用水量平衡法計算渠道輸水過程中各個時段的水量損失ΔW和單位渠長流量損失率σ,計算公式見式(6)、式(7):

式中:ΔT為輸水時段時間長度,s;Q人和Q出分別為流入渠道和流出渠道的流量,m3/s;L為渠道長度,km;然后假定不同的K,Wam組合,用式(1)計算渠道輸水過程中各個時段對應的前期影響含水量以及渠床土壤相對濕潤度θ;再以θ為參數,點繪某一θ下的單位渠長流量損失率與渠道凈流量的關系曲線,即得到渠道的σ~θ~Qn關系圖;最后用式(3)對關系圖中的曲線進行擬合,優選出使擬合優度最大的K、Wam,同時得到不同θ下的A與m值,采用式(4)和式(5)擬合即得到函數關系A(θt)、m(θt)。

渠道輸水損失動態計算模型率定完成后,利用該模型計算渠道輸水的水量損失,與采用水量平衡法計算的水量損失結果進行對比,驗證模型的可靠性。模型檢驗的評價指標選擇平均相對誤差(MRE)和相關系數(R):

式中:n為樣本數;Pi為動態計算模型第i個模擬值;Oi為水量平衡法第i個計算值;Pave和Oave分別為動態計算模型模擬和水量平衡法計算的平均值。

2.2 渠系優化配水模型灌溉渠系輸配水計劃的編制一般要經歷自上而下——自下而上——再自上而下3個編制過程[16],即上級渠道的管理者發布整個灌溉區域的用水計劃,下級灌溉子區用水戶以此為基礎各自制定子區用水計劃,最后再由上級管理者統一協調,確定盡量滿足各子區的需水要求等多個目標的計劃。計劃編制建模過程中通常將其概化為一個多級多水源的多目標決策問題。多級指灌溉渠系是從水源取水、輸送、分配到田間的各級渠道網絡,如總干渠、干渠、分干渠、支渠等骨干渠道和斗渠、農渠、毛渠等田間渠道,多水源包括河流、水庫等多種水源,目標包括滲漏損失最小、產量最大等。為簡化計算,本模型僅以單水源雙目標的兩級骨干渠道組成的渠系優化配水問題建立模型。

假設上級渠道共有M個渠段U1、U2、…、UM,連接有N條下級渠道D1、D2、…、DM,在渠段末端分水給各下級渠道再由下級渠道將水量分配到各灌溉子區C1、C2、…、CZ,配水時段總數為T,如圖1所示。

圖1 灌區渠系示意(U1—UM,上級渠道渠段編號;D1—DN,下級渠道編號;C1—CZ,灌溉子區編號)

2.2.1 目標函數 模型以N條下級渠道在T個時段的渠首配水流量為決策變量,將目標函數設置為下級用水戶缺水量和上級渠道輸水水量損失最小,即:

式中:Wl為上級渠道輸水的水量損失,m3;Ws為下級用水戶各時段的缺水總量,m3;Qjtd為下級渠道Dj渠首第t時段的輸水流量,m3/s;Δt為時段長度,s;Wjt為第j條下級渠道第t時段的需水量,m3;qiut為上級渠道Ui第t時段的損失流量,m3/s。

渠道流量損失分別采用考斯加科夫經驗公式和上節中的動態計算模型計算,以對比分析兩種優化配水結果的差異。采用考斯加科夫經驗公式計算時:

式中:Aiu、miu分別為上級渠道Ui的渠床土壤透水系數和指數;Qitu為上級渠道Ui末端第t時段的輸水流量,m3/s;Liu為上級渠道Ui的長度,km。

采用動態計算模型計算時:

2.2.2 約束條件

(1)渠道輸水能力約束:為了保證渠道安全穩定運行,渠道輸水流量不得超過設計流量。

式中:Qius、Qjds分別為上級渠道Ui、下級渠道Dj的設計流量,m3/s。

(2)灌溉可引水量約束:受來水條件限制,上級渠道渠首引水流量不得超過其最大可引流量。

式中Qtumax為上級渠道第t時段的最大可引流量,m3/s。

(3)渠段水量平衡約束:上級渠道第i渠段的末端流量應等于該渠段末端所連的各下級渠道渠首流量與上級渠道第i+1渠段渠首流量之和。

(4)非負約束:各決策變量的技術性約束。

(5)公平性約束:受到上級渠道設計流量和來水流量的限制,各下級渠道灌溉子區可能出現缺水情況。為了保證配水的公平性,需將各灌溉子區的缺水比例差值控制在合理的范圍內。

式中:Rj為第j條下級渠道所在的灌溉子區的缺水比例;R max、R min分別為各灌溉子區中的最大缺水比例、最小缺水比例;r為允許的缺水比例差值。

3 實例研究

3.1 研究區概況內蒙古自治區河套灌區的地理區域為106°20′—109°19′E、40°19′—41°18′N之間,總土地面積119 hm2,是全國特大型灌區之一,灌區自西向東分為烏蘭布和、解放閘、永濟、義長和烏拉特五大灌域。灌區氣候干旱,屬于沒有灌溉便沒有農業的地區[28]。總干渠是河套灌區一首制引黃灌溉的輸水大動脈,位于河套平原南緣,共連接有14條干渠、9條分干渠和2條支渠,其輸水運行狀況對灌區渠系水的利用效率有著重大影響。為簡化起見,按照所屬灌域及其引水口位置,本文將總干渠上的干渠、分干渠、支渠和直口渠等多級渠道概化為7條干渠,如表1和圖2所示,只研究總干渠-干渠組成的兩級灌溉渠系優化配水問題。

圖2 河套灌區總干渠系示意

總干渠各渠段以及各下級干渠的主要技術參數如表1所示。

表1 總干渠系渠道參數

河套灌區目前渠系輸配水總干渠以續灌為主,而永濟等灌域采取以輪灌為主,大輪灌、小續灌、輪續結合的方式。受地域差異大、“三秋”(秋收、秋翻、秋澆)矛盾突出等問題影響,難以做到科學合理配水和集中統一調度,存在渠系滿負荷運行時間短、半渠水運行和開關口時間不一等配水流量合理控制問題,輸配水效率尚待提升。從圖3顯示的2002—2019年總干渠輸水流量概率分布來看,總干渠整體輸配水流量偏小,小流量和半渠水運行的比例超過了60%,不利于渠道水利用效率的提高。

圖3 總干渠2002—2019年輸水流量的概率分布(流量數據為總干渠渠首2002—2019年的逐日引水流量,m3/s;設計流量為565m3/s;概率密度、累積概率分別采用MATLAB中的ksdensity函數、ecdf函數計算)

圖4為總干渠2002—2019年的典型年開關閘時間變化,可以看出,總干渠開關閘時間和引水進度年際變化大。在開閘階段,部分年份總干渠開閘后流量增加迅速,渠道很快達到大流量運行狀態如2009年;相反,部分年份總干渠開閘后流量增加緩慢,渠道較長時間處于小流量運行甚至停水狀態,如2003年。這種年際變化主要是由不同年份作物種植結構、氣象條件、渠道可引水量等因素的差異導致的,但從渠道輸配水效率來看,總干渠2003年這種輸水方式易造成較大的水量損失,是需要改進的。

圖4 總干渠2002—2019年典型年份的開關閘時間變化分析

3.2 動態滲漏損失模型參數率定與驗證利用河套灌區總干渠及其下級渠道2002—2006年的實測輸水資料對模型參數進行率定。圖5為總干渠各渠段的σ~θ~Qn關系,表2給出了總干渠各渠段渠道輸水損失動態計算模型率定后的參數。由圖5可以看出,渠床土壤相對濕潤度相同時,單位渠長流量損失率與流量呈反比關系;θ越小,對應的σ~Qn曲線位置越高,表明渠道流量相同時,渠床土壤相對濕潤度越大,單位渠長流量損失率越小,這些規律都是與實踐經驗相符的。另外,與考斯加科夫公式的擬合結果相比,渠道輸水損失動態計算模型的擬合優度普遍更大,表明該模型具有更好的擬合效果。

表2 總干渠各渠段渠道輸水損失動態計算模型參數

圖5 渠道的σ~θ~Qn 關系(θ,渠床土壤相對濕潤度;同一顏色的點代表某一θ下單位渠長流量損失率與流量的數據分布;無k的曲線為相應顏色數據散點用動態計算模型擬合的趨勢線,并給出了擬合表達式和擬合優度;標k的曲線是用考斯加科夫公式對所有數據點整體擬合而成)

采用河套灌區總干渠及其下級渠道2007—2009年實測輸水資料對渠道輸水損失動態計算模型和考斯加科夫公式進行驗證,結果如表4所示。可以看出,動態計算模型驗證的相關系數大于0.7,平均相對誤差小于20%,是較為可靠的;同時,與考斯加科夫公式相比,動態計算模型驗證計算的相關系數更大,平均相對誤差更小,表明動態計算模型在計算渠道輸水損失時更為準確。

表4 模型驗證結果

3.3 不同輸配水方案滲漏損失對比分析利用前文建立的渠系優化配水模型,計算求解研究區域2012年夏灌期間(4—6月)總干渠-干渠兩級灌溉渠系優化配水的優化配水過程,以每5日為一個配水時段,共19個時段。根據灌區氣象資料、作物種植情況計算7條概化的干渠需水過程,依據劉家峽水庫下泄流量資料制定灌區可引水量過程。

將基于渠道輸水損失動態計算模型(DCLM)的渠道優化配水結果與基于考斯加科夫公式(Kostia?kov)的渠道配水結果進行對比,兩種渠系優化配水方式下總干渠各渠段的水量損失情況見表5。

由表5可知,在總干渠夏灌期間的渠首引水量皆為15.552億m3的條件下,與基于Kostiakov 公式的渠系優化配水模型相比,基于DCLM 模型配水方案顯著地減少了各個渠段的滲漏損失,滲漏水量損失減少達1560萬m3,渠道水利用效率由0.919提高到0.929,特別是總干渠三閘以下渠段。

表5 總干渠各渠段輸水及損失水量統計

兩種渠系優化配水方式下總干渠各渠段的流量損失率變化狀況如圖6所示。可以看出,隨著輸水時間的延長,兩種配水方式下各渠段的流量損失率均呈現出由初始較大值逐漸減小,最后趨于穩定的規律。然而,在4月上旬至5月上旬的渠道開閘階段,兩種配水方式對應的渠道流量損失率變化過程有顯著的差異。與基于Kostiakov 公式的配水方式相比,基于DCLM 模型的配水方式中,各渠段單位渠長流量損失率的初始值更小,且更快趨于穩定水平:基于Kostiakov 公式的配水方式中,U1—U4單位渠長流量損失率的初始值分別為0.0030、0.0047、0.0130、0.0165(1/km),基于DCLM 模型的配水方式中相應初始值分別為0.0023、0.0034、0.0083、0.0139(1/km),減小幅度為16%~36%;基于Kostiakov 公式的配水方式中,U1—U4單位渠長流量損失率減小到穩定水平的日期分別為5月6日(0.00039)、5月11日(0.00065)、5月6日(0.00143)、5月16日(0.00349),基于DCLM 模型的配水方式中相應日期分別為5月1日(0.00036)、5月1日(0.00067)、5月1日(0.00143)、5月10日(0.00326),相應提前5~10日。5月上旬后,兩種配水方式下各渠段的流量損失率均減小到較低水平,差異不大。

圖6 總干渠單位渠長流量損失率變化對比

由表5和圖6還可以看出,總干渠各渠段的渠道水利用效率差別較大,在兩種優化配水模型中均呈現出從U1至U4依次降低的規律。

3.4 不同輸配水方案配水進度對比分析將基于渠道輸水損失動態計算模型(DCLM)的渠道配水過程與基于考斯加科夫公式(Kostiakov)的渠道配水過程進行對比,總干渠及7條概化干渠輸配水流量過程分別如圖7和圖8所示。

圖7 兩種配水方式對應的總干渠流量過程對比

圖8 兩種配水方式對應的干渠流量過程對比

由圖7、圖8可以看出,兩種配水方式對應的渠道流量過程之間存在明顯的差異。與基于Kostia?kov 公式的渠道輸配水流量過程相比,基于DCLM 模型的優化配水在4月10日至5月5日的渠道開閘階段配水的流量更大,配水中后期運行的流量更小。圖7(b)中,在總干渠開閘階段,基于DCLM模型的配水流量比基于Kostiakov 公式的配水流量大15%~55%,隨著配水的進行,兩種模型的配水流量差距逐漸縮小,在圖7(c)顯示的關閘階段,基于Kostiakov 公式的配水流量比基于DCLM模型的配水流量普遍大10%左右。圖8顯示的兩種配水方式下各干渠的流量過程對比也呈現出與總干渠相似的規律。產生這種區別的原因在于:DCLM模型在計算渠道的輸水損失時考慮了渠床土壤干濕程度的影響,開閘階段渠床土壤較干燥,加大流量運行有利于渠床土壤盡快達到濕潤狀態,為減少后續時段的損失水量打下了良好基礎;配水中后期,渠床土壤已充分濕潤,無需通過加大流量以減少滲漏損失,這時減少缺水量是配水調度的主要目標。

兩種優化模型在渠道流量過程方面的差異導致了渠床土壤干濕程度變化情況亦有較大不同。如圖9所示,基于考斯加科夫公式的配水方式,在渠道輸水初期行水流量較小,導致渠床土壤濕潤較為緩慢,U1—U4的渠床土壤達到完全濕潤狀態的日期分別為5月6日、5月16日、5月11日、5月16日;基于DCLM 模型的配水方式,在渠道輸水初期加大流量運行,使得渠床土壤濕潤更為迅速,U1—U4的渠床土壤達到完全濕潤狀態的日期分別為5月1日、5月1日、5月1日、5月11日,相應提前了5~15日。

綜合圖9中渠床土壤干濕程度的變化情況以及圖6中單位渠長流量損失率的變化過程,可將渠道滲漏損失隨輸配水時間的變化過程劃分為初滲和穩滲兩個階段:初滲階段,渠床土壤由干燥狀態逐漸濕潤,渠道單位渠長流量損失率由初始較大值逐漸減小;穩滲階段,渠床土壤已達到完全濕潤狀態,渠道單位渠長流量損失率趨于穩定較低水平。在渠道初滲階段適當加大流量運行,可使渠道更快進入穩滲階段,減小渠道輸配水過程中的滲漏水量損失。

圖9 總干渠渠床土壤相對濕潤度變化對比

4 結論

本文通過引入渠床土壤前期影響含水量概念來描述渠床土壤濕潤程度動態變化,建立渠床土壤透水性參數與相對濕潤度之間的函數關系,在此基礎上提出了一種新的考慮渠系滲漏損失動態特征的渠系優化配水模型。與考斯加科夫公式等經驗公式相比,該模型能夠反映渠床土壤透水性在渠道輸水過程中的動態變化,對渠道輸水損失機理的描述更加科學準確,配水過程更為合理。在河套灌區總干渠系的夏灌輸配水實例研究中發現:與基于Kostiakov公式的配水方式相比,基于動態計算模型的渠系優化配水方式在渠道輸水初期加大了15%~55%的配水流量,使渠道更快地由初滲階段進入到穩滲階段,渠床土壤達到完全濕潤狀態和單位渠長流量損失率減小到穩定水平的日期分別提前了5~15日、5~10日,從而在不增加總干渠引水量的情況下減少渠道滲漏損失水量1560萬m3,將渠道水利用效率提高1%。

本文提出的渠系優化配水模型可為灌區制定合理的渠系輸配水計劃提供決策依據,以達到減少滲漏損失、提高渠道水利用效率的目的。但該模型對渠道滲漏損失動態變化特性的機理描述還不夠全面,僅考慮了渠床土壤透水性和渠道流量在時間上的動態變化,未能反映這些因素沿渠道流程變化對滲漏損失的影響。進一步的研究將圍繞這方面展開,以期更加全面準確地描述渠道輸水損失機理。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲无卡视频| 亚洲 欧美 偷自乱 图片| 无码综合天天久久综合网| 蜜桃视频一区| 久久久噜噜噜久久中文字幕色伊伊| 色播五月婷婷| 国产精品欧美亚洲韩国日本不卡| 国产三级毛片| 中文字幕精品一区二区三区视频| 日韩一级毛一欧美一国产| 久久窝窝国产精品午夜看片| 精品久久久久成人码免费动漫| 欧美黄色a| 国产精品免费福利久久播放| 国产精品无码一区二区桃花视频| 午夜不卡视频| 国产美女丝袜高潮| 爽爽影院十八禁在线观看| 精品人妻无码中字系列| 成人亚洲天堂| 国产精品白浆在线播放| 毛片视频网| 无码一区二区波多野结衣播放搜索| 伊人狠狠丁香婷婷综合色| 毛片最新网址| 国产最爽的乱婬视频国语对白| 国产综合亚洲欧洲区精品无码| 国产成人精品无码一区二| 国产欧美视频一区二区三区| 国产H片无码不卡在线视频| 91www在线观看| 91国内外精品自在线播放| 中文国产成人精品久久| 国产理论一区| 亚洲无码高清一区二区| 最新日韩AV网址在线观看| 国产人成在线观看| 亚洲无码日韩一区| 任我操在线视频| 国产呦精品一区二区三区下载| 无码精油按摩潮喷在线播放 | 114级毛片免费观看| 热99精品视频| 五月婷婷导航| 国产精品手机视频一区二区| 亚洲精品欧美日本中文字幕| 茄子视频毛片免费观看| 91视频青青草| 日韩精品一区二区三区大桥未久| 久久无码高潮喷水| 无码电影在线观看| 国产新AV天堂| 91欧美亚洲国产五月天| 97国产精品视频自在拍| 成人综合在线观看| 97视频在线精品国自产拍| 亚洲香蕉在线| 狠狠干综合| 激情视频综合网| 国产地址二永久伊甸园| 欧美不卡视频在线| 欧美h在线观看| yy6080理论大片一级久久| 日本黄色a视频| 青青草原国产av福利网站| 亚洲AⅤ波多系列中文字幕| 农村乱人伦一区二区| 国产在线日本| a级毛片免费看| 国产精品视频导航| 无码专区在线观看| 日本在线免费网站| 中国国产高清免费AV片| 国产极品粉嫩小泬免费看| 亚洲男人天堂久久| 亚洲国产精品国自产拍A| 一级毛片在线播放| 久久综合亚洲鲁鲁九月天| 久久精品国产精品青草app| 一级片免费网站| 一本色道久久88亚洲综合| 国产jizzjizz视频|