• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    泵站水力瞬變的二階Godunov格式模型構(gòu)建

    2022-02-06 01:02:40胡安妮吳金遠(yuǎn)
    關(guān)鍵詞:水錘二階水力

    周 領(lǐng),胡安妮,吳金遠(yuǎn)

    泵站水力瞬變的二階Godunov格式模型構(gòu)建

    周 領(lǐng)1,2,胡安妮1,吳金遠(yuǎn)1

    (1. 河海大學(xué)水利水電學(xué)院,南京 210098;2. 河海大學(xué)長江保護(hù)與綠色發(fā)展研究院,南京 210098)

    為了提高復(fù)雜泵站系統(tǒng)水力瞬變數(shù)值模擬的高效性和穩(wěn)定性,該研究基于泵站系統(tǒng)水力瞬變問題,建立有限體積法Godunov格式的數(shù)學(xué)模型,對簡單管道系統(tǒng)和復(fù)雜泵站系統(tǒng)進(jìn)行模擬研究。與常用的特征線法求解泵站水力模型方程不同,該模型引進(jìn)有限體積法二階Godunov格式對模型進(jìn)行離散,用Riemann求解器對離散通量進(jìn)行求解。使用MUSCL-Hancock方法進(jìn)行界面數(shù)值重構(gòu),采用MINMOD斜率限制器避免虛假震蕩。提出雙虛擬單元邊界處理方法,實(shí)現(xiàn)計(jì)算區(qū)域與邊界同時達(dá)到二階精度。將所建模型計(jì)算結(jié)果與精確解、經(jīng)典算例數(shù)據(jù)進(jìn)行對比,并針對庫朗數(shù)取值和計(jì)算網(wǎng)格數(shù)進(jìn)行敏感性分析。結(jié)果表明:所建模型模擬結(jié)果與精確解、經(jīng)典算例數(shù)據(jù)吻合較好;與特征線法相比,二階Godunov格式更加準(zhǔn)確、穩(wěn)定且高效。對于簡單管道系統(tǒng),特征線法計(jì)算耗時0.227 s,二階Godunov格式計(jì)算耗時0.017 s。對于實(shí)際泵站系統(tǒng),由于存在多特性的管道結(jié)構(gòu),二階Godunov格式模擬時需稍微降低庫朗數(shù)。而采用特征線法進(jìn)行泵站水力過渡過程計(jì)算時,若不調(diào)整管道長度或者波速,管道中庫朗數(shù)會小于1,在該文算例中,庫朗數(shù)為0.72~0.76,模擬計(jì)算結(jié)果偏差很大。所以需要調(diào)整局部管道長度或波速,以達(dá)到庫朗數(shù)為1的條件,這樣處理因改變管道特性而引入計(jì)算誤差。綜上,二階Godunov格式模擬方法可以更有效提高傳統(tǒng)泵站系統(tǒng)水力瞬變模擬的高效性、穩(wěn)定性以及準(zhǔn)確性。

    泵站;模型;有限體積法;二階Godunov格式;水力過渡過程;特征線法;一階Godunov格式

    0 引 言

    泵站輸水系統(tǒng)在啟、停機(jī)過程因流量變化而發(fā)生水力瞬變(或水錘),如果水力控制元件操作或水錘防護(hù)設(shè)施設(shè)置不當(dāng),會引起異常的壓力波動,從而導(dǎo)致機(jī)組損壞、爆管等重大安全事故。準(zhǔn)確模擬泵站系統(tǒng)水力瞬變對保證泵站系統(tǒng)安全穩(wěn)定運(yùn)行、實(shí)現(xiàn)智慧精準(zhǔn)化調(diào)度運(yùn)行極為重要[1-4]。

    常采用特征線法(Method of Charicteristic, MOC)對泵站系統(tǒng)水力瞬變進(jìn)行建模模擬[5]。對于簡單管道系統(tǒng),該方法具有簡單、準(zhǔn)確的優(yōu)勢。然而,在實(shí)際工程建設(shè)中,泵站系統(tǒng)會存在復(fù)雜多特性管道結(jié)構(gòu),水力計(jì)算時,為了滿足庫朗數(shù)為1的穩(wěn)定計(jì)算條件,MOC需要改變局部管段波速、網(wǎng)格長度(管道長度),產(chǎn)生計(jì)算誤差[5]。

    近年來,有限體積法逐漸被用于有壓管道水力瞬變計(jì)算。Guinot等[6]將有限體積法一階和二階格式數(shù)值解與精確解進(jìn)行比對,指出有限體積法可以準(zhǔn)確模擬簡單管路中水錘壓力波動。Zhao等[7]基于Godunov格式(有限體積法的一種數(shù)值格式)與Riemann求解器,得到了一階和二階的水錘求解格式,并指出Godunov格式可以用于水錘計(jì)算,具有較好的準(zhǔn)確性和穩(wěn)定性,但是在邊界處理時需要與特征線法相結(jié)合。趙越等[8]針對有壓管道內(nèi)水力瞬變問題,采用有限體積法Godunov格式進(jìn)行一維數(shù)學(xué)建模,發(fā)現(xiàn)一階Godunov格式與特征線法的計(jì)算精度一致,庫朗數(shù)變化會影響計(jì)算精度。Zhou等[9]提出一種基于Godunov格式的二階顯式有限體積法,并將其應(yīng)用基于特征線法的經(jīng)典離散空穴模型,發(fā)現(xiàn)有限體積法模擬值與試驗(yàn)數(shù)據(jù)吻合較好,并且可以避免虛假數(shù)值震蕩。上述研究主要針對簡單管道水力瞬變,而實(shí)際泵站系統(tǒng)含水泵、閥門、各種特性管道等,水力計(jì)算更為復(fù)雜。

    本文采用有限體積法二階Godunov格式[10-16]對泵站系統(tǒng)水力過渡過程進(jìn)行建模,提出將虛擬單元體與水庫、水泵等邊界元件相結(jié)合,通過通量求解模擬管道水錘變化。詳細(xì)分析了庫朗數(shù)和計(jì)算網(wǎng)格數(shù)對MOC方法和有限體積法Godunov格式的數(shù)值耗散影響,并采用Chaudhry泵站系統(tǒng)經(jīng)典算例數(shù)據(jù)[17]與本文有限體積二階Godunov格式計(jì)算值進(jìn)行對比。針對某實(shí)際泵站系統(tǒng),對比分析MOC方法和二階Godunov格式模擬結(jié)果。

    1 數(shù)學(xué)模型

    1.1 水錘基本控制方程

    泵站輸水系統(tǒng)中,管道內(nèi)水錘的基本控制方程[18-22]可由以下連續(xù)方程和動量方程進(jìn)行描述:

    式中是測壓管水頭,m;為沿管軸線方向的長度,m;為時間,s;是水錘波速,m/s;是重力加速度,m/s2;是斷面平均流速,m/s;是管道直徑,m;是按照Darcy-Weisbach方程求得的摩阻系數(shù)。

    將水錘基本控制方程(1)~(2)以矩陣形式表示:

    采用黎曼求解格式進(jìn)行求解,將式(3)可改寫為常系數(shù)雙曲方程(4):

    1.2 水泵控制方程

    1)機(jī)組轉(zhuǎn)動方程

    式中是機(jī)組轉(zhuǎn)動部分加上進(jìn)入該部分液體的轉(zhuǎn)化慣量,kg·m2;M是電機(jī)轉(zhuǎn)矩,N·m;是泵的軸力矩,N·m。

    根據(jù)機(jī)組轉(zhuǎn)動方程式,用泰勒展開作近似變換得到泵的轉(zhuǎn)速平衡方程為

    式中、是本時刻機(jī)組的轉(zhuǎn)速和力矩相對值;0、0是上一時刻機(jī)組的轉(zhuǎn)速和力矩相對值;T為機(jī)組的慣性時間常數(shù)(s),表示機(jī)組在額定軸力矩作用下,額定轉(zhuǎn)速減小到0所需要的時間,表征機(jī)組慣性的大小。

    2)水頭平衡方程

    聯(lián)立機(jī)組轉(zhuǎn)動方程(6)和水頭平衡方程(8),利用牛頓迭代法可求得水力瞬變過程中每一時刻機(jī)組的水頭、流量、轉(zhuǎn)速、力矩等參數(shù)。

    1.3 Godunov求解格式

    1.3.1 通量計(jì)算

    有限體積法將管道離散成若干單元體,每個控制單元內(nèi)的數(shù)值是連續(xù)的,但是在兩個單元控制單元相交的各個界面處的數(shù)值是間斷的。采用Godunov格式進(jìn)行求解時,可通過黎曼問題計(jì)算出每個界面上的離散通量,具體如下:

    1.3.2 一階Godunov格式

    一階Godunov格式通過計(jì)算單元內(nèi)變量的平均值來表示界面左右的值:

    方程(12)中邊界通量由黎曼解確定,先根據(jù)黎曼不變量方程對邊界處的離散通量進(jìn)行計(jì)算,然后再計(jì)算控制單元內(nèi)的離散通量。通過虛擬邊界處理方法,對計(jì)算區(qū)域內(nèi)邊界和控制單元的離散通量統(tǒng)一模擬。

    本文為統(tǒng)一計(jì)算整個計(jì)算區(qū)域的通量值,分別在管道上游邊界和管道下游邊界加入1個虛擬單元0和虛擬單元+1。

    由式(12)以及負(fù)向特征線方程可推導(dǎo)出如下的黎曼不變量方程:

    由式(12)以及正向特征線方程可推導(dǎo)出如下的黎曼不變量方程:

    1.3.3 二階Godunov格式

    1)數(shù)據(jù)重構(gòu)

    將式(18)代入(17)可以得到:

    2)推進(jìn)時間計(jì)算

    3)Riemann問題的近似求解

    將式(23)~(24)代入式(12)即可得到具有二階精度的Godunov格式每個計(jì)算單元界面+1/2處的數(shù)值通量。

    1.3.4 虛擬邊界

    一階Godunov格式為實(shí)現(xiàn)全瞬變區(qū)域內(nèi)的統(tǒng)一模擬,在上下游邊界處各加入一個虛擬單元。由于計(jì)算精度的提高,二階Godunov格式需在上游邊界處左邊和下游邊界處的右邊分別加入2個虛擬單元0、-1和虛擬單元+1、+2。

    對于上游水庫:

    對于下游水庫:

    對于泵站機(jī)組,將式(29)~(32)代入式(8)中即可求解虛擬邊界下的水頭平衡方程:

    1.3.5 時間積分

    一階Godunov格式通過歐拉法進(jìn)行,二階Godunov格式采用具有二階精度龍格庫塔法進(jìn)行計(jì)算,如下:

    2 算例分析

    2.1 模型驗(yàn)證

    2.1.1 簡單管道水錘問題

    設(shè)置一上游為水庫,下游為閥門的簡單管道。管道長20.0 m,水錘波速為1 319 m/s,上游水庫水位為30.0 m,初始流速為0.1 m/s,重力加速度為9.8 m/s2,計(jì)算網(wǎng)格數(shù)為16,下游閥門設(shè)置為瞬時關(guān)閉,光滑管道(無摩阻)。圖1和圖2中精確解是根據(jù)水錘基本控制方程[5]求解得到,在上文式(1)~(2)中,無摩阻項(xiàng)時,基本方程變成線性偏微分方程,理論上解析解即精確解[23-27]。則結(jié)果中的所有壓力衰減均是由于數(shù)值耗散引起的。

    圖1和圖2給出了庫朗數(shù)=0.2、0.6、1.0時,MOC、一階Godunov格式(First-order Godunov scheme, 1st-order GTS)、二階Godunov格式(Second-order Godunov scheme, 2nd-order GTS)模擬閥門處水錘計(jì)算結(jié)果。

    圖1 不同庫朗數(shù)條件下水錘壓力計(jì)算結(jié)果(計(jì)算網(wǎng)格數(shù)Ns=16)

    如圖1a所示,在=1條件下,MOC方法與一階Godunov格式模擬結(jié)果與精確值完全吻合。當(dāng)1時,MOC方法與一階Godunov格式模擬結(jié)果精度完全一致,且越小,數(shù)值衰減越嚴(yán)重。從數(shù)值模型方面分析,兩種方法是相似的。MOC模型可寫成如下形式:

    式中下標(biāo)代表待求節(jié)點(diǎn),-1代表左邊節(jié)點(diǎn),+1代表左邊節(jié)點(diǎn)。

    式(37)~(39)也可寫作式(40):

    如圖1b所示,在=1的條件下,二階Godunov格式方法模擬結(jié)果與精確值完全吻合。在<1時,與MOC和一階Godunov格式方法相比,二階Godunov格式能有效抑制數(shù)值耗散。

    此外,考慮計(jì)算網(wǎng)格數(shù)對數(shù)值模擬結(jié)果的影響,如圖2所示,增加計(jì)算網(wǎng)格數(shù)量,使=256,對比不同庫朗數(shù)條件下MOC方法和二階Godunov格式數(shù)值耗散程度。

    圖2 不同庫朗數(shù)條件下水錘壓力計(jì)算結(jié)果(計(jì)算網(wǎng)格數(shù)Ns=256)

    如圖2所示,對于MOC方法,適當(dāng)增加計(jì)算網(wǎng)格數(shù)對計(jì)算結(jié)果影響很大,網(wǎng)格數(shù)的增加會使數(shù)值耗散減小。對于二階Godunov格式,計(jì)算網(wǎng)格數(shù)的影響較小。同樣庫朗數(shù)情況下,二階Godunov格式下稀疏網(wǎng)格(=16)的模擬結(jié)果與MOC方法網(wǎng)格加密后(=256)的模擬結(jié)果基本一致。在相同精度要求下,=256時,MOC的計(jì)算耗時0.227 s,二階Godunov格式(=16)計(jì)算耗時0.017 s。這說明二階Godunov格式效率更高。

    2.1.2 泵站系統(tǒng)水力瞬變模擬驗(yàn)證

    為驗(yàn)證水泵模型的準(zhǔn)確性,本文模型計(jì)算值與Chaudhry[17]經(jīng)典數(shù)據(jù)(被學(xué)者們廣泛用于仿真模型驗(yàn)證)進(jìn)行對比。水泵系統(tǒng)包括2臺水泵,閥門,上下游水庫。水泵采用并聯(lián),水泵直接與上游水庫相連。2臺泵的參數(shù)一致:水泵中液體的綜合轉(zhuǎn)動慣量2=16.85,額定效率η=0.84,額定轉(zhuǎn)速N=1 100 r/min,額定流量Q=0.25 m3/s,額定揚(yáng)程H=60 m。管道參數(shù)如表1所示。

    表1 Chaudhry算例的管道參數(shù)[17]

    采用特征線法對該系統(tǒng)進(jìn)行水力分析,水泵轉(zhuǎn)速、流量以及水泵出口水頭計(jì)算結(jié)果如表2所示。

    如表2所示,本文提出的二階Godunov格式和已有常用的特征線法進(jìn)行水泵失電工況計(jì)算。與Chaudhry[17]計(jì)算結(jié)果對比(圖3),水泵流量峰值均為0.25 m3/s;水泵轉(zhuǎn)速峰值均為1 100 r/min;水泵揚(yáng)程峰值均為87 m。泵轉(zhuǎn)速、流量、壓力的最小值也基本吻合。因此,本文所建立的二階Godunov格式模型能準(zhǔn)確模擬泵站系統(tǒng)水力瞬變。

    2.2 復(fù)雜泵站系統(tǒng)瞬態(tài)計(jì)算

    2.2.1 泵站參數(shù)

    某泵站系統(tǒng)設(shè)置4臺設(shè)計(jì)軸功率為3 500 kW的水泵機(jī)組,采用“2用2備”運(yùn)行方式,水庫1水位6.4 m,水庫2水位3.15 m。2臺泵的參數(shù)一致:設(shè)計(jì)揚(yáng)程=52.7 m,設(shè)計(jì)流量=5.0 m3/s,額定轉(zhuǎn)速=500 r/min,轉(zhuǎn)動慣量2=14.5 t·m2,設(shè)計(jì)軸功率=3 500 kW。布置方式如圖4所示,管道參數(shù)如表3所示。

    圖3 Chaudhry泵站算例計(jì)算結(jié)果

    圖4 某泵站系統(tǒng)示意圖

    表3 泵站系統(tǒng)管道參數(shù)

    2.2.2 計(jì)算結(jié)果

    針對泵機(jī)組失電工況進(jìn)行計(jì)算分析,采用兩階段關(guān)閥規(guī)律[28-30],快關(guān)時間30 s,總關(guān)閥時間90 s,拐點(diǎn)處閥門開度為0.166 7。實(shí)際泵站系統(tǒng)水力計(jì)算時,MOC一般采用調(diào)整各管道水錘波速,以保證庫朗數(shù)為1(計(jì)算方案1:Case1);另外也可采用保持各管道水錘波速不變,庫朗數(shù)小于1(計(jì)算方案2:Case2)。

    本文主要計(jì)算庫朗數(shù)變化對MOC方法、二階Godunov格式數(shù)值耗散的影響,以及2種計(jì)算方案對計(jì)算結(jié)果的影響,相應(yīng)庫朗數(shù)取值見表4。

    表4 兩種計(jì)算方案Cr取值

    如表5所示,方案1中MOC方法對管道波速進(jìn)行調(diào)整,二階Godunov格式方法管道波速取原值。方案1管道系統(tǒng)參數(shù)見表5。

    2種方法模擬泵站系統(tǒng)發(fā)生失電后120 s水力瞬變過程。整個計(jì)算時段內(nèi)水泵最大反轉(zhuǎn)轉(zhuǎn)速、最大倒流流量、水泵出口閥后最小壓力、水泵出口閥后最大壓力見表6所示。

    表5 計(jì)算方案1中管道系統(tǒng)參數(shù)處理

    表6 兩種計(jì)算方案下不同計(jì)算方法的計(jì)算結(jié)果

    如表6所示,在實(shí)際泵站系統(tǒng)中,庫朗數(shù)條件的變化對于二階Godunov格式模擬結(jié)果影響很小,水泵最大反轉(zhuǎn)轉(zhuǎn)速穩(wěn)定在?417 r/min,水泵最大倒流流量穩(wěn)定在?3.27 m3/s,揚(yáng)程峰值變化范圍很小,這說明二階Godunov格式模擬結(jié)果穩(wěn)定,這與簡單管道結(jié)論一致。然而,對于MOC方法,方案1的管道波速經(jīng)調(diào)整后,均高于實(shí)際管道波速,故水泵出口閥后壓力峰值增大,而方案2采用管道實(shí)際波速,揚(yáng)程峰值減小3.13 m,倒流流量變化至?3.88 m3/s,庫朗數(shù)條件的變化使水泵最大反轉(zhuǎn)轉(zhuǎn)速、最大倒流流量、水泵出口閥后最大、最小壓力的峰值和周期均出現(xiàn)較大誤差。因此對于水泵系統(tǒng)水力瞬變模擬,本文采用的二階Godunov格式更具有穩(wěn)定性、準(zhǔn)確性。機(jī)組轉(zhuǎn)速、流量、水泵出口閥后壓力變化過程線見圖5。

    圖5 泵站系統(tǒng)計(jì)算結(jié)果

    為考慮網(wǎng)格數(shù)對計(jì)算結(jié)果的影響,在上文2個算例基礎(chǔ)上設(shè)置方案3和方案 4,參數(shù)設(shè)置如表7所示。

    表7 方案3、4的參數(shù)設(shè)置

    如圖6a所示,對于MOC方法,庫朗數(shù)小于1時,增加網(wǎng)格數(shù)可以適當(dāng)提高計(jì)算結(jié)果精確度,但是計(jì)算耗時將成倍增加(MOC中Case2和Case4的計(jì)算耗時分別為0.063和0.094 s)。如圖6b所示,對于二階Godunov格式,庫朗數(shù)等于或小于1,網(wǎng)格數(shù)對計(jì)算結(jié)果精確度影響均很小,但計(jì)算耗時會增大(二階Godunov中Case2和Case4的計(jì)算耗時分別為0.077 和0.109 s)。因此,當(dāng)模擬水泵系統(tǒng)瞬變過程面臨庫朗數(shù)小于1情況,為達(dá)到相同精度要求,MOC需要更精細(xì)網(wǎng)格,其計(jì)算耗時將大于二階Godunov格式。

    圖6 泵站系統(tǒng)水泵轉(zhuǎn)速模擬結(jié)果

    3 結(jié) 論

    本文采用有限體積法二階Godunov格式對泵站系統(tǒng)水力瞬變進(jìn)行建模模擬,管道邊界采用虛擬邊界的處理方法,水力元件邊界采用傳統(tǒng)的特征線法進(jìn)行處理。將數(shù)值模擬結(jié)果與MOC計(jì)算所得結(jié)果進(jìn)行對比分析。主要結(jié)論:

    1)本文建立的二階Godunov格式可以準(zhǔn)確模擬水泵系統(tǒng)水力瞬變過程中水泵轉(zhuǎn)速、流量、壓力等參數(shù)變化過程線。通過與已有文獻(xiàn)計(jì)算結(jié)果對比可知,本文建立模型模擬結(jié)果精確,水泵流量峰值均為0.25 m3/s;水泵轉(zhuǎn)速峰值均為1 100 r/min;水泵揚(yáng)程峰值均為87 m。泵轉(zhuǎn)速、流量、壓力的最小值也基本吻合。

    2)與MOC模型相比,二階Godunov模型具有更好的穩(wěn)定性、準(zhǔn)確性以及高效性。當(dāng)庫朗數(shù)<1時,二階Godunov格式能保持較高的計(jì)算精度、結(jié)果穩(wěn)定性。簡單管道系統(tǒng)中,在相同精度要求下,同一臺計(jì)算機(jī)中特征線法計(jì)算耗時是0.227 s,二階Godunov格式計(jì)算耗時是0.017 s。復(fù)雜泵站系統(tǒng),相同精度要求下,MOC計(jì)算耗時在二階Godunov格式基礎(chǔ)上成倍增加,二階Godunov格式計(jì)算效率更高。

    3)在實(shí)際泵站系統(tǒng)中,為避免數(shù)值耗散,MOC模型需要對管道波速和管道分段進(jìn)行前處理,例如調(diào)整波速或者網(wǎng)格長度,或者直接忽略短管,這樣處理不僅繁瑣,而且因改變管道特性而引入計(jì)算誤差。對于有限體積法二階Godunov格式,管道波速或者管道分段無需調(diào)整,因?yàn)閹炖蕯?shù)變化對于二階Godunov格式模擬結(jié)果影響很小。在本文實(shí)際算例中,水泵最大反轉(zhuǎn)轉(zhuǎn)速穩(wěn)定在?417 r/min,水泵最大倒流流量穩(wěn)定在?3.27 m3/s,揚(yáng)程峰值變化范圍很小,這說明二階Godunov格式模擬結(jié)果穩(wěn)定。既簡單方便,又可保證計(jì)算精度。

    [1] 劉金昊,吳建華. 西山供水工程事故停泵水力過渡過程計(jì)算及水錘防護(hù)[J]. 水電能源科學(xué),2021,39(7):113-116.

    Liu Jinhao, Wu Jianhua. Calculation of hydraulic transition process and water hammer protection for accidental pump-stop in xishan water supply project [J]. Water Resources and Power, 2021, 39(7): 113-116. (in Chinese with English abstract)

    [2] 黃偉,楊開林,郭新蕾,等. 抽水蓄能電站極端甩負(fù)荷工況球閥協(xié)同調(diào)節(jié)[J]. 清華大學(xué)學(xué)報(bào)(自然科學(xué)版),2019,59(8):635-644.

    Huang Wei, Yang Kailin, Guo Xinlei, et al. Coordinated regulation of ball valves in pumped storage stations for extreme conditions[J]. Journal of Tsinghua University (Science & Technology), 2019, 59(8): 635-644. (in Chinese with English abstract)

    [3] 楊玉思,徐艷艷,羨巨智. 長距離高揚(yáng)程多起伏輸水管道水錘防護(hù)的研究[J]. 給水排水,2009,45(4):108-211.

    Yang Yusi, Xu Yanyan, Xian Juzhi. Research on water hammer prevention in high lift, hilly and long distance water transmission pipeline[J]. Water & Wastewater Engineering, 2009, 45(4): 108-211. (in Chinese with English abstract)

    [4] 方永,周宏. 管道系統(tǒng)中的水錘及其防護(hù)研究[J]. 石化技術(shù),2016,23(2):115.

    Fang Yong, Zhou Hong. Water hammer and preventive measures in pipeline[J]. Petrochemical Industry Technology, 2016, 23(2): 115. (in Chinese with English abstract)

    [5] Wylie E B, Streeter V L, Suo L. Fluid transients in systems[M]. Englewood Cliffs, NJ: Prentice Hall, 1993.

    [6] Guinot V. Riemann solvers for water hammer simulations by Godunov method[J]. International Journal for Numerical Methods in Engineering, 2000, 49(7): 851-870.

    [7] Zhao M, Ghidaoui M S. Godunov-type solutions for water hammer flows[J]. Journal of Hydraulic Engineering-ASCE, 2004, 130(4): 341-348.

    [8] 趙越,周領(lǐng),劉德有,等. 基于有限體積法Godunov格式的水錘計(jì)算模型[J]. 水利水電科技進(jìn)展,2019,39(1):76-81.

    Zhao Yue, Zhou Ling, Liu Deyou, et al. Water hammer model based on finite volume method and Godunov-type scheme[J]. Advances in Science and Technology of Water Resources. 2019, 39(1): 76-81. (in Chinese with English abstract)

    [9] Zhou L, Wang H, Liu D Y, et al. A second-order finite volume method for pipe flow with water column separation [J]. Journal of Hydro-environment Research, 2016, 17: 47-55.

    [10] 向小華,吳曉玲,牛帥,等. 基于顯式有限體積法的一維河網(wǎng)模型[J]. 水利水電科技進(jìn)展,2015,35(4):6-9,43.

    Xiang Xiaohua, Wu Xiaoling, Niu Shuai, et al. Construction of one-dimensional river network model based on explicit finite volume method[J]. Advances in Science and Technology of Water Resources, 2015, 35(4): 6-9, 43. (in Chinese with English abstract)

    [11] Hu Y Y, Zhou L, Pan T W, et al. Godunov-type solutions for free surface transient flow in pipeline incorporating unsteady friction[J]. Journal of Water Supply: Research and Technology-Aqua, 2022, 71(4): 546-562.

    [12] 鄭劼恒,蔣明,郭芮,等.順序輸送管道水力瞬變模擬的有限體積法[J]. 計(jì)算力學(xué)學(xué)報(bào),2015,32(3):418-422,428.

    Zheng Jieheng, Jiang Ming, Guo Rui, et al. Finite volume method for hydraulic transient simulation of batching pipeline[J]. Chinese Journal of Computational Mechanics, 2015, 32(3): 418-422, 428. (in Chinese with English abstract)

    [13] 楊開林. 電站與泵站中的水力瞬變及調(diào)節(jié)[M]. 北京:中國水利水電出版社,2000.

    [14] Toro E F. Riemann solvers and numerical methods for fluid dynamics[M]. Berlin: Verlag Berlin Heidelberg, 2009: 115-125.

    [15] 劉錦濤,章少輝,許迪,等. 灌溉輸配水系統(tǒng)明滿流的全隱式耦合模擬及驗(yàn)證[J]. 農(nóng)業(yè)工程學(xué)報(bào),2017,33(19):124-130.

    Liu Jintao, Zhang Shaohui, Xu Di, et al. Coupled simulation and validation with fully implicit time scheme for free-surface-pressurized water flow in pipe/channel[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2017, 33(19): 124-130. (in Chinese with English abstract)

    [16] 王振華,馬習(xí)賀,李文昊,等.基于改進(jìn)4-方程摩擦模型的輸水管道水錘壓力計(jì)算[J]. 農(nóng)業(yè)工程學(xué)報(bào),2018,34(7):114-120.

    Wang Zhenhua, Ma Xihe, Li Wenhao, et al. Calculation of water hammer pressure of flow pipeline based on modified four-equation friction model[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2018, 34(7): 114-120. (in Chinese with English abstract)

    [17] Chaudhry M H. 實(shí)用水力瞬變過程[M]. 北京:中國水利水電出版社,2015.

    [18] Lyu J W, Zhang J, Wang T Y. Study on water hammer protection of the siphon breaking structure in the water supply system[J]. AQUA-Water Infrastructure, Ecosystems and Society, 2022, 71(3): 478-489.

    [19] Abdel-Fatah M A, Amin A, Al Bazedi G A. Model and protected design of water piping system to minimize the water hammer effect[J]. Egyptian Journal of Chemistry, 2022, 65(1): 1-2.

    [20] Lupa S I, Gagnon M, Muntean S, et al. The impact of water hammer on hydraulic power units[J]. Energies, 2022, 15(4): 1526.

    [21] 富友,蔣勁,李燕輝,等. 改進(jìn)雙流體模型計(jì)算有液柱分離的管路水錘壓力[J]. 農(nóng)業(yè)工程學(xué)報(bào),2018,34(15):58-65.

    Fu You, Jiang Jing, Li Yanhui, et al. Calculation of pipe water hammer pressure with liquid column separation by improved two-fluid model[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2018, 34(15): 58-65. (in Chinese with English abstract)

    [22] 周領(lǐng),陸燕清.排水管道瞬變流的SWMM模擬能力研究[J]. 中國給水排水,2022,38(5):108-115.

    Zhou Ling, Lu Yanqing. SWMM simulation capability for transient fow in drainage pipe[J]. Water & Wastewater Engineering, 2022, 38(5): 108-115. (in Chinese with English abstract)

    [23] 吳峰. 輸水管道中含氣水流瞬變特性試驗(yàn)研究及數(shù)值模擬[D]. 西安:西安理工大學(xué),2021.

    Wu Feng. Experimental Research and Numerical Simulation on Transient Characteristics of Air-Bearing Water Flow in Water Pipeline[D]. Xi’an: Xi’an University of Technology, 2021. (in Chinese with English abstract)

    [24] 丁澤. 有壓管道中閥門工作狀態(tài)檢測研究[D]. 哈爾濱:哈爾濱商業(yè)大學(xué),2019.

    Ding Ze. Research on Valve Working State Detection in Pressurized Pipeline[D]. Harbin: Harbin University of Commerce, 2019. (in Chinese with English abstract)

    [25] 趙思茂. 輸水管道內(nèi)瞬變流特性的三維數(shù)值模擬研究[D]. 西安:西安理工大學(xué),2020.

    Zhao Simao. Three Dimensional Numerical Simulation of Transient Flow in Water Pipeline[D]. Xi’an: Xi’an University of Technology, 2020. (in Chinese with English abstract)

    [26] Al-Khayat R H, Kadhim A A, Al-Baghdadi M A R S, et al. Flow parameters effect on water hammer stability in hydraulic system by using state-space method[J]. Open Engineering, 2022, 12(1): 215-226.

    [27] Waqar M, Louati M, Ghidaoui M S. Time-reversal of water-hammer waves[J]. Journal of Hydraulic Research, 2022, 60(1): 25-45.

    [28] Li Y J, Hu X D, Zhou F J, et al. A new comprehensive filtering model for pump shut-in water hammer pressure wave signals during hydraulic fracturing[J]. Journal of Petroleum Science and Engineering, 2022, 208: 109796.

    [29] Kim S G, Lee K B, Kim K Y. Water hammer in the pump-rising pipeline system with an air chamber[J]. Journal of Hydrodynamics, 2014, 26(6): 960-964.

    [30] Wan W Y, Li F Q. Sensitivity analysis of operational time differences for a pump-valve system on a water hammer response[J]. Journal of Pressure Vessel Technology, 2016, 138(1): 011303.

    Construction of the second-order Godunov scheme model for hydraulic transients in pumping stations

    Zhou Ling1,2, Hu Anni1, Wu Jinyuan1

    (1.,,210098,; 2.,,210098,)

    This study aims to implement a more efficient and stable numerical simulation of the hydraulic transient in a complex pumping station system. A finite volume method (FVM) Godunov scheme was established to simulate the simple pipeline and complex pumping station system. The FVM was then introduced to discretize the mathematical models, while the Riemann solver was selected to solve the discrete flux. The MUSCL-Hancock method was utilized to reconstruct the numerical data at the interface of control volumes. The higher numerical accuracy and stability were realized in the Godunov scheme, compared with the frequently-used method of characteristics. Meanwhile, the MINMOD slope limiter was used to avoid false oscillation. The boundary processing of the dual virtual unit was then presented for the second-order accuracy of both the computational region and the boundary, particularly for the simpler computation. The simulation of the improved model was in good agreement with the exact solution and the classical examples. The sensitivity analysis was also performed on the Courant and grid number. Furthermore, a more accurate, stable, and efficient performance was achieved in the second-order Godunov scheme, compared with the method of characteristics. More importantly, there was more outstanding attenuation with the decrease of the Courant number for a simple pipeline system. The computation time of the second-order Godunov scheme was 0.017 s at the same accuracy, compared with the method of characteristics (0.227 s). Consequently, a more stable and efficient performance was achieved in the second-order Godunov scheme. In the actual pumping system with the multiple-characteristics pipe structure, the second-order Godunov scheme required only a slight reduction in the Courant numbers, indicating the simple and convenient way for high numerical accuracy. Once the method of characteristics was used to calculate the hydraulic transition of the pumping station, the Courant number in the pipeline was less than 1 at the same length or wave velocity of the pipeline. By contrast, the Courant number was 0.72-0.76 in this case, indicating a very different simulation from the actual. Therefore, it is necessary to adjust the local pipeline length or wave velocity for the condition that the Courant number was 1. The tedious operation can lead to calculation errors, due to the change in pipeline characteristics. The accuracy can be improved but with less computational efficiency, if the wave velocity remained unchanged to increase the number of computational grids. In the method of characteristics, the number of grids can properly improve the accuracy of the calculation but with the doubled computation time, when the Courant number was less than 1. In the second-order Godunov scheme, there was little effect of grid number on the accuracy of the calculation but with the longer calculation time, whether the Courant number was equal to or less than 1. Therefore, a finer grid was preferred in the method of characteristics for the same accuracy requirements, when the Courant number was less than 1 in the transient process of the simulated pump system. Therefore, the second-order Godunov scheme can accurately simulate the process lines of rotational speed, discharge, and outlet pressure parameters during the hydraulic transient of the pump system. Anyway, the second-order Godunov scheme can be expected to effectively improve the efficiency, stability, and accuracy of hydraulic transient simulation of traditional pumping station systems.

    pump station; models; finite volume method; second-order Godunov scheme; hydraulic transition process; method of characteristics; first-order Godunov scheme

    10.11975/j.issn.1002-6819.2022.19.005

    TV143.1

    A

    1002-6819(2022)-19-0042-09

    周領(lǐng),胡安妮,吳金遠(yuǎn). 泵站水力瞬變的二階Godunov格式模型構(gòu)建[J]. 農(nóng)業(yè)工程學(xué)報(bào),2022,38(19):42-50.doi:10.11975/j.issn.1002-6819.2022.19.005 http://www.tcsae.org

    Zhou Ling, Hu Anni, Wu Jinyuan. Construction of the second-order Godunov scheme model for hydraulic transients in pumping stations[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(19): 42-50. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2022.19.005 http://www.tcsae.org

    2022-05-12

    2022-08-24

    國家自然科學(xué)基金項(xiàng)目(51679066,51839008);霍英東教育基金會青年教師基金(161068)

    周領(lǐng),博士,教授,研究方向?yàn)樗娬尽⒈谜?、輸水工程水力學(xué)。Email:zlhhu@163.com

    猜你喜歡
    水錘二階水力
    水力全開
    高水頭短距離泵站水錘計(jì)算分析
    一類二階迭代泛函微分方程的周期解
    一類二階中立隨機(jī)偏微分方程的吸引集和擬不變集
    二階線性微分方程的解法
    一類二階中立隨機(jī)偏微分方程的吸引集和擬不變集
    水錘防護(hù)措施在某長距離供水系統(tǒng)中的應(yīng)用
    球墨鑄鐵管的水力計(jì)算
    輸水管道水錘防護(hù)方案對比分析
    水力噴射壓裂中環(huán)空水力封隔全尺寸實(shí)驗(yàn)
    成年女人毛片免费观看观看9| 久久精品人妻少妇| 欧美一区二区精品小视频在线| 久久精品国产亚洲av香蕉五月| 嫩草影院精品99| 十八禁人妻一区二区| 又爽又黄无遮挡网站| 欧美最黄视频在线播放免费| 久久久久久九九精品二区国产| 久久中文看片网| 91在线精品国自产拍蜜月 | 激情在线观看视频在线高清| 日本与韩国留学比较| 色综合站精品国产| 别揉我奶头~嗯~啊~动态视频| 观看免费一级毛片| 色哟哟哟哟哟哟| 亚洲av美国av| av女优亚洲男人天堂| 国产色爽女视频免费观看| av福利片在线观看| 国产免费av片在线观看野外av| 久久久久久久精品吃奶| 成年女人毛片免费观看观看9| 国产亚洲精品一区二区www| 嫩草影视91久久| 九色成人免费人妻av| 日韩欧美精品免费久久 | 国语自产精品视频在线第100页| 亚洲国产精品sss在线观看| 亚洲国产欧美人成| av片东京热男人的天堂| 搡老熟女国产l中国老女人| 欧美+亚洲+日韩+国产| 国产成人福利小说| 精品欧美国产一区二区三| 日韩欧美国产在线观看| 在线观看66精品国产| 国产高清激情床上av| 欧美激情在线99| 亚洲欧美精品综合久久99| 国产伦一二天堂av在线观看| 久久精品国产自在天天线| 亚洲欧美一区二区三区黑人| 亚洲精品乱码久久久v下载方式 | x7x7x7水蜜桃| 99国产精品一区二区三区| 国产男靠女视频免费网站| 噜噜噜噜噜久久久久久91| 亚洲成人中文字幕在线播放| 日韩欧美免费精品| 精品日产1卡2卡| 内射极品少妇av片p| 村上凉子中文字幕在线| 99国产精品一区二区蜜桃av| 日韩欧美国产在线观看| 99久久久亚洲精品蜜臀av| 亚洲国产欧洲综合997久久,| 国内精品久久久久久久电影| 熟妇人妻久久中文字幕3abv| 草草在线视频免费看| 免费人成视频x8x8入口观看| 精品久久久久久成人av| 国产毛片a区久久久久| 美女免费视频网站| 亚洲av熟女| 亚洲成av人片免费观看| 国产成人福利小说| 欧美日韩中文字幕国产精品一区二区三区| 在线观看av片永久免费下载| 色在线成人网| 国产精品久久久人人做人人爽| 久久久久久久久久黄片| 有码 亚洲区| 亚洲五月婷婷丁香| 一进一出抽搐动态| 日本撒尿小便嘘嘘汇集6| 99久国产av精品| 亚洲在线观看片| 欧美区成人在线视频| 此物有八面人人有两片| 国产91精品成人一区二区三区| 亚洲内射少妇av| 久久久国产精品麻豆| 可以在线观看的亚洲视频| 日日摸夜夜添夜夜添小说| 悠悠久久av| 久久久国产精品麻豆| 亚洲成a人片在线一区二区| 亚洲一区二区三区不卡视频| 九九久久精品国产亚洲av麻豆| 一区二区三区免费毛片| a级一级毛片免费在线观看| 日韩人妻高清精品专区| 综合色av麻豆| 精品99又大又爽又粗少妇毛片 | 国产成人啪精品午夜网站| 免费av不卡在线播放| 国产欧美日韩一区二区三| 18禁裸乳无遮挡免费网站照片| 黑人欧美特级aaaaaa片| 日韩免费av在线播放| 欧美性猛交黑人性爽| 在线看三级毛片| 国产精品亚洲美女久久久| 99国产极品粉嫩在线观看| 制服丝袜大香蕉在线| 99热6这里只有精品| 亚洲无线观看免费| 色尼玛亚洲综合影院| 真人做人爱边吃奶动态| 99国产极品粉嫩在线观看| av在线天堂中文字幕| 91麻豆av在线| av在线蜜桃| 九九热线精品视视频播放| 国产成人欧美在线观看| 女同久久另类99精品国产91| 亚洲无线在线观看| 99久国产av精品| 一级黄色大片毛片| 午夜福利免费观看在线| 成年女人毛片免费观看观看9| 国产黄片美女视频| 欧美黑人欧美精品刺激| 美女 人体艺术 gogo| 国产午夜精品论理片| 一区福利在线观看| a在线观看视频网站| 一级毛片女人18水好多| 人人妻,人人澡人人爽秒播| 制服丝袜大香蕉在线| 日本三级黄在线观看| 国产99白浆流出| 成人特级av手机在线观看| 亚洲一区高清亚洲精品| 亚洲中文字幕日韩| 我要搜黄色片| 国产av一区在线观看免费| 两个人的视频大全免费| 亚洲欧美一区二区三区黑人| 男女那种视频在线观看| ponron亚洲| 搡女人真爽免费视频火全软件 | 美女高潮的动态| 99久久九九国产精品国产免费| 久久久成人免费电影| 给我免费播放毛片高清在线观看| 好看av亚洲va欧美ⅴa在| 最新中文字幕久久久久| 成年人黄色毛片网站| 精品久久久久久久人妻蜜臀av| 国产视频一区二区在线看| 国产主播在线观看一区二区| 在线观看美女被高潮喷水网站 | 日本免费a在线| 亚洲人与动物交配视频| 内地一区二区视频在线| 国产伦精品一区二区三区四那| 在线a可以看的网站| 日韩欧美精品v在线| 国产一区二区激情短视频| www.www免费av| 狂野欧美激情性xxxx| 小蜜桃在线观看免费完整版高清| 成年女人看的毛片在线观看| 在线观看一区二区三区| 中文字幕熟女人妻在线| 在线观看日韩欧美| 日本一本二区三区精品| 白带黄色成豆腐渣| 五月玫瑰六月丁香| 首页视频小说图片口味搜索| 久久香蕉精品热| 国产一区二区亚洲精品在线观看| 99久国产av精品| 亚洲电影在线观看av| 亚洲性夜色夜夜综合| 少妇熟女aⅴ在线视频| 日韩欧美 国产精品| 精品久久久久久久毛片微露脸| 亚洲国产中文字幕在线视频| 日韩精品青青久久久久久| 18禁美女被吸乳视频| 观看免费一级毛片| 村上凉子中文字幕在线| 1024手机看黄色片| 亚洲av五月六月丁香网| 嫩草影院精品99| 国产av一区在线观看免费| 大型黄色视频在线免费观看| 18美女黄网站色大片免费观看| av专区在线播放| 亚洲成人中文字幕在线播放| 三级毛片av免费| 少妇高潮的动态图| 欧美一区二区国产精品久久精品| 制服人妻中文乱码| 亚洲精品在线观看二区| 毛片女人毛片| 成人国产综合亚洲| 欧美性猛交╳xxx乱大交人| 日本三级黄在线观看| 亚洲 国产 在线| 亚洲aⅴ乱码一区二区在线播放| 夜夜躁狠狠躁天天躁| 久久国产乱子伦精品免费另类| 99riav亚洲国产免费| 中文资源天堂在线| www.999成人在线观看| 美女高潮的动态| 欧美三级亚洲精品| 一个人观看的视频www高清免费观看| 人妻丰满熟妇av一区二区三区| 成人三级黄色视频| 高清毛片免费观看视频网站| 欧美日韩瑟瑟在线播放| 国产成人福利小说| 国产 一区 欧美 日韩| av天堂在线播放| 高清毛片免费观看视频网站| 成人国产综合亚洲| 18禁在线播放成人免费| 九色国产91popny在线| 午夜影院日韩av| 美女cb高潮喷水在线观看| 熟妇人妻久久中文字幕3abv| 久久久国产精品麻豆| 精品久久久久久久久久久久久| 亚洲精品亚洲一区二区| 日韩高清综合在线| 757午夜福利合集在线观看| 亚洲精品粉嫩美女一区| 亚洲天堂国产精品一区在线| av欧美777| 狂野欧美白嫩少妇大欣赏| 免费av不卡在线播放| 99国产综合亚洲精品| 久久国产精品人妻蜜桃| 久久这里只有精品中国| 成年女人毛片免费观看观看9| 麻豆久久精品国产亚洲av| 亚洲av熟女| 少妇的逼水好多| 91麻豆精品激情在线观看国产| 两人在一起打扑克的视频| 一区二区三区激情视频| 一个人观看的视频www高清免费观看| 色综合婷婷激情| 亚洲激情在线av| 欧美3d第一页| 国产三级中文精品| 国产淫片久久久久久久久 | 狂野欧美白嫩少妇大欣赏| 麻豆成人av在线观看| 香蕉丝袜av| 天堂网av新在线| avwww免费| 免费在线观看成人毛片| 日韩免费av在线播放| 国模一区二区三区四区视频| 国产精品国产高清国产av| 色综合站精品国产| 窝窝影院91人妻| 亚洲激情在线av| 亚洲第一欧美日韩一区二区三区| 精品人妻偷拍中文字幕| 亚洲av免费高清在线观看| 少妇高潮的动态图| 两个人视频免费观看高清| 小蜜桃在线观看免费完整版高清| 国产成人福利小说| 观看免费一级毛片| 日韩欧美一区二区三区在线观看| 精品99又大又爽又粗少妇毛片 | 一区二区三区高清视频在线| 十八禁人妻一区二区| 国内精品一区二区在线观看| 欧美一区二区亚洲| 中文字幕av在线有码专区| 国产精品一及| 级片在线观看| 久久婷婷人人爽人人干人人爱| 国产成人欧美在线观看| 国产午夜精品久久久久久一区二区三区 | 国产真实乱freesex| 波野结衣二区三区在线 | 国产激情偷乱视频一区二区| 日本精品一区二区三区蜜桃| 国产成人啪精品午夜网站| 亚洲成人精品中文字幕电影| 成人性生交大片免费视频hd| 国产视频一区二区在线看| 久久久久久久精品吃奶| 久久久久国产精品人妻aⅴ院| 亚洲无线在线观看| 99久久精品一区二区三区| 久久中文看片网| 亚洲国产精品久久男人天堂| 日本 av在线| 免费看日本二区| 此物有八面人人有两片| 麻豆成人午夜福利视频| 国产亚洲精品久久久久久毛片| 日韩av在线大香蕉| 欧美精品啪啪一区二区三区| 观看美女的网站| 人人妻人人看人人澡| 99精品欧美一区二区三区四区| 法律面前人人平等表现在哪些方面| 国产野战对白在线观看| 最新在线观看一区二区三区| 人人妻人人澡欧美一区二区| 美女 人体艺术 gogo| 免费在线观看影片大全网站| 一本一本综合久久| 动漫黄色视频在线观看| 午夜老司机福利剧场| 国产精品国产高清国产av| 在线观看美女被高潮喷水网站 | 国产乱人伦免费视频| 日本a在线网址| 亚洲国产日韩欧美精品在线观看 | 久久久久国产精品人妻aⅴ院| 偷拍熟女少妇极品色| 日韩欧美在线乱码| 淫秽高清视频在线观看| 午夜日韩欧美国产| 国产精品自产拍在线观看55亚洲| 色综合婷婷激情| 亚洲国产欧美人成| 在线国产一区二区在线| 19禁男女啪啪无遮挡网站| 真实男女啪啪啪动态图| 88av欧美| 黑人欧美特级aaaaaa片| 亚洲精品美女久久久久99蜜臀| 尤物成人国产欧美一区二区三区| 免费av毛片视频| 欧美日韩亚洲国产一区二区在线观看| 亚洲成人中文字幕在线播放| 欧美另类亚洲清纯唯美| 搡女人真爽免费视频火全软件 | av天堂在线播放| 51国产日韩欧美| 国内少妇人妻偷人精品xxx网站| 国产亚洲av嫩草精品影院| 变态另类成人亚洲欧美熟女| 18禁裸乳无遮挡免费网站照片| 黄色成人免费大全| 一个人观看的视频www高清免费观看| 久久久久国产精品人妻aⅴ院| 51午夜福利影视在线观看| 国产精品99久久久久久久久| 精品免费久久久久久久清纯| 手机成人av网站| 岛国在线免费视频观看| 白带黄色成豆腐渣| 在线观看午夜福利视频| 全区人妻精品视频| 免费观看的影片在线观看| 一区二区三区激情视频| 最近在线观看免费完整版| 国产99白浆流出| 亚洲最大成人手机在线| 色哟哟哟哟哟哟| 真人一进一出gif抽搐免费| 精品一区二区三区人妻视频| 国产淫片久久久久久久久 | 18禁黄网站禁片免费观看直播| www国产在线视频色| 精品久久久久久,| 免费看日本二区| 亚洲成人中文字幕在线播放| 欧洲精品卡2卡3卡4卡5卡区| 宅男免费午夜| 99热精品在线国产| 国产一区二区在线av高清观看| av黄色大香蕉| 丰满的人妻完整版| 国产又黄又爽又无遮挡在线| 免费高清视频大片| 日日摸夜夜添夜夜添小说| 亚洲精品国产精品久久久不卡| 国产激情欧美一区二区| 亚洲狠狠婷婷综合久久图片| 在线观看舔阴道视频| 国内毛片毛片毛片毛片毛片| 日韩国内少妇激情av| 欧美成人性av电影在线观看| 国产精品99久久99久久久不卡| 欧美xxxx黑人xx丫x性爽| 在线a可以看的网站| 俺也久久电影网| 999久久久精品免费观看国产| 久久九九热精品免费| 色哟哟哟哟哟哟| 成人性生交大片免费视频hd| 婷婷精品国产亚洲av在线| xxxwww97欧美| 欧美bdsm另类| 精品久久久久久久久久久久久| 色吧在线观看| 可以在线观看毛片的网站| 欧美一区二区国产精品久久精品| 国产伦一二天堂av在线观看| 99久久无色码亚洲精品果冻| 日本黄色视频三级网站网址| 亚洲国产色片| 亚洲精品在线观看二区| 久久精品国产99精品国产亚洲性色| av国产免费在线观看| 亚洲男人的天堂狠狠| 色综合亚洲欧美另类图片| 真实男女啪啪啪动态图| 国产黄a三级三级三级人| av专区在线播放| 亚洲av电影在线进入| 好看av亚洲va欧美ⅴa在| 日本a在线网址| 亚洲av二区三区四区| www.www免费av| 国内精品美女久久久久久| 亚洲成人久久爱视频| 婷婷精品国产亚洲av在线| 又爽又黄无遮挡网站| 亚洲美女黄片视频| 毛片女人毛片| 黄色日韩在线| 少妇人妻一区二区三区视频| 18禁国产床啪视频网站| 国产精品永久免费网站| 免费在线观看亚洲国产| 国产探花极品一区二区| 亚洲精品国产精品久久久不卡| 国产精品一区二区免费欧美| 男人舔女人下体高潮全视频| 久久午夜亚洲精品久久| 精品一区二区三区视频在线观看免费| 最近最新中文字幕大全免费视频| 麻豆一二三区av精品| 男女床上黄色一级片免费看| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 久久久国产成人精品二区| 日韩中文字幕欧美一区二区| 九九热线精品视视频播放| 国产av麻豆久久久久久久| 波多野结衣高清作品| 国产色婷婷99| 美女高潮喷水抽搐中文字幕| 欧美成人免费av一区二区三区| 亚洲电影在线观看av| 国产精品久久视频播放| 国产探花极品一区二区| 两个人视频免费观看高清| 制服人妻中文乱码| 91麻豆精品激情在线观看国产| 欧美国产日韩亚洲一区| 亚洲精品日韩av片在线观看 | 成年女人看的毛片在线观看| а√天堂www在线а√下载| www日本在线高清视频| 亚洲精华国产精华精| 三级国产精品欧美在线观看| 最新在线观看一区二区三区| www.色视频.com| 欧美一级a爱片免费观看看| 九九热线精品视视频播放| 天天一区二区日本电影三级| 日韩欧美精品免费久久 | 综合色av麻豆| 亚洲av五月六月丁香网| 亚洲欧美日韩高清专用| 无遮挡黄片免费观看| 国产av一区在线观看免费| 床上黄色一级片| 又黄又粗又硬又大视频| 午夜免费激情av| 成人国产一区最新在线观看| 欧美成人免费av一区二区三区| 美女黄网站色视频| 91在线观看av| 久久天躁狠狠躁夜夜2o2o| 麻豆久久精品国产亚洲av| 成年免费大片在线观看| 91麻豆精品激情在线观看国产| av专区在线播放| 亚洲av第一区精品v没综合| 午夜两性在线视频| aaaaa片日本免费| 搡女人真爽免费视频火全软件 | 丰满人妻一区二区三区视频av | 亚洲av美国av| 国产午夜精品论理片| 中国美女看黄片| 大型黄色视频在线免费观看| 精品99又大又爽又粗少妇毛片 | 香蕉av资源在线| 午夜视频国产福利| 免费搜索国产男女视频| ponron亚洲| 色综合欧美亚洲国产小说| 亚洲精品色激情综合| 久久久久久久亚洲中文字幕 | 久久精品国产99精品国产亚洲性色| 国产欧美日韩一区二区三| 少妇的逼水好多| 男女视频在线观看网站免费| 每晚都被弄得嗷嗷叫到高潮| 老司机午夜福利在线观看视频| 精品乱码久久久久久99久播| 国产亚洲精品综合一区在线观看| 天天躁日日操中文字幕| 成人国产一区最新在线观看| 久久伊人香网站| 国产一区二区在线av高清观看| 757午夜福利合集在线观看| tocl精华| 亚洲自拍偷在线| 久久人妻av系列| 叶爱在线成人免费视频播放| 欧美日韩综合久久久久久 | 国产一区二区三区在线臀色熟女| 变态另类成人亚洲欧美熟女| 很黄的视频免费| 国产欧美日韩一区二区三| 欧美高清成人免费视频www| 最近最新免费中文字幕在线| 日本免费一区二区三区高清不卡| 国产成人欧美在线观看| 国产高清视频在线观看网站| 熟妇人妻久久中文字幕3abv| tocl精华| 久久精品国产自在天天线| 99精品在免费线老司机午夜| 狠狠狠狠99中文字幕| 12—13女人毛片做爰片一| 老司机福利观看| 中文字幕人妻丝袜一区二区| 757午夜福利合集在线观看| 麻豆成人午夜福利视频| avwww免费| 麻豆一二三区av精品| 日本 欧美在线| 国产精品影院久久| 中文字幕人妻熟人妻熟丝袜美 | 欧美日韩综合久久久久久 | 男女之事视频高清在线观看| 麻豆久久精品国产亚洲av| 特大巨黑吊av在线直播| 亚洲av成人av| 天天一区二区日本电影三级| 午夜精品久久久久久毛片777| 国产真人三级小视频在线观看| 成人一区二区视频在线观看| 夜夜爽天天搞| 国产精品影院久久| xxx96com| 国产高清有码在线观看视频| 91久久精品电影网| 欧美丝袜亚洲另类 | 国产精品久久久久久久久免 | 白带黄色成豆腐渣| 亚洲精华国产精华精| 日本黄色视频三级网站网址| 老鸭窝网址在线观看| 日韩中文字幕欧美一区二区| 好男人在线观看高清免费视频| 国产精品乱码一区二三区的特点| 久久久国产精品麻豆| 精品久久久久久久人妻蜜臀av| 51午夜福利影视在线观看| 日本黄色视频三级网站网址| 在线看三级毛片| 美女大奶头视频| 校园春色视频在线观看| 久久久久久久久大av| 国产激情欧美一区二区| 综合色av麻豆| 欧美精品啪啪一区二区三区| 国产亚洲精品久久久com| 每晚都被弄得嗷嗷叫到高潮| 99热6这里只有精品| 欧美bdsm另类| 国内精品久久久久精免费| 欧美日韩中文字幕国产精品一区二区三区| 亚洲欧美日韩高清在线视频| a级毛片a级免费在线| 搡老熟女国产l中国老女人| 国产亚洲精品综合一区在线观看| 黄色视频,在线免费观看| 久久久久国产精品人妻aⅴ院| 国产成人影院久久av| 亚洲专区中文字幕在线| 亚洲av电影不卡..在线观看| 久久久久久久久大av| 亚洲av成人精品一区久久| 有码 亚洲区| 老司机深夜福利视频在线观看| 丝袜美腿在线中文| 国产高清视频在线播放一区| 中出人妻视频一区二区| 午夜福利欧美成人| 国产综合懂色| 哪里可以看免费的av片| 色综合婷婷激情| 狠狠狠狠99中文字幕| 熟女电影av网| av片东京热男人的天堂| 精品国产三级普通话版| 久久久精品大字幕| 极品教师在线免费播放| 欧美最新免费一区二区三区 | 国产亚洲欧美在线一区二区| 特级一级黄色大片|