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

旋轉(zhuǎn)葉片-機匣碰摩振動響應分析

2017-07-19 12:37:07太興宇郝玉明聞邦椿
振動與沖擊 2017年14期
關鍵詞:有限元振動

馬 輝, 孫 祺, 太興宇, 郝玉明, 聞邦椿

(1. 東北大學 機械工程與自動化學院,沈陽 110819;2. 機械結(jié)構(gòu)強度與振動國家重點實驗室,西安 710049; (3. 沈陽鼓風機集團股份有限公司,沈陽 110869)

旋轉(zhuǎn)葉片-機匣碰摩振動響應分析

馬 輝1,2, 孫 祺1, 太興宇3, 郝玉明3, 聞邦椿1

(1. 東北大學 機械工程與自動化學院,沈陽 110819;2. 機械結(jié)構(gòu)強度與振動國家重點實驗室,西安 710049; (3. 沈陽鼓風機集團股份有限公司,沈陽 110869)

以航空發(fā)動機壓氣機葉片為研究對象,基于ANSYS有限元軟件,考慮離心剛化、旋轉(zhuǎn)軟化以及科氏力對葉片的影響,采用變厚度殼單元建立了旋轉(zhuǎn)葉片的有限元模型。假定圓盤與機匣存在靜態(tài)不對中,葉片由于離心力的作用產(chǎn)生徑向伸長,從而使葉尖與機匣內(nèi)壁發(fā)生碰摩。根據(jù)葉尖處每個節(jié)點徑向位移與間隙之間的位置關系,對不同轉(zhuǎn)速下葉片-機匣碰摩過程進行數(shù)值仿真。研究結(jié)果表明:對于扭形葉片,葉尖只有部分節(jié)點與機匣發(fā)生碰摩;碰摩結(jié)束后,葉片振動響應以低階動頻形式進行衰減,連續(xù)的周期碰摩產(chǎn)生的倍頻成分在與之接近的動頻處會出現(xiàn)幅值放大現(xiàn)象。研究結(jié)果可為葉片結(jié)構(gòu)設計及葉片-機匣碰摩故障診斷提供理論依據(jù)。

變厚度殼單元;葉片-機匣碰摩;高頻振動;動頻;幅值放大

為提高發(fā)動機的性能和效率,就要使轉(zhuǎn)子葉尖與機匣的間隙盡可能小,以減少因工作介質(zhì)泄漏而造成的損失。然而間隙過小,會導致高速旋轉(zhuǎn)的葉片在離心載荷的作用下與機匣發(fā)生碰撞。由于葉尖處有著較大的線速度,碰撞能量大,所以碰摩一旦發(fā)生,就會造成嚴重的后果。據(jù)美國運輸部報道:在1962—1976年間的417百萬飛行小時中,10.2%的發(fā)動機轉(zhuǎn)子事故是由于轉(zhuǎn)靜子之間碰摩引起的[1]。

Padovan等[2]將葉片簡化為懸臂梁,推導了法向接觸力和葉片徑向變形之間的關系,分析了在單葉片和多葉片碰摩情況下,不平衡量、葉片/轉(zhuǎn)子剛度、系統(tǒng)阻尼和摩擦對系統(tǒng)非線性動力學特性的影響。Sinha[3]研究了帶葉片的柔性轉(zhuǎn)子與剛性機匣的碰撞動力學特性,碰摩通過附加的剛度矩陣和阻尼矩陣來模擬,研究發(fā)現(xiàn)當葉尖與機匣發(fā)生接觸時,接觸力類似于Hertzian接觸產(chǎn)生的力。Lesaffre等[4-5]基于Sinha提出的轉(zhuǎn)子葉片模型,并將葉片進一步簡化成具有兩自由度的集中質(zhì)量塊模型,機匣簡化成柔性環(huán),研究了葉片轉(zhuǎn)子轉(zhuǎn)速大于臨界轉(zhuǎn)速時引起機匣-葉片系統(tǒng)的不穩(wěn)定現(xiàn)象。太興宇等[6]建立旋轉(zhuǎn)葉片的連續(xù)體解析模型,采用間隙函數(shù)判定,確定碰摩是否發(fā)生,分析了科氏力對碰摩響應的耦合作用。

根據(jù)試驗測試葉片-機匣接觸力結(jié)果,發(fā)現(xiàn)對于葉片單點以及局部碰摩,其接觸力類似于周期性脈沖力[7-9],根據(jù)這一特定碰摩情況,一些學者提出了基于脈沖力模型的碰摩故障模擬。Sinha[10]提出基于旋轉(zhuǎn)的Timoshenko懸臂梁模型理論,應用數(shù)值方法分析了不同的摩擦參數(shù)和接觸時間對于葉片頻率的影響以及葉片的非線性動力學響應特性。Turner等[11-12]將碰摩力簡化為脈沖力,分析了葉片機匣連續(xù)撞擊情況下的瞬態(tài)動力學特性,確定了葉尖接觸力分布和葉片在不同碰摩情況下的振動形式。太興宇等[13]采用Timoshenko梁單元建立了旋轉(zhuǎn)葉片在脈沖力作用下的動力學方程,并分析了不同激振頻率下葉片的振動響應。

為了更加精確地描述葉片-機匣之間的碰摩特征,一些學者采用接觸動力學理論來模擬葉尖與機匣的接觸特性。Legrand等[14]采用數(shù)值計算以及軟件分析輔助的方法研究了內(nèi)側(cè)帶有阻尼涂層的機匣與葉尖接觸的特性。Roques等[15]以某核電廠汽輪發(fā)電機組的轉(zhuǎn)子系統(tǒng)為研究對象,建立了相應的轉(zhuǎn)子-靜子有限元模型,基于接觸動力學理論分析了轉(zhuǎn)軸與靜子隔板間的碰摩現(xiàn)象。劉書國等[16]對航空發(fā)動機葉片-機匣碰摩過程進行了數(shù)值模擬,考慮實際葉片的葉形特征,分析葉片頂部受瞬時碰撞與摩擦載荷(碰摩載荷)作用下的動力學影響特征。

通過對上述文獻的分析發(fā)現(xiàn),解析模型通常采用較為簡單的懸臂梁和懸臂板模型,并且大多采用脈沖加載的方式來模擬碰摩過程,并不能真實的反映實際的碰摩情況;而基于接觸動力學的方法,可以較為真實的模擬碰摩過程,但是需要耗費大量的計算時間,效率較低。考慮離心剛化、旋轉(zhuǎn)軟化和科氏力的影響,本文采用變厚度殼單元對真實葉片進行有限元建模,建立笛卡爾坐標系下的間隙函數(shù),對葉尖各個節(jié)點進行間隙判定,確定碰摩發(fā)生的位置,并討論了不同轉(zhuǎn)速下葉片的碰摩響應。

1 旋轉(zhuǎn)葉片振動的有限元動力學方程

旋轉(zhuǎn)葉片振動的動力學方程可表示為

(1)

令C1=C+G,K1=K+S-Kspin,旋轉(zhuǎn)葉片的動力學方程可轉(zhuǎn)化為

(2)

2 葉片-機匣碰摩動力學仿真

2.1 變厚度殼葉片的有限元模型

由于真實葉片為變截面扭形葉片,本文基于ANSYS有限元軟件,通過抽取葉片中面數(shù)據(jù),根據(jù)葉片型線數(shù)據(jù),將葉片的厚度ε擬合為隨空間位置變化的變量,進而將葉片簡化為變厚度殼單元。本文采用Shell181單元來模擬變厚度扭形葉片,此單元具有四個節(jié)點,每個節(jié)點有六個自由度,分別是沿X、Y、Z方向的平動自由度和繞X、Y、Z軸的轉(zhuǎn)動自由度。葉尖沿寬度方向等分成20份,21個位置。圖1為變厚度殼葉片有限元模型,給出了葉尖每個位置的節(jié)點號。變厚度殼葉片的有限元模型適合分析薄及中等厚度的板殼結(jié)構(gòu)零件,并支持線性、大扭轉(zhuǎn)和大應變,變厚度非線性等分析。將ε擬合為如下表達式[17]

ε=a0+a1x+a2y+a3z+a4xy+a5xz+a6yz+a7x2+a8y2+a9z2+a10x2y+a11x2z+a12y2x+a13y2z+a14z2y+a15z2x+a16x3+a17y3+a18z3

(3)

2.2 葉片的固有特性

葉片動頻是指葉片繞盤軸旋轉(zhuǎn)時的自振頻率,它隨轉(zhuǎn)速的變化而變化,即是轉(zhuǎn)速的函數(shù)。圖2為葉片前5階動頻隨轉(zhuǎn)速的變化曲線。由圖可知隨著轉(zhuǎn)速的增加各階動頻隨之增大。

圖1 變厚度殼葉片有限元模型Fig. 1 Finite element model of blade

選用的葉片參數(shù)如下:葉片彈性模量E=1.25×1011Pa,密度ρ=4 370 kg/m3,泊松比υ=0.3,圓盤半徑Rd=216.52 mm。忽略葉根-圓盤榫連結(jié)構(gòu)的影響,這里葉根的約束形式為完全固支。

圖2 葉片動頻Fig. 2 Dynamic frequencies of the blade

2.3 間隙函數(shù)

假設機匣與圓盤存在靜態(tài)平行不對中,并且葉片具有安裝角,則葉尖與機匣內(nèi)壁之間的間隙呈不均勻分布。圖3為不對中葉片-機匣間隙示意圖。圖3(b)中O和O1分別為圓盤中心和機匣中心;L為葉片長度;Rc為機匣半徑;rg為葉尖軌跡半徑,rg=L+Rd;Rd為圓盤半徑。

圖3 葉片-機匣間隙示意圖Fig. 3 Schematic diagram of the clearance between blade and casing

ANSYS在模擬葉片-機匣的碰摩過程中,在葉片上施加的是轉(zhuǎn)速效應,葉片并不轉(zhuǎn)動,因此在笛卡爾坐標系下,通過建立時變的間隙函數(shù)來模擬葉尖徑向和機匣之間的間隙。由于葉片存在安裝角,葉尖各節(jié)點在XOZ平面上的投影的位置不同,導致在葉尖各節(jié)點處存在不同的初始相位φi,從而出現(xiàn)不同的初始間隙ci。在圖3(b)中l(wèi)O1B=Rc,lOiO=Rc-cmin-rg,lOA=rg,lAB=ci,∠O1OA=π-Ωt,根據(jù)正弦定理建立笛卡爾坐標系下的葉片-機匣的初始間隙函數(shù)如下

ci=Rcsin(π/2-θi(t))/sin (π/2+φi(t))-rg

(4)

(5)

為了更加真實的反映碰摩力在葉尖上的分布情況,葉尖處的每個節(jié)點都要判定是否侵入,從而確定是否發(fā)生碰摩。葉尖節(jié)點i的侵入量δi的表達式如下

(6)

葉尖各個節(jié)點的碰摩力值需要根據(jù)對應節(jié)點位置的侵入量來分配,其表達式為

(7)

(8)

fn表達式如下

(9)

式中:Γ1=Γ0/kc,

Γ0=3EI/L3+ρAΩ2(81L/280+3Rd/8),

α=(Rd+L)/L

其中L為葉片長度;EI為葉片抗彎剛度;A為葉片橫截面積;Ω為轉(zhuǎn)速。由于航空發(fā)動機機匣的厚度一般為1~2 mm,所以這里設機匣的當量剛度kc=5×106N/m,摩擦因數(shù)μ=0.3。根據(jù)以上分析,得出模擬碰摩過程的流程圖,如圖4所示。

圖4 葉片-機匣碰摩模擬過程流程圖Fig. 4 Flowchart of blade-casing rubbing simulation process

2.4 穩(wěn)態(tài)響應分析

在碰摩模擬過程中,先使葉片在離心力作用下運動,葉片振動響應穩(wěn)定后,再施加碰摩載荷。本文對葉片在5 000 r/min,10 000 r/min和15 000 r/min三種轉(zhuǎn)速下的穩(wěn)態(tài)響應進行分析。

圖5為不同轉(zhuǎn)速下葉片節(jié)點1、節(jié)點13和節(jié)點

圖5 不同轉(zhuǎn)速下葉尖碰摩位移響應Fig. 5 Tangential rubbing-induced displacement responses under different rotational speeds

21的碰摩位移響應。從圖中可以看到,轉(zhuǎn)速為5 000 r/min時,葉尖各個節(jié)點的Z向位移均沒有超過最小間隙值,此時只有葉片在離心力作用下產(chǎn)生的位移;當轉(zhuǎn)速達到10 000 r/min時,節(jié)點21的Z向位移超過了最小間隙,于是在節(jié)點21處發(fā)生碰摩,而節(jié)點1和節(jié)點13的Z向位移都沒有超過最小間隙值,但是依然在碰摩時產(chǎn)生相應的位移。表明扭形葉片在節(jié)點21處發(fā)生碰摩時,會連帶葉尖其他節(jié)點處產(chǎn)生位移。而當轉(zhuǎn)速達到15 000 r/min時,節(jié)點13的Z向位移也超過了最小間隙,故節(jié)點13處除了有碰摩位移也會出現(xiàn)其他節(jié)點的連帶位移。圖中,節(jié)點1處實際上并沒有發(fā)生碰摩,并且Z向位移隨著轉(zhuǎn)速的增加而減小,說明扭形葉片自身的形狀對節(jié)點1是否發(fā)生碰摩有很大的影響。

由圖5可知,在5 000 r/min時葉尖各節(jié)點都不與機匣發(fā)生碰摩,所以下面只給出10 000 r/min和15 000 r/min下不同時刻各節(jié)點的Z向碰摩力,如圖6所示。從圖中可以更加清楚的分析葉尖各節(jié)點處發(fā)生碰摩的情況。在圖6(a)中,葉尖處只有沿Y軸負方向的5個節(jié)點處(節(jié)點17、18、19、20和21)發(fā)生碰摩,碰摩接觸時間越長,碰摩力越大。而圖6(b)中發(fā)生碰摩的節(jié)點數(shù)要多于圖6(a),同樣為沿葉片Y軸負方向的部分節(jié)點發(fā)生碰摩。結(jié)果表明,變厚度扭形葉片在離心力作用下,葉尖與機匣會發(fā)生局部碰摩。

圖6 不同轉(zhuǎn)速下不同時刻各節(jié)點的Z向碰摩力Fig. 6 Forces of each node at different moments in Z direction under different rotational speeds

由圖6可知,節(jié)點21處的碰摩最嚴重,因此這里分析節(jié)點21處位移響應特征。圖7為10 000 r/min和15 000 r/min轉(zhuǎn)速時節(jié)點21處的X向位移響應。單周期碰摩響應指僅碰摩1次后的自由振動響應,其頻譜圖的峰值對應的頻率為葉片的動頻。由圖7可得以下特征:單周期碰摩時節(jié)點21的X向位移幅值隨著時間逐漸衰減;在多周期碰摩發(fā)生之后,葉尖節(jié)點21逐漸穩(wěn)定為周期運動;10 000 r/min下葉片碰摩的振動形式主要由1階振型(1階彎曲)、2階振型(1階扭轉(zhuǎn))、3階振型(2階彎曲)組成,位移響應的頻率特征表現(xiàn)為在1階動頻處 (5×)、2階動頻處 (12×)和3階動頻處 (16×)出現(xiàn)幅值放大現(xiàn)象;15 000 r/min下碰摩發(fā)生時葉尖的振動形式主要由1階振型(1階彎曲)、2階振型(1階扭轉(zhuǎn))、3階振型(2階彎曲)和4階振型(彎曲與扭轉(zhuǎn)耦合)組成,其頻率特征為在1階動頻處 (3×)、2階動頻處 (8×) 、3階動頻處 (11×)和4階動頻處(17×)出現(xiàn)幅值放大現(xiàn)象;10 000 r/min和15 000 r/min的單周期碰摩響應頻譜圖的峰值所對應的葉片動頻與ANSYS中求出的葉片動頻比較吻合。

3 結(jié) 論

考慮了葉片的旋轉(zhuǎn)效應的影響,采用變厚度殼單元建立了葉片的有限元模型,建立了笛卡爾坐標系下的間隙函數(shù),對葉尖各個節(jié)點進行間隙判定,從而確定碰摩發(fā)生的位置。通過對不同轉(zhuǎn)速下的穩(wěn)態(tài)碰摩響應分析,得到以下結(jié)論:

(1) 對于真實的葉片,由于其復雜的曲面特征,葉尖上一側(cè)點的徑向位移隨著轉(zhuǎn)速的增加而增加,另一側(cè)由于彎扭運動,徑向位移減小。當葉尖一側(cè)節(jié)點的徑向位移超過最小間隙時,只會在部分節(jié)點處發(fā)生局部碰摩,另一側(cè)節(jié)點處則產(chǎn)生連帶運動。

圖7 10 000 r/min和15 000 r/min時節(jié)點21的X向位移響應Fig. 7 Displacement response of node 21 in X direction at 10 000 r/min and 15 000 r/min

(2) 葉片在碰摩載荷的作用下葉片振動在葉片彎曲、扭轉(zhuǎn)和彎扭耦合動頻處會出現(xiàn)較大的幅值,而葉片的彎曲振動更容易被激發(fā)。并且在高轉(zhuǎn)速情況下葉片的振動形式會更為復雜,更容易激發(fā)更高階的動頻響應。

[1] IKEDA T, ISHIDA Y, YAMAMOTO T, et al. Nonlinear forced oscillations of an unsymmetrical shaft and an unsymmetrical rotor with quartic nonlinearity [J]. Bulletin of the JSME, 1988, 31(3): 530-538.

[2] PADOVAN J, CHOY F K. Nonlinear dynamics of rotor/blade/casing rub interactions[J]. Journal of Turbomachinery, 1987, 109: 527-534.

[3] SINHA S K. Dynamic characteristics of a flexible bladed-rotor with Coulomb damping due to tip-rub[J]. Journal of Sound and Vibration, 2004, 273: 875-919.

[4] LESAFFRE N, SINOU J J, THOUVEREZ F. Contact analysis of a flexible bladed-rotor[J]. European Journal of Mechanics A/Solids, 2007, 26: 541-557.

[5] LESAFFRE N, SINOU J J, THOUVEREZ F. Stability analysis of rotating beams rubbing on an elastic circular structure[J]. Journal of Sound and Vibration, 2007, 299(1): 1005-1032.

[6] 太興宇, 馬輝, 譚禎, 等. 基于連續(xù)體旋轉(zhuǎn)梁模型的碰摩故障動力學特征分析[J]. 振動與沖擊, 2013, 32(18): 43-48. TAI Xingyu, MA Hui, TAN Zhen, et al. Dynamic characteristics of a continuous rotating beam model with a rubbing fault[J]. Journal of Vibration and Shock, 2013, 32(18): 43-48.

[7] CHEN G. Study on the recognition of aero-engine blade-casing rubbing fault based on the casing vibration acceleration [J]. Measurement, 2015, 65: 71-80.

[8] CHEN G. Simulation of casing vibration resulting from blade-casing rubbing and its verifications [J]. Journal of Sound and Vibration, 2016, 361: 190-209.

[9] MA H, TAI X, HAN Q, et al. A revised model for rubbing between rotating blade and elastic casing [J]. Journal of Sound and Vibration, 2015, 337: 301-320.

[10] SINHA S K. Non-linear dynamic response of a rotating radial Timoshenko beam with periodic pulse loading at the free-end[J]. International Journal of Nonlinear Mechanics, 2005, 40: 113-149.

[11] TURNER K, ADAMS M, DUNN M. Simulation of engine blade tip-rub induced vibration[C]∥Proceedings of GT2005, Ren-Tahoe. Nevada, USA, 2005.

[12] TURNER K, DUNN M, PADOVA C. Airfoil deflection characteristics during rub events[J]. Journal of Turbomachinery, 2012, 134: 011018-1-7.

[13] 太興宇, 馬輝, 譚禎, 等. 脈沖力加載下的葉片-機匣動力學特性研究[J]. 東北大學學報(自然科學版), 2012, 33(12): 1759-1799. TAI Xingyu, MA Hui, TAN Zhen, et al. Research on dynamic characteristics of blade-casing with impulse loading [J]. Journal of Northeastern University (Natural Science), 2012, 33(12): 1759-1799.

[14] LEGRAND M, BATAILLY A, PIERRE C. Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law[J]. Journal of Computational and Nonlinear Dynamics, 2012, 7: 011010-1-11.

[15] ROQUES S, LEGRAND M, CARTRAUD P, et al. Modeling of a rotor speed transient response with radial rubbing [J]. Journal of Sound and Vibration, 2010, 329: 527-546.

[16] 劉書國,洪杰,陳萌. 航空發(fā)動機葉片-機匣碰摩過程的數(shù)值模擬[J]. 航空動力學報, 2011, 26(6): 1282-1288. LIU Shuguo, HONG Jie, CHEN Meng, et al. Numerical simulation of the dynamic process of aero-engine blade-to-case rub-impact [J]. Journal of Aerospace Power, 2011, 26(6): 1282-1288.

[17] 馬輝,太興宇,張志,等. 根部聯(lián)接剛度對旋轉(zhuǎn)葉片振動特性的影響[J]. 東北大學學報(自然科學版),2013, 34(2): 265-270. MA Hui, TAI Xingyu, ZHANG Zhi, et al. Influence of root connecting stiffness on vibration responses of rotating blades [J]. Journal of Northeastern University (Natural Science), 2013, 34(2): 265-270.

Vibration response analysis on the rotating blade-casing rubbing

MA Hui1,2,SUN Qi1,TAI Xingyu3,HAO Yuming3,WEN Bangchun1

(1.School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819,China; 2.State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University, Xi’an 710049, China;3.Shengu Group Compressor Co., Ltd., Shenyang 110869, China)

Taking a compressor blade as the research object and considering the influences of centrifugal stiffening, spin softening and Coriolis force, a finite element model of the rotating blade was established by using the variable thickness shell element based on ANSYS software. Assuming that there is a static misalignment between the disk and casing, the blade will generate a radial elongation under the impact of centrifugal force, moreover, the blade-casing rubbing will appear. The numerical simulation on rubbing processes under different rotational speeds was carried out according to the location relation between the radial displacements and the clearances at each top node of the blade. The results show that the rubbing only occurs on partial nodes for twisting blades, the vibration response of the blade will present a decay with a low-order dynamic frequency, the amplitude amplification phenomenon will appear when the continual periodic rubbing caused multiple frequency components are close to corresponding dynamic frequencies during rubbing. The simulation results provide a theoretical foundation to the blade structure design and the fault diagnosis on blade-casing rubbing.

variable thickness shell element; blade-casing rubbing; high-frequency vibration; dynamic frequency; amplitude amplification

國家自然科學基金委員會與中國民用航空局聯(lián)合資助項目(U1433109);中央高校基本科研業(yè)務費專項資金(N140301001; N150305001)

2016-01-19 修改稿收到日期: 2016-06-17

馬輝 男,博士,教授,1978年生

TH212;TH213.3

A

10.13465/j.cnki.jvs.2017.14.004

猜你喜歡
有限元振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
This “Singing Highway”plays music
新型有機玻璃在站臺門的應用及有限元分析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
磨削淬硬殘余應力的有限元分析
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: A级毛片无码久久精品免费| 亚洲三级视频在线观看| 尤物成AV人片在线观看| 亚洲国产欧洲精品路线久久| 东京热av无码电影一区二区| igao国产精品| 亚洲毛片网站| 重口调教一区二区视频| 婷婷午夜天| 国产福利一区视频| 国产精品xxx| 97综合久久| 91青青草视频在线观看的| 亚洲综合第一区| 久久精品国产一区二区小说| 成人福利免费在线观看| 国产自在线播放| 麻豆精品在线播放| 国产清纯在线一区二区WWW| 国产自视频| 亚洲三级a| 在线精品亚洲一区二区古装| 日韩无码黄色网站| 久久精品午夜视频| 91在线国内在线播放老师| 欧美激情综合| 国产精品高清国产三级囯产AV| 国产尤物在线播放| 国产噜噜在线视频观看| 极品国产在线| 成人福利在线视频| 婷婷六月激情综合一区| jijzzizz老师出水喷水喷出| 国产亚洲精品97AA片在线播放| 亚洲国产成人麻豆精品| 国产在线一二三区| 制服丝袜一区| 凹凸精品免费精品视频| 成人国产一区二区三区| 青青国产成人免费精品视频| 国产黄在线免费观看| 国产情侣一区二区三区| 亚洲二三区| 男女性午夜福利网站| 欧美国产综合色视频| 一区二区日韩国产精久久| 国产成人一级| 久久成人国产精品免费软件| 91在线日韩在线播放| 国产精品美女网站| 一级香蕉人体视频| 自拍偷拍欧美日韩| 无码中文字幕乱码免费2| 伊人查蕉在线观看国产精品| 亚洲无码久久久久| 免费一级毛片在线观看| 国产69囗曝护士吞精在线视频| 亚洲二区视频| 欧美日韩一区二区三区在线视频| 久久精品波多野结衣| 毛片国产精品完整版| 成人国产精品网站在线看| 国产成年无码AⅤ片在线| 色偷偷一区二区三区| 日韩毛片免费观看| 91综合色区亚洲熟妇p| 亚洲色图欧美在线| 国产精品视频免费网站| 国产成人精品高清在线| 波多野结衣一区二区三区四区视频 | 四虎成人在线视频| 亚洲精品成人片在线观看| 99久久精品国产自免费| 精品欧美日韩国产日漫一区不卡| 欧美精品高清| 国产精品无码制服丝袜| 在线a网站| 一本大道香蕉久中文在线播放| 欧美一区中文字幕| 日韩精品无码不卡无码| 99资源在线| 国产91高跟丝袜|