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

航空發動機雙轉子系統的模態及其表達方法

2022-11-21 03:39:36廖明夫程榮輝叢佩紅黃江博
振動與沖擊 2022年21期
關鍵詞:表達方法模態振動

王 瑞, 廖明夫, 程榮輝, 叢佩紅, 黃江博

(1.西北工業大學 動力與能源學院,西安 710072;2.航發沈陽發動機研究所,沈陽 110015)

航空發動機雙轉子-支承系統含有2個轉子,即高壓轉子和低壓轉子。通常,高、低壓轉子之間設置有中介軸承,用于支承高壓轉子的渦輪端,從而形成結構耦合。雙轉子系統為一多源激勵系統。由于存在高壓激勵源和低壓激勵源,因此,雙轉子系統的模態分為高壓轉子激勵的模態和低壓轉子激勵的模態兩類模態[1]。另外,雙轉子系統的模態和振動響應還與高、低壓轉子的轉速和轉速控制律,包括轉動方向緊密相關[2]。轉子系統模態的確定和清晰的表達是轉子動力學設計的核心任務,也是轉子故障診斷的基礎。目前,在轉子動力學設計時,一般會根據轉子系統的結構,建立有限元模型,進行數值計算,確定轉子的模態。在此過程中,需要判斷高壓轉子主激勵的正進動模態和反進動模態,以及低壓轉子主激勵的正進動模態和反進動模態,并予以清晰地表達。但對于對轉的雙轉子系統,在動力學設計時,模態的確定和表達往往會出現混淆的情況,造成判斷失誤。因此,建立準確和清晰的模態確定和表達方法,意義重大。

針對航空發動機雙轉子系統的模態特性,國內外許多學者進行了理論分析、計算和實驗驗證。傳遞矩陣法[3-5]和有限元法[6-9]是研究轉子模態特性的基礎,隨著計算機的發展,兩者的求解精度和計算效率均達到了較高水平,然而對雙轉子系統的雙源激勵和正、反進動的判別依然處于發展之中。繆輝等[10-11]利用商業軟件分析了輪盤的陀螺力矩對對轉雙轉子系統模態的影響,指出當輪盤作反進動時,陀螺力矩使轉軸當量剛度減小,轉子臨界轉速降低,并基于以上結論提出了判別高、低壓激勵的振型篩選法。張大義等[12-13]給出了基于振型篩選法的對轉雙轉子坎貝爾圖繪制方法,其根據轉子系統振型中占主導地位的轉子對正、反進動曲線進行篩選,當高低壓轉子耦合較小時可取得較好的效果,然而在耦合較強時則會造成較大的誤差。王存[14]基于商業軟件提出了可區分雙轉子坎貝爾圖中正、反進動曲線的完全法,以及在雙轉子系統具有簡單轉速關系時,無需迭代即可求出全部臨界轉速的直接法,最后以諧響應分析結果進行了驗證。上述文獻多是利用商業軟件對雙轉子動力學特性進行求解,并根據振型和坎貝爾圖中臨界轉速變化情況進行正、反進動判斷,未給出經過理論證明的統一判斷方法。

Gasch等[15]在其轉子動力學著作中,引入了正進動和反進動概念。Liao等[16-18]在研究帶裂紋轉子振動時,建立了進動分析方法,用于裂紋故障和不平衡故障的診斷,實驗結果證實了進動分析方法的有效性。針對支承各向異性的轉子,Liao等[19-23]又提出了進動比函數的概念,并利用進動比函數的特性來診斷轉子支座松動、轉軸裂紋和轉子不平衡故障,建立了利用進動比函數預估轉子特性參數的方法,提出了轉子進動的4個定理。上述文獻建立了較為完善的進動理論,但未系統性地提出雙轉子系統正、反進動模態的判斷方法。

Ferraris等[24]首次提出了計算坎貝爾圖的方法,之后一直應用于轉子模態特性的表達中,在單轉子或同轉雙轉子系統中,其斜率為正的曲線一般代表正進動,斜率為負的曲線一般代表反進動,但在對轉雙轉子系統中上述規律不再適用。張大義等繪制同轉雙轉子坎貝爾圖時,將高、低壓激起的各階臨界轉速繪制在一起,去掉反進動曲線,并同時代入高、低壓轉速關系曲線,兩條轉速線與正進動曲線的所有交點均為臨界轉速。由于未對高、低壓激勵進行區分,因此坎貝爾圖表達較為擁擠,且不適用于對轉雙轉子。于超平等[25]繪制了碰摩約束轉子模態頻率區間隨轉速變化的坎貝爾圖,可觀察發生碰摩故障時各轉速下的臨界轉速區間,但沒有對各階正、反進動進行區分。羅貴火等[26-28]在研究雙轉子系統動力學特性時,將雙轉子系統不平衡響應分別以低壓轉子和高壓轉子的進動頻率進行了分解,但表達在同一平面內,容易混淆。因此,判斷雙轉子系統正、反進動,對其響應進行分解,不僅需要理論分析,也需要清晰、直觀的表達方法,現階段的研究較少,模態判斷以及響應表達不直觀。

