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

機動飛行下磁流變阻尼器-轉子系統動力學特性

2024-06-15 06:31:08王俊劉云飛秦朝燁馬梁洪芳芳褚福磊
振動工程學報 2024年5期

王俊 劉云飛 秦朝燁 馬梁 洪芳芳 褚福磊

摘要: 航空發動機在機動飛行過程中,工作條件非常惡劣,飛行過程會產生不規則的瞬態振動,易引發故障。采用有限元法建立機動飛行下基于雙線性本構方程的磁流變阻尼器?轉子系統有限元模型,并利用Newmark?β數值方法進行求解,研究轉子系統在機動飛行過程中的動態特性。在此基礎上,考慮磁流變阻尼器作用,研究其對沖擊載荷下轉子系統瞬態及穩態響應的影響。結果表明,機動飛行開始和結束瞬間會產生瞬態沖擊,激發轉子系統一階模態響應。在合適的電流作用下,變阻尼器可以有效抑制機動飛行過程中轉子系統瞬態及穩態響應。此外,在機動飛行下,由于軸頸離心率較大,易導致磁流變阻尼器產生非線性行為。

關鍵詞: 航空發動機; 機動飛行; 磁流變阻尼器; 轉子系統; 瞬態沖擊

中圖分類號: V231.96??? 文獻標志碼: A??? 文章編號: 1004-4523(2024)05-0747-09

DOI:10.16385/j.cnki.issn.1004-4523.2024.05.003

引? 言

機動飛行包括瞬時轉彎、爬坡?下降、加速、持續轉彎和翻滾等,是飛機特技飛行和空戰中的常見動作,對航空發動機轉子系統的動力學特性有著重要影響[1?2]。

近年來,國內外學者對機動飛行環境下轉子系統的動力學特性進行了大量的研究。El?Saeidy等[3]考慮基礎激勵和不平衡質量,建立了剛性轉子系統的動力學方程,利用解析方法得到了系統的時域分析結果。祝長生等[4]利用Lagrange方程建立了飛機在任意機動飛行條件下的轉子系統動力學統一模型,討論了飛機的典型機動飛行動作對轉子系統動力學特性的影響。Andrés等[5]實驗研究了機動載荷下氣體軸承?轉子系統動力學響應。Yu等[6]研究了不同機動載荷對轉子?聯軸器系統的動力學特性的影響,發現增大機動速度會增強聯軸器的非線性行為,使轉子系統的振動由單周期運動轉變為多周期、分岔或混沌運動。Gao等[7?9]對機動飛行下轉子系統動力學特性進行了較為詳細的研究,通過理論及實驗分析了機動載荷對軸承、擠壓油膜阻尼器支承下正常及故障轉子系統非線性行為的影響。Chen等[10]研究了擠壓油膜阻尼器?轉子系統動力學特性,結果表明對于較大的機動飛行動作,需要增加油膜的間隙,以降低油膜的非線性水平。

機動飛行時,發動機轉子系統將承受很大的附加離心力和陀螺力矩,并且這些載荷隨飛行狀態和時間發生變化,導致系統運行不穩定。特別在機動載荷作用下,發動機轉子發生較大的撓曲變形,擠壓油膜阻尼器(Squeeze Film Damper, SFD)會因為油膜過度被擠壓而非線性顯著增強,產生雙穩態跳躍、鎖死和非協調進動等嚴重的有害現象,進而惡化轉子系統的振動狀態,甚至可能造成轉子失穩。且SFD作為被動式支承結構無法根據實際工況做出反饋調節,改善轉子系統的振動特性。

本文以機動飛行下的轉子系統為研究對象,考慮具有變剛度和變阻尼特性的磁流變阻尼器,基于有限元方法建立機動飛行下磁流變阻尼器?轉子耦合系統動力學模型,采用Newmark?β數值方法對動力學方程進行求解,研究磁流變阻尼器對沖擊載荷下轉子系統瞬態及穩態動力學特性的影響。

1 機動飛行下磁流變阻尼器?轉子系統有限元模型

如圖1所示為磁流變阻尼器支承下轉子系統的有限元模型。轉軸考慮為Timoshenko梁,并被劃分為10個梁單元和11個節點,每個節點具有四個自由度,即xoy平面內的平動和繞x軸、y軸的轉動。剛性圓盤作為一個集中質量單元疊加在相應節點上。Dd表示圓盤直徑,Bd表示圓盤厚度,Ds表示轉軸直徑,Ls表示轉軸長度。

