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

袋型阻尼密封動力學(xué)特性雙控制體Bulk Flow模型

2024-07-17 00:00:00桂佳強(qiáng)李志剛李軍
西安交通大學(xué)學(xué)報 2024年7期

收稿日期:2023-12-08。

作者簡介:桂佳強(qiáng)(1998—),男,碩士生;李志剛(通信作者),男,教授,博士生導(dǎo)師。

基金項目:國家自然科學(xué)基金資助項目(52106153)。

網(wǎng)絡(luò)出版時間:2024-02-29""" 網(wǎng)絡(luò)出版地址:https:∥link.cnki.net/urlid/61.1069.T.20240228.1705.002

摘要:為快速準(zhǔn)確預(yù)測袋型阻尼密封泄漏特性和動力學(xué)特性,針對傳統(tǒng)單控制體Bulk Flow模型預(yù)測精度低、無法預(yù)測交叉動力系數(shù)的問題,提出了袋型阻尼密封雙控制體Bulk Flow模型和動力學(xué)特性數(shù)值預(yù)測方法,并開發(fā)了計算程序。首先,依據(jù)邊界層理論,將袋型密封腔室劃分為兩個控制體,推導(dǎo)了控制體的連續(xù)性、周向動量和能量方程,引入Swamee-Jain和Takahashi方程,計算流體-壁面間和流體-流體間的周向黏性摩擦力;其次,采用牛頓-拉夫森算法和攝動分析法分別求解0階和1階控制方程,獲得各剛度、阻尼動力特性系數(shù);然后,通過與袋型阻尼密封泄漏量和動力特性系數(shù)的實驗值、單控制體Bulk Flow模型和非定常計算流體動力學(xué)(CFD)數(shù)值結(jié)果進(jìn)行比較,驗證了模型和方法的準(zhǔn)確性和可靠性;最后,研究了轉(zhuǎn)子轉(zhuǎn)速(10000、15000、20000r/min)和預(yù)旋比(0.067、0.724、0.997)對袋型阻尼密封動力學(xué)特性的影響。結(jié)果表明:所發(fā)展的模型和方法具有計算速度快、預(yù)測精度高(泄漏量預(yù)測誤差小于6%,動力特性系數(shù)預(yù)測誤差小于38%)的優(yōu)點;轉(zhuǎn)子轉(zhuǎn)速和進(jìn)口預(yù)旋的增大均會導(dǎo)致袋型阻尼密封有效阻尼顯著減小,穿越頻率顯著增大,易誘發(fā)軸系失穩(wěn)。

關(guān)鍵詞:袋型阻尼密封;泄漏特性;動力學(xué)特性;雙控制體;Bulk Flow模型

中圖分類號:TK262" 文獻(xiàn)標(biāo)志碼:A

DOI:10.7652/xjtuxb202407003" 文章編號:0253-987X(2024)07-0026-13

A Two-Control-Volume Bulk Flow Model for Rotordynamic

Characteristics of Pocket Damper Seals

GUI Jiaqiang, LI Zhigang, LI Jun

(Institute of Turbomachinery, Xi’an Jiaotong University, Xi’an 710049, China)

Abstract:To predict the leakage and rotordynamic characteristics of pocket damper seals (PDS) quickly and accurately, a two-control-volume (CV) Bulk Flow model and numerical method are proposed for the prediction of the rotordynamic characteristics of pocket damper seals, and the calculation program is also developed. This new model is expected to solve the problems of the traditional one-control-volume Bulk Flow model, such as the low accuracy and the inability to predict cross-coupled coefficients. Firstly, based on the boundary layer theory, each pocket cavity of the pocket damper seal is divided into two control volumes, and the continuity equations, circumferential momentum equations, and energy equations are derived for the two control volumes. Swamee-Jain and Takahashi formulas are introduced to calculate the circumferential friction between the fluid and the wall and between the fluid and the fluid. Secondly, Newton-Raphson method and perturbation analysis method are used to solve the zeroth-order and first-order governing equations, and the stiffness and damping rotordynamic characteristic coefficients are obtained. Then, the accuracy and reliability of the proposed model and method are verified by comparing with the experimental data of leakage and dynamic characteristic coefficients of pocket damper seals, along with the numeric results of the one-control-volume Bulk Flow model and the unsteady CFD (Computational Fluid Dynamics). Finally, the influences of rotor speed (10000, 15000, and 20000r/min) and inlet preswirl ratio (0.067, 0.724, and 0.997) on the dynamic characteristics of pocket damper seals are calculated and analyzed. Results show that the model and method proposed have the advantages of fast calculation speed and high prediction accuracy (leakage prediction error less than 6%, and dynamic coefficients prediction error" less than 38%); The increasing rotor speed and inlet preswirl result in a significant decrease in the effective damping and an obvious increase in the crossover frequency, increasing the serious risk of rotor instability.

Keywords:pocket damper seals; leakage characteristics; rotordynamic characteristics; two control volume; Bulk Flow model

旋轉(zhuǎn)動密封是現(xiàn)代透平機(jī)械的關(guān)鍵部件之一。先進(jìn)的動密封技術(shù)不僅能顯著地降低透平機(jī)械的泄漏損失從而提高透平機(jī)械的運行效率,而且還能顯著地提高透平機(jī)械的軸系穩(wěn)定性。迷宮密封因其結(jié)構(gòu)簡單、技術(shù)成熟和制造維護(hù)成本低的特點,成為透平機(jī)械中應(yīng)用最為廣泛的一種旋轉(zhuǎn)動密封。但是,迷宮密封因其周向貫通的腔室結(jié)構(gòu),密封腔內(nèi)部周向旋流速度大,會產(chǎn)生顯著的流體激振力,存在著轉(zhuǎn)子穩(wěn)定性差的隱患[1]。Vance和Schultz[2]提出了一種名為TAMSEAL的新型阻尼密封,后改名為袋型阻尼密封(PDS)。針對袋型阻尼密封泄漏特性的實驗和數(shù)值研究[3-5]已證明,袋型阻尼密封腔室內(nèi)的周向擋板能有效地減小周向旋流速度。文獻(xiàn)[6-8]研究表明:相比于迷宮密封,袋型阻尼密封不僅能更有效地減少,泄漏量,還能增大阻尼,提高轉(zhuǎn)子的穩(wěn)定性。近年來,袋型阻尼密封在高參數(shù)(高壓、高轉(zhuǎn)速)透平機(jī)械中受到了越來越多的關(guān)注和應(yīng)用,因此發(fā)展能快速準(zhǔn)確預(yù)測袋型阻尼密封泄漏量和動力特性的方法對于高性能袋型阻尼密封設(shè)計和先進(jìn)透平機(jī)械動力學(xué)分析具有重要的工程應(yīng)用價值。

目前,針對旋轉(zhuǎn)動密封泄漏特性和動力學(xué)特性的研究方法主要有:①實驗測量法;②計算流體動力學(xué)(CFD)數(shù)值研究法;③Bulk Flow模型分析法。實驗測量法存在著實驗條件難以實現(xiàn)和耗費資金大的問題。雖然非定常CFD數(shù)值研究法具有預(yù)測精度高的優(yōu)勢,但存在計算耗時長、成本高的缺點。Bulk Flow模型分析法具有預(yù)測速度極快(lt;3min)、對計算資源要求低且預(yù)測精度滿足工程要求的優(yōu)點,常用來進(jìn)行密封快速選型和結(jié)構(gòu)優(yōu)化設(shè)計,是一種重要的密封工業(yè)設(shè)計和性能分析方法。為了滿足現(xiàn)代透平機(jī)械日益提高的運行參數(shù)、效率與軸系穩(wěn)定性需求,推廣先進(jìn)袋型阻尼密封的應(yīng)用,開展更精確、快速的袋型阻尼密封Bulk Flow模型分析法研究至關(guān)重要。

Iwataubo[9]于1980年首次提出了適用于迷宮密封的單控制體(1-CV)Bulk Flow模型,模型包括連續(xù)性方程和周向動量方程。Childs等[10]在Iwataubo的基礎(chǔ)上考慮了轉(zhuǎn)子渦動引起的腔室周向截面積變化,并由此預(yù)測了直通式迷宮密封的動力學(xué)特性。Scharrer[11]改進(jìn)了適用于直通式迷宮密封的雙控制體Bulk Flow模型,考慮了渦流速度對動力學(xué)特性的影響,相較于Childs的單控制體Bulk Flow模型,交叉剛度和直接阻尼的預(yù)測精度更高。Picardo和Childs[12]將單控制體和雙控制體模型的預(yù)測結(jié)果與迷宮密封高壓實驗結(jié)果進(jìn)行對比發(fā)現(xiàn):對于迷宮密封,單控制體Bulk Flow模型的總體預(yù)測精度比雙控制體Bulk Flow模型高。Cangioli等[13]在Childs的基礎(chǔ)上引入了能量方程,通過假設(shè)泄漏過程絕熱來分析流體焓變的影響,提高了Bulk Flow模型在負(fù)預(yù)旋工況下的迷宮密封動力特性系數(shù)預(yù)測精度。王天昊等[14]在Cangioli的基礎(chǔ)上,對不同經(jīng)驗公式組合成的72種泄漏模型進(jìn)行了適用性分析,獲得了最佳泄漏模型,并利用該模型研究了壓比和預(yù)旋比對直通式迷宮密封動力學(xué)特性的影響規(guī)律。胡樂豪等[15]基于數(shù)值計算結(jié)果,針對迷宮密封提出了一種新的摩擦系數(shù)預(yù)測模型,能夠提高Bulk Flow對黏性摩擦力的預(yù)測精度。

不同于迷宮密封周向貫通的密封腔室,袋型阻尼密封的密封腔內(nèi)均勻布置有周向擋板結(jié)構(gòu),形成周向孤立的袋型密封腔。針對袋型阻尼密封,Ertas[16]在Iwataubo[9]的迷宮密封Bulk Flow模型基礎(chǔ)上,結(jié)合了Alford[17]提出的一維軸向流動模型,提出了適用于袋型阻尼密封的單控制體Bulk Flow模型來預(yù)測動力學(xué)特性,但該模型忽略了流體的周向流速,不考慮周向動量方程,導(dǎo)致無法預(yù)測交叉剛度系數(shù)。李志剛[18]在文獻(xiàn)[16]的基礎(chǔ)上對不同經(jīng)驗公式組合成的9種泄漏模型進(jìn)行了適用性分析,獲得了最佳泄漏模型,并利用該模型研究了幾何參數(shù)對袋型阻尼密封動力學(xué)特性的影響規(guī)律。由于袋型阻尼密封的傳統(tǒng)單控制體Bulk Flow模型忽略了密封腔內(nèi)的周向流動,無法評估轉(zhuǎn)速和進(jìn)口預(yù)旋等關(guān)鍵運行參數(shù)的影響,所以對密封動力特性系數(shù)預(yù)測誤差大,無法滿足袋型阻尼密封結(jié)構(gòu)設(shè)計和性能評估要求。李志剛等[19]針對袋型阻尼密封動力特性預(yù)測,提出了基于多頻渦動模型的非定常CFD數(shù)值預(yù)測方法,提高了動力特性系數(shù)的預(yù)測精度,但非定常CFD數(shù)值計算成本極高,無法滿足工程應(yīng)用需求。

為解決傳統(tǒng)袋型阻尼密封單控制體Bulk Flow模型預(yù)測精度低、不能預(yù)測交叉動力系數(shù),以及無法評估轉(zhuǎn)速和預(yù)旋等關(guān)鍵運行參數(shù)影響的問題,本文提出了一種新型的袋型阻尼密封雙控制體Bulk Flow模型和動力學(xué)特性數(shù)值預(yù)測方法,并開發(fā)了計算程序,不僅考慮了由于周向擋板導(dǎo)致的腔室內(nèi)流體自循環(huán)流動,也考慮了周向擋板間隙處的流體繞轉(zhuǎn)子的周向流動。采用該模型的預(yù)測程序研究了轉(zhuǎn)子轉(zhuǎn)速和預(yù)旋比對袋型阻尼密封動力學(xué)特性的影響規(guī)律。

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

1.1" 密封腔控制體

圖1是文獻(xiàn)[18]用CFD模擬出的袋型阻尼密封腔內(nèi)三維流線分布圖,圖2和圖3分別是其對應(yīng)的軸向和周向截面流線圖。流體穿過密封軸向齒和周向擋板與轉(zhuǎn)子面的徑向間隙后形成射流,從而將腔室劃分為兩個區(qū)域,一個是貼近轉(zhuǎn)子表面的射流區(qū)域,另一個是腔室內(nèi)的自循環(huán)流動區(qū)域。本文采用兩個控制體(CV1和CV2)分別對兩個區(qū)域流動進(jìn)行數(shù)學(xué)表征,采用平均射流厚度定義控制體CV1與CV2的分界面。由于軸向的射流厚度遠(yuǎn)遠(yuǎn)小于周向的射流厚度,所以忽略軸向射流對控制體劃分界面位置的影響,最終劃分結(jié)果的周向截面示意圖和三維示意圖如圖4和圖5所示。

假設(shè)射流厚度沿周向的變化遵循在半不定平板上的湍流邊界層公式[20],控制體CV1在徑向上的高度可表示為

H1i=A1iRα=∫Rα0δ(ξ)dξRα=

∫Rα0(Hdi+0.37ξ/Re1/5ξ)dξRα(1)

Reξ=V1iξvi (2)

式中:H1i為軸向第i個腔室的控制體CV1的徑向高度;A1i為軸向第i個腔室的控制體CV1的周向截面積;R為轉(zhuǎn)子半徑;α為單個腔室所占的周向角度;δ(ξ)為射流在坐標(biāo)ξ處的徑向厚度;Hdi為軸向第i個腔室的周向擋板間隙;Reξ為射流在坐標(biāo)ξ處的雷諾數(shù);V1i為軸向第i個腔室的控制體CV1內(nèi)流體的周向流速;vi為軸向第i個腔室內(nèi)流體的運動黏度。

1.2" 模型假設(shè)

為簡化控制方程和求解過程,需對密封內(nèi)的流動做如下假設(shè):

(1)每個控制體周向截面上流體壓力均勻分布;

(2)控制體CV1與CV2之間無質(zhì)量交換;

(3)控制體CV2內(nèi)流體軸向自循環(huán)速度U2i和周向自循環(huán)速度V2i為恒定值,不受轉(zhuǎn)子渦動影響;

(4)控制體CV2內(nèi)流體壓力P2i為恒定值,不受轉(zhuǎn)子渦動影響,且與控制體CV1內(nèi)流體壓力穩(wěn)態(tài)值相等;

(5)每個腔室的控制體CV1和CV2內(nèi)流體的焓值相等,并且為恒定值,不受轉(zhuǎn)子渦動影響;

(6)不考慮流體的慣性;

(7)每個腔室的聲共振頻率遠(yuǎn)高于轉(zhuǎn)子轉(zhuǎn)速的頻率;

(8)轉(zhuǎn)子渦動的軌道偏心距比密封齒徑向間隙小得多,流體激振力和渦動位移滿足線性關(guān)系;

(9)計算中采用實際氣體性質(zhì),流動過程為絕熱過程;

(10)工質(zhì)為單相氣體。

1.3 "控制方程

基于Childs[10]的迷宮密封單控制體Bulk Flow模型,考慮動力特性系數(shù)的頻率相關(guān)性,推導(dǎo)出圖5中控制體CV1的連續(xù)方程和周向動量方程如下

t(ρ1iA1i)+ρ1iA1iV1iR+i+1-i=0(3)

t(ρ1iA1iV1i)+ρ1iA1iV21iR+i+1V1i-iV1i-1=

-A1iRP1i+τr1iar1i-τs1ias1i-τf1iaf1i(4)

式中:ρ1i和P1i分別為軸向第i個腔室的控制體CV1內(nèi)流體的密度和壓力;i為軸向第i個密封齒間隙處單位周長上的質(zhì)量流率;τf1i為軸向第i個腔室的控制體CV1內(nèi)流體作用于軸向第i個腔室的控制體CV2內(nèi)流體的周向切應(yīng)力;τr1i和τs1i分別為轉(zhuǎn)子面和定子面作用于軸向第i個腔室的控制體CV1內(nèi)流體的周向切應(yīng)力;af1i、ar1i和as1i分別為周向截面上τf1i、τr1i和τs1i作用于流體的長度,可表示為

af1i=Li (5)

ar1i=Li (6)

as1i=2(H1i-Hi) (7)

式中:Li為軸向第i個腔室的軸向長度;Hi為軸向第i個密封齒間隙。

假設(shè)控制體CV1與CV2之間無質(zhì)量交換,控制體CV2內(nèi)連續(xù)性方程(質(zhì)量守恒)是恒成立的,控制體CV2內(nèi)流體的物性不受轉(zhuǎn)子渦動影響,推導(dǎo)出圖5中控制體CV2的周向動量方程如下

τf1iaf1i-2τf2iaf2i-τs2ias2i=0 (8)

式中:τf2i為軸向第i個腔室的控制體CV2內(nèi)流體自循環(huán)產(chǎn)生的周向切應(yīng)力;τs2i為定子面作用于軸向第i個腔室的控制體CV2內(nèi)流體的周向切應(yīng)力;af2i和as2i分別為周向截面上τf2i和τs2i作用于流體的長度,計算方式如下

af2i=Li (9)

as2i=Li+2B-as1i (10)

其中B為腔室深度。

由于流體的焓不受轉(zhuǎn)子渦動的影響[13],因此能量方程僅考慮零階形式。假設(shè)控制體CV1與CV2內(nèi)流體壓力與焓的穩(wěn)態(tài)值相等,可基于文獻(xiàn)[13]的迷宮密封單控制體Bulk Flow模型,針對整個腔室推導(dǎo)出能量方程如下

ihi+V21i2-i+1hi-1+V21i-12=τr1iar1iRΩ (11)

式中:hi為軸向第i個腔室內(nèi)流體的焓;Ω為轉(zhuǎn)子轉(zhuǎn)動角速度。

1.4" 泄漏模型

準(zhǔn)確可靠的泄漏模型不僅可以幫助更準(zhǔn)確地計算出穩(wěn)態(tài)下各腔室的壓力,而且還可以提高密封流體激振力和動力特性系數(shù)的預(yù)測精度。文獻(xiàn)[18]證明了迷宮密封與袋型阻尼密封具有相同的泄漏特性,適用于迷宮密封的泄漏預(yù)測模型同樣適用于袋型阻尼密封。本文采用適用于實際氣體的廣義Neumann迷宮密封泄漏方程[21]作為袋型阻尼密封的泄漏模型,表達(dá)式如下

i=μiCfiHiP1i-1ρ1i-1-P1iρ1i (12)

式中:μi為Neumann[22]動能輸運系數(shù),在軸向第1個密封齒間隙處值為1,在其他密封齒間隙處表達(dá)式如下

μi=N(1-Ji)N+Ji (13)

Ji=1-1+16.6HiLi-1-2 (14)

式中:N為密封軸向齒個數(shù)。

式(12)中的Cfi為Chaplygin[23]流量系數(shù),其表達(dá)式如下

Cfi=ππ+2-5si+2s2i (15)

si=P1i-1P1iγ-1γ-1 (16)

式中:γ為工質(zhì)的絕熱指數(shù)。

當(dāng)密封進(jìn)出口壓比增大到一定值時,流體在軸向最后一密封齒間隙處出口速度會達(dá)到聲速(ζ),形成阻塞流,此時泄漏量不受密封出口側(cè)壓力影響,式(12)在i=N時需替換成如下

i=CfNHNρ1N-1ζ (17)

式中:ζ參考流體性質(zhì)數(shù)據(jù)庫取值[24]。流體在第i個密封齒間隙處的出口速度可表示為如下

U1i=iCfiHiρ1i-1 (18)

1.5" 切應(yīng)力求解

轉(zhuǎn)子和定子作用在流體上的流體-壁面黏性摩擦周向切應(yīng)力可表示為

τr1i=ρ1i2fr1i(RΩ-V1i)U21i+(RΩ-V1i)2 (19)

τs1i=ρ1i2fs1iV1iU21i+V21i (20)

τs2i=ρ2i2fs2iV2iU22i+V22i (21)

式中:ρ2i為軸向第i個腔室的控制體CV2內(nèi)流體的密度;fr1i、fs1i和fs2i為達(dá)西摩擦系數(shù);U2i為軸向第i個腔室的控制體CV2內(nèi)流體軸向自循環(huán)速度,Scharrer[25]用邊界層理論對其進(jìn)行過估計,表達(dá)式如下

U2i=0.206U1i (22)

本文采用Takahashi切應(yīng)力公式[26]定義流體-流體之間的黏性摩擦力,表達(dá)式如下

τf1i=12λfρ1i(U1i-U2i)2+(V1i-V2i)2(V1i-V2i)(23)

τf2i=12λfρ2i4U22i+4V22i2V2i (24)

λf=0.054 (25)

本文采用Swamee-Jain公式[27]對fr1i、fs1i、fs2i進(jìn)行估計,表達(dá)式如下

fr1i=116lger3.7Dh1i+5.74Re0.9r1i-2 (26)

fs1i=116lges3.7Dh1i+5.74Re0.9s1i-2 (27)

fs2i=116lges3.7Dh2i+5.74Re0.9s2i-2 (28)

式中:er、es分別為轉(zhuǎn)子、定子表面粗糙度;Dh1i、Dh2i分別為周向流動的水力直徑,Rer1i、Res1i、Res2i為雷諾數(shù),表示為如下

Dh1i=4H1iLi2(H1i+Li) (29)

Dh2i=4H2iLi2(H2i+Li) (30)

Rer1i=Dh1iU21i+(RΩ-V1i)2νi (31)

Res1i=Dh1iU21i+V21iνi (32)

Res2i=Dh2iU22i+V22iνi (33)

其中H2i為軸向第i個腔室的控制體CV2的徑向高度,可表示為如下

H2i=Hi+B-H1i (34)

1.6" 零階方程求解

在給定進(jìn)出口壓力、溫度和轉(zhuǎn)速邊界條件下,假設(shè)轉(zhuǎn)子軸心與密封中心重合且轉(zhuǎn)子不渦動(穩(wěn)態(tài)),可計算出各控制體CV1和CV2內(nèi)流體各參數(shù)的穩(wěn)態(tài)值,這些穩(wěn)態(tài)值可用于一階分析中。基于上述假設(shè),可將式(3)、(4)、(7)、(11)簡化為如下形式

01=…=0i=…=0N=0 (35)

0V01i-0V01i-1=τ0r1iar1i-τ0s1ias1i-τ0f1iaf1i (36)

τ0f1iaf1i-2τ0f2iaf2i-τ0s2ias2i=0 (37)

0h0i+V201i2-0h0i-1+V201i-12=τ0r1iar1iRΩ (38)

式中:下標(biāo)0代表穩(wěn)態(tài)值。

上述運算得到的方程為非線性隱式方程,采用牛頓-拉夫森算法迭代求解,采用相對容差為10-6的收斂準(zhǔn)則。

1.7" 一階方程求解

為了求解式(3)、(4)組成的非線性多元偏微分方程組,采用攝動分析法將方程組線性化。如圖6所示,假設(shè)轉(zhuǎn)子軸心繞定子中心沿橢圓軌跡渦動,此時熱力學(xué)變量和運動學(xué)變量(這里用符號φ表示)可近似為泰勒級數(shù)的前兩項,表示如下

φi=φ0i+φ1i(t,) (39)

式中:φ0i為穩(wěn)態(tài)項;φ1i(t,)為攝動(一階)項。

根據(jù)上述物理假設(shè),Bulk Flow變量的表達(dá)式可總結(jié)為

Pi=P0i+P1i(t,)

Vi=V0i+V1i(t,)

hi=h0i

Hi=H0i+H1(t,) (40)

將式(40)代入式(3)、(4)并去除穩(wěn)態(tài)項,可以得到一階攝動方程組。

一階攝動方程組是描述轉(zhuǎn)子渦動時密封腔室內(nèi)流體的壓力攝動項、周向旋流速度攝動項和密封間隙攝動項之間關(guān)系的二維非穩(wěn)態(tài)方程組,需要采用以下變換進(jìn)行求解。

如圖6所示,轉(zhuǎn)子軸心渦動軌跡表達(dá)式如下

x=acosωt; y=bsin ωt (41)

式中:a和b分別為橢圓渦動軌跡的長半軸和短半軸長度;ω為轉(zhuǎn)子渦動角速度。

根據(jù)轉(zhuǎn)子運動規(guī)律,密封間隙攝動函數(shù)可表示為如下復(fù)數(shù)形式

H1(t,)=-acosωtcos-bsin ωtsin =

Re-a+b2ej(-ωt)-a-b2ej(+ωt) (42)

控制體CV1內(nèi)流體的壓力和周向流速的攝動函數(shù)可表示為相似的形式

P11i(t,)=ReP+11iej(+ωt)+P-11iej(-ωt) (43)

V11i(t,)=ReV+11iej(+ωt)+V-11iej(-ωt) (44)

式中:P+11i、P-11i、V+11i和V-11i是常量復(fù)數(shù)系數(shù)。

將式(42)~(44)代入一階攝動方程組并消去ej(+ωt)和ej(-ωt)項,可以得到復(fù)變量線性方程組如下

A-1iP-1i-1P+1i-1V-1i-1V+1i-1+A0iP-1iP+1iV-1iV+1i+A+1iP-1i+1P+1i+1V-1i+1V+1i+1=Bi ""(45)

式中:A-1i、A0i和A+1i為四階復(fù)常數(shù)系數(shù)矩陣;Bi為四維復(fù)常數(shù)向量。上述復(fù)常數(shù)均由零階方程的解決定。

對于具有N個密封齒的袋型阻尼密封,聯(lián)立N-1個密封腔室的復(fù)變量線性方程組式(45),得到4(N-1)階復(fù)變量線性方程組并求解,將結(jié)果代入式(43)、(44),可以得到各控制體CV1內(nèi)流體的壓力和周向流速的攝動函數(shù)。

1.8" 動力特性系數(shù)的提取

根據(jù)動密封小位移渦動理論,轉(zhuǎn)子做小位移渦動時,袋型阻尼密封轉(zhuǎn)子面受到的流體激振力(Fx,F(xiàn)y)、轉(zhuǎn)子渦動位移(x,y)和渦動速度(,)滿足力-位移線性方程[28]

-Fx(t)Fy(t)=Kk-kKx(t)y(t)+Cc-cC(t)(t) (46)

式中:K、k、C、c分別為直接剛度系數(shù)、交叉剛度系數(shù)、直接阻尼系數(shù)和交叉阻尼系數(shù)。

將式(41)代入式(46)并寫作復(fù)數(shù)形式

F=Fx+jFy=-a+b2[(K+cω)-j(k-

Cω)]ejωt-a-b2[(K-cω)-j(k+Cω)]e-jωt(47)

作用于轉(zhuǎn)子的流體激振力F還可通過將壓力攝動函數(shù)和表面切應(yīng)力沿轉(zhuǎn)子周向表面積分,得到

F=Fx+jFy=-∑N-1i=1LiR∫2π0(P11i-jariτ1r1i)ejd(48)

式中:τ1r1i為作用于轉(zhuǎn)子表面的切應(yīng)力攝動量。

將壓力攝動函數(shù)P11i和切應(yīng)力攝動量τ1r1i的表達(dá)式代入式(48),整理得

F=-a+b2Z-ejωt-a-b2Z+e-jωt (49)

Z+=∑N-1i=1πLiR1-jariτr1iP1i2+11ia-b-

jariτr1iV1i2+11ia-b+jariτr1iHi (50)

Z-=∑N-1i=1πLiR1-jariτr1iP1i2-11ia+b-

jariτr1iV1i2-11ia+b+jariτr1iHi (51)

式中:+11i、-11i、+11i、-11i分別與P+11i、P-11i、V+11i、V-11i共軛。

結(jié)合式(47)和式(49),可得袋型阻尼密封動力特性系數(shù)如下

K=12ReZ++Z-; k=-12ImZ++Z-

C=-12ωImZ+-Z-; c=-12ωReZ+-Z-(52)

2" 模型驗證

本文將Delgado[29]8齒貫穿擋板袋型阻尼密封實驗件的密封泄漏量、動力特性實驗結(jié)果作為驗證算例,并將本文雙控制體Bulk Flow模型計算結(jié)果與Ertas[16]的單控制體Bulk Flow模型和Delgado[29]的非定常CFD計算結(jié)果進(jìn)行比較。表1和表2分別給出了袋型阻尼密封的實驗幾何參數(shù)和運行工況參數(shù)。圖7給出了袋型阻尼密封的定子面結(jié)構(gòu)及其上半部分的正視圖,上面標(biāo)注了表1和表2工況2中的一些數(shù)據(jù)。

由于傳統(tǒng)單控制體Bulk Flow模型僅能預(yù)測袋型阻尼密封的直接剛度系數(shù)K和直接阻尼系數(shù)C,所以僅能將Ertas的單控制體Bulk Flow模型與本文雙控制體Bulk Flow模型的K、C進(jìn)行對比。

圖8和圖9給出了在表2所示工況2下,袋型阻尼密封的直接剛度系數(shù)K和直接阻尼系數(shù)C隨渦動頻率f的變化曲線。對于袋型阻尼密封的K和C,本文雙控制體Bulk Flow模型預(yù)測值與實驗值符合良好,具有與非定常CFD數(shù)值方法相近的、明顯優(yōu)于傳統(tǒng)單控制體Bulk Flow模型的預(yù)測精度。

如圖8所示,本文雙控制體Bulk Flow模型能夠準(zhǔn)確預(yù)測袋型阻尼密封K隨渦動頻率f的變化趨勢。在低頻區(qū)(flt;150Hz),預(yù)測值偏大;在高頻區(qū)(fgt;150Hz),預(yù)測值與實驗值吻合良好,平均預(yù)測誤差小于38%,明顯優(yōu)于非定常CFD預(yù)測結(jié)果。如圖9所示,本文雙控制體Bulk Flow模型的C預(yù)測值在整個頻率范圍內(nèi)均與實驗值吻合良好,平均預(yù)測誤差小于12%,與非定常CFD的平均預(yù)測誤差(小于14%)相近。

圖10給出了在表2所示工況2下,袋型阻尼密封k隨f的變化曲線。本文雙控制體Bulk Flow模型的交叉剛度系數(shù)k預(yù)測值在整個頻率范圍內(nèi)均與實驗值吻合良好,平均預(yù)測誤差小于22%,在實驗測量誤差范圍內(nèi),優(yōu)于非定常CFD平均預(yù)測誤差(小于42%)。

為評估袋型阻尼密封動力特性對轉(zhuǎn)子系統(tǒng)穩(wěn)定性的影響,引入了動密封有效阻尼系數(shù)Ceff,其表達(dá)式如下

Ceff=C-kω (54)

圖11給出了在表2所示工況2下,袋型阻尼密封有效阻尼系數(shù)Ceff隨渦動頻率f的變化曲線。本文雙控制體Bulk Flow模型能夠準(zhǔn)確預(yù)測袋型阻尼密封有效阻尼系數(shù)Ceff隨渦動頻率f的變化趨勢,總體上預(yù)測值偏大,平均預(yù)測誤差小于17%,與非定常CFD的平均預(yù)測誤差小于23%相近。

表3給出了表2所示3個工況下的袋型阻尼密封泄漏量。本文雙控制體Bulk Flow模型的泄漏量預(yù)測值與實驗值吻合良好,平均預(yù)測誤差小于6%,與非定常CFD的平均預(yù)測誤差小于1%相近。

綜合以上分析,本文所發(fā)展的雙控制體Bulk Flow模型及動力特性數(shù)值預(yù)測方法能夠快速、準(zhǔn)確預(yù)測袋型阻尼密封的泄漏量和動力特性系數(shù),具有與非定常CFD數(shù)值預(yù)測方法相近的預(yù)測精度(泄漏量預(yù)測誤差小于6%,動力特性系數(shù)預(yù)測誤差小于38%),同時具有與傳統(tǒng)Bulk Flow模型同等的計算速度(小于2min)。

3" 動力學(xué)特性影響因素研究

傳統(tǒng)單控制體Bulk Flow模型無法評估轉(zhuǎn)速和預(yù)旋這兩個關(guān)鍵運行參數(shù)對袋型阻尼密封動力學(xué)特性的影響。本文采用所發(fā)展的雙控制體Bulk Flow模型和計算程序,分析了表1所示的袋型阻尼密封在不同工況下的動力特性系數(shù)對渦動頻率的變化情況,并且與Thiele[30]的實驗值進(jìn)行了比較。

3.1" 轉(zhuǎn)子轉(zhuǎn)速的影響

為了研究轉(zhuǎn)子轉(zhuǎn)速對袋型阻尼密封動力特性的影響,本節(jié)計算了袋型阻尼密封在3種轉(zhuǎn)子轉(zhuǎn)速(10000、15000、20000r/min)下的動力特性系數(shù)。計算采用的進(jìn)口壓力Pin=70bar、出口壓力Pout=35bar、預(yù)旋比λ=0.060。

圖12給出了不同轉(zhuǎn)速下袋型阻尼密封動力特性系數(shù)的雙控制體Bulk Flow模型預(yù)測值與實驗測量結(jié)果隨渦動頻率f的變化。如圖12所示,通過與實驗值的比較可知,本文雙控制體Bulk Flow模型能準(zhǔn)確預(yù)測轉(zhuǎn)子轉(zhuǎn)速對各動力特性系數(shù)以及隨渦動頻率f的變化規(guī)律的影響。袋型阻尼密封直接剛度隨轉(zhuǎn)子轉(zhuǎn)速的增大而增大(轉(zhuǎn)速從10000r/min增大到20000r/min,K增大了46%~92%),在高頻區(qū)的增長更為明顯,隨渦動頻率的增大而增大,在高頻區(qū)的頻率相關(guān)性更為顯著;交叉剛度和直接阻尼均隨轉(zhuǎn)子轉(zhuǎn)速的增大而增大(轉(zhuǎn)速從10000r/min增大到20000r/min,k增大了204%~296%,C增大了12%~19%),二者均與頻率弱相關(guān);直接阻尼在整個渦動頻率范圍內(nèi)均為正值,且隨轉(zhuǎn)子轉(zhuǎn)速的增大略有增大,但有效阻尼隨轉(zhuǎn)子轉(zhuǎn)速的增大而顯著減小,尤其是低頻區(qū)(flt;80Hz),這主要是由于轉(zhuǎn)子轉(zhuǎn)速顯著地增大了密封交叉剛度;有效阻尼的穿越頻率fc隨轉(zhuǎn)子轉(zhuǎn)速的增大而增大(轉(zhuǎn)速從10000r/min 增大到20000r/min,fc從20Hz增大到80Hz),轉(zhuǎn)子系統(tǒng)穩(wěn)定工作的渦動頻率范圍隨轉(zhuǎn)子轉(zhuǎn)速的增大而降低。因此,轉(zhuǎn)子轉(zhuǎn)速導(dǎo)致袋型阻尼密封性能惡化,易誘發(fā)轉(zhuǎn)子失穩(wěn)。

根據(jù)周向動力方程式(42),控制體CV1內(nèi)的周向速度對腔室動態(tài)壓力具有顯著影響。轉(zhuǎn)子轉(zhuǎn)速的變化改變了袋型阻尼密封各控制體CV1內(nèi)的周向旋流速度,進(jìn)而影響到密封的動力學(xué)特性。圖13給出了不同轉(zhuǎn)子轉(zhuǎn)速下袋型阻尼密封的控制體CV1內(nèi)周向旋流速度沿軸向的分布。轉(zhuǎn)子轉(zhuǎn)速增大導(dǎo)致控制體CV1內(nèi)的周向旋流速度等比例增大(轉(zhuǎn)速從10000r/min增大到20000r/min,周向旋流速度增大了120% ~134%)。

3.2" 預(yù)旋比的影響

為研究預(yù)旋比對袋型阻尼密封動力特性的影響,計算分析了袋型阻尼密封在3種預(yù)旋比(λ=0.067,0.724,0.997)下的動力特性系數(shù)。密封進(jìn)口壓力Pin=70bar、出口壓力Pout=45.5bar、轉(zhuǎn)子轉(zhuǎn)速n=10000r/min。

預(yù)旋比λ的定義為氣體進(jìn)口預(yù)旋速度Uin與轉(zhuǎn)子表面線速度之比,表達(dá)式如下

λ=60Uin2πnR (54)

式中:n為轉(zhuǎn)子轉(zhuǎn)速。

圖14給出了不同預(yù)旋比下本文模型預(yù)測值與實驗測量結(jié)果隨渦動頻率f的變化。由圖可知,本文模型能準(zhǔn)確預(yù)測預(yù)旋比對各動力特性系數(shù)以及隨渦動頻率f的變化規(guī)律。袋型阻尼密封直接剛度隨預(yù)旋比的增大而增大(預(yù)旋比從0.067增大到0.997時,K增大了22%~87%),在高頻區(qū)的增長更為明顯,與頻率弱相關(guān)。交叉剛度和直接阻尼均隨預(yù)旋比的增大而增大(預(yù)旋比從0.067增大到0.997 時,k增大了258%~354%,C增大了5%~7%),二者均與頻率弱相關(guān);直接阻尼在整個渦動頻率范圍內(nèi)均為正值,且隨預(yù)旋比的增大略有增大,但有效阻尼隨預(yù)旋比的增大而顯著減小,尤其是低頻區(qū)(flt;120Hz),這主要是由于進(jìn)口預(yù)旋顯著增大了密封交叉剛度;有效阻尼的穿越頻率fc隨預(yù)旋比

的增大而增大(預(yù)旋比從0.067增大到0.997時,fc從30Hz增大到120Hz),轉(zhuǎn)子系統(tǒng)穩(wěn)定工作的渦動頻率范圍隨預(yù)旋比的增大而顯著減小。因此,進(jìn)口預(yù)旋導(dǎo)致袋型阻尼密封阻尼性能惡化,易誘發(fā)轉(zhuǎn)子失穩(wěn),需采用止旋裝置有效抑制密封進(jìn)口預(yù)旋。

預(yù)旋比的變化改變了袋型阻尼密封各控制體CV1內(nèi)的周向旋流速度,進(jìn)而影響到密封的動力學(xué)特性。圖15給出了不同預(yù)旋比下袋型阻尼密封的控制體CV1內(nèi)周向旋流速度沿軸向的分布。進(jìn)口預(yù)旋顯著增大上游腔室周向旋流,而對下游腔室影響微弱;因此,密封腔止旋裝置安裝于上游密封腔或密封進(jìn)口止旋效果最優(yōu)。

3" 結(jié)" 論

本文提出了一種新型的袋型阻尼密封雙控制體Bulk Flow模型和動力學(xué)特性數(shù)值預(yù)測方法,并開發(fā)了計算程序。相比于傳統(tǒng)的單控制體Bulk Flow模型和非定常CFD數(shù)值預(yù)測方法,本文提出的雙控制體Bulk Flow模型和數(shù)值預(yù)測方法具有計算速度快、預(yù)測精度高(泄漏量預(yù)測誤差小于6%,動力特性系數(shù)預(yù)測誤差小于38%,計算時長小于2min)的優(yōu)點。

轉(zhuǎn)子轉(zhuǎn)速和進(jìn)口預(yù)旋的增大均會導(dǎo)致袋型阻尼密封有效阻尼顯著減小,穿越頻率顯著增大,易誘發(fā)軸系失穩(wěn)。

轉(zhuǎn)子轉(zhuǎn)速增大使各密封腔室周向旋流等比例顯著增大;進(jìn)口預(yù)旋顯著增大上游腔室周向旋流,而對下游腔室影響微弱;因此,密封腔止旋裝置安裝于上游密封腔或密封進(jìn)口止旋效果最優(yōu)。

所發(fā)展的雙控制體Bulk Flow和動力特性數(shù)值預(yù)測方法可為袋型阻尼密封動力學(xué)特性的快速評估和結(jié)構(gòu)優(yōu)化設(shè)計提供可靠的技術(shù)手段。

參考文獻(xiàn):

[1]CHILDS D W. Turbomachinery rotordynamics: phenomena, modeling, and analysis [M]. New York, USA: Wiley, 1993.

[2]VANCE J M, SCHULTZ R R. A new damper seal for turbomachinery [C]//ASME 1993 Design Technical Conferences. New York, USA: ASME, 1993: 139-148.

[3]LI Zhigang, LI Jun, FENG Zhenping, et al. Numerical investigations on the leakage flow characteristics of pocket damper seals [C]//ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition. New York, USA: ASME, 2011: 701-712.

[4]李志剛, 李軍, 豐鎮(zhèn)平. 高轉(zhuǎn)速袋型阻尼密封泄漏的特性 [J]. 航空動力學(xué)報, 2012, 27(12): 2828-2835.

LI Zhigang, LI Jun, FENG Zhenping. Leakage flow of pocket damper seal at high rotational speed [J]. Journal of Aerospace Power, 2012, 27(12): 2828-2835.

[5]LI Zhigang, LI Jun, YAN Xin, et al. Numerical investigations on the leakage flow characteristics of pocket damper labyrinth seals [J]. Proceedings of the Institution of Mechanical Engineers: Part A" Journal of Power and Energy, 2012, 226(7): 932-948.

[6]LI Zhigang, LI Jun, FENG Zhenping. Numerical investigations on the leakage and rotordynamic characteristics of pocket damper seals: part Ⅰ" effects of pressure ratio, rotational speed, and inlet preswirl [J]. Journal of Engineering for Gas Turbines and Power, 2015, 137(3): 032503.

[7]LI Zhigang, LI Jun, FENG Zhenping. Numerical investigations on the leakage and rotordynamic characteristics of pocket damper seals: part Ⅱ" effects of partition wall type, partition wall number, and cavity depth [J]. Journal of Engineering for Gas Turbines and Power, 2015, 137(3): 032504.

[8]LI Zhigang, LI Jun, FENG Zhenping. Numerical comparison of rotordynamic characteristics for a fully partitioned pocket damper seal and a labyrinth seal with high positive and negative inlet preswirl [J]. Journal of Engineering for Gas Turbines and Power, 2016, 138(4): 042505.

[9]IWATSUBO T. Evaluation of instability forces of labyrinth seals in turbines or compressors[EB/OL]. (1980-01-01)[2023-12-15]. https://ntrs.nasa.gov/citations/19800021214.

[10]CHILDS D W, SCHARRER J K.An Iwatsubo-based solution for labyrinth seals: comparison to experimental results [J]. Journal of Engineering for Gas Turbines and Power, 1986, 108(2): 325-331.

[11]SCHARRER J K. Theory versus experiment for the rotordynamic coefficients of labyrinth gas seals: part Ⅰ" a two control volume model [J]. Journal of Vibration, Acoustics, Stress, and Reliability in Design, 1988, 110(3): 270-280.

[12]PICARDO A, CHILDS D W.Rotordynamic coefficients for a tooth-on-stator labyrinth seal at 70 bar supply pressures: measurements versus theory and comparisons to a hole-pattern stator seal [J]. Journal of Engineering for Gas Turbines and Power, 2005, 127(4): 843-855.

[13]CANGIOLI F, PENNACCHI P, VANNINI G, et al. Effect of energy equation in one control-volume bulk-flow model for the prediction of labyrinth seal dynamic coefficients [J]. Mechanical Systems and Signal Processing, 2018, 98: 594-612.

[14]王天昊, 李志剛, 李軍. 采用Bulk-Flow模型的直通式迷宮密封轉(zhuǎn)子動力特性研究 [J]. 西安交通大學(xué)學(xué)報, 2021, 55(5): 25-33.

WANG Tianhao, LI Zhigang, LI Jun. Investigation on the rotordynamic characteristics of straight-through labyrinth seal using bulk-flow model [J]. Journal of Xi’an Jiaotong University, 2021, 55(5): 25-33.

[15]胡樂豪, 鄧清華, 劉洲洋, 等. 超臨界二氧化碳迷宮密封內(nèi)摩擦損失數(shù)值研究及流動特性分析 [J]. 西安交通大學(xué)學(xué)報, 2023, 57(3): 68-78.

HU Lehao, DENG Qinghua, LIU Zhouyang, et al. Numerical investigation for windage loss and flow characteristics analysis with supercritical CO2 in labyrinth seal [J]. Journal of Xi’an Jiaotong University, 2023, 57(3): 68-78.

[16]ERTAS B H.Rotordynamic force coefficients of pocket damper seals [D]. College Station: Texas Aamp;M University, 2005.

[17]ALFORD J S. Protecting turbomachinery from self-excited rotor whirl [J]. Journal of Engineering for Power, 1965, 87(4): 333-343.

[18]李志剛. 袋型阻尼密封泄漏特性和轉(zhuǎn)子動力特性的研究 [D]. 西安: 西安交通大學(xué), 2013.

[19]李志剛, 李軍, 豐鎮(zhèn)平. 袋型阻尼密封轉(zhuǎn)子動力特性的多頻單向渦動預(yù)測模型 [J]. 西安交通大學(xué)學(xué)報, 2012, 46(5): 13-18.

LI Zhigang, LI Jun, FENG Zhenping. Multiple frequencies One-Dimensional prediction model of rotordynamic characteristics for pocket damper seals [J]. Journal of Xi’an Jiaotong University, 2012, 46(5): 13-18.

[20]SCHLICHTING H. Boundary-layer theory [M]. 7th ed. New York, USA: McGraw-Hill, 1979.

[21]CANGIOLI F, PENNACCHI P, RIBONI G, et al. Sensitivity analysis of the one-control volume bulk-flow model for a 14 teeth-on-stator straight-through labyrinth seal [C]//ASME Turbo Expo 2017: Turbomachinery Technical Conference and Exposition. New York, NY, USA: ASME, 2017: V07AT34A002.

[22]NEOMAHN K. Zur Frage der verwendung von durchblickdichtungen im dampfturbinenbau [J]. Maschinenbautechnik, 1964, 12(4): 188-195.

[23]GUREVICH M I. The theory of jets in an ideal fluid [M]. Oxford, UK: Pergamon Press, 1966.

[24]BELL I H, QUOILIN S, WRONSKI J, et al. CoolProp: an open-source reference-quality thermophysical property library [EB/OL]. [2023-11-20]. https://so1.linfen3.top/scholar?q=Coolprop%3A+An+open-source+reference-quality+thermophysical+property+library.

[25]SCHARRER J K. Theory versus experiment for the rotordynamic coefficients of labyrinth gas seals: part Ⅰa two control volume model [J]. Journal of Vibration, Acoustics, Stress, and Reliability in Design, 1988, 110(3): 270-280.

[26]TAKAHASHI N, MIURA H, NARITA M, et al. Development of scallop cut type damper seal for centrifugal compressors [J]. Journal of Engineering for Gas Turbines and Power, 2015, 137(3): 032509.

[27]CANGIOLI F, PENNACCHI P, VANNINI G, et al. On the thermodynamic process in the bulk-flow model for the estimation of the dynamic coefficients of labyrinth seals [J]. Journal of Engineering for Gas Turbines and Power, 2018, 140(3): 032502.

[28]晏鑫, 李軍, 豐鎮(zhèn)平. Bulk Flow方法分析孔型密封轉(zhuǎn)子動力特性的有效性 [J]. 西安交通大學(xué)學(xué)報, 2009, 43(1): 24-28.

YAN Xin, LI Jun, FENG Zhenping. Validation of two-control-volume bulk flow method for rotordynamic characteristics of hole-pattern seals [J]. Journal of Xi’an Jiaotong University, 2009, 43(1): 24-28.

[29]DELGADO A,SAN ANDRS L, YANG Jing, et al. Experimental force coefficients for a fully-partitioned pocket damper seal and comparison to other two seal types [J]. Journal of Engineering for Gas Turbines and Power, 2023, 145(5): 051019.

[30]THIELE J J.Dynamic characterization of fully partitioned damper seals [D]. College Station: Texas Aamp;M University, 2019.

(編輯" 杜秀杰)

主站蜘蛛池模板: 婷婷丁香色| 日韩天堂视频| 日韩 欧美 小说 综合网 另类| 四虎影视国产精品| 欧美中文字幕在线二区| 日韩毛片视频| 玖玖精品视频在线观看| 91久久精品日日躁夜夜躁欧美 | 亚洲成人www| 亚欧成人无码AV在线播放| 99久久精品久久久久久婷婷| 色成人综合| 视频二区国产精品职场同事| 国内精品免费| 亚洲爱婷婷色69堂| 亚洲香蕉伊综合在人在线| 国产成人久视频免费| 2022国产91精品久久久久久| 国产毛片基地| 国产激情影院| 99久久精品免费观看国产| 在线国产毛片| 欧美精品伊人久久| 国产女人18水真多毛片18精品 | 97在线公开视频| 呦女精品网站| 亚洲V日韩V无码一区二区| 亚洲男人天堂久久| 在线亚洲小视频| 国产国模一区二区三区四区| 99久久精品国产综合婷婷| www亚洲天堂| 国产精品视频系列专区| 亚洲精品视频网| 波多野结衣久久精品| 日韩欧美中文字幕在线韩免费| 青草国产在线视频| 无码精品福利一区二区三区| 一级毛片在线播放免费| 亚洲精品日产精品乱码不卡| 中文字幕人成人乱码亚洲电影| P尤物久久99国产综合精品| 亚洲天堂在线免费| 亚洲精品无码AⅤ片青青在线观看| 日韩精品久久无码中文字幕色欲| AV网站中文| 亚洲第一色视频| 欧美一级黄色影院| 中文字幕久久亚洲一区| 欧美一级高清片欧美国产欧美| 五月综合色婷婷| 热思思久久免费视频| 夜夜拍夜夜爽| 亚洲国产精品国自产拍A| 99这里只有精品在线| 午夜视频免费试看| 亚洲无码日韩一区| 午夜精品一区二区蜜桃| 国产老女人精品免费视频| 一级全免费视频播放| 国产成人综合日韩精品无码不卡| 99视频在线观看免费| 99在线视频精品| 黄片在线永久| 看国产一级毛片| 欧美国产综合视频| 国产精品色婷婷在线观看| 亚洲成人手机在线| 国产美女丝袜高潮| 国产亚洲高清视频| 色婷婷综合激情视频免费看| 精品人妻系列无码专区久久| 日本高清在线看免费观看| 国产精品亚洲αv天堂无码| 欧美综合在线观看| 国产毛片高清一级国语| 欧美国产成人在线| 精品国产中文一级毛片在线看| 人妻21p大胆| 成人一级免费视频| 999国内精品久久免费视频| 青青操国产视频|