本文以典型的發動機雙轉子為對象,建立有限元模型,基于進動分析理論,建立雙轉子系統正、反進動模態的判斷方法,同時,建立雙轉子自振頻率和臨界轉速的“兩面三維”表達方法,既包含轉子正、反進動自振頻率,也包含高、低壓轉子轉速控制律,即共同工作線。另外,還建立了雙轉子系統不平衡響應的“三面三維”表達方法,在2個立面上,分別表示在某一測點測到的轉子振動信號中,高壓頻率分量和低壓頻率分量隨轉速的變化,平面上則表示轉速控制律,即共同工作線。這樣的表達形式,清晰地將雙轉子的模態與轉速、轉向和轉速控制律關聯在一起,為理解和分析轉子動力學特性提供了直觀和準確的圖景。

1 雙轉子模型及其模態特性

圖1表示典型的航空發動機雙轉子結構。低壓轉子有3個支承,高壓轉子則有2個支承,其中后支承為中介支承。

圖1 典型的航空發動機雙轉子結構Fig.1 Typical dual-rotor structure of aero-engine

以該雙轉子結構為原型,基于動力學相似,搭建了一套雙轉子模擬實驗器系統,如圖2所示。雙轉子模擬實驗器高、低壓轉子均通過電機同軸直驅,轉速調節范圍為0~6 000 r/min。低壓轉子為1-1-1支承方案,由低壓風扇盤、低壓軸以及低壓渦輪盤組成;高壓轉子為1-0-1支承方案,考慮到高壓電機占據了實驗器軸向空間,而高壓轉子為剛性轉子,在保證雙轉子系統動力學相似的基礎上,將發動機原型中高壓轉子的四級壓氣機盤整合為三級,即高壓轉子由高壓一級壓氣機盤、二級壓氣機盤、三級壓氣機盤、高壓軸以及高壓渦輪盤組成。實驗器的1,3和5支點為鼠籠式彈性支承結構,帶有擠壓油膜阻尼器。4支點為中介軸承,為避免高轉速下軸承游隙增大,其內環安裝在高壓轉子上,外環的軸承座安裝在低壓轉子上。雙轉子模擬實驗器系統的詳細參數已在文獻[29]中給出,本文不再贅述。

圖2 雙轉子模擬實驗器系統Fig.2 Dual-rotor test rig system

目前,在役和在研的有高、低壓轉子同轉型的發動機,也有對轉型的發動機。根據典型雙轉子發動機結構,建立有限元模型,可得到轉子如下的振動微分方程

(1)

式中:M為轉子系統質量矩陣;C為系統阻尼矩陣;GL為低壓轉子系統陀螺力矩矩陣;Gh為高壓轉子系統陀螺力矩矩陣;K為系統剛度矩陣;ΩL為低壓轉子轉速;Ωh為高壓轉子轉速;Q為系統外力向量;q為系統廣義位移向量,由所有節點的廣義位移向量依次排列,表示為

(2)

式中,下標為節點編號,N為節點總數。

求解轉子系統運動微分方程式(1)對應的齊次方程,便可得到雙轉子系統的自振頻率及對應的振型。

(3)

設h=h0eλt,則方程式(3)變為如下特征方程

(4)

式中:[I]為4N×4N階單位矩陣;[0]為4N×4N階零矩陣。

在一定的主激勵轉速下,解方程式(4)得到左邊矩陣的復特征根λ=α+jΩr,則帶阻尼進動頻率ωr和該模態下的阻尼比D為

(5)

特征值λ的實部α為衰減指數,虛部Ωr為轉子進動頻率。它們在復平面中的幾何關系,如圖3所示。

圖3 復平面中帶阻尼進動頻率和無阻尼進動頻率的幾何關系Fig.3 Geometric relationship between damped precession frequency and undamped precession frequency in complex plane