基于Lagrange方程,轉子系統動力學方程可表示為[11]:

式中? ,,其中,x和y分別表示單元節點在x和y方向的位移,和分別表示單元節點繞x軸和y軸轉動的角度;ω為轉軸角速度;M,J和K分別為系統質量矩陣、陀螺矩陣和剛度矩陣;Q1和Q2為一般外力。

1.1 磁流變阻尼器動力學模型

坐標系XYZ(如圖2所示)用于定義流體流動的位置和速度分量。eB為內環偏心量,θ為相對于最小油膜位置處的油膜方位角,h為相應方位角下油膜的厚度,其可被近似表達為:

(2)

式中? c為初始徑向油膜間隙。

基于雙線性本構方程,磁流變阻尼器油膜壓力可表示為[12?13]:

當0≤Z≤Zc時:

當Zc≤Z≤L時:

式中? 定義αη=ηc/η為黏度比,η表示未施加電流作用時的磁流變液黏度,ηc表示硬核處磁流變液黏度;為油膜厚度對時間t的一階導數;Zc為硬核在Y方向上充滿整個油膜間隙的起始點軸向坐標位置(如圖3所示);為油膜壓力對變量Z的一階偏導;pc為油膜壓力在軸向Zc位置處的壓力值;C1和C2為積分常數;τc表示硬核邊界處磁流變液剪切應力;τy為磁流變液屈服應力,其與磁場強度之間的關系表達式為[14]:

式中? H=IN/(2h)為磁場強度,I為電流強度,N為線圈匝數;a1=0.0297;a2=19.75;a3=1.102×104;a4=-2.482×104;c1=-22.29;c2=2.601×104為常系數。

通過求解式(4)可獲得油膜沿軸向的壓力梯度,對沿軸向進行積分可獲得油膜壓力分布:

式中? pA為大氣壓。

則磁流變阻尼器的徑向油膜力和切向油膜力可表示為:

式中? R為磁流變阻尼器半徑;L為阻尼器軸向長度。

磁流變阻尼器油膜力在x和y方向上的分力可表示為:

式中? xe和ye分別為磁流變阻尼器內環中心在x和y方向上的位移。

1.2 轉子支承系統動力學模型

滾動軸承示意圖如圖4所示。軸承外圈與鼠籠彈性支承相連,內圈與轉軸相連,并隨轉軸旋轉。假設滾珠在軸承保持架上等間距排列,并做純滾動。

在圖4中,Ri為內滾道半徑,Ro為外滾道半徑,φi為第i個滾珠在t時刻的角位置,即:

式中? ωo為滾珠中心角速度;Nb為滾珠數量。

第i個滾珠與滾道的法向接觸變形δi可表示為:

(13)

式中? xb和yb分別為軸承內圈中心在x和y方向上的位移;r0為軸承游隙。

基于非線性Hertz接觸理論,軸承力模型可表示為[15]:

式中? Cb為Hertz接觸剛度;H(·)為Heaviside函數。

轉子系統支承結構如圖5所示。根據牛頓第二定律,磁流變阻尼器支承系統動力學方程可表示為[14]:

式中? mb為軸頸質量;mm為軸承外圈質量;ce為軸承處阻尼系數;ka為鼠籠剛度。

1.3 機動飛行動力學模型

如圖6所示為描述機動飛行的模型和坐標系。基于拉格朗日方程,祝長生等[4]提出了一種在任意機動飛行條件下,帶有不平衡的多剛性盤、多集總質量和多軸承的線性和非線性轉子系統一般動力學模型。參照[4]的方法,機動飛行下轉子?軸承系統的動力學方程可以表示為:

式中? 為轉軸單元節點的位移向量;,,,分別為轉軸單元的質量矩陣、剛度矩陣、阻尼矩陣和陀螺矩陣;為施加在單元節點上的一般外力。CB,i,KB,i和FB,i分別表示機動飛行引起的等效附加阻尼矩陣、剛度矩陣和附加力向量:

式中? 和分別為機動飛行所施加單元節點的直徑轉動質量和極轉動慣量;mi為機動飛行所施加單元節點的質量;X,Y以及Z分別表示機動飛行水平、垂直和移動方向;θX,θY和θZ分別表示繞X,Y和Z方向的旋轉角度;vX和aX分別為X方向上的速度和加速度;vθX和aθX分別為繞X軸旋轉的角速度和角加速度,它們在Y和Z軸上具有相似的表達方式。

2 結果與討論

