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

彈性/黏彈性混合管道瞬變流參數(shù)校核方法

2023-02-25 13:46:16孫強王祉皓伍悅濱徐瑩
科學(xué)技術(shù)與工程 2023年2期
關(guān)鍵詞:實驗

孫強, 王祉皓, 伍悅濱, 徐瑩

(1.東北林業(yè)大學(xué)土木工程學(xué)院人工環(huán)境與能源應(yīng)用研究所, 哈爾濱 150040; 2.哈爾濱工業(yè)大學(xué)建筑學(xué)院, 哈爾濱 150090; 3.哈爾濱工業(yè)大學(xué)寒地建筑科學(xué)與工程研究中心, 哈爾濱 150090; 4.哈爾濱商業(yè)大學(xué)能源與建筑工程學(xué)院, 哈爾濱 150028)

有壓管道經(jīng)常會因突然停泵或關(guān)閥引起的壓力突增或驟降的水力現(xiàn)象,可能導(dǎo)致設(shè)備損壞或管道炸裂,這種管內(nèi)流體的瞬變流動現(xiàn)象,在工程中稱之為水錘現(xiàn)象[1-2]。為保證輸水系統(tǒng)安全運行,對管道瞬變流動進(jìn)行分析和控制是學(xué)者們重點研究的方向[3-4]。

近些年,隨著高分子材料技術(shù)的迅猛發(fā)展,高分子聚合物管道因其價格低、抗腐蝕性高、節(jié)能環(huán)保等特點,廣泛應(yīng)用于供水管道系統(tǒng)中。現(xiàn)如今,塑料管道常用于替換傳統(tǒng)金屬管道,使得現(xiàn)階段管道系統(tǒng)中經(jīng)常出現(xiàn)金屬管道與塑料管道混合連接的情況。由于塑料管道的力學(xué)特性不同于傳統(tǒng)金屬管道的彈性力學(xué)特性,其同時表現(xiàn)出彈性和黏性兩種力學(xué)行為,因而也被稱為黏彈性管道。黏彈性管道在瞬變應(yīng)力的激勵下,所產(chǎn)生的應(yīng)變可以表示為彈性應(yīng)變和延遲應(yīng)變的疊加。這種黏彈力學(xué)行為還會影響瞬變流動中壓力波動的峰值、相位和衰減。因黏彈性管道力學(xué)行為的復(fù)雜性,傳統(tǒng)彈性管道瞬變流的數(shù)學(xué)模型已不再適用于黏彈性管道瞬變流的數(shù)值模擬。近些年來,學(xué)者開始對黏彈性管道瞬變流壓力波動的影響機(jī)制進(jìn)行了理論和實驗研究[5-6],并先后提出并優(yōu)化了黏彈性管道瞬變流計算模型[7-9],證明了這些模型的準(zhǔn)確性。

在黏彈性管道瞬變流模型中,管道本構(gòu)模型主要采用廣義開爾文-福伊特(Kelvin-Voigt,KV)力學(xué)模型,該模型中包含蠕變?nèi)崃亢脱舆t時間兩個參數(shù),蠕變參數(shù)對瞬變流壓力波動的影響規(guī)律分析和參數(shù)在本構(gòu)模型中的選取一直是學(xué)者們重點關(guān)注的問題所在。Urbanowicz等[10]分析了不同水溫下蠕變參數(shù)對黏彈性管道瞬變流壓力波曲線的影響規(guī)律。Javadi等[11]分析了黏彈性管道瞬變流壓力波動周期對蠕變參數(shù)的敏感性。為探究蠕變參數(shù)的校核方法,Covas等[5]利用力學(xué)拉伸實驗,測得高密度聚乙烯(high denisity polyethylene, HDPE)管道的本構(gòu)曲線,并將其用于數(shù)值模擬,結(jié)果表明模擬結(jié)果與實驗結(jié)果偏差較大。Pezzinga[12]和Pan等[13-14]分別應(yīng)用遺傳算法和頻域分析方法對黏彈性管道瞬變流模型中的蠕變參數(shù)進(jìn)行校核。但這些參數(shù)校核方法存在校核程序復(fù)雜,耗時長等問題。