圖3中2個點表示復特征根總是以共軛形式存在的。

求解各個自轉角速度下的進動頻率時,取主激勵轉速(進動轉速)和高、低壓轉子轉速比,計算方程式(4)的特征值,就是該自轉角速度下的進動頻率。一個自轉角速度下有8N個特征值,即8N個進動頻率,連接各個自轉角速度下相同階進動頻率,得到雙轉子坎貝爾圖,如圖4所示。其橫坐標為主激勵自轉轉速Ω,縱坐標為計算出的各階進動頻率Ωr。由不平衡力激起共振時,其特點是激振力的頻率等于自轉角速度,在坎貝爾圖上作出Ωr=Ω,即斜率為1的直線,該直線與各進動頻率線的交點即為各階臨界轉速。它們對應的模態振型就是臨界轉速下的振型。在圖4(a)中Ω=ΩL,即低壓轉子主激勵下的坎貝爾圖,從下往上第1個、第3個、第4個和第5個交點對應同步反進動,第2個、第6個和第7個交點對應同步正進動;圖3(b)中Ω=Ωh,即高壓轉子主激勵下的坎貝爾圖,從下往上第1個、第3個、第4個和第7個交點對應同步反進動,第2個、第5個和第6個交點對應同步正進動。

(a) 低壓轉子主激勵

(b) 高壓轉子主激勵圖4 雙轉子坎貝爾圖Fig.4 Dual-rotor Campbell diagram

雙轉子系統中的2個轉子是以各自的轉速旋轉的,而轉速又是變化的。有時,高壓轉子和低壓轉子的轉速控制律不能用簡單的表達式表示,在求解頻率方程時將遇到困難。在這種情況下,可采用如下的計算方法。

在求解臨界轉速時,先給定低壓轉子的轉速,計算得到高壓轉子同步正進動的各階臨界轉速;改變低壓轉子的轉速,又可得到一組高壓轉子不平衡激起的各階臨界轉速。依此方法,求得在低壓轉子若干不同的轉速下,高壓轉子不平衡激起的臨界轉速。將同一階臨界轉速點連成曲線,得到該雙轉子系統由高壓轉子不平衡激起的臨界轉速,隨低壓轉子轉速的變化曲線,即HRE(high pressure rotor excitation)曲線。依照相同的方法,求得該雙轉子系統低壓轉子不平衡激起的臨界轉速,隨高壓轉子轉速的變化曲線,即LRE(low pressure rotor excitation)曲線。在這個圖譜中,繪出2個轉子的共同工作線,工作線與臨界轉速圖譜中各線的交點便是實際的各階臨界轉速。圖5顯示了雙轉子臨界轉速圖譜。點劃線是一條斜率為45°的直線,它與兩組曲線相交,且交點兩兩重合。這是因為,這條直線表示高壓轉子轉速等于低壓轉子轉速。這時,低壓主激勵和高壓主激勵是一致的。這條直線可以用來檢驗所計算的雙轉子臨界轉速圖譜是否正確。

圖5 雙轉子臨界轉速圖譜Fig.5 Dual-rotor critical speed diagram

圖5為雙轉子系統的臨界轉速圖譜,其中共同工作線與LRE曲線的交點為低壓激勵臨界轉速,共同工作線與HRE曲線的交點為高壓激勵臨界轉速。如果低壓轉子和高壓轉子的共同工作線發生變化,圖5中相應的各階臨界轉速也會發生變化。

圖5所示的結果只包含雙轉子系統正進動模態。在轉子動力學設計中還需確定雙轉子系統的反進動模態。

2 正、反進動模態的判斷

如第1章所述,由狀態方程對應的特征方程求得的特征根以復共軛的形式成對出現。此時需要判斷,對于給定的高、低壓轉速和轉速比,若以高壓轉速或以低壓轉速為參考基準,所求得的自振頻率對應的是正進動還是反進動。

根據轉子自轉轉速與公轉(進動)轉速的方向,可以將轉子進動分為正進動和反進動。利用有限元法計算出進動頻率后,需要根據對應的振型,即以上特征值問題的特征向量來判斷轉子的進動形式。

對于同轉雙轉子系統,正、反進動自振頻率比較容易判斷。一般情況下,高、低壓轉子同轉時,正進動自振頻率隨轉速升高而升高,反進動自振頻率隨轉速升高而降低。