為了研究磁流變阻尼器對機動飛行下轉子系統動力學行為的影響,利用Newmark?β法對動力學方程進行求解,可獲得系統動力學響應。仿真中,轉軸、磁流變阻尼器、滾動軸承以及機動飛行參數如表1~4所示。

2.1 機動飛行下轉子系統動力學響應分析

無機動飛行時轉子系統在圓盤位置處的幅頻特性曲線如圖7所示。振動幅值隨著轉速的增大而逐漸增大,并在ω=580 rad/s (92.31 Hz)時達到峰值0.65 mm。當轉速ω超過580 rad/s時,隨著轉速的增大,振動幅值逐漸減小。可以發現,580 rad/s是轉子系統的一階臨界轉速。

根據一階臨界轉速為92.31 Hz,分別設置亞臨界和超臨界狀態下的轉速為55.39 Hz(轉速比λ=0.6)和129.23 Hz(轉速比λ=1.4)。接下來詳細分析機動飛行過程中轉子系統在亞臨界和超臨界狀態下的動力學特性。

圖8所示為λ=0.6和λ=1.4時,機動飛行下轉子系統的時域響應。由圖8(a)可知,在機動飛行開始時,轉子系統在豎直方向上的位移瞬間增加到一個峰值,然后經過幾個周期的衰減后穩定到新的平衡位置。當機動飛行結束時,轉子系統再次經歷該過程,并返回到機動飛行前的穩態響應。這些現象表明機動飛行會對轉子系統產生瞬態沖擊效應。定義沖擊系數IFb和IFe來定量地描述機動飛行開始時和機動飛行結束時的沖擊效應:

(18)

式中? Amb和Ame分別為機動飛行開始時和機動飛行結束時轉子系統在豎直方向上的瞬時最大位移。Ad和Aa分別為機動飛行期間和機動飛行之后轉子系統在豎直方向上的穩態響應幅值。

通過對比分析圖8可獲得以下結論:

(1) 在機動飛行過程中,轉子系統在豎直方向上穩態響應的平衡位置從0 mm變為-0.50 mm,水平方向上穩態響應的平衡位置保持不變。這意味著由機動飛行引起的附加離心力僅影響轉子系統在豎直方向上的振動,而沒有耦合作用。

(2) 除平衡位置偏移以外,機動飛行產生的附加阻尼會降低轉子系統在機動飛行過程中的穩態響應幅值。豎直方向上的振幅從0.11 mm減小到0.06 mm,水平方向上的振幅從0.11 mm減小到0.08 mm,降幅分別約為45.5%和27.3%,這意味著與附加離心力的影響不同,附加阻尼效應具有耦合作用。

(3) 機動飛行期間轉子系統在豎直方向上和水平方向上的穩態響應幅值均高于機動飛行前的穩態響應幅值,這表明與亞臨界狀態不同,超臨界狀態下的附加阻尼會增大轉子系統的振動。

(4) 超臨界狀態下的IFb和IFe分別為3.62和2.77,低于亞臨界狀態下的13.00和4.00。這意味著在亞臨界狀態下運行的轉子系統對機動飛行引起的沖擊載荷更加敏感。

2.2 磁流變阻尼器支承下轉子系統動力學響應分析

圖9~11分別為亞臨界狀態下(λ=0.6),考慮與不考慮磁流變阻尼器作用時,轉子系統動力學響應的時間歷程圖、頻譜圖以及時頻圖。

(1) 機動飛行前,由時間歷程圖可知,與無磁流變阻尼器相比,磁流變阻尼器能為轉子系統提供有效的阻尼,使得轉子系統在y方向上的振動幅值降低。由頻譜圖可知,磁流變阻尼器作用下轉子系統主要表現為基頻成分,且與無磁流變阻尼器作用時相比,基頻幅值由0.104 mm降低為0.080 mm,降幅約為23.1%。

(2) 機動飛行中,由時間歷程圖可知,與無磁流變阻尼器相比,磁流變阻尼器作用下轉子系統在y方向上的振動幅值基本保持不變。由頻譜圖可知,考慮磁流變阻尼器作用時,頻譜圖中出現2×,3×等倍頻成分及非協調頻率成分,其中2倍頻最為明顯。其原因是在機動飛行中,由于附加離心力的影響,使得轉子系統在y方向上振動的平衡位置發生大幅偏移,油膜受過度受擠壓而非線性增強。