對于管道瞬變流動的研究,前人主要針對單一管材的管道系統(tǒng),而針對彈性/黏彈性混合管道系統(tǒng)的研究較少。為闡明混合管道的瞬變流壓力波動影響規(guī)律,Garg等[15]進(jìn)行了鋼管/玻璃纖維塑料混合管道的瞬變流實驗,在實驗中改變塑料管段長度在整個管道系統(tǒng)中的占比,從而分析混合管道瞬變流壓力波動影響規(guī)律。Trabelsi等[16]在彈性管道系統(tǒng)上游設(shè)置塑料短管,研究塑料短管對水錘抑制作用。但目前還未提出適用于彈性/黏彈性混合管道的高效且準(zhǔn)確的參數(shù)校核方法。

因此,現(xiàn)搭建彈性/黏彈性混合管道瞬變流實驗臺,并進(jìn)行了鋼管、鋼管/PPR混合管道、鋼管/HDPE混合管道的瞬變流實驗,分析混合管道瞬變流壓力波動的影響機(jī)制,闡明了波速和蠕變參數(shù)對壓力波動的影響規(guī)律,進(jìn)而提出一種簡單且準(zhǔn)確的混合管道波速和蠕變參數(shù)的校核方法,基于實驗數(shù)據(jù)驗證該方法的準(zhǔn)確性,為工程中混合管道的模擬計算和參數(shù)校核提供參考和依據(jù)。

1 數(shù)學(xué)模型及求解方法

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

經(jīng)典的彈性管道瞬變流控制方程包括連續(xù)性方程和動量方程[17],分別表示為

(1)

(2)

式中:H為壓力水頭;V為平均流速;a為波速;g為重力加速度;x為距離;t為時間;JQ為擬穩(wěn)態(tài)摩阻引起的水頭損失,可以使用達(dá)西-魏斯巴赫(Darcy-Weisbach)公式為

(3)

式(3)中:Q為流量;f為沿程阻力系數(shù);D為管徑;A為管道斷面截面積。

管道瞬變流壓力波的波速a計算公式為

(4)

式(4)中:E為楊氏模量;KL為水的體積模量;ρ為水的密度;e為管道壁厚;C0為約束系數(shù),其計算公式為

(5)

式(5)中:μ為管道泊松比。

相比彈性管道,黏彈性管道瞬變流控制方程中的連續(xù)性方程中增加了延遲應(yīng)變項,用以模擬管道黏彈性效應(yīng),其連續(xù)性方程為

(6)

式(6)中:εr為管道的延遲應(yīng)變。

黏彈性管道的總應(yīng)變ε為

ε=εe+εr

(7)

式(7)中:εe為黏彈性材料的瞬時應(yīng)變。

根據(jù)玻爾茲曼疊加原理,總應(yīng)變ε可表示為

(8)

式(8)中:J0為瞬時蠕變?nèi)崃浚籎(t′)為t′時刻的蠕變?nèi)崃浚沪?t)為t時刻的應(yīng)力。

黏彈性材料的蠕變函數(shù)使用廣義KV模型表式為

(9)

式(9)中:J0為瞬時蠕變?nèi)崃浚籎K為第K個KV元件的蠕變?nèi)崃浚沪觡為延遲時間;n為蠕變元件個數(shù)。

1.2 模型求解

(1)彈性管道瞬變流特征線法求解。使用特征線法,將連續(xù)性方程[式(1)]和動量方程[式(2)]轉(zhuǎn)化為常微分方程,稱之為相容性方程,即

(10)

將式(10)分別沿正負(fù)特征線(圖1中C+和C-特征線)對相容性方程進(jìn)行積分,得到代數(shù)方程為

圖1 差分網(wǎng)格Fig.1 Grid of characteristics

(11)