但對于對轉雙轉子系統,上述變化規律并不完全成立。這是由于高、低壓轉子旋轉方向相反,所產生的陀螺效應也是相互反向的。例如,轉子系統某階模態是以高壓轉速反進動為主的模態,此時,高壓轉子的運動形式為協調反進動。而低壓轉子則發生的是非協調正進動。若此模態下,低壓轉子應變能占優,則隨轉速升高,對應的以高壓轉速反進動的自振頻率會升高。因此,對于對轉雙轉子,正、反進動自振頻率的判斷易出現混淆。為此,建立如下統一的判斷方法。

在第i階模態處,轉子的振動(狀態向量)則為

hi=a0i(h0eλit+h0i+1eλi+1t)

(6)

式中,a0i為常數,由初始條件決定。

將特征根代入式(6),可得

(7)

式中,括號外的時變系數a0eαt決定了轉子系統的穩定性。而模態的進動特性則由特征向量來確定。

特征向量h0i包含了轉子所有節點的振動速度和振動位移。設特征向量在轉子某一截面沿x和y方向的分量分別為

xi=xi,re+jxi,im

(8)

yi=yi,re+jyi,im

(9)

轉子的進動為

(10)

將式(8)和式(9)代入式(10),并將單位進動量ejΩrit和e-jΩrit的系數合并,得到如下的結果

(11)

式中:括號內第一項為正進動分量;第二項為反進動分量。

當如下條件成立時

(xi,re-yi,im)2+(xi,im+yi,re)2>

(xi,re+yi,im)2+(yi,re-xi,im)2

(12)

或進一步,當

xi,imyi,re-xi,reyi,im>0

(13)

轉子發生正進動,第i階模態為正進動模態。反之,當

xi,imyi,re-xi,reyi,im<0

(14)

轉子進動為反進動,第i階模態則為反進動模態。

上述判據中的量xi,re,xi,im,yi,re和yi,im在求解模態時可直接獲得,利用式(13)或式(14)就可很容易地判斷出模態的進動性質。

對于正進動模態,振型由轉子不同截面處正進動幅值和相位確定;而對于反進動模態,振型則由反進動幅值和相位確定。

3 正、反進動自振頻率的“兩面三維”表達方法

如第2章所述,雙轉子系統的模態包含高壓轉子主激勵的正進動模態和反進動模態,以及低壓轉子主激勵的正進動模態和反進動模態。每一模態不僅與主激勵頻率相關,還與高、低壓轉子轉速比,即共同工作線相關。圖4 表示的雙轉子坎貝爾圖描述了轉子自振頻率與主激勵轉子轉速的關系。但圖4無法確定每一個自振頻率所對應的轉速比,或另一個轉子的轉速。轉子共振發生的條件是:激振頻率與當前高、低壓轉子轉速所對應的自振頻率相等。另外,圖4將轉子正、反進動自振頻率隨主激勵轉速變化的曲線畫在同一象限,與主激勵頻率均有交點,容易混淆臨界轉速。圖5只描述了雙轉子系統正進動模態對應的自振頻率。

為更清晰地表達雙轉子系統正、反進動自振頻率和臨界轉速,建立如圖6所示的“兩面三維”表達方法。

(a) 低壓轉子主激勵

(b) 高壓轉子主激勵圖6 雙轉子自振頻率和臨界轉速的“兩面三維”表達方法Fig.6 “2 plane 3 dimension” interpretation method of the dual-rotor self-vibration frequency and the critical speed

圖6(a)中,圓圈“°”為低壓主激勵的臨界轉速及其對應的高壓轉子轉速;圖6(b)中,圓圈“°”為高壓主激勵的臨界轉速及其對應的低壓轉子轉速

圖6中,在垂直面(面1)上同時畫出主激勵下的正進動和反進動自振頻率,垂直軸為自振頻率,水平軸為主激勵轉子的轉速,同時,它將正進動和反進動自振頻率上下分開。在水平面(面2)上,畫出另一轉子的轉速坐標和共同工作線。這種“兩面三維”的表達方法物理意義非常清晰。圖6中自振頻率曲線上的任一點,都清晰地與高、低壓轉子轉速和共同工作線相關聯。不妨以低壓轉子主激勵的自振頻率曲線來說明。給定共同工作線,選擇一低壓轉速,就可通過共同工作線得到一高壓轉速。在此條件下,可計算得到轉子的各階自振頻率。以此類推,逐點求解,直到低壓轉子最高轉速,就可得到圖6中的曲線。而臨界轉速直接由直線ωL=ΩL與正進動自振頻率曲線的交點確定。直線ωL=ΩL與反進動自振頻率曲線無交點。自振頻率曲線是與給定的共同工作線關聯的。當共同工作線變化時,自振頻率曲線會隨之變化。