(3) 機動飛行開始和結束瞬態過程中,由時間歷程圖可知,考慮磁流變阻尼器作用時,機動飛行開始和結束時轉子系統在豎直方向上的瞬時最大位移Amb和Ame均減小,Amb由0.78 mm減小為0.73 mm,Ame由0.44 mm減小為0.41 mm。此外,由時頻圖可知,在機動飛行開始瞬間和機動飛行結束瞬間,除基頻外,時頻圖中還出現了較寬的頻率分量,且頻率越高衰減越快,沖擊引起的寬頻分量的能量主要集中在93.30 Hz,接近轉子系統一階固有頻率92.31 Hz。原因是機動飛行引起的瞬態沖擊負載具有較寬的頻帶,激發了轉子系統的一階固有頻率。考慮磁流變阻尼器作用時,沖擊引起的瞬態響應幅值減小,且瞬態響應時間縮短,以上現象說明磁流變阻尼器對機動飛行引起的瞬態沖擊效應具有一定的緩解作用。

圖12~14分別為超臨界狀態下(λ=1.4),考慮與不考慮磁流變阻尼器作用時,轉子系統動力學響應的時間歷程圖、頻譜圖以及時頻圖:

(1) 機動飛行前,由時間歷程圖可知,考慮磁流變阻尼器作用時,轉子系統在y方向上穩態響應的振動幅值變化微弱。頻譜圖主要表現為基頻成分,且幅值由0.196 mm增大到0.203 mm,增幅約為3.45%。

(2) 機動飛行中,由時間歷程圖可知,磁流變阻尼器對轉子系統振動幅值的影響較弱。頻譜圖中主要表現為基頻成分,且與不考慮磁流變阻尼器作用時相比,基頻幅值減小,由0.236 mm減小為0.226 mm,降幅約為4.24%。此外,頻譜圖中還出現微弱的2倍頻成分。

(3) 機動飛行開始和結束瞬態過程中,由時頻圖可知,超臨界轉速下,機動飛行引起的瞬態沖擊效應會激發轉子系統的一階模態響應,沖擊引起的寬頻分量的能量主要集中在93.65 Hz,接近轉子系統一階固有頻率92.31 Hz,考慮磁流變阻尼器作用時,其響應的幅值變化不明顯,但瞬態響應的時間縮短。

如圖15所示為亞臨界狀態下(λ=0.6)機動飛行前和機動飛行中轉子系統在豎直方向上穩態響應的級聯圖。隨著電流的增大,磁流變阻尼器對機動飛行前轉子系統穩態響應的幅頻特性影響較小,而對機動飛行中轉子系統穩態響應的幅頻特性影響較大。與機動飛行前穩態響應相比,機動飛行中轉子系統穩態響應的非線性明顯增強。由圖15(b)可知,當I≤0.4 A時,隨著電流的增大,基頻幅值增大,且非線性逐漸增強。這是因為,當電流較小時,磁流變液屈服應力較小,機動飛行產生的附加載荷使得軸頸離心率過大,油膜受過度擠壓而產生較強的非線性。當I>0.4 A時,磁流變液屈服應力增大,使得軸頸離心率減小,磁流變阻尼器的阻尼效應起主要作用。隨著電流的增大,阻尼力增強,基頻幅值逐漸減小,在I=1.0 A時,磁流變阻尼器的阻尼效應最為明顯。隨后,隨著電流的增大,基頻振動幅值逐漸增大,非線性逐漸增強。這是因為,當電流較大時,磁流變液屈服應力較大,磁流變阻尼器剛度效應增強并起主要作用,從而限制了油膜的擠壓作用,削弱了磁流變阻尼器的阻尼效應,使得振動幅值逐漸增大。同時,在較大電流作用下,磁流變液硬核體積增大,油膜變薄,致使油膜力非線性增強。當I>1.6 A時,磁流變阻尼器對轉子系統動力學特性的影響基本保持不變,這是因為,在大的電流作用下,磁流變阻尼器處于準剛性支承狀態。

如圖16所示為超臨界狀態下(λ=1.4)機動飛行前和機動飛行中轉子系統在豎直方向上穩態響應的級聯圖。隨著電流的增大,機動飛行前轉子系統穩態響應的基頻幅值逐漸增大,并在I=1.2 A時趨于穩定。與機動飛行前轉子系統穩態響應的級聯圖相比,機動飛行中轉子系統穩態響應的非線性明顯增強。當電流I≤1.2 A時,隨著電流的增大,基頻幅值逐漸增大,系統非線性逐漸增強,級聯圖中出現較為明顯的連續非協調頻率成分。當I>1.2 A時,系統處于準剛性支撐,振動狀態基本保持不變。通過對比圖15可以發現,亞臨界狀態下,磁流變阻尼產生的非線性頻率成分主要在基頻以上,且出現較為明顯的倍頻成分。而在超臨界狀態下,磁流變阻尼產生的非線性頻率成分主要在基頻以下,且表現為連續頻譜現象。