(2)黏彈性管道瞬變流特征線法求解。使用特征線法,將黏彈性管道連續(xù)性方程[式(6)]和動量方程[式(2)]轉(zhuǎn)化為常微分方程,即

(12)

分別沿正負(fù)特征線對方程進(jìn)行積分,得到代數(shù)方程為

(13)

2 實驗裝置與結(jié)果分析

2.1 實驗設(shè)備及方法

如圖2所示,彈性/黏彈性混合管道瞬變流實驗臺是由供水水箱、管道、氣動閥、替換活接、下游水槽等組成。實驗管段為總長為15 m的鍍鋅鋼管,管道的公稱直徑為DN25、壁厚為3 mm,彈性模量為1.4×1014Pa,泊松比μ為0.23。上游水箱的水位為2.7 m。在管道末端設(shè)置氣動球閥作為管道瞬變流的激勵方式。管道中段兩處活接的間距為5.7 m,活接部分管道占整個管道的38%,關(guān)閥時間0.15 s,線性關(guān)閥,通過實驗數(shù)據(jù)分析和理論計算確定管道沿程阻力系數(shù)f為0.036。實驗設(shè)備的具體型號參數(shù),如表1所示。實驗裝置的局部實物圖如圖3所示。

圖2 實驗裝置示意圖Fig.2 Schematic diagram of experimental setup

表1 設(shè)備及主要參數(shù)Table 1 Main parameters of equipment

本次實驗分兩部分進(jìn)行,第一部分為彈性管道瞬變流實驗,通過關(guān)閉管道下游的氣動閥門產(chǎn)生管道瞬變流動;第二部分為彈性與黏彈性混合管道瞬變流實驗,使用管道中預(yù)先設(shè)置好的活接,將管道中段鋼管替換為高密度聚乙烯(HDPE)和聚丙烯(PPR)管道,從而進(jìn)行彈性/黏彈性混合管道瞬變流實驗,所替換的黏彈性管段的管徑和壁厚與彈性管道相同。所進(jìn)行的實驗工況,如表2所示。

本次實驗所采用的壓力傳感器的精度為±0.25% FS,在本實驗中測得壓力波的峰值為12 m,滿度相對誤差為1.4%,在可接受的范圍內(nèi)。另外,工況2的水溫為24 ℃,其他工況的水溫均為15 ℃。

圖3 實驗裝置局部實物圖Fig.3 Partial physical picture of experimental device

表2 實驗工況基本參數(shù)Table 2 Parameters of experimental tests

2.2 實驗數(shù)據(jù)分析

對上游水箱作為邊界條件、末端關(guān)閥作為瞬變流激勵方式的實驗進(jìn)行分析。在相同實驗參數(shù)(表2)下,分別進(jìn)行鍍鋅鋼管道(工況1)、鋼管/PPR管道(工況4)和鋼管/HDPE(工況5)管道瞬變流壓力波動實驗數(shù)據(jù),如圖4所示。

從圖4可見,彈性管道瞬變流壓力波的最大波峰值為11.3 m,而鋼管/PPR混合管道瞬變流壓力波的最大波峰值為10.2 m,鋼管/HDPE混合管道瞬變流壓力波最大波峰值為8.8 m。與彈性管道瞬變流壓力波曲線相比,鋼管/PPR混合管道和鋼管/HDPE混合管道的壓力波最大峰值分別下降9%和22%。另外,在部分管段替換為HDPE管材后,其產(chǎn)生的瞬變流壓力波曲線相比彈性管道會出現(xiàn)更加嚴(yán)重的衰減現(xiàn)象,但在替換為PPR管道后衰減幅度與彈性管道相似。但是,鋼管/PPR混合管道相比鋼管/HDPE混合管道產(chǎn)生了更大的相位延遲現(xiàn)象。

圖4 工況1、工況4和工況5的實驗數(shù)據(jù)比較Fig.4 Comparison of experimental data for Cases 1, 4 and 5

3 參數(shù)校核方法及驗證

