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

    基于伴隨方程復雜管網(wǎng)泄漏反向?qū)ぴ吹臄?shù)值模擬研究

    2022-09-06 08:08:08韓宇波洪文鵬
    東北電力大學學報 2022年1期
    關鍵詞:管網(wǎng)測點方程

    韓宇波,常 暢,宋 琪,洪文鵬

    (1.珠海橫琴能源發(fā)展有限公司,廣東 珠海 519015;2.東北電力大學能源與動力工程學院,吉林 吉林 132012)

    目前大多比較先進成熟的泄漏檢測定位方法和技術均是針對單一的長直管道,例如長輸油、輸天然氣管道[1-3].然而,對于拓撲結構復雜的多分支型管網(wǎng)來說,例如供熱管網(wǎng)、供水管網(wǎng),這些方法將不能直接被利用.復雜供水管網(wǎng)廣泛存在破損、泄漏等問題.給社會生產(chǎn)和人民生活造成巨大經(jīng)濟損失以及諸多不便,因此,為了及時的發(fā)現(xiàn),定位,并修補泄漏,消除水泄漏造成的潛在危險,本文針對供水管網(wǎng)泄漏檢測定位,提出了基于伴隨方程的反向?qū)ぴ捶椒?

    1994年,靳世久[4]提出基于負壓波的管道泄漏檢測定位方法.1996年,楊開林、郭宗周[5]提出了基于恒定流動模擬的靜態(tài)泄漏檢測法和基于水力瞬態(tài)模擬的瞬態(tài)泄漏檢測法.2001年,Verde和Mpesha[6-7]分別提出了最小非線性觀測器和頻率響應法來確定泄漏點的位置.2002年,Costa[8]運用參數(shù)估計法來對管網(wǎng)進行泄漏檢測和泄漏點定位.2003年,F(xiàn)errante[9]提出了基于小波的分析方法,根據(jù)上、下游傳感器的時間差進行泄漏點定位.2007年,劉暢等[10]利用多元統(tǒng)計學上的聚類分析方法,對管網(wǎng)故障模擬的數(shù)據(jù)進行聚類分析,并利用多元統(tǒng)計學上的判別分析方法,對泄漏情況進行故障定位.2010年,楊紅英等[11]研究了基于Duffing振子的管道泄漏檢測方法.2014年,呂玉坤等[12]對于單點泄漏的枝狀供熱管網(wǎng),引入插入虛擬節(jié)點的方法,來判斷管網(wǎng)是否泄漏與泄漏可能發(fā)生的位置.2015年,Wachla等[13]采用人工神經(jīng)模糊模型對供水管道泄漏進行定位.2015年,燕山大學黃滿義等[14]針對單一參數(shù)研究輸油管道泄漏,引入混沌特性分析及遞歸圖方法對管道泄漏檢測與定位方法進行分析.2016年,劉煒,劉宏昭[15]從時間序列角度出發(fā),提出一種基于模糊集理論的輸油管道泄漏檢測定位方法.

    Enting[16]將反問題用于大氣污染物擴散研究.Bagtzoglou[17]將反演方法用于地下水污染定位研究.郭少東等[18]采用計算流體力學的方法求解伴隨方程和MCMC方法的室內(nèi)污染源反演結果的研究.Zhang[19]將源項反演方法用在室內(nèi)污染源定位研究.

    Liu和Zhai[20-22]使用伴隨方法實現(xiàn)了室內(nèi)污染源的快速辨識,但他們并沒有將其作為最優(yōu)化方法使用,而是利用伴隨方法的特性,將在流場固定的基礎上實現(xiàn)了時間和對流項的取反,以反向方法的方式實現(xiàn)了污染源的辨識.

    總結以上國內(nèi)外研究現(xiàn)狀發(fā)現(xiàn),現(xiàn)有的管網(wǎng)漏損檢測定位方法并不能滿足日常管網(wǎng)漏損檢測的需求,在精度或者經(jīng)濟性上均存在一定的缺點.目前,基于流場求解的反向?qū)ぴ醇夹g能夠?qū)崿F(xiàn)準確地檢測和定位,在診斷空氣污染和水污染源頭的方向上取得了較好的效果,因此,本文將基于伴隨方程的反向?qū)ぴ醇夹g應用于供水管網(wǎng)中,以更快的實現(xiàn)復雜管網(wǎng)的泄漏檢測和定位.

    1 數(shù)學模型的建立

    1.1 單相流瞬態(tài)伴隨方程的推導

    建立管網(wǎng)流場瞬態(tài)方程,找出與流場瞬態(tài)壓力耦合的參數(shù),保留這些參數(shù),對流場瞬態(tài)方程進行簡化.

    不考慮摩擦、阻尼等因素的影響,一維半無限區(qū)間下壓力瞬態(tài)流動方程為

    (1)

    初始條件為

    當t=0時,p(x,t)=0,v(x,t)=0;

    邊界條件為

    當x=0時,p(x,t)=p0(t),v(x,t)=0;當x→+∞時,p(x,t)=0;

    公式中:p為壓力,Pa;v為液體流速,m3/s;a為壓力波速,m/s;ρ為液體密度,kg/m3.

    (1)以管網(wǎng)泄漏對流場形成的壓力波影響為依據(jù),選取適當敏感度函數(shù).所選取的壓力波敏感函數(shù)

    h(α,p)=p(x,t)δ(x-x*)δ(t-t*),

    (2)

    公式中:p(x,t)為壓力波傳播函數(shù);α為狄克拉函數(shù),α與α為泄漏或堵塞故障發(fā)生的位置與時間;α為敏感度參數(shù)(過程變量);p為壓力.

    (2)將敏感度函數(shù)對壓強求導,并引入伴隨算子,推導管網(wǎng)流場瞬態(tài)伴隨方程.利用高斯散度定理處理,設置邊界條件,得到流場瞬態(tài)伴隨方程.

    量化系統(tǒng)狀態(tài)量的度量目標函數(shù)

    (3)

    公式中:h(α,p)為系統(tǒng)狀態(tài)函數(shù);α為系統(tǒng)矢量參數(shù)(例如[v,D,a];p為壓力,積分區(qū)間是給定的空間-時間區(qū)間,系統(tǒng)關于矢量參數(shù)αk的臨界靈敏度可以通過對參數(shù)αk求微分得到

    (4)

    (5)

    初始條件為

    當t=0時,Φ(x,t)=0;當t=0時,λ(x,t)=0;

    邊界條件為

    當x→+∞時,φ(x,t)=0;

    當x=0時,λ(x,t)=0,φ(x,t)=0;當x→+∞時,λ(x,t)=0.

    (6)

    初始條件為

    φ*(x,T)=0,λ*(x,T)=0;

    邊界條件為

    φ*(0,t)=p0(t),λ*(0,t)=0.

    這樣伴隨方程(6)中未定義函數(shù)僅剩余h(α,p)項,它的形式還需要根據(jù)所要研究的具體工況條件確定,將方程(2)中選取的適當函數(shù)兩邊取偏導可得

    (7)

    將方程(7)帶入到伴隨方程(6)中,經(jīng)過整理得

    (8)

    初始條件為

    φ*(x,T)=0,λ*(x,T)=0;

    邊界條件為

    φ*(0,t)=p0(t),λ*(0,t)=0.

    考慮流動過程中摩擦阻力,最終得出管網(wǎng)流場瞬態(tài)伴隨方程

    (9)

    邊界條件為

    φ*(0,t)=p0(t),φ*(L,t)=0,λ*(0,t)=0,λ*(L,t)=0;

    初始條件為

    φ*(x,T)=0,λ*(x,T)=0.

    公式中:t為時間變量;x為距離變量;L為距離常數(shù);T為時間常數(shù).

    (3)將伴隨方程進行空間與時間的反向處理,應用拉普拉斯變換和傅里葉變換求解析解,得到反向?qū)ぴ茨P?

    做反向處理,τ=td-t,x為-x,可得反向伴隨方程

    (10)

    公式中:τ=td-t,td為檢測的時間;xd為檢測的位置.

    1.2 模型驗證

    傳統(tǒng)的正向求解方法是通過對每個可能的泄漏位置進行正向求解,判斷泄漏點的位置,但是這種方法因為需要多次求解從而變的低效.因此,本文利用伴隨理論和數(shù)學方法結合流體力學中的單相流體控制方程推導伴隨方程,建立逆向模型,如圖1所示.

    圖1 泄漏位置檢測原理圖

    從檢測點開始傳播的壓力變化曲線,如圖2所示.從圖2中可以看到明顯的突升變化過程,該曲線是通過在假定的壓力測定點出檢測到負壓波信號,并設置為xd=0,即為原點,得到反向傳播時間,τ=5 000 m/1 300 m/s,且泄漏點出的壓力變化為突升25 kPa.

    圖2 從檢測點開始傳播的壓力變化曲線圖3 檢測點開始傳播的壓力曲線

    假設在xd=0監(jiān)測點處開始檢測壓力波,如圖3所示,模擬了泄漏時間τ=5.36 s時,壓力波反向傳播變化曲線,壓力值從25 kPa突變到50 kPa,泄漏突變的位置在x=6 968 m處.

    綜上所述,基于壓力輸運控制方程,應用伴隨理論及一些數(shù)學方法推導出相應的管網(wǎng)瞬態(tài)伴隨方程,通過MATLAB模擬了檢測點處逆向壓力波變化規(guī)律,明確了由壓力輸運方程推導出來反向伴隨方程可以應用于管道泄漏檢測與定位.

    2 實驗結構及數(shù)值模擬計算

    2.1 實驗結構

    以實驗為基準驗證使用Open FOAM進行數(shù)值模擬的結果,實驗中對管道泄漏前后以及泄漏瞬間的壓力進行了測量.實驗系統(tǒng)圖如圖4所示,主要由一個實驗段,五個壓力變送器,兩個電磁流量計,兩個過濾器,一個水箱,一個水泵,一個數(shù)據(jù)采集儀和一個電腦組成,整體分為復雜管網(wǎng)微縮實驗系統(tǒng)、信號傳輸系統(tǒng)和數(shù)據(jù)處理系統(tǒng)三個基礎部分.

    圖4 實驗系統(tǒng)圖

    由于實際供水管網(wǎng)雷諾數(shù)大致范圍在5×104到4×105之間,且管網(wǎng)穩(wěn)定運行時管內(nèi)流態(tài)為阻力平方區(qū),結合模擬供水管網(wǎng)、運行參數(shù)及實驗室空間,根據(jù)相似理論進行了管網(wǎng)模型設計.選擇水為管網(wǎng)運行介質(zhì),為了滿足可視性,使用強度大、安全性高、安裝簡單等要求,選擇有機玻璃為管網(wǎng)主體管道,U-PVC為閥門和四通材料.

    管網(wǎng)主體為4×4復雜環(huán)狀管網(wǎng),如圖5所示.管材為內(nèi)徑12 mm外徑20 mm的有機玻璃,泄漏內(nèi)徑4 mm,外徑6 mm.管網(wǎng)各環(huán)管段尺寸相同,長均為0.85 m,寬均為0.38 m,管段中共有25個節(jié)點,其中5個節(jié)點為壓力測點.共設置20個泄漏點,其中10個位于管段節(jié)點,另外10個為管段中點,分別對節(jié)點和管段中泄漏進行研究.

    圖5 實驗臺壓力測點和泄漏點的分布

    2.2 數(shù)值模擬計算

    2.2.1 網(wǎng)格的劃分

    根據(jù)實驗段尺寸,建立模擬模型,忽略壁厚,管道直徑12 mm,泄漏孔直徑4 mm,整體采用三維六面體結構化網(wǎng)格,管道交叉處的與泄漏口處分別采用Y-Block和O-block,實驗段為4×4管網(wǎng),網(wǎng)格數(shù)達9 061 188,為了便于計算,將其簡化為2×2管網(wǎng),網(wǎng)格數(shù)減少到4 420 092,計算效率大大提高.

    圖6 4*4模型網(wǎng)格劃分圖7 2*2模型網(wǎng)格劃分

    本文選擇三維六面體結構化網(wǎng)格,如表1所示.表1做了網(wǎng)格無關性驗證,本次模擬采用的網(wǎng)格數(shù)為4 420 092.

    表1 網(wǎng)格無關性驗證

    由表1可知,以泄漏點1為例,當網(wǎng)格數(shù)達到400萬,監(jiān)測點壓力的壓降幾乎不變,為了提高管網(wǎng)泄漏模擬計算的精度和效率,對管道相貫處,相貫處的泄漏點,以及管道進出口等進行了網(wǎng)格加密處理.

    2.2.2 復雜管網(wǎng)泄漏檢測及相關性分析

    同一泄漏源產(chǎn)生的壓力波,其傳播到不同傳感器所需的時間不同,因此不同傳感器上所接收到的信號記錄也不會相同,但是當把兩個傳感器記錄其中之一給予一定的時移后,會發(fā)現(xiàn)此時兩個信號序列之間具有很大的相似性.根據(jù)壓力波傳到兩個監(jiān)測點的時間差,壓力波波速即可確定泄漏點的位置.

    為提高時延估計的精度,采用了希爾伯特變換和二次相關法相結合的新時延估計方法,其基本思路在于結合兩種分析方法,將其與相關函數(shù)的絕對值相減,相關峰值點被保留,其附近的點的相關值相應的減小,主峰值點被銳化,精度大大提高,在如上所述的二次相關法的基礎上,加入了希爾伯特變換,并結合了二次相關方法,原理如圖8所示.

    圖8 新時延估計法原理圖

    銳化了主峰值點后,利用互相關系數(shù)法判斷信號相似程度,其原理是兩列信號的相關系數(shù)在泄漏前通常在零值附近,出現(xiàn)泄漏后,信號的相關系數(shù)在零附近有所偏移,信號峰值點所對應的時間就是相關系數(shù)最大時所對應的時間點.假定兩列不同隨機信號s1(t)、s2(t),互相關函數(shù)為

    (11)

    公式中:RS1S2(τ)可以計算出s1(t)和s2(t)信號的相似程度,它是信號的相關函數(shù),其中,τ為時移.

    以泄漏孔1、泄漏孔2、泄漏孔3為例,壓力波的傳播速度1 483.2 m/s,設置時間步長為5e-7s,采用PISO算法,自定義求解器,管網(wǎng)內(nèi)流體流動穩(wěn)定到2 s,出現(xiàn)泄漏.

    簡化4×4管網(wǎng)之后的管網(wǎng)系統(tǒng)圖,如圖9所示.設置兩個壓力測點,測點1(左上方測點)和測點2(右下方測點),監(jiān)測泄漏前后的壓力波信號.模擬泄漏前和整個過程測點1和測點2通過新時延估計法和互相關分析法處理后得到的信號相關性圖,如圖10~圖12所示.相關系數(shù)的最大值在泄漏前通常在零值附近,出現(xiàn)泄漏后,信號相關系數(shù)的最大值在零附近有所偏移,通過相關性對比圖可知,管道泄漏的出現(xiàn),最大相關系數(shù)的值下降,我們通過對比管道泄漏前后的閾值,判斷管道是否出現(xiàn)了泄漏.

    圖9 管網(wǎng)布置示意圖

    圖10 不同測點,1點泄漏前后壓力及相關性對比圖

    圖11 不同測點,2點泄漏前后壓力及相關性對比圖

    圖12 不同測點,3點泄漏前后壓力及相關性對比圖

    從上圖10至圖12中可知,管網(wǎng)出現(xiàn)泄漏后,壓力瞬間增大,與實際實驗數(shù)據(jù)相反,這是由于模擬過程添加了逆向伴隨模型,大大加快了計算效率.管網(wǎng)流動穩(wěn)定后,2 s時出現(xiàn)泄漏,壓力波迅速向管道上下游進行傳播,在泄漏后的0.2 s內(nèi)壓力變化明顯,出現(xiàn)上升的趨勢,管道內(nèi)部的沿程壓力也升高,略微波動后穩(wěn)定于略微高于管道正常壓力的值.管道入口我們采用速度進口,入口速度一定即流量一定,泄漏后短時間內(nèi),泄漏點上游的流速下降后又恢復到了泄漏前的穩(wěn)定值,左上方測點距離泄漏孔較近,右下方測點較遠,其接收到壓力波信號時間大于左上方測點,壓力值大于左上方測點,總體趨勢相同,壓力變化范圍在10 kPa~28 kPa之間,泄漏后相關系數(shù)減小,最大相關系數(shù)對應的時間有延遲,符合相關性驗證.

    管道內(nèi)流體流動2 s穩(wěn)定后,1,2兩點同時發(fā)生泄漏的模擬圖,相對于一點泄漏,兩點泄漏的峰值差更小,壓力變化范圍也較小,但總體趨勢一致,符合逆向模型得出的結果,可以利用上述方法對管網(wǎng)實現(xiàn)在線監(jiān)測,如圖13所示.

    圖13 不同測點,1,2兩點同時泄漏前后壓力對比圖

    綜上所述,管道突然出現(xiàn)泄漏后,相比于泄漏前,管道首末兩端的壓力差更大,壓降作用也更明顯.加入了反向模型的管網(wǎng),泄漏后管網(wǎng)各個監(jiān)測點的壓力迅速上升后略有降低,之后維持在一個較高的值,比較不同泄漏位置的壓力值,泄漏位置靠近管道進口,下游壓力降低幅度越高,越靠近管道出口,上游壓力下降幅度越小,根據(jù)管道沿程各監(jiān)測點收集到的壓力值梯度,可知泄漏情況的發(fā)生.

    2.2.3 復雜管網(wǎng)泄漏位置的確定

    管網(wǎng)泄漏檢測及定位流程圖,如圖14所示.當管網(wǎng)出現(xiàn)泄漏時,我們首先通過上述方法進行泄漏檢測,確定可能發(fā)生泄漏的位置,判斷各部分泄漏量的多少,然后利用全局搜索方式,最終確定泄漏點.

    結合模擬結果,我們分析了兩列信號s1(t)和s2(t)的相關性(以1點泄漏采集為例),計算可知壓力波傳到兩個監(jiān)測點的時間差為△t=1.15×10-3s,根據(jù)已知的壓力波波速信息就能計算出泄漏孔距離兩個測點的距離,兩者距離差約五倍,由于管網(wǎng)線是同程的,可知可能出現(xiàn)泄漏的位置如下圖所示的A-D,根據(jù)圖15所示,測點1測得的壓力值較大,且有所延遲,說明泄漏點在測點2附近,泄漏位置則是泄漏點A或B.

    圖15 泄漏可能出現(xiàn)位置示意圖

    因此,為了識別泄漏位置,提出了一種搜索算法.為了排出其他泄漏點,我們采取了新的方案,該方案包括兩個特點:它包括環(huán)形網(wǎng)絡;測量點的排列可以符合試驗的要求.這個管道泄漏產(chǎn)生的負壓波可以到達壓力測量點通過不同的方式,我們可以利用這個特性實驗誤差分析.具體地:

    步驟1:搜索最接近泄漏位置的節(jié)點

    (12)

    公式中:tj、tk為在tj和tk時刻,在節(jié)點j和k檢測到泄漏瞬變;Si為較小的殘值表示泄漏發(fā)生在i節(jié)點的概率較高;τik為讓我們把節(jié)點i到k的最短旅行時間表示出來.

    假如泄漏在節(jié)點i發(fā)生,當i=1……N(N=管網(wǎng)中的節(jié)點數(shù)目),從而

    (tj-tk)-(τij-τik)=0.

    (13)

    然而,由于時間、測量和其他誤差,這個等式的左邊永遠不會為零.

    步驟2:沿著連接到最近節(jié)點的管道段搜索泄漏位置.

    在這一步驟中,沿著步驟1確定的與節(jié)點ni相連的管道部分(即沿著圖模型的邊緣)放置一組新的虛擬節(jié)點.

    第一步是對網(wǎng)絡中所有節(jié)點進行粗略的全局搜索;

    第二步在最近的節(jié)點估計附近執(zhí)行本地搜索,以確定沿著管道段的最可能泄漏位置.

    此外,管道泄漏產(chǎn)生的壓力波傳播將通過不同的路徑到達測量點,我們可以以A點泄漏時,M1和M2點為例.

    圖16 A點泄漏時,傳遞到測點M2的路徑

    圖17 A點泄漏時,傳遞到測點M1的路徑

    如圖16和17所示,假設A點泄漏,則最符合本次模擬的路線是圖16中的(a)和圖17中的(b).

    表2為根據(jù)圖14所示的方法計算出的的每個節(jié)點的Ei值,由表中可以看出節(jié)點3的Ei值最小,節(jié)點4的Ei值次之,故離泄漏位置最近的節(jié)點為3,并且泄漏位置介于節(jié)點3和4之間,泄漏點為點A.根據(jù)△t=1.15×10-3s,已知壓力波的波速,我們可以求出泄漏點最終確定泄漏位置414 mm,實際泄漏點在431 mm處(泄漏孔1),在誤差允許的范圍內(nèi),所以該方法可行.

    表2 不同節(jié)點Ei計算結果

    同時,模擬了不同泄漏孔大小和不同進口壓力下的5個泄漏點的定位情況,分別進行了30次和20次模擬計算,為了便于分析,取以測點2為起點,傳播路徑最短的一組數(shù)據(jù),現(xiàn)以泄漏孔3和泄漏孔5為例,定位情況如表3、表4所示,誤差均在允許的范圍內(nèi),說明該方法可行.

    表3 以泄漏孔3為例的不同泄漏孔的定位情況

    表4 以泄漏孔5為例的不同進口壓力的定位情況

    2.2.4 不同方法及算法對比

    2.2.4.1 不同方法對比

    由表5可知,近些年現(xiàn)有的檢測方法如ITA(Inverse Transient Analysis),Model Falsification & EPANET,PSO-SVM(Particle Swarm Optimization-Support Vector Machine)等方法的定位能力較高,均在90%以上.本文提到的于反向?qū)ぴ捶ǘㄎ坏臏蚀_度在92%~96%之間,在復雜管網(wǎng)的泄漏定位中有一定的應用價值.

    表5 不同泄漏檢測及定位方法對比

    2.2.4.2 算法對比

    如圖18所示,利用表1中經(jīng)過網(wǎng)格無關性驗證的網(wǎng)格數(shù)為442 009的管網(wǎng)模型,在相同的邊界條件和參數(shù)的設置下,對比了在50 000、100 000、150 000、200 000計算步數(shù)下,SIMPLE、PISO及加入伴隨方程的PISO算法所需的迭代時間,隨著計算步數(shù)的增多,三種方式的計算時間均增長.結果表明,SIMPLE算法應用于管網(wǎng)泄漏模擬中,計算速度最慢,伴隨方程的PISO算法計算速度最快.

    圖18 不同計算方法計算速度對比

    我們利用了進行了網(wǎng)格無關性驗證的管網(wǎng)模型了,進行了一段時間的計算,對比了SIMPLE算法,PISO算法和本文加入了伴隨方程的PISO算法,在相同的時間步長下,本文的計算速度明顯加快.

    2.3 實驗驗證

    在此部分通過對不同漏點進行實驗,分析管網(wǎng)中的壓力變化規(guī)律,將實驗數(shù)據(jù)同模擬結果進行對比分析.以泄漏點1為例,模擬和實驗壓力變化規(guī)律同模擬結果符合度較高.

    如圖19所示為以泄漏孔1為例的泄漏孔處的壓力變化曲線,圖中藍色線條代表壓力變化的總體趨勢,黑色線條分別為試驗和模擬條件下壓力變化趨勢.從壓力變化曲線中可以發(fā)現(xiàn)兩者變化規(guī)律相同,實驗數(shù)據(jù)在泄漏前壓力值保持在一個高位水平,泄漏瞬間壓力值突降,在到達最底端后,恢復值一個較低的壓力值,并逐漸穩(wěn)定.而模擬數(shù)據(jù)與實驗變化規(guī)律相反,這是由于添加了逆向模型的緣故.對比實驗和模擬數(shù)據(jù),模擬過程阻力小于實驗,無噪聲的影響,因此模擬壓力的最大值高于實驗,而泄漏前后的壓降大小基本相同.總的來說,從橫軸觀察泄漏瞬間對應的時間差值可以發(fā)現(xiàn),泄漏是在一瞬間發(fā)生的壓力變化時間極端,這與模擬的結果相吻合,并且正向壓力傳播曲線與模擬所得反向壓力傳播曲線呈現(xiàn)的壓力變化趨勢相吻合.

    圖19 實驗模擬對比圖

    接下來為了進一步驗證方法的合理性,對與模擬工況對應的多組實驗數(shù)據(jù)做了處理,將50組實驗和模擬的工況進行了計算,對比了各自的定位準確度,如圖20所示.

    圖20 不同泄漏孔直徑下實驗模擬定位準確度對比圖圖21 不同泄漏點管道泄漏定位準確度變化曲線圖

    不同泄漏孔直徑下實驗模擬定位準確度對比圖,如圖20所示.可知隨著泄漏孔直徑增大實驗和模擬點位準確度均增大,但當泄漏孔直徑達到12 mm時,再增大泄漏孔,由于泄漏孔徑很大,泄漏產(chǎn)生的壓力波經(jīng)過兩個測點時間很短,時間差計算不準確,進而導致定位位置計算不準確,準確度有所下降.總的來說,實驗結果與模擬結果相近,定位準確度在2%左右波動.

    不同泄漏點管道泄漏定位準確度變化曲線圖,如圖21所示.實驗管網(wǎng)的形狀以及五個泄漏點的位置一定,泄漏點1和泄漏點5距離兩個壓力測點路徑長度一致,定位準確度基本一致,泄漏點2距離測1距離相對較遠,誤差比1、5大,泄漏點3的定位準確度最低,3點泄漏信號傳遞到兩個壓力測點,無論怎樣傳遞都要兩次經(jīng)過四通,并且距離兩個測點的長度較長,定位準確度比較低.

    總體來說,實驗結果與模擬結果基本吻合,均在誤差允許的范圍內(nèi),該方法可行.另外本研究是基于實際管網(wǎng)進行的簡化,實際地下供熱管網(wǎng)沒有四通,多用二個臨近三通來代替,壓力波的反射情況完全不同,給復雜實際管網(wǎng)檢測中帶來一定的干擾,但是這種干擾只是增加了一部分計算過程中的局部阻力,壓力波傳遞情況稍有改變,對整體壓力變化趨勢以及泄漏檢測和定位精度影響較小.

    3 結 論

    本文的主要結合反向?qū)ぴ丛砑鞍殡S方法,推導瞬態(tài)流場反向伴隨方程.借助MATLAB軟件和開源有限元軟件OpenFOAM,構建單管和復雜管網(wǎng)流體流動數(shù)值模擬模型,在OpenFOAM自定義求解器,將反向?qū)ぴ磫栴}轉化為伴隨方程的求解問題.根據(jù)得到的圖形分析不同測點的壓力等參數(shù)的變化,結合新時延估計方法和互相關方法,對管網(wǎng)的運行狀態(tài)進行分析檢測,完成對伴隨方程的驗證和求解,確定泄漏源的可能泄漏位置,最終利用全局搜索方法完成了泄漏源的定位.得出以下結論:

    (1)基于正向壓力輸運模型,引入目標函數(shù),將伴隨理論與靈敏度分析方法相結合推導了正向伴隨方程.使敏感度函數(shù)對壓強求導,并引入伴隨算子,推導了瞬態(tài)流場反向伴隨方程,求得了其解析解,并用MATLAB軟件,完成了伴隨模型的驗證,結果證明反向?qū)ぴ捶椒捎糜诠芫W(wǎng)泄漏檢測.

    (2)在OpenFOAM中完成了復雜管網(wǎng)流體流動數(shù)值模擬模型的構建和伴隨方程的求解,并監(jiān)測了泄漏孔泄漏瞬間壓力變化,運用新時延估計方法銳化了峰值點,最后用互相關法進行了相關性分析,發(fā)現(xiàn)管道泄漏后,信號相關系數(shù)的最大值在零附近有所偏移,通過相關性對比圖可知,管道泄漏的出現(xiàn),最大相關系數(shù)的值下降,以泄漏點1為例證明了方法的可行性.

    (3)比較了不同泄漏孔大小及不同進口壓力下,模擬泄漏定位的準確度,結果表明,運用基于反向?qū)ぴ捶ǘㄎ坏臏蚀_度在92%~96%之間,說明此方法可以應用在復雜管網(wǎng)泄漏檢測及定位中.

    (4)本文在不同計算步數(shù)時,對比了單一的SIMPLE、PISO及加入伴隨方程的PISO算法所需的迭代時間,結果表明,結合伴隨方程的PISO算法的計算速度要明顯優(yōu)于單一的SIMPLE和PISO算法,以此證明了伴隨方法的優(yōu)越性.

    在本文中,基于復雜管網(wǎng)泄漏檢測及定位模型,未考慮同樣引起泄漏的堵塞、變形及污染等問題,應定義相應的求解器,設置相應的參數(shù)及邊界條件,盡可能將定位方法轉化成OpenFOAM語言,便可通過模擬體系,直接得到不同情況下的診斷結果.此外,實驗系統(tǒng)存在一定的不穩(wěn)定性和局限性,與實際復雜管網(wǎng)還有一定差距,只是在宏觀上完成了伴隨方程的求解和驗證,需進一步結合實際工況強化實驗和模擬.

    猜你喜歡
    管網(wǎng)測點方程
    液壓支架整機靜強度試驗及等效應力分析
    方程的再認識
    方程(組)的由來
    基于CATIA的汽車測點批量開發(fā)的研究與應用
    圓的方程
    管網(wǎng)獨立是妥協(xié)還是改革
    能源(2018年8期)2018-09-21 07:57:20
    從管網(wǎng)獨立看國企改革
    能源(2018年8期)2018-09-21 07:57:20
    管網(wǎng)改革虛實
    能源(2018年8期)2018-09-21 07:57:18
    織起一張共管網(wǎng)
    中國公路(2017年8期)2017-07-21 14:26:20
    拱壩結構損傷的多測點R/S分析
    国产精品久久电影中文字幕| 黄色成人免费大全| 久久狼人影院| 交换朋友夫妻互换小说| 色综合欧美亚洲国产小说| 在线国产一区二区在线| 99国产极品粉嫩在线观看| 在线天堂中文资源库| 国产精华一区二区三区| 久久人妻av系列| 在线十欧美十亚洲十日本专区| 亚洲精品一区av在线观看| 亚洲精品国产精品久久久不卡| a级毛片在线看网站| 91麻豆av在线| 亚洲一区中文字幕在线| 日韩 欧美 亚洲 中文字幕| 女性生殖器流出的白浆| 免费在线观看亚洲国产| 日韩免费高清中文字幕av| 欧美+亚洲+日韩+国产| 亚洲av成人av| 黄色片一级片一级黄色片| 日韩大尺度精品在线看网址 | 男女下面进入的视频免费午夜 | 涩涩av久久男人的天堂| 欧美日韩乱码在线| 欧美另类亚洲清纯唯美| 国产精品1区2区在线观看.| 少妇粗大呻吟视频| 在线播放国产精品三级| 日韩免费高清中文字幕av| 精品国产亚洲在线| 亚洲国产精品999在线| 91麻豆精品激情在线观看国产 | 国产1区2区3区精品| 女性被躁到高潮视频| 亚洲精品中文字幕在线视频| 亚洲精品一区av在线观看| 国产精品免费一区二区三区在线| 日韩三级视频一区二区三区| 欧美日韩乱码在线| 亚洲五月婷婷丁香| 性少妇av在线| 国内久久婷婷六月综合欲色啪| 窝窝影院91人妻| 不卡av一区二区三区| 中文亚洲av片在线观看爽| 新久久久久国产一级毛片| 久久精品国产清高在天天线| 俄罗斯特黄特色一大片| 在线观看免费日韩欧美大片| 亚洲欧美精品综合久久99| 亚洲在线自拍视频| 超碰成人久久| 精品国产美女av久久久久小说| 一级,二级,三级黄色视频| 多毛熟女@视频| 国产伦一二天堂av在线观看| 少妇 在线观看| 日韩欧美国产一区二区入口| 亚洲三区欧美一区| 久久久久亚洲av毛片大全| 丰满的人妻完整版| 国产成人精品在线电影| 中国美女看黄片| 性少妇av在线| 精品第一国产精品| 超碰成人久久| 一进一出抽搐gif免费好疼 | 免费久久久久久久精品成人欧美视频| 日日夜夜操网爽| 色老头精品视频在线观看| 伊人久久大香线蕉亚洲五| 两个人免费观看高清视频| 亚洲av日韩精品久久久久久密| 麻豆av在线久日| av片东京热男人的天堂| 最新美女视频免费是黄的| 亚洲一区二区三区不卡视频| 老司机靠b影院| av在线播放免费不卡| 美女高潮到喷水免费观看| 午夜福利在线观看吧| 久久草成人影院| 中文字幕高清在线视频| 色在线成人网| 日韩人妻精品一区2区三区| 国产精品 欧美亚洲| 亚洲专区国产一区二区| 亚洲视频免费观看视频| 欧美日韩一级在线毛片| 国产1区2区3区精品| 亚洲伊人色综图| 亚洲九九香蕉| 久久伊人香网站| 亚洲精品成人av观看孕妇| 免费在线观看完整版高清| 制服诱惑二区| 久久久久久久久中文| av天堂在线播放| 精品乱码久久久久久99久播| av超薄肉色丝袜交足视频| 亚洲国产中文字幕在线视频| 国产亚洲精品综合一区在线观看 | 亚洲avbb在线观看| 久久精品91蜜桃| 中文字幕人妻丝袜一区二区| 桃色一区二区三区在线观看| 精品福利观看| a级片在线免费高清观看视频| 青草久久国产| 在线观看舔阴道视频| 亚洲国产毛片av蜜桃av| 十八禁人妻一区二区| 久久人人爽av亚洲精品天堂| 亚洲精品在线观看二区| 色精品久久人妻99蜜桃| 黄色毛片三级朝国网站| 精品国产一区二区三区四区第35| 中亚洲国语对白在线视频| 久久影院123| 午夜视频精品福利| 日韩视频一区二区在线观看| 欧美中文综合在线视频| 亚洲专区中文字幕在线| 久久精品国产亚洲av高清一级| www.精华液| 精品一区二区三区视频在线观看免费 | 亚洲男人天堂网一区| 欧美日韩一级在线毛片| 日韩精品中文字幕看吧| 日本黄色日本黄色录像| 99久久综合精品五月天人人| 亚洲精品一卡2卡三卡4卡5卡| 亚洲第一青青草原| 黄色片一级片一级黄色片| 一本大道久久a久久精品| 男女下面进入的视频免费午夜 | 99精品欧美一区二区三区四区| 99精品久久久久人妻精品| 精品卡一卡二卡四卡免费| 欧美在线黄色| 亚洲一卡2卡3卡4卡5卡精品中文| 激情在线观看视频在线高清| 成年人免费黄色播放视频| av片东京热男人的天堂| 亚洲三区欧美一区| 国产一区二区三区在线臀色熟女 | 免费在线观看完整版高清| 日韩有码中文字幕| 热99国产精品久久久久久7| 免费av毛片视频| 亚洲熟妇中文字幕五十中出 | 啦啦啦免费观看视频1| 天堂中文最新版在线下载| 欧洲精品卡2卡3卡4卡5卡区| 18禁观看日本| 在线观看午夜福利视频| 午夜福利影视在线免费观看| 另类亚洲欧美激情| 日韩人妻精品一区2区三区| 国产成人欧美在线观看| 91精品国产国语对白视频| 精品国产超薄肉色丝袜足j| av视频免费观看在线观看| 亚洲av熟女| 成年女人毛片免费观看观看9| 一区二区三区国产精品乱码| 国产无遮挡羞羞视频在线观看| 精品国产美女av久久久久小说| 欧美在线一区亚洲| 最新在线观看一区二区三区| 精品高清国产在线一区| 18禁观看日本| 国产一区二区在线av高清观看| 欧美人与性动交α欧美精品济南到| 黄色怎么调成土黄色| 黄色 视频免费看| 久久精品影院6| 18禁黄网站禁片午夜丰满| 黄片播放在线免费| 丰满迷人的少妇在线观看| 国产精华一区二区三区| 大码成人一级视频| 久久国产亚洲av麻豆专区| 丝袜在线中文字幕| 欧美老熟妇乱子伦牲交| 免费人成视频x8x8入口观看| 精品欧美一区二区三区在线| 两性午夜刺激爽爽歪歪视频在线观看 | 精品久久久精品久久久| 91字幕亚洲| 国产一区二区三区综合在线观看| 波多野结衣一区麻豆| 中文字幕av电影在线播放| 少妇粗大呻吟视频| 国产欧美日韩精品亚洲av| 精品免费久久久久久久清纯| 91成人精品电影| 国产91精品成人一区二区三区| 99国产综合亚洲精品| 免费在线观看亚洲国产| 真人做人爱边吃奶动态| 麻豆久久精品国产亚洲av | 亚洲一区高清亚洲精品| 国产亚洲精品一区二区www| 欧美日本亚洲视频在线播放| 极品教师在线免费播放| 精品久久久久久电影网| 91在线观看av| 免费看a级黄色片| 极品教师在线免费播放| 成年女人毛片免费观看观看9| 欧美激情久久久久久爽电影 | 18禁国产床啪视频网站| 久久精品国产99精品国产亚洲性色 | 日韩免费高清中文字幕av| 免费在线观看完整版高清| 亚洲成国产人片在线观看| 超碰97精品在线观看| 少妇裸体淫交视频免费看高清 | 日本 av在线| 免费人成视频x8x8入口观看| 最新美女视频免费是黄的| 色在线成人网| 亚洲一区二区三区欧美精品| 久久午夜综合久久蜜桃| 亚洲成人免费av在线播放| 午夜精品在线福利| 亚洲五月色婷婷综合| 天天添夜夜摸| 精品国内亚洲2022精品成人| 久久天堂一区二区三区四区| 欧美另类亚洲清纯唯美| 亚洲自拍偷在线| tocl精华| 亚洲人成网站在线播放欧美日韩| 性色av乱码一区二区三区2| 午夜免费鲁丝| a级片在线免费高清观看视频| 亚洲成国产人片在线观看| 老司机靠b影院| 免费人成视频x8x8入口观看| 亚洲自偷自拍图片 自拍| www.999成人在线观看| 成熟少妇高潮喷水视频| 亚洲美女黄片视频| 无人区码免费观看不卡| 亚洲精品中文字幕在线视频| 十八禁人妻一区二区| 超碰成人久久| 久久中文字幕一级| 一区二区三区国产精品乱码| 成年人免费黄色播放视频| 色尼玛亚洲综合影院| 99国产极品粉嫩在线观看| 成年人黄色毛片网站| 欧美激情久久久久久爽电影 | 十分钟在线观看高清视频www| 丝袜在线中文字幕| av超薄肉色丝袜交足视频| 欧美日韩国产mv在线观看视频| 免费在线观看完整版高清| 人妻久久中文字幕网| 日韩精品青青久久久久久| 十八禁人妻一区二区| 在线十欧美十亚洲十日本专区| 色综合站精品国产| 亚洲精品久久午夜乱码| 国产精品免费一区二区三区在线| 国产av一区在线观看免费| 色在线成人网| 桃色一区二区三区在线观看| 成人av一区二区三区在线看| 久久性视频一级片| 国产乱人伦免费视频| 亚洲九九香蕉| 真人做人爱边吃奶动态| 欧美激情极品国产一区二区三区| 两个人看的免费小视频| 黑人巨大精品欧美一区二区蜜桃| 最近最新中文字幕大全免费视频| 免费在线观看亚洲国产| 亚洲五月婷婷丁香| 真人做人爱边吃奶动态| 欧美午夜高清在线| 久久精品国产亚洲av高清一级| 亚洲色图综合在线观看| 神马国产精品三级电影在线观看 | 日本免费一区二区三区高清不卡 | 一级片免费观看大全| 人人妻人人添人人爽欧美一区卜| 999久久久国产精品视频| 国产有黄有色有爽视频| 久久久国产一区二区| 高清欧美精品videossex| 久久精品国产亚洲av香蕉五月| 18禁黄网站禁片午夜丰满| 成人永久免费在线观看视频| 日韩欧美一区视频在线观看| 嫁个100分男人电影在线观看| 亚洲中文日韩欧美视频| 亚洲美女黄片视频| 日韩av在线大香蕉| 久久久久精品国产欧美久久久| 悠悠久久av| 大陆偷拍与自拍| 精品午夜福利视频在线观看一区| 不卡一级毛片| 亚洲七黄色美女视频| 精品久久久久久久久久免费视频 | 久久人妻熟女aⅴ| 欧美中文综合在线视频| 身体一侧抽搐| 俄罗斯特黄特色一大片| 国产成人啪精品午夜网站| 老汉色∧v一级毛片| 俄罗斯特黄特色一大片| 国产片内射在线| 大型黄色视频在线免费观看| 精品免费久久久久久久清纯| 国产日韩一区二区三区精品不卡| 久久欧美精品欧美久久欧美| 久久亚洲真实| 黄色毛片三级朝国网站| 伊人久久大香线蕉亚洲五| 欧美+亚洲+日韩+国产| 亚洲国产欧美一区二区综合| 精品国产国语对白av| 久热这里只有精品99| 欧美日韩福利视频一区二区| 美女国产高潮福利片在线看| 国产精品二区激情视频| 女人被狂操c到高潮| 精品少妇一区二区三区视频日本电影| 国产精品自产拍在线观看55亚洲| 丰满人妻熟妇乱又伦精品不卡| 乱人伦中国视频| 国产成人影院久久av| 老司机深夜福利视频在线观看| 女生性感内裤真人,穿戴方法视频| 亚洲男人的天堂狠狠| 三上悠亚av全集在线观看| av超薄肉色丝袜交足视频| 色综合婷婷激情| 99国产精品免费福利视频| 老司机午夜十八禁免费视频| 水蜜桃什么品种好| 国产午夜精品久久久久久| 婷婷丁香在线五月| 国产精品影院久久| 桃红色精品国产亚洲av| 美国免费a级毛片| 美女高潮喷水抽搐中文字幕| 大码成人一级视频| 亚洲精品在线美女| 国产精品日韩av在线免费观看 | 99热国产这里只有精品6| 美女国产高潮福利片在线看| 国产黄a三级三级三级人| 狂野欧美激情性xxxx| 老汉色∧v一级毛片| 欧美成人性av电影在线观看| 国产99久久九九免费精品| 精品一区二区三卡| 久久香蕉激情| 日韩精品中文字幕看吧| 国产亚洲欧美在线一区二区| 熟女少妇亚洲综合色aaa.| 91字幕亚洲| 亚洲人成电影观看| 日韩精品中文字幕看吧| 五月开心婷婷网| 日韩欧美在线二视频| 麻豆成人av在线观看| 久久久久久大精品| 久久久久久久久中文| 国产欧美日韩精品亚洲av| 人妻久久中文字幕网| 亚洲少妇的诱惑av| 欧美成人免费av一区二区三区| 久久国产精品男人的天堂亚洲| 自拍欧美九色日韩亚洲蝌蚪91| 巨乳人妻的诱惑在线观看| 成熟少妇高潮喷水视频| 他把我摸到了高潮在线观看| 黄片播放在线免费| 久久久久久久午夜电影 | 久久99一区二区三区| 日韩精品中文字幕看吧| 在线观看午夜福利视频| 黄网站色视频无遮挡免费观看| 一区二区三区精品91| 欧美日韩黄片免| 久9热在线精品视频| 亚洲全国av大片| 免费在线观看亚洲国产| 久久国产精品男人的天堂亚洲| 欧美另类亚洲清纯唯美| 18禁国产床啪视频网站| 国产激情欧美一区二区| 天天添夜夜摸| 女性被躁到高潮视频| 亚洲国产欧美一区二区综合| 精品国产亚洲在线| 女警被强在线播放| 欧美一区二区精品小视频在线| 人人妻人人澡人人看| 神马国产精品三级电影在线观看 | 宅男免费午夜| 啦啦啦在线免费观看视频4| 国产一区二区激情短视频| 久久久久久久久中文| 国产有黄有色有爽视频| 午夜激情av网站| 又大又爽又粗| 黄色 视频免费看| 热99re8久久精品国产| 婷婷精品国产亚洲av在线| 国产成人欧美在线观看| 国产免费av片在线观看野外av| 丰满迷人的少妇在线观看| 国产一区二区三区综合在线观看| 黑人操中国人逼视频| 一夜夜www| 老司机在亚洲福利影院| 两性夫妻黄色片| 麻豆久久精品国产亚洲av | 精品一品国产午夜福利视频| 久久影院123| 久久久水蜜桃国产精品网| 99久久99久久久精品蜜桃| 精品免费久久久久久久清纯| 亚洲精品成人av观看孕妇| 欧美色视频一区免费| 国产深夜福利视频在线观看| www国产在线视频色| 久久精品aⅴ一区二区三区四区| 欧美大码av| 在线天堂中文资源库| 9191精品国产免费久久| 亚洲一码二码三码区别大吗| 日本黄色视频三级网站网址| 国产精品秋霞免费鲁丝片| 天天躁狠狠躁夜夜躁狠狠躁| 国产精品久久电影中文字幕| 国产xxxxx性猛交| 日韩欧美在线二视频| 人人妻人人爽人人添夜夜欢视频| 色综合婷婷激情| 亚洲欧美一区二区三区黑人| 精品少妇一区二区三区视频日本电影| 亚洲人成电影免费在线| 美女大奶头视频| 黄色丝袜av网址大全| 亚洲av成人av| 一本大道久久a久久精品| 最近最新免费中文字幕在线| 亚洲一区二区三区色噜噜 | 亚洲成人久久性| 免费不卡黄色视频| 啦啦啦免费观看视频1| 亚洲avbb在线观看| 免费在线观看完整版高清| 99国产精品99久久久久| 日本免费a在线| 少妇 在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 久久久久亚洲av毛片大全| 久久久国产一区二区| 日韩欧美三级三区| 美女福利国产在线| 久久精品国产亚洲av高清一级| 国产视频一区二区在线看| av福利片在线| 在线视频色国产色| 国产高清videossex| 国产一区二区在线av高清观看| 中文字幕精品免费在线观看视频| 老鸭窝网址在线观看| 妹子高潮喷水视频| 极品人妻少妇av视频| 国产又色又爽无遮挡免费看| 成人18禁高潮啪啪吃奶动态图| 国产精品日韩av在线免费观看 | 国产高清激情床上av| 人人澡人人妻人| 99精国产麻豆久久婷婷| 亚洲欧美精品综合一区二区三区| 如日韩欧美国产精品一区二区三区| 神马国产精品三级电影在线观看 | 黄色视频不卡| 免费在线观看视频国产中文字幕亚洲| 最近最新免费中文字幕在线| 久久久久久人人人人人| 国产单亲对白刺激| 久久中文看片网| 亚洲国产欧美网| 日韩人妻精品一区2区三区| 国产高清视频在线播放一区| 中出人妻视频一区二区| 欧美黄色片欧美黄色片| 亚洲av片天天在线观看| 一级片免费观看大全| 老熟妇仑乱视频hdxx| 精品久久蜜臀av无| 亚洲av美国av| 午夜福利影视在线免费观看| 这个男人来自地球电影免费观看| 欧美+亚洲+日韩+国产| 国产激情久久老熟女| 在线天堂中文资源库| 色婷婷久久久亚洲欧美| 国产av又大| 日韩成人在线观看一区二区三区| av天堂久久9| 亚洲狠狠婷婷综合久久图片| 亚洲欧美精品综合一区二区三区| 国产精品美女特级片免费视频播放器 | 男女做爰动态图高潮gif福利片 | 欧美日韩瑟瑟在线播放| 国产av又大| 天天躁狠狠躁夜夜躁狠狠躁| 两个人看的免费小视频| 丝袜美腿诱惑在线| 国产精品亚洲一级av第二区| 日韩一卡2卡3卡4卡2021年| 国产av一区二区精品久久| 国产成人免费无遮挡视频| 亚洲黑人精品在线| 纯流量卡能插随身wifi吗| 人妻久久中文字幕网| 日韩三级视频一区二区三区| 国产精品国产高清国产av| 99在线视频只有这里精品首页| 搡老岳熟女国产| 国产成人av激情在线播放| 日本免费a在线| 不卡一级毛片| 99久久综合精品五月天人人| 久久精品影院6| 欧美激情久久久久久爽电影 | 精品高清国产在线一区| 中亚洲国语对白在线视频| 99久久人妻综合| 男女高潮啪啪啪动态图| 午夜福利在线免费观看网站| 一二三四在线观看免费中文在| 日本一区二区免费在线视频| 国产精品香港三级国产av潘金莲| 日日摸夜夜添夜夜添小说| 老鸭窝网址在线观看| 黑人巨大精品欧美一区二区蜜桃| 国产片内射在线| 俄罗斯特黄特色一大片| 在线观看www视频免费| 首页视频小说图片口味搜索| 午夜福利,免费看| 不卡av一区二区三区| 别揉我奶头~嗯~啊~动态视频| 一区福利在线观看| 免费一级毛片在线播放高清视频 | 搡老乐熟女国产| 亚洲性夜色夜夜综合| 久久精品91蜜桃| 亚洲第一欧美日韩一区二区三区| 国产精品一区二区三区四区久久 | 免费av中文字幕在线| 嫩草影院精品99| 久久精品91蜜桃| 欧美成人免费av一区二区三区| 国产午夜精品久久久久久| 国产精品1区2区在线观看.| 99久久精品国产亚洲精品| 亚洲,欧美精品.| 免费日韩欧美在线观看| 国产欧美日韩综合在线一区二区| 婷婷六月久久综合丁香| 精品久久久精品久久久| 亚洲色图 男人天堂 中文字幕| 亚洲国产看品久久| 亚洲中文日韩欧美视频| 一区二区三区精品91| 国产xxxxx性猛交| 美女福利国产在线| 黄片大片在线免费观看| 欧美在线黄色| 国产成人精品久久二区二区免费| 国产成人影院久久av| 国产精品久久久av美女十八| 免费一级毛片在线播放高清视频 | 久久久久久久久免费视频了| 男女下面插进去视频免费观看| 精品电影一区二区在线| 国产aⅴ精品一区二区三区波| 亚洲精品中文字幕一二三四区| 91字幕亚洲| 久久热在线av| 久久99一区二区三区| 神马国产精品三级电影在线观看 | 免费女性裸体啪啪无遮挡网站| 老熟妇仑乱视频hdxx| av天堂久久9| 在线看a的网站| 乱人伦中国视频| 人妻丰满熟妇av一区二区三区| 日韩人妻精品一区2区三区| 久久久国产欧美日韩av| 国产精品电影一区二区三区| 最近最新免费中文字幕在线| 精品久久久久久,| 看片在线看免费视频| 久久精品亚洲av国产电影网| 国产xxxxx性猛交|