3 結? 論

(1) 機動載荷會對轉子系統產生沖擊效應,使轉子系統在機動飛行開始瞬間和機動飛行結束瞬間的振動幅值大幅增大。此外,機動飛行引起的沖擊載荷具有較寬的頻帶,會激發轉子系統振動的一階自然模態。

(2) 在亞臨界狀態下,磁流變阻尼器能夠有效抑制機動飛行前的穩態響應幅值,且對機動飛行引起的瞬態沖擊效應具有一定的緩解作用,能夠降低瞬態沖擊幅值,縮短響應時間。在超臨界狀態下,磁流變阻尼器對轉子系統動力學特性的影響較小。

(3) 在整個機動飛行過程中,通過施加合適的電流能夠有效抑制轉子系統的振動幅值及緩解機動飛行過程中的瞬態沖擊效應,但較大的電流反而會使磁流變阻尼器產生較強的非線性,造成轉子系統的失穩,且在機動飛行中,由于油膜被過度擠壓,磁流變阻尼器的非線性更為明顯。

參考文獻:

[1]????? Han B, Ding Q. Forced responses analysis of a rotor system with squeeze film damper during flight maneuvers using finite element method[J]. Mechanism and Machine Theory, 2018, 122: 233-251.

[2]????? Lin F, Meng G. Study on the dynamics of a rotor in a maneuvering aircraft[J]. Journal of Vibration and Acoustics, 2003, 125(3): 324-327.

[3]????? El-Saeidy F M A, Sticher F. Dynamics of a rigid rotor linear/nonlinear bearings system subject to rotating unbalance and base excitations[J]. Journal of Vibration and Control, 2010, 16(3): 403-438.

[4]????? 祝長生, 陳擁軍. 機動飛行時發動機轉子系統動力學統一模型[J]. 航空動力學報, 2009, 24(2): 371-377.

Zhu C S, Chen Y J. General dynamic model of aeroengine's rotor system during maneuvering flight[J]. Journal of Aerospace Power, 2009, 24(2): 371-377.

[5]????? Andrés L S, Rodríguez B. Experiments with a rotor-hybrid gas bearing system undergoing maneuver loads from its base support[J]. Journal of Engineering for Gas Turbines and Power, 2020, 142(11): 111004.

[6]????? Yu Y, Ding K, Zhao T, et al. Nonlinear dynamics of flexible diaphragm coupling's rotor system during maneuvering flight[J]. The Journal of Strain Analysis for Engineering Design, 2023, 58(3): 236-254.

[7]????? Gao T, Cao S, Sun Y. Nonlinear dynamic behavior of a flexible asymmetric aero-engine rotor system in maneuvering flight[J]. Chinese Journal of Aeronautics, 2020, 33(10): 2633-2648.

[8]????? Gao T, Cao S, Hou L, et al. An experimental study on the nonlinear vibration phenomenon of a rotor system subjected to barrel roll flight and coupled rub-impact faults[J]. Measurement, 2020, 153: 107406.

[9]????? Gao T, Yuan S, Liu P, et al. Vibration behavior of dual-rotor caused by maneuver load and intershaft bearing defect[J]. AIAA Journal, 2023, 61(3): 1396-1410.

[10]??? Chen X, Gan X, Ren G. Dynamic modeling and nonlinear analysis of a rotor system supported by squeeze film damper with variable static eccentricity under aircraft turning maneuver[J]. Journal of Sound and Vibration, 2020, 485: 115551.

[11]??? Wang J, Liu Y, Qin Z, et al. Dynamic performance of a novel integral magnetorheological damper-rotor system[J]. Mechanical Systems and Signal Processing, 2022, 172: 109004.

[12]??? Wang J, Liu Y, Qin Z, et al. Nonlinear characteristic investigation of magnetorheological damper-rotor system with local nonlinearity[J]. Chinese Journal of Aeronautics, 2023: 36(2): 111-126..

[13]??? Wang J, Zhang X, Liu Y, et al. Rotor vibration control via integral magnetorheological damper[J]. International Journal of Mechanical Sciences, 2023, 252: 108362.

[14]??? Wang J, Ma L, Zhang J, et al. Mitigation of nonlinear rub-impact of a rotor system with magnetorheological damper[J]. Journal of Intelligent Material Systems and Structures, 2020, 31(3): 321-338.