提出了一種合理且簡便的彈性/黏彈性混合管道瞬變流參數(shù)校核方法,并基于不同工況下的實驗數(shù)據(jù),利用黏彈性管道瞬變流模型驗證該校核方法的準(zhǔn)確性。對于黏彈性管道,非穩(wěn)態(tài)摩阻對壓力波動的影響較小,而黏彈性效應(yīng)對管道瞬變流壓力波動的影響占主導(dǎo)地位[18-19],因此本文研究中將使用黏彈性管道瞬變流模型結(jié)合一維擬穩(wěn)態(tài)摩阻模型進(jìn)行校核。

3.1 彈性管道瞬變流波速校核

以工況1的實驗數(shù)據(jù)為例,對彈性管道瞬變流波速進(jìn)行校核。通過實驗數(shù)據(jù)計算波速的范圍區(qū)間,基于實驗數(shù)據(jù)計算波速的方法有以下兩種。

(1)在下游氣動閥瞬時關(guān)閉后,壓力波依次通過壓力傳感器2和壓力傳感器1,分別記錄兩傳感器首次出現(xiàn)壓升的時間、壓力首次達(dá)到波峰15%的時間、壓力首次達(dá)到波峰的時間。根據(jù)公式a=(L1-L2)/(T1-T2),可計算出壓力波波速。

(2)基于壓力波每次到達(dá)穩(wěn)態(tài)的時間計算出壓力波的周期值,根據(jù)公式a=4L/T,計算出壓力波波速。

應(yīng)用所述兩種方法計算出工況1的波速均為350 m/s左右,考慮到實驗數(shù)據(jù)存在誤差,仍需要應(yīng)用管道瞬變流模型對波速進(jìn)行進(jìn)一步校核。

根據(jù)式(4)計算彈性管道的瞬變流壓力波波速,水體積模量KL取2.1×109Pa、空氣體積模量Kg取 1×105Pa、鋼管彈性模量E取1.4×1011Pa。泊松比取0.23。根據(jù)式(5)計算出約束系數(shù)C0為1.17,計算瞬變流波速為1 100 m/s。理論計算結(jié)果同實驗結(jié)果差距過大,推測為管道含氣的原因。根據(jù)文獻(xiàn)[20]中的含氣管道波速計算公式,當(dāng)波速為350 m/s時,管道含氣量為0.000 74。

通過實驗數(shù)據(jù)計算出波速為350 m/s后,以實驗波速的±5%作為取值范圍,通過彈性管道瞬變流模型進(jìn)一步校核波速。在彈性管道瞬變流動過程中,波速是影響壓力波曲線變化的決定性因素,直接影響壓力波的周期和首個波峰值。在取值范圍內(nèi)改變波速,找到壓力波首個波峰和周期與實驗曲線契合程度最好的波速值即為校核結(jié)果。波速值為341 m/s為工況1的校核結(jié)果,模擬結(jié)果如圖5所示。結(jié)果表明,使用一維擬穩(wěn)態(tài)摩阻模型可以較為準(zhǔn)確的模擬出壓力波動的峰值和周期情況,但無法準(zhǔn)確表達(dá)壓力波動的衰減情況,這是由于彈性管道的瞬變流壓力波動的衰減很大程度上受到非穩(wěn)態(tài)摩阻的影響,而一維擬穩(wěn)態(tài)摩阻模型低估了瞬變流壓力波動過程的摩阻耗散[21]。

圖5 工況1的實驗數(shù)據(jù)及模擬結(jié)果對比Fig.5 Comparison of experimental and simulation results in case 1

3.2 混合管道水錘波速校核

用3.1節(jié)相同的方法對工況2的實驗數(shù)據(jù)進(jìn)行計算,計算得出工況2的實驗波速的估算值為230 m/s。與彈性管道相同,實驗計算波速并非準(zhǔn)確的結(jié)果,需要應(yīng)用黏彈性管道瞬變流模型進(jìn)行進(jìn)一步校核。