4 穩態不平衡響應求解與“三面三維”表達方法

在雙轉子系統中,高壓轉子、低壓轉子均存在不平衡激振力。由于氣動性能要求,高壓轉子/低壓轉子轉速必須符合特定的關系,即轉速控制律,或稱共同工作線。兩轉子轉速不同,因而高、低壓轉子不平衡力激振頻率不相等。兩轉子的不平衡力都能激起轉子系統的共振,共振時的轉速都是轉子系統的臨界轉速。通過中介軸承以及支承和機匣,載荷將由一個轉子傳遞給另一個轉子。在發動機上測得的振動包含低壓和高壓轉子的頻率成分(主要分量)。低壓轉子的不平衡響應是由低壓轉子不平衡量引起,高壓轉子的不平衡響應是由高壓轉子不平衡量引起,其運動形式有所不同。當低壓轉子不平衡力為主激勵時,低壓轉子作同步正進動,即自轉與公轉同步,并強迫高壓轉子以低壓轉速作公轉運動,使得高壓轉子公轉與自轉轉速不一樣,因此,高壓轉子作非同步正進動。如果高壓轉子與低壓轉子旋轉方向相反,高壓轉子則作非同步反進動。反之,高壓轉子不平衡量為主激勵激起的振動時,高壓轉子作同步正進動,低壓轉子作非同步正進動(或非同步反進動)。發動機上測得的響應是這兩種不平衡響應的疊加。

主激勵轉速(選定高壓轉速或者低壓轉速)為Ω,低壓轉速和高壓轉速分別為ΩL和Ωh,轉子系統的不平衡力可以表示為

Q=QccosΩt+QssinΩt

(15)

則雙轉子運動微分方程的穩態解可設為

q=qccosΩt+qssinΩt

(16)

代入式(1),并推導、整理得到

(17)

從而

q=qcos(Ωt-?)

(18)

其中

(19)

由以上推導可以看到,每給定一個主激勵轉速以及高壓轉速和低壓轉速的關系,即共同工作線,可以求出在給定不平衡量作用下雙轉子系統各節點的穩態響應。

如上所述,雙轉子具有高、低壓轉子主激勵的臨界轉速。在實際中,高、低壓轉子會同時存在不平衡量。在機匣上某一測點所測得的振動為高、低壓轉子不平衡激振力共同作用下的響應,同時還與轉速控制律相關。分析結果中,不易辨識高、低壓轉子主激勵的臨界峰值。為此,建立圖7所示的“三面三維”表達方法來描述雙轉子的不平衡響應。

圖7 雙轉子不平衡響應的“三面三維”表達方法Fig.7 “3 plane 3 dimension” interpretation method of the dual-rotor unbalanced response

圖7表征的是雙轉子模擬實驗器低壓渦輪測點的不平衡響應隨轉速和共同工作線的變化。實驗之前在高壓轉子和低壓轉子上分別添加了20 g·cm的不平衡量。雙轉子實驗器按照共同工作線增速,設計點1的低壓轉速為2 600 r/min,高壓轉速為4 810 r/min,設計點2的低壓轉速為4 000 r/min,高壓轉速為6 000 r/min。

當高、低壓轉子轉速按照共同工作線變化時,可從振動信號中分解出高壓基頻和低壓基頻分量,分別描繪在2個立面上,得到高壓轉子不平衡激勵下的響應曲線和低壓轉子不平衡激勵下的響應曲線。響應曲線上的峰值就分別對應高、低壓轉子主激勵的臨界峰值。圖7的水平面(面3)上畫出了轉速控制律,即共同工作線。還可以在水平面上將振動總量隨轉速的變化表達出來,即圖中的黑色曲線。“三面三維”表達方法,將雙轉子的振動與轉速、轉向和轉速控制律清晰地描繪在一起,便于理解和分析轉子的振動特性。

5 結 論

本文以典型的發動機雙轉子系統為對象,分析了雙轉子系統的模態特性,理論上證明了判斷正、反進動模態的條件,建立了雙轉子自振頻率和臨界轉速的“兩面三維”表達方法,以及雙轉子不平衡響應的“三面三維”表達方法,并應用于雙轉子模擬實驗器,效果顯著。本文結論如下:

(1) 雙轉子系統高、低壓轉子主激勵的正、反進動模態除與主激勵轉速相關外,還與高、低壓轉子轉速控制律(共同工作線)有關。利用常規的方法不能清晰表達雙轉子系統的模態特征。

(2) 本文建立的雙轉子自振頻率和臨界轉速“兩面三維”表達方法,是在垂直面(面1)上同時畫出主激勵下的正進動和反進動自振頻率,垂直軸為自振頻率,水平軸表示主激勵轉子的轉速,同時,它將正進動和反進動自振頻率上下分開。在水平面(面2)上,畫出另一轉子的轉速坐標和共同工作線。這樣的表達方法能將自振頻率與高、低壓轉子轉速和共同工作線清晰地關聯在一起。

(3) 雙轉子系統的不平衡響應是高、低壓轉子不平衡激勵的共同作用結果,同時還與轉速控制律相關。在對振動信號的分析中,不易辨識高、低壓轉子主激勵的幅頻特性。

(4) 本文建立的雙轉子不平衡響應“三面三維”表達方法利用3個垂直面分別表示高、低壓轉子主激勵下的幅頻特性,利用水平面表示高、低壓轉子共同工作線和振動總量隨高、低壓轉子轉速變化。“三面三維”表達方法,將雙轉子的振動與轉速、轉向和轉速控制律清晰地描繪在一起,便于理解和分析轉子的振動特性。

猜你喜歡
表達方法模態振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
如果藝術有一萬種表達方法
藝術啟蒙(2022年11期)2022-12-06 09:34:38
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
中立型Emden-Fowler微分方程的振動性
英語中序數詞的表達方法
系統醫學(2016年8期)2016-02-20 02:55:08
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
人腸系膜血管平滑肌細胞BKCa 通道在HEK293 細胞上的表達方法研究
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 午夜精品福利影院| 欧美a在线看| 99久视频| 日本a∨在线观看| 欧美在线网| 久久国产精品电影| 国产精品嫩草影院av| 91欧美亚洲国产五月天| 欧美自拍另类欧美综合图区| 日日噜噜夜夜狠狠视频| 在线日韩日本国产亚洲| 日本欧美精品| 欧美国产另类| 精品少妇人妻无码久久| 亚洲欧美日韩天堂| 午夜免费小视频| 欧美日韩在线亚洲国产人| 国产成人精品第一区二区| 成人福利在线观看| 日韩中文字幕亚洲无线码| 久久国产高清视频| 日韩精品无码不卡无码| 欧美精品亚洲精品日韩专| 超清无码熟妇人妻AV在线绿巨人 | 亚洲全网成人资源在线观看| 久久99这里精品8国产| 国产激情无码一区二区三区免费| 国产亚洲精久久久久久久91| 99久视频| 国产精品大白天新婚身材| 好紧好深好大乳无码中文字幕| 亚洲日韩精品综合在线一区二区| 性69交片免费看| 亚瑟天堂久久一区二区影院| 日韩在线第三页| 久久综合色播五月男人的天堂| 久久黄色影院| 伊在人亞洲香蕉精品區| 国产91精品久久| 国产亚洲男人的天堂在线观看| 自拍欧美亚洲| 99re这里只有国产中文精品国产精品| 99久久99视频| 999国内精品视频免费| 好吊妞欧美视频免费| 婷婷成人综合| 国产成人精品亚洲77美色| 天天躁夜夜躁狠狠躁躁88| AV不卡无码免费一区二区三区| 亚洲va在线观看| 毛片久久网站小视频| 欧美精品一二三区| 任我操在线视频| 成人免费午间影院在线观看| 欧美日韩综合网| 国产成人亚洲精品色欲AV| 思思99热精品在线| 国产精品无码AⅤ在线观看播放| 欧美区一区| 国产在线视频自拍| 久久精品午夜视频| 久久久久免费精品国产| 国产特级毛片aaaaaa| 97精品久久久大香线焦| 中文无码影院| 日韩福利视频导航| 午夜无码一区二区三区| 免费一级无码在线网站| 91久久国产综合精品女同我| 亚洲成肉网| a级免费视频| 男女男免费视频网站国产| a国产精品| 精品人妻无码区在线视频| 欧美一级专区免费大片| 性视频一区| 最新加勒比隔壁人妻| 精品视频一区在线观看| 一级毛片免费高清视频| 日本精品一在线观看视频| 欧美激情首页| 国产精品久久久久久久久久98|