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

    考慮氫氣實際氣體效應(yīng)和阻塞流效應(yīng)的螺旋槽干氣密封動態(tài)特性分析

    2017-12-22 05:37:12許恒杰宋鵬云毛文元鄧強國
    化工學(xué)報 2017年12期
    關(guān)鍵詞:干氣氣膜端面

    許恒杰,宋鵬云,毛文元,鄧強國

    (1昆明理工大學(xué)化學(xué)工程學(xué)院,云南 昆明 650500;2昆明理工大學(xué)環(huán)境科學(xué)與工程學(xué)院,云南 昆明 650500)

    考慮氫氣實際氣體效應(yīng)和阻塞流效應(yīng)的螺旋槽干氣密封動態(tài)特性分析

    許恒杰1,2,宋鵬云1,毛文元1,鄧強國1

    (1昆明理工大學(xué)化學(xué)工程學(xué)院,云南 昆明 650500;2昆明理工大學(xué)環(huán)境科學(xué)與工程學(xué)院,云南 昆明 650500)

    以 Chen等提出的氫氣實際氣體狀態(tài)方程描述氫氣的實際氣體行為,以氣體出口速度達到音速作為產(chǎn)生阻塞效應(yīng)的條件,確定出口壓力邊界條件,采用小擾動法分析了干氣密封操作參數(shù)對螺旋槽干氣密封動態(tài)特性的影響規(guī)律,并與理想氣體、強制出口壓力邊界模型中的動態(tài)特性系數(shù)進行了對比。結(jié)果表明:研究高壓螺旋槽干氣密封的動態(tài)特性時應(yīng)當(dāng)考慮實際氣體效應(yīng)和阻塞流效應(yīng)。兩種效應(yīng)使氫氣螺旋槽干氣密封的直接動態(tài)氣膜剛度減小,使直接動態(tài)氣膜阻尼增大。隨著壓縮數(shù)、進口壓力的增大,兩種效應(yīng)對動態(tài)氣膜剛度的影響逐漸增強。以頻率比為變量時,兩種效應(yīng)主要影響氣膜剛度,對氣膜阻尼的影響作用較小。針對所研究的工況,與理想氣體和強制壓力出口邊界條件相比,考慮實際氣體效應(yīng)和阻塞流效應(yīng),以壓縮數(shù)為變量時,動態(tài)氣膜阻尼(Czz、Cαα、Cαβ)的平均偏差分別為2.28%、1.93%、2.79%;以進口壓力為變量時,3種氣膜阻尼的平均偏差分別達到4.08%、2.07%、1.82%。

    螺旋槽干氣密封;動態(tài)特性;實際氣體;阻塞流;數(shù)值分析

    引 言

    在石油煉制、石油化工和煤化工行業(yè)中,加氫裝置中的氫氣壓縮機是一種關(guān)鍵設(shè)備,而泵入式螺旋槽干氣密封是氫氣壓縮機里應(yīng)用較多的一種干氣密封[1],靠近介質(zhì)側(cè)干氣密封的沖洗潤滑介質(zhì)一般為氫氣。泵入式螺旋槽干氣密封的高壓側(cè)位于密封端面外徑處,低壓側(cè)位于密封端面內(nèi)徑處,一般與緩沖氣體或大氣相通。在徑向方向上,由于密封間隙進出口之間的壓差作用,潤滑氣體在端面外徑處被泵入密封間隙,在內(nèi)徑處流出密封間隙。此過程中氣體的流通面積是逐漸減小的,且潤滑氣體在密封壩區(qū)的流動是一個降壓過程,從質(zhì)量守恒和能量守恒兩種角度出發(fā),潤滑氣體在密封間隙內(nèi)流動時速度會逐漸增大,一旦此速度達到音速,會使密封端面出口發(fā)生阻塞,導(dǎo)致出口壓力大于環(huán)境壓力,改變密封端面氣膜壓力分布規(guī)律[2],從而會對泵入式螺旋槽干氣密封性能產(chǎn)生一定的影響。

    目前在干氣密封性能的研究中,大多數(shù)學(xué)者均將氣體的出口壓力視為環(huán)境壓力[3-5](稱為強制出口壓力邊界),考慮阻塞流效應(yīng)的干氣密封性能研究較少。1971年,針對光滑平面的泵出式端面密封,Zuk等[6]給出了密封發(fā)生阻塞時的出口壓力解析表達式,在進口壓力一定的前提下分析了進出口壓比對氣體泄漏率的影響,同時通過實驗驗證了其理論分析的正確性。Arghir等[7]對環(huán)形密封中的阻塞流提出一種數(shù)值計算方法,研究了存在阻塞流時環(huán)形密封的靜態(tài)穩(wěn)定性,指出若在出口處發(fā)生阻塞,環(huán)形密封可能發(fā)生失穩(wěn)。Thomas等[8]認(rèn)為氣體在密封間隙出口發(fā)生阻塞時端面泄漏將達到最大,并將這一條件作為阻塞壓力出口條件對端面光滑的錐形干氣密封進行了熱彈流潤滑分析。馬春紅等[9]求解雷諾方程時將Zuk阻塞壓力表達作為出口邊界條件,研究了高速螺旋槽干氣密封的穩(wěn)態(tài)性能。同樣,Thatte等[10]采用Zuk阻塞壓力邊界研究了螺旋槽干氣密封壩區(qū)的Mach數(shù)和壓力分布規(guī)律。以上研究報道均是針對阻塞流效應(yīng)對干氣密封穩(wěn)態(tài)特性的影響。而工程應(yīng)用已表明,干氣密封氣膜的動態(tài)特性不良是密封喪失穩(wěn)定性的原因之一[11],表征氣膜動特性的兩個系數(shù)(動態(tài)氣膜剛度、阻尼)反映了氣膜的動態(tài)穩(wěn)定性能[12]。因此,開展干氣密封氣膜動態(tài)特性的研究對設(shè)計、評價干氣密封的動態(tài)穩(wěn)定性能具有重要意義。干氣密封的動態(tài)性能分析是在其穩(wěn)態(tài)性能的基礎(chǔ)上展開的,發(fā)生阻塞流效應(yīng)后,其動態(tài)特性必然會與亞音速流動的情況存在差異。迄今為止,尚未發(fā)現(xiàn)考慮出口阻塞流影響的干氣密封動態(tài)特性研究,而開展干氣密封動態(tài)特性研究是干氣密封穩(wěn)定性研究至關(guān)重要的部分。此外,大多數(shù)研究沒有考慮潤滑氣體的實際氣體行為,根據(jù)已有報道[13-14],實際氣體效應(yīng)對干氣密封的動態(tài)特性存在著一定的影響。綜合考慮,忽略實際氣體效應(yīng)和阻塞流效應(yīng)會影響干氣密封性能的準(zhǔn)確分析與評價。

    本文以氫氣為例,采用不同的實際氣體狀態(tài)方程計算氫氣壓縮因子,對密封出口速度進行判定,建立了考慮實際氣體效應(yīng)和阻塞流效應(yīng)的泵入式螺旋槽干氣密封動態(tài)特性分析數(shù)學(xué)模型,并通過算例驗證了所建模型的正確性,研究了兩種效應(yīng)對螺旋槽干氣密封動態(tài)特性的影響,并與理想氣體、強制出口壓力邊界算例進行了對比,為螺旋槽干氣密封的工程設(shè)計提供一定的理論依據(jù)。

    1 幾何模型

    泵入式螺旋槽干氣密封動態(tài)特性的分析模型如圖1所示,動環(huán)端面外側(cè)周向均勻開有12個螺旋槽,潤滑氣體被泵入密封間隙后,由于槽結(jié)構(gòu)的特殊性和一定的轉(zhuǎn)速,在槽根處會因壓縮產(chǎn)生氣體動壓使干氣密封動靜環(huán)之間形成一定厚度的氣膜,最后經(jīng)由密封間隙內(nèi)徑流出。

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

    2.1 壓縮因子

    壓縮因子Z實質(zhì)上表示實際氣體比容Vre與相同溫度、相同壓力下理想氣體比容Vid之比,按理想氣體狀態(tài)方程的形式pVid=RT,實際氣體狀態(tài)方程可以寫成

    明顯地,當(dāng)Z=1時,式(1) 為理想氣體狀態(tài)方程;Z<1時,可認(rèn)為是同溫同壓下,實際氣體的比容小于理想氣體的比容,即實際氣體更易于壓縮;當(dāng)Z>1時情況則恰好相反。

    圖1 螺旋槽干氣密封動態(tài)特性分析模型[15]Fig.1 Dynamic analysis model of spiral groove dry gas seal[15]

    建立一個既適用于所有氣體,而且壓力、溫度又包括寬廣的范圍、同時能確切反映實際氣體行為的通用方程是很困難的。因此,實際氣體的壓縮因子需要根據(jù)氣體的種類選用合理的狀態(tài)方程進行計算。目前,計算壓縮因子的手段有多種,主要包括Redlich-Kwong(R-K)方程、維里(Virial)方程和擬合NIST數(shù)據(jù)等,維里方程為二階項維里方程,Chen等[16]通過對NIST數(shù)據(jù)進行擬合,在溫度為173~393 K、壓力為0.1~100 MPa范圍內(nèi),得出了一個氫氣的實際氣體狀態(tài)方程。

    以氫氣為例,采用以上3種方式計算的壓縮因子Z隨壓力的變化趨勢如圖2所示??梢钥闯觯S著壓力的增大,R-K方程、Chen等提出的氫氣狀態(tài)方程與NIST數(shù)據(jù)均保持良好的吻合性,相對于其他兩種方程,二階項維里方程與NIST數(shù)據(jù)之間的偏差略大。為了確定精確的氫氣壓縮因子表達方式,計算3種壓縮因子方程與NIST數(shù)據(jù)之間的平均誤差分別為0.272%、0.263%、1.985%。因此,本文選用Chen等提出的氫氣狀態(tài)方程表達氫氣的實際行為。

    圖2 壓縮因子隨壓力變化的趨勢(T=303.15 K)Fig.2 Compressibility factor versus pressure (T=303.15 K)

    Chen的氫氣狀態(tài)方程中的壓縮因子為

    式中,B為常數(shù),其值為1.9155×10-6K·Pa-1。

    2.2 實際氣體的音速

    式中,kv=Z/(Z-RZT2/cp)為實際氣體的容積絕熱指數(shù),其中ZT為導(dǎo)數(shù)壓縮因子,cp為實際氣體的比定壓熱容??梢钥闯觯?dāng)氣體為理想氣體時,kv等于理想氣體的比熱比γ,式(3)可以化為理想氣體的音速公式c=(γRgT)0.5。導(dǎo)數(shù)壓縮因子ZT的計算如下[18]

    實際氣體比定壓熱容cp的計算按文獻[19]的方法進行。

    2.3 壓力控制方程

    螺旋槽干氣密封的端面氣膜厚度一般為微米級,與端面氣膜的徑向和周向尺度相比,其軸向方向上的壓力梯度可以忽略不計,則密封端面間的氣膜壓力選用二維柱坐標(biāo)下可壓縮流體的雷諾方程進行描述

    式中,η為潤滑氣體的黏度;ρ為潤滑氣體的密度,ω為干氣密封的轉(zhuǎn)速。

    文獻[9]對絕熱流動和等溫流動兩種假設(shè)下螺旋槽干氣密封發(fā)生阻塞時的質(zhì)量流量進行了計算,發(fā)現(xiàn)等溫流動假設(shè)下其理論計算結(jié)果更接近Zuk的實驗結(jié)果。本文假設(shè)密封端面間的氣體流動為等溫流動,潤滑氣體為牛頓流體。

    聯(lián)立式(1)和式(3)并代入式(5)中即可得考慮實際氣體效應(yīng)的螺旋槽干氣密封壓力控制方程。在干氣密封的實際運行過程中,旋轉(zhuǎn)動環(huán)往往會發(fā)生軸向微振和角向偏擺,氣膜厚度在軸向和角向方向上會產(chǎn)生微小的波動。因此,將干氣密封受到外部激勵ω1而發(fā)生的微小運動視為在平衡運動狀態(tài)上的疊加,采用小擾動法可推導(dǎo)出考慮實際氣體效應(yīng)的量綱1穩(wěn)態(tài)、動態(tài)雷諾方程,具體的推導(dǎo)過程可參考文獻[13,15,20]。

    式中,量綱1變量定義如下

    其中,p0為穩(wěn)態(tài)壓力,ppi為參考壓力,文中取為5 MPa,hg為槽深,h0為平衡膜厚,ro為密封環(huán)外徑,Λ為壓縮數(shù),Γ為頻率比,σ為頻率數(shù),pzr、pzi為軸向微擾壓力。在圖1中,定義α向為動環(huán)繞x軸偏擺的方向,則pαr、pαi為使動環(huán)繞x軸偏擺的角向微擾壓力,β向為動環(huán)繞y軸偏擺的方向,則pβr、pβi為使動環(huán)繞y軸偏擺的角向微擾壓力。

    2.4 定義端面氣膜動態(tài)特性系數(shù)

    如圖1所示,密封環(huán)受3個自由度的微擾,動態(tài)特性參數(shù)有18個。已有的三自由度微擾下干氣密封動態(tài)特性研究表明[21],干氣密封受三自由度微擾的運動可以簡化為兩種相互獨立的微擾運動,即軸向的微擾運動和沿兩個正交軸作角向擺動的微擾運動,因此軸向微擾引起的角向動特性系數(shù)或角向微擾引起的軸向動特性系數(shù)可以忽略不計。同時,由于端面氣膜的軸對稱性,α向和β向微擾引起的動態(tài)特性系數(shù)滿足以下關(guān)系:Kαα=Kββ、Kαβ=-Kβα、Cαα=Cββ、Cαβ=-Cβα。這里將由j向微擾引起的j向動態(tài)剛度、阻尼稱為直接動態(tài)剛度、阻尼,其余的稱為交叉耦合動態(tài)剛度、阻尼,其中j=α、β。即干擾自身方向的動態(tài)剛度和阻尼稱為直接動態(tài)剛度、直接動態(tài)阻尼;干擾引起的另外兩個方向的剛度和阻尼,稱為交叉耦合剛度和交叉耦合阻尼。因此,本文具體以軸向直接動態(tài)剛度、阻尼系數(shù)(Kzz、Czz),角向直接動態(tài)剛度、阻尼系數(shù)(Kαα、Cαα),角向交叉耦合動態(tài)剛度、阻尼系數(shù)(Kαβ、Cαβ)為研究對象,其余12個動特性系數(shù)不予以討論。定義端面

    氣膜動態(tài)剛度和動態(tài)阻尼系數(shù)計算式如下[22]

    因此,有量綱動態(tài)特性系數(shù)與量綱1動態(tài)特性系數(shù)之間的關(guān)系可根據(jù)量綱1變量的定義導(dǎo)出

    動態(tài)剛度反映了受到微擾時氣膜抵抗變形的能力,其包含穩(wěn)態(tài)氣膜剛度;動態(tài)阻尼反映了受到微擾時氣膜抑制振動發(fā)散的能力,即阻尼越大,氣膜越易趨于穩(wěn)定。

    3 數(shù)值求解

    3.1 計算方法

    采用向后和中心差分,可以將式(6)~式(12)離散,具體的差分離散格式可以參考文獻[23]??紤]到密封間隙出口處可能會發(fā)生阻塞,在迭代計算過程中需對氣體的出口速度進行判定,從而選用合適的出口壓力邊界。由于在端面氣膜密封中,對于密封壩區(qū)的氣體泄漏通常采用經(jīng)典的黏性可壓流理論(即與膜厚3次方的關(guān)系)來分析,且文獻[15]已指出徑向泄漏率主要由壓差流引起,轉(zhuǎn)速引起的剪切流對徑向泄漏率的影響不大,即旋轉(zhuǎn)效應(yīng)對徑向速度和徑向泄漏率的影響可以忽略,因此上述的出口速度指的是潤滑氣體的徑向速度。值得注意的是,這里僅是指氣膜厚度確定后,計算泄漏率和出口速度時不考慮旋轉(zhuǎn)速度的影響。實際上,旋轉(zhuǎn)速度對端面氣膜厚度,即密封間隙的形成有重要貢獻,在計算氣膜厚度時考慮了轉(zhuǎn)速的影響。此外,在目前干氣密封端面間隙出口處考慮阻塞流效應(yīng)的研究文獻中,都有考慮入口壓力損失,但是通過分析文獻[9]中的數(shù)據(jù)發(fā)現(xiàn),針對泵入式螺旋槽干氣密封,考慮入口損失后的進口壓力與原進口壓力僅僅相差 1%左右,因此為了簡化計算,本文忽略入口壓力損失。

    在穩(wěn)態(tài)壓力的迭代計算過程中,首先對出口壓力邊界賦一初值(出口處環(huán)境壓力),迭代計算后獲得出口速度uexit,然后與此出口壓力下的實際氣體音速us進行比較,若uexit小于或等于us,則以此出口壓力作為出口的壓力邊界條件進行后續(xù)計算;若uexit大于us,則需對出口壓力邊界重新賦值,直至uexit等于us為止[24]后才進行穩(wěn)態(tài)、動態(tài)壓力迭代計算。以軸向動態(tài)剛度與阻尼為例,具體的計算流程如圖3所示,角向動態(tài)剛度與阻尼的計算過程與之類似。

    圖3 計算流程Fig.3 Calculation procedure

    3.2 邊界條件

    在迭代計算穩(wěn)態(tài)、動態(tài)氣膜壓力時,需要給定邊界條件,在迭代中邊界條件分為兩類。

    (1)在密封環(huán)的外、內(nèi)側(cè),即潤滑氣體進、出口處,軸向、角向微擾壓力和穩(wěn)態(tài)進口壓力分別滿足[25]Pjr=0,Pji=0,Po=po/ppi

    (2)周期性邊界條件

    這里P為量綱1穩(wěn)態(tài)壓力、微擾壓力的統(tǒng)稱,包含P0、Pjr、Pji,其中j=z、α、β。

    4 算例驗證

    4.1 出口阻塞流算法驗證

    在密封環(huán)內(nèi)、外徑分別為70 mm和90 mm,進口壓力為 20 MPa,密封間隙為 3.17 μm,轉(zhuǎn)速為261.8 rad·s-1的工況下,計算了光滑平行端面干氣密封的壓力分布,并與文獻[8]的數(shù)據(jù)進行了對比,對比結(jié)果如圖4所示。可以看出本文算法所得的結(jié)果與Thomas等[8]的計算結(jié)果趨勢一致,出口壓力約相差4%。說明本文的出口阻塞流算法是合理的。

    圖4 阻塞出口邊界壓力計算值與文獻值對比Fig.4 Comparison of steady pressure distribution between current study and Thomas’ at choked boundary

    4.2 密封動態(tài)特性算法驗證

    在潤滑氣體為理想氣體、密封間隙進出口均采用強制壓力出口邊界的前提下,Ruan[26]采用有限元法給出了三自由度微擾下螺旋槽干氣密封的動態(tài)剛度和阻尼數(shù)據(jù)。本節(jié)選用 Ruan的計算參數(shù),將壓縮因子取為 1,即B=0,同時不對密封間隙的出口速度進行判定,選定大氣壓力作為干氣密封出口壓力,選取軸向動態(tài)氣膜剛度和阻尼作為驗證參數(shù),對本文的動態(tài)特性算法進行驗證。

    案例計算參數(shù):密封環(huán)內(nèi)徑ri=30 mm;外徑ro=42 mm;螺旋槽槽根半徑rg=33.6 mm;螺旋角α=20°;槽深hg=5 μm;槽臺比為 1;氣體黏度η=1.79×10-5Pa·s;大氣壓力pa=0.101 MPa;進口壓力po=0.303 MPa。

    將計算值轉(zhuǎn)化成有量綱的剛度與阻尼,本文算法所計算的螺旋槽干氣密封軸向動態(tài)氣膜剛度和阻尼系數(shù)與文獻值的對比如圖5所示,可以看出,隨著頻率比Γ的增大,本文的剛度系數(shù)和阻尼系數(shù)計算結(jié)果與文獻值趨勢一致,數(shù)值相差不大,這說明本文的差分格式、網(wǎng)格密度、迭代精度控制都是合理的,軸向動態(tài)剛度系數(shù)和阻尼系數(shù)的計算程序是有效的。

    圖5 動態(tài)特性系數(shù)計算值與文獻值對比Fig.5 Comparison of dynamic force coefficients between current study and Ruan’s

    5 結(jié)果分析與討論

    在此節(jié)中,螺旋槽干氣密封的槽深hg=3 μm,平衡膜厚h0=3 μm,其余的結(jié)構(gòu)參數(shù)沿用4.2節(jié)中的數(shù)據(jù)。由于螺旋槽干氣密封的動特性系數(shù)隨其結(jié)構(gòu)參數(shù)的變化規(guī)律已有了充分的研究[27-29],因此本文側(cè)重于密封操作參數(shù),選取壓縮數(shù)Λ、頻率比Γ和進口壓力po作為變量(其中壓縮數(shù)可反映干氣密封的轉(zhuǎn)速),探討實際氣體效應(yīng)和阻塞流效應(yīng)對螺旋槽干氣密封動特性系數(shù)的影響。潤滑氣體選用氫氣,其溫度為 298.15 K,黏度為η=8.9135×10-6Pa·s,氣體常數(shù)R=8.314 J·(mol·K)-1。需要說明的是,文中的算法是經(jīng)過量綱1化的,因此若沒有特別說明,后續(xù)分析中的參數(shù)均為量綱1值。

    此外,定義了4種計算模型:考慮實際氣體效應(yīng)和阻塞流效應(yīng)模型(CaseⅠ);實際氣體和強制壓力出口邊界模型(Case Ⅱ);理想氣體和阻塞壓力出口邊界模型(Case Ⅲ);理想氣體和強制壓力出口邊界模型(Case Ⅳ)。

    5.1 壓縮數(shù)Λ的影響

    選定頻率數(shù)σ=20,進口壓力為20 MPa??紤]實際氣體效應(yīng)、阻塞流效應(yīng)和理想氣體、強制壓力出口邊界兩種模型中的動態(tài)剛度隨壓縮數(shù)的變化,結(jié)果如圖6所示。從壓縮數(shù)的定義可以得知,當(dāng)其他參數(shù)保持不變,壓縮數(shù)可以表征干氣密封的轉(zhuǎn)速大小。從圖6可以看出,就絕對值而言,在兩種模型中,軸向、角向直接動態(tài)剛度(Kzz、Kαα)、角向交叉耦合動態(tài)剛度Kαβ均隨壓縮數(shù)的增大而近似呈線性增長??梢哉J(rèn)為這是由較高的轉(zhuǎn)速導(dǎo)致密封端面間出現(xiàn)了較強的動壓效應(yīng),從而使得密封端面間壓力增大所致。對比兩種模型中的動態(tài)剛度計算結(jié)果,可以發(fā)現(xiàn)考慮實際氣體效應(yīng)和阻塞流效應(yīng)時的動態(tài)剛度小于理想氣體、強制壓力出口邊界時的動態(tài)剛度,動態(tài)剛度在兩種模型之間的偏離程度隨壓縮數(shù)的增大而逐漸變得顯著,在所研究的壓縮數(shù)范圍內(nèi),3種動態(tài)剛度在兩種模型之間的平均誤差分別約為 9.17%、11.6%和 0.958%,最大偏差分別為9.52%、16.8%、1.99%。

    圖6 壓縮數(shù)對動態(tài)剛度的影響Fig.6 Dynamic stiffness coefficients change with squeeze number

    圖7 壓縮數(shù)對動態(tài)阻尼的影響Fig.7 Dynamic damping coefficients change with squeeze number

    壓縮數(shù)對動態(tài)阻尼系數(shù)的影響如圖7所示,兩種模型中的量綱1動態(tài)阻尼系數(shù)隨著壓縮數(shù)的增大而增大,這與文獻[15]中螺旋槽順流泵端面氣膜密封阻尼系數(shù)隨壓縮數(shù)的變化規(guī)律相似,但是有量綱直接動態(tài)阻尼系數(shù)的變化趨勢恰好與之相反,具體數(shù)據(jù)如表1所示。壓縮數(shù)增大,有量綱直接動態(tài)阻尼減小意味著高轉(zhuǎn)速下氣膜的穩(wěn)定性變差。同時從表1可以看出,有量綱直接動態(tài)阻尼系數(shù)czz的值要遠遠大于cαα的值,這說明高轉(zhuǎn)速工況下密封容易失穩(wěn),且相較于軸向運動,隨著轉(zhuǎn)速的增大端面氣膜密封的角向運動更容易發(fā)生失穩(wěn)。這些狀況與端面氣膜密封穩(wěn)定性分析和實踐運行的結(jié)果吻合。另外,CaseⅠ中的動態(tài)阻尼大于 Case Ⅳ中的動態(tài)阻尼。兩種模型之間的平均動態(tài)阻尼偏差分別為2.28%、1.93%和2.79%,最大偏差分別為5.77%、5%、3.91%。

    表1 有量綱直接動態(tài)氣膜阻尼隨壓縮數(shù)的變化Table 1 Dimension direct dynamic damping coefficients change with squeeze number

    實際氣體效應(yīng)和阻塞流效應(yīng)分別對干氣密封氣膜的動態(tài)特性系數(shù)產(chǎn)生何種具體的影響很難從圖6和圖7看出。因此,以軸向直接動態(tài)剛度和阻尼(均為有量綱值)作為分析對象,比較了4種計算模型中的氣膜動特性系數(shù),如圖8所示。根據(jù)前述對4種計算模型的定義,CaseⅠ和CaseⅢ、CaseⅡ和 CaseⅣ中動態(tài)特性參數(shù)之間的差異表示實際氣體效應(yīng)對干氣密封動特性的影響;CaseⅠ和 CaseⅡ、CaseⅢ和CaseⅣ中的動態(tài)特性參數(shù)之間的差異表示阻塞流效應(yīng)對干氣密封動特性的影響??梢钥闯?,4種模型中計算所得的氣膜軸向直接動態(tài)剛度滿足以下關(guān)系:CaseⅠ<Case Ⅲ、CaseⅡ<Case Ⅳ;CaseⅠ< CaseⅡ、Case Ⅲ<Case Ⅳ,這說明當(dāng)潤滑氣體為氫氣時,不管密封出口是否發(fā)生阻塞,實際氣體效應(yīng)均使動態(tài)氣膜剛度減??;不管潤滑氣體被視為實際氣體還是理想氣體,阻塞流效應(yīng)同樣使動態(tài)氣膜剛度減小。而氣膜軸向直接動態(tài)阻尼則滿足:CaseⅠ> Case Ⅲ、CaseⅡ>Case Ⅳ;CaseⅠ>CaseⅡ、Case Ⅲ>Case Ⅳ,這說明實際氣體效應(yīng)和阻塞流效應(yīng)均使動態(tài)氣膜阻尼增大。對于兩種效應(yīng)使動態(tài)剛度減小、使動態(tài)阻尼增大這一現(xiàn)象,這是由于氫氣的壓縮因子大于 1,較理想氣體而言,氫氣實際氣體具有更大的比體積;再者阻塞流效應(yīng)會提升密封間隙的出口壓力,導(dǎo)致干氣密封進出口壓差減小。在氣膜厚度一定的情況下,兩種效應(yīng)均使得密封的泄漏率減小,此時相較于理想氣體,氫氣實際氣體密度較小,氣膜吸收外界干擾能量的能力增強,抑制振動發(fā)散的能力提高,因此將潤滑氣體視為理想氣體和考慮阻塞流效應(yīng)時,氣膜具有較小的動態(tài)剛度和較大的動態(tài)阻尼。

    圖8 4種計算模型中有量綱直接動態(tài)剛度、阻尼對比Fig.8 Comparison of dimension dynamic stiffness and damping coefficients in four calculated cases

    5.2 頻率比Γ的影響

    選定壓縮數(shù)Λ=50,進口壓力為20 MPa。螺旋槽干氣密封動特性系數(shù)隨頻率比Γ的變化規(guī)律分別如圖 9、圖 10所示。可以看出,在低頻率比范圍(Γ=0.01~0.5)內(nèi)動態(tài)剛度和動態(tài)阻尼幾乎保持恒定。隨著頻率比的增大,動態(tài)剛度(Kzz、Kαα、Kαβ)均增大,動態(tài)阻尼(Czz、Cαα、Cαβ)均減小,其中動態(tài)剛度和動態(tài)阻尼在Γ=0.5~10范圍內(nèi)變化幅度非常明顯,在Γ>10以后動態(tài)剛度和動態(tài)阻尼的變化幅度趨于平緩。文獻[30]指出,對于密封-轉(zhuǎn)子系統(tǒng),直接阻尼越大越有利于密封的穩(wěn)定。因此從阻尼的角度出發(fā),外界擾動頻率越低,干氣密封越容易趨于穩(wěn)定,這同樣符合端面氣膜密封穩(wěn)定性分析和實際運行的結(jié)果。

    圖9 頻率比對動態(tài)剛度的影響Fig.9 Dynamic stiffness coefficients vary with frequency ratio

    圖10 頻率比對動態(tài)阻尼的影響Fig.10 Dynamic damping coefficients vary with frequency ratio

    從以上兩圖中還可以看出,不同頻率比下,考慮實際氣體、阻塞流效應(yīng)模型中的氣膜動態(tài)剛度系數(shù)比理想氣體、強制出口壓力邊界模型中的動態(tài)剛度小,動態(tài)阻尼則表現(xiàn)出相反的變化趨勢。在所研究的頻率比范圍內(nèi),對于動態(tài)氣膜剛度(Kzz、Kαα、Kαβ),兩種模型之間的平均誤差分別為 7.04%、8.64%和7.23%,最大誤差分別為9.03%、11.7%和18.8%;對于動態(tài)氣膜阻尼(Czz、Cαα、Cαβ),兩種模型之間的平均誤差分別為 1.87%、1.29%和2.41%,最大誤差分別為4.45%、3.55%和6.26%。對比以上兩組誤差,顯然實際氣體效應(yīng)和阻塞流效應(yīng)對動態(tài)剛度的影響更為顯著。

    5.3 進口壓力po的影響

    選定頻率數(shù)σ=100,壓縮數(shù)Λ=50。不同進口壓力下的動態(tài)氣膜剛度變化如圖11所示??梢钥闯鲭S著進口壓力的增大,就絕對值而言,3種動態(tài)氣膜剛度均持續(xù)增大,且在高壓范圍內(nèi),3種動態(tài)氣膜剛度的增大幅度均變緩,其中角向直接動態(tài)剛度Kαα曲線近似水平。三者的增長幅度并不相同,其中以Kzz最大,Kαβ次之,Kαα最小。相較于理想氣體、強制壓力出口邊界條件,實際氣體和阻塞流效應(yīng)使氣膜動態(tài)剛度減小。此外,隨著進口壓力的增大,實際氣體效應(yīng)和阻塞流效應(yīng)對氣膜動態(tài)剛度的影響越發(fā)明顯,即考慮實際氣體效應(yīng)、阻塞流效應(yīng)和理想氣體、強制壓力出口邊界兩種模型間的動態(tài)剛度之差隨進口壓力的增大而增大,說明進口壓力增大會導(dǎo)致實際氣體效應(yīng)、阻塞流效應(yīng)對氣膜動態(tài)剛度的影響逐漸變強。當(dāng)進口壓力為20 MPa時,兩種模型之間動態(tài)氣膜剛度(Kzz、Kαα、Kαβ)的誤差分別為8.46%、10.7%和0.25%。

    圖11 進口壓力對動態(tài)剛度的影響Fig.11 Dynamic stiffness coefficients vary with sealed pressure

    圖12 進口壓力對動態(tài)阻尼的影響Fig.12 Dynamic damping coefficients vary with sealed pressure

    不同進口壓力下的動態(tài)氣膜阻尼變化如圖 12所示。從圖中可以看出直接動態(tài)阻尼(Czz、Cαα)隨進口壓力的增大而增加,其增大幅度逐漸變緩;在1~10 MPa范圍內(nèi)角向交叉耦合阻尼系數(shù)Cαβ表現(xiàn)出與直接動態(tài)阻尼相同的趨勢,隨著進口壓力的進一步增加,Cαβ逐漸減小。相較于理想氣體、強制壓力出口邊界條件,實際氣體效應(yīng)和阻塞流效應(yīng)使直接氣膜動態(tài)阻尼略微變大,但是對于角向交叉阻尼,則呈現(xiàn)出較復(fù)雜的變化規(guī)律,這可能與外界干擾頻率較大有關(guān),具體原因還須進一步的探討。兩種模型之間動態(tài)氣膜阻尼(Czz、Cαα、Cαβ)的平均誤差分別為4.08%、2.07%和1.82%,最大偏差分別為12.6%、3.65%、3.25%。

    對于某一特定的動特性系數(shù),在1~3 MPa范圍內(nèi)兩種模型中的動特性系數(shù)曲線幾乎重合,這是由于在此進口壓力范圍內(nèi),密封間隙的出口并沒有阻塞流出現(xiàn),且實際氣體效應(yīng)在此壓力范圍內(nèi)可以忽略的緣故。

    6 結(jié) 論

    (1)相較于理想氣體和強制出口壓力邊界模型,實際氣體效應(yīng)和阻塞流效應(yīng)使氫氣螺旋槽干氣密封的直接動態(tài)氣膜剛度減小,使直接動態(tài)氣膜阻尼增大。針對算例,頻率比Γ在0.01~20范圍內(nèi),兩種效應(yīng)使軸向動態(tài)氣膜剛度Kzz平均減小了7.04%,最大減少了 9.03%;使軸向動態(tài)氣膜阻尼Czz平均增大了約1.87%,最大增幅達到4.45%。

    (2)實際氣體效應(yīng)和阻塞流效應(yīng)對動態(tài)氣膜剛度的影響隨壓縮數(shù)、進口壓力的增大而增強。針對計算案例,軸向動態(tài)氣膜剛度隨壓縮數(shù)變化的最大增幅達到 9.52%、隨進口壓力增加的最大增幅達8.46%。在低壓范圍(1~3 MPa)內(nèi),兩種模型中的剛度曲線幾乎重合,可以忽略兩種效應(yīng)的影響。

    (3)針對所研究工況,與理想氣體和強制邊界條件相比,實際氣體效應(yīng)和阻塞流效應(yīng)使動態(tài)氣膜阻尼發(fā)生改變。壓縮數(shù)Λ在20~120范圍內(nèi),動態(tài)氣膜阻尼(Czz、Cαα、Cαβ)的平均偏差分別為 2.28%、1.93%、2.79%;最大偏差分別為5.77%、5%、3.91%。進口壓力1~20 MPa范圍內(nèi),3種氣膜阻尼的平均偏差分別達到4.08%、2.07%、1.82%,最大偏差分別為12.6%、3.65%、3.25%。

    [1]程海峰,王為民,王艷東,等.TM02D型串聯(lián)式干氣密封在重整循環(huán)氫氣壓縮機上的應(yīng)用[J].潤滑與密封,2010,35(12):119-122.CHENG H F,WANG W M,WANG Y D,et al.Application of TM02D-tandeming dry gas seal on reforming circulation hydrogen compressor[J].Lubrication Engineering,2010,35(12):119-122.

    [2]THOMAS S,BRUNETIERE N,TOUMERIE B.Numerical modeling of high pressure gas face seals[J].Journal of Tribology,2006,128(2):241-242.

    [3]產(chǎn)文,宋鵬云,毛文元,等.螺旋槽干氣密封端面氣膜溫度場的數(shù)值分析[J].排灌機械工程學(xué)報,2015,33(5):422-428.CHAN W,SONG P Y,MAO W Y,et al.Numerical an alysis of temperature field of gas film in spiral groove dry gas seal[J].Journal of Drainage and Irrigation Machinery Engineering,2015,33(5):422-428.

    [4]宋鵬云,產(chǎn)文,毛文元,等.實際氣體效應(yīng)對螺旋槽干氣密封性能影響的數(shù)值分析[J].排灌機械工程學(xué)報,2015,33(10):874-881.SONG P Y,CHAN W,MAO W Y,et al.Numerical analysis on effect of real gas on spiral groove dry gas seal performance[J].Journal of Drainage and Irrigation Machinery Engineering,2015,33(10):874-881.

    [5]王衍,孫見君,陶凱,等.T型槽干氣密封數(shù)值分析及槽型優(yōu)化[J].摩擦學(xué)學(xué)報,2014,34(4):420-427.WANG Y,SUN J J,TAO K,et al.Numerical analysis of T-groove dry gas seal and groove optimization[J].Tribology,2014,34(4):420-427.

    [6]ZUK J,LUDWIG L P,JOHNSON R L.Compressible flow across shaft face seals[C]//Fifth International Conference on Fluid Sealing.Coventry.England:BHRA,1971:Paper H6.

    [7]ARGHIR M,DEFAYE C,FRENE J.The Lomakin effect in annular gas seals under choked flow conditions[J].Journal of Engineering for Gas Turbines & Power,2007,129(4):1028-1034.

    [8]THOMAS S,BRUNETIERE N,TOUMERIE B.Thermoelastohydrodynamic behavior of mechanical gas face seals operating at high pressure[J].Journal of Tribology,2007,129(4):841-850.

    [9]馬春紅,白少先,彭旭東,等.螺旋槽端面微間隙高速氣流潤滑密封特性[J].摩擦學(xué)學(xué)報,2015,35(6):699-706.MA C H,BAI S X,PENG X D,et al.Properties of high speed airflow lubrication in micro-clearance of spiral-groove face seals[J].Tribology,2015,35(6):699-706.

    [10]THATTE A,ZHENG X.Hydrodynamics and sonic flow transition in dry gas seals[C]//Proceedings of the ASME Turbo Expo.Fairfield,United States:American Society of Mechanical Engineers,2014.

    [11]江錦波.高速干氣密封端面型槽仿生設(shè)計與實驗研究[D].杭州:浙江工業(yè)大學(xué),2016.JIANG J B.Theoretical and experimental study of the bionic design of grooved surface of a high speed dry gas seal[D].Hangzhou:Zhejiang University of Technology,2016.

    [12]MILLER B A,GREEN I.Semi-analytical dynamic analysis of spiral-grooved mechanical gas face seals[J].Journal of Tribology,2003,125(2):403-413.

    [13]宋鵬云,胡曉鵬,許恒杰,等.實際氣體對T槽干氣密封動態(tài)特性的影響[J].化工學(xué)報,2014,65 (4):1344-1352.SONG P Y,HU X P,XU H J,et al.Effect of real gas on dynamic T-groove dry gas seal[J].CIESC Journal,2014,65 (4):1344-1352.

    [14]WANG W,LIU Y,JIANG P.Numerical investigation on influence of real gas properties on nonlinear behavior of labyrinth seal-rotor system[J].Applied Mathematics & Computation,2015,263:12-24.

    [15]劉雨川.端面氣膜密封特性研究[D].北京:北京航天航空大學(xué),2000.LIU Y C.Study on the characteristics of gas film face seal[D].Beijing:Beihang University,2000.

    [16]CHEN H,ZHENG J,XU P,et al.Study on real-gas equations of high pressure hydrogen[J].International Journal of Hydrogen Energy,2010,35(7):3100-3104.

    [17]蘇長蓀.高等工程熱力學(xué)[M].北京:高等教育出版社,1987:361-362.SU C S.Advanced Thermodynamics[M].Beijing:Higher Education Press,1987:361-362.

    [18]劉暉.實際氣體溫度絕熱指數(shù)和容積絕熱指數(shù)的計算[J].石油化工高等學(xué)校學(xué)報,2000,13(4):42-45.LIU H.Calculation of the isentropic temperature change exponent and isentropic volume change exponent of real gas[J].Journal of Petrochemical Universities,2000,13(4):42-45.

    [19]鄧成香,宋鵬云,馬愛琳.干氣密封的實際氣體焦耳-湯姆遜效應(yīng)分析[J].化工學(xué)報,2016,67(9):3833-3842.DENG C X,SONG P Y,MA A L.Analysis of Joule-Thomson effect of real gas system sealed by dry gas[J].CIESC Journal,2016,67(9):3833-3842.

    [20]許恒杰,宋鵬云.三自由度微擾下的靜壓干氣密封動態(tài)特性分析[J].排灌機械工程學(xué)報,2017,35(1):56-64.XU H J,SONG P Y.Analysis on dynamic characteristics of aerostatic dry gas seals with three degrees of freedom perturbation[J].Journal of Drainage and Irrigation Machinery Engineering,2017,35(1):56-64.

    [21]LIU Y,SHEN X,XU W.Numerical analysis of dynamic coefficients for gas film face seals[J].Journal of Tribology,2002,124(4):743-754.

    [22]劉雨川,徐萬孚,王之櫟,等.端面氣膜密封動力特性系數(shù)的計算[J].清華大學(xué)學(xué)報(自然科學(xué)版),2002,42(2):185-189.LIU Y C,XU W F,WANG Z L,et al.Dynamic coefficients for gas film face seal[J].J.Tsinghua Univ.(Sci.& Technol.),2002,42(2):185-189.

    [23]黃平,許蘭貴,孟永鋼,等.求解磁頭/磁盤超薄氣膜潤滑性能的有效有限差分算法[J].機械工程學(xué)報,2007,43(3):43-49.HUANG P,XU L G,MENG Y G,et al.Effective finite difference method to calculate lubricating performances of ultra-thin gas film of magnetic head/disk[J].Chinese Journal of Mechanical Engineering,2007,43(3):43-49.

    [24]LEBECK A O.Principles and Design of Mechanical Face Seals[M].John Wiley & Sons,Inc.,1991:371.

    [25]RUAN B.Numerical modeling of dynamic sealing behaviors of spiral groove gas face seals[J].Journal of Tribology,2002,124(1):186-195.

    [26]RUAN B.A semi-analytical solution to the dynamic tracking of non-contacting gas face seals[J].Journal of Tribology,2002,124(1):196-202.

    [27]陳源,彭旭東,李紀(jì)云,等.螺旋槽結(jié)構(gòu)參數(shù)對干氣密封動態(tài)特性的影響研究[J].摩擦學(xué)學(xué)報,2016,36(4):397-405.CHEN Y,PENG X D,LI J Y,et al.The influence of structure parameters of spiral groove on dynamic characteristics of dry gas seal[J].Tribology,2016,36(4):397-405.

    [28]劉向鋒,徐辰,黃偉峰.基于半解析法的極端工況干氣密封動態(tài)特性研究與參數(shù)設(shè)計[J].清華大學(xué)學(xué)報(自然科學(xué)版),2014,54(2):223-228.LIU X F,XU C,HUANG W F.Analysis and parametric design of the dynamics of a dry gas seal for extreme operating conditions using a semi-analytical method[J].J.Tsinghua Univ.(Sci.& Technol.),2014,54(2):223-228.

    [29]ZIRKELBACK N.Parametric study of spiral groove gas face seals[J].Tribology Transactions,2000,43(2):337-343.

    [30]晏鑫,蔣玉娥,李軍,等.迷宮密封轉(zhuǎn)子動力學(xué)特性的數(shù)值模擬[J].熱能動力工程,2009,24(5):566-570.YAN X,JIANG Y E,LI J,et al.Numerical calculation of dynamic coefficients for gas film cylinder seal[J].Journal of Engineering for Thermal Energy and Power,2009,24(5):566-570.

    date:2017-05-31.

    Prof.SONG Pengyun,songpengyunkm@sina.com

    supported by the National Natural Science Foundation of China (51465026).

    Dynamic characteristics of spiral groove dry gas seals with consideration of hydrogen real gas and choked flow effects

    XU Hengjie1,2,SONG Pengyun1,MAO Wenyuan1,DENG Qiangguo1
    (1Faculty of Chemical Engineering,Kunming University of Science and Technology,Kunming650500,Yunnan,China;2Faculty of Environmental Science and Engineering,Kunming University of Science and Technology,Kunming650500,Yunnan,China)

    The exit pressure boundary was determined by Chen’s real gas equation of hydrogen and choked flow condition of gas exit speed reaching to the sound speed.Then,dynamic characteristics of spiral groove dry gas seal (S-DGS) at various operating parameters was analyzed by perturbation method,and compared to those of ideal gas and coercive exit pressure boundary models.The results show that real gas and choked flow effects should be taken into account for studying dynamic characteristics of S-DGS at high pressure.Both effects reduced direct stiffness coefficients of hydrogen S-DGS but increased direct damping coefficients.The influence of these two effects on dynamic gas film stiffness gradually enhanced with the increase of squeeze number and inlet pressure.In addition,when frequency ratio was varied,these two effects had a significant influence on gas film dynamic stiffness but minimal influence on gas film dynamic damping.Compared to models of ideal gas and coercive pressure boundary condition at the studied operating circumstance,these two effects caused three gas film dynamic damping coefficients (Czz,Cαα,Cαβ) by mean standard deviations of 2.28%,1.93%,2.79% when squeeze number is variable and 4.08%,2.07%,1.82% when inlet pressure is variable.

    spiral groove dry gas seal; dynamic characteristics; real gas; choked flow; numerical analysis

    S 277.9;TH 136

    A

    0438—1157(2017)12—4675—10

    10.11949/j.issn.0438-1157.20170704

    2017-05-31收到初稿,2017-07-07收到修改稿。

    聯(lián)系人:宋鵬云。

    許恒杰(1989—),男,博士研究生。

    國家自然科學(xué)基金項目(51465026)。

    猜你喜歡
    干氣氣膜端面
    KDF3E成型機濾棒端面觸頭的原因及排除方法
    T 型槽柱面氣膜密封穩(wěn)態(tài)性能數(shù)值計算研究
    高溫熔鹽泵干氣螺旋密封性能的研究
    氣膜孔堵塞對葉片吸力面氣膜冷卻的影響
    靜葉柵上游端壁雙射流氣膜冷卻特性實驗
    火箭推進(2020年2期)2020-05-06 02:53:56
    銅基合金襯套端面鍍鉻質(zhì)量的改善
    優(yōu)化吸收穩(wěn)定單元操作
    化工管理(2017年36期)2018-01-04 03:26:13
    躲避霧霾天氣的氣膜館
    老舊端面磨齒機故障處理
    降低干氣中C3含量的技術(shù)措施
    化工管理(2015年21期)2015-05-28 12:12:56
    嫩草影院精品99| 麻豆精品久久久久久蜜桃| 白带黄色成豆腐渣| 午夜老司机福利剧场| 亚洲中文字幕一区二区三区有码在线看| 亚洲婷婷狠狠爱综合网| 在线播放无遮挡| 欧美3d第一页| 国产私拍福利视频在线观看| 天堂av国产一区二区熟女人妻| 亚洲国产最新在线播放| 国产午夜精品久久久久久一区二区三区| 欧美丝袜亚洲另类| 亚洲综合色惰| 午夜精品一区二区三区免费看| 嫩草影院精品99| 麻豆久久精品国产亚洲av| 色视频www国产| 插逼视频在线观看| 亚洲成人精品中文字幕电影| 啦啦啦啦在线视频资源| 国产男人的电影天堂91| 三级国产精品片| 少妇猛男粗大的猛烈进出视频 | 亚洲18禁久久av| 日韩欧美在线乱码| 六月丁香七月| 免费黄网站久久成人精品| 国产视频内射| 青青草视频在线视频观看| www.色视频.com| 91午夜精品亚洲一区二区三区| 免费人成在线观看视频色| 欧美丝袜亚洲另类| 欧美成人一区二区免费高清观看| 亚洲国产日韩欧美精品在线观看| 国内精品宾馆在线| 国产私拍福利视频在线观看| 国产毛片a区久久久久| 最近的中文字幕免费完整| 国产高清视频在线观看网站| 亚洲经典国产精华液单| 午夜激情欧美在线| 色吧在线观看| 国产精品国产高清国产av| 一区二区三区免费毛片| 国产一区二区三区av在线| 99热这里只有是精品50| 亚洲高清免费不卡视频| 哪个播放器可以免费观看大片| 十八禁国产超污无遮挡网站| 中国国产av一级| 亚洲色图av天堂| 精品人妻一区二区三区麻豆| 一级毛片久久久久久久久女| 精品人妻视频免费看| 在线播放国产精品三级| 一边摸一边抽搐一进一小说| 成人毛片a级毛片在线播放| 99久久人妻综合| 国内精品宾馆在线| 男人舔奶头视频| 久久久久久久久中文| 国产亚洲av嫩草精品影院| 99热网站在线观看| 欧美不卡视频在线免费观看| 少妇人妻精品综合一区二区| 亚洲自偷自拍三级| 女人久久www免费人成看片 | 大话2 男鬼变身卡| 国产精品永久免费网站| 菩萨蛮人人尽说江南好唐韦庄 | 日韩人妻高清精品专区| 成人综合一区亚洲| 日韩亚洲欧美综合| 99久久中文字幕三级久久日本| 黄色配什么色好看| 精品午夜福利在线看| 国产精品一区二区在线观看99 | 欧美高清性xxxxhd video| 五月玫瑰六月丁香| 在线播放无遮挡| 一级黄片播放器| 综合色av麻豆| 欧美日韩综合久久久久久| 免费看日本二区| 一级二级三级毛片免费看| 亚洲精品,欧美精品| 国产一区亚洲一区在线观看| 国产老妇女一区| 熟妇人妻久久中文字幕3abv| 国产成人精品一,二区| 免费看日本二区| 精品一区二区三区视频在线| 69av精品久久久久久| 亚洲欧洲日产国产| 免费av毛片视频| av免费观看日本| 日日啪夜夜撸| 1000部很黄的大片| 精品国产露脸久久av麻豆 | 99热网站在线观看| 亚洲精品日韩在线中文字幕| 亚洲av二区三区四区| 成人av在线播放网站| 亚洲性久久影院| 久久精品国产鲁丝片午夜精品| 水蜜桃什么品种好| 丰满少妇做爰视频| 深爱激情五月婷婷| 九九热线精品视视频播放| 国产欧美另类精品又又久久亚洲欧美| 成人国产麻豆网| 成人国产麻豆网| 高清日韩中文字幕在线| 成人三级黄色视频| 国产精品福利在线免费观看| 国产老妇女一区| 精品少妇黑人巨大在线播放 | 亚洲经典国产精华液单| 精品少妇黑人巨大在线播放 | 高清av免费在线| 人人妻人人澡欧美一区二区| 熟女电影av网| 夫妻性生交免费视频一级片| 特级一级黄色大片| 亚洲成av人片在线播放无| 看片在线看免费视频| 国产在线一区二区三区精 | 尤物成人国产欧美一区二区三区| 在线免费观看不下载黄p国产| 最近手机中文字幕大全| 免费搜索国产男女视频| 舔av片在线| 国产精品一区二区三区四区久久| 男女下面进入的视频免费午夜| 九九久久精品国产亚洲av麻豆| 国产精品一及| 在线免费观看的www视频| www.色视频.com| 亚洲欧美成人综合另类久久久 | 最近中文字幕高清免费大全6| 超碰av人人做人人爽久久| 视频中文字幕在线观看| 长腿黑丝高跟| 成人无遮挡网站| 女的被弄到高潮叫床怎么办| 亚洲av中文字字幕乱码综合| 成人亚洲欧美一区二区av| 汤姆久久久久久久影院中文字幕 | 国产精品伦人一区二区| 日本黄色片子视频| 精品99又大又爽又粗少妇毛片| 日韩一区二区视频免费看| 欧美最新免费一区二区三区| 1000部很黄的大片| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲国产精品成人综合色| 在线免费十八禁| 亚洲欧美日韩高清专用| 婷婷色av中文字幕| 久久久久久久久久久丰满| 高清毛片免费看| 国产精品麻豆人妻色哟哟久久 | 国产精品一二三区在线看| 成年免费大片在线观看| 老司机福利观看| 亚洲av不卡在线观看| 一级黄片播放器| 秋霞在线观看毛片| 亚洲欧美一区二区三区国产| 亚洲国产精品成人久久小说| 99在线视频只有这里精品首页| 亚洲精品国产av成人精品| 能在线免费看毛片的网站| 免费看美女性在线毛片视频| 一夜夜www| 久久6这里有精品| 超碰av人人做人人爽久久| 欧美zozozo另类| 久久久久九九精品影院| 亚洲国产最新在线播放| 午夜老司机福利剧场| 欧美成人一区二区免费高清观看| av卡一久久| 好男人视频免费观看在线| 久久精品综合一区二区三区| 国产av不卡久久| 久久精品夜夜夜夜夜久久蜜豆| 91午夜精品亚洲一区二区三区| 国产毛片a区久久久久| 在线播放国产精品三级| 国产av一区在线观看免费| 亚洲人成网站高清观看| 亚洲人成网站高清观看| 国产色婷婷99| 大又大粗又爽又黄少妇毛片口| 国产精品野战在线观看| 久久久久久久久久久免费av| 日本免费在线观看一区| 汤姆久久久久久久影院中文字幕 | 成人美女网站在线观看视频| 国产成人一区二区在线| 熟女电影av网| 欧美一区二区精品小视频在线| 欧美日本亚洲视频在线播放| 热99在线观看视频| 国产精品人妻久久久影院| 99国产精品一区二区蜜桃av| av国产久精品久网站免费入址| 国产乱来视频区| 亚洲人成网站在线播| 黄色欧美视频在线观看| 日韩,欧美,国产一区二区三区 | 久久久久久久久大av| 天堂√8在线中文| 国国产精品蜜臀av免费| 国产免费视频播放在线视频 | 日韩亚洲欧美综合| 99热这里只有是精品在线观看| 精品人妻一区二区三区麻豆| 亚洲丝袜综合中文字幕| 国产精品人妻久久久久久| 国内精品一区二区在线观看| 欧美区成人在线视频| 汤姆久久久久久久影院中文字幕 | 一级爰片在线观看| 国产激情偷乱视频一区二区| 伦理电影大哥的女人| 在线观看av片永久免费下载| 丰满人妻一区二区三区视频av| av视频在线观看入口| 日韩精品青青久久久久久| 日产精品乱码卡一卡2卡三| 在线观看一区二区三区| 黑人高潮一二区| 丰满乱子伦码专区| 午夜亚洲福利在线播放| 91久久精品国产一区二区成人| 毛片女人毛片| 国产美女午夜福利| 国产一级毛片在线| 色尼玛亚洲综合影院| 男女那种视频在线观看| 欧美人与善性xxx| 插逼视频在线观看| 一区二区三区高清视频在线| 自拍偷自拍亚洲精品老妇| 亚洲美女视频黄频| 99在线人妻在线中文字幕| 少妇裸体淫交视频免费看高清| 97热精品久久久久久| 非洲黑人性xxxx精品又粗又长| 极品教师在线视频| 日韩视频在线欧美| 岛国毛片在线播放| 99视频精品全部免费 在线| 天天躁夜夜躁狠狠久久av| 热99在线观看视频| 久久久久久国产a免费观看| 日本免费a在线| 禁无遮挡网站| 中文字幕精品亚洲无线码一区| av在线老鸭窝| 看十八女毛片水多多多| 国产精品.久久久| 神马国产精品三级电影在线观看| av免费在线看不卡| 日本av手机在线免费观看| 午夜爱爱视频在线播放| 最近最新中文字幕大全电影3| 国产精品一区二区三区四区免费观看| 亚洲av日韩在线播放| 亚洲美女视频黄频| 日韩强制内射视频| 日韩av在线大香蕉| 成人鲁丝片一二三区免费| 可以在线观看毛片的网站| 在线a可以看的网站| 五月伊人婷婷丁香| 亚洲av一区综合| 美女cb高潮喷水在线观看| 纵有疾风起免费观看全集完整版 | 五月玫瑰六月丁香| 有码 亚洲区| 国产亚洲精品久久久com| 国产又黄又爽又无遮挡在线| 伊人久久精品亚洲午夜| 国产免费福利视频在线观看| 久久精品夜夜夜夜夜久久蜜豆| 狂野欧美激情性xxxx在线观看| 嫩草影院入口| 国产片特级美女逼逼视频| 丰满少妇做爰视频| 丝袜喷水一区| 精品国产露脸久久av麻豆 | 国产伦精品一区二区三区四那| 免费观看在线日韩| 少妇丰满av| 国产高清有码在线观看视频| 免费在线观看成人毛片| 免费人成在线观看视频色| 亚洲第一区二区三区不卡| 亚洲一区高清亚洲精品| 精品久久久久久久人妻蜜臀av| 国产精品一区二区性色av| 日韩中字成人| 超碰av人人做人人爽久久| 婷婷色综合大香蕉| 久久国产乱子免费精品| 国产精品三级大全| 亚洲国产最新在线播放| 老司机福利观看| 国产精品人妻久久久久久| 精品一区二区三区人妻视频| 亚洲怡红院男人天堂| 天美传媒精品一区二区| 最近2019中文字幕mv第一页| 日本一本二区三区精品| .国产精品久久| 久久久成人免费电影| 国语对白做爰xxxⅹ性视频网站| 91久久精品国产一区二区成人| 色播亚洲综合网| 久久久国产成人精品二区| 日韩在线高清观看一区二区三区| 日韩人妻高清精品专区| or卡值多少钱| 超碰97精品在线观看| 久久久欧美国产精品| 永久免费av网站大全| 国产私拍福利视频在线观看| 成人午夜高清在线视频| 亚洲精品一区蜜桃| 网址你懂的国产日韩在线| 国内揄拍国产精品人妻在线| av在线播放精品| 欧美另类亚洲清纯唯美| 一级黄色大片毛片| 我的女老师完整版在线观看| 少妇猛男粗大的猛烈进出视频 | 亚洲三级黄色毛片| 白带黄色成豆腐渣| 国产三级在线视频| 亚洲av成人av| 好男人视频免费观看在线| 国产久久久一区二区三区| 免费无遮挡裸体视频| 国产真实伦视频高清在线观看| 搡女人真爽免费视频火全软件| 国产高清有码在线观看视频| 一级黄片播放器| 欧美日本视频| 中文欧美无线码| 听说在线观看完整版免费高清| 水蜜桃什么品种好| 国产成人午夜福利电影在线观看| 亚洲精品aⅴ在线观看| 天堂√8在线中文| 午夜免费激情av| av黄色大香蕉| 黄色日韩在线| av播播在线观看一区| 久久久精品欧美日韩精品| 国产精品嫩草影院av在线观看| 最后的刺客免费高清国语| 99久国产av精品| 国产一区有黄有色的免费视频 | 日韩欧美 国产精品| 日韩欧美精品v在线| 亚洲av电影在线观看一区二区三区 | 久久99热这里只有精品18| 亚洲综合精品二区| 国产91av在线免费观看| 日本黄大片高清| 美女黄网站色视频| 午夜福利视频1000在线观看| 国产色爽女视频免费观看| 国产熟女欧美一区二区| 中文字幕精品亚洲无线码一区| 18禁动态无遮挡网站| 高清视频免费观看一区二区 | 成年版毛片免费区| 中文天堂在线官网| АⅤ资源中文在线天堂| 亚洲国产色片| 免费观看精品视频网站| 成人毛片a级毛片在线播放| 亚洲精品国产成人久久av| 免费看日本二区| h日本视频在线播放| 国产成人91sexporn| 亚洲综合精品二区| 黄色欧美视频在线观看| 久久欧美精品欧美久久欧美| 国产精品女同一区二区软件| 亚洲一区高清亚洲精品| 精品酒店卫生间| 免费不卡的大黄色大毛片视频在线观看 | 黄色日韩在线| 亚洲熟妇中文字幕五十中出| 午夜爱爱视频在线播放| 中文字幕久久专区| 精品人妻偷拍中文字幕| av在线天堂中文字幕| 欧美极品一区二区三区四区| 女人久久www免费人成看片 | 久久人人爽人人片av| 亚洲,欧美,日韩| 久久久久精品久久久久真实原创| 偷拍熟女少妇极品色| 久久精品国产99精品国产亚洲性色| 内射极品少妇av片p| 26uuu在线亚洲综合色| 亚洲va在线va天堂va国产| 99热这里只有是精品50| 国产黄色小视频在线观看| 精品欧美国产一区二区三| 国产精品电影一区二区三区| 日本-黄色视频高清免费观看| 日日撸夜夜添| 青春草国产在线视频| 99久久九九国产精品国产免费| 日韩精品有码人妻一区| 午夜福利视频1000在线观看| 久久人妻av系列| 日本黄大片高清| 嫩草影院精品99| 欧美成人免费av一区二区三区| 99视频精品全部免费 在线| 九九热线精品视视频播放| 韩国av在线不卡| 日韩国内少妇激情av| 深夜a级毛片| 好男人视频免费观看在线| 国产日韩欧美在线精品| 国产黄片美女视频| 中文精品一卡2卡3卡4更新| 亚洲国产最新在线播放| 国产乱人偷精品视频| 全区人妻精品视频| 精品久久国产蜜桃| 中文字幕精品亚洲无线码一区| 少妇熟女aⅴ在线视频| 日本wwww免费看| 国产人妻一区二区三区在| 人人妻人人看人人澡| 国产 一区精品| 精品欧美国产一区二区三| 国产精品一区二区在线观看99 | 99热这里只有精品一区| 亚洲国产精品合色在线| 日本五十路高清| 插逼视频在线观看| 97热精品久久久久久| 免费在线观看成人毛片| 黄色一级大片看看| 亚洲精品色激情综合| 人妻系列 视频| 精品国产一区二区三区久久久樱花 | 能在线免费看毛片的网站| 国产淫语在线视频| 哪个播放器可以免费观看大片| 国产精华一区二区三区| 在线播放无遮挡| 精品一区二区免费观看| 精品人妻熟女av久视频| 看免费成人av毛片| 国产欧美另类精品又又久久亚洲欧美| 免费黄色在线免费观看| 成人高潮视频无遮挡免费网站| 亚洲av.av天堂| 99热全是精品| 小说图片视频综合网站| 亚洲欧美精品自产自拍| 日日撸夜夜添| 男人狂女人下面高潮的视频| 成人毛片a级毛片在线播放| 两个人视频免费观看高清| 午夜精品一区二区三区免费看| 亚洲欧美精品专区久久| 欧美成人一区二区免费高清观看| 成人午夜精彩视频在线观看| 精品国产三级普通话版| 国产乱来视频区| 国语自产精品视频在线第100页| 一边亲一边摸免费视频| 晚上一个人看的免费电影| 国产av码专区亚洲av| 嫩草影院精品99| 国产亚洲av片在线观看秒播厂 | 国产精品久久电影中文字幕| 欧美一区二区精品小视频在线| 一个人观看的视频www高清免费观看| 成人美女网站在线观看视频| 深夜a级毛片| 国产片特级美女逼逼视频| 亚洲国产精品合色在线| 九色成人免费人妻av| 欧美高清性xxxxhd video| 纵有疾风起免费观看全集完整版 | 日本免费一区二区三区高清不卡| 99久久无色码亚洲精品果冻| 赤兔流量卡办理| 亚洲综合色惰| 女人被狂操c到高潮| 永久免费av网站大全| 亚洲高清免费不卡视频| 日日摸夜夜添夜夜添av毛片| 免费电影在线观看免费观看| 欧美3d第一页| 国产免费福利视频在线观看| 丝袜美腿在线中文| 小蜜桃在线观看免费完整版高清| 色综合亚洲欧美另类图片| 桃色一区二区三区在线观看| 精品国产露脸久久av麻豆 | 中文字幕av成人在线电影| av国产久精品久网站免费入址| 中文资源天堂在线| av又黄又爽大尺度在线免费看 | 人人妻人人看人人澡| 99久久中文字幕三级久久日本| 国产黄色视频一区二区在线观看 | 国产精品女同一区二区软件| 国产日韩欧美在线精品| 啦啦啦啦在线视频资源| 日韩欧美 国产精品| 国模一区二区三区四区视频| 亚洲最大成人中文| 国产一区二区亚洲精品在线观看| 中文字幕人妻熟人妻熟丝袜美| 日本色播在线视频| 精品无人区乱码1区二区| 99视频精品全部免费 在线| 日本av手机在线免费观看| 麻豆一二三区av精品| 国产精品一及| 亚洲av电影不卡..在线观看| 只有这里有精品99| 熟女人妻精品中文字幕| 国产精品久久久久久久电影| av在线播放精品| 色5月婷婷丁香| a级毛片免费高清观看在线播放| 真实男女啪啪啪动态图| 久久精品国产自在天天线| 99热这里只有是精品在线观看| 午夜福利在线观看免费完整高清在| av视频在线观看入口| 国产亚洲精品av在线| 变态另类丝袜制服| 国产乱人偷精品视频| 日本午夜av视频| 日韩一区二区三区影片| 男的添女的下面高潮视频| 亚洲国产最新在线播放| 精品酒店卫生间| 我要搜黄色片| 亚洲精华国产精华液的使用体验| 亚洲av福利一区| 亚洲天堂国产精品一区在线| 日韩精品有码人妻一区| 欧美3d第一页| 国产av在哪里看| 简卡轻食公司| 国产免费福利视频在线观看| 最后的刺客免费高清国语| 免费无遮挡裸体视频| 午夜亚洲福利在线播放| 日日摸夜夜添夜夜添av毛片| 九色成人免费人妻av| 男女那种视频在线观看| 国产精品一及| 亚洲国产精品国产精品| 最新中文字幕久久久久| 一级av片app| 久久99蜜桃精品久久| 久久久久性生活片| 亚洲精品日韩在线中文字幕| 2021少妇久久久久久久久久久| 精品不卡国产一区二区三区| 两个人视频免费观看高清| 高清视频免费观看一区二区 | 九九久久精品国产亚洲av麻豆| 18禁动态无遮挡网站| 亚洲自拍偷在线| 免费在线观看成人毛片| 中文字幕av成人在线电影| av女优亚洲男人天堂| 亚洲av日韩在线播放| 国产在视频线精品| 人妻夜夜爽99麻豆av| 成人性生交大片免费视频hd| 大话2 男鬼变身卡| 男人的好看免费观看在线视频| 欧美日韩综合久久久久久| 亚洲高清免费不卡视频| 亚洲av二区三区四区| 国产精品久久视频播放| 国产精品1区2区在线观看.| 一边摸一边抽搐一进一小说| 全区人妻精品视频| 国产亚洲av嫩草精品影院| 在线a可以看的网站| 亚洲五月天丁香| 欧美日韩综合久久久久久| 大香蕉97超碰在线| 国产精品一区二区性色av| 国内少妇人妻偷人精品xxx网站| 国内精品宾馆在线| 国产精品一区二区性色av| 非洲黑人性xxxx精品又粗又长| 亚洲真实伦在线观看| 18禁动态无遮挡网站|