在參數(shù)校核的過程中,瞬變流波速大小的影響主要體現(xiàn)在壓力波的峰值和周期。但壓力波的峰值和周期同時受到波速和蠕變參數(shù)的影響,僅僅使用波峰作為波速的校核基準(zhǔn)無法得到準(zhǔn)確的結(jié)果。在瞬變流壓力波動過程中前3個周期,管壁的黏彈性對壓力波周期的影響較小,可以忽略不計,但在第4個壓力波動周期以后,受管壁的黏彈性影響,壓力波會產(chǎn)生少量的相位延遲現(xiàn)象。所以,使用壓力波動的前3個周期作為基準(zhǔn)可以較準(zhǔn)確校核出壓力波速值。以本實驗為例,在黏彈性管道瞬變流模型中選取波速值為251 m/s時,模擬結(jié)果的前3個周期與實驗曲線的契合程度最高。

3.3 本構(gòu)參數(shù)校核

對于黏彈性管道的數(shù)值模擬,僅僅通過波速不能充分模擬出壓力波動情況,還需要在本構(gòu)模型(即KV模型)中確定蠕變參數(shù)。并且,由于混合管道中管材占比不同,通過力學(xué)實驗直接測量管道的蠕變參數(shù)的方法是不可行的,只能借助實驗數(shù)據(jù)使用參數(shù)校核的方式得到混合管道的蠕變參數(shù)。

黏彈性管道本構(gòu)模型(即KV模型)常用于表達(dá)黏彈性管道的蠕變動力學(xué)行為。廣義的KV模型是多個KV元件組成的,每個KV元件是由一個彈簧和粘壺并聯(lián)組成的。KV元件的個數(shù)影響蠕變曲線的準(zhǔn)確性,選擇2個及以上的KV元件即可準(zhǔn)確模擬出管道的黏彈性行為[22]。蠕變參數(shù)的影響體現(xiàn)在壓力波的衰減和延遲上。在模型中,蠕變參數(shù)由蠕變?nèi)崃亢脱舆t時間表達(dá)。在參數(shù)校核時將延遲時間固定在0.05 s和0.25 s,通過改變?nèi)渥內(nèi)崃渴沟媚M曲線接近實驗曲線,確定出使模擬結(jié)果最接近實驗曲線的蠕變參數(shù)。

表3列出3組不同大小的蠕變參數(shù),蠕變參數(shù)的大小對壓力波曲線的影響如圖6所示,其中a組的蠕變參數(shù)為工況2校核結(jié)果,模擬曲線同實驗曲線契合程度最高。使用該組蠕變參數(shù)以及3.2節(jié)校核得出得瞬變流波速進(jìn)行模擬計算,模擬結(jié)果如圖7所示。使用黏彈性管道瞬變流模型結(jié)合一維擬穩(wěn)態(tài)摩阻模型可以較好地描述黏彈性管道瞬變流壓力波動過程,只有第2、3個周期壓力波的衰減會有所差異,這是由于在模型中沒有考慮非穩(wěn)態(tài)摩阻的影響。

表3 三組蠕變參數(shù)Table 3 Three sets of creep parameters for simulation

圖6 不同蠕變參數(shù)下PPR混合管道的模擬結(jié)果Fig.6 Simulation results of PPR mixed pipeline under different creep functions

圖7 工況2的實驗數(shù)據(jù)及模擬結(jié)果對比Fig.7 Comparison of experimental and simulation results in case 2

如圖6所示,在延遲時間固定的情況下,b組中使用更大的蠕變參數(shù)會導(dǎo)致更大的壓力波衰減,使得壓力波的過早得趨于平緩,同時相位也會產(chǎn)生更大的延遲現(xiàn)象,c組同理。在校核過程中,第一個波峰值也會在一定程度上受到蠕變參數(shù)的影響,尤其在關(guān)閥時間較長的系統(tǒng)中,蠕變參數(shù)對第一個峰值的影響更為明顯。所以,在校核蠕變參數(shù)的過程中應(yīng)綜合考慮波峰值、衰減幅度和相位延遲情況。