[15]??? 羅躍綱, 王鵬飛, 王晨勇, 等. 迷宮密封-滾動軸承-懸臂轉子系統非線性動力學特性分析[J]. 振動工程學報, 2020, 33(2): 256-264.

Luo Y G, Wang P F, Wang C Y, et al. Nonlinear dynamic characteristics of labyrinth seal-rolling bearing-cantilever rotor system[J]. Journal of Vibration Engineering, 2020, 33(2): 256-264.

Dynamic characteristics of rotor system with magnetorheological damper during maneuvering flight

Abstract: In the process of maneuvering flight, the aeroengine will bear very harsh working conditions, leading to irregular transient vibrations that can result in failure. In this paper, the effect of a semi-active magnetorheological damper (MR damper) on the dynamic characteristics of a rotor system under maneuvering flight is investigated. The finite element model of the rotor system with MR damper under maneuvering flight is established using the finite element method. The Newmark-β numerical method is used to solve the dynamic equations, and the dynamic characteristics of the rotor system during maneuvering flight are studied. On this basis, considering the effects of MR damper on the transient, the steady state responses of the rotor system under maneuvering flight are analyzed. The results show that transient impact is caused at the beginning and the end of maneuvering flight, which stimulates the first order modal response of the rotor system. The MR damper with suitable current can effectively suppress the amplitudes of transient and steady-state responses of the rotor system during maneuvering flight. In addition, due to the large eccentricity of the journal in maneuvering flight, the MR damper is prone to produce nonlinearity.

Key words: aeroengine; maneuvering flight;magnetorheological damper;rotor system;transient impact

主站蜘蛛池模板: 免费三A级毛片视频| 99热国产这里只有精品9九| 国产超碰一区二区三区| 国产亚洲男人的天堂在线观看| 日韩a级毛片| 久久亚洲精少妇毛片午夜无码| 91福利国产成人精品导航| 免费在线a视频| 国产乱人伦偷精品视频AAA| 国产亚洲精| 久久免费视频播放| 亚洲天堂免费| 日韩经典精品无码一区二区| 97综合久久| 人妻丰满熟妇啪啪| 亚洲一区第一页| а∨天堂一区中文字幕| 乱码国产乱码精品精在线播放| 国产成人AV大片大片在线播放 | 亚洲视频色图| 国产一区二区丝袜高跟鞋| 亚洲成人黄色在线观看| 午夜精品一区二区蜜桃| 国产精品视频白浆免费视频| 国产小视频网站| 亚洲h视频在线| 白丝美女办公室高潮喷水视频| 四虎在线高清无码| 女人爽到高潮免费视频大全| 亚洲精品在线91| 美女一级毛片无遮挡内谢| 亚洲精品色AV无码看| 亚洲第一区精品日韩在线播放| 极品尤物av美乳在线观看| 免费视频在线2021入口| 就去色综合| 日韩无码精品人妻| 中文字幕 欧美日韩| 欧美一级特黄aaaaaa在线看片| 国产无码网站在线观看| 亚洲欧美自拍中文| 色综合激情网| 九九热精品免费视频| 欧美在线黄| 欧美三级视频在线播放| 一区二区三区国产精品视频| 不卡国产视频第一页| av一区二区人妻无码| 啪啪国产视频| 乱人伦中文视频在线观看免费| 91黄视频在线观看| 亚洲大学生视频在线播放| 有专无码视频| 亚洲人成网站观看在线观看| 日韩精品一区二区三区swag| 久久青草免费91线频观看不卡| 日本AⅤ精品一区二区三区日| 午夜日韩久久影院| 亚洲综合欧美在线一区在线播放| 国产精品自在自线免费观看| 无码在线激情片| 免费网站成人亚洲| 欧美另类视频一区二区三区| 一本一道波多野结衣一区二区| 国产成人永久免费视频| 四虎国产在线观看| 国产成人亚洲综合A∨在线播放| 国产尤物在线播放| 77777亚洲午夜久久多人| 夜夜操天天摸| 日本成人一区| 国产主播在线一区| 亚洲不卡av中文在线| 国产麻豆福利av在线播放| 55夜色66夜色国产精品视频| 99久久精品免费看国产免费软件 | P尤物久久99国产综合精品| 国产日韩欧美在线视频免费观看| 欧美伦理一区| 国产v精品成人免费视频71pao| 国产精品99r8在线观看| 亚洲最新地址|