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

    碰摩約束下柔性轉(zhuǎn)子模態(tài)特性及其計(jì)算方法

    2020-12-28 08:35:20于平超陳果王存楊默晗
    航空學(xué)報(bào) 2020年12期
    關(guān)鍵詞:模態(tài)系統(tǒng)

    于平超,陳果,王存,楊默晗

    1. 南京航空航天大學(xué) 民航學(xué)院,南京 211106 2. 北京動(dòng)力機(jī)械研究所,北京 100074

    由于航空發(fā)動(dòng)機(jī)對(duì)氣動(dòng)效率極高要求,其轉(zhuǎn)靜件間隙不斷減小,同時(shí)航空發(fā)動(dòng)機(jī)在全壽命周期內(nèi)遭遇的載荷工況又極為復(fù)雜,導(dǎo)致碰摩一直是航空發(fā)動(dòng)機(jī)中的常見(jiàn)現(xiàn)象。轉(zhuǎn)靜件碰摩會(huì)引發(fā)轉(zhuǎn)子劇烈振動(dòng)以及轉(zhuǎn)子軸系內(nèi)部的交變應(yīng)力,進(jìn)而導(dǎo)致轉(zhuǎn)子軸系的疲勞斷裂,危及轉(zhuǎn)子乃至整機(jī)的安全運(yùn)行,是影響航空發(fā)動(dòng)機(jī)結(jié)構(gòu)可靠性的關(guān)鍵問(wèn)題之一。

    當(dāng)碰摩發(fā)生時(shí),會(huì)對(duì)轉(zhuǎn)子產(chǎn)生沖擊、摩擦以及約束等多種物理效應(yīng)[1-2],進(jìn)而帶來(lái)轉(zhuǎn)子復(fù)雜的物理現(xiàn)象,為此國(guó)內(nèi)外學(xué)者開(kāi)展了廣泛的研究。Muszynska[3]對(duì)轉(zhuǎn)靜子局部碰摩問(wèn)題進(jìn)行了研究,表明系統(tǒng)存在階次為1/2、1/3、1/4的次諧波振動(dòng)。Chu和Zhang[4]對(duì)碰摩轉(zhuǎn)子的分岔和穩(wěn)定性進(jìn)行了分析,揭示了轉(zhuǎn)子動(dòng)力響應(yīng)從穩(wěn)定周期運(yùn)動(dòng)經(jīng)過(guò)倍周期分岔、擦邊分岔等變?yōu)閿M周期和混沌響應(yīng)的過(guò)程。Chen等[5]針對(duì)航空發(fā)動(dòng)機(jī)整機(jī)系統(tǒng)碰摩,建立了梁?jiǎn)卧麢C(jī)模型,分析了碰摩剛度比、不平衡量等參數(shù)對(duì)轉(zhuǎn)子振動(dòng)響應(yīng)、分岔以及混沌特性的影響。Ma等[6-8]建立了單點(diǎn)、多點(diǎn)以及全周碰摩形式的轉(zhuǎn)子模型,研究了不同碰摩形式對(duì)轉(zhuǎn)子系統(tǒng)時(shí)頻響應(yīng)和軸心軌跡的影響。

    碰摩除對(duì)轉(zhuǎn)子產(chǎn)生非線性激勵(lì)作用外,還對(duì)轉(zhuǎn)子產(chǎn)生約束效應(yīng)。約束效應(yīng),即轉(zhuǎn)子與靜子一直接觸或者持續(xù)間斷接觸而導(dǎo)致轉(zhuǎn)子系統(tǒng)剛度增加的現(xiàn)象,也稱剛化效應(yīng)[9]。利用該現(xiàn)象,Chu和Lu[10]建立了碰摩過(guò)程中轉(zhuǎn)子系統(tǒng)剛度的參數(shù)識(shí)別方法,并將其應(yīng)用于碰摩轉(zhuǎn)子的故障診斷。Wang等[11]分析了碰摩約束對(duì)轉(zhuǎn)子動(dòng)力學(xué)特性的影響,指出附加約束剛度使轉(zhuǎn)子共振區(qū)間擴(kuò)張,并出現(xiàn)不穩(wěn)定接觸區(qū)域。

    由于碰摩約束導(dǎo)致的轉(zhuǎn)子系統(tǒng)剛度變化,會(huì)使轉(zhuǎn)子系統(tǒng)模態(tài)頻率/振型改變,然而以往碰摩研究中很少關(guān)注碰摩轉(zhuǎn)子這類強(qiáng)非線性系統(tǒng)的模態(tài)特性。Hong等[12]結(jié)合Floquet理論和Hill方法建立了間歇碰摩轉(zhuǎn)子模態(tài)特性的計(jì)算方法,分析了碰摩約束剛度時(shí)變對(duì)轉(zhuǎn)子模態(tài)頻率和穩(wěn)定性的影響,但該方法需事先假設(shè)碰摩過(guò)程,局限性較大。陳艷華和江俊[13]首次將非線性模態(tài)理論引入到碰摩轉(zhuǎn)子中,通過(guò)解析方法推導(dǎo)獲得了碰摩轉(zhuǎn)子的模態(tài)頻率,其分析思路為研究碰摩約束對(duì)轉(zhuǎn)子模態(tài)的影響提供了借鑒。在此之上,Hong等[14]通過(guò)半解析的諧波平衡法建立了碰摩轉(zhuǎn)子模態(tài)的計(jì)算方法,分析了碰摩約束轉(zhuǎn)子的模態(tài)特征和穩(wěn)定性演化過(guò)程。然而,目前關(guān)于碰摩轉(zhuǎn)子模態(tài)的研究?jī)H局限于單盤(pán)Jeffcott轉(zhuǎn)子這類簡(jiǎn)化系統(tǒng),對(duì)于以航空發(fā)動(dòng)機(jī)轉(zhuǎn)子為代表的這類工程復(fù)雜柔性轉(zhuǎn)子,其不僅具有復(fù)雜質(zhì)量/剛度分布、陀螺效應(yīng)影響,同時(shí)工作轉(zhuǎn)速范圍內(nèi)含多階彎曲模態(tài)[15],碰摩約束對(duì)轉(zhuǎn)子模態(tài)特性的影響更加復(fù)雜。另一方面,工程中的轉(zhuǎn)子系統(tǒng)往往需要大量自由度進(jìn)行描述[16],這使得現(xiàn)有求解碰摩轉(zhuǎn)子非線性模態(tài)的解析或半解析方法無(wú)法勝任。

    鑒于此,本文在非線性模態(tài)理論框架下,研究具有真實(shí)結(jié)構(gòu)特征的航空發(fā)動(dòng)機(jī)柔性轉(zhuǎn)子在碰摩約束下的模態(tài)特性。首先,針對(duì)柔性轉(zhuǎn)子復(fù)雜結(jié)構(gòu)特征,基于梁?jiǎn)卧邢拊ń⒖紤]碰摩約束的轉(zhuǎn)子動(dòng)力學(xué)模型;其次,結(jié)合頻域中自由度壓縮、諧波平衡思想以及時(shí)頻轉(zhuǎn)換技術(shù)等提出一套適用于含大規(guī)模自由度復(fù)雜轉(zhuǎn)子非線性模態(tài)的求解方法;最后,詳細(xì)分析碰摩約束對(duì)實(shí)際航空發(fā)動(dòng)機(jī)低壓柔性轉(zhuǎn)子模態(tài)特性的影響規(guī)律。

    1 考慮碰摩約束的柔性轉(zhuǎn)子動(dòng)力學(xué)模型

    圖1所示為典型高涵道比渦扇發(fā)動(dòng)機(jī)低壓柔性轉(zhuǎn)子示意圖,轉(zhuǎn)子包含風(fēng)扇、多級(jí)增壓級(jí)及渦輪葉盤(pán),并采用0-2-1三支點(diǎn)支承方案。由于葉盤(pán)、錐殼等結(jié)構(gòu)使用,轉(zhuǎn)子具有較為復(fù)雜的質(zhì)量/剛度分布。同時(shí),此類轉(zhuǎn)子具有細(xì)長(zhǎng)軸段、懸臂大質(zhì)量/轉(zhuǎn)動(dòng)慣量、長(zhǎng)跨度支承等結(jié)構(gòu)特點(diǎn),為典型弱剛度柔性轉(zhuǎn)子,工作轉(zhuǎn)速范圍內(nèi)通常會(huì)有多階彎曲模態(tài)被激起。當(dāng)碰摩發(fā)生時(shí),碰摩產(chǎn)生的約束剛度往往與轉(zhuǎn)子自身的橫向等效剛度相當(dāng),這將導(dǎo)致柔性轉(zhuǎn)子多階模態(tài)均可能發(fā)生顯著變化。

    圖1 航空發(fā)動(dòng)機(jī)低壓柔性轉(zhuǎn)子結(jié)構(gòu)示意圖Fig.1 Schematic diagram of low-pressure flexible rotor structure in aero engines

    1.1 柔性轉(zhuǎn)子的有限元建模方法

    針對(duì)圖1所示的柔性轉(zhuǎn)子,建模的關(guān)鍵是準(zhǔn)確模擬其質(zhì)量/轉(zhuǎn)動(dòng)慣量、橫向剛度沿軸向的分布特征,為此采用模化能力較強(qiáng)的有限元法進(jìn)行動(dòng)力學(xué)建模。

    本文基于梁?jiǎn)卧⒋祟愞D(zhuǎn)子的模化方法,基本流程如下:首先基于低壓轉(zhuǎn)子的CAD模型進(jìn)行結(jié)構(gòu)簡(jiǎn)化,簡(jiǎn)化過(guò)程如圖2所示,忽略螺栓、篦齒等局部結(jié)構(gòu),采用等截面直梁模擬軸段和鼓筒,變截面梁模擬錐殼,彈簧單元模擬支承。對(duì)于葉盤(pán)結(jié)構(gòu),采用質(zhì)量單元模擬,保證質(zhì)量元的質(zhì)量和轉(zhuǎn)動(dòng)慣量與實(shí)際葉盤(pán)結(jié)構(gòu)的質(zhì)量和轉(zhuǎn)動(dòng)慣量一致;值得說(shuō)明的是,葉盤(pán)結(jié)構(gòu)也可采用梁?jiǎn)卧M,通過(guò)仔細(xì)劃分模擬葉盤(pán)結(jié)構(gòu)的梁?jiǎn)卧孛媸沟迷摬糠至簡(jiǎn)卧馁|(zhì)量和轉(zhuǎn)動(dòng)慣量與實(shí)際葉盤(pán)結(jié)構(gòu)一致即可。本文中,風(fēng)扇葉盤(pán)采用梁?jiǎn)卧M,其他葉盤(pán)結(jié)構(gòu)采用質(zhì)量元,最終建立的簡(jiǎn)化模型如圖3(a)所示。

    圖2 轉(zhuǎn)子模型簡(jiǎn)化示意圖Fig.2 Simplification process of rotor model

    其次,根據(jù)上述簡(jiǎn)化模型,在ANSYS中建立其有限元模型,梁?jiǎn)卧捎肂eam188單元,質(zhì)量單元為Mass21,彈簧單元為Combin14。梁?jiǎn)卧孛鎸傩院唾|(zhì)量單元實(shí)常數(shù)通過(guò)上述簡(jiǎn)化模型獲得,彈簧單元的實(shí)常數(shù)基于發(fā)動(dòng)機(jī)實(shí)際支承剛度獲得。所建立的有限元模型如圖3(b)所示,模型包含46個(gè)單元和40個(gè)節(jié)點(diǎn),需要說(shuō)明的是,梁?jiǎn)卧獙?shí)際為一維的中心線模型,圖3(b)是ANSYS中為了視覺(jué)效果而設(shè)置的截面顯示功能。

    為了保證模型的精度,一方面在進(jìn)行結(jié)構(gòu)簡(jiǎn)化和網(wǎng)格劃分時(shí),需要校核關(guān)鍵部件質(zhì)量、轉(zhuǎn)動(dòng)慣量與真實(shí)值的誤差,以及軸段的等效剛度與真實(shí)值的誤差,以此為基礎(chǔ)進(jìn)行適當(dāng)修正;另一方面,對(duì)于建立好的模型,也可以通過(guò)模態(tài)分析結(jié)果與試驗(yàn)?zāi)B(tài)結(jié)果或者精細(xì)的實(shí)體有限元模型模態(tài)分析結(jié)果進(jìn)行對(duì)比,驗(yàn)證模型的精度。本文所建立模型與所方提供的實(shí)體模型模態(tài)計(jì)算結(jié)果對(duì)比,所關(guān)心模態(tài)的誤差均在5%以下,滿足精度要求。

    圖3 低壓柔性轉(zhuǎn)子的梁?jiǎn)卧邢拊P虵ig.3 Beam element model of low-pressure flexible rotors

    1.2 轉(zhuǎn)靜件碰摩約束模型

    圖4所示為轉(zhuǎn)子系統(tǒng)碰摩約束模型,當(dāng)轉(zhuǎn)子振幅超過(guò)轉(zhuǎn)靜間隙r0時(shí),轉(zhuǎn)子與機(jī)匣發(fā)生碰摩,假設(shè)機(jī)匣剛度為kc,轉(zhuǎn)靜件之間的摩擦系數(shù)為μ,該截面p處轉(zhuǎn)子位移為xp和yp,令zp=xp+iyp,i為復(fù)數(shù),則碰摩對(duì)轉(zhuǎn)子產(chǎn)生的附加剛度為

    k′=H(|zp|-r0)kc(1+isign(vrel)μ)·

    (1-r0/|zp|)

    (1)

    式中:vrel=ω|zp|+Ω·rp,bd為碰摩點(diǎn)相對(duì)速度,ω為進(jìn)動(dòng)角速度,Ω為轉(zhuǎn)子旋轉(zhuǎn)角速度;rp,bd為截面p處的葉盤(pán)半徑;H(·)為Heaviside函數(shù);sign(·)為符號(hào)函數(shù),表達(dá)式為

    (2)

    (3)

    將碰摩剛度分解到x和y方向,可以獲得截面p處由碰摩產(chǎn)生的附加剛度矩陣:

    (4)

    可以看出,碰摩附加剛度矩陣為節(jié)點(diǎn)自由度xp和yp的非線性函數(shù),且其在非對(duì)角線元素上存在非0項(xiàng),即碰摩約束引起x和y向的交叉耦合。根據(jù)節(jié)點(diǎn)自由度編號(hào)將各截面處由碰摩產(chǎn)生的附加剛度矩陣進(jìn)行組集,便可得到整體的碰摩附加剛度矩陣K′。

    圖4 轉(zhuǎn)子系統(tǒng)的碰摩約束模型Fig.4 Rubbing constraint model of rotor system

    1.3 考慮碰摩約束的轉(zhuǎn)子動(dòng)力學(xué)方程

    根據(jù)1.1節(jié)柔性轉(zhuǎn)子動(dòng)力學(xué)建模方法,可在ANSYS獲得柔性轉(zhuǎn)子有限元模型,并將其質(zhì)量矩陣、剛度矩陣、阻尼矩陣及陀螺矩陣導(dǎo)出。進(jìn)一步在對(duì)應(yīng)自由度位置處引入1.2節(jié)的轉(zhuǎn)靜碰摩約束模型,可獲得含碰摩約束的柔性轉(zhuǎn)子動(dòng)力學(xué)方程:

    K′(u)u(t)=0

    (5)

    式中:M、C、G、K、u分別為質(zhì)量、阻尼、陀螺、剛度矩陣和位移向量;K′(u)為附加剛度矩陣,其為轉(zhuǎn)子自由度的非線性函數(shù)矩陣。

    2 基于諧波平衡的非線性模態(tài)求解方法

    2.1 諧波平衡法原理

    非線性模態(tài)由Rosenberg首次提出,將其定義為無(wú)阻尼非線性系統(tǒng)的一種同步周期振子[17],此后為了能夠考慮系統(tǒng)阻尼耗散,Krack[18]、Laxalde[19]等學(xué)者將線性系統(tǒng)復(fù)模態(tài)概念進(jìn)一步擴(kuò)展,提出了復(fù)非線性模態(tài)概念。根據(jù)復(fù)非線性模態(tài)概念,動(dòng)力學(xué)方程式(5)的特征值為

    (6)

    式中:ω0為特征頻率;ξ為模態(tài)阻尼比。此處假設(shè)模態(tài)阻尼比是與頻率無(wú)關(guān)的。

    根據(jù)非線性模態(tài)定義,非線性模態(tài)具有周期運(yùn)動(dòng)特性,故可以展開(kāi)成傅里葉級(jí)數(shù)形式,但由于系統(tǒng)存在阻尼,需引入慢變的衰減項(xiàng)e-kβt,其中,k=1,2,…,l表示諧波階次,l為諧波系數(shù)總數(shù),t為時(shí)間,則非線性模態(tài)運(yùn)動(dòng)表示為

    (7)

    該非線性模態(tài)運(yùn)動(dòng)的基頻即為λ=-β+iω,A0為靜力項(xiàng),Ak、Bk分別為正弦和余弦展開(kāi)項(xiàng)的系數(shù)向量。若忽略靜力項(xiàng)并僅保留k=1項(xiàng),上述解便退化為線性阻尼系統(tǒng)的指數(shù)解假設(shè);若系統(tǒng)為保守系統(tǒng),即不存在阻尼耗散,則慢變衰減項(xiàng)e-kβt=1。

    對(duì)于碰摩約束項(xiàng)g(u)=K′(u)u(t),由于K′(u)為位移u(t)的非線性函數(shù),故該項(xiàng)也可以展開(kāi)為傅里葉形式,如式(8)所示:

    Pksinkωt)

    (8)

    式中:P0為碰摩載荷的傅里葉展開(kāi)常數(shù)項(xiàng);Pk和Qk為碰摩載荷的傅里葉展開(kāi)正弦和余弦分量。

    將式(6)~式(8)代入到式(5)中,并分別令方程的e-kβtsinkωt和e-kβtcoskωt等于0,整理得到代數(shù)方程:

    (9)

    KA0+P0=0

    (10)

    其中:D=C+G。將各個(gè)階次的平衡方程進(jìn)行組裝,得到如下的代數(shù)方程組:

    (11)

    其中:

    Λk=

    H(Z,ω,β)=Λ(ω,β)Z+b(Z,ω)=0

    (12)

    其中:

    2.2 自由度縮減方法

    由于柔性轉(zhuǎn)子系統(tǒng)具有較多的自由度,因此式(12)所示的非線性代數(shù)方程通常具有較高規(guī)模,這將導(dǎo)致代數(shù)方程在求解過(guò)程中存在求解效率低、難于收斂問(wèn)題。為此,給出如下的自由度縮減方法。假設(shè)系統(tǒng)非線性自由度個(gè)數(shù)為m,線性自由度個(gè)數(shù)為s(m+s=N,N為轉(zhuǎn)子系統(tǒng)自由度總數(shù))。對(duì)于式(13)所示的非線性代數(shù)方程組,按線性和非線性自由度進(jìn)行分塊:

    (13)

    其中:

    式(13)進(jìn)行展開(kāi),可得:

    ΛqqZq+ΛqpZp+Rq=0

    (14)

    ΛpqZq+ΛppZp=0

    (15)

    根據(jù)式(15),得到Zp的表達(dá)式為

    (16)

    將式(16)代入(14),得到:

    (17)

    2.3 數(shù)值求解方法

    對(duì)前述非線性代數(shù)方程組進(jìn)行求解,還需解決如下兩方面問(wèn)題,其一是由于特征值λ=-β+iω的引入,使得非線性代數(shù)方程的待求變量比方程數(shù)量多2,為此,可借鑒線性系統(tǒng)模態(tài)歸一化的思想,引入附加方程,從而使代數(shù)方程組靜定。選取第m個(gè)自由度的一次諧波系數(shù)進(jìn)行歸一化處理,即令:

    (18)

    第2個(gè)問(wèn)題是如何根據(jù)位移u(t)的傅里葉系數(shù)獲得碰摩項(xiàng)g(u)=K′(u)u(t)的傅里葉系數(shù),此處通過(guò)時(shí)頻轉(zhuǎn)換技術(shù)予以計(jì)算,如式(19)所示,通過(guò)位移的各項(xiàng)諧波系數(shù),可以確定位移的時(shí)域表達(dá),根據(jù)時(shí)域中碰摩約束項(xiàng)與位移間的顯示關(guān)系,確定時(shí)域中的碰摩約束項(xiàng)表達(dá),便可以得到碰摩約束項(xiàng)的各項(xiàng)諧波系數(shù)。變換時(shí),傅里葉系數(shù)與其時(shí)域間的關(guān)系通過(guò)離散傅里葉變換與逆變換實(shí)現(xiàn)。

    (19)

    最終,獲得如下的非線性代數(shù)方程:

    (20)

    X(k+1)=X(k)-F′(X(k))-1F(X(k))

    (21)

    計(jì)算時(shí),取線性模態(tài)分析結(jié)果作為初值X(0),F(xiàn)′(X(k))為向量函數(shù)的Jacobi矩陣,通過(guò)有限差分法計(jì)算得到。

    2.4 計(jì)算流程

    結(jié)合第1節(jié)建模及本節(jié)的計(jì)算方法,本文提出求解航空發(fā)動(dòng)機(jī)轉(zhuǎn)子系統(tǒng)非線性模態(tài)的一般流程,如圖5所示。

    圖5 航空發(fā)動(dòng)機(jī)轉(zhuǎn)子非線性模態(tài)分析流程Fig.5 Modal analysis process for rotor system in aero-engines

    首先結(jié)合CAD簡(jiǎn)化模型獲得轉(zhuǎn)子截面參數(shù)、質(zhì)量轉(zhuǎn)動(dòng)慣量參數(shù)等參數(shù)后,基于ANSYS軟件建立其梁?jiǎn)卧P?,而后?dǎo)出HBMAT格式的動(dòng)力學(xué)矩陣(質(zhì)量矩陣、剛度矩陣等)以及Mapping映射文件(含自由度編號(hào)信息);在MATLAB文件中,首先根據(jù)HBMAT格式的動(dòng)力學(xué)矩陣和Mapping文件獲得滿秩格式的動(dòng)力學(xué)矩陣,根據(jù)節(jié)點(diǎn)自由度編號(hào)引入碰摩約束項(xiàng),建立考慮碰摩約束的轉(zhuǎn)子動(dòng)力學(xué)方程,進(jìn)一步結(jié)合本文建立的轉(zhuǎn)子非線性模態(tài)計(jì)算方法,獲得碰摩約束轉(zhuǎn)子的模態(tài)特性。

    3 計(jì)算結(jié)果與討論

    3.1 無(wú)碰摩時(shí)轉(zhuǎn)子的模態(tài)特性

    針對(duì)圖3的柔性轉(zhuǎn)子系統(tǒng),本節(jié)分析其在碰摩約束下的非線性模態(tài)特性。首先了解該轉(zhuǎn)子系統(tǒng)在無(wú)碰摩狀態(tài)下的固有特性,根據(jù)第1節(jié)建立的轉(zhuǎn)子模型,基于ANSYS軟件對(duì)不同轉(zhuǎn)速下轉(zhuǎn)子進(jìn)行線性特征值計(jì)算,得到各轉(zhuǎn)速下轉(zhuǎn)子線性模態(tài)頻率,圖6和圖7所示分別為相應(yīng)的坎貝爾圖及兩階正進(jìn)動(dòng)臨界轉(zhuǎn)速下的振型,結(jié)果表明:該轉(zhuǎn)子在工作轉(zhuǎn)速范圍內(nèi)存在兩階正/反進(jìn)動(dòng)臨界轉(zhuǎn)速,一階正/反進(jìn)動(dòng)臨界轉(zhuǎn)速為2 720 r/min和1 480 r/min;二階正/反進(jìn)動(dòng)臨界轉(zhuǎn)速為4 660 r/min 和2 210 r/min。對(duì)應(yīng)振型分別為轉(zhuǎn)子一階彎曲和渦輪平動(dòng)振型。

    圖6 無(wú)碰摩時(shí)柔性轉(zhuǎn)子的坎貝爾圖Fig.6 Campbell diagram of flexible rotor without rub-impact

    圖7 兩階臨界轉(zhuǎn)速下轉(zhuǎn)子振型Fig.7 Mode shapes at two order critical speeds

    3.2 碰摩影響下轉(zhuǎn)子模態(tài)特性

    考慮風(fēng)扇碰摩對(duì)轉(zhuǎn)子模態(tài)特性的影響,給定如下計(jì)算參數(shù):機(jī)匣剛度kc=107N/m,摩擦系數(shù)μ=0.2,轉(zhuǎn)靜間隙r0=0.2 mm,葉盤(pán)直徑800 mm,轉(zhuǎn)子1#支點(diǎn)存在阻尼,值為10 000 N·s/m,同時(shí)假設(shè)轉(zhuǎn)子工作于巡航轉(zhuǎn)速3 600 r/min。由于非線性模態(tài)具有能量相關(guān)性,需考慮不同振動(dòng)狀態(tài)下轉(zhuǎn)子模態(tài),本文通過(guò)風(fēng)扇幅值大小來(lái)反映轉(zhuǎn)子振動(dòng)水平,在2.3節(jié)的模態(tài)歸一化中給定不同的風(fēng)扇幅值,而后基于2.4節(jié)計(jì)算流程獲得不同風(fēng)扇幅值下轉(zhuǎn)子模態(tài)頻率,結(jié)果如圖8所示。本文在分析中主要考慮工作轉(zhuǎn)速范圍內(nèi)的兩階模態(tài),圖中結(jié)果可得出如下結(jié)論:

    1) 轉(zhuǎn)子振動(dòng)水平較小時(shí),風(fēng)扇幅值小于轉(zhuǎn)靜間隙,未發(fā)生碰摩,通過(guò)本文非線性模態(tài)求解方法計(jì)算得到的模態(tài)頻率與圖6中ANSYS線性特征值計(jì)算結(jié)果進(jìn)行對(duì)比,2種方法計(jì)算結(jié)果相同,側(cè)面驗(yàn)證了本文建立的非線性模態(tài)計(jì)算方法的有效性。

    2) 當(dāng)風(fēng)扇幅值超過(guò)轉(zhuǎn)靜間隙時(shí),風(fēng)扇機(jī)匣約束使得轉(zhuǎn)子各階正/反進(jìn)動(dòng)模態(tài)頻率增加,增加趨勢(shì)呈現(xiàn)逐漸變緩趨勢(shì),當(dāng)風(fēng)扇振幅達(dá)到一定程度,即轉(zhuǎn)子振動(dòng)能量足夠高時(shí),轉(zhuǎn)子各階模特頻率趨近于定值??傮w而言,考慮碰摩約束影響時(shí)轉(zhuǎn)子模特頻率具有區(qū)間分布特征,不同振動(dòng)狀態(tài)下轉(zhuǎn)子模態(tài)頻率均位于該區(qū)間內(nèi)。

    圖8 風(fēng)扇碰摩對(duì)兩階模態(tài)頻率的影響(巡航轉(zhuǎn)速3 600 r/min)Fig.8 Influence of fan rubbing on two order modal frequencies (Operating speed 3 600 r/min)

    表1 不同風(fēng)扇相對(duì)幅值下轉(zhuǎn)子模態(tài)頻率的變化率

    以上分析了碰摩約束轉(zhuǎn)子在特定轉(zhuǎn)速下的模態(tài),針對(duì)不同轉(zhuǎn)速,分別計(jì)算轉(zhuǎn)子模態(tài)頻率隨振幅變化曲線,提取各轉(zhuǎn)速下轉(zhuǎn)子模態(tài)頻率區(qū)間最小值和最大值,如圖9所示,進(jìn)而可獲得碰摩約束轉(zhuǎn)子的各階臨界轉(zhuǎn)速。結(jié)果表明:風(fēng)扇碰摩下,轉(zhuǎn)子一階正/反進(jìn)動(dòng)臨界轉(zhuǎn)速區(qū)間分別為[2 720,3 460]r/min和[1 480,1 810]r/min;二階正/反進(jìn)動(dòng)臨界轉(zhuǎn)速區(qū)間分別為[4 660,4 683]r/min和[2 210,2 260]r/min,風(fēng)扇振動(dòng)為主的一階臨界轉(zhuǎn)速受風(fēng)扇碰摩約束改變較為顯著,使得該階臨界轉(zhuǎn)速相對(duì)巡航轉(zhuǎn)速裕度顯著降低。航空發(fā)動(dòng)機(jī)對(duì)氣動(dòng)效率的高要求,葉尖間隙通常接近零間隙設(shè)計(jì),工作時(shí)碰摩時(shí)常發(fā)生,臨界轉(zhuǎn)速分析結(jié)果說(shuō)明實(shí)際工程中考慮碰摩約束帶來(lái)的轉(zhuǎn)子臨界轉(zhuǎn)速及其裕度是極其必要的。

    圖9 風(fēng)扇碰摩約束下轉(zhuǎn)子系統(tǒng)坎貝爾圖Fig.9 Campbell diagram of rotor system with fan rub-impact

    圖10所示為碰摩約束對(duì)轉(zhuǎn)子振型的影響,圖中碰摩約束下的轉(zhuǎn)子振型均取模態(tài)頻率接近極限值(風(fēng)扇振幅2 mm)時(shí)的結(jié)果,同時(shí)為便于對(duì)比,將無(wú)碰摩時(shí)和碰摩時(shí)的振型均作振幅歸一化,兩階振型均取遠(yuǎn)離風(fēng)扇碰摩位置的渦輪作為基準(zhǔn)進(jìn)行歸一化。結(jié)果表明:碰摩對(duì)轉(zhuǎn)子振型產(chǎn)生一定影響,對(duì)于風(fēng)扇振動(dòng)為主的一階模態(tài),風(fēng)扇碰摩約束可使得風(fēng)扇位置的相對(duì)振動(dòng)幅值有所降低,即抑制了風(fēng)扇振動(dòng);而對(duì)于渦輪振動(dòng)為主的二階模態(tài),風(fēng)扇碰摩約束使得風(fēng)扇相對(duì)振動(dòng)略有增加,二階振型存在由渦輪振動(dòng)向低壓轉(zhuǎn)子二階彎曲演化的趨勢(shì),但總體而言風(fēng)扇碰摩約束對(duì)二階振型的影響相對(duì)較小。

    圖10 風(fēng)扇碰摩約束對(duì)轉(zhuǎn)子振型的影響Fig.10 Influence of fan rub-impact on rotor modal shapes

    根據(jù)2.1節(jié),轉(zhuǎn)子模態(tài)阻尼對(duì)應(yīng)著特征值實(shí)部,其決定著模態(tài)運(yùn)動(dòng)隨時(shí)間的衰減速率。對(duì)于大多數(shù)線性系統(tǒng),模態(tài)阻尼通常為正,系統(tǒng)模態(tài)運(yùn)動(dòng)隨時(shí)間衰減,而對(duì)于非線性系統(tǒng),非線性力的引入可能導(dǎo)致模態(tài)阻尼為負(fù),因而其模態(tài)并不總是穩(wěn)定的。圖11所示為轉(zhuǎn)子各階模態(tài)阻尼比隨轉(zhuǎn)子振動(dòng)幅值的變化曲線,結(jié)果表明:未碰摩時(shí),轉(zhuǎn)子各階模態(tài)阻尼比為正,其原因是由于1#支點(diǎn)存在阻尼;當(dāng)風(fēng)扇振幅超過(guò)轉(zhuǎn)靜間隙時(shí),轉(zhuǎn)靜件發(fā)生碰摩,各階模態(tài)阻尼比將發(fā)生變化,對(duì)于正進(jìn)動(dòng)模態(tài),其模態(tài)阻尼比隨振幅增加,而反進(jìn)動(dòng)模態(tài)阻尼則隨振幅減小,當(dāng)碰摩足夠劇烈時(shí),轉(zhuǎn)子系統(tǒng)的

    反進(jìn)動(dòng)模態(tài)阻尼比為負(fù),此時(shí)反進(jìn)動(dòng)模態(tài)不穩(wěn)定。正反進(jìn)動(dòng)模態(tài)阻尼的變化機(jī)理主要與碰摩接觸點(diǎn)處摩擦力做功有關(guān)[14],對(duì)于反進(jìn)動(dòng)模態(tài),摩擦力對(duì)轉(zhuǎn)子系統(tǒng)的模態(tài)運(yùn)動(dòng)做正功,增加轉(zhuǎn)子能量,導(dǎo)致反進(jìn)動(dòng)模態(tài)阻尼比下降,當(dāng)摩擦力輸入能量超過(guò)轉(zhuǎn)子支點(diǎn)阻尼耗散能力時(shí),轉(zhuǎn)子等效模態(tài)阻尼比為負(fù),此時(shí)模態(tài)失穩(wěn)。需要說(shuō)明的是,由反進(jìn)動(dòng)模態(tài)不穩(wěn)定誘發(fā)的轉(zhuǎn)子振動(dòng)響應(yīng)失穩(wěn)現(xiàn)象已在簡(jiǎn)單轉(zhuǎn)子試驗(yàn)器上得到了證實(shí)[3]。然而,在實(shí)際航空發(fā)動(dòng)機(jī)轉(zhuǎn)子中,反進(jìn)動(dòng)模態(tài)不穩(wěn)定導(dǎo)致的轉(zhuǎn)子振動(dòng)失穩(wěn)卻鮮有出現(xiàn),主要原因可能在于航空發(fā)動(dòng)機(jī)轉(zhuǎn)子中通常具有擠壓油膜阻尼器,同時(shí)轉(zhuǎn)子中大量連接界面的存在等,使得轉(zhuǎn)子系統(tǒng)具有足夠的相對(duì)阻尼,因此轉(zhuǎn)子系統(tǒng)的振動(dòng)能量能夠得到有效的耗散。

    3.3 碰摩參數(shù)的影響規(guī)律

    實(shí)際工程中,航空發(fā)動(dòng)機(jī)風(fēng)扇機(jī)匣可能采用不同材質(zhì)并內(nèi)涂不同材質(zhì)/硬度的涂層[20],這對(duì)機(jī)匣自身剛度以及葉片/機(jī)匣之間的摩擦系數(shù)有明顯影響。因此,本節(jié)重點(diǎn)分析機(jī)匣剛度和摩擦系數(shù)對(duì)碰摩轉(zhuǎn)子模態(tài)特性的影響。

    3.3.1 機(jī)匣剛度

    機(jī)匣剛度分別取5×106N/m、1×107N/m和2×107N/m,其他參數(shù)與3.1節(jié)相同。計(jì)算得到不同機(jī)匣剛度下,轉(zhuǎn)子兩階正/反進(jìn)動(dòng)模態(tài)頻率隨振幅變化曲線,如圖12所示??梢钥闯?,隨機(jī)匣剛度增加,轉(zhuǎn)子兩階正/反進(jìn)動(dòng)模態(tài)頻率均增加,其原因是在相同振動(dòng)狀態(tài)下,機(jī)匣剛度越高,其對(duì)

    圖12 機(jī)匣剛度對(duì)模態(tài)頻率的影響Fig.12 Influence of casing stiffness on modal frequency

    轉(zhuǎn)子系統(tǒng)的約束作用越強(qiáng),導(dǎo)致模態(tài)頻率增加的越高。另外,從機(jī)匣剛度對(duì)兩階模態(tài)頻率影響的相對(duì)變化而言,如表2所示,一階模態(tài)頻率對(duì)碰摩剛度變化更為敏感,而二階模態(tài)頻率對(duì)碰摩剛度變化則極不敏感,以兩階正進(jìn)動(dòng)模態(tài)頻率為例,當(dāng)碰摩剛度由5×106N/m增加至2×107N/m,一階正進(jìn)動(dòng)

    表2 不同碰摩剛度下柔性轉(zhuǎn)子系統(tǒng)模態(tài)頻率變化率(風(fēng)扇幅值2 mm)

    模態(tài)頻率變化率可由8.94%增加至31.49%;而二階正進(jìn)動(dòng)模態(tài)頻率變化率由0.06%增加至0.19%。

    由于反進(jìn)動(dòng)模態(tài)存在失穩(wěn)的可能,因此分析機(jī)匣剛度對(duì)反進(jìn)動(dòng)模態(tài)阻尼的影響,如圖13所示,可以看出:隨機(jī)匣剛度增加,反進(jìn)動(dòng)模態(tài)阻尼比等于0的臨界振幅逐漸減小,表明機(jī)匣剛度越高,碰摩時(shí)反進(jìn)動(dòng)模態(tài)越容易失穩(wěn)。

    圖14所示為不同機(jī)匣剛度下低壓轉(zhuǎn)子的振型,不同機(jī)匣剛度下轉(zhuǎn)子振型均相對(duì)渦輪進(jìn)行了振幅歸一化處理。結(jié)果表明:① 對(duì)于風(fēng)扇振動(dòng)為主的一階彎曲振型,隨著機(jī)匣剛度增加,風(fēng)扇位置的相對(duì)振幅逐漸增加,表明機(jī)匣剛度越高,對(duì)風(fēng)扇振動(dòng)的抑制作用愈發(fā)明顯;② 對(duì)于渦輪振動(dòng)為主的二階振型,隨著機(jī)匣剛度增加,風(fēng)扇相對(duì)振幅有增加的趨勢(shì),低壓轉(zhuǎn)子向二階彎曲振型演化的趨勢(shì)愈發(fā)顯著。

    圖14 機(jī)匣剛度對(duì)模態(tài)振型的影響Fig.14 Influence of casing stiffness on modal shapes

    3.3.2 摩擦系數(shù)

    本節(jié)分析摩擦系數(shù)對(duì)轉(zhuǎn)子模態(tài)特性的影響,如前所述,風(fēng)扇碰摩對(duì)轉(zhuǎn)子一階模態(tài)頻率的影響顯著高于二階模態(tài),故此處僅一階反進(jìn)動(dòng)模態(tài)為例,給出計(jì)算結(jié)果。圖15所示為摩擦系數(shù)取0.1、 0.2和0.4時(shí),一階反進(jìn)動(dòng)模態(tài)頻率和阻尼變化曲線。結(jié)果表明,摩擦系數(shù)對(duì)轉(zhuǎn)子模態(tài)頻率的影響極小,表明轉(zhuǎn)子模態(tài)頻率主要取決于接觸點(diǎn)處法向約束,而對(duì)切向的摩擦作用極不敏感;而對(duì)于反進(jìn)動(dòng)模態(tài)阻尼,隨摩擦系數(shù)增加,轉(zhuǎn)子反進(jìn)動(dòng)模態(tài)阻尼比顯著下降,同時(shí)模態(tài)阻尼比等于0的臨界振幅逐漸減小,表明摩擦系數(shù)越高,反進(jìn)動(dòng)模態(tài)越容易失穩(wěn)。

    圖15 摩擦系數(shù)對(duì)轉(zhuǎn)子一階反進(jìn)動(dòng)模態(tài)的影響Fig.15 Influence of friction coefficient on first order backward whirl mode

    4 結(jié) 論

    1) 本文針對(duì)航空發(fā)動(dòng)機(jī)中具有復(fù)雜結(jié)構(gòu)和質(zhì)量/剛度分布的柔性轉(zhuǎn)子系統(tǒng),結(jié)合梁?jiǎn)卧ㄌ岢隽丝煽紤]碰摩約束的復(fù)雜轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)建模方法;將諧波平衡思想與頻域的自由度縮減技術(shù)結(jié)合,建立了適用于復(fù)雜非線性轉(zhuǎn)子系統(tǒng)非線性模態(tài)特性的求解方法;基于ANSYS和MATLAB平臺(tái),將建模方法與求解方法融合,提出了可適用于工程中含大規(guī)模自由度的非線性轉(zhuǎn)子系統(tǒng)模態(tài)特性的計(jì)算分析流程。建模及求解方法的適用性在某航空發(fā)動(dòng)機(jī)柔性轉(zhuǎn)子系統(tǒng)的非線性模態(tài)分析中得到了驗(yàn)證。

    2) 碰摩約束使轉(zhuǎn)子系統(tǒng)模態(tài)頻率/臨界轉(zhuǎn)速顯著增加,且轉(zhuǎn)子振幅越高,機(jī)匣對(duì)轉(zhuǎn)子約束作用越強(qiáng),模態(tài)頻率也越高,針對(duì)文中所分析的發(fā)動(dòng)機(jī)柔性轉(zhuǎn)子對(duì)象及參數(shù),風(fēng)扇碰摩可使一階彎曲型的正/反進(jìn)動(dòng)模態(tài)頻率最高增加約16%和29%。但需要注意的是,碰摩約束轉(zhuǎn)子的模態(tài)頻率始終位于特定區(qū)間范圍內(nèi)。此外,風(fēng)扇碰摩同時(shí)對(duì)轉(zhuǎn)子一階彎曲型的模態(tài)振型產(chǎn)生一定抑制作用。陀螺效應(yīng)、轉(zhuǎn)子振型以及機(jī)匣剛度對(duì)碰摩約束帶來(lái)的轉(zhuǎn)子模態(tài)頻率變化有較為顯著的影響,但摩擦系數(shù)的影響可以忽略。

    3) 由于碰摩點(diǎn)處切向摩擦力做功影響,碰摩對(duì)轉(zhuǎn)子系統(tǒng)模態(tài)阻尼亦會(huì)產(chǎn)生顯著影響,轉(zhuǎn)子各階正進(jìn)動(dòng)模態(tài)阻尼始終大于0,而反進(jìn)動(dòng)模態(tài)阻尼在碰摩較為劇烈時(shí)可能小于0,表明碰摩轉(zhuǎn)子的正進(jìn)動(dòng)模態(tài)始終穩(wěn)定,而反進(jìn)動(dòng)模態(tài)存在失穩(wěn)的可能。結(jié)合模態(tài)阻尼變化機(jī)理及關(guān)鍵參數(shù)影響規(guī)律,實(shí)際中通過(guò)擠壓油膜阻尼器提升轉(zhuǎn)子系統(tǒng)的阻尼耗散能力、抑制碰摩程度等方式均可降低反進(jìn)動(dòng)模態(tài)失穩(wěn)的可能性。

    與線性模態(tài)類似,非線性模態(tài)為非線性轉(zhuǎn)子系統(tǒng)的一種固有特性,必然在理解轉(zhuǎn)子非線性振動(dòng)響應(yīng)方面能發(fā)揮重要作用。本文揭示了碰摩約束對(duì)航空發(fā)動(dòng)機(jī)低壓轉(zhuǎn)子非線性模態(tài)的影響,進(jìn)一步,碰摩轉(zhuǎn)子的非線性模態(tài)與其振動(dòng)響應(yīng)之間的影響關(guān)系如何,例如非線性模態(tài)在何種運(yùn)動(dòng)狀態(tài)下被激起,不同階次的非線性模態(tài)對(duì)碰摩轉(zhuǎn)子動(dòng)力學(xué)行為的演化有何貢獻(xiàn),這些理論問(wèn)題將在未來(lái)予以重點(diǎn)關(guān)注。

    猜你喜歡
    模態(tài)系統(tǒng)
    Smartflower POP 一體式光伏系統(tǒng)
    WJ-700無(wú)人機(jī)系統(tǒng)
    ZC系列無(wú)人機(jī)遙感系統(tǒng)
    基于PowerPC+FPGA顯示系統(tǒng)
    半沸制皂系統(tǒng)(下)
    連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
    車輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對(duì)比
    國(guó)內(nèi)多模態(tài)教學(xué)研究回顧與展望
    高速顫振模型設(shè)計(jì)中顫振主要模態(tài)的判斷
    基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
    久久热在线av| 视频区欧美日本亚洲| 看黄色毛片网站| 99riav亚洲国产免费| 午夜福利高清视频| 一本综合久久免费| 一个人观看的视频www高清免费观看 | 精品久久久久久久人妻蜜臀av | 制服丝袜大香蕉在线| 精品国产国语对白av| 操出白浆在线播放| 欧美日韩中文字幕国产精品一区二区三区 | 中亚洲国语对白在线视频| av中文乱码字幕在线| 成熟少妇高潮喷水视频| 亚洲中文av在线| 午夜精品国产一区二区电影| 亚洲精品av麻豆狂野| 国产精品美女特级片免费视频播放器 | 国产精品久久久久久人妻精品电影| 18禁观看日本| 超碰成人久久| 精品欧美一区二区三区在线| 亚洲七黄色美女视频| 日本免费一区二区三区高清不卡 | 久久精品国产综合久久久| 美女 人体艺术 gogo| 久久久久久国产a免费观看| 久久精品91蜜桃| 制服诱惑二区| 91成年电影在线观看| 午夜a级毛片| 曰老女人黄片| 国产精品一区二区在线不卡| 国产麻豆成人av免费视频| 亚洲欧美一区二区三区黑人| 精品久久久久久久人妻蜜臀av | 欧美激情极品国产一区二区三区| 97人妻精品一区二区三区麻豆 | 亚洲精品一区av在线观看| 亚洲国产日韩欧美精品在线观看 | 免费在线观看视频国产中文字幕亚洲| www.www免费av| 国产精品久久久久久精品电影 | 女警被强在线播放| 免费av毛片视频| 亚洲 欧美 日韩 在线 免费| 亚洲在线自拍视频| 可以在线观看的亚洲视频| 国产av精品麻豆| av中文乱码字幕在线| 女人爽到高潮嗷嗷叫在线视频| 91国产中文字幕| 国产亚洲欧美在线一区二区| 午夜福利欧美成人| 亚洲欧美日韩另类电影网站| 色播在线永久视频| 亚洲 欧美一区二区三区| 欧美另类亚洲清纯唯美| 久久久水蜜桃国产精品网| 日韩欧美一区视频在线观看| 国产精品电影一区二区三区| 午夜免费激情av| 欧美日韩中文字幕国产精品一区二区三区 | 中文字幕av电影在线播放| 日本撒尿小便嘘嘘汇集6| 久久精品91无色码中文字幕| 色综合站精品国产| 久久狼人影院| 18禁观看日本| 免费在线观看影片大全网站| 精品高清国产在线一区| 亚洲成国产人片在线观看| 香蕉国产在线看| 久久精品成人免费网站| 国产三级黄色录像| 国产免费av片在线观看野外av| 97人妻精品一区二区三区麻豆 | 一本精品99久久精品77| 欧美一级a爱片免费观看看| 国产久久久一区二区三区| 欧美色欧美亚洲另类二区| 制服丝袜大香蕉在线| 日本 av在线| 最好的美女福利视频网| 悠悠久久av| 赤兔流量卡办理| 午夜福利成人在线免费观看| 老司机福利观看| 观看美女的网站| 亚洲无线在线观看| 成年女人永久免费观看视频| 色综合亚洲欧美另类图片| 一进一出抽搐gif免费好疼| 日日夜夜操网爽| 99热这里只有精品一区| 国产淫片久久久久久久久| 99riav亚洲国产免费| 两个人视频免费观看高清| 国产高清不卡午夜福利| 极品教师在线视频| 毛片女人毛片| 69人妻影院| 亚洲精品色激情综合| 国产精品久久久久久精品电影| 精品人妻1区二区| 伦理电影大哥的女人| 少妇熟女aⅴ在线视频| 免费在线观看日本一区| 国产精品伦人一区二区| 99精品久久久久人妻精品| 女人十人毛片免费观看3o分钟| 久久人人爽人人爽人人片va| 国产精品av视频在线免费观看| 成人午夜高清在线视频| 91av网一区二区| 成人国产一区最新在线观看| 国内精品美女久久久久久| 看免费成人av毛片| 老司机福利观看| 性欧美人与动物交配| 一本一本综合久久| 日本一二三区视频观看| 有码 亚洲区| 99在线人妻在线中文字幕| 久久99热这里只有精品18| 国产av在哪里看| 不卡一级毛片| 久久精品91蜜桃| 国产伦一二天堂av在线观看| 国产精品人妻久久久久久| 亚洲在线自拍视频| 联通29元200g的流量卡| 精品久久久久久久人妻蜜臀av| 亚洲欧美日韩高清在线视频| 天堂影院成人在线观看| 韩国av在线不卡| 欧美3d第一页| 如何舔出高潮| 日本在线视频免费播放| 老司机午夜福利在线观看视频| 成年女人毛片免费观看观看9| bbb黄色大片| 午夜爱爱视频在线播放| 中国美女看黄片| 嫁个100分男人电影在线观看| 男插女下体视频免费在线播放| 性插视频无遮挡在线免费观看| 美女黄网站色视频| 亚州av有码| 日韩欧美精品v在线| 久久精品夜夜夜夜夜久久蜜豆| 又粗又爽又猛毛片免费看| 亚洲精品亚洲一区二区| 久久精品久久久久久噜噜老黄 | 成人一区二区视频在线观看| 午夜免费男女啪啪视频观看 | 一夜夜www| 99九九线精品视频在线观看视频| 亚洲av一区综合| a级毛片免费高清观看在线播放| 欧美日韩乱码在线| 男女之事视频高清在线观看| 亚洲精品久久国产高清桃花| 国内精品宾馆在线| 久久午夜福利片| 欧美黑人巨大hd| 男人和女人高潮做爰伦理| 国产精品嫩草影院av在线观看 | 两人在一起打扑克的视频| 乱码一卡2卡4卡精品| 日韩一区二区视频免费看| 老司机深夜福利视频在线观看| 亚洲一区二区三区色噜噜| 国产三级在线视频| 国产精品免费一区二区三区在线| 国产黄色小视频在线观看| 桃色一区二区三区在线观看| 免费av观看视频| av福利片在线观看| 日本五十路高清| 我的女老师完整版在线观看| 听说在线观看完整版免费高清| 亚洲精华国产精华精| 99在线视频只有这里精品首页| 国产视频内射| 亚洲在线观看片| 日韩精品中文字幕看吧| 我的老师免费观看完整版| 色噜噜av男人的天堂激情| 俄罗斯特黄特色一大片| x7x7x7水蜜桃| 国产精品人妻久久久久久| 日韩一本色道免费dvd| 国产精品国产三级国产av玫瑰| 色播亚洲综合网| 在现免费观看毛片| 色av中文字幕| 亚洲精品乱码久久久v下载方式| 3wmmmm亚洲av在线观看| 男女下面进入的视频免费午夜| 十八禁网站免费在线| 校园人妻丝袜中文字幕| 成人av在线播放网站| 亚洲欧美精品综合久久99| 国产精品久久视频播放| 桃红色精品国产亚洲av| 亚洲成人久久爱视频| 人妻丰满熟妇av一区二区三区| 亚洲国产日韩欧美精品在线观看| 亚洲一区二区三区色噜噜| 亚洲av中文av极速乱 | 一边摸一边抽搐一进一小说| 午夜福利成人在线免费观看| 免费av观看视频| 亚洲色图av天堂| 99国产极品粉嫩在线观看| 久久亚洲精品不卡| 成人特级av手机在线观看| 性色avwww在线观看| 两人在一起打扑克的视频| 亚洲三级黄色毛片| 99视频精品全部免费 在线| 无遮挡黄片免费观看| 日韩欧美三级三区| 日本免费a在线| 一本久久中文字幕| 精品无人区乱码1区二区| 中文在线观看免费www的网站| 国产伦一二天堂av在线观看| 久久精品91蜜桃| 一级a爱片免费观看的视频| 国产伦精品一区二区三区视频9| 久久久久久久午夜电影| 99国产精品一区二区蜜桃av| 三级毛片av免费| 亚洲欧美精品综合久久99| 国产aⅴ精品一区二区三区波| 久久精品国产亚洲av涩爱 | 五月伊人婷婷丁香| 欧美一区二区国产精品久久精品| 99国产精品一区二区蜜桃av| 午夜精品在线福利| 女人被狂操c到高潮| 可以在线观看毛片的网站| 国产视频一区二区在线看| 久久久午夜欧美精品| 午夜福利欧美成人| 床上黄色一级片| 嫩草影院新地址| 久久久国产成人免费| 亚洲avbb在线观看| 特级一级黄色大片| 老熟妇仑乱视频hdxx| 日韩欧美三级三区| 国产精品自产拍在线观看55亚洲| 性色avwww在线观看| 国产69精品久久久久777片| 日本在线视频免费播放| 免费av观看视频| 欧美三级亚洲精品| 国产精品av视频在线免费观看| 国产精品三级大全| 国产精品久久视频播放| 性欧美人与动物交配| 亚洲在线自拍视频| 国产高清不卡午夜福利| 精华霜和精华液先用哪个| 精品乱码久久久久久99久播| 69av精品久久久久久| 老司机午夜福利在线观看视频| 国产高清不卡午夜福利| 欧美日韩中文字幕国产精品一区二区三区| 99在线视频只有这里精品首页| 国产色婷婷99| 日韩中文字幕欧美一区二区| 欧美日韩精品成人综合77777| 狠狠狠狠99中文字幕| 夜夜夜夜夜久久久久| 一个人看视频在线观看www免费| 99国产极品粉嫩在线观看| 长腿黑丝高跟| 亚洲精品久久国产高清桃花| 91麻豆av在线| 成人永久免费在线观看视频| 亚洲真实伦在线观看| 国产亚洲91精品色在线| 麻豆国产av国片精品| 人妻丰满熟妇av一区二区三区| av.在线天堂| 国产伦精品一区二区三区视频9| 99riav亚洲国产免费| 成人一区二区视频在线观看| 1000部很黄的大片| АⅤ资源中文在线天堂| 欧美日韩乱码在线| 夜夜夜夜夜久久久久| netflix在线观看网站| 国产精品久久久久久久电影| 婷婷色综合大香蕉| 免费观看人在逋| 欧美色欧美亚洲另类二区| 亚洲av成人精品一区久久| 欧美日韩瑟瑟在线播放| 国产 一区精品| 琪琪午夜伦伦电影理论片6080| 日本欧美国产在线视频| 高清在线国产一区| 亚洲国产精品合色在线| 日本 av在线| 亚洲欧美日韩卡通动漫| 男人狂女人下面高潮的视频| 黄色一级大片看看| 最新中文字幕久久久久| 日韩欧美精品免费久久| 国产大屁股一区二区在线视频| 变态另类成人亚洲欧美熟女| 99久久精品热视频| 1024手机看黄色片| 久久久久久久亚洲中文字幕| 日本黄大片高清| 国产男人的电影天堂91| 久久午夜亚洲精品久久| 欧美日韩亚洲国产一区二区在线观看| 禁无遮挡网站| 人妻久久中文字幕网| 亚洲性久久影院| 国内少妇人妻偷人精品xxx网站| 91麻豆精品激情在线观看国产| 国产精品嫩草影院av在线观看 | 熟女人妻精品中文字幕| 亚洲av成人精品一区久久| 一进一出抽搐gif免费好疼| 99在线人妻在线中文字幕| 亚洲国产色片| 国产精品电影一区二区三区| 婷婷精品国产亚洲av在线| 九九热线精品视视频播放| 亚洲人与动物交配视频| 国产激情偷乱视频一区二区| 亚洲中文字幕一区二区三区有码在线看| 国产成年人精品一区二区| 国产精华一区二区三区| 乱码一卡2卡4卡精品| 人妻夜夜爽99麻豆av| 长腿黑丝高跟| 久久精品国产亚洲网站| 赤兔流量卡办理| 此物有八面人人有两片| 欧美性猛交黑人性爽| 亚洲成a人片在线一区二区| 人妻制服诱惑在线中文字幕| 久久国产乱子免费精品| 天美传媒精品一区二区| 国产精品伦人一区二区| 看免费成人av毛片| 亚洲av免费在线观看| 99久久精品国产国产毛片| 大又大粗又爽又黄少妇毛片口| 91久久精品电影网| 综合色av麻豆| 国产单亲对白刺激| 老熟妇乱子伦视频在线观看| 国产免费av片在线观看野外av| 日本 欧美在线| 午夜精品一区二区三区免费看| 久久精品91蜜桃| 床上黄色一级片| 狂野欧美白嫩少妇大欣赏| 18+在线观看网站| 亚洲精品成人久久久久久| 99久久无色码亚洲精品果冻| 国产精品永久免费网站| 在线天堂最新版资源| 亚洲一级一片aⅴ在线观看| 欧美xxxx性猛交bbbb| 两性午夜刺激爽爽歪歪视频在线观看| 成人国产麻豆网| 18禁裸乳无遮挡免费网站照片| 亚洲美女搞黄在线观看 | 窝窝影院91人妻| 美女免费视频网站| 欧美性猛交╳xxx乱大交人| 国产日本99.免费观看| 亚洲无线观看免费| 欧美日韩精品成人综合77777| videossex国产| 天堂网av新在线| 精品久久久久久,| 亚洲黑人精品在线| 亚洲美女搞黄在线观看 | 天堂av国产一区二区熟女人妻| 22中文网久久字幕| 成人午夜高清在线视频| 别揉我奶头 嗯啊视频| 三级国产精品欧美在线观看| 狂野欧美激情性xxxx在线观看| 最新在线观看一区二区三区| 日韩,欧美,国产一区二区三区 | 国产精品福利在线免费观看| 九九热线精品视视频播放| 别揉我奶头~嗯~啊~动态视频| 日日摸夜夜添夜夜添av毛片 | 日本欧美国产在线视频| 免费人成在线观看视频色| 亚洲国产精品sss在线观看| 一区二区三区四区激情视频 | 欧美一区二区精品小视频在线| or卡值多少钱| 亚洲欧美精品综合久久99| 午夜福利在线观看免费完整高清在 | 亚洲内射少妇av| 三级毛片av免费| 国产精品99久久久久久久久| 日韩精品有码人妻一区| 久久午夜亚洲精品久久| 免费一级毛片在线播放高清视频| 午夜福利高清视频| 国产精品免费一区二区三区在线| 毛片女人毛片| 亚洲午夜理论影院| 男人舔女人下体高潮全视频| 亚洲国产色片| 亚洲国产精品合色在线| 欧美黑人欧美精品刺激| 午夜福利成人在线免费观看| 成年女人看的毛片在线观看| 国产精品精品国产色婷婷| 国产精品久久久久久久久免| 免费观看在线日韩| 色综合站精品国产| 精品午夜福利视频在线观看一区| 黄色日韩在线| 国产精品精品国产色婷婷| 国产精品永久免费网站| 亚洲av中文av极速乱 | 色视频www国产| 国产精品不卡视频一区二区| 国产主播在线观看一区二区| 国产精品嫩草影院av在线观看 | 动漫黄色视频在线观看| 亚洲最大成人av| 女的被弄到高潮叫床怎么办 | 亚洲最大成人av| 乱码一卡2卡4卡精品| 男人舔奶头视频| videossex国产| 中文在线观看免费www的网站| 国内精品美女久久久久久| 一区二区三区四区激情视频 | 婷婷色综合大香蕉| 久9热在线精品视频| 91狼人影院| 男女做爰动态图高潮gif福利片| 色5月婷婷丁香| 国内精品久久久久精免费| 美女免费视频网站| 乱系列少妇在线播放| 成人国产综合亚洲| 色哟哟·www| 黄色丝袜av网址大全| 久久久午夜欧美精品| 午夜老司机福利剧场| 成人永久免费在线观看视频| 国产精品人妻久久久影院| 人妻久久中文字幕网| 国产精品一区二区性色av| 嫩草影院新地址| 国产探花在线观看一区二区| 婷婷色综合大香蕉| 一进一出抽搐动态| 九九在线视频观看精品| 男女之事视频高清在线观看| 精品久久久久久久人妻蜜臀av| 免费av毛片视频| 蜜桃亚洲精品一区二区三区| 在线观看av片永久免费下载| 性色avwww在线观看| 天天一区二区日本电影三级| a级毛片a级免费在线| 一本精品99久久精品77| 婷婷精品国产亚洲av在线| 亚洲欧美日韩高清专用| .国产精品久久| 18禁在线播放成人免费| 亚洲三级黄色毛片| 亚洲美女视频黄频| 国产成人av教育| 99视频精品全部免费 在线| 亚洲欧美日韩卡通动漫| 中文字幕高清在线视频| 免费av毛片视频| 久久久国产成人精品二区| 久久精品91蜜桃| 亚洲精品影视一区二区三区av| 免费一级毛片在线播放高清视频| 免费av毛片视频| 麻豆精品久久久久久蜜桃| 亚洲电影在线观看av| 亚洲aⅴ乱码一区二区在线播放| 伦精品一区二区三区| 舔av片在线| 可以在线观看毛片的网站| 国产精品人妻久久久影院| 欧美另类亚洲清纯唯美| 欧美日本亚洲视频在线播放| 国产人妻一区二区三区在| 国产午夜福利久久久久久| 性插视频无遮挡在线免费观看| 搞女人的毛片| 丰满乱子伦码专区| 99热这里只有是精品50| 少妇丰满av| 男人狂女人下面高潮的视频| 在线观看av片永久免费下载| 亚洲av免费高清在线观看| 一个人看的www免费观看视频| 欧美日韩乱码在线| 国产亚洲精品综合一区在线观看| 国产男人的电影天堂91| 午夜福利成人在线免费观看| 精品久久久久久,| 欧美3d第一页| 亚洲真实伦在线观看| 韩国av在线不卡| 欧美一区二区精品小视频在线| 免费观看的影片在线观看| 国产视频一区二区在线看| 亚洲av熟女| 久久久久免费精品人妻一区二区| 变态另类丝袜制服| 日本一本二区三区精品| 国产成人a区在线观看| 国产三级中文精品| 99精品久久久久人妻精品| 亚洲一区二区三区色噜噜| 成人无遮挡网站| 搞女人的毛片| 淫秽高清视频在线观看| 欧美潮喷喷水| 99在线人妻在线中文字幕| 99热精品在线国产| 国产国拍精品亚洲av在线观看| 亚洲男人的天堂狠狠| 男女下面进入的视频免费午夜| 人妻久久中文字幕网| 色噜噜av男人的天堂激情| 特大巨黑吊av在线直播| 亚洲精品在线观看二区| 美女 人体艺术 gogo| 22中文网久久字幕| 人妻少妇偷人精品九色| 男女啪啪激烈高潮av片| 51国产日韩欧美| 久久精品人妻少妇| 狂野欧美白嫩少妇大欣赏| 此物有八面人人有两片| 国内精品久久久久久久电影| 日韩精品青青久久久久久| 能在线免费观看的黄片| 国产高潮美女av| 精品一区二区三区视频在线| 亚洲第一区二区三区不卡| 18禁裸乳无遮挡免费网站照片| av天堂中文字幕网| 精品久久久久久久末码| 老熟妇乱子伦视频在线观看| 欧美在线一区亚洲| 精品福利观看| 中亚洲国语对白在线视频| 亚洲专区国产一区二区| 欧美成人a在线观看| 精品福利观看| 欧美成人a在线观看| 国产精品1区2区在线观看.| 精品久久久久久久末码| 欧美日韩综合久久久久久 | 在线观看午夜福利视频| 中文亚洲av片在线观看爽| 国模一区二区三区四区视频| 999久久久精品免费观看国产| 久久久色成人| 色哟哟哟哟哟哟| 1024手机看黄色片| 三级国产精品欧美在线观看| 精品欧美国产一区二区三| 精品人妻偷拍中文字幕| 很黄的视频免费| 免费看a级黄色片| 干丝袜人妻中文字幕| 在线观看一区二区三区| 日韩欧美国产在线观看| 欧美黑人欧美精品刺激| 日日夜夜操网爽| 久久久久久久午夜电影| 亚洲国产日韩欧美精品在线观看| 超碰av人人做人人爽久久| 少妇熟女aⅴ在线视频| 夜夜爽天天搞| 国产视频一区二区在线看| 美女cb高潮喷水在线观看| 成人av一区二区三区在线看| 无人区码免费观看不卡| 国产午夜精品久久久久久一区二区三区 | 日本色播在线视频| 桃色一区二区三区在线观看| 美女大奶头视频| 一卡2卡三卡四卡精品乱码亚洲| 丰满的人妻完整版| av在线老鸭窝| 国产一区二区亚洲精品在线观看| 尾随美女入室| 亚洲七黄色美女视频| .国产精品久久| 男人和女人高潮做爰伦理|