3.4 不同流量的鋼管/PPR混合管道校核模擬

以工況3和工況4的實驗數(shù)據(jù)為基礎(chǔ),對不同流量的工況進(jìn)行模擬校核并驗證不同流量下數(shù)學(xué)模型的準(zhǔn)確性。該兩組實驗在同溫度、同約束條件下進(jìn)行,工況3、工況4的初始流量分別為0.48、0.58 m3/h。

利用前文所提到的主要參數(shù)校核方法,校核得出該兩組實驗的波速為303 m/s,延遲時間τ1和τ2分別為0.05、0.25 s,蠕變?nèi)崃縅1和J2分別為0.19×10-9、0.6×10-9Pa-1。工況3和工況4的波速與本構(gòu)參數(shù)同工況2下稍有不同的原因是因為受到了水溫等外界條件影響。更高的水溫會導(dǎo)致更大的蠕變參數(shù)和更小的波速。

使用校核得出的主要參數(shù)進(jìn)行模擬的結(jié)果,如圖8所示。從圖8可知,使用黏彈性管道瞬變流模型結(jié)合一維擬穩(wěn)態(tài)摩阻模型可以在一定程度上較準(zhǔn)確地模擬出不同初始流量下的混合管道瞬變流壓力波動,模擬結(jié)果與實驗結(jié)果產(chǎn)生一定相位偏差是在可接受的范圍之內(nèi)。

圖8 工況3和工況4的實驗數(shù)據(jù)及模擬結(jié)果對比Fig.8 Comparison of experimental and simulation results in case 3 and case 4

3.5 鋼管/HDPE混合管道的校核模擬

基于工況5的鋼管/HDPE混合管道瞬變流實驗數(shù)據(jù),利用所提到的參數(shù)校核方法進(jìn)行數(shù)值模擬。校核得出該工況下的波速為290 m/s、延遲時間τ1和τ2分別為0.05 s和0.5 s,蠕變?nèi)崃縅1和J2分別為0.13×10-9Pa-1和0.5×10-9Pa-1。模擬結(jié)果與實驗結(jié)果的對比,如圖9所示。使用黏彈性管道瞬變流模型結(jié)合一維擬穩(wěn)態(tài)摩阻模型可以較準(zhǔn)確地模擬出混合管道瞬變流壓力波動前期的峰值和衰減情況,誤差可以保持在4%以下。但模擬結(jié)果與實驗結(jié)果的相位出現(xiàn)一定程度的偏差,達(dá)到10%,在可以接受的范圍之內(nèi)。

4 結(jié)論

通過設(shè)計并搭建彈性/黏彈性混合管道瞬變流實驗臺,分別進(jìn)行了鍍鋅鋼管、鋼管/PPR混合管道以及鋼管/HDPE混合管道的瞬變流實驗,分析彈性/黏彈性混合管道的瞬變流壓力波動機(jī)制。同時提出了一種簡單準(zhǔn)確的波速和蠕變參數(shù)校核方法,并應(yīng)用黏彈性管道瞬變流模型結(jié)合一維擬穩(wěn)態(tài)摩阻模型對該方法的準(zhǔn)確性進(jìn)行驗證。主要得到如下結(jié)論。

圖9 工況5的實驗數(shù)據(jù)及模擬結(jié)果對比Fig.9 Comparison of experimental and simulation results in case 5

(1)在彈性管道系統(tǒng)中部分替換為黏彈性管道可以使得瞬變流壓力波波速變小、峰值下降及相位延遲。這一現(xiàn)象在鋼管/PPR混合管道中相比鋼管/HDPE混合管道更為明顯,在本文的實驗中,對于黏彈性管道占比為38%的混合管道,鋼管/PPR混合管道和鋼管/HDPE混合管道的壓力波最大峰值分別下降9%和22%。

(2)使用彈性管道一維擬穩(wěn)態(tài)摩阻模型無法準(zhǔn)確地模擬出彈性管道瞬變流壓力波動衰減情況。但是不同于彈性管道,使用黏彈性管道瞬變流模型結(jié)合一維擬穩(wěn)態(tài)摩阻模型可以較準(zhǔn)確地對混合管道瞬變流壓力波動進(jìn)行的數(shù)值模擬,盡管在衰減和相位上有一定的偏差,但均在可接受的范圍之內(nèi)。

(3)提出了黏彈性管道瞬變流壓力波波速和蠕變參數(shù)的校核方法。該方法首先通過實驗數(shù)據(jù)計算波速的估算區(qū)間。其次,在此區(qū)間內(nèi)以壓力波前3個周期為基準(zhǔn),校核出與實驗曲線的前3個周期契合程度最高的波速值。最后,以壓力波的首個波峰和周期作為基準(zhǔn),校核出與實驗曲線的峰值和相位契合度最高的蠕變參數(shù)值。

猜你喜歡
實驗
我做了一項小實驗
記住“三個字”,寫好小實驗
我做了一項小實驗
我做了一項小實驗
記一次有趣的實驗
有趣的實驗
小主人報(2022年4期)2022-08-09 08:52:06
微型實驗里看“燃燒”
做個怪怪長實驗
NO與NO2相互轉(zhuǎn)化實驗的改進(jìn)
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
主站蜘蛛池模板: 黄片在线永久| 国产成人精品在线| 国产第一页免费浮力影院| 日韩欧美中文| 国产精品美女自慰喷水| 二级特黄绝大片免费视频大片| 欧美成人亚洲综合精品欧美激情| 在线看片中文字幕| 精品色综合| 99视频免费观看| 亚洲国产成熟视频在线多多 | 国产乱子伦一区二区=| 久久精品国产精品青草app| 欧美日韩一区二区在线播放| 国产精品三级av及在线观看| 综合色在线| 亚洲综合色婷婷中文字幕| 国产拍在线| 欧美在线国产| 免费不卡在线观看av| 午夜不卡视频| 国产欧美亚洲精品第3页在线| 性69交片免费看| 久久久久人妻一区精品色奶水| 97超碰精品成人国产| 亚洲成人播放| 亚洲一区二区三区国产精品| 尤物成AV人片在线观看| 精品人妻无码区在线视频| 99久久99这里只有免费的精品| 日本五区在线不卡精品| 欧美特级AAAAAA视频免费观看| 91精品国产自产91精品资源| 在线观看欧美国产| a毛片免费观看| 国产色伊人| 色婷婷亚洲十月十月色天| 亚洲精品成人福利在线电影| 91极品美女高潮叫床在线观看| 欧美三级自拍| 中文字幕色站| 久草视频精品| 国禁国产you女视频网站| 91网在线| 精品国产美女福到在线不卡f| av手机版在线播放| 亚洲第一视频区| 国内精品91| 国产欧美亚洲精品第3页在线| 大陆国产精品视频| 亚洲中文在线看视频一区| 国产一区二区三区在线精品专区| 秋霞午夜国产精品成人片| 欧美激情视频一区| 国产熟睡乱子伦视频网站| 欧美日在线观看| 麻豆国产精品| 亚洲色图综合在线| 青青草国产在线视频| 真实国产乱子伦视频| 国产成本人片免费a∨短片| 伊伊人成亚洲综合人网7777 | 国产无码精品在线| 亚洲一级毛片在线观播放| 福利国产微拍广场一区视频在线| 欧美成人精品一区二区| 午夜a视频| 久久久久免费精品国产| 欧美日韩中文字幕在线| 色哟哟精品无码网站在线播放视频| 国产91导航| 欧美在线网| 91伊人国产| 久久中文字幕av不卡一区二区| 欧美午夜小视频| 亚洲码一区二区三区| 成人综合久久综合| 亚洲制服丝袜第一页| 欧美午夜在线观看| 国产爽歪歪免费视频在线观看 | 亚洲欧美成人综合| 色吊丝av中文字幕|